地质统计学方法在地下水水位估值中的应用

1、应用协同-泛克立格法进行潜水位估值陈家军* 通讯作者：Fax: +86-10-62200397 E-mail address: jeffchen (Jiajun Chen)基金项目: 国家重点基础研究发展规划(973)项目(编号:G19990436)Foundation Item: the National Development Plan for Basic Studies (973 Project: G19990436); 教育部青年骨干教师资助计划项目the Excellent Young Teachers Program of Ministry of Education 郭乔羽 王金生

2、（北京师范大学环境科学研究所，环境模拟与污染控制国家重点联合实验室，北京100875）摘 要 在许多区域水资源水环境问题中，地下水流或水质模拟往往采用数值方法，需给出各个节点的初始水位值以反映流场的初始状态。本文阐述了，以容易获取的节点地表高程信息，作为协同区域化变量，用协同-泛克立格法进行地下水位估值的原理，以河南省焦作市修武段地下水数值模拟分析区为例，给出了利用该方法时的实验变差函数获取方法，分析了这一方法的估值情况，并与未加入协同区域化变量地表高程的泛克立格法进行了比较，指出协同-泛克立格法在进行松散多孔介质区域地下水位估值时能取得更为良好的效果。关键词 地下水，水位估值，协同-泛克立格

3、法，变差函数，协同区域化变量.1. 概述在区域水资源、水环境问题中，地下水水流模拟和水质模拟往往采用数值方法，需要在每个节点上给出初始流场值作为初始条件，反映流场的初始状态。通常情况下，剖分节点数总是多于统一测量地下水位的观测孔数，因而大部分节点上的初始地下水位要根据统一测量点的地下水位，通过插值求得。地质统计学方法克立格法是一种最佳空间估计法，其本质是最佳无偏估计，最初起源于地质矿产储量预测并在这一领域得到迅猛发展，近年来也广泛地应用于水资源、水环境问题的相关研究领域，如沼泽中化学物质循环的空间变异规律1（Marnette，Stein，1993），水域污染空间模式，海洋微生物2(Beliae

4、ff，Marie-laure，1995)，含水层导水率的分布3 (Chu et al., 1982)，污染物的浓度4 (Dingman et al. 1988)，区域降水量估值5 (Geren et al., 1994)，降水网络设计6 (Bastin et al., 1984)，流量估值7 (Huang et al., 1998)。利用协同-克立格在地下水位估值中也有一定的研究8-12 (Hoeksema et al.,1984, 1989; Yeh et al., 1995; Parks et al., 1996)，Hoeksema et al.(1989) 10在地下水位估值研究中所用的

5、协同-克立格方法是协同-普通克立格方法，其中忽略了地下水位变化的趋势特征。本文在泛克立格法的基础上，加入协同区域化变量地表高程，用协同-泛克立格法估计地下水水位，因为对于区域地下水(潜水)其从水位高处流向水位低处且往往与地形标高相一致，这种方法在可获得的国内外文献中尚未见到。提出了通过多次迭代获取剩余的稳健实验变差函数，从而获取区域化变量的稳健实验变差函数的方法。对协同-泛克立格法和单纯的泛克立格法在估值上的准确性进行了比较，表明采用协同-泛克立格法进行地下水位估值能取得更为准确的结果。2. 地下水水位估值的协同-泛克立格法在所研究的地下水水流系统中遍布观测孔是不可能的。进行水流数值模拟时，节

6、点数总是多于统一测量水位的观测孔数，而且还有一些统一测量水位点不能与节点重合，因而大多数节点上的初始水位要根据统一测量水位通过插值求得。在节点水位值信息不足的同时，节点的地表高程信息是容易从地形图等途径获得的。利用这些补充信息，使信息量增大，会使地下水位估值的准确性提高。对于松散多孔介质中有潜水分布的区域，地表高程与地下水位具有相关性，即地下水位往往随着地形的高低起伏而同步变化，尽管这不是绝对的，但这一变化趋势是存在的，因而可采用协同-泛克立格法进行地下水位估值。在用这一方法进行估值时，地下水位、地表高程作为协同区域化变量，其中包括结构部分和随机部分，而地下水位作为主估计量。通常区域流场和区域

