王子思禹 師立勤 劉四清 鐘秋珍 陳艷紅 閆曉輝 石育榕 何欣燃

1(中國科學(xué)院國家空間科學(xué)中心 北京 100190)

2(中國科學(xué)院大學(xué) 北京 100049)

3(中國科學(xué)院空間環(huán)境態(tài)勢感知技術(shù)重點(diǎn)實(shí)驗(yàn)室 北京 100190)

0 引言

太陽風(fēng)是太陽高層大氣向外流動(dòng)所形成的超聲速等離子體流。當(dāng)由日冕物質(zhì)拋射事件或冕洞產(chǎn)生的太陽風(fēng)到達(dá)地球時(shí)，會(huì)向磁層注入更多能量和粒子，可能引起磁層擾動(dòng)，從而引發(fā)地磁暴[1]。地磁暴會(huì)嚴(yán)重影響衛(wèi)星等飛行器的性能和安全，從而對(duì)整個(gè)空間和地面的技術(shù)系統(tǒng)的正常運(yùn)行造成威脅。因此，對(duì)太陽風(fēng)造成的地磁擾動(dòng)進(jìn)行及時(shí)預(yù)報(bào)具有重要的實(shí)際意義。

隨著機(jī)器學(xué)習(xí)技術(shù)的不斷突破，其自適應(yīng)、自學(xué)習(xí)的強(qiáng)大并行信息處理能力，在很多方面取得了突破性進(jìn)展。特別是在目標(biāo)推薦系統(tǒng)中的應(yīng)用：通過分析用戶歷史行為，提取特征并進(jìn)行相似度計(jì)算，可以解決用戶從大量信息中挑選目標(biāo)信息這樣復(fù)雜且耗時(shí)的問題[2]。同時(shí)，在空間環(huán)境監(jiān)測方面，作為衡量近地空間全球磁擾強(qiáng)度的重要指標(biāo)之一的地磁Kp指數(shù)，已有長達(dá)80年無間斷的數(shù)據(jù)[3]；太陽風(fēng)方面也積累了幾十年的衛(wèi)星監(jiān)測數(shù)據(jù)。因此，將機(jī)器學(xué)習(xí)技術(shù)應(yīng)用于空間環(huán)境監(jiān)測信息處理，可以更有效地進(jìn)行地磁擾動(dòng)預(yù)報(bào)。

目前已有的機(jī)器學(xué)習(xí)地磁Kp指數(shù)預(yù)報(bào)模型有：Rice 大學(xué)Costello 模型[4]、Lund 天文臺(tái)的Boberg 模型[5]、USAF WingKp模型（JHU/APL 模型）[6]、Rice大學(xué)的Bela 模型[7]等，但這些方法都屬于基于人工神經(jīng)網(wǎng)絡(luò)的回歸問題，需要較多的模型參數(shù)實(shí)現(xiàn)地磁Kp指數(shù)的預(yù)報(bào)。同時(shí)，預(yù)報(bào)員也因?yàn)槭芟抻跈C(jī)器學(xué)習(xí)的“黑盒”過程，不能充分結(jié)合自身經(jīng)驗(yàn)判斷進(jìn)行預(yù)報(bào)。本文通過機(jī)器學(xué)習(xí)的聚類算法和相似度算法，建立了太陽風(fēng)推薦模型。利用該模型可以計(jì)算并推薦出與輸入太陽風(fēng)事例相似的多個(gè)歷史太陽風(fēng)事例。通過對(duì)模型推薦的太陽風(fēng)事例的后續(xù)Kp指數(shù)進(jìn)行分析，以及推薦太陽風(fēng)事例與輸入的太陽風(fēng)事例參數(shù)變化的對(duì)照分析，預(yù)報(bào)員可以更加有效地將自身經(jīng)驗(yàn)融入到地磁擾動(dòng)的預(yù)報(bào)當(dāng)中。同時(shí)，根據(jù)歷史相似太陽風(fēng)事例的后續(xù)影響，預(yù)報(bào)員可以更早地對(duì)當(dāng)前太陽風(fēng)事例后續(xù)變化做出判斷和預(yù)報(bào)。

1 數(shù)據(jù)

