安徽省CORS坐標時間序列噪聲和速度場分析

何 蓉 陶庭葉 丁 鑫 陶征廣

1 合肥工業大學土木與水利工程學院，合肥市屯溪路193號，230009

全球導航衛星系統(global navigation satellite system, GNSS)技術的不斷發展為精確探究地殼運動規律開辟了途徑。通過對長期觀測得到的數據進行精密處理，可以測定板塊運動的速度和方向，進而得到地殼運動形變規律[1-2]。連續運行參考站(continuously operating reference stations, CORS)是采集地理空間信息的重要基礎設施，其能夠提供連續穩定的GNSS觀測數據，為研究區域三維運動狀態提供重要依據。

安徽省衛星定位綜合服務系統(AHCORS)建成于2011年，其包含50余個均勻分布于全省的參考站，是建設“數字安徽”的基礎框架，同時也為研究安徽省地殼運動形變規律積累了數據。安徽省地處我國華東地區，郯廬斷裂帶將其一分為二，以東屬于下揚子斷塊，以西屬于華北斷塊。后者南部的大別山區屬于秦嶺-大別山斷塊，該區域地質構造較為復雜，是地震活躍地區[3-4]，對該區域的地殼運動狀態進行研究具有重要意義。袁鵬等[5]獲取了AHCORS 2011-11～2012-09的數據，首次定量分析安徽省地殼運動狀態，并揭示華北平原的地面沉降狀況。本文在前人研究的基礎上，解算AHCORS 2013-01～2018-06期間的觀測數據，得到安徽省CORS在ITRF2008框架下的坐標時間序列；通過模型擬合得到噪聲序列，并計算噪聲序列的譜指數，從而確定安徽省CORS的最優噪聲模型；同時，顧及有色噪聲的影響，解算AHCORS在ITRF2008框架下的水平速度場和垂直速度場，并給出其相對于歐亞板塊的運動狀態。

1 GNSS數據處理

1.1 數據處理

采用GAMIT/GLOBK軟件對AHCORS 19個測站2013-01～2018-06的觀測數據進行解算。為了提高數據解算的精度，同時解算我國及周邊地區5個IGS站(CHAN、DAEJ、URUM、TCMS、LHAZ)的同期觀測數據。數據處理主要分為2個部分：1)運用GAMIT進行基線解算，得到AHCORS和IGS 測站的松弛單日解；2)通過GLOBK 進行濾波平差，同時引入由SOPAC處理分析的全球子網IGS1-7解算，進而得到AHCORS在ITRF2008框架下的坐標時間序列。

單日解解算時，主要策略包括：1)采用雙差相位解算方法；2)同時解算衛星軌道和測站坐標；3)對流層延遲改正采用維也納映射函數(Vienna mapping function 1, VMF1)，每2 h估計1個天頂濕延遲參數；4)海潮改正采用有限元解2004(finite element solutions 2004, FES2004)網格模型；5)使用IGS精密星歷和地球自轉參數，并給予適當約束，對極移和極移速率分別給予3″/d、0.3″/d的約束，對UT1變化和UT1變化速率分別給予0.2″/d和0.02″/d的約束；6)選擇我國及周邊地區的5個IGS站作為基準站，解算時IGS站點的3個坐標分量均給予0.05 m的約束。

運用GLOBK獲取測站坐標時間序列時，主要解算方案為：1)為了進行全球網平差，引入由SOPAC處理分析的全球子網IGS1-7。在平差解算時將AHCORS單日松弛解與IGS1-7松弛解合并，并將解方案中協方差矩陣的相對權重因子取1.0；2)利用GLOBK綜合各網單日解，并以IGS站在ITRF2008 框架下的坐標和速度為基準，解算AHCORS的坐標時間序列。

圖1給出了滁州來安(CZLA)站和阜陽界首(FYJS)站的坐標時間序列圖，其中2015-02～04由于設備原因出現數據缺失。由圖可見，CZLA站和FYJS站的東分量和北分量均具有較為明顯的線性趨勢，相比之下，CZLA站高程方向的運動趨勢不如FYJS站的明顯，后者垂直分量上出現非常明顯的下降趨勢。

圖1 CZLA站和FYJS站坐標時間序列Fig.1 Coordinate time series of CZLA and FYJS stations

1.2 坐標時間序列分析

GNSS坐標時間序列一般根據如下模型進行擬合，以便更好地研究其線性項、周期項和噪聲項：

