基于响应面法的蓝莓离散元参数标定

张刘炜 ，王海滨 ，郭艳玲 ，李 森 ，王枫辉

（东北林业大学 a.工程技术学院，b.机电工程学院，哈尔滨 150040）

蓝莓与沙棘、山葡萄和树梅一同被联合国粮农组织(FAO)确定为人类21 世纪健康食品，被誉为第三代水果中的黄金水果和浆果之王[1-3]。 蓝莓含有丰富的花青素，在保护视力，增强免疫力上有很大的功效，国际国内消费市场越来越好。 因此，国内外对蓝莓采摘装备的研究也日益增加[4-5]。 针对现有的龙门式蓝莓采摘机，在采摘的过程中发现，蓝莓果实运动轨迹有一定的规律，采摘机的接果板尺寸、机体宽度等会对蓝莓的运动轨迹产生影响。近些年来，作为新的分析算法，离散元法（discrete element method，DHM）在农业装备研究上有着越来越多的应用[6-11]。

进行蓝莓的离散元仿真分析，需要确定蓝莓果实离散元仿真的材料参数和接触参数，可靠的参数是离散元结果准确的保证，有助于离散元法在蓝莓采摘装备研发中的应用。 查阅国内外的文献，目前有两种方法得到这些参数，直接测量和参数标定法。直接测量中，有COHTZHH 等[12]采用剪切试验和压缩试验测出玉米的内摩擦角和刚度，BALHVIĈIUS 等[13]使用豌豆和有机玻璃进行了滑动试验，测量了两者之间的静摩擦系数，鲍玉冬[14]通过静准态压缩等试验和间接计算的方式，估算出了蓝莓的离散元参数，用于蓝莓采摘机的收获过程仿真。 由于农业物料的生物学特性，形状不规则，成熟度的差异也会带来参数的差异，故农业物料的离散元参数更为复杂，颗粒与颗粒间的参数值相差较大，直接测量的参数值不总是可行，代表性差。 很多的学者选择了使用参数标定法来获取农业物料的离散元参数，主要步骤是通过不断地改变参数值，直到模拟结果与相应的实体试验之间取得最佳匹配。 王云霞等[15]进行了休止角的试验，并生成数学模型，以此确定了玉米种子间静摩擦系数和滚动摩擦系数；刘凡一等[16]基于响应曲面法，应用Plackett-Burman 和Box Behnken 试验，标定了小麦的离散元模型参数；冯俊小等[17]进行了堆积仿真试验、滚筒内颗粒的仿真试验和实体试验，通过对比，标定了试验颗粒的离散元参数；武涛等[18]进行了Box Behnken 试验和堆积角仿真试验，进而标定了土壤离散元模型的参数。

本研究以蓝莓果实为研究对象，以休止角为响应值，使用了响应面法，建立了回归模型，对蓝莓的离散元参数进行标定，且分析了参数间的交互作用。 使用标定后的值进行仿真，与实体试验进行对比发现无显著性差异。 故能将标定的参数用于蓝莓采摘机工作过程中的蓝莓运动轨迹分析，优化采摘机的结构参数。

1 材料与方法

1.1 蓝莓基本参数

在辽宁省丹东市振安区五龙背镇新康蓝莓种植基地进行蓝莓采摘机的田间收获作业，收获的过程和采收的蓝莓如图1，使用该批蓝莓进行试验。 采摘的蓝莓品种为蓝丰，根据目测，单颗的蓝莓果实为扁圆形。 由于蓝莓果实的生物学特性，果实的大小具有随机性，符合正态分布，故随机选取采摘的部分蓝莓，首先依据目测按大小分成大、中、小3 组。 对各组蓝莓进行计数，并在每组蓝莓中随机挑选5 颗蓝莓，使用游标卡尺测量每组蓝莓的直径（d）和高(h)（图2），每个数据测量5 次后取均值，记在表1 中。 通过上海恒平电子天平JA5003 测量单颗蓝莓的平均质量在 0．0020～0．0035kg，通过排水法测量体积，间接测量出它的密度为 1．018×103kg·m-3。

图1 田间作业Figure 1 Field operation

1.2 标定方法

