董 慶，李本威，錢仁軍，滕懷亮，張 赟

(1.海軍裝備部駐蘇州地區(qū)軍事代表室， 江蘇 蘇州 215000；2.海軍航空大學 航空基礎學院， 山東 煙臺 264000)

航空發(fā)動機是一個復雜的非線性時變系統(tǒng)[1]。在發(fā)動機故障診斷[2]、性能分析與仿真及其控制規(guī)律設計[3]等多個領域，發(fā)動機模型都發(fā)揮著無可代替的作用。針對裝配某型渦軸發(fā)動機的直升機將在不同大氣環(huán)境(高原、高空及高低溫等)下進行飛行試驗的相關工作，為實現(xiàn)對發(fā)動機狀態(tài)的實時監(jiān)控和對發(fā)動機穩(wěn)態(tài)控制規(guī)律的進一步研究，建立精度和實時性均滿足實際要求的發(fā)動機穩(wěn)態(tài)模型成為相關工作關鍵的一環(huán)。

目前，航空發(fā)動機模型建立的方法主要有兩種：一種是分析部件氣動特性、機械關系的解析法，另一種是系統(tǒng)辨識的方法。解析法通過假設、近似處理和大量的迭代運算來獲取發(fā)動機的性能數(shù)據，過程復雜。而系統(tǒng)辨識的方法無需了解發(fā)動機復雜的部件特性，只需足夠的發(fā)動機輸入輸出參數(shù)即可快速建立精度滿足要求的發(fā)動機模型。目前，基于數(shù)據驅動的模型辨識方法得到快速的發(fā)展，常用的辨識方法有最小二乘法[4]和子空間狀態(tài)辨識法[5]等。它們假設發(fā)動機是線性時不變系統(tǒng)，導致建立的發(fā)動機動態(tài)模型精度不能達到很好的效果。近幾年，人工神經網絡[6]。(artificial neural networks，ANN)、支持向量機[7](support vector machine，SVM)和極限學習機[8]。(extreme learning machine，ELM)等機器學習算法快速發(fā)展，具有了很好的非線性學習能力，很多人也將它們應用到航空發(fā)動機模型的建立。陳超等[9]利用RBF神經網絡對航空發(fā)動機起動模型進行辨識與仿真；王海濤等[10]利用稀疏最小支持向量機建立發(fā)動機動態(tài)過程模型。這些傳統(tǒng)的神經網絡存在初始參數(shù)選擇無依據、無反饋單元等各種問題。耿宏等[11]利用NARX網絡建立發(fā)動機參數(shù)動態(tài)辨識模型；Pogorelv等[12]利用動態(tài)神經網絡對雙軸燃氣輪機差分形式的非線性動態(tài)模型進行辨識，并應用于啟動控制。NARX網絡具有反饋延時單元，具有很好的動態(tài)特性，且對非線性動態(tài)系統(tǒng)具有良好的非線性逼近能力，效果明顯優(yōu)于前饋神經網絡，但是網絡初始特性參數(shù)由隨機函數(shù)產生影響辨識模型精度。

粒子群優(yōu)化算法(particle swarm optimization,PSO)穩(wěn)定性好且具有較好的全局搜索和尋優(yōu)能力，利用PSO對NARX網絡的初始參數(shù)進行尋優(yōu)，使得NARX網絡取得更好的效果。文獻[13]驗證了粒子群優(yōu)化算法的優(yōu)越性。文獻[14]驗證了NARX網絡對非線性動態(tài)系統(tǒng)具有良好的非線性逼近能力。

針對以上分析，本文以直升機飛參數(shù)據為驅動，結合提出的三點自適應判斷法從飛參數(shù)據中提取發(fā)動機穩(wěn)態(tài)數(shù)據樣本，采用PSO-NARX方法對渦軸發(fā)動機穩(wěn)態(tài)模型進行辨識。利用樣本之外架次的飛參數(shù)據對辨識模型的性能進行驗證并與其他不同辨識算法相比較。結果表明本文PSO-NARX方法對渦軸發(fā)動機穩(wěn)態(tài)模型辨識的可行性與優(yōu)越性，為后續(xù)在高原、高空及高低溫環(huán)境條件下進行飛行試驗的此型渦軸發(fā)動機狀態(tài)實時監(jiān)控和控制規(guī)律的優(yōu)化奠定基礎。同時，也為工廠對此型渦軸發(fā)動機即將開展的試車試驗提供相關支持。

