CN102779332A - 时域卡尔曼滤波修正的非线性拟合红外非均匀校正方法 - Google Patents

Abstract

本发明提供一种时域卡尔曼滤波修正的非线性拟合红外非均匀校正方法。技术方案是首先假设红外焦平面探测器各阵列单元的响应变化曲线在时间上是连续的，采用高阶多参数非线性多项式可对其进行拟合和描述；然后采集四个不同温度点各阵列单元的响应输出值，利用非线性方程求解方法确定各阵列单元响应表达式中的多个参数；为了解决探测器响应随时间漂移问题，利用时域Kalman滤波对响应表达式中各参数进行时漂修正，得到探测器阵列单元响应曲线的解析表示；最后利用该解析表示计算任意时间任意温度条件下探测器各阵列单元响应输出。本发明实现对探测器响应不一致的校正，解决红外图像非均匀性问题。

Description

时域卡尔曼滤波修正的非线性拟合红外非均匀校正方法

技术领域

本发明属于红外图像预处理技术领域，涉及一种对红外图像进行非均匀校正的方法。

背景技术

红外图像非均匀校正是图像处理领域里的一个重要研究方向。红外图像非均匀校正的主要目的是利用图像处理手段，解决由于材料缺陷、电路的稳定性以及集成工艺水平的限制，导致的探测器输出响应不一致的问题，使校正后的图像便于后续处理。目前，红外图像非均匀校正的主要方法可分为两类：基于参考源的温度定标校正法和基于场景的自适应校正法。这两类方法的基本原理和技术特点如下：

一、基于参考源的温度定标校正方法

基于参考源的温度定标校正方法主要包括单点温度定标法，两点温度定标法和多点温度定标法等。这类方法的设计思想是：利用参考辐射源给红外焦平面阵列提供均匀辐照度，对每个探测器单元的响应输出进行测量，由此计算得出个探测器单元的校正参数。

二、基于场景的自适应校正法

基于场景的自适应校正方法主要包括时域高通滤波法、神经网络校正法、常量统计法、线性滤波校正法和场景匹配法等。这类算法的基本原理是计算增益系数和偏移量的数据不是取自参考辐射源，而是全部或部分来自于场景的估计。

非均匀性是红外焦平面探测器的固有属性，在理想情况下，红外焦平面面阵探测器受均匀入射辐射时，其各个像元的信号输出值应该完全一致；但实际上在制作器件的半导体材料不均匀（杂质浓度、晶体缺陷、内部结构的不均匀性等）、器件工作状态、生产工艺过程以及外界输入等的综合影响下，其输出幅度并不相同，这就是所谓的红外图像非均匀性（Non-Uniformity，NU）。对于单点扫描方式的探测器来说不存在非均匀性问题，线阵扫描方式的探测器中的非均匀性存在于线阵方向，而红外焦平面探测器的非均匀性存在于整个焦平面上，愈是大规模的器件，非均匀性问题就愈突出。这种非均匀性会导致系统的温度分辨率下降，使目标图像的质量受到严重影响，从而限制了其在高灵敏度检测方面的应用。

现有的图像校正方法在解决红外图像非均匀性问题时存在不足，主要表现在以下几个方面：

基于参考源的温度定标校正方法：校正参数是固定不变的，而实际上，随着器件工作温度和环境温度的变化，器件的工作状态会发生变化，如果还采用原先计算出的校正参数进行校正，就会使校正效果变差；基于场景的自适应校正法：便于处理存在运动目标的红外图像，而且计算量大，实时处理需要先进的多处理器机构。

综上所述，针对红外图像非均匀性问题设计适于工程实现的、适用图像范围广的非均匀校正方法是一个急需解决的工程技术问题。目前尚未发现有关这个问题的公开研究资料。

发明内容

本发明的目的是提供一种红外图像非均匀校正方法，解决红外焦平面探测器输出响应不一致所导致的图像降质问题。

