衛(wèi)星遙測數(shù)據(jù)相關性知識發(fā)現(xiàn)方法*

楊甲森，孟 新，王春梅

(1. 中國科學院 國家空間科學中心 復雜航天系統(tǒng)電子信息技術重點實驗室， 北京 100190； 2. 中國科學院大學， 北京 100049)

航天器遙測數(shù)據(jù)是地面運管系統(tǒng)判斷其在軌運行狀態(tài)的唯一依據(jù)[1]。作為典型復雜系統(tǒng)，航天器所包含的供配電、姿控、軌控、熱控、有效載荷等分系統(tǒng)[2]之間、分系統(tǒng)內部各模塊之間，均存在著大量電氣、數(shù)據(jù)、熱控接口以及復雜的系統(tǒng)交互，這就決定了反映星載設備狀態(tài)的遙測數(shù)據(jù)之間普遍存在著相關關系。掌握這些相關關系對于實現(xiàn)多元遙測數(shù)據(jù)的關聯(lián)檢測、航天器異常事件原因的深層次挖掘以及特征選擇和數(shù)據(jù)降維等都具有十分重要的意義。

圍繞變量相關性，國內外學者提出了多種有效的相關性度量方法。其中Pearson相關系數(shù)[3]、最大信息壓縮指數(shù)[4]、最小平方回歸誤差等被廣泛用于度量變量間的線性相關，但難以刻畫遙測變量間普遍存在的非線性相關；互信息[5]、信息增益比等測度雖然能夠同時對線性相關和非線性相關進行度量，但這些測度依賴的概率密度函數(shù)的估計較為困難[5-6]，且不具有普適性和等價性[7]。2011年，Reshef等[7]通過基于網(wǎng)格劃分的互信息估計思想提出的最大信息系數(shù)(Maximal Information Coefficient, MIC)，不但能夠同時刻畫變量間線性相關和非線性相關，而且對函數(shù)、超函數(shù)以及分段函數(shù)關系的度量都十分有效，具有較好的普適性和等價性，被廣泛應用于復雜裝備監(jiān)測數(shù)據(jù)[8]、航天器遙測數(shù)據(jù)[9]的相關性分析。然而，在將MIC方法應用于高維度、大規(guī)模衛(wèi)星遙測數(shù)據(jù)的相關性分析時，存在以下問題：①算法時間復雜度高，即使基于動態(tài)規(guī)劃的近似算法，時間復雜度也達到O(n2.4)[10]；②MIC測度偏向于多值變量，取值較少變量的得分往往偏小。提高MIC算法的精度和性能歷來是相關性研究的兩大熱點。文獻[11]基于二次優(yōu)化過程提高了MIC精度但卻再度提高了復雜度；文獻[12]采用圖形處理單元和現(xiàn)場可編程門陣列組成的異構加速器提高了MIC計算的性能；文獻[13]提供了一種基于并行處理的MIC快速計算跨平臺工具；文獻[14]基于MapReduce模型對MIC算法進行了并行化設計；文獻[15]采用C語言重新實現(xiàn)了MINE套件。上述文獻針對MIC算法性能提升的研究，多是圍繞硬件環(huán)境、并行運算等算法實現(xiàn)的方式展開，并未改變其等深劃分后基于動態(tài)規(guī)劃算法尋優(yōu)的本質。此外，已有文獻尚無關于MIC測度偏向多值變量問題的討論。

1 最大信息系數(shù)

1.1 MIC算法原理

MIC是以互信息為基礎的，其主要思想是基于一種認識：如果兩個變量存在相關關系，那么可以在兩個變量構成的散點圖上繪制網(wǎng)格，數(shù)據(jù)在網(wǎng)格中的分布情況可以反映二者之間的相關關系。網(wǎng)格繪制中，MIC方法考慮兩個因素：①網(wǎng)格劃分的數(shù)量；②網(wǎng)格劃分的方法。其原理可通過以下定義[7]來說明。

定義1(特定劃分下最大信息系數(shù)I(D,s,t))對二維有序對數(shù)據(jù)集D(X,Y)，分別在X、Y方向上進行s、t段劃分。定義該劃分對應的網(wǎng)格為G，G劃分下D的概率分布為D|G、互信息為I(D|G)，那么此劃分下的最大信息系數(shù)為：

I(D,s,t)=maxI(D|G)

(1)

式中max的含義是指：將二維變量劃分為s、t段的方法有多種，如等寬、等深劃分[16]，每一種方法對應不同的互信息值，I(D,s,t)為其中的最大值。I(D|G)計算方法如下：

(2)

