鐘智超 肖雄武 涂建光
(1. 武漢大學 遙感信息工程學院, 湖北 武漢 430079;2. 武漢大學 測繪遙感信息工程國家重點實驗室, 湖北 武漢 430079)
隨著傳感器與控制技術不斷發展,無人機已經成為測繪遙感領域獲取數據的重要工具[1-2]。目前常見的無人機航線規劃算法致力于求解旅行客問題[3](traveling salesman problem,TSP)。常用的求解法有動態規劃[4]、遺傳算法、蟻群算法、概率模型估計[5]等,但這類算法很少考慮具體的飛行任務,例如,建筑物測圖[6]、道路勘探、電路巡檢[7]、農業植保[8]等。另外,無人機在建筑物測圖領域的應用場景越來越多[9]。例如,在城市規劃中,無人機結合傾斜攝影測量技術可以實現城市實景三維模型的構建;在地圖生產中,無人機高精度的航片可以用于正射影像生產、數字線劃圖生產[10]等。
針對常見航線規劃算法對建筑物測圖不具備針對性的現狀,提出一種針對建筑物測圖的無人機航線自動規劃算法。針對建筑物與植被粘連的問題,結合影像的光譜特性分離植被,縮小作業范圍;設計基于無人機多視圖立體幾何的地表地勢的估計方法,分析出飛行區域內粗略的高區和低區;結合植被信息與地勢信息完成對建筑物區域的分割;設計基于邊緣檢測的建筑物主方向識別,優化航線角度;最后利用霍普菲爾德網絡求解覆蓋多個建筑物區域的最優路徑。
建筑物測圖任務中,研究對象是建筑物,不包含其他地物。而常規的無人機規劃線路往往對整個測區進行“條帶狀”的航帶劃分,如圖1所示,這會導致大量的無用航點和無用航片,尤其是在建筑物集群分布而且群體之間的距離較遠的情況。

圖1 “條帶狀”航線和常見“無效航片”示例
植被是建筑物測圖中“無效航片”的重要來源,建筑物集群分布的特性導致存在大量僅被植被覆蓋的區域,可以利用植被特有的光譜特性來區分其與別的地物類型,如式(1)所示。
(1)
式中,(x,y)為像素;NIR(x,y)為近紅外波段值;R(x,y)為紅色波段的灰度值。歸一化植被系數(normalized differential vegetation index,NDVI)可以用于區分常見的綠色植物。綠色植物對近紅外波段的反射率較高,對可見光紅色波段的吸收率較高,因此將NDVI(x,y)大于閾值的像素標記為植被。
除了植被之外,“裸地”也是非建筑物地物中占地面積較大的部分,尤其對比較干旱的測區,“裸地”的影響程度更大。區分建筑物與“裸地”的重要方式是通過地勢高低差來判斷。在無人機作業中,很難直接獲取精確的地勢高程數據,如數字表面模型(digital surface model,DSM)等,因此可依靠無人機多視角測量來估計地表[11]。
1.2.1無人機多視立體幾何
第一步,利用無人機在起飛點附近俯拍整個測區(根據測區大小調整飛行高度),記錄拍照瞬間無人機傳感器的位姿信息,具體包括外參矩陣R1以及根據全球定位系統(global positioning system,GPS)定位獲取的攝影中心初始位置信息,用矩陣t1表示。值得注意的是,初始時姿態角可以直接利用無人機萬向角信息進行初始化。
第二步,移動無人機并調整鏡頭位姿,在另一位置拍攝,同樣覆蓋整體測區。其間通過無人機慣導數據結合初始化的位姿推算第二次拍攝瞬間的位姿R2和t2。
第三步,根據灰度特征匹配兩張影像,計算出匹配點對的深度,從而獲取對應3D點的估計坐標。
圖2中,地表上不同地物同時被兩次拍攝記錄到,對于地表的一點P(X,Y,Z),其對應兩張影像上的兩個像素點p1(u1,v1)和p2(u2,v2),在已知兩次攝影無人機位姿的估計值,那么有式(2)的關系。
(2)
式中,K為相機的內參矩陣,由于未知量是地表點P(X,Y,Z)坐標,因此至少需要兩組同名點對才能求解該問題,因此需要無人機在不同的位姿對同一地表進行拍攝。

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

