基于SMA-GBDT-3D CA 的土壤污染擴(kuò)散模擬方法*

于 淼 王占剛

（北京信息科技大學(xué)信息與通信工程學(xué)院 北京 100101）

1 引言

土壤污染擴(kuò)散是一個(gè)復(fù)雜的理化過程，針對其擴(kuò)散過程的研究多是利用數(shù)值計(jì)算的方法來構(gòu)建數(shù)學(xué)模型，且大部分以城市表層土壤為研究目標(biāo)［1～3］。郭鴻等［4］根據(jù)重金屬污染物在土壤中的擴(kuò)散規(guī)律建立了重金屬污染物擴(kuò)散方程；王建等［5］基于菲克定律及對流—彌散模型構(gòu)建土壤重金屬污染擴(kuò)散模型，并利用GIS平臺進(jìn)行模擬可視化；張靜［6］考慮大氣沉降對土壤中重金屬污染物傳播的影響建立了點(diǎn)源擴(kuò)散的高斯模型。但數(shù)學(xué)模型結(jié)構(gòu)復(fù)雜且求解較為困難，鑒于土壤污染是一個(gè)復(fù)雜、動(dòng)態(tài)的過程，而元胞自動(dòng)機(jī)是一個(gè)通過局部轉(zhuǎn)換規(guī)則來實(shí)現(xiàn)整體變化的動(dòng)力學(xué)模型，能夠清楚、準(zhǔn)確、完整地模擬復(fù)雜的變化過程。已有學(xué)者結(jié)合邏輯回歸與元胞自動(dòng)機(jī)預(yù)測土壤重金屬污染情況［7］，得到較為可靠的空間分布，但影響因素冗余造成模型效率低并且容易出現(xiàn)過擬合問題。此外，針對區(qū)域土壤污染擴(kuò)散過程的研究多是在二維空間上進(jìn)行模擬的。

特征選擇可以在減少特征數(shù)量的同時(shí)保證不變甚至獲得更優(yōu)異的學(xué)習(xí)性能。梯度提升決策樹（GBDT）［8］的開發(fā)是為了識別主要特征并且有效地進(jìn)行特征選擇。Wang H 等［9］利用GBDT 構(gòu)建特征選擇的集成模型用于質(zhì)量審查檢測。但是，GBDT在構(gòu)建決策樹時(shí)需要多次計(jì)算特征分支點(diǎn)的信息增益，冗余的初始輸入會(huì)導(dǎo)致時(shí)間和空間效率低下。然而，黏菌優(yōu)化算法（Slime Mould Algorithm，SMA）［10］的收斂效率高、適應(yīng)性強(qiáng)的特點(diǎn)可以優(yōu)化GBDT模型。

綜上所述，本文提出一種基于黏菌優(yōu)化算法（SMA）優(yōu)化梯度提升決策樹（GBDT）并結(jié)合三維元胞自動(dòng)機(jī)（3D Cellular Automat，3D CA）的土壤重金屬點(diǎn)源污染擴(kuò)散模型。為優(yōu)化GBDT 的初始輸入，利用其分類精度來評價(jià)輸入各個(gè)特征的質(zhì)量，由此減少SMA 陷入局部最優(yōu)的可能，使得算法實(shí)現(xiàn)全局優(yōu)化并刪除相關(guān)性低的特征，在提高運(yùn)算效率的同時(shí)保證模型的準(zhǔn)確性。為了提高模擬算法的效率及精度，在傳統(tǒng)三維元胞自動(dòng)機(jī)的基礎(chǔ)上增加自適應(yīng)調(diào)節(jié)機(jī)制，在模擬擴(kuò)散的過程中確定有限的擴(kuò)散空間從而優(yōu)化SMA-GBDT-3D CA算法。

2 研究方法

2.1 梯度提升決策樹

梯度提升決策樹（GBDT）是一種迭代決策樹算法，它是由多個(gè)決策樹組成的加法模型。構(gòu)建決策樹使得殘差向梯度方向遞減［11］，GBDT 模型可以表示為式（1）：

