CN103871031A - 基于核回归的sar图像相干斑抑制 - Google Patents

Abstract

本发明公开了一种基于核回归的SAR图像相干斑抑制方法，主要是解决现有的核回归去噪方法不能应用于噪声模型为乘性的SAR图像中。其实现过程是：对SAR图像选取邻域窗；用索贝尔sobel算子计算图像块的水平梯度与垂直梯度；依据得到的梯度求得局部自适应协方差矩阵，计算平滑矩阵；通过权影响估计该像素点的理想值得到基本的去噪结果；在基本的滤波结果上重复该过程，得到满意的滤波效果。本发明在抑制相干斑的同时保持点目标和边缘纹理细节信息，提高了SAR图像的相干斑抑制效果，可用于SAR图像目标识别和地物分类。

Description

基于核回归的SAR图像相干斑抑制

技术领域

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

背景技术

合成孔径雷达SAR是一种高分辨率成像雷达。SAR不受气候、昼夜的影响，具备全天时、全天候的特点，在对地观测系统中，装备高分辨率的合成孔径雷达的航空系统机动性好、响应速度快，在灾害应急处理方面优势显著。作为卫星数据的有效补充，雷达在地形测绘、灾害监测评估、环境与国土资源监测等领域具有广泛的应用价值。特别是在灾害应急处理方面，雷达系统能及时为灾害应急处理提供信息丰富、清晰准确的图像，对及时开展救援、进行受灾评估、减轻灾害损失具有重要意义。但由于SAR系统采用的是相干微波源进行成像，回波的振幅在相干叠加时会相互叠加或抵消，得到的总回波强度与子回波平均强度之间会存在一定的偏差，在SAR图像中这种随机偏差以强噪声的形式呈现，这种噪声通常被称为相干斑。因为原理的缺陷性，SAR图像中相干斑影响了后续的图像解译，因此如何抑制SAR图像中的相干斑，提高图像的解读能力以及获得更多的信息成为一个重要的问题。SAR图像相干斑抑制的目标是在去除图像中的相干斑及尽可能的保持图像的细节信息。

SAR图像的相干斑是一种复杂的乘性噪声模型，对于斑点噪声这种特殊的性质，在过去的近二十年中，人们提出了很多经典的SAR图像相干斑抑制方法。如Lee滤波，增强Lee滤波，Kuan滤波等。这些方法是用一个已定义好的滤波器窗来估计局域相干斑的方差，并进行滤波处理，其结果通常过分的平滑了图像的边缘细节信息，在一定程度上取得了较好的效果，但存在模糊了边缘和点目标的缺陷。除了在空域进行图像相干斑抑制外，1995年，美国学者Donoho把小波理论引入到图像去噪中，提出了小波软阈值方法。小波软阈值法开创了变换域进行图像去噪的先河，之后涌现了许多变换域去噪的优秀成果，包括多尺度变换的图像去噪方法。小波软阈值方法是一种非线性的算法，依然存在破坏图像细节信息的问题，对图像的辐射特性保持也不好。

核回归方法本质上是一种空域局部平均的去噪方法，不同的是它依赖于被处理数据本身来指定这个模型的结构，这个模型是基于泰勒局部展开，采用的是非参数化方法，模型建立在图像的数据特征上。核回归分析应用于图像去噪最早是由H.Takeda和S.Farsiu等人于2005年提出，随后又得到了进一步的发展。其滤波效果与目前非常优秀的非局部均值方法NLM和三维块匹配的联合滤波BM3D方法相差无几，但是以上核回归算法目前只应用于噪声模型为加性的自然图像中，不能用到噪声模型为乘性的SAR图像中。

发明内容

本发明的目的在于提出了一种基于核回归的SAR图像相干斑抑制方法，将性能优异的核回归去噪方法扩展到噪声模型为乘性的SAR图像中，实现对SAR图像充分滤波的同时保持点目标和边缘细节信息。

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

(1)读入一幅SAR图像，并进行边界镜像对称扩展，扩展半径为N，取图像中的一个像素点x，以该像素点为中心取(2×N+1)大小的邻域窗w;

(2)用索贝尔（Sobel）算子计算邻域窗内w每一个像素点的水平梯度与垂直梯度；

(3)根据得到的水平梯度与垂直梯度计算以像素点x为中心的邻域窗w自适应协方差矩阵C；

