基于近似多位姿的轮系式钵苗移栽机构运动综合

孙 良 胡艺翔 邢子勤 俞高红 俞亚新

(1.浙江理工大学机械与自动控制学院， 杭州 310018； 2.浙江省种植装备技术重点实验室， 杭州 310018)

0 引言

相比传统水稻毯面的机插技术，水稻钵苗移栽具有不伤根、无缓苗期、可缩短生长期等优点，对促进我国水稻产量增长及扩大可适种面积具有重要意义[1-3]。移栽机构是水稻钵苗移栽的核心工作部件，移栽机构创新设计及其设计理论一直是众多学者研究的重点。不等速轮系式机构具有结构对称、运动平稳等特点，学者将其应用在移栽机构设计中。

目前，轮系式水稻钵苗移栽机构的设计方法主要有运动分析与运动综合。运动分析是选定机构、建立机构运动学模型、确定机构关键参数，并优化机构参数，以获得理想的轨迹和姿态。俞高红等[4-5]在研究蔬菜钵苗取苗机构的基础上，采用运动分析方法提出一种取栽一体式的水稻钵苗移栽机构。孙良等[6]建立了夹秧片姿态和尖点轨迹的目标函数，获得一组满足水稻钵苗移栽要求的结构参数。叶秉良等[7]提出一种包含1个不完全非圆齿轮的非圆齿轮行星轮系水稻钵苗移栽机构，得到一组实现水稻钵苗移栽较优工作轨迹和姿态的机构参数。左彦军等[8]提出基于三次非均匀有理B样条曲线的非圆齿轮行星轮系钵苗移栽机构，并通过运动学分析建立了夹秧片姿态和尖点轨迹的目标函数，得到一组满足水稻钵苗移栽要求的结构参数。

在运动分析中对机构初始参数选择的盲目性会影响轨迹形状与姿态优化过程的收敛性，而运动综合可以在保证轨迹形状与移栽姿态的基础上，通过优化传动齿轮凸性来达到目标要求，便于初始参数的选取。赵雄等[9]采用封闭节曲线拟合技术，通过给定型值点确定机构作业要求的运动轨迹，建立机构逆运动学模型，求解非圆齿轮非匀速传动比函数，设计了非圆齿轮节曲线。SUN等[10]提出一种基于球面曲线的空间行星轮系机构运动综合方法，实现了非圆齿轮-非圆锥齿轮行星轮系宽窄行分插机构的参数反求。值得注意的是，现有轮系式移栽机构的运动综合中主要针对给定形状的轨迹求解机构参数，并获得满足凸性要求的传动比，求解模型中未考虑移栽臂的姿态因素。

考虑若干位姿(位置与姿态)信息的机构运动综合问题，即刚体导引问题，在连杆机构上已有许多研究。BURMESTER[11]从纯几何角度出发，研究了平面连杆的3、4、5个位姿的精确综合方法。SUH等[12]使用矩阵法建立了一系列非线性方程，求解平面四杆以及平面五杆最多五位姿精确综合，并对该连杆机构进行设计。SUBBIAN等[13]使用连续法对平面四杆6、7位置进行轨迹综合求解。文献[14-15]采用运动映射方法进行了平面四连杆机构的多位姿精确综合研究。ZHAO等[16-18]进行了连杆机构的混合多位姿综合研究，在五位姿无精确解的情况下，通过扩大误差的方法获得近似解。虽然目前基于多位姿的运动综合在连杆机构中应用较多，但应用在轮系机构中却鲜见相关报道。

本文基于运动学映射理论对不等速轮系机构进行运动综合，提出一种可兼顾轨迹形状和多位姿的水稻钵苗移栽机构设计方法。在求解非圆齿轮节曲线时采用遗传算法优化其圆度性。

1 水稻钵苗移栽轨迹

