高頻GPS雙差殘差模型監(jiān)測強(qiáng)震地表運(yùn)動(dòng)

馮 威，黃丁發(fā)，李 萌，張 熙，嚴(yán) 麗

西南交通大學(xué)測量工程系，成都 610031

1 引 言

隨著GPS接收機(jī)技術(shù)的發(fā)展和高頻接收機(jī)的普及，高頻GPS技術(shù)在強(qiáng)震領(lǐng)域得到了越來越廣泛的研究和應(yīng)用［1－5］.且與傳統(tǒng)地震儀相比，由于不受振幅限制，GPS能夠有效監(jiān)測大幅度的地面震動(dòng)，其可為地震預(yù)警、救援和相關(guān)理論研究提供可靠的信息.

為了滿足強(qiáng)震地表運(yùn)動(dòng)監(jiān)測的高精度定位和單歷元解算的要求，常用精密單點(diǎn)定位（PPP）模式和差分定位模式來實(shí)現(xiàn)高精度的位移解算.由于地震影響的范圍較大，采用PPP技術(shù)進(jìn)行地震的地表位移解算［6］具有一定的優(yōu)勢.但由于PPP技術(shù)對(duì)衛(wèi)星鐘差的精度要求很高，而目前還難以實(shí)時(shí)獲取到精密的衛(wèi)星鐘差信息，因此PPP常用于事后的位移解算.差分定位模式雖然還無需精密衛(wèi)星鐘，但在地震監(jiān)測應(yīng)用中需在附近選擇一個(gè)參考站，而選擇較遠(yuǎn)的測站作為參考站時(shí)，數(shù)據(jù)處理過程變得相對(duì)復(fù)雜.目前廣泛使用的單歷元數(shù)據(jù)處理軟件有GAMIT／GLOBK軟 件 的 TRACK 模 塊［7］，BERNESE 軟 件的運(yùn)動(dòng)學(xué)坐標(biāo)解算模塊［8］，GIPSY軟件的高精度運(yùn)動(dòng)學(xué)數(shù)據(jù)處理模塊［9］，以及Geodetics公司研制的商業(yè)軟件RTD［10］.雖然這些軟件能夠?qū)崿F(xiàn)單歷元高精度的位移解算，但它們難以同時(shí)在實(shí)時(shí)數(shù)據(jù)處理以及長基線高精度解算方面滿足地震地表實(shí)時(shí)監(jiān)測的要求.

本文基于震時(shí)地表震動(dòng)持續(xù)時(shí)間短的特點(diǎn)，以及短時(shí)間內(nèi)雙差殘差的強(qiáng)相關(guān)性，提出了雙差殘差預(yù)報(bào)的GPS位移監(jiān)測方法，對(duì)各雙差殘差進(jìn)行建模，再結(jié)合預(yù)報(bào)殘差實(shí)現(xiàn)短時(shí)間尺度內(nèi)的位移監(jiān)測.最后本文將通過長約1100km的靜態(tài)基線數(shù)據(jù)和El Mayor－Cucapah 7.2級(jí)地震時(shí)94個(gè)測站的數(shù)據(jù)來驗(yàn)證本文方法的有效性.

2 雙差殘差建模

長基線GPS雙差觀測方程可表達(dá)成如下形式：其中Φ為經(jīng)天線相位中心、對(duì)流層模型等改正后的相位觀測值，λ為對(duì)應(yīng)的波長，N為整周模糊度，T為對(duì)流層殘差，I為電離層延遲，O為軌道誤差，M為多路徑延遲，ε為觀測噪聲，i（i＝1，2）代表不同頻率，ΔΔ代表雙差運(yùn)算.與常用的測站天頂延遲估計(jì)加映射函數(shù)的模型不同，對(duì)流層延遲先通過經(jīng)驗(yàn)?zāi)Ｐ透恼?1］，再對(duì)各衛(wèi)星的雙差對(duì)流層殘差分別進(jìn)行建模.

