一种新颖的亚像素级SAR图像水陆分割方法

地表水是地球水资源的重要组成部分，在地球的水文和生物化学循环中具有举足轻重的作用。利用遥感图像实现精确的水陆分割是开展基于水体各项研究的工作基础，有着重要的现实意义。合成孔径雷达(Synthetic Aperture Radar, SAR)是一种主动微波遥感技术，可穿透云层、降雨、雾霾等实现全天时、全天候对地观测，且具有响应速度快、覆盖范围广、周期性观测等特点，目前已广泛应用于大尺度地表调查与研究中[1]。

SAR图像水陆分割研究发展至今，经典的阈值、聚类、主动轮廓模型(Active Contour Model，ACM)等方法，以及新兴的深度学习算法都有较多的研究与应用[2-6]。这些方法虽然能够有效地实现SAR图像水陆分割，但分割尺度多停留在像素级别。对于星载SAR系统而言，由于受最小天线孔径积的限制，常见工作模式的图像空间分辨率通常在几米至几十米量级[7]。例如Sentinel-1卫星TOPS模式下的空间分辨率为20 m，国产高分三号(GF-3)卫星标准条带模式下的空间分辨率为8 m。这些中低分辨率的图像中，水域边界上会存在大量既包含陆地又包含水域的混合象元，传统算法无法在混合象元内确定真实水域轮廓位置，只能将整个象元划分为陆地或者水域，极大影响了水陆分割结果的精度。同时，由于SAR图像固有的相干斑噪声在水域边界附近分布杂乱，会破坏真实的水域边界，这也是导致水陆分割精度难以提高的原因之一。而SAR图像水陆分割中几个像素的误差对应于实际地物往往是数十米甚至上百米的差异。因此，设计一套性能良好且可实现亚像素级水陆分割的方法，对于更加有效地利用中低分辨率SAR图像具有重要意义。

针对上述问题，本文提出了一种新颖的亚像素级SAR图像水陆分割方法，该方法以改进非局域相干斑滤波算法与模糊C均值(Fuzzy C-Mean,FCM)聚类算法完成像素级水陆粗分割，在此基础上根据水域特征，消除FCM聚类粗分割结果中的噪声、阴影等干扰，利用双三次样条插值算法对获得的像素级水域轮廓邻域进行上采样，并在该邻域内应用几何主动轮廓(Geometric Active Contour, GAC)模型获取亚像素级水陆分割位置点，最终通过分段曲线拟合方法获得完整的亚像素级水域轮廓。为了验证本文所提方法的有效性，应用国产GF-3卫星的多模式SAR图像进行验证实验，结果表明，与传统方法相比，本文所提方法平均像素偏移精度提升了一个数量级，具有良好的水陆分割效果。

1 亚像素级水陆分割方法

本文提出的亚像素级SAR图像水陆分割方法流程如图1所示。该方法的第一步，使用改进的非局域滤波算法抑制图像中的相干斑噪声；第二步，应用FCM算法对水陆边界进行像素级的粗分割；第三步，结合双三次样条插值和GAC模型，在像素级水域轮廓的基础上实现亚像素级水陆分割，通过以上步骤可以获取到更精确的亚像素级水域轮廓。

1.1 非局域相干斑滤波算法

相干斑是一种SAR图像特有的乘性噪声，当电磁波照射粗糙表面时，电磁波到达的相位不同，导致回波产生干涉现象，使得SAR图像中分布着大量由干涉引起的斑点噪声，这些噪声会使SAR图像处理的难度大大增加。近年来，在抑制SAR图像相干斑噪声方面，非局域滤波算法取得了较好的效果。

非局域滤波算法中采用像素相关性反映图像相邻位置像素值的相关程度。文献[8]证明了比值距离(Ratio Distance, RD)是一种鲁棒性强的特征，可以用来描述SAR图像中像素之间的相关性。本节将介绍一种基于RD的非局域SAR图像滤波算法。

SAR图像中每个像素的测量值Y，真实值X以及噪声值Z三者的关系可以被描述为Y=X×Z。对于一景L视SAR幅度图像，测量值通常用Goodman模型中的幅度模型描述，该模型又被称为均方根伽马(Nagakami-Rayleigh)分布[8]，如式(1)所示。

式中，A为每个像素的观测值，σ为待估计的真实强度值，则真实的幅度值X为σ的均方根。

