陕西省黄土高原地表水体遥感监测与时空演变特征分析

李辉, 高学睿, 谢治国, 赵春明, 连迎馨, 王纪超

(1.西北农林科技大学 水利与建筑工程学院,陕西 杨凌 712100; 2.西北农林科技大学 水土保持研究所,陕西 杨凌 712100; 3.陕西省林业科学院湿地与草原研究所,陕西 西安 710082;4.安康市白河县农业技术推广站,陕西 安康 725899)

水是生命之源、生产之要、生态之基,是自然和社会经济系统可持续健康发展的基本要素之一[1]。作为区域水资源量的重要组成部分,地表水资源在流域水循环和水平衡中扮演着重要角色[2-3]。然而,由于地表水直接参与陆-气交互过程,其对环境变化十分敏感。在全球气候变暖与下垫面强人类活动叠加作用影响下,流域地表水资源演变规律发生显著变化,对区域生态系统保护与修复和国民经济高质量发展产生重要影响。因此,开展区域地表水资源动态监测,获取地表水资源时空变化信息对制定水资源保护规划、优化水资源配置格局和提升区域水资源承载力具有十分重要的意义。

地表水体动态监测是气候水文与生态系统变化研究领域的热点问题。及时准确地监测下垫面水体时空变化特征可为生态保护、经济发展和政府决策提供重要的科学依据。遥感技术具有重返周期短、效率高、覆盖面广等特点,为大尺度区域水体变化监测提供了有效的技术手段。目前,许多学者基于多源卫星数据对不同区域水体的分布和时空演变规律进行了监测和分析。梁益同等[4]基于Landsat卫星数据采用特征指数与最大似然法分类相结合的方法对洪湖进行了遥感监测,分析了洪湖水体面积的长系列时变特征及驱动因素;张行清等[5]利用HJ-1卫星数据采用单波段阈值法与目视解译相结合的手段，识别了广西壮族自治区主要水库水面动态变化特征及其成因。在更大尺度上,FENG M等[6]利用Landsat影像生成了2000年全球内陆地表水体数据集;PRIGENT T等[7]利用多源卫星观测数据重建了1993—2007年全球0.25°空间分辨率的月尺度地表水体变化数据集,并分析了水体变化与区域人口之间的关系。已有研究虽然在特定区域或者全球范围内开展了地表水体信息的提取和识别工作,但受影像空间分辨率、卫星重返周期和数据可用时段的影响,在小型内陆地表水体识别的精细度以及水体长时间序列演变规律解译分析方面还存在一些不足。随着遥感技术的快速发展,各种高时空分辨率遥感数据和产品的出现为流域地表水体精细化识别和长时序演变规律解析提供了极大的便利[8]。

当前,黄河流域生态保护与高质量发展已经上升为国家战略[9]。毫无疑问,水资源是黄河流域生态保护和高质量发展的核心要素,也是流域未来发展前景的控制性要素[10]。陕西省黄土高原是黄河流域中部的重要组成部分。近些年,退耕还林(草)工程的实施显著改变了当地的下垫面条件,加之气候变化与社会经济发展等因素的共同影响,使该区地表水体和水资源量发生了显著变化[11]。准确识别和提取陕西省黄土高原地表水体时空分布信息,解析地表水体长时序演变规律及其驱动因子对实施黄河流域生态保护与高质量发展战略,促进该区生态建设和产业可持续发展具有十分重要的意义。

基于此,本文利用GEE(Google Earth Engine)平台对陕西省黄土高原长时间序列(2000—2018年)的Landsat系列遥感影像进行解译和分析。基于多波段谱间关系法的构建思想[12],结合水体指数和植被指数建立了一种简单可靠的地表水体识别算法模型,并对研究区永久性水体和季节性水体进行了识别,揭示了新世纪以来陕西省黄土高原地表水体的时空演变特征。与此同时,为探究研究区地表水体面积变化的主要驱动因子,采用累积量斜率变化分析法量化了自然和人类活动两因素对水体面积动态变化的贡献率。研究结果可为黄河流域生态保护和区域水资源可持续利用提供科学参考。

