基于網格劃分的空間關聯區域VOCs濃度預測研究

陸秋琴， 蘭 瓊， 黃光球

(西安建筑科技大學 管理學院， 陜西 西安 710055)

隨著工業化的快速發展,區域性大氣污染日益突出，VOCs的大幅度排放引發了諸多環境問題。作為PM2.5和臭氧等大氣污染物的重要前體物，VOCs能發生光化學反應并生成有害的二次有機氣溶膠等物質[1-2]。部分VOCs易燃易爆，部分VOCs有毒，可以致癌、引起病變，嚴重危害人體健康[3]，所以“十三五”生態環境保護規劃將VOCs納入大氣污染防治的重要模塊[4]。因此，對VOCs濃度進行預測研究，有助于掌握其發展和變化規律，對制定有效的污染防治對策具有重要意義。不同的研究方法拓展和推動了預測理論的發展，為其他行業的預測研究提供了參考。同時，該預測研究可為環境保護規劃提供重要的數據積累，對開展污染控制有著積極的參考意義，也促進了公眾參與和居民環保意識的提高。

當前，對VOCs等大氣污染物濃度的預測研究主要是在其排放清單的基礎上展開的，通過建立基準年的污染物排放清單，來實現其他時段的預測[5]。國內外學者還利用大氣排放因子S型曲線預測大氣污染物的未來排放趨勢[6-7]；除此之外，還有基于情景分析法的污染物濃度預測，通過識別關鍵不確定因素，構建幾種可能出現的情景并分析內容[8]；優化模型也是污染物濃度預測的常見方法[9-10]。已經提出的大氣污染濃度預測模型主要有回歸分析、灰色模型[11]、神經網絡模型[12]、混沌模型[13]、基于時間序列的模型等[14]，以及他們的組合和改進模型。最優定權組合法大氣污染物濃度預測是基于多個空氣質量模式，以各單項空氣質量模式的組合預測誤差平方和最小為原則，構建出針對大氣污染的預測模型[15]。模糊綜合評價方法一般都是結合預測模型來使用。通過模糊聚類分析，將影響環境質量的各因素按主次區分，預測時考慮主要因素[16]。

以上研究還存在一些不足：①由于資金、地理條件等限制，對VOCs并不能做到全方位監測，所獲取的數據和信息不太完整；②研究主要集中在數量預測方面，較少通過劃分區域精細到每一個網格進行研究；③預測過程中較少考慮氣象指標等因素對預測結果的影響。為了解決上述問題，本文提出基于網格劃分的空間關聯區域VOCs濃度預測方法，以實現區域內VOCs精細化預測研究。

1 網格劃分與編號

1.1 區域坐標集合

根據選定區域建立相應的坐標系，建立原則為其中的每一點都能用坐標表示，可以取所選范圍比例尺為坐標刻度，獲取不同地方的坐標，形成區域坐標集合Rc：

Rc={(x1，y1)，(x2，y2)，…，(xn，yn)}

(1)

式中：(xi，yi)表示選定區域中的第i個坐標，用二維平面坐標表示，其中i=1，2，…，n；n表示區域坐標點總個數。

1.2 點云網格劃分算法與編號

點云網格劃分算法是利用點與點之間的距離關系來實現網格劃分，基于一點搜索臨近點形成線段，根據線段中點臨近檢索第三點，連接三點形成一個三角網格。對其新邊進行中點臨近檢索，依次形成網格體系，具體步驟如下。

1) 獲取區域坐標點集合Rc，初始化一個種子網格。基于點p1=(xm1,ym1)進行臨近檢索到第二個坐標點p2=(xm2,ym2)，連接兩點形成線段L(p1,p2),再基于線段L的中點臨近檢索第三點p3=(xm3,ym3)，連接點p3形成第一個三角網格，如圖1所示。將網格形成過程中產生的每條邊存入集合El，開始時El=?。

圖1 種子網格

El=El∪(p1,p2)∪(p4,p3)∪…∪(pi,pj)

i,j=1,2,…,n

(2)

2) 在種子網格的基礎上進行網格擴充，利用中點檢索，形成原始網格。從邊集合El中獲取未進行中點檢索的邊Lh(h=1,2,…,l；l為邊的數量)，其端點坐標為pi=(xmi,ymi)、pj=(xmj,ymj),計算其中點坐標Ci,j；從集合Rc檢索距離點Ci,j最近且未形成邊的點，中點邊與新點構造出兩條新邊，形成一個新的三角網格,并將新產生的邊存入集合El中。重復該步驟，直到邊集合El中不再提供外邊中點檢索為止。

(3)

(4)