(4)依据步骤(4)得到的局部自适应协方差矩阵C，计算以像素点x为中心的邻域窗w平滑矩阵H；

(5)依据计算得到的平滑矩阵H，利用权重影响，求得该像素点步骤1中的像素点x的相干斑抑制结果；

(6)对SAR图像每一个像素点按步骤(1)—(5)进行处理，得到每一个像素点的相干斑抑制结果，从而得到所述SAR图像的相干斑抑制结果。

(7)将步骤6得到的结果，按照步骤(1)—步骤(6)进行迭代处理，得到最终的所述SAR图像的相干斑抑制结果。

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

1、本发明是基于SAR图像的乘性噪声模型而提出的，更加符合SAR图像的相干斑抑制，对同质区域的滤波更加平滑；

2、本发明充分利用了图像的结构信息，克服了增强Lee的滤波等方法无法表示图像大量丰富的局部结构；

3、本发明用索贝尔Sobel算子计算图像的初始水平梯度与垂直梯度，减少计算量，能够快速有效的处理大幅面的SAR图像；

附图说明

图1是本发明的流程图；

图2是本发明仿真使用的barbara自然图像模拟的SAR图像；

图3是本发明仿真使用的光学图像模拟的SAR图像；

图4是本发明仿真使用的真实的幅度SAR图像；

图5是本发明仿真使用的真实的强度SAR图像；

图6是用现有方法及本发明对barbara自然图像模拟的SAR图像滤波后的结果图；

图7是用现有方法及本发明对光学图像模拟的SAR图像滤波后的结果图；

图8是用现有方法及本发明对真实的幅度SAR图像滤波后的结果图；

图9是用现有方法及本发明对真实的强度SAR图像滤波后的结果图。

具体实施步骤

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

步骤1，对一幅SAR图像，进行边界镜像对称扩展，扩展半径为N，取图像中的一个像素点x，以该像素点为中心取(2×N+1)大小的邻域窗w。

邻域窗太小会出现过平滑现象，太大会导致降斑效果不明显。优选地，N=3。

步骤2，用索贝尔（Sobel）算子计算邻域窗内每一个像素点的水平梯度与垂直梯度；

具体过程如下：

2a)设置水平梯度算子模板与垂直梯度算子模板分别用来检测图像的垂直边缘及水平边缘，模板内的数字为模板系数，水平梯度算子模板的方向与垂直边缘的方向相垂直，垂直梯度算子模板的方向与水平边缘的方向相垂直。

$\frac{1}{4}*\left[\begin{array}{ccc}-1& -2& -1\\ 0& 0& 0\\ 1& 2& 1\end{array}\right],\frac{1}{4}*\left[\begin{array}{ccc}-1& 0& 1\\ -2& 0& 2\\ -1& 0& 1\end{array}\right]$

水平梯度算子模板     垂直梯度算子模板

将水平梯度算子模板与垂直梯度算子模板表达为：

${\left[M\right]}^i=\left[\begin{array}{ccc}{M}_{-1,-1}^i& M_{-\mathrm{1,0}}^i& M_{-\mathrm{1,1}}^i\\ M_{0,-1}^i& M_{\mathrm{0,0}}^i& M_{\mathrm{0,1}}^i\\ M_{1,-1}^i& M_{1.0}^i& M_{\mathrm{1,1}}^i\end{array}\right]$

式中，i=1,2分别代表水平梯度算子模板，垂直梯度算子模板。

代表水平梯度算子模板或垂直梯度算子模板内的中心位置值。当i=1时，

${\left[M\right]}^1=\left[\begin{array}{ccc}{M}_{-1,-1}^1& M_{-\mathrm{1,0}}^1& M_{-\mathrm{1,1}}^1\\ M_{0,-1}^1& M_{\mathrm{0,0}}^1& M_{\mathrm{0,1}}^1\\ M_{1,-1}^1& M_{1.0}^1& M_{\mathrm{1,1}}^1\end{array}\right]=\frac{1}{4}*\left[\begin{array}{ccc}-1& -2& -1\\ 0& 0& 0\\ 1& 2& 1\end{array}\right].$

2b)步骤1中的邻域窗w内每一个像素点的垂直梯度与水平梯度的计算公式为：

$z_i\left(x_{j,k}\right)=\underset{m=-1}{\overset{1}{\mathrm{\Σ}}}\underset{n=-1}{\overset{1}{\mathrm{\Σ}}}{x}_{j,k}{M}_{m,n}^i$

