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.