陳向陽,于金池,葛 建,侯勇濤

(1. 南通職業(yè)大學(xué)建筑工程學(xué)院，江蘇 南通 226007； 2. 南通大學(xué)交通與土木學(xué)院，江蘇 南通 226019； 3. 上海華測導(dǎo)航技術(shù)股份有限公司，上海 201702)

近些年來，GPS被廣泛地運用于長周期地殼形變、火山監(jiān)測和臺風(fēng)預(yù)警[1-7]。這些監(jiān)測主要依靠于毫米級的GPS坐標時間序列，與此同時，GPS站的采樣間隔通常設(shè)置為30 s[8]。地震及其引起的火山爆發(fā)、滑坡和泥石流等對人類的生產(chǎn)和生活影響巨大，因此，快速準確地監(jiān)測地震具有尤為重要的意義。隨著GPS定位精度與方法的成熟，高采樣率的GPS數(shù)據(jù)(如1、10、50 Hz)已經(jīng)被用于監(jiān)測由地震引起的瞬時地殼形變。高頻GPS被成功運用于監(jiān)測2002年Denali地震、2008年汶川地震和2011年日本宮城9.0級地震等[9-11]。之后，高頻GPS結(jié)合其他手段也被成功應(yīng)用于2003年Miyagi地震、2009年L’Aquila地震和2015年尼泊爾地震[12-13]。鑒于高頻GPS與強震儀互有優(yōu)缺點，因此可以融合GPS與強震儀的數(shù)據(jù)來實時監(jiān)測地震引起的動態(tài)形變。本文以2016年意大利中部的Mw 6.0地震為例，探討融合GPS數(shù)據(jù)和強震儀數(shù)據(jù)實時監(jiān)測中等強度地震引起瞬時地殼形變的可行性與優(yōu)勢。

1 數(shù)據(jù)與方法

1.1 GPS數(shù)據(jù)

2016年8月24號，意大利中部發(fā)生了Mw 6.0級地震，震中位于Perugia，震源深度為8.1 km。目前，意大利境內(nèi)已有完善的高密度GNSS觀測網(wǎng)絡(luò)，主要包括RING、ISPRA、DPC、Regione Lazio、Regione Abruzzo、Leica ITALPOS和Topcon NETGEO，使探測中等強度引起的同震動態(tài)形變成為可能。此次地震引起的同震形變成功地被近場的10 Hz高頻GPS和200 Hz強震儀所記錄。為了融合GPS與強震儀數(shù)據(jù)實時監(jiān)測地震引起的地殼形變，筆者挑選了3組位置近似相同的GPS/強震儀組合進行分析(即RM33/MTER、GSA/INFN和AMT/AMAT)，GPS與強震儀之間的相對距離少于1.1 km。此次所用的高頻GPS數(shù)據(jù)由INGV提供，強震儀數(shù)據(jù)由美國強震運動監(jiān)測中心(CESMD)提供。如圖1所示。

1.2 數(shù)據(jù)融合方法

鑒于強震儀記錄的加速度二次積分求解的位移波形不理想，而高頻GPS由于采樣率的問題不能夠更為細致地描述位移波形，許多學(xué)者嘗試將GPS與強震儀數(shù)據(jù)結(jié)合來監(jiān)測地震的同震動態(tài)形變。文獻[14]于2011年提出了利用Kalman濾波將GPS與強震儀數(shù)據(jù)融合來獲取準確和細致的位移波形。

將以上Kalman濾波應(yīng)用于測站地表位移波形融合時，假設(shè)其在N、E、U 3個方向上均為一維的運動形式，可用一階線性微分方程來描述其運動狀態(tài)，即

(1)

對上述一階線性微分方程離散化，其狀態(tài)方程系數(shù)陣會發(fā)生變化，若離散數(shù)據(jù)采樣間隔為τ，采用離散GPS位移觀測量zk和強震儀加速度ak進行Kalman濾波融合，可建立如下的離散狀態(tài)方程與觀測方程

(2)

(3)

1.3 交叉小波變換

首先介紹連續(xù)小波變換(CWT)，一組時間序列Xn的連續(xù)小波變換如下

(4)

WXY=WXWY*

(5)

式中，Y*為WY的共軛矩陣。平均相位角可以更好地描述兩個時間序列的相位關(guān)系。一組相位角(ai,i=1,2，…，n)的平均相位角可以表示為

(6)

對am進一步處理可以更直觀地描述兩組時間序列之間的相關(guān)性。交叉小波相似值可以用下式來計算

ρ=cosαm

(7)

式中，ρ=1表示兩組時間序列完全相關(guān)；ρ=0表示兩組時間序列沒有明顯的相關(guān)性；ρ=-1表示兩組時間序列呈負相關(guān)。

2 結(jié)果與分析

2.1 GPS與強震儀位移波形對比

