基于球諧變換及負荷響應解算地表三維形變

王希禾， 陳 雨， 劉慧玲

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

陸地水儲量的大幅度變化會使其周邊地表產生相應的形變，嚴重時甚至會引起滑坡、地震等地質災害的發生，因此對地表形變的監測及預測是至關重要的課題。由負荷響應理論可知，地表質量載荷的變化會使地球表面產生相應的位移。地球表層大氣壓力、海洋中的潮汐和非潮汐波動以及陸地上水、雪和冰的分布都是顯著的地表質量載荷，其中，對于陸地水儲量變化模型的研究較少[1]。

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

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

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

1 方法原理

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

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

1.2 球諧變換及負荷形變

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

(1)

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

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

(2)

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

地球實質上是一個黏彈體，地表質量載荷發生變化時，會引起相應的地表彈性負荷形變[22]。在展開球諧系數計算地表水文負荷形變時，地表位移與地表質量負荷之間存在一定的關聯，兩者相關公式為

(3)

(4)

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

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

(5)

(6)

(7)

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

1.3 計算時間序列

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

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

2 三峽庫區實例驗證

2.1 實驗區域及數據

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

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

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

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

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

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

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

黑色不規則形狀為當月水位數據下三峽水庫的流域橫截面圖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 解算水平形變及其時間序列

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

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

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

黑色不規則形狀為當月水位數據下三峽水庫的流域橫截面圖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

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

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

黑色不規則形狀為當月水位數據下三峽水庫的流域橫截面圖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 結論

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

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

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