于文強(qiáng) 李厚樸 秦清亮 宋立忠 王志遠(yuǎn)

1 海軍工程大學(xué)電氣工程學(xué)院，武漢市解放大道717號(hào)，430033 2 海洋環(huán)境專項(xiàng)辦公室，北京市，100081 3 91208部隊(duì)，青島市，266000

地磁變化場(chǎng)是地磁學(xué)和空間物理學(xué)的研究重點(diǎn)[1]。疊加在地球穩(wěn)定磁場(chǎng)上的變化磁場(chǎng)可以分為2類：1)隨時(shí)間緩慢變化的長(zhǎng)期變化地磁場(chǎng)，其變化周期一般為幾十年甚至更長(zhǎng)；2)受到固體地球外部各種空間電流體系影響產(chǎn)生的短期變化地磁場(chǎng)[2]。短期變化地磁場(chǎng)又可分為平靜變化和擾動(dòng)變化2類：地磁平靜變化比較平緩，有一定的周期性，分為太陽(yáng)靜日變化和太陰日變化；擾動(dòng)變化的發(fā)生較為偶然，持續(xù)一段時(shí)間后消失，變化短暫且不規(guī)律[2]。

變化磁場(chǎng)在地球總磁場(chǎng)中占比僅1%左右，但卻對(duì)地磁導(dǎo)航影響巨大[3]。通常情況下，地磁平均日變幅可達(dá)幾十nT，因此對(duì)于地磁導(dǎo)航而言，地磁圖除了需要高精度的穩(wěn)定磁場(chǎng)外，還需要高時(shí)空分辨率的變化磁場(chǎng)，而現(xiàn)有的地磁變化場(chǎng)模型的時(shí)間分辨率不足以滿足地磁導(dǎo)航的需求。針對(duì)變化磁場(chǎng)，傳統(tǒng)的建模方法是根據(jù)其物理機(jī)制、觀測(cè)數(shù)據(jù)、形成猜想以及物理方程建立空間電流體系模型，但其物理過(guò)程極為復(fù)雜、理論化難度較大，無(wú)法達(dá)到相應(yīng)的導(dǎo)航精度[4]。隨著計(jì)算機(jī)技術(shù)和人工智能的發(fā)展，在變化磁場(chǎng)建模的過(guò)程中可以依靠神經(jīng)網(wǎng)絡(luò)等方法避開復(fù)雜的物理機(jī)制，僅在數(shù)據(jù)驅(qū)動(dòng)下開展研究。Sutcliffe[5]首次提出將人工智能應(yīng)用于地磁日變化模型的開發(fā)；Yi等[6]建立地磁變化場(chǎng)多時(shí)間尺度組合模型，該方法結(jié)合可變地磁場(chǎng)的物理特征，融合多個(gè)模型，預(yù)測(cè)時(shí)間跨度較長(zhǎng)；Moghadam等[7]提出一種ENN神經(jīng)網(wǎng)絡(luò)算法，并對(duì)Kp和Dst等地磁活動(dòng)指數(shù)進(jìn)行預(yù)測(cè)；程文凱等[8]提出基于XGBoost的機(jī)器學(xué)習(xí)方法重構(gòu)地磁日變數(shù)據(jù)，相較于反向傳播神經(jīng)網(wǎng)絡(luò)，其精度更高。

本文提出融合混沌理論和RBF神經(jīng)網(wǎng)絡(luò)的預(yù)測(cè)模型，在變化磁場(chǎng)的重構(gòu)相空間中，利用RBF神經(jīng)網(wǎng)絡(luò)逼近相點(diǎn)的演化規(guī)律，對(duì)變化磁場(chǎng)進(jìn)行預(yù)測(cè)。

1 理論與方法

1.1 混沌理論

混沌是非線性動(dòng)力系統(tǒng)的固有特性，廣泛存在于自然系統(tǒng)和社會(huì)系統(tǒng)中[9]?；煦缦到y(tǒng)是指在一個(gè)確定的系統(tǒng)里，存在看似隨機(jī)發(fā)展的不規(guī)則運(yùn)動(dòng)，這種運(yùn)動(dòng)無(wú)周期性、不收斂，具有不確定性和不可重復(fù)性，并且對(duì)于初始值極為敏感。混沌理論是非線性系統(tǒng)在一定的參數(shù)條件下，發(fā)生分叉、周期與非周期運(yùn)動(dòng)相互糾纏而產(chǎn)生的某種非周期有序運(yùn)動(dòng)的理論。

