一种基于非负约束2D变分模态分解的红外弱小目标检测方法与流程

一种基于非负约束2D变分模态分解的红外弱小目标检测方法，用于红外弱小目标检测，属于红外弱小目标检测技术领域。

背景技术：

红外成像技术在现代军事、天基探测、安防监控等领域得到了广泛的应用。其中，红外目标检测与跟踪系统是其关键技术和核心模块。由于多数应用中，红外成像距离远，导致目标往往只有几个像素到几十个像素，属于小目标。此外，由于红外探测系统的成像质量不高，以及大气系统对成像的干扰，目标的信杂比低，属于弱目标。红外弱小目标具有面积小、信杂比低、无形状及纹理特征以及快速运动等特点，其检测和跟踪具有很大难度。设计合适的红外弱小目标检测算法，可以大大提升红外目标检测与跟踪系统的有效性。

对于一幅含有红外弱小目标的图像而言，其低频成分主要由平滑或模糊的背景所构成。相反地，成像过程中的噪声构成了红外图像的高频成分。弱小目标由于面积小，缺乏纹理特征，一般处于红外图像平滑背景和高频噪声之间的中间频段。红外图像成分之间的这种频率差异性常常被用来进行图像去噪以及弱小目标检测。典型的方法包括基于小波变换的红外弱小目标检测方法，基于剪切波变换的检测方法以及基于轮廓波变换的方法等。这些方法利用了小波变换以及超小波变换的频率特性，在变换域中去除噪声及背景杂波。这些方法往往不能自适应地找到目标所在频带，难以去除与目标接近的频率成分，且容易受噪声干扰。

频域显著性也广泛地应用于红外弱小目标检测中。Qi等(2013)提出了一种基于四元组傅里叶变换的红外图像小目标增强方法，通过构建二阶方向导数滤波器，对输入图像进行四元组傅里叶变换，其相位谱包含了红外弱小目标的信息，同时抑制了背景杂波。此后，Han(2015)等提出了基于改进的四元组傅里叶变换的红外弱小目标检测方法，结合灰度及运动特性，在很多场景下都取得了良好的检测结果。这些方法依赖于红外弱小目标的运动特性及其高斯分布特性。当红外弱小目标的灰度分布不再接近高斯分布时，会影响检测效果。

在红外图像处理中，由于事先无法知道红外图像及弱小目标所处的准确频带，需要对红外图像的频率成分进行自适应估计。邓鹤、刘建国等(2010)将经验模态分解(Empirical Mode Decomposition,EMD)方法引入红外小目标图像处理中，进行图像的二维经验模态分解。经验模态分解(1998)是一种传统的信号分解方法，在每一次迭代过程中，重复地检测信号的局部最大/最小值，再通过插值的方式求得信号的上下包络。对上下包络求均值，可以得到当前信号的“中心信号”。用原信号减去该中心信号，所得到的差作为下一次迭代的输入。该方法在某些情况下，可以将信号分解为其本征分量的集合。但由于其缺少理论支撑，且算法效果极易受到极值点搜索结果的影响，该方法是不稳定的。在红外弱小目标检测中，该方法难以处理变化的背景，容易产生较大的虚警。

技术实现要素：

本发明针对上述不足之处提供了一种基于非负约束2D变分模态分解的红外弱小目标检测方法，解决现有技术中红外弱小目标所在的频带难以被准确估计、检测结果易受噪声及背景杂波干扰以及检测效率低的问题。

为了实现上述目的，本发明采用的技术方案如下：

一种基于非负约束2D变分模态分解的红外弱小目标检测方法，其特征在于，还包括以下步骤：

步骤1：获得一幅待处理的红外图像f(x,y)；

步骤2：采用带通滤波器进行红外图像f(x,y)预处理，即得到去除部分高频噪声及平滑的背景的红外图像f′(x,y)∈R^m×n；

