基于數據挖掘方法的空間大氣模型修正

廖川， 白雪， 徐明

(北京航空航天大學宇航學院, 北京 100083)

數據挖掘(又稱從數據中發現知識)起源于20世紀80年代后期，在20世紀90年代有了突飛猛進的發展，數據挖掘提供了發現隱藏在大型數據集中的模式的技術，關注可行性、有用性、有效性和可伸縮性問題[1]。大數據指的是以不同形式存在于數據庫、網絡等媒介上蘊含豐富信息的規模巨大的數據。大數據是一個寬泛的概念，其基本特征可以用4個V(Volume、Variety、Value和Velocity)來總結，隨著大數據的發展，也有演變成5個甚至6個V(Volume、Variety、Value、Velocity、Veracity和Variability)的趨勢[2]。

近年來，隨著各行業從業人員和產品的開發應用，數據挖掘技術得到了長足的發展。例如：金融證券行業利用大數據預測股票價格的波動和走勢，分析用戶投資習慣；電力行業利用大數據進行電力負載情況的統計，加強分時分流管理；互聯網行業利用大數據挖掘上網用戶的使用習慣，為用戶推薦相應優質內容。而在航天領域，國內外也有相當多的學者進行了深入研究。

Chen等[3]站在數據庫研究人員的角度，對數據挖掘技術做出了綜合性闡述，調查各類數據挖掘新技術，包括挖掘關聯規則的Apriori算法、數據立方體方法、基于決策樹(Decision Tree，DT)的數據分類方法等內容。Tanner等[4]提出了星上數據挖掘的概念，將傳感器數據進行星上處理，可以提升網絡通信能力，并降低成本，可適用于無人自主航天器、火星表面測繪衛星、生物識別系統等。Sánchez-Sánchez等[5]通過大數據和神經網絡(Neural Network,NN)學習連續確定非線性系統的狀態反饋最優控制問題，利用深度神經網絡對最優控制問題生成的軌道進行有監督學習。最終利用得到的神經網絡求解Hamilton-Jacobi-Belman方程，借此直接在星上實時生成近優控制行為。Hennes等[6]考慮航天器在主帶小行星間的小推力最優轉移問題，針對相位值、最大初始質量和最大最終質量，采用計算智能技術(包括分解框架的多目標進化算法、超體積指標和機器學習回歸)，分別對3個目標函數進行最優設計。

李德仁等[7]分析了空間數據挖掘和知識發現(Spatial Data Mining and Knowledge Discovery，SDMKD)的內涵和外延，利用SDMKD從空間數據庫中自動或半自動地挖掘事先未知卻潛在有用的空間模式。宮輝力等[8]針對國內缺乏對海量衛星遙感數據的有效數據挖掘和知識發現技術，設計了一種面向地學應用的衛星多源遙感圖像數據挖掘系統框架和原型。Peng和Bai[9]針對空間目標編錄中通常缺失的面質比信息，利用隨機森林(Random Forest，RF)方法學習一致性誤差和面質比間的關系，從而確定空間目標的面質比類型。

異常檢測與模式挖掘作為數據挖掘的一個重要分支，也逐漸成為研究熱點方向。胡小平等[10]分析了液體火箭發動機的工作和故障特點，提出利用數據挖掘方法從數據倉庫的角度對液體火箭發動機進行故障檢測和診斷的策略。徐宇航和皮德常[11]為了及時發現衛星故障情況并提取異常模式，提出一種基于PrefixSpan算法，對衛星異常數據集進行模式挖掘。肇剛和李言俊[12]研究時間序列數據挖掘技術現狀，提出一種基于時間序列數據挖掘的航天器故障診斷方法，包括7個步驟：確定診斷對象、遙測數據選取、數據再處理、時間序列特征表示、數據挖掘方法執行、挖掘結果的解釋與評估、軟件系統開發。對促進航天器故障診斷技術的發展和完善有重要意義。

某型號衛星是一顆典型的低軌衛星，在其軌道預報和實際運行中，大氣阻力是最主要的非保守力。本文基于以上現狀，將利用基于決策樹的隨機森林方法、神經網絡、支持向量機(Support Vector Machine，SVM)這3類數據挖掘的方法，學習大氣阻力在不同誤差情況下導致的軌道預報的偏差，借此對所使用的大氣模型進行修正，并評估3種方法的修正性能，以及給出方法間的對比。本文的創新點如下：1)針對有監督學習方法需要先驗類標簽的問題，利用2種不同精度的軌道預報模型作為仿真環境，產生訓練數據和測試數據。2)采用了大數據處理技巧中的3種數據挖掘方法與傳統數理統計方法相比，具有快速處理大數據集(本文中數據點達到18萬左右，甚至更多)、能夠挖掘隱藏在軌道預報微小誤差中的潛在信息等優勢，驗證了數據挖掘方法的應用可行性，為未來在實際運行環境下進一步修正大氣模型提供參考。

