CN103413279A - 基于ad-nsct算法的sar图像去噪方法 - Google Patents

Abstract

本发明公开了一种基于AD-NSCT算法的SAR图像去噪方法，主要改善现有图像去噪方法边缘模糊和噪声去除不充分的现象。其实现过程是：(1)输入一幅SAR图像；(2)对其按四个尺度、四个方向的子带进行非下采样contourlet变换，生成系数矩阵；(3)计算系数矩阵的噪声标准差；(4)计算图像的掩膜矩阵；(5)计算图像的Affinity矩阵，并结合系数矩阵，计算Alpha矩阵；(6)用Alpha矩阵修正对应的掩膜矩阵，并计算后验比；(7)用修正后的掩膜矩阵，计算先验比；(8)用先验比和后验比，计算图像的缩放因子；(9)对系数矩阵进行缩放处理，对其结果进行逆变换，实现SAR图像的去噪。本发明在保证较高平滑度的同时能很好地保持SAR图像的边缘细节信息，使去噪后的SAR图像在视觉效果上得到了显著的改善。

Description

基于AD-NSCT算法的SAR图像去噪方法

技术领域

本发明属于图像处理技术领域，涉及SAR图像去噪方法，可应用于图像处理和识别的预处理等领域，是图像后续更高层次分析、处理的基础。

背景技术

合成孔径雷达(Synthetic Aperture Rader，SAR)图像中由于存在固有的斑点噪声，故不利于图像中场景的自动分析和SAR图像的理解。因此斑点噪声的去除对SAR图像的后续处理例如边缘检测、图像分割等是非常重要的。SAR和光学成像相比有其独特优点，即不受时间和云雾限制。SAR作为地貌和地物传感器，涉及面十分广泛，其高分辨、大区域和全天候工作的优点，使它被广泛地用于地球遥感、海洋研究、资源勘探、实情预报和军事侦察等领域。SAR图像自身固有相干斑噪声，其产生的原因：雷达波束照射在一个相对于雷达波长而言比较粗糙的表面时，而回波信号是由许多散射体所反射的电磁波的合成，由于表面粗糙的原因，在分辨单元内的散射体与接收机之间的距离是变化的，所以接收到的电磁波在频率上是相关的，而在相位上是不相关的，回波信号表现的时强时弱。SAR图像是对连续回波信号进行相关处理而得到的。因此使得SAR图像出现了很大的点与点之间的灰度级波动，这种波动在视觉上呈现出一种类似颗粒的噪声，因此被称为相关斑噪声。简称斑点噪声，它是SAR图像固有的。SAR图像的斑点噪声严重影响了图像的质量，因而大大增加了SAR图像处理和识别的难度。因此滤除相干斑噪声成为SAR图像处理研究领域的一个重要的研究课题。

2008年，Synthetic aperture radar image despeckling via spatially adaptiveshrinkage in the nonsubsampled contourlet transform domain.Journal of ElectronicImaging，17，pp.013013(1-13)这篇文章，提出一种基于非下采样Contourlet变换（nonsubsampled Contourlet transform，NSCT）域中的空间自适应缩放SAR降斑方法。该方法首先对SAR图像进行下采样Contourlet变换，这种变换采用四个尺度，每个尺度上四个方向；然后，对各尺度上的每一方向子带的进行非下采样contourlet变换，并求解方向子带的初始化掩膜矩阵，再利用求得的掩膜矩阵求出图像中像素的先验比和似然比；最后，用似然比对先验比加权，并将加权结果作为NSCT变换方向子带系数的缩放因子，对系数进行缩放，将缩放后的变换系数进行NSCT逆变换处理，实现SAR图像的降斑。此方法中，求解似然比时，利用的是像素为重要信息或非重要信息的数据分布模型。但这种数据分布模型只是一种合理的假设，并没有充分利用图像的冗余信息，求解似然比的方法也有待改善，去噪结果也没达到最佳的平滑度和边缘保持度。

发明内容

本发明的目的在于根据上述已有技术的不足，提出一种基于AD-NSCT算法的SAR图像去噪方法，以提高去噪后SAR图像的平滑度和边缘保持度，实现SAR图像降斑。

