徐克科,伍吉倉

(1.河南理工大學 測繪與國土信息工程學院，河南 焦作 454000;2.同濟大學 測繪與地理信息學院，上海 200092)

GNSS測站速度估計模型重建

徐克科1,伍吉倉2

(1.河南理工大學 測繪與國土信息工程學院，河南 焦作 454000;2.同濟大學 測繪與地理信息學院，上海 200092)

獲取準確的GNSS測站速度對于研究全球板塊運動、地殼形變、地震活動及其地球動力學過程至關重要。為此，以GNSS基線解算后的基線向量作為觀測值，重建了最小二乘綜合速度解算模型和卡爾曼濾波速度估計模型。模型中，考慮了測站坐標、運動速度、年、半年周期項一同作為參數(shù)，在基線解網(wǎng)平差的同時，一并求解獲取速度估值。同時，利用坐標時序分析的方法，顧及白噪聲和冪律噪聲的影響，對單日解坐標時間序列重建了時序速度擬合估計模型，以獲取長期趨勢項作為速度值?；?種模型，對川滇地區(qū)中國地殼運動觀測網(wǎng)絡2010—2014年21個GPS基準站的觀測數(shù)據(jù)進行了速度估計與對比分析。結果表明：3種模型所估計測站速度非常接近，差異基本處于0～1 mm/a范圍之內(nèi)；速度估值中誤差均在亞毫米水平。由此得出，3種速度估計模型具有本質(zhì)的一致性，均可正確估計測站運動速度，能夠滿足高精度地殼形變研究的需要。

GNSS；速度；基線向量；卡爾曼濾波；時序分析

要正確理解全球板塊運動、地殼形變構造運動過程及其動力學機制，首先必須先獲取高精度、可靠的測站運動速度[1-6]。如果只觀測兩期，通過GNSS高精度處理軟件經(jīng)赫爾模特坐標變換可得到這兩期的ITRF框架下的坐標，兩期的坐標相減除以兩期間隔時間便可得到測站實際運動速度。然而，地殼形變的研究通常需要經(jīng)過長時間連續(xù)站連續(xù)觀測或區(qū)域站多期的持續(xù)觀測。隨著GNSS觀測資料的持續(xù)積累，通過GNSS長期的觀測數(shù)據(jù)解算的單日解坐標時間序列可以估計出更高精度的測站運動速度。目前，GNSS測站速度估計的方法主要可歸結于兩大類。一是將測站坐標和速度作為參數(shù)，建立關于基線向量的觀測方程，在對基線解網(wǎng)平差的同時一起求解[7-8]。二是對網(wǎng)平差后單日解測站坐標時間序列進行時序分析，擬合出線性趨勢項作為速度[9-12]。有關研究表明，當GNSS坐標時間序列跨度大于4.5 a時，擬合中是否有季節(jié)項的加入對速度估計誤差影響可以忽略不計；低于2.5 a，對速度估計影響較大，不可接受；當時間跨度是整數(shù)年加上半年時，季節(jié)項對速度場的影響最小[13-14]。GAMIT軟件的輸出結果文件H文件提供了單日解基線向量參數(shù)估計結果及其協(xié)方差陣[9，15]，這為利用合并單天解估計速度提供了有利條件。對此，以基線向量作為觀測值，引入測站坐標、運動速度、年、半年周期項一同作為參數(shù)，重建了基線向量最小二乘綜合解算模型、基線向量卡爾曼濾波解算模型；顧及白噪聲和冪律噪聲的影響構建了坐標時序分析速度擬合估計模型。利用川滇地區(qū)中國地殼運動觀測網(wǎng)絡2010—2014年GPS基準站觀測數(shù)據(jù)，分別基于這3種模型進行了GNSS測站運動速度的估計和對比分析。

1 基線向量最小二乘綜合解算模型

(1)

(2)

(3)

上述的觀測方程中，待求參數(shù)為測站i和j在參考時刻m的坐標、測站i，j的運動速度和年、半年周期項系數(shù)，觀測值為基線三維向量，權陣為GNSS基線解算的協(xié)方差陣的逆陣。這只是基于其中一條基線所列出的誤差方程。同理，可列出所有時段整個GNSS網(wǎng)所有基線的觀測方程并聯(lián)立求解。因基線解算結果單天解是松弛解，聯(lián)合解的站坐標組成的網(wǎng)沒有固定的參考基準，因此需要確定參考框架。采用若干穩(wěn)定的IGS站的站坐標和速度進行擬穩(wěn)平差或重心基準平差。

2 基線向量卡爾曼濾波估計模型

