郝國成，鍋 娟，譚淞元，曾佐勛

(1.中國地質(zhì)大學(武漢) 機械與電子信息學院，湖北 武漢 430074； 2.中國科學院 測量與地球物理研究所 大地測量與地球動力學國家重點實驗室，湖北 武漢 430077； 3.中國地質(zhì)大學(武漢)復雜系統(tǒng)先進控制與智能自動化湖北省重點實驗室，湖北 武漢 430074；4.中國地質(zhì)大學(武漢) 智能地學信息處理湖北省重點實驗室，湖北 武漢 430074； 5.中國地質(zhì)大學(武漢) 地球科學學院，湖北 武漢 430074)

地球天然脈沖電磁場(the Earth’s natural pulse electromagnetic field, ENPEMF)是指可在地球表面接收的由天然場源產(chǎn)生的綜合電磁總場[1-2].地震、滑坡等地質(zhì)災害現(xiàn)象可在地表產(chǎn)生甚低頻(very low frequency, VLF)信號脈沖波動，“微破裂機-電轉(zhuǎn)換”機制和“地殼波導”是上述電磁現(xiàn)象的機理之一[3].ENPEMF信號中的電磁異常信息具有潛在變化趨勢及典型的非平穩(wěn)特征，可反映地質(zhì)活動的孕育發(fā)展趨勢，可用于震前電磁異常監(jiān)測分析[4].

隨著信息技術(shù)、人工智能及機器學習理論的不斷發(fā)展，電磁信息預測模型逐漸成為前沿熱點.由于影響ENPEMF信號的場源較多，具有非平穩(wěn)信號的特點，很難根據(jù)脈沖波形歸納其規(guī)律.非線性模型中的神經(jīng)網(wǎng)絡方法是基于經(jīng)驗風險最大化原則的機器學習算法，非線性擬合能力較強，可以對采集到的震前ENPEMF信號進行建模，并擬合其強度趨勢的變化特點.

徑向基函數(shù)(radial basis function, RBF)神經(jīng)網(wǎng)絡可逼近任意的非線性函數(shù)，其優(yōu)點是學習速度快、非線性逼近能力強，具有良好泛化能力，可為ENPEMF信號構(gòu)建非線性預測模型.本文提出基于混沌參數(shù)優(yōu)化RBF算法的預測模型對其強度趨勢進行預測，為數(shù)據(jù)分析和災害監(jiān)測提供支持.首先采用混沌理論對實測ENPEMF數(shù)據(jù)進行分析，其中假鄰近法(false nearest neighbor, FNN)及自相關(guān)函數(shù)法分別求得嵌入維數(shù)和延遲時間等混沌特征參數(shù)，并優(yōu)化RBF神經(jīng)網(wǎng)絡；采用訓練完成的混沌參數(shù)優(yōu)化RBF神經(jīng)網(wǎng)絡模型對ENPEMF數(shù)據(jù)進行預測，并與傳統(tǒng)RBF神經(jīng)網(wǎng)絡進行比較.結(jié)果表明，混沌參數(shù)RBF算法可基于混沌理論確定動態(tài)系統(tǒng)的混沌特性，從采集到的震前ENPEMF信號強度數(shù)據(jù)中找到其變化態(tài)勢，并預測14 d(4月7～20日)的ENPEMF數(shù)據(jù)強度趨勢，且預測效果及精度均優(yōu)于傳統(tǒng)RBF預測模型.

1 混沌理論

由于震前ENPEMF信號的產(chǎn)生機理復雜、孕育過程非線性，其信號強度數(shù)據(jù)具有非平穩(wěn)特點和混沌特性.因此，本文引入混沌理論對其數(shù)據(jù)內(nèi)部特征進行挖掘，找到隱藏的混沌特點.假設(shè)一段時間內(nèi)采集的震前ENPEMF信號數(shù)據(jù)為{x(tj),j=1,2,…,n}，其中n表示采集的數(shù)據(jù)點數(shù)，通過混沌理論中的假鄰近法及自相關(guān)函數(shù)法對數(shù)據(jù)進行處理，得到數(shù)據(jù)變化形式：

X(t)=[x(t),x(t+τ),…,x(t+(m-1)τ)] .

(1)

式中，τ和m分別表示震前ENPEMF信號數(shù)據(jù)的延遲時間和嵌入維數(shù)，用于描述該信號隱藏的混沌特征并為RBF神經(jīng)網(wǎng)絡輸入節(jié)點提供判斷依據(jù).

1.1 假鄰近法