1 研究区概况

本文以陕西省黄土高原为研究区(图1)。研究区总面积约13.4万km2,其国土面积、人口数量、经济总量分别占陕西省相应总量的65%、76%和87%;空间范围包括榆林市、延安市、宝鸡市、咸阳市、铜川市、西安市、渭南市的整个区域和商洛市的部分区域;区域地势西高东低,按地貌类型可划分为3个大类6个亚类,即陕北黄土高原区(长城沿线风沙滩地区、陕北黄土丘陵沟壑区、陕北黄土高原沟壑区、陕北土石低山区)、关中盆地区(关中黄土台塬区)和秦岭—关山地区(秦岭北坡—关山山地)。研究区基本属于暖温带半干旱气候,仅北部边缘为中温带半干旱气候，年均气温8～12 ℃,年平均降水量150～250 mm,且降水时空分布不均。

图1 研究区范围

2 数据与方法

2.1 数据类型及来源

2.1.1 遥感数据

1)Landsat系列遥感影像数据。本文使用的Landsat系列遥感影像数据包括Landsat5、Landsat7、Landsat8,研究时段为2000—2018年。研究区在Landsat系列遥感影像中的条带号为126～128,列编号为33～36。该系列遥感影像的时间分辨率为16 d,空间分辨率为30 m。Landsat系列遥感影像对应波段的具体信息见表1。

表1 Landsat卫星对应的各波段信息

2)JRC(Joint Research Centre)水体数据集。JRC水体数据集制图的初衷是为了厘清全球地表水的时空变化趋势,以期为水管理部门提供决策依据。JRC水体数据集是由300多万景长时间序列遥感影像生成的[13]。该数据集中的每个像素被分为水体像素和非水体像素,以月为时间尺度记录了两个时期(1984—1999年、2000—2019年)全球地表水体的范围和变化情况。本文利用JRC数据集优化水体提取结果。

3)Sentinel-2MSI遥感数据。Sentinel-2MSI是高分辨率多光谱成像卫星,卫星的重访周期为10 d。其遥感数据的空间分辨率为10 m,精度较Landsat系列卫星的高,且成像时间与Landsat系列卫星的间隔较小。本文在GEE云平台上通过筛选操作获得云量低于10%的高分辨率影像，以此对水体提取结果进行精度验证。

2.1.2 气象数据

本文采用的气象数据来源于中国气象驱动数据集CMFD(China Meteorological Forcing Dataset)。该数据集专门为研究中国陆面过程而制作,是被网格化且具有高时空分辨率的近地表气象数据集[14]。CMFD数据集由遥感产品、再分析数据集和气象站现场观测数据融合而成,时间分辨率为3 h,空间分辨率为0.1°,记录了从1979年1月开始,一直延续到2018年12月的7个近地表气象要素的数据,包括降雨量，距地面2 m处的气温、地面气压、空气湿度，距地面10 m处的风速、向下短波辐射、向下长波辐射。本文从CMFD数据集中获取了研究区2000—2018年的降雨量数据。

2.2 研究方法

2.2.1 数据预处理及地表水体提取

GEE遥感云计算平台是谷歌公司联合其他部门开发的快速处理遥感影像的云平台。该平台包括算法函数区、代码编写区、文件输出区以及地图展示区4个部分。用户在平台上的代码编写区编写代码,形成处理操作命令,通过API接口发送处理操作命令至远程服务器,最后将结果反馈至平台的文件输出区和地图展示区。开发GEE遥感云计算平台的初衷是为了处理海量的遥感数据。其最大的优势在于处理数据的时效性，用户可以根据自身需要编写相关代码快速处理长时间序列、大范围尺度遥感数据。本文利用GEE遥感云计算平台,通过去云算法,完成对Landsat5、Landsat7、Landsat8地表反射率遥感影像数据集的预处理。

