盧兆興 呂志峰 李 婷 張金生 姚 垚

1 中國(guó)人民解放軍75838部隊(duì)，廣州市，510515 2 火箭軍士官學(xué)校，山東省青州市，262500 3 火箭軍工程大學(xué)，西安市，710025

地磁場(chǎng)是一種隨時(shí)間和空間變化的矢量場(chǎng)，在多個(gè)領(lǐng)域均有廣泛應(yīng)用，如空間天氣監(jiān)測(cè)和預(yù)報(bào)[1]、地震預(yù)報(bào)[2]、地磁導(dǎo)航與地磁尋的制導(dǎo)技術(shù)[3]等。地磁場(chǎng)的時(shí)空變化對(duì)于地磁場(chǎng)應(yīng)用的精度和范圍具有重要影響，因此對(duì)地磁變化場(chǎng)進(jìn)行建模和預(yù)測(cè)在科學(xué)研究中具有重大意義[4]。現(xiàn)有資料表明[5-8]，地磁場(chǎng)的成因機(jī)理極其復(fù)雜，很難通過(guò)機(jī)理建模的方法對(duì)其時(shí)空變化規(guī)律進(jìn)行準(zhǔn)確描述；部分針對(duì)地磁變化場(chǎng)的高精度預(yù)測(cè)方法，由于受地磁數(shù)據(jù)量限制，缺乏時(shí)間和空間上的廣泛適用性。

目前，全球多個(gè)國(guó)家已建立地磁觀測(cè)臺(tái)站，以研究地磁場(chǎng)的時(shí)空變化，臺(tái)站總數(shù)已達(dá)百余個(gè)。這些地磁臺(tái)站分布在全球各地，24 h不間斷測(cè)量，采樣頻率達(dá)1 kHz，觀測(cè)數(shù)據(jù)量極大，其中必然蘊(yùn)含著地球磁場(chǎng)變化在時(shí)間和空間上的分布規(guī)律。本文以現(xiàn)有地磁臺(tái)站的觀測(cè)數(shù)據(jù)為基礎(chǔ)，分析地磁變化場(chǎng)數(shù)據(jù)在時(shí)空上的關(guān)聯(lián)性，研究數(shù)據(jù)中蘊(yùn)含的規(guī)律信息，并利用BP神經(jīng)網(wǎng)絡(luò)建立地磁變化場(chǎng)預(yù)測(cè)模型，以實(shí)現(xiàn)地磁變化場(chǎng)的時(shí)空預(yù)測(cè)。

1 地磁變化場(chǎng)時(shí)空關(guān)聯(lián)性分析

為說(shuō)明地磁變化場(chǎng)在空間和時(shí)間上的關(guān)聯(lián)性，以長(zhǎng)春、烏魯木齊、成都和蘭州4個(gè)臺(tái)站的數(shù)據(jù)為例，分析相同緯度不同經(jīng)度、相同經(jīng)度不同緯度、不同月同日及同月不同日的地磁數(shù)據(jù)，并通過(guò)定性(作圖觀察)和定量(計(jì)算相關(guān)系數(shù))方法探討其在時(shí)空上的關(guān)聯(lián)性。

采用相關(guān)分析法對(duì)關(guān)聯(lián)性進(jìn)行定量計(jì)算。相關(guān)分析法是研究2個(gè)或2個(gè)以上變量之間的相關(guān)程度并用一定函數(shù)來(lái)表達(dá)相互關(guān)系的方法，相關(guān)變量之間不存在確定性關(guān)系，通常用相關(guān)系數(shù)R表述2個(gè)變量的密切程度，R絕對(duì)值越接近1表明相關(guān)性越好，其計(jì)算公式為：

1.1 地磁變化場(chǎng)時(shí)均值在空間上的相關(guān)性分析

