Structural Similarity Fusion of Infrared and Visible Image in the M-SWF Domain
-
摘要: 为了解决常规滤波器组在红外与可见光图像融合领域中存在提取结构信息不充分和融合视觉效果不佳的问题，本文提出了一种基于多尺度滑动窗口滤波器（Multi-scale Sliding Window Filter，M-SWF）图像融合方法。首先，提出一种基于SWF的多尺度图像分解方法实现对源图像的结构细节层和基础层提取；其次，采用L1范数融合规则（L1-Fusion, L1F）整合结构细节层，提取图像的结构信息；然后，利用一种图像能量贡献融合规则（Energy Attribute-Fusion, EAF）整合基础层，突出显著性目标；最后，融合图像通过叠加整合后的多尺度结构细节层和基础层得到。实验首先通过分析能量贡献系数，从主客观方面得到M-SWF域内红外与可见光图像融合较为适宜的能量贡献系数；其次，在该取值下，本文提出的M-SWF融合模型与其他的融合方法相比，不仅提高了对源图像结构信息的提取能力，而且通过整合图像的能量属性，改善了融合效果不佳问题，有效地突出了显著性目标。Abstract: This study introduces a multiscale sliding window filter (M-SWF) image fusion method to address issues with traditional filter banks in infrared and visible image fusion. First, a multiscale image decomposition method based on SWF is proposed to extract the structural detail layers and base layers of the source image. Second, the L1 norm fusion rule (L1-Fusion, L1F) is used to integrate the structural detail layers, which can extract the structure of the image. Then, to highlight the salient objects, energy attribute fusion (EAF), which is a rule for fusing image energy contributions, is used to integrate the base layers, and the fusion results are obtained by stacking the integrated multiscale structure detail layers and base layers. The energy contribution coefficient was analyzed, and a suitable energy contribution coefficient was obtained for the fusion of infrared and visible images in the M-SWF domain from subjective and objective perspectives. Compared with other fusion methods, the M-SWF not only improves the ability to extract the structural information of the source image but also improves the poor fusion effect and effectively highlights salient targets by integrating the energy attributes of the image.
-
Keywords:
- infrared image /
- visible image /
- image fusion /
- side window filter /
- structural similarity
-
0. 引言
地球静止轨道凝视相机可实现固定位置的高帧频成像，具有高空间分辨率、高时间分辨率和高灵敏度探测等优点，在红外目标探测领域发挥着重要作用。由于背景辐射的空间分布不均匀性，相机视线抖动导致落在探测器像元上的辐射能量发生变化，产生视线抖动引起的杂波。抖动引起杂波的强度与背景特征、相机参数、视线稳定度和背景抑制算法等因素有关。由于探测灵敏度和空间分辨率需求的不断提升，抖动引起杂波对相机探测性能的影响逐渐显著。定量化描述和分析视线抖动引起的杂波，有助于系统总噪声的计算、视线稳定度的控制和相机参数的设定与调整。
国内外许多学者对相机视线抖动引起的杂波进行过相关研究。胡方明[1]利用运动传递函数，分析了线性、正弦和随机3种运动形式对红外成像系统性能的影响，建立了成像系统动态性能理论模型，并结合背景杂波对约翰逊准则的修正，对动态性能模型进行了仿真计算。Pohlman[2-3]针对星载凝视相机，度量和分析了系统噪声和杂波，主要包括：探测器和电路的噪声、背景辐射引起的光子噪声、相机视线不稳定性引起的杂波、云的运动和太阳光线变化引起的杂波。Rapier[4-5]提出了杂波等效目标（Clutter-equivalent Target，CET）的快速估计近似方法，在N阶差分算法条件下，将CET表示为背景功率谱密度、光学系统参数、探测距离和相机视线运动特征的函数。Myers[6]针对一阶和二阶差分算法，通过对短波和中波红外谱段下的强反射背景进行模拟，分析了CET与相机空间分辨率、视线运动特征和抖动幅度的关系。Fraedrich[7]评估了帧间差分算法的杂波抑制能力，分析了杂波抑制因子与相机视线特征、探测距离、瞬时视场和帧周期等参数的关系。Casey[8]提出了通过系统建模和权衡分析向系统设计提供反馈的方法，例如：对于采用帧间差分算法检测运动目标的凝视系统，可通过减少抖动而提高系统性能，但需要权衡相应的成本增加。Schroeder[9]和Lee[10]分析了点源探测器在功率谱描述背景中的一维抖动。Liang[11]分析了阶跃辐射背景中的一维抖动，提出了定量评估视线抖动对相机性能影响（信噪比下降）的方法。Hu[12]考虑了二维抖动和相机空间分辨率的影响，但仅使用功率谱模型描述背景，没有考虑抖动引起杂波的空间分布。
本文首先分析了探测器像元上的辐通量梯度与背景辐射强度梯度的关系，推导了焦平面上的辐通量分布与背景辐射强度分布的关联方程，然后综合考虑抖动引起杂波的时间相关影响因素，推导了与背景无关、适用于各个探测器像元的抖动等效角模型，将平均抖动杂波强度表示为抖动等效角和均方根背景辐射强度梯度统计量的乘积，最后通过MATLAB仿真模拟，对抖动引起杂波模型进行了验证。
1. 探测器像元上的辐通量梯度
本文将抖动描述为相机视线与基准位置θx=θy=0的二维时变角度偏差(θx(t), θy(t))。在焦平面坐标(x, y)下，(θx(t), θy(t))对应于焦平面上的图像偏移(Δx(t), Δy(t))，其中Δx(t)=f0θx(t)，Δy(t)=f0θy(t)，f0为光学系统焦距。设E(x, y, t)为t时刻焦平面上点(x, y)处的辐照度，则以(x, y)为中心的单个探测器像元上的辐通量P(x, y, t)为：
$$ P\left(x,y,t\right)={\displaystyle \underset{\begin{array}{l}像元\\ 面积\end{array}}{\iint }E\left(x,y,t\right)\text{d}x\text{d}y} $$ (1) 对于Δx(t0)=Δy(t0)=0的基准时刻t0，则(x′, y′)处的辐照度为E(x′, y′, t0)。当抖动在t时刻将图像平移(Δx, Δy)时，辐照度E(x′, y′, t)为：
$$ E(x′, y′, t)=E(x′+Δx, y′+Δy, t_{0})$$ (2) 因此，焦平面上固定点(x′, y′)处辐照度的时间变化可与t0时刻辐照度的空间变化相关联。由公式(1)和(2)，探测器像元上的辐通量为：
$$ \begin{array}{l}P\left(x,y,t\right)={\displaystyle \underset{\begin{array}{l}像元\\ 面积\end{array}}{\iint }E\left({x}^{\prime }+\Delta x,{y}^{\prime }+\Delta y,{t}_{0}\right)\text{d}{x}^{\prime }\text{d}{y}^{\prime }}\\ \;\;\;\;\;\;\;\;\;=P\left(x+\Delta x,y+\Delta y,{t}_{0}\right)\end{array} $$ (3) 则单个探测器像元上辐通量的时间变化可与t0时刻辐通量的空间变化相关联。对于较小的抖动位移，考虑辐通量的空间变化，上式右边可用泰勒级数线性近似为[11]：
$$ P\left( {x + \Delta x,y + \Delta y,{t_0}} \right) \cong P\left( {x,y,{t_0}} \right) + \frac{{\partial P}}{{\partial x}} \cdot \Delta x + \frac{{\partial P}}{{\partial y}} \cdot \Delta y $$ (4) 因此，探测器像元上辐通量的时间变化为：
$$ P\left( {x,y,t} \right) - P\left( {x,y,{t_0}} \right) = {\left. {\frac{{\partial P}}{{\partial x}}} \right|_{t = {t_0}}} \cdot \Delta x\left( t \right) + {\left. {\frac{{\partial P}}{{\partial y}}} \right|_{t = {t_0}}} \cdot \Delta y\left( t \right) $$ (5) 或用向量表示为：
$$ P\left(x,y,t\right)-P\left(x,y,{t}_{0}\right)=\stackrel{\rightharpoonup }{\nabla }P\left(x,y,{t}_{0}\right)\cdot \stackrel{\rightharpoonup }{\Delta }r\left(t\right) $$ (6) 式中：$ \stackrel{\rightharpoonup }{\Delta }r\left(t\right)=\left(\Delta x\left(t\right),\Delta y\left(t\right)\right) $，$ \stackrel{\rightharpoonup }{\nabla }P $为辐通量P的空间梯度。
辐通量梯度$ \stackrel{\rightharpoonup }{\nabla }P $方向上抖动位移导致的辐通量变化与梯度幅度|$ \stackrel{\rightharpoonup }{\nabla }P $|成正比，而垂直于梯度方向上的抖动位移不会导致辐照度变化。因此，对于二维抖动，只有平行于$ \stackrel{\rightharpoonup }{\nabla }P $的抖动分量是重要的。对于各向同性二维抖动，辐通量变化与|$ \stackrel{\rightharpoonup }{\nabla }P $|成正比。虽然抖动偏移(Δx(t), Δy(t))对于整个焦平面是相同的，但每个探测器像元分别对应具有特定辐通量梯度的不同背景区域，因此，不同探测器像元的响应波动将不同。
2. 背景辐射强度梯度模型
对于地球静止轨道对地凝视红外相机所成图像，设背景辐亮度为L(xe, ye)，其为地面位置坐标(xe, ye)的函数，下标e表示地球坐标系。对一个小矩形区域(Δxe, Δye)上的辐亮度进行积分，得到辐射强度为L(xe, ye)ΔxeΔye，则这个小区域所成图像为以(x′, y′)为中心的焦平面上面积为Δx′Δy′的矩形，其平均辐照度为E(x′, y′)。在没有抖动的基准时刻t0，焦平面坐标与地球坐标的关系为：
$$ x' = \frac{{{f_{\text{o}}}}}{R}{x_{\text{e}}},y' = \frac{{{f_{\text{o}}}}}{R}{y_{\text{e}}} $$ (7) 式中：fo为光学系统焦距；R为相机到地球的距离。
类似地，所成图像区域的尺寸为：
$$ \Delta x' = \frac{{{f_{\text{o}}}}}{R}\Delta {x_{\text{e}}},\Delta y' = \frac{{{f_{\text{o}}}}}{R}\Delta {y_{\text{e}}} $$ (8) 小区域的平均辐照度为收集功率除以成像面积，则：
$$ \begin{gathered} E\left( {x',y',{t_0}} \right) = \frac{{A \cdot \tau }}{{{R^2}}}L\left( {{x_{\text{e}}},{y_{\text{e}}}} \right)\frac{{\Delta {x_{\text{e}}}\Delta {y_{\text{e}}}}}{{\Delta x'\Delta y'}} \hfill \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\; = \frac{{A \cdot \tau }}{{f_{\text{o}}^{\text{2}}}}L\left( {{x_{\text{e}}},{y_{\text{e}}}} \right) \hfill \\ \end{gathered} $$ (9) 式中：A为相机入瞳面积；τ为系统透过率。
由公式(1)、(7)和(9)可知，探测器像元上的辐通量与背景辐亮度的关系为：
$$ \begin{array}{l}P\left(x,y,{t}_{0}\right)={\displaystyle \underset{\begin{array}{l}像元\\ 面积\end{array}}{\iint }E\left({x}^{\prime },{y}^{\prime },t\right)\text{d}{x}^{\prime }\text{d}{y}^{\prime }}\\ \;\;\;\;\;\;\;\;\;\;\;={\displaystyle \underset{\begin{array}{l}像元\\ 面积\end{array}}{\iint }\frac{A\cdot \tau }{{f}_{\text{o}}^{2}}L\left({x}_{\text{e}},{y}_{\text{e}}\right)\text{d}{x}^{\prime }\text{d}{y}^{\prime }}\\ \;\;\;\;\;\;\;\;\;\;\;=\frac{A\cdot \tau }{{R}^{2}}{\displaystyle \underset{\begin{array}{l}瞬时视场\\ 范围\end{array}}{\iint }L\left({x}_{\text{e}},{y}_{\text{e}}\right)\text{d}{x}_{\text{e}}\text{d}{y}_{\text{e}}}\end{array} $$ (10) 又因为梯度幅度$ \left|\stackrel{\rightharpoonup }{\nabla }P\right| $的定义为：
$$ \left|\stackrel{\rightharpoonup }{\nabla }P\right|={\left({\left(\frac{\partial P}{\partial x}\right)}^{2}+{\left(\frac{\partial P}{\partial y}\right)}^{2}\right)}^{1/2} $$ (11) 所以，需要对上式求导。而地面分辨率对x和y的约束为：
$$ \frac{R}{{{f_{\text{o}}}}}x - \frac{{\Delta {x_{\text{c}}}}}{2} < {x_{\text{e}}} < \frac{R}{{{f_{\text{o}}}}}x + \frac{{\Delta {x_{\text{c}}}}}{2} $$ (12) $$ \frac{R}{{{f_{\text{o}}}}}y - \frac{{\Delta {y_{\text{c}}}}}{2} < {y_{\text{e}}} < \frac{R}{{{f_{\text{o}}}}}y + \frac{{\Delta {y_{\text{c}}}}}{2} $$ (13) 式中：ΔxcΔyc为地面坐标系中地面分辨率大小。将上式对x和y分别求导，从积分极限中提取出R/F因子，则：
$$ \frac{{\partial P}}{{\partial x}} = \frac{{A \cdot \tau }}{{R \cdot {f_{\text{o}}}}}\int\limits_{{y_{\text{c}}} - \frac{{\Delta {y_{\text{c}}}}}{2}}^{{y_{\text{c}}} + \frac{{\Delta {y_{\text{c}}}}}{2}} {\left[ {N\left( {{x_{\text{c}}} + \frac{{\Delta {x_{\text{c}}}}}{2},{y_{\text{e}}}} \right) - N\left( {{x_{\text{c}}} - \frac{{\Delta {x_{\text{c}}}}}{2},{y_{\text{e}}}} \right)} \right]{\text{d}}{y_{\text{e}}}} $$ (14) $$ \frac{{\partial P}}{{\partial y}} = \frac{{A \cdot \tau }}{{R \cdot {f_{\text{o}}}}}\int\limits_{{x_{\text{c}}} - \frac{{\Delta {x_{\text{c}}}}}{2}}^{{x_{\text{c}}} + \frac{{\Delta {x_{\text{c}}}}}{2}} {\left[ {N\left( {{x_{\text{e}}},{y_{\text{c}}} + \frac{{\Delta {y_{\text{c}}}}}{2}} \right) - N\left( {{x_{\text{e}}},{y_{\text{c}}} - \frac{{\Delta {y_{\text{c}}}}}{2}} \right)} \right]{\text{d}}{x_{\text{e}}}} $$ (15) 式中：(xc, yc)为瞬时视场范围的中心。
上面两式右边的积分值完全取决于背景和地面分辨率，可从二维图像数据中进行评估。将背景表示为像素大小为ΔxΔy的二维阵列，Lm, n为第m行和第n列像素的辐亮度。设瞬时视场范围对应k×k像素的方形，(i, j)为k×k像素方形角落处的图像像素，Li, j为其辐亮度，则x方向上的梯度分量近似为：
$$ \frac{{\partial P}}{{\partial x}} = \frac{{A \cdot \tau }}{{R \cdot {f_{\text{o}}}}}\sum\limits_{p = i}^{i + k - 1} {\left( {{L_{p,j + k}} - {L_{p,j}}} \right)} \Delta y $$ (16) 式中：∂P/∂x的单位为W/(sr⋅m)，可通过下式计算每个像素上的辐射强度bm, n得到：
$$ b_{m,n}=L_{m,n}ΔxΔy $$ (17) 将∂P/∂x重新表示为：
$$ \frac{{\partial P}}{{\partial x}} = \frac{{A \cdot \tau }}{R}\frac{1}{{\Delta x}}\sum\limits_{p = i}^{i + k - 1} {\left( {{b_{p,j + k}} - {b_{p,j}}} \right)} $$ (18) 类似地，y方向上的梯度分量为：
$$ \frac{{\partial P}}{{\partial y}} = \frac{{A \cdot \tau }}{R}\frac{1}{{\Delta y}}\sum\limits_{p = j}^{j + k - 1} {\left( {{b_{i + k,p}} - {b_{i,p}}} \right)} $$ (19) 公式(18)和(19)结合梯度的定义式，给出了求$ \left|\stackrel{\rightharpoonup }{\nabla }P\right| $的所有信息。
$ \left|\stackrel{\rightharpoonup }{\nabla }P\right| $的计算结果可由背景辐射强度梯度表示。定义背景辐射强度梯度统计量$ \left|\stackrel{\rightharpoonup }{\nabla }I\right| $，其为瞬时视场范围内背景辐射强度随单位位移的变化量。当k=3时，瞬时视场范围的位移如图 1所示。因此，$ \left|\stackrel{\rightharpoonup }{\nabla }I\right| $的计算式为：
$$ \begin{array}{l}{\left|\stackrel{\rightharpoonup }{\nabla }I\right|}^{2}={\left[\frac{1}{\Delta x}{\displaystyle \sum _{p=i}^{i+k-1}\left({b}_{p,j+k}-{b}_{p,j}\right)}\right]}^{2}+\\ \;\;\;\;\;\;\;\;{\left[\frac{1}{\Delta y}{\displaystyle \sum _{p=j}^{j+k-1}\left({b}_{i+k,p}-{b}_{i,p}\right)}\right]}^{2}\end{array} $$ (20) ![]()
图 1 瞬时视场范围（3×3）移动一个像素前后示意图Figure 1. Schematic diagram of the IFOV (3×3) before and after moving by one pixel因此，可建立辐通量梯度幅度$ \left|\stackrel{\rightharpoonup }{\nabla }P\right| $与背景辐射强度梯度统计量$ \left|\stackrel{\rightharpoonup }{\nabla }I\right| $之间的关系。由于抖动是在角坐标系下描述的，因此角坐标下的辐通量梯度幅度为$ \left|\stackrel{\rightharpoonup }{\nabla }P\right|\cdot {f}_{\text{o}} $。由公式(11)、(18)、(19)和(20)可知，$ \left|\stackrel{\rightharpoonup }{\nabla }P\right| $可由$ \left|\stackrel{\rightharpoonup }{\nabla }I\right| $表示为：
$$ \left|\stackrel{\rightharpoonup }{\nabla }P\right|\cdot {f}_{o}=\frac{A\cdot \tau }{R}\cdot \left|\stackrel{\rightharpoonup }{\nabla }I\right| $$ (21) 由于相机所成图像特征与背景辐射强度空间梯度相关，本节推导了将相机焦平面上的辐通量分布与背景辐射强度分布关联起来的方程，分析了探测器像元上的辐通量梯度$ \stackrel{\rightharpoonup }{\nabla }P $与瞬时视场范围内背景辐射强度梯度$ \left|\stackrel{\rightharpoonup }{\nabla }I\right| $的关系。
3. 抖动等效角模型
对于以一定背景辐射强度梯度为中心的单个探测器像元，抖动杂波与抖动频谱、探测器积分时间、帧周期和背景抑制算法有关。这些因素可以合并为一个对所有探测器像元都相同的性能参数σJ，即系统的抖动等效角。
对于单个探测器像元，在积分时间Tint内积分得到电子数n(t) [2]：
$$ n\left( t \right) = \int\limits_0^{{T_{\rm{int} }}} {\frac{{\lambda \eta }}{{hc}}P\left( {x,y,t'} \right){\text{d}}t'} $$ (22) 式中：λ为波长；η为量子效率；h为普朗克常量；c为光速。分析可知，电子数n(t)波动的统计特性取决于抖动的统计特性。设θ(t)为梯度方向上的时变抖动角度偏差；Ψ(τ)为其相关函数，称为抖动相关函数，定义为：
$$ \varPsi (τ)=＜θ(t)θ(t－τ)＞$$ (23) 在各向同性二维抖动的情况下，Ψ(τ)也可以描述任何其他方向的抖动特性。
由所建立的辐通量梯度模型可知：
$$ P(x,y,t)=P(x,y,t_{0})+|\nabla P| \cdot f_{\text{o}} \cdot θ(t) $$ (24) 因此，辐通量变化的相关函数为：
$$ < P\left( t \right)P\left( {t - \tau } \right) > = {\left| {\nabla P} \right|^2} \cdot f_{\text{o}}^2 \cdot \psi \left( \tau \right) $$ (25) 对于帧间差分背景抑制算法，抖动杂波的强度可由帧间差分法的均方输出（以电子数形式）来描述。若相机帧周期为s，则t时刻帧间差分法的均方输出为：
$$ \begin{align}＜n(s,t)^{2}＞= &＜[n(t)－n(t－s)]^{2}＞= \\&2[＜n^{2}(t)＞－＜n(t)n(t－s)＞]\end{align}$$ (26) 式中：＜n＞表示均值；＜n2＞表示均方值。
上式的最后一步是假设抖动为平稳随机过程，即n(t)的变化是平稳随机过程，因此，＜n(t－s)2＞=＜n2(t)＞。
采样相关函数R(τ)为：
$$ R(τ)=＜n(t)n(t－τ)＞ $$ (27) 则帧间差分法的均方输出为：
$$ ＜n(s,t)^{2}＞=2R(0)－2R(s) $$ (28) 经变量变换t′=t－Tint+v，式(22)可表示为：
$$ n\left( t \right) = \int\limits_0^{{T_{\rm{int} }}} {\frac{{\lambda \eta }}{{hc}}P\left( {t - {T_{\rm{int} }} + v} \right){\text{d}}v} $$ (29) 因此，R(τ)可由二重积分表示为：
$$ R\left( \tau \right) = {\left( {\frac{{\lambda \eta }}{{hc}}} \right)^2}\int\limits_0^{{T_{\rm{int} }}} {\int\limits_0^{{T_{\rm{int} }}} { < P\left( {t - {T_{\rm{int} }} + v} \right)P\left( {t - \tau - {T_{\rm{int} }} + w} \right) > {\text{d}}v} {\text{d}}w} $$ (30) 根据上述方程，上式可以用抖动相关函数表示为：
$$ \begin{gathered} R\left( \tau \right) = {\left( {\frac{{\lambda \eta }}{{hc}}\left| {\nabla P} \right|{f_{\text{o}}}} \right)^2}\int\limits_0^{{T_{\rm{int} }}} {\int\limits_0^{{T_{\rm{int} }}} {\psi \left( {v - w + \tau } \right){\text{d}}v{\text{d}}w} } \hfill \\ = {\left( {\frac{{\lambda \eta }}{{hc}}\left| {\nabla P} \right|{f_{\text{o}}}} \right)^2}\int\limits_0^{{T_{\rm{int} }}} {\left( {{T_{\rm{int} }} - u} \right)\left[ {\psi \left( {\tau + u} \right) + \psi \left( {\tau - u} \right)} \right]{\text{d}}u} \hfill \\ \end{gathered} $$ (31) 则采样相关函数R(0)和R(s)可用抖动相关函数表示为[11]：
$$ R\left( 0 \right) = {\left( {\frac{{\lambda \eta }}{{hc}}\left| {\nabla P} \right|{f_{\text{o}}}} \right)^2}\int\limits_0^{{T_{\rm{int} }}} {\left( {{T_{\rm{int} }} - u} \right)\left[ {\psi \left( u \right) + \psi \left( u \right)} \right]{\text{d}}u} $$ (32) $$ R\left( s \right) = {\left( {\frac{{\lambda \eta }}{{hc}}\left| {\nabla P} \right|{f_{\text{o}}}} \right)^2}\int\limits_0^{{T_{\rm{int} }}} {\left( {{T_{\rm{int} }} - u} \right)\left[ {\psi \left( {s + u} \right) + \psi \left( {s - u} \right)} \right]{\text{d}}u} $$ (33) 抖动统计量通常用功率谱密度ω(f)表示，而功率谱密度与抖动相关函数之间存在以下关系式[15]：
$$ \omega \left( f \right) = 4\int\limits_0^\infty {\psi \left( \tau \right)\cos 2{\text{π }}f\tau {\text{d}}\tau } $$ (34) $$ \psi \left( \tau \right) = \int\limits_0^\infty {\omega \left( f \right)\cos 2{\text{π }}f\tau {\text{d}}f} $$ (35) 则将采样相关函数用功率谱密度表示为：
$$ R\left( 0 \right) = {\left( {\frac{{\lambda \eta }}{{hc}}\left| {\nabla P} \right|{f_{\text{o}}}} \right)^2}\int\limits_0^\infty {2\omega \left( f \right)\frac{1}{{{{\left( {2{\text{π }}f} \right)}^2}}} \cdot \left[ {1 - \cos 2{\text{π }}f{T_{\rm{int} }}} \right]{\text{d}}f} $$ (36) $$ \begin{array}{l} R\left( s \right) = {\left( {\frac{{\lambda \eta }}{{hc}}\left| {\nabla P} \right|{f_{\text{o}}}} \right)^2} \hfill \\ \int\limits_0^\infty {2\omega \left( f \right)\cos \left( {2{\text{π }}fs} \right)\frac{1}{{{{\left( {2{\text{π }}f} \right)}^2}}} \cdot \left[ {1 - \cos 2{\text{π }}f{T_{\rm{int} }}} \right]{\text{d}}f} \hfill \\ \end{array} $$ (37) 最后，将公式(33) (34)代入公式(28)，得到帧间差分法的均方输出的表达式为：
$$ < n{\left( {s,t} \right)^2} > = {\left( {\frac{{\lambda \eta }}{{hc}}\left| {\nabla P} \right|{f_{\text{o}}}} \right)^2}T_{\rm{int} }^2\int\limits_0^\infty {\omega \left( f \right)S\left( {f,s,{T_{\rm{int} }}} \right){\text{d}}f} $$ (38) 其中：
$$ S\left( {f,s,{T_{\rm{int} }}} \right) = \frac{1}{{{{\left( {{\text{π }}f{T_{\rm{int} }}} \right)}^2}}}\left[ {1 - \cos \left( {2{\text{π }}f{T_{\rm{int} }}} \right)} \right]\left[ {1 - \cos \left( {2{\text{π }}fs} \right)} \right] 。 $$ 因此，抖动等效角σJ可由抖动功率谱密度ω(f)和帧间差分算法传递函数S表示为：
$$ \sigma _{\text{J}}^2 = \int\limits_0^\infty {\omega \left( f \right)S\left( {f,s,{T_{\rm{int} }}} \right){\text{d}}f} $$ (39) 由于ω(f)的单位为μrad2/Hz，S无量纲，因此σJ的单位为μrad。
由以上推导过程可知，抖动等效角实际上描述了在探测器积分时间和帧间差分背景抑制算法的通频带内抖动功率大小，表征了抖动频谱、探测器积分时间、帧周期和帧间差分背景抑制算法对抖动杂波的影响，对焦平面上每个探测器像元都相同。
4. 抖动引起杂波模型
由公式(21)(38)和(39)，对于单个探测器像元，帧间差分法的均方根输出为：
$$ \begin{array}{l}\sqrt{ < n{\left(s,t\right)}^{2}{ > }_{i,j}}=\frac{\lambda \cdot \eta \cdot A\cdot \tau \cdot {\sigma }_{\text{J}}\cdot \left|\stackrel{\rightharpoonup }{\nabla }I\right|\cdot {T}_{\mathrm{int}}}{h\cdot c\cdot R}\\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;=\frac{\text{π}\lambda \cdot \eta \cdot {D}^{2}\cdot \tau \cdot {\sigma }_{\text{J}}\cdot \left|\stackrel{\rightharpoonup }{\nabla }I\right|\cdot {T}_{\mathrm{int}}}{4h\cdot c\cdot R}\end{array} $$ (40) 由于上式中除$ \left|\stackrel{\rightharpoonup }{\nabla }I\right| $之外的所有因子对于每个探测器像元都相同，因此，通过求$ \left|\stackrel{\rightharpoonup }{\nabla }I\right| $的均方根，可得到对整个焦平面计算的帧间差分法均方根输出为：
$$ \sqrt{ < {n}^{2} > }=\frac{\text{π}\lambda \cdot \eta \cdot {D}^{2}\cdot \tau \cdot {\sigma }_{\text{J}}\cdot \sqrt{ < \stackrel{\rightharpoonup }{\nabla }{I}^{2} > }\cdot {T}_{\mathrm{int}}}{4h\cdot c\cdot R} $$ (41) 根据辐亮度与电子数的关系，则整个焦平面的平均抖动引起杂波强度为：
$$ \begin{array}{l}\text{CET}=\frac{4h\cdot c\cdot {f}_{\text{o}}^{2}\cdot \sqrt{ < {n}^{2} > }}{\text{π}\lambda \cdot \eta \cdot \tau \cdot {D}^{2}\cdot {x}_{\text{d}}^{2}\cdot {T}_{\mathrm{int}}}\cdot \Delta {x}_{\text{c}}\cdot \Delta {y}_{\text{c}}\\ \;\;\;\;\;\;={\sigma }_{\text{J}}\cdot \sqrt{ < \stackrel{\rightharpoonup }{\nabla }{I}^{2} > }\cdot R\end{array} $$ (42) 由公式(42)式可知，平均抖动杂波强度与抖动等效角σJ和均方根背景辐射强度梯度统计量$ \sqrt{ < \stackrel{\rightharpoonup }{\nabla }{I}^{2} > } $成正比。式(42)显示了如何将抖动杂波强度模型分解为几个独立的分量。其中，抖动等效角包含了相机的时间相关特性，而均方根背景辐射强度梯度统计量表征了背景的空间结构。在相同观测条件和系统条件下，当采用同一相机观测不同背景时，式(42)中唯一变化项为均方根背景辐射强度梯度统计量；因此，其为描述抖动杂波强度的关键统计量。
5. 抖动引起杂波模型验证
针对地球静止轨道凝视红外相机，为了验证相机视线抖动引起的杂波模型，考虑在轨影响因素的多样性和测量仪器设备的有限性，首先开展了相机视线抖动特性地面测量实验，基于角位移传感器数据和实际相机参数，分析了相机视线抖动频谱，计算了抖动等效角。然后根据所建立的抖动引起杂波模型，计算了不同背景下的抖动杂波强度。最后以在轨实测图像为背景，通过MATLAB进行抖动条件下的图像序列仿真模拟，将抖动杂波强度的仿真模拟与模型计算结果进行了对比和分析。本文所研究相机的具体参数如表 1所示。
表 1 相机参数Table 1. The sensor parametersItem Value Average wavelength/μm 2.85 Quantum efficiency 0.6 Planck constant/J⋅s 6.626×10-34 Speed of light/m⋅s－1 3×108 System transmittance 0.6 Ground sample distance/m 600 Detector size/μm 15×15 Optical aperture/m 0.35 Distance from sensor to Earth/km 35786 Integration time/s 0.4 Frame period/s 1 5.1 视线抖动频谱分析和抖动等效角计算
相对于低轨卫星相机，地球静止轨道（高轨）凝视红外相机成像积分时间长，更容易受到抖动的影响。根据卫星设计方案，星上活动部件主要包括：数传天线、太阳电池阵驱动机构（Solar Array Drive Assembly，SADA）、动量轮和脉冲管制冷机。数传天线在相机成像时不运动，SADA的基频较低，能量主要集中在0.2~20 Hz的低频段内[16]。通过对SADA的扰动测试可知，相比于制冷机和动量轮，SADA的扰动影响较小[17]。因此，本文将制冷机和动量轮作为地球静止轨道凝视红外相机视线抖动的主要扰动源[18]，基于相机主镜位置的角位移传感器数据，分析了制冷机和动量轮引起的相机视线抖动频谱。
在轨工作期间，制冷机对红外焦平面进行低温制冷，从而保证探测器的成像任务。当制冷机工作时，高压气体压力的波动和活塞等部件动量的不平衡，均会导致相机视线抖动[19]。在与相机视线垂直的x和y方向上，制冷机引起的相机视线抖动频谱如图 2和3所示。
![]()
图 2 制冷机引起相机视线在x方向上的抖动频谱Figure 2. Dither spectra of the camera line-of-sight in the x direction caused by cryocooler![]()
图 3 制冷机引起相机视线在y方向上的抖动频谱Figure 3. Dither spectra of the camera line-of-sight in the y direction caused by cryocooler由图 2和图 3可知，制冷机引起的相机视线抖动表现为典型的谐波特性，频率为制冷机工作频率和工作频率的倍频。
动量轮包含高速转动的惯性组件，由于加工误差等原因，惯性组件的质心与实际转动轴不重合，会导致转动不平衡，从而引起相机视线抖动[20]。本文实验中的卫星配备了0.5 Nm和0.1 Nm两种规格的动量轮，工况为多台动量轮同时开机。在x和y方向上，动量轮引起的相机视线抖动频谱如图 4和图 5所示。
![]()
图 4 动量轮引起相机视线在x方向上的抖动频谱Figure 4. Dither spectra of the camera line-of-sight in the x direction caused by momentum wheels![]()
图 5 动量轮引起相机视线在y方向上的抖动频谱Figure 5. Dither spectra of the camera line-of-sight in the y direction caused by momentum wheels由图 4和5可知，动量轮引起的相机视线抖动表现为典型的谐波和宽频噪声特性。
根据相机视线抖动频谱，结合实际相机参数，可计算不同扰动源引起的相机视线抖动等效角，如表 2所示。
表 2 不同扰动源引起的相机视线抖动等效角Table 2. The dither-equivalent angle caused by different disturbance sourcesx Direction y Direction Total angle Jitter-equivalent angle caused by cryocooler/μrad 0.00027 0.0016 0.00163 Jitter-equivalent angle caused by momentum wheel/μrad 0.0747 0.0829 0.112 由表 2可知，制冷机引起的相机视线抖动等效角比动量轮引起的相机视线抖动等效角小近两个数量级。根据所建立的抖动等效角和抖动引起杂波模型，相比于制冷机，动量轮的扰动为影响抖动杂波强度的主要因素。
5.2 仿真模拟和模型计算结果对比
根据上文建立的抖动引起杂波模型，即公式(42)，通过抖动等效角σJ、均方根背景辐射强度梯度统计量$ \sqrt{ < \stackrel{\rightharpoonup }{\nabla }{I}^{2} > } $和相机到地球的距离R，可计算不同背景下的平均抖动杂波强度，作为抖动杂波强度的模型计算结果。
在相同观测条件和系统条件下，根据地球静止轨道凝视红外相机在轨所成图像，选取典型的海面、陆地和云层背景区域，根据5.1节中的视线抖动频谱，生成帧间偏移量，在MATLAB中模拟生成抖动条件下的图像序列；对图像序列中相邻两帧图像做差分运算，统计每个像素灰度值的均方根和整个背景区域的均方根值，并转换为等效辐射强度，以此作为抖动杂波强度的仿真模拟结果。
图 6显示了典型的海面、陆地和云层背景图像，图 7显示了不同背景的辐射强度梯度分布（单位：W/sr/m），图 8显示了不同背景下仿真模拟的抖动杂波强度分布（单位：W/sr）。为了清晰地表征物理量的数值对比和空间分布特征，图 7和8采用了二维彩色云图的形式。
![]()
图 7 不同背景的辐射强度梯度分布Figure 7. Radiation intensity gradient distribution of different backgrounds![]()
图 8 不同背景下仿真模拟的抖动杂波强度分布Figure 8. Simulated dither-caused clutter intensity distribution of different backgrounds表 3显示了不同背景下平均抖动杂波强度的模型计算与仿真模拟结果。表 3中均方根背景辐射强度梯度统计量来自公式(20)和图 7；模型计算的平均抖动杂波强度来自公式(42)，即通过抖动等效角σJ、均方根背景辐射强度梯度统计量$ \sqrt{ < \stackrel{\rightharpoonup }{\nabla }{I}^{2} > } $和相机到地球的距离R计算得到；仿真模拟的平均抖动杂波强度来自图 8。
表 3 不同背景下的平均抖动杂波强度Table 3. Average clutter intensity of different backgroundsBackground RMS background radiation
intensity gradient statistics/(W/sr/m)CET calculated from
the model/(W/sr)CET simulated/
(W/sr)Relative deviation Sea 0.0624 0.25 0.23 8.7% Land 0.2671 1.07 1.01 5.9% Clouds 0.6382 2.57 2.86 10.1% 由表 3可知，不同背景下平均抖动杂波强度的模型计算与仿真模拟结果相近，相对偏差小于15%。在相同系统条件和不同背景下，计算结果符合背景辐射强度梯度越大则抖动杂波强度越大的特性，间接证明了所建立的抖动引起杂波模型的合理性。
所建立的抖动引起杂波模型给出了较合理的抖动杂波强度，与仿真模拟结果具有一致性，可简化描述背景特征、相机参数、相机视线抖动特性与抖动引起杂波的关系，而无需依赖抖动条件下的图像序列仿真模拟，便于系统设计中抖动杂波强度的估计、相机参数的设定和视线稳定度的控制。
6. 结论
对于地球静止轨道凝视红外相机，相机视线抖动是时域杂波的重要来源。针对使用功率谱描述背景的杂波模型会导致边缘信息缺失问题，本文探讨了探测器像元上的辐通量梯度与背景辐射强度梯度的关系，建立了相机视线抖动等效角和抖动引起杂波模型。该模型表明：探测器像元的抖动杂波强度与它对应的背景辐射强度梯度成正比，整个焦平面的平均抖动杂波强度与均方根背景辐射强度梯度统计量成正比；因此，抖动杂波强度的分布与背景辐射强度梯度分布相关。本文所建立的抖动引起杂波模型可以计算焦平面上的抖动杂波强度，具有较好的通用性和高效性，可为优化凝视红外相机设计、改进目标检测算法和提高系统探测性能提供依据。
-
![]()
图 5 不同方法在“Camp_1811”图像下的融合结果
Figure 5. The fusion results by different methods under "Camp_1811"
![]()
图 6 不同方法在“Trees_4917”图像下的融合结果
Figure 6. The fusion results by different methods under "Trees_4917"
![]()
图 7 不同方法在“Kaptein_1123”图像下的融合结果
Figure 7. The fusion results by different methods under "Kaptein_1123"
表 1 不同P值下的融合定量数据
Table 1 Fusion of quantitative data under different P
Qabf MIP SCD SSIM MSSIM VIFF RT/s Fig.4(c) 0.5007 0.9249 1.5446 0.7640 0.9071 0.3044 0.510 Fig.4(d) 0.5303 0.9329 1.7225 0.7637 0.9287 0.3311 0.513 Fig.4(e) 0.5293 0.9364 1.7823 0.7604 0.9287 0.3323 0.508 Fig.4(f) 0.5233 0.9380 1.7936 0.7565 0.9350 0.3407 0.504 Fig.4(g) 0.5163 0.9382 1.7969 0.7527 0.9181 0.3512 0.518 Fig.4(h) 0.5096 0.9369 1.7989 0.7490 0.9135 0.3626 0.506 
下载: 导出CSV
表 2 客观评价数据
Table 2 Objective evaluation data
Qabf MIP SCD SSIM MSSIM VIFF RT/s Camp_1811 GFF 0.4241 0.8979 1.1453 0.6664 0.7186 0.2724 0.0990 CBF 0.3920 0.8701 1.3044 0.6238 0.7438 0.2839 3.4876 FFIF 0.4562 0.9104 1.2300 0.7257 0.8442 0.2933 0.3520 Resnet 0.3392 0.8823 1.4860 0.7777 0.8707 0.3203 12.3780 MLEPF 0.3640 0.8463 1.2691 0.6230 0.6878 0.2791 8.9630 M-SWF 0.4802 0.8857 1.7003 0.7476 0.9249 0.4410 0.2520 Trees_4917 GFF 0.2526 0.8807 0.9299 0.7559 0.8560 0.2651 0.1160 CBF 0.3296 0.8086 1.2718 0.5963 0.6850 0.5599 3.1175 FFIF 0.2271 0.8928 0.6808 0.7627 0.7782 0.0877 0.4260 Resnet 0.3847 0.8492 1.2224 0.8462 0.9249 0.2975 11.8230 MLEPF 0.4963 0.7815 0.4175 0.5972 0.6365 0.3880 9.1730 M-SWF 0.5431 0.8626 1.6195 0.8095 0.9597 0.4528 0.2370 Kaptein_1123 GFF 0.1079 0.8486 0.4842 0.4425 0.2967 0.3722 0.2830 CBF 0.4315 0.8645 1.4597 0.5370 0.6509 0.3280 8.9544 FFIF 0.3053 0.8954 1.2110 0.6758 0.8472 0.3671 0.4770 Resnet 0.3303 0.8911 1.6111 0.7467 0.8582 0.3024 40.4150 MLEPF 0.4110 0.8171 1.3662 0.5523 0.5958 0.1951 20.7530 M-SWF 0.5778 0.8955 1.8601 0.7183 0.9423 0.4522 0.4840
下载: 导出CSV
表 3 TNO数据集上的平均客观评价数据
Table 3 Average objective evaluation data on TNO
Qabf MIP SCD SSIM MSSIM VIFF RT/s GFF 0.2899 0.9020 1.1450 0.6736 0.7554 0.2488 0.2118 CBF 0.6054 0.8897 1.1488 0.6679 0.7864 0.2895 9.7087 FFIF 0.2881 0.9016 1.1914 0.7098 0.8050 0.2271 0.4828 Resnet 0.4093 0.8683 1.2537 0.6625 0.7516 0.3145 27.2588 MLEPF 0.4325 0.8514 1.3164 0.6006 0.6582 0.2634 20.2127 M-SWF 0.5691 0.9084 1.8071 0.7476 0.9317 0.3827 0.3944
下载: 导出CSV
-
[1] MA Jiayi, MA Yong, LI Chang. Infrared and visible image fusion methods and applications: a survey [J]. Information Fusion, 2019, 45: 153-178. DOI: 10.1016/j.inffus.2018.02.004。
[2] ZHAO Z, XU S, ZHANG C, et al. Bayesian fusion for infrared and visible images[J]. Signal Processing, 2020, 177: 165-168. DOI: 10.1016/ j.sigpro.2020.107734
[3] HUAN Kewei, LI Xiangyang, CAO Yutong, et al. Infrared and visible image fusion of convolutional neural network and NSST[J]. Infrared and Laser Engineering, 2022, 51(3): 20210139. DOI: 10.3788/IRLA20210139.
[4] CHENG Boyang, LI Ting, WANG Yulin. Fusion of infrared and visible light images based on visual saliency weighting and maximum gradient singular value[J]. Chinese Optics, 2022, 15(4): 675-688. DOI: 10.37188/CO.2022-0124
[5] 李威, 李忠民. 一种基于EASSF的红外与可见光图像视觉保真度融合[J]. 红外技术, 2022, 44(7): 686-692. http://hwjs.nvir.cn/article/id/2c60451f-34b1-47a3-bab6-629774a73a0c LI Wei, LI Zhongmin. Visual fidelity fusion of infrared and visible image using edge-aware smoothing-sharpening filter[J]. Infrared Technology, 2022, 44(7): 686-692. http://hwjs.nvir.cn/article/id/2c60451f-34b1-47a3-bab6-629774a73a0c
[6] 李永萍, 杨艳春, 党建武, 等. 基于变换域VGGNet19的红外与可见光图像融合[J]. 红外技术, 2022, 44(12): 1293-1300. http://hwjs.nvir.cn/article/id/48843dec-a48b-468f-b743-4a00c345f406 LI Yongping, YANG Yanchun, DANG Jianwu, et al. Infrared and visible image fusion based on transform domain VGGNet19[J]. Infrared Technology, 2022, 44(12): 1293-1300. http://hwjs.nvir.cn/article/id/48843dec-a48b-468f-b743-4a00c345f406
[7] 雷大江, 杜加浩, 张莉萍, 等. 联合多流融合和多尺度学习的卷积神经网络遥感图像融合方法[J]. 电子与信息学报, 2022, 44(1): 237-244. Doi: 10.11999/JEIT200792. LEI Dajiang, DU Jiahao, ZHANG Liping, et al. Multi-stream architecture and multi-scale convolutional neural network for remote sensing image fusion [J]. Journal of Electronics & Information Technology, 2022, 44(1): 237-244. Doi: 10.11999/JEIT200792.
[8] 马梁, 苟于涛, 雷涛, 等. 基于多尺度特征融合的遥感图像小目标检测[J]. 光电工程, 2022, 49(4): 49-65. https://www.cnki.com.cn/Article/CJFDTOTAL-GDGC202204005.htm MA Liang, GOU Yutao, LEI Tao, et al. Small object detection based on multi-scale feature fusion using remote sensing images[J]. Opto-Electron Eng, 2022, 49(4): 49-65. https://www.cnki.com.cn/Article/CJFDTOTAL-GDGC202204005.htm
[9] 钱金卓, 马骏, 李峰, 等. 面向CMOS遥感相机的多曝光图像融合方法[J]. 遥感信息, 2022, 37(4): 51-57. https://www.cnki.com.cn/Article/CJFDTOTAL-YGXX202204008.htm QIAN Jinzhuo, MA Jun, LI Feng, et al. Multi-exposure image fusion method for CMOS remote sensing camera [J]. Remote Sensing Information, 2022, 37(4): 51-57. https://www.cnki.com.cn/Article/CJFDTOTAL-YGXX202204008.htm
[10] LI Shutao, KANG Xudong, HU Jianwen. Image fusion with guided filtering[J]. IEEE Transactions on Image Processing, 2013, 22(7): 2864-2875. Doi: 10.1109/TIP.2013.2244222.
[11] Shreyamsha Kumar B K. Image fusion based on pixel significance using cross bilateral filter[J]. Signal Image Video Process, 2015, 9(5): 1193-1204. Doi: 10.1007/s11760-013-0556-9.
[12] ZHAN K, XIE Yuange, MIN Yufang. Fast filtering image fusion[J]. J. Electron. Imaging, 2017, 26(6): 063004. Doi: 10.1117/1.JEI.26.6.063004.
[13] LI Hui, WU Xiaojun, Tariq S Durrani. Infrared and visible image fusion with ResNet and zero-phase component analysis[J]. Infrared Physics & Technology, 2019, 102: 1030390. Doi: 10.1016/j.infrared.2019.-103039.
[14] TAN Wei, Thitn W, XIANG P, et al. Multi-modal brain image fusion based on multi-level edge-preserving filtering[J]. Biomedical Signal Processing and Control, 2021, 64: 102280. Doi: 10.1016/j. bspc.2020.102280.
[15] YIN Hui, GONG Yuanhao, QIU Guoping. Side window guided filtering[J]. Signal Process, 2019, 165: 315-330. Doi: 10.1016/j.sigpro. 2019.07.026.
[16] LI Xiaosong, ZHOU Fuqiang, TAN Haishu, et al. Multimodal medical image fusion based on joint bilateral filter and local gradient energy[J]. Information Sciences, 2021, 569: 302-325.
[17] LI Hui, WU Xiaojun. DenseFuse: a fusion approach to infrared and visible images [J]. IEEE Trans. Image Process, 2019, 28(5): 2614-2623. Doi: 10.1109/TIP.2018.2887342.
[18] Toet A. TNO Image Fusion Dataset[EB/OL]. [2022-10-30]. http://figshare.com/articles/-TNO_Image_Fusion_Dataset/1008029.
[19] Xydeas C S, Petrovic V. Objective image fusion performance measure [J]. Electron. Lett, 2000, 36(4): 308-309. Doi: 10.1109/ICCV.2005.175.
[20] Aslantas V, Bendes E. A new image quality metric for image fusion: the sum of the correlations of differences[J]. AEU- Int. J. Electron. Commune, 2015, 69(12): 1890-1896. Doi: 10.1016/j.aeue.2015.09.004.
计量
- 文章访问数:
- HTML全文浏览量:
- PDF下载量:
