陳 超,鄒 蓉,曹家銘,李 瑜,梁 宏,方智偉

1. 中國地質(zhì)大學(xué)(武漢)地球物理與空間信息學(xué)院,湖北 武漢 430074; 2. 中國地震局中國地震臺網(wǎng)中心,北京 100045; 3. 中國氣象局氣象探測中心,北京 100081

隨著地球氣候以前所未有的速度變暖,全球范圍內(nèi)降水和蒸散的規(guī)模、頻率、時(shí)間和類型正發(fā)生著變化[1-4]。預(yù)計(jì)這些變化將進(jìn)一步加劇,影響到全球水文循環(huán),因此監(jiān)測水資源變化變得更加重要。近年來,GNSS作為一種新的監(jiān)測陸地水儲量(terrestrial water storage,TWS)的大地測量技術(shù)[1-2,4]應(yīng)運(yùn)而生,由此誕生了GNSS水文學(xué)這個(gè)交叉學(xué)科。GNSS監(jiān)測TWS原理與GRACE(gravity recovery and climate experiment)有所不同,GRACE衛(wèi)星重力觀測能直接感應(yīng)TWS(質(zhì)量)變化,而GNSS地殼垂直位移(vertical crustal displacements,VCD)記錄了由于海平面變化、大氣環(huán)流、冰川和冰蓋演變、雪和降雨以及河流、湖泊、土壤和地下水含水層儲水變化等不同表面荷載作用下地球表面變形而引起的地球表面微小但可測量的荷載引起的位移[5-6]。相比GRACE的天基技術(shù),地基觀測的GNSS可借助持續(xù)發(fā)展、通用性強(qiáng)的連續(xù)運(yùn)行參考站網(wǎng)系統(tǒng)(continuously operating reference system,CORS)(例如美國西部板塊邊界觀測實(shí)驗(yàn)室(Plate Boundary Observatory,PBO)、日本地球觀測網(wǎng)絡(luò)(GPS Earth Observation Net work of Japan,GEONET)、歐洲EPN(European Permanent Network)等),無須高額投入。研究表明,GNSS監(jiān)測TWS與GRACE結(jié)果具有廣泛的一致性[7]。密集的GNSS臺網(wǎng)可以得到更高時(shí)空分辨率的TWS變化圖像。例如PBO可將美國西部地區(qū)TWS空間分辨率提升至50 km[1-2],利用GNSS技術(shù)得出的洛基山脈及太平洋西海岸TWS分布狀況圖像可清楚顯示俄勒岡和華盛頓州的季節(jié)性降水集中于Cascade和Olympic山區(qū),而相鄰的Columbia和Harney盆地降水不多,盆山地形過渡帶明顯對應(yīng)了一條干濕分界線,與氣象觀測十分接近[3]。相反,該地區(qū)GRACE圖像很難分辨以上細(xì)微特征[7]。在四川地區(qū)開展的TWS研究也得出類似結(jié)論,GNSS圖像揭示出TWS分布在成都平原、川西高原及川南地區(qū)間的顯著差異,山區(qū)TWS展示出更大的季節(jié)性變化,尤其受西南印度季風(fēng)控制的季節(jié)性降水導(dǎo)致川南一帶TWS的大幅變動(dòng)[8-11]。以上研究結(jié)果表明,GRACE能夠量化300 km范圍內(nèi)(空間尺度)的TWS每月(時(shí)間尺度)變化,而GNSS能反映具有更高空間分辨率(50 km)和時(shí)間分辨率(1 d,甚至數(shù)小時(shí))的TWS精細(xì)變化,與GRACE得到的總體結(jié)果形成互補(bǔ)。

