地下水观测网优化设计的基本原理

发布网友 发布时间：2022-05-01 00:18

我来回答

共1个回答

热心网友 时间：2022-06-21 07:16

目前，地下水观测网优化设计主要采用的方法有时间序列法、水文地质学法、地质统计法以及一些最优化方法，这些方法在实际中已取得良好效果，促进了地下水观测网优化设计这一新兴交叉学科的发展。下面具体介绍地质统计法，包括地质统计学基础、普通克立格法、正克立格法、改进的克立格法及克立格法涉及的球状模型的拟合技术问题的研究，及其在地下水观测网优化设计中的应用。

3.2.1 地质统计学

地质统计学是一门新兴边缘学科。它已广泛地应用于地质勘探、煤田地质、石油地质、水文地质、工程地质、环境地质等地质学领域。在地质学领域的应用，包括矿体变化性估计、取样最优化、合理勘探方案的选择、资源评价的丰度估计、矿产资源的最优化估计、地下水位、地下水中化学组分浓度及含水层厚度等值线描述、含水层参数估计、地下水数值模拟的逆问题、地下水观测网合理布局等问题。

3.2.1.1 地质统计学基础

应用地质统计学方法设计地下水观测网时，要应用到两个重要的概念，即区域化变量和变差函数。区域化变量理论是地质统计学的核心，从随机变量的概念引申而来。区域化变量具有部分随机性、部分确定性特征。例如，地下水位、地下水质、含水层的渗透系数、给水度、孔隙度以及厚度等，都可看做区域化变量。传统的观点认为，地下水的特性（如水位、水质等）是“确定性”或完全给定的问题，即完全给定几何形状、参数、边界条件和初始条件。只要给定这些信息，地下水的特性可应用连续性方程和达西定律唯一求得。可是，在水文地质应用上，几乎不能精确地给出上述问题：诸如边界、初始条件和所有输入项。可现有信息的不完备性，以及采用带有测量误差的测量值来研究地下水问题，区域化变量理论为研究地下水问题提供了可能。地质统计学以区域化变量理论为基础，以半变差函数为主要工具，研究那些在空间分布上既有随机性又有结构性的自然现象。

（1）区域化变量

1）定义。以空间点x的三个直角坐标x、y、z为自变量的随机函数Z（x，y，z）=Z（x）称为一个区域化变量。区域化变量Z（x）的含义具有两重性，即观测前，把Z（x）看做随机函数；观测后，把Z（x）看做一个普通的三元实值函数（或空间点函数）。

2）特点。由于区域化变量是一种随机函数，因而能同时反映观测变量的结构性和随机性。例如，地下水某种化学组分的浓度分布，具有结构性和随机性的特征。结构性是指在地下水系统内，空间两个不同点x及（x+h）（此处h是三维向量（hx、hy、hz）T，它的模‖h‖=

，表示x与（x+h）点的距离）处的地下水水化学浓度Z（x）和Z（x+h）具有某种程度的相关性。随机性是指在地下水系统内，任意空间点x处，其地下水化学组分浓度Z（x）是一个随机变量。这就体现了区域化变量Z（x）的随机性特征。

（2）变差函数

变差函数是地质统计学分析的主要内容，它既能描述区域化变量的结构性变化，又能描述其随机性变化，同时它的计算又是地质统计学计算的基础。

1）定义。假定地下水系统内x处和（x+h）处的地下水水位分别为Z（x）和Z（x+h），是两个区域变量，h是两点间的距离，这两个区域化变量是相关的，这两个区域化变量之间变化性用变差函数来描述。变差函数可定义为x和（x+h）处区域化变量Z（x）和Z（x+h）之差的平均的数学期望，即

2γ（x，h）=E[Z（x）-Z（x+h）]²或γ（x，h）=1/2E[Z（X）-Z（x+h）]²（3.1）

在地质统计学的观测值中，同一点取样品只能取一个样品，只能得到一对Z（x）和Z（x+h），因而E[Z（x）-Z（x+h）]²、E[Z（x）-Z（x+h）]值无法求得。必须对区域化变量提出一些假设。

二阶平稳假设满足下列条件：

a.在整个研究区域内E[Z（x）]=m。

