反演地表質量變化的附有方差約束的徑向點質量方法

尹恒， 游為， 范東明， 方偉浩， 萬祥禹， 宋夢芝

西南交通大學地球科學與環境工程學院， 成都 611756

0 引言

2002年以來，由美德聯合實施的GRACE重力衛星以前所未有的精度和時空分辨率觀測時變地球重力場信號，提供了研究地球表層質量變化的直接觀測手段(Wahr et al., 2004; Tapley et al., 2004b).近20年來，基于GRACE時變重力場進行的全球或局部陸地水儲量變化研究(Rodell et al., 2007；馮偉, 2013；Zhong et al., 2018)、極地冰蓋和高山冰川質量變化研究(Luthcke et al., 2006b, 2013；Baur, 2013；朱傳東等, 2015；高春春等, 2016)、海平面變化研究(Chambers, 2006；Uebbing et al., 2019)等取得了豐碩成果，促進了大地測量學、地球物理學、海洋學和冰川學等領域的發展，因此利用重力衛星觀測數據研究全球質量變化具有重大意義.

目前利用時變地球重力場表達地表物質遷移的方法主要分為球諧位系數法(Wahr et al., 1998)、Mascon法(Rowlands et al., 2005；Luthcke et al., 2006a)、點質量法(Baur and Sneeuw, 2011)和徑向基函數(楊帆, 2017)等方法.由于GRACE衛星本身的儀器載荷誤差、近極近圓飛行軌道和高頻背景攝動力模型誤差等問題，使得高階次球諧位系數被噪聲嚴重污染，導致直接使用球諧位系數反演的地表質量變化呈現嚴重的南北條帶誤差(Tapley et al., 2004a).許多學者研究了眾多濾波平滑技術，如常用的多項式擬合去相關算法和高斯濾波方法(Wahr et al., 1998；Swenson and Wahr, 2006).雖然濾波方法的實現使得利用球諧位系數法計算地表質量變化被廣泛應用，但濾波在過濾噪聲的同時不可避免地會過濾掉有用的信號甚至改變信號的空間分布，相比于對球諧位系數施加平滑濾波，對地球物理信號施加時間和空間約束能得到更高精度和空間分辨率的地表質量變化反演結果(Luthcke et al., 2006a)，因此有學者提出了另一種計算地表質量遷移的Mascon方法.

Mascon方法最早被用于月球重力場模型的確定(Muller and Sjogren, 1968)，其基本思路是將研究區域劃分為若干質量塊并假設每個質量塊內部質量變化相等，直接以每個質量塊的質量變化作為未知參數，建立衛星軌道和星間距離變率等原始觀測值與未知參數的函數關系.Rowlands等(2005)、Luthcke等(2006a)和Rodell等(2007)將地表質量塊與GRACE衛星K波段星間距離變率建立聯系，利用時空約束方程求解地表質量變化，通過與水文模型和衛星測高對比，驗證了該方法的有效性.Watkins等(2015)利用星間距離變率作為觀測值，引入GLDAS (Global Land Data Assimilation System)、衛星測高和ECCO2(Estimating the Circulation and Climate of the Ocean Phase 2)海洋模型等外部物理信號建立時空約束方程反演等效水高變化，Save等(2016)利用球諧系數法反演結果作為先驗信息構建約束矩陣反演全球質量變化.

有學者發現也可直接建立地表質量塊與重力衛星所受攝動加速度(或重力位)的函數關系來求解地表質量變化，該類方法被稱為點質量方法.目前點質量方法主要有擾動位點質量方法(Han et al., 2005)、徑向點質量方法(Baur and Sneeuw, 2011)以及三維加速度點質量方法(蘇勇等, 2017)，其觀測值分別為擾動位變化、徑向加速度變化以及三維加速度變化.Han等(2005)使用擾動位點質量方法反演了亞馬遜流域和奧里諾科河流域的質量變化，表明其空間分辨率優于球諧位系數法結果.Baur和Sneeuw(2011)、Baur(2013)利用徑向點質量方法計算了格陵蘭島冰蓋質量變化，驗證了該方法的有效性.郭飛霄等(2018)、蘇勇等(2019)和魏偉等(2020)分別利用附有空間約束的徑向點質量方法和三維點質量方法反演了亞馬遜流域和華北平原陸地水儲量變化，對比分析了不同空間約束的效果，表明采用的各類空間約束均能提高反演結果精度，其中線性形式的空間約束效果最好.Chen等(2016,2020)利用徑向點質量方法通過截斷奇異值分解和Tikhonov正則化方法計算青藏高原冰川凍土消融情況，提高了計算結果的信噪比.Forsberg等(2017)利用徑向點質量方法計算了格陵蘭島和南極大陸冰蓋融化情況對海平面上升的貢獻.Ran等(2018)以地表質量變化先驗誤差信息為約束條件，利用基于特征值分解的徑向點質量方法計算了格陵蘭島冰蓋消融情況.大量研究結果已驗證了點質量方法的可行性和有效性，但由于約束矩陣的缺陷，其反演陸地水儲量變化結果仍受條帶誤差和海陸信號泄漏的影響.

本文仿照CSR Mascon正則化矩陣的構造思路，幾乎不加入任何外部數據的質量變化信息，直接以GRACE時變重力場模型的計算結果作為先驗質量變化信息，對CSR Mascon正則化矩陣的構造思路進行若干改進，推導出了附有方差約束的徑向點質量方法.具體改進的地方包括：(1)以單位矩陣為正則化矩陣的徑向點質量方法計算的質量變化的RMS作為先驗質量變化信息；(2)考慮各種質量變化信號的正則化參數和正則化矩陣不唯一，構造了迭代正則化的處理策略.