目前,GNSS反演TWS主要有格林函數(shù)[5,8,11]和Slepian基函數(shù)[10,12-13]兩種方法;基于Farrell[6]的空間域卷積格林函數(shù)方法可用于確定50～100 km分辨率的陸地水存儲變化,明顯優(yōu)于GRACE提供的分辨率;當(dāng)研究區(qū)域相對較小(區(qū)域、省份)且GNSS測站覆蓋足夠密集時(shí),格林函數(shù)反演方法較合適。當(dāng)研究范圍較大(國家、大洲)、GNSS網(wǎng)絡(luò)的空間覆蓋不太密集時(shí),格林函數(shù)方法則受到限制[14];Slepian基函數(shù)方法類似于GRACE的大范圍尺度的反演技術(shù)[15]。GRACE重力數(shù)據(jù)通常使用球面諧波(SH)基函數(shù)(整個(gè)球面內(nèi)正交)進(jìn)行分析,而陸地范圍只占整個(gè)球面29.2%,基于局部SH基函數(shù)或球面Slepian基函數(shù)對GNSS數(shù)據(jù)進(jìn)行參數(shù)化和反演TWS,以獲得區(qū)域(而非全球)表面質(zhì)量變化。

理論上,格林函數(shù)與球諧函數(shù)在數(shù)學(xué)模型上是等價(jià)的[16],而在實(shí)際應(yīng)用上卻存在差異。目前,眾多學(xué)者研究GNSS數(shù)據(jù)反演陸地水儲量均采用單一方法,但是在同一地區(qū)不同方法的反演結(jié)果究竟存在多大差異,尚缺少定量分析。因此,本文在同一研究區(qū)域,使用相同數(shù)據(jù),分別采用上述兩種方法反演GNSS等效水高(equivalent water height,EWH),并進(jìn)行結(jié)果對比研究,定量分析格林函數(shù)和Slepian基函數(shù)反演方法的差異。

本文研究區(qū)域選擇在云南地區(qū)(97°E—107°E、20°N—30°N),首先介紹格林函數(shù)和Slepian基函數(shù)反演陸地水的方法和數(shù)學(xué)模型,然后基于模擬數(shù)據(jù)與真實(shí)GNSS連續(xù)觀測站坐標(biāo)時(shí)間序列詳細(xì)對比分析格林函數(shù)和Slepian基函數(shù)反演效果,同時(shí)與GRACE球諧系數(shù)產(chǎn)品、GLDAS水文模型和降水?dāng)?shù)據(jù)進(jìn)行了相關(guān)性分析。本文的研究結(jié)果可以為近實(shí)時(shí)觀測的GNSS在流域尺度研究水文荷載提供參考。

1 原理與方法

1.1 格林函數(shù)