为实现上述目的，本发明包括如下步骤：

（1）输入一幅SAR图像，并设置分解级数为n_levels=[2222]；

（2）对输入的SAR图像，按四个尺度、四个方向的子带进行非下采样contourlet变换，生成系数矩阵c^(l，d)，其中l表示尺度，d表示方向子带；

（3）对系数矩阵c^(l，d)，计算各尺度各方向子带的噪声标准差σ_(l，d)；

（4）根据系数矩阵c^(l，d)和噪声标准差σ_(l，d)，对输入的图像，计算各尺度各方向子带的的掩膜矩阵x^(l，d)；

（5）根据系数矩阵c^(l，d)，计算输入图像在各尺度各方向子带的Affinity矩阵W^(l，d)，并利用Affinity矩阵W^(l，d)，计算其对应尺度对应方向子带上的拉普拉斯矩阵L^(l，d)；

（6）根据拉普拉斯矩阵L^(l，d)和系数矩阵c^(l，d)，计算输入图像在各尺度各方向子带的Alpha矩阵α^(l，d)；

（7）用Alpha矩阵α^(l，d)修正对应尺度对应方向子带的掩膜矩阵x^(l，d)，得到修正后的掩膜矩阵x'^(l，d)；

（8）根据Alpha矩阵α^(l，d)，计算输入图像在各尺度各方向子带的后验比ζ^(l，d)：

${\mathrm{\ζ}}^{(l,d)}=\frac{{\mathrm{\α}}^{(l,d)}}{1-{\mathrm{\α}}^{(l,d)}};$

（9）根据修正后的掩膜矩阵x'^(l，d)，计算输入图像在各尺度各方向子带的先验比η^(l，d)；

（10）根据先验比η^(l，d)和后验比ζ^(l，d)，计算输入的图像在各尺度各方向子带的缩放因子ρ^(l，d)：

ρ^(l，d)＝η^(l，d)ζ^(l，d)/(1+η^(l，d)ζ^(l，d))；

（11）根据缩放因子ρ^(l，d)，对各尺度各方向子带的系数矩阵c^(l，d)进行缩放处理，得到缩放后的系数

${\hat{y}}^{(l,d)}={\mathrm{\ρ}}^{(l,d)}{c}^{(l,d)};$

（12）对缩放后的系数矩阵

进行非下采样Contourlet逆变换，逆变换后的结果即为去噪后的SAR图像数据。

本发明与现有的技术相比具有以下优点：

1.本发明采用的Alpha矩阵求解数据分布特征，比起用数据分布特征的分布函数，在SAR去噪方面更加能准确的描述像素点的分布信息，而且Alpha矩阵对SAR图像特征的描述有独特的优势；同时由于用Alpha矩阵修正其对应的掩膜矩阵，这样就使基于掩膜矩阵求得的先验比得到了优化，可以得到更好的像素空间分布特征；

2.本发明利用Alpha矩阵来求解后验比，更加凸显出像素为重要信息的概率，并用后验比对先验比加权，再对系数进行收缩处理，因此对SAR图像噪声有较强的抑制能力，在保证有较高平滑度的同时能很好地保持SAR图像的边缘细节信息，使得去噪后的SAR图像在视觉效果上得到了显著的改善。

附图说明

图1是本发明的实现流程图；

图2是用本发明与现有两种方法对一幅2视数X波段幅度SAR图像进行去噪的比较图；

图3是用本发明与现有两种方法对一幅6视数X波段幅度SAR图像去噪的比较图；

图4是用本发明与现有两种方法对一幅4视数Ku波段强度SAR图像去噪的比较图。

具体实施方式

参照图1，本发明的具体实施过程如下：

步骤1.输入一幅SAR图像，对其进行非下采样Contourlet变换：

设置分解级数为n_levels=[2222]，采用金字塔滤波器和方向滤波器组的非下采样Contourlet变换，对输入SAR图像进行分解；即对输入的SAR图像，按四个尺度、四个方向的子带进行非下采样contourlet变换，生成系数矩阵c^(l，d)，其中l表示尺度，d表示方向子带。

