鐘波， 李賢炮， 李建成*， 汪海洪,2， 丁劍

1 武漢大學(xué)測繪學(xué)院， 武漢 430079 2 地球空間環(huán)境與大地測量教育部重點實驗室， 武漢 430079 3 湖北珞珈實驗室， 武漢 430079 4 中國測繪科學(xué)研究院， 北京 100830

0 引言

全球定位系統(tǒng)(Global Positioning System, GPS)作為現(xiàn)代大地測量的重要觀測手段，被廣泛應(yīng)用于參考框架確定、板塊構(gòu)造動力學(xué)機制研究和地表形變監(jiān)測等科學(xué)問題.根據(jù)質(zhì)量負(fù)荷格林函數(shù)理論(Farrell, 1972)，地表質(zhì)量負(fù)荷(如地表水、地下水和冰川等)的遷移與重新分布會導(dǎo)致地殼發(fā)生形變.例如，Heki(2001)研究發(fā)現(xiàn)日本東北地區(qū)地表積雪的季節(jié)性變化會使GPS站點產(chǎn)生抬升和沉降.由于地表形變可以被GPS以毫米級精度連續(xù)地進(jìn)行監(jiān)測，因此可以利用GPS形變觀測數(shù)據(jù)獨立地反演地表質(zhì)量變化(Blewitt et al., 2001).特別是，很多國家和地區(qū)建立了密集的GPS連續(xù)運行臺網(wǎng)對區(qū)域地表形變進(jìn)行監(jiān)測(姜衛(wèi)平，2017)，這為反演精細(xì)的區(qū)域地表質(zhì)量變化提供了有力的數(shù)據(jù)保障.

扣除大氣和海洋質(zhì)量變化以后，地表質(zhì)量變化主要為陸地水儲量變化(Terrestrial Water Storage Change，TWSC).TWSC是指儲存在地表以及地下的全部水分，包括地表水、地下水、土壤水、樹冠水以及積雪等，是水循環(huán)的重要組成部分(Rodell and Famiglietti, 1999; Cazenave and Chen, 2010; Rodell et al., 2018).由于GPS垂直位移相比水平位移能更準(zhǔn)確地反映質(zhì)量負(fù)荷變化信息(Wahr et al., 2013)，當(dāng)前大多數(shù)研究是利用GPS垂直位移對區(qū)域TWSC進(jìn)行反演.Argus 等(2014)利用密集的GPS垂直位移數(shù)據(jù)反演了美國加州地區(qū)的TWSC，結(jié)果表明其空間分辨率可達(dá)到～50 km，并且與重力恢復(fù)與氣候?qū)嶒炐l(wèi)星(Gravity Recovery and Climate Experiment，GRACE)觀測結(jié)果及水文模型符合較好.隨后，GPS被廣泛用于反演不同地區(qū)的TWSC，如美國中西部地區(qū)(Borsa et al., 2014; Fu et al., 2015; Jin and Zhang, 2016; Argus et al., 2017)、中國西南地區(qū)(何思源等, 2018; Zhong et al., 2020; 成帥等, 2021; Jiang et al., 2021a,b; Shen et al., 2021)以及南美洲(Ferreira et al., 2019)等典型區(qū)域.其中，Borsa 等(2014)利用GPS垂直位移反演的TWSC研究了美國西部的干旱情況，結(jié)果表明GPS反演的水量流失總量與水文氣象觀測結(jié)果一致；Fu等(2015)利用GPS垂直位移反演了美國華盛頓和俄勒岡州的TWSC，結(jié)果顯示GPS反演結(jié)果與GRACE及水文模型結(jié)果具有高度相關(guān)性；Jin和Zhang(2016)利用GPS垂直位移反演的TWSC有效監(jiān)測了2012年美國西南部的干旱事件，表明GPS形變觀測具有干旱指示器的作用；Argus等(2017)利用GPS垂直位移估算了2012—2015年美國加州山區(qū)持續(xù)干旱期間的水儲量流失情況；何思源等(2018)利用GPS垂直位移反演了云南省TWSC的時空分布；Ferreira等(2019)探討了利用每日間隔的GPS垂直位移反演南美地區(qū)TWSC的空間分辨率；Jiang等(2021a)利用每天的GPS垂直位移反演了云南省日間隔的TWSC，并追蹤了極端天氣事件導(dǎo)致的水儲量變化；Shen 等(2021)提出了將隨機森林算法和地殼負(fù)荷模型相結(jié)合的機器學(xué)習(xí)負(fù)荷反演方法，有效提高了利用稀疏GPS測站數(shù)據(jù)反演TWSC的精度.此外，利用GPS反演的TWSC還可作為GRACE及其繼任者GFO(GRACE Follow-on)觀測結(jié)果的有益補充，特別是可以有效彌補GRACE/GFO缺失月份和任務(wù)間隔期的TWSC信息(Fu et al., 2015；Zhong et al., 2020).

利用GPS垂直位移反演區(qū)域TWSC屬于典型的病態(tài)問題，常規(guī)方法是基于Tikhonov正則化方法(Tikhonov, 1963)并利用二階拉普拉斯算子作為正則化約束矩陣對病態(tài)問題進(jìn)行求解(Argus et al., 2014; Fu et al., 2015; Jin and Zhang, 2016; 何思源等, 2018; Zhong et al., 2020; Jiang et al., 2021a,b; Shen et al., 2021)，或者是基于截斷奇異值分解(Truncated Singular Value Decomposition，TSVD)的正則化方法求解病態(tài)問題(Lai et al., 2020).但這兩種正則化方法各有優(yōu)勢，其中Tikhonov正則化是通過引入先驗約束矩陣并結(jié)合正則化參數(shù)選取方法對病態(tài)問題進(jìn)行求解，有效提高了反演結(jié)果的穩(wěn)定性；而TSVD正則化是通過舍棄設(shè)計矩陣中貢獻(xiàn)較小的奇異值并僅由觀測數(shù)據(jù)進(jìn)行獨立求解，其反演結(jié)果不依賴于先驗信息且計算簡單.因此，通過組合TSVD和Tikhonov正則化方法(即TSVD-Tikhonov)可實現(xiàn)兩者優(yōu)勢互補，即利用TSVD正則化舍棄小奇異值對反演結(jié)果的影響，并利用Tikhonov正則化引入先驗約束以提高解的穩(wěn)定性，進(jìn)而得到比單獨正則化解更穩(wěn)定和可靠的反演結(jié)果.例如，Chen T Y等(2020)基于TSVD-Tikhonov正則化方法利用GRACE時變重力場模型位系數(shù)對青藏高原地區(qū)的地表質(zhì)量變化進(jìn)行估計，結(jié)果顯示TSVD-Tikhonov正則化比單獨利用TSVD或Tikhonov正則化的反演結(jié)果更優(yōu).