衛(wèi)星軌道誤差所導(dǎo)致的單差測距誤差可表示為［12］：

式（1）中，電離層延遲可通過雙頻消電離層組合消除其影響，而多路徑延遲則可通過改進(jìn)的恒心日濾波的方法來建模［13－14］.經(jīng)過消電離層和恒心日濾波處理后，式（1）可寫成如下形式：

式（3）中的消電離層雙差模糊度不再具有整周特性，但是在沒有周跳發(fā)生的情況下，其是不隨時(shí)間變化的量.一般情況下，天頂對(duì)流層延遲量約為2.4m，可以通過模型有效削除其中90%以上的延遲影響.圖1顯示了不同衛(wèi)星高度角下映射函數(shù)的值.當(dāng)高度角高于25°時(shí)（圖中的虛線處），映射函數(shù)值隨高度角的變化十分緩慢；當(dāng)衛(wèi)星高度角低于20°時(shí)，隨著高度角的降低，映射函數(shù)值急劇增大，因此式（3）中不同衛(wèi)星高度角的雙差對(duì)流層殘差需采用不同的模型.

圖1 不同衛(wèi)星高度角的映射函數(shù)值Fig.1 The values of mapping function of different elevation

經(jīng)過雙差處理后，對(duì)流層殘差將得到進(jìn)一步削弱.對(duì)流層殘差具有極強(qiáng)的時(shí)空相關(guān)性，測站天頂延遲變化常是以小時(shí)為時(shí)間尺度來估計(jì)［15］，同時(shí)考慮到地震發(fā)生時(shí)地表震動(dòng)時(shí)間相對(duì)較短，衛(wèi)星高度角變化也較小，由映射函數(shù)所引起的誤差可近似表達(dá)成線性函數(shù).結(jié)合圖1，雙差殘差一般可簡化成如下模型：

式中t為時(shí)間，θ為衛(wèi)星高度角，θ0為不同模型的分界值，建議值為25°，εr為模型噪聲.由于模糊度是一個(gè)定值，將模糊度合并到式（4）后有

稱a和b為雙差殘差參數(shù)，f（t）為雙差殘差模型.式（3）代入式（5）可得：

地震發(fā)生前，測站精確坐標(biāo)已知，即式（6）中的ΔΔ ρ可視為已知量.通過地震發(fā)生前的若干歷元，結(jié)合最小二乘法即可解算出雙差殘差參數(shù)a、b.

3 基于雙差殘差預(yù)報(bào)的地表運(yùn)動(dòng)監(jiān)測

當(dāng)?shù)卣鸢l(fā)生時(shí)，式（6）中的ΔΔ ρ將會(huì)發(fā)生變化，可以認(rèn)為距震中較遠(yuǎn)（上千公里）的參考站其空間位置保持不動(dòng)，ΔΔ ρ的變化是由于離震中較近的測站震動(dòng)而引起.由于天頂延遲通常是以小時(shí)為間隔來估計(jì)其變化量，而地震地表震動(dòng)時(shí)間短，可忽略其變化量.再結(jié)合第2節(jié)映射函數(shù)相關(guān)影響的分析，可以認(rèn)為地表震動(dòng)過程中雙差殘差不發(fā)生突變，各雙差殘差模型f（t）保持不變.圖2為4.1節(jié)中靜態(tài)試驗(yàn)部分某兩顆衛(wèi)星的雙差殘差序列，可以看出在整個(gè)過程中雙差殘差量變化較為平緩.當(dāng)然在大氣變化異常時(shí)也可能會(huì)出現(xiàn)異常情況，此時(shí)可以在建模和定位過程中進(jìn)行適當(dāng)?shù)呐袛啵綔y是否存在異常.

根據(jù)上面分析，式（6）可表達(dá)為：

f0（t）為震前所確定的各雙差殘差模型f（t）的預(yù)報(bào)值.當(dāng)一個(gè)歷元有n顆衛(wèi)星（參考星除外）可用時(shí)，可得如下觀測方程：

