CN107808380B - 一种基于G0和Gamma联合分布的多尺度SAR影像水体分割方法 - Google Patents

Abstract

本发明在图像金字塔结合水平集理论的多尺度迭代处理框架内，采用G⁰和Gamma联合分布构建全局能量泛函，并在金字塔顶层SAR图像处理时引入快速模糊C均值(fast fuzzy C means，FFCM)算法初始化水平集函数，提出一种多尺度的基于G⁰和Gamma联合分布的水体分割新方法。本发明基于G⁰和Gamma联合分布的变分水平集分割，对其中涉及的水平集初始化方法进行了进一步的改进。

Description

一种基于G0和Gamma联合分布的多尺度SAR影像水体分割方法

技术领域

本发明涉及SAR遥感领域，尤其涉及一种基于G⁰和Gamma联合分布的多尺度高分辨率SAR 影像水体分割方法

背景技术

合成孔径雷达(synthetic aperture radar,SAR)具备全天时、全天候的优势，能实现大范围、实时、动态的地物监测。与光学传感器相比，其无惧风雨，穿云透雾的全天候特点，更适用于洪涝灾害信息的监测和提取。

由于SAR图像的相干斑噪声特点，在分割方法中引入SAR统计特性的基于统计模型的分割方法正受到越来越多的重视。传统的基于区域分裂/合并的分割方法虽然可以和SAR统计模型结合，但其只能基于SAR局域操作，无法获取全局最优分割结果。基于全局的分割方法一般将SAR图像分割看做一个全局的问题，通过优化一个全局的目标函数来实现分割。基于 MRF模型的方法和变分水平集方法是其中的代表。变分水平集方法可以很好地处理拓扑结构的变化，提供统一的框架结合SAR统计信息和先验信息，与MRF方法相比，变分水平集方法采用的变分模型基于连续图像域，可以避免MRF模型离散网格出现的网格偏向误差。

2005年Ben在水平集框架内引入描述相干斑噪声的Gamma模型，在不需要相干斑预处理的情况下，获得SAR图像比较准确的分割结果。2014年徐川等人在将多尺度分析技术与变分水平集方法结合，在多尺度水平集框架内，采用otsu初始化初次迭代的水平集函数，利用 Gamma统计模型构造能量泛函，实现SAR影像水体的分割。

SAR多视强度图像的统计特性一般会建模成Gamma模型。但是，随着SAR影像分辨率的提高，地物在SAR图像上呈现出越来越多的细节，许多在低分辨率SAR图像上展现的较为均质的地物在高分辨SAR图像(分辨率在10m以内)上呈现出不均质的状态，其统计特性已不再完全符合Gamma模型。分辨率较低时，SAR图像单位像元尺寸较大，一个像元可以包含大量的随机散射点，使单位像元后向散射特征展现为较均质的状态。而高分辨SAR图像相干斑发育不充分，像元尺寸小于强散射体的尺寸，纹理特征趋于明显，强散射点更突出，图像特征趋于不均质。为了描述高分辨率SAR图像的统计特性，许多统计模型被提出，如k 分布、广义Gamma Reyleigh分布等。其中，Frery基于乘性噪声模型，在相干斑分量服从单位均值Gamma分布，背景雷达散射截面积(Radar Cross Section，RCS)服从广义高斯分布的假设下推导出G⁰模型。2014年Cui等人对目前常用的SAR统计模型在同一场景不同分辨率和不同场景同一分辨率两种情况下建模能力的对比实验揭示了G⁰分布模型不但具有广泛均匀度变化下强大的建模能力，也具有较好的模型向下兼容性。

水体在高分辨率SAR图像上表现为质地均匀的面状区域，相干斑发育充分，可以用Gamma分布拟合。而水体周围通常杂糅分布着人工建筑物和混合植被等质地复杂的地物，水体周边地形也不尽是平坦，其统计特性已不符合Gamma模型。

