高禮彬，陳金浩，張怡卓，王克奇

(東北林業(yè)大學(xué)機(jī)電工程學(xué)院，哈爾濱 150040)

木材順紋抗壓強(qiáng)度是結(jié)構(gòu)用材重要的力學(xué)指標(biāo)之一，而傳統(tǒng)的破壞性檢測(cè)方法時(shí)間長(zhǎng)、成本高，且由于木材的各向異性及不均勻性，造成操作復(fù)雜且準(zhǔn)確度低。近紅外光譜技術(shù)(NIR)具有簡(jiǎn)單便捷、快速有效的特點(diǎn)，通過建立NIR定量分析模型，可實(shí)現(xiàn)快速無損測(cè)定[1]。黃安民等[2]指出近紅外光譜技術(shù)在木材物理力學(xué)研究方面具有極大的潛力；虞華強(qiáng)等[3]將近紅外光譜作為初步檢測(cè)杉木木材力學(xué)性質(zhì)的工具；趙榮軍等[4]利用近紅外光譜方法預(yù)測(cè)人工林粗皮桉木材力學(xué)性質(zhì)效果較好。因此，近紅外光譜技術(shù)在木材力學(xué)性能檢測(cè)中被證明是一種理想可靠的方法。

在近紅外光譜建模過程中，原始光譜會(huì)包含大量無用信息，要先進(jìn)行特征選擇，實(shí)現(xiàn)降維處理。云永歡等[5]和Yun等[6]提出模型集群分析(MPA)算法框架運(yùn)用于光譜特征優(yōu)選，打破了傳統(tǒng)一次性建模思路，力求最大限度地利用已有樣本集的信息，通過隨機(jī)采樣，從大量的子模型中提取統(tǒng)計(jì)信息并進(jìn)行選擇分析?；贛PA算法框架開發(fā)的算法常見的有：變量組合集群分析(VCPA)[7-8]算法，其根據(jù)優(yōu)秀子集中變量出現(xiàn)頻率作為變量重要性的評(píng)價(jià)標(biāo)準(zhǔn)選取特征變量組合；競(jìng)爭(zhēng)性自適應(yīng)重加權(quán)(CARS)[9-10]算法，其通過搜索方法及評(píng)價(jià)變量重要性程度的標(biāo)準(zhǔn)以得到最優(yōu)的變量子集；變量迭代空間收縮(VISSA)[11-12]算法，其在變量空間中迭代收縮，保證新的變量空間優(yōu)于前一個(gè)，得到最優(yōu)變量組合。其中，VCPA和CARS的弊端在于當(dāng)校正集樣本出現(xiàn)波動(dòng)時(shí)，對(duì)于變量的重要性判斷會(huì)受到一定程度的影響，會(huì)降低模型的預(yù)測(cè)能力。VISSA避免了這種缺點(diǎn)，但因其在每輪迭代子數(shù)據(jù)集時(shí)，得到的是特征波長(zhǎng)數(shù)量比例不同的子數(shù)據(jù)集，可能無法準(zhǔn)確評(píng)價(jià)子模型的優(yōu)劣，導(dǎo)致有效信息丟失，增加不穩(wěn)定性。因而引入采樣誤差分布分析(SEPA)[13-14]與VISSA相結(jié)合進(jìn)行改進(jìn)，相比于原VISSA能以一種更溫和的方式消除非信息變量。

合適的建模方法是建立光譜數(shù)據(jù)與真值之間聯(lián)系的關(guān)鍵，由于木材的各向異性及不均勻性，使高維度特征變量之間的相關(guān)性非常敏感，采用線性建模會(huì)導(dǎo)致模型的精度大大降低，如偏最小二乘回歸法(PLS)[15-16]。而非線性建模能有效提升預(yù)測(cè)精度，有較好的泛化能力且不會(huì)出現(xiàn)過擬合的問題，如支持向量回歸(SVR)[17-18]，但其不提供后驗(yàn)概率且參數(shù)不易設(shè)置。因此，在SVR基礎(chǔ)上提出基于貝葉斯框架的稀疏概率學(xué)習(xí)模型——相關(guān)向量機(jī)(RVM)[19]，其不僅具有SVR模型的優(yōu)點(diǎn)，且核函數(shù)可以任意指定，不是必須正定，此外調(diào)節(jié)參數(shù)較少，便于模型優(yōu)化，相較于SVR模型能夠完成更精確的預(yù)測(cè)。

