基于定量結構-性質相關性的烴類物質爆炸下限預測*

潘 勇,蔣軍成,王 睿

(1.中國科學技術大學火災科學國家重點實驗室,安徽 合肥230026;2.南京工業(yè)大學城市建設與安全工程學院,江蘇 南京210009)

1 引 言

爆炸極限是評價可燃氣體或蒸氣爆炸危險性的重要參數之一,在防爆技術中應用廣泛。它是風險評估以及確定可燃物質儲存、運輸、生產、使用安全性的依據之一;在監(jiān)測監(jiān)控技術中,它是一個重要的爆炸指示參量。因此,掌握各類可燃物質的爆炸極限數據具有重要的現實意義。

通過實驗測定是目前獲取爆炸極限數據的有效方法。但實驗方法對設備要求高、工作量大;實驗結果影響因素眾多,差別較大。同時,對于測量上有困難或尚未合成的物質,也無法基于實驗進行測定。因此,有必要開發(fā)簡便可靠的爆炸極限理論預測模型,以彌補實驗方法的不足。已有的爆炸極限理論預測方法主要包括經驗關聯法和基團貢獻法2 大類[1],這些方法具有計算簡單、使用簡便等優(yōu)點。然而,在實際使用過程中也存在明顯的缺陷:經驗關聯法需要使用臨界壓力、等張比容等不常見的物化參數,而這些參數的實驗數據往往本身就比較缺乏;基團貢獻法的應用范圍則受研究體系的影響較大,如果某個基團不在建模時選取的那組基團范圍之內,那么對于含有該基團的化合物就無法應用該模型進行預測。因此,上述方法的應用范圍都具有較大的局限性,限制了在實際工程中的進一步應用。

相關研究[2-3]表明,定量結構-性質相關性研究是一種能夠根據分子結構實現有機物理化性質預測的有效方法,它根據化合物性能與分子結構密切相關的原理,尋求分子結構與物質性質之間的內在定量關系。其基本假設是有機物的性能與分子結構密切相關,分子結構不同,性能就不同。而分子結構可以用反映分子結構特征的各種參數來描述,即有機物的理化性質可用化學結構的函數來表示。根據對分子結構參數和實驗數據的內在定量關系,采用合適的統計建模方法進行關聯,建立分子結構參數與理化性質之間的關系模型。一旦建立了可靠的QSPR 模型,僅需分子的結構信息,就可預測新的或尚未合成的有機物的各種性質。目前,該方法已被廣泛應用于有機物理化性質及生物活性的預測研究[3]。

本文中根據QSPR 研究基本原理,從分子結構角度出發(fā),對354 種烴類物質的爆炸下限與其分子結構間的內在定量關系進行研究,建立根據分子結構預測烴類物質爆炸下限的理論模型。

2 研究方法

2.1 數據來源

實驗數據的準確度直接影響到所建QSPR 模型的預測精度。爆炸下限實驗值與其測定方法和測定裝置有關,由于不同研究者采用的實驗裝置和方法不同,測定的爆炸下限數據之間存在一定的差異。為排除這些因素可能對預測結果帶來的影響,確保實驗樣本的可靠性和標準化,本文中統一選用美國化學工程師協會下屬DIPPR 數據庫中的354 種烴類化合物的爆炸下限數據作為實驗樣本。這些化合物涵蓋廣泛的化學多樣性空間,為建立有效的預測模型奠定了基礎。

2.2 結構參數的計算和刪減

要想尋求化合物性能與其分子結構間的相關性,必須首先對反映化合物分子結構特征和信息的各種結構參數(如拓撲參數、組成參數、電性參數以及幾何參數等)進行計算,以獲取分子的結構信息,實現分子結構的參數化描述。本文中首先將所有樣本化合物結構輸入Hyperchem 7.0 程序中,采用分子力學方法MM+進行初步優(yōu)化,在此基礎上用半經驗方法AM1 進一步進行幾何優(yōu)化。隨后將優(yōu)化完的分子結構輸入Dragon 2.1 程序中,對化合物的各種結構參數進行計算。通過計算,每個化合物均得到1 481 種結構參數。

