基于GF-1影像的山区河道丰、枯水期水体提取方法改进

赵程铭，董晓华，薄会娟，章程焱，张庆玉

（1.三峡大学水利与环境学院，湖北宜昌 443002；2.三峡库区生态环境教育部工程研究中心，湖北宜昌 443002；3.水资源安全保障湖北省协同创新中心，武汉 430072）

0 引 言

从遥感影像中准确提取水体信息对于水资源的调查、规划和保护具有重要意义［1］。黄柏河东支流域位于长江北岸，是湖北省宜昌市的重要水源地，但该流域同时拥有全国储量第二的磷矿资源，磷矿开采导致位于黄柏河东支上游的玄庙观水库与天福庙水库出现了不同程度的水华，且两座水库都处于中等营养水平［2］，为遏制污染物向下游水库的进一步恶化蔓延，需要对该流域进行严格监测；而天福庙水库与下游西北口水库之间河道水体地处山区，地形蜿蜒曲折，河道细小狭长，水体实地采样难度增大。当前，遥感是一项迅速发展的技术，可以实现在人迹罕至地区对大范围的实时数据进行长时间跨度连续收集，为地方、区域以及全球尺度的水环境变化提供低成本且可靠的信息［3］。利用遥感影像对流域水体进行精确提取，建立水质反演模型，对监测流域水体富营养化情况具有明显优势［4］；因此基于遥感影像消除噪声干扰，准确提取山区河道水体信息，是合理、正确地应用遥感技术分析水环境区域变化的重要基础。

目前，利用多光谱遥感影像对大面积水体进行自动、快速提取已被许多国内外学者广泛研究。例如：Work 和Gilmer 等［5］利用单波段阈值法从遥感影像中提取水体；McFeeters［6］提出归一化差分水体指数（英文全称，NDWI）算法，利用绿波段和近红外波段来描述开阔水体；徐涵秋等［7］采用改进的NDWI（英文全称，MNDWI）算法，用短波红外波段（英文全称，SWIR）代替近红外波段增强水体反射率特征；Qiao 等［8］提出了一种自适应的水体提取方法，该方法将水体特征进行分层提取；Feyisa 等［9］提出自动水体提取指数（英文全称，AWEI），旨在通过波段相加、差分运算，提高水体像元与非水体像元之间的可分性；陈文倩等［10］提出阴影水体指数（英文全称，SWI），使山体阴影与水体之间的区分度显著提高。上述水体提取方法被广泛应用于提取遥感影像中如水库、湖泊等面积较大水域，然而对于狭长且处于复杂地形山区环境下的水体，仅利用光谱特性构建的水体提取方法并不能全面地避免环境噪声，具有一定的局限性［11］。针对山区狭长水体提取，Goumehei 等［12］利用研究区DEM 高程数据消除高海拔山体阴影和暗面影响，相较于仅使用MNDWI 水体指数精度提高4.88%；朱长明等［13］在DEM 生成水系数据辅助下，采用Landsat 遥感数据，提出一种基于自适应阈值水体指数空间优化迭代算法，水体提取正确率达到90%左右；李艳华等［14］基于国产GF-1 遥感数据应用基于规则的面向对象方法，辅以DEM 数据使山区细小水体提取总体精度达到93.5%，Kappa系数达到0.87；薛源等［15］利用GF-1 遥感影像与DEM 高程数据相结合，以总体判别精度为89.5%对山区细小河流边界进行识别。通过DEM 数据进行水文分析对山区细小水体提取边界进行修正的方法虽以较为成熟，但并未确定山体阴影分布，无法消除与水体相邻近的阴影噪声；同时，乔丹玉等人研究表明，当水体处于不同地物背景条件下，提取背景的差异决定了水体提取方法的效率［16］；地物特征的变化是造成提取背景差异的主要原因，在不同水文期典型月，丰、枯水期地物特征呈明显季节性变化，如丰水期典型月多处于夏、秋雨水充沛季节，植被覆盖率高，地物多以茂盛植被为主；枯水期典型月则多处于春、冬少雨季节，森林植被稀疏，地物多由裸露山体组成。

