復(fù)雜氣象條件下考慮結(jié)冰風(fēng)險(xiǎn)的無人機(jī)飛行策略

郭琪磊，桑為民，?？〗?，袁燁

1．西北工業(yè)大學(xué) 航空學(xué)院，西安 710072

2．中國民用航空飛行學(xué)院 工程技術(shù)訓(xùn)練中心，廣漢 618307

3．中國空氣動力研究與發(fā)展中心 結(jié)冰與防除冰重點(diǎn)實(shí)驗(yàn)室，綿陽 621000

4．中國商用飛機(jī)有限責(zé)任公司 上海飛機(jī)設(shè)計(jì)研究院，上海 201210

近年來無人機(jī)在物流運(yùn)輸、災(zāi)后救援、遙感測繪、航拍監(jiān)測、電力巡檢等諸多領(lǐng)域發(fā)揮著重要作用［1-2］。然而結(jié)冰問題日益成為威脅無人機(jī)飛行安全的重要因素。結(jié)冰問題是指航空器在積冰氣象環(huán)境下飛行時(shí)，由于云層中過冷水滴撞擊在其表面形成積冰而導(dǎo)致氣動性能下降的現(xiàn)象［3］。根據(jù)美國聯(lián)邦航空局（FAA）統(tǒng)計(jì)，2003─2008年間，有380余起飛機(jī)結(jié)冰相關(guān)的飛行事件［4］，而在無人飛機(jī)中，由于結(jié)冰造成的事故約占總事故的25%［5］。

無人機(jī)相比有人機(jī)而言遭遇結(jié)冰后更易導(dǎo)致嚴(yán)重后果，原因在于無人機(jī)抗結(jié)冰能力差，失速響應(yīng)時(shí)間更短，且可用于防除冰的能量更加有限，限制了防除冰裝置效能。因此，若能根據(jù)無人機(jī)實(shí)際應(yīng)用場景和結(jié)冰氣象環(huán)境，快速準(zhǔn)確地預(yù)測無人機(jī)的結(jié)冰情況，進(jìn)而對無人機(jī)航跡進(jìn)行優(yōu)化就顯得尤為重要。無人機(jī)日益深化的應(yīng)用場景要求其能夠在復(fù)雜的氣象條件下執(zhí)行飛行任務(wù)。就結(jié)冰問題而言，美國聯(lián)邦航空條例（FAR）25部附錄C中定義了大氣結(jié)冰條件，指出航空器結(jié)冰主要受結(jié)冰氣象環(huán)境影響，如液態(tài)水含量（Liquid Water Content，LWC）、平均有效水滴 直 徑（droplets Median Volumetric Diameter，MVD）、相對濕度（Relative Humidity， RH）、環(huán)境溫度、壓力以及云層范圍等［6］。

目前主要是通過考慮上述與結(jié)冰相關(guān)物理量的演化預(yù)報(bào)和診斷計(jì)算，利用結(jié)冰預(yù)報(bào)算法診斷飛機(jī)結(jié)冰情況。國際民航組織（International Civil Aviation Organization， ICAO）根據(jù)產(chǎn)生積冰時(shí)周圍大氣的溫度、濕度等條件，構(gòu)建了可預(yù)測結(jié)冰趨勢的IC積冰算法［7］。美國國家大氣研究中心（Na?tional Center for Atmospheric Research， NCAR）考慮到大氣的環(huán)境溫度、云中過冷水含量以及云滴的中位數(shù)體積直徑等參數(shù)而提出了積冰嚴(yán)重性指數(shù)［8］。美國空軍全球天氣中心開發(fā)的RAOB積冰算法［9］，將溫度劃分為不同的區(qū)間，并結(jié)合相應(yīng)溫度區(qū)間內(nèi)不同的溫露差和溫度垂直遞減率，進(jìn)而區(qū)分積冰強(qiáng)度和積冰類型。美國國家大氣研究中心提出的RAP積冰算法［10］，該算法根據(jù)某一高度層的溫度與相對濕度，將積冰定為4種類型。McDonough等［11］基于模糊邏輯和積冰情景決策樹分類方法建立了潛在積冰預(yù)報(bào)算法（Forecast Ic?ing Potential algorithm， FIP）。Bernstein等［12］將與結(jié)冰有關(guān)的濕度、溫度、云量、云水等要素與結(jié)冰可能性和結(jié)冰強(qiáng)度相聯(lián)系，從而提出當(dāng)前結(jié)冰潛勢（Current Icing Potential，CIP）算法。隨后在CIP算法基礎(chǔ)上結(jié)合地面觀測和探空資料，提出改進(jìn)CIP算法，并計(jì)算了全球各地的結(jié)冰及過冷大水滴的氣候分布［13］。王洪芳等［14］對比多種積冰算法，采用中尺度數(shù)值模式和預(yù)報(bào)產(chǎn)品，建立了飛機(jī)積冰預(yù)報(bào)模型。李佰平等［15］在CIP方法基礎(chǔ)上提出了綜合溫度、濕度、云頂溫度等要素的改進(jìn)結(jié)冰潛勢模糊邏輯診斷方法。然而上述各結(jié)冰預(yù)報(bào)算法均根據(jù)溫度、相對濕度、云量等結(jié)冰氣象參數(shù)對結(jié)冰潛勢做出模糊診斷，進(jìn)而預(yù)測出航空器結(jié)冰的可能性與結(jié)冰強(qiáng)度，但沒有對航空器結(jié)冰風(fēng)險(xiǎn)做較為準(zhǔn)確的量化評價(jià)。因此規(guī)劃出滿足結(jié)冰影響和其他飛行約束條件及性能指標(biāo)的最佳飛行航跡，進(jìn)而在復(fù)雜的氣象環(huán)境中能夠有效保證無人機(jī)的可靠性和行進(jìn)安全就顯得尤為重要。

