基于Rényi熵的q-指數(shù)分布及其可靠性分析應(yīng)用

謝 暄，王敏夷，白穎利，汪東敏，李西峰*，謝永樂(lè)

(1. 電子科技大學(xué)自動(dòng)化工程學(xué)院 成都 611731；2. 中國(guó)空間技術(shù)研究院通信衛(wèi)星事業(yè)部 北京 海淀區(qū) 100094；3. 四川慧龍科技有限責(zé)任公司 成都 610041)

隨著集成系統(tǒng)復(fù)雜性的日益增加，系統(tǒng)建模和可靠性保證面臨越來(lái)越大的挑戰(zhàn)。在可靠性分析的各種場(chǎng)景中，雖然指數(shù)分布模型已被較廣泛地用于建立系統(tǒng)壽命模型[1]，但由于指數(shù)分布的風(fēng)險(xiǎn)函數(shù)是恒定的，因此直接將指數(shù)分布用于描述系統(tǒng)壽命，存在無(wú)法描述損傷過(guò)程和無(wú)法準(zhǔn)確反映故障累積效果的問(wèn)題，最典型的例子是將指數(shù)分布簡(jiǎn)單應(yīng)用于描述人類死亡率和電子設(shè)備生命周期，效果不夠理想。對(duì)這類過(guò)程，通常需要采用具備浴盆特征的風(fēng)險(xiǎn)函數(shù)所對(duì)應(yīng)的壽命分布來(lái)準(zhǔn)確描述。韋伯分布作為指數(shù)分布的概括，以及它的帶有浴盆型風(fēng)險(xiǎn)函數(shù)的擴(kuò)展得到了重視，并被廣泛應(yīng)用于許多領(lǐng)域[2-4]。

目前關(guān)于指數(shù)分布推廣的研究，主要集中在廣義韋伯分布的累積分布函數(shù)的參數(shù)加法上。大多數(shù)情況下，這些韋伯類型的一般化是通過(guò)簡(jiǎn)單的參數(shù)加法技術(shù)來(lái)實(shí)現(xiàn)的。此方法一般除保留了原有的參數(shù)，還引入了一些新的參數(shù)，使得模型對(duì)壽命數(shù)據(jù)的擬合，多數(shù)情況下優(yōu)于沒有新參數(shù)的模型。

基于參數(shù)加法，文獻(xiàn)[5]提出了四參數(shù)廣義指數(shù)模型，稱為廣義指數(shù)化線性指數(shù)分布。結(jié)果表明，這一擴(kuò)展可以導(dǎo)出一系列指數(shù)型分布，如指數(shù)化韋伯分布。結(jié)合韋伯分布和改進(jìn)的五參數(shù)韋伯分布，文獻(xiàn)[6]提出了一種修正的韋伯分布擴(kuò)展?；诖?，文獻(xiàn)[7]又提出了一個(gè)離散五參數(shù)修正型的韋伯分布，并發(fā)現(xiàn)基于這種分布的離散數(shù)據(jù)擬合勝過(guò)其他3個(gè)修正的韋伯模型。

盡管韋伯推廣在壽命模型中具有良好的效果，但由于參數(shù)估計(jì)過(guò)程復(fù)雜，參數(shù)物理意義不明確(至少這些參數(shù)在這些廣義的分布中，不能直接指示系統(tǒng)的當(dāng)前狀態(tài))，限制了韋伯推廣的應(yīng)用[8-10]。

本文提出了一種廣義雙參數(shù)指數(shù)分布，它可以作為可靠性分析的基礎(chǔ)，例如浴盆型風(fēng)險(xiǎn)函數(shù)的構(gòu)造和有用壽命預(yù)測(cè)。與傳統(tǒng)的參數(shù)加法不同，本文用最大熵原理得到廣義指數(shù)分布。從物理角度看，這種廣義分布的主要優(yōu)點(diǎn)是兩個(gè)參數(shù)都具有明顯的物理意義。一個(gè)參數(shù)具有分形意義，表示系統(tǒng)的穩(wěn)定程度；另一個(gè)參數(shù)則與系統(tǒng)的平均行為密切相關(guān)。兩個(gè)參數(shù)均可以用來(lái)刻畫系統(tǒng)的性能。

