實時高頻GPS在地震學中的應用研究*

殷海濤，劉希強，甘衛軍

0 引言

GPS技術作為地學研究的一個有效途徑，具有全天候、高精度、自動化、高效益等顯著特點，在地殼運動監測和地球動力學等領域發揮著越來越重要的作用 (Wang et al，2001;Gan et al，2007;Zhang，Gan，2008)。目前GPS通常利用日均值來分析和研究長周期的地殼形變信息，可為中長期地震預測提供可靠的依據，但它不能捕捉內地殼瞬時形變信息。

隨著接收機技術和存儲能力的提高，我們可以獲取高頻 (1 Hz)和超高頻 (20～50 Hz)GPS觀測數據，而單歷元GPS處理技術的逐漸成熟，也使高頻GPS有能力觀測到周期大于1 s的真實地表位移。目前國內外學者已經利用高頻GPS觀測網成功的獲取了幾次大地震所引起的瞬時地表運動狀態，如2003年加州San Simeon 6.5級地震 (Ji et al，2004)，2004年 Sumatra-Andaman 9.3級巨震 (Ohta et al，2006)，2008年汶川8.0級地震(Yin et al，2010)和 2011年日本 9.0級地震(Grapenthin，Freymueller，2011)等，這為我們更好地理解地震破裂過程以及地震波在地表的傳播方式開辟了新的途徑。

通訊技術的飛速發展使高頻GPS數據的實時傳輸成為了可能，可獲取實時地表位移狀態，彌補現有技術手段的不足。當前實時高頻GPS技術成為國際研究的新熱點 (Allen，Ziv，2011)。本文介紹了實時高頻GPS技術及其數據處理方法，并與傳統地震 (強震)儀進行了對比分析，論述其如何應用于地震預警等實時地震學領域。

1 實時高頻GPS及處理方法

GPS全面建成至今只有不到20年的時間，但卻廣泛地應用于導航、大地測量、攝影測量、地震學等多個領域。隨著技術和方法的不斷發展，其定位精度和采樣頻率有了明顯提高，目前GPS定位精度可達亞毫米級，相對定位精度可達10－9。而數據采樣率從30 s提升到了0.02～1 s，甚至0.01 s，極大拓展了觀測的時間尺度。實時高頻GPS技術 (簡稱RH GPS)是目前GPS發展的最新階段，主要是通過接收機和4顆以上的衛星來實時獲取高頻數據 (采樣率高于1 s)，并通過通訊網絡和處理中心實時解算出地面點亞厘米級的位移。其數據處理方法主要有差分模式和非差模式兩種，差分模式必須要選擇一個參考站點，然后對其他站點進行相對定位，非差模式主要依賴于高精度的衛星軌道信息。兩種方法在后處理情況下精度相當，但在高精度的實時應用中，非差模式由于無法及時獲取可靠的軌道信息而受限。

目前國際上的GAMIT/GLOBK、RTD和GIPSY軟件均可做實時高頻GPS數據處理，前兩種軟件為差分模式定位，GIPSY為非差模式定位。本文主要詳細闡述GAMIT/GLOBK的實時處理方法，該軟件中的TrackRT模塊是由麻省理工學院的Thom-as Herring教授于2010年開發出來的，主要用來進行實時GPS數據處理。該模塊可直接得出站點三維位移時間序列和對流層延遲，其水平位移精度優于1 cm，垂直方向精度優于2 cm。目前該軟件已應用于UNAVCO，USGS，PBO等GPS觀測網的實時解算和發布。雖然實時高頻GPS數據的解算方法與高頻GPS類似，但在其軟件安裝、數據獲取、衛星軌道的選擇等方面有很大差異。

1.1 實時數據獲取

在安裝、運行TrackRT模塊前，必須要安裝BKG NTRIP Client(BNC)和QT軟件的libraries和include文件。BNC軟件主要是利用互聯網協議和客戶端來傳輸RTCM數據，可通過互聯網端口為trackRT模塊提供GPS原始數據。BNC軟件除了可以實時傳輸GPS原始數據外，還具有將原始數據轉換為RINEX格式數據、存儲星歷文件、發布數據、精密單點定位等功能，可通過 http://igs.bkg.bund.de/ntrip/download下載。

