CN109544567A - 基于两阶段区域合并策略的极化合成孔径雷达图像分割方法 - Google Patents

Abstract

本发明公开了一种基于两阶段区域合并策略的极化合成孔径雷达图像分割方法，目的是提升非均匀区域的分割精度，避免过多的复杂运算，提高计算效率。技术方案是首先对PolSAR图像进行超像素分割，然后依次进行Wishart合并和KummerU合并。Wishart合并阶段对明确来自相同地物的原始超像素进行快速合并，而KummerU合并阶段对剩余的区域进行迭代合并，得到最终的分割结果。采用本发明可提升合并的准确性和计算效率，且将KummerU能量损失、边缘惩罚函数值以及匀质度惩罚因子相结合，得到新的度量标准，可以进一步提升分割精度。

Description

基于两阶段区域合并策略的极化合成孔径雷达图像分割方法

技术领域

本发明属于雷达图像处理领域，涉及极化合成孔径雷达(PolarimetricSynthetic Aperture Radar，简称PolSAR)图像分割方法，尤其是使用区域合并框架的PolSAR图像分割方法。

背景技术

PolSAR是一种先进的遥感观测***，能够获取比单极化SAR更加丰富、全面的地物和目标散射信息，在国民经济建设和国防军事领域取得了非常广泛的应用。在PolSAR图像解译中，很多情况下需要首先将图像分割为局部同质区域，然后基于这些区域(也称为对象)进行处理，这对于提升后续处理效率、降低相干斑噪声影响以及获得均匀稳定的结果具有重要意义。理想的分割结果需要与图中地物和目标的轮廓保持一致，同时各分割区域的大小和形状应该随不同地物的大小和形状而变化，即分割结果既能使用单个分割区域来描述大块目标，也能保留细小目标的信息。单尺度分割往往会导致过分割或欠分割的现象，过分割即分割尺度较小，大块地物类型或目标也被分割为小的区域；而欠分割指分割尺度较大，不能有效地将相邻的地物目标分割开来，导致很多细节信息的丢失。为了解决该问题，一种比较好的思路是在过分割的基础上，对来自同一地物的相邻小区域进行合并，进而得到最终的分割结果。

在基于区域合并(Region Merging)的PolSAR图像分割方法中，通常依据相邻区域之间的距离度量或者能量损失函数来决定是否对该相邻区域进行合并，现有文献中的距离度量或者能量损失函数大都建立在PolSAR数据服从Wishart分布的基础之上。例如文献1：Beaulieu J M，Touzi R，“Segmentation of textured polarimetric SAR scenes bylikelihood approximation”，IEEE Transactions on Geoscience and Remote Sensing，2004，42(10)：2063-2072(Beaulieu J M等人于2004年在《电气与电子工程师协会地理科学与遥感会刊》第42卷第10期发表的“通过似然近似方法对极化SAR纹理区域进行分割”)基于Wishart能量损失函数来确定逐次被合并的相邻区域。Wishart分布是PolSAR图像处理中应用最广泛的统计分布，基于Wishart分布推导出的Wishart能量损失函数具有形式简单、运算高效等优点，但是Wishart分布比较适合描述均匀区域的PolSAR数据，随着图像分辨率的提高，PolSAR数据包含更加丰富的细节信息，地物后向散射特性的纹理引起PolSAR数据的非高斯特性，导致数据统计分布偏离Wishart分布。因此，文献1中方法在非均匀区域(森林、城区等)不能有效保持区域的细节信息。为了改善在非均匀区域的分割效果，需要结合PolSAR数据在复杂纹理区域的统计特性来构建区域之间新的距离度量或者能量损失函数。一些方法利用PolSAR数据的K分布、KummerU分布等来改善图像分割性能，例如文献2：Bombrun L，Beaulieu J M，“Fisher distribution for texture modeling ofPolarimetric SAR data”，IEEE Geoscience and Remote Sensing Letters，2008，5(3)：512-516(Bombrun L等人于2008年在《电气与电子工程师协会地理科学与遥感快报》第5卷第3期发表的“将Fisher分布用于PolSAR数据的纹理建模”)提出PolSAR数据的KummerU分布，并推导出KummerU能量损失函数，用于指导相邻区域之间的合并。KummerU分布能够更好地区分具有不同纹理的区域，但是KummerU统计分布的表达式以及参数估计非常复杂，导致整个分割处理的效率很低。因此，在实施区域合并过程中，需要同时兼顾分割精度和运算效率，从而提升分割方法的可靠性和实用性。

另外，区域合并依赖于图像过分割的结果。近些年来，超像素分割在图像处理领域引起广泛关注，是一种非常流行的图像过分割方法。超像素分割指把图像分割为尺寸相似、形状规则的匀质区域，并且这些区域能保留图中目标的信息，同时与目标边界相吻合。因此，所得到的超像素可以作为后续合并操作的基本处理单元。在PolSAR图像超像素分割方面，研究人员提出了多种方法，例如文献3：Wang W，Xiang D，et al，“SuperpixelSegmentation of Polarimetric SAR Data Based on Integrated Distance Measureand Entropy Rate Method”，IEEE Journal of Selected Topics in Applied EarthObservations and Remote Sensing，2017，10(9)，4045-4058(Wang W等人于2017年在《电气与电子工程师协会在应用地球观测与遥感领域主题期刊》第10卷第9期发表的“基于统一距离度量和熵率方法的PolSAR数据超像素分割方法”)提出了一种新的PolSAR图像超像素分割方法，在均匀区域和非均匀区域都可以取得较为理想的超像素结果。

