基于不同地質(zhì)統(tǒng)計方法的滲透系數(shù)場對污染物運移的影響

岳松梅，楊 蘊，吳劍鋒，吳吉春

(1.表生地球化學教育部重點實驗室，南京大學地球科學與工程學院，江蘇 南京 210023;2.河海大學地球科學與工程學院，江蘇 南京210098)

研究大尺度地下水流、污染物濃度和運動規(guī)律問題時，區(qū)域性模型構(gòu)建中的一個突出問題是勘探資料分布不均或缺乏，以及含水系統(tǒng)本身的非均質(zhì)性，造成了水文地質(zhì)參數(shù)的空間變異性。在通常的地下水數(shù)值模擬中只能用簡單的參數(shù)分區(qū)來描述參數(shù)的非均質(zhì)性，但由于分區(qū)過大而忽略分區(qū)內(nèi)部參數(shù)的差異性，往往導致模擬結(jié)果的不確定性。尋求一種盡可能利用有限的勘探資料，對未知區(qū)域內(nèi)的含水層參數(shù)進行合理估值的方法，是目前大區(qū)域地下水流模擬中的關(guān)鍵問題［1－2］。

在影響地下水系統(tǒng)的所有不確定因素中，滲透系數(shù)是表征含水層特性的一個關(guān)鍵要素［3－5］，它在空間上的分布變化及其復雜，其不確定性主要是由自然界含水介質(zhì)的非均質(zhì)各向異性和取樣、測試中的失真以及試驗量測誤差造成的。這種空間分布的不確定性使我們將其以隨機變量來處理，在不可能測得含水層每一點的介質(zhì)特征和水力特征的前提下，引入含水層空間隨機場，將實際含水層視為隨機場的一個實現(xiàn)。但實際情況中，滲透系數(shù)不僅具有隨機性，也具有一定的結(jié)構(gòu)性［6－9］。地質(zhì)統(tǒng)計學就是研究這種具有“二重性”的區(qū)域化變量的數(shù)學工具，其以變差函數(shù)為基本工具來研究分布于空間并呈現(xiàn)一定結(jié)構(gòu)性與隨機性的自然現(xiàn)象［10－11］，因而用地質(zhì)統(tǒng)計學方法來研究水文水資源水環(huán)境系統(tǒng)參變量的空間變異性，可以定量地揭示水文水資源水環(huán)境系統(tǒng)參變量在空間不同方向上的變化規(guī)律，辨識各參數(shù)的最大變異方向，尤其是它能給出區(qū)域化變量的最優(yōu)估計值及估計方差［10－11］，而這種特性是其他任何一種估計方法所不曾擁有的。

本文在前人研究的基礎上［12］，將地質(zhì)統(tǒng)計學推廣到大尺度的地下水污染運移的研究工作中，結(jié)合GSLIB軟件庫［18］，利用美國麻省軍事保留區(qū)(Massachusetts Military Reservation，MMR)的場地實測數(shù)據(jù)，采用普通克里格法和指示克里格法、順序高斯模擬法和順序指示模擬法等四種地質(zhì)統(tǒng)計方法，插值估測和模擬再現(xiàn)含水層滲透系數(shù)隨機場，進而對比研究四種滲透系數(shù)場對大尺度污染物運移的影響。

1 四種地質(zhì)統(tǒng)計方法簡述

1.1 普通克里格和指示克里格方法

普通克里格法屬線性平穩(wěn)地質(zhì)統(tǒng)計學范疇，是地質(zhì)統(tǒng)計學中最基本的方法，也是最常用的方法之一。設區(qū)域化變量Z(u)滿足二階平穩(wěn)假設，其數(shù)學期望為未知常數(shù)，協(xié)方差函數(shù)和變異函數(shù)存在且平穩(wěn)，以任一未知點為中心點的某一鄰域范圍內(nèi)n個樣本點，賦予相應權(quán)重計算未知點的屬性值，形成估計領域內(nèi)n個信息值的線性組合。通過求得一個無偏的最優(yōu)估計，并給出每個估值的誤差方差，該方法即為普通克里格方法(ordinary kriging，OK)［10－11］。

