楊少龍, 黃 金, 向先波, 李偉超

(華中科技大學船舶與海洋工程學院, 湖北 武漢 430074)

0 引 言

隨著一帶一路建設,海上經(jīng)貿(mào)活動不斷增加,隨之海上事故頻發(fā)。在發(fā)生海難事故時為了最大限度減少生命財產(chǎn)損失,海上應急救援需要準確、快速地估計搜尋區(qū)域,設計合適區(qū)域搜尋規(guī)劃路徑提高海上搜尋成功率(probability of success,POS)并減少搜尋時間,使搜尋力量(如無人艇)以最短時間完成對高概率區(qū)域覆蓋式搜尋。

海上遇險目標搜尋區(qū)域確定是成功搜尋的重要前提。Breivik等人采用普林斯頓海洋模型分析海洋環(huán)境,建立了海上搜救模型,并利用蒙特卡羅模擬確定搜尋區(qū)域。肖方兵等人建立包含環(huán)境擾動的遇險目標漂移模型,通過蒙特卡羅隨機粒子仿真法結合最優(yōu)搜尋矩形規(guī)則確定搜尋區(qū)域。但蒙特卡羅隨機粒子仿真法產(chǎn)生搜尋區(qū)域中少量邊緣粒子遠離中心高密度粒子區(qū)域,數(shù)量稀少且分布發(fā)散。隨著時間推移,邊緣粒子分布會更加分散,顯著擴大待搜尋矩形區(qū)域面積。這對搜尋力量投入帶來巨大負擔,對于分秒必爭的救援任務而言更需要對粒子高密度聚集區(qū)域以最短時間完成搜尋覆蓋任務。

搜救力量是海上遇險目標搜尋規(guī)劃執(zhí)行主體,但有人搜尋力量難以勝任大面積、長時間、惡劣海況的區(qū)域搜尋任務。隨著無人系統(tǒng)快速發(fā)展,無人艇、無人機等新生力量為海上遇險目標搜尋帶來了新的機遇。以無人機為例, Bouzid等人、Martins等人研究無人機海上搜尋,但由于惡劣海況、惡劣天氣、續(xù)航等因素難以長時間工作。相比之下,無人艇可在水面航行,對惡劣海況適應能力更強且續(xù)航滿足大范圍的區(qū)域搜尋。李亞南等人研究了無人艇集群對包含障礙的區(qū)域進行動態(tài)規(guī)劃,較好地適應了動態(tài)環(huán)境,但具體到海上搜救應用場景,搜尋代價、區(qū)域覆蓋程度等是更加重要的關注點。

區(qū)域覆蓋搜尋問題在各領域廣泛存在,受到國內(nèi)外學者普遍關注。Choset等人首次利用boustrophedon細胞分解實現(xiàn)了機器人對區(qū)域的完全覆蓋,但梳狀規(guī)劃路徑未考慮機器人自身的操縱性限制。李亞南等人采用模仿獨居動物社會行為的分布式反集群研究了無人艇集群區(qū)域全覆蓋方法,同樣忽略了欠驅(qū)動無人艇操縱性導致的轉向路徑跟蹤問題。針對區(qū)域覆蓋路徑尋優(yōu),Li等人將航行距離與最大航行距離多目標優(yōu)化問題轉化為單目標優(yōu)化。彭輝等人將多目標航路代價估計值和有效時間估計值加權求和轉化為單目標優(yōu)化。但這類多目標加權轉化為單目標做法,在實際搜尋規(guī)劃任務中卻存在權值系數(shù)難以確定的問題。

