陳昱君，孫樊榮，許學(xué)吉，沐 瑤

（南京航空航天大學(xué)民航學(xué)院，南京 211100）

在民航業(yè)蓬勃發(fā)展的背景下，中國航空公司的盈利情況卻不甚良好，其中，燃油正是制約航空公司（簡稱航司）盈利的重要因素之一。為了提高燃油利用效率，促進航司盈利，需要更加精準地計算飛機攜帶燃油量。在雷達管制十分普及的情況下，航空器不再按照傳統(tǒng)的進離場程序飛行，而是在一定空域范圍內(nèi)受管制員指揮，不同管制員的管制方式會形成不同的飛行路徑，造成了以往按照程序計算燃油消耗（簡稱油耗）的方式不再精準。綜上所述，基于航跡數(shù)據(jù)進行油耗分析與研究對提升航空公司的運行效益有著重要意義。

在飛機油耗的計算方面，國內(nèi)外學(xué)者已經(jīng)進行了大量研究。在國外，Khadilkar 等[1]利用飛行數(shù)據(jù)記錄器的數(shù)據(jù)，基于最小二乘回歸法以滑行時間、停車次數(shù)、轉(zhuǎn)彎次數(shù)和加速活動4 個變量建立了兩個線性模型對航空器滑行階段的油耗進行預(yù)測；Lawrance 等[2]通過對高斯過程回歸、K-近鄰回歸、簡單線性、新的基于空氣動力學(xué)的4 種模型采用高頻數(shù)據(jù)計算不同階段的K值并進行對比，結(jié)果表明新的基于空氣動力學(xué)的模型計算油耗更加精確；Oh 等[3]將地面燃油消耗與發(fā)動機推力狀態(tài)相關(guān)聯(lián)，且發(fā)動機推力狀態(tài)根據(jù)地面停車狀況的不同可分為峰值狀態(tài)和常規(guī)狀態(tài)，因此可通過國際民用航空組織（ICAO，International Civil Avi-ation Organization）數(shù)據(jù)庫建立基于停車轉(zhuǎn)彎的地面油耗模型與基于加速的地面油耗模型，并將這兩種模型的計算結(jié)果與QAR（quick access recorder）數(shù)據(jù)對比，結(jié)果顯示加速運動才是影響飛機滑行階段油耗的因素。在國內(nèi)：趙冬林[4]利用Pearson 系數(shù)進行相關(guān)性分析，結(jié)合真實的QAR 數(shù)據(jù)分析了航行滑出、滑入、爬升、巡航、下降5 個階段油耗與影響因素的關(guān)系，并建立了基于支持向量機的油耗模型，粗略地計算了下降階段的油耗，但忽略了進近階段的油耗；錢宇等[5]利用一元線性回歸分析方法對QAR 滑行時間數(shù)據(jù)進行擬合，提高了滑行時間的預(yù)測精度，從而對油耗計算進行優(yōu)化，但其只考慮了滑行時間一個影響因素，對于滑行時具體油耗估算并不準確；黃倩文等[6]基于BADA（user manual for the based of aircraft data）油耗模型，結(jié)合QAR 數(shù)據(jù)采用最小二乘法擬合氣動參數(shù)中的升力系數(shù)和阻力系數(shù)，并同時考慮側(cè)風(fēng)的影響對離場油耗模型進行了優(yōu)化，其計算精度高于風(fēng)洞計算。