同時,為避免建模過程中“機會相關”現象的發(fā)生,必須預先對眾多的結構參數進行刪減,以剔除那些不能為模型提供有用信息的參數。首先,對于那些對所有化合物來說數值為常數或者近似常數(變化較小)的結構參數,由于無法對導致性質差別的結構差異性進行有效表征,因此被刪除;同時,對于兩者之間相關系數大于0.96 的任意2 個結構參數,由于存在“共線性”,因此其中之一也被刪除。經過刪減,共剩下445 種結構參數進行下一步的特征變量篩選。

2.3 基于遺傳算法的特征變量選擇

如何從眾多的結構參數中應用特征變量選擇方法篩選出與目標性質最密切相關的結構參數,是QSPR 研究非常關鍵的問題。目前常用的特征變量選擇方法主要包括3 類:(1)基于多元線性回歸(M LR)的變量篩選方法,如向前/后選擇變量法、逐步回歸法等;(2)基于偏最小二乘(PLS)的變量篩選法,包括修正PLS 權重或系數以消除模型中無用變量的方法等;(3)基于搜索算法的變量篩選方法,如模擬退火法、遺傳算法(GA)等搜索算法和M LR、PLS 等多種建模方法相結合的變量篩選方法。其中,第1 類方法適用于變量間不存在多重共線性數據的變量篩選,優(yōu)點是方法簡單直觀,缺點是它們不能遍歷所有的變量組合,也就不能保證尋找到變量空間里的最優(yōu)解。第2 類方法與第1 類方法類似,僅搜索變量空間的某些范圍,而不具備全局搜索能力,因此他們得到的常常是變量空間的局部最優(yōu)解,而非全局最優(yōu)解。

遺傳算法是J.H.Holalnd[3]于1975 年最早提出的一種模擬生物在自然環(huán)境下的遺傳和進化過程而形成的一種自適應全局優(yōu)化概率搜索方法。它對需要解決問題的參量進行編碼運算,并沿多種路線進行平行搜索,一般不會陷入局部最優(yōu)的陷阱中,能夠在多個局部較優(yōu)中找到全局最優(yōu)解。由于GA 具有相當強的搜索能力,當它與MLR、PLS 等建模方法相結合后,在一定條件下能夠在有限的時間內搜尋到變量空間的最佳模型。因此,近年來GA 在QSPR 研究中得到了較為廣泛的應用[5]。

本研究將GA 與PLS 相結合(GA-PLS),對剩余的445 種結構參數進行優(yōu)化篩選,找出與爆炸下限最密切相關的結構參數作為分子描述符,表征化合物的結構特征。GA-PLS 程序通過Matlab 軟件實現,選用“留1/10 法”交互驗證的均方根誤差(RMS)作為適應度函數。“留1/10 法”交互驗證是指從訓練集中每次篩除訓練樣本數的1/10 個化合物,用其余的化合物建模,來預測所篩除化合物的性質,這樣得到一個交互驗證的RMS 來評價模型性能的好壞。GA-PLS 所選用的其他相關參數為:種群規(guī)模30;初始染色體串長度5;最大染色體串長度30;突變概率1%;最大主成分數由對包含所有自變量的模型進行交互驗證確定(≤15);最大進化代數100。

2.4 預測模型的建立

最優(yōu)分子描述符確定以后,以它為輸入變量,分別采用多元線性回歸和支持向量機2 種方法,對烴類物質爆炸下限與分子描述符間的內在定量關系進行關聯,建立相應的預測模型。

多元線性回歸分析是一種較為常用的多元統計方法,它以線性最小二乘為基礎,通過對輸入變量和輸出變量間的內在定量關系進行統計學習,建立相應的線性模型。本文中MLR 分析采用統計分析軟件SPSS 11.0 實現,在95%置信區(qū)間,采用全回歸算法對烴類物質爆炸下限數據與其分子描述符數據之間的線性關系進行擬合。

