辛 平， 李志華

(江南大學(xué) 物聯(lián)網(wǎng)工程學(xué)院 計(jì)算機(jī)科學(xué)系，江蘇 無錫 214122)

0 引 言

在地震震相檢測(cè)中，首先要做的是精確拾取震相到時(shí)，即精確拾取P波到時(shí)。地震波的準(zhǔn)確拾取在許多與地震相關(guān)的研究中起著關(guān)鍵作用。目前，自動(dòng)拾取方法已占據(jù)主導(dǎo)地位[1]。文獻(xiàn)[2]用簡(jiǎn)單的振幅和能量閾值來確定P波到時(shí);由于該算法只用了振幅和能量是否大于閾值判定P波到時(shí)，所以其拾取的精確度不高；文獻(xiàn)[3]采用長(zhǎng)短時(shí)平均(short term average/long term average,STA/LTA)方法確定地震震相到時(shí),可以快速拾取震相到時(shí)，但其拾取震相到時(shí)的精確度受噪聲影響比較大，高噪音時(shí)出現(xiàn)誤拾甚至拾取不到P波。此外經(jīng)驗(yàn)閾值的選取直接影響STA/LTA拾取的準(zhǔn)確度；文獻(xiàn)[4]采用最大似然(maximum likelihood,ML)估計(jì)方法來確定地震P波到時(shí)，雖然能夠比較精確確定P波到時(shí)，但ML方法認(rèn)為噪聲干擾是一成不變的，忽略了突發(fā)人為或者其他非地震因素所引起的噪聲，這顯然是不合理的；文獻(xiàn)[5]通過計(jì)算P波偏振度來確定P波到時(shí),該方法提高了識(shí)別可靠度，但計(jì)算量大；文獻(xiàn)[6]將STA/LTA、偏振分析、AR分析以及瞬時(shí)頻率分析這4種方法結(jié)合起來確定P波到時(shí)，雖然精確度上有了極大的提高，但算法的計(jì)算量很大；文獻(xiàn)[7]采用自回歸方法來拾取P波到時(shí)，能夠比較精確確定P波到時(shí)；文獻(xiàn)[8]提出了一種基于局部極大值分布的P波到時(shí)自動(dòng)識(shí)別方法，算法將能量和頻率特征結(jié)合起來，能夠在高噪音時(shí)比較精確拾取P波。

本文提取地震信號(hào)的15維特征作為數(shù)據(jù)樣本，并通過特征選擇來對(duì)數(shù)據(jù)樣本進(jìn)行降維處理，最后利用支持向量機(jī)(support vector machine,SVM)進(jìn)行分類預(yù)測(cè)，確定P波到時(shí)。

1 基于改進(jìn)的局部線性嵌入和核分類方法的P波拾取

1.1 多維特征的構(gòu)造

為了獲取豐富的地震信號(hào)特征，本文通過對(duì)實(shí)測(cè)地震信號(hào)進(jìn)行變換和處理，從時(shí)域和頻域兩方面綜合選取15維特征利于數(shù)據(jù)深度挖掘。15維特征包括震相、振幅、偏振度、分形維數(shù)、STA/LTA，ML、振幅比、能量比、曲線長(zhǎng)度比、峭度、偏斜度、均方根振幅、平均能量、振幅峰態(tài)，頻域參數(shù)中選取頻率。具體如下:

1)震相到時(shí)由STA/LTA確定。

2)振幅從信號(hào)中直接讀取。

3)偏振度F采用文獻(xiàn)[5]中的方法計(jì)算，F(xiàn)=m1/m2,m1和m2分別為由地震信號(hào)數(shù)據(jù)形成的3×3階矩陣的最大特征值和次最大特征值。其矩陣構(gòu)造為

(1)

式中x，y，z為地震信號(hào)3分量數(shù)據(jù)，cov(x,y)為x和y的協(xié)方差

(2)

式中μx和μy為x和y的平均值。

4)分形維數(shù)D的計(jì)算為D=1-S,S為不同步長(zhǎng)和其對(duì)應(yīng)的曲線長(zhǎng)度擬合的直線斜率。

5)i時(shí)刻振幅比為

(3)

