樊桂菊 李 釗 毛文華 婁 偉 梁 昭 姜紅花
(1.山東農業大學機械與電子工程學院, 泰安 271018; 2.山東省農業裝備智能化工程實驗室, 泰安 271018;3.山東省園藝機械與裝備重點實驗室, 泰安 271018; 4.中國農業機械化科學研究院, 北京 100083;5.山東農業大學信息科學與工程學院, 泰安 271018)
水果種植屬于勞動密集型產業。為降低剪枝、疏花疏果、采摘等作業環節的勞動強度,提高工作效率,學者們研制了多種類型輔助人工操作的果園作業平臺[1-4],根據升降機構的不同,作業平臺可分為剪叉式和臂式。水果種類繁多,各地種植模式不同[5],果園作業平臺的工作空間直接影響其適用性。因此,開展基于工作空間的果園作業平臺結構參數優化研究對其推廣應用具有重要意義。
針對考慮工作空間的結構參數優化問題,吳超宇等[6]以提出的全局混合性能指標為目標函數,利用粒子群算法優化得到直線驅動型并聯機器人的最優尺寸參數。權龍哲等[7]以包容目標工作空間及桿長因素為目標,采用遺傳算法優化立體苗盤管理機器人的結構參數。丁淵明等[8]提出以工作空間和能量消耗綜合指標為目標函數,使用遺傳算法優化機械臂結構參數。孫小勇等[9]通過分析并聯機構結構尺寸對工作空間的影響確定了結構參數的優化范圍,利用fminimax函數求解得到優化后的結構參數。
基于工作空間的結構參數優化主要集中在機器人或機械臂等方面,針對果園機械的相關研究較少。本文基于已研制的果園作業平臺,以喬砧密植的紡錘形蘋果種植模式確定其目標工作空間,采用網格迭代算法計算可達工作空間體積,分析平臺結構參數對其空間性能的影響規律,建立以空間性能和結構緊湊為指標的多目標優化函數,利用遺傳算法求解最優參數,并通過仿真與試驗進行驗證,以期提高果園作業平臺的適用性和可操作性。
前期研制的果園作業平臺主要由動力裝置、行走機構、回轉機構、升降機構、調平機構、工作臺和控制系統等組成,實現工作臺升降、旋轉和坡地調平等功能。其結構、工作原理和運動學方程見文獻[10],但未考慮操作人員在工作臺的左右和前后移動。為更好地描述平臺可操作性和各機構運動,將原參考點在豎直方向上移動至成人肩高[11]位置作為平臺執行末端參考點,建立D-H坐標系[12-13],如圖1所示。其中,l1~l5為各桿長度,分別為964、812、210、900、680 mm。各參數如表1所示,表中αi、ai、di、θi(i為關節編號,i=1,2,…,6)分別表示連桿轉角、連桿長度、連桿偏移和關節角。d5min為工作臺鉸接點安裝距離,為70 mm;d6max與d6min互為相反數,使工作臺關于橫梁對稱以保持平衡。其運動學方程為
(1)
其中
P=[pxpypz]T
式中R——姿態矩陣P——位置矩陣
n、o、a、p——關于θi與di的函數
O——零矩陣

