近50 a萨雷兹堰塞湖水域时序变化及其驱动因素研究

洪永欣,邓晚倩,张 新,

(1. 河北工程大学 地球科学与工程学院, 河北 邯郸 056038; 2. 中国科学院大学, 北京 100049;3. 中国科学院空天信息创新研究院 遥感科学国家重点实验室, 北京 100101)

0 引言

萨雷兹堰塞湖位于亚洲中部的帕米尔高原上,是世界上最深的堰塞湖,于1911年因地震致使山体大面积崩塌而形成。形成至今,萨雷兹湖水位呈波动性显著上涨趋势[1],乌索伊大坝坝体一旦崩解将会给坝体下游中亚人民带来巨大灾难[2]。近年来,萨雷兹湖所在地区引起了国际社会的高度关注[3],加强灾害的识别与评估,可有效提高灾害风险防治的水平[4]。堰塞湖作为洪水灾害的潜在因素,已有多学科学者对其形态[5]、稳定性[6]、崩溃预测[7]做相关研究;对于萨雷兹堰塞湖研究的主要集中在坝体监测预警[8-10]、水量监测[1]、风险评估[9]、防灾减灾[3,11]、地震对萨雷兹湖的影响[12-15],洪水演算模拟[2]等方面。因此,监测萨雷兹湖水资源动态变化,研究湖泊水位变化特征规律对于萨雷兹湖的综合治理、水量调控等具有重要的指导作用以及科学价值。

国内外学者对萨雷兹湖水资源多从成因灾害评估、溃坝模拟、影响等角度对其进行了阐述,而对萨雷兹湖水资源水位水量积蓄过程中的变化特点和发展趋势、驱动因子分析等相关的内容研究还比较缺乏。然而了解萨雷兹堰塞湖积蓄过程,及影响该过程的相关因素,可以达到直接和间接的观测效果,实现预测萨雷兹湖水位水量趋势、溃坝风险,最大程度减少或者避免损失。归一化水体指数NDWI(normalized water index, NDWI)是遥感水体最为有效的提取方法之一[16]。MCFEETERS等[17]发现利用NDWI方法在近红外波段水体吸收最强,而植被和土壤在近红外波段反射率较高,李飞等[18]在研究复杂水体信息提取方法研究中表示,NDWI能够很大程度上去除植被信息,强调水体。已有的图像分割方法有:基于区域[19]、边缘[20]、学习分类[21]和阈值[22]的方法等,其中采用阈值分割对研究区内水体遥感影像信息进行提取的方法是最为广泛,吕争等[23]提取水体时,在初提取阶段灵活运用全局阈值分割法中的最大类间方差算法,在精处理时采用局部阈值分割法,全局-部分相结合精度高,工程量小。文中遂以萨雷兹堰塞湖为研究对象,从区域、季节等多个角度对萨雷兹堰塞湖水域水量变化进行监测,并结合遥感数据以及气温、降水、蒸散发等气象因素,对萨雷兹湖水资源变化的驱动因素进行探讨。

1 研究区域及数据来源

1.1 研究区概况

萨雷兹湖位于塔吉克斯坦的帕米尔高原,位于72°32′30″～73°11′21″E,38°9′21″～38°20′20″N,萨雷兹湖湖区海拔高度约为3 263 m,湖集水面积16 506 km2,萨雷兹湖年平均长度约 55.8 km,平均宽度为 1.44 km;流域整体的平均海拔为3 263 m,堰塞湖总体积约17 km3,最深点在500 m,平均深度为202 m。巴尔坦格河为萨雷兹湖的下游,穆尔加布河为上游。萨雷兹湖年平均环境气温为1 ℃,近年来年平均降水量为108 mm。

1.2 数据资料

1.2.1 Landsat数据

Landsat是美国航天局和美国地质调查局共同管理维护的对地卫星,主要用于探究地下矿物资源、海洋能源和地下水资源,研究农作物生长规律、预测农作物产量,调查预测自然灾害,制作专题地图等。考虑时间节点、云层干扰、湖面结冰情况,最终选取了1972—2019年间的Landsat分辨率多光谱数据,用于分析萨雷兹湖面积变化长时间序列。

1.2.2 气温与降水数据

气温数据来源于全球历史气候学网络数据库(www.ncdc.noaa.gov/cdo-web)。文中选择了距离萨雷兹湖最近、数据覆盖最完整的霍罗格气象台站来研究气温数据。降水数据来自中亚生态与环境研究中心,采用了萨雷兹湖水文监测站点的1972—2018年萨雷兹湖日观测水位数据以及2005—2019年萨雷兹湖入湖流量日观测数据,用以分析萨雷兹湖水资源变化及其变化的驱动因素。

