基于時間序列的GPS 速率可靠性分析*

丁曉光 甘衛(wèi)軍 肖根如 張 藝

1)地震動力學國家重點實驗室，中國地震局地質(zhì)研究所，北京 100029

2)陜西省地震局，西安 710068

3)東華理工大學，撫州344000

1 引言

中國地殼運動觀測網(wǎng)絡”(簡稱“網(wǎng)絡工程”)［1］和“中國大陸構(gòu)造環(huán)境監(jiān)測網(wǎng)絡”(簡稱“陸態(tài)網(wǎng)絡”)［2］的GPS 觀測站既有長年不間斷連續(xù)觀測的基準站，也有每隔1 ～2年臨時連續(xù)觀測3 ～4天的流動站。截至2013年，網(wǎng)絡工程的27 個GPS 基準站已積累了約14年跨度的連續(xù)觀測數(shù)據(jù)，1 055 個流動站至少已有7 ～8 期的重復觀測;而陸態(tài)網(wǎng)絡建成于2012年，233 個新建連續(xù)站僅有2 ～3年的觀測時長，1 000 個新建流動站只有2009年和2011年兩期觀測。由于兩個項目的GPS 觀測時間跨度和重復觀測期次相差懸殊，那么基于時間跨度較短的連續(xù)GPS 觀測和非連續(xù)GPS 觀測所得到的站點速度可靠性如何?該問題也即，要獲取可靠有效的GPS 站點速度，連續(xù)站需要累積多長時間的觀測，流動站需要幾期的觀測和多長的時間跨度?針對這個疑問，本文將以“網(wǎng)絡工程”基準站長時間跨度的連續(xù)GPS 觀測數(shù)據(jù)和速度結(jié)果為基礎，采用模擬觀測(數(shù)據(jù)分段模擬連續(xù)觀測，以及隨機抽取連續(xù)數(shù)據(jù)模擬流動觀測)和數(shù)值計算方法給出了全面的統(tǒng)計分析。

2 連續(xù)GPS 數(shù)據(jù)及數(shù)據(jù)處理

利用GIPSY-OASIS(Ver.6.0)軟件［3，4］精密單點定位模式對“網(wǎng)絡工程”站點1999—2012年的觀測數(shù)據(jù)進行計算，獲得單日松弛約束解，然后采用聯(lián)合平差軟件QOCA［5］對所有站點的單日松弛約束解嚴密平差，獲得各站點的坐標變化序列。數(shù)據(jù)處理中具體要點為:以消除電離層影響的線性組合作為觀測量，衛(wèi)星截止高度角為15°，應用JPL 精密星歷和時鐘產(chǎn)品并估計鐘差參數(shù);先驗對流層延遲模型采用全球大氣壓溫度模型(GPT，Global Pressure and Temperature)，投影函數(shù)為GMF，再用參數(shù)估計對流層延遲剩余部分;海洋潮汐改正由FES2004 模型基于格林函數(shù)在線計算獲得;進行接收機天線和衛(wèi)星天線絕對相位中心改正;采用Ambizap 模塊進行整周模糊度解算［6］，該模塊的處理算法是利用固定點法則來確定觀測網(wǎng)參數(shù)的各種線性組合，并為整個觀測網(wǎng)生成消除整周模糊度且唯一、自洽的單日解;最后利用QOCA 軟件與全球271 個IGS 站的單日松弛約束解進行聯(lián)合平差，通過選取的IGS 站進行7參數(shù)變換將單日解統(tǒng)一轉(zhuǎn)換到ITRF2008 框架下。

從解算得到的各站點時間序列中選取數(shù)據(jù)較為完整、結(jié)果較為穩(wěn)定的XNIN(西寧)站作為實驗對象。通過QOCA 軟件analys_timseri 模塊對原始時間序列進行分析，刪除個別偏離線性過大的突跳點，去除2001年11月14日昆侖山地震的同震影響，以及其他幾次由于更換天線造成的大幅階躍(圖1)。為便于后續(xù)研究，截取了時序年變穩(wěn)定的一段數(shù)據(jù)進行分析(2001年年積日100天至2010年年積日220天)，利用最小二乘算法求得該段數(shù)據(jù)N、E、U分量的年平均速率分別為:－3.80 mm/a、37.94 mm/a 和0.02 mm/a。

