葉 松，袁艷艷

（北京航天自動控制研究所，北京，100854）

0 引 言

慣性器件利用硬件冗余的方式，可以提高運載火箭制導(dǎo)系統(tǒng)的可靠性，如果某個慣性器件出現(xiàn)故障，需要對冗余慣性測量組合進行故障檢測與隔離。

基于主成分析法（Principal Component Analysis，PCA）的故障檢測與隔離技術(shù)已經(jīng)得到了國內(nèi)外學(xué)者的廣泛關(guān)注，其實質(zhì)是一種向量空間變換，將多變量空間轉(zhuǎn)化為子空間，該子空間可以在最小維數(shù)下保持原空間的最大方差。從變量優(yōu)化的角度會發(fā)現(xiàn)子空間維數(shù)小于原空間的維數(shù)，且相互獨立。Feng Xu等人提出了一種基于區(qū)間觀測器和不變集的傳感器故障檢測與隔離方法[1]，故障檢測用于實時故障報警，故障隔離用于檢測瞬態(tài)性能，仿真結(jié)果表明了該方法的有效性。Mania Navi等人研究了主成分析法用于自主水下航行器傳感器的故障檢測與隔離[2]，仿真結(jié)果表明，該方法能用于非線性系統(tǒng)。Walton R. Williamson等人提出了一種用在深空探測衛(wèi)星上的故障檢測與隔離算[3]法，該方法需要通過 Shiryayev序列概率比來測試并計算故障發(fā)生的可能性。Kim等人做了開創(chuàng)性研究，提出了兩種冗余系統(tǒng)的故障檢測方法[4]，一種是結(jié)合小波變換的奇偶空間法，另一種是利用閾值確定故障時間間隔并基于模型的故障檢測與隔離方法，通過擴展卡爾曼濾波器來估計殘差作為檢測故障的參考值，仿真結(jié)果表明：這兩種方法都能用于冗余系統(tǒng)的故障診斷。Alvarez利用主成分析法對廢水處理廠的機械進行故障診斷[5]。Daniel Skoogh提出了一種能同時檢測不同振幅故障的檢測與隔離方法[6]，傳感器信號的測量矢量被投影到配置矩陣的正交補空間，通過觀察投影點的分布來進行故障診斷。Oliveira在故障檢測的數(shù)據(jù)處理過程中加入中值濾波來消除噪聲的影響[7]，仿真結(jié)果表明該方法能夠檢測出微小的階躍故障信號。Sérgio Brás提出利用狀態(tài)觀測器提供的狀態(tài)信息來確定不同模型的故障[8]，然而該方法需要箭載計算機具有強大的數(shù)據(jù)處理能力。張昱君和劉愛倫針對化工生產(chǎn)過程中故障提出了一種改進的主元分析方法[9]，可以檢測出對于過程影響較小的故障。考慮到傳統(tǒng)的主成分析法對訓(xùn)練數(shù)據(jù)集之外的數(shù)據(jù)異常敏感，Yvon Tharrault建立了一個穩(wěn)健的主成分析模型并使用多故障檢測和隔離的結(jié)構(gòu)化殘差[10]。主成分析法應(yīng)用在冗余慣組故障檢測研究較少。近年來，Lee and Park提出了一種結(jié)合了奇偶向量法的改進主成分析法[11]，該方法能夠檢測動態(tài)環(huán)境下的故障信息，但不能進行故障隔離。

本文提出了一種基于廣義似然比的改進主成分析法，廣義似然算法用來進行故障檢測，主成分析法用來進行故障隔離?？紤]到加入濾波器能夠降低系統(tǒng)的噪聲，因此在故障檢測前先對奇偶向量進行濾波預(yù)處理，目的是能夠檢測到小幅值故障并具有故障隔離能力。數(shù)學(xué)仿真結(jié)果表明：該方法能夠檢測小幅值故障并提高故障檢測與隔離的效果。

1 經(jīng)典主成分析法

主成分析法可以被定義為原始數(shù)據(jù)到新數(shù)據(jù)的線性變換，并且這些新的數(shù)據(jù)集相互線性無關(guān)。眾所周知，傳統(tǒng)的主成分析法進行故障檢測的方法是將統(tǒng)計量Q和 T2分別與其閾值進行比較，當(dāng)統(tǒng)計量的數(shù)值超過了控制限，就可以判定過程中出現(xiàn)了異常情況。然而這樣僅僅能判定在這個時間點出現(xiàn)了故障，而無法判定是哪個部分的異常導(dǎo)致。能夠幫助確定哪個部分出現(xiàn)故障工具就是貢獻圖，它通過每個過程變量對統(tǒng)計量的貢獻值來確定。

