應(yīng)用GPS數(shù)據(jù)和Slepian基函數(shù)反演川云渝地區(qū)陸地水儲量變化

成帥， 袁林果， 姜中山， 劉中冠， 張迪， 徐小鳳

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

0 引言

陸地水包括了地表水(湖泊、河流和水庫蓄水)、土壤水、地下水、冰雪以及樹冠水等(Famiglietti et al., 2011)，它們在人類的日常生活、農(nóng)業(yè)灌溉、工業(yè)生產(chǎn)等領(lǐng)域都扮演著至關(guān)重要的角色.中國是一個極度缺水的國家，水資源短缺已經(jīng)成為我們國家社會經(jīng)濟(jì)發(fā)展面臨的一個重大資源問題(孫冬營等，2020).川云渝地區(qū)作為南水北調(diào)西線工程的重要源頭(張金良等，2020)，監(jiān)測和量化當(dāng)?shù)仃懙厮淖兓瘜Φ胤秸M(jìn)行干旱與洪澇的預(yù)防、氣候變化分析甚至國家水資源管理政策的制定等都著有重大意義.目前多種現(xiàn)代空間大地測量技術(shù)和水文同化建模已廣泛應(yīng)用于監(jiān)測陸地水的變化.

全球尺度下的水文模型主要包含陸面模型和水文水資源模型，當(dāng)前應(yīng)用最為廣泛的是美國國家航空航天局NASA(Nation Aeronautics and Space Administration)發(fā)布的全球陸面數(shù)據(jù)同化系統(tǒng)GLDAS(Global Land Data Assimilation System)，該模型的優(yōu)點是能提供全球高時空分辨率的多種陸地表狀態(tài)和通量場，例如降水、蒸散、地表徑流和地下徑流等(Rodell et al., 2004).Wang和Zeng(2012)利用青藏高原地區(qū)的63個氣象站點觀測資料驗證了GLDAS能提供較為準(zhǔn)確降水?dāng)?shù)據(jù)，同時氣溫的描述也較好.李霞等(2014)利用水文站觀測數(shù)據(jù)驗證了模型在黃河源區(qū)的適用性，發(fā)現(xiàn)GLDAS數(shù)據(jù)能較好地體現(xiàn)黃河源區(qū)的流量特征.但GLDAS數(shù)據(jù)也存在著一些缺陷，例如Noah模型缺乏部分地表水和深層地下水?dāng)?shù)據(jù)(Jin and Zhang, 2016)，不能完整地監(jiān)測區(qū)域陸地水變化.同時部分地區(qū)存在徑流與蒸散結(jié)果不夠精確.例如，王文等(2017)研究表明GLDAS低估了中國外流區(qū)的徑流以及高估了蒸散作用的影響，且該模型簡化了物理過程，只考慮了氣候因素作用，沒有考慮流域內(nèi)人為因素對水儲量造成的影響(Tatsumi and Yamashiki, 2015).因此，GLDAS數(shù)據(jù)很難以精確反映實際的陸地水變化.

重力恢復(fù)與氣候試驗GRACE(Gravity Recovery and Climate Experiment)是美國國家航空航天局和德國航空航天中心于2002年3月17日發(fā)射的重力衛(wèi)星.其主要目標(biāo)是獲取每月的時變地球重力場來探究水圈、冰層和海洋中的質(zhì)量變化(Tapley et al., 2004).當(dāng)前GRACE提供的時變重力場數(shù)據(jù)已廣泛應(yīng)用于研究全球范圍內(nèi)季節(jié)和年際尺度上的陸地水儲量變化、地下水儲量變化、冰層質(zhì)量再分布等(Seo et al., 2006; Tiwari et al., 2009; Velicogna, 2009).除此之外，GRACE對局部地區(qū)較為明顯的質(zhì)量變化也有良好的響應(yīng)，例如有學(xué)者研究利用GRACE探測三峽水壩以及巴貢大壩的蓄水量變化(Wang et al., 2011; Tangdamrongsub et al., 2019).但GRACE一方面由于條帶誤差和高頻誤差的原因解算時需要進(jìn)行濾波處理(Swenson and Wahr, 2006)，這在一定程度上抑制真實信號特征，另一方面，GRACE在全球范圍內(nèi)只能保持一個月的時間分辨率和約300 km的空間分辨率，導(dǎo)致小范圍內(nèi)的高頻質(zhì)量變化的監(jiān)測受到限制.此外，在GRACE與后期的GRACE Follow-on任務(wù)銜接期間有著接近一年的觀測空缺，因此需要一種代替手段來彌補這部分空缺.

