基于克里金的粵港澳大灣區(qū)典型化工廠地塊地下水流場(chǎng)分布空間分析

陸海建，莫健瑩，2，鄧一榮*，張曉露，吳 儉，劉麗麗

（1.廣東省環(huán)境科學(xué)研究院 廣東省污染場(chǎng)地環(huán)境管理與修復(fù)重點(diǎn)實(shí)驗(yàn)室／粵港澳環(huán)境質(zhì)量協(xié)同創(chuàng)新聯(lián)合實(shí)驗(yàn)室，廣東 廣州 510045；2.暨南大學(xué) 環(huán)境學(xué)院，廣東 廣州 510006）

隨著中國(guó)城市化進(jìn)程的加快與“退二進(jìn)三”政策的實(shí)施，大量有機(jī)化工等企業(yè)關(guān)閉、搬遷，遺留地塊污染問(wèn)題逐步凸顯，有些污染地塊甚至存在嚴(yán)重的地下水污染。據(jù)《全國(guó)地下水污染防治規(guī)劃（2011—2020）》數(shù)據(jù)，我國(guó)長(zhǎng)三角、珠三角、淮河流域等地區(qū)主要城市地下水中普遍存在有毒有機(jī)污染物檢出，部分地區(qū)檢出率超80%，而且越是城市化發(fā)展迅速、經(jīng)濟(jì)建設(shè)好的地區(qū)，其地下水中有毒物質(zhì)的種類(lèi)和數(shù)量也就越多［1-2］。

隨著地理信息系統(tǒng)的發(fā)展和對(duì)空間數(shù)據(jù)高質(zhì)量的要求，GIS 中的空間數(shù)據(jù)插值方法及其應(yīng)用也發(fā)揮著越來(lái)越重要的影響。在對(duì)地塊地下水流場(chǎng)研究時(shí)，通過(guò)有限的井位信息獲得整個(gè)地塊的地下水流場(chǎng)情況，這時(shí)候就可以利用ArcGIS 中的空間插值法來(lái)計(jì)算［3］。李俊曉等［4］運(yùn)用ArcGIS 的克里金插值方法對(duì)高程點(diǎn)數(shù)據(jù)進(jìn)行試驗(yàn)分析，給出克里金插值方法的精度分析及應(yīng)用優(yōu)勢(shì)；吳學(xué)文等［5］利用普通Kriging法對(duì)溫度的空間分布進(jìn)行模擬，給出參數(shù)的設(shè)置及變異函數(shù)模型選擇的合理性；楊陽(yáng)等［6］利用半變異函數(shù)及取樣間距對(duì)克里金法在海洋地層分析中的影響開(kāi)展了研究，通過(guò)對(duì)海底地層厚度進(jìn)行空間插值得到地層分布特征驗(yàn)證其插值效果。SHYAMALA 等［7］通過(guò)單變量統(tǒng)計(jì)分析，得到克里金法是插值曲面的最佳方法。LI 等［8］利用小波分析方法研究了西安地下水流場(chǎng)演變特征的地下水位時(shí)空變化。劉偉朋等［9］應(yīng)用ArcGIS 插值分析、柵格代數(shù)運(yùn)算、對(duì)比分析等方法，闡明了不同影響因素對(duì)地下水流場(chǎng)演化的控制作用。

空間插值分析算法是一種應(yīng)用于將離散點(diǎn)的測(cè)量數(shù)據(jù)轉(zhuǎn)換為連續(xù)數(shù)據(jù)表面的算法，能夠?qū)⑦B續(xù)數(shù)據(jù)曲面與其他空間現(xiàn)象的分布情況進(jìn)行比較，它在空間信息方面具有廣泛的應(yīng)用場(chǎng)景，尤其是地理信息方面［10］。其主要研究對(duì)象是人類(lèi)生活環(huán)境，其主要應(yīng)用場(chǎng)景包括遙感（RS）、全球定位系統(tǒng)（GPS）、地理信息系統(tǒng)（GIS）等技術(shù)［11］。ArcGIS 統(tǒng)計(jì)模塊中的空間插值法能根據(jù)某區(qū)域已獲得的部分空間信息來(lái)對(duì)未知空間屬性信息等進(jìn)行計(jì)算與分析，從而做到節(jié)約勘探成本。同時(shí)能較準(zhǔn)確的繪制出調(diào)查地塊地下水流場(chǎng)。文章就基于ArcGIS中克里金插值法的分析和應(yīng)用，利用克里金插值法來(lái)計(jì)算和分析粵港澳大灣區(qū)典型化工廠地塊的地下水流場(chǎng)，由此可分析地塊地下水流向的空間分布情況。

