基于AR-MED的闊葉材原木聲信號特征參數提取及原木質量分等

瞿玉瑩， 楊 揚， 徐 鋒

(南京林業大學 信息科學技術學院，南京 210037)

闊葉材原木的質量隨樹種、產地等變化較大?？斩?、節子、外傷或其他生長缺陷降低了原木的強度并影響其外觀，因此降低了其產品的價值[1-2]。通常一塊較好的單板原木，其價格是工廠普通一級鋸木價格的1.5～6倍[3]，如未能及時檢測到木材內部隱藏的缺陷，往往會給木材制造商或企業造成重大的經濟損失。闊葉材原木早期缺陷檢測的另一重要目的是從待加工的原木流中去除無利潤收益的原木(加工成本大于收益)。理想上，為實現利潤的最大化，這些原木應出售給更低加工成本的企業，制作成如枕木、木托盤、紙漿和燃料等產品。

聲波技術在木材檢測領域的研究已有幾十年歷史。研究表明聲學分析能夠被用來預測原木結構材的硬度和強度[4-6]。這些研究為聲學技術應用于針葉材原木的質量分等和分級開辟了道路。目前，一些商用聲學儀器已被部分林產品企業用于針葉材原木的在線質量控制和現場分等[7-9]，但由于闊葉材缺陷種類繁多，聲波法在其內部缺陷探測的研究還非常有限。Wang等[10-11]應用聲共振法測試了幾種闊葉材原木并基于聲速實現了粗略地原木質量分等，研究指出聲速與原木板材的彈性模量、等級和缺陷等存在相關性。盡管聲速簡單分等已被一些試驗證明是有效的方法，但聲速與原木性質之間的直接相關性并不顯著，其決定系數一般位于0.38～0.54[12]，原木分等遠沒有達到期望的精度。遺憾的是除了聲速，目前尚未發現其他聲參數應用于闊葉材原木內部缺陷的檢測或板材質量的評估[13]。因此，進一步挖據數據信息，提取聲信號特征參數是聲波法應用于闊葉材原木質量評估或分等的重要任務之一。

峭度(kurtosis)是反映信號分布特性的數值統計量，對信號方差和幅值的變化非常敏感，現已被逐漸運用到動態信號的特征提取中[14]。傳輸于健全木材中的聲波，信號的幅值近似于正態分布，峭度值波動較小，約為某一常數；當原木內部存在缺陷時(如，空洞、腐朽、開裂等)，因疊加了缺陷信號分量，信號幅值分布發生偏離，峭度值隨之改變。因而從峭度值的變化可以判定是否有隨機成分被引入到信號中。原木聲振動信號由于傳播路徑、木材自身變異性、缺陷和噪聲的影響，特征較為復雜且有很強的隨機性，很多情況下若直接計算峭度值進行缺陷判斷，效果并不理想。

自回歸(Autoregressive, AR)模型是時序分析方法中最基本、實際應用最廣泛的時序模型。研究表明，AR模型的自回歸參數對狀態變化規律反映最為敏感[15-16]。Wang等[17]建立了一種AR模型線性濾波器，用于分離裂紋故障振動信號，并將預測誤差信號的峭度值作為評判閾值診斷齒輪故障。田艷等[18]對塔機在沖擊載荷作用下的狀態位移數據建立AR模型，以方差作為特征參數對塔機損傷狀態進行分析，研究指出方差的大小跟損傷程度呈正相關。但AR模型對信號相位不敏感，因此難以區分噪聲和突變信號。

最小熵反褶積(Minimum Entropy Deconvolution, MED)方法最早由Wiggins[19]于1978年首次提出，并成功應用于地震波處理中，是一種不需要任何先驗假設的信號時域盲解卷積技術[20]。該方法可以提高被分析信號的信噪比，增強信號中的突變成分，使得采集信號盡可能還原損傷或缺陷特征成分，非常適用于早期的缺陷診斷。

因此，本文提出一種基于AR-MED的聲信號特征參數提取方法。首先對采集的聲信號進行AR模型預測，濾除確定性信號，減小其對缺陷信號的影響；其次對包含缺陷信號和噪聲的殘差信號進行MED去噪并增強信號突變成分；然后計算解卷后的聲信號峭度，并依據峭度進行原木質量粗略分等；最后與傳統的聲速分等方法進行比較，論證本文方法的可行性和有效性。

