馬長剛，李 青，陳 明，孫大林

（1.空軍勤務學院航空彈藥系，江蘇 徐州 221000；2.解放軍94676 部隊，上海 202150）

0 引言

通過某裝備檢測參數(shù)可以判斷出該裝備的質(zhì)量狀態(tài)，但測試參數(shù)種類眾多，數(shù)量龐大，其中包含了眾多表征該裝備質(zhì)量狀態(tài)的冗余測試參數(shù)。在進行其質(zhì)量狀態(tài)評估、故障預測、檢測周期優(yōu)化等內(nèi)容時，無論是準確度還是計算的效率，均受到了數(shù)量龐大的測試參數(shù)群的影響。完全可從大量參數(shù)群中篩選出具有代表性的，同時與該裝備質(zhì)量狀態(tài)相關(guān)性較強的一部分測試參數(shù)，去除大部分冗余檢測參數(shù)，作為評估裝備質(zhì)量狀態(tài)、預測裝備故障等的狀態(tài)特征量主成分，方便后續(xù)基于狀態(tài)檢測參數(shù)的問題研究。

1 某裝備質(zhì)量狀態(tài)參數(shù)時間序列分析

假設(shè)某裝備在k 個測試時刻其質(zhì)量狀態(tài)測試的測試參數(shù)進行收集，共有m 種狀態(tài)測試參數(shù)，用Zt表示t 時刻檢測參數(shù)值向量，zit表示在時刻t 質(zhì)量狀態(tài)測試參數(shù)的測試值，則Zt=（z1t，z2t，…，zmk）是一個k 維的隨機列向量，而且假設(shè)每種質(zhì)量狀態(tài)檢測參數(shù)所有的歷史測試數(shù)據(jù)個數(shù)為n。由于該裝備檢測參數(shù)數(shù)據(jù)量比較大，如果考慮到當前t 時刻之前該裝備質(zhì)量狀態(tài)檢測所有歷史檢測時刻的檢測參數(shù)，龐大的數(shù)據(jù)量會給計算帶來很大麻煩，而且由于疊加效應，有可能會導致計算結(jié)果偏離實際值的情況。所以選取該裝備質(zhì)量狀態(tài)檢測特征量參數(shù)的關(guān)鍵是確定延遲期p，解決到底將歷史檢測時刻往后推到哪個時刻比較合理這一問題。利用時序分析法中多維自回歸（AR）［1］確定該裝備狀態(tài)檢測信息在時序上的延遲期p，具體步驟如下。

1.1 數(shù)據(jù)預處理

要利用自回歸模型確定延遲期p，首先需要保證不同測試參數(shù)項在同一時間序列中數(shù)據(jù)值的均值為零，對均值不為零的原始質(zhì)量狀態(tài)檢測數(shù)據(jù)做如下處理：

一階差分之后，利用自協(xié)方差式（3）判斷是否滿足平穩(wěn)序列要求：

按照一階跟二階的差分方法繼續(xù)對非平穩(wěn)序列進行差分處理，直至其滿足平穩(wěn)化條件。平穩(wěn)化處理的目的就是消除質(zhì)量狀態(tài)檢測參數(shù)在時間序列上的規(guī)律性和線性特點［3］。

1.2 確定AR 模型的階數(shù)

由式（1）變換計算后得到該裝備質(zhì)量狀態(tài)檢測參數(shù)的平穩(wěn)時間序列，建立此時k 維質(zhì)量狀態(tài)參數(shù)變量的AR 模型［4］。

根據(jù)以上定義，建立該裝備質(zhì)量狀態(tài)檢測參數(shù)的AR（p）模型：

式（6）可以表示為：

式（8）可以用下面的式子來表示：

要求解式（9）中的最小值，借助Q 關(guān)于Φ 的偏導來求解，令偏導為0，即：

即：

通過該裝備質(zhì)量狀態(tài)檢測參數(shù)的零均值處理、平穩(wěn)化處理確定了時間序列延遲期p，也確定了其AR（p）模型，并對模型中的參數(shù)進行了估計，但確定的時間序列延遲期p 是否滿足要求，是否精確，估計值是否合理，需要對其進行檢驗。在這里，根據(jù)估計值，取顯著性水平α 進行檢驗。

可通過分析殘差序列來檢驗模型的合適與否，采用Ljung 和Box 改進并提出的LB 統(tǒng)計量檢驗法進行檢驗：

假定原假設(shè)H0：模型延遲期p 符合要求；

H1：模型延遲期p 不符合要求；

2 某裝備質(zhì)量狀態(tài)特征量選取

