杜楠楠，陳建，馬奔，王術(shù)波, 3，張自超

1. 中國農(nóng)業(yè)大學 工學院，北京 100083

2. 北京航空航天大學 自動化科學與電氣工程學院，北京 100191

3. 武漢大學 測繪遙感信息工程國家重點實驗室，武漢 430079

近年來，隨著無人機及其相關(guān)技術(shù)的發(fā)展，無人機已經(jīng)廣泛應用于國防軍事[1]、災后搜救[2]以及遙感[3-4]等領(lǐng)域。隨著太陽能電池效率的提高，環(huán)保、可以實現(xiàn)長距離航行、并且能夠執(zhí)行多種任務的太陽能無人機進入了人們的視野[5]。1974年，第1架太陽能無人機“Sunrise I”試飛成功[6]。2003年，瑞士蘇黎世聯(lián)邦理工大學開始研制小型無人機“Sky-Sailor”，并于2008年6月在瑞士完成了27 h的低空飛行[7]。2005年，美國AC公司成功測試了主要應用于民用遙感測繪的電動滑翔太陽能無人機“SoLong”[8]。2017年，同樣是由瑞士蘇黎世聯(lián)邦理工大學研發(fā)的太陽能無人機“Atlantik Solar”以連續(xù)飛行81.5 h打破了輕小型太陽能無人機的航時記錄[9-10]。

無人機作業(yè)中最重要的問題之一是路徑規(guī)劃。傳統(tǒng)的路徑規(guī)劃是點對點的，即規(guī)劃出無人機和終點之間的一條路徑[11]。在某些特定情況下，無人機需要獲得某一區(qū)域的全面信息，或在該地區(qū)進行全覆蓋作業(yè)，因此需要進行覆蓋路徑規(guī)劃?！芭８ā笔且环N典型的覆蓋路徑規(guī)劃方法，即通過往復的直線運動來將區(qū)域覆蓋完全[12]。在“牛耕法”中，一般選擇使轉(zhuǎn)彎次數(shù)最小的方向作為覆蓋飛行方向[13-14]。對于不含障礙物的凸多邊形區(qū)域，Gustavo等提出了基于車輛路徑規(guī)劃問題的、以作業(yè)時間最短為求解目標的多無人機覆蓋路徑規(guī)劃模型，并在該問題中考慮了無人機的起飛操作時間對總作業(yè)時間的影響[13]。對于存在障礙物的區(qū)域，Latombe首先提出了梯形分割法，即通過和覆蓋方向平行的割線首先將區(qū)域分為多個含障礙物或不含障礙物的子區(qū)域，再對每個不含障礙物的子區(qū)域進行牛耕式覆蓋[15]。之后，Choset提出了Boustrophedon細胞分解方法[16]，該方法可將區(qū)域分割為面積更小的單元。之后，Huang提出了一種分解區(qū)域的割線可變的“基于線掃描的分解方法”，該方法可將區(qū)域分為數(shù)量更少的子區(qū)域[17]。在此基礎(chǔ)上，陳海等考慮了無人機的轉(zhuǎn)向時間較長，通過最小化覆蓋作業(yè)的轉(zhuǎn)彎次數(shù)來確定覆蓋方向[18]。此外，Nedjati等首先將整個區(qū)域分為需要覆蓋、不需要覆蓋和無人機準備起飛的區(qū)域，之后再將整個區(qū)域劃分為均勻網(wǎng)格，最終求解出一條覆蓋所有不含障礙物的網(wǎng)格路徑[19]。

在覆蓋作業(yè)中，傳統(tǒng)的電動無人機存在續(xù)航時間短的缺點，針對該問題，本文提出了應用多架太陽能無人機的覆蓋路徑規(guī)劃問題。之后，基于文獻[13]中針對凸多邊形區(qū)域的覆蓋路徑規(guī)劃模型，提出針對凹多邊形以及內(nèi)部含障礙物的待覆蓋區(qū)域、以總作業(yè)時間和能量利用情況作為優(yōu)化目標、基于無向圖搜索的覆蓋路徑優(yōu)化模型，通過定義約束方程來求解每架無人機的最優(yōu)或相對最優(yōu)的覆蓋飛行方向和飛行路徑。在實際應用無人機進行覆蓋作業(yè)的任務中，一般以總作業(yè)時間為優(yōu)化目標。例如在應用無人機植保作業(yè)中，越短的作業(yè)時間意味著越高的作業(yè)效率，同時意味著更高的經(jīng)濟效益。而在應用太陽能無人機作業(yè)中，還需要考慮對太陽能的利用情況，當無人機機載設備耗能較大時，需要同時考慮能量利用情況和總作業(yè)時間這2個優(yōu)化目標，通過2個優(yōu)化目標在具體任務中的重要性占比來確定各自的優(yōu)化權(quán)重。本文創(chuàng)新點可總結(jié)如下：

1) 針對邊界存在障礙物的凹多邊形區(qū)域或內(nèi)部含障礙物的多邊形區(qū)域，提出基于無向圖搜索的覆蓋路徑優(yōu)化模型，以總作業(yè)完成時間最短為優(yōu)化目標，通過定義約束條件來求解無人機的最優(yōu)飛行路徑。

2) 將總作業(yè)完成時間最短和能量流動效率最高同時作為優(yōu)化目標，建立雙目標優(yōu)化方程。基于總作業(yè)完成時間最短，通過有限遍歷的方式來求解在2個優(yōu)化目標的權(quán)重系數(shù)不同時，使作業(yè)時間和能流效率相對最優(yōu)的覆蓋飛行方向和飛行路徑。

1 太陽能無人機能量模型的建立

