沈 淳 李健兵* 高 航 陳柏緯 韓啟光 王雪松

①(國防科技大學(xué)電子科學(xué)學(xué)院 長沙 410073)

②(香港天文臺 香港 999077)

1 引言

飛機(jī)尾流是機(jī)翼上下表面壓力差而在其后方形成的一種反向旋轉(zhuǎn)的強(qiáng)烈氣流(如圖1所示)，是飛機(jī)飛行時產(chǎn)生升力的必然產(chǎn)物，具有空間尺度大、持續(xù)時間長、旋轉(zhuǎn)強(qiáng)烈等特點。在航空安全管制方面，飛機(jī)尾流會對后續(xù)進(jìn)入起飛/降落滑道或預(yù)定航線的飛機(jī)造成影響，可能發(fā)生顛簸、橫滾，乃至失去控制。在大氣物理方面，飛機(jī)尾流的存在會改變大氣中水的相態(tài)，產(chǎn)生凝結(jié)尾跡、管道云等自然現(xiàn)象，甚至可能影響地球熱輻射和小范圍的天氣[1]。飛機(jī)尾流的特性、探測和演化規(guī)律研究已成為當(dāng)前空中交通安全管制、大氣物理等領(lǐng)域共同關(guān)注的前沿科學(xué)問題。

圖1 飛機(jī)尾流形成示意圖Fig.1 Illustration of aircraft wake vortex

在空中交通管制領(lǐng)域，為了避免飛機(jī)因遭遇前機(jī)產(chǎn)生的飛機(jī)尾流而發(fā)生事故，國際民航組織(International Civil Aviation Organization,ICAO)于20世紀(jì)70年代針對不同飛機(jī)的起飛重量規(guī)定了相鄰兩架飛機(jī)起飛/降落的最小安全距離，但這些規(guī)則比較保守，很大程度上限制了機(jī)場的起降容量[2–4]，難以滿足航空業(yè)快速發(fā)展的需求。歐美等國家在此基礎(chǔ)上開展了飛機(jī)重分類(Recategorization,RECAT)[5,6]研究并形成適合當(dāng)?shù)氐匦魏蜌庀蟓h(huán)境條件的分類體系。中國民航引入了歐美飛機(jī)重分類并開始在廣州和深圳機(jī)場試運行，但歐美RECAT規(guī)則主要針對歐美機(jī)場運行需求開發(fā)，具有地形和氣象環(huán)境條件約束，我們對該規(guī)則背后物理規(guī)律的認(rèn)識還不充分。我國正處于自主創(chuàng)新、快速發(fā)展的歷史機(jī)遇期，開展飛機(jī)尾流短時行為預(yù)測，為后續(xù)飛機(jī)尾流危害評估和構(gòu)建適合中國民航的尾流動態(tài)間隔標(biāo)準(zhǔn)奠定基礎(chǔ)，對我國掌握相關(guān)核心技術(shù)和突破國際壟斷，具有重大的應(yīng)用價值和科學(xué)意義。

飛機(jī)尾流行為預(yù)測是根據(jù)飛機(jī)參數(shù)(重量、翼展、飛行速度、飛行高度等)以及氣象環(huán)境參數(shù)(橫向風(fēng)、風(fēng)切變、大氣湍流、大氣分層等)對其飛行時產(chǎn)生的左右渦旋的演化和消散進(jìn)行反演，預(yù)知飛機(jī)尾流的運動路線和強(qiáng)度變化情況，提前確定飛機(jī)尾流危害區(qū)域，達(dá)到有效縮短飛機(jī)起降安全間隔，提升航空效率和機(jī)場容量的目的。

基于上述目的，歐美等國家開展了長期研究并建立了部分尾流行為預(yù)測模型，歐洲的德國宇航中心(Deutsches Zentrum furLuft-und Raumfahrt,DLR)的Deterministic/Probabilistic Two-Phase wake vortex decay and transport model (D2P/P2P)[7,8]模型，將尾流演化分成初始和快速消散兩個階段，考慮了飛機(jī)重量、翼展、飛行速度、位置、飛行航跡傾角等飛行參數(shù)和迎風(fēng)速度、橫向風(fēng)、風(fēng)切變、湍流消散率等天氣參數(shù)。比利時魯汶天主教大學(xué)(Universite Catholique de Louvain,UCL)[9]建立了DVM/PVM(Deterministic/Probabilistic wake Vortex Model)模型，該模型基于離散尾渦粒子系統(tǒng)并結(jié)合Monte-Carlo仿真生成尾流渦旋和近地鏡面渦旋，考慮了飛機(jī)位置、翼展、重量、飛行速度等參數(shù)。美國的Sparpkaya,Robins,Delisi等研究團(tuán)隊[10]聯(lián)合構(gòu)建了APA(AVOSS Prediction Algorithm)模型，該模型考慮了橫向風(fēng)、氣溫、湍流強(qiáng)度等天氣參數(shù)并利用鏡像渦理論模擬近地面效應(yīng)。美國西北研究協(xié)會(Northwest Research Associates,NWRA)[11]建立了VIPER (Vortex algorithm Including Parameterized Entrainment Results)模型，利用從大渦模擬(Large Eddy Simulation,LES)中參數(shù)化學(xué)習(xí)建立的方程分別對左旋和右旋兩個尾渦進(jìn)行預(yù)測。TDAWP(TASS Driven Algorithm for Wake Prediction)[12]模型由美國國家航空航天局(National Aeronautics and Space Administration,NASA)建立，該模型基于控制容積法和流體動量、角動量守恒原理，可在不同天氣條件和飛機(jī)飛行參數(shù)下進(jìn)行預(yù)測。

