基于傅里葉描述子和共振頻率的不規(guī)則果形梨果硬度評(píng)估指標(biāo)構(gòu)建

王 鵬 ，吳 杰 ,2※

（1.石河子大學(xué)機(jī)械電氣工程學(xué)院，石河子 832000；2.綠洲特色經(jīng)濟(jì)作物生產(chǎn)機(jī)械化教育部工程研究中心，石河子 832003）

0 引 言

果蔬硬度是其重要的內(nèi)部品質(zhì)之一，硬度與果蔬的口感、成熟度等密切相關(guān)，是梨果品質(zhì)檢測(cè)重要指標(biāo)。國(guó)內(nèi)外學(xué)者研究發(fā)現(xiàn)共振頻率是評(píng)估硬度變化的重要參數(shù)。通常果蔬硬度采用頻率指標(biāo)f2m2/3（f為赤道部感測(cè)頻率，m為質(zhì)量）評(píng)價(jià)[1]，但已有研究表明，果形對(duì)果蔬共振頻率有影響，Chen等[2]研究指出，評(píng)價(jià)指標(biāo)f2m2/3適于近球形果蔬的硬度評(píng)價(jià)，但對(duì)非球形果蔬硬度評(píng)價(jià)誤差較大。為此，Cherng等[3]提出了一種適于評(píng)價(jià)橢球形果蔬硬度的雙頻指標(biāo)SCherng=()2/3m2/3ρ1/3（f1為赤道部感測(cè)頻率，f2為萼端或梗端感測(cè)頻率，ρ為果蔬密度），筆者團(tuán)隊(duì)針對(duì)多為卵圓形的庫爾勒香梨，從香梨赤道部和萼端分別感測(cè)共振頻率f1和f2，構(gòu)建了一種雙頻硬度指標(biāo)Sf1f2= ( 0.4f1+ 0.6f2)2m2/3，提高了香梨硬度評(píng)價(jià)性能[4]。盡管如此，梨果種類豐富，果形復(fù)雜，其中啤梨為典型的果形不規(guī)則梨果，現(xiàn)有的各種基于規(guī)則果形構(gòu)建的頻率指標(biāo)，可能都不適于啤梨硬度的精確評(píng)價(jià)。因此，有必要基于共振頻率構(gòu)建一種適于評(píng)價(jià)果形不規(guī)則梨果硬度的指標(biāo)。

針對(duì)復(fù)雜果形果蔬，了解果形對(duì)共振頻率的影響規(guī)律，是構(gòu)建消除果形影響的新硬度評(píng)價(jià)指標(biāo)的關(guān)鍵。大多研究采用縱橫比、圓度比、圓形度、偏心度、離心率等描述果形[5-9]，李凡等[5]采用縱橫比和圓度比對(duì)香梨形狀定量描述獲得5種果形，申翠香等[7]采用圓形度描述蘋果果形。這些描述方法多用于蘋果、桃子、柑橘等“比較端正”果形的描述[8-11]，但卻很難全面描述果形復(fù)雜不規(guī)則的梨果。近年來，郝敏等[12]采用具有平移伸縮和旋轉(zhuǎn)不變性且有多特征描述能力的Zernike矩描述馬鈴薯，但Zernike矩描述果形時(shí)，數(shù)據(jù)處理量大，計(jì)算耗時(shí)長(zhǎng)[13]。傅立葉描述子具有很強(qiáng)的果蔬外形重建功能且對(duì)復(fù)雜果形有較強(qiáng)的描述能力[14-16]，這一方法在描述木瓜、畸形馬鈴薯、楊桃等果蔬外形輪廓時(shí)，能夠兼顧全局及局部特征[17-19]，其中 Limsiroratana等[17]研究木瓜果形描述方法時(shí)，發(fā)現(xiàn)傅里葉描述子計(jì)算快速且抗噪性強(qiáng)；EIMasry等[18]采用傅里葉描述子描述畸形馬鈴薯，并選取4個(gè)傅里葉描述子識(shí)別畸形馬鈴薯，識(shí)別率可達(dá) 96.2%；Abdullah等[19]采用傅里葉描述子描述楊桃外形，對(duì)3種變形果的識(shí)別率最高可達(dá) 100%。上述研究表明傅里葉描述子可以描述果形復(fù)雜多變的果蔬外形特征。

