利用質量遷移指紋和GRACE數據解算地心運動

(1. 福州大學 空間數據挖掘與信息共享教育部重點實驗室，福建 福州350108；2. 閩江學院 海洋學院，福建 福州 350108；3.代爾夫特理工大學 地球科學與遙感學院，荷蘭 南荷蘭省 代爾夫 2600AA)

地表質量，包括大氣、海洋、冰川、陸地水等，不斷交換、遷移，致使固體地球形心(center of figure， CF)與地球質心(center of mass， CM)時刻處于相對運動中。采用國際地球自轉和參考框架服務機構(the international earth rotation and reference system service，IERS)的建議，將CM定義為地心，CM相對于CF的運動即為地心運動[1]。

地心運動分量ΔX,ΔY,ΔZ與地球重力場模型的一階項系數ΔC10,ΔC11,ΔS11相對應，即[2]：

(1)

其中：a為地球平均半徑；k1為一階彈性負載勒夫數，其值取決于坐標系的選取，例如，在CF坐標系中，k1=0.021[2]。以下為了簡潔，省去公式中的Δ符號，但均應理解為時變信號。

重力恢復與氣候實驗(gravity recovery and climate experiment，GRACE)重力衛星的成功，空前提高了地球時變重力場模型的解算精度和分辨率，首次實現人類對地球系統中大規模質量遷移的直接觀測，加深了對氣候變化影響的理解[3]。GRACE任務采用低-低衛星跟蹤衛星方式獲取觀測數據，兩顆衛星始終環繞CM運動，因而無法獲取CM相對于CF的運動信息，即無法獲得地心運動數據。因此，基于GRACE數據解算的地球重力場模型不包含一階項系數。

為填補GRACE重力場模型的一階項系數，可采用包含該信息的其他衛星數據，如衛星激光測距(satellite laser ranging，SLR)數據。但由于目前SLR地面觀測站數量稀疏，且集中于北半球中緯度地區，所得結果誤差較大。Blewitt等[4]利用全球分布的GPS地面站點探測地表形變，反演地表質量遷移，解算地心運動。該方法雖得到一定推廣[5]，但由于其無法解算地心運動的趨勢，在GRACE重力場模型用戶中未得到廣泛應用。目前GRACE數據處理中使用最為廣泛的方法是由Swenson等[6]首先提出的基于GRACE數據和洋底壓力(ocean bottom pressure，OBP)模型的GRACE-OBP方法。其基本思路是利用洋底壓力模型獲取海洋各點在CF坐標系中的質量變化，利用GRACE數據獲取這些點在CM坐標系中的變化，認為二者差值是由于一階項缺失所導致的，從而可以利用最小二乘法反算一階項。隨后，Sun等[7]對該方法進行了改進，使之可以同時解算C20；并利用模擬數據確定了該方法的最佳參數組合以便獲得最優的地心運動周期項估計[1]；通過結合輸入數據的誤差信息得到統計最優的地心運動[8]。GRACE-OBP方法依賴于精度未知的OBP，最終得到的地心運動結果精度難以評定。

近年來，指紋法在全球及區域海平面變化研究中得到初步應用。Rietbroek等[9]通過建立指紋數據集，結合GRACE和衛星測高數據，解算了全球及區域海平面變化。Sun等[10]利用GRACE數據和指紋數據集解釋了C20變化。與GRACE-OBP方法相比，指紋法不直接依賴地球物理模型，且可以給出地球系統各個分量的貢獻。本研究旨在利用GRACE和指紋數據集解算地心運動。

1 算法介紹

由于自吸引效應(self-attraction and loading effect)，地表任一點的質量變化都會導致全球性的質量遷移，從而造成全球重力場變化。本研究將點P(θ，φ)處的單位質量變化所導致的全球性質量重新分布模式定義為指紋fM，可以在頻率域內展開，以球諧系數形式表達。

假設某點或某區域發生單位質量變化，可利用海平面方程(參考Tamisiea等[11]公式(1)～(5))計算相應的fM。當該點或該區域質量變化為αfM而非單位質量時，所導致的全球質量變化可表示為αfMfM。將各點質量遷移所導致的全球質量變化疊加即可得到GM：

(2)

其中：FM為fM的集合；m為指紋個數；n為每個指紋展開至某一階次所包含的球諧系數的個數；α為指紋的系數向量。由于質量系數可由重力場模型的斯托克斯(Stokes)系數通過簡單轉換來求解[13]，即：

(3)

