利用GNSS垂直位移研究長江流域陸地水儲量變化的時(shí)空分布特征

湯苗， 鐘萍*， 姜中山， 劉晚純， 成帥， 楊興海

1 西南交通大學(xué)地球科學(xué)與環(huán)境工程學(xué)院， 成都 611756 2 湖南省交通規(guī)劃勘察設(shè)計(jì)院有限公司， 長沙 410219

0 引言

陸地水儲量(Terrestrial Water Storage, TWS)通常定義為存儲在地表、土壤、地下、冰雪、樹冠以及生物當(dāng)中水分的總和(Famiglietti, 2004; Wahr et al., 2004).與大氣、海洋負(fù)荷不同，陸地水儲量往往因缺少可靠的觀測數(shù)據(jù)而難以量化(Dai and Trenberth, 2002).因此，傳統(tǒng)的陸地水儲量研究一般通過氣象數(shù)據(jù)結(jié)合地表模型建立的水文模型進(jìn)行，但這些數(shù)學(xué)模型未對所有水成分完整建模，并且低估了人為影響以及氣候變遷導(dǎo)致的陸地水異常情況(Scanlon et al., 2018).現(xiàn)代衛(wèi)星大地測量技術(shù)的發(fā)展為陸地水儲量變化研究提供了革命性的工具，例如衛(wèi)星重力觀測技術(shù)、全球衛(wèi)星定位系統(tǒng)等長達(dá)數(shù)十年的連續(xù)觀測為水儲量變化研究提供了寶貴的基礎(chǔ)數(shù)據(jù)(重力場、地表變形數(shù)據(jù))(Tapley et al., 2004; Argus et al., 2014)，這些數(shù)據(jù)能夠反映所有水成分的總體變化情況，彌補(bǔ)了水文學(xué)研究中地表及地下水實(shí)測值不足的缺陷.

重力恢復(fù)與氣候試驗(yàn)GRACE (Gravity Recovery and Climate Experiment)衛(wèi)星通過主動(dòng)測距方式獲取的全球時(shí)變重力場被廣泛用于全球大尺度陸地水儲量變化研究(Tapley et al., 2004; Wahr et al., 2004).近年來，基于GRACE數(shù)據(jù)對長江流域內(nèi)水資源及極端天氣事件的研究也較多，例如Zhang等(2015)采用GRACE數(shù)據(jù)研究了長江流域干旱事件，并探索了陸地水異常與厄爾尼諾的可能聯(lián)系；Zhang等(2016, 2019)利用GRACE數(shù)據(jù)對長江流域總體及子流域水文干旱分別進(jìn)行了評估；Wang等(2020)提出了一種加權(quán)水儲量虧損估計(jì)方法并結(jié)合GRACE數(shù)據(jù)來表征長江流域干旱事件.但GRACE數(shù)據(jù)主要存在以下局限：(1)受限于軌道與運(yùn)行機(jī)制，其時(shí)空分辨率較低(時(shí)間分辨率為1月，空間分辨率～350 km)；(2)高階噪聲與條帶誤差(Swenson and Wahr, 2006)需要進(jìn)行濾波處理，進(jìn)一步損耗和壓制了原始信號；(3)數(shù)據(jù)質(zhì)量不高并伴隨著數(shù)據(jù)缺失，插值等后處理容易引入人為誤差.這些因素導(dǎo)致GRACE衛(wèi)星對高頻和區(qū)域尺度的陸地水儲量變化不敏感(Enzminger et al., 2018; Knappe et al., 2019)，難以滿足區(qū)域水資源管理利用、自然災(zāi)害防范的高時(shí)空分辨率需求(Feng et al., 2018).

全球定位系統(tǒng)(Global Positioning System)能夠?qū)崟r(shí)監(jiān)測固體地球表面由于質(zhì)量遷移產(chǎn)生的彈性形變，為實(shí)時(shí)監(jiān)測陸地水儲量變化提供了一種全新的解決方案(Argus et al., 2014; Xue et al., 2021).Farrell(1972)提出的質(zhì)量負(fù)荷理論詳細(xì)描述了負(fù)荷變遷與其造成的地表彈性變形之間的關(guān)系，目前根據(jù)質(zhì)量負(fù)荷計(jì)算地表響應(yīng)一般采用負(fù)荷格林函數(shù)法或球諧函數(shù)法(沈迎春等, 2017).鑒于無法獲得全球性的(海洋及陸地)負(fù)荷位移，不能利用傳統(tǒng)球諧函數(shù)將其轉(zhuǎn)化至譜域內(nèi)，大多數(shù)學(xué)者常用格林函數(shù)法對陸地水儲量變化進(jìn)行反演.例如，Argus等(2014)利用GPS位移數(shù)據(jù)結(jié)合負(fù)荷格林函數(shù)方法估算了美國西部的陸地水儲量變化；Borsa等(2014)使用該方法研究了美國西部干旱期間的地下水流失量；Fu等(2015)計(jì)算了美國西北各州的陸地水儲量異常；Johnson等(2017)利用GPS垂向形變時(shí)間序列來約束小地震斷層面上每月的水文負(fù)荷和由此產(chǎn)生的應(yīng)力變化模型，推斷加州地震活動(dòng)受水文循環(huán)的適度調(diào)節(jié)；Enzminger等(2018)則利用該反演模型估算了美國西部山區(qū)的雪水當(dāng)量.利用GPS位移時(shí)間序列結(jié)合負(fù)荷格林函數(shù)方法能夠獲取50～100 km內(nèi)天或天以下的負(fù)荷分布情況(Han and Razeghi, 2017)，明顯優(yōu)于GRACE所能提供的時(shí)空分辨率.

相關(guān)學(xué)者試圖采用Slepian函數(shù)法從譜域中獲取大尺度的陸地水儲量變化(Han and Razeghi, 2017; Jiang et al., 2021a; 成帥等, 2021).Slepian函數(shù)最早用于解決一維信號在時(shí)頻域的最優(yōu)集中問題(Slepian, 1964).為了解決大地測量學(xué)中存在的“極空白”問題，Sneeuw和Van Gelderen(1997)及Albertella等(1999)在球體上構(gòu)建了一組新的Slepian基函數(shù)，使大部分能量最佳地集中在由整個(gè)地球減去極地間隙組成的緯向帶內(nèi).在此基礎(chǔ)上，Simons等(2006)又實(shí)現(xiàn)了球面空間任意區(qū)域的譜集中.后來的學(xué)者基于Slepian函數(shù)對冰蓋質(zhì)量變化進(jìn)行了研究，例如Harig和Simons(2012)利用Slepian方法對GRACE球諧解進(jìn)行譜集中，以此研究格陵蘭島質(zhì)量時(shí)空變化；高春春等(2019)基于此提出了一種Slepian空域反演的方法，估算了南極冰蓋子流域的質(zhì)量變化.國內(nèi)學(xué)者又將這一方法引入到重力場研究領(lǐng)域當(dāng)中，陳石等(2017)利用其解算了中國大陸重力場的球諧模型；韓建成等(2021)基于地面重力觀測建立了華北地區(qū)120階的時(shí)變重力場模型.而Han和Razeghi(2017)利用Slepian基函數(shù)對GPS數(shù)據(jù)進(jìn)行參數(shù)化和反演，獲得了澳大利亞水文和大氣質(zhì)量的變化情況，使Slepian方法逐漸應(yīng)用于陸地水儲量反演；Jiang等(2021a)和成帥等(2021)采用同樣的方法分別對中國地區(qū)與川云渝地區(qū)的陸地水儲量變化進(jìn)行了研究.