步骤2.对系数矩阵c^(l，d)，计算各尺度各方向子带的噪声标准差σ_(l，d)：

σ_(l，d)＝Median(|c^(l，d)-Median(c^(l，d))|)/0.6745，

其中，Median()为求取中值函数，c^(l，d)为系数矩阵，l∈{1，2，3，4}，d∈{1，2，3，4}。

步骤3.根据系数矩阵c^(l，d)和噪声标准差σ_(l，d)，对输入的SAR图像，计算各尺度各方向子带的掩膜矩阵x^(l，d)。

本步骤的实现可采用硬阈值法或软阈值法，本实例采用硬阈值法，其按如下步骤进行：

3a）计算l=4时的掩膜矩阵，即第4个尺度上第d个子带的掩膜矩阵x^(4，d)：

$x_k^{(4,d)}=\left\{\begin{array}{c}0,|m_k^{(4,d)}\&CenterDot;c_k^{(1,d)}|<{\mathrm{\&sigma;}}_{(4,d)}^2\\ 1,|m_k^{(4,d)}\&CenterDot;c_k^{(1,d)}|\&GreaterEqual;{\mathrm{\&sigma;}}_{(4,d)}^2\end{array}\right.,$

其中，

为掩膜矩阵x^(4，d)中第k个像素点的元素，d∈{1，2，3，4}，为第4个尺度上第d个方向子带系数的噪声系数方差，

为系数矩阵c^(4，d)中第k个像素点的元素；

3b）计算当l=3时的掩膜矩阵，即第3个尺度上第d个子带的掩膜矩阵x^(3，d)：

$x_k^{(3,d)}=\left\{\begin{array}{c}0,|m_k^{(3,d)}\&CenterDot;{\hat{y}}_k^{(4,d)}|<{\mathrm{\&sigma;}}_{(3,d)}^2\\ 1,|m_k^{(3,d)}\&CenterDot;{\hat{y}}_k^{(4,d)}|\&GreaterEqual;{\mathrm{\&sigma;}}_{(3,d)}^2\end{array}\right.,$

其中，

为掩膜矩阵x^(3，d)中第k个像素点的元素，d∈{1，2，3，4}，

为第3个尺度上第d个方向子带的噪声系数方差，为系数矩阵c^(3，d)中第k个像素点的元素，

是缩放后的系数矩阵

中第k个像素点的元素；

3c）计算当l=2时的掩膜矩阵，即第2个尺度上第d个子带的膜矩阵x^(2，d)：

$x_k^{(2,d)}=\left\{\begin{array}{c}0,|m_k^{(2,d)}\&CenterDot;{\hat{y}}_k^{(3,d)}|<{\mathrm{\&sigma;}}_{(2,d)}^2\\ 1,|m_k^{(2,d)}\&CenterDot;{\hat{y}}_k^{(3,d)}|\&GreaterEqual;{\mathrm{\&sigma;}}_{(2,d)}^2\end{array}\right.,$

其中，为掩膜矩阵x^(2，d)中第k个像素点的元素，d∈{1，2，3，4}，

为第2个尺度上第d个方向子带的噪声系数方差，为系数矩阵c^(2，d)中第k个像素点的元素，

是缩放后的系数矩阵

中第k个像素点的元素；

3d）计算当l=1时的掩膜矩阵，即第1个尺度上第d个子带的掩膜矩阵x^(1，d)：

$x_k^{(1,d)}=\left\{\begin{array}{c}0,|m_k^{(1,d)}\&CenterDot;{\hat{y}}_k^{(2,d)}|<{\mathrm{\&sigma;}}_{(1,d)}^2\\ 1,|m_k^{(1,d)}\&CenterDot;{\hat{y}}_k^{(2,d)}|\&GreaterEqual;{\mathrm{\&sigma;}}_{(1,d)}^2\end{array}\right.,$

其中，

为掩膜矩阵x^(1，d)中第k个像素点的元素，d∈{1，2，3，4}，为第1个尺度上第d个方向子带的噪声系数方差，

为系数矩阵c^(1，d)中第k个像素点的元素，

是缩放后的系数矩阵

中第k个像素点的元素。

步骤4.根根据系数矩阵c^(l，d)，计算输入图像在各尺度各方向子带的Affinity矩阵W^(l，d)，并利用Affinity矩阵W^(l，d)，计算其对应尺度对应方向子带上的拉普拉斯矩阵L^(l，d)。

4a）对第l尺度上第d方向子带的掩膜矩阵x^(l，d)进行腐蚀处理；

4b）计算第l个尺度上第d个子带方向的Affinity矩阵W^(l，d)：

$W^{(l,d)}=\frac{1}{\left|{\mathrm{\&omega;}}_k\right|}\mathrm{\&Sigma;}(1+\frac{{(m_i^{(l,d)}-{\mathrm{\&mu;}}_k)}^T(m_j^{(l,d)}-{\mathrm{\&mu;}}_k)}{{\mathrm{\&sigma;}}_k^2+\mathrm{\&epsiv;}}),$

其中，

为系数矩阵c^(l，d)中第i个像素点的元素， 为系数矩阵c^(l，d)中第j个像素点的元素，μ_k和

为k邻域窗口的均值和方差，|ω_k|为k邻域窗口的像素个数，k为系数矩阵c^(l，d)中第i个像素点的邻域窗口，ε为平滑参数；

4c）根据Affinity矩阵W^(l，d)，计算拉普拉斯矩阵L^(l，d)：

L^(l，d)＝D^(l，d)-W^(l，d)，

其中，D^(l，d是一个对角矩阵，该对角矩阵中的第i行第i列为： $D^{(l,d)}(i,i)={\mathrm{\&Sigma;}}_j{W}^{(l,d)}(i,,W^{(l,d)}(i,j)$ 表示W^(l，d)矩阵中的第i行第j列。

步骤5.根据拉普拉斯矩阵L^(l，d)和系数矩阵c^(l，d)，计算输入图像在各尺度各方向子带的Alpha矩阵α^(l，d)。

5a）建立α^(l，d)矩阵的代价函数J(α^(l，d))：

J(α^(l，d))＝α^(l，d)TL^(l，d)α^(l，d)+λ(α^(l，d)-E^T)^TL_F(α^(l，d)-E^T)+λα^(l，d)TL_Bα^(l，d)，

其中，α^(l，d)T表示α^(l，d)的转置，L_F表示腐蚀后的掩膜矩阵为1的元素即代表前景先验信息，L_B表示腐蚀后的掩膜矩阵为0的元素即代表背景先验信息，E＝[1，1，１，...，1]为N维的向量，N为图像的像素点数，λ为控制平滑的权重值；

5b）计算代价函数J(α^(l，d))的最小值，得到Alpha矩阵α^(l，d)：

α^(l，d)＝[L^(l，d)+λ(L_F+L_B)]^-1λL_FE^T，

其中，E^T表示E的转置。

步骤6.用Alpha矩阵α^(l，d)修正对应尺度对应方向子带的掩膜矩阵x^(l，d)。

6a）设定门限值为0.5；

6b）将Alpha矩阵α^(l，d)中元素值与门限值0.5做比较，若Alpha矩阵α^(l，d)元素值大于等于0.5，则x^(l，d)矩阵中对应位置的元素值就设置为1，否则为0，得到修正后的掩膜矩阵x'^(l，d)。

步骤7.根据Alpha矩阵α^(l，d)，计算输入图像在各尺度各方向子带的后验比ζ^(l，d)：

${\mathrm{\&zeta;}}^{(l,d)}=\frac{{\mathrm{\&alpha;}}^{(l,d)}}{1-{\mathrm{\&alpha;}}^{(l,d)}}.$

步骤8.根据修正后的掩膜矩阵x'^(l，d)，计算输入图像在各尺度各方向子带的先验比η^(l，d)：

${\mathrm{\&eta;}}_k^{(l,d)}=\exp [\mathrm{\&delta;}\&CenterDot;{\mathrm{\&Sigma;}}_{s\&Element;N_k}({2x\&prime;}_s^{(l,d)}-1)],$

其中，δ为平滑系数，其最佳取值为0.5，N_k为第k个像素点的8邻域，s为k的8邻域索引号，

为修正后掩膜矩阵x'^(l，d)第k个像素点的邻域元素。

步骤9.根据先验比η^(l，d)和后验比ζ^(l，d)，计算输入的图像在各尺度各方向子带的缩放因子ρ^(l，d)：

${\mathrm{\&rho;}}^{(l,d)}={\mathrm{\&eta;}}^{(l,d)}{\mathrm{\&zeta;}}^{(l,d)}/(1+{\mathrm{\&eta;}}^{(l,d)}{\mathrm{\&zeta;}}^{(l,d)}).$

步骤10.根据缩放因子ρ^(l，d)，对各尺度各方向子带的系数矩阵c^(l，d)进行缩放处理，得到缩放后的系数

${\hat{y}}^{(l,d)}={\mathrm{\&rho;}}^{(l,d)}{c}^{(l,d)}.$

步骤11.对缩放后的系数矩阵

进行非下采样Contourlet逆变换，逆变换后的结果即为去噪后的SAR图像数据。

本发明的效果通过以下实验进一步说明：

1.实验条件：

仿真实验环境为：MATLAB R2011a，CPU intel Pentium Dual-Core E53002.60GHz，内存2G，Windows XP SP3。

2.实验内容：

用基于非下采样Contourlet变换域中的空间自适应缩放SAR降斑算法和基于线性最小化均方误差小波缩放的非局部SAR图像降斑算法分别对SAR图像进行去噪处理，并用本发明和上述两种去噪方法进行比较。其中基于非下采样Contourlet变换域中的空间自适应缩放SAR降斑算法D-NSCT称为方法1；基于线性最小化均方误差小波缩放的非局部SAR图像降斑算法NLWS称为方法2。

实验1：用本发明和上述两种去噪方法对一幅大小为256×256的2视数X波段幅度SAR1图像进行去噪，结果如图2所示。其中图2(a)所示是原始的SAR1图像进，图2(b)为基于非下采样Contourlet变换域中的空间自适应缩放SAR降斑算法的去噪图像，图2(c)为基于线性最小化均方误差小波缩放的非局部SAR图像降斑算法的去噪图像，图2(d)为本发明去噪后的图像。

实验2：用本发明和上述两种去噪方法对大小为256×256的一幅6视数X波段幅度SAR2图像进行去噪，结果如图3所示。其中图3(a)所示是原始的SAR2图像，图3(b)为基于非下采样Contourlet变换域中的空间自适应缩放SAR降斑算法的去噪图像，图3(c)为基于线性最小化均方误差小波缩放的非局部SAR图像降斑算法的去噪图像，图3(d)为本发明去噪后的图像。

实验3：用本发明和上述两种去噪方法对大小为256×256的一幅4视数Ku波段强度SAR3图像进行去噪，结果如图4所示。其中图4(a)所示是原始的SAR3图像，图4(b)为基于非下采样Contourlet变换域中的空间自适应缩放SAR降斑算法的去噪图像，图4(c)为基于线性最小化均方误差小波缩放的非局部SAR图像降斑算法的去噪图像，图4(d)为本发明去噪后的图像。

3.实验结果与分析：

从图2(b)、图3(b)和图4(b)均可以看出，对于SAR1、SAR2和SAR3图像，基于非下采样Contourlet变换域中的空间自适应缩放SAR降斑算法得到的去噪后的图像，在平滑区域的平滑度和边缘边缘信息保持度方面未能达到最优，还有待提高；

从图2(c)可以看出，基于线性最小化均方误差小波缩放的非局部SAR图像降斑算法得到的去噪后的图像，在某些直线边缘信息有较好的保持度，但整体上在很多重要的边缘信息出现过模糊现象，使视觉观察度很差；

从图2(d)可以看出，本发明得到的去噪后的图像平滑度较高，能有效的抑制SAR图像中噪声，更好地保留了细节信息和边缘信息。

采用客观评价指标equivalent number of looks(ENL)和edge sustain index(ESI)，对本发明与上述两种方法的去噪后的图像性能进行比较，其结果如表1所示。其中，在计算ENL值时，图像中所选区域如图2(a)、图3(a)和图4(a)的方框所示，图2(a)区域1大小50×50，区域2大小40×40，图3(a)和图4(a)所的选区域1和区域2大小均为40×40；在计算ESI值时，所选区域为整幅图像。

表1.本发明与上述两种方法得到的去噪后图像的ENL值和ESI值

从上表中可以看出，与现有的两种去噪方法相比，用本发明去噪后的SAR图像的ENL值和ESI值均有较大幅度的提高。

综上，本发明具有较强的抑制能力，在保证有较高平滑度的同时能很好地保持SAR图像的边缘细节信息，使得去噪后的SAR图像在视觉效果上得到了显著的改善。

Claims (8)

1.一种基于AD-NSCT算法的SAR图像去噪方法，包含以下步骤： 

（1）输入一幅SAR图像，并设置分解级数为n_levels=[2222]； 

（2）对输入的SAR图像，按四个尺度、四个方向的子带进行非下采样contourlet变换，生成系数矩阵c^(l，d)，其中l表示尺度，d表示方向子带； 

（3）对系数矩阵c^(l，d)，计算各尺度各方向子带的噪声标准差σ_(l，d)； 

（4）根据系数矩阵c^(l，d)和噪声标准差σ_(l，d)，对输入的图像，计算各尺度各方向子带的的掩膜矩阵x^(l，d)； 

（5）根据系数矩阵c^(l，d)，计算输入图像在各尺度各方向子带的Affinity矩阵W^(l，d)，并利用Affinity矩阵W^(l，d)，计算其对应尺度对应方向子带上的拉普拉斯矩阵L^(l，d)； 

（6）根据拉普拉斯矩阵L^(l，d)和系数矩阵c^(l，d)，计算输入图像在各尺度各方向子带的Alpha矩阵α^(l，d)； 

（7）用Alpha矩阵α^(l，d)修正对应尺度对应方向子带的掩膜矩阵x^(l，d)，得到修正后的掩膜矩阵x'^(l，d)； 

（8）根据Alpha矩阵α^(l，d)，计算输入图像在各尺度各方向子带的后验比ζ^(l，d)： 

（9）根据修正后的掩膜矩阵x'^(l，d)，计算输入图像在各尺度各方向子带的先 

验比η^(l，d)； 

（10）根据先验比η^(l，d)和后验比ζ^(l，d)，计算输入的图像在各尺度各方向子带的缩放因子ρ^(l，d)： 

ρ^(l，d)＝η^(l，d)ζ^(l，d)/(1+η^(l，d)ζ^(l，d))； 

（11）根据缩放因子ρ^(l，d)，对各尺度各方向子带的系数矩阵c^(l，d)进行缩放处 理，得到缩放后的系数

（12）对缩放后的系数矩阵

进行非下采样Contourlet逆变换，逆变换后的结果即为去噪后的SAR图像数据。 

2.根据权利要求1所述的方法，其中步骤(3)所述的计算各尺度各方向子带的噪声标准差σ_(l，d)，通过如下公式进行： 

σ_(l，d)＝Median(|c^(l，d)-Median(c^(l，d))|)/0.6745， 

其中，Median()为求取中值函数，c^(l，d)为系数矩阵，l∈{1，2，3，4}，d∈{1，2，3，4}。 

3.根据权利要求1所述的方法，其中步骤(4)所述的计算各尺度各方向子带的掩膜矩阵x^(l，d)，有硬阈值法和软阈值法，本步骤采用硬阈值法，通过如下公式进行： 

4a）计算l=4时的掩膜矩阵，即第4个尺度上第d个子带的掩膜矩阵x^(4，d)： 

其中，

为掩膜矩阵x^(4，d)中第k个像素点的元素，d∈{1，2，3，4}，

为第4个尺度上第d个方向子带系数的噪声系数方差，

为系数矩阵c^(4，d)中第k个像素点的元素； 

4b）计算当l=3时的掩膜矩阵，即第3个尺度上第d个子带的掩膜矩阵x^(3，d)： 

其中，

为掩膜矩阵x^(3，d)中第k个像素点的元素，d∈{1，2，3，4}，

为第3个尺度上第d个方向子带的噪声系数方差，

为系数矩阵c^(3，d)中第k个像素点的元素，

是缩放后的系数矩阵

中第k个像素点的元素； 

4c）计算当l=2时的掩膜矩阵，即第2个尺度上第d个子带的膜矩阵x^(2，d)： 

其中，

为掩膜矩阵x^(2，d)中第k个像素点的元素，d∈{1，2，3，4}，

为第2个尺度上第d个方向子带的噪声系数方差，

为系数矩阵c^(2，d)中第k个像素点的元素，

是缩放后的系数矩阵

中第k个像素点的元素； 

d）计算当l=1时的掩膜矩阵，即第1个尺度上第d个子带的掩膜矩阵x^(1，d)： 

其中，

为掩膜矩阵x^(1，d)中第k个像素点的元素，d∈{1，2，3，4}，

为第1个尺度上第d个方向子带的噪声系数方差，

为系数矩阵c^(1，d)中第k个像素点的元素，

是缩放后的系数矩阵

中第k个像素点的元素。 

4.根据权利要求1所述的方法，其中步骤(5)所述的计算输入图像在各尺度各方向子带的Affinity矩阵W^(l，d)，通过如下公式进行： 

其中，为系数矩阵c^(l，d)中第i个像素点的元素，

为系数矩阵c^(l，d)中第j个像素点的元素，μ_k和为k邻域窗口的均值和方差，|ω_k|为k邻域窗口的像素个数，k为系数矩阵c^(l，d)中第i个像素点的邻域窗口，ε为平滑参数。 

5.根据权利要求1所述的方法，其中步骤(5)所述的计算拉普拉斯矩阵L^(l，d)，通过如下公式进行： 

L^(l，d)＝D^(l，d)-W^(l，d)， 

其中，D^(l，d)是一个对角矩阵，且D^(l，d)(i，i)表示D(^l，d)矩 阵中的第i行第i列，W^(l，d)(i，j)表示W^(l，d)矩阵中的第i行第j列。 

6.根据权利要求1所述的方法，其中步骤(6)所述的计算各尺度各方向子带的Alpha矩阵α^(l，d)，通过如下公式进行： 

6a）建立α^(l，d)矩阵的代价函数J(α^(l，d))： 

J(α^(l，d))＝α^(l，d)TL^(l，d)α^(l，d)+λ(α^(l，d)-E^T)^TL_F(α^(l，d)-E^T)+λα^(l，d)TL_Bα^(l，d)， 

其中，α^(l，d)T表示α^(l，d)的转置，L_F表示腐蚀后的掩膜矩阵为1的元素即代表前景先验信息，L_B表示腐蚀后的掩膜矩阵为0的元素即代表背景先验信息，E＝[1，1，1，...，1]为N维的向量，N为图像的像素点数，λ为控制平滑的权重值； 

6b）计算代价函数J(α^(l，d))的最小值，得到Alpha矩阵α^(l，d)： 

α^(l，d)＝[L^(l，d)+λ(L_F+L_B)]^-1λL_FE^T， 

其中，E^T表示E的转置。 

7.根据权利要求1所述的方法，其中步骤(7)所述的用Alpha矩阵α^(l，d)修正对应尺度对应方向子带的掩膜矩阵x^(l，d)，是根据设定门限值为0.5，将Alpha矩阵α^(l，d)中元素值与门限值0.5做比较，如果Alpha矩阵α^(l，d)元素值大于等于0.5，则x^(l，d)矩阵中对应的行列元素的值就设置为1，否则为0，得到修正后的掩膜矩阵x'^(l，d)。 

8.根据权利要求1所述的方法，其中步骤(9)所述的计算各尺度各方向子带的先验比通过如下公式进行： 

其中，δ为平滑系数，其最佳取值为0.5，N_k为第k个像素点的8邻域，s为k的8邻域索引号，

为修正后掩膜矩阵x'^(l，d)第k个像素点的邻域元素。 

Images