基于協(xié)同克里金的飲用水硝酸鹽濃度預測研究

曹文翰，張 強，羅孝芹，時 曉，李紅葉

(1.成都理工大學國家重點實驗室，成都 610059；2.中國電建集團貴陽勘測設計研究院有限公司，貴陽 550081)

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

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

1 材料與方法

1.1 研究區(qū)域

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

圖1 地下水樣采集點分布圖Fig.1 The distribution of groundwater sampling points

1.2 相關性分析

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

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

(1)

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

1.3 地統(tǒng)計學方法

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

(2)

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

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

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

(3)

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

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

(4)

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

1.4 模型評價指標

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

(5)

2 結果和分析

2.1 數(shù)據分析

表1 Distance模型分析結果表Tab.1 Expressions of distance model

表2 偏相關分析結果表Tab.2 Expressions of particle model

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

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

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

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

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

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

2.3 預測結果的交叉論證

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

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

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

圖濃度預測表面圖Fig.3 The prediction surface map of Nitrate

3 結 語

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

(3)交叉論證結果表明，與普通克里金相比(MS=0.019 7,RMSS=0.979)，協(xié)同克里金的精度更高(MS=0.018 8,RMSS=0.988)，其所生成的硝酸鹽濃度預測表面圖能夠為污染物的治理提供一定的指導意義。

參考文獻：

[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] 王仁鐸，胡光道． 線性地質統(tǒng)計學[M]．北京: 地質出版社，1987．

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

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

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

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

[10] 李俊曉,李朝奎,殷智慧. 基于ArcGIS的克里金插值方法及其應用[J]. 測繪通報,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.