水稻钵苗移栽机构根据钵苗的取出方式可分为顶出式和夹取式，而夹取式水稻钵苗移栽机构根据夹取对象又分为夹苗式和夹钵式两种，其中夹苗式水稻钵苗移栽机构是主要研究对象[19-21]。夹苗式水稻钵苗移栽机构首先利用夹秧片在取苗点处夹取钵苗的底部茎秆，然后夹秧片沿着拔苗段轨迹将钵苗从穴盘中拔出，秧苗随着机构调整姿态并运动到移栽轨迹的最低点(植苗点)，松开夹秧片并通过推秧器将钵苗推入水田中以完成一次移栽。如图1所示，其轨迹是一条带环口的封闭曲线，而曲线的形状以及移栽臂的姿态影响秧苗的栽植效果[22]。

图1 “8字”型移栽轨迹位姿点选取示意图Fig.1 Schematic of selection of position points of “8-word” type transplanting trajectory

根据夹苗式水稻钵苗移栽机构的农艺要求以及机构运动特点，选取的轨迹姿态点一般要满足以下要求[23]：①取苗角α1大于-5°且小于15°。②拔苗段距离l1大于40 mm。③轨迹与秧箱之间的距离l2大于20 mm。④推苗角α2大于55°且小于75°。⑤取苗角和推苗角的差为55°。⑥取苗段轨迹应与钵盘近乎垂直。⑦动轨迹允许的植苗高度h不小于10 cm。⑧移栽机构齿轮箱回转最低点高于地面不小于2 cm。

结合上述要求，在轨迹的规划中考虑取苗环扣处移栽臂的位置和姿态对取苗成功率以及秧苗损伤率的影响[24]，故取苗环口处预设5个位姿点以约束取苗环扣的轨迹形状和相应的取苗姿态；持苗段的约束条件主要是控制轨迹与秧箱的合理间距，避免干涉，预设1个位姿点；植苗段的位姿点影响立苗效果，预设2个位姿点以控制植苗角以及轨迹的整体高度。即，共用8个位姿点约束整个轨迹的各关键部位(图1)。合理的移栽轨迹形状和姿态是指满足移栽要求的一个特定的数据区间[25]，并非是严格约束，为此，本文利用运动学映射理论开展近似8位姿的移栽机构运动综合。

2 多位姿近似综合方法

单行星架的轮系式移栽机构若不考虑齿轮约束，可以简化为如图2a所示的开链2R机构。开链2R移栽机构中移栽臂末端点的位置和姿态即为移动坐标系下的刚体在固定坐标系下的运动关系，如图2b所示。

图2 刚体位移在两个坐标系间的转换Fig.2 Transformation of planar displacement between two different coordinate systems

将刚体的位姿信息表示为(L1,L2,θ)，M表示与刚体相连的动坐标系，F表示固定坐标系，运动刚体内点P的位移可以表示为点坐标或线坐标从M到F的转换[26]，即

(1)

式中 (L1,L2)——动坐标系M的原点在固定坐标系F下的坐标

(x,y)——点P在动坐标系M下的坐标

(X,Y)——点P在固定坐标系F下的坐标

θ——动坐标系M相对于固定坐标系F的旋转角度，逆时针为正向

给定一组笛卡尔空间的平面位姿参数(L1,L2,θ)，用平面运动学映射理论可以将其转换为图像空间坐标B=(B1,B2,B3,B4)表示，其中

(2)

将式(1)代入式(2)，可得

(3)

如图2所示，两个平面刚体之间的相对运动关系可以看成是一个平面开链2R机构运动模型，动坐标系xy下的点P满足固定坐标系XY下圆方程的一般形式，即参数(X,Y)满足圆方程

2b1X+2b2Y+b3=b0(X2+Y2)

(4)

将式(3)代入式(4)可得

(5)

其中，p1=-b0,p2=b0y,p3=b0y,p4=b1,p5=b2,p6=-b1y+b2x,p7=-0.5(b1x+b2y),p8=0.25[b3-b0(x2+y2)]，(b0,b1,b2,b3)是约束圆的齐次坐标。另外，p1～p88个系数须满足附加方程组

(6)

根据给定的8位姿条件，式(5)可转换为

(7)

其中

通过对矩阵BTB进行奇异值分解，得到8个奇异向量，一般对于精确位姿求解可以选取3个对应的特征向量vα、vβ和vγ构成解空间。在解空间中，向量可以表示为

