CN102968776B - 一种线性与非线性滤波器相结合的非均匀性校正方法 - Google Patents

Abstract

本发明公开了一种线性与非线性滤波器相结合的非均匀性校正方法，在对线性校正图像进行空间平滑滤波得到真实场景的估计图像时，如果非均匀性水平较高，使用模板为11x11的均值滤波，否则，使用模板为5x5的非线性平滑滤波，在实现边缘保持滤波的同时，计算校正参数的自适应迭代步长，进行异常像素（坏元、冲激噪声）检测与替换，以滤波输出图像作为对真实场景的估计图像，利用变化参考图像实现变化检测，只有当校正图像与变化参考图像的差异大于变化阈值时才对非均匀性校正系数进行迭代更新。本发明既具有线性平滑滤波器平滑作用大、收敛速度快的优点，又具有非线性平滑滤波的边缘保持特性，而且计算复杂度较低，适合硬件电路实现。

Description

一种线性与非线性滤波器相结合的非均匀性校正方法

技术领域

本发明属于图像探测处理领域，更具体地，涉及一种线性与非线性滤波器相结合的非均匀性校正方法。

背景技术

由于制造过程中的差异，非制冷红外焦平面阵列(Focal Plane Array，简称FPA)探测器中的每个探测元的光谱响应特性存在差异，称其为非均匀性。非均匀性表现为探测器输出图像中叠加在场景上的固定图案噪声(Fixed Pattern Noise，简称FPN)。该噪声严重影响图像质量，不利于目标检测、识别等应用，需要通过校正进行抑制甚至消除。由于非均匀性随时间以及环境条件的变化而漂移，仅在使用前进行一次基于标准黑体的标定是不够的。但在使用过程中周期性进行基于黑体的标定不但需要中断工作状态，而且需要额外的光学、机械设备。基于场景的方法可利用图像中的场景或目标运动实现非均匀性校正，不需要额外设备，不需要中断工作状态，因而得到广泛研究与应用。

在基于场景的非均匀校正方法中，最小均方(Least Mean Square，简称LMS)方法计算复杂度低、获得非均匀参数的速度快、并可持续跟踪非均匀性参数变化。与常统计量(Constant Statistics，简称CS)方法(J.G.Harris and Y.-M.Chiang,“Minimizing the Ghosting Artifact inScene-Based Nonuniformity Correction,”in SPIE Conference onInfrared Imaging Systems:Design Analysis,Modeling,and Testing IX,vol.3377,Orlando,Florida,1998)或基于图像配准的方法(Chao Zuo,“Scene-based nonuniformity correction algorithm based oninterframe registration”,J.Opt.Soc.Am.A/Vol.28,No.6/June2011)相比，LMS非均匀性校正方法在硬件电路实现实时校正中更具优势。

LMS非均匀性校正方法最早由Scribner提出，其基本思想是用空间平滑滤波器对当前校正图像进行滤波，作为对真实场景的估计。将误差函数定义为当前校正图像与估计图像的平方差，并采用最陡下降法迭代求解使误差函数最小的非均匀性校正系数(Adaptive Nonuniformity Correctionfor IR Focal Plane Arrays using Neural Networks：Proceedings of theSPIE:Infrared Sensors:Detectors,Electronics,and SignalProcessing,T.S.Jayadev,ed.,vol.1541,pp.100–109)。Scribner方法中先后使用过4邻域均值滤波器以及21x21均值滤波。Vera在2005年提出利用未进行非均匀性校正图像的局部方差自适应地控制LMS算法的迭代步长，从而在边缘区域减小迭代步长，在平坦区域增大迭代步长(FastAdaptive Nonuniformity Correction for Infrared Focal-Plane ArrayDetectors,”URASIP Journal on Applied Signal Processing 13,1994–2004)。Vera采用3x3均值滤波作为空间平滑滤波器。Hardie于2009年进一步改进Vera的方法，加入变化检测门控功能，只有当像素值发生明显变化时才执行非均匀性参数迭代，从而抑制静止目标的退化与鬼影问题(Scene-Based Non-uniformity Correction with Reduced Ghosting Usinga Gated LMS Algorithm”，OPTICS EXPRESS,Vol.17,No.17,17August2009,14918)。Hardie使用的空间平滑滤波器为21x21高斯滤波器。Rossi在2009年提出使用7x7双边滤波器作为空间平滑滤波器(Bilateral filter-based adaptive nonuniformity correction for infraredfocal-plane array systems”，Optical Engineering 49(5),057003)。该方法能够克服线性滤波器在边缘区域估计误差较大，易造成鬼影的问题，但无法从根本上抑制“鬼影”的出现。