在航跡數(shù)據(jù)處理方面：Gariel 等[7]利用基于雷達數(shù)據(jù)的K-means 算法，對終端區(qū)航跡進行聚類分析，有效地實現(xiàn)了對航空器的實時監(jiān)控；王莉莉等[8]利用局部異常因子（LOF，local outlier factor）算法對航跡點的局部可達性密度進行計算，對大于1 的航跡點進行野點甄別和排除；陳勇[9]通過line-segment-clustering 算法對每條航跡鄰域范圍內(nèi)航跡數(shù)量小于閾值的航跡進行剔除；韓云祥[10]基于航空器物理運動的6 個參量建立連續(xù)動態(tài)模型，并結(jié)合飛行管理系統(tǒng)航路點的8個離散變量建立離散切換模型，然后通過Petri 網(wǎng)進行融合，并加入風(fēng)場變量插值，從而擬合出航跡飛行的高度剖面和水平剖面；何艷等[11]利用均值及合成航跡線的思想對現(xiàn)有的航跡數(shù)據(jù)特征進行總結(jié)，并提出基于航跡形狀約束的半監(jiān)督K-means 算法對航跡進行聚類。

通過對上述研究總結(jié)發(fā)現(xiàn)：①計算航空器油耗的方式多種多樣，而基于BADA 數(shù)據(jù)庫進行油耗計算是被廣泛采用和認可的一種方法；②利用QAR 數(shù)據(jù)分析航行各個階段油耗已日益完善和精細，能夠?qū)叫懈鱾€階段進行詳細的動態(tài)分析計算，但下降進近階段由于雷達管制的普及，各個管制員的指揮方式相差較大，飛機的飛行軌跡差異大，對飛機的油耗影響不可忽略；③現(xiàn)階段研究都是從單個航空公司的QAR 數(shù)據(jù)庫進行分析，缺少來自空管雷達全部數(shù)據(jù)的綜合分析，存在一定的單一性和片面性。

針對上述研究的不足，提出采用BADA 數(shù)據(jù)庫提供的燃油計算方法，搭建下降進近階段的油耗模型，并通過大量的空管航跡數(shù)據(jù)計算油耗，對計算結(jié)果進行假設(shè)檢驗，分析燃油分布特點并計算出航司燃油攜帶建議值。

1 航跡數(shù)據(jù)處理

空中交通雷達航跡數(shù)據(jù)記錄了航空器航行的諸項元素，但由于獲取到的數(shù)據(jù)雜亂且存在異常數(shù)據(jù)以及缺失數(shù)據(jù)，為了能夠精準地計算油耗從而為航司提供科學(xué)的決策量，在利用雷達航跡數(shù)據(jù)進行油耗計算前，要對紛繁復(fù)雜的航跡數(shù)據(jù)進行分析和處理[12]。

1.1 航跡數(shù)據(jù)分析

一條完整的雷達航跡由若干個航跡段組成，而每個航跡段又是由單個孤立的航跡點構(gòu)成。從集合的角度看一條航跡表示為D={R1，R2，…，Rn}，其中，n 為航跡段個數(shù)，Ri為第i 個航跡段，Ri={p1，p2，…，pm}，表示該條航跡段由m 個航跡點組成，每個航跡點按照時間序列依次排列，由于航跡段長度不同，航跡點的個數(shù)也不相同。pj={t，la，lo，h，f，v}，表示第j 個航跡點是由以下6 項航行元素定義的：t 表示航空器位于該點的時間，la表示該點的緯度，lo表示該點的經(jīng)度，h 表示高度，f 表示航向，v 表示航空器的速度。

總而言之，航跡由航跡段構(gòu)成，航跡段是計算燃油的基礎(chǔ)單位，航跡段由若干個航跡點組成，每個航跡點的時間間隔為雷達波束的掃描間隔即4 s，而每個航跡點又是由6 項航行元素表征的。

1.2 航跡數(shù)據(jù)處理

由于獲取到的雷達數(shù)據(jù)是雜亂的，要想在海量的數(shù)據(jù)中提取有效因素進行油耗分析，首先要對數(shù)據(jù)進行規(guī)范化處理，即將原始數(shù)據(jù)可視化；由于研究范圍為進近階段，所以對原始數(shù)據(jù)進行范圍提取；由于雷達有一定容錯率，難免會有噪聲和缺失，所以對數(shù)據(jù)離群點進行剔除，并為了計算的連貫性利用插值法補全數(shù)據(jù)；最后進行特征點提取，即針對每條航跡進行拐點識別，從而對航跡分段來計算攜帶燃油消耗量[13]。航跡數(shù)據(jù)處理流程，如圖1 所示。

