鐘智超 肖雄武 涂建光

(1. 武漢大學(xué) 遙感信息工程學(xué)院, 湖北 武漢 430079;2. 武漢大學(xué) 測繪遙感信息工程國家重點實驗室, 湖北 武漢 430079)

0 引言

隨著傳感器與控制技術(shù)不斷發(fā)展,無人機已經(jīng)成為測繪遙感領(lǐng)域獲取數(shù)據(jù)的重要工具[1-2]。目前常見的無人機航線規(guī)劃算法致力于求解旅行客問題[3](traveling salesman problem,TSP)。常用的求解法有動態(tài)規(guī)劃[4]、遺傳算法、蟻群算法、概率模型估計[5]等,但這類算法很少考慮具體的飛行任務(wù),例如，建筑物測圖[6]、道路勘探、電路巡檢[7]、農(nóng)業(yè)植保[8]等。另外,無人機在建筑物測圖領(lǐng)域的應(yīng)用場景越來越多[9]。例如，在城市規(guī)劃中,無人機結(jié)合傾斜攝影測量技術(shù)可以實現(xiàn)城市實景三維模型的構(gòu)建;在地圖生產(chǎn)中,無人機高精度的航片可以用于正射影像生產(chǎn)、數(shù)字線劃圖生產(chǎn)[10]等。

針對常見航線規(guī)劃算法對建筑物測圖不具備針對性的現(xiàn)狀,提出一種針對建筑物測圖的無人機航線自動規(guī)劃算法。針對建筑物與植被粘連的問題,結(jié)合影像的光譜特性分離植被,縮小作業(yè)范圍;設(shè)計基于無人機多視圖立體幾何的地表地勢的估計方法,分析出飛行區(qū)域內(nèi)粗略的高區(qū)和低區(qū);結(jié)合植被信息與地勢信息完成對建筑物區(qū)域的分割;設(shè)計基于邊緣檢測的建筑物主方向識別,優(yōu)化航線角度;最后利用霍普菲爾德網(wǎng)絡(luò)求解覆蓋多個建筑物區(qū)域的最優(yōu)路徑。

1 測區(qū)內(nèi)建筑物區(qū)域定位

建筑物測圖任務(wù)中,研究對象是建筑物，不包含其他地物。而常規(guī)的無人機規(guī)劃線路往往對整個測區(qū)進行“條帶狀”的航帶劃分,如圖1所示,這會導(dǎo)致大量的無用航點和無用航片,尤其是在建筑物集群分布而且群體之間的距離較遠的情況。

圖1 “條帶狀”航線和常見“無效航片”示例

1.1 植被區(qū)域分離

植被是建筑物測圖中“無效航片”的重要來源,建筑物集群分布的特性導(dǎo)致存在大量僅被植被覆蓋的區(qū)域,可以利用植被特有的光譜特性來區(qū)分其與別的地物類型,如式(1)所示。

(1)

式中,(x,y)為像素;NIR(x,y)為近紅外波段值;R(x,y)為紅色波段的灰度值。歸一化植被系數(shù)(normalized differential vegetation index,NDVI)可以用于區(qū)分常見的綠色植物。綠色植物對近紅外波段的反射率較高,對可見光紅色波段的吸收率較高,因此將NDVI(x,y)大于閾值的像素標(biāo)記為植被。

1.2 地表地勢估計

除了植被之外,“裸地”也是非建筑物地物中占地面積較大的部分,尤其對比較干旱的測區(qū),“裸地”的影響程度更大。區(qū)分建筑物與“裸地”的重要方式是通過地勢高低差來判斷。在無人機作業(yè)中,很難直接獲取精確的地勢高程數(shù)據(jù),如數(shù)字表面模型(digital surface model,DSM)等,因此可依靠無人機多視角測量來估計地表[11]。

1.2.1無人機多視立體幾何

第一步,利用無人機在起飛點附近俯拍整個測區(qū)(根據(jù)測區(qū)大小調(diào)整飛行高度),記錄拍照瞬間無人機傳感器的位姿信息,具體包括外參矩陣R1以及根據(jù)全球定位系統(tǒng)(global positioning system,GPS)定位獲取的攝影中心初始位置信息,用矩陣t1表示。值得注意的是,初始時姿態(tài)角可以直接利用無人機萬向角信息進行初始化。

第二步,移動無人機并調(diào)整鏡頭位姿,在另一位置拍攝,同樣覆蓋整體測區(qū)。其間通過無人機慣導(dǎo)數(shù)據(jù)結(jié)合初始化的位姿推算第二次拍攝瞬間的位姿R2和t2。

第三步,根據(jù)灰度特征匹配兩張影像,計算出匹配點對的深度,從而獲取對應(yīng)3D點的估計坐標(biāo)。