其中：ρ為地球平均密度；l為階數；kl為l階彈性勒夫數。因此，公式(2)可轉化為：

(4)

其中，F為f的集合，G由GRACE月重力場模型提供，因此一階項缺失，f則為一套包含一階項的球諧系數。另外，由于GRACE重力場模型的C20不準確，本研究假設其未知，與一階項同時解算，因此公式(4)應寫成：

(5)

其中T為截斷矩陣：

(6)

其作用是消掉f中的C10，C11，S11和C20，從而與G中的系數保持一致。為便于理解，假設只需要5個指紋：1個ANT(南極，Antarctica)指紋，1個GRE(格陵蘭島，Greenland)指紋，1個GLA(陸地冰川，glaciers)指紋，1個TWS(陸地水，terrestrial water storage)指紋和1個GIA(冰川均衡調整，glacial isostatic adjustment)指紋，則公式(5)可以進一步展開為：

(7)

當殘差G-TFα較小時，可以認為對地球系統的區塊或模式劃分比較合理。注意，α通過最小二乘法求解，理論上G中只需包含4個系數即可(僅需解算4個未知數)，但使用較多的系數有利于減少系數誤差對于結果的影響，使得求解穩定精確。同時應注意，某一時刻f中的所有系數均對應同一個αf，也即利用C21及以上系數所解算的αf也適用于一階項和C20。因此：

(8)

利用GRACE重力場模型月解(不包含一階項，去除C20)，可以得到該模型的一階項和C20。利用該方法，每月解算一組，最終形成時間序列。不僅可以得到總的一階項時序，還可以得到各個貢獻源的分量，即：

(9)

(10)

式中，mk表示第k個月份，前后兩個月份分別表示為mk-1和mk+1。

圖1展示了對GIA相關指紋所采用的時域約束。注意圖1中貢獻源僅為4個，共包含1個ANT，1個GRE、1個TWS和2個GIA共計5個指紋。在實際解算過程中，為了使所添加的約束有效，需要令其有較高的權重。由于指紋系數一般小于10-3量級，因此采用圖1所示約束即可。

該示意圖僅展示當GRACE月重力場模型為4個，貢獻源數量為4個，指紋數據量為5個(包含2個GIA指紋)的情況。對GIA貢獻的約束通過設計矩陣下方增加行來進行，所對應的數值根據公式(10)在mk-mk-1=mk+1-mk的假設下得到圖1 對GIA貢獻進行時間域約束Fig. 1 Scheme of implementing the time-domain constraint on GIA contributions

2 輸入數據

本節運用由美國空間研究中心(Center for Space Research，CSR)解算的RL06 GRACE月重力場模型(GRACE-Level2數據GSM)，包括從2003年1月至2016年8月共計149個月重力場模型。統一去除該時段內的平均重力場模型得到時變月重力場模型。由于地心運動主要受大規模質量移動影響，高階項系數表達的中小尺度質量遷移對其影響較小，可以忽略，因此采用解算至60階次的模型。另外，由于解算α所使用的系數達3 000多個，系數中的各種誤差可以得到很大緩解[10]，對月重力場模型無需采用任何去相關、濾波處理，從而可以避免由此所導致的信號畸變。

參考現有指紋數據集建立方案[12]，建立了一套新的共包含161個指紋數據的集合。其中，將南極地區按照Zwally等[14-15]劃分的匯水區分成27個區塊；格陵蘭島則首先按照匯水區劃分為8個區塊，再利用2 000 m高程線將每個區塊一分為二，形成16個區塊。假設區塊內各點的質量變化均一，區塊內總質量變化均為1 Gt，分別計算其指紋。陸地冰川一般為點狀分布，位置信息由世界冰山清單/空基全球陸地冰雪觀測(world glacier inventory/global land ice measurement from space，WGI/GLIMS)[16]數據集提供。將聚集在一起的冰山群劃歸為一個區塊，每個冰山群的總質量變化統一為1 Gt，但各點質量變化不同，與各點的冰山數量成正比。極地地區和陸地冰川的區塊劃分方案見圖2。由于TWS變化比較復雜，很難用區塊劃分的方式進行解構，因此利用Water GAP全球水文模型(water GAP global hydrology model，WGHM)[17-18]解析出TWS分量的前60個經驗正交分解函數(experiment orthogonal function，EOF)，針對每個EOF解算一個指紋。由于水文模型不能預測人為因素導致的質量變化，因此另外增加了印度地區指紋以便捕捉到該地區的地下水減少信號。注意，WGHM中與陸地冰川重疊的格網點需要事先去除。GIA指紋則是基于ICE-6G冰史解算的Fennoscandia 2個、Laurentide 3個以及南極地區1個共計6個區域GIA模型。

