基于高斯混合模型的發(fā)動機(jī)穩(wěn)態(tài)數(shù)據(jù)特征值提取方法

哈 圣，徐 昊，唐 震，朱赤洲

（中國航發(fā)沈陽發(fā)動機(jī)研究所，沈陽 110015）

0 引言

航空發(fā)動機(jī)穩(wěn)態(tài)性能參數(shù)是表征發(fā)動機(jī)技術(shù)狀態(tài)的重要參量，對這些參數(shù)的準(zhǔn)確測量是發(fā)動機(jī)研制和使用過程都必須關(guān)注的重要問題之一。數(shù)據(jù)測量不可避免地受到環(huán)境因素和測試系統(tǒng)本身的影響，測得的數(shù)據(jù)中會存在噪聲[1-2]。隨著測量設(shè)備采樣率的提升，穩(wěn)態(tài)數(shù)據(jù)樣本點(diǎn)數(shù)量往往數(shù)以千計(jì)。有時(shí)試驗(yàn)需要的穩(wěn)態(tài)采樣時(shí)間較長，其樣本數(shù)量成倍增加，噪聲數(shù)據(jù)的樣本點(diǎn)數(shù)量也隨之增加。數(shù)據(jù)的充分利用與數(shù)據(jù)量增加帶來篩選難度增大的矛盾問題較為尖銳。由數(shù)據(jù)噪聲引入的異常值無論其數(shù)據(jù)顯著性強(qiáng)弱與否，在大量數(shù)據(jù)的累積效應(yīng)下將會影響性能評估結(jié)果的準(zhǔn)確性。因此，數(shù)據(jù)中的噪聲直接影響發(fā)動機(jī)穩(wěn)態(tài)性能的考核結(jié)果，不利于發(fā)動機(jī)研制，并且如何合理利用發(fā)動機(jī)穩(wěn)態(tài)數(shù)據(jù)進(jìn)行融合以反映發(fā)動機(jī)實(shí)際穩(wěn)態(tài)工作特性成為難點(diǎn)。

針對數(shù)據(jù)噪聲的剔除與融合方法在數(shù)據(jù)處理領(lǐng)域開展了大量應(yīng)用性研究，如模糊指數(shù)滑動平均濾波法對數(shù)據(jù)進(jìn)行濾波[3]以及基于模糊一致性矩陣的多源多模型加權(quán)決策融合診斷方法[4]。葉川等[1]對工程中常用的萊茵達(dá)準(zhǔn)則、狄克遜準(zhǔn)則等異常值剔除方法進(jìn)行較為詳細(xì)介紹，表明其評判標(biāo)準(zhǔn)在于置信區(qū)間的選擇；陳震宇等[5]基于分布圖法對航空發(fā)動機(jī)穩(wěn)態(tài)數(shù)據(jù)異常值進(jìn)行剔除，并采用分批估計(jì)方法進(jìn)行數(shù)據(jù)融合，提升了穩(wěn)態(tài)數(shù)據(jù)計(jì)算結(jié)果的可靠性。由于此類方法大多建立在數(shù)據(jù)源為單一正態(tài)分布特性基礎(chǔ)上，其理論依據(jù)主要是試驗(yàn)參數(shù)采用等精度傳感器測量，且測量形式符合中心極限定理使用情況[6-7]，通過對置信區(qū)間的選定來剔除異常值，適用于小樣本空間，若置信區(qū)間劃分不當(dāng)容易刪除有用的原始數(shù)據(jù)，從而影響數(shù)據(jù)信息整合。隨著機(jī)器學(xué)習(xí)的興起，許多專家與學(xué)者采用智能分類算法進(jìn)行異常值剔除，常見的有K-means聚類算法[8-10]、LOF離異值檢測算法等[11]。章永來等[12]介紹了關(guān)于K-means算法自20世紀(jì)60年代以來經(jīng)Bradley與Berkhin等在原算法的基礎(chǔ)上不斷改進(jìn)與優(yōu)化后，提升算法本身收斂速度并擴(kuò)展到了分布式聚類領(lǐng)域，使其成為應(yīng)用較廣、較為高效的聚類方法，廣泛應(yīng)用于對數(shù)據(jù)異常值的剔除。雖然該聚類方法屬于無監(jiān)督學(xué)習(xí)方法，但劃分方式大多依靠閔可夫斯基距離，且屬于硬聚類，發(fā)動機(jī)穩(wěn)態(tài)試驗(yàn)數(shù)據(jù)測量值的正態(tài)性無法很好地以概率方式表征，使得該方法對發(fā)動機(jī)穩(wěn)態(tài)數(shù)據(jù)異常值剔除的適用性不強(qiáng)。

