CN102073989A - 融入pca的非局部均值的极化sar数据相干斑抑制方法 - Google Patents

Abstract

本发明公开了一种基于融入PCA的非局部均值的极化SAR数据相干斑噪声抑制的方法，主要解决现有极化SAR滤波方法不能很好的滤除同质区域的相干斑噪声和不能有效保持边缘细节信息的问题。其实现过程为：(1)输入极化SAR数据的相干矩阵T；(2)对相干矩阵T进行亮目标保持；(3)使用PCA的方法求得span数据的特征向量；(4)对相干矩阵T的元素进行非局部均值滤波，滤波权值由span数据的特征向量计算所得；(5)通过Pauli向量方法将滤波后的相干矩阵T生成伪彩图。本发明与现有技术相比显著提高了极化SAR数据的相干斑噪声抑制的能力，能够有效的平滑同质区域和保持边缘细节信息，可用于极化SAR数据的预处理过程。

Description

融入PCA的非局部均值的极化SAR数据相干斑抑制方法

技术领域

本发明属于图像数据处理技术领域，具体地说是一种相干斑抑制方法，该方法可用于对极化SAR数据的去噪。

背景技术

随着雷达技术的发展，极化SAR已成为SAR的发展趋势，极化SAR能够得到更丰富的目标信息，有利于提高目标检测，辨别和分类能力等等的特性体现了极化SAR系统的优势，但是它和SAR一样，受着相干斑噪声的严重干扰。因此，相干斑的抑制成为一个经久不衰的研究课题。对于极化SAR数据，抑斑的目的在于能够在抑制相干斑的同时且保持数据的极化特性，边缘细节以及纹理信息。现有对极化SAR数据相干斑抑制的方法很多，其中：

1)极化白化滤波PWF是最早的一个滤波方法，该方法通过对极化SAR数据散射矩阵元素的优化组合来完成对span数据的相干斑抑制，但是该方法的缺点在于它只对极化SAR数据中的span数据进行相干斑抑制，而其余极化SAR数据的各元素并没有进行相干斑的抑制。

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

3)最近新提出的改进的sigma滤波，它解决了原始sigma滤波的暗像素不被滤波和滤波数据存在误差等缺点，并有效的保持了亮目标像素，该方法无论在边缘的保持上还是同质区域的平滑上都优于精致极化Lee滤波方法，但是在边缘和纹理的处理上，由于相干斑噪声的影响，这种滤波还是不能最好的区分相干斑噪声和边缘纹理信息，使有用的边缘纹理信息不能被完整的保留。

发明内容

本发明的目的在于克服上述已有技术的缺点，提出一种基于融入PCA的非局部均值的极化SAR数据的相干斑抑制方法，以在抑制相干斑的同时保持亮目标和边缘纹理细节信息，提高极化SAR数据的相干斑抑制效果。

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

(1)将一组极化SAR数据表示为含有9个元素的3×3相干矩阵T，并使用T矩阵中的元素T11和T22对T矩阵进行亮目标检测和保留；

(2)取T矩阵的元素T11、T22和T33获得span数据，其中span＝T11+T22+T33，使用主成分分析方法求出span数据中最大的d个特征向量w₁～w_d，d取6；

(3)对相干矩阵T各元素的非亮目标像素进行如下基于融入PCA的非局部均值滤波：

3a)对T矩阵的任意元素取一个待滤波像素点x，以x为中心确定7×7的局部邻域Φ

和21×21的搜索窗Ω，对待滤波像素点x的局部邻域Φ进行对数变换，将对数变换后的局部邻域拉伸成49×1的待滤波向量Lx；

3b)在搜索窗内取一像素y，取y的7×7的局部邻域Ψ，对该局部邻域Ψ先进行对

数变换，再将其拉伸成49×1的搜索窗向量Ly；

3c)用span数据的d个最大的特征向量w₁～w_d与待滤波向量Lx进行内积，得到待滤波内积结果fx(m，1)，再用span数据的d个最大的特征向量w₁～w_d与搜索窗向量Ly进行内积得到搜索窗内积结果fy(m，1)，m为1～d的整数；

3d)计算待滤波像素x与搜索窗内像素y的相似度d(x，y)：

$d(x,y)=\underset{m=1}{\overset{d}{\mathrm{\Σ}}}{[\mathrm{fx}(m,1)-\mathrm{fy}(m,1)]}^2$