步骤3：将二维变分模态分解方法与非负约束相结合，构建目标函数，再将步骤2得到的红外图像f′(x,y)∈R^m×n输入目标函数，并根据目标函数求解输出分解结果，即得到K个非负窄带子信号；

步骤4：提取步骤3得到的非负窄带子信号中的红外弱小目标的窄带子信号，忽略其余不含有红外弱小目标的窄带子信号，得到目标子信号；

步骤5：将步骤4中提取的目标子信号进行自适应阈值分割，确定红外弱小目标位置及大小，并输出检测结果。

进一步，所述步骤2中采用的带通滤波器为高斯差分滤波器或巴特沃斯带通滤波器。

进一步，所述步骤3的具体步骤如下：

(31)将二维变分模态分解方法与非负约束相结合，构建目标函数；具体步骤如下：

(311)假设输入的预处理后的红外图像f′(x,y)∈R^m×n将被分解为K个窄带子信号的集合，窄带子信号表示为u_k(x,y),k＝1,...,K；将空间域二维坐标(x,y)简记为x，计算窄带子信号u_k(x)的解析信号：其中δ表示二维脉冲信号，ω_k表示第k个子信号对应的频率域二维坐标，j为虚数单位，ω_k,⊥表示与ω_k相垂直的频率域二维坐标，*表示卷积运算，<·>代表内积运算；

(312)将解析信号进行频移，移动到预先估计出的图像基准频带φ_k，频移后的信号记为其中，u_AS,k(x)为解析信号，j为虚数单位，<·>代表内积运算；

(313)求解频移后的解析信号的带宽B_k，

(314)根据步骤(311)-步骤(313)，对窄带子信号u_k(x,y)进行非负约束以及数据保真约束，得到最终的目标函数，目标函数的公式如下：

其中，f′(x)表示步骤2得到的红外图像f′(x,y)∈R^m×n，u_k(x)≥0表示在最终的结果中，只有窄带子信号中的正数部分被视为有效信号，从而得以保留，负数部分被去掉，是确保分解出的窄带子信号不偏离预处理后的红外图像f′(x,y)∈R^m×n，ε为误差容许因子；

(32)将步骤2得到的预处理后的红外图像f′(x,y)∈R^m×n输入目标函数，采用交迭方向乘子方法(Alternative directional method of multipliers,ADMM)求解目标函数，并输出非负窄带子信号，其中，在第n+1次迭代过程中，窄带子信号的更新采用维纳滤波，傅里叶变换表示为图像基准频带的更新通过计算当前窄带子信号的功率谱重心实现；目标函数的具体求解过程如下：

(321)初始化迭代次数n＝0、窄带子信号图像基准频带傅里叶变换分别表示为以及其中k＝1,...,K；

(322)根据步骤(321)初始化的结果计算Hilbert模板：其中，表示当前的图像基准频带，ω表示频率域二维坐标，sgn表示符号函数；

(323)利用Hilbert模板，在频率域中计算第k个窄带子信号在第n+1次迭代中的解析信号其中表示Hilbert模板，ω表示频率域二维坐标，表示f′(x)的傅里叶变换结果，λ为拉格朗日乘子，α是计算中引入的权重因子；

(324)将解析信号从频率域反变换回空间域，公式如下：

(325)在反变换回空间域的窄带子信号中引入非负约束，得到非负窄带子信号，公式如下：

(326)更新迭代次数：n＝n+1；

(327)判断当前迭代次数是否大于最大迭代次数，如果是，则停止迭代，转到步骤(328)；如果不是，更新图像基准带值和频率域拉格朗日乘子返回步骤(322)；

(328)输出步骤(325)的分解结果，即非负窄带子信号u_k(x,y),k＝1,...,K。

进一步，所述步骤(327)中更新图像基准带值和频率域拉格朗日乘子的公式分别为：

其中，λ为拉格朗日乘子。

进一步，所述步骤4具体步骤如下：

