曹文翰，張 強(qiáng)，羅孝芹，時(shí) 曉，李紅葉

(1.成都理工大學(xué)國(guó)家重點(diǎn)實(shí)驗(yàn)室，成都 610059；2.中國(guó)電建集團(tuán)貴陽(yáng)勘測(cè)設(shè)計(jì)研究院有限公司，貴陽(yáng) 550081)

傳統(tǒng)的水文地質(zhì)調(diào)查方法需要對(duì)每戶人家的飲用水點(diǎn)進(jìn)行取樣，這種方法成本高，在調(diào)查點(diǎn)數(shù)目眾多且極為分散的情況下難度極大。地統(tǒng)計(jì)學(xué)可以通過少量樣本點(diǎn)的硝酸鹽濃度獲得大尺度下硝酸鹽濃度的分布情況[5]。其中克里金估值法相比于傳統(tǒng)空間的確定性插值具有較高的可靠性，不僅可以進(jìn)行插值，還可以定量計(jì)算誤差和不確定度。協(xié)同克里金法把區(qū)域化變量之間的最佳估值方法從單一屬性發(fā)展到了二個(gè)以上的協(xié)同區(qū)域化屬性，充分考慮變量之間的統(tǒng)計(jì)相關(guān)性和空間關(guān)系，提高了主變量的估計(jì)精度[6,7]。

協(xié)同克里金法要求協(xié)變量和主變量之間存在較強(qiáng)的聯(lián)系，變量之間聯(lián)系越強(qiáng)模型就越穩(wěn)定。因此，在眾多變量之中找出同硝酸鹽相關(guān)性最強(qiáng)的變量是進(jìn)行協(xié)同克里金插值的基礎(chǔ)。Distance模型可以計(jì)算觀測(cè)值之間的距離大小即變量之間的相似或不相似程度的測(cè)度，但無法給出顯著性P值；偏相關(guān)分析可以在多個(gè)變量存在時(shí)，控制其他多個(gè)變量的影響，計(jì)算兩個(gè)變量的相關(guān)性[8,9]。本文借助Particle模型和Distance模型判斷出相關(guān)性最強(qiáng)的變量并作為協(xié)變量，利用協(xié)同克里金方法反演研究區(qū)內(nèi)硝酸鹽的濃度，為后期的污染治理提供依據(jù)。

1 材料與方法

1.1 研究區(qū)域

研究區(qū)位于位于貴陽(yáng)市東郊龍洞堡地區(qū)，研究區(qū)內(nèi)居民分布松散，均是從山上引泉水下來生活飲用，飲用水源數(shù)目眾多且分散。在易采樣點(diǎn)地區(qū)收集24個(gè)水樣點(diǎn)進(jìn)行簡(jiǎn)分析后發(fā)現(xiàn)有超過42%的樣品中硝酸鹽含量超過美國(guó)環(huán)境保護(hù)署的最大污染物水平(EPA-MCL)指標(biāo)值10 mg/L。因此查明該地區(qū)的硝酸鹽分布情況顯得尤為重要(見圖1)。

圖1 地下水樣采集點(diǎn)分布圖Fig.1 The distribution of groundwater sampling points

1.2 相關(guān)性分析

偏相關(guān)分析在數(shù)理統(tǒng)計(jì)學(xué)中應(yīng)用廣泛，同傳統(tǒng)的雙變量分析相比，它能夠在多個(gè)相互影響的變量之間剔除多余變量的影響，只分析兩個(gè)變量之間的相關(guān)程度。該方法在處理高達(dá)12個(gè)變量的簡(jiǎn)分析數(shù)據(jù)時(shí)，擁有無可比擬的優(yōu)勢(shì)。

Distnce模型能將多個(gè)變量之間的相似程度用距離的大小表示出來，其結(jié)果在聚類分析、因子分析等多元分析方法之中應(yīng)用廣泛。在本文中選用歐幾里得距離公式，它能夠處理測(cè)量值相差懸殊的數(shù)據(jù)。

(1)

式中：i,j分別為數(shù)據(jù)的行列號(hào)；p為顯著性值。

1.3 地統(tǒng)計(jì)學(xué)方法

變異函數(shù)是地統(tǒng)計(jì)學(xué)所特有的基本工具，它既能表征區(qū)間變量的空間變異結(jié)構(gòu)，又能描述其隨機(jī)性變化，是很多統(tǒng)計(jì)學(xué)計(jì)算的基礎(chǔ)。樣本的變異函數(shù)可用下式計(jì)算：

(2)

式中：N(h)為樣點(diǎn)對(duì)的個(gè)數(shù)；Z為樣本硝酸鹽濃度；xi為樣本點(diǎn)的位置。