(a)南極冰蓋的劃分方式(27個區塊)；(b)格陵蘭島冰蓋的劃分方式(16個區塊)；(c)為陸地冰川集群，同時也展示了南極(10-20)和格陵蘭島(24-28)周邊的冰川集群圖2 極地冰蓋及陸地冰川區塊劃分

3 結果與分析

將指紋乘以相應的系數并累加(即TFα)應能恢復實測全球GRACE全球重力場變化。圖3(a)～(d)，展示了由原始GRACE數據解算得到的大地水準面變化，選取2003年1月、2006年6月、2009年9月和2015年12月作為示例進行展示(其他月份結果類似)。圖3(e)～3(h)為TFα得到的結果，可見指紋法反演結果與原始GRACE數據的結果特征一致，且無明顯條帶誤差。子圖3(i)～3(l)和3(m)～3(p)分別與子圖3(a)～3(d)和3(e)～3(h)類似，但展示內容為質量遷移。因此，基于該數據集所反演的地球重力場變化及地球系統質量遷移能夠捕捉到大尺度變化特征，可以用于反演地球重力場低階項。

(a)～(d)中分別展示了2003年1月、2006年6月、2009年9月和2015年12月由原始GRACE月重力場解算的大地水準面變化(未加入一階項和C20)；(e)～(h)中分別展示了上述4個月份基于指紋法反演的大地水準面變化(TFα)；(i)～(l)展示了上述月份由原始GRACE月重力場解算的地表質量遷移(未加入一階項和C20)；(m)～(p)中分別展示了上述月份基于指紋法反演的質量遷移

由于解算GRACE月重力場模型之前，大氣及海洋去噪1B級(atmosphere and ocean dealiasing level 1B，AOD1B)產品所模擬的高頻大氣和海洋動力信號已經扣除，因此GRACE月重力場模型中不包含相應信號，所得一階項結果中也不包含大氣和海洋動力信號的貢獻，可以直接用于補全GRACE月重力場模型。AOD1B產品的月平均值存儲于GRACE Level-2產品GAC中，因此只需加上GAC產品的一階項數據，即可得到地心運動總信號。圖4展示了由指紋法得到的地心運動時序及其誤差，并將其與基于GRACE-OBP法的結果進行比較。由于二者對應的GAC產品相同，即大氣與海洋動力信號的貢獻相同，圖中所展示結果均未加入GAC產品的一階項。

圖中展示由指紋法解算的一階項時序及其誤差，并與由GRACE-OBP法解算的一階項時序(條帶的寬度表示其誤差)進行對比。注意，所得結果未加入大氣與海洋動力信號的貢獻。為了易于辨識、在圖中對所示時間序列進行了上下平移，由上至下以此為C10，C11，S11

表1給出了由指紋法解算的地心運動與其他方法結果的年周期振幅、相位及趨勢信息。注意，為了與其他結果進行相比，表中基于指紋法和GRACE-OBP方法的結果均加上了大氣和海洋動力信號的貢獻。由表1可知，各種方法所解算的C10年周期振幅相差較大。Swenson方法解算的年周期振幅僅為1.9±0.1 mm，較其他方法小近1 mm。這是由于Swenson方法認為陸海質量交換時，流入、流出海洋的質量平均分布在整個海洋上。Sun等[1]利用模擬數據證明了該假設會嚴重低估地心運動C10(Z軸分量)的年周期振幅，考慮自引力效應對流入、流出海洋質量的分布影響之后，C10的年周期振幅明顯增大，達到2.7±0.1 mm，也即GRACE-OBP方法給出的結果。基于指紋法的地心運動時序的年周期振幅達3.0±0.1 mm 與GRACE-OBP方法和Rietbroek方法接近。而基于GPS和SLR數據的方法所得到的年周期振幅普遍較大，相較基于GRACE數據的方法，差值超過0.5 mm。由于C10主要反映的是地心在軸方向的運動，因此與質量在南北半球之間的運動較為相關，南北極地區質量變化觀測的精確程度成為關鍵因素，由于GPS和SLR在南北極地區的站點都相對稀疏，因此所得到的C10結果可能誤差較大。變換解算方式發現所得C10的年周期振幅變化幅度可達0.6 mm，也反映了該方法所得結果的不確定性。C11和S11的年周期振幅(2.0±0.1 mm)介于Swenson方法和GRACE-OBP方法之間。Sun等[1]的模擬并未發現Swenson方法在C11和S11的解算中會造成較大偏差，因此仍具有可比性。表1中列出的各種方法所得到的C11年周期振幅一致，相差不超過0.2 mm。各種方法所得S11的年周期振幅也較為一致，只有Rietbroek方法的結果偏大，原因不明。

