CN104715458A - 一种双模非局部均值滤波方法 - Google Patents

Abstract

本发明公开了一种双模非局部均值滤波方法，该方法首先利用目标先验知识确定滤波模型参数，根据图像中待处理像素，定位相似窗以及下一帧图像中的干涉区域，然后遍历干涉区域像素，定位匹配窗，对于每一个匹配窗，首先提取相似窗向量和匹配窗向量，分别计算相似窗和匹配窗的目标模型相似度以及背景模型相似度，然后利用目标模型相似度和背景模型相似度计算高斯模型和复合指数模型权值，根据干涉区域内所有双模权值得到相似窗中心像素的背景估计值，最后用原始图像和估计背景图像作差，获取背景抑制结果。本发明有效解决了由于复杂背景造成的强杂波干扰问题，并克服弱小目标造成的有效的几何轮廓特征缺乏问题，适用于复杂背景弱小运动目标的提取。

Description

一种双模非局部均值滤波方法

技术领域

本发明属于图像处理与目标检测技术领域，更具体地，是一种用于弱小运动目标提取的双模非局部均值滤波方法。

背景技术

在复杂背景且信噪比较低的情况下，图像序列中弱小运动目标的检测一直是图像处理领域中比较困难的问题。当目标处于远距离时，它们在像面上仅占几个像素，无形状与结构特征，且强度较弱，并处于复杂背景下，被看作是复杂环境下的弱小目标。弱小运动目标处于复杂背景中时，由于目标和背景的双重复杂性，使得目标检测的难度大为增加。因此，对复杂背景下弱小运动目标的检测技术的研究，对于提高现代化高技术武器的作战距离及反应速度具有重要意义。

目标检测算法按照目标增强和背景估计两个方面的思路进行分类。目标增强的本质是充分利用目标本身所固有的、规律性的多种特性，使变换后的目标特性得到有效的增强，改进目标特性的显著性和与背景的分离性；而背景估计的本质是充分利用背景本身所固有的、规律性的变化特性，使变换后的背景特性具有某种相对的稳定性，改进背景估计的精度和与目标特性的分离性。在目标的可利用特征有限的情况下，背景估计方法更有优势。图像中的任何一个像素都不是孤立的，而是与其周围的像素点结合在一起，共同构成图像的几何结构。以某一像素点为中心的窗口邻域(或称图像块)，可以很好地描述像素点的结构特征。任何一个像素点的图像块的所有集合，可以看作是图像的一种过完备表示。另外，图像还具有自相似性质，也就是说，处于图像中不同位置上的像素点常常表现出明显的强相关性。

对于背景估计算法，国内外学者已做了很多富有成效的工作。简单的时空域滤波方法如高通滤波，中值滤波，匹配滤波等无法满足复杂背景下弱小目标的提取需求，基于小波变换的滤波方法对小波基的选取要求较高，简单的小波基对弱小目标提取的效果并不理想，而复杂的小波基又在计算复杂度和计算量上有待改进。固定权值背景预测算法是最简单的线性低通滤波背景预测算法，计算简单快捷，但不适用于预测强起伏的背景；基于形态学滤波的背景预测利用形态学中腐蚀与膨胀操作的对称性可以很好地估计出背景。膨胀运算使得在比结构元素小的区域中的黑色细节将被减少或去除，腐蚀运算使得在比结构元素小的区域中的明亮细节经腐蚀后将被减少或去除，其结果与高通滤波器是相似的，虽然在抑制噪声方面有一些优点，但对低信噪比图像的小目标检测性能有限，而且形态学滤波中对结构元的选取要求很高。非局部均值滤波方法是由Buades提出的一种基于单一模型的空域滤波方法，利用对含噪声的图像的背景估计得到去噪图像，针对高斯噪声非常有效，但用于弱小目标检测中的背景抑制会导致目标丢失。

基于背景估计原理的目标检测算法，其难点在于需要在目标存在的情况下精确地估计出不含目标的背景，由于背景估计过程建立在原始图像之上，往往会受到目标的影响而无法估计出正确背景，导致目标残留，得出的残差图像中目标丢失，检测率降低。而且预警系统中的目标检测算法受到星上运算系统的限制，计算量和存储量都有一定要求，因此实时鲁棒的小目标自动检测技术尚未完全突破。

发明内容

针对现有技术的以上缺陷或改进需求，本发明提供了一种双模非局部均值滤波方法，能对图像原始背景进行估计，从而对背景抑制结果进行目标提取，有效地解决了由于复杂背景造成的强杂波干扰以及弱小目标可利用特征有限的问题，提高了复杂背景下弱小目标的检测概率。

为了实现上述目的，本发明提供了一种双模非局部均值滤波方法，包括如下步骤：

(1)图像数据输入：

输入待处理的图像序列；

(2)滤波模型参数确定：

根据目标先验知识确定滤波模型参数；

(3)相似窗和干涉区域定位：

根据图像中的待处理像素为中心划定相似窗，并在下一帧图像中确定干涉区域范围；

(4)匹配向量提取：

以干涉区域中当前像素确定匹配窗，根据(3)中的相似窗和当前匹配窗的像素灰度值得到相似窗向量和匹配窗向量；

(5)背景模型和目标模型相似度测量：

利用步骤(4)中的相似窗向量和匹配窗向量进行背景模型和目标模型的相似度测量。

(6)双模权值计算

根据步骤(5)中的背景模型和目标模型相似度测量值计算双模权值。

(7)背景估计值计算：

判断当前处理的匹配窗是否是以干涉区域内最后一个像素为中心的匹配窗，如果满足条件，则利用步骤(4)-(6)得到的所有匹配窗中心像素值和该匹配窗的权值获得待处理像素对应的估计值然后进行步骤(8)，否则返回步骤(4)继续处理干涉区域下一个像素；

(8)残差获取：

判断当前处理像素是否为图像中最后一个像素，如果满足条件，则将所得估计背景图像与原始图像作差，得到目标得到增强而背景得到抑制的图像，否则转到步骤(3)处理图像中的下一像素。

