彭 藝,馮小虎,賈樹澤,韓 琦
(國家衛星氣象中心,北京 100081)
衛星狀態的確定,在當前條件下,一般還是采用技術人員手動操作的方式,來確認衛星是否處于異常狀態,這些方法過于依賴地面工作人員的經驗,難以實現異常檢測的智能化。另外,一旦遙測數據在閾值范圍內出現異常,閾值法在使用過程中就會出現問題,其結果就是出現故障漏報。
此外,這些方法只考慮單個變量的變化情況,未考慮相互之間的耦合關系,即某個變量的異常可能造成其它變量的異常,甚至造成重大故障的發生。
當前常見的故障檢測方法有三種:一是數據驅動法,二是知識方法,三是解析模型法。受到廣泛研究的數據驅動法包括:一是以信號處理為基礎的方法,二是以多元統計為基礎的方法,三是以機器學習為基礎的方法[1-3]。其中,第二種方法適用于復雜的故障檢測,一般經常采用方法包括:一是獨立成分分析(Independent Component Analysis,ICA);二是偏最小二乘(Partial Least Square,PLS);三是主元分析(Principal Components Analysis,PCA)[4-6]。但是,上述這些方法在建模方法上無一例外是線性化的,就實踐過程中使用的數據來看,絕大部分是非線性的。就如何解決非線性這個難題,很多研究者提出改進算法,其中得到普遍認可的是以核函數為基礎的算法。如Lee[7]提出新的方法KPCA,可以監測非線性連續過程。降維過程中,如果使用ICA、PLS、PCA這些以全局為基礎的方法,會導致數據點之間原有近鄰關系出現扭曲,就會出現低維數據中保留的局部關系無效的情況[8]。近幾年,流行學習方法在機器學習、數據挖掘和模式識別領域顯得日益重要,此方法特點為:以原始數據樣本的局部信息作為基礎;從中找出低維表示,它的要求是能夠維持局部結構,而且是最大限度的。它與提取全局信息的PCA等方法存在本質區別。當前普遍使用的方法有四種:一近鄰保持嵌入(Neighborhood preserving embedding,NPE)[9];二局部保持投影(Locality preserving projections,LPP)[10];三等距映射(Isometric mapping,Isomap)[11];四局部線性嵌入(Locally linear embedding,LLE)[12]。LPP、NPE等為線性流行學習方法,不適用于非線性過程。Isomap、LLE等非線性流行學習算法,在對訓練數據降維過程中,僅可以使用隱晦的非線性計算方法,而無法在降維過程中采取新的數據樣本,其原因是無法獲得輸入空間和特征空間之間的簡單直接的映射關系[13]。基于多元統計分析的傳統PCA、流行學習等方法不能有效地監測數據空間局部結構特征變化,所以利用貢獻圖法進行故障識別的精度就會差一些。
針對以上問題,本文提出PCA-DNMFSC算法構建衛星故障檢測模型,遙測數據分解的基矩陣的稀疏化,能使異常局部特征凸顯,有效監測特征的變化;動態化表示當前時刻的觀測變量以表達觀測值間的時序相關性,能充分考慮前后數據樣本的特征;采用PCA對DNMFSC進行初始化處理,有效提高算法穩定性;將某段時期內貢獻值進行權加和來確定異常發生是由哪些變量產生的以避免單點的統計量貢獻值無法很好地描述系統的異常。
以局部表示為基礎的非負矩陣分解(Non-negative Matrix Factorization,NMF)屬于一種稀疏表示方法,但其存在三個問題:一是控制困難;二是稀疏程度較弱;三是稀疏能力較差[14]。然而,稀疏化非負矩陣分解(Non-negative Matrix Factorization with Sparseness Constraints,NMFSC)的適用性就更好,主要體現在:一是提取的成分向量可以很好地反映樣本的局部特性;二是可以控制矩陣的稀疏性;三是可以控制矩陣的非負性。另外,考慮到衛星遙測觀測數據序列相關性,本文提出利用前l時刻的觀測數據動態表示當前時刻t的樣本數據,即動態稀疏化非負矩陣分解(Dynamic NMFSC,DNMFSC)。
NMFSC[15]通過對系數矩陣或基矩陣進行稀疏性約束而發展的一種非負矩陣分解算法。
“稀疏”是基于大量數據,而只采用少量數據單位對典型的數據向量進行有效表示。這表征具有稀疏性的數據大部分為零。基于L1和L2范數,NMFSC的稀疏因子定義為