為研究地磁變化場(chǎng)在相同緯度不同經(jīng)度的相關(guān)性，選取長(zhǎng)春臺(tái)(44.0°N、125.2°E)及烏魯木齊臺(tái)(43.8°N、87.8°E)數(shù)據(jù)，對(duì)其進(jìn)行由宏觀特性到個(gè)例的分析。圖1為2012年長(zhǎng)春臺(tái)和烏魯木齊臺(tái)全年的地磁變化場(chǎng)時(shí)均值，由圖可知：1)從宏觀上看，2個(gè)臺(tái)站的地磁變化場(chǎng)時(shí)均值曲線的形態(tài)非常相似；2)2個(gè)臺(tái)站尖峰值的出現(xiàn)均是由于磁暴，且出現(xiàn)尖峰值的時(shí)間幾乎一致；3)2個(gè)臺(tái)站數(shù)據(jù)的相關(guān)系數(shù)R=0.710 8。綜上可知，相同緯度不同經(jīng)度的地磁變化場(chǎng)具有較強(qiáng)的相關(guān)性。為更加精細(xì)地體現(xiàn)其相關(guān)性，任意抽取4 d的數(shù)據(jù)進(jìn)行分析。圖2為2012-01-01～01-04長(zhǎng)春臺(tái)與烏魯木齊臺(tái)地磁變化場(chǎng)時(shí)均值變化曲線，經(jīng)計(jì)算可知，2個(gè)臺(tái)站每日的地磁變化場(chǎng)相關(guān)系數(shù)分別為0.871 5、0.903 3、0.924 8和0.963 3，相關(guān)系數(shù)較大，說(shuō)明相同緯度不同經(jīng)度的地磁變化場(chǎng)具有較強(qiáng)的相關(guān)性。

圖1 2012年長(zhǎng)春臺(tái)和烏魯木齊臺(tái)全年的地磁變化場(chǎng)時(shí)均值對(duì)比Fig.1 Hourly mean variable geomagnetic field during 2012 in Changchun station and Urumqi station

圖2 2012-01-01～01-04長(zhǎng)春臺(tái)和烏魯木齊臺(tái)地磁變化場(chǎng)時(shí)均值曲線Fig.2 Hourly mean variable geomagnetic field from January 1, 2012 to January 4, 2012 in Changchun station and Urumqi station

為研究地磁變化場(chǎng)在相同經(jīng)度不同緯度的相關(guān)性，選取成都臺(tái)(103.7°E、31.0°N)及蘭州臺(tái)(103.8°E、36.1°N)數(shù)據(jù)進(jìn)行分析。通過(guò)計(jì)算可知，2個(gè)臺(tái)站2012年全年數(shù)據(jù)的相關(guān)系數(shù)R=0.999 8，2012-01-01～01-04每日的地磁變化場(chǎng)相關(guān)系數(shù)分別為0.945 5、0.926 6、0.975 7和0.928 6，相關(guān)系數(shù)極大，說(shuō)明相同經(jīng)度不同緯度的地磁變化場(chǎng)也具有較強(qiáng)的相關(guān)性。

1.2 地磁變化場(chǎng)時(shí)均值在時(shí)間上的相關(guān)性分析

選用長(zhǎng)春臺(tái)2012-01-01、04-01、07-01和10-01的地磁變化場(chǎng)時(shí)均值數(shù)據(jù)進(jìn)行分析，結(jié)果見圖3。從圖3可以看出，不同月同日的地磁變化場(chǎng)時(shí)均值的變化趨勢(shì)一致，僅幅值大小存在差別，這些差別在很大程度上由不同月份所造成。假設(shè)1月份為1，4月份為2，7月份為3，10月份為4，Rij為i與j的相關(guān)系數(shù)，經(jīng)計(jì)算可得，R12=0.885 9，R13=0.699 0，R14=0.901 6，R23=0.771 5，R24=0.882 0，R34=0.828 8，說(shuō)明不同月同日的地磁變化場(chǎng)具有較強(qiáng)的相關(guān)性。

圖3 長(zhǎng)春臺(tái)2012-01-01、04-01、07-01、10-01地磁變化場(chǎng)時(shí)均值Fig.3 Hourly mean variable geomagnetic field on January 1, April 1, July 1, October 1, 2012 in Changchun station

選取長(zhǎng)春臺(tái)2012-01-01～01-04的地磁變化場(chǎng)時(shí)均值數(shù)據(jù)進(jìn)行分析，結(jié)果見圖4。假設(shè)第1日為1，第2日為2，第3日為3，第4日為4，經(jīng)計(jì)算可得，R12=0.921 9，R13=0.825 6，R14=0.885 8，R23=0.960 8，R24=0.923 0，R34=0.844 9，說(shuō)明同一個(gè)區(qū)域同月不同日的地磁變化場(chǎng)也具有較強(qiáng)的相關(guān)性。