3e)根据相似度d(x，y)计算搜索窗内像素y与待滤波像素x的滤波权值w(x，y)：

$w(x,y)=e^{(-\frac{d(x,y)}{h^2})}$

其中h为一个滤波参数，这里h取5倍的噪声标准差；

3f)逐个扫描搜索窗内的像素，重复执行步骤3b)-3e)，得到最终待滤波像素x的滤波结果

$\hat{z}\left(x\right)=\frac{1}{C\left(x\right)}\underset{y\∈\mathrm{\Ω}}{\mathrm{\Σ}}w(x,y)z\left(y\right)$

其中C(x)为归一化函数，表示为

z(y)为搜索窗Ω内像素y的值；

3g)对T矩阵元素的逐个像素进行上述3a)-3f)的步骤，完成T矩阵所有元素的滤波；

(4)用Pauli向量法对滤波后的T矩阵合成伪彩图。

本发明具有如下优点：

A)本发明由于将PCA融入到非局部均值滤波的思想中，因而比现有的精制极化Lee滤波和改进的sigma滤波，在同质区域的滤波结果上更为平滑；

B)本发明由于使用了非局部均值滤波的思想，比起精制极化Lee滤波和改进的sigma滤波，能够更好的保持数据的边缘、纹理等细节和特征信息；

C)本发明由于使用span数据提取极化SAR数据的特征来对T矩阵元素进行滤波，因而很好的保持了极化SAR数据的极化相关性，保证滤波后T矩阵合成的伪彩图没有失真。

附图说明

图1是本发明的流程图；

图2是本发明使用的两组原始极化SAR数据的图像；

图3是用本发明与现有精致极化Lee滤波和改进的sigma滤波对第一组原始极化SAR数据的滤波结果图；

图4是用本发明与精致极化Lee滤波和改进的sigma滤波对第二组原始极化SAR数据的滤波结果图。

具体实施方式

参照图1，本发明的具体实现步骤如下：

步骤1，对极化SAR数据相干矩阵T的亮目标进行检测并保留。

1a)将极化SAR数据相干矩阵T表示为：

$\left[T\right]=\left[\begin{array}{ccc}{\mathrm{AA}}^{*}& {\mathrm{AB}}^{*}& {\mathrm{AC}}^{*}\\ {\mathrm{BA}}^{*}& {\mathrm{BB}}^{*}& {\mathrm{BC}}^{*}\\ {\mathrm{CA}}^{*}& {\mathrm{CB}}^{*}& {\mathrm{CC}}^{*}\end{array}\right]$

其中[A B C]＝[S_hh+S_vv S_hh-S_vv 2S_hv]，这里S_hh表示h向发射和h向接收的回波数据，S_vv表示v向发射和v向接收的回波数据，S_hv表示h向发射v向接收的回波数据，AA^*到CC^*分别代表T11到T33九个元素，即AA^*代表T11，AB^*代表T12，AC^*代表T13，BA^*代表T21，BB^*代表T22，BC^*代表T23，CA^*代表T31，CB^*代表T32，CC^*代表T33，

由T11＝AA^*＝|S_hh+S_vv|²，T22＝BB^*＝|S_hh-S_vv|²可知，通常它们有着较强的回波值，而T33为|S_hv|²，它的回波值通常很小，因此本实施例只用T11和T22来检测亮目标；

1b)将T11的所有像素从小到大排列，取出第t个像素，得到该像素值K，

为取整符号，n为T11像素总数；

1c)使用3×3滑窗对T11的逐个像素进行扫描，当滑窗中9个像素中大于K的个数超过Tk时，将此3×3区域视为亮目标区域，Tk取5或6；

1d)对T22做与上述步骤1a)-1c)相同的处理，得到T22的亮目标区域；

1e)将T11和T22得到的亮目标区域的位置一起作为整个极化SAR数据T矩阵的亮目标，并保留这些亮目标不被滤除。

步骤2，使用主成分分析方法求出span＝T11+T22+T33数据中最大的d个特征向量w₁～w_d。

2a)对span数据做对数变换，得到变换后的span数据，对span数据做对数变换是为了把乘性噪声模型转化为加性噪声模型；

2b)求对数变换后的span数据的均值

2c)求对数变换后的span的特征矩阵Cy：

$\mathrm{Cy}=\frac{1}{\left|\mathrm{\Ψ}\right|}\underset{i\∈\mathrm{\Ψ}}{\mathrm{\Σ}}(y\left(i\right)-\stackrel{\‾}{y}){(y\left(i\right)-\stackrel{\‾}{y})}^T$

