王希禾， 陳 雨， 劉慧玲

(四川大學(xué)電子信息學(xué)院， 成都 610065)

陸地水儲量的大幅度變化會使其周邊地表產(chǎn)生相應(yīng)的形變，嚴重時甚至?xí)鸹?、地震等地質(zhì)災(zāi)害的發(fā)生，因此對地表形變的監(jiān)測及預(yù)測是至關(guān)重要的課題。由負荷響應(yīng)理論可知，地表質(zhì)量載荷的變化會使地球表面產(chǎn)生相應(yīng)的位移。地球表層大氣壓力、海洋中的潮汐和非潮汐波動以及陸地上水、雪和冰的分布都是顯著的地表質(zhì)量載荷，其中，對于陸地水儲量變化模型的研究較少[1]。

重力衛(wèi)星技術(shù)是研究水儲量變化的一項重要技術(shù)，自重力恢復(fù)與氣候試驗(gravity recovery and climate experiment，GRACE)重力衛(wèi)星任務(wù)發(fā)布以來，吸引了眾多中外學(xué)者利用其重力場數(shù)據(jù)進行大地測量及地球物理相關(guān)學(xué)科的研究[2-6]。Wahr等[7]利用GRACE衛(wèi)星重力數(shù)據(jù)獲取地球時變重力場信息，并推演出陸地水變化的理論及方法，為后續(xù)相關(guān)研究奠定了基礎(chǔ)。但是使用GRACE衛(wèi)星數(shù)據(jù)反演水儲量變化的總體趨勢時，空間尺度范圍在300 km左右，存在信號平滑作用，導(dǎo)致反演精度降低[8-9]。GRACE全球覆蓋的特性使其在陸地水儲量較少的區(qū)域及陸地水較為豐富的區(qū)域均可以得到較好的應(yīng)用，但其產(chǎn)品的球諧系數(shù)通常僅有60階或96階，丟失了大部分的高階信息。

相較于GRACE來說，全球?qū)Ш叫l(wèi)星系統(tǒng)(global navigation satellite system，GNSS)測量的空間分辨率更高，可以通過連續(xù)監(jiān)測站點坐標的變化，反映陸地水和大氣等負荷遷移引起的地表形變[10-11]。全球定位系統(tǒng)(global positioning system，GPS)作為GNSS包含的四大導(dǎo)航系統(tǒng)之一，其數(shù)據(jù)常被用于反演地表質(zhì)量分布[12-14]。Kusche等[15]對全球GPS形變數(shù)據(jù)采用了物理激勵的正則化方法，保證了更高程度的穩(wěn)定反演。但GPS系統(tǒng)的安置及運行費用較高，從而限制了GPS網(wǎng)的分布密度，其觀測站點的數(shù)目通常與地區(qū)發(fā)展程度成正比，這種現(xiàn)象導(dǎo)致許多地區(qū)的水文變化無法監(jiān)測[16]。

計算負荷地表形變時，相同分辨率情況下負荷格林函數(shù)與球諧函數(shù)計算結(jié)果的精度一致，但在多站點大區(qū)域計算時，球諧函數(shù)的計算效率比負荷格林函數(shù)快100倍左右[17]?，F(xiàn)提出一種基于球諧變換解算地表三維形變(垂直形變、東西向形變、南北向形變)及其時間序列的方法，在僅有數(shù)字高程模型(digital elevation model，DEM)數(shù)據(jù)以及水位數(shù)據(jù)的情況下便可對相應(yīng)水體負荷造成的地表三維形變進行預(yù)測，從而提升全球柵格矩陣的空間分辨率。同時，以三峽水庫作為實例進行方法可行性的驗證，可將研究方法進一步推廣于GPS觀測站點稀少或不存在的地區(qū)，尤其是地表水文大幅遷移的區(qū)域中。例如，正在建設(shè)中的埃塞俄比亞文藝復(fù)興大壩，可在已知大壩注水速率的情況下推導(dǎo)其水位變化，繼而使用該方法預(yù)測地表水儲量變化引起的地表三維形變。

