基于衛星圖像的平原植保作業區域檢測研究

邵長超,王思捷

（河南科技學院信息工程學院,河南 新鄉 453003）

植保無人機作業主要受非植保作業時間、能量消耗、路程、植保作業的覆蓋率和遺漏率等約束條件影響[1].目前植保無人機航線主要靠人工規劃,在進行植保作業區域,尤其是遇到不規則作業區域時容易產生較大誤差,致使錯誤估算施藥量,增加環境污染和植保作業時間.因此在植保作業前需要對作業區域感知,獲取作業區域的形狀和面積等信息,為作業規劃提供二維數學計算環境[2].

目前,已有國內外學者圍繞相關問題進行了研究.針對作業區域邊緣問題,Umaselvi 等[3]基于聚類的彩色圖像無監督分割方法,對衛星圖像中的植被和市區進行分類.張明杰等[4]利用大飛機平臺獲取的高分辨率航空影像對平坦地區較規整的田埂進行研究,以提取田埂界線信息.針對區域面積計算問題,楊靜靜等[5]采用經典數學定理鞋帶公式（Shoelace Formula）及計算機圖形學中向量積法計算區域面積.孫藝哲等[6]針對小塊農田和不規則農田使用基于改進后的Alpha Shapes 算法農機作業面積測量方法,結果誤差控制在3.5%以內.

可以看到,目前研究大都建立在高分辨率的遙感影像基礎之上,民用獲取尚有難度,使用分辨率有限的民用衛星圖像,可有效地降低作業規劃的前期準備難度與成本.本文基于分辨率有限的民用衛星圖像,圍繞作業區域邊緣檢測和面積計算問題進行研究,為作業規劃和導航提供研究基礎.

1 研究方法設計

總體設計是獲取到衛星地圖后對圖像進行處理,獲得作業區域的二值圖像,最后運用像素比例計算實際作業面積,組織結構見圖1.研究具體過程為：首先利用高德地圖[7]獲取衛星圖像,對獲取到的衛星圖像進行預處理;然后使用邊緣檢測算子檢測作業區域邊界,獲得邊緣檢測圖像后疊加掩膜與運算得到帶有綠色域的圖像,之后將綠色色域轉換為背景色;再后對處理圖像進行形態學處理與端點檢測修補,連接斷點以獲得完整連通區域,對連通區域進行輪廓繪制和填充;最后將二值圖像進行像素比例計算,依據顏色占比與實際圖像總面積換算可獲得作業區域的實際面積.

圖1 組織結構Fig.1 The structure of organization

為模擬真實植保無人機作業規劃軟件,運用高德地圖開放平臺JS API[8]開發了多邊形區域面積簡易測量工具,用以與圖像處理得到的面積進行擬合度對比.區域面積測量工具采用Shoelace 算法（也稱高斯面積公式）,公式以n×2 的矩陣形式表示多邊形上順序排列的頂點,行列式的計算又存在錯位,形如所系的“鞋帶”而得名[9].其缺點是無法直接應用在區域內有弧度的節點上,即使方法改造后也會喪失一部分精度,加之為手工固定測量點,測量精度無法保證[10].

1.1 邊緣檢測

本文使用Python 環境下Open C V 庫對圖像進行處理.為了得到最佳的結果,獲取圖像后使用了多種邊緣檢測算子進行對比檢測.主要有一階Roberts Cross 算子、Sobel 算子等,二階Laplacian 算子、Canny 算子等.Sobel 算子在邊緣定位方面有較好的表現,缺點是容易產生過多的像素邊緣[11];Roberts 算子雖然在定位邊緣的準確度上也相對較好,但在抗噪聲方面表現稍有遜色.經過驗證對比最終選擇了具有較高的檢測準確率以及信噪比,且圖像邊緣的連續性和完整性較好的Canny 算子[12].Canny 算子的檢測過程為[13-14]：

過程一：將彩色圖像轉換為灰度圖像;

過程二：對圖像進行高斯濾波去除圖像噪聲;

過程三：使用一階偏導的有限差分計算梯度的幅值和方向;

過程四：對梯度幅值進行非極大值抑制,尋找像素點局部最大值;

過程五：雙閾值算法檢測和連接邊緣.