時間序列的本質(zhì)是將系統(tǒng)高維空間坐標的運動軌跡投影到低維空間.當嵌入維數(shù)較小時，系統(tǒng)空間軌道中本來相距很遠的相點相互擠壓折疊，未能充分展開，這些點為假鄰近點[5].隨著嵌入空間維數(shù)的增加，軌道逐漸展開，投影到低維空間的假鄰近點隨之分離.當所有的假鄰近點消失時所對應的最小嵌入空間維數(shù)即為最佳嵌入維數(shù).給定正整數(shù)m，可構(gòu)造m維重構(gòu)向量：

ym(n)=(x(n),x(n+τ),…,
x(n+(m-1)τ))T.

(2)

(3)

將維數(shù)從1維增加到m+1維.m+1維空間重構(gòu)向量如式(4)所示：

ym+1(n)=(x(n),x(n+τ),…,
x(n+(m-1)τ),x(n+mτ))T.

(4)

(5)

(6)

(7)

(8)

對所有的重構(gòu)向量，利用判據(jù)找出鄰近點中的假鄰近點，并記錄下所有假鄰近點的數(shù)目FN(m).繼續(xù)增加維數(shù)，當找到一個整數(shù)mε使得FN(mε)=0時，mε即為所求嵌入維數(shù).當假鄰近點所占比率即假鄰近率隨著嵌入維數(shù)的增加趨于平穩(wěn)不再降低時，所對應的嵌入維數(shù)m為最佳嵌入維數(shù).本文通過統(tǒng)計假鄰近點數(shù)的比率隨嵌入維數(shù)升高逐漸減小，最后維持不變的情況，確定最優(yōu)嵌入維數(shù).

1.2 自相關(guān)函數(shù)法

由于實際時間序列長度有限且存在噪聲，選取合適的延遲時間至關(guān)重要.延遲時間τ過小，將使重構(gòu)的系統(tǒng)由于相關(guān)性較強造成相空間的擠壓，不能充分展示系統(tǒng)的動力特征；延遲時間τ太大，會造成相鄰兩時刻的動力學形態(tài)劇烈變化，使構(gòu)造的相空間比實際空間復雜.

自相關(guān)函數(shù)法可在降低相關(guān)性的同時保證原動力學的系統(tǒng)信息不丟失，使重構(gòu)相空間能充分展現(xiàn)系統(tǒng)拓撲性質(zhì)和幾何性質(zhì).首先寫出時間序列的自相關(guān)函數(shù)，然后做出自相關(guān)函數(shù)隨時間變化的函數(shù)圖，找到自相關(guān)函數(shù)首次達到零點時對應的時間即為延遲時間τ.自相關(guān)函數(shù)定義為

(9)

自相關(guān)函數(shù)值隨時間變化逐步下降，當其下降到初始值的(1-1/e)時對應的時間為所求延遲時間.

2 RBF神經(jīng)網(wǎng)絡算法

徑向基函數(shù)(RBF)神經(jīng)網(wǎng)絡是一種包含輸入層、隱含層和輸出層的三層神經(jīng)網(wǎng)絡[6].其中輸入層與隱含層之間為非線性變換，從隱含層到輸出層為線性變換.

在RBF神經(jīng)網(wǎng)絡中，輸入層僅作為通道傳輸信號.隱含層中神經(jīng)元的變換函數(shù)為徑向基函數(shù)，通過非線性變換可將信號從輸入層傳遞到隱含層.輸出層是對輸入信號的響應.RBF神經(jīng)網(wǎng)絡結(jié)構(gòu)可根據(jù)具體問題在訓練階段進行自適應調(diào)整.

在網(wǎng)絡中，X=(x1,x2,…,xn)T為輸入樣本，Y=(y1,y2,…,yn)T為輸出響應.RBF神經(jīng)網(wǎng)絡算法需要求解基函數(shù)的中心Ci，基函數(shù)的寬度Di以及隱含層到輸出層的權(quán)值Wi三個參數(shù).

RBF神經(jīng)網(wǎng)絡的訓練過程分為兩步：首先進行無監(jiān)督學習，計算輸入層與隱含層之間的和，得到隱含層輸出：

(10)

其中：i=1,2,…,N;‖X-Ci‖為歐氏范數(shù).

然后采用最小二乘法求隱含層與輸出層之間的權(quán)值ωi.最終得到RBF神經(jīng)網(wǎng)絡的輸出y：

(11)