本文提出一種復(fù)雜氣象條件下考慮結(jié)冰風(fēng)險(xiǎn)的無人機(jī)航跡規(guī)劃方法。首先，構(gòu)建基于中尺度WRF （Weather Research and Forecasting）模式的結(jié)冰氣象預(yù)測模型，通過基于最佳參數(shù)化方案組合的結(jié)冰氣象模擬獲得模擬時(shí)段內(nèi)目標(biāo)區(qū)域的溫度、壓力、LWC空間分布及時(shí)序變化。其次，構(gòu)建基于代理模型的水滴收集質(zhì)量快速預(yù)測方法。在獲取FAR 25部附錄C中連續(xù)最大結(jié)冰條件下40個(gè)采樣點(diǎn)處水滴收集質(zhì)量分布的基礎(chǔ)上，基于本征正交分解（Proper Orthogonal Decomposition，POD）降階模型和Kriging插值算法，建立溫度、壓力、LWC、MVD等結(jié)冰氣象參數(shù)與水滴收集質(zhì)量之間的代理模型，可快速預(yù)測出目標(biāo)區(qū)域水滴收集質(zhì)量的空間分布與時(shí)序變化。最后，針對現(xiàn)有結(jié)冰預(yù)報(bào)算法未能量化結(jié)冰風(fēng)險(xiǎn)的缺點(diǎn)，根據(jù)飛機(jī)結(jié)冰強(qiáng)度劃分等級，以不同結(jié)冰強(qiáng)度下水滴收集質(zhì)量閾值為結(jié)冰安全約束，利用基于粒子群優(yōu)化（Particle Swarm Optimization， PSO）的結(jié)冰容限航跡規(guī)劃方法進(jìn)行考慮結(jié)冰風(fēng)險(xiǎn)的無人機(jī)飛行策略研究。本文所提出的復(fù)雜氣象條件下考慮結(jié)冰風(fēng)險(xiǎn)的無人機(jī)飛行策略方法可根據(jù)數(shù)值預(yù)報(bào)的實(shí)時(shí)結(jié)冰氣象環(huán)境，快速地定量預(yù)測出目標(biāo)區(qū)域的結(jié)冰風(fēng)險(xiǎn)，并為無人機(jī)規(guī)劃出安全合理的航跡路線。

1 理論方法與數(shù)值模型

1.1 WRF模式配置

本文中結(jié)冰氣象預(yù)測采用WRF模式，版本為4.1.2。WRF模式是由美國國家環(huán)境預(yù)測中心（National Centers for Environmental Prediction，NCEP）和美國國家大氣研究中心等聯(lián)合開發(fā)的中尺度天氣數(shù)值預(yù)報(bào)系統(tǒng)。該模式采用可壓非靜力動力框架和地形追隨坐標(biāo)系，水平方向采用Arakawa-C型網(wǎng)格劃分計(jì)算域，具有多重網(wǎng)格嵌套功能和豐富的大氣物理過程參數(shù)化方案。結(jié)合各國氣象機(jī)構(gòu)提供的全球氣象數(shù)據(jù)和資料同化技術(shù)，WRF模式廣泛應(yīng)用于不同時(shí)間/區(qū)域尺度天氣現(xiàn)象模擬、空氣質(zhì)量建模、大氣海洋和理想化模擬等領(lǐng)域［16］，用于結(jié)冰氣象分析也有優(yōu)異表現(xiàn)［17-18］。

云層中上升氣流的強(qiáng)度是決定水凝體微觀物理特性的主要因素，已有研究表明在熱帶地區(qū)溫度低至?18 ℃的強(qiáng)上升氣流中仍然發(fā)現(xiàn)有明顯的過冷液滴存在［19］。此外，由于熱帶地區(qū)云層的云底溫度較高，而暖空氣比冷空氣需要更大的飽和含水量，因此相較于中/高緯度地區(qū)，熱帶地區(qū)云層中水汽含量也會更大［20］。海南位于中國南海北部，地處熱帶北緣，四面環(huán)海，海陸效應(yīng)明顯，強(qiáng)對流活動強(qiáng)烈且頻繁［21］。因此海南地區(qū)高空云層中的航空器飛行安全也容易受到結(jié)冰問題威脅。

