于海英 徐 旋 張同宇

?

強(qiáng)震動記錄的低頻截止頻率自動搜索算法1

于海英1，2）徐 旋1，2）張同宇1，2）

1）中國地震局工程力學(xué)研究所，哈爾濱 150080 2）中國地震局地震工程與工程振動重點(diǎn)實(shí)驗(yàn)室，哈爾濱 150080

本文針對傳統(tǒng)方法在強(qiáng)震動記錄處理中確定濾波低頻截止頻率效率較低的問題，提出一種確定濾波低頻截止頻率的自動搜索模型，并利用統(tǒng)計(jì)學(xué)習(xí)方法中的損失函數(shù)確定自動搜索模型流程的結(jié)束條件?；?008年汶川和2013年蘆山兩次地震主震和余震獲得的強(qiáng)震動記錄，將自動搜索算法得出的結(jié)果與傳統(tǒng)方法確定的低頻截止頻率進(jìn)行比較，分析自動搜索算法產(chǎn)生誤差的原因，進(jìn)而提出自動搜索算法應(yīng)遵循的原則和適用條件。結(jié)果表明該算法計(jì)算效率比傳統(tǒng)方法有很大提升，特別適合海量強(qiáng)震動記錄的批量處理。

強(qiáng)震動記錄 濾波 低頻截止頻率 自動搜索模型 損失函數(shù)

引言

強(qiáng)震動記錄處理是地震工程與工程地震相關(guān)研究的基礎(chǔ)工作，可以為抗震設(shè)防、地震烈度速報和預(yù)警以及重要工程強(qiáng)震動觀測等工作提供高質(zhì)量的科學(xué)數(shù)據(jù)。然而，受環(huán)境振動干擾等因素影響，強(qiáng)震動記錄不可避免地隱含噪聲。因此，研究如何對強(qiáng)震動記錄中的噪聲進(jìn)行識別尤為重要。而無論是過去還是將來，頻率域去噪都是一項(xiàng)主要的去噪技術(shù)（張軍華等，2006），所以濾波技術(shù)的研究在強(qiáng)震動記錄處理中占有非常重要的位置。

濾波器的選擇是強(qiáng)震動記錄處理中的核心問題之一。不同強(qiáng)震動記錄處理機(jī)構(gòu)選擇的濾波器有所不同。USGS采用雙向巴特沃斯非因果濾波器；CSMIP使用Ormsby或低通巴特沃斯濾波器進(jìn)行非因果濾波；COSMOS采用因果與非因果結(jié)合的多種濾波器類型。根據(jù)張同宇（2016）的研究，對于強(qiáng)震動記錄事件分析推薦使用巴特沃斯（IIR）非因果濾波，對于強(qiáng)震動記錄實(shí)時處理推薦使用巴特沃斯（IIR）因果濾波，在分析蘆山地震強(qiáng)震動記錄時推薦使用4階巴特沃斯（IIR）濾波。本文采用雙向巴特沃斯非因果濾波器，濾波階數(shù)為4階。

對于數(shù)字強(qiáng)震動記錄，濾波通帶的高通低頻截止頻率和低通高頻截止頻率一直是濾波器中的重要參數(shù)。根據(jù)奈奎斯特采樣定理，強(qiáng)震動記錄的高頻截止頻率最大值不得大于采樣頻率的二分之一，低頻截止頻率的最小值（低頻截止頻率的下限）不得小于記錄時間長度倒數(shù)的2倍。CSMIP一般選23Hz作為高頻濾波的截止頻率，25Hz作為終止頻率，選0.05—0.07Hz作為低頻濾波的截止頻率。Syun'itiro等（1988）根據(jù)擬速度譜（阻尼比為5%）在長周期段的特點(diǎn)來確定截止頻率。Kenzo等（1988）研究了SMAC強(qiáng)震儀所獲得的強(qiáng)震記錄，截止頻率由1/譜確定。根據(jù)謝禮立等（1983）和周寶峰（2012）的研究，高頻截止頻率對于強(qiáng)震動記錄位移時程基線偏移的影響很小，考慮到抗震結(jié)構(gòu)設(shè)計(jì)中比較感興趣的頻帶范圍，本文高頻截止頻率取35Hz。

