合成孔径雷达(synthetic aperture radar，SAR)是一种借助于机载或星载平台实现高分辨率成像的遥感系统，通过雷达波束向目标地物发射相位调制脉冲，整合地面反射回来的后向散射能量实现成像和对地观测^[1]。由于雷达自带照射源，因此它可全天时、全天候工作，且可几乎无失真地穿过可透射电磁波的水汽和云层，具备一定的地下穿透能力。凭借以上独特优势，SAR广泛应用于环境保护^[2-3]、目标检测^[4-5]、海洋监测^[6]、军事领域^[7]、地质和矿物质勘探^[8]等方面。

水体分割就是从SAR图像数据中将水体和陆地等区分开。通常，水体的均匀性被认为是表征水面的最佳特征，但SAR图像有着很强的斑点响应特性，其中，斑点噪声的存在也使得难以准确地提取水体^[4]，这阻碍了对图像的精准分析。常见的分割水体方法主要包括：灰度阈值法、区域生长方法和主动轮廓模型等。

基于灰度阈值的分割方法因其具备有效性及高效率被广泛应用于图像分割中，较常见的算法有大津法(Otsu)^[9]、自适应阈值法^[10]等。该方法原理简单，利用图像全局信息快速分割水体，但通常需要人工干预才能精准定位最优分割阈值，易被建筑物、树木、山体阴影及斑点噪声影响，导致获得的阈值与真实阈值相差甚远，并且若用此阈值分割图像会造成虚警率增大。

基于区域生长的方法^[11]常用于提取单块水体的范围，并不适用于复杂环境下水体的检测。主要原因是在复杂环境下，初始种子点的自动选取较为困难，一旦选择不恰当则会导致目标轮廓分割不完整或者过分割^[12]。

基于主动轮廓模型(active contour model, ACM)的方法在图像分割中也较为常见，它对初始位置有很强的依赖性，一般先经过阈值分割等传统方法进行初始分割，后续利用ACM修正分类结果。若是数据周围环境复杂，则对初分割边界精度要求更高。此外，由于ACM迭代耗时长且复杂度较高，需要多次尝试选择优化参数^[13]，所以效率相对较低。冷英等^[14]首先利用两级Otsu获取水体的初始轮廓线，基于区域ACM方法在轮廓线附近的窄带内迭代优化提取鄱阳湖边缘。但此方法更适用于独立的大面积湖泊提取，对于湖泊群效果则较差，边界易出现交叉现象。

鉴于上述传统方法存在诸多问题，本文提出一种基于改进的高斯混合模型(Gaussian mixture model, GMM)和图割模型(graph cut model, GCM)的水体图像提取算法。先是对原始SAR图像进行一定的预处理，包括下采样及滤波；其后采用两级Otsu初步分割图像，将基于输出结果计算得到的水体和非水体的统计数据作为GMM的初始参数，应用Xu等^[15]在2013年提出的最大期望算法(expectation-maximization algorithm, EM)迭代获取GMM的最优参数；最后，将获得的GMM中像素为水体和非水体的概率作为基于改进的图割模型中t-link边的权重，而n-link边的权重则是与像素邻域间的灰度和欧式距离有关，同时重新定义了图割模型的能量函数，通过优化能量函数实现湖泊分割。具体的流程如图 1所示。

1 图像预处理 1.1 下采样

通常，SAR图像数据尺寸较大，为提高分割速度，本文通过最近邻插值的方法对SAR图像数据进行下采样以等比例缩小图像尺寸，减小计算量的同时增强对水体信息的感受范围。对尺寸为I_M×I_N的SAR图像进行下采样，如下式所示

$ I_{\mathrm{m}}=\left(I_M / m\right) \times\left(I_N / m\right) . $

(1)

其中: m为下采样倍数。通常，针对大尺寸SAR图像采取下采样方法能够使水体特征更集中，提高运行效率。若是小尺寸图像则进行较小倍数的下采样，避免过度下采样导致有用信息大幅度减少。

1.2 高斯滤波

由于SAR图像受到相干斑点噪声的高度抑制，因此本文采用高斯滤波器^[16]降低相干斑噪声的影响，并将其与下采样相结合，使其在研究区域水体相对较少的条件下，直方图仍能呈现出明显的波峰和波谷^[17]。高斯窗函数^[16]被定义为

