鄭洪艷， 田 曉， 蘇廣利， 張 超

(中國地震局第一監(jiān)測中心，天津 300180)

基于粒子濾波的唐山臺跨斷層數(shù)據(jù)分析

鄭洪艷， 田 曉， 蘇廣利， 張 超

(中國地震局第一監(jiān)測中心，天津 300180)

基于游麗蘭等給出的跨斷層單測線動態(tài)系統(tǒng)模型，采用殘差重采樣粒子濾波算法對唐山臺跨斷層水準(zhǔn)(基線)觀測數(shù)據(jù)和降水量、氣溫等輔助觀測數(shù)據(jù)進行綜合處理。假設(shè)跨斷層形變觀測量僅受輔助觀測項的影響，可以根據(jù)殘差重采樣粒子濾波結(jié)果獲得環(huán)境因子(降水量、氣溫)對唐山臺跨斷層形變觀測量的影響和斷層運動情況。

粒子濾波；唐山臺；輔助觀測；跨斷層形變

0 引言

跨斷層形變測量通過直接測定活動斷層兩側(cè)參考點水平距離和相對高差的微小變化推斷斷層兩側(cè)介質(zhì)的三維運動方式[1]，是地震預(yù)測和地學(xué)研究的有效途徑。實際情況下，跨斷層形變監(jiān)測數(shù)據(jù)不僅包含斷層兩盤相對構(gòu)造運動信息，也包括大量氣象、地下水等因素造成的非構(gòu)造活動信息[2]。前人進行了大量探索，試圖利用多道維納濾波[3]、多元線性回歸[4]等方法排除干擾因素影響，提取跨斷層形變監(jiān)測中的地震前兆信息。游麗蘭等[5]結(jié)合跨斷層場地布設(shè)情況，系統(tǒng)建立了利用卡爾曼濾波處理跨斷層觀測數(shù)據(jù)的函數(shù)模型和隨機模型。

粒子濾波是基于貝葉斯采樣估計的重要性采樣思想發(fā)展起來的一種近似數(shù)值解方法。它的基本思想是利用狀態(tài)空間一組加權(quán)隨機樣本粒子逼近狀態(tài)的后驗概率分布[6-7]，不受模型狀態(tài)量和誤差高斯分布假設(shè)的約束，適用于任意非線性非高斯動態(tài)系統(tǒng)。為解決粒子濾波算法中存在的粒子退化問題，Gordon等[8]引入了重采樣方法，Douc等[9]對比并指出殘差重采樣算法的優(yōu)越性。

本文基于跨斷層單測線函數(shù)模型[5]，采用殘差重采樣粒子濾波算法對唐山臺的水準(zhǔn)(基線)觀測數(shù)據(jù)與降水量、氣溫等輔助觀測數(shù)據(jù)進行綜合處理。假設(shè)跨斷層形變觀測量僅受輔助觀測項的影響，可以根據(jù)殘差重采樣粒子濾波結(jié)果獲得環(huán)境因子對唐山臺跨斷層形變觀測量的影響和斷層運動情況。

1 跨斷層單測線動態(tài)系統(tǒng)模型

本文采用游麗蘭等[5]推導(dǎo)的跨斷層單測線動態(tài)系統(tǒng)模型，如式(1)、式(2)所示。

(1)

(2)

2 粒子濾波算法

2.1 貝葉斯估計

貝葉斯估計是在系統(tǒng)初始狀態(tài)、噪聲特征和測量信息已知的條件下，遞推獲得k時刻的狀態(tài)向量，是一個將先驗知識與測量數(shù)據(jù)加以綜合的過程，主要包括預(yù)測和更新2個步驟，如式(3)和式(4)所示：

p(xk|y1:k-1)=∫p(xk|xk-1)p(xk-1|y1:k-1)dxk-1，

(3)

(4)

其中，xk和yk分別表示到k時刻的狀態(tài)變量和觀測變量；p(xk-1|y1:k-1)為k-1時刻狀態(tài)變量的后驗概率分布；p(xk|xk-1)為狀態(tài)轉(zhuǎn)移概率；p(xk|y1:k-1)為k時刻預(yù)測得到的狀態(tài)變量先驗概率分布；p(yk|xk)為根據(jù)新的觀測數(shù)據(jù)計算得到的似然函數(shù)。實際工作中很難使用這種遞歸運算處理問題，式(3)和式(4)的閉環(huán)解通常無法獲得，因此應(yīng)用十分有限。

