基于田间摄像的多参数水稻土深松扰动行为与效应研究

丁启朔 李 杨 BELAL Eisa Adam 梁 磊 何瑞银 汪小旵

(1.南京农业大学工学院， 南京 210031； 2.江苏省智能化农业装备重点实验室， 南京 210031； 3.扎林盖大学农业工程系， 扎林盖 06)

0 引言

土壤深松已成为保护性农业的一项重要耕作技术[1]，通过加深耕作层、改善土壤孔隙、提高土壤透水和透气性实现作物增产[2-6]。在生产实践中，土壤深松与秸秆还田[7-9]、深松与轮耕系统设计[10-12]、深松保障作物增产等研究内容不断拓展[13-15]，兼顾各种复合效益的深松机具设计也在不断完善[16-19]。在深松理论上，有限元[20-21]和离散元[22-25]等仿真模型的应用也深化了人们对深松等土壤耕作力学机理的认识。尽管如此，不同田间条件下深松机具在作业过程中的一些扰动行为及其效应仍有待揭示，例如，深松机具在作业过程中的碎土和土壤抬升等行为需要准确描述，碎土与抬土行为造成的微地貌及土层紧实度变化更需要直观表达，只有如此才能定量表述基于触土部件结构设计和作业参数的土壤深松扰动行为，进而优化耕作系统。

深松过程的碎土与土垡抬升所造成的微地貌及土壤紧实度变化也进一步决定着是否需要配合二次整地等后续作业，因此关系到保护性耕作系统的设计。可见，揭示土壤耕作的过程行为也是当前的一个技术难题，涉及能否简化耕作流程、实现复式作业，甚至整合深松、播种施肥及秸秆还田的技术可行性，更关系到能否解决好保护性农业需要最大程度疏松土层，但同时要求节约能耗的问题，以及需要加深耕层、但同时要求不扰乱土层和减少微地貌破坏的现实性矛盾。

本研究采用土壤耕作原位综合测试平台进行单铲深松试验，从5个方位协同录制深松作业视频，以获取关键的水稻土深松扰动行为指标，结合深松前后的土壤微地貌和紧实度变化、耕作阻力测试与EDEM仿真模拟，定量表达不同入土角和耕深条件下的水稻土深松扰动行为及效应。

1 材料与方法

1.1 试验地概况

试验地位于南京市六合区八百桥(118°59′E，31°98′N)，该地区常年稻麦轮作，土壤类型为壤质粘土，0～30 cm土层土壤物理参数见表1，试验前清除地表留茬。

表1 0～30 cm土层基础物理参数Tab.1 Basic physical properties in 0～30 cm soil layer

1.2 原位深松摄像试验设计

田间试验使用土壤耕作原位综合测试平台[29]，挂接单铲，牵引速度为0.1 m/s(图1b)。试验前沿深松铲行进方向(x轴)建立坐标系(图1a)，5台摄像机按表2所示空间坐标布局。

图1 水稻土原位深松试验机位布置图Fig.1 Camera positions of subsoiling experiment

鉴于深松铲入土角α及耕深D对作业性能的显著影响[24,30-34]，参照文献[35]，确定α为15°、22.5°、30°；深松或锄耕(Chiseling，多见于西方保护性耕作机具)的耕深范围在10～30 cm[36-39]。另外生产中通常考虑到稻作制条件大于30 cm的深松泥脚过深，影响机具作业且破坏阻渗层，故确定D为10、20、30 cm，本试验双因素共9种处理(见表3)，参照标准JB/T 9788—1999《深松铲和深松铲柄》仿制一组深松铲用于试验。

表2 摄像机位空间坐标Tab.2 Camera space coordinates cm

表3 深松过程参数Tab.3 On-site parameters of subsoiling

注: 同列相同字母上标表示不同处理之间差异不显著(P>0.05); 小写字母表示入土角的显著性，大写字母表示耕深的显著性。

1.3 试验操作和测试方法

1.3.1深松过程的摄像法测试

深松过程的土壤扰动行为由5个机位协同拍摄。鉴于深松铲的抬土过程是周期性发生，仅以摄像中的抬土最高点静态图像提取深松过程行为指标，具体包括x、y、z方向上的纵向碎土距离R、扰动宽度W和抬土高度H。其中，H和R从机位3、4的视频中提取(图2a)，W从机位5的视频中提取(图2b)，视频处理使用Digimizer软件。每个处理重复6次，采用Excel 2010和SPSS 19.0软件对试验数据进行统计分析。图2上部组图为6次重复样本，下部组图为一次完整深松过程的动态展示，其中黄色边框为标注大图。

