丁 鴻，徐芳芳，杜 慧，張 欣,3，徐 冰，吳 云,3，王振中,3，肖 偉,3*

1.南京中醫(yī)藥大學(xué)，江蘇 南京 210023

2.江蘇康緣藥業(yè)股份有限公司，江蘇 連云港 222001

3.中藥制藥過(guò)程新技術(shù)國(guó)家重點(diǎn)實(shí)驗(yàn)室，江蘇 連云港 222001

4.北京中醫(yī)藥大學(xué) 中藥信息學(xué)系，北京 102400

熱毒寧注射液是由金銀花、青蒿和梔子3 種中藥提取精制而成[1-2]，其前提取精制過(guò)程包括金銀花提取濃縮、青蒿提取濃縮、梔子提取濃縮、金銀花青蒿（金青）醇沉、金青萃取、梔子萃取等工序，均在中藥數(shù)字化提取精制工廠中完成[3]。中藥數(shù)字化提取工廠配備集散控制系統(tǒng)（distributed control system，DCS）、車間制造執(zhí)行系統(tǒng)（manufacturing execution system，ΜES）等，可實(shí)現(xiàn)對(duì)中藥前提取設(shè)備傳感器實(shí)時(shí)數(shù)據(jù)以及生產(chǎn)過(guò)程數(shù)據(jù)的抓取。在工業(yè)數(shù)據(jù)的挖掘與分析中，機(jī)器學(xué)習(xí)是常用的方法之一[4]。杜慧等[5]以金銀花和青蒿濃縮浸膏醇沉工序?yàn)檠芯繉?duì)象，采用偏最小二乘法（partial least squares，PLS）進(jìn)行數(shù)據(jù)挖掘，PLS 需要對(duì)預(yù)測(cè)變量和目標(biāo)變量進(jìn)行降維投影處理[6]，不利于直觀解釋生產(chǎn)規(guī)律，該研究篩選了關(guān)鍵變量但未劃分變量的控制范圍，不便于生產(chǎn)的實(shí)際調(diào)控。決策樹(shù)算法是機(jī)器學(xué)習(xí)中一個(gè)樹(shù)狀預(yù)測(cè)模型，適用于多階段決策、表查找和最優(yōu)化等問(wèn)題，能夠直觀解釋決策過(guò)程，目前已成功應(yīng)用于醫(yī)療、廣告、交通和金融等多個(gè)領(lǐng)域[7-8]。本實(shí)驗(yàn)采用決策樹(shù)及其相關(guān)算法對(duì)歷史數(shù)據(jù)進(jìn)行分析，深入挖掘并尋找關(guān)鍵變量的優(yōu)化控制范圍，期望提升熱毒寧注射液的質(zhì)量控制水平。

1 材料

建模數(shù)據(jù)由江蘇康緣藥業(yè)股份有限公司數(shù)字化提取工廠提供，金青中間體批次為170105 至180905共205 批；數(shù)據(jù)的訓(xùn)練集和測(cè)試集的劃分應(yīng)用Python 的Scikit-Learn 工具實(shí)現(xiàn)。建模算法均在Salford Predictive Μodeler 8.3（美國(guó)Μinitab 公司）實(shí)現(xiàn)，單值過(guò)程控制圖在Μinitab 19（美國(guó)Μinitab公司）繪制。

2 方法

2.1 目標(biāo)變量的確定

