郝美仙 王鑫 張珂 劉穎 尹戰(zhàn)軍 張建中 鄭鈺

1)內(nèi)蒙古自治區(qū)地震局，呼和浩特 010010 2)中國地震局地球物理研究所，北京 100081

0 引言

地震是地球上主要的自然災(zāi)害之一，其突發(fā)性極強(qiáng)，破壞性巨大，一旦在城市或人口稠密地區(qū)附近發(fā)生破壞性地震，將會對重大基礎(chǔ)設(shè)施、生命線工程以及人民的生命財產(chǎn)造成嚴(yán)重危害。因此，科學(xué)有效地做好地震監(jiān)測和災(zāi)害防御工作，可以最大限度地減輕地震災(zāi)害所造成的損失(蔣長勝等，2016)。

地震預(yù)警是有效減輕地震災(zāi)害損失的有效手段之一，其原理是在地震發(fā)生后，搶在破壞性地震波到達(dá)之前，向可能遭受地震影響的目標(biāo)區(qū)發(fā)出警報(Zhang et al，2016；馬強(qiáng)，2008)，為相關(guān)行業(yè)緊急處置和人們逃生避險預(yù)留時間。地震預(yù)警信息發(fā)布的核心之一就是在破壞性地震發(fā)生過程中，能夠快速準(zhǔn)確地測定出潛在破壞區(qū)的范圍(張紅才等，2013)。

關(guān)于破壞性地震發(fā)生后如何快速估算預(yù)警目標(biāo)區(qū)的潛在破壞區(qū)，國內(nèi)外多名學(xué)者已基于地震動參數(shù)展開了相關(guān)研究(Teng et al，1997；馬強(qiáng)，2008；馬東，2013；宋晉東等，2017；王玉石等，2010；陳鯤等，2012)。地震動參數(shù)主要包括峰值加速度PGA、峰值位移PGD和峰值速度PGV等，由這些參數(shù)產(chǎn)出的潛在破壞區(qū)空間分布圖依賴于地震基本參數(shù)和地震事件持續(xù)時間，過程耗時幾分鐘到十幾分鐘，時效性不強(qiáng)，因此有必要采用一種時效性強(qiáng)的方法來快速獲取破壞性地震的潛在破壞區(qū)范圍，以便更好地服務(wù)于社會公眾及城市重要設(shè)施。

目前我國多地已建立了地震預(yù)警示范區(qū)，這些示范區(qū)依托于現(xiàn)有的地震觀測臺網(wǎng)，采用將實(shí)時數(shù)據(jù)流傳輸?shù)脚_網(wǎng)中心進(jìn)行統(tǒng)一處理的方式，對于地震預(yù)警算法的研究主要集中在震級估算與目標(biāo)區(qū)的破壞情況研究2個方面。對于震級估算，目前主要使用2種預(yù)警參數(shù)(位移幅值Pd和特征周期τc)，國內(nèi)外眾多學(xué)者對此展開了研究，并獲取到了適應(yīng)于特定地區(qū)的經(jīng)驗(yàn)關(guān)系式。對于目標(biāo)區(qū)的破壞情況研究，Zollo等(2010)曾提出一種綜合的地震預(yù)警方法，該方法實(shí)時聯(lián)合測量P波到達(dá)后3s窗內(nèi)的Pd與τc值，通過與預(yù)設(shè)的閾值相比較，達(dá)到快速區(qū)分是否為破壞性地震的目的；彭朝勇等(2019)利用該方法對四川廬山和云南魯?shù)?次地震的潛在破壞區(qū)范圍進(jìn)行計算，計算結(jié)果與實(shí)際烈度調(diào)查結(jié)果基本吻合。地震預(yù)警參數(shù)的經(jīng)驗(yàn)關(guān)系式具有區(qū)域相依性，與內(nèi)蒙古地區(qū)地殼速度結(jié)構(gòu)、斷裂帶分布及發(fā)震地區(qū)附近的臺站分布情況等因素有關(guān)。至今適合內(nèi)蒙古地區(qū)的地震潛在破壞區(qū)預(yù)警參數(shù)的經(jīng)驗(yàn)關(guān)系式的相關(guān)研究工作還未開展。