发明内容

本发明要解决的技术问题是提出一种基于两阶段区域合并策略的极化合成孔径雷达图像分割方法，使用KummerU分布来提升非均匀区域的分割精度，同时避免过多的复杂运算，从而提高计算效率。

技术方案是首先对PolSAR图像进行超像素分割，然后依次进行Wishart合并(对应Wishart合并阶段)和KummerU合并(对应KummerU合并阶段)。Wishart合并阶段(WishartMerging Stage，简称WMS)对明确来自相同地物的原始超像素进行快速合并，而KummerU合并阶段(KummerU Merging Stage，简称KUMS)对剩余的区域进行迭代合并，得到最终的分割结果。

本发明主要包括以下步骤：

步骤1：对PolSAR图像进行超像素分割，方法是：

采用文献3中所述PolSAR图像超像素分割方法对原始PolSAR图像进行超像素分割，得到超像素集合V＝{v₁,v₂,...,v_i,...,v_j,...,v_N}，其中N为超像素个数，v_i表示第i个超像素，1≤i≤N，1≤j≤N。超像素集合V及集合中元素的关系可以用区域邻接图(RegionAdjacency Graph，简称RAG)G＝(V,A)来表示，其中A为超像素之间相邻关系的集合，a_ij∈A表示超像素v_i和v_j在几何位置上相邻，令N_a表示A中元素的个数，N_a为正整数。另外，用b_ij来表示相邻超像素v_i和v_j之间的共享边界，为该边界上像素的集合。

步骤2：对V进行Wishart合并，方法是：

2.1对V中所有超像素v₁,v₂,...,v_i,...,v_j,...,v_N计算分布参数Z₁,Z₂,...,Z_i,...,Z_j,...,Z_N，分布参数为某区域(此处是超像素，即合并的最小元素；后面步骤中所述区域指超像素合并后的结果)中协方差矩阵C的均值，v_i的分布参数为Z_i，按公式(1)计算：

其中，|v_i|_c表示超像素v_i的基数，即v_i中像素的个数，下标c指代基数。s用来表示一个像素，C_s为像素s的协方差矩阵。

2.2根据A计算所有相邻超像素之间的Wishart度量标准，构建第一三元组集合D，具体步骤如下：

2.2.1令k＝1，集合D为空集；

2.2.2设A中第k个元素为a_ij，1≤i≤N，1≤j≤N，a_ij所对应的两个相邻超像素为v_i和v_j；

2.2.3计算超像素v_i和v_j合并前后产生的Wishart能量损失计算式为：

其中，v_ij＝v_i∪v_j表示超像素v_i和v_j合并之后的区域，|v_ij|_c为v_ij的基数，|v_i|_c和|v_j|_c分别为v_i和v_j的基数，Z_ij表示v_ij的分布参数，由公式(3)计算得出。|Z_ij|表示分布参数Z_ij的行列式值，|Z_i|和|Z_j|分别为分布参数Z_i和Z_j的行列式值；

2.2.4计算超像素v_i和v_j之间的边缘惩罚函数值计算式为：

其中，V_s表示像素s的边缘强度值，该值介于0到1之间，K为决定边缘惩罚强度的参数，取值为0.3；

2.2.5将2.2.3步计算得到的Wishart能量损失和2.2.4步得到的边缘惩罚函数值相加，得到超像素v_i和v_j在Wishart合并阶段的度量标准按公式(5)计算：

其中，参数β是控制相对于的权重，取值设为5；

2.2.6令度量标准元素将d_k以及所对应的i和j构成三元组(d_k,i,j)，放入集合D中；

2.2.7令k＝k+1，判断k是否小于等于N_a，N_a为A中元素的个数，若满足，转2.2.2步；否则，得到了所有相邻超像素在Wishart合并阶段的度量标准以及在此基础上构成的三元组的集合D，D中元素个数为N_a，转2.3步。

2.3对D中三元组按三元组中第一个数值的大小进行排序，排列顺序为从小到大的顺序，取前M个三元组，构成集合D_M，1≤M≤N_a，M的取值根据超像素个数N以及百分比参数a来设定，M＝N×a％，a取值范围为40～60。

2.4根据D_M中三元组，在A中找出所对应的M个元素，构成集合A_M，具体步骤如下：

2.4.1令k₁＝1，A_M为空集；

2.4.2取出D_M中第k₁个三元组，设该三元组第二和第三个数值为i和j，1≤i≤N，1≤j≤N，将a_ij放入集合A_M；

2.4.3令k₁＝k₁+1，判断k₁是否小于等于M，若满足，转2.4.2步；否则，得到包含M个元素的集合A_M，转2.5步。

2.5基于A_M对超像素集合V进行合并，得到合并后区域的集合V⁽⁰⁾和区域邻接图G⁽⁰⁾＝(V⁽⁰⁾,A⁽⁰⁾)，具体步骤如下：

2.5.1令k₁＝1；

2.5.2设A_M中第k₁个元素为a_ij，1≤i≤N，1≤j≤N，a_ij所对应的两个相邻超像素为v_i和v_j，则将超像素v_i和v_j合并；