1.2.3 潜在蒸散

潜在蒸散数据来自东英吉利大学气候研究格网单元与英国国家大气研究中心联合发布的1901—2019年期间全球陆面月平均数据集(CRU-TS-4.04)。其空间分辨率为0.5°×0.5°(https://catalogue.ceda.ac.uk/)。文中采用的是萨雷兹湖湖面1972—2019年间潜在蒸散量的月平均值。

2 研究方法

综合考虑影像的时间节点、云层覆盖以及湖面结冰状况,以及萨雷兹湖周围山体众多,存在影像被干扰的情况,选用了归一化水体指数与全局-局部阈值分割技术相结合的自适应水体提取算法[24]对萨雷兹湖的湖岸廓线进行提取。

2.1 归一化水体指数提取

本研究采用的是归一化水体指数NDWI,可以实现有效地增强水体光谱特征并抑制植被等环境信息。MCFEETERS等[17]给出的NDWI公式如式(1)所示:

(1)

式中:ρGreen为遥感影像波段中大气反射率的绿波段,对应Landsat的MSS中的波段1;ρNIR为遥感影像中表观反射率的近红外波段,对应于Landsat中的波段3,在TM/ETM+中对应的波段为波段2与波段4,因此可以灵活地获取长时间序列的提取。由于不同传感器的标定系数和太阳光谱辐照度不同,NDWI的计算是基于大气表观反射率。

2.2 全局-局部分割阈值选择

通过对NDWI的物理特征的了解,从中选取全局阈值,其阈值的范围通常为[-0.1,0.1]。第1步,全局分割整幅图像。首先通过设定初始全局阈值(如T0=-0.1)对影像进行分割,对分割结果向外扩充若干单元设立相应的缓冲区,其中缓冲区的像元总数与湖泊像元总数相同,得到初分割湖泊。第2步,针对萨雷兹湖的特征作局部分割,其周围的背景像素作为局部分析单元。将上述包含缓冲区的局部单元进行直方图统计,若频率分布直方图呈双峰分布,则通过式(2)阈值将湖泊与缓冲区进一步分割,得到新湖岸廓线;对该湖泊矢量建立新的缓冲区,按照下述阈值进行再次分割迭代,至2次湖泊矢量的像元个数差值小于给定阈值。

(2)

式中:μwater为水体的平均值;μland为陆地像元的平均值;σwater、σland分别为水体与陆地像元的标准差。

第3步,混合像素。由于山体阴影与湖泊水体具有相似的光谱特征,因此仅通过影像中的光谱信息难以对其进行有效的区分。通过DEM数据生成的坡度图进行分割,将坡度图与湖泊矢量进行叠加,剔除山体阴影的像素点。最后,为保证提取矢量图层的高精度,针对每期提取后的结果进行目视判读,迭代检验矫正,最终可以得出萨雷兹湖矢量数据。

2.3 精度检验

对于遥感数据预处理后所得萨雷兹湖岸线精度,本研究聚焦于萨雷兹湖的动态变化,因此采用目视判读的方法。凭借实践经验和专业知识在萨雷兹湖DEM遥感影像上获取地物特征,实践表明判读的结果满足萨雷兹湖水体研究动态变化的条件。在分析水文参数趋势变化检验方面,使用Mann-Kendall非参数统计法[25-26](M-K)来研究萨雷兹湖动态变化趋势。文中统计了水域各季节的平均面积变化,运用M-K以四季为单位,统计构建得到萨雷兹堰塞湖水体面积的变化趋势与突变时间拐点,具体计算方法如式(3),针对水体面积时间序列:x1,x2,x3,…,xn,构造一个秩序列ri,ri表示xi>xj(1≤j≤i)的样本累积数。定义Sk:

(3)

式中,Sk的均值E(Sk)及方差var(Sk)定义如式(4):

(4)

水域面积时间序列是随机独立分布为假设的前提条件,那么定义统计量UF表示为:

(5)

式中,UFk为标准正态分布,其显著性水平为0.05时的临界值UFα=±1.96,当|UFk|>UFα时,测得的萨雷兹湖水域面积序列表现出显著的增长趋势或减少的趋势;其中,把UFk绘构为曲线c1。

以上述的计算方法作用于反序列上,迭代计算,并将结果乘以-1,得到UBk;所有的UBk在图中构成曲线c2。通过绘制c1、c2曲线图,若UFk或UBk的值大于0,则表明萨雷兹湖水域面积时间序列中的该时段上为上升趋势,反之则为下降;当曲线超过临界值时,表明该时间段内存在剧烈的上升或下降幅度;若c1、c2两曲线交点在临界线区间内,表示在该交点处,为突变点的起始时间可能性大。

3 结果与分析

3.1 萨雷兹湖水域面积时序变化

经过对Landsat遥感影像数据进行NWDI和全局-局部分割阈值处理后的萨雷兹湖水域面积时序变化进行具体分析,如图1所示。在1972—1998年间,水域面积呈缓慢上升趋势,上升速率达到0.189 km2/a;1999—2004年间,水域面积呈缓慢下降趋势,下降速率为-0.293 km2/a;2005年,萨雷兹湖水域面积急速扩张,上升速率高达到4.038 km2/a; 2006—2009年水域面积呈现出急速下降的趋势,下降速率为-0.410 km2/a;2010—2019年水域面积始终保持周期性上升,上升速率为0.221 km2/a;将萨雷兹湖的面积变化时间序列按照逐年季节进行平均计算,春夏秋季节萨雷兹湖水域面积整体变化呈振荡上升趋势,年均的变化速率略有不同,春季水域面积的年平均增长率为0.157 km2/a,为一年中增长最快的季节;秋季在1972—1988年之间变化的水域面积增加幅度高达4.77 km2;四季的水域面积变化率,春夏秋按顺序逐渐减少,冬季保持波动性平稳。

图1 萨雷兹湖时间-面积变化序列图Fig. 1 Time-area change sequences of Sarez Lake

将萨雷兹湖水域面积时序进行年平均结合Mann-Kendall非参数检验分析,如图2。由于数据缺失的原因,暂不考虑1976年,近40 a 的UF统计量的值均为正值,发现在1991—1993年这3 a之间的UF统计量在临界值之上,其显著水平均大于0.05临界值,则通过具体分析得到萨雷兹湖水表面积具有显著的上升趋势;具有同样突出特征的年份有:1996—2001年以及2013年至今。然而,在2012—2013年,UF与UB两统计量曲线有交点,则表明在该年份间萨雷兹湖湖面的年均面积变化剧烈。

图2 萨雷兹湖年平均水域面积突变检验分析Fig. 2 Mutation test of annual average water areas of Sarez Lake

在春夏秋冬四季中水域面积的变化有所不同,经M-K突变检测得春夏两季表现规律相似,UF统计量的值始终为正值。2006年以后春季平均水域面积上升趋势显著,2012年以后水域面积具有显著上升特征的季节为夏季,其UF统计量表现为在0.05临界值线之上,此外春夏两季的临界线之间的交点分别位于2005年与2011年;秋季在1972—1983年UF统计量为负值,但UB统计量仍为正值,水体表面积始终处于上升趋势,其余年份UF统计值均为正值,自2015年以来,其UF统计值均超过0.05显著水平的临界点,说明水域面积有明显上升趋势,水域面积变化的突变点与夏季相近,位于2011—2012年间;因此,1997—2000年及2003—2013年水域面积呈上升趋势,且上升趋势较平缓,在1972—2019年间上升下降趋势变化幅度较大。因此,萨雷兹湖水域面积变化在2010—2013年间存在突变拐点,且变化的趋势成剧烈状态。

3.2 萨雷兹湖水域体量变化

本研究获取了萨雷兹湖的实测站点水位数据,获取了1972—1991年、2005—2011年、2014—2018年3个时段水位高程,分析水位在该3段时间序列里的变化,得出1971—2018年水位总体的趋势为波动性上升。继而,根据现有的观测数据,主要包括5个时间段的结果:1972—1980年、1981—1985年、1986—1991年、2005—2011年、2014—2018年,表现出来的变化规律略有区别,对这5个阶段的水位时间序列分别建立时间-水位变化函数,整体上1972—2018年间水位呈0.147 m/a的速率周期性上升,在1972—1980年间水位上升速率为0.560 m/a,而1981—1985年水位快速上涨速率达到0.784 m/a,在1986—1991年间水位呈0.356 m/a的波动型缓慢上涨;2005—2011年为唯一时间段水位呈周期性下降期速率达到-0.473 m/a,随后在2014—2018年,水位再度呈现上涨趋势速率达到0.444 m/a。为研究萨雷兹湖水域面积与同期水位之间的相关性,建立了两者的关系曲线,通过时间匹配选取了133组面积水位数据点,二者显示出了良好的相关性,并建立面积-水位关系模型:

S=-0.036 55×dh2+4.762 07×dh-68.717 94

(6)

式中:dh为相对水位(m),dh=h-3 200;S为水域面积(km2)。

根据湖泊水位时间序列与面积时间序列的计算,可以按梯形台锥台体积式(7)估算湖泊的体积变化,从而获得萨雷兹湖体积变化:

(7)

式中:H1为最低水位高度;H2为最高水位高度;A1为高水位时湖泊的面积;A2为低水位时湖泊的面积。

由于水域面积与水位序列的时间匹配不佳的问题,文中通过面积-水位关系模型获取与水位数据相匹配的面积序列,并结合水位时间序列来进行湖泊体积变化的计算。1972—2018年萨雷兹湖量总体表现为波动增长的趋势,其中主要包括6个不同时间段的结果:1972—1980年、1980—1983年、1983—1985年、1985—1991年、2005—2011年、2014—2018年,对这6个阶段的水量时间序列分别建立时间-水量变化函数,整体上1972—2018年间水量呈0.012 km3/a的速率周期性上升。

3.3 萨雷兹湖水域变化驱动力因素

湖泊水量变化的驱动因素众多,比如有降水量、蒸散发、冰川补给等。为了深入了解萨雷兹堰塞湖水资源变化的驱动因素,文中结合气象数据、估测的萨雷兹湖流域冰川变化以及萨雷兹湖体积变化,并进行了驱动力分析。

3.3.1 气象因素

对萨雷兹湖水量变化的驱动力进行气象因素方面研究,分析降水量的日数据变化发现,最大的降水量出现在春季(3—5月),但从7—9月几乎没有降雨。比较萨雷兹湖水量时序数据与降水量数据,发现两者的月变化趋势没有较大关联性,通过时间序列的匹配与筛选,将二者时间序列对齐进行相关系数的计算后结果为-0.265,这与常规的湖泊水量与总降水量呈正相关不符。

将萨雷兹湖月平均水量与气温变化比较,见图3,发现两者并无直接关系。采用全球陆面月平均地面数据集的潜在蒸散量对水量变化进行了进一步的相关分析,如图4所示,1972—2019年间萨雷兹湖所在区域潜在蒸散量总体上呈相对平稳的趋势,但水量与潜在蒸散相关系数为-0.103,整体上两者之间表现出负相关性,并且相关性较弱。

图3 萨雷兹湖月平均水量变化与月平均气温变化比较 图4 萨雷兹湖月平均水量变化与潜在蒸散量变化比较Fig. 3 Comparison of the temperature data near Sarez Lake and the water volume change Fig. 4 Comparison of the potential evapotranspiration data near Sarez Lake and the water volume change

虽然月平均气温与萨雷兹湖月平均水量变化并无直接关系,但是萨雷兹湖水量的年际变化趋势与年平均气温、年均最大气温(年内每月最大气温的平均值)、 年均最小气温(年内每月最小气温的平均值)的变化相比较,发现其趋势一致(如图5所示)。多年来,萨雷兹湖气温年际变化呈季节性波动, 并缓慢上升, 对齐二者时间序列相关分析结果分别为0.579、0.601、0.439,表明萨雷兹湖年际水量变化与区域年均温度具有正相关性。推测其随着气温的升高影响萨雷兹湖的水文特征变化状况,使萨雷兹湖水资源得到了补给。

图5 萨雷兹湖年平均水量变化与多年年均气温、多年年均最高气温、多年年均最低气温变化比较Fig. 5 Annual temperature data, the maximum temperature data, and the minimum temperature data near Sarez Lake and comparison with annual water volume change

3.3.2 入湖补给

萨雷兹湖水量近年来始终保持周期性上涨,为更好地显示水资源变化的影响因素,本研究给出了萨雷兹湖水位变化与穆尔加布河入湖流量的变化关系,研究发现水位变化与湖泊入湖流量变化的趋势在大部分时间内表现一致,两者的相关系数为0.538,相关性良好,因此可以认为萨雷兹湖的水位上涨很大程度上源于穆尔加布河对其的补给。由于萨雷兹湖的水位变化呈现明显的季节性周期变化,这表明穆尔加布河的补给作用可能受到其他因素的影响。研究表明,在1972—2000年间,由于气候变暖导致帕米尔高原东部的冰川融溶体量增大[27],流域内的冰川冻土消融,融水一部分可能通过穆尔加布河补给汇入萨雷兹湖,导致穆尔加布河入流量的上升,最终给萨雷兹湖的水位变化特征造成影响。

3.3.3 流域冰川消融

分析流域冰川的变化对萨雷兹湖水量变化的影响,对入湖流量与流域冰川变化进行了相关分析。由于穆尔加布河上的产流区主要包括融雪产流与冰川产流,其中融雪产流主要集中在3—7月,冰川产流集中在6—9月,如图6(a)为入湖流量的月平均变化,可以看出萨雷兹湖入湖流量较大的月份在7—9月,属于穆尔加布河的冰川产流阶段之内。此外,2003—2009年间,萨雷兹湖入湖流量呈季节性减少趋势;在2010—2018年间,入湖流量呈增加趋势。图6(b)所示,显示了萨雷兹湖流域冰川面积与年平均入湖流量变化的关系,可以看出冰川面积变化与萨雷兹湖流量变化趋势呈负的相关性,年平均入湖流量与北部冰川面积、南部冰川面积的相关系数分别为-0.531、-0.462。图6(c)中冰川体积变化的趋势与流域内的河流流量的变化趋势变现出显著负相关,入湖流量与流域总体冰川体积变化,北部山脉冰川体积变化,南部山脉冰川体积变化之间的相关系数分别为-0.748、-0.781、-0.612,可以看出,无论是南北部山脉冰川变化还是整体流域山脉的冰川与流域内的水资源补给都有着密切的联系。图6(d)显示了流域冰川体积变化对萨雷兹湖水量变化造成的影响,萨雷兹湖在2003—2009年间年平均水量变化呈减少趋势,在2010—2018年呈缓慢上升趋势,与流域冰川体积在2003—2009年增加,2010—2019年间减少的趋势相反。

图6 萨雷兹湖流域冰川变化与入湖流量变化比较Fig. 6 Comparison of the inflow discharges of Sarez Lake and glacier changes

可以看出,萨雷兹湖的水资源变化与流域冰川积雪消融密切相关,2003—2009年间受多种气候因素影响流域冰川的消融过程减弱,此时,积雪消融的径流为主要入湖补给,入湖流量与萨雷兹湖坝体渗透量相比可能相对较小,因此,总体上2003—2009年间萨雷兹湖水量呈下降趋势;2010—2018年间流域冰川受到同期温度的影响,冰川消融呈缓慢加剧趋势,冰川融水导致穆尔加布河补给量增加,最终影响萨雷兹湖水位上升。假设在冰川系统和传输损耗(蒸发、渗透损失)过程中的蒸发和升华过程引起的质量损失可忽略不计,则根据入湖流量年平均变化与冰川体积年平均变化,计算2010—2018年间入湖流量在损耗的冰川体积中的年均占比可知:2010—2018年间萨雷兹湖流域冰川质量损失中有超过15.4%的比例被汇入穆尔加布河的河道产流,最终汇入萨雷兹湖,导致湖泊水量增多。

4 结论与讨论

文中针对萨雷兹湖流域受地形因素和天气要素的影响,传统方法难以高效观测的问题,采用光学影像、气象数据等对萨雷兹湖水资源变化进行监测,并探讨湖水变化的相关性因素,得到的主要结论如下:

1)结合长时间序列Landsat影像与萨雷兹湖实测水位数据,监测了萨雷兹湖1972—2019年间的湖岸廓线并计算了水域面积时间序列,结果表明萨雷兹湖水域面积变化在2010—2013年间存在突变拐点,变化趋势加剧,萨雷兹湖呈现出一定的季节性周期变化规律,并且四季水域面积变化各有所区别,1972—2019年水域面积整体上呈现了逐渐增加的趋势。