因此，本文采用傅里葉描述子對(duì)啤梨果形進(jìn)行定量描述，并結(jié)合試驗(yàn)?zāi)B(tài)和有限元模態(tài)分析方法，了解啤梨果形變化對(duì)其頻率的影響規(guī)律，構(gòu)建適于不規(guī)則果形的梨果硬度準(zhǔn)確評(píng)價(jià)指標(biāo)，為復(fù)雜果形梨果硬度振動(dòng)法無損檢測(cè)提供研究基礎(chǔ)。

1 材料與方法

1.1 啤梨試樣

啤梨試樣于2019年10月于石河子水果批發(fā)市場(chǎng)采購，剔除損傷及病蟲害試樣，立即貯藏于?2～0 ℃、相對(duì)濕度85%～95%的冷庫備用。試驗(yàn)前將試樣在(23±1)℃室溫下回溫24 h。

1.2 研究思路

首先結(jié)合試驗(yàn)?zāi)B(tài)和有限元模態(tài)分析結(jié)果確定啤梨試樣適宜的振型和對(duì)應(yīng)的模態(tài)頻率，采用傅里葉描述子和縱橫比描述啤梨試樣果形；然后利用Pearson線性相關(guān)性分析果形對(duì)啤梨振動(dòng)頻率的影響規(guī)律，進(jìn)而構(gòu)建含有果形描述子和共振頻率的新硬度評(píng)估指標(biāo)以消除果形對(duì)硬度評(píng)估的影響；最后結(jié)合M-T穿刺硬度實(shí)測(cè)值，采用普通最小二乘法（ordinary least square, OLS）分別對(duì)基于頻率的無損硬度與M-T穿刺硬度SMT進(jìn)行線性回歸分析。研究方案流程如圖1所示。

1.3 啤梨試驗(yàn)?zāi)B(tài)分析

試驗(yàn)?zāi)B(tài)分析系統(tǒng)如圖2所示，選擇單點(diǎn)感測(cè)多點(diǎn)激勵(lì)方式采集振動(dòng)信號(hào)構(gòu)建頻率響應(yīng)函數(shù)[4]，啤梨試樣用橡皮筋懸掛使其處于接近自由振動(dòng)狀態(tài)。根據(jù)前期試驗(yàn)在試樣表面均布56個(gè)點(diǎn)，取赤道部的32號(hào)點(diǎn)（見圖3）為感測(cè)點(diǎn)，LC1305型力錘（河北朗斯有限責(zé)任公司）以9～12 N的力從1號(hào)點(diǎn)依次開始敲擊各點(diǎn)[4]，固結(jié)在感測(cè)點(diǎn)的LC0408T型加速度傳感器（河北朗斯有限責(zé)任公司）感測(cè)信號(hào)，隨機(jī)選取10個(gè)試樣進(jìn)行試驗(yàn)?zāi)B(tài)分析。

1.4 啤梨有限元模態(tài)分析

如圖4所示，采用Rigelscan Elite型三維掃描儀（精度0.0185 mm，武漢中觀自動(dòng)化科技有限公司）獲得啤梨試樣.asc格式點(diǎn)云模型，然后將其導(dǎo)入正逆向混合設(shè)計(jì)軟件Geomagic Wrap 2017（美國(guó)，Geomagic公司）進(jìn)行逆向處理，獲得.stp格式實(shí)體模型，最后將實(shí)體模型導(dǎo)入有限元分析軟件Abaqus 2017（法國(guó)，達(dá)索SIMULIA公司）中進(jìn)行有限元分析。