2.5.3令k₁＝k₁+1，判断k₁是否小于等于M，若满足，转2.5.2步；否则，即完成了Wishart合并，得到合并后区域的集合V⁽⁰⁾，V⁽⁰⁾中元素为区域，设区域个数为N⁽⁰⁾，区域表示为1≤m≤N⁽⁰⁾，1≤n≤N⁽⁰⁾。V⁽⁰⁾及V⁽⁰⁾中元素的关系用区域邻接图G⁽⁰⁾＝(V⁽⁰⁾,A⁽⁰⁾)来表示，其中A⁽⁰⁾为Wishart合并后区域之间相邻关系的集合，令表示A⁽⁰⁾中元素的个数，为正整数，用b_mn来表示相邻区域r_m和r_n之间的共享边界。

步骤3：对V⁽⁰⁾进行KummerU迭代合并，方法是：

3.1对V⁽⁰⁾中每个区域，进行参数估计，并计算区域的KummerU能量，具体步骤如下：

3.1.1令k₂＝1；

3.1.2V⁽⁰⁾中第k₂个元素为r_k，根据式(6)计算r_k的分布参数Z_k：

其中，|r_k|_c表示区域r_k的基数，即r_k中像素的个数；

3.1.3基于区域r_k中所有像素的协方差矩阵，估计区域r_k中KummerU分布的纹理形状参数ξ_k和ζ_k。参数估计方法采用文献4：Anfinsen S N，Eltoft T，“Application of theMatrix-Variate Mellin Transform to Analysis of Polarimetric Radar Images”，IEEE Transactions on Geoscience and Remote Sensing，2011，49(6)：2281-2295(Anfinsen S N等人于2011年在《电气与电子工程师协会地理科学与遥感会刊》第49卷第6期发表的“矩阵变量Mellin变换在PolSAR图像分析中的应用”)中的矩阵对数累计方法；

3.1.4计算区域r_k的KummerU能量计算式为：

其中，L为图像视数，d等于3，为矩阵的迹，Γ(x)是x的伽马函数值，x为自变量，为以Ld+ζ_k、Ld-ξ_k+1、为参数的第二类合流超几何函数；

3.1.5令k₂＝k₂+1，判断k₂是否小于等于N⁽⁰⁾，若满足，转3.1.2步；否则，得到V⁽⁰⁾中所有区域的KummerU能量，转3.2步。

3.2计算V⁽⁰⁾中所有区域的KummerU能量总和SE⁽⁰⁾：

3.3根据A⁽⁰⁾计算所有相邻区域之间的KummerU度量标准，构建第二三元组集合D⁽⁰⁾，具体步骤如下：

3.3.1令k₃＝1，集合D⁽⁰⁾为空集；

3.3.2设A⁽⁰⁾中第k₃个元素为a_mn，1≤m≤N⁽⁰⁾，1≤n≤N⁽⁰⁾，a_mn所对应的两个相邻区域为r_m和r_n；

3.3.3计算区域r_m和r_n合并前后产生的KummerU能量损失计算式为：

其中，r_mn＝r_m∪r_n表示区域r_m和r_n合并之后的区域，Z_mn为区域r_mn的分布参数，由公式(11)计算得出，ξ_mn和ζ_mn为区域r_mn的纹理形状参数，参照步骤3.1.3基于区域r_mn中所有像素的协方差矩阵估计得到。由式(10)可以看出，该式前三项与式(2)中的Wishart能量损失相似，式(10)中第二行可以看作是由PolSAR数据的纹理模型引入的修正部分。

3.3.4计算区域r_m和r_n之间的边缘惩罚函数值计算式为：

3.3.5分别计算区域r_m、区域r_n以及合并后区域r_mn的匀质度H_m，H_n和H_mn，计算采用变化系数(Coefficient of Variation，CoV)的方法，具体方法参考文献5：Lopes A，Touzi R，Nezry E，“Adaptive speckle filters and scene heterogeneity”IEEE Transactionson Geoscience and Remote Sensing，1990，28(6)：992-1000(Lopes A等人于1990年在《电气与电子工程师协会地理科学与遥感会刊》第28卷第6期发表的“自适应相干斑滤波器和场景异质性”)。在此基础上，计算区域合并的匀质度惩罚因子F_mn，计算式为：

3.3.6将3.3.3步计算得到的KummerU能量损失3.3.4步得到的边缘惩罚函数值以及3.3.5步得到的匀质度惩罚因子F_mn相结合，得到区域r_m和r_n在KummerU合并阶段的度量标准按公式(14)计算：

3.3.7令第二度量标准元素将以及所对应的m和n构成三元组放入集合D⁽⁰⁾中；

3.3.8令k₃＝k₃+1，判断k₃是否小于等于 为A⁽⁰⁾中元素的个数，若满足，转3.3.2步；否则，得到所有相邻区域在KummerU合并阶段的度量标准以及在此基础上构成的三元组的集合D⁽⁰⁾，D⁽⁰⁾中元素个数为转3.4步。

3.4令l表示KummerU迭代合并的次数，且l＝0，表明此时尚未对区域进行合并。

3.5找出D^(l)里所有三元组中第一个数值最小的三元组，设该三元组第二和第三个数值为m和n，则对区域r_m和r_n进行合并，合并后区域仍记为区域r_m(相当于r_n合并进了r_m)。同时，将D^(l)中与区域r_m和r_n相关的三元组从集合中删除，将删除元素之后的D^(l)记为D^(l+1)。