1 峭 度

峭度K又稱峰態系數，表征概率密度分布曲線在平均值處峰值高低的特征數，是歸一化的4階中心矩。其定義為

(1)

2 AR模型

AR模型是用自身作回歸變量的過程，即利用前期若干時刻的隨機變量的線性組合來描述以后某時刻隨機變量的線性回歸模型。

設{xi}(i=1,2,…,N)是長度為n的零均值平穩信號序列，p是AR模型的階次，則xk的AR模型可以表示為

(2)

AR模型的階數選擇是一個很重要的問題，階數越低，估計譜越平滑、譜峰越少，越能反映信號的總體情況，但也造成較大的預測偏差；階數越高，估計譜波動越劇烈、譜峰越多，越能突出信號的細節信息，但會導致虛假譜峰和統計的不穩定性[21]。本文使用赤池信息量準則(Akaike Information Criterion, AIC)來確定AR模型的最優階數。AIC建立在熵的概念上，提供了權衡估計模型復雜度和擬合數據優良性的標準。其定義為

(3)

3 最小熵解卷積

最小熵解卷積 (MED) 是一種自適應濾波器設計方法，反作用信號的傳遞路徑，以恢復輸入信號的突變成分。最小熵解卷積濾波器的參數是通過對信號峭度取最大化而獲得，根據Wiggins的定義，信號的峭度最大，則其熵最小，所以稱這種方法為最小熵解卷積。

設濾除平穩成分的原木聲振動信號為εk(k=1,2,…,N)，則由信號與系統理論可知，其可表示為

εk=(gk+nk)*hk

(4)

式中：gk為缺陷信號；nk為混入的背景噪聲；hk為傳輸路徑的沖擊響應函數；*為卷積運算。MED的基本目的就是求解最優逆濾波器fk系數解集以消除傳輸路徑和噪聲的影響，使濾波輸出信號yk不斷逼近缺陷信號。

MED逆濾波器的實現利用了目標函數法(Objective Function Method, OFM)，即通過不斷修正濾波器的系數來優化目標函數，實現MED輸出信號的峭度值達到最大。當系數在設定的容限范圍內收斂時，優化過程結束。MED逆濾波器的實現過程如下：

步驟1設逆濾波器fk的長度為L，輸出為yk，輸入為εk，則

(5)

步驟2求濾波器輸出信號的峭度值

(6)

步驟3為求解滿足式(6)取最大值時的濾波器系數，分別對式(6)兩邊求關于濾波器f的1階偏導數，并令其為零

(7)

步驟4式(5)兩邊對濾波器f求偏導，得

(8)

步驟5將式(5)和式(8)代入式(7)，整理可得

(9)

式(9)可寫成矩陣形式

b=Af

(10)

式中：列向量b為濾波器f的輸入與輸出信號的互相關加權和；A(L×L)為輸入信號的Toeplitz自相關矩陣；f為所求逆濾波器的系數向量。

步驟6求解濾波器f系數。因式(10)的高度非線性，直接求解很難進行，因此可通過迭代計算實現最優解的求取。迭代算法步驟歸納如下：

步驟a計算自相關矩陣A；

步驟b初始化逆濾波器系數f(0)，一般可設其為單位時延脈沖濾波器，即f(0)=[0,…,0,1,0,…,0]T；

步驟c由輸入信號εk和初始濾波器系數f(0)按式(5)計算輸出信號y(0)；

步驟d由輸入信號εk和輸出信號y(0)按式(10)計算列向量b(1)；

步驟e求解新濾波器系數f(1)=A-1b(1)；

步驟f計算迭代終止條件

(11)

當計算誤差的期望值E(err)小于設定的容限時，迭代結束；否則更新逆濾波器的系數f(0)=f(1)并返回步驟c，重新下一輪的迭代循環。

4 基于AR-MED的特征參數提取

