朱傳東，劉金釗，王同慶，陳兆輝，張 品，張雙喜

中國地震局第一監(jiān)測中心，天津 300180

地球重力場是地球固有的物理特性之一，它反映了地球內(nèi)部物質(zhì)分布、運動及其變化的狀態(tài)。因此，準確確定地球重力場及其變化歷來是大地測量學、地球物理學、海洋學、空間科學以及地球動力學中的重要研究內(nèi)容，這對于深入研究地球形狀與內(nèi)部構造、防震減災、勘查礦產(chǎn)資源、保障國防安全等諸多領域具有非常重要的意義。

地球重力場時刻隨著地球質(zhì)量的重新再分布而變化。由于難以獲取瞬時的地球重力場模型，人們通常采用一段時期內(nèi)的觀測數(shù)據(jù)求取其平均地球重力場模型。如果觀測時間足夠短，則可將解算的模型看作反映該段時間重力場變化的“時變重力場模型”。重力恢復與氣候?qū)嶒炐l(wèi)星GRACE(gravity recovery and climate experiment)自2002年發(fā)射以來，已累積獲取了近15年的時變重力場模型序列，這些數(shù)據(jù)一般按月或周發(fā)布，并以離散的引力位球諧系數(shù)(Clm,Slm)表達，以代表觀測期間全球重力場的平均變化[1]。但事實上，地球重力場是隨時間連續(xù)變化的，為更好描述其動態(tài)變化過程，就需為離散的球諧系數(shù)構建函數(shù)模型，即在傳統(tǒng)的描述三維空間變化的重力場模型基礎上引入時間變量，以構造出能同時描述重力場時空變化的四維結構模型(圖1)，這能為各種時變重力場參量特征分析及其應用研究提供有意義的基礎參考。

圖1 GRACE時變重力場模型Fig.1 GRACE time variable gravity model

現(xiàn)有的GRACE時變重力場在其動態(tài)變化模型的提取及構建過程中仍面臨兩個主要問題。一方面，受儀器測量誤差、相關誤差、截斷誤差和混頻誤差等影響，GRACE解算結果的中高階球諧系數(shù)會存在較大誤差。目前，聯(lián)合以高斯濾波為代表的濾波因子和去相關濾波方法雖然可有效去除部分誤差，但不同濾波方法的參數(shù)選取不同，則時變重力場模型的中短波部分的改正量也會有所不同，其相應的濾波效果(信噪比)也會有所差異[2-6]。另一方面，在年際、季節(jié)性或更短的時間尺度上，時變重力場變化主要產(chǎn)生于地球表面圈層的水循環(huán)、大氣及海洋質(zhì)量變化等[7]。因此，多種時空尺度地球物理場過程的相互復雜作用，本質(zhì)上決定了時變重力場的非線性變化特征。例如，以往的研究結果表明，基于不同時間跨度的GRACE數(shù)據(jù)反演的極地冰蓋物質(zhì)平衡在整體和局部上呈現(xiàn)出顯著的年際變化特征，可導致不同時間段的冰雪質(zhì)量變化趨勢存在顯著差異[8]。即使在全球范圍內(nèi)，GRACE反演的陸地水質(zhì)量變化也存在顯著的加速度信號，這會對趨勢信號的評估造成較大干擾[9-10]。在此背景下，僅顧及趨勢項和季節(jié)項的影響且以時間序列必須線性演變或諧波型周期變化為前提的假設方法，就很難準確描述時變重力場的非線性動態(tài)變化過程[11-12]。