圖2中,地表上不同地物同時被兩次拍攝記錄到,對于地表的一點P(X,Y,Z),其對應(yīng)兩張影像上的兩個像素點p1(u1,v1)和p2(u2,v2),在已知兩次攝影無人機位姿的估計值,那么有式(2)的關(guān)系。

(2)

式中,K為相機的內(nèi)參矩陣，由于未知量是地表點P(X,Y,Z)坐標(biāo),因此至少需要兩組同名點對才能求解該問題,因此需要無人機在不同的位姿對同一地表進行拍攝。

圖2 利用無人機多視立體幾何估算地表示意圖

1.2.2估計地形表面

通過無人機多視立體幾何可以求解地面點的三維坐標(biāo)[12],因為位姿不夠精確,所以這種方式求解的地面點是不夠精確的。但是這一步只需要對地表地勢有一個相對估計,即能夠區(qū)分地表上哪些位置地勢更高(如建筑物、樹木等),哪些位置地勢更低(如裸地、坑洼等)[13],因為地表本身是有坡度的,首先需要對地表面進行擬合。建筑物或樹木產(chǎn)生的高度差是突變的,而地表面是緩慢地升高或降低,因此可以設(shè)定一個權(quán)重函數(shù)限制這種突變,如Kraus和Pfeifer提出的分段式權(quán)重函數(shù),如圖3所示,根據(jù)高度變化殘差的大小來重新加權(quán)計算高程[14]。

圖3 用于估計地形表面的權(quán)重函數(shù)

圖3中,p(r)為權(quán)重關(guān)于殘差變化的函數(shù);r為殘差大小;g和t為分段權(quán)重的殘差閾值。首先用等權(quán)重來加權(quán)這些點的坐標(biāo)Z值,擬合一個初始的地表面,然后逐點比較坐標(biāo)Z值與擬合表面的估計值的差異,根據(jù)差異的大小來重新加權(quán)Z值。該分段權(quán)重函數(shù)把殘差大于閾值t的點直接設(shè)置為“外點”,權(quán)重設(shè)置為0。

把分段擬合后的地表面記為GS,那么區(qū)分某點P(X,Y,Z)是否屬于建筑物的方法如式(3)所示。其中(x,y)表示與點P(X,Y,Z)對應(yīng)的像素坐標(biāo)。

(3)

2 單區(qū)域航線規(guī)劃與全局最優(yōu)求解

2.1 單區(qū)域航線規(guī)劃

邊緣特征描述了像素之間的梯度關(guān)系。數(shù)學(xué)上,用一階導(dǎo)數(shù)來表達像素梯度被稱為一階微分算子,如Sobel算子、Canny算子等。利用二階導(dǎo)數(shù)來檢測梯度的算子被稱為二階微分算子,如LoG算子等。為了提高建筑物邊界檢測結(jié)果的完整性,防止出現(xiàn)對不明顯的邊界的漏檢,設(shè)計了一種8方向的增強一階微分算子,算子的結(jié)構(gòu)如圖4所示。該算子對比普通的一階算子的改進包括,第一,從普通一階算子只檢測2個方向的梯度變化,拓展到檢測8個方向;第二,擴大了檢測算子的鄰域大小,從3×3拓展到5×5,增大梯度檢測的感知范圍。

圖4 8方向增強邊緣檢測算子

為了驗證該改進算子的提取能力,設(shè)計了與經(jīng)典的Canny算子、LoG算子的對比實驗,驗證結(jié)果如圖5所示?？梢钥闯?虛線框中目標(biāo)建筑物被陰影遮擋的邊緣梯度很小,提取難度較大,該8方向增強算子能保持較好的檢出率。

圖5 8方向增強算子提取建筑物邊緣的對比結(jié)果

建筑物的主方向可以表示該建筑物的朝向,假設(shè)用方向向量d表示主方向,首先利用Hough變換對2.1中獲取的建筑物邊緣進行線段化,得到所有表示建筑物邊界的直線段集合L。根據(jù)直線段判斷建筑物區(qū)域的主方向步驟如下,如圖6所示。

圖6 確定主方向與單區(qū)域航線規(guī)劃

第1步,對于任意L中的每一條直線段,先計算其長度,對于長度小于閾值(根據(jù)經(jīng)驗一般設(shè)定為真實尺寸1.5 m)的設(shè)定為噪聲,對于長度大于閾值的保留為有效線段。

第2步,對于每一條有效線段,計算其方向向量。按照每30度為一個區(qū)間劃分6個方向角區(qū)間,根據(jù)方向向量將這些有效線段分到這6個區(qū)間內(nèi)。

