聯合GNSS和GRACE數據分析南極跟蹤站的徑向形變

李偉偉 鞠曉蕾 沈云中 張子占

（1同濟大學測繪與地理信息學院，上海200092；2同濟大學空間信息科學與及可持續發展應用中心，上海200092；3中國科學院測量與地球物理研究所大地測量與地球動力學國家重點實驗室，湖北武漢430077）

0 引言

利用GNSS跟蹤站長期觀測數據確定的地殼運動信息（如速度場）對地球動力學的研究發揮了重要的作用。影響GNSS跟蹤站速度場估計精度的主要因素包括跟蹤站坐標序列的長度，周期信號的合理建模，有色噪聲正確確定。通常要求坐標序列的時間跨度大于2.5年，考慮周年與半年周期項［1］，顧及閃爍噪聲和隨機游走噪聲兩類有色噪聲才能準確估計GNSS跟蹤站的速度及其精度［2-6］。近年來，衛星重力測量（如GRACE）從質量變化的角度分析徑向形變也得到了應用。Davis等［7］首次分析比較了亞馬遜流域GRACE時變重力資料與GNSS坐標序列，兩者具有較好的一致性。Van Dam等［8］比較分析了歐洲區域GNSS與GRACE 34個月估計得到的地表垂直位移，并將兩者呈現較低相關性的原因歸結于GNSS數據處理模型的不準確，且在沿海區域的站點表現更明顯。隨后，Tregoning等［9］在 Van Dam研究基礎上，重新處理GNSS與GRACE數據，兩者的相關性有了明顯提高，接近25%的站點相關系數都大于 0.5。Tesmer等［10］分析全球區域 115個站GNSS站觀測坐標和GRACE估算得到的徑向位移序列，60%站相關系數達到0.5以上，34%站相關系數在0和0.5之間，還有6%站呈現負相關。可見，不同區域GNSS與GRACE相關程度不同。

基準站速度主要包括大尺度板塊運動，局部構造形變和冰后回彈三種地殼形變因素的綜合影響。基準站位于基巖上且在板塊內部相對穩定，就可消除局部構造形變的影響。而板塊運動一般沿水平方向，此時的徑向變化則主要是冰后回彈引起的地殼形變［11］。Bevis等［12］利用西南極 GNSS觀測網坐標序列獲得徑向形變值，并比較了 ICE-5G（VM2），IJ05（6A），RF3S20（β＝0.2）以及 HUY09m四種GIA模型的預測徑向形變值。觀測徑向形變值與預測徑向形變值之差的經驗分布函數表示四種GIA模型均存在低估形變或高估形變的現象。埃爾斯沃斯山脈和東南極沿海區域的GNSS站徑向形變結果與前兩種GIA模型預測值比較同樣有一定的偏差［12］，并按站點分析了其偏差的可能原因。

本文主要利用GNSS數據分析南極區域的噪聲特性和形變速度，并與GRACE時變重力場結果進行比較分析，研究兩者在徑向位移的相關性，且比較了GNSS估計形變與GIA模型（W12a）預測站點形變。

1 數學模型

1.1 GNSS估計形變

GNSS跟蹤站坐標序列通常可以用線性趨勢項、周年和半年項、階躍變化項和噪聲項表示［14］，即

式中ti是以年為單位的GNSS站點單日歷元，y0為常數項，r為線性速度，cf和df表示周期性運動（f＝1表示年周期項，f＝2表示半年周期項），表示發生在歷元Tgj大小為gj的階躍，H表示階梯函數。按照最小二乘參數估計的函數模型，式（1）可以改寫為

其中y為觀測值向量，A為設計矩陣，x＝［y0，r，c1，d1，c2，d2，g1，…，gng］為待估的模型參數，e為觀測噪聲向量。若Σy為觀測值的協方差矩陣，則基于最小二乘準則求得的參數估值及其協方差陣為

考慮到白噪聲、閃爍噪聲和隨機游走噪聲的影響，觀測值的協方差矩陣Σy可表示為

其中U1，U2，U3分別為上述三種噪聲的協因數陣，θ1，θ2，θ3為對應的噪聲分量［3］，采用方差-協方差估計公式可計算各種噪聲分量的估值［15］。

1.2 GRACE估計形變

由GRACE估計的徑向位移可用（6）式表示，

