多無人機協同目標分配與航跡規劃

許先靜，姜海燕

(福州大學 電氣工程與自動化學院，福州 350000)

0 引言

無人機具有體積小、重量輕、成本低等諸多優點，被廣泛運用于農業、軍事、飛行表演等多個行業。在很多實際場景中，單個無人機無法滿足作業需求，需要多無人機的協同作業，而協同目標分配是多無人機協同作業的基礎和保障。無人機的目標分配問題是一個約束性多且復雜的優化問題，其解空間隨無人機和目標總數的增加呈現指數級增加，求解最優解的時間也呈指數級增加。因此，無人機的目標分配問題是一個多參數、多約束的非線性多項式完全(Non-deterministic polynomial Complete，NPC)問題[1]。在進行目標分配時，需要考慮的因素很多，如無人機數量、任務類別、飛行環境等，同時還需要兼顧航行代價、分配算法以及各種協同約束條件等。

目前針對無人機的協同目標分配主要有集中式和分配式。集中式系統可以統一控制，缺點是擴展性差；而分布式系統的擴展性強，但對通信系統要求較高[2]。目前無人機的協同目標分配的解決方案主要有3種類型，即傳統的數學規劃法、市場機制法和智能優化算法。傳統的數學規劃法應用于集中式系統中，具體可分為匈牙利法[3]、分支界定法、隱枚舉法、混合整數線性規劃法、動態規劃法等；市場機制法應用于分布式系統，具體可分為合同網協議[4]、拍賣法等市場機制法[5]；智能優化算法既適用于集中式也適用于分布式，具體包括人工神經網絡算法、模擬退火算法[6]、遺傳算法[7-9]、差分進化算法、和聲搜索算法[10]、蟻群算法、A_算法[11-12]、粒子群優化算法[13-16]等，以及由若干算法組合而成的混合優化算法[1，17]。

無人機的編隊切換作為無人機的協同目標分配問題的等效(將N個目標分配給N個無人機)，大部分可以歸納為：給定無人機的起點和終點，賦予一定的目標任務、設定障礙或者航行代價等條件。目前許多的無人機的研究多基于二維平面，而無人機實際是在三維環境中執行各種協同任務的，所以在三維空間中研究無人機的協同目標分配，將會更具有理論研究意義和實際應用價值。對于多無人機隊形變換的評價標準比較復雜，因為縮小某一單架無人機的路徑長度，必然也會增加其它無人機的飛行路徑長度。在單個無人機的最大飛行路徑問題上，文獻[18]旋翼無人機協同任務指派問題研究與算法改進中，通過將匈牙利算法中代價矩陣各元素值替換為各自值的m次方，有效縮小了單個無人機的最大飛行距離，但m取值不同，其得到的結果也是不同的；文獻[19]多旋翼無人機編隊隊形變換算法研究中，運用粒子群算法實現隊形變換時間最小，縮小了單個無人機最大飛行距離，但該研究缺乏飛行距離整體性的考慮。

本文從協同目標分配出發,研究無人機編隊切換問題，發現標準的匈牙利算法在多無人機編隊切換應用中，存在部分無人機分配的飛行路徑過長的問題，同時由于無人飛行的能量消耗大于無人機懸停時能量消耗。從而導致個別無人機的電量下降迅速，因此相較于其他無人機會提前降落，而且會導致其他無人機在編隊切換過程中等待時間過長，影響整體編隊飛行的時長和編隊切換時長。因此縮小單個無人機的最大飛行路徑，既可以縮短編隊切換的時間，還可以平衡無人機的耗能，延長整體的編隊時長。利用匈牙利算法本身的特性，求解總的飛行距離最短，將目標函數改為求解在滿足約束條件下，單個無人機的最大飛行距離最小，對匈牙利算法進行改進。改進的匈牙利算法可以使得在總的飛行距離盡量小的前提下，盡可能縮減單個無人機的飛行距離，從而有效地降低單個無人機的最大能耗，減少無人機在編隊切換時懸停等待的時間和編隊切換的時間，有效延長整體編隊飛行的時間。適用于每架無人機到達各自目標空域點后需要懸停，等待所有無人機就位后再開始執行任務，比如多無人機定點拍照、監控或測量、協同打擊等時間協同性強的應用場景。而且算法簡單，沒有過多的參數設置，無需占用太多的計算資源，適合無人機嵌入式系統的實時計算。