第3步,統(tǒng)計每個區(qū)間內(nèi)有效線段的總長度,挑選出總長度最長的一個區(qū)間,將區(qū)間內(nèi)的所有線段的方向向量按照長度進行加權(quán)求和,得到該區(qū)域內(nèi)建筑物的主方向。

2.2 全局最優(yōu)航線求解

無人機的單次任務(wù)是從起飛到降落的全過程,需在盡可能短的時間內(nèi),完成對所有建筑物區(qū)域的飛行任務(wù),即找到從起飛點到完成所有建筑物拍攝并返回起飛點的最短路徑。Hopfield神經(jīng)網(wǎng)絡(luò)是一種結(jié)合存儲系統(tǒng)和二元系統(tǒng)的神經(jīng)網(wǎng)絡(luò)?；谶B續(xù)Hopfield的特點,如果把一個最優(yōu)化問題的目標(biāo)函數(shù)轉(zhuǎn)換成網(wǎng)絡(luò)的能量函數(shù),把問題的變量對應(yīng)于網(wǎng)絡(luò)的狀態(tài),那么Hopfield神經(jīng)網(wǎng)絡(luò)就能夠用于解決優(yōu)化組合問題[15]。

2.2.1能量函數(shù)構(gòu)建

首先構(gòu)建無人機飛行的Hopfield次序矩陣e,次序矩陣中的元素eij代表無人機在第i時刻執(zhí)行j建筑物區(qū)域的飛行任務(wù)。假設(shè)所有的目標(biāo)建筑物區(qū)域總數(shù)為N,那么將產(chǎn)生如下約束條件。約束1,次序矩陣的每一行只有一個神經(jīng)元被激活,即每個區(qū)域有且只有1次飛行任務(wù);約束2,次序矩陣的每一列只有一個神經(jīng)元被激活,即無人機不能同時執(zhí)行2個或2個以上區(qū)域的任務(wù);約束3,次序矩陣所有激活的神經(jīng)元數(shù)量之和等于所有的建筑物區(qū)域數(shù)量總和,即無人機需要執(zhí)行完所有建筑物區(qū)域的飛行任務(wù)。三個約束條件如式(4)所示。

(4)

無人機整體測區(qū)路徑優(yōu)化的目標(biāo)是所有路徑長度加起來最小,如式(5)所示,其中d表示距離矩陣,dik代表第i個飛行區(qū)域到第k個飛行區(qū)域的距離,默認在初始時對這些待飛區(qū)域進行了編號。

(5)

因此,將三個約束條件以及最小化總路徑長度的目標(biāo)函數(shù)轉(zhuǎn)化為能量函數(shù)的懲罰項,可以得到總體能量函數(shù)的表達式，如式(6)所示。

(6)

式中,前3項是約束項；第4項為目標(biāo)函數(shù)；A、B、C、D為懲罰參數(shù),為固定值,這些值反映了這幾項在最小化過程中的相對重要性。

然后將式(6)描述的航線優(yōu)化問題按照Hopfield網(wǎng)絡(luò)標(biāo)準(zhǔn)能量函數(shù)對應(yīng)起來,其標(biāo)準(zhǔn)形式如式(7)所示。

(7)

式中,權(quán)重和偏重Wijkl、Iij展開表示如式(8)所示。

(8)

式中,δij代表指示函數(shù),其在i=j時取值為1,i≠j時取值為0。

2.2.2迭代求解最優(yōu)航線

確定航線優(yōu)化問題的約束條件和標(biāo)準(zhǔn)能量函數(shù)之后,需要不斷迭代修正權(quán)重矩陣,直到狀態(tài)趨于穩(wěn)定,具體的步驟如下。

第1步,設(shè)定懲罰系數(shù)A、B、C、D和最大迭代次數(shù)Tmax,隨機初始化次序矩陣e。

第2步,計算各個建筑物區(qū)域之間的距離,形成距離矩陣d。

第3步,根據(jù)初始化的輸入狀態(tài)U0,確定狀態(tài)轉(zhuǎn)移函數(shù),如式(9)所示。

(9)

第4步,計算輸入狀態(tài)的增量,如式(10)所示。

(10)

第5步,利用動態(tài)方程和一階歐拉法更新下一刻的輸入狀態(tài)和輸出狀態(tài),如式(11)所示。

(11)

第6步,檢查當(dāng)前的輸出狀態(tài)e是否滿足約束,判斷迭代是否結(jié)束,否則返回第4步,直到迭代次數(shù)達到設(shè)定的最大次數(shù)。

3 實驗分析