3) 原始網格擴展，形成新網格。第二步結束形成一個原始網格，檢索集合Rc是否存在未形成邊的點，如果存在，則尋找新的種子網格重復第一、第二步，直到集合Rc不再有未形成邊的點為止，如圖2所示，此種情況下所選區域中存在大量的坐標點。在形成網格過程中，如果出現中斷現象，只需重復上述第一、第二步形成新的網格即可。

4) 編制網格順序碼，標識網格信息。在初始化種子網格時，將初始化的第一個三角網格編號為001，表示該區域的第一個網格。在網格擴充時，根據網格劃分步驟以及檢索點算法，對形成的新網格依次編號，最后輸出編號后的區域網格以及網格編號信息 [(pi，pj，pk)，Num](k=1,2,…,n)，如圖2所示。其中(pi，pj，pk)表示形成該網格的三個坐標點，即pi=(xmi,ymi)、pj=(xmj,ymj)、pk=(xmk,ymk)，Num表示網格編號，其編號值范圍為0～999的整數。

圖2 原始網格及網格編號圖

1.3 點云網格劃分的優點

1) 不規則劃分。根據所取點不規則形成大小不一的三角網格。

2) 劃分區域選點靈活。根據劃分需求可以隨意選取點，選點過程能有效避免山川、河流等地理條件的限制。

3) 自動編碼。在劃分過程中自動編碼表示網格，達到網格唯一性和明確性的要求。

4) 點利用率高。在網格劃分中采取三點為一的原因是可以將研究區域中所有的點全部劃分完，不會遺留未劃分的點。

2 空間關聯區域數據預估與收集

2.1 空間關聯區域數據預估

1) 網格數據預估原理

在實現VOCs精細化監管的過程中，將區域劃分成網格，在網格內設置監測點，監測設備在固定時段對網格內VOCs污染物進行監測，能夠準確地標識該網格內VOCs污染物的監測濃度值。但由于網格數眾多，并不是每一個網格都會設置監測點，為了收集和計算未設置監測點的網格數據，以及預估其污染物發展態勢，采取克里金插值法，通過已知網格數據及其與未知網格之間的空間關聯性來預估未知網格數據。

2) 克里金插值法預估過程

克里金插值被稱為空間最優無偏估計器，它是以變異函數理論和結構分析為基礎[17]，所選變異函數由數學期望、隨機場內特定點的數學期望、方差運算組成。克里金插值法會根據所選的變異函數模型進行模擬，最終對待估點進行預估。

設區域網格坐標點pi處設有監測點，監測值為V(pi)，i=1，2，…，n，則未設置監測點p0的估計值可以通過周圍n個監測點的監測值V(pi)求得，即

(5)

式中：λi為監測點pi的權重，λi的取值不僅要考慮監測點與預測點之間的距離，而且需結合二者的空間分布關系來確定,樣點分布如圖3所示。

圖3 樣點分布圖

設p0為待估計點，已知其鄰域內有p1，p2，…，p8共8個采樣點，其位置如圖3所示，各點的權重分別是λ1，λ2，…，λ8，由于圖中p1、p2、p3、p6到p0的距離相同，并且有p2與p3、p1與p6關于p0對稱，則有λ2=λ3，但由于樣點p5、p7、p8與p6叢聚在一起，這種叢聚作用降低了樣點p6對待估計點p0的影響，p1是一個單獨的樣點不存在叢聚影響，而且點p6與p0之間存在點p4，由于點p4距離點p0更近，對p6存在屏蔽效應，所以λ1>λ6。

要得到無偏最優估計值，必須滿足下面兩個條件：

a) 無偏估計，即E=[V(p0)-V*(p0)]=0

b) 估計方差最小，即

Var[V(p0)-V*(p0)]=min

則要求權重λi滿足下列方程：

(6)

2.2 數據收集及預處理

1) 數據收集

現有的監測設備不僅可以監測到VOCs的濃度(即單位體積排放量)，而且可以分析出該區域內VOCs不同組成成分的含量，并將監測數據上傳至服務器進行存儲，對于設有監測點的網格，通過監測設備獲取到VOCs監測值，并按照統一格式處理。已知監測點的監測數據，通過克里金插值法計算未設有監測點網格的VOCs組成成分預估值。將網格監測數據與網格預估數據合并，得到區域網格的VOCs污染物濃度值，如表1所示。

表1 區域網格VOCs污染物濃度值

表1對VOCs主要成分依次劃分了編號：苯為1號、甲苯為2號、……、苯乙烯為12號，并結合單元網格順序碼，描述不同網格中不同成分的監測濃度值，如V001(1)表示001號網格中苯的濃度值、V028(12)表示028號網格中苯乙烯的濃度值，依次收集得到區域網格VOCs污染物的濃度值。

2) 數據預處理