现有的SAR图像的基于统计模型的水平集水体分割方法多采用Gamma建模，不能适应高分辨率SAR数据处理的需要，分割精度较低。其次，现有的方法多采用otsu等二值分割结果初始化水平集函数，使用该种方法初始化水平集函数，极容易使最终分割结果受灰度值与水体接近的阴影区域干扰，导致分割精度低。

发明内容

本发明在图像金字塔结合水平集理论的多尺度迭代处理框架内，采用G⁰和Gamma联合分布构建全局能量泛函，并在金字塔顶层SAR图像处理时引入快速模糊C均值(fastfuzzy C means，FFCM)算法初始化水平集函数，提出一种多尺度的基于G⁰和Gamma联合分布的水体分割新方法。

该方法首先对SAR图像进行金字塔分解，并采用G⁰分布拟合SAR水体周边背景区，Gamma分布模型拟合水体区，对SAR图像整体综合建模并构建全局能量泛函；然后在顶层SAR图像上，用FFCM聚类初始化水平集函数；接着基于区域竞争模型，在能量泛函的梯度下降流推动下，迭代实现能量泛函的最小化，得到首层分割结果；该层分割结果反演至下一层，初始化下一层分割。这样，分割结果在多尺度迭代链中不断精化，最终得到轮廓细节完整的水体区域。

本发明所提出的基于G⁰和Gamma联合分布的多尺度高分辨率SAR影像水体分割方法的总体框架。分为以下步骤：1)金字塔分解SAR影像至K层，形成0,1,2,…,i,…,K层图像序列，其中第0层为原始图像，分辨率最高，从第0层至第K层分辨率逐步下降，第K层为最粗尺度图像，分辨率最低；2)取第i层SAR影像，i的初始取值为K；3)对第i层影像执行 G⁰和Gamma联合分割算法，得到分割结果；4)若i＞0,则将第i层影像的分割结果反演至第i-1层，然后将i的值赋值为i-1，转回第3)步继续执行分割；否则，执行第5)步；5) 输出分割结果并评价其精度。

优选的是，步骤3)中，采用G⁰分布拟合SAR水体周边背景区，Gamma分布模型拟合水体区，采用G⁰和Gamma联合分布对SAR图像整体综合建模并构建全局能量泛函；然后在K层SAR图像上，用FFCM聚类初始化水平集函数，在i≠K层，由第i+1层分割结果作为第i 层的初始化；估算水平集函数确定的水体内、外区域满足的Gamma、G⁰模型的参数；接着在能量泛函的梯度下降流推动下，迭代实现能量泛函的最小化，得到水体区域的分割结果。进一步优选的是，采用G⁰和Gamma联合分布对SAR图像整体综合建模并构建全局能量泛函具体包括：

假定SAR图像Ω被划分为水体内外两个区域，即为水体内目标区域Ω_F和水体外背景区域Ω_B。

令I(x，y)表示SAR图像强度数据，目标区域Ω_F满足的Gamma分布概率密度公式为：

其中，参数n为等效视数，u为均值,Γ是Gamma函数；

背景区域Ω_B满足的G⁰分布概率密度公式[2]为：

(x，y)∈Ω_B

-α，γ，n＞0，(x，y)∈Ω_B

其中，Γ是Gamma函数；参数n为等效视数；参数α为形状参数，反映被测区的均匀度，α∈(-∞，0)；参数γ为尺度参数，反映被测区的综合能量；

令P_f＝P_f(I(x，y)，n，u)，(x，y)∈Ω_F，P_b＝P_b(I(x，y)，n，u)，(x，y)∈Ω_B，则SAR图像强度数据可以用两个模型的混合概率模型表示：

P＝ω_fP_f+ω_bP_b (3)

ω_f和ω_b分别为目标区域和背景区域的先验概率，满足ω_f+ω_f＝1；

SAR图像概率分布的似然函数对数形式为：

利用-log(l_FB)构建能量泛函：

其中，η，μ，λ≥0；式(5)第一项是距离正则项，该项系数η经验值为0.04；第二项为曲线长度项，曲线c是水体边缘轮廓线，噪声越多的图像，系数μ越大；第三项为似然函数项，用统计模型来拟合图像，λ为系数。

