許林云，劉冠華，周 杰，周宏平，蔣雪松，宣 言

（南京林業(yè)大學(xué)機械電子工程學(xué)院，南京 210037）

0 引 言

現(xiàn)有林果采收的主要方式為通過林果采收設(shè)備對果樹主干或側(cè)枝處施加振動激勵，振動經(jīng)過樹枝傳遞到果實并引起果實振動響應(yīng)，當(dāng)果實的慣性力大于果實的果柄結(jié)合力時果實脫落[1]。研究果樹的振動采收機理有助于開發(fā)林果振動采收裝備。在針對果樹的動態(tài)響應(yīng)研究中，最常見的研究方法為動力學(xué)建模[2-5]和試驗研究分析[6-10]。果樹動力學(xué)建模方面，Miller[11]將樹體簡化為集中質(zhì)量-連桿模型，并將杜芬方程與樹木振動模型結(jié)合起來，對樹木的非線性響應(yīng)進行分析。Kenneth 等[12]通過多自由度的質(zhì)量-彈簧-阻尼系統(tǒng)構(gòu)建果樹振動模型，將一級側(cè)枝模型與主干模型相連接，二級側(cè)枝模型與對應(yīng)的一級側(cè)枝模型相連接，研究了多級樹枝振動。Murphy 等[13]構(gòu)建了多級連桿模型，將多連桿模型果樹振動方程在振動平衡位置附近簡化為線性方程，求解出了多連桿果樹振動模型的固有頻率和振型，并利用模態(tài)分析方法計算了果樹的響應(yīng)。王琳[14]根據(jù)云杉的生態(tài)特性將其簡化為一端固定、一端自由的變截面彈性桿模型，將樹干的部分質(zhì)量與樹冠質(zhì)量看作一個集中質(zhì)量球固定在桿的頂端，求得了云杉的自由振動響應(yīng)形式。翁凌云等[15]把Y 型果樹簡化為變截面楔形梁，將梁頂端施加一個分段樹干質(zhì)量與樹冠質(zhì)量的等效質(zhì)量塊，求得了其受迫振動下的響應(yīng)。

近年來許多學(xué)者利用三維建模與有限元法對果樹進行模態(tài)和振動響應(yīng)分析[16-19]。Bentaher 等[20]通過對樹木的樹形結(jié)構(gòu)測量，利用COSMOSM-(Geostar)數(shù)值模擬軟件對果樹進行了三維重構(gòu)，并求解了果樹在不平衡偏心塊式激振下的響應(yīng)情況。Peng 等[21]用三維軟件構(gòu)建了冬棗果樹的三維模型，實測了冬棗樹的密度、彈性模量和剪切模量，利用ANSYS 軟件對冬棗果樹模型進行了模態(tài)分析和諧響應(yīng)分析。Burt 等[22]通過QSM 點云處理軟件對威薩姆森林的點云進行處理，利用Abaqus軟件構(gòu)建樹木風(fēng)振有限元模型，通過數(shù)值模擬得到了在風(fēng)振時樹木的應(yīng)力分布。賀磊盈[23]采用2 張不同角度的核桃樹圖像對核桃樹進行了三維重建，并對重建的三維模型模擬果樹側(cè)枝修剪，用有限元軟件對不同程度修剪的果樹進行模態(tài)計算，結(jié)果發(fā)現(xiàn)隨著側(cè)枝樹木數(shù)量的增加，固有頻率出現(xiàn)下降趨勢。王冬等[24]通過Pro/E 軟件建立了3 種常見形態(tài)的果樹有限元模型，計算了3 種不同激勵下的果樹響應(yīng)，確定了3 種樹型最合適的激振采收方式。

綜上所述，傳統(tǒng)的建模方式只適合將復(fù)雜多形態(tài)的果樹極大簡化后建模，有限元建模方法雖可構(gòu)建與復(fù)雜形態(tài)果樹的模型，但現(xiàn)有的建模方法均忽略了樹葉和果實對果樹的影響，而理論模型應(yīng)對與實際復(fù)雜形態(tài)果樹較一致的有果有葉果樹建模才有實際應(yīng)用價值及指導(dǎo)意義。本研究提出了將激光點云信息與有限元相結(jié)合，提取果樹的骨架點及節(jié)點樹枝半徑，構(gòu)建果樹空間有限元模型。基于同種果樹在其收獲期間果實與樹葉的分布規(guī)律基本一致的假設(shè)，通過對果實和樹葉在果樹上的分布進行統(tǒng)計分析，確定果實和樹葉在果樹上的分布規(guī)律，繼而對模型節(jié)點上附加果實質(zhì)量和樹葉質(zhì)量構(gòu)建有果有葉的果樹空間振動理論模型，并用銀杏樹對模型進行試驗驗證及分析。

1 果樹空間有限元振動模型

1.1 有果有葉果樹振動模型構(gòu)建方法

將樹枝看作空間圓柱梁，采用空間6 自由度圓柱梁進行三維建模。假設(shè)果樹在小變形振動下為線性振動，可將果樹的振動視為空間剛體振動，并將機械振動理論中的三維空間圓柱梁結(jié)構(gòu)有阻尼一般振動微分方程應(yīng)用于果樹，形成果樹振動微分方程如式（1）所示

