魏新華 張 敏 劉青山 李 林

(江蘇大學江蘇省農業(yè)裝備與智能化高技術研究重點實驗室， 鎮(zhèn)江 212013)

0 引言

農業(yè)機械無人駕駛可以減少勞動力需求，提高作業(yè)精度和效率[1-2]，基于聯(lián)合收獲機的無人駕駛技術的研究尤為重要。傳統(tǒng)的聯(lián)合收獲機無人駕駛技術主要是基于衛(wèi)星導航信號實現車輛自主行走[3-4]。在收獲機作業(yè)時，需要實時監(jiān)測作物高度，為撥禾輪高度實時調控提供依據。當無人收獲機處于田地邊緣作物邊界參差不齊、樹木遮擋無法接收衛(wèi)星信號以及衛(wèi)星信號短暫丟失等場景時，僅依靠衛(wèi)星導航不能有效作業(yè)，需要對收割邊界進行實時監(jiān)測，實現輔助導航。

目前關于作物高度檢測的研究方法主要有視覺和傳統(tǒng)測距兩種。SRITARAPIPAT等[5]通過圖像處理分析稻田中固定標記棒和水稻的相對高度，以獲取實際水稻高度。文獻[6-8]通過深度相機進行作物高度信息提取。文獻[9-10]使用雙目立體相機的視差圖實現作物高度計算。此外，文獻[11-12]利用激光雷達測距系統(tǒng)進行作物高度測量。以上測高方式或固定測量、或需要同時檢測作物的底點與頂點、或僅適用于有壟種植方式，均不適用于稻麥作物的實際收獲場景。

邊界檢測方法主要有視覺和激光兩種。文獻[13-16]通過圖像處理的方法進行邊界提取，文獻[17-18]利用激光測距方法，根據收割邊界兩側高度差異來檢測邊界。以上視覺方法均基于單目相機，在谷物成熟度不同、背景對比度不同等多樣化的收獲場景下，算法缺乏穩(wěn)定適用性。激光雷達不適用于充滿灰塵和秸稈碎末的收獲場景，且在作物生長稀疏的場景無法有效提取收割邊界線。

為滿足無人駕駛收獲機在動態(tài)收割場景下對作物高度和收割邊界的實時檢測，本文提出一種基于雙目視覺的稻麥聯(lián)合收獲機田間作物高度和收割邊界的提取方法。利用雙目相機獲取被測物體的三維信息，將其轉換為實際高度，結合彩色圖像對作物高度進行計算，同時根據實際高度數據檢測收割邊界，實現對田間作物高度和收割邊界的動態(tài)提取。

1 高度數據獲取

1.1 數據采集

雙目視覺技術是將2個普通相機相距固定的間隔安裝，先進行相機標定獲得內外參數，校正畸變圖像，然后對左右相機拍攝的2幅圖像進行立體匹配，確定被測物體分別在2幅圖像上的位置，最后根據被測物體在2幅圖像上的位置差異和幾何關系，來計算其深度信息。本文使用ZED雙目相機，該相機可輸出左右2幅彩色圖像、深度圖、相機坐標系下三維坐標等。

數據采集系統(tǒng)包括計算機端、相機和IMU(慣性測量單元)，相機安裝在中聯(lián)重科PL60 4LZT-6.0ZE型稻麥聯(lián)合收獲機上，使用Inter i5處理器、運行Visual Studio 2015軟件開發(fā)工具，相機安裝位置和數據采集系統(tǒng)如圖1a所示。對水稻進行收獲作業(yè)，過程中獲取左右2個相機分別拍攝的彩色圖像、左圖像素點對應的左相機坐標系下的三維相機坐標，彩色圖像分辨率為1 280像素×720像素，圖1b、1c為一組數據示意圖。

圖1 數據采集系統(tǒng)和數據示意圖Fig.1 Data acquisition system and data diagram1.計算機端 2.相機+IMU

1.2 實際高度獲取

