沈曉波 張智瑋 汪圣華 李學(xué)盛 牟杰 包其富

（1.華東理工大學(xué)資源與環(huán)境工程學(xué)院，上海 200237；2.浙江省應(yīng)急管理科學(xué)研究院浙江省安全工程與技術(shù)研究重點(diǎn)實(shí)驗(yàn)室，浙江 杭州 310012）

0 引言

粉塵爆炸是工業(yè)生產(chǎn)過程中較容易出現(xiàn)的危險(xiǎn)情況之一［1］。了解粉塵以及其他易燃材料的性質(zhì)對(duì)于粉塵事故防控具有重要意義。然而，通過實(shí)驗(yàn)方法來了解材料的特性不僅存在一定的危險(xiǎn)性，且成本較高。隨著醫(yī)藥、化工和新材料行業(yè)的發(fā)展，新的化學(xué)品被不斷地合成和生產(chǎn)出來，面對(duì)這些未知化學(xué)品的粉塵爆炸危險(xiǎn)性，尋找一種快速對(duì)樣品的燃爆敏感度進(jìn)行預(yù)測(cè)，實(shí)現(xiàn)未知化學(xué)品粉塵爆炸風(fēng)險(xiǎn)的快速評(píng)估，同時(shí)也能作為實(shí)驗(yàn)測(cè)試前粉塵危險(xiǎn)特性預(yù)篩的方法非常重要。定量構(gòu)效性質(zhì)關(guān)系QSPR（Quantitative Structure Property Relationship）描述了化學(xué)或生物特性與結(jié)構(gòu)屬性之間的數(shù)學(xué)關(guān)系。在安全和爆炸領(lǐng)域，QSPR 模型被證明對(duì)目標(biāo)特性（如分子單點(diǎn)能、偶極矩和分子極化指數(shù)等）具有良好的預(yù)測(cè)能力。QSPR 方法也經(jīng)常被用于藥物、材料和表面活性劑的制造以及安全性調(diào)查和其他的研究之中。

SHEN S 等［2］通過遺傳算法以及多元線性回歸方法建立了燃料混合物可燃性上限的QSPR 模型，該模型可以穩(wěn)定地預(yù)測(cè)燃料混合物的可燃性上限；朱曉等［3］基于QSPR 原理，利用遺傳-偏最小二乘（GA-PLS）方法篩選出相關(guān)描述符，采用多元線性回歸建立了鏈烷烴辛烷值的預(yù)測(cè)模型，為有效預(yù)測(cè)辛烷值提供了一種新方法。

針對(duì)可燃性粉塵，REYES O J 等［4］收集了31 種不同的化學(xué)粉塵，利用Material Studio 5.0 軟件以及遺傳算法（GFA）建立了最大超壓（Pmax）以及粉塵燃爆指數(shù)（Kst）的數(shù)學(xué)預(yù)測(cè)模型；CHAUDHARI P 等［5］利用隨機(jī)森林、人工神經(jīng)網(wǎng)絡(luò)和遺傳函數(shù)等數(shù)學(xué)方法建立了最小點(diǎn)火能預(yù)測(cè)模型并進(jìn)行優(yōu)化，使得隨機(jī)森林算法的測(cè)試集R2達(dá)到0.85，人工神經(jīng)網(wǎng)絡(luò)和遺傳函數(shù)逼近模型的測(cè)試集R2則為0.79 和0.71。

在國內(nèi)，趙琳［6］首次采用線性回歸方法計(jì)算了可燃?xì)怏w及蒸氣的分子描述符，建立了最佳多元線性回歸與啟發(fā)式回歸方法的最小點(diǎn)火能預(yù)測(cè)模型，并取得了預(yù)期的研究成果；國際上，CHAUDHARI P 等［7］首次利用隨機(jī)森林和決策樹算法建立了芳香族化合物的最小點(diǎn)火能二分類預(yù)測(cè)模型，并證明隨機(jī)森林是應(yīng)用于QSPR 中分子描述符約簡(jiǎn)的重要工具。

