CN103514589A - 图像恢复中的滤波方法 - Google Patents

Abstract

本发明公开了一种图像恢复中的滤波方法，包括以下步骤：1）对待处理图像进行分块预处理；2）建立图像在格林空间的状态空间模型；3）生成图像的固定滞后平滑器并对图像进行滤波处理。本发明将格林空间的滤波计算方法应用于图像处理中，在计算量上相比增广卡尔曼滤波算法大大降低，在未知噪声统计特性的情况下实现图像恢复，有效的提高了图像噪声滤波处理的适用性以及图像恢复的鲁棒性。

Description

图像恢复中的滤波方法

【技术领域】

本发明属于计算机图像处理领域，具体涉及一种图像恢复中的滤波方法。其中，图像恢复是指通过计算机处理，对质量下降的图像加以重建或恢复的处理过程。

【背景技术】

图像滤波是图像恢复中的一个重要步骤。现有技术中比较常用的是Hilbert空间的卡尔曼（Kalman）滤波。众所周知，卡尔曼滤波是一个不断地预测、修正的递推滤波过程，由于其在求解时不需要存储大量的观测数据，并且当得到新的观测数据时，可随时算得新的参数滤波值，便于实时的处理观测结果，因此被越来越多地应用于动态定位数据处理中。然而，在Hilbert空间使用Kalman滤波对系统进行滤波处理，我们需要事先知道系统的准确模型及噪声的统计特性，这在实际生产生活中都是很困难的。

同时，卡尔曼滤波方法在初始迭代过程中，滤波器所得结果是与初值的选取密切相关的，一般认为该滤波方法在滤波的初始阶段误差较大、效果也差；只有当时间充分大时，滤波器的效果才与初值的选取无关，而滤波器趋于稳定状态。

【发明内容】

本发明所要解决的技术问题是：提出一种图像恢复中的滤波方法，解决卡尔曼滤波要求事先已知噪声的统计特性的问题，用来解决图像恢复问题。

本发明的技术问题通过以下的技术方案予以解决：

一种图像恢复中的滤波方法，包括以下步骤：1）对待处理图像进行分块预处理；2）建立图像在格林空间的状态空间模型；3）生成图像的固定滞后平滑器并对图像滤波。

本发明与现有技术对比的有益效果是：

本发明将格林空间的H∞滤波（即固定滞后平滑器）理论应用到图像恢复中，研究证实，H∞滤波在一定条件下就是格林空间卡尔曼滤波,因此可将Hilbert空间的滤波问题转化为格林空间的卡尔曼滤波问题，用来解决图像恢复问题。这种格林空间的滤波理论可以在没有事先知道噪声统计特性的情况下，对图像进行滤波并实现图像恢复。

另外，采用H∞滤波器的滤波方法使得滤波误差始终维持在一个固定范围内，能够保证在整个图像范围内滤波效果都是较好的，包括初始阶段，有效的提高了滤波器对图像噪声的广泛适用性以及图像恢复的鲁棒性。

【附图说明】

图1是本发明实施例的技术方案流程框图。

图2是本发明实施例对待处理图像分块的原理图。

【具体实施方式】

下面结合具体实施方式并对照附图对本发明做进一步详细说明。

如图1所示，为本具体实施方式中的图像恢复中基于格林空间的固定滞后平滑器的滤波方法。利用格林空间中滤波器进行图像恢复的模块流程图，包括以下步骤：

U1）：对待处理图像进行分块预处理；

此步骤中，包含对待处理图像进行分块，以及确定分块后图像的每个块所代表的物理信息及状态向量的关系。

具体地，对图像的分块预处理过程（如图2所示）可以表述为：每4个相邻的像素点按照如下规则整理成大小是4×1的列向量（X_1,1(i),BX_3,3(i)），然后选取这样的9个相邻向量组合得到一个36×1的状态向量X(i)。（上述数字4和9可以改为其它数，是数字n的平方即可，n≥2,但n最好不要太大超过7,如果太大会影响滤波效果。）

