龔 圣，朱雅寧，曾陳娟，馬秀英，彭 成，郭 力*

1.成都中醫(yī)藥大學(xué)藥學(xué)院西南特色中藥資源國家重點(diǎn)實(shí)驗(yàn)室，四川 成都 611137 2.雅安三九制藥有限公司，四川 雅安 625000 3.四川佳能達(dá)攀西藥業(yè)有限公司，四川 布拖 616350

引 言

附子為毛茛科植物烏頭(AconitumcarmichaeliiDebx.)膨大側(cè)根加工品，具有回陽救逆，補(bǔ)火助陽，散寒止痛之功效，被譽(yù)為“回陽救逆第一要藥”[1]，是一種重要的、高價(jià)值的川產(chǎn)道地藥材。獨(dú)特的地理?xiàng)l件、生態(tài)環(huán)境和栽培管理技術(shù)，“江油附子”具有獨(dú)特的品質(zhì)獨(dú)特和顯著的療效，為附子上品，具有“道地藥材”和“地理標(biāo)志產(chǎn)品”的雙重美譽(yù)，其療效與經(jīng)濟(jì)效益均高于其他地區(qū)附子[2]。在過去的10多年，由于附子種植帶來良好的經(jīng)濟(jì)效益，全國多地都進(jìn)行了引種，四川布拖、陜西漢中、云南祿勸等地種植規(guī)模較大，對江油附子市場形成了沖擊。

產(chǎn)地造假是中藥材質(zhì)量監(jiān)管面臨的難題之一，特別是那些質(zhì)量與療效和產(chǎn)地有明顯關(guān)聯(lián)性的藥材被造假的風(fēng)險(xiǎn)更大。附子的療效和毒性與產(chǎn)地具有明顯的關(guān)聯(lián)性，不同產(chǎn)地附子的毒性成分相差巨大[3]。研究報(bào)道，江油附子毒性低于布拖附子與巍山附子[4]，云南附子的毒性是其18倍[4-5]。市場上仍然存在用其他產(chǎn)地附子冒充江油附子的現(xiàn)象?！暗赖厮幉摹焙汀暗乩順?biāo)志產(chǎn)品”標(biāo)簽的隨意使用不僅損害了消費(fèi)者的利益，更不利于市場監(jiān)管，甚至有可能導(dǎo)致嚴(yán)重的用藥事故。為了對道地藥材及地理標(biāo)志產(chǎn)品的聲譽(yù)保護(hù)，利于市場監(jiān)管，對附子產(chǎn)地進(jìn)行快速、有效的識別非常必要。

近紅外光譜(near-infrared spectroscopy, NIR)作為一種綠色分析技術(shù)，具有方便、快速、高效、準(zhǔn)確、低損耗等其他分析方法無法比擬的優(yōu)勢，廣泛應(yīng)用于中藥材的產(chǎn)地鑒定與質(zhì)量評價(jià)中，例如川貝母[6]、人參[7]、太子參[8]、丹參[9]、余甘子[10]、葛根[11]、三七[12]等。然而，近紅外區(qū)[780～2 500 nm(12 000～4 000 cm-1)]譜帶較寬，中藥樣品的吸收帶重疊嚴(yán)重，使得用常規(guī)方法進(jìn)行近紅外光譜分析非常困難，因此在近紅外光譜分析中，變量優(yōu)選是一個關(guān)鍵步驟。隨機(jī)森林算法[10]、遺傳算法[13]、模擬退火法[14]、區(qū)間偏最小二乘法[15]、競爭自適應(yīng)加權(quán)抽樣[16]、蒙特卡羅無信息變量消去法[17]等算法是目前常用的近紅外光譜變量篩選手段。隨機(jī)森林是處理大型和背景嘈雜數(shù)據(jù)集最有效的學(xué)習(xí)算法之一[18]。