1 q-指數(shù)分布模型

1.1 一般化的指數(shù)分布

基于Rényi熵，可以通過(guò)最大熵方法，得到廣義指數(shù)分布：

式中，q是分形參數(shù)；λ是分布的期望值；q和λ都是非負(fù)的。根據(jù)參數(shù)q和λ的值，概率密度函數(shù)fq會(huì)有不同的形態(tài)。圖1所示為不同的參數(shù)q時(shí)固定均值λ下的fq形態(tài)?？梢钥闯鰍-指數(shù)分布是單峰分布，也是右偏分布。因此，q-指數(shù)分布可以被考慮用于對(duì)壽命相關(guān)的可靠性分析。

圖1 q-指數(shù)分布fq(x)

計(jì)算可得，q-指數(shù)分布的累積分布函數(shù)是：

式中，2F1是高斯超幾何函數(shù)。

1.2 特殊衍生分布

下面將給出q-指數(shù)分布的兩個(gè)主要的衍生分布，它們?cè)谏娣治?、壓縮感知、剩余使用壽命預(yù)測(cè)和其他的可靠性相關(guān)領(lǐng)域有潛在實(shí)用價(jià)值。

首先，當(dāng)q→1時(shí)，該分布簡(jiǎn)化為均值為λ的指數(shù)分布：

其次，使用式(1)可以得到一個(gè)廣義的q-拉普拉斯分布：

圖2顯示了不同參數(shù)值q和固定λ值(λ=10)的廣義拉普拉斯分布。拉普拉斯分布在壓縮感知領(lǐng)域起到了重要作用[11-13]，這里衍生出的拉普拉斯分布在統(tǒng)計(jì)信號(hào)處理和機(jī)器學(xué)習(xí)領(lǐng)域具有潛在應(yīng)用價(jià)值。

圖2 不同參數(shù)q下的q-拉普拉斯分布gq(x)

2 模型性質(zhì)

2.1 q-指數(shù)分布的風(fēng)險(xiǎn)函數(shù)

生存函數(shù)是表示一系列事件的隨機(jī)變量函數(shù)，通常用于表示一些基于時(shí)間的系統(tǒng)失敗或死亡概率。假設(shè)T表示產(chǎn)品使用壽命，其分布函數(shù)為F(t)，那么該產(chǎn)品壽命大于t的概率為：

S(t)=P(T＞t)=1-F(t)

式中，S(t)被稱為生存函數(shù)。在此基礎(chǔ)上，可用風(fēng)險(xiǎn)函數(shù)刻畫已有效使用到t時(shí)刻的產(chǎn)品，在[t,t+Δt]極短時(shí)間內(nèi)“死亡”的風(fēng)險(xiǎn)：

據(jù)此計(jì)算可得q-指數(shù)分布的生存函數(shù)：

風(fēng)險(xiǎn)函數(shù)：

圖3為具有不同參數(shù)值的風(fēng)險(xiǎn)函數(shù)圖像，可以看到風(fēng)險(xiǎn)函數(shù)具有許多不同的形狀，q-指數(shù)分布中當(dāng)q＞1，顯示遞增特性；當(dāng)q＜1，顯示遞減性。

圖3 q-指數(shù)分布的風(fēng)險(xiǎn)函數(shù)

注意到q-指數(shù)分布的風(fēng)險(xiǎn)函數(shù)呈現(xiàn)多態(tài)性，其中部分具有澡盆特征。可根據(jù)實(shí)際使用環(huán)境，選擇恰當(dāng)?shù)膓-指數(shù)建模風(fēng)險(xiǎn)過(guò)程，從而可以提高壽命估計(jì)精度。

2.2 q-指數(shù)分布的矩

q-指數(shù)分布fq的均值是λ，方差是：

可以得到分布的k階矩為：

式中，k表示分布的k階矩。

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