1 無人機編隊切換問題描述和算法改進

1.1 隊形切換的數學模型及匈牙利算法

無人機的編隊切換問題，主要指的是在已知起點和目標位置的情況下進行位置移動，形成新的隊形，主要影響編隊切換效率的因素是每個無人機目標位置的選取。標準的匈牙利算法和其他算法大多數都只考慮無人機總飛行路徑最短，并在此基礎上去選取每個無人機的目標位置。然而考慮無人機總飛行路徑最短時，就會存在部分無人機分配的飛行路徑過長的問題，導致個別無人機電量迅速下降，提前降落且造成其他無人機等待時間過長，影響整體編隊的飛行時長。因此以單個無人機最大的飛行距離最小為優化目標，并滿足多約束條件，保證總移動距離盡可能小的情況下，可以延長整體編隊的飛行時長。

為了簡化問題的研究，假設每個無人機都是勻速飛行，忽略轉彎以及外界環境的影響并不考慮碰撞。記無人機編隊中各無人機的編號為i(i=1,2,3,...,n)，目標位置編號為j(j=1,2,3,...,n)；無人機的初始位置坐標(Xa,Ya,Za),無人機目標的坐標(Xs,Ys,Zs)?？梢杂嬎愠龅趇架無人機飛往第j個目標的路徑長度：

(1)

匈牙利算法是基于Hall定理中充分性證明的思想，它是部圖匹配最常見的算法，該算法的核心就是尋找增廣路徑，它是一種用增廣路徑求二分圖最大匹配的算法[20-21]。匈牙利算法在多無人機編隊切換應用中是將對于目標位置的分配抽象成指派問題。即：

1)有n個飛機飛往n個目標。

2)每個飛機只能飛往一個目標，每個目標只能有一個飛機飛往。

3)已知每個飛機飛往目標的能耗值。

4)求如何分配使得總飛行距離最短。

將上述問題轉換成數學模型可以得到：

(1)無人機/目標機對定義一個二值函數Xij，

(2)則總飛行距離最短的性能指標函數為：

(i=1,2,3,...,n;j=1,2,3,...,n)

(2)

(3)需要滿足的約束條件有：

(3)

(4)

Xij=0,1

(5)

滿足上述約束條件去求性能指標函數最優解，庫恩(W.W.Kuhn)于1955年提出了指派問題的一種解法，他引用了匈牙利數學家康尼格(D.konig)一個關于矩陣中0元素的定理：系數矩陣中獨立0元素的最多個數等于能覆蓋所有0元素的最少直線數。這解法稱為匈牙利算法。將其運用到多無人機編隊切換應用中，其求解的過程是：

1)先將每個飛機飛往每個目標的飛行距離求出并賦值為代價矩陣H(n×n)。行代表飛機編號，列代表目標編號。

(i=1,2,3，...，n;j=1,2,3，...，n)

(6)

2)把代價矩陣進行行歸約和列歸約，即每行減去該行的最小值，再把每列減去該列的最小值。

3)進行試指派，用盡可能少的直線去覆蓋矩陣中所有的0，如果所用的直線為n條則停止，得到最優匹配結果A，如果少于n，則繼續下一步。

4)在第三步中得到的代價函數中找尋最小的元素X，然后將直線未覆蓋的元素都減去該元素，并在直線相交處加上該元素。再更新代價矩陣后重復第三步。

圖1為算法流程圖：

圖1 標準匈牙利算法流程圖

圖2是一個4個飛機飛往4個目標的指派問題中匈牙利算法的過程。左1為代價矩陣，中間的圖是完成求解過程第二步和第三步后的矩陣，由于圖中只有3條直線，而n等于4，故執行第4步，減去最小元素2，并在直線交點處加上2，更新矩陣，得到右1圖。然后再從該矩陣中利用匈牙利算法找出最優匹配。

圖2 匈牙利求解過程

1.2 改進的匈牙利算法設計