(41)给出目标相似度因子(Target similarity factor,TSF)的公式：TSF_k＝(m_k-μ_k)/μ_k，m_k为u_k(x,y)的最大值，μ_k为u_k(x,y)的均值；

(42)将目标相似度因子的公式用于衡量当前非负窄带子信号为目标子信号的可能性，即具有最大的TSF的窄宽子信号被认为是含有目标的信号，其余窄宽子信号被认为是不含有目标的信息。

进一步，所述步骤5中，采用Otsu方法或最大熵分割方法对选择出的目标子信号进行分割，得到红外弱小目标位置及大小的估计，输出检测结果，即输出最终的红外弱小目标检测结果u(x,y)，其中，分割的公式如下：

其中，T为自适应阈值，u_t(x,y)为目标所在子信号，t∈R^K。

综上所述，由于采用了上述技术方案，本发明的有益效果是：

一、本发明将红外弱小目标检测问题转化为目标函数求解问题，自适应地将输入信号分解为数个非负窄带子信号的集合，可以高效、准确的找到目标所对应的窄带子信号，且不需要事先估计出目标所对应的频率；

二、本发明考虑了红外图像及红外目标的特点，利用了红外背景、目标及噪声之间的频率成分差异，同时考虑了目标成分的非负特性，去除计算过程中产生的干扰，同时抑制了背景杂波；

三、本发明具有较高的计算效率，迭代次数在5次以内，可以满足实时性要求。

附图说明

图1为本发明的流程图；

图2为一幅含有红外弱小目标的图像；

图3为带通滤波器传递函数示意图，其中，包括三维图及其二维投影；

图4为图2的带通滤波结果；

图5为图4所示的滤波结果经过本发明所述的方法处理后得到的结果图，及结果图所对应的三维投影图；

图6为本发明中选择出的目标子信号经过自适应阈值分割的结果。

具体实施方式

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

一种基于非负约束2D变分模态分解的红外弱小目标检测方法，先采用带通滤波器对红外图像进行预处理，再结合红外图像及目标的特点，对预处理后的红外图像进行二维变分模态分解，对待分解子信号施加非负约束，将窄带子信号分解为窄带信号的集合，得到K个窄带子信号。选择目标子信号，对目标子信号进行自适应阈值分割，得到红外弱小目标检测结果。具体如下：

一种基于非负约束2D变分模态分解的红外弱小目标检测方法，还包括以下步骤：

步骤1：输入一幅待处理的红外图像f(x,y)，其大小为m×n，如图2所示，其中，方框所示即为弱小目标所在位置；

步骤2：采用带通滤波器进行红外图像f(x,y)预处理，即得到去除部分高频噪声及平滑的背景的红外图像f′(x,y)∈R^m×n，带通滤波器可以为高斯差分滤波器或巴特沃斯带通滤波器，如图3、图4所示，；

步骤3：将二维变分模态分解方法与非负约束相结合，构建目标函数，再将步骤2得到的红外图像f′(x,y)∈R^m×n输入目标函数，并根据目标函数求解输出分解结果，即得到K个非负窄带子信号。具体步骤如下：

(31)将二维变分模态分解方法与非负约束相结合，构建目标函数；具体步骤如下：

(311)假设输入的预处理后的红外图像f′(x,y)∈R^m×n将被分解为K个窄带子信号的集合，窄带子信号表示为u_k(x,y),k＝1,...,K；将空间域二维坐标(x,y)简记为x，计算窄带子信号u_k(x)的解析信号：其中δ表示二维脉冲信号，ω_k表示第k个子信号对应的频率域二维坐标，j为虚数单位，ω_k,⊥表示与ω_k相垂直的频率域二维坐标，*表示卷积运算，<·>代表内积运算；

(312)将解析信号进行频移，移动到预先估计出的图像基准频带φ_k，频移后的信号记为u_shift(x)，其中，u_AS,k(x)为解析信号，j为虚数单位，<·>代表内积运算；

(313)求解频移后的解析信号的带宽B_k，