上述尾流行為預(yù)測模型雖然考慮了多種氣象條件因素和飛機(jī)參數(shù)的影響，但是依賴初始風(fēng)場氣象條件的設(shè)置，在預(yù)測尾流演化過程中風(fēng)場演化非常復(fù)雜。隨著時間的變化，初始?xì)庀髼l件和飛機(jī)參數(shù)的設(shè)定與實際氣象條件等存在誤差，導(dǎo)致渦心位置和速度環(huán)量等趨勢預(yù)測與尾流實際演化趨勢存在偏差，進(jìn)而影響了飛機(jī)尾流行為預(yù)測的準(zhǔn)確性和魯棒性[13,14]。在此基礎(chǔ)上，Pruis等人[15]提出綜合運用多個預(yù)測模型，給出尾流演化的概率性預(yù)測結(jié)果，但這種綜合方法仍然依賴于初始?xì)庀髼l件的設(shè)置，未能解決氣象條件實時變化的問題。Sch?enhals等人[16]提出應(yīng)用數(shù)據(jù)融合的方法，利用線性卡爾曼濾波或擴(kuò)展卡爾曼濾波對尾流演化過程中預(yù)測和實測數(shù)據(jù)之間的誤差進(jìn)行估計并對預(yù)測數(shù)據(jù)進(jìn)行修正，但其仍以經(jīng)典D2P模型為尾流預(yù)測模型，并未根據(jù)實際風(fēng)場演化來實時調(diào)整預(yù)測模型氣象參數(shù)，修正魯棒性和后續(xù)預(yù)測效果不佳。

本文設(shè)計了非線性卡爾曼濾波數(shù)據(jù)同化預(yù)測模型，對尾流預(yù)測模型氣象環(huán)境參數(shù)進(jìn)行實時擬合修正，能根據(jù)實時探測數(shù)據(jù)對代表飛機(jī)尾流演化趨勢的渦心位置和速度環(huán)量等進(jìn)行實時修正并重新預(yù)測。

針對經(jīng)典尾流演化預(yù)測模型比較依賴初始?xì)庀髼l件和飛機(jī)飛行參數(shù)等初始條件的問題，結(jié)合復(fù)雜風(fēng)場非線性演化的特點，本文提出了基于無跡卡爾曼濾波(Unscented Kalman Filter,UKF)數(shù)據(jù)同化的飛機(jī)尾流行為預(yù)測方法，實現(xiàn)對飛機(jī)尾流行為演化的實時修正。針對經(jīng)典尾流演化預(yù)測模型初始?xì)庀髼l件和飛機(jī)飛行參數(shù)設(shè)置后未根據(jù)實際風(fēng)場演化實時調(diào)整的問題，應(yīng)用風(fēng)場線性切變和最小二乘曲線擬合方法，根據(jù)實時探測獲得的尾流風(fēng)場數(shù)據(jù)計算得到預(yù)測模型的氣象環(huán)境參數(shù)，構(gòu)建了參數(shù)化尾流實時演化預(yù)測模型。

本文的組織結(jié)構(gòu)如下：第2節(jié)對基于數(shù)據(jù)同化的飛機(jī)尾流行為預(yù)測思路進(jìn)行描述；第3節(jié)根據(jù)風(fēng)場線性切變和最小二乘擬合理論建立參數(shù)化飛機(jī)尾流演化預(yù)測模型；第4節(jié)設(shè)計數(shù)據(jù)同化預(yù)測模型，用于尾流行為預(yù)測的修正；第5節(jié)通過仿真和實測數(shù)據(jù)驗證基于數(shù)據(jù)同化的飛機(jī)尾流行為預(yù)測方法的有效性；第6節(jié)給出本文的結(jié)論。

2 飛機(jī)尾流行為預(yù)測方法

2.1 激光雷達(dá)探測尾流場景設(shè)置

在晴空條件下，飛機(jī)尾流回波散射主要由尾流內(nèi)部的浮塵等微粒決定，激光雷達(dá)波束的衰減較小，探測距離較遠(yuǎn)且時間和角度分辨率較高，所以激光雷達(dá)是當(dāng)前探測晴空飛機(jī)尾流廣泛使用的傳感器[17]。激光雷達(dá)常用的掃描方式為距離高度指示器(Range-Height Indicator,RHI)掃描，通過上下掃描飛機(jī)尾流可以得到掃描平面上各個探測單元在不同掃描時刻的飛機(jī)尾流回波信號(如圖2所示)。激光雷達(dá)探測飛機(jī)尾流一般采用側(cè)視方式，雷達(dá)設(shè)置于機(jī)場跑道一側(cè)，在與跑道正交的掃描平面上上下來回掃描飛機(jī)尾流，激光波束的仰角范圍為[αmin-αmax]。當(dāng)飛機(jī)在跑道上起飛或降落，將會在機(jī)場跑道上方的起飛/降落滑道上產(chǎn)生一對旋轉(zhuǎn)方向相反的漩渦，它們的渦心位置分別為Oc1和Oc2。