支持向量機(SVM)算法是V.N.Vapnik 等[6]在統計學習理論的基礎上,于1995 年提出的一種新型機器學習方法。其基本思想是通過事先確定的非線性映射將輸入向量x 映射到一個高維特征空間中,然后在此高維空間中構建最優(yōu)超平面,將問題轉化為2 次規(guī)劃。在2 次規(guī)劃中,目標函數只涉及內積運算。如果采用核函數就可以避免在高維空間進行復雜運算,而通過原空間的函數來實現內積運算。因此,選擇適當的內積函數就可以實現某一非線性變換后的線性計算,而計算復雜度卻沒有增加。由于SVM 具有嚴格的理論基礎,能較好地解決小樣本、非線性、高維數和局部最小等實際問題,因此在QSPR 研究領域也開始得到應用[7]。本文中SVM 算法采用Libsvm 軟件運行。為獲得最佳的泛化能力,在建模過程中需要調節(jié)相應的參數組合,即選擇合適的核函數、確定核函數的參數、懲罰系數C 以及ε-不敏感損失函數中ε的大小。目前實際應用中最多的是徑向基(RBF)形式的核函數,它具有較高的學習效率和學習速率。對于RBF核函數,最重要的參數是核函數的寬度γ,它與C 及ε一起,決定了SVM 的泛化能力及預測性能。由于這幾個參數之間有很大的相關性,因此采用格點搜索方法尋找最優(yōu)的參數組合。采用“留一法”交互驗證方式進行模擬,最優(yōu)模型所對應的那組參數即被確定為最優(yōu)參數。

2.5 模型驗證

模型驗證對于QSPR 研究至關重要,在QSPR 模型建立以后,必須對模型的穩(wěn)定性、預測能力及泛化能力進行檢驗。常用的驗證方法有內部驗證與外部驗證2 種。“留一法”交互驗證是常用的內部驗證方法,其結果以“留一法”交互驗證的復相關系數表示。內部交互驗證主要是為了檢驗模型的穩(wěn)定性以及其內部預測能力,而相關研究表明[8],內部驗證結果的好壞并不能說明其外部預測能力的大小,對模型預測能力的評價必須通過對未參與訓練的物質進行預測,即外部驗證的方式來進行。通過對預測集化合物進行外部驗證,既能夠體現模型對未參與訓練的樣本的真實預測能力,又能夠反映模型的泛化能力。模型的外部預測能力常用交互驗證的外部復相關系數進行衡量。因此,將樣本集隨機劃分為訓練集(樣本總數80%)和預測集(樣本總數20%)2 部分,其中,訓練集主要用于變量篩選和構建預測模型,預測集則主要用于對所建模型進行外部驗證。

此外,均方根誤差RMS 和平均絕對誤差AAE 作為重要的統計學參數,也被用來衡量所建模型的估計能力及預測性能。

3 結果與討論

3.1 GA-PLS結果

運用GA-PLS 方法進行特征變量篩選,通過比較各模型預測RMS 的大小,確定本研究中與烴類物質爆炸下限最密切相關的4 個結構參數列于表1。表中,x1、x 2 和x3均屬于拓撲描述符,主要描述分子中原子的連接信息。其中,x1與分子的相對頂點距離復雜度有關;x2則主要描述分子的形狀,與分子中元素的差異度密切相關;x3則主要反映整個分子的拓撲結構。x4屬于RDF 描述符,主要描述整個分子中的原子間距離信息,除此之外,它還與分子中鍵的距離、環(huán)的類型、平面和非平面體系及原子類型等信息有關。

表1 GA-PLS篩選出的結構參數Table 1 The descriptors selected by GA-PLS

從以上討論可以看出,4 個結構參數主要與分子的一些簡單結構特征有關,如分子的大小和形狀、分子中原子的分布等。因此,可以將影響烴類物質爆炸下限的主要結構特征歸納為分子的大小、形狀以及分子中原子的分布信息等。

對于任意已知分子結構的化合物,上述4 個結構參數均可以通過Dragon 軟件進行計算,各參數的具體計算公式及過程可參見文獻[9]。

3.2 MLR 模型結果

以GA-PLS 提取的4 個結構參數作為輸入變量,針對訓練集樣本,運用M LR 方法建立爆炸下限線性預測模型

式中:φ為爆炸下限,n 為樣本數,R2為復相關系數,F 為F 檢驗值,S 為標準誤差,P 為方程的顯著性概率。由式(1)可知,模型具有較高的相關系數和較低的標準誤差,說明模型是可靠的;顯著性概率遠小于0.05,表明回歸方程具有統計學意義。

