一种新的基于卡尔曼滤波的点扩散函数估计方法

摘要

本发明为一种用于图像复原的基于Kalman模型的点扩散函数估计方法。本发明首先建立点扩散函数估计的对称全平面Kalman状态方程和观测方程；其次根据选择图像平稳区域，确定点扩散函数模型的观测矩阵；最后根据平稳区域进行Kalman迭代，估计出图像的点扩散函数。本发明有效改善了点扩散函数估计不准确和非自适应的缺点，可以根据图像自适应的估计点扩散函数，估计的点扩散函数数值更准确。为后续图像复原提供了更准确的点扩散函数。

说明书

技术领域

本发明涉及图像处理领域，具体涉及一种基于Kalman模型的点扩散函 数估计方法。

背景技术

图像在获取过程中，受噪声和光照等因素的影响，存在散焦现象。在数 学上，散焦过程可以用点扩散函数PSF英文解释来描述，即未散焦图像f(x, y)和点扩散函数h(x,y)相卷积，得到了模糊图像g(x,y)，g(x,y)即为实际中 获取的图像，这里仅考虑散焦模糊的影响。图像复原的方法已有不少文献报 道，目前备受关注的有两类方法：基于参数估计的方法和基于自适应滤波的 方法，前者根据某种估计准则，比如最大似然准则和贝叶斯准则，寻求最大 似然函数，从而得到的图像和PSF就是复原的图像和估计的PSF，但是PSF 估计的结果较为均匀，不能反映实际的PSF。后者采用状态空间描述其数学 公式，利用观测数据求状态向量各个分量的最小二乘估计，通过递推计算复 原的图像，但是PSF必须是已知的。

分析已有的图像复原方法，发现PSF的准确估计大多需要基于某种先验分 布，比如高斯分布，student-t分布等，才能够估计出比较准确的PSF数值， 进而获得较准确的复原图像。然而受噪声和光照等因素的影响，PSF不一定 服从事先假设的高斯等分布，这就存在估计的PSF不准确问题。因此基于先 验分布PSF估计方法的自适应性和抗噪声能不强。

发明内容

为避免以上现有技术的不足，本发明提出一种新的基于卡尔曼滤波的点 扩散函数估计方法，以解决图像在散焦情况下，实现引起散焦的点扩散函数 的自适应估计，提高点扩散函数估计的准确性，增强图像复原的准确性。

本发明的目的通过以下技术方案来实现：

本发明的目的通过以下技术方案来实现：

一种新的基于卡尔曼滤波的点扩散函数估计方法,该方法包括如下步 骤：

1）选择点扩散函数模糊后的图像平稳区域作为点扩散函数估计区域；

2）对模糊后的图像建立用于点扩散函数估计的Kalman状态方程和观测 方程；

3）根据所述模糊后的图像、状态方程和观测方程，计算出图像的观测矩 阵；

4）对所述状态方程和观测方程中的参数进行初始化；

5）根据步骤3估计的观测矩阵、状态方程和观测方程，进行Kalman迭 代计算，得到引起图像模糊的点扩散函数PSF；

6）对所述点扩散函数PSF数值进行调整。

进一步，所述步骤1）对图像平稳区域的选取方式如下：

经过大量实验，当灰度变化不超过图像灰度最大值15%时，方差不超过 30，该区域Θ可以近似为平稳随机过程。15%和30等具体数值可根据不同 类型的图像进行调整。不局限于于15%和30等具体数值。

进一步，所述步骤2）对点扩散函数估计的Kalman状态方程和观测方程 的建立包括如下步骤：

201）建立Kalman状态方程p(k+1)=A(k+1,k)p(k)+ξ(k)；该方程不局限于 线性方程。其中，p(k)表示点扩散函数，A(k+1,k)表示状态转移矩阵，ξ(k)表 示状态噪声。

202）建立Kalman观测方程g(k)=G(p(k))+η(k)，G(p(k))表示对p(k)的非 线性运算；该方程是非线性方程，在一定条件下，可以近似为线性方程 g(k)=C(k)p(k)+η(k)，其中C(k)是观测矩阵；g(k)的观测的散焦图像，为已 知量；η(k)是测量噪声。

进一步，所述步骤3）对观测矩阵的计算包括如下步骤：