低頻截止頻率對于強(qiáng)震動記錄位移時程基線偏移的影響很大。而目前已有的低頻截止頻率確定方法，如經(jīng)驗(yàn)法、擬速度譜法、震源譜理論法、傅里葉幅值譜法等，都不能夠有效且準(zhǔn)確地濾除低頻噪聲。比較好的方法是通過濾波后積分速度和位移的效果來判定選擇的低頻截止頻率的合理性和準(zhǔn)確性，也就是傳統(tǒng)方法，但是這種方法的效率很低，不適用于我國強(qiáng)震動臺網(wǎng)海量強(qiáng)震動記錄的處理。本文針對傳統(tǒng)方法的不足，開展自動搜索算法的研究，以滿足目前海量強(qiáng)震動記錄處理的需要，并把傳統(tǒng)方法獲得的結(jié)果作為自動搜索算法準(zhǔn)確性的參考標(biāo)準(zhǔn)。

1 基礎(chǔ)數(shù)據(jù)

本文選擇的強(qiáng)震動記錄為2008年汶川和2013年蘆山兩次地震主震和余震獲得的強(qiáng)震動加速度記錄。在這兩次地震中，國家強(qiáng)震動臺網(wǎng)獲得了大批完整的強(qiáng)震動記錄。從汶川和蘆山地震加速度記錄中篩選出符合條件的作為研究的基礎(chǔ)，篩選的條件如下：

（1）三通道加速度峰值都大于等于20cm/s2。

（2）避免選取具有奇異波形特征（周寶峰等，2014a；Zhou等，2014b）的強(qiáng)震動記錄。

由此，可以從汶川地震數(shù)據(jù)（中國地震局震害防御司，2008；于海英等，2009a）和蘆山地震數(shù)據(jù)（國家強(qiáng)震動臺網(wǎng)中心，2014；于海英等，2014）中，挑選符合條件的主余震強(qiáng)震動加速度記錄共609條。所選強(qiáng)震動記錄震級和強(qiáng)震動記錄數(shù)目的關(guān)系如圖1所示。

圖1 各震級的強(qiáng)震動數(shù)目

2 低頻截止頻率自動搜索算法

2.1 算法原理與模型

傳統(tǒng)方法通過觀察強(qiáng)震動記錄濾波后的位移時程來確定低頻截止頻率。在臺站不發(fā)生永久位移的情況下，合理濾波后的位移時程末尾段均值接近于零并且末尾段擬合直線保持水平狀態(tài)。類比傳統(tǒng)方法提出確定低頻截止頻率的自動搜索算法模型如下：

（1）設(shè)置一系列頻率點(diǎn)存儲于一維數(shù)組LPS中，并將LPS按照升序進(jìn)行排列，作為低頻截止頻率搜索的范圍。頻率點(diǎn)取的越密集，計(jì)算結(jié)果越精確，但是將增加數(shù)據(jù)的計(jì)算量。

（2）依次從LPS中取出一個低頻截止頻率進(jìn)行濾波。因?yàn)長PS中的頻率點(diǎn)已經(jīng)升序排列，所以可以盡可能保留強(qiáng)震動記錄的低頻信息。

（3）根據(jù)濾波后的加速度時程計(jì)算位移時程。

（4）判斷位移時程末尾段均值是否小于該位移時程峰值（PGD）的1/1并且位移時程末尾段斜率的絕對值小于對應(yīng)位移峰值的1/2。若不滿足條件，跳轉(zhuǎn)至（2）；若滿足條件，即確定低頻截止頻率。本文1和2為正數(shù)，根據(jù)傳統(tǒng)方法的經(jīng)驗(yàn)，1的范圍為1—5，每間隔1取5個測試點(diǎn)；2的范圍為50—800，每間隔10取1個測試點(diǎn)。