為了在不同場(chǎng)景下選擇恰當(dāng)?shù)膮?shù)，這里給出兩種估計(jì)q-指數(shù)分布中參數(shù)的方法：最大似然估計(jì)法(maximum likelihood estimation,MLE)和信息似然估計(jì)法。最大似然估計(jì)法適用于滿足高斯分布的數(shù)據(jù)集，信息似然估計(jì)法適用于高斯和非高斯分布的數(shù)據(jù)集。

3.1 最大似然估計(jì)

根據(jù)式(1)，建立對(duì)數(shù)似然函數(shù)：

然后有：

q和λ可以通過(guò)求解方程組(6)獲得：

3.2 信息似然估計(jì)

另一個(gè)估計(jì)Rényi信息未知參數(shù)的方法，將信息論和譜估計(jì)相結(jié)合，可在最小先驗(yàn)的條件下，取得經(jīng)驗(yàn)風(fēng)險(xiǎn)最小的參數(shù)估計(jì)值[14]。具體步驟如下：

首先，已知Rényi信息頻譜定義為：

當(dāng)λ →1，

Rényi信息頻譜趨向于香農(nóng)熵。

同時(shí)，Rényi信息的頻譜梯度被定義為：

式中，D(·)是隨機(jī)變量x的方差。

令φf(shuō):=-2LR(1)

一方面，將式(1)中的fq代 入式(9)計(jì)算φf(shuō)，根據(jù)文獻(xiàn)[13]可得:

另一方面，使用核方法估計(jì) φf(shuō)。具體而言，φf(shuō)可以通過(guò)核方法計(jì)算如下：

式中，

fn(x)是為無(wú)參數(shù)內(nèi)核密度估計(jì)量，定義為：

假設(shè)K(·)為有界變分的概率密度函數(shù)(核函數(shù))，其支撐集位于部分有限區(qū)間。Xi是隨機(jī)變量X的樣本。設(shè)bn為滿足以下收斂條件的序列：

比較式(10)和式(11)，式(1)中的參數(shù)q可以通過(guò)求解式(13)得到：

4 數(shù)據(jù)分析

本文使用了3個(gè)與可靠性相關(guān)的實(shí)驗(yàn)去評(píng)估q-指數(shù)分布在可靠性分析中的有效性，包括白血病病人生存期分析、設(shè)備壽命預(yù)測(cè)以及鋰電池壽命預(yù)測(cè)。

第一個(gè)數(shù)據(jù)集是40個(gè)白血病病人的生存期數(shù)據(jù)[15]，它具有不斷增加的風(fēng)險(xiǎn)率。第二個(gè)數(shù)據(jù)集來(lái)自50個(gè)組件的樣本[16]，具有浴盆型風(fēng)險(xiǎn)的風(fēng)險(xiǎn)率。利用K-S統(tǒng)計(jì)量評(píng)估了q-指數(shù)與上述兩個(gè)著名數(shù)據(jù)集的擬合結(jié)果，并計(jì)算出赤池信息準(zhǔn)則(AIC)值和貝葉斯信息準(zhǔn)則(BIC)值，用于比較它和其他廣義指數(shù)分布的擬合優(yōu)越度。第三個(gè)實(shí)驗(yàn)建立了關(guān)于3組鋰電池的容量退化模型，與普通指數(shù)分布的結(jié)果相比，利用q-指數(shù)分布可以得到更準(zhǔn)確的剩余壽命預(yù)測(cè)結(jié)果。

4.1 白血病病人生存期

該數(shù)據(jù)集由文獻(xiàn)[15]給出，數(shù)據(jù)源自沙特阿拉伯衛(wèi)生醫(yī)院部門，如表1所示。它記錄了40名白血病患者的生命周期。

表1 40名白血病患者的生存期

