內河船舶尾氣監測的多無人機路徑規劃研究

王廣生，孫祎崢，孫海軍，魚 童，張 鋮，周云鵬

（1.中國港灣工程有限責任公司，北京 100027；2.河海大學港口海岸與近海工程學院，江蘇南京 210098）

引言

近年來，盡管國際海事組織（IMO）及各國政府通過降低船用燃料油中的硫含量來控制船舶大氣污染物排放，但由于高、低硫油價格存在較大差距，低硫油高昂的使用成本使得高硫燃油使用的違規率依然較高。為解決此問題，已有眾多研究嘗試通過移動式監測系統[1]、岸基遙感監測[2]、無人機嗅探檢測[3]等方式開展船舶尾氣污染防治工作。其中，無人機嗅探檢測方式憑借在機動性強和執法效率高等優勢，受到了行業內學者的廣泛關注，在海事管理領域的可行性得到了充分的證明[4-6]。

然而，在船舶流量密度較大的情況下，僅依靠人工操控無人機完成檢測工作效率低、危險性也較高，因此研究續航時間有限的條件下，多無人機協同是實現無人機自動化監測的關鍵問題。

多無人機協同類屬于任務分配問題[7]，即在已知無人機數量和任務數量的條件下，考慮到環境威脅、無人機的續航能力、載荷能力約束以及任務的特殊需求，進行合理的任務分配，使得各無人機均能夠以最優效率協同完成工作。路徑尋優問題一般可將其歸為經典的旅行商問題[8]，多無人機協同的情景涉及到多旅行商問題（Multiple Traveling Salesman Problem,MTSP）和移動目標旅行商問題（Moving Target Traveling Salesman Problem,MTTSP）。針對MTSP 問題，目前主要采用啟發式算法進行求解，如雜草入侵算法[9]、人工蜂群算 法[10]、遺傳算法[11]、粒子群算法[12]等。例如，Bourjolly等[13]建立了具有時間相關性的MTTSP 問題模型，并使用禁忌搜索算法求解航天器在軌服務問題；Moraes 等[14]通過對比遺傳算法、蟻群算法、模擬退火算法在解決MTTSP 時的表現，得出遺傳算法在求解質量和求解效率等方面表現最佳。

本文基于MTTSP 和MTSP 模型，構建了多無人機巡檢路徑規劃模型（Multiple UAVs’ Path Planning,MUPP）。針對多無人機巡檢任務分配的需求，提出兩階段算法（Two Phase Algorithm,TPA）以提升算法的求解性能，并以長江干線航道斷面為例進行仿真，驗證模型的科學性和算法的有效性。

1 模型構建

1.1 模型假設

為了將實際問題理論化，便于模型構建，需對模型做出以下假設：

1）無人機前往檢測時的速度恒定且快于船舶；在檢測時間內，船舶按當前速度和航向行駛；無人機到達目標點后能夠完成船舶燃油品質的檢測；

2）無人機跟隨檢測時的速度與被檢測船舶速度一致，檢測所需時間不隨船舶變化；

3）旋翼無人機不同于固定翼無人機，其轉彎半徑較小，故不考慮最小轉彎半徑約束；不考慮飛行過程中威脅、地形和氣候等因素的影響。

1.2 模型建立

基于上述假設，多無人機巡檢路徑規劃（Multiple UAVs’ Path Planning,MUPP）問題的表達如下：

式中：

L—巡航路徑總長度；

K—無人機的集合；

M—無人機的數量；

dij—船舶si與船舶sj的距離；

T—檢測總時間；

Q—無人機續航能力；

式（1）為目標函數，表示優化的目的為多無人機總飛行長度最小。約束條件（2）表示對巡檢總時間T 進行限制。約束條件（3）表示每艘船舶僅能被一架無人機檢測一次。約束條件（4）、（5）表示路徑中每個船舶只能分別與兩個船舶相連，保證路徑的連續性。

1.3 數據準備

MUPP 模型中的L、T值按照圖1 方法計算，圖中tr、dr為無人機返回基站的時間和距離；tj為無人機飛往船舶的時間；td為無人機跟隨檢測的時間；vd為無人機速度；vH(i)為船舶H(i)的速度；tH(i)為檢測時間；dH(i)為飛行距離。

圖1 無人機飛行總長度和時間的計算方法Fig.1 Calculation method of total flight length and time of UAVs

其中，飛行時間tf的計算原理為：已知無人機檢測完H(i)船后的位置D 點坐標為(x1,y1)，待檢船舶H(j)位置O 點坐標為(x2,y2)。以(x2,y2)為原點，H(j)船舶航行方向為y 軸建立直角坐標系，船舶H(j)勻速前行至檢測點M 后，無人機飛行方向與y 軸夾角為θ，得到方程式（6），化簡可得飛行時間tf的計算式（7）。計算簡圖如圖2 所示。