式中，x是輸入特征，h(x)為基本學(xué)習(xí)器，M為其數(shù)量，γm是每個(gè)學(xué)習(xí)器的參數(shù)。

1）初始化分類回歸樹：

2）迭代m=1:M次，計(jì)算殘差的梯度方向：

3）根據(jù)最小二乘法計(jì)算模型參數(shù)：

4）計(jì)算當(dāng)前決策樹的權(quán)重：

5）更新模型：

2.2 黏菌優(yōu)化算法

黏菌優(yōu)化算法（SMA）是通過模擬粘菌覓食行為進(jìn)行搜索的元啟發(fā)算法［12］。設(shè)置權(quán)重來模擬粘菌在覓食過程中產(chǎn)生的正負(fù)反饋，從而形成三種不同的形態(tài)。黏菌的位置更新規(guī)則如下：

其中，LB和UB分別為搜索的上、下界，rand和r為[0,1]區(qū)間的隨機(jī)數(shù)，是范圍[-a,a]間的隨機(jī)參數(shù)，從1到0線性下降，t為當(dāng)前迭代，表示當(dāng)前最佳個(gè)體位置，和為黏菌中兩個(gè)隨機(jī)選擇的個(gè)體，W→是黏菌重量，S(i)為適應(yīng)度（i∈1,2,…,N），DF代表最佳適應(yīng)度。參數(shù)a及權(quán)重表達(dá)式如下：

式中，condition表示S(i)排前一半的菌群，bF和wF分別表示當(dāng)前迭代最優(yōu)和最差適應(yīng)度，SI代表適應(yīng)度序列。

2.3 三維元胞自動(dòng)機(jī)模型

元胞自動(dòng)機(jī)（Cellular Automata，CA）是一種時(shí)間維度、空間維度以及狀態(tài)都離散的動(dòng)力學(xué)模型。每個(gè)離散狀態(tài)的元胞都分散在空間網(wǎng)格中并遵循著一定的規(guī)律進(jìn)行同步更新［13］。

每個(gè)元胞的位置用(i,j,k)表示，以表示其在t時(shí)步的元胞狀態(tài)。在水平方向上采用周期型邊界條件；垂直方向采取反射型邊界條件。每個(gè)時(shí)步t 內(nèi)，中心元胞與其鄰居進(jìn)行污染物傳輸，t+1時(shí)步元胞的狀態(tài)計(jì)算公式如式（10）：

其中，為t+1 時(shí)刻(i,j,k)位置土壤中重金屬濃度；和分別表示水平方向和垂直方向上的污染物擴(kuò)散系數(shù)。

3 SMA-GBDT-3D CA模型及其優(yōu)化

3.1 SMA-GBDT-3D CA模型

不同于液體與空氣，土壤在空間分布上具有不連續(xù)性，因此選定任一空間立方體為中心元胞，與之相鄰的10個(gè)元胞為其鄰居，如圖1所示。

圖1 三維元胞鄰居

由于大氣中的大多數(shù)重金屬污染物會(huì)沉降在土壤中，因此，針對區(qū)域表層土壤，在考慮污染物自身擴(kuò)散作用的基礎(chǔ)上進(jìn)一步考慮風(fēng)場對污染擴(kuò)散的影響則有：

式中，為自身濃度擴(kuò)散系數(shù)，為風(fēng)場影響的擴(kuò)散系數(shù)。

理想狀態(tài)下污染物的傳播由污染源向四周近似圓形地緩慢擴(kuò)散［14］。然而污染擴(kuò)散受到各種自然、人為因素的影響空間形狀會(huì)產(chǎn)生偏移。為了獲得風(fēng)向、風(fēng)速、高程等因素對各個(gè)方向擴(kuò)散的影響程度，引入SMA優(yōu)化的GBDT算法。

假設(shè)D 為初始數(shù)據(jù)集的維度，即屬性數(shù)量，可行解集為xi=(xi1,xi2,…,xiD),i={1,2,…,N}，其中xij(j∈1,2,…,N)對應(yīng)解空間中原始數(shù)據(jù)集的特征?？尚薪獗硎具x擇的特征子集，以“1”和“0”表示選擇或不選擇該特征。

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

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