對于固定翼無人機來說，無人機機翼上方覆蓋有太陽能電池板，將吸收來的太陽能轉(zhuǎn)化為電能。應用太陽能無人機作業(yè)時，最重要的問題之一就是太陽能無人機的能量模型對作業(yè)效果的影響。本節(jié)在文獻[20]的基礎(chǔ)上建立了太陽能無人機的能量模型，進一步提出能量流動效率這一評價指標。

1.1 能量輸入模型

已有太陽輻射模型主要分為兩種：一種是考慮氣象條件的全球輻射模型；一種是不考慮天氣影響的晴空輻射模型。由于前一種模型的計算需要具備多種氣象參數(shù)、計算復雜且準確度低，因此本文選擇晴空輻射模型。在晴空輻射模型中，有ASHRAE、Badescu、Bashasu等模型，其區(qū)別在于模型的輸入?yún)?shù)不同。其中，ASHRAE晴空輻射模型的相關(guān)參數(shù)齊全、計算相對簡單且精度較高，該模型最受全球研究工作者的認可[20]。因此本文選擇ASHRAE晴空輻射模型來建立太陽能無人機的能量輸入估計模型。

太陽能電池板表面吸收的太陽輻射強度I與太陽光線在機翼表面的入射角呈線性關(guān)系?？捎嬎愠龃怪庇谒矫娴奶栞椛?，即

(1)

式中：太陽常數(shù)I0=1 353(±21) W/m2[21]；nday為從1月1日開始算起的太陽天數(shù)。

進一步可以計算出太陽能電池板水平面的單位面積所吸收的太陽直接輻射強度Ib和太陽直接散射強度Id，如式(2)～式(4)所示:

(2)

(3)

(4)

式中：b和d分別為直射和散射空氣質(zhì)量指數(shù)；τb和τd分別為直射和散射光學深度，可通過查表獲得；mr為空氣質(zhì)量比；αe為太陽高度角。

建立能量輸入模型，首先應建立太陽光線與無人機之間的位置關(guān)系，如圖1所示。圖1中：αs表示太陽方位角，是目標物體的北方向與太陽光線入射方向的順時針夾角；αe表示太陽高度角，是地平面與太陽入射光線之間的夾角；λ表示太陽入射角，是太陽入射光線與翼面法線之間的夾角。由文獻[22]可計算出某一時刻和地點的αs和αe，如式(5)～式(8)所示：

圖1 太陽光線與無人機之間的位置關(guān)系

sinαe=sinnlatsinδ+cosnlatcosδcosω(t)

(5)

(6)

δ=0.409 3sin[2π(284+nday)/365]

(7)

ω(t)=0.261 8×(12-tlocal)

(8)

式中：nlat為作業(yè)地點的緯度；ω(t)為太陽時；δ為太陽赤緯角；tlocal為作業(yè)時刻。

進而計算太陽入射角λ，即

(9)

式中：Spb為太陽光線的單位向量在機體坐標系中的表示，即

(10)

式中：φ、θ、ψ分別為無人機的滾轉(zhuǎn)角、俯仰角和偏航角；R(φ)、R(θ)、R(ψ)分別為在這3個方向的導航坐標系到機體坐標系的轉(zhuǎn)移矩陣。

基于以上，可計算出太陽入射角的余弦，如式(11)所示：

cosλ=cosαesinαs(cosψsinφ-cosφ·

sinψsinθ)-cosαecosαs(sinφsinψ+

cosφcosψsinθ)+cosφsinαecosθ

(11)

根據(jù)ASHRAE晴空輻射模型，可計算單位機翼面積吸收的太陽輻射強度PS(λ)及太陽能無人機的能量輸入功率Pin，即

(12)

Pin=ηsolSPS

(13)

式中：S為電池板的面積；ηsol為其光電轉(zhuǎn)化效率。由此可計算出太陽能無人機在[t1,t2]時間內(nèi)吸收的總能量，即太陽能無人機的輸入能量Ein，即

(14)

1.2 能量輸出模型

無人機的輸出能量，主要指無人機克服空氣阻力做功和機載作業(yè)消耗。可計算出無人機定高飛行時的能量輸出功率Pout，即

(15)

式中：v為飛行速度；T為無人機螺旋槳拉力；ηprop為螺旋槳效率；PA為無人機機載作業(yè)消耗功率。由于無人機定高飛行時重力mg與升力L相等，因此可求出飛行速度v，如式(16)～式(18)所示：

Lcosφ=mg

(16)

(17)

(18)

式中：CL0為零升力系數(shù)；CLα為斜線升力系數(shù)；ρ為空氣密度；α為無人機的迎角。進而可計算出螺旋槳的拉力T，如式(19)～式(21)所示：

T=D

(19)

(20)

(21)

式中：Ra為無人機的展弦比；CD為飛行阻力系數(shù)；K為空氣動力系數(shù)；ε為奧斯瓦爾德效率因子。

基于式(1)～式(21)，可以計算出無人機在[t1,t2]時間內(nèi)的輸出能量Eout，即

(22)

選擇文獻[7]中的“Sky-Sailor”太陽能無人機作為本文的無人機模型，其主要參數(shù)如表1所示。

表1 太陽能無人機模型參數(shù)

1.3 能量流動效率

為了表示太陽能無人機在飛行作業(yè)過程中對太陽能的利用率，本文提出能量流動效率ηflow這一評價指標，即

(23)

