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

西南交通大學測量工程系，成都 610031

1 引 言

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

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

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

2 雙差殘差建模

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

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

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

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

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

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

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

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

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

3 基于雙差殘差預報的地表運動監(jiān)測

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

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

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

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

4 試驗結果分析

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

4.1 靜態(tài)試驗

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

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

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

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

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

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

4.2 El Mayor－Cucapah地震位移解算

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

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

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

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

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

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

5 結 論

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

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

（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級地震震時近場地表變形過程.科學通報，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軟件測試結果與驗證.武漢大學學報 （信息科學版），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ā)，丁曉利等.虛擬參考站技術中對流層誤差建模方法研究.測繪學報，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