筆者以樺木木材為試驗(yàn)對(duì)象，使用近紅外光譜儀采集光譜數(shù)據(jù)并測(cè)定順紋抗壓強(qiáng)度真值，然后以SEPA作為VISSA的改進(jìn)策略進(jìn)行特征波長(zhǎng)篩選，在此基礎(chǔ)上采用RVM非線性建模方法進(jìn)行木材抗壓強(qiáng)度預(yù)測(cè)，并與其他波長(zhǎng)優(yōu)選算法和建模方法效果進(jìn)行比較。

1 材料與方法

1.1 試驗(yàn)材料

樺木木材淡褐色至紅褐色，可用作地板、家具、內(nèi)部裝飾材料、車船設(shè)備、膠合板等。所制家具光滑耐磨，花紋明晰，如今多用于結(jié)構(gòu)、鑲花木細(xì)工和內(nèi)部框架的制作。試驗(yàn)樺木樣本來自黑龍江帶嶺林業(yè)局林場(chǎng)，取5株樺樹樣木，樹齡均在20 a以上，在每株樣木的胸高(1.3 m)附近截取約50 cm的圓盤后去皮，氣干后鋸解成力學(xué)試材毛坯條。參照GB/T 1935—2009《木材順紋抗壓強(qiáng)度試驗(yàn)方法》，制取尺寸為30 mm×20 mm×20 mm的順紋方向抗壓力學(xué)試樣，并挑選出無疵試樣100條，按1～100編號(hào)，將試樣放入恒溫恒濕箱內(nèi)調(diào)至含水率為12%，保持室溫為(20±2)℃、相對(duì)濕度為(65±3)%，然后進(jìn)行近紅外光譜掃描和力學(xué)性能測(cè)試。

1.2 近紅外光譜采集與預(yù)處理

研究表明，木材近紅外光譜在1 000～1 600 nm波長(zhǎng)范圍內(nèi)攜有重要信息，能夠較好地預(yù)測(cè)木材力學(xué)強(qiáng)度等性質(zhì)[3,16]。因此，采用INSION公司近紅外光纖光譜儀進(jìn)行光譜掃描，采集波長(zhǎng)范圍900～1 700 nm，環(huán)境溫度20 ℃，平均相對(duì)濕度50%。使用兩分叉光纖探頭及鎢光源，光纖探頭垂直固定在金屬支架上，試樣置于支架底面，非接觸測(cè)量距1 mm，光斑直徑5 mm，采用聚四氟乙烯白板標(biāo)定后進(jìn)行光譜采集。在每個(gè)試件的徑、弦切面上分別均勻掃描采集4個(gè)樣本點(diǎn)，每點(diǎn)掃描30次自動(dòng)平均為1個(gè)光譜，記錄保存。將每個(gè)試件的徑、弦切面共8條光譜數(shù)據(jù)平均成1條光譜，代表該編號(hào)試件整體的近紅外吸收光譜。

本試驗(yàn)利用INSION公司開發(fā)的SPEC view 7.1 軟件進(jìn)行數(shù)據(jù)采集，由于原始光譜存在基線漂移、樣本顆粒大小不一和光散射等噪聲信息影響因素而混入非信息變量，因此需對(duì)所采集數(shù)據(jù)做預(yù)處理后再建模。分別采用多元散射校正(MSC)、卷積平滑(SG)和MSC-SG這3種方法對(duì)原始光譜數(shù)據(jù)進(jìn)行預(yù)處理，然后針對(duì)各處理后的全波段光譜建立PLS模型。SG-PLS模型效果最優(yōu)，校正集和預(yù)測(cè)集的決定系數(shù)分別為0.908 8和0.881 4，相應(yīng)的均方根誤差分別為3.765 3和3.086 5，如表1所示。由表1可知，經(jīng)過SG預(yù)處理后能夠?qū)υ脊庾V信息進(jìn)行去噪，所建立的抗壓強(qiáng)度PLS模型更加穩(wěn)定和精確。

