1 引言

撞击感度是衡量含能材料安全性能的一个重要参数^[1]，多硝基化合物在民用生产和军事上应用广泛，因此获得多硝基化合物的撞击感度数据具有一定实用价值。由于实验测试费时费力，且测试结果易受环境、操作人员、设备等因素的影响，此外，对于正处于探索阶段尚未合成出的物质而言，也无法通过实验获得测试值，因此有必要对理论预测撞击感度的方法进行研究。定量结构-性质相关性研究(Quantitative Structure Property Relationship，QSPR)是一种能根据分子结构实现有机物理化性质预测的有效方法，已广泛应用于有机物的各类理化性质及生物活性的预测中^[2]。

近年来，国内外学者已经对撞击感度与其分子结构间的相关关系进行了广泛的探讨。赵俊^[3]等人对36种硝基化合物和33种硝胺化合物的撞击感度和其分子特征量之间的关联顺序进行了研究。王睿^[4-5]等人应用电性拓扑状态指数对41种多硝基苯类化合物和46种多硝基脂肪族化合物分别进行了撞击感度的定量构效关系研究。葛素红^[6]等人计算并分析了18种三硝基芳香族炸药分子的原子化能以及原子化能与分子的电子结构能的比值与撞击感度的关联关系。Guillaume Fayet^[7]等人用最佳多元线性回归(B-MLR)研究60个硝胺类化合物的撞击感度与Codessa计算所得64个分子描述符间的线性关系。N. R. Badders^[8]等人研究了22个多硝基化合物的撞击感度与6个量子化学参数间的相关性。这些研究只专注于特定种类的样本集，所用样本较少，所得模型的适用范围有很大的局限性。因此，也有学者针对较大样本集对硝基化合物的撞击感度进行了进一步研究。袁方强^[9]等人根据分子内氢键、分子结构、对称性、氧平衡OB₁₀₀、活性指数F、“拥挤性”等分子结构参数与对156个硝基类炸药的撞击感度进行了研究。房伟^[10]等人选取氧系数、对称性、—COOR、氧杂环、苯环、α-H、α-OH、α-CH、硝基(—NO₂)、氨基(—NH₂)等作分子结构描述符，用多元线性回归计算了123个硝基化合物的撞击感度。而这些研究是凭借分子结构与撞击感度之间的定性关联，人为地发掘数个对撞击感度有影响的结构参数，因此往往构建的模型所用变量较多，或者可能存在选用参数不足以很好得表征研究的性质的问题。因此，有必要在较大样本集条件下，从分子结构角度出发对结构参数进行筛选，从而建立多硝基化合物撞击感度和其分子结构间的定量关系模型。

为了实现上述目标，这里根据QSPR基本原理，采用遗传算法^[11](Genetic Algorithm，GA)对计算所得的描述符进行筛选，分别建立多元线性回归(MLR)模型和人工神经网络^[12](ANN)模型对149种多硝基化合物的撞击感度定量构效关系进行研究，并进行模型比较。

2 样本与方法 2.1 试验样本

实验样本的数量和准确度直接影响撞击模型的适用性和预测精度。这里选取了149种多硝基化合物的撞击感度^[13]作为实验样本，其中包含47个芳香族类化合物、70个脂肪族化合物以及32个硝胺类化合物。按照4:1的比例随机选择120种作为训练集用于建立QSPR模型，剩余的29种作为测试集用于对所建模型进行外部评估。H₅₀为2.5 kg落锤实验中测得的样品爆炸几率为50%时的落高，所涉及的149种多硝基化合物的撞击感度实验值(lgH₅₀)及QSPR计算得到的预测值见附表(附表可查阅本刊网站)。由于篇幅原因，表 1只给出了部分多硝基化合物的撞击感度实验值和预测值。

2.2 研究方法

首先，借助化学软件HyperChem7.5^[14]进行分子结构的输入。在采用分子力学MM+初步优化的基础上，采用量子化学半经验方法PM3对分子结构进一步优化以获得能量低的稳定构型。所有计算均限制在Hartree-Fock能级，采用Ploak-Ribiere方法直至均方根误差梯度达到0.001 kJ/mol。然后将优化的分子结构导入Dragon2.1软件^[15]，计算获得1481种分子描述符，将近似常数或者相关系数大于0.95的描述符删除，初步筛选获得561个分子描述符。采用Material Studio 7.0软件^[16]实现遗传算法对分子描述符进行进一步筛选。继而分别采用MLR和ANN研究所得描述符与撞击感度间的定量关系，多元线性回归由SPSS软件^[17]实现，人工神经网络建模由MATLAB R2014a^[18]实现。试验步骤流程如图 1所示。