3.6更新3.5步区域合并之后的区域邻接图，用G^(l+1)＝(V^(l+1),A^(l+1))来表示，并更新合并后的能量总和SE^(l+1)和D^(l+1)：

3.6.1区域合并后的能量总和SE^(l+1)更新为：

其中，SE^(l)为区域合并前的KummerU能量总和，为第3.5步合并过程中产生的KummerU能量损失，由3.3.3步计算得出；

3.6.2基于区域r_m中所有像素的协方差矩阵，对r_m的分布参数Z_m以及纹理形状参数ξ_m和ζ_m进行估计，估计方法参照第3.1.2步和第3.1.3步；

3.6.3将A^(l+1)中与区域r_m相关的相邻关系元素构成设集合中元素个数为对中元素逐个计算相应的第二度量标准，并对D^(l+1)进行更新，具体步骤如下：

3.6.3.1令k₄＝1；

3.6.3.2设中第k₄个元素为a_mn，计算KummerU能量损失边缘惩罚函数和匀质度惩罚因子F_mn，并按公式(14)计算得到度量标准计算方法参照第3.3.3步至3.3.6步；

3.6.3.3令第二度量标准元素将以及所对应的m和n构成三元组放入集合D^(l+1)中；

3.6.3.4令k₄＝k₄+1，判断k₄是否小于等于若满足，转3.6.3.2步；否则，实现了对D^(l+1)的更新，转3.6.4步。

3.6.4对KummerU迭代合并次数l进行更新：l＝l+1，判断l是否小于等于N⁽⁰⁾-1，若满足，转3.5步；否则，此时所有区域被合并为一个整体，转步骤4。

在每一步迭代合并过程中，步骤3.6.1得到新的KummerU能量总和SE，从而得到KummerU能量总和SE序列：(SE⁽⁰⁾由3.2步得到，由3.6.1得到)

步骤4：确定PolSAR图像最终分割结果的区域个数

基于SE序列采用文献6：Salvador S，Chan P，“Determining the number of clusters/segments in hierarchical clustering/segmentation algorithms”Proc.16th IEEE ICTAI，2004，578–579(Salvador S等人于2004年在《第16届电气与电子工程师协会人工智能工具会议论文集》发表的“在层次分割方法中确定分割区域的个数”)中的方法来确定最终分割结果所对应的区域个数，得到区域个数为N_r。该方法构建一种评估曲线图，图的x轴为区域个数，y轴表示某一种评估指标(在本发明中为KummerU能量总和SE)，曲线图的拐点，即最大曲率所对应的点则为最终分割结果的区域个数。基于区域个数N_r确定最终的分割结果为

采用本发明可以达到以下技术效果：

1.本发明在第2.2.5步和3.3.6步中使用边缘惩罚函数，当两个相邻超像素或区域之间存在明显边缘时，边缘惩罚函数值增大，将会阻止这两个区域之间的合并，有助于提升合并的准确性。

2.本发明在步骤2(Wishart合并阶段)直接确定在本阶段需要合并的所有超像素对，即集合A_M，不需要对度量标准进行更新，可以显著提升计算效率。

3.本发明在步骤3(KummerU合并阶段)将KummerU能量损失、边缘惩罚函数值以及匀质度惩罚因子相结合，得到新的度量标准，可以进一步提升分割精度。

附图说明

图1本发明总体流程图；

图2PolSAR图像合并过程对应的精度-回召率曲线示意图。

具体实施方式

下面结合附图对本发明进行详细解释和说明。

图1是本发明总体流程图，本发明包括以下步骤：

步骤1：对PolSAR图像进行超像素分割，方法是：

采用PolSAR图像超像素分割方法对原始PolSAR图像进行超像素分割，得到超像素集合V＝{v₁,v₂,...,v_i,...,v_j,...,v_N}，其中N为超像素个数，v_i表示第i个超像素，1≤i≤N，1≤j≤N。用区域邻接图G＝(V,A)表示超像素集合V及集合中元素的关系，其中A为超像素之间相邻关系的集合，a_ij∈A表示超像素v_i和v_j在几何位置上相邻，令N_a表示A中元素的个数，N_a为正整数，用b_ij来表示相邻超像素v_i和v_j之间的共享边界，为该边界上像素的集合。

步骤2：对V进行Wishart合并，方法是：

2.1对V中所有超像素v₁,v₂,...,v_i,...,v_j,...,v_N计算分布参数Z₁,Z₂,...,Z_i,...,Z_j,...,Z_N，按公式(1)计算v_i的分布参数Z_i。

2.2根据A计算所有相邻超像素之间的Wishart度量标准，构建第一三元组集合D，具体步骤如下：

2.2.1令k＝1，集合D为空集；

2.2.2设A中第k个元素为a_ij，1≤i≤N，1≤j≤N，a_ij所对应的两个相邻超像素为v_i和v_j；

2.2.3按公式(2)、(3)计算超像素v_i和v_j合并前后产生的Wishart能量损失

2.2.4按公式(4)计算超像素v_i和v_j之间的边缘惩罚函数值

2.2.5按公式(5)计算超像素v_i和v_j在Wishart合并阶段的度量标准

2.2.6令度量标准元素将d_k以及所对应的i和j构成三元组(d_k,i,j)，放入集合D中；

2.2.7令k＝k+1，判断k是否小于等于N_a，N_a为A中元素的个数，若满足，转2.2.2步；否则，得到了所有相邻超像素在Wishart合并阶段的度量标准以及在此基础上构成的三元组的集合D，D中元素个数为N_a，转2.3步。

