用于強震記錄仿真地動位移的高精度累積數值積分方法研究1

曲保安 劉希強 蔡 寅 范曉勇于慶民 趙銀剛 張 明 李 峰 孫 豪

1）山東省地震局泰安基準地震臺，泰安 271000

2）山東省地震局，濟南 250014

3）中國地震局地球物理研究所，北京 100081

4）山東省地震局安丘地震臺，濰坊 262100

5）山東省地震局煙臺地震臺，煙臺 264000

引言

地震預警用P波前3s記錄波形估算震級和烈度需要用到卓越周期和峰值位移等參數（Wu等，2006；2008），卓越周期的計算和峰值位移的量取都需要位移記錄，而目前的數字地震計為速度或加速度記錄，如何將速度記錄和加速度記錄積分為位移記錄是一個必須面對的難題，尤其是對于需要兩次積分的加速度記錄，其積分后的位移記錄的正確性和精確性一直被懷疑（Chen等，2007；Emore等，2007）。

關于利用強震加速度記錄積分為位移記錄的研究，早在20世紀40年代就已經開始。自1933年在加州長灘地震中記錄到第一條強震動加速度記錄以來，關于用強震加速度記錄估算近斷層同震位移的研究得到了深入的發展。但是，由于強震動記錄受到地面傾斜、環境噪音等各種因素的影響，直接使用強震動加速度記錄積分計算得到的速度和位移通常會出現嚴重的漂移失真現象，因此需要對加速度記錄進行專門的校正處理以便估算比較真實的同震位移，對于此方面的研究已經得到深入而廣泛的認同和長足的發展（金星等，2005；于海英等，2009；彭小波等，2011）。但是關于用于地震記錄數值積分方法本身的研究，卻未受到足夠的關注。

本文從數值積分方法本身的研究著手，首先對記錄進行分段三次Hermite插值，客觀上相當于用類似上采樣的方式將采樣率變為原始數據采樣率的4倍，之后采用積分區間四等分的Newton-Cotes積分公式進行積分（李慶揚等，2000），并對實際強震動加速度記錄進行了應用驗證。

1 積分方法

在對強震動記錄波形特征進行幅頻分析并對數值計算方法和數值積分以及累積數值積分的多種方法進行學習分析的基礎上，在目前對于累積數值積分只有梯形積分算法被廣泛應用的現狀，以及 Simpson、Romberg等多種數值積分算法和復合型數值積分方法得到深入研究的推動下，作者首先考慮采用n=4時的Newton-Cotes公式作為計算累積數值積分的方法，但是考慮到強震儀頻帶的高頻截至 80Hz以及近場強震動記錄信號覆蓋高頻部分，所以直接對原始強震記錄進行累積數值積分計算可能導致信號高頻損失過多，并且根據強震動記錄波形的性質（鄭君里等，2011），選擇分段三次Hermite插值（李慶楊等，2000；2008；劉彬清等，2002；龍愛芳等，2013）。

1.1 插值

目前應用較為普遍的一維插值方法有四種：最近鄰插值；線性插值；三次樣條插值；分段三次Hermite插值（李慶揚等，2008）。作者用周期為2π的正弦函數，選擇從0到10間隔為1的11個點，分別采用四種插值方法，插值后從0到10之間為41個點，繪制了四種插值方法的效果圖（圖1）。其中，圓圈為原始的11個點，菱形點為插值后的 41個點，線是為了直觀展示而將41個點連接的連線。圖1（a）為最近鄰插值；（b）為線性插值；（c）為三次樣條插值；（d）為分段三次Hermite插值。

根據強震動記錄波形的性質（鄭君里等，2011），選擇對原始數據保真性較強而且平滑較好的分段三次Hermite插值。

圖1 四種插值算法插值效果圖Fig. 1 Comparison of four interpolation algorithms

1.2 積分公式（李慶揚等，2008）

根據Newton-Cotes公式：

從計算精度和計算量兩方面綜合分析，本文選擇n=4時的Newton-Cotes公式：

1.3 誤差分析

當n=4時，Newton-Cotes公式計算結果具有5次代數精度，截斷誤差為：

而當n=1時，Newton-Cotes公式為：

即為梯形積分公式，計算結果具有1次代數精度，截斷誤差為：

