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

劉慧玲 陳 雨 王希禾

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

水壩通過調(diào)節(jié)水流，為人類和環(huán)境的各種需求提供了可靠水源[1]。埃塞俄比亞的文藝復興大壩(Grand Ethiopian Renaissance Dam，GERD)是現(xiàn)今非洲最大的大壩，預計其有效庫容可達74 km3[2]。GERD水庫蓄水將引發(fā)大量的水文遷移，可能導致周邊地區(qū)發(fā)生大規(guī)模的地表形變[3]，地表形變量超出一定界限會演變成為地質(zhì)災害，對人類生命和財產(chǎn)安全帶來嚴重損失[4]。

地表形變監(jiān)測網(wǎng)絡是監(jiān)測區(qū)域地表形變的有效方法之一[5]，但該方法不適用于GERD水庫區(qū)域，因為GERD附近240 km范圍內(nèi)沒有可用的GPS站點，距離GERD 142 km的ASOS站點已于2017年停止工作。衛(wèi)星遙感可以進行長時間、高覆蓋率的觀測，能連續(xù)監(jiān)測水庫蓄水和地表覆蓋變化，許多研究者利用多種遙感數(shù)據(jù)相結(jié)合成功估算了特定地區(qū)的地表形變[6-7]。GRACE以及最新的GRACE-FO三大數(shù)據(jù)處理機構(gòu)(JPL、CSR、GFZ)將地表引力場變化通過球諧變換，得到隨時間變化的GRACE重力數(shù)據(jù)的球諧系數(shù)(spherical harmonic coefficients，簡稱SHC)形式(即Level 2產(chǎn)品)[8]，利用該產(chǎn)品并引入負荷勒夫數(shù)(Love numbers)可反演得到相應重力變化引起的地表負荷形變[9]。但由于GRACE的分辨率只有200～300 km，其時變重力場的SHC的階數(shù)只能展開到60～96階[10]，受制于分辨率的影響，大部分研究中基于GRACE反演得到的地下水儲量變化與實際結(jié)果的相關系數(shù)僅為0.6～0.7[11]。所以，通常利用GPS數(shù)據(jù)或其他水文模型對GRACE反演后的結(jié)果進行校正。

根據(jù)以上分析，本文基于球諧變換和Love numbers提出一種預測GERD水庫蓄水引起的地表形變的方法。首先假設3種蓄水方案：短期5 a、中期10 a、長期15 a內(nèi)注滿74 km3的水庫，并根據(jù)數(shù)字高程模型(DEM)分別對3種蓄水方案中每月蓄水表面積、體積和平均水高進行數(shù)字擬合。然后以平均水高為基礎創(chuàng)建一個7 200×3 600的柵格網(wǎng)絡，該柵格水庫區(qū)域的各點設為平均水高值，水庫區(qū)域外各點設為0。通過對比等效水高(equivalent water height，EWH)和平均水高分析高階球諧截斷所造成的信號能量損失，再以2 300階的球諧變換系數(shù)為基礎，通過引入Love numbers計算水庫及周邊區(qū)域的垂直和水平形變。

1 數(shù)據(jù)準備和研究方法

1.1 研究區(qū)域及數(shù)據(jù)準備

GERD位于埃塞俄比亞的藍色尼羅河上，地處蘇丹和埃及的上游，于2011-04開始建設，計劃正常水高為640 m，水庫蓄水體積為74 km3，約為大壩所在地年均流量的1.5倍[12]。

采用美國宇航局的SRTM-DEM數(shù)據(jù)(https:∥dwtkns.com/srtm30m/)，空間分辨率為30 m。根據(jù)GERD的實際位置先選取32°50′～48°9′E、3°21′～ 15°1′N的區(qū)域，然后將GERD的位置作為傾瀉點，從上述區(qū)域中劃分出本文的研究區(qū)域——上藍色尼羅河流域。

1.2 研究方法

本文的數(shù)據(jù)處理流程分為5個步驟(圖1)：1)利用SRTM-DEM數(shù)據(jù)對研究區(qū)域在不同蓄水策略下提取的水體進行分析，得出蓄水過程中每月水庫的蓄水表面積、體積和平均水高；2)構(gòu)建全局柵格圖，將水庫區(qū)域內(nèi)的柵格值設為當月平均水高，其余區(qū)域設為0；3)將此柵格網(wǎng)絡球諧擴展為最大諧波度(Nmax)分別為60、700、1 500、2 300的SHC；4)利用SHC計算EWH，先通過比較EWH和平均水高，分析高階SHC對信號能量損失的影響，再利用Nmax=2 300階時的SHC計算地表形變；5)提取垂直、東、西、南和北向形變的時間序列圖。