（5）將第4步確定的頻率與根據(jù)奈奎斯特采樣定理確定的低頻截止頻率的下限進(jìn)行比較，取較大值作為低頻截止頻率。

算法模型流程圖如圖2所示。

圖2 自動搜索算法模型流程圖

2.2 參數(shù)確定

在已有理論和研究的基礎(chǔ)上，先對自動搜索模型中涉及到的相關(guān)參數(shù)做出如下約定：

（1）從0.04—1Hz之間間隔0.01Hz抽取97個頻率點(diǎn)作為低頻截止頻率的搜索范圍。

（2）采用雙向巴特沃斯非因果濾波器，濾波階數(shù)為4階。

（3）濾波器的高頻截止頻率為35Hz。

（4）位移時程末尾段定義為強(qiáng)震動記錄末尾占記錄總長度四分之一的部分。

影響傳統(tǒng)方法獲取濾波低頻截止頻率的主要因素是對強(qiáng)震動記錄位移時程情況的主觀判斷，因此自動搜索模型流程的結(jié)束條件是影響自動搜索算法精度的主要因素。從算法模型流程圖（圖2）中可知，確定算法模型流程的結(jié)束條件是確定1和2兩個重要參數(shù)的值。在確定1和2之前，根據(jù)統(tǒng)計(jì)學(xué)習(xí)方法（李航，2012）中的思想，作如下說明：

（1）參與統(tǒng)計(jì)的所有強(qiáng)震動記錄存儲于X中。

（2）將自動搜索算法模型稱為決策函數(shù)fn(X)，其中n＝(1，2)T，T表示轉(zhuǎn)置。決策函數(shù)fn(X)表示利用n計(jì)算每一個強(qiáng)震動記錄的低頻截止頻率，即自動搜索算法模型。

（3）根據(jù)傳統(tǒng)方法確定強(qiáng)震動記錄的低頻截止頻率存儲于Y中。

（4）X、fn(X)和Y為列向量。

作出如上說明后，需要制定判斷n優(yōu)劣的準(zhǔn)則，這樣自然地引入損失函數(shù)來度量算法預(yù)測的好壞。損失函數(shù)是fn(X)和Y的非負(fù)實(shí)值函數(shù)，記作(Y，fn(X))。損失函數(shù)的種類比較多，本文采用的損失函數(shù)為：

損失函數(shù)L(Y，fn(X))越小，決策函數(shù)fn(X)所采用的參數(shù)n就越好。根據(jù)傳統(tǒng)方法的經(jīng)驗(yàn)，自動搜索算法模型在n1、n2的取值范圍進(jìn)行搜索測試點(diǎn)。所有測試點(diǎn)計(jì)算出的損失函數(shù)如圖3所示。其中n1＝4和n2＝440時，損失函數(shù)L(Y，fn(X))＝0.58為最小值。因此自動搜索算法的結(jié)束條件為：位移時程末尾段均值小于對應(yīng)位移峰值的1/4并且末尾段線性回歸斜率的絕對值小于對應(yīng)位移峰值的1/440。

2.3 實(shí)際震例的應(yīng)用

本文以傳統(tǒng)方法確定的低頻截止頻率為參考標(biāo)準(zhǔn)，基于汶川和蘆山兩次地震主震和余震獲得的609條強(qiáng)震動加速度記錄，利用自動搜索算法計(jì)算低頻截止頻率。傳統(tǒng)方法和自動搜索算法的結(jié)果如圖4所示。

圖4 傳統(tǒng)方法和自動搜索算法的結(jié)果

以記錄051YAL130420113402為例，依據(jù)自動搜索算法確定低頻截止頻率。該記錄的時程長度為80.625s，根據(jù)奈奎斯特采樣定理，該記錄的低頻截止頻率下限為0.025Hz。該記錄濾波后加速度時程、濾波后位移時程、PGD統(tǒng)計(jì)曲線（圖中fPGD是不同頻率濾波后的PGD，rPGD是未濾波的PGD）、位移時程末尾段均值統(tǒng)計(jì)曲線以及末尾段線性回歸斜率的統(tǒng)計(jì)曲線如圖5所示。