圖1 航跡數(shù)據(jù)處理流程Fig.1 Trajectory data processing process

1）規(guī)范化處理

用Python 對原始數(shù)據(jù)進行規(guī)范化處理，以期得到規(guī)整的.csv 格式文件，讀取excel 形式的.csv 數(shù)據(jù)文件，查看整體數(shù)據(jù)。

2）研究范圍截取

由于獲取到的航跡數(shù)據(jù)是整個管制區(qū)域的全部航跡數(shù)據(jù)，為了只保留計算所需數(shù)值，截取起始值為開始下降的高度，終止值為QFE（query field elevation）100 m。

3）離群數(shù)據(jù)和缺失數(shù)據(jù)處理

雷達數(shù)據(jù)中常含有噪聲和錯誤數(shù)據(jù)，常見的錯誤類型有數(shù)據(jù)突變和數(shù)據(jù)缺失。若時間、速度和高度有異常或缺失值則刪除整行數(shù)據(jù)。

航跡截取和異常數(shù)據(jù)處理流程如圖2 所示。

圖2 航跡截取和異常數(shù)據(jù)處理流程圖Fig.2 Flow chart of trajectory interception and abnormal data processing

圖2 中，航跡中航跡點最小Nmin=500，即當航跡點數(shù)少于500 時，將該條航跡作為錯誤數(shù)據(jù)刪除，航跡點最大高度Hmax=8 000 m，即文中只考慮8 000 m以下的進離場航班。

4）特征點提取

特征點是指高度改變的點。由于航空器在空中飛行難免遇到顛簸對高度產(chǎn)生影響，航空器處于平飛階段期間，雷達數(shù)據(jù)的高度顯示很難一成不變，所以通過觀察20 組進近航跡數(shù)據(jù)平飛階段各個點的高度誤差，計算得出平飛階段每相鄰兩個航跡點允許10 m 以內(nèi)的高度容差。在特征點的提取中也基于10 m 容差范圍對數(shù)據(jù)抖動情況進行排除，即滿足if data[i+1]-data[i]<-10 or if data[i]-data[i+1]<-10，則可視為data[i+1]=data[i]。對每條航跡進行分析，采用相鄰3點檢查，檢查某個點與其相鄰前后兩點高度的大小關(guān)系，判斷該點是否為拐點。若檢查結(jié)果滿足data[i-1]=data[i]and data[i+1] data [i] and data[i+1]=data[i]，則數(shù)據(jù)i 為下降轉(zhuǎn)平飛的拐點。如此反復(fù)可得到所有特征點。

2 基于BADA 數(shù)據(jù)庫的油耗計算模型

BADA 數(shù)據(jù)庫關(guān)于某機型發(fā)動機的油耗模型是利用機型性能參數(shù)表獲得航空器在平飛階段和下降階段的燃油消耗系數(shù)及推力系數(shù)，從而獲得單位時間內(nèi)的燃油消耗率[14]，再結(jié)合雷達數(shù)據(jù)給出的航跡點時間信息得到航跡數(shù)據(jù)對應(yīng)的燃油消耗量。

2.1 推力計算方法

2.1.1 最大爬升推力模型

1）標準大氣條件下的計算模型

最大爬升推力是計算各個階段推力的基礎(chǔ)，而最大爬升推力是根據(jù)真空速、標準海平面氣壓高度和溫差計算的。以下給出各類發(fā)動機在標準大氣（ISA，international standard atmosphere）條件下的最大爬升推力模型，分別給出各類型單臺發(fā)動機（ISA 條件下）最大爬升推力計算公式。

渦噴式

渦槳式

活塞式

式中：CTC，1、CTC，2、CTC，3為推力系數(shù)，由機型決定；Hp為標準海平面氣壓高度（1 ft=0.340 8 m），取航段開始高度；VTAS為真空速（1 kn=1.852 km/h），取航段平均值。