p=αvα+βvβ+γvγ

(8)

式中α、β、γ——解空间系数

然而，在许多实际应用等式中往往无法找到精确解，故允许使用更多近似解决方案进一步求解，即一种逐渐提高误差容忍度并扩展解空间的策略。如式(8)，这个近似的零空间是具有最小二乘代数拟合误差的局部最优解空间，为了允许更大的容错能力，可使用更高维度的解空间，即用4个奇异向量来构造一般解空间，以将一般解空间扩展到四维，表达式为

p=αvα+βvβ+γvγ+μvμ

(9)

式中μ——解空间系数

在四维解空间中，它包含4个要确定的齐次变量α、β、γ和μ，以及2个齐次二次方程。由于通解中有一个未知数没有确定，通常可以获得无穷多个解。为了满足实际的工程要求，还可施加其他具体的限制，如，可以将p第一个元素设置为0，以便获得圆柱副连接来提供更好的支撑性，但是最终的二杆组会落入RR或者RP类型中。因此，可以考虑使固定铰链点位于某一条直线上，即固定铰链点(Xf，Yf)在直线L1Xf+L2Yf+L3=0上。根据式(5)还可以得到

(10)

由式(10)可得

-L1p4-L2p5+L3p1=0

(11)

将式(11)代入式(9)可得

Cαα+Cββ+Cγγ+Cμμ=0

(12)

其中

(13)

式(13)是(α，β，γ，μ)的线性约束方程，它定义了所得的RR型二杆组固定铰链点所在的直线。该附加的线性约束方程可以与方程式(6)中的2个二次约束方程组合，可得

(14)

其中Kij是vα、vβ、vγ和vμ组成的表达式，其中8个约束移栽轨迹的位姿点见表1，求解得到解空间的4个特征向量见表2，最后利用求解二元二次方程组通法解得向量p的4种结果见表3。

表1 8个位姿点的参数Tab.1 Parameters of eight pose points

表2 解空间的4个特征向量Tab.2 Four eigenvectors of solution space

表3 向量p的4种结果Tab.3 Four results of vector p

将向量p标准化并附加实际的约束条件，即固定铰链点位于Y=0水平线上，然后结合2个齐次方程(式(6))，可以求解出未知参数α、β、γ和μ。对于N个平面姿态，可以由3种平面二杆开链类型中的一种进行近似综合：RR、PR或RP(图3)。

图3 平面二杆开链类型Fig.3 Type of planar two-bar open chain

依据向量p可以确定平面二杆开链的类型，4种模式为：①如果p1=p2=p3=p4=p5=0，得到的平面二杆组类型是PP型。②如果p1=p2=p3=0，得到的平面二杆组类型是PR型。③如果p1=p4=p5=0，得到的平面二杆组类型是RP型。④如无上述情况，生成的平面二杆组类型是RR型。

表3列出了p的4种结果，其中应排除p1=p4=p5=0(如表3中第2组解)的特殊情况，因为这会使得p始终满足式(11)。从而，可获得3个RR型二杆组，可依次近似经过给定的8个位姿。考虑到水稻钵苗移栽机构工作范围限制，选择第1、3二杆组构造机构。

3 平面开链2R机构及轨迹

将动铰链点P的坐标(x,y)转换成固定坐标系XY下的坐标(X,Y)，通过选取合适的二杆组，得到拥有2个固定铰链点的一组四杆机构。表4中，输入杆的长度是l1，输出杆的长度是l2，(Xf，Yf)和(Xm，Ym)分别是固定铰链点和动铰链点的坐标。图4中的四杆机构可以近似地通过8个给定的位姿点，其中实线部分的开链2R机构即为与轮系式移栽机构对应的行星架和移栽臂。

表4 四杆机构参数Tab.4 Parameters of four-bar mechanism mm

图4 四杆机构及其轨迹Fig.4 Four-bar mechanism and its trajectory

图5 原轨迹以及截取的部分轨迹段Fig.5 Original trajectory and part of intercepted trajectory