將經(jīng)過參數(shù)訓練的RBF神經(jīng)網(wǎng)絡用于預測混沌時間序列.可利用混沌時間序列數(shù)據(jù)校正上述權(quán)值參數(shù)，提高神經(jīng)網(wǎng)絡的非線性泛化能力.其中，當輸入層節(jié)點個數(shù)為混沌時間序列的嵌入維數(shù)m時，預測結(jié)果較好；隱含層節(jié)點數(shù)目根據(jù)實驗實時調(diào)整確定：預先設(shè)定RBF神經(jīng)網(wǎng)絡的精度值，隱含層節(jié)點個數(shù)遞增，當神經(jīng)網(wǎng)絡達到預設(shè)精度時，該節(jié)點個數(shù)即為神經(jīng)網(wǎng)絡的隱含層節(jié)點數(shù).

本文提出的基于混沌參數(shù)優(yōu)化RBF神經(jīng)網(wǎng)絡算法的預測模型工作流程如圖1所示.

圖1 本文算法預測模型工作流程

首先對實測ENPEMF數(shù)據(jù)進行混沌分析，然后分別采用假鄰近法和自相關(guān)函數(shù)法求得最優(yōu)嵌入維數(shù)和延遲時間等混沌特征參數(shù)，將得到的參數(shù)作為確定輸入層節(jié)點個數(shù)的依據(jù)并優(yōu)化RBF神經(jīng)網(wǎng)絡；對參數(shù)優(yōu)化的RBF神經(jīng)網(wǎng)絡進行訓練，學習其內(nèi)部混沌特征，最后用訓練完成的混沌參數(shù)優(yōu)化RBF算法預測ENPEMF信號的強度趨勢.

3 實驗仿真

隱藏在混沌時間序列內(nèi)部的某些特征信息可以通過延遲時間及嵌入維數(shù)表現(xiàn)出來.本文用經(jīng)典的Rossler混沌時間序列驗證假鄰近法及自相關(guān)函數(shù)法獲得嵌入維數(shù)和延遲時間等參數(shù)的可行性.

Rossler系統(tǒng)可用式(12)微分方程組進行描述：

(12)

選取參數(shù)a=b=0.2，c=5，初值x(0)=-1，y(0)=0，z(0)=1，積分時間步長h=0.05，生成長度為3 000的連續(xù)混沌時間序列如圖2所示.

圖2 Rossler混沌時間序列

Rossler混沌時間序列在有限區(qū)域內(nèi)運動時趨向于一個穩(wěn)定的點，完全展開系統(tǒng)內(nèi)部的混沌特性.嵌入維數(shù)的取值范圍為[1,8]，閾值的判別門限范圍為[2,15]，Rossler時間序列長度為3 000，采用假鄰近法計算時間序列x分量的嵌入維數(shù)m，結(jié)果如圖3所示.

當嵌入維數(shù)從1增加到4時假鄰近率急速下降；當嵌入維數(shù)達到5時，假鄰近率趨于平緩，此時的嵌入維數(shù)達到理想值，為所求最佳嵌入維數(shù)，即Rossler混沌時間序列的嵌入維數(shù)為5.

采用自相關(guān)函數(shù)法求延遲時間τ，選擇x分量序列計算延遲時間τ，仿真結(jié)果如圖4所示.

圖4中，直線為初始值的(1-1/e)，曲線為自相關(guān)函數(shù)曲線，選取值為自相關(guān)函數(shù)下降到初始值的(1-1/e)時所對應的時間，即Rossler混沌時間序列的延遲時間τ為16 s.

圖3 假鄰近法求嵌入維數(shù)

圖4 自相關(guān)函數(shù)法求延遲時間

取迭代后長度為2 000的時間序列，其中前500個數(shù)據(jù)進行訓練，后1 500個數(shù)據(jù)進行預測.以嵌入維數(shù)m=5作為RBF神經(jīng)網(wǎng)絡的輸入節(jié)點判斷依據(jù)，采用訓練完成的混沌參數(shù)優(yōu)化RBF算法對數(shù)據(jù)進行預測.圖5為預測模型的仿真結(jié)果.“×”表示真實值，“-”表示預測值.圖6為預測結(jié)果的絕對誤差.

圖5 Rossler混沌時間序列的預測結(jié)果

圖6 預測結(jié)果的絕對誤差

圖5和圖6中，混沌參數(shù)優(yōu)化RBF算法可以較好擬合經(jīng)典的Rossler混沌時間序列，且預測誤差較小.

4 混沌參數(shù)優(yōu)化RBF算法的預測研究

