基于STA/LTA算法的地脈動有效信號自動識別研究

朱 勝，李 平，李玉影，劉應(yīng)慈，周 楷

（1.防災(zāi)科技學(xué)院地質(zhì)工程學(xué)院，河北三河 065201；2.河北省地震災(zāi)害防御與風(fēng)險評價重點實驗室，河北三河 065201）

引言

地脈動是地面上一種平穩(wěn)的非重復(fù)性隨機(jī)波動，主要是由風(fēng)、海浪活動等自然振源或場地周圍的機(jī)械振動、交通工具等人工振源激發(fā)，經(jīng)過在不同性質(zhì)的巖土層界多次反射和折射，疊加后傳播到地面的振動。地脈動具有頻率低、振幅小等特點，其頻率范圍一般為0.1～10 Hz，而振幅一般在0.1～1μm。地脈動的不同振幅和頻率的變化，既能反映測試場地的土層特征，又能反映工程場地的動力特性。早在1907年，學(xué)者Omori［1］觀測到了一種振幅很小的振動，他將這種微小的振動稱之為地動脈（Microtremor），并發(fā)表了《論地脈動》一文，在文中初步描述了地脈動的振動特征，此后日本和美國的學(xué)者對地脈動這一熱點課題進(jìn)行了大量的研究。但是很多學(xué)者對地脈動產(chǎn)生的機(jī)理不明確，所以地脈動的研究一度陷入沉寂，直到19世紀(jì)60年代日本學(xué)者Kiyoshi等［2］提出將地脈動測試應(yīng)用于場地評價中去，日本學(xué)者Nakamura提出用單點譜比法處理地脈動數(shù)據(jù)，地脈動研究又重新引起國內(nèi)外學(xué)者的興趣，隨后美國、歐盟、中國等各國學(xué)者也相繼開展大量的地脈動測試研究，促進(jìn)了地脈動的在工程中的應(yīng)用和發(fā)展；但是Lermo等［3］一直認(rèn)為地脈動產(chǎn)生的機(jī)理不明確，對地脈動用于場地評價仍一直持有否定的態(tài)度；地脈動的研究至今已有百年的歷史，我國學(xué)者對地脈動的研究起步于20世紀(jì)60年代，郭明珠等［4］，師黎靜等［5］，劉紅彪等［6］，分別地脈動測試應(yīng)用有在震害分析、建筑結(jié)構(gòu)設(shè)計、場地類別劃分和地基場地處理效果評價等工作中取得了一些的工程應(yīng)用成果；師黎靜等［7］通過選取46處場地的地脈動數(shù)據(jù)，采用單點譜比法計算出場地卓越頻率，探究地脈動測試法應(yīng)用于場地特征參數(shù)測試的適用性，并且給出了覆蓋層厚度、不同深度等效剪切波速以及特征周期與卓越頻率之間的擬合曲線，驗證其相關(guān)性；王偉君等［8］利用在北京五棵松地區(qū)進(jìn)行地脈動測試，采用單點譜比法分析場地卓越頻率，臺陣地脈動測試F-K法，反演土層的速度結(jié)構(gòu)。研究成果表明：地脈動能夠反映出場地土層的構(gòu)造信息，相比于目前的場地動力特性測試方法，具有經(jīng)濟(jì)、便捷和高效的特點，有廣闊的工程應(yīng)用前景。就目前而言，學(xué)者們已經(jīng)對地脈動測試的現(xiàn)場采集方法達(dá)成了共識，但是在后期數(shù)據(jù)的處理仍存在爭議。鑒于此，文中以松原5.7級地震砂土液化場地的地脈動測試數(shù)據(jù)為研究對象，施加反STA/LTA算法，探究數(shù)據(jù)處理中有效信號時間窗的選取問題，為以后規(guī)范化處理數(shù)據(jù)提供依據(jù)。

1 地脈動數(shù)據(jù)來源