然而，在多數(shù)文獻(xiàn)中針對(duì)燃爆參數(shù)特性預(yù)測(cè)的研究對(duì)象多采用氣體、烴類化合物及芳香族化合物，并且多局限于使用回歸方法進(jìn)行預(yù)測(cè)，本文采用鏈?zhǔn)接袡C(jī)化合物，由于鏈?zhǔn)交衔锒酁橹咀寤衔?，并且大多?shù)情況下是由碳?xì)浣Y(jié)構(gòu)組成，體現(xiàn)出了極高的易燃性。通常情況下對(duì)于脂肪族化合物的研究只考慮了其毒性而未考慮其粉塵爆炸危險(xiǎn)性，本文利用機(jī)器學(xué)習(xí)算法進(jìn)行二分類預(yù)測(cè)，建立此類化合物分子結(jié)構(gòu)特征參數(shù)與最小點(diǎn)火能（MIE）之間的構(gòu)效關(guān)系模型，并開展模型預(yù)測(cè)能力評(píng)估和驗(yàn)證，以彌補(bǔ)此方面研究的不足，研究成果可為同類化合物MIE的快速篩選提供有力工具。

1 研究方法

1.1 數(shù)據(jù)集

各類可燃性有機(jī)化合物粉塵的燃爆參數(shù)的實(shí)驗(yàn)數(shù)據(jù)從GESTIS-DUST-EX 數(shù)據(jù)庫［8］中獲得。本文從數(shù)據(jù)庫6 973 條現(xiàn)有數(shù)據(jù)中選擇出了有最小點(diǎn)火能（MIE）記錄的37 種鏈?zhǔn)接袡C(jī)化合物進(jìn)行研究，以上述37 種鏈?zhǔn)接袡C(jī)化合物粉塵的最小點(diǎn)火能（MIE）作為目標(biāo)特性進(jìn)行預(yù)測(cè)。表1 為化合物的最小點(diǎn)火能（MIE）分布。MIE 較小的物質(zhì)點(diǎn)火敏感性較低，相對(duì)來說更加易燃甚至易爆。根據(jù)美國消防協(xié)會(huì)行業(yè)標(biāo)準(zhǔn)NFPA 652 和NFPA 77 規(guī)定，由于人體是良好的電導(dǎo)體，在正?；顒?dòng)中，人體電位可達(dá)10 ～15 kV，產(chǎn)生的火花能量可達(dá)20 ～30 mJ，從而導(dǎo)致靜電放電事件的發(fā)生，因此點(diǎn)火敏感材料的MIE 為30 mJ或更低。根據(jù)美國勞工部所著的《職業(yè)安全與健康管理局技術(shù)手冊(cè)》的第6 部分《可燃性粉塵》規(guī)定，低MIE 可燃粉塵是MIE≤30 mJ 的粉塵，如果任何化學(xué)品制造商的數(shù)據(jù)、雇主委托的MIE 測(cè)試結(jié)果或公布的數(shù)據(jù)顯示，材料的MIE為≤30 mJ或可能為≤30 mJ，CSHO 應(yīng)根據(jù)本協(xié)議將粉塵視為危險(xiǎn)采樣。所以通常情況下30 mJ 被認(rèn)為是判斷粉塵點(diǎn)火敏感性與危險(xiǎn)性的重要指標(biāo)。因此本文以30 mJ 作為MIE的分類界限，將上述37 種粉塵分成“最小點(diǎn)火能小于30 mJ”和“最小點(diǎn)火能大于30 mJ”二大類，模型中的輸出代號(hào)分別為“0”和“1”。為開展機(jī)器學(xué)習(xí)建模和驗(yàn)證，隨機(jī)選擇其中的70%用于模型開發(fā)（即訓(xùn)練集），30%用于模型驗(yàn)證（即測(cè)試集）。

表1 鏈?zhǔn)交衔镒钚↑c(diǎn)火能分布

1.2 幾何優(yōu)化及分子描述符篩選

通過Pubchem 數(shù)據(jù)庫獲取了本文所考察的37種鏈?zhǔn)交衔锏姆肿咏Y(jié)構(gòu)。在計(jì)算分子描述符之前應(yīng)進(jìn)行幾何優(yōu)化以獲得具有最低能量的分子結(jié)構(gòu)，以保證模型開發(fā)的可靠性。本文選擇在B3LYP 密度泛函理論以及6-31G（d）基組的條件下，利用Guassian 16（G16）程序?qū)Ψ肿咏Y(jié)構(gòu)進(jìn)行優(yōu)化。該方法被證實(shí)具有較佳的優(yōu)化效率［9］。隨后，將優(yōu)化后的分子結(jié)構(gòu)導(dǎo)入Material Studio 2019 軟件進(jìn)行分子描述符計(jì)算，初步篩選出352 個(gè)分子描述符。