301）将步骤1）确定的平稳区域设为Θ

302）令Θ是表示步骤1）所确定的平稳区域，B(k)∈Θ，B(k)可以看作 是Θ的一个样本；

303）B(k)在Θ中取值，B(k)维数和p(k)维数相同，观测矩阵C(k)是B(k) 中元素的组合，其选择规则根据卷积定理实现。

进一步，所述步骤4）中对参数进行初始化,具体是：误差协方差矩阵P (0|0)=0，0表示全0矩阵，状态估计表示全1矩阵；R_ξ=10^-4，R_η=10， R_ξ和R_η分别为状态误差和测量误差的方差。设定迭代次数K。

进一步，所述步骤5）对点扩散函数Kalman迭代计算包括如下步骤：

501）根据步骤4）中的初始化参数，计算k+1时刻的状态估计

502）根据当前时刻的测量矩阵C(k+1)、R_η和P(k+1|k)，计算增益矩阵H(k+1)

503）计算k+1时刻的PSF估计值

504）由P(k+1|k)和H(k+1)计算校正的误差协方差矩阵P(k+1|k+1)

505）如果没有达到迭代次数K，返回501）重复以上步骤；如果达到 迭代次数K，就是估计的点扩散函数PSF。

进一步，所述步骤6）对点扩散函数数值的调整方式如下：对估计得到 的点扩散函数PSF，对其数值进行筛选排序，去除出现的过大或者过小值， 将其变为从内至外逐渐衰减的排列。本发明的优点在于：

1.采用Kalman滤波的方法估计点扩散函数，具有自适应的优点，可以 只适应的根据图像散焦的程度估计点扩散函数

2.不受点扩散函数先验分布的约束

3.Kalman模型可以根据实际需要进行更改，灵活性较大。

附图说明

图1：图像的获取和复原降噪过程；

图2：本发明语音转换合成方法的具体步骤流程图；

图3：点扩散函数PSF图示；

图4：观测矩阵B(k)的变化过程图示。

具体实施方式

基于AR模型的Kalman图像复原方法将图像视为平稳随机场，而实际 上，图像是非平稳的，这样的近似带给图像复原更多的噪声。因此本文采取 无AR模型的Kalman复原方法复原图像，有效的降低了图像在复原过程中 的噪声，并有效保持了图像的边缘信息。

散焦的过程可以看作是未散焦图像f(x,y)和点扩散函数h(x,y)的卷积， 如图1所示，卷积后的图像，受加性噪声n(x,y)影响，形成最终的观测图像 g(x,y)。

聚焦的过程就是估计出PSF，采用估计的PSF复原图像并去除噪声，得 到f(x,y)，该过程就是图像复原。其中模糊包括运动模糊和散焦模糊，摄像 机位置固定，仅考虑散焦模糊。n(x,y)是加性高斯白噪声。复原系统p(x,y) 对模糊图像g(x,y)进行聚焦和降噪处理，得到估计的原始图像f(x,y)。该过 程可以通过式(1)表示

$g_{i,j}=\underset{l}{\Sigma }\underset{m}{\Sigma }{h}_{i-l,j-m}{f}_{l,m}+n_{i,j}$(1)

${\hat{f}}_{i,j}=F\left(g_{i,j}\right)$

其中g_i，j,h_i，j,f_i，j和n_i，j分别表示观测图像、模糊系统、原始图像和噪声。 F是复原降噪系统。本文的目的就是仅从散焦和噪声污染图像g_i，j中复原图 像f(x,y)。

如图2所示为本发明语音转换合成方法的具体步骤流程图。下面对本发 明语音转换合成方法的具体步骤进行详细说明。

步骤一：选择图像平稳区域进行点扩散函数估计

将获得的观测图像看作非因果、对称全平面、高斯马尔科夫随机场的一 个采样，设图像大小为M×M。假设点扩散函数PSF的大小是L=3×3(J×J， J＜M)，如图3所示。进行PSF估计时，需要选取图像中变化均匀的区域， 近似为平稳区域。经过大量实验，当灰度变化不超过图像灰度最大值的15% 时，方差不超过30，该区域Θ可以近似为平稳随机过程，假设该区域大小 为J×N，J＜N＜M。

步骤二：建立用于点扩散函数估计的Kalman状态方程和观测方程