圖3 用于估計地形表面的權重函數
圖3中,p(r)為權重關于殘差變化的函數;r為殘差大小;g和t為分段權重的殘差閾值。首先用等權重來加權這些點的坐標Z值,擬合一個初始的地表面,然后逐點比較坐標Z值與擬合表面的估計值的差異,根據差異的大小來重新加權Z值。該分段權重函數把殘差大于閾值t的點直接設置為“外點”,權重設置為0。
把分段擬合后的地表面記為GS,那么區分某點P(X,Y,Z)是否屬于建筑物的方法如式(3)所示。其中(x,y)表示與點P(X,Y,Z)對應的像素坐標。
(3)
邊緣特征描述了像素之間的梯度關系。數學上,用一階導數來表達像素梯度被稱為一階微分算子,如Sobel算子、Canny算子等。利用二階導數來檢測梯度的算子被稱為二階微分算子,如LoG算子等。為了提高建筑物邊界檢測結果的完整性,防止出現對不明顯的邊界的漏檢,設計了一種8方向的增強一階微分算子,算子的結構如圖4所示。該算子對比普通的一階算子的改進包括,第一,從普通一階算子只檢測2個方向的梯度變化,拓展到檢測8個方向;第二,擴大了檢測算子的鄰域大小,從3×3拓展到5×5,增大梯度檢測的感知范圍。

圖4 8方向增強邊緣檢測算子
為了驗證該改進算子的提取能力,設計了與經典的Canny算子、LoG算子的對比實驗,驗證結果如圖5所示。可以看出,虛線框中目標建筑物被陰影遮擋的邊緣梯度很小,提取難度較大,該8方向增強算子能保持較好的檢出率。

圖5 8方向增強算子提取建筑物邊緣的對比結果
建筑物的主方向可以表示該建筑物的朝向,假設用方向向量d表示主方向,首先利用Hough變換對2.1中獲取的建筑物邊緣進行線段化,得到所有表示建筑物邊界的直線段集合L。根據直線段判斷建筑物區域的主方向步驟如下,如圖6所示。

圖6 確定主方向與單區域航線規劃
第1步,對于任意L中的每一條直線段,先計算其長度,對于長度小于閾值(根據經驗一般設定為真實尺寸1.5 m)的設定為噪聲,對于長度大于閾值的保留為有效線段。
第2步,對于每一條有效線段,計算其方向向量。按照每30度為一個區間劃分6個方向角區間,根據方向向量將這些有效線段分到這6個區間內。
第3步,統計每個區間內有效線段的總長度,挑選出總長度最長的一個區間,將區間內的所有線段的方向向量按照長度進行加權求和,得到該區域內建筑物的主方向。
無人機的單次任務是從起飛到降落的全過程,需在盡可能短的時間內,完成對所有建筑物區域的飛行任務,即找到從起飛點到完成所有建筑物拍攝并返回起飛點的最短路徑。Hopfield神經網絡是一種結合存儲系統和二元系統的神經網絡。基于連續Hopfield的特點,如果把一個最優化問題的目標函數轉換成網絡的能量函數,把問題的變量對應于網絡的狀態,那么Hopfield神經網絡就能夠用于解決優化組合問題[15]。
2.2.1能量函數構建
首先構建無人機飛行的Hopfield次序矩陣e,次序矩陣中的元素eij代表無人機在第i時刻執行j建筑物區域的飛行任務。假設所有的目標建筑物區域總數為N,那么將產生如下約束條件。約束1,次序矩陣的每一行只有一個神經元被激活,即每個區域有且只有1次飛行任務;約束2,次序矩陣的每一列只有一個神經元被激活,即無人機不能同時執行2個或2個以上區域的任務;約束3,次序矩陣所有激活的神經元數量之和等于所有的建筑物區域數量總和,即無人機需要執行完所有建筑物區域的飛行任務。三個約束條件如式(4)所示。
(4)
無人機整體測區路徑優化的目標是所有路徑長度加起來最小,如式(5)所示,其中d表示距離矩陣,dik代表第i個飛行區域到第k個飛行區域的距離,默認在初始時對這些待飛區域進行了編號。
(5)
因此,將三個約束條件以及最小化總路徑長度的目標函數轉化為能量函數的懲罰項,可以得到總體能量函數的表達式,如式(6)所示。
(6)
式中,前3項是約束項;第4項為目標函數;A、B、C、D為懲罰參數,為固定值,這些值反映了這幾項在最小化過程中的相對重要性。
然后將式(6)描述的航線優化問題按照Hopfield網絡標準能量函數對應起來,其標準形式如式(7)所示。
(7)
式中,權重和偏重Wijkl、Iij展開表示如式(8)所示。
(8)
式中,δij代表指示函數,其在i=j時取值為1,i≠j時取值為0。
2.2.2迭代求解最優航線
確定航線優化問題的約束條件和標準能量函數之后,需要不斷迭代修正權重矩陣,直到狀態趨于穩定,具體的步驟如下。
第1步,設定懲罰系數A、B、C、D和最大迭代次數Tmax,隨機初始化次序矩陣e。
第2步,計算各個建筑物區域之間的距離,形成距離矩陣d。
第3步,根據初始化的輸入狀態U0,確定狀態轉移函數,如式(9)所示。
(9)
第4步,計算輸入狀態的增量,如式(10)所示。
(10)
第5步,利用動態方程和一階歐拉法更新下一刻的輸入狀態和輸出狀態,如式(11)所示。
(11)
第6步,檢查當前的輸出狀態e是否滿足約束,判斷迭代是否結束,否則返回第4步,直到迭代次數達到設定的最大次數。
為了驗證該航線規劃算法的有效性,利用國際攝影測量與遙感學會(International Society for Photogrammetry and Remote Sensing,ISPRS)提供的公開數據進行模擬實驗,該地區的建筑物占地面為4 099.5 m2。該地區建筑物面積較大,分散分布且樹木環繞,首先對其中的建筑物區域進行定位,結果如圖7所示。

