基于Volterra 級數(shù)的核爆地震參數(shù)化非線性特征提取方法*

韓紹卿，李夕海，伍海軍，劉代志

(第二炮兵工程學(xué)院，陜西 西安710025)

地震波方法是監(jiān)測核實驗的重要方法［1］。分布于全球范圍內(nèi)的世界標(biāo)準(zhǔn)地震臺網(wǎng)以及各個國家和地區(qū)的地震臺網(wǎng)可以及時有效地監(jiān)測到某一地區(qū)的地震事件。由地震儀獲取的地震波形信號為一時間序列，該時間序列中包含了有關(guān)震源的豐富信息。地下核爆炸的自動監(jiān)測問題是一個模式識別問題，解決該問題的關(guān)鍵在于如何對地震時間序列提取出有效的特征，將核爆炸事件同其他干擾事件(如天然地震、化學(xué)爆炸)區(qū)分開來。

在對時間序列信號的分析和特征提取中，目前常用的方法可分為兩類［2-3］:一類是非參數(shù)信號表示技術(shù)，比如傅立葉變換，各種時頻分析技術(shù);另一類是參數(shù)化的信號建模技術(shù)。對于地震波時間序列的特征提取，采用的非參數(shù)化分析方法主要有:傅立葉變換［4］、短時傅立葉變換［5］、魏格納分布［6］、小波變換［7］等;參數(shù)化信號建模則主要采用ARMA(autoregressive moving average)模型［8］。參數(shù)化核爆地震特征提取的基本思想是:首先對地震波時間序列建立ARMA 模型，然后提取模型的參數(shù)作為特征進行分類判別。對于線性時不變最小相位系統(tǒng)，ARMA 模型是在實際中獲得廣泛應(yīng)用的時間序列信號模型。該模型是建立在觀測信號為平穩(wěn)、高斯時間序列的基礎(chǔ)之上的。但是對于地震波信號的分析和處理來說，大地系統(tǒng)不是一個線性系統(tǒng)，時不變和最小相位特性通常是不滿足的，地震波信號也并非是平穩(wěn)的滿足高斯分布的信號［9］。因此基于ARMA 模型的核爆地震特征的識別率不高。研究顯示，地震波信號是一種非線性信號，且具有混沌特性［10-11］，所以采用非線性方法對地震波信號進行建模才更合理。Volterra 級數(shù)是一種泛函級數(shù)，大多數(shù)非線性動態(tài)系統(tǒng)都可以用Volterra 級數(shù)逼近到任意精確的程度［12］。

本文中，基于混沌時間序列分析中的相空間重構(gòu)理論，采用Volterra 自適應(yīng)預(yù)測模型對地震波時間序列進行建模，進而提取模型的參數(shù)作為地震波信號的特征進行分類識別，獲得比ARMA 模型參數(shù)更好的分類效果。研究表明，地震波信號中的線性、非線性和高階統(tǒng)計信息在核爆地震分類中發(fā)揮著重要作用。

1 基于模型的地震波非線性特征提取方法

提出的核爆地震非線性特征提取的核心思想是:基于地震波信號的混沌特性，對地震時間序列進行相空間重構(gòu)，然后采用非線性級數(shù)Volterra 級數(shù)建立自適應(yīng)預(yù)測模型，利用模型參數(shù)的信息凝聚性獲得地震波信號的低維非線性特征。該方法的關(guān)鍵是地震波信號的相空間重構(gòu)和重構(gòu)相空間內(nèi)Volterra 自適應(yīng)預(yù)測模型的建立。

1.1 混沌時間序列的Volterra 自適應(yīng)預(yù)測模型

混沌時間序列預(yù)測的基礎(chǔ)是時間延遲相空間重構(gòu)理論，假設(shè)觀測到的混沌時間序列為x(n)，n=1，2，…，N，通過延遲坐標(biāo)重構(gòu)可獲得如下的延遲坐標(biāo)向量

式中:m 為嵌入維數(shù)，τ 為延遲時間。鑒于實際中大量的非線性系統(tǒng)都可用Volterra 級數(shù)得到很好的表征，采用Volterra 級數(shù)展開式來構(gòu)造混沌時間序列的非線性預(yù)測模型［13］。設(shè)非線性離散動力系統(tǒng)的輸入為X(n)，輸出為y(n)=^x(n+1)，則該非線性系統(tǒng)函數(shù)的Volterra 級數(shù)展開式為