其中，Δr為地球表面徑向位移變化，R為地球的平均半徑（6 378.136 km），θ、λ分別為余緯和經度，Nmax為最大階數為完全規格化的勒讓德函數，ΔClm和ΔSlm為每月Stokes系數與系數平均值之差，h′l及 k′l為 l階負荷勒夫數，Wl為高斯濾波函數，主要用來減小高階系數本身誤差導致的影響。本文采用Han等［16］提供的負荷勒夫數，取Nmax為60，高斯濾波半徑為 500 km［17］。

2 數據處理與結果分析

2.1 數據介紹

2.1.1 GNSS數據

由于南極區域各跟蹤站建站時間不同，因此各站采集的數據時間段也不同。為了與GRACE數據時間統一，本文選取了南極大陸12個時間跨度為2004—2012年的GNSS站點作為分析對象，每站數據長度不少于3年。各站經緯度和時間跨度如表1所示。圖1為選取的GNSS站點分布，除CRAR站和MCM4站處于火山地帶外，其他站均建立在基巖上。

表1 南極GNSS站數據概況Table 1.Data overview of GNSS stations in Antarctica

圖1 南極GNSS站點分布Fig.1.Distribution of GNSS stations in Antarctica

2.1.2 GRACE數據

本文采用的GRACE數據為2012年3月CSR（Centre for Space Research）公布的GSM RL05數據，其時間跨度為2004年1月到2012年6月，共99個月的數據。該數據采用了新的地球物理模型（海洋、大氣、潮汐等），其空間分辨率、精度和周期性變化特性等都優于RL04數據［18］。由于GRACE GSM模型已經扣除海洋大氣等的影響，因此在與GNSS數據進行對比時，需要加回level 2數據中反映海洋和大氣質量變化引起的時變重力場信息的GAC產品系數；此外，為了使GRACE數據能夠與GNSS數據的參考框架統一［19］，需要考慮GRACE數據的一階項［20］。

2.2 結果分析

2.2.1 噪聲分量

表2列出了各站估計的白噪聲、閃爍噪聲和隨機游走噪聲分量，由此可明顯看出南極GNSS跟蹤站的觀測噪聲中白噪聲并不是主要的噪聲，因此對其形變速度進行估計時只顧及白噪聲的影響是不符合實際的。另外，除MCM4站以外，大部分站點均不存在隨機游走噪聲或其分量較小，這可能是由于MCM4站附近的火山運動引起基巖不穩定造成的［21］。

表2 噪聲分量Table 2.Noise component

2.2.2 相關性比較

如2.1節所述，GNSS與GRACE數據的時間分辨率不同，為了比較GNSS和GRACE去除線性趨勢項的形變數據，本文對GNSS數據進行了月平均。圖2為12個站點的徑向月變化比較圖（藍色為GNSS，紅色為GRACE），可見GNSS和GRACE確定的形變總體趨勢一致，但前者的變化振幅比后者大。為了研究GNSS和GRACE的相關性，計算了兩者相關系數（如表3，顯著性水平為0.05）。結果顯示兩者的相關系數最大可達0.663，總體上呈現較強的相關性。

2.2.3 徑向形變速度比較

表4列出12個站點GNSS估計的徑向形變速度和W12a模型預測的形變速度。GNSS結果顯示位于南極半島的OHI3和PALM站的徑向形變速度最大，12個站中4個站下降，但變化趨勢均比較緩慢。由于CRAR站和MCM4站位于火山活動地區，其形變速度不能準確反映形變趨勢。其他10個站形變方向基本一致，且6個站與W12a模型吻合較好，W12a模型在DEVI站和DUM1站較高估計形變，而在南極半島較低估計形變。

圖2 GNSS與GRACE數據求得的徑向坐標月變化Fig.2.Vertical displacements derived from GNSS and GRACE data

表3 GNSS與GRACE數據求得的徑向坐標相關系數Table 3.Correlation coefficient of GNSS and GRACE

表4 站點速度（mm·a－1）Table 4.Velocities of GNSS stations（mm·a－1）

3 結論

本文計算了GNSS與GRACE月平均坐標的相關性，顧及有色噪聲估計南極GNSS站點的徑向形變，并比較分析了GNSS站估計形變和GIA模型預測形變。主要得出以下結論：

（1）南極GNSS坐標序列的觀測噪聲主要是白噪聲和閃爍噪聲，然而部分站（如MCM4站）受當地環境影響隨機游走噪聲更明顯，因此在估計形變速度時，只考慮白噪聲是不準確的。

