红崖山水库近20年面积变化遥感调查及驱动力分析

郝固状，甘甫平，闫柏琨，李贤庆，胡辉东

(1.中国矿业大学(北京)煤炭资源与安全开采国家重点实验室，北京 100083；2.中国矿业大学(北京)地球科学与测绘工程学院，北京 100083；3.中国自然资源航空物探遥感中心，北京 100083)

0 引言

水资源是人类生存和社会发展不可缺少的生命之源和物质基础，而水库作为水资源中不可或缺的淡水资源，为人类的生产、生活提供了重要的保障。水库的分布状况在一定程度上反映了水资源的区域变化和时空差异。水库水量变化在工业、农业、城市以及防汛抗旱等方面发挥着重要的调节作用[1]，而水库的水域面积变化直接反映着水库水量的增减状况。红崖山水库是位于石羊河流域的大型水库，是上游冰川融水、降水水源补给的主要汇集地与下游工、农、生活用水的主要补给与调节地，其水量、水面分布变化受气候、人为因素的双重作用与影响。

随着遥感技术的不断成熟，充分发挥遥感技术宏观性、动态性和实时性的优势，通过遥感技术可以提供大范围、长时间序列的动态监测，可以节省大量的人力、物力、财力，并快速准确地获取湖泊和水库面积信息。李均力等[2]通过对柴窝湖面积近50 a的时序变化进行分析，得出面积变化分为3个阶段，面积呈现先缓慢增加后快速退缩的趋势；魏善蓉等[3]利用MODIS09卫星数据通过对比6种水体指数法的湖泊提取结果，对柴达木盆地连续14 a的湖泊面积进行解译，湖泊呈现先增加后减少的趋势变化；张文等[4]利用Landsat和GF影像对鄱阳湖进行了长时间序列监测，研究结果表明在空间上，鄱阳湖水体面积呈缓慢收缩的趋势，在时间上，除冬季面积有明显下降趋势外，其他季节整体趋势变化不大，不同季节鄱阳湖的湖面范围有明显的差异，湖泊面积具有显著的季节性；张兆鹏等[5]以95景多时相Landsat卫星影像为主，还原了1987—2016年年际变化(以8月水面分布为对比基准)及部分年内月度变化，并定性分析了水量变化的自然、人为因素影响，认为年际水面变化呈“W”形，总体为增加趋势，人工增雨增雪、农业灌溉、水库建设是主要的影响因素。该研究对红崖山水库的水面面积变化及驱动力得出了基本与初步结论。

本文利用谷歌地球引擎(Google Earth Engine，GEE)平台进行水域面积的提取，这是一种更为便捷快速的研究方法；利用更高分辨率的遥感影像来验证真实的水库水域面积变化，并进行精度对比，提取出更加真实的水域面积；利用气象数据，综合气温、降水、蒸散发的数据以及多角度影响因素对变化原因进行分析。在全球气候显著变化的背景下，还原监测其水量、水面历史分布变化，分析解剖变化的驱动力，对于研究我国西北内陆河流域和水库水资源变化、保障、气候变化影响与适应对策等具有重要的代表意义。

1 研究区概况与数据源

1.1 研究区概况

亚洲最大的人工沙漠水库——红崖山水库位于西北地区甘肃省民勤县，被腾格里和巴丹吉林两大沙漠包围，属于典型沙漠平原水库,因位于石羊河下游，其水源补给主要来自石羊河流域(图1)，红崖山水库是沙漠地区的一座中型洼地蓄水工程，也是唯一的水利调蓄工程，西面依红崖山而建，其他三面都是人工所筑，而且又修建在沙漠中，这在全国甚至全世界都是罕见的[6-7]。

图1 研究区地理位置Fig.1 Location of the study area

1.2 数据源

1.2.1 遥感数据及预处理

