趙 越, 李言鋒
(上海理工大學管理學院, 上海 200093)
近年來,無人機在農業領域的應用日益廣泛,植保無人機憑借其靈活性,可以在農田中進行植物檢測、施肥、噴灑農藥等操作。 中國農業土地遼闊且分布不均,不僅有集中的大塊田地、也有較為分散的小塊農田,不只有單一作物種植、也存在種植作物分散的多樣化種植區[1],因此對多區域無人機航跡規劃進行研究具有實際應用價值,可以減少不必要的航行和能耗,降低任務成本。
針對單架無人機的路徑規劃問題,已經得到了廣泛的研究和探討。 武錦龍[2]和李靖等學者[3]在避開障礙物的同時對多區域噴灑路徑進行全局規劃。 當單架次無人機對多區域執行遍歷任務時,王嘉琪[4]和Liu 等學者[5]分別使用改進編碼的差分進化算法和基于蟻群二元迭代優化的動態遺傳算法對轉化后的TSP 問題進行求解。 也有路徑規劃研究以多目標來確定任務收益,如飛行總路程、多任務區的偵查時間[6]、目標摧毀概率[7]等。 當無人機被分配到一個區域執行任務時,提出了全區域覆蓋搜索路徑規劃問題[8-9],大多數研究方法基于無人機飛行航向對無人機工作區間進行劃分與路徑規劃。 多架無人機的路徑規劃也得到了研究。 杜芳芳等學者[10]等和張哲等學者[11]研究了多無人機以最小成本完成任務目標的MTSP 問題,并泛化出異構無人機執行不同任務的多車型VRP 問題[12]等。
上述研究方法大多是將區域看作一點,將路徑規劃問題轉化為旅行商問題或多車輛路徑問題,或者對單一區域內的路徑覆蓋問題展開研究,但在農業應用中,將區域轉化為一點并不實際,無人機在區域內的覆蓋路徑將直接影響區域間的飛行距離。 針對上述問題,本文從超級網絡的角度,對多作業地塊之間的調度問題以及地塊內覆蓋路徑選擇問題進行描述,建立更具一般化的數學模型,并設計出一套高效求解的算法,再與將區域轉化為一點的處理方法做結果對比,證明本方法的可行性。
假設現有若干不連續地塊,使用無人機對其執行噴灑任務,使得路徑從出發點開始遍歷所有地塊并回到出發點后無人機能耗成本最小。 將地塊間的飛行任務看作旅行商問題,將地塊內的飛行任務看作覆蓋路徑規劃問題,因此,無人機執行地塊噴灑任務為旅行商問題與覆蓋路徑規劃問題的結合。
機械覆蓋路徑方式主要有牛耕往復法和內外螺旋法,在相同大小地塊內作業時,內外螺旋法的飛行長度大于牛耕往復法,因此本文選擇牛耕往復法作為區域內無人機作業方式。 牛耕往復法在矩形區域作業方式有2 種,如圖1 所示。 圖1 中,a、b為地塊的長度與寬度,d為無人機作業直徑,圖1 (a)、(b)兩種方式的作業長度相同,均為但由于轉彎會產生一定的非作業長度,將其稱為“轉彎長度”,分別為可以看出(b)方式的轉彎長度更小,因此在區域內使用(b)方式遍歷路徑可使飛行能耗更小。 但是如果只追求區域內路徑長度最小,可能會導致區域間飛行長度或轉彎角度增加,考慮到地塊的飛入點與飛出點會直接影響區域間的飛行路徑,因此要獲得最優的多區域覆蓋路徑,應合理選擇地塊飛入點,使得整體飛行能耗最小。

圖1 無人機覆蓋路徑方式示意圖Fig. 1 UAV coverage path planning diagram
使用牛耕往復法對區域內進行遍歷,一個地塊內有8 個可行的飛入點,假設區域長度與寬度均為作業直徑的整數倍,那么8 個飛入點對應8 個可行的飛出點,將其進行編號,如圖2 所示。

圖2 區域進入點示意圖Fig. 2 Regional gateway diagram
超級網絡是指通過添加虛擬節點和弧段形成一個網絡拓撲結構,對研究對象的行為與特征進行抽象化描述的一種研究方法[13]。 在地塊規模較大時,利用超級網絡的虛擬節點與虛擬弧段代替無人機真實飛行路徑,使得問題與數據信息更加簡潔清晰。下面利用超級網絡對無人機噴灑農藥過程進行描述,如圖3 所示。 地塊飛入的8 個點對應了8 條不同的覆蓋路徑,將無人機在地塊內的8 條覆蓋路徑轉化為8 條虛擬弧段,該虛擬弧段包含區域內飛行距離、覆蓋線路起點坐標、終點坐標等已知信息。 地塊間的飛行路徑實際為地塊內覆蓋線路終點至下一地塊覆蓋線路起點的弧段,由于無人機實際以哪一條線路覆蓋地塊未知,因此將地塊看作一矩形點用虛擬弧段將其連接起來,該矩形點包含地塊排列序號與覆蓋線路序號的未知信息,該虛擬弧段包含地塊間飛行距離的未知信息。