1.2 衛星星歷文件

實時GPS數據處理時所需的衛星軌道文件也必須滿足實時的要求。目前國際GPS服務機構(IGS)提供了3種衛星星歷:最終星歷 (IGF)、快速星歷 (IGR)和預報星歷 (IGU)。對于IGS最終星歷來說，精度雖然很高，但時間延遲12～18 d，這限制了實時或準實時用戶的使用。IGS提供的快速精密星歷 (時間延遲17～41 h)，也不能滿足實時的要求。為了實現實時應用，IGS機構于2000年3月5日 (GPS周1 052周)提供了超快星歷。超快星歷和超快衛星鐘差的更新率為6 h，在每天 UTC(CoordinateUniversalTime)3:00、9:00、15:00和21:00各發布一次。超快星歷的歷元間隔為15 min，共包括48 h的軌道信息，前24 h是基于跟蹤站的觀測值計算得到，后24 h是外推預報得到的。實時GPS數據處理只能利用后24 h外推結果來進行計算。除了IGS發布的精密星歷外，觀測數據本身也可實時獲取導航星歷，但其精度較低，不建議使用。表1對現有GPS衛星星歷文件進行了比較，發現只有預報星歷中外推的24 h結果在精度和時間延遲上滿足實時高頻GPS數據處理的要求。

1.3 實時處理方法

實時處理GPS數據時，TrackRT模塊中的部分命令與Track模塊相同，需要在命令文件中設置一個參考站，然后進行差分處理，但由于實時數據流中并不包括站點的位置坐標，所以必須要預先給定站點坐標值。在命令文件中還需利用SITE_STATS和ATM_STATS命令對站點的坐標噪聲和大氣延遲進行先驗約束，確保其解算精度。參考站的選擇方法可參考殷海濤等 (2012)，為了使參考站“固定不動”，通常將參考站的這兩組值設置為0。在命令文件中還需設定站點的天線和接收機型號，以及更新數據碼偏心文件 (DCB)，這樣不僅可以提高解算精度，還可減少計算量，提高計算效率。

表1 GPS星歷文件對比Tab.1 Comparison of GPS satellite ephemeris files

實時GPS數據處理的核心問題也是模糊度的快速解算，TrackRT模塊使用Melbourne-Wubbena寬相 (MW－WL)，甚寬相 (EX－WL)和消電離層浮點估計 (LC)來計算整周模糊度。設L1和L2的模糊度整周數為N1，N2，則:

(1)MW－WL是基于相位和距離數據來估計N1－N2;

(2)EX－WL=N1－ (f1/f2)*N2是L1的整周數，但L2的周數變為了1.283*N2;

(3)LC=2.546*N1－1.984*N2。

EX－WL不受幾何距離變化的影響，但是依賴于電離層延遲。對于短基線來說，EX－WL應該接近于0才符合N1和N2的選擇。當N1和N2正確的時，LC殘差也應該接近于0。

在定位精度方面，目前實時高頻GPS的定位精度在亞厘米級，這表示它可以為大地震的震級估計提供約束。在解算時間方面，實時高頻GPS的數據處理在理論上可達到實時解算，但在實際操作過程中，其數據采集、傳輸和處理受到通訊設備、處理機配置等硬件因素的影響，時間延遲約為2～5 s，可為地震預警等實時地震學應用提供幫助。

2 與地震 (強震)儀的比較

震時地表位移是估計地震破裂過程和地震震級的一項重要參量，傳統的方法是對地震儀或強震儀所測速度 (加速度)進行積分得到 (殷海濤等，2009)。雖然地震儀 (強震儀)和高頻GPS都可以獲得震時的地表位移信息，但兩種技術在觀測方法上有著本質的區別:慣性地震儀的初始狀態是平靜的 (不工作的)，當地面震動時則產生測量值，震動結束后恢復平靜。數字地震儀提供的是一系列速度或加速度觀測值，其值是在一個慣性框架下產生的，然后再經過積分 (二次積分)得到站點位移。相對而言，GPS記錄的是由GPS衛星來確定的GPS天線的位置，其位置時間序列是在地球參考框架下計算出來的。表2列舉了兩種技術在性能上的差異 (殷海濤等，2012)。