大數(shù)據(jù)技術(shù)在航空發(fā)動機(jī)研發(fā)領(lǐng)域也愈發(fā)體現(xiàn)其價(jià)值性[13]。為此，本文采用期望極大（Exceptation Maximization algorithm，EM）算法[14-15]來求解高斯混合模型（Gaussian mivtnre model,GMM）。利用混合模型自身良好的回歸特性表征數(shù)據(jù)結(jié)構(gòu)，通過EM算法在逐步迭代過程中能夠近似實(shí)現(xiàn)對觀測數(shù)據(jù)的極大似然估計(jì)特性[16-17]，提出了發(fā)動機(jī)穩(wěn)態(tài)試驗(yàn)數(shù)據(jù)特征值的提取方法。

1 高斯混合模型及模型求解方法

1.1 模型介紹

高斯混合模型是指模型本身概率分布由多個(gè)高斯分布概率密度函數(shù)疊加而成的一種混合模型，其概率分布形式為

式中：K為分模型總數(shù)；y為樣本點(diǎn)；αk為第k個(gè)分模型系數(shù)為高斯分布密度函數(shù)，為 均值為 方 差 參 數(shù)，φ(y|θk)=稱為第k個(gè)分模型。

1.2 高斯混合模型參數(shù)估計(jì)的EM算法

EM算法是一種迭代算法，常用于含有隱參數(shù)變量的概率模型極大似然估計(jì)，或極大后驗(yàn)概率估計(jì)。EM算法在每次迭代過程中包含：E步，求解模型期望值；M步，求解模型極大似然估計(jì)值[16]。EM算法在工程上是一種求解高斯混合模型的有效方法。具體求解過程如下：

（1）選取模型參數(shù)的初始值進(jìn)行迭代；

（2）E步：依據(jù)當(dāng)前模型參數(shù)，計(jì)算分模型k對觀測數(shù)據(jù)yj的響應(yīng)度

（3）M步：計(jì)算新一輪迭代的模型參數(shù)

（4）重復(fù)（2）、（3）直到收斂，收斂后可認(rèn)為M步最終計(jì)算各分模型參數(shù)以及近似于實(shí)際分模型參數(shù)與αk。

2 聚類數(shù)的確定方法

由第1章可知，為了充分利用混合模型良好的回歸特性，需要確定高斯混合模型分模型數(shù)目，即確定數(shù)據(jù)聚類數(shù)。

本文采用赤池信息準(zhǔn)則（Akaike Information Criterion，AIC）準(zhǔn)則與貝葉斯信息準(zhǔn)則（Bayesian Information Criterion,）準(zhǔn)則的評估準(zhǔn)則來確定數(shù)據(jù)樣本的聚類數(shù)。用于評估真實(shí)模型與估計(jì)模型之間的K-L散度（Kullback-Leibler Divergence，KLD）中無偏估計(jì)項(xiàng)與GMM中參數(shù)數(shù)量漸近相等[18]，由此可分別利用AIC準(zhǔn)則在小樣本與BIC準(zhǔn)則在大樣本空間的模型定階優(yōu)勢作為GMM聚類數(shù)的確定方法。

2.1 赤池信息準(zhǔn)則（AIC）

AIC信息準(zhǔn)則[19-20]又稱赤池信息準(zhǔn)則，是日本統(tǒng)計(jì)學(xué)家赤池弘次（H.Akaike）在解決時(shí)間序列模型定階問題時(shí)，從隨機(jī)建模觀點(diǎn)出發(fā)針對信息論的研究提出來的基本信息量定階準(zhǔn)則。AIC準(zhǔn)則在統(tǒng)計(jì)分析特別是在統(tǒng)計(jì)模型中應(yīng)用廣泛，其計(jì)算結(jié)果越小，說明模型接近程度越高

式中：K為模型參數(shù)數(shù)量；L為模型極大似然函數(shù)值。

2.2 貝葉斯信息準(zhǔn)則（BIC）

與AIC準(zhǔn)則相同，BIC準(zhǔn)則在模型復(fù)雜性與模型對數(shù)據(jù)集的描述能力二者之間尋求最佳平衡。當(dāng)數(shù)據(jù)的樣本空間較大時(shí)，AIC準(zhǔn)則由于似然函數(shù)值過大，會削弱模型參數(shù)的影響；為解決大樣本空間下AIC準(zhǔn)則的不足問題，Schwarz[21]基于貝葉斯理論提出了BIC準(zhǔn)則，定義為