令p(k)是点扩散函数PSF在k时刻的估计，p(k)是L×1的向量，p(k)= [p₁₁ p₂₁ p₃₁ p₁₂ p₂₂...p₃₃]^T，其中T表示转置，p(k)初始值设置为1/L，以保证 p(k)矩阵中所有元素相加的结果为1。建立如式(2)所示的状态方程，其中 A(k+1,k)为状态转移矩阵，表示状态向量p(k)从k时刻到k+1时刻的转移， 由于p(k)在转移过程中仅受噪声影响，因此A(k+1,k)为单位阵；ξ(k)为状态 误差，是对状态p(k)的修正。

p(k+1)=A(k+1,k)p(k)+ξ(k)    (2)

g(k)=G(p(k))+η(k)    (3)

建立如式(3)所示的观测方程，g(k)表示在k时刻的观测值，η(k)为测量 噪声，ξ(k)和η(k)可以认为满足以下的统计特性：

E(ξ(k))=0  E(ξ(k)ξ*(j))=R_ξδ_kj

E(η(k))=0  E(η(k)η*(j))=R_ηδ_kj    (4)

E(ξ(k)η*(j))=0

其中E(·)表示数学期望，R_ξ和R_η分别为状态误差和测量误差的方差。 它们在Kalman估计过程中是统计后的单一数值，经过大量实验，对于PSF 估计，R_ξ=10^-4，R_η=10。

${\delta }_{\mathrm{kj}}=\left(\begin{array}{cc}1& k=j\\ 0& k\ne j\end{array}\right)---\left(5\right)$

步骤三：观测矩阵的计算

式(3)是一个非线性方程，其中G(p(k))=C(k)p(k)，观测矩阵C(k)是L×L矩 阵，B(k)是空变的L×1向量，B(k)=[b₁₁ b₂₁ b₃₁ b₁₂ b₂₂…b₃₃]^T，B(k)∈Θ，Θ是观 测图像中的某块区域，可近似为平稳随机过程。B(k)可以看作是Θ的一个样本， 因此G(p(k))可以近似认为是线性的。

图4中的区域Θ是平稳随机过程，k时刻B(k)为左边黑色方框的像素点， 在k+1时刻，B(k)变化为右边黑色方框的像素点，新的像素点b₁₄ b₂₄ b₃₄进入 到方框中，形成B(k+1)，设矩阵c是L×1向量，c=[c₁₁ c₂₁ c₃₁ c₁₂ c₂₂…c₃₃]^T。 令c=B(k+1)，

[c₁₁,c₂₁,c₃₁,c₁₂,c₂₂,c₃₂,c₃₁,c₃₂,c₃₃]^T=[b₁₂,b₂₂,b₃₂,b₃₁,b₃₂,b₃₃,b₁₄,b₂₄,b₃₄]^T    (6)

观测矩阵C(k)如式(7)所示：

$C\left(k\right)=\left(\begin{array}{ccccccccc}0& 0& 0& 0& c_{11}& c_{21}& 0& c_{12}& c_{22}\\ 0& 0& 0& c_{11}& c_{21}& c_{31}& c_{12}& c_{22}& c_{32}\\ 0& 0& 0& c_{21}& c_{31}& 0& c_{22}& c_{32}& 0\\ 0& c_{11}& c_{21}& 0& c_{12}& c_{22}& 0& c_{13}& c_{23}\\ c_{11}& c_{21}& c_{31}& c_{12}& c_{22}& c_{32}& c_{13}& c_{23}& c_{33}\\ c_{21}& c_{31}& 0& c_{22}& c_{32}& 0& c_{23}& c_{33}& 0\\ 0& c_{12}& c_{22}& 0& c_{13}& c_{23}& 0& 0& 0\\ c_{12}& c_{22}& c_{32}& c_{13}& c_{23}& c_{33}& 0& 0& 0\\ c_{22}& c_{32}& 0& c_{23}& c_{33}& 0& 0& 0& 0\end{array}\right)---\left(7\right)$

随着k的更新，观测矩阵C(k)也在不断更新。

步骤四：初始化参数