圖2 激光雷達(dá)探測飛機(jī)尾流場景設(shè)置Fig.2 Geometry configuration for Lidar detection of wake vortex

激光雷達(dá)掃描飛機(jī)尾流時，尾流中懸浮的粒子會受到激光的照射并散射能量。激光雷達(dá)接收粒子散射回來的光信號后，將其轉(zhuǎn)化為電信號，再通過相干處理，得到尾流區(qū)域的多普勒信息(如圖3所示)。通過處理多普勒信息，可以反演飛機(jī)尾流的渦心位置和速度環(huán)量等特征參數(shù)，為飛機(jī)尾流行為預(yù)測提供數(shù)據(jù)支撐。

圖3 飛機(jī)尾流多普勒速度RHI回波分布圖Fig.3 Doppler velocity distribution of wake vortex in an RHI

2.2 基于數(shù)據(jù)同化的飛機(jī)尾流行為預(yù)測方法

基于數(shù)據(jù)同化的飛機(jī)尾流行為預(yù)測(Data Assimilation,DA)方法如圖4所示。應(yīng)用激光雷達(dá)探測數(shù)據(jù)反演得到的尾流渦心位置和速度環(huán)量等特征參數(shù)，基于經(jīng)典尾流預(yù)測模型和非線性卡爾曼濾波，構(gòu)建參數(shù)化尾流預(yù)測模型進(jìn)行尾流演化行為預(yù)測，同時利用演化過程中實時探測到的尾流特征數(shù)據(jù)進(jìn)行非線性卡爾曼濾波，對預(yù)測的渦心位置和速度環(huán)量等特征參數(shù)進(jìn)行修正并重新進(jìn)行行為預(yù)測，同時利用實測數(shù)據(jù)對預(yù)測模型的氣象環(huán)境參數(shù)進(jìn)行曲線擬合修正，使尾流行為預(yù)測更加貼近真實。具體預(yù)測步驟如下：

圖4 飛機(jī)尾流行為預(yù)測流程圖Fig.4 Wake vortex behavior prediction flow chart

(1) 根據(jù)激光雷達(dá)實測獲得的尾流初始渦心位置和速度環(huán)量等特征參數(shù)，以及雷達(dá)和風(fēng)速儀等設(shè)備測量反演得到的風(fēng)速、風(fēng)切變系數(shù)、湍流耗散率等風(fēng)場數(shù)據(jù)，應(yīng)用參數(shù)化尾流演化預(yù)測模型設(shè)計非線性狀態(tài)傳輸函數(shù)，預(yù)測得到飛機(jī)尾流在一定時間段內(nèi)的渦心演化和強(qiáng)度變化趨勢。

(2) 在飛機(jī)尾流演化過程中，根據(jù)接收到的實時測量數(shù)據(jù)，將速度環(huán)量和渦心水平位置、高度等作為預(yù)估量，應(yīng)用無跡卡爾曼濾波，將非線性狀態(tài)傳輸?shù)玫降奈擦魈卣鲄?shù)預(yù)測值與實測數(shù)據(jù)進(jìn)行數(shù)據(jù)同化，計算對應(yīng)的狀態(tài)估計值、估計方差以及卡爾曼濾波增益系數(shù)，對代表尾流演化趨勢的速度環(huán)量、左右渦心位置等特征參數(shù)進(jìn)行修正。同時，根據(jù)實測數(shù)據(jù)進(jìn)行最小二乘曲線擬合更新修正速度環(huán)量演化模型的氣象環(huán)境參數(shù)，根據(jù)實測風(fēng)場數(shù)據(jù)進(jìn)行渦心演化模型的背景切向風(fēng)場速度和尾渦下降速度等參數(shù)的實時估計，利用修正后的氣象環(huán)境參數(shù)和速度參數(shù)對數(shù)據(jù)同化的非狀態(tài)傳遞函數(shù)進(jìn)行修正。

(3) 根據(jù)修正后的數(shù)據(jù)同化模型和尾流渦心位置、速度環(huán)量等特征參數(shù)對尾流演化趨勢進(jìn)行更新預(yù)測得到更加貼近真實的預(yù)測數(shù)據(jù)，演化過程中，一旦得到新的尾流實測數(shù)據(jù)，就重復(fù)第(2)步和第(3)步。

3 參數(shù)化尾流演化預(yù)測模型

參數(shù)化飛機(jī)尾流演化預(yù)測模型在經(jīng)典尾流演化模型基礎(chǔ)上，根據(jù)雷達(dá)連續(xù)探測反演得到的數(shù)個渦心位置和速度環(huán)量等特征參數(shù)進(jìn)行最小二乘擬合得到速度環(huán)量預(yù)測模型參數(shù)，根據(jù)實時探測到的風(fēng)場數(shù)據(jù)并利用線性切變得到渦心演化模型參數(shù)，進(jìn)而對飛機(jī)尾流的強(qiáng)度(速度環(huán)量)和渦心位置的演化趨勢進(jìn)行估計，并在演化過程中利用數(shù)據(jù)同化對渦心位置、速度環(huán)量以及預(yù)測模型的氣象環(huán)境參數(shù)進(jìn)行不斷修正，實現(xiàn)速度環(huán)量和渦心位置的實時行為預(yù)測。