非局域滤波算法中，定义以每个像素为中心的小正方形窗为一个相似块。假设两个相似块的真实幅度值相同，观测值的RD将从YNi/YNj转变为ZNi/ZNj。将新的RD带入到块相似性测量方程S(i,j)以及两个相似块的联合概率密度函数p(i,j)可得式(2)和式(3)。

式中：j是非局域窗中的像素，非局域窗是以像素i为中心、边长为r1的正方形窗；G为标准高斯核函数；ZNk是相似块Nk的观测值之和，Nk是以像素k为中心、边长为r2的正方形相似块；p(ZNi/ZNj)是ZNi/ZNj的概率密度函数。

滤波后的SAR图像I可由新定义的S(i,j)计算得到，如式(4)所示。

1.2 基于FCM算法的像素级水陆分割方法

FCM是一种无监督的聚类算法。一幅图像可描述为集合I={i1,i2,…,iN}，传统FCM算法通过设置c个不同的聚类中心，计算集合I中每个像素与各聚类中心的隶属度，完成对图像的分类。FCM算法执行时，需遍历图像中每一个像素计算隶属度，当处理大尺度SAR图像时，非常耗时。针对该问题，文献[4]提出了一种基于灰度分布的FCM算法，采用如式(5)所示的目标函数Jm，将与聚类中心隶属度的计算由像素转变为灰度级，在灰度级上完成对图像的分类，从而大大提高了算法的运行效率。

式中：Q是图像I的灰度级数量，

是具有相同灰度值i的像素总数; vj表示第j个聚类中心;m为模糊指数;uji表示第i个灰度值与第j个聚类中心的隶属度，隶属度求解方程如式(6)和式(7)所示。

本文用于提取像素级水陆分割的FCM算法执行步骤如下：

(1)设置聚类中心数c=4，并随机设置聚类中心数值；设置模糊指数m=2；设置迭代次数m1=15。

(2)根据式(6)计算样本隶属度，即SAR图像中各灰度级到各聚类中心的隶属度。

(3)根据式(7)更新各聚类中心。

(4)判断迭代是否完成，未完成则循环执行步骤(2),步骤(3)。

(5)迭代完成后，若‖Ii-vj‖2<‖Ii-vk‖2,(j,k)∈[1,c]，则将灰度值等于i的像素划分为第j类。

1.3 亚像素级水陆分割方法

本节提出了一种新颖的亚像素级水陆分割方法，该方案以像素级水陆分割结果为基础，利用双三次样条插值算法提高其周边区域的清晰度，并结合GAC模型求解梯度最大位置来确定亚像素级水陆边界位置，方法流程如图1中步骤3所示。

该方法是在以像素级水域边界线的像素点为中心，半径为r3的邻域内进行水陆精细分割的，其有效性建立在两个假设之上：第一，邻域内需同时包含水域与陆地，如图2(a)所示；第二，当雷达后向散射系数梯度值最大时，邻域内的陆地与水域产生分界[9]，如图2(b)所示。

亚像素级水陆分割过程包括以下5个步骤：

(1)自然水域周边地形一般较为复杂，陡峭山脉等地物会因起伏地形未受到电磁波照射，没有回波，会产生较多的阴影现象，在SAR图像中呈现为黑色区域。而水面对电磁波产生表面散射，后向散射系数较低，平静水域在SAR图像中也同样呈现为与阴影相似的黑色区域。为实现对二者的区分，需进一步提取感兴趣区域(Region of Interest, ROI)。在粗分割获取的像素级水域中，记最大连通区域面积为Smax，设置面积阈值Ta=μ·Smax。当连通区域面积大于Ta时，被记为ROI，否则丢弃。同时，将提取到的ROI区域边界信息记为像素级水域轮廓。本文根据经验将μ设置为0.2。

(2)以步长k遍历像素级水域轮廓，以访问到的每一个像素为中心，建立边长为r3的正方形亚像素处理窗(Subpixel Processing Window，SPW)，记录该SPW中心坐标为(x1,y1)，计算得到第一个像素坐标为(x2,y2)=(x1-r3,y1-r3)，如图2(a)中虚线框所示。在一个SPW中，以第一个像素为坐标原点，建立坐标系，则该SPW中所包含的像素级水域轮廓点可通过多项式拟合，表示为曲线方程：fw l(x,y)。方程阶数取决于SPW的大小，SPW越大，包含的像素级水域轮廓点越多，多项式阶数越高。通常，多项式阶数过高可能会导致最终结果不具备鲁棒性，相反小的SPW则会对高阶多项式要求较小。

