GNSS實時衛(wèi)星鐘差估計在地震監(jiān)測中的應用

王浩浩,郝 明,莊文泉

(中國地震局第二監(jiān)測中心,西安 710054)

0 引言

近年來,隨著全球導航衛(wèi)星系統(tǒng)(global navigation satellite system, GNSS)高頻定位技術的不斷發(fā)展,其在地殼運動監(jiān)測領域發(fā)揮著越來越重要的作用,為地震監(jiān)測提供了一種有效的新型技術途徑[1-2]。其中,利用差分相對定位技術可以計算得到監(jiān)測站相對于某一固定參考站毫米級的動態(tài)位移,但監(jiān)測區(qū)域內(nèi)有時因地質條件等客觀因素難以布設觀測環(huán)境良好、穩(wěn)定的基準站,而且坐標精確已知的基準站在強震發(fā)生階段可能會發(fā)生移動,導致動態(tài)解算的定位精度顯著降低[3-4]。而采用國際GNSS服務(international GNSS service, IGS)組織發(fā)布的精密衛(wèi)星軌道和鐘差等產(chǎn)品的精密單點定位(precise point positioning, PPP)技術,具有不依賴于某一特定參考基準站、實時性強等優(yōu)勢,僅利用單臺GNSS接收機即可獲得國際地球參考框架下高精度的“絕對位置”,具備準確捕捉地震位移波形的能力,更適合長距離、大范圍的地殼形變監(jiān)測[3-4]。文獻[5]實驗結果表明,衛(wèi)星鐘差的采樣率越高,利用PPP動態(tài)解捕捉遠場站點形變信息的優(yōu)勢越明顯。文獻[6]利用武漢大學自主研發(fā)的PANDA軟件對高頻GNSS觀測信號進行PPP后處理,能夠很好地獲取2016-11-13新西蘭Mw7.8地震產(chǎn)生的動態(tài)位移特征。

地震瞬時地表動態(tài)位移的實時高可靠性監(jiān)測對地震預警系統(tǒng)而言至關重要,能夠為震中以及震級的快速確定等研究工作提供關鍵信息[3]。然而,實時PPP高精度動態(tài)解的實現(xiàn)取決于精密衛(wèi)星軌道和鐘差的實時估計性能。目前超快速實時預報軌道產(chǎn)品已經(jīng)能夠滿足其時效性和可靠性等要求,但衛(wèi)星鐘差的精確預報極易受到自身時頻特性以及復雜太空環(huán)境的影響,高可靠的厘米級GNSS實時衛(wèi)星鐘差則需利用地面站觀測數(shù)據(jù)進行實時估計得到[7]。文獻[8]實時估計的多模GNSS衛(wèi)星鐘差與武漢大學最終精密鐘差互差優(yōu)于0.2 ns。文獻[9-10]利用均方根信息濾波實現(xiàn)GNSS四系統(tǒng)實時衛(wèi)星鐘差的快速估計,且與最終精密鐘差產(chǎn)品符合性較好。文獻[11]提出了一種基于序貫最小二乘的在線質量控制方法,GNSS四系統(tǒng)實時估計衛(wèi)星鐘差的標準差均值優(yōu)于0.1 ns。

鑒于此,本文以長安大學北斗分析中心為平臺,基于多模全球導航衛(wèi)星系統(tǒng)實驗跟蹤網(wǎng)(multi-GNSS experiment, MGEX)監(jiān)測站1 Hz的GNSS觀測數(shù)據(jù),采用模型嚴密、精度較高的非差估計法[12],實現(xiàn)了多模GNSS衛(wèi)星鐘差的實時估計和性能評估,并將其用于高頻動態(tài)PPP實時獲取2021年漾濞Mw6.4地震和瑪多Mw7.4地震發(fā)生時段的地表形變波形,具體分析震時點位的運動變化情況。

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

1.1 鐘差估計

本文采用無電離層組合觀測值進行GNSS衛(wèi)星鐘差的實時估計,非差偽距和載波相位的觀測方程為

(1)

由于公式中的偏差參數(shù)具有很強甚至完全的線性相關性,將其作為未知數(shù)進行估計,不僅會引入大量待估參數(shù),還會減弱衛(wèi)星鐘差的估計精度。其中,衛(wèi)星端的碼偏差參數(shù)具有較高的時間穩(wěn)定性,可以在觀測方程中將其合并到衛(wèi)星的鐘誤差參數(shù)中,而且相位延遲偏差可以被相應的相位模糊度參數(shù)吸收。進而可得觀測模型

(2)