格林函數(shù)(Green's function)是一個(gè)數(shù)學(xué)物理方程[6,17],它表示一種特定的“場”和產(chǎn)生這種場的“源”之間的關(guān)系。對地球物理大地測量領(lǐng)域里的反演而言,線性反演的首要工作是建立線性的函數(shù)模型(線性觀測方程組)[18]。在GNSS水文研究中,格林函數(shù)為地表水單位負(fù)荷對GNSS觀測點(diǎn)的形變貢獻(xiàn)量。線性反演的數(shù)學(xué)模型表述如下[2,17,19]

(1)

根據(jù)Tikhonov正則化原理,建立GNSS反演陸地水儲量估計(jì)準(zhǔn)則[5,8,20-22]

((Ax-u)/σ)2+β2(T)2→min

(2)

1.2 Slepian基函數(shù)

格林函數(shù)反演結(jié)果受GNSS測站數(shù)量和空間分布影響程度大,研究區(qū)域較大時(shí),劃定的格網(wǎng)面積過大會(huì)導(dǎo)致反演效果不佳,因此在大尺度陸地水儲量的反演研究還是采用GRACE技術(shù)為主。然而如前所述,球諧函數(shù)在全球范圍內(nèi)滿足正交性,在局部研究區(qū)域內(nèi)會(huì)存在“泄漏”,球面Slepian基函數(shù)可以最大限度地減少信號從研究區(qū)域“泄漏”[14]。因此,研究者開始研究使用Slepian基函數(shù)反演GRACE和GNSS數(shù)據(jù)以獲得大尺度等效水高[12-14]。

(3)

1.3 數(shù)據(jù)和處理策略

定量對比研究GNSS兩種反演EWH方法,本文收集了兩類實(shí)測空間大地測量數(shù)據(jù)進(jìn)行反演計(jì)算。

(1) 使用美國CSR(Center for Space Research)提供的GRACE/GRACE-FO RL06 Mascon數(shù)據(jù)(圖1(a)),并使用計(jì)算機(jī)隨機(jī)函數(shù)rand()在整個(gè)研究區(qū)域內(nèi)隨機(jī)生成100個(gè)點(diǎn)位作為隨機(jī)分布GNSS模擬站(圖1(a)藍(lán)色)和將1600個(gè)格網(wǎng)平均分成100個(gè)4×4的區(qū)域,每個(gè)區(qū)域內(nèi)隨機(jī)生成1個(gè)點(diǎn)位(共100點(diǎn))作為平均分布GNSS模擬站(圖1(a)洋紅色);

圖1 CSR提供的GRACE-Mascon等效水高以及模擬GNSS測站(洋紅表示平均分布、藍(lán)色表示隨機(jī)分布)Fig.1 GRACE-Mascon equivalent water height provided by CSR and simulated GNSS stations (magenta indicates average distribution and blue indicates random distribution)

(2) 采用中國大陸構(gòu)造環(huán)境監(jiān)測網(wǎng)絡(luò)簡稱“陸態(tài)網(wǎng)絡(luò)”(Crustal Motion Observation Network of China,CMONOC)和中國氣象局(China Meteorological Administration,CMA)建設(shè)的54個(gè)GNSS連續(xù)觀測站(圖1(b)),由圖1(b)地形中清晰看到,云南地區(qū)西北海拔高(最大海拔7600 m),東南海拔低(最低海拔僅76 m),呈現(xiàn)明顯的西北-東南的海拔降低梯度。云南全省25°坡度以下區(qū)域占比56.46%,文獻(xiàn)[27]研究表明,顯著的地形坡度會(huì)影響GNSS反演陸地水的結(jié)果,忽略地形改正使得GNSS反演陸地水的振幅偏大。

GNSS日解坐標(biāo)時(shí)間序列數(shù)據(jù)來源于中國地震數(shù)據(jù)中心(China Earthquake Data Center,CEDC,http:∥www.eqdsc.com/)和中國氣象局。為了更好地反映真正的地表垂直彈性變形,需要對GNSS位置時(shí)間序列中的非水文負(fù)載信號進(jìn)行識別和剔除,本文采用德國地學(xué)中心(Geo Forschungs Zentrum,GFZ)的地球系統(tǒng)建模小組(Earth System Modeling Group,ESM)發(fā)布的地球物理流體荷載產(chǎn)品(http:∥esmdata.gfz-potsdam.de:8080/repository)去除GNSS時(shí)間序列數(shù)據(jù)包含的固體地球表面幾何中心框架(CF)中非潮汐海洋和大氣荷載位移;冰后回彈效應(yīng)導(dǎo)致中國各地地殼抬升0.3～0.6 mm/a[11],本文使用ICE-6G_D模型進(jìn)行校正[11]。GNSS數(shù)據(jù)經(jīng)過原始數(shù)據(jù)解算和后處理改正后得到主要包含由陸地水文變化造成的地表位移波動(dòng)序列。

對于每個(gè)GNSS時(shí)間序列u(t),采用截距b,速率v和余弦、正弦表示的周年和半周年的季節(jié)性項(xiàng)的線性模型

(4)

式中,時(shí)間t以年為單位;周年項(xiàng)的振幅Amp和相位day(最大值出現(xiàn)在一年中的第幾天)計(jì)算方法為

(5)

本文數(shù)據(jù)處理策略如圖2。

圖2 本文數(shù)據(jù)處理策略和流程Fig.2 Data processing strategy and process

(1) 圖2(a) 數(shù)據(jù)集:①利用CSR提供的GRACE Mascons格網(wǎng)數(shù)據(jù),在研究區(qū)域內(nèi),按照隨機(jī)方式和平均方式各生成100個(gè)GNSS測站;②實(shí)測GNSS連續(xù)站數(shù)據(jù);③GRACE球諧系數(shù);④GLDAS-Noah產(chǎn)品。

(2) 圖2(b)方法模型:①模擬數(shù)據(jù)采用格林函數(shù)方法和Slepian基函數(shù)方法反演等效水高;②實(shí)測GNSS連續(xù)站數(shù)據(jù)分別采用兩種反演方法計(jì)算等效水高;③GRACE球諧函數(shù)法;④GLDAS-Noah水文模型。

(3) 圖2(c)水文信號分析:①GNSS反演的等效水高和GRACE Mascons格網(wǎng)數(shù)據(jù)進(jìn)行定量對比分析;②定量比較兩種反演方法計(jì)算的等效水高時(shí)間序列(振幅、相位);③將GNSS反演結(jié)果和GRACE球諧函數(shù)法及GLDAS-Noah水文模型進(jìn)行對比。

2 結(jié)果與分析

本文采用棋盤測試[8]來評估格林函數(shù)反演結(jié)果,圖3顯示了3種GNSS分布策略的棋盤測試結(jié)果(β=0.01);由圖3可以清晰地看出,測站分布越均勻(圖3(b)和(c)對比)、測站數(shù)量越多(圖3(d)和(b)、(c)對比),反演效果越好。

圖3 棋盤測試評估格林函數(shù)反演結(jié)果Fig.3 Green's function inversion checkerboard test

Slepian基函數(shù)空間分布呈現(xiàn)明顯“空限”特征(圖4),本文展開到60階,一共計(jì)算了(L+1)2=3721個(gè)Slepian基函數(shù),其中聚集度大于0.01閾值的共有21個(gè)基函數(shù)(圖5展示了前8個(gè)Slepian基函數(shù)空間分布),其貢獻(xiàn)度達(dá)到99.74%,余下的3700個(gè)Slepian基函數(shù),占總數(shù)的99.44%,其貢獻(xiàn)度不足0.3%。

圖4 展開L=60階,3721個(gè)Slepian基函數(shù)的聚集度Fig.4 Expanded L=60 order, the aggregation degree of 3721 Slepian basis functions

圖5 前8個(gè)Slepian基函數(shù)的空間分布Fig.5 Spatial distribution of the top 8 Slepian basis functions

基于GRACE Mascon的2015年1月“l(fā)we_thickness”格網(wǎng)數(shù)據(jù),從模擬隨機(jī)分布的GNSS測站(圖1(a)、圖3(b)中的藍(lán)色點(diǎn))中選取其中的30、60、100個(gè)站點(diǎn),利用格林函數(shù)和Slepian基函數(shù)分別進(jìn)行等效水高EWH的反演。格林函數(shù)方法計(jì)算結(jié)果如圖6(a)—(c)所示,Slepian基函數(shù)計(jì)算結(jié)果如圖6(d)—(f)所示。

圖6每個(gè)子圖右下角展示了不同GNSS測站數(shù)反演結(jié)果與圖1(a)所示的GRACE Mascons格網(wǎng)數(shù)據(jù)差值統(tǒng)計(jì)的直方分布,可直觀反映反演結(jié)果的好壞。在研究區(qū)域內(nèi)格林函數(shù)方法受測站數(shù)和站點(diǎn)空間分布影響較大(圖6(a)—(c)),如圖6所示,選取60個(gè)GNSS測站(圖6(b)),格林函數(shù)反演即可獲得較好結(jié)果。從差值直方統(tǒng)計(jì)可以發(fā)現(xiàn),Slepian基函數(shù)反演效果比格林函數(shù)反演結(jié)果略差;同時(shí),Slepian基函數(shù)反演結(jié)果對測點(diǎn)空間分布差異和測站點(diǎn)數(shù)多寡敏感程度較弱。

更進(jìn)一步,定量分析GNSS測站數(shù)和Slepian基函數(shù)最大截?cái)鄶?shù)對反演結(jié)果的影響。圖7紅線為對應(yīng)模擬隨機(jī)選取不同測站數(shù)后反演結(jié)果RMS(root mean square)值的變化,多于42個(gè)測站后,RMS值變化相對平緩。圖7說明在研究區(qū)域,超過42個(gè)分布合理的GNSS測站能獲得較為理想的反演結(jié)果。圖7藍(lán)線為Slepian基函數(shù)選取不同截?cái)鄶?shù)J反演結(jié)果RMS值,J=9時(shí),可以獲得最佳效果。

