張 健,羅鑫悅,黎宗孝,趙嶷飛
(中國民航大學 空中交通管理學院,天津 300300)
當前,無人駕駛航空業(yè)正不斷快速發(fā)展、迭代演進。在人工智能等全新技術革命為主的第四次工業(yè)革命浪潮中,無人駕駛航空應運而生并蓬勃發(fā)展,已成為先進生產(chǎn)力的重要載體。隨著在農(nóng)業(yè)、物流、航拍乃至載人等通用和運輸航空領域的推廣應用,無人駕駛航空正在逐步融入國家空域系統(tǒng)。從行業(yè)細分市場來看,我國消費級無人機預計2024年市場規(guī)模將達到568億元,工業(yè)級無人機預計2024年市場規(guī)模達到1 508億元[1]。
然而,與無人機數(shù)量急劇增長的態(tài)勢相比,無人機管理規(guī)范與法律建設明顯滯后。中國民航局、JARUS、FAA、EASA相繼發(fā)布監(jiān)管規(guī)則[2-6]。將小型無人機的運行空域限制為120 m以下,以避免干擾有人機運行。該規(guī)則要求小型無人機必須在飛手視距內(nèi)或在具有有效通信能力的觀察員視距內(nèi)運行。Thomas等[7]、Jeffrey等[8]分別介紹以部分國家為代表的無人機主流運營機構(gòu),設想通過空域和航路設計等實現(xiàn)無人機超視距運行,構(gòu)建無人機空管系統(tǒng)。參照有人機在空中運行安全領域的發(fā)展路徑,對無人機安全管理而言,無人機空管系統(tǒng)建設的關鍵和難點在于開發(fā)類似于有人機配備的機載防撞系統(tǒng)的無人機沖突探測與解脫(conflict detection &resolution,CDR)技術;作為CDR的前提,航跡預測模型研究成為核心課題[9]。
CDR技術在有人機領域開展比較早,很多學者均在這方面做出過很多研究。James[10]結(jié)合美國林肯實驗室在機載防撞系統(tǒng)研發(fā)方面積累的經(jīng)驗,對有人機航跡預測及CDR建模方法進行了系統(tǒng)總結(jié)。然而對于小型無人機航跡預測的研究,國內(nèi)外均處于起步階段,主流方法包括基于運動學模型以及機器學習2種方法。
1)基于運動學模型
Vinay等[11]使用適用于有人機的BADA(base of aircraft data)模型對無人機航跡建模,生成無人機航跡預測的性能參數(shù),發(fā)現(xiàn)輕型、電動或垂直飛行無人機對模型適用性不足。Chen等[12]針對多無人機協(xié)同跟蹤多無人機編隊問題,提出了1種分散式協(xié)同跟蹤控制方法。James等[13]證明了無人機差異性對航跡建模誤差具有非常大的影響。D’Souza等[14-15]提出1種廣義六自由度多旋翼無人機航跡模型,用于有風狀態(tài)下無人機性能識別。
2)基于機器學習
Alligier等[16]提出1種基于機器學習的方法來學習無人機的運行因素,改進了無人機爬升預測的準確性。Wang等[17]提出1種基于飛行模式識別的混合方法。James等[18]利用神經(jīng)網(wǎng)絡學習無人機視覺特征,提出1種新的無人機探測方法。趙嶷飛等[19]從無人機運行狀態(tài)角度,實現(xiàn)航跡的狀態(tài)預測。李亦同[20]利用馬爾可夫過程,基于歷史數(shù)據(jù)進行航跡預測方法研究。陸佳歡等[21]使用長短期記憶神經(jīng)網(wǎng)絡算法和歷史數(shù)據(jù)進行訓練,提出無人機實時航跡預測算法。
綜上可知,針對無人機的航跡預測問題,已有研究均存在不足之處。無人機的機型種類繁多,且存在大量自行開發(fā)組裝,性能參數(shù)各不相同,加之制造商不愿公布參數(shù),采用質(zhì)點運動模型進行無人機航跡預測,研究中缺乏完整的無人機性能手冊?;跓o人機歷史航跡數(shù)據(jù),采用機器學習類方法進行航跡預測成為必然之舉。然而已有聚焦于歷史航跡數(shù)據(jù)的模型,存在預測長度和精度之間難以平衡的問題,僅用于短期預測。
鑒于此,本文采用貝葉斯網(wǎng)絡作為底層方法。首先,基于經(jīng)驗模型構(gòu)建了無人機航跡變量之間的靜態(tài)貝葉斯網(wǎng)絡和動態(tài)貝葉斯網(wǎng)絡。之后,統(tǒng)計歷史航跡離散化數(shù)據(jù),確定貝葉斯網(wǎng)絡中各個子節(jié)點變量的條件概率矩陣并進行降維存儲。最后,基于動態(tài)貝葉斯網(wǎng)絡,編程調(diào)用條件概率矩陣并反復迭代,最終實現(xiàn)航跡預測。
航跡數(shù)據(jù)包含的變量很多,且變量間存在多重耦合,例如,城市物流無人機低空飛行比高空飛行時速度要小得多、高速飛行比低速飛行時轉(zhuǎn)彎率要小得多。結(jié)合航跡預測功能的需要,涉及的主要變量是:高度數(shù)據(jù)、飛行速度、加速度、升降率和轉(zhuǎn)彎率?;诋斍俺鞘形锪鳠o人機航跡數(shù)據(jù),一般可以直接得到高度、速度,而加速度、轉(zhuǎn)彎率和升降率則需要推導獲取。
1)加速度變量。速度可以直接獲得或加速度分量整合獲得,通過對速度變量做差得到加速度數(shù)值。
2)升降率變量。通過無人機前后時刻高度變化獲得。
3)轉(zhuǎn)彎率變量。通過無人機航向前后時刻變化得到,或通過速度分量比值推導間接獲得。
貝葉斯網(wǎng)絡作為1種概率圖模型,已在民航安全評估領域進行不少應用。本文將貝葉斯網(wǎng)絡模型引入無人機航跡預測問題研究。首先,基于無人機性能和專家經(jīng)驗,構(gòu)建靜態(tài)貝葉斯網(wǎng)絡和動態(tài)貝葉斯網(wǎng)絡。接著,對無人機歷史航跡數(shù)據(jù)進行離散化處理,以訓練和構(gòu)建變量間條件概率矩陣。最后,通過預測數(shù)據(jù)的連續(xù)化和可視化處理,實現(xiàn)航跡預測功能。
貝葉斯網(wǎng)絡(Bayesian network,BN)被稱為因果概率網(wǎng)絡,基礎是貝葉斯公式。貝葉斯網(wǎng)絡各節(jié)點變量的條件概率是在貝葉斯公式與統(tǒng)計學的基礎上發(fā)展而來。設x1,x2,x3,…,xn,是變量X的1個互不相容的完備事件組,且p(xi)>0,Y是X引發(fā)的事件,貝葉斯公式如式(1)所示:
(1)
式中:xi為變量X的某個事件;xj為變量X的某個事件;p(xi)為事件xi的發(fā)生概率;p(Y|xi)是xi引發(fā)Y事件的條件概率。
動態(tài)貝葉斯網(wǎng)絡模型以離散變量集為研究對象,滿足馬爾科夫特性和平穩(wěn)特性。動態(tài)貝葉斯網(wǎng)絡可以看作是應用貝葉斯統(tǒng)計思想結(jié)合動態(tài)結(jié)構(gòu)的模型,既能夠體現(xiàn)變量之間的概率依存關系,又能描述這一系列變量隨時間變化的情況,是貝葉斯網(wǎng)絡在時間變化過程上的擴展。動態(tài)貝葉斯網(wǎng)絡作為1種機器學習方法,已被廣泛應用于語音識別、地面交通軌跡數(shù)據(jù)挖掘等領域。
貝葉斯網(wǎng)絡將每個變量視為節(jié)點,用有向邊表示它們之間的條件概率依賴關系,如圖1所示。
注:L(t)為t時刻高度,m;V(t)為t時刻速度,m/s;A(t)為t時刻加速度,m/s2;Ah(t)為t時刻升降率,m/s;RATE(t)為t時刻轉(zhuǎn)彎率,(°/s)。圖1 航跡預測中的靜態(tài)貝葉斯網(wǎng)絡遞推關系Fig.1 Recursive relationship of static Bayesian network in track prediction
由圖1可知,該網(wǎng)絡模型基于航空器性能參數(shù)模型和經(jīng)驗模型匯總得來,同時借鑒有人機在航跡預測模型中的參數(shù)關系。以升降率Ah(t)為例,影響其取值的父節(jié)點包括高度L(t)、速度V(t)、加速度A(t)。
動態(tài)貝葉斯網(wǎng)絡背后的假設是,系統(tǒng)的下一個狀態(tài)只取決于當前狀態(tài),也被稱為馬爾可夫假設?;诤娇掌餍阅芎蛯<医?jīng)驗,同時也借鑒有人機在航跡預測模型中的參數(shù)關系,得出如下動態(tài)貝葉斯網(wǎng)絡航跡模型,如圖2所示。
圖2 航跡預測中的動態(tài)貝葉斯網(wǎng)絡遞推關系Fig.2 Recursive relationship of dynamic Bayesian network in track prediction
由圖2可知,給定當前時刻高度、加速度、升降率、轉(zhuǎn)彎率4個變量值,根據(jù)動態(tài)貝葉斯網(wǎng)絡及各個變量的條件概率矩陣,即可實現(xiàn)下一時刻航跡預測和生成。
歷史航跡數(shù)據(jù)是變量取值的出處,也是進行模型訓練的基礎。本文基于貝葉斯網(wǎng)絡的模型求解數(shù)據(jù),來自2021年5月10日-2021年5月16日期間,順豐無人機于江西省贛州市某配送航線運行記錄。經(jīng)統(tǒng)計,日約26 000條,合計182 511條航跡數(shù)據(jù)。其中,每天存在約15航次飛行活動,飛行活動涵蓋離場上升、巡航、下降著陸3個階段。執(zhí)行每次飛行活動的無人機機型為Ark方舟無人機或H4四旋翼無人機2款。論文所采用的某城市物流無人機航跡數(shù)據(jù)格式如表1所示。
表1 無人機航跡數(shù)據(jù)格式Table 1 UAV track data format
基于前文航跡變量介紹,高度變量可以直接獲取,速度、加速度、轉(zhuǎn)彎率和升降率則通過推導獲得。
1)速度和加速度變量
速度通過已知速度分量整合獲得,如式(2)所示:
(2)
式中:VE(t)為東向速度,m/s;VN(t)為北向速度,m/s;VG(t)為地向速度,m/s。
對前后2個時刻速度做差獲取加速度數(shù)值,如式(3)所示:
A(t)=V(t+1)-V(t)
(3)
2)升降率變量
通過對無人機前后2個時刻高度數(shù)值做差獲得,如式(4)所示:
Ah(t)=L(t+1)-L(t)
(4)
3)航向角和轉(zhuǎn)彎率變量
將東向速度與北向速度做比,得到航向參考數(shù)據(jù)。設定北為Y軸正方向,東為X軸正方向,將將航向參考數(shù)據(jù)轉(zhuǎn)化為大地坐標系下航向角,如式(5)所示:
(5)
式中:θ(t)為無人機在大地坐標系航向角,(°)。
若VE(t)>0,VN(t)>0,則θ(t)=θ(t);若VE(t)<0,VN(t)>0,則θ(t)=360°-θ(t);若VE(t)>0,VN(t)<0,則θ(t)=180°-θ(t);若VE(t)<0,VN(t)<0,則θ(t)=180°+θ(t)。
進而計算轉(zhuǎn)彎率,如式(6)所示:
RATE(t)=θ(t+1)-θ(t)
(6)
3.2.1 條件概率矩陣的訓練
基于貝葉斯網(wǎng)絡模型要求,對變量歷史數(shù)據(jù)取值進行離散化,并將這些量化區(qū)間用序號命名。例如,將高度劃分為[100,140)、[140,180)、[180,220)、[220,260]4個區(qū)間。針對高度、速度、加速度、升降率、轉(zhuǎn)彎率離散值,運用歷史數(shù)據(jù)統(tǒng)計和訓練,生成各節(jié)點變量的條件概率矩陣。過程為:
1)針對不同機型,分別篩選出各自飛行活動所產(chǎn)生的航跡數(shù)據(jù),用于后續(xù)分別針對特定機型的航跡數(shù)據(jù)訓練。
2)針對貝葉斯網(wǎng)絡中的每個變量取值,基于數(shù)據(jù)分布情況分別定義一系列的分割點c1c2…cn對變量數(shù)據(jù)進行量化。數(shù)值小于c1的量化到第1個區(qū)間,大于cn的量化到第n+1個區(qū)間,以此類推,各個變量分別被量化為1至n+1的離散值。
3)針對某型無人機每次航線飛行活動航跡數(shù)據(jù),基于貝葉斯網(wǎng)絡中變量間約束條件,分別統(tǒng)計各種約束關系下的出現(xiàn)次數(shù)。以轉(zhuǎn)彎率為例,動態(tài)貝葉斯網(wǎng)絡中的當前時刻轉(zhuǎn)彎率受到上一時刻高度(4個離散值)、當前時刻升降率(4個離散值)、上一時刻轉(zhuǎn)彎率(5個離散值)3個變量約束和決定,要分別針對各組合情況統(tǒng)計當前時刻轉(zhuǎn)彎率取值概率情況,生成條件概率矩陣。
4)針對該型無人機的每次航線飛行活動,分別進行第3步操作,并以之對同型航空器當天其他飛行航次進行預測。選取預測效果最匹配的1次條件概率矩陣,作為當天的條件概率矩陣。
動態(tài)貝葉斯網(wǎng)絡里的轉(zhuǎn)彎率概率矩陣部分樣例展示如表2所示。
表2 轉(zhuǎn)彎率條件概率矩陣Table 2 Conditional probability matrix of RATE
3.2.2 條件概率矩陣的存儲
多重變量條件概率矩陣非常龐雜,為有效存儲和調(diào)用,提出了多維矩陣降維方法。首先,建立空白矩陣,列數(shù)由該變量離散量化區(qū)間數(shù)決定,行數(shù)由其父節(jié)點個數(shù)和各自離散量化區(qū)間數(shù)決定。以轉(zhuǎn)彎率的條件概率矩陣存儲為例,行數(shù)設定時,將上一時刻高度取值定義為百位,將當前時刻升降率取值定義為十位,將上一時刻轉(zhuǎn)彎率取值定義為個位。結(jié)合各節(jié)點區(qū)間個數(shù),確定行數(shù)為445行,列數(shù)為5列。之后將條件概率矩陣中的每一行概率值導入特定行中,取值范圍之外數(shù)據(jù)做空白處理,如表3所示。
表3 轉(zhuǎn)彎率條件概率矩陣存儲格式Table 3 Storage format for conditional probability matrix of RATE
針對無人機升降率、轉(zhuǎn)彎率和加速度,均可建立類似條件概率矩陣。在之后迭代過程中,通過調(diào)用變量條件概率矩陣,依照子節(jié)點取值概率分布生成預測值。
通過初始貝葉斯網(wǎng)絡和高度初始數(shù)據(jù),動態(tài)貝葉斯網(wǎng)絡中的高度、加速度、升降率、轉(zhuǎn)彎率均獲得初始數(shù)據(jù),之后在每輪迭代中生成加速度、升降率和轉(zhuǎn)彎率的預測值。在此過程中,高度沒有參與預測,需要由升降率新值連續(xù)化,更新高度原值生成高度新值,再對高度信息進行離散化,以驅(qū)動下一輪預測過程。
預測結(jié)束后,需要將升降率、加速度、轉(zhuǎn)彎率的離散值進行逆向連續(xù)化處理。直接進行連續(xù)化處理,會導致隨機取值的波動性過大,毛刺化嚴重。為此,本文提出1種連續(xù)化重復采樣方法,即如果離散值在一段時間內(nèi)始終不變,則只取其對應量化區(qū)間內(nèi)同一隨機數(shù)值。其中,若量化區(qū)間內(nèi)存在0值,則取0作為該時段數(shù)值。以升降率逆向連續(xù)化為例進行展示,如圖3所示。
圖3 升降率連續(xù)化重復采樣Fig.3 Repeated continuous sampling of Ah
由圖3可知,連續(xù)化處理可以最大程度降低毛刺現(xiàn)象。上述處理后,將升降率、加速度依次疊加到高度和速度變量上,間接完成高度和速度的連續(xù)化,如式(7)~(8)所示:
L(t+1)=L(t)+Ah(t+1)
(7)
V(t+1)=V(t)+A(t+1)
(8)
針對預測得到的無人機航跡數(shù)據(jù),進行可視化處理,建立三維坐標系。其中,Z軸取值地面垂直向上為正,X軸取值向東為正,Y軸取值向北為正。在XY平面上,設定航向初始值,將預測出來的轉(zhuǎn)彎率變化量不斷疊加,迭代生成航向新值?;谶\動學方程,生成無人機在X軸、Y軸、Z軸的實時位置,如式(9)~(12)所示:
x(t+1)=x(t)+V(t)cosθ(t)
(9)
y(t+1)=y(t)+V(t)sinθ(t)
(10)
z(t)=L(t)
(11)
(12)
式中:x(t)為無人機在X軸坐標;y(t)為無人機在Y軸坐標;z(t)為無人機在Z軸坐標;θ0為無人機初始航向。
預測航跡與歷史數(shù)據(jù)對比如圖4所示。
圖4 預測航跡和真實航跡對比Fig.4 Comparison between predicted track and real track
從圖4可知,預測航跡與真實航跡相似度很高。對預測效果進行分析,如表4所示。
表4 航跡模型的預測性能Table 4 Prediction performance of track model
由表4可知,基于動態(tài)貝葉斯網(wǎng)絡的無人機航跡預測數(shù)值與真實航跡數(shù)據(jù)總體符合性較好,總體預測精度接近94%。但航跡對比圖顯示部分航段符合度不夠高,進一步數(shù)據(jù)分析也驗證這個觀察。
1)本文首次將動態(tài)貝葉斯網(wǎng)絡引入城市物流無人機航跡預測研究中,構(gòu)建無人機航跡模型,完成無人機航跡變量多個離散取值的動態(tài)迭代和預測,最終實現(xiàn)無人機航跡預測功能。
2)數(shù)據(jù)量化處理和逆向連續(xù)化采樣存在近似處理,數(shù)據(jù)離散化處理時存在量化區(qū)間數(shù)量確定上的不足,后續(xù)通過收集更多的無人機運行數(shù)據(jù)加密量化區(qū)間。此外,由于每次對變量離散預測值的連續(xù)化采樣具有隨機性誤差,進而導致推導出的高度、速度、航向存在誤差,后續(xù)考慮通過優(yōu)化重復采樣方法降低該誤差。
3)不同飛行階段和不同型號無人機預測效果存在差異。因無人機不同飛行階段展示的機動差異,使得預測航跡在上升轉(zhuǎn)入巡航階段存在較大誤差,后續(xù)需要針對城市物流無人機不同飛行航段,分別進行精細化預測。此外,文中航跡預測模型是基于同型無人機構(gòu)建,換言之,不同型無人機需要分別構(gòu)建預測模型。當前物流運輸無人機航線飛行是在飛控系統(tǒng)控制下自動飛行,因此可基本排除人為操作習慣引起的差異。大風等天氣也會對無人機飛行航跡保持造成影響,因此,論文中針對不同自然日,分別進行模型訓練和構(gòu)建。后續(xù)研究可融入天氣因素,完善模型預測效果。
4)本文提出的無人機航跡模型,作為1種普適性航跡預測方法,可以實現(xiàn)針對任意航線、任意機型的航跡預測。隨著飛行活動增加,必然會打破當前無人機航線隔離劃設的現(xiàn)狀。因此,實現(xiàn)特定航線、特定機型的航跡預測,可為不同型號無人機之間交叉運行時沖突探測和解脫技術提供基礎。