最终，将采用相关系数(R²)、均方根误差(RMSE)和平均绝对误差(MAE)等统计参数对所建立的模型性能进行全面的评估，并采用Williams图对模型的应用域进行分析。

3 结果与讨论 3.1 GA筛选结果

基于Accelrys公司的Material Studio 7.0软件，采用遗传算法(GA)对561个分子描述符进行筛选。设置方程式的初始长度6，最大方程式长度为10，种群数为50，最大代数为500，变异概率为0.1。选用Friedman提出的拟合缺失分数(Lack-Of-Fit，LOF)作为其适应度函数，比例控制LOF的平滑参数α为0.5。在遗传过程中，选择在整个种群中出现次数最多的变量为适应度函数的最优解。最终确定的6个分子描述符的名称、类型和定义见表 2。

3.2 GA-MLR模型的建立

首先，运用MLR研究GA筛选出的描述符和129种多硝基化合物撞击感度之间可能存在的线性关系。基于训练集所建立的撞击感度多元线性回归模型如下：

$ \begin{array}{l} {\rm{lg}}{H_{50}} = 16.268 - 0.285 \times IC1 - 10.656 \times ATS1{\rm{e}} - \\ \;\;\;\;\;\;\;\;\;\;\;1.578 \times MATS1{\rm{e}} + 0.031 \times RDF050v - \\ \;\;\;\;\;\;\;\;\;\;\;0.037 \times RTu + 0.326 \times R6{m_ + }\\ \;\;\;\;\;\;\;\;\;\;\;{\rm{n}} = 129,{R^2} = 0.854,{Q^2}_{{\rm{LOO}}} = 0.808,RMSE = 0.195,\\ \;\;\;\;\;\;\;\;\;\;\;{S_{\rm{D}}} = 0.200,F = 110.514,P < 0.001 \end{array} $

(1)

式中，IC1、ATS1e、MATS1e、RDF050v、RTu和R6m₊是由GA筛选获得的描述符，具体含义见表 2; n为样本数; R²为相关系数; Q_LOO²为“留一法”交互验证的相关系数; RMSE为均方根误差; S_D为标准偏差; F为F检验值; P为模型显著性概率。

该模型有较高的相关系数和交互验证系数，均方根误差和标准偏差均较小，表明该模型是可靠的。模型的F检验值远大于F_理论(95%置信度)=2.198，因此可以认为该回归方程及所选变量的影响均是显著的。表 3所示是模型的显著性检验结果，标准系数可以确定不同分子描述符对模型的影响程度。从表 3可以看出，6个描述符对撞击感度的影响大小排序为ATS1e＞MATS1e＞RTu＞RDF050v＞IC1＞R6m₊，所有描述符与撞击感度lgH₅₀均呈负相关性。

为检验模型的外部预测能力，应用测试集对所建模型进行外部验证。计算得到，训练集的相关系数R²为0.854，均方根误差RMSE为0.195，平均绝对误差MAE为0.157;测试集的相关系数R_EXT²为0.790，均方根误差RMSE_EXT为0.215，平均绝对误差MAE为0.183。模型的主要性能参数见表 4，可以看出，无论对于训练集还是测试集，模型均有较小的均方根误差和平均绝对误差，说明该线性模型有着较好的内部稳定性。

训练集和测试集撞击感度预测值与实验值的比较及残差分布分别见图 2和图 3。如图 2所示，所有的数据点均随机分布在对角线附近，表明此线性模型具有一定的预测能力，但是模型的精度不高，预测结果不是很好。由此猜测，多硝基类化合物的撞击感度和分子结构间可能存在非线性关系，因此将尝试建立GA-ANN非线性模型进行研究。而由图 3可知，模型的预测残差均匀且随机分布于基准线的两侧。由此可以认为，GA-MLR模型在建立过程中未产生系统误差。

3.3 GA-ANN模型的建立

在用神经网络研究化学结构性质的领域，三层反向传播的前馈神经网络应用最为广泛。理论上已经证明：具有偏差和至少一个S型隐含层加上一个线性输出层的网络，能够逼近任何有理函数^[12]。因此，这里建立包含输入层、输出层和一层输隐含层的三层BP人工神经网络用于非线性模型研究。BP神经网络的构建在MATLAB R2014a环境中实现，经过多次训练试验，最终确定6-8-1的神经网络结构的精度最好。具体参数设定如下：输入层节点数为6，隐含层节点数为8，输出层节点数为1，选择最小训练速率0.1，最大迭代次数为1000，隐含层激活转移函数采用双曲正切S型函数tansig，输出层激活转移函数采用线性函数purelin，训练函数采用trainlm，网络学习函数采用带动量项的BP学习规则learngdm。