通過對Canny 算子分析,Canny 算子圖像最終得到的為灰度圖像,無法區分檢測到的邊緣為田埂或者綠色作物,容易影響到邊緣的精度.針對該問題,本文在Canny 算子的基礎上對其做出了適應性改進,主要思想是將原始圖像和Canny 結果疊加掩膜后進行and 運算,得到帶有綠色域的圖像,之后將圖像模式轉換為HSV 顏色模型,并將綠色域轉換為黑色背景色,最后得到無綠色域作物的圖像,見圖2.其具體的實現過程為：

圖2 原始圖像與Canny 圖像Fig.2 Original image and Canny image

過程一：將原始圖像進行灰度處理;

過程二：為了減少噪聲對邊緣檢測結果的影響,采用自適應濾波核對圖像去除噪聲;

過程三：使用一階偏導的有限差分計算梯度的幅值和方向;

過程四：對梯度幅值進行非極大值抑制,尋找像素點局部最大值;

過程五：雙閾值檢測,閾值設定為最小,保留最多細節;

過程六：Canny 灰度圖像與原始圖像利用掩膜（mask）進行二進制and 運算;

過程七：圖像轉為HSV 格式,將圖像綠色域顏色轉換為背景色;

1.2 斷點修補

得到色域轉換后的灰度圖像,發現存在邊緣像素斷點問題.為獲得圖像封閉區域,需要對圖像進行斷點檢測和斷點缺陷修補.為減少算法開銷,本文先使用形態學閉運算進行斷點修補,閉運算過程為圖像先膨脹后腐蝕,其作用是填充物體內細小空洞,連接鄰近物體和平滑邊界[15].經過驗證,本文方案使用矩形膨脹方法,結構元2×2,迭代次數2 時可以將大部分斷點封閉.算法閉運算公式

形態學處理后大部分的斷點已經封閉,但也暴露出其缺點,如圖3 所示.

圖3 閉運算連通Fig.3 Closed operation

當兩條輪廓像素距離相距特別近時,在進行閉運算操作時容易使兩個輪廓黏連在一起,對輪廓的提取會產生一定影響.所以要使斷點連接只在輪廓端點進行,需要對其進行端點檢測,結果見圖4.

圖4 端點檢測Fig.4 Endpoint detection

由圖4 可知,以P1 為中心的八鄰域,P1 周邊數值的和為1 時為端點,和為0 時為孤點.得到端點后,斷點修補的方法將非常明確,只需要判斷兩端點的距離設置一定的閾值,小于等于該閾值即可將斷點用直線連接.

1.3 輪廓繪制填充

植保作業區域斷點修補后,形成完整連通區域,再通過輪廓提取獲得可參考的植保作業區域,而進行輪廓填充的主要目的為統計二值圖像顏色像素比例.以圖5 為例,主要過程為首先使用cv.findContours（）函數提取全部圖像輪廓;再將提取到的全部輪廓使用cv.contourArea（）函數進行面積計算,通過計算只保留并繪制出的最大面積輪廓,也就是植保作業區域;最后將最大輪廓進行填充并統計二值圖像顏色像素顏色個數與比例.

圖5 輪廓提取與填充Fig.5 Contour Extraction and Contour Filling

1.4 作業區域面積

通過對面積計算函數分析,使用opencv 提供的面積計算函數contourArea（）時得到的值并非輪廓內亮點的像素值.根據官方文檔提供的信息可知,contourArea（）使用格林公式計算面積,其復雜度比求非0像素數低,在像素分辨率有限且相對較低時,兩種算法會出現不同結果：對于大或正方形的輪廓時,誤差最小;對于較小輪廓或者橢圓形等不規則輪廓時,誤差會較大.像素輪廓見圖6.

圖6 5×4 像素輪廓Fig.65×4 Image boundary

驗證圖6 中的5×4 像素輪廓,contourArea（）算法取輪廓像素中心點計算像素數,其面積ContourArea=12,計算輪廓內非0 像素數面積NonZero=20.結果可知,兩種算法會產生不同的結果且誤差較大.為減少兩者誤差,本文使用求非0 像素數獲得作業區域面積,主要思想是二值圖像黑白兩色像素占比與實際總面積的比例換算,求得植保作業區域實際面積.具體算法過程如下

過程一 求出圖像的實際長寬,公式為

式（2）、（3）中：L 表示實際長度,H 表示實際高度,l 表示寬像素數、h 表示高像素數,d 表示DPI,2.54表示英寸與厘米的比例關系,b 表示地圖圖像比例尺.

過程二 計算植保作業區域實際面積,公式為

式（4）中：S 表示作業區域面積,a 表示輪廓非0 像素數量,lh 表示圖像總像素數,LH 表示土地面積.