表2 高頻GPS與地震儀 (強震儀)性能對比Tab.2 Performance comparison of high rate GPS and seismometer(strong motion seismography)

通過分析對比高頻GPS與地震 (強震)儀的工作原理及性能，可以發現高頻GPS技術可以為傳統地震學技術進行有效補充:

(1)位移時間序列是GPS的直接觀測量，但是對于地震 (強震)儀來說，其觀測量為速度(加速度)，所以必須要對其進行積分處理才能得到位移，而積分處理經常伴隨著傳感器旋轉、傾斜、滯后造成的漂移和積分過程中出現的不確定性，很容易產生放大噪聲和扭曲真實信號的情況(Boore et al，2002);

(2)當大地震發生時，地震儀能夠產生振幅飽和現象，為了不記錄到滿幅的速度和加速度，采取限幅的方法，而GPS在振幅方面不會產生飽和，它沒有儀器響應來限制接收機的觀測能力，相對而言，位移越大其定位相對精度會越高，所以GPS技術更適用于大震造成的地表形變;

(3)地震儀只有在地面震動后才能產生觀測值，但很多地震事件如:斷層滑移，巖石粘性變化，流性變化和震后形變等，雖然產生了地殼形變，但是并未產生地震波。在火山活動和非構造性變形過程等，由于形變的時間尺度較大而未必產生地震波。利用高頻GPS監測這些事件，可以較好地獲得低頻形變信息，進而了解慢形變過程與機理;

(4)目前的地震 (強震)儀在大震中很難提取出長周期位移。而Allen和Ziv(2011)，Bock等(2011)認為:實時GPS定位得到的最穩定的信息就是水平位移，可有效彌補傳統地震技術的不足，為大地震的震級估計和地震破裂過程研究提供非常重要的信息。

雖然高頻GPS在地震學應用上具有很多的優勢，在大地震發生時可以記錄到明顯的位移等，但它也存在一些限制因素:(1)相對于地震儀而言，高頻GPS的噪聲基底較大，Ji等 (2004)分析認為，1～Hz GPS數據計算的EW向，NS向和垂直方向位置標準差可達到3 mm，7 mm和11 mm。Langbein和 Bock(2004)，Bock等 (2011)研究認為1 Hz GPS觀測數據可以在99%的置信水平，估計振幅超過6 mm的水平位移，因此目前高頻GPS還不足以記錄小震所造成的微弱地表形變。(2)目前GPS接收機僅能獲得1 Hz或者超高頻(20～50 Hz)的數據，但是對于更高頻 (約200 Hz)的信號，它最終還是要受限于接收機的頻帶寬度。

因此，結合GPS和地震儀資料各自發揮其優勢，可起到相互補充的作用，目前GPS觀測技術與地震學的觀測譜范圍逐漸合并，逐步發展成了一門新興的學科——GPS地震學 (Nikolaidis et al，2001)，將會把傳統地震學向前推進一大步。

3 實時地震學應用

目前地震預警系統 (Earthquake Early Warning，簡稱EEW)等實時地震學應用大都只基于地震儀觀測網，利用P波檢測來估計震級，這對快速估計中小地震震級是非常有效的。但目前有可靠的證據證明 (Yamada et al，2008)，在利用前幾秒的P波信息估計震級時，當M＞7地震發生時，震級估計會出現飽和。也就是傳統的P波檢測方法無法區分所探測到的地震是6.8級還是7.8級，雖然說這并不影響地震預警的發布，但是以2011年3月11日日本的9.0級巨大地震為例，就是因為低估了震級，才導致了對地震和海嘯沒有做出相應的防御措施，因此無論從科學研究還是從實際應用的角度來講，都需要繼續開展這方面的研究。

3.1 高頻GPS與強震儀的結合

地震震級可以通過震時位移的振幅反演來進行估計，在獲取震時地表位移的方法中，只有有效結合地震儀和高頻GPS才能得到可靠的高精度強震動地表位移。Bock等 (2011)引入了連續Kalman濾波 (卡爾曼濾波)和最優平滑濾波等方法來合并強震儀與高頻GPS結果。