盡管TSVD-Tikhonov組合正則化方法是一種更優(yōu)的病態(tài)問題求解方案，但其并未有效用于解決GPS反演區(qū)域TWSC的病態(tài)問題.本文以四川省為實驗區(qū)，利用中國地殼運動觀測網(wǎng)絡(luò)(Crustal Movement Observation Network of China，CMONOC)發(fā)布的72個GPS測站的垂直位移數(shù)據(jù)，基于TSVD-Tikhonov正則化方法反演2010年12月至2021年2月的TWSC時間序列.首先，通過數(shù)值模擬對TSVD、Tikhonov和TSVD-Tikhonov正則化方法采用不同正則化參數(shù)選取策略進(jìn)行TWSC反演分析，以驗證TSVD-Tikhonov正則化方法的有效性.其次，基于TSVD-Tikhonov正則化方法利用實測的GPS垂直位移反演四川省的TWSC時間序列，并與JPL(Jet Propulsion Laboratory)、CSR(Center for Space Research)和GSFC(Goddard Space Flight Center)發(fā)布的GRACE/GFO RL06 版本的Mascon模型進(jìn)行對比.最后，利用廣義三角帽方法(Generalized Three-Cornered Hat，GTCH)融合不同類型的水文氣象資料得到更為可靠的降水、蒸散發(fā)和徑流數(shù)據(jù)，并根據(jù)水量平衡方程計算的四川省dTWSC/dt時間序列對GPS和GRACE/GFO反演結(jié)果進(jìn)行驗證.

1 理論與方法

1.1 負(fù)荷格林函數(shù)理論

根據(jù)彈性負(fù)荷理論(Farrell，1972)，地表質(zhì)量負(fù)荷變化(如TWSC)引起的地球形變可以近似表示為彈性形變(包括垂直和水平方向).其中，垂直形變可以表示為質(zhì)量負(fù)荷與垂直位移格林函數(shù)的卷積(Wahr et al., 2013; Jin and Zhang, 2016; Zhong et al., 2020)：

(1)

式中，U(φ)為垂直形變；ΔM為點質(zhì)量負(fù)荷；R和Me分別為地球的半徑和質(zhì)量；hl為負(fù)荷Love數(shù)，本文采用Wang等(2012)的計算結(jié)果；Pl為l階勒讓德多項式；φ為點質(zhì)量負(fù)荷與GPS測站之間的角距.

根據(jù)式(1)，GPS垂直位移與TWSC之間的觀測方程可以表示為：

y=Ax+e,e～(0,σ2I),

(2)

式中，y為GPS垂直位移觀測向量；A為格林函數(shù)設(shè)計矩陣；x為待求的TWSC向量，通常采用等效水高(Equivalent Water Height，EWH)表示；e和σ2分別是GPS垂直位移y的觀測值殘差向量和誤差方差；I是與y相對應(yīng)的單位矩陣.

1.2 正則化方法

由于相鄰GPS測站的垂直位移觀測值之間存在相關(guān)性，以及GPS測站個數(shù)通常遠(yuǎn)小于待估參數(shù)個數(shù)，這些都會導(dǎo)致式(2)求解的法方程矩陣條件數(shù)過大，因此利用GPS垂直位移反演區(qū)域TWSC屬于典型的離散病態(tài)問題.以下簡單介紹病態(tài)問題求解的Tikhonov和TSVD正則化方法，并給出兩者組合的TSVD-Tikhonov正則化方法.

1.2.1 Tikhonov正則化

Tikhonov正則化方法是通過引入先驗正則化約束矩陣對病態(tài)問題進(jìn)行穩(wěn)定求解.根據(jù)式(2)的觀測方程，正則化約束求解TWSC向量x的目標(biāo)函數(shù)為(Argus et al., 2014)：

(3)

式中，λ>0為正則化參數(shù)，其選取方法主要有L-curve法(Hansen, 1992)、廣義交叉檢驗法(Generalized Cross Validation，GCV)(Golub et al., 1979)以及基于先驗信息的均方根誤差(Root Mean Square Error，RMSE)最小準(zhǔn)則等；L為正則化約束矩陣，本文采用改進(jìn)的正則化拉普拉斯矩陣(李賢炮等，2022)進(jìn)行求解，該矩陣通過考慮邊緣格網(wǎng)點與內(nèi)部格網(wǎng)點之間的相關(guān)關(guān)系，可更好地抑制邊緣效應(yīng)的影響.于是，式(2)中TWSC參數(shù)向量的估值為：

(4)

需要說明的是，式(4)中GPS垂直位移的中誤差σ通常難以準(zhǔn)確獲得，解算中可將σλ作為正則化參數(shù)進(jìn)行估計.

1.2.2 TSVD正則化

TSVD正則化方法是通過對觀測方程的設(shè)計矩陣進(jìn)行奇異值分解，并利用貢獻(xiàn)較大的奇異值及相應(yīng)奇異向量對設(shè)計矩陣進(jìn)行重構(gòu)，再根據(jù)重構(gòu)的觀測方程獨立求解待估參數(shù).

對式(2)中格林函數(shù)設(shè)計矩陣A進(jìn)行奇異值分解，可得(嵇昆浦和沈云中，2020)：

(5)

式中，Q=(q1,q2,…,qm)和V=(v1,v2,…,vn)分別為設(shè)計矩陣A的左、右奇異向量，m和n分別為觀測值個數(shù)和待估參數(shù)個數(shù)；Σ=diag(σ1,σ2,…,σm)，σ1>σ2>…>σm>0為設(shè)計矩陣A的奇異值.

于是，由奇異值分解得到式(2)中TWSC參數(shù)向量的最小二乘估計解為：

(6)

(7)

1.2.3 TSVD-Tikhonov正則化

(8)

1.3 廣義三角帽方法

廣義三角帽方法(GTCH)最大的優(yōu)點是可在缺少數(shù)據(jù)先驗信息的情況下對多組數(shù)據(jù)集的不確定性進(jìn)行評估(Tavella and Premoli, 1994; Galindo et al., 2001)，已被廣泛用于估計不同GRACE時變重力場模型及陸面模型的不確定度(Long et al., 2017; 姚朝龍等，2019；Zhao et al., 2019).因此，本文采用GTCH方法對不同數(shù)據(jù)集(或模型)的不確定度進(jìn)行估計，并根據(jù)估計的方差進(jìn)行最優(yōu)融合，進(jìn)而得到更加可靠的結(jié)果.

針對不同的數(shù)據(jù)集{Xi,i=1,2,…,N}，其加權(quán)融合的結(jié)果可以表示為：

(9)

其中，pi為數(shù)據(jù)集Xi對應(yīng)的權(quán)重，具體計算方法為：

(10)

式中，ωi為GTCH方法估計得到的數(shù)據(jù)集Xi的方差，具體計算方法參見姚朝龍等(2019)和Zhao 等(2019).

2 數(shù)值模擬

2.1 實驗區(qū)概況