由于此次地震持續(xù)時間較短，僅選取震前10 s到震后60 s這個時間段的坐標時間序列進行分析。與GPS站不同的是，強震儀記錄的往往是地震引起的加速度變化值，需要一個二次積分的過程才能夠得到震時的位移波形。對強震儀記錄的加速度進行二次積分后，利用SeismoSignal軟件進行基線校正和低通濾波。圖2為GPS位移波形和強震儀二次積分得到的位移波形的對比。由圖2可知，GPS/強震儀組合RM33/MTER在水平分量上符合較好，兩種手段獲取的位移波形差別在2 cm以內(nèi)；由于垂直方向的精度問題，GPS/強震儀組合在垂直方向上位移波形符合度較差。

2.2 GPS與強震儀數(shù)據(jù)融合

Kalman濾波可以利用GPS和強震儀的優(yōu)點，進而補充GPS和強震儀各自的缺陷。融合后的位移波形可以更為細致地描述地震引起的動態(tài)地表形變，同時也可以探測地震引起的同震永久階躍。以GPS/強震儀組合RM33/MTER為例，利用Kalman濾波進行數(shù)據(jù)的融合。從圖3可以清楚看到，Kalman濾波獲取的位移波形與GPS位移波形符合很好，與此同時，地震引起的同震階躍在N方向上可以明顯看到。GPS位移數(shù)據(jù)可以作為融合數(shù)據(jù)的整體趨勢，而細節(jié)部分則用高采樣率的強震儀數(shù)據(jù)進行描述。Kalman濾波獲取的位移波形相對于GPS位移波形更為平滑，說明融合GPS與強震儀數(shù)據(jù)能夠獲取更為精確的位移波形。盡管GPS垂向方向精度較差，但是在強震儀數(shù)據(jù)的輔助下Kalman濾波仍然能夠獲取精確的位移波形。

2.3 交叉小波分析

為了進一步驗證Kalman濾波位移波形與GPS、強震儀位移波形的相關(guān)性，筆者采用交叉小波變換進行任意2個位移波形之間的相關(guān)性分析。在交叉小波譜中，黑色粗線包圍的區(qū)域表示通過了95%的顯著性檢驗。兩組時間序列的相對相位關(guān)系如交叉小波譜中的箭頭所示(向右表示相位關(guān)系相同，向左表示相位關(guān)系相反，向上表示相位落后90°變化，向下表示相位超前90°變化)。

圖4為GPS/強震儀組合RM33/MTER的強震儀與Kalman濾波位移波形、強震儀與GPS位移波形、GPS與Kalman濾波位移波形的交叉小波譜。對于GPS與強震儀位移波形的交叉小波譜(如圖4(b)所示)，共振區(qū)域內(nèi)的箭頭幾乎全部朝右，表明GPS與強震儀位移波形在地震發(fā)生期間呈明顯的正相位關(guān)系。GPS與強震儀在共振區(qū)域的交叉小波相似值為0.99，表明二者在地震時呈明顯的正相關(guān)性。對于強震儀、GPS與Kalman濾波位移波形的交叉小波譜(如圖4(a)和(b)所示)，共振區(qū)域內(nèi)的部分箭頭向下有一定的偏移，說明Kalman濾波位移波形比GPS、強震儀位移波形的相位有一定的超前。圖4(a)和(b)的共振區(qū)域的交叉小波相似值均在0.8以上，說明Kalman濾波波形與GPS、強震儀位移波形均呈較強的正相關(guān)性。對于3個組合，共振區(qū)域的小波相似值，GPS與強震儀位移波形的交叉小波相似值均在0.9以上，而Kalman濾波與GPS、強震儀位移波形的交叉小波相似值均在0.8以上。交叉小波分析表明，Kalman濾波獲取的位移波形與GPS、強震儀位移波形均呈較強的正相關(guān)性，表明Kalman濾波可能獲取準確的位移波形。

3 結(jié) 語

本文以2016年意大利Mw 6.0級地震為例，基于3組位置近似相同的GPS/強震儀組合，探討了利用Kalman濾波融合GPS與強震儀數(shù)據(jù)監(jiān)測中等強度地震引起的瞬時地殼形變的可行性與優(yōu)勢。結(jié)果表明Kalman濾波能夠有機地結(jié)合GPS與強震儀數(shù)據(jù)的優(yōu)點，融合后的位移波形精度更高，并能夠?qū)Φ卣鹨鸬膭討B(tài)地殼形變的細節(jié)部分進行較好地描述。然后利用交叉小波變換，對融合后的位移波形與GPS、強震儀的位移波形在時頻域的相關(guān)性進行了分析。交叉小波分析表明，Kalman濾波位移波形與GPS、強震儀位移波形呈較強的正相關(guān)性，說明Kalman濾波確實獲取了比較穩(wěn)定可靠的位移波形。由于GPS與強震儀位置不完全一致，Kalman濾波獲取的位移波形與GPS、強震儀位移波形會有略微的相位差異。