秸秆旋埋还田后空间分布效果仿真与试验

周 华 李 栋 刘政源 李支轶 骆双成 夏俊芳

(1.华中农业大学工学院， 武汉 430070； 2.农业农村部长江中下游农业装备重点实验室， 武汉 430070)

0 引言

秸秆还田是保护性耕作的一种方式，秸秆还田作业一方面可解决秸秆焚烧带来的环境污染问题[1-3]，另一方面秸秆在土壤中腐解可为作物提供丰富的氮、磷、钾及有机物等，可有效减少化肥的使用量，提高作物产量[4-9]。

秸秆还田作业从秸秆与土壤的位置关系可分为覆盖还田(秸秆在土壤表面)[10-12]、开沟深埋还田(秸秆埋在土壤下面)[13-15]及秸秆土壤混埋还田(秸秆与土壤相对均匀混合)[16-18]等。不同的秸秆还田方法各有优缺点。秸秆覆盖还田是将秸秆覆盖在土壤表面，可有效降低土壤表面水分蒸发，对土壤的蓄水保墒具有一定效果，但秸秆都留在地表，不但秸秆腐解速度较慢，还对后茬作物的播种作业产生不利影响。开沟深埋还田是通过特定机具开出一条深沟，然后将秸秆聚集起来埋入沟底，再回填土壤覆盖，作业过程较为复杂，虽然可以将大部分秸秆埋入沟底，但秸秆在沟底聚集成堆，也严重影响秸秆的腐解速度。秸秆土壤混埋还田是对秸秆和土壤进行一定的扰动，达到秸秆和土壤混合的目的，秸秆和土壤混合得越均匀则秸秆与土壤的接触面积越大，不但有利于提高土壤中微生物群的代谢活动，还有助于提高秸秆的腐解速度[19-20]。

针对秸秆土壤混埋还田(简称秸秆还田)的研究，前期主要集中在秸秆埋覆率[21]，即仅以秸秆混埋入土壤的百分比作为评价指标，而没有深入研究秸秆埋入土壤中的真实状态。离散元法作为数值模拟方法已广泛应用于农业机械作业过程的研究，国内外学者利用离散元软件对不同的耕作装备进行了仿真和试验，基本证实了离散元仿真可较为准确地模拟实际田间作业[22-26]。应用离散元软件进行秸秆还田空间分布研究可为实际田间秸秆还田效果的研究提供理论支撑。

为了实现秸秆在土壤中空间分布的量化及可视化，本文利用设计的秸秆三维坐标测量装置对3种耕作装备的秸秆还田效果进行测量，并与离散元仿真进行对比分析，以期为耕作装备的秸秆还田效果研究提供参考。

1 材料与方法

1.1 试验地概况

2018年11月在华中农业大学现代农业科技试验基地进行了秸秆还田现场试验。试验现场为秋季稻茬田，秸秆残茬立于田中。土壤类型为粉质黏壤土，黏粒、粉粒和砂粒质量分数分别为31%、63%和6%。在试验前，对土壤和秸秆残留物特征进行测量。所有初始土壤和秸秆残留物参数均采用五点取样法在整个试验田进行测量。土壤坚实度由TJSD-750型土壤紧实度测定仪(浙江托普云农科技股份有限公司，精度：±0.05%)测量。利用环刀采集土壤样品，将土壤样品称量，在干燥箱105℃下干燥24 h，并再次称量以确定土壤含水率和干密度。秸秆残茬长度通过钢直尺直接在田中测量。将1 m×1 m的钢制取样框放置在田中，通过剪刀剪取取样框内所有秸秆残留物并称量以确定秸秆残留量。耕前土壤和秸秆参数如表1所示。

表1 耕前土壤和秸秆参数Tab.1 Soil and straw parameters before tillage

1.2 试验方法

1.2.1耕作装备与试验设计

