CN105184740A - 一种红外焦平面图像非均匀性条带校正方法 - Google Patents

Abstract

本发明提出了一种红外焦平面图像非均匀性条带校正方法，具体为：1)采集含条带噪声的红外焦平面图像的像素点位置；2)对含条带噪声的红外焦平面图像建立单方向总变分去条带能量泛函模型；3)引入辅助变量替换能量泛函模型中的不可微项,进而得到新的能量泛函；4)依据新的能量泛将求解问题拆解为三个子问题，对所述三个子问题进行交替迭代求解，输出得到校正后图像f。本方法利用了条带只会对垂直于其的梯度产生影响，而不影响沿其方向的梯度这一特性知识，能快速有效地去除单方向条带噪声的方法，具有兼顾条带去除和图像细节保存的优点，适应性强，计算复杂度低。

Description

一种红外焦平面图像非均匀性条带校正方法

技术领域

本发明属于红外图像处理领域，更具体地，涉及到一种红外焦平面图像非均匀性条带校正方法，适用于快速去除多传感元成像数据中的单方向条带噪声。

背景技术

由于红外传感器材料、工艺水平限制和外界环境的影响，红外焦平面阵列中各探测元间响应不一致，导致图像数据中出现非均匀性条带噪声，严重影响了后续的数据处理。条带噪声校正的困难在于有效地去除各种不同类型的条带噪声的同时较完整地保存图像原有的结构信息。条带校正算法可以大致分为三大类：

第一大类基于统计匹配的方法，以直方图匹配和矩匹配方法为代表。统计匹配的方法简单快速，易于实现。缺点是需要假设图像的行或列的统计特性一致，并且还需事先已知条带的位置，不符合实际情况，去条带效果差。

第二大类基于图像滤波技术，以低通滤波、功率滤波器为代表，优点在于只对条带进行滤波处理，而不处理非条带内容，简单易于实现。但，严重依赖条带检测的准确性，漏检的条带被忽略不被处理，错检则会将图像内容也滤除掉，导致信息丢失。实际情况中很难准确检测条带。

第三大类基于变分正则化方法，将条带校正问题看作估计不含条带图像的逆问题，建立包含数据项和正则化项的能量函数，通过最小化能量函数迭代求解得到不含条带图像。Bouali等人的文章《TowardoptimaldestripingofMODISdatausingaunidirectionalvariationalmodel》提出了一种遥感卫星图像条带校正的变分方法，其能量泛函只包含x方向的梯度域数据保真项TV_x(f)和y方向上的数据约束项αTV_y(f-g)的总变分能量泛函：

$\begin{array}{c}E\left(f\right)={\mathrm{TV}}_x\left(f\right)+\mathrm{\α}{\mathrm{TV}}_y(f-g)\\ =\underset{\mathrm{\Ω}}{\∫}\left|\frac{\∂\left(f\right)}{\∂x}\right|\mathrm{dxdy}+\mathrm{\α}\underset{\mathrm{\Ω}}{\∫}\left|\frac{\∂(f-g)}{\∂y}\right|\mathrm{dxdy}\end{array}$

其中TV_x和TV_y分别表示对图像x与y方向求总变分。第二项表示我们寻找一个和条带图像在竖直方向有几乎一样变化的图像，因为条带噪声几乎只对水平方向的图像梯度产生影响，而第一项表示只对条带图像在水平方向进行平滑，该模型有效的反映了条带的单方向性质，但是该模型结果图像与原始图像在灰度值上整体出现较大偏差。

发明内容

本发明的目的在于提出一种红外焦平面图像非均匀性条带校正方法，将条带校正问题表示为从条带图像中估计真实图像的逆问题，利用了条带只会对垂直于其的梯度产生影响，而不影响沿其方向的梯度这一特性知识，能快速有效地去除单方向条带噪声的方法，具有兼顾条带去除和图像细节保存的优点，适应性强，计算复杂度低。

本发明与现有技术相比，具有如下优点：

第一，兼顾条带去除和图像细节保存。本发明将条带校正问题表示为从条带图像中估计真实图像的逆问题，利用了条带只会对垂直于其的梯度产生影响，而不影响沿其方向的梯度这一特性知识。既能有效地去除条带噪声，又能较好的保存图像细节。

第二，计算复杂度低。本发明利用分裂Bregman方法进行数值优化求解，有效地规避了不可微项。同时采用Guass_Sidel迭代进行求解，具有更快的收敛速度。