1 數據與方法

1.1 數據處理與驗證數據

本文選用CSR發布的2003-01—2014-11共133個月的Level-2 Release 06 GSM時變重力場模型數據反演地表質量變化, 并將球諧位系數截斷至60階.由于GRACE時變重力場解算的一階項為0(C10,C11,S10)，采用Sun等(2016)解算的一階項替換原有一階項(Swenson et al., 2008； Sun et al., 2016).由于GRACE解算的C20項精度較低，采用SLR (Satellite Laser Ranging,衛星激光測距)解算的C20替代GRACE球諧位系數中的C20項(Cheng et al., 2013； Loomis et al., 2020).

采用最新的GLDAS_NOAH_M.2.1水文模型(簡稱GLDAS2.1)用于評估本文反演的陸地水儲量變化精度.該模型是一組地表水和能量場的數據集，由GSFC(Goddard Space Flight Center,戈達德航天中心)和NCEP(National Centers for Environmental Prediction,美國國家環境預報中心)聯合研制(Rodell et al., 2004)，本文采用其提供的空間分辨率為0.25°×0.25°，時間分辨率為1個月的模式數據.

為更好地評價反演結果，本文采用CSR和JPL RL06 Mascon(以下簡稱CSRM和JPLM)模型進行全球、陸地、冰凍圈信號綜合對比.由于上述Mascon模型均扣除了2004年1月至2009年12月的均值，為與MASCON模型保持一致，本文所用球諧位系數也扣除相同時段的球諧位系數均值(Watkins et al., 2015；Save et al., 2016).

1.2 空間格網

徑向點質量方法通過將地表點質量變化與空間點重力擾動建立函數關系求解區域或全球質量變化，因此需要事先將地球表面劃分為不同質量塊，并確保所有質量塊均勻地覆蓋在地球表面.為確保劃分地表質量塊面積的均勻性，本文采用基于二十面體劃分的等面積格網，格網劃分方法可參考Majewski等(2002).為了有效減少海陸泄漏誤差，本文按上述方法將地球表面劃分為40962個等面積格網，每個格網面積約為12450 km2，格網點之間距離約為120 km，即格網點空間分辨率約為1°.圖1為本文用于徑向點質量方法的等面積海陸格網空間分布圖，其中紫色表示陸地格網，藍色表示海洋格網(圖中白色部分為地圖投影問題，真實格網已完全覆蓋地球表面).

1.3 徑向點質量模型方法基本原理

地球表層質量變化引起的重力擾動徑向分量可由球諧位系數變化量 ΔClm,ΔSlm表示為(Baur and Sneeuw, 2011)：

×[ΔClmcos(mλ)+ΔSl msin(mλ)],(1)

圖2 徑向點質量方法幾何示意圖Fig.2 Radial point-mass modeling geometry

如圖2所示，在地心坐標系中有質量變化為 δmj的地面點Pj(λj,φj,a)，根據牛頓第二運動定律，其引起的外部空間點Si(λi,φi,r)的徑向重力擾動δgij可表示為：

(2)

cosψij=sinφisinφj+cosφicosφjcos(λi-λj),(3)

式中，lij是空間點Si到地面點Pj的距離，ψij為空間點Si到地面點Pj的球面角距.根據圖2中的幾何關系，特定區域內所有地表點質量變化引起的某空間點徑向重力擾動可表示為：

(4)

式(4)左端為時變重力場球諧位系數計算的徑向重力擾動偽觀測值，式(4)右端δmj為待求的單個地表質量變化未知數.當利用等效水高變化來表示地表質量變化時，可設δmj=ρw×sj×Δhj，其中ρw=1025 kg·m-3為水的密度，sj為單個格網面積，Δhj為等效水高變化，將單位換算為cm后，式(4)表示為：

i=1,…,s.

(5)

聯合式(1)和式(5)即可得利用徑向點質量方法計算地球表面等效水高變化的函數模型，其中i=1,…,s為空間點的個數，j=1,…,p為研究區域地面格網點的個數，本文令空間點個數s和地面格網點個數p均為40962.進一步，將式(5)線性化可得線性方程：

y=Ax+e,(6)

式中y(s×1)表示觀測值向量，A(s×j)表示設計矩陣，x(j×1)為待求未知數向量，e(s×1)為觀測值擬合殘差.

1.4 參數估計方法

由于利用衛星高度處的重力擾動觀測值計算地表質量變化本質上是一個重力向下延拓問題，會使設計矩陣A呈現病態性，因此直接用最小二乘法則求解式(6)不能獲得穩定解(Kusche and Klees, 2002).為了得到穩定的估值，本文利用Tikhonov正則化原理，引入正則化參數和正則化矩陣作為約束條件，構造以下估計準則(王振杰, 2003)：

min{‖y-Ax‖2+α‖Mx‖2},(7)

將式(7)對x求偏導數并聯合式(6)可得待估參數的解為：

(8)

式中α是正則化參數，用于平衡觀測值殘差和約束條件；M為正則化矩陣，為方程提供約束條件，‖·‖是歐式2-范數.