主成分分析法（PCA）利用降維的思想，把多指標轉(zhuǎn)化為少數(shù)幾個綜合指標，即主成分，其中每個主成分都能夠反映原始變量的大部分信息，且所含信息互不重復［6］。PCA 方法中僅僅考慮到了狀態(tài)監(jiān)測信息間的互相關(guān)關(guān)系，而動態(tài)主成分DPCA 不僅考慮到了狀態(tài)檢測信息間的互相關(guān)而且考慮到同一狀態(tài)檢測信息不同時序時刻間的自相關(guān)［7］。

在求得模型階次p 的基礎(chǔ)上，考慮到該裝備狀態(tài)檢測信息間既互相關(guān)又在不同時序上自相關(guān)，采用動態(tài)主成分分析法提取狀態(tài)檢測信息的主成分。步驟如下：

1）計算協(xié)方差矩陣

通過時序模型的計算檢驗，確定模型的階次為p，狀態(tài)檢測信息間的協(xié)方差矩陣由互相關(guān)協(xié)方差矩陣C（0），以及自相關(guān)協(xié)方差矩陣C（1），C（2），…，C（p）組成，此時的協(xié)方差矩陣為一個由（p+1）×（p+1）個m×m 矩陣塊組成的協(xié)方差矩陣，即C。每一矩陣塊表示為：

i，j=1，2，…，p+1，Zit表示t 時刻狀態(tài)檢測信息i的檢測值，表示其檢測值的平均值。

2）確定相關(guān)矩陣

3）計算相關(guān)矩陣的特征值

4）提取主成分

主成分的提取，通過主成分的綜合貢獻率進行確定，根據(jù)步驟3）中的特征值的大小排序，前k 個主成分的綜合貢獻率［8］為：

根據(jù)綜合貢獻率要求就可確定出k 個主成分。

3 實例分析

本章前面部分已經(jīng)詳細地介紹了該裝備質(zhì)量狀態(tài)特征量選取的方法，方法合適與否，都需要根據(jù)真實的數(shù)據(jù)進行驗證。從一線保障部隊收集到了部分該裝備壽命周期內(nèi)完整的測試數(shù)據(jù)。用本文提出來的將時間序列和主成分分析法結(jié)合起來的動態(tài)主成分分析法來提取這幾枚該裝備質(zhì)量狀態(tài)特征量。選擇收集到的某型號該裝備完整歷史測試數(shù)據(jù)，進行特征量的選取，該裝備總共8 次歷史測試次數(shù)的原始測試數(shù)據(jù)如表1 所示，表中C1～C90分別對應該裝備各項參數(shù)。

表1 測試參數(shù)原始數(shù)據(jù)

根據(jù)表1 中該裝備原始的8 次歷史狀態(tài)檢測數(shù)據(jù)，下面利用前面一節(jié)中提到的該裝備質(zhì)量狀態(tài)特征量參數(shù)選取方法提取該裝備的狀態(tài)特征量。具體實現(xiàn)步驟如下。

3.1 原始數(shù)據(jù)預處理

表1 中收集了該裝備從出廠后到最近一次測試的數(shù)據(jù)，數(shù)據(jù)鏈比較完整，首先可以選擇其中的一部分參數(shù)，從數(shù)據(jù)層面來了解一下隨著檢測次數(shù)的增加每一項測試參數(shù)數(shù)據(jù)值的變化趨勢，在這里以表1 中的部分測試參數(shù)在8 次測試中數(shù)據(jù)值的變化情況為例，畫出其變化曲線圖，如圖1 所示：

圖1 原始參數(shù)數(shù)據(jù)變化趨勢

圖1 中選擇了90 多項測試參數(shù)中的3 項，以它們?yōu)槔?，分析其原始?shù)據(jù)的變化趨勢，從圖中可以看出，每項參數(shù)總共8 次測試的均值都不為零，數(shù)據(jù)波動起伏比較大，并沒有嚴格圍繞0 值線規(guī)律地上下波動，因此，數(shù)據(jù)值不滿足模型求解的零均值要求，需要對其進行零均值處理。

收集到的該裝備的測試參數(shù)項目總共90 項，如果在后面故障預測以及質(zhì)量狀態(tài)評估的過程中，考慮所有測試參數(shù)，將這么多的數(shù)據(jù)代入進行計算的話，計算起來比較困難，而且計算的準確度也受到很大的影響，前面已經(jīng)分析過了，這90 項測試參數(shù)中，有很大一部分參數(shù)之間是相互關(guān)聯(lián)的，它們的變化具有同步性。由圖1 中每一項參數(shù)變化情況的分析中，其變化不具備零均值和平穩(wěn)化的要求。根據(jù)特征量選取模型計算的要求，首先，需要對原始數(shù)據(jù)零均值、平穩(wěn)化地處理。

