基于累積絕對位移值的震級估算方法

馬美帥，王延偉，汪祚賺，趙慶旭，潘代洪

（1.桂林理工大學(xué)廣西巖土力學(xué)與工程重點實驗室，廣西 桂林 541004；2.重慶南江工程勘察設(shè)計集團有限公司，重慶 401121）

引言

地震是一種無法避免的自然現(xiàn)象，具有極大破壞力，當(dāng)前科技無法準確地預(yù)測地震的發(fā)生，而地震預(yù)警（EEW）系統(tǒng)被認為是一種減輕地震災(zāi)害的有效的方法［1-2］。目前，世界上已有多個國家已經(jīng)或正在建立地震預(yù)警系統(tǒng)，如日本［3］，美國［4］，墨西哥［5］，土耳其［6］，羅馬尼亞［7］，意大利［8］，中國［9］等。地震預(yù)警系統(tǒng)是在地震發(fā)生后，利用先到達的P 波確定地震信息，在破壞性S 波到達前幾秒至數(shù)十秒向目標區(qū)域發(fā)出警報，使其能夠采取緊急措施，盡量減少地震破壞。震級作為反映地震大小的地震學(xué)基本參數(shù)，也是關(guān)鍵的地震預(yù)警信息，地震預(yù)警系統(tǒng)能否快速準確地估算震級是其成敗的關(guān)鍵。

地震預(yù)警中的震級估算，主要是通過初至P波的預(yù)警參數(shù)與震級的經(jīng)驗關(guān)系來實現(xiàn)的。目前，被廣泛研究和應(yīng)用的預(yù)警參數(shù)，可以分為三類。一類是與周期有關(guān)的參數(shù)，這類參數(shù)是根據(jù)震級越大地震波的周期越長提出的，主要有最大卓越周期τpmax［10］、平均卓越周期τc［11］、平均對數(shù)周期τlog［12］和基于小波變化的周期λlog［13］等。另一類是與地震波波形幅值有關(guān)的參數(shù)，這類參數(shù)主要依據(jù)震級越大相同震源距測點的地表變形越大，以Wu和Zhao［14］提出的位移幅值Pd的應(yīng)用和研究最為廣泛。還有一類是與地震輻射能量有關(guān)的參數(shù)，其依據(jù)是震級越大地震輻射能量越大，通過對地震波進行積分（求和）來近似求得地震輻射能量，主要有速度平方積分IV2［15］和累積絕對速度值CAV［16］。這三類預(yù)警參數(shù)都從不同方面代表了初至P 波中與震級相關(guān)的信息，與震級相關(guān)性越高的預(yù)警參數(shù)，越有利于震級的估算。為了從前述各預(yù)警參數(shù)中尋找估算震級準確性最好的參數(shù)，多位學(xué)者進行了對比分析，結(jié)果表明Pd參數(shù)與震級有最好的相關(guān)性，估算震級的準確性最高［17-22］。此外，多個已建地震預(yù)警系統(tǒng)的震級估算方法也主要采用Pd估算震級，例如中國的臺灣［23］和福建［24］，日本［25］，美國［26］和伊朗［27］。然而，從信息量上看，Pd參數(shù)僅保留了初至P波的最大變形，舍棄了初至P波中與斷層破裂過程有關(guān)的變形過程信息，并不利于震級的估算，仍有改進的空間。

為了更充分利用初至P波中與震級相關(guān)的信息，提高估算震級的準確性，本文提出利用累積絕對位移值（CAD）參數(shù)來估算震級。為驗證CAD參數(shù)的震級估算效果，利用大量的日本強震記錄（28018 條豎向記錄）的初至3 s ～6 s地震波，分析CAD參數(shù)與震級的相關(guān)性，以及CAD參數(shù)估算震級的誤差，并與Pd參數(shù)進行對比。結(jié)果表明，CAD參數(shù)估算震級的準確性高于Pd參數(shù)。

1 數(shù)據(jù)選擇及處理

