水体在地球陆地面积中占据了较大的比例，在生态环境保护、旱涝灾害预防等方面发挥着重要作用。快速大面积地获取水体分布信息对于国土资源管理、用水安全保障、灾害快速评估等方面具有重要意义(Lv 等，2010；Yang 等，2001；王小标 等，2018；沈占锋 等，2013)。光学影像中的水体呈现出不同的光谱特征，增大了基于光学影像的水体分布制图难度。合成孔径雷达(SAR)影像以对云雨雾的不敏感性，反映了地物的后向散射信息，这种性质使得以镜面散射为主的水体呈现出不同于其他地物的后向散射能力，有效弥补了光学影像的局限性(Wang 等，2010)。高分三号(GF-3)卫星的成功发射，大大缓解了国内SAR数据源的长期缺乏问题，为了尽快将GF-3多极化、全极化数据服务于国家水资源宏观管理业务，并为未来国内多频率、多极化SAR遥感技术的快速发展提供技术支撑，研究一种高自动化、高精度的多极化SAR水体信息提取方法具有重要的应用价值。

SAR图像后向散射强度的大小与地物的表面粗糙度息息相关，表面光滑的水体多产生镜面散射，对雷达信号的后向散射能力较弱(Mason 等，2014)，因此基于阈值分割的SAR图像水体信息提取方法发展较早。对单极化影像中水体和非水体类别的后向散射强度进行直方图统计分析，获得类别分割阈值的方法简单高效(孙亚勇 等，2014；安成锦 等，2010；李智慧 等，2013)，已发展了较多的阈值自动、半自动选择方法。Hong等(2015)利用Landsat影像进行ISODATA分类，从而对Radarsat-1强度影像中的水体类别和其他类别进行统计分析，根据概率密度函数相等原则选择分割阈值，并结合地面高程信息获得了较高精度的洪水分布图，该方法要求在研究区有高质量的光学影像覆盖，且分割阈值不具有唯一性。邓滢等(2016)综合利用了全极化影像的Freeman极化特征和多尺度分割特征，使水体提取精度优于单一方法，但是需要不断尝试分割尺度和合并尺度。Huang等(2018)以地表覆盖类型产品为样本数据，利用随机森林的方法提取哨兵影像中的水体信息，采用随机选点的方式进行精度验证，实现了较高的自动化程度。以上研究多根据图像的统计特征分析地物的后向散射强度，以影像目视解译结果或已有的水体—非水体分类产品为先验知识，主要针对地表覆盖类型识别及水域面积提取，要求有一定的人机交互操作，并不适合于紧急灾害下的水体信息快速提取。

为了更精细高效地进行水域资源的监测管理，对水体信息的提取精度和快速响应能力要求越来越高，国内外学者在水体与非水体的阈值自动化选择方面进行了大量研究，多采用的是国外发展较早、图像质量较好的星载SAR影像。陈玲艳等(2014)利用水体的GAMMA分布和背景值的高斯分布特点，求得最佳阈值分割Radarsat-2影像，得到的水体提取结果优于传统的Otsu和纹理信息提取法，该方法仅需优化图像的分块数量，最大程度地降低了人工参与度。Xie等(2016)首先对全极化Radarsat-2进行wishart监督分类，然后基于河流、湖泊、池塘的不同形态特征，利用SVM、CART、TB、RF(随机森林)等分类器进行水体类型的识别，RF分类器获得了95%的水体类型分类精度。透视收缩和阴影现象的存在降低了SAR图像在地形复杂地区的成像质量(Sohn 等，2005)，目前多采用基于R-D定位模型的正射纠正方法，在发生透视收缩和山体阴影等像素失真区域采用了图像重采样、空间插值等方法弥补缺失数据，这给图像统计分析带来一定误差。国内外研究多采用地面高程信息和SAR轨道参数计算图像失真区域，从而在水体粗提取结果中剔除，得到精度更高的水体精提取图(Musa 等，2015；Hahmann 等，2008；庞科臣 等，2016)，具有计算量少，可操作性强的优势。