從公式（5）和（6）來看，n=4時，Newton-Cotes公式積分方法的誤差比梯形積分精度高4個數量級。所以，梯形公式代數精度較低，Newton-Cotes公式的代數精度比梯形公式代數精度高，能更好地近似積分值（王少英，2012)

2 方法驗證

為了驗證本文提出方法的積分效果，采用構造的正弦信號和實際的強震動加速度記錄分別進行積分效果分析。

2.1 構造正弦信號

構造周期為2π，采樣率為10，長度為4π的正弦信號，按照本文方法的處理步驟，首先對構造信號進行插值，見圖 2。其中，藍色為原始構造信號，紅色為插值后的信號，為了能夠清楚分辨插值前后的變化，圖2（b）對橫軸［1，2］區間的信號進行了放大。

圖3為構造信號經過一次和兩次積分后產生的信號對比。其中，藍色信號為構造信號，紅色信號為經過一次積分后的信號，綠色為經過兩次積分后的信號。對正弦信號進行一次積分，積分后信號為余弦信號，進行兩次積分，積分后信號仍為正弦信號，但是變為負值。圖3的結果與此吻合。

2.2 強震動加速度記錄

為了驗證本文累積數值積分方法對實際強震記錄的適用性，選取2013年4月20日雅安蘆山地震和2008年5月12日汶川地震的強震動加速度記錄，采用原始梯形積分方法和本文方法分別對強震動加速度記錄進行積分，分別比較二者兩次積分的結果。

圖4為寶興地辦強震臺記錄的2013年4月20日雅安蘆山地震P波初至的前3s波形，對此強震加速度記錄先去除均值，再插值，見圖5。其中，藍色為原始強震記錄信號，紅色為插值后的信號。圖5（b）為更加清楚分辨插值前后變化，截取了圖5（a）前0.2s的數據。

圖6為分別用本文方法和梯形積分方法對強震記錄進行一次積分的結果對比。藍色為本文方法，紅色為梯形積分方法。圖6（b）為選取圖6（a）末尾0.2s數據進行了局部放大。

圖7為分別用本文方法和梯形積分方法對強震記錄進行兩次積分的結果對比。藍色為本文方法，紅色為梯形積分方法。圖7（b）為選取圖7（a）末尾0.2s數據進行了局部放大。

圖2 構造信號插值前后對比Fig. 2 Comparison of constructed signal before and after interpolation

圖3 構造信號與一次和兩次積分后信號Fig. 3 Original structured signal and its integration signals after 1st and 2nd integration

圖4 寶興地辦強震臺記錄的2013年4月20日雅安蘆山地震P波初至的前3s垂直向波形Fig. 4 The first three-second vertical component record of Yaan Lushan earthquake on Apr. 20, 2013 recorded at BXDB strong motion station

圖5 雅安蘆山強震加速度記錄插值前后對比Fig. 5 Comparison of Yaan Lushan strong motion record before and after interpolation

圖6 兩種積分方法對雅安蘆山強震記錄一次積分結果Fig. 6 The 1st integral results of strong motion records of Yaan Lushan earthquake with two different integration methods

圖7 兩種積分方法對雅安蘆山強震記錄兩次積分結果Fig. 7 The 2nd integral results of strong motion records of Yaan Lushan earthquake with two different integration methods

圖8 為汶川臥龍強震臺記錄的2008年5月12日汶川地震P波初至的前3s波形，對此強震加速度記錄先去除均值，再插值，見圖 9。其中藍色為原始強震記錄信號，紅色為插值后的信號。圖9（b）為更加清楚分辨插值前后變化，截取了圖9（a）前0.2s的數據。

圖8 汶川臥龍強震臺記錄的2008年5月12日汶川地震P波初至的前3s東西向波形Fig. 8 The first three-second horizontal east-west component record of Wenchuan earthquake on May. 12, 2008 recorded at WCW strong motion station

圖9 汶川強震加速度記錄插值前后對比Fig. 9 Comparison of strong motion record of Wenchuan earthquake before and after interpolation

圖10 為分別用本文方法和梯形積分方法對強震記錄進行一次積分的結果對比。藍色為本文方法，紅色為梯形積分方法。圖10（b）為選取圖10（a）末尾0.2s數據進行了局部放大。

