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

1)地震動(dòng)力學(xué)國家重點(diǎn)實(shí)驗(yàn)室，中國地震局地質(zhì)研究所，北京 100029

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

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

1 引言

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

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

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

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

3 GPS 速率模擬計(jì)算與可靠性分析

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

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

圖1 XNIN 站坐標(biāo)時(shí)間序列Fig.1 Coordinate time series of XNIN station

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

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

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

3.2 流動(dòng)站速率可靠性分析

首先對(duì)XNIN 站時(shí)間序列中所有坐標(biāo)點(diǎn)進(jìn)行加權(quán)滑動(dòng)平均，滑動(dòng)窗口為4天。對(duì)原坐標(biāo)序列所有數(shù)據(jù)進(jìn)行滑動(dòng)平均后，得到了新的“流動(dòng)GPS”坐標(biāo)時(shí)間序列(簡稱“新序列”)，每一時(shí)間點(diǎn)均對(duì)應(yīng)4天流動(dòng)測(cè)量。

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

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

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

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

方案一:每1年觀測(cè)1 期。觀測(cè)總時(shí)長為3、4、5、6、7 期;

方案二:每2年觀測(cè)1 期。觀測(cè)總時(shí)長為2、3、4 期;

方案三:每3年觀測(cè)1 期。觀測(cè)總時(shí)長為2、3期。

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

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

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

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

表1 各方案得到的速率序列的標(biāo)準(zhǔn)差(單位: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 時(shí)間跨度為兩年的速率曲線比較Fig.3 Comparation of the velocity curves in 2 years

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

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

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

流動(dòng)觀測(cè)中兩期觀測(cè)時(shí)刻分別為ti和tj，相應(yīng)坐標(biāo)為x(ti)和x(tj)，則速度

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

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

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

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

時(shí)間序列中任意時(shí)刻t 點(diǎn)均可表示為

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

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

5 結(jié)論

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

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

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

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

1 甘衛(wèi)軍，等.中國地殼運(yùn)動(dòng)觀測(cè)網(wǎng)絡(luò)的建設(shè)及應(yīng)用［J］.國際地震動(dòng)態(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)測(cè)網(wǎng)絡(luò)的建設(shè)與應(yī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ù)處理策略及應(yīng)用［J］.地球物理學(xué)進(jìn)展，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 位置時(shí)間序列中的中長期誤差研究［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 流動(dòng)觀測(cè)水平速度精度評(píng)估［J］.大地測(cè)量與地球動(dòng)力學(xué)，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)造形變對(duì)GPS 連續(xù)站位置時(shí)間序列的影響和修正［J］.地球物理學(xué)報(bào)，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)