四川省位于我國西南地區(qū)，地處長江上游，界于北緯26°03′—34°19′和東經(jīng)97°21′—108°12′之間，全省面積約為4.86×105km2，自然資源豐富，地形復(fù)雜多樣，高差懸殊較大，地勢呈西高東低的特點(見圖1).四川省內(nèi)氣候溫暖濕潤，冬暖夏熱，包含亞熱帶季風(fēng)氣候、高原山地氣候以及高原高寒氣候等多種氣候類型，年平均氣溫約為16.4 ℃，多數(shù)地區(qū)的年降水量在900～1200 mm之間.在全球氣候變暖的背景下，監(jiān)測四川省TWSC的時空變化對于該地區(qū)的水循環(huán)研究及自然資源配置和經(jīng)濟發(fā)展具有重要意義.

為了更好的反演TWSC，并降低邊緣效應(yīng)的影響，本文將解算區(qū)域從四川省的邊界向外擴展了2°，并使用了整個區(qū)域內(nèi)共72個連續(xù)運行的GPS測站垂直位移數(shù)據(jù)進(jìn)行計算.這些GPS測站的位置分布情況如圖1所示，其中南部和東北部的測站相對密集，而西北部和東南部的測站相對稀疏，整個計算區(qū)域平均每10000 km2約有38.6個GPS測站，測站之間的平均距離約為169.6 km.

圖1 實驗區(qū)的地形及GPS測站位置分布Fig.1 Topography of the test area and distribution of GPS stations

2.2 不同正則化方法的模擬結(jié)果

為了驗證TSVD-Tikhonov正則化方法的有效性，利用模擬的GPS垂直位移數(shù)據(jù)，基于TSVD、Tikhonov和TSVD-Tikhonov正則化方法采用不同的正則化參數(shù)選取策略進(jìn)行反演.數(shù)值模擬中，將研究區(qū)域劃分為0.5°×0.5°格網(wǎng)，以GLDAS(Global Land Data Assimilation Systems)水文模型(Rodell et al., 2004)計算的2005年1—12月的TWSC作為原始信號，根據(jù)式(1)模擬研究區(qū)域內(nèi)72個GPS測站的垂直位移觀測值，并在垂直位移數(shù)據(jù)中加入均值為0 mm、標(biāo)準(zhǔn)差為1 mm的高斯白噪聲，然后分別由式(4)、式(7)和式(8)采用不同的正則化參數(shù)選取策略(RMSE最小準(zhǔn)則、GCV法和L-curve法)反演TWSC，并將其與原始信號進(jìn)行對比分析.

圖2a和2b分別為基于Tikhonov和TSVD正則化方法采用不同的正則化參數(shù)選取策略反演的2005年1—12月的TWSC及與原始信號的區(qū)域平均時間序列比較.可以看出，Tikhonov和TSVD正則化方法的反演結(jié)果都能夠與原始信號保持較好的一致性，但Tikhonov正則化對應(yīng)的反演結(jié)果與原始信號更為接近，而TSVD的反演結(jié)果與原始信號的整體差異相對較大(如2005年9月)，這表明Tikhonov正則化相比TSVD正則化的反演結(jié)果更為穩(wěn)定和可靠.

圖2 基于TSVD和Tikhonov正則化由不同正則化參數(shù)選取方法反演的TWSC(a和b)及其與原始信號的差值標(biāo)準(zhǔn)差(c和d)和相應(yīng)的最優(yōu)正則化參數(shù)與截斷階數(shù)(e和f)Fig.2 TWSC solved by Tikhonov and TSVD regularizations through different regularization parameter selection methods (a and b), STD of differences between the estimated TWSC and original signal (c and d), and the corresponding optimal regularization parameters and truncation orders (e and f)

圖2c和2d分別為基于Tikhonov和TSVD正則化方法由不同的正則化參數(shù)選取策略反演的TWSC與原始信號差值的標(biāo)準(zhǔn)差(Standard Deviation，STD).可以看出，在Tikhonov正則化的反演結(jié)果中，RMSE最小準(zhǔn)則對應(yīng)的STD最小，GCV法對應(yīng)的STD次之，而L-curve法對應(yīng)的STD最大；在TSVD正則化的反演結(jié)果中，RMSE最小準(zhǔn)則對應(yīng)的STD同樣為最小，而L-curve法對應(yīng)的STD在大多數(shù)月份小于GCV法對應(yīng)的STD，但GCV法相比L-curve法對應(yīng)的STD變化更為穩(wěn)定.總體而言，TSVD正則化對應(yīng)的STD均要大于Tikhonov正則化對應(yīng)的STD，這也表明Tikhonov正則化比TSVD正則化的反演結(jié)果精度更高.

圖2e和2f分別為Tikhonov和TSVD正則化采用RMSE最小準(zhǔn)則、GCV法和L-curve法反演TWSC所選取的最優(yōu)正則化參數(shù)和截斷階數(shù).對于Tikhonov正則化而言，GCV法與RMSE最小準(zhǔn)則確定的正則化參數(shù)更為接近，而L-curve法確定的正則化參數(shù)在8—12月出現(xiàn)較大差異，其對應(yīng)的誤差STD也相對較大(見圖2c)，這是因為L-curve法有時難以準(zhǔn)確獲取最優(yōu)的正則化參數(shù)；對于TSVD正則化，L-curve法與RMSE最小準(zhǔn)則確定的最優(yōu)截斷階數(shù)相近的月份更多，但L-curve法確定的最優(yōu)截斷階數(shù)逐月變化差異較大，而GCV法確定的最優(yōu)截斷階數(shù)的變化更為平穩(wěn).

為了對TSVD-Tikhonov正則化的截斷階數(shù)進(jìn)行分析，圖3a和3b分別給出了格林函數(shù)設(shè)計矩陣A的奇異值和不同截斷階數(shù)對應(yīng)的累積貢獻(xiàn)率，可知設(shè)計矩陣A的奇異值在40～50階以后的變化趨于緩慢，其對應(yīng)的累積貢獻(xiàn)率在80%以上.例如，截斷至44階對應(yīng)的累積貢獻(xiàn)率達(dá)到了87.75%.

圖3 設(shè)計矩陣奇異值(a)和不同截斷階數(shù)對應(yīng)的累積貢獻(xiàn)率(b)Fig.3 Singular values of design matrix (a) and their cumulative contribution ratio (b) with different truncation orders