式中K、C 和M 分別為果樹的剛度矩陣、阻尼矩陣和質(zhì)量矩陣，其中質(zhì)量矩陣M 由樹枝質(zhì)量矩陣Mz、樹葉質(zhì)量矩陣My和果實質(zhì)量矩陣Mg三部分構(gòu)成，果樹剛度矩陣K 由樹枝剛度構(gòu)成；、和u 分別為節(jié)點加速度向量、速度向量和位移向量；p(t)為激勵力向量。

依據(jù)有限元理論[25]可知6 自由度空間圓柱梁單元的質(zhì)量矩陣和剛度矩陣分別如對稱矩陣（2）和對稱矩陣（3）所示

式中r 為對應(yīng)節(jié)點圓柱梁的半徑，mm，其所對應(yīng)單元的截面積為A=πr2，mm2；m 為對應(yīng)單元的單元質(zhì)量，kg； L 是單元長度，mm；梁的拉伸剛度Hm=EA，N，其中E為果樹的彈性模量，MPa；梁的扭轉(zhuǎn)剛度為Ht=Gmr2/2，N·m2，其中G 為果樹的剪切模量，MPa；繞Y 軸的抗彎剛度 Hy=Em(3r2+L2)/12，N·m2；繞 Z 軸的抗彎剛度Hz=Em(3r2+L2)/12，N·m2；Sym 表示省略部分關(guān)于對角線對稱。

通過果樹激光點云提取果樹的骨架節(jié)點構(gòu)建的果樹振動模型，較人工測量手繪模型可以更好地重現(xiàn)果樹的枝桿姿態(tài)。以激光掃描提取的骨架點為基礎(chǔ)，將2個具有相鄰?fù)負潢P(guān)系的骨架點作為空間梁單元的2 個節(jié)點構(gòu)成一個空間梁單元，樹枝振動模型由已構(gòu)建的單個單元按照拓撲關(guān)系依次連接而成，連接關(guān)系如圖1所示。

圖1 果樹骨架模型示意圖 Fig.1 Schematic diagram of fruit tree skeleton model

實際的果實采收期間，果樹枝桿上有大量的樹葉和果實，樹葉和果實的質(zhì)量會影響果樹整體的質(zhì)量分布，只使用果樹枝干模型計算得到的模態(tài)參數(shù)與有果有葉的果樹模型計算的結(jié)果存在很大差異。本研究提出根據(jù)果樹的果實和樹葉分布規(guī)律構(gòu)建果樹有果有葉振動模型?，F(xiàn)有的研究中，已經(jīng)有部分學(xué)者采用統(tǒng)計的方式對植物的生態(tài)特性進行了研究[26-27]。本研究設(shè)單位長度的樹枝上果實和樹葉的平均數(shù)量分別為ng和ny。單個果實的平均質(zhì)量為mg，kg；單個樹葉的平均質(zhì)量為my，kg。則附加果實與樹葉質(zhì)量矩陣Mg和My如式（4）和式（5）所示

附加樹葉與果實質(zhì)量后的總質(zhì)量矩陣如式（6）所示

式（1）中的阻尼矩陣表達了樹對運動能量的耗散，樹木的阻尼是由樹木振動過程中樹木內(nèi)部和外部摩擦共同作用的結(jié)果。Milne[28]對云杉的阻尼進行了分析，其發(fā)現(xiàn)樹冠與鄰近樹木的碰撞占總阻尼的50%，其次為葉子運動時對應(yīng)的空氣阻尼占總阻尼的40%，樹體自身的黏塑性阻尼占總阻尼的10%。Sellier 等[29]在研究樹木風(fēng)振時采用瑞利阻尼構(gòu)建果樹有阻尼振動模型。Moore[30]認為樹體阻尼與枝條和樹葉在樹枝上的分布有關(guān)，且枝條和樹葉實際上影響了果樹樹枝的質(zhì)量分布，因此采用瑞利阻尼進行計算是合理的。綜上所述，本研究采用瑞利阻尼構(gòu)建果樹振動模型。瑞利阻尼C 一般式如式（7）所示

式中α 為質(zhì)量阻尼系數(shù)，β 為剛度阻尼系數(shù)。果樹的瑞利阻尼系數(shù)可通過實測獲取果樹的多階固有頻率及對應(yīng)的黏性阻尼系數(shù)，采用最小二乘法進行擬合，即可由式（8）確定α、β。

式中ωi為第i 階固有圓頻率，rad/s；ξi為第i 階固有圓頻率對應(yīng)的黏性阻尼系數(shù)。

則第i 階固有頻率如式（9）所示

式中fi表示第i 階固有頻率，Hz。

1.2 計算有果有葉果樹固有頻率

由機械振動原理可知，系統(tǒng)有阻尼固有頻率ωd和無阻尼固有頻率ωn的關(guān)系如式（10）所示

式中ωd為有阻尼固有圓頻率，rad/s；ωn為無阻尼固有頻率，rad/s；ξ 為黏性阻尼系數(shù)。