相空間重構(gòu)是混沌時(shí)間序列預(yù)測(cè)的重要基礎(chǔ)理論，可以從單一時(shí)間序列中提取和恢復(fù)出原系統(tǒng)的規(guī)律[10]。假設(shè)變化磁場(chǎng)時(shí)間序列為{z}j,j=1,2,…,n，用時(shí)間延遲的方法構(gòu)造M=N-(m-1)τ個(gè)m維相空間矢量：

Z(ti)=(zi,zi+τ,…,zi+(m-1)τ)∈Rn

(1)

重構(gòu)軌道為：

Z=[Z1,Z2,…,Zm]T

(2)

式中，m和τ分別為嵌入維數(shù)和延遲時(shí)間，Z為M×m維矩陣，Zj由m個(gè)分量組成，n個(gè)相點(diǎn)在m維相空間中構(gòu)成一個(gè)相型，各相點(diǎn)連線即為描述動(dòng)力系統(tǒng)在m維相空間里的演化軌跡[10]。在相空間重構(gòu)過(guò)程中，合適的m和τ對(duì)于地磁變化時(shí)間序列中的信息挖掘具有重要意義。

地磁變化場(chǎng)起源于磁層和電流層體系，是一個(gè)極為復(fù)雜的系統(tǒng)，難以用理論模型進(jìn)行精準(zhǔn)建模。但變化磁場(chǎng)數(shù)據(jù)是一個(gè)非線性、非平穩(wěn)的時(shí)間序列，具有混沌特性[11]，因此本文引入混沌理論來(lái)挖掘數(shù)據(jù)深層特征，利用數(shù)據(jù)驅(qū)動(dòng)尋找隱藏規(guī)律，從而完成預(yù)測(cè)。

1.2 RBF神經(jīng)網(wǎng)絡(luò)

RBF神經(jīng)網(wǎng)絡(luò)是一種使用徑向基函數(shù)作為激活函數(shù)的人工神經(jīng)網(wǎng)絡(luò)，其輸出是輸入的徑向基函數(shù)和神經(jīng)元參數(shù)的線性組合。RBF神經(jīng)網(wǎng)絡(luò)由3層結(jié)構(gòu)組成，網(wǎng)絡(luò)拓?fù)浣Y(jié)構(gòu)見圖1。第1層是輸入層，其節(jié)點(diǎn)個(gè)數(shù)等于輸入向量維數(shù)；第2層是隱含層，該層是網(wǎng)絡(luò)的核心部分，決定網(wǎng)絡(luò)性能，其變換函數(shù)是局部分布的非負(fù)非線性函數(shù)，對(duì)中心點(diǎn)徑向?qū)ΨQ且逐漸減弱；第3層是輸出層，該層是隱含層輸出的線性加權(quán)求和[10]。RBF神經(jīng)網(wǎng)絡(luò)有較強(qiáng)的線性擬合能力，收斂速度較快，有較好的泛化能力，可用于構(gòu)建變化磁場(chǎng)預(yù)測(cè)模型。

圖1 RBF神經(jīng)網(wǎng)絡(luò)拓?fù)浣Y(jié)構(gòu)Fig.1 RBF neural network topology

2 地磁變化預(yù)測(cè)模型構(gòu)建

2.1 地磁變化場(chǎng)混沌特性檢驗(yàn)

首先使用李雅普諾夫特性指數(shù)(Lyapunov characteristic exponents，LCE)判斷地磁變化時(shí)間序列是否具有混沌特性。

假設(shè)某一段時(shí)間內(nèi)采集的地磁變化場(chǎng)時(shí)間序列為{zj},j=1,2,…,n，具體判斷步驟如下：

1)通過(guò)FFT變換計(jì)算地磁變化時(shí)間序列的延遲時(shí)間τ和平均周期p。

2)根據(jù)m≥2D+1的原則確定嵌入維數(shù)m，根據(jù)τ和m重構(gòu)相空間：

Z(ti)=(zi,zi+τ,…,zi+(m-1)τ)∈Rn

(3)

式中，i=1,2,…,M，N=M+(m+1)τ。將新產(chǎn)生的序列定義為{Zj},j=1,2,…,M。

(4)

4)計(jì)算中心點(diǎn)Z(t)的dt(i)：

(5)

5)計(jì)算平均z(i)：

(6)

