申世才，雷 杰，郝曉樂

(中國飛行試驗研究院發(fā)動機所，西安，710089)

喘振是一種危害性極大的發(fā)動機壓縮系統(tǒng)氣動失穩(wěn)現(xiàn)象，輕則造成發(fā)動機工況急劇惡化,重則導致發(fā)動機機械損傷，進而引發(fā)嚴重飛行事故，引起人員和財產(chǎn)的巨大損失[1]。因此，尋求一種快速、準確檢測發(fā)動機喘振故障的方法是發(fā)動機穩(wěn)定性研究領(lǐng)域的熱點問題，對于保障航空器飛行安全具有重要的現(xiàn)實意義和價值。

當前，國內(nèi)外的喘振檢測方法基本是以壓氣機出口流場壓力為對象，根據(jù)壓力信號時域或頻域特征進行喘振算法研究。其中基于壓力信號時域特征的喘振檢測算法主要有：自相關(guān)函數(shù)法[2]、短時能量法[3]、變化率法[4-5]、壓差法[6]、方差分析法[7]、統(tǒng)計特征法[8]等，基于壓力信號頻域特征的喘振檢測算法主要有：頻譜分析法[9-11]、小波分析法[12-13]、頻域幅值法[14-15]等。上述時域和頻域特征喘振檢測算法多數(shù)是以壓氣機部件級試驗為支撐，取得了很多較好的喘振檢測效果，但由于受到試驗對象和條件限制，應(yīng)用在發(fā)動機整機喘振檢測的算法相對較少，有必要尋求一種簡單、可靠、準確的發(fā)動機整機喘振檢測方法。

1 壓氣機出口脈動壓力統(tǒng)計特征分析

發(fā)動機壓氣機出口動態(tài)壓力p31包含平均壓力成分p31.a和脈動壓力成分Δp31，在發(fā)動機穩(wěn)定狀態(tài)下，壓力機內(nèi)部流場穩(wěn)定，其出口脈動壓力Δp31相對穩(wěn)定，而在發(fā)動機出現(xiàn)喘振時，壓氣機出口流場Δp31呈不穩(wěn)定狀態(tài)[4]。所以當發(fā)動機穩(wěn)定工作時，如果Δp31服從某種概率分布規(guī)律，根據(jù)其概率分布函數(shù)，就可以通過Δp31的統(tǒng)計特征設(shè)置閾值來檢測發(fā)動機喘振故障的發(fā)生。

正態(tài)分布是統(tǒng)計學中最重要的一種分布規(guī)律，大量隨機現(xiàn)象可以用正態(tài)分布規(guī)律來描述或近似，同時正態(tài)分布具有很多優(yōu)良的性質(zhì)，所以不論是在理論研究還是工程實踐中，正態(tài)分布具有廣泛的應(yīng)用。這里首先假設(shè)Δp31測量數(shù)值服從正態(tài)分布，下面采用某發(fā)動機實測數(shù)據(jù)檢驗壓氣機出口脈動壓力時間序列數(shù)據(jù)是否符合正態(tài)概率分布規(guī)律。

某發(fā)動機試飛期間，測量了不同飛行高度、速度條件下發(fā)動機穩(wěn)態(tài)及瞬態(tài)的p31時間序列數(shù)據(jù)，獲取了相應(yīng)的Δp31時間序列數(shù)據(jù)。將這些不同飛行條件下發(fā)動機穩(wěn)態(tài)和瞬態(tài)Δp31時間序列數(shù)據(jù)作為總體統(tǒng)計量，采用簡單隨機抽樣原則從總體統(tǒng)計量中抽取發(fā)動機穩(wěn)態(tài)和瞬態(tài)Δp31時間序列數(shù)據(jù)作為檢驗樣本，進行Δp31時間序列正態(tài)符合性檢驗分析。