i=1,2，z₁(x_j,k),z₂(x_j,k)分别是邻域窗w内当前处理的像素点x_j,k的垂直梯度和水平梯度，其中j,k的取值范围均为[1,2×N+1]。

步骤3，根据步骤2得到的水平梯度与垂直梯度计算以像素点x为中心的邻域窗w的自适应协方差矩阵C：

$C\≈\left[\begin{array}{cc}\underset{x_{j,k}\∈w}{\mathrm{\Σ}}{z}_1\left(x_{j,k}\right)z_1\left(x_{j,k}\right)& \underset{x_{j,k}\∈w}{\mathrm{\Σ}}{z}_1\left(x_{j,k}\right)z_2\left(x_{j,k}\right)\\ \underset{x_{j,k}\∈w}{\mathrm{\Σ}}{z}_1\left(x_{j,k}\right)z_2\left(x_{j,k}\right)& \underset{x_{j,k}\∈w}{\mathrm{\Σ}}{z}_2\left(x_{j,k}\right)z_2\left(x_{j,k}\right)\end{array}\right]$

z₁(·)表示垂直梯度，z₂(·)表示水平梯度。

步骤4，依据步骤3得到的局部自适应协方差矩阵C，计算以像素点x为中心的邻域窗w平滑矩阵H，

其中，det(C)表示自适应协方差矩阵C的行列式，h为一平滑参数，通常取值为[4.0,4.5]，通过实验验证，在本发明中h＝4.3时得到的实验结果最理想。

步骤5，依据计算得到的平滑矩阵H，利用权重影响，求得步骤1中的像素点x的相干斑抑制结果

表示邻域窗w内所有像素点的值的均值，k＝1/k₁是一个权重数值，

表示邻域窗w内像素点的值的方差，

表示SAR图像的相干斑方差，对于给定的强度格式的SAR图像，该SAR图像的相干斑方差为：

而对于给定的幅度格式的SAR图像，SAR图像的相干斑方差为：L表示图像的视数。

步骤6，对步骤1中所述的SAR图像每一个像素点按步骤1—5进行处理，得到每一个像素点的相干斑抑制结果，从而得到所述SAR图像的相干斑抑制结果。

步骤7，将步骤6得到的结果，按照步骤1—步骤6进行迭代处理，迭代次数一般为[8,10]次，得到最终的所述SAR图像的相干斑抑制结果。通过实验验证，在本发明中迭代9次时得到的实验结果最理想。

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

1.实验条件

实验仿真环境：MATLAB2010a，Intel(R)Pentium(R)1CPU2.4GHz,Window XPProfessional。

实验仿真图像：barbara自然图像模拟的SAR图像，添加了视数为4的幅度相干斑，图像大小为512×512；光学图像模拟的SAR图像，添加了视数为8的幅度相干斑，图像大小为512×512；真实的幅度SAR图像是来自英国国家防务局SAR，X波段,3m分辨率，英格兰Bedfordshire的一处农田，图像视数为3，图像大小为256×256；真实的强度SAR图像是来自英国国家防务局SAR，X波段，3m分辨率，英国Bedford附近，图像视数为6，图像大小为256×256。

实验对比方法：增强Lee滤波，小波软阈值滤波以及西安电子科技大学智能感知与图像理解教育部重点实验室的凤宏晓提出的改进非局部均值滤波，前两种方法是SAR图像空域滤波和频域滤波中最具代表性的方法，第三种方法代表了目前SAR图像相干斑抑制的前沿技术。

实验结果评价：对于模拟的SAR图像用峰值性噪比PSNR和结构相似性指数SSIM，两个指标反映噪声抑制的程度和对图像结构的保持，PSNR越大对相干斑抑制效果越好，SSIM的理想值是1，越接近1对图像结构的保持性越好，对于真实的SAR图像还可以用平滑区域的等效视数ENL作为评价标准，ENL越大对区域的平滑效果越好，用主观的视觉效果也可以对相干斑抑制效果的好坏进行判断。

2.实验内容与结果

仿真1，利用现有的增强Lee方法、小波软阈值方法、改进非局部均值方法及本发明方法，对图2所示的barbara自然图像模拟的SAR图像进行滤波，滤波后的结果如图6所示，其中，图6(a)是增强Lee的滤波结果图，图6(b)是小波软阈值的滤波结果图，图6(c)改进非局部均值的滤波结果图，图6(d)是本发明的滤波结果图。

