基于RELAX 算法的飛機(jī)油耗性能估計(jì)方法

劉家學(xué)， 馬 濤， 陳靜杰

(中國(guó)民航大學(xué)航空自動(dòng)化學(xué)院，天津 300300)

0 引言

2013 年，我國(guó)民航業(yè)運(yùn)輸飛機(jī)起降731．5 萬(wàn)架次，運(yùn)輸周轉(zhuǎn)量達(dá)671．7 億噸千米，航空煤油消耗近2000萬(wàn)噸，碳排放量近6200 萬(wàn)噸［1］，運(yùn)輸飛機(jī)燃油成本目前已經(jīng)占到航空公司運(yùn)行成本的60%以上［2］。對(duì)航空運(yùn)輸業(yè)而言，良好的油耗性能同時(shí)也意味著較低的運(yùn)營(yíng)成本和碳排放量，航空公司在運(yùn)力調(diào)配上更為關(guān)注機(jī)型的油耗性能指標(biāo)評(píng)估。因此，利用實(shí)際飛行數(shù)據(jù)建立機(jī)型油耗性能模型并提取指標(biāo)參數(shù)越來(lái)越引起業(yè)界的關(guān)注和重視。目前關(guān)于飛機(jī)油耗性能的模型中，文獻(xiàn)［3］使用了神經(jīng)網(wǎng)絡(luò)的方法從飛機(jī)水平面運(yùn)動(dòng)、加速度、風(fēng)向風(fēng)速、傾斜角、縱向加速度等方面估計(jì)計(jì)算飛機(jī)瞬時(shí)燃油流量，偏重微觀分析;文獻(xiàn)［4］側(cè)重于從效率角度評(píng)估航線運(yùn)輸，提出了優(yōu)化飛機(jī)運(yùn)行航路、改善機(jī)務(wù)維修方法、提高飛行操作能力等方法來(lái)提高飛機(jī)燃油效率，偏重于管理科學(xué)領(lǐng)域;文獻(xiàn)［5］綜述了飛機(jī)燃油消耗中可用的幾類機(jī)器學(xué)習(xí)模型，包括使用線性回歸、模型樹(shù)、MLP 多層感知、SVM 支持向量機(jī)、自組織映射與Hopfield 神經(jīng)網(wǎng)絡(luò)等，并分析了幾類方法各自的特點(diǎn)。關(guān)于飛機(jī)飛行整體油耗性能評(píng)估領(lǐng)域研究較少，而按照飛行小時(shí)與飛行千米評(píng)估飛機(jī)油耗性能應(yīng)區(qū)分多個(gè)運(yùn)行模式，不同模式表現(xiàn)不同的指標(biāo)模型，區(qū)分模型需要一種合適的分離理論與估計(jì)方法。

松弛譜估計(jì)(RELAX)算法是一種信號(hào)分離算法，通過(guò)對(duì)信號(hào)分量的反復(fù)迭代估計(jì)信號(hào)參數(shù)。算法在變換域中運(yùn)算，不改變信號(hào)的相對(duì)強(qiáng)度，信號(hào)分量在整體范圍內(nèi)表現(xiàn)相似特性，容易分離信號(hào)，但無(wú)法很好地估計(jì)不同時(shí)間呈現(xiàn)不同模式的混合信號(hào)［5］，基于廣義內(nèi)積(GIP)的非一致性檢測(cè)(NHD)方法能夠檢測(cè)出信號(hào)分量的界限。針對(duì)不同時(shí)間區(qū)間上不同信號(hào)強(qiáng)度作用這個(gè)特性進(jìn)行模式區(qū)分。

1 飛機(jī)油耗性能指標(biāo)模型

民航飛機(jī)飛行主要分為起飛爬升段，空中巡航段與下降進(jìn)近段3 個(gè)階段，各段的油耗性能指標(biāo)相差迥異，如圖1 所示。

圖1 Trent700 發(fā)動(dòng)機(jī)典型燃油流量曲線Fig．1 Typical fuel flow chart of Trent700 aero engine

起飛爬升階段總體燃油流量平均處于高值震蕩區(qū)間。3000 s 后進(jìn)入巡航階段，燃油流量趨于平穩(wěn)，均值相比爬升段減小明顯;在保持巡航階段基本平穩(wěn)燃油流量后飛機(jī)進(jìn)入下降階段，下降段燃油流量呈現(xiàn)低值震蕩。3 段在燃油流量分布上十分明顯。設(shè)起飛爬升段、巡航飛行段、下降進(jìn)近段飛行油耗指標(biāo)樣本集為{ycl，ycr，yde}，ycl＞ycr＞yde，若按照飛行距離分類有如下的兩個(gè)特性。