與最小二乘解相比, 式(8)右端求逆部分加入了αMTM一項，由于它的作用抑制了法方程的病態性, 使得(ATA+αMTM)的求逆變得穩定, 因此式(8)可以得到穩定的估值.式(8)的結果與正則化參數α有關,過小的正則化參數會導致解仍然不穩定，過大的參數會使反演結果過度平滑.常見選取正則化參數的方法有L-曲線法(Hansen, 1992)、U-曲線法(Krawczyk-Stańdo and Rudnicki,2007)、GCV(Generalized Cross-Validation，廣義交叉驗證)方法(Xu, 2009)和RMSE最小法(Xu et al., 2006)等.本文使用Hansen(1992)針對線性病態問題提出的L曲線法估計正則化參數，以‖y-Ax‖為橫坐標，‖x‖為縱坐標獲得L曲線圖，L曲線上曲率最大的點對應參數即為L曲線法確定的正則化參數，但需要注意的是此方法確定的正則化參數僅是近似最優的.

作為約束條件的正則化矩陣M有多種構造方法，最初的方法是直接將M設為單位矩陣(Baur and Sneeuw, 2011; 郭飛霄等, 2017)，后來有學者利用空間平滑條件設計正則化矩陣，以此對地表質量變化點進行空間平滑約束(郭飛霄等, 2018; 魏偉等, 2020).直接利用單位矩陣或空間平滑條件的約束矩陣均能有效去除部分條帶誤差，但其結果仍受部分條帶誤差影響(蘇勇等, 2019).

1.5 無空間約束徑向點質量方法反演地表質量變化

本節利用CSR RL06 GSM數據通過正則化矩陣為單位矩陣的無空間約束徑向點質量方法(以下簡稱PM)計算2003-01—2014-11共133個月的全球等效水高，圖3為PM反演各月份等效水高的L曲線圖，不同月份對應的L曲線有所差異，取所有L曲線均值作為平均正則化曲線.

圖3 2003-01—2014-11不同月份L曲線(綠)和平均L曲線(紅)Fig.3 L curve (green) and average L curve (red) in different months from January 2003 to November 2014

圖3中的紅色曲線為所有月份平均L曲線，當正則化參數為1×10-19時，L曲線曲率最大，因此將其作為所有月份的正則化參數.為了比較不同正則化參數的效果，同時選取了大于和小于最佳正則化參數反演結果作為對比，繪制了如圖4的2010-01等效水高分布圖，為了更直觀地對比反演效果，同時繪制了濾波半徑為300 km的高斯濾波反演結果分布圖.

對比圖4b和4d可以看出，當正則化參數為1×10-19(最優正則化參數)時，PM反演結果仍受條帶誤差和泄漏誤差影響，其條帶誤差與濾波半徑為300 km的高斯濾波反演結果相當.當增大正則化參數至3×10-19時(圖4c)，PM反演結果信號強度隨之減小，說明過大的正則化參數會抑制反演結果的信號強度，但與此同時條帶誤差也在減小.當減小正則化參數為5×10-20時(圖4a)，信號強度和條帶誤差都顯著增大.從以上分析可以看出PM仍受條帶誤差和泄漏誤差影響，通過增大正則化參數可減少條帶誤差，但也會嚴重削弱信號.

2 附有方差約束的徑向點質量方法

針對無空間約束的PM仍受條帶誤差和海陸信號泄漏影響的缺陷，本文利用無空間約束的PM反演結果設計正則化矩陣，構建附有方差約束的徑向點質量方法，旨在通過先驗信息的約束，消除條帶誤差和海陸泄漏誤差.

2.1 正則化矩陣

VCPM需要借助質量塊先驗方差信息設計正則化矩陣，Watkins等(2015)利用陸地水文模型、衛星測高和海底壓力模型等多源數據的方差信息構造約束矩陣，Save等(2016)利用GRACE球諧位系數法反演結果估計質量塊的方差信息構建正則化矩陣.式(9)為本文所用的正則化矩陣結構，由于Save等(2016)指出質量塊之間的協方差信息會給反演結果帶來橫向誤差，因此僅取質量塊方差信息設計正則化矩陣.

(9)

圖4 不同正則化參數的點質量方法(a,b,c)和高斯濾波半徑為300 km的球諧位系數法(d)反演的2010-01等效水高分布圖Fig.4 Equivalent water height estimated by (a,b,c) point-mass method with different regularization parameters or (d) spherical harmonic potential coefficient method with Gaussian filter radius of 300 km in January 2010

式中σi,i=1,2,…,n為第i個質量塊的方差信息，通過計算第i個質量塊多個月的RMS值得到.Save等(2016)設計了兩種正則化矩陣，一種是直接以2004—2009年等效水高RMS設計的0409F，另一種是在0409F基礎上再結合其他單個月的等效水高信息設計的附有時間相關性的0409V.結果表明：由于0409V包含了2004—2009年以外的陸地水儲量變化信息，因此以0409V反演的2004—2009年時間段以外的等效水高解效果優于0409F反演結果，而在2004—2009年時間段以內，兩個正則化矩陣反演結果相同.本文直接以2003—2014年全時段PM反演結果構建正則化矩陣，包含了時間段內的所有地表水儲量變化信號特征，因此不再構建時間相關的正則化矩陣.

本文的正則化矩陣分兩步設計，第一步：設計初始正則化矩陣，獲得中間解，這一步主要是為了消除海陸泄漏誤差，同時獲得海洋區域質量塊真實方差信息；第二步：利用第一步獲得的真實海洋方差信息和初始陸地質量變化信息構建迭代正則化矩陣，通過迭代的方法逐步提取質量變化信號.對于第一步的初始正則化矩陣，本文與Save等(2016)的方法完全相同，第二步中Save等(2016)采用扣除了陸地振幅信號和高緯冰川趨勢信號的剩余信號構建正則化矩陣，利用該矩陣一次性反演全球等效水高信號，本文為充分提取觀測值中的物理信號設計了三個不同的正則化矩陣迭代地反演全球等效水高.

2.2 初始正則化矩陣設計