$ I_G(x, y)=\frac{1}{\sqrt{2 \pi} \sigma_x \sqrt{2 \pi} \sigma_y} \exp \left(-\left(\frac{x^2}{2 \sigma_x^2}+\frac{y^2}{2 \sigma_y^2}\right)\right) . $

(2)

其中: σ_x和σ_y分别控制窗的长和宽，这里σ_x=σ_y=3, (x, y)为图像I_m中像素坐标。

2 基于SAR图像的水体提取 2.1 GMM的建立

本文假设像元为水体和非水体的分布符合高斯分布，因此利用GMM拟合SAR图像直方图分布。GMM为概率分布模型，通常用所有类的高斯密度的加权叠加来表示。因此，模型可以表示为

$ P_i(x \mid \lambda)=\sum\limits_{i=1}^N w_i F\left(x \mid u_i, \sigma_i^2\right) . $

(3)

其中: P_i(x|λ)为概率密度，F为概率密度函数，x={x₁, x₂, …, x_M}为遥感图像中所有像元，M为SAR图像中所有像素数量。参数集λ={w_i, u_i, σ_i²; (i=1, 2)}, i=1为像元被标记为水体，i=2则被标记为非水体，其中，w_i为权重，u_i为均值，σ_i为标准差，N为分类数目，在本文中N=2。

像元属于每一类的概率密度函数如下所示

$ \begin{gathered} F\left(x \mid u_i, \sigma_i\right)= \\ \frac{1}{(2 \pi)^{\frac{1}{2}}\left|\sigma_i\right|} \exp \left\{\frac{\left(x-u_i\right)^T \Sigma^{-1}\left(x-u_i\right)}{2}\right\} . \end{gathered} $

(4)

从式(3)、式(4)可知，求解GMM的关键在于获得权重系数、方差、标准差的最优值，本文利用最大似然函数法来估计GMM的参数，似然函数表达式为

$ L(\theta)=\prod\limits_{j=1}^M \sum\limits_{i=1}^N w_i F\left(x_j \mid u_i, \sigma_i^2\right) . $

(5)

由于L(θ)的高度非线性，所以一般用EM算法迭代求解最佳参数。

2.2 EM算法优化模型参数

EM^{[15, 18]}算法获得GMM最优的参数集λ的步骤如下所示：

首先，输入$\lambda_0=\left\{w_{00}, u_{00}, \sigma_{00}^2 ; w_{11}, u_{11}, \sigma_{11}^2\right\}$，λ₀为初始参数集。本文采用基于Otsu的两级阈值分割算法^[9]获取λ₀。通常，水体分布相对集中，且集中在灰度值偏低的波峰区域，而这部分仍存在噪声以及植被等灰度值与水体特性相近的区域，因此，在对整体SAR图像进行初次Otsu阈值分割后，对灰度值和较小的那类所有像素再进行一次Otsu分割，以降低湿地、道路等部分对结果的干预，计算最终得到的二分类图像中水体和非水体的各自所占比重、均值以及方差，初步估计每一类的密度分布。

其次，应用EM算法使似然函数最大化，将L(θ)取对数为

$ \ln L(\theta)=\sum\limits_{i=1}^N \ln \sum\limits_{j=1}^M F\left(x_i ; \lambda\right) . $

(6)

整理后，基于式(6)分别对参数求偏导取极值，比重参数简化为

$ w_j^{\prime}=\frac{1}{N} \sum\limits_{i=1}^N P\left(I_i=j \mid x_i, \lambda\right) . $

(7)

均值为

$ u_j^{\prime}=\frac{\sum\limits_{j=1}^M P\left(I_i=j \mid x_i, \lambda\right) \cdot x_i}{\sum\limits_{j=1}^M P\left(I_i=j \mid x_i, \lambda\right)} $

(8)

方差为

$ \sigma_j^{2 \prime}=\frac{\sum\limits_{j=1}^M P\left(I_i=j \mid x_i, \lambda\right) \cdot x_i^2}{\sum\limits_{j=1}^M P\left(I_i=j \mid x_i, \lambda\right)}-\left(u_j^{\prime}\right)^2 . $

(9)

将初始参数λ₀代入式(7)~式(9)中，得到新的参数集合λ，不断重复上述步骤更新各参数，直到达到对数函数lnL(θ)收敛的终止条件，EM算法迭代结束，便得到最终参数结果λ={w_k, u_k, σ_k²; (k=1, 2)}。后续将依据确定的双高斯混合模型进一步设定GCM权重。

