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

（1.廣東省環(huán)境科學(xué)研究院 廣東省污染場地環(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）

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

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

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

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

1.1 研究區(qū)概況

調(diào)查地塊為粵港澳大灣區(qū)一關(guān)停的化工廠，工廠于1966 年搬遷至本地塊，主要車間有乙酸乙酯車間、乙酸丁酯車間、甲醛車間、鍋爐房、機(jī)修車間等，面積為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；場地中部地勢低于南面和北面，且呈中間高東西向低，地塊中部最大標(biāo)高為24.82 m，東面最低標(biāo)高為18.54 m，西面最低標(biāo)高為15.53 m。工廠生產(chǎn)運(yùn)營時主要從事化工原料的生產(chǎn)加工。生產(chǎn)過程涉及的特征污染物主要有苯系物、多環(huán)芳烴、石油烴等有機(jī)物。

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

調(diào)查地塊內(nèi)地下水按其賦存條件、水理性質(zhì)和水力特征，劃分為松散巖類孔隙水和基巖裂隙水2 大類型。松散巖類孔隙水主要賦存于第四系沖殘積土層中，基巖裂隙水賦存于全風(fēng)化泥質(zhì)粉砂巖和基巖中。地塊內(nèi)揭露的地層有填土（雜填土、素填土）、粉質(zhì)粘土、砂質(zhì)粘土、砂土（主要是中砂）、全風(fēng)化泥質(zhì)粉砂巖和砂巖，其中填土層埋深厚度基本分布在0.5~4.0 m范圍，填土層土質(zhì)松散，滲透性較大，地塊部分區(qū)域地下水位分布在填土層中，該土層是本地塊一個富水性較好的含水層；粉質(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)測井測得地下水位埋深在1.21~5.55 m，地下水位標(biāo)高在11.81~20.25 m。

1.3 數(shù)據(jù)來源

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

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

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

2.1 普通克里金法

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

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

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

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

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

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

2.3 精度評價

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

采用平均值誤差（ME）、均方根誤差（RMSE）、標(biāo)準(zhǔn)平均值誤差（SEM）、標(biāo)準(zhǔn)均方根誤差（RMSSE）和平均標(biāo)準(zhǔn)誤差（ASE）來評價插值結(jié)果的精度。對于1 個精確的模型，其平均誤差和標(biāo)準(zhǔn)平均誤差的絕對值應(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 軟件對所得數(shù)據(jù)進(jìn)行正態(tài)分布和對數(shù)正態(tài)分布檢驗(yàn)，分析結(jié)果如圖2、表2、表3。檢驗(yàn)結(jié)果表明，該地區(qū)地下水水位高程H的范圍在11.81~20.25 m，經(jīng)自然對數(shù)變換后ln H范圍為2.47~3.01 m。由圖2 對比結(jié)果顯示，地下水水位高程H 和ln H 的標(biāo)準(zhǔn)Q-Q分析圖2散點(diǎn)分布近似一條直線并沿對角線分布，數(shù)據(jù)符合正態(tài)分布的要求；根據(jù)表3，2組數(shù)據(jù)的偏度與峰度約為0，同時偏度與峰度的Z-score 絕對值均小于1.96，近似服從正態(tài)分布；本研究有效樣本數(shù)為48，采用S-W 檢驗(yàn)的P值分別為0.70 和0.25，均大于0.05，通過5%的置信區(qū)間檢驗(yàn)，說明該地塊地下水水位高程數(shù)據(jù)同時服從正態(tài)分布和對數(shù)正態(tài)分布，有利于提高插值預(yù)測精度［19］。由于S-W檢驗(yàn)結(jié)果中，正態(tài)分布結(jié)果優(yōu)于對數(shù)正態(tài)分布，所以將地下水水位高程H作為研究對象進(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ù)的變程均大于步長間隔，表明在步長間隔范圍內(nèi)3 種變異函數(shù)均存在空間自相關(guān)性［13］；從結(jié)構(gòu)性角度來看，根據(jù)區(qū)域化變量空間相關(guān)性程度分級標(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%，說明該地塊地下水水位高程H具有強(qiáng)烈的空間相關(guān)性。

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

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

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

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

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

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

4 地下水流場插值效果討論與分析

本次調(diào)查地塊南北面均為山體且地塊內(nèi)局部存在一地勢較高的區(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)確性，利用地下水流速流向儀對研究區(qū)GW33點(diǎn)位的流速和流向進(jìn)行了現(xiàn)場測定，如圖5，GW33點(diǎn)位地下水流向?yàn)槲鞅毕驏|南流（流向角度為150°）、測定的平均流速約為0.45 m/d。球面函數(shù)模型與指數(shù)函數(shù)模型空間插值的結(jié)果顯示GW33 點(diǎn)處地下水流向?yàn)槲鞅毕驏|南流（流向角度約為150°）與實(shí)測流向高度吻合，而高斯函數(shù)模型空間插值的效果與實(shí)測值吻合度較弱。

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

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

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

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

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

5 結(jié)論

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

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

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