(314)根据步骤(311)-步骤(313)，对窄带子信号u_k(x,y)进行非负约束以及数据保真约束，得到最终的目标函数，目标函数的公式如下：

其中，f′(x)表示步骤2得到的红外图像f′(x,y)∈R^m×n，u_k(x)≥0表示在最终的结果中，只有窄带子信号中的正数部分被视为有效信号，从而得以保留，负数部分被去掉，是确保分解出的窄带子信号不偏离预处理后的红外图像f′(x,y)∈R^m×n，ε为误差容许因子；

(32)将步骤2得到的预处理后的红外图像f′(x,y)∈R^m×n输入目标函数，采用交迭方向乘子方法(Alternative directional method of multipliers,ADMM)求解目标函数，并输出非负窄带子信号，其中，在第n+1次迭代过程中，窄带子信号的更新采用维纳滤波，傅里叶变换表示为图像基准频带的更新通过计算当前窄带子信号的功率谱重心实现，如图5所示，图5中从左到右，从上到下依次为k＝1，k＝2以及k＝3的窄带子信号分量，及其对应的三维投影，可以看出，在该场景下，目标分量主要包含在k＝1的分量中，k＝2以及k＝3的分量主要包含了信号的震荡以及方向信息；目标函数的具体求解过程如下：

(321)初始化迭代次数n＝0、窄带子信号图像基准频带傅里叶变换分别表示为以及其中k＝1,...,K；

(322)根据步骤(321)初始化的结果计算Hilbert模板：其中，表示当前的图像基准频带，ω表示频率域二维坐标，sgn表示符号函数；

(323)利用Hilbert模板，在频率域中计算第k个窄带子信号在第n+1次迭代中的解析信号其中表示Hilbert模板，ω表示频率域二维坐标，表示f′(x)的傅里叶变换结果，λ为拉格朗日乘子，α是计算中引入的权重因子；

(324)将解析信号从频率域反变换回空间域，公式如下：

(325)在反变换回空间域的窄带子信号中引入非负约束，得到非负窄带子信号，公式如下：

(326)更新迭代次数：n＝n+1；

(327)判断当前迭代次数是否大于最大迭代次数，如果是，则停止迭代，转到步骤(328)；如果不是，更新图像基准带值和频率域拉格朗日乘子返回步骤(322)；其中，更新图像基准带值和频率域拉格朗日乘子的公式分别为：

其中，λ为拉格朗日乘子。

(328)输出步骤(325)的分解结果，即非负窄带子信号u_k(x,y),k＝1,...,K。

步骤4：提取非负窄带子信号对应于红外弱小目标的某个窄带子信号，忽略其余不含有目标信息的窄带子信号，得到目标子信号；具体步骤如下：

(41)给出目标相似度因子(Target similarity factor,TSF)的公式：TSF_k＝(m_k-μ_k)/μ_k，其中，m_k为u_k(x,y)的最大值，μ_k为u_k(x,y)的均值，计算所有窄带子信号的TSF值，具有最大的TSF值的子信号被认为是目标信号；

(42)将目标相似度因子的公式用于衡量当前非负窄带子信号为目标子信号的可能性，即具有最大的TSF的窄宽子信号被认为是含有目标的信号，其余窄宽子信号被认为是不含有目标的信息。

步骤5：将步骤4中提取的目标子信号进行自适应阈值分割，确定目标子信号位置及大小，并输出检测结果。采用Otsu方法或最大熵分割方法对选择出的目标子信号进行分割，得到红外弱小目标位置及红外弱小目标大小的估计，输出检测结果，即输出最终的红外弱小目标检测结果u(x,y)，如图6所示，假设目标所在的窄带子信号为u_t(x,y),t∈R^K，计算得到的自适应阈值为T，则分割的公式如下：

其中，T为自适应阈值，u_t(x,y)为目标所在子信号，t∈R^K，R代表实数。进行自适应阈值处理，可进一步抑制杂波。

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