2.2 殘差重采樣粒子濾波算法

粒子濾波算法依據(jù)大數(shù)定理，采用蒙特卡洛方法求解貝葉斯估計中的積分運算。其核心思想是，通過尋找一組在狀態(tài)空間傳播的加權(quán)粒子來近似狀態(tài)變量的后驗概率分布，當(dāng)粒子數(shù)量N→∞時，可以逼近任何形式的概率密度分布。

(5)

其中，δ為Dirac函數(shù)。由于狀態(tài)變量后驗概率分布p(xk|y1:k)通常未知，因此采用容易抽樣的重要性分布函數(shù)q(xk|y1:k)進行抽樣。通常采用先驗概率分布函數(shù)作為重要性分布函數(shù)，重要性權(quán)值為：

(6)

(7)

為解決粒子濾波算法中存在的粒子退化問題引入了重采樣方法。本文采用殘差重采樣算法解決粒子退化問題，具體方法參考Zhang等[10]。

3 數(shù)據(jù)處理方法

本文采用粒子濾波進行數(shù)據(jù)處理的主要步驟包括：

6)判斷是否結(jié)束循環(huán)，若是則退出循環(huán)；否則令k=k+1，接收下一時刻測量值并執(zhí)行步驟2)。

4 處理結(jié)果與分析

唐山跨斷層臺站位于1976年唐山7.8級地震的主破裂帶上，其下與開平礦的5號斷層連接。2條水準(zhǔn)基線同樁觀測測線：水準(zhǔn)測線3-2(基線測線I-II)長約24 m，水準(zhǔn)測線4-1(基線測線III-IV)長約48 m(圖1)。臺站每天進行跨斷層短水準(zhǔn)、短基線和氣溫、降水量觀測。

4.1 唐山水準(zhǔn)3-2和基線I-II

唐山水準(zhǔn)3-2近年來變化形態(tài)呈現(xiàn)一定周期性，但每年周期形態(tài)略有差異：以2013年為界，2010—2016年年變周期幅度經(jīng)歷了不斷減小后逐步增大的過程，且在2013年以后年周變規(guī)律更加明顯并伴隨趨勢下降；2016年中旬曲線持續(xù)下降，破年變(圖2)。唐山基線I-II近年來形態(tài)多變：2012年以前呈波動狀態(tài)并具有一定年周期性，2012年中旬至2015年年變幅度明顯減小，2015年后年變消失，2016年中旬大幅下降后有所回升(圖3)。

圖1 唐山臺測線布設(shè)圖

對獲得的2010年2月至2016年10月近7年的唐山水準(zhǔn)3-2和基線I-II觀測數(shù)據(jù)與輔助觀測(降水量和氣溫)數(shù)據(jù)進行殘差重采樣粒子濾波處理，結(jié)果分別如圖2和圖3所示。

注：Yr 水準(zhǔn)3-2原始觀測值；PF 粒子濾波圖2 唐山水準(zhǔn)3-2濾波值時序圖

注：Yr 基線I-II原始觀測值；PF 粒子濾波圖3 唐山基線I-II濾波值時序圖

由于系統(tǒng)模型狀態(tài)向量同時包含斷層運動速率和環(huán)境因子影響參數(shù)，假設(shè)形變觀測量僅受輔助觀測項干擾，可以根據(jù)殘差重采樣粒子濾波結(jié)果分析環(huán)境因素(氣溫、降水量)對唐山水準(zhǔn)3-2和基線I-II跨斷層形變觀測的影響和斷層運動情況(圖4～7)。

a 氣溫及其形變觀測影響量

b 降水量及其形變觀測影響量圖4 唐山水準(zhǔn)3-2環(huán)境因子及其形變影響量

a 氣溫及其形變觀測影響量

b 降水量及形變觀測影響量圖5 唐山基線I-II環(huán)境因子及其形變影響量

圖6 唐山水準(zhǔn)3-2構(gòu)造形變位移累積量

圖7 唐山基線I-II構(gòu)造形變位移累積量

