CN104463805B - 一种基于同质性显著度和方向选择的极化sar相干斑噪声抑制方法 - Google Patents

Abstract

本发明具体提供了一种基于同质性显著度和方向选择的极化SAR相干斑噪声抑制方法，实现步骤为：1)读入相干矩阵格式的极化SAR数据得到图像O，计算O的SPAN系数得到SPAN图像S；2)计算同质性显著度图，得到显著图Z；3)判断O中每个点的方向；4)确定每一点的相似邻域D；5)在D内对极化SAR图像O降斑处理；6)将搜索窗尺寸扩大到17×17，对上一步得到的滤波结果重复3)～5)的操作再进行一次滤波，得到最终滤波结果；7)用Pauli向量法将滤波后的相干矩阵合成伪彩图。本发明相比现有技术显著提高了极化SAR图像相干斑的抑制能力，同时很好的保护了边缘与纹理细节信息，可用于极化SAR图像的预处理过程。

Description

一种基于同质性显著度和方向选择的极化SAR相干斑噪声抑 制方法

技术领域

本发明属于图像处理技术领域，涉及极化SAR图像斑点噪声抑制，可用于对极化SAR图像进行滤波，具体是一种基于同质性显著度和方向选择的极化SAR相干斑噪声抑制方法。

背景技术

极化合成孔径雷达(Polarization SAR,POL-SAR)是用来测量地表回波信号的新型成像雷达,它可以在发射回波信号极化脉冲时近乎同时接收回波信号极化脉冲,各极化脉冲之间是相干的。与单极化数据相比,全极化SAR在对地观测和遥感应用方面有着显著的优势。POL-SAR数据记录了地物任意极化状态下的散射回波,既有振幅信息,也有相位信息,极化SAR包含更多的地物信息将大大提高对地物的识别能力。可是极化SAR图像同样受到相干斑噪声影响,增加了图像解译的难度。因此，极化SAR图像降斑对后续图像处理应用的意义重大。

1)精致极化Lee滤波是一个经典的滤波方法，它通过使用方向窗口进行滤波，滤波后的数据在边缘的特性保持方面效果显著，但是，在结构细节信息的保持上，滤波效果并非十分理想，因此在相干斑的抑制中，原始数据的一些特性无法很好地保留。

2)Lee提出的改进的Sigma滤波，它解决了原始Sigma滤波的暗像素点不被滤波和存在滤波数据误差等缺点，并有效保持了亮目标像素点，在边缘保持和细节信息保留上都要优于精致极化Lee滤波方法，但是在边缘和纹理的处理上，由于斑点噪声的影响，这种方法还是不能最好的区分相干斑噪声与边缘纹理信息，滤波后会出现一些“毛刺”，且同质区域的平滑程度也未达到理想效果。

3)非局部方法的极化SAR滤波。Buades等人提出的非局部均值滤波器使用两个图像块的相似程度来作为滤波权值，利用这种思想，Deledalle C.等人提出了PPB滤波器，并应用于了SAR与极化SAR图像的滤波；J.Chen等人提出了Pretest的极化SAR图像滤波方法并都取得了很好的滤波效果。相比于以上两种基于统计特性的滤波方法，非局部的方法在同质区域的平滑程度上取得了很大提高。但与此同时，PPB的方法没有充分利用极化SAR数据的散射信息，Pretest的方法在判别相似点的时候忽略了一些细节，所以仍有需要提高的地方。

发明内容

本发明的目的在于克服上述已有技术的缺点，提出一种结合局部同质显著度与图像块相似性的极化SAR图像降斑方法，以实现在有效抑制相干斑的同时保持点目标和纹理结构的细节信息，提高极化SAR图像的相干斑抑制效果。

实现本发明的技术方案是：一种基于同质性显著度和方向选择的极化SAR相干斑噪声抑制方法，包括如下步骤：

(1)读入相干矩阵格式的极化SAR数据得到图像O，计算O的SPAN系数图像S，计算方法如下:

SPAN＝|S_hh|²+2|S_hv|²+|S_vv|²，

