文章快速检索     高级检索
  大地测量与地球动力学  2021, Vol. 41 Issue (12): 1282-1287  DOI: 10.14075/j.jgg.2021.12.015

引用本文  

刘慧玲, 陈雨, 王希禾. 基于球谐变换的GERD水库形变预测[J]. 大地测量与地球动力学, 2021, 41(12): 1282-1287.
LIU Huiling, CHEN Yu, WANG Xihe. Deformation Prediction of GERD Reservoir Based on Spherical Harmonic Transform[J]. Journal of Geodesy and Geodynamics, 2021, 41(12): 1282-1287.

项目来源

四川大学社会科学研究项目(skbm201915)。

Foundation support

Social Sciences Research Project of Sichuan University, No.skbm201915.

通讯作者

陈雨,博士,副教授,主要研究方向为GRACE重力卫星数据应用、遥感图像处理,E-mail: ychen@scu.edu.cn

Corresponding author

CHEN Yu, PhD, associate professor, majors in GRACE gravity satellite data application, remote sensing image processing, E-mail: ychen@scu.edu.cn.

第一作者简介

刘慧玲,硕士生,主要研究方向为信号与信息处理,E-mail: liuhuil@sina.cn

About the first author

LIU Huiling, postgraduate, majors in signal and information processing, E-mail: liuhuil@sina.cn.

文章历史

收稿日期:2021-02-23
基于球谐变换的GERD水库形变预测
刘慧玲1     陈雨1     王希禾1     
1. 四川大学电子信息学院,成都市一环路南一段24号,610065
摘要:提出一种基于球谐变换预测GERD水库蓄水过程中引起的地表形变的新方法。首先假定GERD水库分别在5 a、10 a、15 a间被注满,利用SRTM-DEM数据求取蓄水过程中每月水库的面积、体积和平均水高。经球谐展开后计算对应蓄水过程中引起的地表形变,并通过提高最大谐波度的上限来增加结果的精确度。3种蓄水策略下引起的垂向位移速率分别约为11.8 mm/a、5.9 mm/a、3.9 mm/a;引发的东向和西向位移速率分别约为1.8 mm/a和1.4 mm/a、0.9 mm/a和0.7 mm/a、0.6 mm/a和0.4 mm/a;引发的南向和北向位移速率分别约为2.6 mm/a和4.4 mm/a、1.3 mm/a和2.2 mm/a、0.8 mm/a和1.5 mm/a。
关键词球谐变换GERD垂直形变水平形变时间序列

水坝通过调节水流,为人类和环境的各种需求提供了可靠水源[1]。埃塞俄比亚的文艺复兴大坝(Grand Ethiopian Renaissance Dam,GERD)是现今非洲最大的大坝,预计其有效库容可达74 km3[2]。GERD水库蓄水将引发大量的水文迁移,可能导致周边地区发生大规模的地表形变[3],地表形变量超出一定界限会演变成为地质灾害,对人类生命和财产安全带来严重损失[4]

地表形变监测网络是监测区域地表形变的有效方法之一[5],但该方法不适用于GERD水库区域,因为GERD附近240 km范围内没有可用的GPS站点,距离GERD 142 km的ASOS站点已于2017年停止工作。卫星遥感可以进行长时间、高覆盖率的观测,能连续监测水库蓄水和地表覆盖变化,许多研究者利用多种遥感数据相结合成功估算了特定地区的地表形变[6-7]。GRACE以及最新的GRACE-FO三大数据处理机构(JPL、CSR、GFZ)将地表引力场变化通过球谐变换,得到随时间变化的GRACE重力数据的球谐系数(spherical harmonic coefficients,简称SHC)形式(即Level 2产品)[8],利用该产品并引入负荷勒夫数(Love numbers)可反演得到相应重力变化引起的地表负荷形变[9]。但由于GRACE的分辨率只有200~300 km,其时变重力场的SHC的阶数只能展开到60~96阶[10],受制于分辨率的影响,大部分研究中基于GRACE反演得到的地下水储量变化与实际结果的相关系数仅为0.6~0.7[11]。所以,通常利用GPS数据或其他水文模型对GRACE反演后的结果进行校正。