3 GPS 速率模擬計算與可靠性分析

3.1 連續(xù)站速率可靠性分析

為研究獲得穩(wěn)定可靠的速率所需最短觀測時長，采用了以下計算方法:以第1天至360天數(shù)據(jù)得到的速率值為起始點，每增加1天，重新計算一次速率，這樣便得到了361，362，…直至全部數(shù)據(jù)得到的速率，將其連成一條曲線，當曲線波動小于給定閾值時，即表示速率值穩(wěn)定。再以第10天至370天的速率為起始點，重復以上計算，又會得到一條曲線…以此類推，分析起算時間位于坐標序列不同相位時的速率變化。

圖1 XNIN 站坐標時間序列Fig.1 Coordinate time series of XNIN station

利用QOCA 軟件時序分析模塊計算XNIN 站時間序列N、E、U 分量的周年項相位，分別為347.83°、224.39°、316.76°，與田云鋒［7］計算結(jié)果(351.69°、215.15°、337.68°)基本相同。由此推算N 方向的0 相位處于每年年初10 ～15天。圖2 為N 分量起算速率的相位大致位于0°、90°、180°、270°時的速率變化曲線。

從圖2 來看，水平分量速率曲線均呈現(xiàn)振幅遞減的波動狀態(tài)，振幅約1 ～3 mm，在持續(xù)3 ～4年后，波動變化小于0.5 mm，速率趨于穩(wěn)定;在不同初始相位時，各曲線變化趨勢各不相同，N 分量變化范圍在－1 ～－5 mm 之間，E 分量則達到35 ～41 mm。垂向分量由于振幅較大，而線性速率很小，故擬合速率呈直線變化，受突跳點影響出現(xiàn)毛刺，未得到可用于分析的結(jié)果。

速率曲線上下波動逐漸趨近于穩(wěn)定值，這一現(xiàn)象說明，在計算所用數(shù)據(jù)量不能使速率穩(wěn)定之前，增加觀測天數(shù)未必能得到更真實的速率值;要獲得穩(wěn)定可靠的結(jié)果需3 ～4年的觀測時長。值得指出的是，高精度GPS 軟件在處理長時間觀測數(shù)據(jù)時，除對測站坐標線性擬合出速率外，還要在穩(wěn)定參考框架下進行。與已知準確速率的測站組綁定一起平差，能夠在相對較短的時間內(nèi)得到穩(wěn)定的速率，這與本文討論的利用時間序列擬合獲得速率的方法有很大不同，關(guān)于流動站速率的可靠性的分析與之類似。

3.2 流動站速率可靠性分析

首先對XNIN 站時間序列中所有坐標點進行加權(quán)滑動平均，滑動窗口為4天。對原坐標序列所有數(shù)據(jù)進行滑動平均后，得到了新的“流動GPS”坐標時間序列(簡稱“新序列”)，每一時間點均對應4天流動測量。

圖2 以N 分量坐標序列不同相位起算得到的速率曲線Fig.2 The rate curves calculated by different phases of coordinate series of component N

以新序列中第1 組數(shù)據(jù)為一期觀測，間隔N年后再選1 組數(shù)據(jù)作為第二期觀測，每期觀測時間點的選取方法為:

Δt 為在1 至100 間隨機抽取的整數(shù)，模擬每期觀測時間相差的天數(shù)。在自動選取得到多期數(shù)據(jù)之后，擬合獲得第一組速率值。然后再以新序列中第2 組數(shù)據(jù)為第一期觀測，按照同樣方法選取多期模擬觀測值，擬合得到第二組速率值…以此類推，可以得到多組速率值，約定每組速率值與各自第一期的觀測時間形成對應關(guān)系，構(gòu)成一條速率-時間曲線。

根據(jù)每期觀測時間間隔和總時長的不同，設計了三種流動模擬觀測方案。

方案一:每1年觀測1 期。觀測總時長為3、4、5、6、7 期;

方案二:每2年觀測1 期。觀測總時長為2、3、4 期;

方案三:每3年觀測1 期。觀測總時長為2、3期。

每種方案都計算了新序列的前720 組數(shù)據(jù)，從而得到720 個速率結(jié)果組成的序列。將各方案得到的速率序列與擬合的連續(xù)站速率相比較，計算其差異的標準差，作為外符合精度;再與本速率序列的平均值比較，計算差異的標準差作為內(nèi)符合精度(表1)。

