基于SMA-GBDT-3D CA 的土壤污染扩散模拟方法*

于 淼 王占刚

（北京信息科技大学信息与通信工程学院 北京 100101）

1 引言

土壤污染扩散是一个复杂的理化过程，针对其扩散过程的研究多是利用数值计算的方法来构建数学模型，且大部分以城市表层土壤为研究目标［1～3］。郭鸿等［4］根据重金属污染物在土壤中的扩散规律建立了重金属污染物扩散方程；王建等［5］基于菲克定律及对流—弥散模型构建土壤重金属污染扩散模型，并利用GIS平台进行模拟可视化；张静［6］考虑大气沉降对土壤中重金属污染物传播的影响建立了点源扩散的高斯模型。但数学模型结构复杂且求解较为困难，鉴于土壤污染是一个复杂、动态的过程，而元胞自动机是一个通过局部转换规则来实现整体变化的动力学模型，能够清楚、准确、完整地模拟复杂的变化过程。已有学者结合逻辑回归与元胞自动机预测土壤重金属污染情况［7］，得到较为可靠的空间分布，但影响因素冗余造成模型效率低并且容易出现过拟合问题。此外，针对区域土壤污染扩散过程的研究多是在二维空间上进行模拟的。

特征选择可以在减少特征数量的同时保证不变甚至获得更优异的学习性能。梯度提升决策树（GBDT）［8］的开发是为了识别主要特征并且有效地进行特征选择。Wang H 等［9］利用GBDT 构建特征选择的集成模型用于质量审查检测。但是，GBDT在构建决策树时需要多次计算特征分支点的信息增益，冗余的初始输入会导致时间和空间效率低下。然而，黏菌优化算法（Slime Mould Algorithm，SMA）［10］的收敛效率高、适应性强的特点可以优化GBDT模型。

综上所述，本文提出一种基于黏菌优化算法（SMA）优化梯度提升决策树（GBDT）并结合三维元胞自动机（3D Cellular Automat，3D CA）的土壤重金属点源污染扩散模型。为优化GBDT 的初始输入，利用其分类精度来评价输入各个特征的质量，由此减少SMA 陷入局部最优的可能，使得算法实现全局优化并删除相关性低的特征，在提高运算效率的同时保证模型的准确性。为了提高模拟算法的效率及精度，在传统三维元胞自动机的基础上增加自适应调节机制，在模拟扩散的过程中确定有限的扩散空间从而优化SMA-GBDT-3D CA算法。

2 研究方法

2.1 梯度提升决策树

梯度提升决策树（GBDT）是一种迭代决策树算法，它是由多个决策树组成的加法模型。构建决策树使得残差向梯度方向递减［11］，GBDT 模型可以表示为式（1）：

式中，x是输入特征，h(x)为基本学习器，M为其数量，γm是每个学习器的参数。

1）初始化分类回归树：

2）迭代m=1:M次，计算残差的梯度方向：

3）根据最小二乘法计算模型参数：

4）计算当前决策树的权重：

5）更新模型：

2.2 黏菌优化算法

黏菌优化算法（SMA）是通过模拟粘菌觅食行为进行搜索的元启发算法［12］。设置权重来模拟粘菌在觅食过程中产生的正负反馈，从而形成三种不同的形态。黏菌的位置更新规则如下：

其中，LB和UB分别为搜索的上、下界，rand和r为[0,1]区间的随机数，是范围[-a,a]间的随机参数，从1到0线性下降，t为当前迭代，表示当前最佳个体位置，和为黏菌中两个随机选择的个体，W→是黏菌重量，S(i)为适应度（i∈1,2,…,N），DF代表最佳适应度。参数a及权重表达式如下：

式中，condition表示S(i)排前一半的菌群，bF和wF分别表示当前迭代最优和最差适应度，SI代表适应度序列。

2.3 三维元胞自动机模型

元胞自动机（Cellular Automata，CA）是一种时间维度、空间维度以及状态都离散的动力学模型。每个离散状态的元胞都分散在空间网格中并遵循着一定的规律进行同步更新［13］。

每个元胞的位置用(i,j,k)表示，以表示其在t时步的元胞状态。在水平方向上采用周期型边界条件；垂直方向采取反射型边界条件。每个时步t 内，中心元胞与其邻居进行污染物传输，t+1时步元胞的状态计算公式如式（10）：

其中，为t+1 时刻(i,j,k)位置土壤中重金属浓度；和分别表示水平方向和垂直方向上的污染物扩散系数。

3 SMA-GBDT-3D CA模型及其优化

3.1 SMA-GBDT-3D CA模型