2018年5月28日1點50分，在吉林松原市寧江區(qū)發(fā)生了5.7級地震，震源深度為13 km。地震發(fā)生后，防災(zāi)科技學(xué)院迅速組織人員，對震中受災(zāi)區(qū)域開展震害調(diào)查工作。選取了26處典型的場地進(jìn)行地脈動測試數(shù)據(jù)采集，其中液化場地為19處，未發(fā)生液化場地7處，先后2年時間一共采集108條地脈動數(shù)據(jù)記錄。

本次測試設(shè)備采用美國Kinemet-rics公司生產(chǎn)的Basalt型一體式強(qiáng)震儀，其內(nèi)部設(shè)置三通道加速度計，能夠同時記錄2個水平方向和垂直方向的加速度，儀器的動態(tài)范圍為155 dB。儀器布置時利用地質(zhì)羅盤定位，按南北、東西和垂直3個正交方向放置。采樣率設(shè)為200 sps，電壓靈敏度為2.5 V/g，每個測點每次測試連續(xù)采樣300 s，每個測點數(shù)據(jù)采樣3次［19］，地脈動的采集結(jié)果如圖1所示。

圖1 測點M1地脈動三分量記錄結(jié)果Fig.1 Microtremors 3-components time-histories at site M1

地脈動測試是一種在無強(qiáng)震源的情況下，僅僅是依靠自然震源或人為震源產(chǎn)生的地表微弱的振動，需要使用靈敏度較高的加速度傳感器進(jìn)行監(jiān)測得到。因為振動幅值比較小，所以在數(shù)據(jù)采集的過程中，需在觀測環(huán)境中在一定范圍內(nèi)無特定的振動源，如汽車、機(jī)械運(yùn)轉(zhuǎn)的影響等。地脈動的測點一般選擇自由場地，并且遠(yuǎn)離道路。選擇的數(shù)據(jù)采集的時間一般選擇晚上，但是考慮到實際的測試情況進(jìn)而可以考慮選擇白天進(jìn)行地脈動測試，從早上8點到下午6點，進(jìn)行地脈動數(shù)據(jù)觀測。數(shù)據(jù)采集時一般在風(fēng)和日麗的天氣進(jìn)行。刮風(fēng)下雨天氣對測試結(jié)果也有一定的影響，下雨之后，雨水滲透作用導(dǎo)致覆蓋土層抗剪強(qiáng)度下降，對高頻波有抑制作用，對低頻波有放大效應(yīng)。

本次地脈動的測試采集應(yīng)為是地震應(yīng)急科考期間時間比較急迫，所以選擇白天開展地脈動測試。在地脈動數(shù)據(jù)處理的過程中，發(fā)現(xiàn)了地脈動數(shù)據(jù)有一定的“毛刺”和脈沖信號干擾，采用帶通濾波仍然無法有效的去除，如圖2所示。因此，文中提出對地脈動信號施加反STA/LTA算法，提取地脈動的有效信號段，剔除“毛刺”和脈沖信號段以達(dá)到濾波的效果。

圖2 測點M 13地脈動三分量記錄結(jié)果Fig.2 Microtremors 3-components time-histories at site M13

2 STA/LTA法簡介

震相識別是地震學(xué)研究中一個重要課題分支，它為地球內(nèi)部構(gòu)造的研究、振源的定位、震源機(jī)制的推演、地震速報和預(yù)警等一系列研究工作提供了基礎(chǔ)性的研究資料。根據(jù)識別的方法，可將震相識別劃分為人工識別和自動識別。傳統(tǒng)的人工識別方法主要有：直觀檢驗法、走時表檢驗法、綜合檢驗法和合成地震圖檢驗法等等，這些方法都是基于識別者的工作經(jīng)驗來加以判斷、辨識，主要的缺點在于判別效率低，并且由于人為的誤差、讀數(shù)錯誤，可能導(dǎo)致不同的人識別的精度相差較大［9-10］。隨著2008年以后我國數(shù)字化地震臺網(wǎng)的大規(guī)模建立，地震觀測數(shù)據(jù)量巨大，單單依靠人工識別是不切實際的，因此應(yīng)運(yùn)而生了穩(wěn)定、快速和高效的自動識別方法。目前，自動識別方法主要分為：時域分析法、頻域分析法、時頻分析法、模式識別法和綜合分析法［9-10］。文中著重介紹時域分析法中的STA/LTA法。