目前，以經(jīng)驗正交函數(shù)分解(EOF)及其擴展形式多通道奇異譜分析(MSSA)為代表的統(tǒng)計分析方法，是構建GRACE離散時變重力場動態(tài)變化模型的有效技術手段[13-14]。其中，EOF又稱為主成分分析，在以往的研究中，可在頻譜域和空間域上應用于GRACE矩陣數(shù)據(jù)的數(shù)據(jù)壓縮降維、噪聲消除[15-16]，以及主要地球物理場模態(tài)的提取分析解釋[17-19]。不過，EOF在將GRACE矩陣數(shù)據(jù)變換為統(tǒng)計不相關的新變量過程中會存在模態(tài)混疊問題，由此不能準確解釋各分量的真實物理意義[14,20]。有別于EOF方法，MSSA可進一步利用遲后協(xié)方差矩陣，通過遲后排列的空間型表達原場中主要空間結構的時間演變特征。對于事先未知物理本質(zhì)的系統(tǒng)，MSSA尤其不需要假定時間序列必須線性演變或諧波型周期變化，而是根據(jù)資料本身最優(yōu)確定的，從而識別原序列的傾向、周期和噪聲信號[21]。有鑒于此，MSSA特別適用于提取GRACE矩陣數(shù)據(jù)的非線性變化特征[11,22]，尤其是在局部空間域上分析、解釋與多種氣候模態(tài)相關的GRACE陸地水文變化過程[23-24]。不過，MSSA直接應用于GRACE月球諧系數(shù)矩陣的信號提取過程中會受到系統(tǒng)相關誤差和計算效率的顯著影響，由此不能有效識別其準確的物理變化特征[25]。針對GRACE月時變重力場模型的誤差特征和演變特性，本文結合EOF和MSSA兩種方法的互補性，在利用EOF對球諧系數(shù)矩陣實現(xiàn)降噪和去冗余的基礎上，進一步利用MSSA提取隱藏在EOF混疊模態(tài)中的非線性傾向和周期信號特征，以期有效構建可以準確描述GRACE離散時變重力場的大尺度動態(tài)變化模型。

1 數(shù)據(jù)處理

本文采用的GRACE數(shù)據(jù)為CSR(center for space research)提供的2002年8月至2016年8月共153個月的RL05版本時變重力場模型，球諧系數(shù)截斷到60階，其中潮汐、大氣和海水質(zhì)量變化的影響均已采用模型扣除。對GRACE月時變重力場模型作如下處理：采用衛(wèi)星激光測距測得的更為準確的C20項對其作替換[26]；采用GRACE和海洋模型數(shù)據(jù)解算的一階項加回到球諧系數(shù)中[27]；對球諧系數(shù)時間序列中16個月的缺測值采用最小二乘擬合作插補；從模型中扣除時間段內(nèi)球諧系數(shù)的平均值；為削弱模型中高階項球諧系數(shù)噪聲的影響，采用300 km高斯濾波對其做平滑處理[2]；經(jīng)過高斯濾波處理后，某些固定次的球諧系數(shù)的偶(奇)數(shù)階之間仍存在明顯的系統(tǒng)相關誤差，這會導致GRACE反演結果在空間上呈現(xiàn)出顯著的南北條帶噪聲。為此，對以上處理得到的169個月的時變重力場模型構成的去中心化的球諧系數(shù)矩陣，見式(1)

(1)

采用EOF分別將其分解成相互正交的空間特征向量和空間主成分的乘積之和，見式(2)

(2)

通過顯著性檢驗選取其前幾個空間主成分即可有效反映時變重力場模型的主要信號變化特征。這可起到兩方面的作用：①有效消除球諧系數(shù)中高階項系數(shù)存在的顯著系統(tǒng)相關誤差，提取其更優(yōu)信噪比的信號[15-16]；②有效降低球諧系數(shù)矩陣空間維數(shù)的同時可顯著減少后續(xù)MSSA的計算量。

對于式(2)中通過顯著性檢驗的S-PC構成的多維時間序列xlt，其中l(wèi)=1,2,…,L是通道序號，t=1,2,…,T是時間序號，包含了多種模態(tài)的混疊信號。為此，進一步利用MSSA提取隱藏在xlt中的非線性傾向和周期信號[28]。首先，把去中心化的xlt在時間遲后相空間排列，構建得到軌跡矩陣X，其第i+1列，也就是第i+1個狀態(tài)向量為

Xi=(x1i+1,…,x1i+M,…,xLi+1,…,xLi+M)

(3)

式中，i=0,1,…,N-M；后遲量M稱為窗口長度。

然后，對矩陣X作EOF分解，其分量形式為

(4)

(5)

MSSA最重要的應用功能是通過重建成分RC(reconstructed component)實現(xiàn)的，并由此從被分析對象中提取出與某個或某些特征成分相聯(lián)系的部分[29]。根據(jù)最小二乘原理，由第k個特征成分得到的重建成分為