式中，q為dt(i)的數(shù)量。由最小二乘法計(jì)算式(5)中dt(i)的偏微分，從而求得上式斜率，即最大李雅普諾夫指數(shù)λ[12]。當(dāng)λ<0時(shí)，系統(tǒng)有穩(wěn)定的不動(dòng)點(diǎn)；當(dāng)λ=0時(shí)，對(duì)應(yīng)分叉點(diǎn)或系統(tǒng)的周期解，即系統(tǒng)出現(xiàn)周期現(xiàn)象；當(dāng)λ>0時(shí)，系統(tǒng)具有混沌特性。

2.2 偽最近鄰域法確定嵌入維數(shù)

混沌時(shí)間序列是高維相空間混沌運(yùn)動(dòng)的軌跡在低維空間上的投影，在投影過(guò)程中，混沌軌跡會(huì)發(fā)生一定程度的扭曲，可能導(dǎo)致高維空間中不相鄰的2個(gè)相點(diǎn)投影到低維空間后成為相鄰相點(diǎn)，即偽最近鄰點(diǎn)，導(dǎo)致混沌時(shí)間序列無(wú)規(guī)律[13]。偽最近鄰域法確定嵌入維數(shù)的方法如下：

若Rm+1(t)與Rm(t)差別較大，則可以認(rèn)為是高維混沌吸引子2個(gè)點(diǎn)在投影到低維空間上時(shí)變?yōu)閭巫罱忺c(diǎn)。令

(7)

若Sm>ST，則判定ZF(t)是Z(t)的虛假最近鄰點(diǎn)。ST為閾值，可在[10, 50]之間選擇。

從最小嵌入維數(shù)m0開始計(jì)算偽最近鄰點(diǎn)比值，逐漸增大嵌入維數(shù)m，當(dāng)m增大到偽最近鄰點(diǎn)比值小于5%或偽最近鄰點(diǎn)不隨嵌入維數(shù)的增加而減少時(shí)，可以認(rèn)為奇異吸引子已經(jīng)完全展開，此時(shí)的m即為最佳嵌入維數(shù)[13]。

2.3 自相關(guān)函數(shù)法確定延遲時(shí)間

對(duì)于地磁變化序列{zj},j=1,2,…,n，其時(shí)間跨度為jτ的自相關(guān)函數(shù)定義為：

(8)

固定j，繪制自相關(guān)函數(shù)關(guān)于時(shí)間τ的函數(shù)圖像，當(dāng)自相關(guān)函數(shù)下降到初始值的1-1/e倍時(shí)，所得時(shí)間τ即為延遲時(shí)間[13]。

2.4 混沌RBF神經(jīng)網(wǎng)絡(luò)組合預(yù)測(cè)

在RBF神經(jīng)網(wǎng)絡(luò)的基礎(chǔ)上，融合混沌理論對(duì)樣本集進(jìn)行優(yōu)化，得到更適用于預(yù)測(cè)地磁變化的混沌RBF神經(jīng)網(wǎng)絡(luò)模型，具體步驟如下。

1)將已知的地磁變化垂向分量數(shù)據(jù)作為預(yù)測(cè)的初始時(shí)間序列，根據(jù)混沌理論中的偽最近鄰域法和自相關(guān)函數(shù)法計(jì)算最優(yōu)嵌入維數(shù)m和延遲時(shí)間τ，重構(gòu)相空間，得到經(jīng)混沌參數(shù)優(yōu)化的數(shù)據(jù)集：

Z=[Z1,Z2,…,Zm]T

(9)

Z(ti)=(zi,zi+τ,…,zi+(m-1)τ)∈Rn

(10)

2)將經(jīng)相空間重構(gòu)后的地磁變化垂向分量數(shù)據(jù)作為RBF神經(jīng)網(wǎng)絡(luò)輸入數(shù)據(jù)，一部分?jǐn)?shù)據(jù)作為訓(xùn)練樣本集，另一部分?jǐn)?shù)據(jù)作為測(cè)試樣本集。

3)利用優(yōu)化后的訓(xùn)練樣本對(duì)神經(jīng)網(wǎng)絡(luò)進(jìn)行訓(xùn)練，并通過(guò)調(diào)試參數(shù)構(gòu)建合理網(wǎng)絡(luò)，最終通過(guò)測(cè)試樣本得到地磁變化垂向分量的預(yù)測(cè)數(shù)據(jù)。

3 實(shí)驗(yàn)驗(yàn)證