STA/LTA法是由Stevenson［11］最早提出，后經(jīng)過Allen等［12］的不斷改進(jìn)，用來檢測地震事件和拾取震相。其原理是在原信號上施加兩個時間窗，分別為短時窗和長時窗，用短時窗的平均值（short-term average，簡稱STA）和長時窗的時長（long-term average，簡稱LTA）的比值（R）來反映信號能量的變化程度，當(dāng)?shù)卣鹦盘柕竭_(dá)時，STA值比LTA值變化的快，STA值快速變大，R值會有一個顯著的增大趨勢，當(dāng)R值超過設(shè)定的閾值（THR）時，則此點被判斷為初動，地震記錄系統(tǒng)開始記錄有效信號。

式中：X(i),(i=1,2,3···,N)，表示短時間窗內(nèi)記錄數(shù)據(jù)；Y(j),(j=1,2,3···,M)，表示長時窗內(nèi)記錄數(shù)據(jù)，M和N分別表示長短窗內(nèi)的樣本個數(shù)。

STA/LTA方法計算量小、運(yùn)算程序簡單，但是在實際應(yīng)用過程中，其計算結(jié)果常常出現(xiàn)計算偏差，導(dǎo)致誤判的現(xiàn)象。為了提高信號識別的準(zhǔn)確率，增加信號到達(dá)時的靈敏度，所以學(xué)者Allen［10］提出基于原始的地震記錄，重新構(gòu)造能夠靈敏反應(yīng)信號時間序列CF（t），并將新的時間序列作為計算的輸入?yún)?shù)。新的時間序列C F（t）被稱為特征函數(shù)。常用的特征函數(shù)有3種：

隨著研究的不斷深入，很多學(xué)者提出改進(jìn)的特征函數(shù)形式，國內(nèi)學(xué)者武東坡和高淑芳［9，17］利用三角函數(shù)推導(dǎo)，提出了一種改進(jìn)的STA/LTA震相自動識別的特征函數(shù)：

劉晗通過大量微震監(jiān)測信號分析發(fā)現(xiàn)，STA/LTA算法特征函數(shù)取式（3）和式（5），對微震記錄信號識別響應(yīng)比較靈敏，并且式（5）對頻率變化也比較敏感［13］。根據(jù)短時窗與長時窗的關(guān)系，可以將STA/LTA計算方法分為標(biāo)準(zhǔn)STA/LTA法和延遲STA/LTA法。標(biāo)準(zhǔn)STA/LTA法是短時窗與長時窗之間沒有重疊或延遲，進(jìn)行截取計算；而延遲STA/LTA法，顧名思義即在短時窗和長時窗之間有一段延遲。實際應(yīng)用中，為了提高計算效率，在標(biāo)準(zhǔn)法的基礎(chǔ)上，逐漸發(fā)展了遞歸STA/LTA法。遞歸法相比于標(biāo)準(zhǔn)法計算更加快速，拾取信號的靈敏度更高，對脈沖信號有一定的抑制作用，減小誤觸的風(fēng)險。

遞歸法的計算公式如下：

式中：STAi和LTAi分別為信號在i時刻短時窗與長時窗的平均值；C F（i）為信號在i時刻特征函數(shù)值；Nsta和Nlta分別為信號短時窗和長時窗的信號記錄點數(shù)。

3 基于地脈動數(shù)據(jù)有效信號自動識別

目前隨著各個學(xué)科之間不斷交融與發(fā)展，STA/LTA算法的應(yīng)用不僅僅局限于震相識別領(lǐng)域的，在其他信號“初動”識別領(lǐng)域也被廣泛應(yīng)用，如煤礦沖擊地壓的預(yù)測、煤與瓦斯突出的預(yù)防、深部礦山評估巖體穩(wěn)定性、預(yù)測巖爆等動力型破壞等方面［13-14］。