作物高度計算以及根據高度數據檢測收割邊界，需對三維相機坐標進行處理，獲得其對應的實際高度。圖2為作業(yè)場景示意圖，(O,X,Y,Z)為左相機坐標系，(x,y,z)為相機坐標，h為點(x,y,z)實際高度。計算地面平面方程后，根據點到平面距離公式計算實際高度。

圖2 作業(yè)場景示意圖Fig.2 Schematic of operation scene1.相機 2.作物 3.地面所在平面

1.2.1初始平面標定

為計算作業(yè)前的初始平面，在地面平坦且沒有作物的區(qū)域采集數據，獲取該場景三維相機坐標。對三維數據進行濾波處理，先進行統(tǒng)計濾波，去除明顯離群點，以提高平面計算的準確性；然后通過直通濾波選取地面感興趣區(qū)域，選取坐標值屬于中間部分的數據，去掉處于邊緣及距離較遠等精度相對較低的數據，同時去除收獲機割臺部分及周圍無關環(huán)境。采用隨機抽樣一致算法(Random sample consensus， RANSAC)對濾波后的數據進行平面擬合[19]，擬合結果視為初始平面，如圖3所示。

圖3 初始平面標定示意圖Fig.3 Initial plane calibration diagram

1.2.2實時平面獲取

由于田間地面起伏不平并非理想平面，收獲機作業(yè)時車身會發(fā)生俯仰、橫滾以及地盤凹陷等情況，導致相機產生位姿變化，如果按照初始平面方程作為整個作業(yè)過程的基準平面，會對后續(xù)作物信息的計算造成誤差。因此，使用WT901C型IMU獲取相機的位姿變化，對地面平面進行實時校正，減少計算物體實際高度數據時產生的誤差。

安裝相機和IMU時，使兩者各軸方向一致，根據IMU狀態(tài)獲得相機坐標系的變化。相機和IMU坐標軸方向如圖4a所示，具體安裝方式如圖4b所示，相機靜置在一平面板上，相機底面與平面板貼合，調整IMU位置，使其頂面與相機背面貼合，同時側邊與平面板貼合，圖4c為安裝方式右視圖，此時相機與IMU各軸方向一致，將相機和IMU按照該狀態(tài)固定。

圖4 相機和IMU安裝示意圖Fig.4 Installation diagram of camera and IMU

IMU數據基于預積分與相機數據融合，第k、k+1幀圖像數據分別對應第i、j時刻的IMU數據，對第i、j時刻之間測量的IMU數據進行預積分，計算旋轉、速度和位移增量ΔRij、Δvij和Δpij。

(1)

式中 ΔRij——i時刻到j時刻的旋轉增量

Δvij——i時刻到j時刻的速度增量

Δpij——i時刻到j時刻的位移增量

ΔRin——i時刻到n時刻的旋轉增量

Δvin——i時刻到n時刻的速度增量

nαn、ngn——IMU噪聲

bαn、bgn——IMU零偏

Δt——IMU測量值時間間隔

結合預積分增量和前一幀圖像時刻對應的IMU狀態(tài)，更新IMU狀態(tài)

(2)

式中Rk、vk、pk——k幀圖像對應的IMU狀態(tài)

Rk+1、vk+1、pk+1——k+1幀圖像對應的IMU狀態(tài)

gw——世界坐標系下重力加速度

Δtij——圖像幀時間間隔

選取初始平面上3點，計算作業(yè)過程中IMU位姿變化后對應的新坐標，根據新坐標計算此時作業(yè)的實時地面平面。相機坐標值變換前后關系為

(3)

式中 (x,y,z)——初始狀態(tài)相機坐標

(x′,y′,z′)——當前狀態(tài)相機坐標

R′——兩狀態(tài)之間的旋轉矩陣

p′——兩狀態(tài)之間的平移矩陣

O——零矩陣

R′和p′分別由IMU此刻旋轉狀態(tài)和平移狀態(tài)獲得。

1.2.3實際高度計算

被測物體的實際高度即為相機坐標(x,y,z)到地面平面的距離，根據點到平面距離計算公式有

(4)

式中Ac、Bc、Cc、Dc——相機坐標系下地面平面方程參數