日本Kiban-Kyoshin Net（KiK-Net）臺網(wǎng)記錄的大量歷史強震數(shù)據(jù)，被廣泛應(yīng)用于地震領(lǐng)域的各項研究。為保證數(shù)據(jù)規(guī)模，從數(shù)據(jù)庫中下載了1997年到2019年的地表地震強震記錄，并參考相關(guān)研究和地震預(yù)警系統(tǒng)的特點，按條件篩選強震記錄：震級≥4 級［20］，這里的震級為日本氣象廳（JMA）給出的震級，該震級在4～8級范圍內(nèi)與Ms 震級的差異可以忽略不計［28］；震源距范圍25 km～200 km，保證地震記錄至少有3 s P 波且盡可能包含近海地震事件；信噪比（SNR）大于10 dB，減小噪聲對震級估算的影響［12，29］；三分向合成加速度的峰值大于2 Gal［30］。此外，為保證數(shù)據(jù)質(zhì)量和使用，對下載的記錄進行了常規(guī)處理：檢查記錄的基線并進行校正；采樣頻率統(tǒng)一為100 Hz；自動撿拾P 波到時并人工校驗。經(jīng)過數(shù)據(jù)篩選和處理后，共獲得28 018 條豎向地震加速度記錄，作為驗證和對比的數(shù)據(jù)集，涉及3 519 次4～9 級地震事件，不同記錄數(shù)量隨震源距和震級的分布如圖1所示。

圖1 地震記錄數(shù)量隨震源距和震級的分布Fig.1 Distribution of seismic record number with source distance and magnitude

2 累積絕對位移值估算震級方法

2.1 預(yù)警參數(shù)CAD

地表的變形在一定程度上可以反映地震的大小，預(yù)警參數(shù)Pd是P波到達后幾秒內(nèi)位移的最大幅值，早在Wu等［14］就通過研究發(fā)現(xiàn)當(dāng)震源距已知時，可以通過Pd建立估算震級的經(jīng)驗公式。此后，多數(shù)研究也表明Pd估算震級的準確性要比其他參數(shù)好得多［17-22］，也被廣泛應(yīng)用于各個已建地震預(yù)警系統(tǒng)中［23-27］。然而，Pd參數(shù)僅能代表最大變形與震級的關(guān)系，并不能反映出變形過程（斷層破裂過程）與震級的相關(guān)性，為此，本文提出累積絕對位移值CAD，CAD與Pd參數(shù)相比，既包含了Pd參數(shù)所代表的最大變形信息，也包含了變形的過程信息，具有更多的與震級相關(guān)的信息。

在利用強震記錄計算CAD時，需要對加速度記錄進行2次積分得到位移記錄，再累積取絕對值的位移記錄，CAD的計算公式如式（1）：

式中，a（t）代表t時刻加速度（Gal），T代表了初至地震波的時長（s）。在CAD的計算過程中，為避免積分計算引入長周期誤差，采用與Pd參數(shù)計算過程中相同的數(shù)據(jù)處理方法［12］，即在每次積分后對記錄進行0.075 Hz的高通濾波（2階Butterworth）。

2.2 估算震級

震級的估算是通過建立CAD、震源距和震級的經(jīng)驗公式來實現(xiàn)的，經(jīng)驗公式的形式與Pd是相同的，如公式（2）所示：

式中，Para是預(yù)警參數(shù)（CAD或Pd）；M是震級；R是震源距；a，b，c是待定的回歸系數(shù)?；貧w系數(shù)是通過大量地震記錄的回歸分析得到，在利用不同時長初至地震波估算震級時，需要采用相應(yīng)時長初至地震波的回歸系數(shù)。

3 計算結(jié)果

3.1 初至3 s時P波的震級估算