標準的匈牙利算法只解決了在滿足約束條件的情況下性能指標函數Z的最優解，得到的分配結果A：

A=[S1J,S2J,S3J,...,S18J]，(j=1,2,3,...,n)

(7)

注：S1J代表編號為1的無人機飛往編號為j位置。

其中A集合中最大值即單個飛機飛行的最大飛行距離。

Smax=max(A)

(8)

考慮到時間的協同性，當無人機到達既定位置時需要等待其他無人機到達后統一開始執行任務，則完成任務的總時間T為：

T=Smax/v

(9)

而無人機無論是處于懸停還是飛行都需要損耗能量，則整個過程的能量損耗為：

(10)

注：k1代表無人機飛行能耗，k2代表無人機懸停能耗。

顯然無人機的總的耗能、所有無人機飛行的總距離、單個無人機最大的飛行距離存在關聯關系。而且由于無人機的飛行耗能要大于懸停耗能，因此減小單個無人機的最大飛行距離，不僅可以減小編隊切換的時間，還可以減小總的能耗及減小單個無人機的最大能耗，延長無人機整體編隊的時間。

針對標準的匈牙利算法，需要解決的問題就變成在滿足約束條件下，如何減小Smax的值，并保證性能指標Z的值盡量小(無人機飛行的總飛行距離盡量短)。此時優化的性能指標變為:

f=min[max(A)]

(11)

由上述分析可知，標準匈牙利算法在多無人機編隊切換應用中只考慮了總飛行距離最短，而沒有考慮在任務分配后單個飛機飛行的最大飛行距離。因此造成個別無人機的電量下降迅速相較于其他無人機提前降落且會造成其他無人機等待時間過長，影響整體編隊飛行的時長。所以在算法的改進過程中將單個飛機飛行的最大飛行距離考慮在內就可以解決這個問題。

改進算法的核心在于如何縮減單機最大飛行距離，本文通過標準的匈牙利算隨后將代價矩陣中比該值大的元素都賦值一個更大的元素M。在矩陣更新后，再通過標準的匈牙利算法進行求解。重復上述步驟，直至分配結果中出現元素M，接著返回上一次的分配結果，即可得出最優解。該算法可以保證在總的飛行距離盡可能小的情況下，縮小單個飛機飛行的最大飛行距離。算法具體求解的偽代碼如下：

begin

輸入：無人機的初始位置坐標(Xs,Ys,Zs),無人機目標的坐標(Xa,Ya,Zs)，

初始化：無人機數量n，矩陣H為n×n的零矩陣，最大元素M=0，矩陣A為1×n的零矩陣

For i=1 to n:

For j=1 to n:

Hij←Sij//將得到的距離值賦值到代價矩陣中

End for

End for

M←10×max(H)//M賦值為矩陣元素最大值的10倍

While Ture：

A←proceduce(Hungarian method)//進行標準匈牙利算法，并將結果賦值給A

If (max(A)=M)： //判斷單個無人機飛行最大距離是否等于M

Break //等于則返回上一次結果，結束循環

End if

For i->1 to n:

For j->1 to n:

If (Hij≥max(A)),

Hij←M//將大于單個無人機最大飛行距離的值都賦值為M，更新代價矩陣

End if

End for

End for

a=A

End while

輸出：最優分配結果a

End

注：M必須大于1倍max(H)，這樣才能保證第一次循環不會報錯，若小于等于1倍的maxH，第一次循環的分配結果的maxA有可能等于M導致循環報錯。

1.3 算法說明

設代價矩陣H為：

(i=1,2,3,...,n;j=1,2,3,...,n)

(12)

由匈牙利算法計算結果為：

A=[S1J,S2J,S3J,...,S18J]，(j=1,2,3,...,n)

(13)

由標準的匈牙利算法可知，最大元素Smax=max(A)出現在最優解的結果中時，是因為在步驟1進行行歸約、列歸約和步驟3更新代價矩陣時該元素Sij被賦值為0。

即在步驟2進行試指派的代價矩陣中：Sij←0將大于等于Sij的元素都賦值為M記此時代價矩陣為H′，再進行匈牙利算法計算。

