復雜氣象條件下考慮結冰風險的無人機飛行策略

郭琪磊，桑為民，牛俊杰，袁燁

1．西北工業大學 航空學院，西安 710072

2．中國民用航空飛行學院 工程技術訓練中心，廣漢 618307

3．中國空氣動力研究與發展中心 結冰與防除冰重點實驗室，綿陽 621000

4．中國商用飛機有限責任公司 上海飛機設計研究院，上海 201210

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

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

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

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

1 理論方法與數值模型

1.1 WRF模式配置

本文中結冰氣象預測采用WRF模式，版本為4.1.2。WRF模式是由美國國家環境預測中心（National Centers for Environmental Prediction，NCEP）和美國國家大氣研究中心等聯合開發的中尺度天氣數值預報系統。該模式采用可壓非靜力動力框架和地形追隨坐標系，水平方向采用Arakawa-C型網格劃分計算域，具有多重網格嵌套功能和豐富的大氣物理過程參數化方案。結合各國氣象機構提供的全球氣象數據和資料同化技術，WRF模式廣泛應用于不同時間/區域尺度天氣現象模擬、空氣質量建模、大氣海洋和理想化模擬等領域［16］，用于結冰氣象分析也有優異表現［17-18］。

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

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

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

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

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

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

1.2 基于POD和Kriging的代理模型

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

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

其中：( · ， · )表示內積運算。

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

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

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

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

給定n個輸入參數樣本點及其相對應的系統響應值其中X(i)為m維向量。Kriging模型假設目標函數值與設計變量之間的真實關系可以表達為［36］

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

其中：Z(X(i)，X(j))為輸入參數變量X(i)和X(j)的相關函數，其與兩點間空間位置密切相關，并隨著空間距離增大而減小。

式中：θq為第q個輸入參數分量xq的距離權值。

構建Kriging代理模型的關鍵在于尋找最佳的θ向量，通過最大似然法將該問題轉化為如下的帶約束優化問題：

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

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

在求得最優θ向量后，對任意一個新的輸入參量向量X(0)，對應的Kriging模型的估計值為

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

1.3 基于PSO的結冰容限航跡規劃方法

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

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

式中：w為粒子的慣性權重系數；c1和c2為正實數，分別表示粒子的認知和社會加速度系數；r1和r2為分布于［0，1］區間的隨機數。

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

圖1 所構建的二維航跡規劃空間Fig. 1 Constructed 2D path planning workspace

本文中，所規劃航跡需滿足以下條件：① 所規劃航跡航程最短；② 所規劃航跡需滿足指定的結冰風險。因此單機航跡規劃問題的數學模型表示為

式中：ωi為任意待選航跡ph中的導航節點；JL為任意待選航跡的航程長度；Lmax為無人機的最大飛行航程；dis(ωi，ωi+1)為待選航跡ph中導航節點ωi與ωi+1之間的歐式距離；workplace為整個規劃空間；semi-free workplace為規劃空間中滿足結冰風險的半自由規劃空間。

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

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

基于PSO的結冰容限航跡規劃算法流程如算法1所示。

2 數值模型驗證

2.1 基于WRF模式的結冰氣象預測

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

算法1 基于PSO的結冰容限航跡規劃Algorithm 1 Ice tolerance route planning based on PSO

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

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

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

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

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

圖4為LWC時序變化與文獻［41］中結果對比。圖中黑色圓點為實際觀測值，藍色線為本文中WRF模擬結果，紅色線為文獻［41］中模擬結果。可知本文中數值模擬結果雖然在數值精度上與實際觀測值存在偏差，但能夠較好反映出實際觀測值的時序變化。此外本文和文獻［41］中LWC的MAE分別0.084 g/m3和0.099 g/m3，本文預測LWC的表現優于文獻［41］中結果。因此本文所構建WRF模式參數化方案組合能夠準確預測出云層中的液態水含量。

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

2.2 基于代理模型的水滴收集質量快速預測

根據1.2節中所述的代理模型，構建水滴收集質量快速預測方法，相應的壁面水滴收集質量計算式為

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

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

圖5 連續最大結冰條件下采樣點分布示意圖Fig. 5 Schematic diagram of sampling point distribution under continuous maximum icing conditions

表3 部分采樣點工況Table 3 Working conditions of partial sampling points