(3)利用双三次样条插值方法，对SPW内图像进行λ倍的上采样，增加局部像素数目，减小后续在亚像素尺度下调整水陆边界时对最终边界的影响，如图2(c)与图2(d)所示。

(4)插值后的SPW中，梯度最大的位置表示水域轮廓，再利用GAC模型对水域轮廓进行精细化处理，初始轮廓由曲线方程fw l(x,y)获得。如图2(d)中所示，红线位置表示梯度最大位置，同时也表示亚像素级水域轮廓所在位置，由此可得到亚像素水域轮廓点的坐标集合

(5)根据步骤(1)与步骤(3)中记录的坐标，还原SPW中获取到的亚像素级水域轮廓到原始图像中，得到新的亚像素级水域轮廓点的坐标集合，可表示为

值得注意的是，通过像素级水域轮廓计算得到的新的水域轮廓可能会存在多个位置，如图3中黄点所示。造成这种情况的原因是距离较近的像素点产生的SPW会有部分重叠，从而使重叠区域被多次处理，产生多个亚像素水域轮廓位置。通过分段曲线拟合方法可以解决该问题，由此拟合得到的一条连续曲线即为亚像素级水域轮廓，结果如图3中红线所示。

1.4 基于GAC模型的亚像素级水陆分割线获取方法

GAC模型是基于Snake模型和水平集函数(Level-Set Function，LSF)发展而来的，它将能量最小化问题转化为在黎曼空间中求最小曲线路径的问题[5]。GAC模型对应的LSF可以表示为

式中，∂φ为初始LSF， width=17,height=17,dpi=110

为梯度算子，g(x)为基于梯度信息的边缘停止函数，α为控制轮廓膨胀或收缩的气球力参数。GAC算法的执行步骤如下：

(1)按如式(9)所示方法初始化LSF。

式中，ψ为SPW内的图像，ψ0为SPW中的水域区域，σψ0为由SPW中拟合的fwl函数获得的曲线，ρ为常数参数，本文设置ρ=1。

(2)设置气球力参数α，并根据式(8)演化LSF。本文令α=5。

(3)使用高斯滤波器调整LSF，如φ=φ×Gτ。τ为高斯函数的标准差，范围为0.8到1.5。本文设置高斯核的边长h=5，标准差τ=1。

(4)设置最大演化次数m2=30，若未达到演化次数，返回步骤(2)。

2 实验结果与分析

2.1 研究区域及数据

实验数据采用国产GF-3号卫星在丹江口水库主库区观测得到的两景SAR图像，范围如图4中黑色框所示，图像详细信息如表1所示。国产GF-3号卫星是一颗C波段全极化SAR卫星，适用于多种工作场景，包括海洋监测、灾害监测等[10]。

2.2 实验结果

除本文上述步骤中已给出的实验参数外，本文方法其他的实验参数如下：(1)非局域相干斑滤波步骤中，非局域窗边长r1=15，相似块边长r2=7，视数L=1。(2)亚像素水陆分割线提取方案中，步长k=2，SPW边长r3=7，上采样倍数λ=5。

图5展示了非局域滤波算法对SAR图像相干斑噪声抑制的结果。两景原始SAR图像如图5(a)和图5(f)所示。为了更好地展示细节，选取区域1至区域4放大，原始图像如图5(b)、(d)、(g)和(i)所示，滤波后图像如图5(c)、(e)、(h)和(j)所示。由图示结果可知，该滤波方法具有的良好平滑性能和边缘保持效果，陆地与水域上的斑点噪声被平滑，陆地呈均匀的灰色，水域呈均匀的黑色，水域边界也得到了保持和加强，整体上提高了SAR图像的均匀性，为后续处理创造了良好的条件。

图6展示了像素级水陆分割线提取的结果和亚像素级水陆分割线提取的结果。为了更为清晰地观察，选择区域A到D放大，分割结果如图6(b)、(c)、(e)和(f)所示。红线代表亚像素级水域轮廓线，蓝线代表像素级水域轮廓线。两景SAR图像的标称分辨率分别为5 m和10 m，因此两景SAR图像中的一个像素分别代表25 m2和100 m2。结果的目视分析可知，亚像素级水域轮廓(红线)比像素级(蓝线)更加平滑，更加接近真实水陆边界。因此，在水域SAR遥感监测中，亚像素级水陆分割具有更好的性能与应用前景。

2.3 对比实验与定量分析