為此,針對海上搜救區(qū)域搜尋規(guī)劃面臨的目標包含概率(probability of containment, POC)、搜尋探測概率(probability of detection, POD)表示搜救力量探測到搜尋目標的概率,以及搜尋代價等優(yōu)化需求,開展基于無人艇的多目標優(yōu)化搜尋規(guī)劃方法研究,在最少搜尋路徑代價下實現(xiàn)最大化搜尋POS=POC×POD,表示搜尋行動成功找到目標的概率。一方面,搜尋區(qū)域確定聚焦高密度區(qū)域及最佳POC;另一方面,搜尋規(guī)劃優(yōu)化實現(xiàn)總路徑最短且POD最大。首先,建立海上遇險目標漂移粒子分布位置,然后依次對搜尋區(qū)域和搜尋路徑進行優(yōu)化求解。對于搜尋區(qū)域,設計一種基于腐蝕膨脹的置信橢圓搜尋區(qū)域確定算法,實現(xiàn)區(qū)域內(nèi)粒子POC和單位面積粒子數(shù)達到全局最優(yōu)。對于搜尋路徑,為實現(xiàn)優(yōu)化搜尋線間距、遍歷次序等從而獲取全局最優(yōu)總路徑和POD采用非支配排序的遺傳算法(nondominated sorting genetic algorithm II, NSGA-II),該算法使用帶精英策略快速非支配排序從而實現(xiàn)快速、準確的搜索性能。本文最后通過對比常規(guī)矩形區(qū)域的梳狀路徑搜尋規(guī)劃,揭示所提方法在搜尋工作路徑等指標的顯著優(yōu)勢。

1 無人艇區(qū)域搜尋規(guī)劃總體思路

無人艇區(qū)域搜尋規(guī)劃總體思路如圖1所示。搜尋任務分為兩個階段:一是根據(jù)已知海況信息以及船只失事前最后一次上傳信息利用蒙特卡羅隨機粒子仿真獲得粒子初始位置分布,然后根據(jù)粒子分布特點采用聚類算法將粒子分類,最后采用規(guī)則區(qū)域確定算法鎖定高密度區(qū)域,使得搜尋力量盡可能投入到高概率區(qū)域;二是根據(jù)搜尋區(qū)域位置和形狀信息結合無人艇轉向數(shù)學模型,進行多目標搜尋路徑規(guī)劃,同時滿足最大化POD和最短總路徑。

圖1 水上遇險目標搜尋整體框架Fig.1 Framework of distress targets overboard searching

無人艇執(zhí)行搜尋路徑跟蹤任務,通過艇載視覺、雷達、紅外探測等傳感器探測落水人員。單次任務結束后根據(jù)搜尋結果重復上述步驟,直至滿足要求。

2 基于腐蝕膨脹的置信橢圓搜尋區(qū)域確定

2.1 海上遇險目標概率分布

海上遇險目標位置分布中包含多種不確定因素,導致目標漂移軌跡呈現(xiàn)多樣性。為此,在文獻[3]基礎上采用蒙特卡羅隨機粒子仿真法建立海上遇險人員漂移位置分布。預測時長0.5 h,風場、流場數(shù)據(jù)以及落水人員風壓屬性見文獻。如圖2所示,海上遇險目標概率分布呈現(xiàn)較為明顯的類簇特征,且邊緣粒子分布相對分散。

圖2 海上遇險目標概率分布Fig.2 Probability distribution of targets in distress at sea

2.2 海上遇險目標概率分布聚類

如圖2所示,若直接利用矩形區(qū)域框選包括所有粒子來確定搜尋區(qū)域,則會顯著增大搜尋代價,且在粒子低密度區(qū)域浪費搜尋力量。因此,采用聚類算法將所有粒子分簇,再結合分簇結果設計合理的搜尋區(qū)域。傳統(tǒng)K均值聚類算法是不考慮概率特性的計算點與類簇中心的歐式距離來判定歸屬,不適應本場景應用,因此采用高斯混合模型(Gaussian mixture model,GMM)進行聚類,假設GMM模型由組高斯模型組成,對應GMM概率密度函數(shù)如下:

(1)

其中,第組分模型概率密度公式為

(2)

式中:為權值;為均值;為方差;為數(shù)據(jù)維度(=2);為高斯模型數(shù);為描述二維坐標的樣本向量。

首先確定模型組數(shù)值,結合圖2分布特征設定=2,然后使用最大期望算法迭代求解直至產(chǎn)生的新聚類中心與上一次聚類中心差值小于預設常數(shù)0005,結束后將獲得所有粒子分簇的劃分結果。如圖3所示,粒子分布呈現(xiàn)兩個類簇。此外,少量粒子顯著遠離類簇中心且呈相對分散的特征,從搜尋角度而言需要投入相當大的搜尋代價才能完全覆蓋這些概率密度較低的分散區(qū)域,不利于提高搜尋效率。