金青醇沉工段由金青浸膏分配、金青浸膏醇沉、上清液暫存和上清液回收4 個(gè)過(guò)程組成。金銀花和青蒿的水提濃縮浸膏通過(guò)稱量罐稱定質(zhì)量后分配至2 個(gè)醇沉罐（醇沉罐1 和醇沉罐2），在醇沉過(guò)程中，向分配浸膏的2 個(gè)醇沉罐中分3 個(gè)階段按不同的體積流量加入一定量的乙醇并不斷攪拌，攪拌停止靜置一段時(shí)間并將醇沉上清液通過(guò)真空泵轉(zhuǎn)入濃縮罐，刮板濃縮后得金青醇沉浸膏。金青醇沉工段簡(jiǎn)圖如圖1所示，各過(guò)程主要參數(shù)如表1所示。前期辨識(shí)出金青醇沉工段的關(guān)鍵因素為金青醇沉上清液和上清液傳料體積、金青總質(zhì)量、醇沉罐1 及醇沉罐2 金銀花浸膏質(zhì)量[5]。上清液傳料體積為上清液暫存罐傳送至濃縮罐的上清液體積；醇沉上清液為上清液暫存罐收集的上清液體積；金青浸膏總質(zhì)量為2 個(gè)醇提罐加入金銀花和青蒿總質(zhì)量；醇沉罐1和醇沉罐2 金銀花質(zhì)量為分配至2 個(gè)醇沉罐的金銀花質(zhì)量。所篩選的參數(shù)均為金青浸膏醇沉過(guò)程輸入物料及輸出物料的質(zhì)量和體積參數(shù)，醇沉過(guò)程工藝較復(fù)雜，涉及影響因素較多，為進(jìn)一步挖掘潛在規(guī)律，本研究選擇金青浸膏醇沉過(guò)程為研究目標(biāo)建立過(guò)程模型。

表1 金青醇沉4 個(gè)過(guò)程的主要參數(shù)Table 1 Main parameters of four processes of LA alcohol precipitation

圖1 金青醇沉回收工段流程簡(jiǎn)化圖Fig.1 Data map for LA alcohol precipitation process

為消除過(guò)程物料質(zhì)量體積相關(guān)性對(duì)模型的影響，本研究擬以浸膏醇沉過(guò)程收率為目標(biāo)變量。將金青浸膏醇沉收率和金青醇沉工段收率進(jìn)行相關(guān)性分析，相關(guān)圖如圖2所示，相關(guān)系數(shù)（r）為0.653。控制浸膏醇沉過(guò)程對(duì)提高金青醇沉工段的生產(chǎn)效率十分重要。綜上，本實(shí)驗(yàn)以金青浸膏醇沉收率作為目標(biāo)變量，建立浸膏醇沉過(guò)程模型。

圖2 金青浸膏醇沉收率與金青醇沉工段收率相關(guān)圖Fig.2 Correlation diagram between yield of LA extract alcohol precipitation and yield of LA alcohol precipitation section

2.2 建模數(shù)據(jù)準(zhǔn)備

本實(shí)驗(yàn)建模數(shù)據(jù)主要來(lái)源于批記錄及數(shù)字化提取工廠數(shù)據(jù)庫(kù)集成數(shù)據(jù)。批記錄為手動(dòng)錄入，主要為原料藥、中間體物料相關(guān)數(shù)據(jù)；數(shù)據(jù)庫(kù)收集的數(shù)據(jù)主要從DCS 系統(tǒng)中收集各個(gè)設(shè)備傳感器實(shí)時(shí)生產(chǎn)過(guò)程中抓取的時(shí)序數(shù)據(jù)，抓取頻率為每分鐘1 次。中藥生產(chǎn)環(huán)境復(fù)雜，不同時(shí)間點(diǎn)內(nèi)時(shí)序參數(shù)相關(guān)性強(qiáng)，前期對(duì)DCS 數(shù)據(jù)進(jìn)行數(shù)據(jù)清洗和特征提取后得到建模特征參數(shù)[6]。本實(shí)驗(yàn)在建模特征參數(shù)中選擇30 個(gè)醇沉罐傳感器數(shù)據(jù)參數(shù)納入建模數(shù)據(jù)集，考慮到浸膏質(zhì)量和加醇體積間存在較強(qiáng)的相關(guān)性，本實(shí)驗(yàn)將金銀花浸膏質(zhì)量、青蒿浸膏質(zhì)量和加醇量轉(zhuǎn)化為金青浸膏質(zhì)量比和加醇料液比。在中藥制劑的過(guò)程中，中藥材的質(zhì)量差異會(huì)傳遞至處方藥味、中間體及成品，直接影響中藥制劑批間質(zhì)量的穩(wěn)定[9]。為考察藥材及金銀花浸膏及青蒿浸膏的各生產(chǎn)階段可能對(duì)醇沉過(guò)程的影響，收集藥材質(zhì)量參數(shù)、藥材廠商及藥材至濃縮浸膏各過(guò)程工序的收率作為預(yù)測(cè)變量，樣本收集時(shí)間為2017年1月—2018年9月，共205 批歷史數(shù)據(jù)。以批號(hào)關(guān)聯(lián)建模參數(shù)，建模參數(shù)共49 個(gè)，最終得到205×49 的建模數(shù)據(jù)集，建模變量如表2所示。