1 研究區(qū)域概況及數(shù)據(jù)來(lái)源

1.1 研究區(qū)概況

調(diào)查地塊為粵港澳大灣區(qū)一關(guān)停的化工廠，工廠于1966 年搬遷至本地塊，主要車(chē)間有乙酸乙酯車(chē)間、乙酸丁酯車(chē)間、甲醛車(chē)間、鍋爐房、機(jī)修車(chē)間等，面積為170 000 m2。2005年關(guān)停，廠區(qū)設(shè)備及廠房已全部關(guān)閉拆除。調(diào)查地塊位于兩山之間的山谷中，地貌為剝蝕殘丘與山間凹地地貌單元，地塊平均標(biāo)高為21.34 m。地塊南面和北面為山，北面調(diào)查最大標(biāo)高為29.40 m，南面調(diào)查最大標(biāo)高為26.90 m；場(chǎng)地中部地勢(shì)低于南面和北面，且呈中間高東西向低，地塊中部最大標(biāo)高為24.82 m，東面最低標(biāo)高為18.54 m，西面最低標(biāo)高為15.53 m。工廠生產(chǎn)運(yùn)營(yíng)時(shí)主要從事化工原料的生產(chǎn)加工。生產(chǎn)過(guò)程涉及的特征污染物主要有苯系物、多環(huán)芳烴、石油烴等有機(jī)物。

1.2 研究區(qū)水文地質(zhì)特征

調(diào)查地塊內(nèi)地下水按其賦存條件、水理性質(zhì)和水力特征，劃分為松散巖類(lèi)孔隙水和基巖裂隙水2 大類(lèi)型。松散巖類(lèi)孔隙水主要賦存于第四系沖殘積土層中，基巖裂隙水賦存于全風(fēng)化泥質(zhì)粉砂巖和基巖中。地塊內(nèi)揭露的地層有填土（雜填土、素填土）、粉質(zhì)粘土、砂質(zhì)粘土、砂土（主要是中砂）、全風(fēng)化泥質(zhì)粉砂巖和砂巖，其中填土層埋深厚度基本分布在0.5~4.0 m范圍，填土層土質(zhì)松散，滲透性較大，地塊部分區(qū)域地下水位分布在填土層中，該土層是本地塊一個(gè)富水性較好的含水層；粉質(zhì)粘土、砂質(zhì)粘土在地塊內(nèi)基本連續(xù)分布，埋深一般分布在填土層之下，滲透系數(shù)在2.3×10-6~4.4×10-5cm/s 范圍，滲透性較小，可概化為弱透水層；強(qiáng)風(fēng)化的泥質(zhì)粉砂巖，主要成份為砂質(zhì)粘土局部含有粉質(zhì)粘土，土層較密實(shí)，滲透性大小與上述粉質(zhì)粘土差不多，可概化為弱透水層。調(diào)查地塊監(jiān)測(cè)井測(cè)得地下水位埋深在1.21~5.55 m，地下水位標(biāo)高在11.81~20.25 m。

1.3 數(shù)據(jù)來(lái)源

文章用于地下水流場(chǎng)插值數(shù)據(jù)主要是通過(guò)現(xiàn)場(chǎng)實(shí)際測(cè)試獲取，首先通過(guò)RTK測(cè)量獲取每個(gè)地下水監(jiān)測(cè)井孔口的空間信息（X、Y、Z），然后利用地下水水位測(cè)試儀測(cè)量每個(gè)孔的地下水埋深，根據(jù)上述獲取的基礎(chǔ)數(shù)據(jù)換算每個(gè)孔地下水潛水面的空間信息（Xi、Yi、Zi），利用上述獲取的離散測(cè)量點(diǎn)位數(shù)據(jù)導(dǎo)入到ArcGIS軟件中［3-4］，文章獲取了地塊內(nèi)48個(gè)水位的空間信息，在地質(zhì)統(tǒng)計(jì)學(xué)中屬于小樣本案例［12］，各點(diǎn)位的分布情況如圖1。

