不確定風(fēng)場下平流層浮空器全局路徑規(guī)劃

翟嘉琪，楊希祥，鄧小龍，龍遠(yuǎn)，張經(jīng)倫，柏方超

(國防科技大學(xué) 空天科學(xué)學(xué)院，長沙 410073)

平流層浮空器是指依靠輕于空氣的氣體提供升力，工作于平流層并可實現(xiàn)持久駐空的浮空類飛行器[1-2]，具有成本低、載重能力強(qiáng)、可長期駐空等優(yōu)點(diǎn)。作為在臨近空間開展應(yīng)用的主要平臺，平流層浮空器通過攜帶任務(wù)載荷，具備長時間、全方位、全天時和實時信息獲取能力，可以為高分辨率對地觀測、預(yù)警探測、通信中繼、防災(zāi)減災(zāi)、環(huán)境監(jiān)測等提供技術(shù)途徑[3]。

目前平流層浮空器的研究主要關(guān)注長時駐空、區(qū)域駐留等方向[4-8]。例如，Loon 氣球已經(jīng)完成了單球連續(xù)駐空335 天23 時的總飛行時間，并能提供網(wǎng)絡(luò)通信服務(wù)的飛行試驗[9]，且Loon 團(tuán)隊首次將強(qiáng)化學(xué)習(xí)應(yīng)用于氣球的飛控系統(tǒng)，并在太平洋赤道附近進(jìn)行了為期39 天，保持在定點(diǎn)距離50 km 以內(nèi)區(qū)域駐留的驗證試驗[10]；國內(nèi)多個機(jī)構(gòu)，如北京航空航天大學(xué)、中國科學(xué)院空天信息創(chuàng)新研究院、中國電子科技集團(tuán)第三十八研究所等研究的平流層浮空器也基本具備長航時駐空和區(qū)域駐留的能力。針對未來戰(zhàn)場的快速變化及平流層浮空器的實際應(yīng)用場景需求，如何實現(xiàn)從放飛點(diǎn)到目標(biāo)區(qū)域的快速部署，即浮空器的路徑規(guī)劃，將是制約浮空器未來大規(guī)模應(yīng)用的關(guān)鍵問題。

浮空器路徑規(guī)劃是指依靠風(fēng)場環(huán)境和浮空器設(shè)計參數(shù)，在約束條件下（如能源、動力、時間及距離等），規(guī)劃從起始點(diǎn)到目標(biāo)點(diǎn)的最優(yōu)路徑，不同約束條件可以規(guī)劃出不同的飛行路徑。

目前的路徑規(guī)劃大多是針對無人機(jī)、水下航行器、機(jī)器人及無人駕駛車輛等對象[11-14]，常見的路徑規(guī)劃方法主要分為2 類，即傳統(tǒng)路徑規(guī)劃和智能路徑規(guī)劃。其中，傳統(tǒng)路徑規(guī)劃中的路線圖構(gòu)建方法簡單、易于實現(xiàn)，主要能用于二維路徑規(guī)劃中；單元分解法能完整描述環(huán)境信息，但是規(guī)劃速度與分解的單元個數(shù)息息相關(guān)；這2 種方法都需要通過搜索方法來尋找最優(yōu)路徑，常見的有A*、D*方法；人工勢場法是一種虛擬算法，智能體在環(huán)境中受到來自目標(biāo)點(diǎn)的吸引力勢場和來自障礙物的排斥力勢場組成的合勢場來決定其運(yùn)動信息，該方法實時性好，計算簡單，但是容易出現(xiàn)局部鎖死、路徑振蕩等問題。針對傳統(tǒng)路徑規(guī)劃自適應(yīng)性較差，提出智能路徑規(guī)劃方法，因其源于模擬人的經(jīng)驗或生物的行為，具有自組織、自學(xué)習(xí)及一定的容錯能力，將其應(yīng)用于路徑規(guī)劃中，可以提高系統(tǒng)的自主能力，使得系統(tǒng)更加靈活，自主性和適應(yīng)性更強(qiáng)，主要包括遺傳算法、神經(jīng)網(wǎng)絡(luò)、強(qiáng)化學(xué)習(xí)及深度強(qiáng)化學(xué)習(xí)等。

目前針對平流層浮空器的路徑規(guī)劃研究較少，多為無人機(jī)等研究對象，通過對比，大尺寸、低動態(tài)的平流層浮空器在飛行過程中更易受風(fēng)場影響，因此，在浮空器的路徑規(guī)劃中可充分利用風(fēng)場提供推力，同時引入一定的水平弱動力，在風(fēng)場中可通過規(guī)劃動作序列，避免部分不利風(fēng)場影響，快速部署至目標(biāo)區(qū)域。