基于上述試驗(yàn)結(jié)果,采用云南地區(qū)陸態(tài)網(wǎng)絡(luò)和氣象局實(shí)測的54個(gè)GNSS垂直位移數(shù)據(jù),(時(shí)間跨度2010—2022年),進(jìn)行GNSS水文研究。格林函數(shù)反演采用β=0.03,Slepian基函數(shù)與上面采用同樣的展開階數(shù)L=60和截?cái)鄶?shù)J=9。在研究范圍內(nèi)不同區(qū)域(滇中、滇東南、滇東北、滇西、滇西南、滇西北)選取6個(gè)測站(圖1(b))反演的等效水高時(shí)間序列(圖8右側(cè));對比分析結(jié)果表明,格林函數(shù)和Slepian基函數(shù)反演的等效水高時(shí)間序列一致性非常高,統(tǒng)計(jì)所有測站兩種方法結(jié)果的相關(guān)系數(shù)達(dá)到0.98(表1)。GRACE球諧函數(shù)方法和GLDAS水文模型與GNSS反演的等效水高時(shí)間序列相關(guān)性均大于0.65,差異體現(xiàn)在GNSS反演的EWH空間分布更精細(xì)。同時(shí),可以發(fā)現(xiàn)EWH時(shí)間序列波峰均出現(xiàn)在降雨量峰值之后,等效水高時(shí)間序列波動(dòng)和降水?dāng)?shù)據(jù)(圖8 Preci)具有明顯的一致性,波峰出現(xiàn)的時(shí)間比最大降雨滯后1～2個(gè)月,可能的原因是該地區(qū)地殼對荷載變化響應(yīng)存在滯后,而且荷載變化被GNSS監(jiān)測并在位移中體現(xiàn)出來需要時(shí)間。

