基于L1-L2聯(lián)合范數(shù)約束的中藥近紅外光譜波長(zhǎng)選擇

任 真 李四海

(甘肅中醫(yī)藥大學(xué)信息工程學(xué)院 甘肅 蘭州 730000)

0 引 言

中藥具有物質(zhì)成分多、作用機(jī)理復(fù)雜等特點(diǎn)，保證中藥產(chǎn)品質(zhì)量的可靠穩(wěn)定是實(shí)現(xiàn)中藥現(xiàn)代化必須要解決的關(guān)鍵問(wèn)題。近紅外光譜分析技術(shù)是20世紀(jì)90年代以后迅速發(fā)展起來(lái)的一種新型在線檢測(cè)技術(shù)，具有簡(jiǎn)便、低成本、不破壞樣品等優(yōu)點(diǎn)，已在農(nóng)業(yè)、食品、石油化工、藥物質(zhì)量控制等領(lǐng)域得到廣泛應(yīng)用。在中藥質(zhì)量控制方面，近紅外光譜能夠快速、準(zhǔn)確地鑒別中藥材的真?zhèn)?、種類(lèi)和產(chǎn)地，并且能夠快速測(cè)定中藥材中有效成分的含量以及中藥輔料的品質(zhì)[1-2]。

近紅外光譜信號(hào)具有維度高、變量多重共線性嚴(yán)重等特點(diǎn)[3]。文獻(xiàn)[4]用隨機(jī)蛙跳算法(Random-frog)對(duì)生物柴油近紅外光譜進(jìn)行波長(zhǎng)選擇，分別建立了生物柴油含水量的偏最小二乘和支持向量機(jī)預(yù)測(cè)模型，取得了較好的效果。文獻(xiàn)[5]建立了金銀花醇沉過(guò)程中綠原酸含量的偏最小二乘回歸模型，通過(guò)競(jìng)爭(zhēng)自適應(yīng)抽樣CARS變量選擇方法，提高了模型的預(yù)測(cè)精度。以上研究表明，基于波長(zhǎng)變量選擇的近紅外光譜建模方法能夠有效提高模型的預(yù)測(cè)能力。

現(xiàn)有的波長(zhǎng)變量選擇方法大多是封裝式的：通過(guò)遍歷搜索，找到變量空間的最優(yōu)特征子集，使模型的均方根誤差最小。由于特征子集的搜索采用前向、后向或蒙特卡洛隨機(jī)搜索策略，對(duì)每一個(gè)特征子集的評(píng)價(jià)都需要重新訓(xùn)練模型，計(jì)算開(kāi)銷(xiāo)較大。嵌入式特征選擇是機(jī)器學(xué)習(xí)中一類(lèi)重要的特征選擇方法，其通過(guò)最小化目標(biāo)函數(shù)，使得特征選擇和模型訓(xùn)練融為一體，同步完成，由于L1范數(shù)的稀疏特性，嵌入式特征選擇能夠自動(dòng)實(shí)現(xiàn)變量選擇，降低模型的過(guò)擬合風(fēng)險(xiǎn)，提高模型的可解釋性。

本文通過(guò)在偏最小二乘回歸模型中同時(shí)引入L1和L2正則項(xiàng)，建立一種嵌入式的光譜特征變量選擇方法。將光譜變量選擇和預(yù)測(cè)模型的建立融合在一起，解決光譜數(shù)據(jù)存在的多重共線性問(wèn)題，提高了偏最小二乘回歸模型的可解釋性，實(shí)現(xiàn)了對(duì)當(dāng)歸中藁本內(nèi)酯含量的快速、準(zhǔn)確檢測(cè)。

1 正則偏最小二乘回歸算法設(shè)計(jì)

1.1 L1正則與稀疏性

正則化具有產(chǎn)生稀疏模型的能力[6-9]，Tibshirani[10]提出的線性回歸Lasso模型通過(guò)將嶺回歸中的L2正則項(xiàng)替換為L(zhǎng)1正則項(xiàng)，使模型具有變量選擇和數(shù)據(jù)降維能力。

假設(shè)自變量矩陣X∈Rn×p，因變量Y∈Rn×1，線性回歸的Lasso模型為：

(1)

式中：‖·‖2表示L2范數(shù)，β為回歸系數(shù)向量，‖β‖1表示所有回歸系數(shù)的絕對(duì)值之和，為L(zhǎng)1范數(shù)罰，λ為懲罰系數(shù)。Lasso模型的解可通過(guò)坐標(biāo)下降算法求得[11]：

