基于各向异性时空域全变分非均匀性校正的方法与流程

本发明涉及红外视频图像的非均匀校正技术，具体涉及一种基于各向异性时空域全变分非均匀性校正的方法。

背景技术：

红外焦平面阵列器件(IRFPA)是当前最主流的红外探测器，广泛应用于各个领域，但由于制造工艺以及工作环境等因素影响，IRFPA的探测元的响应度不一致，具体表现为输出图像存在一定的非均匀固定条纹噪声(NUC-FPN)，严重影响成像质量。

目前，国内外的非均匀校正(NUC)方法大致可分为两大类：基于定标的方法和基于场景的方法。其中基于定标的方法利用温度分布均匀的红外热源(诸如黑体)对探测器进行标定来解决非均匀性的问题，该类方法的准确度比较高，但是在发生温漂时需要重新标定来不断的修正校正参数，操作十分麻烦。基于场景的方法主要有基于统计的和基于配准的两大类。基于统计的方法通常对焦平面接收到的辐射量进行时间上以及空间上的统计假设，在此假设的基础上不断修正校正参数，完成非均匀性校正。其中代表性的有时域高通法、统计恒定法、神经网络法、恒定范围法及其相应的扩展形式。然而，由于图像场景千变万化，上述算法的假设并不能在任何场合均得到满足，因此这类校正算法通常伴随较为严重的鬼影。另一类是基于配准的方法，其假设前提为，在较短的时间间隔内，对相同的图像场景，如果不存在非均匀性，则每个像元的响应应该是相同的，因此这类技术需要对红外图像序列进行精准的配准。其中比较有代表性的有全景图积累法，代数校正法等。这类算法存在计算量与存储量大，校正误差易逐级累积并传播，无法处理旋转和缩放情况下相邻帧间配准等问题。

基于全变分红外图像非均匀性校正方法最初由Esteban Vera和Pablo Meza提出(V.Esteban,M.Pablo,“Total variation approach for adaptive nonuniformity correction in focal-plane arrays,”Optics Letters,36,172-174(2011))，通过将非均匀噪声视为固定噪声的一种，从单纯的图像角度进行去噪处理。具有实时性好，适用范围广，校正精度高等优点。缺点是论文中采用的变分模型是同时沿着x轴和y轴进行的，忽视了非均匀噪声在结构上表现为条纹噪声的特点。同时文中采用了L2范数进行正则约束，在滤去噪声的同时，也模糊了图像的细节，也没有考虑到视频图像帧间的相关性。因此，本专利在该篇论文的基础上提出了一种基于各向异性时空域全变分非均匀性校正的新方法，能在校正非均匀性条纹噪声的同时最大限度的保护图像细节。

技术实现要素：

本发明的目的在于提供一种基于各向异性时空域全变分非均匀性校正的方法，可以在无需预先标定的情况下实时的进行非均匀性条纹噪声的去噪，在有效去除图像的固定条纹噪声的同时，最大程度的保留图像的细节信息，提高图像的质量。该方法具有除噪效果好，处理速度快等优点，能满足了视频处理的实时性要求。

实现本发明目的的技术解决方案为：一种基于各向异性时空域全变分非均匀性校正的方法，包括以下步骤：

步骤1、采集含非均匀性噪声的红外图像序列f(i,j,t)＝u(i,j,t)+n(i,j,t)，其中f(i,j,t)表示实际观察到的图像，u(i,j,t)表示理想的无噪图像，n(i,j,t)表示由红外探测器成像单元的暗电流和响应率不同引起的固定条纹噪声，(i,j)表示图像像素点的位置，t表示视频序列的第t帧图像。

步骤2、采用全变分模型对上述含非均匀性噪声的红外图像序列f(i,j,t)建立图像的全变分模型J(u)，并对图像的全变分模型J(u)进行更新：

其中f为f(i,j,t)的通用表示，u为u(i,j,t)的通用表示，H(u)为正则项；J(u)为图像能量泛函；为保真项；λ为平滑系数；Ω表示单帧图像像素点构成的空间。

又对正则项H(u)取H(u)＝λ₁||u_x||₁+λ₂||u_t||₁ (2)

d_x为图像时空域上对于像素点的x轴的一阶差分，d_t为图像时空域上对于像素点的帧间的一阶差分；λ₁为图像时空域上对于像素点的x轴的一阶差分的平滑系数，λ₂为时空域上对于像素点的帧间的一阶差分的平滑系数；其具体的表达式如下：