1 基于PSO的NARX網絡原理

1.1 NARX網絡模型

英文縮寫詞NARX(nonlinear auto regressive with exogenous inputs model)的意思為有外部輸入的非線性自回歸模型，是線性外生輸入的自回歸(ARX)模型改進后在非線性情況下的應用。其通過添加延時單元使過去某些時刻的輸入、輸出成為模型的輸入，再利用多層感知器的非線性映射能力，來實現(xiàn)時間序列的動態(tài)預測[11]。NARX模型可由帶有外部輸入的非線性差分方程表示為

y(t+1)=f[y(t),y(t-1),…,y(t-k),

x(t),x(t-1),…,x(t-m)]

(1)

式中：f表示自變量的一個非線性函數(shù) ;y(t-k)表示輸出參數(shù)相應時延為k的數(shù)值；x(t-m)表示輸入參數(shù)相應時延為m的數(shù)值；y(t+1)表示得到的下一時刻的輸出數(shù)值。NARX網絡模型見圖1。

圖1 NARX網絡模型示意圖

某t時刻其第i個隱含層節(jié)點的輸出di(t)為

(2)

式中：fsig為隱層節(jié)點使用的雙曲線正切sigmoid函數(shù)；wir(t)、wil(t)為t時刻第i個隱層節(jié)點與輸入u(t-r)、輸出y(t-l)之間的權值；bi為第i個隱層節(jié)點的偏置。第j個輸出層節(jié)點輸出y(t)為

(3)

式中：wij(t)為t時刻第i個隱層節(jié)點到第j個輸出層節(jié)點的權值；Tj為第j個輸出層節(jié)點偏置；N為隱含層節(jié)點個數(shù)。

1.2 粒子群優(yōu)化算法(Particle Swarm Optimization)

粒子群優(yōu)化算法屬于進化算法的一種，源于對鳥群捕食的行為研究。該算法從隨機解出發(fā)，通過迭代找到最優(yōu)解，通過適應度來評價解的品質，并通過追隨當前搜索到的最優(yōu)值來尋找全局最優(yōu)[16-17]。

每個粒子都有位置xi和速度vi兩個屬性，在每次的迭代中，粒子通過跟蹤個體極值和全局極值來更新自己。假設在一個D維的目標搜索空間中，即每個粒子(解)都是一個D維的向量，粒子群由N個粒子構成，則其中第i個粒子位置可以表示為

xi=(xi1,xi2,…,xiD)，i=1,2,…,N

第i個粒子的速度也是一個D維的向量，記為

vi=(vi1,vi2,…,viD)，i=1,2,…,N

并設定第i個粒子迄今為止搜索到的最優(yōu)位置為個體極值，記為pi，而整個粒子群迄今為止搜索到的最優(yōu)位置為全局極值，記為pg。同樣，pi、pg均為D維向量。在找到這兩個最優(yōu)值后，粒子通過下面的公式(4)來更新自己的速度和位置。

(4)

式中：c1、c2為慣性因子和學習因子，也稱加速常數(shù)，反映粒子間信息交換的強度，通常取c1=c2=2，r1、r2為[0,1]范圍內的均勻隨機數(shù)，pi為第i個粒子迄今為止搜索到的最優(yōu)位置，即個體極值，pg為整個粒子群迄今為止搜索到的最優(yōu)位置，即全局極值。

目前多采用的使線性遞減權值策略，動態(tài)w能夠獲得比固定值更好的結果。w的引入使得PSO算法性能有了很大的提高，提高了算法的全局和局部的搜索能力，w的值為

w=(wini-wend)(Gk-g)/Gk-wend