由于風(fēng)場環(huán)境的不確定性和時空復(fù)雜性，基于隨機(jī)概率的方法是研究浮空器的路徑規(guī)劃問題的合理手段，為此，研究人員針對馬爾可夫決策過程（Markov decision process, MDP）方法進(jìn)行了一定的研究。MDP 是一個描述系統(tǒng)與外部環(huán)境之間相互作用的模型，并嘗試選擇合理的動作使得累積報酬最大，不需要給定任何狀態(tài)下的監(jiān)督信號，可以用來處理不確定情況下的規(guī)劃問題，在路徑規(guī)劃中應(yīng)用廣泛。Wolf 等[15]提出一種基于MDP 的概率運(yùn)動規(guī)劃方法，其研究對象是運(yùn)行在土衛(wèi)六的Montgolfieré氣球（熱氣球的一種），仿真結(jié)果指出該方法可實現(xiàn)Montgolfieré氣球在土衛(wèi)六的全局路徑規(guī)劃。陳魁[16]基于MDP 模型設(shè)計了碼垛機(jī)器人的路徑規(guī)劃，采用最小二乘策略迭代方法設(shè)定基函數(shù)，求取最優(yōu)狀態(tài)-動作值函數(shù)，并且基于激光測距儀設(shè)計了基于部分可觀MDP（partially observable Markov decision process，POMDP）的局部路徑規(guī)劃模型，使機(jī)器人具有避障能力。Yu 等[17]針對無人機(jī)在惡劣環(huán)境下的避碰問題，提出一種基于MDP 的方法，規(guī)劃出來的無人機(jī)路徑安全無碰，仿真結(jié)果表明，設(shè)計出來的路徑可有效避免障礙物的攻擊，但是該方法并沒有將環(huán)境的不確定性引入其中，環(huán)境信息提前確定，不能很好地滿足實際情況。Kularatne等[18]基于MDP 計算在不確定海洋環(huán)境下的水下航行器路徑規(guī)劃問題，約束條件是最低能耗和有效避障，仿真結(jié)果表明：MDP 方法更適用于起始點(diǎn)和目標(biāo)點(diǎn)位置固定且不確定性已知的情況。Nanaz 等[19]通過對氣球的動力學(xué)模型進(jìn)行簡化和解耦，將全局路徑規(guī)劃問題轉(zhuǎn)化為圖形搜索問題，使用Dijkstra方法計算給定目標(biāo)點(diǎn)與土衛(wèi)六上任意位置的最短時間路徑，但是該方法不適用于不確定情況下的路徑規(guī)劃。

本文以具有弱動力的平流層浮空器為研究對象，提出一種基于MDP 的全局路徑規(guī)劃方法。所提方法充分考慮了風(fēng)場預(yù)測模型的誤差，并將誤差引起的風(fēng)場不確定性引入MDP 模型中，建立浮空器的二維全局路徑規(guī)劃模型。

1 問題描述及基本假設(shè)

由于平流層浮空器在長期駐空過程中的運(yùn)動多為二維平面運(yùn)動，本文對基于MDP 的浮空器全局路徑規(guī)劃方法有以下4 個假設(shè)。

1）不考慮浮空器垂直方向的運(yùn)動。

2）浮空器吊艙四周安裝4 個相互獨(dú)立的動力設(shè)備（電機(jī)和螺旋槳）為其提供水平驅(qū)動力，單個動力設(shè)備能提供的最大速度umax=5 m/s，浮空器最多允許2 個動力設(shè)備同時工作。

3）浮空器無動力時的速度與風(fēng)速一致，即處于隨風(fēng)飄飛的狀態(tài)。

4）不考慮浮空器的動力學(xué)特性，只考慮其運(yùn)動學(xué)特性。

本文研究區(qū)域在我國某地區(qū)，浮空器在該區(qū)域內(nèi)的位置用r?[x,y]表 示，其中x為 經(jīng)度，y為緯度，風(fēng) 速 用w(r,t)?[wx(r,t),wy(r,t)]表 示，其 中wx(r,t)和wy(r,t) 分 別表示在t時刻x軸 （東西方向）和y軸（南北方向）的風(fēng)速大小。在給定區(qū)域內(nèi)，以?x=0.3?,?y=0.3?為精度進(jìn)行離散化，得到如圖1所示的21×21柵 格區(qū)域。針對柵格i，si表示浮空器的狀態(tài)信息，ri表示浮空器在狀態(tài)si下 的位置，w(ri,t)表示風(fēng)場信息。

圖1 區(qū)域離散化示意圖Fig.1 Diagram of regional discretization

