近红外光谱信息筛选在玛咖产地鉴别中的应用

王元忠，赵艳丽，张 霁，金 航

云南省农业科学院药用植物研究所，云南 昆明 650200

近红外光谱信息筛选在玛咖产地鉴别中的应用

王元忠，赵艳丽，张 霁，金 航*

云南省农业科学院药用植物研究所，云南 昆明 650200

食药植物玛咖富含多种营养成分，极具药用价值。采用近红外漫反射光谱，对采自秘鲁及云南共139份玛咖样品进行产地鉴别。采用多元信号校正结合二阶导数和Norris平滑预处理光谱，利用光谱标准偏差初步选择光谱波段(7 500～4 061 cm-1)，结合主成分-马氏距离(principal component analysis-mahalanobis distance，PCA-MD)筛选出适宜的主成分数为5。基于所筛选的光谱波段及主成分数，利用“模群迭代奇异样本诊断”方法剔除2个异常样品后，分别采用竞争自适应重加权法(competitive adaptive reweighted sampling，CARS)、蒙特卡洛-无信息变量消除法(monte carlo-uninformative variable elimination，MC-UVE)、遗传算法(genetic algorithm，GA)和子窗口重排(subwindow permutation analysis，SPA)四种方法筛选光谱变量信息，利用模型集群分析(model population analysis，MPA)思想对所筛选的光谱变量信息进行评价。结果显示，RMSECV(SPA)>RMSECV(CARS)>RMSECV(MC-UVE)>RMSECV(GA)，分别为2.14，2.05，2.02，1.98，光谱变量数分别为250，240，250和70。采用偏最小二乘判别分析法(partial least squares discriminant analysis，PLS-DA)对四种方法筛选的光谱变量建立判别模型，随机选择97份样品作为建模集，其余40份样品作为验证集。通过R2，RMSEC和RMSEP分析可知，R2： GA>MC-UVE>CARS>SPA，RMSEC和RMSEP： GA

玛咖； 近红外光谱； 鉴别； 光谱信息筛选； 模型集群分析

引 言

玛咖(LepidiummeyeniiWalp.)为十字花科(Cruciferae)独荇菜属(Lepidium)植物，常用名Maca(音译“玛咖”)，又名Maka，Maca-maca，Peruvian ginseng，Maino，Ayakwillku，Ayakchichira等，原产于海拔3 500～4 500 m的南美洲安第斯山区，现主要分布于秘鲁中部的Puna生态区和秘鲁东南部城市Puno[1-2]。玛咖块根可药食两用，极具营养价值和经济价值，富含蛋白质、人体必需氨基酸、矿物质和多种维生素，在安第斯山区已有悠久的使用历史。研究发现玛咖中除含有多种营养成分外[3]，还含有玛咖酰胺、玛咖烯、芥子油苷及其衍生物、甾醇等多种活性物质[4-6]。现代药理学研究表明玛咖具有抗氧化[7]、增强生育力[8]、改善性功能[9]、抗疲劳[10]、调节内分泌[11]、增强免疫力[12]、抗压力[13]、缓解更年期综合症[14]等多种功效。近年来，玛咖因其独特功效而倍受关注。我国云南、新疆、吉林、西藏和湖南等省(区)均有引种栽培，其中以云南省种植规模最大、品种最多。据统计2012年，云南全省玛咖种植面积达22 000 hm2，已成为云南省高寒山区的重要经济作物[15]。然而，不同产地种植的玛咖其营养价值和药用价值存在一定差异，且随着市场需求逐年增加，市售玛咖品质差异较大，粉末样品掺杂现象突出。因此，不同产地玛咖样品的鉴别，为后续指导玛咖GAP种植、品质及安全评价具有重要的作用。

近红外光谱(NIRS)对复杂体系中C—H，O—H，N—H等含氢基团振动的倍频和合频吸收极为敏感，全光谱具有同时反映物料的物化信息等特点。但近红外全光谱变量中，很多是无效信息、甚至是干扰变量，引入此类变量建立判别模型增加了模型的复杂性，同时降低了模型的判别精度，进而对样本的鉴别和变量的计算产生重要影响[16]。为建立较为准确的判别模型，需对光谱信息进行选择及优化，使所筛选的光谱信息具有全面性和客观性。收集了139份不同产区种植的玛咖样品，拟采用竞争自适应重加权法(CARS)[17]、蒙特卡洛-无信息变量消除法(MC-UVE)[18-19]、遗传算法(GA)[20]和子窗口重排(SPA)[21]筛选玛咖样品的近红外光谱变量后，利用模型集群分析(MPA)[22]思想对所筛选的光谱变量信息进行评价，并确定筛选的光谱波数，进一步采用偏最小二乘判别分析法(PLS-DA)建立各方法的判别模型，用于预测验证样品。该方法旨在为探析构建近红外光谱判别模型的新方法提供参考，通过集群分析结合PLS-DA建立不同产地玛咖样品的判别模型，以期为近红外光谱鉴别和评价中药材产地、品种等提供依据。