u_x(i,j,t)＝u(i+1,j,t)-u(i,j,t) (3)

u_t(i,j,t)＝u(i,j,t)-u(i,j,t+1) (4)

将式(2)～(4)带入式(1)，更新图像的全变分模型J(u)：

步骤3：通过Split Bregman方法，对式(5)进行迭代求最小化，获得最终的无噪图像u_final。

上述步骤2中，采用全变分模型对上述含非均匀性噪声的红外图像序列f(i,j,t)建立图像的全变分模型，并对全变分模型进行更新，具体的步骤为：

步骤2-1)输入视频图像序列f(i,j,t)＝u(i,j,t)+n(i,j,t)，读取相邻两帧图像f_t1(M,N,t₁)和f_t2(M,N,t₂)，其中M为图像的总行数，N为图像的总列数，t₁和t₂均代表帧数。

步骤2-2)对于t₁时刻的图像，对于像素点(i,j)建立图像噪声模型f_t1(i,j,t₁)＝u(i,j,t₁)+n(i,j,t₁)，其时空域x轴方向和帧间方向上t的一阶差分表示如下：

对于x轴方向，当j＜M时，取u_x(i,j,t₁)＝u(i,j+1,t₁)-u(i,j,t₁)_；

对帧间方向：u_t(i,j,t)＝u(i,j,t₂)-u(i,j,t₁)；

对于图像边界处j＝M，取u_x(i,M,t₁)＝0；u_x(i,M,t₂)＝0；

步骤2-3)将式(2)～(4)带入式(1)，得到全变分模型J(u)离散化表示：

设置帧间阈值T，判断f(i,j,t₂)-f(i,j,t₁)是否大于T：

当f(i,j,t₂)-f(i,j,t₁)≤T，保留式(5)中的u_t(i,j,t₁)项；

当f(i,j,t₂)-f(i,j,t₁)>T，删除式(5)中的u_t(i,j,t)项，此时J(u)退化为一个仅对图像空间域竖直方向进行校正的单帧模型：

。

进一步的，上述步骤3中，通过Split Bregman方法，对式(5)进行迭代求最小化，获得最终的无噪图像u_final，具体的步骤为：

3-1)加入二次惩罚项，通过拉格朗日乘子法，将式(5)中转化为一个无约束最小化问题，得到式(7)；

其中d_x是图像水平方向的二次惩罚项，d_t是时间方向的二次惩罚项，b_x是引入图像水平方向的辅助变量，b_t是引入的时间方向的辅助变量，α和β是方程的Bregman惩罚系数；

3-2)将式(7)中u，d_x，d_y这三个未知变量进行分离，得到两个最优化子问题；

对于变量u，分离后最优化子方程如下：

对于变量d_x和d_t，分离后最优化子方程组如下：

3-3)式(8)等价于一个基于线性算子的迭代过程，即式(10)：

通过快速傅里叶变换，得到u的封闭解如下：

其中F是快速傅里叶变换，F^-1为傅里叶逆变换；

3-4)采用shrink算子对式(9)的迭代最小化，得到；

其中λ为收缩系数，取0.05，收缩算子shrink，对于d_x，当取当取同理，对d_t也采用如上处理方式；

3-5)对于余下的未知数b_x和b_y，采用代入法迭代求解：

3-6)使用Split Bregman进行最优化时，采用改变一个变量并固定其他变量的策略交替迭代，当迭代到第n次时，满足uⁿ⁺¹-uⁿ＜10^-3，此时迭代停止，得到的变量u即为最终的无噪图像u_final。

本发明与现有技术相比，其显著优点为：1、利用视频图像序列中相邻帧图像之间的相关性扩展了全变差去噪模型，相比较单纯的空域模型能有更好的平滑条纹噪声的效果。

2、利用红外非均匀性固定条纹噪声的几何性质，在图像行方向最小化变分，在列方向则尽量保存图像信息。

3、引入Split Bregman法来求解全变分方程，解决了传统的最陡梯度法速度过慢的缺陷，将整个模型的最优化过程分离为数个子最优化过程，大大提升了速度，使其满足红外视频图像的非均匀校正实时性要求。

4、引入帧间阈值T来防止红外视频序列帧间差别过大的情况，提高模型的精确性。

附图说明