通過SMA 來優(yōu)化GBDT的輸入，以減少不相關(guān)的特征提高模型效率，再使用GBDT 進(jìn)一步減少數(shù)據(jù)集的特征并獲取特征重要性作為各屬性及方向的擴(kuò)散系數(shù)，根據(jù)規(guī)則進(jìn)行擴(kuò)散過程模擬。算法具體步驟如下：

1）獲取研究區(qū)相應(yīng)數(shù)據(jù)并進(jìn)行數(shù)據(jù)處理；

2）建立三維元胞空間，初始化元胞狀態(tài)；

3）初始化黏菌種群，計(jì)算適應(yīng)度，更新參數(shù)；

4）更新個(gè)體位置及全局最優(yōu)解，判斷是否達(dá)到最大迭代，如果沒達(dá)到返回計(jì)算適應(yīng)度；

5）獲取特征子集作為GBDT 的輸入，使用CART 決策樹及平方損失函數(shù)。初始化GBDT 模型；

6）由式（3）計(jì)算殘差梯度方向，擬合模型；

7）在不斷迭代過程中構(gòu)建完整的梯度提升決策樹模型；

9）分析變量對于各個(gè)方向擴(kuò)散的貢獻(xiàn)度并確定擴(kuò)散系數(shù)；

10）根據(jù)狀態(tài)轉(zhuǎn)換規(guī)則更新元胞狀態(tài)；

11）利用循環(huán)函數(shù)實(shí)現(xiàn)污染擴(kuò)散過程的動(dòng)態(tài)模擬與顯示；

12）判斷元胞狀態(tài)是否達(dá)到閾值，如果沒達(dá)到則繼續(xù)更新元胞狀態(tài)，否則停止模擬。

3.2 SMA-GBDT-3D CA模型邊界反彈機(jī)制

元胞自動(dòng)機(jī)是一種并行的動(dòng)力學(xué)模型，每一時(shí)步，元胞空間中的所有元胞同時(shí)更新狀態(tài)。當(dāng)污染未擴(kuò)散至整個(gè)元胞空間時(shí)，每一時(shí)步擴(kuò)散半徑增長1個(gè)元胞；當(dāng)所有元胞的狀態(tài)都不為0時(shí)，即污染物遍布元胞空間，此時(shí)所有元胞繼續(xù)同步更新元胞狀態(tài)。此種更新機(jī)制每一次需要更新整個(gè)元胞空間，當(dāng)研究區(qū)域較大時(shí)，模型實(shí)現(xiàn)的效率較低且不符合重金屬污染物易吸附于土壤的實(shí)際情況。為了增加模型模擬的效率及準(zhǔn)確率，提出模型擴(kuò)散范圍邊界反彈機(jī)制。具體流程如圖3所示。

圖3 模型邊界反彈機(jī)制

首先每次更新狀態(tài)前確定當(dāng)前元胞邊界及范圍，再根據(jù)算法獲得的規(guī)則更新當(dāng)前范圍內(nèi)的元胞空間。設(shè)定邊界元胞擴(kuò)散閾值，當(dāng)邊界元胞最大值超過閾值則增加擴(kuò)散范圍，即擴(kuò)散半徑加一。由污染擴(kuò)散分布規(guī)律可知，一般污染源周邊污染物濃度呈正態(tài)分布，因此，當(dāng)max(邊界元胞)未超出閾值時(shí)，按照與污染源距離的正態(tài)分布分配超出邊界閾值的部分，距離污染源越近的值越大。由此，可減少每次更新的元胞數(shù)量，提高模型效率及模擬精度。

4 實(shí)驗(yàn)及結(jié)果分析

4.1 數(shù)據(jù)來源及預(yù)處理