X_1,1(i)=[x_1,1(i) x_2,1(i) x_1,2(i) x_2,2(i)]^T

X_2,1(i)=[x_3,1(i) x_4,1(i) x_3,2(i) x_4,2(i)]^T

X_3,1(i)=[x_5,1(i) x_6,1(i) x_5,2(i) x_6,2(i)]^T

X_1,2(i)=[x_1,3(i) x_2,3(i) x_1,4(i) x_2,4(i)]^T

X_2,2(i)=[x_3,3(i) x_4,3(i) x_3,4(i) x_4,4(i)]^T

C

X_3,3(i)=[x_5,5(i) x_6,5(i) x_5,6(i) x_6,6(i)]^T

$X\left(i\right)={\left[\begin{array}{ccccccc}{X}_{\mathrm{1,1}}^T\left(i\right)& X_{\mathrm{2,1}}^T\left(i\right)& X_{\mathrm{3,1}}^T\left(i\right)& X_{\mathrm{1,2}}^T\left(i\right)& X_{\mathrm{2,2}}^T\left(i\right)& B& X_{\mathrm{3,3}}^T\left(i\right)\end{array}\right]}^T$

其中，列向量X_m,n(i)(1≤m,n≤3)是由第m,n个小方块内所包含的4个像素点的灰度值所组成的，而列向量X(i)是由9个小方块对应的列向量X_m,n(i)组成的。

至此，依照该方法实现了对图像的分块，将大型的图像数据矩阵分成若干个小型矩阵，使得我们对图像的处理更加容易实现，也可在提高运算速度的同时减少计算量。

U2）：建立图像在格林空间的状态空间模型；

此步骤中，要建立图像在格林空间的状态空间模型。首先应依据U1）中对图像的分块结果以及每个向量所含的物理意义，列写系统模型的基本方程；其次再确定模型中各系数矩阵的具体表达式。

根据格林空间系统模型的一般表达式，我们令图像具有如下形式的观测方程和量测方程：

X(i+1)=FX(i)+Gu(i) (1)

Y(i)=HX(i)+v(i) (2)

其中，F,G分别具有如下形式：

$F={\left[\begin{array}{ccc}0& I& 0\\ 0& 0& I\\ 0& 0& 0\end{array}\right]}_{9\×9},$ $G={\left[\begin{array}{ccc}0& 0& 0\\ 0& 0& 0\\ G_1& G_2& G_3\end{array}\right]}_{9\×9}$

一般在实际图像中，通过图像的拍摄、传输等过程都存在对图像模糊化的现象。因此，矩阵H是参考点扩散函数（PSF）的特点来进行选取的。通常，矩阵H是由点扩散系数h组成的，而h满足下式：