耕作装备包括3种类型：传统旋耕机(TR)、秸秆旋埋还田机(SR)及深松+秸秆旋埋还田机(SSR)。耕作装备均通过三点悬挂连接在拖拉机上。SR和SSR的主要工作部件旋埋刀辊的结构如图1所示。与传统旋耕机相比其增加了螺旋横刀、弯刀和刀盘。旋耕刀的回转半径为245 mm，螺旋横刀的回转半径略小于旋耕刀，为210 mm。螺旋横刀的主要作用为通过滑切作用将秸秆切断并压埋入土。SSR与SR相比在机架前方加装了4个深松铲，形成了深松+秸秆还田联合耕整方式。该耕作方式通过深松使土壤疏松，有助于后续的秸秆入土。

图1 SR和SSR旋埋刀辊结构图Fig.1 Structure diagram of rotary burying blade roller of SR and SSR1.螺旋横刀 2.旋耕刀 3.弯刀 4.刀盘 5.刀轴

由于拖拉机的前进速度、动力输出轴的转速及耕深等均会对秸秆的埋覆效果产生影响。因此，对3种耕作装备设定相同的作业参数，作业幅宽为2 000 mm，拖拉机作业速度为1.55 km/h，刀辊转速为327 r/min，秸秆还田深度为150 mm，松土深度为250 mm。

1.2.2秸秆空间分布测量工具及取样方法

为实现秸秆空间分布的测量，设计了一种秸秆三维坐标测量装置，如图2所示。该装置由铝型材搭建，其主要测量部分包括1个数显角度尺、1个由游标卡尺改装的数显深度尺和1个数显标尺。

图2 秸秆三维坐标测量装置Fig.2 Measuring device of straw three-dimensional coordinate1.数显深度尺 2.数显角度尺 3.数显标尺 4.取样框

对3种耕作方式作业后的土壤进行取样。取样方式为将一个500 mm×500 mm×150 mm的钢制取样框(无顶面和底面)沿耕作前进方向缓慢压入耕后土壤，压入过程中需借助剪刀将位于取样框边框上的秸秆剪断，以免边框对秸秆位置产生影响。待取样框顶部与地表齐平，将取样框外围土壤挖开，并在前进方向挖开适当长度，然后将一个面积略大于取样框底面的钢制托盘插入取样框底部，插入过程借助铁锤敲入，直至托盘完全将取样框托起。此时，一个土壤样品完全取出。待所有土壤取样完毕，将取样土壤带回实验室进行下一步的测量工作。

1.2.3秸秆分布的测量方法

由于水稻秸秆量较大，且有的秸秆主茎秆和外皮在耕作过程中已被分离，又外皮较杂乱，难以测量，因此只以主茎秆作为测量对象。主茎秆含水率(68.15%)较高，不易发生弯折，基本可视为直线段，秸秆三维坐标测量装置可通过测量主茎秆两个端点的空间坐标来确定整根秸秆的空间位置。针对一些发生弯折的特殊主茎秆，除测量其两个端点空间坐标外，再对其弯折的最大曲率点处测量一个空间坐标，将其近似视为两条直线段，由此可基本确定发生弯折秸秆的空间位置。

测量秸秆坐标时先通过调整角度尺的角度和数显标尺的滑动将深度尺调整到秸秆一端点附近，向下滑动深度尺的探针，同时对角度尺和数显标尺进行微调，直至深度尺的探针刚好落在该端点位置，此时，分别记录数显标尺、角度尺和深度尺的读数L1、θ1和H1。以相同的方法测量另一端点的L2、θ2和H2。由此，该秸秆的空间位置完全确定。

测量过程中先依次测量露出土壤部分的秸秆端点坐标，然后在尽量不扰动秸秆位置的前提下慢慢移出表层土壤，继续测量另一端坐标，直至土样中所有秸秆测量完毕。

该测量装置初始测量值为极坐标，通过换算转换为笛卡尔坐标。将所获得的所有秸秆坐标导入三维绘图软件Pro/Engineer中生成秸秆的三维立体图(以线段的形式表示秸秆)。土壤中秸秆分布状态在绘图软件中可视化显示，如图3a所示。通过在三维图中对立方体进行分层处理(图3b)，以100 mm为单位进行横向划分(沿耕作装备幅宽方向，图3c)和纵向划分(沿耕作装备前进方向，图3d)等，并借助三维软件对各处理区域内秸秆长度进行数值统计，可较为准确地判断秸秆在整个立体空间中分布均匀性。

2 离散元仿真分析

