CN102903080B - 合成孔径雷达图像相干斑噪声抑制性能的非监督评估方法 - Google Patents

Abstract

本发明公开的合成孔径雷达图像相干斑噪声抑制性能的非监督评估方法，具体按照以下步骤实施：1）利用各项异性高斯核平行窗通过比率运算生成合成孔径雷达SAR图像的ESM与方向信息；2）获得合成孔径雷达SAR图像的ESM阈值估计；3）根据ESM阈值实现合成孔径雷达SAR图像边缘与非边缘区域的自动划分；4）利用非极小值抑制方法从合成孔径雷达SAR图像边缘区域提取细化边缘；5）获得等效视数ENL估计和局部边缘保持指数EKI估计；6）完成SAR图像基于ENL与EKI参数的相干斑抑制性能的评估。实现SAR图像相干斑抑制性能的全程自适应非监督评价，并有效克服传统方法评估不全面、受区域选择影响大的问题，使相干斑抑制性能得到更加客观稳健的全面评估。

Description

合成孔径雷达图像相干斑噪声抑制性能的非监督评估方法

技术领域

本发明属于合成孔径雷达图像处理技术领域，涉及一种合成孔径雷达图像相干斑噪声抑制性能的非监督评估方法。

背景技术

合成孔径雷达借助大量随机分布的散射体反射的雷达回波相干叠加成像，从而不可避免地在合成孔径雷达SAR图像中产生一种称为相干斑的乘性噪声。相干斑噪声的存在严重降低了合成孔径雷达SAR图像的视觉质量，限制了特征提取、目标识别等后续解译处理技术的有效性。因此，相干斑抑制对于改进合成孔径雷达SAR图像成像质量，提高后续解译处理效果都具有重要意义。

当前，合成孔径雷达SAR图像相干斑抑制方法已形成空域滤波、变换域滤波与偏微分扩散滤波三大类，而合成孔径雷达SAR图像抑斑性能的有效评估不仅是评价抑斑算法优劣的主要手段，而且也为抑斑合成孔径雷达SAR图像的后续解释处理提供了必要的先验信息。

在合成孔径雷达SAR图像抑斑性能的评价参数中，反映抑斑程度的等效视数(the Equivalent Number of Looks，ENL)与衡量边缘保持能力的边缘保持指数(Edge Keeping Index，EKI)是最重要的两类指标。但传统方法在利用ENL与EKI两类指标评价合成孔径雷达SAR图像抑斑性能时，往往采用监督评估方法，通过人为选择几块合成孔径雷达SAR图像同质区来估计ENL参数，而通过选择几段边缘区域来估计EKI参数。这类方法由于需要人为选择评估区域，势必造成评估结果主观性较强且受评估区域选择影响大，同时，由于人为选择大量的评估区麻烦且不现实，势必存在评估不全面的问题，而对不同合成孔径雷达SAR图像不能进行自适应的自动评估，也为评估工作的实施带来了不便。

发明内容

本发明的目的在于提供一种合成孔径雷达图像相干斑噪声抑制性能的非监督评估方法，实现SAR图像相干斑抑制性能的全程自适应非监督评价，并有效克服传统方法评估不全面、受区域选择影响大的问题，使相干斑抑制性能得到更加客观稳健的全面评估。

本发明所采用的技术方案是，合成孔径雷达图像相干斑噪声抑制性能的非监督评估方法，体按照以下步骤实施：

步骤1、利用各项异性高斯核平行窗通过比率运算生成合成孔径雷达SAR图像的边缘强度映射dge strength map简称ESM与方向信息：

采用带有θ方向的各向异性高斯核平行窗按照公式（1）计算各向异性高斯平行窗权重系数，然后生成合成孔径雷达SAR图像的边缘强度映射ESM与方向信息：

带有θ方向的各向异性高斯核平行窗计算公式如下：

其中：r_x＝xcosθ-ysinθ，r_y＝xsinθ+ycosθ，x与y分别表示水平与垂直方向坐标，θ表示方向，k表示两个平行窗的间距，σ_x与σ_y分别表示x与y方向的标准差，r_x,r_y表示对x，y按照θ方向旋转后形成的新坐标；

步骤2、获得合成孔径雷达SAR图像的边缘强度映射ESM阈值T：

经步骤1获得合成孔径雷达SAR图像的边缘强度映射ESM与方向信息后，开始对边缘强度映射ESM图像进行多次重复阈值递减的阈值化处理，统计阈值化处理后的边缘强度映射ESM图像中连通像素点集合的个数，确定合成孔径雷达SAR图像区域划分的阈值T；

步骤3、根据ESM阈值实现合成孔径雷达SAR图像边缘与非边缘区域的自动划分：

将步骤2中获得的边缘强度映射ESM阈值T代入公式（2）中，即可获得合成孔径雷达SAR图像边缘区域与同质区域的自动划分；

其中，ψ若为非零的区域为即为边缘区域，ψ若为零的区域就视为为同质区；

步骤4、利用非极小值抑制方法从合成孔径雷达SAR图像边缘区域提取细化边缘；

步骤5、获得等效视数ENL估计和局部边缘保持指数EKI估计：