采用有限元軟件Abaqus 2017（法國(guó)，達(dá)索SIMULIA公司）劃分網(wǎng)格時(shí)，由于啤梨試樣小且形狀不規(guī)則，因此選用10節(jié)點(diǎn)二次四面體結(jié)構(gòu)實(shí)體單元C3D10進(jìn)行自由網(wǎng)格劃分，避免大量分割操作導(dǎo)致工作量增加和人為誤差。為減少計(jì)算時(shí)間，在保證精度的前提下將網(wǎng)格大小設(shè)定為2 mm。假設(shè)啤梨為均質(zhì)、各向同性的線彈性體，模型泊松比取0.3[20]，采用排水法測(cè)得啤梨密度。采用有限元分析軟件 Abaqus進(jìn)行啤梨有限元模態(tài)分析時(shí)，在Frequency分析步中選擇子空間迭代算法（subspace）求解器，求解運(yùn)動(dòng)方程+Ka(t) = 0（M為質(zhì)量矩陣，為加速度矢量，K為剛度矩陣，a(t)為速度矢量），設(shè)置最大頻率為1 000 Hz，頻移為100 Hz，特征值提取數(shù)量為20，每次循環(huán)矢量數(shù)為28次，最大循環(huán)次數(shù)為30次，獲得啤梨前20階模態(tài)振型和對(duì)應(yīng)的固有頻率進(jìn)行分析。

1.5 啤梨固有頻率歸一化

由于質(zhì)量對(duì)模態(tài)固有頻率影響很大[21-22]，采用模態(tài)分析啤梨果形對(duì)頻率影響時(shí)，用式（1）對(duì)各階模態(tài)固有頻率進(jìn)行歸一化處理以消除質(zhì)量對(duì)固有頻率的影響[20]：

式中m為試樣質(zhì)量，kg；m0為常數(shù)，根據(jù)前期試驗(yàn)取0.1 kg；f為某一階模態(tài)固有頻率，Hz；fn為歸一化的固有頻率，Hz。

1.6 啤梨果形描述方法

1.6.1 果形輪廓的提取

采用EOS 750D型佳能數(shù)碼相機(jī)（日本，佳能有限公司）采集啤梨圖像，采用Matlab 2014（美國(guó)，MathWorks公司）軟件中圖像處理工具箱對(duì)啤梨圖像進(jìn)行灰度變換、閾值分割、邊界提取，獲得試樣外形輪廓，如圖5所示。

1.6.2 試樣輪廓半徑序列的計(jì)算

如圖6，計(jì)算試樣輪廓曲線的形心坐標(biāo)x0和y0[23]：

式中n為試樣邊界輪廓曲線上的總樣點(diǎn)數(shù)；k為邊界輪廓曲線上點(diǎn)的序數(shù)。求得形心坐標(biāo)x0和y0后，以形心坐標(biāo)為圓心，果梗與果體交界處P1點(diǎn)作為起始點(diǎn)，按逆時(shí)針方向求半徑序列r(k)：

1.6.3 半徑序列的傅里葉變換

為避免試樣大小對(duì)果形分析的影響，將半徑序列r(k)歸一化處理[24]：

式中rav為試樣的平均半徑，Pixels；n為試樣邊界輪廓的樣點(diǎn)數(shù)，rg(k)為歸一化半徑序列。對(duì)歸一化半徑序列rg(k)進(jìn)行離散傅里葉變換：

式中Fh為傅里葉描述子，h為傅里葉變換后的自變量，2π/n為相鄰兩樣點(diǎn)間夾角，F(xiàn)h可理解為半徑序列變化h次的程度[23]。當(dāng)h=0時(shí)，F(xiàn)0表征試樣平均半徑信息，F(xiàn)1表征試樣外形的彎曲程度信息，F(xiàn)2表征試樣外形的伸長(zhǎng)度信息，F(xiàn)3表征試樣外形的三角度信息，F(xiàn)4表征試樣外形的方形度信息[24]。

1.6.4 縱橫比q的計(jì)算

采用縱橫比描述啤梨果形時(shí)，測(cè)量如圖7所示啤梨試樣肩高的最大值Hmax(mm)和最小值Mmin(mm)以及赤道部最大直徑Dmax(mm)和最小直徑Dmin(mm)，縱橫比q由公式（9）計(jì)算求得[4]。