對(duì)于i=1,2,…,p

1.2 偏最小二乘回歸

偏最小二乘回歸(PLSR)是主成分分析(PCA)和典型相關(guān)分析(CCA)的有效結(jié)合，其對(duì)CCA方法進(jìn)行了進(jìn)一步拓展。PLSR能有效解決高維變量之間的多重共線性問(wèn)題，在近紅外光譜的定量分析中得到了廣泛應(yīng)用[14]。

假設(shè)X和Y分別為光譜數(shù)據(jù)矩陣和待測(cè)含量矩陣，X∈Rn×p，Y∈Rn×q，X和Y的每一列均為零均值且標(biāo)準(zhǔn)差為1。

PLSR首先從X和Y中提取第一對(duì)主成分u1和t1，滿足：Var(u1)→max，Var(t1)→max，Corr(u1,t1)→max。即u1和t1最大化且二者的相關(guān)性最大化。然后計(jì)算X和Y殘差并根據(jù)殘差繼續(xù)提取下一個(gè)主成分。由于各投影方向之間相互正交，抽取的特征位于不同的投影方向，因此偏最小二乘能夠有效去除變量之間的多元共線性。PLSR的最優(yōu)主成分個(gè)數(shù)一般通過(guò)交叉驗(yàn)證方法確定。

1.3 正則偏最小二乘回歸

近紅外光譜數(shù)據(jù)具有高維度、小樣本的特點(diǎn)。PLSR的投影向量是所有原始波長(zhǎng)變量的線性組合，一方面，將部分噪聲變量也納入投影向量參與模型擬合會(huì)使有效變量的回歸系數(shù)產(chǎn)生衰減，導(dǎo)致預(yù)測(cè)精度下降；另一方面，當(dāng)p>>n，即p遠(yuǎn)大于n時(shí)，投影向量包含所有的原始波長(zhǎng)變量導(dǎo)致模型可解釋性不強(qiáng)。針對(duì)上述問(wèn)題，在偏最小二乘回歸模型中引入正則項(xiàng)，建立正則偏最小二乘回歸RPLS算法，RPLS可形式化為如下的最優(yōu)化問(wèn)題：

(2)

s.t.αTα=1

式中：η、λ1和λ2均為常量，用于初始化算法，c為主成分載荷系數(shù)向量α的副本且與α取值相近。 式(2)可通過(guò)交替迭代算法求解[15]：

1) 固定c，求解α。對(duì)固定的c，式(2)變?yōu)椋?/p>

(3)

s.t.αTα=1

令Z=XTY,η′=(1-η)/(1-2η)，則式(3)可重寫(xiě)為：

(4)

s.t.αTα=1

式中的α可通過(guò)拉格朗日乘子法求解。

2) 固定α，求解c。對(duì)固定的α，式(2)變?yōu)椋?/p>

(5)

問(wèn)題轉(zhuǎn)化為因變量為ZTα的彈性網(wǎng)問(wèn)題[16]，c可以通過(guò)LARS算法求解[17]。

式中：Z1=XTY/‖XTY‖為第一投影方向單位向量。

Step2 當(dāng)k≤K時(shí)

1.4 算法參數(shù)選擇

對(duì)于單因變量回歸問(wèn)題，一般取η=1/2,λ2→∞。因此，算法的關(guān)鍵參數(shù)有兩個(gè)：L1正則項(xiàng)系數(shù)λ1和最優(yōu)主成分個(gè)數(shù)k。

λ1用于控制選擇的波長(zhǎng)變量個(gè)數(shù)，其值越大，選擇的波長(zhǎng)變量數(shù)越少，模型可解釋性越好，但參與建模的變量數(shù)過(guò)少會(huì)導(dǎo)致模型預(yù)測(cè)能力下降。因此，λ1的選取要權(quán)衡稀疏度和模型預(yù)測(cè)能力，通過(guò)實(shí)驗(yàn)選擇最優(yōu)的λ1值。

最優(yōu)主成分的個(gè)數(shù)k通過(guò)留一交叉驗(yàn)證法計(jì)算得到的Qk2值來(lái)確定[18]。

Qk2= 1-PRESSk/SS(k-1)

(6)

式中：PRESSk為使用前k個(gè)主成分對(duì)預(yù)留樣本預(yù)測(cè)誤差的平方和，SS(k-1)為使用前k-1個(gè)主成分對(duì)所有樣本擬合誤差的平方和。當(dāng)Qk2≥0.097 5時(shí)，認(rèn)為第k個(gè)主成分作用顯著。