當原木截面受到單位縱向沖擊力作用時，產生的聲振動信號反應了原木內部的質量狀態。若原木處于均勻連續的無缺陷狀態，則由一維桿振動理論可知其脈沖響應為無窮多個固有頻率主振型的疊加。而當原木內部存在缺陷時(如空洞、節子、腐朽和裂紋等)，缺陷所導致的一些突變、畸變和衰減等成分必然使振動響應信號發生偏斜或調制等狀態。因此，實際采集的振動信號xk可以描述成原木共振信號dk(確定性信號)、缺陷信號gk和背景噪聲nk與傳輸路徑響應函數hk的卷積，即

xk=[dk+gk+nk]*hk

(12)

AR-MED過程的目的就是要從檢測信號xk中分離出缺陷信號gk，并提取聲特征參數用于闊葉材原木質量的分等。圖1所示為AR-MED過程，具體步驟如下：

步驟1按式(2)AR模型預測聲信號確定性成分，并將其從原信號中減去，得到包含缺陷成分和噪聲的殘差信號；

步驟2將殘差信號作為MED濾波器的輸入，經多次迭代運算使濾波輸出信號逼近缺陷信號成分；

步驟3計算輸出信號的峭度值，并將其作為聲特征參數用于原木質量分等。

圖1 AR-MED解卷積過程Fig.1 AR-MED filtering deconvolution process

5 算例仿真

為驗證AR-MED算法增強突變信號的有效性，并考慮原木聲信號的成分多樣性、瞬時性以及噪聲的干擾，構建如下信號并加入噪聲予以仿真

(13)

式中：n(t)為高斯白噪聲；a=0.05為阻尼系數；參數f1=260 Hz，f2=510 Hz，f3=3 kHz。信號x(t)的波形如圖2(a)所示，采樣頻率為20 kHz，采樣點數為2 048點。

為說明AR-MED對信號突變成分的增強能力，分別計算信號x(t)、殘差ε(t)和輸出信號y(t)(如圖2所示)的峭度。由圖2可知，原信號的峭度值為2.89，經AR濾波器去除周期平穩成分后的殘差信號(包括突變成分和噪聲)峭度值上升到3.26。圖2(b)為AIC熵值與AR模型階次的對應關系。由AIC準則可知，階數90為AR預測濾波器實現平穩信號與突變成分分離的最佳階數。但由于AR模型基于最小相位假設，相應的預測譜中只包含了信號的幅度信息，不包含相位信息，因此，不能很好地分離突變成分與噪聲，即殘差信號中既包含了突變成分也包括了噪聲，因此峭度值提升并不明顯。圖2(d)為MED濾波器的輸出信號。經濾波器進一步壓制噪聲，輸出信號越加逼近突變成分，峭度值達到了4.70，較初始信號峭度值有較明顯的提高，實現了信號突變成分的增強。

通過仿真信號的分析可知，本文提出的方法能有效地增強信號的突變成分，提高其峭度值，因此能夠用作原木的缺陷識別，并對其質量進行分等。

圖2 仿真信號及AR-MED濾波Fig.2 Simulated signal and AR-MED filtering process

6 材料與實驗

6.1 樣品材料

測試樣本為美國西弗吉尼亞州阿巴拉契亞中部山區Mead West Vaco森林的黃楊樹木(Liriodendron tulipifera)，隨機選取15棵伐倒。每棵樹被鋸截成3～5根商用長度的原木，共52根。觀察發現這些原木的質量等級差異較大，一些原木鋸截后存在較明顯的腐爛痕跡，一些原木則有較深的外部損傷，另外也有部分高質量原木。測試實驗在砍伐后一個月，時間為三月初。原木樣品的含水量為45%～60%，均在纖維飽和點(30%)以上。因此，本文討論的所有聲參數均指生材參數。

6.2 實驗過程

6.2.1 原木基本參數測量

為了獲取原木的基本物理參數，首先對原木進行高分辨率三維激光掃描。掃描沿原木長度方向進行，采樣步進為1.59 mm，依據原木周長差異，每周像素點數為250～450點，平均分辨率為3 pixel/cm。圖3為原木掃描實驗現場。隨后用量程為2 000 kg，精度為0.5 kg的吊秤(LHS4000a, ADAME quipment, Inc., Oxford,CT)對每根原木質量進行測定。最后根據原木的總質量和三維激光掃描數據中獲得的精確體積計算原木密度。所有上述參數均列于表1。