地脈動是一種平穩(wěn)、隨機(jī)波動，其最大的特點就是振幅小、頻率低。地脈動測試受環(huán)境影響較大，因此測試最理想的時間段應(yīng)選擇在夜深人靜的深夜，但是在實際測量中卻做不到這一點。而白天地脈動測試時，拾震器往往接收到一些振幅較大的干擾信號，干擾信號直接影響數(shù)據(jù)處理的結(jié)果，對結(jié)果造成不利的影響，因此需要對這些干擾信號進(jìn)行篩選剔除。文中基于震相自動識別思想，對地脈動數(shù)據(jù)施加反STA/LTA算法，快速有效地對地脈動數(shù)據(jù)進(jìn)行篩選，自動識別出異常信號，獲得有價值的信號，為后期數(shù)據(jù)處理奠定了基礎(chǔ)。

3.1 反STA/LTA算法識別流程

地脈動數(shù)據(jù)施加反STA/LTA算法的識別過程為：信號輸入→計算→判斷→拾取。輸入?yún)?shù)有地脈動測試信號Y（t）、數(shù)據(jù)采樣頻率fs、短時窗的時長STA、長時窗的時長LTA、對應(yīng)的判斷閾值Rd和切片時間窗時長Tw。識別算法具體流程如下：

（1）計算地脈動數(shù)據(jù)t時刻的C F（t）值；

（2）計算STA/LTA法的t時刻的R（t）值；

（3）判別t時刻的R（t）值是否小于Rd值，若判斷結(jié)果為是，則進(jìn)行下一步判斷；若判斷結(jié)果為否，則進(jìn)入下一時刻并返回（2）；

（4）將地脈動時程進(jìn)行切片為t時刻到t+Tw時刻，分別計算出切片段的R（t）值，判別切片段的最大值是否小于Rd值，若判斷結(jié)果為是，則進(jìn)行下一步；若判斷結(jié)果為否，則進(jìn)入下一時刻并返回（2）；

（5）記錄下T0=t時刻和T1=t+Tw時刻值，生成時間窗矩陣，并進(jìn)行下一步判別；

（6）判別t時刻是否到達(dá)結(jié)束時刻t0，若判斷結(jié)果為是，則繪制時間窗圖形；若判斷結(jié)果為否，則將進(jìn)入t+Tw時刻，并返回（2）。

（7）如圖3所示，為算法程序的流程圖。

圖3 STA/LTA算法流程圖Fig.3 Flow chart of STA/LTA algorithm

3.2 反STA/LTA算法的影響因素探究

由于STA/LTA算法常常被用于地震震相的識別，震相識別與地脈動干擾信號識別的相似之處在于：兩者都是在平穩(wěn)信號中，識別突變的異常信號；而地震信號與地脈動信號的最大的區(qū)別在于：當(dāng)?shù)卣鹦盘杹砼R時，振幅迅速增大，短時窗平均值會迅速變化，STA/LTA值短時間內(nèi)增大，超過閾值，地震記錄系統(tǒng)開始工作，記錄地震信號；而地脈動最大的特點是振幅小、頻率低，并且干擾信號的幅值，相較于地脈動信號的幅值相差不大。文中將從STA/LTA計算方法、特征函數(shù)和長短時窗的時長的選取等方面，討論STA/LTA法對于地脈動測試數(shù)據(jù)處理的適用性。文中選取了一條在2018年5月27日松原發(fā)生5.7級地震后采集的典型地脈動測試記錄作為研究對象。此地脈動記錄點為姜家圍子村委會前測試點，截取的時程為5 min，如圖4所示。

圖4 測點M1地脈動三分量記錄結(jié)果Fig.4 Microtremors 3-components time-histories at site M1

3.2.1 不同的計算方法對算法的影響