第三，适应性强。本发明中提供算法参数调节接口，适应多种条带噪声。

附图说明

图1为本发明方法流程图；

图2为真实红外图像轻度条带去除效果，图2a为原始的噪声图像；图2b为非均匀性校正后的图像；图2c为原始噪声图像与非均匀性校正后图像差值；

图3为真实红外图像中度条带去除结果，图3a为原始的噪声图像；图3b为非均匀性校正后的图像；图3c为原始噪声图像与非均匀性校正后图像差值；

图4为真实红外图像重度条带去除效果，图4a为原始的噪声图像；图4b为非均匀性校正后的图像；图4c为原始噪声图像与非均匀性校正后图像差值。

具体实施方式

为了使本发明的目的、技术方案及优点更加清楚明白，以下结合附图及实施例，对本发明进行进一步详细说明。

本发明根据单方向条带噪声的特点，在提出包括灰度域数据保真项y方向的梯度域数据保真项βTV_y(f-g)和x方向上的数据约束项αTV_x(f)的能量泛函，通过分裂Bregman方法求解和离散化实现。

实现本发明的方法是：首先利用单方向条带噪声的性质提出合理的能量泛函模型

$\begin{array}{c}E\left(f\right)=\frac{1}{2}{\left|\right|g-f\left|\right|}_2^2+\mathrm{\α}{\mathrm{TV}}_x\left(f\right)+\mathrm{\β}{\mathrm{TV}}_y(f-g)\\ =\frac{1}{2}\underset{\mathrm{\Ω}}{\∫}{(g-f)}^2\mathrm{dxdy}+\mathrm{\α}\underset{\mathrm{\Ω}}{\∫}\left|\frac{\∂\left(f\right)}{\∂x}\right|\mathrm{dxdy}+\mathrm{\β}\underset{\mathrm{\Ω}}{\∫}\left|\frac{\∂(f-g)}{\∂y}\right|\mathrm{dxdy}\end{array}$

g(x,y)为探测器输出的条带图像，f(x,y)为原始无条带图像，E(f)表示关于待估计图像f的能量泛函；能量泛函包括灰度域数据保真项y方向梯度数据保真项βTV_y(f-g)和x方向梯度惩罚项αTV_x(f)；其中TV_x和TV_y分别表示对图像x与y方向求总变分；α，β是正则化参数，用于调节正则化强度。(x、y)表示图像中像素点的位置，Ω表示图像像素坐标集合。第二行为具体展开表达式。

分裂Bregman方法求解引入辅助变量d_x,d_y,b_x,b_y,建立新的能量泛函:

$E(f,\mathrm{dx},\mathrm{dy},b_x,b_y)=\frac{1}{2}{\left|\right|g-f\left|\right|}_2^2+\mathrm{\α}{d}_x|+\mathrm{\β}{d}_y|+\frac{{\mathrm{\λ}}_1}{2}{\left|\right|d_x-{\▿}_{x}f-b_x\left|\right|}_2^2\frac{{\mathrm{\λ}}_2}{2}{\left|\right|d_y-{\▿}_y(f-g)-b_y\left|\right|}_2^2$

d_x,d_y是替代项，d_x替代αTV_x(f)，d_y替代βTV_y(f-g)，b_x,b_y则是数据找回项，使 $\frac{{\mathrm{\λ}}_1}{2}{\left|\right|d_x-{\▿}_{x}f-b_x\left|\right|}_2^2\frac{{\mathrm{\λ}}_2}{2}\left|\right|d_y-{\▿}_y(f-g)-b_y\left|\right|$ 更快的接近于0，加快迭代收敛速度。α，β是正则化参数，调节正则化约束强度。λ₁，λ₂为惩罚参数，用于约束引入的辅助变量。

对单方向总变分去条带能量泛函模型进行迭代求解，输出得到校正后图像f。依据新的能量泛将求解问题拆解为三个子问题，对所述三个子问题进行交替迭代求解，输出得到校正后图像f；其中，

第一子问题为：固定d，b，求解f，得到

$E\left(f\right)=\frac{1}{2}{\left|\right|g-f\left|\right|}_2^2+\frac{{\mathrm{\λ}}_1}{2}{\left|\right|d_x-{\▿}_{x}f-b_x\left|\right|}_2^2\frac{{\mathrm{\λ}}_2}{2}{\left|\right|d_y-{\▿}_y(f-g)-b_y\left|\right|}_2^2$

第二子问题为：固定f，b，求解d，得到