地表季節(jié)性的陸地水變化會導(dǎo)致彈性地球由于負(fù)荷作用產(chǎn)生周期性的形變.Van Dam and Wahr(1998)發(fā)現(xiàn)單個GPS站點的位移時間序列包含了地表負(fù)荷對地球造成的顯著位移，Bevis等(2005)認(rèn)為GPS觀測技術(shù)能對周圍100 km內(nèi)的負(fù)荷變化有著較好的靈敏度，后續(xù)也有許多學(xué)者將GPS位移時間序列與GRACE計算的位移數(shù)據(jù)進(jìn)行對比，發(fā)現(xiàn)GPS時間序列中的季節(jié)性成分主要與陸地水變化有關(guān)(Van Dam et al., 2001; Davis et al., 2004; 王林松等，2014; Gu et al., 2017; 丁一航等，2018).于是有學(xué)者探索利用GPS垂直位移時間序列來研究陸地水儲量變化，例如Argus等(2014)與Fu等(2015)基于GPS站垂直位移時間序列，應(yīng)用圓盤負(fù)荷格林函數(shù)的方法分別反演計算了美國加利福利亞地區(qū)以及美國西部地區(qū)的陸地水儲量變化，Borsa等(2014)基于GPS監(jiān)測的地表隆起計算了美國西部2013年干旱時期流失的水質(zhì)量.國內(nèi)也有學(xué)者根據(jù)GPS監(jiān)測的地表形變來研究了云南境內(nèi)的水儲量變化以及干旱特征(何思源等，2018; Jiang et al., 2017).Hsu等(2020)基于GPS垂直形變反演了臺灣地區(qū)2005至2016年間的陸地水儲量變化，并分析了降雨模式、入滲率、土壤飽和度和徑流等多種因素對陸地蓄水的影響.他們的結(jié)論表明除了GRACE與水文模型以外，GPS可以作為一種獨立的數(shù)據(jù)源來監(jiān)測陸地水變化.且相較于GRACE而言，GPS時間序列時間分辨率高(達(dá)到天或天以下)，監(jiān)測時間長(部分站點觀測時長達(dá)二十余年)，使得GPS在小尺度上的反演計算上有著獨特的優(yōu)勢.

然而格林函數(shù)方法對區(qū)域內(nèi)GPS站點的分布和密度有著較高的要求，不適合用于大尺度陸地水儲量的反演工作(Han and Razeghi, 2017).為了能夠使GNSS用于大尺度水儲量的研究工作，將其轉(zhuǎn)為與GRACE類似的空間譜是一種有效的分析方法.本文探究應(yīng)用局部Slepian函數(shù)將負(fù)荷位移信號轉(zhuǎn)換為Slepian空間譜，并以此反演地表質(zhì)量變化.Slepian問題最早的提出，主要用于解決有限信號截斷問題和提高信號功率譜密度的估算等應(yīng)用(Slepian, 1983)，后來進(jìn)一步完善后逐漸開始應(yīng)用于地球科學(xué)等各個領(lǐng)域.Slepian函數(shù)與GRACE所使用的球諧系數(shù)法類似，是用一定階數(shù)的Slepian基函數(shù)表示局部區(qū)域的地球物理信號場.Simons等(1997)基于該方法研究了金星的重力地形構(gòu)造.Han等用該方法研究了2004年蘇門答臘地震的震后局部重力場變化以及局部月球重力場的優(yōu)化(Han et al., 2008; Han, 2008).Harig和Simons(2012)將該方法應(yīng)用于監(jiān)測冰川消融現(xiàn)象，利用重力衛(wèi)星跟蹤數(shù)據(jù)計算了格陵蘭島地區(qū)2002到2011年間的質(zhì)量變化趨勢.在國內(nèi)該方法也已普遍應(yīng)用于地月重力場方面的相關(guān)研究(朱廣彬等，2012；孫雪梅等，2015；陳石等，2017).在反演陸地水方面，Han和Razeghi(2017)使用Slepian函數(shù)反演了澳大利亞大陸地區(qū)質(zhì)量變化，結(jié)果顯示該方法可以適當(dāng)降低區(qū)域內(nèi)站點空間覆蓋率的要求，實現(xiàn)更大尺度的陸地水儲量變化的反演計算.

川云渝三省位于我國西南地區(qū)，地處青藏高原東南部，地勢西高東低，是青藏高原向東擠壓的構(gòu)造單元中活動性最強(qiáng)的一個，其內(nèi)部結(jié)構(gòu)具有明顯的空間變異性(Jin et al., 2019).氣候由西部青藏高原的高原季風(fēng)氣候向東逐漸過渡為亞熱帶季風(fēng)氣候，區(qū)域內(nèi)水資源豐富，有著較強(qiáng)的季節(jié)性現(xiàn)象且空間差異明顯.本文目標(biāo)是開發(fā)基于Slepian基函數(shù)的GNSS陸地水反演方法并驗證該方法在川云渝地區(qū)的適用性.文中將GPS垂直位移信號應(yīng)用局部Slepian基函數(shù)結(jié)合地表質(zhì)量負(fù)荷理論反演得到2011—2015年川云渝三省陸地水儲量變化，并在空間和時間上與GRACE、水文模型和降雨數(shù)據(jù)進(jìn)行對比來分析該方法的可行性.

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