2）基于溫度修正的計算模型

在航空器實際運行過程中，由于高度或氣候等因素改變，運行外界溫度不可能一成不變。而溫度變化對發(fā)動機推力也會產(chǎn)生影響，因此給出了在實際大氣溫度條件下3 類發(fā)動機的最大爬升推力的溫度修正表達式。

實際大氣條件下的單臺發(fā)動機最大爬升推力為

式中（ΔT）eff為基于溫差的修正量，其修正方式為

最大爬升推力模型的約束條件為

式中：CTC，4、CTC，5為推力系數(shù)，由機型決定；ΔT 為實際大氣溫度與ISA 的差值（℃）。

2.1.2 平飛階段推力模型

根據(jù)定義，平飛階段推力等于阻力。然而真實運行過程中可用的最大推力是有限的。因此，實際運行過程中的最大平飛推力以最大爬升推力和相應(yīng)的推力系數(shù)得出。

平飛階段的推力模型為

式中：（Tcr）max為平飛階段中單臺發(fā)動機的最大推力（kN）；CTr為最大平飛推力系數(shù)，由機型決定，無量綱。

2.1.3 下降階段推力模型

下降階段的推力表達式根據(jù)高高度下降、低高度下降、進近形態(tài)和著陸形態(tài)采用不同的修正系數(shù)，但所用的計算方式相同，在公式中統(tǒng)一用j 指代飛行階段，其表達式如下

式中：Tdes，j為第j 下降階段推力（kN）；CTdes，j為第j 下降階段推力系數(shù)。

由于不同下降階段對應(yīng)的推力大小相差很大，所以BADA 數(shù)據(jù)庫關(guān)于CTdes，j的取值根據(jù)高度和飛行階段的不同給出以下4 種情況，如表1 所示，其中：Hdes為下降階段的參考高度（ft），Ht為航空器當前飛行高度（ft）。

表1 下降時推力系數(shù)與推力對應(yīng)表Tab.1 Correspondence table of thrust and its co-efficiency when descending

2.2 平飛階段燃油消耗率

平飛階段的燃油消耗率采用的是標準燃油流量，是關(guān)于推力油耗比和推力的表達，適用于除了巡航和慢車推力狀態(tài)的各個階段[15]。其中，單位時間單位推力的油耗是關(guān)于VTAS的函數(shù)，根據(jù)發(fā)動機類型的不同其表達式也不同。

1）渦噴式、渦槳式發(fā)動機燃油消耗率

式中：fnom為標稱燃油消耗率（kg·min-1）；THR為平飛階段的推力（kN）；η 為單位時間單位推力的燃油消耗量（kg·min-1·kN-1），分為以下兩種計算方式。

a）渦噴式發(fā)動機

式中Cf1，Cf2為第一單位推力燃油消耗系數(shù)。

b）渦槳式發(fā)動機

2）活塞式發(fā)動機燃油消耗率

2.3 下降階段燃油消耗率

1）渦噴式、渦槳式發(fā)動機燃油消耗率

根據(jù)BADA 數(shù)據(jù)庫，對于噴氣發(fā)動機和渦輪螺旋槳發(fā)動機，當飛機切換到下降形態(tài)，尤其是進近和著陸時，下降的推力部分停止[16]。因此，進近和著陸階段的燃料流量計算應(yīng)以標準燃料流量為基礎(chǔ)，并在必要時限于最小燃料流量。也就是說，下降階段的燃油消耗率fdes應(yīng)當以標準燃油消耗率和最小燃油消耗率的最大值為準，即

式中：fnom與上一節(jié)提到的計算方法一致；fmin為最小燃油消耗率，計算方式如下

2）活塞式發(fā)動機燃油消耗率

與活塞式發(fā)動機的下降條件相對應(yīng)的燃油流量規(guī)定為該機型的最小燃油流量也即fmin，此參數(shù)規(guī)定為一個常數(shù)，該常數(shù)由飛機的機型決定，即