其中，Acx+Bcy+Ccz+Dc=0為相機坐標系下地面平面方程。

驗證基于相機測高的精度，在3種不同場景下，對與相機相隔不同距離的物體進行高度測量，每種場景測10組數據，與被測物體的實際高度比較，結果如表1所示，距離越遠，誤差越大，平均絕對誤差為0.006 8 m。感興趣區(qū)域內作物與相機距離小于5 m，本文測高方法滿足實際測高需求。由于地面和相機呈一定角度，相機坐標轉換到高度數據時，相機深度測量帶來的誤差會一定程度上減小。

表1 相機測高精度Tab.1 Camera height measurement accuracy

將相機坐標轉換為實際高度數據后，獲取和左彩色圖像像素點一一對應的高度灰度圖，像素值代表實際高度，單位為m。選取最近待收割區(qū)域中干擾最少的區(qū)域作為感興趣區(qū)域，在此區(qū)域內進行作物信息計算，如圖5a紅色矩形框所示區(qū)域，圖5b為感興趣區(qū)域對應的高度灰度圖，圖5c為感興趣區(qū)域高度三維示意圖。

圖5 獲取感興趣區(qū)域Fig.5 Obtaining region of interest

2 作物高度和收割邊界提取

2.1 作物高度提取

本文將感興趣區(qū)域分為：水稻上部、水稻中下部和已收割區(qū)域。水稻上部為作物高度計算的目標區(qū)域，水稻中下部和已收割區(qū)域統(tǒng)稱為背景。水稻上部高度大于背景高度，但是水稻上部和中下部之間存在高度的連續(xù)變化，兩者邊界處高度差值較小，因此僅依靠高度數據無法將其準確區(qū)分。彩色圖像水稻上部與水稻中下部具有明顯邊界特征，但水稻收獲場景多樣化、作物成熟度不同等，難以找出統(tǒng)一的特征區(qū)分水稻上部和已收割區(qū)域。因此，本文提出一種融合高度數據和彩色圖像的方法計算水稻高度，圖6為水稻高度計算流程圖。

圖6 水稻高度計算流程圖Fig.6 Flow chart of rice height calculation

首先采取聚類方法對高度數據進行劃分，將高度平均值最大的類作為水稻高度計算的候選區(qū)域；同時對彩色圖像進行處理，將水稻上部目標和中下部區(qū)分出來；然后融合兩者結果，得到水稻高度計算目標區(qū)域，計算該區(qū)域內高度數據均值，獲取水稻高度。

2.1.1高度數據聚類

K均值算法應用廣泛，計算速度快，但要隨機選取初始聚類中心，聚類個數需要預先指定，因此聚類結果受初始值影響較大且不穩(wěn)定。密度峰聚類算法(Density peaks clustering，DPC)是由RODRIGUEZ等[21]提出的算法，主要思想是結合密度及距離2種特征，選取聚類中心進行聚類，但密度峰聚類算法需要計算每個數據點和其余所有點的距離，不適合數據樣本大的場合，且需要人工根據決策圖選取聚類中心，不適用實時處理的場景。針對以上2種算法優(yōu)缺點，本文提出一種基于改進的結合密度峰聚類和K均值聚類的方法對高度數據聚類。

改進后的聚類算法主要流程為：作出高度灰度圖的直方圖對應的曲線，將曲線極大值點作為初始聚類中心候選點，引入權值綜合描述初始聚類中心候選點的密度和距離，規(guī)定權值判定準則自適應選取初始聚類中心，將高度數據進行區(qū)間劃分犧牲了一定精確性，改進算法獲得的初始聚類中心與密度峰聚類算法結果會有較小差值，最后用K均值聚類迭代計算達到最佳聚類結果。

改進算法具體步驟如下：

(1)將高度數據劃分到m個區(qū)間，作出直方圖，計算直方圖對應的曲線f(x)的n個極大值點xi={x1,x2,…,xn}，圖7a為極大值示意圖，將極大值點當作高局部密度候選點，即候選聚類中心，對應的極大值為密度ρi，即