式中:hp(m1，m2，…，mp)為p 階Volterra 核。

這種無窮級數(shù)展開式在實際應(yīng)用中無法實現(xiàn)，必須采用有限截斷和有限求和的形式。在二階截斷求和的情況下Volterra 級數(shù)展開式為

實際應(yīng)用中，濾波器的長度N1、N2必須取有限長。根據(jù)Takens 嵌入定理，可將N1、N2均取為N1=N2=m。式(3)為一非線性模型，線性擴展后可用線性自適應(yīng)FIR 濾波器實現(xiàn)。定義濾波器的輸入信號矢量為

以及濾波系數(shù)矢量

則式(3)可表示為

對于式(6)這種自適應(yīng)濾波器，可采用最小均方誤差自適應(yīng)算法來求解。

以上簡要闡述了二階截斷情況下的Volterra 自適應(yīng)預(yù)測模型，同理可將其擴展到更高階，在此不再贅述。

1.2 地震波非線性特征提取方法

核爆地震和天然地震信號是混沌信號，但其混沌參數(shù)(如分維數(shù))不適合作模式識別的特征參數(shù)［11］。本文中則根據(jù)地震波的混沌特性，利用Volterra 級數(shù)建立地震波時間序列的自適應(yīng)預(yù)測模型，在此基礎(chǔ)上提出一種參數(shù)化的核爆地震非線性特征提取算法，該算法的流程如圖1 所示。

圖1 核爆地震非線性特征提取算法流程圖Fig.1 Flow chart of nonlinear feature extraction algorithm of nuclear explosions and earthquakes

在特征提取過程中，非常關(guān)鍵的一個步驟是混沌系統(tǒng)參數(shù)——時間延遲和嵌入維的確定。在一般的混沌時間序列(特別是基于混沌動力學(xué)模型的時間序列)分析中，已經(jīng)發(fā)展了多種估計時間延遲和嵌入維的方法。但是在核爆地震特征提取中，這些方法不太適用，原因如下:(1)不同的參數(shù)估計方法獲得的參數(shù)往往不同，實踐中到底該采用哪種方法尚無嚴(yán)格標(biāo)準(zhǔn)。(2)在進行參數(shù)估計時，大部分情況下數(shù)據(jù)量要求都比較大，而核爆炸是在瞬間完成的，一般情況下，天然地震的持續(xù)時間也不會很長，再加上地震儀器采樣率的限制，實際中獲得地震波信號長度都很有限。(3)實際中獲得的地震波信號都會受到噪聲的污染，但很多情況下對噪聲的統(tǒng)計分布和強度情況了解較少，此外，即使獲得了噪聲的先驗信息，但降噪會給原動力系統(tǒng)的信息造成損失。(4)核爆地震的分類必然涉及大量的樣本，即使是同類樣本之間的混沌參數(shù)也會有一定的差異。特征提取的目的在于找到同類樣本的共性和不同類樣本的差異，而沒有必要一定要獲得對于每個單個樣本最佳的混沌參數(shù)。此外，樣本的混沌參數(shù)不同會造成樣本特征空間維數(shù)的不同，給后續(xù)樣本的分類帶來很大麻煩。這正是模式識別過程中特征提取與一般信號分析的不同點。為解決這一問題，采用如下的方法:首先采用自相關(guān)和去偏復(fù)自相關(guān)法確定每一個地震波樣本的時間延遲，經(jīng)過綜合分析得出地震波樣本的時間延遲的分布范圍，采用Grassberger 和Procaccia 提出的G-P 法和Cao 氏法確定每一個地震波樣本的嵌入維，經(jīng)過綜合分析得出地震波樣本的嵌入維的分布范圍，然后再將二者的取值范圍適當(dāng)增大，在增大的取值范圍內(nèi)對所有的核爆和地震樣本都分別采用相同的時間延遲和相同的嵌入維進行自適應(yīng)預(yù)測建模，最后進行特征提取和分類判別，根據(jù)識別率確定最佳的時間延遲和嵌入維。這樣處理的依據(jù)和好處是:(1)Takens 定理表明，只要相空間嵌入維大于飽和嵌入維時，都可以得到原動力系統(tǒng)的重構(gòu)。此外定理中嵌入維應(yīng)當(dāng)滿足的條件只是一個充分條件，而不是必要條件，并不能排除在比這一要求更小的嵌入維時得到較好的重構(gòu)。這說明滿足重構(gòu)的嵌入維是在一個范圍內(nèi)變化的。(2)延遲時間的確定以相關(guān)性為準(zhǔn)則，希望相鄰點間既不能不相關(guān)，又不能相關(guān)性太強。這種矛盾表明延遲時間也是在一定范圍內(nèi)滿足要求。(3)現(xiàn)有計算方法獲得的時間延遲和嵌入維相對于相空間重構(gòu)來說是最佳的，但對于特征提取和分類而言不一定最佳。而本文中提出的在一定范圍內(nèi)經(jīng)過比較后獲得最佳時間延遲和嵌入維，其判別標(biāo)準(zhǔn)是識別率，這樣能保證獲得對于分類而言最佳的參數(shù)。