由式(23)可知，能量流動效率是無人機飛行消耗和機載作業(yè)消耗的總能量與吸收的太陽能的比值。該值的大小體現(xiàn)了太陽能無人機是否及時有效地將吸收來的太陽能轉(zhuǎn)化為作業(yè)或飛行消耗所用的能量。由式(23)可知，能量流動效率ηflow必然小于1，并且該值越大，在飛行作業(yè)過程中對太陽能的利用率越高。此外，對于同一架太陽能無人機，能量流動效率越高，意味著電池的負載對于該作業(yè)來說越??；反而當能量流動效率過低時，電池的負載對于該作業(yè)任務來說則是一種“負擔”。因此在基于太陽能無人機的覆蓋路徑規(guī)劃問題中，可在以覆蓋作業(yè)時間為優(yōu)化目標的基礎(chǔ)上，將能量流動效率同時作為優(yōu)化目標來求解無人機的飛行路徑。

2 含障礙區(qū)域覆蓋路徑單目標優(yōu)化

在無人機覆蓋作業(yè)中，最主要的優(yōu)化目標是覆蓋作業(yè)時間。因此本節(jié)研究的問題是，有M架相同的無人機需要覆蓋一個邊界存在障礙物的凹多邊形區(qū)域或內(nèi)部含障礙物的多邊形區(qū)域。假設無人機的飛行高度恒定，不考慮續(xù)航時間對作業(yè)時間的限制。需要求解的是當以總作業(yè)完成時間最短為優(yōu)化目標時每架無人機的最優(yōu)路徑。為解決該問題，則需建立針對含障礙物區(qū)域的覆蓋路徑優(yōu)化模型。

2.1 覆蓋飛行方向的確定

本文基于牛耕法進行覆蓋路徑規(guī)劃，因此首先需要確定覆蓋飛行方向。為保證覆蓋作業(yè)的質(zhì)量，無人機常需要飛出區(qū)域進行轉(zhuǎn)彎，因此轉(zhuǎn)彎次數(shù)的增加往往對應著總路徑長度的增加[13]。為了減少無人機作業(yè)時間則應盡量減少無人機的轉(zhuǎn)彎次數(shù)。由于僅考慮總作業(yè)完成時間這一優(yōu)化目標，因此選擇能使轉(zhuǎn)彎次數(shù)最少的方向作為最佳覆蓋飛行方向。

首先定義一組坐標點來表示含障礙物待覆蓋區(qū)域的邊界，假設該組坐標點的集合為B={bk=(xk,yk),k=1,2,…,n}。進一步可以得到待覆蓋區(qū)域的邊集為L={lk=[(y-yk)(xk-xk+1)=(x-xk)(yk-yk+1),x∈(xk,xk+1),k=1,2,…,n]}。由于待覆蓋區(qū)域為凹多邊形或內(nèi)部含障礙物，因此將邊集L分為凹邊或內(nèi)部障礙物邊界集合L1={lk,k=1,2,…,h}和凸邊或外部非障礙物邊界集合L2={lk,k=1,2,…,n-h}，如圖2所示。

圖2 待覆蓋區(qū)域邊集的劃分

在凸多邊形區(qū)域中，垂直于多邊形最小寬度的方向是轉(zhuǎn)彎次數(shù)最少的方向[14]。而當待覆蓋區(qū)域存在凹邊或內(nèi)部含有障礙物時，由于無人機在凹邊處或障礙物邊界處也需要轉(zhuǎn)彎，因此轉(zhuǎn)彎次數(shù)不僅和待覆蓋區(qū)域?qū)挾扔嘘P(guān)，還與凹邊或障礙物的位置和形狀有關(guān)。

因此，針對含障礙物區(qū)域的最佳覆蓋飛行方向的選取原則如下。首先在無人機的覆蓋作業(yè)寬度確定之后，用一組間距為單架無人機覆蓋作業(yè)寬度的平行直線以不同的角度與待覆蓋區(qū)域相交，通過遍歷的方式找到相交節(jié)點數(shù)量最小的角度范圍。之后再在這組角度范圍中選取出待覆蓋區(qū)域最小寬度方向的垂直方向，作為無人機的最佳覆蓋飛行方向。其中，最小寬度方向是指間距最短的包容該多邊形的平行線的方向。覆蓋飛行方向的確定過程如圖3所示，其中無人機的最佳覆蓋飛行方向如圖中箭頭所示。

圖3 覆蓋飛行方向的確定

一旦覆蓋飛行方向被確定，就可以對待覆蓋區(qū)域進行行劃分，為后續(xù)的覆蓋任務分配做準備。這里假設在最佳覆蓋飛行方向的垂直方向上，待覆蓋區(qū)域的寬度為h，假設單架無人機的覆蓋作業(yè)寬度為L，進而可求得覆蓋作業(yè)的行數(shù)Nl，即

(24)

式中：為保證對待覆蓋區(qū)域的完全覆蓋，即兩個覆蓋行之間不存在未經(jīng)覆蓋的空白區(qū)域，因此對Nl進行向上取整。之后修正每一行的實際覆蓋寬度dl，即

(25)

2.2 覆蓋路徑優(yōu)化模型

在覆蓋作業(yè)的基本條件被確定之后，基于文獻[13]提出的針對凸多邊形待覆蓋區(qū)域的覆蓋路徑規(guī)劃模型，建立本文中針對含障礙物區(qū)域的基于無向圖搜索方法的覆蓋路徑優(yōu)化模型。首先，用一組平行且間距為dl的直線沿著最佳覆蓋飛行方向與待覆蓋區(qū)域邊界相交，得到一組交點P={pi=(xi,yi),i=1,2,…,N-1}。之后這組交點和無人機的起飛節(jié)點共同組成無向圖G=(V,E)中的節(jié)點集合V，如圖4所示。

圖4 無向圖節(jié)點的確定