表1 幾種方法反演等效水高時(shí)間序列相關(guān)性系數(shù)

圖8 GNSS垂直位移和反演EWH時(shí)間序列Fig.8 GNSS vertical displacement and inversion EWH time series

降雨數(shù)據(jù)表明,云南地區(qū)常年平均5月最為干燥,9月最為濕潤。通過GNSS的2015年5月和2015年9月的垂直位移和反演等效水高(圖9)進(jìn)一步分析發(fā)現(xiàn),Slepian基函數(shù)反演結(jié)果(圖9(c)、(f))比格林函數(shù)反演結(jié)果(圖9(b)、(e))振幅略大,但兩者空間分布趨勢整體較為一致,相位幾乎一致(圖8(b)、(d)、(f)、(h)、(j)、(l))。

圖9 GNSS垂直位移((a)、(d))、格林函數(shù)反演結(jié)果((b)、(e))、Slepian基函數(shù)反演結(jié)果((c)、(f))Fig.9 GNSS vertical displacement ((a), (d)), Green's function inversion results ((b), (e)), Slepian basis function inversion results ((c), (f))

統(tǒng)計(jì)研究區(qū)域內(nèi)GNSS測站反演的結(jié)果表明,Slepian基函數(shù)反演的振幅比格林函數(shù)反演的振幅平均大25%(圖10、圖11)。筆者推測可能的原因是Slepian基函數(shù)反演GNSS等效水高,存在系統(tǒng)性偏差;GNSS站點(diǎn)的水文荷載導(dǎo)致的垂直位移,是周圍半徑50 km范圍內(nèi)[28]所有水文荷載貢獻(xiàn),而從Slepian基函數(shù)反演數(shù)學(xué)模型能看出,其將周圍所有的水文荷載貢獻(xiàn)都?xì)w集于GNSS站點(diǎn)處的水文荷載,這顯然會(huì)導(dǎo)致Slepian基函數(shù)反演結(jié)果偏大;不同半徑(范圍)的水文荷載對GNSS站點(diǎn)垂直位移貢獻(xiàn)程度除了與距離有關(guān),還與荷載大小有關(guān)系。

圖10 3種不同的格網(wǎng)圓盤及地表荷載形變響應(yīng)Fig.10 Three different grid disks, modeled vertical displacement under a disk load of radius 15.7 km and an equivalent water depth of 1 m