長江流域地域遼闊，地形起伏大，同時(shí)受東亞季風(fēng)、印度季風(fēng)影響，區(qū)域內(nèi)降水時(shí)空分布不均，是旱澇災(zāi)害的多發(fā)地.近年來，全球變暖導(dǎo)致了水汽循環(huán)加速，降水的時(shí)空分布和強(qiáng)度受到了極大的影響(Cao et al., 2011)，長江流域極端天氣事件也隨之增加.例如2010年發(fā)生特大洪水(Sun et al., 2017)、2011年春季出現(xiàn)嚴(yán)重干旱等(Zhang et al., 2015; Sun et al., 2018; Wang et al., 2020)，這些自然災(zāi)害給當(dāng)?shù)厮Y源、農(nóng)業(yè)、自然生態(tài)系統(tǒng)和社會造成了重大影響.因此，量化長江流域陸地水儲量變化的時(shí)空分布特征將有助于理解區(qū)域內(nèi)水文循環(huán)模式以及水文與氣象相互作用的機(jī)制(Wang et al., 2020)，對旱澇災(zāi)害的預(yù)防、水資源管理利用及區(qū)域可持續(xù)發(fā)展有著重要意義(Jiang et al., 2021a).

以上研究說明目前迫切需要掌握長江流域水儲量變化及其時(shí)空格局，但關(guān)于利用GNSS研究長江流域陸地水儲量變化還較少.因此，本文利用長江流域GNSS垂直坐標(biāo)時(shí)間序列，采用Slepian函數(shù)法估計(jì)長江流域陸地水儲量變化，并探索了流域內(nèi)水儲量變化的時(shí)空分布模式.將估算結(jié)果與其他大地測量和水文數(shù)據(jù)集進(jìn)行對比分析，驗(yàn)證了長江流域水儲量變化的時(shí)空分布特征，并量化了各子區(qū)域(金沙江流域、傳統(tǒng)上游和長江中下游)不同數(shù)據(jù)集結(jié)果的相關(guān)程度和時(shí)間滯后關(guān)系，同時(shí)對反演模型的精度進(jìn)行了評定.最后利用GNSS估計(jì)的水儲量異常對極端天氣事件進(jìn)行追蹤，并結(jié)合歷史極端天氣事件進(jìn)行驗(yàn)證，為長江流域極端天氣事件的監(jiān)測提供了新的觀測依據(jù).

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

1.1 GNSS數(shù)據(jù)及處理

本文采用的GNSS時(shí)序產(chǎn)品從國家地震科學(xué)數(shù)據(jù)中心-東部形變數(shù)據(jù)服務(wù)中心(http:∥www.eqdsc.com)獲取，該產(chǎn)品利用中國地殼運(yùn)動(dòng)觀測網(wǎng)絡(luò)(Crustal Movement Observation Network of China, CMONOC)獲取的觀測數(shù)據(jù)，采用麻省理工學(xué)院開發(fā)的GAMIT/GLOBK 10.40 (Herring et al., 2018)軟件結(jié)合QOCA軟件解算得到.本文選取了長江流域98個(gè)站點(diǎn)2011年1月—2020年12月單天解坐標(biāo)時(shí)間序列，所選站點(diǎn)空間分布如圖1所示.由于垂向位移對水文質(zhì)量變化最為敏感(Dong et al., 2002; Kraner et al., 2018; Chanard et al., 2018; Panda et al., 2018)，因此本文僅使用垂直位移進(jìn)行反演建模.為了獲取水文負(fù)荷位移變形，我們還需扣除其他非水文的影響因素，例如：線性的構(gòu)造運(yùn)動(dòng)、非潮汐海洋與大氣負(fù)荷變化(姜衛(wèi)平等, 2014)、強(qiáng)震影響或儀器更新等引起的階躍現(xiàn)象、以及熱膨脹變形和冰后回彈效應(yīng)等.原始GNSS垂直位移時(shí)序處理步驟如下：

(1)扣除非潮汐海洋與大氣負(fù)荷位移：位移時(shí)序中包含的非潮汐海洋與大氣負(fù)荷位移均采用德國地學(xué)研究中心(Geo Forschungs Zentrum, GFZ)提供的產(chǎn)品進(jìn)行扣除(http:∥rz-vm115.gfz-potsdam.de:8080/repository)(Dill and Dobslaw, 2013).

(2)改正熱膨脹變形：熱膨脹效應(yīng)對GNSS臺站垂直位移的影響分為兩部分，一是安裝GNSS天線的水泥墩隨溫度變化熱脹冷縮，另外地表溫度變化傳遞到基巖也會引起GNSS臺站發(fā)生垂向位移(閆昊明等, 2010; 姜衛(wèi)平等, 2015).本文從中國氣象局獲取每日地表溫度(http:∥data.cma.cn)，并按閆昊明等(2010)提出的經(jīng)驗(yàn)?zāi)Ｐ凸浪銦崤蛎涀冃?，其中模型參?shù)值采用賈路路等(2018)給出的相關(guān)參數(shù).

(3)提取水文負(fù)荷形變：扣除上述因素引起的垂向位移后，將剩余的位移時(shí)間序列按公式(1)進(jìn)行最小二乘擬合去除長期線性趨勢、偏移量和震后變形(Jiang et al., 2017; 姜衛(wèi)平等, 2018)：

(1)

式中，y(ti)為ti時(shí)刻位置，y0為站點(diǎn)參考位置，v為線性項(xiàng)速率，Aj與Bj分別為周期項(xiàng)對應(yīng)的振幅(j=1時(shí)為年信號，j=2時(shí)為半年信號)，gk為階躍項(xiàng)，Tk為階躍發(fā)生的時(shí)刻，H為階梯函數(shù)，ε(ti)為隨機(jī)過程.

圖1 GNSS站點(diǎn)分布及長江流域的地理范圍藍(lán)線為長江干流，黃線為長江流域邊界，Ⅰ為金沙江流域，Ⅱ?yàn)閭鹘y(tǒng)上游地區(qū)，Ⅲ為長江中下游地區(qū).Fig.1 GNSS station distribution and geographical range of the Yangtze River Basin Blue line represents the mainstream of the Yangtze River, and yellow line is the boundary of the Yangtze River Basin. Regions Ⅰ—Ⅲ are the Jinsha River Basin, the traditional upper reaches, and the middle and lower reaches of the Yangtze River Basin, respectively.