之后，對該集合中的每個節(jié)點編號，所有節(jié)點的數(shù)量為N。其中無人機起飛節(jié)點的編號為1，與第1個覆蓋行相關(guān)的節(jié)點編號為2和3，與第2個覆蓋行相關(guān)的節(jié)點編號為4和5，以此類推。最終每個覆蓋行都與其對應的偶數(shù)和奇數(shù)節(jié)點相關(guān)。無向圖G=(V,E)中的邊集E由兩兩連接圖中的N個節(jié)點的所有線段組成，從而形成一個完整的無向圖。

在數(shù)學上，無向圖G可以用N×N的成本矩陣C來表示，其中元素Cij表示節(jié)點i和節(jié)點j之間的歐幾里得距離。無向圖G中的邊集表示為E={lij=[(y-yi)(xi-xj)=(x-xi)(yi-yj),x∈(xi,xj)],i,j=1,2,…,N}。假設凹邊及內(nèi)部障礙物區(qū)域由于存在障礙物而無法穿越，因此與凹邊及障礙物邊界相交的路徑將不被選取。為方便后續(xù)約束方程的定義，此處將成本矩陣中表示與凹邊相交的邊的元素置為零，即

if ? ?(xa,ya),

xa∈(xk,xk+1)∪(xi,xj),

(ya-yk)(xk-xk+1)=(xa-xk)(yk-yk+1),

xa∈(xk,xk+1),k=1,2,…,n;

(ya-yi)(xi-xj)=(xa-xi)(yi-yj),

xa∈(xi,xj),i,j=1,2,…,N,Cij=0

(26)

將無向圖中和凹邊及內(nèi)部障礙物邊界相交的邊去除后，可將針對凹多邊形或內(nèi)部含障礙物的多邊形區(qū)域的覆蓋路徑優(yōu)化模型的無向圖表示如圖5所示。

圖5 表示含障礙區(qū)域覆蓋路徑優(yōu)化問題的無向圖

在基于無向圖搜索的覆蓋路徑優(yōu)化模型中，位于起飛節(jié)點處的多架無人機的任務是遍歷劃分為行的待覆蓋區(qū)域，進一步可以將該任務描述為每架無人機須對無向圖中的邊進行搜索，保證其沿著覆蓋飛行方向訪問無向圖中的所有節(jié)點，即必須訪問無向圖中所有表示待覆蓋行的邊。假設覆蓋任務為無人機遙感作業(yè)，并且假設采用幾何遙感模型。根據(jù)相機的成像原理，單架無人機的遙感區(qū)域可建模為矩形或梯形[14]。因此，在本文中，單架無人機的覆蓋寬度即在覆蓋飛行方向的垂直方向上的遙感矩形的寬度。由上可知，在本文的覆蓋作業(yè)中，該覆蓋路徑優(yōu)化模型能保證無人機覆蓋每一個待覆蓋行，并在覆蓋行和邊界相交的節(jié)點處轉(zhuǎn)彎，因此能保證對待覆蓋區(qū)域的完全覆蓋。

根據(jù)優(yōu)化目標來定義目標方程，并通過定義約束條件來強制無人機在路線中按照一定規(guī)則飛行，最終保證每個待覆蓋行都有一架無人機“負責”。通過求解表示該優(yōu)化問題的方程組，最終可得到每架無人機訪問節(jié)點的最佳順序，即每架無人機的最優(yōu)路線，該求解過程如圖6所示。

圖6 針對含障礙區(qū)域的無向圖搜索方法的求解示意圖

基于以上變量的定義，可求出第k架無人機的作業(yè)時間Tk，即

(27)

在本節(jié)中，將總作業(yè)完成時間最短這一單一目標作為優(yōu)化目標，即飛行路徑最長的無人機的作業(yè)時間最短。進而，該優(yōu)化問題的目標方程可以表示為

min(max(Tk))

(28)

如式(28)所示，該目標函數(shù)可表示為一個線性的min-max問題。下面通過定義約束方程，來限制無人機按照一定的規(guī)則飛行。

約束1式(29)中的約束條件表示除起飛節(jié)點外的每個節(jié)點只能被一架無人機訪問一次。

(29)

約束2式(30)中的約束條件保證了到達某個節(jié)點的無人機和離開該節(jié)點的無人機是同一架。

p=2,3,…,N；k=1,2,…,m

(30)

約束3為了強制每架無人機在起飛節(jié)點1處起飛，并回到起飛節(jié)點，確保該路徑中沒有內(nèi)部子循環(huán)路徑，式(31)表示了標準的子循環(huán)消除約束。

ui,uj∈Z;i,j=2,3,…,N

(31)

式中：ui和uj代表任意實數(shù)[23]。

約束4為了保證無人機能夠按照最佳覆蓋方向飛行，確保每一個覆蓋行都能被無人機覆蓋，還需要定義約束方程：

i=2,4,…,N

(32)

式(32)的約束條件保證了訪問某一個覆蓋行，其中一個節(jié)點的無人機必將訪問該覆蓋行的另一個節(jié)點?？紤]到節(jié)點的編號順序，該約束條件通過強制訪問某個偶數(shù)節(jié)點的無人機必須訪問其對應的下一個奇數(shù)節(jié)點以及訪問某個奇數(shù)節(jié)點的無人機必須訪問其對應的上一個偶數(shù)節(jié)點來實現(xiàn)。因此，該約束對于使該方法成為解決覆蓋路徑優(yōu)化問題的可行方案來說是至關(guān)重要的。

約束5為了避免無人機沿著與覆蓋方向不平行的邊穿過待覆蓋區(qū)域，防止無人機急轉(zhuǎn)彎，則需要在優(yōu)化模型中添加兩個可選的約束，即

