Vegetation coverage extraction model of winter wheat based on pixel dichotomy
-
摘要:目的
快速准确提取冬小麦返青期植被覆盖度信息。
方法利用无人机获取田间冬小麦可见光图像，提取图像中4种常见可见光植被指数；在像元二分法原理的基础上，分别构建基于差异植被指数(Visible-band difference vegetation index，VDVI)、过绿指数(Excess green，EXG)、归一化绿蓝差异指数(Normalized green-blue difference index，NGBDI)和归一化绿红差异指数(Normalized green-red difference index，NGRDI)的植被覆盖度提取模型，采用支持向量机(Support vector machine，SVM)监督分类结果作为真值对各模型进行精度验证。
结果4种模型中，利用VDVI植被覆盖度提取模型获取的植被覆盖度精度最高，提取效果较好。与监督分类结果对比，4种植被覆盖度提取模型的提取误差(EF)分别为3.36%、15.68%、8.74%和15.46%，R2分别为0.946 1、0.934 4、0.695 3和0.746 0，均方根误差(RMSE)分别为0.021 9、0.059 5、0.042 0和0.055 9。
结论采用可见光植被指数结合像元二分法构建植被覆盖度提取模型实现了冬小麦返青期植被覆盖度准确快速提取，为植被覆盖度提取提供了一种新途径，可为无人机遥感监测提供参考。
Abstract:ObjectiveTo quickly and accurately extract vegetation coverage information of winter wheat in turning green period.
MethodThe field visible light images of winter wheat were obtained by UAV, and four common visible light vegetation indices were extracted. According to the principle of pixel dichotomy model, the vegetation coverage extraction models were established based on visible-band difference vegetation index(VDVI), excess green (EXG), normalized green-blue difference index (NGBDI) and normalized green-red difference index (NGRDI) respectively. The accuracies of four models were verified using the support vector machine (SVM) supervised classification results as the truth values.
ResultThe VDVI vegetation coverage extraction model had the highest accuracy and the best extraction effect to extract vegetation coverage among the four models. Compared with the supervised classification results, the extraction errors (EF) of four vegetation coverage extraction models were 3.36%, 15.68%, 8.74% and 15.46% respectively, the values of R2 were 0.946 1, 0.934 4, 0.695 3 and 0.746 0 respectively, and the values of root mean square error (RMSE) were 0.021 9, 0.059 5, 0.042 0 and 0.055 9 respectively.
ConclusionThe vegetation coverage extraction model based on visible vegetation index and pixel dichotomy has realized accurate and rapid extraction of vegetation coverage of winter wheat in turning green period, which provides a new way to extract vegetation coverage and a reference for UAV remote sensing monitoring vegetation coverage information.
-
Keywords:
- UAV remote sensing /
- winter wheat /
- pixel dichotomy /
- visible vegetation index /
- vegetation coverage
-
稻田杂草抑制水稻生长[1]，增加水稻病虫害发生概率[2]，是影响水稻产量和品质的主要原因[3-4]。因此，有效的稻田杂草防控一直是水稻安全生产的重要保障。近年来，机械除草作为一种环境友好型的除草技术被广泛地应用于水稻生产中，多用于去除行间杂草[5]。与行间杂草相比，水稻株间杂草多为近株杂草，对水稻生长的影响更大，但由于作物行的干扰，如何避开水稻苗去除株间杂草是水稻机械化除草技术的难点。
国内外现有的水稻株间机械除草装置多利用移栽稻和杂草的生物力学特性、苗草根系深浅差异设计柔性除草部件或控制除草深度实现株间除草。日本三菱农机株式会社研制的水田除草机[6]，采用柔性毛刷除去株间杂草，田间试验表明该类型除草机能有效提高水稻产量和品质[7]。欧莱克株式会社研制了一种四轮乘坐式除草机，行间采用滚轮除草，株间则采用弹性钢丝除草[8]。美善株式会社研制了一种独轮的手扶式除草机，行间采用滚轮除草，株间则使用从动转动式的羽轮除草，兼备除草和开沟的功能[9]。王金武等[10]研制的3ZS-150型水稻中耕除草机采用耕耘锄去除行间杂草；采用软轴驱动立式除草装置，利用弹性钢丝旋转去除株间杂草。以上水稻株间除草装置均不具备水稻植株和杂草的识别与定位功能，株间除草率仅为60%左右[11]，需多次作业才能完成除草要求，但多次作业不仅降低了生产效率，而且还增加了机具对秧苗的碾压概率。
为了降低伤苗率、提高除草效果，株间除草技术常采用传感器定位作物或杂草位置，选择性地避开秧苗，从而去除杂草。有研究采用实时差分GPS技术定位作物位置并精确地控制除草装置避开作物[12-13]，但该技术的硬件成本相对较高，而且不适用于水稻等小株距作物。日本石井农机株式会社研制的水稻株间除草机采用距离传感器识别稻株并控制电磁阀通断，喷射高压液体除草[14]。采用机器视觉识别植株的避苗技术[15-16]具有成本较低和控制系统相对简单的优点，蒋郁等[17]和Jiang等[18]根据水田环境特点设计了基于机器视觉的水稻植株定位系统，用以引导除草机构避开水稻苗实现株间除草。
除了与作物或杂草的定位精度有关，株间机械除草效果还和除草机械装置的结构参数和运动参数关系密切[19]。牛春亮等[20]对水田株间除草部件的运动轨迹进行分析，为稻田机械除草整机的研制提供了技术参考。贾洪雷等[21]通过对玉米中耕避苗除草装置的运动速度与避苗区域半径的关系进行分析，验证控制系统响应时间是否满足要求。在运动学分析基础上进行机械与水田土壤耦合的动力学分析[22-23]，能够直观地揭示各种不同结构及参数对除草部件工作性能的影响，从而对除草部件的局部结构进行更有针对性的优化。
在前期水稻植株定位方法研究的基础上，本研究拟设计一种气动式水稻株间机械除草装置，并对其进行运动学和动力学分析，优化除草装置结构和工作参数，以期提高水稻株间的除草效果。
1. 气动式株间除草装置工作原理及总体结构
1.1 株间机械除草装置工作原理
气动式水稻株间机械除草装置根据秧苗动态图像的定位结果作出相应的动作决策控制。图像采集和处理系统的主要硬件有：智能彩色相机(ISM1400C，美国Cognex)、工业镜头(M0814-MP，日本Computer)和上位机(B470，中国Lenovo)。信号转换及控制系统的主要硬件为I/O模块(CIO-MICRO-CC，美国Cognex)和可编程逻辑控制器PLC(FX2N-32MT，日本Mitsubishi)。
机具前进时，装置处于除草状态，当水稻植株图像进入检苗区域后，采用机器视觉定位[17]对秧苗进行动态定位；上位机发送定位信息的数字信号至智能相机，智能相机将数字信号转换成编码电信号后传输到I/O模块，I/O模块解码信号并输出信号至PLC，PLC运行内部程序触发输出端口；除草刀齿进入秧苗保护区后，发出切换电磁阀位置的指令，气缸杆伸出带动除草部件摆动，从而实现避苗动作；绕过稻株后，PLC再次发出切换电磁阀位置的指令，收回气缸杆，装置重新恢复至除草状态。控制指令将一直循环，直至工作结束。株间除草作业流程如图1所示。
1.2 株间除草装置总体结构
气动摆动式株间除草装置由气源(ZB-0.13/8, 中国浙江日金实业有限公司)、电磁阀(4v210-08, 中国Airtac)、气缸(SC50×25)、连杆、摆杆、除草刀齿和控制系统组成，具体组成如图2所示。除草装置随牵引机具向垂直纸面方向移动，安装在摆杆上的除草刀齿沿前进方向在株间入土破坏杂草，达到株间除草的效果。当除草刀齿到达秧苗根部保护区域内时，电磁阀动作使气路换向，气缸杆伸出推动连杆带动摆杆转动，安装在摆杆末端的除草刀齿随摆杆转动离开株间土壤完成避苗动作。绕过稻株后，电磁阀再次使气路换向，气缸杆回缩拉动连杆带动摆杆使除草刀齿入土执行株间除草动作。
![图 2 气动式株间除草装置示意图]()
图 2 气动式株间除草装置示意图1：气源；2：电磁换向阀；3：节流调速阀；4：气缸；5：节流调速阀；6：连杆；7：摆杆；8：除草部件Figure 2. Sketch map of pneumatic intra-row weeding device1: Air supply; 2: Solenoid directional valve; 3: Throttle; 4: Cylinder; 5: Throttle; 6: Connecting rod; 7: Oscillating rod; 8: Weeding component2. 气动式株间除草机构的运动分析
气动式株间除草机构运动简图如图3所示，根据机构运动原理，机构平面自由度(F)的计算如下：
![图 3 气动式株间除草机构运动简图]()
图 3 气动式株间除草机构运动简图v为气缸伸出速度(m/s)；a为连杆长(mm)；b为摆杆长(mm)；c为除草部件到回转中心水平距离(mm)；d为除草部件到回转中心垂直距离(mm)；w为除草部件横向摆动距离(mm)；h为除草部件提升高度(mm)；α为摆杆转角(°)；β为辅助角度(°)Figure 3. Kinematic diagram of mechanism of pneumatic intra-row weeding devicev is velocity of cylinder(m/s); a is length of connecting rod(mm); b is length of oscillating rod(mm); c is horizontal distance from weeding blade to rotation center(mm); d is vertical distance from weeding blade to rotation center(mm); w is oscillating distance of weeding blade at horizontal direction(mm); h is oscillating distance of weeding blade at vertical direction(mm); α is angle of oscillating rod(°); β is auxiliary angle(°)$$F = 3n - 2P,$$ (1) 式中，n为杆件和滑块数量，n=3，P为铰链转动副和移动副的数量，P=4。根据式(1)计算出F=1，说明机构在平面内有唯一确定的运动轨迹。除草刀齿中心在平面内的运动轨迹为绕点O作圆周摆动，摆动的幅度与角速度和气缸的运动行程与速度相关。假设气缸伸出时，推杆滑块作速度为v的匀速直线运动，在气缸伸出的任意时刻(t)，连杆长(a)、摆杆长(b)、摆动角度(α)和辅助角度(β)的几何关系根据余弦定理有：
$$\begin{split}{a^2} =& {b^2} + {\left[\sqrt {{b^2} + {{(a - vt)}^2}} \right]^2} -\\ &2b\sqrt {{b^2} + {{(a - vt)}^2}} \cos (\alpha + \beta ),\end{split}$$ (2) 其中，根据直角三角形的正切值关系有：
$$\beta = \arctan \frac{{a - vt}}{b},$$ (3) 将式(3)代入式(2)，则摆动角度(α)与时间(t)的关系可求得为：
$$\alpha = \arccos \left[\frac{{2{b^2} + {v^2}{t^2} - 2avt}}{{2b\sqrt {{b^2} + {{(a - vt)}^2}} }}\right] - \arctan \frac{{a - vt}}{b},$$ (4) 用摆动角度对时间求导，可得任意时刻除草刀齿的摆动角速度(ω)：
$$\omega = \frac{{{\rm{d}}\alpha }}{{{\rm{d}}t}},$$ (5) 通过式(2)~(5)可得连杆长和摆杆长与摆动杆的摆动角度和角速度的关系。根据几何关系，可得任意时刻除草刀齿中心提升高度(h)和横向摆动距离(w)：
$$h = d - \sqrt {{c^2} + {d^2}} \cos \left(\alpha + \arctan \frac{c}{d}\right),$$ (6) $$w = \sqrt {{c^2} + {d^2}} \sin \left(\alpha + \arctan \frac{c}{d}\right) - c{\text{。}}$$ (7) 气动式株间除草机构要有效地绕过秧苗，依照入土深度最大为30 mm计算，在避苗状态时规定除草部件的中心提升高度要≥40 mm，此时除草部件完全离开土壤，可避免损伤秧苗根部。茎基部保护区直径为60 mm的情况下，规定除草部件横向移动距离≥60 mm，此时可以避免碰伤秧苗的茎部。同时，为了提升传动的效率，摆杆转角≤20°[24]；为避免机架碰撞秧苗茎部，除草部件到回转中心垂直距离≥170 mm。
本研究中，先设定摆杆转角为20°，由于气缸行程为25 mm，则有vt=25 mm，根据式(4)可得连杆长(a)和摆杆长(b)的关系如图4所示，图中b对a的变化不敏感(a增大20 mm，b只增大0.5 mm)，故只需考虑a。为使机具整体结构紧凑和杆件轻量化，a应当尽量小，同时为避免杆件过小导致铰链与摆动杆干涉，设计a=35.00 mm，此时可求得b=72.24 mm。
![图 4 连杆长与摆杆长的关系]()
图 4 连杆长与摆杆长的关系Figure 4. Relationship of length of connecting rod and length of oscillating rod将式(6)和式(7)移项后求平方和，得到除草部件到回转中心水平距离(c)和除草部件到回转中心垂直距离(d)关于除草部件横向摆动距离(w)和提升高度(h)的关系式：
$${h^2} + {w^2} = 2dh - 2wc,$$ (8) 暂设h=40 mm，w=60 mm，可得到：
$$d = \frac{{130 + 3c}}{2},$$ (9) 将式(9)分别代入式(6)和式(7)可知，当c=84 mm，d=191 mm时，h=40.25 mm，w=60.26 mm，均符合设计要求，且此时除草部件到回转中心水平距离和除草部件到回转中心垂直距离处于最小状态，满足紧凑和轻量化的设计需求。根据计算，本研究设计气动式株间除草机构a=35 mm，b=72.24 mm，c=84 mm，d=191 mm。
为研究设计出的气动式株间除草机构的运动轨迹是否满足规定需要，将气动式株间除草机构的三维模型导入Pro-Engineer软件，运用机构运动仿真功能得到除草刀齿中心的运动轨迹。输入推杆滑块速度0.25 m/s，机具前进速度为0.25 m/s，可得除草刀齿中心的运动轨迹(图5)，除草刀齿能绕过护苗区，设计参数符合要求。
3. 除草刀齿设计与仿真优化
3.1 除草刀齿结构设计
除草刀齿摆动入土后通过破坏杂草或扰动杂草根系附近的土壤达到除草目的，故除草刀齿要求入土阻力小、入土后对土壤扰动大且不能出现壅土。根据上述要求设计出直齿刀和弯齿刀2种除草刀齿，如图6所示。为覆盖秧苗苗带区域，除草刀齿的作业宽度设计为60 mm。根据经验，刀齿数过多易引起壅土，过少则导致除草率低，综合考虑刀齿数量设计为4片。2种刀齿入土过程中先接触土壤的部分设计为锐角尖端，使其碰触土壤时更容易侵入破坏土壤表面从而更容易入土，刀片外形设计能减小刀齿平移方向的土壤阻力，同时依靠刀片侧面黏附和扰动土壤达到除草的目的。直齿刀的刃部从根部到尖端过渡形状为直线，弯齿刀的过渡形状则为曲线。2种刀刃的设计具有如下特点：直齿刀的直线触土面容易产生较大的黏附效果，除草效果好但阻力较大；而弯齿刀的曲线触土面对土壤具有一定的翻转作用且拥有较优的力学性能[25]。刀刃根部轮廓切线与垂直平面逆时针所成的角度称为倾角。为确保土壤不沿刀刃向上积累引起壅土，应倾角≥0°，同时应避免倾角过大导致刀片尺寸过长，故设计倾角的范围为0~20°。根据水田土壤高湿高阻的工作环境，关键除草部件的材料选用抗氧化性、耐腐蚀性和耐磨性的20CrMnTi号钢，采用数控电火花机床进行加工，考虑除草刀片要有足够的机械强度，刀片厚2 mm。
3.2 除草刀齿−水田土壤相互作用的仿真分析
为了模拟并分析2种除草刀结构(直齿刀和弯齿刀)和刀齿倾角(0、5、10、15和20°)对前进阻力和土壤的扰动程度的影响，得到较小前进阻力和较大土壤扰动程度的因素组合，运用离散元动力学仿真方法[26-27]，并参考相关文献[28]的方法，在EDEM软件中建立2种尺寸的球形颗粒模型分别代表水和土壤，将水颗粒覆盖在土壤颗粒上并使其自然向下渗透，模拟自然泥水的形成过程。颗粒模型采用Hertz-Mindlin with JKR Cohesion模型建立，该模型能够很好地模拟水田土壤颗粒间因静电和水分等原因发生明显黏结和团聚的物理特性。所建立的水田土壤模型，上部分为水颗粒层，中间为水颗粒与土壤颗粒混合的饱和土壤层，下层为有少量水颗粒的土壤层。广东省的水田土壤以砂质黏土为主，根据相关文献[23,28-29]确定水田土壤参数与水颗粒参数，通过EDEM官方材料数据库(GEMM数据库)获取除草刀具材料参数。土壤与水离散元颗粒和除草刀模型的材料参数如表1所示，材料间接触参数如表2所示。
表 1 材料参数Table 1. Material parameters材料
Material泊松比
Poisson’s ratio剪切模量/MPa
Shear modulus密度/(kg·m−3)
Density颗粒半径/mm
Particle radius颗粒接触半径/mm
Contact radius土壤 Soil 0.500 11.5 1830 0.5 0.7 水 Water 0.500 100.0 1000 0.2 除草刀 Weeding blade 0.288 2.8×104 7860 表 2 材料间接触参数Table 2. Contact parameters between materials材料
Material恢复系数
Coefficient of restitution静摩擦系数
Coefficient of static friction动摩擦系数
Coefficient of rolling friction水−水Water-water 0.05 0.05 0.01 土壤−土壤 Soil-soil 0.05 0.90 0.01 土壤−水 Soil-water 0.05 0.05 0.01 土壤−除草刀 Soil-weeding blade 0.10 0.20 0.20 水−除草刀 Water-weeding blade 0.10 0.20 0.20 对不同结构与不同倾角的除草刀齿与水田土壤的相互作用进行离散元动力学仿真。除草刀的前进速度设置为0.2 m/s，最大入土深度不宜超过移栽秧苗的根系深度[30]，刀齿最大入土深度为3 cm。土壤阻力通过土壤阻力−时间折线图(图7)表明，规定除草刀齿与土壤接触的时刻为0，仿真时间共0.5 s。
根据仿真结果，不同倾角(0、5、10、15和20°)的直齿刀对应的阻力最大值分别为24.18、19.38、15.55、21.89和19.20 N，对应的阻力平均值为14.30、13.53、12.94、15.82和15.49 N；而不同倾角(0、5、10、15和20°)弯齿刀对应的阻力最大值为5.78、5.69、4.98、5.20和5.73 N，对应的阻力均值为4.04、3.63、3.12、3.40和3.54 N。直齿刀所受的土壤阻力大于弯齿刀，不同类型的刀齿受土壤的阻力最大值与平均值均随着倾角增大呈先减后增的趋势。综上，倾角为10°的弯齿刀破坏土壤时受到土壤的阻力较小。
利用EDEM的后处理器，得到土壤颗粒的运动速度分布，土壤颗粒受除草刀扰动会产生运动速度，可根据产生运动速度的颗粒数量判定土壤受扰动程度的大小。垂直于刀齿前进方向，在刀齿尖端处做截面，观察同一时刻(0.5 s)截面上颗粒的运动速度分布，如图8所示。颗粒的运动速度范围为0~0.100 m/s，在刀齿附近的土壤颗粒为扰动土壤颗粒，速度分布在0.033~0.100 m/s，因此规定速度大于0.033 m/s的颗粒为受扰动的颗粒。使用Photoshop图像处理软件统计出截面内受扰动颗粒像素所占的面积，面积越大，说明土壤受扰动程度越大。统计结果见表3，根据统计结果，倾角为10°的弯齿刀产生的土壤扰动程度较大，被扰动的土壤区域面积为149.69 cm2。
![图 8 不同除草刀齿不同倾角作用下水田土壤颗粒运动速度分布]()
图 8 不同除草刀齿不同倾角作用下水田土壤颗粒运动速度分布Figure 8. Distribution of particle motion velocity under different weeding blade claw with different inclinations in paddy soil表 3 受扰动颗粒统计结果Table 3. Statistical result of particles under turbulence除草刀 Weeding blade 倾角/(°) Tilt angle 像素数量 Number of pixels 图像面积/cm2 Image area 直齿刀
Straight-blade0 70201 87.35 5 67096 83.49 10 58381 72.65 15 54773 68.16 20 48267 60.07 弯齿刀
Bending-blade0 82182 102.30 5 108747 135.37 10 120277 149.69 15 111040 138.22 20 102145 127.14 综合仿真结果可知，倾角为10°的弯齿刀与土壤相互作用时产生的阻力较小、土壤扰动程度较大，是较优的设计。
3.3 仿真结果验证
为验证离散元动力学仿真结果是否与实际情况相符，于2018年11月在华南农业大学工程学院土槽试验平台对弯齿刀进行测试试验，每种除草刀进行3次重复试验，每次试验后将土槽泥土恢复平整，测试不同倾角弯齿刀除草部件作业时的土壤阻力。具体设计如图9所示，土槽试验平台由土槽、滑轨支架和移动平台组成。土槽在试验前进行泡水、搅拌、混合等处理，并静置2 d，形成上部覆有水层的仿饱和水田土壤；移动平台由130BYG350D 步进电机(步距角1.2°，步距精度5%) 提供前进动力，测试过程中移动平台的速度为0.2 m/s，除草刀的入土深度为30 mm；除草刀和DS2-500N-S推拉力传感器(测量范围0~500 N，精度为±0.1%F.S，最小读数达到0.01 N)安装在移动平台上，拉力传感器的采集频率为50 Hz。由于平台匀速移动，拉力传感器测得的牵引力与土壤阻力大小相等、方向相反，所测的牵引力经转换后以数字信号的形式存储于计算机。
![图 9 仿真结果土槽验证试验]()
图 9 仿真结果土槽验证试验1：电控箱；2：计算机；3：测力计；4：步进电机；5：连接装置；6：拉力传感器；7：除草刀Figure 9. Soil bin validation experiment of simulation result1: Electronic control cabinet；2: Computer；3: Dynamometer；4: Stepper motor；5: Mounting device；6: Sensor of traction；7: Weeding blade土壤阻力实测值结果如图10所示，弯齿刀接触土壤后所受的阻力呈平稳上升趋势。在刀齿完全入土后土壤阻力在2~7 N稳定波动，且较平稳，数据变化趋势与仿真结果相似。试验过程中不同倾角弯齿刀所受最大阻力和平均阻力的实测值见表4，结果(表4)显示，弯齿刀的实测值大于仿真值，验证结果与仿真值误差小于20%，土壤阻力仿真结果可靠。
表 4 弯齿刀土槽试验实测数据与仿真结果Table 4. Measured data of bending blade experiment in soil bin and the simulation result倾角/(°)
Inclination最大阻力/N Maximum soil resistance 平均阻力/N Average soil resistance 实测值
Test data仿真值
Simulation data相对误差/%
Relative error实测值
Test data仿真值
Simulation data相对误差/%
Relative error0 6.80 5.78 17.65 4.62 4.06 14.36 5 6.30 5.69 10.72 4.31 3.67 18.73 10 5.00 4.98 0.40 3.72 3.16 19.23 15 6.17 5.20 18.65 3.98 3.40 17.06 20 5.20 5.68 9.25 4.20 3.54 18.64 ![图 10 不同倾角弯齿刀的土壤阻力实测值]()
图 10 不同倾角弯齿刀的土壤阻力实测值Figure 10. Measured soil resistance to bending-blade claw in variant inclination4. 田间试验
4.1 试验材料
株间除草装置结构优化(选用弯刀式、倾角10°的除草刀齿)后，为验证在不同工作参数下的综合工作性能(除草率及伤苗率)，并选出较优的工作条件，2019年4月5日在广东省肇庆市农业科学研究所试验田中进行机械除草单体搭载的田间性能试验。试验于插秧后15 d进行，试验使用的水稻品种为‘恒丰优3512’，秧苗平均高度为270 mm、株距约为150 mm。试验田内杂草以稗草、千金子、节节菜和陌上菜为主，杂草的平均根系深度为22 mm，株间杂草平均密度为每平方米160株。田块于试验前1 d蓄水30 mm，试验当天放水，试验时泥脚深度约为200 mm。田间试验现场如图11所示。
![图 11 株间除草装置田间试验现场]()
图 11 株间除草装置田间试验现场Figure 11. Scene of field experiment of paddy intra-row weeding device4.2 试验方法
株间除草装置单体的田间性能试验以气缸伸缩速度、机具前进速度和除草深度3个工作参数作为试验因子，采用三因素五水平二次旋转正交试验设计，分析各工作参数对除草率(y1)和伤苗率(y2)的影响，除草率与伤苗率定义如下：
$$ {y_1} = \frac{{Z - S}}{Z} \times 100{\text{%}} , $$ (10) $$ {y_2} = \frac{I}{M} \times 100{\text{%}} , $$ (11) 式中，Z为试验区域内水稻行间杂草总株数，S为除草作业后试验区域内水稻行间剩余杂草株数，I为除草作业后试验区域内被压折、连根拔起和倒伏的损伤秧苗数，M为试验区域内的总秧苗数。试验因素水平编码表如表5所示。
表 5 试验因素水平编码表Table 5. Factor level coding table for experiment编码值
Coded
value气缸伸缩
速度(x1)/(m·s−1)
Cylinder
protruding speed机具前进
速度(x2)/(m·s−1)
Machine
forward speed除草
深度(x3)/cm
Weeding
depth+γ 0.43 0.28 2.84 +1 0.40 0.25 2.50 0 0.35 0.20 2.00 −1 0.30 0.15 1.50 −γ 0.27 0.12 1.16 选取23块30 m长的地块作为试验区域，每个试验区域前后各5 m长的区域作为加速和减速的缓冲区，以保证机器能以稳定的工作状态通过试验区域。由2名计数员分别统计各自测试区域内杂草数和伤苗数，对测得的试验数据取平均值并转化为相应的除草率与伤苗率。
4.3 试验数据及分析
试验结果如表6所示。根据表6的试验数据，分析各因子对除草率的影响结果的方差分析如表7所示，比较各因子的F值可知，除草率显著性的影响因子依次为除草深度(x3)、机具前进速度(x2)、气缸伸缩速度(x1)。将影响不显著的因子(x1)剔除，得到除草率(y1)的回归模型方程：
表 6 试验方案及结果Table 6. Experimental scheme and results试验编号
Test number气缸伸缩速度/(m·s−1)
Cylinder protruding speed机具前进速度/(m·s−1)
Machine forward speed除草深度/cm
Weeding depth除草率/%
Weeding rate伤苗率/%
Seedling injury rate1 0.30 0.30 1.50 61.24 2.40 2 0.40 0.30 1.50 60.55 2.02 3 0.30 0.25 1.50 66.89 4.55 4 0.40 0.25 1.50 67.94 2.76 5 0.30 0.30 2.50 72.98 3.77 6 0.40 0.30 2.50 72.16 3.35 7 0.30 0.25 2.50 81.13 4.04 8 0.40 0.25 2.50 81.17 3.85 9 0.27 0.000 2.00 75.39 4.23 10 0.43 0.000 2.00 74.03 2.90 11 0.35 0.12 2.00 59.12 2.49 12 0.35 0.28 2.00 80.59 4.25 13 0.35 0.20 1.16 51.76 1.84 14 0.35 0.20 2.84 78.93 3.93 15 0.35 0.20 2.00 72.92 3.27 16 0.35 0.20 2.00 75.77 3.43 17 0.35 0.20 2.00 76.14 3.03 18 0.35 0.20 2.00 72.37 3.11 19 0.35 0.20 2.00 73.71 3.27 20 0.35 0.20 2.00 73.96 2.93 21 0.35 0.20 2.00 72.24 3.52 22 0.35 0.20 2.00 71.90 3.17 23 0.35 0.20 2.00 74.06 3.48 表 7 各因子对除草率影响的方差分析表Table 7. Variance analysis of influence from each factor to weeding rate来源 Source 平方和 Sum of square 自由度 Degree of freedom 均方和 Mean square F P 模型 1157.79 9 128.64 31.59 <0.0001 x1 0.54 1 0.54 0.1318 0.7224 x2 321.95 1 321.95 79.05 <0.0001 x3 682.08 1 682.08 167.48 <0.0001 x1x2 0.85 1 0.85 0.2075 0.6563 x1x3 0.16 1 0.16 0.0399 0.8448 x2x3 2.12 1 2.12 0.5210 0.4832 $ x_1^2$ 4.37 1 4.37 1.07 0.3194 $ x_2^2$ 22.59 1 22.59 5.55 0.0349 $ x_3^2$ 123.41 1 123.41 30.30 0.0001 残差 Residual 52.94 13 4.07 失拟项 Lack of fit 34.77 5 6.95 3.06 0.0776 纯误差 Pure error 18.18 8 2.27 总和 Sum 1210.74 22 $$ {y_1} = 73.66 + 4.86{x_2} + 7.07{x_3} - 1.19x_2^2 - 2.79x_3^2{\text{。}} $$ (12) 根据表6的试验数据，各因子对伤苗率的影响结果的方差分析如表8所示，模型的F值为16.85，P<0.0001，表明该模型可信；失拟项F值为2.48，P=0.1215，差异不显著，表明该模型的拟合程度良好；因子x1、x2、x3影响极显著，因子x2x3在水平α=0.05时影响显著，因子x1x3、x12在水平α=0.1时影响显著，其他因子影响不显著；比较各因子的F值可知，影响伤苗率显著性的因子依次为除草深度、机具前进速度和气缸伸缩速度。将影响不显著的因子剔除，得到伤苗率(y2)的回归模型方程为：
表 8 各因子对伤苗率影响的方差分析表Table 8. Variance analysis of influence from each factor to seedling injury rate来源 Source 平方和 Sum of square 自由度 Degree of freedom 均方 Mean square F P 模型 Model 9.9900 9 1.1100 16.85 <0.0001 x1 1.8400 1 1.8400 27.98 0.0001 x2 3.2100 1 3.2100 48.71 <0.0001 x3 3.3800 1 3.3800 51.32 <0.0001 x1x2 0.1740 1 0.1740 2.64 0.1280 x1x3 0.3042 1 0.3042 4.62 0.0511 x2x3 0.5618 1 0.5618 8.53 0.0119 $ x_1^2$ 0.2560 1 0.2560 3.89 0.0703 $ x_2^2$ 0.0534 1 0.0534 0.81 0.3842 $ x_3^2$ 0.2046 1 0.2046 3.11 0.1015 残差 Residual 0.8563 13 0.0659 失拟项 Lack of fit 0.5207 5 0.1041 2.48 0.1215 纯误差 Pure error 0.3356 8 0.0420 总和 Sum 10.8500 22 $$\begin{split} {y_2} = &3.24 - 0.37{x_1} + 0.48{x_2} + 0.50{x_3} +\\ &0.20{x_1}{x_3} - 0.27{x_2}{x_3} + 0.13x_1^2{\text{。}} \end{split}$$ (13) 根据回归方程(12)，气缸伸缩速度对除草率的影响不显著，因此将气缸伸缩速度设定在零水平，即0.35 m/s，得到的响应曲面如图12所示。分析可知，除草率最低为60%，除草率随着机具前进速度和除草深度的增加而增加，两者达到最大时，除草率最大(82%)。通过比较发现，除草深度对除草率的影响更加显著。
![图 12 机具前进速度和除草深度对除草率影响的响应曲面]()
图 12 机具前进速度和除草深度对除草率影响的响应曲面Figure 12. Response surface of machine forward speed and weeding depth to weeding rate将除草深度设定为零水平，即2.0 cm，得到气缸伸缩速度与机具前进速度的伤苗率响应曲面(图13a)。分析图13a可知，伤苗率最低为2.73%，最高为4.39%，且伤苗率随着气缸伸缩速度的升高和机具前进速度的降低而降低。气缸伸缩速度升高会导致机具前进速度对伤苗率的影响减弱，同时机具前进速度的降低也会导致气缸伸缩速度对伤苗率的影响减弱。当气缸伸缩速度最高而机具前进速度最低时，伤苗率达到最小值。
![图 13 不同影响因子对伤苗率影响的响应曲面]()
图 13 不同影响因子对伤苗率影响的响应曲面Figure 13. Response surface of influence from each factor to seedling injury rate将机具前进速度设置为零水平，即0.20 m/s，得到气缸伸缩速度与除草深度的伤苗率响应曲面(图13b)。分析图13b可知，伤苗率最低为2.25%，最高为3.91%，且伤苗率随着气缸伸缩速度的升高和除草深度的降低而降低。气缸伸缩速度升高会导致除草深度对伤苗率的影响增强，同时增加除草深度也会减弱气缸伸缩速度对伤苗率的影响。当气缸伸缩速度最高而除草深度最低时，伤苗率达到最小值。
将气缸伸缩速度设置为零水平，即0.35 m/s，得到机具前进速度与除草深度的伤苗率响应曲面(图13c)。分析图13c可知，伤苗率最低为1.99%，最高为3.89%，且伤苗率随着机具前进速度和除草深度的增加而增加。当机具前进速度和除草深度同时最小时，伤苗率达到最小值。
根据以上的分析，结合机具实际工作的要求，设定机具前进速度为0.25 m/s，气缸伸缩速度为0.45 m/s，除草深度为1.0~3.0 cm。以除草率>80%，伤苗率<4%为目标，利用Design-Expert 11软件进行参数优化分析，在符合优化目标的参数中，除草深度的均值为2.42 cm，因此除草深度选择2.50 cm。
在国家水稻产业技术体系肇庆市农业科学研究所开展水稻株间机械除草性能多点重复田间试验，采用优化的工作参数，气缸伸缩速度0.45 m/s、机具前进速度0.25 m/s和除草深度2.5 cm。每次试验分别选取30 m长的地块作为试验区域，重复5次，试验结果取平均值。试验结果表明肇庆不同水稻品种的平均除草率为83.91%，伤苗率为3.63%，该除草装置能够满足水稻株间除草的性能要求。
5. 结论
本研究在稻株机器视觉识别定位技术研究的基础上，采用机械设计理论、离散元动力学仿真方法结合田间试验研制出气动式水稻株间机械除草装置，结论如下：
1)运用机械设计原理设计气动式株间除草机构及其关键参数，得到气动式株间除草机构的连杆长35 mm，摆杆长72.24 mm，除草部件到回转中心水平距离84 mm，垂直距离191 mm。通过理论计算及Proe运动学仿真验证其运动轨迹并测得除草部件提升高度40.25 mm>40 mm，除草部件横向摆动距离60.26 mm>50 mm，符合株间避苗除草的要求。
2)利用EDEM软件建立水田土壤泥水混合的颗粒模型，在此基础上进行不同结构和倾角组合的除草刀齿与水田土壤相互作用的离散元动力学仿真分析。根据仿真结果可得：关键除草部件形态为弯齿刀、倾角为10°时，关键除草部件与土壤接触时阻力较小，为3.12 N，且对土壤的扰动程度较大，受影响的面积达149.69 cm2。土槽试验验证结果表明土壤阻力实测值与仿真值结果一致，误差基本小于20%，仿真结果可靠。
3)对株间除草装置整机进行田间试验，试验采用二次正交旋转组合设计，得到气缸伸缩速度、机具前进速度和除草深度对除草率和伤苗率的影响规律，优化结果表明：在机具前进速度为0.25 m/s，气缸伸缩速度为0.45 m/s，除草深度为2.5 cm的条件下，平均除草率为83.91%，伤苗率为3.63%，满足株间除草率>80%，伤苗率<4%的性能要求。
-
![]()
图 5 不同植被覆盖度提取模型的线性拟合结果
Figure 5. The linear fitting results of different vegetation coverage extraction models
表 1 冬小麦分类精度评价
Table 1 Evaluation of classification accuracy of winter wheat
类别
Category小麦/像素
Wheat土壤/像素
Soil样本总数/像素
Total sample size用户精度/%
User accuracy小麦/像素 Wheat 20 779 3 20 782 99.99 土壤/像素 Soil 92 25 806 25 898 99.64 样本总数/像素 Total sample size 20 871 25 809 46 680 用户精度/% User accuracy 99.56 99.99 
下载: 导出CSV
表 2 植被覆盖度精度评价
Table 2 Accuracy evaluation of vegetation coverage
植被指数模型
Vegetation index model覆盖度 Vegetation coverage 提取误差/%
Extraction error像元二分法
Pixel dichotomy监督分类
Supervised classification差值
Difference valueVDVI 0.321 452 0.332 623 0.011 171 3.36 NGBDI 0.303 564 0.332 623 0.029 059 8.74 NGRDI 0.384 043 0.332 623 0.051 420 15.46 EXG 0.384 780 0.332 623 0.052 157 15.68
下载: 导出CSV
表 3 植被覆盖度等级分布统计结果
Table 3 Statistical results of vegetation coverage grade distribution
植被覆盖度
Vegetation coverage像元数
Pixel number百分比/%
Percentage0~0.10 1 960 500 8.15 0.10~0.30 12 265 775 51.01 0.30~0.45 4 166 478 17.33 0.45~0.60 2 460 841 10.23 0.60~1.00 3 192 675 13.28 总计 Total 24 046 269 100
下载: 导出CSV
-
[1] 赵广才, 常旭虹, 王德梅, 等. 中国小麦生产发展潜力研究报告[J]. 作物杂志, 2012(3): 1-5. [2] 刘治开, 牛亚晓, 王毅, 等. 基于无人机可见光遥感的冬小麦株高估算[J]. 麦类作物学报, 2019, 39(7): 859-866. [3] PUREVDORJ T, ISHIYAMA T, HONDA Y, et al. Relationships between percent vegetation cover and vegetation indices[J]. Int J Remote Sens, 1998, 19(18): 3519-3535. doi: 10.1080/014311698213795
[4] GITELSON A A, KAUFMAN Y J, STARK R, et al. Novel algorithms for remote estimation of vegetation fraction[J]. Remote Sens Environ, 2002, 80(1): 76-87. doi: 10.1016/S0034-4257(01)00289-9
[5] 牛亚晓, 张立元, 韩文霆, 等. 基于无人机遥感与植被指数的冬小麦覆盖度提取方法 [J]. 农业机械学报, 2018, 49(4): 212-221. [6] SONG W, MU X, RUAN G, et al. Estimating fractional vegetation cover and the vegetation index of bare soil and highly dense vegetation with a physically based method[J]. Int J Appl Earth Obs, 2017, 58: 168-176. doi: 10.1016/j.jag.2017.01.015
[7] RUNDQUIST B C. The influence of canopy green vegetation fraction on spectral measurements over native tallgrass prairie[J]. Remote Sens Environ, 2002, 81(1): 129-135. doi: 10.1016/S0034-4257(01)00339-X
[8] ZHOU Q, ROBSON M, PILESJO P. On the ground estimation of vegetation cover in Australian Rangelands[J]. Int J Remote Sens, 1998, 19(9): 1815-1820. doi: 10.1080/014311698215261
[9] 陈虹兵, 黄贝贝, 彭道黎. 基于混合像元分解的植被覆盖度模型比较研究[J]. 西北林学院学报, 2018, 33(3): 203-207. doi: 10.3969/j.issn.1001-7461.2018.03.31 [10] 佟斯琴, 包玉海, 张巧凤, 等. 基于像元二分法和强度分析方法的内蒙古植被覆盖度时空变化规律分析[J]. 生态环境学报, 2016, 25(5): 737-743. [11] 李苗苗, 吴炳方, 颜长珍, 等. 密云水库上游植被覆盖度的遥感估算[J]. 资源科学, 2004, 26(4): 153-159. doi: 10.3321/j.issn:1007-7588.2004.04.022 [12] 何勇, 岑海燕, 李艺健, 等. 基于无人机遥感的农田信息获取技术及其机载装备[J]. 农业工程技术: 农产品加工, 2018, 38(9): 28-32. [13] 汪小钦, 王苗苗, 王绍强, 等. 基于可见光波段无人机遥感的植被信息提取[J]. 农业工程学报, 2015, 31(5): 152-159. doi: 10.3969/j.issn.1002-6819.2015.05.022 [14] 赵静, 杨焕波, 兰玉彬, 等. 基于无人机可见光图像的夏季玉米植被覆盖度提取方法[J]. 农业机械学报, 2019, 50(5): 232-240. doi: 10.6041/j.issn.1000-1298.2019.05.027 [15] TORRES-SÁNCHEZ J, PEÑA J M, DE CASTRO A I, et al. Multi-temporal mapping of the vegetation fraction in early-season wheat fields using images from UAV[J]. Comput Electron Agr, 2014, 103: 104-113. doi: 10.1016/j.compag.2014.02.009
[16] 龙满生, 何东健. 玉米苗期杂草的计算机识别技术研究[J]. 农业工程学报, 2007, 23(7): 139-144. doi: 10.3321/j.issn:1002-6819.2007.07.027 [17] HUNT E R, CAVIGELLII M, DAUGHTRY C S T, et al. Evaluation of digital photography from model aricraft for remote sensing of crop biomass and nitrogen status[J]. Precis Agric, 2005, 6(4): 359-378. doi: 10.1007/s11119-005-2324-5
[18] MEYER G E, NETO J C. Verification of color vegetation indices of automated crop image application[J]. Comput Electron Agr, 2008, 63(2): 282-293. doi: 10.1016/j.compag.2008.03.009
[19] 陈晋, 陈云浩, 何春阳, 等. 基于土地覆盖分类的植被覆盖率估算亚像元模型与应用[J]. 遥感学报, 2001, 5(6): 416-422. doi: 10.11834/jrs.20010603 [20] ZRIBI M, HÉGARAT-MASCLE S L, TACONET O, et al. Derivation of wild vegetation cover density in semi-arid regions: ERS2/SAR evaluation[J]. Int J Remote Sens, 2003, 24(6): 1335-1352. doi: 10.1080/01431160210146668
[21] 张军国, 韩欢庆, 胡春鹤, 等. 基于无人机多光谱图像的云南松虫害区域识别方法[J]. 农业机械学报, 2018, 49(5): 249-255. doi: 10.6041/j.issn.1000-1298.2018.05.029