1）數(shù)據(jù)零均值處理

2）數(shù)據(jù)平穩(wěn)化處理

在滿足零均值的情況下，根據(jù)平穩(wěn)化處理的方法，結(jié)合式（2），對1）中零均值后的數(shù)據(jù)進行平穩(wěn)化處理，并利用式（3）對其平穩(wěn)化效果進行檢驗，最終可以得到滿足后面質(zhì)量狀態(tài)特征量選取模型的零均值、平穩(wěn)化后的數(shù)據(jù)。經(jīng)過編寫MATLAB 程序進行計算，平穩(wěn)化處理后的數(shù)據(jù)滿足平穩(wěn)化要求。圖1中3 項參數(shù)在8 次測試中的參數(shù)值經(jīng)平穩(wěn)后處理后的變化趨勢圖如圖2 所示。

圖2 平穩(wěn)化后參數(shù)數(shù)據(jù)變化趨勢

3.2 延遲時序P 的確定

根據(jù)建立的該裝備質(zhì)量狀態(tài)檢測參數(shù)AR（p）的模型，要確定AR 模型的階數(shù)p，需要求出式（7）、式（10）中的所有測試參數(shù)在P 個時刻內(nèi)完整的測試數(shù)據(jù)矩陣X 和當前時刻質(zhì)量狀態(tài)參數(shù)測試數(shù)據(jù)矩陣Z。算法實現(xiàn)步驟如下：

Step1：首先假定延遲時序P=1；

Step2：根據(jù)p 的值，確定出矩陣X、Z；

Step6：如果通過檢驗則算法結(jié)束，認為p 值符合要求，若檢驗不通過，則令p=p+1，則跳轉(zhuǎn)至Step2進入循環(huán)，直至p 的值通過檢驗，算法結(jié)束。

根據(jù)上面的算法步驟，編寫MATLAB 代碼進行計算，經(jīng)計算，當P=2 時，通過檢驗，因此，最終確定p 的值為2。

3.3 某裝備質(zhì)量狀態(tài)特征量的選取

確定出了時間序列模型的延遲時序P，并檢驗了P 值的合理性，運用主分量分析的方法提取出該裝備的質(zhì)量狀態(tài)特征量，提取步驟及計算過程如下。

首先，根據(jù)前面所確定的延遲時序P 為2，所以在考慮同一參數(shù)不同檢測時刻測試值間相互關(guān)系時考慮到往前推移2 個時刻，則本例中，根據(jù)式（14）計算協(xié)方差矩陣，這里計算出的協(xié)方差矩陣是由9 個90×90 維的矩陣組成的高維矩陣，即：

C11，…，C33分別為上面提到的9 個90×90 維的協(xié)方差矩陣，然后分別求出協(xié)方差矩陣所對應的相關(guān)矩陣R，并求解出相關(guān)矩陣R 的特征值和特征向量。由于協(xié)方差陣和相關(guān)陣維數(shù)較高，具體矩陣在這里不作展示。本例中，最終求解出的大部分特征值極小，這里僅列出部分特征值，如表2 所示。

表2 部分特征值

對所有特征值進行排序，并選取前100 項作出特征值散點圖進行分析特征值分布規(guī)律，如3圖所示。

圖3 特征值分布規(guī)律

從圖3 中可以看出來，特征值從第7 項開始已經(jīng)趨近于零，直觀判斷，可以將前6 項特征值所對應的質(zhì)量狀態(tài)測試參數(shù)作為我們要選取的狀態(tài)特征量。但根據(jù)前面建模過程中的方法，應從其特征值綜合貢獻率出發(fā)來計算，前7 項特征值綜合貢獻率的和為96.12%，滿足對該裝備質(zhì)量狀態(tài)特征量選取綜合貢獻率大于95%的要求。因此，最終選取出來的該裝備質(zhì)量狀態(tài)特征量為6 項，分別是特征值排序中前6 項所對應的參數(shù)項。這6 項參數(shù)分別對應表1 中的C76、C77、C78、C79、C80、C81。

4 結(jié)論

針對某裝備測試參數(shù)眾多，冗余參數(shù)較多，在利用參數(shù)進行該裝備故障預測或者狀態(tài)評估時計算很不方便，預測和評估的準確度受該裝備狀態(tài)參數(shù)維度高的影響較大。因此，借助時間序列分析確定該裝備質(zhì)量狀態(tài)不同時序條件下的自相關(guān)、互相關(guān)關(guān)系，又利用動態(tài)主分量分析法從眾多測試參數(shù)中選取出決定該裝備質(zhì)量狀態(tài)的部分主分量，選取效果非常明顯。