i=2,4,6,…,N

(33)

i=3,5,7,…,N

(34)

在針對含障礙區(qū)域的覆蓋路徑優(yōu)化模型中，為了使問題盡量有解，即路徑不會陷在凹邊或者障礙物邊界處，約束5有時不需要添加進來。如圖7，當添加了該約束條件后，無人機只能沿著與凹邊或障礙物邊界相交的邊l35、l37、l39飛行，因此問題就會變得無解。當不添加此約束時，無人機則可從節(jié)點3沿著l36、l38飛至節(jié)點6或8，從而飛出凹邊或障礙物邊界區(qū)域。當然無解只可能出現(xiàn)在待覆蓋行數(shù)量極少的情況下，而當待覆蓋行數(shù)量較多時，可選擇添加該約束條件。對于該約束條件對覆蓋路徑求解的影響將在2.3節(jié)中通過對比實驗來說明。

圖7 是否添加約束5對路徑求解的影響

約束6為了保證無人機不沿著穿越凹邊或障礙物邊界的路徑飛行，需要添加約束方程：

(35)

通過求解式(28)中的目標函數(shù)和式(29)～式(35)中的約束方程，即可求解出給定數(shù)量的無人機在最短的時間內(nèi)遍歷建模為無向圖G的待覆蓋區(qū)域的最優(yōu)路徑。

此外，目標函數(shù)(28)的求解目標是使得多架無人機中最長路徑的那架無人機的路徑最短，因此，一旦找到最長的無人機路徑，就無法保證其他無人機的路徑也最小化。為了保證所有無人機的路徑長度都最小，采用迭代的解決方案，通過多次迭代來保證每架無人機的路徑最優(yōu)。首先，在第1次迭代中，使用原始的無向圖G，解決的優(yōu)化問題正是前面提出的問題。在第2次迭代中，首先移除在第1次迭代中分配給最長路徑的無人機以及除起飛節(jié)點外所有屬于該路徑的節(jié)點和邊。更新參數(shù)之后，再用第1次迭代中的方法求解。以此類推，直至所有節(jié)點都被遍歷。通過迭代的方法，可以保證每架無人機的路徑都是最優(yōu)的，即在目標函數(shù)不變的前提下，最大可能地將覆蓋任務的總成本降到最低。該覆蓋路徑優(yōu)化模型的算法流程如圖8所示。

圖8 覆蓋路徑優(yōu)化模型的算法流程圖

2.3 仿真實驗

本文中的仿真實驗在運行內(nèi)存為4 GB RAM、CPU為Intel Core i3的計算機的MATLAB 2018b中進行，并通過Yalmip工具箱求解。

首先假設作業(yè)地點為(125°E,50°N)，作業(yè)時間為2020年5月1日上午09：00。假設無人機的覆蓋作業(yè)寬度為L=130 m。當無人機攻角α=2°時，可求得其速度為v=10.778 4 m/s。

針對凹多邊形區(qū)域和內(nèi)部含障礙物的多邊形區(qū)域，根據(jù)待覆蓋區(qū)域的形狀分別將實驗分為兩組，每組中假設可用的無人機數(shù)量m分別為2架、3架 兩種情況。仿真得到每種情況下的最佳覆蓋方向δ、修正后的覆蓋寬度dl、每架無人機的作業(yè)時間Tk以及總作業(yè)完成時間t。之后根據(jù)太陽能無人機的能量模型計算出其能量輸入功率Pin、能量輸出功率Pout，進而求出該覆蓋作業(yè)中的能量流動效率ηflow，用來評價太陽能無人機在該覆蓋作業(yè)中對能量的利用率。

2.3.1 針對凹多邊形區(qū)域

在本節(jié)的第1組實驗中，定義待覆蓋作業(yè)區(qū)域邊界的坐標分別為(-2 555,-986)、(-2 770,-99)、(341,-49)、(345,-961)、(-833,-1 073)、(-1 754,-683)、(2 074,-895) m，無人機起飛節(jié)點的坐標為(-3 335,-289) m。

當無人機數(shù)量分別為m=2，3時，每架無人機的最優(yōu)路徑如圖9所示。

圖9 每架無人機的最優(yōu)路徑(凹多邊形區(qū)域，第1組實驗)

由仿真求解可得，在該組實驗下，無人機的最佳覆蓋飛行方向為δ=1.55 rad，修正后的覆蓋寬度為dl=125.96 m。當無人機數(shù)量不同時，每架無人機的作業(yè)時間和總作業(yè)完成時間如表2所示。

表2 每架無人機的作業(yè)時間和總作業(yè)完成時間(凹多邊形區(qū)域，第1組實驗)

在本節(jié)的第2組實驗中，定義待覆蓋作業(yè)區(qū)域邊界的坐標分別為(-2 707,-1 152)、(-2 702,-38)、 (440,119)、(553,-932)、(-279,-1 111)、(-969,-1 140)、(-986,-791)、(-1 592,-433)、 (-2 071,-475)、(-2 058,-1152) m，無人機起飛節(jié)點的坐標為(-3 334,-47) m。

當無人機數(shù)量分別為m=2，3時，每架無人機的最優(yōu)路徑如圖10所示。

圖10 第架無人機的最優(yōu)路徑(凹多邊形區(qū)域，第2組實驗)

由仿真求解可得，在該組實驗下，無人機的最佳覆蓋飛行方向為δ=1.52 rad，修正后的覆蓋寬度為dl=119.83 m。當無人機數(shù)量不同時，每架無人機的作業(yè)時間和總作業(yè)完成時間如表3所示。