式中fdes為下降階段燃油消耗率（kg·min-1）。

2.4 油耗計算模型

在假定各個航跡段燃油消耗率恒定的情況下，利用上述燃油消耗率模型及雷達數(shù)據(jù)中的時間、速度和高度信息對某條航跡的油耗進行計算。

由圖3 可見，一條下降航跡由若干個平飛階段和下降階段組成。每個階段的油耗等于燃油消耗率與時間的乘積，其中不同機型不同階段的燃油消耗率的計算方式已在第2 節(jié)的前述部分給出，每個航段的時間在識別航跡數(shù)據(jù)的特征點后也可直接得出，最后將平飛階段和下降階段的油耗疊加組合計算得出航空器下降進近階段所需的油耗總量。

圖3 下降進近航跡Fig.3 Descent approach trajectory

若某一條下降進近航跡命名為d，其由x 個平飛階段和y 個下降階段組成，每個航段的燃油消耗率保持不變，結(jié)合2.2 和2.3 中提出的計算方式，其油耗總量計算模型如下

式中：Ne為航空器發(fā)動機數(shù)量；ti、tj為第i 平飛階段和第j 下降階段的持續(xù)時間（min）為第i 平飛階段燃油消耗率（kg/min）為第j 下降階段燃油消耗率（kg/min）。

3 算例分析

3.1 基于BADA 數(shù)據(jù)庫的油耗計算

3.1.1 數(shù)據(jù)介紹

1）航跡數(shù)據(jù)

采用中國民用航空西南地區(qū)空中交通管理局云南分局提供的30 d 中同一時段11：00—12：00 之間進近階段的航跡數(shù)據(jù)作為研究對象，驗證所提出的進近階段的油耗計算模型及油耗數(shù)據(jù)分析方法的可行性和有效性。

2）機型數(shù)據(jù)

將依據(jù)采集的雷達數(shù)據(jù)進行油耗計算，經(jīng)數(shù)據(jù)篩選共涉及4 種機型，在BADA 數(shù)據(jù)庫中找到這些機型的性能參數(shù)，并摘錄其相應(yīng)的計算推力和油耗的參數(shù)，如表2 所示。

表2 各機型油耗參數(shù)表Tab.2 Fuel consumption of each type of plane

3.1.2 單條航跡算例驗證

圖4 為某架B737-800 的進近航跡示意圖，其進近開始高度為10 394 ft，進近結(jié)束高度為2 469 ft，共分為2 個平飛階段和3 個下降階段，公式中的高度采用各個航段的起始高度，真空速采用航段平均速度。按照第2 節(jié)給出的計算模型進行計算如表3 所示。

圖4 某架B737-800 進近航跡圖Fig.4 Approach trajectory chart of a B737-800

表3 B737-800 進近航跡油耗計算Tab.3 Fuel consumption calculation of an approaching B737-800

3.1.3 科學(xué)燃油決策量

利用SPSS 對各機型的油耗數(shù)據(jù)進行分析，計算自由度df 和顯著性P 值，如表4 所示，觀察其顯著性檢測的P 值均大于0.05，可以發(fā)現(xiàn)各機型油耗均服從正態(tài)分布。再以某架B737-800 為例進行進一步驗證，其Q-Q 圖如圖5 所示。由圖5 可見，油耗數(shù)據(jù)的分布非常接近于一條直線，因此確定B737-800 機型在同一時段的進近階段的油耗服從正態(tài)分布，假設(shè)成立。通過分析其余各個機型Q-Q 圖也均滿足條件，不再贅述。

圖5 B737-800 油耗Q-Q 圖Fig.5 Q-Q diagram of B737-800 fuel consumption

表4 各機型SPSS 正態(tài)性檢驗表Tab.4 SPSS normal test table of each type of plane

在驗證飛機油耗符合正態(tài)分布后，利用SPSS 進行雙邊95%置信度檢驗，得到各個機型的均值的95%置信區(qū)間的上限值，該值為滿足昆明長水國際機場航班進近階段95%情況下足夠用油的科學(xué)燃油決策量。其計算結(jié)果如表5 所示。