根据以上分析,本文基于球谐变换和Love numbers提出一种预测GERD水库蓄水引起的地表形变的方法。首先假设3种蓄水方案:短期5 a、中期10 a、长期15 a内注满74 km3的水库,并根据数字高程模型(DEM)分别对3种蓄水方案中每月蓄水表面积、体积和平均水高进行数字拟合。然后以平均水高为基础创建一个7 200×3 600的栅格网络,该栅格水库区域的各点设为平均水高值,水库区域外各点设为0。通过对比等效水高(equivalent water height,EWH)和平均水高分析高阶球谐截断所造成的信号能量损失,再以2 300阶的球谐变换系数为基础,通过引入Love numbers计算水库及周边区域的垂直和水平形变。

1 数据准备和研究方法 1.1 研究区域及数据准备

GERD位于埃塞俄比亚的蓝色尼罗河上,地处苏丹和埃及的上游,于2011-04开始建设,计划正常水高为640 m,水库蓄水体积为74 km3,约为大坝所在地年均流量的1.5倍[12]

采用美国宇航局的SRTM-DEM数据(https://dwtkns.com/srtm30m/),空间分辨率为30 m。根据GERD的实际位置先选取32°50′~48°9′E、3°21′~ 15°1′N的区域,然后将GERD的位置作为倾泻点,从上述区域中划分出本文的研究区域——上蓝色尼罗河流域。

1.2 研究方法

本文的数据处理流程分为5个步骤(图 1):1)利用SRTM-DEM数据对研究区域在不同蓄水策略下提取的水体进行分析,得出蓄水过程中每月水库的蓄水表面积、体积和平均水高;2)构建全局栅格图,将水库区域内的栅格值设为当月平均水高,其余区域设为0;3)将此栅格网络球谐扩展为最大谐波度(Nmax)分别为60、700、1 500、2 300的SHC;4)利用SHC计算EWH,先通过比较EWH和平均水高,分析高阶SHC对信号能量损失的影响,再利用Nmax=2 300阶时的SHC计算地表形变;5)提取垂直、东、西、南和北向形变的时间序列图。

图 1 数据处理过程 Fig. 1 Data processing
1.2.1 计算水库蓄水表面积、体积和平均水高

首先使用一组高程值间隔1 m的水平表面截取研究区域的SRTM-DEM栅格图,获得指定存储层中每个水平表面下的体积,再构建以体积为自变量的体积/水平面高程值拟合函数。如果以74 km3作为水库额定容量,将蓄水年限分为5 a、10 a和15 a分别进行计算,对于3种年限,每月分别需蓄水1.23 km3、0.61 km3和0.41 km3。利用体积/水平面高程值拟合函数计算不同蓄水策略下每月蓄水体积对应的水平表面高程值;然后使用ArcGIS工具计算水库存储层中每月水平表面高程值下的蓄水表面面积;最后用月蓄水体积除以蓄水表面积得到每月平均水高。

1.2.2 设计平均水高栅格及球谐展开

本文构建了一个分辨率为0.05°×0.05°的包含全球经纬度的栅格网络。水库范围内的栅格值设为当月平均水高,其余栅格值设为0;然后将此全局栅格网络球谐展开至指定Nmax的SHC,并根据Wahr等[13]的方法求出每月变化量。

1.2.3 计算EWH和地表形变

利用球谐函数计算EWH和地表变形时需引入Love numbers,每月SHC变化量结合Love numbers可计算由水文质量载荷变化引起的EWH变化和3D位移。ΔEWH的具体计算公式可参考文献[14],3D位移的具体计算公式可参考文献[15-16]。

2 结果分析 2.1 每月水库蓄水表面积、体积和平均水高

图 2是5 a、10 a和15 a蓄水策略的模拟实验中水库蓄水表面积和平均水高随时间变化曲线。由图可见,3种蓄水策略下GERD水库平均水高的月平均变化量分别为1.63 m、0.88 m和0.61 m,水库蓄水表面积的月平均变化量分别为28.66 km2、14.59 km2和9.72 km2