卡爾曼濾波器包括兩個主要過程:預估與校正。預估過程主要利用時間更新方程建立對當前狀態的先驗估計，及時向前推算當前狀態變量和誤差協方差估計的值，以便為下一個時間狀態構造先驗估計值;校正過程負責反饋，利用測量更新方程在預估過程的先驗估計值及當前測量變量的基礎上建立起對當前狀態改進的后驗估計。

我們假定每個站點的坐標方向 (N，E，U)都是獨立的一維運動，將每個方向的運動作為一階線性差分等式，利用前向卡爾曼濾波來融合GPS和強震儀數據并監測狀態變量。在每個時間步長中，每個坐標方向的位移和速度可表示為

時間更新 (先驗估計)為

測量值更新 (后驗估計)為

其中，X為狀態變量，k為時間歷元，P－為先驗估計誤差協方差矩陣，P為后驗估計誤差協方差矩陣;A、B和H為狀態變換矩陣，是狀態變換過程中的調整系統，數值是從建立的系統數學模型中導出來的;Q為過程激勵噪聲協方差矩陣，R為觀測噪聲協方差矩陣;觀測變量Z表示GPS位移，U為系統控制輸入，表示強震儀數據。

強震儀時間序列提供了系統輸入，其采樣率通常為80～250 Hz;而GPS位移控制測量過程，其采樣率通常為1～5 Hz。這樣式(3)在每一個時間步長內進行時間和測量值更新，如果強震儀采樣率是GPS采樣率的整數倍，則可以簡化公式。如果采樣頻率自始至終保持不變，其噪聲特性也不改變，則式(3)中的各項參數A，B，Q和R均保持不變。此方法的特點是矩陣的維數較少，可以使數值計算更加簡潔。此外，卡爾曼濾波只需要當前的數據，這樣可以實現實時應用。

應用前向卡爾曼濾波可以解決預測問題，但如果已經存在了一部分數據，在計算中會用到所有存在的和未來的數據，在任意點的估計值可以進行優化。雖然平滑濾波通常應用于非實時處理，但可以通過約束平滑的時間間隔來實現近實時的處理。固定時間間隔的平滑濾波包括3步:(1)前向卡爾曼濾波;(2)后向卡爾曼濾波，將時間順序轉換到整個時間序列中。前兩步提供的分散信息流在第三步中結合，并計算出最終狀態值。結果如圖1所示，通過濾波處理，兩種技術手段可以進行有效結合，獲取更精確的寬相位移，而且也減少了處理過程中的主觀因素。

3.2 地震矩震級估計

快速的矩張量 (Moment Tensor)估計為地震響應、海嘯預警等提供了有價值的信息。目前所有的方法基本都使用地震儀資料，利用時間域的波形匹配反演方法來實時估計 (Dreger，2003)，這對中小地震來說非常有效，但對較大地震來說，仍然是研究的熱點。

Kanamori和Rivera(2008)發展了W相位反演方法，這使大地震CMT(Centroid moment tensor)的計算前進了一大步，目前該方法正應用于USGS等機構。雖然W相位反演方法非常穩定，但也需要長周期位移資料，而目前的地震 (強震)儀很難在大震中提取出長周期位移，限制了其進一步發展。Diego等 (2012)提出了一種利用實時高頻GPS反演地震矩張量的方法fastCMT，該反演方法主要利用地表的同震位移來反演深部的震源參數，主要由以下4部分組成:

(1)在指定區域構建反演節點，其節點格網的構建方法有兩種:一是需要預先設定斷層面模型，這就要求對該地區的地質條件有很深入的了解;二是不需要已知斷層模型，直接根據格網搜索法來進行估計，這對斷層復雜的地區非常有幫助，但是必須要有足夠密度的GPS觀測網才能實現;

圖1 卡爾曼濾波前后位移時間序列①汶川地震對XYAN站點觀測數據.(a)南北向分量;(b)東西向分量;(c)高程分量Fig.1 Displacement time series before and after Kalman filter(a)NS-component;(b)EW-component;(c)evaluation component

