普通克里金法與泰森多邊形法在土壤污染模擬中的應(yīng)用對比

周宜一 湯傳棟 陳國偉 崔志強(qiáng)

(銘瑞環(huán)境科技(南京)有限公司 江蘇南京 210000)

隨著我國對土壤環(huán)境的重視，近年來出臺(tái)了一系列土壤污染防治相關(guān)規(guī)定與技術(shù)規(guī)范，其中針對土壤污染狀況調(diào)查表明：污染物含量超過土壤污染風(fēng)險(xiǎn)管控標(biāo)準(zhǔn)的地塊，要求開展土壤污染風(fēng)險(xiǎn)評估。風(fēng)險(xiǎn)評估工作需在土壤污染狀況調(diào)查階段得到的分散采樣點(diǎn)檢測結(jié)果的基礎(chǔ)上，通過相關(guān)模型進(jìn)行統(tǒng)計(jì)分析，從而得到污染物空間分布，確定土壤污染范圍，估算修復(fù)面積與方量，為后續(xù)修復(fù)技術(shù)方案編制與土壤污染修復(fù)工程實(shí)施提供依據(jù)。

我國相關(guān)土壤技術(shù)導(dǎo)則尚未規(guī)定土壤污染范圍模擬方法，目前已有較多研究采用克里金插值法[1-6]、反距離加權(quán)法[7-8]、BP神經(jīng)網(wǎng)絡(luò)模型法[9-10]、無污染點(diǎn)位連線法等方法模擬污染范圍，近年來泰森多邊形法也出現(xiàn)了較多實(shí)際應(yīng)用。目前，多項(xiàng)研究表明克里金法較其他方法模擬精度更高，但前提是數(shù)據(jù)需服從或近似服從正態(tài)分布，否則對其模擬精度產(chǎn)生較大影響，目前尚無明確的結(jié)論表明何種模擬方法最優(yōu)[11-14]。

本次選取當(dāng)前常用的普通克里金插值法與泰森多邊形法模擬某地塊土壤污染范圍，研究模擬結(jié)果差異的影響因素，并從實(shí)際操作應(yīng)用角度對比分析兩種方法的優(yōu)缺點(diǎn)。

1 研究方法

1.1 研究地塊概況

本次研究選取的地塊位于江蘇省某地區(qū)，地塊總面積約為4.4萬m2，地塊歷史上存在碳素廠、金屬加工廠等多家工業(yè)企業(yè)，地塊周邊曾建有焦化廠等企業(yè)。由于該地塊未來用地性質(zhì)將轉(zhuǎn)換為中小學(xué)用地，按照《中華人民共和國土壤污染防治法》等相關(guān)法律法規(guī)要求，由土壤使用權(quán)人組織開展了地塊土壤污染狀況調(diào)查、風(fēng)險(xiǎn)評估等工作。

根據(jù)土壤污染狀況調(diào)查結(jié)論，地塊土壤污染物含量超過《土壤環(huán)境質(zhì)量建設(shè)用地土壤污染風(fēng)險(xiǎn)管控標(biāo)準(zhǔn)（試行）》（GB 36600-2018）規(guī)定的第一類建設(shè)用地土壤污染風(fēng)險(xiǎn)篩選值，屬于污染地塊。

地塊初步調(diào)查、詳細(xì)調(diào)查以及補(bǔ)充調(diào)查共布設(shè)了140個(gè)土壤采樣點(diǎn)，布點(diǎn)密度可滿足每400 m2不少于一個(gè)點(diǎn)位的要求，每個(gè)點(diǎn)位按照0.5 m 以內(nèi)取一個(gè)表層樣、0.5 m 以下每2 m 取一個(gè)樣的原則，分別取0.5 m、1.5 m、2.5 m、4 m、6 m的樣品進(jìn)行檢測分析。根據(jù)調(diào)查結(jié)果，地塊最大污染深度為1.5 m。本次選取地塊主要污染物苯并[a]芘作為模擬因子，基于地塊所有采樣點(diǎn)位不同深度樣品的檢測結(jié)果，將地塊分為0～＜0.5 m、0.5～1.5 m 兩層，分別采用普通克里金插值法與泰森多邊形法，模擬地塊污染范圍，并分析兩種方法的差異。研究地塊苯并[a]芘檢測結(jié)果基本情況如表1所示。