圖2 飛行時間計算示意圖Fig.2 Schematic diagram of flight time calculation

1.4 模型求解

MUPP 模型中的路徑規劃問題屬于NP-hard 問題，本文采用融合了遺傳算法（GA）和模擬退火算法（SA）兩種算法的混合模擬退火遺傳算法[15]（Simulated Annealing Genetic Algorithm,SAGA）作為基本算法。該改進遺傳算法同時具有遺傳算法的全局搜索特性和模擬退火算法的局部搜索特性。

此外，由于MUPP 模型需要多無人機協同完成巡檢任務，故涉及到任務分配的問題。基于改進的聚類算法，將該任務分配問題抽象為MTSP 后，再將MTSP 轉化為多個TSP 進行求解。

考慮到船舶運動特點及無人機任務分配的均衡，在上下行船舶數量較為接近的條件下，提出一種基于上下行船舶分組和任務均分的K-means++算法。其主要思想為將船舶數據集分為上下行兩組數據集，再引入容量限制要求，使得聚類形成的集合中船舶數量較為均勻。GT-K-means++算法的具體步驟如下：

設V=｛v1|i=1,2,…,n｝表示待檢測船舶節點集，初始簇集為s=(c1,c2,…,cj,…,ck)，k 為類簇數量。初始類簇cj的中心為μj,j=1,2,…,k，每個類簇中節點的個數qj=0。

1）根據船舶航向數據，將船舶節點集分為上行船舶節點集Vup和下行船舶節點集Vdown；接著初始化類簇容量，即每個類簇中船舶數量上限為Q=ceil(n/k)。

2）對上行船舶節點集Vup，初始聚類中心是該集合中隨機選取的一個船舶點；用最短距離D(x)表示每個船舶點與已有聚類中心之間的距離，然后用輪盤賭法選出下個聚類中心點；重復以上操作，直到選出上行船舶的kup個聚類中心點。

3）對下行船舶節點集Vdown，按照步驟2 中的方法選出下行船舶的kdown個聚類中心點。

4）針對各上行船舶點vi?Vup，按照公式（8）計算各船舶點vi到其類簇中心μj的歐幾里得距離ρ(vi,μj)。

接著將所有船舶點按照距離值從小到大排序。

5）聚類。取到類簇中心μj距離值最小的船舶點vi，并令qj=qj+1，判斷類簇cj當前船舶點個數是否滿足qj≤Q。若滿足，則將該船舶點vi分配給類簇cj；若不滿足，則將該船舶點vi到類簇cj的距離值設為無窮大，重新排序并分配下一船舶點，直到分配完所有船舶點。

6）根據公式更新聚類中心點jμ的坐標，類簇cj中船舶點數量是，若其中心點坐標沒有發生變化，則聚類完成收斂，輸出聚類結果；否則，轉4）。

7）針對各下行船舶點vi?Vdown，按照4）至 6）的方法完成聚類，并輸出聚類結果。

綜上，本文采用了GT-K-means++算法對水域內待檢測船舶點進行聚類，將船舶分成多個子集后，運用SAGA 對每一個船舶子集進行無人機巡檢路徑規劃，從而構成“GT-K-means++”+“SAGA”的二階段算法（Two Phase Algorithm,TPA），以此來分段求解MUPP。

2 實驗分析

2.1 研究范圍

本文選擇中國長三角船舶排放控制區范圍內長江干線航道南京段進行仿真實驗（如圖3 所示），并作參數影響分析。研究范圍中頂點經緯度為：32.25241 °N，119.156 °E；32.23642 °N，119.22114 °E；32.21797 °N，119.21249 °E；32.2372 °N，119.15147 °E。研究范圍的斷面長度約4 km，寬度約1.3 km。船舶航行數據均源自中國交通運輸部海事局的AIS 信息服務平臺。圖中圓形點表示設定的無人機基站位置，黑色線段為航道邊界線，船形標志表示船舶分布位置。

圖3 基于AIS 的研究范圍和船舶實時位置Fig.3 The research scope based on AIS and the real-time positions of the ships

2.2 仿真結果分析

算法采用Python 語言編寫，相關參數與模型參數設置如表1 所示，結合船舶AIS 航行數據，給范圍內所有船舶設定相應的編號、速度以及坐標。

表1 仿真參數設置Tab.1 Simulation parameter settings