式中x為地震數(shù)據(jù)，m為地震數(shù)據(jù)時(shí)窗的起點(diǎn)，n為地震數(shù)據(jù)時(shí)窗的終點(diǎn)。

6)i時(shí)刻能量比為

(4)

7)i時(shí)刻曲線長(zhǎng)度比為

(5)

式中 Δt為地震數(shù)據(jù)的采樣間隔。

8)峭度為

(6)

式中s為時(shí)窗大小。

1.1.1 改進(jìn)的STA/LTA方法

傳統(tǒng)STA/LTA計(jì)算為

(7)

式中m為STA窗口大小,n為L(zhǎng)TA窗口大小。

為提高震相拾取的準(zhǔn)確度，對(duì)STA/LTA方法進(jìn)行改進(jìn),即

(8)

1.1.2 改進(jìn)的ML方法

由于震相初至前后地震信號(hào)的振幅有很大的變化，因此可以用ML估計(jì)來確定震相的初至?xí)r間。由于地震信號(hào)呈隨機(jī)分布，所以本文假設(shè)地震信號(hào)變化前后的分布都服從正態(tài)分布，則ML計(jì)算為

(9)

式中σ1為震相初至前噪音方差，σ2為震相初至后方差。

(10)

改進(jìn)的ML計(jì)算方法概括如下：

3)采用式(10)計(jì)算窗口內(nèi)的似然估計(jì)值L。

1.2 基于改進(jìn)的LLE特征抽取方法

LLE計(jì)算序列間距離采用歐氏距離，考慮每個(gè)區(qū)域上各個(gè)點(diǎn)分布是不均勻的，為了使每個(gè)區(qū)域上各點(diǎn)的分布整體呈均勻化，LLE中距離的計(jì)算改進(jìn)為

(11)

式中M(i)，M(j)分別為與其點(diǎn)距離的均值。

改進(jìn)的LLE算法具體描述如下:

1)計(jì)算每個(gè)樣本點(diǎn)xi的k個(gè)近鄰點(diǎn)xij(j=1,2,…,k)。

(12)

3)計(jì)算xi和xj的低維映射值yi和yj，使得ε(Y)最小

(13)

式中ε(Y)為損失函數(shù)，式(14)滿足

(14)

M=(I-W)T(I-W)

(15)

1.3 面向時(shí)間序列的核分類方法

由于傳統(tǒng)的SVM[10]距離度量采用歐氏距離，計(jì)算時(shí)間復(fù)雜度低(O(n2))，但其無量綱且只適用于等長(zhǎng)的時(shí)間序列[11]。由于所研究的地震信號(hào)是時(shí)間序列，具有多維、離散的特性，且每一維都有其特定的量綱，為了使SVM更好地適用于地震時(shí)間序列數(shù)據(jù)，本文提出面向時(shí)間序列的SVM(time series SVM,TS-SVM)分類方法，在TS-SVM方法中，時(shí)間序列間相似性度量采用動(dòng)態(tài)時(shí)間彎曲距離[12]。描述如下:

假設(shè)有時(shí)間序列X和Y，則其相似性度量測(cè)度為

(16)

(17)

式中 < >為空序列,Rest(X)={x2,x3,…，xn},Rest(Y)={y2,y3,…,yn}。

TS-SVM方法在核函數(shù)選取上采用高斯核函數(shù)，用DTW距離代替歐氏距離，核函數(shù)定義為

(18)

得到最優(yōu)分類函數(shù)為

f(x)=sgn{(w*,x)+b*}

(19)

新提出的TS-SVM方法描述如下：

1)由計(jì)算時(shí)間序列X的動(dòng)態(tài)彎曲距離DTW。

2)計(jì)算時(shí)間序列的內(nèi)核函數(shù)，求解二次規(guī)劃問題

得到最優(yōu)解a*={a1,a2,…，ak}。

4)計(jì)算f(x)的值，實(shí)現(xiàn)對(duì)時(shí)間序列X的分類。

動(dòng)態(tài)彎曲距離的時(shí)間復(fù)雜度為O(n2)，歐氏距離的時(shí)間復(fù)雜度也為O(n2)，TS-SVM并未改變傳統(tǒng)SVM的時(shí)間復(fù)雜度。