不同于液体与空气，土壤在空间分布上具有不连续性，因此选定任一空间立方体为中心元胞，与之相邻的10个元胞为其邻居，如图1所示。

图1 三维元胞邻居

由于大气中的大多数重金属污染物会沉降在土壤中，因此，针对区域表层土壤，在考虑污染物自身扩散作用的基础上进一步考虑风场对污染扩散的影响则有：

式中，为自身浓度扩散系数，为风场影响的扩散系数。

理想状态下污染物的传播由污染源向四周近似圆形地缓慢扩散［14］。然而污染扩散受到各种自然、人为因素的影响空间形状会产生偏移。为了获得风向、风速、高程等因素对各个方向扩散的影响程度，引入SMA优化的GBDT算法。

假设D 为初始数据集的维度，即属性数量，可行解集为xi=(xi1,xi2,…,xiD),i={1,2,…,N}，其中xij(j∈1,2,…,N)对应解空间中原始数据集的特征。可行解表示选择的特征子集，以“1”和“0”表示选择或不选择该特征。

SMA-GBDT-3D CA算法流程如图2所示。

图2 SMA-GBDT-3D CA算法流程图

通过SMA 来优化GBDT的输入，以减少不相关的特征提高模型效率，再使用GBDT 进一步减少数据集的特征并获取特征重要性作为各属性及方向的扩散系数，根据规则进行扩散过程模拟。算法具体步骤如下：

1）获取研究区相应数据并进行数据处理；

2）建立三维元胞空间，初始化元胞状态；

3）初始化黏菌种群，计算适应度，更新参数；

4）更新个体位置及全局最优解，判断是否达到最大迭代，如果没达到返回计算适应度；

5）获取特征子集作为GBDT 的输入，使用CART 决策树及平方损失函数。初始化GBDT 模型；

6）由式（3）计算残差梯度方向，拟合模型；

7）在不断迭代过程中构建完整的梯度提升决策树模型；

9）分析变量对于各个方向扩散的贡献度并确定扩散系数；

10）根据状态转换规则更新元胞状态；

11）利用循环函数实现污染扩散过程的动态模拟与显示；

12）判断元胞状态是否达到阈值，如果没达到则继续更新元胞状态，否则停止模拟。

3.2 SMA-GBDT-3D CA模型边界反弹机制

元胞自动机是一种并行的动力学模型，每一时步，元胞空间中的所有元胞同时更新状态。当污染未扩散至整个元胞空间时，每一时步扩散半径增长1个元胞；当所有元胞的状态都不为0时，即污染物遍布元胞空间，此时所有元胞继续同步更新元胞状态。此种更新机制每一次需要更新整个元胞空间，当研究区域较大时，模型实现的效率较低且不符合重金属污染物易吸附于土壤的实际情况。为了增加模型模拟的效率及准确率，提出模型扩散范围边界反弹机制。具体流程如图3所示。

图3 模型边界反弹机制

首先每次更新状态前确定当前元胞边界及范围，再根据算法获得的规则更新当前范围内的元胞空间。设定边界元胞扩散阈值，当边界元胞最大值超过阈值则增加扩散范围，即扩散半径加一。由污染扩散分布规律可知，一般污染源周边污染物浓度呈正态分布，因此，当max(边界元胞)未超出阈值时，按照与污染源距离的正态分布分配超出边界阈值的部分，距离污染源越近的值越大。由此，可减少每次更新的元胞数量，提高模型效率及模拟精度。

4 实验及结果分析

4.1 数据来源及预处理

本实验的数据来源于华北某区的土壤重金属采样数据，研究区采样点分布如图4 所示。研究区域土壤中As、Cr、Cu、Cd、Hg、Ni、Pb、Zn 的含量如表1。

图4 研究区域采样点分布

表1 研究区域土壤中重金属含量

利用SPSS 软件对研究区域采样的八种土壤重金属进行主成分分析，得到了两种主成分，在第一主成分中Cr、Pb、Hg、Ni、As 有较高载荷，表明这五种重金属污染因子含有原始数据较多的信息量，其对土壤污染结果影响较大。区域土壤重金属点源污染扩散模拟以重金属Hg为例。

对研究区域的145个采样点进行克里金插值［15］，根据污染核密度图选定点源污染区域并获得数据矩阵。同时，将研究区域划分为同样大小的空间网格并用三维元胞表示。将研究区域的风场、高程等数据同样进行插值处理，获得经纬度、风向、风速等矢量数据。在构建模型前需要将风向场、含量值等数据归一化处理以便模型计算，公式如式（12）。