因此，本文首先对比传统高效的单波段阈值法、NDWI水体指数和SWI 阴影水体指数在丰、枯水期不同地物特征背景下水体提取的应用效果，确定丰、枯水期典型月水体提取的最适方法；其次对丰、枯水期水体提取结果进行综合改进处理，针对山体阴影、水系周围零散水体小斑块以及水体提取结果断线现象等问题，本文引入遥感成像时太阳高度和太阳辐射角度，对实验区中逐个像元进行亮度识别，得到实验区山体阴影图层，确定实验区山区阴影准确位置信息，并使用DEM数据计算流域坡度和坡向生成流域河网水系对水体边界进行校正，最后利用膨胀滤波和Pavlidis 异步细化填充水体断线现象处，获到丰、枯水期山区河道水体准确分布信息，以期为后续流域水资源监控与管理提供可靠支撑。

1 数据与方法

1.1 研究区概况

黄柏河流域位于宜昌市中部，地处东经111°04′～111°30′，北纬30°43′～31°29′之间典型峡谷型河流。黄柏河由东、西两支组成，其中东支发祥于宜昌市夷陵区黑良山山脉，是黄柏河的主流，其长度为130 km；西支发源于夷陵区的五郎寨，长度为78 km［17］。两支在夷陵区两河口汇合后在葛洲坝坝址上游流入长江。黄柏河流域总面积1 902 km2，年平均气温16.9 ℃，属于亚热带季风气候，降水分布具有明显的季节性，暴雨、大暴雨集中在6-8月份出现，其3个月降水量几乎占全年降水的一半［18］，根据黄柏河东支流域2016年的月降雨量及月流量情况（图2），选取7月为丰水期典型月，2月为枯水期典型月。黄柏河东支流域内水体主要以水库与河道水体为主，依次建有玄庙观水库、天福庙水库、西北口水库、尚家河水库四座梯级水库。本文选取天福庙水库与西北口水库之间山区河道水体为实验对象，研究区地理位置及实验区遥感影像如图1所示。

图1 研究区概况图Fig.1 Overview of the study area

图2 黄柏河东支流域2016年降雨-流量变化图Fig.2 Rain-flow variation in the East Branch of Huangbai River basin in 2016

1.2 数据来源

国产GF-1 卫星于2013年发射，是中国高分辨率对地观测系统中的第一颗卫星；GF-1 号配有两台2 m 分辨率全色/8 m 分辨率多光谱复合PMS（Panchromatic and Multispectral）传感器，四台16 m 分辨率多光谱WFV（Wide Field of View）传感器，重返周期4 d。本研究采用GF-1（WFV）遥感影像为数据源，数据来源于中国资源卫星应用中心（http：//www.cresda.com/CN/），WFV传感器光谱波段信息如表1 所示；选择成像时间为2016年7月14日与2016年2月18日的两景遥感影像分别作为黄柏河东支流域丰水期和枯水期进行水体提取实验。同时使用到分辨率为30 m×30 m 的数字高程数据（ASTER GDEM），数据来源于地理空间数据云（http：//www.gscloud.cn/search）。

表1 GF-1卫星WFV传感器光谱波段信息Tab.1 Spectral band information of GF-1 satellite WFV sensor

1.3 数据预处理

GF-1 遥感影像预处理主要包括正射校正、辐射定标、大气校正和影像裁剪。正射校正选择与待校正高分影像时相相近的Landsat8-OLI 全色影像为参照，采用ENVI 软件中RPC（有理多项系数）正交校正模块进行；继而将完成正射校正的遥感影像进行辐射定标，将像元亮度DN 值（Digital Number）转换为具有实际物理意义的表观反射率；最后使用ENVI 中FLAASH（Fast Line-of-Sight Atmospheric Analysis of Spectural Hypercubes）模块进行大气校正。同时裁剪出研究区域相应的ASTER GDEM 30m数据。

1.4 水体提取方法

水体提取方法通常可以分为单波段阈值法和多波段算法两种方法。单波段阈值法通常是依据多光谱图像选择出水体的特征波段，然后确定该波段的合适阈值，以区分其他典型地物；多波段组合算法则是利用不同波段之间的反射差异，通过不同的波段组合来增强水体与其他典型地物之间的差异。目前，常用的波段组合水体指数计算方法有NDWI、改进NDWI（MNDWI）、EWI、SWI 和AWEI 等，由于高分一号卫星数据的波段数量有限，因此本文选择性能较好的单波段阈值法、NDWI水体指数和SWI阴影水体指数。水体提取流程如图3所示。

图3 水体提取流程图Fig.3 Flow chart of water extraction