本文中結(jié)冰氣象預(yù)測區(qū)域位于中國海南樂東地區(qū)（18.69N， 109.16E），模擬時(shí)段為UTC 2021-05-01T00：00至UTC 2021-07-31T00：00。WRF模式輸入資料采用NCEP發(fā)布的全球預(yù)報(bào)系統(tǒng)（Global Forecast System， GFS）生成的實(shí)時(shí)預(yù)報(bào)數(shù)據(jù)，其水平精度為0.25°×0.25°，時(shí)間分辨率為6 h。出于節(jié)省計(jì)算資源的考慮，模式采用2層嵌套方式，并遵循雙向嵌套策略，即粗糙的外部模型域（D01）為精細(xì)的內(nèi)部嵌套域（D02）提供邊界條件，而嵌套域在計(jì)算后將信息反饋到模型域。此外，嵌套區(qū)域水平網(wǎng)格分辨率分別為3 km（網(wǎng)格數(shù)為200×200）、1 km（網(wǎng)格數(shù)為190×190），垂直方向均勻劃分為34層。

基于前期大量數(shù)值試驗(yàn)并參考相關(guān)文獻(xiàn)，所確定WRF模式參數(shù)化方案如下所述：邊界層參數(shù)化采用Yonsei University（YSU）方案［22］，積云參數(shù)化采用Kain-Fritsch（K-F）方案［23］，長波/短波輻射分別采用RRTM （Rapid Radiative Transfer Model）方案［24］與Dudhia方案［25］，近地面層采用MM5相似方案［26］，陸面過程采用Noah陸面模型［27］。為優(yōu)化WRF模式參數(shù)化方案，采用4種微物理過程進(jìn)行敏感性分析工作，分別為WSM6方案［28］、Purdue-Lin方案［29］、Thompson方案［30］、Morrison方案［31］，具體4種微物理過程方案特征如表1所述。

表1 4種微物理過程方案特征Table 1 Characteristics of four microphysics schemes

基于WRF的結(jié)冰氣象預(yù)測模型的評價(jià)指標(biāo)分別采用絕對平均誤差（Mean Absolute Error，MAE）、相對平均誤差（Root Mean Absolute Er?ror， RMAE）以及相關(guān)系數(shù)（Pearson correlation coefficient，r），具體定義為

式 中：ymeas和ypred分 別 為 觀 測 值 與 預(yù) 測 值；yˉmeas和yˉpred分別為觀測值與預(yù)測值的均值；Np為樣本數(shù)量。

1.2 基于POD和Kriging的代理模型

本征正交分解（Proper Orthogonal Decomposition，POD）方法的核心是使用特殊正交基函數(shù)的疊加來實(shí)現(xiàn)未知參數(shù)的近似［32］。基本理論為假設(shè)某一空間Ω中存在物理場y=y(x，t)，可表示為一系列特征基函數(shù)與相關(guān)時(shí)間系數(shù)的線性組合：

式中：x和t為參數(shù)變量；αk(tk)為經(jīng)驗(yàn)系數(shù)；Φk(x)為特征基函數(shù)。借助Sirovich［33-35］引入的方法，取N個(gè)樣本集合為一組快照。找到這組快照所張成空間Ψ中的一組規(guī)范正交基使得集合中的元素在這組基上投影最大，即

其中：( · ， · )表示內(nèi)積運(yùn)算。

將式（6）代入式（5）可以得到一個(gè)新的特征值問題：

式中：A為N×N的協(xié)方差矩陣，其中各元素為Ai，j=(U(i)，U(j))；λ(1)，λ(2)，…，λ(N)為矩陣A從大到

矩陣A特征值的大小可以表征其對應(yīng)正交基在整個(gè)樣本空間中所占能量的比重。因此進(jìn)行物理場重構(gòu)時(shí)便可以僅保留部分正交基函數(shù)就能夠表征原物理場的主要特征。假設(shè)選用M個(gè)正交基函數(shù)，且有M≤N。在進(jìn)行矩陣分解得到最佳正交基后，對選取的正交基與樣本進(jìn)行內(nèi)積即可得到對應(yīng)的系數(shù)向量：

因此，重構(gòu)的物理場可以使用正交基的線性疊加表示為

給定n個(gè)輸入?yún)?shù)樣本點(diǎn)及其相對應(yīng)的系統(tǒng)響應(yīng)值其中X(i)為m維向量。Kriging模型假設(shè)目標(biāo)函數(shù)值與設(shè)計(jì)變量之間的真實(shí)關(guān)系可以表達(dá)為［36］

式中：F(β，X)為全局近似模型，為確定性部分；H(X)為一維隨機(jī)過程，其均值為0，方差為σ2，協(xié)方差表示為

其中：Z(X(i)，X(j))為輸入?yún)?shù)變量X(i)和X(j)的相關(guān)函數(shù)，其與兩點(diǎn)間空間位置密切相關(guān)，并隨著空間距離增大而減小。

式中：θq為第q個(gè)輸入?yún)?shù)分量xq的距離權(quán)值。

構(gòu)建Kriging代理模型的關(guān)鍵在于尋找最佳的θ向量，通過最大似然法將該問題轉(zhuǎn)化為如下的帶約束優(yōu)化問題：