对所述状态方程和观测方程中的参数进行初始化，初始化状态估计 在0时刻的状态估计值表示全1矩阵，初始化误差协方差 矩阵P(0|0)=0，0表示全0矩阵。及对状态误差的方差R_ξ，测量误差的方差 R_η进行初始化，R_ξ=10^-4，R_η=10，并设定迭代次数K。

步骤五：根据步骤3得到的观测矩阵，进行Kalman迭代计算

基于Kalman的PSF估计步骤如下所示：

假定k时刻状态估计值为其误差协方差矩阵为P(k|k)，k+1时刻 的状态估计为

其中是根据状态方程从k时刻的估计结果外推得到

$\hat{p}(k+1|k)=A(k+1,k)\hat{p}\left(k\right|k)---\left(9\right)$

W(k+1)是第k个观测量包含的新息，即实际观测值和模型方程计算值的 差：

$W(k+1)=g(k+1)-C(k+1)\hat{p}(k+1|k)---\left(10\right)$

新息W(k+1)与k时刻之前所有时刻的观测数据g(1)，...g(k)正交，且不同 时刻的新息也是正交的。即新息是具有白噪声性质的随机过程，提供了观测 数据g(k)所包含的新的信息。式(6)中的H(k+1)是增益矩阵，其计算也基于 正交原理：

H(k+1)={P(k+1|k)C(k+1)^T}{C(k+1)P(k+1|k)C(k+1)^T+R_η}^-1    (11)

其中P(k+1|k)是k时刻未经校正的误差协方差矩阵，

P(k+1|k)=A(k+1,k)P(k|k)A^T(k+1,k)+R_ξ    (12)

校正的误差协方差矩阵P(k+1|k+1)可以表示为：

P(k+1|k+1)={I-H(k+1)C(k+1)}P(k+1|k)    (13)

在进行PSF估计前，首先初始化在0时刻的状态估计值1表示全1矩阵，误差协方差矩阵P(0|0)=0，0表示全0矩阵。

所述基于Kalman的PSF估计步骤实现方法具体如下：

(1)初始化矩阵：如步骤四所示；

(2)根据P(0|0)=0和公式(12)，计算未经校正的误差协方差矩阵P (k+1|k)；

(3)根据当前时刻的测量矩阵C(k+1)、R_η和P(k+1|k)，由公式(11)计算 增益矩阵H(k+1)；

(4)根据公式(9)和公式(8)计算k+1时刻的PSF估计值

(5)由P(k+1|k)、H(k+1)和公式(13)计算校正的误差协方差矩阵P (k+1|k+1)；

(6)k＝k+1，返回(1)；当k达到设定的最大值，完成Kalman估计。

通过上述过程，获得点扩散函数

步骤六：点扩散函数的整理

在某些情况下，Kalman估计过程中，点扩散函数的估计值中有可能出 现估计的过大或者过小值。按照点扩散函数数值从内至外逐渐减小的特性， 每一圈数值去掉最大值和最小值，剩余值取平均，实现Kalman估计的点扩 散函数的稳定性。

得到的点扩散函数是一个J×J的矩阵，1＜J。点扩散函数矩阵的数值从 内至外可能有大有小，但是理论上，PSF从内至外，其数值应当呈衰减的趋 势。所以对估计得到的点扩散函数PSF，对其数值进行筛选排序，去除可能 出现的过大或者过小值，将其变为从内至外逐渐衰减的排列。以图3为例进 行说明，p₂₂是中心点，保持不变；距p₂₂最近的四个点p₁₂，p₂₁,p₃₂和p₂₃，去 掉这四个点中的最大值和最小值后，对剩余的两个值取平均作为距p₂₂最近 的四个点的值，即这四个点的值是相同的；最后距中心点p₂₂较远的四个点 p₁₁，p₁₃,p₃₁和p₃₃，去掉其中的最大值和最小值，剩余值取平均后作为最终的 p₁₁，p,₁₃，p₃₁和p₃₃。当点扩散函数矩阵变大时，同样采用这种方法对PSF矩阵 进行排序。

应当理解，以上借助优选实施例对本发明的技术方案进行的详细说明是 示意性的而非限制性的。本领域的普通技术人员在阅读本发明说明书的基础 上可以对各实施例所记载的技术方案进行修改，或者对其中部分技术特征进 行等同替换；而这些修改或者替换，并不使相应技术方案的本质脱离本发明 各实施例技术方案的精神和范围。