基于GEE平台的1990—2019年松涛水库水体面积时空变化

陆天启， 邵长高*， 任旭光， 马洪波， 王陆一

(1.广州海洋地质调查局自然资源部海底矿产资源重点实验室， 广州 510075； 2.广州海洋地质调查局三亚南海地质研究所， 三亚 572025; 3.吉林省气象台， 长春 130062)

水库是地表水资源的主要载体，是陆地水圈的重要组成部分，在淡水供用和水资源调控过程中发挥着不可或缺的作用[1-3]。水库的功能多样，其不仅具有调节区域水资源和气候的作用，还可用来灌溉、发电、养殖、旅游，为人类的生活、生产活动提供水源保障[4-5]。水库筑成后，在自然因素的影响下不断地演变，气候变化会导致水库面积发生不同程度的萎缩或扩张，直接影响着人类社会的生存与发展。水库萎缩会导致一系列生态环境问题的出现，如水质恶化、生态系统失衡、旱涝灾害加剧等[6-7]，进而影响区域生产、生活用水安全，对社会经济发展带来一定的危害。

受气候因素的影响，水库蓄水量具有短期快速波动、长期演变的特征，而对水库面积长期的动态监测是研究水库蓄水量的重点[8-9]。对水库水体面积的动态监测有助于了解水库演变和自然气候之间的相关性，以及评估人类活动对生态环境影响的强度，同时也是进行水资源可持续管理的必要条件。遥感技术具有大面积、动态观测等多方面的优势，自20世纪70年代以来，已被逐渐地应用在水体识别和监测领域[10]。基于水体的光谱特征，构建的归一化差异水体指数(normalized difference water index，NDWI)[11]、改进归一化差异水体指数(modified normalized difference water index，MNDWI)[12-13]、自动提水指数(automated water extraction index，AWEI)[6]、穗帽湿度变换指数(tasseled cap wetness，TCW)[14-15]等水体指数已被广泛地应用在水体提取研究中。针对不同传感器和研究区的特征，中外学者已应用不同的水体指数对不同类型的水体进行了提取分析。基于Landsat TM和ETM+数据，Ouma等[16]利用构建的新型水体指数对肯尼亚5个裂谷湖的变化进行了定量分析，Sarp等[17]使用多种水体指数对土耳其布尔杜尔湖3个年份的时空变化进行了评价分析。利用MODIS数据，万华伟等[18]采用NDWI提取了2000—2013年呼伦湖的水体，在此基础上对其动态变化和驱动力进行了分析。基于Landsat TM和ALOS数据，李志刚等[19]在比较不同水体指数提取海域面积精度的基础上，选用MNDWI研究了1988、1998、2006年锦州湾海域变化情况。

