采樣誤差分布分析法建立氯丙醇含量的穩健近紅外光譜分析模型

陳萬超 陶鑫 范長春 張飛宇 杜一平

摘?要?利用近紅外光譜技術結合采樣誤差分布分析（SEPA）方法建立了二氯丙醇產品生產過程的氯化液雜質3-氯-1，2-丙二醇濃度的分析模型。對樣本數據進行1000次隨機劃分，建立1000個子模型，獲得多個潛變量數下的交互檢驗誤差，進行統計分析。繪制了誤差分布圖，計算其中位數、標準偏差、偏斜度和分布峰度等統計指標，通過這些指標的綜合分析對近紅外光譜分析模型進行條件優化、建模和模型評價等。4種光譜處理方法顯示出比較理想的模型性能，作為候選與不同波長區域的選擇相結合，繼續運用SEPA運算，進一步優化模型。最終優化的建模條件為：一階導數結合標準正態變換; 6931～6017 vcm1波數區間;使用5個偏最小二乘潛變量。校正、交互檢驗和獨立驗證誤差分別為0.881%、1.282%和1.167%。所選擇的波長具有可解釋性，模型的各項統計參數合理、可信。研究結果表明，SEPA能全面、合理地考察多項統計指標，可以建立實用、穩健的近紅外光譜分析模型。

關鍵詞?采樣誤差分布分析; 近紅外光譜; 建模; 穩健

1?引 言

近紅外光譜分析技術具有簡單、快速、無損、準確等優點，近年來應用廣泛。應用近紅外光譜法時，通常需要借助化學計量學建立光譜多元分析模型，實現定性和定量分析，在實際應用時，建模常是難點，也容易出現錯誤。近紅外光譜分析建模一般包括條件優化（光譜預處理方法選擇和波長選擇等）、建模和模型評價等部分，這些工作均離不開數據集合的劃分，如將數據分為校正集、交互檢驗集和獨立驗證集等，不同的集合劃分常獲得不同的建模結果。近年來，集合劃分方法及其相關研究備受關注[1～5]。劃分方法大體可分為人為劃分和隨機劃分兩大類，前者根據指標性質數值或光譜之間距離[6]，或者從指標和光譜相關關系考慮樣品之間的差異[2]劃分樣本集;后者通過多次隨機采樣，達到選擇具有代表性樣本數據集的目的[1]。對于復雜的近紅外光譜數據，經常包含特殊數據，如指標性質數值相差不大，或光譜卻有較大差異。這種數據采用人為劃分方法難以得到理想的結果，而隨機方法，只要樣本采樣次數足夠大，各種情況的數據都能被考慮到，因此隨機劃分方法極具競爭力。Xu等[1]提出的蒙特卡羅交互檢驗法（MCCV）很好地利用了這種隨機采樣的特性，多次隨機劃分數據集，建立多個子模型，對這些子模型的誤差取平均，用于交互檢驗。Li等[7]提出的模型集群分析方法（MPA）采用類似的策略建立多個子模型，并考慮了集群的子模型多種統計性質用于建模。本研究組在MPA算法的基礎上，近期提出了采樣誤差分布分析（SEPA）算法[8]，利用多種統計指標綜合分析隨機采樣誤差，如誤差分布的中位數和標準偏差、分布偏斜度和分布峰度等，用于近紅外光譜分析中的異常點篩查、光譜預處理方法和波長選擇、交互檢驗、模型評價，還用于模型轉移[9]和波長選擇[10]等。多次隨機采樣結合多種統計指標的綜合考慮，能更全面地考察數據的各種狀況，在此基礎之上的條件優化和模型建立及評價可以更加合理，所建立的模型更加穩健。

二氯丙醇化工生產中的關鍵副產物一氯丙二醇對產品質量控制至關重要。由于主產物和被測物質結構非常相似，同時樣品中還含有含量較小的其它化合物，所測的近紅外光譜中難以獲得被測組分清晰的特征譜峰，光譜干擾嚴重。本研究以此副產物為研究目標，采用近紅外光譜進行分析研究， 利用SEPA方法對近紅外光譜數據進行了細致研究，建立了合理、穩健的分析模型。

2?實驗方法與理論

2.1?儀器與數據

樣本取自江蘇省某化工企業二氯丙醇生產工段的氯化液，樣品為液體。用Antaris II近紅外光譜儀（美國賽默飛公司）測定樣品的透射光譜，比色皿光程為8 mm，波數范圍3799～11999 cm1，分辨率4 cm1，去除吸光度飽和的數據，實際使用的光譜范圍為5361～11999 cm1。掃描次數為16，取平均光譜。3-氯-1，2-丙二醇的濃度采用GC7820氣相色譜儀（美國安捷倫公司）進行測定，方法參考國標法GB/T 13097-2007[11]。