3.1 速度環(huán)量演化預(yù)測

經(jīng)典D2P模型[7]將尾流速度環(huán)量演化分為穩(wěn)定段和快速消散段兩個階段，模型如式(1)

尾流在實際演化過程中，受低空風(fēng)場各種因素的疊加影響，強(qiáng)度變化趨勢復(fù)雜，但經(jīng)典D2P模型的參數(shù)設(shè)置并未考慮風(fēng)場演化因素的實時影響，無法真實反映出尾流強(qiáng)度的變化趨勢，需根據(jù)實時探測反演特征值對速度環(huán)量演化模型進(jìn)行修正。本參數(shù)化預(yù)測模型將等模型參數(shù)作為未知量，通過雷達(dá)連續(xù)探測反演得到的N(N ≥4)個初始環(huán)量值，結(jié)合式(1)進(jìn)行最小二乘曲線擬合求得，同時在演化過程中結(jié)合后續(xù)探測反演得到的環(huán)量值，不斷對等模型參數(shù)進(jìn)行最小二乘擬合修正。則由湍流耗散率、浮力頻率、大氣分層系數(shù)等實時探測得到的氣象條件參數(shù)決定。

3.2 左右渦心位置演化預(yù)測

飛機(jī)尾流受左右渦相互作用Vd不斷下降，同時受背景切向風(fēng)場Vc的影響進(jìn)行漂移。左右渦心相互作用的下降速度分別為

其中，b0為左右渦心間距，ΓL和ΓR分別為左渦和右渦速度環(huán)量值。

背景切向風(fēng)場Vc應(yīng)用線性切變理論進(jìn)行求解，在一個RHI掃描周期中，背景風(fēng)場由x和y方向的兩部分分量組成：

其中，x方向的分量可假設(shè)為線性切變：

圖5 估計背景風(fēng)場的非尾流區(qū)域Fig.5 Regions free of wake vortex that was used to estimate the background wind

根據(jù)求解得到的背景切向風(fēng)場Vc和渦心相互作用下降速度Vd，垂直方向的尾流渦心軌跡演化模型如式(5)

水平方向的尾流渦心軌跡演化模型如式(6)

在尾流演化過程中，使用上述演化模型來反演渦心軌跡垂直方向和水平方向的位置。在數(shù)據(jù)同化過程中，當(dāng)探測到新的雷達(dá)數(shù)據(jù)時，實時根據(jù)新的風(fēng)場數(shù)據(jù)進(jìn)行等尾流背景風(fēng)場參數(shù)的估計并重新進(jìn)行預(yù)測。

4 數(shù)據(jù)同化模型

數(shù)據(jù)同化模型基于無跡卡爾曼濾波，利用雷達(dá)探測得到的飛機(jī)尾流渦心位置和速度環(huán)量等特征參數(shù)對尾流行為預(yù)測進(jìn)行修正，達(dá)到根據(jù)實測數(shù)據(jù)進(jìn)行實時修正的目的。Julier等人[18]根據(jù)近似非線性函數(shù)的均值和方差遠(yuǎn)比近似非線性函數(shù)本身更容易的特點，提出了基于確定性采樣的UKF算法。

本文基于無跡卡爾曼濾波，針對尾流演化特點和雷達(dá)探測實際，設(shè)計了非線性傳遞函數(shù)和量測函數(shù)，采用無損(Unscented Transformation,UT)變換[19],利用比例修正對稱采樣策略獲取的Sigma采樣點通過非線性函數(shù)的傳遞進(jìn)行尾流行為預(yù)測，進(jìn)一步結(jié)合實測數(shù)據(jù)并利用量測函數(shù)實時修正行為預(yù)測結(jié)果。

4.1 基于無跡卡爾曼濾波的尾流特征參數(shù)數(shù)據(jù)同化

相比于線性卡爾曼濾波，無跡卡爾曼濾波能夠模擬非線性系統(tǒng)變化且估計精度能夠達(dá)到泰勒級數(shù)展開的2階精度，可應(yīng)用于復(fù)雜風(fēng)場演化預(yù)測。飛機(jī)尾流演化過程中受復(fù)雜風(fēng)場影響，其演化規(guī)律呈現(xiàn)非線性特點，其非線性演化傳遞和量測函數(shù)構(gòu)建如式(7)其中，xk=[Γk,1,Γk,2,xk,1,xk,2,yk,1,yk,2]為傳遞狀態(tài)矢量，Γk,1和Γk,2分別為左渦和右渦速度環(huán)量值，xk,1和xk,2分別為左右渦心水平位置，yk,1和yk,2分別為左右渦心垂直位置，zk=vD為量測矢量，vD為多普勒速度矢量，f(·)和h(·)分別為系統(tǒng)的非線性狀態(tài)函數(shù)和量測函數(shù)。