休止角是表征散粒体物料颗粒流动、摩擦等特性的参数[19-20]，休止角试验常作为农业物料颗粒离散元模型的验证试验[21-22]。 用休止角的值作为表征蓝莓果实离散元参数的目标响应值，通过实体的圆筒提升试验与离散元仿真相结合的方式， 对蓝莓的离散元模型参数进行标定。 实体试验中，使用内直径为16mm 的无底圆筒，放在边为60mm×60mm 的钢板中心，将蓝莓随机落入圆筒中，高度约250mm，上表面铺平。 用电机带动圆筒以0．05m·s-1的速度匀速上升，由于重力的作用，会逐渐形成蓝莓锥形堆（图3），待蓝莓锥形堆稳定后使用直尺测量出蓝莓堆的高(h)和底面直径(d)，可以计算出休止角。 共进行5组重复试验， 测量的休止角平均值为 24．52°， 作为后续仿真试验的参考角度。

在三维建模软件Pro/H 中建立试验的模型，导入软件HDHM 2018 中， 然后建立蓝莓仿真模型。 通过Design Hxpert 8 软件进行Plackett-Burman 试验设计，筛选出显著性的参数，再通过最陡爬坡试验，确定中心点，作为响应面法的中心，然后进行中心组合试验，寻找合适的回归模型，确定蓝莓离散元模型的最优值。

1.3 仿真模型

1.3.1 仿真的接触模型 目前常见的颗粒模型有软球模型和硬球模型两种。硬球模型，不考虑颗粒接触的过程，假设颗粒间的接触和碰撞都是瞬时完成，颗粒接触后直接得到颗粒速度，通过接触力对时间的积分来描述接触的整个过程。 软球模型，会引入弹簧、阻尼和滑动器，不用考虑接触力作用的整个过程，直接计算两个颗粒间的变形重叠部分。 采摘的蓝莓果实硬度低，故碰撞的接触过程不能看作瞬时的，适合采用软球模型。在软球模型中，会将蓝莓间的接触过程简化为弹簧振子的阻尼振动(图4)，运动方程为：

式中：m 为振子质量(kg)；x 为振子距离平衡位置的位移量（m）；c 为弹黃阻尼系数(Ns·m-1)；k 为弹性系数(N·m-1)。

Hertz-Mindlin（no slip）模型是目前较为常用的软球模型。所以本研究的离散元仿真应用Hertz-Mindlin（no slip）模型表征颗粒与钢板、颗粒与颗粒间的接触碰撞。

图3 蓝莓堆Figure 3 Pile of blueberry

图2 蓝莓果实Figure 2 Blueberriey fruit

表1 蓝莓的尺寸Table 1 The size of blueberries

图4 弹簧阻尼系统Figure 4 Spring damping system

1.3.2 蓝莓果实模型 蓝莓果实为扁球形， 在软件HDHM 2018 中， 使用四球模型来作为蓝莓的仿真模型 （图5），通过依次指定每个圆球模型的相对位置和直径来生成模型。 根据前面测量的蓝莓尺寸数据，蓝莓的仿真模型也分为 3 组，标准直径依次为 20．61，17．48，15．32mm，对应的标准高为 12．83，11．38，10．70mm。 仿真时，蓝莓依照标准尺寸按照正态分布生成 （均值为1， 标准差为0．05），每组蓝莓的数目比为 1∶2∶1，与实测数据保持一致。

本文以晶体管直流增益作为中子辐射损伤效应的宏观表征参数，采用电压补偿方法解决了远程监测中因晶体管工作电压损耗引起的测量误差，并通过模块化软件架构以及电压回读技术，建立了晶体管直流增益在线测试系统，实现了不同中子注量辐照下晶体管直流增益的实时监测，获得了辐照期间晶体管直流增益随不同中子注量的变化规律，为晶体管的中子辐射损伤效应评估提供了重要的测试依据。

图5 蓝莓模型Figure 5 Blueberry model

1.3.3 圆筒提升模型 离散元仿真中，圆筒的内径和高度、底部钢板的尺寸都与实物试验中采用的一致，在Pro/H 中完成建模，导入到HDHM 软件中。设定蓝莓颗粒在圆筒顶面的面型粒子工厂中自动随机生成，生成速率为1000 个·s-1，总生成时间设置为2s，生成后的颗粒在重力作用下做自由落体运动，填充圆筒。待系统稳定后，起初在软件中给圆筒添加的0．05m·s-1向上的线性速度会出现，缓慢的，圆筒会与底部钢板分离，由于重力作用，蓝莓果实模型将从圆筒底端缓慢流出，圆筒完全与蓝莓颗粒分离后，蓝莓会在钢板上形成锥形堆。

将Rayleigh 时间步长设为15%，数据保存时间间隔为0．05s，网格大小为2．5 倍最小模型颗粒直径，进行仿真，仿真总时长设定为8s。