水体面积是衡量水库储水量的一个重要指标[20-21]。联合国可持续发展目标六“确保安全的饮水和卫生设施”明确表示将“与水有关的生态系统空间变化”作为重要指标之一[22]。在研究水体面积长时间序列变化的过程中，GEE(Google Earth Engine)平台具有无可比拟的优势。GEE是由谷歌公司推出的针对地球科学数据进行在线可视化计算分析处理的平台(https://code.earthengine.google. com)，它是由谷歌公司、美国地质调查局和卡内基梅隆大学共同研发，包含了40多年的历史遥感影像、气候、高程等数据集，并且每日进行更新和扩展，其强大的遥感数据集成和计算能力，为进行大范围、长时间序列的遥感分析处理提供了新途径[23-24]。围绕着水体面积这一指标，以海南省松涛水库为研究区，基于GEE云计算平台，选取了Landsat系列遥感影像数据，提取1990—2019年研究区的水体并计算其面积，分析了水体面积变化的时空特征，在此基础上讨论了降水和温度对水体面积变化的影响。

1 研究区及数据源

1.1 研究区概况

松涛水库位于海南省儋州市南部，库区地跨儋州市、白沙黎族自治县两市县，其经纬度范围为19°12′～19°25′N，109°27′～109°41′E(图1)。水库所属区域为热带季风性气候，夏秋季雨水充沛，冬春季雨量稀少。松涛水库流域面积达1 503.27 km2，总库容33.45亿m3，是海南省第一大淡水湖泊。

图1 松涛水库地理位置Fig.1 Location of Songtao Reservoir

松涛水库是一座集灌溉、防洪、发电、养殖、旅游等多种功能于一体的大型水库，同时还是海南省最大的饮用水源保护区，为周边多个市县的生活、生产活动提供水源支撑[25]。因此，松涛水库对于海南省社会经济发展以及生态环境保护具有重大意义。本文中选取松涛水库为研究对象。

1.2 数据源

研究使用的数据为美国Landsat系列遥感数据，包括Landsat 5、7、8三代探测器获取的影像，其部分波段信息及空间分辨率见表1。基于GEE平台，选用Landsat 5、7、8的L1TP级地表反射率数据。L1TP级地表反射率数据是进行过地形精校正后的数据，它通过综合利用地面控制点和数字高程模型进行辐射和正射校正，以此来修正地形位移，是最适合像元级时间序列分析的、质量最高的1级产品[26-27]。选取1990年1月1日—2019年12月31日共30年的848景Landsat 5、7、8遥感影像(图2)。

图2 1990—2019年松涛水库Landsat卫星影像数量Fig.2 Number of Landsat images in different years of Songtao Reservoir from 1990 to 2019

表1 Landsat系列卫星遥感数据参数Table 1 Remote sensing data parameters of Landsat series satellites

为了探究气候因素对松涛水库水体面积变化的影响，研究使用了国家气象信息中心提供的中国地面降水日值0.5°×0.5°格点数据集和中国地面气温日值0.5°×0.5°格点数据集。该数据集是基于2 472个国家级气象观测站的基本气象要素资料，利用薄板样条法结合三维地理空间信息进行空间插值而得到。在剔除异常值后，计算出了1990—2019年松涛水库区域年际和月际的降水量和平均温度。

2 研究方法

2.1 水体提取

传统的NDWI是利用遥感影像的绿波段和近红外波段进行运算[11]，通过在近红外波段水体吸收强而植被反射率高的特征来提取水体信息。但是在提取城镇附近水体的过程中，NDWI计算的结果往往会因混有建筑用地信息而使得提取的水体面积比实际值偏高[12]。与NDWI相比，MNDWI很大程度上降低了建筑物对水体提取结果的影响，同时还可以区分阴影和水体，解决了水体信息中掺杂有阴影的问题[28]。MNDWI是目前应用最为有效、最为广泛的水体提取方法之一，采取MNDWI对松涛水库水体信息进行提取，其计算公式为

(1)

式(1)中：Green、MIR分别为遥感影像的绿波段和短波红外波段。Landsat 5、Landsat 7和Landsat 8传感器的绿波段和短波红外波段所在的波段编号有所差异。其中，Landsat 5、Landsat 7为波段2(绿波段)和波段5(短波红外波段)，Landsat 8为波段3(绿波段)和波段6(短波红外波段)。

2.2 水体面积计算

基于GEE平台对松涛水库1990—2019年的年际和月际水体信息进行提取并计算水体面积。水体面积计算的具体流程如图3所示。首先，选取GEE平台集成的Landsat 5、Landsat 7、Landsat 8地表反射率Tier 1数据集，分别对其进行去云处理。然后，对Landsat 8数据的波段号进行重新编号，以便其绿波段和短波红外波段的波段编号与Landsat 5、Landsat 7的波段编号一致。接着将Landsat 5、Landsat 7、Landsat 8数据集合并成一个新的数据集，筛选出新数据集中松涛水库区域1990—2019年的遥感影像，然后逐年、逐月对遥感影像的每个像素进行取均值处理，得到年际合成数据和月际合成数据。进而计算每幅合成数据的MNDWI，使用大津算法确定的阈值对水体和非水体进行分割处理。大津算法是一种使用最大类间方差自动确定阈值的方法，它根据图像的灰度特征将图像分为前景和背景两个部分，是计算机视觉和图像处理中常用的一种阈值分割方法[29-30]。当像元的MNDWI值大于分割阈值时，将像元划分为水体，反之，则将像元划分为非水体。最后，统计水体像元总个数，计算出水体面积。整个计算过程在GEE平台中通过代码实现，最后由GEE平台在线显示水体面积。

图3 水体面积计算流程图Fig.3 Flow chart of water area calculation

2.3 水体变化分析

淹没频率是指像元在一定时间范围内被观测为水体的次数占总观测次数的比例，可用来描述水体分布的空间格局及区域变化特征[31-32]。本文中采用淹没频率来分析松涛水库的水体变化特征，其计算公式为

(2)

式(2)中：F为像元的淹没频率；N为观测的总次数；Wi代表第i次观测的像元是否为水体，其中i=1表示为水体，i=0表示为非水体。

3 结果与分析

3.1 松涛水库水体时间变化

1990—2019年期间，松涛水库水体面积的年际变化波动较大(图4)。其中，1990年水体面积最大，为130.83 km2，2005年水体面积最小，为69.54 km2，分别比平均水体面积(99.19 km2)高31.90%和低29.89%。

图4 1990—2019年松涛水库水体面积变化Fig.4 Change of Songtao Reservoir area from 1990 to 2019

1990—2019年松涛水库的水体面积动态变化大致可分为5个阶段(图4)。1990—1999年的总体下降阶段，该阶段，虽然松涛水库水体面积在1996年和1997年相对于上一年有一定幅度的增长，但水体面积整体呈减少趋势，减少面积57.45 km2；1999—2002年的逐步上升阶段，2002年相对1999年水体面积扩大了40.33 km2，年平均增长面积达到13.44 km2；2002—2005年的逐步下降阶段，2005年相对2002年水体面积减少了44.18 km2，年平均减少面积14.73 km2，且在2005年松涛水库的水体面积达到了30年内的最低值。2005—2014年的总体上升阶段，期间，虽然2007年和2013年的松涛水库水体面积相对于上一年有所减少，但整体上水体面积呈增加趋势，增加面积达52.06 km2；2014—2019年的波动下降阶段，该阶段水体面积变化规律性不显著，但整体上面积减少，2019年相对2014年水体面积减少11.47 km2。

基于1990—2019年水体面积月际合成数据的计算结果，得到松涛水库水体面积月际变化曲线如图5所示。由图5可以看出，松涛水库水体面积的月际变化特征较为明显，呈先下降后上升的趋势。

图5 松涛水库水体面积月际变化Fig.5 Monthly change of Songtao Reservoir area

在整个研究时期，松涛水库的月平均水体面积最大值在1月份，为116.74 km2，最小值在7月份，为93.26 km2，平均值为105.49 km2，其动态变化过程可以分为3个阶段。1—6月份的大幅减少阶段，月平均减少面积4.57 km2；6—8月份的稳定阶段，该阶段的3个月份水体面积均处在95.00 km2以下，且彼此相差最大值仅0.65 km2；8—12月份的总体上升阶段，12月份相对8月份水体面积增加了17.86 km2，月平均增加面积达4.46 km2。

3.2 松涛水库水体空间变化

图6(a)和图6(b)分别显示了1990—2019年松涛水库年际和月际的水体淹没频率分布情况。

从空间分布来看，水库的年际淹没频率总体上呈现“内高外低”的分布格局，且由内向外淹没频率逐渐减小[图6(a)]。水库及主渠道中心位置的年际淹没频率最高，表明其在研究期间拥有永久性水体，而水库边缘的年际淹没频率较低，表明该区域在研究期间多数期次为非水体，仅有少数期次为水体。年际淹没频率较低的区域主要分布在南丰镇、马荷村、油文村、皇帝殿、总了岭、南吉、番打村、峨亮来岭、志针村以及大沙岭附近，其中番打村附近尤为显著。松涛水库月际淹没频率的分布特征[图6(b)]与年际淹没频率的分布特征[图6(a)]总体上具有较好的一致性，区别在于前者具有更大面积的永久性水体。

图6 1990—2019年松涛水库水体淹没频率分布Fig.6 Inundation frequency distribution of Songtao Reservoir from 1990 to 2019

根据淹没频率计算结果，将其划分为6个等级，分别是0

表2 1990—2019年松涛水库年际水体淹没频率统计表Table 2 Statistics of the annual inundation frequency of Songtao Reservoir from 1990 to 2019

表3 1990—2019年松涛水库月际水体淹没频率统计表Table 3 Statistics of the monthly inundation frequency of Songtao Reservoir from 1990 to 2019

年际淹没频率的统计结果表明(表2)，在近30年时间里，松涛水库拥有66.43 km2的永久性水体(F=100%)，这部分区域不存在水路更替的现象，所占比例仅有50.77%，表明松涛水库在研究期间水体面积年际变化幅度较大。其余64.41 km2(总面积的49.23%)的区域经历了水陆更替的过程，这部分区域在未来存在一定的旱涝风险，应是重点防控区域。其中，约25.67 km2(总面积的19.63%)的区域属于高淹没频率区(60%

根据月际淹没频率的统计结果可知(表3)，松涛水库拥有90.39 km2的永久性水体(F=100%)，比例达到了77.43%，占有绝对的优势。高淹没频率(60%

3.3 气候因素对水体面积变化影响分析

降水和温度可能是影响水库面积变化的气候因素。一方面，降水量增加可直接为水库提供水量，扩大水体面积；另一方面，气温的升高会促使水库表面水的蒸发，导致水面积的减小。

1990—2019年松涛水库区域的年降水量和年平均温度波动较大，但总体上均呈上升趋势[图7(a)]。其中，年降水量上升趋势不明显，而年平均温度上升趋势略微明显，研究期间共上升0.63 ℃，部分年份的年降水量与年平均气温呈负相关关系。相关性分析结果表明，年降水量、年平均温度在统计上与松涛水库水体面积无显著相关关系[图7(b)、图7(c)]。湖泊面积对气候变化响应存在一定的延迟效应[33]，因此对应年份的气候因素对当年湖泊的面积影响不显著。通过分析对应年份的降水量和平均温度与下一年度松涛水库水体面积的相关性，可看出虽然水体面积与平均温度的相关性依旧很低[图8(a)]，但与降水量表现出较强的正相关关系[图8(b)]，皮尔孙相关系数达到了0.64，表明松涛水库区域的年降水量对其下一年度的水体面积变化有着一定的影响。

图7 1990—2019年松涛水库主要气候因素年变化及与 水体面积相关性分析Fig.7 Annual variation of main climate factors and their correlation with water area in Songtao Reservoir from 1990 to 2019

图8 松涛水库主要气候因素与下一年度水体面积 相关性分析Fig.8 Correlation analysis ofmain climatic factors and water area of the next year in Songtao Reservoir

与年降水量和年平均温度变化趋势相比，松涛水库区域的月降水量和月平均温度变化特征明显，均呈先上升后下降的趋势[图9(a)]。5—10月的降水量较多，均超过了150 mm，而其中5—9月的平均温度也较高，均大于25 ℃。分别对月降水量、月平均温度与水体面积进行相关性分析，结果表明月降水量、月平均温度与水体面积之间存在着显著的负相关关系[图9(b)、图9(c)]，皮尔逊相关系数分别为-0.72和-0.88，这说明月降水量和月平均温度对松涛水库水体面积的影响显著，随着月降水量和月平均温度的增加，水体面积明显减少。

图9 1990—2019年松涛水库主要气候因素月变化及与 水体面积相关性分析Fig.9 Monthly variation of main climate factors and their correlation with water area in Songtao Reservoir from 1990 to 2019

4 结论与讨论

基于GEE云端运算平台对1990—2019年松涛水库的水体变化进行研究，分别从时间和空间角度分析了水体变化的特征，在此基础上，分析了气候因素对水体面积变化的影响，得到以下结论。

(1)1990—2019年松涛水库水体面积的年际变化波动较大，最大年份1990年与最小年份2005年相差61.29 km2，整体上呈先下降后上升的趋势。水体面积月际变化表现出明显的先下降后上升特征，从1月开始逐渐减小，7月达到最小值，具有显著的季节性变化规律。

(2)研究期间松涛水库水体表现出内部稳定、周围变化频率高的空间特征，特别在年际尺度上表现得更为明显。近30年松涛水库拥有约66.43 km2(总面积的50.77%)的永久性水体，其余64.41 km2(总面积的49.23%)的区域存在水陆更替的现象，是旱涝风险防控的重点区域。

(3)松涛水库区域月降水量、月平均温度对水体面积有着显著的影响，随着月降水量和月平均温度的增加，水体面积明显减少。在年际尺度上，区域降水量、平均温度与水体面积无明显相关性，但年降水量与下一年度的水体面积表现出较强的正相关关系。

松涛水库库区降水主要集中在5—10月，形成了7、8月份的夏汛和9、10月份的秋汛。在汛期期间，当水库蓄水达到警戒值，水库将加大泄洪量。同时，在夏秋季节由于月平均温度的增高，使得农业生产用水量增大，对水库面积变化造成影响。至春冬季节的枯水期，为了农业生产生活用水的安全，水库加大了蓄水量。人为的泄洪和蓄水是松涛水库水体面积与气候因素相关性低的主要原因。通过对松涛水库水体面积监测的研究可以为其在不同季节的泄洪量和蓄水量做出估算，为维持水库的合理水量提供参考。另外，水库的消涨区主要集中在边缘地带，这些重点区域应该合理用水、加强生态保护，为维持松涛水库的优良水质提供保障。