地震臺站獲得的地震波信號是震源(爆炸、地震)發(fā)出的地震波經(jīng)過大地系統(tǒng)多次濾波后產(chǎn)生的時間序列。由于震源激發(fā)地震波的過程本身是一個非線性過程，且大地系統(tǒng)是一個非線性系統(tǒng)，因此地震波信號是非線性信號。有關(guān)震源的豐富信息蘊含在地震波時間序列中，當(dāng)建立地震波時間序列的非線性模型之后，震源信息也同樣蘊含在模型參數(shù)之中。參數(shù)模型的一個最大特點是信息的凝聚性，即大量數(shù)據(jù)所蘊含的信息凝聚成少數(shù)幾個模型參數(shù)。而特征提取的目的之一也是為了實現(xiàn)信息的凝聚，從而能夠使分類判決在維數(shù)較低的空間進行。這是本文中提出的參數(shù)化非線性特征提取方法的理論依據(jù)。

2 實驗結(jié)果

實驗所用的數(shù)據(jù)庫為北京、蘭州、恩施、昆明、瓊中、上海、烏魯木齊、海拉爾、牡丹江等9 個中周期地震臺站所采集的核爆炸和天然地震所產(chǎn)生的地震波數(shù)據(jù)，各有54 個樣本，樣本的采樣率為20 Hz。為了計算方便和便于對比，對每個樣本從地震波初至點開始截取1 024 個數(shù)據(jù)點。圖2 給出了一次典型的核爆炸樣本和一次典型的天然地震樣本的時域波形圖?？梢钥闯觯瑢τ诘湫偷暮吮ê吞烊坏卣?，二者時域波形有較為顯著的差異。但是通過對數(shù)據(jù)庫中樣本的統(tǒng)計分析發(fā)現(xiàn)，并非所有的核爆炸樣本和天然地震樣本之間都具有這樣的差異。因此必須通過特征提取，將核爆炸和天然地震之間的差異量化。此外，不同樣本的幅度有很大差異。造成這種差異的主要原因是:不同的地震事件震級不同;不同的地震事件距離接收臺站的距離不同。為了消除波形幅度對判據(jù)的影響，特征提取算法中首先需要對數(shù)據(jù)進行歸一化處理。

圖2 核爆炸和天然地震產(chǎn)生的地震波形Fig.2 Seismic waveforms from nuclear explosion and natural earthquakes