技术方案是首先假设红外焦平面探测器各阵列单元的响应变化曲线在时间上是连续的，采用高阶多参数非线性多项式可对其进行拟合和描述；然后采集四个不同温度点各阵列单元的响应输出值，利用非线性方程求解方法确定各阵列单元响应表达式中的多个参数；第三步为了解决探测器响应随时间漂移问题，利用时域Kalman滤波对响应表达式中各参数进行时漂修正，得到探测器阵列单元响应曲线的解析表示；最后利用该解析表示计算任意时间任意温度条件下探测器各阵列单元响应输出，从而实现对探测器响应不一致的校正，解决红外图像非均匀性问题。

本发明的技术方案包括以下步骤：

第一步：参数估计

在红外探测器常规工作温度10℃到70℃之间，任意选择未经校正的四个不同温度点T₁,T₂,T₃,T₄的红外图像，f_i,j(T₁),f_i,j(T₂),f_i,j(T₃),f_i,j(T₄)分别为这四幅红外图像像素点的灰度值，其中下标i,j表示像素点行和列，将其分别代入公式一，求解探测器响应方程：

$\left\{\begin{array}{c}\stackrel{\‾}{F}\left(T_1\right)=A_{i,j}\·{\left(f_{i,j}\left(T_1\right)\right)}^3+B_{i,j}\·{\left(f_{i,j}\left(T_1\right)\right)}^2+C_{i,j}\·\left(f_{i,j}\left(T_1\right)\right)+D_{\mathrm{ij}}\\ \stackrel{\‾}{F}\left(T_2\right)=A_{i,j}\·{\left(f_{i,j}\left(T_2\right)\right)}^3+B_{i,j}\·{\left(f_{i,j}\left(T_2\right)\right)}^2+C_{i,j}\·\left(f_{i,j}\left(T_2\right)\right)+D_{\mathrm{ij}}\\ \stackrel{\‾}{F}\left(T_3\right)=A_{i,j}\·{\left(f_{i,j}\left(T_3\right)\right)}^3+B_{i,j}\·{\left(f_{i,j}\left(T_3\right)\right)}^2+C_{i,j}\·\left(f_{i,j}\left(T_3\right)\right)+D_{\mathrm{ij}}\\ \stackrel{\‾}{F}\left(T_4\right)=A_{i,j}\·{\left(f_{i,j}\left(T_4\right)\right)}^3+B_{i,j}\·{\left(f_{i,j}\left(T_4\right)\right)}^2+C_{i,j}\left(f_{i,j}\left(T_4\right)\right)+D_{\mathrm{ij}}\end{array}\right.$ （公式一）

其中，

n＝1,2,3,4为上述四幅红外图像灰度均值，A_i，j，B_i，j,C_i，j，D_i,j为待估计的非线性曲线各阶描述参数。利用上述方程，求解出A_i，j,B_i，j，C_i,j，D_i，j的估计值，记为

第二步：实现校正

将任意温度点T下需校正的红外图像像素点的灰度值g_i,j(T)代入公式二，得到校正后红外图像像素点的灰度值

计算公式如下：

${\hat{g}}_{i,j}\left(T\right)={\hat{A}}_{i,j}\·{\left(g_{i,j}\left(T\right)\right)}^3+\·{\hat{B}}_{i,j}{\left(g_{i,j}\left(T\right)\right)}^2+{\hat{C}}_{i,j}\·\left(g_{i,j}\left(T\right)\right)+{\hat{D}}_{\mathrm{ij}}$ （公式二）

第三步：时域Kalman滤波修正

校正后的红外探测器在工作一段时间后，会出现焦平面阵列输出响应随工作时长产生波动的现象，即发生非均匀性时漂，时漂的出现严重影响探测器的工作性能，因此需要对校正后红外图像像素点的灰度值

进行时漂处理，对探测器响应曲线参数进行修正，本发明利用时域Kalman滤波修正方式解决时漂问题。具体实现如下：

第1步：依据图像NU(T)（Non-uniformity，非均匀度）值确定是否进行时漂处理，并确定α_k和β_k两个时漂因子值，NU(T)值计算方法如公式三所示：

$\mathrm{NU}\left(T\right)=\frac{\underset{i,j}{\mathrm{\Σ}}|{\hat{g}}_{i,j}\left(T\right)-\stackrel{\‾}{G}\left(T\right)|}{N\stackrel{\‾}{G}\left(T\right)}$ （公式三）

其中：

为红外图像灰度值

的均值，N为该红外图像像素数目， $\stackrel{\‾}{G}\left(T\right)=\frac{1}{N}\underset{i,j}{\mathrm{\Σ}}{\hat{g}}_{i,j}\left(T\right).$

当NU(T)≤3.0‰时不需要进行时漂处理，则技术方案结束；当NU(T)＞3.0‰，取α_k=α∈[0.99,1]，β_k=β∈[0.99,1]，并进行时漂处理；当NU(T)＞4.0‰，则令α_k=α∈[0.9,0.99]，β_k=β∈[0.9,0.99]，并进行时漂处理；如果时漂更大，即当NU(T)＞4.5‰，令α_k=α∈[0.8,0.9]，β_k=β∈[0.8,0.9]，并进行时漂处理。

时漂处理是指下述第2步至第4步：

第2步：建立Kalman滤波的状态方程和观测方程：

状态方程：X_i,j(k+1)＝Φ_kX_i,j(k)+M_k+W_k    （公式四）

观测方程：Y_i,j(k)＝H_kX_i,j(k)+V_k    （公式五）

其中，k表示当前时刻，k+1表示下一时刻，k＝0,1,2,…，并且k＝0对应于得到需要进行时漂处理的红外图像的时刻T；状态向量X_i,j(k)定义为

分别表示k时刻非线性曲线描述参数，

${\mathrm{\Φ}}_k=\left[\begin{array}{cc}{\mathrm{\α}}_k& 0\\ 0& {\mathrm{\β}}_k\end{array}\right]$ 为状态转移矩阵，驱动噪声均值定义为 $M_k=\left[\begin{array}{cc}1-{\mathrm{\α}}_k& 0\\ 0& 1-{\mathrm{\β}}_k\end{array}\right]{\hat{X}}_{i,j}\left(0\right),$ 其中

Y_i,j(k)表示k时刻观测到的红外图像像素点的灰度值，

为观测矩阵；W_k和V_k分别为噪声干扰，其协方差Q_k和R_k分别为；

$Q_k=\left[\begin{array}{cc}(1-{\mathrm{\α}}_k^2){\mathrm{\σ}}_{\mathrm{\α}0}^2& 0\\ 0& (1-{\mathrm{\β}}_k^2){\mathrm{\σ}}_{\mathrm{\β}0}^2\end{array}\right],$ $R_k={\mathrm{I\σ}}_{\mathrm{vk}}^2.$ (公式六)

根据实践经验取 ${\mathrm{\σ}}_{\mathrm{\α}0}^2\∈\left[\mathrm{0.2,0.25}\right],$ ${\mathrm{\σ}}_{\mathrm{\β}0}^2\∈\left[\mathrm{0.05,0.1}\right],$ ${\mathrm{\σ}}_{\mathrm{vk}}^2\∈\left[\mathrm{0.05,0.1}\right].$

第3步：将α_k和β_k的取值代入Kalman滤波的状态方程和观测方程，迭代得到最终的滤波估计值

即得到

为经过时域卡尔曼滤波修正后的非线性曲线描述参数。

第4步：实现时漂修正

将

代入公式七中，对

进行时漂修正，得到经过时漂处理的红外图像像素点的灰度值

${\hat{\hat{g}}}_{i,j}\left(T\right)=(1-{\mathrm{\α}}_k)\·{\hat{\hat{A}}}_{i,j}\left(k\right)\·\stackrel{\‾}{G}\left(T\right)+(1-{\mathrm{\β}}_k)\·{\hat{\hat{B}}}_{i,j}\left(k\right)$ (公式七)

采用本发明可以取得以下技术效果：

本发明能够稳定可靠地实现红外图像非均匀性校正，明显改善由于探测器各像元响应的不一致所导致的图像非均匀性噪声，提高系统温度分辨率、提升图像预处理质量，为后续的图像检测、识别、跟踪提供良好数据源。本发明所提出的红外图像非均匀校正方法，具有以下比较明显的特点和优势：

1.根据本发明的实验结果以及和其它常用红外图像非均匀校正方法的性能对比结果可见：本发明能够有效地克服非均匀噪声的干扰，具有简单易行、准确率高、鲁棒性强、易于FPGA（Field Programmable Gate Array，现场可编程门阵列）；

2.依据本发明第一步和第二步，可很好地表征红外探测器输出响应随温度的非线性变化规律，解决了传统校正方法难于处理接近温度饱和区域的红外图像校正问题；

3.依据本发明第三步“时域Kalman滤波修正”，可动态完成探测器响应时漂修正，解决了校正后的红外探测器在工作一段时间后，阵列输出响应随工作时长产生波动的问题，极大提高了校正方法的场景适应性。

附图说明

图1是本发明的原理流程图；

图2是20℃温度点时利用本发明进行的仿真实验结果；

图3是50℃温度点时利用本发明进行的仿真实验结果；

图4是利用本发明和两点法进行对比实验的结果；

图5是利用本发明和稳态卡尔曼滤波法进行对比实验的结果。

具体实施方式

下面结合附图对本发明做进一步说明。

图2至图5进行实验的结果，横坐标均为图像帧号，纵坐标为不同的图像非均匀性描述参数。实验中利用本发明的第一步进行参数估计时，选择四个不同温度点T₁,T₂,T₃,T₄分别是15℃、25℃、40℃、55℃。对温度点T下需校正的红外图像进行处理，温度T的取值分别是20℃、50℃和40℃。

图2（a）为20℃温度点红外图像序列灰度均方差变化曲线图，左图为校正前结果，右图为校正后结果；图2（b）为20℃温度点红外图像序列非均匀度变化曲线图，左图为校正前结果，右图为校正后结果；从实验结果可以看出，本发明对20℃温度点红外图像的非均匀性具有明显的校正效果，图像灰度均方差由校正前的171左右降低到校正后的1.4左右，图像非均匀度由校正前的9.5‰左右，降低到校正后的0.14‰左右，结果均证明了本发明的有效性。

图3（a）为50℃温度点红外图像序列灰度均方差变化曲线图，左图为校正前结果，右图为校正后结果；图3（b）为50℃温度点红外图像序列非均匀度变化曲线图，左图为校正前结果，右图为校正后结果；从实验结果可以看出，本发明对50℃温度点红外图像的非均匀性具有明显的校正效果，图像灰度均方差由校正前的384-385降低到校正后的1.3-2.0，图像非均匀度由校正前的13.08‰-13.10‰，降低到校正后的0.6‰-0.10‰，结果均证明了本发明的有效性。

图4为针对50℃温度点红外图像序列，本发明与两点定标校正法的校正效果对比图，本发明对应曲线用“----”表示，两点定标校正法对应曲线用“——”表示；图4（a）为灰度均方差变化曲线图，左图为校正前对比结果，右图为校正后对比结果；图4（b）为非均匀度变化曲线图，左图为校正前对比结果，右图为校正后对比结果；从实验结果可以看出，与在工程实践过程中广泛用到的两点定标校正法比较，本发明的校正效果明显优于两点定标校正法，图像灰度均方差和非均匀度两个图像非均匀性描述参数明显小于两点定标法相应的参数指标，特别是在接近温度饱和的情况下，由于本发明采用了非线性拟合方法实现对探测器响应参数的估计，避免了线性拟合方法带来的失真问题，保证了校正效果的鲁棒性。

图5为针对40℃温度点红外图像序列，本发明与稳态卡尔曼滤波校正方法的校正效果对比图，本发明对应曲线用“----”表示，稳态卡尔曼滤波校正方法对应曲线用“——”表示；图5（a）为灰度均方差变化曲线图，左图为校正前对比结果，右图为校正后对比结果；图5（b）为非均匀度变化曲线图，左图为校正前对比结果，右图为校正后对比结果；从实验结果可以看出，与典型的基于场景的红外非均匀校正方法—稳态卡尔曼滤波校正方法相比，利用图像灰度均方差和非均匀度两个评价指标进行非均匀性校正效果评价，由于本发明充分利用了实时成像特征数据，因此在校正性能明显优于稳态卡尔曼滤波法。

Claims (1)

1.一种时域卡尔曼滤波修正的非线性拟合红外非均匀校正方法，其特征在于，包括下述步骤：

第一步：参数估计：

在10℃到70℃之间，任意选择未经校正的四个不同温度点T₁,T₂,T₃,T₄的红外图像，f_i,j(T₁),f_i,j(T₂),f_i,j(T₃),f_i,j(T₄)分别为这四幅红外图像像素点的灰度值，其中下标i，j表示像素点行和列，将其分别代入公式一，求解探测器响应方程：

$\left\{\begin{array}{c}\stackrel{\‾}{F}\left(T_1\right)=A_{i,j}\·{\left(f_{i,j}\left(T_1\right)\right)}^3+B_{i,j}\·{\left(f_{i,j}\left(T_1\right)\right)}^2+C_{i,j}\·\left(f_{i,j}\left(T_1\right)\right)+D_{\mathrm{ij}}\\ \stackrel{\‾}{F}\left(T_2\right)=A_{i,j}\·{\left(f_{i,j}\left(T_2\right)\right)}^3+B_{i,j}\·{\left(f_{i,j}\left(T_2\right)\right)}^2+C_{i,j}\·\left(f_{i,j}\left(T_2\right)\right)+D_{\mathrm{ij}}\\ \stackrel{\‾}{F}\left(T_3\right)=A_{i,j}\·{\left(f_{i,j}\left(T_3\right)\right)}^3+B_{i,j}\·{\left(f_{i,j}\left(T_3\right)\right)}^2+C_{i,j}\·\left(f_{i,j}\left(T_3\right)\right)+D_{\mathrm{ij}}\\ \stackrel{\‾}{F}\left(T_4\right)=A_{i,j}\·{\left(f_{i,j}\left(T_4\right)\right)}^3+B_{i,j}\·{\left(f_{i,j}\left(T_4\right)\right)}^2+C_{i,j}\left(f_{i,j}\left(T_4\right)\right)+D_{\mathrm{ij}}\end{array}\right.$ （公式一）

上式中，

n＝1,2,3,4为上述四幅红外图像灰度均值，A_i，j，B_i，j，C_i，j，D_i,j为待估计的非线性曲线各阶描述参数；利用公式一，求解出A_i，j，B_i，j，C_i，j，D_i,j的估计值，记为

第二步：实现校正：

将任意温度点T下需校正的红外图像像素点的灰度值g_i,j(T)代入公式二，得到校正后红外图像像素点的灰度值

${\hat{g}}_{i,j}\left(T\right)={\hat{A}}_{i,j}\·{\left(g_{i,j}\left(T\right)\right)}^3+\·{\hat{B}}_{i,j}{\left(g_{i,j}\left(T\right)\right)}^2+{\hat{C}}_{i,j}\·\left(g_{i,j}\left(T\right)\right)+{\hat{D}}_{\mathrm{ij}}$ （公式二）

第三步：时域Kalman滤波修正：

第1步：计算图像NU(T)：

$\mathrm{NU}\left(T\right)=\frac{\underset{i,j}{\mathrm{\Σ}}|{\hat{g}}_{i,j}\left(T\right)-\stackrel{\‾}{G}\left(T\right)|}{N\stackrel{\‾}{G}\left(T\right)}$ （公式三）

其中：

为红外图像灰度值

的均值，N为该红外图像像素数目， $\stackrel{\‾}{G}\left(T\right)=\frac{1}{N}\underset{i,j}{\mathrm{\Σ}}{\hat{g}}_{i,j}\left(T\right);$

当NU(T)≤3.0‰时不需要进行时漂处理，则技术方案结束；当NU(T)＞3.0‰，取α_k=α∈[0.99,1]，β_k=β∈[0.99,1]，并进行时漂处理；当NU(T)＞4.0‰，则令α_k=α∈[0.9,0.99]，β_k=β∈[0.9,0.99]，并进行时漂处理；如果时漂更大，即当NU(T)＞4.5‰，令α_k=α∈[0.8,0.9]，β_k=β∈[0.8,0.9]，并进行时漂处理；

时漂处理是指下述第2步至第4步：

第2步：建立Kalman滤波的状态方程和观测方程：

状态方程：X_i,j(k+1)＝Φ_kX_i,j(k)+M_k+W_k    （公式四）

观测方程：Y_i,j(k)＝H_kX_i,j(k)+V_k    （公式五）

其中，k表示当前时刻，k+1表示下一时刻，k＝0,1,2,…，并且k＝0对应于得到需要进行时漂处理的红外图像的时刻T；状态向量X_i，j(k)定义为

分别表示k时刻非线性曲线描述参数，

${\mathrm{\Φ}}_k=\left[\begin{array}{cc}{\mathrm{\α}}_k& 0\\ 0& {\mathrm{\β}}_k\end{array}\right]$ 为状态转移矩阵，驱动噪声均值定义为 $M_k=\left[\begin{array}{cc}1-{\mathrm{\α}}_k& 0\\ 0& 1-{\mathrm{\β}}_k\end{array}\right]{\hat{X}}_{i,j}\left(0\right),$ 其中

Y_i,j(k)表示k时刻观测到的红外图像像素点的灰度值，

为观测矩阵；W_k和V_k分别为噪声干扰，其协方差Q_k和R_k分别为；

$Q_k=\left[\begin{array}{cc}(1-{\mathrm{\α}}_k^2){\mathrm{\σ}}_{\mathrm{\α}0}^2& 0\\ 0& (1-{\mathrm{\β}}_k^2){\mathrm{\σ}}_{\mathrm{\β}0}^2\end{array}\right],$ $R_k={\mathrm{I\σ}}_{\mathrm{vk}}^2;$ (公式六)

上式中，取 ${\mathrm{\σ}}_{\mathrm{\α}0}^2\∈\left[\mathrm{0.2,0.25}\right],$ ${\mathrm{\σ}}_{\mathrm{\β}0}^2\∈\left[\mathrm{0.05,0.1}\right],$ ${\mathrm{\σ}}_{\mathrm{vk}}^2\∈\left[\mathrm{0.05,0.1}\right];$

第3步：将α_k和β_k的取值代入Kalman滤波的状态方程和观测方程，迭代得到最终的滤波估计值

即得到

为经过时域卡尔曼滤波修正后的非线性曲线描述参数；

第4步：实现时漂修正：

将代入公式七中，对进行时漂修正，得到经过时漂处理的红外图像像素点的灰度值

${\hat{\hat{g}}}_{i,j}\left(T\right)=(1-{\mathrm{\α}}_k)\·{\hat{\hat{A}}}_{i,j}\left(k\right)\·\stackrel{\‾}{G}\left(T\right)+(1-{\mathrm{\β}}_k)\·{\hat{\hat{B}}}_{i,j}\left(k\right)$ (公式七)。

Images