经步骤4在合成孔径雷达SAR图像区域划分基础上，非边缘区域内各像素局部ENL估计要通过加窗实现，对像素C(m,n)施加大小为N的窗R，计算像素C(m,n)的局部ENL估计；

利用非极小值抑制方法获得合成孔径雷达SAR图像细化边缘基础上，采用公式对坐标为(x,y)的任意一个边缘像素估计其局部EKI；

步骤6、利用统计方法将各像素局部ENL与EKI估计形成对比曲线，完成SAR图像基于ENL与EKI参数的相干斑抑制性能的评估：

经步骤5获得的大量局部ENL与EKI估计值，利用统计方法获得不同抑斑算法下抑斑SAR图像ENL与EKI参数的对比分布，进而对这两类指标进行全面且直观的总体评估。

本发明的特点还在于，

步骤1具体按照以下步骤实施：

采用带有θ方向的各向异性高斯核平行窗，各向异性高斯核平行窗的方向分别：方向为θ=0°的各向异性高斯平行窗，方向为θ=45°的各向异性高斯平行窗，方向为θ=90°的各向异性高斯平行窗，方向为θ=135°的各向异性高斯平行窗；

利用公式（1）生成的各向异性高斯平行窗权重系数，再将各向异性高斯核平行窗内各像素与待估计中心像素的距离作为权重分配的依据，实现加权平均，即可获得合成孔径雷达SAR图像的边缘强度映射ESM与方向信息，获得更好的SAR图像ESM估计与区域划分。

步骤2具体按照以下步骤实施：

1）经步骤1获得合成孔径雷达SAR图像的边缘强度映射ESM与方向信息之后，开始对边缘强度映射ESM图像进行多次重复阈值递减的阈值化处理，即阈值T从最大值1开始递减，每次减少量为0.01，直至阈值T为0.01为止；

2）每次阈值化处理后，统计阈值化处理后的边缘强度映射ESM图像中连通像素点集合的个数，即可获得不同阈值T下阈值化边缘强度映射ESM图像中连通像素点集合个数曲线，当边缘强度映射ESM图像中连通像素点集合个数曲线由最大值迅速降低直至出现第一次极小值时，此时具有最好的图像区域划分性能，这个极小值就是合成孔径雷达SAR图像区域划分的阈值T。

骤4具体按照以下步骤实施：

利用非极小值抑制方法从合成孔径雷达SAR图像边缘区域提取细化边缘：

1）若ψ(x，y)＞0，则在ψ域中以位置(x,y)为中心的3×3窗内，沿与θ(x,y)垂直的方向取两相邻值ψ(x₁,y₁)与ψ(x₂,y₂)；

2）若设ψ(x,y)＜ψ(x₂,y₂)与ψ(x,y)＜ψ(x₁,y₁)同时成立为条件Ω，则边缘抽取可表示为：

其中，Edge(x,y)＝1的像素为边缘区域中的边缘像素，否则为边缘区域的非边缘像素。

步骤5具体按照以下步骤实施：

经步骤4在合成孔径雷达SAR图像区域划分基础上，非边缘区域内各像素局部ENL估计要通过加窗实现，对像素C(m,n)施加大小为N的窗R，像素C(m,n)的局部ENL估计采用公式（4）计算，窗R相对于整个图像是局部窗，计算出的是窗R中各元素的平均值，称为窗内局部均值E，具体的计算公式如下：

$E(m,n)=\frac{1}{N}\underset{(i,j)\∈R}{\mathrm{\Σ}}I(m+i,n+j),---\left(4\right)$

窗中各元素的方差称为窗内局部方差D可表示为：

$E(m,n)=\frac{1}{N-1}\underset{(i,j)\∈R}{\mathrm{\Σ}}{[I(m+i,n+j)-E(m,n)]}^2,---\left(5\right)$

经步骤4中由阈值化边缘强度映射ESM生成的区域划分信息ψ，对矩形窗R内各像素施加约束，将矩形窗R中与窗中心像素C(m,n)不同质的像素剔除，形成不规则窗；再通过给矩形窗R内各元素加权的方法间接实现不规则加窗运算；对于矩形窗R，窗中心像素C(m,n)属于某同质区，则像素I(i,j)∈R对应的权值δ(i,j)可表示为：

可将公式（4）调整为：

$E(m,n)=\frac{\underset{(i,j)\∈R}{\mathrm{\Σ}}I(m+i,n+j)\mathrm{\δ}(m+i,n+j)}{\underset{(i,j)\∈R}{\mathrm{\Σ}}\mathrm{\δ}(i,j)},---\left(7\right)$

将公式（5）调整为：

$D(m,n)=\frac{\underset{(i,j)\∈R}{\mathrm{\Σ}}{[I(m+i,n+j)-E(m,n)]}^2\mathrm{\δ}(m+i,n+j)}{\underset{(i,j)\∈R}{\mathrm{\Σ}}\mathrm{\δ}(i,j)-1},---\left(8\right)$

利用非极小值抑制方法获得合成孔径雷达SAR图像细化边缘基础上，对坐标为(x,y)的任意一个边缘像素估计其局部EKI时，采用如下的公式（9）估计：