由于Tikhonov正則化比TSVD正則化的反演結(jié)果更可靠(見圖2)，因此圖4a和4b以Tikhonov正則化反演結(jié)果作為參考，分別給出了Tikhonov正則化和采用不同截斷階數(shù)的TSVD-Tikhonov組合正則化反演結(jié)果對應(yīng)的誤差STD.圖4a中矩形點線是Tikhonov正則化采用RMSE最小準(zhǔn)則確定最優(yōu)正則化參數(shù)的反演結(jié)果與原始信號差值的STD；圓形點線是TSVD-Tikhonov正則化采用RMSE最小準(zhǔn)則選取TSVD截斷階數(shù)、GCV法確定Tikhonov最優(yōu)正則化參數(shù)所對應(yīng)的反演結(jié)果與原始信號差值的STD.可以看出，除個別月份(7月)以外，TSVD-Tikhonov正則化比Tikhonov正則化反演結(jié)果對應(yīng)的STD更小，說明TSVD-Tikhonov正則化方法的反演結(jié)果更優(yōu).

圖4 Tikhonov正則化和不同截斷階數(shù)K的TSVD-Tikhonov正則化反演結(jié)果的誤差STD(a)和月平均誤差STD(b)Fig.4 STD (a) and monthly mean STD (b) of the errors of the inversion results based on Tikhonov regularization and TSVD-Tikhonov regularization with different truncation order K

圖4b中虛線和實線分別為圖4a中Tikhonov和TSVD-Tikhonov正則化反演結(jié)果對應(yīng)STD的月平均值，而圓形點線則為不同截斷階數(shù)的TSVD-Tikhonov正則化反演結(jié)果對應(yīng)的月平均STD值.可以看出，在截斷階數(shù)為42階以后，TSVD-Tikhonov正則化比Tikhonov正則化反演結(jié)果對應(yīng)的月平均STD更小，并且截斷階數(shù)為44階時達(dá)到最優(yōu)，相應(yīng)的累積貢獻(xiàn)率為87.75%(見圖3b).而在TSVD-Tikhonov組合正則化過程中，每個月都以RMSE最小準(zhǔn)則選取TSVD最優(yōu)截斷階數(shù)，并采用GCV法確定最優(yōu)正則化參數(shù)，其反演結(jié)果對應(yīng)的月平均STD為最小(圖4b中實線)，可認(rèn)為是理論最優(yōu)結(jié)果.

為了進(jìn)一步評估TSVD-Tikhonov正則化方法的有效性，圖5給出了TSVD、Tikhonov和TSVD-Tikhonov正則化反演的2005年1月至12月TWSC序列與原始信號差值STD的空間分布.可以看出，在整個四川省區(qū)域內(nèi)，TSVD正則化反演結(jié)果的STD最大，平均STD為14.97 mm；Tikhonov正則化與TSVD-Tikhonov正則化對應(yīng)的STD相對較小且總體接近，但在四川省東部Tikhonov正則化對應(yīng)的STD比TSVD-Tikhonov正則化對應(yīng)的STD更大，其相應(yīng)的平均STD分別為7.03 mm和5.04 mm.也就是說，TSVD-Tikhonov正則化對應(yīng)的平均STD在三者之中最小，即TSVD-Tikhonov正則化方法反演TWSC的精度最高.

圖5 TSVD (a)、Tikhonov (b)和TSVD-Tikhonov (c)正則化反演的TWSC序列與原始信號差值STD的空間分布Fig.5 Spatial distributions of STD of the differences between the estimated TWSC series and original signal based on TSVD (a), Tikhonov (b), and TSVD-Tikhonov (c) regularizations

以上數(shù)值模擬結(jié)果表明，TSVD-Tikhonov正則化方法比單獨利用TSVD或Tikhonov正則化方法反演的四川省TWSC的精度更高、穩(wěn)定性也更好，初步驗證了TSVD-Tikhonov組合正則化方法的有效性.需要指出的是，實測數(shù)據(jù)處理中無法獲得可靠的先驗信息由RMSE最小準(zhǔn)則選取TSVD對應(yīng)的最優(yōu)截斷階數(shù)，因此在后續(xù)實測數(shù)據(jù)反演中以數(shù)值模擬得出的最優(yōu)階數(shù)(44階)作為TSVD-Tikhonov正則化的截斷階數(shù).同時，綜合考慮GCV法和L-curve法的數(shù)值穩(wěn)定性及反演精度(見圖2)，在后續(xù)實測數(shù)據(jù)計算中采用GCV法確定最優(yōu)正則化參數(shù).

3 反演建模與結(jié)果分析

3.1 數(shù)據(jù)來源及預(yù)處理

3.1.1 GPS垂直位移數(shù)據(jù)

利用CMONOC的72個GPS測站(位置分布見圖1)獲取的2010年12月至2021年2月期間的垂直位移時間序列進(jìn)行TWSC反演(數(shù)據(jù)下載自國家地震科學(xué)數(shù)據(jù)中心http:∥www.eqdsc.com).CMONOC提供的GPS單日位移時間序列相對于ITRF2008參考框架，是由GAMIT/GLOBK軟件(10.4版)基于雙差模型解算得到(詳細(xì)說明參見官網(wǎng)數(shù)據(jù)處理手冊).由于CMONOC發(fā)布的GPS位移時間序列已經(jīng)扣除了固體潮、海潮以及極潮等的影響，因此獲取的GPS垂直位移主要反映的是陸地水、非潮汐大氣和海洋質(zhì)量變化引起的垂直形變.此外，由于人類活動(Yao et al., 2020)、設(shè)備更換以及自然環(huán)境變化等因素的影響，GPS位移時間序列會存在數(shù)據(jù)跳躍和間斷，因此需要對GPS垂直位移數(shù)據(jù)進(jìn)行預(yù)處理，主要包括粗差探測和數(shù)據(jù)插值等.本文利用四分位粗差探測方法(Nikolaidis, 2012)剔除粗差，并采用克里金-卡爾曼濾波方法(Liu et al., 2018)對缺失的垂直位移數(shù)據(jù)進(jìn)行插值補齊.圖6給出了四川巴中(SCBZ)和四川鹽源(SCYY)兩個測站(位置見圖1)的垂直位移時間序列在濾波前后的對比，圖中顯示濾波處理可以很好地去除高頻噪聲并填補缺失的數(shù)據(jù).

圖6 SCBZ和SCYY測站的垂直位移時間序列濾波前后比較Fig.6 Comparisons of the vertical displacement time series at SCBZ and SCYY stations before and after filtering

另外，GPS垂直位移時間序列的長期趨勢中包含了由地表質(zhì)量變化、冰川均衡調(diào)整(Glacial Isostatic Adjustment，GIA)和板塊運動等地球物理現(xiàn)象的綜合貢獻(xiàn)，但卻難以將GIA和板塊運動等引起的長期趨勢準(zhǔn)確地分離出來.考慮到研究區(qū)域的地表質(zhì)量變化以季節(jié)性特征為主，因此在GPS垂直位移時間序列的預(yù)處理過程中扣除了長期趨勢項，保留了周年和半周年變化信號.此外，為了最終反演得到TWSC序列，利用非潮汐大氣與海洋去混頻產(chǎn)品(Atmosphere and Ocean De-aliasing Level-1B, AOD1B)RL06版本扣除了GPS垂直位移中由非潮汐大氣和海洋效應(yīng)導(dǎo)致的地表形變影響.同時，為了與GRACE/GFO Mascon模型的時間尺度保持一致，本文將GPS垂直位移數(shù)據(jù)取月平均，然后用于TWSC的反演.