其中S_hh，S_hv和S_vv分别代表极化SAR数据三个通道HH、HV、VV的强度；

(2)利用Wishart距离计算O中以每一像素点为中心的局部匀质度M₁，利用局部方差系数计算S上每一像素点的系数M₂，联合M₁与M₂得到极化SAR图像对应的同质性显著图Z；

(3)判断O中每个像素点的方向：

3a)定义8个3×3的方向窗d1～d8：

3b)对于O上某一像素点x，取SPAN系数图像S上对应的以x为中心的3×3图像块并与8个方向窗d1～d8分别进行卷积，得到8个卷积结果矩阵；

3c)在8个卷积结果矩阵中，选择矩阵中心点值最大的结果对应的方向窗为该像素点x的方向d；

(4)在S中，以x为中心取一个5×5的搜索窗，以窗内每一像素点为中心取图像块与d进行卷积，选择结果大于0的像素点组成x的相似邻域D；

(5)对极化SAR图像O进行双边滤波：

5a)在相似邻域D中对x进行滤波，用x对应的步骤(2)中得到的显著图Z上的值来自适应的调节双边滤波的空域参数h_s与极化域参数h_d，即对像素点x进行了降斑处理；

5b)遍历所有像素点，得到滤波结果；

(6)将搜索窗尺寸扩大到17×17，对上一步得到的滤波结果重复步骤(3)～(5)的操作再进行一次滤波，得到最终滤波结果；

(7)用Pauli向量法将滤波后的相干矩阵合成伪彩图进行输出，得到滤波结果。

上述步骤(2)中的极化SAR图像对应的同质性显著图Z，具体按如下步骤进行计算：

2a)以像素点x为中心取一个3×3大小的图像块，计算x为中心的局部匀质度，即图像块中每一像素点的相干矩阵与x的相似度之和M₁：

其中，T_x表示像素点x的相干矩阵，T_y表示图像块中第y个像素点的相干矩阵，k＝9；同时，在SPAN系数图像S上同样取x的3×3大小的图像块，计算该图像块的方差M₂，即x的局部方差系数：

其中，μ_x表示图像块的均值，P_y表示图像块中第y个点的值，k＝9；

2b)联合M₁与M₂，得到像素点x的局域同质性显著度:

遍历每一像素点后，得到一幅与原始极化SAR同样尺寸大小的同质性显著图Z。

上述步骤5a)所述的在相似邻域D中对x进行双边滤波，具体计算公式如下：

其中，表示滤波后像素点x的相干矩阵，表示x的相似邻域D中第j个点的相干矩阵，分别表示空间域和极化散射域的权重，它们的表达式如下：

这里(m,n)和(p,q)分别表示点x_j和点x的空间坐标；是对相干矩阵中每个元素计算对数操作，其结果为与同大小的矩阵；h_s和h_d分别是空间域和极化域的滤波参数，控制着和的衰减程度；计算h_s和h_d的表达式为：

Z_p,q代表点x在显著图Z上对应点的显著度值；K₁是一个手动设定的参数，这里取K₁＝2。

上述步骤(7)所述的用Pauli向量法将滤波后的相干矩阵合成伪彩图，按如下步骤进行：

7a)将滤波后的相干矩阵的第二行第二列元素T22＝|S_hh-S_vv|²作为待合成伪彩图的红色分量R；

7b)将滤波后的相干矩阵的第三行第三列元素T33＝|S_hv|²作为待合成伪彩图的绿色分量G；

7c)将滤波后的相干矩阵的第一行第一列元素T11＝|S_hh+S_vv|²作为待合成伪彩图的蓝色分量B；

7d)用R、G、B三个分量合成伪彩图。

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

1.本发明利用了自适应搜索窗，在区域边界处能找到更多的相似点，使得滤波更加充分并且不会破坏纹理结构，克服了非局部滤波中依赖于搜索窗盲目选择相似点的缺点；