$V_{\mathrm{EKI}}(x,y)=G_{{\mathrm{\θ}}^{\⊥}}^{\′}(x,y)/G_{{\mathrm{\θ}}^{\⊥}}(x,y)---\left(9\right)$

其中，和分别表示抑斑前、抑斑后合成孔径雷达SAR图像位于(x,y)处的边缘像素沿与该像素方向θ(x,y)垂直方向的梯度模值。

步骤6具体按照以下步骤实施：

经步骤5获得的大量局部ENL与EKI估计值，利用统计方法获得不同抑斑算法下抑斑SAR图像ENL与EKI参数的对比分布，进而对这两类指标进行全面且直观的总体评估：

ENL参数的统计方式是：

1）将各抑斑算法对应的抑斑合成孔径雷达SAR图像局部ENL，在分贝单位下，先进行比率运算，具体是利用公式（1）生成的各向异性高斯平行窗权重系数，计算两个平行窗内的局部均值，然后将这两个局部均值直接进行比率运算；

2）再按各像素在合成孔径雷达SAR图像中的实际位置顺序，统计其对比分布；

EKI参数的统计方式是：

将各抑斑算法对应的抑斑合成孔径雷达SAR图像局部EKI估计值分别利用直方图统计各像素局部EKI在不同取值区域的对比分布。

本发明的有益效果在于，

本发明与现有方法比较具有以下优点：

1）采用各向异性高斯核平行窗代替传统矩形平行窗生成基于比率边缘强度映射（edge strength map，ESM），从而获得更有利于合成孔径雷达SAR图像区域划分的比率边缘强度映射ESM与方向信息，合成孔径雷达SAR图像区域划分的更加精确，使边缘强度映射ESM与EKI参数的评估更准确。

2）针对比率ESM阈值与阈值化ESM图像中连通像素点集合个数之间的对应关系，提出了一种获得ESM阈值的非监督估计方法，使本发明可无需人为参与设置阈值参数，利用ESM阈值化处理完成对合成孔径雷达SAR图像边缘区域与非边缘区域的有效自适应划分。

3）在合成孔径雷达SAR图像区域化分基础上，对边缘区域采用非极小值抑制方法获得SAR图像较为精确的细化边缘信息，从而自适应的完成了对估计区域更加精细有效的划分与选择。

4）利用ESM约束对所有非边缘区域像素施加大尺度不规则窗完成局部ENL估计，利用ESM方向信息计算合成孔径雷达SAR图像抑斑前后对应边缘像素的梯度比值获得局部EKI估计，获得比传统矩形窗统计更加优异的局部ENL与EKI估计性能。

5）利用二维直方图统计法将合成孔径雷达SAR图像各像素局部ENL与EKI估计形成对比曲线，可更全面直观的展示不同抑斑图像的相干斑抑制性能。

附图说明

图1是本发明方法的流程图；

图2是传统ESM估计方法中所采用的一种带有θ方向的矩形平行窗的结构图；

图3是本发明ESM估计方法中所采用的各向异性高斯窗中方向为0°的各向异性高斯平行窗的结构图；

图4是本发明ESM估计方法中所采用的各向异性高斯窗中方向为45°的各向异性高斯平行窗的结构图；

图5是本发明ESM估计方法中所采用的各向异性高斯窗中方向为90°的各向异性高斯平行窗的结构图；

图6是本发明ESM估计方法中所采用的各向异性高斯窗中方向为135°的各向异性高斯平行窗的结构图；

图7是以图9中ESM图像为基础获得的阈值化ESM图像中连通像素点集合个数随阈值T变化关系曲线图；

图8是一幅测试用256×256大小的幅度格式六视仿真SAR图像；

图9是对图8中的六视仿真SAR图像采用传统ESM估计方法中的四方向矩形平行窗所获得的ESM图像；

图10是对图8中的六视仿真SAR图像采用本发明ESM估计方法中的四方向各向异性高斯平行窗所获得的ESM图像；

图11是对图9中的ESM图像进行阈值化操作获得的阈值化ESM图像；

图12是对图10中的ESM图像进行阈值化操作获得的阈值化ESM图像；

图13是对图8中的六视仿真SAR图像利用经典的空域Kuan滤波抑制相干斑后的抑斑SAR图像；

图14是对图8中的六视仿真SAR图像利用各向异性扩散滤波抑制相干斑后的抑斑SAR图像；

图15是对图13与图14中两类抑斑SAR图像所对应的全部非边缘像素的局部ENL估计值的比值曲线图；

图16是对图13与图14中两类抑斑SAR图像所对应的绝大部分边缘像素的局部EKI估计值的对比直方图；

图17是一幅测试用400×400大小的幅度格式五视城镇郊区真实SAR图像；

图18是对图17中的五视城镇郊区真实SAR图像利用经典的空域Kuan滤波抑制相干斑后的抑斑SAR图像；

图19是对图17中的五视城镇郊区真实SAR图像利用非下采样Contourlet变换域滤波（简记为：NSCT滤波）抑制相干斑后的抑斑SAR图像；