1)飛機(jī)在本場(chǎng)起降訓(xùn)練或近場(chǎng)調(diào)機(jī)等短程航線飛行時(shí)，多處于爬升-下降或爬升-低高度巡航-下降兩種模式。兩種模式中起飛爬升段距離占飛行全程的絕大部分，表現(xiàn)在油耗指標(biāo)上則ycl占據(jù)全程油耗指標(biāo)的主導(dǎo)地位，小部分的ycr與yde無(wú)法抵消影響，總體指標(biāo)很高。若根據(jù)能量守恒，短程航線中，爬升段的主導(dǎo)使飛機(jī)動(dòng)能源源不斷轉(zhuǎn)化為勢(shì)能，同時(shí)還要保證一定的分型速度產(chǎn)生足夠升力。如圖2 所示，地速與高度的快速增加使對(duì)應(yīng)動(dòng)能與勢(shì)能都在提升，表現(xiàn)在油耗指標(biāo)水平平均較高。

2)隨著飛行距離的增加，飛機(jī)巡航階段所占比重逐步加大。表現(xiàn)在油耗上ycr在3 段指標(biāo)中所占比例也逐步加大，成為全程航段的主導(dǎo)油耗指標(biāo)。如圖2所示，在巡航期間，飛機(jī)動(dòng)能與勢(shì)能基本保持不變，不產(chǎn)生額外的轉(zhuǎn)化關(guān)系，油耗指標(biāo)ycr穩(wěn)定偏小。此時(shí)可認(rèn)為隨著巡航距離增加時(shí)，占據(jù)主導(dǎo)的ycr不斷中和ycl與yde的高低指標(biāo)影響，使整體指標(biāo)回歸平穩(wěn)，回歸以指數(shù)形式漸漸逼近飛機(jī)最佳油耗指標(biāo)，并且可以以距離為界分類油耗性能模式。

圖2 某航段飛行高度與地速歸一化曲線Fig．2 Normalized flight altitude and ground speed diagram

在航空公司指標(biāo)評(píng)價(jià)體系中，飛機(jī)千米油耗yfkm多用來(lái)測(cè)算某機(jī)型與機(jī)隊(duì)本身燃油利用效率，本文使用此指標(biāo)評(píng)價(jià)上述指標(biāo)模型較為準(zhǔn)確。圖3 所示為某飛機(jī)千米油耗與飛行距離的樣本采樣圖，可看到隨著距離增加，千米油耗yfkm遞減，遞減趨勢(shì)明顯呈現(xiàn)不同模式。指數(shù)下降型函數(shù)可以有效地描述這一遞減過(guò)程，并且針對(duì)不同下降模式有著良好的參數(shù)描述。

模型為

式中:Ki表示指數(shù)遞減項(xiàng)系數(shù);n 代表主衰減分量個(gè)數(shù)，根據(jù)上述分析，n=2，分類兩種模式;λi表示趨勢(shì)變化實(shí)指數(shù)，根據(jù)回歸趨勢(shì)從大到小緩慢逼近且兩種模式相差很大的特點(diǎn)，認(rèn)為λ2＜λ1＜0;C 為常數(shù)項(xiàng)，代表當(dāng)飛行距離x 趨近于無(wú)窮時(shí)油耗性能回歸的極限值;en代表噪聲信號(hào)即在其他因素作用下的油耗性能參數(shù)的噪聲起伏。通過(guò)不同的實(shí)指數(shù)與衰減系數(shù)搭配固定油耗常數(shù)可以估計(jì)某機(jī)型的油耗指標(biāo)。

圖3 A330 飛機(jī)飛行距離與千米油耗數(shù)據(jù)Fig．3 A330 flight distance vs fuel mileage

2 油耗模型的參數(shù)估計(jì)

根據(jù)收集到的油耗樣本參數(shù)相應(yīng)的飛行距離數(shù)據(jù)，數(shù)據(jù)分布與航線航路距離吻合，采樣間隔并不平均，在某些區(qū)間范圍內(nèi)采樣點(diǎn)比較分散，需要對(duì)傳統(tǒng)的RELAX 算法進(jìn)行一定的算法適應(yīng)改變。因?yàn)槟Ｐ蜔o(wú)法采用FFT 方法進(jìn)行運(yùn)算，本文采用RELAX 算法的信號(hào)分離與參數(shù)迭代的思想進(jìn)行運(yùn)算。

根據(jù)距離油耗性能指標(biāo)模型

取