其中上標(biāo)代表不同的衛(wèi)星，利用最小二乘法就能解算出震時(shí)各個(gè)歷元的地表位移.

4 試驗(yàn)結(jié)果分析

為檢驗(yàn)本文方法的有效性，分別采用靜態(tài)數(shù)據(jù)和真實(shí)地震數(shù)據(jù)，利用本文方法進(jìn)行動(dòng)態(tài)定位.試驗(yàn)過程中，采用5min的數(shù)據(jù)來估計(jì)各雙差殘差模型，并用各模型預(yù)報(bào)后5min的雙差殘差，利用預(yù)報(bào)的雙差殘差來解算測站的動(dòng)態(tài)位移，解算時(shí)采用IGS的預(yù)報(bào)精密星歷.

4.1 靜態(tài)試驗(yàn)

靜態(tài)數(shù)據(jù)采用兩個(gè)CORS站的觀測數(shù)據(jù)，基線長約1100km，采樣間隔為1s，時(shí)長為24h，衛(wèi)星截止高度角設(shè)置為5°.以每10min為一個(gè)時(shí)段進(jìn)行分段處理，并剔除衛(wèi)星數(shù)少于5的時(shí)段.

圖3顯示了根據(jù)預(yù)報(bào)的雙差參數(shù)模型解算出的測站動(dòng)態(tài)位移的分布情況，共約解算了38000個(gè)歷元的位移.統(tǒng)計(jì)結(jié)果表明，N、E、U三個(gè)方向的中誤差分別為6mm，6mm，13mm.

圖2 不同衛(wèi)星雙差殘差序列圖Fig.2 Time series of the double differenced residual value for different satellites

進(jìn)一步分析雙差殘差模型，圖4顯示了兩個(gè)具有代表性的雙差參數(shù)的計(jì)算值和模型值，其中，雙差殘差模型是根據(jù)前300s的數(shù)據(jù)解算得出，后300s的模型值為預(yù)報(bào)值.可以看出，當(dāng)衛(wèi)星高度角較低時(shí)，雙差殘差變化顯著，設(shè)置為常數(shù)顯然不合理.另外為了防止高次模型因噪聲影響而出現(xiàn)震蕩現(xiàn)象，應(yīng)盡量采用低次的模型.因此本文根據(jù)不同的衛(wèi)星高度角分別用一次模型和常數(shù)模型來對(duì)雙差殘差建模.

隨著預(yù)報(bào)時(shí)間的增加，雙差殘差的相關(guān)性將會(huì)減弱.對(duì)于一般地震來說，地面震動(dòng)時(shí)候一般不超過5min，因此本文將只分析在殘差預(yù)報(bào)的時(shí)間達(dá)到5min時(shí)，本文方法進(jìn)行位移監(jiān)測所能達(dá)到的精度.為了削弱隨機(jī)噪聲的影響，將根據(jù)最后20個(gè)歷元的位移信息取平均來計(jì)算最終的位移量.本試驗(yàn)共解算了126個(gè)時(shí)段的數(shù)據(jù)，結(jié)果如圖5所示.

可以看出，水平方向的最大值不超過20mm，豎直方向基本不超過40mm.統(tǒng)計(jì)結(jié)果表明，N、E方向超過92.0%的位移小于10mm，U方向超過91.2%的位移小于20mm.

4.2 El Mayor－Cucapah地震位移解算

2010年4月4日，在墨西哥Baja地區(qū)發(fā)生了7.2級(jí)強(qiáng)烈地震，主震持續(xù)時(shí)間超過了40s，地表破裂長度達(dá)120km.本文利用加州實(shí)時(shí)觀測網(wǎng)絡(luò)（CRTN）的1Hz GPS數(shù)據(jù)［16］，通過本文的方法解算該地震各測站的動(dòng)態(tài)位移和地震同震位移場，參數(shù)設(shè)置與靜態(tài)試驗(yàn)相同.