图 2 3种策略下每月平均水高与水库表面积 Fig. 2 Monthly average water height and area under three strategies

图 3为3种蓄水策略在不同蓄水阶段的体积、蓄水表面积和水位的变化。分析可知,5 a策略的第1个月蓄水体积为1.23 km3,其水面面积为82.52 km2,平均水高为14.96 m;10 a和15 a蓄水策略中对应的数据为0.61 km3、46.05 km2、13.32 m和0.41 km3、34.11 km2、12.06 m。当5 a蓄水策略完成时,10 a蓄水策略完成50%,15 a蓄水策略完成约33%;当10 a蓄水策略完成时,15 a蓄水策略完成约66%;当整个蓄水完成时,会形成一个体积约为73.99 km3、面积约为1 773.2 km2、平均水高约为41.72 m的水库。

图 3 3种策略模拟水库容量、水位、面积 Fig. 3 Simulated reservoir volume, area and water height using three strategies
2.2 基于不同Nmax的EWH对比与分析

图 4(a)为5 a蓄水策略中第30个月的平均水高图,水库内的栅格值为34 m(蓝色区域),水库外栅格值为0(黄色区域)。在此蓄水阶段中,根据不同Nmax计算各自的EWH,结果如图 4(b)~4(e)所示。

图 4 5 a蓄水策略下第30个月的平均水高与EWH Fig. 4 Average water height and EWH of the 30th month under 5 a filling strategy

Parseval能量守恒定理表明,信号在空间域的平均功率等于时域中各次谐波分量的平均功率之和[17]。虽然Parseval定理通常用于描述傅里叶变换的统一性,但该特性的最一般形式应更恰当地被称为Plancheral定理。在数学中,Plancheral定理是谐波分析的结果,所以能量定理也适用于谐波分析。当研究区域为面积较大的水库或湖时,区域地表质量变化的主要影响因素是水库或湖中水的储量变化。研究中常用ΔEWH(θ, φ)来表示陆地水储量的变化[18]

由于球谐扩展受Nmax的限制,高阶谐波被截断导致信号在空间域被平滑和扩展,所以利用球谐函数计算ΔEWH和3D位移是一种数值再分配模式下的计算,计算后的区域比研究的水库区域要大得多。图 4(b)显示了由于高于60阶的谐波信号被截断导致信号能量被平滑、衰减和扩展至水库边界之外的现象,其栅格单元中EWH的范围为-0.164~0.123 m。

随着Nmax的逐渐增大,受截断效应影响的信号逐渐减少,EWH信号向GERD水库中心集中,呈现出水库中心区域的EWH值明显大于周边区域的特性。当Nmax=700,水库中心区域的EWH最大值为16.5 m,EWH的范围为-3.2~16.5 m;当Nmax=1 500时,水库中心区域的EWH最大值增加至30.03 m,EWH的范围增加至-3.83~30.03 m;当Nmax=2 300时,已显示出了EWH沿GERD水库区域分布的特性,EWH的范围为-8.28~38.3 m(由于吉布斯现象会出现少数大于34 m和小于0的栅格值)。

利用Nmax=2 300时的SHC计算5 a蓄水策略下的平均水高与EWH的时间序列。经误差分析,两者的平均相对误差仅有1.6%,因此选定Nmax=2 300阶计算后续地表形变。

2.3 地表形变及其时间序列 2.3.1 垂直形变

计算5 a策略中第30个月的垂直形变,结果如图 5(a)所示。从图中可以看出,栅格中所有质量载荷点的垂直位移分量都向水库质心运动,为负值。该月垂直形变位移的范围为-91.13~-0.49 mm, 且越靠近水库区域垂直形变位移量越大,水库中心区域的垂直形变位移范围为-91.13~-55.04 mm。

图 5 5 a蓄水策略下第30个月的形变结果 Fig. 5 Deformation results of the 30th month under 5 a filling strategy
2.3.2 水平形变