樣本數(shù)據(jù)x的正態(tài)性檢驗方法有很多，主要分為統(tǒng)計圖法和統(tǒng)計指標法。統(tǒng)計圖法包括直方圖、P-P圖、箱式圖、概率圖、分位數(shù)圖等，均能為樣本正態(tài)性提供一個粗略估計，本文采用分位數(shù)圖對樣本Δp31正態(tài)性進行初步估計；統(tǒng)計指標法可對樣本正態(tài)性進行定量檢驗，主要包括偏度峰度聯(lián)合檢驗、S-W(Shapiro-Wilk)檢驗、A-D(Anderson-Darling)檢驗、K-S(Kolmogorov-Smirnov)檢驗等，本文選擇K-S檢驗作為樣本Δp31的正態(tài)性定量檢驗方法。

圖1為高壓轉(zhuǎn)速為n2時發(fā)動機穩(wěn)態(tài)過程Δp31時間序列隨機樣本分位數(shù)圖，可見該發(fā)動機穩(wěn)態(tài)輸入樣本散點不完全分布在直線附近，故而不能完全為正態(tài)分布提供粗略支持，進一步采用K-S檢驗。

圖1 穩(wěn)態(tài)過程時間序列樣本數(shù)據(jù)分位數(shù)圖

K-S檢驗統(tǒng)計量為：

Dn=max{|Fn(xi)-F0(xi)|,|Fn(xi+1)-F0(xi)|}

(1)

式中：Fn(x)為樣本的概率分布函數(shù)；F0(x)為標準正態(tài)分布函數(shù)。如果檢驗統(tǒng)計量Dn大于給定顯著水平α和樣本容量n確定的檢驗臨界值Dn(α)，則拒絕零假設(shè)，否則接受零假設(shè)，認為Δp31符合正態(tài)分布。

對圖1中發(fā)動機穩(wěn)態(tài)過程時間序列樣本Δp31進行計算，得到Dn=0.008 3，取顯著性水平α為0.01時，查表得到臨界值Dn(α)=0.002 5，Dn>Dn(α)，故拒絕零假設(shè)，所以發(fā)動機穩(wěn)態(tài)過程Δp31不服從正態(tài)分布。

進一步分析發(fā)動機瞬態(tài)過程Δp31是否服從正態(tài)分布，圖2為發(fā)動機瞬態(tài)過程Δp31隨機樣本分位數(shù)圖，顯然發(fā)動機瞬態(tài)過程Δp31輸入樣本散點同樣不完全分布在直線附近，進一步進行K-S檢驗：Dn等于0.014 7，取顯著性水平α為0.01時，查表Dn(α)為0.004，Dn>Dn(α)，拒絕零假設(shè)，故發(fā)動機瞬態(tài)過程Δp31同樣也不服從正態(tài)分布。

圖2 瞬態(tài)過程時間序列樣本數(shù)據(jù)分位數(shù)圖

進行全部Δp31時間序列樣本正態(tài)符合性檢驗后，發(fā)現(xiàn)Δp31不服從正態(tài)分布，這是因為發(fā)動機壓氣機是高溫、高速旋轉(zhuǎn)機械部件，其Δp31信號包含了發(fā)動機轉(zhuǎn)子頻率、葉片頻率的基頻和倍頻成分及受發(fā)動機瞬態(tài)過程影響，其測量數(shù)據(jù)隨時間統(tǒng)計序列不服從正態(tài)分布，所以對于Δp31這種未知分布規(guī)律，難以根據(jù)其統(tǒng)計特征設(shè)置喘振檢測閾值，想要Δp31按正態(tài)分布設(shè)置喘振檢測閾值，必須將Δp31時間序列樣本數(shù)據(jù)進行正態(tài)轉(zhuǎn)換。

2 壓氣機出口脈動壓力正態(tài)轉(zhuǎn)換

2.1 Johnson正態(tài)轉(zhuǎn)換方法