2.本发明利用了每一点的局部同质显著度，在滤波的过程中自适应的采用了不同的滤波参数，使得在同质程度高的区域尽量平滑，提高了降斑方法的平滑能力；在同质程度低的区域(纹理、结构区域)尽可能的在滤波同时保持细节信息，提高了降斑方法的细节保护能力；

3.本发明采用了双边滤波的框架，联合了空间域与极化域的信息，充分利用了极化SAR数据特征，更为合理。

附图说明

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

图2是本发明仿真采用的四组极化SAR图像,其中区域A、区域B、区域C和区域D是计算等效视数(ENL)时的所选取的区域，区域1、区域2、区域3和区域4是计算基于比值的边缘保持指数(EPD-ROA)时所选取的区域；

图3是本发明对四组图像的同质性显著度结果图像；

图4是本发明与现有的精制极化Lee滤波、改进的Sigma滤波和Pretest滤波对第一组极化SAR图像的滤波结果对比图；

图5是本发明与现有的精制极化Lee滤波、改进的Sigma滤波和Pretest滤波对第二组极化SAR图像的滤波结果对比图；

图6是本发明与现有的精制极化Lee滤波、改进的Sigma滤波和Pretest滤波对第三组极化SAR图像的滤波结果对比图；

图7是本发明与现有的精制极化Lee滤波、改进的Sigma滤波和Pretest滤波对第四组极化SAR图像的滤波结果对比图。

具体实施方式

参照图1，本发明的具体实施步骤包括：

步骤1，读入一幅相干矩阵格式的极化SAR图像O，图像O中的每一像素点都代表着一个3×3的相干矩阵T：

其中A＝S_hh+S_vv，B＝S_hh-S_vv，C＝2S_hv，S_hh代表h方向发射和h方向接受的回波数据，S_vv代表v方向发射和v方向接受的回波数据，S_hv代h方向发射和v方向接受的回波数据，^*表示共轭。同时计算图像O的SPAN系数：

SPAN＝|S_hh|²+2|S_hv|²+|S_vv|²，

得到SPAN系数图像S，其中||表示求模运算。

步骤2，计算极化SAR图像对应的同质性显著图Z：

2.1)在O上以像素点x为中心取一个3×3大小的图像块，计算块中每一像素点的相干矩阵与x的相似度之和M₁：其中T_x表示像素点x的相干矩阵，T_y表示图像块中第y个像素点的相干矩阵，k＝9。同时，在SPAN图像S上同样取x的3×3的邻域，计算该图像块的方差M₂：其中μ_x表示图像块的均值，P_y表示图像块中第y个像素点的值，k＝9。

2.2)联合M₁与M₂，可以得到像素点x的局域同质性显著度:

遍历每一像素点后，可以得到一幅与原始极化SAR同样尺寸大小的同质性显著图Z。

步骤3，判断像素点x的方向：

3.1)定义八个3×3的方向窗d₁～d₂：

3.2)取SPAN系数图像S中对应的以x为中心的3×3图像块并与八个方向窗d₁～d₈分别进行卷积，得到8个卷积结果矩阵；

3.3)在8个卷积结果矩阵中，选择矩阵中心点值最大的结果对应的方向窗为像素点x的方向d。

步骤4，确定像素点x的相似邻域：

4.1)以像素点x为中心，取一个5×5大小的邻域搜索窗；

4.2)在S中以邻域窗内每个像素点为中心取3×3的图像块与3c)中得到的当前像素点x的方向d进行卷积，选出结果矩阵中心点值大于0的像素点组成x的相似邻域D；

步骤5，对极化SAR图像O进行滤波：

5.1)在x的相似邻域D内，对x进行双边滤波估计，计算表达式如下：

其中，表示滤波后像素点x的相干矩阵，表示x的相似邻域D中第j个点的相干矩阵。分别表示空间域和极化散射域的权重，它们的表达式如下：

这里(m,n)和(p,q)分别表示像素点x_j和像素点x的空间坐标。h_s和h_d分别是控制空间域和极化散射域权重衰减程度的参数。这里计算h_s和h_d的计算表达式为：