1 方法原理

1.1 平均水厚柵格矩陣的建立

通過創(chuàng)建平均水厚柵格矩陣模擬了一個理想化的以水體為質(zhì)量中心的載荷模型，用于反演水體載荷引起的地表形變。該模型首先創(chuàng)建一個分辨率為0.01°，尺寸為36 000×18 000的空白柵格矩陣。此空白矩陣的橫縱坐標分別對應(yīng)地理坐標系的經(jīng)緯度，每一柵格長度反映的地理距離約為1.1 km?；谘芯繀^(qū)域的DEM數(shù)據(jù)計算不同水平面以下的蓄水體積和各點水深，并求出研究區(qū)域不同水位所對應(yīng)的水體平均水厚度。在前述創(chuàng)建的空白矩陣中將研究區(qū)域水體所對應(yīng)的各點值設(shè)置為平均水厚度值，在此區(qū)域外的各點值均設(shè)為零，平均水厚柵格矩陣創(chuàng)建完畢。

1.2 球諧變換及負荷形變

將上述創(chuàng)建的平均水厚柵格矩陣擴展為球面諧波，采用了Wahr等[7]于1998年提出的球諧系數(shù)表達式，其計算公式為

(1)

負荷勒夫數(shù)被用來表示重力位對負荷引力位的響應(yīng)，簡稱負荷數(shù)。負荷勒夫數(shù)與地球自身內(nèi)部結(jié)構(gòu)及物質(zhì)的密度和彈性分布有關(guān)，一般表示為k、l、h，可用于表示地球在負荷作用下發(fā)生的形態(tài)變化[18]。地球表層大氣、陸地水以及海洋等負荷的質(zhì)量發(fā)生變化，均會使地球重力場產(chǎn)生相應(yīng)變化，形成負荷引力位，進而引起地球內(nèi)部質(zhì)量的重新分布。反演陸地水儲量變化以及地表形變時，由于球諧系數(shù)階次的增加，需要加載對應(yīng)階數(shù)的負荷勒夫數(shù)共同參與計算，一般采用的數(shù)據(jù)為Farrell[19]于1972年推導(dǎo)的經(jīng)典解。研究方法采用郭俊義等[20]求解的負荷勒夫數(shù)漸進值，該值比Farrell[19]的結(jié)論在精度上高出一個數(shù)量級。通常以等效水高的形式表示陸地水儲量變化，其計算公式為[21]

[ΔClmcos(mφ)+ΔSlmsin(mφ)]

(2)

式(2)中：ΔEWH(θ,φ)為等效水高；ρE為固體地球的平均密度，ρE=5 517 kg/m3；kl為l階處的負荷勒夫數(shù)。

地球?qū)嵸|(zhì)上是一個黏彈體，地表質(zhì)量載荷發(fā)生變化時，會引起相應(yīng)的地表彈性負荷形變[22]。在展開球諧系數(shù)計算地表水文負荷形變時，地表位移與地表質(zhì)量負荷之間存在一定的關(guān)聯(lián)，兩者相關(guān)公式為

(3)

(4)

式中：u(θ)為地表垂直方向的位移形變；v(θ)為地表水平方向的位移形變；me為負荷質(zhì)量；Pl(cosθ)為l階勒讓德函數(shù)；hl和ll為第l階的負荷勒夫數(shù)。

利用時變重力場反演地表質(zhì)量變化以及由此引起的負荷形變時，地表水文負荷導(dǎo)致的形變計算公式為[23]

(5)

(6)

(7)

式中：Δh(θ,φ)為垂直方向位移形變；Δn(θ,φ)、Δe(θ,φ)分別為北向以及東向位移形變，當?shù)貧し謩e向北、向東移動時，兩者為正值，反之為負。

1.3 計算時間序列