本文使用俄羅斯科學院托木斯克分院GR-01型設(shè)備接收ENPEMF信號.在武漢九峰地震臺放置了3臺設(shè)備，接收方向為W-E和N-S，3個通道為CN1，CN2，CN3.設(shè)備記錄了ENPEMF信號的AH數(shù)據(jù)(超過設(shè)定閾值的脈沖幅度)及NH數(shù)據(jù)(超過設(shè)定閾值的脈沖個數(shù))，可表征地表天然磁場的強弱[7].設(shè)備的工作頻率為甚低頻段：5～25 kHz，在武漢設(shè)置的接收頻率為14.5 kHz.

ENPEMF信號為非周期、非平穩(wěn)信號，具有明顯混沌特性，設(shè)備輸出為數(shù)字量化后的信號，數(shù)據(jù)存儲格式為：時間-幅度-脈沖數(shù)(t-AH-NH)[8].幅度單位僅為信號包絡大小變化的參照量，已不具有原來量綱的直接意義.

對采集的ENPEMF數(shù)據(jù)進行混沌特性分析，找到其內(nèi)部隱藏的混沌特征及趨勢變化特點，結(jié)合RBF神經(jīng)網(wǎng)絡算法對其信號強度進行預測，識別孕震信息[9].本文采用Grassberger和Procaccia提出的關(guān)聯(lián)維數(shù)算法(G-P算法)求解ENPEMF信號關(guān)聯(lián)維數(shù)，可判斷信號是否具有混沌特性[10].

對于時間序列x(1),x(2),…,x(t)，其長度為M，對其進行相空間重構(gòu)，得到向量：X(t)=[x(t),x(t+τ),…,x(t+(m-1)τ)],其中t=1,2,…,N,N=M-(m-1)τ.給定正數(shù)ε足夠小，當空間向量間的距離小于ε時，向量關(guān)聯(lián).關(guān)聯(lián)向量的關(guān)聯(lián)積分表達式為

(13)

其中θ(·)為Heaviside階躍函數(shù)，滿足式(14)：

(14)

當時間序列的長度N→∞，半徑ε→0時，關(guān)聯(lián)積分與半徑的關(guān)系為

(15)

其中，D為所求關(guān)聯(lián)維數(shù)，變形后得到式(16)：

(16)

給定一系列半徑ε和嵌入維數(shù)m，作半徑隨嵌入維數(shù)變化的關(guān)聯(lián)積分圖組，用最小二乘法對圖中l(wèi)bC(N,ε)～lbε最接近直線的一段擬合最佳直線，該直線斜率即所求關(guān)聯(lián)維數(shù)D.

本文數(shù)據(jù)于2013年4月20日在中國四川省蘆山地震期間收集，地震的位置和地震臺的位置如圖7所示[11].

圖7 蘆山地震的位置和地震臺的位置

圖8為4月10～20日通道2的AH數(shù)據(jù)，箭頭指向中國四川省7.0級蘆山地震發(fā)生時間.

圖8 4月10～20日通道2的AH數(shù)據(jù)

地震發(fā)生前的11 d內(nèi)，ENPEMF數(shù)據(jù)從14日到15日有較大的峰值變化，在16日回落至正常水平.在17日和18日觀測到顯著的峰值變化，在19日信號又跌落至正常.因此，在地震發(fā)生前ENPEMF信號脈沖強度會發(fā)生劇烈變化.

本文選擇經(jīng)過平滑及歸一化的20 d ENPEMF數(shù)據(jù)作為實驗數(shù)據(jù)(4月1～20日).其中前6 d數(shù)據(jù)(4月1～6日)為模型的訓練樣本，后14 d數(shù)據(jù)(4月7～20日)作為模型的預測樣本.利用G-P算法計算4月1～6日ENPEMF信號的關(guān)聯(lián)維數(shù).信號的關(guān)聯(lián)積分組在一定范圍內(nèi)呈近似直線分布；隨著嵌入維數(shù)的增加，直線斜率增大，且最后關(guān)聯(lián)維數(shù)趨于穩(wěn)定，說明ENPEMF信號具有混沌特性.

對選取的ENPEMF信號中前6 d數(shù)據(jù)(4月1～6日)進行預處理，采用假鄰近法計算嵌入維數(shù)m，自相關(guān)函數(shù)法計算延遲時間τ，得到圖9和10，從中了解震前ENPEMF信號數(shù)據(jù)的混沌特性.

圖9 假鄰近法求嵌入維數(shù)

圖10 自相關(guān)函數(shù)法計算延遲時間

圖9中，隨著嵌入維數(shù)的增加，當嵌入維數(shù)為4時，假鄰近率趨于平穩(wěn)，此時的嵌入維數(shù)即ENPEMF信號數(shù)據(jù)的最佳嵌入維數(shù).圖10中，信號自相關(guān)函數(shù)達到初始值的(1-1/e)時τ=5 s，即ENPEMF信號數(shù)據(jù)的延遲時間為5 s.用得到的參數(shù)確定RBF神經(jīng)網(wǎng)絡的輸入節(jié)點個數(shù)為4，進而將訓練完成的混沌參數(shù)優(yōu)化RBF算法用于數(shù)據(jù)預測.