本次試驗(yàn)共解算了94個(gè)測站的位移，p090離震中約920km，將其作為參考站.主震之后的部分位移場圖如圖6所示.

圖7畫出了其中離震中較近和較遠(yuǎn)的幾個(gè)測站水平方向的實(shí)時(shí)動(dòng)態(tài)位移信息，可以看出，利用本文的方法解算出的GPS時(shí)間序列可清楚記錄各個(gè)測站的地表震動(dòng).

為驗(yàn)證解算結(jié)果，將本方法解算的同震位移與GPS Explorer提供的高精度后處理結(jié)果［16］進(jìn)行對(duì)比，統(tǒng)計(jì)了其中66個(gè)測站的同震位移的差值，統(tǒng)計(jì)結(jié)果如圖8.可以看出水平方向偏差基本小于20mm，豎直方向偏差基本小于30mm.統(tǒng)計(jì)結(jié)果表明，與后處理結(jié)果相比，N、E方向分別有79%和91%的測站相差小于10mm，U方向有82%的測站相差小于20mm.

進(jìn)一步分析N方向差值最大的4個(gè)測站，偏差從大到小依次為p494、p744、p497和p498.從圖6中可以看出這4個(gè)點(diǎn)距震中較近，且位置相對(duì)集中.由于事后處理的結(jié)果是根據(jù)震后幾天的數(shù)據(jù)解算得出，其位移場是震后幾天的位移場，而本文解算的是震后幾分鐘的位移場.兩種解算結(jié)果在4個(gè)測站N、E方向的位移差分別為：－34、－28、－25、－19mm和9、－3、3、0mm.可以看出它們主要呈現(xiàn)出往南偏的趨勢，這與同震位移的方向大體保持一致，初步分析兩種方法解算結(jié)果中較大的差異主要是由于震后滑移所引起.

SOPAC提供的p494和p497兩個(gè)測站（p744、p498缺少）的快速靜態(tài)解和本文的解算結(jié)果吻合較好，兩測站N方向相差約為8mm和14mm，E方向相差約2mm和16mm，進(jìn)一步證明了該方法的有效性.

5 結(jié) 論

根據(jù)地震地表震動(dòng)時(shí)間短的特點(diǎn)、以及短時(shí)間尺度內(nèi)雙差殘差具有極強(qiáng)的相關(guān)性，提出了基于雙差殘差預(yù)報(bào)的強(qiáng)震地表運(yùn)動(dòng)實(shí)時(shí)監(jiān)測方法.與傳統(tǒng)的基線解算方法不同，該方法不單獨(dú)解算模糊度信息，而是直接對(duì)各雙差殘差建模，并利用殘差模型預(yù)報(bào)實(shí)現(xiàn)上千公里基線的解算.試驗(yàn)結(jié)果表明，在5min的預(yù)報(bào)時(shí)間內(nèi)，1Hz的靜態(tài)數(shù)據(jù)該方法N、E、U三個(gè)方向的中誤差分別為6mm、6mm和13mm，El Mayor－Cucapah地震數(shù)據(jù)的解算結(jié)果也與實(shí)際情況有較好的一致性，試驗(yàn)結(jié)果證明了本文方法的可行性.本方法只需要預(yù)報(bào)精密星歷，且計(jì)算效率高，在有實(shí)時(shí)觀測數(shù)據(jù)的情況下，可實(shí)現(xiàn)強(qiáng)震時(shí)測站位移的實(shí)時(shí)監(jiān)測，這對(duì)地震預(yù)警和地震快速救援具有一定的實(shí)用和研究價(jià)值.