如圖5（a）和圖5（b）所示，STA/LTA的算法分別選擇的是遞歸法和標(biāo)準(zhǔn)法，遞歸法選擇的是式（6）、式（7），特征函數(shù)選擇的是式（5），短時窗的時長為4 s，長時窗的時長為20.48 s，閾值設(shè)定為3，時間窗長為20.48 s。在數(shù)據(jù)處理的過程中發(fā)現(xiàn)，地脈動的三分量中的垂直分量受脈沖影響較大，在施加STA/LTA算法的計算過程中，以垂直分量作為計算的首選計算樣本，截取出數(shù)據(jù)窗之后再檢驗水平分量的R值是否超過閾值，若超過閾值則將該數(shù)據(jù)窗舍去，反之沒有超過則進(jìn)行下一步的計算。為了減小之后對地脈動數(shù)據(jù)分析處理的誤差，相鄰的時間窗設(shè)定重合50%，對比2種方法，結(jié)果表明：采用遞歸法，對干擾信號的識別的靈敏度較高，而傳統(tǒng)的標(biāo)準(zhǔn)法，對干擾信號計算R值有一定的平滑作用，靈敏度較低。

如圖5（a）、（c）和（d）所示，探究特征函數(shù)的選取對STA/LTA算法的影響。特征函數(shù)分別選擇的是式（5）、式（2）和式（3），STA/LTA計算方法選擇遞歸法，其他參數(shù)都設(shè)置一樣的，短時窗的時長為4 s，長時窗的時長為20.48 s，閾值設(shè)定為3，時間窗長為20.48 s，相鄰的時間窗重合50%；對比3種不同特征函數(shù)，如圖5中的相同段e與f段，結(jié)果表明：特征函數(shù)選擇平方法和高淑芳提出的改進(jìn)法，對于干擾信號識別的靈敏度較高，而采用絕對值法的對干擾信號的識別較差。其中原因可能是地脈動的幅值較小，利用絕對值法計算出的特征函數(shù)，對R值的計算結(jié)果的放大作用不明顯。

圖5 不同特征函數(shù)選取對時間窗的影響Fig.5 The influence of different characteristic function selection on time window

3.2.2 不同參數(shù)選取對算法的影響

如圖5（e）～（h）所示，探究短時窗和長時窗的時長選取值對STA/LTA算法結(jié)果的影響。其中計算方法選擇遞歸法，特征函數(shù)選擇的是式（5），短時窗的時長分別選取了2、4、8 s，長時窗的時長為10.24、20.48、20.72 s，閾值仍然設(shè)定為3，時間窗時長為20.48 s，相鄰的時間窗重合50%。經(jīng)過對比，如圖6中的g與g′段，結(jié)果表明：隨著短時窗的時長的增加，STA/LTA算法對干擾識別的靈敏度是降低，而長時窗的時長變化對識別信號的靈敏度影響較小，尤其是持續(xù)信號干擾段，算法的識別能力較差，可以通過減小短時窗時長或適當(dāng)增大長時窗時長的方法提高精度。

圖6 不同參數(shù)選取對的時間窗選擇Fig.6 The influence of different parameter selection on time window selection

文中沒有討論反STA/LTA算法的閾值（THR）設(shè)定的問題，本篇文章中提及的所有對比試驗閾值都取值3。其主要原因是：在地脈動分析數(shù)據(jù)時施加反STA/LTA算法的主要目的是對時程曲線進(jìn)行快速切片提取、過濾干擾信號，以便于下一步的數(shù)據(jù)分析處理。在實際分析地脈動數(shù)據(jù)的工作時，由于地脈動是幅值比較小的是平穩(wěn)信號，而干擾信號多為脈沖信號，易于判斷出干擾信號，從而可根據(jù)時程曲線圖來初步判斷出閾值，并調(diào)整閾值的取值。