對(duì)于多個(gè)變量，常需要計(jì)算協(xié)變量函數(shù)γij(h)。

Zi(xi′+h)] [Zj(xi′)-Zj(xi′+h)]

(3)

式中：N′(h)為分離距離h的Zi(x)和Zj(x)的樣本對(duì)數(shù)。

如果變異函數(shù)表明變量具有空間相關(guān)性，克里金插值法就可以在變異函數(shù)理論及結(jié)構(gòu)分析基礎(chǔ)上，在有限區(qū)域內(nèi)對(duì)區(qū)域化變量的取值進(jìn)行無偏最優(yōu)估計(jì)并通過分析有限的樣本間的相關(guān)性，探索樣本間的分布關(guān)系，并進(jìn)行預(yù)測(cè)。協(xié)同克里金法把區(qū)域化變量的最佳估計(jì)方法從單一屬性發(fā)展到二個(gè)以上的協(xié)同區(qū)域化屬性，充分考慮變量之間的統(tǒng)計(jì)相關(guān)性和空間關(guān)系，提高了主變量的估計(jì)精度[10]。在二階平穩(wěn)假設(shè)條件下，其協(xié)方差公式為：

(4)

式中：Z1(xi)為主變量的測(cè)量值；Z2(xj)、λ1i、λ2j分別是Z1、Z2的權(quán)重，∑λ1j=1， ∑λ2j=1；m、n分別為主變量和協(xié)變量的測(cè)量點(diǎn)數(shù)量。

1.4 模型評(píng)價(jià)指標(biāo)

模型的評(píng)價(jià)指標(biāo)主要有標(biāo)準(zhǔn)平均值MS(Mean Standardized)以及均方根預(yù)測(cè)誤差RMSS(Root Mean Square Standardized)。如果預(yù)測(cè)誤差是無偏的，MS應(yīng)接近于0。RMSS是用來考量誤差的不確定性，以考察預(yù)測(cè)誤差是否是最優(yōu)的，RMSS越接近1精度越高，其公式如下：

(5)

2 結(jié)果和分析

2.1 數(shù)據(jù)分析

表1 Distance模型分析結(jié)果表Tab.1 Expressions of distance model

表2 偏相關(guān)分析結(jié)果表Tab.2 Expressions of particle model

2.2 半變異函數(shù)模型

克里金插值是在ArcGIS10.1的地統(tǒng)計(jì)分析模塊(Geostatistical Analyst,GA模塊)進(jìn)行的。同其他的插值軟件相比，ArcGIS最大的優(yōu)點(diǎn)在于它會(huì)自動(dòng)將一部分?jǐn)?shù)據(jù)作為訓(xùn)練樣本，另一部分?jǐn)?shù)據(jù)作為檢驗(yàn)樣本以在交叉驗(yàn)證中保證結(jié)果的精度[11]。

克里金插值要求數(shù)據(jù)本身或變換后滿足正態(tài)分布以及內(nèi)蘊(yùn)假設(shè)，因此在半變異函數(shù)建模之前，需對(duì)數(shù)據(jù)進(jìn)行探索性分析(見圖2)。發(fā)現(xiàn)數(shù)據(jù)進(jìn)行l(wèi)og變換后會(huì)更加符合正態(tài)分布且需剔除全局趨勢(shì)以排除變量的自相關(guān)性對(duì)模型精度產(chǎn)生的影響。

圖2 正態(tài)Q～Q plot分布圖Fig.2 Normal Q～Q plot of nitrate

地統(tǒng)計(jì)分析之后可以通過半變異函數(shù)建模來表示表面的非隨機(jī)即確定性的組成部分，即從數(shù)據(jù)中剔除趨勢(shì)后剩余的殘差繼續(xù)建模以表示表面局部變化的組成部分。塊金值C0與基臺(tái)值(C0+C)的比值表示可度量空間自相關(guān)的變異所占的比例，表明系統(tǒng)變量的空間自相關(guān)程度，比值越小，說明樣本間的變異受隨機(jī)變量影響越小，模型越穩(wěn)??；模型的變程反映了簡(jiǎn)分析指標(biāo)的空間連續(xù)性范圍，變程以外的范圍不存在相關(guān)性[12,13]。結(jié)果表明，協(xié)同克里金的塊金值C0與基臺(tái)值(C0+C)的比值為34.01%，優(yōu)于協(xié)同克里金的塊金值C0與基臺(tái)值(C0+C)的比值58.04%，但兩者均在25%～75%之間，表明系統(tǒng)具有中等的空間相關(guān)性，這與該地區(qū)地下水的補(bǔ)給來源較多且所取水樣成分復(fù)雜有關(guān)。兩個(gè)模型的變程均較大，這是由于該研究區(qū)域?qū)儆谕凰牡刭|(zhì)單元，空間相關(guān)的范圍較大，同時(shí)也與研究區(qū)的空間尺度較大有關(guān)(見表3)。