無跡卡爾曼數(shù)據(jù)同化的飛機(jī)尾流行為預(yù)測方法以比例修正對稱采樣策略選取Sigma點，利用Sigma點通過UT變換來近似尾流軌跡和速度環(huán)量非線性變化函數(shù)以及量測函數(shù)的均值和方差，其狀態(tài)演化傳遞流程如下。

假設(shè)k時刻的狀態(tài)估值和估計方差Pk|k，由比例修正對稱采樣策略，可得到2n+1個Sigma采樣矢量點以及對應(yīng)的權(quán)值然后將采樣點進(jìn)行非線性狀態(tài)函數(shù)傳遞得：對應(yīng)演化傳遞的均值和方差如式(8)

當(dāng)k+1時刻沒有雷達(dá)測量得到的量測值時，傳遞狀態(tài)矢量作為數(shù)據(jù)同化結(jié)果輸出為尾流演化預(yù)測結(jié)果；當(dāng)k+1時刻存在雷達(dá)測量得到的量測值時，根據(jù)狀態(tài)演化得到的及采樣策略公式可得2n+1個Sigma采樣矢量點和相應(yīng)的權(quán)值經(jīng)過非線性量測函數(shù)傳遞得量測變量的均值、方差以及協(xié)方差為

根據(jù)該時刻的實測向量zk+1可求出濾波增益Kk+1以及狀態(tài)估計和估計方差

4.2 非線性狀態(tài)傳遞函數(shù)和量測函數(shù)

4.2.1 狀態(tài)傳遞函數(shù)構(gòu)建

飛機(jī)尾流行為預(yù)測主要是對代表飛機(jī)尾流強(qiáng)度和位置的速度環(huán)量與渦心位置等特征參數(shù)演化進(jìn)行估計，系統(tǒng)非線性狀態(tài)傳遞函數(shù)f(·)為狀態(tài)向量的步進(jìn)預(yù)測

根據(jù)尾流預(yù)測模型，環(huán)量Γk，渦心水平位置xk和渦心垂直位置yk關(guān)于時間t*的導(dǎo)數(shù)為

4.2.2 量測函數(shù)構(gòu)建

排除背景風(fēng)的情況下，雷達(dá)探測可得到掃描平面的風(fēng)場多普勒速度分布，該多普勒速度為飛機(jī)尾流渦旋切向速度在雷達(dá)掃描波束上的投影(如圖6所示)，使用正弦定理可得多普勒速度與切向速度的關(guān)系式為

圖6 切向速度與多普勒速度的幾何關(guān)系Fig.6 Relationship between the tangential velocity and Doppler velocity

經(jīng)典尾流速度模型有Hallock-Burnham速度模型、Rankine速度模型和Lamb-Oseen速度模型等[20]，已證實Hallock-Burnham速度模型與實測數(shù)據(jù)吻合度較高[21]，而Rankine模型和Lamb-Oseen模型存在對環(huán)量的低估[22]，故采用經(jīng)典尾流速度模型Hallock-Burnham速度模型[23]來計算尾流的切向速度

其中，rc=0.052b0,b0為渦核大小。

將切向速度關(guān)系式(17)代入多普勒速度關(guān)系式(15)中，可得多普勒速度VD(R,α)的表達(dá)式為

在雷達(dá)實際探測過程中得到的多普勒速度為尾流速度分量和背景風(fēng)場分量的疊加，設(shè)計量測函數(shù)時需將背景風(fēng)場分量去除，因此量測向量zk設(shè)置為多個空間點去除背景風(fēng)后的多普勒速度矢量，其中一個空間點(R,α)的量測輸出為

5 實驗驗證分析

本節(jié)通過仿真實驗和實測數(shù)據(jù)對基于數(shù)據(jù)同化的飛機(jī)尾流演化趨勢預(yù)測方法的有效性進(jìn)行驗證分析。首先構(gòu)建飛機(jī)尾流仿真數(shù)據(jù)模塊，隨后使用兩步定位[24]和路徑積分方法[25]對預(yù)測需要用到的渦心位置和速度環(huán)量等尾流特征參數(shù)進(jìn)行估計，進(jìn)而對預(yù)測方法從仿真數(shù)據(jù)和實測數(shù)據(jù)兩個方面進(jìn)行驗證分析。

5.1 飛機(jī)尾流仿真數(shù)據(jù)模塊構(gòu)建

飛機(jī)尾流仿真數(shù)據(jù)模塊根據(jù)預(yù)先設(shè)定的飛機(jī)參數(shù)模擬飛機(jī)尾流的形成以及在復(fù)雜大氣環(huán)境中受風(fēng)場影響演化運動的過程，主要包含復(fù)雜風(fēng)場模擬模型、飛機(jī)尾流速度和演化模型以及激光雷達(dá)回波模型。

