師利中,孫天宇,蔡舒妤
(1.中國民航大學航空工程學院,天津300300;2.航天工程大學宇航科學與技術系,北京101416)
隨著航空發動機智能維修技術的發展,發動機性能監控對于跟蹤其性能變化、分析其運行狀態、實時預警性能異常具有重要意義。快速存取記錄器(Quick Access Recorder, QAR)是飛機上主要的性能數據記錄設備,其采集性能數據量大、特征維數高,使得數據分析處理困難,目前對QAR數據的應用還不到其潛在價值的10%。
近年來,針對發動機性能數據處理方法的理論和應用研究不斷深入,成果顯著。Zhang等[1]提出關聯事件QAR超限綜合嚴重程度的GREEN模型,將QAR超限與其并發相關事件相關聯;Wang等[2]提出利用QAR數據對典型機載系統進行故障診斷和狀態評估,挖掘相應的QAR參數與失效模式之間的關系;楊慧等[3]針對QAR數據集提出1種改進的小波聚類算法CWave Cluster,劃分非均勻網格,進一步細化邊界網格,對不滿足密度閾值的網格進行處理,最終形成聚類;張春曉等[4]利用HOLT雙參數指數平滑方法,建立基于機載QAR數據的對稱發動機性能參數的差異監控模型,并給出實現對稱發動機差異監控的算法;張鵬等[5]提出1種融合卷積神經網絡與長短時記憶網絡的雙通道融合模型CNN-LSTM,使模型能同時表達數據在空間維度和時間維度上的特征。在應用研究方面,Li等[6]將性能監控數據應用于波音飛機空調系統的溫度狀態預測中;Liang等[7]以民機引氣系統為對象,提出引氣系統故障檢測算法,設計實現了基于QAR數據的民機引氣系統健康監測方法;Chen等[8]提出基于QAR數據的飛機縱向飛行控制律辨識機制;LYU等[9]提出1種基于大規模QAR數據的超限風險評估方法;Wang等[10]提出基于QAR數據的航空發動機潤滑油消耗率在線監測方法;Gao等[11]提出基于QAR數據的民機發動機起動系統健康監測方法;Sun等[12]提出1種基于時態地理信息系統的FQA方法,對不同時空飛行數據的有效相關性進行分析;孫見忠等[13]提出1種基于渦輪葉片外場故障數據及性能數據的渦輪葉片剩余壽命評估方法,為民航發動機在翼壽命評估及送修方案的制定提供決策支持;谷潤平等[14]引入外界環境因素變量,分析發動機燃油流量的影響因素,建立發動機性能異常檢測模型。
本文通過引入變差函數研究發動機性能數據的圖像轉化方法,以融合不同時刻不同參數的性能數據,提升特征表示的關聯性和運算效率,并提出了基于變差函數的發動機性能數據異常判別方法,實現發動機性能圖像中異常的直觀顯示和高效判別。
為實現發動機性能數據的降維,減小多參數關聯性分析的復雜度,需要對性能圖像轉化方法進行研究。
發動機性能數據采集于不同的飛行時刻,其高度、氣壓、溫度等環境因素存在較大差異,使得性能數據不具備可比性。因此,需通過標準化方法將性能數據換算到標準大氣狀態下。
假設性能數據中某性能參數“時間-參數”數據分布圖像如圖1所示。

圖1 “時間-參數”數據分布
可依據對應的高度值對性能參數的分布進行統計,并對其進行曲線擬合,如圖2所示。

圖2 “高度-參數”數據統計分布
對于性能參數α,假設曲線L為高度-性能參數α擬合曲線,曲線方程為y=f(x)。H0為設定的標準高度,則a0=f(H0)。當α=a1時,對應的高度值為H1,將L平移到a1點的高度,則有曲線L'的方程為

式中:y為曲線方程因變量;x為曲線方程自變量;f(x)為擬合曲線方程;a1為H1高度對應的性能參數取值。
性能參數α的標準化轉換公式為

式中:Y'為標準化轉換公式。
采用性能圖像轉化方法將實現性能參數連續值域空間映射到RGB離散彩色空間,如圖3所示。

圖3 性能數據圖像化映射
假設αmax、αmin為標準化修正后性能參數α的最大值、最小值,當α=a時,則性能圖像轉化公式為

