基于地统计学的伊犁-巩乃斯河谷渗透系数空间变异性研究

宋 浩 ， 张 琛， 刘 明, 曾 磊， 李 瑛， 马稚桐

(1. 长安大学 旱区地下水与生态效应教育部重点实验室， 陕西 西安 710061； 2.长安大学 环境科学与工程学院， 陕西 西安 710061； 3.中国地质调查局 西安地质调查中心， 陕西 西安 710054)

1 研究背景

水文地质参数的空间变异性是影响饱和及非饱和带地下水水流和溶质运移不确定性的主要因素[1]，是开展区域数值模拟、地下水资源形成演化以及水资源评价的基础，而这种空间变异性大多都具有很强的结构性[2]。20世纪70年代以来，自DAGAN率先提出了随机理论研究方法，使得空间变异性得以在数学模型的基础上进行解析，随后，空间变异性便成了研究领域的热点问题。国内外众多学者[3-16]注重于研究河床渗透系数、表土饱和渗透系数以及农业小流域的渗透系数，如Fu Tonggang等[7]和Botros等[8]运用实验变差函数拟合了喀斯特地貌下土壤饱和渗透系数的最优模型来表征该地区的空间变异性，并发现最优模型均是高斯模型。黄冠华等[11]通过测量试验地土壤相关的水文地质参数，发现浅层土壤的渗透系数大致呈对数正态分布。而施小清等[12]和商玮麟等[13]通过研究发现，承压含水层的渗透系数样本数据大多服从正态分布或者对数正态分布。

已有的研究中鲜有使用正态性更好的Box-Cox变换，承压含水层的空间变异性研究也比较少，并且传统的水文地质参数分区大多是以地形地貌作为参考依据。本文旨在对伊犁-巩乃斯河谷基于钻孔的渗透系数进行相关分析，探讨伊犁河谷潜水含水层渗透系数的空间分布规律，寻求一种利用有限水文地质资料获取区域渗透系数的方法，并为野外水文地质工作提供一个可供参考的渗透系数分区。

2 研究区概况

伊犁河谷区隶属于伊犁哈萨克自治州。地理坐标范围为80°9′42″～91°01′45″E，40°14′16″～49°10′45″N，总面积为19 900 km2。本文主要以伊犁-巩乃斯河谷作为研究区，研究区总面积约为10 000 km2，见图1。

研究区西端总宽度约为90 km，向东延伸，宽度逐渐减小。地层结构呈现单层、多层交替变化，岩层颗粒也粗细交替，并在最东端形成断陷谷地。在这些地层中，赋存有大量的浅层潜水以及深层的承压水，通过大范围的水文地质调查，山前冲洪积扇区地下水位大于100m，而河谷区两侧地下水位只有不足15 m。伊犁河谷区的地下水补给来源主要为北部以及南部山区的冰雪消融，整个谷地地下水总流向与伊犁河一致，为自东向西径流。

研究区的第四层地质大体可以描述为下更新统西域砾岩主要分布于南部库克苏河源头，中更新统冰水沉积主要分布在巩乃斯谷地两侧，上更新统的风积层主要分布于霍城以北，冲洪积层主要分布于乌孙山北麓以及特克斯以东，全新统风积层分布于伊宁市以西的伊犁河两侧。

图1 研究区位置与抽水试验点分布图

3 数据来源及方法

通过收集整理伊犁河谷地区的抽水试验资料，共获取52个潜水含水层和31个承压含水层的渗透系数数据。抽水试验点分布图见图1。本文用到的主要理论及方法包括Box-Cox变化、变差函数理论等。

3.1 Box-Cox变换

本文中由于潜水含水层渗透系数既不满足正态分布，也不满足对数正态分布，所以特引入Box-Cox变换从而提高了样本的正态性。 变换[17]就是对反应变量y进行变换,变换的公式为：

(1)

式中：λ为待定参数,λ>0;y为反应变量。

由公式(1)可以看出对于不同的λ会产生不一样的变换，很明显，λ=0就是对数变换，可见Box-Cox变换就是一种比对数变换更高级的变换。λ一般采用最大似然估计求得。

3.2 变差函数理论

变差函数是一个分析区域变量随机性与结构性特别有效的地质统计学方法。该方法假设区域变化量Z(x)是二阶平稳[18]的，则其变差函数定义为：

(2)

假设Z(xi)是平稳区域化变量Z(x)在n个点上的观测值，数据对{Z(xi),Z(xi+h)}为在某一方向上相隔|h|的点对(xi,xi+h)上的观测值，则定义Z(x)的实验变差函数为：

(3)