復(fù)雜風(fēng)場模擬模型包含橫向風(fēng)、風(fēng)切變和大氣湍流等，其中平均橫向風(fēng)是指風(fēng)速的基準(zhǔn)值，通常用來表明特定時間的風(fēng)速平均值；風(fēng)切變是指空間中任意兩點間風(fēng)向或風(fēng)速的突然變化，低空情況(距離地面高度小于600 m)下可劃分為水平風(fēng)的水平切變、水平風(fēng)的垂直切變和垂直風(fēng)的切變[26]。大氣湍流是介質(zhì)的不規(guī)則運動，其參數(shù)極其復(fù)雜且無規(guī)律變化，可采用基于馮·卡門(Von Karman)頻譜的擬小波方法來模擬隨機(jī)湍流場[27,28]。

飛機(jī)尾流穩(wěn)定段的速度可以簡化為左右兩個漩渦產(chǎn)生的速度場疊加的合速度，其演化情況受大氣溫度、湍流能量耗散率、風(fēng)切變以及尾流的初始強(qiáng)度、初始位置等多方因素影響[8]，演化模型采用經(jīng)典兩階段(Probabilistic Two-Phase,P2P)尾流耗散模型[7]來模擬飛機(jī)尾流的強(qiáng)度和位置在復(fù)雜大氣環(huán)境中的演化過程。

激光雷達(dá)接收的是探測單元內(nèi)大量隨機(jī)分布的氣溶膠的所有后向散射信號，忽略粒子之間的互耦，可認(rèn)為每個探測單元的后向散射信號為所有氣溶膠后向散射信號的相干疊加[29]。激光雷達(dá)信號處理把相鄰N個采樣回波信號補零到M點后再作M點的快速傅里葉變換(Fast Fourier Transform,FFT)得到多普勒頻譜。為提升回波的信噪比，可以對多個脈沖重復(fù)周期(Pulse Repetition Time,PRT)之間的多普勒頻譜進(jìn)行累積，以提高激光雷達(dá)的探測能力[30]，進(jìn)一步用矩估計方法可從頻譜的1階矩和2階矩得到回波信號的多普勒速度和多普勒譜寬等信息。更多關(guān)于飛機(jī)尾流仿真模塊構(gòu)建的細(xì)節(jié)可參考文獻(xiàn)[31]。

5.2 仿真數(shù)據(jù)驗證分析

應(yīng)用5.1節(jié)構(gòu)建的激光雷達(dá)探測飛機(jī)尾流仿真數(shù)據(jù)模塊，仿真總時長設(shè)置為90 s，左右尾流渦的初始速度環(huán)量都設(shè)置為Γ0=400 m2/s，左右渦心的初始位置分別為Oc1(t0)=(550,107)m和Oc2(t0)=(610,105)m，背景風(fēng)場包含橫向風(fēng)Vc0=-5 m/s、風(fēng)切變變化率β=0.05、縱向風(fēng)Vcy=-0.3 m/s和大氣湍流 EDR=0.05 m2/s3。激光雷達(dá)探測參數(shù)設(shè)置如表1所示。

表1 激光雷達(dá)探測參數(shù)設(shè)置Tab.1 Detection parameters of the Lidar

對仿真飛機(jī)尾流回波數(shù)據(jù)應(yīng)用兩部渦心定位方法和路徑積分方法進(jìn)行飛機(jī)尾流渦心定位和速度環(huán)量的計算，將定位和環(huán)量計算結(jié)果作為激光雷達(dá)探測數(shù)據(jù)輸入，對DA方法進(jìn)行驗證，同時加入數(shù)據(jù)融合[16](Data Fusion,DS)方法預(yù)測結(jié)果作為對比。DS方法采用經(jīng)典D2P模型對尾流進(jìn)行初步短時演化，同時根據(jù)實際探測值應(yīng)用線性卡爾曼濾波方法，將尾流速度環(huán)量、垂直距離和水平距離的估計誤差(實際探測值與經(jīng)典D2P模型預(yù)測值之差)作為修正量，建立關(guān)于該修正量的傳遞矩陣并進(jìn)行傳遞預(yù)測，采用泰勒級數(shù)展開方法計算誤差轉(zhuǎn)移矩陣，當(dāng)有測量數(shù)據(jù)時計算相應(yīng)的卡爾曼增益系數(shù)，通過迭代方法將模型預(yù)測數(shù)據(jù)與實測數(shù)據(jù)進(jìn)行融合得到修正量，應(yīng)用修正量對經(jīng)典D2P理論的預(yù)測結(jié)果和誤差傳遞矩陣進(jìn)行修正，數(shù)據(jù)融合流程如圖7所示，具體方法可見文獻(xiàn)[16]。

圖7 數(shù)據(jù)融合方法流程Fig.7 Flowchart of data fusion method

經(jīng)典飛機(jī)尾流預(yù)測模型根據(jù)初始設(shè)置氣象環(huán)境參數(shù)和飛機(jī)參數(shù)進(jìn)行尾流行為預(yù)測，為檢驗本文提出的數(shù)據(jù)同化方法對行為預(yù)測氣象環(huán)境參數(shù)估計錯誤修正的有效性和魯棒性，分別對經(jīng)典飛機(jī)尾流預(yù)測模型初始?xì)庀蟓h(huán)境參數(shù)的影響進(jìn)行分析。

(1) 渦心位置水平位置演化分析

圖8 飛機(jī)尾流行為預(yù)測方法對比(水平方向軌跡)Fig.8 Comparison between different wake vortex behavior prediction (Horizontal trajectories)