圖1 研究區(qū)地下水監(jiān)測(cè)井分布圖

2 普通克里金法與變異函數(shù)理論基礎(chǔ)

2.1 普通克里金法

普通克里金法是以區(qū)域化變量理論和空間連續(xù)性理論為基礎(chǔ)，以變異函數(shù)為基本工具，研究在空間分布上既有結(jié)構(gòu)性又有隨機(jī)性的地質(zhì)變量，建立符合地質(zhì)規(guī)律的統(tǒng)計(jì)模型來(lái)反映變量的變化規(guī)律，并對(duì)其空間分布進(jìn)行預(yù)測(cè)［13］。普通克里金插值法對(duì)于任意待估點(diǎn)的估計(jì)值是通過(guò)該點(diǎn)影響范圍內(nèi)的n個(gè)有效樣品值的線性組合得到，表達(dá)式如下所示。

其中：Z（x0）為區(qū)域內(nèi)x0位置的預(yù)測(cè)值；Z（xi）為區(qū)域內(nèi)xi位置的觀測(cè)值；λi為克里金權(quán)重系數(shù)；n為監(jiān)測(cè)點(diǎn)個(gè)數(shù)。

2.2 變異函數(shù)理論基礎(chǔ)

由于在普通克里金插值法中λ值是通過(guò)樣本變異函數(shù)求得的，通過(guò)變異函數(shù)的計(jì)算結(jié)果，能夠反映區(qū)域化變量的結(jié)構(gòu)性［14］，因此，選擇一個(gè)合適的變異函數(shù)來(lái)擬合該地塊地下水流場(chǎng)的變異性，有助于提高克里金插值的效果。

在本次調(diào)查地塊地下水流場(chǎng)空間插值研究中，選擇ArcGIS中普通克里金法進(jìn)行空間插值研究，并對(duì)普通克里金法插值中的球面、指數(shù)、高斯變異函數(shù)模型進(jìn)行進(jìn)一步的優(yōu)選，優(yōu)選出最佳的模型繪制出最準(zhǔn)確的地下水流場(chǎng)圖。其中球面函數(shù)，指數(shù)函數(shù)、高斯的變異函數(shù)模型表達(dá)式如表1。

表1 變異函數(shù)模型表達(dá)式

2.3 精度評(píng)價(jià)

為進(jìn)一步優(yōu)選出最佳的變異函數(shù)模型對(duì)插值的結(jié)果進(jìn)行交叉驗(yàn)證分析。交叉驗(yàn)證法能夠評(píng)價(jià)和比較不同插值方法，是插值精度評(píng)價(jià)的重要方法［15-16］。交叉驗(yàn)證基本思路是對(duì)每一個(gè)有觀察值的地方，去除原有的觀察值運(yùn)用其他點(diǎn)觀察值對(duì)其進(jìn)行估算，如此循環(huán)對(duì)所有的觀察值進(jìn)行逐點(diǎn)估算。然后用原有的觀察值與該點(diǎn)的估算值進(jìn)行對(duì)比，從而比較選出最佳的變異函數(shù)模型［5-6］。

采用平均值誤差（ME）、均方根誤差（RMSE）、標(biāo)準(zhǔn)平均值誤差（SEM）、標(biāo)準(zhǔn)均方根誤差（RMSSE）和平均標(biāo)準(zhǔn)誤差（ASE）來(lái)評(píng)價(jià)插值結(jié)果的精度。對(duì)于1 個(gè)精確的模型，其平均誤差和標(biāo)準(zhǔn)平均誤差的絕對(duì)值應(yīng)接近于0，均方根誤差和平均標(biāo)準(zhǔn)差誤差最接近，其標(biāo)準(zhǔn)均方根誤差應(yīng)該接近于1［17-18］。

3 研究數(shù)據(jù)的空間分布分析

3.1 數(shù)據(jù)分布狀態(tài)