1)水体提取规则。水体提取的关键在于增大水体与其他地物的光谱特征差异,差异越明显,识别效果越好。本文以归一化植被指数(Normalized Difference Vegetation Index,NDVI)[15]、改进的归一化水体指数(Modified Normalized Difference Water Index,MNDWI)[16]、增强植被指数(Extreme Value Index,EVI)[17]之间的大小关系构建水体提取模型,完成了对研究区地表水体的提取。即其数值符合MNDWI>EVI或MNDWI>NDVI、EVI<0.1规则的像元为水体像元,否则为非水体像元。该水体提取模型涉及的波段数多,指数间比较关系显著增强了水体与其他地物间的光谱特征差异[18]。归一化植被指数(NDVI)、改进的归一化水体指数(MNDWI)、增强植被指数(EVI)的计算公式分别为:

(1)

(2)

EVI=2.5(ρnir-ρred)-(ρnir+6ρred-7.5ρblue+1)。

(3)

式中ρblue、ρgreen、ρred、ρnir、ρswir1分别为蓝波段、绿波段、红波段、近红外波段、短波红外波段的反射率。

2)不同类型水体的判定。噪声的存在会对水体提取结果产生影响。本文通过计算水体频率去除了由低质量遥感影像产生的噪声,具体计算流程为:①将水体像元赋值为1,非水体像元赋值为0。②计算1年中水体像元数与总体能观测到的良好像元数的比值,将该比值定义为水体频率。③若水体频率大于0.75则定义该像元为永久性水体像元；若水体频率的变化范围为0.25～0.75，则定义该像元为季节性水体像元；若水体频率小于0.25，则定义该像元为非水体像元[19-21]。水体频率Wfreq的计算公式为:

(4)

式中:S为观测到的像元总数;mk为第k个像元的属性。当该像元为水体像元时，mk=1；当该像元为非水体像元时，mk=0。

2.2.2 混淆矩阵参数计算

本文采用混淆矩阵对研究区地表水体提取结果进行精度验证。具体的指标包括生产者精度(Pay)、用户精度(Uay)、总体精度(Oay)和Kappa系数(Kc)4个参数[22]。其中,生产者精度(Pay)表示在分类中,水体类别的地面真实参考数据被正确分类的概率;用户精度(Uay)表示分类器将像元归类为水体像元时,参考数据(Sentinel-2)的真实类别为水体的概率;总体精度(Oay)指所有被正确分类的样本点占总样本点数的百分比;Kappa系数(Kc)考虑的是分类结果和参考数据间的一致性,以及取样结果和分类结果的一致性,一般而言,Kappa系数的变化范围为0～1,其值越大表示分类精度越高。生产者精度(Pay)、用户精度(Uay)、总体精度(Oay)和Kappa系数(Kc)的计算公式分别为:

(5)

(6)

(7)

(8)

式中:Sij为第i行第j列样本点的观测值;Si为第i行的样本点总数;Sj为第j列的样本点总数;Stotal为正确分类的样本点总数;n为地面验证样本点总数;r为行数。

2.2.3 贡献率计算

本文采用累积量斜率变化分析法识别自然和人类活动因素对研究区地表水体面积变化的影响。首先,根据累积距平法得到拐点年份。然后,根据拐点年份左右时间段的年累积量斜率变化率,计算影响因素的贡献率值。贡献率的计算流程如图2所示,具体的计算步骤如下:

步骤1 计算累积距平值,生成累积距平曲线。在因变量y随时间变化的序列值中,第t个时段的累积距平值Dt可以表示为:

(9)

式中:yi为第i个时段的因变量值；y为因变量在一定时间范围内的平均值。同理可以得到自变量x的累积距平曲线。累积距平值的变化反映了一系列离散数据相对于其均值的大小关系，若累积距平值增大,表明离散数据大于其均值,反之则小于其均值。通过累积距平曲线可以直观地看到一系列累积距平值随时间的变化趋势,并得到变化趋势突变的拐点[23]。

步骤2 计算累积量斜率变化率。累积距平曲线的拐点将因变量累积量划分为变化前后两个时期,分别计算这两个时期因变量累积量拟合直线的斜率ka和kb,由此计算出因变量y的累积量斜率变化率Rkw。具体计算公式为:

(10)

同理，可以计算出自变量x的累积量斜率变化率Rkp。