2)结合水文站点实测水位数据建模,同时进行了湖泊体积变化的计算,构建了萨雷兹湖1972—2018年的水量变化时间序列。表明1972—2018年萨雷兹湖水量整体以0.012 km3/a的变化速率波动上升。依据变化趋势不同,分为了6个不同时间段的结果水量并呈现出“缓慢上升—缓慢下降—快速上升—缓慢上升—缓慢下降—缓慢上升”的变化规律。

3)结合气温、降水、潜在蒸散等气象因素以及萨雷兹湖入湖流量、流域冰川变化,分析了萨雷兹湖水资源变化的驱动因素。分析结果表明萨雷兹湖水资源变化与降水补给、蒸散发在统计意义上均无显著的相关关系,但与年平均温度变化呈显著正相关,还与穆尔加布河径流补给量密切相关。综合萨雷兹湖流域冰川的变化,表明,萨雷兹湖流域冰川的消融与其水文变化特征具有直接联系,根据入湖流量年均变化与冰川体积年均变化可知2010—2019年间,萨雷兹湖流域冰川质量损失的15.4%多在穆尔加布河的河道产流中。

下一步研究需要考虑短时间内的降雪与融溶易对实验造成影响以及对驱动因素研究需要增加大量验证数据等问题,进而为全球变暖的大背景下萨雷兹堰塞湖溃坝的风险评估提供科学支撑。