$E\left(d_x\right)=\mathrm{\α}\left|d_x\right|+\frac{{\mathrm{\λ}}_1}{2}{\left|\right|d_x-{\▿}_{x}f-b_x\left|\right|}_2^2$

$E\left(d_y\right)=\mathrm{\β}\left|d_y\right|+\frac{{\mathrm{\λ}}_2}{2}{\left|\right|d_y-{\▿}_y(f-g)-b_y\left|\right|}_2^2$

第三子问题为：固定f，d，求解b，得到

$E\left(b_x\right)=\frac{{\mathrm{\λ}}_1}{2}{\left|\right|d_x-{\▿}_{x}f-b_x\left|\right|}_2^2$

$E\left(b_y\right)=\frac{{\mathrm{\λ}}_2}{2}{\left|\right|d_y-{\▿}_y(f-g)-b_y\left|\right|}_2^2$

为了方便硬件实现，采用离散化实现。如图1所示，本发明迭代求解的具体流程如下：

1)输入含条带噪声的红外焦平面图像g，图像大小为M×N；初始化迭代次数k＝1，图像f^k-1＝g，辅助变量为M×N大小的零矩阵。

2)将辅助变量代入能量泛函迭代目标函数求解得到校正后图像f^k；具体求解过程为：

固定b，d求解f的数值实现

$\begin{array}{c}{f}^k(i,j)=({\mathrm{\λ}}_1(f^k(i-1,j)+f^{k-1}(i+1,j)+(d_x^{k-1}-b_x^{k-1})(i,j)-(d_x^{k-1}-b_x^{k-1})(i+1,j))\\ +{\mathrm{\λ}}_2(f^k(i,j-1)+f^{k-1}(i,j+1)+(d_y^{k-1}-b_y^{k-1})(i,j)-(d_y^{k-1}-b_y^{k-1})(i,j+1)-g(i.j-1)\\ -g(i,j+1))+g(i,j)(1+2{\mathrm{\λ}}_2))/(1+2{\mathrm{\λ}}_1+2{\mathrm{\λ}}_2)\end{array}$

3)将校正后图像f^k代入第二子问题求解d和第三子问题求解b对应的迭代目标函数得到辅助变量

固定b，f求解d的数值实现

$\begin{array}{c}{d}_x^k=\mathrm{shrink}({\▿}_x{f}^k+b_x^{k-1},\frac{\mathrm{\α}}{{\mathrm{\λ}}_1})\\ =\max \left\{\right|{\▿}_x{f}^k+b_x^{k-1}|-\frac{\mathrm{\α}}{{\mathrm{\λ}}_1},0\}\frac{{\▿}_x{f}^k+b_x^{k-1}}{|{\▿}_x{f}^k+b_x^{k-1}|}\end{array}$

$\begin{array}{c}{d}_y^k=\mathrm{shrink}({\▿}_y(f^k-g)+b_y^{k-1},\frac{\mathrm{\β}}{{\mathrm{\λ}}_2})\\ =\max \left\{\right|{\▿}_y(f^k-g)+b_y^{k-1}|-\frac{\mathrm{\β}}{{\mathrm{\λ}}_2},0\}\frac{{\▿}_y(f^k-g)+b_y^{k-1}}{|{\▿}_y(f^k-g)+b_y^{k-1}|}\end{array}$

固定d，f求解b的数值实现

$b_x^k=b_x^{k-1}+({\▿}_x{f}^k-d_x^k)$

$b_y^k=b_y^{k-1}+({\▿}_y(f^k-g)-d_y^k)$

4)如果‖f^k+1-f^k‖₂/‖f^k+1‖₂＜ε或者迭代次数大于要求的迭代次数，ε为预定阈值，则输出f^k+1作为最终处理图像，结束；否则，进入步骤5)。

5)更新k＝k+1，返回步骤3)。

在本方法中引入了正则化参数α和β，在分裂Bregman优化求解过程中又引入惩罚参数λ₁和λ₂。本发明针对不同强度的条带示例列出多组与之对应的参数，其仅仅作为示例，不作为对本发明限制。迭代次数固定在40次，在轻条带下的α＝20,β＝350,λ₁＝20，λ₂＝90；在中度条带下的α＝80,β＝800,λ₁＝30，λ₂＝110；在重度条带下的α＝200,β＝800,λ₁＝15，λ₂＝130。