(1)
n為非負向量x的維數。當且僅當所有元素都相等時,稀疏因子為0;當x中只包含一個非零分量時,稀疏因子為1。稀疏因子在[0,1]區間內滑動,稀疏度通過它來表示。sparseness(x)越大,意味著向量越稀疏,反之越稠密。
NMFSC定義為:假定非負矩陣V∈Rm×n,將其分解為兩個非負子矩陣W∈Rm×r和H∈Rr×n,使得目標函數值
E(W,H)=‖V-WH‖2
(2)
最小,并滿足如下約束條件
sparseness(wi)=Sw,?i
(3)
sparseness(hi)=Sh,?i
(4)
wi、hi分別為W的第i列、H的第i行。Sw為基矩陣W期望的稀疏度,Sh為系數矩陣H期望的稀疏度。
基于L1和L2范數,NMFSC控制稀疏度,即在超平面和超球面上的投影,實現對矩陣非負性和稀疏性的控制。
觀測變量存在序列相關性,是衛星遙測數據的特點之一,也就是說當前時刻與過去若干時刻存在關系。考慮到數據的時序相關性,使用前l時刻的觀測數據對此刻t的樣本數據進行擴充,調整如下
V=[VtVt-1…Vt-l]T
(5)
但是這種方法造成輸入數據的維數增加了原來的l倍。因此,引入指數加權滑動平均對衛星數據進行處理,具體為

(6)
式中,wi表示權重,對應變量Vt-(i-1);Vt表示觀測變量V的第t時刻采樣點;N表示滑動窗口的寬度,本文實驗中經驗取值為3。
為了提高DNMFSC算法穩定性,本文提出對DNMFSC進行PCA初始化處理。通過構建統計量H2和SPE確定衛星運行是否有異常。基于PCA-DNMFSC對衛星故障進行在線檢測,其框架如圖1所示,其流程如下:首先利用指數加權滑動平均對衛星正常數據進行動態化,利用PCA初始化DNMFSC,計算得到W和H;其次計算正常數據的統計量H2和SPE,確定統計量的控制限;然后對采集的當前時刻的衛星數據進行動態化表示,利用離線建模中的基矩陣W計算當前時刻的系統矩陣H,并由此計算當前時刻統計量H2和SPE;最后判斷當前時刻統計量是否超出控制限。如果有超出,表明有故障發生。

圖1 基于PCA-DNMFSC的異常檢測框架
與NMFSC類似,初始化H以及W時,DNMFSC也存在隨機性,結果可能各不相同,影響算法的穩定性[16]。本文提出采用PCA對DNMFSC進行初始化,首先使用PCA分解數據集,然后將分解矩陣中負元素調整為零,作為初始化矩陣代入DNMFSC。
標準化Xn*m

(7)


(8)
本文實驗中,設α為0.85。
得分矩陣如下

(9)
Tr初始化NMFSC分解矩陣H。
由于DNMFSC算法的非負性約束,為了能用Tr和Pr對W和H進行初始化,具體為:
P(W)=max(0,W)
(10)
P(H)=max(0,H)
(11)
之后,P(H)初始化NMFSC分解矩陣H,P(W)初始化分解矩陣W。
利用PCA對DNMFSC初始化后,構造基于DNMFSC的故障檢測模型。通過衛星遙測離線正常數據構建的統計量H2和SPE,在線監測當前時刻遙測數據,確定衛星運行是否異常。
統計量H2表征模型內部變化的測度,定義如下

(12)

統計量SPE表征模型外部變化的測度,定義如下

(13)

3.2.1 統計量控制限
DNMFSC檢測中,經過分解所得到的信號不一定服從正態分布,所以統計量的控制限不能通過正態分布的置信區間來確定。針對這種情況,目前研究者大多采用非參數化方法。文獻[17]使用非參數化方法-核密度估計(Kernel Density Estimation,KDE)來確定過程的控制限。本文借鑒此思想來求取控制限,核密度估計為

(14)
x為待求控制限值,xi表示觀測值,n表示樣本數,h表示平滑參數。K為核函數,本文選用高斯核

(15)
則核密度估計為

(16)
3.2.2 統計量累計貢獻
如果檢測到系統發生異常,就要對此異常進行分析,識別故障發生的原因。但當異常發生時,對于時序性數據,單點SPE統計量貢獻值無法良好地描述系統的異常,原因在于監測數據的異常變化不太明顯。基于此,本文將某段時期內貢獻值進行加和,并加上一定的權值,及早確定異常的發生是由哪個變量產生的。SPE統計量累計貢獻定義如下