式中:γ*(h)为h点的实验变异函数;N(h)为数据对{Z(xi),Z(xi+h}的对数。

3.3 变差函数的理论模型

变差函数的理论模型[19]主要分为有基台值和无基台值模型两大类，本文主要研究的是有基台值的情况，所以就只罗列有基台值的三大模型，见表1。

表1中，各公式中a，C，C0分别为各自变差函数模型中的数值。

4 结果分析与讨论

4.1 渗透系数特征统计分析

对获取的样本数据进行统计分析，结果分别见表2以及图2、3。

从表2可以看出，该地区的潜水含水层渗透系数分布极度不均，渗透系数值分布在0.46～40.68 m/d之间，进一步表明研究区的潜水含水介质变化较大，既有类似砂砾等粗颗粒介质，也有类似于黏土等截水能力较强的细颗粒介质。承压含水层的渗透系数统计分析结果表明，渗透系数分布在0.54～20.40 m/d之间，相比于潜水含水层，这个值较稳定，主要原因是承压含水介质较潜水含水层来说，在整个地区的更替变化较小，地层结构更加稳定。从计算结果来看，潜水和承压水的原始K值以及变换后的K值都在0.1～1.0范围，但是很明显的区别是进行Box-Cox变换后，变异系数达到了最小，可见Box-Cox变换可以减小样本的离散度。

由图2、3可以很直接地看出来，首先是不论潜水还是承压水，原始K值直方图拟合的钟型曲线效果都不好，正态性比较差，主要还是由于K的样本离散度太大。此外，由Q-Q图(描述正态分布的散点图)也可以很明显的看出来，K值的正态概率分布稳定性很差，体现在样本点在斜率为标准差、截距为均值的直线上比较分散。而对K值进行对数转换和Box-Cox转换后，很明显地可以看出来钟型曲线拟合效果有所提升，尤其是Box-Cox变换，达到最优效果，这正体现出了Box-Cox变换对于样本正态性提升的优越性。由公式(1)来看，对数转换是Box-Cox变换的一种特殊形式，对于离散度较高的样本数据，Box-Cox变换则能显示出它的优势，这点从图2可以直接观察出来。但是在类似于承压含水层渗透系数稳定性较好情况下，对数变换能达到最优转换的效果。从表2可以看出，对数变换和Box-Cox变换的结果基本上是一致的，由图3也可以看出，两者的正态性相似。

运用K-S检验和S-W检验对K、lnK、以及进行Box-Cox变换后的K值进行正态性检验，结果见表3。

由表3可知，潜水含水层渗透系数只有进行Box-Cox变换后，样本显著性Sig值才大于0.05，说明潜水含水层渗透系数只有进行Box-Cox变换后，才会在95%的置信水平下服从于正态分布；如前面所言，由于承压含水层渗透系数离散度低，对数变换就可以服从正态分布。因此，对于潜水含水层渗透系数，分析过程中需要对原有数据进行Box-Cox变换，而承压含水层渗透系数，分析过程中只需要对其做对数变换即可。

4.2 实验变差函数的计算与分析

前面的研究发现伊犁-巩乃斯河谷的潜水含水层渗透系数服从于Box-Cox变换正态分布，承压含水层渗透系数服从对数正态分布。运用公式(3)计算实验变差函数[20]，以Box-Cox变换后的潜水含水层K值和对数变换后的承压含水层K值作为区域变化量，在进行基本的统计之后，发现样本值较为平均，符合平稳假设的条件。运用GS+软件进行拟合，潜水含水层样本数据选用滞后距为143 412 m，间距为13 200 m，承压含水层样本数据选用滞后距为123 987 m，间距为10 300 m，方向容限都为22.5°的情况下，计算该地区的潜水和承压含水层样本值的实验变差函数，得到变差函数特征参数的计算结果与理论模型拟合结果，见表4。

表2 渗透系数K统计分析结果

图2 潜水含水层渗透系数直方图与Q-Q图

图3 承压含水层渗透系数直方图与Q-Q图

从表4可以看出来，潜水含水层和承压含水层渗透系数的拟合结果均表明球状、指数和高斯模型的拟合优度相差不大，但是三者的变程相差比较大，很明显，高斯模型的变程a都是最小的，并且C0/(C0+C)都是最高的，所以选择高斯模型作为拟合结果的最优模型。依据Iqbal等[21]所阐述的观点，如果C0/(C0+C)<25%，说明系统的空间相关性很强烈。对比计算结果，可以看出，潜水含水层和承压含水层渗透系数的空间相关性都很强烈，依据相关研究[22]可知，这些变异主要由结构性因素所引起，也就是受地层结构的影响较大。图4给出了实验变差函数及理论模型的拟合图。

表3 渗透系数正态性检验结果

表4 潜水与承压水K值变差函数特征参数的计算与拟合

图4 实验变差函数及理论模型拟合图

结合表4与图4来看，潜水含水层的高斯模型拟合优度达到了0.378，C0/(C0+C)达到了18.85%，说明潜水含水层K值的Box-Cox变换空间相关性较高；变程a=10 046 m，所以当两个样本点的距离小于10 046 m时是空间相关的。承压含水层的高斯模型拟合优度达0.441，C0/(C0+C)达到了9.20%，说明承压含水层K值的对数变换空间相关性较好；变程a=9 180 m，所以说当两个样本点的距离小于9 180 m时是空间相关的。

4.3 渗透系数参数分区

4.3.1 潜水含水层渗透系数参数分区 利用克里金插值法可以得到该地区潜水含水层K值的等值线图，插值结果见图5。

图5中由于潜水含水层只有Box-Cox值服从正态分布，所以潜水含水层插值数据采用的是渗透系数的Box-Cox变换值，相应地，承压含水层采用的是K值的对数变换值。

由图5的插值结果可以看出来，研究区潜水含水层的渗透系数呈现出霍城县以西的渗透系数较高，而在伊宁县附近渗透系数突然变得很小，产生了整个研究区的最小的渗透系数区。察县和喀什河与伊犁河交界点偏南部的渗透系数又很高，而在整个巩乃斯河附近的渗透系数变化较稳定，渗透系数呈现出西部偏大而东部偏小的特点。结合研究区的地貌类型见图6。

由地貌类型图6可以看出，该地区的地层岩性从西向东主要是呈现出颗粒越来越细的特点。有研究表明，大尺度的渗透系数的空间变异性都是以地貌类型条件作为控制性因素，正是由于地貌类型的千差万别，导致河流的水动力条件发生了改变，从而使得沉积物岩性在空间上分布不均，直接导致渗透系数的大小也分布不均。研究区的砾质细土平原区占地面积较大，也就是说，较粗颗粒的地貌类型占研究区比例较大，正因为如此，研究区总体应该呈现出西部地区的渗透系数大于东部地区的渗透系数的特点，这与克里金插值结果一致。图5中渗透系数最小的区域为图6中的山前黄土丘陵区，图5中的3个渗透系数最大值，全都在砾质平原区。东部地区由于主要地貌类型为山前黄土丘陵以及含砾细土平原区，所以渗透系数相比于西部地区普遍较小。结合图5、6，对研究区的潜水含水层的渗透系数进行参数的分区，分区结果见图7。

图7的分区主要是依据图5所示的渗透系数的插值结果以及研究区的地形地貌和水文地质条件，R1～R8表示潜水含水层渗透系数各个不同的分区，各分区的渗透系数值是按照插值结果的平均值进行确定的，结果见表5，由于研究区的样本数据比较少，只能分出一个大概的框架。如果要进行比较细致的水文地质工作，需要收集更详细的资料来对研究区进行更加细致地划分。

表5 潜水含水层渗透系数K分区结果 m/d

4.3.2 承压含水层渗透系数参数分区 利用克里金插值法可以得到该地区承压含水层K值的等值线图，见图8。从图8可以看出，整个区域渗透系数最大值点分布在察县东部以及新源县北部偏西处，最小值区域分布在霍城县西部。承压含水层的渗透系数主要受地区深部地层结构以及沉积物类型的影响较大。经过调查得知，研究区内上游主要为卵砾石层，山前洪积扇中部以下和各大河流中下游逐渐变为砂砾石、黏土等互层结构。由此可见，该地区的深层沉积物主要表现为从东至西颗粒逐渐变细，且沉积物厚度越来越大，从而导致了该地区的含水介质渗透系数呈现出东部大于西部的特点。对研究区的承压含水层渗透系数进行参数分区，分区结果见图9。

图5 潜水含水层渗透系数插值结果图 图6 研究区主要地貌类型图

图7 潜水含水层渗透系数分区图 图8 承压水含水层渗透系数插值结果图

由于承压含水层的沉积物类型没有可供参考的资料，所以承压含水层的分区主要是依据插值结果图以及渗透系数的大小进行分区，R1-R6表示承压含水层渗透系数不同的分区，分区结果见表6。

图9 承压含水层渗透系数分区图

表6 承压含水层渗透系数K分区结果m/d

5 结 论

采用地质统计学对伊犁河谷的潜水和承压水含水层渗透系数进行了正态性检验，用实验变差函数理论进行了模型拟合，结合区域地貌地形以及克里金插值结果对该地区的渗透系数进行了分区，取得了如下认识：

(1)伊犁-巩乃斯河谷地区渗透系数样本数据经Box-Cox变换后符合正态性分布；而承压含水层渗透系数样本经对数变换即可满足正态性分布。

(2)伊犁-巩乃斯河谷地区的渗透系数显示出了强烈的空间相关性，并且该地区渗透系数的空间变异性主要是由岩层的结构本身引起的。

(3)基于克里金插值方法，结合研究区地形地貌特征，开展了研究区渗透系数分区工作，分区结果表明：研究区总体应该呈现出西部地区的渗透系数大于东部地区的渗透系数的特点，与克里金插值结果一致。分区结果显示该地区总体渗透系数较大，最大值在R3区，达到了31 m/d；最小值在R5区，仅为0.7 m/d。

(4)根据调查结果，承压含水层的地层岩性基本表现为从东至西逐渐变细的特点，相应地其渗透系数也在逐渐变小，这与插值结果也比较符合。分区结果显示承压含水层渗透系数总体偏小，最大值在R6区，可以达到13.5 m/d；最小值在R1区，仅为1.1 m/d。