式中，p(s,t)為聯(lián)合概率密度，p(s)、p(t)為邊緣概率密度。根據(jù)大數(shù)定理，當觀測數(shù)據(jù)量足夠大時，p(s,t)的計算可用落入網(wǎng)格中的樣本數(shù)量占比來近似，p(s)、p(t)分別用落入(k,k+1)和(l,l+1)網(wǎng)格的樣本數(shù)量占比來近似，其中k∈[0,s-1],l∈[0,t-1]。

定義2(特征矩陣M(D))特征矩陣M(D)第s行、t列元素的定義為：在s,t段網(wǎng)格劃分下的最大信息系數(shù)(定義1)的歸一化矯正，如式(3)。s,t為任意正整數(shù)，即M(D)是無限維矩陣。

(3)

定義3(最大信息系數(shù)MIC(D))設D的樣本容量為n，那么：

MIC(D)=maxst≤B(n){M(D)s,t}

(4)

式中，B(n)為網(wǎng)格劃分數(shù)量的上限，它將特征矩陣M(D)限定為有限維，文獻[7]給出了B(n)的推薦值為n0.6。

1.2 基于動態(tài)規(guī)劃的MIC近似算法

在計算I(D,s,t)時，為了避免網(wǎng)格窮舉切割，進行遍歷尋優(yōu)，降低MIC方法計算的復雜度，文獻[7]提出了一種基于動態(tài)規(guī)劃算法來近似求解MIC的方法。詳細步驟如下：

步驟1：在X方向，以等值樣本落入同一網(wǎng)格為原則，確定X方向s=2段的等深劃分Q，如圖1所示。等深的含義是指，劃分后每一個分段的樣本點數(shù)相差不多。

圖1 X方向等深劃分Fig.1 X direction equal partition

圖2 Y方向候選劃分Fig.2 Y direction candidate partition

步驟2：在Y方向，以等值樣本落入同一網(wǎng)格為原則，確定Y方向候選劃分P，如圖2所示。其方法是以X方向分割線與數(shù)據(jù)曲線(按Y值升序順序，逐點連接而成的線)的交點作Y軸的垂線，Y方向等值樣本若不在同一X劃分Q中，則需增加一個Y方向的候選劃分。

F(cr,2)

=max1≤i

(5)

F(cr,l)

(6)

F(ck,l)是所有樣本l分段下，H(P)-H(Q,P)的最大值。由于互信息I(Q;P)=H(Q)+H(P)-H(Q,P)，在X方向Q劃分不變的情況下，H(Q)為定值，因此F(ck,l)+H(Q)即為所有樣本l分段下的最大互信息值。

步驟7：將縱坐標與橫坐標交換，重復步驟1～6，求得Y方向等深劃分下的特征矩陣MY(D)。

步驟8：以兩個特征矩陣MX(D)、MY(D)中的最大元素作為兩個變量的MIC值。

近似算法計算量主要集中在步驟4，即動態(tài)規(guī)劃算法求F值部分，該部分時間復雜度為O(k2st)[7]，等深劃分循環(huán)次數(shù)為B/2，因此其整體復雜度為O(k2stB)=O(k2B2)=O(n2.4)。

2 Mini Batch K-Means算法

Mini BatchK-Means算法是K-Means算法的變種。與傳統(tǒng)K-Means方法相比，Mini BatchK-Means在每次迭代更新質心的過程中，采用隨機抽樣獲得的數(shù)據(jù)子集進行更新。實踐證明Mini BatchK-Means在數(shù)據(jù)量較大時，可以有效地減少算法收斂時間，但其準確度稍有下降[17]。

給定樣本集Z={z1,…,zn}，n為樣本容量。Z將被劃分為k個簇，簇的中心為C={c1,…,ck}。Mini BatchK-Means聚類的迭代步驟如下。

步驟1：從Z中隨機選擇k個樣本作為初始中心C。

步驟2：從Z中隨機抽取容量為b的樣本子集L={l1,…,lj,…,lb}，組成一個Batch。

步驟3：對L中每一個樣本點lj，計算其與k個簇類中心的相似度，將樣本點lj劃入相似度最大的簇。

步驟4：L中所有樣本經(jīng)過步驟3后，根據(jù)各樣本的簇標號重新計算聚類中心。

步驟5：判斷是否滿足聚類結束的條件(如達到迭代次數(shù)t)，若未滿足，回到步驟2，否則進入步驟6。

步驟6：對Z中的每個樣本點，根據(jù)其與k個聚類中心的相似程度，將其劃分給相似度最大的簇。

Mini BatchK-Means算法的主要計算量集中在步驟3，即計算每個樣本點的相似度并確定歸屬，其時間復雜度為O(knt)=O(n)。

