CN104103046A - 一种极化sar图像的快速滤波处理方法 - Google Patents

Abstract

本发明涉及一种极化SAR图像快速滤波处理方法，可有效提高极化SAR图像滤波处理速度，节约费用的问题，技术方案是，包括以下步骤：1)载入极化SAR图像数据，并进行扩展，为滤波做准备；2)利用扩展后的极化SAR图像计算得到总功率图像；3)基于非局部均值中搜索窗Ω，利用积分图像计算实现快速滤波处理；4)使用Pauli向量法对滤波后的整个极化协方差矩阵C合成伪彩图，输出滤波结果；本发明大大降低了遍历计算范围，有效减少了计算开销，提高了极化SAR滤波处理速度；扩展了非局部均值滤波处理的适用范围，利用积分图，一方面是相似性度量的计算次数，另一方面是矩形窗口积分计算，对SAR图像的相干斑抑制能得到比现有技术更好的效果。

Description

一种极化SAR图像的快速滤波处理方法

技术领域

本发明涉及图像处理技术领域，特别是一种极化SAR图像快速滤波处理方法，可用于极化SAR图像的场景分析和图像理解的预处理。 

背景技术

合成孔径雷达(Synthetic Aperture Radar,SAR)就是利用雷达与目标的相对运动把尺寸较小的真实天线孔径用数据处理的方法合成一较大的等效天线孔径的雷达，也称综合孔径雷达。SAR不受天气、光照等条件的限制，能够连续、实时地进行对地观测。极化SAR为遥感图像中的信息获取和处理提供了更加丰富的目标信息，使其在采集地表或地面覆盖物的物理和电磁结构信息的应用中起着越来越重要的作用。由于极化SAR受相干斑噪声严重的干扰，相干斑噪声不仅导致数据质量下降，同时给极化SAR数据分类、目标检测等后处理带来困难。因此相干斑噪声的滤波处理成为极化SAR信息处理中的首要问题，相干斑噪声如果能够有效的抑制，可以显著提高极化SAR数据的可读性，有利于图像的处理、解译和信息提取。 

非局部(non local means，NLM)算法在近年来受到广泛关注和应用，其基本原理是利用图像块的相似性作为均值滤波的加权权重。这种算法对于自然图像的去噪结果十分理想，尤其是在边缘纹理等细节信息的保持方面明显优于其它现有的方法。目前该去噪方法正逐步推广到极化SAR相干斑噪声滤波处理领域。但是非局部均值算法的一个重要问题是计算量问题，在NLM类算法中像素点间相似性度量的计算最为耗时，因为任意两个图像块之间的相似性，都要通过计算对应像素之间的差异，然后进行累积获得。可以想象，随着滤波窗口的滑动，像素之间相似性度量运算存在大量的重复操作，使现有非局部均值类算法很难应用于工程实践，适用面小，计算速率低，费用高，为此需要提供极化SAR的非局部均值滤波类快速算法，但是今未见有公开报导。 

发明内容

针对上述情况，为克服现有技术之缺陷，本发明之目的就是提供一种极化SAR图像快速滤波处理方法，以拓展非局部均值滤波方法的适用范围，有效提高极化SAR图像滤波处理速度，节约费用的问题 

本发明解决的技术方案是，首先使用不局限于高斯分布的相似性度量计算滤波权重，然后利用积分图在极化SAR功率图上实现相似性度量的快速计算，最后用该相似度分别对极化散射矩阵或者协方差矩阵进行滤波处理，包括以下步骤： 

1)载入极化SAR图像数据，并进行扩展，为滤波做准备； 

2)利用扩展后的极化SAR图像计算得到总功率图像； 

3)记Δt表示目标像素t所在的邻域窗，Δs表示中心像素s所在的邻域窗，基于非局部均值中搜索窗Ω，利用积分图像计算实现快速滤波处理，步骤如下： 

3a)初始化邻域窗Δt在搜索窗内的当前位置 

3b)在扩展的总功率图像上，分别记A为邻域窗Δs中心像素构成的全局图像，B为邻域窗Δt中心像素构成的全局图像，按照下式计算全局图像A和B对应像素的相似性度量DS，其中I_A为图像A内像素的灰度值，I_B为图像B内像素的灰度值，(x,y)为像素坐标： 