對于全球定位系統(tǒng)(global positioning system, GPS)、格洛納斯衛(wèi)星導航系統(tǒng)(global navigation satellite system, GLONASS)、伽利略衛(wèi)星導航系統(tǒng)(Galileo satellite navigation system, Galileo)以及我國的北斗衛(wèi)星導航系統(tǒng)(BeiDou navigation satellite system, BDS)而言,不同衛(wèi)星導航系統(tǒng)所采用的時空基準以及信號體制不一致,導致衛(wèi)星信號在GNSS接收機內(nèi)部產(chǎn)生系統(tǒng)間偏差(inter system bias, ISB)以及GLONASS所特有的頻率間偏差(inter frequency bias, IFB)[11]。若以GPS時間系統(tǒng)為基準,則顧及ISB/IFB參數(shù)的GNSS衛(wèi)星鐘差實時估計模型為

(3)

1.2 處理策略

無電離層組合模型是GNSS精密衛(wèi)星鐘差估計、精密單點定位等數(shù)據(jù)處理常用的函數(shù)模型,表1總結了GNSS衛(wèi)星鐘差實時估計以及實時動態(tài)PPP基于該模型的數(shù)據(jù)處理策略。對于GNSS衛(wèi)星鐘差的實時估計,則是在觀測方程中將衛(wèi)星鐘誤差作為白噪聲進行估計以免受到鐘跳的影響,將衛(wèi)星軌道和測站坐標作為已知值進行改正。對于GNSS實時動態(tài)PPP,則是將測站的位置坐標作為待估參數(shù)進行估計,將精密衛(wèi)星軌道和衛(wèi)星鐘差作為已知值進行改正。

表1 數(shù)據(jù)處理策略Tab.1 Data processing strategy

2 多模GNSS實時衛(wèi)星鐘差估計分析

在MGEX中選取60個均勻分布在全球的連續(xù)運行跟蹤站。從IGS下載2021年5月19日—2021年5月22日連續(xù)4 d所選測站1 s采樣間隔的GNSS觀測數(shù)據(jù)用于衛(wèi)星鐘差的實時估計。在實時估計多模GNSS衛(wèi)星鐘差的過程中,測站始終處于靜止狀態(tài)且位置已知,可以將測站坐標固定到IGS周解?，F(xiàn)階段,超快速實時預報軌道產(chǎn)品與最終精密軌道產(chǎn)品的精度量級相當,衛(wèi)星軌道可以固定為武漢大學提供的包含G/R/E/C四系統(tǒng)軌道的6 h超快速解,以減少待估參數(shù)的個數(shù)。同時,選擇某一測站的接收機鐘差作為參考基準鐘進行約束[15],基于雙頻非差消電離層組合的載波相位和偽距觀測值,采用序貫最小二乘平差和驗后殘差質量控制算法,最終估計得到歷元間隔為5 s的多模GNSS實時衛(wèi)星鐘差。

2.1 估計時間分析

當利用60個測站構成的地面跟蹤站網(wǎng)進行多模GNSS衛(wèi)星鐘差實時估計時,每個歷元大約需要處理3 000個觀測值并對將近2 000維的矩陣進行求逆。每個歷元需要處理的觀測值數(shù)目和法方程維數(shù)如圖1所示。為加快高維矩陣的運算速度,提高鐘差估計算法的計算效率,本文采用LAPACK函數(shù)庫中相關的矩陣運算解算法方程,以實現(xiàn)多模GNSS實時衛(wèi)星鐘差的快速估計[16]。圖2統(tǒng)計了LAPACK方法下,實時估計2021年年積日第141天G/R/E/C四系統(tǒng)衛(wèi)星鐘差時每個歷元的計算時間?？梢缘贸?在全球分布60個連續(xù)運行跟蹤站的情況下,所有歷元的鐘差估計耗時均小于5 s,單歷元鐘差估計的平均計算時間為3.4 s,可見上述算法對高頻觀測數(shù)據(jù)的處理效率足以實時估計5 s歷元間隔的多模GNSS衛(wèi)星鐘差。

圖1 每個歷元的觀測值數(shù)目和法方程維數(shù)Fig.1 The number of observations and the dimension of the normal equation system at each epoch

圖2 多模GNSS實時鐘差單歷元估計時間Fig.2 The computation time for the multi-GNSS real-time clock offset estimation at each epoch

2.2 鐘差精度分析