(6)

(7)

2 計算結果

2.1 EOF分解

對球諧系數(shù)矩陣ZC和ZS分別作EOF分解，以對系數(shù)矩陣進行降維和除噪。圖2給出了前50個S-EOF的方差貢獻，可以看出方差較大的值主要集中在前面幾個且遞減速度很快。依次對S-PC進行柯爾莫諾夫-斯米爾諾夫檢驗，結果表明在95%的置信水平下，ZC和ZS前7個特征值中除第6個以外均通過了顯著性檢驗，相應的累積方差分別達到92.8%和91.0%。

圖2 前50個S-PC的方差貢獻Fig.2 The eigen value spectra of the first 50 S-PC

選取通過顯著性檢驗的S-PC和S-EOF重建統(tǒng)計量ZC和ZS，并由此反演全球陸表網(wǎng)格的等效水柱高變化時間序列。為了說明EOF的除噪效果，圖3比較分析了經(jīng)高斯濾波和“高斯濾波&EOF”處理后GRACE反演的2008年6月的全球陸表水儲量變化結果?？梢钥闯?，只經(jīng)過高斯濾波處理后的反演結果在兩極地區(qū)的條帶噪聲盡管有所減弱，但在赤道兩側地區(qū)依舊比較明顯；而經(jīng)過“高斯濾波&EOF”處理后的條帶噪聲明顯減少，且局部地區(qū)的信號變化得到有效呈現(xiàn)。另外，本文進一步比較分析了經(jīng)“高斯濾波&EOF”和經(jīng)驗性的“高斯濾波&P5M11”[3-4]處理后的GRACE反演結果。圖4對應給出了兩種反演結果與GLDAS Noah水文模式(經(jīng)高斯濾波處理)全球陸表之差的標準差序列。整體上看，“高斯濾波&EOF”要小于“高斯濾波& P5M11”處理后的標準差結果，相應的平均標準差分別為6.1和6.3 cm，從以上比較可以看出，EOF分解可有效消除GRACE球諧系數(shù)中系統(tǒng)相關誤差的影響，提取其更優(yōu)信噪比的信號。

圖3 GRACE 反演的2008年6月全球陸表水儲量變化Fig.3 The terrestrial water change recovered from GRACE in June 2008

圖4 GRACE反演結果與GLDAS Noah之差的標準差Fig.4 The standard deviation of the differences between GRACE and GLDAS Noah

2.2 季節(jié)性變化

圖5 前50個ST-PC的方差貢獻Fig.5 The eigenvalue spectra of the first 50 ST-PC

圖6給出了第2個和第3個特征值對應的ST-EOF和ST-PC。據(jù)圖5，特征值2與3接近相等，同時滿足|λk-λk+1|≤min{σλk,σλk+1}，且兩個特征值對應的ST-EOF和ST-PC呈現(xiàn)出顯著的周期性變化并相互正交。根據(jù)文獻[29]，ST-EOF2&3可構成一對潛在的周期振蕩模態(tài)，兩者的總方差貢獻可達到50.6%。

圖6 第2個和第3個周期振蕩模態(tài)ST-EOF及其對應的ST-PCFig.6 The second and third ST-EOF and ST-PC

采用ST-EOF2&3和ST-PC2&3重建MSSA分析對象RC，即對S-PC序列的擬合。選取RC方差最大的通道序列做最大熵譜估計，結果表明ST-EOF2&3為周期12個月的周年振蕩模態(tài)。進一步由RC計算GRACE球諧系數(shù)的周年變化成分STRC(spatial-temporal RC)，并由此反演全球陸表網(wǎng)格的等效水柱高變化時間序列。圖7給出了全球陸表水各網(wǎng)格的周年振蕩最高相位與最低相位對應的差值結果?？梢钥闯?，陸表水周年振蕩變化較大的區(qū)域主要位于赤道兩側，其中亞馬孫流域最大可達99.7 cm。