b.在整个研究区域内，Z（x）的协方差存在且相同。cov｛Z（x），Z（x+h）｝=E[Z（x）Z（x+h）]-E[Z（x）]E[Z（x+h）]=E[Z（x）Z（x+h）]-m²=C（h）

在实际工作中，二阶平稳假设不能满足，故提出本征假设。

本征假设满足下列条件：

a.在整个研究区域内E[Z（x）-Z（x+h）]=0；

b.在整个研究区内，增量的方差存在且平稳，即

Var｛Z（x）-Z（x+h）｝=E[Z（x）-Z（x+h）]²-｛E[Z（x）]-E[Z（x+h）]｝²=E[Z（x）-Z（x+h）]²=C（h）（3.2）

在本征假设下

γ（x，h）=1/2E[Z（X）-Z（x+h）]²

2）变差函数的性质。设Z（x）满足二阶平稳假设，则有γ（h）存在且平稳。于是，变差函数有以下的性质：

a.γ（0）=0；

b.γ（h）≥0；

c.γ（-h）=γ（h）；

d.[-γ（h）]必须是条件非负定函数（即-γ（x_i-x_j）构成的矩阵是条件非负矩阵）。

e.γ（h）与协方差C（h）的关系为γ（h）=C（0）-C（h）

式中：C（0）为验前方差。

3）变差函数理论模型。设区域化变量Z（x）满足本征假设，此时，变差函数理论模型有：

a.球状模型，一般表达式为

岩溶地区地下水与环境的特殊性研究

式中：a为变程；a₁为褴；a₀为块金常数。

b.指数模型

一般表达式为

岩溶地区地下水与环境的特殊性研究

式中：a为不是变程，3a是变程。

c.高斯模型

岩溶地区地下水与环境的特殊性研究

式中：a不是变程，

是变程。

d.幂函数模型

一般表达式为γ（h）=α₀+α₁h^λ，0 ＜λ＜2，当λ=1 时，γ（h）=α₀+α₁h，α₀为直线截距，α₁为直线斜率，为线性模型。

4）变差函数的计算。应用地质统计学进行实际计算时，重要的一步就是计算变差函数。一般根据地下水系统内实测点水文地质变量（水位、水质、导水系数等）值、计算实验变差函数值、绘制变差函数图、选取合适的理论变差函数，采取最佳的拟合技术，确定变差函数的表达式^[2]。

常用的公式是Matheron依据本征假设，提出的传统实验变差函数计算公式。

岩溶地区地下水与环境的特殊性研究

式中：N_h=N（N-1）/2，N为地下水系统内观测点数目，N_h为被距离矢量h分隔的区域化变量Z（x）和Z（x+h）数值对的数目；

（h）是实验变差函数。

Matheron指出，这种计算方法变差函数的可靠性随距离增大而减少，当区域化变量为非正态时，其估计效果明显下降；当采样点的分布极不规则时也就不稳定了。第二种方法是Cressie和Hawkins（1980）基于正态性假设提出的实验变差函数计算公式。第三种方法Henning More提出的稳健变差函数估计法，即

岩溶地区地下水与环境的特殊性研究

岩溶地区地下水与环境的特殊性研究

岩溶地区地下水与环境的特殊性研究

式中：â_i为待定的未知权系数。

在水文地质学应用中，一般采用第一种方法。

3.2.1.2 克立格法

（1）普通克立格法

克立格法是一种对时空分布变量求最优、线性、无偏内插估计量的方法^[3]。在水文地质方面，它可根据已知观测点上的实测数据，对观测变量进行结构性分析（变差函数确定）之后，对周围已知点的测量值赋予一定的权系数，进行加权平均估计待估点观测变量。

若Z（x）满足本征假设，则变差函数存在，在不考虑估值权系数非负约束的条件下，普通克立格法是通过如下的方程组来获得权系数及估计误差标准差的。

岩溶地区地下水与环境的特殊性研究

岩溶地区地下水与环境的特殊性研究

岩溶地区地下水与环境的特殊性研究

若令γ_ij=γ（x_i，x_j），且γ（x_i，x_i）=0，将式（3.9）写成矩阵形式为

γ·λ_μ=γ₀（3.11）