2.3对D中三元组按三元组中第一个数值的大小进行排序，排列顺序为从小到大的顺序，取前M个三元组，构成集合D_M，1≤M≤N_a，M的取值根据超像素个数N以及百分比参数a来设定，M＝N×a％，a取值范围为40～60。

2.4根据D_M中三元组，在A中找出所对应的M个元素，构成集合A_M，具体步骤如下：

2.4.1令k₁＝1，A_M为空集；

2.4.2取出D_M中第k₁个三元组，设该三元组第二和第三个数值为i和j，1≤i≤N，1≤j≤N，将a_ij放入集合A_M；

2.4.3令k₁＝k₁+1，判断k₁是否小于等于M，若满足，转2.4.2步；否则，得到包含M个元素的集合A_M，转2.5步。

2.5基于A_M对超像素集合V进行合并，具体步骤如下：

2.5.1令k₁＝1；

2.5.2设A_M中第k₁个元素为a_ij，1≤i≤N，1≤j≤N，a_ij所对应的两个相邻超像素为v_i和v_j，则将超像素v_i和v_j合并；

2.5.3令k₁＝k₁+1，判断k₁是否小于等于M，若满足，转2.5.2步；否则，即完成了Wishart合并，得到合并后区域的集合V⁽⁰⁾，V⁽⁰⁾中元素为区域，设区域个数为N⁽⁰⁾，区域表示为1≤m≤N⁽⁰⁾，1≤n≤N⁽⁰⁾。V⁽⁰⁾及V⁽⁰⁾中元素的关系用区域邻接图G⁽⁰⁾＝(V⁽⁰⁾,A⁽⁰⁾)来表示，其中A⁽⁰⁾为Wishart合并后区域之间相邻关系的集合，令表示A⁽⁰⁾中元素的个数，为正整数。另外，用b_mn来表示相邻区域r_m和r_n之间的共享边界。

步骤3：对V⁽⁰⁾进行KummerU迭代合并，方法是：

3.1对V⁽⁰⁾中每个区域，进行参数估计，并计算区域的KummerU能量，具体步骤如下：

3.1.1令k₂＝1；

3.1.2V⁽⁰⁾中第k₂个元素为r_k，根据式(6)计算r_k的分布参数Z_k；

3.1.3基于区域r_k中所有像素的协方差矩阵，估计区域r_k中KummerU分布的纹理形状参数ξ_k和ζ_k。

3.1.4根据公式(7)、(8)计算区域r_k的KummerU能量

3.1.5令k₂＝k₂+1，判断k₂是否小于等于N⁽⁰⁾，若满足，转3.1.2步；否则，得到V⁽⁰⁾中所有区域的KummerU能量，转3.2步。

3.2按公式(9)计算V⁽⁰⁾中所有区域的KummerU能量总和SE⁽⁰⁾；

3.3根据A⁽⁰⁾计算所有相邻区域之间的KummerU度量标准，构建第二三元组集合D⁽⁰⁾，具体步骤如下：

3.3.1令k₃＝1，集合D⁽⁰⁾为空集；

3.3.2设A⁽⁰⁾中第k₃个元素为a_mn，1≤m≤N⁽⁰⁾，1≤n≤N⁽⁰⁾，a_mn所对应的两个相邻区域为r_m和r_n；

3.3.3根据公式(10)、(11)计算区域r_m和r_n合并前后产生的KummerU能量损失

3.3.4根据公式(12)计算区域r_m和r_n之间的边缘惩罚函数值

3.3.5分别计算区域r_m、区域r_n以及合并后区域r_mn的匀质度H_m，H_n和H_mn并按公式(13)计算区域合并的匀质度惩罚因子F_mn；

3.3.6按公式(14)计算区域r_m和r_n在KummerU合并阶段的度量标准

3.3.7令第二度量标准元素将以及所对应的m和n构成三元组放入集合D⁽⁰⁾中；

3.3.8令k₃＝k₃+1，判断k₃是否小于等于 为A⁽⁰⁾中元素的个数，若满足，转3.3.2步；否则，得到所有相邻区域在KummerU合并阶段的度量标准以及在此基础上构成的三元组的集合D⁽⁰⁾，D⁽⁰⁾中元素个数为转3.4步。

3.4令l表示KummerU迭代合并的次数，且l＝0，表明此时尚未对区域进行合并。

3.5找出D^(l)里所有三元组中第一个数值最小的三元组，设该三元组第二和第三个数值为m和n，则对区域r_m和r_n进行合并，合并后区域仍记为区域r_m(相当于r_n合并进了r_m)。同时，将D^(l)中与区域r_m和r_n相关的三元组从集合中删除，将删除元素之后的D^(l)记为D^(l+1)。

3.6更新3.5步区域合并之后的区域邻接图，用G^(l+1)＝(V^(l+1),A^(l+1))来表示，并更新合并后的能量总和SE^(l+1)和D^(l+1)：

3.6.1按公式(15)更新区域合并后的能量总和SE^(l+1)；

3.6.2基于区域r_m中所有像素的协方差矩阵，对r_m的分布参数Z_m以及纹理形状参数ξ_m和ζ_m进行估计；

3.6.3将A^(l+1)中与区域r_m相关的相邻关系元素构成设集合中元素个数为对中元素逐个计算相应的第二度量标准，并对D^(l+1)进行更新，具体步骤如下：