基于此，本文采用內(nèi)蒙古地區(qū)測震臺站記錄到的地震事件數(shù)據(jù)，擬合出適應(yīng)于內(nèi)蒙古地區(qū)的預(yù)警參數(shù)關(guān)系式，包括特征周期τc與震級ML、位移幅值Pd與τc和震源距R的統(tǒng)計關(guān)系；根據(jù)評估內(nèi)蒙古地區(qū)地震潛在破壞區(qū)范圍時定義的相應(yīng)預(yù)警參數(shù)閾值，利用2次已有地震開展線下模擬，對擬合關(guān)系式進(jìn)行驗(yàn)證。

1 數(shù)據(jù)選取與預(yù)處理

本研究選取2016—2019年內(nèi)蒙古及周邊地區(qū)ML≥2.0地震，共計251次，所選地震事件參數(shù)信息直接采用全國統(tǒng)一地震編目(1)中國地震臺網(wǎng)中心地震編目系統(tǒng). 內(nèi)部資料.結(jié)果，波形數(shù)據(jù)來源于內(nèi)蒙古測震臺網(wǎng)波形數(shù)據(jù)庫，文中所用震級標(biāo)度均為ML，所選事件和記錄臺站如圖1 所示。

圖1 本文選取的地震事件與臺站分布

由于所選該時間段內(nèi)發(fā)生的中強(qiáng)震事件較少，僅有1次ML5.4地震和1次ML6.1地震，因此在對獲取到的內(nèi)蒙古地區(qū)的測震波形數(shù)據(jù)進(jìn)行篩選處理時，采用了與彭朝勇等(2019)不同的規(guī)則，限制了震級下限、記錄臺站個數(shù)和震源距范圍。具體篩選條件為：①所選地震事件震級ML≥2.0；②有3個以上臺站可以清晰記錄到該事件；③震源距在200km范圍內(nèi)。通過以上條件共篩選出251個地震事件和3075條臺站記錄，本文所用震例分布情況見圖2。

圖2 本文所用震例分布情況統(tǒng)計(a)研究所用事件的震級-頻度分布；(b)研究所用記錄的震級-震源距分布

對于篩選后的臺站記錄數(shù)據(jù)，首先使用Matlab逐一打開各條記錄并手動選取P波和S波到時，保存至文本文件中，然后對各臺站垂直向速度記錄進(jìn)行積分后得到位移值，再利用Butterworth高通濾波器對位移值進(jìn)行連續(xù)濾波處理，最后分別根據(jù)P波觸發(fā)后3s窗段和全P波段垂直向位移時程和速度時程計算得出τc和Pd值，以進(jìn)行對比分析。另外，在處理過程中分別使用了2階和4階Butterworth高通濾波器進(jìn)行濾波處理，以便于考量不同階數(shù)濾波器對結(jié)果的影響。

2 相關(guān)參數(shù)的擬合關(guān)系式

2.1 τc與ML的統(tǒng)計關(guān)系式擬合

利用選取的所有地震事件記錄，采用下式進(jìn)行線性擬合

lg(τc)=A×ML+B

(1)

式中，A和B為待擬合系數(shù)，τc的單位為s。