y(ti)=a+bti+csin(2πti)+dcos(2πti)+

esin(4πti)+fcos(4πti)+

(1)

式中，ti為時間(單位年)；待定系數a為截距，b為線性速率，c和d為周年項系數，e和f為半周年項系數；gj為遠場地震后產生的同震位移；hj為震后地殼運動速度的改變量；系數kj為震后變形的指數衰減；H(t)為階躍函數；τi為松弛時間；vi為殘差。根據上述參數模型對各個參數進行擬合，可以有效地估計坐標時間序列中的位置和速度信息及其季節性變化。

2 噪聲分析

外界環境變化(如季節性變化)使觀測墩收縮膨脹、解算GNSS觀測數據時引入的校正模型不正確等因素，均會導致GNSS坐標時間序列存在噪聲。研究噪聲之前需要按照以上模型扣除線性和非線性項，得到殘差序列。在討論噪聲序列時，最初認為其中僅包含與時間無關的噪聲，即白噪聲，后經研究認為，其中還含有與時間相關的噪聲。

噪聲的冪律性質可以表示為：

P(f)=P0fα

(2)

式中，P(f)為噪聲序列的功率譜密度；f為該類噪聲對應的頻率；P0和α為待定參數。將式(2)兩邊取對數得：

lnP(f)=lnP0+αlnf

(3)

由式(3)可以看出，譜指數α是功率譜密度對數與噪聲頻率對數分布圖的斜率。譜指數為-3～1之間的實數，不同的譜指數對應的噪聲不同，α=0時表示噪聲類型為白噪聲(white noise, WN)；α=-1時噪聲類型為閃爍噪聲(flick noise, FN)；α=-2時噪聲類型為隨機游走噪聲(random walk noise, RWN)。除了白噪聲，其余噪聲均屬于有色噪聲。

2.1 頻譜分析

采用周期圖法求得參考站各方向噪聲序列的功率譜密度，得到功率譜密度與頻率的雙對數圖，再利用最小二乘擬合截距和斜率得到P0和α。其中CZLA站的功率譜密度與頻率雙對數分布圖如圖2所示。經計算，19個站點的3個坐標分量的噪聲譜指數均在-0.7～0之間，結合文獻[6]可以判斷，噪聲項具有白噪聲和閃爍噪聲結合(WN+FN)的特點。

圖2 CZLA站各分量噪聲序列功率譜密度與頻率雙對數分布Fig.2 The double logarithmic graphics of power spectral density and frequency noise sequence of each component of CZLA station

2.2 顧及有色噪聲對速度估計的影響

GNSS坐標時間序列中含有有色噪聲，僅考慮白噪聲可能會低估速度估計的不確定度[7]。本文利用CATS軟件[8]分別計算在白噪聲(WN)和白噪聲+閃爍噪(WN+FN)2種模型下的速度項并進行對比。表1列出測站坐標時間序列各分量在不同噪聲模型下的速度估計值，其中第1行為WN模型，第2行為WN+FN模型；倍數列是WN+FN模型下各分量速度估計值的精度與WN模型下的比值。由表1可知，在分析GNSS坐標時間序列時，噪聲序列中僅含有白噪聲會導致速度項的精度被明顯高估，且對速度值的估計產生影響。袁鵬等[5]僅采用白噪聲假設得到的速度不確定度可能會過于樂觀。

3 速度場

AHCORS在ITRF2008框架下的水平速度是以地球質心為基準的地殼水平運動速度，這其中包含了歐亞板塊相對于ITRF全球參考框架的運動，且一般情況下板塊運動速度要比變形速度大。因此，要得到安徽省CORS的水平差異性運動，需要將其在ITRF2008 框架下的水平速度場轉化到相對于歐亞板塊的運動速度。根據剛性塊體的旋轉運動模型，在地心坐標系中，如果有一個塊體的絕對歐拉矢量為Ω(ωx,ωy,ωz)，在該塊體上某點(λ,φ)的矢徑為r(x,y,z)，其東向速度和北向速度可表示為[9]：

(4)

3.1 水平速度場