特征分子描述符的篩選對(duì)減小模型開發(fā)規(guī)模和提高計(jì)算效率非常關(guān)鍵。首先進(jìn)行預(yù)篩選，篩除恒為0 的常數(shù)及描述符之間相關(guān)系數(shù)大于0.95 的分子描述符，初步篩選后分子描述符數(shù)量降至108 個(gè)。然后應(yīng)用隨機(jī)森林算法對(duì)分子描述符開展進(jìn)一步篩選。當(dāng)隨機(jī)森林算法應(yīng)用于特征選擇時(shí)，會(huì)在每個(gè)決策樹中計(jì)算每個(gè)特征的置換特征重要性得分。在計(jì)算出每個(gè)特征的重要性得分后，隨機(jī)森林算法會(huì)根據(jù)重要性得分從大到小對(duì)所有特征進(jìn)行排序，刪除重要性得分較低的特征。由此，最終得到9 個(gè)特征重要度較高的分子描述符，如表2 所示。

表2 分子描述符特征重要度

1.3 模型可靠性指標(biāo)

對(duì)于二分類問題，一般用精確率、準(zhǔn)確率、召回率、F1值以及PR曲線來評(píng)估二分類模型的預(yù)測(cè)結(jié)果質(zhì)量。

TP（True Positive）指正類預(yù)測(cè)為正類，F(xiàn)P（False Positive）指負(fù)類預(yù)測(cè)為正類，F(xiàn)N（False Negative）指正類預(yù)測(cè)為負(fù)類，TN（TrueNegative）指負(fù)類預(yù)測(cè)為負(fù)類。

PR 曲線是指在二分類任務(wù)中，根據(jù)學(xué)習(xí)器（分類器）的預(yù)測(cè)結(jié)果對(duì)樣例進(jìn)行排序，按此順序逐個(gè)將樣本作為正例進(jìn)行預(yù)測(cè)，每次計(jì)算出準(zhǔn)確率作為縱坐標(biāo)，召回率作為橫坐標(biāo)，得到PR 曲線（Precision-Recall Curve）。

2 結(jié)果與討論

2.1 二分類算法優(yōu)異性對(duì)比

2.1.1 不同算法的PR 曲線對(duì)比

本文中將最小點(diǎn)火能以30 mJ 為界限，對(duì)37種鏈?zhǔn)交衔锏淖钚↑c(diǎn)火能為對(duì)象進(jìn)行二分類預(yù)測(cè)。為了選出較適合粉塵燃爆數(shù)據(jù)集二元分類的算法，首先將10 種常用的二分類算法用于數(shù)據(jù)集的學(xué)習(xí)和預(yù)測(cè)，以PR 曲線結(jié)果作為初步評(píng)價(jià)標(biāo)準(zhǔn)（圖1）。

圖1 二分類模型PR 曲線

在10 種不同算法中，邏輯回歸、支持向量機(jī)、隨機(jī)森林以及梯度提升決策樹的PR 曲線的區(qū)域面積在95%以上，而其他的6 種算法在該數(shù)據(jù)集上的表現(xiàn)欠佳。

2.1.2 模型評(píng)估指標(biāo)對(duì)比

上述4 類算法PR 曲線評(píng)估較佳的二分類模型詳細(xì)比較如表3 所示。

表3 二分類模型評(píng)估分?jǐn)?shù)

二分類模型的混淆矩陣見圖2。

圖2 二分類模型的混淆矩陣

從上述評(píng)估結(jié)果可以看出，以線性函數(shù)擬合決策邊界的邏輯回歸在多特征的情況下表現(xiàn)欠佳，而核函數(shù)為（高斯）徑向基函數(shù)核（Radial Basis Function Kernel）的支持向量機(jī)在多特征的情況下預(yù)測(cè)性能良好，它與集成算法中的隨機(jī)森林和梯度提升決策樹的評(píng)估分?jǐn)?shù)都在90%以上，表現(xiàn)出了極好的預(yù)測(cè)能力與穩(wěn)定性。

從上述9 類特征中選出重要度最高的分子柔性和脂水分配系數(shù)作為二特征二分類模型評(píng)估的輸入數(shù)據(jù)，結(jié)果如表4 所示。

表4 二分類模型評(píng)估分?jǐn)?shù)

不同模型的二分類結(jié)果見圖3。

圖3 不同模型的二分類結(jié)果