根據(jù)傳統(tǒng)方法，利用不同的低頻截止頻率繪制濾波后速度時程、位移時程如圖6所示。由圖6可知低頻截止頻率對速度時程的影響很小，但是，對位移時程的影響很大。通過對位移時程尾部的觀察，對于記錄051YAL130420113402，想得到合理的位移時程，低頻截止頻率應(yīng)為0.08Hz，和頻率搜索法確定的0.08Hz一致。

從圖4可以看出無論是傳統(tǒng)方法還是自動搜索算法，確定大部分強(qiáng)震動記錄的低頻截止頻率都在0.1Hz左右，自動搜索算法得出的結(jié)果與傳統(tǒng)較為接近。從圖5和圖6可知，由自動搜索算法確定低頻截止頻率的過程以及計(jì)算效率都是可行的。在實(shí)際海量強(qiáng)震動記錄處理工作中，自動搜索算法的計(jì)算效率比傳統(tǒng)方法有顯著提高，節(jié)省大量人力。原則上只需要調(diào)整本算法的結(jié)束條件就可以適用于不同地震的強(qiáng)震動記錄。

圖5 記錄051YAL130420113402頻率搜索

2.4 誤差統(tǒng)計(jì)與分析

自動搜索算法主要在張同宇（2016）關(guān)于低頻截止頻率的研究基礎(chǔ)上發(fā)展而來，以下將張同宇關(guān)于低頻截止頻率的研究稱為張同宇算法。自動搜索算法和張同宇算法主要區(qū)別是對于數(shù)據(jù)位移時程末尾段的定義和算法結(jié)束條件的不同。依據(jù)張同宇的研究，位移時程末尾段被定義為強(qiáng)震動記錄末尾10s的長度。算法的結(jié)束條件是：

（1）位移時程末尾段10s的均值小于0.1cm；

（2）位移時程末尾段10s的均值小于0.04Hz濾波后位移時程末尾段10s均值的1/10；

（3）位移時程末尾段10s的均值小于對應(yīng)位移時程峰值的1/10。

為評估自動搜索算法的有效性，利用自動搜索算法和張同宇算法對挑選出的609條強(qiáng)震動記錄進(jìn)行處理，確定每條記錄的低頻截止頻率。為比較準(zhǔn)確度，以傳統(tǒng)方法為標(biāo)準(zhǔn)，將自動搜索算法和張同宇算法分別與傳統(tǒng)方法獲得的低頻截止頻率作差，得到誤差統(tǒng)計(jì)結(jié)果。自動搜索算法和張同宇算法的誤差統(tǒng)計(jì)圖如圖7所示，誤差統(tǒng)計(jì)表如表1所示。可以看出，相對于傳統(tǒng)方法，自動搜索算法所確定的低頻截止頻率誤差基本上都在0.03Hz以內(nèi)。

圖6 記錄051YAL130420113402濾波后速度、位移時程曲線

圖7 自動確定低頻截止頻率誤差統(tǒng)計(jì)

2.5 本文算法適用條件

從圖7可知，有兩條強(qiáng)震動記錄低頻截止頻率的誤差值在0.2Hz以上，這兩條記錄分別為051AXT080512142801和051CDZ130420080203，其加速度時程如圖8所示。