定義y1(x)為弱目標(biāo)信號(hào)，y2(x)為強(qiáng)目標(biāo)信號(hào)，不妨假設(shè)當(dāng)距離x 增大到3 倍的距離常數(shù)，即x ＞3λ-1i 時(shí)，可認(rèn)為信號(hào)已衰減完畢。估計(jì)時(shí)確定初始常數(shù)信號(hào)值，從弱目標(biāo)信號(hào)到強(qiáng)目標(biāo)信號(hào)，利用差分得到修正常數(shù)信號(hào)，整體不斷交叉迭代修正所有分量信號(hào)，直到達(dá)到收斂條件為止。具體算法過(guò)程如下所述。

1)設(shè)定迭代次數(shù)，迭代次數(shù)應(yīng)隨著信號(hào)信噪比增加而減小，為了盡量達(dá)到良好性能，可取迭代次數(shù)n =100。并初始化常數(shù)信號(hào)，因?yàn)榭赡艹霈F(xiàn)的噪聲干擾信號(hào)，無(wú)法準(zhǔn)確得知當(dāng)飛行距離x 趨向于無(wú)窮時(shí)油耗性能常數(shù)C0的值，故取信號(hào)后部平均，根據(jù)噪聲大小可動(dòng)態(tài)選取平均值范圍，一般認(rèn)為后5%區(qū)間項(xiàng)均值作為初始C^0較為合適。設(shè)z1i為去除初始常數(shù)分量后的信號(hào)，對(duì)于每點(diǎn)信號(hào)yi有

2)進(jìn)行弱分量與強(qiáng)分量的分層，分層原理借鑒了空時(shí)信號(hào)自適應(yīng)處理中常用的非均勻檢測(cè)(NHD)中廣義內(nèi)積(GIP)的方法。GIP 可以計(jì)算信號(hào)的馬氏距離，分類不同信號(hào)的作用區(qū)間［6］。對(duì)z1i信號(hào)做相鄰點(diǎn)差分，形成新的序列s，并且將新序列進(jìn)行長(zhǎng)度為K的多點(diǎn)分層形成序列R^。通過(guò)GIP 后的數(shù)據(jù)GIP(x)可明顯看出信號(hào)中兩種樣本集的相似度差異。計(jì)算GIP 的變異系數(shù)，根據(jù)變異系數(shù)大小取平均值與系數(shù)標(biāo)準(zhǔn)差組合進(jìn)行樣本選取判斷。

在樣本選取界定范圍時(shí)需要考慮噪聲影響，在有噪聲信號(hào)附加下GIP 值會(huì)在同頻率處出現(xiàn)抖動(dòng)，影響樣本容量的選取，此時(shí)，可考慮使用二次GIP 運(yùn)算或GIP 性能改進(jìn)算法等類似濾波算法［7］，再一次細(xì)化可用信號(hào)區(qū)間，在確定后的區(qū)間內(nèi)進(jìn)行參數(shù)估計(jì)［8-9］。

樣本的NHD 分類如圖4 所示。

圖4 樣本的NHD 分類Fig．4 NHD classification of the sample

3)先在弱分量區(qū)間進(jìn)行參數(shù)估計(jì)，使用差分方法來(lái)平均化消除噪聲信號(hào)，并且考慮采樣點(diǎn)的非周期性，對(duì)于差分采用每?jī)牲c(diǎn)組合方式，共C2N 種，可得

4)用原始信號(hào)減去已估計(jì)出的弱信號(hào)分量區(qū)間形成新的信號(hào)，并在強(qiáng)分量區(qū)間進(jìn)行參數(shù)估計(jì)，可得

式中，iNHD1，iNHD2，NNHD1，NNHD2分別為對(duì)弱信號(hào)和強(qiáng)信號(hào)使用GIP 計(jì)算后確定的數(shù)據(jù)分界點(diǎn)和數(shù)據(jù)量。

5)設(shè)ej為誤差向量，每點(diǎn)數(shù)據(jù)為xj，則

使用該誤差向量的Frobenius 范數(shù)‖E‖F(xiàn)與門限值ε 進(jìn)行參數(shù)估計(jì)結(jié)果的評(píng)估，若在n 次迭代過(guò)程中各參數(shù)估計(jì)結(jié)果的變化量低于ε，并且此時(shí)‖E‖F(xiàn)處于多次迭代過(guò)程的極小值區(qū)間，即可認(rèn)為已達(dá)到“實(shí)際收斂”。

6)從估計(jì)出的各參數(shù)出發(fā)重新估計(jì)常數(shù)項(xiàng)C，去除估計(jì)信號(hào)后可再做一次GIP 運(yùn)算得到新的iNHD3與NNHD3，挑選適合重新計(jì)算常數(shù)C 的區(qū)間，即

