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

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

0 引言

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

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

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

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

1.1 模型介紹

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

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

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

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

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

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

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

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

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

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

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

2.1 赤池信息準則（AIC）

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

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

2.2 貝葉斯信息準則（BIC）

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

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

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

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

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

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

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

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

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

卡方擬合優(yōu)度檢驗以及J-B檢驗均在5%顯著水平下進行正態(tài)性檢驗，檢驗結(jié)果為0表示數(shù)據(jù)符合正態(tài)分布特性，檢驗結(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)度檢驗結(jié)果分別如圖1～4所示。

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

圖1 Pjl3檢驗結(jié)果

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

圖2 T6檢驗結(jié)果

圖3 Ps16檢驗結(jié)果

圖4 Wf檢驗結(jié)果

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

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

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

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

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

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

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

為了檢驗高斯混合模型聚類效果以及穩(wěn)態(tài)數(shù)據(jù)特征值計算結(jié)果的準確性，添 加300個N（1000，4），400個N（980，3）以及300個N（1030，2）共計1000個數(shù)據(jù)點的樣本作為發(fā)動機穩(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ù)。

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

6 結(jié)論

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

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

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

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

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