隨機(jī)森林(random forest，RF)是由Leo Breiman提出的一種集成算法[19]，多應(yīng)用于分類問題。RF是一種由多棵弱決策樹為基分類器的集成算法(圖1)，通過組合多棵獨(dú)立的弱分類決策樹后根據(jù)投票或取均值的方式得到最終預(yù)測結(jié)果的機(jī)器學(xué)習(xí)方法，往往比單棵決策樹具有更高的準(zhǔn)確率和更強(qiáng)的穩(wěn)定性。隨機(jī)抽取樣本和樣本特征以及算法集成使RF具有出色的性能，隨機(jī)性使RF具有更穩(wěn)定的抗過擬合能力，多棵決策樹集成使RF具有更高的準(zhǔn)確率。RF隨機(jī)列抽樣的特性使RF能處理高維數(shù)據(jù)，對離散型和連續(xù)型數(shù)據(jù)能同時(shí)處理，也無需對數(shù)據(jù)做標(biāo)準(zhǔn)化處理，可以將數(shù)據(jù)缺失樣本單獨(dú)作為一類處理，因此RF對數(shù)據(jù)的格式要求較低。然而RF算法本身較為復(fù)雜，建模速度偏慢，并且隨著決策樹數(shù)目增多，訓(xùn)練模型需要更長的時(shí)間。簡而言之，RF具有處理樣本類不平衡數(shù)據(jù)、特征遺失數(shù)據(jù)和高維數(shù)據(jù)等優(yōu)勢，基于RF構(gòu)建的模型具有優(yōu)秀的抗過擬合和抗噪聲的性能[20]。這種機(jī)器學(xué)習(xí)方法已被用于疾病診斷[21]、愈后分析[22]，氣候分析[23]，種群識別[24]，公共衛(wèi)生[25]，食品摻假[26]，產(chǎn)地溯源[27]等問題，目前文獻(xiàn)中沒有關(guān)于NIR光譜采用RF在附子產(chǎn)地鑒定中的研究。本研究采集了附子三個主要栽培區(qū)域共255份樣本的NIR光譜數(shù)據(jù)，按照1∶3的比例將數(shù)據(jù)集隨機(jī)劃分為測試集和訓(xùn)練集，并分別訓(xùn)練了省級預(yù)測模型和市/縣級預(yù)測模型，采用測試隊(duì)列評估了模型的實(shí)用性和可靠性。

圖1 隨機(jī)森林算法原理

1 實(shí)驗(yàn)部分

1.1 藥材

附子樣本分別來源于四川(江油、布拖、鹽源、安縣)，云南(德欽、祿勸、虎跳峽)和陜西(南鄭、洋縣、眉縣)3省10個栽培地共255份樣本，為了確保樣品的地理真實(shí)性，研究中親自前往各栽培地采集樣本，詳細(xì)信息見表1。所有樣本經(jīng)成都中醫(yī)藥大學(xué)藥學(xué)院高繼海副教授鑒定為毛茛科烏頭屬植物AconitumcarmichaeliiDebx.的側(cè)根。樣本采集后使用超純水清洗，陶瓷刀切片，45 ℃烘箱干燥，水分被控制在10%左右，恒溫干燥箱中保存?zhèn)溆谩Ｓ糜诤罄m(xù)測定的樣本，干燥后粉碎，過80目篩，存放于干燥箱，備用。所有實(shí)驗(yàn)用的樣本存放于成都中醫(yī)藥大學(xué)藥學(xué)院郭力教授實(shí)驗(yàn)室。

表1 樣品信息

1.2 儀器與軟件

PerkinElmer Frontier傅里葉變換近紅外光譜儀(PerkinElmer公司，美國)；Mili-Q去離子水發(fā)生器(Millipore公司，美國)；101-2AB電熱鼓風(fēng)干燥箱(北京中興偉業(yè)儀器有限公司，中國)；HX高速粉碎機(jī)(中國永康溪岸五金藥具廠，中國)；80目不銹鋼篩網(wǎng)(北中西泰安公司，中國)。