(17)
w1+w2+w3=1
(18)
其中,a為變量序號,et,a為異常時刻t的預測誤差et的第a變量對應的預測誤差,et-1,a為t前一時刻第a變量對應的預測誤差,et+1,a為t后一時刻第a變量對應的預測誤差。
通過計算系統運行過程中每一個變量的統計量貢獻值并進行作圖,觀測每一個變量的運行情況,從而確定異常由哪些變量引起的。
為了驗證本文所提出的PCA-DNMFSC故障檢測方法的性能,選取兩種在軌衛星實際發生故障進行驗證。
1) 極軌氣象衛星輻射計發生的故障。表1列出與此故障可能相關的遙測數據,用TMC進行標示定義,在某時刻觀測變量可表示為Z=[TMC1 TMC2… TMC19]。微波溫度計某年某月某日世界時 21:21:33轉動機構突然發生卡滯,本來做勻速/變速掃描的探測機構無法正常轉動掃描,定標源溫度無法保持恒定。

表1 衛星遙測變量
2) 靜止氣象衛星輻射計發生的故障。與此故障可能相關的遙測數據有26個,如表2所示,用TMC進行標示定義,在某時刻觀測變量可表示為L=[TMC1 TMC2… TMC26]。輻射計某年某月某日世界時 00:55:15南北鏡正常掃描,東西鏡異常,其位置緩慢移動。

表2 衛星遙測變量
實驗從有效性和合理性兩方面驗證本文所提出的算法。硬件運行環境為AMD Atlonm(tm) X4 750 Quad Core Processor 3.40GHZ 4.00GB,軟件運行環境為Matlab R2015a。本文樣本選取如下:采用衛星載荷未發生故障前一天以世界時 00:00時刻開始的500個正常數據作為訓練樣本,發生故障前后的1000個數據作為測試樣本。極軌衛星微波溫度計測試樣本為世界時 21:10:05時刻開始的1000個樣本數據,靜止衛星輻射計測試樣本為故障發生世界時 00:55時刻開始的1000個數據。
采用本文所提的PCA-DNMFSC方法對以上兩種故障進行檢測,檢測結果如圖2-5所示。圖中,實線表示測試樣本的統計量,虛線是統計量控制限,超過控制限的統計量表明此時刻檢測到衛星故障。圖2-3為極軌衛星微波溫度計故障的檢測結果,圖4-5為靜止衛星輻射計故障的檢測結果。圖2圖4為H2統計量,圖3圖5為SPE統計量。從圖2-3可以看出,極軌衛星微波溫度計從第200個樣本,統計量開始逐漸超過閾值,而第200個樣本時刻為世界時 21:20:35,實際發生故障的時間大約為世界時 21:21:33,符合實際情況。從圖4-5可以看出,靜止衛星輻射計故障從第300個樣本,統計量大多超過閾值,而第300個樣本時刻為世界時00:54:59,這與發生故障的時間大約世界時00:55:15相吻合。這些表明本文所提出的方法能準確地檢測出衛星所發生的故障。

圖2 微波溫度計H2統計量的檢測結果

圖3 微波溫度計SPE統計量的檢測結果

圖4 輻射計H2統計量的檢測結果
為了對異常進行識別并進一步驗證所提出方法的有效性,對變量的累計貢獻率進行分析。圖6-7分別表示極軌衛星微波溫度計和靜止衛星輻射計遙測變量對SPE統計量的累計貢獻程度,其中橫軸表示變量,縱軸表示累計貢獻率大小。由圖6可知,各遙測變量對微波溫度計故障累計貢獻率從大到小依次是:變量12定標源源體溫度、變量6接收機主備機狀態、變量5綜合處理器主備機狀態、變量17冷空觀測開始角度高字節、變量18冷空觀測開始角度低字節。這說明載荷微波溫度計的異常產生的原因,同時也說明了最相關的遙測變量依次為變量12、6、5、17、18。圖7說明了輻射計異常最相關的遙測變量為:輻射計激磁+5V、輻射計激磁-5V、輻射計掃描機構狀態、輻射計EW角度、輻射計NS角度、輻射計感應同步器+15V、輻射計感應同步器-15V、輻射計電機控制+5V以及輻射計電機驅動+8V。這些都與實際情況相符,進一步驗證算法的有效性。

圖5 輻射計SPE統計量的檢測結果

圖6 微波溫度計故障貢獻圖

圖7 輻射計故障貢獻圖
為了實現衛星的智能化健康管理及避免故障漏報,本文提出一種采用PCA-DNMFSC算法進行衛星遙測故障檢測。通過不同衛星故障數據開展方法驗證,實驗證實了以NMFSC為基礎進行故障檢測使局部特征提取更加有效,樣本數據的動態化表示能夠表達變量間的相關性,使用PCA進行初始化解決DNMFSC算法不穩定問題,結論表明算法可以及時、有效、智能地檢測出故障。并且通過統計量累計貢獻率確定產生故障的相關變量與實際相符,這進一步驗證算法的有效性。進一步開展算法研究,使算法更加完善,力求應用到實際業務中。