7、地形都具有一定的趋势，也就是说其期望值不是常数而是空间坐标的函数，协同区域化变量地下水位和地表高程是非平稳的，存在漂移。假定协同区域化变量地下水位、地表高程满足： (1)式中：地下水位；地表高程；空间坐标（向量）；的数学期望值；漂移函数；方差；变差函数；k=1代表水位，k=2代表地表高程。设在流场内有个统一测量地下水位点，个统一测量地表高程点，则任一点地下水位的估计值由下式计算： (2)式中：任一点水位的估计值；协同-泛克立格权系数；统一测点地下水位值；统一测点地表高程值.漂移可写成13： (3)通常可表示成多项式。由于地下水初始流场在一般情况下水力坡度以及地形变化较为平缓，因而在本文中假定平

8、面二维流场，线性漂移，则可写成： (4)即，.根据无偏性条件：，经推导可得下面的无偏条件组： (5)而估计方差为： (6)要在无偏条件组下使估计方差达到极小，可用求条件极值的拉格朗日乘数法，即要求使： (7)达到极小时的诸、和。上式中、为拉格朗日乘数。将对诸、和求偏导，并令其为0，经推导得协同-泛克立格方程组： (8)式中K=M=在求得了诸、和（；）情况下，利用协同-泛克立格方程组式（8），则前面任一估值点的协同-泛克立格估计方差的表达式可进一步写成：（13） (9)只有得到变差函数才能进行后续的计算，因而在实际计算中，就必须获取稳健实验变差函数。然而泛克立格条件下稳健实验变差函数的获得是比较

9、困难的。但是我们知道当成立时，的变差函数就等于其剩余的变差函数，因为： (10)而剩余的变差函数应是稳健的。这样，求取的变差函数就转变为求取其剩余的变差函数。已有原始观测数据，若再能得到漂移，就可以求出剩余，从而求得其变差函数。因而我们通过估计漂移的泛克立格法2，多次迭代来获取稳健的实验变差函数。对于松散多孔介质，如岩性变化不大则可取成各向同性变差函数，用球状模型来拟合实验变差值。一般球状模型的公式如下： 0， ， (11) ，式中：为步长，为变程、为块金值、为拱高。 估计漂移时也需要变差函数为已知，为此初次计算时，变差函数取标准球状模型2，即、分别为0、1。根据球状模型的理论基础，如假定整个

10、区域的地下水位相互影响，地表高程相互联系，则初次计算时变程可由区域的最大宽度来确定1。 在初次计算中把通过标准球状模型计算得到的变差函数代入估计的泛克立格方程组2，求得漂移权系数并进而得到各个观测点的漂移估计值。由此可求得剩余： (12)再由下式计算得到剩余的实验变差： (13)通过对实验变差值的曲线拟合得到剩余的变差函数，即为实验变差函数，以此代替标准球状模型进行再次计算得到新的实验变差。不断重复这一过程，直到实验变差稳定下来，从而确定最终的稳健实验变差函数。地下水位与地表高程的最终稳健实验变差都确定之后，由下式得到交叉实验变差： (14)通过对交叉实验变差值的曲线拟合得到最终的交叉实验变差

11、函数。协同-泛克立格方程组中矩阵只与信息样品点的空间配置和变差函数有关，可首先由矩阵解出其逆矩阵，则根据协同-泛克立格方程组式（8），任一待估点的协同-泛克立格权系数和拉格朗日乘数由下式计算： (15)由此据式（2）可计算出任一点的地下水位估计值，由式（7）则可给出任一点的估计方差。 Nmm3. 实例分析以河南省焦作市修武段浅层地下水数值模拟分析区为例。修武段地下水数值模拟分析区处于一个相对完整的水文地质单元上山前冲洪积扇。模拟分析区浅层地下水为松散介质孔隙潜水。总体流向为西北东南方向。数值模拟分析区的范围为西北东南方向延伸1618公里，西南东北方向延伸约32公里，面积约图1 修武数值模拟分析