圖3 基于GMM聚類粒子分布圖Fig.3 Distribution of GMM-based clustering particles

為此,進一步設計基于腐蝕膨脹的橢圓搜尋區(qū)域確定算法,聚焦高密度區(qū)域提高搜尋效率。

2.3 基于腐蝕膨脹的橢圓搜尋區(qū)域確定算法

基于腐蝕膨脹的置信橢圓搜尋區(qū)域確定算法分為3個步驟。

根據(jù)聚類結果分別求出各類數(shù)據(jù)均值以及數(shù)據(jù)協(xié)方差矩陣對應的特征值和特征向量。

獲取最大特征值并求出對應的特征向量,求出最大特征向量與軸夾角。

根據(jù)橢圓公式構造卡方檢驗,其中自由度為2,選取置信度為95%,可得橢圓長短半軸公式為

(3)

(4)

式中:,分別為橢圓長半軸和短半軸;為置信度;,分別為協(xié)方差矩陣的最大和最小特征值;,為腐蝕膨脹系數(shù)。

通過上述3步獲得置信橢圓中心、長軸和短軸,從而確定橢圓面積和邊界。通過調(diào)節(jié)腐蝕膨脹系數(shù)可以改變區(qū)域大小從而改變單位面積粒子數(shù)和POC。為了獲取全局最優(yōu)單位面積粒子數(shù)和POC建立了如下多目標優(yōu)化模型:

max POC

(5)

(6)

式(5)與式(6)分別表示優(yōu)化目標POC(橢圓內(nèi)粒子數(shù)除以隨機粒子總數(shù))和單位面積的粒子數(shù),其中表示產(chǎn)生隨機粒子總數(shù)。

采用基于非支配排序并結合擁擠度算子求解該多目標優(yōu)化問題,確定最佳腐蝕膨脹系數(shù)。根據(jù)橢圓區(qū)域特點,腐蝕膨脹系數(shù)范圍為08～12,步長001。以圖3中紅色粒子區(qū)域1為例,對不同腐蝕膨脹系數(shù)下POC和單位面積粒子數(shù)進行非支配排序并結合擁擠度算子得到優(yōu)化結果如圖4所示。

圖4 不同腐蝕系數(shù)下POC和單位面積粒子數(shù)Fig.4 POC and number of particles per unit area under different corrosion coefficients

基于最優(yōu)前沿面采用擁擠度算子對結果排序得到擁擠度最小值所對應的腐蝕膨脹系數(shù)即全局最優(yōu)腐蝕膨脹系數(shù),即圖4中點取值。選擇圖4中最優(yōu)前沿面上、、點繪制橢圓區(qū)域,如圖5對比可見,紅色橢圓區(qū)域相比于藍色區(qū)域單位面積粒子密度更高,可縮小待搜尋區(qū)域;與青色區(qū)域相比,紅色橢圓區(qū)域POC概率更高。

圖5 不同腐蝕膨脹系數(shù)下橢圓區(qū)域Fig.5 Elliptical region under different corrosion expansion coefficients

基于此方法優(yōu)化得出兩組聚類粒子的最佳置信橢圓區(qū)域如圖6所示。區(qū)域1長軸腐蝕膨脹系數(shù)=093,短軸腐蝕膨脹系數(shù)=099,POC概率為0469。區(qū)域2長軸腐蝕膨脹系數(shù)=092,短軸腐蝕膨脹系數(shù)=101,對應POC為0473。兩個區(qū)域整體POC概率為0942。

圖6 基于腐蝕膨脹的置信橢圓區(qū)域Fig.6 Confidence ellipsoid region based on corrosion expansion

3 無人艇區(qū)域覆蓋搜尋路徑規(guī)劃

3.1 問題描述