在風(fēng)場的作用下，浮空器的速度v(si)=w(ri)+u(si)，其中u(si)表 示浮空器在狀態(tài)si下自身所能提供的速度。針對柵格i，si表 示當(dāng)前狀態(tài)，sj表示下一狀態(tài)，風(fēng)場信息從狀態(tài)si到 達(dá)臨近狀態(tài)sj的過程中始終保持為w(ri)，到達(dá)狀態(tài)sj時才更新為w(rj)。在目標(biāo)點(diǎn)確定時，位于給定區(qū)域內(nèi)任意位置的浮空器在風(fēng)場和自身動力的共同作用下向目標(biāo)點(diǎn)靠近，但是有的位置在共同作用下始終無法抵達(dá)目標(biāo)點(diǎn)，而某些位置可以在很短的時間內(nèi)到達(dá)目標(biāo)點(diǎn)，如圖1 所示，即本文重點(diǎn)研究的區(qū)域可達(dá)性和最短時間的最優(yōu)路徑問題。

定義浮空器的移動方向為正北、東北、正東、東南、正南、西南、正西、西北8 個方向，將360°平均分為8 個部分，如圖2 所示，其中，θi j為由狀態(tài)si進(jìn) 入相鄰狀態(tài)sj將空間分成8 個部分的較小邊界角度。

圖2 浮空器狀態(tài)轉(zhuǎn)移方向Fig.2 Transition direction of a stratospheric aerostat

2 相關(guān)方法模型

2.1 馬爾可夫決策過程

MDP 是一種描述智能體和外部環(huán)境之間相互作用的模型，如圖3 所示，常被用來作為序列決策問題的框架，可以用一個五元組M=(S,A,P,R,γ)來表示[14]，S代表環(huán)境狀態(tài)的集合，A代表智能體能執(zhí)行的基本動作集合，Pa(s,s′)表 示智能體在狀態(tài)s下執(zhí)行動作a到 達(dá)狀態(tài)s′的 概率，Ra(s,s′)表示智能體在狀態(tài)s下執(zhí)行動作a到 達(dá)狀態(tài)s′所得到的即時獎勵，γ表示折扣因子，取值范圍為 [0,1]，γ越大表示后續(xù)獎勵對當(dāng)前獎勵的影響越大。

圖3 MDP 模型Fig.3 Model of Markov decision processes

當(dāng)MDP 用來表示序列決策問題時，決策問題的解即為策略（π），是從狀態(tài)集合到動作集合的映射，即 π:S→A，求解MDP 即可表示為求解最優(yōu)策略（π?）的問題，即浮空器在所處狀態(tài)采取的最優(yōu)動作，使得所獲得的累積獎勵最大。累積獎勵可形式化表示為浮空器所處狀態(tài)的價值（狀態(tài)-值函數(shù)）：

式中：Vπ(s) 為 執(zhí)行策略 π后 在狀態(tài)s下 的價值；t為決策時刻；Rt為t時刻下的即時獎勵。

對于任意狀態(tài)，使得值函數(shù)最優(yōu)的充要條件是Bellman 方程，最優(yōu)策略 π?是使得值函數(shù)取得最大值時的策略：

式中：π?(s)為 智能體在狀態(tài)s下的最優(yōu)策略。

求解最優(yōu)策略常見的方法有策略迭代和值迭代[20]，其中策略迭代需在每個狀態(tài)-值收斂后再進(jìn)行策略提升，不斷地迭代直到策略收斂，而值迭代在狀態(tài)-值每更新一次之后進(jìn)行策略提升，再重復(fù)進(jìn)行狀態(tài)-值更新和策略提升直至狀態(tài)-值收斂[21]。因值迭代方法在本文研究的問題中具有更高的計算效率，本文采用值迭代方法進(jìn)行求解。

2.2 二維風(fēng)場分布

在實際情況下，風(fēng)場數(shù)據(jù)為該區(qū)域內(nèi)的各個離散氣象站提供的預(yù)測數(shù)據(jù)，在考慮風(fēng)場因素時直接把某個氣象站的數(shù)據(jù)作為該區(qū)域內(nèi)的整體風(fēng)場數(shù)據(jù)，風(fēng)向、風(fēng)速只隨高度變化，不會對該區(qū)內(nèi)各個高度的二維風(fēng)場分布進(jìn)行處理，但是在實際路徑規(guī)劃中，必須考慮二維風(fēng)場的空間分布。本文根據(jù)各個氣象站的風(fēng)場預(yù)測數(shù)據(jù)，采用雙線性差值法如圖4 所示，根據(jù)各點(diǎn)距離氣象站的距離確定權(quán)值，可以得到該區(qū)域內(nèi)各個高度的二維風(fēng)場分布，具體步驟如下：