(4)補(bǔ)齊數(shù)據(jù)：本文獲取的站點(diǎn)時(shí)間序列長度均在1000天以上，平均數(shù)據(jù)量達(dá)到3446天，占研究時(shí)間總長度的94%，具有較高的數(shù)據(jù)完整性.為保證具有相同采樣率的連續(xù)時(shí)間序列進(jìn)行時(shí)變反演，缺失的數(shù)據(jù)通過使用GMIS軟件(GNSS Missing Data Interpolation Software)進(jìn)行補(bǔ)齊(Liu et al., 2018).這種動(dòng)態(tài)時(shí)空插值方法采用克里金-卡爾曼濾波(Kriging Kalman Filter, KKF)，在處理基于網(wǎng)絡(luò)坐標(biāo)時(shí)間序列的連續(xù)和隨機(jī)數(shù)據(jù)間斷方面具有優(yōu)異的性能.以部分站點(diǎn)為例，圖2展示了JSYC、HNLY、QHBM、CQCS、SCLT、YNXP六個(gè)站點(diǎn)GNSS水文負(fù)荷位移時(shí)序經(jīng)過補(bǔ)齊前后的結(jié)果.

圖2 GMIS軟件對水文負(fù)荷位移時(shí)序進(jìn)行補(bǔ)齊后的濾波結(jié)果OBS表示原始時(shí)間序列，INP表示補(bǔ)齊后的結(jié)果.Fig.2 The filtered results of hydrological loading displacement time series that filled by GMIS software OBS represents the original time series and INP represents the result after filling.

1.2 GRACE數(shù)據(jù)

本文采用的GRACE數(shù)據(jù)為德克薩斯大學(xué)空間研究中心(Center for Space Research, CSR)發(fā)布的RL06 Mascon完全校準(zhǔn)模型第二版本(http:∥www2.csr.utexas.edu/grace/RL06_mascons.html)，研究時(shí)段為2011年1月至2020年12月，期間共有31個(gè)月的數(shù)據(jù)缺失.本文未對GRACE與GRACE Follow-on的銜接空缺(11個(gè)月)進(jìn)行處理，其余缺失的數(shù)據(jù)采用三次樣條插值進(jìn)行補(bǔ)齊，最后將數(shù)據(jù)重采樣至0.5°×0.5°格網(wǎng).該模型解直接通過GRACE衛(wèi)星Level-1B數(shù)據(jù)進(jìn)行解算，在時(shí)間和空間上采用先驗(yàn)信息進(jìn)行約束，以減少南北條帶誤差和測量誤差(Save et al., 2016)，因此數(shù)據(jù)使用過程中不需要任何濾波或縮放因子，可以省略平滑或去條帶的后處理過程.數(shù)據(jù)發(fā)布前修正了C20、C30項(xiàng)，并對一階地心項(xiàng)進(jìn)行了改正，同時(shí)扣除了冰后回彈效應(yīng)以及非潮汐大氣、海洋負(fù)荷.最后扣除了2004—2009年的平均值，直接發(fā)布陸地水儲量的每月變化值.

1.3 水文數(shù)據(jù)

水文模型采用美國航天航空局提供的全球陸地?cái)?shù)據(jù)同化模型(Global Land Data Assimilation System, GLDAS) Noah2.1同化數(shù)據(jù)(https:∥disc.gsfc.nasa.gov/datasets)，數(shù)據(jù)以0.25°×0.25°的格網(wǎng)形式發(fā)布，時(shí)間分辨率為1個(gè)月.該模型依賴于多個(gè)陸地表面模型和氣象驅(qū)動(dòng)數(shù)據(jù)集的組合，綜合了降水、太陽輻射、氣溫和其他氣象因素的影響(Rodell et al., 2004)，最終提供從地表到2 m深度的土壤水估計(jì)值、雪水當(dāng)量和冠層水分量.本文對各水成分求和后經(jīng)過去平均處理獲得GLDAS水儲量變化，同樣將其重采樣至0.5°×0.5°格網(wǎng).

降水作為陸地水的來源，是直接影響陸地水儲量變化的一個(gè)重要因素，因此本文結(jié)合中國氣象局氣象數(shù)據(jù)服務(wù)中心的逐月降水資料嘗試探索降水量與陸地水儲量變化的關(guān)系.該數(shù)據(jù)集整合了中國大陸均勻分布的2472個(gè)地面氣象臺站的降水資料，處理后經(jīng)交叉驗(yàn)證對質(zhì)量狀況進(jìn)行檢驗(yàn)，最后通過薄板樣條的空間插值方法給定中國大陸所有格網(wǎng)點(diǎn)的降水值，并以0.5°×0.5°的格網(wǎng)形式發(fā)布.

2 原理與方法

2.1 Slepian基函數(shù)

理想情況下，全球性的物理場量在譜域內(nèi)具有無限帶寬(Mitrovica et al., 1994;朱廣彬等, 2012)，以球諧函數(shù)作為基函數(shù)可以恢復(fù)球面上任意位置或任意空間尺度下的物理信號.由于參數(shù)估計(jì)的有限性，估計(jì)量具有帶限特征，常用截?cái)嗟侥畴A次的球諧位系數(shù)結(jié)合基函數(shù)來表征這個(gè)物理場量.在譜域內(nèi)任意時(shí)變的球面信號σ(θ,λ,t)可近似表達(dá)為：

(2)

式中，a為地球半徑，l,m分別為展開的階數(shù)和次數(shù)，θ,λ分別為經(jīng)度和余緯，cl m為l階m次面球諧系數(shù)，Yl m為面球諧函數(shù)，L為帶限寬度.

Slepian函數(shù)法能夠構(gòu)建既在全球又在感興趣區(qū)域正交的一組帶限基函數(shù)(Simons and Dahlen, 2006).對于局部而非全球的信號而言，用這組新的基函數(shù)的線性組合能夠近似地表征該信號，使信號能量最佳地集中在研究區(qū)域內(nèi)，物理場量σ可通過新的基函數(shù)表示為：

(3)

式中g(shù)α表示第α項(xiàng)Slepian基函數(shù)，cα表示第α項(xiàng)相應(yīng)的Slepian系數(shù).其中g(shù)α通過如下定義：

(4)

式中g(shù)α,l m表示第α項(xiàng)Slepian基函數(shù)的l階m次膨脹系數(shù)，其本質(zhì)是一組能將原始球諧函數(shù)轉(zhuǎn)化為雙正交(既在全球又在目標(biāo)區(qū)域正交)的Slepian基函數(shù)的線性組合系數(shù).膨脹系數(shù)的最優(yōu)解則當(dāng)基函數(shù)gα在給定封閉區(qū)域R內(nèi)的能量與在整個(gè)地球Ω的能量之比達(dá)到最大值時(shí)取得，即：

(5)

