文广超, 李 兴, 吴冰洁, 王晓鹤, 谢洪波
(河南理工大学资源环境学院,河南焦作 454000)
湖泊参与流域水循环,对维持区域水量平衡[1]、满足生产生活用水[2]、河流流量调节[3]发挥着重要作用。柴达木盆地地处青藏高原东北部,盆地内湖泊众多,且多位于盆地中心的低洼地带[4],为地表与地下水的汇集区和河流水系的尾闾,是气候变化与人类活动的敏感指示器[5],湖泊及其引起的生态环境的变化在一定程度上可以反映区域乃至全球的气候事件[6],在干旱的内陆盆地区,湖泊是维持区域生态系统的重要支撑,同时,盆地内的盐湖作为资源也是当地的经济支柱[7]。近年来,随着全球气候变化和人类活动干扰的加剧,盆地内湖泊出现了面积萎缩乃至干涸的现象[8-9]。湖泊面积的变化不仅影响当地的社会经济和生态环境的可持续发展,而且对区域乃至青藏高原气候变化均具有重要影响[10]。因此,精准监测柴达木盆地湖泊水体动态,揭示自然与人类活动对湖泊的影响规律,对合理开发利用和保护湖泊水域,维系区域生态平衡具有重要的意义[11]。
卫星遥感技术以其探测范围广、信息容量大、实时性与动态性强、免费共享等优势,在湖泊水体信息提取中得到了国内外学者的青睐[12-16],湖泊水体的提取最早主要依靠目视解译,目视解译提取水体的效果好坏受人为主观因素影响,而且耗费的时间成本,人力成本高。用计算机自动提取水体的早期研究主要是利用水体在单个波段上表现出的与非水体截然不同的特性,如1985 年,Jupp 等[17]通过对TM7(TM 影像的第七波段)设定阈值识别水体信息;1992年陆家驹等[18]通过对TM5设定阈值识别水体信息,受限于分辨率,对于1000~4000 m2之间的水体识别起来还有困难。随后的研究利用了水体在多个波段上的组合特征,产生了多种水体指数方法,以1996 年McFeeters[19]提出的归一化差异水体指数NDWI(Normalized Difference Water Index)为代表,通过抑制植被和土壤信息,提取了水体信息,但其对建筑物、阴影、云影、细小河流的区分度不高;杨存建等[20]利用水体在TM2 波段与TM3 波段的亮度值之和大于TM4 波段与TM5 波段的亮度值之和这一特征提取水体;徐涵秋[21]提出了改进的归一化差异水体指数法MNDWI(Modified NDWI),较好地提取了城镇水体信息,该方法仅适用于拥有SWIR的传感器,易受冰雪的影响;闫霈等[22]提出增强型水体指数EWI(Enhanced Water Index),在提取水体时区分了半干涸河道与背景噪声,但对细小水体分离效果差,出现湖泊漏提现象;丁凤[23]将TM7 波段应用于构建新型水体指数NWI(New Water Index)上,能够清晰提取水体信息,但无法区分山体;2014年王晴晴等[24]提出简单比值型水体指数SRWI(Simple Ratio of Water Index)用来增强水体与植被、建筑物、裸地和阴影的差异;丁占峰等[25]基于对鄱阳湖水体信息提取的研究提出ONDWI(Override NDWI)用来去除水体中含有的沼泽湿地,同时抑制植被和土壤信息,但对冰雪和阴影的区分度不高。此外,在沙漠水体信息提取[26-27]、城市水体信息提取[28-29]、水域边界或河流信息提取[30-31]、背景噪音去除[32-33]及提高土壤、植被与水体信息区分度[34-35]等方面,学者们也做了较为深入的研究,建立了一系列有效实用的水体指数。针对柴达木盆地湖泊水体动态变化,前人主要利用目视解译、单波段阈值法、NDWI、MNDWI 等方法开展研究,取得了一系列成果[36-38]。以上水体提取方法中,单波段阈值法与比值法依赖动态的阈值对湖泊水体信息进行提取,且对阴影、冰雪噪声的去除效果不理想;基于各种水体指数法提取水体信息的阈值基本稳定为正值,因此将经水体模型运算后结果大于0的区域判别为水体来提取信息,然而试验表明NDWI 在提取西部湖泊水体信息时会将沼泽湿地与冰雪一起提取出来,MNDWI、EWI、NWI、ONDWI 等方法则存在湖泊漏提和误提冰雪的情况。
随着信息处理、定量遥感等相关技术的发展,遥感解译正由半自动化向全自动智能化解译方向发展[39],对遥感解译的速度、分类精度与适用性提出了更高的要求。为了解决上述问题,本文以青海省德令哈市境内的可鲁克湖流域为试验区,以可鲁克湖、托素湖、尕海为试验对象建立水体识别模型,为了验证模型的识别精度,选择了柴达木盆地中的大柴旦湖、小柴旦湖、尕斯库勒湖和察尔汗盐湖4个典型湖泊作为模型的验证区域。同时,为了验证模型在柴达木盆地以外区域的适用性,将该模型应用于邻近柴达木盆地的青海湖、尕海湖,以及太湖、天池和滇池的水体提取中。基于Landsat 系列卫星影像,进行湖泊水体信息自动提取方法研究,以精度高、速度快、适用性强为目标,旨在提出一种适合柴达木盆地区域的自动提取湖泊水体信息的方法,进而为湖泊水体的动态监测提供可靠的技术支撑。
柴达木盆地位于青藏高原东北缘,地处青海省西北部(34°40′~39°20′N、90°00′~99°20′E,图1)。中部的平原区海拔2676~3000 m,地势低洼,四周高山环绕,西南暖湿气流难以进入,降水稀少、蒸发量大、日照充足,属于内陆极端干旱高寒气候。盆地内共有大小河流79 条,这些河流以盆地为中心,呈聚合状分布,河水通过地表及地下径流最终注入各自的汇水中心,形成尾闾湖[40]。
图1 研究区位置及采样点分布Fig.1 Location of the study area and the distribution of sampling points
1.2.1 遥感数据 选用Landsat 系列卫星影像数据作为数据源,数据从美国地质勘探局网站(https://earthexplorer.usgs.gov/)下载,空间分辨率30 m,影像级别为Level 1TP(地形校正)级。除了柴达木盆地内的影像外,还收集了柴达木盆地内外其他区域不同时期、不同传感器的共计90 个Landsat 图像场景数据,旨在验证本文提出的湖泊水体提取方法在不同时间、空间以及传感器的普遍适用性。本文选取自2000—2020年可鲁克湖、小柴旦湖、尕斯库勒湖、察尔汗盐湖、青海湖等区域的部分影像进行论述,数据基本信息见表1。
表1 所用部分遥感影像基本信息Tab.1 Basic information of some remote sensing images used
1.2.2 遥感数据预处理 由于影像上的DN(digital number)是由辐射亮度值经过线性变换拉伸到(0~255)范围内后得到的灰度值,不能准确表征地物在每个波段真实的反射率,因此在ENVI 软件中利用产品提供的影像定标文件将DN值转换成TOA(Topof-atmosphere)反射率。TOA反射率相对于DN值有两大优点:一是消除了传感器在不同时间获取数据时太阳天顶角带来的余弦效应;二是补偿了外大气太阳辐照度的不同值[41]。
1.2.3 采样点布设 为了获取不同地物的波谱曲线,在对影像进行图像增强的基础上,结合谷歌地球上同区域高分辨率影像,采用目视判读法共收集了55615 个样本点,其中包括26676 个水体样点,28939个非水体样点。
此外,在青海湖区域2018 年7 月的影像发布后,于当月在距青海湖东北3 km 的尕海湖周边,对尕海湖水体和非水体之间的边界线采集了8406 个实测定位点数据,这些采样点组成的虚线代表实际的湖岸线,用于湖泊水体信息提取结果的视觉评估以及精度评价参考。
1.2.4 技术路线 模型建立分三步完成:第一步,基于下载的Landsat 系列影像,借助ENVI 软件,利用Landsat 产品提供的元数据文件,将DN 值转换为TOA反射率值;第二步,通过实地调研,结合目视解译,构建不同地物的波谱曲线,分析湖泊水体与干扰地物要素之间的差异,初步提出湖泊水体信息提取模型;第三步,对模型提取精度进行验证,检验精度超过95%时,确定建立湖泊水体信息提取模型,检验精度较低时,重复第二和第三步操作,湖泊水体信息提取模型建立的思路见图2。
图2 湖泊水体信息提取模型建立的技术路线Fig.2 Technical route for establishing lake water body information extraction model
1.2.5 精度评价方法 为了验证湖泊水体信息提取的精度,引入视觉评估法和混淆矩阵[42]2种方法。
混淆矩阵[43-44]中的总体分类精度OA(公式1),Kappa 系数(公式2),以及用户精度User accuracy(公式3)等参数可以定量反映分类效果。
式中:n是总采样测试点;nkk是第k类地物样本点的总数;nk+是从分类数据中导出第k类样点的总数;n+k是从数据集导出的第k类地物样本点的总数;q是样本类别的总数。
总体分类精度(OA)、Kappa 系数和用户精度的值越大,说明分类效果越好。
不同类型的地物对一定波长的电磁波表现出不同的反射特性,地物反射电磁波特性在遥感影像上的定量表现就是TOA反射率数值不同,通常用波谱曲线来刻画地物在不同波段的TOA 反射率表现。根据不同地物TOA反射率的差异构建模型,即通过增强目标地物抑制非目标地物的方式,将目标地物与其他地物区分开。
在柴达木盆地内,湖泊水体的矿化度存在较大的差异,有淡水湖、半咸水湖、咸水湖和盐湖4类,平均高程约2788 m[40]。在可鲁克湖流域,有可鲁克湖(淡水湖)、托素湖(咸水湖)与尕海(盐湖)3个湖泊,矿化度分别为:0.80~0.88 g·L-1、29.62~30.20 g·L-1、103.6~156 g·L-1,流域平均高程2980 m[45]。无论从地形地貌还是矿化度上看,可鲁克湖流域地物要素的光谱特征在柴达木盆地均具有典型性。为此,本文选择可鲁克湖流域作为试验区开展湖泊水体信息提取方法研究。
基于上述思路,在可鲁克湖周边区域,在实地调研的基础上,通过目视解译选择对湖泊水体信息提取干扰较大的沼泽湿地、冰雪、阴影、裸地、建筑物等5类典型要素,通过样本点(图1)获取其在不同波段上的TOA 反射率,分别求取其平均TOA 反射率,以波段为横坐标,TOA 反射率为纵坐标,绘制5类典型地物和湖泊水体的波谱曲线(图3)。
图3 研究区6种地类的波谱曲线Fig.3 Spectral signatures of six major land use classes in the study area
从图3 可以看出,湖泊水体的TOA 反射率与其他5类地物有明显的不同,表现在以下2个方面:一是湖泊是唯一一类波谱曲线整体呈下降趋势的地物;二是湖泊在Blue、Green 波段的反射率值接近其他地物,而在Red、NIR、SWIR1、SWIR2 波段明显低于其他地物。基于上述差异,在研究区内,可以根据TOA反射率来区分湖泊水体和其他地物。
对于清水,在Blue-Green波段反射率为4%~5%,波长0.6 μm以上的Red波段反射率降到2%~3%,在NIR波段上几乎吸收全部的入射能量[46]。对于柴达木盆地的湖水来说,随着悬浮物浓度的增加,与清水相比,反射率数值会发生一定的影响[47],但整体的反射率峰值仍然处于Blue、Green 波段范围内,且反射率随波长的增加而减小[48-49]。基于此,对1.2.3小节中的采样结果,取TOA 反射率,执行(Blue+Green-Red-NIR-SWIR1-SWIR2)运算,结果表明,湖泊水体中,所有湖泊样本点的TOA反射率运算结果均大于0,而几乎所有的其他地物样本点表现出相反的运算结果。由此可以初步得出以下结论:(Blue+Green-Red-NIR-SWIR1-SWIR2)运算可以将湖泊水体信息增强为大于0 的值,将其他非水体信息抑制为小于0的值。据此,提出了基于Landsat系列卫星影像的湖泊水体信息提取模型——湖泊水体差分模型(Lake Water Differential Model,LWDM)。
为了验证模型的可行性,对比不同模型的提取效果,本文选择位于柴达木盆地东部的可鲁克湖、尕海和托素湖,中部的大柴旦湖、小柴旦湖,西部的尕斯库勒湖,以及察尔汗盐湖为验证对象,利用不同传感器数据,对LWDM 模型进行验证。同时,为证明模型在更大空间尺度的普遍适用性,还选择了临近柴达木盆地东部边缘的青海湖、尕海湖作为验证对象。
采用LWDM、NDWI和MNDWI方法部分提取结果如下图所示(图4~图5)。图4a 采用真彩色影像合成并进行平方根拉伸,图5a采用真彩色影像合成并进行高斯增强,旨在突出湖泊在原影像上的位置。
图4 可鲁克湖区域湖泊水体提取结果Fig.4 Results of lake water extraction in Keluke Lake
图5 青海湖区域湖泊水体提取结果Fig.5 Results of lake water extraction in Qinghai Lake
在可鲁克湖流域内,由图4 可知,3 种方法均可以将湖泊水体信息提取出来,但NDWI 方法在提取湖泊信息(图4c)的同时提取了河流水体、沼泽湿地和部分阴影的信息,MNDWI 方法同时提取了河流水体、沼泽湿地、阴影和冰雪的信息(图4d);这种情况在小柴旦湖周边区域以及尕斯库勒湖区域同样存在。
在青海湖区域,选择同一时间段(2018年7月),对比野外实测湖岸线与基于LWDM 的湖泊提取结果(图6),可以看出,提取结果与实测湖岸线误差在一个像素内,能够反映实际湖水岸线的曲折情况。
图6 青海湖区域LWDM方法水体局部提取结果Fig.6 Results of lake water extraction by LWDM in Qinghai Lake
虽然,在模型建立之初就已经考虑到了柴达木盆地不同湖泊水体矿化度的差异,作为测试对象的尕海,其矿化度已经达到了卤水的标准,但与察尔汗盐湖矿化度高达400 g·L-1相比,还存在较大的差距,为了进一步验证模型的适用性,运用本文提出的方法,提取了察尔汗盐湖及其周边的湖泊水体信息(图7),由图7 可知,LWDM 较为完整的提取了湖泊水体信息,有效去除了区内的道路、建筑物等干扰信息,总体分类精度为97.98%,Kappa 系数为0.9339,用户精度99.94%,说明LWDM 模型在柴达木盆地内具有较好的适用性。
图7 察尔汗盐湖水体提取结果Fig.7 Results of lake water extraction by LWDM in Chaerhan Salt Lake
通过上述对比可以得出如下结论:LWDM 在提取湖泊水体信息的同时,不仅可以去除植被、裸地、建筑物等干扰信息,也可以较好地去除绝大部分雪、阴影、沼泽湿地和地表河流等干扰因素,与NDWI、MNDWI 等水体提取方法相比,LWDM 方法可以更好地提取柴达木盆地内的湖泊水体信息。
根据原始的混淆矩阵,经过公式(1)~(3)的运算,得到精度评价结果(表2)。
由表2可知,湖泊水体信息的提取精度表现为:LWDM>NDWI>MNDWI。在青海湖区域附近,3 种方法提取湖泊水体信息的精度接近,而在柴达木盆地内的可鲁克湖区域、小柴旦湖区域、尕斯库勒湖区域、察尔汗盐湖区域,3种水体提取方法在湖泊提取精度上表现出显著的差异,其原因是NDWI、MNDWI 在提取湖泊水体时,会同时提取地表河流、冰雪、沼泽湿地等非湖泊水体信息,从而降低了用户精度。
表2 湖泊水体信息提取精度评价结果Tab.2 Evaluation result of information extraction accuracy of lake water
此外,在柴达木盆地的其他区域不同时期、不同传感器共计90个Landsat图像场景数据中,OA达到99.48%,Kappa值达到0.9877。
由此可以得出如下结论:基于Landsat 系列影像,LWDM 可以较好地提取柴达木盆地湖泊水体信息,相对于NDWI、MNDWI等水体信息提取方法,去噪能力强,提取精度显著提高。对于柴达木盆地之外的其他区域,针对不同清洁度的湖泊,LWDM 方法的适用性有待进一步的研究。
以Landsat TM、ETM+和OLI 多光谱遥感影像为数据源,通过分析可鲁克湖流域内湖泊水体与冰雪、阴影、沼泽湿地、裸地、建筑物等对象的TOA 反射率差异,根据湖泊水体在Blue、Green 波段的TOA反射率高于其余波段,且不同波段TOA反射率值呈整体下降趋势的特征,利用蓝光、绿光、红光、近红外与短波红外波段,构建了柴达木盆地湖泊水体信息提取模型——湖泊水体差分模型(LWDM)。将该模型应用于柴达木盆地的可鲁克湖、托素湖、尕海、尕斯库勒湖、小柴旦湖以及盆地附近青海湖地区的尕海湖等湖泊水体信息的提取,与Mcfeeters[19]的NDWI和徐涵秋[21]的MNDWI相比,LWDM不仅可以较为准确地提取出湖泊水体信息,而且可以有效去除冰雪、阴影、沼泽湿地、建筑物等干扰信息。模型提出后,利用Landsat(Level-2 产品)的地表反射率作为输入数据对模型进行进一步测试,结果表明其提取精度与TOA 反射率作为输入数据的提取精度基本一致。因此,本文提出的湖泊水体差分模型既可以将TOA反射率作为输入数据,也可以将地表反射率作为输入数据;同时应用不同时段、不同季节的Landsat 系列影像数据进行了测试,结果表明,该模型所反映的湖泊水体规律在不同年份、不同季节是一致的。此外,NDWI和MNDWI提出的目的是提取地表水体信息[19,21],并未考虑湖泊水体信息与河流水体信息的差异,因此,在提取湖泊水体信息的同时,也混入了地表河流水体信息(图4c~图4d)。而LWDM建立的目标就是提取湖泊水体信息,在前期试验研究过程中充分考虑了湖泊水体与地表水体信息之间的差异,在提取湖泊水体信息的同时,较好地剔除地表河流的水体信息(图4b)。
从提取精度评价的结果来看,LWDM 方法的提取精度要高于NDWI 和MNDWI,原因有三:第一,NDWI是在分析植被与水体在可见光波段和近红外波段的反射强度特征后提出的[19],目的是最大程度地抑制植被信息,增强水体信息。而MNDWI 是针对NDWI不能很好区分城市范围内的水体而改进的一种水体指数[21],目的是增加建筑物和水体的差异,进而更精确的提取城市范围内的水体,但是不能很好地抑制含水量较大的土壤干扰信息。二者以提取地表水体(包括河流、湖泊、水库等)为目标,建立时更多地关注植被、建筑物对水体信息提取的干扰。第二,虽然徐涵秋[21]在提出MNDWI时,选择河流、湖泊、海洋3 种水体对NDWI 和MNDWI 进行了水体信息提取试验,实验结果证明,MNDWI 具有一定的阴影去除效果,但未考虑像柴达木盆地这样复杂的环境:一是盆地属于高原型内陆盆地,极高山、高山、戈壁、丘陵、平原、盐沼等地貌并存;二是区内湖泊、河流密布,沼泽湿地普遍存在,不同湖泊矿化度差异较大;三是盆地周边的极高山、高山区有冰雪覆盖区。这在一定程度上降低了MNDWI模型提取柴达木盆地湖泊水体信息的精度。第三,LWDM 模型以提取柴达木盆地湖泊水体信息为目标,建立时不仅考虑了上述背景条件,而且借助波谱曲线捕获到了湖泊水体与其他地物要素的TOA反射率差异,模型在精确刻画地物差异的同时,提高了提取结果的可靠性。
LWDM 在建立时既考虑了可鲁克湖流域内水体在可见光波段,特别是蓝光与绿光波段反射率相对较强,在近红外与短波红外强吸收的特征,又考虑了DN值在表征地物时失真的不足。模型建立过程中采用差值运算方法并引入较多的波段,一方面可以有效避免由于比值运算造成局部信息丢失的现象;另一方面引入多波段可以探测到更多的水体微细变化信息,对与柴达木盆地相似的不同矿化度湖泊水体信息并存的区域,可取得比NDWI、MNDWI更好的提取效果。在阈值选择方面,有学者采用神经网络技术自动计算样本阈值,得到了水体信息自动提取模型,但该方法容易得到局部极小值而延长计算时间,而本文提出的LWDM模型将提取湖泊水体信息的阈值设定为0,与NDWI 和MNDWI 在局部地区信息提取时需要调整阈值[21]操作及采用神经网络技术导致计算时间较长相比,LWDM 提取湖泊水体信息操作更加便捷。
为了考证LWDM的适用性,选择柴达木盆地以外位于不同地区的天池、太湖、滇池作为验证对象,利用LWDM 对其周边区域的湖泊水体信息进行了提取,借助总体分类精度OA、Kappa 系数及用户精度对提取结果进行了精度评价(表3),由表3 可知,天池的提取结果较好,太湖和滇池的效果稍差,原因是太湖和滇池的水质相对较差,说明LWDM模型对清洁水体的适用性较好,而对受污染水体的适用性相对较差,前期笔者对太湖含蓝藻水体进行了专题研究,提出了基于Landsat影像的含蓝藻水体提取方法(Blue+Green-Red-SWIR1-SWIR2>0)[50],通过对比发现,提取结果优于NDWI、MNDWI、MBWI、NWI、WI2015,后续针对不同清洁度的湖泊水体信息,还需要对模型进行完善,以提高其适用性。利用MODIS数据(成像时间,2005年5月25日)进一步对LWDM 的适用性进行分析,选择盆地内外的哈拉湖、青海湖作为试验对象,对提出的水体提取指数进行了试验验证,采用LWDM提取了2个湖泊水体,提取结果的总体分类精度为99.81%,Kaapa 系数为0.9973,用户精度99.78%,说明对于不同的传感器,在青藏高原地区该方法也有较好的适用性。
表3 滇池、太湖、天池湖泊水体信息提取精度评价结果Tab.3 Evaluation results of water body information extraction accuracy of Dianchi Lake,Taihu Lake and Tianchi Lake
由于湖泊水体信息受制于多种因素,当前的研究方法尚不能完成全面、系统的湖泊水体信息提取,特别是对较差水质的湖泊水体信息,模型还有待进一步完善。本文基于前人的研究成果,对柴达木盆地的湖泊水体分布自动提取方法进行了探索性研究,后续还要在以下2个方面加强研究:一是引入机器学习等人工智能技术,自动识别青藏高原湖泊水体与冰雪、沼泽湿地、阴影等要素的TOA 反射率差异,完善湖泊水体信息提取模型,将其应用到塔里木盆地、云贵高原、四川盆地等不同区域,在应用中不断改进模型,增强模型的普适性;二是考虑SPOT、QuickBird、ASTER、IKONOS 及国产高分辨率对地观测等卫星遥感数据,采用多源遥感影像数据融合技术,关注气象变化与人类活动的影响,提高湖泊水体分布信息的提取精度,进而为湖泊水体动态变化及其周边区域生态环境演化提供理论与技术支撑。
本文基于Landsat系列遥感影像,通过对比分析可鲁克湖流域不同地物要素的TOA反射率差异,建立了基于Landsat 影像的湖泊水体分布自动提取方法——湖泊水体差分模型(LWDM),将该方法应用于柴达木盆地、青海湖、尕海湖、太湖、天池、滇池等湖泊水体信息提取,并分析了同一卫星不同传感器、不同成像时间、不同卫星不同传感器的提取结果。主要结论如下:
(1)在统计分析的基础上,根据地物在遥感影像不同波段的TOA反射率差异特征,通过增强目标地物信息,抑制非目标地物信息,可以有效提取所需信息。
(2)基于TOA反射率差异特征提出的湖泊水体差分模型(LWDM)可以用稳定的阈值(0)实现区域湖泊水体信息的快速提取,与NDWI、MNDWI 等方法相比,在柴达木盆地内,该方法能有效抑制地表河流、冰雪、沼泽湿地、阴影等干扰因素,显著提高湖泊水体信息提取的精度。
(3)基于Landsat 系列卫星影像,LWDM 可用于柴达木盆地湖泊水体分布快速识别,基于不同传感器不同时期的遥感影像,通过对比不同湖泊提取结果,OA达到99.48%,Kappa值达到0.9877,用户精度达到99.66%。基于不同卫星传感器,利用MODIS遥感影像数据,提取了青海湖和哈拉湖水体信息,OA达到99.81%,Kappa 值达到0.9973,用户精度达到99.78%,结果具有较高的可信度,可为柴达木盆地湖泊水体动态监测及湖泊周边生态环境保护提供技术支持。