(5)

式中：Gk為最大迭代次數(shù)；wini為初始慣性權值；wend為迭代至最大進化代數(shù)時的慣性權值，典型權值wini取0.9，wend取0.4[17]。

1.3 基于PSO的NARX網絡特征參數(shù)優(yōu)化流程

NARX網絡結構將隱含層個數(shù)設定為1，隱含層神經元個數(shù)q值的經驗公式為

(6)

式中：d為輸入層節(jié)點數(shù)；l為輸出層節(jié)點數(shù)；a為10～14的常數(shù)。

利用PSO優(yōu)化NARX網絡中輸入層到隱含層的權值vih,隱含層到輸出層的權值whj,隱含層神經元的偏置bh，輸出層神經元的偏置bj，i=1,2,…,d；h=1,2,…,q；j=1,2,…,l。以下是具體的優(yōu)化流程：

步驟1：初始化。初始化粒子群群體數(shù)量N，種群數(shù)量和最大迭代步數(shù)。每個粒子mi由vih、whj、bh、bj組成：

mi=(v1h,v2h,…,vih,w1j,w2j,…,whj,b1,b2,…,bj,b1,b2,…,bh)

利用rands函數(shù)初始化每個粒子mi的xi和vi，xi和vi中的每個元素的取值范圍都在[0,1]之間。

步驟2：計算每個粒子適應度的值。將模型預測輸出與目標輸出的均方誤差作為粒子的適應度函數(shù)值，如式(7)所示，均方誤差用來評價模型與實測數(shù)據的吻合程度，PSO算法尋優(yōu)的目標就是使粒子的適應度函數(shù)值即均方誤差最小。

(7)

式中：Emse(y,yd)為均方誤差；y(t)為模型的預測值；yd(t)為輸出實測值。

步驟3：計算每個粒子，用它的適應度值ffitness(mi)和個體極值pi比較，如果ffitness(mi)

步驟4：重復步驟1～3直到最大迭代步數(shù)，利用PSO算法得到NARX網絡的特征參數(shù)。

2 渦軸發(fā)動機穩(wěn)態(tài)模型辨識方法

2.1 建立渦軸發(fā)動機穩(wěn)態(tài)模型

航空發(fā)動機是一個高度非線性的時變系統(tǒng)。大氣靜溫T0、絕對氣壓高度Hp、飛行馬赫數(shù)Ma、燃油流量Wf都會影響發(fā)動機工作過程。渦軸發(fā)動機主要工作狀態(tài)參數(shù)為燃氣發(fā)生器轉速ng、發(fā)動機排氣溫度EGT和發(fā)動機扭矩T。

發(fā)動機當前時刻的工作狀態(tài)與過去時刻的工作狀態(tài)密切相關。渦軸發(fā)動機穩(wěn)態(tài)模型的輸入參數(shù)設置為T0、Hp、Ma、Wf的當前時刻值與過去時刻的值以及ng、EGT和T過去時刻的值，輸出參數(shù)設置為ng、EGT和T的當前時刻值。渦軸發(fā)動機穩(wěn)態(tài)模型輸入輸出參數(shù)如表1所示。

綜上，渦軸發(fā)動機穩(wěn)態(tài)模型為：

(8)

式中：fss為需要辨識的發(fā)動機模型；k、m為輸入輸出延時階次。文獻[18]提出：通常，發(fā)動機模型輸入輸出參數(shù)延時階次設置為2為最佳方案。因此，將k和m設置為2。

表1 發(fā)動機穩(wěn)態(tài)模型輸入輸出參數(shù)

2.2 基于PSO-NARX網絡的穩(wěn)態(tài)模型辨識方法

根據表(1)知輸入層節(jié)點數(shù)d為4，輸出層節(jié)點數(shù)l為3，根據式(6)得到本文NARX網絡隱含層神經元個數(shù)q選為15。

基于PSO-NARX網絡的渦軸發(fā)動機穩(wěn)態(tài)模型辨識的具體流程如圖2所示。