Save等(2016)利用Save等(2012)正則化的時變重力場模型經200 km高斯濾波反演的等效水高結果構建初始正則化矩陣，本文使用信噪比更高的PM計算結果構建初始正則化矩陣.圖5a為利用PM反演的133個月等效水高計算的均方根值(RMS)分布圖，從中可以直觀地看到許多代表質量變化的物理信號：如陸地冰川冰蓋融化信號，各大流域季節性水儲量變化信號，印度西北部、華北平原、美國加利福利亞地區地下水消減信號等；除了陸地信號之外，海洋區域也能看出部分質量變化較強的信號，其中蘇門答臘、東京和智利等海底地震信號最為明顯，里海、黑海、紅海、澳大利亞卡奔塔利亞灣等封閉半封閉海域質量變化也較為明顯.同時能看到在長期趨勢較大的冰川冰蓋融化區域，海陸泄漏誤差極為嚴重，如南極、格陵蘭島、阿拉斯加、巴塔哥尼亞冰原等(谷延超, 2018).

為了消除海陸泄漏誤差，同時保證封閉、半封閉海域信號和地震海域信號的恢復，對圖5a的RMS信號做如下特殊處理：(1)保持陸地區域RMS信號不變；(2)保留質量變化明顯的海底地震區域RMS(蘇門答臘地震、東京地震和智利地震)以及封閉、半封閉海域RMS(里海、紅海、黑海、哈德孫灣、福克斯灣、澳大利亞卡奔塔利亞灣、泰國灣、北冰洋部分信號明顯海域)；(3)將其余海洋區域RMS均設置為3 cm(因為大多數海洋區域RMS小于3 cm)(Save et al., 2016; 谷延超, 2018).

圖5b為經以上特殊處理后的RMS信號分布圖，將此RMS格網點元素作為式(9)M-1的對角線元素，即獲得初始正則化矩陣.將以上初始化矩陣代入式(8)，用L曲線法尋找最佳正則化參數，即可獲得利用初始正則化矩陣求得的地表質量變化，稱為中間解.圖6a和6b分別為2010-01的無空間約束點質量解和利用初始正則化矩陣求得的中間解，通過對比可以看出：中間解在初始正則化矩陣的約束作用下基本消除海陸泄漏誤差，條帶誤差也大幅減弱；此外，中間解等效水高信號強度和空間分辨率也大大提高.信號增強現象在質量變化較大的區域最為明顯，如：亞馬遜流域、贊比西河流域以及南北極冰川冰蓋融化區域；但也能看出在質量變化較小的區域，信號幾乎沒有得到增強，如：青藏高原、華北平原和美國加利福利亞區域等.出現上述信號恢復情況不一致的原因是：全球質量變化信號復雜，僅使用單一的正則化參數不能同時兼顧質量變化強的信號和弱的信號的特征.

圖5 利用PM(Point-Mass，點質量)方法反演的2003-01—2014-11共133個月等效水高計算的(a)RMS分布圖和(b)初始正則化矩陣RMS分布圖Fig.5 (a) RMS distribution and (b) initial regularization matrix RMS distribution of equivalent water height calculated from January 2003 to November 2014 by PM

2.3 迭代正則化矩陣設計

圖5b的初始正則化矩陣RMS分布圖反映的是全球各陸地區域質量變化大小的分布，從圖中可以直觀地看出，全球質量變化大小分布極度不均勻，要消除質量變化較強區域的條帶誤差，需要較大的正則化參數進行約束，但對于信號較小的區域，過大的正則化參數會削弱其真實信號(Save et al., 2016)，因此單一的正則化參數很難平衡全球質量變化差異很大的情況.為兼顧全球質量變化信號較強區域的條帶誤差消除，以及質量變化較小區域的信號提取，可構建迭代的附有方差約束的徑向點質量方法，逐步提取時變重力場模型信號(Luthcke et al., 2013； Arendt et al., 2013).本文第一次迭代的正則化矩陣為利用133個月的中間解直接計算的RMS，其包含所有的質量變化信號；第二次用于迭代的正則化矩陣為扣除全球陸地區域等效水高季節性信號和高緯度冰蓋區域及海底地震區域長期趨勢信號后剩余信號計算的RMS，其主要包括低緯度地區趨勢信號較弱的質量變化信號；第三次迭代的是扣除所有陸地區域季節性信號和所有陸地及地震海域長期趨勢信號后剩余信號計算的RMS，其主要包括海洋信號和黑海、紅海等內陸海域信號計算的RMS.圖7a、7b和7c分別為三次迭代的正則化矩陣RMS分布圖.

圖6 分別利用單位矩陣和初始正則化矩陣的點質量方法解算的2010-01(a)無約束解和(b)中間解等效水高分布圖Fig.6 Equivalent water height distribution of January 2010 (a) unconstrained solution and (b) intermediate solution calculated by point-mass method using unit matrix and initial regularization matrix respectively

圖7 用于VCPM方法的不同迭代次數的正則化矩陣RMS分布圖Fig.7 Regularized matrix RMS distribution for different iterations of VCPM

其中第一次迭代所用的觀測值是式(1)直接計算出的重力擾動；第二次迭代的觀測值是利用式(6)計算出的第一次迭代剩余的重力擾動殘差，其包含的是第一次迭代未提取的重力擾動信號；第三次迭代的觀測值是前兩次迭代剩余的重力擾動殘差，包含的是前兩次迭代剩余的微小重力擾動信號.