圖4 雙線性差值法示意圖Fig.4 Diagram of bilinear interpolation method

1）二維線性差值。按照雙線性差值法如式（4）所示，可以得到二維區(qū)域內(nèi)任意位置的數(shù)值。

式中：（x1,y1）、（x1,y2）、（x2,y1）、（x2,y2）、（x,y）分別為Q11、Q12、Q21、Q22、P的位置坐標(biāo)。

2）確定權(quán)值。由于區(qū)域中各個點(diǎn)的風(fēng)速風(fēng)向與每個氣象站的實際數(shù)據(jù)有關(guān)，即目標(biāo)點(diǎn)距離某個氣象站越近，該點(diǎn)的風(fēng)速風(fēng)向越接近該氣象站的實際數(shù)據(jù)，如圖5 所示，因此，目標(biāo)點(diǎn)的氣象數(shù)據(jù)與各個氣象站在該點(diǎn)的權(quán)值有關(guān)，將其設(shè)置為與氣象站的距離成反比，為目標(biāo)點(diǎn)與氣象站的水平距離，權(quán)值按照式（5）取值。

圖5 二維平面權(quán)值確定示意圖Fig.5 Diagram of determination of two-dimensional plane weights

式中：η為 權(quán)值；(xk,yk)為 待求目標(biāo)的位置坐標(biāo)；d為待求目標(biāo)與已知點(diǎn)之間的距離；N為氣象站點(diǎn)數(shù)量。

3）二維風(fēng)場可視化。二維風(fēng)場速度的大小通過箭頭的長度和顏色表示，箭頭越短，顏色越藍(lán)，風(fēng)速越?。伙L(fēng)向通過箭頭的指向來表示，箭頭的方向即為風(fēng)的來向。為了和浮空器的狀態(tài)數(shù)量保持一致，二維風(fēng)場分布也劃分為 21×21個柵格。如圖6所示。

圖6 不同高度的二維風(fēng)場分布對比Fig.6 Comparison of two-dimensional wind field distribution at different heights

2.3 不確定風(fēng)場模型

由于風(fēng)場模型可能存在的不準(zhǔn)確性，風(fēng)場數(shù)據(jù)存在一定的誤差，為了使得路徑規(guī)劃更加符合實際情況，本文將不確定性引入風(fēng)場中，建立不確定風(fēng)場模型。風(fēng)場信息包括風(fēng)向和風(fēng)速，分別用 θi和wi表示，兩者具有相互獨(dú)立的概率分布，風(fēng)速和風(fēng)向由w(ri)表示[15]。

1）風(fēng)向概率分布。風(fēng)向概率分布采用von Mises 分布，該分布是圓上的連續(xù)概率分布，近似于包裹正態(tài)分布，是正態(tài)分布的圓形模擬，概率密度函數(shù)為

式中：θi=∠w(ri)為可獲得的風(fēng)場數(shù)據(jù)中的實際風(fēng)向，整個分布圍繞它展開，取值范圍為 [0?,360?]，其中正北風(fēng)向為 0?（3 60?），正東方向為 90?，以此類推可以得到所有風(fēng)向?qū)?yīng)的角度；κ為風(fēng)向分布的集中程度，κ＞0，且 κ →0，分布越分散，κ →∞，分布越集中；I0(κ)為0 階Bessel 修正函數(shù)。

2）風(fēng)速概率分布。風(fēng)速概率分布采用Gaussian分布，其概率密度函數(shù)為

式中：σi=ρwi，wi=//w(ri)//，w(ri)為可獲得的風(fēng)場數(shù)據(jù)中的實際風(fēng)速，整個風(fēng)速分布圍繞它展開；ρ為分布幅度的參數(shù)，可按實際情況取值。

根據(jù)實際風(fēng)場數(shù)據(jù)，在19 000 m 高度、37°經(jīng)度、86°緯度處，風(fēng)向 θi=251.916?，κ=1/400，風(fēng)速wi=[0.344 6,?3.559]，ρ=0.2，可以得到該點(diǎn)處的風(fēng)向、風(fēng)速的概率分布，如圖7 所示。

圖7 不確定模型下風(fēng)向和風(fēng)速概率分布Fig.7 Probability distribution of wind direction and speed under uncertain model

3 基于MDP 的全局路徑規(guī)劃

3.1 MDP 模型中關(guān)鍵參數(shù)設(shè)計