式中：Z(θ)為相關(guān)函數(shù)矩陣，其第i行第j列元素為Z(X(i)，X(j))；σ2為方差估計(jì)值，其表達(dá)式為

其中：F為所有元素均為1的向量。

在求得最優(yōu)θ向量后，對任意一個(gè)新的輸入?yún)⒘肯蛄縓(0)，對應(yīng)的Kriging模型的估計(jì)值為

式中：z為n維相關(guān)向量，其第i個(gè)元素為Z(X(0)，X(i))。

1.3 基于PSO的結(jié)冰容限航跡規(guī)劃方法

本文中對無人機(jī)飛行策略的考察主要從航跡規(guī)劃角度出發(fā)。無人機(jī)航跡規(guī)劃是指在綜合考慮無人機(jī)自身物理約束及飛行區(qū)域中威脅和障礙物分布等諸多約束條件下，為無人機(jī)在給定的飛行區(qū)域中規(guī)劃出一條連接飛行起點(diǎn)和目標(biāo)終點(diǎn)的滿足所有約束條件的最優(yōu)或次優(yōu)的可飛路徑［37］，其中成功求解航跡規(guī)劃問題的關(guān)鍵在于航跡規(guī)劃算法研究。

粒子群優(yōu)化（Particle Swarm Optimization，PSO）算法是通過模擬鳥群覓食過程中遷徙和群聚行為的一種全局隨機(jī)優(yōu)化技術(shù)［38］，由于其基于種群的特性、易于實(shí)現(xiàn)和快速收斂等優(yōu)勢，因此被廣泛應(yīng)用于求解航跡規(guī)劃問題。PSO算法的基本原理為：假設(shè)D維搜索空間內(nèi)，K個(gè)粒子均被賦予位置和速度向量，即pi=[pi1，pi2，…，piD]Τ和Vi=[vi1，vi2，…，viD]Τ。尋優(yōu)過程中每個(gè)粒子都通過跟蹤個(gè)體極值Pi和全局極值Pg來更新位置。在第k次迭代過程中，可由已知的和可對粒子狀態(tài)進(jìn)行更新，即

式中：w為粒子的慣性權(quán)重系數(shù)；c1和c2為正實(shí)數(shù)，分別表示粒子的認(rèn)知和社會加速度系數(shù)；r1和r2為分布于［0，1］區(qū)間的隨機(jī)數(shù)。

在基于PSO的結(jié)冰容限航跡規(guī)劃問題中，采用文獻(xiàn)［39］中提出的方法建立規(guī)劃空間模型。如圖1所示，構(gòu)建全局笛卡兒坐標(biāo)系Oxy表示無人機(jī)航跡的規(guī)劃空間，其中點(diǎn)st和ta分別表示所規(guī)劃無人機(jī)航跡的起始位置和終點(diǎn)位置，陰影區(qū)域表示規(guī)劃空間中的禁飛區(qū)、障礙物、威脅源等。在所構(gòu)建的坐標(biāo)系中，起始位置st和終點(diǎn)位置ta沿x方向可以等分為ns+1等分，其中ns為規(guī)劃航跡導(dǎo)航節(jié)點(diǎn)數(shù)，為決策者預(yù)定義的常量。

圖1 所構(gòu)建的二維航跡規(guī)劃空間Fig. 1 Constructed 2D path planning workspace

本文中，所規(guī)劃航跡需滿足以下條件：① 所規(guī)劃航跡航程最短；② 所規(guī)劃航跡需滿足指定的結(jié)冰風(fēng)險(xiǎn)。因此單機(jī)航跡規(guī)劃問題的數(shù)學(xué)模型表示為

式中：ωi為任意待選航跡ph中的導(dǎo)航節(jié)點(diǎn)；JL為任意待選航跡的航程長度；Lmax為無人機(jī)的最大飛行航程；dis(ωi，ωi+1)為待選航跡ph中導(dǎo)航節(jié)點(diǎn)ωi與ωi+1之間的歐式距離；workplace為整個(gè)規(guī)劃空間；semi-free workplace為規(guī)劃空間中滿足結(jié)冰風(fēng)險(xiǎn)的半自由規(guī)劃空間。