太陽風(fēng)及行星際數(shù)據(jù)來源于OMNI 網(wǎng)站*https://omniweb.gsfc.nasa.gov/。OMNI網(wǎng)站整合了來自IMP 8、Geotail、Wind 和 Ace 衛(wèi)星的太陽風(fēng)數(shù)據(jù)。本文選取ACE 衛(wèi)星1998－2019年的太陽風(fēng)數(shù)據(jù)，數(shù)據(jù)時(shí)頻為60 s?？紤]到衛(wèi)星太陽風(fēng)數(shù)據(jù)可能存在缺失，對(duì)于缺失數(shù)據(jù)首先進(jìn)行數(shù)據(jù)清洗：將個(gè)別的缺省數(shù)據(jù)近似取值為其前60 s 對(duì)應(yīng)的數(shù)據(jù)，保證了數(shù)據(jù)的完整性。OMNI 網(wǎng)站提供的太陽風(fēng)參數(shù)包括太陽風(fēng)速度v、太陽風(fēng)質(zhì)子密度N、太陽風(fēng)溫度T、GSM 坐標(biāo)下行星際磁場By和Bz分量的值以及行星際磁場強(qiáng)度B。為了便于太陽風(fēng)推薦模型的訓(xùn)練和學(xué)習(xí)，對(duì)太陽風(fēng)參數(shù)數(shù)據(jù)均進(jìn)行了歸一化處理，消除不同量綱對(duì)計(jì)算的影響。

地磁Kp指數(shù)來源于OMNI 網(wǎng)站整合自德國地球科學(xué)研究中心（GFZ）的結(jié)果。根據(jù)機(jī)器學(xué)習(xí)模型輸入要求，將Kp指數(shù)進(jìn)行臨近值預(yù)處理：即Kp指數(shù)30轉(zhuǎn)化為3，3–轉(zhuǎn)化為2.7，3+轉(zhuǎn)化為3.3，以此類推。根據(jù)當(dāng)前空間天氣預(yù)報(bào)業(yè)務(wù)一般規(guī)范，把Kp指數(shù)介于5–和6+（包含5–和6+）定義為小磁暴，把Kp指數(shù)介于7–和90（包含7–和90）定義為大磁暴。

本文將世界時(shí)整點(diǎn)時(shí)刻后3 h 內(nèi)的太陽風(fēng)數(shù)據(jù)作為一個(gè)太陽風(fēng)事例，進(jìn)而將歷史數(shù)據(jù)集劃分成約18 萬個(gè)太陽風(fēng)事例，其中造成小磁暴影響的事例約4800 個(gè)，造成大磁暴影響的事例約450 個(gè)。另外，將時(shí)間上緊隨其后的第一個(gè)Kp指數(shù)定義為與這一太陽風(fēng)事例對(duì)應(yīng)的Kp指數(shù)，太陽風(fēng)事例領(lǐng)先其對(duì)應(yīng)的Kp指數(shù)1～3 h，即預(yù)報(bào)時(shí)間提前量為1～3 h。

2 模型方法

太陽風(fēng)是造成地磁暴的直接原因，因此具有一致性特征的太陽風(fēng)引起的地磁擾動(dòng)也應(yīng)該具有相似的特征。在歷史數(shù)據(jù)足夠多的條件下，可以在歷史數(shù)據(jù)中尋找到與當(dāng)前太陽風(fēng)特征相似的歷史太陽風(fēng)事例。歷史太陽風(fēng)事例引起的地磁變化，就可以作為當(dāng)前Kp指數(shù)預(yù)報(bào)的參考。

本文主要利用機(jī)器學(xué)習(xí)中的相似度算法來構(gòu)建尋找相似特征太陽風(fēng)事例的推薦模型。具體的步驟分為3 步。第1 步，太陽風(fēng)特征參數(shù)選取?？晒┻x取的太陽風(fēng)特征參數(shù)有速度、密度、溫度、行星際磁場、行星際磁場By分量、行星際磁場Bz分量；在相似度計(jì)算中，過多的特征參數(shù)可能會(huì)導(dǎo)致推薦耗時(shí)倍增，且推薦結(jié)果會(huì)被無影響或極低影響的參數(shù)干擾，需要合理選取對(duì)地磁擾動(dòng)起主要作用的參數(shù)作為相似度計(jì)算的特征參數(shù)。第2 步，參數(shù)權(quán)重計(jì)算。在地磁擾動(dòng)中，太陽風(fēng)各種特征參數(shù)的影響作用不同，需要對(duì)所選取的特征參數(shù)根據(jù)其對(duì)地磁擾動(dòng)的影響計(jì)算相應(yīng)的權(quán)重值。第3 步，推薦模型構(gòu)建。本文以選定的太陽風(fēng)特征參數(shù)作為輸入，基于機(jī)器學(xué)習(xí)中的相似度算法，構(gòu)建篩選相似太陽風(fēng)事例的推薦模型。