表3 每架無人機的作業(yè)時間和總作業(yè)完成時間(凹多邊形區(qū)域，第2組實驗)

2.3.2 針對內(nèi)部含障礙物區(qū)域

在本節(jié)的第1組實驗中，定義待覆蓋作業(yè)區(qū)域外邊界的坐標分別為(-2 780,-1 495)、(-2 780,-305)、(-550,-305)、(-550,-1 495) m，定義待覆蓋區(qū)域內(nèi)障礙物邊界的坐標為(-2 283, -958)、(-2 210,-525)、(-1 838,-482)、 (-1 799,-910) m，無人機起飛節(jié)點的坐標為(-3 424,-300) m。

當無人機數(shù)量分別為m=2，3時，每架無人機的最優(yōu)路徑如圖11所示。

圖11 每架無人機的最優(yōu)路徑(含障礙物區(qū)域，第1組實驗)

由仿真求解可得，在該組實驗中，無人機的最佳覆蓋飛行方向為δ=1.57 rad，修正后的覆蓋寬度為dl=119 m。當無人機數(shù)量不同時，每架無人機的作業(yè)時間和總作業(yè)完成時間如表4所示。

表4 每架無人機的作業(yè)時間和總作業(yè)完成時間(含障礙物區(qū)域，第1組實驗)

在本節(jié)的第2組實驗中，定義待覆蓋作業(yè)區(qū)域邊界的坐標分別為(-2 204,-170)、(-533,-1 227)、(92,-453)、(11,-283)、(72,76)、(-1 691, 390) m，定義待覆蓋區(qū)域內(nèi)障礙物邊界的坐標為(-1 118,-642)、(-1 260,-626)、(-1 352,-497)、(-1 304,-364)、(-961, -299)、(-921,-396) m，無人機起飛節(jié)點的坐標為(641,-5) m。

當無人機數(shù)量分別為m=2，3時，每架無人機的最優(yōu)路徑如圖12所示。

圖12 每架無人機的最優(yōu)路徑(含障礙物區(qū)域,第2組實驗)

由仿真求解可得，在該組實驗中，無人機的最佳覆蓋飛行方向為δ=1.75 rad，修正后的覆蓋寬度為dl=124.38 m。當無人機數(shù)量不同時，每架無人機的作業(yè)時間和總作業(yè)完成時間如表5所示。

表5 每架無人機的作業(yè)時間和總作業(yè)完成時間(含障礙物區(qū)域,第2組實驗)

在針對含障礙物區(qū)域的多智能體覆蓋路徑規(guī)劃問題中，文獻[24]提出一種多農(nóng)機編隊的基于旋轉(zhuǎn)光束與改進概率路線圖的路徑覆蓋算法。在該方法中，3架農(nóng)機組成倒“V”形編隊，等效成一架農(nóng)機后進行覆蓋作業(yè)。針對含障礙物待覆蓋區(qū)域的第1組和第2組實驗，實驗條件相同，應用該算法的實驗結(jié)果如圖13所示。

由圖13可看出，采用編隊作業(yè)進行覆蓋路徑規(guī)劃時，路徑中多次出現(xiàn)重復飛行的路徑，導致總飛行路徑增加。該方法與本文方法在作業(yè)時間上的對比如表6所示。其中，t1表示應用文獻[24]中方法的作業(yè)時間，t2表示應用本文方法的作業(yè)時間，a表示[(t2-t1)/t1]×100%，即本文方法在作業(yè)時間上提高的百分比。

圖13 其他覆蓋路徑規(guī)劃算法的對比實驗結(jié)果

表6 在覆蓋作業(yè)時間上的對比實驗結(jié)果

由表6可以看出，本文方法在路徑規(guī)劃時間上較文獻[24]中的方法提高了10%以上，且文獻[24]中的方法沒有考慮智能體往返于出發(fā)點的路徑長度。因此由仿真實驗可看出，本文采用多無人機分散作業(yè)且基于無向圖搜索的覆蓋路徑優(yōu)化模型可行且較優(yōu)。

最后，假設本文中用于覆蓋作業(yè)的無人機的機載作業(yè)消耗功率為PA=60 W，可計算出每架太陽能無人機的能量輸入功率Pin和輸出功率Pout為

Pin=72.462 6 W,Pout=68.305 7 W

根據(jù)Pin和Pout以及每架無人機的作業(yè)時間可求得作業(yè)中每架無人機的能量生產(chǎn)總值和能量消耗總值。進一步可計算出每架無人機的能量流動效率為

2.3.3 添加約束5的對比實驗

當考慮添加約束5 時，無人機將不被允許沿著與覆蓋行不平行的邊飛行。此處選取內(nèi)部含障礙物區(qū)域的第1組實驗中無人機數(shù)量m=3的情況，分別添加和不添加約束5進行對比。圖14為添加該約束條件的實驗結(jié)果，表7為兩種情況下每架無人機的作業(yè)時間和總作業(yè)完成時間的對比實驗結(jié)果。

圖14 添加約束5的對比實驗

表7 添加約束5的對比實驗結(jié)果

通過圖11和圖14的對比可以看出，未添加該約束條件的無人機路徑中出現(xiàn)了3條與待覆蓋行不平行的路徑，而在添加該約束條件的實驗中，無人機均沿著與覆蓋行平行的方向飛行。根據(jù)表7 中的實驗結(jié)果可以看出，當未添加該約束條件時，每架無人機的飛行作業(yè)時間比較平均，總作業(yè)完成時間也較短。而當添加了該約束條件時，為了使問題有解，每架無人機的路徑長度不夠平均，總作業(yè)完成時間也較長。