以上文献表明，基于SAR影像的水体信息提取研究需要光学影像或者SAR图像的目视解译结果作为先验知识，以阈值分割或监督分类的方法得到水体—非水体分布图，该类方法需要目视解译等人机交互过程，无法实现较高的自动化。国外中、高分辨率星载SAR传感器对我国的拍摄时间及拍摄地区有限，尽管哨兵1卫星可实现大范围、短周期地对地监测，但是较低的空间分辨率难以满足小区域水体识别。随着中国GF-3数据的不断丰富，研究一种基于GF-3影像的高自动化、高精度的水体信息提取方法对国产遥感影像的快速业务化应用具有重要推动作用。

实验区位于湖南省东北部(图1)，经纬度范围是26°47′—29°6′ N，111°51′—113°27′ E，覆盖了长沙市、株洲市、湘潭市、益阳市等地级市。海拔高度110—150 m，亚热带季风气候，雨量充沛，年均降水量达到1361 mm，内部有湘江及众多支流穿过，有较多的山地分布，山间有较小的湖泊和水洼，北部有东烂泥湖、横岭湖两个面积较大的湖泊。雨季的大量降雨使得当地受到严重的洪涝灾害。本文研究的目标是水体分布信息的提取，包括河流、湖泊、水库、坑塘、内陆滩涂、沟渠、水工建设用地等。

图 1 研究区位置 Figure 1 The location of study area

采用GF-3 L1级数据，包括1景精细成像模式(FSII)的HH、HV双极化影像和1景全极化模式(QPSI)影像，成像时间分别是2017年7月9日、2017年5月24日。FSII影像的距离向分辨率为1.1 m，方位向分辨率是5.8 m，入射角为27.8°，南北长156 km，东西宽146 km。QPSI影像的距离向分辨率为4.5 m，方位向分辨率是4.7 m，入射角为49.4°，南北长40 km，东西宽34 km。采用ASTER DEM作为辅助数据。本文采用两种成像模式的GF-3影像覆盖范围(图2)，底图是HV极化强度影像。

图 2 GF-3影像覆盖范围 Figure 2 The coverage of GF-3 image

为定量地评价本研究所发展的水体信息提取方法的精度，获取了实验区内与GF-3成像时间接近的Landsat 8 OLI影像和高分一号(GF-1)WFV模式影像，经过正射校正和融合处理后，随机选择检验样点，目视解译样点的水体—非水体类别属性，在光学影像难以辨识的地区还会结合GF-3假彩色合成影像和Google Earth影像进行仔细识别，作为后期实验结果的精度评价参考数据。比较样点位置在实验结果和参考影像上的水体—非水体类别信息，若类别一致则定义为正确的检测样点，不一致定义为错误的检测样点。在QPSI模式影像上随机选择了100个检验样点，FSII模式影像的覆盖范围较大，随机选择了200个检验样点。图3中黄色多边形分别是QPSI模式和FSII模式的GF-3影像在Landsat 8影像和GF-1影像上的覆盖范围。

图 3 水体提取参考图(R: NIR, G: R, B: G) Figure 3 The reference map of water extraction

采用正确率指标(TDA)，即检测到的正确水体样点和正确的非水体样点总数量占总样点数量的百分比，定量评价水体信息提取的精度。

本文首先对GF-3数据进行预处理(图4)，转换到地理坐标系下。由于在山体阴影产生了辐射失真，严重影响了图像的概率分布特征，已有的研究是在获得水体信息分布图以后，利用DEM和SAR轨道参数生成的阴影掩模文件掩模水体分布图，并不能消除失真像元对未失真像元的概率分布影响。本研究首先利用阴影掩模文件去除SAR图像的辐射失真区域，然后利用Otsu、KI二值化方法进行单极化阈值分割，得到初始的水体分布图。Fisher变换和基于高斯分布的最大期望—马尔可夫随机场(GD-MRF)算法以阈值分割得到的初始水体分布图为先验知识，综合了多极化图像的特征，通过不断的循环迭代，获得水体精提取分布图。最后以成像时间接近的光学影像为参考，随机选取样点，通过目视解译的方法定性和定量地评价水体信息提取的精度。

图 4 水体提取技术流程 Figure 4 Technique flow chart of water extraction

GF-3原始数据采用两个波段将实部与虚部存储在一个文件中，本文按照式(1)，从元数据中获取定标参数 $Qu$ 和 $Ci$ ，分别对实部和虚部进行辐射定标，然后组合为单视复数(SLC)数据。结合雷达入射角和距离向、方位向的几何关系，确定FSII模式采用距离向5视、方位向2视，QPSI模式采用距离向3视、方位向2视进行多视化处理。采用空间分辨率为30 m的ASTER DEM数据将GF-3影像从雷达坐标系转换到地理坐标系，并经过Lee滤波和取对数运算，生成空间分辨率为15 m的SAR正射影像产品和阴影掩模文件Mask，完成自动化预处理过程。

