基于衛(wèi)星圖像的平原植保作業(yè)區(qū)域檢測(cè)研究

邵長超,王思捷

（河南科技學(xué)院信息工程學(xué)院,河南 新鄉(xiāng) 453003）

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

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

可以看到,目前研究大都建立在高分辨率的遙感影像基礎(chǔ)之上,民用獲取尚有難度,使用分辨率有限的民用衛(wèi)星圖像,可有效地降低作業(yè)規(guī)劃的前期準(zhǔn)備難度與成本.本文基于分辨率有限的民用衛(wèi)星圖像,圍繞作業(yè)區(qū)域邊緣檢測(cè)和面積計(jì)算問題進(jìn)行研究,為作業(yè)規(guī)劃和導(dǎo)航提供研究基礎(chǔ).

1 研究方法設(shè)計(jì)

總體設(shè)計(jì)是獲取到衛(wèi)星地圖后對(duì)圖像進(jìn)行處理,獲得作業(yè)區(qū)域的二值圖像,最后運(yùn)用像素比例計(jì)算實(shí)際作業(yè)面積,組織結(jié)構(gòu)見圖1.研究具體過程為：首先利用高德地圖[7]獲取衛(wèi)星圖像,對(duì)獲取到的衛(wèi)星圖像進(jìn)行預(yù)處理;然后使用邊緣檢測(cè)算子檢測(cè)作業(yè)區(qū)域邊界,獲得邊緣檢測(cè)圖像后疊加掩膜與運(yùn)算得到帶有綠色域的圖像,之后將綠色色域轉(zhuǎn)換為背景色;再后對(duì)處理圖像進(jìn)行形態(tài)學(xué)處理與端點(diǎn)檢測(cè)修補(bǔ),連接斷點(diǎn)以獲得完整連通區(qū)域,對(duì)連通區(qū)域進(jìn)行輪廓繪制和填充;最后將二值圖像進(jìn)行像素比例計(jì)算,依據(jù)顏色占比與實(shí)際圖像總面積換算可獲得作業(yè)區(qū)域的實(shí)際面積.

圖1 組織結(jié)構(gòu)Fig.1 The structure of organization

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

1.1 邊緣檢測(cè)

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

過程一：將彩色圖像轉(zhuǎn)換為灰度圖像;

過程二：對(duì)圖像進(jìn)行高斯濾波去除圖像噪聲;

過程三：使用一階偏導(dǎo)的有限差分計(jì)算梯度的幅值和方向;

過程四：對(duì)梯度幅值進(jìn)行非極大值抑制,尋找像素點(diǎn)局部最大值;

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

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

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

過程一：將原始圖像進(jìn)行灰度處理;

過程二：為了減少噪聲對(duì)邊緣檢測(cè)結(jié)果的影響,采用自適應(yīng)濾波核對(duì)圖像去除噪聲;

過程三：使用一階偏導(dǎo)的有限差分計(jì)算梯度的幅值和方向;

過程四：對(duì)梯度幅值進(jìn)行非極大值抑制,尋找像素點(diǎn)局部最大值;

過程五：雙閾值檢測(cè),閾值設(shè)定為最小,保留最多細(xì)節(jié);

過程六：Canny 灰度圖像與原始圖像利用掩膜（mask）進(jìn)行二進(jìn)制and 運(yùn)算;

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

1.2 斷點(diǎn)修補(bǔ)

得到色域轉(zhuǎn)換后的灰度圖像,發(fā)現(xiàn)存在邊緣像素?cái)帱c(diǎn)問題.為獲得圖像封閉區(qū)域,需要對(duì)圖像進(jìn)行斷點(diǎn)檢測(cè)和斷點(diǎn)缺陷修補(bǔ).為減少算法開銷,本文先使用形態(tài)學(xué)閉運(yùn)算進(jìn)行斷點(diǎn)修補(bǔ),閉運(yùn)算過程為圖像先膨脹后腐蝕,其作用是填充物體內(nèi)細(xì)小空洞,連接鄰近物體和平滑邊界[15].經(jīng)過驗(yàn)證,本文方案使用矩形膨脹方法,結(jié)構(gòu)元2×2,迭代次數(shù)2 時(shí)可以將大部分?jǐn)帱c(diǎn)封閉.算法閉運(yùn)算公式

形態(tài)學(xué)處理后大部分的斷點(diǎn)已經(jīng)封閉,但也暴露出其缺點(diǎn),如圖3 所示.

圖3 閉運(yùn)算連通Fig.3 Closed operation

當(dāng)兩條輪廓像素距離相距特別近時(shí),在進(jìn)行閉運(yùn)算操作時(shí)容易使兩個(gè)輪廓黏連在一起,對(duì)輪廓的提取會(huì)產(chǎn)生一定影響.所以要使斷點(diǎn)連接只在輪廓端點(diǎn)進(jìn)行,需要對(duì)其進(jìn)行端點(diǎn)檢測(cè),結(jié)果見圖4.

圖4 端點(diǎn)檢測(cè)Fig.4 Endpoint detection

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

1.3 輪廓繪制填充

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

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

1.4 作業(yè)區(qū)域面積