1 軌道動力學模型

高精度軌道預報需要對空間環境進行精確建模，軌道預報的主要攝動因素包括地球引力、太陽引力、太陽光壓和大氣阻力。首先，對各攝動因素進行了精確的數學建模；然后，提出2種軌道動力學模型，作為軌道預報和產生訓練數據的基礎。

1.1 高階地球引力場

根據參考文獻[13]給出的相關內容，地球引力勢可以寫成

(1)

式中：

(2)

GMe為地球引力常數，G為萬有引力常數，Me為地球質量；r為地心固連坐標系中的航天器位置矢量；φ和λ分別為航天器的地心緯度和地心經度；Re為地球的赤道半徑；Pnm為n階m次締合勒讓德多項式；Cnm和Snm分別為相應的重力勢系數。

這種形式的引力場模型可以擴展至任意階次，而不僅僅限制在帶諧項，因此可以用以精確描述地球引力場。

(3)

式中：

(4)

式(1)～式(4)得到的是地心固連坐標系中的加速度，為統一各攝動加速度，后文所有其他加速度也將轉換到地心固連坐標系下，如式(5)所示：

(5)

式中：U(t)為坐標轉換矩陣，用以刻畫地球自轉，實際計算中還需考慮歲差和章動效應；r為航天器位置矢量，下標e和s分別表示地心固連坐標系和空間固定坐標系。

1.2 非地球引力攝動

太陽引力引起的攝動加速度可在地心慣性坐標系下表示為

(6)

式中：ri和s分別為航天器和太陽的地心坐標；MS為太陽的質量。

由于太陽輻射光壓導致的攝動加速度可以表示為

(7)

式中：PS為太陽輻射壓；1 AU為一個天文單位，1 AU=1.496×108km；rS為太陽地心矢量；n為航天器表面的單位法線矢量，A為航天器的迎風面積；eS為太陽方向單位矢量；θ為矢量n和矢量eS間的夾角；ε為航天器所用材料的反射率。衛星在運行過程中，由于存在太陽光遮擋的問題，本文在進行軌道預報的過程中，將采用參考文獻[14]中的地影預報算法，計算衛星當前是否處于地影區，從而判斷當前狀態是否存在太陽光壓攝動。

對于低軌衛星而言，大氣阻力是最主要的非保守力。大氣阻力引起的攝動加速度可表示為

(8)

式中：CD為阻力系數；m為航天器質量；ρ為航天器所處位置的大氣密度，本文采用MSISE86大氣密度模型[15]；vr為航天器相對于大氣的速度；ev為相對速度的單位矢量，即ev=vr/vr。

1.3 精確模型和誤差簡化模型

1.1節和1.2節已經給出了航天器軌道動力學模型，該模型可以用于航天器軌道預報。為產生分類器的訓練數據和測試數據，本節將給出2種模型，精確模型將作為航天器軌道預報的基準，模擬航天器在軌“真實”情況；不同的誤差模型將人為對大氣模型施加誤差，從工程實踐角度出發，CD(A/m)ρ可以被視作一個阻力系數CS，因此本文實際上是對CS進行修正，因此式(8)在誤差簡化模型中可寫成如下形式：

(9)

通過精確模型產生的軌道，與誤差簡化模型產生的軌道的各項數據差值，即可作為分類器的訓練數據和測試數據。2種模型具體的參數設置由表1所示。

表1 精確模型與誤差簡化模型參數Table 1 Parameters of accurate model and error simplified model

2 數據挖掘方法

數據挖掘能夠從大量雜亂的數據中挖掘出有價值的信息，但是其過程是相當復雜的。針對本文的研究內容及數據的特性，本節將介紹3種分類方法，利用這3種分類方法就能夠學習第1節得到的訓練數據，并通過測試數據對3種分類方法進行準確度檢驗。

2.1 決策樹與隨機森林

決策樹又稱分類樹，是應用最為廣泛的分類方法之一。如圖1所示，決策樹是一種監督學習方法，通過不斷地劃分數據，使依賴變量的差別最大，最終目的是將數據分類到不同的組織或不同的分枝，在依賴變量的基礎上建立最強的歸類，其訓練結果是類似流程圖的樹結構。

但是決策樹是一種“貪心算法”，在應用過程中，每一步的判斷，僅是針對當前測試做出最優選擇，并不考慮全局結果。如果將多棵樹以某種關系結合起來，對數據進行分類，以解決單一決策樹泛化能力弱的缺點，這就是隨機森林。

圖1 決策樹Fig.1 Decision tree

隨機森林實際上是套袋法與決策樹的結合，如圖2所示，隨機森林的隨機性體現在：

