許錫文 張志偉 李辰風(fēng) 席瑞杰

1 江西應(yīng)用技術(shù)職業(yè)學(xué)院測繪地理信息學(xué)院，江西省贛州市文峰路9號，341000 2 長江三峽勘測研究院有限公司(武漢)，武漢市光谷創(chuàng)業(yè)街99號，430074 3 武漢大學(xué)衛(wèi)星導(dǎo)航定位技術(shù)研究中心，武漢市珞喻路129號，430079

GNSS長期坐標(biāo)時間序列可以反映測站的運(yùn)動趨勢，對于地殼運(yùn)動及動力學(xué)研究至關(guān)重要。國內(nèi)外眾多學(xué)者已對GNSS坐標(biāo)時間序列的噪聲特性進(jìn)行過分析。Mao等[1]、Williams等[2]和Zhang等[3]認(rèn)為，不論是全球還是區(qū)域分布的測站，其坐標(biāo)時間序列的噪聲模型主要表現(xiàn)為閃爍噪聲+白噪聲(FN+WN)的組合。李昭等[4]對11個IGS基準(zhǔn)站近15 a的坐標(biāo)序列進(jìn)行分析，證明中國區(qū)域IGS基準(zhǔn)站的噪聲模型存在多樣性，主要表現(xiàn)為FN+WN和帶通冪律噪聲加白噪聲(BPPL+WN)。張風(fēng)霜[5]的研究表明，F(xiàn)N+WN或冪律噪聲+白噪聲(PL+WN)為云南地區(qū)GNSS各坐標(biāo)分量的最優(yōu)噪聲模型。

雖然以上學(xué)者對GNSS坐標(biāo)時間序列的噪聲模型進(jìn)行了詳細(xì)分析，但是并未研究環(huán)境負(fù)載對噪聲特性的影響，而環(huán)境負(fù)載是GNSS信號呈現(xiàn)周期性變化的主要原因[6]。因此，本文選取云南省25個GNSS觀測站長時間的坐標(biāo)序列，計(jì)算相應(yīng)的環(huán)境負(fù)載位移，采用多種不同的噪聲組合模型來分析環(huán)境負(fù)載改正前后云南省區(qū)域GNSS坐標(biāo)時間序列的噪聲特性變化。

1 數(shù)據(jù)來源

1.1 GNSS數(shù)據(jù)

本文采用中國大陸構(gòu)造環(huán)境監(jiān)測網(wǎng)絡(luò)(簡稱陸態(tài)網(wǎng)絡(luò))云南省25個GNSS測站2010～2017年的坐標(biāo)時間序列，站點(diǎn)的空間分布如圖1所示。采用GAMIT/GLOBK 10.4同時解算所有GNSS觀測數(shù)據(jù)，獲得每日松弛解，然后通過GLOBK采用7參數(shù)轉(zhuǎn)換的方法將松弛解約束到ITRF框架下，具體解算策略可參考陸態(tài)網(wǎng)絡(luò)數(shù)據(jù)處理手冊。

圖1 測站分布Fig.1 Distribution of the stations

1.2 環(huán)境負(fù)載相關(guān)數(shù)據(jù)

環(huán)境負(fù)載主要指水文負(fù)載(HYDL)、非潮汐大氣負(fù)載(NTAL)和非潮汐海洋負(fù)載(NTOL)等，根據(jù)彈性負(fù)載理論[7]可以計(jì)算地球上任意點(diǎn)環(huán)境負(fù)載所引起的地殼形變。本文采用美國宇航局(National Aeronautics and Space Administration, NASA)全球模型與同化辦公室(Global Modeling and Assimilation Office, GMAO)提供的環(huán)境負(fù)載模型產(chǎn)品進(jìn)行雙三次插值得到所有站點(diǎn)的負(fù)載形變，其中水文負(fù)載模型產(chǎn)品采用MERRA-2(modern-era retrospective analysis for research and applications, version 2)，模型分辨率為0.625°×0.5°×1 h；非潮汐大氣負(fù)載模型產(chǎn)品采用GEOS-FPIT，模型分辨率為0.625°×0.5°×3 h×72層；非潮汐海洋負(fù)載模型產(chǎn)品采用由GFZ(German Research Centre for Geosciences)提供數(shù)據(jù)計(jì)算得到的OMCT05(ocean model for circulation and tides)。插值后平均，得到日均負(fù)載位移量，與GNSS時間序列的時間分辨率統(tǒng)一。