其中，Ψ是搜索窗的子集，大小为7×7，|Ψ|是区域Ψ的像素个数，y(i)表示区域Ψ做对数变化后的第i个像素，特征矩阵Cy的大小由局部区域大小的平方而定，例如是7×7的局部区域时，此时Ψ区域的大小取7×7，获得的特征矩阵为49×49；

2d)计算特征矩阵Cy的特征向量ω，对特征向量ω按照特征值从小到大的顺序排列，取其ω中最后d个特征向量作为span数据最大的d个特征向量w₁～w_d。

步骤3，对相干矩阵T各元素的非亮目标像素进行融入PCA的非局部均值滤波。

3a)对T矩阵的任意元素取一个待滤波像素点x，以x为中心确定7×7的局部邻域Φ和21×21的搜索窗Ω，将x的局部邻域Φ做对数变换，将变换后的局部邻域拉伸成49×1的待滤波向量Lx；

3b)在搜索窗内取一像素y，取y的7×7局部邻域Ψ，将此局部邻域Ψ做对数变换，然后拉伸成49×1的搜索窗向量Ly；

3c)用span数据的d个最大的特征向量w₁～w_d与待滤波向量Lx进行内积，即特征向量中的元素与待滤波向量中的对应元素相乘求和，得到待滤波内积结果fx(m，1)，再用span数据的d个最大的特征向量w₁～w_d与搜索窗向量Ly进行内积，得到搜索窗内积结果fy(m，1)，m为1～d的整数；

3d)计算待滤波像素x与搜索窗像素y的相似度d(x，y)：

$d(x,y)=\underset{m=1}{\overset{d}{\mathrm{\Σ}}}{[\mathrm{fx}(m,1)-\mathrm{fy}(m,1)]}^2$

3e)根据相似度d(x，y)计算搜索窗内像素y与待滤波像素x的滤波权值w(x，y)：

$w(x,y)=e^{(-\frac{d(x,y)}{h^2})}$

其中h为一个修正系数，h取值为5倍的噪声标准差；

3f)逐个扫描搜索窗内的像素，重复执行步骤3b)-3e)，得到最终待滤波像素x的滤波结果

$\hat{z}\left(x\right)=\frac{1}{C\left(x\right)}\underset{y\∈\mathrm{\Ω}}{\mathrm{\Σ}}w(x,y)z\left(y\right)$

其中C(x)为归一化函数，表示为C(x)＝∑_y∈Ωw(x，y)，z(y)为搜索Ω内局部邻域Ψ的第y个像素；

3g)对T矩阵元素的逐个像素进行上述3a)-3f)的步骤，完成T矩阵所有元素的滤波，最终得到整个滤波后的相干矩阵T。

步骤4，使用Pauli向量法对滤波后的整个相干矩阵T合成伪彩图，以观察滤波的效果。

Pauli向量法主要是使用相干矩阵T中的T11，T22和T33三个元素来合成伪彩图，具体步骤如下：

4a)对滤波后的相干矩阵T的元素T22进行开平方处理：

将|S_hh-S_vv|作为待合成伪彩图的红色分量R；

4b)对滤波后的相干矩阵T的元素T33进行开平方处理：

将|S_hv|作为待合成伪彩图的绿色分量G；

4c)对滤波后的相干矩阵T的元素T11进行开平方处理：

将|S_hh+S_vv|作为待合成伪彩图的蓝色分量B；

4d)用R、G、B三个颜色分量作为三基色，根据三基色原理得到最终伪彩图。

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

1.实验条件和内容

实验仿真环境为：MATLAB 7.0.4，VC++6.0，Intel(R)Pentium(R)1CPU 2.4GHz，Window XP Professional。

实验内容包括：本发明分别使用两组极化SAR数据做测试实验。第一组数据是加拿大Ottawa的区域，视数为四，来源于CONVAIR。第二组极化SAR数据是荷兰Flevoland省的区域，视数为四，来源于AIRSAR。

实验方法分别为现有的精致极化Lee滤波、改进sigma滤波和本发明，其中现有的这两种方法都是极化SAR数据空域滤波中最经典效果最显著的两种方法。

评价结果分别用其细节纹理边缘信息的保持以及同质区域的等效视数ENL、均值mean以及标准差std的大小来客观衡量滤波结果的好坏。

2.实验结果