在获得的四杆机构轨迹中截取由8个给定位姿控制的曲线作为已知段轨迹，如图5所示，并结合选定的开链2R机构确定2R机构的相对角位移。为了得到整个运动周期的角位移曲线，先在求解的已知角位移中间隔选取若干数值点(本文中选取19个点，为确定的数值点)；再参考回程段轨迹(图5)对应的角位移选定若干数值点(本文选取5个点，为可变数值点)；然后，基于设定的24个数值点，利用三次非均匀B样条曲线拟合获得完整的角位移曲线(图6)。此时，再利用开链2R机构可复演出一条新的完整的封闭轨迹[27]。

图6 角位移曲线Fig.6 Angular displacement curve

4 移栽机构设计

4.1 传动比分配

根据运动综合得到的开链2R机构，在其第一杆上添加一定的传动齿轮后，可获得移栽机构传动所需的单自由行星轮系机构(图7a)。机构以与机架固定的太阳轮为中心，对称布置两套两级非圆齿轮传动，由行星架动力输入，通过齿轮传动，实现两移栽臂的复合运动，即移栽运动[28]。如图7b所示，将行星架输入点记为A0，行星架输出点记为Ac，行星架记为杆G1，移栽臂记为杆G2，移栽臂端点记为C。

图7 七齿行星轮系水稻钵苗移栽机构简化模型Fig.7 Transplanting mechanism and simplified model of seven-tooth planetary gear system rice bowl seedlings1.中心轴 2.太阳轮 3.中间轮A 4.中间轮B 5.行星轮 6.移栽臂部件 7.行星轴 8.中间轴 9.齿轮箱

显然，轮系机构的总传动比以及各级齿轮传动比可依据上一节拟合的角位移曲线(图6)获得。

对于两级非圆齿轮轮系机构，总传动比为

(15)

式中i——机构总传动比

w0——太阳轮角速度

w1——一杆(输入杆)角速度

w2——二杆(输出杆)角速度

确保两对齿轮的传动比相近，分配的传动比计算式为[29]

(16)

i2=i/i1

(17)

式中i1——第1级非圆齿轮的传动比

i2——第2级非圆齿轮的传动比

k——峰谷值调整系数

图8为总传动比、第1级和第2级传动比曲线，曲线峰值位置和整体波动趋势相近，有利于获得圆度较好的非圆齿轮节曲线。图9为根据初始分配的传动比求解的两对非圆齿轮节曲线。显然，从节曲线形状上看，非圆齿轮整体光滑性有待进一步改善以利于齿廓设计。

图8 传动比分配结果Fig.8 Distribution of transmission ratio

图9 未经优化的非圆齿轮节曲线Fig.9 Unoptimized non-circular gear pitch curves

4.2 非圆齿轮的凹凸性优化

针对图6所示的角位移曲线，可通过改变点17、18之间的5个可变的角位移参数来调节齿轮节曲线以及非关键轨迹段的形状。将非圆齿轮节曲线的凹凸性优化转变为目标函数，变量为图6中5个控制点的横坐标，同时满足条件

λ17≤X1≤X2≤X3≤X4≤X5≤λ18

其中，独立变量X1～X5分别表示可变型值点1～5(图6)的横坐标，λ17、λ18表示固定型值点17、18(图6中实心点)的横坐标。

建立一种通用的目标函数或评价函数，其目标是优化非圆齿轮节曲线的圆度性，以提升齿轮设计可行性与运动学性能。具体函数建立方式如下：若太阳轮节曲线中一个拟合点的向径小于相邻两个拟合点的向径时，函数Y1表示为

(18)

式中r1i——太阳轮节曲线上第i个拟合点的向径

r1(i+20)、r1(i-20)——第i个拟合点两侧拟合点的向径

N——一条节曲线上拟合点的总数

第二级非圆齿轮节曲线函数Y2形式为

(19)

式中r3j——中间轮3上第j个拟合点的向径

r3(j+20)、r3(j-20)——中间轮3上第j个拟合点两侧拟合点的向径

r4j——行星轮4上第j个拟合点的向径