3.1.2 GRACE/GFO Mascon模型

為了對GPS反演結(jié)果進(jìn)行比較，采用CSR(Save et al., 2016)、JPL(Watkins et al., 2015; Wiese et al., 2016)和GSFC(Loomis et al., 2019)發(fā)布的GRACE/GFO RL06版本的Mascon模型，分別計算了四川省2010年12月至2017年6月和2018年6月至2021年2月期間的TWSC時間序列.需要說明的是，Mascon模型通過采用約束方法對GRACE/GFO觀測中含有的南北條帶誤差進(jìn)行抑制，具有更高的空間分辨率，因此可以使用Mascon模型直接計算TWSC，并不需要對其進(jìn)行濾波處理和信號泄漏改正.此外，由于GRACE/GFO反演的低階項系數(shù)存在較大的不確定性，Mascon模型中對應(yīng)的C20和C30(只針對GFO解)項已使用衛(wèi)星激光測距(Satellite Laser Ranging，SLR)的估計結(jié)果進(jìn)行了替換，并且地心運動導(dǎo)致的一階項變化和GIA等效應(yīng)也進(jìn)行了相應(yīng)改正，參見Save 等(2016)和Watkins 等(2015).

3.1.3 降水、蒸散發(fā)和徑流數(shù)據(jù)

利用不同類型水文氣象資料得到的降水、蒸散發(fā)和徑流數(shù)據(jù)對GPS和GRACE/GFO反演結(jié)果進(jìn)行檢驗.這些資料包括國家氣象信息中心NMIC(National Meteorological Information Center)提供的地面降水月值0.5°×0.5°格網(wǎng)數(shù)據(jù)集V2.0(http:∥data.cma.cn/)；熱帶降雨測量任務(wù)TRMM(Tropical Rainfall Measurement Mission)(Huffman et al., 2007)提供的月值0.25°×0.25°格網(wǎng)降水?dāng)?shù)據(jù)(https:∥disc.gsfc.nasa.gov/datasets/TRMM_3B43_7/summary)；全球陸地蒸散發(fā)模型GLEAM(Global Land Evaporation Amsterdam Model)(Miralles et al., 2011; Martens et al., 2017)提供的月值0.25°×0.25°格網(wǎng)蒸散發(fā)數(shù)據(jù)(https:∥www.gleam.eu/)；英國國家大氣科學(xué)中心NCAS(National Centre for Atmospheric Science)研制的氣候產(chǎn)品CRU(Climatic Research Unit)計算的月值0.5°×0.5°格網(wǎng)降水和蒸散發(fā)數(shù)據(jù)(https:∥crudata.uea.ac.uk/cru/data/hrg/)；以及GLDAS的三個陸面模型Noah、CLSM和VIC(Rodell et al., 2004)計算的月尺度降水、蒸散發(fā)和徑流數(shù)據(jù)(https:∥earthdata.nasa.gov/search?q=GLDAS)，其中Noah數(shù)據(jù)為0.25°×0.25°格網(wǎng)，CLSM和VIC數(shù)據(jù)為0.5°×0.5°格網(wǎng).

對于一個特定的區(qū)域，其水循環(huán)主要是由降水、蒸散發(fā)和徑流控制(Chen J L et al., 2020).于是，區(qū)域內(nèi)的水量平衡方程可以表示為(Famiglietti et al., 2011; 李瓊等，2013)：

(11)

式中，P為降水，ET為蒸散發(fā)，R為徑流，dS/dt為TWSC關(guān)于時間t的一階導(dǎo)數(shù)，即為dTWSC/dt.同時，由GPS或GRACE/GFO反演的TWSC序列計算dS/dt的公式可表示為：

(12)

式中，dSt為利用GPS或GRACE/GFO反演的連續(xù)兩個月之間的TWSC，這里dt為1個月.

因此，可以采用水量平衡方程計算的月平均dTWSC/dt對GPS和GRACE/GFO計算的月平均dTWSC/dt進(jìn)行檢驗.但水量平衡方程(式(11))計算的月平均dS/dt對應(yīng)的時間為月中，而式(12)采用中心差分計算的GPS和GRACE/GFO月平均dSt/dt對應(yīng)的時間為月初或月末.為了使兩者在時間上保持一致，本文首先采用Landerer 等(2010)提出的時間序列濾波方法分別對水量平衡方程中的降水、蒸散發(fā)和徑流數(shù)據(jù)進(jìn)行濾波處理，即：

(13)

3.2 棋盤測試結(jié)果

利用GPS垂直位移反演TWSC的分辨率和可靠性還與計算區(qū)域內(nèi)測站的數(shù)量及空間分布有關(guān).因此，以下通過棋盤測試分析現(xiàn)有72個GPS測站分布對研究區(qū)內(nèi)陸地水儲量的恢復(fù)能力.將研究區(qū)域劃分為0.5°×0.5°規(guī)則格網(wǎng)，棋盤信號采用2°×2°的EWH表示，如圖7a所示，其中白色和棕色棋盤分別表示0 mm和300 mm的EWH信號.利用圖7a中的棋盤信號模擬了72個GPS測站的垂直位移，然后將其反演得到如圖7b所示的EWH空間分布.從圖7b可以看出，川東盆地和川南地區(qū)的GPS測站相對密集，其反演結(jié)果對原始信號的恢復(fù)較好；而川西北地區(qū)的GPS測站分布相對稀疏，其反演結(jié)果對原始信號的恢復(fù)較差.總體而言，本文采用的72個GPS測站數(shù)據(jù)，它們之間的平均測站距離約為169.6 km，可以較好地恢復(fù)四川省內(nèi)約2°×2°空間分辨率的TWSC信號.

圖7 72個GPS測站恢復(fù)EWH信號的棋盤測試(a) 原始信號； (b) 反演結(jié)果.Fig.7 Checkerboard test of EWH signal recovered from 72 GPS stations(a) Original signal; (b) Inversion result.

3.3 與GRACE/GFO反演結(jié)果對比

圖8為官方機構(gòu)(JPL、CSR和GSFC)發(fā)布的GRACE/GFO Mascon模型與GPS反演的TWSC在2010年12月至2021年2月期間的周年振幅空間分布.其中，圖8a—c中的周年振幅為同時擬合GRACE和GFO數(shù)據(jù)的結(jié)果，圖8d中GPS反演的TWSC采用了截斷階數(shù)為44階的TSVD-Tikhonov組合正則化方法，并采用GCV方法確定最優(yōu)正則化參數(shù).

圖8 GRACE/GFO Mascon模型(JPL、CSR和GSFC)與GPS垂直位移反演的TWSC周年振幅空間分布Fig.8 Spatial patterns of annual amplitudes for TWSC derived from GRACE/GFO Mascons (JPL, CSR, and GSFC) and GPS vertical displacements