許多學(xué)者的試驗研究表明樹木的阻尼比大約在0.2左右[28,31]在低頻區(qū)有阻尼固有頻率與無阻尼固有頻率十分接近。由于果樹實際阻尼很小，在計算果樹固有頻率時，可以忽略果樹的阻尼項，在不施加外部激勵時，則式（1）可簡化為果樹無阻尼自由振動微分方程如式（11）所示

根據(jù)線性振動理論中的模態(tài)振型定義，當(dāng)果樹處于某一模態(tài)振型? 時，果樹各點的振幅之比是唯一的。設(shè)果樹在某一模態(tài)振型下振動時各點的位移振動響應(yīng)如式（12）所示

式中ω 表示響應(yīng)頻率，rad/s；θ 表示響應(yīng)與激勵之間的相位差，rad；t 表示響應(yīng)時間，s。

將式（12）帶入式（11）中約去共同因子后得到式（13）

求解式（13）需給定邊界條件，果樹的邊界條件主要為根部對果樹的約束。Láng[32]在進行果樹有限元分析時對樹木的根部簡化為彈簧-阻尼系統(tǒng)。Dupuy 等[33]構(gòu)建了包括樹根在內(nèi)的有限元模型，要構(gòu)建此類模型需要根部的幾何形狀，要將樹根挖掘出來并進行測量，這樣獲取根部的形狀是很困難的。Sillier 等[29]將根部簡化為2 個扭轉(zhuǎn)彈簧，保證果樹能夠繞水平面的2 個坐標(biāo)軸轉(zhuǎn)動。Wood 等[34]認為根部雖不能使樹木非常牢固地固定在地面上，但這些學(xué)者認為這種現(xiàn)象對樹木的固有頻率影響不大。綜上所述，構(gòu)建樹根的約束需在根部節(jié)點施加對應(yīng)的約束條件，本研究考慮了Wood 等[34]的結(jié)論與實際計算的便捷性后，將根部約束簡化為固定端約束。消去果樹底部節(jié)點的6 個自由度后求解式（13）的特征行列式即可得到6 自由度三維空間梁單元構(gòu)建的果樹模型的固有頻率。

2 銀杏樹三維模型構(gòu)建

為驗證上節(jié)所構(gòu)建的果樹振動模型的準(zhǔn)確性，本研究以銀杏樹為研究對象，在南京林業(yè)大學(xué)校園內(nèi)采伐一棵有果有葉的適用于室內(nèi)試驗的小型銀杏樹（圖2a）進行模型構(gòu)建及相關(guān)試驗測試，采伐時間為2019 年9 月，正值果實成熟時期。采伐后的果樹用地鉗將果樹根部固定在地面上。為研究有果有葉、無果有葉和無果無葉不同狀態(tài)果樹的頻譜特性，采用人工干預(yù)的方式，先從有果有葉狀態(tài)時的果樹進行測試，然后人為摘去果實形成無果有葉果樹，最后摘去樹葉形成無果無葉果樹（圖2b）進行各項試驗測試，所有室內(nèi)試驗在果樹采伐后的3 d 內(nèi)完成。

要精準(zhǔn)確定果樹骨架點，構(gòu)建果樹骨架振動模型，只能先對無果無葉銀杏樹通過激光掃描獲取骨架點云，然后再在骨架振動模型基礎(chǔ)上，通過分別添加果實與樹葉分布陣，獲得有果有葉的果樹振動模型。

2.1 銀杏樹樹枝點云處理及模型構(gòu)建

本研究采用二維激光雷達掃描圖2b 所示的果樹，獲取果樹骨架點云數(shù)據(jù)，點云采集系統(tǒng)（圖3a）包括二維激 光 傳 感 器（ UTM-30LX ） 、 2.5 m 絲 杠 滑 臺（FLS40L100010C7）、步進電機及驅(qū)動器、STC51 單片機、24 V 直流電源和PC 電腦。通過單片機控制滑臺以3 mm/s 帶動激光掃描儀，將絲杠導(dǎo)軌放置在距果樹中心5 m 的圓周上，以120°夾角分3 次對果樹進行掃描，可全面獲取果樹枝桿的形態(tài)信息。此外，在距果樹中心1 m范圍內(nèi)還需將3 個直徑不同的靶標(biāo)球布置在3 個不同高度的位置上，且要求放置各靶標(biāo)球時對掃描不會產(chǎn)生遮擋果樹或相互遮擋問題。

圖2 固定于室內(nèi)地鉗上的銀杏樹 Fig.2 Ginkgo tree fixed on indoor ground clamp

圖3 點云掃描裝置及掃描方法 Fig.3 Point cloud scanning device and scanning method

3 個位置點3 次掃描果樹樹枝點云圖可通過3 個靶標(biāo)球進行拼接，拼接后的原始點云通過半徑法對果樹原始點云進行去噪處理，對去噪后的點云采用基于點云主方向移動圓柱體的方法獲取果樹的骨架點，由于受拼接精度和激光掃描精度的影響，通過點云提取的果樹骨架點截面半徑會存在某些粗大誤差，對誤差較大的數(shù)據(jù)采用人工測量實際果樹的方式進行修正，最終構(gòu)建果樹空間梁單元振動模型（圖4）。

