基于GSA的厭氧發(fā)酵原料碳氮比NIRS快速檢測

劉金明 程秋爽 甄 峰 許永花 李文哲,5 孫 勇,5

(1.東北農(nóng)業(yè)大學(xué)工程學(xué)院， 哈爾濱 150030； 2.黑龍江八一農(nóng)墾大學(xué)電氣與信息學(xué)院， 大慶 163319；3.中國科學(xué)院廣州能源研究所中國科學(xué)院可再生能源重點實驗室， 廣州 510640；4.東北農(nóng)業(yè)大學(xué)電氣與信息學(xué)院， 哈爾濱 150030；5.黑龍江省寒地農(nóng)業(yè)可再生資源利用技術(shù)及裝備重點實驗室， 哈爾濱 150030)

0 引言

碳氮比(碳、氮元素質(zhì)量的比值)對厭氧發(fā)酵微生物的生長繁殖和產(chǎn)物合成具有重要影響，測定厭氧發(fā)酵底物的碳氮比已成為優(yōu)化厭氧發(fā)酵原料配比的重要環(huán)節(jié)[1-2]。在以玉米秸稈為主要原料進行厭氧發(fā)酵生產(chǎn)沼氣時，通過預(yù)處理打破玉米秸稈自身緊密的木質(zhì)纖維素結(jié)構(gòu)，能夠有效提高玉米秸稈的生物轉(zhuǎn)化利用率[3-4]。同時，由于玉米秸稈碳氮比過高，常與碳氮比較低的畜禽糞便混合共同發(fā)酵，以提高厭氧發(fā)酵產(chǎn)沼氣的效率和能力[5-6]。為了分析碳氮比對厭氧發(fā)酵的影響，進而對預(yù)處理后玉米秸稈及秸稈和糞便混合物的厭氧發(fā)酵過程進行有效調(diào)控，有必要對發(fā)酵原料的碳氮比進行快速、準(zhǔn)確的測定，采用傳統(tǒng)化學(xué)方法測定其碳氮比時存在測試速度慢、成本高的問題。

近紅外光譜(Near infrared spectroscopy, NIRS)分析技術(shù)具有快速、低成本及多組分同步檢測等優(yōu)點[7-8]，已廣泛用于有機動植物廢棄物碳氮含量的快速測定[9-10]。當(dāng)使用NIRS對畜禽糞便中的碳、氮含量進行檢測時，氮含量的檢測精度明顯優(yōu)于碳含量[11-12]。在使用NIRS對植物中所含碳、氮成分進行檢測時，能夠獲得較高的氮含量檢測精度[13-14]，但碳含量檢測精度較低[15-16]。針對厭氧發(fā)酵過程中原料碳氮比的快速檢測需求，以及NIRS在秸稈和糞便碳氮含量檢測方面的優(yōu)勢與不足，本文提出使用NIRS對厭氧發(fā)酵原料直接進行碳氮比的快速檢測。

針對以全譜建模時冗余波長點嚴(yán)重影響模型檢測精度和效率的問題，相關(guān)學(xué)者提出應(yīng)用遺傳算法(Genetic algorithm，GA)[17-18]進行NIRS波長變量組合優(yōu)化。GA因其具有較強的魯棒性和全局搜索能力，在NIRS特征波長優(yōu)選方面得到了廣泛應(yīng)用，其隨機搜索能力能夠有效解決光譜波長點之間的共線性問題[19]，還可以與其他光譜譜區(qū)優(yōu)化算法相結(jié)合進行特征波長點的優(yōu)選[20-21]。但GA存在早熟問題，且進化后期搜索效率較低。遺傳模擬退火算法(Genetic simulated annealing algorithm, GSA)[22]通過將GA與模擬退火算法(Simulated annealing algorithm, SA)相結(jié)合，引入溫度參數(shù)對適應(yīng)度函數(shù)進行改進設(shè)計，在充分發(fā)揮算法強大搜索能力的同時，有效解決了GA存在的不足。

