考慮隨機(jī)漂移-擴(kuò)散的非線性維納過(guò)程可靠性分析

中圖分類號(hào)：TB9;TH133.33

文獻(xiàn)標(biāo)志碼：A 文章編號(hào)：1674-5124（2025）06-0056-11

Nonlinear Wiener process degradation modeling considering random drift-diffusion

XIAO Meng1， SHEN Ao12， XU Zhibiaol， SHAN Susu1， XIN Mingjiang1 （1. School ofRail Transit， Wuyi University， Jiangmen 529020，China; 2. Tongcheng County Metrology Verification and Testing Institute， Xianning 437499， China）

Abstract： Individual difference in the process of degradation is one of the difficulties in degradation modeling and the main source of prediction uncertainty. The common random drift model considers the variability of degradation rate，butregards the diffusionparameter asa constant，and ignores the actual individual diferences in the degradation fluctuation，which is diffcult to adequately describe the heterogeneity in the degradation process，and thus affects the reliability evaluationofthe equipment.To solve this problem，a reliabilityanalysis method based on random drift-difusion Wiener process is proposed， which can consider the heterogeneity of degradation rate and degradation fluctuation at the same time，and can select the appropriate distribution according to the information in the actual degradation data，which improves the versatility and flexibility of modeling. Aiming at the increasing complexity of the proposed model， a Bayesian estimation algorithm based on Markov chain Monte Carlo（MCMC） method is proposed， which can estimate the model parameters at one time.The accuracy of the proposed algorithm is verified by data simulation experiments，and the efectiveness of the proposed model is verified by comparison experiments between turbofan engine data and laser data. Keywords：Wiener process;MCMC; reliability： heterogeneity

0 引言

