聶兆生 祝芙英 付寧波

1)中國(guó)地震局地震研究所，武漢4300712)地殼運(yùn)動(dòng)與地球觀測(cè)實(shí)驗(yàn)室，武漢4300713)江西省測(cè)繪產(chǎn)品質(zhì)量監(jiān)督檢驗(yàn)站，南昌330046

Kalman濾波在地震電離層TEC異常探測(cè)中的應(yīng)用*

聶兆生1，2)祝芙英1，2)付寧波3)

1)中國(guó)地震局地震研究所，武漢430071

2)地殼運(yùn)動(dòng)與地球觀測(cè)實(shí)驗(yàn)室，武漢430071
3)江西省測(cè)繪產(chǎn)品質(zhì)量監(jiān)督檢驗(yàn)站，南昌330046

基于中國(guó)地殼運(yùn)動(dòng)觀測(cè)網(wǎng)絡(luò)GPS觀測(cè)資料解算的電離層TEC數(shù)據(jù)，利用Kalman濾波方法對(duì)2008年5月12日汶川Ms8.0地震前的電離層TEC進(jìn)行異常探測(cè)研究，并與以前的處理結(jié)果進(jìn)行對(duì)比，對(duì)比結(jié)果表明:利用Kalman濾波方法能夠有效地探測(cè)到震前的電離層TEC異常擾動(dòng)。

Kalman濾波;電離層;總電子含量;異常;探測(cè)

AbstractIt has been verified by plentiful existing observations that the ionospheric disturbance do appear prior to earthquake.In this paper，on the basis of the ionospheric TEC derived from the GPS observation data from the reference stations of Crustal Movement Observational Network of China，we processed and analyzed the ionospheric TEC data prior to the WENCHUAN Ms8.0 earthquake by Kalman filtering and then we compared the results with the previous conclusion，the results indicate that the establishment of the model of the Kalman filter is reasonable and reliable in the detection of the ionospheric TEC anomaly prior to earthquake effectively.

Key words:Kalman filtering;ionosphere;TEC(Total Electron Content);anomaly;detecting

1 引言

研究表明，由地震引起的電離層擾動(dòng)的確存在，而且在震級(jí)大于5級(jí)的地震發(fā)生前的幾天到幾個(gè)小時(shí)內(nèi)會(huì)發(fā)生電離層擾動(dòng)［1－8］。目前，電離層TEC異常探測(cè)方法主要有滑動(dòng)均值法、中位數(shù)法、四分位距法［9］等，雖然對(duì)一些顯著的TEC異常擾動(dòng)也能探測(cè)出來，但由于異常探測(cè)的背景參考值大多為當(dāng)月所有的觀測(cè)資料，即對(duì)照及探測(cè)的觀測(cè)資料在時(shí)間上已經(jīng)涵蓋了震發(fā)時(shí)刻，從而使得探測(cè)結(jié)果的可靠性有所降低。

Kalman濾波作為一種新的重要的最優(yōu)估計(jì)理論已被廣泛地應(yīng)用于各種測(cè)量數(shù)據(jù)的處理應(yīng)用中［10－13］。鑒于Kalman濾波的特點(diǎn)，基于中國(guó)地殼運(yùn)動(dòng)觀測(cè)網(wǎng)絡(luò)的GPS觀測(cè)數(shù)據(jù)解算的電離層總電子密度(TEC，Total Electron Content)資料，本文嘗試用該方法對(duì)汶川Ms8.0地震前的電離層TEC資料進(jìn)行處理和分析，以期能夠探測(cè)到震前的電離層TEC異常擾動(dòng)。

2 電離層TEC的解算

GPS反演的方法具有臺(tái)站分布廣、數(shù)量多，利用一臺(tái)GPS接收機(jī)可以同時(shí)獲得至少4個(gè)電離層星下點(diǎn)處的垂直TEC等優(yōu)點(diǎn)，可以獲得大尺度范圍內(nèi)的電離層狀態(tài)信息。在GPS反演電離層時(shí)，通常使用的是單層球殼模式，即把電離層近似成某一個(gè)高度的無限薄的球殼，TEC就是球殼上的電子濃度，數(shù)據(jù)處理中，我們選取的球殼高度為350km，觀測(cè)資料為偽距觀測(cè)量與載波相位觀測(cè)量，取球諧函數(shù)作為電離層TEC的擬合模型，考慮了硬件延遲誤差的影響［14，15］?；谥袊?guó)地殼運(yùn)動(dòng)監(jiān)測(cè)網(wǎng)絡(luò)的近30個(gè)GPS跟蹤站的觀測(cè)數(shù)據(jù)，我們獲得了中國(guó)區(qū)域上空的電離層地圖(GIM，Global Ionospheric Map)。為檢驗(yàn)TEC解算結(jié)果的可靠性，在GIM中任取一點(diǎn)(35°N，105°E)，以歐洲定軌中心(CODE)發(fā)布的最終電離層TEC產(chǎn)品為真值，將本文解算的TEC值與真值進(jìn)行比對(duì)，我們統(tǒng)計(jì)了4月1—30日的誤差RMS分布情況(圖1)，由圖1可知，本文解算的TEC誤差RMS多在2 TECu(1 TECu=1016el/m2)以內(nèi)。