后向散射能力的大小与地物的表面粗糙度及影像入射角大小息息相关。陆地水域以镜面散射为主，后向散射能力很弱，而植被、城镇等非水体表面粗糙，对雷达波束具有较强的后向散射能力，因此可采用一定的阈值分割方法，当图像的后向散射强度小于阈值时定义为水体，大于阈值定义为非水体，从而实现水体信息的自动提取。为寻找最佳的分割阈值和极化方式，本研究在FSII模式图像上分别选择了水体和非水体的感兴趣区ROI，其中水体ROI包括了河流、湖泊、坑塘等地物类型，共173113个像元，非水体ROI包括了耕地、城镇、森林、裸地等地物类型，共184346个像元。分别统计FSII模式下HH和HV极化的水体非水体两个类别的像元后向散射强度范围，构建HH、HV两幅影像的统计直方图(图5)。

图 5 FSII模式下HH、HV极化的水体—非水体直方图 Figure 5 The water-nowater histogram of HH, HV polarization under FSII mode

在QPSI模式的影像中选择水体和非水体的感兴趣区，分别包含了15031个水体像元和65423个非水体像元，构建HH、HV、VH、VV 4幅影像的统计直方图(图6)。

图 6 QPSI模式下HH、HV、VH、VV极化的水体非水体直方图 Figure 6 The water-nowater histogram of HH, HV, VH, VV polarization under QPSI mode

从以上水体—非水体的统计直方图中可以看出，水体的后向散射强度明显小于非水体，水体—非水体类别呈现出两个波峰分布的特点，具有较好的可分性，交叉极化(HV，VH)和同极化(HH、VV)数据相比，水体与非水体的直方图重叠区域更小，在水体与非水体的分类方面具有更强的可分离性，原因在于交叉极化图像信息主要由体散射所引起，对镜面反射的敏感程度较低，水体表面光滑，均一性较强，在交叉极化图像上产生的噪声水平更低，具有较小的类内方差，更有利于水体信息的识别。通过对交叉极化概率密度曲线的目视解译发现，FSII模式HV图像以小于–27 dB为水体，QPSI模式HV图像以小于–22.5 dB为水体，可获得最好的水体—非水体的分类效果。这是一种以目视解译为主要手段的半自动化水体提取方法，该方法类似于监督分类，通过对训练样本的后向散射强度分析，总结图像在不同地类上的分布特征，从而较准确地提取出水体目标，已被许多学者广泛使用。为提高阈值选择的自动化程度，本文利用Otsu、KI二值化的方法进行图像分割，并将自动选择的阈值和目视解译的阈值对比分析，通过两者之间的接近程度判断自动阈值分割法的有效性。

Otsu方法是一种发展较早的图像二值化方法，该方法采用图像中灰度级循环迭代的方法，最佳阈值应满足类内方差最小，类间方差最大的准则，当分割阈值为T时，两个类别的像元个数分别是N₁和N₂，均值分别是m₁和m₂，方差分别为σ₁和σ₂，通过一系列的公式推导，当使得式(2)取得最大值时的T^*即为最佳的分割阈值。

虽然对数强度图并不完全符合高斯分布，但是已有的经验统计模型中，高斯分布具有最好的拟合效果。KI方法以图像的概率分布为特征，在代价函数的辅助下，当满足总体分类代价最小时即可获得具有最小错分率的结果。式(3)是高斯分布下的判别准则，Pa是先验概率，当J(T)取得最小值时的T^*即为最佳的分割阈值。

Otsu、KI等阈值分割法只是采用简单的阈值将图像分成两个类别，在实验流程中并没有加入抑制斑点噪声的机制，多视、滤波等预处理过程虽然会降低图像的噪声水平，但是也会导致图像空间分辨率的降低和水体边缘的模糊化，损失了图像的部分有用信息。应当发展一种既保持图像的细节信息，又对斑点噪声具有较好抑制作用的水体提取方法。