然而，在實際研究中還會碰到采樣數(shù)據(jù)中存在特異值得問題，所謂特異值是指那些比全部數(shù)值的均值或中位數(shù)高很多的值，其既非分析誤差所致，也非采樣方法等人為誤差引起，指示克里格法(indicator kriging，IK)就是為解決上述問題而發(fā)展起來的一種非參數(shù)地統(tǒng)計學法。該方法在給定條件數(shù)據(jù)的前下可以估計任意位置上的條件累積分布函數(shù)(conditional cumulative distribution function，CCDF)，適用于連續(xù)型和離散型的信息預測［10－11］。

相比于普通克里格，指示克里格法是將對區(qū)域化變量Z(u)的研究轉(zhuǎn)化為對其指示函數(shù)的研究，它無需假設數(shù)值來自某種特定分布的總體，也無需對原始數(shù)據(jù)進行變換，可以估計落于某些區(qū)間的每一部分的值，也可以估計比某一門坎值高(或低)的那部分值。

1.2 順序高斯模擬和順序指示模擬方法

順序高斯模擬(sequential Gaussian simulation，SGSIM)也稱序貫高斯模擬，是將順序模擬思想和高斯模擬相結(jié)合的一種隨機模擬方法。它假設一個高斯隨機域，基于變異函數(shù)特征，對每一個待模擬的空間位置利用克里格方法求構(gòu)建高斯分布函數(shù)，之后隨機在其中抽取一個數(shù)值作為該點的模擬值。

順序指示模擬(sequential indicator simulation，SISIM)將指示克里格方法和順序模擬算法相結(jié)合，實現(xiàn)了非參數(shù)的連續(xù)型與離散型的分布模擬。該方法不需要對原始條件數(shù)據(jù)分布的參數(shù)形式做任何假設，而是在現(xiàn)有資料的基礎上，通過一系列的門檻值把條件數(shù)據(jù)轉(zhuǎn)化成指示數(shù)據(jù)，相應的變差函數(shù)模擬轉(zhuǎn)換為指示變差函數(shù)模型，并采用指示克里格法對每個網(wǎng)格點處的條件概率分布進行估計［14－16］。

2 污染物運移及其結(jié)果評估

研究地下水污染問題必須用兩類方程來描述:第一類方程用來描述液體的流動;第二類方程用來描述溶質(zhì)即污染物質(zhì)在地下水中的運移。在一般或者較復雜的水文地質(zhì)條件下，污染物運移只有通過數(shù)值方法來實現(xiàn)，而且無論是采用有限單元或者有限差分方法，都要通過地下水流模擬地下水的流速場才能由運移模型得到污染物運移的模擬結(jié)果。在這里，分別采用常用的軟件水流模擬模型MODFLOW和溶質(zhì)運移模型MT3DMS來進行模擬。

2.1 滲透系數(shù)場的產(chǎn)生

2.1.1 原始數(shù)據(jù)的基本統(tǒng)計特征

包括均值、方差、標準差、變異系數(shù)、偏態(tài)系數(shù)和峰態(tài)系數(shù)的計算與分析;畫出原始數(shù)據(jù)分布直方圖并分析數(shù)據(jù)是否滿足正態(tài)分布特征。

根據(jù)原始數(shù)據(jù)分布直方圖，依據(jù)“偏度、峰度檢驗法”對原始數(shù)據(jù)進行正態(tài)分布檢驗，進而判斷原始數(shù)據(jù)的正態(tài)性。

2.1.2 原始數(shù)據(jù)的穩(wěn)健性處理

地質(zhì)統(tǒng)計學必須滿足數(shù)據(jù)為正態(tài)分布這一要求，如不滿足正態(tài)分布，直接計算的傳統(tǒng)地質(zhì)統(tǒng)計學變差函數(shù)必然有偏，進而影響估計的效果，故需對數(shù)據(jù)進行穩(wěn)健性處理。分別對數(shù)據(jù)先后進行去類分析、特異值處理和正態(tài)變換，使數(shù)據(jù)接近正態(tài)分布。

2.1.3 滲透系數(shù)分布場的產(chǎn)生

本文利用GSLIB軟件庫，對滲透系數(shù)進行原始數(shù)據(jù)基本參數(shù)統(tǒng)計和數(shù)據(jù)穩(wěn)健性處理，然后采用普通克里格法和指示克里格法、高斯序列模擬法和指示序列模擬法分別對數(shù)據(jù)進行插值和條件模擬，得到相應的滲透系數(shù)(lnK)分布場。