表1 不同預(yù)處理方法的PLS模型對(duì)比Table 1 PLS model comparisons of different pretreatment methods

1.3 真值測(cè)試

參照GB/T 1935—2009測(cè)定樺木的無疵試樣抗壓強(qiáng)度真值。對(duì)100組樺木樣本，以校正集與預(yù)測(cè)集之比為3∶1的比例，采用光譜-理化值共生距離(SPXY)算法[20]對(duì)其進(jìn)行分割，得到校正集樣本75個(gè)，預(yù)測(cè)集樣本25個(gè)，如表2所示。由表2可知，抗壓強(qiáng)度范圍為60.822 1～107.476 2 MPa，預(yù)測(cè)集樣品信息被校正集樣品信息所覆蓋。

表2 樣本校正集和預(yù)測(cè)集抗壓強(qiáng)度的測(cè)定結(jié)果Table 2 Test results of compressive strengths of sample calibration set and prediction set

1.4 近紅外光譜SEPA-VISSA-RVM建模方法

1.4.1 SEPA-VISSA算法

VISSA的核心在于通過加權(quán)二進(jìn)制采樣(WBMS)構(gòu)建子模型，然后提取預(yù)測(cè)誤差最小的子模型，統(tǒng)計(jì)各變量出現(xiàn)的頻率，分配變量權(quán)重來實(shí)現(xiàn)變量空間收縮。其中，WBMS主要根據(jù)變量當(dāng)前權(quán)重來分配子數(shù)據(jù)集中的變量數(shù)，確保權(quán)值較大的變量有較高被選擇概率，使每個(gè)變量有不同的采樣頻率。變量的權(quán)重公式為：

ωi=fi/mbest

(1)

式中：fi為變量在最佳子模型中的頻率；mbest為最佳模型的數(shù)量；ωi為包含變量的子模型比例。

SEPA是將子模型均方根誤差(RMSE，公式中記為RMSE)的中位數(shù)和標(biāo)準(zhǔn)差相加作為RMSE的一種估計(jì)。使用中位數(shù)而不是平均值和最小值，是因?yàn)橹形粩?shù)更穩(wěn)健。在VISSA算法中通過WBMS建立m個(gè)子模型，最大主成分?jǐn)?shù)為n，然后根據(jù)SEPA原理，以其中一個(gè)子模型的n個(gè)交叉驗(yàn)證均方根誤差(RMSECV，公式中記為RMSECV)進(jìn)行計(jì)算，取RMSECV的中位數(shù)和計(jì)算n個(gè)RMSECV的標(biāo)準(zhǔn)差后相加，公式如下：

(2)

式中：RMSECV為子模型的交叉驗(yàn)證均方根誤差；n為RMSECV的個(gè)數(shù)；STD為各子模型RMSE的標(biāo)準(zhǔn)差。

根據(jù)上述原理，設(shè)光譜數(shù)據(jù)為X，抗壓強(qiáng)度真值為Y，特征波長(zhǎng)個(gè)數(shù)為p，WBMS采樣數(shù)為m，則SEPA-VISSA算法過程如下：

步驟1，賦予特征波長(zhǎng)p值均為0.5的初始權(quán)重。

步驟2，采用WBMS對(duì)光譜數(shù)據(jù)X生成m個(gè)子數(shù)據(jù)集，分別建立PLS子模型。

1.4.2 RVM建模原理

RVM是建立在貝葉斯稀疏框架下的稀疏概率模型。與SVR相比，滿足Merrcer核正定的條件不再是必要的，且在計(jì)算負(fù)擔(dān)方面也更小。在波長(zhǎng)優(yōu)選的結(jié)果下，令校正集為{xn,tn|n=1,2,…,N}，建立木材抗壓強(qiáng)度預(yù)測(cè)模型：