本文提出基于GSA構(gòu)建遺傳模擬退火區(qū)間偏最小二乘算法(Genetic simulated annealing interval partial least squares algorithm, GSA-iPLS)和雙重遺傳模擬退火偏最小二乘算法(Double genetic simulated annealing partial least squares algorithm, DGSA-PLS)分別用于特征譜區(qū)優(yōu)選和特征波長點優(yōu)選，進而獲取與碳氮比相關(guān)性高的有效特征波長變量，建立厭氧發(fā)酵原料碳氮比的NIRS快速檢測模型。

1 材料與方法

1.1 樣品采集與制備

實驗用玉米秸稈取自東北農(nóng)業(yè)大學(xué)校內(nèi)實驗田，豬糞取自哈爾濱市三元畜產(chǎn)實業(yè)公司，牛糞取自哈爾濱市宇峰奶牛養(yǎng)殖農(nóng)民專業(yè)合作社，羊糞取自東北農(nóng)業(yè)大學(xué)阿城實驗實習(xí)基地。玉米秸稈、豬糞、牛糞和羊糞各采集1個樣品，用于后續(xù)秸稈預(yù)處理和秸稈糞便混合樣品制備。采集的玉米秸稈樣品自然風(fēng)干后裝袋保存，采集的豬糞、牛糞和羊糞樣品于-18℃冷凍保存。在厭氧發(fā)酵樣品制備過程中，先將風(fēng)干玉米秸稈切成10 mm長的秸稈段，再將秸稈段、豬糞、牛糞、羊糞干燥、粉碎，過40目篩后裝袋備用。依據(jù)秸稈厭氧發(fā)酵過程中堿性預(yù)處理的高效性和生物預(yù)處理的環(huán)境友好性，分別采用地衣芽孢桿菌(生物方法)、NaOH溶液(堿性試劑)、豬糞沼液(富含微生物的弱堿性試劑)和NaOH沼液溶液(富含微生物的強堿性試劑)對玉米秸稈進行預(yù)處理實驗。生物方法預(yù)處理秸稈樣品取自本課題組進行的地衣芽孢桿菌秸稈降解實驗[23]。按最優(yōu)培養(yǎng)條件活化培養(yǎng)地衣芽孢桿菌后，將其液體菌種接種于秸稈粉末固體培養(yǎng)基上，進行為期10 d的降解實驗，每2 d采樣1次，共計采樣5個。其它方法預(yù)處理實驗過程中，將10 mm長的秸稈段浸泡于處理液中3 s后，撈出擠壓排水并裝入自封袋進行密封處理。定時采樣并用蒸餾水充分洗滌樣品5次后，將樣品干燥、粉碎過40目篩后裝袋密封保存，制備預(yù)處理秸稈樣品45個。NaOH、沼液預(yù)處理秸稈實驗方案如表1所示。

將粉碎玉米秸稈按比例與粉碎后的豬糞、牛糞、羊糞粉末進行混合，制備秸稈和糞便混合發(fā)酵原料樣品36個?；旌媳壤?質(zhì)量比)為9∶1、8∶2、7∶3、6∶4、5∶5、4∶6、3∶7、2∶8、1∶9和3個隨機比例。連同玉米秸稈、豬糞、羊糞、牛糞樣品4個，共計采集與制備樣品90個。

1.2 碳氮比的測定

樣品碳、氮含量的測定按照干燒法的原理[24]，采用EURO EA3000型元素分析儀測定，測試模式為碳氫氮模式，測試溫度為980℃，運行時間320 s，樣品杯為5 mm×9 mm錫囊，裝樣質(zhì)量1～3 mg，氦氣為載氣，反應(yīng)管型號為E13041。標(biāo)樣為琥乙紅霉素標(biāo)準(zhǔn)品(C43H75NO16)，其碳、氫、氮、氧質(zhì)量分?jǐn)?shù)分別為59.911%、8.769%、1.625%、29.695%，與玉米秸稈碳氮比接近，適于測試玉米秸稈碳氫氮含量。每個樣品測試3次，取3次的平均值作為待測樣品的碳、氮含量值，然后通過計算可得樣品的碳氮比。