2.3 基于GCM提取水体

本文中的水体分割是通过将GCM中有向加权图的能量函数最小化实现的。有学者选择采用模拟退火、动态规划得到全局最优化解决方案，但是模拟退火方法受限于其较为低下的工作效率，动态规划则受限于巨大的空间复杂度^[19]。而GCM的高效、精确完全满足图像分割的需求^[20]，因此本文利用此模型对SAR图像数据进行分类，最终分为2类：前景和背景，分别表示水体和非水体^[21]。

GCM的有向图G^[22]如图 2所示，通常由一组节点V和一组连接它们的有向边集E^[19]构成，通过s/t的割(cut)将图中的所有像素节点划分为2个不相交的子集S和T，节点s在S中，节点t在T中，任何s/t割都将节点划分为不相交的组，而每个组仅包含一个终端。一般地，将从S到T的切割成本$C=\{\mathcal{S}, \mathcal{T}\}$定义为“边界”，也是边缘所有像元对(p, q)的成本总和，其中$p \in \mathcal{S}, q \in \mathcal{T}$，其汇总了有向边的权重。

边集E存在着2种边，邻域连接边n-link和终端连接边t-link，像素的8邻域像素(p, q)构成n-link，每个像素都与前景s、背景t相连接构成的边(s, p)和(p, t)为t-link。图的所有边都被赋予一定的权重，n-link的权重是源于像素之间的不连续性的损失，t-link的权重相当于将相应标签分配给像素的代价。

图中的每条边被赋予非负的权重，图的每个割成本等于能量函数中的所有项，因此要想实现水体和非水体的分割就要求取最小的成本^[22]。提出的算法对图割模型的能量函数进行改进，定义如下

$ E=\lambda E_{\text {region }}+(1-\lambda) E_{\text {image }} . $

(10)

其中: λ∈[0, 1]进一步明确区域项E_region和图像特征项E_image在能量函数中占据的比重关系。区域项能量函数E_region指式(4)求得的像元为前景和背景的概率，作为区分像元属性的依据，更多地关注潜在的目标像元，而能量函数E_image强调像元邻域强度及空间位置的关系。

定义的区域项能量函数E_region为

$ E_{\text {region }}=\sum\limits_{p \in P} D_p\left(l_p\right), $

(11)

$ D_p\left(l_p\right)= \begin{cases}-P_0, & \text { if } l_p=\mathrm{bg}, \\ -P_1, & \text { if } l_p=\mathrm{fg} .\end{cases} $

(12)

其中: P代表图像中所有像素，l_p则表示SAR图像中像元p被分配的标签，前景(fg)或背景(bg)，区域项能量函数为由式(3)确定的P₁(非水体)和P₂(水体)概率组成，D_p(l_p)则表示将标签分配给像元的惩罚。

除此之外，将图像特征项融入到能量函数中，图像特征项函数E_image^[23]公式如下

$ E_{\text {image }}=\sum\limits_{\{p, q\} \in N} V_{\{p, q\}}, $

(13)

$ V_{\{p, q\}}=\exp \left(-\frac{\left\|I_p-I_q\right\|^2}{2 \sigma^2}\right) \cdot \frac{1}{d}, $

(14)