表2 金青浸膏醇沉過(guò)程49 建模變量Table 2 Fourty-nine modeling variables in LA extract alcohol precipitation process

2.3 建模算法

決策樹(shù)（decision tree）是一種基本的分類和回歸方法，該學(xué)習(xí)方法主要包括特征選擇、決策樹(shù)生成和修剪，目前決策樹(shù)學(xué)習(xí)常用的算法有ID3、C4.5和分類與回歸樹(shù)（classification and regression tree，CART）[10]。當(dāng)單一決策樹(shù)學(xué)習(xí)器不能很好地處理問(wèn)題時(shí)，可嘗試集成學(xué)習(xí)方法（ensemble learning）提高模型預(yù)測(cè)性能。集成學(xué)習(xí)通過(guò)訓(xùn)練多個(gè)學(xué)習(xí)器并進(jìn)行組合來(lái)解決同一問(wèn)題[11]。Bagging 算法和Boosting 算法是2 種經(jīng)典的集成學(xué)習(xí)算法[12]。本實(shí)驗(yàn)采用CART 算法以及隨機(jī)森林、TreeNet 2 種集成樹(shù)算法進(jìn)行建模。

2.3.1 CART CART 模型是應(yīng)用廣泛的決策樹(shù)算法[13]，它是由特征選擇、樹(shù)的生成及剪枝組成，既可以用于分類，也可以用于回歸[10]。相較于ID3，CART 算法選擇基尼不純度作為度量標(biāo)準(zhǔn)，并采用二分遞歸分割方法劃分?jǐn)?shù)據(jù)集，最終生成二叉決策樹(shù)，二叉劃分適用更加廣泛，也更易處理數(shù)值型數(shù)據(jù)[6]。

2.3.2 隨機(jī)森林（random forests，RF） RF 是由多個(gè)Bagging 算法訓(xùn)練的決策樹(shù)組合而成[14]，在解決回歸問(wèn)題時(shí)，其以多個(gè)決策樹(shù)的平均值為輸出結(jié)果[12,15]。它是通過(guò)利用bootstrap 重抽樣方法從原始樣本中抽取多個(gè)樣本，沒(méi)有抽到的樣本被稱為帶外數(shù)據(jù)（out-of-bag，OOB）[16]。OOB 誤差估計(jì)是RF算法的一個(gè)無(wú)偏估計(jì)，可替代數(shù)據(jù)集的交叉驗(yàn)證，明顯降低計(jì)算復(fù)雜性[12,17]。

2.3.3 TreeNet TreeNet 被稱為隨機(jī)梯度提升（stochastic gradient boosting）[18]。梯度提升是一種實(shí)現(xiàn)Boosting 的方法，其主要思想是每次建模都在前一模型損失函數(shù)（loss function）梯度下降方向進(jìn)行[19]。TreeNet 在訓(xùn)練的過(guò)程中首先生成1 棵小樹(shù)并計(jì)算該模型的預(yù)測(cè)殘差，之后建立第2 棵樹(shù)對(duì)第1 棵樹(shù)的殘差進(jìn)行預(yù)測(cè)，第3 棵樹(shù)對(duì)前2 棵樹(shù)的殘差進(jìn)行預(yù)測(cè)。上步驟遞歸，新產(chǎn)生的小樹(shù)對(duì)前面的模型進(jìn)行不斷修正，最終得到最佳模型。