表1 果園作業平臺D-H參數Tab.1 D-H parameter of orchard platform
果園作業平臺的目標工作空間是其結構參數設計和優化的重要依據,與果園地形、果樹種類和種植模式密切相關。而我國果樹種類繁多,種植模式各異,為使作業平臺能夠靈活、高效地在復雜環境中進行升降、旋轉和調平等動作,實現輔助人工完成剪枝、疏花疏果、套袋及果實采摘等作業環節,結合作業平臺的通用性,本文以種植面積位居前列的蘋果為研究對象,目前其種植模式仍以喬砧密植為主,樹形多為紡錘形[14-15],株距為3~4 m,行距4~5 m,株高為2.8~3.5 m,主干高度為0.6~1.0 m,冠徑為1.4~3.0 m。以O0為坐標原點,X軸為行距方向,Y軸為平臺前進方向,Z軸為株高方向,平臺參考點的目標工作空間如圖2所示。其中,陰影部分為平臺參考點的目標工作空間;Dtx、Dty、Dtz分別表示目標工作空間在X、Y、Z方向范圍,分別為[-((L0-Ld)/2-Lr),(L0-Ld)/2-Lr]、(0,D0/2-Lr]、[Hr+Lr,Hs-Lr],其中L0為行距,Ld為冠徑,D0為株距,Hr為成人肩高,Hs為株高,Lr為成人掌心與肩膀距離。
為使平臺最大限度地滿足作業需求,各參數取值為:L0、D0、Hs分別取最大值,冠徑取最小值,由文獻[11]可知,成人掌心與肩膀距離取0.5 m,肩高取1.4 m。得目標工作空間在X、Y、Z方向范圍分別為:Dtx為[-1 300 mm,1 300 mm]、Dty為(0,1 500 mm]、Dtz為[1 900 mm,3 000 mm]。
平臺參考點的可達工作空間指在各關節約束下參考點運動學方程所有解的集合。本文采用改進的蒙特卡洛法[10]得到其工作空間點云圖,如圖3所示。由圖3可知,參考點可達工作空間在X、Y、Z方向范圍分別為:Drx為[-1 781 mm,1 802 mm]、Dry為[-730 mm,1 819 mm]、Drz為[1 190 mm,2 917 mm]。
由上文可知,作業平臺參考點的可達工作空間與目標工作空間存在偏差,X方向和Y方向偏差可以通過調整平臺底盤位置消除,因此本文僅研究Z方向偏差。
對比圖2和圖3,為保證操作人員對兩行果樹有效作業,以μz1和μz2描述平臺工作空間尺寸偏差,計算式為
(2)
式中Zr1max、Zr1min——X為-1 300 mm時可達工作空間中Z的最大值、最小值,mm
Zr2max、Zr2min——X為1 300 mm時可達工作空間中Z的最大值、最小值,mm
Ztmax、Ztmin——目標工作空間在Z方向的最大值、最小值,mm
工作空間體積[16]是表述作業平臺工作能力的重要參數。為得到較高精度的空間體積,本文采用網格迭代算法[17]。設某子網格編號為Numi,其相鄰網格編號如圖4所示,算法流程如圖5所示。
通過多次迭代計算平臺參考點的工作空間體積,體積迭代差計算式為
(3)
式中Vj、Vj-1——第j、j-1次迭代工作空間體積,m3
以體積迭代差ε≤0.1%為體積收斂標準。圖6為工作空間體積迭代變化曲線,第16次迭代時體積迭代差為0.045%,故取工作空間體積為4.26 m3。
為量化作業平臺在X、Y、Z方向的運動能力,以平臺可操作性作為評價參數。目前關于機器人可操作性的研究方法主要有:基于雅可比矩陣、基于Hessian矩陣和基于剛度矩陣[18-20],其中由于雅可比矩陣代表了機器關節空間的微分運動向操作空間的微分運動之間的映射轉換,而被廣泛應用。因此,本文在平臺參考點運動學模型基礎上推導雅可比矩陣,計算式為
(4)
其中
D=[dxdydzδxδyδz]TJ=[J1J2J3J4J5J6]Θ=[θ1θ2θ3θ4d5d6]T
(5)
式中D——參考點沿X、Y、Z軸的微分平移和繞X、Y、Z軸的微分旋轉矩陣
J——雅可比矩陣Θ——關節空間
基于雅可比矩陣有多種指標可以描述平臺可操作性,其中可操作度κ考慮了機構執行末端的各方向運動能力而優于其他指標[21],計算式為
(6)
式中λi——矩陣JJT的特征值
由式(4)、(6)可知,當平臺結構參數給定后,可操作度隨平臺參考點空間位置變化而不同,為描述作業平臺在工作空間的全局操作性能,引入平均可操作度τ,計算式為
(7)
其中Ω={(x,y,z)|x∈Drx,y∈Dry,z∈Drz}
式中Vi——平臺可達工作空間子網格空間體積,m3
V——平臺可達工作空間體積,m3
Ω——平臺可達工作空間范圍,m
通過分析果園作業平臺運動學模型和其工作空間,影響其空間性能的結構參數主要是桿長和關節變量(關節角和連桿偏移)。根據其結構特點及作業要求,主要分析桿2、4長度,關節1、3關節角,以及關節5、6連桿偏移對空間性能的影響規律。
2.4.1桿長
保持關節變量和桿1、3、5長度不變,桿2、4長度在0~1 600 mm之間變化,通過蒙特卡洛方法多次模擬可達工作空間,由2.1~2.3節分別得桿2、4對可達工作空間尺寸偏差μz1和μz2、體積V及平均可操作度τ的影響,如圖7所示。
由圖7可知,桿2、4對μz1和μz2影響趨勢一致,桿2影響更大:當桿2長度約為950、960 mm時,μz1和μz2最小,分別為42.82、40.05 mm;當桿4長度約為800、810 mm時,μz1和μz2最小,分別為260.49、258.97 mm;桿2對V和τ影響較小,桿4對V和τ影響較大且呈正相關。
2.4.2關節變量
保持桿長和關節2、4關節角不變,關節1、3的關節角分別在[60°,300°]和[-130°,-50°]范圍內變化;由表1可知,d5min為70 mm,|d6max|=|d6min|,故僅考慮d5max、d6max對工作空間的影響,關節5、6的連桿偏移最大值在[0,1 200 mm]范圍內變化。得關節1、3的關節角和關節5、6的連桿偏移對可達工作空間性能的影響,如圖8~10所示。
由圖8~10可知,θ1對μz1、μz2影響較小,當θ3為-111°~-63°時,μz1和μz2最小,分別為62.18、54.74 mm;θ1與θ3對V影響較大,對τ影響較小但均呈正相關。當d5max大于700 mm時,μz1和μz2基本不再變化;當d6max大于800 mm時,μz1和μz2基本不再變化;d5max對V和τ均影響較大呈正相關;d6max對V影響較大呈正相關,對τ影響較小呈負相關。
針對復雜的果園作業環境空間,為使作業平臺能夠安全可靠地完成動作,在能夠滿足目標工作空間的作業要求下,執行升降、旋轉和調平任務時更加靈活、高效,可達工作空間與目標工作空間的尺寸偏差應盡量小,空間體積和平均可操作度應盡量大。因此可達工作空間性能目標優化函數為
(8)
式中w——待優化參數矩陣
同時,為使作業平臺結構更加緊湊,在滿足目標工作空間作業任務前提下,平臺結構尺寸和關節變量應盡量小,即總桿長和關節變量總和應盡量小。因此結構緊湊目標優化函數為
(9)
式中ζ1、ζ2——量綱統一系數
根據2.4節的分析和上述目標優化函數涉及的變量可得出平臺待優化的參數為桿2、4長度,關節1、3關節角范圍和關節5、6連桿偏移范圍,用矩陣w表示為
w=[l2l4θ1maxθ1minθ3maxθ3mind5maxd6max]T
(10)
綜上所述,作業平臺結構參數優化函數共有6個子目標,8個待優化參數,屬于多元多目標函數優化問題,引入比例因子和權重,將其變為單目標優化函數;由圖7~10得到待優化結構參數范圍作為約束條件,優化模型為
(11)
其中
∑ηi=1
式中fi(w)——各子目標函數
F(w)——總目標函數
si——比例因子ηi——權重
根據平臺原有結構和實際作業需求,對優化模型的各指標設定適當的比例因子和權重。
3.2.1比例因子
工作空間尺寸偏差越小越好,考慮優化精度和實際作業情況,其比例因子取10 mm;工作空間體積、平均可操作度及結構參數指標的比例因子均以原樣機結構為基礎,分別按照前述對應內容計算求解,結果如表2所示。