利用本文的結(jié)論，對地脈動采集數(shù)據(jù)施加反STA/LTA算法快速提取出有效信號段，進(jìn)而利用單點譜比法求取場地的卓越頻率。單點譜比法也稱為H/V譜比法，是日本學(xué)者Nakamura［18］提出的一種場地周期計算的方法，在文獻(xiàn)［17］中對單點譜比法的計算也有探討，文中不再贅述。將提取出來的時間窗分別計算H/V譜，并計算其平均值和標(biāo)準(zhǔn)差，最終求得卓越頻率。由于篇幅有限文中只給出了測點M1的結(jié)果，如圖7所示。由于干擾信號都為脈沖信號，低頻成分豐富，相比于施加反STA/LTA算法之后的結(jié)果，卓越頻率略有減小。

圖7 H/V譜比法結(jié)果對比Fig.7 Comparison of results of H/V spectra

對松原地區(qū)采集的108條地脈動數(shù)據(jù)記錄施加反STA/LTA算法，參數(shù)設(shè)置為：STA分別為1 s和2 s，LTA為20.48 s，R分別設(shè)置為1.5和2.5。如圖8所示，第1組（STA=1 s，LTA=20.48 s，R=1.5）的靈敏度較高，提取有效時窗個數(shù)平均值在10個；第2組（STA=1 s，LTA=20.48 s，R=2.5），提取的有效時窗個數(shù)平均值在15個；第3組（STA=2 s，LTA=20.48 s，R=2.5），提取的有效時窗個數(shù)平均值在20個。經(jīng)過對比分析，結(jié)果表明：（1）當(dāng)短時窗時長較小，施加反STA/LTA算法敏靈度較高，提取時窗個數(shù)較少；（2）當(dāng)提取有效時窗個數(shù)小于10個時，計算得到的場地卓越頻率有減小的趨勢；（3）也存在當(dāng)提取有效時窗個數(shù)達(dá)到15個時，計算得到的場地卓越頻率呈現(xiàn)穩(wěn)定的趨勢。

圖8 H/V譜比法結(jié)果對比Fig.8 Comparison of results of H/V spectra

4 結(jié)論與討論

文中以地脈動測試數(shù)據(jù)為研究對象，施加反STA/LTA算法，提取地脈動信號的平穩(wěn)段。通過改變計算方法、不同的特征函數(shù)和長短時窗的時長的選取，探究其對于地脈動測試數(shù)據(jù)分析的可行性。對108條地脈動數(shù)據(jù)記錄分析，得到了如下結(jié)論：（1）反STA/LTA算法程序?qū)Φ孛}動測試數(shù)據(jù)的干擾信號，能夠成功的識別；（2）使用遞歸法STA/LTA，特征函數(shù)選用平方法和高淑芳提出的改進(jìn)法靈敏度較高；（3）短時窗的時長對結(jié)果影響較大，而長時窗的時長取值對結(jié)果影響較小；（4）利用算法提取的時窗數(shù)量越多，計算求得的卓越頻率越趨于穩(wěn)定。因此，推薦使用遞歸STA/LTA法參數(shù)設(shè)置為：短時窗（STA）的時長設(shè)定為1～2 s，長時窗（LTA）的時長設(shè)定為20.48 s，閾值（THR）設(shè)定為2.5，對地脈動測試數(shù)據(jù)進(jìn)行處理。

在地脈動數(shù)據(jù)的處理時，需要人工截取數(shù)據(jù)段進(jìn)行分析，存在偶然性。文中主要的探討問題是：如何規(guī)范化處理地脈動數(shù)據(jù)。進(jìn)而文中論證了將STA/LTA算法應(yīng)用于地脈動數(shù)據(jù)處理的可行性，還要在今后的實際工程應(yīng)用中驗證其適用性。對于干擾信號的有效識別取決于算法的靈敏度和閾值的設(shè)定值。特征函數(shù)的選取也間接的影響算法的靈敏度，所以對特征函數(shù)的優(yōu)化選取將是下一步工作的重點；同時在數(shù)據(jù)處理過程中，當(dāng)采集的地脈動數(shù)據(jù)如圖9所示時，測試過程中一直有脈沖荷載干擾時，算法無法有效識別出有效信號，只能通過增加短時窗時長或增大閾值。

圖9 測點M13計算值Fig.9 Calculated value of M13