數據預處理是對收集到的網格數據進行整理的過程，通過研究區域每個網格的VOCs污染物濃度數據，形成區域VOCs污染物數據集合：

(7)

式中：D表示整個研究區域網格VOCs組成成分濃度集合；vij表示第i網格內第j類污染物的濃度值。

3 基于隨機森林算法的VOCs預測模型

3.1 數據集及模型結構

1) VOCs預測模型特征

VOCs預測模型特征可分為兩大類型，VOCs污染物和氣象指標，具體特征如表2所示。

表2中VOCs污染物特征是指VOCs污染物的烷類、烴類、酯類、醇類、苯系物等具體監測成分；氣象指標是指監測當天的氣象特征。表2中所有特征形成特征向量集合F。

表2 VOCs特征表

2) VOCs預測模型原始訓練樣本數據集合

基于研究區域VOCs污染物數據特征以及時間維度，形成區域VOCs數據集VD：

(8)

式中：αt1，αt2，…，αti和βt1，βt2，…，βti是時序特征向量，分別表示某一時間段內區域VOCs污染物濃度集合和區域VOCs總濃度序列數據；γ1，γ2，…，γm是非時序特征向量，包含氣象指標參數值、VOCs污染物特征。

3) VOCs預測模型構建

在上述數據處理的基礎上，運用隨機森林算法對研究區域VOCs濃度進行預測建模，建模過程如圖4所示。

圖4 基于隨機森林的區域VOCs預測建模過程

首先，利用Bootstrap方法從原始訓練樣本集VD中隨機抽取多個訓練樣本子集，對每個子集分別進行決策樹建模，然后利用測試集對各決策樹進行測試，綜合多棵決策樹測試結果，通過投票得出最終的預測模型。

3.2 訓練樣本子集的隨機選取

原始訓練樣本子集由兩部分構成：一類為VD中區域VOCs總量數據集合βti，將其作為預測模型的輸出；另一類為對應的區域網格VOCs污染物平均濃度集合αti和非時序特征數據集合γm，將其作為預測模型輸入。

利用Bootstrap方法從VD隨機選取w個訓練樣本子集V1，V2，…，Vw，用于構建w棵分類回歸樹(CART)。由于訓練樣本集的選取采用有放回的采樣方法，在采樣過程中會有36.8%的原始樣本不會出現在采集的樣本集合中，這些數據稱為袋數(out-of-bag，OOB)，對CART決策樹的誤差進行估計。對誤差估計取平均，便可得到隨機森林的泛化誤差估計值，由此可以對VOCs濃度預測模型的精度進行量化度量[18]。

3.3 CART決策樹的構建

對每個訓練樣本子集，采用CART算法生成一棵決策樹，共生成w棵決策樹。為保證決策樹構建的隨機性，采用隨機子空間思想，從VOCs特征集合F中隨機選取m個特征作為隨機特征變量，參與決策樹節點分裂過程，其中m≤log2(M+1)，而M表示特征集合F的集合長度。此外，整個隨機森林中決策樹的棵數w需根據預測結果來調整。

3.4 VOCs濃度預測結果投票及性能評價

1) VOCs濃度預測結果

當w棵樹構建完成后，利用測試集對數據進行仿真。將測試集數據Vk作為輸入，得到各決策樹模型預測的結果序列{fk1(V1)，fk2(V2)，…，fkw(Vw)}，基于隨機森林算法的預測模型最終預測輸出的VOCs濃度采用投票方式產生：

k=1,2,…,n

(9)

式中：Fk為組合預測模型；fki為單棵決策樹預測模型；I為示性函數；Yk為各決策樹預測的結果序列。將預測模型進行線性組合，即可得到區域VOCs濃度預測模型。

2) 性能評價指標

采用通用的模型誤差、擬合程度、效率作為度量指標，進行多模型量化評估，如平均相對誤差(MRE)和決定系數(R2)。其中R2表示模型輸入變量對輸出變量的解釋程度，也稱為擬合優度，取值在0到1之間。MRE越小，R2越接近于1，說明模型準確度越高。

(10)

(11)

4 案例分析

4.1 數據源

以西安市某區域涉及VOCs排放的企業為研究對象，企業清單來源于北極星網站，時間跨度為2018年6月至2018年12月。VOCs具體濃度數據通過企業年報、地方統計年鑒以及天氣后報網站獲得。將研究區域劃分成不同大小的網格，收集設有監測設備網格的污染物數據，通過克里金插值估計法計算出未設監測設備網格的污染物數據，形成VOCs數據集VD。

4.2 研究區域網格劃分及數據收集

1) 網格劃分

通過點云網格算法對西安市某區進行網格劃分并且對網格進行編號。首先獲取該區的坐標點集合，初始化種子網格，然后在種子網格的基礎上繼續擴充，形成新的網格，以此類推，將整個區域的網格劃分完畢，并編制網格順序碼，標識網格信息，結果如圖5所示。