1.3 光譜數(shù)據(jù)采集

對發(fā)酵原料樣品使用Bruker TANGO型近紅外光譜儀進行積分球漫反射光譜掃描，光譜采集范圍3 946～11 542 cm-1(波長866～2 534 nm)，分辨率為8 cm-1，樣品掃描32次，裝樣方式為50 mm樣品杯旋轉(zhuǎn)臺掃描，裝樣質(zhì)量約7 g，將所測吸光度進行保存。在光譜掃描時，使用鍍金樣品杯上蓋壓實樣品粉末，并采用3次掃描的平均值作為樣品的原始光譜。每個樣品原始光譜的波長點數(shù)為1 845個。

1.4 波長優(yōu)選方法

1.4.1GSA算法

GSA算法[25]融合了GA[26]的高效遺傳操作和SA[27]的退溫策略，采用二進制編碼方案，以偏最小二乘(Partial least squares，PLS)回歸模型的交叉驗證均方根誤差(Root mean squared error of cross-validation，RMSECV)為目標(biāo)函數(shù)，通過結(jié)合溫度參數(shù)設(shè)計適應(yīng)度函數(shù)，基于Metropolis判別準(zhǔn)則實現(xiàn)擾動解的選擇復(fù)制，有效解決了GA算法存在的早熟收斂和進化后期搜索效率低的兩點不足，又克服了SA算法進化速度慢的問題。

1.4.2GSA-iPLS算法

針對GSA以全譜波長點個數(shù)為碼長進行二進制編碼時，碼長過長容易導(dǎo)致解空間發(fā)散的問題，將GSA與區(qū)間偏最小二乘法(Interval PLS，iPLS)[28]相結(jié)合構(gòu)建GSA-iPLS算法，通過選取多個有效的光譜子區(qū)間參與建模，能夠有效提高模型的性能。

GSA-iPLS基于iPLS的思想，將NIRS劃分為n個等寬子區(qū)間，然后使用GSA優(yōu)選出有效的特征譜區(qū)建模，以提高模型精度。GSA-iPLS采用二進制編碼方式，以子區(qū)間個數(shù)為碼長，進行GSA的種群初始化?！?”和“0”分別表示對應(yīng)子區(qū)間所包含波長點對應(yīng)的數(shù)據(jù)“是”、“否”選中參與運算。根據(jù)種群初始化結(jié)果計算各染色體的目標(biāo)函數(shù)值，確定初始溫度和降溫操作，并計算各染色體的適應(yīng)度函數(shù)值。然后依據(jù)適應(yīng)度函數(shù)值對種群中的染色體依次執(zhí)行帶最優(yōu)保留策略的賭輪選擇、離散重組交叉、離散變異和Metropolis選擇復(fù)制操作，完成一輪次的GSA種群進化過程。經(jīng)過多個輪次的種群進化，滿足設(shè)定的算法終止條件后，即完成NIRS特征譜區(qū)優(yōu)選。按如上方法，執(zhí)行多次特征譜區(qū)優(yōu)選算法，求出不同子區(qū)間個數(shù)下碳氮比對應(yīng)的多個備選特征子區(qū)間組合，通過綜合評測模型性能后，基于RMSECV確定發(fā)酵原料碳氮比對應(yīng)的最佳子區(qū)間個數(shù)和最佳特征譜區(qū)。

1.4.3DGSA-PLS算法

GSA-iPLS優(yōu)選的光譜子區(qū)間內(nèi)部仍然可能存在不相關(guān)的波長點和波長點之間的共線性問題。DGSA-PLS算法利用GSA對GSA-iPLS優(yōu)選的光譜子區(qū)間進行二次優(yōu)化，能夠有效去除譜區(qū)內(nèi)的不相關(guān)波長點，解決波長點之間的共線性問題。