若分配結果中max(A′)=M時，說明在步驟1進行行歸約、列歸約和步驟3更新代價矩陣時某個為M的元素被賦值為0。此時在滿足指派問題的約束條件下，必須將代價矩陣H′中等于M的元素添加到最優解中。

又因為已知元素等于M的值在原代價矩陣H中的值都是大于Sij的；所以上一次的分配結果A即為滿足指派問題約束條件下的最優解。

2 實驗結果與分析

2.1 實驗背景

本次仿真模擬18架無人機的編隊切換，每架無人機到達各自目標空域點后需要懸停，等待所有無人機就位后再開始執行任務，比如多無人機定點拍照、監控或測量、協同打擊等。本文提出的方法是在已知無人機自身的位置和周圍環境信息下的目標分配與航跡規劃方法。無人機的初始位置坐標由自身攜帶的定位系統確定，目標位置坐標由無人機控制系統給定。

假定無人機初始位置坐標：

A_x=[100,140,180,220,260,300,100,140,180,140,260,260,260,140,140,140,180,220]

A_y=[400,400,400,400,400,400,100,100,100,350,300,250,200,150,200,300,250,250]

A_z=[500,500,500,500,500,500,200,200,200,450,400,350,300,250,300,400,350,350]

將上述對應初始位置的無人機進行編號為1-18。

無人機目標點位置坐標：

B_x=[100,140,180,220,260,300,100,140,180,100,260,260,260,220,120,160,200,240]

B_y=[400,400,400,400,400,400,100,100,100,350,350,150,100,100,150,190,260,315]

B_z=[500,500,500,500,500,500,200,200,200,450,450,250,200,200,250,290,360,415]

將上述對應目標點位置編號為任務1～18。

編隊飛行實驗從隊形F到隊形Z，如圖3～4所示。

圖3 編隊F 圖4 編隊Z

為了驗證算法的性能，針對上述問題設計了粒子群算法并且分別用標準的匈牙利算法、m次方的匈牙利算法[18]、和本文中改進的匈牙利算法進行目標位置的分配，無人機編隊切換航線由各算法的目標分配結果確定，以二維空間的下匈牙利算法為例，其坐標值為上述坐標中去除z軸坐標，將初始位置坐標和目標點坐標賦予匈牙利算法，得到目標分配結果。算法輸出結果如下：任務1～18依次分配的無人機編號為：

[1,2,3,4,5,6,7,8,9,10,11,13,12,15,16,18,17,14]。根據算法分配結果的飛行路線如圖5所示，圖中圓表示無人機的初始位置，*號表示無人機的目標位置。對應的飛行距離為：

[0,0,0,0,0,0,0,0,0,40,50,150,50,20,22,101,22,150]

圖5 無人機編隊切換示例

2.2 粒子群算法

粒子群算法的思想是模擬鳥群捕食的行為演化而來，捕食的鳥群都是通過各自的搜索與群體的合作最終發現食物所在的位置。通過一種無質量的粒子來模擬鳥群中的鳥，每個粒子含有兩種屬性，速度和位置，速度表示移動的快慢，位置表示移動的方向。算法中將無人機當作粒子，首先，每個無人機在空間中搜索當前個體極值，每個個體極值通過比較，以此找到當前最優的個體極值作為當前全局最優解，根據當前全局最優，粒子來調整自己的速度和位置，往最優化方向靠近來求得最優解。

而本次任務要求分配方案要滿足無人機飛行路徑相近，因此，在粒子搜索最優解(即路徑和最小)的過程中加入第一優先級：滿足相近條件。在相近的條件下，再來找到路徑和最小。PSO算法流程圖，如圖6所示。

圖6 PSO算法流程圖

2.3 實驗結果

將無人機的初始位置和目標位置分別賦予給匈牙利算法、m次方的匈牙利算法、粒子群算法和本文中改進的匈牙利算法，通過各算法給出的目標分配結果指定無人機的航線并計算無人機的飛行距離。文獻[18]中比較了當m分別等于2、3、4情況下算法的匹配結果，當m越大時，其結果變化趨勢與之前相同，綜合考慮優化參數和實際中對路線的要求，最終將m的值設定為2。本次實驗中引入該算法，并且分別比較當m等于2和3兩種情況下的匹配結果。如表1所示為上述各算法的目標分配結果，可見各算法的匹配結果不見相同。