表1 地心運動時序周期及趨勢統計

注：表中所統計的地心運動時序均包含大氣和海洋動力信號的貢獻

年周期相位方面，各種方法所得到的C11和S11結果相符，差別一般不超過2周。但是C10的年周期相位差別較大。基于GRACE數據的方法(包括GRACE-OBP[8，10]和指紋法)所得到的相位相較其他方法，如Wu等[20]使用的GPS反演法，差別較大(超過一個月)，其原因可能是GPS反演地面觀測數據分布不均勻。

線性趨勢方面，不同方法得到的結果差別較大。C10的趨勢較C11和S11明顯。指紋法所得到的趨勢為-0.66±0.22 mmyr-1，與其他兩種方法(Wu等[21]和Rietbroek等[19])的結果-0.88±0.09 mmyr-1和-1.08 mmyr-1符號一致但有一定差異。當單獨考慮線性趨勢的兩個分量(PDMT和GIA)時，差別也較大。例如，指紋法得到的GIA分量只有其他方法所得結果的一半，這可能與目前GIA模型較大的不確定性有關。

與其他方法相比，指紋法的優勢在于可以同時給出地心運動各個分量的貢獻，有利于解釋地心運動的驅動機制。圖5展示了各信號源對地心運動的貢獻。其中大氣和海洋(atmosphere and ocean，AO)由AOD1B產品提供。可以看到，地心運動的三個分量的周期性變化幾乎全部由TWS和AO引起。趨勢項則主要存在于C10中，主要由格陵蘭地區質量遷移信號所致。

由指紋法解算的地球系統各個分量(由上到下依次為ANT，GRE，GLA，TWS，GIA和AO)對地心運動的貢獻，其中AO分量由AOD1B產品提供。子圖(a)～(c)分別展示了C10, C11和S11的結果圖5 不同地心運動分量的貢獻

4 討論與結論

該研究表明指紋法可有效用于解算地心運動，也即地球重力場模型的一階項時變。其中C11和S11分量的年周期振幅、相位均與GPS反演法和GRACE-OPB法所得結果吻合。C10的年周期振幅為3.0±0.1 mm與GRACE-OBP接近，但與GPS反演法的差別為0.3～0.9 mm。考慮到GPS反演法所得結果自身具有較大的不確定性(可達0.6 mm)，可以認為結果合理。Sun等[1]利用端到端的數值模擬證明了在參數選擇合理的情況下GRACE-OBP方法可以準確解算地心運動；另外，Sun等[8]將各種方法所得一階項結果用于補全GRACE重力場模型，并解算沙漠和南極東部地區質量遷移信號微弱的區域，發現使用基于GRACE數據的一階項結果后可以有效消除研究區域質量遷移時序中的周期信號，這從兩方面印證了GRACE-OBP和指紋法的結果較為合理。

與其他方法相比，該方法可以獲取分布均勻且密集的數據，并同時給出地球系統各個分量對于地心運動的貢獻，有利于理解地心運動的驅動機制。本研究結果表明，地心運動的季節性周期運動主要由AO和TWS兩個分量驅動，且二者貢獻相當；趨勢項則主要由極地冰蓋消融驅動。值得注意的是，格陵蘭島和南極冰蓋的質量遷移對于地心運動的影響方向相反，但格陵蘭島作用明顯強于南極作用。GIA也會驅動地心運動的長期趨勢變化，本研究結果顯示GIA對于地心運動的貢獻很小，需要進一步研究確認。

本研究所使用的指紋法不僅是一種可行的地心運動計算方法，還具有其他應用潛力。例如，該方法所恢復的地球重力場沒有明顯條帶誤差的干擾，隨著人們對質量遷移的認識不斷加深，區塊劃分方案更加符合實際，利用指紋法恢復的月地球重力場模型將更加準確，其優勢在于可以同時補充一階項、改正C20和去除條帶誤差。另外，傳統方法利用GRACE計算某一區塊的質量變化時，由于其空間分辨率低，需要考慮信號泄露問題。而指紋法則可直接利用相應區塊已知的質量分布結合通過最小二乘方法解算得到的每個指紋(區塊)所對應的權重系數來計算其質量變化，避免了討論信號泄露問題。