利用球諧函數(shù)將等效水高及3種地表形變展開之后，得到與平均水厚柵格矩陣同尺寸的柵格數(shù)據(jù)，并利用該柵格數(shù)據(jù)求解時間序列。依據(jù)Parseval定理可知，在空間域中，信號的平均功率等于其各個諧波分量的平均功率之和。該定理特性的最一般形式也可用于球諧系數(shù)的分析中，稱為Plancherel定理。利用球諧函數(shù)反演變換時，高次諧波的丟失將導(dǎo)致信號在空間域變得平滑和延展。

研究方法的重點不在于每個柵格的值或其平均值，而是以信號變換的總能量為著手點，盡可能減少球諧變換前后能量的損失，得到相對精確的時間序列。Δσ(θ,φ)以及ΔEWH(θ,φ)均表示水高，不同點在于Δσ(θ,φ)表示的質(zhì)量載荷為平均水厚，在其柵格數(shù)據(jù)中，只有水庫范圍內(nèi)的柵格單元取值為水庫的平均水厚度，水庫外的其余單元填充值為零；而在ΔEWH(θ,φ)的網(wǎng)格數(shù)據(jù)中，由于球面諧波被截斷造成的延展效應(yīng)，數(shù)值非零的單元分布在遠大于水庫實際范圍的區(qū)域中。經(jīng)實驗發(fā)現(xiàn)，在研究區(qū)域足夠廣泛的情況下，即信號能量的求和區(qū)域增大時，Δσ(θ,φ)與ΔEWH(θ,φ)的平均值更加接近，從而可以由Δσ(θ,φ)變換估算出相對精確的ΔEWH(θ,φ)數(shù)值。劃定研究區(qū)域的范圍后，通過對網(wǎng)格數(shù)據(jù)進行球諧域的變換，結(jié)合式(1)以及形變式(5)～式(7)，即可得到由于地表水體變化而引起的地表三維形變時間序列。

2 三峽庫區(qū)實例驗證

2.1 實驗區(qū)域及數(shù)據(jù)

為驗證研究方法的可靠性，選取三峽水庫作為實例驗證。如圖1所示，三峽水庫位于湖北省宜昌市三斗坪境內(nèi)，庫區(qū)蓄水工作于2003年開始，正常蓄水情況下庫區(qū)水位大約保持在145～175 m，總庫容量約為393×108m3[24]。三峽工程的建成帶來了巨大的防洪效益以及航運效益，水電站的建設(shè)更是體現(xiàn)了清潔能源開發(fā)的里程碑式發(fā)展。

圖1 三峽水庫位置示意圖Fig.1 The location of the Three Gorges Reservoir

使用的數(shù)字高程模型數(shù)據(jù)為NASA提供的SRTM-DEM數(shù)據(jù)(https://dwtkns.com/srtm30m/)，分辨率為30 m，獲取的DEM數(shù)據(jù)塊包含湖北宜昌到重慶的長江主干流。采用的水位數(shù)據(jù)來源于中國長江三峽集團有限公司(https://www.ctg.com.cn/)發(fā)布的水情信息，時間跨度為2011年1月—2018年12月。以2011年1月為基準按月獲取，共計96個水位數(shù)據(jù)。

2.2 等效水高及垂直形變驗證

GRACE重力衛(wèi)星發(fā)布的產(chǎn)品在球諧系數(shù)上僅有60階或96階，缺少中高階系數(shù)包含的細節(jié)信息，反映在空間圖像中是十分模糊的。利用球諧函數(shù)反變換解算等效水高及地表形變時，信號能量主要集中在水庫附近，但依舊存在一部分能量向水庫邊緣不斷延伸。當球諧階數(shù)增加時，信號能量會越來越集中，水庫位置包含的信息增加，從而使反變換得到的數(shù)據(jù)更精確。選取球諧截斷階數(shù)時，以反演的等效水高與輸入的等效水厚度的誤差作為判別依據(jù)。結(jié)果表明，球諧階數(shù)為800階時，96個數(shù)據(jù)中，誤差率最大為0.815%，最小為0.612%，平均誤差為0.742%。由于球諧系數(shù)的增加會降低計算效率，且800階時誤差率已經(jīng)達到較好的期望效果，因此后續(xù)計算時間序列時，選取的球諧截斷系數(shù)均為800階。

利用式(1)將創(chuàng)建好的平均水厚網(wǎng)格數(shù)據(jù)擴展至有限次的球面諧波，獲取相應(yīng)的球諧系數(shù)，對地表垂直形變進行反演。均以2011年1月為基準，對比計算2018年6月三峽庫區(qū)發(fā)生的地表形變。將球諧系數(shù)的截斷階次分別設(shè)置為60階、200階、600階、800階，觀測不同系數(shù)下垂直形變的空間分布狀況。圖2中，相較于2011年1月，該月水位高度下降，地殼垂直向上抬升，垂直形變?yōu)檎?。球諧截斷階數(shù)為60階時，由于信號能量的大量缺失，整個三峽庫區(qū)范圍的形變僅有0.6 mm左右。截斷階數(shù)逐漸增加時，高階系數(shù)的能量信息也被計算在內(nèi)，形變值的“分層”現(xiàn)象越發(fā)明顯，不同形變值的邊界也越發(fā)清晰，800階時三峽庫區(qū)中心形變值達到6.938 44 mm左右。