3.6.3.1令k₄＝1；

3.6.3.2设中第k₄个元素为a_mn，计算KummerU能量损失边缘惩罚函数和匀质度惩罚因子F_mn，并按公式(14)计算得到度量标准

3.6.3.3令第二度量标准元素将以及所对应的m和n构成三元组放入集合D^(l+1)中；

3.6.3.4令k₄＝k₄+1，判断k₄是否小于等于若满足，转3.6.3.2步；否则，实现了对D^(l+1)的更新，转3.6.4步。

3.6.4对KummerU迭代合并次数l进行更新：l＝l+1，判断l是否小于等于N⁽⁰⁾-1，若满足，转3.5步；否则，此时所有区域被合并为一个整体，转步骤4。

步骤4：确定PolSAR图像最终分割结果的区域个数

基于SE序列采用文献6即Salvador S等人于2004年在《第16届电气与电子工程师协会人工智能工具会议论文集》发表的“在层次分割方法中确定分割区域的个数”中的方法来确定理想分割结果所对应的区域个数，得到区域个数为N_r。该方法构建一种评估曲线图，图的x轴为区域个数，y轴表示某一种评估指标(在本发明中为KummerU能量总和SE)，曲线图的拐点，即最大曲率所对应的点则为最终分割结果的区域个数。基于区域个数N_r确定最终的分割结果为

实施例：为说明本发明的有效性，使用一副实测PolSAR图像开展分割实验。同时，使用经典的层次分割方法进行实验对比，包括文献1中基于Wishart度量标准的迭代合并方法(简称IM-Wishart)和文献2中基于KummerU度量标准的迭代合并方法(简称IM-KummerU)。实验在通用计算机平台上运行，采用Matlab语言进行程序编程。

PolSAR数据为来自于德国宇航局(DLR)的ESAR数据，图像尺寸大小为510×530，图中包含均匀区域和非均匀区域，地物类型包含人造建筑物、森林、农田等。首先对PolSAR图像进行超像素分割，超像素个数为706，即N＝706。接着对超像素分割结果进行Wishart合并(步骤2)，合并后剩余区域的个数为353，即N⁽⁰⁾＝353。Wishart合并阶段将均匀区域中的大部分超像素进行合并，剩余的超像素基本位于非均匀区域，例如复杂建筑物区域以及图像边缘附近。对剩余区域进行KummerU迭代合并(步骤3)，在每一个迭代步骤中，将最小度量标准所对应的两个区域合并。基于KummerU能量总和SE的变化序列，步骤4计算出最终分割结果的区域个数N_r＝37。N⁽⁰⁾-N_r＝353-37＝316，由此得到最终的PolSAR图像分割结果为V⁽³¹⁶⁾。

对实施例的分割结果基于地面真值(ground truth)数据进行数值评估，一般使用精度和回召率两个指标来表征分割的准确性。精度(precision)指分割结果中正确边缘点的个数与所有区域边缘点的比值，而回召率(recall)表示分割结果中正确边缘点的个数与图像所有真实边缘点个数的比值。图2给出了本发明以及IM-Wishart方法和IM-KummerU方法在合并过程中得到的精度—回召率曲线，圆圈标记处为最终分割结果在曲线中的位置，可以看出本发明方法能取得更高的精度和回召率，表明该方法在合并过程中可以更好地保持图像细节和真实边缘。

表1列出了不同分割方法的运行时间。计算效率不仅取决于度量标准所采用的统计模型，同时也与区域合并的策略有关。本发明中，Wishart合并阶段对均匀区域的超像素进行快速合并，该处理过程只需要2到3秒。IM-Wishart和IM-KummerU方法从最原始的超像素分割结果开始进行迭代合并，IM-KummerU方法的运算时间超过700秒；而本发明中的KUMS阶段在Wishart合并所得结果的基础上执行迭代合并，时间消耗显著降低，仅为256秒。因此，基于两阶段的合并策略，本发明的计算效率得到明显提高。

表1不同方法的计算时间对比

综上，上述实验不仅验证了本发明中PolSAR图像分割方法的可行性和实用性，同时还表明了该方法能够有效提升分割精度和计算效率。

虽然参照上述实施例详细描述了本发明，但是应该理解本发明并不限于所公开的实施例。对于本专业领域的技术人员来说，可以对其形式和细节进行各种改变。本发明意欲涵盖所附权利要求书的精神和范围内的各种变形。

Claims (8)

1.一种基于两阶段区域合并策略的极化合成孔径雷达图像分割方法，首先对极化合成孔径雷达图像即PolSAR图像进行超像素分割，然后进行区域合并；其特征在于区域合并时先进行Wishart合并，再进行KummerU迭代合并，包括以下步骤：

步骤1：对PolSAR图像进行超像素分割，方法是：

采用PolSAR图像超像素分割方法对原始PolSAR图像进行超像素分割，得到超像素集合V＝{v₁,v₂,...,v_i,...,v_j,...,v_N}，其中N为超像素个数，v_i表示第i个超像素，1≤i≤N，1≤j≤N；超像素集合V及集合中元素的关系用区域邻接图G＝(V,A)来表示，其中A为超像素之间相邻关系的集合，a_ij∈A表示超像素v_i和v_j在几何位置上相邻，令N_a表示A中元素的个数，N_a为正整数；用b_ij来表示相邻超像素v_i和v_j之间的共享边界，为该边界上像素的集合；