2.1 太陽風(fēng)特征參數(shù)選取

采用XGBoost 算法模型進(jìn)行太陽風(fēng)特征參數(shù)選取。XGBoost[8]是一個(gè)決策樹訓(xùn)練模型，該模型依據(jù)輸入的特征參數(shù)在迭代建立決策樹中的作用和價(jià)值，即每個(gè)特征在所有樹中作為劃分屬性的次數(shù)判斷特征重要性，繼而通過每個(gè)屬性分割點(diǎn)改進(jìn)性能度量的量來計(jì)算單個(gè)決策樹的重要性，并由節(jié)點(diǎn)負(fù)責(zé)的觀察數(shù)量加權(quán)，最終將一個(gè)屬性在所有提升樹中的結(jié)果進(jìn)行加權(quán)求和后然后平均，得到重要性得分[9]。

以1998－2019年歸一化后的太陽風(fēng)速度v、太陽風(fēng)質(zhì)子密度N、太陽風(fēng)溫度T、行星際磁場強(qiáng)度B、行星際磁場By、Bz分量1 h 的平均值作為輸入，以相應(yīng)的Kp指數(shù)作為標(biāo)簽，通過XGBoost 模型迭代建立決策樹，對(duì)輸入的六個(gè)特征參數(shù)進(jìn)行了評(píng)分。評(píng)分高低即為XGBoost 模型認(rèn)為不同參數(shù)對(duì)地磁擾動(dòng)（Kp指數(shù)）影響的重要程度。評(píng)分結(jié)果如表1。該結(jié)果也比較符合當(dāng)前的研究結(jié)論，即磁暴的主相發(fā)展與行星際磁場的南向分量和正值行星際電場密切相關(guān)[10]。為了避免過多特征參數(shù)造成推薦結(jié)果的混亂，提高模型運(yùn)行效率，經(jīng)過隨機(jī)的部分樣本的測試，選擇對(duì)地磁擾動(dòng)影響最重要的兩個(gè)參數(shù)：太陽風(fēng)風(fēng)速v和行星際磁場Bz分量作為太陽風(fēng)事例推薦模型的輸入?yún)?shù)時(shí)，推薦序列相似性和推薦耗時(shí)情況最好，所以本文最終選擇v和Bz作為相似度計(jì)算輸入的特征。

表1 XGBoost 對(duì)太陽風(fēng)特征參數(shù)評(píng)分結(jié)果Table 1 Scores of feature parameters using XGBoost

2.2 參數(shù)權(quán)重計(jì)算

為了有效地計(jì)算篩選后的特征的權(quán)重，避免多余特征的影響，本文采用單層全連接神經(jīng)網(wǎng)絡(luò)對(duì)篩選后保留的太陽風(fēng)風(fēng)速v，行星際磁場Bz分量進(jìn)行訓(xùn)練，并評(píng)估它們各自對(duì)地磁擾動(dòng)影響的權(quán)重。以歸一化后的太陽風(fēng)風(fēng)速v，行星際磁場Bz分量作為單層神經(jīng)網(wǎng)絡(luò)的輸入層參數(shù)，時(shí)間分辨率是1 h，輸出層為每小時(shí)相應(yīng)的Kp指數(shù)，損失函數(shù)為預(yù)測值和實(shí)際值的均方差函數(shù)，優(yōu)化器采用Adam（Adaptive moment estimation）算法[11]。通過單層全連接神經(jīng)網(wǎng)絡(luò)訓(xùn)練后，最終得到的權(quán)重比為|Wv|∶|WBz|=1.41∶1。

2.3 推薦模型構(gòu)建