采用二次差的方法計算實時估計衛(wèi)星鐘差與武漢大學最終精密衛(wèi)星鐘差的差異[17],利用該差異統(tǒng)計其標準差(standard deviation, STD),并與法國空間研究中心(Centre National D’études Spatiales,CNES)實時播發(fā)的衛(wèi)星鐘差產(chǎn)品進行對比,對多模GNSS實時衛(wèi)星鐘差估計算法的有效性進行評估?？紤]到衛(wèi)星鐘差實時估計的收斂時間,僅對2021年5月20日—2021年5月22日(年積日第140—142天)的實時衛(wèi)星鐘差進行精度分析,并在圖3中展示了單顆衛(wèi)星連續(xù)3 d的鐘誤差時間序列,BDS衛(wèi)星鐘差的整體穩(wěn)定性要弱于其他導航衛(wèi)星系統(tǒng)。同時,GNSS實時估計鐘差與CNES實時鐘差產(chǎn)品的精度對比如表2和圖4所示。其中,表2為GPS、GLONASS、Galileo、BDS各系統(tǒng)單天的平均STD對比,圖4為GNSS每顆衛(wèi)星連續(xù)3 d的平均STD對比。

圖3 GPS/GLONASS/Galileo/BDS衛(wèi)星鐘誤差的時間序列Fig.3 Time series of GPS/GLONASS/Galileo/BDS satellite clock error

圖4 GPS/GLONASS/Galileo/BDS衛(wèi)星鐘差精度對比Fig.4 Accuracy comparison of GPS/GLONASS/Galileo/BDS satellite clock

表2 衛(wèi)星鐘差精度統(tǒng)計Tab.2 Accuracy statistics of satellites clock offset ns

由于衛(wèi)星鐘差與軌道誤差的耦合性,在實時估計多模GNSS衛(wèi)星鐘差時,90%的軌道徑向誤差會被鐘差所吸收,從而導致實時估計鐘差和最終精密鐘差兩者之間并不完全吻合[18]。結合表2和圖4可以得出,實時估計GPS衛(wèi)星鐘差的STD數(shù)值區(qū)間為0.044～0.094 ns,每顆GPS衛(wèi)星的鐘差精度均優(yōu)于0.1 ns,整體穩(wěn)定性較好,所有衛(wèi)星的STD均值為0.064 ns。除G18衛(wèi)星的鐘差STD值為0.215 ns外,CNES實時播發(fā)的GPS衛(wèi)星鐘差產(chǎn)品與自主估計的GPS實時鐘差精度基本一致,這可能與G18衛(wèi)星的型號為Block Ⅲ,發(fā)射時間較晚,CNES的軌道模型還未對其完全精化有關。對于Galileo系統(tǒng)而言,實時估計衛(wèi)星鐘差與CNES實時衛(wèi)星鐘差二次差結果序列相差不大,所有Galileo衛(wèi)星的STD均值分別為0.058 ns和0.080 ns,且每顆衛(wèi)星的鐘差精度基本上均能優(yōu)于0.1 ns,能夠達到與GPS衛(wèi)星鐘差相當?shù)木人?這可能與Galileo衛(wèi)星所搭載的高精度氫原子鐘有關。CNES實時播發(fā)的GLONASS鐘差產(chǎn)品中,R09衛(wèi)星連續(xù)3 d的鐘差精度均將達到1.5 ns,這是因為該衛(wèi)星的鐘差序列出現(xiàn)多次中斷,導致精度異常。對于其他GLONASS衛(wèi)星,二次差結果序列的STD數(shù)值區(qū)間為0.087～0.598 ns,STD均值為0.232 ns,稍差于自主估計的GLONASS實時衛(wèi)星鐘差精度,所有衛(wèi)星的STD均值為0.163 ns。對于CNES實時BDS衛(wèi)星鐘差產(chǎn)品,地球靜止軌道(geostationary Earth orbit,GEO)衛(wèi)星的鐘差精度明顯差于中圓地球軌道(medium Earth orbit,MEO)衛(wèi)星和傾斜地球同步軌道(inclined geosynchronous orbit,IGSO)衛(wèi)星,導致BDS衛(wèi)星鐘差整體精度較差。除GEO衛(wèi)星外,其余MEO/IGSO衛(wèi)星的STD均值為0.294 ns,大約是自主估計BDS MEO/IGSO實時衛(wèi)星鐘差STD均值的2倍。從整體上看,自主估計的GPS、GLONASS、Galileo、BDS GEO、BDS IGSO和BDS MEO實時衛(wèi)星鐘差STD均值分別為0.064、0.163、0.058、0.232、0.154和0.179 ns,要稍優(yōu)于CNES實時播發(fā)的衛(wèi)星鐘差產(chǎn)品,兩者G/R/E/C四系統(tǒng)的STD均值分別為0.142 ns和0.347 ns,這是由于在實時估計GNSS衛(wèi)星鐘差時,二者選取測站的數(shù)據(jù)質量和解算策略存在差異所致。