运用离散元仿真软件EDEM 2018建立耕作装备-土壤-秸秆间作用模型，在设定工况下对秸秆还田田间作业状态进行模拟，分析耕作装备对秸秆还田效果的影响。

2.1 耕作模型建立

利用三维绘图软件Pro/Engineer分别对3种耕作装备TR、SR和SSR进行建模，为简化整体结构和运算，对机架等非触土部件进行简化，以.igs格式文件分别将其三维模型导入离散元软件中，仿真模型如图4所示。

图4 SSR离散元仿真模型Fig.4 Discrete element simulation model of SSR

2.2 土壤-秸秆颗粒模型

采用EDEM 2018进行仿真，仿真关键是土壤接触模型和基本仿真参数的确定。选用球形颗粒模拟土壤颗粒，由于离散元仿真中颗粒越小，仿真计算速度越慢，占用计算机存储空间越大，因此仿真颗粒一般都远大于实际土壤颗粒[27-29]，本研究选用土壤颗粒半径为8 mm。为防止仿真过程中虚拟土槽侧面及底面对耕作装备作业的影响，建立长×宽×高为2 500 mm×2 500 mm×400 mm的虚拟土槽。

接触模型是离散元法的重要基础，其实质是准静态下颗粒固体的接触力学弹塑性分析结果[30]。接触模型的分析计算直接决定了颗粒受力和力矩的大小，不同的仿真对象须建立不同的接触模型，以提高仿真结果的准确性。本区域水稻土为黏壤土，因此，设定土壤颗粒间的接触模型为Hertz-Mindlin with JKR。JKR模型引入了颗粒间表面能的概念，考虑了湿颗粒间黏结力对颗粒运动规律的影响，是一个凝聚力接触模型，适用于模拟颗粒间因静电、水分等原因发生明显黏结和团聚的物料，如农作物和土壤颗粒间的黏聚作用[31]。秸秆由于其特殊的性质，在离散元仿真中一直难以模拟，试验中水稻秸秆易弯折性和空心结构等离散元都难以做到，考虑到仿真量较大，为了简化秸秆及减少仿真时间，秸秆颗粒采用刚性模型，以16个直径为10 mm、球心间距为6 mm，总长度为100 mm的线性模型来模拟秸秆。

在EDEM软件中，模型参数主要包括材料参数和接触参数。材料参数主要通过文献和试验获得，模型参数主要通过文献、试验和参数标定获得。选取土壤、秸秆和钢的泊松比分别为0.38、0.4和0.3，密度分别为2 680、227、7 865 kg/m3，剪切模量分别为1×108、1×106、7.9×1010GPa[32-33]。通过文献、试验测量及参数标定确定模型的接触参数如表2所示[32，34-35]。

表2 模型的接触参数Tab.2 Contact parameters of model

2.3 仿真过程

为减少仿真过程中土壤和秸秆颗粒随机生成时位置不同对耕作装备作业效果的影响，仿真时先生成土壤和秸秆颗粒，再替换不同的耕作装备进行仿真作业。仿真过程中，设置2层Polygon虚拟平面分别生成土壤和秸秆颗粒，在前3.4 s内先累计生成661 515个土壤颗粒，再生成2 500个秸秆颗粒。为便于区分耕层土壤和底层土壤，0～150 mm和150～400 mm土壤颗粒以不同的颜色区分。由于秸秆难以实际模拟田间秸秆状况，因此将秸秆颗粒均匀分布在土壤表面。耕作装备作业过程中设定与田间试验相同的作业参数，设置Rayleigh时间步长为10%(1.41×10-5s)，网格尺寸设置为最小颗粒半径的2倍，仿真总时间为12 s，其中0～3.4 s为颗粒生成时间，3.4～10 s为作业时间，10～12 s为耕作装备离开和颗粒稳定时间。

在仿真作业完成后，通过设置计算区域Geometry Bin来计算区域内秸秆数量。对每种耕作装备作业后随机选取3个位置(以3个位置的平均值作为结果)，每个位置设置3层计算区域，分别对应实际田间作业的垂直分层处理，每层计算区域尺寸为500 mm×500 mm×50 mm。隐藏土壤颗粒，保留秸秆颗粒，效果如图5所示。图中处于上、中、下层计算区域的秸秆显示为橙色，黑色为计算区域外的秸秆，若秸秆颗粒的长度超过一半处于计算区域，则该秸秆算作区域内秸秆，否则不在计算范围内。同理，通过设置不同形式和尺寸的计算区域分别对应实际田间作业的横向和纵向划分。

