NIR光譜的kernel Isomap-PLS建模方法及其在血府逐瘀口服液提取過程分析的應用研究

申永祥， 楊輝華 ， 覃 鋒， 宋昊鯤， 羅國安

（1.桂林電子科技大學計算機與控制學院，廣西桂林 541004;2.清華大學分析中心，北京 100084;3.永州職業技術學院，湖南永州 425006）

1 引言

血府逐瘀湯［1，2］由桃仁、紅花、當歸、川芎、生地、赤芍、牛膝、柴胡等藥材組成，具有活血化瘀、行氣止痛的功效，用于瘀血內阻、頭痛或胸痛、內熱瞀悶、失眠多夢、心悸怔忡、急躁善怒等疾病。由于藥材成分的復雜性及批次間的差異等因素的影響，其產品的質量必定會存在差異，從而無法確保產品質量的穩定性和均一性。要解決這一問題，則必須實現對藥物體系中化學成分的直接控制。

過程分析技術（PAT）是實現工業生產過程在線、實時分析并實現反饋控制，從而保證生產過程可控、最優的重要技術［3］。近紅外光譜在線分析是一種最近廣泛應用的過程分析技術。它是利用物質對近紅外線有特征吸收的原理，并結合標準或認可方法測得的組成或性質數據建立分析模型，從而實現間接測量其組成或性質的一種分析方法。

降維是一種有效消除近紅外光譜數據噪聲并提取有用信息的方法［4，5］，PCA、PLS、ANN、SVR 已廣泛應用于近紅外光譜處理中。等距映射（Isomap）［6］是 Tenenbaum 等人于 2000年在 Science上提出的一種新的非線性降維方法。文獻［7］結合Isomap與PLS，提出了Isomap－PLS算法，建立的校正模型在準確性方面優于PLS。為了進一步提高校正模型的準確性，本文在Isomap－PLS算法的基礎上提出了一種新的算法 kernel Isomap－PLS，即先用 kernel Isomap［8］對近紅外光譜進行非線性降維，再用PLS建立校正模型（kernel Isomap－PLS算法）。實驗結果表明:kernel Isomap－PLS算法的建模效果優于 Isomap－PLS 與 PLS。

2 算法

2.1 Isomap算法

Isomap算法建立在多維尺度變換（MDS）的基礎上，力求保持數據點的內在幾何性質，即保持2點間的測地距離。算法描述如下［7］:

輸入:樣本x1，x2，…，xq，xi∈RP，q為樣本本真維數，k為鄰域參數。

第一步:構建k-鄰域圖G。計算每個樣本點的近鄰點，并在樣本集上定義一個賦權無向圖G。若xi與xj互為近鄰點，則邊的權值賦為

第二步:計算圖G中兩點間的最短路徑。用迪杰斯特拉（Dijkstra）算法或弗洛伊德（Floyd）算法求得最短路徑距離矩陣，記為DG。

第三步:計算q維嵌入。用MDS求低維嵌入流形。

① 代價函數為Dy=dy（i，j）。

2.2 kernel Isomap算法

kernel Isomap 算法描述如下［7］:

算法1:kernel Isomap算法

輸入:樣本x1，x2，…，xq，xi∈RP，q為樣本本真維數，k為鄰域參數。

第一步:構建k-鄰域圖G（同Isomap算法第一步）。

第二步:計算最短路徑（同 Isomap算法第二步）。

其中H=I－（1/n）eeT，e=［1，…，1］T∈Rn

第四步:計算矩陣

的最大特征值c*。構建Mercer核矩陣:

第五步:計算ˉK的最大d個特征值和對應的特征向量，其中特征值矩陣Λ∈Rd×d，特征向量V∈Rn×d。

第六步:相應的在d維空間的n個點以M=VΛ1/2得出。

因為矩陣ˉK是Mercer核矩陣，所以它的（i，j）元素可以被下式取代。

其中φ（g）是到特征空間的非線性映射。利用核技巧，可以計算出測試數據xk在低維空間的投影。

其中，［g］i是向量的第i個元素，vi是ˉK的第i個特征向量。投影一個新的變量算法總結如下:

算法2:投影一個測試樣本

輸入:Mercer核矩陣ˉK和它的特征向量V，加法常數c，測試樣本xk。

輸出:xk在d維空間的投影mk。

第一步:計算xk與原始輸入空間的任一xj樣本點的最短路徑dG（k，j）。

第二步:計算測試數據xk的核:

其中，ˉdG（k，j）=dG（k，j）+c，j=1，…，n。

第三步:用（5）式計算出相應的投影mk。

2.3 kernel Isomap－PLS 算法

第二步:用上面kernel Isomap算法2求出測試樣本xk的低維嵌入rk。

第三步:進行PLS分解。M=TS+E，Y=UQ+F

其中，T和U分別為M和Y的得分矩陣;S和Q分別為M和Y的載荷矩陣;E和F分別為M和Y的PLS擬和殘差矩陣。