在τc與ML的關(guān)系式擬合過程中，根據(jù)震級范圍和選取的P波不同窗長分別使用2階和4階Butterworth高通濾波器進(jìn)行數(shù)據(jù)擬合計算。對比結(jié)果顯示，ML2.0～6.1、ML2.5～6.1和ML3.0～6.1的記錄數(shù)據(jù)擬合獲得的結(jié)果，無論使用的濾波器階數(shù)和P波窗長為多少，標(biāo)準(zhǔn)偏差相對較大，相關(guān)性微弱，甚至出現(xiàn)了負(fù)相關(guān)的情況。通過對ML<3.5的地震事件記錄進(jìn)行分析發(fā)現(xiàn)，部分記錄在P波到達(dá)前存在較大的波形擾動情況，從而對計算獲得的τc值造成了極大的影響，由個別記錄算得的τc值超過了5s。這說明在利用τc參數(shù)開展震級估算時，記錄的數(shù)據(jù)需要具有較高的質(zhì)量。

震級選取在ML3.5～6.1和ML4.0～6.1的記錄數(shù)據(jù)，無論使用的濾波器階數(shù)和P波窗長為多少，均具備較好的相關(guān)性，尤其是對于震級范圍在ML4.0～6.1的記錄，其相關(guān)系數(shù)均在0.9以上。此外，使用4階Butterworth高通濾波器處理的結(jié)果要優(yōu)于使用2階濾波器，這主要?dú)w因于高通濾波器的過濾帶更窄，可以更好地濾除濾波頻帶范圍外的能量。圖3 為震級范圍為ML4.0～6.1時，不同窗長的P波段經(jīng)4階濾波器處理后所對應(yīng)的擬合曲線。

圖3 特征周期τc與震級ML的擬合曲線實(shí)心圓為每次事件所對應(yīng)各個臺站τc的平均值，空心圓為每次事件中各記錄臺站計算得出的τc值；黑色實(shí)線為τc與ML擬合曲線；黑色虛線為1倍標(biāo)準(zhǔn)方差；SDV為標(biāo)準(zhǔn)差；R為相關(guān)系數(shù)；Pole為所用的濾波器階數(shù)

2.2 Pd與τc和震源距R的統(tǒng)計關(guān)系

由于較小的地震不會引起破壞，一般不需要對其進(jìn)行預(yù)警處理，因此本文只分析ML4.0～6.1地震事件波形。對于所選用的震級范圍ML4.0～6.1的記錄數(shù)據(jù)，采用不同的垂直向P波窗長和4階Butterworth高通濾波器進(jìn)行處理后，經(jīng)過多變量線性擬合，所獲得的Pd與τc和震源距R的相關(guān)性分別為

lg(Pd)=1.5122+2.4665×lg(τc)-2.0336×lg(R)

(2)

lg(Pd)=0.3338+2.2531×lg(τc)-1.4228×lg(R)

(3)

式(2)、(3)分別為采用垂直向P波3s窗和全P波窗得出的擬合關(guān)系式，式中Pd的單位為cm；τc的單位為s；R為震源距。關(guān)系式對應(yīng)的擬合曲線如圖4，從圖中可以看出，Pd值明顯呈現(xiàn)出隨τc值增加而穩(wěn)步增加的趨勢，且利用垂直向全P波段數(shù)據(jù)得到的Pd值變化范圍更為穩(wěn)定。

圖4 Pd與τc和震源距R的相關(guān)性擬合曲線(a)P波3s窗；(b)全P波窗

3 臺站地震預(yù)警閾值設(shè)定

在開展?jié)撛谄茐膮^(qū)范圍估算時，需要定義相應(yīng)的預(yù)警參數(shù)閾值，以便實(shí)時跟蹤每個觸發(fā)臺站的變化情況。按照Wu等(2008)、Zollo等(2010)、彭朝勇等(2019)的閾值定義模式，需要同時獲取Pd和τc的閾值。對于Pd閾值，一般根據(jù)Pd與PGV的相關(guān)性以及PGV與儀器地震烈度的對應(yīng)關(guān)系來獲得。采用震級范圍ML4.0～6.1的記錄數(shù)據(jù)，同時采用4階Butterworth濾波器，針對不同的垂直向P波窗長擬合得到的相關(guān)性分別為：

垂直向P波3s窗