圖4 長(zhǎng)春臺(tái)2012-01-01～01-04地磁變化場(chǎng)時(shí)均值Fig.4 Hourly mean variable geomagnetic field from January 1, 2012 to January 4, 2012 in Changchun station

綜合分析可知，地磁變化場(chǎng)的觀測(cè)數(shù)據(jù)包含地磁變化場(chǎng)在空間和時(shí)間上的變化規(guī)律，且地磁變化場(chǎng)與其相對(duì)應(yīng)的地理經(jīng)度、緯度及時(shí)間等要素之間具有復(fù)雜的非線性函數(shù)關(guān)系，屬于高維問(wèn)題。因此可將空間因素和時(shí)間因素作為神經(jīng)網(wǎng)絡(luò)的輸入，通過(guò)數(shù)據(jù)驅(qū)動(dòng)固化網(wǎng)絡(luò)參數(shù)，進(jìn)而表達(dá)地磁變化場(chǎng)在時(shí)間和空間上的復(fù)雜非線性映射關(guān)系。

2 基于神經(jīng)網(wǎng)絡(luò)的地磁變化場(chǎng)預(yù)測(cè)

BP神經(jīng)網(wǎng)絡(luò)也稱反向傳播前饋性神經(jīng)網(wǎng)絡(luò)，具有強(qiáng)大的計(jì)算能力，可表達(dá)各種復(fù)雜映射，適用于反映地磁場(chǎng)的復(fù)雜非線性變化規(guī)律[9]，且在地磁空間環(huán)境的預(yù)測(cè)預(yù)報(bào)方面具有較好的應(yīng)用[10]。本文選用BP神經(jīng)網(wǎng)絡(luò)對(duì)地磁變化場(chǎng)進(jìn)行時(shí)空預(yù)測(cè)，確定神經(jīng)網(wǎng)絡(luò)的輸入量和參考輸出，實(shí)現(xiàn)未知區(qū)域(無(wú)地磁臺(tái)站)地磁變化場(chǎng)的預(yù)測(cè)。

2.1 神經(jīng)網(wǎng)絡(luò)結(jié)構(gòu)

神經(jīng)網(wǎng)絡(luò)中隱含層神經(jīng)元的傳遞函數(shù)設(shè)定為正切S型傳遞函數(shù)tansig，輸出層神經(jīng)元的傳遞函數(shù)為線性傳遞函數(shù)purelin，訓(xùn)練函數(shù)設(shè)定為L(zhǎng)evenberg-Marquardt算法訓(xùn)練函數(shù)trainlm。神經(jīng)網(wǎng)絡(luò)參數(shù)設(shè)定后，可得到區(qū)域地磁變化場(chǎng)時(shí)空預(yù)測(cè)的BP神經(jīng)網(wǎng)絡(luò)結(jié)構(gòu)示意圖(圖5)，其中所建網(wǎng)絡(luò)的輸入量個(gè)數(shù)為i，中間隱含層的節(jié)點(diǎn)數(shù)為m，輸出量個(gè)數(shù)為n。

圖5 BP神經(jīng)網(wǎng)絡(luò)示意圖Fig.5 Schematic diagram of BP neural network

由于獲取的數(shù)據(jù)信息量過(guò)大，因此首先需要分析并確定神經(jīng)網(wǎng)絡(luò)的輸入和輸出，然后將輸入輸出作為訓(xùn)練數(shù)據(jù)輸入到神經(jīng)網(wǎng)絡(luò)，通過(guò)訓(xùn)練函數(shù)trainlm對(duì)神經(jīng)網(wǎng)絡(luò)各層之間的權(quán)系數(shù)進(jìn)行調(diào)整。訓(xùn)練過(guò)程中權(quán)系數(shù)w可反映輸入輸出之間的映射關(guān)系，訓(xùn)練結(jié)束后可得到符合精度要求的神經(jīng)網(wǎng)絡(luò)結(jié)構(gòu)，從而實(shí)現(xiàn)對(duì)地磁變化場(chǎng)進(jìn)行時(shí)空預(yù)測(cè)。具體流程如圖6所示。