图20是对图18与图19中两类抑斑SAR图像所对应的全部非边缘像素的局部ENL估计值的比值曲线图；

图21是对图18与图19中两类抑斑SAR图像所对应的绝大部分边缘像素的局部EKI估计值的对比直方图；

图中，1.方向为0°的各向异性高斯平行窗，2.方向为45°的各向异性高斯平行窗，3.方向为90°的各向异性高斯平行窗，4.方向为135°的各向异性高斯平行窗，5.仿真SAR图像中用于局部ENL参数估计的同质区一，6.仿真SAR图像中用于局部ENL参数估计的同质区二，7.仿真SAR图像中用于局部ENL参数估计的同质区三，8.仿真SAR图像中用于局部EKI参数估计的边缘区一，9.仿真SAR图像中用于局部EKI参数估计的边缘区二，10.仿真SAR图像中用于局部EKI参数估计的边缘区三，11.城镇郊区真实SAR图像中用于局部ENL参数估计的同质区一，12.城镇郊区真实SAR图像中用于局部ENL参数估计的同质区二，13.城镇郊区真实SAR图像中用于局部ENL参数估计的同质区三，14.城镇郊区真实SAR图像中用于局部EKI参数估计的边缘区一，15.城镇郊区真实SAR图像中用于局部EKI参数估计的边缘区二，16.城镇郊区真实SAR图像中用于局部EKI参数估计的边缘区三。

具体实施方式

下面结合具体实施方式和附图对本发明进行详细说明。

本发明合成孔径雷达图像相干斑噪声抑制性能的非监督评估方法，具体按照以下步骤实施：

步骤1、利用各项异性高斯核平行窗通过比率运算生成合成孔径雷达SAR图像的边缘强度映射dge strength map简称ESM与方向信息：

根据合成孔径雷达SAR图像相干斑的乘性噪声特性，利用图像处理中常用的“梯度计算方法”，这里是指将平行窗两个窗内像素的平均值进行减运算，对合成孔径雷达SAR图像进行区域划分、边缘检测时存在对噪声敏感边缘检测非恒虚警的问题；为此，在合成孔径雷达SAR图像中常利用带方向平行窗通过比率运算获取边缘强度映射ESM，实现对噪声具有较强鲁棒性的边缘恒虚警检测。

比率边缘强度映射ESM的计算需先在如图2所示的矩形平行窗内进行：

计算中心像素C(x,y)在两个平行窗内的局部均值m₁(x,y|θ)与m₂(x,y|θ)，然后通过比较全部θ方向局部均值估计的比值获得ESM估计值与方向信息θ，计算公式为：