2.4 數(shù)據(jù)集劃分與訓(xùn)練

數(shù)據(jù)按3∶1 的比例隨機(jī)劃分為訓(xùn)練集和測(cè)試集。通過(guò)訓(xùn)練集的數(shù)據(jù)來(lái)訓(xùn)練數(shù)據(jù)，用測(cè)試集來(lái)評(píng)估模型。在訓(xùn)練模型的過(guò)程中需要訓(xùn)練不同模型參數(shù)的模型來(lái)篩選最佳模型，將保留訓(xùn)練集的一部分來(lái)評(píng)估各個(gè)模型的性能，新的保留集稱為驗(yàn)證集。本實(shí)驗(yàn)采用10 折交叉驗(yàn)證來(lái)評(píng)估CART 和TreeNet的訓(xùn)練模型，RF 為OOB 誤差估計(jì)。本實(shí)驗(yàn)以均方根誤差（root mean square error，RΜSE）、平均絕對(duì)偏差（mean absolute difference，ΜAD）以及決定系數(shù)（R2）作為模型評(píng)價(jià)指標(biāo)。RΜSE 及ΜAD 越接近0，R2越接近1 表示模型誤差越小，性能越好。

n為校正集或驗(yàn)證集的樣本數(shù)，t∈[1，n]，y為參考值，yi為預(yù)測(cè)值，為所有樣品參考值的平均值

3 結(jié)果與分析

3.1 各模型性能評(píng)價(jià)

基于訓(xùn)練集數(shù)據(jù)分別建立CART、RF和TreeNet模型，并采用10 折交叉驗(yàn)證來(lái)評(píng)估不同超參數(shù)的CART 和TreeNet 訓(xùn)練模型性能，OOB 誤差估計(jì)評(píng)估RF 模型性能。最終得到各模型最佳的超參數(shù)設(shè)置如表3所示，最優(yōu)模型的性能評(píng)價(jià)指標(biāo)如表4所示，各模型的預(yù)測(cè)變量重要性如表5所示。由表5中可知，與CART 算法相比，RF 和TreeNet 2 種集成樹(shù)算法的驗(yàn)證集和預(yù)測(cè)集均方根更小，決定系數(shù)更大，模型效果更好。說(shuō)明集成學(xué)習(xí)算法有效地提高了單一決策樹(shù)算法模型的性能。

表3 3 種模型的超參數(shù)設(shè)置Table 3 Hyperparameter setting table of three models

表4 3 種最優(yōu)模型的性能指標(biāo)評(píng)價(jià)Table 4 Performance index evaluation of three models

由表5可知，3 種模型的共同重要變量為X16（罐1 料液比）、X34（罐2 料液比）和X5（金銀花提取液濃縮收率）。以這3 個(gè)變量建立模型，其模型效果如表6所示。

表5 CART、RF 和TreeNet 模型變量重要性Table 5 Importance chart of predicted table of CART,RF and TreeNet

表6 篩選后3 種模型的性能指標(biāo)評(píng)價(jià)Table 6 Performance index evaluation of three models after screening

由變量篩選后模型性能觀察比較可知，變量篩選后測(cè)試集和預(yù)測(cè)集的模型誤差及R2變化較小，均保持較好的模型性能。說(shuō)明醇沉罐1 料液比、醇沉罐2 料液比以及金銀花提取液濃縮收率為金青浸膏醇沉過(guò)程的關(guān)鍵變量。

3.2 關(guān)鍵變量依存性分析

為進(jìn)一步描述3 個(gè)關(guān)鍵變量與目標(biāo)變量的潛在規(guī)律。對(duì)TreeNet 模型的3 個(gè)關(guān)鍵變量的變量依存度圖（partial dependency plots，PDPs）進(jìn)行分析，單變量依存性散點(diǎn)圖如圖3所示。