1.7 啤梨硬度的M-T穿刺法測(cè)定

采用TA.XP plus質(zhì)構(gòu)儀（英國(guó)，Stable Micro System公司）在啤梨梗端、赤道部和萼端 3個(gè)區(qū)域內(nèi)每個(gè)區(qū)域隨機(jī)選取4個(gè)點(diǎn)進(jìn)行帶皮M-T穿刺試驗(yàn)[25-26]。采用直徑5 mm 的圓柱探頭，穿刺速度為 1 mm/s，穿刺深度為8 mm。穿刺硬度由穿刺前曲線的斜率計(jì)算得出，然后取12次穿刺試驗(yàn)平均值作為梨整體硬度。

1.8 硬度評(píng)估指標(biāo)的評(píng)價(jià)方法

在室溫(23±1) ℃、相對(duì)濕度 20%環(huán)境下，隨機(jī)選取140個(gè)啤梨試樣等分為7組，每2 d進(jìn)行一次測(cè)試試驗(yàn)。試驗(yàn)時(shí)隨機(jī)選取一組啤梨試樣，每個(gè)試樣在15 min內(nèi)依次進(jìn)行激振試驗(yàn)、圖像采集和M-T穿刺試驗(yàn)，獲得啤梨基于頻率的無損硬度和M-T穿刺硬度。

采用普通最小二乘法（Ordinary Least Square, OLS）分別對(duì)基于頻率的無損硬度與M-T穿刺硬度SMT（N/mm）進(jìn)行線性回歸分析[26]，回歸方程為

式中SZ為各頻率指標(biāo)中的某一個(gè)；α和β分別為回歸方程的回歸系數(shù)和截距。

2 結(jié)果與分析

2.1 試驗(yàn)?zāi)B(tài)對(duì)有限元模態(tài)的驗(yàn)證

2.1.1 試驗(yàn)?zāi)B(tài)與有限元模態(tài)的分析結(jié)果

瑞士有很多傳統(tǒng)的奶酪飲食，其中奶酪火鍋就是其中最著名的一款。如果說瑞士有什么美食的話，那一定首推奶酪火鍋。瑞士盛產(chǎn)奶酪，而且法令規(guī)定一個(gè)地區(qū)只能生產(chǎn)該地區(qū)產(chǎn)的奶酪，因此瑞士的奶酪直至今日仍能夠保持極其傳統(tǒng)的風(fēng)味，且不同地區(qū)的奶酪味道各不相同。

通過啤梨試驗(yàn)?zāi)B(tài)分析獲得 3個(gè)模態(tài)，分別為彎曲模態(tài)、擠壓模態(tài)和呼吸模態(tài)，如圖8所示?？梢钥吹?，彎曲模態(tài)在啤梨梗端發(fā)生較大變形，在萼端幾乎不發(fā)生振動(dòng)變形；擠壓模態(tài)在赤道部發(fā)生最大變形，呼吸模態(tài)主要沿啤梨徑向發(fā)生擴(kuò)張和收縮的振動(dòng)變形。

啤梨進(jìn)行有限元模態(tài)分析所獲得的20階有限元模態(tài)振型中，第1～6階頻率近似為0，為剛體模態(tài)；第16～20階振型變形復(fù)雜，為二階模態(tài)，本研究不予考慮。在第7～15階振型中，第8階與第7階變形相同，均為彎曲模態(tài)，第10階為擠壓模態(tài)，第13～15階模態(tài)變形相同，均為呼吸模態(tài)。

2.1.2 試驗(yàn)?zāi)B(tài)與有限元模態(tài)頻率對(duì)比