迭代正則化矩陣按照上述方法設計的原因是：第一次迭代的正則化矩陣包含了所有質量變化特征，因此需要較大的正則化參數來抑制質量變化信號較大區域(如格陵蘭島)的條帶誤差，因此第一次迭代主要提取的是質量變化較大的信號，而質量變化較小的信號會留在重力擾動觀測值內，留待下一次迭代提取；第二次迭代的正則化矩陣由于扣除了高緯冰川融化趨勢和陸地振幅等信號，因此僅需要較小的正則化參數就能抑制條帶誤差，本次迭代就可以很容易的提取出重力擾動殘差中剩余的質量變化較小的信號(如華北平原地下水信號)；第三次迭代的正則化矩陣扣除了所有的陸地振幅和趨勢信號，因此可以使用更小的正則化參數提取前兩次未完全提取的海洋信號和內陸海洋信號，因此三次迭代提取的信號累加即可得到觀測值中所有的質量變化信號.

Save等(2016)利用扣除高緯冰川趨勢信號和陸地振幅信號的剩余信號RMS構造正則化矩陣，以此正則化矩陣一次性提取觀測值中蘊含的質量變化信號，因為這樣可以較好解決單一正則化參數不能平衡正則化矩陣中質量變化大小極不均勻的問題，但是這樣做也會導致先驗方差信息不真實的問題；而本文利用迭代的方法解決單一正則化參數的適應性問題，避免了先驗方差信息不真實的情況.

2.4 附有方差約束的徑向點質量方法反演結果

根據上述正則化矩陣的設計思路，附有方差約束的徑向點質量方法的計算流程圖如圖8所示，利用VCPM方法計算了2003-01—2014-11的全球等效水高變化，以2010-01為例，分析此方法單個月的反演效果.圖9為不同迭代次數的等效水高變化分布圖，其中圖9a，9b，9c分別是第一次到第三次迭代的反演結果，圖9d為三次迭代結果之和，也即VCPM計算的2010-01等效水高變化結果.為了反映不同迭代次數提取的重力擾動物理信號，圖10給出了2010-01重力擾動初始觀測值和不同迭代次數后的重力擾動殘差.可以看到:第一次迭代提取了重力擾動觀測值中包含的大多數質量變化信號，其中質量變化較大的信號(如南極、亞馬遜、奧卡萬戈三角洲等信號)已經基本提取完畢，但是從圖10b第一次迭代的重力擾動殘差中可以看出，質量變化較小的信號(如華北平原、青藏高原、堪察加半島和黑海)仍有部分未被提取；因此第二次迭代提取的主要是質量變化較小的信號，如：黑海、華北平原、青藏高原、非洲西部熱帶雨林等；第三次迭代主要提取的是前兩次迭代剩余的微小海洋信號和紅海、黑海等內陸海信號，此時圖10d第三次迭代后的重力擾動殘差中已基本不可見質量變化的物理信號，剩余的信號主要是條帶誤差，說明VCPM經過三次迭代能夠充分提取重力擾動觀測值中蘊含的質量變化信號.

圖8 附有方差約束的徑向點質量方法的計算流程圖Fig.8 Flow chart of variance-constrained radial point-mass method

為驗證三次迭代相對于一次迭代的優勢以及三次迭代提取的信號是否為真實信號，如圖11所示，分別利用一次迭代和三次迭代VCPM的2003-01—2014-11的全球質量變化計算了相對于CSRM和JPLM兩種Mascon數據的RMSE在20個流域的按面積加權平均值，另外也計算了CSRM和JPLM之間的RMSE在各流域的加權平均值，并將其作為參考值.從圖中可以看出，在除塔里木河流域之外的絕大多數流域，三次迭代的VCPM與兩種Mascon數據的RMSE平均值均要小于一次迭代的VCPM與兩種Mascon數據的RMSE平均值，這說明三次迭代的VCPM相對于一次迭代的VCPM要更加接近CSRM和JPLM這兩種Mascon數據，因此也能說明后兩次迭代提取了更多代表質量變化的真實信號，三次迭代能夠有效解決單一的正則化參數無法適應全球質量變化不均勻的問題.其次還可以看出三次迭代VCPM與CSRM的RMSE在各流域的加權平均值均要小于CSRM與JPLM的相應RMSE加權平均值，這說明VCPM與CSRM在這些流域之間的差異要小于CSRM與JPLM之間的差異，也能說明三次迭代的VCPM能夠達到與這兩種Mascon同等的精度.

圖9 利用VCPM方法反演的2010-01不同迭代次數和三次迭代之和的等效水高分布圖Fig.9 Equivalent water height distribution of different iterations and the sum of three iterations estimated by VCPM in January 2010

圖10 2010-01重力擾動初始值和不同迭代次數的重力擾動殘差Fig.10 Initial observation of gravity disturbation and residual of gravity disturbation with different iterations in January 2010

圖11 不同方法反演的2003-01—2014-11等效水高之間的RMSE在20個流域的按面積加權平均值柱狀圖Fig.11 Area-weighted mean values of RMSE between equivalent water heights estimated by different methods in 20 river basins from January 2003 to November 2014

圖12 VCPM與CSRM分別反演的2008年1、3、5月等效水高分布圖Fig.12 Equivalent water height distribution estimated by VCPM and CSRM methods in January, March, May 2008 respectively

為直觀展示VCPM反演質量變化的去條帶和泄漏誤差效果，將VCPM反演的等效水高與CSRM進行對比.圖12為 VCPM與CSRM分別反演的2008年1、3、5月等效水高分布圖，可以看出，兩種方法均已扣除條帶誤差及海陸泄漏誤差的影響，提高了信號精度和空間分辨率，并且反演的等效水高分布基本一致，信號強度在大多數區域基本相同.以上對比分析表明本文設計的VCPM方法相對于傳統濾波方法和PM方法更有優勢，在抑制重力場模型噪聲和獲取高精度高分辨率等效水高變化效果上與CSRM模型方法相當，后文將繼續在多個尺度上比較各反演結果優缺點.