1 实验部分

1.1 材料

2014年采集15个不同产地种植的药用植物玛咖139份，经鉴定，样品来源见表1。

Table 1 Source and number of materials

1.2 仪器

AntarisⅡ近红外光谱仪(美国赛默飞世尔科技，配备InGaAs检测器的漫反射模块，Result 2.1软件采集光谱图、TQ 8.6软件定性分析)； DFT-100型中药粉碎机(浙江温岭市林大机械有限公司)，80目不锈钢筛盘(北京，中西泰安)； SIMCA-P+11.0软件(瑞典，UMETRICS)； MATLAB R2010a分析软件，代码来源于http://code.google.com/p/carspls/和http://www.mathworks.cn/。

1.3 样品制备

样品采集后用自来水洗净根茎部位，再用蒸馏水冲洗3次，将其切成薄片，置于干净的白纸上在常温下阴干。用中药粉碎机粉碎，过80目筛，样品不少于20.0 g，存放于自封袋中，备用。

1.4 近红外光谱采集

将样品置于平衡箱中平衡2 h，使样品水分含量在10%～12%。准确称取平衡后的样品20.0 g，充分混匀，置于近红外采样杯中，压紧。在仪器预热2 h后，使用Result 2.1软件漫反射模块采集近红外光谱图。测试条件： 扫描次数： 64次，分辨率： 4 cm-1，扫描范围： 10 000～4 000 cm-1。样品平行测定3次，取平均光谱。随机选择99份样品作为训练集，其余40份样品作为验证集。采用TQ 8.6软件对光谱进行优化处理，消除基线漂移和噪音，结果见图1和图2。

Fig.1 Original NIRS spectra of Maca

Fig.2 Second derivative spectra of NIRS of Maca

2 数据处理

2.1 近红外光谱信息初步筛选[23]

利用TQ 8.6软件，使用MSC+SD+ND(13∶3)方法初步优化所采集的光谱，对不同产地样品类型分别赋值： 香格里拉1#“1”，香格里拉2#“2”，香格里拉3#“3”至轿子雪山13#“13”，大山包14#“14”，秘鲁15#“15”。通过光谱标准偏差选择光谱波段，利用主成分-马氏距离(PCA-MD)初步建立分类模型，并对主成分数进行选择，结果见图3和表2。

2.2 近红外光谱集群分析

2.2.1 奇异样品诊断

由于光谱数据(X)和量测指标(Y)的采集和测定中，仪器的稳定性及异常波动会使部分X或Y偏离整体分布的数据点，导致所建立的判别模型存在奇异点。为提高判别模型的预测精度，建立稳健的判别模型，采用SIMCA-P+11.0软件将2.1中所筛选的光谱波段(7 500～4 061 cm-1)转置为数据格式，并利用“模群迭代奇异样本诊断”方法通过多次循环校正，将所选择的光谱数据(X=7 500～4 061 cm-1)与量测指标(Y=产区分类赋值1～15)建立判别模型，并计算模型的偏度和峰度，结果见图4。

Table 2 Contribution of principal component

Fig.3 Spectral range selected by standard deviation

Fig. 4 Singular sample diagnosis

2.2.2 集群分析优化光谱变量

剔除2.2.1中异常样本1—2和12—6后，选择主成分数为5，对光谱数据与产地分类指标进行判别分析。将建模集样品和验证集样品按0.8的比例(即Ratio=0.8)，分别利用CARS(运行次数为N=100次，主成分数A=5，交叉验证数K=10，统计每个光谱波数选中频率)、MC-UVE、GA(MC-UVE运行次数为N=100次，主成分数A=5； GA运行次数为N=100次，Autoscaling=2，Deletion groups=5，Chromosomes=30，Variables as a maximum=30，Probability of mutation=0.01, Probability of cross-over=0.5，统计每个光谱波数选中频率)和SPA(主成分数p.component=5，显著性水平p.sig=0.01，载荷p.Q=0.3)优化光谱波数，计算标表征变量重要性的指标，CARS，MC-UVE，GA和SPA分别为频率(Frequency，F)、变量稳定指数(reliability index，RI)[17]、频率和条件协同得分(conditional synergetic score，COSS)[22]，按降序排列，以蒙特卡洛采样技术，采用向前推移选择变量重要性指标引入变量(最大变量数为250)重复建模，使用蒙特卡洛交互验证评价预测性能[22]，以RMSECV值最小或变化趋是较小时确定最佳变量数，结果见图5。