1.1 T 2統(tǒng)計量

T2統(tǒng)計量反映的是主成分在變化趨勢和變化幅值上與模型的偏差程度，計算公式如下：

式中 ti為測量樣本的第i個主元得分；λi為協(xié)方差矩陣的第i個特征值；Xi為第i時刻的采樣值；為主元模型的載荷矩陣T2為統(tǒng)計量，是多個變量共同累加的標量，表示數(shù)據(jù)點到模型中心的距離（建模是經(jīng)過標準化后中心在原點），因此它的大小反映了新數(shù)據(jù)偏離正常狀態(tài)的程度。

如果實際協(xié)方差未知，則T2統(tǒng)計量服從F分布，T2統(tǒng)計量閾值為

當(dāng)一個過程變量的信息沒有被很好的反映在主元模型中，這個變量出現(xiàn)數(shù)據(jù)異常情況后通過 T2統(tǒng)計量無法檢測到故障，因此在這種情況下需要通過分析測量數(shù)據(jù)的殘差來計算Q統(tǒng)計量。

1.2 Q統(tǒng)計量

Q統(tǒng)計量是平方預(yù)測誤差，它代表測量值與主元模型的偏差程度。Q統(tǒng)計量的分布如下：

式中 ?P為由最小二乘法計算得到的載荷矩陣；e為由剩余載荷矩陣生成的誤差矢量。

1.3 載荷矩陣計算

載荷矩陣算法也稱非線性迭代部分最小二乘算法（Non-linear Iterative Partial Least Squares，NIPALS）。

a）從X中任選一列 Xj，并記為t1，即t1= Xj；

b）計算 P1： P

c）將 P1的長度歸一化；

e）將步驟a中的t1與步驟d作比較，如果他們幾乎一樣，則算法收斂，計算停止，如果它們不一樣，回到步驟a；

f） X = X - t1p1T，得到新的誤差矩陣，計算其他主元即可。

1.4 主元個數(shù)選擇方法

采用基于平均特征值的累計方差貢獻率法來確定主元個數(shù)。累計方差貢獻率的定義：矩陣X的協(xié)方差矩陣的前k個特征值的和除以它的所有特征值的和（特征值要按由大到小順序排列），它表示了前k個主元所解釋的數(shù)據(jù)變化占全部數(shù)據(jù)變換的比例。

式中 control為人為設(shè)定的控制限，一般選取 85%，此時滿足條件的最大主元個數(shù)記為 k1。

接著，計算相關(guān)系數(shù)矩陣所有特征值的均值，選取大于平均特征值作為主元特征值，同時舍棄掉那些小于均值的特征值，對應(yīng)的最小的主元特征值的序號記為 k2。

最后，取 k1， k2兩者中的最大值為主元個數(shù)，即

2 改進的主成分析法

傳統(tǒng)主成分析法只能檢測出較大幅值故障，本文提出將最大廣義似然法應(yīng)用到傳統(tǒng)的主成分析法中，目的是檢測小幅值故障。方法的具體流程如圖 1所示（以典型的正六面體結(jié)構(gòu)為例）。

圖1 故障檢測與隔離流程Fig.1 Flowchart of Fault Detection and Isolation

步驟1：通過安裝測量矩陣計算解耦矩陣V，采集正常工作訓(xùn)練數(shù)據(jù)TZ，得到訓(xùn)練等價向量TP，并標準化后計算計算投影向量W。

步驟 2：應(yīng)用是將傳感器采集得到的數(shù)據(jù)Z，計算得到等價向量P，再通過投影向量得到Ty。

步驟3：對Ty進行濾波處理。

步驟 4：計算統(tǒng)計量，根據(jù)統(tǒng)計量的貢獻值進行故障檢測。

2.1 廣義似然法

冗余捷聯(lián)慣組含有n個傳感器，測量輸出只包含測量噪聲，那么其測量方程為

因此，利用矩陣 V可得到等價向量P = V Z = Vε+Vf， P ∈ Rn-3。即等價向量與被測量X是線性無關(guān)的，只與誤差ε和故障f相關(guān)。

根據(jù)等價空間原理，解耦矩陣V位于安裝測量矩陣H的左零空間，V的秩為 n -3。為了確定V的元素，Potter和Suman建議將V陣選擇為具有正對角元的上三角陣，然后通過正交化，便可確定V陣的元素。具體的算法為