本研究主要采用的是Landsat系列影像数据(数据来源：https://developers.google.com/earth-engine/datasets/catalog/landsat)和高分二号(GF-2)卫星影像(数据来源：http://geocloud.cgs.gov.cn)，其中Landsat包括Landsat5 TM，Landsat7 ETM+和Landsat8 OLI数据，影像获取的时间是2000—2019年，空间分辨率为30 m，其中Landsat5有254景，Landsat7有34景，Landsat8有171景(表1)，GF-2包括多光谱(PMS)和全色(PAN)，分辨率分别为4 m和1 m，共12景(表2)。本文对Landsat5/7/8的表面反射率(surface reflectance，SR)数据进行了几何纠正、云雪去除、影像拼接、影像融合等预处理，为后续水体的提取奠定了基础；同时对GF-2数据进行了辐射定标、大气校正、正射校正、影像配准、影像融合和影像裁剪等预处理。本文采用Landsat影像数据与采集时间基本一致的6景GF-2卫星遥感影像作水体提取精度对比。

表1 Landsat遥感影像数据Tab.1 Landsat remote sensing image data

表2 高分二号(GF-2)遥感影像数据Tab.2 GF-2 remote sensing image data

1.2.2 气象数据

获取研究区周边5个气象站点的2000—2018年的气温、降水数据(数据来源：中国气象数据网http://data.cma.cn)。本次研究的气象站点的数量有限，为获取区域内各点的气温降水值，采用克里金插值法得到气象要素栅格图[8]，并提取红崖山水库质心处气温降水数据。

2 研究方法

本文红崖山水库水域面积的提取方法，精度验证以及气象因子、植被覆盖度和人工自然等其他因素对水域面积变化驱动力的分析流程如图2所示。

图2 总体技术路线Fig.2 Overall technical route

2.1 水体指数法提取水体

基于遥感影像数据的水体提取方法主要包括光谱特性提取和纹理特征提取2类。基于光谱信息的提取方法有单波段阈值法[9]、光谱指数法、分类器法；基于纹理信息的提取方法主要是水体特殊的纹理。本文Landsat遥感数据主要对改进归一化差异水体指数(modified normalized difference water index，MNDWI)、归一化差异水体指数(normalized difference water index，NDWI)、增强型水体指数(enhance water index，EWI)、新水体指数(new water index，NWI)4种光谱指数法进行对比分析。

运用水体指数方法的关键在于阈值的选择，本文选用Ostu算法[10]进行阈值的自动分割。Ostu是一种确定图像分割阈值的算法，按其求得的阈值进行图像二值化分割后，前景与背景图像的类间方差最大；其被认为是图像分割中阈值选取的最佳算法，计算简单，不受图像亮度和对比度的影响，因此在数字图像处理上得到了广泛的应用。

1)McFeeters[9]提出的NDWI主要是根据水体的光谱特性利用水体在绿波段(Green)和近红外波段(NIR)的反射特性，使得影像中的水体和其他地物区分开来，同时抑制植被和土壤信息。其计算公式如下：

(1)

式中：Green为TM/ETM+/OLI遥感数据中绿光波段数据；NIR为近红外波段数据。

2)MNDWI[11]水体指数是对NDWI指数的进一步修正，削弱了提取水体时由于山体阴影的影响，可有效抑制建筑物信息，能够更好地提取城市范围内的水体信息。其计算公式如下：

(2)

式中：MIR为TM/ETM/OLI遥感数据中红外波段数据。

3)在GEE平台使用NWI[12]指数时，Landsat数据的表面反射率(surface reflectance,SR)和大气表观反射率(top of atmosphere,TOA)产品的水体提取效果是不同的，经过试验，NWI(TOA)只要大于0就可以很好地提取出水体，而NWI(SR)则要自己调整阈值，根据普适性，NWI(SR)的组合更能满足多地区的水体提取。其计算公式如下：

(3)

式中：B1,B4,B5,B7分别对应TM/ETM+遥感数据中的第1,4,5,7波段数据(OLI对应遥感数据中的第2,5,6,7波段)；C为常数，本文中取200，也可取100或其他正值。

4)EWI[13]主要是将MNDWI和NDWI相结合，利用绿波段，红外波段和短波红外波段进行水体提取，其计算公式如下：

(4)

2.2 植被覆盖度的计算

归一化植被指数(normalized difference vegetation index,NDVI)计算公式如下：

NDVI=(NIR-R)/(NIR+R)。

(5)

NDVI与植被覆盖度之间有紧密相关性，综合大量研究，把NDVI与像元二分模型结合起来，建立了基于NDVI的植被覆盖度估算模型[14]，计算公式如下：

(6)

式中：NDVI为混合像元的值；NDVIveg为纯植被像元的值；NDVIsoil为纯非植被像元的值，理论上NDVIveg=1，NDVIsoil=0。

由于在现实情况中，植被覆盖受到多种气象因子的影响，一般不可能取到理论值，又因为缺乏大面积地表实测数据，经大量研究取定置信区间内的NDVI最大值(NDVImax)和最小值(NDVImin)来代替[15]，因此公式变为：