表3 半變異函數(shù)模型分析表Tab.3 Expression of semi-variogram model

2.3 預(yù)測(cè)結(jié)果的交叉論證

表4 交叉論證結(jié)果分析表Tab.4 Cross-validation results of model

2.4 研究區(qū)硝酸鹽濃度分析

從硝酸鹽濃度預(yù)測(cè)表面圖圖3可以看出，硝酸鹽濃度整體上西比東高，南北比中間高并在最北端達(dá)到最高值22.01 mg/L。在區(qū)域內(nèi)的敏感地帶趙家灣人口聚集區(qū)和黃泥哨水庫(kù)雖然硝酸鹽濃度相對(duì)較低，主要原因?yàn)樵摰貐^(qū)使用氮肥較少，但濃度仍超過10 mg/L，可能會(huì)對(duì)居民身體健康造成不利影響。

圖濃度預(yù)測(cè)表面圖Fig.3 The prediction surface map of Nitrate

3 結(jié) 語

(2)對(duì)原始數(shù)據(jù)進(jìn)行l(wèi)og變換并剔除趨勢(shì)面后，進(jìn)行半變異函數(shù)建模。模型的空間相關(guān)度均在25%～75%之間，表明系統(tǒng)具有中等的空間相關(guān)性，推測(cè)為地下水補(bǔ)給來源較多及所取水樣成分復(fù)雜的原因。相比較而言，協(xié)同克里金模型比普通克里金模型更加穩(wěn)健[C0/(C0+C)=34.01%，A0=1 435]。

(3)交叉論證結(jié)果表明，與普通克里金相比(MS=0.019 7,RMSS=0.979)，協(xié)同克里金的精度更高(MS=0.018 8,RMSS=0.988)，其所生成的硝酸鹽濃度預(yù)測(cè)表面圖能夠?yàn)槲廴疚锏闹卫硖峁┮欢ǖ闹笇?dǎo)意義。

參考文獻(xiàn)：

[1] Galloway J N, Aber J D, Erisman J W, et al. The nitrogen cascade [J]. Bioscience, 2003,53(4):341-356.

[2] Galloway J N, Townsend A R, Erisman J W, et al. Transformation of the nitrogen cycle: recent trends, questions and potential solutions[J].Science, 2008,320(5878):889-892.

[3] Vitousek P M, Naylor R, Grews T, et al. Nutrient imbalances in agricutural development [J]. Science, 2009,324(5934):1 519-1 520.

[4] Nie S, Gao Y. Review of current status and research approaches to nitrogen pollution in farmlands[J]. Agriculture Sciences in China, 2009,8(7):843-849.

[5] 王仁鐸，胡光道． 線性地質(zhì)統(tǒng)計(jì)學(xué)[M]．北京: 地質(zhì)出版社，1987．

[6] 劉治政,吳曉東,林洪孝. Kriging插值模型在地下水位監(jiān)測(cè)網(wǎng)優(yōu)化中的應(yīng)用[J]. 人民長(zhǎng)江,2010,41(9):14-17.

[7] 徐 馳,曾文治,黃介生,等. 基于高光譜與協(xié)同克里金的土壤耕作層含水率反演[J]. 農(nóng)業(yè)工程學(xué)報(bào),2014,30(13):94-103.

[8] 謝志英,劉 浩,唐新明. 北京市MODIS氣溶膠光學(xué)厚度與PM10質(zhì)量濃度的相關(guān)性分析[J]. 環(huán)境科學(xué)學(xué)報(bào),2015,35(10):3 292-3 299.

[9] 成玉婷,李 鵬,徐國(guó)策,等. 丹江流域氫氧同位素變化特征[J]. 水土保持學(xué)報(bào),2014,28(5):129-133.

[10] 李俊曉,李朝奎,殷智慧. 基于ArcGIS的克里金插值方法及其應(yīng)用[J]. 測(cè)繪通報(bào),2013,(9):87-90,97.

[11] Zhang Renduo. Applied geostatistics in environmental science[M]. Beijing: Science Press, 2005.

[12] Liu X, Wu J, Xu J. Characterizing the risk assessment of heavy metals and sampling uncertainty analysis in paddy field by geostatistics and GIS[J]. Environmental Pollution, 2006,141(2):257-264.

[13] Calzolari C, Ungaro F. Predicting shallow water table depth at regional scale from rainfall and soil data[J]. Journal of Hydrology, 2012,s414-415(2):374-387.