式中，γα表示能量集中度(0≤γα≤1)，γα趨近于1的基函數(shù)能量集中在R內(nèi)，而趨近于0的基函數(shù)能量集中在R外.(5)式是Slepian問題在二維球面內(nèi)的表達(dá)，引入核矩陣D對該問題進(jìn)行求解：

(6)

Dgα,l m=γαgα,l m.

(7)

矩陣中每一項(xiàng)為各階次兩個(gè)球諧函數(shù)在區(qū)域R上的乘積積分，對D矩陣求解特征向量和特征值分別得到gα,l m與γα(式(7)).D矩陣的構(gòu)建由區(qū)域R及帶限階數(shù)L共同確定，D矩陣確定的同時(shí)膨脹系數(shù)gα,l m及能量集中度γα被唯一確定，然后利用(4)式便可確定Slepian基函數(shù).

圖3a給出了Slepian基函數(shù)的能量空間分布，圖中黃色實(shí)線表示長江流域邊界，黑色虛線表示研究區(qū)擴(kuò)充3°后的邊界擴(kuò)充線(Harig and Simons, 2012; 成帥等, 2021).由圖3a可以看出，Slepian基函數(shù)能量通常從中心向外平滑擴(kuò)散，但對邊緣信號響應(yīng)不敏感，邊界擴(kuò)充在一定程度上能夠削弱這一影響.此外，站點(diǎn)水文負(fù)荷位移是站點(diǎn)附近區(qū)域內(nèi)所有水文質(zhì)量負(fù)荷變遷的綜合體現(xiàn)，進(jìn)行邊界擴(kuò)充可以更加準(zhǔn)確地約束研究區(qū)邊緣水文質(zhì)量變化.

圖3 基函數(shù)能量空間分布及能量集中度變化黃色實(shí)線為長江流域邊界，黑色虛線為邊界擴(kuò)充線.α表示Slepian基函數(shù)按能量大小排序后的序號，γ表示能量集中度.Fig.3 Energy spatial distribution and energy concentration change of basis function Solid yellow line is the boundary of the Yangtze River Basin and dashed black line is the expansion line of the boundary. α represents the serial number of the Slepian basis function sorted by energy, and γ represents the degree of energy concentration.

通常情況下，加入運(yùn)算的Slepian基函數(shù)個(gè)數(shù)由香農(nóng)數(shù)N=(L+1)2A/4π確定.香農(nóng)數(shù)作為局部譜優(yōu)化的度量，定義為能量最佳集中在單位球上面積為A/4π的空間區(qū)域內(nèi)最優(yōu)基函數(shù)個(gè)數(shù)(Kennedy and Sadeghi, 2013; von Hippel and Harig, 2019).但在實(shí)際情況中，采用的基函數(shù)個(gè)數(shù)應(yīng)小于GNSS站點(diǎn)數(shù)，否則求解Slepian系數(shù)會出現(xiàn)欠定情況，同時(shí)也要防止基函數(shù)太少和站點(diǎn)過多導(dǎo)致能量損失過多和模型過擬合的現(xiàn)象.本文采用的最大階數(shù)為60階，計(jì)算的香農(nóng)數(shù)值為34，第34個(gè)基函數(shù)對應(yīng)的能量集中度為0.51(如圖3b).結(jié)合香農(nóng)數(shù)與所選站點(diǎn)個(gè)數(shù)，并考慮保留足夠的多余觀測，最終選取前34個(gè)基函數(shù)對研究區(qū)域內(nèi)水文信號進(jìn)行恢復(fù)，其總能量恢復(fù)比例為87.8%.

2.2 反演原理

GNSS垂向位移時(shí)間序列經(jīng)1.1節(jié)步驟處理后，將區(qū)域內(nèi)GNSS站點(diǎn)垂直位移時(shí)間序列作為位移信號，結(jié)合(4)式求解的基函數(shù)并按(8)式通過最小二乘求解位移Slepian系數(shù).

(8)

式中，d,n,J分別表示歷元數(shù)、站點(diǎn)數(shù)和所選基函數(shù)個(gè)數(shù)，u表示水文負(fù)荷位移信號，Uα(td)為td時(shí)刻第α個(gè)位移Slepian系數(shù)，gi(θn,λn)為第n個(gè)站點(diǎn)第i個(gè)Slepian基函數(shù).

根據(jù)質(zhì)量負(fù)荷理論，質(zhì)量負(fù)荷和地表響應(yīng)可通過負(fù)荷勒夫數(shù)在譜域內(nèi)進(jìn)行線性轉(zhuǎn)化(Farrell, 1972)，負(fù)荷位移Slepian系數(shù)與陸地水儲量Slepian系數(shù)在局部地區(qū)譜域內(nèi)的轉(zhuǎn)換關(guān)系如下：

(9)

式中，Sβ(t)表示陸地水儲量Slepian系數(shù)，ρw表示水密度(1000 kg·m-3)，ρe為地球平均密度(5517 kg·m-3)，h′l為l階垂向負(fù)荷勒夫數(shù)(Wang et al., 2012)，兩者的詳細(xì)轉(zhuǎn)化過程可參見Han和Razeghi (2017).

最后，將求解的陸地水儲量Slepian系數(shù)利用Slepian基函數(shù)轉(zhuǎn)化到空域內(nèi)，并按能量集中度加權(quán)即可求解研究區(qū)域內(nèi)陸地水儲量變化.

(10)

3 結(jié)果與分析

3.1 垂向運(yùn)動(dòng)的周年變化

準(zhǔn)確提取陸地水變遷引起的形變位移是有效約束區(qū)域內(nèi)水儲量變化的關(guān)鍵，本節(jié)對非潮汐海洋與大氣負(fù)荷位移以及熱膨脹變形的周年振幅進(jìn)行了詳細(xì)說明，并對水文負(fù)荷引起的地殼垂向季節(jié)性振蕩進(jìn)行了分析.非潮汐大氣、海洋負(fù)荷位移以及熱膨脹變形的周年振幅空間變化模式與地理區(qū)域具有很強(qiáng)的相關(guān)性.如圖4a所示，青藏高原東部地區(qū)非潮汐大氣負(fù)荷位移周年振幅為1～3 mm，四川盆地為3～4 mm，東部平原地區(qū)～5 mm.從沿海到內(nèi)陸，非潮汐海洋負(fù)荷位移周年振幅表現(xiàn)出明顯的遞減趨勢(圖4b)，近海站點(diǎn)的非潮汐海洋負(fù)荷位移周年振幅最高可達(dá)1.2 mm，但遠(yuǎn)海站點(diǎn)通常小于0.6 mm，空間分布特征與van Dam等(2012)的結(jié)果較為一致.熱膨脹變形的周年振幅表現(xiàn)出更為復(fù)雜的分布模式，除海拔因素以外，與緯度也具有強(qiáng)相關(guān)性(圖4c).總的來說，非潮汐大氣、海洋負(fù)荷以及熱膨脹效應(yīng)引起的最大周年振幅分別可達(dá)5 mm、1.2 mm和1.2 mm，對站點(diǎn)垂向位移的貢獻(xiàn)不容忽視.