在本发明的一个实施例中，在步骤(2)滤波模型参数确定步骤中，根据输入的目标尺寸的先验知识确定滤波模型的尺寸，滤波模版分别为相似窗，匹配窗和干涉区域，相似窗和匹配窗尺寸相同，都为R×R，干涉区域大小为L×L，且满足R≤L＜2R，窗口模版分为目标模型和背景模型区域，目标模型为以待处理像素为中心的中央小型窗r×r大小的区域B，背景模型区域为剩下部分组成的环形区域A，目标尺寸若不超过t×t像素块大小，则r,t,R的关系满足R＞r＝t，一般取值为R＝r+2,L＝R+2。

目标模型权值采用高斯模型ω＝exp(-DIS_B/h₂)计算权值，背景模型采用复合指数模型ω＝exp(h₁/(DIS_A+1))计算权值。其中DIS_A和DIS_B分别为背景模型和目标模型的相似度，h₁,h₂为滤波参数，50＜h₂＜150，且需要满足

$\left\{\begin{array}{c}{h}_1{h}_2>({\mathrm{\&sigma;}}_1^2+1)({\mathrm{\&sigma;}}_1^2+{\mathrm{\&sigma;}}_3^2)\\ h_1{h}_2>\frac{(2{\mathrm{\&sigma;}}_3^2-{\mathrm{\&sigma;}}_2^2)({\mathrm{\&sigma;}}_1^2+1)({\mathrm{\&sigma;}}_1^2+{\mathrm{\&sigma;}}_2^2+1)}{{\mathrm{\&sigma;}}_2^2}\\ h_1{h}_2<({\mathrm{\&sigma;}}_1^2+{\mathrm{\&sigma;}}_2^2)({\mathrm{\&sigma;}}_1^2+{\mathrm{\&sigma;}}_2^2+1)\end{array}\right.$

其中，为噪声标准差，为背景标准差，为目标标准差。

在本发明的一个实施例中，在步骤(3)相似窗和干涉区域定位步骤中，相似窗和干涉区域的位置和大小按照根据如下方法确定：

相似窗以待处理像素为中心，选取R×R大小的图像块，干涉区域以下一帧图像中与待处理像素同一位置的像素为中心选取L×L大小的图像块。

在本发明的一个实施例中，所述步骤(4)具体包括以下子步骤：

(4-1)以当前像素为中心选取R×R大小的图像块为匹配窗；

(4-2)将以像素i为中心的相似窗N_i划分为两个区域，分别是以待处理像素为中心的中央小型窗r×r大小的区域B(目标模型区域)和剩下部分组成的环形区域A(背景模型区域)，对以像素j为中心的匹配窗N_j进行同样处理；

(4-3)将相似窗N_i与匹配窗N_j中的区域A中的像素灰度值按同样的先后顺序分别组成向量v(N_iA)和v(N_jA)；

(4-4)将相似窗N_i与匹配窗N_j中的区域B中的像素灰度值按同样的先后顺序分别组成向量v(N_iB)和v(N_jB)；

在本发明的一个实施例中，所述步骤(5)具体包括以下子步骤：

(5-1)根据(4-3)中得到的相似窗与匹配窗区域A对应的向量v(N_iA)和v(N_jA)按照如下公式求取两个向量的欧氏距离，得到背景模型相似度DIS_A：

${\mathrm{DIS}}_A={\left|\right|v\left(N_{\mathrm{iA}}\right)-v\left(N_{\mathrm{jA}}\right)\left|\right|}_2^2$

(5-2)根据(4-4)中得到的相似窗与匹配窗区域B对应的向量v(N_iB)和v(N_jB)，并按照如下公式求取两个向量的欧氏距离，得到目标模型相似度DIS_B：

${\mathrm{DIS}}_B={\left|\right|v\left(N_{\mathrm{iB}}\right)-v\left(N_{\mathrm{jB}}\right)\left|\right|}_2^2$

在本发明的一个实施例中，所述步骤(6)具体包括以下子步骤：

(6-1)利用相似窗N_i与匹配窗N_j对应的背景模型相似度DIS_A和目标模型相似度DIS_B根据如下公式来计算与匹配窗中心像素j对应的权值：

$\mathrm{\&omega;}(i,j)=\exp (\frac{h_1}{{\mathrm{DIS}}_A+1}-\frac{{\mathrm{DIS}}_B}{h_2})$

其中,DIS_A和DIS_B分别为背景模型和目标模型的相似度，h₁,h₂为滤波参数；

(6-2)将匹配窗N_j对应的权值累加至相似窗N_i对应的所有匹配窗权值累加和中，用于权值归一化计算；

在本发明的一个实施例中，所述步骤(7)具体包括以下子步骤：

(7-1)判断当前处理的匹配窗是否为像素i的干涉区域中最后一个像素对应的匹配窗，如果满足条件，则按(7-2)所述计算背景估计值，否则继续处理干涉区域中下一像素；

(7-2)按照如下公式计算待处理像素i的背景估计值

$\tilde{f}\left(i\right)=\underset{j\&Element;I}{\mathrm{\&Sigma;}}{\mathrm{\&omega;}}^{\&prime;}(i,j)f\left(j\right)=\underset{j\&Element;I}{\mathrm{\&Sigma;}}\frac{\mathrm{\&omega;}(i,j)}{Z\left(i\right)}f\left(j\right)$

$Z\left(i\right)=\underset{j\&Element;I}{\mathrm{\&Sigma;}}\mathrm{\&omega;}(i,j)$

其中,I为待处理像素i对应的干涉区域，j为干涉区域I中所有匹配窗中心像素，f(j)为其像素灰度值，ω(i,j)为匹配窗N_j对应的权值，Z(i)为相似窗N_i的干涉区域内所有匹配窗权值累加和。

总体而言，通过本发明所构思的以上技术方案与现有技术相比，能够取得下列有益效果：

由于将相似窗和匹配窗分为背景模型区域和目标模型区域，分别统计背景模型相似度和目标模型相似度，并利用高斯模型和复合指数模型进行双模权值计算，使得背景估计对目标的存在不敏感，能有效估计出原始背景像素值，抑制复杂背景造成的强杂波干扰，有利于目标提取，对于弱小目标仍然具有鲁棒性，有效提高了弱小目标检测概率，降低了虚警率。

附图说明

图1为本发明双模非局部均值滤波方法的处理流程图；

图2为本发明实施例中相似窗和干涉区域定位示意图；

图3为本发明实施例中滤波参数曲线示意图；

图4为本发明实施例中匹配窗定位示意图；

图5为本发明实施例中相似窗和匹配窗分区示意图；

图6为本发明实施例中像素转化向量示意图；

图7为本发明实施例中待处理原始图像第一帧；

图8为本发明实施例中待处理原始图像第二帧；

图9为本发明实施例中背景估计结果图像；

图10为本发明实施例中背景抑制结果图像；

图11为本发明实施例中背景抑制图像三维图像；

图12为本发明实施例中SNR_in＜2时各算法性能比较ROC曲线；

图13为本发明实施例中2＜SNR_in＜5时各算法性能比较ROC曲线；

图14为本发明实施例中SNR_in＞5时各算法性能比较ROC曲线；

具体实施方式

为了使本发明的目的、技术方案及优点更加清楚明白，以下结合附图及实施例，对本发明进行进一步详细说明。应当理解，此处所描述的具体实施例仅用以解释本发明，并不用于限定本发明。此外，下面所描述的本发明各个实施方式中所涉及到的技术特征只要彼此之间未构成冲突就可以相互组合。

本发明利用图像的自相关性，通过对像素块及其邻域进行目标模型相似度和背景模型相似度测量，根据相似度测量值计算像素块对应的双模权值，对邻域像素进行加权平均，用该平均值估计像素灰度值，从而得到图像的背景估计图像，通过与原始图像作差使背景杂波得到抑制。具体为：遍历干涉区域根据每一个像素依次定位匹配窗，对于每一个匹配窗，首先提取相似窗向量和匹配窗向量，根据相似窗向量和匹配窗向量分别计算相似窗和匹配窗的目标模型相似度以及背景模型相似度，然后利用目标模型相似度测量值和背景模型相似度测量值进行高斯模型和复合指数模型权值计算，根据干涉区域内所有双模权值得到相似窗中心像素的背景估计值，最后用原始图像和估计背景图像作差，获取背景抑制结果。

如图1所示，本发明一种双模非局部均值滤波方法，包括以下步骤：

(1)输入图像序列

本发明适用于复杂背景下弱小运动目标的检测，尤其是低信噪比输入图像。本实施例中待处理图像大小为M×N＝256×256，叠加高斯噪声参数大小为μ＝0,σ＝2。由于相似度大小受到目标位置的影响，目标运动速度如果大于2像素/帧，算法性能会比目标运动速度小于2像素/帧的图像序列更好。

(2)确定滤波参数大小

高斯点扩展标准差σ＝0.3，在此条件下，弱小目标在相面上不超过3×3像素块大小。相似窗和干涉区域的尺寸确定方法为：根据关系式r＝t,R＝r+2,L＝R+2。令r＝3,R＝5,L＝7。其中r为区域B(目标模型区域)尺寸，R为相似窗尺寸，L为干涉区域尺寸。根据像素分别在相邻两帧图像中确定相似窗和干涉区域，如图2所示。对于滤波参数h₁,h₂的确定，推导如下：

权值计算公式为

$\mathrm{\&omega;}(i,j)=\mathrm{eap}(\frac{h_1}{{\mathrm{DIS}}_A+1}-\frac{{\mathrm{DIS}}_B}{h_2})---\left(1\right)$

其中，DIS_A和DIS_B分别为背景模型和目标模型的相似度测量值，或称为不相似因子，h₁,h₂为滤波参数；

令为噪声方差，为背景方差，为目标方差，并将式(1)简化为讨论幂次方的大小，即

$\mathrm{\&omega;}=\frac{h_1}{{\mathrm{DIS}}_A+1}-\frac{{\mathrm{DIS}}_B}{h_2}---\left(2\right)$

DIS_A可以近似为相似窗和匹配窗中背景模型对应像素值的方差，DIS_B可以近似为相似窗和匹配窗中目标模型对应像素值的方差，由于目标是运动的，不可能出现目标周围背景相同的情况，将其排除之后，权值计算时会出现以下几种情况：

①相似窗和匹配窗背景相同，一个其中有目标一个没有，则两者背景模型中不相似因子为噪声，目标模型中不相似因子为目标和噪声，那么权值计算如下：

${\mathrm{\&omega;}}_1=\frac{h_1}{{\mathrm{\&sigma;}}_1^2+1}-\frac{{\mathrm{\&sigma;}}_1^2+{\mathrm{\&sigma;}}_3^2}{h_2}---\left(3\right)$

②相似窗和匹配窗背景相同，两者都不含目标，则两者背景模型中不相似因子为噪声，目标模型中不相似因子为噪声，那么权值计算如下：

${\mathrm{\&omega;}}_2=\frac{h_1}{{\mathrm{\&sigma;}}_1^2+1}-\frac{{\mathrm{\&sigma;}}_1^2}{h_2}---\left(4\right)$

③相似窗和匹配窗背景不同，一个其中有目标一个没有，则两者背景模型中不相似因子为背景和噪声，目标模型中不相似因子为目标、背景和噪声，那么权值计算如下：

${\mathrm{\&omega;}}_3=\frac{h_1}{{\mathrm{\&sigma;}}_1^2+{\mathrm{\&sigma;}}_2^2+1}-\frac{{\mathrm{\&sigma;}}_1^2+{\mathrm{\&sigma;}}_2^2+{\mathrm{\&sigma;}}_3^2}{h_2}---\left(5\right)$

④相似窗和匹配窗背景不同，两者都不含目标，则两者背景模型中不相似因子为背景和噪声，目标模型中不相似因子为背景和噪声，那么权值计算如下：

${\mathrm{\&omega;}}_4=\frac{h_1}{{\mathrm{\&sigma;}}_1^2+{\mathrm{\&sigma;}}_2^2+1}-\frac{{\mathrm{\&sigma;}}_1^2+{\mathrm{\&sigma;}}_2^2}{h_2}---\left(6\right)$

⑤相似窗和匹配窗背景不同，两者都含有目标，则两者背景模型中不相似因子为背景和噪声，目标模型中不相似因子为背景和噪声，那么权值计算如下：

${\mathrm{\&omega;}}_4=\frac{h_1}{{\mathrm{\&sigma;}}_1^2+{\mathrm{\&sigma;}}_2^2+1}-\frac{{\mathrm{\&sigma;}}_1^2+{\mathrm{\&sigma;}}_2^2}{h_2}---\left(7\right)$

由于④⑤都为不匹配情况，且权值计算结果相同，因此将其合并讨论。

由于同一相似窗对应的所有匹配窗的权值需要累加进行归一化计算，而最终估计值的结果是由匹配窗权值与匹配窗中心像素值的乘积累加得到，每一个匹配窗权值大小决定了该匹配窗中心像素值占最后的估计结果的比例大小，因此正确匹配的情况下权值应远大于不匹配情况下的权值大小，为了达到更好的背景抑制效果，我们有如下期望：

(1)对于ω的绝对大小，我们希望尽量有①、②这种正确匹配的情况对应的权值ω₁,ω₂＞0，而③、④、⑤这种不匹配的情况对应的权值ω₃,ω₄＜0。

(2)对于ω的相对大小，我们希望对于①、②这种正确匹配的情况对应的权值远大于③、④、⑤这种不匹配的情况对应的权值

(3)对于ω的相对大小，我们希望对于①，②这两种正确匹配的情况对应的权值之差远大于它们与③、④、⑤这种不匹配的情况对应的权值之差。

即拉开匹配与不匹配情况下对应权值的范围的距离大小，可看做(2)的的补充。

(1)由 ${\mathrm{\&omega;}}_1=\frac{h_1}{{\mathrm{\&sigma;}}_1^2+1}-\frac{{\mathrm{\&sigma;}}_1^2+{\mathrm{\&sigma;}}_3^2}{h_2}>0$ 得到

$h_1{h}_2>({\mathrm{\&sigma;}}_1^2+1)({\mathrm{\&sigma;}}_1^2+{\mathrm{\&sigma;}}_3^2)---\left(8\right)$

由 ${\mathrm{\&omega;}}_2=\frac{h_1}{{\mathrm{\&sigma;}}_1^2+1}-\frac{{\mathrm{\&sigma;}}_1^2}{h_2}>0$ 得到

$h_1{h}_2>{\mathrm{\&sigma;}}_1^2({\mathrm{\&sigma;}}_1^2+1)---\left(9\right)$

由 ${\mathrm{\&omega;}}_3=\frac{h_1}{{\mathrm{\&sigma;}}_1^2+{\mathrm{\&sigma;}}_2^2+1}-\frac{{\mathrm{\&sigma;}}_1^2+{\mathrm{\&sigma;}}_2^2+{\mathrm{\&sigma;}}_3^2}{h_2}<0$ 得到

$h_1{h}_2<({\mathrm{\&sigma;}}_1^2+{\mathrm{\&sigma;}}_2^2+1)({\mathrm{\&sigma;}}_1^2+{\mathrm{\&sigma;}}_2^2+{\mathrm{\&sigma;}}_3^2)---\left(10\right)$

由 ${\mathrm{\&omega;}}_4=\frac{h_1}{{\mathrm{\&sigma;}}_1^2+{\mathrm{\&sigma;}}_2^2+1}-\frac{{\mathrm{\&sigma;}}_1^2+{\mathrm{\&sigma;}}_2^2}{h_2}<0$ 得到

$h_1{h}_2<({\mathrm{\&sigma;}}_1^2+{\mathrm{\&sigma;}}_2^2)({\mathrm{\&sigma;}}_1^2+{\mathrm{\&sigma;}}_2^2+1)---\left(11\right)$

由于 $({\mathrm{\&sigma;}}_1^2+1)({\mathrm{\&sigma;}}_1^2+{\mathrm{\&sigma;}}_3^2)>{\mathrm{\&sigma;}}_1^2({\mathrm{\&sigma;}}_1^2+1),({\mathrm{\&sigma;}}_1^2+{\mathrm{\&sigma;}}_2^2)({\mathrm{\&sigma;}}_1^2+{\mathrm{\&sigma;}}_2^2+1)<({\mathrm{\&sigma;}}_1^2+{\mathrm{\&sigma;}}_2^2+1)({\mathrm{\&sigma;}}_1^2+{\mathrm{\&sigma;}}_2^2+{\mathrm{\&sigma;}}_3^2)$ 要使式(8)-(11)成立，需有式(8)、式(11)成立，即

$({\mathrm{\&sigma;}}_1^2+1)({\mathrm{\&sigma;}}_1^2+{\mathrm{\&sigma;}}_3^2)<({\mathrm{\&sigma;}}_1^2+{\mathrm{\&sigma;}}_2^2)({\mathrm{\&sigma;}}_1^2+{\mathrm{\&sigma;}}_2^2+1)---\left(12\right)$

化简为

${\mathrm{\&sigma;}}_3^2<\frac{2{\mathrm{\&sigma;}}_1^2+{\mathrm{\&sigma;}}_2^2+1}{{\mathrm{\&sigma;}}_1^2+1}{\mathrm{\&sigma;}}_2^2---\left(13\right)$

(2)要使ω₁,ω₂＞＞ω₃,ω₄，由于

${\mathrm{\&omega;}}_3=\frac{h_1}{{\mathrm{\&sigma;}}_1^2+{\mathrm{\&sigma;}}_2^2+1}-\frac{{\mathrm{\&sigma;}}_1^2+{\mathrm{\&sigma;}}_2^2+{\mathrm{\&sigma;}}_3^2}{h_2}<{\mathrm{\&omega;}}_4=\frac{h_1}{{\mathrm{\&sigma;}}_1^2+{\mathrm{\&sigma;}}_2^2+1}-\frac{{\mathrm{\&sigma;}}_1^2+{\mathrm{\&sigma;}}_2^2}{h_2},$ 只需ω₁,ω₂＞＞ω₄

$\begin{array}{c}{\mathrm{\&omega;}}_1-{\mathrm{\&omega;}}_4=\frac{h_1}{{\mathrm{\&sigma;}}_1^2+1}-\frac{{\mathrm{\&sigma;}}_1^2+{\mathrm{\&sigma;}}_3^2}{h_2}-\frac{h_1}{{\mathrm{\&sigma;}}_1^2+{\mathrm{\&sigma;}}_2^2+1}+\frac{{\mathrm{\&sigma;}}_1^2+{\mathrm{\&sigma;}}_2^2}{h_2}\\ =\frac{({\mathrm{\&sigma;}}_1^2+{\mathrm{\&sigma;}}_2^2+1)h_1-({\mathrm{\&sigma;}}_1^2+1)h_1}{({\mathrm{\&sigma;}}_1^2+1)({\mathrm{\&sigma;}}_1^2+{\mathrm{\&sigma;}}_2^2+1)}+\frac{{\mathrm{\&sigma;}}_2^2-{\mathrm{\&sigma;}}_3^2}{h_2}\\ =\frac{{\mathrm{\&sigma;}}_2^2{h}_1}{({\mathrm{\&sigma;}}_1^2+1)({\mathrm{\&sigma;}}_1^2+{\mathrm{\&sigma;}}_2^2+1)}+\frac{{\mathrm{\&sigma;}}_2^2-{\mathrm{\&sigma;}}_3^2}{h_2}\\ ={\mathrm{\&sigma;}}_2^2(\frac{h_1}{({\mathrm{\&sigma;}}_1^2+1)({\mathrm{\&sigma;}}_1^2+{\mathrm{\&sigma;}}_2^2+1)}+\frac{1}{h_2})-\frac{{\mathrm{\&sigma;}}_3^2}{h_2}\\ ={\mathrm{\&sigma;}}_2^2\left(\frac{h_1}{({\mathrm{\&sigma;}}_1^2+1)({\mathrm{\&sigma;}}_1^2+{\mathrm{\&sigma;}}_2^2+1)}\right)-\frac{{\mathrm{\&sigma;}}_3^2-{\mathrm{\&sigma;}}_2^2}{h_2}\end{array}---\left(14\right)$

$\begin{array}{c}{\mathrm{\&omega;}}_2-{\mathrm{\&omega;}}_4=\frac{h_1}{{\mathrm{\&sigma;}}_1^2+1}-\frac{{\mathrm{\&sigma;}}_1^2}{h_2}-\frac{h_1}{{\mathrm{\&sigma;}}_1^2+{\mathrm{\&sigma;}}_2^2+1}+\frac{{\mathrm{\&sigma;}}_1^2+{\mathrm{\&sigma;}}_2^2}{h_2}\\ =\frac{({\mathrm{\&sigma;}}_1^2+{\mathrm{\&sigma;}}_2^2+1)h_1-({\mathrm{\&sigma;}}_1^2+1)h_1}{({\mathrm{\&sigma;}}_1^2+1)({\mathrm{\&sigma;}}_1^2+{\mathrm{\&sigma;}}_2^2+1)}+\frac{{\mathrm{\&sigma;}}_2^2}{h_2}\\ =\frac{{\mathrm{\&sigma;}}_2^2{h}_1}{({\mathrm{\&sigma;}}_1^2+1)({\mathrm{\&sigma;}}_1^2+{\mathrm{\&sigma;}}_2^2+1)}+\frac{{\mathrm{\&sigma;}}_2^2}{h_2}\\ ={\mathrm{\&sigma;}}_2^2(\frac{h_1}{({\mathrm{\&sigma;}}_1^2+1)({\mathrm{\&sigma;}}_1^2+{\mathrm{\&sigma;}}_2^2+1)}+\frac{1}{h_2})\end{array}---\left(15\right)$

在 $t=\frac{{\mathrm{\&sigma;}}_2^2{h}_1}{({\mathrm{\&sigma;}}_1^2+1)({\mathrm{\&sigma;}}_1^2+{\mathrm{\&sigma;}}_2^2+1)}$ 一定的情况下，有

${\mathrm{\&omega;}}_1-{\mathrm{\&omega;}}_4=t-\frac{{\mathrm{\&sigma;}}_3^2-{\mathrm{\&sigma;}}_2^2}{h_2}---\left(16\right)$

${\mathrm{\&omega;}}_2-{\mathrm{\&omega;}}_4=t+\frac{{\mathrm{\&sigma;}}_2^2}{h_2}---\left(17\right)$

(3)由于 ${\mathrm{\&sigma;}}_3^2>0,$ 因此 ${\mathrm{\&omega;}}_1=\frac{h_1}{{\mathrm{\&sigma;}}_1^2+1}-\frac{{\mathrm{\&sigma;}}_1^2+{\mathrm{\&sigma;}}_3^2}{h_2}<{\mathrm{\&omega;}}_2=\frac{h_1}{{\mathrm{\&sigma;}}_1^2+1}-\frac{{\mathrm{\&sigma;}}_1^2}{h_2},$ 且有

${\mathrm{\&omega;}}_3=\frac{h_1}{{\mathrm{\&sigma;}}_1^2+{\mathrm{\&sigma;}}_3^2+1}-\frac{{\mathrm{\&sigma;}}_1^2+{\mathrm{\&sigma;}}_2^2+{\mathrm{\&sigma;}}_3^2}{h_2}<{\mathrm{\&omega;}}_4=\frac{h_1}{{\mathrm{\&sigma;}}_1^2+{\mathrm{\&sigma;}}_2^2+1}-\frac{{\mathrm{\&sigma;}}_1^2+{\mathrm{\&sigma;}}_2^2}{h_2},$ 则只需ω₂-ω₁＜＜ω₁-ω₄

即 $\frac{{\mathrm{\&sigma;}}_3^2}{h_2}<{\mathrm{\&sigma;}}_2^2\left(\frac{h_1}{({\mathrm{\&sigma;}}_1^2+1)({\mathrm{\&sigma;}}_1^2+{\mathrm{\&sigma;}}_2^2+1)}\right)-\frac{{\mathrm{\&sigma;}}_3^2-{\mathrm{\&sigma;}}_2^2}{h_2}$

化简为

$h_1{h}_2>\frac{({2\mathrm{\&sigma;}}_3^2-{\mathrm{\&sigma;}}_2^2)({\mathrm{\&sigma;}}_1^2+1)({\mathrm{\&sigma;}}_1^2+{\mathrm{\&sigma;}}_2^2+1)}{{\mathrm{\&sigma;}}_2^2}---\left(18\right)$

结合式(11)可得 $\frac{({2\mathrm{\&sigma;}}_3^2-{\mathrm{\&sigma;}}_2^2)({\mathrm{\&sigma;}}_1^2+1)({\mathrm{\&sigma;}}_1^2+{\mathrm{\&sigma;}}_2^2+1)}{{\mathrm{\&sigma;}}_2^2}<({\mathrm{\&sigma;}}_1^2+{\mathrm{\&sigma;}}_2^2)({\mathrm{\&sigma;}}_1^2+{\mathrm{\&sigma;}}_2^2+1),$ 化简得

${\mathrm{\&sigma;}}_3^2<\frac{{\mathrm{\&sigma;}}_2^2({\mathrm{\&sigma;}}_1^2+{\mathrm{\&sigma;}}_2^2)}{2({\mathrm{\&sigma;}}_1^2+1)}+\frac{{\mathrm{\&sigma;}}_2^2}{2}=\frac{{2\mathrm{\&sigma;}}_1^2+{\mathrm{\&sigma;}}_2^2+1}{2}\&CenterDot;\frac{{\mathrm{\&sigma;}}_2^2}{{\mathrm{\&sigma;}}_1^2+1}---\left(19\right)$

而式(13)有 ${\mathrm{\&sigma;}}_3^2<\frac{{2\mathrm{\&sigma;}}_1^2+{\mathrm{\&sigma;}}_2^2+1}{{\mathrm{\&sigma;}}_1^2+1}{\mathrm{\&sigma;}}_2^2=({2\mathrm{\&sigma;}}_1^2+{\mathrm{\&sigma;}}_2^2+1)\&CenterDot;\frac{{\mathrm{\&sigma;}}_2^2}{{\mathrm{\&sigma;}}_1^2+1}$

合并得 ${\mathrm{\&sigma;}}_3^2<\frac{{2\mathrm{\&sigma;}}_1^2+{\mathrm{\&sigma;}}_2^2+1}{2}\&CenterDot;\frac{{\mathrm{\&sigma;}}_2^2}{{\mathrm{\&sigma;}}_1^2+1}$

综上所述，得到滤波参数h₁h₂的约束条件

$\left\{\begin{array}{c}{h}_1{h}_2>({\mathrm{\&sigma;}}_1^2+1)({\mathrm{\&sigma;}}_1^2+{\mathrm{\&sigma;}}_3^2)\\ h_1{h}_2>\frac{(2{\mathrm{\&sigma;}}_3^2-{\mathrm{\&sigma;}}_2^2)({\mathrm{\&sigma;}}_1^2+1)({\mathrm{\&sigma;}}_1^2+{\mathrm{\&sigma;}}_2^2+1)}{{\mathrm{\&sigma;}}_2^2}\\ h_1{h}_2<({\mathrm{\&sigma;}}_1^2+{\mathrm{\&sigma;}}_2^2)({\mathrm{\&sigma;}}_1^2+{\mathrm{\&sigma;}}_2^2+1)\end{array}\right.---\left(20\right)$

在 ${\mathrm{\&sigma;}}_3^2<\frac{{2\mathrm{\&sigma;}}_1^2+{\mathrm{\&sigma;}}_2^2+1}{2}\&CenterDot;\frac{{\mathrm{\&sigma;}}_2^2}{{\mathrm{\&sigma;}}_1^2+1}$ 情况下滤波效果最好。

一般而言，有因此且有h₁＞0,h₂＞0，根据式(16)，式(17)可知ω₁-ω₄和ω₂-ω₄分别为从正方向和负方向逼近t的曲线。

本实施例中有则由式(19)可得

${\mathrm{\&sigma;}}_3^2<\frac{{\mathrm{\&sigma;}}_2^2({\mathrm{\&sigma;}}_1^2+{\mathrm{\&sigma;}}_2^2)}{2({\mathrm{\&sigma;}}_1^2+1)}+\frac{{\mathrm{\&sigma;}}_2^2}{2}=\frac{49\&times;53}{2\&times;5}+\frac{49}{2}=284.2$

本实施例中取 ${\mathrm{\&sigma;}}_3^2=196.$ 由于 ${\mathrm{\&omega;}}_1-{\mathrm{\&omega;}}_4=t-\frac{{\mathrm{\&sigma;}}_3^2-{\mathrm{\&sigma;}}_2^2}{h_2},{\mathrm{\&omega;}}_2-{\mathrm{\&omega;}}_4=t+\frac{{\mathrm{\&sigma;}}_2^2}{h_2},$ 其中 $t=\frac{{\mathrm{\&sigma;}}_2^2{h}_1}{({\mathrm{\&sigma;}}_1^2+1)({\mathrm{\&sigma;}}_1^2+{\mathrm{\&sigma;}}_2^2+1)},$ 当t一定时，令 $k_1={\mathrm{\&omega;}}_1-{\mathrm{\&omega;}}_4-t=-\frac{{\mathrm{\&sigma;}}_3^2-{\mathrm{\&sigma;}}_2^2}{h_2},$ 则k₁,k₂的值随h₂变化的曲线如图3所示。

由图可知当h₂＜50时，k₁比较小，使得ω₁-ω₄的值减小；而当h₂＞150，k₁，k₂的变化不明显，需要ω₁-ω₄和ω₂-ω₄都尽量大，则h₂应控制在50-150范围内。

由式(16)得到

$\left\{\begin{array}{c}{h}_1{h}_2>({\mathrm{\&sigma;}}_1^2+1)({\mathrm{\&sigma;}}_1^2+{\mathrm{\&sigma;}}_3^2)=1145\\ h_1{h}_2>\frac{(2{\mathrm{\&sigma;}}_3^2-{\mathrm{\&sigma;}}_2^2)({\mathrm{\&sigma;}}_1^2+1)({\mathrm{\&sigma;}}_1^2+{\mathrm{\&sigma;}}_2^2+1)}{{\mathrm{\&sigma;}}_2^2}=1890\\ h_1{h}_2<({\mathrm{\&sigma;}}_1^2+{\mathrm{\&sigma;}}_2^2)({\mathrm{\&sigma;}}_1^2+{\mathrm{\&sigma;}}_2^2+1)=2862\end{array}\right.$

因此有1890＜h₁h₂＜2862。本实施例中我们取h₁＝20,h₂＝100。

(3)相似窗和干涉区域定位

在待处理图像上以待处理像素为中心划分出5×5大小的相似窗和7×7大小的干涉区域，如图2所示。

(4)相似窗与匹配窗向量提取，具体计算包括以下子步骤：

(4-1)以干涉区域内当前像素为中心确定5×5大小的匹配窗；如图4所示为下一帧图像中干涉区域内的匹配窗。

(4-2)将以像素i为中心的相似窗N_i划分为两个区域，分别是以待处理像素为中心的中央小型窗r×r大小的区域B和剩下部分组成的环形区域A，对以像素j为中心的匹配窗N_j进行同样处理，如步骤(2)中所述，本实施例中弱小目标在相面上不超过3×3像素块大小，因此令r＝3，相似窗和匹配窗分区示意图如图5所示；

(4-3)将相似窗N_i与匹配窗N_j中的区域A(背景模型区域)中的像素灰度值按同样的先后顺序分别组成向量v(N_iA)和v(N_jA)，根据像素灰度值得到向量的示意图如图6所示

(4-4)将相似窗N_i与匹配窗N_j中的区域B(目标模型区域)中的像素灰度值按同样的先后顺序分别组成向量v(N_iB)和v(N_jB)，根据像素灰度值得到向量的示意图如图6所示

(5)相似度测量，具体包括以下子步骤：

(5-1)根据(4-3)中得到的相似窗与匹配窗区域A对应的向量v(N_iA)和v(N_jA)按照如下公式求取两个向量的欧氏距离，得到背景模型相似度DIS_A：

${\mathrm{DIS}}_A={\left|\right|v\left(N_{\mathrm{iA}}\right)-v\left(N_{\mathrm{jA}}\right)\left|\right|}_2^2$

(5-2)根据(4-4)中得到的相似窗与匹配窗区域B对应的向量v(N_iB)和v(N_jB)，并按照如下公式求取两个向量的欧氏距离，得到目标模型相似度DIS_B：

${\mathrm{DIS}}_B={\left|\right|v\left(N_{\mathrm{iB}}\right)-v\left(N_{\mathrm{jB}}\right)\left|\right|}_2^2$

(6)双模权值计算，具体计算包括以下子步骤：

(6-1)利用相似窗N_i与匹配窗N_j对应的背景模型相似度DIS_A和目标模型相似度DIS_B根据如下公式来计算与匹配窗中心像素j对应的权值：

$\mathrm{\&omega;}(i,j)=\exp (\frac{h_1}{{\mathrm{DIS}}_A+1}-\frac{{\mathrm{DIS}}_B}{h_2})$

如步骤(2)中所述，本实施例中令h₁＝20,h₂＝100；

(6-2)将匹配窗N_j对应的权值累加至相似窗N_i对应的所有匹配窗权值累加和Z(i)中，用于权值归一化计算；

(7)背景像素值估计，包括以下子步骤：

(7-1)判断当前处理的匹配窗是否为像素i的干涉区域中最后一个像素对应的匹配窗，如果满足条件，则按(7-2)计算估计值，否则返回(4)继续处理干涉区域的下一像素；

(7-2)按照如下公式计算待处理像素i的背景估计值

$\tilde{f}\left(i\right)=\underset{j\&Element;I}{\mathrm{\&Sigma;}}{\mathrm{\&omega;}}^{\&prime;}(i,j)f\left(j\right)=\underset{j\&Element;I}{\mathrm{\&Sigma;}}\frac{\mathrm{\&omega;}(i,j)}{Z\left(i\right)}f\left(j\right)$

$Z\left(i\right)=\underset{j\&Element;I}{\mathrm{\&Sigma;}}\mathrm{\&omega;}(i,j)$

其中,I为待处理像素i对应的干涉区域，j为干涉区域I中所有匹配窗中心像素，f(j)为其像素灰度值，ω(i,j)为匹配窗N_j对应的权值，Z(i)为相似窗N_i的干涉区域内所有匹配窗权值累加和。

(8)残差获取，包括以下子步骤：

判断当前处理像素是否为图像中最后一个像素，如果满足条件，则将所得估计背景图像与原始图像相减得到目标增强图像，否则后移一个像素转到步骤(3)对其进行背景值计算。

原始图像第一帧如图7所示，第二帧如图8所示，全部像素值处理完后得到原始图像的背景估计图像，如图9所示，背景抑制结果图像如图10所示，背景抑制结果的三维图像如图11所示。

本实施例中累积处理图像1250帧，并将本算法与多种算法比较，通过统计检测概率和虚警率得到各算法性能比较的ROC曲线，分别对输入信噪比的三个范围进行统计，其中满足SNR_in＜2的图像有172帧，统计结果如图12所示，满足2＜SNR_in＜5的有955帧，统计结果如图13所示，满足SNR_in＞5的有123帧，统计结果如图14所示。

本发明对复杂背景下弱小目标检测起着至关重要的影响，在强杂波干扰和目标可利用特征有限的情况下，本发明利用图像的自相关性，通过对像素块及其邻域进行目标模型相似度和背景模型相似度测量，根据相似度测量值计算像素块对应的双模权值，用加权平均值对原始图像进行背景估计，有效恢复出不含目标的背景图像，并且对目标的存在不敏感，通过与原始图像作差得到背景抑制目标增强的图像，提高了弱小目标的检测概率，降低了虚警率，使其适用于目标检测系统的图像滤波处理。

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

Claims (8)

1.一种双模非局部均值滤波方法，其特征在于，包括以下步骤：

(1)输入待处理的图像序列；

(2)根据目标先验知识确定滤波模型参数；

(3)以图像中的待处理像素为中心，划定相似窗，并在下一帧图像中确定干涉区域范围；

(4)根据干涉区域中当前像素确定匹配窗，利用(3)中的相似窗和当前匹配窗的像素灰度值得到相似窗向量和匹配窗向量；

(5)利用步骤(4)中的相似窗向量和匹配窗向量进行背景模型和目标模型的相似度测量；

(6)根据步骤(5)中的背景模型和目标模型相似度测量值计算双模权值；

(7)判断当前处理的匹配窗是否是以干涉区域内最后一个像素为中心的匹配窗，如果满足该条件，则利用步骤(4)-(6)得到的所有匹配窗中心像素值和该匹配窗的权值获得待处理像素对应的估计值，然后进行步骤(8)，否则返回步骤(4)继续处理干涉区域下一个像素；

(8)判断当前处理像素是否为图像中最后一个像素，如果满足该条件，则将所得估计背景图像与原始图像作差，获得目标得到增强而背景得到抑制的图像，否则转到步骤(3)处理图像中的下一像素。

2.如权利要求1所述的双模非局部均值滤波方法，其特征在于，所述步骤(2)具体为，根据输入的目标尺寸的先验知识确定滤波模型的尺寸，滤波模版分别为相似窗，匹配窗和干涉区域，相似窗和匹配窗尺寸相同，都为R×R，干涉区域大小为L×L，且满足R≤L＜2R，窗口模版分为目标模型和背景模型区域，目标模型为以待处理像素为中心的中央小型窗r×r大小的区域B，背景模型区域为剩下部分组成的环形区域A，其中目标尺寸不 超过t×t像素块大小，r,t,R的关系满足R＞r＝t；

采用高斯模型ω＝exp(-DIS_B/h₂)计算目标模型权值，采用复合指数模型ω＝exp(h₁/(DIS_A+1))计算背景模型权值，其中DIS_A和DIS_B分别为背景模型和目标模型的相似度，h₁,h₂为滤波参数，50＜h₂＜150，且有

其中，为噪声标准差，为背景标准差，为目标标准差。

3.如权利要求2所述的双模非局部均值滤波方法，其特征在于，满足关系R＝r+2,L＝R+2。

4.如权利要求1或2所述的双模非局部均值滤波方法，其特征在于，所述步骤(3)中，相似窗和干涉区域的位置和大小按照根据如下方法确定：

相似窗以待处理像素为中心，选取R×R大小的图像块，干涉区域以下一帧图像中与待处理像素同一位置的像素为中心选取L×L大小的图像块。

5.如权利要求1或2所述的双模非局部均值滤波方法，其特征在于，所述步骤(4)包括以下子步骤：

(4-1)以当前像素为中心选取R×R大小的图像块为匹配窗；

(4-2)将以像素i为中心的相似窗N_i划分为两个区域，分别是以待处理像素为中心的中央小型窗r×r大小的区域B(目标模型区域)和剩下部分组成的环形区域A(背景模型区域)，对以像素j为中心的匹配窗N_j进行同样处理；

(4-3)将相似窗N_i与匹配窗N_j中的区域A中的像素灰度值按同样的先后顺序分别组成向量v(N_iA)和v(N_jA)；

(4-4)将相似窗N_i与匹配窗N_j中的区域B中的像素灰度值按同样的 先后顺序分别组成向量v(N_iB)和v(N_jB)。

6.如权利要求5所述的双模非局部均值滤波方法，其特征在于，所述步骤(5)包括以下子步骤：

(5-1)根据(4-3)中得到的相似窗与匹配窗区域A对应的向量v(N_iA)和v(N_jA)按照如下公式求取两个向量的欧氏距离DIS_A：

(5-2)根据(4-4)中得到的相似窗与匹配窗区域B对应的向量v(N_iB)和v(N_jB)，并按照如下公式求取两个向量的欧氏距离DIS_B：

。

7.如权利要求1或2所述的双模非局部均值滤波方法，其特征在于，所述步骤(6)包括以下子步骤：

(6-1)利用相似窗N_i与匹配窗N_j对应的背景模型相似度DIS_A和目标模型相似度DIS_B根据如下公式来计算与匹配窗中心像素j对应的权值：

其中,DIS_A和DIS_B分别为背景模型和目标模型的相似度，h₁,h₂为滤波参数；

(6-2)将匹配窗N_j对应的权值累加至相似窗N_i对应的所有匹配窗权值累加和中，用于权值归一化计算。

8.如权利要求1或2所述的双模非局部均值滤波方法，其特征在于，所述步骤(7)包括以下子步骤：

(7-1)判断当前处理的匹配窗是否为像素i的干涉区域中最后一个像素对应的匹配窗，如果满足该条件，则按(7-2)所述计算背景估计值，否则继续处理干涉区域中下一像素；

(7-2)按照如下公式计算待处理像素i的背景估计值

其中,I为待处理像素i对应的干涉区域，j为干涉区域I中所有匹配窗中心像素，f(j)为其像素灰度值，ω(i,j)为匹配窗N_j对应的权值，Z(i)为相似窗N_i的干涉区域内所有匹配窗权值累加和。