DGSA-PLS以GSA-iPLS優(yōu)選后特征譜區(qū)包含的特征波長點數(shù)為碼長，進行二進制編碼和種群初始化?！?”和“0”分別表示該波長點對應(yīng)的數(shù)據(jù)“是”、“否”選中參與運算。在確定初始溫度、降溫操作，計算適應(yīng)度函數(shù)值后，執(zhí)行多個輪次的GSA選擇、交叉、變異和Metropolis選擇復(fù)制進化操作，完成NIRS特征波長點的優(yōu)選。針對GSA優(yōu)化結(jié)果的隨機性問題，多次執(zhí)行特征波長點優(yōu)選算法，并選擇多次重復(fù)選中的波長點為特征波長變量建立PLS回歸模型，能夠得到較高的回歸模型性能。

1.5 回歸模型建立及評價

本文算法(包括光譜預(yù)處理、樣品集劃分、GSA-iPLS算法、DGSA-PLS算法及回歸模型構(gòu)建等)全部在Matlab R2012b軟件平臺中實現(xiàn)。

2 結(jié)果與分析

2.1 采集數(shù)據(jù)分析

圖1 樣品光譜數(shù)據(jù)Fig.1 Spectroscopic data of samples

對90個樣品的原始光譜經(jīng)一階導(dǎo)數(shù)預(yù)處理后，使用KS法按2∶1的比例進行樣本劃分，得到校正集樣本60個、驗證集樣本30個，對應(yīng)的碳氮比如表2所示。

表2 厭氧發(fā)酵原料碳氮比Tab.2 C/N ratio of anaerobic fermentation materials

對預(yù)處理后的NIRS進行主成分(Principal components, PCs)分析，第1、第2和第3主成分的貢獻率分別為71.136%、8.596%和5.327%，前3個PCs的累積貢獻率達85.059%。校正集和驗證集的三維主成分空間分布情況如圖2所示。在樣本主成分空間分布圖中，左側(cè)為秸稈糞便混合物樣本對應(yīng)數(shù)據(jù)點，右側(cè)為預(yù)處理秸稈樣本對應(yīng)數(shù)據(jù)點，產(chǎn)生如此清晰分類的結(jié)果與樣品性狀差異、原始光譜數(shù)據(jù)分布吻合。

圖2 樣本主成分空間分布Fig.2 Distribution of samples in PCs space

由表2和圖2可知，校正集樣本碳氮比基本涵蓋了驗證集，且校正集和驗證集樣本在主成分空間上分布比較均勻，可以使用該樣本劃分方法進行NIRS分析。

2.2 特征波長優(yōu)選

2.2.1GSA-iPLS特征譜區(qū)優(yōu)選

GSA-iPLS先按照iPLS將全譜劃分成多個均勻的子區(qū)間，再以子區(qū)間個數(shù)為碼長、以RMSECV為目標(biāo)函數(shù)運行GSA算法，優(yōu)選特定子區(qū)間數(shù)下的特征譜區(qū)。為考察分割波長點個數(shù)對波長選擇及模型預(yù)測性能的影響，分別按約30、40、50、60、80、100、120個波長點劃分子區(qū)間，依次將預(yù)處理后的一階導(dǎo)數(shù)光譜劃分為61、46、37、31、23、18、15個子區(qū)間，依據(jù)RMSECV優(yōu)選有效的子區(qū)間組合作為GSA-iPLS優(yōu)選的特征譜區(qū)。為解決GSA優(yōu)選結(jié)果的隨機性問題，在每個子區(qū)間劃分個數(shù)下，執(zhí)行10次GSA-iPLS算法，并選定回歸模型性能最佳的子區(qū)間組合作為該子區(qū)間數(shù)下的碳氮比特征譜區(qū)。在進行GSA-iPLS特征譜區(qū)優(yōu)選時，種群規(guī)模設(shè)為100，初溫確定系數(shù)取200，降溫系數(shù)取0.950，進化代數(shù)取200，交叉概率取0.700，變異概率取0.010，鄰域解擾動位數(shù)取碼長的1/10向上取整。不同子區(qū)間數(shù)下優(yōu)選的碳氮比特征譜區(qū)信息如表3所示。