$\left\{\begin{array}{c}{V}_{\mathrm{ESM}}(x,y)=\underset{\mathrm{\θ}\∈[0,\mathrm{\π})}{\min }\left[\min (\frac{m_1(x,y|\mathrm{\θ})}{m_2(x,y|\mathrm{\θ})},\frac{m_2(x,y|\mathrm{\θ})}{m_1(x,y|\mathrm{\θ})})\right]\\ \mathrm{\θ}(x,y)=\underset{\mathrm{\θ}}{\arg }\left\{\underset{\mathrm{\θ}\∈[0,\mathrm{\π})}{\min }(\frac{m_1(x,y|\mathrm{\θ})}{m_2(x,y|\mathrm{\θ})},\frac{m_2(x,y|\mathrm{\θ})}{m_1(x,y|\mathrm{\θ})})\right\}\end{array}\right.$

传统比率ESM估计时所采用的图2所示的矩形平行窗，在进行窗内像素平均时，仅采用了简单的算术平均。

本发明方法采用如图3~图6所示的各向异性高斯核平行窗，各向异性高斯核平行窗的方向分别：方向为θ=0°的各向异性高斯平行窗，方向为θ=45°的各向异性高斯平行窗，方向为θ=90°的各向异性高斯平行窗，方向为θ=135°的各向异性高斯平行窗；

带有θ方向的各向异性高斯核平行窗，计算公式如下：

其中：r_x＝xcosθ-ysinθ，r_y＝xsinθ+ycosθ，x与y分别表示水平与垂直方向坐标，θ表示方向，k表示两个平行窗的间距，σ_x与σ_y分别表示x与y方向的标准差，r_x,r_y表示对x，y按照θ方向旋转后形成的新坐标；

利用公式（1）生成的各向异性高斯平行窗权重系数，再将各向异性高斯核平行窗内各像素与待估计中心像素的距离作为权重分配的依据，实现加权平均，即可获得合成孔径雷达SAR图像的边缘强度映射ESM与方向信息，获得更好的SAR图像ESM估计与区域划分；加权平均其实就是权重系数数组，该系数数组中的各元素随着与窗中心像素距离的不同而大小各异。

步骤2，获得合成孔径雷达SAR图像的边缘强度映射ESM阈值T：

在已有基于比率边缘强度映射ESM的合成孔径雷达SAR图像区域划分或边缘检测算法中，阈值T是通过人为手动进行设置，本发明方法利用ESM阈值与阈值化ESM图像中连通像素点集合个数之间的对应关系，给出了一种ESM阈值的非监督估计方法：

1）经步骤1获得了合成孔径雷达SAR图像的边缘强度映射ESM与方向信息后，开始对边缘强度映射ESM图像进行多次重复阈值递减的阈值化处理，即阈值T从最大值1开始递减，每次减少量为0.01，直到阈值T为0.01为止；

2）每次阈值化处理后，统计阈值化处理后的边缘强度映射ESM图像中连通像素点集合的个数，可以获得不同阈值T下阈值化边缘强度映射ESM图像中连通像素点集合个数曲线，其曲线如图7所示，当边缘强度映射ESM图像中连通像素点集合个数曲线由最大值迅速降低，直到出现第一次极小值时，此时具有最好的图像区域划分性能，为此，可以利用阈值化处理边缘强度映射ESM图像中连通像素点集合个数曲线，方便有效的找到合成孔径雷达SAR图像区域划分的阈值T；

如图7所示，其中阈值T=0.84时，该阈值化边缘强度映射ESM图像虚假边缘与边缘断裂现象都最少，此时具有最好的图像区域划分性能。

步骤3、根据ESM阈值实现合成孔径雷达SAR图像边缘与非边缘区域的自动划分：

经步骤2可获得合成孔径雷达SAR图像的边缘强度映射ESM阈值T，实现对图像进行区域划分，阈值化比率边缘强度映射ESM操作是划分合成孔径雷达SAR图像边缘与非边缘区域的一种有效方法，边缘强度映射ESM阈值化可表示为如下公式：

其中，T为阈值，ψ为非零的区域为即为边缘区域，若为零的区域就视为为同质区；

利用步骤2中自适应获得的边缘强度映射ESM阈值T，将阈值T代入公式（2）中，即可获得合成孔径雷达SAR图像边缘区域与同质区域的自动划分。

步骤4、利用非极小值抑制方法从合成孔径雷达SAR图像边缘区域提取细化边缘：

由于局部边缘保持指数EKI估计时需要提供边缘信息，因此需要对合成孔径雷达SAR图像区域划分得到的边缘区域进行细化处理，进而从边缘区域中抽取精确的边缘信息；传统方法借助Canny边缘检测算子，利用非极大值抑制方法对边缘区域进行细化处理，同时考虑到ESM估计值，其取值越小说明其为边缘的概率越大。

本发明方法采用非极小值抑制方法来提取细化边缘，具体处理方法为：

1）若ψ(x,y)＞0，则在ψ域中以位置(x,y)为中心的3×3窗内，沿与θ(x,y)垂直的方向取两相邻值ψ(x₁,y₁)与ψ(x₂,y₂)；

2）若设ψ(x,y)＜ψ(x₂,y₂)与ψ(x,y)＜ψ(x₁,y₁)同时成立为条件Ω，则边缘抽取可表示为

其中，Edge(x,y)＝1的像素为边缘区域中的边缘像素，否则为边缘区域的非边缘像素。

步骤5、获得等效视数ENL估计和局部边缘保持指数EKI估计：

由边缘强度映射ESM约束生成的大尺度不规则窗获得非边缘区域各像素局部等效视数ENL估计，并由抑斑前、后合成孔径雷达SAR图像各边缘像素沿垂直于自身边缘强度映射ESM方向计算梯度模值的比值获得局部EKI估计，具体方法如下：

经步骤4在合成孔径雷达SAR图像区域划分基础上，非边缘区域内各像素局部ENL估计需通过加窗实现，若对像素C(m,n)施加大小为N的窗R，则像素C(m,n)的局部ENL估计可采用如下公式计算，由于窗R相对于整个图像是局部窗，因此所计算的是窗R中各元素的平均值，称为窗内局部均值E，具体的计算公式如下：

$E(m,n)=\frac{1}{N}\underset{(i,j)\∈R}{\mathrm{\Σ}}I(m+i,n+j),---\left(4\right)$

窗中各元素的方差称为窗内局部方差D可表示为：

$E(m,n)=\frac{1}{N-1}\underset{(i,j)\∈R}{\mathrm{\Σ}}{[I(m+i,n+j)-E(m,n)]}^2,---\left(5\right)$

由于ENL的估计需在合成孔径雷达SAR图像同质区进行，若窗R为传统矩形窗，则当窗R内中心像素C靠近图像边缘区域时，窗R会包含边缘区域像素并覆盖不同灰度值水平的两块或多块同质区，造成局部统计量E与D出现较大偏差，使局部ENL估计失效，因此对合成孔径雷达SAR图像非边缘区域各像素加窗时，窗内像素应与窗中心像素C处于一个同质区内，显然传统矩形窗不能实现这一要求。另外，在进行局部ENL估计时，窗R尺度越大将会有更多像素参与到局部ENL估计中，使得局部ENL估计更加平稳有效；但另一方面，直接对各像素施加大尺度窗又使得窗R内容易包含不同区域的像素，造成局部ENL估计偏差较大；为此，对各像素施加ESM约束下的大尺度不规则窗，从窗内参与局部ENL估计运算的各像素中，剔除灰度值起伏较大的像素，避免这些像素参与窗内统计运算，既可以增加窗内参与局部ENL估计的像素个数，同时又避免了不同区域像素的引入而降低局部统计的偏差。

在具体实现时，本发明方法利用步骤4中由阈值化边缘强度映射ESM生成的区域划分信息ψ，对传统矩形窗内各像素施加约束，将传统矩形窗中与窗中心像素C(m,n)不同质的像素剔除，形成不规则窗，为避免直接利用不规则窗而给后续运算带来不便，又通过给传统矩形窗内各元素加权的方法间接实现不规则加窗运算；对于传统矩形窗R，设窗中心像素C(m,n)属于某同质区，则像素I(i,j)∈R对应的权值δ(i,j)可表示为：

将公式（4）调整为：

$E(m,n)=\frac{\underset{(i,j)\∈R}{\mathrm{\Σ}}I(m+i,n+j)\mathrm{\δ}(m+i,n+j)}{\underset{(i,j)\∈R}{\mathrm{\Σ}}\mathrm{\δ}(i,j)},---\left(7\right)$

将公式（5）调整为：

$D(m,n)=\frac{\underset{(i,j)\∈R}{\mathrm{\Σ}}{[I(m+i,n+j)-E(m,n)]}^2\mathrm{\δ}(m+i,n+j)}{\underset{(i,j)\∈R}{\mathrm{\Σ}}\mathrm{\δ}(i,j)-1},---\left(8\right)$

在利用非极小值抑制方法获得合成孔径雷达SAR图像细化边缘基础上，对坐标为(x,y)的任意一个边缘像素估计其局部EKI时，采

用如下的公式估计：

$V_{\mathrm{EKI}}(x,y)=G_{{\mathrm{\θ}}^{\⊥}}^{\′}(x,y)/G_{{\mathrm{\θ}}^{\⊥}}(x,y)---\left(9\right)$

其中，和分别表示抑斑前、抑斑后合成孔径雷达SAR图像位于(x,y)处的边缘像素沿与该像素方向θ(x,y)垂直方向的梯度模值。

步骤6、利用统计方法将各像素局部ENL与EKI估计形成对比曲线，完成SAR图像基于ENL与EKI参数的相干斑抑制性能的评估：

经步骤5获得的大量局部ENL与EKI估计值，利用统计方法获得不同抑斑算法下抑斑SAR图像ENL与EKI参数的对比分布，进而对这两类指标进行全面且直观的总体评估：

对ENL参数的统计方式是：

1）将各抑斑算法对应的抑斑合成孔径雷达SAR图像局部ENL，在分贝单位下，先进行比率运算，具体是利用公式（1）生成的各向异性高斯平行窗权重系数，计算两个平行窗内的局部均值，然后将这两个局部均值直接进行比率运算；

2）再按各像素在合成孔径雷达SAR图像中的实际位置顺序，统计其对比分布。

EKI参数的统计方式是：

将各抑斑算法对应的抑斑合成孔径雷达SAR图像局部EKI估计值分别利用直方图统计各像素局部EKI在不同取值区域的对比分布。

如图8所示，是六视仿真合成孔径雷达SAR图像，如图9及图10所示，是对图8所示的六视仿真合成孔径雷达SAR图像在相同窗尺度条件下，分别采用传统的矩形平行窗与各向异性高斯核平行窗所获得的ESM图，观察图9及图10可看出：图10的ESM图像中的边缘区域要比图9的ESM图像中的边缘区域更清晰，也更光滑；而图10所示同质区由相干斑引起的条纹也比图9所示同质区由相干斑引起的条纹起伏更小；对图9及图10采用相同的阈值化处理获得图11及图12，由图11及图12中很明显可以看出：图12获得的边缘区域远比图11获得的边缘区域光滑且边缘区域断裂现象更少，这说明：各向异性高斯核平行窗可以获得比传统矩形平行窗更好的SAR图像ESM估计与区域划分性能。

为验证本发明方法在利用ENL与EKI两类指标评价合成孔径雷达SAR图像抑斑性能时，有效克服了传统监督评估采用人为选择评估区域所存在的评估不全面且受区域选择影响大的缺陷，设计了两个实验：

实验一：选择一幅256×256幅度格式的六视仿真SAR图像，如图8所示，然后利用经典的空域Kuan滤波与各向异性扩散滤波（简记为PM滤波），分别对图8所示SAR图像进行相干斑抑制，获得如图13与图14所示的抑斑图像，然后再利用本发明方法与传统方法分别对两幅抑斑图像进行相干斑抑制性能评估：其中，图15与图16为本发明方法获得的评估结果，而表1为传统方法在选定区域获得的评估结果。

表1传统方法对仿真SAR图像抑制性能的实验评估结果

在实验一中，从表1所示的由传统方法获得的评估结果可以发现：在图8中的区域5～区域7，Kuan滤波抑斑图像的ENL指标均优于各向异性扩散滤波抑斑图像；而在图8中的区域8～区域10，各向异性扩散滤波抑斑图像的EKI指标在两个区域优于Kuan滤波抑斑图像，因此，传统方法由这些结果可以得出：Kuan滤波抑斑图像在表征相干斑抑制程度的ENL指标上优于各向异性扩散滤波抑斑图像，但在表征SAR图像边缘保持程度的EKI指标上却不如各向异性扩散滤波抑斑图像。另一方面，从图15与图16所示的本发明方法获得的评估结果可以发现：图15所示的两抑斑图像非边缘像素局部ENL比值曲线取值大于0dB部分明显要远远多于小于0dB部分，这充分说明Kuan滤波抑斑图像虽在部分区域相干斑抑制程度优于各向异性扩散滤波抑斑图像，但就整体而言各向异性扩散滤波抑斑图像相干斑抑制程度却优于Kuan滤波抑斑图；在图16所示的两抑斑图像边缘像素局部EKI直方图对比曲线中，Kuan滤波抑斑图像边缘像素局部EKI估计值明显要比各向异性扩散滤波抑斑图像更加集中分布于取值1的周围，而且其局部EKI估计值大于10以上的边缘保持较差像素个数也要远远少于各向异性扩散滤波抑斑图像，这充分说明Kuan滤波抑斑图像的边缘保持程度不仅在整体上好于各向异性扩散滤波抑斑图像，而且边缘像素保持较差的像素个数也明显少于各向异性扩散滤波抑斑图像。显然，对同样的两幅抑斑图像，传统方法与本文发明方法得到了刚好相反的相干斑抑制性能评价结果。

实验二：选择一幅400×400幅度格式的五视城镇郊区SAR图像，如图17所示，然后利用空域Kuan滤波与非下采样Contourlet变换域滤波（简记为NSCT滤波），分别对图17所示SAR图像进行相干斑抑制，获得如图18与图19所示的抑斑图像，然后再利用本发明方法与传统方法分别对两幅抑斑图像进行相干斑抑制性能评估：其中，图20与图21为本文方法获得的评估结果，而表2为传统方法在选定区域获得的评估结果。

表2传统方法对真实SAR图像抑制性能的实验评估结果

在真实城镇郊区SAR图像的评估实验二中，表2所示的由传统方法获得的评估结果与图20以及图21所示的本发明方法获得的评估结果也有类似的矛盾结论。

在上述两个实验中，导致传统方法与本发明方法得出矛盾评估结果，显然是由传统方法采用监督评估方式由人为主观的选择SAR图像部分区域参与评估造成的。因为依靠人为参与选择大量的评估区域来获取充分的局部统计样本麻烦且不现实，而利用有限的ENL与EKI局部统计样本来对SAR图像整体抑斑性能进行评估时，就可能造成评估结果因区域选择不同而大相径庭。另外，由于本文方法采用非监督方式，自适应的选择SAR图像中的非边缘像素与边缘像素，分别用于估计局部ENL与EKI参数，不仅排除了人为主观因素对评估结果的影响，而且可以获取充分的局部统计样本，为获得更全面的评估结果奠定基础。

Claims (1)

1.合成孔径雷达图像相干斑噪声抑制性能的非监督评估方法，其特征在于，具体按照以下步骤实施：

步骤1、利用各项异性高斯核平行窗通过比率运算生成合成孔径雷达SAR图像的边缘强度映射edge strength map简称ESM与方向信息：

采用带有θ方向的各向异性高斯核平行窗按照公式(1)计算各向异性高斯平行窗权重系数，然后生成合成孔径雷达SAR图像的边缘强度映射ESM与方向信息：

带有θ方向的各向异性高斯核平行窗计算公式如下：

其中：r_x＝xcosθ-ysinθ，r_y＝xsinθ+ycosθ，x与y分别表示水平与垂直方向坐标，θ表示方向，k表示两个平行窗的间距，σ_x与σ_y分别表示x与y方向的标准差，r_x,r_y表示对x，y按照θ方向旋转后形成的新坐标；

步骤1具体按照以下步骤实施：

采用带有θ方向的各向异性高斯核平行窗，各向异性高斯核平行窗的方向分别：方向为θ＝0°的各向异性高斯平行窗，方向为θ＝45°的各向异性高斯平行窗，方向为θ＝90°的各向异性高斯平行窗，方向为θ＝135°的各向异性高斯平行窗；

利用公式(1)生成的各向异性高斯平行窗权重系数，再将各向异性高斯核平行窗内各像素与待估计中心像素的距离作为权重分配的依据，实现加权平均，即可获得合成孔径雷达SAR图像的边缘强度映射ESM与方向信息，获得更好的SAR图像ESM估计与区域划分；

步骤2、获得合成孔径雷达SAR图像的边缘强度映射ESM阈值T，具体按照以下步骤实施：

步骤2.1、经步骤1获得合成孔径雷达SAR图像的边缘强度映射ESM与方向信息之后，开始对边缘强度映射ESM图像进行多次重复阈值递减的阈值化处理，即阈值T从最大值1开始递减，每次减少量为0.01，直至阈值T为0.01为止；

步骤2.2、每次阈值化处理后，统计阈值化处理后的边缘强度映射ESM图像中连通像素点集合的个数，即可获得不同阈值T下阈值化边缘强度映射ESM图像中连通像素点集合个数曲线，当边缘强度映射ESM图像中连通像素点集合个数曲线由最大值迅速降低直至出现第一次极小值时，此时具有最好的图像区域划分性能，这个极小值就是合成孔径雷达SAR图像区域划分的阈值T；

步骤3、根据ESM阈值实现合成孔径雷达SAR图像边缘与非边缘区域的自动划分：

将步骤2中获得的边缘强度映射ESM阈值T代入公式(2)中，即可获得合成孔径雷达SAR图像边缘区域与同质区域的自动划分；

其中，ψ若为非零的区域为即为边缘区域，ψ若为零的区域就视为为同质区；

步骤4、利用非极小值抑制方法从合成孔径雷达SAR图像边缘区域提取细化边缘，具体按照以下步骤实施：

利用非极小值抑制方法从合成孔径雷达SAR图像边缘区域提取细化边缘：

步骤4.1、若ψ(x,y)>0，则在ψ域中以位置(x,y)为中心的3×3窗内，沿与θ(x,y)垂直的方向取两相邻值ψ(x₁,y₁)与ψ(x₂,y₂)；

步骤4.2、若设ψ(x,y)<ψ(x₂,y₂)与ψ(x,y)<ψ(x₁,y₁)同时成立为条件Ω，则边缘抽取可表示为：

其中，Edge(x,y)＝1的像素为边缘区域中的边缘像素，否则为边缘区域的非边缘像素；

步骤5、获得等效视数ENL估计和局部边缘保持指数EKI估计，具体方法如下：

经步骤4在合成孔径雷达SAR图像区域划分基础上，非边缘区域内各像素局部ENL估计要通过加窗实现，对像素C(m,n)施加大小为N的窗R，像素C(m,n)的局部ENL估计采用公式(4)计算，窗R相对于整个图像是局部窗，计算出的是窗R中各元素的平均值，称为窗内局部均值E，具体的计算公式如下：

$E(m,n)=\frac{1}{N}\underset{(i,j)\&Element;R}{\mathrm{\&Sigma;}}I(m+i,n+j)---\left(4\right);$

窗中各元素的方差称为窗内局部方差D可表示为：

$D(m,n)=\frac{1}{N-1}\underset{(i,j)\&Element;R}{\mathrm{\&Sigma;}}{[I(m+i,n+j)-E(m,n)]}^2---\left(5\right);$

经步骤4中由阈值化边缘强度映射ESM生成的区域划分信息ψ，对矩形窗R内各像素施加约束，将矩形窗R中与窗中心像素C(m,n)不同质的像素剔除，形成不规则窗；再通过给矩形窗R内各元素加权的方法间接实现不规则加窗运算；对于矩形窗R，窗中心像素C(m,n)属于某同质区，则像素I(i,j)∈R对应的权值δ(i,j)可表示为：

可将公式(4)调整为：

$E(m,n)=\frac{\underset{(i,j)\&Element;R}{\mathrm{\&Sigma;}}I(m+i,n+j)\mathrm{\&delta;}(m+i,n+j)}{\underset{(i,j)\&Element;R}{\mathrm{\&Sigma;}}\mathrm{\&delta;}(i,j)}---\left(7\right);$

将公式(5)调整为：

$D(m,n)=\frac{\underset{(i,j)\&Element;R}{\mathrm{\&Sigma;}}{[I(m+i,n+j)-E(m,n)]}^2\mathrm{\&delta;}(m+i,n+j)}{\underset{(i,j)\&Element;R}{\mathrm{\&Sigma;}}\mathrm{\&delta;}(i,j)-1}---\left(8\right);$

利用非极小值抑制方法获得合成孔径雷达SAR图像细化边缘基础上，对坐标为(x,y)的任意一个边缘像素估计其局部EKI时，采用如下的公式(9)估计：

$V_{\mathrm{EKI}}(x,y)=G_{{\mathrm{\&theta;}}^{\&perp;}}^{\&prime;}(x,y)/G_{{\mathrm{\&theta;}}^{\&perp;}}(x,y)---\left(9\right);$

其中，和分别表示抑斑前、抑斑后合成孔径雷达SAR图像位于(x,y)处的边缘像素沿与该像素方向θ(x,y)垂直方向的梯度模值；

步骤6、利用统计方法将各像素局部ENL与EKI估计形成对比曲线，完成SAR图像基于ENL与EKI参数的相干斑抑制性能的评估，具体按照以下步骤实施：

经步骤5获得的大量局部ENL与EKI估计值，利用统计方法获得不同抑斑算法下抑斑SAR图像ENL与EKI参数的对比分布，进而对这两类指标进行全面且直观的总体评估：

ENL参数的统计方式是：

1)将各抑斑算法对应的抑斑合成孔径雷达SAR图像局部ENL，在分贝单位下，先进行比率运算，具体是利用公式(1)生成的各向异性高斯平行窗权重系数，计算两个平行窗内的局部均值，然后将这两个局部均值直接进行比率运算；

2)再按各像素在合成孔径雷达SAR图像中的实际位置顺序，统计其对比分布；

EKI参数的统计方式是：

将各抑斑算法对应的抑斑合成孔径雷达SAR图像局部EKI估计值分别利用直方图统计各像素局部EKI在不同取值区域的对比分布。