Fig.5 RMSECV change trend chart

采用上述四种方法对所筛选的光谱波数进行分析，结果见图6—图10。由图6—图10可知，CARS法选中的频率大于1%的光谱波数分布于少数几个波段，其中主要分布于4 072～4 539，4 917～4 933，5 191～5 993，7 046～7 116 cm-1等波段，选中的光谱变量数为240个； GA法选中的频率大于4%的光谱波数分布均匀，主要集中于4 315～7 254 cm-1，选中的光谱变量数为70个； MC-UVE法选中的变量稳定指数(RI)大于3的光谱波数主要分布于4 396～7 170 cm-1，选中的光谱变量数为250个； SPA法选中的条件协同得分(COSS)大于1的光谱波数主要分布于4 072～6 981 cm-1，选中的光谱变量数为240个。

Fig.6 Spectrum wave number selected by CARS

Fig.7 Spectrum wave number selected by GA

Fig.8 Spectrum wave number selected by MC-UVE

Fig.9 Spectrum wave number selected by SPA

2.3 PLS-DA判别模型的建立

采用SIMCA-P+11.0软件对2.2.2中四种方法优化的光谱波数及其分类赋值建立判别模型，剔除奇异样品后，训练集为97份，验证集为40份，分类赋值见2.2.1。以GA法为例，建立判别模型，用于预测验证集样品。采用决定系数(R2)、校准均方根误差(RMSEC)和预测均方根误差(RMSEP)三个参数评估预测模型的准确性[23]，结果见图10和表3。

Fig.10 3D plot of PLS-DA by GA analysis

Table 3 Results of PLS-DA models for discrimination of different geographical Maca samples by GA, MC-UVE, CARS and SPA analysis

续表3

11-4111110.5720.30311UI11.5730.405111110.8470.108111110.8270.12211-9111111.0000.000111111.4070.288111110.4730.373111110.6530.24612-2121211.8600.099121212.4210.297121211.9260.052121212.1020.07212-5121212.4780.338121212.4320.30612UI12.5690.40212DE12.7960.56313-2131313.1500.106131313.2440.172131312.5470.320131313.0680.04813-6131312.6860.222131313.2230.15813UI12.3070.490131312.6690.23413-9131312.6330.260131313.2280.161131312.6680.235131312.8670.09414-2141413.5980.284141414.2670.18914DE13.2380.539141413.4540.38614-4141414.0000.000141413.7910.148141413.6320.26014UI13.3640.45014-8141414.1290.091141413.5330.330141413.5250.336141413.7220.19715-1151514.8850.082151514.9680.023151514.6440.252151514.8910.07715-5151514.5720.302151515.2650.187151515.0630.045151515.1280.09115-8151514.8020.140151515.1600.113151514.5430.323151514.5860.293正确率/%95.0092.5090.0085.00R20.99620.99380.98900.9864RMSEC0.180.230.310.34RMSEP0.280.350.470.52

Note： AC: actual class; CC: calculated class; Ypre: predicted value; Ydev: deviation; UI: unidentified; DE: discriminant error

3 结果与讨论

3.1 近红外光谱初步筛选分析

由图3和表2可知，15个不同产区139份样品建立的PCA-MD分类模型，光谱波段选择7 500～4 061 cm-1，主成分数为5时，提取的光谱信息较为充分，提取的总光谱贡献率为98.192 7%，选择的7 500～4 061 cm-1光谱波段提取贡献率为98.9995%，且随主成分数的增加，光谱信息数据量无明显变化。因此，选择主成分数为5较为适宜。

3.2 近红外光谱集群分析

由图4分析可知，共剔除奇异点(outlier)数2个，即样品1—2和12—6，计算出模型的偏度和峰度分别为-0.31和-0.52。结果表明，样品1—2和12—6可能在所筛选的光谱波段与其余样品存在较大差异，致使样品出现奇异。通过模型偏度和峰度分析可知，在所选择的光谱波段范围内，分析样品较为集中，说明该光谱波段对模型的建立具有重要影响。