圖7 建筑物區域定位實驗結果
如圖7(a)所示,目標測區內建筑物分散分布且在朝向上具有明顯差異;圖7(b)中亮區代表植被,暗區代表非植被;圖7(c)中亮區代表地勢高的區域,暗區代表低地勢;結合圖7(b)和圖7(c)可以定位具體的建筑物區域,利用增強8方向的檢測算子提取目標建筑物的邊界,如圖7(d)所示。
根據建筑物區域的定位結果與邊緣提取的結果,然后計算每個區域內建筑物的主方向,根據主方向確定每個小區域內的航線。最后將所有小區域的航線全部生成后,通過全局航線優化確定無人機執行各區域飛行任務的順序。在初始化Hopfield距離矩陣時需考慮到建筑物區域本身具有的空間寬度,距離是兩個區域最小外接矩形幾何中心之間的距離。確定初始距離矩陣之后,設置初始系數A=B=D=50,C=20進行Hopfield網絡的迭代趨近,結果如圖8(b)所示,輸出狀態在迭代20次左右逐漸穩定,最小損失穩定在315.654,并得到如圖8(a)所示的整體航線優化結果。

圖8 全局航線優化結果與Hopfield網絡迭代收斂
為了驗證該方法規劃的航線與不考慮建筑物分布的“條帶狀”航線的效率差異,設計了同區域的對比實驗。設置相同的重疊率使得航帶間距為12 m,固定無人機飛行速度為6 m/s,忽略無人機拍照時間與等待,生成的航線結果如圖9所示。為了合理地評價無人機執行飛行任務的效率,這里給出一些指標。
有效航線距離(d1):屬于航線規劃的航線,在建筑物區域(參考圖8(a)所示的建筑物區域定義)內的航線長度。
通勤距離(d2):從當前航線飛往下一航線的中間過程飛行距離。
無效航線距離(d3):屬于航線規劃的航線,但是不屬于建筑物區域內部,在無效航線上拍攝的航片稱為無效航片。
航線效率(r1):有效航線距離與航線規劃的總長度的比值(不包括通勤距離)。
作業效率(r2):有效航線距離與無人機作業全流程飛行總距離的比值(包括通勤距離)。

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

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