本文以機(jī)器學(xué)習(xí)中的相似度算法構(gòu)建太陽風(fēng)推薦模型，主要采用了歐氏距離和動(dòng)態(tài)時(shí)間規(guī)整（Dynamic Time Warping，DTW）兩種算法。歐式距離是指m維空間中兩個(gè)點(diǎn)之間的真實(shí)距離，或者向量的自然長度，在二維和三維空間中的歐氏距離就是兩點(diǎn)之間的實(shí)際距離[12]。設(shè)空間中兩點(diǎn)X、Y的坐標(biāo)為X(X1,X2,···,Xn),Y(Y1,Y2,···,Yn)，則X點(diǎn)與Y點(diǎn)間的歐式距離為

當(dāng)其輸入為兩組特征參量時(shí)，可以反映兩組特征參量之間的綜合差異。

動(dòng)態(tài)時(shí)間規(guī)整算法（Dynamic Time Warping，簡稱DTW 算法）是尋找時(shí)間序列相似的算法之一，其原理是給定兩個(gè)序列Q={Q1,Q2,···,Qi,···,Qn}和C={C1,C2,···,Cj,···,Cm}，其長度分別是n和m，構(gòu)造一個(gè)n×m的矩陣網(wǎng)格，如圖1（a）所示，用紅線和藍(lán)線分別表示兩個(gè)序列，以歐氏距離D為標(biāo)準(zhǔn)，矩陣元素(i,j)為Qi和Cj兩個(gè)點(diǎn)的距離D(Qi,Cj)。從序列起始段所在的矩陣角為邊界條件，在滿足連續(xù)性和單調(diào)性約束的同時(shí)，通過動(dòng)態(tài)規(guī)劃求得距離累積值最小的路徑即為最佳路徑[13]。

如圖1(b)所示，動(dòng)態(tài)時(shí)間規(guī)整算法認(rèn)為，通過對(duì)時(shí)間序列進(jìn)行局部非線性縮放，使兩個(gè)序列形態(tài)盡可能對(duì)齊后的歐氏距離累計(jì)值才是兩個(gè)時(shí)間序列的最小距離。

圖1 動(dòng)態(tài)時(shí)間規(guī)整算法原理Fig.1 Principle of DTW

太陽風(fēng)推薦模型包含初篩和精篩兩個(gè)部分。初篩采用歐式距離方法，輸入為一個(gè)太陽風(fēng)事例中太陽風(fēng)風(fēng)速v與行星際磁場Bz的1 h 平均值，3 h 即為維度3×2 的矩陣。初篩通過逐個(gè)計(jì)算歷史事例的3×2 矩陣與輸入事例對(duì)應(yīng)的3×2 矩陣的賦權(quán)歐式距離進(jìn)行篩選，找到3 個(gè)距離最近的歷史事例，作為推薦事例。精篩采用動(dòng)態(tài)時(shí)間規(guī)整算法，輸入為一個(gè)太陽風(fēng)事例中太陽風(fēng)風(fēng)速v和行星際磁場Bz的60 s值，3 h 即為維度180×2 的矩陣。精篩計(jì)算歷史上每個(gè)造成地磁暴的太陽風(fēng)事例與輸入事例之間的局部非線性縮放的歐氏距離，通過篩選找到最相似的3 個(gè)歷史太陽風(fēng)事例作為推薦的太陽風(fēng)事例。太陽風(fēng)推薦模型實(shí)際運(yùn)行時(shí)，首先對(duì)輸入事例進(jìn)行初篩，若初篩結(jié)果為輸入太陽風(fēng)事例不會(huì)造成地磁暴，則不繼續(xù)進(jìn)行更高精度的推薦，僅選用初篩的結(jié)果作為預(yù)報(bào)參考。若初篩結(jié)果認(rèn)為可能發(fā)生磁暴，則以每分鐘時(shí)間間隔的特征參數(shù)作為輸入，繼續(xù)通過精篩在引發(fā)磁暴的歷史太陽風(fēng)事例中尋找相似事例，得到參數(shù)時(shí)間序列變化相似的歷史太陽風(fēng)事例作為模型精篩的推薦結(jié)果。空間天氣業(yè)務(wù)預(yù)報(bào)中，一般對(duì)可能發(fā)生磁暴的情況更為關(guān)注，因此推薦模型的精篩主要針對(duì)磁暴事例。推薦模型中沒有全部采用精篩方式，也考慮到在大量歷史數(shù)據(jù)中進(jìn)行動(dòng)態(tài)時(shí)間規(guī)整算法計(jì)算時(shí)，用時(shí)可能較長，降低預(yù)報(bào)實(shí)時(shí)性。