图2 土壤动态扰动参数示意图Fig.2 Soil dynamic disturbance parameters

1.3.2深松前、后微地貌测试

本文使用小蚁4K+运动摄像机(900万像素，视频帧率30 f/s)以0.1 m/s的速度分别对深松前、后土壤微地貌进行3次不同角度和高度的拍摄，利用Adobe-Premiere软件提取高精度视频帧(图3a)用于Agisoft-Photoscan软件三维重建，结合地面控制点坐标后获取数字化微地貌模型，利用Geomagic-Studio软件的网格修补模块对微地貌三维点云模型进行修复(图3b)，使用曲线模块提取土壤表面模型的样线作为地表线[40]，并在深松后地表线上过最高点做一条水平直线为基准线，在整个深松铲深松宽度范围内以固定值(30 mm)等分标记测点，在深松作业前、后根据比例采样获取各测点至水平基准线的垂直距离，计算出平均值和标准差，并以标准差表示地表平整度，上述试验分别测量5次取平均值[25]。摄像法测试完成后，采用刘俊安等[26]的方法测得深松扰动剖面的土壤紧实度等值线图。耕作阻力由土壤耕作原位综合测试平台测试记录[29,36]。

图3 耕后土壤微地貌Fig.3 Disturbed microrelief

垂直距离与平整度计算式为

(1)

式中ak——第k次测量的深松前、后地表至水平基准线的垂直距离

akj——第k次测量中第j个等分点的深松前、后地表至水平基准线的垂直距离

nk——第k次测量中的测试点数

Sk——第k次测量的地表平整度

2 结果与讨论

2.1 摄像法所得深松过程参数与效应

2.1.1摄像法所得深松过程参数

表3表明，入土角α和耕深D对纵向碎土距离R、抬土高度H和扰动宽度W均影响显著，表4汇总了水稻土深松的扰动行为指标与α和D间的拟合关系，表明H随α呈线性增加，而LI等[24]和张锐等[41]研究成果表明随着α增大，机具推土能力增强，而抬土能力减弱。视频的回放进一步表明水稻土的深松过程不以脆性断裂为主，而是绕着O点(图2a)，沿着曲柄段向上攀升，至直柄段被剪切回落。

表3、4所示α和D对R的影响规律与TONG等[42]和GODWIN等[43]的研究结果一致,R随对应耕深范围内铲的纵向水平长度L[44]线性增加。不仅如此，α和D对扰动宽度W也有显著影响。W随α线性增加，因为铲前方土壤位移场变大，扰动程度增强，新月形失效面宽度[43]随α增加，这与张锐等[41]和TONG等[42]的研究结果一致。

表4 水稻土深松作业性能指标与深松铲入土角、耕深的关系Tab.4 Regressional analysis on subsoiler performance parameters with respect to subsoiling parameters

注：H10、H20、H30表示耕深D为10、20、30 cm的抬土高度，H15、H22.5、H30表示入土角α为15°、22.5°、30°的抬土高度。其余类似。

2.1.2水稻土深松的微地貌特征

通过摄像法获得重构数字云微地貌，并从中提取地表平整度参数与针板法测量结果的平均相对误差为8.6%(表5)，表明数字云重构三维土壤模型能够高效、准确地提供微地貌信息。

表5 摄像法与针板法测试微地貌误差分析Tab.5 Error analysis on video measurement and pin-board measurement

地表平整度S是评价耕作质量的重要指标[25,28,45]，表3、5表明D对S有显著影响，且S随D先增大后减小，这一变化规律与以往研究所得结论一致[36]。

摄像法所得抬土高度H随D显著减小，但地表平整度变化趋势与此不同。视频显示在10 cm耕深时土壤发生新月形失效，碎土充分，土块间的粘接呈完全断裂破坏，土垡回落轨迹为从铲两侧翻落至地表，回落距离最大；随着耕深增加，耕作土层越加深厚致密，土壤失效模式转变为新月形失效和侧向挤压失效[43]，土壤剪切破碎不充分，土垡与两侧土壤粘连，土垡回落距离显著减小，这也进一步支撑了针对土层状态来设计深松铲的理论假说[46]。

2.1.3摄像法所得土壤剖面扰动过程

图4 土壤剖面扰动过程Fig.4 Soil disturbance process recorded at stations 1 and 2