lg(PGV)=0.5977×lg(Pd)+0.6000±0.3717

(4)

垂直向全P波窗

lg(PGV)=0.7491×lg(Pd)+0.9736±0.3137

(5)

其中，PGV表示臺站三分向合成的地震動速度峰值，其單位為cm/s。

圖 5為Pd與PGV的相關(guān)性擬合曲線，由圖可以看出，無論是在標(biāo)準(zhǔn)偏差還是相關(guān)系數(shù)方面，利用垂直向全P波窗獲得的Pd與PGV的相關(guān)性要明顯優(yōu)于P波3s窗。因此，在后續(xù)的數(shù)據(jù)處理與分析中，所有的擬合參數(shù)關(guān)系式均選取震級范圍ML4.0～6.1，且采用4階Butterworth高通濾波器對垂直向全P波窗進(jìn)行分析。

圖5 Pd與PGV的相關(guān)性曲線(a)P波3s窗；(b)全P波窗

根據(jù)《儀器地震烈度計算暫行規(guī)程》(2)《儀器地震烈度計算暫行規(guī)程》(中震測發(fā)[2015]18號).內(nèi)部資料.，同時考慮到現(xiàn)有地震事件對應(yīng)的最大震級較低(ML6.1)，因此設(shè)定Ⅳ度對應(yīng)的三分向合成PGV下限值(1.77cm/s)作為對應(yīng)的閾值。根據(jù)式(5)可得到相對應(yīng)的Pd閾值為0.1075cm。當(dāng)設(shè)定ML=5時，得出相應(yīng)的τc值為0.686s。依據(jù)Pd和τc的閾值，可定義4類臺站地震預(yù)警級別：①3級：Pd≥0.1075cm且τc≥0.686s，說明正在發(fā)生的地震震級可能很大，并且發(fā)生在記錄臺站附近，預(yù)計會對記錄臺站附近和較遠(yuǎn)地區(qū)造成較大破壞，此為最高警報級別；②2級：Pd≥0.1075cm且τc<0.686s，說明地震發(fā)生在記錄臺站附近且震級較小，只可能對臺站周邊地區(qū)造成影響；③1級：Pd<0.1075cm 且τc≥0.686s，說明地震震級大但震中遠(yuǎn)離記錄臺站，僅在震中附近區(qū)域會引起破壞；④0級：Pd<0.1075cm且τc<0.686s，地震不會造成破壞。

在地震發(fā)生過程中，獲取到多個觸發(fā)臺站所對應(yīng)的Pd和τc實(shí)際測量值后，就可以采用與彭朝勇等(2013、2019)一致的方法，通過對各個臺站的τc值取平均值，同時將Pd值設(shè)定在閾值0.1075cm，根據(jù)式(3)可以計算出潛在破壞區(qū)的半徑。

在實(shí)際應(yīng)用時，首先對臺站覆蓋區(qū)進(jìn)行固定大小的單元網(wǎng)格劃分，單元格內(nèi)如果有觸發(fā)臺站，則直接利用記錄臺站Pd和τc的實(shí)測值；如果沒有觸發(fā)臺站，則使用式(3)和各臺站τc測量平均值來推算Pd值，對震中區(qū)域進(jìn)行內(nèi)插，以填補(bǔ)沒有記錄臺站的區(qū)域。生成的等值線Pd=0.1075cm即為評估破壞程度大小的臨界線。在中強(qiáng)地震發(fā)生時，利用該方法可以根據(jù)一定的時間間隔來實(shí)時估算破壞區(qū)范圍，并可隨著記錄臺站的增多逐步更新。

