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

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

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

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

0 引言

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

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

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

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ù)時頻為60 s?？紤]到衛(wèi)星太陽風(fēng)數(shù)據(jù)可能存在缺失，對于缺失數(shù)據(jù)首先進行數(shù)據(jù)清洗：將個別的缺省數(shù)據(jù)近似取值為其前60 s 對應(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分量的值以及行星際磁場強度B。為了便于太陽風(fēng)推薦模型的訓(xùn)練和學(xué)習(xí)，對太陽風(fēng)參數(shù)數(shù)據(jù)均進行了歸一化處理，消除不同量綱對計算的影響。

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

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

2 模型方法

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

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

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

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

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

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

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

為了有效地計算篩選后的特征的權(quán)重，避免多余特征的影響，本文采用單層全連接神經(jīng)網(wǎng)絡(luò)對篩選后保留的太陽風(fēng)風(fēng)速v，行星際磁場Bz分量進行訓(xùn)練，并評估它們各自對地磁擾動影響的權(quán)重。以歸一化后的太陽風(fēng)風(fēng)速v，行星際磁場Bz分量作為單層神經(jīng)網(wǎng)絡(luò)的輸入層參數(shù)，時間分辨率是1 h，輸出層為每小時相應(yīng)的Kp指數(shù)，損失函數(shù)為預(yù)測值和實際值的均方差函數(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)建

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

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

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

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

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

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

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

3 模型測試結(jié)果

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

選取測試用例的最優(yōu)推薦事例（非自身的top1）對應(yīng)Kp指數(shù)與輸入事例對應(yīng)的實際Kp指數(shù)比較，測試結(jié)果如圖2。圖中橫軸是輸入太陽風(fēng)事例對應(yīng)的實際地磁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平均絕對誤差為0.65，Kp指數(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ù)報模型與本文模型的比較。對比結(jié)果表明，本文的推薦模型優(yōu)化了參數(shù)輸入量的選取，Kp預(yù)報結(jié)果的相關(guān)性更好，預(yù)報結(jié)果誤差較小。

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

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

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

4 討論

由于太陽風(fēng)對地球磁層的能量輸入方式極其復(fù)雜且高度動態(tài)，進入地球磁層的能量可能快速耗散，也可能暫時儲存在磁層內(nèi)部緩慢地釋放，所以現(xiàn)有各種模型很難給出精準(zhǔn)預(yù)報。本文的推薦模型可以推薦出3 個最相似事例，預(yù)報員在實際應(yīng)用中，可以直接采用最優(yōu)推薦結(jié)果，也可以通過比對分析，優(yōu)選其中之一作為預(yù)報參考。推薦結(jié)果不光可以為預(yù)報員提供直接的Kp指數(shù)預(yù)報值，還可以提供推薦事例與輸入事例太陽風(fēng)參數(shù)的時序變化對比，供預(yù)報員結(jié)合自身經(jīng)驗進行預(yù)報調(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 的時間序列。由圖可知，輸入的太陽風(fēng)事例風(fēng)速v更大，根據(jù)人工預(yù)報經(jīng)驗，當(dāng)前輸入的太陽風(fēng)可能在地球磁層發(fā)生更多的磁重聯(lián)，單位時間內(nèi)的磁通量可能更大，因此輸入到地球的能量可能更多，則對推薦事例的結(jié)果可以向上微調(diào)。而實際結(jié)果正印證了本文的推測，即實際輸入的太陽風(fēng)事例對應(yīng)相鄰的地磁Kp指數(shù)為9–，推薦的太陽風(fēng)事例對應(yīng)相鄰的地磁Kp指數(shù)為8–。

圖3 2001年第90 天03:00 UT 太陽風(fēng)與2003年第324 天11:00 UT 太陽風(fēng)參數(shù)變化的對比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ù)，利用機器學(xué)習(xí)相似度算法計算輸入事例的參數(shù)矩陣與歷史事例的參數(shù)矩陣的距離，并根據(jù)距離排序結(jié)果進行推薦。與傳統(tǒng)基于人工神經(jīng)網(wǎng)絡(luò)模型相比，本文的推薦模型基于XGBoost 算法剔除了其它特征參數(shù)對地磁Kp指數(shù)影響的非線性關(guān)系，直接通過相似的太陽風(fēng)歷史事例的地磁擾動影響指導(dǎo)當(dāng)前太陽風(fēng)的預(yù)報，測試結(jié)果表明該模型具有較好的地磁Kp指數(shù)預(yù)報效果。

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