2 模型的建立

2.1 實(shí)驗(yàn)數(shù)據(jù)

采集甘肅不同產(chǎn)地的當(dāng)歸樣本76個(gè)，使用美國(guó)Thermo公司的Nicolet-6700型近紅外光譜儀掃描得到所有樣本的近紅外光譜，光譜范圍為4 000 cm-1～10 000 cm-1，全譜共包括1 557個(gè)波數(shù)變量。76個(gè)當(dāng)歸樣本的近紅外光譜如圖1所示。

圖1 當(dāng)歸樣本的近紅外光譜

為消除基線漂移并提高光譜信號(hào)的信噪比，對(duì)原始光譜進(jìn)行一階導(dǎo)數(shù)結(jié)合正交信號(hào)校正預(yù)處理。當(dāng)歸中藁本內(nèi)酯的含量采用高效液相色譜法(HPLC)測(cè)定[19]。

2.2 波長(zhǎng)變量選擇

將樣本劃分為訓(xùn)練集和測(cè)試集，訓(xùn)練集樣本56個(gè)，測(cè)試集樣本20個(gè)。使用訓(xùn)練集建立正則偏最小二乘回歸模型對(duì)當(dāng)歸近紅外光譜進(jìn)行波長(zhǎng)變量選擇，因變量為當(dāng)歸中的藁本內(nèi)酯含量。訓(xùn)練集和測(cè)試集中藁本內(nèi)酯含量的分布情況見(jiàn)表1。

表1 訓(xùn)練集和測(cè)試集中藁本內(nèi)酯含量分布

實(shí)驗(yàn)中用eta參數(shù)控制變量選擇個(gè)數(shù)和模型稀疏度，對(duì)eta和k的最優(yōu)組合采用網(wǎng)格尋優(yōu)法確定。設(shè)定主成分k的搜索范圍為[2,8]，步長(zhǎng)為1，eta的搜索范圍為[0.7,0.9]，步長(zhǎng)為0.1。最優(yōu)波長(zhǎng)變量子集通過(guò)權(quán)衡訓(xùn)練集上的預(yù)測(cè)均方根誤差和波長(zhǎng)變量個(gè)數(shù)來(lái)最終確定。

3 實(shí)驗(yàn)結(jié)果與分析

3.1 波長(zhǎng)變量選擇結(jié)果

根據(jù)主成分?jǐn)?shù)k和參數(shù)eta的搜索范圍，共得到21個(gè)波長(zhǎng)變量子集，分別建立以上變量子集在56個(gè)訓(xùn)練樣本上的偏最小二乘回歸模型。不同變量子集在測(cè)試集上的預(yù)測(cè)結(jié)果對(duì)比如圖2所示。

圖2 k和eta對(duì)選擇波長(zhǎng)變量個(gè)數(shù)的影響

從圖2可知，主成分個(gè)數(shù)k和eta對(duì)選擇的波長(zhǎng)變量個(gè)數(shù)都會(huì)產(chǎn)生影響，變量個(gè)數(shù)隨著k的增加而增加，其中eta參數(shù)對(duì)變量個(gè)數(shù)的影響較大，eta值越大，選擇的波長(zhǎng)變量個(gè)數(shù)越少。當(dāng)eta為0.7時(shí)，選擇變量個(gè)數(shù)最少為100個(gè)，最多為119個(gè)，模型的擬合效果很好，但參與建模的變量過(guò)多，模型的可解釋性不強(qiáng)。當(dāng)eta為0.8時(shí)，選擇變量個(gè)數(shù)最少為42個(gè)，最多為67個(gè)，此時(shí)選擇的變量仍然過(guò)多。當(dāng)eta為0.9時(shí)，選擇變量個(gè)數(shù)最少為15個(gè)，最多為22個(gè)。綜合考慮選擇變量個(gè)數(shù)和模型的預(yù)測(cè)能力，本文選擇eta=0.9,k=3為最優(yōu)參數(shù)，此時(shí)共選擇16個(gè)波長(zhǎng)變量。分別為5 164.4 cm-1、5 218.4 cm-1、5 222.3 cm-1、5 226.1 cm-1、5 230.0 cm-1、5 233.9 cm-1、5 237.7 cm-1、5 241.6 cm-1、5 245.4 cm-1、5 249.3 cm-1、5 253.1 cm-1、5 257.0 cm-1、5 260.9 cm-1、5 704.4 cm-1、5 708.3 cm-1、 5 712.1 cm-1。根據(jù)近紅外光譜的特征峰理論，所選擇的波數(shù)大多位于C=O基團(tuán)的伸縮振動(dòng)吸收峰附近，與當(dāng)歸中的內(nèi)酯類(lèi)化合物有關(guān)，這說(shuō)明所選擇的波長(zhǎng)變量具有較好的化學(xué)意義。