圖10顯示,本文研究區(qū)域內(nèi)采用0.25°×0.25°格網(wǎng),每個(gè)格網(wǎng)邊長約為28 km,3種不同的格網(wǎng)圓盤(圖10(a)內(nèi)切圓半徑14 km、圖10(b)外接圓半徑19.8 km和圖10(c)等面積圓半徑15.7 km);等面積圓能較好地代表每個(gè)格網(wǎng)的信息,本文采用格網(wǎng)等面積圓半徑作為圓盤半徑,圖10(d)為等效水高1 m、圓半徑15.7 km的圓盤造成的地表形變響應(yīng),圖10(e)展示圖10(d)二維剖面,其中藍(lán)色面積為±25 km2范圍內(nèi)荷載響應(yīng)占整體比值0.705。

將Slepian基函數(shù)反演式(3)修改如下

(6)

式中,ω表示Slepian反演方法“系統(tǒng)性偏差”改正參數(shù),其根據(jù)荷載格林函數(shù)矩陣特性計(jì)算;若GNSS測站反映周圍半徑50 km范圍內(nèi)的陸地水荷載(空間分辨率50 km),根據(jù)格林函數(shù)矩陣可以獲得不同格網(wǎng)圓盤半徑的荷載-位移曲線,定義積分計(jì)算曲線在±25 km2以內(nèi)的面積占整體面積的比值為參數(shù)ω,本文研究中ω=0.705。圖11顯示了Slepian反演方法系統(tǒng)性偏差改正效果;ω參數(shù)改正前(圖11(c)、(e))后(圖11(d)、(e)),格林函數(shù)反演結(jié)果振幅與Slepian反演結(jié)果振幅相關(guān)性系數(shù)從0.48上升到0.88。

綜合而言,在本文研究范圍,實(shí)測的GNSS垂直時(shí)間序列反演等效水高,采用格林函數(shù)法和Slepian基函數(shù)法效果相當(dāng),兩種方法均可用于該區(qū)域水文研究。另一方面,統(tǒng)計(jì)Slepian反演結(jié)果的振幅比格林函數(shù)反演結(jié)果的振幅平均大25%,Slepian反演方法存在系統(tǒng)性偏差,應(yīng)根據(jù)荷載格林函數(shù)響應(yīng)特性進(jìn)行改正?；赟lepian反演GNSS陸地水在其他更大空間尺度區(qū)域上(空間分辨率>50 km)是否同樣存在系統(tǒng)性偏差,將在未來研究中更詳細(xì)驗(yàn)證與量化計(jì)算。

3 結(jié) 語

大地測量水文學(xué)研究通常探索水文荷載作用下地球變形與氣候驅(qū)動(dòng)因素之間的經(jīng)驗(yàn)關(guān)系,以改善洪水和干旱等水文事件的預(yù)報(bào)?，F(xiàn)如今,更高空間和時(shí)間分辨率的GNSS技術(shù),能夠研究近乎實(shí)時(shí)觀測的流域尺度水文荷載。本文定量分析格林函數(shù)和Slepian基函數(shù)反演陸地水差異,在研究區(qū)域得到如下結(jié)論:①基于模擬數(shù)據(jù),格林函數(shù)反演結(jié)果受GNSS測站數(shù)量和空間展布影響程度比Slepian基函數(shù)反演結(jié)果更大,而Slepian基函數(shù)反演結(jié)果受最大截?cái)嚯A數(shù)的影響較大;同等情況下(GNSS測站數(shù)≥42),格林函數(shù)反演結(jié)果總體精度比Slepian基函數(shù)反演結(jié)果更優(yōu)。②基于實(shí)測的“陸態(tài)網(wǎng)絡(luò)”和中國氣象局的GNSS連續(xù)測站垂直時(shí)間序列數(shù)據(jù),兩種方法反演等效水高結(jié)果相關(guān)性達(dá)到0.98,Slepian基函數(shù)反演結(jié)果振幅比格林函數(shù)反演結(jié)果振幅平均大25%。③GNSS數(shù)據(jù)反演的陸地水和GRACE、GLDAS推斷的陸地水相關(guān)性均大于0.65,并且與該地區(qū)月度降水?dāng)?shù)據(jù)一致性很好;GNSS反演等效水高序列波峰出現(xiàn)的時(shí)間比最大降雨滯后1～2個(gè)月。