Fisher变换是一种图像线性变换方法，在先验知识的指导下，将多维数据通过线性变换投影到一维数据中，投影的准则是使初始图像中相同类别的像元实现最大程度地聚集，不同类别的像元实现最大程度地分离。杨国鹏等(2008)详细介绍了Fisher变换的实现过程，经过一系列的数理推导，该方法利用式(4)计算出多维影像的变换系数 ${\alpha} $ ，将多极化对数强度影像转换为一维影像 $g $ (式(5))，QPSI模式的影像需要计算 ${\alpha} $ 的4个分量，FSII模式的影像只计算 ${\alpha} $ 的前2个分量。

式中， ${{{S}}_{\rm{W}}}$ 是在KI二值化分割得到的水体分布图的指导下，计算得到的水体和非水体两个类别的类内散布矩阵之和(式(6))， ${{x}} = {[HH,HV,VH,VV]^{\rm{T}}}$ ， ${{M}}_i = {[mi,mi,mi,mi]^{\rm{T}}}$ i是水体和非水体两个类别序号。

Fisher变换将多极化的GF-3影像融合在一起，产生了水体—非水体类别具有更高可分离度的一维影像，更有利于水体信息的提取。GD-MRF以高斯分布下的最大期望算法为框架，在利用单一像元后向散射强度的基础上加入了像元间的空间上下文信息，使相同类别的像元更连贯地聚集在一起。张斌(2013)已详细描述了MRF的推导过程，本研究通过KI阈值法得到的初始水体分布图计算二阶邻域内的水体—非水体像元个数f_i，结合像元本身的类条件概率构建水体和非水体两类的能量函数U(式(7))，进而通过式(8)计算每个像元的后验概率，以最大后验概率准则确定像元的水体—非水体类别属性。式中，f是二阶邻域内与中心像元的类别属性不同的像元个数， $\lambda $ 是f的权重系数。

在Fisher变换和GD-MRF每一次的迭代计算中，不断更新初始的水体分布图、均值m、方差 $\sigma $ 、权重系数 ${\textit{λ}} $ ， ${\textit{λ}} $ 初始值设置为0.3，在实际运用中可采用牛顿—拉夫逊方法计算 $\lambda $ 的近似值。当满足相邻两次后验概率的差值小于0.001，或者迭代次数达到30次时停止迭代，输出为最终的水体信息分布图。该方法既融合了GF-3数据的多极化信息，在保持影像空间分辨率的前提下对固有的斑点噪声也具有较好的抑制作用，有效提高了算法的自动化程度。

利用Otsu、KI两种阈值分割方法对两种成像模式下的HV极化影像进行阈值分割，为研究雷达辐射失真像素对水体提取精度的影响程度，本研究首先对整幅图像进行Otsu和KI阈值分割，然后利用Mask文件剔除辐射失真区域，重新进行KI自动化阈值分割，此种方法表示为Mask-KI。对比目视解译阈值及自动化确定阈值的大小(表1)，可以看出Otsu方法的阈值与目视解译的真实阈值差别很大，难以达到理想的水体提取效果，而采用基于高斯分布的KI方法自动选择的阈值与真实阈值非常接近，可以达到较好的水体提取精度。经过Mask文件掩模以后，Mask-KI方法自动确定的阈值与真实阈值更加接近，误差在0.5 dB以内，显示了本文所发展的阈值分割方法的有效性。

表 1(Table 1)

表 1 不同阈值分割方法比较 Table 1 The comparation of threshold using different threshold segmentation method /dB 图像 目视解译阈值 Otsu KI Mask-KI FSII HV –27 2.35 –26.65 –27.01 QPSI HV –22.5 1.70 –23.69 –22.96

表 1 不同阈值分割方法比较 Table 1 The comparation of threshold using different threshold segmentation method

以Mask-KI的结果为初始的水体分布图，利用Fisher变换将Mask掩模后的双极化和全极化数据分别变换为单一图像，并输入到中GD-MRF模型中，通过不断地循环迭代实现对水体信息的进一步精细提取，得到最终水体信息分布图(图7)。显然，GD-MRF提取的水体信息分布位置和面积与图2原始GF-3影像及图3参考影像非常接近，具有较好的识别效果。

图 7 GD-MRF方法的水体分布图 Figure 7 Water distribution map based on GD-MRF method