1）環(huán)境狀態(tài)集合S。本文浮空器所處區(qū)域，經(jīng)過離散化得到 21×21個柵格，MDP 模型中的狀態(tài)S即為 21×21個柵格中的風(fēng)速和風(fēng)向。由于風(fēng)向的范圍為 [0?,360?]，按照圖3 中對角度和方向的劃分，風(fēng)向可分為如下8個方向，分別為：①正北：[0?,22.5?]∪[337.5?,360?]；②東北：[22.5?,67.5?]；③正東：[67.5?,112.5?]；④東 南：[112.5?,157.5?]；⑤正 南：[157.5?,202.5?]；⑥西 南：[202.5?,247.5?]；⑦正 西：[247.5?,292.5?]；⑧西北：[292.5?,337.5?]；

當(dāng)風(fēng)場中的風(fēng)向數(shù)值確定之后，浮空器在無動力情況下的移動方向就確定為這8 個方向之一。

2）基本動作集合A。浮空器自身的速度由4 個相互獨(dú)立的動力設(shè)備提供，且最多允許2 個設(shè)備同時工作，通過控制動力設(shè)備的開啟和關(guān)閉，可以得到9 個基本動作，正東方向為x軸正方向，正北方向為y軸正方向，浮空器基本動作方向如圖8 所示。動力設(shè)備均不工作，浮空器無動力：[0,0]m/s。

圖8 浮空器基本動作方向Fig.8 Basic action direction of stratospheric aerostat

只有1 個動力設(shè)備工作，浮空器分別具有向正北、正東、正南、正西方向的速度，通過調(diào)整動力設(shè)備使得浮空器可具備5 m/s 和2 m/s 的速度，速度大小為：[0,5]，[5,0]，[0,?5]，[?5,0]，[0,2]，[2,0]，[0,?2]，[?2,0]m/s。

2 個動力設(shè)備同時工作（本文假設(shè)2 個設(shè)備所能提供的速度一致），浮空器分別具有向東北、東南、西南、西北的速度，通過調(diào)整動力設(shè)備使得浮空器在單個方向可具備5 m/s 和2 m/s 的速度，此時速度大小為：[5,5]，[5,?5][?5,?5]，[?5,5]，[2,2]，[2,?2]，[?2,?2]，[?2,2] m/s。

3）轉(zhuǎn) 移 概率P。Pa(si,sj)表 示浮 空 器 在狀 態(tài)si下 通過采取某個動作a進(jìn) 入到狀態(tài)sj的概率。轉(zhuǎn)移過程中包括2 種模型：一種是確定模型，即浮空器在風(fēng)速和自身驅(qū)動速度的共同作用下，從初始狀態(tài)si轉(zhuǎn) 移到的唯一確定的臨近狀態(tài)sj；另一種是不確定模型，即由于風(fēng)場的不確定性，導(dǎo)致浮空器在不確定的風(fēng)速和自身驅(qū)動速度的共同作用下，從初始狀態(tài)si轉(zhuǎn) 移到的臨近狀態(tài)sj是不確定的，是具有一定的概率的。在浮空器無動力的情況下，風(fēng)向分布決定浮空器的轉(zhuǎn)移概率為

浮空器在無動力、速度為 [5,5]m /s 及 [?5,?5]m/s情況下的轉(zhuǎn)移概率分布情況如圖9 所示。

圖9 不確定風(fēng)場下浮空器轉(zhuǎn)移概率分布Fig.9 Transition probability distribution of stratospheric aerostate under uncertain model

4）獎勵函數(shù)R。為描述浮空器狀態(tài)、動作和目標(biāo)之間的關(guān)系，獎勵函數(shù)用Ra(si,sj)表示，具體含義為浮空器在選擇動作a的情況下由狀態(tài)si轉(zhuǎn)移到狀態(tài)sj所 消耗時間的負(fù)值，獎勵函數(shù)R為即時獎勵，且共有 21×21×9個 數(shù)值。浮空器在由狀態(tài)si移動到sj過 程中所消耗時間為 δt，則Ra(si,sj)=?δt，如下：

式中：dist(·,·)為 2 個柵格間的實際距離；vi,a為浮空器實際速度矢量，由風(fēng)速和浮空器自身速度合成，vi,a=w(ri)+ua；w(ri)為 風(fēng)場在位置ri處的矢量風(fēng)速；u(si)為 浮空器在狀態(tài)si下所能提供的速度。

路徑規(guī)劃目的是使浮空器在給定區(qū)域內(nèi)找到一條從起始點(diǎn)到目標(biāo)點(diǎn)飛行時間最短路徑，因此，獎勵函數(shù)需考慮到距離區(qū)域邊界、遠(yuǎn)離目標(biāo)點(diǎn)、趨向目標(biāo)點(diǎn)等所有情況，在確定風(fēng)場下的獎勵函數(shù)為