图1是本发明一种基于各向异性时空域全变分模型的非均匀校正方法算法流程图。

图2是视频序列f(i,j,t)恒满足f(i,j,t+1)-f(i,j,t)≤T时t₁的情况；其中图(a)为t₁时刻的图像原始的视频输入图像，图(b)为对原始视频图像做本专利提出算法进行非均匀校正后的图像。

图3是视频序列f(i,j,t)满足f(i,j,t+1)-f(i,j,t)＞T的情况；其中图(a)为t₁时刻的图像原始的视频输入图像，图(b)为对原始视频图像做本专利提出算法进行非均匀校正后的图像。

具体实施方式

下面结合附图对本发明作进一步详细描述。

本发明是一种基于各向异性时空域全变分非均匀性校正的方法。其原理为：将图像的看做一个有界变差函数空间，即BV空间，通过建立一个最小化保真项和一个全变分正则项H(u)的过程来实现对图像的复原和去噪，同时考虑到视频处理时帧间的相关性和红外非均匀性校正固定条纹噪声的几何特性，我们对全变分正则项进行改进，竖直方向上的处理项与前面保真项结合在一起，水平和时间方处理项仍用全变分正则项进行平滑，从而在实现去除非均匀性条纹噪声的同时，不会破坏图像的边缘和细节。

一幅含噪图像可以表示为：f(i,j,t)＝u(i,j,t)+n(i,j,t)，其中f(i,j,t)表示实际观察到的图像，u(i,j,t)表示理想的无噪图像，n(i,j,t)表示条纹噪声。我们的目的是根据f(i,j,t)最大限度地还原出u(i,j,t)，从概率的角度来看，即最大化概率P(u/f)，根据贝叶斯原理P(u/f)＝P(f/u)*P(u)，说明P(f/u)最大，同时P(u)最大时，P(u/f)能最大。因此，去噪模型可以写成第一项代表P(f/u)，即保真项，第二项H(u)代表P(u)，即正则项。因为两项之间不是完全独立，a代表权衡系数，一般来说a大的时候图像更平滑，a小的时候图像细节保留的更好。

对于红外图像的非均匀条纹噪声去噪过程，我们取H(u)项为一个全变分模型，这个模型最大的好处就是能对x轴和y轴进行分离的变换，还能沿不同的方向进行扩展，很适合非均匀性噪声(大多数表现为固定的条纹噪声)具有几何稳定性的特定。又由于在处理竖纹噪声时，我们主要要进行水平方向的平滑，垂直方向尽量要宝真，于是得到本专利所提出的去除非均匀性噪声的模型：

其中u为无噪图像，f为含非均匀噪声图像，u_x为空间域水平方向一阶差分，u_t为时间域上的一阶差分。通过分离布拉格曼(Spilt Bregman)法对(1)式求最小值既可以得到所需要的无噪图像u。

结合图1，本发明一种基于各向异性时空域全变分非均匀性校正的方法，包括以下步骤：

步骤1、采集含非均匀性噪声的红外图像序列f(i,j,t)＝u(i,j,t)+n(i,j,t)，其中f(i,j,t)表示实际观察到的图像，u(i,j,t)表示理想的无噪图像，n(i,j,t)表示由红外探测器成像单元的暗电流和响应率不同引起的固定条纹噪声，(i,j)表示图像像素点的位置，t表示视频序列的第t帧图像；同时设每帧图像的总行数为M，总列数为N。

步骤2、采用全变分模型对上述含非均匀性噪声的红外图像序列f(i,j,t)建立图像的全变分模型J(u)，并对图像的全变分模型J(u)进行更新：

其中正则项H(u)取H(u)＝λ₁||u_x||₁+λ₂||u_t||₁ (2)

u_x为图像时空域上对于像素点的x轴的一阶差分，u_t为图像时空域上对于像素点的帧间的一阶差分；λ₁为图像时空域上对于像素点的x轴的一阶差分的平滑系数，λ₂为时空域上对于像素点的帧间的一阶差分的平滑系数。其具体的表达式如下：

对于t₁时刻的图像，对于x轴方向，当j＜M时，取

u_x(i,j,t₁)＝u(i,j+1,t₁)-u(i,j,t₁) (3)对于图像边界处当j＝M，取u_x(i,M,t₁)＝0；u_x(i,M,t₂)＝0；

对于t₁时刻的图像，对于帧间方向，当f(i,j,t₂)-f(i,j,t₁)≤T：