從上述結(jié)果可以看出，在二特征二分類的情況下，不同算法因決策邊界不同而對(duì)預(yù)測(cè)結(jié)果產(chǎn)生影響。決策邊界是一個(gè)將屬于不同類別標(biāo)簽的數(shù)據(jù)點(diǎn)分開的表面，決策邊界對(duì)整個(gè)數(shù)據(jù)集的特征空間進(jìn)行了劃分。決策邊界可以表現(xiàn)出每個(gè)特定模型對(duì)數(shù)據(jù)集中每個(gè)數(shù)據(jù)點(diǎn)的敏感程度。

在上述4 個(gè)PR曲線評(píng)分較高的模型中，邏輯回歸和核函數(shù)為線性核（Linear kernel）的支持向量機(jī)的決策邊界是一條直線，這是因?yàn)槎卣鞯臄?shù)據(jù)分布呈現(xiàn)出良好的線性可分性，所以在此種情況下由線性函數(shù)創(chuàng)建的線性決策邊界進(jìn)行預(yù)測(cè)點(diǎn)分割體現(xiàn)出了良好的性能。而在進(jìn)行多模型比較時(shí)由于數(shù)據(jù)特征參數(shù)存在多維分布且具有非線性特征，使得邏輯回歸表現(xiàn)較差，而核函數(shù)為（高斯）徑向基函數(shù)核（Radial basis function kernel）的支持向量機(jī)在非線性數(shù)據(jù)分布中同樣體現(xiàn)出較好的穩(wěn)定性。

以決策樹為基礎(chǔ)的隨機(jī)森林和梯度提升樹的決策邊界，則更加靈活。由于決策樹的決策邊界是由訓(xùn)練樣本推斷的簡(jiǎn)單閾值規(guī)則組合構(gòu)成的，因此經(jīng)常是以平行于x 軸和y 軸的線條通過決策樹不同分枝上的隨機(jī)屬性進(jìn)行投票而劃分出的區(qū)域。隨機(jī)森林和梯度提升樹的決策邊界則是集成學(xué)習(xí)下不同決策樹疊加而產(chǎn)生的分割區(qū)域。而相比梯度提升樹，隨機(jī)森林的優(yōu)越性在于，隨機(jī)森林僅通過隨機(jī)選擇的決策樹，利用隨機(jī)特征進(jìn)行多數(shù)投票來劃分類別。每個(gè)決策樹互相獨(dú)立且互不影響，而梯度提升樹的決策樹是串行迭代過程，每一次的提升都是靠上次的預(yù)測(cè)結(jié)果與訓(xùn)練數(shù)據(jù)的差值作為新的訓(xùn)練數(shù)據(jù)進(jìn)行重新訓(xùn)練，預(yù)測(cè)結(jié)果為所有決策樹的累加值，因此決策樹之間互相影響，并且因?yàn)榉诸悩?biāo)簽數(shù)據(jù)是離散的，梯度提升樹無法直接從輸出類別去擬合類別輸出的誤差，所以分類結(jié)果會(huì)略差于隨機(jī)森林。

綜上，支持向量機(jī)和隨機(jī)森林在對(duì)上述37 種鏈?zhǔn)交衔镒钚↑c(diǎn)火能分類預(yù)測(cè)中表現(xiàn)出較優(yōu)異的預(yù)測(cè)效能和穩(wěn)定性。

2.2 分子描述符與MIE 相關(guān)性的機(jī)理分析

在所篩選出的2 種算法中，脂水分配系數(shù)和分子柔性具有較高的相關(guān)性，見圖4。脂水分配系數(shù)（LogP）描述了化合物在脂質(zhì)和水不相容的雙相體系中的溶解傾向，它可以衡量溶質(zhì)在水與脂質(zhì)中的溶解程度。主要溶于水層的溶質(zhì)表現(xiàn)為親水性，而主要溶于脂質(zhì)的溶質(zhì)表現(xiàn)為疏水性。極性大、親水性化合物的LogP較低（可以為負(fù)），非極性疏水化合物則具有更高的LogP。典型值范圍從-3（極性）到7（非極性）。因此，LogP 可以用于反映分子的極性。隨著極性的不斷增強(qiáng)，分子間作用力亦不斷增強(qiáng)。分子間作用力將物質(zhì)緊密地結(jié)合在一起。當(dāng)外部受熱平均動(dòng)能不斷增加超過分子間作用力的束縛物質(zhì)即會(huì)發(fā)生相變，所以分子間作用力更低的非極性物質(zhì)會(huì)更容易發(fā)生相變，其物質(zhì)也更易點(diǎn)燃。