圖5 西安市某區網格劃分及編號圖

2) 數據集

研究區域中有部分網格設有監測點，由監測點獲取到網格VOCs監測數據，包括VOCs污染物組分中的甲苯、乙烯、苯乙烯等12種物質，具體監測數值如表3所示。

表3 監測點VOCs污染物濃度值

4.3 空間關聯區域數據預估

根據網格坐標點及VOCs污染物濃度值，構建一個40×40的網格，標注范圍為1～40，即使網格間距為1。創建矩陣S和Y分別存儲坐標值和觀測值(即VOCs污染物濃度值)用于預測，根據其預估點和已知數值網格坐標點的空間位置，形成預測值表面，如圖6所示。

注：黑色點表示原始散點數據

根據圖6中預測值表面，結合每個點的擬合誤差值，求解出待估點的預估值，擬合誤差值如圖7所示。

圖7 擬合誤差值

在λi滿足式(6)的條件下，將其相關數值代入式(5)計算出未設有監測點網格的VOCs污染物預估值，具體數值如表4所示。

表4 預估點VOCs污染物濃度值

4.4 隨機森林模型預測

1) 模型構建及變量相關性分析

通過上述數據收集，獲得1 237組VOCs濃度數據，按式(8)處理得到數據集VD形成原始訓練樣本集，將其劃分為訓練集和驗證集，構建隨機森林回歸模型預測VOCs污染物濃度。VOCs特征集合F作為變量參與決策樹的分裂，模型預測中每個特征所起的作用不同，其相關系數如表5所示。

表5 VOCs部分特征相關系數表

根據相關系數表，VOCs與異丁烷以及環戊烷的線性相關性最大，相關系數達到了0.8以上，但是異戊烷與丙烯、甲苯之間的相關系數也達到了0.8以上，即各因素之間存在多重共線性，不滿足相互獨立條件，不能直接進行線性回歸，所以采用隨機森林預測。

2) 模型訓練、驗證和評估

將原始數據集合分為訓練集和驗證集，由式(10)、(11)分別進行模型的訓練和驗證，并對模型訓練和驗證結果進行評估，如表6所示。

表6 模型評估參數表

表6中訓練集和驗證集的相關評估參數值相差很小，其決定系數R2以及解釋度均達到了98%以上，表明模型在自變量不發生變化的情況下，因變量的變異概率極小。模型訓練過程中，各特征參數的重要性如圖8所示。

圖8 VOCs特征影響系數表

圖8表明， VOCs污染物的預測中，烷烴類污染物重要性比較強，相對而言溫度及壓強作用比較小。

3) VOCs污染物濃度預測

從設有監測點網格中選取19組數據作為預測集輸入模型，得到各決策樹的預測結果序列，再根據式(9)投票篩選出最優預測結果，預測結果如表7所示。

表7 VOCs污染物濃度預測結果

4) 模型比較

本文是基于網格空間特性以及隨機森林回歸模型實現VOCs污染物濃度預測，現將預測結果與常用的BP神經網絡預測結果進行比較，如表8所示。

表8給出了不同網格在兩種預測模型下的VOCs污染物預測值，未設置監測點的網格VOCs污染物實際值用克里金插值預測結果代替；分別采用相對誤差和平均相對誤差對兩種模型進行分析。由表8可知，隨機森林模型和BP神經網絡模型的VOCs總量預測值的平均誤差分別是3.15%和13.36%，由此可見，隨機森林回歸模型誤差更小。

表8 不同預測方法的結果對比

5 結 論

本次預測是根據區域空間關聯性以及VOCs污染物特征，對其濃度進行的精細化預測，意在解決監測設備不能普及部署以及區域之間污染物的流動影響問題。

1) 各區域之間的污染物存在相互影響。克里金插值法通過網格的空間地理位置來預估未設置監測點的網格數據；隨機森林模型基于污染物特征之間的相關關系預測污染物的濃度，隨機森林模型預測的結果更加精準。

2) 和BP神經網絡模型相比，隨機森林模型誤差更小，其VOCs總濃度預測值的平均誤差為3.15%。模型構建過程考慮了氣象指標對預測結果的影響，更能體現出VOCs特征之間的關聯性及相互影響作用。

3) 運用基于隨機森林算法的預測模型預測區域VOCs總濃度,同時也可以預測其組成成分的濃度(如苯、甲苯、苯乙烯等)，將其與國家VOCs排放控制標準限值進行對比，當超出限值時，結合區域網格編號信息[(pi，pj，pk)，Num]獲得其坐標信息(pi，pj，pk)，而坐標定位位置可為管理者超前管控提供依據。