2.2 污染羽的空間矩評估

污染物在空間上的變化特征通?？梢酝ㄟ^刻畫污染羽隨時間變化的空間矩來評價。污染羽在空間上的零階矩表示污染物的總質(zhì)量;一階矩表示污染羽的質(zhì)心位置;二階矩表示污染羽圍繞質(zhì)心的分布范圍。

污染羽在任一時刻t關(guān)于原點的零階矩，即污染物總質(zhì)量，可以表示為

污染羽在空間上關(guān)于質(zhì)心的二階矩可以表示為協(xié)方差變量U:

本文以MMR為實例，建立水流模型和溶質(zhì)運移模型，利用空間矩評估污染羽的變化范圍。由于模型考慮化學反應的兩種基本類型:固液表面反應(吸附作用)和一階速率反應，污染物的質(zhì)量(零階矩)在運移過程中會發(fā)生改變，因此我們考慮污染羽的零階矩、一階矩和二階矩隨滲透系數(shù)場的空間變異特征。

3 實例研究

3.1 算例概況

美國麻省軍事保護區(qū)位于科德角(Cape Cod)的法爾茅斯鎮(zhèn)(Town of Falmouth)附近，始建于20世紀初，面積大約為89 km2。該區(qū)為第四紀冰川沉積物覆蓋，這些沉積物組成的非承壓含水層為科德角幾個社區(qū)的主要供水水源。位于MMR的東南部的美國空軍的給養(yǎng)基地，由于長期的軍事活動，出現(xiàn)了以三氯乙烯(Trichloroethylene，TCE)為主要污染物的化學物溢漏區(qū)(Chemical－Spill 10，CS－10)，對整個區(qū)域的環(huán)境造成嚴重破壞，直接影響了當?shù)鼐用裾５纳a(chǎn)生活。因此，本次研究就以CS－10區(qū)為模擬計算區(qū)域(圖1)。

圖1 研究區(qū)位置、CS－10模擬區(qū)域及TCE污染治理系統(tǒng)的平面分布［19］

3.2 模型的建立

本文建立了描述研究區(qū)含水層中TCE污染物運移的數(shù)學模型，采用有限差分法求解。其中只有含水層滲透系數(shù)為一隨機變量，其他模型參數(shù)均為常數(shù)。水流模型中的滲透系數(shù)K是通過地質(zhì)統(tǒng)計學對實測的值進行插值和條件模擬得到的，變化范圍較大，在泥沙含水層中可以小于3 m/d，粗砂含水層中能大于91 m/d。在假定研究區(qū)水流呈穩(wěn)定態(tài)，而溶質(zhì)運移呈非穩(wěn)定態(tài)的基礎上，整個研究區(qū)域剖分為159列、161行和21層，總模擬面積為57 km2，總模擬時間為30 a，1個應力期。水平方向步長為34m，越接近邊界步長會越大，垂直方向上每一層的厚度變化較大，最小可以小于1.5 m，最大可以大于15 m。側(cè)向邊界為定水頭邊界，由區(qū)域水流模型計算然后插值得到;上部為定流量邊界，補給量變化范圍為41～86 cm/a;下部隔水為零通量邊界。該區(qū)內(nèi)的地下水水流主要是沿水平方向從南向西南流動，水平的水力坡度平均值為0.001。

溶質(zhì)運移模型是建立在MT3DMS的基礎上的。該模型采用和MODFLOW水流模型相同的網(wǎng)格剖分，側(cè)向和垂向上都是定通量邊界。在邊界上，不考慮彌散項而只考慮對流項，如果流量為0，則溶質(zhì)運移也等于0。其余參數(shù)如，有效孔隙度為0.3，橫向彌散度、縱向彌散度和垂向彌散度分別為11 m、1.1 m 和 0.11 m［19］。TCE 濃度超過 5μgL－1的污染羽三維分布圖見圖2所示，

圖2 溶質(zhì)運移模型運行前CS－10 TCE污染羽濃度的三維分布圖