太陽風(fēng)事例領(lǐng)先其對(duì)應(yīng)的Kp指數(shù)1～3 h，因此推薦模型的預(yù)報(bào)提前量為1～3 h。

3 模型測試結(jié)果

以1998－2019年的太陽風(fēng)數(shù)據(jù)作為歷史太陽風(fēng)數(shù)據(jù)對(duì)推薦模型進(jìn)行測試。測試中，從1998－2019年歷史太陽風(fēng)事例中隨機(jī)選擇120 個(gè)太陽風(fēng)事例用作測試用例，這些測試事例覆蓋了1998－2019年的全部年份，并且覆蓋了從0～9 全部Kp指數(shù)的等級(jí)，其中沒有造成地磁暴影響的、造成小磁暴影響的和造成大磁暴影響的比例為1∶1∶1，各40 個(gè)。

選取測試用例的最優(yōu)推薦事例（非自身的top1）對(duì)應(yīng)Kp指數(shù)與輸入事例對(duì)應(yīng)的實(shí)際Kp指數(shù)比較，測試結(jié)果如圖2。圖中橫軸是輸入太陽風(fēng)事例對(duì)應(yīng)的實(shí)際地磁Kp指數(shù)，縱軸是推薦的歷史太陽風(fēng)事例的地磁Kp指數(shù)；I 區(qū)是無磁暴影響的太陽風(fēng)，II 區(qū)是造成小磁暴影響的太陽風(fēng)，III 區(qū)是造成大磁暴影響的太陽風(fēng)。由圖可知，推薦太陽風(fēng)事例的地磁Kp指數(shù)與輸入太陽風(fēng)事例的地磁Kp指數(shù)具有較好的一致性。Kp平均絕對(duì)誤差為0.65，Kp指數(shù)推薦值與實(shí)際值的均方根誤差為0.79，相關(guān)系數(shù)為0.93。

圖2 推薦模型測試結(jié)果Fig.2 Test results of the recommanded model

表2 列出了典型的基于人工神經(jīng)網(wǎng)絡(luò)方法的已有Kp指數(shù)預(yù)報(bào)模型與本文模型的比較。對(duì)比結(jié)果表明，本文的推薦模型優(yōu)化了參數(shù)輸入量的選取，Kp預(yù)報(bào)結(jié)果的相關(guān)性更好，預(yù)報(bào)結(jié)果誤差較小。

表2 本文模型與現(xiàn)有典型模型比較Table 2 Comparison between the model proposed in this paper and the existing typical models

本文模型可以推薦出3 個(gè)最相似的事例，測試中在3 個(gè)最相似的事例中總能找出較為理想的推薦結(jié)果（推薦事例與實(shí)際Kp的偏差在1 以內(nèi)）。表3 給出了推薦的三個(gè)最相似的事例的距離與Kp平均絕對(duì)誤差（保留兩位小數(shù)）的統(tǒng)計(jì)結(jié)果。由表3 可知，盡管存在一些距離與Kp指數(shù)最小變化精度不同導(dǎo)致的例外，但距離越相近的事例，Kp指數(shù)預(yù)報(bào)值越可能接近實(shí)測值。

表3 不同Kp 指數(shù)下三個(gè)最相似事例的距離與平均Kp 誤差Table 3 Distances of 3 recommanded cases and Kp error in different Kp indexes

4 討論