表1 研究地塊土壤中苯并[a]芘含量基本情況表

1.2 插值方法介紹

1.2.1 普通克里金插值法

普通克里金插值法是通過一組具有z值的分散點(diǎn)生成估計(jì)表面的統(tǒng)計(jì)過程。該方法基于點(diǎn)與點(diǎn)之間的統(tǒng)計(jì)關(guān)系進(jìn)行內(nèi)插或外推，利用區(qū)域化變量的原始數(shù)據(jù)和變異函數(shù)的結(jié)構(gòu)特點(diǎn)，對未知樣點(diǎn)進(jìn)行線性無偏、最優(yōu)估計(jì)。目前已有較多文獻(xiàn)介紹了利用Surfer 等軟件通過普通克里金法模擬污染范圍的操作方法與案例[15-17]。

本次研究將140個(gè)采樣點(diǎn)位0～＜0.5 m、0.5～1.5 m苯并[a]芘檢測值分別在Surfer軟件中進(jìn)行數(shù)據(jù)網(wǎng)格化，由于研究地塊數(shù)據(jù)不符合正態(tài)分布，因此對數(shù)據(jù)采用對數(shù)形式進(jìn)行變異函數(shù)理論模型轉(zhuǎn)換，轉(zhuǎn)換后0～＜0.5 m、0.5～1.5 m 數(shù)據(jù)峰度分別為1.26、3.19，數(shù)據(jù)偏度分別為1.29、1.94，較原始數(shù)據(jù)更接近正態(tài)分布，基于對數(shù)轉(zhuǎn)換后的數(shù)據(jù)建立線性白化網(wǎng)格，并生成等值線圖層，分別按篩選值設(shè)置最小等值線值，從而得到0～＜0.5 m、0.5～1.5 m污染范圍圖。

1.2.2 泰森多邊形法

泰森多邊形法構(gòu)造方法為在所有點(diǎn)中劃分出符合Delaunay準(zhǔn)則的不規(guī)則三角網(wǎng)（TIN），三角形各邊的垂直平分線即可形成泰森多邊形的邊，各平分線的交點(diǎn)決定泰森多邊形折點(diǎn)的位置。各個(gè)點(diǎn)覆蓋的區(qū)域中任何位置到其關(guān)聯(lián)點(diǎn)的距離都比到任何其他點(diǎn)輸入要素的距離近。

將140個(gè)采樣點(diǎn)位的苯并[a]芘按照0～＜0.5 m、0.5～1.5 m樣品檢測值分別列出并導(dǎo)入ArcGIS，創(chuàng)建泰森多邊形，并裁切得到地塊紅線范圍內(nèi)所有采樣點(diǎn)的泰森多邊形，利用按屬性選擇功能，輸入檢測值大于篩選值的條件表達(dá)式，可得到0～＜0.5 m、0.5～1.5 m深度污染范圍圖。

2 結(jié)果與分析

分別采用普通克里金法及泰森多邊形法插值得到0～＜0.5 m、0.5～1.5 m污染范圍圖，并將不同深度的各方法模擬結(jié)果疊加進(jìn)行對比，如圖1所示。

圖1 不同方法污染范圍插值結(jié)果疊加對比圖