污染物在含水層中運移前初始污染物濃度均值為16 μg/L，最高濃度值遠超過 4 690 μg/L［13］。在美國，建立在保護人類健康和環(huán)境的安全飲用水法規(guī)的條件下，TCE的最大污染標準(MCL)是5 μg/L，這說明研究區(qū)域被高濃度的TCE污染。

3.3 滲透系數(shù)分布場的產(chǎn)生

在對研究區(qū)已知的375個觀測數(shù)據(jù)進行穩(wěn)健統(tǒng)計分析，并在此基礎上構(gòu)造該區(qū)域化變量的幾何各向異性套合模型，最后與GSLIB數(shù)據(jù)庫相結(jié)合，采用四種不同的地質(zhì)統(tǒng)計學方法分別對數(shù)據(jù)進行插值和條件模擬，得到相應的滲透系數(shù)(lnK)分布場。對于克里格插值來說，一旦各種參數(shù)確定(變差函數(shù)模型)以后，滲透系數(shù)(lnK)分布場就確定了，對應的污染物運移結(jié)果也確定了;而對于條件模擬來說，則每一次產(chǎn)生的滲透系數(shù)(lnK)分布場都是隨機的，相應的污染物運移結(jié)果也是隨機的。在這里，分別對克里格插值和條件模擬得出的結(jié)果做出對比，條件模擬模擬次數(shù)采用300次。

3.4 模擬結(jié)果與分析討論

普通克里格和指示克里格法所得的滲透系數(shù)(lnK)均值在3.7左右，范圍在－3～7之間，標準偏差在1左右，但指示克里格值(0.607 6)比普通格里格值(1.104 9)偏小。四種地質(zhì)統(tǒng)計方法得到的滲透系數(shù)場(lnK)估計方差的統(tǒng)計參數(shù)詳見表1。

表1 四種地質(zhì)統(tǒng)計方法得到滲透系數(shù)場(lnK)估計方差的統(tǒng)計參數(shù)

由表1可以得出，兩種克里格插值估計方差均值在1左右，指示克里格法略小于普通克里格法，但標準偏差比普通克里格法略大。圖3分別為污染運移5 475 d(15 a)和10 950 d(30 a)后，克里格插值下污染羽濃度的三維分布圖。

圖3 克里格插值生成滲透系數(shù)場所對應的15 a和30 a后污染羽分布圖(μgL－1)

從圖3中不難看出，普通克里格法的污染羽擴散更加均勻，指示克里格所對應的污染物濃度極值范圍更大。對于污染羽的二階空間矩(表2)評估，普通克里格插值對應的污染羽在空間上的展布范圍(二階矩)明顯大于和指示克里格，這主要受滲透系數(shù)的空間變異方差的影響，滲透系數(shù)方差的偏大，污染羽在空間上的二階矩隨之升高，表明污染物圍繞中心的擴散范圍不斷加大。這些不同之處是由于指示克里格法是用個別范圍個別變異函數(shù)來刻畫該區(qū)域化變量的關(guān)聯(lián)性的特性，無需假設數(shù)值來自某種特定分布(如正態(tài)分布)的總體，也無需對原始數(shù)據(jù)進行變換(如對數(shù)變換)，因此指示克里格法不必去掉重要而實際存在的高值數(shù)據(jù)即可處理各種不同現(xiàn)象［12］，而普通克里格法表現(xiàn)得更圓滑效應，滲透系數(shù)估計方差大且平滑，進而污染羽擴散更加均勻。

估計法是用加權(quán)平均法對研究區(qū)域的未知量求得線性、無偏和最小方差的內(nèi)插估計量及其相應的估計精度，而順序模擬法是通過系列隨機模擬對地質(zhì)變量進行局部估計，克服了估計法的平滑，能夠更好的再現(xiàn)真實曲線的波動性。對于條件模擬每一次產(chǎn)生的滲透系數(shù)(lnK)分布場都是隨機的，相應的污染物運移結(jié)果也是隨機的，這里采用300次模擬次數(shù)。

順序高斯模擬和順序指示模擬(模擬次數(shù)300次)所得的滲透系數(shù)(lnK)均值相同，即3.7997，范圍均在－3～7之間，但標準偏差順序指示模擬值(0.5941)比順序高斯模擬值(0.758 2)偏小。這兩種方法得到滲透系數(shù)場(lnK)估計方差均值在1左右，標準偏差順序高斯模擬值略小，詳見表1。