2 驗證與分析

獲得到輪廓填充圖像后最重要的一步是計算圖像面積.因為衛星地圖圖像面積精確程度直接影響到最后的植保作業區域面積精度,所以首先需要求得衛星圖像長寬與真實土地的比例關系.高德地圖提供了一種地圖縮放級別（zoom）,在高德開放平臺的開發文檔中,地圖顯示縮放級別范圍默認為3～18,取值范圍3～18[16].因開發文檔中沒有給出縮放等級與實際距離的比,不過通過驗證計算可以得知,圖像分辨率在96 DPI 環境下,每減小一個縮放等級距離像素比增大一倍,計算公式如下

式（5）中：b 表示比例尺,z 表示縮放等級.

根據計算公式可得比例尺與地圖中的縮放級別的映射關系表1.

表1 縮放等級與映射關系Tab.1 Zoom level and mapping relationship

2.1 驗證結果

地圖驗證采取高德地圖縮放等級zoom=18（驗證結果表明縮放等級越小,誤差越大）下的衛星圖像作為實驗數據,地圖采集地點在河南省新鄉市紅旗區某村農田,采集的模擬作業區域數為4 塊.地圖圖像采集完成后,使用本文建構的方案對作業區域輪廓和面積計算方法加以驗證,結果見圖7 至圖10.

圖7 區域1 驗證結果Fig.7 Verification results of region 1

圖8 區域2 驗證結果Fig.8 Verification results of region 2

圖9 區域3 驗證結果Fig.9 Verification results of region 3

圖10 區域4 驗證結果Fig.10 Verification results of region 4

以上驗證的4 塊區域圖像中a 表示原始衛星圖像,b 表示模擬作業區域輪廓提取,b 表示經過填充的區域,d 表示程序測量的結果.測量面積對比見表2.

表2 測量面積對比Tab.2 Results comparison

驗證結果表明,本文方案對衛星圖像處理后獲得的輪廓信息,一定程度上抑制了邊緣檢測中的偽邊緣,有效地清除掉了農田中的非種植區域,顯示出農田真實輪廓.通過對表2 計算面積與程序測算面積的比對驗證,兩種方法在以上驗證地塊土地中的擬合程度最高達到了99.53%,從均值來看兩種方式擬合度均值達到了95.78%.由本文方案得到的農田輪廓與面積信息表明具有較高的精確度,可以為植保作業提供數據參考.

2.2 誤差分析

從驗證結果上看,本文方案所得到的結果與手工測算具有些微誤差,在此對誤差產生原因進行分析：

（1）由于面積計算程序采用的算法產生的誤差問題.由于鞋帶公式的原理決定了無法避免測量誤差,加之區域節點為手工添加也是圖像擬合度降低的一個原因.

（2）衛星遙感影像微分糾正誤差問題.衛星遙感影像偏斜是指由于地球自轉,造成掃描行西移,使衛星影像呈平行四邊形的偏斜現象[17].數字微分糾正指根據有關的參數與數字地面模型,利用相應的構像方程式,或按一定的數學模型用控制點解算,從原始非正射投影的數字影像獲取正射影像[18].比較有趣的新聞就是“三峽歪了”[19],該新聞描述的內容就是微分糾正傾誤差時出現了問題.本文方案利用高德獲取到的地圖都是經過微分糾正算法糾正過的圖像.如果要進一步減少誤差,可使用農田檢測無人機提前獲取作業區域的正視投影,但這無疑增加了部署時間,浪費人力物力.

（3）獲取到的衛星圖像分辨率低與比例尺誤差的影響.由于獲取到的衛星分辨率有限,圖像內部分細節丟失,加上民用地圖的比例尺精度要求有限等因素,所以無法避免誤差.

（4）衛星地圖投影產生的失真問題.高德地圖使用的Web Mercator 投影,也是誤差產生的重要原因.

3 小結

本文主要將提取的衛星地圖進行圖像處理,采用改進的邊緣檢測算子檢測了植保作業區域邊緣,對圖像進行形態學處理與端點檢測修補后得到封閉作業區域,最后基于像素比例計算得到作業區域的面積數據.經過對本文方法驗證,使用分辨率有限的民用衛星圖像,其結果也可較好地得到植保作業區域的輪廓范圍,能一定程度的減少植保作業前期準備工作.再通過對本文方法與面積計算工具的對比,結果得到兩者具有較高的擬合度,可以為植保無人機的作業規劃提供較為精確的二維數學計算環境.