第四步:將T和U作線性回歸:U=TB，B=（T′T）－1T′Y。

第五步:計算預測值。由S和測試集的低維嵌入求出測試集的得分Tscore，然后根據Tscore和回歸系數B，可以計算出測試集的預測值:Y=TscoreBQ。

3 實驗

3.1 數據集描述

從吉林敖東延邊藥業服份有限公司的血府逐瘀口服液的提取生產過程中，利用北京英賢儀器有限公司的INCE9500MT近紅外光譜儀在線采集其近紅外光譜，采用透射方式測量，光程2 mm，光譜范圍1 000～2 200 cm，波長間隔4.8 cm。測光譜的同時抓取樣本共163個，所得提取液光譜如圖1所示。

采用HPLC檢測提取液中羥基紅花素－A和甘草酸銨的含量。儀器使用SH MADZU－LC－2－10A 高效液相色譜儀。在樣品中，測得羥基紅花素－A的最大含量為0.195 mg/mL，最小含量為0.015 mg/mL，甘草酸銨的最大含量為0.192 mg/mL，最小含量為0.066 mg/mL。

圖1 INCE9500MT近紅外光譜儀在線采集血府逐瘀口服液提取液163個樣本光譜圖

3.2數據處理

對所得光譜數據，分別用kernel Isomap－PLS與PLS建立回歸模型，模型采用留一法交叉驗證均方根誤差（RMSECV）或測試集均方根誤差與決定系數（R2）來評價:若RMSECV或RMSEP越小，R2越大，則模型越優。kernel Isomap算法采用Choi開發的程序包［8］，在MatlabTM 7.1上自編程實現了參數優化與其它相關程序。

4 結果與討論

4.1 ernel Isomap算法與Isomap一樣有兩個可調整的參數:鄰域參數k和樣本本真維數d、k和d均取整數，對k和d采用網絡搜索法進行優化。在數據集中，對樣本集中所有樣本均選取d從3到80，k從3到70，用kernel Isomap方法進行降維。降維后的數據再用PLS建立回歸模型，模型評價采用留一法交叉驗證均方根誤差（RMSECV）與決定系數（R2），RMSECV越小，R2越大，模型越優。當d=60，k=52時，羥基紅花素－A具有最小的RMSECV值與最大的R2值，RMSECV=0.038 9，R2=0.935 2;當d=32，k=43時，甘草酸銨的RMSECV值最小，為0.034 4。R2值最大，為0.957 1。

4.2 kernel Isomap－PLS算法和PLS算法在校正集上建模效果與測試集上預測結果的比較

在樣本集中，隨機抽取15個樣本作為測試集，其余樣本作為校正集。分別用kernel Isomap－PLS、Isomap－PLS與PLS在校正集上采用留一法交叉驗證建立回歸模型。在校正集中采用留一法交叉驗證均方根誤差（RMSECV）與決定系數（R2）來評價模型，RMSECV越小，R2越大，模型越優;在測試集中采用測試集均方根誤差與決定系數（R2）來評價模型，RMSEP越小，R2越大，預測結果越準確。模型結果如表1所示。從表中可以看出，kernel Isomap－PLS建立的校正模型及預測結果均優于PLS與Isomap－PLS。

表1 kernel Isomap-PLS與Isomap-PLS、PLS在校正集上建模效果和測試集上預測結果的比較

5 結論

kernel Isomap能夠對NIR光譜進行流形降維，剔除噪聲，提取有用信息，與PLS結合后可實現NIR光譜的非線性建模。實驗結果表明kernel Isomap－PLS算法比單獨的PLS算法與Isomap－PLS算法建模與預測效果更好。該方法可實現血府逐瘀口服液提取過程的在線檢測和質量控制。

［1］丁志山，高承賢，程東慶，等.血府逐瘀湯對牛內皮細胞增殖和遷移的影響［J］.中成藥，2003，25（5）:423－424.

［2］王 雁，瑞 龍.不同提取方法對血府逐瘀湯中芍藥苷含量的影響［J］.中成藥，1998，20（3）:11－12.

［3］陸婉珍，褚小立.近紅外光譜（NIR）和過程分析技術（PAT）［J］.現代科學儀器，2007，4:13－17.

［4］楊輝華，覃 鋒，王 勇，等.NIR光譜的LLE－PLS非線性建模方法及應用［J］.光譜學與光譜分析，2007，10（10）:1955－1958.

［5］覃 鋒，楊輝華，呂琳昂，等.NIR光譜結合LLE－PLS建模用于安神補腦液提取過程分析的研究［J］.中成藥，2008，30（10）:1465－1468.

［6］Tenenbaum J B，Silva V D，Langford J C.A global geometric framework for nonlinear dimensionality reduction［J］.Science，2000，290（22）:2319－2323.

［7］覃 鋒.Isomap用于中藥生產過程近紅外光譜在線檢測研究［D］，桂林電子科技大學，2008.

［8］Choi H，Choi S.Kernel Isomap ［J］.Electron Lett，2004，40（25）:1612－1613.