（1）单波段阈值法。单波段阈值法主要是依据不同地物反射率差异设定合适的划分阈值对水体进行分割。通过对实验区中水体、山体阴影、裸露山体和植被等典型地物进行光谱特征分析（图4），可以发现裸露山体、植被和山体阴影在Band1、Band2、Band3、Band4 具有相似的变化趋势，同时水体在Band2、Band3、Band4 与其他典型地物无交集，并且其表观反射率值普遍低于其他地物，呈逐渐下降趋势，在Band4（NIR）处达到反射率最低值，这一光谱特点是区分水体与其他典型地物的关键。因此，本文在两幅GF-1遥感影像数据的近红外波段（NIR）建立一个阈值进行水体提取。

图4 水体及典型地物光谱曲线图Fig.4 Spectral curves of water bodies and typical ground objects

（2）NDWI（Normalized Difference Water Index）水体指数阈值法。

式中：ρGreen和ρNIR分别为绿波段和近红外波段的反射率数值。

（3）SWI（Shadow Water Index）阴影水体指数阈值法。

式中：ρGreen为蓝波段反射率数值。

1.5 水体提取结果综合改进处理

（1）水体边界修正。通过上述3 种水体提取方法划定合适的阈值，确保水体得到完整提取，获得的最佳水体提取结果不可避免会存在少量阴影和零散水体小斑块等噪声干扰；因此本文引入遥感成像时太阳辐射角度和太阳高度（太阳高出地平线的角度或坡度），结合DEM 高程数据（ASTER GDEM 30 m）计算得到的实验区坡度和坡向，对实验区中逐个像元进行亮度识别，得到包含山坡阴影和斜坡暗影等信息的山体阴影图层（图5），确定实验区山体阴影准确分布信息，使山体阴影得到全面可视化，并将其与上述提取的水体相叠加，去除阴影噪声干扰；并依据DEM高程数据获得的实验区河网水系，确定河道水体整体分布，通过高程落差限定，去除实验区河道干流水体周边零散水体小斑块，对经过最佳阈值划分的河道干流水体提取边界进行修正。

图5 实验区山体阴影与河网水系叠加效果图Fig.5 Overlay effect of mountain shadow and river network in the experimental area

（2）水体断线处填充。由于实验区水体形态分布呈狭小细长状，同时受丰、枯不同水文期水量影响，在提取过程中部分水体存在明显的断线现象。膨胀滤波算法被广泛应用于填充比自身结构元素小的像元结构中［19］，通过统计出断线处像元数量级，设定相同数量级的结构元素，即可通过公式（3）将断线处加以膨胀连接。Pavlidis 异步细化算法可以在保持原影像的像元结构基础上，识别并标记像元骨架，将影像上多像元宽度的直线、曲线沿中心轴线细化成适当宽度［20］。使经过填充后的膨胀水体细化为实际宽度的水体。

式中：f(s，t)为图像输入函数；b(x，y)为结构元素；Df、Db分别是f和b的定义域。

1.6 精度验证

总体精度［21］（Overall accuracy，OA）表示所涉及到的所有像元分类的正确性；Kappa系数［22］可以从整体上对两幅图中像元点的一致性进行评价，针对实验区水体信息提取精度进行全面定量评估，本文通过从丰枯水期水体提取结果影像中随机抽取100 个像元点与参照样本之间建立误差矩阵，计算得到总体精度（OA）和Kappa系数；参照样本数据来源于同时期Google Earth 提供的研究区高空间分辨率影像。Kappa系数数值与质量分级关系见表2［23］。

表2 Kappa系数数值与质量分级关系表Tab.2 Relationship between Kappa coefficient value and quality grade

式中：Oreal为真实水体信息提取区域像元数量；Ototal为标准参照图像中的真实水体信息像元数量；K为Kappa系数；N为总像元数量；xi+、x+j分别为第i行和第j列总像元数量；xij为第i行第j列单元上的像元数量。

2 结果与分析