2.2?采樣誤差分布分析（SEPA）

利用SEPA[8]對近紅外光譜數據進行光譜預處理、波長選擇、建模和模型評價，均以交互檢驗的形式進行工作。因此，在此以交互檢驗過程為例說明SEPA。用X表示光譜，維度nm，y表示指標濃度，維度n1，n和m分別表示樣本數和波長數，隨機采樣次數為N，最高潛變量數為Ncv，建模方法為偏最小二乘法（PLS）。首先通過隨機采樣將樣本集合劃分為校正集和交互驗證集，采樣N次，獲得N個不同的樣本集合劃分;然后分別對每個樣本集合用校正集建立校正子模型，并針對各個潛變量數（nLVs）計算交互驗證預測均方根誤差（RMSECV）;對所有得到的N個子模型的RMSECV 誤差數據（共計NNcv）進行統計分析，在每個潛變量數下，繪制誤差分布圖，并計算其誤差中位數、標準偏差（STD）、分布偏斜度和分布峰度;最終綜合考慮上述統計結果， 進行條件選擇和模型評價。

采用中位數而非平均值估計RMSECV，是因為前者更加穩健，標準偏差STD能反映誤差分布寬度，其值越小， 說明采樣偏差越小，對建模越有利;而分布偏斜度和分布峰度則是判斷誤差分布的合理性，前者反映分布的對稱性，理想值為0，后者表示分布的峰形態（簡單理解就是峰太瘦或太胖的程度），其理想值為3。由于采用隨機選擇樣本集劃分，并多指標綜合考慮誤差分布，因此所選擇的參數和建立的模型更加合理和穩健。

2.3?數據分析

樣本總數為83，被測指標為3-氯-1，2-丙二醇的含量，光譜有1722個波長點。采用Kennard -Stone（K-S）方法[7]將83個樣本劃分為校正集和驗證集，前者用于條件優化和建模，后者用來對模型進行獨立的檢驗，以此設置校正集68個樣本，驗證集15個樣本。條件優化過程均采用交互檢驗的方法，在建模集68個樣本中，隨機選取20個樣本組成交互檢驗集，其它樣本進行建模（稱為校正集）。在進行SEPA運算時，隨機采樣次數為1000。表1列出了被測指標3-氯-1，2-丙二醇含量的相關數據。所有計算均在Matlab（R2010a）環境下自編程序完成。

3?結果與討論

3.1?近紅外光譜

在所測得的近紅外光譜中，低波數的吸光度較高，但噪聲也高，而高波數位置，除了在8600 cm1附近有一個小峰外，基本無吸收。明顯的吸收峰包括8600、6900、6300、6000和5800 cm1（這些譜峰在建模中的作用參見3.5節）。

3.2?采用原始數據初步建模

SEPA的核心是多次隨機采集校正集和交互檢驗集，建立多個子模型，通過對子模型的綜合分析進行條件優化和建模。為了說明SEPA的核心工作過程和數據分析方法，用建模集的68個樣本的原始光譜數據，采用SEPA方法初步建模。交互檢驗集樣本數為20，校正集樣本數為48，隨機采集這些樣本次數800，即可建立800個子模型，每個子模型都產生各個潛變量數（nLVs=1～15）的一組交互檢驗誤差，對每個潛變量數的800個誤差進行統計分析，繪制誤差分布圖，在每個潛變量數下，取其中位數作為交互檢驗誤差，對nLVs作圖。

在潛變量數為8時交互檢驗誤差最小，其標準偏差也最小。此時的誤差分布均勻，相對比潛變量數少于8時，更能體現出正態分布的輪廓，作為比較給出4個潛變量的誤差分布，其分布不如8個潛變量的誤差分布。但8個潛變量的誤差分布左右略不對稱，計算發現其偏斜度為1.479，峰度為12.32，距的理想值0和3相差較遠。這些不足可以通過光譜預處理和波長選擇進一步得到改善（參見3.4和3.5節）。

3.3?隨機采樣次數的選擇

數據集劃分通過隨機采樣完成，隨機采樣次數分別選擇600、800、1000、1200，按3.2節的方法進行交互檢驗，重復計算5次，考察結果的穩定性。結果表明，采樣次數為600時，潛變量數有時選擇8，有時選擇9;采樣次數為800時，潛變量數多選擇8，偶爾選擇9;當采樣次數為1000或更高時，潛變量數均選擇8。而重復計算誤差的相對標準偏差RSD（%）分別為0.7578、0.8845、0.1361和0.2366。因此，選擇隨機采樣次數為1000，計算時間約20 s。