圖4 果樹建模過程 Fig.4 Fruit tree modeling process

2.2 銀杏樹果實分布特性

銀杏樹的果實和樹葉只會生長在腋芽處，本研究對室內(nèi)試驗銀杏樹上的樹葉和果實進行統(tǒng)計分析，樹葉、果實及腋芽分布狀態(tài)如圖5 所示。

圖5 樹葉、果實及腋芽局部分布圖 Fig.5 Distribution of leaves, fruits and axillary buds

果樹上共計335 個腋芽，其中長果實的腋芽共有242個，約占腋芽總數(shù)的72.2%。由于構(gòu)建有果的果樹模型需要獲得單位長度樹枝上的果實的平均質(zhì)量與平均數(shù)量，因此統(tǒng)計內(nèi)容包括相鄰長果的腋芽分布間距、同一長果腋芽處果實生長的數(shù)量和同一腋芽上果實的平均質(zhì)量。并對統(tǒng)計出的各參量分布規(guī)律用正態(tài)分布進行擬合，擬合公式如式（14）所示

式中a 和b 為擬合常系數(shù)，u 為分布期望值即均值，σ 為分布標(biāo)準(zhǔn)差。對應(yīng)各參數(shù)分布如圖6 所示。

生長果實的腋芽間距分布如圖6a 所示，生長果實的腋芽間距主要集中在3～6 cm 之間，占比為81.4%，整體的分布規(guī)律為均值4.24 cm，標(biāo)準(zhǔn)差1.93 cm 的正態(tài)分布。由圖6b 可知，腋芽處生長的果實數(shù)量主要集中在4～10個之間，其中生長有4、6 和8 個果實的腋芽數(shù)量明顯高于其他果實個數(shù)的腋芽數(shù)量，其原因是絕大多數(shù)銀杏果在一個果柄上會并列生長出2 個果實，所以果實數(shù)量出現(xiàn)雙數(shù)的概率明顯大于奇數(shù)，因此總體分布規(guī)律擬合曲線相決定系數(shù)較低，只有0.757，單個腋芽處生長果實數(shù)量的分布為均值5.69 個，標(biāo)準(zhǔn)差為4.84 個的正態(tài)分布。由圖6c 可知，生長果實的腋芽主要集中在直徑為3～10 mm 的樹枝上，占比為93.0%，總體分布規(guī)律為均值為6.09 mm，標(biāo)準(zhǔn)差4.09 mm 的正態(tài)分布。由圖6d 可以發(fā)現(xiàn)腋芽處的果實質(zhì)量主要集中在5.9～6.5 g 之間，占比為81.4%，其分布為均值6.15 g，標(biāo)準(zhǔn)差0.45 g 的正態(tài)分布。除了腋芽處生長的果實數(shù)量擬合正態(tài)曲線的相關(guān)系數(shù)較低，其余3 種情況擬合正態(tài)曲線時的決定系數(shù)均高達0.922 以上，說明各分布規(guī)律均較好地符合正態(tài)分布。

圖6 果實分布規(guī)律 Fig.6 Fruits distribution law

2.3 樹葉分布特性

對試驗銀杏樹上的樹葉分布進行研究，共計317 個腋芽處生長樹葉，約占腋芽總數(shù)的94.6%。對生長樹葉腋芽的腋芽相鄰間距、同一腋芽處生長樹葉的數(shù)量、同一腋芽處生長樹葉的平均重量和生長樹葉腋芽處的樹枝直徑進行統(tǒng)計，各參數(shù)分布如圖7 所示。

圖7 樹葉分布規(guī)律 Fig.7 Leaves distribution law

由圖7a 可知，生長樹葉的腋芽間距主要分布在2～6 cm 之間，占比為86.1%，腋芽間距分布規(guī)律為均值4.04 cm，標(biāo)準(zhǔn)差2.26 cm 的正態(tài)分布。由圖7b 可知，腋芽處的樹葉數(shù)量主要集中在3～7 片之間，數(shù)量分布最多的為5 片，樹葉數(shù)量分布符合均值為4.85，標(biāo)準(zhǔn)差為2.84的正態(tài)分布。由圖7c 可以看出，能夠生長出樹葉的腋芽主要分布在直徑2～10 mm 的樹枝上，占比為85.8%，有極少數(shù)的直徑為2 mm 以下、10 mm 以上的樹枝也存在少量的樹葉，其分布總體符合均值4.95 mm，標(biāo)準(zhǔn)差4.18 mm的正態(tài)分布。根據(jù)圖7d 所示，葉片的平均質(zhì)量主要分布在0.3～0.6 g，占比為84.9%，其分布符合均值0.40 g，標(biāo)準(zhǔn)差0.28 g 的正態(tài)分布?？傮w來說，腋芽處樹枝直徑按正態(tài)分布曲線進行擬合時的決定系數(shù)最低，但也達到0.852，其他3 種參數(shù)擬合時的決定系數(shù)均高達0.969 以上，說明具有較好的正態(tài)分布規(guī)律。

2.4 有果有葉模型