通過對(duì)面積計(jì)算函數(shù)分析,使用opencv 提供的面積計(jì)算函數(shù)contourArea（）時(shí)得到的值并非輪廓內(nèi)亮點(diǎn)的像素值.根據(jù)官方文檔提供的信息可知,contourArea（）使用格林公式計(jì)算面積,其復(fù)雜度比求非0像素?cái)?shù)低,在像素分辨率有限且相對(duì)較低時(shí),兩種算法會(huì)出現(xiàn)不同結(jié)果：對(duì)于大或正方形的輪廓時(shí),誤差最小;對(duì)于較小輪廓或者橢圓形等不規(guī)則輪廓時(shí),誤差會(huì)較大.像素輪廓見圖6.

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

驗(yàn)證圖6 中的5×4 像素輪廓,contourArea（）算法取輪廓像素中心點(diǎn)計(jì)算像素?cái)?shù),其面積ContourArea=12,計(jì)算輪廓內(nèi)非0 像素?cái)?shù)面積NonZero=20.結(jié)果可知,兩種算法會(huì)產(chǎn)生不同的結(jié)果且誤差較大.為減少兩者誤差,本文使用求非0 像素?cái)?shù)獲得作業(yè)區(qū)域面積,主要思想是二值圖像黑白兩色像素占比與實(shí)際總面積的比例換算,求得植保作業(yè)區(qū)域?qū)嶋H面積.具體算法過程如下

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

式（2）、（3）中：L 表示實(shí)際長度,H 表示實(shí)際高度,l 表示寬像素?cái)?shù)、h 表示高像素?cái)?shù),d 表示DPI,2.54表示英寸與厘米的比例關(guān)系,b 表示地圖圖像比例尺.

過程二 計(jì)算植保作業(yè)區(qū)域?qū)嶋H面積,公式為

式（4）中：S 表示作業(yè)區(qū)域面積,a 表示輪廓非0 像素?cái)?shù)量,lh 表示圖像總像素?cái)?shù),LH 表示土地面積.

2 驗(yàn)證與分析

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

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

根據(jù)計(jì)算公式可得比例尺與地圖中的縮放級(jí)別的映射關(guān)系表1.

表1 縮放等級(jí)與映射關(guān)系Tab.1 Zoom level and mapping relationship

2.1 驗(yàn)證結(jié)果

地圖驗(yàn)證采取高德地圖縮放等級(jí)zoom=18（驗(yàn)證結(jié)果表明縮放等級(jí)越小,誤差越大）下的衛(wèi)星圖像作為實(shí)驗(yàn)數(shù)據(jù),地圖采集地點(diǎn)在河南省新鄉(xiāng)市紅旗區(qū)某村農(nóng)田,采集的模擬作業(yè)區(qū)域數(shù)為4 塊.地圖圖像采集完成后,使用本文建構(gòu)的方案對(duì)作業(yè)區(qū)域輪廓和面積計(jì)算方法加以驗(yàn)證,結(jié)果見圖7 至圖10.

圖7 區(qū)域1 驗(yàn)證結(jié)果Fig.7 Verification results of region 1

圖8 區(qū)域2 驗(yàn)證結(jié)果Fig.8 Verification results of region 2

圖9 區(qū)域3 驗(yàn)證結(jié)果Fig.9 Verification results of region 3

圖10 區(qū)域4 驗(yàn)證結(jié)果Fig.10 Verification results of region 4

以上驗(yàn)證的4 塊區(qū)域圖像中a 表示原始衛(wèi)星圖像,b 表示模擬作業(yè)區(qū)域輪廓提取,b 表示經(jīng)過填充的區(qū)域,d 表示程序測(cè)量的結(jié)果.測(cè)量面積對(duì)比見表2.

表2 測(cè)量面積對(duì)比Tab.2 Results comparison

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

2.2 誤差分析

從驗(yàn)證結(jié)果上看,本文方案所得到的結(jié)果與手工測(cè)算具有些微誤差,在此對(duì)誤差產(chǎn)生原因進(jìn)行分析：

（1）由于面積計(jì)算程序采用的算法產(chǎn)生的誤差問題.由于鞋帶公式的原理決定了無法避免測(cè)量誤差,加之區(qū)域節(jié)點(diǎn)為手工添加也是圖像擬合度降低的一個(gè)原因.

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

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

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

3 小結(jié)

本文主要將提取的衛(wèi)星地圖進(jìn)行圖像處理,采用改進(jìn)的邊緣檢測(cè)算子檢測(cè)了植保作業(yè)區(qū)域邊緣,對(duì)圖像進(jìn)行形態(tài)學(xué)處理與端點(diǎn)檢測(cè)修補(bǔ)后得到封閉作業(yè)區(qū)域,最后基于像素比例計(jì)算得到作業(yè)區(qū)域的面積數(shù)據(jù).經(jīng)過對(duì)本文方法驗(yàn)證,使用分辨率有限的民用衛(wèi)星圖像,其結(jié)果也可較好地得到植保作業(yè)區(qū)域的輪廓范圍,能一定程度的減少植保作業(yè)前期準(zhǔn)備工作.再通過對(duì)本文方法與面積計(jì)算工具的對(duì)比,結(jié)果得到兩者具有較高的擬合度,可以為植保無人機(jī)的作業(yè)規(guī)劃提供較為精確的二維數(shù)學(xué)計(jì)算環(huán)境.