通過第2節(jié)確定最佳搜尋區(qū)域后,欠驅(qū)動無人艇需要設計合適的搜尋規(guī)劃路徑完成區(qū)域覆蓋任務。采用NSGA-II算法同時考慮搜尋總路徑與POD雙目標優(yōu)化問題,結合非支配排序和擁擠度算子求得非支配解更接近真實Pareto最優(yōu)解集,適用于本場景問題。

3.2 問題建模

欠驅(qū)動無人艇搜尋區(qū)域路徑規(guī)劃主要關注搜尋總路徑和POD兩個目標。其中,搜尋總路徑包含工作路徑和非工作路徑,工作路徑為無人艇在搜尋區(qū)域內(nèi)所有直線作業(yè)行路徑總和,非工作路徑為各作業(yè)行轉向路徑總和。執(zhí)行任務過程中需要滿足約束條件包括:① 轉向路徑最短;② 轉向次數(shù)最少;③ 每條作業(yè)行均被遍歷一次且不回起始行。約束條件中轉向路徑需依據(jù)無人艇操縱特性建立轉向數(shù)學模型,根據(jù)欠驅(qū)動無人艇轉向特性和作業(yè)行間距建立了U型和Ω型兩種轉向類型。如圖7所示,其中P1B1、P2B2為工作路徑,B1B2之間圓弧段路徑為非工作路徑,為當前作業(yè)行與下一條作業(yè)行之間的間距,兩個圓心距表示無人艇最小回轉直徑,=arccos()。

圖7 兩種類型轉向模式Fig.7 Two types of steering modes

受橢圓搜尋區(qū)域邊界特征影響,搜尋區(qū)域內(nèi)各條直線作業(yè)行長度不一(如和點),這會對無人艇的轉向路徑?jīng)Q策帶來影響。為此,建立適應不同高度差的無人艇轉向數(shù)學模型為

=|·tan-Δ|+

(7)

(8)

式中:、表示作業(yè)行路徑編號≠;表示從當前作業(yè)行駛出到駛入下一條作業(yè)行的轉向路徑總長度;=(π2+);Δ=|-|表示第作業(yè)行和第作業(yè)行高度差,其中式(7)和式(8)分別對應于圖7中Ω型轉彎模式和U型轉彎模式。

將作業(yè)行總路徑和POD作為優(yōu)化目標,以起始作業(yè)行、作業(yè)行遍歷順序和搜尋線間距為變量,以橢圓長軸方向為搜尋方向,可建立描述該問題數(shù)學模型。其中作業(yè)總行數(shù)計算式為

(9)

式中:為相鄰兩條搜尋作業(yè)行間距,即搜尋線間距;start為搜尋作業(yè)起始行,總行數(shù)計算為向下取整;為決策變量,表示是否包含從行到行的轉向路徑,其中=0表示不包含,=1表示包含;為無人艇掃海寬度;表示集合{1,2,…,};||表示集合中元素個數(shù),則整體優(yōu)化模型表述如下:

(10)

s.t.

=, ?,∈

(11)

(12)

∈{1,0},,∈;≠

(13)

式(10)為優(yōu)化目標函數(shù),要求總路徑最短且POD最大;式(11)表示轉向路徑具有對稱性即從行到行與從行到行轉向路徑長度相等;式(12)與式(13)表示每條作業(yè)行有且僅被遍歷一次且最終不回到起始行。

3.3 NSGA-II優(yōu)化算法設計

利用NSGA-II算法獲得最優(yōu)起始行、作業(yè)行遍歷順序和搜尋線間距,實現(xiàn)總距離和POD值最優(yōu),整個過程分為以下3步。

設定群體規(guī)模、編碼長度、交叉概率、變異概率以及最大遺傳代數(shù)等參數(shù),根據(jù)無人艇航行性能確定搜尋線間距范圍和間距步長。

利用遺傳算法產(chǎn)生初始起始行和作業(yè)行遍歷順序初始值,然后按照確定搜尋線間距計算在起始行確定條件下,遺傳算法確定最佳遍歷順序,由此確定每個子代總路徑和POD,采用非支配排序和擁擠度算子獲取父代和子代合并后新子代作為下一次計算的父代,反復迭代直至達到終止條件并保存最優(yōu)總路徑和POD。