3.4?光譜預處理方法的選擇

光譜預處理可以改善模型精度，而不同的方法常會導致預測誤差和最佳潛變量數均發生明顯的變化。常見的光譜預處理方法[12]包括Savitzky-Golay平滑（用S表示）、多元散射校正（MSC）、標準正態變換（SNV）、一階導數（1D）、二階導數（2D）以及它們的組合。考察了25種有效的單一方法及其組合，分別進行SEPA計算。不同光譜預處理方法的交互檢驗誤差隨潛變量數變化圖，因為有些方法得到的變化圖非常相似，圖中只繪制了其中的一種。有4種方法（SNV、1D、1D+S和1D+SNV）的誤差較低，而且潛變量數也較低。它們的交互檢驗誤差圖及標準偏差（誤差棒），SNV的預測誤差最低，為1.391%，但潛變量數較高（nLVs=7），另外3種方法的預測誤差稍高，但可以選擇更低的潛變量數，如果選擇4，預測誤差分別為1.568%、1.552%和1.495%。 本研究將這4種方法作為候選方法，通過波長選擇進一步優化模型。

3.5?波長選擇

近年來，發展了很多波長選擇方法[13～20]，如移動窗口偏最小二乘法[13]、無信息變量剔除法[14]、競爭性自適應加權抽樣法[15]和基于集群分析變量提取法[16]等，各有優缺點。為了簡化，本研究采用最簡單的相關系數法選擇候選波長區域，即計算濃度與各波數處吸光度的相關系數R。各個波數吸光度與濃度的相關系數圖，單一波數的相關系數R都不高，按R>0.4選擇3個波數區間： R1=6931～6017 cm1、R2=5959～5897 cm1、R3=5758～5523 cm1。光譜中8600 cm1峰因為強度太低，在波長選擇中沒有體現。5800 cm1位置因為干擾物也同樣具有強吸收，在波長選擇中也沒有體現。波長選擇的R2和R3譜帶為5800 cm1兩邊區域，避開了強干擾區域。R1主要反映6900、6300和6000 cm1 3個譜帶。

被測組分3-氯-1，2-丙二醇在4769、4421和3995 cm1具有明顯的吸收峰[21]，但樣品中含有其它含氯丙醇，含量較高的是二氯丙醇，對被測組分構成比較嚴重的光譜干擾。3-氯-1，2-丙二醇在7000～5700 cm1范圍有一個很寬的吸收帶，而主要干擾物二氯丙醇除在7000和5900 cm1處和5700～5800 cm1范圍有明顯吸收峰，其它區域吸收較弱[21]。3個譜帶正是反映了被測組分和干擾物之間的光譜差異。在R1（綠色）、R2（橙色）和R3（藍色）波數范圍，3-氯-1，2-丙二醇的吸收比較強，而二氯丙醇的吸收比較弱，尤其是R1帶。

將3.4節中的4種光譜預處理方法與所選擇的波長區間組合，利用SEPA綜合考察模型性能。表2給出了模型的各項指標，以及在各種優化條件下，建模后采用驗證集計算預測誤差RMSEP。

與原始光譜相比，光譜預處理后模型得到明顯改善，最佳方法為SNV方法，在7個潛變量下，模型的建模集和驗證集誤差RMSEC和RMSEP為0.985%和1.195%。進一步考慮波長選擇，當使用3個波長區段建模時，各種誤差并未降低，原因可能是在R2和R3段光譜干擾依然很嚴重。R1區域的光譜干擾相對比較弱，如果只選擇R1進行建模，模型誤差進一步降低。最佳結果是1D+SNV結合R1，在5個潛變量下，RMSEC和RMSEP降為0.881%和1.167%。1000次計算交互檢驗誤差的分布圖，4種方法的分布都比較均勻，分布的對稱性也得到顯著改善。最佳方法1D+SNV+R1的偏斜度為0.4726，峰度為3.9573，與原始光譜建模的1.4789和12.3246（見3.2節）相比，更接近理想值的0和3。另一外值得指出的是，標準偏差STD在選用不同的處理方法和波長時變化不大，表中STD均在0.186～0.270之間，因此STD指標在條件優化時作用不大。

3.6?最終模型

在最佳條件下，即光譜處理方法為一階導數1D結合SNV，波長區間為R1（6931～6017 cm1），潛變量數為5，用建模集建立模型，對驗證集進行預測。建模集和驗證集樣本的預測結果。各個樣本基本上均勻分布在斜率為1的直線兩側，雖然誤差并不太小，但可能是這套數據本身的特征，所建模型反映的正是數據的本質。 RMSEC和RMSEP為0.881%和1.167%，對應的相關系數分別為0.95和 0.88。