1.4 基于TS-SVM的 P波拾取方法AMPAT

AMPAT流程如圖1所示。

圖1 基于TS-SVM的P波拾取方法流程

AMPAT方法概括如下：

1)數(shù)據(jù)預(yù)處理：計(jì)算15維特征構(gòu)造特征空間集，考慮到特征空間的每一維特征數(shù)量級(jí)不同，需要進(jìn)行歸一化處理。對(duì)歸一化處理的特征空間用LLE進(jìn)行降維，選取最佳特征降維后的維數(shù)。選取合理數(shù)量的樣本數(shù)據(jù)，50 %樣本數(shù)據(jù)作為訓(xùn)練集，50 %作為測(cè)試集。

2)TS-SVM參數(shù)尋優(yōu):在訓(xùn)練樣本集上進(jìn)行有監(jiān)督學(xué)習(xí)的訓(xùn)練，確定最佳參數(shù)。

3)構(gòu)造P波拾取模型：對(duì)LLE降維后的特征數(shù)據(jù)，采用合理的TS-SVM尋優(yōu)參數(shù)進(jìn)行測(cè)試，拾取P波到時(shí)。

2 實(shí)驗(yàn)與結(jié)果分析

2.1 實(shí)驗(yàn)數(shù)據(jù)

本文數(shù)據(jù)由IRIS(Incorporated Research Institutions for Seismology)提供，數(shù)據(jù)格式為SAC。以美國(guó)東部海岸 2016年 3月份地震數(shù)據(jù)為樣本，利用MATLAB進(jìn)行仿真實(shí)驗(yàn)來驗(yàn)證本文提出的方法的準(zhǔn)確性。樣本的采樣頻率均為40Hz，記錄了地震信號(hào)在某一時(shí)間段內(nèi)的加速度變化情況。

2.2 結(jié)果分析

根據(jù)專家經(jīng)驗(yàn)對(duì)樣本信號(hào)進(jìn)行人工拾取，其中P波到時(shí)如圖2中直線。拾取得P波到時(shí)在[161.65,161.675]s之間，精確P波到時(shí)為161.652 8 s。

圖2 人工拾取P波

對(duì)樣本信號(hào)采用STA/LTA和ML方法處理,P波到時(shí)分別為161.618 7 s和161.631 5 s,如圖3。

圖3 2種方法處理結(jié)果

AMPAT方法得到P波到時(shí)為161.632 6 s。另取IRIS 2016年3月份美國(guó)中部海岸的10組地震數(shù)據(jù)，用這4種方法拾取P波到時(shí)如表1?？芍?，以人工拾取結(jié)果為參考，STA/LTA方法拾取誤差比較大，最大為0.338 8 s，ML拾取結(jié)果比STA/LTA方法好，最大誤差為0.061 4 s,AMPAT拾取結(jié)果比STA/LTA和ML結(jié)果要好，最大誤差為0.055 4 s。進(jìn)一步比較，以人工拾取方法結(jié)果為參照，其他3種方法計(jì)算平均誤差為：STA/LTA，0.073 7 s；ML，0.0331 s;AMPAT，0.025 1 s。可知，AMPAT的平均誤差要小于ML方法和STA/LTA方法的，ML拾取誤差較好于STA/LTA方法。

表1 拾取10組地震數(shù)據(jù)的P波初至?xí)r間 s

綜上，AMPAT在最大誤差和平均誤差都優(yōu)于STA/LTA和ML方法，其對(duì)P波拾取的精度高于這2種方法。

3 結(jié) 論

針對(duì)地震數(shù)據(jù)單個(gè)特征受噪聲等其他因素影響準(zhǔn)確拾取P波到時(shí)難度之大這一問題，本文提取地震數(shù)據(jù)的15維特征并通過基于LLE降維和SVM拾取P波，實(shí)驗(yàn)對(duì)比證明，AMPAT方法獲取的P波到時(shí)更加精確。由于AMPAT方法選取了15維特征，導(dǎo)致計(jì)算量較多，如何有效降低計(jì)算復(fù)雜度是下一階段的研究目標(biāo)。