翟涌光, 张 鑫, 冀鸿兰, 牟献友, 张宝森
(1.内蒙古农业大学水利与土木建筑工程学院,内蒙古 呼和浩特 010018;2.黄河水利委员会黄河水利科学研究院,河南 郑州 450003)
冰凌由河流中的河水冻结、堆积而产生,在河道中堆积易引发冰塞、冰坝等现象,轻则河冰堆积上岸,重则造成漫堤、溃堤决口[1]。黄河内蒙古段河道形态多样,河势复杂,易产生冰凌灾害[2]。因此,准确提取河冰、水体信息对研究河冰生消过程、河冰空间变化以及冰凌洪水的研判至关重要。传统野外观测局限于“点对点”式的局部观测,无法获得完整的河冰信息,冰情观测受到限制。卫星遥感则具有观测范围大、时效性强、数据内容客观等优点,近些年在河冰研究领域已被广泛应用[3]。
目前,利用遥感技术提取河冰、水体信息主要以目视解译、单波段反射率阈值、归一化积雪指数(NDSI)、归一化差异水体指数(NDWI)以及改进的归一化积雪指数(MNDSI)与单波段阈值相结合的方法为主[4]。杨中华等[5]基于“四星三源”模式,在国内最早将遥感技术应用于黄河冰情监测当中,限于当时卫星影像空间分辨率,作者仅进行大尺度宏观观测,并未进行局部河段微观分析;赵水霞等[6]基于Landsat 8 彩色合成影像对黄河什四份子弯道河冰过程及分布进行目视解译,但作者并未进行提取精度验证;Chu 等[7]基于MODIS 数据近红外波段监测了加拿大奴隶河河冰过程,作者并未确定河冰、水体在近红外波段的反射率值,而是假设当近红外波段的反射率最大时河道内均为河冰,而反射率最小时河道内均为水体;Chaouch 等[8]利用MODIS 可见光及近红外波段的阈值构建决策树,提取河冰信息并获取河冰范围;Dozier[9]于1989 年提出的用于监测积雪变化的雪信息提取NDSI指数,该指数对积雪与土壤的区分度较高,目前也被应用于大尺度冰川、湖冰等领域[10-13];Li 等[14]利用NDSI 对1999—2018年青藏高原八宝河河冰分布进行监测,而八宝河流域冬季河水完全封冻并不存在清沟(狭长型不结冰的裸露水面),提取任务主要以区分河冰与土壤像元为主,NDSI 分类精度令人满意;牟献友等[15]根据黄河水体的光谱特征,以蓝色及绿色波段构建归一化差异未封冻水体指数(NDUWI),在利用NDSI 排除土壤、树木等其他地物像元的基础上,进而利用NDUWI进行河冰、水体提取,但以水体作为切入点的提取方法的精度有待进一步验证;勾鹏等[16]发现近红外波段与短波红外波段构建的归一化指数直方图中两峰值间距大于NDSI 直方图中的两峰间距,有利于区分河冰、水体,由此建立改进的归一化差异积雪指数(MNDSI)对纳木错湖冰进行识别,取得很好的效果。综合来看,对于黄河河冰、水体的监测,当前遥感技术和卫星影像的时空分辨率已基本满足冰期黄河冰凌监测的能力,但单一使用遥感指数模型存在难以区分河冰、水体与土壤的问题,河道边界难以提取;同时,各指数模型对于黄河河冰、水体监测的精度还有待验证。
黄河内蒙古段由于其自身河道水力条件的特殊性,凌汛期易发生卡冰结坝现象,增加“武开河”的风险[17]。鉴于此,本文选取黄河内蒙古段5 个试验河段为研究对象,基于Landsat 8 OLI 遥感影像数据,在利用Google Earth 历史高清影像结合NDSI 指数模型排除河道外无关地物的基础上,使用近红外波段反射率值、NDWI、NDSI、MNDSI 以及NDUWI,利用各指数阈值对河冰、水体进行分类,进而计算监测精度,对阈值的稳定性及方法的适用性进行评价,以期为遥感技术监测河流冰情及冰凌洪水的研判提供参考。
研究区始于黄河内蒙古段海勃湾水库,终于万家寨水库库尾区,全长约为650 km,整体呈“几”字形分布[18](图1)。分别在河段上、中、下游选取试验河段,河段类型包括过渡型、游荡型及弯曲型;过渡型河段浅滩多且河道摆动幅度大;游荡型河段淤积严重且经常发生河冰漫滩、偎堤等现象[19]。实验区内地物类型包括主槽冰、河滩冰、清沟及土壤,河冰主要以平滑、冲击岸冰为主要类型。河段每年12月初开始流凌,至翌年3月末河段冰凌全部消融,受来水条件、河道特性等因素影响,水体在河床内及河滩上结冰;在流速较大的急流、浅滩处或排水口处常出现狭长型不结冰的裸露水面(清沟)。
图1 研究区概况Fig.1 Overview of the study area
1.2.1 Landsat 8 遥感影像河冰信息提取基于Landsat 8 OLI影像(后文简称L8影像)数据,数据获取自美国地质调查局官网(United States Geological Survey,简称USGS,https://earthexplorer.usgs.gov/),所用影像均为L1T 级产品,该级别影像是使用地面控制点和数字高程模型数据进行精确校正的数据产品,误差小于0.4 个像元。试验所用使用的影像为无云或河道及河道附近无云的影像数据。影像信息见表1。
1.2.2 验证数据试验使用Google Earth 历史高清影像作为参考影像,选取与L8影像拍摄日期同日或相近日期的历史高清影像,由于拍摄日期处在河冰稳定封冻的1—2月,L8影像与参考影像成像时间差内的河冰覆盖基本不变。基于参考影像目视解译出的各类主要地物矢量数据将用于分类精度评价。因参考影像的空间分辨率在亚米级,故目视解译的误差将忽略不计。
1.2.3 数据预处理L8影像的预处理过程为辐射定标及大气校正。预处理过程利用集成在ENVI 5.5版本中的辐射定标及FLAASH(Fast line-of-sight atmospheric analysis of spectral hypercubes)大气校正工具获得地表反射率数据。其中,FLAASH 中大气模型参数根据影像拍摄日上午的大气柱水汽量确定,该数据获取自MODIS标准大气产品中的MOD05数据集(下载地址:http://ladsweb.nascom.nasa.gov)。大气校正参数见表1。
表1 Landsat 8 OLI遥感影像信息及大气校正参数Tab.1 Landsat 8 OLI remote sensing image information and atmospheric correction parameters
所有L8影像均经过堤防矢量数据进行裁剪,只保留两岸堤防间的部分。首先基于Google Earth 历史高清影像目视解译绘制的河道矢量边界进行裁剪,排除大部分河道外无关地物,其次利用NDSI 对河岸边缘土壤及其他无关地物进行二次排除。因实验区位于两岸堤防内的河道范围,相关法律规定堤防内部不得设立建筑、道路等设施且河滩地不得占用,故经堤防矢量裁剪后的影像数据可有效排除附近民居、道路、工厂等其他无关地物;又因河道内树木极少且处于冬季,故假定河道内仅有河冰、水体、土壤3种地物。而NDSI对土壤和冰雪的分类精度已有前人证实在此不再赘述[9]。为尽可能保留全部河冰、水体像元,故以0 作为土壤与河冰、水体的分类阈值,影像预处理及分类流程见图2,流程化处理后待提取像元仅为河冰与水体,几乎无其他无关地物。
图2 影像预处理及分类实验流程图Fig.2 Flowchart of image preprocessing and classification experiment
1.3.1 5 种遥感指数模型选取5 种遥感指数模型对黄河内蒙古段5个子段的河冰、水体进行分类,分别为:
(1)近红外波段反射率阈值法[7]。水体在近红外波段反射率值低,而冰雪的反射率较高,故适用于冰雪与水体的分类。分类算法如下:
式中:ρNIR为近红外波段反射率值;pixel 为遥感影像像元点;threshold 为阈值,下同。为节省储存空间,USGS 发布L8 影像数据的像元反射率值均放大10000倍,计算时除以10000将取值范围限制在0到1之间以便于与其他指数值进行比较。
(2)归一化差异水体指数(NDWI)阈值法[20]。NDWI利用绿色波段与近红外波段构建归一化差异指数,是由McFeeters[20]提出的目前最为常用的陆表水体指数。因河道内并不涉及建筑用地,NDWI 的缺点不会影响其对河冰、水体提取。分类算法如下:
式中:ρSWIR为短波红外波段反射率值,下同。
(4)改进的归一化差异积雪指数(MNDSI)[23]。MNDSI 利用近红外及短波红外波段构建归一化差异指数,而其计算方法与Ouma 等[23]提出的系列水体指数中的NDWI3计算方法一致,该方法对水体边缘提取效果较好。分类算法如下:
式中:ρBlue为蓝色波段反射率值。
1.3.2 精度评价为比较各算法分类精度,利用总体分类精度(p0,指对每一个随机样本,所分类的结果与检验数据类型相一致的概率,单位为%)和Kappa 系数(κ,基于混淆矩阵衡量河冰、水体分类精度的指标,介于[-1,1])评价各分类方法的提取精度(pe为中间变量,表示偶然一致性,单位为%),其计算过程见式(6)~(8)。根据各方法的提取结果列出河冰、水体混淆矩阵(表2),并分别计算各指数分类的p0及κ。Landis 等[24]提出κ>0.81 时分类几乎完全吻合,但为达到精确提取的目的,选取κ>0.9作为分类达到几乎完全吻合的标准。
表2 河冰、水体混淆矩阵Tab.2 Confusion matrix of river ice and water body
分别采用5种遥感指数模型对黄河内蒙古段河冰、水体信息进行提取,限于篇幅,仅展示2016年和2018 年2 期影像提取结果(图3、图4)。总体来看,基于NDSI结合Google Earth历史高清影像提取到的河道边界排除了土壤及河岸滩地上冰雪等地物,降低了无关地物对于后续总体分类精度及Kappa系数计算精度的影响,从分类效果来看,5种方法分类效果由高到低依次为:归一化未封冻水体指数、归一化差异水体指数、近红外波段反射率值、改进的归一化差异积雪指数、归一化差异积雪指数。在2018年影像中,NDSI 较其余4 种方法在河冰、水体信息提取上欠佳(图3),误分现象严重;NDWI与NDUWI在遥感影像中分类效果最佳,其中NDUWI 总体分类精度和Kappa系数均达到95.00%和0.90以上。
图3 2016年2月1日5种遥感指数模型河冰、水体分类结果Fig.3 Classification results of ice-water in five remote sensing models on February 1,2016
图4 2018年2月22日5种遥感指数模型河冰、水体分类结果Fig.4 Classification results of ice-water in five remote sensing models on February 22,2018
从5种遥感指数模型所对应的最大(最小)阈值来看(表3),河冰、水体分类结果与阈值大小并无直接关系,NDSI 阈值范围最大,但其分类结果在5 个试验区并不稳定,NDWI 与NDUWI 阈值范围适中,但分类效果在5个试验区内分类结果较好。进一步分析发现,NDUWI 在5 期遥感影像中的河冰、水体最优区分阈值大体分布于阈值中值附近,如在2016年遥感影像中NDUWI的阈值为(0~0.46),最优阈值为0.21。
表3 5种遥感指数模型对应最大(最小)阈值范围Tab.3 Maximum(minimum)threshold range of the five remote sensing models
2.2.1 敏感性分析图5 给出了试验区5 期遥感影像河冰、水体分类精度对于阈值变化的响应。从各期影像中各方法的分类精度来看(表4),在2020 年的2期影像中近红外波段反射率阈值法在相较其他方法展现出最高的提取精度,p0及κ分别为99.18%、0.99及99.98%、0.99;2014、2016年及2020年影像中NDUWI 阈值法相较其他方法展现出最高的分类精度,p0及κ分别为99.38%、98.90%、99.80%及0.98、0.97、0.99。从各方法在各期影像中的分类表现来看,近红外波段反射率阈值法在2018年1期和2020年2 期共3 期影像中p0及κ均达到并超过了90.00%及0.90,2014年和2016年影像中p0未达到90.00%,κ未超过0.7;与近红外波段反射率阈值法相似,NDWI在前2期影像中分类效果不理想,在后3期影像中p0及κ均达到并超过90.00%及0.90;NDSI 在2014 及2016 年2 期影像中p0及κ均达到并超过90.00%及0.90,而在后3 期影像中κ出现为0 的情况,说明误分、漏分现象严重;MNDSI 在5 期影像中p0及κ总体保持在85.00%及0.80左右;NDUWI阈值法在5期遥感影像中的p0及κ均达到并超过了90.00%及0.90。
表4 5种遥感指数模型对应最高分类精度Tab.4 Highest classification accuracy of the remote sensing index models
进一步分析发现(图5),当采用适当阈值时,对于上述5 种遥感指数模型,所能达到的最高精度从低到高排序依次为:NDSI、MNDSI、近红外波段反射率阈值法、NDWI、NDUWI。此外,从分类精度对阈值变化的响应机制来看,近红外波段反射率阈值法变化幅度最大,即在分类时只要阈值有微小的变动,就会造成提取结果的显著变化,总体精度差值为50.10%;相比较而言,NDWI阈值变化对分类精度敏感程度较低,选取恰当阈值,便可得到较好的精度;对于其他3种指数模型,阈值对分类精度影响的敏感程度由低到高分别为:NDSI、MNDSI、NDUWI。总体来看,NDUWI 对阈值变化的敏感程度较小,总体p0及κ较MNDSI 高,在适当阈值内p0及κ更接近100%和1,分类精度高。
图5 试验区p0和κ对阈值的响应Fig.5 Response of p0 and κ in the test area on the threshold
综合试验区5 期提取结果发现,L8 影像不同指数模型在同一研究区内阈值对河冰、水体提取精度的影响程度不同;相同指数模型在不同研究区内所能达到的提取精度不同,阈值对河冰、水体提取精度影响程度也不同;例如,在2016年遥感影像中,当选用一定阈值时,NDUWI精度最高,MNDSI次之;而在2018年影像中,NDUWI在一定阈值内精度最高,近红外波段反射率阈值法次之。
2.2.2 不确定性分析基于遥感影像数据对黄河内蒙古段冰水信息提取中,由于主观和客观因素的存在,往往影响到最后的提取精度。例如,NDSI 在2014年和2016年2期影像中分类精度较高,而在后3期影像的分类上差强人意。本研究的不确定性来源于以下几个方面:遥感影像的选取、人工目视解译经验及研究区环境变化。冬季,由于河道上空雾、云及霾的存在,光学数据在冰情监测中可用数据量有限,每个试验区仅能甄选1~2 景可用于试验的影像[25]。本研究所采用的L8 影像均为无云或河道及河道附近无云的影像数据,尽量保证试验数据的质量及可靠性。受客观因素影像,目前尚未引入其他卫星影像作为补充数据源。人工目视解译的误差主要来自混合像元,利用阈值模型计算的栅格数据一般包含多种地物信息,除河道内部少量平滑冰外[26](即静态岸冰,指平均流速低于河道临界流速的水力条件下,在河道内同时沿纵向和横向发展的表面较为光滑的河冰),少有纯净像元,特别是河道边界河冰、水体、积雪以及土壤之间界限难于准确拟定,在人为判别时易影响分类精度[27]。在变化环境方面,大气中水汽、气溶胶及臭氧、河道清沟内水体的含沙量、河岸边界土壤湿度等造成“同谱异质、同质异普”等现象,对边界提取、冰水情监测造成影响[28-29]。为了最大程度上降低河道边界提取中的误判,本研究设计了NDSI 结合Google Earth 历史高清影像对边界提取的流程化提取操作,使得土壤等无关地物最大程度得以排除。
黄河内蒙古段全长约为650 km,东西跨度较大,整体呈“几”字形分布,野外冰情受到限制,基于遥感影像在冰情监测方面效果较好[18]。同时,黄河冰情信息监测受限于数据的时相分辨率以及天气条件,因此,常在研究中出现某一稳封期数据缺失的情况[6]。与传统的大尺度海冰(湖冰)冰情监测不同,由于黄河自身的泥沙特性以及河道淤积蔓延情况,本文研究基于L8影像,综合Google Earth历史高清影像和NDSI首先排除了河道外土壤、建筑及植被等无关地物,仅保留冰水信息以保证被研究像元的真实性,解决了河道边界河冰、水体信息难以识别的问题[30-33]。
本研究发现,遥感影像选取、河道上空气象条件以及河流水质条件是影响冰情监测提取精度的重要因素。在数据选取方面,由于冬季河道上空雾、云及霾的存在,不同指数分类阈值难以拟定,存在较大阈值波动情况,故应选用稳封期(即河冰不再沿河道横向发展)无云、雪且气象条件良好时期的卫星影像[34]。在冰水信息提取实验设计方面,NDSI 虽不能精确区分河道内部冰情、水情信息,但该指数可结合历史高清影像作为分类化流程实验的第一步,即提取河道边界,排除研究对象外的无关地物,使真实河冰、水体像元尽最大限度得以保留[9]。
黄河内蒙古段河冰信息提取与前人大尺度海冰(湖冰)提取的研究相比,既具统一性,也具差异性。例如,海勃湾库区是典型的库区型河段,水流流速相对较小,河冰信息提取方法同海冰(湖冰)提取大体相同;巴彦高勒河段则为典型的游荡型河段,河段淤积严重,河冰漫滩现象严重,河冰信息提取受不确定因素增加。5种遥感指数模型在黄河内蒙古段冰情信息提取精度对阈值的响应不尽相同,近红外波段反射率阈值法提取精度对阈值响应较为敏感,即阈值的少量波动就会对提取精度产生较大影响;NDSI、MNDSI 以及NDWI 在流程化实验中提取精度对阈值响应适中,但精度略小于近红外波段反射率阈值法;NDUWI 在多时相、多河段研究中表现出较高提取精度,同时发现其河冰、水体最优区分阈值大体分布于阈值中值附近。同时发现,凹岸侧模型指数值略大于凸岸侧,这可能与弯道处离心力的存在导致凹岸侧堆冰增加有关,故河道凹岸侧应加强堤岸防护[35-36]。
黄河冰情遥感监测不仅受数据质量、气象条件以及光照等的影响,并且也受水位,流量以及上游水库防凌调度的影响[37-39]。例如:气温越高,河冰分布范围越小,稳封期清沟所占比例越大,已发生凌汛灾害;水位越高,河道与大气表面接触面积就越大,河冰冻结所需要的热量条件也随之增大,河冰分布范围增大,河冰漫堤风险随之增加等[40]。因此,在未来研究中,如何在模型中结合气象条件、水利条件等数据以及自动化提取冰水像元,同时结合水文站野外观测资料做对比分析,是一个值得考虑的方向,这可能会进一步提高凌汛期黄河河冰、水体情监测的时效性和准确性。
黄河内蒙古段河势复杂、地形地貌多样,基于遥感影像模型监测河冰、水体情分布不易获取且各指数模型的适用性尚未有人做出评估。针对上述问题,本文运用L8 影像数据,通过设计实验方案提取河道边界,利用近红外波段反射率阈值、NDWI、NDSI、MNDSI 及NDUWI 5 种遥感模型对研究区河冰、水体分布进行监测,进而分析各指数的阈值稳定性及不确定性,得出以下结论:
(1)基于NDSI 与Google Earth 历史高清影像结合的二次排除无关地物的实验方案在河冰、水体情监测时能较好的去除周围无关地物。
(2)NDUWI 在黄河内蒙古段影像中的分类精度均达到90.00%及0.90 以上,总体分类精度较高,适用性较强,且NDUWI 在5 期遥感影像中的河冰、水体最优区分阈值大体分布于阈值中值附近;但NDUWI无法将土壤与河冰进行有效区分,使用前应先提取河道矢量边界,排除无关地物。