3 實(shí)際數(shù)據(jù)驗(yàn)證分析

驗(yàn)證算法采用兩種方式進(jìn)行:一種方式使用固定信號(hào)加油耗指標(biāo)偏移噪聲信號(hào)，噪聲為高斯白噪聲，信噪比步進(jìn)5 dB，區(qū)間［-20 dB，20 dB］，再考慮一種信噪比無(wú)窮大時(shí)即沒(méi)有噪聲情況，評(píng)估考察算法誤差的最佳性能，噪聲信號(hào)疊加于所有飛行距離;另一種方式采用實(shí)際數(shù)據(jù)進(jìn)行，采集了某航空公司2012 年10 月至2013 年8 月A330 飛機(jī)飛行千米油耗樣本共8430項(xiàng)，對(duì)所有的樣本進(jìn)行原始點(diǎn)分布估計(jì)，使用本文方法進(jìn)行仿真驗(yàn)證。

3．1 算法的噪聲運(yùn)算性能

使用標(biāo)準(zhǔn)信號(hào)作為固定信號(hào)加噪聲影響驗(yàn)證此方法的通用性，對(duì)給定的最小信噪比RS，N，白噪聲en～N(0，σ2)，標(biāo)準(zhǔn)差σ=K1e-3×10(-RS，N/20)。在不同SNR條件下各參數(shù)50 組估計(jì)平均值與克拉美-勞(CRB)對(duì)比如圖5 所示。

圖5 參數(shù)估計(jì)精度與CRB 對(duì)比圖Fig．5 Parameter estimation precision and CRB under the condition of standard signal

從圖5 可看出，估計(jì)性能相對(duì)CRB 持續(xù)下降，算法效果顯著，從50 組平均均方根誤差與CRB 對(duì)比亦可看出上述關(guān)系，如圖6 所示。

由圖5、圖6 可看出，算法的運(yùn)算性能可以在不同信噪比條件下很好地估計(jì)出模型的各參數(shù)，在強(qiáng)噪聲環(huán)境下，因?yàn)槎啻蔚虶IP 的使用，也可以很好地排除噪聲的影響，逼近最佳性能。

當(dāng)信號(hào)為理想形式即不疊加噪聲時(shí)，算法收斂速度平滑減慢，震蕩逼近正確信號(hào)形式，圖7 所示為理想信號(hào)條件下誤差均方根在100 次迭代后震蕩收斂到原信號(hào)值，初始參數(shù)取值如表1 原始信號(hào)估計(jì)值，最佳收斂性能‖E‖F(xiàn)=0．0091。

圖6 不同信噪比下誤差均方根與CRB 對(duì)比Fig 6 RMSE under various SNR comparied with CRB

圖7 理想信號(hào)條件下迭代次數(shù)與誤差均方根變化趨勢(shì)Fig．7 The RMSE trends under constant iteration

3．2 飛行數(shù)據(jù)試驗(yàn)結(jié)果

如圖3 所示，A330 飛機(jī)所飛航段跨越了短中長(zhǎng)各個(gè)階段，從200 ～10000 km 均有涉及，總體符合前文模型所預(yù)示結(jié)果，其數(shù)據(jù)偏差似噪聲擾動(dòng)，并且噪聲功率較強(qiáng)，使用本文方法在經(jīng)過(guò)300 次迭代運(yùn)算后，選取誤差均方最小值點(diǎn)。各參數(shù)計(jì)算結(jié)果如圖8 與表1 所示。

圖8 A330 飛機(jī)8429 組數(shù)據(jù)實(shí)驗(yàn)結(jié)果Fig．8 8429 set of simulation data of A330

對(duì)所有原始數(shù)據(jù)進(jìn)行128 點(diǎn)交叉平均采樣處理，迭代100 次后可滿足所有收斂條件，數(shù)據(jù)趨于平滑，平滑噪聲影響后進(jìn)行估計(jì)，兩者估計(jì)參數(shù)與變化百分比計(jì)算結(jié)果如表1 所示。

表1 參數(shù)估計(jì)結(jié)果Table 1 Result of parameter estimation

可以看出，當(dāng)數(shù)據(jù)進(jìn)行多點(diǎn)平均后，對(duì)于原始參數(shù)變化沒(méi)有超出10%，說(shuō)明算法在平均化輸入后亦可得到較好的結(jié)果。

3．3 模型仿真結(jié)果推論