圖6 基于神經(jīng)網(wǎng)絡(luò)的地磁變化場(chǎng)預(yù)測(cè)示意圖Fig.6 Schematic diagram of variable geomagnetic field forecasting based on neural network

2.2 神經(jīng)網(wǎng)絡(luò)輸入輸出

通過(guò)分析可知，地磁變化場(chǎng)在時(shí)間和空間上具有較好的相關(guān)性，因此可將時(shí)間因素和空間因素作為神經(jīng)網(wǎng)絡(luò)的輸入。為了便于表達(dá)和計(jì)算，將1 a中每日及1 d中每小時(shí)進(jìn)行編號(hào)，小時(shí)編號(hào)用t1表示，每日編號(hào)用t2表示，1 d中00:00～23:00的編號(hào)t1分別為1、2、3、…、24，1 a中01-01～12-31的編號(hào)t2分別為1、2、3、…、365(閏年為366)。盡管t1和t2具有不同取值，但不能將這些不同值直接作為神經(jīng)網(wǎng)絡(luò)的輸入，因?yàn)楫?dāng)日23:00與次日00:00非常接近，當(dāng)年12-31與次年01-01也非常接近，而神經(jīng)網(wǎng)絡(luò)并不能對(duì)其進(jìn)行區(qū)分。因此，將24個(gè)世界時(shí)向圓周投影，并按照式(1)將地方時(shí)t1投影為x1和x2，將日期t2投影為x3和x4，從而完成時(shí)間的圓周投影：

(1)

神經(jīng)網(wǎng)絡(luò)輸入的空間信息主要用經(jīng)度λ和緯度φ來(lái)表示，但λ和φ實(shí)際對(duì)應(yīng)一個(gè)圓心角，對(duì)其進(jìn)行圓周投影就相當(dāng)于將球面位置坐標(biāo)放到平面上，其方法與世界時(shí)的圓周化相同，具體投影公式為：

(2)

2.3 神經(jīng)網(wǎng)絡(luò)訓(xùn)練、驗(yàn)證與預(yù)測(cè)

確定神經(jīng)網(wǎng)絡(luò)的輸入量為8，輸出量為1。為更直觀地說(shuō)明輸入輸出數(shù)據(jù)的形式，統(tǒng)一數(shù)據(jù)的轉(zhuǎn)化順序?yàn)榻?jīng)度、緯度、日、時(shí)，假設(shè)位置為120°E、30°N，時(shí)間為03-01 06:00，變化磁場(chǎng)為20 nT，經(jīng)過(guò)式(1)和式(2)轉(zhuǎn)換，得到輸入量和輸出量分別為X=[-0.5，0.866，0.866，0.5，0.512 4，0.858 8，0，1]T，y=20。

將現(xiàn)有觀測(cè)臺(tái)站的可用數(shù)據(jù)按照上述輸入量和輸出量的形式進(jìn)行表示，并將其分為訓(xùn)練數(shù)據(jù)和驗(yàn)證數(shù)據(jù)2類，訓(xùn)練數(shù)據(jù)(占總數(shù)量90%以上)用于對(duì)神經(jīng)網(wǎng)絡(luò)進(jìn)行訓(xùn)練，驗(yàn)證數(shù)據(jù)用于驗(yàn)證訓(xùn)練后神經(jīng)網(wǎng)絡(luò)的正確性。最后將未知區(qū)域的輸入量加入到已經(jīng)訓(xùn)練好的神經(jīng)網(wǎng)絡(luò)中，實(shí)現(xiàn)對(duì)任意未知區(qū)域地磁變化場(chǎng)的預(yù)測(cè)。

3 地磁臺(tái)站數(shù)據(jù)驗(yàn)證

由于獲取的數(shù)據(jù)量有限，無(wú)法對(duì)全球所有臺(tái)站的地磁場(chǎng)數(shù)據(jù)進(jìn)行收集整理，僅以我國(guó)地磁臺(tái)站的數(shù)據(jù)信息為例，本文共收集10個(gè)地磁臺(tái)站2008～2013年共6 a的地磁總場(chǎng)分鐘值數(shù)據(jù)及磁暴數(shù)據(jù)(表1)。

表1 地磁臺(tái)站信息

