賀志遠(yuǎn)，呂衛(wèi)民，胡文林

(海軍航空大學(xué)， 山東 煙臺(tái) 264001)

速率陀螺儀是導(dǎo)彈姿態(tài)控制系統(tǒng)的重要組成部分，它控制著飛行過程角速度的輸出，一旦出現(xiàn)故障將影響飛行精度，導(dǎo)致飛行失控甚至彈體自毀。速率陀螺儀作為一種高精密的設(shè)備，內(nèi)部包括眾多電子元器件、傳感器以及機(jī)械連接部分。它的性能隨著工作時(shí)間逐漸改變，同時(shí)，受到環(huán)境因素(溫度、濕度、振動(dòng))的影響，可能出現(xiàn)電路短路、斷路、傳感器失靈以及機(jī)械部件變形等失效。因此，對(duì)速率陀螺儀進(jìn)行可靠性評(píng)估，合理預(yù)測(cè)剩余壽命，對(duì)于保證導(dǎo)彈的安全性具有重大意義。

傳統(tǒng)的可靠性評(píng)估方法主要是根據(jù)失效機(jī)理建立模型進(jìn)行失效分析，但是隨著設(shè)備復(fù)雜程度的提高，很難建立起可靠的的失效物理模型[1]。針對(duì)這一問題，基于性能退化數(shù)據(jù)的分析方法出現(xiàn)，該方法利用設(shè)備性能隨時(shí)間變化過程中的數(shù)據(jù)，建立性能退化模型，計(jì)算過程較為簡(jiǎn)便且結(jié)果具有較高的可靠度，已成為可靠性評(píng)估及壽命預(yù)測(cè)的重要手段之一[2-4]。目前，性能退化模型主要包括隨機(jī)系數(shù)模型和隨機(jī)過程模型。Lu等[5]最早使用隨機(jī)系數(shù)模型研究退化過程，文獻(xiàn)中假定模型的參數(shù)服從某種隨機(jī)分布。這種方法的弊端在于，研究對(duì)象的參數(shù)一旦確定，其退化軌跡就確定，體現(xiàn)不出退化過程的隨機(jī)性。而基于隨機(jī)過程模型的方法解決了這一問題，該模型很好的體現(xiàn)了退化過程的不確定性，更加符合設(shè)備的實(shí)際退化規(guī)律。隨機(jī)過程模型中主要包括Gamma過程和Wiener過程等， Yuan[6]建立Gamma過程退化模型對(duì)核電站某些元件的性能進(jìn)行可靠性評(píng)估。王浩偉等[7]提出利用Gamma過程參數(shù)的非共扼先驗(yàn)分布進(jìn)行Bayesian統(tǒng)計(jì)推斷的剩余壽命預(yù)測(cè)方法。Gamma過程的退化增量是非負(fù)的，且嚴(yán)格單調(diào)遞增，適用于單調(diào)變化的退化過程。而Wiener過程具有不嚴(yán)格單調(diào)的特性，計(jì)算分析特性較好，適用于退化過程具有波動(dòng)性的產(chǎn)品。Park等[8]利用Wiener過程建立模型，對(duì)在加速應(yīng)力條件下得到的退化數(shù)據(jù)進(jìn)行處理。王小林[9]通過Wiener過程隨機(jī)模型預(yù)測(cè)了金屬化膜電容器的剩余壽命。