式中：λ_μ=（λ₁，λ₂，…，λ_N，μ）^T；γ₀=（γ₁₀，γ₂₀，…，γ_N0，1）^T

岩溶地区地下水与环境的特殊性研究

式（3.10）改写成

岩溶地区地下水与环境的特殊性研究

（2）正克立格法

在运用普通克立格方程计算克立格系数γ_i（i=1，2，…，N）时，往往计算出的克立格系数中会出现一些负值，没有要求λ_i（i=1，2，…，N）为正值。负的权系数的存在给出了一些不能令人满足的计算结果，造成地下水位（或地下水化学组成的浓度）的估计在某些点出现较大的偏差。因此，在应用普通克立格法时，当出现负权系数数值时，要对该方法进行改进，这里介绍一种改进普通克立格法的正克立格法。

若令基本的λ≥0，非基本的λ=0，基本的α=0，非基本的α≥0，即λ向量与α向量分别写成为

岩溶地区地下水与环境的特殊性研究

式中：下标“a”和“b”指基本的和非基本的向量。则正克立格方程组为

岩溶地区地下水与环境的特殊性研究

这一方程组具有N×1个方程和N×N未知数，可以直接求解。

正克立格方程求解步骤：

第一步，用普通克立格法求解，若求得的全部权系数λ_i（i=1，2，…，N）是非负的权系数，则终止算法，所求λ_i为最佳值。

第二步，若用普通克立格法求得的λ_i（i=1，2，…，N）中部分值为负，必须进行附加计算，按式计算，采用迭代法，在求解的过程中找出最小的权系数。

第三步，若得到的最小权系数是非负的，则算法把当前的解作为最佳解而终止。若最小的权系数是负的，则从当前的一组λ_i的基本值中消去对应的λ_i，这时的解（矩阵的矩和后来的权）不断被修改。

第四步，反复迭代，在解的过程中找到最小的权系数值，不断修改，直至所有的λ_i全部为非负为止。

（3）改进克立格法

在普通克立格方程组中，如果考虑权系数的非负性，则令：

岩溶地区地下水与环境的特殊性研究

这样，求解普通克立格方程组的问题就可转化成如下的一个线性规划问题：

岩溶地区地下水与环境的特殊性研究

令T=（t₁，t₂，…，t_n）^T，L=（λ₁，λ₂，λ_n）^T，M=（μ₁，μ₂）^T，Y=（L，M）^T，U为元素全为1的n维横向量，θ为元素全为零的n+2维横向量，I为n阶单位矩阵，A、B分别为

岩溶地区地下水与环境的特殊性研究

则线性规划可写为如下矩阵形式：

岩溶地区地下水与环境的特殊性研究

IT+AY≥B

IT-AY≥-B（3.18）

UL=1

岩溶地区地下水与环境的特殊性研究

（4）克立格法使用范围

克立格法在水文地质学研究中，具有广泛的应用前景，在用于地下水观测网优化设计方面，它具有以下几个特点：

1）应用变差函数描述水文地质变量的时、空的结构性和随机性，避免应用复杂的地下水系统数学模型。因此，该方法较简单，计算量小。

2）不需要更深入了解水文地质物理背景，可应用已有的地下水观测信息进行结构性分析。因而这种方法不涉及水文地质初始、边界条件及地下水系统的确定性参数和随机性参数。

3）应用实测数据确定变差函数时，需要较多的观测点（一般大于30个），以保证计算精度。

4）对于人工干扰较为强烈的地下水系统，克立格法难于考虑输入项的时空变化。

基于上述特征，克立格法的应用范围是：

1）可用于大区域优化设计地下水观测网密度的优化设计。

2）可用于人工干扰不大的水文地质区地下水观测网密度的优化设计。

3）克立格法结合数学规划法，可以确定给定经费约束下的地下水观测网的优化设计问题。

4）克立格法结合地下水系统数值模拟方法，可以研究任何条件下的地下水观测网优化设计。

3.2.2 加权线性规划法在球状模型拟合中的应用