式中:R/G/B為計算所得顏色分量值。
由此形成的性能圖像橫軸為時間域的表示,縱軸為各參數取值,以實現對不同時刻不同參數的性能數據的融合。
變差函數常用于紋理結構分析和紋理信息的相關性統計[15-16]。不僅可以用于描述低頻紋理特征在小尺度上的相關性,還能反映高頻紋理信息在大尺度上的結構性。因此,引入變差函數理論對發動機性能圖像進行分析,并對關鍵特征點進行提取,以表示發動機性能圖像的特征信息,實現發動機性能多參數關聯分析。
假設性能圖像大小為M·N,其中任意像素點P坐標為(x,y),對應圖像特征值為 I(x,y)。 在尺寸為2d+1的滑動窗口中,窗口中心點設為(x0,y0),變差函數定義為

式中:h為滑動窗口內兩點的距離;N(h)為滑動窗口內所有距離為h的點對個數;I(P)為滑動窗口內任意像素點的圖像特征值;I(P+h)為滑動窗口內與點P距離為h的像素點的圖像特征值。
細化性能圖像在 0°、45°、90°、135°4 個方向上變差函數定義為

以 0°、45°、90°、135°4 個方向上的變差函數均值為窗口中心點變差函數值,即該窗口中心像素點(x0,y0)的性能特征值為

對于性能圖像,計算其性能特征值,并采用SIFT算法提取關鍵特征點。利用SIFT方法的多量性,保證在較少的性能圖像信息中仍能夠提取較為充分的關鍵特征點。
設第i幅圖像Pi尺寸為Mi·Ni,關鍵特征點共計ni個,第i幅圖像第j個關鍵特征點表示為rij(j=1,2,…,ni),其坐標表示為(xij,yij)。
由于不同性能數據樣本數有所差異,生成的性能圖像大小有所不同,不同圖像間提取的關鍵特征點坐標不具備可比性,因此需對關鍵特征點坐標進行處理,轉化為相對坐標(x'ij,y'ij)。

式中:xij、yij和 x'ij、y'ij分別為第 i幅圖像第 j個關鍵特征點的橫、縱坐標值和橫、縱坐標相對坐標值;Mi、Ni分別為第i幅圖像長度和寬度。
對于任意若干發動機性能圖像Pi,選取圖像Pk,滿足nk=min{ni}。
對于性能圖像Pk中的關鍵特征點r'kl,于其他性能圖像Pi中,尋找匹配關鍵特征點r'ij,滿足
發動機性能圖像差異距離定義為

發動機性能圖像差異距離是識別發動機性能圖像差異的標準。
基于變差函數的發動機性能異常判別方法具體流程如圖4所示。

圖4 方法流程
設N表示發動機性能數據樣本數量,αi表示發動機性能數據中第i個性能參數,其中α0為高度參數,αij表示第i個性能參數第j個數據的取值。
(1)按高度值排序。對高度參數α0進行排序,并擴展至其他性能參數αi。
(2)按高度值分組。按高度參數α0將發動機性能數據樣本分為n組,計算每個性能參數分組后取值。表示第i個性能參數的第s組。對應取值為=
(3)性能參數曲線擬合。對每個性能參數αi分組后取值進行曲線擬合,得到α0-αi曲線Li:αi=f(α0)。
(4)性能數據標準化修正。對每個性能參數αi,將擬合曲線方程Li:αi=f(α0)代入式(2),得到性能參數αi數據標準化修正方程,進行性能數據標準化修正。αij表示標準化修正后第i個性能參數第j個數據的取值。
(5)獲取性能數據變化率。對每個性能參數αi,獲取其性能變化率 δij=αij+1-αij。
(6)數據歸一化。對每個性能參數δij,進行數據歸
(7)數據圖像化。對每個性能參數 αi,依據式(3),得到發動機性能圖像。對于性能參數數量較少的情況,形成圖像時可放大映射。
(8)性能圖像特征表示。對性能圖像中所有像素點 P(x,y),依據式(5)~(11)計算得到各像素變差函數值r(x,y)。
(9)性能關鍵特征點提取。采用SIFT算法提取關鍵特征點,獲取關鍵特征點坐標,并進行相對坐標處理。
(10)對所有性能數據進行性能圖像轉化及關鍵特征點提取。
(11)發動機性能圖像差異計算。選取性能圖像中,關鍵特征點數最少的圖像;對于其他性能圖像,根據式(12),計算關鍵特征點差異距離。
(12)不同航段性能圖像差異識別。比較不同性能圖像差異距離值,在發動機不同運行狀態下的差異距離值有明顯變化,以此識別發動機性能差異。
飛機的起飛、進近和著陸階段約占總飛行時間的17%,但其事故發生概率高達78%。為驗證基于圖像化變差函數的發動機性能數據異常判別方法,選取某機型飛機起飛階段性能數據為實例,共計4組數據樣本量見表1。性能數據包含21個性能參數,分別為:2臺發動機的增壓比REP(engine pressure ratio)、低壓轉子轉速N1(low pressure compressor rotor speed)、高壓轉子轉速N2(high pressure compressor rotor speed)、排氣出口溫度TEG(exhaust gas temperature)、燃油流量FF(fuel flow)、油門桿角度ATL(throttle lever angle)、滑油流量Oq(oil quantity)、滑油壓力Op(oil pressure)、低壓轉子振動值N1V(low pressure compressor rotor vibrate)、高壓轉子振動值 N2V(high pressure compressor rotor vibrate)。在4組數據樣本量中,3組為正常起飛狀態數據;1組為異常起飛狀態數據,起飛階段推油門至40%N1時,出現警告終止起飛。