2.2 數(shù)字濾波

流程中對Ty進行濾波預(yù)處理。眾所周知低通濾波器可以降低系統(tǒng)噪聲帶來的不利影響，而且可以對外來數(shù)據(jù)不敏感，因此設(shè)計的低通濾波器的傳遞函數(shù)為

式中1T為濾波器的時間常數(shù)。

將式（5）進行Z變換可得：

對式（6）左右進行簡化，可以得到一階濾波器[16]的最終表達式為

濾波器的參數(shù)影響濾波器性能的關(guān)鍵因素，如何選取一個合適的濾波器參數(shù)是非常重要的。一方面濾波器參數(shù)反映了系統(tǒng)的響應(yīng)速度，另一方面反映的是系統(tǒng)對于高頻信號的衰減，因此濾波器參數(shù)的選取是一個需要折中的問題。文獻[16]里面的相關(guān)研究，選取的一階濾波器的參數(shù)為 0.434。此時濾波器的差分方程如下：

3 仿真分析

3.1 仿真條件

分別對兩種配置結(jié)構(gòu)（正六面體結(jié)構(gòu)和三正交二斜置[15]結(jié)構(gòu)）進行仿真分析。

陀螺儀的誤差參數(shù)如表1所示。

表1 慣組誤差Tab.1 Error Parameter

仿真條件如下：在400 s的時候給第5支陀螺儀加入1（°）/h的常值故障。

3.2 PCA方法仿真

3.2.1 正六面體

利用圖 1的流程，對正六面體冗余慣組進行仿真分析，結(jié)果如圖2、圖3所示。

圖2 特征平面的投影點的分布Fig.2 Projection of Feature Plane

圖3 各個變量對統(tǒng)計量Q的貢獻值Fig.3 Contribution of Each Variable to Statistics

3.2.2 三正交二斜置

同樣三正交二斜置的慣組配置進行了仿真，結(jié)果如圖4、圖5所示。

圖4 特征平面的投影點的分布Fig.4 Projection of Feature Plane

圖5 各個變量對統(tǒng)計量Q的貢獻值Fig.5 Contribution of Each Variable to Statistics

從圖2和圖4可以看出：待檢測數(shù)據(jù)和正常數(shù)據(jù)在特征平面上的投影沒有完全分開，說明原PCA方法檢測出1（°）/h幅值的故障困難。并且投影點在X軸的位置和Y軸的位置隨時間變化曲線可以看出有故障發(fā)生，但是不明顯，然而圖3和圖5表明故障隔離效果較好。

3.3 改進PCA算法仿真

3.3.1 正六面體

為了檢驗加入濾波器之后的故障檢測性能，得到的仿真結(jié)果如圖6、圖7所示。

圖6 特征平面上的投影Fig.6 Projection of Feature Plane

圖7 各個變量對統(tǒng)計量Q的貢獻值Fig.7 Contribution of Each Variable to Statistics

3.3.2 三正二斜

三正交二斜置的慣組配置的仿真結(jié)果如圖8、圖9所示。

圖8 特征平面上的投影Fig.8 Projection of Feature Plane

圖9 各個變量對統(tǒng)計量Q的貢獻值Fig.9 Contribution of Each Variable to Statistics

圖6 代表的是投影位置，可以看出過程中是否發(fā)生了故障，可看出正常數(shù)據(jù)的投影點和待檢測數(shù)據(jù)的投影點分開，說明已經(jīng)檢測到了故障，沒有完全分開是因為故障只發(fā)生了一段時間而非整個過程（在400 s處發(fā)生故障）。從圖7可以看出第5個傳感器對統(tǒng)計量Q的貢獻值最大，因此可以斷定第5個傳感器發(fā)生了故障。圖3和圖7有細微的差異（各個變量對統(tǒng)計量Q的貢獻量數(shù)值不完全一致），這是由于陀螺儀誤差項包含有隨機誤差，而故障隔離的結(jié)果是正確的。同理，對于三正交二斜置的安裝方式也有相同的結(jié)論。

4 結(jié) 論

傳統(tǒng)的主成分析法應(yīng)用在動態(tài)系統(tǒng)的故障診斷時，會存在故障檢測效果略差但故障隔離效果很好的問題。本文提出了一種基于廣義似然比算法的改進主成分分析法，并使用濾波器進行預(yù)處理，新方法能夠克服傳統(tǒng)主成分析法只能檢測大幅值故障的不足。