12、区剖分节点440平方公里。潜水含水介质主要为中砂、细砂，潜水含水层厚度为4080m，渗透系数为16.1321.33m/d，给水度为0.1620.221。可见区域潜水含水层岩性变化不大，可近似采用球状模型变差函数。区域地下水的补给主要来自降水补给和西北山区的侧向补给。地下水排泄的主要方式为农灌开采和侧向排泄。在数值模拟计算中对该区进行四边形和三角形混合剖分，共形成节点n=581个，见图1。1996年4月进行了一次地下水位统一观测，观测节点个数n=30个，见表1。进行估值的第一步是得到地下水位、地表高程的变差函数，以及两者的交叉变差函数。把所有30个观测点两两组合成点对，计算出其距离值。因为不同的

13、距离点对之间的相互影响力是有差别的，故按距离值的不同分为10组，见表2。表1 修武统一观测点节点号观测地下水位（海拔高度m）观测地表高程（海拔高度m）节点号观测地下水位（海拔高度m）观测地表高程（海拔高度m）1191.2099.9028283.5092.5044101.20105.0028683.0084.00102100.00117.8529280.5082.6010396.20114.0531379.0081.0010692.7095.0033479.1581.9012687.3089.0034480.0083.4016879.0086.0035683.10108.6018191.0092

14、.9037693.00110.6018385.0085.4538177.0080.1018875.0085.0045688.85100.0520780.00110.0048388.85126.0021973.6084.5050288.40116.0022383.00121.0053588.00123.0022882.7083.9055389.20118.0023182.7083.9056888.70125.00表2 修武统一观测点分组组别点对数平均距离(m)r10km23816607.21510kmr9km279543.56889kmr8km298480.88028kmr7km327450.7

15、7527kmr6km246517.59426kmr5km215519.09535kmr4km274486.17724kmr3km133601.70703kmr2km172564.4438r2km71504.6504以每组的平均距离作为该组的代表数据点，代入标准球状模型得出变差函数，进而用前面所述的方法求得实验变差，用一般球状模型式(11)拟合实验变差值得到实验变差函数。拟合即为寻找与实验变差值之间偏差最小的球状模型曲线，即分别对、求偏导，并令其为零，从而分别得到、的值，确定拟合曲线。以所得的球状模型曲线代替标准球状模型进行再次计算，得到第二次的实验变差值，与第一次实验变差值比较，如果任一点的两

16、者之差在以上，就继续拟合第二次球状模型曲线，再代替第一次的模型进行计算。直到任一点两次的实验变差值之差小于，这时认为实验变差值已稳定，从而拟合得到最终的实验变差函数。实验变差值稳定后，用前述方法得到交叉实验变差值。对交叉实验变差值的拟合与前面有所不同的是，从交叉实验变差值的图形看，认为其遵从指数函数模型，而非球状模型13： 0， ， (16)m故用指数函数模型拟合交叉变差函数。实验变差值拟合结果见图2。其中地下水位实验变差函数，即球状模型的参数为=10000；=0；=46.27024；地表高程实验变差函数，即球状模型的参数为=10311.75；=0；=170.3306；交叉实验变差函数，即指数

17、函数模型的参数为=49907.97；=0；=282.1266(见式16)。在地下水位实验变差值中，第十个点偏差比较大，认为是异常点而去除。mm图2 实验变差函数曲线拟合把拟合所得到的变差函数应用于前述的协同-泛克立格法中，即可计算估值点的地下水位。利用前面的方法原理编制计算程序10，用已知的30个地下水位观测点和地表高程观测点资料，对区域内581个节点进行估值，得到地下水位估计值、估计方差并绘制了相应的等水位线图，见图3、4。从图3可以看到，协同-泛克立格法地下水位估计方差的最大值在区域的东北角，这与区域内东北角的信息点最少，难以对估计值控制相符；除此以外，区域内西南角与西北边界上的估计方差也