常用的非正態(tài)數(shù)據(jù)正態(tài)轉(zhuǎn)換方法主要包括：倒數(shù)轉(zhuǎn)換、對數(shù)轉(zhuǎn)換、平方根轉(zhuǎn)換、平方根反正旋轉(zhuǎn)換、平方根后取倒數(shù)轉(zhuǎn)換、Johnson轉(zhuǎn)換等，其中前5種方法對樣本數(shù)據(jù)的要求較高，要求樣本x為正值，由于Δp31存在負值數(shù)據(jù)，故難以采用前5種轉(zhuǎn)換方法對Δp31數(shù)據(jù)進行正態(tài)轉(zhuǎn)換。

Johnson轉(zhuǎn)換條件相對寬松，允許樣本x存在負數(shù)的情況，具有良好適應(yīng)性且Johnson轉(zhuǎn)換可將隨機樣本x直接轉(zhuǎn)換為標準正態(tài)分布N(0,1)。本文采用Johnson轉(zhuǎn)換方法對Δp31時間序列樣本數(shù)據(jù)進行正態(tài)變換。樣本x進行Johnson正態(tài)轉(zhuǎn)換時，其基本轉(zhuǎn)換公式為：

z=γ+ηf(x,ε,λ)

(2)

式(2)中f(x,ε,λ)為特定Johnson分布曲線函數(shù)，具體為SB(有界)、SL(對數(shù)正態(tài))、SU(無界)3種類型，η、γ、λ、ε等通過對應(yīng)轉(zhuǎn)化類型計算公式得到，3種類型函數(shù)及參數(shù)約束條件見表1。

確定Johnson轉(zhuǎn)換函數(shù)的方法有很多，其中樣本百分位數(shù)法得到了廣泛應(yīng)用，本文采用百分位數(shù)法進行正態(tài)轉(zhuǎn)換。

首先選擇一個合適的r,在標準正態(tài)分布表中查找(-sr，-r，r，sr)的分布概率P-sr、P-r、Pr、Psr，通過樣本獲取對應(yīng)分位數(shù)x-sr，x-r，xr，xsr后，計算m=xsr-xr、k=x-r-x-sr、l=xr-x-r，得到分位數(shù)比mn/l2。

當mk/l2<1時為SB分布，mk/l2=1時為SL分布，mk/l2>1時為SU分布。由于分位數(shù)比與s和r取值有關(guān)，需要確定合適的s和r，研究表明s取值為3及r取值范圍為(0.25，0.26，…，1.25)時，可獲得較好的轉(zhuǎn)換效果[16-17]。

本文取r值為1.2(在0.25～1.25范圍內(nèi)，r取值1.2時正態(tài)轉(zhuǎn)換效率最高)，對于圖1和圖2中發(fā)動機穩(wěn)態(tài)和瞬態(tài)過程隨機樣本，計算得到mk/l2值分別為0.721和0.921 5，即mk/l2<1，Δp31的Johnson轉(zhuǎn)換函數(shù)為SB型，對應(yīng)的η、γ、λ、ε計算公式為:

(3)

(4)

(5)

(6)

2.2 Johnson正態(tài)轉(zhuǎn)換結(jié)果

采用K-S方法對正態(tài)轉(zhuǎn)換后的脈動壓力時間序列樣本Δp31.norm正態(tài)性符合性進行檢驗。其中圖1所示發(fā)動機穩(wěn)態(tài)過程Δp31.norm檢測結(jié)果為：Dn=0.001 1，Dn(α)=0.002 5，DnDn(α)，拒絕了零假設(shè)，故Johnson正態(tài)轉(zhuǎn)換后的瞬態(tài)Δp31.norm不服從正態(tài)分布。

對其他隨機樣本采用Johnson轉(zhuǎn)換并檢測轉(zhuǎn)換后Δp31.norm的正態(tài)性，發(fā)現(xiàn)Johnson轉(zhuǎn)換方法對于穩(wěn)態(tài)過程具有較高的轉(zhuǎn)換效率，穩(wěn)態(tài)過程Δp31.norm正態(tài)轉(zhuǎn)化成功率P可達到90%以上，但是對于瞬態(tài)過程隨機樣本正態(tài)轉(zhuǎn)換效果較差，瞬態(tài)過程Δp31.norm正態(tài)轉(zhuǎn)化成功率P不足50%。