进一步优选的是，在水平集函数初始化的基础上，迭代实现能量泛函的最小化，得到K 层分割结果的具体步骤包括：

引入Heaviside函数的近似形式H_ε以及δ_ε＝H′_ε函数，描述为以下式子(7)形式，ε一般取 1.5：

水平集函数经过初始化，在水体内区域为正值，在水体外区域为负值，如式子(8)所示：

其中，inside(c)为水体内区域，outside(c)为水体外区域，φ(x，y)为水平集函数；

令φ＝φ(x，y)，I＝I(x，y)，将式子(7)、(8)带入表示能量泛函的式子(5)中，式子(5)变为：

是梯度运算符；式子(9)中第一项距离正则项描述为以下形式：

在变分水平集的框架下，能量泛函的梯度下降流如式子(12)所示：

将式子(9)、(10)、(11)带入式子(12)中，得到梯度下降流为：

div(·)指散度；在初始条件φ(x，y，t)|_t＝0＝φ₀(x，y)下，解方程(13)得到解φ，其对应的φ(x，y)≥0区域就是所要提取的水体区域。

进一步优选的是，在实际计算时，用迭代的方式来逐渐逼近能量泛函的最小值，即设k 为迭代次数，Δt为时间间隔，则

将式子(13)代入式子(14)得到实现能量泛函最小化的迭代求解公式。

优选的是，步骤3中用FFCM聚类水平集函数初始化方法具体如下：

将SAR图像映射到灰度空间，再对灰度级运用FCM聚类；其目标函数为：

其中n为图像实际灰度级个数，c为类别数，h_k为第k级灰度，n(h_k)为图像中灰度级为 h_k的像素个数；平滑参数m＝2，μ_ik为灰度级为h_k的样本x(h_k)到类别i的隶属度；d_ik为样本到类别i的距离，这里取欧几里德距离，即：

d_ik＝||x(h_k)-P_i)^T||₂ (16)

变分法使目标函数取极小值的必要条件为：

进一步优选的是，具体的聚类步骤如下：

4)输入类别c，2≤c≤n，迭代停止阈值ε，聚类初始中心P_i ⁽⁰⁾，1≤i≤c；

5)根据式子(17)计算所有样本对当前聚类中心的隶属度；

6)根据步骤2)计算出隶属度利用式子(18)重新计算当前聚类中心；

重复步骤2)，3)，直到相邻两部分聚类中心小于阈值ε；

得到聚类结果或分割结果后，按式子(19)初始化水平集函数：

其中，曲线c为FFCM聚类结果确定的水体轮廓线，d_xy为像素点(x,y)距离最近的非0值的欧几里德距离；inside(c)指水体轮廓内区域，outside(c)指水体轮廓外区域。

优选的是，对于G⁰模型待估计的形状参数α、尺度参数γ和等效视数n，采用基于Mellin变换的参数估计方法，包括：

从Mellin变换定义函数的第二类第一特征函数中导出第二类第二特征函数和前三阶第二特征对数累积量，取对数累积量表达式组成如下非线性方程组：

其中，Ψ(z)＝dlgΓ(z)/dz为Digamma函数，Ψ_p(z)＝d^p+1lgΓ(z)/dz^p+1为Polygamma函数； c_p，p＝1,2,3为p阶对数累积量，可通过如下表达式求得：

x₁，x₂，x₃...x_N为N个用于估算G⁰分布参数的图像区域像素值；利用数值迭代的方法可以很快收敛，解出该非线性方程组的唯一解。

Gamma模型待估计的等效视数n和图像均值u。其中视数n由先验知识获取，均值u由最大似然估计算法得到如下计算方式：

x₁，x₂，x₃...x_N为N个用于估算均值u的图像区域像素值。

本发明利用G⁰和Gamma联合分布拟合高分辨率SAR影像统计特性并构建全局能量泛函，在多尺度水平集方法框架下实现能量泛函的最小化，最终实现水体的分割。引入FFCM聚类改进了水平集函数的初始化，降低了统计特性与水体类似的阴影区域对分割结果的干扰；引入 G⁰和Gamma联合分布改进了全局能量泛函，改善了场景复杂区域的水体分割结果。与传统分割方法相比，本发明在水体分割精度上有较好的提高，是一种较为实用的高分辨率SAR影像水体分割方法。