圖4 非潮汐大氣負(fù)荷、非潮汐海洋負(fù)荷及熱膨脹效應(yīng)引起的垂向位移周年振幅站點(diǎn)顏色表示位移的周年振幅大小. (a)—(c)分別表示非潮汐大氣負(fù)荷位移、非潮汐海洋負(fù)荷位移以及熱膨脹變形的周年振幅.Fig.4 Annual amplitudes of vertical displacement caused by non-tidal atmosphere loading, non-tidal ocean loading, and thermal expansionThe color of site represents the annual amplitude of the displacement. (a)—(c) represent the annual amplitudes of non-tidal atmospheric loading displacement, non-tidal ocean loading displacement, and thermal expansion deformation, respectively.

為了驗(yàn)證提取的GNSS站點(diǎn)水文負(fù)荷位移的準(zhǔn)確性，我們對GNSS及GRACE數(shù)據(jù)獲取的負(fù)荷位移周年振幅進(jìn)行了比較.圖5展示了由全球Mascon模型解結(jié)合負(fù)荷格林函數(shù)正演得到的站點(diǎn)水文負(fù)荷位移的周年振幅GRACE-VCD(Vertical Crustal Deformation, VCD)以及提取的GNSS水文負(fù)荷位移周年振幅GNSS-VCD.從圖5a與圖5b可以看出，GNSS與GRACE獲取的站點(diǎn)水文負(fù)荷位移周年振幅在空間分布上較為相似，但兩者在量級上具有一定的差異.例如，西南地區(qū)GNSS-VCD普遍大于GRACE-VCD，其中GNSS-VCD最大可達(dá)9.7 mm，而GRACE-VCD僅有7.5 mm.此外，JXHK站點(diǎn)GNSS-VCD為8.3 mm，其附近站點(diǎn)GNSS-VCD也相對較大，但是在GRACE-VCD的分布中卻沒有體現(xiàn).這可能與GRACE的低時(shí)空分辨率有關(guān)，另外，GNSS負(fù)荷位移時(shí)序中包含的其他未剔除干凈的地球物理信號也可能會導(dǎo)致這一現(xiàn)象.圖5c展示了GRACE-VCD與GNSS-VCD兩組數(shù)據(jù)之間的相似性，散點(diǎn)坐標(biāo)由各站點(diǎn)GNSS和GRACE獲取的周年振幅確定.所有散點(diǎn)經(jīng)過一次多項(xiàng)式擬合得到的直線斜率為0.54，其中68%的站點(diǎn)GNSS-VCD大于GRACE-VCD，這一結(jié)果說明在長江流域內(nèi)GNSS獲取的水文負(fù)荷位移大于GRACE，其原因可以分為兩個(gè)方面.一方面，從GNSS時(shí)間序列中提取的水文負(fù)荷位移信號可能包含其他非水文信號的影響，例如交年項(xiàng)誤差、非線性的地殼構(gòu)造運(yùn)動(dòng)和觀測噪聲等.Fu等(2015)認(rèn)為GNSS偏大的結(jié)果可能是因?yàn)榻荒觏?xiàng)誤差所引起，并指出交年項(xiàng)振幅占周年項(xiàng)振幅的平均比例為21.24%±13.56%；另一方面，由于空間分辨率的原因，GRACE Mascon模型解中仍然存在一定的泄漏誤差，局部地區(qū)的Mascon結(jié)果不可避免會受到壓制.

圖5 GNSS與GRACE獲取的站點(diǎn)水文負(fù)荷位移周年振幅及其比較Fig.5 Annual amplitudes of hydrological loading displacements obtained by GNSS and GRACE and their comparison

為進(jìn)一步揭示GNSS垂向位移周年相位包含的信息，厘清長江流域在一周年內(nèi)水儲量變化的演變過程，本文結(jié)合降水的周年相位對GNSS站點(diǎn)反映的同期平均周年相位進(jìn)行了討論.圖6展示的GNSS站點(diǎn)水文負(fù)荷變形的周年相位與降水周年相位具有相似的空間分布.從圖6b可以看出，6月由于南北暖冷空氣相遇，形成梅雨(Ye et al., 2019)，長江中下游受梅雨影響率先達(dá)到降水峰值；隨著亞洲夏季風(fēng)的推進(jìn)，其余地區(qū)相繼在7月內(nèi)達(dá)到降水峰值，云南地區(qū)氣候主要受西南季風(fēng)支配，降水峰值稍晚到達(dá).圖6a展示了周年內(nèi)水文負(fù)荷引起的地殼垂向振蕩首次達(dá)到低谷的時(shí)間(即水文負(fù)荷變化首次達(dá)到峰值的時(shí)間)，相比降水周年相位有明顯的延遲.結(jié)合GNSS周年相位分布情況，推測各區(qū)域水文負(fù)荷變化首次達(dá)到峰值的時(shí)間具有以下先后順序：首先，6至7月長江中下游率先達(dá)到水文負(fù)荷變化峰值；結(jié)合降水周年相位結(jié)果，長江流域除中下游以外地區(qū)達(dá)到水文變化峰值的時(shí)間相近，但在8月中旬左右，GNSS結(jié)果顯示青藏高原東邊緣以上區(qū)域提前達(dá)到峰值，這可能是青藏高原東部冰川由于升溫融化提供大量水分導(dǎo)致；最后，9月初至10月初，四川盆地及西南地區(qū)相繼達(dá)到峰值(Jiang et al., 2021b)，相較于降水晚2～3個(gè)月.

圖6 水文負(fù)荷變形及降水在各GNSS臺站的周年相位Fig.6 Annual phases of hydrological loading deformation and precipitation at GNSS stations

3.2 反演模型的不確定性

邊界擴(kuò)充是反演模型提升邊緣信號響應(yīng)、充分約束邊緣信號的重要手段，本節(jié)對不同擴(kuò)充邊界的反演模型估算結(jié)果進(jìn)行了討論，并通過模擬實(shí)驗(yàn)對反演模型的截?cái)嗾`差進(jìn)行了探究.將2011—2020年全球GRACE Mascon模型解作為模擬值，圖7a展示了同期長江流域水文負(fù)荷變化周年振幅的空間分布.采用負(fù)荷格林函數(shù)方法獲得98個(gè)站點(diǎn)的水文負(fù)荷位移時(shí)間序列，再使用Slepian函數(shù)方法結(jié)合不同擴(kuò)充邊界進(jìn)行反演，恢復(fù)了不同緩沖區(qū)下的水文負(fù)荷變化信號(圖7b—7f).由圖7可見，擴(kuò)充邊界為0°或1°時(shí)，信號恢復(fù)的結(jié)果與模擬結(jié)果差異較大，且邊界處信號響應(yīng)不明顯，而當(dāng)擴(kuò)充邊界達(dá)到3°時(shí)，信號恢復(fù)結(jié)果趨于穩(wěn)定.通過穩(wěn)定的恢復(fù)結(jié)果與模擬結(jié)果對比，兩者具有相似的空間分布，但西南邊界處與長江中下游東部平原的恢復(fù)信號小于模擬信號，這可能與反演模型的截?cái)嗾`差有關(guān)，即選取基函數(shù)的過程類似于抑制高頻信號，導(dǎo)致局部地區(qū)恢復(fù)結(jié)果更加平滑，造成區(qū)域內(nèi)部的信號泄漏.實(shí)驗(yàn)結(jié)果表明，當(dāng)擴(kuò)充邊界為3°時(shí)，研究區(qū)域內(nèi)水文負(fù)荷變化能夠得到有效約束，即300 km左右空間分辨率下3°以外質(zhì)量負(fù)荷的影響幾乎可以忽略，同時(shí)該模擬實(shí)驗(yàn)也驗(yàn)證了反演模型的有效性.