3 不同方法反演結果比較與分析

為了具體評價VCPM反演地表質量變化的效果，引入CSRM、JPLM和GLDAS2.1水文模型進行對比分析.同時為了比較VCPM方法相對于傳統濾波方法的優勢，引入300 km高斯和去相關組合濾波算法(簡稱DS300)(Chen et al., 2007).由于下載的CSRM和JPLM數據中包含GAD產品，因此需要將其扣除；同時CSRM和JPLM均已經扣除冰川均衡調整GIA模型ICE-6G_D(VM5a)，因此VCPM和DS300也扣除該模型(Peltier et al., 2018).

3.1 全球尺度信號對比

為比較不同方法反演等效水高在多個月上的一致性，分別計算了2003-01—2014-11的VCPM與CSRM、JPLM、DS300數據和CSRM與JPLM、PM數據的均方根誤差(RMSE)，如圖13所示.同時表1給出了相應在全球、陸地和海洋區域的RMSE按面積加權平均值.從圖13a中可以看出，VCPM與CSRM一致性較好，海洋區域相差較小，除蘇門答臘、東京地震海域外的其余海域RMSE均很小；陸地區域中歐亞大陸、北美洲、澳大利亞和非洲北部等大部分地區RMSE很小，一致性較高；RMSE較大的區域主要分布在南極、格陵蘭島、阿拉斯加和巴塔哥尼亞等趨勢信號較大的區域，除此之外，恒河流域、印度河流域和亞馬遜流域北部RMSE也較大.CSRM與JPLM的RMSE分布圖和VCPM與CSRM的RMSE分布圖基本一致，RMSE較大區域與上述列舉的位置基本相同.另外從圖13c和13e中可以看出，DS300和PM在格陵蘭島等冰川冰蓋融化區域存在巨大泄漏誤差.

從表1的RMSE加權平均值來看，VCPM與CSRM的RMSE加權平均值在全球和陸地區域均最小(2.12 cm和4.16 cm)，優于CSRM與JPLM在全球和陸地區域的RMSE加權平均值(2.22 cm和4.62 cm).這可能是因為VCPM和CSRM均只采用GRACE觀測數據導出的正則化矩陣反演等效水高，沒有利用任何外部數據，而JPLM采用了GLDAS、衛星測高等外部物理模型設計正則化矩陣，另外VCPM與CSRM均為1°分辨率而JPLM為3°分辨率，因此在細節信號上VCPM與CSRM更為一致.在海洋區域RMSE加權平均值最小的是CSRM與JPLM，說明這兩種數據在海洋區域反演結果最接近；在海洋區域RMSE最大的是CSRM與PM，這是因為PM在最優正則化參數下條帶誤差仍然較大且存在泄漏誤差.VCPM與DS300在陸地區域的RMSE加權平均值最大，這是因為DS300由于平滑濾波導致了大量的信號泄漏.從以上對比分析可以看出，VCPM在利用先驗方差信息設計的正則化矩陣約束下能有效去除條帶誤差和抑制泄漏誤差.

表1 不同數據之間均方根誤差按面積加權的均值(cm)Table 1 Area-weighted mean value of root mean square error between different data (cm)

本文也利用不同數據計算了全球長時間序列等效水高的周年振幅和長期趨勢，如圖14和圖15所示.從周年振幅上來看，VCPM與CSRM和JPLM在空間分布和信號強度上表現出極好的一致性，尤其在亞馬遜流域、贊比西河流域、印度北部區域、湄公河流域、格陵蘭島、阿拉斯加等區域信號相似.DS300由于濾波平滑作用，其周年振幅被嚴重削弱，PM較DS300稍好，但其高緯度區域振幅信號也被嚴重削弱.從長期趨勢來看，DS300信號泄漏最為嚴重，很難分辨出信號較小的趨勢信號.PM除冰川冰蓋融化區域外，其余區域信號泄漏均較小，能清晰分辨出各種長期趨勢變化的細節信息，這反映出PM方法優于DS300方法.VCPM、CSRM和JPLM三種數據在長期趨勢上基本一致，其中VCPM和CSRM一致性略優于JPLM，這是因為地球質量變化信號中長期趨勢信號更加復雜，而JPLM空間分辨率為3°，VCPM和CSRM空間分辨率均為1°，能夠反映更多細節信號.從整體來看三者均能夠有效提取出地球上主要的質量變化信號，如南極、格陵蘭島、阿拉斯加、巴塔哥尼亞及亞洲高山等地區的冰雪消融信號；青藏高原降雨增加信號、印度西北部地下水枯竭、中國華北平原地下水消耗、美國加利福利亞地區地下水消耗和德克薩斯州持續干旱信號、里海海平面下降、非洲奧卡萬戈三角洲降雨增加信號等水儲量變化信號；2004年蘇門答臘地震、2010年智利地震和2011年東京地震信號(Chen et al., 2014; Scanlon et al., 2016; Peltier et al., 2018; Rodell et al., 2018;谷延超, 2018).從以上分析可以看出，在反演等效水高變化的周年振幅和長期趨勢方面，VCPM方法明顯優于DS300和PM方法，能夠達到官方發布的Mascon數據的效果.

圖13 不同方法反演的2003-01—2014-11全球等效水高之間的均方根誤差Fig.13 Root mean square errors between global equivalent water heights estimated from January 2003 to November 2014 by different methods