圖3 關(guān)鍵變量為金銀花提取液濃縮收率 (A)、罐1 料液比 (B) 及罐2 料液比 (C) 單變量依存度散點(diǎn)圖Fig.3 PDPs of univariate dependency with key variable were concentration yield of honeysuckle extract liquid (A) and solidliquid ratios in extraction tank 1 (B) and that in tank 2 (C)

通過(guò)對(duì)3 個(gè)關(guān)鍵變量的PDPs 分析，各變量與其部分依存度不是簡(jiǎn)單的線性關(guān)系。變量依存度分析以響應(yīng)值的正負(fù)表示對(duì)轉(zhuǎn)化率影響的好壞，將工藝優(yōu)化問(wèn)題轉(zhuǎn)化為計(jì)算問(wèn)題。當(dāng)金銀花提取液濃縮收率小于0.039 kg/L 時(shí)，目標(biāo)變量預(yù)測(cè)響應(yīng)值為正值。醇沉罐1 和醇沉罐2 料液比的依存度隨料液比的不斷增加的變化可分為3 個(gè)范圍：料液比較低，當(dāng)料液比較小時(shí)依存度為負(fù)值，較小的料液比將降低目標(biāo)變量的預(yù)測(cè)響應(yīng)值；料液比居中，隨著料液比的不斷增加依存度也逐漸提高并轉(zhuǎn)為正值；料液比較高，當(dāng)增加到一定范圍時(shí)，料液比的增加對(duì)目標(biāo)變量預(yù)測(cè)響應(yīng)值的影響不變。將料液比控制在合適區(qū)間，在提高金青浸膏醇沉的收率同時(shí)可避免加入過(guò)量乙醇，降低生產(chǎn)成本。

模型結(jié)果中醇沉罐1 的料液比控制范圍大約為5.54～5.91 L/kg，醇沉罐2 的范圍大約為5.58～5.91 L/kg。對(duì)2 罐的范圍進(jìn)行中位數(shù)檢驗(yàn)，分析2 罐間是否存在差異，以便于確定是否需要分罐控制。結(jié)果如表7所示。

表7 不同提取罐料液比中位數(shù)Table 7 Median of solid-liquid ratios in different extraction tanks

醇沉罐1 和醇沉罐2 相比，P值＞0.05，總體中位數(shù)相等，醇沉罐1 和醇沉罐2 間料液比優(yōu)化范圍無(wú)顯著差異。綜上，醇沉罐1 和醇沉罐2 可統(tǒng)一控制，為控制金青浸膏醇沉過(guò)程的收率，需主要對(duì)金銀花提取液濃縮收率和各罐料液比進(jìn)行控制。

3.3 關(guān)鍵變量過(guò)程控制

上文篩選影響金青浸膏醇沉的關(guān)鍵變量為金銀花提取液濃縮收率和料液比，金銀花醇沉上清液濃縮收率的影響因素較多，涉及上一工段各工藝參數(shù)及設(shè)備過(guò)程參數(shù)，需另外建模研究。本實(shí)驗(yàn)主要對(duì)料液比的過(guò)程控制進(jìn)行研究。

控制料液比主要控制金銀花浸膏、青蒿浸膏和加醇量。通過(guò)收集落入優(yōu)化范圍內(nèi)外的金銀花浸膏、青蒿浸膏和加醇量，進(jìn)行中位數(shù)檢驗(yàn)分析差異。結(jié)果如表8所示。統(tǒng)計(jì)分析表明，影響料液比的變化主要是由于金銀花浸膏質(zhì)量和加醇量的波動(dòng)引起的。企業(yè)內(nèi)部控制標(biāo)準(zhǔn)中要求金銀花浸膏濃縮至1.10～1.12 g/cm3（70～80 ℃），并分別對(duì)密度為1.10、1.11、1.12 g/cm3的浸膏的出膏量給出不同的控制范圍。依存度分析優(yōu)化范圍內(nèi)的金銀花浸膏出膏溫度均符合要求，浸膏密度分為1.11、1.12 g/cm3。分別對(duì)密度為1.11、1.12 g/cm3的金銀花浸膏做控制圖，控制圖如圖4所示。