步骤2：对V进行Wishart合并，方法是：

2.1对V中所有超像素v₁,v₂,...,v_i,...,v_j,…,v_N计算分布参数Z₁,Z₂,...,Z_i,...,Z_j,...,Z_N，分布参数为超像素中协方差矩阵C的均值，v_i的分布参数为Z_i，按公式(1)计算：

其中，|v_i|_c表示超像素v_i的基数，即v_i中像素的个数，下标c指代基数，s表示一个像素，C_s为像素s的协方差矩阵；

2.2根据A计算所有相邻超像素之间的Wishart度量标准，构建第一三元组集合D，具体步骤如下：

2.2.1令k＝1，集合D为空集；

2.2.2设A中第k个元素为a_ij，1≤i≤N，1≤j≤N，a_ij所对应的两个相邻超像素为v_i和v_j；

2.2.3计算超像素v_i和v_j合并前后产生的Wishart能量损失计算式为：

其中，v_ij＝v_i∪v_j表示超像素v_i和v_j合并之后的区域，|v_ij|_c为v_ij的基数，|v_i|_c和|v_j|_c分别为v_i和v_j的基数，Z_ij表示v_ij的分布参数，由公式(3)计算得出，|Z_ij|表示分布参数Z_ij的行列式值，|Z_i|和|Z_j|分别为分布参数Z_i和Z_j的行列式值；

2.2.4计算超像素v_i和v_j之间的边缘惩罚函数值计算式为：

其中，V_s表示像素s的边缘强度值，K为决定边缘惩罚强度的参数；

2.2.5将2.2.3步计算得到的Wishart能量损失和2.2.4步得到的边缘惩罚函数值相加，得到超像素v_i和v_j在Wishart合并阶段的度量标准按公式(5)计算：

其中，参数β是控制相对于的权重；

2.2.6令度量标准元素将d_k以及所对应的i和j构成三元组(d_k,i,j)，放入集合D中；

2.2.7令k＝k+1，判断k是否小于等于N_a，N_a为A中元素的个数，若满足，转2.2.2步；否则，得到了所有相邻超像素在Wishart合并阶段的度量标准以及在此基础上构成的三元组的集合D，D中元素个数为N_a，转2.3步；

2.3对D中三元组按三元组中第一个数值的大小进行排序，排列顺序为从小到大的顺序，取前M个三元组，构成集合D_M，1≤M≤N_a；

2.4根据D_M中元素，在A中找出所对应的M个元素，构成集合A_M，转2.5步；

2.5基于A_M对超像素集合V进行合并，得到合并后区域的集合V⁽⁰⁾和区域邻接图G⁽⁰⁾＝(V⁽⁰⁾,A⁽⁰⁾)，V⁽⁰⁾中元素为区域，设区域个数为N⁽⁰⁾，区域表示为1≤m≤N⁽⁰⁾，1≤n≤N⁽⁰⁾；V⁽⁰⁾及V⁽⁰⁾中元素的关系用区域邻接图G⁽⁰⁾＝(V⁽⁰⁾,A⁽⁰⁾)来表示，其中A⁽⁰⁾为Wishart合并后区域之间相邻关系的集合，令表示A⁽⁰⁾中元素的个数，为正整数，用b_mn来表示相邻区域r_m和r_n之间的共享边界；

步骤3：对V⁽⁰⁾进行KummerU迭代合并，方法是：

3.1对V⁽⁰⁾中每个区域，进行参数估计，并计算区域的KummerU能量，具体步骤如下：

3.1.1令k₂＝1；

3.1.2 V⁽⁰⁾中第k₂个元素为r_k，根据式(6)计算r_k的分布参数Z_k：

其中，|r_k|_c表示区域r_k的基数，即r_k中像素的个数；

3.1.3基于区域r_k中所有像素的协方差矩阵，估计区域r_k中KummerU分布的纹理形状参数ξ_k和ζ_k；

3.1.4计算区域r_k的KummerU能量计算式为：

其中，L为图像视数，d等于3，为矩阵的迹，Γ(x)是x的伽马函数值，x为自变量，为以Ld+ζ_k、Ld-ξ_k+1、为参数的第二类合流超几何函数；

3.1.5令k₂＝k₂+1，判断k₂是否小于等于N⁽⁰⁾，若满足，转3.1.2步；否则，得到V⁽⁰⁾中所有区域的KummerU能量，转3.2步；

3.2计算V⁽⁰⁾中所有区域的KummerU能量总和SE⁽⁰⁾：

3.3根据A⁽⁰⁾计算所有相邻区域之间的KummerU度量标准，构建第二三元组集合D⁽⁰⁾，具体步骤如下：

3.3.1令k₃＝1，集合D⁽⁰⁾为空集；

3.3.2设A⁽⁰⁾中第k₃个元素为a_mn，1≤m≤N⁽⁰⁾，1≤n≤N⁽⁰⁾，a_mn所对应的两个相邻区域为r_m和r_n；

3.3.3计算区域r_m和r_n合并前后产生的KummerU能量损失计算式为：

其中，r_mn＝r_m∪r_n表示区域r_m和r_n合并之后的区域，Z_mn为区域r_mn的分布参数，由公式(11)计算得出，ξ_mn和ζ_mn为区域r_mn的纹理形状参数，基于区域r_mn中所有像素的协方差矩阵估计得到；

3.3.4计算区域r_m和r_n之间的边缘惩罚函数值计算式为：

