区域浅层地下水硝酸盐含量评价的指示克立格法

2001年3月水 利 学 报SHUILI XUEBAO 第3期收稿日期:2000 01 10基金项目:国家重点基础研究发展规划项目(G1999011810);中澳合作ACIAR 项目(LWR1 96 164).作者简介:李保国(1964-),男,山西襄汾人,教授、博士,主要从事水土资源管理和资源环境信息技术的研究.文章编号:0559 9350(2001)03 0001 05区域浅层地下水硝酸盐含量评价的指示克立格法李保国1,胡克林1,黄元仿1,张凤荣1(1 中国农业大学土壤和水科学系,北京 100094)摘 要:在进行区域水土资源特性的分析中经常会发现一些特异值,用普通克立格法来求解,通常是不稳健的.为此,本文引入了地统计学方法中的指示克立格法,详细介绍了其基本原理和计算方法.通过对华北冲积平原区曲周县667km 2的139个点的浅层地下水井取样分析,发现硝酸盐含量中存在有特异值,应用指示克立格法对其进行了分析,并绘制了硝酸盐含量的分布图及其含量超过最大允许值的概率图.这为区域水土资源利用的风险性评价提供了一种新的方法,同时也为区域水土资源质量的决策管理提供了指导.关键词:指示克立格法;地下水;硝酸盐;地统计学中图分类号:P641 文献标识码:A通常在进行水土资源特性的分析中会出现一些特异值.这些特异值既非分析误差所致,也非采样方法等人为误差引起,而是实际存在于所研究的母体之中,虽然只占全部数据的极少部分,但作用却不容忽视.为此,地质学家和统计学家进行了多种有成效的研究[1~3],产生了非参数统计方法,又称无分布法,指示克立格法是其中的一种.它是在不必去掉重要而实际存在的高值数据的条件下来处理各种不同的现象,并且给出在一定风险条件下未知量的估计值及概率[4,5],侯景儒利用该方法估算了矿床金品位的储量[6].但目前将其用到区域地下水的环境评价中的研究少有报道.本文利用指示克立格法研究了区域浅层地下水中硝酸盐含量的空间分布,并且给出了其含量超过国家规定的最大允许值的概率分布图.1 基本原理1 1 指示函数[7] 设在区域D 上取样并测定了某观测项目,若区域D 上的阈值或边界值为Z ,则在D 上的每个样品点x D 上定义1个Z 的如下阶梯函数:i (x ;Z )=1当x 点上的观测值Z (x ) Z 0当x 点上的观测值Z (x )>Z (1)在D 上的任一区域A D 内,低于阈值Z 的观测值Z (x )所占区域A 的比例表示为: (A ;Z )=1A !A i (x ;Z )d x [0,1](2)(A ;Z )是关于Z (x )和Z 的二元函数,即小于阈值Z 的全部i (x ;Z )(x A )的平均值.1 2 指示半变异函数 在给定的阈值Z 的条件下,随机函数I (x ;Z )服从二项分布,其期望值是:E {I (x ;Z )}=1∀Prob{Z (x ) Z }+0∀Prob{Z (x )>Z }=Prob{Z (x ) Z }=F (Z )(3)式中的F (Z )为观测值Z (x )的分布函数在Z (阈值)处的值.当I (x +h ;Z )和I (x ;Z )为被矢量h 分隔的两个随机变量时,则指示半变异函数可定义如下:1(h ;Z )=12E {[I (x +h ;Z )-I (x ;Z )]2}(4)1 3 指示克立格方程组及空间分布的指示克立格估计1 3 1 (A ;Z )的估计值 设在区域D 上有N 个有效数据,在D 上的一个域A (A D )内有n 个有效数据,{Z (x ),x A , =1,2,∃n },在给定阈值之后,得到样品的指示函数空间;{i (x ;Z ), =1,2,∃,n },则 (A ;Z )的估计值可表示如下:*(A ;Z )=%n =1! (Z )∀i (x ;Z )(5)1 3 2 指示克立格方程组 如同普通克立格法一样,为了求得 *(A ;Z )必须在无偏和估计方差极小的条件下求式(5)中的权系数! (Z )(=1,2,∃,n ).根据式(2)、式(3),偏差的期望可写成:E { *(A ;Z )- (A ;Z )}=%n =1! ∀E {I (x ;Z )}-F (Z )=%n =1! -1∀F (Z )=0(6) 当%n =1! =1时,满足上式,所以,无偏条件是:%n =1! =1.与普通克立格法相似[2,3],也可以给出指示克立格法的估计方差∀2E I :∀2E I =E {[ *(A ;Z )- (A ;Z )]2}=2%n =1 i (x ,A ;Z )- i (A ,A ;Z )-2%n =1%n #=1! !# i (x ,x #;Z )(7)要使估计方差在无偏条件下变为最小,这是一个求条件极值的问题,同样,采用标准拉格朗日乘子法,最后得到对于n +1个未知数(n 个!和1个∃)的n +1个方程的指示克立格方程组:%n #=1!# i (x ,x #;Z )+∃= i (x ,A ;Z )%n =1! =1( =1,2,∃,n )(8)式(7)和式(8)中:∃为拉格朗日乘子; i (x ,x #;Z )为在给定的阈值Z 条件下,矢量h 的两个端点分别在信息域x 、x #内的所有对点的平均指示半变异函数值; i (A ,A ;Z )为在给定阈值Z 条件下,矢量h 的两个端点在待估域A 内所有对点的平均指示半变异函数值; i (x ,A ;Z )为在给定阈值Z 条件下,矢量h 的一个端点在信息域x 内,另一端点在待估域A 内所有对点的平均指示半变异函数值.1 3 3 条件概率解释 可用条件概率将指示值i (x ;Z )表示如下:i (x ;Z )=Prob{Z (x ) Z |Z (x )=Z }(9)这时,某点x 的估计值i *(x ;Z )可表示为:i *(x ;Z )=Prob{Z (x ) Z |周围数据}(10)而某待估域A 的 (A ;Z )的估值 *(A ;Z )可表示如下: *(A ;Z )=1A !A Prob{Z (x ) Z }|周围数据}d x(11)1 4 待估域A 平均值的指示克立格估计 指示克立格法是应用某种克立格法(如普通克立格法)求得 (A ;Z )的线性估算值 *(A ;Z ),最后得到待估域A 的平均估计值[Z (x )]*[1].[Z (x )]*=%L l =1[I l (x )]*∀[Z (x )|x Z l ]*(12)1 5 指示克立格法的计算步骤 (1)对原始数据进行统计分析:包括样品分布特征的研究,绘制样品累积分布曲线等;(2)将原始数据变换成指示值:将原始数据变换成顺序统计量,做累积直方图,以10%,20%,∃,90%9个分位数(根据实际情况而定)对应的样品观测值{Z(x10),Z(x20),∃,Z(x90)}为确定的阈值Z l=(l=1,2,∃,9),求出每一样品对应于Z l(l=1,2,∃,9)的指示值并列成表格.这时,指示值的定义为:i(x;Z l)=1, 当Z(x)Z l(l=1,2,∃,9) 0; 当Z(x)>Z l(=1,2,∃,n)(3)计算每个阈值下的指示半变异函数 *l(h;Z l)(l=1,2,∃,9).为了减少工作量,在条件允许时,可用中位数(50%分位数)指示半变异函数进行研究,或只选那些有明显差异的指示半变异函数;(4)进行结构分析,给出不同阈值下的指示半变异函数模型;(5)求待估域的平均估计值和条件概率.2 实验方案2 1 实验方案 研究区设在河北省曲周县境内,面积为667km2,位于漳河冲积扇、漳河滏阳河冲积平原和黄河冲积平原的交汇处.经度36&35∋43(~36&57∋56(,纬度114&50∋22 3(~115&13∋27 4(,属暖温带、半湿润、大陆性季风气候区.年平均气温为13 1),平均降水量为556 2mm,降水主要集中在7 ~9月份,占全年降水量的2 3,年平均蒸发量是年平均降水量的3 3倍,地下水埋深一般在4m以下.作物种植制度主要为:小麦-玉米,小麦-大豆,小麦-棉花.小麦生育期内主要依赖地下水灌溉,林地、荒地均很少,有部分菜地.在该县境内农田灌溉的浅层水井里随机、但尽可能均一地提取地下水样,共进行139个点的采样,其地理位置(经度、纬度)用GPS定位,采样时间为1999年5月.2 2 测定项目及测定方法 (1)地下水位埋深:采取水样时用井绳标记地下水水面与地面的距离,用刻度尺量得水位.(2)地下水含盐量:使用DDS-11A T电导率仪测定.(3)地下水硝酸盐含量:取样后加1ml500ppm P MA抑制微生物活性,在冰箱中保存,用自动流动分析仪测定.3 结果分析3 1 观测项目的统计特征值 从表1、图1及图2可以看出:硝酸盐含量既不服从正态分布也不服从对数正态分布,其变异系数为2 25,属于高变异强度[8],并且含有特异值.因此,适合用指示克立格法进行分析.表1 曲周县地下水埋深及水质统计特征值观测项目样点数分布类型均值标准差最小值最大值变异系数地下水埋深 m139正态9 814 072 5518 970 414电导率 (mS c m)139对数正态2 941 430 766 640 486硝酸盐含量 (mg l)139偏态2 395 370 0136 282 253 2 指示值的确定 由于地下水中硝酸盐含量超过了一定的浓度就会威胁人体健康,世界卫生组织(W HO)颁布的饮用水质标准规定硝酸盐含量的最大允许浓度为10mg l.我国规定为20mg l,因而取20mg l作为阈值.根据图1试验累积频率分布图,可知阈值20mg l所对应的百分数为0 97,中位数0 50所对应的浓度值为0 3mg l.故选取0 3mg l和20mg l作为阈值.求出每个样品值所对应的阈值的指示值并列成表格,以便计算各个阈值的半变异函数.3 3 指示半变异函数的计算及结果分析 通过计算硝酸盐原始数据不同方向(0&30&60&90&4个方向)的半变异函数,得到了如图3的结果.从图3可以看出,曲周县浅层地下水中硝酸盐含量的分布在不同的方向上的差异不明显(不存在带状和几何各向异性),故可以认为其分布是各向同性的.图1 硝酸盐含量(mg l)累积频率分布图2 硝酸盐含量(mg l)频率分布图3 原始数据不同方向半变异函数图4 指示半变异函数(a阈值为0 3mg l b阈值为20mg l)根据式(4)可以得到不同阈值下的指示半变异函数值,用不同的模型来进行拟合,最后得到如图5的半变异函数模型.图5(a)的半变异函数模型为纯块金形式,表示在所取的采样密度下,各个点之间是相互独立的.图5(b)表示用球状模型拟合了阈值为20mg l的半变异函数,并给出了理论模型的3个参数.3 4 待估域的平均估计值及其条件概率 根据式(12)可以得到不同待估域的指示克立格的平均估计值(见图5).根据式(11)可以得到超过最大允许浓度20mg l的条件概率(见图6).从图5可以看出,地下水硝酸盐含量高的地方主要集中在曲周县的中西部和北部的一小块区域,这主要是因为中西部处于县城周围一带,这里容纳了大量富含氮素的生活废水,种植蔬菜的面积居全县之首,并且此区域也是流经曲周县境内几大河流(滏阳河,漳河)汇集之处,而这些河流的水质已遭到明显的污染,所有这些因素都会造成地下水中硝酸盐含量升高.北部的一小块区域为曲周县的第二大镇###河南疃镇,其生活污水,蔬菜种植面积居全县次位,故也造成了其硝酸盐含量较高.从图6可以看出,其浓度超过最大允许值20mg l的概率达到95%的地方,也主要是这两个区域.因此,就目前而言,应在这两个区域采取有效管理措施(如节水灌溉、科学施用氮肥等),对浅层地下水质量进行保护.将图5和图6与用传统的普通克立格法得到的区域硝酸盐含量分布图[9]对比分析,可知,在普通克立格法得到的曲周县浅层地下水硝酸盐含量分布图中,河南疃镇受污染的程度没有反映出来,这说明了普通克立格法具有明显的平滑效应,而采用指示克立格法更能准确地反映实际情况,且得到的信息更为丰富.图5 硝酸盐含量(mg l)分布图6 硝酸盐含量超过20mg l的概率4 结语在水土资源特性的数据分析中经常会出现特异值,这些值不可任意地舍弃,因为它们所起的作用往往是不可忽视的.普通克立格方法对此处理是不稳健的,用非参数的地统计学的指示克立格法就能够很好地处理这类问题,如本文对区域浅层地下水中硝酸盐含量的分析.应用指示克立格法不仅可以得到未知区域的估值,而且可以得到未知区域的条件概率,这就为进行区域水土资源利用的风险性评价提供了新的方法,也为区域水土资源质量的决策管理提供了指导.预计在环境污染指标的确定,区域水环境和土壤质量的评估等研究和实践领域都具有十分广阔的应用前景.致谢:本文得到了北京科技大学侯景儒先生的指教,特此感谢!参 考 文 献:[1] Journel A G.Non parametric estimation of spatial distribu tion[J].Math.Geol.1983,15:445-468.[2] 儒尔奈耳A J,等(侯景儒、黄竞先编译).地质统计学的理论及方法[M].北京:地质出版社,1989.[3] 侯景儒,黄竞先.地质统计学及其在矿产储量计算中的应用[M].北京:地质出版社,1982.[4] Journel A G.Nonparametric geostatistics for risk and additional sampling assessment[A].In Principles of Environ mental Sampling[M].Keith L H(ed.)American Chemical Society,Washington,DC,1988,45-72.[5] Goovaerts parative performance of indicator algorithms for modeling conditional probability distribution functions[J].Math.Geol.1994,26:389-411.[6] 侯景儒.指示克立格的方法及原理[J].地质与勘探,1990,26(3):28-36.[7] Henley Stephen.Nonparametric Geostatistics[M].Applies Science Publishers LTD,1981.[8] 雷志栋,杨诗秀,谢森传.土壤水动力学[M].北京:清华大学出版社,1988.[9] 胡克林,李保国,陈德立.区域浅层地下水埋深和水质的空间变异性特征[J].水科学进展,2000,11(4):409-415.Application of indicator Kriging method for assessing nitratecontent of regional shallow groundwaterLI Bao guo1,H U Ke lin1,HUANG Yuan fang1,Z HANG Feng rong1(1 China Agric ultural Unive rsity,Bei j ing 100094,China)Abstract:Some outliers usually are found in the process of analyzing the c haracteristics of soil and water resources.The universal Kriging method is not suitable to such data processing to obtain a steady solu tion.For this reason,a geostatistical method indicator Kriging method(a nonparamtric method)is intro duced in this paper.Samples collected from139observation points of667Km2in Quzhou County,North China Plain,are analyzed.Outliers are found in the nitrate content of groundwater.The indicator Kriging method is applied to analyze the data and form the distribution map of nitrate content and the probability map of its content exceeded a maximum allowed value.The proposed method provides a ne w way for as sessing the risk of utilizing soil and water resources in regional scale.It also can be used as the guide for decision making and management of regional soil and water resources.Key words:indicator Kriging method;groundwater;nitrate;geostatistics。