从机位1、2的视频按时间序列提取土壤扰动过程(图4)，其中2 s时铲尖入土，3.7 s时铲尖“挤”进土层且底层土沿水平方向产生扩展裂纹，同时铲尖上方的土壤开始向上运动，铲尖下方的土壤则因为深松铲的挤压向下运动。4.6 s时，剪切失效的土体沿深松铲曲柄段抬升，被抬升的土体向前和向上运动，同时土壤剖面的裂纹继续向侧斜上方延伸,摄像显示的土壤剖面失效过程与已有的研究成果一致[21,23,27,47-49]；铲尖和曲柄附近的底层土受到挤压剪切形成沟槽。5.3 s时被抬升的土体受深松铲直柄和曲柄段剪切产生土垡和裂纹，但摄像显示仅有铲右侧的土垡被继续剪切抬起，土体失效过程呈现出非对称性，前人研究成果表明凿形铲等窄齿耕作部件作业形成的地表裂纹呈左右交替性分布[30,36],这与本文研究结果一致。6.3 s时深松铲已经穿越了抬升的土体，大部分土垡完成回落，仅有少量土垡黏附在深松铲上继续前移。可见，摄像法直观展示了水稻土深松过程中，深松铲结构与土壤扰动行为的紧密联系，可为今后深松铲的设计和优化提供决策依据。

图5 土壤扰动过程分析Fig.5 Simulation of soil disturbance processes

2.2 水稻土深松的EDEM触土过程分析

近年来大量的离散元(Discrete element method,DEM)模型被用于解析土壤深松的发生机理[23-25,27,41]，本文沿用经过试验标定和检验的湿粘水稻土EDEM模型仿真方法[49]，获取水稻土深松触土过程时间序列图(图5)，分别从纵向和横向剖面2个视窗同步解析摄像法记录的结果。

前人的研究成果证明离散元仿真能够准确预测耕作阻力[23-25]，本文耕作阻力F的仿真值变化趋势与试验值一致，平均相对误差为6.6%,EDEM仿真所得耕作阻力符合实测结果(表6)。但更关键的是DEM便于从微观颗粒层面解析耕作时土壤的运动与受力[41]。由图5可知，3.7 s时铲尖附近的底层土壤因为被挤压获得速度，耕作层部分土壤随底层土抬升获得速度，铲尖上方的土壤向上运动，铲尖下方的土壤则向下运动，抬升区和压实区开始形成，应力集中区主要存在于铲尖周围。

图6 深松扰动剖面紧实度等值线图Fig.6 Contours of soil cone index of subsoiling trench profile

入土角α/(°)耕深D/cm耕作阻力试验值/kN耕作阻力仿真值/kN耕作阻力相对误差/%100.540.599.2615.0201.902.1312.10304.534.744.64100.670.681.4922.5202.212.439.95305.175.302.51100.840.884.7630.0202.632.889.51305.645.935.14

4.6 s时主要为曲柄作用阶段，临界深度(15 cm)以上疏松耕作层发生新月形失效，以下的致密底层土发生侧向挤压失效[30,36,43,50-51]。红色土壤颗粒所呈现的扰动边界轮廓线与摄像法实测相符。铲周边土壤速度最大，等速度土壤颗粒分布曲线与铲柄的弧线基本吻合，这与黄玉祥等[27]的研究结果一致。其中，新月形失效区的土壤速度方向主要为向上、向前和向两侧，侧向失效区的土壤速度方向主要为向前和向两侧。应力集中区从铲尖附近沿着曲柄破土刃[52]向上延伸，因为此时曲柄正剪切抬升土层。

5.3 s时铲右侧的失效土壤继续向上运动且获得了更大的速度，而左侧土壤的运动速度与4.6 s无差异(图5b)。横向扰动过程显示的土壤非对称性失效进一步表明离散元模型用于解析土壤耕作力学的科学性。另外此时应力集中区并没有继续上移，反而下移至铲尖周围的土壤，这是因为待剪切破碎的土体经过曲柄的抬升已较为松散，剪切力较小，而铲尖此时仍在挤压底部的致密土体。

6.3 s时深松铲已经完成了对土壤的局部扰动，剪切形成的土垡大部分回落至地表和深松沟，少量黏附在深松铲上，铲尖此时已进入到下一周期的“挤压入土”阶段，形成摄像法观测到的循环往复抬土和碎土现象。

