基于球諧變換的GERD水庫形變預測

劉慧玲 陳 雨 王希禾

1 四川大學電子信息學院，成都市一環路南一段24號，610065

水壩通過調節水流，為人類和環境的各種需求提供了可靠水源[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個月的平均水高與EWHFig.4 Average water height and EWH of the 30thmonth 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 不同蓄水策略下的年位移變化、總位移變化和最大位移變化

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%的誤差率均證明了本實驗的合理性。