1.1 GPS位移時間序列及預(yù)處理

本文采用國家地震科學(xué)數(shù)據(jù)共享中心提供的單天解坐標(biāo)時間序列(http:∥www.eqdsc.com/data/pgv-sjxl.htm)，選取了川云渝地區(qū)分布較均勻的72個GPS站，站點分布如圖1所示.本文選取的數(shù)據(jù)時間跨度為2011—2015年，且在研究時段內(nèi)具有較好的完整性.數(shù)據(jù)解算使用麻省理工學(xué)院提供的GAMIT軟件，采用IGS提供的軌道和鐘差信息，先驗對流層延遲采用GMF模型改正，固體潮和極潮改正使用IERS03模型，海潮使用FES2004模型進(jìn)行改正.為了能精確地提取陸地水變化所造成的地表形變信號，本實驗在后續(xù)的時間序列處理中，使用國際質(zhì)量負(fù)荷服務(wù)(http:∥massloading.net/)所提供的大氣負(fù)荷(Merra2模型)和非潮汐海洋負(fù)荷(OMTC05模型)的垂向位移改正產(chǎn)品.圖2為兩者產(chǎn)生的垂向位移在站點處的周年振幅大小.由圖可知，川云渝地區(qū)受到大氣負(fù)荷的影響較為明顯，垂向周年振幅可達(dá)6.5 mm，受到非潮汐海洋負(fù)荷的影響較小，最大約為0.4 mm.在站點原始坐標(biāo)時間序列中將這兩種負(fù)荷信號扣除得到主要由陸地水變化所產(chǎn)生的位移信號.最后為了準(zhǔn)確獲取季節(jié)性位移，本文擬合并去除了時間序列中的線性趨勢以及跳變，以減小構(gòu)造運動，地震以及儀器更換等其他因素的影響.擬合公式如下：

(1)

式中，ti為觀測時間，y為觀測序列，a0為線性項，H(ti-Tgj)為地震與儀器更換所造成的階躍項，Ck

圖1 西南地區(qū)GPS站分布圖Fig.1 Distribution of GPS stations in Southwest China

圖2 站點處大氣負(fù)荷(Non-tidal atmospheric loading: NTAL)與非潮汐海洋負(fù)荷(Non-tidal ocean loading: NTOL)垂直位移周年振幅Fig.2 Annual vertical displacement amplitude of atmospheric loading and non-tidal oceanic loading at GPS station

與Sk表示待估的正余弦振幅，e(ti)為殘差.當(dāng)k為1，2時，時間序列分別擬合周年變化與半周年變化.

為了保證反演計算時數(shù)據(jù)的完整性，本文采用了Kriging-Kalman-Filter模型補齊缺失的GPS數(shù)據(jù)(Liu et al., 2018)，該方法插值時考慮了GPS站之間的空間相關(guān)性，可以更好地恢復(fù)GPS時間序列中的地球物理信號.

1.2 GRACE和GLDAS數(shù)據(jù)處理

本文使用的GRACE時變重力場數(shù)據(jù)來自美國德克薩斯大學(xué)空間研究中心(Center of Space Research，CSR)提供的RL06 v02 Mascon數(shù)據(jù)，本文選取數(shù)據(jù)的時間跨度為2011年至2015年，共47個月.這些Mascon解純粹由GRACE與GRACE-FO信息進(jìn)行正則化約束計算得到，為了便于分析重采樣為0.25°×0.25°的格網(wǎng)形式.數(shù)據(jù)解算時采用衛(wèi)星激光測距的估計值替代C20項和C30項(僅GRACE-FO解)，并根據(jù)Swenson等(2008)計算的估計值校正了地心改正項以及基于ICE6G-D模型修正冰川均衡調(diào)整.Mascon數(shù)據(jù)較好地解決了條帶誤差，不需要再附加濾波以及縮放因子等處理過程，相較于GRACE球諧解能更好地保留真實的質(zhì)量變化特征.

水文模型使用美國航天航空局提供的全球陸地數(shù)據(jù)同化模型GLDAS-Noah(http:∥daac.gsfc.nasa.gocv)，數(shù)據(jù)為格網(wǎng)形式.文中選擇空間分辨率為1°×1°，時間分辨率為一個月的L4表面陸地模型，GLDAS中主要體現(xiàn)了淺層土壤水以及冰雪等成分的變化，將模型中的積雪量、樹冠總蓄水和四層(0～2 m)土壤含水量累加就可得到陸地總水儲量.再去均值后即可得到研究時段內(nèi)的陸地水儲量變化.

2 原理和方法

