陸地水GNSS反演的格林函數和Slepian基函數比較分析

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

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

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

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

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

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

1 原理與方法

1.1 格林函數

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

(1)

根據Tikhonov正則化原理,建立GNSS反演陸地水儲量估計準則[5,8,20-22]

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

(2)

1.2 Slepian基函數

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

(3)

1.3 數據和處理策略

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

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

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

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

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

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

(4)

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

(5)

本文數據處理策略如圖2。

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

(1) 圖2(a) 數據集:①利用CSR提供的GRACE Mascons格網數據,在研究區域內,按照隨機方式和平均方式各生成100個GNSS測站;②實測GNSS連續站數據;③GRACE球諧系數;④GLDAS-Noah產品。

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

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

2 結果與分析

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

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

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

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

圖5 前8個Slepian基函數的空間分布Fig.5 Spatial distribution of the top 8 Slepian basis functions

基于GRACE Mascon的2015年1月“lwe_thickness”格網數據,從模擬隨機分布的GNSS測站(圖1(a)、圖3(b)中的藍色點)中選取其中的30、60、100個站點,利用格林函數和Slepian基函數分別進行等效水高EWH的反演。格林函數方法計算結果如圖6(a)—(c)所示,Slepian基函數計算結果如圖6(d)—(f)所示。

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

更進一步,定量分析GNSS測站數和Slepian基函數最大截斷數對反演結果的影響。圖7紅線為對應模擬隨機選取不同測站數后反演結果RMS(root mean square)值的變化,多于42個測站后,RMS值變化相對平緩。圖7說明在研究區域,超過42個分布合理的GNSS測站能獲得較為理想的反演結果。圖7藍線為Slepian基函數選取不同截斷數J反演結果RMS值,J=9時,可以獲得最佳效果。

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

表1 幾種方法反演等效水高時間序列相關性系數

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

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

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

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

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

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

(6)

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

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

3 結 語

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