對比克里格插值得出的結(jié)果，四種地質(zhì)統(tǒng)計法所得到的滲透系數(shù)(lnK)均值在3.7左右，所對應污染羽擴散范圍基本一致，污染羽的質(zhì)心坐標也基本相同，驗證了質(zhì)心位置是由滲透系數(shù)的平均值來決定的理論。圖4分別為污染運移5 475 d(15 a)和10 950 d(30 a)后，條件模擬(模擬次數(shù)300次)下污染羽的濃度分布。表3為條件模擬模擬次數(shù)300次所對應的15 a和30 a后污染羽在空間上的分布特征。

綜合圖4和表3，不難看出，順序指示模擬的高濃度范圍明顯高于順序高斯模擬法，但污染羽的二階空間矩略小于后者。

表2 克里格插值方法所對應的15 a和30 a后污染羽在空間上的分布特征

圖4 條件模擬次數(shù)300次模擬生成滲透系數(shù)場所對應的15 a和30 a后污染羽分布圖

兩種模擬方法在模擬次數(shù)增加的情況下，污染羽的擴散范圍與質(zhì)心坐標沒有明顯差別，這有力的驗證了質(zhì)心位置是由滲透系數(shù)的平均值來決定的理論;但污染羽的空間二階矩在某一方向上卻有明顯減小的現(xiàn)象，這與二階空間矩主要受滲透系數(shù)的空間變異方差影響的結(jié)論相呼應。單獨對比順序高斯模擬與順序指示模擬這兩種方法，在模擬次數(shù)增加的情況下，滲透系數(shù)估計方差與污染物空間矩評估值這兩方面，后者的變化程度明顯高于前者。這是由于順序指示模擬采用指示克里金來估計，指示克里金不同于其他克里金法，它并不依賴于變量的平穩(wěn)性和要求變量服從某種分布。順序指示模擬法可以模擬復雜各向異性的地質(zhì)現(xiàn)象及連續(xù)性分布的極值，對于具有不同連續(xù)性分布的變量，可利用不同的變差函數(shù)進行表征，建立各向異性的模擬結(jié)果，進而更能體現(xiàn)數(shù)據(jù)的差異性。所以，采用順序指示模擬方法模擬的污染羽分布，其高濃度范圍明顯高于順序高斯模擬得出的結(jié)果，并且空間二階矩的變化程度明顯更大。

表3 條件模擬模擬次數(shù)300次所對應的15 a和30 a后污染羽在空間上的分布特征

4 結(jié)語

(1)本文利用野外場地實測數(shù)據(jù)，采用四種地質(zhì)統(tǒng)計方法，插值估測和模擬再現(xiàn)隨機滲透系數(shù)場，進而對比研究四種滲透系數(shù)場對大尺度污染物運移的影響。研究結(jié)果表明，污染羽一階空間矩證明了污染羽的質(zhì)心坐標基本相同，而四種地質(zhì)統(tǒng)計法所得到的滲透系數(shù)(lnK)范圍均在－3～7之間，均值在3.7左右，這有力的驗證了質(zhì)心位置是由滲透系數(shù)的平均值來決定的理論;通過指示克里格和順序指示模擬法得到的滲透系數(shù)(lnK)標準偏差比其余兩種更小，污染羽空間二階矩偏小，這也與二階矩主要受滲透系數(shù)的空間變異方差影響的結(jié)論相呼應。

(2)四種地質(zhì)統(tǒng)計方法在處理樣本數(shù)據(jù)時各有利有弊，應根據(jù)樣本的特點來作出相應對策。指示克里格與普通克里格得出的結(jié)果最大不同之處在于，前者污染羽的極值特點更加明顯，后者污染羽擴散更加均勻。條件模擬追求的是模擬的真實性，克服了估計法的平滑效果，能較好地再現(xiàn)真實曲線的波動性;隨著模擬次數(shù)的增大，滲透系數(shù)(lnK)估計方差與污染羽空間二階矩呈現(xiàn)變小的趨勢，并且順序指示模擬程度更加明顯，這是由于順序指示模擬在模擬復雜各向異性的地質(zhì)現(xiàn)象及連續(xù)性分布的極值上更加突出。