地震預(yù)警系統(tǒng)屬于秒級數(shù)據(jù)處理和信息發(fā)布系統(tǒng)，在P波震相到達后，越快速地給出具有一定準確性的震級，就可以爭取到越多的預(yù)警時間。多數(shù)研究結(jié)果表明［11-12，14，31］，采用初至3 s P 波估算的震級較好地兼顧了地震預(yù)警的時效性和準確性，可以用于發(fā)布預(yù)警信息。因此，提高初至3 s P 波時估算震級的準確性尤為重要。利用前述歷史地震記錄數(shù)據(jù)集的初至3 s P波回歸分析公式（2），分別得到CAD和Pd與震級、震源距的經(jīng)驗公式（3）和（4）。為便于圖形展示和分析震級與CAD、Pd的相關(guān)性，利用震源距將CAD、Pd歸一化到震源距10 km 處［32］，震源距歸一化后，Pd用Pd10表示，CAD用CAD10表示。由圖2 可以看出，取對數(shù)后的CAD10和Pd10均與震級呈現(xiàn)出了較好的線性關(guān)系，與震級的相關(guān)系數(shù)CAD10略高于Pd10（CAD10為0.80 和Pd10為0.76）；CAD10和Pd10與震級的線性關(guān)系在不同的震級范圍內(nèi)存在差異，對于約6.5級以上地震事件，CAD10和Pd10總體上隨著震級的增大而減小，特別是7級以上地震尤為明顯，表現(xiàn)出了飽和的情況，這與多數(shù)學(xué)者所發(fā)現(xiàn)的“震級飽和”問題是一致的［11，14，22，32］。

圖2 CAD10、Pd10與震級的相關(guān)性Fig.2 Relationship between CAD10、Pd10 and magnitude

利用式（3）和式（4）估算震級，對比CAD和Pd估算震級的效果。圖3 為真實震級與估算震級的誤差分布：由圖可以看出用CAD估算震級的誤差的離散性略小，且誤差標準差更?。–AD為0.58，Pd為0.66）；CAD和Pd估算震級的誤差分布在不同震級范圍內(nèi)是不一致的，對于6.5級以上地震事件，也表現(xiàn)出了圖2中的“震級飽和”的問題。

圖3 估算震級的誤差分布Fig.3 Estimated magnitude error distribution

為了更直觀的展現(xiàn)CAD和Pd估算震級的誤差分布情況，作震級誤差分布直方圖（圖4）。由圖明顯看出CAD和Pd估算的震級誤差極少數(shù)分布在±1 個震級范圍外，少部分分布在（0.5，1）和（-1，-0.5）震級范圍之內(nèi)，大多數(shù)分布在［-0.5，0.5］震級范圍內(nèi)，且在［-0.5，0.5］震級范圍內(nèi)CAD的占比（誤差在震級范圍內(nèi)的記錄數(shù)量/總記錄數(shù)量*100%）約是Pd的1.2倍（CAD為63%，Pd為55%）。

圖4 震級誤差分布（百分數(shù)為占比）Fig.4 Magnitude error distribution

進一步對比CAD和Pd估算震級誤差隨震級的變化，圖5 為平均絕對誤差隨震級的變化情況：可明顯看出，兩參數(shù)的平均絕對誤差均隨震級的增大而增大，CAD在［4，6］級范圍內(nèi)估算震級的平均絕對誤差明顯小于Pd，在6 級以上時估算震級的平均絕對誤差與Pd差別很小。值得注意的是，7.3 級以上地震CAD和Pd估算震級的平均絕對誤差過大，超過了1個震級誤差。

圖5 平均絕對誤差Fig.5 Mean absolute error

3.2 在不同時長初至地震波時的震級持續(xù)估算

考慮到地震預(yù)警中的震級估算是隨著初至地震波時長的增加而不斷更新的過程，一般來說地震波時長越長越有利于震級估算的準確性。在前述結(jié)果的基礎(chǔ)上，進一步持續(xù)增加初至地震波時長檢驗CAD方法和Pd方法的持續(xù)震級估算效果。利用P波到達后的4 s、5 s和6 s 地震波求得參數(shù)CAD和Pd，通過最小二乘回歸分析式（2），分別得到CAD和Pd與震級的經(jīng)驗公式并估算震級，其回歸系數(shù)如表1所示。