對于強(qiáng)震動記錄051AXT080512142801，根據(jù)于海英等（2009a）的研究可知，本記錄的強(qiáng)震動臺站受汶川主震的影響發(fā)生了永久位移，在汶川8.0級地震中有21組近場強(qiáng)震動臺站都有不同程度的永久位移發(fā)生，由于本文篇幅限制，對以上21組三分量近場強(qiáng)震動記錄的積分計(jì)算和濾波結(jié)果不一一列出。如圖9左圖所示，利用自動搜索算法確定的低頻截止頻率對該強(qiáng)震動記錄進(jìn)行帶通濾波（低頻截止頻率0.29Hz）以后，兩次積分得到的位移時程就會損失掉永久位移信息。對大地震近場獲得的強(qiáng)震動記錄經(jīng)濾波處理會濾掉真實(shí)的地震動記錄低頻成分（于海英等，2009b）。所以，不僅該強(qiáng)震動記錄不適合濾波，凡是有永久位移的強(qiáng)震動記錄都不適合采用濾波進(jìn)行零線校正處理，即不適合采用自動搜索算法確定低頻截止頻率。

對于強(qiáng)震動記錄051CDZ130420080203，從圖8右圖可知，本記錄數(shù)據(jù)不完整。記錄前10s的最大值為0.055cm/s2，后10s的最大值為5.190cm/s2，震前和震后的地面震動差別比較大，所以判斷本記錄缺少末尾段。從圖9右圖可知，利用自動搜索算法確定的低頻截止頻率對該強(qiáng)震動記錄進(jìn)行帶通濾波（低頻截止頻率0.29Hz）以后，兩次積分得到的位移時程末尾段的均值較大。所以，不完整的強(qiáng)震動記錄也不適合采用自動搜索算法確定低頻截止頻率。

綜上分析，除了要避免選取具有奇異波形特征的強(qiáng)震動記錄以外，還要求強(qiáng)震動記錄滿足以下條件時自動搜索算法才能準(zhǔn)確確定低頻截止頻率：

（1）獲得強(qiáng)震動記錄的臺站沒有發(fā)生永久位移；

（2）強(qiáng)震動記錄要完整，沒有缺少末尾段的情況。

圖8 051AXT080512142801和051CDZ130420080203加速度時程

圖9 051AXT080512142801和051CDZ130420080203位移時程

3 結(jié)語

本文通過分析選擇低頻截止頻率應(yīng)遵循的原則，以2008年汶川地震、2013年蘆山地震獲得的強(qiáng)震動加速度記錄為例，提出一種低頻截止頻率自動搜索算法，利用統(tǒng)計(jì)學(xué)習(xí)方法中的損失函數(shù)來確定該算法流程的結(jié)束條件。該算法和傳統(tǒng)方法確定的低頻截止頻率非常接近。目前，雖然確定地震動的低頻截止頻率只能對比傳統(tǒng)方法，尚無標(biāo)準(zhǔn)可以參考，然而，本文為自動確定低頻截止頻率提供了新的量化途徑。通過本文研究得出以下結(jié)論：

（1）本文提出的低頻截止頻率自動搜索算法準(zhǔn)確性接近傳統(tǒng)方法，計(jì)算效率比傳統(tǒng)方法有很大提升，特別適合批量處理海量強(qiáng)震動記錄；

（2）該算法適用于強(qiáng)震動臺站產(chǎn)出記錄完整且無永久位移的強(qiáng)震動記錄；

（3）原則上只需要調(diào)整本算法的結(jié)束條件就可以適用于不同地震的強(qiáng)震動記錄。

國家強(qiáng)震動臺網(wǎng)中心，2014．中國強(qiáng)震記錄匯報，第十七集，第一卷，蘆山7.0級地震及余震未校正加速度記錄．北京：地震出版社．

李航，2012．統(tǒng)計(jì)學(xué)習(xí)方法．北京：清華大學(xué)出版社．

謝禮立，李沙白，錢渠炕等，1983．我國強(qiáng)震記錄處理和分析方法的若干特點(diǎn)．地震工程與工程振動，3（1）：1—14．

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

于海英，王棟，楊永強(qiáng)，等，2009b．汶川8.0級地震強(qiáng)震動加速度記錄的初步分析．地震工程與工程振動，29（1）：1—13．