(3)需要說明的是，本文進行的污染物運移模擬是建立在四種地質(zhì)統(tǒng)計方法估計模擬滲透系數(shù)場的前提下實現(xiàn)的，不可避免得會出現(xiàn)估計誤差，除此之外，在變差函數(shù)理論模型擬合中存在必然的人為誤差和模型誤差。用數(shù)值法計算中也存在誤差，模型只能計算位于網(wǎng)格單元中心的結(jié)點處的濃度，而監(jiān)測孔的實際位置大多不是位于有限差分網(wǎng)格單元的中心。

［1］Kai－ Wei Juang，Dar－ Yuan Lee，Timothy R.Ellsworth.Using rank－order geostatistics for spatial interpolation of highly skewed data in a heavy － metal contaminated site［J］.Journal of Environmental Quality，2001，30:894－903.

［2］Bastante F G，Ordóez C，Taboada J，et al.Comparison of indicator kriging，conditional indicator simulation and multiple－point statistics used to model slate deposits［J］.Engineering Geology，2008，98:50－59.

［3］Diana dell＇Arciprete，Riccardo Bersezio，F(xiàn)abrizio Felletti，et al.Comparison of three geostatistical methods for hydrofacies simulation:a test on alluvial sediments［J］.Hydrogeology Journal，2012，20:299－311.

［4］Freeze R A.A stochastic conceptual analysis of one－dimensional groundwater flow in non － uniform homogeneous media［J］.Water Resources Research.1975，11(5):725－741.

［5］施小清，吳吉春，袁永生.滲透系數(shù)空間變異性研究［J］.水科學進展.2005，16(2):210－215.

［6］陳彥，吳吉春.含水層滲透系數(shù)空間變異性對地下水數(shù)值模擬的影響［J］.水科學進展.2005，16(4):482－48.

［7］劉文婷，朝倫巴根，劉艷偉，王力，宋君.地質(zhì)統(tǒng)計學方法在滲透系數(shù)空間變異性研究中的應用［J］.水利科技與經(jīng)濟.2010，16(4):364－366.

［8］李少龍，張家發(fā)，楊金忠.滲透系數(shù)空間相關(guān)性對滲流場統(tǒng)計特征影響分析［J］.長江科學院院報.2009，26(10):17－20.

［9］閆婷婷，吳劍鋒.滲透系數(shù)的空間變異性對污染物運移的影響研究［J］.水科學進展.2006，17(1):29－36.

［10］張仁鐸.空間變異理論及應用［M］.北京:科學出版社.2004.

［11］劉愛利，王培法，丁園園.地統(tǒng)計學概論［M］.北京:科學出版社.2012.

［12］劉玲玲，吳劍鋒，吳吉春.不同地質(zhì)統(tǒng)計方法在確定滲透系數(shù)場中的對比研究［J］.水文地質(zhì)工程地質(zhì).2009，5:66－71.

［13］吳劍鋒，彭偉，錢家忠.基于INPGA的地下水污染治理多目標優(yōu)化管理模型［J］.地質(zhì)論評.2011，57(3):437－443.

［14］李少華，張昌民，王振奇.利用順序指示模擬方法預測儲集層巖性［J］.新疆石油地質(zhì).2002，23(1):59－62.

［15］宋永忠，于晶，年靜波.隨機順序指示模擬技術(shù)與應用［J］.大慶石油地質(zhì)與開發(fā).2003，1:52－54.

［16］余振，何靜，魏福吉.序貫指示模擬和序貫高斯模擬在某地區(qū)精細流體預測中的聯(lián)合應用［J］.天然氣地球科學.2012，23(6):1170－1174.

［17］施小清，吳吉春，吳劍鋒.多個相關(guān)隨機參數(shù)的空間變異性對溶質(zhì)運移的影響［J］.水科學進展.2012，23(4):509－515.

［18］Deutsch C V，Journel A G.GSLIB－Geostatistcal Software Library and User’s Guide［M］.2nd ed.New York:Oxford University Press，1998.

［19］Zheng Chunmiao，Wang P P.2002.A field demonstration of the simulation－ optimization approach for remediation system design.Ground Water.40(3):258 ～265.