圖1 與CODE結(jié)果比對(duì)的RMS時(shí)序圖Fig.1Comparing RMS time serial with the TEC from CODE

3 基于Kalman濾波的電離層TEC異常探測(cè)

設(shè)離散Kalman濾波狀態(tài)方程［16］為:

式中Xk+1為狀態(tài)向量，Φk+1/k為轉(zhuǎn)移矩陣，Γk為系統(tǒng)驅(qū)動(dòng)矩陣，Yk是觀測(cè)向量，Hk是觀測(cè)系數(shù)向量，vk為觀測(cè)量噪聲，uk為狀態(tài)噪聲序列，均值為0，且為互不相關(guān)的白噪聲。由相應(yīng)的狀態(tài)方程和噪聲的統(tǒng)計(jì)特性求tk+1時(shí)刻狀態(tài)向量Xk+1的最佳估值，其算法詳見文獻(xiàn)［15］。

要實(shí)現(xiàn)Kalman濾波法逐步判別異常值算法，首先，應(yīng)確定初始狀態(tài)初始值及其方差矩陣，本文中初始值選取為觀測(cè)弧段的前10天TEC值，即用4月12—21日共10天的TEC數(shù)據(jù)作為Kalman濾波初值，然后進(jìn)行Kalman濾波逐步遞推，在遞推的每一步利用預(yù)報(bào)殘差對(duì)該觀測(cè)量進(jìn)行異常判別。記k+ 1時(shí)刻的觀測(cè)量為Yk+1，狀態(tài)量為Xk+1，觀測(cè)方程為成立時(shí)，對(duì)異常值Yk+1(i)的判別概率為95%，或者說對(duì)異常值Yk+1(i)的遺漏概率為5%。

4 汶川地震前電離層TEC異常的探測(cè)

選取臨震前一個(gè)月的電離層TEC資料進(jìn)行處理，結(jié)果顯示在震前數(shù)天多次出現(xiàn)了電離層TEC的異常擾動(dòng)。異常分布圖的映射方式參照文獻(xiàn)［5］，作為代表，圖2(a)和圖2(b)分別給出了震前第3天即5月9日08:00 UT和10:00 UT時(shí)中國(guó)區(qū)域上空附近的電離層TEC的異常分布。由圖2可以看出:5月9日在震中上空周圍的區(qū)域出現(xiàn)了明顯的電離層TEC異常擾動(dòng)現(xiàn)象，電離層TEC為異常增大，最大異常近5TECu，而且，隨著時(shí)間的推移，TEC異常擾動(dòng)的空間分布由東向西漂移。這與之前的結(jié)論相一致［6，7］。

圖2 中國(guó)區(qū)域上空的ΔTEC分布(單位:TECu)Fig.2Distribution of the ΔTEC value over China(unit:TECu)

為了更加清楚地了解孕震區(qū)上空電離層TEC異常擾動(dòng)隨時(shí)間變化的一般規(guī)律，我們給出了震中上空電離層ΔTEC隨時(shí)間變化的分布圖(圖3)，選取的考察時(shí)段為震前第20天至震后第10天，即2008年的4月22日—5月22日共31天，圖中橫軸代表時(shí)間，縱軸代表時(shí)刻，負(fù)號(hào)代表震前，0代表地震發(fā)生的當(dāng)天。

圖3 汶川Ms8.0地震震中上空電離層ΔTEC的時(shí)間分布圖(單位:TECu)Fig.3Distribution of the ΔTEC value with time over the epicenter of Wenchuan Ms8.0 earthquake(unit: TECu)