3 震時地表形變分析

3.1 2021年云南漾濞Mw 6.4地震

為進一步評估實時估計多模GNSS衛(wèi)星鐘差的精度和穩(wěn)定性,利用中國大陸構造環(huán)境監(jiān)測網(wǎng)絡(簡稱“陸態(tài)網(wǎng)絡”)基準站2021年5月21日21:45—21:59時段內(nèi)的高頻GNSS觀測數(shù)據(jù),基準站分布情況如圖5所示,采用表1中實時動態(tài)PPP的數(shù)據(jù)處理策略,基于5 s歷元間隔的GNSS實時衛(wèi)星鐘差,在實時處理的定位模式下,通過序貫最小二乘估計方法,歷元之間不繼承坐標信息,分別進行GPS單系統(tǒng)、GNSS多系統(tǒng)的高頻動態(tài)PPP,并與武漢大學最終精密衛(wèi)星鐘差的GNSS多系統(tǒng)組合動態(tài)PPP結果進行對比(表3),從而得到不同實驗下2021年云南漾濞Mw6.4地震發(fā)震時段內(nèi)單歷元高頻動態(tài)解的時間序列。上述3種高頻動態(tài)PPP實驗分別對應實驗1、實驗2和實驗3,在不同定位實驗下,云南下關站(XIAG)、云南云龍站(YNYL)和云南麗江站(YNLJ)在東(East,E)、北(North,N)和高程(Up,U)方向上的動態(tài)位移時間序列如圖6所示,圖中垂直于橫軸的紅色實線為漾濞地震的發(fā)震時刻。

圖5 漾濞地震震中周邊GNSS基準站分布Fig.5 Distribution of GNSS reference station around the epicenter of Yangbi earthquake

圖6 漾濞地震發(fā)生時段的動態(tài)位移時間序列Fig.6 The time series of dynamic displacement during the Yangbi earthquake

表3 不同測站高頻動態(tài)PPP的標準差Tab.3 Standard deviation of high-rate kinematic PPP at different stations

由圖6和表3可見,隨著其他系統(tǒng)的加入,基于實時估計衛(wèi)星鐘差的GNSS多系統(tǒng)融合精密單點定位的動態(tài)位移波形更為穩(wěn)定,其高頻動態(tài)解在E、N、U方向上的平均標準差分別為0.4、0.3、1.0 cm,相對于GPS單系統(tǒng)動態(tài)定位結果的平均標準差,在三個方向上分別提升了33%、25%和50%,印證了多系統(tǒng)融合定位的優(yōu)越性[19]。同時,自主估計的GNSS實時衛(wèi)星鐘差與武漢大學最終精密衛(wèi)星鐘差的定位性能相當,其GNSS多系統(tǒng)組合PPP動態(tài)解的平均標準差能夠達到水平方向0.5 cm左右,垂直方向1.0 cm左右的精度水平。對于2021年5月21日云南漾濞地震,GNSS高頻動態(tài)解可以觀測到距震中39.7 km的XIAG站在水平方向上具有相對明顯的波動,這與文獻[20]漾濞地震的水平位移主要發(fā)生在震中距50 km范圍內(nèi)相吻合?？梢?基于上述數(shù)學模型估計的多模GNSS衛(wèi)星鐘差能夠用于實時監(jiān)測地震發(fā)生階段地表的動態(tài)形變特征。

3.2 2021年青?，敹郙w 7.4地震

