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

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

①(國(guó)防科技大學(xué)電子科學(xué)學(xué)院 長(zhǎng)沙 410073)

②(香港天文臺(tái) 香港 999077)

1 引言

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

相比于線性卡爾曼濾波，無(wú)跡卡爾曼濾波能夠模擬非線性系統(tǒng)變化且估計(jì)精度能夠達(dá)到泰勒級(jí)數(shù)展開(kāi)的2階精度，可應(yīng)用于復(fù)雜風(fēng)場(chǎng)演化預(yù)測(cè)。飛機(jī)尾流演化過(guò)程中受復(fù)雜風(fēng)場(chǎng)影響，其演化規(guī)律呈現(xiàn)非線性特點(diǎn)，其非線性演化傳遞和量測(cè)函數(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為量測(cè)矢量，vD為多普勒速度矢量，f(·)和h(·)分別為系統(tǒng)的非線性狀態(tài)函數(shù)和量測(cè)函數(shù)。

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

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

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

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

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

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

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

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

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

排除背景風(fēng)的情況下，雷達(dá)探測(cè)可得到掃描平面的風(fēng)場(chǎ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]，已證實(shí)Hallock-Burnham速度模型與實(shí)測(cè)數(shù)據(jù)吻合度較高[21]，而Rankine模型和Lamb-Oseen模型存在對(duì)環(huán)量的低估[22]，故采用經(jīng)典尾流速度模型Hallock-Burnham速度模型[23]來(lái)計(jì)算尾流的切向速度

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

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

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

5 實(shí)驗(yàn)驗(yàn)證分析

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

5.3 實(shí)測(cè)數(shù)據(jù)驗(yàn)證分析

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

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

(1) 探測(cè)場(chǎng)景1

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

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

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

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

(2) 探測(cè)場(chǎng)景2

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

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

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

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

6 結(jié)論

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

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