2.1 Slepian函數(shù)原理

球諧函數(shù)法是在保證全局正交性的條件下將球面上全球分布的地球物理信號用一系列不同階次的球諧系數(shù)表示.但研究時常有目標(biāo)信號只集中于某一特定區(qū)域內(nèi)的情況，Slepian函數(shù)原理是尋找?guī)н吔缦拗频木植壳蛑C函數(shù)，函數(shù)的信號能最佳地集中在單位球表面上的某個封閉部分，進(jìn)而可以表示特定區(qū)域內(nèi)的局部物理場(Simons and Dahlen, 2006)，Slepian函數(shù)只與研究區(qū)域與選取的最大階數(shù)有關(guān).基于最大階數(shù)L解算第a項Slepian基函數(shù)ga的公式為

(2)

θ為余緯，λ為經(jīng)度，Yl m為l階m次的表面球諧函數(shù)，gα,l m為第α項Slepian函數(shù)的l階m次球諧膨脹系數(shù)，它可以由球面Ω上某個封閉區(qū)域R內(nèi)的能量集中度γ來確定.

(3)

Slepian函數(shù)的目的在于得到γ接近于1的一系列g(shù)α,l m值來求得信號集中于R內(nèi)的基函數(shù)gα.方法是用各階次兩個球諧函數(shù)積的積分構(gòu)造矩陣Dl m,l′m′

(4)

解上式中的特征值與特征向量即可得到能量集中度γ與gα,l m.本文選擇封閉積分區(qū)域為西南地區(qū)川云渝三省邊界，最大階數(shù)由積分區(qū)域內(nèi)GPS站點密度決定，試驗后本文選擇60階，計算得到特征值和特征向量分別為γ與gα,l m.圖3為60階Slepian函數(shù)下求得的3721個能量集中度γ.

2008年奧運會之際，峨眉武術(shù)被國務(wù)院列為第二批國家級非物質(zhì)文化遺產(chǎn)。成為中央電視臺《華夏之魂——中國非物質(zhì)文化遺產(chǎn)大觀》專題片重要選題之一。習(xí)近平總書記在相關(guān)報告中也強(qiáng)調(diào)：道路自信心，理論上的自信心和制度上的自信心，最根本的是文化自信，文化自信，是一種更基本、更廣、更深的自信。

圖3 各階次下Slepian函數(shù)能量集中度Fig.3 Concentration ratio of each Slepian function

從圖3可以得出越靠前的項Slepian函數(shù)能量集中度越接近于1，即信號在指定區(qū)域所占的比重越大.圖4為前六項基函數(shù)信號在積分區(qū)域的空間分布.其中紅色表示信號為正，藍(lán)色表示信號為負(fù)，顏色的深淺表示信號的大小，內(nèi)輪廓為實際研究區(qū)域邊界，外輪廓是為了能讓信號能充分覆蓋區(qū)域周圍，將研究區(qū)域向外擴(kuò)充一定緩沖區(qū)后的積分邊界.結(jié)合圖3和圖4也可以得到γ越接近于1時，信號越集中于積分區(qū)域內(nèi).隨著項數(shù)α的增大，基函數(shù)信號逐漸往積分區(qū)域外擴(kuò)散，當(dāng)能量集中度減少至0時，表明信號主要集中于實驗區(qū)域以外.因此一般只需要選擇γ大于一定閾值的前J項基函數(shù)，使其能充分覆蓋積分區(qū)域即可有效轉(zhuǎn)換內(nèi)部所有地球物理信號.本文中的實驗選擇γ大于0.1以上的Slepian基函數(shù).

2.2 Slepian函數(shù)計算等效水高

根據(jù)上節(jié)中Slepian函數(shù)的特性，我們可以將地球表面局部位移場信號用一組Slepian基函數(shù)來表示.各個站點處的垂向位移信號u可以近似為

圖4 前6項Slepian基函數(shù)信號分布α為基函數(shù)項數(shù)，γ表示能量集中度.內(nèi)輪廓為川云渝三省邊界，外輪廓為擴(kuò)充后積分邊界.Fig.4 Signal distribution of the first 6 Slepian basis functionsα is the number of base function items, γ represents the energy concentration.The inner contour is the boundary of Sichuan, Yunnan and Chongqing, and the outer contour is the extended boundary for integral calculation.

(5)

(6)

ρe為地球密度，ρw為水密度，h′l為負(fù)荷勒夫數(shù)，最后得到區(qū)域內(nèi)(θ,λ)處等效水高形式的質(zhì)量負(fù)荷為

(7)

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

3.1 邊界敏感性測試