3 改進算法

3.1 改進的歸一化因子

式(3)中，MIC采用最大熵log2(min(s,t))對互信息進行歸一化矯正。當變量取值較少，其數(shù)據(jù)分布只集中在少量網(wǎng)格中時，式(3)計算的MIC測度往往偏小，而這種取值較少的變量，在航天器遙測領域普遍存在。以量子衛(wèi)星遙測參數(shù)“延時遙控指令數(shù)”“注入軌道剩余點數(shù)”為例，其取值如圖3所示，二者采用式(3)方法計算的MIC值為0.41，而皮爾遜相關系數(shù)達到了0.95。

分析MIC取值偏小的原因為：歸一化因子采用的最大熵，對應的是特殊的均勻分布，而不是變量的實際分布。本文以信息熵代替原有最大熵，作為歸一化因子對互信息進行矯正，如式(7)所示，式中變量分布對分子、分母具有等同的貢獻而被相互約減，可有效降低變量分布對MIC測度的影響。其中Q為X方向的s段劃分，P為Y方向的t段劃分。

(7)

(a) 遙測時間序列(a) Time series of telemetry

(b) 二維遙測分布(b) Two-dimensional telemetry distribution圖3 數(shù)據(jù)分布集中在少量網(wǎng)格Fig.3 Data is concentrated in a small number of grids

3.2 改進的MIC算法

3.2.1 相關定義

傳統(tǒng)Mini BatchK-Means算法的初始聚類中心是隨機選擇的，聚類結果具有不可預見性[18]。本文定義4、5，提出一種依據(jù)樣本分布選擇初始聚類中心的方法。

定義5(初始k個聚類中心)S′={〈Value′i,Count′i〉|i=1,…,m}是一維序列統(tǒng)計信息S按Count的降序排列，則初始k個聚類中心C={c1,…,ck}={Value′1,…,Value′k}。

3.2.2 小批量K均值MIC改進算法(MBKM_MIC)

MBKM_MIC由兩個階段組成，如算法1所示，第一階段(行1～2)獲取一維序列的統(tǒng)計信息，并按Count字段降序排列，最頻繁的取值將被作為第二階段初始的聚類中心；第二階段(行3～8)進行MIC值計算，主要包括3個步驟：①利用Mini BatchK-Means算法將X軸劃分為s段、Y軸劃分為t段；②計算在s、t劃分下的互信息值并進行歸一化得到信息系數(shù)；③循環(huán)①、②遍歷st≤B條件下所有劃分的信息系數(shù)，選擇極大值作為MIC值。

算法1 MBKM_MIC算法

算法1調用了算法2～4，描述如下。

算法2 GetStatisticInfo算法

3.2.3 時間復雜度分析

第一階段完成4個快速排序，時間復雜度為O(4nlog2n)；第二階段中，X、Y方向的劃分為Mini BatchK-Means算法，其復雜度為O(n)，互信息計算需要遍歷所有數(shù)據(jù)，其復雜度亦為O(n)，故第二階段復雜度為O(3nB)=O(n1.6)。算法總復雜度為O(4nlog2n)+O(n1.6)=O(n1.6)。

算法3 MBKM算法

算法4 MI算法

4 實驗結果與分析

為驗證所提方法的有效性，進行了三組實驗。第一組實驗采用最大熵、信息熵兩種歸一化因子對遙測數(shù)據(jù)進行相關性分析，用于驗證信息熵因子對取值較少變量的適用性；第二組實驗用于驗證MBKM_MIC方法的效能，分析該方法得分與動態(tài)規(guī)劃方法得分的差異，評估所提方法的工程可用性；第三組實驗，用于驗證MBKM_MIC方法的處理效率，確認該方法對大規(guī)模數(shù)據(jù)相關性分析的適用性。

4.1 改進的歸一化因子實驗驗證

采用量子衛(wèi)星遙測數(shù)據(jù)，對改進前后歸一化因子的應用效果進行試驗。其中最大熵因子算法采用基于動態(tài)規(guī)劃的ApproxMIC[7]算法，信息熵因子算法是在該算法C語言代碼基礎上進行的改進。

實驗選用3組遙測數(shù)據(jù)，如表1所示(變量維數(shù)為所選時段方差不為0的遙測變量個數(shù))。對每組數(shù)據(jù)中兩兩變量的MIC值進行計算，以閾值0.8作為相關性篩選條件。實驗中發(fā)現(xiàn)的相關關系數(shù)量對比如圖4所示。

表1 三組量子衛(wèi)星數(shù)據(jù)

圖4 兩種歸一化因子發(fā)現(xiàn)的相關關系數(shù)量Fig.4 Number of correlation relationships found by two normalization factors