根據(jù)模擬結(jié)果，可以看出泰森多邊形法插值結(jié)果不受點(diǎn)位檢測值影響，其污染范圍邊界固定位于相鄰點(diǎn)位連線的垂直平分線，而普通克里金法插值結(jié)果除布點(diǎn)位置外，還取決于點(diǎn)位污染物含量。選取圖1 中部分污染分區(qū)與相應(yīng)點(diǎn)位進(jìn)一步分析污染物含量對普通克里金法插值結(jié)果的影響，如表2 所示。由此可以看出：超標(biāo)點(diǎn)位檢測值越大，模擬得到污染范圍越大；插值結(jié)果也受到相鄰不超標(biāo)點(diǎn)檢測值影響，超標(biāo)點(diǎn)位檢測值不變的情況下，周邊不超標(biāo)點(diǎn)檢測值越大，污染范圍越大。

表2 污染物含量對普通克里金插值法污染范圍的影響對比

利用ArcGIS統(tǒng)計(jì)兩種方法插值得到的污染面積。泰森多邊形法插值得到的苯并[a]芘0～＜0.5 m深度污染面積為10 641 m2，0.5～1.5 m深度污染面積為7 599 m2；而普通克里金法對數(shù)轉(zhuǎn)換后插值得到的苯并[a]芘0～＜0.5 m 深度污染面積為8 454 m2，0.5～1.5 m 深度污染面積為6 333 m2。因此，基于該地塊苯并[a]芘調(diào)查數(shù)據(jù)，泰森多邊形法插值得到的污染面積大于普通克里金法插值得到的污染面積。

根據(jù)普通克里金法與泰森多邊形法插值原理，在對數(shù)轉(zhuǎn)換的情況下，研究地塊檢測值在5 mg/kg左右時(shí)兩種方法插值結(jié)果接近，如圖1 中G 點(diǎn)位。檢測值為4.2 mg/kg，污染面積分別為403.94 m2、394.9 m2，僅相差2%。分析研究地塊原始數(shù)據(jù)可知：0～＜0.5 m、0.5～1.5 m超標(biāo)點(diǎn)位檢測值中位數(shù)分別為1.45 mg/kg、1.80 mg/kg，上四分位數(shù)則分別為2.38 mg/kg、3.86 mg/kg，總體小于4 mg/kg，因此研究地塊普通克里金法插值得到的污染面積小于泰森多邊形法，如圖1 中A 點(diǎn)位，檢測值為0.6 mg/kg，普通克里金法與泰森多邊形法插值污染面積分別為392.61 m2、1.43 m2，顯示出明顯的差異。

3 討論

根據(jù)模擬過程操作及模擬結(jié)果對比普通克里金法與泰森多邊形法在土壤污染模擬中的應(yīng)用差異如下。

3.1 操作使用對比

根據(jù)相關(guān)研究，Surfer 軟件插值功能較ArcGIS 更優(yōu)[18-20]，且便于設(shè)置任意的最小等值線，從而得到污染范圍圖，但Surfer 中的普通克里金法每次僅能模擬單一深度單一因子污染范圍，不能進(jìn)行批量網(wǎng)格化與等值線設(shè)置，對于污染分層較多、污染因子較多的地塊，模擬工作量較大。同時(shí)，Surfer在坐標(biāo)提取與統(tǒng)計(jì)功能上相對欠缺，需配合CAD 或ArcGIS 等軟件進(jìn)一步分析，操作較為繁瑣。本研究模擬過程中，在Surfer中設(shè)置最大、最小等值線均為篩選值，得到污染范圍線，導(dǎo)出為.shp 文件，再導(dǎo)入ArcGIS 進(jìn)行定義投影、統(tǒng)計(jì)面積、提取坐標(biāo)等操作，由于研究地塊邊界處仍有污染，Surfer白化后使污染范圍線在地塊邊界處斷開，在導(dǎo)入ArcGIS后還需通過編輯器進(jìn)行閉合。若需得到地塊疊加污染范圍，或按重金屬污染、有機(jī)污染等不同污染類別進(jìn)行分類疊加，則需采取更多的操作步驟。

泰森多邊形法一般通過ArcGIS軟件進(jìn)行模擬，可通過按屬性要素選擇功能，輸入條件表達(dá)式，一次性得到所有污染單元，或根據(jù)需求改變條件表達(dá)式，分層、分因子模擬污染范圍，方法較靈活便捷。