圖7 不同邊界擴(kuò)充條件下的水儲量變化空間分布Fig.7 Spatial distribution of water storage variation under different boundary extension conditions

利用Slepian基函數(shù)方法恢復(fù)區(qū)域尺度水儲量變化過程中涉及到兩次截?cái)啵旱谝淮谓財(cái)嗍怯捎赟lepian基函數(shù)的帶限特征所引起，本文選取的帶限寬度為60階，大于60階的高頻信號被舍棄；60階Slepian展開應(yīng)共有(60+1)2個(gè)基函數(shù)，根據(jù)香農(nóng)數(shù)選取前34個(gè)基函數(shù)會造成第二次截?cái)?為分析與評估截?cái)嗾`差對估算結(jié)果的影響，將提取的水文負(fù)荷位移u表示為兩個(gè)分量(截?cái)嗪蟮奈灰菩盘杣res與截?cái)嗾`差δu1)(Han and Razeghi, 2017):

u(θ,λ,t)=ures(θ,λ,t)+δu1(θ,λ,t).

(11)

以2011—2020年全球GRACE Mascon模型解為例，首先利用負(fù)荷格林函數(shù)法計(jì)算站點(diǎn)水文負(fù)荷位移，站點(diǎn)位移時(shí)序RMS值(Root Mean Square)代表模擬值(圖8a).然后利用站點(diǎn)位移時(shí)序獲取的負(fù)荷位移Slepian系數(shù)結(jié)合Slepian基函數(shù)再次獲取站點(diǎn)的水文負(fù)荷位移時(shí)序，模擬時(shí)采用3°擴(kuò)充邊界，截?cái)嗾`差用兩次位移時(shí)序RMS變化量表示(圖8b).從圖8b可以看出，站點(diǎn)模擬值大小范圍為1.5～5.5 mm，兩次截?cái)鄬φ军c(diǎn)位移信號的損耗最大～1 mm，其中93%的站點(diǎn)信號截?cái)嗾`差低于0.8 mm (邊界范圍內(nèi)站點(diǎn)截?cái)嗾`差均小于0.8 mm)，截?cái)嗾`差對站點(diǎn)位移的平均占比為8.1%，對位移信號影響較小.本文利用全球Mascon求解的水文負(fù)荷位移能夠反映短波長和長波長的水文負(fù)荷總體變化情況，而西南地區(qū)及長江中下游的邊界擴(kuò)充范圍內(nèi)站點(diǎn)截?cái)嗾`差較大，結(jié)合上文邊界擴(kuò)充實(shí)驗(yàn)結(jié)果，分析其原因可能是反演模型對邊緣信號不敏感(Harig and Simons, 2012)以及未考慮擴(kuò)充邊界外站點(diǎn)位移的約束，這也能部分說明邊界擴(kuò)充的必要性與有效性.就總體而言，反演模型在研究范圍內(nèi)對原始信號的削弱程度較小，證明該模型能夠完成從離散位移信號到區(qū)域水文負(fù)荷的反演工作.

圖8 站點(diǎn)位移的模擬值與站點(diǎn)位移的截?cái)嗾`差(a) 站點(diǎn)顏色表示站點(diǎn)位移時(shí)序的RMS值； (b) 站點(diǎn)顏色表示截?cái)嗲昂驲MS變化值.Fig.8 Simulation value of station displacement and truncation error of station displacement(a) The color of station represents the RMS value of the station displacement time series; (b) The color of station represents the change of RMS before and after truncation.

3.3 水儲量變化的空間模式

本文利用GNSS站點(diǎn)垂向位移時(shí)間序列，首先反演出各歷元的長江流域水儲量變化，再對得到的各格網(wǎng)點(diǎn)水儲量變化時(shí)序進(jìn)行擬合后獲得格網(wǎng)點(diǎn)的水儲量變化周年振幅大小，所有水儲量變化周年振幅均用等效水高(Equivalent Water Height, EWH)表示.圖9a反映了基于GNSS的長江流域陸地水儲量變化的周年振幅空間分布.從總體上看，長江流域水儲量變化在金沙江流域及長江中下游東部平原較大，在太湖及長江流域中部較小.其中金沙江流域陸地水儲量變化最為明顯，區(qū)域內(nèi)GNSS-EWH最大周年振幅可達(dá)～214 mm，并呈現(xiàn)由西南向東北逐漸減小的趨勢.這一現(xiàn)象可能原因有兩個(gè)：(1)西南地區(qū)氣候受亞洲夏季風(fēng)支配，濕季(5—10月)降水豐沛，為流域內(nèi)補(bǔ)充了大量水分；(2)西南地區(qū)植被覆蓋度高，河流眾多，有利于夏季水的存儲.GRACE-EWH(圖9b)與GLDAS-EWH(圖9c)結(jié)果與GNSS-EWH有較為一致的空間分布特征，三者在金沙江流域變化顯著，且GNSS-EWH與GRACE-EWH在長江中下游表現(xiàn)出較大的變化.降水作為陸地水存儲的主要來源，降水量變化在一定程度上能夠作為陸地水儲量變化的依據(jù)，如圖9d所示，長江流域多年平均降水量呈現(xiàn)由東至西逐步遞減的趨勢.但GLDAS在長江中下游的陸地水儲量變化不明顯，區(qū)域內(nèi)最大周年振幅～60 mm，遠(yuǎn)小于GNSS-EWH的～143 mm及GRACE-EWH的～120 mm.長江中下游眾多支流為地表水的存儲提供了條件，因此這種差異很可能是由于GLDAS模型未對地表水成分進(jìn)行建模所引起(Scanlon et al., 2018).太湖地處亞熱帶季風(fēng)氣候區(qū)，豐富的降水為其補(bǔ)充了大量水分，但由于氣候與所處地理位置的影響，區(qū)域內(nèi)蒸散、徑流量相對較大，水儲量收支較為平衡，這可能是太湖區(qū)域陸地水儲量變化較小的原因.