需要注意的是雙差參數(shù)預(yù)報(bào)的精度會(huì)隨著預(yù)報(bào)時(shí)間的增加而降低，本文方法不適用于長時(shí)間的地表運(yùn)動(dòng)緩慢的監(jiān)測.本文初步測試了5min預(yù)報(bào)時(shí)間內(nèi)方法的有效性，對(duì)于更長時(shí)間的位移監(jiān)測，方法的有效性尚需進(jìn)一步研究.此外，周跳對(duì)本方法的影響不容忽視，其將影響殘差模型和后面的位移解算.

（References）

［1］Shi C，Lou Y D，Zhang H P，et al.Seismic deformation of the Mw8.0Wenchuan earthquake from high－rate GPS observations.Advances in Space Research，2010，46（2）：228－235.

［2］Elósegui P，Davis J L，Oberlander D，et al.Accuracy of high－rate GPS for seismology.Geophys.Res.Lett.，2006，33（11）：L11308，doi：10.1029／2006G1＿026065.

［3］Langbein J，Bock Y.High－rate real－time GPS network at Parkfield： Utility for detecting fault slip and seismic displacements.Geophys.Res.Lett.，2004，31（15）：L15S20，doi：10.1029／2003GL019408.

［4］殷海濤，張培震，甘衛(wèi)軍等.高頻GPS測定的汶川Ms8.0級(jí)地震震時(shí)近場地表變形過程.科學(xué)通報(bào)，2010，55（26）：2621－2626.Yin H T，Zhang P Z，Gan W J，et al.Near－field surface movement during the Wenchuan Ms8.0earthquake measured by high－rate GPS.Chinese Science Bulletin，2010，55（23）：2529－2634.

［5］Crowell B W，Bock Y，Squibb M B.Demonstration of earthquake early warning using total displacement waveforms from real－time GPS networks.Seismological Research Letters，2009，80（5）：772－782.

［6］方榮新，施闖，徐培亮等.GPS地震儀：PANDA軟件測試結(jié)果與驗(yàn)證.武漢大學(xué)學(xué)報(bào) （信息科學(xué)版），2011，36（4）：453－456.Fang R X，Shi C，Xu P L，et al.GPS seismometer：PANDA software testing results and validation.Geomatics and Information Science of Wuhan University （in Chinese），2011，36（4）：453－456.

［7］http：／／geoweb.mit.edu／～tah／track＿example／

［8］Bock Y，Nikolaidis R，de Jonge P J，et al.Instantaneous geodetic positioning at medium distances with the global positioning system.J.Geophys.Res.，2000，105（B12）：28233－28253.

［9］Jet Propulsion Laboratory California Institute of Technology.GIPSY5.0Release Note.https：／／gipsy－oasis.jpl.nasa.gov／gipsy／docs／Release＿Nates＿5.0.pdf，2008

［10］http：／／www.geodetics.com／products／software／rtd＿pro.php

［11］Saastamoinen J.Atmospheric correction for the troposphere and stratosphere in radio ranging of satellites，in the use of artificial satellites for geodesy.Geophysics Monograph，1972，15（16）：247－251.

［12］Han S W.Carrier phase－based long－range GPS kinematic positioning［Ph.D.thesis］.Sydney：The University of New South Wales，1997.

［13］Choi K，Bilich A，Larson K M，et al.Modified sidereal filtering： Implication for high－rate GPS positioning.Geophys.Res.Lett.，2004，31：L24610，doi：10.1029／2004GL021621.

［14］Niell A E.Global mapping functions for the atmosphere delay at radio wavelengths.J.Geophys.Res.，1996，101（B2）：3227－3246.

［15］熊永良，黃丁發(fā)，丁曉利等.虛擬參考站技術(shù)中對(duì)流層誤差建模方法研究.測繪學(xué)報(bào)，2006，35（2）：118－121.Xiong Y L，Huang D F，Ding X L，et al.Research on the modeling of tropospheric delay in virtual reference station.Acta Geodaetica et Cartographica Sinica （in Chinese），2006，35（2）：118－121.

［16］http：／／geoapp03.ucsd.edu／gridsphere／gridsphere？cid＝El＋Mayor＋Cucapah