從圖8可以看出，JPL和CSR Mascon模型得到的TWSC周年振幅的空間分布較為一致，其在川東盆地的信號較強，川西高原的信號較弱，并且JPL比CSR反演結(jié)果的信號更強；而GSFC Mascon模型與GPS反演結(jié)果的周年振幅空間分布在總體上更為接近，它們在川西南山地的信號較強，而在東北部地區(qū)的信號較弱，并且GPS比GSFC反演結(jié)果的信號更強.總體上看，GPS反演結(jié)果與GRACE/GFO Mascon模型的周年振幅保持了較好的一致性，但GPS反演的TWSC信號振幅更強，這是由于基于地基觀測的GPS垂直形變數(shù)據(jù)對局部水質(zhì)量負(fù)荷變化更為敏感，而GRACE/GFO反演結(jié)果的空間分辨率有限(～300 km)，受到信號泄漏及混頻效應(yīng)的影響，其對大尺度的TWSC估計更為準(zhǔn)確.同時，GPS和GRACE/GFO觀測數(shù)據(jù)誤差來源及數(shù)據(jù)處理策略的不同也是導(dǎo)致兩者反演結(jié)果具有差異的一個重要原因.另外，與Jiang 等(2021b)利用GPS反演的四川省TWSC的周年振幅相比，本文反演結(jié)果的空間分布特征與其總體接近，但振幅量級相對偏小.引起該差異的主要原因可能有兩方面：一方面是由于不同的GPS解算軟件和解算策略會引起反演結(jié)果的差異，另一方面是本文反演TWSC所采用的四川省內(nèi)的GPS測站數(shù)有限且分布更為稀疏，進(jìn)而影響了TWSC信號的恢復(fù)性能.

圖9為GPS反演的TWSC、GRACE/GFO Mascon模型(JPL、CSR和GSFC)及GTCH方法融合的Mascon模型(定義為GTCH-TWSC)在2010年12月至2021年2月之間的時間序列變化.可以看出，GPS和GRACE/GFO反演結(jié)果總體保持較好的一致性，并且都很好地反映了四川省TWSC明顯的季節(jié)性變化特征.此外，利用GPS反演的TWSC時間序列的振幅總體大于GRACE/GFO反演結(jié)果，并且能夠很好地填補GRACE和GFO任務(wù)間隔期的數(shù)據(jù)空白.需要指出的是，圖9還顯示出GPS反演結(jié)果在2011、2012和2017等年份的年初或年終存在異常信號，且比GRACE/GFO反演結(jié)果更敏感，這可能是由于GPS時間序列中仍然含有未被模型化的虛假信號所引起，如交點年誤差(Fu et al., 2015; 何思源等，2018).

圖9 GRACE/GFO Mascon模型(JPL、CSR和GSFC)及其融合模型(GTCH-TWSC)與GPS垂直位移反演的TWSC時間序列比較Fig.9 Comparisons of TWSC time series derived from GRACE/GFO Mascons (JPL, CSR, and GSFC) and their fusion model (GTCH-TWSC), and GPS vertical displacements

由于GPS垂直位移時間序列在預(yù)處理中扣除了長期趨勢項，為了與其保持一致性，將GRACE/GFO反演結(jié)果同樣扣除趨勢項，然后統(tǒng)計兩者的相關(guān)性和差值STD.表1為去掉趨勢后GRACE/GFO Mascon模型(JPL、CSR和GSFC)及其GTCH融合模型(GTCH-TWSC)與GPS反演的TWSC時間序列的相關(guān)系數(shù)和差值STD統(tǒng)計.對比表1中三個Mascon模型的統(tǒng)計結(jié)果可知，JPL與GPS反演結(jié)果的相關(guān)性最好且STD最小，GSFC次之，而CSR對應(yīng)的相關(guān)性和STD略差.而GTCH融合模型GTCH-TWS與GPS反演結(jié)果的相關(guān)性和差值STD均為最好，這說明通過GTCH方法估計的誤差方差對不同模型進(jìn)行加權(quán)得到的融合模型GTCH-TWS更加可靠，其與GPS反演結(jié)果的一致性也更好.

表1 GRACE/GFO Mascon模型(JPL、CSR和GSFC)及其融合模型(GTCH-TWSC)與GPS反演的TWSC時間序列去除線性趨勢后的相關(guān)系數(shù)和差值STD統(tǒng)計Table 1 Correlation coefficients and STD of the differences between TWSC time series derived from GPS and GRACE/GFO Mascons (JPL, CSR, and GSFC) and their fusion model (GTCH-TWSC) after removing linear trends

表2給出了GPS反演的TWSC、GRACE/GFO Mascon模型(JPL、CSR和GSFC)及其融合模型GTCH-TWSC分別在2010年12月—2021年2月、2010年12月—2017年6月和2018年6月(GPS為2017年7月)—2021年2月等三個時間段對應(yīng)時間序列的周年振幅、相位和線性趨勢統(tǒng)計.由表2可知，GPS反演的TWSC時間序列在三個時段的周年振幅均大于GRACE/GFO結(jié)果，這是因為GPS反演的空間分辨率更高，并且地基GPS對局部范圍的短波信號更為敏感.對于周年相位，GPS和GRACE/GFO反演結(jié)果的一致性較好，不同時段內(nèi)各種結(jié)果的相位差異小于半個月.對于線性趨勢，各種Mascon模型及其融合模型均呈現(xiàn)增長趨勢，其中JPL反演結(jié)果的增長趨勢最為明顯，CSR與GSFC反演結(jié)果的增長趨勢相當(dāng).相反，GPS反演結(jié)果在以上三個時段均呈現(xiàn)下降趨勢.其主要原因是，本文對GPS垂直位移時間序列做了去趨勢處理(用于扣除GIA和板塊構(gòu)造運動導(dǎo)致的長期趨勢)，進(jìn)而影響了GPS反演結(jié)果對應(yīng)線性趨勢的準(zhǔn)確性.需要指出的是，四川省TWSC主要為季節(jié)性變化，其趨勢變化較小，因而各種結(jié)果估計趨勢的不確定度也較大.

表2 GPS、GRACE/GFO Mascon模型(JPL、CSR和GSFC)及其融合模型(GTCH-TWSC)估計的TWSC時間序列的周年振幅、相位和線性趨勢統(tǒng)計Table 2 Annual amplitudes, phases, and linear trends of TWSC time series estimated by GPS, GRACE/GFO Mascons (JPL, CSR, and GSFC) and their fusion model (GTCH-TWSC)

3.4 水文氣象數(shù)據(jù)驗證結(jié)果