利用GAMIT軟件的輸出結果的H文件提供的單日解基線向量參數(shù)估計結果及其協(xié)方差陣，基于動態(tài)卡爾曼濾波理論，將測站坐標、測站速度、年、半年周期項作為具有方差信息的準觀測值，建立了動態(tài)卡爾曼濾波漂移速度估計的狀態(tài)模型。對一系列的單日解進行卡爾曼濾波合并。并通過對全球IGS 核心站的約束來進行參考框架的定義和七參數(shù)轉換，從而獲取了ITRF參考框架下的測站的坐標和運動速度[9，15]。因測站坐標參數(shù)是隨機變化的，采用馬爾科夫過程來描述這些參數(shù)是隨機漫步過程。將測站運動速度、年、半年周期項系數(shù)作為恒定參數(shù)，測站坐標作為時變參數(shù)，列出關于基線向量的卡爾曼濾波模型的觀測方程和狀態(tài)方程。觀測方程同第1節(jié)基線向量最小二乘綜合解算模型中的觀測方程，見式(3)。建立以測站坐標、速度、年、半年周期項系數(shù)作為狀態(tài)向量的狀態(tài)轉移方程，見式(4)。

(4)

Dt=

(5)

從k到k+1時刻的預報公式為

(6)

(7)

從k到k+1時刻的修正公式為

(8)

(9)

(10)

式中，K為卡爾曼濾波增益矩陣。

上述卡爾曼濾波估計的時間順序是逐漸增加的，直到最后參考歷元時刻的參數(shù)最佳估計值和協(xié)方差矩陣，這是前向卡爾曼濾波。還需要再運行后向卡爾曼濾波估計進行回代。后向卡爾曼濾波估計與前向卡爾曼濾波估計公式相同，只是時間順序相反。

3 坐標時序擬合估計模型

對于長時間連續(xù)GNSS坐標時間序列，為獲取GNSS測站精確的速度場信息，可以對此采用時序分析的方法進行速度擬合。GNSS單站、單分量坐標時間序列速度擬合函數(shù)模型可表示為

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

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

(11)

式中:ti為以年為單位的GPS站點單日歷元，H表示Heaviside階梯函數(shù)，參數(shù)a表示站點的起始位置，b表示線速度，c和d組合表示年周期項系數(shù)，e和f表示半年周期項系數(shù)，第七項表示發(fā)生在歷元Tgj處大小為gj的ng個偏移常量，即站點由于地震或者更換接收機天線等原因產(chǎn)生的后續(xù)歷元位置整體偏移。

在速度估計時，考慮到有色噪聲對速度估計的影響，通過構造最大似然函數(shù)估計譜指數(shù)。對于冪律譜噪聲，協(xié)因數(shù)陣Jpl可表示成轉換矩陣與其轉置的乘積[11,16]，即

(12)

轉換矩陣T(a)可表示為

(13)

(14)

式中:σwh，σpl為白噪聲和有色噪聲分量；Rwh，Jpl為白噪聲和有色噪聲協(xié)方差矩陣；Rwh為單位陣。構造噪聲分量和譜指數(shù)的極大似然函數(shù)為

(15)

式中，n表示時間序列的總長度，其對數(shù)形式為

(16)

采用單純形法求解上式中的最大值，得到相應的似然值及白噪聲與有色噪聲分量σwh，σpl。研究結果表明，原始時間序列的最優(yōu)噪聲模型與PL+VW噪聲模型的極大似然值相差甚微，冪律噪聲的水平可以反映時間序列的整體噪聲水平[17]。因為閃爍噪聲和隨機漫步噪聲均屬于冪律噪聲類型，在此，采用PL+VW噪聲模型進行了速度擬合。

4 實例分析

利用川滇地區(qū)中國地殼運動觀測網(wǎng)絡21個基準站(包括框架站)自2010—2014年的觀測數(shù)據(jù)，經(jīng)GAMIT10.4軟件處理后分別得到了單日基線解和單日坐標解，部分站的結果見圖1和圖2。因不同軟件，不同方法的解算結果之間的比較，更能客觀地評價數(shù)據(jù)處理結果的外部一致性與可靠性。為此，同時利用BERNESE5.2軟件也進行了解算，作了比較。由圖1基線結果看，兩種軟件解算的DX、基線長度DL離散度均在5 mm之內(nèi)，DY，DZ離散度在10 mm之內(nèi)；兩者基線時間序列整體吻合度較高，變化趨勢較一致。由圖2坐標解算結果看，兩種軟件網(wǎng)平差后單日解坐標時間序列離散度接近，基本都在10 mm之內(nèi)，中誤差均優(yōu)于3 mm；兩者坐標時序整體趨勢變化一致。可見，數(shù)據(jù)處理過程中所解算基線向量和坐標結果可靠，為后續(xù)測站運動速度估計的準確性提供了保障。

圖1 基線SCLH-SCXJ、SCDF-SCYX三分量與基線長度單日解時間序列

圖2 測站SCLH、SCXJ、SCDF和SCYX單日解坐標時間序列及中誤差