表3 GSA-iPLS優(yōu)選結(jié)果Tab.3 Results optimized by GSA-iPLS

由表3可知，采用子區(qū)間劃分個數(shù)為23，優(yōu)選譜區(qū)的選中子區(qū)間數(shù)為8，波長點數(shù)為641時，回歸模型的性能最佳。GSA-iPLS優(yōu)選譜區(qū)如圖3陰影部分所示。

圖3 GSA-iPLS優(yōu)選譜區(qū)Fig.3 Spectral intervals selected by GSA-iPLS

基于各含氫基團在近紅外譜區(qū)中的分布特性可知，在選中的8個子區(qū)間中，3 950～4 935 cm-1(波長2 026～2 534 nm)對應(yīng)著C—C、—CH、—CH2、—CH3和—NH2基團組合頻，7 242～7 567 cm-1(波長1 322～1 381 nm)、7 901～8 226 cm-1(波長1 215～1 266 nm)和8 560～8 885 cm-1(波長1 125～1 168 nm)對應(yīng)著—CH、—CH2和—CH3基團的二級倍頻，9 219～9 544 cm-1(波長1 048～1 084 nm)對應(yīng)著—CH、—NH2基團的三級倍頻。

當(dāng)特征波長點在整個譜區(qū)中分布比較集中時，GSA-iPLS譜區(qū)優(yōu)選算法的性能優(yōu)越，去除冗余波長點的效果較好。當(dāng)特征波長點的分布比較分散時，GSA-iPLS以子區(qū)間為單位進行特征譜區(qū)優(yōu)選，在去除冗余波長點時會連帶去除部分有效波長點，進而影響回歸模型的性能。此時，需要增大子區(qū)間個數(shù)，減小子區(qū)間內(nèi)波長點的數(shù)量，防止GSA-iPLS算法去除過多的有效波長點。但子區(qū)間數(shù)過多時，編碼碼長過長，影響GSA算法搜索效率的同時還可能導(dǎo)致解空間的發(fā)散問題。因此，在進行問題求解時，需要結(jié)合實際情況，設(shè)置合理的算法參數(shù)，實現(xiàn)算法運行效率和求解精度的統(tǒng)一。

2.2.2DGSA-PLS特征波長優(yōu)選

DGSA-PLS在進行特征波長點優(yōu)選時，以GSA-iPLS優(yōu)選的特征譜區(qū)波長點數(shù)為碼長，隨機生成160個碼長為641的染色體構(gòu)建初始種群，鄰域解擾動位數(shù)取20，其它算法初始參數(shù)與GSA-iPLS一致。為消除GSA算法的隨機性，執(zhí)行算法50次對碳氮比特征波長點進行優(yōu)選。多次執(zhí)行時，每次都選中的波長點代表了染色體的優(yōu)良基因，以這些特征波長點作為特征波長變量建立回歸模型時，可以有效消除GSA算法的隨機性，且能夠得到較高的回歸模型性能。DGSA-PLS波長優(yōu)選結(jié)果與預(yù)處理后光譜的平均值對比如圖4所示。

圖4 DGSA-PLS優(yōu)選波長變量Fig.4 Wavelength variables selected by DGSA-PLS

圖4中，重復(fù)選中次數(shù)為1時，選中628個波長點；重復(fù)選中次數(shù)為50時，選中19個波長點。測試發(fā)現(xiàn)，校正集RMSECV和驗證集RMSEP都隨選中波長點數(shù)的增加呈先減小后增大的趨勢，但兩者的趨勢存在較大差別。為了分析特征波長變量數(shù)目與模型性能的關(guān)系，繪制RMSECV、RMSEP與選中波長點數(shù)的關(guān)系圖，如圖5所示。

圖5 RMSE與選中波長點數(shù)間的關(guān)系Fig.5 Relationship between RMSE and number of variables