由于Slepian基函數(shù)是從中心向外平滑擴(kuò)散的，邊界附近遠(yuǎn)離積分區(qū)域中心的Slepian基函數(shù)靈敏度會有所降低(Harig and Simons, 2012).因此為了保證有足夠的靈敏度來響應(yīng)邊界附近的站點信號，也為了驗證該方法反演水儲量變化的可靠性.本文進(jìn)行了恢復(fù)川云渝三省質(zhì)量變化趨勢的實驗來確定緩沖區(qū)的大小.實驗首先將2012年的GLDAS水儲量變化轉(zhuǎn)化為類似于GRACE的60階球諧系數(shù)，并以此計算了區(qū)域內(nèi)站點處相應(yīng)的垂直位移時間序列，用來進(jìn)行不同緩沖區(qū)下的信號恢復(fù)實驗.實驗選擇最大階數(shù)為60，緩沖區(qū)邊界在0°至3°之間，分別計算了各緩沖區(qū)下水儲量變化恢復(fù)情況如圖5所示.圖中第一幅(Simulation)為GLDAS數(shù)據(jù)2012年水儲量周年振幅，后面為不同緩沖區(qū)下恢復(fù)的區(qū)域內(nèi)水儲量變化圖像.可以發(fā)現(xiàn)當(dāng)緩沖區(qū)小于2°時，由于邊界附近Slepian基函數(shù)信號靈敏度較低，邊界會有明顯的信號缺失現(xiàn)象，同時信號大量集中于區(qū)域中央部分.隨著緩沖區(qū)的增大，邊界信號濃度逐漸提升，當(dāng)緩沖區(qū)為2°至3°時趨勢逐漸穩(wěn)定且較好地恢復(fù)了模擬數(shù)據(jù)的水儲量變化，當(dāng)緩沖區(qū)為3°時可恢復(fù)90%以上的模擬水儲量變化.分析模擬實驗中的誤差來源一是Slepian函數(shù)的截斷誤差的影響，二是由于站點分布不均勻，導(dǎo)致無法完全恢復(fù)部分地區(qū)的質(zhì)量變化信號.針對這個問題，陳石等(2017)的相關(guān)實驗表明在站點覆蓋的空白區(qū)容易產(chǎn)生較大的剩余誤差，而均勻分布的站點對Slepian函數(shù)恢復(fù)信號效果有著較好的提升.

3.2 反演計算結(jié)果

基于積分區(qū)域大小內(nèi)站點分布密度，本文最終選擇最大階數(shù)60，緩沖區(qū)為3°，γ>0.1的前31項Slepian函數(shù)將處理后GPS垂向位移時間序列通過式(5)轉(zhuǎn)換為區(qū)域內(nèi)的Slepian位移場，轉(zhuǎn)換后的GPS空間位移場的周年振幅如圖6所示.圖6a為川云渝區(qū)域內(nèi)空間插值后的垂直位移周年振幅，圖6b為用Slepian函數(shù)轉(zhuǎn)化的空間位移場周年振幅.圖中顯示川云渝地區(qū)陸地水造成的垂向季節(jié)性變化有著明顯的從西南至東北減少的趨勢，具體表現(xiàn)為云南以及四川西南部地區(qū)位移周年變化較為明顯，往四川盆.地方向開始變化逐漸減小，最大位于云南西部山地可達(dá)13 mm，最小位于重慶地區(qū)約為4 mm.圖6a和圖6b比較發(fā)現(xiàn)Slepian函數(shù)能有效地將位移信號展開到了整個區(qū)域，圖6b相比圖6a信號過渡更加平滑.分析認(rèn)為這是由于兩種方法恢復(fù)的信號成分不同，空間插值側(cè)重對空間局部點值變化的擬合，而Slepian函數(shù)是從空間譜頻率上合成地球物理信號，該方法對局部的高頻干擾具有明顯的壓制作用.

圖5 不同緩沖區(qū)下模擬信號恢復(fù)實驗(a) GLDAS模擬的陸地水變化; (b)—(h) 不同擴(kuò)充緩沖區(qū)下恢復(fù)的信號變化.Fig.5 Experiment of signal recovery under different buffersThe first sub-picture is the land water change simulated by GLDAS, and the follow pictures are the signal changes recovered under different expansion buffers.