3.2 模擬結(jié)果對比

泰森多邊形法僅通過點(diǎn)位之間的位置關(guān)系劃定單元，最終確定的污染范圍僅與布點(diǎn)密度有關(guān)，與點(diǎn)位污染程度無關(guān)，因此模擬結(jié)果不確定性較大。

普通克里金法在各土壤點(diǎn)位的位置關(guān)系基礎(chǔ)上，考慮各點(diǎn)位同一深度污染物含量，插值得到未知區(qū)域污染物含量，可直觀反映出地塊某一深度的污染物含量梯度變化情況。但在數(shù)據(jù)不能服從正態(tài)分布的情況下，模擬結(jié)果可能出現(xiàn)較大偏差，在超標(biāo)倍數(shù)較高的點(diǎn)位，因其平滑效應(yīng)使預(yù)測范圍過大，可能導(dǎo)致過度修復(fù)，或在超標(biāo)倍數(shù)較低的點(diǎn)位，可能預(yù)測范圍過小，將污染區(qū)誤判為清潔區(qū)。由于土壤的不均質(zhì)性，調(diào)查得到的污染物數(shù)據(jù)較多情況下不能服從正態(tài)分布，目前雖有較多文獻(xiàn)介紹利用Surfer 軟件普通克里金法進(jìn)行土壤污染模擬的操作方法，但基本直接采用線性插值，而數(shù)據(jù)通過變異函數(shù)進(jìn)行優(yōu)化后模擬結(jié)果可能發(fā)生較大變化。

3.3 劃定修復(fù)范圍應(yīng)用對比

普通克里金法模擬得到的污染范圍為平滑曲線，對于后續(xù)采取開挖方式進(jìn)行異位修復(fù)的情況，需通過裁彎取直進(jìn)一步劃定范圍，從而得到修復(fù)拐點(diǎn)坐標(biāo)，以便土壤開挖的實(shí)施，因此實(shí)際修復(fù)范圍將會(huì)擴(kuò)大，修復(fù)工程量相應(yīng)增加。且不同深度模擬得到的修復(fù)區(qū)邊界垂直投影差異較大，不便于分層施工定界。

利用泰森多邊形法模擬結(jié)果為直線多段線型修復(fù)邊界，可直接提取修復(fù)拐點(diǎn)坐標(biāo)，模擬得到的污染單元即為建議修復(fù)范圍。泰森多邊形法在各采樣點(diǎn)位不同深度均有數(shù)據(jù)的情況下，各層劃分的單元格垂直投影一致，便于實(shí)際施工。李靜文等人[21]采用克里金法模擬無數(shù)據(jù)點(diǎn)位，再利用泰森多邊形法進(jìn)行修復(fù)區(qū)劃定的方法，進(jìn)一步解決了泰森多邊形法中數(shù)據(jù)不足對施工定界的影響。

4 結(jié)語

由于土壤具有不均質(zhì)性，且不同污染源與污染因子在土壤中吸附與遷移過程復(fù)雜，目前雖有多種劃定修復(fù)范圍的方法，但尚無明確結(jié)論何種方法最優(yōu)。泰森多邊形法在操作使用與修復(fù)工程定界方面較普通克里金法有較大的優(yōu)勢，但插值結(jié)果不確定性較普通克里金法更大。

建議在劃定地塊污染范圍時(shí)，可采用兩種方法分別進(jìn)行模擬，對于模擬結(jié)果相近的區(qū)域，優(yōu)先選用泰森多邊形法，便于后續(xù)修復(fù)范圍劃定與修復(fù)施工；對于模擬結(jié)果差別較大的區(qū)域，結(jié)合地塊歷史用途、生產(chǎn)特征、土壤性質(zhì)等因素，具體分析數(shù)據(jù)特征、污染物遷移規(guī)律，并進(jìn)行加密布點(diǎn)和異常值分析，進(jìn)一步減少不確定性。