簡(jiǎn)介剩余使用壽命預(yù)測(cè)是設(shè)備健康管理與預(yù)測(cè)性維護(hù)領(lǐng)域的核心問題其目標(biāo)是在設(shè)備退化過程中動(dòng)態(tài)評(píng)估其距失效閾值的時(shí)間。隨機(jī)過程模型因能同時(shí)刻畫退化趨勢(shì)與隨機(jī)波動(dòng)并提供壽命的概率分布成為工程中兼顧機(jī)理與數(shù)據(jù)的高性價(jià)比方案。Wiener過程作為典型的隨機(jī)退化建模工具通過漂移項(xiàng)與擴(kuò)散項(xiàng)分別表征平均退化速率與不確定性結(jié)合首達(dá)時(shí)間理論可解析推導(dǎo)出剩余壽命的逆高斯分布。極大似然估計(jì)法可簡(jiǎn)便估計(jì)模型參數(shù)使壽命點(diǎn)預(yù)測(cè)與區(qū)間預(yù)測(cè)均可計(jì)算。該方法適用于軸承、鋰電池、LED光源等具有漸進(jìn)退化特征的部件在工業(yè)現(xiàn)場(chǎng)數(shù)據(jù)驅(qū)動(dòng)與機(jī)理建模之間架起實(shí)用橋梁。本文從Wiener過程的數(shù)學(xué)本質(zhì)出發(fā)推導(dǎo)RUL分布公式并給出完整的Matlab仿真、參數(shù)估計(jì)、滾動(dòng)預(yù)測(cè)及調(diào)參避坑指南幫助工程人員快速落地一套可運(yùn)行的剩余壽命預(yù)測(cè)方案。 搞設(shè)備健康管理這些年我自己最大的感受是真正能落地的壽命預(yù)測(cè)方法往往不是那些結(jié)構(gòu)極其復(fù)雜的深度模型而是數(shù)學(xué)機(jī)理清晰、參數(shù)含義明確、工程上能閉環(huán)的經(jīng)典模型。基于Wiener過程的剩余使用壽命RUL預(yù)測(cè)恰好就是這類方法中的代表。這篇博文我打算從原理、推導(dǎo)、Matlab實(shí)現(xiàn)到調(diào)參避坑完整拆解一套可以直接跑通的項(xiàng)目源碼和數(shù)據(jù)。無論你是在校碩博做PHM課題還是在工業(yè)現(xiàn)場(chǎng)做預(yù)測(cè)性維護(hù)這個(gè)模型都值得認(rèn)真吃透——它既能處理帶趨勢(shì)的隨機(jī)退化過程又能給出壽命的點(diǎn)估計(jì)和區(qū)間估計(jì)而且代碼量控制在幾百行以內(nèi)對(duì)工程部署非常友好。1. 項(xiàng)目整體思路RUL預(yù)測(cè)到底在解決什么問題1.1 從壞了再修到提前預(yù)測(cè)先看這個(gè)項(xiàng)目解決的實(shí)際問題。旋轉(zhuǎn)機(jī)械軸承、鋰電池、LED光源、電容器這類設(shè)備/部件失效過程通常不是瞬間發(fā)生的而是伴隨著某種可測(cè)物理量的漸進(jìn)退化軸承的振動(dòng)幅值逐漸增大電池容量逐漸衰減LED光通量逐漸下降。如果有一個(gè)傳感器能持續(xù)記錄退化量那么我們能回答的核心問題就是還要多久這個(gè)量會(huì)超過失效閾值這個(gè)還要多久就是剩余使用壽命Remaining Useful Life, RUL。它和可靠性壽命等概念有本質(zhì)區(qū)別RUL是動(dòng)態(tài)的隨著觀測(cè)數(shù)據(jù)的積累不斷更新。舉個(gè)例子一塊鋰電池當(dāng)前容量是額定容量的85%退役標(biāo)準(zhǔn)是80%。電壓、溫度、循環(huán)次數(shù)都在影響容量衰減的速率。如果簡(jiǎn)單地用當(dāng)前85%除以初始100%再乘設(shè)計(jì)壽命來估算剩余壽命誤差會(huì)很大——因?yàn)樗p路徑不是線性的而且每塊電池的制造偏差也不小。Wiener過程模型的價(jià)值就在于它把退化路徑建模為確定趨勢(shì) 隨機(jī)波動(dòng)的組合既抓住了平均退化速度也量化了不確定性最后輸出的是一個(gè)完整的剩余壽命概率分布而不僅僅是一個(gè)拍腦袋的數(shù)字。1.2 為什么選用Wiener過程來建模退化很多做RUL預(yù)測(cè)的初學(xué)者第一反應(yīng)是直接用回歸模型擬合退化曲線比如指數(shù)衰減或線性擬合外推到閾值點(diǎn)就算出壽命了。這個(gè)思路不是不行但它有一個(gè)致命的盲區(qū)——它只給出一條平均退化軌跡沒有建模樣本間的隨機(jī)波動(dòng)更沒有給出壽命預(yù)測(cè)的置信區(qū)間。工程上做維修決策時(shí)你除了需要知道估計(jì)還能用100個(gè)小時(shí)更需要知道這個(gè)估計(jì)有多可靠。Wiener過程維納過程的建模思路完全不同。它假設(shè)退化量 (X(t)) 滿足隨機(jī)微分方程[ dX(t) \mu dt \sigma dB(t) ]其中 (B(t)) 是標(biāo)準(zhǔn)布朗運(yùn)動(dòng)(\mu) 是漂移系數(shù)(\sigma) 是擴(kuò)散系數(shù)。這個(gè)方程可以直觀理解為(\mu dt)確定性趨勢(shì)項(xiàng)退化量在單位時(shí)間里平均增加 (\mu)(\sigma dB(t))隨機(jī)波動(dòng)項(xiàng)每次觀測(cè)都會(huì)疊加一個(gè)服從正態(tài)分布的隨機(jī)擾動(dòng)這兩個(gè)參數(shù)恰好對(duì)應(yīng)工程中的平均退化速率和退化過程的不確定性物理含義非常清晰。1.3 各主流RUL預(yù)測(cè)方法橫向?qū)Ρ任以趯?shí)際項(xiàng)目中接觸過不少RUL預(yù)測(cè)方案簡(jiǎn)單做個(gè)對(duì)比方法代表模型優(yōu)勢(shì)劣勢(shì)適用場(chǎng)景基于物理模型Paris裂紋擴(kuò)展、Arrhenius方程機(jī)理清晰、外推能力強(qiáng)需要深入理解失效機(jī)理材料疲勞、單點(diǎn)退化機(jī)理明確基于隨機(jī)過程Wiener過程、Gamma過程、Inverse Gaussian過程同時(shí)給出點(diǎn)估計(jì)與區(qū)間估計(jì)、不確定性量化自然線性漂移假設(shè)約束較強(qiáng)退化軌跡近似單調(diào)或帶波動(dòng)基于數(shù)據(jù)驅(qū)動(dòng)LSTM、GRU、Transformer、CNN無需機(jī)理、擬合能力強(qiáng)訓(xùn)練成本高、需大量失效數(shù)據(jù)、不確定性量化困難傳感器數(shù)據(jù)豐富、非線性強(qiáng)退化Wiener過程在其中的定位很明確它介于物理模型與數(shù)據(jù)驅(qū)動(dòng)之間不需要深挖失效機(jī)制也能建模隨機(jī)波動(dòng)同時(shí)數(shù)學(xué)上可解析推導(dǎo)壽命分布。對(duì)大多數(shù)工程場(chǎng)景來說它是性價(jià)比最高的起點(diǎn)。2. Wiener過程模型的核心原理與RUL推導(dǎo)2.1 Wiener過程的數(shù)學(xué)本質(zhì)跑步機(jī)上的隨機(jī)游走Wiener過程的離散化形式對(duì)寫代碼非常重要。當(dāng)采樣間隔為 (\Delta t) 時(shí)[ X_k X_{k-1} \mu \Delta t \sigma \sqrt{\Delta t} \cdot \varepsilon_k ]其中 (\varepsilon_k \sim N(0,1)) 是標(biāo)準(zhǔn)正態(tài)隨機(jī)數(shù)。這里有一個(gè)非常關(guān)鍵的細(xì)節(jié)也是初學(xué)者最容易寫錯(cuò)的地方隨機(jī)項(xiàng)的方差與步長(zhǎng)成正比。也就是說隨機(jī)擾動(dòng)是 (\sigma \sqrt{\Delta t}) 而不是 (\sigma \Delta t)。為什么因?yàn)椴祭蔬\(yùn)動(dòng)的獨(dú)立增量方差是 (\text{Var}[B(t\Delta t)-B(t)] \Delta t)標(biāo)準(zhǔn)差就是 (\sqrt{\Delta t})。如果代碼里寫成sigma * dt * randn那么步長(zhǎng)細(xì)分時(shí)退化過程的隨機(jī)特性會(huì)收斂到零模擬結(jié)果完全錯(cuò)誤。我用一個(gè)類比來幫助理解把Wiener過程想象成一個(gè)人戴著眼罩在跑步機(jī)上走路。跑步機(jī)的速度決定漂移項(xiàng) (\mu)人的左右晃動(dòng)幅度決定擴(kuò)散項(xiàng) (\sigma)。如果我們每秒鐘記錄一次位置位置的變化量里既有向前的固定推進(jìn)(\mu \Delta t)又有隨機(jī)的左右偏移(\sigma \sqrt{\Delta t} \varepsilon)。注意晃動(dòng)幅度和觀察時(shí)間間隔的平方根成正比觀察時(shí)間越長(zhǎng)累計(jì)的隨機(jī)偏移越大。2.2 首達(dá)時(shí)間與逆高斯分布Wiener過程模型做RUL預(yù)測(cè)的核心是首達(dá)時(shí)間First Hitting Time, FHT概念。設(shè)失效閾值為 (\omega)退化過程 (X(t)) 從初始值 (X(0)x_0) 出發(fā)壽命 (T) 定義為退化量第一次超過閾值的時(shí)刻[ T \inf{t \geq 0: X(t) \geq \omega} ]對(duì)于帶正漂移 (\mu0) 的線性Wiener過程首達(dá)時(shí)間是服從逆高斯分布Inverse Gaussian distribution的隨機(jī)變量其概率密度函數(shù)為[ f_T(t) \frac{\omega - x_0}{\sqrt{2\pi\sigma^2 t^3}} \exp\left(-\frac{(\omega - x_0 - \mu t)^2}{2\sigma^2 t}\right) ]這個(gè)公式長(zhǎng)得有點(diǎn)嚇人但它就是整個(gè)項(xiàng)目的核心公式。用它就能計(jì)算在任一時(shí)刻 (t) 失效發(fā)生的概率密度進(jìn)而得到累計(jì)分布函數(shù)、期望壽命、壽命分位數(shù)等全部關(guān)鍵指標(biāo)。再看RUL。假設(shè)當(dāng)前時(shí)刻為 (t_k)觀測(cè)到的退化量為 (x_k)定義剩余退化量 (a_k \omega - x_k)。那么剩余壽命 (L_k) 的分布可以寫成[ f_{L_k}(l) \frac{a_k}{\sqrt{2\pi\sigma^2 l^3}} \exp\left(-\frac{(a_k - \mu l)^2}{2\sigma^2 l}\right) ]注意這里的一個(gè)關(guān)鍵簡(jiǎn)化剩余壽命只取決于當(dāng)前退化量與閾值之間的距離 (a_k)與過去的歷史路徑無關(guān)。這正是馬爾可夫性的體現(xiàn)。這意味著在預(yù)測(cè)時(shí)我們只需要把當(dāng)前觀測(cè)值 (x_k) 代入公式即可不需要回溯整個(gè)退化歷史。點(diǎn)估計(jì)方面Wiener過程和首達(dá)時(shí)間理論給出了一個(gè)非常優(yōu)雅的結(jié)果當(dāng) (\mu 0) 時(shí)期望剩余壽命為[ E[L_k] \frac{a_k}{\mu} ]也就是說把剩余距離除以平均速度就能得到平均剩余壽命。這和初中物理里的速度-時(shí)間公式一模一樣只是這里的失效時(shí)間是隨機(jī)變量需要加一個(gè)分布去描述它的不確定性。2.3 參數(shù)估計(jì)極大似然估計(jì)與貝葉斯估計(jì)的選擇參數(shù)估計(jì)是Wiener過程RUL預(yù)測(cè)中最容易出問題的一環(huán)。核心要估計(jì)的是漂移系數(shù) (\mu) 和擴(kuò)散系數(shù) (\sigma^2)。假設(shè)在時(shí)間 (t_0, t_1, ..., t_n) 觀測(cè)到退化數(shù)據(jù) (x_0, x_1, ..., x_n)記增量 (\Delta x_k x_k - x_{k-1})時(shí)間間隔 (\Delta t_k t_k - t_{k-1})。由于Wiener過程的增量相互獨(dú)立且服從正態(tài)分布[ \Delta x_k \sim N(\mu \Delta t_k, \sigma^2 \Delta t_k) ]因此對(duì)數(shù)似然函數(shù)為[ \ell(\mu, \sigma^2) -\frac{1}{2} \sum_{k1}^{n} \left[\ln(2\pi\sigma^2\Delta t_k) \frac{(\Delta x_k - \mu\Delta t_k)^2}{\sigma^2\Delta t_k}\right] ]對(duì) (\mu) 和 (\sigma^2) 分別求偏導(dǎo)并令其為零得到極大似然估計(jì)MLE[ \hat{\mu} \frac{x_n - x_0}{t_n - t_0} ][ \hat{\sigma}^2 \frac{1}{n} \sum_{k1}^{n} \frac{(\Delta x_k - \hat{\mu}\Delta t_k)^2}{\Delta t_k} ]如果采樣是等間隔的(\Delta t_k) 恒等于 (\Delta t)那么第二個(gè)公式可以簡(jiǎn)化為[ \hat{\sigma}^2 \frac{1}{n \Delta t} \sum_{k1}^{n} (\Delta x_k - \hat{\mu}\Delta t)^2 ]MLE最大的優(yōu)點(diǎn)是無偏且漸進(jìn)有效在樣本量足夠大的情況下表現(xiàn)很好計(jì)算也極其簡(jiǎn)單。但它的缺點(diǎn)是在小樣本情況下估計(jì)不穩(wěn)定尤其是擴(kuò)散系數(shù)容易偏小。這時(shí)候可以考慮貝葉斯估計(jì)給 (\mu) 和 (\sigma^2) 設(shè)置先驗(yàn)分布然后使用MCMC或者共軛先驗(yàn)解析求解后驗(yàn)分布。工程上我一般建議先用MLE快速驗(yàn)證模型如果發(fā)現(xiàn)結(jié)果對(duì)參數(shù)異常敏感再上貝葉斯。我在自己項(xiàng)目中嘗過一個(gè)小甜頭先對(duì)原始退化數(shù)據(jù)進(jìn)行平滑比如5點(diǎn)滑動(dòng)平均再估計(jì)參數(shù)。平滑能顯著降低 (\hat{\sigma}^2) 的波動(dòng)讓RUL區(qū)間估計(jì)更穩(wěn)定。代價(jià)是對(duì)(\sigma)的估計(jì)會(huì)偏小不過在工程上這個(gè)偏小通常在可接受范圍內(nèi)。2.4 線性模型的邊界什么時(shí)候必須升級(jí)線性Wiener過程模型有一個(gè)內(nèi)在假設(shè)——漂移系數(shù) (\mu) 在退化全程保持不變。這個(gè)假設(shè)對(duì)某些退化過程成立但對(duì)很多工程部件并不成立。舉個(gè)典型例子鋰電池衰減。電池在早期循環(huán)中容量衰減較慢后期由于內(nèi)阻增大、活性物質(zhì)損失加速容量衰減明顯加快。這時(shí)候如果強(qiáng)行用線性Wiener過程擬合全部數(shù)據(jù)(\mu) 估計(jì)出來會(huì)是一個(gè)平均速率導(dǎo)致后期預(yù)測(cè)嚴(yán)重滯后壽命高估。面對(duì)這種情況有幾個(gè)升級(jí)路徑非線性漂移函數(shù)把 (\mu t) 改成 (\mu \cdot \Lambda(t))其中 (\Lambda(t)) 是時(shí)間尺度變換函數(shù)比如冪函數(shù) (t^\beta)。這樣能捕捉加速退化。隨機(jī)效應(yīng)模型假設(shè)漂移系數(shù) (\mu) 在不同個(gè)體之間服從正態(tài)分布 (\mu \sim N(\mu_0, \sigma_\mu^2))用于刻畫設(shè)備間的制造差異。退化-沖擊復(fù)合模型在Wiener過程基礎(chǔ)上疊加隨機(jī)沖擊的影響適合同時(shí)存在漸進(jìn)退化和突發(fā)沖擊的場(chǎng)景。這些擴(kuò)展都會(huì)顯著增加數(shù)學(xué)復(fù)雜度但核心思想仍然圍繞趨勢(shì) 波動(dòng)展開理解了基礎(chǔ)版后續(xù)升級(jí)就順理成章。3. Matlab完整實(shí)現(xiàn)從數(shù)據(jù)生成到RUL預(yù)測(cè)3.1 項(xiàng)目文件結(jié)構(gòu)與主程序框架我在設(shè)計(jì)Matlab項(xiàng)目時(shí)傾向于把功能拆分成清晰的小文件便于調(diào)試和復(fù)用。這個(gè)項(xiàng)目的目錄結(jié)構(gòu)如下RUL_Wiener/ ├─ main_rul_prediction.m % 主程序數(shù)據(jù)生成、參數(shù)估計(jì)、RUL預(yù)測(cè)、繪圖 ├─ generate_degradation.m % 生成Wiener退化仿真數(shù)據(jù) ├─ estimate_params.m % MLE參數(shù)估計(jì) ├─ predict_rul.m % RUL概率密度計(jì)算與點(diǎn)估計(jì)/區(qū)間估計(jì) └─ plot_results.m % 可視化輔助函數(shù)雖然你也可以把所有代碼塞進(jìn)一個(gè)腳本里但工程上還是建議拆開。原因很簡(jiǎn)單真實(shí)項(xiàng)目里數(shù)據(jù)通常來自傳感器文件而不是仿真生成函數(shù)generate_degradation會(huì)被替換成load_sensor_data保留清晰的函數(shù)接口會(huì)讓替換成本降到最低。3.2 退化數(shù)據(jù)仿真怎么生成像樣的數(shù)據(jù)先來解決數(shù)據(jù)從哪來的問題。如果你手頭沒有公開數(shù)據(jù)集或者實(shí)驗(yàn)數(shù)據(jù)仿真生成是驗(yàn)證算法正確性的最佳途徑——答案已知跑一遍模型就知道代碼寫對(duì)沒有。這里演示生成一段線性退化趨勢(shì)疊加隨機(jī)波動(dòng)的仿真數(shù)據(jù)function [t, X] generate_degradation(mu, sigma, dt, N, x0) % 生成Wiener過程退化仿真數(shù)據(jù) % 輸入: % mu - 漂移系數(shù)單位/單位時(shí)間 % sigma - 擴(kuò)散系數(shù)單位/sqrt(單位時(shí)間) % dt - 采樣時(shí)間間隔 % N - 采樣點(diǎn)數(shù)量 % x0 - 初始退化量 % 輸出: % t - 時(shí)間向量 (1xN) % X - 退化數(shù)據(jù)向量 (1xN) t (0:N-1) * dt; X zeros(1, N); X(1) x0; for k 2:N X(k) X(k-1) mu * dt sigma * sqrt(dt) * randn(); end end關(guān)鍵點(diǎn)有兩個(gè)一是sigma * sqrt(dt) * randn()這里不能漏掉sqrt。前面說過布朗運(yùn)動(dòng)增量的方差與時(shí)間間隔成正比只要這里的縮放關(guān)系正確無論dt取多大或多小模擬出的退化軌跡在統(tǒng)計(jì)意義上都是一致的。二是固定隨機(jī)數(shù)種子。仿真階段建議在主程序中加上rng(42)否則每次運(yùn)行結(jié)果都不同調(diào)試參數(shù)時(shí)你根本分不清是代碼改對(duì)了還是運(yùn)氣好。我習(xí)慣把隨機(jī)種子作為主程序的一個(gè)可配置參數(shù)方便復(fù)現(xiàn)實(shí)驗(yàn)結(jié)果。主程序里調(diào)用它rng(42); mu_true 0.05; % 每單位時(shí)間退化0.05 sigma_true 0.03; % 隨機(jī)波動(dòng)強(qiáng)度 dt 0.02; % 采樣間隔 N 500; % 采樣點(diǎn)數(shù) x0 0; [t, X] generate_degradation(mu_true, sigma_true, dt, N, x0);這樣一條帶隨機(jī)波動(dòng)的退化曲線就生成好了。如果想模擬加速退化可以取指數(shù)漂移mu*lambda(t)但基礎(chǔ)版先保持線性。3.3 滑動(dòng)窗口分割與訓(xùn)練/測(cè)試策略RUL預(yù)測(cè)和普通回歸預(yù)測(cè)有一個(gè)很大的區(qū)別預(yù)測(cè)時(shí)刻是流動(dòng)的。你不可能在設(shè)備出廠時(shí)預(yù)測(cè)一次就完事而是要在t_k時(shí)刻基于當(dāng)前所有觀測(cè)做一次預(yù)測(cè)隨后每個(gè)采樣點(diǎn)都重新預(yù)測(cè)一次形成一個(gè)滾動(dòng)更新的預(yù)測(cè)序列。實(shí)操中我建議以當(dāng)前時(shí)間 (t_k) 為界將全部觀測(cè)分成兩部分歷史窗口((t_1, x_1), ..., (t_k, x_k))用來估計(jì)模型參數(shù)預(yù)測(cè)目標(biāo)當(dāng)前時(shí)刻剩余的壽命分布主程序里用一個(gè)循環(huán)遍歷預(yù)測(cè)起始點(diǎn)。假設(shè)前 (N_0200) 個(gè)點(diǎn)用于初始訓(xùn)練之后每隔step10個(gè)點(diǎn)做一次RUL預(yù)測(cè)直到數(shù)據(jù)末尾t_predict 200:10:N-1; RUL_mean zeros(1, length(t_predict)); RUL_std zeros(1, length(t_predict)); RUL_lower zeros(1, length(t_predict)); RUL_upper zeros(1, length(t_predict)); for idx 1:length(t_predict) k t_predict(idx); % 當(dāng)前預(yù)測(cè)時(shí)刻的索引 mu_hat (X(k) - X(1)) / (t(k) - t(1)); dX diff(X(1:k)); sigma2_hat mean((dX - mu_hat*dt).^2) / dt; [RUL_mean(idx), RUL_std(idx), RUL_lower(idx), RUL_upper(idx)] ... predict_rul(X(k), omega, mu_hat, sigma2_hat); end這種滾動(dòng)預(yù)測(cè)方式模擬了真實(shí)設(shè)備監(jiān)控的場(chǎng)景隨著設(shè)備運(yùn)行時(shí)間增長(zhǎng)預(yù)測(cè)結(jié)果應(yīng)該越來越準(zhǔn)確不確定區(qū)間應(yīng)該越來越窄。如果畫出RUL預(yù)測(cè)值與真實(shí)剩余壽命隨時(shí)間變化的曲線模型效果一目了然。3.4 參數(shù)估計(jì)與RUL計(jì)算的Matlab實(shí)現(xiàn)參數(shù)估計(jì)函數(shù)實(shí)現(xiàn)非常簡(jiǎn)潔本質(zhì)上就是套用前面推導(dǎo)的MLE公式function [mu_hat, sigma2_hat] estimate_params(t, X) % 基于MLE估計(jì)線性Wiener過程的漂移與擴(kuò)散參數(shù) n length(X); mu_hat (X(n) - X(1)) / (t(n) - t(1)); % 總變化量 / 總時(shí)間 dX diff(X); dt diff(t); sigma2_hat mean((dX - mu_hat .* dt).^2 ./ dt); end這里有一個(gè)容易被忽略的細(xì)節(jié)當(dāng)采樣不均勻即dt不是常數(shù)時(shí)必須用帶權(quán)重的公式也就是每一項(xiàng)除以其對(duì)應(yīng)的時(shí)間間隔。如果忽略這個(gè)細(xì)節(jié)直接用均一的dt會(huì)在采樣稀疏的區(qū)段引入偏差。等間隔采樣是沒問題的但真實(shí)傳感器數(shù)據(jù)經(jīng)常有掉幀、延遲代碼里處理非均勻時(shí)間間隔的寫法屬于防御性編程。預(yù)測(cè)函數(shù)是整段代碼的靈魂它實(shí)現(xiàn)了逆高斯分布的RUL概率密度計(jì)算以及點(diǎn)估計(jì)和區(qū)間估計(jì)function [mean_rul, std_rul, lower, upper] predict_rul(x_current, omega, mu_hat, sigma2_hat, alpha) % 計(jì)算剩余壽命的點(diǎn)估計(jì)和區(qū)間估計(jì) % alpha 為置信水平默認(rèn)0.9 if nargin 5 alpha 0.9; end a omega - x_current; % 剩余退化量 if a 0 error(當(dāng)前退化量已超過閾值RUL為0); end % 點(diǎn)估計(jì)期望RUL 剩余退化量 / 平均漂移率 mean_rul a / mu_hat; % 逆高斯分布方差 var_rul a * sigma2_hat^2 / mu_hat^3; std_rul sqrt(var_rul); % 在時(shí)間網(wǎng)格上計(jì)算PDF和CDF用數(shù)值方法求分位數(shù) L_max max(3 * mean_rul, mean_rul 10 * std_rul); l_grid linspace(1e-6, L_max, 5000); % 逆高斯分布PDF pdf_l a ./ sqrt(2 * pi * sigma2_hat^2 * l_grid.^3) .* ... exp(-(a - mu_hat * l_grid).^2 ./ (2 * sigma2_hat^2 * l_grid)); % 數(shù)值積分求CDF cdf_l cumtrapz(l_grid, pdf_l); cdf_l min(max(cdf_l, 0), 1); % 防止數(shù)值誤差導(dǎo)致的越界 % 分位數(shù) low_q (1 - alpha) / 2; high_q 1 - (1 - alpha) / 2; lower interp1(cdf_l, l_grid, low_q); upper interp1(cdf_l, l_grid, high_q); end這里注意幾個(gè)工程細(xì)節(jié)方差公式是逆高斯分布的解析式(\text{Var}(L) a\sigma^2 / \mu^3)。但用它來算標(biāo)準(zhǔn)差有一個(gè)問題——逆高斯分布是右偏的直接用均值加減標(biāo)準(zhǔn)差構(gòu)造的區(qū)間并非嚴(yán)格置信區(qū)間所以代碼里用數(shù)值積分和分位數(shù)法這才是嚴(yán)格的置信區(qū)間。L_max的選擇要自適應(yīng)。如果固定一個(gè)很大的上限概率密度在網(wǎng)格上的分辨率不足如果太保守又會(huì)截?cái)辔膊繉?dǎo)致分位數(shù)計(jì)算錯(cuò)誤。按mean 10*std取上限在絕大多數(shù)情況下是安全的。分位數(shù)查找用interp1比用find更精確而且不需要預(yù)先知道CDF的單調(diào)區(qū)間端點(diǎn)恰好落在網(wǎng)格點(diǎn)上。3.5 結(jié)果可視化與評(píng)價(jià)指標(biāo)做預(yù)測(cè)不能只輸出一行數(shù)字需要畫圖來判斷模型行為是否合理。我通常畫出三張圖第一張是退化軌跡圖標(biāo)注失效閾值和當(dāng)前預(yù)測(cè)時(shí)刻。figure; plot(t, X, b-, LineWidth, 1.2); hold on; yline(omega, r--, Failure Threshold, LineWidth, 1.5); xlabel(Time); ylabel(Degradation Value); title(Degradation Trajectory and Failure Threshold); legend(Degradation, Threshold); grid on;第二張是RUL概率密度曲線。選取幾個(gè)代表性預(yù)測(cè)時(shí)刻比如早期、中期、臨近失效各自畫一條逆高斯分布的PDF曲線。隨著預(yù)測(cè)時(shí)刻越來越接近真實(shí)失效時(shí)間曲線應(yīng)該越來越窄、峰值越來越靠近真實(shí)RUL。第三張是滾動(dòng)預(yù)測(cè)結(jié)果圖。橫軸是真實(shí)剩余壽命縱軸是預(yù)測(cè)RUL理想情況下預(yù)測(cè)值應(yīng)該沿著45度對(duì)角線分布。把預(yù)測(cè)均值畫成點(diǎn)把置信區(qū)間畫成誤差棒就能直觀看到模型是否過于樂觀或過于悲觀。評(píng)估指標(biāo)方面我最常用的三個(gè)指標(biāo)公式說明MAE(\frac{1}{N}\sum | \hat{L}_i - L_i |)平均絕對(duì)誤差越小越好RMSE(\sqrt{\frac{1}{N}\sum (\hat{L}_i - L_i)^2})均方根誤差對(duì)大的離群誤差敏感區(qū)間覆蓋概率(\frac{1}{N}\sum \mathbb{I}(L_i \in [low_i, up_i]))真實(shí)RUL落在預(yù)測(cè)區(qū)間內(nèi)的比例接近置信水平說明區(qū)間標(biāo)定正確區(qū)間覆蓋概率這個(gè)指標(biāo)特別重要很多人會(huì)忽略。如果你的90%置信區(qū)間實(shí)際只有60%的覆蓋率說明模型對(duì)不確定性估計(jì)過于保守會(huì)誤導(dǎo)維修決策。我跑仿真數(shù)據(jù)時(shí)積累了一個(gè)經(jīng)驗(yàn)如果區(qū)間覆蓋率長(zhǎng)期低于理論值首要懷疑的是擴(kuò)散系數(shù) (\sigma^2) 估計(jì)偏小常見原因包括數(shù)據(jù)采樣率過低、傳感器量化誤差導(dǎo)致信息丟失或者退化過程本身是非線性的線性模型壓縮了波動(dòng)范圍。4. 常見問題與排查實(shí)錄4.1 數(shù)據(jù)預(yù)處理階段容易踩的坑Wiener過程模型假設(shè)退化數(shù)據(jù)滿足獨(dú)立增量正態(tài)分布但真實(shí)傳感器數(shù)據(jù)幾乎不可能直接滿足。我的經(jīng)驗(yàn)是正式建模前必須先做兩件事平滑去噪和異常值剔除。平滑推薦用滑動(dòng)平均或Savitzky-Golay濾波。前者適合噪聲均勻的情況后者在保留退化趨勢(shì)拐點(diǎn)方面效果更好。Matlab里smoothdata和sgolayfilt都可以直接用。注意平滑窗口不要選太大窗口過大會(huì)把真實(shí)的退化趨勢(shì)細(xì)節(jié)抹平導(dǎo)致(\hat{\sigma}^2)估計(jì)嚴(yán)重偏小。我試過用100點(diǎn)窗口平滑500點(diǎn)數(shù)據(jù)結(jié)果RUL區(qū)間窄到不真實(shí)后來改成15點(diǎn)窗口才恢復(fù)正常。異常值的處理更關(guān)鍵。Wiener過程對(duì)異常值非常敏感——一個(gè)漂移過大的離群點(diǎn)會(huì)直接拉高(\hat{\sigma}^2)導(dǎo)致置信區(qū)間異常膨脹。在預(yù)處理階段用3σ原則或者Hampel濾波器剔除異常值能顯著提升估計(jì)穩(wěn)定性。注意如果在預(yù)處理階段使用了平滑操作事后評(píng)估區(qū)間覆蓋率時(shí)要意識(shí)到(\hat{\sigma}^2)已經(jīng)偏小。別等到結(jié)果出來對(duì)不上再回頭找原因。4.2 參數(shù)初值與數(shù)值穩(wěn)定性問題MLE的解是解析的理論上不存在初值問題。但如果你擴(kuò)展模型比如非線性漂移、隨機(jī)效應(yīng)需要數(shù)值優(yōu)化這時(shí)候初值選擇就很重要。我的建議是先用線性MLE結(jié)果作為非線性優(yōu)化的初值。比如冪函數(shù)漂移 (\Lambda(t)t^\beta)先用線性模型估計(jì) (\mu)再固定 (\mu) 用一維搜素估計(jì) (\beta)最后所有參數(shù)聯(lián)合優(yōu)化。這樣能避開局部最優(yōu)。數(shù)值穩(wěn)定性方面最容易出問題的地方是RUL概率密度公式中的指數(shù)項(xiàng)。當(dāng)剩余退化量 (a) 很小而時(shí)間 (l) 也較小時(shí)((a - \mu l)^2 / (2\sigma^2 l)) 可能出現(xiàn)較大的中間值導(dǎo)致指數(shù)下溢。我處理的做法是先用logpdf計(jì)算對(duì)數(shù)密度再取指數(shù)必要的時(shí)候用logsumexp技巧歸一化。Matlab里直接對(duì)概率密度做cumtrapz的前提是PDF數(shù)值上沒有出現(xiàn)NaN或Inf所以在仿真數(shù)據(jù)里我把時(shí)間網(wǎng)格的起點(diǎn)設(shè)成1e-6避免除以零。4.3 預(yù)測(cè)結(jié)果偏保守或偏激進(jìn)的調(diào)參方向如果你跑出來的RUL預(yù)測(cè)系統(tǒng)性地偏離真實(shí)值別急著改代碼——先問自己?jiǎn)栴}出在參數(shù)估計(jì)還是模型假設(shè)本身先說參數(shù)層面。如果預(yù)測(cè)的RUL均值比真實(shí)值偏小即預(yù)測(cè)過早失效通常意味著 (\hat{\mu}) 被高估。常見原因是后期退化加速而MLE用了全時(shí)段平均漂移率。把參數(shù)估計(jì)窗口改成滑動(dòng)窗口只使用最近一段時(shí)間的退化數(shù)據(jù)來估計(jì) (\mu)往往能改善。反過來如果預(yù)測(cè)RUL系統(tǒng)偏大預(yù)測(cè)過晚失效常見原因是數(shù)據(jù)早期退化緩慢拉低了平均漂移率?；蛘唛撝翟O(shè)置不合理實(shí)際失效閾值比設(shè)定的 (\omega) 更低。這時(shí)候要重新審視閾值是怎么來的——如果閾值是拍腦袋定的預(yù)測(cè)結(jié)果很難準(zhǔn)確。還有一類情況是RUL預(yù)測(cè)隨觀測(cè)時(shí)間劇烈波動(dòng)一會(huì)兒偏大一會(huì)兒偏小。這在數(shù)據(jù)噪聲比較強(qiáng)時(shí)很常見。解決辦法是增加觀測(cè)數(shù)據(jù)的平滑、延長(zhǎng)參數(shù)估計(jì)用的時(shí)間窗口或者對(duì)多個(gè)歷史時(shí)刻的預(yù)測(cè)結(jié)果做時(shí)間上的平滑濾波。4.4 模型失效的三大典型特征如果在項(xiàng)目中使用Wiener過程時(shí)遇到以下三種現(xiàn)象基本可以判定模型與數(shù)據(jù)不匹配需要升級(jí)模型架構(gòu)第一種退化增量不是正態(tài)分布。用histogram畫一下增量直方圖如果呈現(xiàn)明顯偏態(tài)或厚尾Wiener過程就不好用了。此時(shí)可以考慮Gamma過程專門建模單調(diào)遞增的退化過程或逆高斯過程。第二種退化路徑存在明顯非線性。把所有樣本歸一化到同一時(shí)間尺度畫出退化軌跡。如果軌跡呈明顯曲率比如下凸的加速退化線性Wiener過程會(huì)產(chǎn)生系統(tǒng)性的壽命預(yù)測(cè)偏差。升級(jí)為指數(shù)漂移Wiener過程或者利用時(shí)間尺度變換效果立竿見影。第三種個(gè)體差異過大。如果是從多個(gè)同類設(shè)備收集的數(shù)據(jù)且不同設(shè)備間平均退化速率差異巨大單一Wiener過程無法刻畫這種heterogeneity。用隨機(jī)效應(yīng)模型讓漂移系數(shù) (\mu) 服從一個(gè)超先驗(yàn)分布能顯著提升預(yù)測(cè)精度。這三種情況在學(xué)術(shù)論文和工業(yè)項(xiàng)目里都很常見。所謂模型的邊界就在這里——不是Wiener過程不好用而是你需要在正確的地圖里使用它。最后再分享一個(gè)小技巧做參數(shù)估計(jì)和RUL預(yù)測(cè)時(shí)建議把隨機(jī)種子固定下來同時(shí)把完整參數(shù)配置寫成一個(gè)結(jié)構(gòu)體存下來。哪怕是仿真數(shù)據(jù)每次運(yùn)行結(jié)果也應(yīng)該嚴(yán)格一致。這樣當(dāng)后續(xù)結(jié)果出現(xiàn)異常時(shí)你能確定是代碼邏輯變化導(dǎo)致的而不是隨機(jī)性造成的。我在做這個(gè)項(xiàng)目的過程中深刻體會(huì)到RUL預(yù)測(cè)的本質(zhì)是在不確定性中給出判斷——分布比單點(diǎn)更重要區(qū)間覆蓋率比平均誤差更需要盯住。希望這份實(shí)現(xiàn)和踩坑記錄能幫你少繞幾段彎路。本文還有配套的精品資源點(diǎn)擊獲取