1) 從樣本中隨機選擇n個樣本。

2) 從所有屬性中隨機選擇k個屬性。

3) 選擇最佳分割屬性作為節點建立決策樹；

4) 以上步驟重復j次，即得到j棵決策樹，進而完成隨機森林訓練。

5) 分類問題中，通過投票決定數據的類別。回歸問題中，由j棵決策樹預測結果的均值作為最后的預測結果。

隨機森林的優點很多，對于多種數據，隨機森林可以產生高精度的分類器；建造森林時，隨機森林可以在內部對于一般化后的誤差產生不偏差的估計；在決定類別時，可以評估變數的重要性等等。當然隨機森林的缺點也是顯而易見的，因為要訓練m棵決策樹，其訓練過程會是訓練單棵決策樹的數倍。

圖2 隨機森林Fig.2 Random forest

2.2 神經網絡

神經網絡是應用類似于人腦神經突觸連接的結構進行信息處理的數學模型。在早期對人腦的理解中，人腦主要由稱為神經元的神經細胞組成，神經元通過軸突的纖維絲連接在一起，當一個神經元受到刺激時，神經脈沖通過軸突從一個神經元傳遞到另一個神經元。一個神經元通過樹突連接到其他神經元的軸突，樹突是神經元的延伸物。樹突與軸突之間的連接點叫做神經鍵。人腦通過在同一脈沖反復刺激下改變神經元之間的神經鍵連接強度來進行學習。類似于人腦的結構，神經網絡也是由一組相互連接的節點和有向鏈構成。

圖3 神經網絡感知器Fig.3 Neural network perceptron

圖3為神經網絡最簡單的模型——感知器，包含多個輸入節點、一個輸出節點，節點即是神經元。在感知器中，每個輸入節點通過一個加權的鏈連接到輸出節點，該加權鏈用以模擬神經元之間神經鍵的連接強度，如人腦一樣，訓練一個感知器就是不斷調整鏈的權值，直到能擬合出輸入輸出間的關系為止。

本文將采用前饋負向傳播神經網絡，前饋是指：每一層的神經元僅與前一層連接，且輸出給下一層，各層之間沒有反饋，不對網絡的參數進行調整。負向傳播是指：誤差會進行方向傳播，用以調整權值和閾值。至少含有一個隱藏層的多層神經網絡可以用來擬合任何目標函數，且可以處理冗余特征，但是神經網絡對于噪聲十分敏感，訓練過程也很耗時。

2.3 支持向量機

支持向量機是建立在統計學習理論基礎上的機器學習方法。支持向量機通過構建分割多類的超平面，在構建過程中，會試圖將多類之間的分割達到最大化(如圖4所示)。

支持向量機的理論包括3個要點：最大化間距、對偶理論和核函數。“最大化間距”實際上是一個優化問題，該優化問題，可以利用拉格朗日對偶理論變換到對偶變量的優化問題。然而在輸入空間中，數據可能是不可分的，支持向量機就需要通過非線性映射，將數據映射到某個其他點積空間。本文所使用的支持向量機將采用高斯核函數：

(10)

式中：X為空間中任一點;Y為核函數中心;σ為函數的寬度參數。高斯核函數對于噪聲有著較好的抗干擾能力，但是其參數決定了函數作用范圍，超過了這個范圍，數據的作用就“基本消失”。

圖4 支持向量機分類方法示意圖Fig.4 Schematic diagram of classification method of SVM

3 分類器訓練與評估

首先，將利用第1節提出的2種軌道動力學模型進行軌道預報，并生成訓練數據；然后，將訓練數據導入第2節中介紹的各分類方法，進行適當參數調整后完成訓練。

某型號衛星是中國第一代傳輸型立體測繪衛星，于2010年8月24日發射成功，其軌道高度為500 km左右，運行過程中，大氣阻力對其造成較大的軌道衰減。因此本節將以該衛星在軌實際狀態作為軌道預報起點，從而生成一系列訓練數據。

3.1 訓練數據生成

仿真起始時刻為2016年5月31日21：07：05，衛星真實在軌狀態如表2所示，其中：a為半長軸；e為偏心率；i為軌道傾角；Ω為升交點赤經；ω為近地點幅角；M為平近點角;rx、ry、rz為衛星在慣性系下的位置三軸分量；vx、vy、vz為衛星在慣性系下的速度三軸分量。

表2 仿真初始狀態Table 2 Initial state of simulation

向后進行1 d軌道預報可以得到多條衛星軌道，預報間隔為10 s，每一個點都對應一定時間后由25個變量組成的衛星在軌狀態Xi：

(11)

式中：Δt為該點的預報時長；前標Δ表示誤差簡化模型與精確模型間的差值；下標i表示該組狀態的編號。

