徐克科,伍吉倉(cāng)

(1.河南理工大學(xué) 測(cè)繪與國(guó)土信息工程學(xué)院，河南 焦作 454000;2.同濟(jì)大學(xué) 測(cè)繪與地理信息學(xué)院，上海 200092)

GNSS測(cè)站速度估計(jì)模型重建

徐克科1,伍吉倉(cāng)2

(1.河南理工大學(xué) 測(cè)繪與國(guó)土信息工程學(xué)院，河南 焦作 454000;2.同濟(jì)大學(xué) 測(cè)繪與地理信息學(xué)院，上海 200092)

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

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

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

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

(1)

(2)

(3)

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

2 基線向量卡爾曼濾波估計(jì)模型

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

(4)

Dt=

(5)

從k到k+1時(shí)刻的預(yù)報(bào)公式為

(6)

(7)

從k到k+1時(shí)刻的修正公式為

(8)

(9)

(10)

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

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

3 坐標(biāo)時(shí)序擬合估計(jì)模型

對(duì)于長(zhǎng)時(shí)間連續(xù)GNSS坐標(biāo)時(shí)間序列，為獲取GNSS測(cè)站精確的速度場(chǎng)信息，可以對(duì)此采用時(shí)序分析的方法進(jìn)行速度擬合。GNSS單站、單分量坐標(biāo)時(shí)間序列速度擬合函數(shù)模型可表示為

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

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

(11)

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

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

(12)

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

(13)

(14)

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

(15)

式中，n表示時(shí)間序列的總長(zhǎng)度，其對(duì)數(shù)形式為

(16)

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

4 實(shí)例分析

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

圖1 基線SCLH-SCXJ、SCDF-SCYX三分量與基線長(zhǎng)度單日解時(shí)間序列

圖2 測(cè)站SCLH、SCXJ、SCDF和SCYX單日解坐標(biāo)時(shí)間序列及中誤差

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

表1 3種速度估計(jì)模型的估計(jì)結(jié)果 mm/a

5 結(jié) 論

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

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

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

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

[1] 趙國(guó)強(qiáng),蘇小寧.基于GPS獲得的中國(guó)大陸現(xiàn)今地殼運(yùn)動(dòng)速度場(chǎng)[J]. 地震,2014,34(1):97-103.

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

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

[4] 魏子卿,劉光明,吳富梅. 2000中國(guó)大地坐標(biāo)系:中國(guó)大陸速度場(chǎng)[J]. 測(cè)繪學(xué)報(bào),2011,40(4):403-410.

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

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

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

[8] 王解先.GPS精密定軌定位[M].上海:同濟(jì)大學(xué)出版社,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坐標(biāo)時(shí)間序列跨度對(duì)噪聲模型建立的影響分析[J]. 中國(guó)科學(xué)(D輯：地球科學(xué)),2014,44(11):2461-2478.

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

[責(zé)任編輯：劉文霞]

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

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

10.19349/j.cnki.issn1006-7949.2017.04.002

2016-02-17

國(guó)家重點(diǎn)973基礎(chǔ)研究發(fā)展計(jì)劃資助項(xiàng)目(2013CB733304);國(guó)家自然科學(xué)基金資助項(xiàng)目(41404023)

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

P228.4

A

1006-7949(2017)04-0006-06