圖14 不同方法反演的2003-01—2014-11全球等效水高周年振幅Fig.14 Annual amplitude of global equivalent water height estimated by different methods from January 2003 to November 2014

圖15 不同方法反演的2003-01—2014-11全球等效水高長期趨勢Fig.15 Trends of global equivalent water height estimated by different methods from January 2003 to November 2014

3.2 陸地流域信號對比

為評價VCPM模型反演區域陸地水儲量變化的可靠性，選取了亞馬遜流域(Amazon)、奧里諾科河流域(Orinoco)、贊比西河流域(Zambeze)、育空河流域(Yukon)、印度河流域(Indus)、華北平原(NCP)等陸地流域和里海(Caspian)、紅海(Redsea)等內陸海域作為研究對象進行對比分析.圖16為不同數據在以上不同區域反演的2003-01—2014-11等效水高時間序列圖，表2給出了將所有數據扣除季節項和趨勢項后得到的CSRM與其余數據的時間序列相關系數，表3和表4分別給出了不同數據在各流域反演結果的長期趨勢和周年振幅.

表2 CSRM與其余不同數據均扣除季節項和趨勢項后的時間序列相關系數Table 2 Correlation coefficients of time series between CSRM and other different data after deducting the seasonal items and trend items

表3 不同方法反演的2003-01—2014-11等效水高在不同區域的長期趨勢(cm·a-1)Table 3 Trends of equivalent water height in different regions estimated by different methods from January 2003 to November 2014 (cm·a-1)

結合上述圖表可以看到：VCPM、CSRM和JPLM三種數據在各流域一致性都比較好，扣除季節項和趨勢項后的時間序列相關系數在除紅海之外的各區域仍能達到0.8以上，紅海區域相關系數較差的原因可能是該區域較為狹小，各方法反演結果的空間分布差異較大；DS300和PM在亞馬遜河、奧里諾科河、贊比西河流域區域與以上三種數據一致性較好，在華北平原、里海、紅海區域一致性較差.這是因為亞馬遜河、奧里諾科河、贊比西河流域均為季節性信號較強而且面積較大的區域，PM和DS300信號雖然會有泄漏誤差，但大多數信號都泄漏在區域內；而紅海區域雖然季節性信號占主要地位，但該區域面積較小且形狀狹長，因此PM和DS300信號極易泄漏到區域外；華北平原和里海為趨勢信號較強而季節性信號較弱區域，由于信號泄漏對趨勢信號影響較大，因此PM和DS300在該區域與VCPM、CSRM、JPLM一致性較差.在育空區域，VCPM與JPLM一致性更好，而CSRM與PM、DS300一致性較好，表現出了較大的負增長趨勢信號，這可能是阿拉斯加冰雪融化信號泄漏至該區域引起的.

圖16 不同方法反演的2003-01—2014-11等效水高在不同區域的時間序列Fig.16 Time series of equivalent water height in different regions estimated by different methods from January 2003 to November 2014

表4 不同方法反演的2003-01—2014-11等效水高在不同區域的周年振幅(cm)Table 4 Annual amplitude of equivalent water height in different regions estimated by different methods from January 2003 to November 2014 (cm)

GLDAS2.1模型反映的振幅信號和長期趨勢信號在多個區域與其他模型有較大差異，這是因為該模型主要反映的是淺層土壤水、地表積雪及植被冠層水等質量變化，而忽略了深層土壤水、地下水和地表湖泊等影響(詹金剛等, 2011).從華北平原、印度河流域的時間序列可以看出，GLDAS2.1時間序列趨勢和其他模型相反，表明這兩個區域存在地下水虧損的現象(Chen et al., 2014;魏偉等, 2020).

3.3 冰凍圈信號對比

由于全球氣候變暖導致南極和格陵蘭島存在巨大的冰雪消融趨勢信號，PM方法和傳統濾波方法在這兩個區域均存在巨大的泄漏誤差，因此本文利用各數據在南極和格陵蘭島進行對比分析以驗證VCPM消除泄漏誤差的能力.由于CSRM和JPLM均選用了ICE6G_D(VM05)模型進行GIA改正，因此本文其他模型也進行了同樣的GIA改正，其中DS300扣除的GIA模型經過了與之相同的截斷與濾波處理(Peltier et al., 2018).經計算，ICE6G_D(VM05)在南極的趨勢改正值約為98.8 Gt/a，在格陵蘭島改正值約為2.5 Gt/a.

圖17為不同方法反演的2003-01—2014-11南極區域質量變化長期趨勢，DS300由于濾波平滑作用存在很大泄漏誤差，其趨勢信號也被極大削弱；PM泄漏誤差也較為嚴重，反映出利用單位矩陣作為正則化矩陣在處理泄漏誤差方面的缺陷.VCPM由于先驗方差的約束作用極大程度地消除了南極區域的泄漏誤差，同時提高了南極區域質量變化的空間分辨率和信號精度，從南極半島、毛德皇后地、恩比德地、威爾克斯地等區域反映的趨勢信號來看，VCPM反演結果與CSRM和JPLM一致性極好.圖18為不同數據計算的南極區域質量變化時間序列，從整體來看VCPM、CSRM、JPLM三種數據時間序列基本一致，VCPM與CSRM的相關系數為0.997，與JPLM的相關系數為0.991，均有很好的相關性.由于DS300和PM兩種數據在南極區域泄漏誤差很大，因此對應時間序列與以上三種數據差異較大.通過對時間序列進行最小二乘擬合，得到了如表5所示的南極質量變化速度與加速度，其中VCPM模型反演的南極冰蓋融化速率為-139.8 Gt/a，與CSRM(-130.4 Gt/a)和JPLM(-123.1 Gt/a)的融化速率差值均在20 Gt/a以內.三種方法均探測到南極冰川呈現加速融化的趨勢，VCPM反演的加速度最大(-9.8 Gt/a2)，CSRM反演的加速度最小(-9.1 Gt/a2)，其差值也在8%以內，可能是因為VCPM利用迭代方法充分提取了觀測值中所包含的信號.