圖4 密度為1.11 g·cm-3 金銀花浸膏質(zhì)量 (A) 與加醇量 (B) 及密度為1.12 g·cm-3 金銀花浸膏質(zhì)量 (C) 與加醇量 (D) 的過(guò)程控制圖Fig.4 Process control diagram of quality of honeysuckle extract with a density of 1.11 g·cm-3 (A),amount of alcohol (B),quality of honeysuckle extract with a density of 1.12 g·cm-3 (C),and amount of alcohol (D)

表8 不同料液比的相關(guān)變量中位數(shù)Table 8 Median of related variables for different solidliquid ratios

如控制圖所示，通過(guò)篩選得到密度為1.11 的金銀花浸膏分配控制范圍為557.92～639.62 kg，加醇量的控制范圍為3.370～3.828 m3；密度為1.12 的金銀花浸膏的控制范圍為540.4～616.9 kg，加醇量的控制范圍為3.317～3.859 m3。綜上分析，優(yōu)化范圍內(nèi)加醇量的差異不大，金銀花浸膏質(zhì)量的波動(dòng)是影響金青醇沉過(guò)程的重要因素。需加強(qiáng)對(duì)金銀花浸膏質(zhì)量的監(jiān)控手段和方法，將金銀花浸膏質(zhì)量控制在穩(wěn)定的優(yōu)化范圍內(nèi)控制加醇量，從而控制醇沉過(guò)程高效穩(wěn)定進(jìn)行。

4 總結(jié)與展望

醇沉工序是中藥提取精制過(guò)程的關(guān)鍵工序之一，在實(shí)際生產(chǎn)控制過(guò)程中對(duì)金青醇沉過(guò)程控制的穩(wěn)定性還有待提高。本實(shí)驗(yàn)依托熱毒寧提取精制生產(chǎn)大數(shù)據(jù)，收集金青醇沉過(guò)程歷史數(shù)據(jù)，以金青浸膏醇沉上清液收率為目標(biāo)變量，基于決策樹(shù)算法進(jìn)行建模分析，綜合3 個(gè)模型建模結(jié)果，篩選出重要變量為金銀花上清液濃縮收率和醇沉罐料液比；通過(guò)變量依存度分析，優(yōu)選批次并作金銀花質(zhì)量及加醇量過(guò)程控制圖優(yōu)化控制范圍。從熱毒寧注射液工業(yè)數(shù)據(jù)中挖掘質(zhì)量傳遞規(guī)律，有利于該產(chǎn)品質(zhì)量的追溯與控制，促進(jìn)制造過(guò)程的優(yōu)化及持續(xù)改進(jìn)[17]。

決策樹(shù)算法在中藥工業(yè)大數(shù)據(jù)的建模應(yīng)用報(bào)道較少，本實(shí)驗(yàn)首次嘗試應(yīng)用于熱毒寧注射液的醇沉過(guò)程，相較于PLS 算法，決策樹(shù)算法在缺失值處理、變量篩選及優(yōu)化參數(shù)范圍的劃分等方面更具優(yōu)勢(shì)，為熱毒寧注射液的數(shù)據(jù)挖掘和質(zhì)量控制提供新路徑。本實(shí)驗(yàn)建模數(shù)據(jù)為生產(chǎn)歷史數(shù)據(jù)，在后續(xù)的工作中，還需結(jié)合近紅外快速檢測(cè)技術(shù)收集中間體含量數(shù)據(jù)，多目標(biāo)決策優(yōu)選工藝參數(shù)控制范圍。此外，還需持續(xù)收集生產(chǎn)數(shù)據(jù)，不斷優(yōu)化及驗(yàn)證工藝參數(shù)控制范圍，提升金青醇沉過(guò)程的控制水平。

利益沖突所有作者均聲明不存在利益沖突