表1 經(jīng)驗公式的回歸系數(shù)Table 1 Regression coefficients of empirical formulas

為了分析參數(shù)CAD和Pd與震級的相關(guān)程度和估算震級誤差的離散程度，作相關(guān)系數(shù)圖和誤差標準差圖進行對比（圖6）。圖6（a）表示震源距歸一化后的CAD和Pd分別與震級的相關(guān)系數(shù)隨時間的變化趨勢：從圖中可以看出，兩參數(shù)的相關(guān)系數(shù)均隨地震波時長的增加而有所提高，而CAD10與震級的相關(guān)系數(shù)整體上明顯高于Pd10，且CAD10在3 s、4 s、5 s和6 s時與震級的相關(guān)系數(shù)都達到0.8以上，呈高度相關(guān)。圖6（b）表示CAD和Pd估算震級的誤差標準差隨時間的變化情況：由圖可看出，兩者估算震級的誤差標準差均隨地震波時長的增加而降低，CAD的誤差標準差明顯小于Pd；總體來看，CAD10與震級的相關(guān)系數(shù)始終大于Pd10，CAD估算震級的誤差標準差始終小于Pd。

圖6 相關(guān)系數(shù)和誤差標準差Fig.6 Correlation coefficient and error standard deviation

進一步討論CAD和Pd估算震級的誤差分布情況。圖7給出了P波到達4 s、5 s和6 s時，CAD和Pd估算震級誤差的條形分布情況：由圖可以明顯看出，不同時長下，CAD和Pd估算的震級誤差少部分分布在（0.5，1）和（-1，-0.5）震級范圍之間，大多數(shù)分布在［-0.5，0.5］震級范圍內(nèi)；在［-0.5，0.5］震級范圍內(nèi)，CAD和Pd的占比隨著地震波時長的增加而增加，CAD的占比始終是Pd的1.1倍。

圖7 初至4 s～6 s地震波估算震級的誤差分布（百分數(shù)為占比）Fig.7 Error distribution of estimated magnitude using initial 4 s～6 s wave

進一步分析不同震級下持續(xù)估算震級的準確性。圖8代表P波初至4 s，5 s和6 s時，CAD和Pd估算震級的平均絕對誤差隨震級的變化情況，由于大震飽和使7.3級以上的估算出現(xiàn)明顯的震級低估，平均絕對誤差相比其他震級來說過大，不進行分析對比，故只選取7.3 級以下的地震數(shù)據(jù)研究。由圖8 可看出：在同一時刻，CAD與Pd的平均絕對誤差都隨震級的增大而增加；不同時刻間，二者的平均絕對誤差都隨地震波時長增加而減??；初至地震波時長4 s 時，CAD在［4-6］級范圍內(nèi)估算震級的平均絕對誤差要明顯比Pd更小，6 級以上兩者的平均絕對誤差相差不大；初至地震波時長5 s 和6 s 時，CAD在［4-6.3］級范圍內(nèi)估算震級的平均絕對誤差要明顯比Pd更小，6.3級以上兩者的平均絕對誤差相差不大。

圖8 初至4 s ～6 s地震波估算震級的平均絕對誤差Fig.8 Mean absolute error of estimated magnitude using initial 4 s～6 s wave

4 震例測試

以2019 年6 月18 日日本新潟發(fā)生的6.7 級地震為實例，評估CAD的震級估算效果。該地震為日本近海海域典型地震，震級為6.7級，震中位置（38.6°N，139.5°E），震源深度為14 km。地震造成了房屋破壞、列車停車、公路塌陷、滑坡、液化等問題。地震發(fā)生時，日本JMA 向公眾發(fā)布了地震預(yù)警信息［33］。圖9 為最先觸發(fā)4 個KiK-net 地表監(jiān)測臺站的CAD和Pd持續(xù)估算震級的結(jié)果（首個監(jiān)測臺站YMTH13 的P 波到時記為橫軸0 s），每個臺站都在P波到達后3 s，開始以1 s間隔持續(xù)估算震級，當(dāng)初至地震波時長為6 s時停止估算震級。由圖9 可以明顯看出，隨著地震時間的增長，CAD和Pd估算的震級整體呈上升趨勢；4 個臺站中CAD估算震級的最大偏差為+0.7個單位震級，Pd估算震級的最大偏差為+1.0個單位震級，且CAD估算震級的離散性比Pd更小。