本實(shí)驗(yàn)的數(shù)據(jù)來源于華北某區(qū)的土壤重金屬采樣數(shù)據(jù)，研究區(qū)采樣點(diǎn)分布如圖4 所示。研究區(qū)域土壤中As、Cr、Cu、Cd、Hg、Ni、Pb、Zn 的含量如表1。

圖4 研究區(qū)域采樣點(diǎn)分布

表1 研究區(qū)域土壤中重金屬含量

利用SPSS 軟件對研究區(qū)域采樣的八種土壤重金屬進(jìn)行主成分分析，得到了兩種主成分，在第一主成分中Cr、Pb、Hg、Ni、As 有較高載荷，表明這五種重金屬污染因子含有原始數(shù)據(jù)較多的信息量，其對土壤污染結(jié)果影響較大。區(qū)域土壤重金屬點(diǎn)源污染擴(kuò)散模擬以重金屬Hg為例。

對研究區(qū)域的145個(gè)采樣點(diǎn)進(jìn)行克里金插值［15］，根據(jù)污染核密度圖選定點(diǎn)源污染區(qū)域并獲得數(shù)據(jù)矩陣。同時(shí)，將研究區(qū)域劃分為同樣大小的空間網(wǎng)格并用三維元胞表示。將研究區(qū)域的風(fēng)場、高程等數(shù)據(jù)同樣進(jìn)行插值處理，獲得經(jīng)緯度、風(fēng)向、風(fēng)速等矢量數(shù)據(jù)。在構(gòu)建模型前需要將風(fēng)向場、含量值等數(shù)據(jù)歸一化處理以便模型計(jì)算，公式如式（12）。

式中，Ni為歸一化后的標(biāo)準(zhǔn)值，xi是原始數(shù)據(jù)，xmax和xmin分別是原始數(shù)據(jù)的最大和最小值。最終獲得經(jīng)緯度、ph、風(fēng)向、風(fēng)速、距離、高程等共9 種屬性及對應(yīng)的Hg濃度值。

4.2 模擬結(jié)果分析

對7056 組數(shù)據(jù)進(jìn)行隨機(jī)分組，令70%的數(shù)據(jù)用于訓(xùn)練，30%用于驗(yàn)證。利用SMA 算法優(yōu)化GBDT 輸入，在保證準(zhǔn)確率的同時(shí)增加效率，智能優(yōu)化算法全部迭代100 次，對比傳統(tǒng)的優(yōu)化算法以及近年提出的平衡優(yōu)化器（EO）［16］算法優(yōu)化GBDT 得到的特征選擇結(jié)果如圖5。

圖5 模型效率及屬性篩選對比

圖中，方法1～6 分別為SMA-GBDT、PSO-GBDT、ABC-GBDT、GWO-GBDT、EO-GBDT 以及原始GBDT算法，將6種方法分別進(jìn)行屬性篩選及擬合，獲得模型效率及篩選后剩余屬性對比。其中，為了更加清晰展示運(yùn)行效率比較，以SMA-GBDT 模型為基準(zhǔn)，根據(jù)公式來進(jìn)行各個(gè)方法的效率比較，式中T為各方法的運(yùn)行時(shí)間。由圖可以看出，在屬性篩選能力上，SMA-GBDT 與其他兩種方法均減少至剩余3 種屬性，在特征選擇方面表現(xiàn)良好。在算法效率方面，相比直接將數(shù)據(jù)集輸入GBDT 模型中，優(yōu)化后的模型均提高了運(yùn)行效率，但很明顯，SMA-GBDT算法的表現(xiàn)最為優(yōu)異。

現(xiàn)工程及研究中很多方法以擬合精度為代價(jià)來達(dá)到減少特征、避免維數(shù)災(zāi)難的目的。為了定量地分析模型的擬合結(jié)果與實(shí)際情況的差距，通過均方根誤差（RMSE）、決定系數(shù)（R2）對結(jié)果進(jìn)行評估，指標(biāo)計(jì)算公式如下：

式中，cr(x,y,z)和cs(x,y,z)分別為同一位置的污染物濃度真實(shí)值和模擬值，為真實(shí)值的平均值，n 是樣本數(shù)據(jù)的數(shù)量。計(jì)算結(jié)果如表2 所示。