對于基線解，分別利用第1節(jié)基線向量最小二乘方法和第2節(jié)基線向量卡爾曼濾波方法進行了速度估計。對于單日解坐標時間序列，利用第3節(jié)時序分析的方法擬合了測站運動速度。3種方法解算的21個臺站的速度結果見表1。由表1所示，3種模型所解算測站速度值非常接近，之差大都小于1 mm/a；速度估值中誤差均在亞毫米水平。這表明，3種速度估計模型均可正確解算測站運動速度，能夠滿足高精度地殼形變研究的需要。最小二乘和卡爾曼濾波估計模型結果十分接近，兩種模型具有本質(zhì)的一致性。時序擬合模型結果表明，所有測站殘差的有色噪聲功率譜指數(shù)基本處于(-1，0)區(qū)間內(nèi)，表現(xiàn)為閃爍噪聲。不考慮時間序列中的有色噪聲會導致擬合速度精度的嚴重高估。顧及有色噪聲的影響，GPS測站擬合速度值變化不大，但是擬合速度中誤差將擴大約2～8倍。同時，對于5 a的GNSS坐標時間序列，是否有季節(jié)項的加入對速度估計影響微乎其微，可以忽略。

表1 3種速度估計模型的估計結果 mm/a

5 結 論

以GAMIT解算的基線向量作為觀測值，引入測站坐標、運動速度、年、半年周期項一同作為參數(shù)，重建了基線向量最小二乘綜合解算模型、基線向量卡爾曼濾波解算模型。并顧及白噪聲和冪律噪聲的影響構建了坐標時序分析速度擬合估計模型。利用川滇地區(qū)中國地殼運動觀測網(wǎng)絡2010—2014年GPS基準站觀測數(shù)據(jù)，分別基于3種模型進行了GNSS測站運動速度估計和對比分析。

1)基線向量最小二乘綜合解算模型是利用三維自由網(wǎng)解得到的所有基線向量納入統(tǒng)一綜合模型作為觀測值，基線向量的驗后方差-協(xié)方差陣則被用來確定觀測值得權陣，將參考時刻的測站坐標、速度、周年、半周年系數(shù)作為未知參數(shù)，把所有時段所有基線聯(lián)合起來利用最小二乘一同求解。同時需附加一定的約束條件，如固定其中的某些點的坐標速度，或進行擬穩(wěn)平差和重心基準平差等，求出最終解，這通常是法方程迭加過程。隨著觀測期數(shù)的增加，A陣變得很大，對計算機存儲空間要求較高，數(shù)據(jù)處理變得復雜。而基線向量卡爾曼濾波估計勿需保留用過的觀測值序列，按照一套遞推算法，把參數(shù)估計和預報有機地結合起來，計算速度較快，節(jié)約內(nèi)存。坐標時序擬合模型是對GNSS單站N，E，U三分量的坐標時間序列進行時序分析，考慮長期線性運動+年/半年周期運動+階躍，基于譜指數(shù)計算或最小范數(shù)二次無偏估計(MINQUE)或極大似然估計，可以顧及白噪聲、閃爍噪聲、隨機游動、冪指數(shù)噪聲和帶通濾波噪聲等不同噪聲組合，采用加權最小二乘法求取測站長時間的運動速度。同時還可以對不同噪聲特性及產(chǎn)生原因進行深入分析。

2)3種模型速度估計結果表明，三者速度值十分接近，差異大都位于0～1 mm/a范圍之內(nèi)，速度值中誤差均在亞毫米水平。由此驗證了3種速度估計模型的一致性和可靠性，這說明3種模型所估計的測站運動速度均可以滿足高精度地殼形變研究的需要。

3)時序分析結果表明，所有測站殘差的有色噪聲功率譜指數(shù)基本處于(-1，0)區(qū)間內(nèi)?？梢?，5 a的數(shù)據(jù)主要表現(xiàn)為閃爍噪聲。不考慮時間序列中的有色噪聲會導致擬合速度精度的嚴重高估。顧及有色噪聲的影響，相比只考慮白噪聲，GPS測站擬合速度值變化不大，但是擬合速度中誤差將擴大約2～8倍。由于隨機漫步噪聲的長延時相關性，而川滇地區(qū)連續(xù)站坐標時間序列的時間過短，僅5 a時間，很可能還不足以準確量化隨機漫步噪聲的特征。同時，對于5 a的GNSS坐標時間序列，是否有季節(jié)項的加入對速度估計影響微乎其微，可以忽略不計。

[1] 趙國強,蘇小寧.基于GPS獲得的中國大陸現(xiàn)今地殼運動速度場[J]. 地震,2014,34(1):97-103.