將觀測(cè)數(shù)據(jù)中的異常值及磁暴數(shù)據(jù)進(jìn)行剔除和修正后，取時(shí)均值用于計(jì)算。選取用于驗(yàn)證神經(jīng)網(wǎng)絡(luò)在時(shí)空域內(nèi)求解地磁變化場(chǎng)有效性的數(shù)據(jù)，本文選取6個(gè)臺(tái)站2010-06-01～06-02的數(shù)據(jù)進(jìn)行分析，驗(yàn)證結(jié)果見圖7。選取用于驗(yàn)證神經(jīng)網(wǎng)絡(luò)仿真計(jì)算精度的數(shù)據(jù)，本文隨機(jī)選取100個(gè)訓(xùn)練數(shù)據(jù)進(jìn)行分析，結(jié)果見圖8。將剩余的數(shù)據(jù)作為訓(xùn)練數(shù)據(jù)，設(shè)定訓(xùn)練函數(shù)為L(zhǎng)evenberg-Marquardt算法訓(xùn)練函數(shù)trainlm，訓(xùn)練100次后對(duì)結(jié)果進(jìn)行驗(yàn)證。

圖7 地磁變化場(chǎng)時(shí)空計(jì)算結(jié)果Fig.7 The temporal and spatial calculation results of the variable geomagnetic field

圖8 地磁變化場(chǎng)時(shí)空計(jì)算精度Fig.8 The temporal and spatial calculation accuracy of the variable geomagnetic field

由圖7可知：1)神經(jīng)網(wǎng)絡(luò)求解的地磁變化場(chǎng)在06-01～06-02的變化情況與實(shí)際地磁變化場(chǎng)基本一致，說(shuō)明其能夠較好地反映地磁變化場(chǎng)在時(shí)間上的變化規(guī)律；2)各個(gè)臺(tái)站06-01～06-02的地磁變化場(chǎng)無(wú)論是幅值還是相位都存在差異，神經(jīng)網(wǎng)絡(luò)計(jì)算對(duì)各臺(tái)站之間的位置差異敏感，求解結(jié)果與實(shí)際變化情況基本一致，說(shuō)明其能夠較好地反映地磁變化場(chǎng)在空間上的分布規(guī)律。

由圖8可知，神經(jīng)網(wǎng)絡(luò)仿真計(jì)算的均方根誤差RMS=4.8 nT，考慮到地磁變化場(chǎng)1 d內(nèi)的波動(dòng)幅度能達(dá)到幾十nT，在一般擾動(dòng)變化時(shí)能達(dá)到幾百nT，且地磁變化場(chǎng)隨時(shí)間和空間的變化規(guī)律非常復(fù)雜，因此本文建立的BP神經(jīng)網(wǎng)絡(luò)模型的計(jì)算精度能夠滿足一般科學(xué)研究對(duì)地磁變化場(chǎng)的精度需求。

4 結(jié) 語(yǔ)

本文綜合分析了影響地磁變化場(chǎng)的時(shí)間和空間因素，基于國(guó)內(nèi)多個(gè)地磁臺(tái)站的觀測(cè)數(shù)據(jù)構(gòu)建用于地磁變化場(chǎng)預(yù)測(cè)的BP神經(jīng)網(wǎng)絡(luò)模型，并驗(yàn)證了模型的有效性，該研究思路和建模方法能夠?yàn)榈卮艌?chǎng)建模和計(jì)算等相關(guān)研究提供借鑒。本文在地磁變化場(chǎng)數(shù)據(jù)處理時(shí)，僅對(duì)磁暴數(shù)據(jù)進(jìn)行修正，但地磁變化場(chǎng)還受太陽(yáng)風(fēng)、空間粒子及其他外部環(huán)境影響，且本文僅基于陸上地磁臺(tái)站的數(shù)據(jù)，對(duì)海上地磁變化場(chǎng)預(yù)測(cè)的適用性有待研究。后期將綜合考慮各類影響因素，并通過(guò)增加各區(qū)域地磁臺(tái)站數(shù)據(jù)等方式對(duì)所建立的神經(jīng)網(wǎng)絡(luò)模型進(jìn)行訓(xùn)練，以提高地磁變化場(chǎng)的預(yù)測(cè)精度和適用范圍。