(2) 渦心位置垂直位置演化分析

圖9 飛機(jī)尾流行為預(yù)測方法對比(垂直方向軌跡)Fig.9 Comparison between different wake vortex behavior prediction (Vertical trajectories)

(3) 速度環(huán)量演化分析

湍流耗散率EDR由正確的0.05m2/s3設(shè)置為錯誤的0.01m2/s3，預(yù)測結(jié)果如圖10所示，因為輸入錯誤的氣象環(huán)境參數(shù)，經(jīng)典飛機(jī)尾流預(yù)測模型對尾流速度環(huán)量的預(yù)測存在一個固定的誤差，DA方法和DS方法都能夠根據(jù)實時探測得到的渦心位置和速度環(huán)量估計進(jìn)行修正，但DA方法能夠更快建立預(yù)測模型，得到的行為預(yù)測曲線相比DS方法更加貼近真實值。

圖10 飛機(jī)尾流行為預(yù)測方法對比(速度環(huán)量)Fig.10 Comparison between different wake vortex behavior prediction (Circulation)

進(jìn)一步對預(yù)測的誤差進(jìn)行對比，定義預(yù)測值與真實值之間的相對誤差，計算公式如下：

其中，Pci是一個RHI的估計位置或環(huán)量值，其對應(yīng)的理論值，N為總的RHI數(shù)量，i=1,2分別對應(yīng)左渦和右渦。當(dāng)計算渦心位置誤差時，Q=b0；當(dāng)計算環(huán)量誤差時，Q=Γ0。

設(shè)置不同的錯誤初始?xì)庀蟓h(huán)境參數(shù)，比較DA方法和DS方法應(yīng)用實測數(shù)據(jù)進(jìn)行實時修正階段的相對誤差，渦心位置相對誤差結(jié)果如表2所示，可以發(fā)現(xiàn)橫向風(fēng)誤差對經(jīng)典尾流預(yù)測模型影響較大，隨著橫向風(fēng)誤差變大，DA方法以及DS方法兩種行為預(yù)測模型誤差也會變大；縱向風(fēng)誤差和湍流耗散率誤差及其誤差大小的變化對行為預(yù)測影響較小，兩種行為預(yù)測模型都能夠達(dá)到較好的修正效果；在不同的誤差設(shè)置條件下，相比DS方法，DA方法能夠根據(jù)雷達(dá)實際探測數(shù)據(jù)進(jìn)行修正并取得更好的效果。

表2 飛機(jī)尾流預(yù)測位置相對誤差Tab.2 Relative error of predict trajectories