1.3.4 参数设置 在软件HDHM2018 中，为了区别不同的物料，仿真前需要设定一些参数。 需要设置的参数包括物料颗粒（蓝莓）和几何体(圆筒、底部钢板)的材料参数和接触参数。 实体试验中，几何体为不锈钢，密度、泊松比和剪切模型参数确定， 而蓝莓果实由于生物学特性，离散元参数的值并不唯一。 蓝莓果实的密度可以通过简单的实验测出，但通过实验测量其他参数的过程特别繁琐，通过查阅国内外文献[14，23-24]，得到蓝莓泊松比、蓝莓剪切模量、蓝莓—蓝莓恢复系数、蓝莓—蓝莓静摩擦系数、蓝莓—蓝莓滚动摩擦系数、蓝莓—钢恢复系数、蓝莓—钢静摩擦系数和蓝莓—钢滚动摩擦系数的变化范围，见表2。

1.4 试验方法设计

1.4.1 Plackett-Burman(PB) 试验设计仿真中用的参数对休止角的影响大小不同，并非所用的参数都对休止角有显著影响[25-27]，有显著影响的参数才能基于响应面法标定出来，能用于仿真分析。 本研究应用软件Design Hxpert8[28-29]首 先 进 行 Plackett-Burman Design[30]，在 多 个参数中筛选出有显著效应的参数。根据HDHM 仿真时要用到的参数，去除可以测量的蓝莓密度，选取了蓝莓泊松比、蓝莓剪切模量、蓝莓—蓝莓恢复系数、蓝莓—蓝莓静摩擦系数、蓝莓—蓝莓滚动摩擦系数、蓝莓—钢恢复系数、蓝莓—钢静摩擦系数和蓝莓—钢滚动摩擦系数等8 个不便直接测量的参数作为变量，由于有8 个变量，预留3 个虚拟变量作误差分析，本研究中采用N=11的Plackett-Burman 设计表。 以蓝莓果实休止角为响应值，每个参数选取高和低2 个水平，记为+1 和-1（表3）。

1.4.2 最陡爬坡试验 通过Plackett-Burman 试验可以筛选出对休止角影响显著的几个变量，但并不能确定变量的具体值。 因次可以利用最陡爬坡试验，来确定变量最优值所在的区间。 影响不显著的变量在最陡爬坡试验中都取平均水平，根据Plackett-Burman 试验结果中显著因素的正负效应来确定变量的增减，按照选定的步长进行试验，并计算仿真休止角与实际休止角的相对误差。

1.4.3 中心点响应面试验 根据最陡爬坡试验的结果，采用Central Composite Design(CCD)试验进行响应面分析，显著性变量在试验设计中取5 个水平，不显著变量取平均水平，试验设计中根据经验采用6 个中心点，设计的试验共进行20 次。 对结果进行分析，通过不断剔除和调整拟合方程中影响不显著的项，优化回归方程，能得到最佳的参数。

表2 仿真参数设置Table 2 The setting of simulation parameters

表3 Plackett-Burman Design 因素水平Table 3 Factors and levels of Plackett-Burman Design

1.4.4 验证试验 用得到的最佳蓝莓果实仿真参数，在HDHM 2018 中进行仿真试验，将仿真试验的结果与实体试验结果进行对比，计算它们的相对误差，来验证得到的参数是否可靠，能否用于蓝莓的离散元仿真中。

1.5 休止角的测量

在仿真试验中，由于HDHM 软件自身测量数值的局限性，故使用计算机图像处理的方法来计算出休止角[31]，准确方便。 得到休止角仿真试验后的图像，利用Python 语言并结合openCV 库和一些科学计算的第三方库，编写程序，对图像进行二值化处理，提取边界线和最小二乘法拟合，所得拟合直线的斜率即所求休止角的正切值，经过转换可以得到休止角。

图 6 图像处理过程Figure 6 Process of image analysis

2 结果与分析

2.1 Plackett-Burman Design 筛选显著性参数

利用软件Design Hxpert 8 对试验结果进行分析，得到选择的8 个参数对休止角的影响效果和显著性，见表5。 由表 5 可知，蓝莓—钢静摩擦系数（G）、蓝莓—蓝莓静摩擦系数（D）和蓝莓—蓝莓恢复系数（C）对休止角的影响显著，而其他因素对休止角的影响非常小，可以忽略。 因为蓝莓为扁球形，非圆球形，相比于有机肥等圆形颗粒，静摩擦系数对休止角的影响比滚动摩擦系数对休止角的影响显著。