使用5個(gè)臺(tái)站的數(shù)據(jù)進(jìn)行實(shí)驗(yàn)驗(yàn)證，分別是北京十三陵(SSL)站、四川郫縣(PXT)站、黑龍江漠河(MHT)站、西藏拉薩(LAT)站、廣州肇慶(ZQT)站，所選臺(tái)站分布較分散，可以代表我國(guó)大部分區(qū)域的地磁變化特征(圖2)。使用磁通門磁力儀測(cè)得秒采樣數(shù)據(jù)，時(shí)間跨度為2020-02-12～03-09。在實(shí)際應(yīng)用中，匹配分量的時(shí)間分辨率越高越好，但是過(guò)高的分辨率會(huì)導(dǎo)致計(jì)算量劇增、預(yù)測(cè)時(shí)間縮短。目前常用的地磁匹配分量為地磁垂向分量Z和磁場(chǎng)總強(qiáng)度F[14]。綜合考慮后，采用Z的30 min均值作為實(shí)驗(yàn)數(shù)據(jù)，并以樣本集前20 d的數(shù)據(jù)訓(xùn)練預(yù)測(cè)模型，選取1周的測(cè)試數(shù)據(jù)作為預(yù)測(cè)檢驗(yàn)樣本集。

圖2 臺(tái)站分布Fig.2 Distribution of the stations

3.1 SSL站建模預(yù)測(cè)

由于混沌動(dòng)力系統(tǒng)中噪聲信號(hào)的不可預(yù)測(cè)性和高破壞性會(huì)掩蓋系統(tǒng)的內(nèi)在動(dòng)力學(xué)特性，影響后續(xù)的預(yù)測(cè)精度以及對(duì)系統(tǒng)的準(zhǔn)確描述，因此首先要對(duì)地磁數(shù)據(jù)進(jìn)行去噪和平滑處理[13]。采用較為簡(jiǎn)單的拉依達(dá)準(zhǔn)則剔除粗大誤差，以SSL站為例進(jìn)行建模及預(yù)測(cè)。

首先使用MATLAB計(jì)算得到地磁變化Z分量數(shù)據(jù)的最大李雅普諾夫指數(shù)λ=0.35，可知地磁變化Z分量具有混沌特性。

然后采用自相關(guān)函數(shù)法和偽最近鄰域法確定延遲時(shí)間τ和嵌入維數(shù)m，結(jié)果見圖3、4。

圖3 自相關(guān)函數(shù)法計(jì)算時(shí)間延遲τFig.3 Calculation of time delay τ by autocorrelation function method

圖3中紅色直線為Rxx=1-1/e，由圖可知，當(dāng)自相關(guān)函數(shù)Rxx下降到初始值的1-1/e倍時(shí)，Z分量的延遲時(shí)間τ=4 s。由圖4可見，偽鄰近率隨m的增大逐漸減小，并在m=14時(shí)成為常數(shù)值，不再隨m發(fā)生變化，說(shuō)明此時(shí)m已經(jīng)能夠還原混沌吸引子的結(jié)構(gòu)，m=14即為模型的最佳嵌入維數(shù)。當(dāng)RBF神經(jīng)網(wǎng)絡(luò)每層神經(jīng)元數(shù)目與混沌時(shí)間序列重構(gòu)相空間的嵌入維數(shù)m相等時(shí)，預(yù)測(cè)效果最佳[15]。因此，本文選取的神經(jīng)元數(shù)目為14。

圖4 偽最近鄰域法計(jì)算嵌入維數(shù)mFig.4 Calculation of embedding dimension m by false nearest neighbor method

將經(jīng)過(guò)平滑和歸一化處理后的研究時(shí)段前20 d共1 296個(gè)Z分量數(shù)據(jù)的前960個(gè)作為訓(xùn)練樣本集，經(jīng)相空間重構(gòu)后分為訓(xùn)練輸入樣本集和訓(xùn)練期望輸出樣本集；將后336個(gè)樣本作為需要預(yù)測(cè)的未知樣本，經(jīng)相空間重構(gòu)后分為測(cè)試輸入樣本集和測(cè)試期望輸出樣本集[10]。選取訓(xùn)練輸入樣本集中最后(m-1)τ+1個(gè)數(shù)據(jù)作為預(yù)測(cè)輸入樣本集進(jìn)行遞推預(yù)測(cè)，將預(yù)測(cè)結(jié)果與傳統(tǒng)RBF神經(jīng)網(wǎng)絡(luò)算法的預(yù)測(cè)結(jié)果進(jìn)行比較(圖5)。其中，測(cè)試輸入樣本集不參與預(yù)測(cè)，期望輸出樣本集作為模型預(yù)測(cè)結(jié)果的參照值[10]。