步骤3 计算贡献率。由Rkp与Rkw的比值得到自变量x对因变量y变化的贡献率Cp。若因变量受多个自变量的影响,则除x外的其他自变量对因变量y的贡献率CH可以表示为:

CH=1-Cp。

(11)

图2 贡献率计算流程图

3 结果分析

3.1 地表水体提取精度分析

由于Sentinel-2MSI遥感影像与Landsat系列遥感影像的成像日期接近,误差在2 d以内,且分辨率(10 m)较高,因此以Sentinel-2MSI遥感影像为地面验证数据,验证陕西省黄土高原的地表水体提取精度。从研究区2015—2018年的地表水体提取结果中,随机选取水体和非水体样本点各500个，并以2015—2018年6—10月的Sentinel-2MSI遥感影像(云覆盖量小于15%)为地面验证数据验证样本点的真实性,采用混淆矩阵表征精度验证结果,详见表2。由表2可知,地表水体的生产者精度(Pay)为96.89%,用户精度(Uay)为93.60%,总体的提取精度(Oay)为95.30%。经计算，一致性检验系数(Kc)为88.20%。这说明地表水提取产品的精度较高,可作为后续研究的基础。

表2 研究区地表水体精度验证混淆矩阵

3.2 陕西省黄土高原地表水体的空间分布

2000—2018年陕西省黄土高原地表水体的空间分布如图3所示。由图3可知:本文采用的水体提取模型对内陆湖泊(红碱淖湿地)、河流(黄河和无定河)、水库(马家庄水库)等类型水体的提取效果良好。

(a)红碱淖湿地;(b)无定河;(c)黄河;(d)马家山水库

图3 2000—2018年陕西省黄土高原地表水空间分布

2000—2018年研究区所辖各市的多年平均地表水体面积数据见表3。由表3可知:①研究区地表水体面积的多年平均值为684.07 km2,空间上看,分布较不均匀,主要集中在关中的渭南市和陕北的榆林市。其中,渭南市地表水体面积占研究区总水体面积的29.7%,榆林市地表水体面积占研究区总水体面积的39.6%。该两市在研究区地表水体中的占比接近70%。因此,在制定水资源保护措施时,应优先考虑榆林市和渭南市。②2000—2018年研究区地表水体中季节性水体和永久性水体的面积基本相当。其中,永久性地表水体面积为335.97 km2,占地表水总水体面积的49.1%;季节性地表水体面积为348.10 km2,占地表水总水体面积的50.9%。③渭南市季节性地表水体分布最多,面积为126.17 km2,是其永久性地表水体面积的1.63倍;榆林市永久性地表水体分布最多,面积为149.48 km2,是其季节性地表水体面积的1.23倍。

表3 2000—2018年研究区所辖各市的多年平均地表水体面积 km2

为了分析研究区地表水体面积的变化趋势,本文以市为单元,线性拟合了2000—2018地表水体面积的变化情况,得到各市的地表水体面积线性拟合方程。对地表水体面积和年份进行了相关性分析和显著性检验，结果见表4。由表4可知:2000—2018年,榆林、铜川的地表水体面积随时间呈减少趋势；宝鸡、咸阳、西安、延安、渭南、商洛的地表水体面积呈增加趋势;显著性检验发现,2000—2018年,宝鸡、咸阳和西安地表水体面积呈显著上升趋势(p值均小于0.05),其原因可能与当地水库工程、湖泊湿地建设有关。

表4 2000—2018年研究区所辖各市地表水体线性拟合方程

3.3 陕西省黄土高原地表水体随时间的变化

2000—2018年研究区的季节性水体面积、永久性水体面积和总水体面积的年际变化如图4所示。由图4可知:①陕西省黄土高原总水体面积呈增加趋势,由636.83 km2增加至760.43 km2,增长率为19.4%。其中,研究区2000年的总水体面积最小,为636.83 km2；2018年的最大,为760.43 km2。②地表水体总面积变化最大的时期为2015—2016年,其值由2015年的648.53 km2增加到2016年的721.44 km2,变化值达72.91 km2。分析其原因为,2015年受厄尔尼诺现象的影响,我国华北、黄淮、西北地区东部出现了不同程度的旱情[24],而2016年陕西发生了400余次洪涝、滑坡等自然灾害。这表明降雨等气象因子会对研究区地表水体的面积变化造成影响,同时也验证了本文方法的适用性。③永久性水体面积呈增加趋势。其中,2014年永久性水体面积最大,为428.44 km2；2001年永久性水体面积最小,为249.08 km2。④季节性水体面积呈递减趋势。其中,2005年季节性水体面积最大,为428.44 km2；2015年季节性水体面积最小,为266.69 km2。