水库区域质量载荷点的位置决定了该点的位移方位(向北或向南),由于水库北部面积大于南部,因此质量载荷在南部单元上产生的拉力大于北部。在5 a策略下第30个月的蓄水阶段中,南向位移的最大值约为9 mm,北向位移的最大值约为37 mm。另外,由于水库北部区域较宽,因此东部和西部的单元也包含了向南移的分量,如图 5(b)所示。

东西向位移与南北向位移相似,水库西岸的载荷点运动趋势向东,表现为正值,水库东岸的载荷点则向西移动,表现为负值,如图 5(c)所示。由于GERD水库东西方向分布较为均匀,所以呈现的东西向位移量也较为均衡,东向位移的最大值为9.31 mm,西向位移的最大值为10.59 mm。

2.3.3 地表形变的时间序列

依据图 5的结果分析垂直向、南北向和东西向的形变区域,并计算3种蓄水策略下5个方向的时间序列,结果如图 6所示,图中,V表示垂直形变,E表示东向形变,W表示西向形变,N表示北向形变,S表示南向形变。

图 6 3种蓄水策略下水平形变的时间序列 Fig. 6 Time series of horizontal deformation under three strategies

使用恒定体积的水蓄满水库时,蓄水初期和中后期将产生较大的垂直位移, 5 a策略下垂直位移变化的最快速率和最慢速率分别为1.2 mm/月和0.4 mm/月;10 a策略下对应速率降为1.0 mm/月和0.35 mm/月;15 a策略下仅为0.45 mm/月和0.2 mm/月。这3种策略的年平均垂直位移量呈递减的趋势,分别约为11.8 mm、5.9 mm、3.9 mm。

图 5(c)所示,5 a蓄水策略下第30个月的东西向位移网格图显示出一定的对称性;这种对称性在东西向位移的时间序列中也有体现,如图 6(b)所示。本文认为,这种对称性是由于水库东西狭窄的形状所引起的。5 a策略的东向年平均位移量约为1.8 mm,最快的变化发生在蓄水的第9个月,达到0.42 mm/月,在蓄水的最后时期,移动速率降到0.1 mm/月;10 a策略下最快和最慢的东向位移速率为0.35 mm/月和0.07 mm/月,分别发生在蓄水期的第19个月和第80个月前后;而在15 a策略下,东向形变最大速率为0.25 mm/月,最小速率为0.03 mm/月。

GERD水库南部和北部质量分布不平衡导致南部和北部的位移分量趋势不同,北部较大的区域可能会分散质量负荷,北部质量载荷单元的移动小于南部。南部和北部网格单元的位移时间序列如图 6(c)所示,负值表示北部的网格单元向南移动,正值表示南部的网格单元向北移动,两者均向水库的质心移动。对于5 a、10 a和15 a蓄水策略,北向位移的年平均位移量分别为4.4 mm、2.2 mm和1.5 mm,南向位移的年平均位移量分别为2.6 mm、1.3 mm和0.8 mm。

表 1为不同蓄水策略下的年位移变化、位移总变化和最大位移变化。可以看出,蓄水的速度越快,年变形量越大。总体而言,在水库蓄水完成后(注满74 km3),存储层将产生约59 mm的垂直位移、22.5 mm的北向位移、13 mm的南向位移、9.2 mm的东向位移和6.7 mm的西向位移。

表 1 不同蓄水策略下的年位移变化、总位移变化和最大位移变化 Tab. 1 Statistics of annual displacements, total displacements and maximum displacements under different strategies
3 结语

本文提出了一种基于球谐变换预测GERD水库地表位移分量的方法。分别假定5 a、10 a、15 a蓄满水库,则由此引发的地表垂直形变速率分别约为11.8 mm/a、5.9 mm/a、3.9 mm/a;引发的南向和北向位移速率分别约为2.6 mm/a和4.4 mm/a、1.3 mm/a和2.2 mm/a、0.8 mm/a和1.5 mm/a;引发的东向和西向位移速率分别约为1.8 mm/a和1.4 mm/a、0.9 mm/a和0.7 mm/a、0.6 mm/a和0.4 mm/a。同时实验发现,蓄水速度越快,月/年形变量越大。