于海英，周寶峰，楊程等，2014．蘆山7.0級地震及余震強(qiáng)震動記錄初步分析．見：第九屆全國地震工程學(xué)術(shù)會議論文集．哈爾濱：中國建筑學(xué)會，中國地震學(xué)會，中國地震工程聯(lián)合會，153—160．

張軍華，呂寧，田連玉等，2006．地震資料去噪方法技術(shù)綜合評述．地球物理學(xué)進(jìn)展，21（2）：546—553．

張同宇，2016．基于蘆山地震強(qiáng)震動記錄的數(shù)據(jù)處理技術(shù)研究．哈爾濱：中國地震局工程力學(xué)研究所．

中國地震局震害防御司，2008．中國強(qiáng)震記錄匯報，第十二集，第一卷，汶川8.0級地震未校正加速度記錄．北京：地震出版社．

周寶峰，2012．強(qiáng)震觀測中的關(guān)鍵技術(shù)研究．哈爾濱：中國地震局工程力學(xué)研究所．

周寶峰，宋廷蘇，于海英等，2014a．蘆山強(qiáng)震記錄中的奇異波形研究．地震工程與工程振動，34（S）：93—99．

周寶峰，溫瑞智，謝禮立，2014b．強(qiáng)震記錄中的“尖刺”現(xiàn)象初步研究．土木工程學(xué)報，47（S）：295—299．

Syun'itiro O., Tokiharu O., Shigeto H., et al., 1988. Data processing method for acceleration records and its application results. In: Proceedings of the Second Workshop on Processing of Seismic Strong Motion Records. Tokyo, 119—135.

Kenzo T., Sumio S., 1988. A method for correction of accelerograms. In: Proceedings of the Second Workshop on Processing of Seismic Strong Motion Records. Tokyo, 183—199.

Zhou B. F., Wang H. Y., Xie L. L., et al., 2015. Bizarre waveforms in strong motion records. Shock and Vibration, 2015: 630362.

Automatic Search Algorithm of Low Cut-off Frequency for Filtering Strong Motion Records

Yu Haiying1, 2）, Xu Xuan1, 2）and Zhang Tongyu1, 2）

1） Institute of Engineering Mechanics, China Earthquake Administration, Harbin 150080, China 2） Key Laboratory of Earthquake Engineering and Engineering Vibration of China Earthquake Administration, Harbin 150080, China

Aiming at the low efficiency of traditional methods to determine the low cut-off frequency for filtering strong-motion records, this paper proposes an automatic search algorithm to determine the cut-off frequency of low-frequency filters and uses the loss function in a statistical learning method to determine the end condition of the automatic search model flow. Based on strong-motion records obtained from the main shock and aftershocks of the 2008 Wenchuan earthquake and 2013 Lushan earthquake, we compare the results of automatic search algorithms and the low cut-off frequency identified by traditional methods, analyze the errors of automatic search algorithms, and propose principles and applicable conditions for their use. The results show that the algorithm improves computational efficiency significantly compared to traditional methods; in particular, the algorithm is suitable for the batch-processing of massive strong-motion records.

Strong motion records; Filter; Low cut-off frequency; Automatic search algorithm; Loss function

于海英，徐旋，張同宇，2018．強(qiáng)震動記錄的低頻截止頻率自動搜索算法．震災(zāi)防御技術(shù)，13（1）：65—74．

10.11899/zzfy20180106

國家重點(diǎn)研發(fā)計(jì)劃（2017YFC1500803），國家自然科學(xué)基金高鐵聯(lián)合基金資助項(xiàng)目（U1534202），黑龍江省自然科學(xué)基金（E2015070），國家自然科學(xué)基金青年項(xiàng)目（51308517）

2017-11-20

于海英，男，生于1962年。研究員。主要從事強(qiáng)震動觀測與數(shù)據(jù)處理技術(shù)方面研究。E-mail：haiyingyu@126.com

徐旋，男，生于1992年。碩士研究生。主要從事強(qiáng)震動數(shù)據(jù)處理技術(shù)方面研究。E-mail：xu_dexuan@outlook.com