圖7 聚類過程參數Fig.7 Clustering process parameter diagram

ρi=f(xi) (i=1,2,…,n)

(5)

(2)計算極大值點xi的最小距離σi，遍歷所有極大值點，找到比極大值點xi密度大且與xi距離最小的極大值點xj，由于本文數據范圍較小，為減小密度過大的影響，xj與xi之間的距離乘以區(qū)間個數m進行修正，σi為

σi=min(|xi-xj|m)
(ρj>ρi;i,j=1,2,…,n)

(6)

得到圖7b所示的決策圖。

(3)計算每個極值點的簇中心權值λi，定義為該點密度和距離的乘積，即

λi=ρiσi(i=1,2,…,n)

(7)

將所有簇中心權值降序排列，計算該點與后一點的差值，拐點定義為首次出現且差值小于閾值的點,拐點之前的簇中心權值對應的極大值點為初始聚類中心。圖7b較大的3個圓點為自適應選取的聚類中心。

(4)確定初始聚類中心和聚類數目后，計算其余高度數據到各個聚類中心的距離,并歸為最小距離所對應的類。

(5)其余的點遍歷完成并歸類后，對每一類計算均值,作為該類的新聚類中心,并重復步驟(4)，直至聚類中心變化小于設定閾值。

根據以上步驟對圖5b的高度數據聚類，得到圖8a所示的聚類結果，由圖8a看出，該樣本被分為3類：①包含大部分水稻上部的數據。②包含水稻中下部和小部分水稻上部的數據。③已收割區(qū)域數據。且水稻上部類均值最大，中下部次之，未收割最小。將均值最大的類作為計算水稻高度的候選區(qū)域，即包含大部分水稻上部數據的類，如圖8b所示。

圖8 聚類結果Fig.8 Clustered result graph

基于本文改進算法和基于密度峰聚類算法選取的初始聚類中心值平均差值為0.010 3 m，占數據總范圍的0.96%，表明本文改進算法計算的初始聚類中心具有較好密度表達能力。在當前數據集下，密度峰聚類算法僅計算所有點兩兩之間距離，耗時大于5 min，且不包括密度計算等操作，改進之后得到聚類中心候選點的平均時間0.131 0 s，整個聚類過程平均時間0.701 4 s。K均值迭代次數，在隨機產生聚類中心時迭代次數平均為18.4次，改進后平均為2.5次。改進后的聚類算法不僅有效自適應找到較為準確的初始聚類中心，且大大縮短了運行時間。

2.1.2基于彩色圖像的水稻目標提取

水稻上部和中下部由于受到光照強度不同，亮度會有明顯差距，而且兩者顏色也有明顯差異，水稻上部顏色偏黃，中下部偏綠。對同時包含亮度和顏色信息的R、G、B分量進行歸一化，由RGB顏色空間轉換到rgb顏色空間，歸一化后的顏色特征因子r、g、b分別為

(8)

圖9a為圖5a中感興趣區(qū)域歸一化的結果。從圖9a可以看到，未收割區(qū)域中歸一化圖像的水稻上部目標區(qū)域比原圖更顯著。對歸一化后的圖像進行灰度化，利用最大類間方差法自適應獲取閾值進行圖像分割，得到二值化結果，如圖9b所示。圖中白色區(qū)域為目標區(qū)域，包括水稻上部區(qū)域，以及已收割區(qū)域中的部分區(qū)域。

圖9 彩色圖像目標提取結果Fig.9 Diagram of color image target extraction results

2.1.3作物高度計算

將高度聚類結果中均值最大的類作為水稻高度計算候選區(qū)域，把該區(qū)域反映到圖像對應位置。候選區(qū)域和彩色圖像處理結果中的目標區(qū)域進行與運算,得到共同目標區(qū)域，即圖10白色區(qū)域。該區(qū)域內的高度數據集為Rh，對屬于Rh的高度數據h計算均值，獲取水稻高度

圖10 水稻高度計算區(qū)域示意圖Fig.10 Diagram of calculated rice height area