式中：N為樣本數(shù)據(jù)數(shù)。

與AIC信息準(zhǔn)則相比，BIC信息準(zhǔn)則左邊第1項(xiàng)用KlnN代替了2K，在大樣本空間中l(wèi)nN遠(yuǎn)大于2，因此在大樣本數(shù)據(jù)中使用極小化BIC定出的模型階數(shù)估計(jì)值要比AIC的低，彌補(bǔ)了AIC準(zhǔn)則在大樣本空間的不足，效果也更明顯。

3 穩(wěn)態(tài)數(shù)據(jù)正態(tài)性檢驗(yàn)

3.1 數(shù)據(jù)正態(tài)性檢驗(yàn)方法

航空發(fā)動機(jī)穩(wěn)態(tài)數(shù)據(jù)分布特性受測量方法以及測量形式影響，數(shù)據(jù)分布特性往往服從正態(tài)分布特性。為了檢驗(yàn)數(shù)據(jù)是否符合正態(tài)分布特性，本文穩(wěn)態(tài)數(shù)據(jù)的正態(tài)性檢驗(yàn)將分別采用主觀判斷法、卡方擬合優(yōu) 度 檢 驗(yàn)（Chi-Square Goodness-of-Fit Test）和（Jarque-Bera，J-B）檢驗(yàn)[22]相結(jié)合的方法。其中，主觀判斷法采用數(shù)據(jù)頻數(shù)直方圖以及Q-Q圖（Quantile-Quantile chart）正態(tài)性檢驗(yàn)法。Q-Q圖正態(tài)性檢驗(yàn)法是將樣本得分位數(shù)與按照正態(tài)分布計(jì)算得到的分位數(shù)作為坐標(biāo)軸，如果樣本服從正態(tài)分布，則樣本應(yīng)呈1條圍繞第一象限對角線的直線型散點(diǎn)?？紤]到穩(wěn)態(tài)數(shù)據(jù)的分布特性整體與局部差異性，需要對穩(wěn)態(tài)數(shù)據(jù)整體性驗(yàn)證后再進(jìn)行局部數(shù)據(jù)檢驗(yàn)。由于Jarque-Bera檢驗(yàn)基于偏度和峰度進(jìn)行正態(tài)性檢驗(yàn)，檢驗(yàn)結(jié)果受異常值影響過大，因此對于局部數(shù)據(jù)正態(tài)性檢驗(yàn)僅采用卡方擬合優(yōu)度檢驗(yàn)法。

3.2 正態(tài)性檢驗(yàn)示例

為了充分體現(xiàn)數(shù)據(jù)采集系統(tǒng)不同來源數(shù)據(jù)特性，分別對電壓信號輸出的單精度壓力傳感器、數(shù)字量數(shù)據(jù)輸出設(shè)備溫度采集模塊DTS3250和壓力采集模塊DSA3217以及頻率量輸出信號的采集數(shù)據(jù)進(jìn)行數(shù)據(jù)特性檢驗(yàn)，選擇穩(wěn)態(tài)數(shù)據(jù)參數(shù)Pjl3、T6、Ps16，Wf作為此次檢驗(yàn)參數(shù)，其中，Pjl3為某加力區(qū)供油壓力、T6為低壓渦輪后溫度、Ps16為外涵出口靜壓、Wf為頻率輸出的流量計(jì)測得流經(jīng)發(fā)動機(jī)燃油流量，為直接測量參數(shù)。為了縮小數(shù)量相對關(guān)系以及消除量綱的影響[23]，對穩(wěn)態(tài)數(shù)據(jù)進(jìn)行標(biāo)準(zhǔn)分?jǐn)?shù)歸一化方法（z-score normalization，ZSN），轉(zhuǎn)化函數(shù)為

式中：x'為轉(zhuǎn)換后數(shù)據(jù)；x、μ、σ分別為轉(zhuǎn)換前原始數(shù)據(jù)、樣本均值、標(biāo)準(zhǔn)差。

卡方擬合優(yōu)度檢驗(yàn)以及J-B檢驗(yàn)均在5%顯著水平下進(jìn)行正態(tài)性檢驗(yàn)，檢驗(yàn)結(jié)果為0表示數(shù)據(jù)符合正態(tài)分布特性，檢驗(yàn)結(jié)果為1表示數(shù)據(jù)不符合正態(tài)分布特性。Pjl3、T6、Ps16、Wf的穩(wěn)態(tài)數(shù)據(jù)整體結(jié)果見表1，各參數(shù)的頻數(shù)直方圖、Q-Q圖以及局部數(shù)據(jù)每組卡方擬合優(yōu)度檢驗(yàn)結(jié)果分別如圖1～4所示。