表1 起飛階段性能數據實例樣本
首先依據性能圖像轉化方法對實例數據進行預處理,將高維度性能數據轉化成2維性能圖像,獲得實例數據對應性能圖像,如圖5所示。

圖5 4個實例性能圖像轉化結果
從圖 5(a)~(c)中可見,對于正常起飛階段,經標準化修正后性能參數較為平穩;轉化的性能圖像各參數對應圖像區域,沒有明顯顏色變化;局部存在由飛行高度引起的輕微波動。
從圖5(d)中可見,對于異常中斷起飛狀態下的性能圖像,PER、N1、N2、TEG、FF等 5 個性能參數所對應的圖像區域色彩波動明顯,Oq、N1V2個參數值明顯高于正常狀態。
依據式(9),對發動機性能圖像計算其性能特征值,實現實例圖像進行特征表示;分別采用Harris方法和SIFT方法對特征表示后的性能圖像提取關鍵特征點,結果圖像如圖6、7所示。

圖6 Harris方法提取關鍵特征點

圖7 SlFT方法提取關鍵特征點
從圖6、7中可見,對于信息并不充分的發動機性能圖像,SIFT方法所提取的特征點明顯比Harris方法數量多,為后續性能特征分析提供了數據保證。
由于實例性能數據樣本量不同,需將關鍵特征點坐標轉化為相對坐標,依據式(10)、(11),對 SIFT 方法所提取的關鍵特征點轉化為相對坐標,散點如圖8所示。
根據關鍵特征點匹配并計算差異距離值,結果見表2,該方法運算時間如圖9所示。

圖8 關鍵特征點相對坐標散點

表2 關鍵特征點差異距離值

圖9 實例運算時間
由關鍵特征點相對坐標散點圖和差異距離值可見,實例a、b在關鍵特征點分布上有較大的相似性,差異距離值極為接近,可由此判斷其運行狀態具有較高的一致性。
實例c與a、b在關鍵特征點數量上有所差異,但三者差異距離值較為接近,因此三者應處于相同運行狀態。
實例d與其他實例關鍵特征點分布存在較大差異,差異距離值偏差較大,由此可判斷其運行狀態與其他實例不同。
假設已知其中1個實例的狀態,即可判斷出其他實例的狀態。
從運行時間來看,圖像化、特征表示、特征提取階段所需計算時間相對較長,異常識別階段均在毫秒級時間范圍內完成運算過程,總體運行效率較高。
本文通過研究發動機性能圖像轉化方法,結合變差函數理論,提出了1種基于圖像化變差函數的發動機性能數據異常判別方法,并以發動機實際性能數據為實例,進行了方法驗證。由驗證結果可知:
(1)從運算速度方面,本文方法能夠在毫秒級時間范圍內完成各階段運算過程,時間復雜度較低。
(2)從異常判別質量方面,可依據關鍵特征點分布和差異距離值,實現對不同性能圖像特征的分類,從而實現對應性能數據運行狀態的自動判別。
(3)本文方法為發動機性能數據的智能監控與機器分析提供了數據基礎,后續將對性能特征的完備表示、機器分析進一步研究,以提高發動機性能監控的智能化程度。