為驗證上述果實與樹葉分布規(guī)律的準(zhǔn)確性。本研究將上述對果實、樹葉所獲得的各分布規(guī)律，與在南京林業(yè)大學(xué)校園內(nèi)生長的10 a 以上的5 棵大型銀杏樹上隨機選取的共計150 個腋芽處的樹葉與果實分布規(guī)律進行對比，各項分布中誤差最大為樹葉分布中的腋芽處樹枝直徑，相對誤差為7.6%，其主要原因可能為室外果樹具有較多生長有樹葉的末端細枝，導(dǎo)致室外果樹分布均值偏小。其余各分布的均值相對誤差均<5%，對比結(jié)果均滿足以上的各自正態(tài)分布規(guī)律，說明采用這些分布規(guī)律對果樹進行構(gòu)建有果有葉的動力學(xué)模型是可行的。

應(yīng)用上述獲取的銀杏樹果實和樹葉的分布規(guī)律，對銀杏樹枝桿模型按以下規(guī)律添加果實質(zhì)量：生長果實的腋芽間距為4.2 cm，并結(jié)合實際果實生長雙數(shù)居多的現(xiàn)象，單個腋芽生長果實數(shù)量為6 個，單個腋芽上生長果實的平均質(zhì)量為6 g，則對應(yīng)1 cm 長度的樹枝上果實數(shù)量密度ng=1.43 個，應(yīng)用在果實質(zhì)量矩陣中的單個果實平均質(zhì)量mg=6×10-3kg。依據(jù)生長果實的腋芽處樹枝直徑分布可知，實際果實大多生長在直徑3～10 mm 的樹枝上，因此不能對果樹整體枝干上進行附加果實質(zhì)量，而只能在直徑符合特定條件的枝干上附加果實質(zhì)量。因此只在銀杏樹樹枝模型中直徑≤10 mm 的單元，應(yīng)用式（4）計算附加果實質(zhì)量矩陣Mg。

基于樹葉的分布規(guī)律，在銀杏樹枝桿模型上按以下規(guī)律添加樹葉質(zhì)量：生長樹葉的腋芽間距為4 cm，單個腋芽生長葉片數(shù)量為5 個，單個腋芽上生長樹葉的平均質(zhì)量為0.4 g。則按1 cm 長度樹枝上樹葉密度ny=1.32 個，應(yīng)用在樹葉質(zhì)量矩陣中的單個樹葉平均質(zhì)量my=4×10-4kg，且只對枝桿模型中直徑≤12 mm 的單元，應(yīng)用式（5）計算附加樹葉質(zhì)量矩陣My。

通過果實和樹葉的分布特征可以看出，銀杏的樹葉質(zhì)量明顯小于果實的質(zhì)量，對銀杏樹進行建模時，可進一步簡化即忽略銀杏樹的樹葉質(zhì)量。但本研究為全面體現(xiàn)有果有葉的建模方法，并不忽略樹葉的質(zhì)量，依然根據(jù)有果有葉的質(zhì)量矩陣進行計算。對其他樹種來說，特別是果實與樹葉質(zhì)量均對果樹質(zhì)量分布有明顯影響的樹種，在進行果實建模時同時附加果實和樹葉質(zhì)量矩陣是十分有必要的。

應(yīng)用式（6）即可計算得到銀杏樹有果有葉振動模型的質(zhì)量矩陣，結(jié)合式（13）即可計算銀杏樹有果有葉模型的固有頻率。

3 銀杏樹的頻譜特性試驗與仿真分析

銀杏樹頻譜特性測試方法采用沖擊力錘法，具體測試方法為在果樹主干上布置一個加速度傳感器測點（測點0），在6 個主要側(cè)枝上布置6 個加速度傳感器測點（測點1～6）。為減小傳感器質(zhì)量對測試結(jié)果的影響，本研究選取的側(cè)枝測點位置均盡可能距離主干較近即具有較好的剛性，同時又能測試出反映該側(cè)枝的頻譜特性，具體布點如圖2 所示。通過在測點0 處施加沖擊信號，測得0～6 測點的加速度響應(yīng)信號，即可由式（15）計算果樹各個測點的的傳遞函數(shù)。

式中H(f)為果樹的傳遞函數(shù)；Y(f)為加速度信號的傅里葉變換；F(f)為力錘沖擊信號的傅里葉變換。

頻譜測試裝置包括沖擊力錘（LC-02A）、三向壓電式加速度傳感器（CA-YD-141）、電荷放大器（YE5853A）、數(shù)據(jù)采集系統(tǒng)（NI cDAQ-9174）和測試分析軟件（CRAS V7.1）。振動采收應(yīng)用最多的為偏心塊式激振方式，偏心塊式激振的工作頻率區(qū)間主要為15～25 Hz[35]，因此本研究對果樹頻譜特性的頻率關(guān)注范圍設(shè)定在0～30 Hz。

3.1 銀杏樹的頻譜特性試驗分析