從表1和圖1～4中可見，發(fā)動機(jī)穩(wěn)態(tài)數(shù)據(jù)大體符合正態(tài)分布特性，雖然局部檢驗(yàn)結(jié)果存在不符合情況，但所占比例較小，其檢驗(yàn)結(jié)果受劃分?jǐn)?shù)量影響，因此可認(rèn)為發(fā)動機(jī)穩(wěn)態(tài)數(shù)據(jù)分布特性符合正態(tài)分布。

圖1 Pjl3檢驗(yàn)結(jié)果

表1 穩(wěn)態(tài)數(shù)據(jù)整體性正態(tài)檢驗(yàn)結(jié)果

圖2 T6檢驗(yàn)結(jié)果

圖3 Ps16檢驗(yàn)結(jié)果

圖4 Wf檢驗(yàn)結(jié)果

4 穩(wěn)態(tài)數(shù)據(jù)特征值計(jì)算方法

理論上在航空發(fā)動機(jī)穩(wěn)定工作狀態(tài)，其轉(zhuǎn)子轉(zhuǎn)速、渦輪前后的燃?xì)鉁囟?、耗油量等參?shù)不隨時(shí)間發(fā)生變化[24]，但實(shí)際工作狀態(tài)受控制規(guī)律、工作環(huán)境等多方面因素影響，其狀態(tài)值會在一定范圍內(nèi)波動，有時(shí)會使數(shù)據(jù)本身不再滿足單一正態(tài)分布特性。從時(shí)間域上看，可將這種波動視為若干個(gè)穩(wěn)定工作點(diǎn)間的數(shù)值切換，即測量數(shù)值為若干正態(tài)分布數(shù)據(jù)值與噪聲的混合疊加。

根據(jù)在試驗(yàn)過程中噪聲數(shù)據(jù)樣本點(diǎn)在全樣本空間中的統(tǒng)計(jì)特性，可假設(shè)在整段穩(wěn)態(tài)測量中出現(xiàn)非高斯噪聲為小概率事件，其高斯混合模型所占權(quán)重則遠(yuǎn)小于各穩(wěn)定工作點(diǎn)的分模型權(quán)重，那么穩(wěn)態(tài)數(shù)據(jù)的特征值可寫成

整個(gè)計(jì)算過程包括對選取穩(wěn)態(tài)數(shù)據(jù)片段在高斯混合模型下最優(yōu)劃分形式的求解以及對篩選數(shù)據(jù)的融合計(jì)算過程，計(jì)算過程如圖5所示。

圖5 穩(wěn)態(tài)數(shù)據(jù)特征值計(jì)算流程

5 數(shù)據(jù)驗(yàn)證與結(jié)果分析

5.1 仿真計(jì)算與結(jié)果分析

為了檢驗(yàn)高斯混合模型聚類效果以及穩(wěn)態(tài)數(shù)據(jù)特征值計(jì)算結(jié)果的準(zhǔn)確性，添 加300個(gè)N（1000，4），400個(gè)N（980，3）以及300個(gè)N（1030，2）共計(jì)1000個(gè)數(shù)據(jù)點(diǎn)的樣本作為發(fā)動機(jī)穩(wěn)態(tài)數(shù)據(jù)原始樣本，在原始樣本中分別添加正弦與脈沖噪聲，添加后的數(shù)據(jù)形式如圖6所示。

圖6 數(shù)據(jù)分布形式

圖中，折線1為在原始數(shù)據(jù)上添加的脈沖信號，折線3、5為在原始數(shù)據(jù)上添加的10sin（300t）+15以及40sin（300t）噪聲數(shù)據(jù)。

通過計(jì)算不同聚類數(shù)的高斯混合模型的AIC與BIC準(zhǔn)則來確定分模型的數(shù)量。聚類數(shù)K與AIC、BIC計(jì)算值的關(guān)系見表2。

表2 不同聚類數(shù)AIC與BIC計(jì)算結(jié)果

根據(jù)大樣本數(shù)據(jù)BIC最小原則，確定聚類數(shù)為5，得到聚類效果散點(diǎn)，如圖7所示。

圖7 數(shù)據(jù)聚類效果散點(diǎn)（K=5）

各分模型計(jì)算結(jié)果及分類情況見表3。

表3 數(shù)據(jù)分類結(jié)果（K=5）