隨后,應用式(1)分別對訓練集樣本和預測集樣本進行預測。對預測集中70 個樣本的預測值見表2,對所有樣本所得的爆炸下限預測值與目標值的比較見圖1,模型的主要性能參數見表3。

表2 預測集樣本爆炸下限預測值Table 2 The predicted values for the test set

表3 MLR 和SVM 模型的主要性能參數比較Table 3 Performance parameters obtained by SVM and MLR models

3.3 SVM 模型結果

為了與MLR 模型具有可比性,SVM 模擬中選取的訓練集和預測集與M LR 相同,采用同樣的4 個結構參數作為模型輸入變量。通過格點搜索方法確定SVM 模型的最優(yōu)參數如下:C=512,ε=0.031 25,γ=0.25,相應的支持向量數為129。應用所建模型分別對訓練集樣本和預測集樣本進行預測。對預測集中70 個樣本的預測值見表2,對所有樣本所得的爆炸下限預測值與目標值的比較見圖2,模型的主要性能參數見表3。

圖1 MLR 模型爆炸下限預測值與目標值的比較Fig.1 Comparison between the predicted and observed LFL by MLR

圖2 SVM 模型爆炸下限預測值與目標值的比較Fig.2 Comparison betw een the predicted and observed LFL by SVM

3.4 模型穩(wěn)定性分析

為進一步排除所建模型可能存在的“機會相關”現象,采用“y-scrambling”方法分別對MLR 和SVM 模型進行驗證,評估2 個模型對“機會相關”現象的依賴程度。首先將訓練集樣本中的自變量x 保持不變,將對應的應變量y 順序隨機打亂,使應變量和自變量不再對應,以消除兩者之間存在的內在定量關系;隨后,針對上述改變序列的樣本集建立新的QSPR 模型,并對其相關性能參數如RMS 等進行計算;將上述過程重復50 ～100 次,獲得某一相對“最優(yōu)”的預測模型,并將其與基于原始樣本數據建立的實際預測模型進行比較。若實際模型的性能參數明顯優(yōu)于“最優(yōu)”模型,則認為原始樣本數據中存在真正的QSPR 關系,所建立的預測模型不存在“機會相關”現象。

將“y-scrambling”方法針對2 個預測模型分別運行100 次,對于MLR 模型,所得最小RMS 為5.95,對于SVM 模型,最小RMS 為4.71,兩者均分別接近原始模型預測誤差的100 倍。由此可見,只有在使用正確的應變量數據時才能產生合理的QSPR 模型,因而本文中所建立的預測模型不存在“機會相關”現象,具備較強的穩(wěn)定性。

此外,還給出了2 種模型的預測殘差分布圖(圖3、4)。從圖中可以看出,2 種模型的預測殘差都均勻且隨機分布于0基準線的兩側,可見預測模型在建立過程中未產生系統誤差。

圖3 M LR 模型爆炸下限預測殘差圖Fig.3 Plot of the residuals versus the observed values for the M LR model

圖4 SVM 模型爆炸下限預測殘差圖Fig.4 Plot of the residuals versus the observed values for the SVM model

3.5 模型的比較與分析

3.5.1 M LR 模型與SVM 模型的比較

從表3 可以看出,2 種預測模型的精度均令人滿意,其中MLR 模型的預測平均絕對誤差為0.041%,S VM 模型的預測平均絕對誤差為0.036%,均在實驗誤差允許范圍之內。這說明本文中基于遺傳算法所篩選出的4 個結構參數能夠有效地表征與烴類物質爆炸下限密切相關的結構特征,同時基于這些參數所建立的預測模型能夠成功地對未知物質的爆炸下限進行預測。