圖3 超級網絡示意圖Fig. 3 Super network diagram
考慮研究問題的情景,對模型建立進行了如下假設:
(1)無人機在飛行時通過改變飛行方向和姿態來完成轉彎和變向操作,而不需要繞著一個固定的曲率半徑進行圓弧形的轉彎,即飛行軌跡為直線或折線;
(2)實施噴灑任務的地塊為矩形區域且不重疊;
(3)起點看作邊界點坐標和飛入飛出點坐標均為(0,0)、區域內各線路距離均為0 的地塊;
(4)無人機有足夠的動力完成該任務,并以恒定的高度飛行。目標函數為:
約束條件為:
模型中有2 類符號:模型參數和決策變量。 在式(1)~式(6)中,模型參數主要有:n為目標點總數,包括1 個出發點和n- 1 個地塊;為無人機在覆蓋地塊i的路線起點,共8 條線路;是無人機在覆蓋地塊i的路線終點,共8 條線路;為無人機以第k條路線覆蓋地塊i的終點與無人機以第l條路線覆蓋地塊j的起點的歐式距離;為無人機以第k條線路覆蓋地塊i時的飛行距離。 決策變量主要有:xij表示無人機是否將從地塊i移動到地塊j,如果是,xij=1,否則xij=0;表示如果無人機以第k條路線覆蓋地塊i,則=1,否則=0;zjl為如果無人機以第l條路線覆蓋地塊j,則=1,否則=0。
具體來說,式(1)的目標函數為最小化飛行距離;式(2)、式(3)表示每個地塊只能進入一次;式(4)、式(5)表示每個地塊只能以一種線路覆蓋;式(6)的約束表示消除子環路、但不消除最大的環路。
遺傳算法是一種啟發式優化算法,通過模擬自然界中的進化過程來解決問題。 將問題的解表示為某個個體的基因型,通過選擇、交叉和變異等基因操作,產生新的解,并逐步優化當前的最優解。 現已廣泛應用于優化問題、機器學習、圖像處理等領域。
本文采用順序表示的編碼方式,每一個數字代表一個地塊,遍歷所有的地塊產生一條染色體、即產生了1 個解,以6 個地塊1 個起點為例,地塊編碼是2~7 的數字,若生成的染色體是1342756,表示從起點1,依次經過地塊3、4、2、7、5、6,最終回到起點1的一條路徑回路。
假設由n個地塊構成的染色體為|C1|C2|…|Cn-1|Cn |,個體的適應度值為f=-F,其中F為從起點出發遍歷所有地塊回到起點所產生的距離與地塊內作業距離的總和。 適應度越大表示染色體越優,反之越劣。 定義了一個計算目標函數值F的函數CalObj,輸入參數為區域信息(包含地塊編號、地塊內線路編號、地塊內遍歷長度)以及出發點坐標,將地塊內線路編號序列映射到實際的路徑上,如Xnow={2,4,6,3,2,8},假設最大地塊數為6,那么進行6 次for 循環,將Xnow的第一個元素對應的區域的信息存儲在一個新的矩陣中,其中新矩陣已包含起點坐標,從Xnow數組中刪除第一個元素,變為5 個元素,再進行下一次循環,根據新生成矩陣中節點之間的距離計算目標函數值F。
交叉操作采用“子路徑互換交叉”的交叉方法,確定待交叉的2 個父代個體,從變量列數中隨機生成2 個位置,確定需要進行交叉的位置區間,即2 個隨機位置之間的所有位置。 將父代1 中位置區間的值與父代2 中相同的值進行交換,將父代2 中的位置區間的值與父代1 中相同的值進行交換。 變異操作采用的是2-opt 組合優化鄰域算子,是一種基于交換2 個城市路徑中間的路徑段來產生新解的變異方法。 從待變異的個體中隨機選擇2 個位置,將這2 個位置之間的路徑段反轉,即路徑上這2 個地塊之間的所有地塊的順序反轉,得到一個新的路徑。
變鄰域深度搜索實現了5 種不同類型的鄰域操作,包括交換鄰域(swap)、插入鄰域(insert)、2-opt鄰域、or-opt 鄰域和反轉加上插入鄰域。 對于每種鄰域操作,該函數都會在當前解的鄰域中進行搜索,生成一個新解,并計算新解的目標函數值。 如果新解的目標函數值比當前最優解更優,則更新最優解,并記錄產生最優解的鄰域操作的具體參數,即交換、插入、反轉等操作的具體位置。
步驟1種群初始化。 根據地塊順序對該問題進行編碼形成初始化種群。
步驟2適應度求解與選擇。 根據染色體基因的排列得到初始解的目標函數值,并計算種群中每個染色體的適應度以及種群中最佳的適應度,同時使用錦標賽選擇法對隨機選取的2 個個體進行比較,選取適應度值最優的個體作為本次選擇的結果。
步驟3染色體交叉操作。 首先通過錦標賽從當前種群中選擇出一組父代個體,然后通過交叉概率來決定是否進行交叉操作。 如果需要進行交叉操作,將這2 個父代個體進行子路徑互換交叉,生成新的子代。
步驟4染色體變異操作。 由交叉操作生成的染色體通過變異概率來決定是否進行變異操作,使用2-opt 交換方法進行變異,從而得到一個新的基因序列。
步驟5局部搜索操作。 將交叉變異后的子代個體作為下一代種群并進行局部搜索操作,得到更新后的解。
重復步驟2~5,直至滿足終止條件,即可輸出當前的迭代次數和無人機地塊訪問順序。
實驗在Matlab R2022a 環境下構建基于遺傳算法的植保無人機路徑規劃模型。 本文實驗地塊數據來源于某種植區地圖信息。 首先選用6 塊區域,區域A、B、C、D、E、F的邊界點坐標與中點坐標見表1,無人機噴灑直徑為5 m。