$h=\left\{\begin{array}{cc}1/\left({\mathrm{\πr}}^2\right)& x^2+y^2\≤r^2\\ 0& x^2+y^2>r^2\end{array}\right.$

其中，r是模糊半径。

在式(1)(2)所表示的系统中，X(i)表示每个块内的所有像素点灰度值，X(0)为初始状态，控制量u用来处理图像中所含的噪声、模糊等干扰因素，可选取它具有如下形式：

u=[u₁ u₂ u₃]^T

其中的v(i)表示能量有界的外部扰动，其统计特性未知。Y表示已知的观测量。

根据以上的系统模型公式(1)(2)，我们可以设计相应的基于格林空间的平滑器来实现对图像的滤波并达到图像恢复的目的。

此处，我们将图像的状态空间模型建立在格林空间内。其突出优点是，在格林空间内，允许噪声v(i)的内积为负定的（即v(i)v^T(i)e0），这在普通的Hilberlt空间内是不被允许的，这样更便于理论上对滤波器的设计。

另外，我们假设控制量u(i)和外部扰动v(i)与零均值的白噪声具有相同的统计特性，即有

$E\left\{\left[\begin{array}{c}u\left(i\right)\\ v\left(i\right)\end{array}\right]\left[\begin{array}{cc}{u}^T\left(j\right)& v^T\left(j\right)\end{array}\right]\right\}=\left[\begin{array}{cc}Q& S\\ S^T& R\end{array}\right]{\mathrm{\δ}}_{i,j}$

同时，关于初始状态X(0)的假设如下：初始状态X(0)独立于u(i)和v(i)，EX(0)=μ₀，并且

E[(X(0)-μ₀)(X(0)-μ₀)^T]=P₀

其中，符号E为数学期望，T为转置号，δ_i,j=1(i=j),δ_i,j=0(i≠j)。

此步骤中，还应确定滤波器（即基于格林空间的平滑器）存在的条件(因为我们设计的滤波器是在一定的性能指标下存在的，所以要确定这种滤波器存在需要满足的条件)。首先，初始条件X(0)记作X₀，控制量u满足关系：u=u₀+u_u，这里u₀是已知量，u_u是未知量。

其次，准则函数的应具有的形式：

X₀已知时满足：

$\underset{u_u,w\&Element;L_2[0,N-1],X_0\&Element;R^n}{\sup }\frac{{\left|\right|u-\hat{u}\left|\right|}_2^2}{{\left|\right|u_u\left|\right|}_2^2+{\left|\right|w\left|\right|}_2^2+{X_0}^T{\mathrm{SX}}_0}<{\mathrm{\&gamma;}}^2$

X₀未知时满足：

其中，表示控制量u的估计值，w表示噪声，γ是滤波器鲁棒性能指标。

滤波器存在的条件是满足如下的Riccati方程：

P_j=FP_jF^T+GG^T-FP_jH^T(I+HP_jH^T)^-1HP_jF^T

且，P₀=Π₀。满足以上条件的图像才存在格林空间滤波器。

U3）：设计图像的固定滞后平滑器对图像进行滤波处理，实现图像恢复。

此步骤中，我们将对U2）中确定的系统方程进行格林空间的滤波器设计。涉及到的方法是基于经典Kalman滤波器和Mendel的输入白噪声估值器，设计相应的有限域固定滞后平滑器，并最终利用设计的平滑器对图像滤波以实现图像恢复。

有限域固定滞后平滑器的问题可描述为：基于观测(y(t+N),y(t+N-1),B,y(1)),N>0,求状态x(t)的固定滞后平滑器， $\hat{x}\left(t\right|t+1),t=\mathrm{0,1},B;$ N固定。

采用有限域固定滞后平滑器，根据后面块的图像信息估计前面块的图像信息。此处，我们可以选取系统（1）（2）的控制量u具有白噪声的形式。因为，系统是完全可观、可控的，且F是稳定矩阵，所以有相应的关于控制量u的稳态输入估值器，其形式如下：

$\hat{u}\left(i\right|i+N)=\hat{u}\left(i\right|i+N-1)+M_N\&CenterDot;\mathrm{\&epsiv;}(i+N)---\left(3\right)$

带初值

且稳态增益M_N为

$M_N={\mathrm{QF}}^T{\left[{(I_N-\mathrm{KH})}^T{F}^T\right]}^T{H}^T{Q}_{\mathrm{\&epsiv;}}^{-1}$

其中新息ε(i)的稳态方差阵为Q_ε=H∑H^T+R

由以上信息，可设计带有白噪声估值器的有限域稳态次优固定滞后平滑器具有如下形式的：

$\hat{X}\left(i\right|i+N)=F\hat{X}(i-1|i-1+N)+G\hat{u}(i-1|i-1+N)+K_N\&CenterDot;\mathrm{\&epsiv;}(i+N)$

其中，i代表图像的第i块，K_N是稳态平滑增益阵，控制量的估值器如上面(3)式所示。根据以上的迭代关系，我们可知平滑器

关于滤波初值

是渐进稳定的。通过以上的步骤，根据第(i+N)块的图像信息观测值，实现了对图像第i块的平滑滤波。

本具体实施方式中，通过步骤U1)-U3）即完成对一幅图像的图像恢复。按照上述方法对图像进行分块预处理，然后设计相应的系统模型，并最终依据设计的系统模型建立对图像的固定滞后平滑器，实现图像恢复。