表4 Plackett-Burman Design 方案及结果Table 4 Scheme and results of Plackett-Burman Design

由图7 可知，各参数对休止角影响的显著性排序，并观察到蓝莓—钢静摩擦系数（G）、蓝莓—蓝莓静摩擦系数（D）为正效应，蓝莓—蓝莓恢复系数（C）为负效应，可以确定在后面的最陡爬坡试验中参数的变化方向。

表5 Plackett-Burman Design 试验参数显著性分析Table 5 Analysis of significance of parameters in Plackett-Burman Design test

图7 PB Design 帕累托图Figure 7 Pareto chart of PB Design

2.2 最陡爬坡试验确定显著性参数中心点

最陡爬坡试验可以较快地确定参数的最优值， 即为中心组合的响应面法提供一个中心点。 根据Plackett-Burman 试验的结果，将蓝莓—钢静摩擦系数（G）、蓝莓—蓝莓静摩擦系数（D）按选定的步长增加（正效应），蓝莓—蓝莓恢复系数（C）按选定的步长减少（负效应），其余的参数取中间的水平（蓝莓泊松比0．35，蓝莓剪切模量 0．03MPa，蓝莓—蓝莓滚动系数 0．05，蓝莓—钢恢复系数 0．18，蓝莓—钢滚动系数 0．05）。 试验设计的方案见表6，共进行6 次试验，计算仿真中蓝莓的休止角与实际休止角的相对误差，相对误差最小的参数组即作为响应面的中心点。

由表6 可知，仿真中的蓝莓休止角与实际休止角的相对误差由最初的95．07%下降到1．59%，随后又上升至20．35%，变化的趋势明显。 第5 组试验的相对误差最小为1．59%，说明各参数的最佳值在第5 组附近，即选取蓝莓—蓝莓恢复系数0．1， 蓝莓—蓝莓静摩擦系数0．5，蓝莓—钢静摩擦系数0．65 作为响应面法的中心点。

2.3 中心组合设计试验寻找参数最优值

进行最陡爬坡试验后， 接着用Central Composite Design 响应面[32]分析试验。和最陡爬坡试验一样，非显著性参数依然选择中间水平（蓝莓泊松比0．35，蓝莓剪切模量 0．03MPa，蓝莓—蓝莓滚动系数 0．05，蓝莓—钢恢复系数 0．18， 蓝莓—钢滚动系数 0．05），3 个显著性的参数的设置情况如表7。 根据经验，误差估计用6 个中心点，共进行20 次试验，试验的方案及结果如表8。

方差分析结果如表9，蓝莓—蓝莓恢复系数(C)、蓝莓—蓝莓静摩擦系数(D)、蓝莓—钢静摩擦系数(G)3 个参数对蓝莓休止角的影响十分显著，p＜0．01，该模型的p 值小于0．0001，说明因变量和自变量之间的关系极显著；失拟项p=0．1818＞0．05， 变异系数 CV 为 2%， 比较低， 说明试验有较好的可靠性； 此外方程的决定系数 R2=0．9667，校正决定系数 R2adj=0．9368，都接近 1，说明所得的拟合方程可靠度高，精密度 Adeq Precision=22．334，回归模型有良好的精确度。

通过去除影响不显著的项，如CD、C2和G2等，完成对回归模型的优化。 优化后的回归模型方差分析如表10，发现优化后的失拟项 p=0．2586＞0．05，大于优化前的失拟项，方程拟合度增强；决定系数 R2=0．9610、校正决定系数R2adj=0．9431，比之前有所改善，可靠性高；变异系数CV=1．89%，低于之前的值，可靠性增强；精确度Adeq Precision=27．919，大于之前的精确度，模型更加精准。优化后的回归方程在各项指标上都优于之前，有更好的拟合度和可靠性，优化后的回归方程为：

表6 最陡爬坡试验方案及结果Table 6 Scheme and results of steepest ascent test

表7 Central Composite Design 因素水平Table 7 Factors and levels table of Central Composite Design

R=6．47367+15．95407×C+36．19523×D+2．19483×G-45．66667×C×G+21．00000×D×G-35．91890×D2

同时，蓝莓—蓝莓恢复系数(C)与蓝莓—钢静摩擦系数(G)的交互项(CG)对休止角的影响显著，p＜0．05。 取蓝莓—蓝莓静摩擦系数(D)为0．56，绘制此时的蓝莓—钢静摩擦系数与蓝莓—蓝莓恢复系数的交互作用的响应曲面（图8）。 由图8 可知，在蓝莓—钢静摩擦系数确定时，休止角随着蓝莓—蓝莓恢复系数增大而减小，在蓝莓—钢静摩擦系数越大时， 这变化的速率越大； 在蓝莓—蓝莓恢复系数确定时，休止角随着蓝莓—钢静摩擦系数增大而增大，在蓝莓—蓝莓恢复系数越小时，这变化的速率越大。