从图6可以看出，增强Lee滤波和小波软阈值滤波的效果都不太理想，尤其对纹理的保持无法令人满意，改进的非局部均值滤波效果大大提升，但纹理保持上仍然不及本发明。

仿真2，利用现有的增强Lee方法、小波软阈值方法、改进非局部均值方法及本发明方法，对图3所示的光学图像模拟的SAR图像进行滤波，滤波后的结果如图7所示，其中，图7(a)是增强Lee的滤波结果图，图7(b)是小波软阈值的滤波结果图，图7(c)改进非局部均值的滤波结果图，图7(d)是本发明的滤波结果图。

从图7可以看出，增强Lee滤波与小波软阈值滤波在纹理细节信息保留上效果不理想并将边缘模糊化，本发明在边缘细节和纹理信息的保持上与改进非局部均值方法相差无几。

仿真3，利用现有的增强Lee方法、小波软阈值方法、改进非局部均值方法及本发明方法，对图4所示的真实的幅度SAR图像进行滤波，滤波后的结果如图8所示，其中，图8(a)是增强Lee的滤波结果图，图8(b)是小波软阈值的滤波结果图，图8(c)是改进非局部均值的滤波结果图，图8(d)是本发明的滤波结果图。

从图8可以看出，增强Lee滤波虽然能起到一定的噪声抑制能力，但是对边缘和纹理细节的保持上不尽人意，小波软阈值的方法存在划痕的问题，边缘保持略好于增强Lee滤波，改进的非局部均值无论是对匀质区域的噪声抑制，还是对边缘的保持都是非常优秀的，对点目标的保持也不错，视觉上本发明的方法与改进的非局部均值相差无几。

仿真4，利用增强Lee,小波软阈值及本发明对图6所示的真实的强度SAR图像进行滤波，滤波后的结果如图9所示，其中，图9(a)是增强Lee的滤波结果图，图9(b)是小波软阈值的滤波结果图，图9(c)是本发明的滤波结果图。

从图9可以看出，用本发明方法滤波后的结果，其纹理、边缘、点目标都清晰可见，对匀质区域的噪声抑制也十分理想，而其它两种方法，总是存在一些令人无法接受的问题。

将不同方法噪声抑制后的模拟SAR图像的峰值性噪比PSNR和结构相似性指数SSIM进行比较，结果如表1。

表1不同方法噪声抑制后的模拟SAR图像的PSNR和SSIM进行比较

表1数据表明，将本发明在应用于不同的图像中峰值性噪比PSNR和结构相似性指数SSIM都是最优的。

将不同噪声抑制方法后的真实幅度SAR图像的平滑区域的均值、方差和等效视数ENL进行比较，结果如表2。

表2不同方法噪声抑制后的真实幅度SAR图像的平滑区域的均值、方差和ENL

表2中区域1，区域2，区域3是图5中标注的三个区域.

从表2可见，对幅度SAR图像的相干斑抑制，增强Lee滤波方法在均值保持上是最好的，而本发明也非常接近，主观视觉效果良好的改进非局部均值在均值保持上明显有较大的偏差，在标准差上本发明在这些方法中一般都是最优的，改进的非局部均值次之，增强Lee滤波和小波软阈值滤波都比较差，在最终的等效视数ENL上本发明是所有方法中最好的。

综上所述，本发明提出的基于核回归的SAR图像相干斑抑制方法，能够很好的保持SAR图像的边缘信息、纹理细节和点目标，并且对平滑区域的平滑效果也很理想，因此本发明对SAR图像相干斑的抑制效果显著。

Claims (9)

1.基于核回归的SAR图像相干斑抑制方法，包括如下步骤：

步骤1，对一幅SAR图像，进行边界镜像对称扩展，扩展半径为N，取图像中的一个像素点x，以该像素点为中心取(2×N+1)大小的邻域窗w；

步骤2，用索贝尔（Sobel）算子计算邻域窗内每一个像素点的水平梯度与垂直梯度；

步骤3，根据步骤2得到的水平梯度与垂直梯度计算以像素点x为中心的邻域窗w的自适应协方差矩阵C；

步骤4，依据步骤3得到的局部自适应协方差矩阵C，计算以像素点x为中心的邻域窗w平滑矩阵H；

步骤5，依据步骤4得到的平滑矩阵H，求得步骤1中的像素点x的相干斑抑制结果；