说明书附图

图1是本发明形成的多尺度迭代链示意图；

图2是本发明所述分割方法的框架示意图；

图3是水平集函数零水平集的示意图；

图4是水平集函数初始化的示意图。

具体实施方式

下面通过实施例，对本发明的技术方案做进一步具体的说明。

本发明在图像金字塔结合水平集理论的多尺度迭代处理框架内，采用G⁰和Gamma联合分布构建全局能量泛函，并在金字塔顶层SAR图像处理时引入快速模糊C均值(fastfuzzy C means，FFCM)算法初始化水平集函数，提出一种多尺度的基于G⁰和Gamma联合分布的水体分割新方法。该方法首先对SAR图像进行金字塔分解，并采用G⁰分布拟合SAR水体周边背景区，Gamma分布模型拟合水体区，对SAR图像整体综合建模并构建全局能量泛函；然后在顶层SAR图像上，用FFCM聚类初始化水平集函数；接着基于区域竞争模型，在能量泛函的梯度下降流推动下，迭代实现能量泛函的最小化，得到首层分割结果；该层分割结果反演至下一层，初始化下一层分割。这样，分割结果在图1所示的多尺度迭代链中不断精化，最终得到轮廓细节完整的水体区域。

图2所示为本发明所提出的基于G⁰和Gamma联合分布的多尺度高分辨率SAR影像水体分割方法的总体框架。分为以下步骤：1)金字塔分解SAR影像至K层，形成0,1,2,…,i,…,K 层图像序列，其中第0层为原始图像，分辨率最高，从第0层至第K层分辨率逐步下降，第 K层为最粗尺度图像，分辨率最低；2)取第i层SAR影像，i的初始取值为K；3)对第i 层影像执行G⁰和Gamma联合分割算法，得到分割结果；4)若i＞0,则将第i层影像的分割结果反演至第i-1层，然后将i的值赋值为i-1，转回第3)步继续执行分割；否则，执行第5) 步；5)输出分割结果并评价其精度。

图2所示的总体框架中提到的基于G⁰和Gamma联合分布的变分水平集分割是本发明所提出方法的核心。其具体实现过程下面会详细阐述。

1.水平集演化理论

水平集方法是20世纪80年代由美国数学家Stanley Osher和James Sethian中提出的。基本思想是：将轮廓线看成是一个函数的零水平集，这样的函数称为水平集函数。由相关准则定义关于水平集函数的能量泛函，在变分水平集的框架下，驱动水平集函数向能量减少的方向演化，当能量最小时，此时对应的水平集函数的零水平集为最终分割轮廓线。

设t时刻轮廓线c(t)是水平集函数φ(x，y，t)的零水平集，则c(t)＝{(x，y)|φ(x，y，t)＝0}。即，轮廓线c是一系列水平集函数值为0的点构成的，如图3所示。令φ＝φ(x，y，t)，能量泛函E(φ)是φ(x，y，t)的函数，不同时刻t对应不同水平集函数φ，不同水平集函数φ对应不同能量E，当能量达到最小时，此时水平集函数φ对应的零水平集c就是水体轮廓线。

从水平集演化过程可以看出，构建出符合SAR水体内外图像特性的能量泛函是实现水体分割的关键，因此，本文将G⁰分布和Gamma分布联合建模，构建出符合水体内外统计特性的能量函数。

2.基于G⁰和Gamma联合分布的变分水平集分割

本发明综合考虑高分辨率SAR图像的散射特点，采用G⁰和Gamma联合分布对SAR图像建模。将相干斑发育充分的水体区域采用Gamma分布建模，对于水体外区域，通常集中或散落分布着人工建筑区、混合植物区，采用G⁰分布建模。虽然G⁰分布对均匀区域同样有很好的建模能力，但考虑到G⁰模型计算较为复杂，参数估算也比Gamma模型繁琐，故采用综合建模的方式。