Z_p,q代表像素点x在步骤2中得到的显著图Z上对应像素点的显著度值。K₁是一个手动设定的参数，这里取K₁＝2。

5.2)计算完每个点后，得到第一次滤波结果图像数据O_p。

步骤6，扩展搜索窗尺寸为17×17，重复步骤3～5的处理，对上一步的结果O_p再进行一次滤波，得到最终滤波结果O_f。

步骤7，利用Pauli向量法将最终滤波结果O_f的相干矩阵合成伪彩图，观察降斑效果。合成方法如下：

7.1)将滤波后的相干矩阵的第二行第二列元素T22＝|S_hh-S_vv|²作为待合成伪彩图的红色分量R；

7.2)将滤波后的相干矩阵的第三行第三列元素T33＝|S_hv|²作为待合成伪彩图的绿色分量G；

7.3)将滤波后的相干矩阵的第一行第一列元素T11＝|S_hh+S_vv|²作为待合成伪彩图的蓝色分量B；

7.4)用R、G、B三个分量合成伪彩图。

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

1.实验条件:

图2是本发明仿真采用的四组极化SAR图像,其中区域A、区域B、区域C和区域D是计算等效视数(ENL)时的所选取的区域，区域1、区域2、区域3和区域4是计算基于比值的边缘保持指数(EPD-ROA)时所选取的区域；本发明的仿真实验使用的是四组真实极化SAR数据。第一组和第二组数据来源于机载CONVAIR SAR的Ottawa地区和Ice地区，如图2(a)、(b)所示；第三组和第四组数据来源于机载AIRSAR的San Francisco地区和Flevoland地区，如图2(c)、(d)所示。

实验的仿真环境为：Matlab R2010a，Intel Core(TM)2Duo2.4GHz,2G RAM，Windows XP Professinoal。

本发明仿真实验对比的方法为精致极化Lee滤波，改进的Sigma滤波，Pretest滤波三种经典方法。

实验的降斑结果衡量使用的是两种客观指标：等效视数(ENL)和基于比值的边缘保持指数(EPD-ROA)。

2.仿真实验内容及结果分析：

实验1，对附图2中四组极化SAR数据用精致极化Lee滤波进行降斑处理。结果分别如图4(a)、图5(a)、图6(a)和图7(a)所示。

实验2，对附图2中四组极化SAR数据用改进的Sigma滤波进行降斑处理。结果分别如图4(b)、图5(b)、图6(b)和图7(b)所示。

实验3，对附图2中四组极化SAR数据用Pretest滤波进行降斑处理。结果分别如图4(c)、图5(c)、图6(c)和图7(c)所示。

实验4，对附图2中四组极化SAR数据用本发明的方法进行降斑处理。结果分别如图4(d)、图5(d)、图6(d)和图7(d)所示。

从图4可以看出，本发明的降斑结果要明显优于精致极化Lee与改进的Sigma的结果，在同质区域达到了较为理想的平滑效果。相比于Pretest的滤波结果，本发明的结果在道路细节的保护、点目标的保持都要更优。表1中区域A的ENL指标也表明了平滑程度上本发明也稍优于Pretest的结果，是四种方法中最好的。表2中区域1的EPD-ROA指数说明了本发明的结果在边缘保持的指标上也是最优的。

表1 各方法结果的等效视数(ENL)评价指标

等效视数	区域A	区域B	区域C	区域D
					Lee	84.1943	69.1889	9.9129	79.1545
Sigma	131.5731	107.2854	12.2525	36.4999
					Pretest	410.6762	209.1253	16.5530	70.2079
本发明	482.9438	239.3284	16.6074	211.6169

表2 各方法结果的EPD-ROA评价指标

从图5的降斑结果中可以看到，精致Lee和Pretest的结果都丢失了很多细节，同时造成了一些纹理区域的模糊化。本发明的结果在尽可能保留细节的同时抑制了斑点噪声，而改进的Sigma滤波结果虽然保留了很多细节，但未能很好地去除斑点噪声。表1中区域B和表2中区域2的评价指标也从客观上说明了本发明的结果是最优的。