根據(jù)表1 平均后的數(shù)據(jù)計(jì)算結(jié)果，A330 飛機(jī)的油耗性能參數(shù)(單位為kg/km)及兩個(gè)距離常數(shù)(單位為km)分別為min{yA330}= 6． 21，max{yA330}= 6． 21 +20．01 + 2．7 = 28．92，RA330，1= 1/0． 006 2 = 161．29，RA330，2=1/0．000 252 =3 968．3。指數(shù)衰減項(xiàng)衰減到離性能參數(shù)下限值10%的飛行距離(單位為km)分別為rA330，1= ln(10 ×20． 01/6． 21)× RA330，1=560，rA330，2=ln(10 ×2．7/6．21)×RA330，2=5832。故根據(jù)模型數(shù)據(jù)可推論，A330 飛機(jī)在航段距離小于560 km 時(shí)起飛帶來(lái)的油耗高性能指標(biāo)占據(jù)主導(dǎo)因素，此距離范圍內(nèi)選擇該機(jī)型執(zhí)行航線油耗性能指標(biāo)較高。在5832 km 后可認(rèn)為已達(dá)到理想油耗性能，大于此距離時(shí)飛機(jī)可以忽略起飛段高指標(biāo)影響，油耗性能指標(biāo)平穩(wěn)偏低。

4 總結(jié)

油耗指標(biāo)模型是航空公司評(píng)價(jià)燃油消耗與節(jié)能減排的重要評(píng)估手段，基于RELAX 思想的信號(hào)分離，迭代估計(jì)算法具有一定的抗噪性能以及通用性，進(jìn)行油耗模型的精細(xì)化估計(jì)是可行的，為航空公司監(jiān)控機(jī)隊(duì)油耗水平等方面提供一個(gè)可行的方法，也可為航線機(jī)型選型提供一定的參考。

［1］ 中國(guó)民用航空局．2013 中國(guó)民航業(yè)發(fā)展公報(bào)［EB/OL］．(2013-06-01)． http://www． caac． gov． cn/I1/K3/201406/P020140623612275082363．pdf．(CAAC．Civil aviation development statistical bulletin of China in 2013［EB/OL］．trol，2006，13(3):11-14，28．)(2013-06-01)． http://www． caac． gov． cn/I1/K3/201406/P0201406236 12275082363．pdf．)

［2］ 孫宏，李峰，黎青松．民用航空航班直接運(yùn)行成本測(cè)算分析［J］．交通運(yùn)輸工程與信息學(xué)報(bào)，2007，5(1):1-5．(SUN H，LI F，LI Q S．Research of civil aviation flight direct operating cost［J］． Journal of Transportation Engineering and Information，2007，5(1):1-5．)

［3］ 劉婧．基于飛行數(shù)據(jù)分析的飛機(jī)燃油估計(jì)模型［D］．南京:南京航空航天大學(xué)，2010．(LIU J．The aircraft fuel estimation model based on flight data analysis［D］． Nanjing:Nanjing University of Aeronautics and Astronautics，2010．)

［4］ 羅杰．航空運(yùn)輸燃油效率影響因素分析［D］．廣漢:中國(guó)民用航空飛行學(xué)院，2010．(LUO J．Analysis of the factors affecting fuel efficiency of air transportation［D］． Guanghan:Civil Aviation Flight University of China，2010．)

［5］ LI G Z． Machine learning in fuel consumption prediction of aircraft［C］//The 9th IEEE International Conference on Cognitive Informatics (ICCI)，2010:358-363．

［6］ LI J，PETRE S． Efficient mixed-spectrum estimation with applications to target feature extraction［J］． IEEE Transactions on Signal Processing，1996，44(2):281-295．

［7］ MELVIN W L，WICKS M C． Improving practical spacetime adaptive radar［C］//IEEE National Radar Conference，Syracuse，New York，USA，1997:48-53．

［8］ 郭佳佳，廖桂生，楊志偉，等． 利用廣義內(nèi)積值迭代加權(quán)的空時(shí)協(xié)方差矩陣估計(jì)方法［J］．電子與信息學(xué)報(bào)，2014，36(2):422-427． (GUO J J，LIAO G S，YANG Z W，et al． Iterative weighted covariance matrix estimation method for STAP based on generalized inner products［J］． Journal of Electronics ＆ Information Technology，2014，36(2):422-427．)

［9］ 李艷，肖懷鐵，付強(qiáng)．Radon-Wigner 變換改進(jìn)算法在多目標(biāo)分辨及參數(shù)估計(jì)中的應(yīng)用［J］．電光與控制，2006，13(3):11-14，28． (LI Y，XIAO H T，F(xiàn)U Q． Application of improved Radon-Wigner transform in multi-target resolving and parameter estimation［J］． Electronics Optics ＆ Con-