两组极化SAR数据如图2所示，其中图2(a)所示的第一组极化数据为加拿大Ottawa区域，图2(b)所示的第二组极化SAR数据为荷兰Flevoland省的一个局部区域。

(1)用本发明与现有精致极化Lee滤波和改进的sigma滤波对第一组极化SAR数据的滤波结果如图3，其中图3(a)为精致极化Lee滤波结果，图3(b)是改进的sigma滤波结果，图3(c)为本发明滤波结果。从图3(a)可见，精致极化Lee滤波在边缘的滤波效果上不理想，边缘十分模糊，同质区域平滑效果也不好。从图3(b)可见，改进的sigma滤波在同质区域和边缘处的滤波效果都要明显优于精致极化Lee滤波，可是边缘处仍然不流畅有很多斑点噪声存在，有些边缘仍是断断续续的并不连贯，存在很严重的毛刺。从图3(c)可见，在区域1中，比起图3(a)中的精致极化Lee滤波结果和图3(b)中改进的sigma滤波结果，本发明的滤波结果在边缘处很清晰，没有粘连与毛刺并且图3(c)在其余边缘处也比前两种方法处理的结果更加平滑和清晰。由视觉效果和下表1看出，本发明在同质区域滤波效果明显优于前两种滤波方法，并且在边缘纹理细节信息的保持也非常显著。

(2)用本发明与现有精致极化Lee滤波和改进的sigma滤波对第二组极化SAR数据的滤波结果如图4，其中图4(a)是精致极化Lee滤波结果，图4(b)是改进sigma滤波结果，图4(c)是本发明滤波结果。从图4(c)可见，在区域1中，比起图4(a)中的精致极化Lee滤波结果和图4(b)中改进的sigma滤波结果，本发明的滤波结果在边缘处清晰而平滑，没有毛刺和粘连，并且图4(c)在其余边缘处也比前面两幅更加平滑和清晰。

(3)对图2中A和B两处同质区域的等效视数ENL，等效视数ENL利用公式：ENL＝(mean/std)²进行计算，其中mean和std分别代表区域的均值和标准差，计算的结果如表1所示。

表1：各滤波结果的等效视数ENL值

从表1可见，本发明滤波后同质区域的等效视数比精致极化Lee滤波和改进的sigma滤波的结果要好，并且均值保持比精致极化Lee滤波和改进的sigma滤波的好。从而说明，本发明能有效的滤除了同质区域的相干斑噪声。

综上所述，本发明提出的对极化SAR数据的相干斑噪声的抑制方法，通过使用融入PCA的非局部均值思想，能够很好的保持极化SAR数据的相关性，并且在同质区域的平滑效果上和边缘纹理细节信息的保持效果上都很理想，因此本发明对极化SAR数据相干斑噪声的抑制效果显著。

Claims (4)

1.一种融入PCA的非局部均值的极化SAR数据的相干斑抑制方法，包括如下步骤：

(1)将一组极化SAR数据表示为含有9个元素的3×3相干矩阵T，并使用T矩阵中的元素T11和T22对T矩阵进行亮目标检测和保留；

(2)取T矩阵的元素T11、T22和T33获得span数据，其中span＝T11+T22+T33，使用主成分分析方法PCA求出span数据中最大的d个特征向量w₁～w_d，d取6；

(3)对相干矩阵T各元素的非亮目标像素进行如下基于融入PCA的非局部均值滤波：

3a)对T矩阵的任意元素取一个待滤波像素点x，以x为中心确定7×7的局

部邻域Φ和21×21的搜索窗Ω，对待滤波像素点x的局部邻域Φ进行对数变换，将对数变换后的局部邻域拉伸成49×1的待滤波向量Lx；

3b)在搜索窗内取一像素y，取y的7×7的局部邻域Ψ，对该局部邻域Ψ先进行对数变换，再将其拉伸成49×1的搜索窗向量Ly；

3c)用span数据的d个最大的特征向量w₁～w_d与待滤波向量Lx进行内积，得到待滤波内积结果fx(m，1)，再用span数据的d个最大的特征向量w₁～w_d与搜索窗向量Ly进行内积得到搜索窗内积结果fy(m，1)，m为1～d的整数；

3d)计算待滤波像素x与搜索窗内像素y的相似度d(x，y)：

$d(x,y)\underset{m=1}{\overset{d}{\mathrm{\Σ}}}{[\mathrm{fx}(m,1)-\mathrm{fy}(m,1)]}^2$