LMS非均匀性校正方法中通常对校正图像进行空间平滑滤波，并将滤波结果作为对真实场景的估计。上述方法大多采用固定模板大小的线性滤波器，如均值滤波器或高斯滤波器。模板尺寸较大的滤波器，非均匀性校正范围较大，收敛速度较快，但估计误差较大。相反地，模板尺寸较小的滤波器估计误差较小，但校正范围较小，收敛速度较慢。线性滤波器在图像边缘处存在边缘模糊，容易出现目标退化及“鬼影”现象。非线性滤波器，如双边滤波器具有边缘保持能力，因此对目标退化及“鬼影”现象具有较强抑制作用。但双边滤波器计算复杂度较高，利用硬件电路实现实时校正比较困难。

发明内容

针对现有技术的缺陷，本发明的目的在于提供一种新的非均匀性校正方法。该方法具备线性平滑滤波器平滑效果明显、收敛速度快的优点，以及非线性平滑滤波器的边缘保持特性，同时能够克服线性平滑滤波器易造成“鬼影”，非线性平滑滤波器平滑效果差、收敛速度慢的缺点。该方法不需要计算局部方差即可实现校正参数迭代步长的自适应调整，而且具有坏元、冲激噪声检测与替换功能。另外，该方法具有较低的计算复杂度，可用硬件电路实现实时非均匀性校正。

为实现上述目的，本发明提供了一种线性与非线性滤波器相结合的非均匀性校正方法，包括以下步骤：

步骤101：逐帧接收红外焦平面输出图像序列，对每一帧当前图像，依据图像中每个像素位置的当前增益、偏移校正参数值进行线性校正，得到校正图像，校正公式如下：

y_ij(n)＝g_ij(n)x_ij(n)+o_ij(n)

其中y_ij(n)为校正图像，x_ij(n)为焦平面输出图像，i和j分别表示像素的行、列位置，n表示当前图像在红外焦平面输出图像序列中的序号，g_ij(n)、o_ij(n)分别为当前增益、偏移校正系数。如果当前图像为第一帧，且已对红外焦平面阵列进行过基于黑体的标定，则以标定得到的增益、偏移校正值作为g_ij(n)和o_ij(n)，否则取g_ij(n)＝1，o_ij(n)＝0。如果当前图像不是第一帧，则采用前一帧的校正参数更新值作为g_ij(n)和o_ij(n)；

步骤102：判断当前图像非均匀性水平是较高还是较低，如果当前图像非均匀性水平较高，则转入步骤103；如果当前图像非均匀性水平较低，则转入步骤104；

步骤103：对校正图像进行模板为11x11的均值滤波，并以滤波后的图像作为真实场景的估计图像，以校正图像与估计图像的差值作为误差图像，并计算误差图像的绝对均值，用当前图像的误差图像绝对均值与前一帧图像的误差图像绝对均值之比衡量误差变化程度，如果该比值接近1，则将当前图像非均匀性水平stuck递增1，否则将stuck清零，具体公式如下：

e_ij(n)＝y_ij(n)-t_ij(n)

$\mathrm{ame}\left(n\right)=\frac{1}{\mathrm{MN}}\underset{i=1}{\overset{M}{\mathrm{\Σ}}}\underset{j=1}{\overset{N}{\mathrm{\Σ}}}\left|e_{\mathrm{ij}}\left(n\right)\right|$

ratio(n)＝mae(n)/mae(n-1)

如果1-α≤ratio(n)≤1+α，则stuck＝stuck+1，否则stuck＝0

其中y_ij(n)为校正图像，t_ij(n)为估计图像，e_ij(n)为误差图像，ame(n)为误差图像绝对均值，ratio(n)为当前误差图像绝对均值与前一帧误差图像绝对均值的比值，α为预先设置的常数，且0＜α＜0.1；

步骤104：对校正图像进行模板为5x5的非线性平滑滤波，以滤波输出图像作为对真实场景的估计图像，并将校正图像与估计图像相减，得到误差图像，具体公式为e_ij(n)＝y_ij(n)-t_ij(n)，其中y_ij(n)为校正图像，t_ij(n)为估计图像，e_ij(n)为误差图像；