環(huán)量估計相對誤差如表3所示，可以發(fā)現(xiàn)不同的湍流誤差條件下，DA方法和DS方法都能夠保持良好的修正效果，這是由于湍流耗散率的誤差會導(dǎo)致經(jīng)典預(yù)測模型的誤差始終保持在恒定的誤差值，經(jīng)過一段時間的迭代，兩種行為預(yù)測方法都能夠?qū)⒃摵愣ㄕ`差進(jìn)行有效的修正；同樣可以發(fā)現(xiàn)DA方法能夠較快修正并得到反映真實反演趨勢的傳遞函數(shù)。綜上分析，DA方法能夠較快建立與尾流真實演化趨勢一致的傳遞函數(shù)，行為預(yù)測結(jié)果魯棒性較高。

表3 飛機(jī)尾流速度環(huán)量相對誤差Tab.3 Relative error of wake vortex circulation

5.3 實測數(shù)據(jù)驗證分析

香港天文臺使用多普勒激光雷達(dá)在香港國際機(jī)場開展了一系列的飛機(jī)尾流探測實驗，通過將雷達(dá)放置在機(jī)場不同位置進(jìn)行長時間的連續(xù)探測，驗證了激光雷達(dá)探測晴空飛機(jī)尾流的有效性。本小節(jié)使用在香港國際機(jī)場應(yīng)用激光雷達(dá)探測得到的多普勒回波數(shù)據(jù)對數(shù)據(jù)同化預(yù)測方法的有效性進(jìn)行分析，探測使用的激光雷達(dá)型號為windcube 200s，雷達(dá)在垂直于跑道的平面上做俯仰向掃描并接收得到該平面的RHI多普勒回波分布圖，探測系統(tǒng)參數(shù)設(shè)置如表4所示，更多探測場景細(xì)節(jié)描述可以參考文獻(xiàn)[32]。

表4 香港機(jī)場激光雷達(dá)探測參數(shù)設(shè)置Tab.4 Detection parameters of the Lidar in Hong Kong field campaigns

(1) 探測場景1

雷達(dá)放置在香港機(jī)場附近的亞洲博覽館天臺，距離地面高度約為19 m，與機(jī)場北跑道垂直距離約為400 m，雷達(dá)波束掃描角度范圍0.83°～10.77°，掃描角速度1.87°/s，探測場景如圖11所示。

圖11 香港機(jī)場北跑道激光雷達(dá)探測場景設(shè)置2014Fig.11 Geometry setup of the observation in north runway of Hong Kong international airport,2014

應(yīng)用兩步渦心定位和路徑積分方法反演得到飛機(jī)尾流渦心位置和速度環(huán)量等特征參數(shù)，并對得到的速度環(huán)量與前一相鄰時間點的速度環(huán)量的變化斜率進(jìn)行判斷，斜率絕對值超過閾值時(大于1)的數(shù)據(jù)點舍棄不用，以濾除應(yīng)用路徑積分方法出現(xiàn)的反演奇異值點。進(jìn)而將得到的特征參數(shù)用于驗證飛機(jī)尾流行為預(yù)測方法，結(jié)果如圖12所示，可以發(fā)現(xiàn)DA方法和DS方法都能夠根據(jù)實測數(shù)據(jù)對尾流演化趨勢進(jìn)行修正，但DS方法基于線性卡爾曼濾波模型，對于非線性雷達(dá)探測數(shù)據(jù)，估計出現(xiàn)非常大的偏差，而DA方法能夠根據(jù)實測數(shù)據(jù)進(jìn)行非線性模型修正，行為預(yù)測演化趨勢與雷達(dá)探測值更加貼近。DS方法修正曲線變化劇烈且偏離雷達(dá)探測值，與其使用經(jīng)典D2P模型進(jìn)行預(yù)測時依賴模型初始?xì)庀蟓h(huán)境參數(shù)設(shè)置，與真實風(fēng)場氣象環(huán)境差異較大，進(jìn)而導(dǎo)致模型的傳遞矩陣變化劇烈，無法利用線性卡爾曼濾波進(jìn)行有效修正有關(guān)。

圖12 飛機(jī)尾流行為預(yù)測方法對比Fig.12 Comparison between different wake vortex behavior prediction

(2) 探測場景2

雷達(dá)放置在機(jī)場跑道附近，距離地面高度約為7 m，與機(jī)場南跑道垂直距離約為240 m，雷達(dá)波束掃描角度范圍0～50°，掃描角速度5°/s，探測場景如圖13所示。

圖13 香港機(jī)場南跑道激光雷達(dá)探測場景設(shè)置2018Fig.13 Geometry setup of the observation in south runway of Hong Kong international airport,2018

同樣應(yīng)用兩步渦心定位和路徑積分方法反演得到飛機(jī)尾流渦心位置和速度環(huán)量等特征參數(shù)，將得到的特征參數(shù)用于驗證飛機(jī)尾流行為預(yù)測方法，結(jié)果如圖14所示，同樣可以發(fā)現(xiàn)DA方法預(yù)測效果遠(yuǎn)遠(yuǎn)好于DS方法且與雷達(dá)探測值接近。與探測場景1得到的結(jié)果類似，DS方法速度環(huán)量值的估計變化非常激烈且偏離探測值，與預(yù)測模型的傳遞矩陣不能匹配實際氣象環(huán)境變化，且探測值較少不能通過多次迭代進(jìn)行有效修正有關(guān)。相比DS方法，DA方法能夠通過少量的觀測值建立穩(wěn)定的尾流演化趨勢模型，特別是在尾流靠近地面時，由于地面效應(yīng)的影響兩個尾渦在垂直方向上出現(xiàn)不同程度的反彈，DA方法能夠擬合其變化趨勢并修正渦心位置，算法的有效性和魯棒性更高。

圖14 飛機(jī)尾流行為預(yù)測方法對比Fig.14 Comparison between different wake vortex behavior prediction

6 結(jié)論

本文通過數(shù)據(jù)同化方法研究了飛機(jī)尾流行為預(yù)測問題，根據(jù)非線性卡爾曼濾波模型建立了飛機(jī)尾流行為預(yù)測模型，通過雷達(dá)實時探測數(shù)據(jù)對飛機(jī)尾流行為預(yù)測軌跡進(jìn)行實時修正。結(jié)合經(jīng)典飛機(jī)尾流預(yù)測模型，根據(jù)風(fēng)場線性切變和最小二乘曲線擬合方法，構(gòu)建了參數(shù)化飛機(jī)尾流預(yù)測模型，解決了經(jīng)典模型初始化氣象參數(shù)未隨時間演化進(jìn)行實時調(diào)整的問題。在數(shù)值仿真數(shù)據(jù)和實測數(shù)據(jù)驗證中，建立了飛機(jī)尾流雷達(dá)回波仿真數(shù)據(jù)模塊，通過誤差對比分析，驗證了數(shù)據(jù)同化飛機(jī)尾流行為預(yù)測方法的有效性和魯棒性。

本文的工作是基于數(shù)據(jù)同化方法實現(xiàn)飛機(jī)尾流行為預(yù)測的嘗試，為針對國內(nèi)實際地形和氣象環(huán)境條件制定適合中國民用機(jī)場的飛機(jī)起降動態(tài)間隔標(biāo)準(zhǔn)提供了依據(jù)，為提高中國機(jī)場航班準(zhǔn)點率，節(jié)省燃油消耗開支，顯著提升社會效益等方面提供了技術(shù)支撐和理論基礎(chǔ)。未來將進(jìn)一步驗證數(shù)據(jù)同化模型在不同雷達(dá)實際工作條件下的有效性，以實現(xiàn)更加復(fù)雜風(fēng)場條件下的飛機(jī)尾流行為預(yù)測分析和危害評估。