圖9 長江流域水儲量變化周年振幅空間分布(a),(b),(c)分別表示基于GNSS、GRACE、GLDAS的周年振幅，(d)表示多年平均降水量.Fig.9 Spatial distribution of annual amplitude of water storage change in the Yangtze River Basin(a), (b), and (c) represent the annual amplitudes based on GNSS, GRACE, and GLDAS, respectively. (d) represents the multiyear mean precipitation.

圖10展示了各數(shù)據(jù)集在長江流域及其子區(qū)域按緯度加權(quán)平均后的水儲量變化時(shí)間序列，其中GNSS-EWH被重新采樣成每月的結(jié)果，以便與其他數(shù)據(jù)集進(jìn)行比較.基于GNSS的子區(qū)域平均水儲量變化量級仍與水儲量變化空間分布相匹配，金沙江流域平均水儲量變化量級最大，傳統(tǒng)上游與中下游相當(dāng).其中金沙江流域平均水儲量變化時(shí)序周年振幅可達(dá)137 mm，分別是中下游的2.3倍、傳統(tǒng)上游的1.9倍. GRACE-EWH在金沙江流域明顯小于GNSS-EWH的結(jié)果，主要可能與GNSS水文負(fù)荷位移時(shí)序中包含的其他未剔除的信號有關(guān)；其次，GRACE低空間分辨率導(dǎo)致的信號泄漏與衰減(Chen et al., 2006)也可能存在一定的影響.GNSS-EWH與GRACE-EWH展示了所有水成分的整體變化情況，而GLDAS模型不包括深部地下水和地表水(河流、湖泊等)成分，在水系較多、徑流量大的長江流域，與前兩者結(jié)果相比水儲量變化必然會被低估.GNSS垂向位移不僅記錄了區(qū)域內(nèi)水文負(fù)荷變化，而且對局部質(zhì)量變化具有很高的敏感性(Knappe et al., 2019)，由GNSS約束的結(jié)果能夠反映更為精細(xì)的水文負(fù)荷效應(yīng)(王偉等, 2017)，考慮采樣頻率、建模方法及觀測手段等因素，幾種數(shù)據(jù)集結(jié)果可能會存在一定的差異.

圖10 長江流域及其子區(qū)域水儲量變化時(shí)間序列(a), (b), (c), (d)分別表示整個(gè)長江流域以及金沙江流域、傳統(tǒng)上游、長江中下游的區(qū)域水儲量變化時(shí)間序列.青色線之間表示GRACE與GRACE-FO的銜接時(shí)段.Fig.10 Time series of water storage change in the Yangtze River Basin and its sub-regions(a), (b), (c), and (d) represent the time series of regional water storage change in the whole Yangtze River Basin,the Jinsha River Basin, the traditional upper reaches, and the middle and lower reaches of the Yangtze River. Cyan lines represent the transitional period of GRACE and GRACE-FO.

3.4 時(shí)變水儲量變化與極端事件

為了深入研究季節(jié)性水儲量變化的時(shí)間特征，本節(jié)對長江流域及其子區(qū)域2011年1月—2020年12月間的水儲量變化時(shí)序進(jìn)行了分析，量化了地表降水與水儲量變化的時(shí)間相關(guān)性和時(shí)滯特征，并對GNSS探測流域尺度極端降水事件的能力進(jìn)行了驗(yàn)證.

從圖10a可以看出，GNSS-EWH能夠表征長江流域水儲量變化的季節(jié)性，與月降水量的相關(guān)性較好(CC=0.59)，但兩者存在明顯的相位滯后關(guān)系.各數(shù)據(jù)集的互相關(guān)系數(shù)及滯后月數(shù)在表1中列出，表中大部分互相關(guān)系數(shù)通過了Monte Carlo模擬一階自回歸時(shí)間序列(AR1)置信度為95%的顯著性檢驗(yàn)(高春春等, 2019).其中，GNSS-GRACE、GNSS-GLDAS、GNSS-PRECI、GRACE-GLDAS、GRACE-PRECI以及GLDAS-PRECI的置信水平值分別為0.45、0.45、0.46、0.46、0.48和0.47.GNSS-EWH出現(xiàn)峰值的時(shí)間較降水峰值出現(xiàn)時(shí)間平均晚2～3個(gè)月，在亞馬遜盆地等熱帶地區(qū)也出現(xiàn)了相同的現(xiàn)象(Humphrey et al., 2016)，Hsu等(2020)認(rèn)為這種現(xiàn)象與水資源儲存和運(yùn)輸?shù)膹?fù)雜動(dòng)態(tài)過程有關(guān).需要指出的是，水儲量與降水量作為不同的物理量(二者為微分關(guān)系)，探索兩者之間的定量關(guān)系還需進(jìn)行后續(xù)的數(shù)據(jù)處理.對比子流域的時(shí)間序列結(jié)果(圖10b—d)，GNSS-EWH與其他數(shù)據(jù)集相關(guān)系數(shù)均大于0.48，在金沙江流域GNSS-EWH與GLDAS-EWH的相關(guān)系數(shù)可達(dá)0.84，說明GNSS-EWH反映的流域水儲量變化情況與其他數(shù)據(jù)集反映的水文情況具有較強(qiáng)的一致性.在金沙江流域與傳統(tǒng)上游，GNSS-EWH出現(xiàn)峰值時(shí)間比降水出現(xiàn)峰值時(shí)間晚2個(gè)月，而長江中下游的時(shí)間滯后不明顯(小于1個(gè)月)，達(dá)到峰值時(shí)間普遍早于上游與金沙江，推測可能是由于梅雨天氣與局部空間變異所導(dǎo)致，整個(gè)長江流域的時(shí)滯性特征主要由金沙江流域與傳統(tǒng)上游貢獻(xiàn).

表1 長江流域及其子區(qū)域不同數(shù)據(jù)集之間的相關(guān)系數(shù)與滯后月數(shù)Table 1 Correlation coefficients and lag months between different datasets in the Yangtze River Basin and its sub-regions