表5 各機型正態(tài)分布分析結(jié)果Tab.5 Analysis results of normal distribution of each type of plane

3.2 基于進近程序的油耗計算

以昆明長水國際機場IDPUG 進近程序為例進行計算。起始進近定位點高度為11 800 ft，中間進近定位點高度為10 800 ft，中間進近至最后進近為平飛階段，進近結(jié)束高度為機場標高之上2 000 ft，即8 901 ft，平均速度采用進近標準速度?；贗DPUG 進近程序設(shè)計的進近軌跡，可知飛機飛行的航段計算高度、航段長度、平均速度和航段持續(xù)時間，結(jié)合飛機飛行的油耗計算模型（式（17））可得出各型號飛機在進近階段需要的油耗總量。

進近程序各階段的基本飛行情況如表6 所示。

將表6 中的各飛行階段的燃油消耗率與航段持續(xù)時間相乘，然后求和得到各機型在進近階段的油耗。分別列出基于IDPUG 進近程序計算的油耗及通過大量航跡數(shù)據(jù)計算燃油攜帶建議值的對比，以及所計算的航跡數(shù)據(jù)的平均進近時間，如表7 所示。

表6 IDPUG 進近程序各階段基本飛行情況Tab.6 Basic flight conditions at each stage of the IDPUG approach procedure

表7 油耗對比Tab.7 Comparison of fuel consumption

3.3 計算結(jié)果對比分析

通過3.1 和3.2 的油耗計算結(jié)果可以直觀地看出，基于航跡數(shù)據(jù)所計算的油耗建議值大大超出基于程序計算的結(jié)果，主要原因如下。

1）進近時間和路徑差異大

由于雷達管制普及，航空器進近軌跡不再完全按照進近程序進場著陸，例如基于程序推算的進近時間大約為7 min，而所有雷達數(shù)據(jù)中的進近航跡的平均進近時間都已經(jīng)全部超過這個數(shù)值，由此可見在雷達管制下路徑的不確定性對飛機油耗的影響較大。

2）終端空域航空器相互之間影響大

隨著航空器數(shù)量的日益增加，排隊和等待的現(xiàn)象也日趨常見，基于航跡數(shù)據(jù)的油耗計算采用的是同一時段內(nèi)連續(xù)進近的航空器油耗數(shù)據(jù)，而基于進近程序計算時卻無法考慮航空器之間的間隔造成的額外航空器油耗。

3）航跡數(shù)據(jù)計算基于實際大氣環(huán)境

風(fēng)向、風(fēng)速和氣溫均是影響航空器油耗的因素，基于航跡數(shù)據(jù)進行油耗計算時，航空器的速度是受風(fēng)影響的真實速度，而基于程序的計算則是采用進近圖中推薦的理想速度所得出的數(shù)據(jù)，忽略了自然因素對航空器油耗造成的影響。

綜上所述，基于進近程序?qū)τ秃倪M行預(yù)測不再適用于當前的運行控制，從而更加凸顯對于同一機場、同一機型基于大量雷達航跡數(shù)據(jù)進行油耗預(yù)估的可信度。

4 結(jié)語

基于航跡數(shù)據(jù)的燃油消耗率計算方法可得到科學(xué)客觀的航空器燃油攜帶量，給航空公司決策航空器燃油攜帶量提供了有力的參考。為保證研究的客觀性，采集實際航跡數(shù)據(jù)，對航跡數(shù)據(jù)進行處理與分析，采用BADA 數(shù)據(jù)庫的油耗計算方式計算得出航空器下降進近階段的實際油耗，并給出建議攜帶燃油量，計算結(jié)果貼合實際，更適用于雷達管制已普及的現(xiàn)代航空運輸規(guī)劃。但目前只研究了下降進近階段的油耗，整個飛行階段更精確的油耗計算需要進一步的研究。