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

（1.空軍勤務(wù)學(xué)院航空彈藥系，江蘇 徐州 221000；2.解放軍94676 部隊(duì)，上海 202150）

0 引言

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

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

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

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

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

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

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

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

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

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

式（6）可以表示為：

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

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

即：

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

可通過分析殘差序列來檢驗(yàn)?zāi)Ｐ偷暮线m與否，采用Ljung 和Box 改進(jìn)并提出的LB 統(tǒng)計(jì)量檢驗(yàn)法進(jìn)行檢驗(yàn)：

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

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

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

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

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

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

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

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

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

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

4）提取主成分

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

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

3 實(shí)例分析

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

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

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

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

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

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

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

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

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

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

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

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

3.2 延遲時(shí)序P 的確定

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

Step1：首先假定延遲時(shí)序P=1；

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

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

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

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

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

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

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

表2 部分特征值

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

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

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

4 結(jié)論

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