圖1 數(shù)據(jù)處理過程Fig.1 Data processing

1.2.1 計算水庫蓄水表面積、體積和平均水高

首先使用一組高程值間隔1 m的水平表面截取研究區(qū)域的SRTM-DEM柵格圖，獲得指定存儲層中每個水平表面下的體積，再構(gòu)建以體積為自變量的體積/水平面高程值擬合函數(shù)。如果以74 km3作為水庫額定容量，將蓄水年限分為5 a、10 a和15 a分別進行計算，對于3種年限，每月分別需蓄水1.23 km3、0.61 km3和0.41 km3。利用體積/水平面高程值擬合函數(shù)計算不同蓄水策略下每月蓄水體積對應的水平表面高程值；然后使用ArcGIS工具計算水庫存儲層中每月水平表面高程值下的蓄水表面面積；最后用月蓄水體積除以蓄水表面積得到每月平均水高。

1.2.2 設計平均水高柵格及球諧展開

本文構(gòu)建了一個分辨率為0.05°×0.05°的包含全球經(jīng)緯度的柵格網(wǎng)絡。水庫范圍內(nèi)的柵格值設為當月平均水高，其余柵格值設為0；然后將此全局柵格網(wǎng)絡球諧展開至指定Nmax的SHC，并根據(jù)Wahr等[13]的方法求出每月變化量。

1.2.3 計算EWH和地表形變

利用球諧函數(shù)計算EWH和地表變形時需引入Love numbers，每月SHC變化量結(jié)合Love numbers可計算由水文質(zhì)量載荷變化引起的EWH變化和3D位移。ΔEWH的具體計算公式可參考文獻[14]，3D位移的具體計算公式可參考文獻[15-16]。

2 結(jié)果分析

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蓄水策略中對應的數(shù)據(jù)為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個月的平均水高圖，水庫內(nèi)的柵格值為34 m(藍色區(qū)域)，水庫外柵格值為0(黃色區(qū)域)。在此蓄水階段中，根據(jù)不同Nmax計算各自的EWH，結(jié)果如圖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定理通常用于描述傅里葉變換的統(tǒng)一性，但該特性的最一般形式應更恰當?shù)乇环Q為Plancheral定理。在數(shù)學中，Plancheral定理是諧波分析的結(jié)果，所以能量定理也適用于諧波分析。當研究區(qū)域為面積較大的水庫或湖時，區(qū)域地表質(zhì)量變化的主要影響因素是水庫或湖中水的儲量變化。研究中常用ΔEWH(θ,φ)來表示陸地水儲量的變化[18]。

由于球諧擴展受Nmax的限制，高階諧波被截斷導致信號在空間域被平滑和擴展，所以利用球諧函數(shù)計算ΔEWH和3D位移是一種數(shù)值再分配模式下的計算，計算后的區(qū)域比研究的水庫區(qū)域要大得多。圖4(b)顯示了由于高于60階的諧波信號被截斷導致信號能量被平滑、衰減和擴展至水庫邊界之外的現(xiàn)象，其柵格單元中EWH的范圍為-0.164～0.123 m。

隨著Nmax的逐漸增大，受截斷效應影響的信號逐漸減少，EWH信號向GERD水庫中心集中，呈現(xiàn)出水庫中心區(qū)域的EWH值明顯大于周邊區(qū)域的特性。當Nmax=700，水庫中心區(qū)域的EWH最大值為16.5 m，EWH的范圍為-3.2～16.5 m；當Nmax=1 500時，水庫中心區(qū)域的EWH最大值增加至30.03 m，EWH的范圍增加至-3.83～30.03 m；當Nmax=2 300時，已顯示出了EWH沿GERD水庫區(qū)域分布的特性，EWH的范圍為-8.28～38.3 m(由于吉布斯現(xiàn)象會出現(xiàn)少數(shù)大于34 m和小于0的柵格值)。

利用Nmax=2 300時的SHC計算5 a蓄水策略下的平均水高與EWH的時間序列。經(jīng)誤差分析，兩者的平均相對誤差僅有1.6%，因此選定Nmax=2 300階計算后續(xù)地表形變。

2.3 地表形變及其時間序列

2.3.1 垂直形變