18、较大，同样是因为信息的缺乏而对估值难以控制。从图4可以看到，水位估值所描述的流场方向为西北东南方向，符合区域地下水径流条件，因而估值所反映的模拟流场是可信的。协同-泛克立格法与单纯的泛克立格法相比较，更大限度的利用了容易获得的相关信息。由于信息量的加大，前者与后者相比，在进行地下水位估值时，应该更为准确。为了说明这一点，分别利用两种方法对区域内各个节点进行估值，比较其估计方差，见表3。此外，在30个已知观测点中抽取5个点作为假想的水位未知点，依据其余25个点，分别利用两种方法对其地下水位估值，再与观测值进行比较，从而说明方法本身的准确程度，见表4。这5个点的选取，主要是在地下水位与地表高程的变

19、化趋势较为一致的区域，以便更好地说明问题。从表3可以看出，采用协同-泛克立格法对修武段全部节点进行估值的平均方差和最大方差都小于采用泛克立格法进行估值的平均方差和最大方差，由此可以看出，协同-泛克立格法多利用了地表高程信息所得出的估值比泛克立格法所得估值准确。在表4中，从求五个已知观测点地下水位的估值数据可以看到，协同-泛克立格法所得的估计值要比泛克立格法所得估值准确，与真实值更为贴近。前者的绝对偏差比后者小一个数量级，估计方差也有所减小。表3 修武段581个节点地下水位估值主要参数表估值方法平均估计方差最大估计方差泛克立格法26.2518793.46481协同-泛克立格法25.2873783

20、.52544 Nmmmm N图4 协同-泛克立格法地下水位估值(m)图3 协同-泛克立格法估值方差 表4 泛克立格与协同-泛克立格地下水位估值比较节点号观测地下水位值(m)泛克立格法协同-泛克立格法水位估值(m)绝对偏差估计方差水位估值(m)绝对偏差估计方差16879.0078.24650.753519.026479.07030.070315.230028683.0081.41471.585320.200682.11180.888216.564131379.0080.18501.185028.069079.15510.155122.946938177.0080.59243.592429.452

21、577.90110.901123.715750288.4088.78890.388916.496888.37810.021913.56364. 结论进行地下水位估值时，可以利用地表高程信息，把地下水位与地表高程作为协同区域化变量，应用协同-泛克立格这一有效方法。通过计算可以看出本文所提出的，采用多次迭代获取剩余的稳健实验变差函数，从而获取区域化变量的稳健实验变差函数的方法有效可行。由对比分析可见，采用协同-泛克立格法进行地下水位估值能取得更为准确的结果。与单纯的泛克立格法相比，其更大程度地利用了容易获得的地表高程数据，使得信息量增大；同时由于信息量的增大也使得其估计结果也更为准确，与真实值更为

22、贴近，绝对偏差比用泛克立格法小了一个数量级，各点的估计方差也都有所减小。因此可以认为协同-泛克立格法是一种具有相对优越性的空间估计方法，在进行松散多孔介质区域地下水位估值时效果良好。参考文献1 Marnette, E.C.L., Stein, A., 1993. Spatial variability of chemical compounds related to S-cycling in two moorland pools. Water Research. 27, 1003-1012.2 Beliaeff, B., Cochard, M.L., 1995. Applying geosta

23、tistics to identification of spatial patterns of fecal contamination in a mussel farming area(havre de la vanlee, France). Water Research 29, 1541-1548.3 Chua, S.H., Bras, R.L., 1982. Optimal estimators of mean areal precipitation in regions of orographic influence. Journal of Hydrology 57, 23-48.4

24、Dingman, S.L., Seely-Renolds, D.M., Renolds, R.C., 1988. Application of kriging to estimating mean annual precipitation in a region of orographic influence. Water Resources Bulletin 24, 329-339.5 Garen, D.C., Johnson, G.L., Hanson, C.L., 1994. Mean areal precipitation for daily hydrologic modeling i

25、n mountainous region. Water Resources Bulletin 30, 481-491.6 Bastin, G.B., Duque, C., Gevers, M., 1984. Optimal estimation of the average areal rainfall and optimal selection of rain gauge locations. Water Resour. Res. 20, 463-470.7 Huang, W.C., Yang, F.T., 1998. Streamflow estimation using kriging. Water Resour. Res. 34, 1599-1608.8 Hoeksema, R.J., Kitanidis, P.K., 1984. An application of the geosta