TPA 算法迭代收斂圖如圖4 所示，從迭代結果可知，約350 次后收斂，表2 為使用TPA 算法運算得到的各無人機的巡檢路徑及所需時間。計算結果顯示，各無人機的任務量均為9，各無人機巡檢時間平均值19.5 min，方差為1.57，任務量較為接近，且滿足無人機的續航時間要求，四架無人機可均衡、高效地完成巡檢工作。

表2 巡航路徑和時間Tab.2 Cruise path and tim

圖4 TPA 迭代結果Fig.4 TPA iteration results

2.3 模型參數影響分析

1）船舶數量及無人機數量對TPA 結果的影響

假設各無人機速度vd均為20 m/s，單船檢測時間td為3 s。分別在不同場景下運行10 次程序，得到最優的路徑總時間T，如表3 所示。結果顯示當船舶數量大于60 艘時，出現2 艘無人機無法滿足巡檢要求的情況，并隨著船舶數量的增加，所需無人機數量不斷增加，當船舶數量達到100 艘時，需4 臺無人機才能完成巡檢任務。

表3 16 種場景下的最優路徑時間Tab.3 Optimal path time under 16 scenarios

可以發現，隨著船舶數量的增加，無人機巡檢總時間T 也呈增長趨勢。雖然場景15 中船舶數量比場景11 中多20 艘，其總時間T 卻比場景11 少3.7 min，表明船舶數量是影響巡檢路徑時間的主要因素。當航道中船舶數量為40 艘時，配備兩架續航時間為40 min 的無人機即可滿足巡檢要求；當船舶數量為60 艘時，應配備3 架續航時間為40 min的無人機或2 架續航時間為60 min 的無人機；而對于待檢測船舶數量較多的情況（80 艘、100 艘），則需要配備4 臺續航時間為40 min 的無人機才能夠順利完成巡檢任務。

當船舶數量為60 和船舶數量為100 時，無人機數量對巡檢時間如圖5 所示。圖中方形、圓形線段為兩種場景下T 值的變化曲線，三角形線段為每增加一架無人機時T 值下降率的變化曲線。

圖5 不同無人機數量下的總巡檢時間Fig.5 Total inspection time under different numbers of UAVs

從圖5 可以看出，無人機數量的增加會減少巡檢路徑總時間，當無人機數量由3 增加至4 時，T值下降率約51 %，變化最為顯著，而在無人機數量由4 增加至8 時，T 值下降率保持在14 %左右，巡檢效率提升不算明顯。由于內河航道中船舶流量有明顯的高峰期與低谷期，配備2 架或3 架無人機不能夠滿足較多船舶數量情況的巡檢需求，故每個無人機基站配備4 臺無人機最為合理。

2）船舶速度和無人機速度對TPA 計算結果的影響

以場景15 為例，在vd為14、18、22、26、30 m/s的五種情景下，分別求出td為3、6、9、12 s 時的路徑總時間T，實驗結果如圖6 所示。

圖6 不同飛行速度下的總檢測時間Fig.6 Total detection time at different flight speeds

圖6 顯示，在vd相同的情況下，td每增加3秒，T 值也隨之增長，從3 s 增長至6 s 時，T 值增長量較小，約6 分鐘；從6 s 增長至9 s 時，T 值的增長量約23 分鐘；從9 s 增長至12 s 時，T 值的增長量更多，約42 分鐘。可見td的增加表示無人機跟隨檢測所需長時間的增加，同時由于其他船舶行駛得更遠，也會導致巡檢時間的延長。當td保持不變時，隨著vd的增加，T 值逐漸減小，這是由于飛行速度提高導致檢測間隙時間縮短所致；此外，在vd由14 m/s 增加至22 m/s 的過程中T 值縮減幅度較大，但隨著vd值的繼續增加，T 值縮減幅度呈現逐漸減小趨勢。綜上分析，為提高巡檢效率，建議執法部門選用速度22 m/s 的無人機，并選配檢測設備響應時間小于6 s 的檢測儀來完成巡檢工作。

3 結語

為解決內河船舶尾氣排放監測工作中無人機自動化巡檢的部署方案問題，本文以無人機總飛行長度最小化為目標函數，提出了一種多無人機動態路徑規劃模型（MUPP），并采用“GT-K-means++” +“SAGA”的二階段算法（TPA）對這一NP-hard 問題進行求解。以長三角船舶排放控制區范圍內長江干線航道南京段為例，基于實時AIS 船舶航行數據，利用TPA 算法進行實例分析，結果顯示：待檢船舶數量的增多會對無人機的數量及續航時長有更高的要求，考慮到航道巡檢需求、無人機的巡檢效率，每個基站配備4 臺速度為22 m/s 的無人機，搭載檢測耗時少于6 s 的高性能檢測儀器較為合理。