(9)

式中s——Rh集合內數據個數

2.2 收割邊界提取

互相關函數在信號分析中表示2個信號之間的相關程度，互相關系數越大，表明相關程度越高，描述2個不同信號的相關性時，既適用于隨機信號，也適用于確知信號。圖11曲線為高度灰度圖中一行的高度數據，已收割區(qū)域和未收割區(qū)域存在著明顯高度差異，呈現變形階躍函數的特征?；谶@一數據分布特點，建立模型函數，根據灰度圖中行高度數據序列和模型函數之間的互相關性提取收割邊界點。

圖11 圖像和高度對照Fig.11 RGB image and height comparison chart

模型函數q(t)定義為

(10)

利用模型函數從左至右檢測行高度序列得到收割邊界點，計算每一時刻高度數據序列與模型函數的互相關系數，互相關系數最大值對應的時刻dmax為收割邊界點對應的位置。在d時刻兩者互相關系數定義為[22-23]

(11)

式中rhq(d)——模型函數d時刻的互相關系數

h(t)——高度數據序列

mh——高度數據序列平均值

mq——模型函數平均值

N——高度數據序列的數據個數

為保證準確度的同時提高計算速度，對高度灰度圖每10行進行一次邊界點提取。利用最小二乘法對邊界點進行直線擬合，得到收割邊界直線。邊界點和擬合直線結果在高度灰度圖中位置如圖12a所示，結果顯示邊界點可被準確識別,最下面的邊界點附近作物輕微向內傾斜，使提取的邊界點偏左，利用最小二乘法進行擬合時減小了該點影響。圖12b為彩色圖中邊界點和邊界線對應位置，圖12c為三維高度圖中邊界位置示意圖，為方便顯示，邊界線高度取為0.8 m。

圖12 收割邊界提取結果Fig.12 Cut-edge extraction result

2個相鄰場景收割邊界線的斜率和邊界點位置變化不會太大，所以可根據當前邊界線對下一幀數據邊界點候選區(qū)域進行預測。收獲機正常作業(yè)時，位置連續(xù)調整的左右范圍在15 cm以內，由相機坐標獲知15 cm對應感興趣區(qū)域底部圖像行的平均寬度為22像素，因此規(guī)定50像素為邊界點候選區(qū)域寬度。在下一幀邊界檢測時，以當前邊界線上的點為中心，邊界點候選區(qū)域為50像素，在候選區(qū)域內從左至右進行邊界點檢測。

計算收割邊界的目的主要是提供導航依據，因此要根據擬合的邊界直線計算航向偏差和橫向偏差。航向偏差為邊界直線和收獲機行駛方向所在直線的夾角。相機安裝于收獲機正前方，其距離地面高度3.00 m，收獲機割臺寬2.10 m，相機支架位于中間偏左0.45 m位置。因此相機坐標系的Z軸正軸方向指向收獲機前進方向，如圖13a所示，將Z軸所在直線L0投影到地面平面得到直線L′0，直線L′0表示收獲機行駛方向所在直線。收割邊界直線方程為L1，L1投影到地面平面得到直線L′1，L′0和L′1同在地面平面，L′0和L′1的夾角為航向偏差。L′1相對于L′0左偏定義為正，右偏定義為負，以此區(qū)分邊界線方向。

橫向偏差定義為割臺右側與收割邊界點之間的距離。相機坐標系的X軸與收獲機前進方向垂直，原點為左相機中心。如圖13b中黑色直線為X=0所在直線，X=0所在直線、收割邊界直線和感興趣區(qū)域底端的交點分別為P0、P1，P0、P1對應X軸坐標值分別為X0、X1，且X0=0。左相機距相機中心0.06 m，根據割臺寬度和相機安裝位置，則左相機中心距離割臺右側邊緣0.66 m，割臺右側X軸坐標值為X2，即X2=0.66 m，則橫向偏差為X2-X1。

圖13 航向偏差和橫向偏差示意圖Fig.13 Diagrams of course deviation and lateral deviation

3 結果分析