式中：Ra(si,sj)表示浮空器分別在轉(zhuǎn)移過程中趨向目標(biāo)點(diǎn)、遠(yuǎn)離目標(biāo)點(diǎn)及超出區(qū)域邊界時的即時獎勵，式（14）情況相同。

在不確定風(fēng)場下，由于風(fēng)場信息的不確定導(dǎo)致浮空器的實際速度也無法確定，本文用均值來表示浮空器在狀態(tài)si下 通過采取動作a得到的平均速度，此時相鄰2 個狀態(tài)的飛行時間為

由于風(fēng)場的不確定性，導(dǎo)致浮空器由狀態(tài)si進(jìn)入到相鄰狀態(tài)的sj也是無法確定的，通過不確定風(fēng)場的轉(zhuǎn)移概率，不確定風(fēng)場下的即時獎勵根據(jù)飛行時間期望值來確定，而飛行時間的期望值則根據(jù)轉(zhuǎn)移概率計算：

3.2 值迭代方法

本文中的風(fēng)場模型包括確定風(fēng)場和不確定風(fēng)場2 種模型，因此，值迭代方法也分為2 種，在確定風(fēng)場下，浮空器選擇動作后到達(dá)下一狀態(tài)是確定的；而在不確定風(fēng)場下，浮空器選擇動作后到達(dá)下一狀態(tài)具有一定的概率，區(qū)別在于狀態(tài)轉(zhuǎn)移概率Pa(si,sj)如何求解，值迭代方法步驟如下：

步驟1輸入環(huán)境狀態(tài)、基本動作集合、即時獎勵、轉(zhuǎn)移概率、目標(biāo)點(diǎn)位置，并初始化狀態(tài)值函數(shù)。

步驟2計算所有環(huán)境狀態(tài)下確定風(fēng)場和不確定風(fēng)場下的狀態(tài)值函數(shù)。

步驟3計算2 次狀態(tài)值函數(shù)的差值。

步驟4若差值達(dá)到值迭代收斂閾值，則末次計算的值函數(shù)即為最優(yōu)狀態(tài)值函數(shù)，若未達(dá)到，則返回步驟2，迭代次數(shù)增加1 次。

步驟5若迭代次數(shù)大于最大迭代次數(shù)，則末次計算的值函數(shù)即為最優(yōu)狀態(tài)值函數(shù)，若未達(dá)到，則返回步驟2，迭代次數(shù)增加1 次。

3.3 路徑規(guī)劃流程

全局路徑規(guī)劃流程見圖10。實現(xiàn)步驟如下：

圖10 全局路徑規(guī)劃流程Fig.10 Flowchart of global path planning

步驟1初始化參數(shù)。

步驟2以風(fēng)場預(yù)測數(shù)據(jù)為輸入，建立二維風(fēng)場模型。

步驟3針對風(fēng)場不確定性，建立二維不確定風(fēng)場模型。

步驟4以狀態(tài)集合、動作集合為輸入，計算狀態(tài)轉(zhuǎn)移概率。

步驟5以狀態(tài)集合、動作集合、目標(biāo)點(diǎn)位置為輸入，計算2 種風(fēng)場模型下的即時獎勵。

步驟6根據(jù)浮空器的即時獎勵和轉(zhuǎn)移概率，用值迭代方法計算浮空器的最優(yōu)狀態(tài)值函數(shù)，并輸出最優(yōu)狀態(tài)值函數(shù)和最優(yōu)策略。

步驟7以飛行區(qū)域內(nèi)任意位置為起始點(diǎn)，根據(jù)最優(yōu)狀態(tài)值函數(shù)，獲得每次轉(zhuǎn)移的最優(yōu)動作，結(jié)合風(fēng)場得到轉(zhuǎn)移方向。

步驟8判斷是否到達(dá)目標(biāo)點(diǎn)，若到達(dá)，則結(jié)束流程，若未到達(dá)，則返回步驟8。

4 仿真結(jié)果分析

將第3 節(jié)中的模型和方法應(yīng)用于浮空器的全局路徑規(guī)劃，通過在MATLAB2020a 環(huán)境下建立的模型進(jìn)行仿真，驗證參數(shù)設(shè)計和方法流程的有效性和工程實用性。

4.1 區(qū)域可達(dá)性仿真結(jié)果分析

在實際任務(wù)中，當(dāng)給定目標(biāo)點(diǎn)時，浮空器需要快速機(jī)動部署至目標(biāo)區(qū)域?qū)嵤┤蝿?wù)，但是由于浮空器自身載荷能力、動力設(shè)備等限制，導(dǎo)致浮空器的驅(qū)動速度不高，抗風(fēng)能力較弱。在沒有動力甚至在有一定動力的情況下，若浮空器自身驅(qū)動能力無法補(bǔ)償風(fēng)場的不利影響，也會導(dǎo)致浮空器在某些位置無論采取何種動作均無法達(dá)到目標(biāo)點(diǎn)。