r4(j+20)、r4(j-20)——行星轮4上第j个拟合点两侧拟合点的向径

由于存在两个需要优化的目标，一般来说，分目标函数越多，对设计方案的评价就越全面，但计算也越复杂、耗时。对于多目标优化问题，通常的处理方法是将几项指标综合进行评价，即构造一个多目标函数

f(X1,X2,…,X5)=K1Y1+K2Y2

(20)

式中K1、K2——权重系数

目标函数f(X1,X2,…,X5)由函数Y1和Y2组成。

表5为利用遗传算法优化得到的一组最优值，其上限值是图6中角位移曲线上的实心点18的横坐标，其下限值为图6中角位移曲线上的实心点17的横坐标。利用遗传算法的局部优化综合方法针对角位移曲线上的部分可变的型值点进行优化，其优化前后角位移曲线的对比如图10所示。优化后的角位移曲线在部分曲线段上相比于初始角位移曲线更为平缓。

表5 角位移的最优值

图10 优化前后的角位移曲线Fig.10 Comparison of angular displacement curves before and after optimization

通过多次优化测试得到多目标函数中的权重系数为0.95和0.05。在优化中，随机产生的初始个体为30个，并包含有50代遗传代数，每代的最优适应度如图11所示，最优值为968.115，平均值为968.402，适应度函数呈收敛状态。对比图9和图12可知，通过遗传算法参数优化明显改善了非圆齿轮节曲线的凹凸性。

图11 不同代数时的最优适应度Fig.11 Optimal fitness values per generation

图13为对非圆齿轮进行参数优化前后的移栽轨迹，平面开链2R机构的轨迹发生了比较明显的变化，但由于变化的位置仅限于可变轨迹段，其关键轨迹段仍然近似经过8个给定的位姿点，所以该轨迹依旧能够满足第2节提出的设计要求。

图12 优化后的非圆齿轮节曲线Fig.12 Optimized non-circular gear pitch curve

图13 优化前后的轨迹对比Fig.13 Trajectory comparison before and after optimization

5 移栽机构仿真与试验

为了进一步验证水稻钵苗移栽机构运动综合的数据合理性，利用软件仿真和高速摄影机(Phantom系列)拍摄，对移栽臂通过的各个关键位姿点及轨迹进行分析。

由图14、15可知，实际物理样机的轨迹与虚拟样机仿真分析的轨迹基本一致。表6为8个设定位姿点相对应的最接近轨迹点的位姿数据。三者姿态角(与x轴方向的夹角)最大误差控制在1°范围内，符合近似运动综合设计要求。

图14 仿真轨迹Fig.14 Simulation trajectory

图15 物理样机轨迹Fig.15 Physical prototype trajectory

表6 姿态角对比Tab.6 Attitude angle comparison (°)

在移栽机构试验台取苗试验中，秧苗采用14格×29格钵盘育苗，平均每穴1～3株秧苗，苗龄为25 d，秧苗平均高度为156 mm。移栽机构在试验台架上以80 r/min的速度运行，每分钟能取苗160次，取苗成功率为94.98%，如图16所示。而后，开展田间移栽试验，如图17所示。栽植的秧苗在水田中的翻倒现象有所改善，相比于现有的轮系式水稻钵苗移栽机构，其移栽后秧苗的立苗角度有所提升。

图16 取苗试验Fig.16 Picking seedling test

图17 田间试验Fig.17 Field test

6 结论

(1)提出一种基于近似多位姿的非圆齿轮行星轮系综合方法，完成了一种7个非圆齿轮水稻钵苗移栽机构的设计，并且实现了满足多个给定位姿点的“8字”型复杂轨迹；提出了基于遗传算法的非圆齿轮节曲线非圆度优化方法，在不改变关键轨迹段的情况下优化得到了合理的非圆齿轮节曲线。

(2)通过仿真与试验，验证了综合方法的正确性以及移栽机构设计的可行性，室内取苗试验成功率达到了预期要求。完成了田间移栽试验，栽植的秧苗在水田中的翻倒现象有所改善。