在使用式(3)對Pd值進(jìn)行計算時，需要先根據(jù)P波信息確定地震震中位置與震源深度，以此來計算震源距R值。目前國際和國內(nèi)用于地震預(yù)警的實(shí)時地震定位技術(shù)已有多種，如“著未著”定位方法(Horiuchi et al，2005)、RTLoc定位方法(Satriano et al，2008)、前三/四臺P波到時信息實(shí)時定位方法(張紅才等，2011)及四臺連續(xù)定位方法(李軍等，2014)等。由于地震定位方法不是本文主要的研究內(nèi)容，因此在進(jìn)行線下模擬時，將直接采用地震目錄中的震中信息和震源深度，以3臺觸發(fā)的1s時間點(diǎn)作為起始時間。

4 線下模擬與結(jié)果分析

利用內(nèi)蒙古測震臺網(wǎng)記錄到的阿拉善左旗2015年4月15日ML6.1和2017年6月3日ML5.4地震開展線下模擬(表1)。對記錄臺站覆蓋區(qū)域按照0.25°大小劃分為均勻的單元網(wǎng)格，每條事件記錄數(shù)據(jù)從P波觸發(fā)后開始，每隔1s計算一次Pd和τc值，直到S波到達(dá)后停止計算，線下模擬的2次地震事件相應(yīng)的潛在破壞區(qū)估算結(jié)果如圖6 所示，每個事件對應(yīng)的4張圖分別為發(fā)震后不同時刻計算得出的Pd分布圖。

表1 線下模擬的阿拉善左旗2次中強(qiáng)地震事件參數(shù)

圖6 2次地震在P波觸發(fā)后不同時間段的潛在破壞區(qū)實(shí)驗(yàn)結(jié)果五角星為震中位置；三角形為觸發(fā)臺站

從模擬結(jié)果可以看出，由于該地區(qū)測震臺站分布非常稀疏，無法獲取到地震發(fā)生后完整的潛在破壞區(qū)范圍。因此，在地震發(fā)生過程中采用該方法實(shí)時估算潛在破壞區(qū)時，需要臺網(wǎng)具有一定密度，如能均勻分布，效果更好，否則就會出現(xiàn)結(jié)果僅由經(jīng)驗(yàn)關(guān)系式?jīng)Q定、而實(shí)際觀測值沒有發(fā)揮作用的問題。

5 討論與結(jié)論

本文利用內(nèi)蒙古測震臺網(wǎng)記錄到的2016—2019年地震事件，按照篩選條件選取波形記錄，人工拾取P波和S波到時后，基于P波參數(shù)閾值開展評估內(nèi)蒙古地區(qū)地震潛在破壞區(qū)的研究。該方法所拾取的預(yù)警參數(shù)包括特征周期τc和位移幅值Pd，首先利用P波3s段和全P波段垂直向位移時程和速度時程計算τc和Pd值，研究中選用不同震級范圍、不同垂直向P波窗長和不同階數(shù)Butterworth高通濾波器分別得到內(nèi)蒙古地區(qū)τc與震級ML、Pd與τc和震源距R的統(tǒng)計關(guān)系式。

結(jié)果表明，在特征周期τc與震級ML的相應(yīng)關(guān)系擬合過程中，震級選取范圍在ML4.0～6.1 的事件記錄，對其P波3s窗和全P波窗波形數(shù)據(jù)進(jìn)行擬合，τc與ML均體現(xiàn)出了較好的相關(guān)性，且采用4階Butterworth高通濾波器處理的結(jié)果要優(yōu)于2階濾波器；而采用小于ML3.5的地震事件記錄擬合結(jié)果標(biāo)準(zhǔn)偏差均較大，相關(guān)性微弱，說明在采用τc參數(shù)開展震級估算時，需要波形記錄具有較高的數(shù)據(jù)質(zhì)量，且利用垂直向全P波窗數(shù)據(jù)獲取到的Pd值變化范圍更為穩(wěn)定。