為進一步分析陸表水的周年振蕩變化特征，選取RC方差最大的通道序列，采用位相劃分討論方法，取各個位相陸表水各網(wǎng)格點的平均值得到周年振蕩第1-8位相合成圖(圖8)。由圖8可知，前后4個位相的正負距平分布基本相反。其中，第2位相和第3位相在中緯度地區(qū)呈現(xiàn)較大正值，在赤道地區(qū)呈現(xiàn)較大負值，而第6位相和第7位相則與之相反，反映出南北半球陸表水的季節(jié)性交替變化特征。

2.3 年際性變化

對矩陣ZC和ZS，通過顯著性檢驗的S-PC時間序列做MSSA分析之前，為避免周年項成分對年際信號提取的干擾，首先將其從S-PC中作扣除。圖9給出了前50個ST-EOF的方差貢獻。可以看出，在誤差范圍內(nèi)前50個特征值有2個達到95%的置信水平，相應的累計方差達到77.4%。

圖10給出了第1個和第2個模態(tài)ST-EOF及其對應的ST-PC。可以看出，前兩個ST-PC表現(xiàn)出一定的傾向特征，其方差貢獻分別達到64.5%和12.9%。由于資料相對較短，目前還很難對兩者長期變化的演變趨勢作出判定。

采用ST-EOF1&2和ST-PC1&2分別重建GRACE球諧系數(shù)的長期趨勢變化成分，并據(jù)此反演全球陸表網(wǎng)格的等效水柱高變化時間序列。圖11對應給出了兩者由一元線性回歸估算的全球陸表水變化趨勢?？梢钥闯觯芙邓?、地下水損耗和冰川消融等的影響，本文給出的兩種趨勢信號的空間分布與以往的研究結果基本類似[30-33]。另外，受冰川均衡調(diào)整的混疊影響，北美哈德遜灣、斯堪的納維亞半島和西南極等地區(qū)呈現(xiàn)出較為顯著的質(zhì)量增加信號[34-35]。

2.4 動態(tài)變化模型

基于EOF&MSSA重建的“周年項分量+趨勢項分量”構建得到GRACE球諧系數(shù)(Clm,Slm)的非線性模型STRCnonlinear。為從頻譜域上說明STRCnonlinear的擬合效果，以球諧系數(shù)(C10,1,S10,1)為例，圖12比較分析了其變化時間序列與非線性模型。結果表明，本文構建的模型能較好擬合(C10,1,S10,1)時間序列的非線性變化特征，擬合標準差分別為1.9×10-12和2.0×10-12。圖13進一步給出了EOF&MSSA得到的GRACE 60階次球諧系數(shù)的擬合標準差分布。整體上該方法得到擬合標準差由低階到高階逐步遞減，結果介于2.5×10-15～5.0×10-11，平均值為4.3×10-13。

圖7 全球陸表水周年振蕩最高與最低相位對應的差值Fig.7 The differences between extreme high and low phases of the annual oscillation

圖8 陸地水周年振蕩第1—8位相合成圖Fig.8 Spatial structures of the terrestrial water annual oscillation keyed to eight phases

圖9 前50個ST-PC的方差貢獻Fig.9 The eigen value spectra of the first 50 ST-PC

為從空間域上說明STRCnonlinear的擬合效果，本文進一步比較分析了EOF&MSSA和最小二乘擬合法(考慮趨勢項和周年項，對只經(jīng)過高斯濾波和EOF處理后的60階次的球諧系數(shù)時間序列逐個做擬合計算) 反演得到的全球等效水柱高變化結果。以陸地水變化較為顯著的亞馬孫流域(2.5°S,58.5°W)和南極半島(75.5°S,98.5°W)為例(圖14)，EOF&MSSA和最小二乘擬合法得到的兩個網(wǎng)格點上等效水柱高的擬合標準差分別為13.6和15.1、9.6和10.1 cm，相比而言，EOF&MSSA構建的非線性模型能更好地反映兩個區(qū)域等效水質(zhì)量的動態(tài)變化特征。圖15進一步給出了全球陸表共21 334個網(wǎng)格點上的擬合統(tǒng)計結果，結果表明，EOF&MSSA得到的擬合標準差介于0.3～14.0 cm，平均標準差為1.8 cm，而最小二乘擬合法的結果介于0.3～15.1 cm，平均標準差為1.9 cm，且在69%的網(wǎng)格點上EOF&MSSA得到的擬合標準差要小于最小二乘擬合法的結果，整體上可以看出EOF&MSSA構建的非線性模型能更好地擬合GRACE離散球諧系數(shù)的動態(tài)變化特征。需要指出的是，以上結果分別是由EOF&MSSA對12個S-PC構成的多維時間序列和最小二乘擬合法對3720個球諧系數(shù)時間序列擬合得到的結果。因此，與最小二乘擬合法相比，EOF&MSSA在顯著減少系統(tǒng)相關誤差和計算量的同時即可取得更好的擬合精度，反映出其在提取和構建GARCE球諧系數(shù)大尺度動態(tài)變化模型方面的優(yōu)越性。