(6)

本文通过此种方法，将2000—2019年中偶数年份的NDVI值进行了统计，根据累积百分比选择累积百分数2%～98%为置信区间，累积百分比小于2%为近似纯非植被覆盖，小于98%为纯植被覆盖，分别将对应的NDVI值代入NDVImin和NDVImax进行计算[16]。

3 结果分析

3.1 GF-2/Landsat影像水体面积提取的精度验证

由于Landsat数据为低分辨率遥感影像，为保证提取的水体精度的可信度，故采用空间分辨率为1 m的GF-2影像数据作为精度验证数据。通过4种水体指数与GF-2影像目视解译法提取的水域面积进行对比分析(表3)，可知MDNWI水体指数法提取的水体面积精度最高，绝对误差仅为0.006 km2，相对误差0.029%，结合4种水体指数提取的效果图，MNDWI指数(图3(b))提取的水体更加符合原始的遥感影像(图3(a))中的水体，边界上的提取更加全面真实。NDWI指数(图3(c))和NWI指数(图3(d))在石羊河下游与水库入口的连接处的水体未提取出来，EWI指数(图3(e))提取效果总体上很符合真实情况，但是在水库边界的提取上有部分水体无法提取，因此综合对比分析4种指数(图3(f))的提取效果，本文将采用MNDWI指数法作为水体提取的方法。

(a)Landsat8原始影像 (b)MNDWI水体指数法 (c)NDWI水体指数法

表3 4种水体指数法与GF-2目视解译法水体面积提取精度对比Tab.3 Comparison of water area extraction accuracy between four water index methods and GF-2 visual interpretation method

由于GF-2卫星影像数据缺少中红外波段，采用GF-2 PMS传感器获取的4 m多光谱数据和1 m全色数据进行图像融合，得到1 m空间分辨率的多光谱数据，由于波段限制，GF-2影像采用NDWI水体指数法进行水体的提取(图4)。根据表4和图5的精度对比，Landsat影像数据和GF-2影像数据的相关系数高到0.99，相对误差最大为6.63%，最小仅为0.33%，绝对误差最大不超过1 km2，明显看出不同分辨率的遥感影像数据对红崖山水库的水域面积的提取影响不大，因此Landsat影像数据对水域面积的提取可以反映红崖山水库的真实性。

(a)GF-2原始影像 (b)GF-2二值影像 (c)GF-2水体和非水体

表4 GF-2/Landsat数据水域面积提取精度对比Tab.4 Comparison of precision of water area extraction in GF-2/Landsat data

图5 GF-2/Landsat水域面积相关性Fig.5 GF-2/Landsat water area correlation

3.2 水库水域面积变化特征

3.2.1 月际面积变化特征

水库的水量和湖泊的水量变化十分相似，受季节性影响也存在枯水期和丰水期，本文选取每4年间隔(2000年、2004年、2008年、2012年、2016年、2019年)的年份进行月际之间的对比分析。综合分析6 a的柱状图和趋势线(图6)，结果表明：2000年、2004年、2008年、2016年和2019年5 a的水域面积月际趋势呈倒“正态曲线”分布，且每年出现两次峰值，主要集中在3月份和9—10月份，最低值出现在6月份，根据丰水期和枯水期的划分，2000年、2004年、2008年、2016年和2019年的丰水期为春秋季节的3月份和9—10月份，枯水期为夏季6月份；但2012年较为特殊，主要表现在枯水期为4月份和7月份，丰水期仍然是春秋季节的3月份和9—10月份。总体上分析2000—2019年月均水域面积变化(图7)，月际面积变化呈倒“正态曲线”分布，且丰水期主要集中在春秋季节的3月份和9—10月份，枯水期主要集中在夏季6月份。

(a)2000年 (b)2004年 (c)2008年

图7 2000—2019年红崖山水库月均水域面积Fig.7 Average monthly water area of Hongyashan Reservoir from 2000 to 2019

3.2.2 年际面积变化特征

由于水库不同季节的水量有很大变化，本研究选取不同季节年际之间的水域面积变化进行分析。春季选取4月份，夏季选取7月份，秋季选取10月份，冬季选取12月份进行年际之间水域面积的比较。根据不同季节年际变化情况的研究(图8)，结果表明：春夏秋冬四季水域面积都是逐年上升的，春秋季节平均年际增长率为5.03%和5.22%，秋季平均年际增长率最低仅为2.42%，夏季水域面积平均年际增长率为22.19%，年际变化幅度最大，秋季变化幅度最小，因此夏季水域面积变化最不稳定，秋季水域面积变化最稳定。夏季水域面积呈“V”字波动上升且变化幅度很大。