圖3 運行中的原木激光掃描系統Fig.3 Laser-based field scanning system for log

6.2.2 聲沖擊試驗

掃描完成后，隨即進行聲沖擊試驗。為采集沖擊響應聲信號，將一對傳感器探針(Fakopp spike sensor, Fakopp Enterprise Bt., Agfalva, Hungary)分別插入到原木兩端靠近中心的等高位置。沖擊信號由5.44 kg的大錘敲擊原木一端產生，響應信號由連接于電腦的數據采集卡(NI 5132)采集，采樣頻率為20 kHz，采樣長度為1 000點。試驗在溫度為21 ℃，相對濕度為50%的條件下進行。聲信號采集示意圖如圖4所示。

圖4 應力波信號采集示意圖Fig.4 Stress wave signal sampling schematic diagram

6.2.3 原木鋸切與板材分等

聲學試驗完畢，基于視覺檢查和基本測試參量，從52根原木中挑選出能代表此批次質量的21根原木進行鋸切，對比驗證本文所提方法的可行性和有效性(之所以僅選擇21根原木，主要為經費與時間所限)。按最大得材率和產值方案將這些原木鋸切成厚度為29 mm的板材，并依據美國闊葉材板材協會(the National Hardwood Lumber Association, NHLA)分等規則(NHLA 2015)[22]對這些板材進行分等。板材等級依次分為高等級、普一級、普二級、普三級和等外級。板材總體積、各等級板材體積和等級率列于表2。

表1 黃楊原木尺寸和物理參數

表2 黃楊原木鋸切結果

7 結果與分析

7.1 峭度提取

對采集的聲信號按節4方法分別進行峭度計算?，F以1123原木為例簡單重述峭度提取過程，其余信號的峭度值直接列于表3，不再贅述。

圖5(a)所示為1123原木時域信號，其波形變化規則，似無雜波，對其進行傅里葉變換，頻譜較為清晰，因此很難確定有無缺陷存在。但從板材鋸切結果來看，原木中存在較大的活節等缺陷。然而由于缺陷信號較弱，完全被主振動信號所淹沒。對圖5(a)信號進行AR預測濾波，獲取包含缺陷成分的殘差信號。AR模型階數與AIC值分布關系及殘差信號分別如圖5(c)和圖5(d)所示。為壓制噪聲，增強缺陷信號成分，進一步MED處理，輸出缺陷信號如圖5(e)所示，計算其峭度值，得13.8。圖5(f)為缺陷信號頻譜，顯然與主振動平穩成分實現了分離，突變分量得到了增強。

表3 黃楊原木聲參數計算值

圖5 1123原木信號的AR-MED缺陷成分提取Fig.5 AR-MED based defect components extraction of acoustic signal from 1123 log

7.2 基于峭度的原木質量分等

對上述獲得的峭度參數，按其大小進行排序，并將其對應的原木分成三個等級(高質量、中等和低質量)。因樣本數量較小，分等閾值的設置原則主要基于參數區間的三等分，再結合樣本數目進行適當微調，確保每個等級的原木在5～9棵。為驗證分等效果，對每個質量等級的原木，根據其切開的板材總體積和各等級板材體積，計算各板材等級率(各等級板材體積/板材總體積×100%)。理想上，高質量原木其高等級板材率應最高，而低質量原木高等級板材率應最低，中等質量原木應介于兩者之間。

傳輸于無缺陷原木中的聲信號，其振動波形規則，幅度分布接近于正態分布，峭度值大小中等。當原木內部存在開裂或活節(群)等缺陷時，信號中較大幅值出現的概率增加，幅值分布偏離正態分布，峭度值隨之增大，且缺陷越嚴重，峭度值越大；另，若缺陷以腐朽或空洞為主導時，信號中較小幅值將顯著增多，其峭度值將減小[23]。因此，基于峭度的原木質量分等具體如下：

G1: 6.5≤K<9.5 (高質量等級)；

G2: 9.5≤K<12.5(中等質量等級)；

G3: 12.5≤K或K<6.5(低質量等級)。