假定SAR图像Ω被划分为水体内外两个区域，即为水体内目标区域Ω_F和水体外背景区域Ω_B。

令I(x，y)表示SAR图像强度数据，目标区域Ω_F满足的Gamma分布概率密度公式为：

其中，参数n为等效视数，u为均值,Γ(n)是Gamma函数，表达式为：

背景区域Ω_B满足的G⁰分布概率密度公式[2]为：

(x，y)∈Ω_B

-α，γ，n＞0，(x，y)∈Ω_B

其中，Γ是Gamma函数；参数n为等效视数；参数α为形状参数，在本质上反映被测区的均匀度，α∈(-∞，0)；被测区越不均匀，α越大,即SAR场景包含植被、建筑物时，参数α往往较大；参数γ为尺度参数，反映被测区的综合能量。

Gamma模型是全发展模型，G⁰模型具有很好的模型向下兼容性。经过金子塔分解后的图像，其统计特性仍可以由分解前的模型拟合。

令P_f＝P_f(I(x，y)，n，u)，(x，y)∈Ω_F，P_b＝P_b(I(x，y)，n，u)，(x，y)∈Ω_B，则SAR图像强度数据可以用两个模型的混合概率模型表示：

P＝ω_fP_f+ω_bP_b (3)

ω_f和ω_b分别为目标区域和背景区域的先验概率，满足ω_f+ω_f＝1。

SAR图像概率分布的似然函数对数形式为：

通过最大化上述似然函数，就可以实现目标区和背景区的分割，最大化似然函数就是最小化-log(l_FB)。在水平集演化理论中，最小化水平集函数的能量泛函，就可以实现分割。因此，这里直接利用-log(l_FB)构建能量泛函，加上长度项，并加上距离正则项，最终能量泛函为式子 (5)所示：

其中，η，μ，λ≥0。第一项是距离正则项，加上此项后，水平集函数在演化过程中可以保持距离符号函数的特性，不需要重新初始化。该项系数η不能选的过大，否则会影响算法稳定性，一般经验值为0.04。第二项为曲线长度项，使演化曲线光滑。曲线c是水体边缘轮廓线，在水平集理论中抽象为零水平集。噪声越多的图像，系数μ要尽量选的大些，以免演化曲线陷入局部极小值。第三项为似然函数项，用统计模型来拟合图像，λ为系数。

引入Heaviside函数：

由于Heaviside函数不连续，在实际计算中，通常使用它的近似形式H_ε以及δ_ε＝H′_ε函数，描述为以下式子(7)形式，ε一般取1.5：

水平集函数经过初始化，在水体内区域为正值，在水体外区域为负值，如图4所示，如式子(8)所示：

其中，inside(c)为水体内区域，outside(c)为水体外区域，φ(x，y)为水平集函数。

令φ＝φ(x，y)，I＝I(x，y)，将式子(7)、(8)带入式子(5)中，式子(5)变为：

是梯度运算符。式子(9)中第一项距离正则项描述为以下形式：

在变分水平集的框架下，能量泛函的梯度下降流如式子(12)所示：

将式子(9)、(10)、(11)带入式子(12)中，得到梯度下降流为：

div(·)指散度。在初始条件φ(x，y，t)|_t＝0＝φ₀(x，y)下，解方程(13)得到解φ，其对应的φ(x，y)≥0区域就是所要提取的水体区域。

在实际计算时，上述公式很难直接求解，一般采用迭代的方式来逐渐逼近能量泛函的最小值。设k为迭代次数，Δt为时间间隔，则

将式子(13)代入式子(14)得到实现能量泛函最小化的迭代求解公式。

3.水平集函数初始化方法

在金字塔分解的最粗尺度上，本文引入FFCM聚类，利用聚三类中水体类初始化水平集函数，削弱阴影对最终分割结果的影响。之所以将SAR图像聚成三类，而不是简单采取二值分割，是考虑到水体周边地形起伏会在SAR图像上形成阴影。阴影在SAR图像上灰度值较低，但与水体仍有区别。如果仅利用二值分割结果初始化水平集函数，会使最终分割结果受灰度值与水体接近的阴影区域干扰，导致分割精度降低。