將時變Slepian位移譜場通過式(6)、(7)轉(zhuǎn)換為時變質(zhì)量負(fù)荷分布，在計算過程中用高斯濾波抑制不穩(wěn)定的高頻噪聲信號，得到區(qū)域等效水高形式的地表質(zhì)量變化.為了更好分析陸地水的季節(jié)性變化，本文利用中國氣象數(shù)據(jù)網(wǎng)(http:∥data.cma.cn/)提供的中國地面降水?dāng)?shù)據(jù)計算了區(qū)域平均降水量.采用該降水量與GPS-Slepian、GRACE和GLDAS計算的川云渝地區(qū)陸地水儲量周年振幅進(jìn)行比較，結(jié)果如圖7所示.從圖中可以發(fā)現(xiàn)在空間分布上GPS-Slepian，GRACE和GLDAS結(jié)果有較好的一致性，三種數(shù)據(jù)在區(qū)域內(nèi)分布都表現(xiàn)為南北分化的趨勢，質(zhì)量變化最大均為云南西部地區(qū)，GPS結(jié)果周年振幅達(dá)300 mm左右，大于GRACE的230 mm和GLDAS的160 mm.往川渝地區(qū)開始逐漸減少，最小在四川北部和重慶地區(qū)GPS結(jié)果為70 mm左右，GRACE和GLDAS分別為40 mm和25 mm，結(jié)合降雨數(shù)據(jù)可以發(fā)現(xiàn)川云渝地區(qū)降雨最大位于云南西南部山區(qū)，最大高達(dá)160 mm，這點與前面數(shù)據(jù)結(jié)果保持一致，表明在云南地區(qū)季節(jié)性降水是陸地水變化的一個重要驅(qū)動因子.而在川西高原處降水豐富卻沒有造成明顯的陸地水周年變化，推測原因是區(qū)域內(nèi)地形及地質(zhì)條件不利于降水的存儲.

圖6 (a)空間插值后的周年振幅與(b)Slepian函數(shù)合成的周年振幅Fig.6 (a) Annual amplitude obtained by spatial interpolation and (b) Slepian function synthesis

圖7 川云渝地區(qū)2011—2015年陸地水儲量和降雨數(shù)據(jù)周年振幅Fig.7 Annual amplitude of terrestrial water storage and precipitation data in Southwest China from 2011 to 2015

為了體現(xiàn)川云渝地區(qū)整體陸地水隨時間的變化，本文分別計算了GPS-Slepian、GRACE和GLDAS結(jié)果在地區(qū)內(nèi)的平均等效水高，并結(jié)合區(qū)域內(nèi)平均降水量進(jìn)行了綜合對比(圖8).圖中可以看出，在隨時間變化特征上GPS與GRACE和GLDAS都有著很強(qiáng)的一致性，表現(xiàn)為在每年夏季的8至10月份出現(xiàn)高峰，在冬季的1月至3月達(dá)到低谷.對比降雨數(shù)據(jù)可知川云渝地區(qū)降水集中于夏季，冬季降雨全年最少，這與前面三種數(shù)據(jù)結(jié)論一致.但在峰值出現(xiàn)時間上，GPS結(jié)果要比GRACE與GLDAS要早一定時間，這點在冬季特別明顯，分析原因為西南地區(qū)冬季降雨較少，陸地水變化主要受區(qū)域內(nèi)徑流和蒸散的影響，而GPS對周圍輕微的水文變化更為靈敏，因此造成GPS結(jié)果與GRACE和GLDAS之間存在延遲.姚朝龍等(2019)在西南地區(qū)干旱的研究中也有類似現(xiàn)象.在時間序列大小上，GPS得到的陸地水變化與降雨數(shù)據(jù)相關(guān)性較好，同時發(fā)現(xiàn)在2015年1—2月GPS-Slepian結(jié)果表現(xiàn)出了小幅度的上升，降雨數(shù)據(jù)在2015年1月也出現(xiàn)了異常增大，經(jīng)查詢資料發(fā)現(xiàn)這是由于該時期南方地區(qū)大范圍強(qiáng)降水過程的影響(李妍和李春，2019).但這點在GRACE與GLDAS上沒有明顯體現(xiàn)，表明GPS-Slepian結(jié)果能反映更為真實的局部陸地水變化，而GRACE與GLDAS都在一定程度上抑制了時空信號特征.

3.3 站點時間序列分析

為了體現(xiàn)各站點處GPS結(jié)果與其他數(shù)據(jù)時間序列上的一致性，本文計算了站點處互相關(guān)系數(shù)來判斷三種數(shù)據(jù)間的相關(guān)性.由于地球物理信號類似于紅噪聲，具有顯著的自相關(guān)現(xiàn)象(Grinsted et al., 2004)，因此文中采用了高春春等(2019)分析南極冰蓋地區(qū)所使用的Monte Carlo模擬法進(jìn)行了互相關(guān)函數(shù)的顯著性檢驗.

圖8 區(qū)域內(nèi)平均陸地水儲量變化及降雨時間序列Fig.8 Variation of average terrestrial water storage and precipitation time series

圖9 站點處等效水高時間序列Fig.9 Time series of equivalent water height at GPS stations