图5 3层计算区域内秸秆模型Fig.5 Straw models in three-layer computing areas

3 结果与分析

3.1 秸秆在土壤中的长度

本文主要研究秸秆在土壤中的分布状况，秸秆计算中仅以埋入土壤中的部分作为研究对象。通过三维软件的测量功能统计，3种耕作装备田间试验后的秸秆参数如表3所示。数据显示SR作业后埋入土壤中的秸秆总长为9 708.5 mm，略大于SSR的8 973.6 mm，二者相差不多，但都远远大于TR作业后的值，约为其2倍，因此SR和SSR的秸秆埋入效果明显优于TR。从图6可以看出，TR的秸秆埋入效果明显不如SR和SSR。作业后秸秆的平均长度表明TR的秸秆切断效果最差，SR的效果次之，SSR的秸秆切断效果最好。其主要原因是SR和SSR增加的螺旋横刀相应增加了秸秆被切割的概率。

表3 3种耕作装备作业后的秸秆参数Tab.3 Straw parameters of three tillage equipments after operation mm

3.2 秸秆在土壤中的垂直分布

分层处理主要用于评价秸秆还田后秸秆在土壤中的垂直分布，秸秆在各层的占比越均匀，则秸秆还田后秸秆在土壤中的垂直分布越理想。3种耕作装备作业后秸秆垂直分布的离散元仿真值与田间试验实测值如表4所示。

图6 TR、SR和SSR的秸秆埋入状况Fig.6 Straw incorporation status of TR, SR and SSR

为便于计算，仿真值以秸秆数量作为测量指标，实测值以秸秆长度作为测量指标。由表4可知，仿真值中SR和SSR埋入土壤中的秸秆数量明显大于TR，尤其是下层秸秆数量更加明显。实测值中埋入3个层次的秸秆长度也均大于TR，SR和SSR埋入下层秸秆长度分别达到了TR的5.4倍和6.2倍。3种耕作装备作业后各层秸秆占比的仿真值与实测值比较接近。其中，TR作业后，仿真值和实测值中分别有88.1%和92.8%的秸秆被埋在上层和中层，上层和中层几乎各占一半，分别仅有11.9%和7.1%的秸秆被埋在下层。仿真值与实测值误差较小，由表4也可以看出，TR主要将秸秆埋覆在上层和中层，很难将秸秆埋入下层。SR作业后，仿真值和实测值中秸秆大部分分布在上层和中层，但所占比例明显低于TR，分别有16.1%和18.1%的秸秆被埋在下层。秸秆在上层和中层之间分布不均匀，仿真值和实测值中分别有46.8%和47.4%的秸秆处于中层，这表明SR作业后有接近一半的秸秆被埋入土壤的中层位置。虽然SR作业后秸秆在各层的比例优于TR，但中层位置秸秆较为集中。SSR作业后，仿真值和实测值中秸秆在上层和中层的占比进一步降低，下层的秸秆占比则相应增加，达到了22.5%和22.3%，虽仍未达到3层秸秆占比的平均值，但已明显优于前两者。

表4 3种耕作装备下秸秆垂直分布的仿真值与实测值Tab.4 Simulated and measured values of straw vertical distribution under three tillage equipments

综合来看，3种耕作装备作业后秸秆在土壤中的分布具有相似的规律，即大部分分布在土壤的中层和上层。变异系数可用来比较数据间离散程度，其定义为数据标准差和平均值的比。各层秸秆占比变异系数如图7所示。

图7 3种耕作装备下各层秸秆占比的变异系数Fig.7 Coefficient of variation of straw proportion in different layers under three tillage equipments

