何 蓉 陶庭葉 丁 鑫 陶征廣

1 合肥工業(yè)大學(xué)土木與水利工程學(xué)院，合肥市屯溪路193號(hào)，230009

全球?qū)Ш叫l(wèi)星系統(tǒng)(global navigation satellite system, GNSS)技術(shù)的不斷發(fā)展為精確探究地殼運(yùn)動(dòng)規(guī)律開(kāi)辟了途徑。通過(guò)對(duì)長(zhǎng)期觀測(cè)得到的數(shù)據(jù)進(jìn)行精密處理，可以測(cè)定板塊運(yùn)動(dòng)的速度和方向，進(jìn)而得到地殼運(yùn)動(dòng)形變規(guī)律[1-2]。連續(xù)運(yùn)行參考站(continuously operating reference stations, CORS)是采集地理空間信息的重要基礎(chǔ)設(shè)施，其能夠提供連續(xù)穩(wěn)定的GNSS觀測(cè)數(shù)據(jù)，為研究區(qū)域三維運(yùn)動(dòng)狀態(tài)提供重要依據(jù)。

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

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

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

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

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

運(yùn)用GLOBK獲取測(cè)站坐標(biāo)時(shí)間序列時(shí)，主要解算方案為：1)為了進(jìn)行全球網(wǎng)平差，引入由SOPAC處理分析的全球子網(wǎng)IGS1-7。在平差解算時(shí)將AHCORS單日松弛解與IGS1-7松弛解合并，并將解方案中協(xié)方差矩陣的相對(duì)權(quán)重因子取1.0；2)利用GLOBK綜合各網(wǎng)單日解，并以IGS站在ITRF2008 框架下的坐標(biāo)和速度為基準(zhǔn)，解算AHCORS的坐標(biāo)時(shí)間序列。

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

圖1 CZLA站和FYJS站坐標(biāo)時(shí)間序列Fig.1 Coordinate time series of CZLA and FYJS stations

1.2 坐標(biāo)時(shí)間序列分析

GNSS坐標(biāo)時(shí)間序列一般根據(jù)如下模型進(jìn)行擬合，以便更好地研究其線性項(xiàng)、周期項(xiàng)和噪聲項(xiàng)：

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

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

(1)

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

2 噪聲分析

外界環(huán)境變化(如季節(jié)性變化)使觀測(cè)墩收縮膨脹、解算GNSS觀測(cè)數(shù)據(jù)時(shí)引入的校正模型不正確等因素，均會(huì)導(dǎo)致GNSS坐標(biāo)時(shí)間序列存在噪聲。研究噪聲之前需要按照以上模型扣除線性和非線性項(xiàng)，得到殘差序列。在討論噪聲序列時(shí)，最初認(rèn)為其中僅包含與時(shí)間無(wú)關(guān)的噪聲，即白噪聲，后經(jīng)研究認(rèn)為，其中還含有與時(shí)間相關(guān)的噪聲。

噪聲的冪律性質(zhì)可以表示為：

P(f)=P0fα

(2)

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

lnP(f)=lnP0+αlnf

(3)

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

2.1 頻譜分析

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

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

2.2 顧及有色噪聲對(duì)速度估計(jì)的影響

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

3 速度場(chǎng)

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

(4)

3.1 水平速度場(chǎng)

表2(單位mm/a)統(tǒng)計(jì)了安徽CORS在ITRF2008框架下的三維速度估值和中誤差(在解算速度時(shí)顧及了有色噪聲的影響)。由表可見(jiàn)，AHCORS在ITRF2008框架下的水平運(yùn)動(dòng)方向多為東南方向，且N、E方向的精度明顯優(yōu)于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。參考站平均水平運(yùn)動(dòng)速率為 31.72 mm/a，方向?yàn)镋22.76°S，與袁鵬等[5]的研究結(jié)果一致。圖3 表示AHCORS在ITRF2008框架下的速度場(chǎng)，其中藍(lán)色虛線為安徽省區(qū)域內(nèi)斷層[10]。

根據(jù)式(4)求得AHCORS以歐亞板塊為運(yùn)動(dòng)背景場(chǎng)框架下的水平速度場(chǎng)為：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。相對(duì)于歐亞板塊運(yùn)動(dòng)的水平速度場(chǎng)如圖4所示。由圖可知，以歐亞板塊作為參考框架，安徽省CORS各站的運(yùn)動(dòng)方向較為一致，均為東南方向，其平均速率為6.28 mm/a，方向?yàn)镋2.65°N。另外，計(jì)算中國(guó)大陸及其周邊5個(gè)IGS站點(diǎn)相對(duì)于亞歐板塊的速度場(chǎng)，得到的速度趨勢(shì)與Zhao等[11]的一致。

3.2 垂直速度場(chǎng)

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

表1 測(cè)站坐標(biāo)時(shí)間序列各分量的速度估計(jì)值及其精度

表2 AHCORS在ITRF2008 框架下三維速度及其中誤差統(tǒng)計(jì)

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

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

圖5 AHCORS垂直速度場(chǎng)Fig.5 Vertical velocity field of AHCORS

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

4 結(jié) 語(yǔ)

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

2)采用WN+FN噪聲模型對(duì)AHCORS坐標(biāo)序列建模，獲取ITRF2008框架下三維方向速度場(chǎng)。結(jié)果表明，參考站水平方向平均運(yùn)動(dòng)速度為31.72 mm/a，方向?yàn)镋22.76°S。為了進(jìn)一步獲取AHCORS水平差異性運(yùn)動(dòng)，將水平運(yùn)動(dòng)速度從ITRF2008框架轉(zhuǎn)化到相對(duì)于歐亞板塊的運(yùn)動(dòng)速度。結(jié)果表明，以歐亞板塊為參考的平均運(yùn)動(dòng)速度為6.28 mm/a，方向?yàn)镋2.65°N。垂直方向上，以淮河為界呈現(xiàn)出南方隆升、北方沉降的趨勢(shì)?；春右员庇写竺娣e沉降，大別山地區(qū)的隆起速率大于江淮丘陵區(qū)和黃山地區(qū)，最大隆起速度為6.96 mm/a；江淮丘陵區(qū)抬升速率相對(duì)平緩，平均抬升速率為2.0 mm/a；黃山地區(qū)的平均隆升速率約為3.6 mm/a?；幢钡貐^(qū)沉降明顯，最大沉降速率為32.82 mm/a，發(fā)生在SZDS站，該地區(qū)的明顯沉降可能是近年來(lái)開(kāi)采礦產(chǎn)資源所致。