图2是真实红外噪声图像(轻度条带噪声)的条带去除结果。图2a为原始的噪声图像，含有轻度的条带噪声；图2b为非均匀性校正后的图像；图2c为原始噪声图像与非均匀性校正后图像差值(用于描述条带噪声去除情况)，α＝20,β＝350,λ₁＝20，λ₂＝90。这组图像清晰地展示了本发明对轻度噪声图像的条带具有良好的去除效果，即保留了图像本身的细节，又最大程度的去除了噪声。

图3是真实红外噪声图像(中度条带噪声)的条带去除结果。图3a为原始的噪声图像，含有中度的条带噪声；图3b为非均匀性校正后的图像；图3c为原始噪声图像与非均匀性校正后图像差值(用于描述条带噪声去除情况)，α＝80,β＝800,lambda1＝30，lambda2＝110。这组图像清晰地展示了本发明对中度噪声图像的条带具有良好的去除效果，即保留了图像本身的细节(图像中的汉字依然很清晰)，又最大程度的去除了噪声。

图4是真实红外噪声图像(重度条带噪声)的条带去除结果。图4a为原始的噪声图像，含有重度的条带噪声；图4b为非均匀性校正后的图像；图4c为原始噪声图像与非均匀性校正后图像差值(用于描述条带噪声去除情况)，α＝200,β＝800,lambda1＝15，lambda2＝130。这组图像清晰地展示了本发明对重度噪声图像的条带具有良好的去除效果，同时保留了图像本身的细节(树叶和建筑的边缘依然清晰)，条带噪声被很好的抑制。

通过图2、图3、图4对比，本专利发明的方法对不同强度的条带噪声红外焦平面图像都有很好的条带去除效果，对细节的保存也是很好的。

Claims (5)

1.一种红外焦平面图像非均匀性条带校正方法，其特征在于，该方法具体为：

1)采集含条带噪声的红外焦平面图像；

2)对含条带噪声的红外焦平面图像建立单方向总变分去条带能量泛函模型：

$\begin{array}{c}E\left(f\right)=\frac{1}{2}{\left|\right|g-f\left|\right|}_2^2+{\mathrm{\αTV}}_x\left(f\right)+{\mathrm{\βTV}}_x(f-g)\\ =\frac{1}{2}\underset{\mathrm{\Ω}}{\∫}{(g-f)}^2\mathrm{dxdy}+\mathrm{\α}\underset{\mathrm{\Ω}}{\∫}\left|\frac{\∂\left(f\right)}{\∂x}\right|\mathrm{dxdy}+\mathrm{\β}\underset{\mathrm{\Ω}}{\∫}\left|\frac{\∂(f-g)}{\∂y}\right|\mathrm{dxdy}\end{array}$

其中，E(f)表示关于校正后图像f的能量泛函，g为含条带噪声的红外焦平面图像；能量泛函包括灰度域数据保真项TV_x(f)和TV_y(f)分别表示对校正后图像f的x和y方向的求总变分；α，β是正则化参数，表示求偏导，(x、y)表示图像中像素点的位置，Ω表示图像像素坐标集合；

3)引入辅助变量d_x,d_y,b_x,b_y替换能量泛函模型中的不可微项,进而得到新的能量泛函：

$\begin{array}{c}E(f,d_x,d_y,b_x,b_y)=\frac{1}{2}{\left|\right|g-f\left|\right|}_2^2+\mathrm{\α}\left|d_x\right|+\mathrm{\β}\left|d_y\right|\\ +\frac{{\mathrm{\λ}}_1}{2}{\left|\right|d_x-{\▿}_{x}f-b_x\left|\right|}_2^2+\frac{{\mathrm{\λ}}_2}{2}{\left|\right|d_y-{\▿}_y(f-g)-b_y\left|\right|}_2^2\end{array}$

惩罚参数λ₁、λ₂为用于约束引入的辅助变量；

4)依据新的能量泛将求解问题拆解为三个子问题，对所述三个子问题进行交替迭代求解，输出得到校正后图像f；其中，

第一子问题为：固定d，b，求解f，得到：

$E\left(f\right)=\frac{1}{2}{\left|\right|g-f\left|\right|}_2^2+\frac{{\mathrm{\λ}}_1}{2}{\left|\right|d_x-{\▿}_{x}f-b_x\left|\right|}_2^2+\frac{{\mathrm{\λ}}_2}{2}{\left|\right|d_y-{\▿}_y(f-g)-b_y\left|\right|}_2^2$

第二子问题为：固定f，b，求解d，得到：