3.3.5分别计算区域r_m、区域r_n以及合并后区域r_mn的匀质度H_m，H_n和H_mn，并计算区域合并的匀质度惩罚因子F_mn，计算式为：

3.3.6将3.3.3步计算得到的KummerU能量损失3.3.4步得到的边缘惩罚函数值以及3.3.5步得到的匀质度惩罚因子F_mn相结合，得到区域r_m和r_n在KummerU合并阶段的度量标准按公式(14)计算：

3.3.7令第二度量标准元素将以及所对应的m和n构成三元组放入集合D⁽⁰⁾中；

3.3.8令k₃＝k₃+1，判断k₃是否小于等于 为A⁽⁰⁾中元素的个数，若满足，转3.3.2步；否则，得到所有相邻区域在KummerU合并阶段的度量标准以及在此基础上构成的三元组的集合D⁽⁰⁾，D⁽⁰⁾中元素个数为转3.4步；

3.4令l表示KummerU迭代合并的次数，且l＝0；

3.5找出D^(l)里所有三元组中第一个数值最小的三元组，设该三元组第二和第三个数值为m和n，则对区域r_m和r_n进行合并，合并后区域仍记为区域r_m；同时，将D^(l)中与区域r_m和r_n相关的三元组从集合中删除，将删除元素之后的D^(l)记为D^(l+1)；

3.6更新3.5步区域合并之后的区域邻接图，用G^(l+1)＝(V^(l+1),A^(l+1))来表示，并更新合并后的能量总和SE^(l+1)和D^(l+1)：

3.6.1区域合并后的能量总和SE^(l+1)更新为：

其中，SE^(l)为区域合并前的KummerU能量总和，为第3.5步合并过程中产生的KummerU能量损失；

3.6.2基于区域r_m中所有像素的协方差矩阵，对r_m的分布参数Z_m以及纹理形状参数ξ_m和ζ_m进行估计；

3.6.3将A^(l+1)中与区域r_m相关的相邻关系元素构成设集合中元素个数为对中元素逐个计算相应的第二度量标准，并对D^(l+1)进行更新；

3.6.4对KummerU迭代合并次数l进行更新：l＝l+1，判断l是否小于等于N⁽⁰⁾-1，若满足，转3.5步；否则，转步骤4；

步骤4：确定PolSAR图像最终分割结果的区域个数：基于SE序列确定理想分割结果所对应的区域个数，得到区域个数为N_r基于区域个数N_r确定最终的分割结果为

2.如权利要求1所述的一种基于两阶段区域合并策略的极化合成孔径雷达图像分割方法，其特征在于所述V_s介于0到1之间，K取值为0.3，β取值为5。

3.如权利要求1所述的一种基于两阶段区域合并策略的极化合成孔径雷达图像分割方法，其特征在于所述M的取值根据超像素个数N以及百分比参数a来设定，M＝N×a％，a取值范围为40～60。

4.如权利要求1所述的一种基于两阶段区域合并策略的极化合成孔径雷达图像分割方法，其特征在于2.4步所述根据D_M中三元组，在A中找出所对应的M个元素，构成集合A_M的方法是：

2.4.1令k₁＝1，A_M为空集；

2.4.2取出D_M中第k₁个三元组，设该三元组第二和第三个数值为i和j，1≤i≤N，1≤j≤N，将a_ij放入集合A_M；

2.4.3令k₁＝k₁+1，判断k₁是否小于等于M，若满足，转2.4.2步；否则，得到包含M个元素的集合A_M。

5.如权利要求1所述的一种基于两阶段区域合并策略的极化合成孔径雷达图像分割方法，其特征在于2.5步所述基于A_M对超像素集合V进行合并的方法是：

2.5.1令k₁＝1；

2.5.2设A_M中第k₁个元素为a_ij，1≤i≤N，1≤j≤N，a_ij所对应的两个相邻超像素为v_i和v_j，则将超像素v_i和v_j合并；

2.5.3令k₁＝k₁+1，判断k₁是否小于等于M，若满足，转2.5.2步；否则，即完成了Wishart合并。

6.如权利要求1所述的一种基于两阶段区域合并策略的极化合成孔径雷达图像分割方法，其特征在于步骤3.3.5计算区域r_m、区域r_n以及合并后区域r_mn的匀质度H_m，H_n和H_mn采用变化系数的方法。

7.如权利要求1所述的一种基于两阶段区域合并策略的极化合成孔径雷达图像分割方法，其特征在于步骤3.6.3所述更新D^(l+1)的方法是：

3.6.3.1令k₄＝1；

3.6.3.2设中第k₄个元素为a_mn，计算KummerU能量损失边缘惩罚函数和匀质度惩罚因子F_mn，并按公式(14)计算得到度量标准

3.6.3.3令第二度量标准元素将以及所对应的m和n构成三元组放入集合D^(l+1)中；

3.6.3.4令k₄＝k₄+1，判断k₄是否小于等于若满足，转3.6.3.2步；否则，完成了对D^(l+1)的更新。

8.如权利要求1所述的一种基于两阶段区域合并策略的极化合成孔径雷达图像分割方法，其特征在于步骤4所述基于SE序列确定PolSAR图像最终分割结果的区域个数的方法是：构建一种评估曲线图，图的x轴为区域个数，y轴表示KummerU能量总和SE，曲线图的拐点，即最大曲率所对应的点则为最终分割结果的区域个数。