速率陀螺儀退化過程具有波動(dòng)性，并非單調(diào)過程，本文采用Wiener過程建立隨機(jī)模型進(jìn)行計(jì)算。但是由于現(xiàn)實(shí)工作環(huán)境及設(shè)備本身的復(fù)雜性，退化建模過程還存在以下問題：① 設(shè)備內(nèi)部組件較多，包含了電子部分和機(jī)械部分，退化規(guī)律不同，導(dǎo)致設(shè)備整體退化軌跡呈現(xiàn)非線性。線性Wiener過程已經(jīng)無法有效的研究設(shè)備的退化過程，有些學(xué)者提出將非線性過程通過某種尺度轉(zhuǎn)化為線性過程的方法[10]，王浩偉等[11]通過時(shí)間尺度轉(zhuǎn)化方法將指數(shù)退化模型轉(zhuǎn)化為線性模型對(duì)電連接器的可靠性進(jìn)行了研究。雖然這種方法簡(jiǎn)化了計(jì)算過程，但結(jié)果的精確性也相應(yīng)降低，而且有些退化模型無法轉(zhuǎn)化為線性過程。② 不同于飛機(jī)等航空裝備高頻使用的特點(diǎn)，導(dǎo)彈裝備的壽命更長(zhǎng)，在全壽命周期內(nèi)要經(jīng)歷不同的任務(wù)階段，為了方便研究，分為運(yùn)輸、貯存、戰(zhàn)備三個(gè)階段，每個(gè)階段的環(huán)境因素不完全相同，內(nèi)部設(shè)備退化的主導(dǎo)機(jī)理也在改變，退化趨勢(shì)呈現(xiàn)階段性。針對(duì)多階段退化問題，部分學(xué)者提出建立多階段退化模型。袁慶祥等[12]引入轉(zhuǎn)換函數(shù)改進(jìn)Wiener過程，建立了多階段退化模型預(yù)測(cè)電機(jī)的剩余壽命。該方法有效解決了多階段退化問題，但是引入了新的函數(shù)，計(jì)算難度增大。

針對(duì)速率陀螺儀退化過程的非線性、多階段的特點(diǎn)，本文采用Wiener過程建立非線性退化模型，利用收集到退化數(shù)據(jù)進(jìn)行模型參數(shù)的估計(jì)，然后根據(jù)Kalman濾波算法對(duì)模型進(jìn)行參數(shù)修正，實(shí)現(xiàn)退化過程的多階段性，使剩余壽命預(yù)測(cè)結(jié)果更加精確。

1 剩余壽命預(yù)測(cè)模型

1.1 非線性退化模型

導(dǎo)彈壽命周期經(jīng)歷的環(huán)境因素復(fù)雜，內(nèi)部設(shè)備的真實(shí)退化軌跡很難估計(jì)，也并非時(shí)間的線性函數(shù)，根據(jù)這一退化特點(diǎn)，本文建立基于Wiener過程的非線性退化模型。假設(shè)該設(shè)備退化過程只有一個(gè)關(guān)鍵性能參數(shù)，則非線性退化模型表示如下：

(1)

式中：X(t)為時(shí)刻t的性能退化量，X(0)為初始時(shí)刻的性能退化量(為了體現(xiàn)模型的一般性，令X(0)=0)；μ(t;θ)為非線性函數(shù)，顯然，當(dāng)μ(t;θ)=μ時(shí)，上式轉(zhuǎn)化為一般Wiener過程；σB為擴(kuò)散系數(shù)；B(t)為標(biāo)準(zhǔn)布朗運(yùn)動(dòng)。

1.2 壽命分布

當(dāng)產(chǎn)品的性能退化量達(dá)到預(yù)設(shè)的臨界水平時(shí)，就認(rèn)為此時(shí)產(chǎn)品失效，這一臨界水平稱為失效閾值。假設(shè)該產(chǎn)品的失效閾值為ω，產(chǎn)品的壽命為T，且產(chǎn)品的性能退化軌跡由非線性Wiener過程描述，則產(chǎn)品的壽命T可定義為產(chǎn)品的性能退化量X(t)首次達(dá)到其失效閾值ω的時(shí)間，即

T={t∶X(t)≥ω|X(0)<ω}

(2)

假定μ(t;θ)和σB為固定未知參數(shù)，μ(t;θ)是關(guān)于t的可導(dǎo)函數(shù)，根據(jù)文獻(xiàn)[13]，壽命T的概率密度函數(shù)為：

(3)

其中：

Ε(t;θ)=Λ(t;θ)-tΛ′(t;θ)

給定一個(gè)漂移系數(shù)λ，令μ(t;θ)=λbtb-1(b為固定參數(shù))，壽命T的概率密度函數(shù)可表示為：

(4)