(a)春季2000—2019年面积变化 (b)夏季2000—2019年面积变化

3.2.3 水库水域面积变化总体特征

通过上述年际变化得出夏季的增长趋势与2000—2019年的总体上(图9)的增长趋势极为相似，也呈现“V”形波动上升。根据图9的研究结果表明，总体上水域面积年际增长的相关系数达到0.78，相关性良好；2000—2019年水域总面积增加8.98 km2，面积变化率高达42.6%，水域面积经历了先增加再减少，后持续波动上升的过程，具体表现为：2000—2001年水域面积减少，且在2001年达到了最低值10.7 km2。2001—2003年水域面积增长迅速，面积变化速率为27.7%，为近20 a水域面积增长幅度最大。2003—2004年水域面积骤降4.33 km2，为2000—2019年降幅最大。2004—2017年水域面积呈“V”字波动上升，总体面积增加6.97 km2，增长率为56.6%，其中2004—2005年,2014—2015年为小幅增长；2005—2006年水域面积变化速率为27.4%，接近最大变化速率；2006—2010年，2010—2012年，2012年—2014年,2015—2017年水域面积分别为3个“V”形增长；2017—2019年水域面积近乎直线增长，最接近趋势线走势。根据年际变化率(图10)可以看出年际之间的变化幅度在逐年下降且趋向于0，说明红崖山水库的水域面积趋向于稳定增长趋势。

图9 红崖山水库年均水域面积变化Fig.9 Annual average water area change of Hongyashan Reservoir

图10 红崖山水库水域面积年际变化率Fig.10 Interannual change rate of water area of Hongyashan Reservoir

3.3 水域面积变化驱动力

3.3.1 气象因子的影响

红崖山水库位于西北干旱地区甘肃省民勤县境内，是最大的人工沙漠水库，选取民勤县周边5个气象站点的气温和降水数据进行研究，分析影响水域面积变化的驱动因子(图11)。

1)气温。红崖山水库的主要水源是祁连山东部的冰雪融水[7]，在降水量不变的情况下，气温的升高会导致冰雪融化，入库水量增加，水域面积增大。根据对近20 a红崖山水库的年均气温(图11(a))研究，结果表明：2000—2005年均低于距平气温6.72 ℃，2006年气温骤升至7.17 ℃，水域面积增加了3.43 km2；2006—2008年气温连续下降了0.89 ℃，此时水域面积也在持续下降，2008—2009年气温上升0.8 ℃，2009—2012年气温降至最低气温6.01 ℃，2013年气温达到最高7.26 ℃，2014—2018气温波动上升，年均气温变化总体来说，研究区气温呈上升趋势，并通过气温和水域面积相关性分析(图11(b))得出，气温对水域面积的增加起到正向作用，因此气温上升是水域面积增加的重要的间接驱动因素。

(a)年平均气温距平变化 (b)水域面积-年均气温相关性

2)降水。西北干旱区降水稀少且分布不均，但却是水库水量增加的最直接的补给水源。根据近20年降水量的年际变化(图11(c))表明：2000—2006年(除2003年)降水量低于距平降水量，2003年、2007年和2018年降水量明显高于相邻年份，其对应的年份水域面积也明显增加，2008—2012年降水量波动增加，2012—2015年出现2次明显下降，一次上升，同年的水域面积也发生同样的变化，2015—2018年降水量连续4年上升，总体上看，降水量整体呈明显增加态势，通过水域面积与降水量进行的相关性分析，相关系数高达0.93(图11(d))，进一步说明二者有很好的相关性，因此降水量的增加是影响水域面积增加的主要驱动因素。

3.3.2 植被盖度的影响

植被不仅对区域气候变化有着重要的影响，植被面积的增加可以减少地表对太阳短波辐射的反射率，缓和地层增温，进而减少地气湍流交换和水分蒸发，同时植被恢复也有利于地表水源涵养功能的增加，减少降水的无效下渗和散失，有助于径流的产生，增加水量。通过对研究区周边的植被指数提取，根据年际之间植被覆盖度(图12)显示，总体上植被覆盖度在逐渐增加，2000—2002年植被盖度增加了10%，但是2002—2008年出现持续下降15%，2008年达到最低值，2008—2018年植被盖度逐年上升，增加26%。通过水域面积和植被覆盖度的相关性分析(图13)，表明植被覆盖度对水域面积的增加起了正向作用，植被盖度的增加是导致水域面积增加的重要的间接驱动力。

