用于強(qiáng)震記錄仿真地動(dòng)位移的高精度累積數(shù)值積分方法研究1

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

1）山東省地震局泰安基準(zhǔn)地震臺(tái)，泰安 271000

2）山東省地震局，濟(jì)南 250014

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

4）山東省地震局安丘地震臺(tái)，濰坊 262100

5）山東省地震局煙臺(tái)地震臺(tái)，煙臺(tái) 264000

引言

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

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

本文從數(shù)值積分方法本身的研究著手，首先對(duì)記錄進(jìn)行分段三次Hermite插值，客觀上相當(dāng)于用類似上采樣的方式將采樣率變?yōu)樵紨?shù)據(jù)采樣率的4倍，之后采用積分區(qū)間四等分的Newton-Cotes積分公式進(jìn)行積分（李慶揚(yáng)等，2000），并對(duì)實(shí)際強(qiáng)震動(dòng)加速度記錄進(jìn)行了應(yīng)用驗(yàn)證。

1 積分方法

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

1.1 插值

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

根據(jù)強(qiáng)震動(dòng)記錄波形的性質(zhì)（鄭君里等，2011），選擇對(duì)原始數(shù)據(jù)保真性較強(qiáng)而且平滑較好的分段三次Hermite插值。

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

1.2 積分公式（李慶揚(yáng)等，2008）

根據(jù)Newton-Cotes公式：

從計(jì)算精度和計(jì)算量?jī)煞矫婢C合分析，本文選擇n=4時(shí)的Newton-Cotes公式：

1.3 誤差分析

當(dāng)n=4時(shí)，Newton-Cotes公式計(jì)算結(jié)果具有5次代數(shù)精度，截?cái)嗾`差為：

而當(dāng)n=1時(shí)，Newton-Cotes公式為：

即為梯形積分公式，計(jì)算結(jié)果具有1次代數(shù)精度，截?cái)嗾`差為：

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

2 方法驗(yàn)證

為了驗(yàn)證本文提出方法的積分效果，采用構(gòu)造的正弦信號(hào)和實(shí)際的強(qiáng)震動(dòng)加速度記錄分別進(jìn)行積分效果分析。

2.1 構(gòu)造正弦信號(hào)

構(gòu)造周期為2π，采樣率為10，長(zhǎng)度為4π的正弦信號(hào)，按照本文方法的處理步驟，首先對(duì)構(gòu)造信號(hào)進(jìn)行插值，見圖 2。其中，藍(lán)色為原始構(gòu)造信號(hào)，紅色為插值后的信號(hào)，為了能夠清楚分辨插值前后的變化，圖2（b）對(duì)橫軸［1，2］區(qū)間的信號(hào)進(jìn)行了放大。

圖3為構(gòu)造信號(hào)經(jīng)過一次和兩次積分后產(chǎn)生的信號(hào)對(duì)比。其中，藍(lán)色信號(hào)為構(gòu)造信號(hào)，紅色信號(hào)為經(jīng)過一次積分后的信號(hào)，綠色為經(jīng)過兩次積分后的信號(hào)。對(duì)正弦信號(hào)進(jìn)行一次積分，積分后信號(hào)為余弦信號(hào)，進(jìn)行兩次積分，積分后信號(hào)仍為正弦信號(hào)，但是變?yōu)樨?fù)值。圖3的結(jié)果與此吻合。

2.2 強(qiáng)震動(dòng)加速度記錄

為了驗(yàn)證本文累積數(shù)值積分方法對(duì)實(shí)際強(qiáng)震記錄的適用性，選取2013年4月20日雅安蘆山地震和2008年5月12日汶川地震的強(qiáng)震動(dòng)加速度記錄，采用原始梯形積分方法和本文方法分別對(duì)強(qiáng)震動(dòng)加速度記錄進(jìn)行積分，分別比較二者兩次積分的結(jié)果。

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

圖6為分別用本文方法和梯形積分方法對(duì)強(qiáng)震記錄進(jìn)行一次積分的結(jié)果對(duì)比。藍(lán)色為本文方法，紅色為梯形積分方法。圖6（b）為選取圖6（a）末尾0.2s數(shù)據(jù)進(jìn)行了局部放大。