为了验证本文所提方法的有效性，选取目前主流，且精度较高的马尔可夫随机场(Markov Random Field，MRF)方法[3]、FCM方法[4]和ACM方法[5]这三种SAR图像水陆分割方法开展对比实验，分割结果如图7所示。同时，选择虚警率(False Alarm Rate)、准确率(ACC)以及平均轮廓偏移距离(Average Contour Offset Distance, ACOD)等量化指标分别评估整体分割效果与局部轮廓质量[3]。参考水域为人工追踪法[2]获取的水域轮廓。评价指标计算方法如式(10)、式(11)和式(12)所示。

(3)平均轮廓偏移距离

式中，AE为算法检测的水域面积，AR为人工检测的参考水域面积，AE&R为AE与AR的交集面积，N是提取水域轮廓上的像素数，d(EinRBi;R)是像素i与实际水域轮廓之间的最短距离。

目视观察各方法的实验结果，本文所提方法与MRF和FCM方法相比，对于水域的提取更为精确，受SAR图像阴影、相干斑噪声以及水面复杂散射影响较少。ACM算法在大尺度SAR图像中演化较慢，当达到最大演化次数时，仍未检测出水库边缘的细小分支。该结果也可由如表2所示的定量分析结果验证，本文方法对两组实验数据的提取准确率均大于99%，虚警率均小于0.5%，并且轮廓局部质量优异，轮廓平均偏移距离均小于0.2个像素，取得了更高精度的结果，优于其他方法。

3 结束语

本文针对中低分辨率SAR图像精确水陆分割的难题，提出了一种新颖的亚像素级SAR图像水陆分割方法。该方法利用改进的非局域滤波算法对SAR图像的相干斑噪声进行抑制；采用FCM算法完成像素级水陆粗分割；并在粗分割基础上，应用基于双三次样条插值算法与GAC模型的方法完成亚像素级的水陆精分割。通过对丹江口水库区域的两组GF-3号卫星多模式SAR图像的实验验证，与目前常见SAR图像水陆分割算法相比，本文方法取得了更高的准确率以及更低的虚警率，平均像素偏移均小于0.2个像素，最小达到0.13个像素，精度提升了一个数量级。同时，该方法具有较强的鲁棒性，可获取更为精确的水域轮廓，对利用星载SAR图像检测水库、湖泊等水域变化有较好的工程使用价值。

[1]牛世林, 郭拯危, 李宁，等.星载SAR水域分割研究进展与趋势分析[J].聊城大学学报(自然科学版), 2018, 31(2):72-86.

[2]LI Ning, NIU Shilin, GUO Zhengwei, et al.Dynamic Waterline Mapping of Inland Great Lakes Using Time-Series SAR Data from GF-3 and S-1A Satellites: A Case Study of DJK Reservoir, China[J].IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing, 2019, 12(11):4297-4314.

[3]LIU Z , LI F , LI N , et al.A Novel Region-Merging Approach for Coastline Extraction from Sentinel-1A IW Mode SAR Imagery[J].IEEE Geoscience and Remote Sensing Letters, 2016, 13(3):324-328.

[4]郭拯危, 王乐, 宋国磊.基于混合模糊的SAR图像水陆分割算法[J].国土资源遥感, 2018, 30(4):62-67.

[5]CHEN Jiaqi, WANG Qingwei, WANG Jian, et al.Change Detection of Water Index in Danjiangkou Reservoir Using Mixed Log-Normal Distribution Based Active Contour Model[J].IEEE Access, 2019, 7: 95430-95442.

[6]GUO Yaru,ZHANG Jixian.A New 2D OTSU for Water Extraction from SAR Image[J].ISPRS-International Archives of the Photogrammetry, Remote Sensing and Spatial Information Sciences, 2017: 733-736.

[7]李宁, 牛世林.基于局部超分辨重建的高精度SAR图像水域分割方法[J].雷达学报, 2020, 9(1): 174-184.

[8]FENG Hongxiao , HOU Biao , GONG Maoguo.SAR Image Despeckling Based on Local Homogeneous-Region Segmentation by Using Pixel-Relativity Measurement[J].IEEE Trans on Geoscience & Remote Sensing, 2011, 49(7):2724-2737.

[9]PARDO-PASCUAL J E , ALMONACID-CABALLER J, RUIZ L A, et al.Automatic Extraction of Shorelines from Landsat TM and ETM+ Multi-Temporal Images with Subpixel Precision[J].Remote Sensing of Environment, 2012, 123(6):1-11.

[10]SUN Jili, YU Weidong, DENG Yunkai.The SAR Payload Design and Performance for the GF-3 Mission[J].Sensors, 2017, 17(10):2419-2431.