[2] 鄒鎮(zhèn)宇,江在森,武艷強,等. 基于GPS速度場變化結果研究汶川地震前后南北地震帶地殼運動動態(tài)特征[J]. 地球物理學報,2015,58(5):1597-1609.

[3] 程鵬飛,文漢江,孫羅慶,等. 中國大陸GPS速度場的球面小波模型及多尺度特征分析[J]. 測繪學報,2015,40(10):1063-1070.

[4] 魏子卿,劉光明,吳富梅. 2000中國大地坐標系:中國大陸速度場[J]. 測繪學報,2011,40(4):403-410.

[5] 張風霜,占偉.利用GNSS連續(xù)觀測資料獲取高精度動態(tài)速度場的研究[J].地震研究,2015,38(1):75-83.

[6] 王琪,張培震,馬宗晉.中國大陸現(xiàn)今構造變形GPS觀測數(shù)據(jù)與速度場[J].地學前緣,2002,9(2):415-429.

[7] 王解先.由GPS基線向量解算地面形變[J].同濟大學學報(自然科學版),2005,33(7):967-970.

[8] 王解先.GPS精密定軌定位[M].上海:同濟大學出版社,1997.

[9] HERRING T A,KING R W,MCCLUSKY S C. GLOBK Reference Manual:Global Kalman filter VLBI and GPS analysis program[R].Release 10.4.Department of Earth,Atmospheric,and Planetary Sciences,Massachusset Institute of Technology. 2010.

[10] WILLIAMS S D P.CATS:GPS coordinate time series analysis software[J].GPS solutions. 2008,12(2):147-153.

[11] LANGBEIN J.Estimating rate uncertainty with maximum likelihood:differences between power-law and flicker-random-walk models[J].Journal of Geodesy,.2012,86(9):775-783.

[12] WILLIAMS S D P.The effect of coloured noise on the uncertainties of rates estimated from geodetic time series[J].Journal of Geodesy,2003,76(9-10):483-494.

[13] BLEWITT G,LAVALLéE D.Effect of annual signals on geodetic velocity[J].Journal of Geophysical Research:Solid Earth (1978-2012),2002,107(B7):ETG 9-1-ETG 9-11.

[14] BOS M S,BASTOS L,FERNANDES R M S.The influence of seasonal signals on the estimation of the tectonic motion in short continuous GPS time-series[J].Journal of Geodynamics, 2010,49(3):205-209.

[15] HERRING T A,KING R W,MCCLUSKY S C. GAMIT Reference Manual:GPS Analysis at MIT[R].Release 10.4.Department of Earth,Atmospheric,and Planetary Sciences,Massachusset Institute of Technology，2010.

[16] 姜衛(wèi)平,周曉慧. 澳大利亞GPS坐標時間序列跨度對噪聲模型建立的影響分析[J]. 中國科學(D輯：地球科學),2014,44(11):2461-2478.

[17] 吳偉偉. 華北地區(qū)GPS連續(xù)站坐標序列特征研究[D].北京：中國地震局地震預測研究所,2014.

[責任編輯：劉文霞]

Velocity estimation model reconstruction for GNSS observation station

XU Keke1,WU Jicang2

(1.School of Surveying and Land Information Engineering,Henan Polytechnic University,Jiaozuo 454000,China；2.College of Surveying and Geo-Information,Tongji University, Shanghai 200092,China)

To acquire reliable velocity of GNSS observation station is very important for the study of global plate movement,crustal deformation,seismic activity and geodynamic process. Therefore, considering the baseline vector after baselines calculation as observation value, introducing the coordinates，speed，annual and semi-annual seasonal items of GNSS stations as parameters together to be estimated, this paper reconstructs the models of least squares estimation and Kalman filter estimation, by which the velocity is determined along with GNSS network adjustment. Meanwhile, considering the effect of the white noise and power law noise, it is same to the velocity fitting estimation model by the coordinate time series using the method of timing series analysis, by which the long-term trend item is estimated and considered as velocity value. Based on three speed estimation models,the GPS data from the Crustal Movement Observation Network of China from 2010 to 2014 in Sichuan-Yunnan area is processed and analyzed.The result shows that the velocity difference between different models is less than 1mm/a and the root mean square error is at a submillimeter level.The three speed estimation models are veritied to be consistent and reliable,and the estimated velocity can satisfy the need the study of high precision crustal deformation..

GNSS; velocity;baseline vector;Kalman filter;timing series analysis

引用著錄：徐克科,伍吉倉.GNSS測站速度估計模型重建[J].測繪工程，2017，26(4)：6-11.

10.19349/j.cnki.issn1006-7949.2017.04.002

2016-02-17

國家重點973基礎研究發(fā)展計劃資助項目(2013CB733304);國家自然科學基金資助項目(41404023)

徐克科(1979-),男，副教授，博士.

P228.4

A

1006-7949(2017)04-0006-06