本文將2011—2015年內(nèi)GPS反演的等效水高進(jìn)行月平均，GRACE缺失的數(shù)據(jù)采用三次多項式插值方法補齊后，分別計算各個站點處GPS與GRACE，GPS與GLDAS，GRACE與GLDAS三種等效水高時間序列間的互相關(guān)系數(shù)，同時擬合每個站點三種數(shù)據(jù)等效水高的一階自相關(guān)系數(shù)，采用一階自回歸模型(Autoregressive model：AR)來模擬地球物理信號中的紅噪聲信號，模擬的時間序列與站點處等效水高的長度一致，Monte Carlo模擬次數(shù)為5000.如圖10，本文計算了在不同AR1系數(shù)下90%、95%和99%置信水平的相關(guān)顯著臨界值進(jìn)行Monte Carlo模擬來計算相關(guān)顯著的臨界值，比較互相關(guān)系數(shù)與臨界值來判斷各個站點處三種數(shù)據(jù)的相關(guān)程度.

圖10 Monte Carlo模擬的不同AR1模型在90%、95%和99%置信度檢驗時相關(guān)顯著的臨界值Fig.10 The cross-correlation 90%, 95% and 99% confidence threshold from different AR1 model derived by Monte Carlo experiments

文中選取90%置信水平的臨界值作為標(biāo)準(zhǔn)，為了對比明顯，將未達(dá)到臨界值的站點互相關(guān)系數(shù)置為0.由此得到各站點處三種數(shù)據(jù)間的互相關(guān)系數(shù)如圖11所示.比較三幅子圖發(fā)現(xiàn)，GPS、GRACE和GLDAS之間互相關(guān)系數(shù)在云南地區(qū)和四川西部地區(qū)大于四川東部以及重慶地區(qū)，分析是因為川西和云南地區(qū)站點密集且有較明顯的季節(jié)性信號.比較GPS與其他兩種數(shù)據(jù)間的相關(guān)程度，發(fā)現(xiàn)GPS與GRACE在整體上相關(guān)性更強(qiáng)，站點互相關(guān)系數(shù)平均為0.77.而GPS與GLDAS之間的相關(guān)系數(shù)普遍要小，平均互相關(guān)系數(shù)為0.72.其中GPS和GRACE主要在四川東部以及重慶地區(qū)相關(guān)性較差，區(qū)域內(nèi)的有10個站沒有通過90%置信度檢驗，推測是由于該區(qū)域陸地水季節(jié)性較弱，GPS結(jié)果反映的是站點附近局部范圍的陸地水變化，而GRACE反映的是分辨率300 km左右的整體變化特征，因此在季節(jié)性變化較小的地區(qū)兩種數(shù)據(jù)反映的站點處陸地水變化會有不同.在該區(qū)域GPS與空間分辨率較高的GLDAS相關(guān)性明顯更好也可以驗證這一點.子圖三中GRACE與GLDAS之間在云南以及四川西部地區(qū)相關(guān)性最好，根據(jù)圖9分析，這是因為兩者的陸地水變化時間序列上都有著相似的相位以及振幅，使得兩者的等效水高變化特征更為相近.比較分析三幅子圖，發(fā)現(xiàn)GLDAS與GPS和GRACE之間在四川以北的陜甘青地區(qū)相關(guān)性都較差，分析原因可能是這些地區(qū)水資源匱乏，徑流與降雨較少導(dǎo)致GLDAS反映的水文變化不夠精確，且GLDAS不能探測到完整的陸地水變化，導(dǎo)致在部分地區(qū)與另外兩者有差異.

圖11 GPS，GRACE與GLDAS數(shù)據(jù)之間的互相關(guān)系數(shù)，圖中值為0表示站點處互相關(guān)系數(shù)小于90%置信水平的相關(guān)顯著臨界值Fig.11 The cross-correlation coefficient between GPS, GRACE and GLDAS, and value is 0 means that the cross-correlation coefficient is less than 90% confidence level

4 結(jié)論

本文應(yīng)用Slepian函數(shù)反演了西南地區(qū)川云渝三省內(nèi)的陸地水儲量變化，基于GPS坐標(biāo)時間序列數(shù)據(jù)反演得到區(qū)域內(nèi)2011年至2015年的等效水高時間序列.實驗中選擇川云渝三省邊界解算60階的Slepian函數(shù)，為了保證對邊界站點有足夠的靈敏性，本文實驗?zāi)M了將邊界擴(kuò)充0°至3°的條件下對Slepian函數(shù)恢復(fù)質(zhì)量變化的影響.最后通過對比區(qū)域內(nèi)GPS、GRACE、GLDAS得到的等效水高和降雨數(shù)據(jù)的空間分布和時間序列變化來分析反演結(jié)果.本文的主要研究結(jié)論如下：

(1)局部Slepian函數(shù)可以將空間上離散的GPS位移信號轉(zhuǎn)化為位移空間譜.實驗表明Slepian函數(shù)合成的空間位移場能很好的反映實際的位移分布，與插值得到的結(jié)果相比之間的空間過渡更加平滑.川云渝地區(qū)位移周年振幅有著從西南至東北減少的趨勢，表現(xiàn)為四川南部青藏高原、云南地區(qū)位移周年變化較為明顯，往四川平原方向開始變化逐漸減小，最大位于云南西部山區(qū)可達(dá)13 mm，最小位于重慶地區(qū)約為4 mm.