由于太陽風(fēng)對(duì)地球磁層的能量輸入方式極其復(fù)雜且高度動(dòng)態(tài)，進(jìn)入地球磁層的能量可能快速耗散，也可能暫時(shí)儲(chǔ)存在磁層內(nèi)部緩慢地釋放，所以現(xiàn)有各種模型很難給出精準(zhǔn)預(yù)報(bào)。本文的推薦模型可以推薦出3 個(gè)最相似事例，預(yù)報(bào)員在實(shí)際應(yīng)用中，可以直接采用最優(yōu)推薦結(jié)果，也可以通過比對(duì)分析，優(yōu)選其中之一作為預(yù)報(bào)參考。推薦結(jié)果不光可以為預(yù)報(bào)員提供直接的Kp指數(shù)預(yù)報(bào)值，還可以提供推薦事例與輸入事例太陽風(fēng)參數(shù)的時(shí)序變化對(duì)比，供預(yù)報(bào)員結(jié)合自身經(jīng)驗(yàn)進(jìn)行預(yù)報(bào)調(diào)整。

圖3 中，紅線表示2001年第90 天03:00 UT 的輸入太陽風(fēng)參數(shù)變化，綠線表示推薦事例之一的2003年第324 天11:00 UT 的太陽風(fēng)參數(shù)變化，其中兩張子圖縱軸分別表示太陽風(fēng)風(fēng)速v和行星際磁場Bz分量，橫軸均為太陽風(fēng)事例的3 h 的時(shí)間序列。由圖可知，輸入的太陽風(fēng)事例風(fēng)速v更大，根據(jù)人工預(yù)報(bào)經(jīng)驗(yàn)，當(dāng)前輸入的太陽風(fēng)可能在地球磁層發(fā)生更多的磁重聯(lián)，單位時(shí)間內(nèi)的磁通量可能更大，因此輸入到地球的能量可能更多，則對(duì)推薦事例的結(jié)果可以向上微調(diào)。而實(shí)際結(jié)果正印證了本文的推測，即實(shí)際輸入的太陽風(fēng)事例對(duì)應(yīng)相鄰的地磁Kp指數(shù)為9–，推薦的太陽風(fēng)事例對(duì)應(yīng)相鄰的地磁Kp指數(shù)為8–。

圖3 2001年第90 天03:00 UT 太陽風(fēng)與2003年第324 天11:00 UT 太陽風(fēng)參數(shù)變化的對(duì)比Fig.3 Comparison of solar wind parameters at 03:00 UT on the 90th day in 2001 and 11:00 UT on the 324th day in 2003

5 結(jié)論

本文以太陽風(fēng)風(fēng)速v和行星際磁場Bz作為輸入特征參數(shù)，利用機(jī)器學(xué)習(xí)相似度算法計(jì)算輸入事例的參數(shù)矩陣與歷史事例的參數(shù)矩陣的距離，并根據(jù)距離排序結(jié)果進(jìn)行推薦。與傳統(tǒng)基于人工神經(jīng)網(wǎng)絡(luò)模型相比，本文的推薦模型基于XGBoost 算法剔除了其它特征參數(shù)對(duì)地磁Kp指數(shù)影響的非線性關(guān)系，直接通過相似的太陽風(fēng)歷史事例的地磁擾動(dòng)影響指導(dǎo)當(dāng)前太陽風(fēng)的預(yù)報(bào)，測試結(jié)果表明該模型具有較好的地磁Kp指數(shù)預(yù)報(bào)效果。

本文的推薦模型避免了人工神經(jīng)網(wǎng)絡(luò)模型預(yù)報(bào)的“黑盒”問題。預(yù)報(bào)員可以通過該模型獲取Kp預(yù)測值，也可以從多個(gè)推薦結(jié)果中看到輸入太陽風(fēng)事例與推薦的歷史事例的參數(shù)變化對(duì)比,使預(yù)報(bào)員可以在預(yù)報(bào)中更多地結(jié)合自身預(yù)報(bào)經(jīng)驗(yàn)，優(yōu)選調(diào)整預(yù)報(bào)結(jié)果。由于太陽風(fēng)與磁層相互作用的復(fù)雜性，推薦模型只是通過歷史相似事件給出當(dāng)前輸入太陽風(fēng)的可能影響，Kp指數(shù)預(yù)測精度的進(jìn)一步提高也有賴于歷史數(shù)據(jù)的積累。綜上，本文的推薦模型可以為預(yù)報(bào)員提供一種更為實(shí)用的短時(shí)預(yù)報(bào)地磁暴的工具，也有利于預(yù)報(bào)員預(yù)報(bào)經(jīng)驗(yàn)的提升。