图4 2000—2018年陕西省黄土高原地表水体的面积变化

3.4 贡献率分析

已有研究结果表明，地表水体面积的时空分布主要受自然因素和人类活动的影响[25]。其中,降雨量和蒸散发量是影响地表水体面积变化的主要自然因素。由孙淼[26]的研究结果可知,1990—2018年黄土高原年均蒸散发量变化较小,因此在分析研究区地表水体面积变化的自然因素时主要考虑降雨量。本文在利用累积量斜率变化分析法计算贡献率前,分析了研究区年降雨量和地表水体面积的变化情况(图5)，并通过皮尔逊相关分析得到地表水体面积和年降雨量的相关系数，其值较高(R2=0.716),表明水体面积与当地降雨量密切相关。

图5 陕西省黄土高原地表水体面积和年降雨量的变化情况

由图5可知:地表水体面积的变化情况和年降雨量的变化情况一致,都呈现出增大的趋势。通过计算年均降雨量对研究区地表水体面积变化的贡献率量化其对地表水体的影响。根据累积距平法得到研究区2000—2018年年降雨量和地表水体面积累积距平特征,如图6所示。由图6可知:①年降雨量累积距平曲线只有一个拐点年份(2010年)。2010年前后年降雨量的累积距平出现先减少后增大的趋势。②地表水体面积累积距平曲线有3个拐点年份(2003年、2010年、2015年)。地表水体面积累积距平曲线在2003年前呈减少趋势,2004—2010年呈缓慢增加的变化趋势,2011—2015年又呈缓慢减少趋势,2015年以后迅速增加。采用累积量斜率变化法计算研究区降雨量对地表水体面积变化的贡献率。线性拟合陕西省黄土高原2000—2010年、2011—2018年2个时期的累积年降雨量和2000—2003年、2004—2010年、2011—2015年和2016—2018年4个时期的累积地表水体面积。由线性拟合方程计算累积年降雨量和累计地表水体面积的斜率变化率，结果见表5和表6。由表5和表6可知:年降雨量和地表水体面积的累积距平斜率变化率分别为5.2%、19.4%。由此得到研究区自然因素对地表水体面积变化的贡献率为26.8%,人类活动的贡献率为73.2%。

图6 2000—2018年陕西省黄土高原年降雨量、地表水体面积累计距平变化特征

表5 陕西省黄土高原2000—2018年累计降雨量的斜率变化率

表6 陕西省黄土高原2000—2018年累积地表水体面积的斜率变化率

4 结语

本文以陕西省黄土高原为研究区,基于GEE遥感云平台通过水体提取模型提取了研究区2000—2018年地表水体的分布信息,采用混淆矩阵验证了地表水体的提取精度，分析了地表水体的时空分布特征,此外利用累积量斜率变化分析法识别了研究区地表水体面积变化的主要影响因素。得到如下结论:

1)利用强植被指数(EVI)、归一化植被指数(NDVI)与改进的归一化水体指数(MNDWI)构建的水体提取模型能够实现对陕西省黄土高原地表水体的准确提取,提取精度达95.30%。

2)研究区的地表水体主要分布在榆林市(271.05 km2)和渭南市(203.42 km2),且研究区的地表水体面积呈逐年增加趋势。

3)人类活动为研究区地表水体面积变化的主要影响因素。根据贡献率的计算结果可知,自然因素对研究区地表水体面积增加的贡献率为26.8%,人类活动的贡献率为73.2%。因此，在进行水资源保护和规划时,应重点考虑人类活动的影响,如合理规划区域水利工程建设等。