根據(jù)圖1 中的特征提取算法流程圖我們進行了核爆地震非線性特征提取實驗。在驗證特征的有效性時，采用k-近鄰分類器進行分類實驗，其中錯誤識別率采用留一法獲得。在對地震波時間序列進行建模時，需要確定的參數(shù)為Volterra 級數(shù)的階數(shù)p、相空間重構(gòu)的時間延遲τ 以及濾波器長度(即嵌入維數(shù)m)。隨著階數(shù)p 的增加，Volterra 級數(shù)濾波器系數(shù)的個數(shù)迅速增加，導(dǎo)致計算量飛速增長。因此僅在Volterra 級數(shù)的階數(shù)p=2，3 時展開研究。對于時間延遲τ，經(jīng)過綜合分析得到的取值范圍為［2，8］，其中在τ 取6 時獲得最佳的識別效果。對于嵌入維數(shù)m，經(jīng)過綜合分析得到的取值范圍為［2，11］。為了考察濾波器長度對特征有效性的影響，同時也為了與已有的ARMA 模型參數(shù)特征進行比較，在τ=6、嵌入維m 取值范圍為［2，11］的情況下進行特征提取和分類判別實驗。從形式上看，本文中采用的Volterra 級數(shù)模型為自回歸模型，因此將ARMA 模型中的自回歸模型(AR 模型)參數(shù)特征作為比較的對象，其中AR 模型的階數(shù)取值范圍同樣為［2，11］。

不同模型參數(shù)特征的分類結(jié)果如圖3 所示。為了便于分析和對比，將AR 模型參數(shù)特征和Volterra自適應(yīng)預(yù)測模型參數(shù)特征的分類識別結(jié)果畫在同一張圖中。圖中縱坐標(biāo)E 代表錯誤識別率。需要注意的是，對于AR 模型，圖中的橫坐標(biāo)代表的是不同的階數(shù)(用l 表示)，而對于Volterra 自適應(yīng)預(yù)測模型，橫坐標(biāo)代表的是不同的嵌入維。從圖中可以看出，AR 模型參數(shù)特征識別效果不太理想，在p=3 時達到最低錯誤識別率，在0.3 左右，并且并非階數(shù)越高識別效果越好。對于本文提出的非線性模型參數(shù)特征，在p=2 時，除嵌入維m=2，5 的情況外，其錯誤識別率都低于AR 模型參數(shù)特征的最低錯誤識別率，在m=9 時達到最低錯誤識別率，為0.213。在p=3 時，非線性模型參數(shù)特征在所有嵌入維情況下的錯誤識別率都低于AR 模型參數(shù)特征的最低錯誤識別率，并且在m=7 時達到最低錯誤識別率，為0.194 4。總體來看，本文中提出的非線性模型參數(shù)特征優(yōu)于線性模型的AR 模型參數(shù)特征。在非線性模型參數(shù)特征中，與階數(shù)為2 情況下獲得的特征相比，p=3 情況下獲得的特征識別效果更好，且對于嵌入維數(shù)更加不敏感。

圖3 不同特征提取方法的分類性能比較Fig.3 Classification performance comparison for different feature extraction methods

3 結(jié)果分析

由實驗結(jié)果可以看出，本文中提出的非線性模型參數(shù)特征的分類效果好于現(xiàn)有的線性AR 模型參數(shù)特征的分類效果，根本原因在于Volterra 自適應(yīng)預(yù)測模型很好地刻畫了地震時間序列非線性本質(zhì)特征以及所具有的混沌特性。具體原因可有如下幾點:(1)地震時間序列本質(zhì)上是非線性的，AR 模型只利用了地震波的線性因素，而Volterra 自適應(yīng)預(yù)測模型綜合利用了線性和非線性因素。(2)地震時間序列是非平穩(wěn)的時間序列，用于刻畫非平穩(wěn)時間序列的模型參數(shù)應(yīng)該是時變的。AR 模型參數(shù)是固定不變的，而Volterra 自適應(yīng)預(yù)測模型能夠根據(jù)地震波統(tǒng)計特性的改變而自適應(yīng)地調(diào)整模型參數(shù)，因此后者對地震波的逼近精度更高。(3)AR 模型只利用了地震波的二階統(tǒng)計信息，而Volterra 自適應(yīng)預(yù)測模型充分利用了地震波的高階統(tǒng)計信息。高階統(tǒng)計信息在核爆地震分類中是有用的，三階Volterra 自適應(yīng)預(yù)測模型參數(shù)特征優(yōu)于二階Volterra 自適應(yīng)預(yù)測模型參數(shù)特征也充分說明了這一點。

4 結(jié) 論

