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

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

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

0 引言

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

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

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

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

1 材料與方法

1.1 樣品采集與制備

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

將粉碎玉米秸稈按比例與粉碎后的豬糞、牛糞、羊糞粉末進行混合，制備秸稈和糞便混合發酵原料樣品36個。混合比例(質量比)為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錫囊，裝樣質量1～3 mg，氦氣為載氣，反應管型號為E13041。標樣為琥乙紅霉素標準品(C43H75NO16)，其碳、氫、氮、氧質量分數分別為59.911%、8.769%、1.625%、29.695%，與玉米秸稈碳氮比接近，適于測試玉米秸稈碳氫氮含量。每個樣品測試3次，取3次的平均值作為待測樣品的碳、氮含量值，然后通過計算可得樣品的碳氮比。

1.3 光譜數據采集

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

1.4 波長優選方法

1.4.1GSA算法

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

1.4.2GSA-iPLS算法

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

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

1.4.3DGSA-PLS算法

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

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

1.5 回歸模型建立及評價

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

2 結果與分析

2.1 采集數據分析

圖1 樣品光譜數據Fig.1 Spectroscopic data of samples

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

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

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

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

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

2.2 特征波長優選

2.2.1GSA-iPLS特征譜區優選

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

表3 GSA-iPLS優選結果Tab.3 Results optimized by GSA-iPLS

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

圖3 GSA-iPLS優選譜區Fig.3 Spectral intervals selected by GSA-iPLS

基于各含氫基團在近紅外譜區中的分布特性可知，在選中的8個子區間中，3 950～4 935 cm-1(波長2 026～2 534 nm)對應著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)對應著—CH、—CH2和—CH3基團的二級倍頻，9 219～9 544 cm-1(波長1 048～1 084 nm)對應著—CH、—NH2基團的三級倍頻。

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

2.2.2DGSA-PLS特征波長優選

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

圖4 DGSA-PLS優選波長變量Fig.4 Wavelength variables selected by DGSA-PLS

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

圖5 RMSE與選中波長點數間的關系Fig.5 Relationship between RMSE and number of variables

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

2.3 回歸模型評價與分析

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

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

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

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

3 結束語

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