計算5 a策略中第30個月的垂直形變，結(jié)果如圖5(a)所示。從圖中可以看出，柵格中所有質(zhì)量載荷點的垂直位移分量都向水庫質(zhì)心運動，為負值。該月垂直形變位移的范圍為-91.13～-0.49 mm,且越靠近水庫區(qū)域垂直形變位移量越大，水庫中心區(qū)域的垂直形變位移范圍為-91.13～-55.04 mm。

圖5 5 a蓄水策略下第30個月的形變結(jié)果Fig.5 Deformation results of the 30th month under 5 a filling strategy

2.3.2 水平形變

水庫區(qū)域質(zhì)量載荷點的位置決定了該點的位移方位(向北或向南)，由于水庫北部面積大于南部，因此質(zhì)量載荷在南部單元上產(chǎn)生的拉力大于北部。在5 a策略下第30個月的蓄水階段中，南向位移的最大值約為9 mm，北向位移的最大值約為37 mm。另外，由于水庫北部區(qū)域較寬，因此東部和西部的單元也包含了向南移的分量，如圖5(b)所示。

東西向位移與南北向位移相似，水庫西岸的載荷點運動趨勢向東，表現(xiàn)為正值，水庫東岸的載荷點則向西移動，表現(xiàn)為負值，如圖5(c)所示。由于GERD水庫東西方向分布較為均勻，所以呈現(xiàn)的東西向位移量也較為均衡，東向位移的最大值為9.31 mm，西向位移的最大值為10.59 mm。

2.3.3 地表形變的時間序列

依據(jù)圖5的結(jié)果分析垂直向、南北向和東西向的形變區(qū)域，并計算3種蓄水策略下5個方向的時間序列，結(jié)果如圖6所示，圖中，V表示垂直形變，E表示東向形變，W表示西向形變，N表示北向形變，S表示南向形變。

圖6 3種蓄水策略下水平形變的時間序列Fig.6 Time series of horizontal deformation under three strategies

使用恒定體積的水蓄滿水庫時，蓄水初期和中后期將產(chǎn)生較大的垂直位移, 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個月的東西向位移網(wǎng)格圖顯示出一定的對稱性；這種對稱性在東西向位移的時間序列中也有體現(xiàn)，如圖6(b)所示。本文認為，這種對稱性是由于水庫東西狹窄的形狀所引起的。5 a策略的東向年平均位移量約為1.8 mm，最快的變化發(fā)生在蓄水的第9個月，達到0.42 mm/月，在蓄水的最后時期，移動速率降到0.1 mm/月；10 a策略下最快和最慢的東向位移速率為0.35 mm/月和0.07 mm/月，分別發(fā)生在蓄水期的第19個月和第80個月前后；而在15 a策略下，東向形變最大速率為0.25 mm/月，最小速率為0.03 mm/月。

GERD水庫南部和北部質(zhì)量分布不平衡導致南部和北部的位移分量趨勢不同，北部較大的區(qū)域可能會分散質(zhì)量負荷，北部質(zhì)量載荷單元的移動小于南部。南部和北部網(wǎng)格單元的位移時間序列如圖6(c)所示，負值表示北部的網(wǎng)格單元向南移動，正值表示南部的網(wǎng)格單元向北移動，兩者均向水庫的質(zhì)心移動。對于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)，存儲層將產(chǎn)生約59 mm的垂直位移、22.5 mm的北向位移、13 mm的南向位移、9.2 mm的東向位移和6.7 mm的西向位移。

表1 不同蓄水策略下的年位移變化、總位移變化和最大位移變化

3 結(jié) 語

本文提出了一種基于球諧變換預測GERD水庫地表位移分量的方法。分別假定5 a、10 a、15 a蓄滿水庫，則由此引發(fā)的地表垂直形變速率分別約為11.8 mm/a、5.9 mm/a、3.9 mm/a；引發(fā)的南向和北向位移速率分別約為2.6 mm/a和4.4 mm/a、1.3 mm/a和2.2 mm/a、0.8 mm/a和1.5 mm/a；引發(fā)的東向和西向位移速率分別約為1.8 mm/a和1.4 mm/a、0.9 mm/a和0.7 mm/a、0.6 mm/a和0.4 mm/a。同時實驗發(fā)現(xiàn)，蓄水速度越快，月/年形變量越大。

其次，通過60階、700階、1 500階和2 300階的球諧展開對比發(fā)現(xiàn)，高階次球諧波的截斷會導致邊緣信息丟失和空間信息的平滑和擴展。為了減少該效應對本文結(jié)果的影響，經(jīng)過對比，最后選定2 300階的截斷來完成實驗。對能量守恒定理的分析和輸入的平均水高與所求EWH僅1.6%的誤差率均證明了本實驗的合理性。