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

Abstract

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

Description

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

技术领域

本发明涉及图像处理领域，具体涉及一种基于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}{\mathrm{\Σ}}\underset{m}{\mathrm{\Σ}}{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。

${\mathrm{\δ}}_{\mathrm{kj}}=\left\{\begin{array}{cc}1& k=j\\ 0& k\≠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矩阵进行排序。

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

Claims (7)

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

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

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

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

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

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

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

2.根据权利要求1所述的一种新的基于卡尔曼滤波的点扩散函数估计方法,其特征在于：所述步骤1）对图像平稳区域的选取方式如下：

经过大量实验，当灰度变化不超过图像灰度最大值15%时，方差不超过30，该区域Θ可以近似为平稳随机过程。

3.根据权利要求1所述的一种新的基于卡尔曼滤波的点扩散函数估计方法,其特征在于：所述步骤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)是测量噪声。

4.根据权利要求1所述的一种新的基于卡尔曼滤波的点扩散函数估计方法，其特征在于：所述步骤3）对观测矩阵的计算包括如下步骤：

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

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

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

5.根据权利要求1所述的一种新的基于卡尔曼滤波的点扩散函数估计方法,其特征在于：所述步骤4）中对参数进行初始化,具体是：误差协方差矩阵P(0|0)=0，0表示全0矩阵，状态估计

表示全1矩阵；R_ξ＝10^-4，R_η=10，R_ξ和R_η分别为状态误差和测量误差的方差，设定迭代次数K。

6.根据权利要求1所述的一种新的基于卡尔曼滤波的点扩散函数估计方法,其特征在于：所述步骤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。

7.根据权利要求1所述的一种新的基于卡尔曼滤波的点扩散函数估计方法,其特征在于：所述步骤6）对点扩散函数数值的调整方式如下：对估计得到的点扩散函数PSF，对其数值进行筛选排序，去除出现的过大或者过小值，将其变为从内至外逐渐衰减的排列。

Images