由圖3可以看出震中上空電離層TEC異常分布的一些主要特征:1)電離層TEC異常擾動(dòng)有正有負(fù)，其中，震前第13天和第6天震中上空的電離層ΔTEC呈現(xiàn)出明顯的異常減少現(xiàn)象，震前第3天TEC為異常增加;負(fù)異常擾動(dòng)的幅度不是很大，一般超出限值約3 TECu，而正異常擾動(dòng)的幅度在5月9日最大達(dá)到5 TECu。2)從電離層TEC異常的時(shí)間分布上來看，異常主要發(fā)生在震前第13天、6天和第3天，異常發(fā)生的時(shí)間段主要分布在下午至黃昏時(shí)段，即12:00－18:00(當(dāng)?shù)貢r(shí))，持續(xù)的時(shí)間約為4個(gè)小時(shí)，4月29日和5月6日TEC出現(xiàn)異常的時(shí)段為12:00－16:00，5月9日異常時(shí)段為14:00－18:00。采用同樣的處理方法，我們處理了同時(shí)段的IGS發(fā)布的全球電離層TEC格網(wǎng)數(shù)據(jù)，結(jié)果顯示在上述日期內(nèi)同樣出現(xiàn)了顯著的異常擾動(dòng)。3)從異常擾動(dòng)的起因來看，4月24日—5月25太陽活動(dòng)水平較低，4月29日、5月1日、2日、6日空間環(huán)境有低水平磁擾活動(dòng)［5］，可能會(huì)對(duì)電離層TEC產(chǎn)生一定的影響，因此，4月29日、5月6日電離層TEC的異常減少可能受地磁活動(dòng)的影響，但6日之后的一段時(shí)間地磁活動(dòng)平靜，空間環(huán)境趨于平靜，因此，排除空間環(huán)境及太陽活動(dòng)的影響，在5月9日震前孕震區(qū)上空附近出現(xiàn)的電離層TEC異常增加可能與此次地震有關(guān)。

為考察濾波時(shí)作為初始值的異常觀測(cè)對(duì)探測(cè)結(jié)果的影響，人為地對(duì)電離層TEC的初始觀測(cè)值做了一些修改，第一種方案為對(duì)某一天的所有時(shí)刻的觀測(cè)值的基礎(chǔ)上縮小為原來的50%，第二種方案為在原來觀測(cè)值的基礎(chǔ)上增大為原來的150%，修改后的觀測(cè)值都可以視為異常觀測(cè)，其探測(cè)方法不變，利用含有異常的觀測(cè)資料進(jìn)行異常探測(cè)，其結(jié)果如圖4所示。由圖4可以看出，震前第13天和第6天以及震前第3天TEC均出現(xiàn)了電離層TEC的異常擾動(dòng)，只是在個(gè)別天TEC異常擾動(dòng)的時(shí)間和強(qiáng)度有所不同，但5月9號(hào)出現(xiàn)的的電離層TEC異常增加的形態(tài)基本一致，由此可以看出:隨著時(shí)間的推移，先驗(yàn)異常值對(duì)探測(cè)結(jié)果的影響逐漸減弱。

圖4 先驗(yàn)異常值對(duì)異常探測(cè)結(jié)果的影響(單位:TECu)Fig.4Influence of abnormal TEC to distribution of the ΔTEC value(unit:TECu)

5 結(jié)果與討論

基于Kalman濾波探測(cè)的震前電離層TEC的異常擾動(dòng)，雖然獲得了很好的結(jié)果，但也有其不足之處:先驗(yàn)異常觀測(cè)值對(duì)鄰近的觀測(cè)值的探測(cè)結(jié)果會(huì)產(chǎn)生一定的影響，而且建模時(shí)需要對(duì)電離層TEC有一定的先驗(yàn)知識(shí)，特別是對(duì)濾波值的估算時(shí)，選取的模型噪聲和觀測(cè)噪聲協(xié)方差矩陣不同時(shí)，獲得的濾波值會(huì)有所差異，進(jìn)而影響到異常探測(cè)的結(jié)果。

致謝感謝中國(guó)地殼運(yùn)動(dòng)觀測(cè)網(wǎng)絡(luò)和IGS數(shù)據(jù)中心提供相關(guān)研究資料。

1卓裕榮.電離層地震前兆之研究［D］.國(guó)立中央大學(xué)太空科學(xué)研究所，2002.(Chuo Y J.A study of seismo－ionospheric precursor［D］.NCU Institute of Space Science，2002)

2張強(qiáng).地基GPS監(jiān)測(cè)電離層TEC及其與地震活動(dòng)關(guān)系的研究［D］.中國(guó)地震局地震研究所，2004.(Zhang Qiang.The ionospheric TEC detected by ground-based GPS and the relationship to the earthquake［D］.Institute of Seismology，CEA.2004)

3吳云，等.利用地基GPS探測(cè)震前電離層TEC異常［J］.大地測(cè)量與地球動(dòng)力學(xué)，2005，(2):36－40.(Wu Yun，et al.Preseismic ionospheric TEC anomaly detected by groundbased GPS［J］.Journal of Geodesy and Geodynamics，2005，(2):36－40)

4楊劍，吳云，周義炎.利用GPS無線電掩星數(shù)據(jù)研究地震前電離層異常［J］.大地測(cè)量與地球動(dòng)力學(xué)，2008，(1): 16－22.(Yang Jian，Wu Yun and Zhou Yiyan.Research on Seismo-ionospheric anomalies using GPS radio occultation data［J］.Journal of Geodesy and Geodynamics，2008，(2): 16－22)