通過測試及數(shù)據(jù)處理可獲取各測點傳遞函數(shù)對應(yīng)的頻譜曲線，如圖8 所示為有果有葉頻譜曲線，由于果樹存在較多固有頻率，圖中僅標(biāo)識較為明顯波峰對應(yīng)的頻率即各固有頻率。各測點頻譜曲線中的基頻均為0.39 Hz，在0～5 Hz 內(nèi)存在大量的低矮密集波峰，即果樹在低頻階段聚集大量的固有頻率，這與De Langre[36]的研究結(jié)果是一致的，即在低頻區(qū)出現(xiàn)大量峰點是由于植物的材料柔性大，導(dǎo)致存在多數(shù)固有頻率集中在10 Hz 以下。在5 Hz 以下的各峰點單位力引起的加速度幅值基本處于1 mm/(s2·N)左右且峰寬較窄；高于5 Hz后頻譜曲線的波峰開始變得稀疏、峰寬變寬，即在有果有葉的情況下果樹的高頻振動加速度幅值變大且對應(yīng)阻尼也增大。

圖8 三種狀態(tài)下的果樹頻譜曲線 Fig.8 Spectrum of fruit tree in three states

去掉果實后，各測點頻譜曲線較去除果實前有了一定變化。果樹在無果有葉狀態(tài)下的基頻為 0.46 Hz，比有果實時的基頻略大，且在0～3 Hz 內(nèi)密集波峰基本消失，而在5～10 Hz 內(nèi)由原平緩上升狀態(tài)變?yōu)槌霈F(xiàn)多峰點。各測點頻譜曲線在10～30 Hz 的波峰突出較為明顯，但單位沖擊力引起的加速度幅值并無明顯變化趨勢。

在無果無葉狀態(tài)下，果樹的基頻明顯增大為2.03 Hz。在低頻0～5 Hz 內(nèi)波峰較少，但在5～10 Hz 內(nèi)出現(xiàn)了較多小型波峰，在10～30 Hz 出現(xiàn)大量明顯的較高波峰。與有果有葉和無果有葉狀態(tài)下相比較，無果無葉狀態(tài)下的各測點頻響函數(shù)的單位力激發(fā)的加速度幅值顯著增加，這主要是由于去除樹葉后，果樹的阻尼系數(shù)明顯下降。

綜上所述，隨著逐漸去除掉果樹的果實和樹葉，果樹的低頻固有頻率集中出現(xiàn)的區(qū)域由0～5 Hz 向5～10 Hz 區(qū)域逐步變化，這是由于隨著果實和樹葉的去除，果樹的整體質(zhì)量在下降，但樹枝的整體剛度基本沒有變化。從固有頻率的基本定義可知，其值會隨著剛度與質(zhì)量比的增大逐漸增大，這與上述測試結(jié)果相一致。此外，各測點頻譜曲線幅值在去掉果實和樹葉后有了明顯的提升，尤其在去掉樹葉后增加特別明顯，說明樹葉是影響果樹振動幅值的主要因素，即樹葉的阻尼效應(yīng)明顯大于果樹果實和自身材料的阻尼效應(yīng)。說明樹葉和果實的質(zhì)量分布對果樹的頻譜特性分布有明顯影響，在研究林果采收時不能忽略果實與樹葉附加質(zhì)量的影響，因此構(gòu)建有果有葉的振動模型才能有效反映實際果樹。

3.2 仿真結(jié)果分析

基于第2 節(jié)構(gòu)建的銀杏樹有果有葉有限元模型，利用MATLAB2017 計算該模型的固有頻率。為確定計算模型的所需果樹的密度、彈性模量和剪切模量參數(shù)，本研究通過在校園內(nèi)其他銀杏樹上進行隨機取樣，將樣本木材加工成27 個300 mm×20 mm×20 mm 的木條樣本，測量木條質(zhì)量并計算得銀杏樹的密度為1 250 kg/m3。使用電子萬能試驗機（CMT6104）測試材料的彈性模量，取27 個樣本的平均彈性模量4 713.9 MPa 作為有限元分析模型的彈性模量，取泊松比為0.3，計算可知剪切模量為1 813.0 MPa。由于仿真計算結(jié)果和試驗測試結(jié)果所得固有頻率較多，為方便比較試驗結(jié)果與仿真結(jié)果，通過2個頻率區(qū)域進行驗證模型的準(zhǔn)確性：1）將低頻區(qū)域前15階固有頻率的計算值與實測值進行比對，是為驗證基頻及低階頻率區(qū)域所構(gòu)建模型的仿真計算準(zhǔn)確性；2）在中高頻15～25 Hz 區(qū)域內(nèi)將兩者進行比對，主要是驗證在常用振動采收頻率范圍15～25 Hz 內(nèi)的仿真模型的計算準(zhǔn)確性。