$\mathrm{DS}(x,y)=\ln (\frac{I_A(x,y)}{I_B(x,y)}+\frac{I_B(x,y)}{I_A(x,y)})$

3c)按照下式构建积分图SI： 

$\mathrm{SI}(x,y)=\underset{u\≤x,v\≤y}{\mathrm{\Σ}}\mathrm{DS}(u,v)$

其中，u、v分别为DS图像中左上角为原点与右下角为点(x,y)区域内点的横纵坐标； 

3d)根据下式计算邻域窗Δt与Δs的相似性度量： 

SSI(x,y)＝SI(x,y)+SI(x+2p+1,y+2p+1)-SI(x+2p+1,y)-SI(x,y+2p+1) 

3e)通过以下计算得到目标像素t对中心像素s的滤波权值w(x,y)，其中h为滤波参数： 

$w(x,y)=\exp \left(\frac{\mathrm{SSI}(x,y)}{h^2}\right)$

3f)对极化协方差矩阵C每个通道分别进行加权滤波，为滤波结果，滤波式为： 

$\hat{C}(x,y)=\underset{x,y\∈\mathrm{\Ω}}{\mathrm{\Σ}}w(x,y)C(x,y)$

3g)利用下式对滤波结果进行归一化，其中为归一化系数： 

$\hat{C}(x,y)=\frac{\hat{C}(x,y)}{Z}$

3h)遍历整个搜索窗，完成滤波处理； 

4)使用Pauli向量法对滤波后的整个极化协方差矩阵C合成伪彩图，输出滤波结果。 

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

(1)将整个图像范围内的遍历计算巧妙转化为搜索窗内的遍历计算，大大降低了遍历计算范围，有效减少了计算开销，提高了极化SAR滤波处理速度； 

(2)使用扩展的相似性度量函数，不再局限于高斯分布的极化SAR图像数据，对其他统计特征的极化SAR图像数据也适用，从而扩展了本发明方法的适用对象； 

(3)利用积分图实现了邻域窗相似性度量函数的快速计算，只需遍历一次原图像，计算开销很小，进一步大大提高了极化SAR图像滤波处理速度； 

(4)本发明采用SPAN数据来计算滤波权值。SPAN数据基本包含了极化SAR数据的所有极化特性，本发明选它来计算滤波权值，能够保持极化SAR数据的极化相关性； 

(5)本发明使用了扩展的相似性度量函数，该相似性度量不再局限于高斯分布的极化SAR图像数据，对其他统计特征的极化SAR图像数据也适用，从而扩展了非局部均值滤波处理的适用范围； 

(6)本发明利用积分图实现了极化SAR图像滤波的快速处理。利用积分图可以方便地计算出一个图像块的灰度值累加和，只需遍历一次原图像，计算开销很小。本发明中利用积分图主要是在两方面减少了重复运算，一方面是相似性度量的计算次数，另一方面是矩形窗口积分计算； 

(7)本发明对SAR图像的相干斑抑制能得到比现有技术更好的效果。 

附图说明：

图1是本发明的流程图； 

图2是本发明使用的搜索窗、邻域窗的相互关系示意图； 

图3是本发明利用积分图计算图像块累加和的示意图； 

图4是本发明与已有方法对实测极化SAR数据的滤波效果对比图，其中：(a)为滤波前的原始数据图；(b)Refine-Lee滤波(精致极化Lee滤波)处理后的数据图；(c)非局部均值滤波处理后的数据图；(d)本发明滤波处理后的数据图； 

图5是本发明滤波处理的实测极化SAR数据的统计特征图。 

具体实施方式

以下结合附图对本发明实施例进一步说明： 

如附图1所示，为本发明对极化SAR图像进行快速滤波处理的流程图，利用本发明，进行极化SAR图像快速滤波的处理流程如下： 

1)载入极化SAR图像数据协方差矩阵C，图像大小为m×n×9(单位为像素)，其中  $C=\left[\begin{array}{ccc}{C}_{11}& C_{12}& C_{13}\\ C_{21}& C_{22}& C_{23}\\ C_{31}& C_{32}& C_{33}\end{array}\right],$ C_ij为单通道图像，i和j表示雷达发射与接收极化组合方式，C₁₁、C₁₂、C₁₃、C₂₁、C₂₂、C₂₃、C₃₁、C₃₂、C₃₃为式中9个单通道图像元素； 