其次,通过60阶、700阶、1 500阶和2 300阶的球谐展开对比发现,高阶次球谐波的截断会导致边缘信息丢失和空间信息的平滑和扩展。为了减少该效应对本文结果的影响,经过对比,最后选定2 300阶的截断来完成实验。对能量守恒定理的分析和输入的平均水高与所求EWH仅1.6%的误差率均证明了本实验的合理性。

参考文献
[1]
Ehsani N, Vörösmarty C J, Fekete B M, et al. Reservoir Operations under Climate Change: Storage Capacity Options to Mitigate Risk[J]. Journal of Hydrology, 2017, 555: 435-446 DOI:10.1016/j.jhydrol.2017.09.008 (0)
[2]
Asim I el M. On the Grand Ethiopian Renaissance Dam (GERD)[J]. International Journal of Environmental Sciences and Natural Resources, 2018, 11(2): 40-44 (0)
[3]
Boy J P, Chao B F. Time-Variable Gravity Signal during the Water Impoundment of China's Three-Gorges Reservoir[J]. Geophysical Research Letters, 2002, 29(24): 53-1-4 (0)
[4]
王天祥. 基于时序InSAR技术的地表形变监测研究——以兰州地区为例[D]. 兰州: 兰州大学, 2015 (Wang Tianxiang. Surface Deformation Monitoring Research Based on Time Series InSAR Technology——A Case Study of Lanzhou[D]. Lanzhou: Lanzhou University, 2015) (0)
[5]
瞿伟, 高源, 陈海禄, 等. 利用GPS高精度监测数据开展青藏高原现今地壳运动与形变特征研究进展[J]. 地球科学与环境学报, 2021, 43(1): 182-204 (Qu Wei, Gao Yuan, Chen Hailu, et al. Review on Characteristics of Present Crustal Tectonic Movement and Deformation in Qinghai-Tibet Plateau, China Using GPS High Precision Monitoring Data[J]. Journal of Earth Sciences and Environment, 2021, 43(1): 182-204) (0)
[6]
Tangdamrongsub N, Han S C, Jasinski M F, et al. Quantifying Water Storage Change and Land Subsidence Induced by Reservoir Impoundment Using GRACE, Landsat, and GPS Data[J]. Remote Sensing of Environment, 2019, 233: 111385 DOI:10.1016/j.rse.2019.111385 (0)
[7]
Wang L S, Chen C, Du J S, et al. Detecting Seasonal and Long-Term Vertical Displacement in the North China Plain Using GRACE and GPS[J]. Hydrology and Earth System Sciences, 2017, 21(6): 2905-2922 DOI:10.5194/hess-21-2905-2017 (0)
[8]
Landerer F W, Flechtner F M, Save H, et al. Extending the Global Mass Change Data Record: GRACE Follow-On Instrument and Science Data Performance[J]. Geophysical Research Letters, 2020, 47(12) (0)
[9]
郭飞霄, 谷延超, 崔阳, 等. 联合GPS和GRACE数据分析陕甘宁地区地表垂直形变[J]. 大地测量与地球动力学, 2020, 40(7): 697-703 (Guo Feixiao, Gu Yanchao, Cui Yang, et al. Analysis of Vertical Deformation in Shaanxi-Gansu-Ningxia Area by GPS and GRACE Data[J]. Journal of Geodesy and Geodynamics, 2020, 40(7): 697-703) (0)
[10]
钟敏, 段建宾, 许厚泽, 等. 利用卫星重力观测研究近5年中国陆地水量中长空间尺度的变化趋势[J]. 科学通报, 2009, 54(9): 1290-1294 (Zhong Min, Duan Jianbin, Xu Houze, et al. Trend of China Land Water Storage Redistribution at Medi-and Large-Spatial Scales in Recent Five Years by Satellite Gravity Observations[J]. Chinese Science Bulletin, 2009, 54(9): 1290-1294) (0)
[11]
涂梦昭, 刘志锋, 何春阳, 等. 基于GRACE卫星数据的中国地下水储量监测进展[J]. 地球科学进展, 2020, 35(6): 643-656 (Tu Mengzhao, Liu Zhifeng, He Chunyang, et al. Research Progress of Groundwater Storage Changes Monitoring in China Based on GRACE Satellite Data[J]. Advances in Earth Science, 2020, 35(6): 643-656) (0)
[12]
Abdelazim N, University C, Bekhit H, et al. Operation of the Grand Ethiopian Renaissance Dam: Potential Risks and Mitigation Measures[J]. Journal of Water Management Modeling, 2020(10) (0)
[13]
Wahr J, Molenaar M, Bryan F. Time Variability of the Earth's Gravity Field: Hydrological and Oceanic Effects and Their Possible Detection Using GRACE[J]. Journal of Geophysical Research: Solid Earth, 1998, 103(B12): 30205-30229 DOI:10.1029/98JB02844 (0)
[14]
Bevis M, Melini D, Spada G. On Computing the Geoelastic Response to a Disk Load[J]. Geophysical Journal International, 2016, 205(3): 1804-1812 DOI:10.1093/gji/ggw115 (0)
[15]
Dam T, Wahr J, Lavallée D. A Comparison of Annual Vertical Crustal Displacements from GPS and Gravity Recovery and Climate Experiment(GRACEover Europe[J]. Journal of Geophysical Research: Solid Earth, 2007, 112(B3) (0)
[16]
Yi S, Song C Q, Wang Q Y, et al. The Potential of GRACE Gravimetry to Detect the Heavy Rainfall-Induced Impoundment of a Small Reservoir in the Upper Yellow River[J]. Water Resources Research, 2017, 53(8): 6562-6578 DOI:10.1002/2017WR020793 (0)
[17]
Sims P H, Lentati L, Alexander P, et al. Contamination of the Epoch of Reionization Power Spectrum in the Presence of Foregrounds[J]. Monthly Notices of the Royal Astronomical Society, 2016, 462(3): 3069-3093 DOI:10.1093/mnras/stw1768 (0)
[18]
尼胜楠, 陈剑利, 李进, 等. 利用GRACE卫星时变重力场监测长江、黄河流域水储量变化[J]. 大地测量与地球动力学, 2014, 34(4): 49-55 (Ni Shengnan, Chen Jianli, Li Jin, et al. Terrestrial Water Storage Change in the Yangtze and Yellow River Basins from GRACE Time-Variable Gravity Measurements[J]. Journal of Geodesy and Geodynamics, 2014, 34(4): 49-55) (0)
Deformation Prediction of GERD Reservoir Based on Spherical Harmonic Transform
LIU Huiling1     CHEN Yu1     WANG Xihe1     
1. College of Electronics and Information Engineering, Sichuan University, 24 South-Section 1 of Yihuan Road, Chengdu 610065, China
Abstract: This study proposes a new method for predicting surface deformation caused during GERD storage based on spherical harmonic transform. Firstly, it is assumed that the GERD reservoir is filled up during 5, 10 and 15 years respectively, and the area, volume and water level of the reservoir are obtained for each month of the storage process using SRTM-DEM data. The surface deformation caused by the corresponding storage process is calculated after spherical harmonic expansion and the accuracy of the results is increased by increasing the upper limit of the maximum harmonic degree. The estimated vertical displacement velocity induced for each of the three filling strategies are 11.8 mm/a, 5.9 mm/a and 3.9 mm/a, respectively. The induced east-west displacement velocity are 1.8 mm/a and 1.4 mm/a, 0.9 mm/a and 0.7 mm/a, 0.6 mm/a and 0.4 mm/a, respectively. The induced south-north displacement velocity are 2.6 mm/a and 4.4 mm/a, 1.3 mm/a and 2.2 mm/a, 0.8 mm/a and 1.5 mm/a, respectively.
Key words: spherical harmonic transform; GERD; vertical deformation; horizontal deformation; time series