EDEM仿真显示在一个完整的深松作业周期内，土壤的应力集中区首先发生在铲尖附近，随后上移至曲柄附近，之后又回到铲尖附近。从“动态设计”的角度上看，触土曲面前方应力集中区的周期性波动对于土壤的深松效果和耕作部件的力学性能的提升是有利的[51]，这同样为振动深松实现节能降阻提供一定的理论支撑[18,53-54]。

2.3 实测深松扰动剖面土壤紧实度

上述摄像法、EDEM仿真及力学测试在用于解析水稻土深松过程行为和土壤失效机理时各有技术优势，但尚无法反映出深松扰动剖面的土壤紧实度分布状况(图6)，表明今后针对田间原位土壤深松研究需要综合运用多种测试方法。

前人研究成果表明深松铲的设计参数和作业参数对土壤孔隙度和紧实度的影响并不一致[21,24,26,55-56]。颜色较浅的区域为深松扰动区，与实测深松沟型相符。图6表明，α对深松扰动剖面紧实度分布无显著影响，这与LI等[24]和刘俊安等[26]的研究成果一致；D对其影响显著，这与LI等[24]和RAHMAN等[31]的研究成果一致。

试验前多次测得0～15 cm耕作土层的紧实度平均值在400 kPa附近，故选取深松后紧实度等值线图中小于400 kPa的区域为有效深松区。分土层分析得知，对于0～10 cm土层，3种耕深的紧实度等值线分布基本一致，这是因为熟土层无法显现出深松的作用[56]，因此尽可能降低深松铲在该土层的扰动是保护性深松系统设计的关键；对于10～20 cm土层，当耕深为20 cm时，15°、22.5°、30°入土角所对应的有效深松区面积依次为15 146.08、16 373.98、16 997.62 mm2，当耕深为30 cm时，3种入土角所对应的有效深松区面积依次为7 488.43、8 223.22、8 874.58 mm2，可见耕深为20 cm时，深松铲提供了关键土层(犁底层)的最大扰动；对于20～30 cm土层，当耕深为30 cm时，3种入土角所对应的有效深松区面积依次为1 480.59、1 488.52、1 947.38 mm2，土壤扰动效果显著下降，这是因为该处土层位于临界深度以下，发生侧向挤压失效，紧实度等值线图呈现出狭窄的沟槽形态(图6d)，这不仅表明铲对此处土层的扰动能力急剧减弱，而且反映出水稻土深松的耕作力学特点。

所有处理的铲尖下方都出现土壤紧实度等值线密集分布的现象，说明此处土壤紧实度变化剧烈，产生土壤压实，这与前人研究成果一致[44,57-58]，且这一现象在EDEM仿真和机位2的视频回放中得到验证。

2.4 结合深松理论的对比检验

2.4.1土壤扰动过程理论检验

结合深松铲功能结构和土壤条件的扰动过程分析是对摄像法、耕作阻力测试、EDEM仿真、扰动剖面土壤紧实度测试的必要检验，可对上述试验结果提供宏观的土壤耕作力学解析。以深松铲为参考系对曲柄圆弧进行微元分析得到抬土高度H(图7a)为

H=L′tanε

(2)

式中L′——扰动土体攀升过程中，任一时刻沿水平方向的运动距离

ε——曲柄圆弧切线与水平线所夹锐角

因本文所用结构参数α1<α2<α3，即ε1<ε2<ε3，故抬土高度H1

图7 不同耕深作业情况Fig.7 Performances of different subsoiling depths

稻作区的浅薄耕层源自长期浅旋耕作、反复的机械压实及稻季淹水沉浆过程，大范围的区域调查表明水稻土的有效耕层仅为11～17 cm[59]，基于本试验地的土壤剖面分析也发现稻季土壤排干后的有效耕层在15 cm左右[49]。当耕深为10 cm时，剪切的土层整体容重和剪切模量较小，土体的自由移动和动力条件优越[60]，易剪切抬升，另外地表以上的曲柄圆弧较长，扰动土体有足够的行程沿曲柄破土刃攀升(图7a)；当耕深增大时，剪切的土层整体容重和剪切模量增大，地表以上曲柄弧长减小，土层抬升难度增大的同时又减小了其攀升行程，抬升的土层很快到达直柄处被剪切回落，因此如摄像所示，抬土高度H随耕深D减小。

基于深松铲-土壤系统相互作用的探讨也能够解释W的变化规律，随着耕深增加到20 cm，深松铲的扰动土体量增加，深松铲的触土范围也随之从铲尖和曲柄下部演变成铲尖和全部的曲柄，铲的扰动能力得以充分发挥。当耕深增大至30 cm时，相当部分的深松铲功能区淹没在致密紧实的底层土(图7c)，耕作阻力急剧增大，铲的扰动能力被大大削弱[49]，W减小。