由图5可知，采用四种方法对光谱波数进行筛选，RMSECV(SPA)>RMSECV(CARS)>RMSECV(MC-UVE)>RMSECV(GA)，分别为2.14，2.05，2.02，1.98； 且变量个数分别为250，240，250和70。分析结果显示： 四种集群分析方法中，通过综合评价选择的光谱波数和RMSECV值，其中GA法最优，CARS法和MC-UVE法次之，SPA法稍差。

由图6—图10分析可知，不同产区采集的玛咖样品在近红外光谱分析中，其光谱差异波段主要集中在4 300～7 200 cm-1之间，且光谱波数具有不连续性，四种方法筛选玛咖样品近红外光谱波数均具有一定的代表性，其中GA法筛选的光谱波数较少，且分布范围较广，具有较好的代表性。

3.3 PLS-DA判别模型分析

由3D图10可知，秘鲁15#、香格里拉2#、昭通5#、大山包14#、轿子雪山13#等能清晰分开，其余产地样品虽集聚在一起，但仍能清晰区分。由表3可知，通过上述四种方法筛选光谱波数，采用PLS-DA对筛选的光谱波数与其分类赋值进行分析，并对验证样品进行预测。结果显示： 四种方法筛选的光谱波数建立的判别模型预测验证样品正确率均大于85%，其中GA方法验证结果较好，达95%。通过R2、RMSEC和RMSEP分析可知，GA(R2)>MC-UVE(R2)>CARS(R2)>SPA(R2)； GA(RMSEC，RMSEP)

4 结 论

通过采集不同产地玛咖样品的近红外光谱，采用TQ 8.6软件，应用光谱标准偏差初步选择7 500～4 061 cm-1光谱波段，应用PCA-MD初步建立产地分类模型，提取光谱信息的主成分贡献率选择所提取的主成分数为5，进一步采用“模群迭代奇异样本诊断”方法通过多次循环校正诊断，剔除异常样品1—2和12—6，使用模型偏度和峰度初步判定所选光谱波段的实用性。

为使鉴定效果更好，采用MatlabR2010a分析软件，分别利用GA，CARS，MC-UVE和SPA四种方法对初步筛选的光谱波段及其分类赋值样品进行分析，优化不同产地咖样品的光谱波数，并通过RMSECV和筛选的光谱波数评估方法的可靠性。采用SIMCA-P+11.0软件结合PLS-DA对四种方法筛选的光谱波数建立判别模型，利用R2，RMSEC和RMSEP评价模型的预测性能。结果显示： 四种方法筛选的光谱波数建立的判别模型预测验证样品正确率均大于85%，其中GA方法验证结果较好，达95%，GA(R2)>MC-UVE(R2)>CARS(R2)>SPA(R2)； GA(RMSEC，RMSEP)

[1] YU Long-jiang, JIN Wen-wen, WU Yuan-xi, et al(余龙江, 金文闻, 吴元喜, 等). Natural Product Research and Development(天然产物研究与开发), 2002, 14(5): 71.

[2] YU Long-jiang, SUN You-ping, CHENG Hua, et al(余龙江, 孙友平, 程 华, 等). Acta Botanica Boreali-Occidentalia Sinica (西北植物学报), 2004, 24(10): 1901.

[3] WANG Yi-qiang, CHEN Zhang-jing, WANG Qi-ye, et al(王义强, 陈章靖, 王启业, 等). Nonwood Forest Research(经济林研究), 2014, 32(2): 167.

[4] Zheng B L, He K, Rogers L, et al. Urology, 2000, 55(4): 598.

[5] Li G, Ammermann U, Quiros C F. Economic Botany, 2001, 55(2): 255.

[6] Piacente S, Carbone V, Plaza A, et al. Journal of Agricultural and Food Chemistry, 2002, 50(20): 5621.

[7] Sandoval M, Okuhama N N, Angeles F M, et al. Food Chemistry, 2002, 79(2): 207.

[8] Onyechi O, Lawrence U S E, Akuoma H O. Nutrition Research, 1999, 19(3): 443.

[9] McKay D. Alternative Medicine Review, 2004, 9(1): 4.

[10] Schroeck F R, Hollingsworth J M, Hollenbeck B K, et al. Urology, 2013, 81(6): 1177.