u_t(i,j,t)＝u(i,j,t₂)-u(i,j,t₁)； (4)

得到全变分模型J(u)离散化表示：

对于t₁时刻的图像，对于帧间方向，当f(i,j,t₂)-f(i,j,t₁)>T：

此时J(u)退化为一个仅对图像空间域竖直方向进行校正的单帧模型：

步骤3：通过Split Bregman方法，对式(5)进行迭代求最小化具体的步骤为：

3-1)加入二次惩罚项，通过拉格朗日乘子法，将式(5)中转化为一个无约束最小化问题，得到式(7)；

其中d_x是图像水平方向的二次惩罚项，d_t是时间方向的二次惩罚项，b_x是引入图像水平方向的辅助变量，b_t是引入的时间方向的辅助变量，α和β是方程的Bregman惩罚系数

3-2)将式(7)中u，d_x，d_y这三个未知变量进行分离，得到两个最优化子问题。

对于变量u，分离后最优化子方程为式(8)；

对于变量d_x和d_t，分离后最优化子方程为式(9)；

3-3)式(8)等价于一个基于线性算子的迭代过程，即式(10)：

该方程可以通过快速傅里叶变换(式11)得到u的封闭解：

其中F是快速傅里叶变换，F^-1为傅里叶逆变换。

3-4)采用shrink算子对式(9)的迭代最小化，得到；

其中λ为收缩系数，取0.05，对于收缩算子shrink，对于d_x，当取当取同理，对d_t也采用如上处理方式；

3-5)对于余下的未知数b_x和b_y，则直接采用代入法迭代求解：

3-6)使用Split Bregman进行最优化时，采用改变一个变量并固定其他变量的策略交替迭代，当迭代到第n次时，满足uⁿ⁺¹-uⁿ＜10^-3，此时迭代停止，得到的变量u即为最终的无噪图像u_final。

实施例1

结合图1和图2，当输入视频序列f(i,j,t)恒满足f(i,j,t+1)-f(i,j,t)≤T时，此时视频序列的全变分模型恒为：其去非均匀性噪声步骤为：

步骤1、采集含非均匀性噪声的红外图像序列f(i,j,t)＝u(i,j,t)+n(i,j,t)，其中f(i,j,t)表示实际观察到的图像，u(i,j,t)表示理想的无噪图像，n(i,j,t)表示由红外探测器成像单元的暗电流和响应率不同引起的固定条纹噪声，(i,j)表示图像像素点的位置，t表示视频序列的第t帧图像；同时设每帧图像的总行数为M，总列数为N。

步骤2、对上述含非均匀性噪声的红外图像序列f(i,j,t)建立图像的全变分模型J(u)；

其中正则项H(u)取H(u)＝λ₁||u_x||₁+λ₂||u_t||₁ (2)

u_x为图像时空域上对于像素点的x轴的一阶差分，u_t为图像时空域上对于像素点的帧间的一阶差分；λ₁为图像时空域上对于像素点的x轴的一阶差分的平滑系数，λ₂为时空域上对于像素点的帧间的一阶差分的平滑系数。其具体的表达式如下：

对于t₁时刻的图像，对于x轴方向，当j＜M时；

u_x(i,j,t₁)＝u(i,j+1,t₁)-u(i,j,t₁) (3)

对于图像边界处当j＝M，取u_x(i,M,t₁)＝0；u_x(i,M,t₂)＝0；

对于t₁时刻的图像，对于帧间方向；

u_t(i,j,t)＝u(i,j,t₂)-u(i,j,t₁) (4)

将(2)～(4)式代入式(1)，得到全变分模型J(u)离散化表示：

步骤3：通过Split Bregman方法，对式(5)进行迭代求最小化具体的步骤为：

3-1)加入二次惩罚项，通过拉格朗日乘子法，将式(5)中转化为一个无约束最小化问题，得到式(6)；

其中d_x是图像水平方向的二次惩罚项，d_t是时间方向的二次惩罚项，b_x是引入图像水平方向的辅助变量，b_t是引入的时间方向的辅助变量，α和β是方程的Bregman惩罚系数

3-2)将式(6)中u，d_x，d_y这三个未知变量进行分离，得到两个最优化子问题。

对于变量u，分离后最优化子方程为式(7)；

对于变量d_x和d_t，分离后最优化子方程为式(8)；