本研究利用SPSS 軟件對(duì)所得數(shù)據(jù)進(jìn)行正態(tài)分布和對(duì)數(shù)正態(tài)分布檢驗(yàn)，分析結(jié)果如圖2、表2、表3。檢驗(yàn)結(jié)果表明，該地區(qū)地下水水位高程H的范圍在11.81~20.25 m，經(jīng)自然對(duì)數(shù)變換后ln H范圍為2.47~3.01 m。由圖2 對(duì)比結(jié)果顯示，地下水水位高程H 和ln H 的標(biāo)準(zhǔn)Q-Q分析圖2散點(diǎn)分布近似一條直線并沿對(duì)角線分布，數(shù)據(jù)符合正態(tài)分布的要求；根據(jù)表3，2組數(shù)據(jù)的偏度與峰度約為0，同時(shí)偏度與峰度的Z-score 絕對(duì)值均小于1.96，近似服從正態(tài)分布；本研究有效樣本數(shù)為48，采用S-W 檢驗(yàn)的P值分別為0.70 和0.25，均大于0.05，通過(guò)5%的置信區(qū)間檢驗(yàn)，說(shuō)明該地塊地下水水位高程數(shù)據(jù)同時(shí)服從正態(tài)分布和對(duì)數(shù)正態(tài)分布，有利于提高插值預(yù)測(cè)精度［19］。由于S-W檢驗(yàn)結(jié)果中，正態(tài)分布結(jié)果優(yōu)于對(duì)數(shù)正態(tài)分布，所以將地下水水位高程H作為研究對(duì)象進(jìn)行空間插值研究。

圖2 H和lnH 的標(biāo)準(zhǔn)Q-Q分析圖

表2 地下水水位高程數(shù)據(jù)統(tǒng)計(jì)表

表3 地下水水位高程數(shù)據(jù)正態(tài)分布分析表

3.2 變異函數(shù)空間比較分析

根據(jù)球面、指數(shù)和高斯變異函數(shù)模型的理論公式計(jì)算，所得變異函數(shù)理論模型擬合參數(shù)結(jié)果如表4。

表4 變異函數(shù)模型參數(shù)表

模型擬合曲線如圖3。3 種變異函數(shù)的變程均大于步長(zhǎng)間隔，表明在步長(zhǎng)間隔范圍內(nèi)3 種變異函數(shù)均存在空間自相關(guān)性［13］；從結(jié)構(gòu)性角度來(lái)看，根據(jù)區(qū)域化變量空間相關(guān)性程度分級(jí)標(biāo)準(zhǔn)［5］，塊金系數(shù)（C0/C0+C）能表明區(qū)域化變量之間的空間相關(guān)程度：塊金系數(shù)＜25%表明區(qū)域化變量之間的空間相關(guān)性強(qiáng)烈；塊金系數(shù)在25%~75%范圍內(nèi)表明區(qū)域化變量之間的空間相關(guān)性中等；塊金系數(shù)＞75%表明區(qū)域化變量之間的空間相關(guān)性很弱。3種變異函數(shù)的塊金系數(shù)均小于25%，說(shuō)明該地塊地下水水位高程H具有強(qiáng)烈的空間相關(guān)性。

圖3 變異函數(shù)模型擬合曲線圖

3.3 交叉驗(yàn)證分析

本次對(duì)調(diào)查地塊地下水流場(chǎng)空間插值研究中，交叉驗(yàn)證預(yù)計(jì)誤差中的指標(biāo)參數(shù)統(tǒng)計(jì)如表5。

表5 各變異函數(shù)模型交叉驗(yàn)證參數(shù)對(duì)比統(tǒng)計(jì)表

其中，球面數(shù)模型的均方根預(yù)測(cè)誤差最小，標(biāo)準(zhǔn)平均值誤差最接近于0，標(biāo)準(zhǔn)均方根預(yù)測(cè)誤差較接近于1，平均標(biāo)準(zhǔn)誤差最接近于均方根預(yù)測(cè)誤差。所以，球面數(shù)模型交叉驗(yàn)證參數(shù)基本滿足最優(yōu)標(biāo)準(zhǔn)的數(shù)值，插值擬合效果可信度最高、效果最佳。球面函數(shù)模模型擬合效果如圖4。

圖4 球面函數(shù)模型擬合參數(shù)圖

4 地下水流場(chǎng)插值效果討論與分析