調(diào)整有限元分析彈性模量值，當(dāng)?shù)?7階有限元彎曲模態(tài)頻率與試驗(yàn)?zāi)B(tài)結(jié)果相同時(shí)，確定彈性模量值Eest，取這一彈性模量值為有限元模態(tài)分析的彈性模量并計(jì)算剩余各階有限元模態(tài)的固有頻率值[4]。計(jì)算有限元模態(tài)頻率與試驗(yàn)?zāi)B(tài)頻率的差值，以差值絕對(duì)值與試驗(yàn)?zāi)B(tài)頻率比值大小判別有限元模態(tài)分析的準(zhǔn)確性。任取10個(gè)啤梨試樣進(jìn)行試驗(yàn)，結(jié)果如表1所示。頻率的有限元結(jié)果與試驗(yàn)結(jié)果的相對(duì)誤差不超過3.25%,表明采用各向同質(zhì)線彈性體的梨實(shí)體建模進(jìn)行有限元模擬結(jié)果較準(zhǔn)確，因此可以采用有限元模態(tài)分析果形對(duì)頻率的影響規(guī)律。

表1 試驗(yàn)?zāi)B(tài)分析與有限元模態(tài)分析頻率比較Table 1 Comparisons between the experimental frequencies and the FEM frequencies

2.2 傅里葉描述子對(duì)啤梨外形描述結(jié)果分析

如圖9所示，前20個(gè)傅里葉描述子與F0的比值隨h增大迅速減小，但F4和F6與F0的比值出現(xiàn)反彈。當(dāng)h大于16時(shí)，F(xiàn)h與F0比值曲線趨于水平且比值小于0.005，即啤梨果形信息主要集中在前16個(gè)傅里葉描述子。

梨試樣輪廓和前10、15、20個(gè)傅里葉描述子反變換圖如圖10所示，采用前10個(gè)傅里葉描述子反變換獲得的梨輪廓圖基本能描述梨果形，但對(duì)啤梨梗端或萼端輪廓重合度較差；采用前15個(gè)傅里葉描述子反變換獲得的輪廓與原輪廓形狀度吻合度高。因此，盡管采用前20個(gè)傅里葉描述子反變換獲得的輪廓仍然與梨輪廓保持很好重合度，但采用前15個(gè)傅里葉描述子已足以描述啤梨試樣輪廓的主要信息。

2.3 果形對(duì)啤梨共振頻率的影響

對(duì) 100個(gè)啤梨進(jìn)行有限元模態(tài)分析，本研究取彎曲模態(tài)、呼吸模態(tài)和擠壓模態(tài)的歸一化固有頻率進(jìn)行分析。采用Pearson相關(guān)性分析方法對(duì)啤梨彎曲、呼吸和擠壓3種模態(tài)的歸一化固有頻率和前15個(gè)傅里葉描述子進(jìn)行線性相關(guān)性分析，與歸一化固有頻率有較強(qiáng)相關(guān)性的 7個(gè)傅里葉描述子的相關(guān)系數(shù)如表2所示。

表2 歸一化固有頻率和傅里葉描述子間的相關(guān)性Table 2 Correlation between the normalized frequencies and the Fourier descriptors

由圖11可知，歸一化固有頻率fn與果形系數(shù)間關(guān)系可用下式表示

式中fn為歸一化固有頻率，S為果形描述子，a和b的取值由振型決定，各階振型對(duì)應(yīng)a和b取值見圖11。

2.4 啤梨硬度評(píng)估指標(biāo)的構(gòu)建

果蔬硬度可用彈性模量E表征，且頻率f的平方與彈性模量E成正比[27-28]，故下式成立：

經(jīng)轉(zhuǎn)換可得

將公式（1）和（11）代入公式（13），得

式中m0為常數(shù)，a、b值由圖11獲得，EFEM為有限元模擬彈性模量值，取2.1.2節(jié)計(jì)算Eest的平均值6.86 MPa，代入公式（14）中得到硬度評(píng)價(jià)指標(biāo)SE：

圖12為啤梨在赤道部感測(cè)的振動(dòng)頻譜曲線，共振頻率fBE為彎曲模態(tài)的頻率，其對(duì)應(yīng)加速度幅值不突出試驗(yàn)時(shí)不易提取。擠壓模態(tài)的共振頻率fC的加速度幅值很高，試驗(yàn)時(shí)非常容易提取。因此，在應(yīng)用新的硬度評(píng)價(jià)指標(biāo)檢測(cè)啤梨時(shí)，評(píng)價(jià)指標(biāo)中a、b分別取圖11中擠壓模態(tài)對(duì)應(yīng)的值。