圖6為表3中采樣點工況下水滴收集質量在翼面處分布對比，其中Y為與弦線垂直的方向，翼型前緣點為坐標原點。可知在翼型前緣處水滴收集質量最大：受來流角度影響，水滴收集質量的最大值出現在翼型前緣下沿。而在前緣點之后沿著翼型表面水滴收集質量逐漸減小，相應地下翼面處水滴收集質量變化相較于上翼面更為緩慢，水滴撞擊極限也相對更大。水滴收集質量峰值與撞擊極限分布受LWC、MVD、溫度T、高度H、來流速度等因素綜合影響，其中水滴收集質量極值分布受LWC和來流速度的影響較大，如Case 1與Case 4。而由于直徑較大、速度較大的液滴具有較大的慣性，在撞擊機翼表面時液滴軌跡不易發生改變。因此水滴撞擊極限范圍通常隨MVD和來流速度增大而增加，如Case 4與Case 5。

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

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

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

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

3 結果與討論

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

無人機結冰主要受氣象因素（如大氣溫度、飛行高度、LWC、MVD）以及飛行狀態（如速度）等因素影響。基于WRF模式的結冰氣象預測可獲得無人機在某一飛行環境中的大氣溫度、氣壓與LWC等參數。需要說明的是在本文中，當大氣溫度<273.15 K時，則認為液態水為過冷態，過冷水溫度與大氣溫度保持一致。根據圖5所示連續最大結冰條件與預測所得的LWC，進行插值則可以獲得目標區域內MVD。圖9~圖11分別為UTC 2021-05-26T21：00、UTC 2021-06-05T01：00以及UTC 2021-07-11T08：00在高度層eta=25處目標區域內壓力、溫度與LWC分布云圖。

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

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

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

以上述3組結冰氣象條件為基礎，假定無人機飛行速度為50 m/s，利用所構建的基于POD與Kriging的水滴收集質量代理模型，分別對目標區域內水滴收集質量進行快速預測，獲得水滴收集質量空間分布如圖12所示。

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

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

表5 飛機結冰強度等級劃分Table 5 Aircraft icing intensity classification

以圖12中目標區域內3組水滴收集質量分布作為無人機航跡規劃空間。由于3組航跡規劃空間內最大水滴收集質量均未超過0.02 kg/（s?m2），意味著最大結冰強度等級為中度結冰，因此分別選取表5中微量結冰與輕度結冰等級上極限，即水 滴 收 集 質 量 為0.009 2 kg/（s?m2）與0.015 3 kg/（s?m2）作為無人機航跡規劃中結冰安全約束。當所規劃航跡點水滴收集質量小于上述上極限時，認為該航跡點滿足結冰安全約束；否則違反結冰安全約束，該航跡點應該規避。

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

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

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

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

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

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

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

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

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

4 結 論

為解決無人機在復雜氣象條件下易受結冰影響，從而威脅其飛行安全的問題，提出了一種考慮結冰風險的無人機航跡規劃方法。基于中尺度天氣預報WRF模式對海南樂東地區結冰氣象環境進行預測，通過參數化方案敏感性分析確定了最佳參數化方案組合，進而獲得該時段內目標區域的溫度、壓力、LWC空間分布及時序變化。與此同時，利用OLHS方法對FAR 25部附錄C中連續最大結冰條件進行采樣，并對40個采樣點進行水滴撞擊特性計算，獲得各采樣點處水滴收集質量分布；基于POD降階模型和Kriging插值方法，構建溫度、壓力、LWC、MVD等氣象參數與水滴收集質量間的代理模型。基于結冰氣象預測參數與代理模型，獲得海南樂東地區水滴收集質量的空間分布與時序變化。最后，分別以不同結冰強度下水滴收集質量閾值作為結冰安全約束，利用基于PSO的結冰容限航跡規劃方法對無人機進行考慮結冰風險的飛行策略研究。研究結果表明：

1） 利用WRF模式可獲得目標區域內溫度、壓力、LWC等結冰氣象參數，且經對比證實其預測值與實際觀測值匹配良好。基于WRF的結冰氣象預測可為航空器結冰影響分析提供實時、準確的結冰參數。

2） 基于POD降階模型和Kriging插值方法，所構建的氣象參數與水滴收集質量間的代理模型可較好地體現溫度、壓力、LWC、MVD等氣象參數對于水滴收集質量的影響。雖然獲取采樣點處樣本數據會耗費一定時間，但建立好代理模型后，可以快速準確地預測目標區域內水滴收集質量的空間分布與時序變化。

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