由图7可知，TR、SR和SSR作业后各层秸秆占比仿真值和实测值的变异系数均呈递减趋势，其中TR仿真值和实测值的变异系数最大，分别为55.8%和68.2%，SSR的变异系数最小，分别为28.8%和28.7%。这说明不论是离散元仿真还是田间试验，SSR秸秆还田后，秸秆在土壤中垂直分布的均匀性均最优。TR、SR和SSR的秸秆在土壤中垂直分布均匀性依次增大的主要原因是由于TR的主要工作部件是旋耕刀，旋耕刀滑切角较大，滑切刃较短，旋耕刀旋转作业过程中，秸秆较容易从滑切刃滑出或被切断，被旋耕刀埋入土壤下层的概率较低，因此，垂直分布均匀性最低。SR和SSR的主要工作部件除旋耕刀外还有螺旋横刀，螺旋横刀的滑切角小于旋耕刀，滑切刃长度大于旋耕刀，作业过程中秸秆并不容易滑出，螺旋横刀针对秸秆一方面起到切碎作用，另一方面起到压埋作用。SR作业条件下，螺旋横刀直接作用于未深松土壤，土壤坚实度较大，土壤在螺旋横刀作业过程中对秸秆起到一定的支撑作用，秸秆偏向于切碎，向下的压埋效果略低。而SSR作业条件下，先深松作业，深松使一定范围内的土壤产生松动，其土壤坚实度大大降低，疏松后的土壤给秸秆提供的支撑力小于前者，螺旋横刀压埋秸秆的入土阻力较小，可以将秸秆压入更深的土层，达到一定支撑力后再切断，因此，SSR的垂直分布均匀性最优。

3种耕作装备下仿真值与实测值的变异系数相差不大，TR、SR和SSR的相对误差分别为22.2%、6.4%和0.3%，其中，TR的实测值变异系数略大，SR和SSR的仿真值变异系数略大。3种耕作装备的仿真值与实测值变异系数的平均误差为9.6%。该误差主要是由于离散元仿真只是使各种作业参数尽量接近真实值，但无法达到所有参数与实际一样，离散元仿真参数的变化会影响最终结果。以上对比分析表明，3种耕作装备秸秆还田后的各层占比变异系数可以量化表达离散元仿真和田间试验下秸秆在土壤中垂直分布的均匀性，且仿真值和实测值的相对误差较小，变化趋势基本一致。

3.3 秸秆在土壤中的水平分布

横向和纵向划分主要用于评价秸秆还田后秸秆在土壤中的水平分布。与秸秆在土壤中的垂直分布类似，以横向和纵向划分后各区域秸秆占比的变异系数表示秸秆在土壤中的水平分布，如图8所示。

图8 3种耕作装备下秸秆水平分布的变异系数Fig.8 Coefficient of variation of straw horizontal distribution under three tillage equipments

由图8可知，在对取样立方体进行横向划分和纵向划分后，TR、SR和SSR的各区域秸秆占比仿真值和实测值的变异系数无绝对规律，但整体而言，SSR的变异系数均最小，说明SSR秸秆还田后，秸秆在土壤中水平分布的均匀性最优。在横向划分中，TR作业后的仿真值和实测值的变异系数间相对误差不大，仅为4.0%，但SR和SSR的相对误差较大，分别为30.1%和42.0%，三者间的平均相对误差为25.4%。在纵向划分中，SR作业后的仿真值和实测值的变异系数间相对误差不大，仅为2.2%，但TR和SSR的相对误差较大，分别为47.1%和39.0%，三者间的平均相对误差为29.4%。

在对取样立方体分层处理下，3种耕作装备作业后各层秸秆占比的仿真值与实测值变异系数的平均误差仅为9.6%，而横向和纵向划分处理下，该值分别达到了25.4%和29.4%。其原因为：3种耕作装备下，无论是仿真作业还是实际田间作业后，秸秆在土壤中的空间状态都有相似的特性，即与水平面呈一个较小的角度，如TR、SR和SSR的田间作业实测值经计算后该角度的平均值分别为29.4°、19.9°和20.1°，仿真作业后该角度同样较小，甚至一些秸秆与水平面近乎平行，如图9a所示。秸秆与水平面间较小的角度令秸秆在垂直方向上的投影长度较小，在仿真作业的分层处理中，秸秆更容易划分到单一的层次中，因此，与实测值误差不大。横向和纵向划分是对水平分布的处理，秸秆与水平面间较小的角度令秸秆在水平方向上的投影长度较大。以纵向划分为例(图9b)，A和B为两个相邻的纵向区域，秸秆L的长度约有60%处于A区域，约有40%处于B区域，而仿真计算中的秸秆占比以秸秆的数量作为计算结果，因此，计算过程中秸秆L属于A区域而不属于B区域。但在田间试验结果处理时，秸秆占比以秸秆的长度作为计算结果，借助三维软件可以计算秸秆在不同区域的长度，在计算图9b中的情况时，秸秆L会分成两段，分别计算A区域和B区域的长度，所得结果也更为精确。而横向划分也是类似，因此，仿真值和实测值的误差大于分层处理。