因此可以得到181 440組衛星的在軌狀態，且每一組狀態都對應一個人為施加在大氣模型上的誤差值εi，即可以得到181440個樣本數據及其類別屬性的組合[Xiεi]作為訓練數據，訓練數據的如表3所示。

由表3可以看出，不同誤差值的模型進行軌道預報，得到的在軌狀態數據相差很小，無法用傳統方法得到這些值與誤差值間的關系，更無法通過這些值準確得到誤差值大小，因此本文采用數據挖掘的技巧，對以上數據進行處理、分類，這也是本文的創新點之一。

3.2 分類器訓練

將3.1節得到的訓練數據導入各分類模型，并適當設置參數，即可完成分類器訓練的初始化。各分類模型的訓練參數設置如表4～表6所示。

3.3 訓練結果與評估

隨機森林方法分類結果如圖5所示，橫坐標為181 440組狀態的編號，縱坐標為該組狀態對應的大氣模型誤差值。因為隨機森林的訓練和預測中采用投票機制決定最終結果，所以圖5(b)的結果是通過選擇圖5(a)的最大概率取值得到的。

表3 部分訓練數據Table 3 Partial training data

表4 隨機森林訓練參數設置Table 4 Training parameter setting of random forest

表5 神經網絡訓練參數設置

表6 支持向量機訓練參數設置Table 6 Training parameter setting of SVM

可以看到，5(a)中顏色越接近紅色，表示概率越大，因為隨機性和投票機制，最終的精確度高達99.99%。

神經網絡訓練結果如圖6所示，圖6(a)為在5 461次循環中，神經網絡通過不斷調整神經鍵上的權值，逐步減小性能函數的取值，于5 451次循環時達到最小值0.014 55，并在第5 461次循環后達到最大驗證數據失敗的次數，停止循環；圖6(b)為在5 461次循環中，性能梯度、動量和驗證數據失敗次數的變化，性能梯度于5461次循環時達到最小值0.022 783，動量于5 461次循環時達到最小值1×10-7，驗證數據失敗次數于5 461次循環時達到最大值10次。

圖7橫坐標為181 440組狀態的編號，縱坐標為該組狀態對應的大氣模型誤差值，可以看到訓練結果相比隨機森林已經有了較大的降幅，精確率僅有0.011 6%。而且通過圖8可以看出，在預報時長較小時，因為不同模型得到的在軌狀態差值很小，這一段內的結果誤差尤其劇烈，甚至超過了100%，隨著預報時長的增加，誤差將會逐步減小。支持向量機分類結果如圖9所示，結果誤差較大，其精確度僅為50.7%。對比圖7，神經網絡分類結果僅在一定范圍內波動，而支持向量機的結果可能在全域內離散波動。

圖5 混淆矩陣和隨機森林分類結果Fig.5 Classification results of confusion matrix and random forest

圖6 神經網絡訓練過程Fig.6 Training process of neural network

圖7 神經網絡分類結果Fig.7 Classification results of neural network

對比3種方法的性能(見表7)，隨機森林訓練時長最短，精確度最高，因此該方法的性能最佳，考慮到森林中包含決策樹的數量為50棵，減少決策樹的數會進一步縮短訓練時長。而神經網絡和支持向量機的應用結果都不理想，原因在于本文中使用的在軌狀態和大氣模型誤差值的訓練數據，對于這2種方法來說學習較為困難，尤其是神經網絡，容易出現過擬合以及不學習，即過早達到最大驗證失敗次數的情況。

圖8 神經網絡分類誤差變化Fig.8 Classification error variation of neural network

圖9 支持向量機分類結果Fig.9 Classification results of SVM

表7 3種分類方法性能對比Table 7 Performance comparison of three classification methods

4 結 論

本文以某型號衛星的軌道預報及數據挖掘方法作為切入點，完成了以下內容的研究分析：

1) 對衛星在軌的空間環境進行數學建模，基于各個攝動力因素，提出了2種軌道動力學模型，精確模型用于模擬衛星在軌真實情況，誤差簡化模型通過人為對大氣模型施加-100%～100%的誤差，作為生成訓練數據的基礎。

2) 對本文利用到的支持向量機、神經網絡、隨機森林3種分類方法進行了概述和總結，闡明了3種方法的原理及特點。

3) 利用2種模型向后進行1 d的軌道預報，預報結果作差，由此產生181 440組訓練數據，并導入3種分類模型訓練。

4) 從訓練結果可以看出，數據挖掘方法能夠對大氣模型的誤差進行較好的反演，尤其是隨機森林的方法，反演的結果高達99.99%。

本文的結果為衛星軌道預報以及大氣模型的建立提供了新思路，結合時下熱門的大數據概念，將數據挖掘方法應用至傳統的航天領域,具有良好的參考意義，為航天領域應用人工智能、大數據等技術奠定了基礎。