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

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

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

(1)

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

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

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

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

1 算法介紹

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

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

(2)

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

(3)

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

(4)

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

(5)

其中T為截斷矩陣：

(6)

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

(7)

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

(8)

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

(9)

(10)

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

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

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

2 輸入數(shù)據(jù)

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

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

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

3 結(jié)果與分析

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

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

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

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

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

表1 地心運(yùn)動時序周期及趨勢統(tǒng)計

注：表中所統(tǒng)計的地心運(yùn)動時序均包含大氣和海洋動力信號的貢獻(xiàn)

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

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

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

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

4 討論與結(jié)論

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

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

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