最后，根据回归模型的方程，将实物试验的休止角回代到方程中， 通过Design Hxpert 8 得到3 个显著性参数的最优值，即蓝莓—蓝莓恢复系数(C)为0．06，蓝莓—蓝莓静摩擦系数(D)为 0．56，蓝莓—钢静摩擦系数(G)为 0．70。

表8 Central Composite Design 方案及结果Table 8 Scheme and results of Central Composite Design

表9 Central Composite Design 模型方差分析Table 9 ANOVA of Central Composite Design model

2.4 仿真验证

将得到的参数蓝莓—蓝莓恢复系数(C)0．06、蓝莓—蓝莓静摩擦系数(D)0．56、蓝莓—钢静摩擦系数(G)0．70 代入到离散元模型中进行仿真试验，进行5 次试验，测得蓝莓仿真休止角a 的平均值为24．66°（图9），与蓝莓实际休止角24．52°的相对误差为0．57%，无显著性差异，说明标定的参数是可靠的，可以应用到关于蓝莓的离散元模型中，用于优化蓝莓采摘机的结构。

表10 Central Composite Design 优化回归模型方差分析Table 10 ANOVA of Central Composite Design modified model

3 讨论与结论

以往在用离散元法对蓝莓采摘过程进行研究时，蓝莓的离散元参数通常是使用直接测量的方法将所需的参数逐一测量出来，但是由于蓝莓品种和成熟度的不同，直接测量的参数准确度和适用性都较差，而测量的过程也较繁琐，对测量仪器的要求高。 相对于之前的直接测量法，通过试验证明，本研究使用的基于响应曲面法的间接标定法是可行的，能用于蓝莓这种农业物料，标定方式简易。 与相似的研究对比发现，如有机肥这种圆球形的物料，有机肥间的滚动摩擦系数对休止角的影响显著， 而蓝莓间滚动摩擦系数对休止角的影响并不显著，是因为近似椭球形的蓝莓不易滚动，而圆球形的有机肥则更易于滚动，所以其滚动摩擦系数对休止角的影响显著；相对于玉米颗粒而言，玉米颗粒的硬度大且玉米表面相对光滑，而蓝莓的表面比玉米软，表皮也带有一定的黏性， 所以蓝莓间的静摩擦系数对休止角的影响显著，而玉米间的动摩擦系数对休止角的影响显著，对比的结果也与日常生活的经验相吻合，也一定程度上证明了本研究结果的正确性。

本研究结果表明，使用Plackett-Burman 试验筛选出对蓝莓休止角有显著影响的3 个参数：蓝莓—蓝莓恢复系数、蓝莓—蓝莓静摩擦系数、蓝莓—钢静摩擦系数；蓝莓泊松比、蓝莓剪切模量、蓝莓—蓝莓滚动摩擦系数、蓝莓—钢恢复系数和蓝莓—钢滚动摩擦系数对休止角的影响不显著。根据最陡爬坡试验的结果，确定3 个参数的最优值，并以此进行Central Composite Design 响应面分析试验，建立参数和蓝莓休止角之间的回归方程，对回归方程进行优化和方差分析，分析表明3 个显著性参数对蓝莓果实休止角影响显著，此外蓝莓-蓝莓静摩擦系数的二次项和蓝莓-蓝莓恢复系数与蓝莓-钢静摩擦系数的交互项对蓝莓果实休止角也有显著影响。 将实际的蓝莓休止角代入所得的回归方程中，通过Design Hxpert 8 的预测，得到参数的最优值：蓝莓-蓝莓恢复系数为0．06、蓝莓-蓝莓静摩擦系数为0．56、蓝莓-钢静摩擦系数为0．70。 进行的仿真验证试验表明，仿真实验的蓝莓休止角与实际的休止角无显著差异，表明得到的参数是可靠的，能用于蓝莓采摘机的离散元仿真中，分析采摘过程中蓝莓的运动轨迹。

图8 蓝莓-钢静摩擦系数与蓝莓-蓝莓恢复系数的交互作用Figure 8 Interaction of Blueberry-steel static friction coefficient and Blueberry-Blueberry restitution coefficient

图9 验证试验Figure 9 Experimental verification