引理如果Z～N(μ,σ2)，A，B，C，ω∈R，那么

(5)

Eλ[fT|λ(t|λ)]

(6)

其中，p(λ)為λ的概率密度函數(shù)。

進(jìn)一步可得，

(7)

可靠度函數(shù)可表達(dá)為：

(8)

1.3 剩余壽命預(yù)測(cè)

假設(shè)特定產(chǎn)品在tn和tn+l時(shí)刻的性能退化量分別為X(tn)和X(tn+l)，按照Wiener過程獨(dú)立增量的性質(zhì)，可得

(9)

令

Z(l)=X(tn+l)-X(tn)=

(10)

顯然式(9)是一個(gè)由布朗運(yùn)動(dòng)B(l)驅(qū)動(dòng)的新的隨機(jī)退化過程。

產(chǎn)品的剩余壽命是指從當(dāng)前時(shí)刻到產(chǎn)品發(fā)生失效時(shí)刻的時(shí)間間隔，根據(jù)壽命的定義可知特定個(gè)體產(chǎn)品在時(shí)刻tn的剩余壽命l可表示為：

L=inf{l∶X(tn+l)≥ω}

(11)

獲得產(chǎn)品剩余壽命的關(guān)鍵是找到剩余壽命的概率密度函數(shù)。按照Wiener過程獨(dú)立增量性質(zhì)，結(jié)合式(9)、式(10)，可得

L=inf{l∶X(tn+l)≥ω}=

inf{l∶X(tn+l)-X(tn)≥ω-X(tn)}=

inf{l∶Z(l)≥ω-X(tn)}

(12)

特定個(gè)體產(chǎn)品在時(shí)刻tn的剩余壽命l，即隨機(jī)過程退化量Z(l)首次達(dá)到閾值ω-X(tn)的時(shí)間。

結(jié)合式(3)、式(11)，可得tn時(shí)刻剩余壽命l的概率密度函數(shù)：

(13)

根據(jù)全概率公式，由式(13)可得，可靠度函數(shù)為式(14)：

(14)

2 參數(shù)估計(jì)

2.1 EM算法

上述非線性退化模型的未知參數(shù)設(shè)為Θ=(μλ,σλ,b,σB)。其中，b和σB為固定系數(shù)，描述各階段退化過程的共性。μλ和σλ為隨機(jī)系數(shù)，通過修正來描述不同階段退化的過程的特性。

假設(shè)可測(cè)退化數(shù)據(jù)的設(shè)備有n臺(tái)，分別在時(shí)刻t1，t2，t3，…，tm對(duì)設(shè)備進(jìn)行測(cè)試，得到第i臺(tái)設(shè)備在時(shí)刻j的性能退化數(shù)據(jù)Xi={X(t1,i),X(t2,i),…,X(tj,i)}，相應(yīng)的退化增量為ΔX(tj,i)，因此

ΔX(tj,i)=X(tj,i)-X(tj-1,i)=

(15)

其中，λi為漂移系數(shù)的真實(shí)值。

未知參數(shù)增加為Θ=(λi,μλ,σλ,b,σB)。

根據(jù)式(16)以及最優(yōu)解搜索算法，可以得到未知參數(shù)的極大似然估計(jì)值，但計(jì)算過程較為復(fù)雜且未知參數(shù)數(shù)量較多，針對(duì)這一問題，本文采用Gibbs抽樣算法來估計(jì)未知參數(shù)。

2.2 Gibbs采樣算法

Gibbs采樣是馬爾科夫蒙特卡洛(Markov Chain Monte Carlo，MCMC)方法的一種，用于在難以直接采樣時(shí)從多變量概率分布中抽取樣本，尤其在高維參數(shù)的情形下，Gibbs采樣有更大的優(yōu)勢(shì)。Gibbs采樣算法實(shí)質(zhì)上是M-H抽樣時(shí)，接受概率等于1的特例，但由于其原理簡(jiǎn)單且適用范圍廣泛，因此被認(rèn)為是一種特殊的抽樣算法。

Gibbs采樣算法的步驟如下：