表2 子目標函數指標權重及比例因子Tab.2 Index weight and scale factor of sub-objective function
3.2.2權重
基于果園種植模式和優化模型,在能夠完成作業任務的前提下,考慮農機農藝融合和人機工程學,各指標權重確定原則如下:①空間性能指標較結構緊湊指標更為重要。②在空間性能指標中,工作空間尺寸偏差最重要,μz1、μz2為同一性質優化指標,二者同等重要。③為保證平臺靈活高效地完成執行動作,平均可操作度較工作空間體積更為重要。④通過分析圖8~10中桿長和關節變量對工作空間尺寸偏差的影響,前者大于后者,因此桿長較關節變量重要。
以工作空間性能和結構緊湊為一級指標,6個子目標函數為二級指標。在上述原則前提下,采用德爾菲法[22]征詢相關領域專家意見,同時結合層次分析法[23-24]的分析計算,最終確定各評價指標的權重。以空間性能(A1)為準則的二級指標μz1(B11)、μz2(B12)、V(B13)和τ(B14)為例,其步驟如下:
(1)將專家對各指標重要性評分匯總,依據判斷矩陣標度[25]對指標B11~B14進行兩兩判斷,構造出判斷矩陣EA1,計算式為
(12)
(2)將上述矩陣每一行元素相乘得積為PA1i,對其歸一化處理,即得權重ηA1i計算式為
(13)
(3)為保證專家對指標B11~B14重要程度的判斷一致,以一致性比率CR檢驗判斷矩陣的一致性,計算式為
(14)
式中RI——隨機指標,判斷矩陣4階時RI取0.89[26]
根據以上步驟,將判斷矩陣EA1和權重ηA1代入式(14),得CR=0.057<0.1,即判斷矩陣通過一致性檢驗,ηA1=[0.355 0.355 0.135 0.155]T。
同理,按上述步驟計算其余指標權重,結果如表2所示。
將表2的數據代入式(11),并將其編譯為適應Matlab遺傳算法工具箱的適應度函數,將優化模型的約束條件編譯為適應于算法工具箱的區域描述器。設置初始種群個體數目為40,遺傳代數為200,交叉概率為0.7,變異概率為0.007。運行算法得到適應度函數進化曲線如圖11所示。
由圖11可知,適應度函數約在第50代收斂。考慮實際加工因素,規整優化參數,如表3所示。
將優化前后參數代入運動學模型,模擬可達工作空間,分別計算其空間性能指標、總桿長和關節變量總和,結果如表4所示。