增大抬土高度有利于充分利用耕作能量疏松土层，但是过大的抬土高度往往会增大耕后地表平整度，也不利于后续播种工作[45]，甚至导致芯土层翻落至地表，这违背了松动土层而非搅浑土层的农艺要求[61]。另外，表3、4也表明α和D的组合影响到纵向碎土距离R的大小，已有研究为深松铲的纵向碎土距离R提供了理论参考[42-43]，其中最易于结合本文的理论分析是郭志军等[44]提出的耕深范围内铲的纵向水平长度L(图7b)。虽然增加L可以有效提升铲对前方土壤的扰动能力，但是这又必然对深松铲的刚度和强度提出更高的要求，因此结合土层物理特点及农艺要求，合理确定深松铲设计参数α和作业参数D是深松铲改善耕作质量的关键，为此仍需要大量的田间原位试验研究。

2.4.2深松铲-土壤系统受力理论检验

深松铲的3个功能区包括铲尖、曲柄和直柄，3个功能区的受力分析受到了普遍关注[28,62-63]。

图8显示扰动土体主要受铲尖上、下表面的正压力Nu、Nd和摩擦力fu、fd作用，各力的水平分力分别为Nux、fux、Ndx、fdx，竖直分力分别为Nuy、fuy、Ndy、fdy，曲柄处土壤主要受到垂直于深松铲刃口的切削力Fn和摩擦力fn，直柄处的土壤主要由铲的挤压力P剪切破碎；由图8受力分析可知，曲柄和铲尖上表面处土壤会沿S3方向向前运动、同时沿S2方向向上运动形成扰动边界，铲尖下表面土壤则沿S5方向向下运动形成新的土壤压实带，此外，随着深松铲的前进，切削力Fn则“挤开”周围土壤形成深松沟槽，其中部分土壤由于铲和周围土体的挤压，沿S1和S4方向向后流动；而其中的耕作阻力F对于机具的优化设计和动力配套最为重要[38]。由表3、4可知，α和D都对耕作阻力F影响显著，F随α线性增加，这与已有的研究结果一致[41]。

图8 深松铲-土壤系统受力分析Fig.8 Force analysis for subsoiler-soil system

对扰动土体进行水平方向上受力分析得

F=P+Fnx+fnx+Nux+fux+Ndx+fdx

(3)

其中

Fnx=Fnsinε

fnx=fncosε=μFncosε
Nux=Nusinα
fux=fucosα=μNucosα
Ndx=Ndsinβ
fdx=fdcosβ=μNdcosβ
β=α-ψ

整理得

F=P+Fn(sinε+μcosε)+Nu(sinα+μcosα)+
Nd(sin(α-ψ)+μcos(α-ψ))

(4)

式中μ——深松铲与土壤之间的动摩擦因数

ψ——深松铲铲尖刃角

β——深松铲铲尖隙角

查阅文献[49]可知，μ≤0.35，故F随α增加。另外试验发现，入土角为15°时，深松铲的入土性能较差，因此在耕作阻力可接受的范围内，建议选择入土角为22.5°的深松铲。

另外F随D近似指数增加，这与前人研究成果一致[30,38,64]，且耕深大于20 cm时，耕作阻力急剧增加，考虑到不同耕深土壤有效深松区域分布，20 cm耕深最佳[36,39]，若考虑深松后最小化微地貌破坏和利于小麦根系下扎等农艺要求[3,65-67]，30 cm耕深最佳。

3 结论

(1)田间原位摄像法能够提供深松过程的土壤扰动行为参数，进而有利于定量深松铲设计参数和作业参数对耕后地表平整性和土壤紧实度的影响。

(2)田间原位摄像法、土壤微地貌三维重构、耕作阻力测试、EDEM仿真、深松扰动剖面土壤紧实度测试等方法分别提供了深松过程及深松前后的多方面作业参数，今后针对土壤耕作力学的原位精准研究需要使用多种测试方法。

(3)针对稻麦轮作区水稻土耕层浅薄的问题，采用22.5°入土角并进行20 cm深松作业能够在节约能耗的情况下实现最大土壤扰动。

(4)从耕后地表平整性方面考虑，30 cm耕深最优，同时可以提供更深的沟槽孔隙增强田间灌排能力。