基于雙參數(shù)模型的生存函數(shù)經(jīng)驗(yàn)估計(jì)如圖4所示，圖5描述了模型的經(jīng)驗(yàn)和擬合風(fēng)險(xiǎn)函數(shù)。同時(shí)，由于該數(shù)據(jù)集所示患者死亡風(fēng)險(xiǎn)是隨時(shí)間上升的，因此，韋伯分布、線性失效率分布(linear failure rate distribution, LFR)和q-指數(shù)分布均可作為數(shù)據(jù)擬合的候選者。為了判斷上述哪個(gè)密度函數(shù)更適合數(shù)據(jù)擬合，用MLE方法對(duì)3種模型參數(shù)進(jìn)行了估計(jì)，如表2所示。然后利用幾種不同的測(cè)試統(tǒng)計(jì)測(cè)度對(duì)擬合結(jié)果進(jìn)行了評(píng)估。

表2 白血病數(shù)據(jù)集各模型參數(shù)的MLE值

圖4 白血病數(shù)據(jù)集的生存函數(shù)

由統(tǒng)計(jì)學(xué)可知，K-S值表示分隔程度，一般大于0.2即表明模型具有良好的分隔性能。由表3可知，q-指數(shù)分布與韋伯分布、指數(shù)分布和LFR分布的K-S均有良好分隔能力。同時(shí)，由表3列出的對(duì)數(shù)似然函數(shù)值可知，q-指數(shù)分布與韋伯分布、指數(shù)分布的最大似然函數(shù)值均在相近水平，所以q-指數(shù)分布可以充分利用先驗(yàn)信息獲得對(duì)未知參數(shù)的最大似然估計(jì)。

表3 白血病數(shù)據(jù)的對(duì)數(shù)似然函數(shù)值和K-S

根據(jù)統(tǒng)計(jì)學(xué)原理，p值是判斷原假設(shè)是否成立的依據(jù)，一般認(rèn)為p＞0.05，說(shuō)明兩組樣本無(wú)統(tǒng)計(jì)學(xué)差異。由表4可知，q-指數(shù)分布與韋伯分布的p值在相近水平，q-指數(shù)分布具備描述壽命分布的能力。同時(shí)，通過(guò)計(jì)算AIC或者BIC值，相比指數(shù)分布、韋伯分布和LFR分布，不難發(fā)現(xiàn)q-指數(shù)分布具有更小的AIC或BIC值，因此q-指數(shù)分布具有更好的壽命數(shù)據(jù)擬合性。

表4 白血病數(shù)據(jù)集各模型的p值、AIC和BIC

4.2 元件壽命數(shù)據(jù)

如表5所示，這是由文獻(xiàn)[17]提供的50臺(tái)設(shè)備的壽命數(shù)據(jù)，已有研究人員利用韋伯分布[18-21]、指數(shù)分布[5]、LFR分布[22]分析了這個(gè)數(shù)據(jù)集。表6給出了所使用的每個(gè)分布參數(shù)的MLE估計(jì)值。圖6和圖7分別給出了設(shè)備數(shù)據(jù)集的經(jīng)驗(yàn)參數(shù)生存函數(shù)以及風(fēng)險(xiǎn)函數(shù)。

圖6 Aarset數(shù)據(jù)集的生存函數(shù)

圖7 Aarset數(shù)據(jù)集的風(fēng)險(xiǎn)函數(shù)

表5 50臺(tái)設(shè)備的生存期

表6 設(shè)備時(shí)間數(shù)據(jù)集各模型參數(shù)的MLE值

已知Aarset數(shù)據(jù)集的風(fēng)險(xiǎn)函數(shù)具有浴盆形狀。不失一般性和為了分析簡(jiǎn)便，本文只專注風(fēng)險(xiǎn)函數(shù)中單調(diào)遞增部分，使用雙參數(shù)分布來(lái)近似風(fēng)險(xiǎn)函數(shù)。為了進(jìn)行參數(shù)比較，使用似然檢驗(yàn)去對(duì)比原假設(shè)和備選假設(shè)。此外，利用AIC[23]在多個(gè)模型中選擇最優(yōu)模型。最適合數(shù)據(jù)擬合的模型應(yīng)具有最低的AIC。表7和表8給出了對(duì)數(shù)似然函數(shù)值、KS值、p值、AIC和BIC值[24]。

表7 設(shè)備時(shí)間數(shù)據(jù)集各模型的對(duì)數(shù)似然函數(shù)值和K-S值