圖3(a)展示了800階球諧截斷階數(shù)下求得的垂直形變時間序列，每年的1—7月，水庫水位持續(xù)下降，地殼垂直向上移動，地表垂直形變呈現(xiàn)持續(xù)增大的趨勢。反之，8—12月份，水庫水位不斷上升時，地表垂直形變呈減小的趨勢。王偉等[25]利用負荷格林函數(shù)積分法計算了三峽庫區(qū)茅坪站水位變化引起的地殼垂直形變，并得出2011年1月—2015年6月的垂直形變時間序列，其觀測結(jié)果如圖3(b)所示。由于該參考文獻中茅坪站的水位數(shù)據(jù)與本文輸入的三峽水位數(shù)據(jù)一致，因此兩者結(jié)果具備可比性。比較圖3(a)、圖3(b)中2011年1月—2015年6月的曲線段，可以發(fā)現(xiàn)兩種方法得到的時間序列變化趨勢整體一致，只是在數(shù)值上存在差距。參考文獻中的地表垂直形變?yōu)?1～10 mm，研究方法得到的結(jié)果為-0.58～13.53 mm。出現(xiàn)該差距與本文方法采用較高的球諧階數(shù)進行反演有關(guān)，球諧系數(shù)的增加使解算結(jié)果中包含了更多高階能量信息，水庫所在網(wǎng)格區(qū)域的累加之和也隨之增大、更為精確。

黑色不規(guī)則形狀為當月水位數(shù)據(jù)下三峽水庫的流域橫截面圖2 三峽水位變化對地表垂直形變的 影響(2018年6月)Fig.2 The impact of water level changes in the Three Gorges on vertical surface deformation (June 2018)

圖3 三峽水位變化引起的地表垂直形變時間序列Fig.3 Time series of vertical surface deformation caused by water level changes in the Three Gorges

2.3 解算水平形變及其時間序列

有關(guān)三峽地區(qū)的水平形變研究相較于垂直形變較少，且多為綜合分析，未分解為東向形變及北向形變單獨研究，因此研究中只解算水平形變，不做對比。

對東向形變進行分析時，依舊以2018年6月的數(shù)據(jù)為例，圖4分別為球諧截斷階數(shù)為60階、200階、600階、800階時計算的東向形變分布圖。當球諧截斷階數(shù)增加時，計算的東向形變更集中于三峽水庫周圍，而不是擴散至更遠的區(qū)域，800階時庫區(qū)最大形變達到0.805 566 mm。水位下降時，東側(cè)網(wǎng)格負荷向東移動，位移值為正，西側(cè)網(wǎng)格負荷向西移動，位移值為負；水體負荷均由流域中心向東西兩側(cè)移動，地殼沿庫區(qū)流域向外拉伸。