由表1 可見，在觀測周期相同時，隨觀測期數(shù)的不斷增加，內(nèi)、外符合精度均不斷減小;在時間跨度相同時，觀測期數(shù)不同，得到的內(nèi)、外符合精度基本相當(圖3 ～5)。

根據(jù)圖3 ～5 并結(jié)合表1，3 種方案得到的速率可靠性相近，當然前提是每期觀測數(shù)據(jù)質(zhì)量都能得到保證，實際上，增加觀測期數(shù)來提高結(jié)果的可信度是必要的，如每期測量時采用不同型號天線、天線高量測、拆裝整平天線等都會影響解算結(jié)果。由圖3 ～5 還可以看出，在速率曲線的三個分量上均存在年周期變化，并且波動幅度隨時間跨度的增長而減小。

結(jié)合表1 所列數(shù)據(jù)，雖然在時間跨度達到5 ～6年時，速率序列的水平分量標準差提高到0.5 mm左右，垂向標準差約為1 mm，但真實的流動測量只可能得到序列中的某一個速率值，由于年周期變化的存在，其計算結(jié)果會大于標準差。為了對流動觀測速率可靠性有完整的認識，表2 列出了各方案速率序列與已知的連續(xù)站速率做差后，差異值的變化范圍。

表1 各方案得到的速率序列的標準差(單位:mm/a)Tab.1 Standard deviations of rate series calculated by different schemes(unit:mm/a)

表2 各方案得到的速率序列與連續(xù)站速率值的差異(單位:mm/a)Tab.2 Differences between the rates calculated by different schemes and the rates directly fitted the data of continuous station(unit:mm/a)

圖3 時間跨度為兩年的速率曲線比較Fig.3 Comparation of the velocity curves in 2 years

圖4 時間跨度為3年的速率曲線比較Fig.4 Comparation of the rate curves in 3 years

式中x0為起算數(shù)據(jù)，v 為擬合線性項(速率)，Si、Ci為時間序列諧波分量振幅，主要為周年和半年項周期變化。Hk為第k 次階躍，Bk表示該次階躍對后續(xù)時間序列影響的函數(shù)，Rj為剩余的殘差項。

同樣XNIN 站的時間序列為例，扣除階躍和強震影響，不考慮第四項;則周期項可以表示為振幅和相位的函數(shù)，即Aisin(ωit + φi)。其中，振幅Ai=，相位

流動觀測中兩期觀測時刻分別為ti和tj，相應坐標為x(ti)和x(tj)，則速度

將式(1)代入式(2)，令ti=0，考慮同一序列不同時刻振幅頻率相等，得

圖5 時間跨度為6年的速率曲線比較Fig.5 Comparation of the velocity curves in 6 years

從表2 中可以看到，當時間跨度僅為2年時，得到的流動站速率與真實值的偏差在水平方向可以達到±3 ～4 mm/a，而垂向分量的速率波動接近±8 mm/a;當時間跨度為4年時，水平方向速率的可靠性在±3 mm/a 以內(nèi)，垂向分量的波動水平在±5 mm/a 以內(nèi);即使時間跨度達到6年，水平速率仍可能有超過±1 mm/a 的偏差，垂向的偏差則可能有±3 mm/a。如果再考慮實際觀測時的復雜影響因素，流動GPS 速率的可靠性還要降低。

4 速率可靠性的影響因素與討論

時間序列中任意時刻t 點均可表示為

當初始相位φi已知時，第二項為正弦函數(shù)。由式(3)可知，流動觀測中影響速率估計的因素有相位差、振幅和殘差，其中相位差由兩期觀測的時間差異決定;振幅取決于站點固有頻率周期;而殘差則決定于站點觀測噪聲的大小，它們都能夠隨觀測間隔t的增加而減少;在相同時間間隔內(nèi)加密觀測期數(shù)能夠提高速率v 的擬合精度，卻無法消除周期誤差和殘差的影響，速率的準確性未得到有效提高，例如本文中利用連續(xù)站(相當于每天1 期流動觀測)得到的水平速率結(jié)果，在觀測時長為3年時仍有1 mm左右波動，直到時長為4 ～5年時速率值才趨于穩(wěn)定。當然，根據(jù)連續(xù)站時間序列估計周期項參數(shù)和噪聲類型，可以縮短獲得可靠速率值所用時間，但是要得到較為準確的周期和噪聲模型，也至少需3年有效觀測。