3e)根据相似度d(x，y)计算搜索窗内像素y与待滤波像素x的滤波权值w(x，y)：

$w(x,y)=e^{(-\frac{d(x,y)}{h^2})}$

其中h为一个滤波参数，这里h取5倍的噪声标准差；

3f)逐个扫描搜索窗内的像素，重复执行步骤3b)-3e)，得到最终待滤波像素x的滤波结果

$\hat{z}\left(x\right)=\frac{1}{C\left(x\right)}\underset{y\∈\mathrm{\Ω}}{\mathrm{\Σ}}w(x,y)z\left(y\right)$

其中C(x)为归一化函数，表示为z(y)为搜索窗Ω内像素y的值；

3g)对T矩阵元素的逐个像素进行上述3a)-3f)的步骤，完成T矩阵所有元素的滤波；

(4)用Pauli向量法将滤波后的相干矩阵T合成伪彩图。

2.根据权利要求书1所述的相干斑抑制方法，其中步骤(1)所述的用T矩阵中的元素T11和T22对T矩阵进行亮目标检测和保留，按如下步骤进行：

2a)将相干矩阵T表示为：

$\left[T\right]=\left[\begin{array}{ccc}{\mathrm{AA}}^{*}& {\mathrm{AB}}^{*}& {\mathrm{AC}}^{*}\\ {\mathrm{BA}}^{*}& {\mathrm{BB}}^{*}& {\mathrm{BC}}^{*}\\ {\mathrm{CA}}^{*}& {\mathrm{CB}}^{*}& {\mathrm{CC}}^{*}\end{array}\right]$

其中[A B C]＝[S_hh+S_vv S_hh-S_vv 2S_hv]，这里S_hh表示h向发射和h向接收的回波数据，S_vv表示v向发射和v向接收的回波数据，S_hv表示h向发射v向接收的回波数据，AA^*到CC^*分别代表T11到T33九个元素，T11＝AA^*＝|S_hh+S_vv|²，T22＝BB^*＝|S_hh-S_vv|²；

2b)将T11的所有像素从小到大排列，取出第t个像素，得到该像素值K，

n为T11像素总数，

为取整符号；

2c)使用3×3滑窗对T11的逐个像素进行扫描，当滑窗中9个像素中大于K的个数超过Tk时将此3×3区域视为亮目标区域，Tk取5或6；

2d)对T22进行与上述步骤2b)-2c)相同的处理，得到T22的亮目标区域；

2e)将T11和T22得到的亮目标区域的位置一起作为整个极化SAR数据T矩阵的亮目标，并保留这些亮目标不被滤除。

3.根据权利要求书1所述的相干斑抑制方法，其中步骤(2)所述的使用主成分分析方法PCA求出span数据中最大的d个特征向量w₁～w_d，按如下步骤进行：

3a)对span数据做对数变换，得到变换后的span数据；

3b)求对数变换后的span数据的均值

3c)求对数变换后的span的特征矩阵Cy：

$\frac{1}{\left|\mathrm{\Ψ}\right|}\underset{i\∈\mathrm{\Ψ}}{\mathrm{\Σ}}(y\left(i\right)-\stackrel{\‾}{y}){(y\left(i\right)-\stackrel{\‾}{y})}^T$

其中，Ψ是搜索窗的子集，大小为7×7，|Ψ|是区域Ψ的像素个数，y(i)表示区域Ψ做对数变化后的第i个像素，特征矩阵Cy的大小由局部区域大小的平方而定；

3d)计算特征矩阵Cy的特征向量ω，对特征向量ω按照特征值从小到大的顺序排列，取其ω中最后d个特征向量作为span数据最大的d个特征向量w₁～w_d。

4.根据权利要求书1所述的相干斑抑制方法，其中步骤(4)所述的用Pauli向量法对滤波后的相干矩阵T合成伪彩图，按如下步骤进行：

4a)对滤波后的相干矩阵T的元素T22进行开平方处理：

将|S_hh-S_vv|作为待合成伪彩图的红色分量R；

4b)对滤波后的相干矩阵T的元素T33进行开平方处理：

将|S_hv|作为待合成伪彩图的绿色分量G；

4c)对滤波后的相干矩阵T的元素T11进行开平方处理：将|S_hh+S_vv|作为待合成伪彩图的蓝色分量B；

4d)用R、G、B三个颜色分量作为三基色，根据三基色原理得到最终伪彩图。

Images