4?結 論

利用采樣誤差分布分析方法對二氯丙醇產品生產過程的氯化液進行近紅外光譜分析，建立了雜質3-氯-1，2-丙二醇濃度的分析模型。SEPA方法通過對樣本數據多次隨機劃分建立多子模型，進而計算中位數、標準偏差、偏斜度和分布峰度等統計指標，通過綜合分析對近紅外光譜分析模型進行條件優化、建模和模型評價等工作。結果表明，在光譜預處理方法選擇和波長選擇中，SEPA能選擇更加合理的結果，波長選擇的結果更具解釋性。用獨立驗證集進行模型評價所建模型的穩健性。本研究為建立實用的近紅外光譜分析模型提供了一種有效的方法。

Reference

1?Xu Q S， Liang Y Z. Chemom. Intell. Lab. Syst.， 2001， 56（1）： 1-11

2?Galvao R， Araujo M， Jose G， Pontes M， Silva E， Saldanha T. Talanta， 2005， 67（4）： 736-740

3?Deng B C， Yun Y H， Liang Y Z， Cao D S， Xu Q S， Yi L Z， Huang X. Anal. Chim. Acta， 2015， 880： 32-41

4?JIN Zhao-Xi， ZHANG Xiu-Juan， LUO Fu-Yi， AN Dong， ZHAO Sheng-Yi， RAN Hang， YAN Yan-Lu. Spectrosc. Spect. Anal.， 2016， 36（12）： 3920-3925

靳召晰， 張秀娟， 羅付義， 安 冬， 趙盛毅， 冉 航， 嚴衍祿. 光譜學與光譜分析， 2016， 36（12）： 3920-3925

5?YUN Yong-Huan， DENG Bai-Chuan， LIANG Yi-Zeng. Chinese J. Anal. Chem.， 2015， 43（11）： 1638-1647

云永歡， 鄧百川， 梁逸曾. 分析化學， 2015， 43（11）： 1638-1647

6?Kennard R W， Stone L A. Technometrics， 1969， 11（1）： 137-148

7?Li H D， Zeng M M， Tan B B， Liang Y Z， Xu Q S， Cao D S. Metabolomics， 2010， 6（3）： 353-361

8?Chen W C， Du Y P， Zhang F Y， Zhang R Q， Ding B Y， Chen Z K， Xiong Q. J. Chemom.， 2017： e2933

9?Zhang F Y， Chen W C， Zhang R Q， Ding B Y， Yao H M， Ge J， Ju L， Yang W Y， Du Y P. Chemom. Intell. Lab. Syst.， 2017， 171： 234-240

10?Zhang R Q， Zhang F Y， Chen W C， Yao H M， Ge J， Wu S C， Wu T， Du Y P. Chemom. Intell. Lab. Syst.， 2018， 175： 47-54

11?GB/T 13097-2007， Determination of Phthalate Esters in Foods. National Standards of the People's Republic of China

工業用環氧氯丙烷. 中華人民共和國國家標準， GB/T 13097-2007

12?Burns D A.， Ciurczak E W. Handbook of Near-Infrared Analysis， 2008： CRC Press

13?Jiang J H， James B R， Siesler Heinz W， Yukihiro O. Anal. Chem.， 2002， 74（14）： 3555-3565

14?Cai W S， Li Y K， Shao X G. Chemom. Intell. Lab. Syst.， 2008， 90（2）： 188-194

15?Zheng K Y， Li Q Q， Wang J J， Geng J P， Cao P， Sui T， Wang X， Du Y P. Chemom. Intell. Lab. Syst.， 2012， 112（6）： 48-54

16?Yun Y H， Wang W T， Deng B C， Lai G B， Liu X B， Ren D B， Liang Y Z， Fan W， Xu Q S. Anal. Chim. Acta， 2015， 862： 14-23

17?Lin Y W， Xiao N， Wang L L， Li C Q， Xu Q S. Chemom. Intell. Lab. Syst.， 2017， 168： 62-71

18?Deng B C， Yun Y H， Liang Y Z， Yi L Z. Analyst， 2014， 139： 4836-4845

19?Deng B C， Yun Y H， Ma P， Lin C C， Ren D B， Liang Y Z. Analyst， 2015， 140（6）： 1876-1885

20?Wang L L， Lin Y W， Wang X F， Xiao N， Xu Y D， Li H D， XU Q S. Chemom. Intell. Lab. Syst.， 2018， 172： 229-240

21?Michael B， Hans P V. FT-NIR Atlas. 1993： VCH publisher. 76， 89