(2)根據地殼速度模型計算格林函數，所用程序為 EDGRN/EDCMP。從Hankel變換得到的部分差分運動等式的結果開始，然后應用Thomson－Haskell傳播矩陣將地面和深部的變形進行關聯，其簡化公式為

結合公式 (4)和 (5)，則得到

(3)高頻GPS數據的獲取及處理，首先利用移動平均法對獲取的高頻GPS時間序列進行濾波處理，得到永久位移，然后再對其結果進行加權處理。在反演過程中，就使用了兩個權重，一是觀測噪聲權重，用來提高觀測值的可靠性。二是震中距權重，為了避免近場最大同震形變主導整個反演過程。其觀測噪聲加權公式為

位移分量震前60 s的觀測噪聲標準差。由于實時GPS位移噪聲在幾分鐘尺度內基本保持穩定，而且Bock等 (2011)利用震動臺試驗表明震動前和震動過程中的噪聲水平并沒有明顯的改變，所以我們可以設定觀測噪聲和權重矩陣在整個反演過程中保持不變。

為了避免地面運動衰減產生的影響，我們利用站點震中距作為另一個權重。Dreger(2003)給出

(4)地震矩震級估計，在實時應用過程中，每獲取一個歷元的觀測值就反演一次，產生新的MT值

然后，可以利用 Hanks和 Kanamori(1979)提出來的關系式計算矩震級

其中，M0單位為dyn·cm。

通過上述計算可以實時估計出地震矩震級，并可通過格網搜索法估計地震矩心的位置，利用方差減少法計算反演差值，將方差減少值VR最大的節點視為矩心，公式為

式中，d為GPS站點到節點的距離。

fastCMT方法的特點是可以在沒有任何斷層模型的情況下實時估計大震矩震級和矩心位置，但是前提條件是近震區裝備有足夠密度的高頻GPS站點。Diego等 (2012)通過2003年日本Tokachioki MW8.3地震和2010年EI Mayor-Cucapah MW7.2地震的實際數據分析認為:在地震發生2～3 min，該方法便可計算出精確的地震矩震級，可為地震快速響應和海嘯預警等提供有力支持。

4 討論與結論

實時高頻GPS是目前GPS技術發展的最新階段，大大拓展了其觀測信息頻域和應用范圍，也將使現有的GPS觀測網發揮更大的作用。但由于其研究剛剛起步，在數據處理和地震應用過程中仍存在很多問題亟待解決:

(1)模糊度的快速解算仍然是實時GPS數據處理的核心問題，也是一直以來國際研究的熱點問題。除此以外，接收機進行觀測時易受外部環境因素干擾，產生突跳或衛星失鎖的情況，如何消除突跳和快速收斂也是需要解決的技術難點;

(2)對于差分模式定位來講，如何選擇參考站是主要問題，而非差模式的問題在于如何獲取高精度的衛星軌道信息。對于地震預警系統來說，兩者各有利弊，差分模式可應用于區域預警，而非差模式更適用于現地預警;

(3)如何消除實時GPS觀測中的噪聲影響，得到更可靠的位移信息，是目前所要面臨的主要問題;

(4)如何更好的結合地 (強)震儀與高頻GPS資料，獲取穩定性高敏感性高的震時地表位移是下一步應用研究的重點;

(5)在獲取地表位移后，如何快速提取地表永久位移是估計地震震級的關鍵因素，也是以后研究的主要問題。

在地震學應用中，與傳統地震 (強震)儀相比，RH GPS可直接測量震時地表位移，增強震級估計 (特別是大地震)的可靠性，避免了傳感器出現的旋轉、傾斜、滯后造成的漂移和積分過程中出現的不確定因素。而且它可以不受振幅限制，有效的彌補了地震 (強震)儀的不足，但由于其定位精度要低于地震 (強震)儀，對微小地震不敏感，所以結合兩種手段才能達到最佳的監測效果。

利用實時高頻GPS技術能快速獲取震時地表位移，可以計算出精確的地震矩震級可以為地震破裂過程研究、地震預警、海嘯預警等地震研究和應用提供一項新的手段，因此在我國及時開展此領域的研究工作是十分必要和迫切的。