此外,從表3 還可以看出,無論是對于訓練集還是預測集,SVM 模型的性能都要明顯優(yōu)于MLR 模型。原因可能主要在于以下2 方面:一是烴類物質爆炸下限與其分子結構之間可能存在著較強的非線性關系,與線性MLR 方法相比,S VM 方法因具有強大的非線性擬合能力而體現出較強的優(yōu)越性;同時,SVM 是一種基于結構風險最小化的機器學習方法,它追求置信范圍值的最小化,而非訓練誤差的最小化,因此與MLR 相比具有更好的泛化性能。但是,SVM 方法建立的是一種“黑箱”模型,不能給出直觀的數學模型,因而無法準確了解各分子描述符對模型的貢獻值,而這些正是MLR模型的優(yōu)勢。因此,這2 種方法各有利弊,并存在一定的互補性。

3.5.2 與其他預測模型的比較

為進一步驗證本文預測模型的優(yōu)越性,將其與T.A.Albahri[10]所建立的烴類物質爆炸下限預測模型進行比較。關于樣本集的大小,Albahri 模型的樣本數(472)大于本文模型(354),然而,Albah ri 模型的外部預測樣本數(18)卻明顯少于本文模型(70),僅占其樣本總數的4%左右,不滿足QSPR 研究的一般要求[8];關于模型的輸入參數,Albahri 模型基于基團貢獻法原理以19 個分子基團作為輸入參數,而本文僅以4 個結構參數作為分子描述符,一方面使所建立的預測模型更為簡便和穩(wěn)固,另一方面因結構參數具有明確的物理意義因而便于對影響爆炸下限的重要結構特征進行解釋和研究;關于模型性能參數,T.A.Albahri 建立的非線性預測模型(0.04%)與本文MLR 模型相比(0.045%)較優(yōu),與SVM 模型(0.035%)相比則較差,這驗證了前文關于爆炸下限與分子結構之間可能存在較強非線性關系的推斷;關于模型的適用范圍,Albahri 模型基于基團貢獻法原理建立,其應用范圍受研究體系的影響較大,如果某個基團不在建模所選取的那組基團范圍之內,那么對于含有該基團的化合物就無法應用該模型進行預測,而本文預測模型則以僅需分子結構就能計算的分子結構參數作為描述符,理論上能夠對所有已知分子結構的烴類物質進行預測,適用范圍更廣。由此可見,與已有的烴類物質爆炸下限預測模型相比,本文模型在模型穩(wěn)定性、模型解釋、模型預測能力和適用范圍等方面都具有一定的優(yōu)越性。

4 結 論

應用遺傳算法從大量分子結構參數中優(yōu)化篩選出與烴類物質爆炸下限最密切相關的一組結構參數,得出了影響爆炸下限的主要結構特征為分子的大小、形狀以及分子中原子的分布信息。同時,以上述結構參數作為分子描述符,分別應用MLR 和SVM 方法建立了相應的爆炸下限預測模型,實現了根據分子結構預測爆炸下限的功能。

[1] Vidal M, Rogers W J, Holste J C,et al.A review of estimation methods for f lash points and f lammability limits[J].Process Safety Progress,2004,23(1):47-55.

[2] Katritzky A R, Lobanov V S, Karelson M.QSPR:The correlation and quantitative prediction of chemical and physical properties from structure[J].Chemical Society Review s,1995,24(4):279-287.

[3] 王連生,韓朔睽.分子結構、性質與活性[M].北京:化學工業(yè)出版社,1997.

[4] H olland J H.Adaptation in natural and artificial systems[M].Ann Arbor:University of Michigan Press,1975.

[5] Niculescu S P.Artificial neural netw orks and genetic algorithms in QSAR[J].Journal of Molecular Structure:TH EOCH EM,2003,622(l-2):71-83.

[6] Vapnik V N.The nature of statistical learning theory[M].New York:Springer,1995.

[7] 劉煥香.基于支持向量機方法的QSAR/QSPR 在化學、生物及環(huán)境科學中的應用研究[D].蘭州:蘭州大學,2005.

[8] T ropsha A,Gramatica P, Gombar V K.The importance of being earnest:Validation is the absolute essential for successful application and interpretation of QSPR models[J].QS AR&Combinatorial Science,2003,22(1):69-77.

[9] Todeschini R,Consonni V.Handbook of molecular descriptors[M].Weinheim:Wiley-VCH, 2000.

[10] Albahri T A.Flammability characteristics of pure hydrocarbons[J].Chemical Engineering Science,2003,58(16):3629-3641.