表1 各地塊邊界點坐標與中點坐標Tab. 1 Coordinates of boundary points and midpoints of each block
首先使用將區域轉化為地塊中點的方法對該案例進行求解,設交叉率crossover_prob=0.6,變異率mutation_prob=0.2,種群大小pop=50,當進化到第3 代時得到最優路徑為A-C-B-D-F-E,確定了地塊順序后,使用Matlab 內置遺傳算法求得地塊內線路序號為6-5-4-2-7-2,得到距離最優值為1 065.4 m。
如果考慮地塊內飛行路線對區域間飛行距離的影響對該案例進行求解,設交叉率crossover_prob=0.6,變異率mutation_prob=0.2,種群大小pop=50,迭代50 次,進行10 次實驗,記錄每次實驗下的區域順序、線路序號、距離最優值等結果,并得到最大值、最小值和方差,見表2、表3。

表2 迭代50 次實驗結果Tab. 2 The experimental results after 50 iterations

表3 各地塊邊界點坐標與中點坐標Tab. 3 Coordinates of boundary points and midpoints of each block
根據實驗結果可看到, 當區域順序為A-B-D-F-C-E,對應線路序號為6-6-1-7-2-2 時,飛行距離最小值為935.15 m,與將區域轉化為地塊中點相比飛行距離減少12.2%。 達到飛行距離最小值的次數為4 次,方差較小,結果更加穩定。 迭代50 次進化過程如圖4 所示,對模型進行50 代循環求解時,可以看出算法收斂較快,下降趨勢明顯,第8 代求得疑似最優解。 航跡規劃示意如圖5 所示。圖5 中,藍色線為無人機在區域內作業路徑,紅色線為區域間飛行路徑,綠色點表示地塊飛入點,黃色點表示地塊飛出點,箭頭為飛行方向。

圖4 迭代50 次進化過程圖Fig. 4 The evolutionary process diagram for 50 iterations

圖5 航跡規劃示意圖Fig. 5 Schematic diagram of path planning
當地塊距離較遠時,同樣選取6 塊區域進行求解,區域A、B、C、D、E、F的邊界點坐標與中點坐標見表3。
若使用將區域轉化為地塊中點的方法對該案例進行求解,當進化到第3 代時得到最優路徑為A-B-C-F-D-E,確定了地塊順序后,使用Matlab內置遺傳算法求得地塊內線路序號為5-2-2-7-8-8,得到距離最小值為10 291 m。 若考慮地塊內飛行路線對區域間飛行距離的影響對該案例進行求解,進行10 次實驗,當區域順序為E-A-D-B-C-F,對應線路序號為5-6-5-6-1-1 時,得到飛行距離最小值為9 673.6 m,與將區域轉化為地塊中點相比飛行距離減少6.0%,因此當地塊距離相距較近時,地塊飛入飛出點對區域間飛行距離的影響更大。
為了驗證該方法在解決更大規模問題的能力,創建了20 個矩形區域。 進行100 代循環求解,種群大小為200,記錄10 次實驗下的區域順序、線路序號、距離最優值、最小值和方差等結果,見表4。

表4 迭代100 次實驗結果Tab. 4 The experimental results after 100 iterations
根據實驗結果可知,當區域順序為7-16-20-13-18-1-4-2-14-5-3-17-15-6-8-10-12-9-11-19,對應線路序號為7-7-3-8-4-6-5-4-2-5-3-8-1-1-2-1-7-2-2-4 時,飛行距離最小值為3 331.49 m。飛行航跡規劃示意圖如圖6 所示(地塊內飛行線路省略)。 對模型進行100 代循環求解時,可以看出算法收斂較快,下降趨勢明顯,第50 代求得疑似最優解,如圖7 所示。

圖6 較大規模的區域航跡規劃示意圖Fig. 6 Large-scale area trajectory planning diagram

圖7 迭代100 次進化過程圖Fig. 7 The evolutionary process diagram for 100 iterations
本文針對無人機農藥噴灑問題,考慮了地塊飛入飛出點對總飛行距離的影響,建立了一般化的區域間覆蓋路徑規劃模型,將無人機多地塊農藥噴灑問題轉化為旅行商問題與覆蓋路徑規劃問題的結合。 使用遺傳算法對2 種區域規模的算例進行求解,驗證了該方法的可行性,得到了最佳迭代次數下最優飛行距離。 未來將對不規則地塊進行覆蓋路徑規劃研究。