圖10 兩種積分方法對汶川強震記錄一次積分結果Fig. 10 The 1st integral results of strong motion records of Wenchuan earthquake with two different integration methods

圖11 為分別用本文方法和梯形積分方法對強震記錄進行兩次積分的結果對比。藍色為本文方法，紅色為梯形積分方法。圖11（b）為選取圖11（a）末尾0.2s數據進行了局部放大。

圖11 兩種積分方法對汶川強震記錄兩次積分結果Fig. 11 The 2nd integral results of strong motion records of Wenchuan earthquake with two different integration methods

從圖6（a）、圖7（a）、圖10（a）和圖11（a）來看，本文積分方法與梯形積分方法積分結果基本一致。但從圖6（b）、圖7（b）、圖10（b）和圖11（b）來看，兩種積分方法的結果有一定的差值，但是差值較小。我們認為差值來自兩個方面：一是插值所引入的；二是積分方法的精度不同造成的。

經過用兩種方法對11個強震臺站記錄的2013年4月20日雅安蘆山地震33個通道的加速度記錄和3個強震臺站記錄的2008年5月12日汶川地震9個通道的加速度記錄的兩次積分的結果對比，我們發現本文的積分方法與梯形積分方法積分結果均有一定差值，但差值較小。

3 討論

本文僅限于積分算法本身的研究，然而從強震加速度記錄積分為位移記錄需要去除漂移，導致零線漂移的原因多種多樣，去除漂移的方法也各有不同，這些內容本文均未涉及。

用本文方法計算的結果與傳統梯形積分方法有一定的差值，但是由兩種積分方法的差值引起的卓越周期和峰值位移的差值大小，以及根據峰值位移和卓越周期計算的震級的差值大小尚未知，需要通過積累大量強震加速度記錄震例進行計算統計分析驗證。

4 結論

本文提出了一種累積數值積分方法，首先對被積信號進行分段三次Hermite插值，將插值后的信號長度變為原信號長度的4倍，之后采用積分區間四等分的Newton-Cotes積分公式進行積分。

經過誤差分析以及對構造信號和實際強震加速度信號的驗證，本文積分方法的精度高于梯形積分方法，而且正確可靠，適用于將強震動加速度記錄積分為位移記錄。

金星，馬強，李山有，2005. 利用數字強震儀記錄實時仿真地動位移. 地震學報，27（1）：79—85.

李慶楊，關治，白峰杉，2000. 數值計算原理. 北京：清華大學出版社，251—282.

李慶揚，王能超，易大義，2008. 數值分析. 北京：清華大學出版社.

劉彬清，任亞娣，2002．Newton-Cotes數值求積公式的漸近性. 上海大學學報（自然科學版），8（6）：503—506.

龍愛芳，胡軍浩，2013. 基于Hermite插值的高精度數值積分公式. 華僑大學學報（自然科學版），34（3）：349—352.

彭小波，李小軍，劉啟方，2011. 基于強震記錄估算同震位移的研究進展及方法. 世界地震工程，27（3）：73—80.

王少英，2012. 數值積分若干方法的比較分析. 德州學院學報，28（6）：14—19.

于海英，江汶鄉，解全才等，2009. 近場數字強震儀記錄誤差分析與零線校正方法. 地震工程與工程振動，29（6）：1—12.

鄭君里，應啟珩，楊為理，2011. 信號與系統. 北京：高等教育出版社.

Chen S.M., Loh C.H., 2007. Estimating permanent ground displacement from near-fault strong-motion accelerograms. Bulletin of the Seismological Society of America, 97 (1): 63—75.

Emore G.L., Haase J.S., Choi K. et a1., 2007. Recovering seismic displacements through combined use of 1-Hz GPS and strong-motion accelerometers. Bulletin of the Seismological Society of America, 97 (2): 357—378.

Wu Y.M., Kanamori H., 2008. Development of an earthquake early warning system using real-time strong motion signals. Sensor, 8:1—9.

Wu Y.M., Zhao L., 2006. Magnitude estimation using the first three seconds P-wave amplitude in earthquake early warning. Geophysical Research Letters, 33 (L16312): 1—4.