針對果樹有限元模型計算所得的固有頻率出現(xiàn)成對現(xiàn)象，即連續(xù)出現(xiàn)相鄰2 個固有頻率值十分接近，這是由于仿真模型采用圓柱梁單元，而對稱的圓柱梁會在同一固有頻率或稱模態(tài)頻率下沿軸向截面2 個正交方向上產(chǎn)生2 個相互垂直的彎曲振型，但由于實際樹枝上存在著各種小側(cè)枝且樹枝本身很難軸向?qū)ΨQ分布，從而會出現(xiàn)一對對的非常接近的固有頻率值。因此，本研究將成對出現(xiàn)的固有頻率以較高值作為同一階固有頻率。一般仿真計算的固有頻率要多于實測的固有頻率，將實測頻率與仿真計算頻率比對，與仿真計算頻率接近的實測頻率作為同一階頻率。仿真果樹模型密度與實際密度存在的差異會導(dǎo)致質(zhì)量矩陣的誤差，枝桿模型的彈性模量與剪切模量與實際的差異會導(dǎo)致剛度矩陣的誤差，因此實際計算出的固有頻率存在一定誤差，Siller 等[29]采用修改模型的彈性模型與等效剛度使基頻誤差從20%相對誤差下降到1%以內(nèi)，本研究采用類似方法將實際樹木密度和彈性模量分別提升了15%和降低了10%，3 種模型的原基頻計算精度及調(diào)整后的計算精度如表1 所示。仿真模型經(jīng)過調(diào)整后，無果無葉狀態(tài)下從最大的相對誤差14.8%，明顯降低為0.5%，且其余2 種狀態(tài)下的仿真基頻與實測基頻完全吻合，即相對誤差為0。

表1 仿真基頻與調(diào)整仿真基頻表 Table 1 Simulated fundamental frequency and adjusted simulated fundamental frequency

將調(diào)整后的前15 階固有頻率的仿真與實測結(jié)果列入表2 中，將15～25 Hz 范圍內(nèi)固有頻率的仿真與實測結(jié)果列于表3 中。實測中的主要固有頻率已在圖8 中表示出來。

由表2 可知，無果無葉果樹模型仿真前15 階固有頻率主要集中在10 Hz 以下，基本呈現(xiàn)均布分隔狀態(tài)；與實測結(jié)果相比對，實測中未出現(xiàn)第2、5、6 階頻率，其余12 階仿真值與實測值非常接近；仿真計算最大相對誤差為6.67%，平均相對誤差為2.03%，且其中4.69、5.54、7.81 Hz 在實測頻譜曲線中幾乎為各個測點均出現(xiàn)的共同固有頻率點。4.69、5.54 和7.81 Hz 振型圖如圖9 所示。可以看出這3 個振型圖中大部分側(cè)枝均有較明顯的振幅，但不同頻率下各枝的振型與振動幅值不同。無果有葉狀態(tài)下，各階仿真頻率與無果無葉模型相比出現(xiàn)明顯下降，致使各階頻率間隔明顯縮小，即前15 階固有頻率集中在3 Hz 以下；仿真前15 階固有頻率中只有9 階固有頻率在實測中出現(xiàn)，仿真計算最大誤差只有3.85%。有果有葉狀態(tài)下，前15 階仿真固有頻率進一步降低，密集處于1.12 Hz 以下，絕大多數(shù)仿真頻率在實測頻譜曲線中未能出現(xiàn)，只出現(xiàn)了5 階固有頻率，因較低頻率區(qū)內(nèi)激振力引起的加速度響應(yīng)靈敏度較低，有些固有頻率較難激發(fā)起來，即在頻譜曲線中難以體現(xiàn)出來，能在實測中體現(xiàn)出來的各階固有頻率最大相對誤差為2.75%，平均相對誤差為1.36%。

表2 前15 階仿真與實測固有頻率 Table 2 Top 15 steps of simulated and measured natural frequency

圖9 部分振型圖 Fig.9 Partial mode shape diagram

表3 15～25 Hz 仿真與實測固有頻率 Table 3 15-25 Hz simulated and measured natural frequency

由表3 可知，在15～25 Hz 頻率區(qū)間內(nèi)，無果無葉狀態(tài)下共有7 階仿真固有頻率，實測頻率對應(yīng)出現(xiàn)了5 階；無果有葉與有果有葉狀態(tài)下，均有10 階仿真頻率，實測頻率分別對應(yīng)出現(xiàn)了7 階與6 階，均比前15 階頻率中對應(yīng)程度高得多。3 種仿真模型的最大相對誤差分別為5.76%、2.06%和1.98%，平均相對誤差分別為2.32%、0.82%和0.95%。

綜上所述，仿真固有頻率多于實測頻譜曲線體現(xiàn)出來的固有頻率，這是因為實測的頻譜曲線上某些部位的曲線波峰非常微弱，只使曲線出現(xiàn)微弱起伏，難以判斷此處是否存在固有頻率，即并未將其標(biāo)注出來，致使實測固有頻率數(shù)會少于仿真頻率數(shù)。因此，在處理仿真頻率與實測頻率對應(yīng)關(guān)系時，是將實測固有頻率與仿真頻率中非常接近的固有頻率值作為對應(yīng)的同等階次。無果無葉果樹的各階固有頻率分布較稀疏，而隨著樹葉與果實的添加，各階固有頻率明顯越來越密集，并在實測中所能體現(xiàn)出來的對應(yīng)固有頻率越來越少。同時仿真計算的固有頻率只要在實測中能體現(xiàn)出來，其相對誤差均較低，最大不超過6.67%，說明所構(gòu)建的有果有葉模型是可行的，具有一定的實用價值。

3.3 計算誤差分析