模糊C均值(fuzzy C means，FCM)聚类是基于目标函数的聚类方法。FFCM与FCM不同之处在于前者是基于像素的，后者是基于灰度级的。即FFCM是将SAR图像映射到灰度空间，再对灰度级运用FCM聚类。其目标函数为：

其中n为图像实际灰度级个数，c为类别数，h_k为第k级灰度，n(h_k)为图像中灰度级为h_k的像素个数。平滑参数m＝2，μ_ik为灰度级为h_k的样本x(h_k)到类别i的隶属度。d_ik为样本到类别i的距离，这里取欧几里德距离，即：

d_ik＝||x(h_k)-P_i)^T||₂ (16)

变分法使目标函数取极小值的必要条件为：

具体的聚类步骤如下：

7)输入类别c，2≤c≤n，迭代停止阈值ε，聚类初始中心P_i ⁽⁰⁾，1≤i≤c；

8)根据式子(17)计算所有样本对当前聚类中心的隶属度；

9)根据步骤2)计算出隶属度利用式子(18)重新计算当前聚类中心；

10)重复步骤(2)(3)，直到相邻两部分聚类中心小于阈值ε。

FFCM聚类利用灰度直方图代替像素点进行聚类，减少了参与聚类的样本数量，降低了隶属度矩阵和聚类中心的计算复杂度，与传统FCM算法相比，极大地缩短了图像的分割时间。同时，FFCM算法也克服了传统FCM算法对噪声敏感的缺点，适用于相干斑噪声显著的SAR影像聚类问题。

在金字塔分解的其他尺度SAR图像上，均由上一尺度得到的分割结果初始化下一层分割。

得到聚类结果或分割结果后，按式子(19)初始化水平集函数，如图4所示：

其中，曲线c为FFCM聚类结果确定的水体轮廓线，d_xy为像素点(x,y)距离最近的非0值的欧几里德距离。inside(c)指水体轮廓内区域，outside(c)指水体轮廓外区域。

4.G⁰和Gamma联合分布参数估计

G⁰模型待估计参数有三个，分别为形状参数α、尺度参数γ和等效视数n。传统的参数估计方法主要有矩估计和最大似然估计两种。矩估计方法虽然计算简便，但只能估算形状参数α≤-2时的SAR场景，不能实现全范围的参数估计。而形状参数与SAR场景均匀度有直接关系，场景越不均匀，α越大。在本文中主要利用G⁰模型对不均匀场景的建模能力，显然矩估计方法不满足本文研究的需要。最大似然估计方法无法求解出被估参数的解析式，甚至无法将被估参数表示成样本均值等样本数字特征的隐函数，且求解方程高度非线性化，这不但造成计算量的大幅度攀升，更严重的是，有时根本无法通过数值迭代的方法得出正确的参数估计值，显然，该方法不可行。

综合考虑，本文选择基于Mellin变换的参数估计方法。该方法中，以Mellin变换为基础的第二类型统计量可以将相干噪声分量视为Mellin卷积。这样，计算结果就被极大地简化了。Tison采用Mellin变换和对数累积量的方法对Fisher分布的参数有了较好的估计。考虑到G⁰分布作变量变换后可以得到Fisher分布的性质，相似地，从Mellin变换定义函数的第二类第一特征函数中导出第二类第二特征函数和前三阶第二特征对数累积量，取对数累积量表达式组成如下非线性方程组：

其中，Ψ(z)＝dlgΓ(z)/dz为Digamma函数，Ψ_p(z)＝d^p+1lgΓ(z)/dz^p+1为Polygamma函数。 c_p，p＝1,2,3为p阶对数累积量，可通过如下表达式求得：

x₁，x₂，x₃...x_N为N个用于估算G⁰分布参数的图像区域像素值。由于Digamma函数和Polygamma函数的单调性，利用数值迭代的方法可以很快收敛，解出该非线性方程组的唯一解。