鑒于本文的單機(jī)航跡規(guī)劃問題為約束優(yōu)化，引入任意待選航跡phm中粒子總的違反約束度TDm，在粒子初始化階段，為增加解的多樣性，TDm只需滿足容限要求即可，而在航跡迭代優(yōu)化階段，則要求待選航跡phm滿足最小航跡要求和結(jié)冰安全約束要求，即TDm=0［40］?？偟倪`反約束度TDm的計(jì)算式為

式中：Ds，m為該粒子違反航跡安全性能約束的程度；Df，m為該粒子違反無人機(jī)最大飛行航程限制的程度；Vm，j為待選航跡phm中第j個(gè)節(jié)點(diǎn)滿足結(jié)冰安全約束情況；Lm為待選航跡phm的航程長度。

基于PSO的結(jié)冰容限航跡規(guī)劃算法流程如算法1所示。

2 數(shù)值模型驗(yàn)證

2.1 基于WRF模式的結(jié)冰氣象預(yù)測

目標(biāo)區(qū)域內(nèi)結(jié)冰氣象預(yù)測分別采用表1中所示4種微物理過程所形成的參數(shù)化方案組合進(jìn)行模擬，并通過誤差分析確定最佳參數(shù)化方案。4種參數(shù)化方案所獲得的壓力和溫度預(yù)測值與東方氣象觀測臺站處的觀測值對比如圖2與圖3所示。圖中黑色線為東方氣象臺站觀測值，藍(lán)色線為Purdue-Lin方案預(yù)測值，洋紅色線為WSM6方案預(yù)測值，紅色線為Thompson方案預(yù)測值，綠色線為Morrison方案預(yù)測值，1 hPa=100 Pa?？芍?，盡管在個(gè)別時(shí)刻點(diǎn)處，4種微物理過程方案的預(yù)測值在數(shù)值精度上與實(shí)際觀測值存在偏差，但總體上4種微物理過程方案均可以十分準(zhǔn)確地預(yù)測東方氣象觀測臺站位置處壓力和溫度的時(shí)序變化趨勢。

算法1 基于PSO的結(jié)冰容限航跡規(guī)劃Algorithm 1 Ice tolerance route planning based on PSO

表2為不同微物理過程方案預(yù)測下誤差分析。由表中所述誤差可知，各參數(shù)化方案組合在壓力和溫度預(yù)測中均有出眾的表現(xiàn)，MAE僅約為0.80 hPa和0.98 K，RMAE也均小于0.08%和0.35%，說明各方案組合在壓力和溫度預(yù)測中具有足夠的數(shù)值精度。而各方案在壓力和溫度預(yù)測中的相關(guān)系數(shù)r也分別高達(dá)0.89和0.91，體現(xiàn)出了預(yù)測值與實(shí)際觀測數(shù)據(jù)極強(qiáng)的相關(guān)性。

圖2 不同微物理過程方案下壓力對比Fig. 2 Comparison of pressure in different microphysics schemes

圖3 不同微物理過程方案下溫度對比Fig. 3 Comparison of temperature in different microphysics schemes

表2 4種微物理過程方案預(yù)測下誤差分析Table 2 Error analysis of four microphysics schemes prediction

上述4種微物理過程方案組合在結(jié)冰氣象預(yù)測中，壓力和溫度的預(yù)測值與實(shí)際觀測值在數(shù)值精度和數(shù)據(jù)相關(guān)性方面均有出色的表現(xiàn)?？紤]相較于其他3種方案，Morrison方案不僅可以較好預(yù)測云水、云冰、雨、雪、霰等5種不同水成物的質(zhì)量濃度，還可以給出云冰、雪、雨和霰等粒子的分布圖譜，這為后續(xù)研究中準(zhǔn)確預(yù)測云中過冷水滴MVD提供了便利。因此，本研究選定Morri?son方案為最佳微物理過程方案，并用于目標(biāo)區(qū)域的結(jié)冰氣象預(yù)測。

圖4為LWC時(shí)序變化與文獻(xiàn)［41］中結(jié)果對比。圖中黑色圓點(diǎn)為實(shí)際觀測值，藍(lán)色線為本文中WRF模擬結(jié)果，紅色線為文獻(xiàn)［41］中模擬結(jié)果?？芍疚闹袛?shù)值模擬結(jié)果雖然在數(shù)值精度上與實(shí)際觀測值存在偏差，但能夠較好反映出實(shí)際觀測值的時(shí)序變化。此外本文和文獻(xiàn)［41］中LWC的MAE分別0.084 g/m3和0.099 g/m3，本文預(yù)測LWC的表現(xiàn)優(yōu)于文獻(xiàn)［41］中結(jié)果。因此本文所構(gòu)建WRF模式參數(shù)化方案組合能夠準(zhǔn)確預(yù)測出云層中的液態(tài)水含量。

圖4 時(shí)序變化下LWC對比Fig. 4 Comparison of LWC variation with time

2.2 基于代理模型的水滴收集質(zhì)量快速預(yù)測

根據(jù)1.2節(jié)中所述的代理模型，構(gòu)建水滴收集質(zhì)量快速預(yù)測方法，相應(yīng)的壁面水滴收集質(zhì)量計(jì)算式為

式中：LWC為云層中液態(tài)水含量；V∞為來流速度；β為局部水收集系數(shù)，且有β=A0Ain，其中A0與Ain分別為自由流場中由4條水滴軌跡所圍成的面積及在翼型表面所圍成的曲面面積。

對FAR 25部附錄C中連續(xù)最大結(jié)冰條件進(jìn)行優(yōu)化拉丁超立方采樣（Optimal Latin Hypercube Sampling， OLHS），獲得40個(gè)工況采樣點(diǎn)，采樣點(diǎn)分布如圖5所示。以NACA 0012翼型為模型，其弦長C為0.533 4 m，攻角AOA為4°，來流速度V∞為50~100 m/s。分別對40個(gè)采樣點(diǎn)處水滴撞擊特性進(jìn)行數(shù)值模擬，獲得不同工況下的水滴收集質(zhì)量分布，其中部分采樣點(diǎn)工況如表3所示，表中T為溫度、H為高度。

圖5 連續(xù)最大結(jié)冰條件下采樣點(diǎn)分布示意圖Fig. 5 Schematic diagram of sampling point distribution under continuous maximum icing conditions

表3 部分采樣點(diǎn)工況Table 3 Working conditions of partial sampling points

圖6為表3中采樣點(diǎn)工況下水滴收集質(zhì)量在翼面處分布對比，其中Y為與弦線垂直的方向，翼型前緣點(diǎn)為坐標(biāo)原點(diǎn)?？芍谝硇颓熬壧幩问占|(zhì)量最大：受來流角度影響，水滴收集質(zhì)量的最大值出現(xiàn)在翼型前緣下沿。而在前緣點(diǎn)之后沿著翼型表面水滴收集質(zhì)量逐漸減小，相應(yīng)地下翼面處水滴收集質(zhì)量變化相較于上翼面更為緩慢，水滴撞擊極限也相對更大。水滴收集質(zhì)量峰值與撞擊極限分布受LWC、MVD、溫度T、高度H、來流速度等因素綜合影響，其中水滴收集質(zhì)量極值分布受LWC和來流速度的影響較大，如Case 1與Case 4。而由于直徑較大、速度較大的液滴具有較大的慣性，在撞擊機(jī)翼表面時(shí)液滴軌跡不易發(fā)生改變。因此水滴撞擊極限范圍通常隨MVD和來流速度增大而增加，如Case 4與Case 5。

圖6 不同工況下水滴收集質(zhì)量沿翼面分布對比Fig. 6 Comparison of water droplet collection distribution along airfoil surface under different conditions

圖7 水滴收集質(zhì)量最大誤差隨積累特征值變化Fig. 7 Variation of the maximum error of water droplet collection with cumulative eigenvalues

圖7給出了壁面處水滴收集質(zhì)量最大誤差隨積累特征值的變化。圖中可以看出隨著積累特征值的增加，最大誤差迅速減小，當(dāng)積累特征值達(dá)到20時(shí)，水滴收集量最大誤差減小至接近于0。說明采用POD方法可利用少量特征基函數(shù)對原物理場進(jìn)行重構(gòu)。因此在本文中，使用前20個(gè)特征基函數(shù)進(jìn)行降階，即可滿足精度要求。而非采樣點(diǎn)（即預(yù)測點(diǎn)）處的特征系數(shù)，可通過對前20個(gè)系數(shù)分別進(jìn)行Kriging多維插值獲得。在建立代理模型后，選用9組工況進(jìn)行驗(yàn)證，分別考慮到影響結(jié)冰的溫度、高度、MVD、LWC和速度等5個(gè)參數(shù)，具體驗(yàn)證算例工況如表4所示。圖8示出了快速預(yù)測結(jié)果與數(shù)值模擬結(jié)果的水滴收集質(zhì)量對比?？梢钥闯鰺o論是水滴收集質(zhì)量的峰值還是分布范圍，快速預(yù)測結(jié)果均與數(shù)值模擬結(jié)果符合良好，說明所構(gòu)建的代理模型可以快速且準(zhǔn)確地預(yù)測出翼型表面水滴收集質(zhì)量。

圖8 數(shù)值模擬與快速預(yù)測水滴收集質(zhì)量對比Fig. 8 Comparison of water droplet collection between fast-prediction and simulation

3 結(jié)果與討論

本文采用WRF模式對海南樂東地區(qū)進(jìn)行了結(jié)冰氣象模擬，模擬時(shí)段為UTC 2021-05-01 T00：00至UTC 2021-07-31T00：00，每次模擬的前12 h作為模式的起轉(zhuǎn)時(shí)間。參數(shù)化方案與模式輸入資料的選擇詳見本文1.1節(jié)。

無人機(jī)結(jié)冰主要受氣象因素（如大氣溫度、飛行高度、LWC、MVD）以及飛行狀態(tài)（如速度）等因素影響?；赪RF模式的結(jié)冰氣象預(yù)測可獲得無人機(jī)在某一飛行環(huán)境中的大氣溫度、氣壓與LWC等參數(shù)。需要說明的是在本文中，當(dāng)大氣溫度<273.15 K時(shí)，則認(rèn)為液態(tài)水為過冷態(tài)，過冷水溫度與大氣溫度保持一致。根據(jù)圖5所示連續(xù)最大結(jié)冰條件與預(yù)測所得的LWC，進(jìn)行插值則可以獲得目標(biāo)區(qū)域內(nèi)MVD。圖9~圖11分別為UTC 2021-05-26T21：00、UTC 2021-06-05T01：00以及UTC 2021-07-11T08：00在高度層eta=25處目標(biāo)區(qū)域內(nèi)壓力、溫度與LWC分布云圖。

圖9 UTC 2021-05-26T21：00（eta=25）氣象參數(shù)分布Fig. 9 Meteorological parameter distribution at UTC 2021-05-26T21：00 （eta=25）

圖10 UTC 2021-06-05T01：00（eta=25）氣象參數(shù)分布Fig. 10 Meteorological parameter distribution at UTC 2021-06-05T01：00 （eta=25）

圖11 UTC 2021-07-11T08：00（eta=25）氣象參數(shù)分布Fig. 11 Meteorological parameter distribution at UTC 2021-07-11T08：00 （eta=25）

以上述3組結(jié)冰氣象條件為基礎(chǔ)，假定無人機(jī)飛行速度為50 m/s，利用所構(gòu)建的基于POD與Kriging的水滴收集質(zhì)量代理模型，分別對目標(biāo)區(qū)域內(nèi)水滴收集質(zhì)量進(jìn)行快速預(yù)測，獲得水滴收集質(zhì)量空間分布如圖12所示。

采用結(jié)冰強(qiáng)度作為考慮無人機(jī)航跡規(guī)劃問題中的結(jié)冰安全約束要求。結(jié)冰強(qiáng)度是指單位時(shí)間內(nèi)航空器表面所形成的冰層厚度，亦稱之為結(jié)冰速率。美國聯(lián)邦航空條例與航空信息手冊（AIM）將結(jié)冰強(qiáng)度劃分為4個(gè)等級，用以說明結(jié)冰情況的嚴(yán)重性。在本文研究中，將結(jié)冰速率乘以冰的密度即可獲得水滴收集質(zhì)量，取冰的密度為917 kg/m3，且認(rèn)為冰的密度始終為常數(shù)。飛機(jī)結(jié)冰強(qiáng)度等級劃分及水滴收集質(zhì)量等效轉(zhuǎn)換結(jié)果如表5所示。

圖12 水滴收集質(zhì)量空間分布Fig. 12 Spatial distribution of water droplet collection

表5 飛機(jī)結(jié)冰強(qiáng)度等級劃分Table 5 Aircraft icing intensity classification

以圖12中目標(biāo)區(qū)域內(nèi)3組水滴收集質(zhì)量分布作為無人機(jī)航跡規(guī)劃空間。由于3組航跡規(guī)劃空間內(nèi)最大水滴收集質(zhì)量均未超過0.02 kg/（s?m2），意味著最大結(jié)冰強(qiáng)度等級為中度結(jié)冰，因此分別選取表5中微量結(jié)冰與輕度結(jié)冰等級上極限，即水 滴 收 集 質(zhì) 量 為0.009 2 kg/（s?m2）與0.015 3 kg/（s?m2）作為無人機(jī)航跡規(guī)劃中結(jié)冰安全約束。當(dāng)所規(guī)劃航跡點(diǎn)水滴收集質(zhì)量小于上述上極限時(shí)，認(rèn)為該航跡點(diǎn)滿足結(jié)冰安全約束；否則違反結(jié)冰安全約束，該航跡點(diǎn)應(yīng)該規(guī)避。

利用所構(gòu)建的基于PSO的結(jié)冰容限航跡規(guī)劃方法分別規(guī)劃出上述3種情況下滿足最小航跡要求和結(jié)冰安全約束要求的離散航跡點(diǎn)，最終得到無人機(jī)可飛路徑。圖13~圖15分別為3種情況下所規(guī)劃航跡示意圖，圖中洋紅色線條表示滿足輕度結(jié)冰強(qiáng)度的航跡路線，綠色線條表示滿足微量結(jié)冰強(qiáng)度的航跡路線，紅色方塊點(diǎn)與綠色五角星點(diǎn)分別代表所規(guī)劃航跡路線的起始點(diǎn)與終止點(diǎn)。從圖中可以看出，在規(guī)劃滿足微量結(jié)冰強(qiáng)度的航跡路線時(shí)，由于結(jié)冰安全約束要求較強(qiáng)，所規(guī)劃航跡可規(guī)避開結(jié)冰可能性高的區(qū)域，然而代價(jià)是航跡距離相應(yīng)增大。而在規(guī)劃滿足輕度結(jié)冰強(qiáng)度的航跡路線時(shí)，由于結(jié)冰容限更高，因此可規(guī)劃出滿足給定結(jié)冰安全約束要求的最短航跡。3種規(guī)劃航跡的航程如表6所示。

圖13 UTC 2021-05-26T21：00（eta=25）航跡規(guī)劃Fig. 13 Route planning at UTC 2021-05-26T21：00（eta=25）

圖14 UTC 2021-06-05T01：00（eta=25）航跡規(guī)劃Fig. 14 Route planning at UTC 2021-06-05T01：00（eta=25）

圖15 UTC 2021-07-11T08：00（eta=25）航跡規(guī)劃Fig. 15 Route planning at UTC 2021-07-11T08：00（eta=25）

由式（23）可知，水滴收集質(zhì)量與氣象條件、飛行條件及幾何外形結(jié)構(gòu)有關(guān)。圖16為UTC 2021-06-05T01：00時(shí)刻不同飛行速度下的水滴收集質(zhì)量分布，與圖12（b）相比較，水滴收集質(zhì)量分布幾乎保持一致；然而隨著飛行速度增大至75 m/s和100 m/s時(shí)，空間分布中最大水滴收集質(zhì)量可達(dá)到0.035 8 kg/（s?m2）和0.052 7 kg/（s?m2），根據(jù)表5所示結(jié)冰強(qiáng)度等級劃分，此時(shí)已出現(xiàn)嚴(yán)重結(jié)冰條件。

表6 3種規(guī)劃航跡的航程Table 6 Voyage of three route trajectories

以圖16（a）中飛行速度為75 m/s目標(biāo)區(qū)域內(nèi)水滴收集質(zhì)量分布作為無人機(jī)航跡規(guī)劃空間。由于此時(shí)航跡規(guī)劃空間內(nèi)最大可能結(jié)冰強(qiáng)度等級為嚴(yán)重結(jié)冰，因此分別選取表5中微量結(jié)冰、輕度結(jié)冰與中度結(jié)冰等級上極限，即水滴收集質(zhì)量 為0.009 2、0.015 3、0.030 6 kg/（s?m2）作 為無人機(jī)航跡規(guī)劃中結(jié)冰安全約束。利用基于PSO的結(jié)冰容限航跡規(guī)劃方法，規(guī)劃得到不同結(jié)冰強(qiáng)度等級下的無人機(jī)可飛路徑，如圖17所示。此時(shí)，滿足微量、輕度與中度結(jié)冰等級的規(guī)劃路徑航程分別為253.1、226.6、204.9 km，可知由于滿足中度結(jié)冰等級時(shí)的結(jié)冰容限更高，因此規(guī)劃出滿足給定結(jié)冰安全約束要求的航跡航程最短。

圖16 不同飛行速度下水滴收集質(zhì)量分布Fig. 16 Distribution of water droplet collection at dif?ferent flight speeds

圖17 飛行速度為75 m/s時(shí)航跡規(guī)劃Fig. 17 Route planning at flight speed of 75 m/s

4 結(jié) 論

為解決無人機(jī)在復(fù)雜氣象條件下易受結(jié)冰影響，從而威脅其飛行安全的問題，提出了一種考慮結(jié)冰風(fēng)險(xiǎn)的無人機(jī)航跡規(guī)劃方法?；谥谐叨忍鞖忸A(yù)報(bào)WRF模式對海南樂東地區(qū)結(jié)冰氣象環(huán)境進(jìn)行預(yù)測，通過參數(shù)化方案敏感性分析確定了最佳參數(shù)化方案組合，進(jìn)而獲得該時(shí)段內(nèi)目標(biāo)區(qū)域的溫度、壓力、LWC空間分布及時(shí)序變化。與此同時(shí)，利用OLHS方法對FAR 25部附錄C中連續(xù)最大結(jié)冰條件進(jìn)行采樣，并對40個(gè)采樣點(diǎn)進(jìn)行水滴撞擊特性計(jì)算，獲得各采樣點(diǎn)處水滴收集質(zhì)量分布；基于POD降階模型和Kriging插值方法，構(gòu)建溫度、壓力、LWC、MVD等氣象參數(shù)與水滴收集質(zhì)量間的代理模型?；诮Y(jié)冰氣象預(yù)測參數(shù)與代理模型，獲得海南樂東地區(qū)水滴收集質(zhì)量的空間分布與時(shí)序變化。最后，分別以不同結(jié)冰強(qiáng)度下水滴收集質(zhì)量閾值作為結(jié)冰安全約束，利用基于PSO的結(jié)冰容限航跡規(guī)劃方法對無人機(jī)進(jìn)行考慮結(jié)冰風(fēng)險(xiǎn)的飛行策略研究。研究結(jié)果表明：

1） 利用WRF模式可獲得目標(biāo)區(qū)域內(nèi)溫度、壓力、LWC等結(jié)冰氣象參數(shù)，且經(jīng)對比證實(shí)其預(yù)測值與實(shí)際觀測值匹配良好?；赪RF的結(jié)冰氣象預(yù)測可為航空器結(jié)冰影響分析提供實(shí)時(shí)、準(zhǔn)確的結(jié)冰參數(shù)。

2） 基于POD降階模型和Kriging插值方法，所構(gòu)建的氣象參數(shù)與水滴收集質(zhì)量間的代理模型可較好地體現(xiàn)溫度、壓力、LWC、MVD等氣象參數(shù)對于水滴收集質(zhì)量的影響。雖然獲取采樣點(diǎn)處樣本數(shù)據(jù)會耗費(fèi)一定時(shí)間，但建立好代理模型后，可以快速準(zhǔn)確地預(yù)測目標(biāo)區(qū)域內(nèi)水滴收集質(zhì)量的空間分布與時(shí)序變化。

3） 基于PSO的結(jié)冰容限航跡規(guī)劃方法可在不同結(jié)冰安全約束條件下，規(guī)劃出無人機(jī)的最優(yōu)路徑。此外提出結(jié)冰容限概念，在粒子初始化階段，可增大搜索空間并增加解的多樣性，使得算法可求解出更為精確的無人機(jī)航跡。