步骤6，对步骤1中所述的SAR图像每一个像素点按步骤1-5进行处理，得到每一个像素点的相干斑抑制结果，从而得到所述SAR图像的相干斑抑制结果；

步骤7，将步骤6得到的结果，按照步骤1-步骤6进行迭代处理，得到最终的所述SAR图像的相干斑抑制结果。

2.如权利要求1所述的SAR图像相干斑抑制方法，其中，N=3。

3.如权利要求1所述的SAR图像相干斑抑制方法，其中，用索贝尔（Sobel）算子计算邻域窗内每一个像素点的水平梯度与垂直梯度具体包括：

2a)设置水平梯度算子模板与垂直梯度算子模板分别用来检测图像的垂直边缘及水平边缘，模板内的数字为模板系数，水平梯度算子模板的方向与垂直边缘的方向相垂直，垂直梯度算子模板的方向与水平边缘的方向相垂直；

$\frac{1}{4}*\left[\begin{array}{ccc}-1& -2& -1\\ 0& 0& 0\\ 1& 2& 1\end{array}\right],\frac{1}{4}*\left[\begin{array}{ccc}-1& 0& 1\\ -2& 0& 2\\ -1& 0& 1\end{array}\right]$

水平梯度算子模板     垂直梯度算子模板

将水平梯度算子模板与垂直梯度算子模板表达为：

${\left[M\right]}^i=\left[\begin{array}{ccc}{M}_{-1,-1}^i& M_{-\mathrm{1,0}}^i& M_{-\mathrm{1,1}}^i\\ M_{0,-1}^i& M_{\mathrm{0,0}}^i& M_{\mathrm{0,1}}^i\\ M_{1,-1}^i& M_{1.0}^i& M_{\mathrm{1,1}}^i\end{array}\right]$

式中，i=1,2分别代表水平梯度算子模板，垂直梯度算子模板；

代表水平梯度算子模板或垂直梯度算子模板内的中心位置值；

2b)步骤1中的邻域窗w内每一个像素点的垂直梯度与水平梯度的计算公式为：

$z_i\left(x_{j,k}\right)=\underset{m=-1}{\overset{1}{\mathrm{\Σ}}}\underset{n=-1}{\overset{1}{\mathrm{\Σ}}}{x}_{j,k}{M}_{m,n}^i$

i=1,2，z₁(x_j,k),z₂(x_j,k)分别是邻域窗w内当前处理的像素点x_j,k的垂直梯度和水平梯度，其中j,k的取值范围均为[1,2×N+1]。

4.如权利要求1所述的SAR图像相干斑抑制方法，其中，

$C\≈\left[\begin{array}{cc}\underset{x_{j,k}\∈w}{\mathrm{\Σ}}{z}_1\left(x_{j,k}\right)z_1\left(x_{j,k}\right)& \underset{x_{j,k}\∈w}{\mathrm{\Σ}}{z}_1\left(x_{j,k}\right)z_2\left(x_{j,k}\right)\\ \underset{x_{j,k}\∈w}{\mathrm{\Σ}}{z}_1\left(x_{j,k}\right)z_2\left(x_{j,k}\right)& \underset{x_{j,k}\∈w}{\mathrm{\Σ}}{z}_2\left(x_{j,k}\right)z_2\left(x_{j,k}\right)\end{array}\right]$

z₁(·)表示垂直梯度，z₂(·)表示水平梯度。

5.如权利要求1所述的SAR图像相干斑抑制方法，其中，

其中，det(C)表示自适应协方差矩阵C的行列式，h为一平滑参数，取值为[4.0,4.5]。

6.如权利要求1所述的SAR图像相干斑抑制方法，其中，步骤5求得的像素点x的相干斑抑制结果为

表示邻域窗w内所有像素点的值的均值，k＝1/k₁是一个权重数值，

表示邻域窗w内像素点的值的方差，

表示SAR图像的相干斑方差。

7.如权利要求1所述的SAR图像相干斑抑制方法，其中，对于给定的强度格式的SAR图像，该SAR图像的相干斑方差为：

或者，

对于给定的幅度格式的SAR图像，SAR图像的相干斑方差为：

表示图像的视数。

8.如权利要求1所述的SAR图像相干斑抑制方法，其中，迭代次数为[8,10]次。

9.如权利要求1所述的SAR图像相干斑抑制方法，其中，h＝4.3。

Images