3.2 不同預(yù)測(cè)方法結(jié)果對(duì)比

使用RPLS選擇的16個(gè)波長(zhǎng)變量，在56個(gè)訓(xùn)練樣本上建立偏最小二乘回歸模型，取前3個(gè)主成分。對(duì)20個(gè)測(cè)試樣本中的藁本內(nèi)酯含量的預(yù)測(cè)結(jié)果見(jiàn)圖3。

圖3 藁本內(nèi)酯預(yù)測(cè)值和實(shí)際值對(duì)比

從3可以看出，預(yù)測(cè)值和真實(shí)值之間非常接近，預(yù)測(cè)均方根誤差RMSEP(root mean square error of prediction)為0.181 6，決定系數(shù)為0.995 7。RMSEP的計(jì)算公式如下：

(7)

為進(jìn)一步說(shuō)明本文波長(zhǎng)變量選擇方法的有效性，將RPLS與目前近紅外光譜分析中常用的CARS、Random-frog等波長(zhǎng)變量選擇方法進(jìn)行了對(duì)比分析。CARS算法迭代次數(shù)為1 200次，5折交叉驗(yàn)證，從1 200個(gè)子模型中選擇交叉驗(yàn)證均方根誤差最小的模型。Random-frog算法變量迭代次數(shù)為1 000次，實(shí)驗(yàn)中迭代達(dá)到200次時(shí)選擇的變量結(jié)果基本穩(wěn)定。根據(jù)三種方法所選擇的變量，分別建立PLSR模型，表2對(duì)比了不同的波長(zhǎng)變量選擇方法的預(yù)測(cè)性能。

表2 不同變量選擇方法的預(yù)測(cè)結(jié)果對(duì)比

從表2可以看出，與全譜建模相比，根據(jù)CARS和隨機(jī)蛙跳選出的重要變量建模，能夠獲得與全譜建模相當(dāng)?shù)念A(yù)測(cè)性能。但隨機(jī)蛙跳需要預(yù)先設(shè)置種子變量，然后采用蒙特卡洛抽樣技術(shù)，將其他波長(zhǎng)變量依次添加到種子變量中形成不同的變量子集，根據(jù)均方根誤差確定最優(yōu)的波長(zhǎng)變量子集，其迭代次數(shù)和變量選擇結(jié)果均依賴(lài)于種子變量初始值，算法的計(jì)算開(kāi)銷(xiāo)較大。與CARS和隨機(jī)蛙跳相比，RPLS波長(zhǎng)變量選擇方法通過(guò)設(shè)置eta和k的值來(lái)控制變量選擇個(gè)數(shù)，并不是以變量子集的預(yù)測(cè)均方根誤差來(lái)衡量變量子集的優(yōu)劣，而是在變量稀疏性和預(yù)測(cè)能力之間取得折中，方法更為穩(wěn)健。RPLS方法選擇的波長(zhǎng)變量數(shù)最少，預(yù)測(cè)性能優(yōu)于CARS和隨機(jī)蛙跳，且所選擇的變量具有較好的化學(xué)意義，模型可解釋性較好。

4 結(jié) 語(yǔ)

通過(guò)在偏最小二乘回歸模型中引入L1和L2范數(shù)罰正則項(xiàng)，建立了正則偏最小二乘波長(zhǎng)選擇方法。該方法能夠?qū)⒔t外光譜中噪聲變量在主成分上的載荷系數(shù)置為0，保留有效變量，達(dá)到選擇重要變量的目的。與CARS和隨機(jī)蛙跳變量選擇方法相比，RPLS變量選擇方法在選擇的波長(zhǎng)數(shù)、模型的預(yù)測(cè)精度及可解釋性等方面均具有一定優(yōu)勢(shì)。

本文提出的正則偏最小二乘波長(zhǎng)選擇方法對(duì)噪聲變量和有效變量施加相同的懲罰，對(duì)主成分載荷系數(shù)的估計(jì)是有偏估計(jì)。如何減弱甚至消除對(duì)有效變量的懲罰，得到載荷系數(shù)的近似無(wú)偏估計(jì)，提高近紅外光譜波長(zhǎng)變量選擇的針對(duì)性將是下一步的研究方向。