图9 仿真后秸秆位置空间状态Fig.9 Spatial state of straw after simulation

3.4 讨论

为便于统计作业后土壤中秸秆数量并节约仿真时间，本研究利用离散元法建立了刚性秸秆模型，模拟不同耕作装备的秸秆还田作业。建立刚性秸秆模型是目前土壤-秸秆-机具联合仿真中较常利用的方法，该方法可较为清晰地观测秸秆运动，定位秸秆位置。目前，已有学者开始利用离散元法建立柔性秸秆颗粒[36-37]，虽然柔性秸秆颗粒物理特性更接近于实际秸秆，但鉴于建立柔性秸秆较为困难并且柔性秸秆仿真后难以测量秸秆数据，秸秆弯曲或断裂后无法计算秸秆数量，更难以分析秸秆的垂直分布和水平分布，因此，本研究更适合使用刚性秸秆模型。土壤颗粒间的接触选用了JKR模型，仅针对本地区土壤而言，不同含水率下，土壤的接触参数会有所不同。因此，不同土壤含水率条件下进行仿真作业时需重新标定土壤参数。

通过实际测量秸秆空间坐标的方法研究旋埋还田作业后的秸秆空间分布效果。该研究开展的关键是秸秆坐标的测量，田间试验的研究对象是含水率相对较高的水稻秸秆，该秸秆被旋埋后基本能通过测量2、3个端点确定其空间位置。该测量方法有其优点也有其局限性，若是油菜等硬质且不易弯曲的秸秆，则更适合该测量方法，若是含水率较低且易弯曲的秸秆或已被粉碎机粉碎过的较杂乱的秸秆，则无法利用此测量方法来判断秸秆旋埋还田的空间分布效果。

4 结论

(1)利用传统旋耕机(TR)、秸秆旋埋还田机(SR)及深松+秸秆旋埋还田机(SSR)测量了秸秆还田后秸秆空间分布效果。设计了秸秆三维坐标测量装置，通过该装置可测量秸秆还田后秸秆在土壤中的空间坐标，将所测得的空间坐标导入三维绘图软件可还原秸秆在土壤中的空间分布，及秸秆在三维图中量化与可视化显示。

(2)利用离散元法建立了相应的仿真模型，并与田间试验设定了相同的作业参数。在仿真作业完成后，通过设置不同形式和尺寸的计算区域Geometry Bin计算区域内秸秆数量，并分别对应实际田间作业的分层、横向和纵向划分。

(3)在分层处理中，仿真与实测结果表明，SR和SSR埋入土壤中的秸秆量都明显大于TR，尤其埋入土壤下层的秸秆量均是TR的数倍。TR、SR和SSR作业后各层秸秆占比仿真值和实测值的变异系数均呈递减趋势，其中TR仿真值和实测值的变异系数最大，分别为55.8%和68.2%，SSR的变异系数最小，分别为28.8%和28.7%。离散元仿真与田间试验结果表明，SSR秸秆还田后，秸秆在土壤中垂直分布的均匀性均最优。3种耕作装备下仿真值与实测值的变异系数相差不大，平均误差为9.6%。在对取样立方体进行横向和纵向划分后，TR、SR和SSR的各区域秸秆占比仿真值和实测值的变异系数无绝对规律，SSR的变异系数均最小，SSR秸秆还田后秸秆在土壤中水平分布的均匀性最优。虽然仿真值和实测值的水平分布误差超过20%，但通过优化计算方式可减小该误差。离散元仿真较好地拟合了实际田间作业后秸秆的空间分布状态，相对误差在可接受范围内。