3 含障礙物區(qū)域考慮能流效率的覆蓋路徑雙目標優(yōu)化

當考慮無人機的姿態(tài)變化時，由于無人機轉(zhuǎn)彎時的姿態(tài)變化會影響其能量輸入和能量輸出，進而影響能量流動效率的大小。因此在本節(jié)中針對含障礙物區(qū)域，考慮了無人機在飛行過程中的姿態(tài)變化，將總作業(yè)完成時間最短和總能量流動效率最高同時作為優(yōu)化目標，在優(yōu)化作業(yè)時間的基礎(chǔ)上，通過有限遍歷的方式找到使作業(yè)時間和能流效率相對最優(yōu)的覆蓋飛行方向和每架無人機的飛行路徑。

3.1 轉(zhuǎn)彎路徑細化策略

無人機在平飛狀態(tài)時姿態(tài)角保持不變，而在轉(zhuǎn)彎狀態(tài)時，需要調(diào)整滾轉(zhuǎn)角來使無人機按照一定的轉(zhuǎn)彎半徑實現(xiàn)轉(zhuǎn)彎飛行。若要考慮無人機在轉(zhuǎn)彎時的姿態(tài)變化對最優(yōu)路徑規(guī)劃的影響，則應首先確定無人機在不同情況下的轉(zhuǎn)彎策略。

在本問題中有兩種情況需要考慮：無人機在兩個覆蓋行之間轉(zhuǎn)換的轉(zhuǎn)彎路徑；無人機從起飛節(jié)點飛至第一個節(jié)點后或從最后一個節(jié)點返回至起飛節(jié)點前的轉(zhuǎn)彎路徑。在轉(zhuǎn)彎半徑確定的前提下，以總路徑長度最短為原則來確定最佳的轉(zhuǎn)彎路徑。首先討論無人機在兩個覆蓋行之間轉(zhuǎn)換時轉(zhuǎn)彎路徑的細化。

首先假設無人機轉(zhuǎn)彎細化之前的路徑如圖15 所示。其中兩節(jié)點之間的距離為dAB，兩覆蓋行之間的距離為d，無人機的轉(zhuǎn)彎半徑為R，圖中箭頭的方向表示其飛行方向。

圖15 無人機轉(zhuǎn)彎細化前的路徑

由文獻[25]可知，無人機的最佳轉(zhuǎn)彎策略需要根據(jù)其轉(zhuǎn)彎半徑R、兩行之間的間距d以及轉(zhuǎn)彎處的折線AB與覆蓋行L1之間的夾角β來進行分類討論，每種情況下最優(yōu)的路徑如圖16所示。

圖16 覆蓋行之間的最優(yōu)轉(zhuǎn)彎策略

之后考慮無人機從起飛節(jié)點飛至第一個節(jié)點后或從最后一個節(jié)點返回至起飛節(jié)點前的轉(zhuǎn)彎路徑細化策略。分別考慮起飛節(jié)點到第一個/最后一個節(jié)點之間的連線與待覆蓋行之間的角度為銳角和鈍角的兩種情況。為保證覆蓋作業(yè)的覆蓋率并盡量使路徑最短，設計最佳飛行路徑如圖17所示。

圖17 起點與覆蓋行之間的最優(yōu)轉(zhuǎn)彎策略

由文獻[25]，可計算出無人機轉(zhuǎn)彎半徑R，即

(36)

由文獻[26]可知，固定翼無人機滾轉(zhuǎn)角的安全范圍為-π/4～π/4 rad，即滾轉(zhuǎn)角最大φmax=π/4 rad。當滾轉(zhuǎn)角取得最大值時，R=Rmin。進而由式(36)可計算出最小轉(zhuǎn)彎半徑為Rmin=7.15 m。

3.2 覆蓋路徑優(yōu)化模型

無人機在轉(zhuǎn)彎時，滾轉(zhuǎn)角會發(fā)生變化，該角度會影響無人機的能量流動效率。當將能量流動效率同樣作為優(yōu)化目標時，不同的轉(zhuǎn)彎路徑占比則會影響能量流動效率的大小，而當選擇不同的覆蓋飛行方向時，轉(zhuǎn)彎路徑的占比將發(fā)生變化。在本節(jié)中，將以能量流動效率最高和總作業(yè)完成時間最短同時作為雙優(yōu)化目標來求解相對最優(yōu)的覆蓋飛行方向和每架無人機的飛行路徑。首先給出雙目標優(yōu)化方程為

k1+k2=1

(37)

式中：k1、k2為總作業(yè)完成時間和總能量流動效率在優(yōu)化方程中的權(quán)重系數(shù)?？筛鶕?jù)總作業(yè)完成時間和總能量流動效率在作業(yè)任務中的重要性占比來確定k1和k2的大小。例如當無人機數(shù)量較少且待覆蓋區(qū)域面積較大時，可使總作業(yè)時間的權(quán)重系數(shù)更大，取k1>0.5；當無人機機載耗能較大時，可使總能量流動效率的權(quán)重系數(shù)更大，取k2>0.5。

由于無人機在平飛和轉(zhuǎn)彎時的姿態(tài)不同，進而會導致其能量流動效率不同。假設無人機在轉(zhuǎn)彎狀態(tài)下的能量輸入和能量輸出功率為Pin1和Pout1；無人機在平飛狀態(tài)下的能量輸入和能量輸出功率為Pin2和Pout2；無人機轉(zhuǎn)彎飛行用時tt，非轉(zhuǎn)彎飛行用時tnt，則可計算出當考慮無人機姿態(tài)變化時的能量流動效率ηflow，即

(38)