2 數(shù)據(jù)分析方法

2.1 噪聲分析方法

采用極大似然估計(jì)分析噪聲時間序列。該方法假設(shè)GNSS坐標(biāo)時間序列由函數(shù)模型和隨機(jī)噪聲組成，其參數(shù)模型為[3]：

x(t)=a+bt+s1sin(2πt)+c1cos(2πt)+

(1)

式中，x(t)為t歷元測站坐標(biāo)序列日解，a為起始時刻測站位置，b為速度項(xiàng)，s1、c1、s2、c2為年周期和半年周期項(xiàng)系數(shù)，hi為由于天線或天線罩更換引起的跳變，nI為跳變的次數(shù)，ti為跳變歷元，H為Heaviside階梯函數(shù)，跳變前其值為0，跳變歷元其值為0.5，跳變后其值為1，為地震跳變的次數(shù)，aj為測站的震后變形幅值，tj、TJ分別為地震歷元和震后弛豫效應(yīng)時間，TJ一般取值為1 a，vt為噪聲項(xiàng)，可表示不同噪聲或者噪聲組合。

一般認(rèn)為，F(xiàn)N+WN是描述GNSS坐標(biāo)時間序列噪聲的最佳模型，但是不少學(xué)者發(fā)現(xiàn)GNSS坐標(biāo)時間序列中也包含有隨機(jī)游走噪聲(RW)、一階高斯-馬爾科夫噪聲(FOGM)等[4,8]。為了更加全面、系統(tǒng)地分析環(huán)境負(fù)載對噪聲特性的影響，本文選取閃爍噪聲加白噪聲(FN+WN)、一階高斯-馬爾科夫噪聲加白噪聲(FOGM+WN)、冪律噪聲加白噪聲(PL+WN)、隨機(jī)游走噪聲、閃爍噪聲加白噪聲(RW+FN+WN)和隨機(jī)游走噪聲加白噪聲(RW+WN)模型，應(yīng)用Hector[9]軟件進(jìn)行分析。

2.2 噪聲模型評價準(zhǔn)則

根據(jù)極大似然估計(jì)準(zhǔn)則，時間序列的對數(shù)值隨著組合噪聲模型v的變化而變化，數(shù)值越大表明模型越優(yōu)。但研究發(fā)現(xiàn)，不能單純地認(rèn)為對數(shù)值越大對應(yīng)的噪聲模型越優(yōu)，因?yàn)槠渲翟酱蟊硎灸Ｐ椭械膮?shù)也越多。為確定最優(yōu)噪聲組合模型，采用貝葉斯信息準(zhǔn)則(BIC)[10]添加懲罰因子，BIC最小表示最優(yōu)。貝葉斯信息準(zhǔn)則定義如下：

BIC=klnN-2lnL

(2)

式中，k為總模型參數(shù)，N為坐標(biāo)時間序列總歷元數(shù)，lnL為模型對應(yīng)的極大似然值。

3 結(jié)果與分析

3.1 環(huán)境負(fù)載改正前的最優(yōu)噪聲模型

采用極大似然估計(jì)及其評價準(zhǔn)則對云南省25個GNSS測站的坐標(biāo)時間序列進(jìn)行上述5種噪聲組合模型的噪聲分析，得到75組環(huán)境負(fù)載改正前的最優(yōu)噪聲組合模型，結(jié)果如表1所示。

表1 環(huán)境負(fù)載改正前后坐標(biāo)時間序列最優(yōu)噪聲模型

3.2 環(huán)境負(fù)載改正后的最優(yōu)噪聲模型

環(huán)境負(fù)載是引起GNSS測站季節(jié)性變化的主要因素?；贜ASA公開的全球環(huán)境負(fù)載格網(wǎng)產(chǎn)品，采用雙三次插值計(jì)算陸地水文負(fù)載、非潮汐大氣負(fù)載及非潮汐海洋負(fù)載對測站E、N、U分量造成的位移，圖2為YNDC站計(jì)算結(jié)果。將相應(yīng)的位移從GNSS坐標(biāo)時間序列中扣除后，再采用極大似然估計(jì)進(jìn)行上述噪聲組合分析，即得到環(huán)境負(fù)載改正后的最優(yōu)噪聲組合模型。