(3)

式中：xn為輸入的變量，即光譜的特征波長(zhǎng)；tn為目標(biāo)輸出，即木材抗壓強(qiáng)度；?為權(quán)重；ωm為?中的第m個(gè)元素；M(1≤m≤M)為xn中光譜特征波長(zhǎng)的數(shù)量；K(x,xn)為核函數(shù)；ω0為基礎(chǔ)量；ξn為殘差且滿足ξn～N(0,σ2)。設(shè)tn為獨(dú)立分布，則似然函數(shù)可表示為：

(4)

式中：σ2為噪聲方差；t為目標(biāo)輸出，是以tn組成的向量(t1,t2,…,tn)T；φ為基函數(shù)，是以φn組成的矩陣[φ1,φ2,…,φM]，其中φ(xn)=[1,k(xn,x1),k(xn,x2),…,k(xn,xM)]。

為了確定先驗(yàn)分布，假設(shè)參數(shù)ωn服從以0為均值、σ2為方差的高斯分布，則有：

(5)

式中，α為決定權(quán)值?先驗(yàn)分布的超參數(shù)，其只與相對(duì)應(yīng)的權(quán)值?相關(guān)。假定α和σ2均服從Gamma先驗(yàn)概率，然后根據(jù)貝葉斯準(zhǔn)則可得?的后驗(yàn)分布為：

(6)

高斯(Gaussian)核函數(shù)(式7)、冪指數(shù)(Exponential)核函數(shù)(式8)和拉普拉斯(Laplacian)核函數(shù)(式9)是3種典型的核函數(shù)，具有良好的分布特性。分別采用這3種核函數(shù)，結(jié)合粒子群算法(PSO)[21]選取其最優(yōu)核參數(shù)并在RVM中建模，最后比較優(yōu)劣。各核函數(shù)表達(dá)式如下所示：

(7)

(8)

(9)

式中：yc為核函數(shù)中心；σ為核寬度。

2 結(jié)果與分析

2.1 特征波長(zhǎng)提取

以SG預(yù)處理后的512個(gè)波長(zhǎng)點(diǎn)作為對(duì)象，采用不同算法設(shè)置參數(shù)，并進(jìn)行降維試驗(yàn)，過程如圖1所示。其中，VCPA 的二進(jìn)制采樣(BMS)采樣數(shù)為500，指數(shù)遞減函數(shù)(EDF)運(yùn)行數(shù)為200，比率為0.5；CARS的蒙特卡洛采樣數(shù)為2 000；VISSA和SEPA-VISSA的加權(quán)二進(jìn)制采樣(WBMS)采樣數(shù)均為2 000。試驗(yàn)均設(shè)最大主成分?jǐn)?shù)n為10，交叉驗(yàn)證為5折，每個(gè)試驗(yàn)各進(jìn)行10次，取試驗(yàn)結(jié)果中值。

不同算法的波長(zhǎng)選擇對(duì)比見圖2。結(jié)合圖1和2可知：①VCPA在EDF運(yùn)行到187次時(shí)，最小RMSECV為3.672 2，此時(shí)得到15個(gè)波長(zhǎng)點(diǎn)，占總波長(zhǎng)的3%；②CARS在蒙特卡洛采樣1 799次時(shí)，最小RMSECV為12.597 0，此時(shí)得到12個(gè)波長(zhǎng)點(diǎn)，占總波長(zhǎng)的2%；③VISSA選取了77個(gè)波長(zhǎng)點(diǎn)，占總波長(zhǎng)的15%，此時(shí)模型的最小RMSECV為3.681 8；④SEPA-VISSA選取了111個(gè)波長(zhǎng)點(diǎn)，占總波長(zhǎng)的22%，此時(shí)最小RMSECV±STD(SEPA)為3.822 0。此處4種算法的采樣方法和判斷最優(yōu)變量特征方式都各不相同，由于是隨機(jī)采樣的子模型，僅表現(xiàn)當(dāng)前對(duì)應(yīng)算法選取特征時(shí)的子模型狀況，相互間無直接聯(lián)系。