為了驗證該航線規(guī)劃算法的有效性,利用國際攝影測量與遙感學(xué)會(International Society for Photogrammetry and Remote Sensing,ISPRS)提供的公開數(shù)據(jù)進行模擬實驗,該地區(qū)的建筑物占地面為4 099.5 m2。該地區(qū)建筑物面積較大,分散分布且樹木環(huán)繞,首先對其中的建筑物區(qū)域進行定位,結(jié)果如圖7所示。

圖7 建筑物區(qū)域定位實驗結(jié)果

如圖7(a)所示,目標(biāo)測區(qū)內(nèi)建筑物分散分布且在朝向上具有明顯差異;圖7(b)中亮區(qū)代表植被,暗區(qū)代表非植被;圖7(c)中亮區(qū)代表地勢高的區(qū)域,暗區(qū)代表低地勢;結(jié)合圖7(b)和圖7(c)可以定位具體的建筑物區(qū)域,利用增強8方向的檢測算子提取目標(biāo)建筑物的邊界,如圖7(d)所示。

根據(jù)建筑物區(qū)域的定位結(jié)果與邊緣提取的結(jié)果,然后計算每個區(qū)域內(nèi)建筑物的主方向,根據(jù)主方向確定每個小區(qū)域內(nèi)的航線。最后將所有小區(qū)域的航線全部生成后,通過全局航線優(yōu)化確定無人機執(zhí)行各區(qū)域飛行任務(wù)的順序。在初始化Hopfield距離矩陣時需考慮到建筑物區(qū)域本身具有的空間寬度,距離是兩個區(qū)域最小外接矩形幾何中心之間的距離。確定初始距離矩陣之后,設(shè)置初始系數(shù)A=B=D=50,C=20進行Hopfield網(wǎng)絡(luò)的迭代趨近,結(jié)果如圖8(b)所示,輸出狀態(tài)在迭代20次左右逐漸穩(wěn)定,最小損失穩(wěn)定在315.654,并得到如圖8(a)所示的整體航線優(yōu)化結(jié)果。

圖8 全局航線優(yōu)化結(jié)果與Hopfield網(wǎng)絡(luò)迭代收斂

為了驗證該方法規(guī)劃的航線與不考慮建筑物分布的“條帶狀”航線的效率差異,設(shè)計了同區(qū)域的對比實驗。設(shè)置相同的重疊率使得航帶間距為12 m,固定無人機飛行速度為6 m/s,忽略無人機拍照時間與等待,生成的航線結(jié)果如圖9所示。為了合理地評價無人機執(zhí)行飛行任務(wù)的效率,這里給出一些指標(biāo)。

有效航線距離(d1)：屬于航線規(guī)劃的航線,在建筑物區(qū)域(參考圖8(a)所示的建筑物區(qū)域定義)內(nèi)的航線長度。

通勤距離(d2)：從當(dāng)前航線飛往下一航線的中間過程飛行距離。

無效航線距離(d3)：屬于航線規(guī)劃的航線,但是不屬于建筑物區(qū)域內(nèi)部,在無效航線上拍攝的航片稱為無效航片。

航線效率(r1)：有效航線距離與航線規(guī)劃的總長度的比值(不包括通勤距離)。

作業(yè)效率(r2)：有效航線距離與無人機作業(yè)全流程飛行總距離的比值(包括通勤距離)。

圖9 兩種方法生成的航線結(jié)果

根據(jù)模擬實驗數(shù)據(jù)得到如表1所示的實驗結(jié)果,可以看出,常規(guī)“條帶狀”航線的總飛行距離為3 447.6 m,本文方法總飛行距離為1 090.6 m,相比“條帶狀”航線縮減了2 357 m,降幅在68.4%左右。從具體指標(biāo)來看,針對本文方法的航線效率為91.2%,“條帶狀”航線為38.5%。整體的作業(yè)效率方面,本文方法作業(yè)效率為71%,傳統(tǒng)方法為20.3%。實驗結(jié)果表明,在建筑物測圖類型的任務(wù)中,本文方法可以明顯提高無人機的作業(yè)效率,降低工作時間,同時保證航片質(zhì)量。

表1 兩種航線規(guī)劃算法對比結(jié)果

4 結(jié)束語

本文提出了一種針對建筑物測圖的無人機全自動航線規(guī)劃方法,首先利用NDVI結(jié)果分離區(qū)域內(nèi)的植被;再利用多視立體幾何方法粗略估計地表地勢,計算出建筑物區(qū)域的位置;并提出了一種增強8方向的邊緣檢測算子來估計建筑物的主方向,調(diào)整航向角度;最后利用Hopfield網(wǎng)絡(luò)迭代求解全局的最佳飛行順序,得到完整航線。實驗表明,對于建筑物航測,該方法比“條帶狀”航線在飛行距離上減少了68.4%,總作業(yè)效率提高了50.7%。