與連續(xù)站相比，流動站僅靠自身每期數(shù)天觀測無法獲知周期振幅和殘差模型，此時只有通過縮小兩期觀測時間差異、增加觀測時間跨度來提高速率結(jié)果的可靠性。另外，隨著連續(xù)站布設密度的逐步提高，流動站的周期性變化模型可以由鄰近連續(xù)站擬合近似得出，從而達到改善流動站速率可靠性的目的。

5 結(jié)論

1)針對XNIN 連續(xù)站，通過擬合時間序列獲得的站點速率，需要至少3 ～4年數(shù)據(jù)累積才能獲得穩(wěn)定可靠的結(jié)果(＜±0.5 mm/a)。當數(shù)據(jù)積累時間較短時，計算速率的起始點位于時間序列的不同相位，會引起不同的速率變化趨勢。

2)利用XNIN 連續(xù)站模擬的流動觀測，當時間跨度為2年時，速率的水平分量波動范圍可達±3 ～4 mm/a;時間跨度達到6年時，速率水平分量可靠性接近±1 mm/a，垂向分量約為±3 mm/a。如果考慮實際觀測時每期儀器更換等因素影響，可靠性還要低于估計值。

3)對于所有連續(xù)站和流動站，獲得穩(wěn)定可靠的速率所需時間長短受站點固有周期和噪聲水平影響。

4)增加觀測時間跨度，縮小流動觀測兩期觀測時間差異可以提高速率結(jié)果的可靠性。

1 甘衛(wèi)軍，等.中國地殼運動觀測網(wǎng)絡的建設及應用［J］.國際地震動態(tài)，2007，343(7):43－52.(Gan weijun，et al.Construction and application of crust movement observation network in China Continant［J］.Development of Internanional Earthquake，2007，343(7):43－52)

2 甘衛(wèi)軍，等.中國大陸構(gòu)造環(huán)境監(jiān)測網(wǎng)絡的建設與應用［J］.工程研究，2012，4(4):324－331.(Gan Weijun，et al.Construction and application of tectonic environment observation network in China Continant［J］.Engineering Studies，2012，4(4):324－331)

3 肖根如，甘衛(wèi)軍，殷海濤.GIPSY 軟件的GPS 數(shù)據(jù)處理策略及應用［J］.地球物理學進展，2010，25(4):1 508－1 515.(Xiao Genru，Gan Weijun and Yin Haitao.GPS data process strategies and application of GIPSY software［J］.Progress in Geophys，2010，25(4):1 508－1 515)

4 Gregorius T.GIPSY-OASIS II:How it Works，(self-published)，Univ.of Newcastle upon Tyne，Newcastle，England，U.k.，1996.

5 Dong Danan，Herring T A and King R W.Estimating regional deformation from a combination of space and terrestrical geodetic data［J］.J Geodesy.，1998，72:200－214.

6 Geoffrey Blewitt.Fixed point theorems of GPS carrier phase ambiguity resolution and their application to massive network processing:Ambizap［J］.Journal of Geophysical Reasearch.2008，113，B12410.

7 田云鋒.GPS 位置時間序列中的中長期誤差研究［D］.中國地震局地質(zhì)研究所，2011.(Tian Yunfeng.Study on intermediate-and long-term errors in GPS position time series［D］.Institute of Geology，CEA，2011)

8 占偉，等.GNSS 流動觀測水平速度精度評估［J］.大地測量與地球動力學，2011，(5):84－87.(Zhan Wei，et al.Accuracy evaluation of horizontal velocity in GNSS mobile observation［J］.Journal of Geodesy and Geodynamics，2011，(5):84－87)

9 王敏，等.非構(gòu)造形變對GPS 連續(xù)站位置時間序列的影響和修正［J］.地球物理學報，2005，48(5):1 045－1 052.(Wang Min，et al.Effects of non-tectonic crustal deformation on continuous GPS position time series and correction［J］.Chinese J Geophys.，2005，48(5):1 045－1 052)