進一步分析瞬態(tài)過程Δp31正態(tài)轉(zhuǎn)化率不足的原因：樣本Δp31是時間序列的統(tǒng)計量，其受發(fā)動機狀態(tài)變化影響很大，當發(fā)動機穩(wěn)態(tài)工作時，Δp31幅值范圍基本穩(wěn)定在同一量級，而不同轉(zhuǎn)速時Δp31幅值范圍差異較大，且在瞬態(tài)變化過程中，隨著發(fā)動機轉(zhuǎn)速的提高，Δp31幅值變化范圍也逐漸升高(具體如圖3所示)，這種Δp31幅值隨轉(zhuǎn)速變化的特性是導致瞬態(tài)過程Δp31隨機樣本Johnson正態(tài)轉(zhuǎn)化成功率P過低的主要原因。

為解決發(fā)動機瞬態(tài)過程Δp31正態(tài)轉(zhuǎn)換成功率P過低問題，必須降低瞬態(tài)過程轉(zhuǎn)速變化對于Δp31幅值變化范圍的影響，考慮到樣本Δp31是時間序列的高頻采集信號，其采樣頻率為5 kHz，減小Δp31樣本容量n就可以直接削弱轉(zhuǎn)速變化對于Δp31幅值變化范圍的影響，當樣本容量n減小到一定程度時(如幾十毫秒量級)，可假設(shè)認為發(fā)動機瞬態(tài)過程時間序列Δp31相當于一個準穩(wěn)態(tài)變化過程。按照上述方法，逐步減小隨機樣本容量n，從總體樣本抽取幾十組發(fā)動機瞬態(tài)過程隨機樣本進行Johnson正態(tài)轉(zhuǎn)換，正態(tài)轉(zhuǎn)化成功率P統(tǒng)計情況見表2所示，可見減小樣本容量n可有效提高瞬態(tài)過程Δp31的正態(tài)轉(zhuǎn)化成功率，當隨機樣本量n降低至500及以下，隨機樣本Δp31的Johnson正態(tài)轉(zhuǎn)換成功率達到90%以上，因此要實現(xiàn)瞬態(tài)過程Δp31有效正態(tài)轉(zhuǎn)換，可將樣本容量n確定為100。

表2 瞬態(tài)過程正態(tài)轉(zhuǎn)化成功率統(tǒng)計表

3 喘振檢測閾值和方法

3.1 喘振檢測量

正態(tài)轉(zhuǎn)換后的Δp31.norm～N(0,1)，根據(jù)其概率分布函數(shù)可知，Δp31.norm時間序列落在±3.5范圍內(nèi)的概率為99.95%，即Δp31.norm有99.95%概率分布在-3.5≤Δp31.norm≤3.5范圍內(nèi)，那么根據(jù)表1中SB型Johnson轉(zhuǎn)換函數(shù)的反函數(shù)z-1，可以得到正態(tài)轉(zhuǎn)換前Δp31的99.95%概率分布范圍為：

(7)

由于瞬態(tài)過程Δp31樣本容量n不宜過大，當n過大時將不能有效實現(xiàn)Δp31正態(tài)轉(zhuǎn)換，所以也就不能通過Johnson轉(zhuǎn)換函數(shù)的反函數(shù)z-1獲取Δp31的99.95%概率分布范圍。此外，不同發(fā)動機狀態(tài)下Δp31幅值范圍存在較大差異，導致其99.95%概率分布范圍不同，進而不能根據(jù)Δp31分布范圍使用固定喘振檢測閾值。為解決這一問題，本文采用一種基于滑動窗口的Johnson轉(zhuǎn)換方法，具體如下。