圖10為采用不同類型水文氣象資料計算的四川省2011—2021年的月平均降水、蒸散發(fā)和徑流時間序列.其中，圖10a為不同數(shù)據(jù)源計算的月平均降水變化時間序列，可以看出不同類型降水?dāng)?shù)據(jù)的一致性較好，其周年振幅約為200 mm/month，但在降水量達(dá)到峰值的月份仍有一定的差異(如2012、2018和2020年的7—9月).圖10b為不同數(shù)據(jù)源計算的月平均蒸散發(fā)變化時間序列，其周年振幅約為100 mm/month，但不同類型的蒸散發(fā)數(shù)據(jù)在相位和振幅上都存在較大的差異.圖10c為GLDAS三個陸面模型(Noah、CLSM和VIC)計算的月平均徑流變化時間序列，相比于該區(qū)域的月平均降水和蒸散發(fā)變化，其月平均徑流變化的量級略小，振幅在100 mm/month以內(nèi)，并且三個陸面模型計算的徑流變化之間也存在較大的差異.

由于各種類型水文氣象資料的精度和可靠性不同，為了降低不同數(shù)據(jù)源計算的降水、蒸散發(fā)和徑流變化的差異影響，本文采用GTCH方法對各類數(shù)據(jù)進(jìn)行融合，得到更為可靠的降水、蒸散發(fā)和徑流數(shù)據(jù)，以便更好地對GPS和GRACE/GFO反演的dTWSC/dt進(jìn)行驗證.

圖11a是利用圖10a—c中不同類型的月平均降水(P)、蒸散發(fā)(ET)和徑流(R)時間序列，采用GTCH方法融合得到的降水(GTCH-P)、蒸散發(fā)(GTCH-ET)和徑流(GTCH-R)數(shù)據(jù).圖11b是利用融合后的GTCH-P、GTCH-ET和GTCH-R由水量平衡方程(式(11))計算的dTWSC/dt序列(記為PER-dS/dt)，以及GPS反演的TWSC和GRACE/GFO融合模型GTCH-TWSC由式(12)計算的dTWSC/dt序列(分別記為GPS-dS/dt和GRACE/GFO-dS/dt).其中，PER-dS/dt序列采用了式(13)進(jìn)行濾波處理，以保持其與GPS-dS/dt和GRACE/GFO-dS/dt序列的時間同步.從圖11b可看出，PER-dS/dt序列與GPS-dS/dt和GRACE/GFO-dS/dt序列的季節(jié)性變化具有很好的一致性，但GPS-dS/dt序列的振幅略大于PER-dS/dt和GRACE/GFO-dS/dt序列，并且GPS-dS/dt和GRACE/GFO-dS/dt序列的逐月變化差異較大.

圖10 不同類型水文氣象數(shù)據(jù)計算的四川省月平均降水、蒸散發(fā)和徑流變化(2011—2021年)Fig.10 Monthly average precipitation, evapotranspiration, and runoff changes of Sichuan Province derived from different hydrometeorology datasets (from 2011 to 2021)

圖11 (a)GTCH方法融合后的降水、蒸散發(fā)和徑流月平均時間序列； (b)和(c) 分別為平滑前后的PER-dS/dt、GPS-dS/dt和GRACE/GFO-dS/dt時間序列Fig.11 (a) Monthly average precipitation, evapotranspiration, and runoff time series fused by GTCH method; (b) and (c) are the PER-dS/dt, GPS-dS/dt, and GRACE/GFO-dS/dt time series before and after smoothing

為了避免各類數(shù)據(jù)中高頻噪聲的影響，并更為直觀的分析以上三種dTWSC/dt時間序列的季節(jié)性變化特征，以及考慮到時間序列中可能存在的半周年信號，分別對PER-dS/dt、GPS-dS/dt和GRACE/GFO-dS/dt序列做5個月的滑動平均處理(Jiang et al., 2017)，其結(jié)果如圖11c所示.可以看出，平滑后的PER-dS/dt、GPS-dS/dt和GRACE/GFO-dS/dt序列具有更好的一致性，但平滑后的GPS-dS/dt序列的振幅略大于平滑后的PER-dS/dt和GRACE/GFO-dS/dt序列.此外，PER-dS/dt序列相比GPS-dS/dt和GRACE/GFO-dS/dt序列有明顯的相位滯后(統(tǒng)計結(jié)果見表3).

表3給出了GTCH方法融合的降水、蒸散發(fā)和徑流月平均序列，GPS反演的TWSC序列(GPS-TWSC)和GRACE/GFO融合模型GTCH-TWSC(GRACE/GFO-TWSC)，以及PER-dS/dt、GPS-dS/dt和GRACE/GFO-dS/dt序列平滑前后計算的周年振幅和相位統(tǒng)計值.圖12給出了平滑前后的GPS-dS/dt和GRACE/GFO-dS/dt序列與PER-dS/dt序列之間的散點圖分布及相關(guān)系數(shù)和差值STD統(tǒng)計.

表3 不同類型時間序列數(shù)據(jù)的周年振幅和相位統(tǒng)計Table 3 Statistics of annual amplitudes and phases for the different types of time series datasets

由表3可知，GPS-TWSC和GRACE/GFO-TWSC序列比降水、蒸散發(fā)和徑流數(shù)據(jù)的相位滯后了約1～2個月，這是因為TWSC是降水、蒸散發(fā)和徑流變化的綜合反映，這三者需要經(jīng)過復(fù)雜的水動力學(xué)過程及時間變化才能轉(zhuǎn)化為TWSC.比較不同的dTWSC/dt序列可以看出，平滑前后GPS-dS/dt的振幅均為最大，PER-dS/dt的振幅次之，而GRACE/GFO-dS/dt的振幅最小.此外，平滑前后各種序列的相位幾乎保持不變，但與PER-dS/dt序列相比，GPS-dS/dt與GRACE/GFO-dS/dt序列的相位更為一致.

從圖12可以看出，對于未平滑的結(jié)果，GPS-dS/dt和GRACE/GFO-dS/dt序列與PER-dS/dt序列的相關(guān)系數(shù)分別為0.53和0.73，差值STD分別為43.82 mm和18.37 mm，如圖12a和b所示.而平滑以后(見圖12c和d)，GPS-dS/dt和GRACE/GFO-dS/dt序列與PER-dS/dt序列的相關(guān)系數(shù)分別為0.78和0.87，差值STD分別為15.01 mm和8.44 mm.也就是說，平滑后GPS-dS/dt和GRACE/GFO-dS/dt序列與PER-dS/dt序列的差值STD明顯減小，相關(guān)性有了顯著提高(特別是GPS-dS/dt序列).但與GRACE/GFO-dS/dt序列相比，平滑前后GPS-dS/dt與PER-dS/dt序列之間的差異更大.

圖12 平滑前后GPS-dS/dt和GRACE/GFO-dS/dt序列與PER-dS/dt序列之間的散點圖及相關(guān)系數(shù)(CC)和差值STD比較Fig.12 Comparisons of scatter plot, correlation coefficient (CC), and STD of the differences between GPS-dS/dt, GRACE/GFO-dS/dt and PER-dS/dt series before and after smoothing

