基于機器學習相似度算法的Kp 指數(shù)預報*

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

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

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

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

0 引言

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

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

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

1 數(shù)據(jù)

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

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

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

2 模型方法

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

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

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

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

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

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

2.2 參數(shù)權重計算

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

2.3 推薦模型構建

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

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

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

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

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

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

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

3 模型測試結果

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

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

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

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

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

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

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

4 討論

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

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

圖3 2001年第90 天03:00 UT 太陽風與2003年第324 天11:00 UT 太陽風參數(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 結論

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

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