分子柔性即分子鏈的鏈柔性，通常指大分子有機(jī)物分子鏈內(nèi)骨架，如C—C鍵發(fā)生內(nèi)旋轉(zhuǎn)使高分子鏈卷曲的能力，這會(huì)對(duì)分子鏈的性質(zhì)產(chǎn)生影響。分子間作用力亦會(huì)對(duì)分子鏈段的柔性產(chǎn)生限制，分子間作用力越強(qiáng)，對(duì)鏈段的柔性影響越大，分子柔韌性不斷減小，其物質(zhì)狀態(tài)也更加穩(wěn)固，從而使物質(zhì)更不容易點(diǎn)燃并具有更高的最小點(diǎn)火能值。

電子鍵反映了所連接的電子鍵基團(tuán)的影響，該電子鍵的分子描述符表現(xiàn)出了所連接的甲基（—CH2）。巴拉班指數(shù)是一類拓?fù)涿枋龇?，它可以用來描述分子的結(jié)構(gòu)與形態(tài)，將鏈?zhǔn)交衔锏姆肿咏Y(jié)構(gòu)看作是一個(gè)有n 個(gè)節(jié)點(diǎn)、m 條邊和c 個(gè)連接部件的圖來說，其定義為：

式中，y=m-n+c 是圖的等級(jí)，E（G）是邊集，Di是圖距離矩陣第i 行（或列）的所有條目之和。維納指數(shù)也是一個(gè)拓?fù)涿枋龇?，將鏈?zhǔn)交衔锟醋饔衝 個(gè)節(jié)點(diǎn)的分子結(jié)構(gòu)，其定義為：

式中，d（ij）是圖的距離矩陣。拓?fù)涿枋龇欢ǔ潭壬戏从沉朔肿咏Y(jié)構(gòu)對(duì)分子反應(yīng)的影響以及和最小點(diǎn)火能之間的相關(guān)性。

分子密度是單位體積分子的個(gè)數(shù)，單位體積內(nèi)的分子個(gè)數(shù)越多，分子間作用力越大，分子結(jié)構(gòu)也越為緊密，則其也越難被點(diǎn)燃，分子密度與最小點(diǎn)火能也表現(xiàn)出了較高的相關(guān)性。最高占有分子軌道（HOMO）的能量、最低未占有分子軌道（LUMO）的能量和也是影響最小點(diǎn)火能預(yù)測(cè)精度的重要因素。由于最高占有分子軌道和最低未占有分子軌道構(gòu)成分子軌道，這些描述符與分子的反應(yīng)性有高度相關(guān)性，因此可能也與可燃粉塵的最小點(diǎn)火能大小密切相關(guān)。

3 結(jié)論

1）本文從GESTIS-DUST-EX燃爆特征數(shù)據(jù)庫中篩選了37 種可燃性鏈?zhǔn)接袡C(jī)化合物粉塵，篩選出與最小點(diǎn)火能相關(guān)的分子描述符并進(jìn)行重要度排序，結(jié)果表明，鏈?zhǔn)接袡C(jī)化合物的脂水分配系數(shù)和分子柔性是2 個(gè)與最小點(diǎn)火能相關(guān)重要度較高的分子描述符。

2）利用常見的10 種二分類機(jī)器學(xué)習(xí)算法，對(duì)37種可燃性鏈?zhǔn)接袡C(jī)化合物粉塵的最小點(diǎn)火能進(jìn)行二分類預(yù)測(cè)應(yīng)用和比較分析發(fā)現(xiàn)，支持向量機(jī)和隨機(jī)森林表現(xiàn)出較優(yōu)異的預(yù)測(cè)效能和穩(wěn)定性。

3）從分子間作用力的角度，對(duì)脂水分配系數(shù)和分子柔性是2 個(gè)重要度較高的分子描述符與MIE之間的影響進(jìn)行機(jī)理分析解釋，為相關(guān)QSPR 模型分子描述符的選擇提供指導(dǎo)。

4）目前同類物質(zhì)的MIE 數(shù)據(jù)還非常有限，因此下一步將進(jìn)一步挖掘文獻(xiàn)數(shù)據(jù)，并開展實(shí)驗(yàn)測(cè)試，獲得更多的MIE 數(shù)據(jù)用于構(gòu)建預(yù)測(cè)模型和模型驗(yàn)證。