圖4中，采用最大熵和信息熵歸一化因子發(fā)現(xiàn)的相關關系總計分別為259、305組。經(jīng)核驗：①采用最大熵因子發(fā)現(xiàn)的相關關系是信息熵發(fā)現(xiàn)關系的子集；②同組變量的相關關系，信息熵因子計算的MIC值均高于最大熵因子；③多出的46組相關關系均為變量取值偏少導致最大熵因子下得分較小、信息熵因子下得分較大的案例，以第3.1節(jié)中“延時遙控指令數(shù)”“注入軌道剩余點數(shù)”為例，其信息熵歸一化MIC取值達到了0.94。因此，較之最大熵，采用信息熵歸一化因子能夠發(fā)現(xiàn)遙測數(shù)據(jù)中更為廣泛的相關關系，對于取值較少變量的相關性分析，信息熵歸一化因子仍具有良好的適用性。

4.2 MBKM_MIC算法效能實驗驗證

分別采用MBKM_MIC、ApproxMIC方法，對表1中量子衛(wèi)星遙測數(shù)據(jù)進行了處理。MBKM_MIC方法的實驗參數(shù)設置為B=n0.6，time=3，bsize=100；ApproxMIC參數(shù)設置為B=n0.6，c=15。

三組數(shù)據(jù)的處理結果如圖5～7所示，圖中橫軸為三組量子數(shù)據(jù)中兩兩組合的變量組索引，圖(a)為ApproxMIC計算結果，圖(b)為MBKM_MIC方法計算結果，圖(c)為兩個測度取值的偏差。三組數(shù)據(jù)兩個測度取值的平均絕對誤差(假設ApproxMIC方法取值為標準值)分別為0.103、0.035、0.034，兩種方法取值趨勢基本一致，誤差也非常小。

從圖5～7可見，部分變量組的MBKM_MIC方法測度取值低于ApproxMIC方法，這是由MBKM_MIC直接網(wǎng)格劃分獲得測度取值，而ApproxMIC經(jīng)歷了動態(tài)規(guī)劃尋優(yōu)過程導致的；另有部分變量組的MBKM_MIC取值高于ApproxMIC方法，分析是由二者分別采用信息熵和最大熵作為歸一化因子引起的。

(a) ApproxMIC

(b) MBKM_MIC

(c) ApproxMIC-MBKM_MIC圖5 軌道實時包處理結果比較Fig.5 Comparison of real-time orbit package processing results

(a) ApproxMIC

(b) MBKM_MIC

(c) ApproxMIC-MBKM_MIC圖6 姿控慢速包處理結果比較Fig.6 Comparison of slow attitude package processing results

(a) ApproxMIC

(b) MBKM_MIC

(c) ApproxMIC-MBKM_MIC圖7 平臺實時包處理結果比較Fig.7 Comparison of real-time platform package processing results

4.3 MBKM_MIC算法效率實驗驗證

采用與第4.2節(jié)實驗相同的參數(shù)設置，對三組量子數(shù)據(jù)的處理時間進行統(tǒng)計，實驗共進行了10次，分別對應樣本容量n為1024～10 240 ，間隔為1024，性能測試結果如圖8所示。

(a) 軌道實時包(a) Real-time orbit package

(b) 姿控慢速包(b) Slow attitude package

(c) 平臺實時包(c) Real-time platform package圖8 處理性能比較Fig.8 Comparison of processing performance

由圖8可知，對于遙測數(shù)據(jù)的相關性分析，MBMK_MIC方法的性能普遍優(yōu)于ApproxMIC，且樣本規(guī)模越大，性能優(yōu)勢越明顯。

5 結論

遙測數(shù)據(jù)相關性發(fā)現(xiàn)在衛(wèi)星數(shù)據(jù)分析、故障診斷中具有關鍵的“導航”作用。從高維的航天器遙測數(shù)據(jù)中挖掘出具有可解釋性的相關性知識，可以快捷、高效地發(fā)現(xiàn)星載設備間內在的關聯(lián)關系。本文針對MIC方法應用于大規(guī)模遙測數(shù)據(jù)相關性分析過程中，處理性能較低、偏向于多值變量的問題，提出了改進歸一化因子、以Mini BatchK-Means聚類為前驅過程的改進MIC方法。試驗結果表明，改進的歸一化因子方法對取值較少的變量也具有良好的適用性；改進的MIC方法，在計算結果與ApproxMIC測度偏差可承受的前提下，顯著提升了數(shù)據(jù)處理的性能，是一種適用于大規(guī)模遙測數(shù)據(jù)相關性分析的有效方法。