表3 優化前后結構參數對比Tab.3 Comparison of structural parameters before and after optimization

表4 優化前后指標對比Tab.4 Comparison of index before and after optimization
由表4可知,優化后可達工作空間最大偏差僅為12.33 mm,平均可操作度提高1.43%,表明優化后更接近目標工作空間,且操作更靈活。
4.1.1試驗準備
基于優化后的果園作業平臺結構參數改進樣機,于2020年7月在山東農業大學園藝實驗基地果園(圖12)進行工作空間試驗。該果園為喬砧密植蘋果園,樹形為自由紡錘形,株行距為2 m×3 m,平均株高約為3.5 m,試驗人員肩高為1.41 m,體質量為65 kg。
試驗儀器包括:直川電子科技有限公司ZCT 230M型傾角儀,Panasonic公司HG-C1100型激光位移傳感器,杰科斯JK-100F系列土壤水分儀,史丹利STHT 36191-23型鋼卷尺和李寧AQAP322型秒表等。
4.1.2試驗設計
以D-H坐標原點為實際坐標系原點,工作臺保持水平,以平臺承載質量、橫坡坡度及縱坡坡度為影響因素進行工作空間尺寸偏差測量試驗,設計三因素三水平試驗,如表5所示。

表5 試驗因素水平Tab.5 Factors and levels


(k=1,2;i=1,2,…,9)
(15)
(16)
式中y′ik——第k項指標第i指標值的評分
yik——第k項指標的第i指標實測值
yMk、ymk——第k項指標最大、最小實測值
Rk——第k項指標的實測值極差
ηk——第k項指標的權重,由表2給出

表6 試驗方案與結果Tab.6 Test scheme and result
(1)在果園作業平臺前期研究基礎上,結合農機農藝融合確定了其目標工作空間,采用網格迭代算法求取可達工作空間體積,以工作空間尺寸偏差、工作空間體積和平均可操作度描述可達工作空間性能,分析了平臺結構參數桿長、關節角和連桿偏移對可達工作空間性能的影響。
(2)建立了以可達工作空間性能和結構緊湊為目標的優化模型,通過引入各指標比例因子和權重轉換為單目標函數,利用遺傳算法求解得出最優結構參數為:桿2和桿4的長度分別為988、879 mm,關節1、3的關節角范圍分別為[107°,256°]、[-118°,-76°],關節5、6的連桿偏移最大值分別為720、340 mm。優化后工作空間尺寸偏差分別減小96.09%、95.60%,平均可操作度增加1.43%,表明優化后更接近目標工作空間,且操作更靈活。
(3)進行了不同承載質量、不同坡度下的優化樣機實地果園工作空間試驗,當平臺承載質量65 kg、橫坡坡度15°、縱坡坡度15°時,實際工作空間尺寸偏差最大值僅為16.7 mm,較原型減小了93.76%。