依據(jù)搜尋線間距步長改變搜尋線間距重復步驟2直至達到搜尋線間距范圍上限。然后利用非支配排序并結合擁擠度算子獲得不同搜尋線間距下?lián)頂D度最小值即全局最優(yōu)總路徑和POD,最后輸出最優(yōu)值所對應的起始行、作業(yè)行遍歷次序和搜尋線間距,上述步驟偽代碼如算法1所示。

算法 1 無人艇路徑優(yōu)化算法01 輸入:N,Gen,Gen1,Pc,Pm,S;∥輸入種群大小N,起始行最大遺傳代數(shù)Gen,遍歷順序最大遺傳代數(shù)Gen1,交叉概率Pc,變異概率Pm,搜尋線間距范圍S02 輸出:s[];∥輸出Pareto非劣解集03 for i=S_min:S_step:S_max04 Chrom=crtbp(N,PRECI);∥得到初始種群05 for j=1:Gen06 MP=Mating(Chrom,FrontValue,CrowdD istance);∥二元錦標賽選擇07 SelCh=recombin(‘xovsp’,Pc);∥單點交叉08 Chrom=mut(SelCh,Pm);∥變異09 Sobj=object_calculation(Chrom);∥獲取每個種群初始行10 for m=1:Gen111 Tobj=TSP(Sobj, PRECI);∥獲取每個種群最優(yōu)遍歷順序所對應總路徑和POD12 end for13 FrontV=NonDominateSort(TObj);∥快速非支配排序

14 CrowdD=CrowdDistances(TObj,Fron tV);∥擁擠度距離計算15 end for16 [distance POD]=find_sort(FrontV,CrowdD);∥該搜尋線間距下最優(yōu)值17 save=[save; distance POD];∥保存最優(yōu)值18 end for19 best_FrontV=NonDominateSort(save_best);20 best_CrowdD=CrowdDistances(save,best_FrontV);21 [d P]=find_sort(best_FrontV, best_CrowdD);∥獲得不同搜尋線間距下最優(yōu)值

4 算法應用實例

算法求解在Matlab 2017b中實現(xiàn),仿真計算機為Intel(R) Core(TM)i5-6200U CPU @ 2.30 GHz、內(nèi)存12 GB、操作系統(tǒng)Windows10 64位系統(tǒng)。將世代距離(generational distance,GD)引入作為整體收斂性評價指標,GD值越小表明越接近真實Pareto最優(yōu)解集。此外,海上搜尋規(guī)劃通常以矩形區(qū)域框選感興趣范圍,且基于最小跨度理論規(guī)劃平行線掃描方向減少無人艇往復轉彎次數(shù)。因此,本文以橢圓最小外接矩形作為對比,長軸為長邊,短軸為短邊,且短軸是最小跨度方向,根據(jù)該外接矩形搜尋區(qū)域進行搜尋規(guī)劃性能比較。

NSGA-II算法中初始參數(shù)設定為:種群大小為100,最大遺傳代數(shù)為Gen為200,交叉概率Pc為0.9,變異概率Pm為0.05,搜尋線間距范圍為0.1～0.2 nmile,掃海寬度為0.2 nmile,區(qū)域1的橢圓長半軸為5.21 nmile,短半軸為2.36 nmile,區(qū)域2的橢圓長半軸為6.31 nmile,短半軸為4.56 nmile,無人艇最小回轉直徑為8 m。

以圖6所示橢圓區(qū)域1為例,區(qū)域覆蓋路徑規(guī)劃結果表明當搜尋線間距步長0.02 nmile時,在搜尋線間距為0.12 nmile時,總路徑和POD值達到最優(yōu),分別為337.04 nmile和0.811,作業(yè)行遍歷順序為按序依次遍歷,Pareto前沿面如圖8所示。

圖8 區(qū)域1優(yōu)化結果的Pareto前沿面及規(guī)劃路徑Fig.8 Pareto front of area 1 optimization result and planned path