本節(jié)主要對浮空器自身驅(qū)動速度與給定目標(biāo)的可達(dá)性之間的關(guān)系進(jìn)行分析，使用的是中國某地區(qū)的6 個氣象站點(diǎn)的氣象數(shù)據(jù)，按照第2 節(jié)中的方法得到19 000 m 高度的二維風(fēng)場分布，給定目標(biāo)點(diǎn)為[88.7?,41.2?]，圖11 中白色方框所在位置，圖11 和圖12 分別為在確定風(fēng)場和不確定風(fēng)場的情況下，在無動力、單個動力設(shè)備提供的最大速度umax=2 m/s、umax=5 m/s時，浮空器在區(qū)域內(nèi)任意位置到達(dá)給定目標(biāo)的期望時間。

圖11 確定風(fēng)場下不同最大抗風(fēng)能力期望到達(dá)時間分布Fig.11 Distribution of expected arrival time at different speeds in certain wind field

圖12 不確定風(fēng)場下不同速度期望到達(dá)時間分布Fig.12 Distribution of expected arrival time at different speeds in uncertain wind field

1）確定風(fēng)場模型。在確定風(fēng)場模型下，浮空器在每個狀態(tài)下的移動方向由該狀態(tài)下的風(fēng)速、風(fēng)向及浮空器自身提供的速度唯一確定。圖11 和圖12中柵格的顏色表示到達(dá)目標(biāo)點(diǎn)的期望時間，顏色越紅到達(dá)時間越長，根據(jù)19 000 m 高度處的風(fēng)場分布可知以北風(fēng)分量為主，西風(fēng)分量較小。當(dāng)浮空器無動力時只能通過風(fēng)場的作用到達(dá)目標(biāo)點(diǎn)，只有在特殊幾個柵格狀態(tài)下可以快速部署到目標(biāo)區(qū)域；當(dāng)浮空器的單個動力設(shè)備提供最大速度為2 m/s 時，浮空器在某些位置下期望到達(dá)時間為100 h，當(dāng)某些位置的期望到達(dá)時間過大時，在實際情況下可將其定義為不可到達(dá)點(diǎn)，即起始點(diǎn)不可選在這些位置；當(dāng)浮空器的單個動力設(shè)備提供最大速度為5 m/s時，浮空器在給定區(qū)域內(nèi)的大部分位置到達(dá)時間在10 h 以內(nèi)，只有極少數(shù)的狀態(tài)下期望到達(dá)時間為14 h，該區(qū)域內(nèi)的任意點(diǎn)均可作為起始點(diǎn)。

2）不確定風(fēng)場模型。在不確定風(fēng)場模型下，浮空器的實際速度和實際轉(zhuǎn)移方向是按照一定的概率進(jìn)行分布的，通過Mont Carlo 方法統(tǒng)計得到風(fēng)場的分布概率，由圖12 可知不確定風(fēng)場與確定風(fēng)場相比，期望時間的分布規(guī)律大致相同，風(fēng)場按照第2 節(jié)中的概率分布，導(dǎo)致浮空器在轉(zhuǎn)移到臨近狀態(tài)時存在不確定性，根據(jù)轉(zhuǎn)移概率和即時獎勵（相鄰2 個狀態(tài)之間的飛行時間）來確定狀態(tài)值函數(shù)。隨著浮空器自身提供的驅(qū)動速度增加，給定區(qū)域內(nèi)的各個位置可到達(dá)給定目標(biāo)點(diǎn)的數(shù)量增加，因此，浮空器的水平驅(qū)動速度可改變區(qū)域內(nèi)各位置相對于給定目標(biāo)點(diǎn)的可達(dá)性。

3）2 種風(fēng)場模型的對比。由圖11 和圖12 可知，當(dāng)浮空器無動力時，對于給定目標(biāo)點(diǎn)，不確定風(fēng)場模型下較為合理期望到達(dá)時間占據(jù)柵格更多（藍(lán)色柵格），確定風(fēng)場模型下不可達(dá)的柵格數(shù)量更多（紅色柵格）；但是隨著浮空器水平驅(qū)動能力的加強(qiáng)，確定風(fēng)場模型下藍(lán)色柵格的分布和數(shù)量與不確定風(fēng)場模型下藍(lán)色柵格的分布和數(shù)量逐漸一致，說明浮空器的水平驅(qū)動能力可以抵消一部分風(fēng)場不確定的影響。