表8 設(shè)備時(shí)間數(shù)據(jù)的P值，AIC和BIC

對(duì)于Aarset數(shù)據(jù)集，由表7可知q-指數(shù)分布的K-S值大于0.2，所以具備良好的分隔能力。由所列對(duì)數(shù)似然函數(shù)值可知，q-指數(shù)分布與韋伯分布、指數(shù)分布的對(duì)數(shù)似然函數(shù)值在相近水平，說(shuō)明q-指數(shù)分布亦可以充分利用先驗(yàn)信息獲得對(duì)未知參數(shù)的最大似然估計(jì)。

由表8知，q-指數(shù)分布的p值大于0.05，因而具備描述Aarset數(shù)據(jù)集的能力。同時(shí)，可以看出q-指數(shù)分布在本文所提到的所有分布中具有最小的AIC和最小的BIC值。這說(shuō)明在所列分布中，q-指數(shù)分布能夠最好地?cái)M合本數(shù)據(jù)集。

4.3 鋰電池壽命預(yù)測(cè)

估計(jì)剩余使用壽命(remaining useful life,RUL)有助于降低實(shí)際系統(tǒng)中發(fā)生災(zāi)難性事件的概率[25]。為了研究所提出的q-指數(shù)分布的有效性，本文采用美國(guó)宇航局(NASA) Prognostics公司(PCoE)的3個(gè)電池?cái)?shù)據(jù)集來(lái)預(yù)測(cè)鋰電池的剩余使用壽命。在該數(shù)據(jù)集中用電池容量來(lái)表征壽命，容量越低，剩余壽命越少。

3組電池(即蓄電池005、蓄電池006和蓄電池018)屬于同一類型，通過(guò)在室溫下工作在3種不同的狀態(tài)下(充電、放電和阻抗)進(jìn)行加速老化試驗(yàn)[26]。這種電池的額定容量是2 Ah，當(dāng)電池容量減少到額定容量的70%(從2 Ah減少到1.4 Ah)時(shí)，電池就會(huì)達(dá)到使用壽命終止(EOL)標(biāo)準(zhǔn)，容量數(shù)據(jù)如圖8所示。

圖8 NASA PCoE的電池容量數(shù)據(jù)

容量退化過(guò)程可用狀態(tài)空間模型來(lái)描述：

式中，xk表示k周期的真實(shí)容量值；yk表示k周期的預(yù)測(cè)值；wk-1表 示環(huán)境的擾動(dòng)；υk表示觀測(cè)噪聲。

分別采用指數(shù)分布和q-指數(shù)分布作為h(xk)來(lái)描述容量狀態(tài)轉(zhuǎn)換，采用粒子濾波算法[27]自適應(yīng)地確定預(yù)測(cè)容量值，并使用最大似然法確定指數(shù)分布和q-指數(shù)分布的相應(yīng)參數(shù)。

由粒子濾波基本原理可知，按照頻率派的概率觀點(diǎn)[24]，這里剩余壽命預(yù)測(cè)的概率分布可以通過(guò)對(duì)粒子濾波每個(gè)周期間隔內(nèi)的粒子個(gè)數(shù)計(jì)數(shù)得到。間隔內(nèi)粒子數(shù)越多，說(shuō)明該區(qū)間代表大多數(shù)情況下的壽命長(zhǎng)度。

在電池005實(shí)驗(yàn)中，共有168個(gè)循環(huán)樣本，分別用前60個(gè)、80個(gè)和100個(gè)樣本點(diǎn)訓(xùn)練粒子濾波器，設(shè)定粒子濾波預(yù)測(cè)的起始時(shí)間分別為T=60、80、100。當(dāng)起始點(diǎn)T=100時(shí)，RUL的分布直方圖和RUL預(yù)測(cè)結(jié)果如圖9、圖10和表9所示。其中，圖9b使用參數(shù)q=1.001。

圖9 電池005的RUL直方圖

此外，為了定量評(píng)價(jià)預(yù)測(cè)精度，將預(yù)測(cè)誤差定義如下：

式中，RULP表示預(yù)測(cè)的周期數(shù)；RULt表示真實(shí)的周期數(shù)。