從試驗結(jié)果中可以看出，對果樹附加質(zhì)量后的計算精度相較無果無葉模型的計算精度提高，可應(yīng)用以下相關(guān)計算公式進行分析論述。為方便理解，只取其中某一樹枝單元基于集中質(zhì)量方程進行論述，假設(shè)構(gòu)建模型時給定該單元樹枝的質(zhì)量和剛度分別為mz和kz，而實際質(zhì)量和剛度分別為mz+Δmz和kz+Δkz，其中Δmz和Δkz分別為樹枝的質(zhì)量誤差和剛度誤差且Δmz和Δkz分別遠小于mz和kz，假設(shè)實際樹枝的無誤差精確振動方程如式（16）所示

式中mz為樹枝質(zhì)量，kg；kz為樹枝剛度，N/m；Δmz為樹枝的質(zhì)量誤差，kg；Δkz為樹枝的剛度誤差，N/m；Δmz和Δkz分別遠小于mz和kz；x 表示位移，m；x˙˙表示加速度，m/s2。

則建模型計算獲得的固有頻率誤差Δf 如式（17）所示

式中Δf 為固有頻率誤差，Hz。

樹體上附加果實和樹葉質(zhì)量后并不會對其剛度產(chǎn)生影響，因此式（17）中僅考慮質(zhì)量誤差Δm。對式（17）進行簡化處理并略去高階無窮小項后得式（18）

附加樹葉和果實質(zhì)量后的計算誤差分別如式（19）和式（20）所示

式中Δ?z+y為無果有葉固有頻率計算誤差，Hz；Δ?z+y+g為有果有葉固有頻率計算誤差，Hz；my為樹葉質(zhì)量，kg；mg為果實質(zhì)量，kg；Δmy為樹葉質(zhì)量誤差，kg；Δmg為果實質(zhì)量誤差；且Δmy和Δmg分別且分別遠小于my和mg。

假設(shè)樹葉質(zhì)量和果實質(zhì)量的相對誤差與樹枝質(zhì)量相對誤差處于同一數(shù)量級，則由式（18）～式（20）可知，附加質(zhì)量后的分母上的項 mz+my+mg與分子上的項Δmz+Δmy+Δmg的冪次數(shù)分別為1.5 和1，誤差中的分母項的增加值遠大于分子項的增加值，因此固有頻率計算誤差隨著附加樹葉質(zhì)量和果實質(zhì)量減小。綜上所述，果樹模型在剛度不變的情況下，隨著附加較準(zhǔn)確的樹葉質(zhì)量與果實質(zhì)量，可以減小單純樹枝建模不準(zhǔn)確所帶來的計算誤差。

4 結(jié) 論

本研究基于對無果無葉果樹激光掃描進行三維重建，建立有果有葉果樹的仿真計算模型，并通過試驗對模型進行驗證。具體結(jié)論為：

1）將激光點云與有限元分析相結(jié)合，通過提取果樹枝桿骨架點和樹枝半徑，將樹枝看作空間圓柱梁，通過分別構(gòu)建果樹的質(zhì)量矩陣、剛度矩陣與阻尼矩陣，結(jié)合果樹果實和樹葉的質(zhì)量矩陣，構(gòu)建有果有葉果樹空間6自由度圓柱梁振動模型。

2）提出了構(gòu)建樹葉、果實與腋芽沿果枝生長方向的正態(tài)分布形態(tài)，找到了銀杏樹果實和/或樹葉的腋芽數(shù)量與腋芽間距、果實或樹葉數(shù)量、樹枝直徑、果實或樹葉平均質(zhì)量的分布規(guī)律，確定了銀杏樹空間振動模型中果實與樹葉的分布質(zhì)量及質(zhì)量矩陣。且果實與樹葉質(zhì)量分布矩陣的構(gòu)建方法可應(yīng)用于其他品種果樹有果有葉振動模型的構(gòu)建。

3）對銀杏樹有果有葉、無果有葉和無果無葉3 種狀態(tài)進行頻譜特性測試，發(fā)現(xiàn)樹葉和果實的質(zhì)量會明顯影響果樹的固有頻率數(shù)量及其分布。說明只構(gòu)建果樹的枝干模型確定果樹的頻譜特性，無法應(yīng)用于果樹實際果實基于頻譜的共振采收。

4）通過對銀杏樹具體的仿真建模和試驗獲得的各階固有頻率進行比對驗證，仿真計算固有頻率數(shù)量要多于實測所體現(xiàn)出來的固有頻率數(shù)量，但實測固有頻率均可在仿真頻率中找到十分相近的值，且仿真模型的果樹各枝振型幅值與實測頻譜曲線中各枝的加速度響應(yīng)幅值之間的對應(yīng)關(guān)系基本吻合。有果有葉仿真模型在常用采收頻率15～25 Hz 內(nèi)最大相對誤差均<6.67%，并對仿真模型的計算誤差進行了分析說明，證明本研究所述果樹建模方法的可行性。

本研究所構(gòu)建的果樹振動模型，可通過激光掃描設(shè)備應(yīng)用于田間果樹固有頻率的快速計算，從而可確定田間采收設(shè)備的實際采收頻率，此為本研究的后續(xù)進展。