對(duì)于木材的主要成分纖維素和木質(zhì)素，其含有大量的含氫基團(tuán)C—H、O—H、N—H等在近紅外光譜區(qū)產(chǎn)生吸收，其中二級(jí)倍頻主要位于1 000～1 400 nm處，一級(jí)倍頻主要位于1 400～1 800 nm處，因此使其在近紅外光譜區(qū)域有豐富的吸收信息。從圖2可以看出，在990，1 210，1 460和1 660 nm附近有主要吸收。根據(jù)吸收峰可得，1 210和1 660 nm 處分別為C—H鍵二級(jí)倍頻伸縮和一級(jí)倍頻伸縮，1 460 nm處為O—H鍵一級(jí)倍頻伸縮振動(dòng)，990 nm處為N—H鍵二級(jí)倍頻對(duì)稱伸縮振動(dòng)。由圖2可知，VCPA和CARS選擇波長(zhǎng)數(shù)量較少，可能會(huì)遺漏信息變量而降低建模準(zhǔn)確性，而SEPA-VISSA相比于VISSA所選特征波長(zhǎng)增加了1 460 nm處C—H鍵一級(jí)倍頻振動(dòng)以及一些非主要的吸收峰，這是VISSA算法在特征波長(zhǎng)多而雜時(shí)，子數(shù)據(jù)集隨機(jī)產(chǎn)生所帶來的不穩(wěn)定性，而SEPA-VISSA算法能更穩(wěn)定包含更多信息變量。

圖1 不同算法的波長(zhǎng)選擇過程Fig. 1 Wavelength selection process of different algorithms

圖2 不同算法的波長(zhǎng)選擇對(duì)比Fig. 2 Comparison of wavelength selections of different algorithms

對(duì)各算法優(yōu)選后的特征波長(zhǎng)進(jìn)行PLS建模，如表3所示。其中，SEPA-VISSA在校正集中決定系數(shù)為0.940 2，校正均方根誤差為3.049 2；在預(yù)測(cè)集中，決定系數(shù)為0.959 3，預(yù)測(cè)均方根誤差為2.899 5，相對(duì)分析誤差為3.025 6。相比之下，驗(yàn)證了SEPA-VISSA在校正集和預(yù)測(cè)集中的PLS模型效果均優(yōu)于其他算法，主要原因如下：①相對(duì)于復(fù)雜的抗壓強(qiáng)度本身而言，影響抗壓強(qiáng)度的因素有很多，如木材的各向異性及不均勻性等特點(diǎn)，但通過VCPA和CARS所提取的特征波長(zhǎng)可能會(huì)排除主要信息變量，導(dǎo)致預(yù)測(cè)結(jié)果不理想；②在SEPA-VISSA算法中，運(yùn)用SEPA中的RMSECV±STD作為子模型的評(píng)判指標(biāo)對(duì)VISSA改進(jìn)，提取出的特征波長(zhǎng)信息與抗壓強(qiáng)度實(shí)際值擬合度更高，除了主要的特征波長(zhǎng)，還包含一些非主要的特征波長(zhǎng)，以更加柔和精確的方式剔除無信息特征波長(zhǎng)，進(jìn)一步提升預(yù)測(cè)模型的準(zhǔn)確性。

表3 不同算法的PLS模型對(duì)比Table 3 PLS model comparisons of different algorithms

2.2 建模方法及比較