步驟1：根據發(fā)動機穩(wěn)態(tài)工作原理，建立渦軸發(fā)動機的穩(wěn)態(tài)模型如式(8)，包括輸入、輸出參數(shù)和函數(shù)關系。

步驟2：利用三點自適應法從大量飛參數(shù)據中提取發(fā)動機穩(wěn)態(tài)時的數(shù)據，建立辨識所需的樣本，并對樣本進行預處理。

步驟3：設置PSO算法的種群粒子數(shù)和最大迭代次數(shù)等基本參數(shù)，利用PSO對NARX網絡的特征參數(shù)進行尋優(yōu)。

步驟4：以提取處理的辨識樣本為驅動，利用優(yōu)選特征參數(shù)后的NARX網絡對所建立的渦軸發(fā)動機穩(wěn)態(tài)模型進行辨識。

步驟5：利用驗證樣本對辨識得到的穩(wěn)態(tài)模型進行驗證。

綜上，步驟1～5構成了基于PSO-NARX網絡的渦軸發(fā)動機穩(wěn)態(tài)模型的辨識方法。

圖2 渦軸發(fā)動機穩(wěn)態(tài)模型辨識流程框圖

3 辨識過程與結果分析

以直升機飛參數(shù)據為驅動，建立基于PSO-NARX網絡的渦軸發(fā)動機穩(wěn)態(tài)辨識模型，對辨識模型的性能進行驗證并與NARX、PSO-BP、PSO-SVM建立的辨識模型進行對比。

3.1 從飛參數(shù)據中提取辨識樣本與數(shù)據預處理

針對從飛參數(shù)據中對發(fā)動機穩(wěn)態(tài)狀態(tài)進行識別的問題。本文在此提出一種穩(wěn)態(tài)狀態(tài)的識別算法——三點自適應判斷法。設換算轉速-時間序列為向量ng,c(t)，穩(wěn)態(tài)轉速限制值采用動態(tài)變化值δn為

δn=0.01ng,c(t)

(9)

采用三點自適應判斷法判斷第i個點是否為穩(wěn)態(tài)點，將三點數(shù)值相減得到：

(10)

若Δ<δn，則第i個轉速點為穩(wěn)態(tài)工作點，并令ng,c(i+1)=ng,c(i)，將時間平移繼續(xù)進行判斷。對符合轉速波動的穩(wěn)態(tài)點進行記錄后通過穩(wěn)態(tài)的時間限制值δtime,確定發(fā)動機的穩(wěn)定工作狀態(tài)。根據渦軸發(fā)動機大量分析數(shù)據、工作經驗及其工作特點，通常選取轉速限制值為1%，穩(wěn)態(tài)時間限制值為10 s。

使用DBSCAN聚類方法對某一架次飛參數(shù)據(飛參數(shù)據的采樣頻率為1s)中1發(fā)位置發(fā)動機轉速去噪處理后[20]，利用三點自適應判斷法得到此架次穩(wěn)定狀態(tài)處理曲線如圖3所示，從中可以得出發(fā)動機的穩(wěn)定狀態(tài)。圖3(a)是利用三點自適應判斷法得出的穩(wěn)定狀態(tài)判斷曲線，圖3(b)是圖(a)中3 000～5 000 s段的局部放大。例如，圖3(b)圖方框內的ng,c=28 997 r/min(即水平線所對應的換算轉速)時的穩(wěn)定狀態(tài)即為發(fā)動機的穩(wěn)態(tài)點，由此得到此架次此渦軸發(fā)動機在裝機過程中轉速從0～32 400 r/min(最大連續(xù)狀態(tài))之間所有穩(wěn)態(tài)點。

圖3 飛參數(shù)據穩(wěn)態(tài)處理曲線