據(jù)中國地震臺網(wǎng)正式測定:北京時間2021年5月22日2時4分,青海省果洛藏族自治州瑪多縣(震中34.59°N,98.34°E)發(fā)生7.4級地震,震源深度17 km[21]。此次地震是中國大陸地區(qū)內(nèi)繼2008年汶川Mw8.0地震后發(fā)生的震級最大的地震,其劇烈的地殼形變信號能夠被區(qū)域內(nèi)連續(xù)運行的GNSS基準站可靠監(jiān)測[22]。為具體分析瑪多Mw7.4地震發(fā)生階段對周圍地表形變的影響,利用中國地震臺網(wǎng)中心提供的震源周邊地區(qū)10個陸態(tài)網(wǎng)絡基準站的高頻GNSS觀測數(shù)據(jù),這10個基準站分別為青?，敹嗾?QHMD)、青?，斍哒?QHMQ)、青海都蘭站(QHDL)、青海班瑪站(QHBM)、青海德令哈站(DLHA)、甘肅瑪曲站(GSMA)、四川甘孜站(SCGZ)、青海格爾木站(QHGE)、青海西寧站(XNIN)和西藏昌都站(XZCD),基準站分布情況如圖7所示。以地震發(fā)生前3 d的精密坐標的平均值為參考值,基于自主估計的2021年年積日第141天實時衛(wèi)星鐘差進行GNSS多系統(tǒng)組合實時動態(tài)PPP,獲取瑪多地震震后5 min時段內(nèi)不同基準站高頻動態(tài)解在E、N、U方向上隨時間變化的位移序列,如圖8所示。同時,對上述10個基準站震前3 d以及震后1 d的低頻GNSS觀測數(shù)據(jù)進行事后靜態(tài)PPP處理得到精密坐標的單天解,然后分別計算3個坐標分量上地震前后基準站位置的坐標平均值,通過差分獲得瑪多地震高可靠性的同震形變[23],并與GNSS多系統(tǒng)組合實時動態(tài)PPP獲得的同震形變進行對比,具體形變結果如表4所示。

圖7 瑪多地震震中周邊GNSS基準站分布Fig.7 Distribution of GNSS reference stations around the epicenter of Madoi earthquake

圖8 不同測站動態(tài)位移時間序列Fig.8 The time series of dynamic displacement at different stations

表4 基于高頻和低頻GNSS觀測的地表形變Tab.4 The surface deformation from 1 s and 30 s sampling GNSS observation data

由圖8可見,當瑪多地震發(fā)生后,10個陸態(tài)網(wǎng)絡基準站在不同時刻均發(fā)生不同程度的波動,且開始波動時刻隨基準站震中距的增大而延遲,影響區(qū)域可達距震中約400 km。其中,距震中約39 km的陸態(tài)網(wǎng)絡GNSS基準站QHMD最先感知到地震波,受瑪多地震的影響也最大,其位置在東西向和南北向均發(fā)生不可逆的永久性位移,東西向為23.3 cm,南北向為8.5 cm。不同于東西方向上的動態(tài)位移特征,QHMQ和GSMA相對于其他基準站在南北向受到瑪多地震的影響最為顯著,其震時最大動態(tài)位移分別約為23.5 cm和15.6 cm,說明地震對震源周圍站點的影響不僅和震中距有關,與其方位角也有一定的相關性[24]。然而,瑪多地震在垂直方向上對震源周邊基準站的影響相對較弱,10個基準站均未出現(xiàn)顯著的位移形變。從表4可以得出,受益于GNSS多系統(tǒng)融合具有更為穩(wěn)健的定位精度,其實時PPP高頻動態(tài)解能夠有效地獲取地震發(fā)生階段各測站在水平方向上較為準確的地表形變量,與低頻GNSS觀測的形變結果基本一致,而兩者在垂直方向上的差值最大為1.1 cm,可靠性相對較差。綜上所述,基于多模GNSS實時衛(wèi)星鐘差的PPP高頻動態(tài)解可以有效地揭示瑪多地震的震時地表形變特征。

4 結論

本文采用MGEX中60個連續(xù)跟蹤站的1 Hz觀測數(shù)據(jù),基于序貫最小二乘算法,實現(xiàn)了多模GNSS衛(wèi)星鐘差的實時估計,并將其用于GNSS多系統(tǒng)組合的高頻動態(tài)PPP,具體分析了2021年云南漾濞Mw6.4地震和瑪多Mw7.4地震的震時地表形變特征,得到以下結論:

1)自主估計的GPS、GLONASS、Galileo、BDS GEO、BDS IGSO和BDS MEO實時衛(wèi)星鐘差STD均值分別為0.064、0.163、0.058、0.232、0.154和0.179 ns,整體上稍優(yōu)于CNES實時播發(fā)的衛(wèi)星鐘差產(chǎn)品。

2)基于實時估計衛(wèi)星鐘差的GNSS多系統(tǒng)組合動態(tài)PPP在E、N、U三個方向上的平均標準差分別為0.4、0.3、1.0 cm,相對于GPS單系統(tǒng)的定位結果分別提升了33%、25%和50%,能夠提供波形更為穩(wěn)定的動態(tài)位移時間序列。

3)GNSS多系統(tǒng)融合的實時PPP高頻動態(tài)解能夠快速有效地獲取地震發(fā)生階段站點的形變特征,為斷層運動特征的初步判定以及地震矩震級的快速確定提供較高可靠性的信息。