Gamma模型参数估计比较简单，只涉及两个参数，一个是等效视数n，一个是图像均值u。其中视数n由先验知识获取，均值u由最大似然估计算法得到如下计算方式：

x₁，x₂，x₃...x_N为N个用于估算均值u的图像区域像素值。

综上所述，基于G⁰和Gamma联合分布的变分水平集分割实现步骤如下：

1)基于G⁰和Gamma联合分布的全局能量泛函构建

2)水平集函数初始化：在第K层，由FFCM快速聚类结果初始化水平集函数；在i≠K层，由第i+1层分割结果作为第i层的初始化

3)联合分布参数估算：分别估算水平集函数确定的水体内外区域满足的Gamma、G⁰模型参数。

4)区域边界演化：由能量泛函的梯度下降流推动水平集函数完成一次演化，更新水平集函数，确立新的水体区域

5)判断Maxlter为最大迭代次数，k为迭代次数。若k≤Maxlter，则执行步骤3)4)；否则，执行步骤6)。

6)输出水体结果。

以上实施例仅用于说明本发明，而并非对本发明的限制，有关技术领域的普通技术人员，在不脱离本发明的精神和范围的情况下，还可以做出各种变化和变型，因此所有等同的技术方案也属于本发明的范畴，本发明的专利保护范围应由权利要求限定。

Claims (7)

1.一种基于G⁰和Gamma联合分布的多尺度高分辨率SAR影像水体分割方法，其特征在于，分为以下步骤：

1)金字塔分解SAR影像至K层，形成0，1，2，…，i，…，K层图像序列，其中第0层为原始图像，分辨率最高，从第0层至第K层分辨率逐步下降，第K层为最粗尺度图像，分辨率最低；

2)取第i层SAR影像，i的初始取值为K；

3)对第i层影像执行G⁰和Gamma联合分割算法，得到分割结果；

4)若i＞0，则将第i层影像的分割结果反演至第i-1层，然后将i的值赋值为i-1，转回第3)步继续执行分割；否则，执行第5)步；

5)输出分割结果并评价其精度；

其中，步骤3)中，采用G⁰分布拟合SAR水体周边背景区，Gamma分布模型拟合水体区，采用G⁰和Gamma联合分布对SAR图像整体综合建模并构建全局能量泛函；然后在K层SAR图像上，用FFCM聚类初始化水平集函数，在i≠K层，由第i+1层分割结果作为第i层的初始化；估算水平集函数确定的水体内、外区域满足的Gamma、G⁰模型的参数；接着在能量泛函的梯度下降流推动下，迭代实现能量泛函的最小化，得到水体区域的分割结果；

采用G⁰和Gamma联合分布对SAR图像整体综合建模并构建全局能量泛函具体包括：

假定SAR图像Ω被划分为水体内外两个区域，即为水体内目标区域Ω_F和水体外背景区域Ω_B；

令I(x，y)表示SAR图像强度数据，目标区域Ω_F满足的Gamma分布概率密度公式为：

其中，参数n为等效视数，u为均值，Γ是Gamma函数；

背景区域Ω_B满足的G⁰分布概率密度公式为：

(x，y)∈Ω_B

-α，γ，n＞0，(x，y)∈Ω_B

其中，Γ是Gamma函数；参数n为等效视数；参数α为形状参数，反映被测区的均匀度，α∈(-∞，0)；参数γ为尺度参数，反映被测区的综合能量；

令P_f＝P_f(I(x，y)，n，u)，(x，y)∈Ω_F，P_b＝P_b(I(x，y)，n，u)，(x，y)∈Ω_B，则SAR图像强度数据可以用两个模型的混合概率模型表示：

P＝ω_fP_f+ω_bP_b (3)

ω_f和ω_b分别为目标区域和背景区域的先验概率，满足ω_f+ω_f＝1；

SAR图像概率分布的似然函数对数形式为：

利用-log(l_FB)构建能量泛函：

其中，η，μ，λ≥0；式(5)第一项是距离正则项，该项系数η经验值为0.04；第二项为曲线长度项，曲线c是水体边缘轮廓线，噪声越多的图像，系数μ越大；第三项为似然函数项，用统计模型来拟合图像，λ为系数。