NSGA-II算法優(yōu)化結果GD值為0.015,接近于0,表明所求得解接近真實Pareto最優(yōu)解集?；赑areto最優(yōu)前沿面解采用擁擠度算子對結果排序從而獲得擁擠度最小值即全局最優(yōu)解,根據(jù)所得最優(yōu)解對應的作業(yè)起始行、作業(yè)遍歷順序和搜尋線間距可得到搜尋規(guī)劃路徑。如圖8左上角區(qū)域1所示規(guī)劃結果,圖8右下角為路徑局部放大圖,可見無人艇轉向過程以U型連接兩條作業(yè)行,符合無人艇轉向特性,且具有較短的非工作路徑。

為了對比常規(guī)矩形區(qū)域和橢圓區(qū)域在搜尋工作路徑和POD等方面的指標差異,以橢圓區(qū)域1的外接矩形區(qū)域進行搜尋規(guī)劃,通過優(yōu)化作業(yè)起始行、作業(yè)行遍歷順序和搜尋線間距獲得全局最優(yōu)總路徑和POD,對比結果如表1所示。橢圓區(qū)域POD比矩形區(qū)域略大,但由于POC小于矩形區(qū)域,導致兩者POS相等。本場景下搜尋線間距遠大于無人艇最小回轉直徑,因此,兩種搜尋區(qū)域均采用U型按序依次遍歷完成路徑規(guī)劃。然而,對比總路徑,橢圓區(qū)域比矩形區(qū)域少23.4%,可以有效減少搜尋路徑代價。

表1 橢圓區(qū)域1和矩形區(qū)域1的搜尋指標對比

按照同樣方式對橢圓區(qū)域2進行計算,結果表明當搜尋線間距為0.12 nmile時,總路徑和POD值最優(yōu)且最優(yōu)值分別為719.98 nmile和0.811,所規(guī)劃路徑順序為按序依次遍歷。優(yōu)化結果所對應的Pareto前沿面如圖9所示,GD值為0.032。區(qū)域1和區(qū)域2總體POS為0.764,對應總路徑為1 050.42 nmile。

圖9 區(qū)域2優(yōu)化結果的Pareto前沿面及規(guī)劃路徑Fig.9 Pareto front of region 2 optimization result and planned path

對比兩個橢圓區(qū)域的外接矩形區(qū)域所規(guī)劃結果可知,矩形區(qū)域總體POS為0.769,總路徑為1 384.13 nmile。由此,相比常規(guī)矩形區(qū)域搜尋規(guī)劃方法,本文所提橢圓區(qū)域搜尋路徑規(guī)劃方法在POS幾乎不變的情況下,縮短搜尋路徑代價達0.318,為實際海上應急搜尋提供了高效的搜尋決策依據(jù)。

5 結 論

本文提出面向遇險目標單向搜尋的無人艇區(qū)域搜尋規(guī)劃方法。首先利用GMM聚類算法對蒙特卡羅隨機粒子仿真法產(chǎn)生初始粒子聚類,然后采用基于腐蝕膨脹的置信橢圓法確定被搜尋區(qū)域,最后采用NSGA-II多目標優(yōu)化算法對搜尋區(qū)域進行合理規(guī)劃。

主要貢獻和創(chuàng)新點有:① 采用GMM聚類算法對粒子分類,將粒子分類成多個區(qū)域,減少后續(xù)搜尋區(qū)域確定中不必要的搜尋面積;② 通過基于腐蝕膨脹的置信橢圓法確定全局下最優(yōu)搜尋區(qū)域,使得區(qū)域內(nèi)POD和單位面積粒子數(shù)全局最優(yōu),保證搜尋聚焦高密度區(qū)域;③ 建立面向橢圓區(qū)域邊界特征的欠驅(qū)動無人艇轉向模型,減少非工作路徑長度,提高搜尋效率;④ 采用NSGA-II多目標優(yōu)化算法對總路徑和POD進行優(yōu)化獲得全局最優(yōu)起始作業(yè)行、搜尋線間距和搜尋作業(yè)行遍歷次序,達到滿足一定搜尋成功率前提下,縮短搜尋路徑代價達31.8%,提高搜尋時效性。