[11] Meissner H O, Reich-Bilinska H, Mrozikiewicz R, et al. Menopause, 2005, 12(6): 813.

[12] ZHANG Yong-zhong, YU Long-jiang, WAN Jun-mei, et al(张永忠, 余龙江, 万军梅, 等). Natural Product Research and Development(天然产物研究与开发), 2007, 19(2): 274.

[13] Rubio J, Riqueros M I, Manuel G, et al. Food and Chemical Toxicology, 2006, 44: 1114.

[14] Brooks N A, Wilcox G, Walker K Z, et al. Menopause, 2008, 15(6): 1.

[15] Yang S H, Li G Z, Xue R G, et al. Agricultural Science and Technology, 2013, 14(12): 1877.

[16] RUAN Zhi-gang, LI Bin(阮治纲，李 彬). Chinese Journal of Pharmaceutical Analysis(药物分析杂志)，2011，31(2)： 408.

[17] Li H D, Liang Y Z, Xu Q S, et al. Analytica Chimica Acta, 2009, 648(1): 77.

[18] Cai W, Li Y, Shao X. Chemometrics and Intelligent Laboratory Systems, 2008, 90(2): 188.

[19] Han Q J, Wu H L, Cai C B, et al. Analytica Chimica Acta, 2008, 612(2): 121.

[20] Leardi R. Journal of Chemometrics, 2000, 14(5-6): 643.

[21] LIANG Yi-zeng, XU Qing-song(梁逸曾, 许青松). Instrumental Analysis of Complex Systems-White, Gray and Black Analytical Systems and Their Multivariate Methods(复杂体系仪器分析-白、灰、黑分析体系及其多变量解析方法). Beijing: Chemical Industry Press(北京: 化学工业出版社), 2012. 487.

[22] Li H D, Liang Y Z, Xu Q S, et al. Journal of Chemometrics, 2009, 24(7-8): 418.

[23] ZHAO Yan-li, ZHANG Ji, YUAN Tian-jun, et al(赵艳丽, 张 霁, 袁天军, 等). Spectroscopy and Spectral Analysis(光谱学与光谱分析), 2014, 34(7): 1831.

Study on Application of NIR Spectral Information Screening in Identification of Maca Origin

WANG Yuan-zhong, ZHAO Yan-li, ZHANG Ji, JIN Hang*

Institute of Medicinal Plants, Yunnan Academy of Agricultural Sciences, Kunming 650200, China

Medicinal and edible plant Maca is rich in various nutrients and owns great medicinal value. Based on near infrared diffuse reflectance spectra, 139 Maca samples collected from Peru and Yunnan were used to identify their geographical origins. Multiplication signal correction (MSC) coupled with second derivative (SD) and Norris derivative filter (ND) was employed in spectral pretreatment. Spectrum range (7 500～4 061 cm-1) was chosen by spectrum standard deviation. Combined with principal component analysis-mahalanobis distance (PCA-MD), the appropriate number of principal components was selected as 5. Based on the spectrum range and the number of principal components selected, two abnormal samples were eliminated by modular group iterative singular sample diagnosis method. Then, four methods were used to filter spectral variable information, competitive adaptive reweighted sampling (CARS), monte carlo-uninformative variable elimination (MC-UVE), genetic algorithm (GA) and subwindow permutation analysis (SPA). The spectral variable information filtered was evaluated by model population analysis (MPA). The results showed that RMSECV(SPA)>RMSECV(CARS)>RMSECV(MC-UVE)>RMSECV(GA), were 2.14, 2.05, 2.02, and 1.98, and the spectral variables were 250, 240, 250 and 70, respectively. According to the spectral variable filtered, partial least squares discriminant analysis (PLS-DA) was used to build the model, with random selection of 97 samples as training set, and the other 40 samples as validation set. The results showed that,R2: GA>MC-UVE>CARS>SPA, RMSEC and RMSEP: GA

Lepidium meyenii Walp.; NIR spectroscopy; Identification; Spectral information screening; Model population analysis

Jul. 11, 2014; accepted Nov. 12, 2014)

2014-07-11，

2014-11-12

国家自然科学基金项目(31460538, 81260608)和云南省自然科学基金项目(2013FD066, 2013FZ150)资助

王元忠，1981年生，云南省农业科学院药用植物研究所助理研究员 e-mail： yzwang1981@126.com *通讯联系人 e-mail： jinhang2009@126.com

O657.3

A

10.3964/j.issn.1000-0593(2016)02-0394-07

*Corresponding author