表2 模擬結(jié)果誤差對比

由上表對比可見，SMA-GBDT 算法的擬合效果要優(yōu)于單一的GBDT 算法，并且SMA 優(yōu)化GBDT獲得的特征子集基本代表了整個(gè)數(shù)據(jù)集的特征，能夠很好地解釋變量。綜上可以看出，SMA-GBDT算法的綜合能力較強(qiáng)，在保證精度甚至增強(qiáng)擬合效果的同時(shí)大幅提高模型效率。

計(jì)算SMA-GBDT 模型的特征重要程度以確定表層土壤污染擴(kuò)散系數(shù)，影響因素及其重要性如表3所示。

表3 影響因素重要性

為清晰看出研究區(qū)域污染分布情況，從SMA-GBDT-3D CA 算法的模擬結(jié)果中隨機(jī)選取3條線上20 個(gè)點(diǎn)位的Hg 含量值與真實(shí)值的克里金插值及傳統(tǒng)元胞自動(dòng)機(jī)的模擬結(jié)果進(jìn)行對比，對比結(jié)果如圖6所示。

圖6 插值與模擬值對比

從圖中可以看出，研究算法的模擬結(jié)果相較于傳統(tǒng)元胞自動(dòng)機(jī)模型更加逼近克里金插值得到的數(shù)據(jù)曲線，模擬的精度更好。每條線可以反映污染擴(kuò)散后重金屬Hg 的分布，圖中三條曲線的變化趨勢基本相同，說明最終模擬的污染物分布情況與真實(shí)分布相符。

為了清晰展示研究區(qū)土壤重金屬污染擴(kuò)散過程，在模型模擬結(jié)果的基礎(chǔ)上，將三維可視化、WebGIS 與土壤重金屬點(diǎn)源污染擴(kuò)散過程結(jié)合起來，在Cesium 平臺中進(jìn)行污染物擴(kuò)散過程的三維動(dòng)態(tài)展示，如圖7所示。

圖7 污染擴(kuò)散三維可視化

將污染分布核密度圖疊加在地圖上，由深色到淺色表示污染物Hg 濃度由高到低的情況，污染核密度圖清晰地描述出污染位置、范圍，三維地形模型中還展示了研究區(qū)線、地形地勢等信息。受到風(fēng)向、位置等影響的點(diǎn)源污染擴(kuò)散過程用箭頭來表示，通過箭頭向四周的移動(dòng)表示污染擴(kuò)散的動(dòng)態(tài)過程，箭頭起始位置為污染源，結(jié)束位置為模擬出的污染范圍，可以看出由于主導(dǎo)風(fēng)向以及地理位置的影響，污染物擴(kuò)散路徑從八個(gè)方向產(chǎn)生偏移。擴(kuò)散過程動(dòng)態(tài)可視化清晰地展示出污染物從污染中心向四周逐漸擴(kuò)散，且距離污染源越遠(yuǎn)污染物濃度越低的擴(kuò)散規(guī)律。

5 結(jié)語

本文提出了邊界反彈機(jī)制的SMA-GBDT-3D CA污染擴(kuò)散模型，針對土壤重金屬污染問題，利用SMA 優(yōu)化的GBDT 算法獲得擴(kuò)散系數(shù)，建立三維元胞自動(dòng)機(jī)模型獲得擴(kuò)散模型。實(shí)驗(yàn)表明，提出的方法在保證模擬擴(kuò)散狀態(tài)及分布精度的同時(shí)增加了運(yùn)行效率，模型的綜合能力表現(xiàn)良好，因此該方法可以作為土壤重金屬點(diǎn)源污染擴(kuò)散三維模擬的有效方法。

土壤中重金屬的污染擴(kuò)散還會(huì)受到其他污染源、人為活動(dòng)、植被條件等因素的影響，在以后的研究中，還需要將這些影響因素考慮到模型中，得到更加準(zhǔn)確的擴(kuò)散規(guī)則來提高模型精度。