在工業(yè)生產(chǎn)應(yīng)用中，設(shè)備和系統(tǒng)的可靠性會(huì)隨著時(shí)間的推移而降低。一些正在服役的大型機(jī)械設(shè)備如果發(fā)生故障將會(huì)帶來(lái)嚴(yán)重的經(jīng)濟(jì)損失甚至人員傷亡[1]。大量研究表明在設(shè)備發(fā)生故障前會(huì)表現(xiàn)出異常的信號(hào)特征，例如軸承的振動(dòng)加速度增大[2-3]，以及鋰電池容量的減少[4-5]等。為了設(shè)備安全生產(chǎn)和預(yù)測(cè)性維護(hù)通常人為的設(shè)定失效閾值，一旦設(shè)備的性能指標(biāo)超過(guò)了閾值則認(rèn)為設(shè)備失效。同時(shí)為了預(yù)測(cè)故障時(shí)間并評(píng)估設(shè)備的可靠性，需要監(jiān)測(cè)這些性能指標(biāo)并進(jìn)行建模。由于同一批產(chǎn)品的制造材料以及外部環(huán)境的不同，導(dǎo)致同一批產(chǎn)品中的不同個(gè)體的退化軌跡各不相同，而呈現(xiàn)明顯的個(gè)體差異。這種退化數(shù)據(jù)中蘊(yùn)含的個(gè)體差異給退化建模帶來(lái)了困難，并會(huì)導(dǎo)致壽命預(yù)測(cè)中的不確定性[。常見(jiàn)的隨機(jī)漂移維納過(guò)程模型難以描述這種異質(zhì)性，在眾多退化模型中維納過(guò)程有著直觀的物理解釋性和優(yōu)良的數(shù)學(xué)性質(zhì)，而廣泛應(yīng)用于退化建模中[7]。一般的非線性維納過(guò)程退化模型為：

式中： X（t_k） ?t_k 時(shí)刻設(shè)備退化狀態(tài)；（204號(hào) k —觀測(cè)次數(shù)；（204號(hào) λ ——漂移參數(shù)，表征退化速率；（204號(hào) μ（τ;θ） ——參數(shù) θ 的時(shí)變非線性函數(shù)，用于描述設(shè)備退化狀態(tài)的非線性；B（t_k） ！ σ ———標(biāo)準(zhǔn)布朗運(yùn)動(dòng)、擴(kuò)散參數(shù)，分別用來(lái)描述退化過(guò)程中的時(shí)變隨機(jī)波動(dòng)及其波動(dòng)程度。

不同函數(shù)形式 μ（τ;θ） 可以描述不同形式的退化過(guò)程。例如當(dāng) 時(shí)退化為線性退化過(guò)程， 時(shí)可用來(lái)表示冪函數(shù)退化過(guò)程。在后續(xù)的研究當(dāng)中采用冪函數(shù)的退化過(guò)程來(lái)表征非線性退化。

在維納退化模型中常見(jiàn)的假設(shè)，將λ服從正態(tài)分布，而將 σ 設(shè)定為常數(shù)以表征個(gè)體的異質(zhì)性[。然而在實(shí)際退化的過(guò)程中， σ 也呈現(xiàn)明顯的個(gè)體差異，忽略 σ 中的異質(zhì)性，將會(huì)影響建模性能[8]。例如，在加速退化壽命實(shí)驗(yàn)中的不同水平下，實(shí)驗(yàn)樣本的擴(kuò)散參數(shù)也會(huì)發(fā)生變化，為了解釋這一實(shí)驗(yàn)現(xiàn)象YE等[在常見(jiàn)的退化模型假設(shè)中進(jìn)行了創(chuàng)新，將λ和 σ 視作確定的比例關(guān)系。在文獻(xiàn)[6]與文獻(xiàn)[9]中都假設(shè)λ服從正態(tài)分布，但利用正態(tài)分布對(duì)非均勻的λ存在著一定的缺點(diǎn)。第一個(gè)不足在于，通常在同一批設(shè)備當(dāng)中，往往都是正值或者負(fù)值，但根據(jù)正態(tài)分布的性質(zhì)，其在全值域都取值[10]。此外正態(tài)分布具有對(duì)稱性，當(dāng)λ服從一個(gè)偏峰分布時(shí)，基于λ服從正態(tài)分布的假設(shè)將會(huì)影響預(yù)測(cè)的準(zhǔn)確性[11]。WANG等[12]假設(shè)漂移參數(shù)服從廣義逆高斯分布（GeneralizedInverseGaussian，GIG）以及與 σ 為確定的比例關(guān)系，克服了服從正態(tài)分布假設(shè)的包含負(fù)值域和對(duì)稱分布的不足。雖然ZHAI等[13]通過(guò)單獨(dú)擬合激光退化數(shù)據(jù)集[14]中的各路徑，發(fā)現(xiàn)λ與σ 存在正相關(guān)，但是并非完全正相關(guān)。本文利用類似的方法對(duì)砷化鎵激光器數(shù)據(jù)集進(jìn)行了相同實(shí)驗(yàn)，估計(jì)的個(gè)體參數(shù)λ和 σ 散點(diǎn)圖如圖1所示。從圖1可以看出， λ 與 σ 并不是完全正相關(guān)，也無(wú)法通過(guò)合適的比例系數(shù)來(lái)描述λ與 σ 之間的相關(guān)性。因此， λ 與 σ 為比例關(guān)系的假設(shè)，可能會(huì)降低模型預(yù)測(cè)的準(zhǔn)確性。

鑒于以上問(wèn)題，本文提出一種考慮隨機(jī)漂移-擴(kuò)散的非線性維納過(guò)程退化分析方法。與確定比例關(guān)系模型不同的是，本文 σ 的隨機(jī)性與 λ 不相關(guān)。從文獻(xiàn)[10當(dāng)中可分析，對(duì)于不同的數(shù)據(jù)集退化速率與退化波動(dòng)可能服從不同的分布。對(duì)此，本文通過(guò)最大似然來(lái)選擇兩個(gè)隨機(jī)參數(shù)的分布，這種考慮有助于提升建模通用性和靈活性。當(dāng)然考慮 σ 的隨機(jī)性，增加了模型的復(fù)雜程度，針對(duì)所提模型的復(fù)雜度增加，提出一種基于MCMC方法的模型參數(shù)貝葉斯估計(jì)算法。通過(guò)數(shù)據(jù)仿真實(shí)驗(yàn)驗(yàn)證了所提算法的準(zhǔn)確性，利用渦扇發(fā)動(dòng)機(jī)、激光數(shù)據(jù)比較實(shí)驗(yàn)，驗(yàn)證了所提模型的有效性。

1模型描述

根據(jù)引言的描述，將式（1）中 λ 和 σ 當(dāng)作隨機(jī)參數(shù)，以描述雙隨機(jī)參數(shù)所導(dǎo)致的退化異質(zhì)性。假設(shè)λ～π（α_λ，β_λ），σ～π（α_σ，β_σ） ，其中 π（?） 為未知分布。令總體參數(shù) ?={α_λ，β_λ，α_σ，β_σ，θ} ，假設(shè)有 N 條退化軌跡，則退化軌跡的個(gè)體退化參數(shù) {λ_i，σ_i} 為來(lái)自總體分布p（λ，σ|?） 的獨(dú)立同分布，即可定義為：

{λ_i，σ_i}～_i.i.d.p（λ，σ|?），i=1，2，…，N

基于上述的描述，令本文模型為 M0 。為了體現(xiàn)所提模型的優(yōu)越性，設(shè)定三組對(duì)比模型，來(lái)比較模型的擬合性能以及預(yù)測(cè)性能。將文獻(xiàn)[12]所提的與 σ 成比例的模型定義為M1，即 λ=κσ 。將λ設(shè)定為隨機(jī)參數(shù)， σ 為固定參數(shù)，設(shè)定為M2。將文獻(xiàn)[6]所提模型定義為M3，即 λ 服從正態(tài)分布以及 σ 為固定參數(shù)。

給定由式（1定義的退化過(guò)程，基于首達(dá)時(shí)間（firsthittingtime，F(xiàn)HT）[15]的概念，設(shè)備的壽命 T 可表示為：

T=inf{t_kgt;0，X（t_k）?w}

其中 w 為失效閾值，基于式（1）在給定 λ 和 σ 的情況下，壽命T的邊緣概率密度可以表示為：

其中 ，由于本文考慮了 λ 與 σ 的隨機(jī)性，故對(duì) λ 和 σ 求雙重積分，壽命 T 的邊緣概率密度可定義為：

因式 （5）中涉及多維積分運(yùn)算，一般情況下很難得到解析解。采用蒙特卡羅積分的方法近似求解，從先驗(yàn)分布 p（λ，σ|?） 中抽取 K 組樣本 （λ^（1），σ^（1）），… （λ^（m），σ^（m）） ，則式（5）可近似為：

根據(jù)式（6）壽命 T 的累計(jì)概率密度可近似為：

2 貝葉斯參數(shù)估計(jì)

概率模型參數(shù)的估計(jì)問(wèn)題可通過(guò)極大似然方法或貝葉斯方法來(lái)實(shí)現(xiàn)，具體選擇涉及模型復(fù)雜度，樣本情況以及基于估計(jì)參數(shù)的后續(xù)統(tǒng)計(jì)推理等因素。相比于極大似然方法，貝葉斯方法能直接提供參數(shù)的后驗(yàn)分布，有利于分析參數(shù)估計(jì)的不確定性和后續(xù)在線觀測(cè)數(shù)據(jù)融合[16]。

為方便表述，全部的退化狀態(tài)為 x_i=[x_i，1，x_i，2，… x_i，niJ^T ，由式（1）中給出的條件高斯性質(zhì)，在給定模型參數(shù) {λ_i，σ_i，θ} 條件下，退化狀態(tài)增量 Δx_i，k= x_i，k-x_i，k-1 服從正態(tài)分布：

其中 Δt_i，k=t_i，k-t_i，k-1 表示相鄰觀測(cè)之間的間隔時(shí)間，k=1，2，…，n_i 。根據(jù)Wiener過(guò)程的獨(dú)立增量性質(zhì)和退化觀測(cè)的條件獨(dú)立假設(shè)，對(duì)于第i個(gè)設(shè)備的退化數(shù)據(jù)，其完全似然函數(shù)表達(dá)式為：

因 N 個(gè)設(shè)備的退化觀測(cè)數(shù)據(jù) {x₁，x₂，…，x_N} 是獨(dú)立樣本，則模型全部未知參數(shù) P={λ_i，σ_i，?} 的后驗(yàn)分布為：

式（10）表示的聯(lián)合后驗(yàn)分布較為復(fù)雜，難以通過(guò)積分方法得到參數(shù) ? 的邊緣后驗(yàn)分布。在后續(xù)的研究當(dāng)中借助MCMC方法來(lái)隨機(jī)模擬其邊緣后驗(yàn)分布。

先驗(yàn)分布設(shè)置是貝葉斯分析的重要步驟，根據(jù)上一節(jié)的模型描述對(duì)于個(gè)體參數(shù) λ_i 和 σ_i 都服從未知分布，對(duì)此先要確定其服從的分布。本文給出了4種不含非負(fù)值域以及偏峰的備選分布，分別為伽瑪、對(duì)數(shù)正態(tài)、韋布爾和指數(shù)分布，根據(jù)最大似然對(duì)備選分布的組合進(jìn)行選擇。進(jìn)一步指定參數(shù) ? 的先驗(yàn)分布為：

其中C為半柯西分布，其概率密度函數(shù)為：

式中： m 分布峰值位置的位置參數(shù);

b ——最大值一半處的一半寬度的尺度參數(shù)。

在貝葉斯當(dāng)中參數(shù)集{n，ξ，max，baχ，mβ，bβ，maσ，b_ασ，m_βσ，b_βσ} 可視為 ? 的先驗(yàn)分布中的超參數(shù)。超參數(shù)通常是基于先驗(yàn)知識(shí)或是以弱信息先驗(yàn)形式予以設(shè)置。由式（1）、（2）和（11）構(gòu)成分層貝葉斯退化過(guò)程模型可概括為圖2。陰影節(jié)點(diǎn)表示可觀測(cè)變量，虛線圈節(jié)點(diǎn)代表超參數(shù)，其他節(jié)點(diǎn)為未知參數(shù)或變量，有向邊代表節(jié)點(diǎn)間的統(tǒng)計(jì)相依關(guān)系，圓角矩形框表示對(duì)參數(shù)或變量的獨(dú)立觀測(cè)數(shù)目。

MCMC方法通過(guò)構(gòu)造一條以聯(lián)合后驗(yàn)分布p（?|x₁，x₂，…，x_N） 為平穩(wěn)分布的遍歷馬爾科夫鏈?_j，j=1，2，… ，并從中模擬出后驗(yàn)分布的隨機(jī)樣本。由于MCMC方法對(duì)后驗(yàn)分布期望的近似具有相合性和漸近性正態(tài)等優(yōu)良性質(zhì)，使其成為對(duì)后驗(yàn)分布總體的合理近似[17]。

在MCMC方法當(dāng)中Gibbs采樣[18]和MetropolisHasting算法[19]最為常用，但兩種算法都存在著局限性。從圖2可見(jiàn)，本文的模型參數(shù)具有較強(qiáng)的耦合性，然而Gibbs采樣將會(huì)割裂這種性質(zhì)，進(jìn)而影響參數(shù)估計(jì)的準(zhǔn)確性。此外Metropolis-Hasting算法的抽樣效率將會(huì)影響模型參數(shù)估計(jì)的運(yùn)行時(shí)間。對(duì)此本文采用HMC（HamiltonianMonteCarlo，HMC）算法框架來(lái)進(jìn)行后驗(yàn)分布的隨機(jī)抽樣。HMC是基于漢密爾頓抽樣的MCMC方法，其原理模擬動(dòng)力學(xué)當(dāng)中的跳點(diǎn)法 （leap frogmethod）更新位置參數(shù)[20]?；贖MC算法計(jì)算模型參數(shù)的迭代運(yùn)算步驟如下。

第一步：根據(jù)聯(lián)合后驗(yàn)分布在漢密爾頓動(dòng)力系統(tǒng)下的勢(shì)能，即

第二步：設(shè)置蛙跳算法迭代過(guò)程的次數(shù)L（l=1，2，…，L） 以及步長(zhǎng) δ ，確定馬爾科夫鏈的長(zhǎng)度為 A 以及預(yù)熱長(zhǎng)度S，確定參數(shù)的初始值 ?⁰= {λ_i⁰，σ_i⁰，α_λ⁰，β_λ⁰，α_α⁰，β_α⁰，θ⁰}_c

第三步：從標(biāo)準(zhǔn)正態(tài)分布 ΔN（0，1） 生成隨機(jī)數(shù)確定動(dòng)量變量 P⁰

第四步：根據(jù)第二步的初始值以及步長(zhǎng)為 δ 進(jìn)行蛙跳算法迭代 L 次，更新得到狀態(tài) ?^* 和 P^* ○

第五步：計(jì)算接受概率：

其中， K（P） 為動(dòng)能對(duì)應(yīng)函數(shù)， K（P）=P²/2 0

第六步：從均勻分布 U（0，1） 產(chǎn)生隨機(jī)數(shù) ζ ，若ζlt;φ ，接受備選狀態(tài)作為下次循環(huán)的初始值，否則舍棄，重復(fù)第三步到第六步，直至 A 次。即可得到一系列參數(shù)值 ?_j 組成的馬爾科夫鏈。

第七步：舍棄預(yù)熱期的抽樣，運(yùn)用蒙特卡羅方法計(jì)算參數(shù)的后驗(yàn)估計(jì)。

HMC算法的流程圖如圖3所示，在本文后續(xù)的研究當(dāng)中使用Rstan軟件包來(lái)實(shí)現(xiàn)HMC算法[21-22]。

基于此，本文的方法結(jié)構(gòu)圖如圖4所示。

3實(shí)驗(yàn)研究

為驗(yàn)證本文隨機(jī)退化建模方法的合理性，本節(jié)先通過(guò)一組仿真退化數(shù)據(jù)來(lái)評(píng)價(jià)參數(shù)估計(jì)算法的準(zhǔn)確性。再利用渦扇發(fā)動(dòng)機(jī)退化數(shù)據(jù)和激光器退化數(shù)據(jù)來(lái)對(duì)比分析本文模型M0以及參照模型M1、M2和M3的數(shù)據(jù)擬合與預(yù)測(cè)性能。

3.1 人工仿真實(shí)驗(yàn)

通過(guò)一組仿真數(shù)據(jù)來(lái)進(jìn)行研究，首先通過(guò)人工指定模型參數(shù)，并隨機(jī)生成不同規(guī)模的退化狀態(tài)序列。然后，使用HMC算法對(duì)參數(shù)進(jìn)行估計(jì)。通過(guò)將參數(shù)估計(jì)與真實(shí)值進(jìn)行比較，來(lái)評(píng)估算法性能[23]。

過(guò)30，因此可設(shè)定的超參數(shù)為：

假設(shè) λ～G（α_λ，β_λ） 和 σ～G（α_σ，β_σ） 。式中G（·）為伽瑪分布， α_σ 為形狀參數(shù)， β_σ 為尺度參數(shù)。指定模型參數(shù) ? 為：

α_λ=5，β_λ=2，α_σ=1，β_σ=0.2，θ=3

隨機(jī)生成具有相同觀測(cè)次數(shù)和相同觀測(cè)時(shí)間間隔的三組不同規(guī)模 N=20，30，50 的退化狀態(tài)及其觀測(cè)數(shù)據(jù)，其中觀測(cè)次數(shù)均為 n_i=100 ，觀測(cè)時(shí)間步長(zhǎng)為 Δt_i，k=0.01s 。 N=50 組的仿真退化數(shù)據(jù)觀測(cè)序列軌跡如圖5所示。通過(guò)觀察退化軌跡，最大值不超

在此次實(shí)驗(yàn)當(dāng)中的迭代次數(shù)為10000，其中5000步迭代為馬爾科夫鏈的預(yù)熱，實(shí)際的后驗(yàn)估計(jì)為后5000步。模型參數(shù) ? 的后驗(yàn)抽樣診斷情況如圖6所示。圖6給出了 ? 在10000次迭代抽樣的樣本軌跡，均值遍歷和自相關(guān)圖。

由圖6可見(jiàn)，在前5000次預(yù)熱抽樣后，模型參數(shù) ? 的抽樣自相關(guān)函數(shù)（autocorrelationfunction，ACF逐漸接近于0，并且遍歷均值以及抽樣軌跡都處于平穩(wěn)狀態(tài)。由此可以判斷馬爾科夫鏈?zhǔn)諗浚貜?fù)此實(shí)驗(yàn)500次，然后取每次實(shí)驗(yàn)當(dāng)中的后5000次抽樣。實(shí)驗(yàn)的統(tǒng)計(jì)結(jié)果如表1所示。

從表1的估計(jì)結(jié)果來(lái)看，不同樣本規(guī)模下參數(shù)估計(jì)的 95% 置信區(qū)間（confidence interval，CI） 內(nèi)都包含了真實(shí)值，因此驗(yàn)證了模型參數(shù)估計(jì)算法的準(zhǔn)確性。隨著樣本數(shù)量的增加，參數(shù)估計(jì)效果得到了總體提升。本小節(jié)的仿真實(shí)驗(yàn)代碼和數(shù)據(jù)可從網(wǎng)頁(yè)（github.com/superAaaaaaa/simulation-code）下載。

3.2 渦扇發(fā)動(dòng)機(jī)數(shù)據(jù)實(shí)驗(yàn)

此小節(jié)應(yīng)用了NASA提供的渦扇發(fā)動(dòng)機(jī)退化數(shù)據(jù)集[24]。該數(shù)據(jù)集是使用美國(guó)宇航局陸軍研究實(shí)驗(yàn)室開(kāi)發(fā)的商業(yè)模塊化航空推進(jìn)系統(tǒng)進(jìn)行的，從系統(tǒng)輸出中選擇21個(gè)特征來(lái)表征發(fā)動(dòng)機(jī)的退化過(guò)程。該數(shù)據(jù)集描述了工程場(chǎng)景中渦扇發(fā)動(dòng)機(jī)在不同工況和故障模式的退化特征，并被廣泛應(yīng)用[25-27]。

為了反映發(fā)動(dòng)機(jī)的退化過(guò)程，采用多傳感器數(shù)據(jù)融合方式來(lái)構(gòu)建健康系數(shù)（healthindicator，HI）。使用FD002訓(xùn)練集作為此次實(shí)驗(yàn)的樣本，此訓(xùn)練集共有260臺(tái)渦扇發(fā)動(dòng)機(jī)的信號(hào)特征，從中選取200臺(tái)作為神經(jīng)網(wǎng)絡(luò)的訓(xùn)練，然后將剩下的60臺(tái)作為神經(jīng)網(wǎng)絡(luò)的測(cè)試數(shù)據(jù)集，其中測(cè)試數(shù)據(jù)集為下文當(dāng)中研究的對(duì)象。首先通過(guò)式（18）構(gòu)建目標(biāo) 值，然后將200臺(tái)訓(xùn)練集的特征作為輸入，目標(biāo)HI值為輸出，訓(xùn)練神經(jīng)網(wǎng)絡(luò)的參數(shù)[25]。有關(guān)神經(jīng)網(wǎng)絡(luò)特征提取的結(jié)構(gòu)見(jiàn)圖7，參數(shù)以及超參數(shù)見(jiàn)表2。

AIC=-2log-LF+2l

BIC的表達(dá)式為：

BIC=-2log-LF+llnn

t_k?T_i5%

其中l(wèi)og-LF（log-likelihood）為對(duì)數(shù)似然，由Rstan軟件包輸出的lp獲得； n 為數(shù)據(jù)量，l為模型參數(shù)的數(shù)目。AIC和BIC準(zhǔn)則通過(guò)參數(shù)數(shù)目懲罰項(xiàng)來(lái)防止建模的過(guò)參數(shù)化，以平衡模型復(fù)雜度與擬合性能，其取值越小代表模型性能越優(yōu)。

與人工仿真實(shí)驗(yàn)不同的是，在實(shí)際數(shù)據(jù)實(shí)驗(yàn)當(dāng)中縮小了 b 的取值，這樣有利于抽樣的效率。對(duì)于MO的超參數(shù)設(shè)置如下：

其中 T_i 為第i個(gè)發(fā)動(dòng)機(jī)的全壽命周期，由于神經(jīng)網(wǎng)絡(luò)輸出的 曲線存在著較大的波動(dòng)，這將影響到剩余壽命預(yù)測(cè)的精度，因此通常采用移動(dòng)平均法來(lái)解決這個(gè)問(wèn)題。其平滑方式如式（18）所示。

其中 s 為滑動(dòng)窗口的長(zhǎng)度，在此次實(shí)驗(yàn)當(dāng)中 s=15 。為了讓所構(gòu)建的 HI_kⁱ 更具有真實(shí)性，將平滑后的 HI_kⁱ 進(jìn)行歸一化處理，即 t_k=0 時(shí) HI_kⁱ 為 0 。不失一般性，將200臺(tái)渦扇發(fā)動(dòng)機(jī)當(dāng)中末端 HI_kⁱ 最小值的前（s+1）/2 項(xiàng)作為失效閾值，即 w=0.627 。圖8為所構(gòu)建的 HI_kⁱ 訓(xùn)練集以及測(cè)試集。

為比較四種模型對(duì)渦扇發(fā)動(dòng)機(jī)退化數(shù)據(jù)的擬合性能，采用赤池信息準(zhǔn)則（Akaikeinformationcriterion，AIC）和貝葉斯信息準(zhǔn)則（Bayesianinformationcriterion，BIC）作為性能測(cè)度，AIC的表達(dá)式為：

在M1當(dāng)中， κ 的先驗(yàn)分布為 κ～C（m_κ，b_κ） ，超參數(shù)設(shè)置如下

在M2當(dāng)中， σ 的先驗(yàn)分布為 σ～C（m_σ，b_σ） ，超參數(shù)設(shè)置如下

在M3當(dāng)中， λ 的先驗(yàn)分布為 λ～N（α_λ，β_λ²） 。在貝葉斯分析中，常將Normal-Gamma分布作為正態(tài)分布的均值和精度參數(shù)的共軛先驗(yàn)，及將伽瑪分布作為正態(tài)分布的精度參數(shù)的共軛先驗(yàn)[28]。故M3超先驗(yàn)分布以及先驗(yàn)分布設(shè)置為 ，超參數(shù)設(shè)置如下：

η=0，ξ=1，m_αλ=0，b_αλ=1，

m_βλ=m_σ=2，b_βλ=b_σ=0.1

采用所構(gòu)建的渦扇發(fā)動(dòng)機(jī) HI_kⁱ 為實(shí)際退化數(shù)據(jù)，分別用模型M1、M2、M3和本文模型M0進(jìn)行退化過(guò)程建模。首先根據(jù)備選分布進(jìn)行模型先驗(yàn)分布的選擇，估計(jì)的結(jié)果如表3所示。

從表3估計(jì)的結(jié)果，對(duì)數(shù)正態(tài)分布的似然最大，對(duì)此將λ與 σ 的先驗(yàn)分布設(shè)置為對(duì)數(shù)正態(tài)。根據(jù)選擇的結(jié)果對(duì)模型參數(shù)進(jìn)行估計(jì)，四種模型的結(jié)果如表4所示。

由表4可見(jiàn)，比較對(duì)數(shù)似然、赤池信息準(zhǔn)則和貝葉斯信息準(zhǔn)則，模型M0對(duì)該實(shí)際退化數(shù)據(jù)的擬合性能優(yōu)于模型M1、M2和M3。同時(shí)比較表4中M1和M0的 {α_λ，β_λ，θ} 參數(shù)估計(jì)差異可見(jiàn)，是否考慮退化波動(dòng)的隨機(jī)性會(huì)對(duì)所建模型產(chǎn)生較大的影響。這種影響會(huì)進(jìn)而導(dǎo)致式（6和式（7）對(duì)失效預(yù)測(cè)的性能差異。

依據(jù)式（6和式（7）以及表4中估計(jì)的結(jié)果可得出4個(gè)模型的概率密度曲線，以及累計(jì)概率密度曲線，如圖9和圖10所示。值得注意的是，在圖10中的經(jīng)驗(yàn)累計(jì)概率密度為60組實(shí)驗(yàn)樣本首次達(dá)到閾值的循環(huán)次數(shù)。由圖9可見(jiàn)對(duì)于不同的模型假設(shè)，其概率密度曲線還是有差異的，差異較小的M0和M2的原因可能在于，構(gòu)建的渦扇發(fā)動(dòng)機(jī) HI_kⁱ 的擴(kuò)散系數(shù)差異不大。導(dǎo)致概率密度曲線和累計(jì)概率密度曲線并沒(méi)有太大的差異，但從表4的擬合程度上來(lái)分析，其擬合程度要優(yōu)于M2。

由圖10對(duì)比M0和M1可見(jiàn)，M0基本完全處于 95% 的執(zhí)行區(qū)間當(dāng)中，相對(duì)于M1在經(jīng)驗(yàn)失效時(shí)間的分布擬合得更優(yōu)。同理對(duì)比M0和M3其擬合精度也優(yōu)于M3。從表4、圖9和圖10當(dāng)中對(duì)比M2和M3。M2的擬合效果和預(yù)測(cè)效果要優(yōu)于M3，因此也能夠證明文獻(xiàn)[12]當(dāng)中所提出的模型假設(shè)是成立的，同時(shí)為了突顯考慮負(fù)值域與不考慮負(fù)值域的差異，由圖9和圖10分別給出了M3的兩種情況。從結(jié)果上看，不考慮負(fù)值域的M3對(duì)于退化數(shù)據(jù)的預(yù)測(cè)效果更好。

3.3 激光數(shù)據(jù)實(shí)驗(yàn)

為了驗(yàn)證所提模型的普適性，此小節(jié)應(yīng)用了一組激光器的退化數(shù)據(jù)。激光器的工作電流隨運(yùn)行時(shí)間增大而逐漸變大，反映了激光器的性能退化程度。該實(shí)際退化數(shù)據(jù)包括15個(gè)激光器的工作電流隨運(yùn)行時(shí)間變化的百分率數(shù)據(jù)，其中觀測(cè)時(shí)間間隔為250h ，觀測(cè)時(shí)長(zhǎng)為 4000h 。全部激光器的退化觀測(cè)軌跡如圖11所示。該退化數(shù)據(jù)集來(lái)自穩(wěn)定工況下激光器的老化試驗(yàn)，并廣泛應(yīng)用于設(shè)備可靠性和剩余壽命預(yù)測(cè)研究領(lǐng)域[8-9，1-13，23]。文獻(xiàn)[14]假設(shè)激光器的失效閾值為10，（即工作電流增加 10% 時(shí)失效）。此實(shí)驗(yàn)亦選用同樣的失效閾值，即 w=10 。

延用3.2節(jié)的先驗(yàn)分布選擇，即 λ 與 σ 的先驗(yàn)分布為對(duì)數(shù)正態(tài)。分別用4種模型進(jìn)行激光器的退化過(guò)程建模。4種模型的估計(jì)結(jié)果如表5所示。從表5的估計(jì)結(jié)果看，所提模型具有更好的擬合性能。

根據(jù)參數(shù)估計(jì)結(jié)果，代入到式（6）和（7）中。4種模型的概率密度曲線還有累計(jì)概率密度曲線如圖12和13所示。

由圖12可見(jiàn)， λ 選擇正態(tài)分布其模型的概率密度曲線與 λ 選擇對(duì)數(shù)正態(tài)分布的概率密度曲線有著明顯的差異，并且從表5上可知模型M3并沒(méi)有M2擬合性能優(yōu)。因此驗(yàn)證了對(duì)于不同的數(shù)據(jù)集λ可能服從不同的分布。

圖13中所有模型的累計(jì)概率密度曲線都位于 95% 置信邊界之內(nèi)，因此驗(yàn)證了算法估計(jì)的有效性。

4結(jié)束語(yǔ)

由于實(shí)際退化過(guò)程的復(fù)雜性和隨機(jī)性，來(lái)自同類型設(shè)備中不同個(gè)體的退化過(guò)程呈現(xiàn)顯著的異質(zhì)性。其中λ和 σ 的變異是導(dǎo)致這種個(gè)體差異的原因之一。雖然已有研究驗(yàn)證了與 σ 成比例關(guān)系能夠描述這種異質(zhì)性，但本文通過(guò)估計(jì)激光數(shù)據(jù)集的個(gè)體參數(shù)發(fā)現(xiàn)，比例關(guān)系模型并不能充分的描述個(gè)體之間的差異。對(duì)此，本文提出了一種考慮隨機(jī)漂移-擴(kuò)散的非線性維納過(guò)程可靠性分析方法。仿真數(shù)據(jù)、渦扇發(fā)動(dòng)機(jī)退化數(shù)據(jù)和激光退化數(shù)據(jù)實(shí)驗(yàn)結(jié)果驗(yàn)證了該方法的可行性和有效性。主要結(jié)論如下：

1）所提模型同時(shí)考慮了λ和 σ 的變異情形，能夠進(jìn)一步表達(dá)和量化多元異質(zhì)性的影響。相比于只考慮λ異質(zhì)因素的模型，它更加通用。

2）基于分層先驗(yàn)結(jié)構(gòu)，通過(guò)HMC算法實(shí)現(xiàn)了模型的參數(shù)估計(jì)。該算法能有效地實(shí)現(xiàn)具有多個(gè)隨機(jī)參數(shù)的復(fù)雜模型的統(tǒng)計(jì)推理。仿真數(shù)據(jù)實(shí)驗(yàn)證明了該算法的有效性。

3）渦扇發(fā)動(dòng)機(jī)以及激光器數(shù)據(jù)實(shí)驗(yàn)表明，相對(duì)于確定比例系數(shù)模型以及單隨機(jī)模型，所提模型具有更好的擬合性能和預(yù)測(cè)精度。

4）對(duì)于不同的數(shù)據(jù)集，λ和 σ 可能服從不同的分布，依據(jù)最大似然來(lái)選擇合適的分布.使得相對(duì)于一般的退化模型，它更加具有靈活性。

觀測(cè)誤差是影響模型的擬合性能以及預(yù)測(cè)精度的一種重要因素，在后續(xù)的研究當(dāng)中引入觀測(cè)誤差，將會(huì)更進(jìn)一步提升模型的性能。通過(guò)實(shí)驗(yàn)?zāi)M或工程現(xiàn)場(chǎng)收集更多的退化數(shù)據(jù)，提升所提模型的工程適用性也是本文后續(xù)研究方向。

參考文獻(xiàn)

[1]金光.基于退化的可靠性技術(shù)：模型、方法及應(yīng)用[M].北 京：國(guó)防工業(yè)出版社，2014.

[2]薛林，王豪，王云森，等.基于Autoformer的滾動(dòng)軸承剩余使 用壽命預(yù)測(cè)[J].電子測(cè)量技術(shù)，2023，46（13）：169-175. XUEL，WANGH，WANGYS，etal.Remaininguseful life prediction of rolling bearing based on Autoformer[J]. Electronic Measurement Technology， 2023， 46（13）： 169-175.

[3]LI Y， HUANG X， DING P，et al. Wiener-based remaining useful life prediction of rolling bearings using improved Kalman filtering and adaptive modification [J].Measurement， 2021，182.

[4]趙沁峰，蔡艷平，王新軍.鋰電池在不同放電區(qū)間下的剩余 壽命預(yù)測(cè)[J].中國(guó)測(cè)試，2023，49（3）：159-165. ZHAO Q F， CAI Y P， WANG X J. Rmaining useful life prediction of lithium battery under different discharge intervals[J].China Measurement amp; Test，2023，49（3）：159- 165.

[5]JIN G， MATTHEWS D， ZHOU Z. A Bayesian framework for on-linedegradationassessment and residual lifeprediction of secondary batteries in spacecraft[J].Reliability Engineeringamp; System Safety，2013，113：7-20.

[6]CHIEN Y， SHENG T.Mis-Specification Analysis of Linear Degradation Models[J]. IEEE Transactions on Reliability， 2009， 58（3）： 444-55.

[7]ZHANG Z， SI X， HU C，et al. Degradation data analysis and remaining useful life estimation： A review on Wiener-processbased methods[J]. European Journal of Operational Research， 2018，271（3）： 775-96.

[8]PANG Z， SI X， HU C， et al. A Bayesian inference for remaining useful life estimation by fusing accelerated degradation data and condition monitoring data[J]. Reliability Engineeringamp; System Safety，2021， 208.

[9] YE Z， CHEN N， SHEN Y. A new class of Wiener process models for degradation analysis[J].Reliability Engineeringamp; System Safety，2015，139： 58-67.

[10] YAN B，MA X， YANG L， et al. A novel degradation-ratevolatility related effect Wiener process model with its extension to accelerated ageing data analysis [J]. Reliability Engineeringamp; System Safety，2020，204.

[11] ZHOU S，XUA，LIANY，etal. Variational Bayesian analysis forWiener degradation model with random effects[J]. Communications in Statistics- Theory and Methods，2020， 50（16）： 3769-89.

[12] WANG Z， ZHAI Q，CHEN P. Degradation modeling considering unit-to-unit heterogeneity-A general model and comparative study [J].Reliability Engineering amp; System Safety，2021，216.

[13] ZHAI Q， CHEN P， HONG L， et al. A random-effects Wiener degradation model based on accelerated failure time[J]. ReliabilityEngineeringamp; SystemSafety，2018，180： 94-103.

[14]WILLIAMQM，LUISAE. Statistical Methods for Reliability Data[M].New York： John Wileyamp; Sons，1998.

[15] SI X，WANG W，HU C，et al. Remaining Useful Life Estimation Based on a Nonlinear Diffusion Degradation Process[J].IEEE Transactions on Reliability，2012，61（1）： 50- 67.

[16] GEBRAEEL NZ，LAWLEYMA，LIR， et al. Residual-life distributions from component degradation signals：ABayesian approach[J]. IIE Transactions，2005，37（6）： 543-57.

[17]CHRISTOPHE A，NANDO D F，ARNAUD D，et al. An Introductionto MCMC forMachineLearning[J].Machine Learning，2003，50：5-43.

[18] HEIMER A， COHEN I. Multichannel seismic deconvolution usingMarkov-Bernoulli random-field modeling[J].IEEET GeosciRem0te，2009，47（7）：2047-2058.

[19]ROBERTS G O，SMITHAFM.Simple conditions for the convergence of the Gibbs sampler and Metropolis Hastings algorithms[J]. Stoch Proc Appl， 1994， 49（2）： 207-216.

[20]王志達(dá).基于漢密爾頓蒙特卡羅算法的數(shù)控刀架貝葉斯可 靠性評(píng)估[D].長(zhǎng)春：吉林大學(xué)，2021. WANGZD.Bayesian Reliability EvaluationofNC turret based on Hamiton Monte Carlo algorithm[D].Changchun： JilinUniversity，2021.

[21]武健.四參數(shù)Logistic模型的貝葉斯估計(jì)方法比較分析 [D].沈陽(yáng)：沈陽(yáng)師范大學(xué)，2020. WU J. Comparative analysis of Bayesian estimation methods for four-parameter Logistic model[D]. Shenyang： Shenyang Normal University，2020.

[22]付志慧，武健，馬明玥.Rstan包在四參數(shù)Logistic模型參數(shù) 估計(jì)中的應(yīng)用[J].沈陽(yáng)師范大學(xué)學(xué)報(bào)（自然科學(xué)版），2019， 37（4）： 309-314. FUZH，WUJ，MA MY.Application ofRstan package in parameter estimation of four-parameter logistic model[J]. Journal of Shenyang Normal University（Natural Science Edition），2019，37（4）： 309-314.

[23]XIAOM，ZHANG Y，LIY， et al.Degradationmodeling based on Wiener process considering multi-source heterogeneity[J]. IEEEAccess，2020，8：160982-160994.

[24]ABHINAV S，KAI G，DON S，et al. Turbofan engine degradation simulation data set，NASA Ames Prognostics DataRepository[EB/OL]. https：/ti.acr.nasa.gov/tech/dash/ groups/pcoe/prognostic-data-repository.

[25]LI N，LEI Y，YANT，etal.AWiener-Process-Model-Based Method for Remaining Useful Life Prediction Considering Unit-to-Unit Variability[J].IEEE Transactionson Industrial Electronics，2019， 66（3）： 2092-101.

[26]張豹，劉瓊，吳細(xì)寶，等.基于集成學(xué)習(xí)的渦扇發(fā)動(dòng)機(jī)剩余壽 命預(yù)測(cè)[J].中國(guó)測(cè)試，2022，48（7）：47-52. ZHANGB，LIUQ，WUXB，etal.Rmainingusefullifetime prediction of turbofan engine based on ensemble learning[J]. ChinaMeasurementamp; Test，2022，48（7）： 47-52.

[27]任子強(qiáng)，司小勝，胡昌華，等.融合多源數(shù)據(jù)的非線性退化建 模與剩余壽命預(yù)測(cè)[J].中國(guó)測(cè)試，2020，46（2）：1-8. RENZQ，SIXS，HUCH，etal.Multi-sourcedatafusionfor nonlinear degradation modeling and remaining usefullife prediction[J]. China Measurement amp; Test， 2020， 46（2）： 1-8.

[28]MURPHYKP. Conjugate Bayesian analysisof the Gaussian distribution[J].Def，2007，1：1-16.

（編輯：譚玉龍）