劉慧玲 陳 雨 王希禾

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

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

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

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

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

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

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

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

1.2 研究方法

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

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

1.2.1 計(jì)算水庫(kù)蓄水表面積、體積和平均水高

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

1.2.2 設(shè)計(jì)平均水高柵格及球諧展開(kāi)

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

1.2.3 計(jì)算EWH和地表形變

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

2 結(jié)果分析

2.1 每月水庫(kù)蓄水表面積、體積和平均水高

圖2是5 a、10 a和15 a蓄水策略的模擬實(shí)驗(yàn)中水庫(kù)蓄水表面積和平均水高隨時(shí)間變化曲線。由圖可見(jiàn)，3種蓄水策略下GERD水庫(kù)平均水高的月平均變化量分別為1.63 m、0.88 m和0.61 m，水庫(kù)蓄水表面積的月平均變化量分別為28.66 km2、14.59 km2和9.72 km2。

圖2 3種策略下每月平均水高與水庫(kù)表面積Fig.2 Monthly average water height and area under three strategies

圖3為3種蓄水策略在不同蓄水階段的體積、蓄水表面積和水位的變化。分析可知，5 a策略的第1個(gè)月蓄水體積為1.23 km3，其水面面積為82.52 km2，平均水高為14.96 m；10 a和15 a蓄水策略中對(duì)應(yīng)的數(shù)據(jù)為0.61 km3、46.05 km2、13.32 m和0.41 km3、34.11 km2、12.06 m。當(dāng)5 a 蓄水策略完成時(shí)，10 a蓄水策略完成50%，15 a蓄水策略完成約33%；當(dāng)10 a蓄水策略完成時(shí)，15 a蓄水策略完成約66%；當(dāng)整個(gè)蓄水完成時(shí)，會(huì)形成一個(gè)體積約為73.99 km3、面積約為1 773.2 km2、平均水高約為41.72 m的水庫(kù)。

圖3 3種策略模擬水庫(kù)容量、水位、面積Fig.3 Simulated reservoir volume, area and water height using three strategies

2.2 基于不同Nmax的EWH對(duì)比與分析

圖4(a)為5 a蓄水策略中第30個(gè)月的平均水高圖，水庫(kù)內(nèi)的柵格值為34 m(藍(lán)色區(qū)域)，水庫(kù)外柵格值為0(黃色區(qū)域)。在此蓄水階段中，根據(jù)不同Nmax計(jì)算各自的EWH，結(jié)果如圖4(b)～4(e)所示。

圖4 5 a蓄水策略下第30個(gè)月的平均水高與EWHFig.4 Average water height and EWH of the 30thmonth under 5 a filling strategy

Parseval能量守恒定理表明，信號(hào)在空間域的平均功率等于時(shí)域中各次諧波分量的平均功率之和[17]。雖然Parseval定理通常用于描述傅里葉變換的統(tǒng)一性，但該特性的最一般形式應(yīng)更恰當(dāng)?shù)乇环Q為Plancheral定理。在數(shù)學(xué)中，Plancheral定理是諧波分析的結(jié)果，所以能量定理也適用于諧波分析。當(dāng)研究區(qū)域?yàn)槊娣e較大的水庫(kù)或湖時(shí)，區(qū)域地表質(zhì)量變化的主要影響因素是水庫(kù)或湖中水的儲(chǔ)量變化。研究中常用ΔEWH(θ,φ)來(lái)表示陸地水儲(chǔ)量的變化[18]。

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

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

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

2.3 地表形變及其時(shí)間序列

2.3.1 垂直形變

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

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

2.3.2 水平形變

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

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

2.3.3 地表形變的時(shí)間序列

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

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

使用恒定體積的水蓄滿水庫(kù)時(shí)，蓄水初期和中后期將產(chǎn)生較大的垂直位移, 5 a策略下垂直位移變化的最快速率和最慢速率分別為1.2 mm/月和0.4 mm/月；10 a策略下對(duì)應(yīng)速率降為1.0 mm/月和0.35 mm/月；15 a策略下僅為0.45 mm/月和0.2 mm/月。這3種策略的年平均垂直位移量呈遞減的趨勢(shì)，分別約為11.8 mm、5.9 mm、3.9 mm。

如圖5(c)所示，5 a蓄水策略下第30個(gè)月的東西向位移網(wǎng)格圖顯示出一定的對(duì)稱性；這種對(duì)稱性在東西向位移的時(shí)間序列中也有體現(xiàn)，如圖6(b)所示。本文認(rèn)為，這種對(duì)稱性是由于水庫(kù)東西狹窄的形狀所引起的。5 a策略的東向年平均位移量約為1.8 mm，最快的變化發(fā)生在蓄水的第9個(gè)月，達(dá)到0.42 mm/月，在蓄水的最后時(shí)期，移動(dòng)速率降到0.1 mm/月；10 a策略下最快和最慢的東向位移速率為0.35 mm/月和0.07 mm/月，分別發(fā)生在蓄水期的第19個(gè)月和第80個(gè)月前后；而在15 a策略下，東向形變最大速率為0.25 mm/月，最小速率為0.03 mm/月。

GERD水庫(kù)南部和北部質(zhì)量分布不平衡導(dǎo)致南部和北部的位移分量趨勢(shì)不同，北部較大的區(qū)域可能會(huì)分散質(zhì)量負(fù)荷，北部質(zhì)量載荷單元的移動(dòng)小于南部。南部和北部網(wǎng)格單元的位移時(shí)間序列如圖6(c)所示，負(fù)值表示北部的網(wǎng)格單元向南移動(dòng)，正值表示南部的網(wǎng)格單元向北移動(dòng)，兩者均向水庫(kù)的質(zhì)心移動(dòng)。對(duì)于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為不同蓄水策略下的年位移變化、位移總變化和最大位移變化?？梢钥闯觯钏乃俣仍娇?，年變形量越大。總體而言，在水庫(kù)蓄水完成后(注滿74 km3)，存儲(chǔ)層將產(chǎn)生約59 mm的垂直位移、22.5 mm的北向位移、13 mm的南向位移、9.2 mm的東向位移和6.7 mm的西向位移。

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

3 結(jié) 語(yǔ)

本文提出了一種基于球諧變換預(yù)測(cè)GERD水庫(kù)地表位移分量的方法。分別假定5 a、10 a、15 a蓄滿水庫(kù)，則由此引發(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。同時(shí)實(shí)驗(yàn)發(fā)現(xiàn)，蓄水速度越快，月/年形變量越大。

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