$E\left(d_x\right)=\mathrm{\α}\left|d_x\right|+\frac{{\mathrm{\λ}}_1}{2}{\left|\right|d_x-{\▿}_{x}f-b_x\left|\right|}_2^2$

$E\left(d_y\right)=\mathrm{\β}\left|d_y\right|+\frac{{\mathrm{\λ}}_2}{2}{\left|\right|d_y-{\▿}_y(f-g)-b_y\left|\right|}_2^2$

第三子问题为：固定f，d，求解b，得到：

$E\left(b_x\right)=\frac{{\mathrm{\λ}}_1}{2}{\left|\right|d_x-{\▿}_{x}f-b_x\left|\right|}_2^2$

$E\left(b_y\right)=\frac{{\mathrm{\λ}}_2}{2}{\left|\right|d_y-{\▿}_y(f-g)-b_y\left|\right|}_2^2.$

2.根据权利要求1所述的红外焦平面图像非均匀性条带校正方法，其特征在于，所述步骤3)的具体实施方式为：

1)初始化辅助变量为零矩阵和迭代次数k＝1；

2)将辅助变量代入第一子问题求解f对应的迭代目标函数，求解得到校正后图像f^k；

3)将校正后图像f^k代入第二子问题求解d和第三子问题求解b对应的迭代目标函数得到辅助变量

4)若或k达到预定迭代次数最大值，ε为预定阈值，则校正后图像f^k即为最终校正后图像f，结束；否则，进入步骤4)；

5)更新k＝k+1，返回步骤3)。

3.根据权利要求2所述的红外焦平面图像非均匀性条带校正方法，其特征在于，所述第一子问题对应后的迭代目标函数为

$\begin{array}{c}{f}^k(i,j)=\left({\mathrm{\λ}}_1(f^k(i-1,j)+f^{k-1}(i+1,j)+(d_x^{k-1}-b_x^{k-1})(i,j)-(d_x^{k-1}-b_x^{k-1})(i+1,j)\right)\\ +{\mathrm{\λ}}_2(f^k(i,j-1)+f^{k-1}(i,j+1)+(d_y^{k-1}-b_y^{k-1})(i,j)-(d_y^{k-1}-b_y^{k-1})(i,j+1)-g(i.j-1)\\ -g(i,j+1))+g(i,j)(1+{2\mathrm{\λ}}_2))/(1+{2\mathrm{\λ}}_1+{2\mathrm{\λ}}_2)\end{array}$

f^k(i,j)是第k次迭代(i,j)位置处的像素值；代表第k-1次迭代时的辅助变量；分别代表第k-1次迭代时的Bregman变量；g代表原始图像。

4.根据权利要求3所述的红外焦平面图像非均匀性条带校正方法，其特征在于，所述第二子问题对应后的迭代目标函数为：

$\begin{array}{c}{d}_x^k=\mathrm{shrink}({\▿}_x{f}^k+b_x^{k-1},\frac{\mathrm{\α}}{{\mathrm{\λ}}_1})\\ =\max \left\{\right|{\▿}_x{f}^k+b_x^{k-1}|-\frac{\mathrm{\α}}{{\mathrm{\λ}}_1},0\}\frac{{\▿}_x{f}^k+b_x^{k-1}}{|{\▿}_x{f}^k+b_x^{k-1}|}\end{array}$

$\begin{array}{c}{d}_y^k=\mathrm{shrink}({\▿}_y(f^k-g)+b_y^{k-1},\frac{\mathrm{\β}}{{\mathrm{\λ}}_2})\\ =\max \left\{\right|{\▿}_y(f^k-g)+b_y^{k-1}|-\frac{\mathrm{\β}}{{\mathrm{\λ}}_2},0\}\frac{{\▿}_y(f^k-g)+b_y^{k-1}}{|{\▿}_y(f^k-g)+b_y^{k-1}|}\end{array}.$

其中，是从到的线性算子，其中表示N²维的实数维空间，类似地,也是从到的线性算子。shrink是软阈值函数。 $\mathrm{shrink}(a,b)=\max \left\{\right|a|-b,0\}\frac{a}{\left|a\right|}.$

5.根据权利要求4所述的红外焦平面图像非均匀性条带校正方法，其特征在于，所述第三子问题对应的迭代目标函数为：

$b_x^k=b_x^{k-1}+({\▿}_x{f}^k-d_x^k)$

$b_y^k=b_y^{k-1}+({\▿}_y(f^k-g)-d_y^k).$