由圖2可知，YNDC站水文負(fù)載產(chǎn)生的位移明顯大于非潮汐大氣負(fù)載和非潮汐海洋負(fù)載，并且在3個方向上都呈現(xiàn)出很強(qiáng)的周期性。

圖2 環(huán)境負(fù)載造成YNDC站的位移Fig.2 The displacements caused by environmental loading at YNDC station

圖3為環(huán)境負(fù)載改正前后最優(yōu)噪聲比例分布圖?？梢钥闯?，E、N、U分量的噪聲模型并不完全相同，5種噪聲模型中，N、U分量分別沒有表現(xiàn)出FOGM+WN、RW+WN和FOGM+WN、RW+FN+WN、RW+WN?？傮w來看，噪聲模型所占比例最大的均為FN+WN，E、N、U分量改正前比例分別為52%、72%、64%；改正后比例分別為56%、64%、72%?？梢姡尤氕h(huán)境負(fù)載改正后并沒有大范圍改變測站的最優(yōu)噪聲模型，僅有部分測站的最優(yōu)噪聲模型有所改變，如YNHZ站E分量的最優(yōu)噪聲模型由PL+WN變?yōu)镕N+WN。究其原因，可能是對于大部分測站而言，環(huán)境負(fù)載改正并沒有掩蓋或者凸顯出某一種噪聲特性，而在某些特定測站的負(fù)載形變卻能夠讓原本被湮沒的信號顯現(xiàn)出來，這與測站周圍的地理環(huán)境有關(guān)。

圖3 環(huán)境負(fù)載改正前后坐標(biāo)時間序列的最優(yōu)噪聲組合模型分布Fig.3 Distribution of the best noise models of coordinate time series before and after environmental loading correction

3.3 環(huán)境負(fù)載對速度不確定度的影響

測站速度及其不確定度對地殼形變的影響至關(guān)重要。以表1計(jì)算的最優(yōu)噪聲模型為準(zhǔn)，以式(1)計(jì)算環(huán)境負(fù)載改正前后各測站的速度及其不確定度。根據(jù)計(jì)算結(jié)果，在環(huán)境負(fù)載改正前后速度項(xiàng)幾乎無差異，故只進(jìn)行速度不確定度差異的討論，速度不確定度計(jì)算結(jié)果如圖4所示。

由圖4可以看出，在正確構(gòu)建噪聲模型的基礎(chǔ)上，環(huán)境負(fù)載改正能夠明顯降低速度不確定度，其中U分量變化最大，E、N分量次之。采用加入環(huán)境負(fù)載改正后建立的最優(yōu)噪聲模型估計(jì)各測站的速度不確定度，對于E、N、U分量，分別有72%、80%、68%的測站的速度不確定度減小。在YNML站，E分量速度不確定度甚至降低2/3，YNMH站N分量及YNYS站U分量均降低超過1/2。但有部分測站的速度不確定度不降反升，如YNZD站U分量變大1/2，這可能與GNSS數(shù)據(jù)預(yù)處理不完全、仍存在系統(tǒng)誤差有關(guān)。由此可以得出，環(huán)境負(fù)載改正能夠影響測站的速度不確定度，在保證噪聲模型正確的前提下提高速度估值的精度。

圖4 環(huán)境負(fù)載改正前后的速度不確定度Fig.4 Velocity uncertainty before and afterenvironmental loading correction

4 結(jié) 語

本文采用極大似然方法及BIC準(zhǔn)則確定了云南省陸態(tài)網(wǎng)絡(luò)25個測站的最優(yōu)噪聲模型，大部分測站3個分量的最優(yōu)噪聲模型表現(xiàn)為FN+WN，其次為PL+WN，少數(shù)測站表現(xiàn)為FOGM+WN、RW+FN+WN及RW+WN。通過計(jì)算分析得出，環(huán)境負(fù)載形變會在一定程度上影響測站坐標(biāo)時間序列的噪聲特性，改變其最優(yōu)噪聲模型。同時，環(huán)境負(fù)載改正能夠減小測站的速度不確定度，提高速度估值的精度，在進(jìn)行高精度地殼形變分析及研究時需顧及此項(xiàng)改正。

需要指出的是，本文采用的GNSS坐標(biāo)時間序列是從2010年開始的，可能無法分析到某些長周期的形變信號，使結(jié)果帶有一定的局限性。