分析其原因，一方面考慮到地基GPS垂直形變數(shù)據(jù)對局部水質(zhì)量負(fù)荷變化更敏感，其反演結(jié)果的幅度變化更大，進(jìn)而引起GPS-dS/dt與PER-dS/dt序列之間的差值STD較大.另一方面，不同數(shù)據(jù)源計算的蒸散發(fā)和徑流數(shù)據(jù)之間存在明顯的差異(見圖10b和c)，其GTCH融合結(jié)果的不確定度也相對較大，導(dǎo)致PER-dS/dt序列也存在較大的不確定性.考慮到不同源的降水?dāng)?shù)據(jù)差異較小(見圖10a)，具有更好的可靠性，為此進(jìn)一步分析了GPS和GRACE/GFO反演結(jié)果與GTCH融合的降水序列(GTCH-P)之間的相關(guān)性，結(jié)果顯示GPS反演的TWSC序列與GTCH-P序列經(jīng)5個月滑動平均處理前后的相關(guān)系數(shù)分別為0.55和0.85，而GRACE/GFO融合序列GTCH-TWSC與GTCH-P序列在平滑前后的相關(guān)系數(shù)分別為0.52和0.72，這說明GPS相比GRACE/GFO觀測結(jié)果對降水變化的響應(yīng)更為敏感，再次印證了其對局部水質(zhì)量負(fù)荷變化的感應(yīng)能力更強.以上分析表明，GPS和GRACE/GFO在監(jiān)測TWSC方面各有優(yōu)勢，但由于GPS垂直位移對局部水質(zhì)量負(fù)荷變化更為敏感，因此可作為GRACE/GFO反演區(qū)域TWSC的有益補充.

4 結(jié)論

本文研究了基于TSVD-Tikhonov組合正則化方法求解GPS垂直位移反演TWSC的病態(tài)問題，并通過數(shù)值模擬驗證了其有效性.利用CMONOC發(fā)布的72個GPS測站的垂直位移數(shù)據(jù)反演了四川省的TWSC時間序列，并與官方機構(gòu)(JPL、CSR和GSFC)發(fā)布的GRACE/GFO Mascon模型進(jìn)行對比.同時，利用水文氣象資料由水量平衡方程計算的dTWSC/dt序列對GPS反演結(jié)果進(jìn)行驗證.主要結(jié)論如下：

(1)利用模擬的72個GPS測站的垂直位移數(shù)據(jù)反演四川省TWSC，結(jié)果表明：采用不同的正則化參數(shù)選取策略(RMSE最小準(zhǔn)則、GCV法和L-curve法)，TSVD-Tikhonov正則化反演的TWSC比單獨使用TSVD或Tikhonov正則化的反演結(jié)果具有更高的精度和穩(wěn)定性；TSVD、Tikhonov和TSVD-Tikhonov正則化方法反演2005年1月至12月的TWSC時間序列與原始信號差值的平均STD分別為14.97 mm、7.03 mm和5.04 mm.

(2)采用棋盤測試分析了現(xiàn)有72個GPS測站恢復(fù)四川省TWSC信號的能力，并利用實測的GPS測站垂直位移數(shù)據(jù)反演了四川省2010年12月至2021年2月的TWSC時間序列，結(jié)果表明：現(xiàn)有72個測站可以較好地恢復(fù)四川省內(nèi)2°×2°空間分辨率的TWSC信號；基于TSVD-Tikhonov組合正則化由GPS反演的TWSC與GRACE/GFO Mascon模型(JPL、CSR和GSFC)及其GTCH融合模型的空間分布特征及季節(jié)性變化符合較好，并且GPS反演結(jié)果與GTCH融合模型更為一致；GPS反演的TWSC信號振幅比GRACE/GFO Mascon模型更強，這是由于地基GPS垂直形變觀測數(shù)據(jù)對局部水質(zhì)量負(fù)荷變化更為敏感，而GRACE/GFO反演的空間分辨率約為300 km，受到信號泄漏的影響，其對大尺度的TWSC估計更為準(zhǔn)確.

(3)采用GTCH方法融合不同類型的降水、蒸散發(fā)和徑流數(shù)據(jù)，根據(jù)水量平衡方程計算的dTWSC/dt序列(PER-dS/dt)對GPS反演結(jié)果(GPS-dS/dt)和GRACE/GFO Mascon模型融合結(jié)果(GRACE/GFO-dS/dt)進(jìn)行驗證，結(jié)果顯示：GPS-dS/dt和GRACE/GFO-dS/dt序列與PER-dS/dt序列具有較好的一致性，平滑后GPS-dS/dt和GRACE/GFO-dS/dt序列與PER-dS/dt序列的相關(guān)系數(shù)分別為0.78和0.87，差值STD分別為15.01 mm和8.44 mm；GPS反演的TWSC序列和GRACE/GFO融合序列(GTCH-TWSC)與降水融合序列(GTCH-P)在平滑后的相關(guān)系數(shù)分別為0.85和0.72，說明GPS相比GRACE/GFO對降水變化的響應(yīng)更為敏感.總體上講，GPS和GRACE/GFO反演TWSC各有優(yōu)勢，但GPS觀測對局部水質(zhì)量負(fù)荷變化的感應(yīng)能力更強，因此可作為GRACE/GFO反演區(qū)域TWSC的有益補充.

(4)TSVD-Tikhonov正則化方法有效提高了利用GPS垂直位移反演區(qū)域TWSC的精度和可靠性，但在實際應(yīng)用中如何快速有效地選取最優(yōu)截斷階數(shù)和最優(yōu)正則化參數(shù)還需要深入探討.同時，本文反演四川省TWSC采用的GPS測站分布較為稀疏，導(dǎo)致TWSC信號的恢復(fù)能力有限，后續(xù)可采用更為密集的GPS測站數(shù)據(jù)進(jìn)行反演與驗證.此外，有效分離GPS垂直位移時間序列中非水文負(fù)荷因素(如人類活動、GIA或板塊構(gòu)造運動等)導(dǎo)致的形變，并準(zhǔn)確估計GPS反演的TWSC的長期趨勢仍是一項具有挑戰(zhàn)性的課題.

致謝感謝國家地震科學(xué)數(shù)據(jù)中心提供的CMONOC GPS垂直位移數(shù)據(jù)，JPL、CSR和GSFC提供的GRACE/GFO Mascon模型，NMIC和TRMM提供的降水?dāng)?shù)據(jù)，GLEAM提供的蒸散發(fā)數(shù)據(jù)，NCAS提供的降水和蒸散發(fā)數(shù)據(jù)，以及GLDAS提供的降水、蒸散發(fā)和徑流數(shù)據(jù).感謝三位匿名審稿專家對本文提出的寶貴修改意見和建議.