与Hilbert空间卡尔曼滤波的局限性相比，具有鲁棒性的H∞估计器（即H∞滤波器、基于格林空间的平滑器、固定滞后平滑器）具有如下的显著优势：1）它对噪声不敏感也因此使得系统更鲁棒；2）系统的外部扰动的统计信息不是必须已知的；3）估计器的唯一要求是扰动是能量有界的。在最近的研究探索中，我们发现鲁棒H∞估计器和控制问题以及几个相关的问题都可以统一到一个与Kalman滤波相接近的问题，不是在通常的希尔伯特（Hilbert）空间，而是一个特殊的称作格林空间的不定度量空间。因此，本发明用基于格林空间的鲁棒H∞估计器来进行图像恢复，可具有上述优点。

而且，与Hilbert空间不同的是，格林空间优点中允许向量内积是负定，即v(i)v^T(i)e0）的，更利于对图像的分析处理。

此外，传统的滤波器判断滤波器存在条件是由线性矩阵不等式（LMI）得到充分条件，然而本发明中采用的格林空间滤波器可使用黎卡提（Riccati）方程确定滤波器存在的充分条件，相比于传统滤波器更具优势。

通过这种设计固定滞后平滑器对图像进行平滑滤波的方法，更能充分利用图像前后块之间的相关信息实现对图像的滤波处理，使得图像恢复的效果比以往方法对噪声更具鲁棒性、图像更清晰。

Claims (5)

1.一种图像恢复中的滤波方法，其特征在于：包括以下步骤：

1）对待处理图像进行分块预处理；

2）建立图像在格林空间的状态空间模型；

3）生成图像的固定滞后平滑器并对图像进行滤波处理，实现图像恢复。

2.根据权利要求1所述的图像恢复中的滤波方法，其特征在于：所述步骤1）中对待处理图像进行分块预处理时，分块的方法如下:每a²个相邻的像素点按照规则整理成大小是a²×1的列向量

然后选取这样的b²个相邻向量组合得到一个b²×1的状态向量X(i)，其中a、b为整数，≧2但<7。

3.根据权利要求1所述的图像恢复中的滤波方法，其特征在于：所述步骤2）中建立图像的状态空间模型的方法包括如下步骤：

1）接收用户根据步骤1）中对图像的分块方式，确定的每个块内所含信息的物理意义，和列写的系统模型的基本方程；

2）接收用户确定的模型中各系数矩阵的具体表达式，确定的状态向量信息，建立图像在格林空间的状态空间模型。

4.根据权利要求1所述的图像恢复中的滤波方法，其特征在于：所述步骤3）中图像的固定滞后平滑器具有如下形式的观测方程和量测方程：

X(i+1)=FX(i)+Gu(i) (1)

Y(i)=HX(i)+v(i) (2)

其中，F,G分别具有如下形式：

$F={\left[\begin{array}{ccc}0& I& 0\\ 0& 0& I\\ 0& 0& 0\end{array}\right]}_{9\&times;9},$ $G={\left[\begin{array}{ccc}0& 0& 0\\ 0& 0& 0\\ G_1& G_2& G_3\end{array}\right]}_{9\&times;9}$

矩阵H是参考点扩散函数PSF的特点来进行选取的由点扩散系数h组成，而h满足下式：

$h=\left\{\begin{array}{cc}1/\left({\mathrm{\&pi;r}}^2\right)& x^2+y^2\&le;r^2\\ 0& x^2+y^2>r^2\end{array}\right.$

其中，r是模糊半径;

X(i)表示每个块内的所有像素点灰度值，X(0)为初始状态，控制量u用来处理图像中所含的噪声、模糊等干扰因素，可选取它具有如下形式：

u=[u₁ u₂ u₃]^T

其中的v(i)表示能量有界的外部扰动，其统计特性未知,Y表示已知的观测量。

5.根据权利要求1所述的图像恢复中的滤波方法，其特征在于：所述步骤3）中生成图像的固定滞后平滑器前，根据黎卡提（Riccati）方程求得固定滞后平滑器存在的充分条件。