在覆蓋路徑雙目標優(yōu)化模型中，首先選取一定的角度為單位角度將覆蓋飛行方向劃分為有限種情況。在每種情況中，以時間最短為優(yōu)化目標，求解出每架無人機的最優(yōu)路徑。之后細化無人機的轉(zhuǎn)彎路徑，計算出每種情況下的轉(zhuǎn)彎用時和非轉(zhuǎn)彎用時，進而計算出每種情況下的總作業(yè)完成時間和總能量流動效率。最后設置兩種優(yōu)化目標的權(quán)重系數(shù)，將每種情況下的總作業(yè)完成時間和總能量流動效率代入雙目標優(yōu)化方程中，求解出使作業(yè)時間和能流效率相對最優(yōu)的覆蓋飛行方向和飛行路徑。由于遍歷求解次數(shù)較少，因此這種通過有限遍歷來優(yōu)化的方式在工程上求解相對快速、可行性高。

3.3 仿真實驗

在仿真實驗中，以10°為單位角度，將覆蓋飛行方向分為18種情況，計算每種情況的作業(yè)完成時間和能量流動效率，設置不同的k1、k2值，根據(jù)式(37)得出在不同k1、k2值下相對最優(yōu)的覆蓋飛行方案。

無人機的速度為v=10.778 4 m/s，假設覆蓋作業(yè)寬度為L=130 m。代入式(36)中可計算出最小轉(zhuǎn)彎半徑為Rmin=11.85 m。

進而假設無人機在飛行中的實際轉(zhuǎn)彎半徑為R=60 m，由式(36)可計算出無人機轉(zhuǎn)彎時的滾轉(zhuǎn)角為φ=11.2°。之后可計算出無人機在兩種狀態(tài)下的能量輸入功率和能量輸出功率如下：

Pin1=83.518 3 W,Pin2=72.462 7 W

Pout1=85.073 3 W,Pout2=68.305 7 W

假設覆蓋作業(yè)地點為(125°E,50°N)，作業(yè)時間為2020年5月1日上午10：00，待覆蓋區(qū)域邊界和起飛節(jié)點條件與含障礙物區(qū)域覆蓋路徑優(yōu)化第2組實驗條件相同。此處假設可用無人機數(shù)量為m=2，且忽略起飛操作時間。可得當覆蓋飛行方向不同時的最優(yōu)無人機路徑，如圖18所示。

圖18 覆蓋飛行方向不同時的最佳路徑

且根據(jù)無人機轉(zhuǎn)彎路徑細化策略，可得每種情況下的轉(zhuǎn)彎路徑細化后的飛行路徑，如圖19所示。

圖19 覆蓋飛行方向不同時轉(zhuǎn)彎細化后的最佳路徑

根據(jù)仿真實驗可得在覆蓋飛行方向不同時的單架無人機的作業(yè)時間tk和總作業(yè)完成時間t，并根據(jù)細化后的無人機路徑得到總的轉(zhuǎn)彎路徑用時tt、非轉(zhuǎn)彎路徑用時tnt，將其代入式(38)中可得每種情況下的總能量流動效率ηflow。將以上實驗結(jié)果整理如表8所示。

由表8可知，當覆蓋飛行方向δ=110°時，總能量流動效率最大，為ηflowmax=95.84%，當覆蓋飛行方向δ=10°時，總作業(yè)完成時間最小，為tmin=19.36 min。將以上數(shù)據(jù)代入式(37)中，可得在不同覆蓋飛行方向下當k1取值不同時的總優(yōu)化目標y。

表8 不同覆蓋飛行方向下的實驗結(jié)果

代入式(37)可計算得，當無人機機載耗能較大、總能量流動效率在優(yōu)化問題中的權(quán)重更大，即當k1k2時，在權(quán)重系數(shù)k1=0.6,0.7,0.8,0.9的情況下，覆蓋飛行方向δ=30°，雙目標優(yōu)化值y最小，min(y)=1.002 4,1.001 8,1.001 2,1.000 6。

4 結(jié) 論

針對普通電動無人機在覆蓋作業(yè)中存在的續(xù)航時間短的問題，提出了應用太陽能無人機進行覆蓋作業(yè)，并提出一種基于無向圖搜索的覆蓋路徑優(yōu)化模型，主要工作總結(jié)如下：

1) 建立了應用于覆蓋作業(yè)的太陽能無人機的能量輸入模型和能量輸出模型，定義了能量流動效率這一評價指標來評價太陽能無人機在作業(yè)過程中對吸收的太陽能的利用率。

2) 針對邊界存在障礙物的凹多邊形區(qū)域或內(nèi)部含障礙物的多邊形區(qū)域，以轉(zhuǎn)彎次數(shù)最少為原則確定覆蓋飛行方向，以總作業(yè)完成時間最短為優(yōu)化目標，建立基于無向圖搜索的覆蓋路徑優(yōu)化模型，可求解出每架無人機的最優(yōu)飛行路徑。

3) 考慮了無人機轉(zhuǎn)彎時的姿態(tài)變化對能量流動效率的影響，同時以總作業(yè)完成時間最短和總能量流動效率最高為優(yōu)化目標，建立雙目標優(yōu)化方程，在優(yōu)化作業(yè)時間的基礎(chǔ)上，通過有限遍歷的方式找到使作業(yè)時間和能量流動效率相對最優(yōu)的覆蓋飛行方向和飛行路徑。

本文提出的覆蓋路徑優(yōu)化模型的適用性強，針對不同形狀的待覆蓋區(qū)域，能夠得到使優(yōu)化目標最優(yōu)或相對最優(yōu)的飛行路徑。該方法在工程應用上適用范圍廣、可行性強。