表2(單位mm/a)統計了安徽CORS在ITRF2008框架下的三維速度估值和中誤差(在解算速度時顧及了有色噪聲的影響)。由表可見，AHCORS在ITRF2008框架下的水平運動方向多為東南方向，且N、E方向的精度明顯優于U方向。N方向速度最大值為-9.40 mm/a，最小值為-15.36 mm/a，均值為-12.27 mm/a；E方向速度最大值為31.41 mm/a，最小值為27.59 mm/a，均值為29.25 mm/a。參考站平均水平運動速率為 31.72 mm/a，方向為E22.76°S，與袁鵬等[5]的研究結果一致。圖3 表示AHCORS在ITRF2008框架下的速度場，其中藍色虛線為安徽省區域內斷層[10]。

根據式(4)求得AHCORS以歐亞板塊為運動背景場框架下的水平速度場為：N方向最大值為3.50 mm/a，最小值為-3.00 mm/a，平均值為0.29 mm/a；E方向最大值為8.31 mm/a，最小值為4.60 mm/a，平均值為6.27 mm/a。相對于歐亞板塊運動的水平速度場如圖4所示。由圖可知，以歐亞板塊作為參考框架，安徽省CORS各站的運動方向較為一致，均為東南方向，其平均速率為6.28 mm/a，方向為E2.65°N。另外，計算中國大陸及其周邊5個IGS站點相對于亞歐板塊的速度場，得到的速度趨勢與Zhao等[11]的一致。

3.2 垂直速度場

圖5表示參考站的垂直速度場。由圖可見，安徽省陸地垂直運動的總體規律為：以淮河為界，以南呈隆起趨勢，以北有大面積沉降，垂直方向整體平均運動速率為-0.71 mm/a。其中，大別山地區的隆起速率大于江淮丘陵區和黃山地區，最大隆起速度為6.96 mm/a，該峰值發生在安慶望江(AQWJ)站；江淮丘陵區抬升相對平緩，平均抬升速率為2.0 mm/a；黃山地區的平均隆升速率約為3.6 mm/a。淮北地區沉降明顯，最大沉降速率為32.82 mm/a，發生在宿州碭山(SZDS)站。

表1 測站坐標時間序列各分量的速度估計值及其精度

表2 AHCORS在ITRF2008 框架下三維速度及其中誤差統計

圖3 AHCORS在ITRF2008框架下水平速度場Fig.3 Horizontal velocity field of AHCORS under ITRF2008 frame

圖4 AHCORS在歐亞框架下水平速度場Fig.4 Horizontal velocity field of AHCORS under Eurasian frame

圖5 AHCORS垂直速度場Fig.5 Vertical velocity field of AHCORS

本文結果與袁鵬等[5]對2011～2013年AHCORS坐標時間序列的速度場分析結果相似。主要區別在于，本文淮北地區臺站沉降的最大值為32.82 mm/a，遠大于2011～2013年計算值13.50 mm/a。其原因可能有：1)本文解算的時間跨度比之前文獻長，而時間跨度對噪聲模型和速度估計均有較大影響，在小于10 a的時間跨度內，隨著時間尺度的增加，參考站的速度以及不確定度會逐漸由發散趨于收斂[12-13]；2)隨著礦產資源的開采，該區域的地表沉降速度有逐漸加速的趨勢。

4 結 語

1)為了得到適用于AHCORS坐標時間序列的最優噪聲組合模型，對周期圖法獲取的各方向噪聲序列的功率譜密度進行分析。結果表明，安徽省境內參考站的最佳噪聲模型為白噪聲+閃爍噪聲(WN+FN)模型。若未考慮時間序列中的有色噪聲，則會明顯低估速度項的誤差。

2)采用WN+FN噪聲模型對AHCORS坐標序列建模，獲取ITRF2008框架下三維方向速度場。結果表明，參考站水平方向平均運動速度為31.72 mm/a，方向為E22.76°S。為了進一步獲取AHCORS水平差異性運動，將水平運動速度從ITRF2008框架轉化到相對于歐亞板塊的運動速度。結果表明，以歐亞板塊為參考的平均運動速度為6.28 mm/a，方向為E2.65°N。垂直方向上，以淮河為界呈現出南方隆升、北方沉降的趨勢。淮河以北有大面積沉降，大別山地區的隆起速率大于江淮丘陵區和黃山地區，最大隆起速度為6.96 mm/a；江淮丘陵區抬升速率相對平緩，平均抬升速率為2.0 mm/a；黃山地區的平均隆升速率約為3.6 mm/a。淮北地區沉降明顯，最大沉降速率為32.82 mm/a，發生在SZDS站，該地區的明顯沉降可能是近年來開采礦產資源所致。