使用同樣的數(shù)據處理方法對5臺直升機2016—2018年的27個架次飛參數(shù)據中1發(fā)位置發(fā)動機轉速進行處理，得到同型號渦軸發(fā)動機轉速從0到最大連續(xù)轉速之間3 564個穩(wěn)態(tài)點，每個穩(wěn)態(tài)點的數(shù)據記為一組樣本數(shù)據，3 564組辨識樣本共129 076個數(shù)據點作為渦軸發(fā)動機穩(wěn)態(tài)模型辨識樣本。此27個架次飛參數(shù)據來源于不同工況下的直升機飛行數(shù)據，包括低溫、高溫等多種飛行工況，發(fā)動機的工作點幾乎涵蓋了整個飛行包線。

根據構建的渦軸發(fā)動機穩(wěn)態(tài)模型如式(7)所示，確定樣本參數(shù)為：T0、Hp、Ma、Wf、ng、EGT和T。對選取的樣本數(shù)據進行濾波和歸一化的處理。將處理后的數(shù)據作為模型的辨識樣本，在模型的訓練與驗證結束后對數(shù)據進行還原。采用辨識模型輸出值與飛參數(shù)據實測值最大相對誤差以及最大相對誤差均值來評價建立的辨識模型性能。相對誤差Eref為

(11)

式中：y(t)為辨識模型的輸出值；yd(t)為飛參數(shù)據實測值。

3.2 模型辨識結果

綜合考慮計算時間和收斂性以及資源占用的最優(yōu)值等因素，將PSO算法的最大迭代步數(shù)設置為50，種群粒子數(shù)量設置為30[21-23]。利用PSO算法對NARX網絡特征參數(shù)尋優(yōu)，尋優(yōu)過程的適應度函數(shù)值曲線如圖4所示。結合從飛參數(shù)據提取的辨識樣本，利用PSO算法優(yōu)化過特征參數(shù)的NARX網絡對建立的渦軸發(fā)動機模型進行回歸辨識。

圖4 尋優(yōu)過程的最優(yōu)個體適應度函數(shù)值曲線

利用樣本外的4個架次不同編號發(fā)動機穩(wěn)態(tài)數(shù)據對得到的辨識模型非樣本點推廣能力進行驗證，分別取轉速ng,c=25 417、28 703、28 997、30 301 r/min時的4個穩(wěn)態(tài)點數(shù)據為4組驗證樣本，分別記為M1、M2、M3和M4。通過將驗證樣本輸入辨識模型得到的輸出值與飛參數(shù)據實測值對比來驗證辨識模型的性能。例如，以M3組為驗證樣本時，取M3組中200個連續(xù)的數(shù)據點，得到辨識模型輸出參數(shù)ng、EGT和T與飛參數(shù)據實測值曲線和相對誤差分布曲線如圖5、圖6、圖7所示。

圖5 M3為驗證樣本ng辨識結果曲線

圖6 M3為驗證樣本EGT辨識結果曲線

4組數(shù)據分別作為驗證樣本時得到ng、EGT和T的最大相對誤差及最大相對誤差均值如圖8所示。從圖8可以得到，不同驗證樣本時，輸出參數(shù)ng、EGT和T的最大相對誤差分別為0.11%、0.26%和0.95%，最大相對誤差均值分別為0.105%、0.193%和0.865%。從圖5、圖6、圖7、圖8可以得到，辨識模型輸出結果很好地逼近了飛參數(shù)據實測值。說明本文提出的基于PSO-NARX網絡的發(fā)動機穩(wěn)態(tài)模型精度滿足實際應用時的精度要求。

圖8 M1 M2 M3 M4為驗證樣本時輸出參數(shù)ng、EGT、T誤差的直方圖

3.3 不同方法辨識效果對比

為了更好的說明本文提出的基于PSO-NARX的渦軸發(fā)動機穩(wěn)態(tài)模型辨識的優(yōu)越性，用相同的訓練樣本和驗證樣本，利用PSO-BP、PSO-SVM和NARX網絡辨識建立的渦軸發(fā)動機穩(wěn)態(tài)模型。