2)将上述9个单通道图像尺寸分别扩展至(m+2W+2P)×(n+2W+2P)，边缘填充最相邻像素 灰度值，得到新图像C′_ij，其中W为搜索窗Ω长度的一半并向下取整数，P为邻域窗Δs长度的一半并向下取整数； 

3)利用扩展后的极化SAR图像计算得到总功率图像：SPAN′＝C′₁₁+C′₂₂+C′₃₃；其中C₁₁′、C₁₁′、C₁₁′为扩展后的极化SAR的3个图像元素； 

4)基于搜索窗Ω，利用积分图像计算实现快速滤波处理，记Δt表示目标像素所在的邻域窗，Δs表示中心像素所在的邻域窗，如图2所示： 

4a)邻域窗Δt(的中心像素)在搜索窗内的当前坐标(dx，dy)初始化为(-W，-W)； 

4b)在SPAN’图像上，选取坐标范围为(W：m+2P+W，W：n+2P+W)的像素构成邻域窗Δs的全局图像A，选取坐标范围为(W+dx：m+2P+W+dx，W+dy：n+2P+W+dy)的像素构成邻域窗Δt的全局图像B，按照下式计算全局图像A和B对应像素的相似性度量图像DS； 

$\mathrm{DS}(x,y)=\ln (\frac{I_A(x,y)}{I_B(x,y)}+\frac{I_B(x,y)}{I_A(x,y)})$

其中，I_A为图像A内像素的灰度值，I_B为图像B内像素的灰度值，(x,y)为像素坐标。 

4c)按照下式构建积分图SI并存储； 

$\mathrm{SI}(x,y)=\underset{u\≤x,v\≤y}{\mathrm{\Σ}}\mathrm{DS}(u,v)$

4d)如图3所示，根据下式计算邻域窗Δt与Δs的相似性度量，即对邻域窗内所有像素相似性度量的求和： 

SSI(x,y)＝SI(x,y)+SI(x+2p+1,y+2p+1)-SI(x+2p+1,y)-SI(x,y+2p+1) 

4e)然后通过以下计算得到目标像素t对中心像素s的滤波权值w(x,y)，为了能够较好地兼顾纹理保持和噪声抑制效果，这里参数h取3倍的SPAN图噪声标准差： 

$w(x,y)=\exp \left(\frac{\mathrm{SSI}(x,y)}{h^2}\right)$

4f)对极化协方差矩阵C的9个单通道图像分别进行加权滤波，为滤波结果，滤波式为： 

${\hat{C}}_i(x,y)=\underset{x,y\∈\mathrm{\Ω}}{\mathrm{\Σ}}w(x,y)C_i(x,y)$

4g)利用下式对滤波结果进行归一化，其中为归一化系数，可通过累加求和取得： 

$\hat{C}(x,y)=\frac{\hat{C}(x,y)}{Z}$

4h)遍历整个搜索窗，直到邻域窗Δt在搜索窗内的坐标(dx，dy)为(W，W)；完成滤波处理； 

5)使用Pauli向量法对滤波后的整个极化协方差矩阵C合成伪彩图，以观察滤波的效果。Pauli向量法主要是使用协方差矩阵C中的C₁₁，C₂₂和C₃₃三个元素来合成伪彩图。 

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

1.仿真条件 

硬件平台：Intel(R)Xeon(R)CPU X56702.93GHz CPU，12GB内存，Windows764位操作系统。 

编程工具为MATLAB7.10(R2010a,32bit)。 

实验数据来源于AIRSAR四视San Francisco区域极化SAR图像，大小为300×300(单位为像素)，原始图像见图4(a)。 

2.仿真内容 

应用本发明和极化Refine-Lee滤波方法和原始NLM滤波方法分别对一幅AIRSAR图像进行去斑实验，从同质区域的内部一致性，边缘信息和纹理信息的完整性以及点目标的清晰性来评价这些方法的性能。 

各算法参数取值如下：Refine-Lee滤波的窗口大小取7×7、衰减系数D取1；非局部均值和本发明方法的搜索窗口大小为21×21、相似性窗口大小为7×7，窗口大小的单位均为像素，平滑参数h取3.5。评价结果分别用其同质区域的等效视数ENL、边缘保持指数EPI以及运算时间的大小来衡量滤波结果的好坏。 