選擇4月1～6日數(shù)據(jù)作為訓練樣本，訓練混沌參數(shù)優(yōu)化RBF神經(jīng)網(wǎng)絡，其中輸入層有4個節(jié)點，輸出層1個節(jié)點，隱含層有6個節(jié)點.選擇徑向基高斯函數(shù)作為隱含層神經(jīng)元傳遞函數(shù)，輸出為線性函數(shù).最后，利用訓練完成的混沌參數(shù)優(yōu)化RBF神經(jīng)網(wǎng)絡預測模型和傳統(tǒng)的RBF神經(jīng)網(wǎng)絡預測模型分別實現(xiàn)對4月7～20日ENPEMF數(shù)據(jù)的單步預測.圖11為所提混沌參數(shù)優(yōu)化RBF神經(jīng)網(wǎng)絡預測模型的結(jié)果，圖12為傳統(tǒng)RBF神經(jīng)網(wǎng)絡預測模型結(jié)果.

圖11 混沌參數(shù)優(yōu)化RBF預測結(jié)果

圖12 傳統(tǒng)RBF預測結(jié)果

圖11和12中，分別采用混沌參數(shù)優(yōu)化RBF預測模型與傳統(tǒng)RBF預測模型對4月7～20日ENPEMF信號強度數(shù)據(jù)進行預測，兩種模型均可以模擬出采集到的地震前14 d(7～20日)實際ENPEMF信號強度的波動.對于整體數(shù)據(jù)范圍，傳統(tǒng) RBF 神經(jīng)網(wǎng)絡模型不能較好地跟蹤實際值的變化，而本文所提預測模型對ENPEMF信號具有較好的跟蹤擬合性能；對于17日的數(shù)據(jù)劇烈波動時刻，混沌參數(shù)優(yōu)化RBF預測模型相較于傳統(tǒng)的RBF預測模型擬合效果更好，具有較好的預測結(jié)果，誤差較小，預測優(yōu)勢明顯.為更精確評估所提預測模型的預測效果，選取絕對誤差作為ENPEMF數(shù)據(jù)預測精度評價指標，結(jié)果如圖13所示.

圖13 兩種算法的絕對誤差值對比

混沌參數(shù)RBF神經(jīng)網(wǎng)絡預測模型僅在17日劇烈波動時段存在預測誤差，整體上的預測誤差均小于傳統(tǒng) RBF 神經(jīng)網(wǎng)絡預測模型.為驗證混沌參數(shù)優(yōu)化RBF算法預測結(jié)果的穩(wěn)健性和可靠性，本文采用互相關(guān)系數(shù)對兩種算法預測值與實際值之間進行量化測量，如式(17)所示：

(17)

得到混沌參數(shù)優(yōu)化RBF神經(jīng)網(wǎng)絡模型的互相關(guān)系數(shù)結(jié)果為r1=0.800 4，略大于傳統(tǒng)RBF神經(jīng)網(wǎng)絡模型的互相關(guān)系數(shù)r2=0.792 6.因此，本文所提優(yōu)化算法的預測效果優(yōu)于傳統(tǒng)RBF神經(jīng)網(wǎng)絡算法.

綜上，混沌參數(shù)優(yōu)化RBF預測模型能夠較好地反映采集到的強震前14 d(4月7～20日)ENPEMF信號強度變化的趨勢和規(guī)律，可以滿足對強震前ENPEMF信號強度趨勢的預測需要，期望為地震和地質(zhì)災害前的電磁預測發(fā)揮積極的作用.

5 結(jié) 論

1)本文提出了一種混沌參數(shù)優(yōu)化RBF神經(jīng)網(wǎng)絡預測模型.通過混沌理論確定的特征參數(shù)優(yōu)化RBF神經(jīng)網(wǎng)絡，進而對震前ENPEMF強度數(shù)據(jù)進行預測，并與傳統(tǒng)RBF神經(jīng)網(wǎng)絡算法進行比較.

2)本文混沌參數(shù)優(yōu)化RBF預測模型可有效地預測強震前14 d(4月7～20日)ENPEMF信號的強度趨勢，且預測效果優(yōu)于傳統(tǒng)RBF神經(jīng)網(wǎng)絡算法，期望為地震及地質(zhì)災害前的電磁監(jiān)測分析提供依據(jù).