2.5 各頻率指標(biāo)對(duì)啤梨硬度的評(píng)價(jià)能力比較

啤梨在室溫條件下硬度隨儲(chǔ)藏時(shí)間變化規(guī)律如圖13所示，由圖可知啤梨M-T穿刺硬度在12 d內(nèi)總體呈現(xiàn)先快后慢的下降趨勢(shì)，M-T穿刺硬度在前6 d內(nèi)急劇下降，6～12 d之間硬度下降速度減緩。

將縱橫比q和F1代入公式（15）獲得含有果形描述子的硬度指標(biāo)SE(q)和SE（F1），由1.8節(jié)方法對(duì)基于頻率的硬度指標(biāo)與M-T穿刺硬度SMT進(jìn)行線性回歸分析，結(jié)果如圖14所示。傳統(tǒng)單頻硬度評(píng)估指標(biāo)Sf與 M-T穿刺硬度SMT的決定系數(shù)較低（R2=0.733），啤梨試樣不規(guī)則果形可能是導(dǎo)致單頻指標(biāo)Sf硬度評(píng)估值與M-T穿刺硬度實(shí)測(cè)值相關(guān)性不高的原因。

Cherng等[3]提出的雙頻硬度評(píng)價(jià)指標(biāo)SCherng和筆者團(tuán)隊(duì)研究提出的雙頻指標(biāo)Sf1f2與 M-T穿刺硬度SMT的決定系數(shù)R2非常接近，分別為0.775和0.765，較傳統(tǒng)單頻硬度評(píng)估指標(biāo)而言，與硬度實(shí)測(cè)值SMT相關(guān)系數(shù)有所提高，但由于這 2種硬度指標(biāo)都只能消除近橢形果形對(duì)頻率的影響，相關(guān)系數(shù)提高效果并不顯著。采用新硬度指標(biāo)SE評(píng)估啤梨硬度時(shí)，以縱橫比q和傅里葉描述子F1作為啤梨果形描述子分別帶入新指標(biāo)SE計(jì)算獲得的指標(biāo)SE(q)與M-T穿刺硬度SMT的決定系數(shù)R2值為0.746，指標(biāo)SE（F1）與 M-T穿刺硬度SMT的決定系數(shù)R2值為 0.892，在各指標(biāo)中最高。這表明傅里葉描述子能準(zhǔn)確描述不規(guī)則果形梨果，本研究提出的硬度評(píng)價(jià)指標(biāo)SE（F1）更適于復(fù)雜且不規(guī)則果形的啤梨硬度評(píng)估。

3 結(jié) 論

1）啤梨振動(dòng)試驗(yàn)?zāi)B(tài)分析獲得彎曲模態(tài)、擠壓模態(tài)和呼吸模態(tài)，彎曲模態(tài)和呼吸模態(tài)在赤道部振動(dòng)變形較小，擠壓模態(tài)振動(dòng)變形在赤道部變形較大，采用加速度傳感器在赤道部感測(cè)時(shí)更易于獲得擠壓模態(tài)的共振頻率。

2）采用前 15個(gè)傅里葉描述子足以描述啤梨外形特征，第1個(gè)傅里葉描述子F1與啤梨彎曲模態(tài)、擠壓模態(tài)和呼吸模態(tài)歸一化固有頻率高度線性相關(guān)，線性相關(guān)系數(shù)r分別為?0.923、?0.922、0.700，較適于描述啤梨復(fù)雜果形的變化。

3）基于擠壓模態(tài)的共振頻率（fC）和傅里葉描述子F1構(gòu)建硬度指標(biāo)SE(F1)，可消除果形對(duì)啤梨硬度評(píng)估的影響。與其他研究報(bào)道的頻率評(píng)估指標(biāo)相比，該硬度指標(biāo)檢測(cè)啤梨硬度值與M-T穿刺法檢測(cè)硬度值的決定系數(shù)R2為0.892，可用于果形復(fù)雜梨果硬度的較精確評(píng)估。