表5 不同方法反演的2003-01—2014-11南極質量變化速度與加速度Table 5 Velocity and acceleration of mass variations in Antarctic estimated by different methods from January 2003 to November 2014

圖19為不同方法反演的格陵蘭島質量變化長期趨勢和加速度空間分布圖，與南極區域類似，DS300和PM在格陵蘭島區域也存在很大泄漏誤差，而VCPM、CSRM和JPLM均有效地去除了泄漏誤差且趨勢空間分布較為一致，三種數據均反映出格陵蘭島外圍海岸線附近存在較大的質量虧損趨勢，且以東南和西北區域最為嚴重，是因為這兩個區域均存在流速較快的外流冰川；格陵蘭島北部、西北和東南部分質量損失趨勢比之較小，是因為這三個區域外流冰川較少，冰雪融化在其質量損失中占主要作用(Rignot and Mouginot, 2012；Velicogna et al., 2014).三種數據在格陵蘭島中部地區均表現出了質量積累趨勢，這是因為該地區地處高緯且海拔較高引起了降雪積累，但其空間分布和量級有所差異(Khan et al., 2015).圖20為不同方法計算的格陵蘭島區域質量變化時間序列，表6給出了相應的不同方法反演的質量變化速度與加速度.VCPM、CSRM和JPLM三種數據的時間序列一致性較好，JPLM質量變化速率最大為-308.3 Gt/a，CSRM和VCPM分別為-297.8 Gt/a和-288.6 Gt/a，三者差距在7%以內，其中VCPM與Velicogna等(2020)結果(-261 Gt/a)更為接近；在加速度方面，三模型反演結果無明顯差別.從以上分析可以看出利用質量塊質量變化的先驗方差對點質量方法進行合理約束可有效改正冰川區域的泄漏誤差，提升反演地表質量變化的效果.

圖17 不同方法反演的2003-01—2014-11南極等效水高長期趨勢Fig.17 Trends of equivalent water height in Antarctic estimated by different methods from January 2003 to November 2014

圖18 不同方法反演的2003-01—2014-11南極質量變化時間序列Fig.18 Time series of mass variations in Antarctic estimated by different methods from January 2003 to November 2014

圖19 不同方法反演的2003-01—2014-11格陵蘭島等效水高長期趨勢Fig.19 Trends of equivalent water height in Greenland estimated by different methods from January 2003 to November 2014

圖20 不同方法反演的2003-01—2014-11格陵蘭島質量變化時間序列Fig.20 Time series of mass variations in Greenland estimated by different methods from January 2003 to November 2014

表6 不同方法反演的2003-01—2014-11格陵蘭島質量變化速度與加速度Table 6 Velocity and acceleration of mass variations in Greenland estimated by different methods from January 2003 to November 2014

4 結論

本文利用附有方差約束的徑向點質量方法反演了2003-01—2014-11全球等效水高變化，結合PM、DS300、CSRM、JPLM和GLDAS2.1模型在多個尺度進行了對比分析，得出如下結論：

(1)PM在最佳正則化參數下的反演結果仍會受到條帶誤差影響，同時存在很大泄漏誤差.本文設計的附有方差約束的徑向點質量方法通過不同正則化矩陣迭代提取重力擾動觀測值中所蘊含的質量變化信號，能夠有效抑制海陸泄漏誤差，消除條帶誤差的影響，同時充分提取觀測信號中不同變化尺度的信號.

(2)在全球尺度的對比中，VCPM與CSRM的RMSE在全球和陸地尺度上最小，海洋區域僅略低于CSRM與JPLM的RMSE；VCPM、CSRM和JPLM反演的周年振幅和長期趨勢信號在空間分布和信號強度上較為一致，均能反映出高精度高分辨率的地球各區域質量變化物理信號，表明VCPM在全球尺度上精度與Mascon數據相當.

(3)在各陸地流域和兩個內陸海的時間序列對比中，VCPM與CSRM、JPLM時間序列一致性均較高，在扣除季節項和趨勢項后，除紅海之外的其他區域相關系數均能達到0.86以上，在以季節性信號為主的區域，相關系數可達0.98以上，紅海區域的相關系數較低可能是該區域較小且較為狹長，導致各方法反演的質量變化空間分布不一致引起的；同時發現泄漏誤差對于長期趨勢信號影響更為明顯，PM和DS300在長期趨勢較大區域信號泄漏更大.

(4)在南極和格陵蘭島兩大冰蓋區域，VCPM能夠有效抑制海陸泄漏誤差，在南極區域VCPM、CSRM和JPLM反演的質量減小趨勢分別為-139.8 Gt/a、-130.4 Gt/a和-123.4 Gt/a，在格陵蘭島區域分別為-288.6 Gt/a、-308.3 Gt/a和-297.8 Gt/a，VCPM在反演冰蓋質量變化方面相對于DS300和PM方法有極大提升，表明本文設計的方差約束可以有效提高點質量方法反演地表質量變化的效果.

致謝感謝CSR提供的GRACE時變重力場數據和Mascon數據，JPL提供的Mascon數據，NASA(National Aeronautics and Space Administration，美國國家航空航天局)提供的GLDAS2.1數據.感謝審稿專家提供的寶貴意見.