$ \left\{\begin{array}{l} d=\sqrt{\left(x_p-x_q\right)^2+\left(y_p-y_q\right)^2} \\ \sigma^2=\mathbb{E}\left[\left\|I_p-I_q\right\|^2\right] \end{array} .\right. $

(15)

其中: I_p和I_q表示相邻像素p和q图像灰度值，d为像素p和q的欧氏距离，σ²表示整幅图像所有邻域像素对N的灰度差的均值，V_{{p, q}}则表示像元间相似度的惩罚项。若V_{{p, q}}较大则说明相邻像素间差异较小，很可能属于同一类别，相反，若V_{{p, q}}偏小，则说明像素间差异大，则像素更倾向于为水体和非水体的边界。

因此，按照图 2建立模型的有向图，包括相应的边集和节点，并依据能量函数E对边的权重进一步分析。表 1给出2种边n-link以及t-link的权重设定规则，边集n-link和t-link则按照GMM求取的像素属于不同端节点的概率密度设定。

3 实验与分析 3.1 实验数据

为验证提出算法的有效性，本文采用真实SAR图像数据对改进的GMM和GCM的分割能力进行评估，如图 3所示。图 3(a)、3(b)、3(c)分别是来自于哨兵(Sentinel-1)卫星的太湖、洞庭湖以及东湖数据，其中，太湖位于江苏省南部^[24]，30°55′40″~31°32′58″N和119°52′32″~120°36′10″E之间。洞庭湖位于湖南省内长江右岸，以湘江、紫江、沅江、溧水河为南、西支流；太平口和大池口作为长江水流的引流口，湖水在城陵矶附近的东北出水口排回长江^[25-26]。东湖同样坐落于湖南省，位于岳阳市华容县。

所采用的数据均源自哨兵1号卫星共有的2颗雷达卫星，Sentinel-1A和Sentinel-1B，分别由欧洲航天局于2014和2016年发射^[27]。哨兵1号卫星作为C波段雷达成像系统，其单卫星在12 d的重复周期内可扫描全球一次，2颗卫星仅需6 d就可完成一次观测^[28]，并且可提供单极化和双极化图像。

本文所采用的3景实验数据为VV+VH极化IW模式的Sentinel-1 GRD数据，数据信息如表 2所示。

3.2 实验结果

为说明所采用的二级Otsu阈值算法真实有效，本文基于图 3(a)数据将其与仅进行一次阈值分割算法进行对比，结果如图 4所示。

从图 4(b)可以明显看出，采用二级阈值较为完整地勾勒了水体的轮廓，较好地排除了山体、建筑、阴影等的干扰。

同时，本文选用GMM拟合SAR图像中像素灰度直方图的分布。图 5是太湖的灰度直方图及其拟合的双高斯混合模型，可以很明显地看出预处理后SAR图像的直方图呈现双波峰，波谷清晰可见。采用多级Otsu阈值和EM算法拟合的GMM趋势近似等同于图 5(a)，可见假设成立。但是，由于GMM未考虑到像素邻域之间的关系，本文后续基于3景Sentinel-1数据利用GCM实现湖泊的自动化识别。所提出算法输入的是SAR图像数据的强度图，输出为黑白二值图，白色为识别的水体，黑色为非水体。

为检验提出算法的准确性、鲁棒性，利用Otsu阈值^[9]、马尔科夫随机场(Markov random field, MRF)^[29]方法以及文献[30]提出的改进图割算法(improved graph cut, IGC)对实验数据进行水体提取，与提出的算法进行对比，对比结果如图 6所示。为保证结果真实可靠，对比算法采取了相同的预处理。

3.3 对比分析

图 6每一列分别给出真值图、提出的方法、Otsu算法、MRF算法以及IGC算法分割结果。本文分别从定性和定量分析的角度比较不同算法之间的差异。通过目视解译很明显地发现，3种方法在视觉效果上存在明显的差别，本文提出的方法对于湖泊和细小河流的检测更加准确，而其余3种对比方法在环境复杂区域识别效果较差。如此大的差异可能是由于水体与周围环境灰度相近，对比方法却无法正确区分导致的。

为进一步探索4种方法检测结果的差异，本文根据6个评价指标定量地分析算法精度，其中包括准确度(O_A)、精准度(precision, Pre)，召回率(recall, Rec)、Kappa系数、F₁分数(F₁-score)以及算法的运行时间，各项指标定义及公式如下所示。

O_A表示正确分类为水体与像素总数的比值

$ O_{\mathrm{A}}=\frac{T_{\mathrm{p}}+T_{\mathrm{N}}}{\text { Pre }} \times 100 \% . $

(16)

Pre为正确分为水体与预测为水体的总数的比值

$ \operatorname{Pre}=\frac{T_{\mathrm{P}}}{T_{\mathrm{P}}+F_{\mathrm{P}}} \times 100 \% . $

(17)

Rec为正确分为水体与参考图像中被标记为水体数目的比值

$ \text { Rec }=\frac{T_{\mathrm{P}}}{T_{\mathrm{P}}+F_{\mathrm{N}}} \times 100 \% . $

(18)

Kappa系数用于衡量分类的效果一致性

$ \text { Карра }=\frac{\left(O_{\mathrm{A}}-P_{\mathrm{e}}\right)}{\left(1-P_{\mathrm{c}}\right)} \times 100 \% . $

(19)

F₁-score为Pre和Rec的调和平均数

$ F_1 \text {-score }=2 \times \frac{\text { Pre } \times \text { Rec }}{\text { Pre }+ \text { Rec }} $

(20)

$ \begin{gathered} P_{\mathrm{e}}=\left(\left(T_{\mathrm{P}}+F_{\mathrm{P}}\right) \times\left(T_{\mathrm{P}}+F_{\mathrm{N}}\right)+\right. \left.\left(F_{\mathrm{N}}+T_{\mathrm{N}}\right) \times\left(F_{\mathrm{P}}+T_{\mathrm{N}}\right)\right) / T^2 . \end{gathered} $

(21)

$ P=T_{\mathrm{P}}+F_{\mathrm{P}}+T_{\mathrm{N}}+F_{\mathrm{N}} . $

(22)

其中：T_P为成功预测为水体的像元数；F_P为将非水体预测为水体的数量；T_N为成功预测为非水体的像元数；F_N为将水体错检为非水体的数量。

表 3反映了包含上述评价指标以及运行效率的所有算法定量比较结果。

对于3景图像的结果来说，传统的Otsu以及MRF与提出的方法相比性能差距较大，将大量的建筑、山体等区域误判为水体，导致整体评价指标的结果均逊色于提出的算法。特别是对于洞庭湖的识别结果，提出的方法的提取精度更是远超过对比方法。并且，从表 3中可以看出，提出的方法拥有较高的总体精度(99.61 %、98.60 % 和99.27 %)，表明其在减少错误预测上有所改善，与真值图具有良好的一致性。定性比较本方法和同样采用GMM及GCM的IGC算法，均较为完整地检测到存在的水体；但定量评价指标表明，IGC算法仍存在一定的误判。同时，由于对比方法随机设置初始值，导致需要大量迭代，效率弱于对比方法。此外，本方法相较于IGC提取水体轮廓更加全面，有效降低了错分率，由此可见提出的算法能够更好地捕捉到太湖、洞庭湖以及东湖的边界轮廓，充分突出了本方法在湖泊提取的优越性和有效性。虽然由于EM算法的循环迭代，导致本算法的运行效率弱于Otsu方法，但各项评价指标的结果都明显优于Otsu方法。

3.4 算法复杂度分析

对于本算法，时间复杂度主要集中在二级Otsu阈值算法、EM算法以及有向加权图的权重计算中。其中，Otsu阈值算法时间复杂度为O(n)^[31]，而本文所采用的算法在第二次分类中仅处理第一次Otsu阈值被标记为水体的像素，假定像素数目为m(m < n)，时间复杂度则为O(n)+O(m)。

EM算法中假定将全部像素迭代N次，则其时间复杂度则为O(n×N)。太湖、洞庭湖以及东湖迭代次数分别为2、3、2。

GCM中主要是中心像素与周围8邻域像素关系的权重计算较为耗时，因此算法中这部分的时间复杂度则为O(8n)。

最终，根据加法原则，算法整体复杂度为O(n)+O(m)+O(n×N) +O(8n)= O((9+N)n+m)。上述n均表示参与算法的像素总数目，通常，在分析算法时间复杂度时，可忽略与最高次项相乘的常数，所以算法的总时间复杂度化简为O(n+m)。

4 结论

传统的SAR图像分割算法针对周围环境较为复杂的湖泊存在检测效果较差、边界模糊等问题，因此，本文提出一种基于改进的GMM和GCM的提取方法。首先对SAR图像进行预处理，通过下采样使水体更加集中的同时减小数据运算量，联合高斯滤波去除相干斑噪声。然后，通过两级Otsu阈值分割初步获取GMM的初始参数，从而大幅度地减少EM算法优化参数的迭代次数，提高运行效率。最后，基于像元属于水体和非水体的高斯概率设定GCM有向图的边集权重，进而最小化改进的能量函数实现图像分割。实验结果表明，与传统Otsu分割方法、MRF方法、IGC算法相比，提出的方法更加完整地解译出对象信息，提取湖泊结果更优，精度更高，具有一定的应用价值。该算法在各类湖泊提取中适用性较好，后续研究中，将引入深度学习技术，进一步提高算法在水体识别中的泛化能力。