(2)川云渝地區(qū)陸地水儲量周年變化在空間分布上體現(xiàn)為南部的云南地區(qū)大于北部的川渝地區(qū).其中云南西部山地陸地水周年變化最大，四川北部以及重慶地區(qū)總體變化最小.GPS得到的陸地水周年振幅約為7～30 cm，而GRACE和GLDAS僅為GPS的76%和53%左右，這是因為GPS結(jié)果分辨率高且在監(jiān)測局部區(qū)域環(huán)境質(zhì)量變化上更為靈敏，而Mascon數(shù)據(jù)實際分辨率僅為300 km左右，解算時的正則化約束也與200 km高斯濾波效果相當(dāng)(Save et al., 2016)，GLDAS則缺乏了部分地表水以及深層地下水信息，因此兩者振幅相較于GPS偏小.在平均水儲量變化和站點處變化趨勢中可以發(fā)現(xiàn)，GPS結(jié)果與降雨數(shù)據(jù)在時間上更為接近，GPS能得到近實時的監(jiān)測陸地水變化，而GRACE與GLDAS在冬季存在一定的延遲.在極端氣候的監(jiān)測上，GPS結(jié)果與降雨數(shù)據(jù)均能體現(xiàn)出2015年1月的強(qiáng)降雨事件，而GRACE與GLDAS則沒有表現(xiàn)出明顯特征.在時間分辨率上，GPS能監(jiān)測以天為單位的陸地水儲量變化，而GRACE與GLDAS只能反映月變化，所以GPS在監(jiān)測高頻的水文特征上有著一定優(yōu)勢.

(3)GPS反演的站點處等效水高時間序列與GRACE和GLDAS結(jié)果具有較好的互相關(guān)性，均有80%以上的站點通過90%的置信度檢驗.在空間上相關(guān)程度為南部云南地區(qū)強(qiáng)于北部川渝地區(qū)，四川西部地區(qū)強(qiáng)于四川東部以及重慶地區(qū)，原因是云南和四川西部地區(qū)站點更為密集且季節(jié)性現(xiàn)象較明顯.整體來看GPS結(jié)果在時空上與GRACE更為相似，GPS-GRACE之間的互相關(guān)系數(shù)平均為0.77，強(qiáng)于GPS-GLDAS之間的0.72.GPS與GRACE之間主要在四川東部和重慶地區(qū)相關(guān)性較差，分析是因為地區(qū)內(nèi)陸地水的季節(jié)性變化較小且GPS結(jié)果主要反映站點周圍局部的陸地水變化，而GRACE反映的是300 km左右的整體變化特征.GRACE-GLDAS整體相關(guān)性最優(yōu)，但在四川以北的陜甘青地區(qū)效果較差，推測是由于地區(qū)內(nèi)降雨與徑流較少，且GLDAS缺乏地下水等因素綜合導(dǎo)致GLDAS反映的水文變化不精確.

實驗證明GPS站網(wǎng)也可以作為獨立觀測量來有效的監(jiān)測陸地水儲量變化，是GRACE與水文模型之外的另一種可行方案.在反演計算中，Slepian函數(shù)在GPS站分布相對均勻的條件下可以適當(dāng)降低站點密集程度且能應(yīng)用于面積更大的地區(qū).反演結(jié)果可以得到區(qū)域內(nèi)隨時間變化的等效水高時間序列，結(jié)合學(xué)者利用GRACE探測干旱的方法(Thomas et al., 2014)，本文認(rèn)為此結(jié)果可同樣用于監(jiān)測西南區(qū)域內(nèi)的干旱事件.在目前GRACE和GRACE Follow-on之間任務(wù)銜接導(dǎo)致出現(xiàn)大量數(shù)據(jù)空缺的情況下，該方法可以對應(yīng)用GPS手段填補數(shù)據(jù)空缺以及聯(lián)合多種大地測量手段監(jiān)測全球各地區(qū)陸地水儲量變化提供一定的研究和參考.

致謝感謝中國陸態(tài)網(wǎng)提供的GPS坐標(biāo)時間序列數(shù)據(jù)，感謝美國德克薩斯大學(xué)空間研究中心提供的時變重力場模型，感謝美國航天航空局提供的全球陸地數(shù)據(jù)同化模型數(shù)據(jù)，感謝中國氣象數(shù)據(jù)網(wǎng)提供的中國地面降水?dāng)?shù)據(jù)，同時感謝Simons及其團(tuán)隊提供的Slepian譜分析函數(shù)相關(guān)源碼.以及本文實驗中畫圖所用的GMT軟件及其相關(guān)開發(fā)者們在此一并感謝.