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

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

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

引言

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

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

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

1 數據選擇及處理

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

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

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

2.1 預警參數CAD

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

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

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

2.2 估算震級

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

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

3 計算結果

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

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

圖2 CAD10、Pd10與震級的相關性Fig.2 Relationship between CAD10、Pd10 and magnitude

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

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

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

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

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

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

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

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

表1 經驗公式的回歸系數Table 1 Regression coefficients of empirical formulas

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

圖6 相關系數和誤差標準差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）震級范圍之間，大多數分布在［-0.5，0.5］震級范圍內；在［-0.5，0.5］震級范圍內，CAD和Pd的占比隨著地震波時長的增加而增加，CAD的占比始終是Pd的1.1倍。

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

進一步分析不同震級下持續估算震級的準確性。圖8代表P波初至4 s，5 s和6 s時，CAD和Pd估算震級的平均絕對誤差隨震級的變化情況，由于大震飽和使7.3級以上的估算出現明顯的震級低估，平均絕對誤差相比其他震級來說過大，不進行分析對比，故只選取7.3 級以下的地震數據研究。由圖8 可看出：在同一時刻，CAD與Pd的平均絕對誤差都隨震級的增大而增加；不同時刻間，二者的平均絕對誤差都隨地震波時長增加而減??；初至地震波時長4 s 時，CAD在［4-6］級范圍內估算震級的平均絕對誤差要明顯比Pd更小，6 級以上兩者的平均絕對誤差相差不大；初至地震波時長5 s 和6 s 時，CAD在［4-6.3］級范圍內估算震級的平均絕對誤差要明顯比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 日日本新潟發生的6.7 級地震為實例，評估CAD的震級估算效果。該地震為日本近海海域典型地震，震級為6.7級，震中位置（38.6°N，139.5°E），震源深度為14 km。地震造成了房屋破壞、列車停車、公路塌陷、滑坡、液化等問題。地震發生時，日本JMA 向公眾發布了地震預警信息［33］。圖9 為最先觸發4 個KiK-net 地表監測臺站的CAD和Pd持續估算震級的結果（首個監測臺站YMTH13 的P 波到時記為橫軸0 s），每個臺站都在P波到達后3 s，開始以1 s間隔持續估算震級，當初至地震波時長為6 s時停止估算震級。由圖9 可以明顯看出，隨著地震時間的增長，CAD和Pd估算的震級整體呈上升趨勢；4 個臺站中CAD估算震級的最大偏差為+0.7個單位震級，Pd估算震級的最大偏差為+1.0個單位震級，且CAD估算震級的離散性比Pd更小。

圖9 持續估算震級Fig.9 Continuous magnitude estimation

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

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

5 結論

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

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

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

致謝：

感謝日本防災科學技術研究所（NIED，https：//www.doi.org/10.17598/NIED.0004）提供了KiK-net 強震數據下載服務。