圖5 算法改進(jìn)前后預(yù)測(cè)結(jié)果及誤差絕對(duì)值Fig.5 Prediction results and absolute values of errors before and after algorithm improvement

由圖5可見，2種模型均可以跟蹤預(yù)測(cè)Z分量的趨勢(shì)變化，但與傳統(tǒng)RBF神經(jīng)網(wǎng)絡(luò)算法相比，混沌RBF神經(jīng)網(wǎng)絡(luò)模型的預(yù)測(cè)結(jié)果明顯更接近于實(shí)際觀測(cè)值，這種優(yōu)勢(shì)在地磁平緩變化時(shí)更加顯著。在地磁場(chǎng)發(fā)生較大變化時(shí)，預(yù)測(cè)擬合程度是評(píng)價(jià)模型性能的重要指標(biāo)?？梢园l(fā)現(xiàn)，相比于RBF神經(jīng)網(wǎng)絡(luò)算法，混沌RBF神經(jīng)網(wǎng)絡(luò)模型對(duì)曲線突變的反映效果更好，但是誤差仍偏大，這可能是因?yàn)槭艿教?yáng)活動(dòng)等因素的影響發(fā)生地磁擾動(dòng)時(shí)，地磁變化場(chǎng)復(fù)雜性增大、可預(yù)測(cè)性降低。從圖5還可以看出，RBF神經(jīng)網(wǎng)絡(luò)模型誤差的絕對(duì)值大約在8 nT以內(nèi)，且離散程度較大；本文改進(jìn)模型誤差的絕對(duì)值大約在4 nT以內(nèi)，且分布較為集中，精度提高約50%。

3.2 模型泛用性評(píng)估

為了評(píng)估改進(jìn)模型在我國(guó)其他區(qū)域的適用性，對(duì)5個(gè)臺(tái)站進(jìn)行建模預(yù)測(cè)，并計(jì)算模型結(jié)果的最大誤差、平均絕對(duì)誤差MAD和均方根誤差RMSE，結(jié)果見表1。

表1 各臺(tái)站誤差結(jié)果Tab.1 Error results of each station

由表1可見，預(yù)測(cè)結(jié)果的MAD均小于4 nT，RMSE均小于5 nT，誤差較小，說(shuō)明基于混沌RBF神經(jīng)網(wǎng)絡(luò)的地磁變化預(yù)測(cè)模型能較好地描述變化磁場(chǎng)的特性。模型對(duì)MHT和SSL站的預(yù)測(cè)效果較好，MAD小于2 nT，RMSE小于3 nT，但對(duì)ZQT和LAT站的預(yù)測(cè)結(jié)果誤差偏大。

為了更準(zhǔn)確地研究模型的預(yù)測(cè)精度，計(jì)算各臺(tái)站誤差范圍所占比例，結(jié)果見表2(單位%)。

表2 絕對(duì)誤差統(tǒng)計(jì)Tab.2 Absolute error statistics

由表2可見，所有臺(tái)站的預(yù)測(cè)結(jié)果絕對(duì)誤差小于6 nT的占比均超過(guò)80%，小于4 nT的占比均超過(guò)70%，除LAT站以外，其他臺(tái)站小于2 nT 的占比均超過(guò)55%；MHT站的預(yù)測(cè)效果最好，絕對(duì)誤差小于4 nT的占比超過(guò)90%。綜合分析可知，本文模型對(duì)中高緯度地區(qū)的預(yù)測(cè)精度略高于低緯度地區(qū)。

4 結(jié) 語(yǔ)

1)混沌RBF神經(jīng)網(wǎng)絡(luò)地磁變化單站預(yù)測(cè)模型能夠有效預(yù)測(cè)地磁變化垂向分量Z的變化趨勢(shì)，相較于傳統(tǒng)的預(yù)測(cè)模型，該模型的預(yù)測(cè)效果更好，精度提高約50%。

2)模型預(yù)測(cè)精度對(duì)緯度變化較為敏感，中高緯度地區(qū)精度較高，MAD不超過(guò)2 nT；低緯度地區(qū)精度相對(duì)較低，MAD約為3 nT。