步骤105：将校正图像与变化参考图像相减后取绝对值，将所得差值图像的每个像素值与变化阈值比较，如果像素值大于变化阈值，则设置迭代步长为预先设定的固定迭代步长或步骤104中得到的自适应迭代步长，同时更新变化参考图像中对应像素位置的值；否则，设置迭代步长为0，同时保持变化参考图像中对应像素位置的值不变。具体公式如下：

${\mathrm{\η}}_{\mathrm{ij}}^{\′}\left(n\right)=\left\{\begin{array}{c}{\mathrm{\η}}_{\mathrm{ij}}\left(n\right),|t_{\mathrm{ij}}\left(n\right)-z_{\mathrm{ij}}\left(n\right)|>T\\ 0,\mathrm{else}\end{array}\right.$

$z_{\mathrm{ij}}(n+1)=\left\{\begin{array}{c}{t}_{\mathrm{ij}}\left(n\right),|t_{\mathrm{ij}}\left(n\right)-z_{\mathrm{ij}}\left(n\right)|>T\\ z_{\mathrm{ij}}\left(n\right),\mathrm{else}\end{array}\right.$

其中η_ij(n)为11x11均值滤波器使用的固定迭代步长或5x5非线性平滑滤波的自适应迭代步长，η′_ij(n)为实际迭代步长，t_ij(n)为校正图像，z_ij(n)为变化参考图像，z_ij(n+1)为用于下一帧的变化参考图像。T为变化阈值。z_ij(n)的初始值设置为大于输入图像最大值；

步骤106：执行非均匀性校正参数的更新，具体公式如下：

g_ij(n+1)＝g_ij(n)-η′_ij(n)e_ij(n)·y_ij(n)

o_ij(n+1)＝o_ij(n)-η′_ij(n)e_ij(n)

其中g_ij(n+1)为更新后的增益校正参数，o_ij(n+1)为更新后的偏移校正参数，η′_ij(n)为步骤105得到的迭代步长，e_ij(n)为步骤103或104得到的误差图像，y_ij(n)为步骤101得到的校正图像；

步骤107：重复以上步骤101至106，直到红外焦平面输出图像的所有帧都处理完毕为止。

步骤102中的具体判决规则为，如果当前图像为第一帧，则默认非均匀性水平较高，如果当前图像不是第一帧，且当前图像非均匀性水平<2，则非均匀性水平较高，如果当前图像不是第一帧，且当前图像非均匀性水平＝＝2，则非均匀性水平较低。

步骤104进一步包括以下子步骤：

步骤201：对校正图像进行非线性平滑滤波，具体为采用5×5滑动窗口按照光栅扫描顺序逐像素对校正图像进行滤波操作，具体公式如下：

$t(i,j)=y(i,j)+\underset{m=-2}{\overset{2}{\mathrm{\&Sigma;}}}\underset{n=-2}{\overset{2}{\mathrm{\&Sigma;}}}{w}_{\mathrm{mn}}{\mathrm{\&Delta;}}_{\mathrm{mn}}(i,j)$

${\mathrm{\&Delta;}}_{\mathrm{mn}}(i,j)=\left\{\begin{array}{c}y(i+m,j+m)-y(i,j),|y(i+m,j+n)-y(i,j)|\&le;T\\ T,y(i+m,j+n)-y(i,j)>T\\ -T,y(i+m,j+n)-y(i,j)<-T\end{array}\right.$

其中t为滤波图像，y为校正图像，T为灰度差阈值；

步骤202：利用非线性平滑滤波器进行异常像素检测，具体公式如下：

$N_{\mathrm{mn}}(i,j)=\left\{\begin{array}{c}1,|y(i+m,j+n)-y(i,j)|\&le;T\\ 0,|y(i+m,j+n)-y(i,j)|>T\end{array}\right.$

$M(i,j)=\underset{m=-2}{\overset{2}{\mathrm{\&Sigma;}}}\underset{n=-2}{\overset{2}{\mathrm{\&Sigma;}}}{N}_{\mathrm{mn}}(i,j)$

若当前像素的有效像素计数值M(i,j)小于等于预设的阈值，即M(i,j)≤T_M，则判断(i,j)为异常像素，该位置的校正图像像素值、目标图像像素值均用最邻近已处理像素的值替换，替换规则如下：

$y(i,j)=\left\{\begin{array}{c}y(i,j-1),j\&NotEqual;1\\ y(i-1,j),j=1\end{array}\right.$

$t(i,j)=\left\{\begin{array}{c}y(i,j-1),j\&NotEqual;1\\ y(i-1,j),j=1\end{array}\right.$

其中y(i,j)为像素位置(i,j)的校正图像，y(i,j-1)为像素位置(i,j-1)的校正图像，y(i-1,j)为像素位置(i-1,j)的校正图像，t(i,j)为像素位置(i,j)的目标图像，i为行序号，j为列序号；

步骤203：利用非线性平滑滤波器计算当前像素的迭代步长，具体公式如下：

$\mathrm{\&eta;}(i,j)=\frac{k\max }{1+256\&times;\mathrm{ratio}(i,j)},$

ratio(i,j)＝1-M(i,j)/25

其中kmax为迭代步长可取的最大值。

通过本发明所构思的以上技术方案，与现有技术相比，本发明具有以下的有益效果：

(1)由于在步骤102到104中采用了线性与非线性空间平滑滤波器相结合的方法，本发明方法摆脱了其他同类方法采用固定空间平滑滤波器的缺陷。在初始校正阶段非均匀性水平较高时，使用11x11均值滤波器在较大范围内平滑非均匀性，提高非均匀性校正参数的收敛速度；在非均匀性水平降低到较低水平后，使用具有边缘保持特性的5x5非线性滤波器克服非线性滤波器在边缘区域模糊效应，从而降低非均匀参数估计误差，提高“鬼影”抑制能力。因此，本发明方法的收敛速度比其他单独使用线性平滑滤波器或非线性平滑滤波器的方法更快，非均匀性参数估计准确性更高，“鬼影”抑制能力更强；

(2)由于采用了步骤202，本发明方法可实现异常像素，如坏元、冲激噪声污染像素坏检测/替换。其他同类方法需要从另外的坏元检测单元或者外部提供的坏元模板获得坏元信息才能实现坏元替换。此外，同类算法不具备冲激噪声检测/替换能力；

(3)由于采用了步骤203，本发明方法可实现非均匀性校正参数迭代步长的自适应调整。其他同类算法通常需要计算图像的局部方差才能实现非均匀性校正参数迭代步长的自适应调整。

附图说明

图1为本发明线性与非线性滤波器相结合的非均匀性校正方法的流程示意图。

图2为本发明非线性滤波器实现平滑滤波、异常像素检测以及迭代步长自适应调整的功能示意图。

图3为本发明与其他代表性方法对模拟数据进行非均匀性校正的绝对误差均值比较。

图4(a)至(h)为本发明与其他代表性方法对模拟数据进行非均匀性校正的“鬼影”现象比较。

图5为本发明与其他代表性方法对实际红外序列进行非均匀性校正的收敛速度比较。

图6(a)至(g)为本发明与其他代表性方法对实际红外序列进行非均匀性校正的“鬼影”现象比较。

图7(a)至(c)为本发明与其他代表性方法的冲激噪声抑制能力比较。

具体实施方式

为了使本发明的目的、技术方案及优点更加清楚明白，以下结合附图及实施例，对本发明进行进一步详细说明。应当理解，此处所描述的具体实施例仅仅用以解释本发明，并不用于限定本发明。

如图1所示，本发明线性与非线性滤波器相结合的非均匀性校正方法包括：

步骤101：逐帧接收红外焦平面输出图像序列，对每一帧当前图像，依据图像中每个像素位置的当前增益、偏移校正参数值进行线性校正，得到校正图像，校正公式如下：

y_ij(n)＝g_ij(n)x_ij(n)+o_ij(n)

其中y_ij(n)为校正图像，x_ij(n)为焦平面输出图像，i和j分别表示像素的行、列位置，n表示当前图像在红外焦平面输出图像序列中的序号，g_ij(n)、o_ij(n)分别为当前增益、偏移校正系数。如果当前图像为第一帧，且已对红外焦平面阵列进行过基于黑体的标定，则以标定得到的增益、偏移校正值作为g_ij(n)和o_ij(n)，否则取g_ij(n)＝1，o_ij(n)＝0。如果当前图像不是第一帧，则采用前一帧的校正参数更新值作为g_ij(n)和o_ij(n)；

步骤102：判断当前图像非均匀性水平是较高还是较低，如果当前图像非均匀性水平较高，则转入步骤103；如果当前图像非均匀性水平较低，则转入步骤104。具体判决规则如下：

规则1：如果当前图像为第一帧(n＝1)，则默认非均匀性水平较高；

规则2：如果当前图像不是第一帧(n≠1)，且当前图像非均匀性水平stuck<2，则非均匀性水平较高；

规则3：如果当前图像不是第一帧(n≠1)，且当前图像非均匀性水平stuck＝＝2，则非均匀性水平较低。

在以上规则中，stuck的初始值为0；

步骤103：对校正图像进行模板为11x11的均值滤波，并以滤波后的图像作为真实场景的估计图像，以校正图像与估计图像的差值作为误差图像，并计算误差图像的绝对均值，用当前图像的误差图像绝对均值与前一帧图像的误差图像绝对均值之比衡量误差变化程度，如果该比值接近1，则将stuck递增1，否则将stuck清零，具体公式如下：

e_ij(n)＝y_ij(n)-t_ij(n)

$\mathrm{ame}\left(n\right)=\frac{1}{\mathrm{MN}}\underset{i=1}{\overset{M}{\mathrm{\&Sigma;}}}\underset{j=1}{\overset{N}{\mathrm{\&Sigma;}}}\left|e_{\mathrm{ij}}\left(n\right)\right|$

ratio(n)＝ame(n)/ame(n-1)

如果1-α≤ratio(n)≤1+α，则stuck＝stuck+1，否则stuck＝0

其中y_ij(n)为校正图像，t_ij(n)为估计图像，e_ij(n)为误差图像，ame(n)为误差图像绝对均值，ratio(n)为当前误差图像绝对均值与前一帧误差图像绝对均值的比值，α为预先设置的常数(通常取0＜α＜0.1)；

步骤104：对校正图像进行模板为5x5的非线性平滑滤波，具体操作见步骤201到203。以滤波输出图像作为对真实场景的估计图像，并将校正图像与估计图像相减，得到误差图像。具体见以下公式：

e_ij(n)＝y_ij(n)-t_ij(n)

其中y_ij(n)为校正图像，t_ij(n)为估计图像，e_ij(n)为误差图像

步骤105：将校正图像与变化参考图像相减后取绝对值，将所得差值图像的每个像素值与变化阈值比较，如果像素值大于变化阈值，则设置迭代步长为预先设定的固定迭代步长或步骤104中得到的自适应迭代步长，同时更新变化参考图像中对应像素位置的值；否则，设置迭代步长为0，同时保持变化参考图像中对应像素位置的值不变。具体公式如下：

${\mathrm{\&eta;}}_{\mathrm{ij}}^{\&prime;}\left(n\right)=\left\{\begin{array}{c}{\mathrm{\&eta;}}_{\mathrm{ij}}\left(n\right),|t_{\mathrm{ij}}\left(n\right)-z_{\mathrm{ij}}\left(n\right)|>T\\ 0,\mathrm{else}\end{array}\right.$

$z_{\mathrm{ij}}(n+1)=\left\{\begin{array}{c}{t}_{\mathrm{ij}}\left(n\right),|t_{\mathrm{ij}}\left(n\right)-z_{\mathrm{ij}}\left(n\right)|>T\\ z_{\mathrm{ij}}\left(n\right),\mathrm{else}\end{array}\right.$

其中η_ij(n)为11x11均值滤波器使用的固定迭代步长或5x5非线性平滑滤波的自适应迭代步长，η′_ij(n)为实际迭代步长，t_ij(n)为校正图像，z_ij(n)为变化参考图像，z_ij(n+1)为用于下一帧的变化参考图像。T为变化阈值。z_ij(n)的初始值设置为大于输入图像最大值；

步骤106：执行非均匀性校正参数的更新。计算公式如下：

g_ij(n+1)＝g_ij(n)-η′_ij(n)e_ij(n)·y_ij(n)

o_ij(n+1)＝o_ij(n)-η′_ij(n)e_ij(n)

其中g_ij(n+1)为更新后的增益校正参数，o_ij(n+1)为更新后的偏移校正参数，η′_ij(n)为步骤105得到的迭代步长，e_ij(n)为步骤103或104得到的误差图像，y_ij(n)为步骤101得到的校正图像。

步骤107：重复以上步骤101至106，直到红外焦平面输出图像的所有帧都处理完毕为止。

如图2所示，本发明方法中步骤104进一步包括以下子步骤：

步骤201：对校正图像进行非线性平滑滤波。具体操作如下：

采用5×5滑动窗口按照光栅扫描顺序逐像素对校正图像进行滤波操作。滤波方法如下：

$t(i,j)=y(i,j)+\underset{m=-2}{\overset{2}{\mathrm{\&Sigma;}}}\underset{n=-2}{\overset{2}{\mathrm{\&Sigma;}}}{w}_{\mathrm{mn}}{\mathrm{\&Delta;}}_{\mathrm{mn}}(i,j)$

${\mathrm{\&Delta;}}_{\mathrm{mn}}(i,j)=\left\{\begin{array}{c}y(i+m,j+m)-y(i,j),|y(i+m,j+n)-y(i,j)|\&le;T\\ T,y(i+m,j+n)-y(i,j)>T\\ -T,y(i+m,j+n)-y(i,j)<-T\end{array}\right.$

其中t为滤波图像，y为校正图像，T为灰度差阈值。

步骤202：利用非线性平滑滤波器进行异常像素检测。具体操作如下：

$N_{\mathrm{mn}}(i,j)=\left\{\begin{array}{c}1,|y(i+m,j+n)-y(i,j)|\&le;T\\ 0,|y(i+m,j+n)-y(i,j)|>T\end{array}\right.$

$M(i,j)=\underset{m=-2}{\overset{2}{\mathrm{\&Sigma;}}}\underset{n=-2}{\overset{2}{\mathrm{\&Sigma;}}}{N}_{\mathrm{mn}}(i,j)$

若当前像素的有效像素计数值M(i,j)小于等于预设的阈值，即：

M(i,j)≤T_M

则判断(i,j)为异常像素，该位置的校正图像像素值、目标图像像素值均用最邻近已处理像素的值替换。替换规则如下：

$y(i,j)=\left\{\begin{array}{c}y(i,j-1),j\&NotEqual;1\\ y(i-1,j),j=1\end{array}\right.$

$t(i,j)=\left\{\begin{array}{c}y(i,j-1),j\&NotEqual;1\\ y(i-1,j),j=1\end{array}\right.$

其中y(i,j)为像素位置(i,j)的校正图像，y(i,j-1)为像素位置(i,j-1)的校正图像，y(i-1,j)为像素位置(i-1,j)的校正图像，t(i,j)为像素位置(i,j)的目标图像，i为行序号，j为列序号。

步骤203：利用非线性平滑滤波器计算当前像素的迭代步长。计算公式如下：

$\mathrm{\&eta;}(i,j)=\frac{k\max }{1+256\&times;\mathrm{ratio}(i,j)},$

ratio(i,j)＝1-M(i,j)/25

其中kmax为迭代步长可取的最大值。通过分析可以看出，在图像中的平坦区域，ratio(i,j)接近0，迭代步长η(i,j)接近kmax；在边缘区域，ratio(i,j)接近1，迭代步长较小。因此η(i,j)可随图像纹理变化自适应变化。

为了评估本发明的性能，将其与几种基于LMS的代表性非均匀性校正方法进行比较。参与比较的方法有：采用4邻域均值滤波的Scribner方法(Scribner2x2)，采用21x21均值滤波的Scribner方法(Scribner21x21)，自适应学习率方法(Vera2005)，门控的自适应学习率方法(Hardie2009)，使用双边滤波器的方法(Rossi2009)。

图3给出不同方法对模拟数据进行非均匀性校正的收敛性能比较。以校正图像序列与参考图像序列的绝对误差均值(MAE)作为性能指标，可以看出本发明的误差收敛速度显著快于其他方法。

图4给出不同方法对模拟数据进行非均匀性校正的结果。其中图4(a)与图4(b)分别为参考图像以及加入模拟非均匀性后的图像。图4(c)至图4(h)分别为Scribner2x2、Scribner21x21、Vera2005、Hardie2009、Rossi2009以及本发明的校正结果。可以看出，其他方法均存在不同程度的“鬼影”现象，只有本发明有效地抑制了“鬼影”。

图5给出不同方法对实际红外图像序列进行非均匀性校正的收敛性能比较。以校正图像序列的校正参数绝对均值(MA)作为性能指标，可以看出本发明的参数估计收敛速度显著快于其他方法。

图6给出不同方法对实际红外数据进行非均匀性校正的结果。其中图6(a)为校正前的图像。图6(b)至图6(g)分别为Scribner2x2、Scribner21x21、Vera2005、Hardie2009、Rossi2009以及本发明的校正结果。可以看出，其他方法均存在不同程度的“鬼影”，只有本发明有效地抑制了“鬼影”。

图7给出冲激噪声对不同方法的影响。其中图7(a)为被1％冲激噪声污染的原始图像，图7(b)为其他方法对该图像的校正结果，图7(c)为本发明的校正结果。可以看出，其他方法受冲激噪声影响很大，而本发明能有效消除冲激噪声。

由以上实验结果可以看出，本发明的校正参数收敛速度明显快于其他代表性算法，“鬼影”现象抑制能力明显强于其他代表性算法。此外，本发明具有异常像素检测与替换能力。更重要的是，本发明只需要进行简单的算术运算，因此计算复杂度低，适合硬件实时实现。

本领域的技术人员容易理解，以上所述仅为本发明的较佳实施例而已，并不用以限制本发明，凡在本发明的精神和原则之内所作的任何修改、等同替换和改进等，均应包含在本发明的保护范围之内。

Claims (2)

1.一种线性与非线性滤波器相结合的非均匀性校正方法，其特征在于，包括以下步骤：

步骤101：逐帧接收红外焦平面输出图像序列，对每一帧当前图像，依据图像中每个像素位置的当前增益、偏移校正参数值进行线性校正，得到校正图像，校正公式如下：

y_ij(n)＝g_ij(n)x_ij(n)+o_ij(n)

其中y_ij(n)为校正图像，x_ij(n)为焦平面输出图像，i和j分别表示像素的行、列位置，n表示当前图像在红外焦平面输出图像序列中的序号，g_ij(n)、o_ij(n)分别为当前增益、偏移校正参数，如果当前图像为第一帧，且已对红外焦平面阵列进行过基于黑体的标定，则以标定得到的增益、偏移校正值作为g_ij(n)和o_ij(n)，否则取g_ij(n)＝1，o_ij(n)＝0，如果当前图像不是第一帧，则采用前一帧的校正参数更新值作为g_ij(n)和o_ij(n)；

步骤102：判断当前图像非均匀性水平是较高还是较低，如果当前图像非均匀性水平较高，则转入步骤103；如果当前图像非均匀性水平较低，则转入步骤104；具体判决规则如下：如果当前图像为第一帧，则默认非均匀性水平较高，如果当前图像不是第一帧，且当前图像非均匀性水平stuck<2，则非均匀性水平较高，如果当前图像不是第一帧，且当前图像非均匀性水平stuck＝＝2，则非均匀性水平较低；

步骤103：对校正图像进行模板为11x11的均值滤波，并以滤波后的图像作为真实场景的估计图像，以校正图像与估计图像的差值作为误差图像，并计算误差图像的绝对均值，用当前图像的误差图像绝对均值与前一帧图像的误差图像绝对均值之比衡量误差变化程度，如果该比值接近1，则将当前图像非均匀性水平stuck递增1，否则将stuck清零，具体公式如下：

e_ij(n)＝y_ij(n)-t_ij(n)

$\mathrm{ame}\left(n\right)=\frac{1}{\mathrm{MN}}\underset{i=1}{\overset{M}{\mathrm{\&Sigma;}}}\underset{j=1}{\overset{N}{\mathrm{\&Sigma;}}}\left|e_{\mathrm{ij}}\left(n\right)\right|$

ratio(n)＝ame(n)/ame(n-1)

如果1-α≤ratio(n)≤1+α，则stuck＝stuck+1，否则stuck＝0

其中y_ij(n)为校正图像，t_ij(n)为估计图像，e_ij(n)为误差图像，ame(n)为误差图像绝对均值，ratio(n)为当前误差图像绝对均值与前一帧误差图像绝对均值的比值，α为预先设置的常数，且0＜α＜0.1；

步骤104：对校正图像进行模板为5x5的非线性平滑滤波，以滤波输出图像作为对真实场景的估计图像，并将校正图像与估计图像相减，得到误差图像，具体公式为e_ij(n)＝y_ij(n)-t_ij(n)，其中y_ij(n)为校正图像，t_ij(n)为估计图像，e_ij(n)为误差图像，非线性滤波器同时实现异常像素检测与替换并计算校正参数的自适应迭代步长；

步骤105：将校正图像与变化参考图像相减后取绝对值，将所得差值图像的每个像素值与变化阈值比较，如果像素值大于变化阈值，则设置迭代步长为预先设定的固定迭代步长或步骤104中得到的自适应迭代步长，同时更新变化参考图像中对应像素位置的值；否则，设置迭代步长为0，同时保持变化参考图像中对应像素位置的值不变，具体公式如下：

${\mathrm{\&eta;}}_{\mathrm{ij}}^{\&prime;}\left(n\right)=\left\{\begin{array}{c}{\mathrm{\&eta;}}_{\mathrm{ij}}\left(n\right),|t_{\mathrm{ij}}\left(n\right)-z_{\mathrm{ij}}\left(n\right)|>T\\ 0,\mathrm{else}\end{array}\right.$

$z_{\mathrm{ij}}(n+1)=\left\{\begin{array}{c}{t}_{\mathrm{ij}}\left(n\right),|t_{\mathrm{ij}}\left(n\right)-z_{\mathrm{ij}}\left(n\right)|>T\\ z_{\mathrm{ij}}\left(n\right),\mathrm{else}\end{array}\right.$

其中η_ij(n)为11x11均值滤波器使用的固定迭代步长或5x5非线性平滑滤波的自适应迭代步长，η′_ij(n)为实际迭代步长，t_ij(n)为目标图像，z_ij(n)为变化参考图像，z_ij(n+1)为用于下一帧的变化参考图像，T为变化阈值，z_ij(n)的初始值设置为大于输入图像最大值；

步骤106：执行非均匀性校正参数的更新，具体公式如下：

g_ij(n+1)＝g_ij(n)-η′_ij(n)e_ij(n)·y_ij(n)

o_ij(n+1)＝o_ij(n)-η′_ij(n)e_ij(n)

其中g_ij(n+1)为更新后的增益校正参数，o_ij(n+1)为更新后的偏移校正参数，η′_ij(n)为步骤105得到的迭代步长，e_ij(n)为步骤103或104得到的误差图像，y_ij(n)为步骤101得到的校正图像；

步骤107：重复以上步骤101至106，直到红外焦平面输出图像的所有帧都处理完毕为止。

2.根据权利要求1所述的非均匀性校正方法，其特征在于，步骤104进一步包括以下子步骤：

步骤201：对校正图像进行非线性平滑滤波，具体为采用5×5滑动窗口按照光栅扫描顺序逐像素对校正图像进行滤波操作，具体公式如下：

$t(i,j)=y(i,j)+\underset{m=-2}{\overset{2}{\mathrm{\&Sigma;}}}\underset{n=-2}{\overset{2}{\mathrm{\&Sigma;}}}{w}_{\mathrm{mn}}{\mathrm{\&Delta;}}_{\mathrm{mn}}(i,j)$

${\mathrm{\&Delta;}}_{\mathrm{mn}}(i,j)=\left\{\begin{array}{c}y(i+m,j+n)-y(i,j),|y(i+m,j+n)-y(i,j)|\&le;T\\ T,y(i+m,j+n)-y(i,j)>T\\ -T,y(i+m,j+n)-y(i,j)<-T\end{array}\right.$

其中t为滤波图像，y为校正图像，T为灰度差阈值；

步骤202：利用非线性平滑滤波器进行异常像素检测，具体公式如下：

$N_{\mathrm{mn}}(i,j)=\left\{\begin{array}{c}1,|y(i+m,j+n)-y(i,j)|\&le;T\\ 0,|y(i+m,j+n)-y(i,j)|>T\end{array}\right.$

$M(i,j)=\underset{m=-2}{\overset{2}{\mathrm{\&Sigma;}}}\underset{n=-2}{\overset{2}{\mathrm{\&Sigma;}}}{N}_{\mathrm{mn}}(i,j)$

若当前像素的有效像素计数值M(i,j)小于等于预设的阈值，即M(i,j)≤T_M，则判断(i,j)为异常像素，该异常像素位置的校正图像像素值、目标图像像素值均用最邻近已处理像素的值替换，替换规则如下：

$y(i,j)=\left\{\begin{array}{c}y(i,j-1),j\&NotEqual;1\\ y(i-1,j),j=1\end{array}\right.$

$t(i,j)=\left\{\begin{array}{c}y(i,j-1),j\&NotEqual;1\\ y(i-1,j),j=1\end{array}\right.$

其中y(i,j)为像素位置(i,j)的校正图像，y(i,j-1)为像素位置(i,j-1)的校正图像，y(i-1,j)为像素位置(i-1,j)的校正图像，t(i,j)为像素位置(i,j)的目标图像；

步骤203：利用非线性平滑滤波器计算当前像素的迭代步长，具体公式如下：

$\mathrm{\&eta;}(i,j)=\frac{k\max }{1+256\&times;\mathrm{ratio}(i,j)},$

ratio(i,j)＝1-M(i,j)/25

其中kmax为迭代步长可取的最大值。