假設(shè)滑動窗口長度為d，將測取的發(fā)動機時間序列Δp31第1點數(shù)據(jù)x1作為滑動窗口起始點，向后截取d個Δp31數(shù)據(jù)。對窗口內(nèi)d個Δp31數(shù)據(jù)進行Johnson正態(tài)轉(zhuǎn)化，得到SB型Johnson正態(tài)轉(zhuǎn)換函數(shù)反函數(shù)的η、γ、λ、ε等參數(shù)，根據(jù)公式(7)計算窗口內(nèi)Δp31數(shù)據(jù)99.95%概率分布上界Δp31.max和下界Δp31.min。將Δp31.max、Δp31.min作為第xd點時間序列Δp31的99.95%概率分布上界Δp31.max(xd)和下界Δp31.min(xd)。按照Δp31時間序列順序滑動計算窗口，重復上述計算方法，可以得到Δp31的99.95%概率分布的自適應(yīng)上界Δp31.max和下界Δp31.min，實際等效于Δp31上下“包絡(luò)線”。

樣本容量n為100時正態(tài)轉(zhuǎn)換成功率最高，取滑動窗口d=100，對于圖3所示Δp31時間序列樣本數(shù)據(jù)，其99.95%概率的Δp31.max、Δp31.min分布情況見圖4所示，局部放大情況見圖5，可見采用滑動窗口的Johnson正態(tài)轉(zhuǎn)換方法，能夠獲取時間序列Δp31的99.95%概率分布自適應(yīng)上界Δp31.max和下界Δp31.min。

上文所述Δp31等效“包絡(luò)線”反映了其幅值變化范圍，可以根據(jù)Δp31上下“包絡(luò)線”間距D檢測發(fā)動機喘振，為消除D受發(fā)動機轉(zhuǎn)速變化的影響，定義一種無量綱的發(fā)動機喘振檢測量T31為：

(8)

圖4 脈動壓力概率邊界時間序列

圖5 脈動壓力概率邊界時間序列局部放大圖

式(8)等效于對Δp31.max和Δp31.min間距D進行歸一化處理，可以消除發(fā)動機不同轉(zhuǎn)速時D幅值范圍的差異，便于設(shè)置固定的喘振檢測閾值。如對于圖4中所示的Δp31時間序列，其對應(yīng)的喘振檢測量T31計算結(jié)果如圖6所示。

圖6 發(fā)動機無量綱喘振檢測量時間序列

3.2 喘振檢測閾值及方法

設(shè)置發(fā)動機喘振檢測閾值時，將上文不同飛行條件下發(fā)動機穩(wěn)定工作時測取的穩(wěn)態(tài)和瞬態(tài)時間序列Δp31作為總體統(tǒng)計量，采用滑動窗口的Johnson轉(zhuǎn)換方法計算對應(yīng)的喘振檢測量T31，將獲取的T31作為設(shè)置喘振檢測閾值A(chǔ)的總體樣本。

采用K-S方法檢驗總體樣本T31正態(tài)性，表明其并不服從正態(tài)分布。對總體樣本T31進行Johnson轉(zhuǎn)換，即T31.norm～N(0,1)。對于轉(zhuǎn)換后的T31.norm，其同樣有99.95%概率分布在±3.5范圍內(nèi)，由公式(7)計算發(fā)動機穩(wěn)定工作時T31變化范圍為：

0.018 5≤T31≤0.079 8

(9)

因為當發(fā)動機喘振時，Δp31脈動幅值會迅速增大，其概率意義“包絡(luò)線”間距D必然會出現(xiàn)劇增，那么T31將不可避免超出0.079 8的范圍，可以設(shè)置固定喘振檢測閾值A(chǔ)=0.079 8。但是發(fā)動機穩(wěn)定工作時，T31還存在0.025%概率小幅超過0.079 8，所以綜合考慮喘振檢測靈敏度和誤報率，引入一個優(yōu)化參數(shù)t(1≤t≤2)，設(shè)置保守的固定喘振檢測閾值為：A=0.079 8t，t根據(jù)實際喘振檢測效果設(shè)定，以最低喘振誤報率作為t取值依據(jù)。

綜上所述，根據(jù)發(fā)動機固定喘振檢測閾值A(chǔ)，本文設(shè)計的喘振檢測方法為：當喘振檢測量T31大于A時，認為發(fā)動機出現(xiàn)喘振；當T31小于等于A時，認為發(fā)動機工作穩(wěn)定。