設(shè)置球諧截斷系數(shù)為800階時，得到的三峽大壩的東向形變時間序列如圖5所示。2011—2018年東向位移的最大絕對差值約為0.77 mm，且每年的形變時間序列都呈先下降后上升的規(guī)律性變化趨勢，該趨勢與三峽水庫水位變化趨勢大致相同。

黑色不規(guī)則形狀為當月水位數(shù)據(jù)下三峽水庫的流域橫截面圖4 三峽水位變化對地表東向形變的影響(2018年6月)Fig.4 The impact of water level changes in the Three Gorges on eastward surface deformation (June 2018)

圖5 三峽水位變化引起的地表東向形變時間序列Fig.5 Time series of eastward surface deformation caused by water level changes in the Three Gorges

受三峽水位變化的影響產(chǎn)生的北向形變?nèi)鐖D6所示，與垂直形變及東向形變相似，隨著球諧截斷階數(shù)的增加，計算的形變分布圖更為精確。概括來說，以三峽水庫水體為中心，庫區(qū)兩岸均呈現(xiàn)出由庫區(qū)中心向外部擴散的趨勢。位于庫區(qū)北部的載荷向北移動，呈現(xiàn)正值，位于南部的載荷向南移動，為負值，且南部位移絕對值略大于北部位移。800階時，庫區(qū)北部載荷最大位移值為0.669 066 mm，南部載荷最大位移為2.091 570 mm。

以800階截斷系數(shù)計算三峽大壩的北向形變時間序列，得到的結(jié)果如圖7所示。2011—2018年計算出的北向形變時間序列范圍為-1.96～0.07 mm，且該時間序列表現(xiàn)出的年度規(guī)律與東向形變時間序列基本相同。

黑色不規(guī)則形狀為當月水位數(shù)據(jù)下三峽水庫的流域橫截面圖6 三峽水位變化對地表北向形變的影響(2018年6月)Fig.6 The impact of water level changes in the Three Gorges on northward surface deformation (June 2018)

圖7 三峽水位變化引起的地表北向形變時間序列Fig.7 Time series of northward surface deformation caused by water level changes in the Three Gorges

3 結(jié)論

提出一種利用球諧函數(shù)及負荷響應(yīng)解算地表三維形變的方法，并將反演過程的空間分辨率提高至0.01°。為驗證該方法的可行性，使用三峽水庫作為實例進行驗證，得出如下結(jié)論。

(1)首先利用反演的等效水高進行驗證，球諧系數(shù)為800階時誤差率約為0.742%，即反演精度能達到99.258%。然后計算三峽水庫的地表垂直形變時間序列，并與格林函數(shù)計算結(jié)果對比，二者在趨勢上一致，且研究方法由于網(wǎng)格數(shù)據(jù)精度提高、球諧截斷系數(shù)增加，使得能量更為集中，計算結(jié)果數(shù)值較大、更精確。從等效水高以及垂直形變兩個角度驗證方法可行性之后，依據(jù)負荷形變原理計算水平形變，由于以往對三峽地區(qū)水平形變的研究較少，且多將北向形變、東向形變綜合分析，因此未作對比。

(2)相較于GRACE重力衛(wèi)星來說，研究方法的局限是在水體負荷變化較為明顯的地區(qū)中應(yīng)用時更加有效，如儲存大量水體的水庫或湖泊，但優(yōu)點是提高了反演的空間分辨率。此外，假設(shè)有一在建水庫，GRACE與GPS均不能預(yù)測水庫蓄水引起的地表形變；但在知曉蓄水速率之后，可以使用該方法，通過對數(shù)字高程模型(digital elevation model，DEM)數(shù)據(jù)處理得到對應(yīng)水位數(shù)據(jù)，進而計算出研究區(qū)域由水體變化產(chǎn)生的形變數(shù)據(jù)。因此本方法對于預(yù)測及計算一些大型水域地表三維形變的研究來說具有一定的實際意義。