BP神經網絡使用MATLAB R2014a中的工具箱，設置為輸出層、單隱含層和輸出層，隱含層神經元個數(shù)依據經驗公式設置為15，各層連接權值及神經元偏置通過PSO尋優(yōu)得到，選用Levenberg-Marquardt算法。SVM對渦軸發(fā)動機穩(wěn)態(tài)模型進行回歸辨識時，利用PSO對SVM特征參數(shù)平衡因子C、不敏感度ε和徑向基核函數(shù)參數(shù)σ進行尋優(yōu)。兩種算法中PSO的最大迭代步數(shù)均設置為50，種群粒子數(shù)量均設置為30。NARX網絡使用MATLAB R2014a中的工具箱，隱含層神經元個數(shù)依據經驗公式(6)設置為15，輸入輸出延時階次均設為2，權值及神經元偏置由rands函數(shù)隨機給定，采用train函數(shù)對NARX網絡進行訓練。

圖9 4種方法誤差曲線

圖9列出了采用4種不同方法辨識模型得到的4組驗證樣本輸出參數(shù)ng、T4和T的最大相對誤差均值。表2為4種不同方法建立渦軸發(fā)動機穩(wěn)態(tài)辨識模型精度，由表可見輸出參數(shù)ng、EGT和T的最大相對誤差與最大相對誤差均值。

表2 4種不同方法辨識模型精度

從圖9和表2可知，PSO-NARX的不同驗證樣本最大相對誤差及最大相對誤差均值均小于PSO-BP、NARX和PSO-SVM。主要是因為PSO-NARX具有延時單元，將過去時刻的值反饋對當前時刻的輸出產生影響，PSO對NARX網絡的特征參數(shù)進行優(yōu)化使得辨識模型的精度更高。BP神經網絡是前饋神經網絡無反饋單元，NARX網絡具有反饋單元，但初始權值、偏置也是由隨機函數(shù)產生，網絡特征參數(shù)沒有進行尋優(yōu)。PSO-SVM精度也較好，但SVM對大規(guī)模訓練樣本會耗費大量的計算機內存和運算時間。由于PSO-NARX方法需要在最大迭代步數(shù)內對全部粒子進行適應度更新與尋優(yōu)，PSO-NARX方法的訓練時間也較長，但平均訓練時間明顯短于PSO-SVM方法。

4 結論

1) 辨識模型輸出參數(shù)ng、EGT和T的辨識結果很好的逼近了飛參數(shù)據實測值，證明了基于PSO-NARX網絡的渦軸發(fā)動機穩(wěn)態(tài)模型辨識方法的可行性。

2) 采用4組驗證樣本對辨識模型的性能進行驗證，得到輸出參數(shù)ng、EGT和T的最大相對誤差分別為0.11%、0.26%和0.95%，最大相對誤差均值分別為0.105%、0.193%和0.865%。表明了基于PSO-NARX網絡的渦軸發(fā)動機穩(wěn)態(tài)模型的精度達滿足實際應用的精度要求。

3) 在相同的訓練樣本和驗證樣本情況下，本文采用的基于PSO-NARX網絡的渦軸發(fā)動機穩(wěn)態(tài)模型精度優(yōu)于PSO-BP、NARX和PSO-SVM方法辨識的渦軸發(fā)動機穩(wěn)態(tài)模型，有效的解決了渦軸發(fā)動機穩(wěn)態(tài)模型建模復雜且精度不高的問題。

4) 若已知某些大氣環(huán)境條件和該條件下燃燒室供油特性，即可通過辨識模型獲取此大氣環(huán)境下的發(fā)動機穩(wěn)態(tài)性能，表明了基于PSO-NARX網絡的渦軸發(fā)動機穩(wěn)態(tài)辨識模型可用于其他任意大氣環(huán)境下發(fā)動機穩(wěn)態(tài)性能的遞推估算。

5) 此辨識模型可用于后續(xù)在高原、高空及高低溫環(huán)境條件下進行飛行試驗的發(fā)動機狀態(tài)的實時監(jiān)控，為渦軸發(fā)動機控制規(guī)律的優(yōu)化奠定基礎，也為工廠對此型渦軸發(fā)動機試車試驗提供相關支持。