图6，从视觉上可以看出，本发明的结果拥有更清晰的效果且点目标亮而明显。Pretest、精致Lee的结果丢失了很多上方海洋区域的细节，且下方的道路线条明显断续不清晰。改进的Sigma的结果仍然噪声保留太多。表1区域C和表2区域3的评价指标可以看出，本发明的结果对于其他三种方法依然具有优势。

表1中的区域D与表2中区域4的指标对应着图7中四组结果的降斑效果。从同质区域D的结果我们可以看出，本发明的结果平滑而无杂点，要远优于其他三种方法。在纹理区域4，本发明的结果线条清晰，很好的保护了细节。

综上，本发明基于同质性显著度和方向选择来进行降斑处理，充分利用极化SAR图像的地表纹理特征，自适应的调节滤波参数，在同质区域的相干斑抑制效果和结构纹理、边缘保持上的效果都很理想，相比于传统经典方法优势明显，因此对于极化SAR数据相干斑噪声抑制的效果显著。

本实施方式中没有详细叙述的部分属本行业的公知的常用手段，这里不一一叙述。以上例举仅仅是对本发明的举例说明，并不构成对本发明的保护范围的限制，凡是与本发明相同或相似的设计均属于本发明的保护范围之内。

Claims (3)

1.一种基于同质性显著度和方向选择的极化SAR相干斑噪声抑制方法，其特征在于，包括如下步骤：

(1)读入相干矩阵格式的极化SAR数据得到图像O，计算O的SPAN系数图像S，计算方法如下:

SPAN＝|S_hh|²+2|S_hv|²+|S_vv|²，

其中S_hh，S_hv和S_vv分别代表极化SAR数据三个通道HH、HV、VV的强度；

(2)利用Wishart距离计算O中以每一像素点为中心的局部匀质度M₁，利用局部方差系数计算S上每一像素点的系数M₂，联合M₁与M₂得到极化SAR图像对应的同质性显著图Z；

该步骤(2)中的极化SAR图像对应的同质性显著图Z，具体按如下步骤进行计算：

2a)以像素点x为中心取一个3×3大小的图像块，计算x为中心的局部匀质度，即图像块中每一像素点的相干矩阵与x的相似度之和M₁：

$M_1=\underset{y=1}{\overset{k}{\Σ}}\left(Tr\right(T_x^{-1}{T}_y)+Tr(T_x{T}_y^{-1})-6),$

其中，T_x表示像素点x的相干矩阵，T_y表示图像块中第y个像素点的相干矩阵，k＝9；同时，在SPAN系数图像S上同样取x的3×3大小的图像块，计算该图像块的方差M₂，即x的局部方差系数：

$M_2=\frac{1}{k}\underset{y=1}{\overset{k}{\Σ}}{(P_y-{\mathrm{\μ}}_x)}^2,$

其中，μ_x表示图像块的均值，P_y表示图像块中第y个点的值，k＝9；

2b)联合M₁与M₂，得到像素点x的局域同质性显著度:

$D_x=\frac{M_1+M_2}{{\mathrm{\μ}}_x},$

遍历每一像素点后，得到一幅与原始极化SAR同样尺寸大小的同质性显著图Z；

(3)判断O中每个像素点的方向：

3a)定义8个3×3的方向窗d1～d8：

$\begin{array}{cccc}{d}_1=\left[\begin{array}{ccc}-1& 0& 1\\ -1& 0& 1\\ -1& 0& 1\end{array}\right]& d_2=\left[\begin{array}{ccc}0& 1& 1\\ -1& 0& 1\\ -1& -1& 0\end{array}\right]& d_3=\left[\begin{array}{ccc}1& 1& 1\\ 0& 0& 0\\ -1& -1& -1\end{array}\right]& d_4=\left[\begin{array}{ccc}1& 1& 0\\ 1& 0& -1\\ 0& -1& -1\end{array}\right]\end{array}$