式中，Ni为归一化后的标准值，xi是原始数据，xmax和xmin分别是原始数据的最大和最小值。最终获得经纬度、ph、风向、风速、距离、高程等共9 种属性及对应的Hg浓度值。

4.2 模拟结果分析

对7056 组数据进行随机分组，令70%的数据用于训练，30%用于验证。利用SMA 算法优化GBDT 输入，在保证准确率的同时增加效率，智能优化算法全部迭代100 次，对比传统的优化算法以及近年提出的平衡优化器（EO）［16］算法优化GBDT 得到的特征选择结果如图5。

图5 模型效率及属性筛选对比

图中，方法1～6 分别为SMA-GBDT、PSO-GBDT、ABC-GBDT、GWO-GBDT、EO-GBDT 以及原始GBDT算法，将6种方法分别进行属性筛选及拟合，获得模型效率及筛选后剩余属性对比。其中，为了更加清晰展示运行效率比较，以SMA-GBDT 模型为基准，根据公式来进行各个方法的效率比较，式中T为各方法的运行时间。由图可以看出，在属性筛选能力上，SMA-GBDT 与其他两种方法均减少至剩余3 种属性，在特征选择方面表现良好。在算法效率方面，相比直接将数据集输入GBDT 模型中，优化后的模型均提高了运行效率，但很明显，SMA-GBDT算法的表现最为优异。

现工程及研究中很多方法以拟合精度为代价来达到减少特征、避免维数灾难的目的。为了定量地分析模型的拟合结果与实际情况的差距，通过均方根误差（RMSE）、决定系数（R2）对结果进行评估，指标计算公式如下：

式中，cr(x,y,z)和cs(x,y,z)分别为同一位置的污染物浓度真实值和模拟值，为真实值的平均值，n 是样本数据的数量。计算结果如表2 所示。

表2 模拟结果误差对比

由上表对比可见，SMA-GBDT 算法的拟合效果要优于单一的GBDT 算法，并且SMA 优化GBDT获得的特征子集基本代表了整个数据集的特征，能够很好地解释变量。综上可以看出，SMA-GBDT算法的综合能力较强，在保证精度甚至增强拟合效果的同时大幅提高模型效率。

计算SMA-GBDT 模型的特征重要程度以确定表层土壤污染扩散系数，影响因素及其重要性如表3所示。

表3 影响因素重要性

为清晰看出研究区域污染分布情况，从SMA-GBDT-3D CA 算法的模拟结果中随机选取3条线上20 个点位的Hg 含量值与真实值的克里金插值及传统元胞自动机的模拟结果进行对比，对比结果如图6所示。

图6 插值与模拟值对比

从图中可以看出，研究算法的模拟结果相较于传统元胞自动机模型更加逼近克里金插值得到的数据曲线，模拟的精度更好。每条线可以反映污染扩散后重金属Hg 的分布，图中三条曲线的变化趋势基本相同，说明最终模拟的污染物分布情况与真实分布相符。

为了清晰展示研究区土壤重金属污染扩散过程，在模型模拟结果的基础上，将三维可视化、WebGIS 与土壤重金属点源污染扩散过程结合起来，在Cesium 平台中进行污染物扩散过程的三维动态展示，如图7所示。

图7 污染扩散三维可视化

将污染分布核密度图叠加在地图上，由深色到浅色表示污染物Hg 浓度由高到低的情况，污染核密度图清晰地描述出污染位置、范围，三维地形模型中还展示了研究区线、地形地势等信息。受到风向、位置等影响的点源污染扩散过程用箭头来表示，通过箭头向四周的移动表示污染扩散的动态过程，箭头起始位置为污染源，结束位置为模拟出的污染范围，可以看出由于主导风向以及地理位置的影响，污染物扩散路径从八个方向产生偏移。扩散过程动态可视化清晰地展示出污染物从污染中心向四周逐渐扩散，且距离污染源越远污染物浓度越低的扩散规律。

5 结语

本文提出了边界反弹机制的SMA-GBDT-3D CA污染扩散模型，针对土壤重金属污染问题，利用SMA 优化的GBDT 算法获得扩散系数，建立三维元胞自动机模型获得扩散模型。实验表明，提出的方法在保证模拟扩散状态及分布精度的同时增加了运行效率，模型的综合能力表现良好，因此该方法可以作为土壤重金属点源污染扩散三维模拟的有效方法。

土壤中重金属的污染扩散还会受到其他污染源、人为活动、植被条件等因素的影响，在以后的研究中，还需要将这些影响因素考虑到模型中，得到更加准确的扩散规则来提高模型精度。