4.2 最優(yōu)路徑規(guī)劃仿真結(jié)果分析

本文規(guī)劃的路徑為最短時間路徑，浮空器在到達(dá)目標(biāo)點(diǎn)的過程中需不斷根據(jù)實際風(fēng)場來選擇最優(yōu)動作序列，使其能夠在最短的時間內(nèi)快速機(jī)動至目標(biāo)區(qū)域。通過對浮空器在區(qū)域內(nèi)各個位置的期望到達(dá)時間的分析，可以找到狀態(tài)轉(zhuǎn)移的最優(yōu)方向，去除風(fēng)場的作用，即可得到浮空器在不同驅(qū)動速度下的最優(yōu)動作序列，浮空器在實際飛行中可按照最優(yōu)動作序列控制執(zhí)行機(jī)構(gòu)作動，為浮空器的飛行決策提供指導(dǎo)。

圖13 和圖14 分別為確定風(fēng)場和不確定風(fēng)場下浮空器在不同起始點(diǎn)、不同速度下的最優(yōu)動作序列和規(guī)劃出來的路徑。目標(biāo)點(diǎn)的經(jīng)緯度坐標(biāo)設(shè)置 為 [88.7?,41.2?]，起 始 點(diǎn) 經(jīng) 緯 度 坐 標(biāo) 設(shè) 置 為：[86.3?,42.4?]、[90.8?,37.3?]，浮空器單個動力設(shè)備提供的最大速 度umax=2 m/s 、umax=5 m/s。圖13 和圖14 中的黑色箭頭為該位置下浮空器的最優(yōu)速度方向，紅色箭頭為風(fēng)向，在不確定風(fēng)場中，紅色箭頭的長短表示風(fēng)場分布在該方向的概率大小。

圖13 確定風(fēng)場下路徑規(guī)劃和最優(yōu)動作序列Fig.13 Path planning and optimal velocity sequence in certain wind field

圖14 不確定風(fēng)場下路徑規(guī)劃和最優(yōu)動作序列Fig.14 Path planning and optimal velocity sequence in uncertain wind field

浮空器自身的驅(qū)動能力、起始點(diǎn)的位置及風(fēng)場模型均會影響浮空器最優(yōu)策略即最優(yōu)動作序列的選擇。圖13（a）和圖14（a）中，當(dāng)浮空器單個動力設(shè)備能提供的最大速度為2 m/s 時，有的位置沒有黑色箭頭，說明在這些位置無論浮空器采取何種動作，都無法到達(dá)靠近目標(biāo)點(diǎn)，因此，在浮空器實際飛行過程中應(yīng)該避免到達(dá)這些位置；但是當(dāng)浮空器單個動力設(shè)備能提供的最大速度為5 m/s 時，浮空器能夠在比較短的時間內(nèi)到達(dá)目標(biāo)點(diǎn)，相同位置下，由于浮空器自身速度的改變導(dǎo)致浮空器在整個區(qū)域的可達(dá)性也發(fā)生改變；對比圖13（b）、圖13（c）和圖14（b）、圖14（c）針對同一起始點(diǎn)和目標(biāo)點(diǎn)規(guī)劃的路徑，在不確定風(fēng)場下的路徑轉(zhuǎn)彎的次數(shù)更少，由于風(fēng)場的概率分布，浮空器的實際速度也存在概率分布，導(dǎo)致在不確定風(fēng)場下浮空器的路徑規(guī)劃更符合實際情況。

本文所提方法是在對整個風(fēng)場環(huán)境和浮空器狀態(tài)空間進(jìn)行有效的離散化的基礎(chǔ)上開展的，與基于圖搜索的路徑規(guī)劃方法有一定的相似之處。Dijkstra 方法是其中比較常用的路徑規(guī)劃方法，主要通過搜索空間來求解最優(yōu)路徑，因此，搜索效率不高，尤其是在不確定風(fēng)場下的路徑規(guī)劃問題，不適用于本文的研究背景。

5 結(jié) 論

1）本文所提方法能夠在不確定二維風(fēng)場環(huán)境中規(guī)劃出一條針對目標(biāo)點(diǎn)的最短時間路徑，并能給出給定區(qū)域內(nèi)各個位置下的最優(yōu)動作序列和期望達(dá)到時間，為浮空器的實際飛行執(zhí)行機(jī)構(gòu)作動策略提供理論依據(jù)。

2）本文所提方法還可提供浮空器在不同水平驅(qū)動能力下針對目標(biāo)點(diǎn)的區(qū)域可達(dá)性，為后續(xù)浮空器的飛行策略提供指導(dǎo)。