圖10 第1個和第2個模態(tài)ST-EOF及其對應的ST-PCFig.10 The first and second ST-EOF and ST-PC

圖11 全球陸表水等效水柱高變化趨勢Fig.11 The terrestrial water change trend

圖12 GRACE球諧系數(shù)及其非線性模型變化時間序列Fig.12 GRACE SHCs time series and its nonlinear model

圖13 GRACE 60階次球諧系數(shù)與其非線性模型之差的標準差Fig.13 The standard deviation of the differences between GRACE SHCs and its nonlinear model

圖14 EOF&MSSA法和最小二乘擬合法反演得到的等效水柱高時間序列Fig.14 The time series of equivalent water heights from EOF&MSSA and least square method

圖15 全球陸表網(wǎng)格點等效水柱高的擬合標準差Fig.15 The fitting standard deviation of equivalent water heights on grids over global land areas

另外，從圖13和圖15全球尺度上的結果來看，本文采用EOF&MSSA構建的非線性模型對GRACE球諧系數(shù)及其反演結果仍存在一定的擬合差異。從局部尺度上看，這種擬合差異可能主要來源于兩個方面：①受多種誤差因素的影響，GRACE高階項球諧系數(shù)觀測精度相對較差，這較大干擾了直接從球諧系數(shù)提取ENSO、北極濤動/北大西洋濤動等氣候模態(tài)引起的區(qū)域陸地水文信號變化[23,36-37]；②受某些區(qū)域極端干旱、洪澇等氣候因素的影響，GRACE球諧系數(shù)會存在一些異常信號，這會進一步干擾對球諧系數(shù)變化的擬合精度[38-39]。由此可以看出，目前GRACE重力衛(wèi)星獲取的離散時變重力場模型的精度未來仍有提升空間。

3 結 論

針對GRACE月時變重力場的誤差特征和演變特性，提出聯(lián)合經(jīng)驗正交分解和多通道奇異譜分析的方法，提取、分析其動態(tài)變化特征并構建了相應的非線性變化模型。相比經(jīng)驗性的去相關濾波方法，EOF對GRACE月球諧系數(shù)實現(xiàn)空間壓縮的同時可有效去除系統(tǒng)相關誤差的影響，提取其更優(yōu)信噪比的信號。多通道奇異譜分析的提取結果表明，GRACE離散球諧系數(shù)存在顯著的長期趨勢項和周年項變化成分，其中，周年項模態(tài)的方差貢獻達到50.6%，扣除周年項成分后長期趨勢項的累積方差貢獻達到77.4%。以EOF&MSSA提取的長期趨勢項和周年項成分構建了GRACE離散球諧系數(shù)的非線性模型，其擬合標準差介于2.5×10-15～5.0×10-11，平均標準差為4.3×10-13，其對全球陸表等效水柱高的擬合標準差介于0.3～14.0 cm，平均標準差為1.8 cm。相比傳統(tǒng)的最小二乘擬合方法，EOF&MSSA在顯著減少系統(tǒng)相關誤差和計算量的前提下即可在頻譜域和空間域上更好地擬合GRACE離散時變重力場模型的大尺度動態(tài)變化特征。

從本文的研究結果看，目前GRACE重力衛(wèi)星獲取的時變重力場的精度未來仍有提升的空間。隨著下一代重力衛(wèi)星的GRACE FO發(fā)射，借助其更高時空分辨率的球諧系數(shù)數(shù)據(jù)，有望為離散時變重力場各種參量的動態(tài)變化模型構建及其在全球氣候變化、地球動力學等中的應用提供重要參考。