綜合圖4～7可以看出，唐山水準(zhǔn)3-2和基線I-II形變觀測量主要反映斷層兩盤相對構(gòu)造運動。環(huán)境因子中，氣溫對水準(zhǔn)3-2形變觀測量的影響幅度通常在0.01～0.02 mm之間，累積影響量絕對值不超過0.1 mm；對基線I-II形變觀測量的影響較為顯著，最大幅度可達0.05 mm，累積影響量絕對值可達0.13 mm。降水量對水準(zhǔn)3-2形變觀測量的影響幅度最大可達0.04 mm，累積影響量絕對值不超過0.05 mm；對基線I-II形變觀測量的影響幅度通常不大于0.01 mm，累積影響量也幾乎可以忽略。而相關(guān)性檢驗表明，唐山基線I-II與氣溫存在中等相關(guān)性(|r|=0.53)，與婁關(guān)壽等[2]、黃建平等[11]分析意見一致。

構(gòu)造形變位移累積量曲線(圖6～7)中顯示，在2012年唐山4.7級地震、2015年昌黎4.2級地震和2016年唐山4.1級地震前1～2年時間內(nèi)出現(xiàn)了破趨勢變化異常，說明唐山水準(zhǔn)3-2和基線I-II對其周圍100 km范圍內(nèi)4級左右地震有很好的響應(yīng)。

4.2 唐山水準(zhǔn)4-1和基線III-IV

與唐山水準(zhǔn)3-2相似，水準(zhǔn)4-1近年來的變化形態(tài)也呈現(xiàn)一定周期性：2012年以前，周期性表現(xiàn)不明顯，2013年以后年周變形態(tài)開始顯現(xiàn)，并伴隨趨勢下降(圖8)。唐山基線III～IV在2011—2013年期間波動幅度相對較大，2014年開始轉(zhuǎn)平轉(zhuǎn)穩(wěn)(圖9)。

注：Yr 水準(zhǔn)4-1原始觀測值；PF 粒子濾波圖8 唐山水準(zhǔn)4-1濾波值時序圖

注：Yr 基線III-IV原始觀測值；PF 粒子濾波圖9 唐山基線III-IV濾波值時序圖

對獲得的2010年2月至2016年10月近7年的唐山水準(zhǔn)4-1和基線III-IV觀測數(shù)據(jù)與輔助觀測(降水量和氣溫)數(shù)據(jù)進行殘差重采樣粒子濾波處理(圖8～9)。同樣假設(shè)形變觀測量僅受輔助觀測項干擾，可以根據(jù)殘差重采樣粒子濾波結(jié)果分析環(huán)境因素(氣溫、降水量)對唐山水準(zhǔn)4-1和基線III-IV跨斷層形變觀測的影響和斷層運動情況(圖10～13)。

a 氣溫及其形變觀測影響量

b 降水量及其形變觀測影響量 圖10 唐山水準(zhǔn)4-1環(huán)境因子及其形變影響量

a 氣溫及其形變觀測影響量

b 降水量及形變觀測影響量 圖11 唐山基線III-IV環(huán)境因子及其形變影響量

圖12 唐山水準(zhǔn)4-1構(gòu)造形變位移累積量

圖13 唐山基線III-IV構(gòu)造形變位移累積量

綜合圖10～13可以看出，唐山水準(zhǔn)4-1和基線III-IV形變觀測量主要反映斷層兩盤相對構(gòu)造運動。環(huán)境因子中，氣溫對水準(zhǔn)4-1形變觀測量的影響幅度通常在0.02～0.03 mm之間，累積影響量絕對值不超過0.15 mm；對基線III-IV形變觀測值的最大影響幅度可達0.3 mm，累積影響量絕對值可達0.6 mm。降水量對水準(zhǔn)4-1和基線III-IV形變觀測值的影響幅度最大可達0.04 mm，累積影響量絕對值不超過0.06 mm。

水準(zhǔn)4-1構(gòu)造形變位移累積量曲線(圖12)顯示，同樣在2012年唐山4.7級地震、2015年昌黎4.2級地震和2016年唐山4.1級地震前1～2年時間內(nèi)出現(xiàn)了破趨勢變化異常，說明唐山水準(zhǔn)4-1對其周圍100 km范圍內(nèi)4級左右地震有很好的響應(yīng)。

5 討論與結(jié)論