1) 給定未知參數(shù)λi,μλ,σλ,b,σB的聯(lián)合分布函數(shù)和先驗(yàn)分布函數(shù)，分別為p(λi,μλ,σλ,b,σB)和p(λi),p(μλ),p(σλ),p(b),p(σB)；

3) 計(jì)算出未知參數(shù)λi的條件概率密度，可表示為

(17)

同式(17)，分別計(jì)算其他未知參數(shù)的條件概率密度；

(18)

6) 使用收斂后的函數(shù)q(λi,μλ,σλ,b,σB)替代p(λi,μλ,σλ,b,σB)，重復(fù)3)、4)過程k次，得到k組Gibbs樣本。

最后，利用k組Gibbs樣本和貝葉斯公式計(jì)算未知參數(shù)的后驗(yàn)分布，進(jìn)一步得到未知參數(shù)的估計(jì)值。

3 基于Kalman濾波算法的多階段參數(shù)修正

Kalman濾波是一種高效的遞歸濾波器，建立在隱馬爾科夫模型的基礎(chǔ)上。它利用目標(biāo)的動(dòng)態(tài)信息，去掉噪聲干擾，得到關(guān)于目標(biāo)位置的估計(jì)。這個(gè)估計(jì)包括：濾波(對(duì)當(dāng)前位置的估計(jì))、預(yù)測(cè)(對(duì)將來位置的估計(jì))和平滑(對(duì)過去位置的估計(jì))[15]。

速率陀螺儀的退化過程具有明顯的多階段性，為了更加準(zhǔn)確地預(yù)測(cè)剩余壽命，通過Kalman濾波算法對(duì)各階段模型參數(shù)進(jìn)行修正[16]，根據(jù)式(4)可知，修正對(duì)象為漂移系數(shù)λ。

假設(shè)在時(shí)刻tn進(jìn)行漂移系數(shù)的修正，根據(jù)式(1)并結(jié)合Kalman濾波算法原理，選取狀態(tài)空間方程如下：

(19)

多種色釉于一坯胎中也是在近現(xiàn)代才流行起來，多種色釉的窯變比單色釉可說難，也可說易。難是難在對(duì)釉料料性的掌握，那種色釉流動(dòng)，那種不流動(dòng)，色釉與色釉之間能否結(jié)合，產(chǎn)生的效果如何，幾種色釉結(jié)合當(dāng)怎樣的比例才能燒制出最佳效果等等，這就是難。說它容易，也只能說出現(xiàn)二次創(chuàng)作以后陶瓷藝術(shù)家才敢如此認(rèn)為。顏色釉基礎(chǔ)上的二次創(chuàng)作是根據(jù)藝術(shù)家對(duì)色釉窯變效果來巧妙裝飾的，主觀性較強(qiáng)，但選取的畫面需要和色釉窯變氣氛相吻合，才能達(dá)到異曲同工之妙。

誤差項(xiàng)νδ可能是正值也可能是負(fù)值，這將影響到λn的正負(fù)，根據(jù)文獻(xiàn)[17]，本文研究對(duì)象的漂移系數(shù)非負(fù)，為此確定了一個(gè)新的狀態(tài)空間方程：

(20)

根據(jù)Kalman濾波算法及文獻(xiàn)[14]，引入衰減系數(shù)γ(ti)。

γ(ti)=max{γ(ti),1}

(21)

且

其中，其中η為遺忘因子，0.95<η<1。

(22)

且

Pi|i-1=γ(ti)Pi-1|i-1

更新方差Pi|i。

(23)

最后，當(dāng)i=n時(shí)，完成時(shí)刻tn漂移系數(shù)的階段更新，修正后的退化模型可表示為

(24)

4 實(shí)例分析