$\begin{array}{cccc}{d}_5=\left[\begin{array}{ccc}1& 0& -1\\ 1& 0& -1\\ 1& 0& -1\end{array}\right]& d_6=\left[\begin{array}{ccc}0& -1& -1\\ 1& 0& -1\\ 1& 1& 0\end{array}\right]& d_7=\left[\begin{array}{ccc}-1& -1& -1\\ 0& 0& 0\\ 1& 1& 1\end{array}\right]& d_8=\left[\begin{array}{ccc}-1& -1& 0\\ -1& 0& 1\\ 0& 1& 1\end{array}\right]\end{array};$

3b)对于O上某一像素点x，取SPAN系数图像S上对应的以x为中心的3×3图像块并与8个方向窗d1～d8分别进行卷积，得到8个卷积结果矩阵；

3c)在8个卷积结果矩阵中，选择矩阵中心点值最大的结果对应的方向窗为该像素点x的方向d；

(4)在S中，以x为中心取一个5×5的搜索窗，以窗内每一像素点为中心取图像块与d进行卷积，选择结果大于0的像素点组成x的相似邻域D；

(5)对极化SAR图像O进行双边滤波：

5a)在相似邻域D中对x进行滤波，用x对应的步骤(2)中得到的显著图Z上的值来自适应的调节双边滤波的空域参数h_s与极化域参数h_d，即对像素点x进行了降斑处理；

5b)遍历所有像素点，得到滤波结果；

(6)将搜索窗尺寸扩大到17×17，对上一步得到的滤波结果重复步骤(3)～(5)的操作再进行一次滤波，得到最终滤波结果；

(7)用Pauli向量法将滤波后的相干矩阵合成伪彩图进行输出，得到滤波结果。

2.根据权利要求1所述的一种基于同质性显著度和方向选择的极化SAR相干斑噪声抑制方法，其特征在于，其中步骤5a)所述的在相似邻域D中对x进行双边滤波，具体计算公式如下：

${\hat{T}}_x=\frac{1}{\underset{x_j\∈D}{\Σ}(w_j^s\·w_j^d)}\·\underset{x_j\∈D}{\Σ}(w_j^s\·w_j^d)T_{x_j},$

其中，表示滤波后像素点x的相干矩阵，T_xj表示x的相似邻域D中第j个点的相干矩阵，分别表示空间域和极化散射域的权重，它们的表达式如下：

$\begin{array}{c}{w}_j^s=\exp (-\frac{1}{h_s}\·({\left(m-p\right)}^2+{\left(n-q\right)}^2\left)\right)\\ w_j^d=\exp (-\frac{1}{h_d}\·|\left|\log \right(T_{x_j})-\log (T_x\left)\right|{|}_F^2)\end{array},$

这里(m,n)和(p,q)分别表示点x_j和点x的空间坐标；h_s和h_d分别是空间域和极化域的滤波参数，控制着和的衰减程度；计算h_s和h_d的表达式为：

$\begin{array}{c}{h}_s={\left(\frac{1}{K_1}\exp \right({\left(1-Z_{p,q}\right)}^2\left)\right)}^2\\ h_d={\left(\exp \right({\left(1-Z_{p,q}\right)}^2-1\left)\right)}^2\end{array},$

Z_p,q代表点x在显著图Z上对应点的显著度值；K₁是一个手动设定的参数，这里取K₁＝2。

3.根据权利要求1所述的一种基于同质性显著度和方向选择的极化SAR相干斑噪声抑制方法，其特征在于，其中步骤(7)所述的用Pauli向量法将滤波后的相干矩阵合成伪彩图，按如下步骤进行：

7a)将滤波后的相干矩阵的第二行第二列元素T22＝|S_hh-S_vv|²作为待合成伪彩图的红色分量R；

7b)将滤波后的相干矩阵的第三行第三列元素T33＝|S_hv|²作为待合成伪彩图的绿色分量G；

7c)将滤波后的相干矩阵的第一行第一列元素T11＝|S_hh+S_vv|²作为待合成伪彩图的蓝色分量B；

7d)用R、G、B三个分量合成伪彩图。