4 喘振檢測方法試驗驗證

某發(fā)動機試飛期間先后發(fā)生3起喘振故障，其中地面發(fā)生喘振1次，空中發(fā)生喘振2次(記為空中喘振1、空中喘振2)，地面發(fā)動機喘振發(fā)生在穩(wěn)態(tài)過程中，空中2起喘振分別發(fā)生在加速、減速過程中，使用該型發(fā)動機喘振試飛數(shù)據(jù)對喘振檢測方法的實際效果進行試驗驗證。

圖7～9為該發(fā)動機3次喘振檢測驗證結(jié)果，3次喘振檢測均發(fā)出了喘振信號。當發(fā)動機穩(wěn)定工作時，不管是穩(wěn)態(tài)過程還是瞬態(tài)過程，喘振檢測量T31變化幅值均相對穩(wěn)定，不受發(fā)動機轉(zhuǎn)速變化影響，T31沒有超過檢測閾值A(chǔ)的范圍，而當發(fā)動機出現(xiàn)喘振時，喘振檢測量T31值會出現(xiàn)劇增并超過檢測閾值A(chǔ)。圖7所示發(fā)動機穩(wěn)態(tài)過程，喘振時，從脈動壓力Δp31出現(xiàn)波動到發(fā)出喘振信號用時約為12 ms，喘振期間T31最大值增長至1.748；喘振消失時，從脈動壓力Δp31停止波動到喘振信號消失用時約為14 ms。圖8、圖9所示發(fā)動機瞬態(tài)過程，喘振時從脈動壓力Δp31出現(xiàn)大幅波動到發(fā)出喘振信號用時分別約為16 ms、15 ms，喘振期間T31最高增長到1.653和1.689；喘振消失時從脈動壓力Δp31停止波動到喘振信號消失用時分別約為17 ms、14 ms。

綜上所述，本文發(fā)動機喘振檢測方法成功檢測出3次喘振故障，其識別發(fā)動機喘振及退出喘振所需時間均小于20 ms，未出現(xiàn)虛警、漏報等異常情況，驗證了檢測方法的有效性和準確性，表明方法具有識別率高、報警遲滯小等優(yōu)點。

圖7 發(fā)動機地面喘振檢測結(jié)果

圖8 發(fā)動機空中喘振1檢測結(jié)果

圖9 發(fā)動機空中喘振2檢測結(jié)果

5 結(jié)論

1)在地面及空中發(fā)動機穩(wěn)定工作時，壓氣機出口脈動壓力不服從正態(tài)分布。對于穩(wěn)態(tài)過程脈動壓力，采用Johnson方法的正態(tài)轉(zhuǎn)換成功率較高，可達90%以上，但是對于瞬態(tài)過程正態(tài)轉(zhuǎn)換效果較差，當樣本容量n降低至500以下時，瞬態(tài)過程脈動壓力正態(tài)轉(zhuǎn)化成功率可達到90%以上。

2)采用滑動窗口的Johnson轉(zhuǎn)換方法，可以獲取壓氣機出口脈動壓力99.95%概率分布的自適應(yīng)上邊界和下邊界。在不同發(fā)動機轉(zhuǎn)速下，該上、下邊界距離的幅值范圍差異較大，且隨發(fā)動機轉(zhuǎn)速增大而增大。

3)在地面及空中發(fā)動機穩(wěn)定工作時，提出的喘振檢測量消除了發(fā)動機狀態(tài)變化帶來的影響，在穩(wěn)態(tài)及瞬態(tài)過程中變化范圍差異很小，即喘振檢測量不受發(fā)動機工作狀態(tài)變化的影響。

4)根據(jù)喘振檢測量能夠設(shè)置適應(yīng)發(fā)動機任意狀態(tài)的固定喘振檢測閾值，通過3次喘振試飛數(shù)據(jù)的驗證，其識別發(fā)動機喘振及退出喘振所需時間均小于20 ms，表明喘振檢測方法具有識別率高、報警遲滯小等優(yōu)點。