图12 植被覆盖度的年际变化Fig.12 Interannual changes in vegetation coverage

图13 水域面积-植被覆盖度相关性Fig.13 Water area-vegetation coverage correlation

3.3.3 其他因素的影响

1)加高扩建工程。为缓解红崖山水库防洪与兴利，来水与需水的主要矛盾，实施了红崖山水库加高扩建工程，库区水域总面积约22.0 km2，和通过遥感影像获取的最大面积基本吻合，同时使得总库容从0.993亿m3增加到1.48亿m3[17-18]，不仅使得西北干旱地区的水源得到充足使用，而且提高了水库蓄水的调节能力，加强了农业灌溉的力度，这是红崖山水库水域面积逐渐增大的主要驱动力。

2)生态环境修复。为了石羊河流域的综合治理，恢复青土湖生态，红崖山水库向青土湖的下泄生态水量逐年增大[19]，使得青土湖水量从无到有，水位有所上升，地表植被覆盖逐渐增加，这是红崖山水库水量逐年增加的重要原因。每年夏天农业灌溉结束之后，红崖山水库将在8—11月份有计划地向青土湖生态补水3 180 m3[20]，这是导致每年11月份红崖山水库水域面积变小的重要原因，因需求量逐渐增加，红崖山水库的水量也在不断增加，这是导致红崖山水库逐渐上升的间接驱动力。

3)入库径流变化。1962—2010年均入库径流最小值为2001年的0.70亿m3，据金彦兆等[21]研究，红崖山水库站径流量年际变化曲线表明，2001—2010年入库径流量逐渐增加，根据红崖山水库站季节径流量变化图发现，春夏秋冬四季的径流量自2000年以后都出现回升，近20 a红崖山水库水域面积不断增加的直接原因是入库径流量自2000年逐渐增加。

4)工、农、生活用水。红崖山水库的最重要的职能是供工业园区、农业灌溉和百姓生活用水。民勤县红沙岗工业园区的需水主要来自红崖山水库，以满足工业用水，保证工业生产顺利进行，为民勤县的经济发展提供便利[22]，为当地的国民经济打下坚实基础。每年的夏季6—8月份为农业、林业灌溉期，气温升高，蒸散发增加，导致该时期红崖山水库出现枯水期，红崖山担负了红崖山灌区13个乡镇，2个国营农林场农田及生态经济林草灌溉任务,实际灌溉面积440.4 km2，是主要的灌溉水来源。由于工农业用水的增加导致水域面积处于波动上升的趋势。2010年全县总人口比2000年减少近2.8万人，人口总量与用水呈正相关，人口的减少意味着用水量的减少，这是水域面积增加的间接驱动力。

4 结论

本文通过对西北地区红崖山水库2000—2019年的水面面积变化遥感调查及驱动力分析，得出以下结论：

1)与传统的遥感监测方法相比，利用GEE云计算平台可以在线处理多种类型的遥感影像，节省存储和运行空间，不受时空的制约，处理速度得到数量级提升，效率得到指数提高，为后续的大尺度区域的遥感影像处理提供了有力支撑。

2)红崖山水库的水域总面积增加8.98 km2，面积变化率高达42.6%，呈波动上升的趋势，年际变化率逐年减小，水域面积整体趋于稳定增加状态。2000—2019年水库水域面积月际变化呈倒“正态曲线”，丰水期主要集中在每年春秋季节的3月份和9—10月份，枯水期主要集中在每年夏季的6月份；春夏秋冬四季年际水域面积都是逐年上升，春秋季节平均年际增长率为5.03%和5.22%，秋季平均年际增长率最低仅为2.42%，夏季水域面积平均年际增长率为22.19%，年际变化幅度最大，秋季变化幅度最小，因此夏季水域面积变化最不稳定，秋季水域面积变化最稳定，夏季水域面积呈“V”字波动上升且变化幅度很大。

3)根据对红崖山水库水域面积增加的驱动力分析，直接驱动因素为降水量的变化，加高扩建工程，以及入库径流变化，其中降水量的增加和加高扩建工程是水域面积增加的主要驱动力；间接驱动因素为气温的变化，植被覆盖度的变化，生态环境修复以及工、农、生活用水，都是十分重要的驱动力。