圖6為各等級板材等級率隨峭度的分布。顯然三個質量等級組中板材等級率與峭度值之間存在較顯著的相關性。隨著峭度從G1增加到G3，對應組所含高等級板材的等級率顯著下降(從G1組的72.5%下降到G2組的38.9%再到G3組的21.8%)，而普一、普二、普三等級板材的等級率明顯上升，其總和占比由G1組的26.2%上升至G2組的59.2%到G3組的77.2%。等級率隨峭度值相反的變化趨勢表明，缺陷少、質量好的原木和有缺陷、質量一般的原木對信號傳輸的影響有明顯的差異。另，三組原木中所有等外級(below grade)板材之和約占4.2%，這意味著此批次原木質量相對較好，含有較大缺陷的原木較少。

圖6 各板材等級率隨峭度分組的分布關系Fig.6 Distribution relationship of the board grade yield as the log kurtosis classes

為進一步證實本文方法的有效性，對這些樣本使用聲速進行質量分等。聲速的測量方法參見Xu等的研究，其測量值列于表3。類似峭度分等，也將聲速進行排序并分成三組，各板材等級率與速度關系如圖7所示。

由圖7可知，中等質量原木組中高等級板材率最高，約53.1%，低質量原木組中次之，為47.7%，高于高質量組中的43.9%。速度分組似乎與板材的等級率沒有明確的預測關系，而僅在等外級才顯示出較明顯的相關性，即隨著三組原木聲速的增大，等外級板材率從3%下降到0.6%。從物理上來說，傳播于原木中的聲速大小依賴于原木的強度和密度。原木內部存在的小缺陷雖并不影響整個原木的強度，但根據NHLA板材分級規則，它們卻顯著影響板材的等級。此外，由于聲波的繞射作用，沿著原木長度方向傳播的縱波對一些小的缺陷并不敏感，因此不會影響波的傳播路徑。這些可能是導致聲速與板材等級不相關的原因。反之，若原木內部存在較嚴重缺陷，而此類缺陷會明顯減小木材密度或影響波的傳輸行為時，則會導致原木異常的低聲速?；谏鲜龇治?，5個原木(1147，1154，1167，1128和1129)被確定為具有異常聲速(2.91～3.05 km/s)。檢查鋸切結果，證實這些原木或存在較大的內部腐朽(1128，1129，1154和1167)或存在致使鋸切板材降等的結構缺陷(1129)。因此，聲速仍然能夠有效識別具有嚴重內部缺陷的原木，此結論與文獻[24]一致。

圖7 各板材等級率隨聲速分組的分布關系Fig.7 Distribution relationship of the board grade yield as the log acoustic velocity classes

8 結 論

針對闊葉材原木缺陷聲信號的非平穩性和特征重疊性，提出了基于AR-MED相結合的聲參數(峭度)提取方法，并對原木質量實現分等，由對應板材等級率進行分等結果評估。

(1)AR-MED方法能夠實現缺陷成分與主振動信號的有效分離，并壓制噪聲，使缺陷成分得到顯著增強，為峭度作為分等參數提供明確的依據。

(2)峭度值與闊葉材原木內部缺陷狀態有較顯著的相關性，基于板材等級率分布而言，峭度能夠實現原木質量較準確的分等。

(3)聲速與闊葉材板材等級率并沒有顯示明晰的相關性，即它只能對具有嚴重缺陷的原木進行較準確的質量評估，而對于存在小缺陷但不影響縱波傳輸特性的原木似乎無能為力。

由于闊葉材內部缺陷的聲檢測與評估國內外尚處于初步研究階段，其檢測方法和評估參數的挖掘還待進一步探索，尤其是聲評估參數的定量指標和普適性還有待于大量實驗數據的積累和歸納。

致謝

本項目是南京林業大學與美國農業部林產品實驗室、明尼蘇達德魯斯大學合作研究項目(14-JV-11111133-089)。美國林產品實驗室的王喜平研究員參與了項目的討論，并給出了指導意見。本實驗樣品的鋸切和板材分等工作由美國農業部林業科學實驗室的工程師Neal Bennett和Deborah Conner完成。沒有他們的參與和幫助，項目研究將很難進行。在此，作者對上述研究人員和工程師表示由衷地感謝。