根據(jù)穩(wěn)態(tài)特性數(shù)據(jù)分布相似的基本原則，選取方差量級相同的高斯分布類數(shù)據(jù)作為原始樣本，按照表3中穩(wěn)態(tài)數(shù)據(jù)特征值計(jì)算方法，選取第1～3類均值與權(quán)重參與計(jì)算，計(jì)算結(jié)果為1001.033，具體計(jì)算過程為

由已知原始分布數(shù)據(jù)可知原始實(shí)際值為1001.000，計(jì)算所得數(shù)據(jù)特征值相對誤差為0.033%。

5.2 穩(wěn)態(tài)數(shù)據(jù)驗(yàn)證結(jié)果與分析

為了驗(yàn)證高斯混合模型對真實(shí)發(fā)動機(jī)穩(wěn)態(tài)數(shù)據(jù)的降噪效果以及穩(wěn)態(tài)數(shù)據(jù)時(shí)間片段選取不同對計(jì)算結(jié)果的差異性進(jìn)行評估，對某一參數(shù)穩(wěn)態(tài)數(shù)據(jù)波動較大的數(shù)據(jù)片段進(jìn)行模型驗(yàn)證，該原始數(shù)據(jù)在不同聚類數(shù)下AIC與BIC值計(jì)算結(jié)果見表4。

表4 原始數(shù)據(jù)聚類結(jié)果

從表中可見，當(dāng)K=2，3時(shí)模型結(jié)果劃分較優(yōu)，原始數(shù)據(jù)劃分類別為2和3的計(jì)算結(jié)果見表5，劃分結(jié)果如圖8所示。

表5 K=2、3的數(shù)據(jù)計(jì)算結(jié)果

圖8 原始數(shù)據(jù)劃分結(jié)果

考慮穩(wěn)態(tài)數(shù)據(jù)片段選擇對數(shù)據(jù)穩(wěn)態(tài)特征值計(jì)算結(jié)果的影響，從總的穩(wěn)態(tài)數(shù)據(jù)中隨機(jī)選取2組不同時(shí)間片段數(shù)據(jù)進(jìn)行模型分類驗(yàn)證，各片段參數(shù)計(jì)算結(jié)果見表6，2組數(shù)據(jù)最優(yōu)劃分結(jié)果如圖9所示。

圖9 不同穩(wěn)態(tài)時(shí)間片段聚類效果

表6 2組數(shù)據(jù)片段在K=2時(shí)的計(jì)算結(jié)果

從圖7中的劃分結(jié)果可知，經(jīng)高斯混合模型進(jìn)行數(shù)據(jù)篩選后的穩(wěn)態(tài)數(shù)據(jù)波動水平降低，可達(dá)到降噪效果。

從第4章中的穩(wěn)態(tài)數(shù)據(jù)特征值計(jì)算方法可知總體穩(wěn)態(tài)數(shù)據(jù)以及2組穩(wěn)態(tài)數(shù)據(jù)片段計(jì)算結(jié)果均為K=2時(shí)權(quán)重最大類的均值，片段1與總體穩(wěn)態(tài)數(shù)據(jù)特征值計(jì)算結(jié)果相差0.0095%，而片段2與總體穩(wěn)態(tài)數(shù)據(jù)特征值計(jì)算結(jié)果相差0.2%，片段1、2的整體差異性不大，計(jì)算結(jié)果與發(fā)動機(jī)穩(wěn)定工作狀態(tài)數(shù)據(jù)性能參數(shù)不隨時(shí)間變化的理論特性相符。

6 結(jié)論

（1）高斯混合模型可以有效篩選出具有正態(tài)分布特性的穩(wěn)態(tài)數(shù)據(jù)；

（2）穩(wěn)態(tài)數(shù)據(jù)提取方法所得計(jì)算結(jié)果具有較高的準(zhǔn)確性；

（3）基于高斯混合模型的數(shù)據(jù)篩選方法可以達(dá)到穩(wěn)態(tài)數(shù)據(jù)降噪的效果；

（4）穩(wěn)態(tài)數(shù)據(jù)片段的選取對穩(wěn)態(tài)數(shù)據(jù)特征值計(jì)算結(jié)果影響不大，該計(jì)算結(jié)果可以作為發(fā)動機(jī)真實(shí)的穩(wěn)態(tài)性能參數(shù)指標(biāo)。

由于模型求解采用EM算法，本身僅具有局部收斂性，在數(shù)據(jù)層次較復(fù)雜模型以及設(shè)定的分模型數(shù)目較多時(shí)，可能很難收斂至實(shí)際值。目前，工程上大多采用在給定分模型數(shù)目時(shí)多次計(jì)算以最小的AIC或BIC值作為最終收斂結(jié)果。