圖7為分別用本文方法和梯形積分方法對(duì)強(qiáng)震記錄進(jìn)行兩次積分的結(jié)果對(duì)比。藍(lán)色為本文方法，紅色為梯形積分方法。圖7（b）為選取圖7（a）末尾0.2s數(shù)據(jù)進(jìn)行了局部放大。

圖2 構(gòu)造信號(hào)插值前后對(duì)比Fig. 2 Comparison of constructed signal before and after interpolation

圖3 構(gòu)造信號(hào)與一次和兩次積分后信號(hào)Fig. 3 Original structured signal and its integration signals after 1st and 2nd integration

圖4 寶興地辦強(qiáng)震臺(tái)記錄的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 雅安蘆山強(qiáng)震加速度記錄插值前后對(duì)比Fig. 5 Comparison of Yaan Lushan strong motion record before and after interpolation

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

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

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

圖8 汶川臥龍強(qiáng)震臺(tái)記錄的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 汶川強(qiáng)震加速度記錄插值前后對(duì)比Fig. 9 Comparison of strong motion record of Wenchuan earthquake before and after interpolation

圖10 為分別用本文方法和梯形積分方法對(duì)強(qiáng)震記錄進(jìn)行一次積分的結(jié)果對(duì)比。藍(lán)色為本文方法，紅色為梯形積分方法。圖10（b）為選取圖10（a）末尾0.2s數(shù)據(jù)進(jìn)行了局部放大。

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

圖11 為分別用本文方法和梯形積分方法對(duì)強(qiáng)震記錄進(jìn)行兩次積分的結(jié)果對(duì)比。藍(lán)色為本文方法，紅色為梯形積分方法。圖11（b）為選取圖11（a）末尾0.2s數(shù)據(jù)進(jìn)行了局部放大。

圖11 兩種積分方法對(duì)汶川強(qiáng)震記錄兩次積分結(jié)果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）來看，本文積分方法與梯形積分方法積分結(jié)果基本一致。但從圖6（b）、圖7（b）、圖10（b）和圖11（b）來看，兩種積分方法的結(jié)果有一定的差值，但是差值較小。我們認(rèn)為差值來自兩個(gè)方面：一是插值所引入的；二是積分方法的精度不同造成的。

經(jīng)過用兩種方法對(duì)11個(gè)強(qiáng)震臺(tái)站記錄的2013年4月20日雅安蘆山地震33個(gè)通道的加速度記錄和3個(gè)強(qiáng)震臺(tái)站記錄的2008年5月12日汶川地震9個(gè)通道的加速度記錄的兩次積分的結(jié)果對(duì)比，我們發(fā)現(xiàn)本文的積分方法與梯形積分方法積分結(jié)果均有一定差值，但差值較小。

3 討論

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

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

4 結(jié)論

本文提出了一種累積數(shù)值積分方法，首先對(duì)被積信號(hào)進(jìn)行分段三次Hermite插值，將插值后的信號(hào)長(zhǎng)度變?yōu)樵盘?hào)長(zhǎng)度的4倍，之后采用積分區(qū)間四等分的Newton-Cotes積分公式進(jìn)行積分。

經(jīng)過誤差分析以及對(duì)構(gòu)造信號(hào)和實(shí)際強(qiáng)震加速度信號(hào)的驗(yàn)證，本文積分方法的精度高于梯形積分方法，而且正確可靠，適用于將強(qiáng)震動(dòng)加速度記錄積分為位移記錄。

金星，馬強(qiáng)，李山有，2005. 利用數(shù)字強(qiáng)震儀記錄實(shí)時(shí)仿真地動(dòng)位移. 地震學(xué)報(bào)，27（1）：79—85.

李慶楊，關(guān)治，白峰杉，2000. 數(shù)值計(jì)算原理. 北京：清華大學(xué)出版社，251—282.

李慶揚(yáng)，王能超，易大義，2008. 數(shù)值分析. 北京：清華大學(xué)出版社.

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

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

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

王少英，2012. 數(shù)值積分若干方法的比較分析. 德州學(xué)院學(xué)報(bào)，28（6）：14—19.

于海英，江汶鄉(xiāng)，解全才等，2009. 近場(chǎng)數(shù)字強(qiáng)震儀記錄誤差分析與零線校正方法. 地震工程與工程振動(dòng)，29（6）：1—12.

鄭君里，應(yīng)啟珩，楊為理，2011. 信號(hào)與系統(tǒng). 北京：高等教育出版社.

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.