单波段阈值法、NDWI 水体指数和SWI 阴影水体指数3 种水体提取方法在实验区丰、枯水期的水体提取结果如图6所示。对比3 种水体提取方法，在丰水期典型月，水体提取结果［图6（a）、6（b）、6（c）］受山体阴影影响较小，归因于在丰水期实验区处于夏季，水量充沛，地物生长茂盛，仅有少量植被暗影会被遥感卫星所捕捉，通过单波段阈值法即可获得良好的提取效果。在枯水期典型月，NDWI 水体指数［图6（e）］和SWI 阴影水体指数［图6（f）］相比单波段阈值法能较好地提取河道水体，提取结果受少量山体阴影和零散水体小斑块影响，而单波段阈值法［图6（d）］则会将大量的山体阴影误认为水体被提取；主要因为在枯水期实验区正直冬季，森林植被较稀疏，地形连绵起伏，在遥感影像中从低山河谷地带到高山区极易存在大量的山体阴影。

图6 丰（上）、枯（下）水期3种方法水体提取结果图Fig.6 Water extraction results of the three methods in the full（top）and dry（bottom）water periods

通过计算OA和Kappa系数（表3）对3 种水体提取方法在丰、枯水期提取结果进行精度评价，表明在丰水期对水体提取的最适方法为单波段阈值法，Kappa系数达到“极好”水平（Kappa系数=0.80）；枯水期的最适方法为SWI 阴影水体指数，Kappa系数同样达到“极好”水平（Kappa系数=0.86）。

表3 丰、枯水期水体提取及综合改进处理结果精度评价结果Tab.3 Precision evaluation results of water extraction and comprehensive improved treatment in wet and dry seasons

为进一步去除噪声干扰，提高山区河道水体提取精度，对丰、枯水期最佳水体提取结果进行综合改进处理，依据山体阴影图层和河网水系对丰枯水期水体提取边界做修正处理，结果分别如图7（a）、7（f），阴影噪声明显改善。然而噪声去除后的山区河道水体，受混合像元效应影响，在河道狭窄、转弯处部分水体断线现象较为严重（红色框线），导致局部水体信息丢失，其中图7（b）、7（g）分别为丰、枯水期水体断线处放大效果图；针对水体断线处本文使用3×3 的变换单位矩阵进行膨胀滤波处理，经过3 次迭代［图7（c）、7（h）］，水体断线处得到连接，在膨胀滤波处理基础上进行Pavlidis 异步细化处理，得到丰枯水期水体提取最终结果［图7（d）、7（i）］，从断线填充效果局部放大图7（e）、7（j）中可以看到断线处连接良好。

图7 丰（上）、枯（下）水期水体提取综合改进处理效果图Fig.7 Effect diagram of comprehensive improved treatment of water extraction during abundant（top）and dry（bottom）periods

经综合改进处理后，丰水期水体提取结果总体精度从89.54%提高到99.52%，Kappa系数从0.80 提高到0.98；枯水期总体精度从91.65%提高到99.27%，Kappa系数从0.86 提高到0.97；水体提取精度达到“极好”最优精度等级。为进一步验证水体提取效果，将综合改进后的面状水体与真彩色原始遥感影像进行叠加，图8两幅影像清楚地显示了丰、枯水期水体提取结果与实际水体边界基本吻合。通过统计丰、枯水期水体提取结果像元个数，得知水面面积呈下降趋势，在丰水期水面面积达到350 km2，当枯水期到来时水面面积减小到194 km2。

图8 丰（左）、枯（右）水期水体提取结果叠加图Fig.8 Superposition diagram of water extraction results in abundant（left）and dry（right）periods

3 结 论

本文基于GF-1遥感数据提取了2016年丰枯水期典型月黄柏河东支流域上游天福庙水库与下游西北口水库之间山区河道水体，针对不同水文期确定了最佳水体提取方法，并对丰、枯水期最佳提取结果进行综合改进处理，得到以下结论。

（1）通过对比单波段阈值法、NDWI 水体指数、SWI 阴影水体指数三种水体提取方法在不同水文期的提取效果，发现在丰水期，水体提取干扰噪声主要以植被暗影为主，单波段阈值法为最佳提取方法，水体提取结果总体精度达到89.54%，Kappa系数为0.80；在枯水期，实验区遥感影像中存在大量山体阴影，SWI 阴影水体指数为最佳提取方法，水体提取结果总体精度达到91.65%，Kappa系数为0.86。

（2）基于山体阴影图层、河网水系、膨胀滤波和Pavlidis 异步细化算法对丰、枯水期水体最佳提取结果做综合改进处理后，丰水期水体提取总体精度提高到99.52%，Kappa系数提高到0.98；枯水期水体提取总体精度提高到99.27%，Kappa系数提高到0.97；水体提取精度均达到“极好”最优精度等级。