为定性评价本文方法的有效性，体现出GF-3影像在水体信息提取方面的优势，根据肖茜等(2018)、Guo等(2017)、Sun等(2017)描述的归一化水体指数(NDWI)、改进归一化水体指数 (MNDWI)、自动提取水体指数(AWEI)3种方法，提取Landsat 8影像中的水体信息，GF-1影像波段较少，只采用了NDWI。从Mask-KI法、GD-MRF、NDWI、MNDWI、AWEI提取的水体分布图与光学影像、GF-3 HV极化影像的对照可以看出(图8)，Mask-KI法在湖泊及河流内部出现了较多的空洞，在非水体部分出现了较多的椒盐噪声，存在一定的水体—非水体混淆现象。经过GD-MRF模型的精提取以后，有效降低了噪声的影响，水体—非水体边界更加平滑，获得更加完整的河流、湖泊分布图。NDWI、MNDWI、AWEI 3种水体指数均具有一定的水体识别能力，在样点一和样点二处获得了与GD-MRF接近的水体识别效果。但是样点四中分布有较多地云层，水体指数无法区分云、云阴影和水体，造成了水体信息的错分。样点5和样点6反应了NDWI指数对城镇和水体的混淆。受制于空间分辨率的限制，基于Landsat 8和GF-1影像的水体信息提取方法无法获得较清晰的水体边界，尤其是在弯曲河流处，难以刻画更详细的细节信息。而基于GD-MRF的GF-3的水体提取结果具有较高的空间分辨率，可获得更清晰的水体信息分布图。

利用在Landsat 8和GF-1光学影像上采集的随机样点和定义的TDA评价指标，定量评价本研究所发展的水体提取方法及NDWI、MNDWI、AWEI 3种水体指数法的精度，表2列出了实验方法、参考影像、采样点数及TDA指标，在GF-3两种成像模式下，Mask-KI阈值分割法正确率为79%，经过GD-MRF精提取后正确率达到85%以上，而基于光学影像的水体指数法实验精度均在83%以下，显示出本文所发展的GD-MRF方法较强的水体信息识别能力。

图 8 水体提取细节展示 Figure 8 The detail display of water extraction

表 2(Table 2)

表 2 精度评价表 Table 2 Accuracy evaluation table 成像模式 参考影像 采样点数量 正确检测样点数量 TDA/% FSII GD-MRF GF-1 200 170 85 QPSI GD-MRF Landsat 8 100 87 87 GF-1 NDWI GF-1 200 166 83 Landsat 8 NDWI Landsat 8 100 82 82 Landsat 8 MNDWI Landsat 8 100 81 81 Landsat 8 AWEI Landsat 8 100 79 79

表 2 精度评价表 Table 2 Accuracy evaluation table

本文发展了一种针对国产GF-3影像的水体自动化提取算法，该方法利用元数据中提供的定标参数，结合DEM数据生成了地理坐标系下的遥感正射影像产品及山体阴影掩模文件，利用Mask-KI阈值分割法初步提取水体分布图，在Fisher变换和GD-MRF循环迭代过程中完成水体分布的精细化提取，在湖南省东北部地区展开实验，主要结论如下：

(1)由于水体对雷达波束的镜面散射机制，在SAR图像上可以应用阈值分割的方法进行水体信息的提取。交叉极化(HV、VH)比同极化(HH、VV)具有更强的水体识别能力。

(2)目前的SAR正射校正算法多采用图像重采样、空间插值等方法降低透视收缩和阴影等现象的影响，产生了辐射失真，对原有的SAR图像概率密度函数分布具有一定的干扰作用，去除辐射失真区域后，有助于获得更准确的水体—非水体分割阈值。

(3)Fisher变换和GD-MRF模型的结合综合利用了多极化SAR影像像元的灰度信息和像元间的空间邻域信息，对斑点噪声具有较好的抑制作用，显示出强大的水体识别能力。

(4)通过GF-3水体提取结果和参考影像的对比分析，本文所发展的基于GF-3影像的水体信息提取方法达到了较高的实验精度，实验结果优于基于光学影像的水体指数法，对洪水等自然灾害监测及GF-3影像的快速业务化应用具有重要意义。

通过对水体分布图的仔细检查发现，当河流上存在较大船只、桥梁等高反射目标时，会干扰水体提取的精度，在今后的研究中应选择更多的图像特征，实现水体信息的连贯提取。