長江中下游作為我國水資源配置的戰(zhàn)略水源地，區(qū)域內(nèi)人口密集、工業(yè)發(fā)達(dá)，但極端水文事件常有發(fā)生，因此本文以長江中下游極端事件為例進(jìn)行研究.圖11b展示了基于GNSS-EWH計(jì)算的長江中下游10年間的水儲量盈虧情況，其中，水儲量盈虧通過圖11a展示的水儲量變化估計(jì)值與氣候?qū)W值(參考值)的差值產(chǎn)生.統(tǒng)計(jì)發(fā)現(xiàn)，長江中下游在2012、2014—2016年均發(fā)生了水儲量盈余的現(xiàn)象，而在2011、2013、2018年水儲量發(fā)生虧損，這些結(jié)果與Zhang等(2015)、Ma等(2018)、陳威等(2020)、Jiang等(2021a)的描述相符.除地形、氣候與人為(Chao et al., 2021)等因素外，相關(guān)資料顯示2011年春季長江中下游的水儲量虧損可能是由于La Nia現(xiàn)象導(dǎo)致的干旱造成(Zhang et al., 2015)，2015—2016年水儲量盈余的現(xiàn)象也與該時(shí)期發(fā)生超強(qiáng)El Nio有關(guān)(陳威等, 2020; Ma et al., 2018).圖11c中展示的降水異常與GNSS估算的水儲量盈虧情況較為吻合，兩者互相關(guān)系數(shù)為0.49，說明長江中下游水儲量變化在一定程度上受降水驅(qū)使，該區(qū)域極端事件的發(fā)生往往與降水量聯(lián)系緊密，且降水異常的反饋在陸地水儲量異常結(jié)果中會延遲顯示.通常情況下，實(shí)測降水不足以描述極端水文事件(Famiglietti and Rodell, 2013)，而扣除了大氣和非潮汐負(fù)荷等影響后，GNSS負(fù)荷垂直位移估算結(jié)果是所有水成分變化的綜合體現(xiàn)，可獨(dú)立表征水文極端情況，無需附加其他觀測約束.總的來說，連續(xù)運(yùn)行GNSS網(wǎng)絡(luò)即使在稀疏的空間覆蓋下，也可以作為識別和量化水文極端情況的補(bǔ)充工具.

圖11 長江中下游水儲量盈虧情況及區(qū)域內(nèi)降水異常Fig.11 Water storage surpluses and deficits in the Yangtze River Basin′s middle and lower reaches, as well as precipitation anomalies in this region

4 結(jié)論

本文基于稀疏的GNSS臺站網(wǎng)絡(luò)數(shù)據(jù)，結(jié)合Slepian基函數(shù)獲取長江流域GNSS網(wǎng)位移的局部譜，并利用質(zhì)量負(fù)荷理論實(shí)現(xiàn)從負(fù)荷位移的空間變化到水儲量的空間變化的轉(zhuǎn)化，得到了長江流域及其子區(qū)域的水儲量變化.通過模擬實(shí)驗(yàn)，對反演模型的不確定性進(jìn)行了評估，驗(yàn)證了反演模型的有效性.將GNSS估算結(jié)果與GRACE、GLDAS及氣象降水三個(gè)數(shù)據(jù)集進(jìn)行對比，討論分析了長江流域及其子區(qū)域水儲量變化的時(shí)空分布模式，展示了GNSS獨(dú)立監(jiān)測極端水文事件的能力.研究結(jié)果表明：

(1)站點(diǎn)水文負(fù)荷位移與其推算的水文負(fù)荷分布較為吻合，在西南地區(qū)與江西湖口站(JXHK)水文負(fù)荷位移周年振幅較大，相對應(yīng)的水儲量變化在該站點(diǎn)附近也較為顯著.此外，利用該反演模型對長江流域水儲量進(jìn)行約束，當(dāng)邊界擴(kuò)充為3°時(shí)能夠獲得穩(wěn)定的結(jié)果，在此條件下，反演模型造成的截?cái)嗾`差對原始位移信號影響較小.本文將Mascon模型解計(jì)算的位移作為模擬值，估算出的截?cái)嗾`差量級～1 mm，其中邊界范圍內(nèi)站點(diǎn)截?cái)嗾`差均小于0.8 mm，截?cái)嗾`差對站點(diǎn)位移的平均占比僅為8.1%.

(2)多種數(shù)據(jù)集的結(jié)果在空間上表現(xiàn)出一致的分布情況，GNSS結(jié)果與GLDAS、GRACE空間相關(guān)度分別為91%和79%.具體表現(xiàn)在：西南地區(qū)降水豐沛、持水能力強(qiáng)，在金沙江流域水儲量變化周年振幅最大，并呈現(xiàn)西南向東北遞減的趨勢；長江中下游水系復(fù)雜，且受梅雨天氣與亞洲夏季風(fēng)影響，中下游東部平原局部變化也較為明顯.由于GNSS位移在包含了更為細(xì)節(jié)的局部質(zhì)量空間分布的同時(shí)也可能包含了其他非線性地球物理信號，而GRACE存在信號泄漏、GLDAS存在建模不完整的問題，幾種數(shù)據(jù)集的結(jié)果也存在一定的差異.GNSS估算的長江流域水儲量變化顯著，其最大周年振幅可達(dá)～214 mm，分別是GRACE-EWH最大周年振幅的1.8倍、GLDAS-EWH最大周年振幅的2倍.

(3)不同數(shù)據(jù)集在長江流域及其子區(qū)域水儲量變化均表現(xiàn)出明顯的季節(jié)變化與年際變化.GNSS與GLDAS、GRACE的時(shí)序相關(guān)系數(shù)均在0.56～0.84之間，在金沙江流域和傳統(tǒng)上游，GNSS結(jié)果較降水出現(xiàn)峰值的時(shí)間晚2月左右，而在長江中下游地區(qū)幾乎沒有時(shí)延現(xiàn)象(時(shí)間滯后小于1個(gè)月).利用GNSS能夠監(jiān)測到降水異常的情況，且對GNSS反演結(jié)果水儲量變化異常情況統(tǒng)計(jì)發(fā)現(xiàn)，長江中下游在2012、2014—2016年均發(fā)生了水儲量盈余的現(xiàn)象，在2011、2013、2018年水儲量發(fā)生虧損，這些結(jié)果在相關(guān)文獻(xiàn)中有相應(yīng)的報(bào)道，證明GNSS具有獨(dú)立監(jiān)測極端水文事件的能力.

GNSS反演結(jié)果與多源數(shù)據(jù)集一致的時(shí)空特征驗(yàn)證了Slepian函數(shù)法反演陸地水儲量變化的可靠性，這也為多源數(shù)據(jù)融合提供了一種新的思路.同時(shí)，基于GNSS的反演結(jié)果能夠反映所有水成分整體的變化或異常情況，這也為區(qū)域干旱、洪水的量化研究提供了新的數(shù)據(jù)集.而Slepian函數(shù)法作為一種局部譜集中的方法，密集的地面觀測能夠更加充分地約束地表負(fù)荷的變化，增加站點(diǎn)覆蓋仍然是提升反演結(jié)果空間分辨率、獲取更加精細(xì)水儲量變化空間分布特征的重要手段.

致謝感謝三位審稿專家的寶貴建議.感謝國家地震科學(xué)數(shù)據(jù)中心-東部形變數(shù)據(jù)服務(wù)中心提供的GNSS坐標(biāo)精密時(shí)序，感謝中國氣象局提供的降水與地表溫度數(shù)據(jù)以及德克薩斯大學(xué)空間研究中心(Center for Space Research, CSR)、德國地學(xué)研究中心(Geo Forschungs Zentrum, GFZ)提供的相關(guān)數(shù)據(jù)產(chǎn)品.并衷心感謝美國普林斯頓大學(xué)Simons教授提供的Slepian相關(guān)源碼和耐心解答以及南陽師范學(xué)院高春春副教授給予的指導(dǎo).