3-3)式(8)等价于一个基于线性算子的迭代过程，即式(9)：

该方程可以通过快速傅里叶变换(式10)得到u的封闭解：

其中F是快速傅里叶变换，F^-1为傅里叶逆变换。

3-4)采用shrink算子对式(8)的迭代最小化，得到；

其中λ为收缩系数，取0.05，对于收缩算子shrink，对于d_x，当取当取同理，对d_t也采用如上处理方式；

3-5)对于余下的未知数b_x和b_y，则直接采用代入法迭代求解：

3-6)使用Split Bregman进行最优化时，采用改变一个变量并固定其他变量的策略交替迭代，当迭代到第n次时，满足uⁿ⁺¹-uⁿ＜10^-3，此时迭代停止，得到的变量u即为最终的无噪图像u_final。

实施例2

结合图1和图3，当输入视频序列f(i,j,t)部分满足f(i,j,t+1)-f(i,j,t)≤T，部分满足f(i,j,t+1)-f(i,j,t)＞T的时候，此时视频序列的全变分模型为：

具体数值求解时，当f(i,j,t+1)-f(i,j,t)≤T时，通过Split Bregman同时迭代最小化u_x和u_t，与实施例1完全相同；当f(i,j,t+1)-f(i,j,t)>T时，只需用Split Bregman迭代最小化u_x，最终得到最优解u，具体步骤如下：

步骤1、采集含非均匀性噪声的红外图像序列f(i,j,t)＝u(i,j,t)+n(i,j,t)，其中f(i,j,t)表示实际观察到的图像，u(i,j,t)表示理想的无噪图像，n(i,j,t)表示由红外探测器成像单元的暗电流和响应率不同引起的固定条纹噪声，(i,j)表示图像像素点的位置，t表示视频序列的第t帧图像；

步骤2对上述含非均匀性噪声的红外图像序列f(i,j,t)建立图像的全变分模型J(u)；

其中正则项H(u)取H(u)＝λ₁||u_x||₁+λ₂||u_t||₁ (2)

u_x为图像时空域上对于像素点的x轴的一阶差分，u_t为图像时空域上对于像素点的帧间的一阶差分；λ₁为图像时空域上对于像素点的x轴的一阶差分的平滑系数，λ₂为时空域上对于像素点的帧间的一阶差分的平滑系数。其具体的表达式如下：

对于t₁时刻的图像，对于x轴方向，当j＜M时；

u_x(i,j,t₁)＝u(i,j+1,t₁)-u(i,j,t₁) (3)

其中u_x(i,j,t₁)为t₁时刻图像时空域上对于像素点的x轴的一阶差分；

对于图像边界处当j＝M，取u_x(i,M,t₁)＝0；u_x(i,M,t₂)＝0；

对于t₁时刻的图像，对于帧间方向不更新；

u_t(i,j,t₁)＝0 (4)

其中u_t(i,j,t₁)为t₁时刻图像时空域上对于像素点的帧间的一阶差分；

将(2)～(4)式代入式(1)，得到全变分模型J(u)离散化表示：

步骤3：通过Split Bregman方法，对式(5)进行迭代求最小化具体的步骤为：

3-1)加入二次惩罚项，通过拉格朗日乘子法，将式(5)中转化为一个无约束最小化问题，得到式(6)；

其中d_x是图像水平方向的二次惩罚项，b_x是引入图像水平方向的辅助变量，α是方程的Bregman惩罚系数。

3-2)将式(6)中u，d_x这两个未知变量进行分离，得到两个最优化子问题。

对于变量u，分离后最优化子方程为式(7)；

对于变量d_x，分离后最优化子方程为式(8)；

3-3)式(6)等价于一个基于线性算子的迭代过程，即式(9)：

该方程可以通过快速傅里叶变换(式10)得到u的封闭解：

其中F是快速傅里叶变换，F^-1为傅里叶逆变换。

3-4)采用shrink算子对式(8)迭代最小化，得到；

其中λ为收缩系数，取0.05，对于收缩算子shrink，对于d_x，当取当取

3-5)对于余下的未知数b_x，则直接采用代入法迭代求解：

3-6)使用Split Bregman进行最优化时，采用改变一个变量并固定其他变量的策略交替迭代，当迭代到第n次时，满足uⁿ⁺¹-uⁿ＜10^-3，此时迭代停止，得到的变量u即为最终的无噪图像u_final。