SIMCA-P 13.0(Umetrics AB, 瑞典)；RStudio 4.1(RStudio Team, 2016; http://www.rstudio.com/)。

1.3 光譜采集

取預(yù)處理好的附子樣品3(±0.1)g，置于儀器配套的樣品杯中平鋪適當(dāng)壓緊。采用儀器配備的PerkinElmer Spectrum 14.0光譜采集軟件采集光譜數(shù)據(jù)，采用積分球漫反射模式，掃描范圍10 000～4 000 cm-1，分辨率4.0 cm-1，能量值2 400，累計(jì)掃描64次。掃描樣品前，儀器預(yù)熱30 min，白板扣除背景干擾，為減小實(shí)驗(yàn)誤差，每份樣品平行掃描3次，分析時(shí)取其平均光譜。在光譜采集過程中，保持室內(nèi)溫度約25 ℃。

1.4 異常值診斷

異常樣本會對模型的準(zhǔn)確性產(chǎn)生負(fù)面影響。因此，本研究基于Hotelling’s T2分布的聚類方法，95%置信限以外的樣本被視為異常值樣本[28]。經(jīng)過檢驗(yàn)，本研究中所有的樣本無異常。

1.5 光譜預(yù)處理

受到附子樣本顆粒大小不均勻、環(huán)境溫度、儀器光程等因素的影響，以及噪聲干擾、基線漂移等問題的存在，需要將得到的近紅外光譜數(shù)據(jù)進(jìn)行預(yù)處理，提高后續(xù)分析結(jié)果的準(zhǔn)確性。本研究中，Savitzky-Golay 11點(diǎn)平滑(SG)用來減低其光譜自身所攜帶的隨機(jī)誤差、提高信噪比，多元散射校正(multiplicative scatter correction, MSC)被用于消除因附子樣本顆粒分布不均勻及顆粒大小產(chǎn)生的散射，標(biāo)準(zhǔn)正態(tài)變量變換(standard normal variate transformation, SNV)用于降低由樣本顆粒大小、表面散射及光程變換對漫反射的影響[29]，小波去噪(wavelet denoising，WDS)用于得到原信號的最佳恢復(fù)，上述處理過程由SIMCA-P 13.0完成。

1.6 模型建立

隨機(jī)森林由多棵分類決策樹組成，在構(gòu)建決策樹的過程中，不進(jìn)行修剪，每棵決策樹均通過隨機(jī)選擇觀測值和變量形成。因此，隨機(jī)森林模型有兩個重要參數(shù)，ntree和mtry，它們分別決定了模型的整體分類性能、規(guī)模和單棵決策樹的情況。減小mtry值，決策樹之間的相關(guān)性和分類能力也會相應(yīng)的降低；增大mtry值，兩者也會隨之增大。ntree值由最優(yōu)mtry值確定，兩者越大模型的分類能力就越強(qiáng)，但是這會增加模型過擬合風(fēng)險(xiǎn)；過小則會降低模型的分類能力。事先通過計(jì)算袋外數(shù)據(jù)(out-of-bag, OOB)錯誤率對模型進(jìn)行內(nèi)部評估，選出最優(yōu)mtry和ntree值。

2 結(jié)果與討論

2.1 附子的原始光譜分析

各產(chǎn)地附子的原始光譜如圖2所示，樣本間的光譜整體表現(xiàn)出類似的趨勢。烏頭類生物堿是附子的主要成分之一，烏頭堿、中烏頭堿、苯甲酰烏頭堿、苯甲酰次烏頭堿等雙/單酯型二萜生物堿是附子主要的“毒-效”成分，其中的含氫基團(tuán)(O—H，N—H，C—H)振動吸收也能一定程度上反映樣本信息。4 650 cm-1附近是N-H振動合頻的吸收帶；6 000 cm-1附近波段主要為C—O和H—O一級倍頻與合頻吸收；5 200和7 000 cm-1是O—H的合頻和二倍頻的吸收帶；其中一些產(chǎn)地樣本在6 956和7 252 cm-1附近具有更強(qiáng)的吸收。就附子原始光譜而言，其譜圖間不存在明顯的區(qū)別，疊加嚴(yán)重，因此采用機(jī)器學(xué)習(xí)來探索這些數(shù)據(jù)是必要的。

圖2 原始光譜

2.2 產(chǎn)地預(yù)測模型

2.2.1 主成分分析

將原始光譜進(jìn)行主成分分析(principal component analysis, PCA)，如圖3所示。以3省[圖3(a)]和10個市/縣級產(chǎn)地[圖3(b)]為觀察樣本分別進(jìn)行PCA分析，其中兩類PCA的主成分(PCs)的貢獻(xiàn)率一致，PC1，PC2和PC3分別解釋了38.8%，26.4%和11.25%方差，前3個主成分的解釋方差達(dá)到76.4%。圖3(a)中顯示，3省的樣本在空間分布上相互疊加嚴(yán)重，且同一省份內(nèi)樣本分布也呈現(xiàn)出明顯的差異，尤其是來自四川的樣本。同樣，圖3(b)中得出結(jié)果也比較類似，各組內(nèi)樣本分布相互疊加嚴(yán)重。

圖3 不同產(chǎn)地附子主成分分析三維分值圖

2.2.2 最佳預(yù)處理方式

不同的預(yù)處理方式對模型的預(yù)測準(zhǔn)確率存在顯著影響，因此需要篩選出最佳的預(yù)處理方式，基于RF模型默認(rèn)參數(shù)(ntree=500，mtry=2)比較了14種預(yù)處理方式下模型的預(yù)測準(zhǔn)確率，如圖4所示，SG+MSC和SG+MSC+SNV模式都具有較佳的準(zhǔn)確率，因此后續(xù)分析選擇SG+MSC處理方式處理的NIR數(shù)據(jù)。需要明確的是MSC和SNV屬于同一類數(shù)據(jù)預(yù)處理方法，同時(shí)使用MSC和SNV模式，有過度處理數(shù)據(jù)的可能性，還需要在后續(xù)的實(shí)際研究中探討。

圖4 不同預(yù)處理方式對模型準(zhǔn)確率的影響

2.2.3 RF模型的建立

在默認(rèn)參數(shù)(ntree=500，mtry=2)條件下，使用訓(xùn)練數(shù)據(jù)集初步探索了RF模型對3省樣本的袋外數(shù)據(jù)錯誤率(OOB Error), 模型的決策樹數(shù)目在230～500時(shí)，模型的袋外數(shù)據(jù)犯錯率趨于穩(wěn)定[圖5(a)]。隨后，使用“randomForest”軟件包中“tuneRF”函數(shù)測試了1～100內(nèi)的最佳mtry值，當(dāng)mtry為4時(shí)，模型整體誤判率為最低2.58%[圖5(b)]。

樣本按照3∶1的比例被隨機(jī)劃分為訓(xùn)練集和測試集兩組，進(jìn)行建模分析(ntree=500，mtry=4)，這個過程主要由RStudio 4.1軟件中的“randomForest”和“caret”包完成。

靈敏度(Sensitivity)、特異度(Specificity)、陽性預(yù)測值(Pos Pred Value)、陰性預(yù)測值(Neg Pred Value)、Prevalence、Detection Rate、Detection Prevalence和平衡精度(Balanced Accuracy)等是用來評判模型的性能常用指標(biāo)，本研究中，上述指標(biāo)值由“caret”包自動計(jì)算得出，針對多分類模型時(shí)，將采用加權(quán)計(jì)算上述參數(shù)。其中前4個指標(biāo)都主要反映了模型對分類樣本的識別能力，平衡精度也是一個重要參數(shù)，通常其值越接近1，模型性能越佳。

圖5 最佳ntree值(a)和mtry值(b)

2.2.4 RF模型對省級區(qū)域的識別結(jié)果

首先，使用訓(xùn)練集樣本建立了四川(Sichuan)、云南(Yunnan)和陜西(Shaanxi)省級識別模型[圖6(a)]，并采用測試集樣本評估了模型性能[圖6(b)]。194份訓(xùn)練隊(duì)列樣本僅5例樣本被錯誤劃分(97.42%)，62份測試樣本有61份被準(zhǔn)確預(yù)測(98.39%)。模型對3省樣本的靈敏度和特異度都達(dá)到了1.000 0，平衡精度也分別達(dá)到了0.982 8, 0.988 4和1.000 0[圖6(c)]。由此可見，僅采用全波長數(shù)據(jù)訓(xùn)練省級預(yù)測模型，仍能取得較高的預(yù)測率，在后續(xù)的應(yīng)用中也推薦使用全波長NIR數(shù)據(jù)對附子的大區(qū)域來源進(jìn)行預(yù)測。事實(shí)上，這3個省的土壤、氣候、地理環(huán)境、海拔等因素差異較大，是導(dǎo)致附子的生長代謝等方面出現(xiàn)差異的主要因素。此外，附子在各地的栽培方法、采收時(shí)間等管理因素也有所差異，也是導(dǎo)致附子各成分的積累出現(xiàn)差異的原因之一。

2.2.5 RF模型對市/縣級區(qū)域的識別結(jié)果

按照2.2.3模型，采用訓(xùn)練集樣本全波長數(shù)據(jù)重新訓(xùn)練了一個新的市/縣級模型，探索了最佳ntree和mtry值，分別為500和4，此時(shí)模型犯錯率最低為2.05%，結(jié)果如圖6(d, e)所示。模型對大部分產(chǎn)地樣本的靈敏度和特異度都達(dá)到了1.000 0，除了來自布拖(0.625 0, 0.962 3)和鹽源(0.666 7, 0.945 5)的樣本，平衡精度僅有0.793 6和0.806 1。然而布拖和鹽源在地理上相隔較近，種植基地都處于高海拔區(qū)域(>2 500 m)，且種植技術(shù)比較一致(定點(diǎn)扶貧縣，得到了四川好醫(yī)生攀西藥業(yè)的技術(shù)支持)，導(dǎo)致兩個栽培地之間的樣本差異降低，這在一定程度上增大了分類難度。

由于NIR數(shù)據(jù)中存在許多無效數(shù)據(jù)，導(dǎo)致RF模型運(yùn)行時(shí)間過長，根據(jù)隨機(jī)森林模型對各波數(shù)的重要性評分[基于Gini指數(shù)，圖7(a)]，優(yōu)選出前100個重要特征波數(shù)[圖7(b)]，重新訓(xùn)練模型。結(jié)果顯示，優(yōu)化后的RF模型對布拖(BT)和鹽源(YY)的靈敏度(0.750 00，1.000 00)、特異度(1.000 00，0.963 64)以及平衡精度(0.875 00，0.981 82)得到顯著提升[圖7(c, d)]，對其余產(chǎn)地附子的判別性能則無明顯改善。

圖6 RF模型對省級樣本訓(xùn)練集(a)和測試集(b)的預(yù)測結(jié)果和識別性能(c)；RF模型對市/縣級測試隊(duì)列樣本(d)的預(yù)測結(jié)果和識別性能(e)

圖7 (a)特征波數(shù)重要性評分；(b)前100個重要波數(shù)；(c)優(yōu)化RF模型對BT和YY樣本的判別性能；(d)優(yōu)化的RF模型

在來自四川的樣本內(nèi)，本研究選取了江油附子的種源地之一安縣的樣本，這兩地在地理間隔上要比布拖和鹽源更近。有趣的是，模型對江油和安縣的樣本表現(xiàn)出了100%的識別能力。除此之外，對于來自云南和陜西各地的樣本，模型同樣表現(xiàn)優(yōu)異，高海拔與低海拔之間的樣本得到了完全劃分。盡管都是作為附子主要產(chǎn)地之一，云南多數(shù)附子種植在1 500 m海拔上的地區(qū)，而來自陜西漢中和寶雞的樣本則多數(shù)來源于低海拔(<1 500 m)的地區(qū)。根據(jù)林俊芝[30]和張定堃[4]等的報(bào)道，安縣附子中新烏頭堿、烏頭堿、展花烏頭寧、塔拉烏頭胺、苯甲酰新烏頭原堿、苯甲酰次烏頭原堿、脫水苯甲酰新烏頭原堿及新烏頭原堿含量較江油附子更高，而江油附子中尼奧林、附子靈、次烏頭堿、塔拉他定、卡米查林、附子亭、宋果靈含量更高，這種差異的形成與兩地的栽培條件差異有關(guān)。采用這些化合物的含量差異還可以用于區(qū)分布拖附子和巍山附子。

盡管NIR數(shù)據(jù)主要反映了化合物官能團(tuán)信息，一定程度上這種生物堿含量差異在光譜信息中也得到體現(xiàn)，利用這些差異使NIR在中藥材質(zhì)量監(jiān)控、成分預(yù)測以及在線評估等領(lǐng)域得到應(yīng)用。本研究只初步探索了NIR集成機(jī)器學(xué)習(xí)在附子主要的3大種植生產(chǎn)區(qū)域的應(yīng)用，而我國多數(shù)地區(qū)都是附子適宜引種區(qū)域[1]。各地氣候、環(huán)境、種植和管理等因素差異較大，在一定程度上影響植物的生長發(fā)育和代謝物積累，同時(shí)也會影響微量元素和元素同位素的富集[31]，從而給全國附子產(chǎn)地識別、精準(zhǔn)溯源等方面提供了可能。

3 結(jié) 論

(1)通過比較不同預(yù)處理方式和優(yōu)化RF模型參數(shù)，建立了兩個針對不同尺度范圍的附子產(chǎn)地識別模型，SG+MSC相結(jié)合為最佳的預(yù)處理模式，ntree=500，mtry=4，RF模型達(dá)到最佳性能。

(2)對于省級樣本，僅采用原始光譜數(shù)據(jù)RF模型同樣表現(xiàn)出了優(yōu)異的性能，訓(xùn)練集和測試集的準(zhǔn)確率分別為92.78%和91.94%，預(yù)處理后模型的準(zhǔn)確率達(dá)到98.39%。對于市/縣級樣本，RF模型不僅能準(zhǔn)確識別道地產(chǎn)區(qū)江油附子，而且對其他栽培地附子的識別率超過了85%。

(3)后續(xù)研究中可以通過測定樣本內(nèi)硫、鍶、鉛、硼元素的穩(wěn)定同位素豐度比值(δ34S，δ87Sr，δ207Pb，δ208Pb和δ11B)，獲得各栽培地樣本間元素同位素豐度差異信息，來增強(qiáng)模型的穩(wěn)定性和可靠性。