圖9 持續(xù)估算震級Fig.9 Continuous magnitude estimation

進一步利用JMA 公布的地震預(yù)警信息對比分析震級估算效果。鑒于JMA 在首個監(jiān)測臺站P 波到時后的5.2 s，7.3 s，8.0 s，9.0 s 和9.9 s 持續(xù)更新了地震預(yù)警信息，取CAD和Pd在首個監(jiān)測臺站P 波到達后5.0 s-10.0 s時，各臺站估算震級的平均值作為最終估算震級進行對比。對比結(jié)果如圖10所示（首個監(jiān)測臺站的P波到時記為橫軸0 s）。從圖10中可以看出，在5 s，7 s和8 s時（6 s時，JMA 未給出預(yù)警信息），CAD估算的震級與真實震級最為接近，其次是Pd方法，JMA 估算的震級偏小了約0.6 個震級；在9 s 和10 s 時，JMA 估算震級偏高了約0.2 個震級，CAD偏高了0.4 個震級，Pd偏高了0.5 個震級。需要注意的是，CAD和Pd僅利用了KiK-net的4 個地表臺站的數(shù)據(jù)估算震級，而JMA 利用更多監(jiān)測臺站（如KiK-net井下臺站、K-net臺站、Hi-net臺站等）的數(shù)據(jù)估算震級?？偟膩砜矗珻AD在本次地震中估算震級的偏差是小于JMA和Pd的估算結(jié)果的。

圖10 持續(xù)估算震級平均值Fig.10 The average of estimated magnitude

5 結(jié)論

（1）為提高地震預(yù)警中震級估算的準確性，提出了一種新的震級估算方法，即利用初至地震波的CAD建立震級的經(jīng)驗公式。CAD的計算過程是，通過對加速度記錄進行兩次積分得到位移記錄，再對位移記錄取絕對值后求累積。在這個過程中，CAD利用了初至地震波引起的地表變形的所有信息，間接地包含了斷層的初始破裂過程信息。相比于廣泛使用的Pd參數(shù)（僅利用地表最大變形信息），CAD在理論上可以保留更多與震級估算相關(guān)的信息。

（2）利用日本KiK-Net 數(shù)據(jù)庫中28018 條地表強震觀測記錄（震級范圍4～9 級）對CAD估算震級效果進行了驗證，并與當(dāng)前普遍采用的Pd方法進行了對比。結(jié)果表明，對于4～7.3 級地震，初至地震波時長為3 s～6 s 時，CAD估算震級的效果隨著地震波時長的增加不斷改善，并且與震級的相關(guān)性比Pd大，估算震級的誤差標準差比Pd小，估算震級誤差在［-0.5，0.5］震級范圍內(nèi)的占比是Pd的1.1～1.2倍。這證明了，利用變形過程信息的CAD比僅利用最大變形信息的Pd更有利于震級的估算，CAD可以為地震預(yù)警系統(tǒng)提供更準確的震級估算。

（3）在日本近海6.7 級典型地震事件的測試中，CAD估算震級的偏差總體上是小于JMA 和Pd的估算結(jié)果，這表明CAD用于實際地震預(yù)警系統(tǒng)具有可行性。

致謝：

感謝日本防災(zāi)科學(xué)技術(shù)研究所（NIED，https：//www.doi.org/10.17598/NIED.0004）提供了KiK-net 強震數(shù)據(jù)下載服務(wù)。