（2）南極區域GNSS站序列與GRACE序列的相關系數最大為0.663，總體呈現較強的相關性。

（3）雖然研究站點形變趨勢總體與W12a模型相吻合，但是需要更多的站點來分析兩者的吻合度。

1 Blewitt G，Lavallée D.Effect of annual signals on geodetic velocity.Journal of Geophysical Research，2002，107（B7）：2145，doi：10.1029／2001JB000570.

2 Zhang J，Bock Y，Johnson H，et al.Southern California Permanent GPSGeodetic Array：Error analysis of daily position estimates and site velocities.Journal of Geophysical Research，1997，102（B8）：18035—18055.

3 Mao A L，Harrision C G A，Dixon T H.Noise in GPS coordinate time series.Journal of Geophysical Research，1999，104（B2）：2797—2816.

4 Teferle F N.Strategies for Long-term Monitoring of Tide Gauges Using GPS.Nottingham：University of Nottingham，2003.

5 Santamaría-Gómez A，Bouin M N，Collilieux X，et al.Correlated errors in GPS position time series：Implications for velocity estimates.Journal of Geophysical Research，2011，116（B1），doi：10.1029／2010JB007701.

6 LiW W，Shen Y Z.Velocity estimation of GPS base stations considering the coloured noises.Journal of Applied Geodesy，2012，6（3-4）：149—157.

7 Davis JL，Elósegui P，Mitrovica JX，et al.Climate-driven deformation of the solid Earth from GRACE and GPS.Geophysical Research Letters，2004，31（24），doi：10.1029／2004GL021435.

8 Van Dam T，Wahr J，Lavallée D.A comparison of annual vertical crustal displacements from GPS and Gravity Recovery and Climate Experiment（GRACE）over Europe.Journal of Geophysical Research，2007，112（B3），doi：10.1029／2006JB004335.

9 Tregoning P，Watson C，Ramillien G，et al.Detecting hydrologic deformation using GRACE and GPS.Geophysical Research Letters，2009，36（L15401），doi：10.1029／2009GL038718.

10 Tesmer V，Steigenberger P，Van Dam T，etal.Vertical deformations from homogeneously processed GRACE and global GPS long-term series.Journal of Geodesy，2011，85（5）：291—310.

11 朱新慧，孫付平.用甚長基線干涉測量數據檢測冰期后地殼回彈.地球物理學報，2005，48（2）：308—313.

12 Bevis M，Kendrick E，Smalley R，etal.Geodeticmeasurements of vertical crustal velocity in West Antarctica and the implications for icemass balance.Geochemistry Geophysics Geosystems，2009，10（10），doi：10.1029／2009GC002642.

13 Argus D F，BlewittG，PeltierW R，etal.Rise of the Ellsworthmountainsand partsof the EastAntarctic coastobserved with GPS.Geophysical Research Letters，2011，38（16），doi：10.1029／2011GL048025.

14 Rosanne N.Observation of Geodetic and Seismic Deformation with the Global Positioning System.San Diego：University of California，2002.

15 Li B F，Shen Y Z，Lou L Z.Efficient estimation of variance and covariance components：A case study for GPS stochasticmodel evaluation.IEEE Transactions on Geoscience and Remote Sensing，2011，49（1）：203—210.

16 Han D Z，Wahr J.The viscoelastic relaxation ofa realistically stratified Earth，and a further analysisof post-glacial rebound.Geophysical Journal International，1995，120（2）：287—311.

17 Jekeli C.Alternative methods to smooth the Earth’s gravity field.Department of Geodetic Science and Surveying.Ohio：Ohio State University，1981.

18 Bettadpur S.Product Specification Document.Center for Space Research，The University of Texas at Austin，2007-b.

19 Swenson S，Chambers D，Wahr J.Estimating geocenter variations from a combination of GRACE and ocean model output.Journal of Geophysical Research，2008，113（B8），doi：10.1029／2007 JB0 05338.

20 Crétaux JF，Soudarin L，Davidson F JM，et al.Seasonal and interannual geocentermotion from SLR and DORISmeasurements：comparison with surface loading data.Journal of Geophysical Research，2002，107（B12）：2374，doi：10.1029／2002JB001820.

21 Raymond C A，Ivins E R，Heflin M B，et al.Quasi-continuous global positioning system measurements of glacial isostatic deformation in the Northern Transantarctic Mountains.Global and Planetary Change，2004，42（1-4）：295—303.