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

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

西南交通大學(xué)地球科學(xué)與環(huán)境工程學(xué)院， 成都 611756

0 引言

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

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

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

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

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

1 數(shù)據(jù)與方法

1.1 數(shù)據(jù)處理與驗證數(shù)據(jù)

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

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

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

1.2 空間格網(wǎng)

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

1.3 徑向點質(zhì)量模型方法基本原理

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

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

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

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

(2)

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

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

(4)

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

i=1,…,s.

(5)

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

y=Ax+e,(6)

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

1.4 參數(shù)估計方法

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

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

將式(7)對x求偏導(dǎo)數(shù)并聯(lián)合式(6)可得待估參數(shù)的解為：

(8)

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

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

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

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

本節(jié)利用CSR RL06 GSM數(shù)據(jù)通過正則化矩陣為單位矩陣的無空間約束徑向點質(zhì)量方法(以下簡稱PM)計算2003-01—2014-11共133個月的全球等效水高，圖3為PM反演各月份等效水高的L曲線圖，不同月份對應(yīng)的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曲線，當(dāng)正則化參數(shù)為1×10-19時，L曲線曲率最大，因此將其作為所有月份的正則化參數(shù).為了比較不同正則化參數(shù)的效果，同時選取了大于和小于最佳正則化參數(shù)反演結(jié)果作為對比，繪制了如圖4的2010-01等效水高分布圖，為了更直觀地對比反演效果，同時繪制了濾波半徑為300 km的高斯濾波反演結(jié)果分布圖.

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

2 附有方差約束的徑向點質(zhì)量方法

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

2.1 正則化矩陣

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

(9)

圖4 不同正則化參數(shù)的點質(zhì)量方法(a,b,c)和高斯濾波半徑為300 km的球諧位系數(shù)法(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個質(zhì)量塊的方差信息，通過計算第i個質(zhì)量塊多個月的RMS值得到.Save等(2016)設(shè)計了兩種正則化矩陣，一種是直接以2004—2009年等效水高RMS設(shè)計的0409F，另一種是在0409F基礎(chǔ)上再結(jié)合其他單個月的等效水高信息設(shè)計的附有時間相關(guān)性的0409V.結(jié)果表明：由于0409V包含了2004—2009年以外的陸地水儲量變化信息，因此以0409V反演的2004—2009年時間段以外的等效水高解效果優(yōu)于0409F反演結(jié)果，而在2004—2009年時間段以內(nèi)，兩個正則化矩陣反演結(jié)果相同.本文直接以2003—2014年全時段PM反演結(jié)果構(gòu)建正則化矩陣，包含了時間段內(nèi)的所有地表水儲量變化信號特征，因此不再構(gòu)建時間相關(guān)的正則化矩陣.

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

2.2 初始正則化矩陣設(shè)計

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

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

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

圖5 利用PM(Point-Mass，點質(zhì)量)方法反演的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 迭代正則化矩陣設(shè)計

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

圖6 分別利用單位矩陣和初始正則化矩陣的點質(zhì)量方法解算的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方法的不同迭代次數(shù)的正則化矩陣RMS分布圖Fig.7 Regularized matrix RMS distribution for different iterations of VCPM

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

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

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

2.4 附有方差約束的徑向點質(zhì)量方法反演結(jié)果

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

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

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

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

圖10 2010-01重力擾動初始值和不同迭代次數(shù)的重力擾動殘差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個流域的按面積加權(quán)平均值柱狀圖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反演質(zhì)量變化的去條帶和泄漏誤差效果，將VCPM反演的等效水高與CSRM進行對比.圖12為 VCPM與CSRM分別反演的2008年1、3、5月等效水高分布圖，可以看出，兩種方法均已扣除條帶誤差及海陸泄漏誤差的影響，提高了信號精度和空間分辨率，并且反演的等效水高分布基本一致，信號強度在大多數(shù)區(qū)域基本相同.以上對比分析表明本文設(shè)計的VCPM方法相對于傳統(tǒng)濾波方法和PM方法更有優(yōu)勢，在抑制重力場模型噪聲和獲取高精度高分辨率等效水高變化效果上與CSRM模型方法相當(dāng)，后文將繼續(xù)在多個尺度上比較各反演結(jié)果優(yōu)缺點.

3 不同方法反演結(jié)果比較與分析

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

3.1 全球尺度信號對比

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

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

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

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

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

表2 CSRM與其余不同數(shù)據(jù)均扣除季節(jié)項和趨勢項后的時間序列相關(guān)系數(shù)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等效水高在不同區(qū)域的長期趨勢(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)

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

圖16 不同方法反演的2003-01—2014-11等效水高在不同區(qū)域的時間序列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等效水高在不同區(qū)域的周年振幅(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模型反映的振幅信號和長期趨勢信號在多個區(qū)域與其他模型有較大差異，這是因為該模型主要反映的是淺層土壤水、地表積雪及植被冠層水等質(zhì)量變化，而忽略了深層土壤水、地下水和地表湖泊等影響(詹金剛等, 2011).從華北平原、印度河流域的時間序列可以看出，GLDAS2.1時間序列趨勢和其他模型相反，表明這兩個區(qū)域存在地下水虧損的現(xiàn)象(Chen et al., 2014;魏偉等, 2020).

3.3 冰凍圈信號對比

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

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

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

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

圖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南極質(zhì)量變化時間序列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格陵蘭島質(zhì)量變化時間序列Fig.20 Time series of mass variations in Greenland estimated by different methods from January 2003 to November 2014

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

4 結(jié)論

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

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

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

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

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

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