2.根据权利要求1所述的分割方法，其特征在于，在水平集函数初始化的基础上，迭代实现能量泛函的最小化，得到K层分割结果的具体步骤包括：

引入Heaviside函数的近似形式H_ε以及δ_ε＝H′_ε函数，描述为以下式子(7)形式，ε一般取1.5：

水平集函数经过初始化，在水体内区域为正值，在水体外区域为负值，如式子(8)所示：

其中，inside(c)为水体内区域，outside(c)为水体外区域，φ(x，y)为水平集函数；

令φ＝φ(x，y)，I＝I(x，y)，将式子(7)、(8)带入表示能量泛函的式子(5)中，式子(5)变为：

是梯度运算符；式子(9)中第一项距离正则项描述为以下形式：

在变分水平集的框架下，能量泛函的梯度下降流如式子(12)所示：

将式子(9)、(10)、(11)带入式子(12)中，得到梯度下降流为：

div(·)指散度；在初始条件φ(x，y，t)|_t＝0＝φ₀(x，y)下，解方程(13)得到解φ，其对应的φ(x，y)≥0区域就是所要提取的水体区域。

3.根据权利要求2所述的分割方法，其特征在于，用迭代的方式来逐渐逼近能量泛函的最小值，即设k为迭代次数，Δt为时间间隔，则

将式子(13)代入式子(14)得到实现能量泛函最小化的迭代求解公式。

4.根据权利要求3所述的分割方法，其特征在于，步骤3中用FFCM聚类水平集函数初始化方法具体如下：

将SAR图像映射到灰度空间，再对灰度级运用FCM聚类；其目标函数为：

其中n为图像实际灰度级个数，c为类别数，h_k为第k级灰度，n(h_k)为图像中灰度级为h_k的像素个数；平滑参数m＝2，μ_ik为灰度级为h_k的样本x(h_k)到类别i的隶属度；d_ik为样本到类别i的距离，这里取欧几里德距离，即：

d_ik＝||(x(h_k)-P_i)^T||₂ (16)

变分法使目标函数取极小值的必要条件为：

5.根据权利要求4所述的分割方法，其特征在于，具体的聚类步骤如下：

1)输入类别c，2≤c≤n，迭代停止阈值ε，聚类初始中心P_i ⁽⁰⁾，1≤i≤c；

2)根据式子(17)计算所有样本对当前聚类中心的隶属度；

3)根据步骤2)计算出隶属度利用式子(18)重新计算当前聚类中心；

重复步骤2)，3)，直到相邻两部分聚类中心小于阈值ε；

得到聚类结果或分割结果后，按式子(19)初始化水平集函数：

其中，曲线c为FFCM聚类结果确定的水体轮廓线，d_xy为像素点(x，y)距离最近的非0值的欧几里德距离；inside(c)指水体轮廓内区域，outside(c)指水体轮廓外区域。

6.根据权利要求5所述的分割方法，其特征在于，对于G⁰模型待估计的形状参数α、尺度参数γ和等效视数n，采用基于Mellin变换的参数估计方法，包括：

从Mellin变换定义函数的第二类第一特征函数中导出第二类第二特征函数和前三阶第二特征对数累积量，取对数累积量表达式组成如下非线性方程组：

其中，Ψ(z)＝dlgΓ(z)/dz为Digamma函数，Ψ_p(z)＝d^p+1lgΓ(z)/dz^p+1为Polygamma函数；c_p，p＝1，2，3为p阶对数累积量，可通过如下表达式求得：

x₁，x₂，x₃...x_N为N个用于估算G⁰分布参数的图像区域像素值；利用数值迭代的方法可以很快收敛，解出该非线性方程组的唯一解。

7.根据权利要求6所述的分割方法，其特征在于，Gamma模型待估计的等效视数n和图像均值u按如下方法获得：

其中视数n由先验知识获取，均值u由最大似然估计算法得到如下计算方式：

x₁，x₂，x₃...x_N为N个用于估算均值u的图像区域像素值。

Images