由圖7可知，無人機按照各算法的匹配結果飛至指點目標位置的飛行距離中，匈牙利算法的單個無人機的飛行距離最長為212.13 m，其次是粒子群算法為156.84 m，m次方的匈牙利算法和本文提出的改進的匈牙利算法均為141.42 m。單個無人機的飛行距離影響著無人機的編隊切換時長和懸停等待能耗，同時也可以降低各無人的能耗差，延長無人機的整體編隊飛行時長，m次方的匈牙利算法和本文提出的改進的匈牙利算法在這方面均表現良好。

圖7 根據各算法匹配結果無人機飛至指定目標的距離

無人機的編隊切換時間取決于最大飛行距離，并且無

表1 各算法的任務匹配結果

人機的能耗取決于懸停等待時間和飛行的距離。由表2可知，粒子群算法相比于匈牙利算法，其最大飛行距離雖然縮短了，但粒子群算法不僅有著大量的參數設置，并且其結果不如其他的改進算法，計算結果和效率明顯欠缺。相對于其他算法，改進的匈牙利算法得到的最大飛行距離最小，且總距離只比匈牙利算法多了11 m；m次方匈牙利算法的匈牙利算法的最大距離雖然和改進的匈牙利算法一致，但其總的飛行距離比較大，影響無人機編隊切換的總體能耗。改進的匈牙利算法不僅縮減了無人機的等待和編隊切換時間，還降低了無人機懸停等待能耗和各無人機之間的能耗差，同時其總的飛行距離和平均距離較小，無人機的飛行能耗也較小，降低了整體編隊的耗能，延長了整體編隊飛行時間。

表2 算法結果對比

經過實驗對比，本文提出的算法在面對每架無人機到達各自目標空域點后需要懸停，等待所有無人機就位后再開始執行任務這類場景下的目標分配表現良好。不僅可以縮短無人機的編隊切換時間，還可以降低無人機的總體能耗和各無人機的能耗差，延長無人機的整體編隊飛行時間。并且該算法簡單，無需占用太多的計算資源，相比于計算智能領域的各種算法如粒子群算法、蟻群算法等，沒有過多的參數設置并且其分配結果不會因為參數的改變而變化，適用于無人機嵌入式系統實時計算。能夠適應時間協同性強的目標分配系統；如無人機編隊表演時能延長整體編隊飛行時長，縮短編隊切換時間；對敵方目標進行協同打擊時能縮短無人機的集結時間，快速完成打擊任務等。

3 結束語

標準匈牙利算法在無人機編隊切換中存在問題，即單個無人機飛行距離過長，造成其他無人機懸停等待過長，影響編隊切換的時長又以及因為各個無人機能耗差別導致個別無人機的電量下降迅速而提前降落，影響整體編隊飛行時間。單個無人機的飛行距離可通過粒子群等智能優化算法有效縮短飛行距離，但存在參數的選取等問題，在一部分系統中不能很好的得到最優解，且只考慮了單個無人機的飛行距離，未將總的飛行距離考慮在內。通過仿真實例，利用匈牙利算法特性，將目標函數改為求解滿足約束條件下單個無人機的最大飛行距離最小。經過改進的匈牙利算法對無人機編隊切換有著很好的表現，雖然標準的匈牙利算法可以得到飛行總距離最小的最優解，但使用改進的匈牙利可以在總距離變化不大的情況下，最大化的減小單個無人機的飛行距離?？梢杂行У慕档蛦蝹€無人機的最大能耗，減少無人機在編隊切換時懸停等待和編隊切換的時間，有效延長整個編隊飛行的時間。本文提出的方法相較于匈牙利算法、粒子群算法及m次方的匈牙利算法更適用于無人機飛行時的編隊切換，能夠縮短無人編隊切換的時間，提升整體飛行的時長；同時適用于具有時間協同性的目標分配系統及性能指標為最小化單個個體指標的協同目標分配系統，算法簡單具有良好的適用性。