5祝芙英，等.汶川地震前電離層VTEC的異常響應(yīng)［J］.地震學(xué)報(bào)，2009，31(2):180－187.(Zhu Fuying，et al.A-nomalous response of ionospheric VTEC before the Wenchuan earthquake［J］.Acta Seismologica Sina，2009，31 (2):180－187)

6趙必強(qiáng)，等.震前電離層擾動(dòng)研究進(jìn)展及汶川地震前電離層變化［J］.科技導(dǎo)報(bào)，2008，26(11):30－33.(Zhao Biqiang，et al.Recent advances on the ionospheric pre-cursors of earthquakes and ionospheric variations prior to Wenchuan earthquake［J］.Science and Technology Review，2008，26(11):30－33)

7Zhou Yiyan，et al.Ionospheric anomalies detected by ground－based GPS before the Mw7.9 Wenchuan earthquake，China［J］.Journal of Atmospheric and Solar-Terrestrial physics，2009，71:959－966.

8Yu Tao，et al.Study of the ionospheric anmaly before the Wenchuan earthquake［J］.Chinese Science Bulletin，2009，54(6):1 080－1 086.

9祝芙英，等.震前電離層TEC異常探測(cè)方法研究［J］.大地測(cè)量與地球動(dòng)力學(xué)，2009，(3):50－54.(Zhu Fuying，et al.Study on method of detecting ionospheric TEC anomaly before earthquake［J］.Journal of Geodesy and Geodynamics，2009，(3):50－54)

10崔希璋，等.廣義測(cè)量平差［M］.武漢:武漢測(cè)繪科技大學(xué)出版社，2005.(Cui Xizhang，et al.Generalized surveying adjustment［M］.Wuhan:Wuhan University of Surveying and Mapping Press，2005)

11張朝玉.卡爾曼濾波在多維AR序列建模中的應(yīng)用［J］.大地測(cè)量與地球動(dòng)力學(xué)，2003，(2):92－95.(Zhang Chaoyu.Application of Kalman filter in modeling multi-dimensional AR series［J］.Journal of Geodesy and Geodynamics，2003，(2):92－95)

12劉國(guó)林，等.卡爾曼濾波在InSAR噪聲消除與相位解纏中的應(yīng)用［J］.大地測(cè)量與地球動(dòng)力學(xué)，2006，(2):66－69.(Liu Guolin，et al.Application of Kalman filter to noise eliminating and phase unwrapping of InSAR［J］.Journal of Geodesy and Geodynamics，2006，(2):66－69)

13齊公玉，邱衛(wèi)寧，花向紅.爾曼濾波粗差修正方法應(yīng)用［J］.測(cè)繪工程，2010，19(2):50－52.(Qi Gongyu，Qiu Weining and Hua Xianghong.The method of gross error modification in Kalman filtering and its application［J］.Engineering of Survering and Mapping，2010，19(2):50－52)

14Wilson B D，et al.Subdaily northern hemisphere ionospheric maps using an extensive network of GPS receivers［J］.Radio Science，1999，30:639.

15Yuan Y B，et al.Preliminary results and analyses of using IGS GPS data to determine global ionospheric TEC［J］.Progress in Natural Science，2003，13(8):446.

16Kalman R E.A new approach to linear filtering and prediction proble［J］.Journal of Basic Engineering，1960，82(series D):35－45.

17付寧波.利用GPS跟蹤站數(shù)據(jù)解算TEC與地震電離層擾動(dòng)研究［D］.中國(guó)地震局地震研究所，2010.(Fu Ningbo.Determination of TEC based on GPS tracking station data and research on seismo-ionospheric disturbance［D］.Institute of Seismology，CEA 2010)

APPLICATION OF KALMAN FILTERING IN DETECTING IONOSPHERIC TEC ANOMALY PRIOR TO EARTHQUAKE

Nie Zhaosheng1，2)，Zhu Fuying1，2)and Fu Ningbo3)
1)Institute of Seismology，CEA，Wuhan430071

2)Crustal Movement Laboratory，Wuhan430071
3)Jiangxi Supervision and Inspection Center for Surveying and Mapping Product Quality，Nanchang330046

P318.2+4

A

1671-5942(2011)03-0047-05

2011-01-07

國(guó)家“863”課題(2007AA12Z169);中國(guó)地震局地震研究所所長(zhǎng)基金(IS200926039)

聶兆生，男，1982年生，碩士，主要從事GPS數(shù)據(jù)處理方面的相關(guān)研究.E－mail:niezhaosheng@126.com