如圖10(其中q=1.001)所示，電池005的真實(shí)壽命周期的結(jié)束數(shù)為124，使用q-指數(shù)分布相對(duì)能更好地?cái)M合實(shí)際值。根據(jù)表9，在指數(shù)分布假設(shè)下，在起始預(yù)測(cè)點(diǎn)T=60時(shí)，粒子濾波算法預(yù)測(cè)的壽命周期數(shù)為152。根據(jù)上述定義，指數(shù)分布模型預(yù)測(cè)誤差為：

表9 不同起始點(diǎn)T下鋰電池的RUL預(yù)測(cè)結(jié)果

圖10 電池005在起始點(diǎn)T為100的RUL預(yù)測(cè)結(jié)果

eRUL=RULp-RULt=152-124=28

同樣以T=60為起始預(yù)測(cè)點(diǎn)，q-指數(shù)分布在q=1.001的情況下，預(yù)測(cè)誤差為23，小于指數(shù)分布假設(shè)的結(jié)果。因此，提出的q-指數(shù)分布輔助粒子濾波算法對(duì)電池005的剩余使用壽命有更準(zhǔn)確的估計(jì)。

同理，觀察圖11(其中q=1.010)與圖12(其中q=0.990)中電池006和電池018的結(jié)果，結(jié)合表9，發(fā)現(xiàn)對(duì)于電池006而言，真實(shí)壽命周期的結(jié)束數(shù)為108，而在指數(shù)分布假設(shè)下，在起始預(yù)測(cè)點(diǎn)T=60、80、100時(shí)，使用q-指數(shù)分布獲得預(yù)測(cè)結(jié)果均與指數(shù)分布的相同，說(shuō)明基于q-指數(shù)分布的粒子濾波方法至少可以獲得和指數(shù)分布假設(shè)下一樣的估計(jì)精度。

圖11 電池006在起始點(diǎn)T=100的RUL預(yù)測(cè)結(jié)果

圖12 電池018在起始點(diǎn)T = 80處的RUL預(yù)測(cè)結(jié)果

對(duì)于電池018而言，真實(shí)壽命周期的結(jié)束數(shù)為96，在指數(shù)分布假設(shè)下，在起始預(yù)測(cè)點(diǎn)T=60時(shí)，預(yù)測(cè)的壽命周期數(shù)為92，預(yù)測(cè)誤差為-4。而q-指數(shù)分布在q=0.990的情況下，如表9所示的同一點(diǎn)的預(yù)測(cè)誤差為-2，誤差預(yù)測(cè)減小了50%。對(duì)于T=80和T=100的起始預(yù)測(cè)點(diǎn)，使用q-指數(shù)分布均可取得更準(zhǔn)確的效果。

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

本文基于最大熵方法，通過(guò)計(jì)算均值約束下最大Rényi熵，得到一種新的廣義指數(shù)分布：q-指數(shù)分布。本文對(duì)q-指數(shù)分布的統(tǒng)計(jì)特性進(jìn)行了分析，并給出了均值和各階矩的解析公式。為了便于應(yīng)用于可靠性分析，給出了基于q-指數(shù)分布可靠性的預(yù)測(cè)模型及對(duì)應(yīng)的生存函數(shù)和風(fēng)險(xiǎn)函數(shù)的解析表達(dá)式。指出可采用了兩種方法：極大似然估計(jì)法和信息似然估計(jì)法進(jìn)行雙參數(shù)估計(jì)。最后，結(jié)合醫(yī)學(xué)白血病患者壽命數(shù)據(jù)集、設(shè)備元件壽命數(shù)據(jù)集及鋰電池剩余壽命數(shù)據(jù)集進(jìn)行了驗(yàn)證，通過(guò)與韋伯分布、指數(shù)分布等常用壽命預(yù)測(cè)分布對(duì)比，驗(yàn)證了q-指數(shù)分布的有效性和估計(jì)精度的優(yōu)良性。下一步將挖掘q-指數(shù)分布在復(fù)雜系統(tǒng)建模中的高效應(yīng)用。