3.仿真结果 

分别计算图4(a)中A处同质区域的等效视数ENL，B处纹理区域的边缘保持指数EPI及统计各滤波算法的运行消耗时间见表1。 

表1 等效视数、边缘保持指数及运行时间的对比 

图4(b)为Refine-Lee滤波结果，图4(c)是非局部均值滤波结果，图4(d)为本发明方法滤波 结果。可以发现图4(b)中精致极化Lee滤波结果在平滑效果上不是很理想，整体上出现了人为的块状结构，使得边缘也表现的比较模糊；在均质海洋区域A对比不同滤波结果，可以看出图4(c)中非局部均值滤波结果要明显优于精致极化Lee滤波，而图4(d)中本发明方法滤波结果视觉平滑效果更为清晰。这与表1中等效视数大小正好吻合。 

在建筑物区域B对比不同滤波结果，可以看出精致极化Lee滤波严重破坏了原来的纹理细节，而非局部均值和本发明方法滤波结果房屋较为清晰。同时，观察区域D可以看出，两条类海堤建筑物在本发明方法处理结果清晰完整的显现出来，而在前两个算法处理结果不是丢失，就是模糊。对照表1中边缘保持指数值，充分反应出本发明方法具有很好纹理结构的保持能力。 

运算效率方面，表1中运行时间表明在该实验条件下，本发明方法相比原始非局部均值算法而言，计算速度提高了2000多倍，极大改善了非局部均值算法在极化SAR中应用效率，为该算法的工程应用提供了基础，具有很强的实用价值。 

为了验证本发明方法的适用性，申请人记录了实验过程中极化SAR功率的直方图，如图5所示。从中可以看出极化功率图的直方图近似服从指数分布，远不是理想的高斯分布，所以本发明方法满足于各类噪声分布条件下滤波，较好解决了原始非局部均值算法中高斯分布假设的局限性，应用面广。 

Claims (1)

1.一种极化SAR图像的快速滤波处理方法，其特征在于，包括以下步骤：

1)载入极化SAR图像数据，并进行扩展，为滤波做准备；

2)利用扩展后的极化SAR图像计算得到总功率图像；

3)记Δt表示目标像素t所在的邻域窗，Δs表示中心像素s所在的邻域窗，基于非局部均值中搜索窗Ω，利用积分图像计算实现快速滤波处理，步骤如下：

3a)初始化邻域窗Δt在搜索窗内的当前位置；

3b)在扩展的总功率图像上，分别记A为邻域窗Δs中心像素构成的全局图像，B为邻域窗Δt中心像素构成的全局图像，按照下式计算全局图像A和B对应像素的相似性度量DS，其中I_A为图像A内像素的灰度值，I_B为图像B内像素的灰度值，(x,y)为像素坐标：

$\mathrm{DS}(x,y)=\ln (\frac{I_A(x,y)}{I_B(x,y)}+\frac{I_B(x,y)}{I_A(x,y)})$

3c)按照下式构建积分图SI：

$\mathrm{SI}(x,y)=\underset{u\≤x,v\≤y}{\mathrm{\Σ}}\mathrm{DS}(u,v)$

其中，u、v分别为DS图像中左上角为原点与右下角为点(x,y)区域内点的横纵坐标；

3d)根据下式计算邻域窗Δt与Δs的相似性度量：

SSI(x,y)＝SI(x,y)+SI(x+2p+1,y+2p+1)-SI(x+2p+1,y)-SI(x,y+2p+1)

3e)通过以下计算得到目标像素t对中心像素s的滤波权值w(x,y)，其中h为滤波参数：

$w(x,y)=\exp \left(\frac{\mathrm{SSI}(x,y)}{h^2}\right)$

3f)对极化协方差矩阵C每个通道分别进行加权滤波，为滤波结果，滤波式为：

$\hat{C}(x,y)=\underset{x,y\∈\mathrm{\Ω}}{\mathrm{\Σ}}w(x,y)C(x,y)$

3g)利用下式对滤波结果进行归一化，其中为归一化系数：

$\hat{C}(x,y)=\frac{\hat{C}(x,y)}{Z}$

3h)遍历整个搜索窗，完成滤波处理；

4)使用Pauli向量法对滤波后的整个极化协方差矩阵C合成伪彩图，输出滤波结果。