本次調(diào)查地塊南北面均為山體且地塊內(nèi)局部存在一地勢(shì)較高的區(qū)域，地塊東南面約1 km 處有1 條地表河流。山體為天然分水嶺，東側(cè)的河流面做為地下水排泄面。根據(jù)調(diào)查地塊所在區(qū)域的地形地貌、地表水系的分布以及水文地質(zhì)特征，可以判定地塊內(nèi)的地下水流向整體為從西向東流。但是在廠區(qū)地塊內(nèi)的南北軸方向上存在一局部的分水嶺，因此地塊內(nèi)分水嶺東側(cè)地下水整體向東流，分水嶺西側(cè)地下水整體向西南方向流。3種變異函數(shù)模型插值的地下水等水位線表明調(diào)查地塊中部的北側(cè)區(qū)域存在地下水水位高值區(qū)域，結(jié)果顯示球面函數(shù)模型與指數(shù)函數(shù)模型空間插值與局部分水嶺的分布吻合，而高斯函數(shù)模型空間插值的效果吻合度較弱。

為進(jìn)一步核實(shí)驗(yàn)證模型插值的準(zhǔn)確性，利用地下水流速流向儀對(duì)研究區(qū)GW33點(diǎn)位的流速和流向進(jìn)行了現(xiàn)場(chǎng)測(cè)定，如圖5，GW33點(diǎn)位地下水流向?yàn)槲鞅毕驏|南流（流向角度為150°）、測(cè)定的平均流速約為0.45 m/d。球面函數(shù)模型與指數(shù)函數(shù)模型空間插值的結(jié)果顯示GW33 點(diǎn)處地下水流向?yàn)槲鞅毕驏|南流（流向角度約為150°）與實(shí)測(cè)流向高度吻合，而高斯函數(shù)模型空間插值的效果與實(shí)測(cè)值吻合度較弱。

圖5 研究區(qū)GW33點(diǎn)位地下水流速流向?qū)崪y(cè)值

根據(jù)3種半變異函數(shù)模型插值擬合后的地下水流場(chǎng)的對(duì)比分析如圖6，球面函數(shù)模型與指數(shù)函數(shù)模型空間插值的結(jié)果與實(shí)際結(jié)果吻合度較高，通過(guò)2 種函數(shù)擬合的交叉驗(yàn)證的參數(shù)對(duì)比分析球面函數(shù)模型插值結(jié)果可信度更高。

圖6 不同模型的空間插值效果圖

高斯函數(shù)模型空間插值的效果與實(shí)際結(jié)果吻合度較弱，同時(shí)交叉驗(yàn)證參數(shù)結(jié)果也顯示其擬合效果較差［20-22］。選用球面函數(shù)作為變異函數(shù)，利用地質(zhì)統(tǒng)計(jì)軟件GS+擬合地下水流場(chǎng)空間模型如圖7。擬合模型結(jié)果的流場(chǎng)方向與地下水實(shí)際分布結(jié)果吻合，可知球面函數(shù)模型空間插值的流場(chǎng)結(jié)果更貼合實(shí)際情況。

圖7 地下水流場(chǎng)克里金（球面函數(shù)）插值空間模型

5 結(jié)論

（1）文章以化工廠地塊內(nèi)48個(gè)地下水水位的空間信息的克里金法插值為例，經(jīng)分析其正態(tài)QQ 圖考察數(shù)據(jù)集符合正態(tài)分布，克里金插值法適用于本研究。研究區(qū)地塊內(nèi)分水嶺東側(cè)地下水整體向東流，分水嶺西側(cè)地下水整體向西南方向流。

（2）通過(guò)對(duì)球面、指數(shù)和高斯變異函數(shù)模型的有關(guān)參數(shù)的比較，采用交叉驗(yàn)證的方法，利用各精度分析指標(biāo)，比選出球面變異函數(shù)模型在本次插值應(yīng)用中數(shù)據(jù)的空間自相關(guān)程度最高，擬合的效果最優(yōu)。

（3）根據(jù)實(shí)測(cè)流場(chǎng)及地形判定流場(chǎng)的結(jié)果與普通克里金球面模型擬合的地下水流場(chǎng)空間分布的對(duì)比驗(yàn)證分析，模型插值擬合結(jié)果與實(shí)測(cè)及判定分析結(jié)果較吻合。