由圖5可知，RMSECV最小值早于RMSEP出現(xiàn)，當(dāng)選中波長點數(shù)為189、重復(fù)選中次數(shù)為27次時，RMSECV最小。RMSEP最小時，對應(yīng)的波長點數(shù)為628，重復(fù)選中次數(shù)為1次。RMSECV與RMSEP隨選中波長點數(shù)變化趨勢差異較大的主要原因在于GSA以校正集的RMSECV為依據(jù)進行特征波長優(yōu)選，驗證集性能拐點出現(xiàn)時表明校正集發(fā)生了過擬合。因此，選取RMSEP最小時對應(yīng)的628個波長點作為DGSA-PLS優(yōu)選的碳氮比特征波長變量。

2.3 回歸模型評價與分析

為評測本文構(gòu)建的兩種波長優(yōu)選算法的建模性能，以GSA-iPLS和DGSA-PLS優(yōu)選后的特征波長變量作為PLS回歸模型的輸入，建立厭氧發(fā)酵原料碳氮比定量回歸模型，并與全譜(Full-PLS)、協(xié)同區(qū)間偏最小二乘(Synergy iPLS，SiPLS)[30]和反向區(qū)間偏最小二乘(Backward iPLS，BiPLS)[31]優(yōu)選特征波長的建模性能進行對比，結(jié)果如表4所示。

表4 不同回歸模型評價指標(biāo)Tab.4 Evaluation indicators of different regression models

由表4可知，SiPLS和BiPLS作為兩種最典型iPLS算法，SiPLS的建模性能弱于全譜建模，而BiPLS的建模性能優(yōu)于全譜建模。主要原因在于：SiPLS選取2～4個固定個數(shù)的子區(qū)間作為備選譜區(qū)，再通過比較RMSECV確定最佳譜區(qū)；而BiPLS通過剔除相關(guān)性較差的子區(qū)間，搜索RMSECV最小的子區(qū)間組合作為特征譜區(qū)；BiPLS比SiPLS更適于特征波長點分布比較分散問題的求解。而碳氮比對應(yīng)著譜區(qū)中所有含碳和含氮基團的吸收峰，這些吸收峰在整個譜區(qū)中分布較廣，適合采用BiPLS進行特征譜區(qū)優(yōu)選。GSA-iPLS作為一種新型近紅外光譜特征譜區(qū)優(yōu)選算法，具有良好的隨機搜索能力。在使用GSA-iPLS算法進行特征譜區(qū)優(yōu)選時，與BiPLS相比擴展了搜索結(jié)果的隨機性。通過多次搜索并選取建模性能最佳的搜索結(jié)果作為GSA-iPLS優(yōu)選特征譜區(qū)，該方式能夠有效提高算法的特征波長優(yōu)選性能。

圖6 碳氮比實測值與預(yù)測值分布Fig.6 Distribution of measured and predicted values for C/N ratio

3 結(jié)束語

探討了采用NIRS技術(shù)結(jié)合化學(xué)計量學(xué)方法進行厭氧發(fā)酵原料碳氮比快速檢測的可行性。為提高NIRS回歸模型的檢測精度和效率，基于GSA算法構(gòu)建了GSA-iPLS和DGSA-PLS兩種算法進行碳氮比特征波長的優(yōu)選。GSA-iPLS將光譜數(shù)據(jù)劃分成多個子區(qū)間后，以子區(qū)間個數(shù)為碼長，搜索有效的特征波長子區(qū)間組合作為特征譜區(qū)，有效減少了建模變量個數(shù)，提高了碳氮比檢測模型的精度和效率。DGSA-PLS在GSA-iPLS優(yōu)選譜區(qū)的基礎(chǔ)上，以波長點個數(shù)為碼長進行特征波長變量優(yōu)選，有效去除不相關(guān)的冗余波長點，得到628個特征波長點，建立的碳氮比檢測模型RMSEP為7.178，RPD為3.805。與全譜建模相比，基于DGSA-PLS建立的回歸模型有效波長點個數(shù)減少了65.94%，RMSEP減小了15.87%，有效地提高了模型的檢測精度。