3.1 作物高度

在鎮(zhèn)江市大港新區(qū)進行現場實驗，對成熟期的南粳9108品種水稻進行收獲，靜態(tài)場景采集10個樣本，每個樣本用卷尺測量20組水稻高度數據，測量高度為地面到穗最高點或葉片最高點的距離。實際高度為20組高度數據去掉最大值和最小值的平均值，將采集的樣本實際高度和計算高度進行比較，結果見表2。

從表2可以得出，高度計算平均絕對誤差為0.043 m，標準差為0.020 m。實際田中水稻的頂部并不是一個理想的平整面，而是分布在一定高度范圍內，所以計算的水稻高度平均誤差小于5 cm，滿足測量需求。由于作物高度計算時對目標區(qū)域的整穗和整個葉片上部進行處理，而實際測量穗最高點或葉片最高點，因此計算結果普遍比測量的高度低。田間地面凹凸不平，水稻底點的測量高度不同，也會造成誤差。本文作物高度提取算法平均每幀耗時0.823 s，收獲機作業(yè)速度約1 m/s，感興趣區(qū)域包含前方2 m范圍內作物，每秒視頻數據取一幀數據進行處理，可應用于實時場景。

表2 高度結果比較Tab.2 Height results comparison m

3.2 收割邊界

為驗證本文收割邊界提取算法在不同場景下的適用性，采集多種場景的樣本進行分析，包括常規(guī)、作物生長稀疏、邊緣作物向外傾斜、作物部分存在漏播、作物被施藥機車輪壓過和倒伏等。圖14為不同場景下邊界檢測結果, 邊界均被有效識別。未收割作物內部出現多種干擾情況，割茬呈現不同顏色和紋理，收割區(qū)域和未收割區(qū)域對比程度不同，表明該方法適應性較高，其中對倒伏場景的邊界檢測僅適用于靠近邊界的作物并未完全倒伏的情況。

圖14 不同場景收割邊界檢測結果 Fig.14 Cut-edge detection results in different scenes

對本文算法計算的航向偏差和橫向偏差進行驗證。人工標記的邊界線作為參考邊界線，計算本文算法與參考邊界線之間的夾角，該夾角為航向偏差的誤差。卷尺測量收獲機割臺右端到其最近的邊界點的距離，作為實際橫向偏差。表3為航向偏差和橫向偏差結果比較。

表3 偏差結果比較Tab.3 Deviation results comparison

由表3可得，航向偏差平均誤差為1.04°，實際和計算的橫向偏差之間的平均絕對誤差為0.084 m，表明本文邊界識別較為精確。邊界線與感興趣區(qū)域底部交點和實際割臺右側最近的邊界點兩者存在一定偏差，以及人工測量誤差等都會造成橫向偏差計算的誤差。作物邊界若超過割臺右端發(fā)生漏割，此時橫向偏差小于0，反之則大于0。其中檢測到圖14b、14c橫向偏差小于0，判為發(fā)生漏割，圖14e橫向偏差較大，表明收割范圍較小。對整個高度數據序列進行計算，平均每幀耗時0.802 s，預測候選區(qū)域后，平均每幀運行時間為0.124 s，提高了運行速度，適用于實時的收割場景。對采集的600幀數據進行邊界提取計算，邊界識別正確率為93.30%，表明本文收割邊界提取算法有較高的準確性。

4 結論

(1)利用雙目視覺實現作物高度和收割邊界信息的動態(tài)提取，作物高度提取平均誤差為0.043 m，邊界識別正確率達93.30%，航向偏差平均誤差為1.04°，橫向偏差平均絕對誤差為0.084 m，可以為聯(lián)合收獲機實時自適應調控提供依據。

(2)基于本文改進算法對高度數據分類，可以自適應選取聚類中心，且提高了運行速度。結合彩色圖像提高了作物高度計算的準確性。

(3)計算高度數據序列和模型函數的互相關系數來獲取收割邊界點，根據前一圖像邊界結果界定后一圖像邊界點候選區(qū)域，提高了運行速度。