殷海濤，甘衛軍，肖根如，等.2009.利用高頻GPS技術進行強震地面運動監測的研究進展［J］．地球物理學進展，24(6):2 012－2 019.doi:10.3969/j.issn.1004－2903.2009.06.011.

殷海濤，肖根如，張磊，等.2012.TRACK高頻GPS定位中震時參考站的選取方法［J］．大地測量與地球動力學，32(4):15－19.

Allen R M，Ziv A.2011.Application of real-time GPS to earthquake early warning［J］．GeophysResLett，38(16):1 － 6，doi:10.1029/2011GL047947.

Bock Y，Diego M，Crowell B W.2011.Real-time strong-motion broadband displacements from collocated GPS and Accelerometers［J］．BSSA，101(6):2 904 －2 925，doi:10.1785/0120110007.

Boore D M，Stephens C D，Joyner W B，et al.2002.Comments on baseline correction of digital strong-motion data:examples from the 1999 Hector Mine，California earthquake［J］．BSSA，92(4):1543 － 1560，doi:10.17850120000926.

Diego M，Bock Y，Crowell B W.2012.Real-time centroid moment tensor determination for large earthquakes from local and regional displacement records［J］．Geophys J Int，188(2):703 － 718，doi:10.1111/j1365－246X.2011.05297.x.

Dreger S D.2003.TDMT_INV:Time domain seismic moment tensor inversion［M］//William H K L，Hiroo K，Paul J，et al.International handbook of earthquake and engineering seismology.American:Academic Press，81B:1 627 －1 631.

Gan W J，Zhang P Z，Zheng K S，et al.2007.Present-day crustal motion within the Tibetan Plateau inferred from GPS measurements［J］．JGR，112(B8):1 －14，doi:10.1029/2005JB004120.

Grapenthin R，Freymueller J T.2011.The dynamics of a seismic wave field:Animation and analysis of kinematic GPS data recorded during the 2011 Tohoku-oki earthquake，Japan［J］．Geophys Res Lett，38(18):1 －5，doi:10.1029/2011GL048405.

Hanks T C，Kanamori H.1979.A moment magnitude scale［J］．JGR，84(B5):2 348－2 350.

Ji C，Larson M，Tan Y，et al.2004.Slip history of the 2003 San Simon earthquake constrained by combining 1 － Hz GPS，strong motion，and teleseismic data［J］．Geophys Res Lett，31(17):1 － 4，doi:10.1029/2004GL020448.

Kanamori H，Rivera L.2008.Source inversion of W phase:speeding up seismic tsunami warning［J］．Geophys J Int，175(1):222 － 238.

Langbein J，Bock Y.2004.High-rate real-time GPS network at Parkfield:Utility for detecting fault slip and seismic displacements.Geophys Res Lett，31(15):1 － 4，doi:10.1029/2003GL019408.

Nikolaidis R，Bock Y，Jonge P J，et al.2001.Seismic wave observation with the Global Positioning System［J］．JGR，106(B10):218 097 －21 916.

Ohta Y，Meiano I，Sagiya T，et al.2006.Large surface wave of the 2004 Sumatra-Andaman earthquake captured by the very long baseline kinematic analysis of 1 － Hz GPS data［J］．Earth Planets Space，58(2):153－157.

Wang Q，Zhang P Z，Jeffrey T F，et al.2001.Present-day crustal deformation in continental China constrained by Global Positioning System measurements［J］．Science，294(5542):574 － 577.

Yamada M，Heaton T，Beck J.2008.Real-time estimation of fault rupture extent using envelopes of acceleration［J］．BSSA，98(2):607－619.

Yin H T，Zhang P Z，Gan W J，et al.2010.Near-field surface movement during the Wenchuan MS8.0 earthquake measured by high-rate GPS［J］．Chinese Sci Bull，55(23):2529 － 2534，doi:10.1007/s11434－010－4026－2.

Zhang P Z，Gan W J.2008.Combined model of rigid-block motion with continuous deformation:Patterns of present-day deformation in Continental China［J］．Geological Society of America Special Paper，444:59 －71，doi:10.1130/2008.2444(04).