由于地震時間序列為非線性時間序列，因此需要采用非線性模型對其進行準(zhǔn)確描述。利用地震時間序列的混沌特性，先進行相空間重構(gòu)，然后在重構(gòu)的相空間內(nèi)建立地震波的Volterra 自適應(yīng)預(yù)測模型，進而提取模型的參數(shù)作為核爆炸和天然地震的特征進行分類實驗，獲得了較好的分類效果。通過研究獲得的主要結(jié)論為:基于Volterra 級數(shù)的非線性模型參數(shù)特征優(yōu)于現(xiàn)有的線性AR 模型參數(shù)特征;在非線性模型參數(shù)特征中，三階Volterra 級數(shù)自適應(yīng)預(yù)測模型參數(shù)特征優(yōu)于二階Volterra 級數(shù)自適應(yīng)預(yù)測模型參數(shù)特征。因此，在核爆地震特征提取中，充分考慮和應(yīng)用地震時間序列的非線性和高階統(tǒng)計信息將有助于獲得更有效的特征。

［1］ Hafemeister D.Progress in CTBT monitoring since its 1999 senate defeat［J］.Science and Global Security，2007，15:151-183.

［2］ Sakkalis V，Cassar T，Zervakis M，et al.Parametric and nonparametric EEG analysis for the evaluation of EEG activity in young children with controlled epilepsy［J］.Computational Intelligence and Neuroscience，2008，doi:10.1155/2008/462593.

［3］ Carden E P，Brownjohn J M W.ARMA modelled time-series classification for structural health monitoring of civil infrastructure［J］.Mechanical Systems and Signal Processing，2008，22(2):295-314.

［4］ 李夕海，劉剛，劉代志，等.基于最近鄰支撐向量特征線融合算法的核爆地震識別［J］.地球物理學(xué)報，2009，52(7):1816-1824.LI Xi-hai，LIU Gang，LIU Dai-zhi，et al.Discrimination of nuclear explosions and earthquakes using the nearest support vector feature line fusion classification algorithm［J］.Chinese Journal of Geophysics，2009，52(7):1816-1824.

［5］ Arrowsmith M D，Arrowsmith S J，Stump B W，et al.A suite of discriminants for ground-truth mining events in the western US and its implications for discrimination capability in Russia［C］∥Proceedings of 2008 Monitoring Research Review:Ground-Based Nuclear Explosion Monitoring Technologies.Portsmouth，Virginia，USA，2008:544-553.

［6］ Benbrahim M，Benjelloun K，Ibenbrahim A，et al.A new approaches for seismic signals discrimination［J］.Proceedings of World Academy of Science，Engineering and Technology，2007，21:183-186.

［7］ Piotr F，Hernando O.Consistent classification of non-stationary time series using stochastic wavelet representations［J］.Journal of the American Statistical Association，2009，104(485):299-312.

［8］ Cercone J A，F(xiàn)oster V S，Clark W M.Application of neural networks to seismic signal discrimination research findings［R］.PL-TR-94-2122，1994.

［9］ 邢貞貞，韓立國，王宇，等.高階統(tǒng)計量方法在地震信號分析中的應(yīng)用［J］.吉林大學(xué)學(xué)報(地球科學(xué)版)，2007，37 增刊:139-142.XING Zhen-zhen，HAN Li-guo，WANG Yu，et al.Application of higher-order statistics in seismic signal analysis［J］.Journal of Jilin University(Earth Science Edition)，2007，37 suppl:139-142.

［10］ Liu D，Red S，Wei Y，et al.Attractor analysis and its applications to seismic pattern recognition of nuclear explosion［C］∥3rd International Conference on Signal Processing.Beijing，China，1996:1324-1329.

［11］ 劉代志，鄒紅星，韋蔭康，等.分形分析與核爆地震模式識別［J］.模式識別與人工智能，1997，10(2):153-158.LIU Dai-zhi，ZOU Hong-xing，WEI Yin-kang，et al.Fractal analysis with applications to seismic pattern recognition of nuclear explosion［J］.Pattern Recognition and Artificial Intelligence，1997，10(2):153-158.

［12］ Kalli R R.Nonlinear modeling of radio frequency circuits to estimate third-order nonlinear distortions［D］.Florida:Florida International University，2008.

［13］ 許小可.海雜波的非線性分析與建模［D］.大連:大連海事大學(xué)，2007.