震源距的計算受到震源深度不確定的影響，進(jìn)而會影響經(jīng)驗(yàn)關(guān)系的可靠性。利用已擬合得到的τc與震級ML、Pd與τc和震源距R的統(tǒng)計關(guān)系式，同時考慮到內(nèi)蒙古地區(qū)現(xiàn)有地震事件最大震級為ML6.1，按照震級ML=5.0，烈度為Ⅳ度定義了4類臺站地震預(yù)警級別，利用式(3)獲得了計算內(nèi)蒙古地區(qū)潛在破壞區(qū)的半徑范圍的關(guān)系式。需要注意的是，這里僅以Ⅳ度對應(yīng)的單水平方向PGV下限值計算Pd閾值，主要還是由內(nèi)蒙古地區(qū)現(xiàn)有地震記錄對應(yīng)的震級較小導(dǎo)致。今后在實(shí)際使用時，如果發(fā)生了震級較高的地震，則應(yīng)將PGV下限值提高，這樣會與實(shí)際計算獲得的潛在破壞區(qū)范圍更加吻合。

在地震發(fā)生過程中，依據(jù)地震預(yù)警的實(shí)時定位算法獲得震中位置、震級和震源深度，并通過獲取到多個觸發(fā)臺站所對應(yīng)的Pd和τc平均值后，計算潛在破壞區(qū)的半徑范圍。在計算時如只選取少量震中附近觸發(fā)臺站信息，則誤差較大，但時效性強(qiáng)；而選取更遠(yuǎn)、更多的觸發(fā)臺站信息，估算結(jié)果會更為準(zhǔn)確，但預(yù)警時間也會增加。

本文利用P波全波段的計算結(jié)果和經(jīng)驗(yàn)關(guān)系僅限于理論研究及線下模擬，效果較好，但在預(yù)警階段的實(shí)際作用有限，可主要用于在烈度速報階段配合實(shí)際儀器估算烈度的應(yīng)用，便于提高預(yù)警參數(shù)估算結(jié)果的準(zhǔn)確性。

本文研究獲得的內(nèi)蒙古地區(qū)預(yù)警參數(shù)經(jīng)驗(yàn)關(guān)系、閾值，由于區(qū)域臺站分布稀疏和選取地震事件的震級分布范圍不均勻且缺少大地震事件記錄數(shù)據(jù)，計算結(jié)果與其他地區(qū)存在差異。

隨著“國家地震烈度速報與預(yù)警工程”內(nèi)蒙古子項(xiàng)目的建設(shè)，在未來2～3年內(nèi)，臺網(wǎng)密度將會得到較大程度的提高，內(nèi)蒙古地震預(yù)警項(xiàng)目新建或改造共計61個基準(zhǔn)站、133個基本站和265個一般站。項(xiàng)目建成后將在呼和浩特、包頭、鄂爾多斯和烏蘭察布地區(qū)形成臺間距為14km的預(yù)警監(jiān)測網(wǎng)，在其他一般預(yù)警區(qū)形成臺站間距為60～70km的監(jiān)測網(wǎng)(圖7)。

圖7 “國家地震烈度速報與預(yù)警工程”內(nèi)蒙古子項(xiàng)目建成后臺網(wǎng)(a)和現(xiàn)有臺站分布(b)對比

將該方法應(yīng)用到內(nèi)蒙古近年來發(fā)生的阿拉善左旗2015年4月15日ML6.1和2017年6月3日ML5.4地震，結(jié)果顯示，由于該地區(qū)臺站分布稀疏，未得到潛在破壞區(qū)估算范圍，結(jié)果僅由經(jīng)驗(yàn)關(guān)系式?jīng)Q定，而實(shí)際觀測值沒有發(fā)揮作用。隨著“國家地震烈度速報與預(yù)警工程”內(nèi)蒙古子項(xiàng)目的建設(shè)，利用該方法實(shí)時開展?jié)撛谄茐膮^(qū)范圍估算所能取得的效果將會大大提高，精度將會得到明顯的改善。

致謝：感謝中國地震局地球物理研究所彭朝勇研究員為本文提供計算程序。