通過對(duì)速率陀螺儀進(jìn)行故障模式、機(jī)理及影響分析，發(fā)現(xiàn)壽命周期內(nèi)其性能指標(biāo)受到溫度、濕度及振動(dòng)等環(huán)境因素的影響。通過查閱相關(guān)文獻(xiàn)，本文選取零偏值(°/h)作為性能退化指標(biāo)。試驗(yàn)對(duì)象為6臺(tái)速率陀螺儀，為了模擬該設(shè)備的工作特點(diǎn)，驗(yàn)證多階段非線性模型的準(zhǔn)確性，選取不同的溫度應(yīng)力(0～1 000 h，20 ℃；1 000～2 000 h，55 ℃；2 000～3 192 h，70 ℃)進(jìn)行退化試驗(yàn)，試驗(yàn)過程共進(jìn)行20次數(shù)據(jù)采集，每次間隔168 h，測(cè)試數(shù)據(jù)如表1所示。其中前4臺(tái)測(cè)試數(shù)據(jù)用于模型參數(shù)估計(jì)，后2臺(tái)測(cè)試數(shù)據(jù)用于多階段參數(shù)的修正。試驗(yàn)的失效閾值設(shè)置為0.5(°/h)。

從圖1可以看出不同溫度應(yīng)力下，零偏測(cè)試數(shù)據(jù)的退化存在著明顯的差異，與上文多階段退化的假設(shè)相符。

利用測(cè)試數(shù)據(jù)，通過Gibbs采樣算法計(jì)算得到非線性退化模型的未知參數(shù)，再將計(jì)算結(jié)果與極大似然估計(jì)方法(EM)結(jié)果一起，列于表2。采用Gibbs采樣算法，模型的擬合殘差均方差(MSE)較小，說明該方法估計(jì)的參數(shù)更加符合實(shí)際數(shù)據(jù)的趨勢(shì)。

表1 速率陀螺儀零偏測(cè)試數(shù)據(jù) (°/h)

圖1 零偏測(cè)試數(shù)據(jù)

參數(shù)獲取方法測(cè)試參數(shù)μλσ2λbσBMSEEM0.0022 90.0007 160.9080.025 61.56Gibbs0.002 370.000 4251.020.018 50.873

再運(yùn)用Kalman濾波方法對(duì)各階段退化模型的參數(shù)進(jìn)行修正，修正后的參數(shù)如表3所示，并通過AIC準(zhǔn)則與單階段非線性退化模型參數(shù)進(jìn)行比較，AIC值越小說明參數(shù)更加符合模型。

AIC=-2(maxl)+2p

(25)

其中，maxl表示最大似然函數(shù)的值；p表示模型未知參數(shù)的個(gè)數(shù)。

表3 不同階段模型參數(shù)的修正值

從表3可以看出多階段退化模型參數(shù)AIC的平均值小于單階段退化模型，證明本文的多階段參數(shù)估計(jì)方法更加精確，退化模型更加可靠。

最后利用修正后的模型進(jìn)行剩余壽命預(yù)測(cè)，如圖2所示。選取工作時(shí)間3 000 h，計(jì)算多階段非線性模型和單階段非線性模型的平均剩余壽命，分別為2 362 h和2 721 h。而在退化試驗(yàn)中，速率陀螺儀的平均失效時(shí)間為5 016 h。通過對(duì)比，本文所采用的基于Wiener過程的多階段非線性退化模型所預(yù)測(cè)的剩余壽命更加接近實(shí)際情況。

圖2 剩余壽命預(yù)測(cè)概率密度

5 結(jié)論

1) 本文采用基于非線性Wiener過程建立速率陀螺儀的退化模型，使用Gibbs采樣算法估計(jì)復(fù)雜模型的多個(gè)未知參數(shù)，并通過Kalman濾波算法進(jìn)行漂移系數(shù)的修正，實(shí)現(xiàn)了退化過程的多階段性。

2) 退化試驗(yàn)結(jié)果與不考慮漂移系數(shù)修正的單階段退化模型相比，多階段退化模型更加符合實(shí)際退化過程，對(duì)于剩余壽命預(yù)測(cè)更加準(zhǔn)確。

3) 由于非實(shí)驗(yàn)條件下采集到的數(shù)據(jù)少且無效數(shù)據(jù)較多，使Gibbs采樣算法計(jì)算精度降低。如何在樣本少、干擾多的情況下估計(jì)模型參數(shù)值，需進(jìn)一步研究。