將SEPA-VISSA特征優(yōu)選后的數(shù)據(jù)作為輸入，以PSO算法分別對(duì)Gaussian、Exponential和Laplacian核函數(shù)的核寬度進(jìn)行優(yōu)化，其中PSO參數(shù)慣性權(quán)重ω和學(xué)習(xí)因子C1、C2為預(yù)置值，設(shè)最大迭代數(shù)為80，上下界為26和2-6，5折交叉驗(yàn)證，各運(yùn)行10次，取最佳結(jié)果。不同核函數(shù)的PSO優(yōu)化過程如圖3所示，3種核函數(shù)參數(shù)優(yōu)化的迭代次數(shù)分別為53,59和38，最佳適應(yīng)度值為3.027 1，3.325 1和3.132 3。不同核函數(shù)的PSO-RVM模型效果對(duì)比見表4，根據(jù)表4，采用Laplacian建立RVM模型時(shí)，其核寬度為10.404 3，預(yù)測(cè)決定系數(shù)為0.944 9，預(yù)測(cè)均方根誤差為2.043 2，均優(yōu)于基于Gaussian和Exponential的RVM模型，因此采用Laplacian進(jìn)行建模，可進(jìn)一步提高抗壓強(qiáng)度預(yù)測(cè)模型的精度。

圖3 不同核函數(shù)的PSO優(yōu)化過程Fig. 3 PSO optimization process of different kernel functions

表4 不同核函數(shù)的PSO-RVM模型效果對(duì)比Table 4 Comparison of PSO-RVM model accuracies of different kernel functions

不同建模方法的模型效果對(duì)比如圖4所示，為SEPA-VISSA提取特征波長(zhǎng)分別采用PLS、SVR和RVM這3種方法進(jìn)行建模。其中PLS模型的預(yù)測(cè)決定系數(shù)為0.959 3，優(yōu)于SVR和RVM，這是因?yàn)榻Ｋ斎氲腟EPA-VISSA特征波長(zhǎng)是在PLS建立子模型的基礎(chǔ)上所篩選的，因此，因變量(抗壓強(qiáng)度)的全部變異通過回歸關(guān)系被自變量(特征波長(zhǎng))解釋的比例，即決定系數(shù)是更高的。但是由于木材的各向異性及不均勻性等特點(diǎn)，特征波長(zhǎng)之間更多是處在非線性不可分的狀態(tài)，因此SVR和RVM非線性模型的預(yù)測(cè)均方根誤差分別為2.235 5和2.043 2，相對(duì)分析誤差分別為3.924 3和4.293 6，均優(yōu)于PLS建模，模型更穩(wěn)定和精確。此外，RVM相比于SVR引入了貝葉斯方法，提供了后驗(yàn)概率的輸出，能產(chǎn)生更稀疏的解，同時(shí)其調(diào)參也更為簡(jiǎn)便。從預(yù)測(cè)的結(jié)果數(shù)據(jù)來看，RVM模型的效果最優(yōu)，建模效果優(yōu)于PLS和SVR。

圖4 不同建模方法的模型效果對(duì)比Fig. 4 Comparison of model accuracies of different modeling methods

3 結(jié) 論

在900～1 700 nm近紅外光譜波段上，以樺木為例，提出SEPA-VISSA-RVM的木材抗壓強(qiáng)度預(yù)測(cè)建模方法。通過實(shí)例和數(shù)據(jù)證明該模型不僅在波長(zhǎng)優(yōu)選上具有優(yōu)勢(shì)，且相比于常用的PLS和SVR建模效果更好，在實(shí)際工程中具有良好應(yīng)用性。

1)利用SEPA作為VISSA的改進(jìn)策略提取出的111個(gè)特征波長(zhǎng)建立PLS模型，得到?jīng)Q定系數(shù)為0.959 3，預(yù)測(cè)均方根誤差為2.899 5，相對(duì)分析誤差為3.025 6，與VCPA、CARS和VISSA相比，在使用PLS線性建模時(shí)預(yù)測(cè)精準(zhǔn)度更高，提升了木材抗壓強(qiáng)度預(yù)測(cè)模型的準(zhǔn)確性和魯棒性。

2)將SPEA-VISSA與RVM建模方法相結(jié)合，通過PSO對(duì)RVM進(jìn)行參數(shù)優(yōu)化，選擇Laplacian核函數(shù)建模，得到?jīng)Q定系數(shù)為0.944 9，預(yù)測(cè)均方根誤差為2.043 2，相對(duì)分析誤差為4.293 6，相比于PLS和SVR，RVM作為非線性建模方法可更加精確地計(jì)算出樺木抗壓強(qiáng)度的預(yù)測(cè)值。