对于变差函数模型的拟合问题一直没有很好的解决方法，以前人们通常采用人工拟合的方法，就是在充分考虑地质因素的基础上，根据实验变差函数曲线的特征，选择一定的理论模型，然后用直观的方法确定变差函数的参数。这种方法耗时、费力，而且缺乏统一客观标准；同时还影响了地质统计学整个计算过程在计算机上的自动化。目前已提出好几种拟合方法，如：非线性回归最小二乘法、加权多项式拟合法、线性规划拟合法、目标规划拟合法、遗传算法等。王仁铎教授等提出用加权多项式方法来拟合变差函数球状模型理论的模型参数，把变差函数自动拟合的问题向前推进了一步，但是没有解决理论模型参数的正负号问题。矫希国等提出用线性规划拟合变差函数球状理论模型的参数，由线性方程组非负解决理论模型参数的正负号问题；但该方法对各实验变差函数值等同对待，没有强调前面的几个数据点，而变差函数值估计可靠性是随距离的增大而减少的。近些时候提出的遗传算法，是求解非线性优化问题的有力工具，它具有全局寻优的特点，利用遗传算法能较好地拟合变差函数，但计算编程较为复杂。为此分别运用结合上述两种方法而提出用加权线性规划法对变差函数进行球状模型进行自动拟合，下面概括介绍这种方法的应用。

空间变量的地质统计学研究表明，在对实验变差函数进行拟合时，球状模型及其套合结构形式是最常用的理论模型。

对球状模型来说，拟合主要是针对0＜h≤a₁范围内的变异函数表达式：

岩溶地区地下水与环境的特殊性研究

进行拟合以求出模型中的参数C₀、C、a并使之满足C₀≥0、C＞0、a＞0的要求。将上式改写成：

岩溶地区地下水与环境的特殊性研究

令

岩溶地区地下水与环境的特殊性研究

则上式又可改写为

y=b₀+b₁x₁+b₂x₂（3.22）

鉴于C₀、C、a均有非负要求，所以也要求b₀、b₁、b₂≥0。如果在计算实验边变差函数时有m对数据：｛h_j，γ^*（h_j）（j=1，2，…，m）｝，可先将这m对数据变换成：｛（y_j，x_1j，x_2j）（j=1，2，…，m）｝，并令

t_j=|y_j-b₀-b₁x_1j-b₂x_2j|（j=1，2，…，m）（3.23）

按照最小二乘法原理，最佳拟合应满足

岩溶地区地下水与环境的特殊性研究

另外，考虑到在计算实验变差函数时，间隔h_j较小时参与计算γ^*（h_j）的数据对数目较多，计算结果有较高的可靠性和重要性，应拟合得好一些，精确些。随着h_j的增大，参与计算γ^*（h_j）的数据对数目相对较少，结果的可靠性则相应降低。因此，在进行理论变差函数拟合时，应使所拟合的理论变差函数在h_j较小时尽可能*近γ^*（h_j），h_j较大时误差可能大些。上述思想可通过对不同的t_j赋予不同的权值来实现。

岩溶地区地下水与环境的特殊性研究

式中：ω_j为权值，其值即可由下式计算，也可通过人机对话方式给出。

岩溶地区地下水与环境的特殊性研究

式中：N_j为h_j时计算γ^*（h_j）的数据对数目；A为放大系数。

再考虑到所拟合变量的非负要求，这样，实验变差函数的拟合就可表达为

岩溶地区地下水与环境的特殊性研究

t_j=|y_j-b₀-b₁x_1j-b₂x_2j|（j=1，2，…，m）

b₀≥0（3.27）

b₁≥0

b₂≥0（j=1，2，…，m）

上式是个线性规划问题，进一步可将其写成：

岩溶地区地下水与环境的特殊性研究

t_j+b₀+b₁x_1j+b₂x_2j≥y_j

-t_j+b₀+b₁x_1j+b₂x_2j≤y_j

b₀≥0（3.28）

b₁≥0

b₂≥0

t_j≥0（j=1，2，…，m）

令T=｛t₁，t₂，…，t_m｝^T，B=（b₀，b₁，b₂）^T，W=｛ω₁，ω₂，…，ω_m｝，U=（0，0，0），I为m阶单位矩阵，Y=（y₁，y₂，…，y_m）^T

岩溶地区地下水与环境的特殊性研究

求出b₀，b₁，b₂后，就可按下式得到模型参数值C、C₀、a。

C₀=b₀

岩溶地区地下水与环境的特殊性研究