1)跨斷層單測線動態(tài)系統(tǒng)模型[5]狀態(tài)向量同時包含斷層運動速率和環(huán)境因子影響參數(shù)。利用跨斷層定點臺站的水準(zhǔn)(基線)觀測和降水量、氣溫等輔助觀測數(shù)據(jù)進行綜合分析，可以獲得環(huán)境因素(氣溫、降水量等)對跨斷層水準(zhǔn)(基線)形變觀測的影響量和斷層運動情況。

2)篩選臺站周圍100 km左右范圍內(nèi)4級以上地震，結(jié)合殘差重采樣粒子濾波得到的構(gòu)造形變位移累積量進行分析，可以看出唐山基線I-II和水準(zhǔn)3-2、水準(zhǔn)4-1對其周圍100 km范圍以內(nèi)4級左右地震也有很好的響應(yīng)。

致謝:衷心感謝所有審稿專家對本文提出的寶貴的修改意見。

[1] 中國地震局監(jiān)測預(yù)報司. 地震監(jiān)測預(yù)報40年[M]. 北京: 地震出版社, 2007.

[2] 樓關(guān)壽, 周偉, 金鵬, 等. 跨斷層形變觀測干擾因素的調(diào)查[J]. 大地測量與地球動力學(xué), 2010, 30(S2): 68-74.

[3] 吳大銘, 韓大宇. 用多道維納濾波方法處理唐山地震前后的大灰廠三種形變資料[J]. 地震學(xué)報, 1983, 5(1): 31-38.

[4] 王金明, 盧良玉. 桃花吐短水準(zhǔn)干擾因素的再研究[J]. 防災(zāi)減災(zāi)學(xué)報, 1998(2): 28-34.

[5] 游麗蘭, 劉大杰, 黃加納, 等. 跨斷層測量資料的卡爾曼濾波數(shù)學(xué)模型[J]. 中國地震, 1992(3): 44-52.

[6] Basar T. A new approach to linear filtering and prediction problems[M]//Basar T. Control Theory: Twenty-Five Seminal Papers (1932-1981). New York: Wiley-IEEE Press, 2001: 167-179.

[7] 李新, 擺玉龍. 順序數(shù)據(jù)同化的Bayes濾波框架[J]. 地球科學(xué)進展, 2010, 25(5): 515-522.

[8] Gordon N J, Salmond D J, Smith A F M. Novel approach to nonlinear/non-Gaussian Bayesian state estimation[J]. IEE Proceedings F (Radar and Signal Processing), 1993, 140(2): 107-113.

[9] Douc R, Cappe O. Comparison of resampling schemes for particle filtering[C]//Proceedings of the 4th International Symposium on Image and Signal Processing and Analysis. Zagreb, Croatia, Croatia: IEEE, 2005: 64-69.

[10] Zhang H J, Qin S X, Ma J W, et al. Using residual resampling and sensitivity analysis to improve particle filter data assimilation accuracy[J]. IEEE Geoscience and Remote Sensing Letters, 2013, 10(6): 1404-1408.

[11] 黃建平, 石耀霖, 孫玉軍, 等. 氣溫變化對唐山地震臺跨斷層形變觀測的影響[J]. 中國地震, 2012, 28(2): 222-230.

Analysis of Tangshan Cross-Fault Deformation Using Particle Filter

ZHENG Hong-yan, TIAN Xiao, SU Guang-li, ZHANG Chao

(First Crust Monitoring and Application Center, CEA, Tianjin 300180, China)

Based on the dynamic system model of across faults deduced by YOU Li-lan，using particle filter based on residual resampling method, we analyze the cross-fault leveling data of Tangshan station combing with auxiliary observations such as temperature, precipitation and so on. Assuming that the deformation observations are influenced only by these observed environmental factors (such as precipitation and the temperature), the influence quantity of these environmental factors and the fault activity can be obtained.

particle filter; Tangshan fixed-point station; auxiliary observation data; cross-fault deformation

2016-11-17

中國地震局監(jiān)測、預(yù)測、科研三結(jié)合課題(163302)；中國地震局第一監(jiān)測中心科技創(chuàng)新主任基金(FMC2016005)

鄭洪艷(1989—)，山東臨沂人，助理工程師，現(xiàn)主要從事跨斷層數(shù)據(jù)處理與分析工作.E-mail：hyzhengzheng@163.com

P315.7

A

1003-1375(2017)03-0039-06

10.3969/j.issn.1003-1375.2017.03.007