模型的预测结果见附表，相关模型性能参数见表 4。计算得到，该模型的相关系数为0.974，交叉验证系数为0.975，训练集和测试集的均方根误差分别为0.071和0.248，平均绝对误差分别为0.051和0.197。从相关系数、交互验证系数和均方根误差来看，GA-ANN模型的精度明显高于GA-MLR模型。

如图 4所示，所有多硝基化合物撞击感度预测值和实验值的比较均未偏离对角线很远，整体基本紧挨对角线，这说明该模型比MLR建立的模型精度高。图 5给出了GA-ANN模型的残差图，预测残差均匀且随机分布于基准线，该模型在建立过程中未产生系统误差。

3.4 模型比较

附表给出了两个模型对所有149种多硝基化合物的预测值。为便于比较分析GA-MLR和GA-ANN模型，计算得到两个模型的主要性能参数见表 4。比较两个模型的相关系数R²、均方根误差RMSE、平均绝对误差MAE、平均绝对误差百分比MAPE和交互验证系数Q_LOO²可以发现，GA-ANN模型的性能参数都要优于GA-MLR模型，这表明非线性模型的精度要明显高于线性模型。然而，比较两个模型的测试集相关系数R_EXT²和外部验证系数Q_F1²、Q_F2²和Q_F3²，GA-MLR模型均略优于GA-ANN模型，这表明非线性模型的外部预测能力相比线性模型并没有得到预期的提升。即，所建的两个模型虽然具有一定的预测能力，但是预测结果并没有取得理想的预期效果。这可能由于两个原因：(1)样本集的数据存在问题; (2)建模所用的参数并不能很好表征多硝基化合物的撞击感度。

这里采用Williams图^[19]进一步分析比较两个模型的应用域。横坐标为臂比值，纵坐标为标准残差。标准残差落在(-3，+3)以外的化合物，认为其实验值是离群点; 当化合物的臂比值h_i大于警戒值h^*时，认为该物质显著影响模型的回归效果^[14]。臂比值的计算式如下：

$ {h_i} = x_i^T\left( {X \cdot {X^T}} \right){x_i}\;\;\;(i = 1, \ldots, n) $

(2)

警戒值h^*的表达式如下：

$ {h^*} = 3 \times \left( {m + 1} \right)/n $

(3)

式中，m为建立模型所用变量的个数。

如图 6、图 7所示的Williams图中，GA-MLR模型和GA-ANN模型分别有2个物质和5个物质的标准残差落在了基准线以外，其相应的化合物和撞击感度值见表 5和表 6。

观察表 5、表 6发现，出现异常值的物质其撞击感度的预测残差均较大，这表明，建模所用的参数没能很好的表征其撞击感度。已有研究表明^[20]，除了分子结构外，含能材料的感度还与炸药晶体特性密切相关，如颗粒内部空洞的微观结构、颗粒表面缺陷、颗粒大小及分布、颗粒形状和表面光滑度对感度影响非常显著。因此，仅用分子结构参数来表征撞击感度还不够完善。而2′, 2′-二硝基丙基-2, 4, 6-三硝基苯甲酸酯的标准残差在两个模型中同时出现异常，可认为主要原因是其实验值数据本身不够准确。在表 7中，六硝基苯的臂比值超出了警戒值，这表明，六硝基苯显著影响了模型的回归效果，可以认为其结构相对于整个样本集而言是比较特殊的。

将所建的两个模型与几种已有的多硝基化合物撞击感度QSPR模型进行比较，相关统计参数见表 8。

本文建模采用的是由GA筛选获得分子描述符，因此能更系统全面地表征所研究的性质，且模型中涉及的变量个数远少于已有模型。比较各模型的相关系数、交互验证系数和外部验证的均方根误差，可以认为这里获得的GA-MLR模型和GA-ANN模型明显优于已有模型。

4 结论

(1) 运用遗传算法筛选出6个分子描述符，分别构建了149种多硝基化合物撞击感度的GA-MLR模型和GA-ANN模型，结果表明两个模型均是稳定可靠的;

(2) 所得GA-ANN模型明显优于GA-MLR模型，说明多硝基化合物的撞击感度与分子结构间存在较强的非线性关系，且两个模型的预测精度均高于已有的QSPR模型;

(3) 由于落锤实验受诸如晶体特性、环境等其他因素的影响，因此所建的两个模型虽然具有一定的预测能力，但是预测效果仍有继续完善的空间。