水文地质学中的反演问题

地下水反演是一个综合性问题，研究重点是对含水层特征进行刻画。许多建模人员已经对概念模型的不确定性（显著的时间和空间变化）、尺度依赖性、许多类型的未知参数（导水系数、补给和边界条件等）、非线性以及对含水层性质经常表现出低灵敏性的状态变量（特别是水头和浓度）进行了研究。由于这些难题的存在，因此不可能把模型校准与模拟过程分割开来，而是应当将其作为认识含水层特征的重要一步。事实上，真实参数的估计方法尽管在计算细节上会有一些差异，但在本质上没有什么区别。就是否还有充足的空间能够改进地下水反演问题存在争议：开发用户友好的代码、通过地质统计学方法调整变异性、将地质学信息和其它类型的数据（温度、同位素浓度和地下水年龄等）相结合、正确计算不确定性等。尽管如此，即使利用现有的代码进行自动校准，也仍然极大地推动了模拟工作的发展。因此，主张将这一方法的应用进行标准化。一、概述广义上，反演模拟是指利用模拟参数的测量值来收集模型信息的一个过程，其中包括两个相关的概念：模型识别和参数估计，在本文中，后者与校准具有相同的含义。“模型识别”是应用一些方法来寻找模型的特征，如控制方程、边界条件、时间范围或非均质性等。而“参数估计”是对那些刻画特征的性质进行赋值。地下水模拟中的“反演”与上述定义大致相同。事实上，一些水文地质学反演问题已经建立了标准规程，例如，水文地质学家采用抽水试验来确定含水层的特征。解释这些试验结果的过程包括：首先确定最合适的模型（根据实际的测量值和地质条件），然后根据曲线拟合来估计模型参数，利用自动程序可以很容易完成这项工作。在解析解很少的领域中，如非饱和带，自动化的校准程序已经标准化了（Kool和Parker，1988；Hollenbeck和Jensen，1998）。在通过适当的测试刻画空间变化的领域中，如医学（扫描或X线断层摄影术）（Rudin等，1999；Rhli等，2002）和地球物理学（Sambridge和Mosegaard，2000）等等，自校准已经成为标准。然而，地下水模拟显示出了能够区别于其它模拟的几个特性：

成本地下水模拟费用相对较高。模型每次运行时，都要求解大型方程组（以获取系统准确而真实的特征）。另外，计算机科学的发展可以大大提高模型的性能；

时间依赖性状态变量，如水位和浓度具有时间依赖性。许多流动问题本质上都处于稳定流条件下，但即使如此，仍然会包含水位随时间波动的信息，这样就需要进行非稳定流的模拟；

非均质性水力传导系数K通常是最重要的水力学指标，其值会相差几个数量级；导水系数也是如此，相当于二维模型中的K（在下文中，这两个参数将会交替使用）；

不同类型的参数。当重点考虑导水系数时，其它参数（补给、边界流量等）同样具有不确定性和相关性；

尺度依赖性野外测量的参数通常只能表示含水层的一小部分，这样，在性质和数量上就与模型所需的参数具有很大的差异；

模型不确定性含水层的几何特征和非均质性受地质条件的控制，而地质条件很难得到确切的认识。

灵敏度低在不同的问题中，状态变量可能会对模拟参数表现出很低的灵敏度（即信息量不足）。特别是水头（最常用的，有时是唯一的测量指标）有时几乎不包含关于水力传导系数的信息。由于以上特征，含水层模型预测具有极高的不确定性。另外，参数估计不像其它领域那样可以明确地表达出来，在解释抽水试验结果时，就可能遇到这种情况。研究人员会将试验获得的水头值直接输入到程序中，经过适当的定性模型分析后，得到参数值。反演是含水层模拟过程中的一步重要工作。地下水的姐妹学科——地表水文学也具有许多类似的特征，其中的很多问题已经得到解决。地下水反演的另一个重要特点是其相对独立性。目前，已经开发出了很多独立于其它领域的反演方法。最早的方法是将简单的等效水头（假定已知）代入流动方程中，这样可以得到关于导水系数的一阶偏微分方程（Stallman，1956）。Neuman（1973）将这一方法称为“直接法”，这种方法相对简单、易于理解，在Nelson（1960，1961）之后，得到了广泛应用。事实上，根据测量水头画出流网后，可以得到导水系数（Bennet和Meyer，1952）。然而，这一方法存在几个不足之处：首先，整个研究范围内（时间和空间）的水头值（以及补给、贮水系数和边界条件）必须是已知的，这样才能进行插值工作，但同时也会引入误差，使最终估计的导水系数受到了一些人为影响；其次，这一方法不稳定，水头微小的误差会造成导水系数很大的误差。为了解决第一个问题，最近大部分的反演方法采用了Neuman（1973）所称为的“间接法”，这一方法承认测量存在误差，并寻找可以将这些误差最小化的水力性质。也就是说，通过对目标函数的最小化来寻找参数，这一工作的计算量会很大。为了克服不稳定性，采取了一系列方法：在目标函数中增加一个正则项（调整项）以消除无根据的摆动、考虑其它类型的数据、减少被估参数的数量（但不丢失表现空间变化的能力）等。这些问题（计算、稳定、不同类型的数据和空间变化等）一直是许多研究的焦点，如Yeh（1986）、Carrera（1987）、Kool等（1987）、McLaughlin和Townley（1996）以及deMarsily（1999）等人的研究。本文的目的是介绍含水层特征的反演模拟研究现状，重点是寻找现有方法的相似点，探讨如何去应用，而非是对这些方法的简单论述。校准是否应当作为模拟过程的必要步骤尚存在争议，一般的模拟过程包括：第一步，收集真实系统所有可用的信息，建立合理的概念模型；第二步，采用控制方程来表示概念模型；第三步，离散并求解控制方程，同时将一些未知的水力参数（即面状补给）离散化作为模型参数的函数；第四步，优化模型参数；第五步，如果需要的话，应当进行误差分析，根据新的试验结果将一些最初的概念模型舍弃，保留有效的模型。本文将在以下章节中进一步分析知识和数据在确定概念模型时的重要性，通常认为这是地下水模拟过程中最常见的误差源（Bredehoeft，2004）。在“需要估计什么”一节中介绍如何利用模型参数来表示含水层性质，下一步工作是优化（校准）模型参数。在“如何估计模型参数”和“代码如何工作”的章节中介绍了校准方法和程序。获得参数之后，就需要对模型的质量进行检验，这是“模型性能如何”一节的主要内容。在“实际做什么”一节中，将讨论实际的应用趋势。文末讨论了将来所面临的挑战。二、概念模型：知识和数据知识和数据是建立概念模型的基础。已经通过许多方法对“知识和数据”进行了定义。通常认为“知识”是对真实事物的认识。科学和哲学所关注的是这一认识是否客观而真实。这样，只有客观而真实的认识才被认为是“知识”。此处的“数据”是指关于含水层的所有信息，而不仅仅是野外测量的数据。地下水建模的基础是对地下水行为的科学认识和详细的场地资料，用数学方程来表达对研究对象的认识，通常会将这些方程处理成一系列可进行数值求解的离散方程。数据既可以定性又可以定量地定义概念模型。表1是一些数据类型及其在模拟过程中的作用。可以根据不同的标准对数据进行分类。其中一种分类方式是将数据划分为定性（“软”）和定量（“硬”）数据；另一种分类方式则考虑了数据的时间行为：静态数据指的是不随时间变化的量，而动态数据指的是随时间变化的量；第三种分类方式是将数据划分为被估参数和状态变量（水头和流量等）的测量值。表1水文地质学常用的数据分类类型实例应用知识达西定律和维度等控制方程软数据定性和定量（间接）地质（非均质性）、地球物理学、遥感、粒度分析、测量误差分布用参数表示、校准硬数据定量抽水试验、示踪试验、水头、化学组分的浓度、其它测量验前信息、校准、确认通过测量获得的定量数据会存在误差。因此，应当采用概率密度函数（下文称之为pdf）而非单一值来定义测量。在研究过程中，经常会采用高斯分布，因为已证实它适合用于描述很多变量（如导水系数的对数，下文简称log-T），另外，高斯分布可以采用平均值和标准差来描述。假设这一分布易于利用，经常要对不服从高斯pdf的数据进行转换，得到的pdf接近于高斯分布（如导水系数）。这是一种常见的误差源，因为log-T的点值服从高斯分布并不表示空间分布也服从多高斯分布（Gómez-Hernández和Wen，1998）。事实上，关于这一问题有很多争议（Meier等，1999），但Log-T的多高斯特征是水文地质学中最常见的假设。利用地质统计学软件来产生另一种可选择的空间分布，例如，采用马尔可夫链（Markov-chain）法（Weissmann等，1999）。状态变量原始测量值的转换也有助于增加数据对参数的灵敏度。在地表水水文学中，这一问题已得到了很多关注（Meixner等，1999）。在地下水研究工作中，相关的一个实例是模拟穿透曲线时，采用排泄点处的观测总量而非单个浓度值。也可以利用峰值浓度和峰值时间（如Woodbury和Rubin，2000）。另外，采用不同类型的数据或利用权重方法也可增加数据的灵敏度（如Wagne和Gorelick，1987；Anderman和Hill，1999）。在利用不同类型的数据时，没有通用的规则来确定数据的信息量（Larocque等，1999）。温度资料可以反映关于垂直流量的信息（Woodbury等，1987）。根据环境同位素和测年资料，可以了解区域流动趋势（Varni和Carrera，1998）；Chen等（2003）采用沉降率来帮助确定隔水层的渗透性。Hill（1992）采用河川径流的收益和损失来进行模型校准。总之，有大量的状态变量资料可以利用。可以采用不同方式将数据应用到反演过程中，通常考虑用包含状态变量的观测值h和模型参数p的矢量ω来表示大部分数据：

（1）

式中，h和p表示实际值，h*和p*表示测量值（或适当的转换值，p*通常称为“验前估计”），和表示误差。状态变量测量值的先验误差独立于那些参数的误差，这样的假设比较合理。事实上，所有测量误差都是独立的，假设通常都比较合理。上文提及的转换会造成数据之间的依赖性增加。对于模型参数更是如此，测量值与模型所需的数据绝不会完全吻合。为了获得p*，需要利用所有可用的参数资料，如参数的点测量值和软数据（如地球物理图像）等。用一个数据矢量zp包括所有这些测量值，之后利用zp，将p的验前估计和协方差矩阵用线性条件期望方程的一般形式来表示。\o"点击图片看全图"（2）\o"点击图片看全图"（3）式中，Qpp、Qzze和Qpze分别表示模型参数的协方差、zp中的测量值、模型参数和zp中测量值的交叉协方差；E[·]表示期望值。实际上，评价这些矩阵的方法与p和zp的性质有关，如果zp表示点测量值，则可以通过协方差函数获得Qzze（即Qzze是简单克里金矩阵）。在这种情况下，根据协方差函数也可以求得Qpze，如果p表示log-T的块值，则Qpze可能是平均值。如果p的先验期望E(p)、zp和未知，则需要扩展这些矩阵，来加强无偏条件。另外，如果zp不仅包括点测量结果，还包括软数据，则在Qzze中应当包括这些数据之间的相关性，Qzze会成为一个协克里金矩阵。由于交叉协方差很难估计，因此，通常采用普通克里金法、外部漂移克里金法或其它非平稳的克里金法。显然，计算这些矩阵非常复杂，详见一些地质统计学的文章，如Kitanidis（1997）和Rubin（2003）。Carrera等（1993）对如何处理大量不同类型的数据进行了逐步介绍。然而，无论如何构建这些矩阵，都必须明确以下两个问题：除了严格定性的数据（如地质描述），还需要提供关于状态变量和模型参数的所有类型的数据；对这些数据进行统计描述，以确定它们在反演过程中的权重。三、需要估计什么控制方程的系数是水力性质，需要对这些值进行估计。除了具有梯度或张量性质，所有这些水力性质都具有空间变异性，其中一些也随时间而变化。获得连续模型范围内每一点的水力性质是不可能的。因此，需要采用离散表示法。根据少量的模型参数来表示水力性质的方法称为参数化。选择模型参数并不是件容易的事，建模者往往采用很多自由度来对变异性进行准确描述。然而，如果模型参数量很大，那么反演可能会成为不适定问题（Hadamard，1902）。这样就需要折衷，这是参数化方法的动机。许多研究人员将含水层性质表示为未知模型参数的线性组合，这样，任何水力性质q（如面状补给）都可以表示为：\o"点击图片看全图"

（4）

式中，N是模型参数pi的数量，q0(x)是水力性质的附加因子（初始模型参数为0时的q的初始值），是将面状补给在时间和空间上进行参数化的内插函数。通过离散化可以将式（4）改写为矩阵形式：q=q0+Aqp（5）式中，Aq是面状补给内插函数矩阵，p是模型参数矢量。这样定义参数化方法的根据是：（1）对模型参数p的估计；（2）计算矩阵A的内插方法；（3）水力性质q0的初始分布。参数化方法不同于定义这些矢量和矩阵的方法，最常用的方法如下：（一）分区法将研究区分为若干个子区进行参数化。一般而言，模型参数的每个矢量组成pi与一个子区相关。这样在每个子区内，假定属性q(x)为常数或者规定以事先确定的方式变化，如果点x落在所研究的子区之外，则认为式（4）中内插函数的值为0。通常时间和空间变化是分开的（Carrera和Neuman，1986），因此：q=(x,t)=pifx(x)ft(t)（6）式中，fx是空间函数（例如，如果pi表示子区的水力传导系统，q表示导水系数，则fx可以表示含水层厚度）；ft是时间函数（例如，如果q是面状补给，ft可以表示根据土壤质量平衡获得的随时间变化的补给量，pi是未知系数）。分区的主要优点是在分析地质信息时的一般性和灵活性，例如，分区可以表示全部或部分地质单元等。然而，需要强调的是，分区并不排斥地质统计学方法的使用。事实上，Clifton和Neuman（1982）综合应用了分区法和克里金法。当子区不必很大时，在保证地质单元连贯性的同时，分区的目的就是降低p的维数（Stallman，1956）。事实上，Carrera等（1993）认为，关于参数变化的可利用的地质信息非常具有说服力（可以归类为确定性数据）时，比常规地质统计学方法要更胜一筹。有些学者认为分区法过于严格，因此，提出了对子区形状进行优化的研究。特别引人关注的是在校准过程中得到的子区形状的“几何变形”（Roggero和Hu，1998）。（二）点估计子区大小趋于0（事实上是趋于元素或栅格大小）的分区方式可作为分区法的特例，此时仍可以采用式（4）（如Meier等，2001）。然而，参数的空间维数就会很大，很难找到更合适的替换公式（Kitanidis和Vomvoris，1983；Dagan，1985；McLaughlin和Townley，1996）。下文会作进一步说明。（三）启发式内插函数可以任意选择式（4）中的内插函数αi。选择不同的类型，包括有限元（Yeh和Yoon，1981）和岭函数（Mantoglou，2003）等。这些方法非常灵活，但是在如何定义模型参数的先验信息方面不是很清楚。（四）试验点法在这种情况下，p表示在一组试验点中，属性q的未知值。这一方法最早由deMarsily等（1984）提出，现在已经得到广泛应用（Ramarao等，1995；Vesselinov等，2001；Hernandez等2003），已经成为非线性地质统计反演的标准方法。在试验点测量值不确定的情况下，可以考虑将式（4）进行一般化：\o"点击图片看全图"（7）

式中，和分别表示测量值和试验点的克里金加权系数。对式（7）和式（4）进行比较，显然，q0是等式右边的第一项，\o"点击图片看全图"与相等。这一方法的主要优点是相对简单而且灵活，然而，直到2004年才认识到需要考虑参数的先验信息。传统上，式（7）的第一项是通过直接测量获得的，但是在反演过程中不加以考虑。Doherty（2003）针对模型参数的不均匀性制定了调整标准，但没有用到先验信息。

（五）条件模拟到目前为止，所提出的方法都是基于寻求某种最佳估计。在以下章节可以看出，有时更倾向于寻求q(x,t)的可能模拟结果。Neuman和Wierenga（2003）提出了一种全面的模拟方法。Sahuquillo等（1992）、Gómez-Hernández等（1997）和Capilla等（1998）提出一种自校准方法，其中用q0(x)表示根据所有可用信息（包括软数据和硬数据）得到的log-T模拟结果，αi(x,t)pi表示由一组主点masterpoint造成的扰动。另一种方法是将q(x)作为随机函数的线性组合（Roggero和Hu，1998；Hu，2002），其中表示条件模拟，模型参数只是那些模拟的加权。总之，首先，需要表示水力性质的变化；其次，最常用的参数化方法的基本表达式如式（4）。现在的关键问题是如何找到模型参数p。四、如何估计模型参数p估计问题涉及到“最佳”模型参数的概念。问题是“最佳”意味着什么？关于这一问题没有理想的答案。事实上，可能没有哪一组模型参数可以代表“最佳”的现实情况，这也是条件模拟方法的动力——寻找可能的参数集合。下文将介绍一些最广泛使用的方法。（一）优化方法在这些方法中，将参数集合定义为目标函数的最小值。确定目标函数的最初目标是保证计算值和测量值之间可以很好地拟合。拟合度可以由计算水头和测量水头之间的方差来表示，在数学术语中，这个方差是一个范数（称为L2范数）。然而，也可以利用其它一些范数，而且有一些已经应用到地下水研究工作中。例如，可以用测量值和计算值之差的绝对值（L1范数）来量化误差。每个范数都有其优点：L2范数比L1范数更敏感，L1范数比L2范数更难应用（Woodbury等，1987；Xiang等，1992）。定义性能指标其实也就是一个最小化问题，即寻求使所选范数最小的模型参数。采用L2范数，将最小化目标函数写为：\o"点击图片看全图"（8）式中，是权重矩阵。当Fh作为目标函数时，问题经常是不适定的。目标函数的解也变得不稳定，即参数集不同，但得到的Fh值却非常相似。有时对初始参数表现出极度的敏感性，这样模拟人员会认为解是不唯一的。为了解决这些问题，一些研究人员选择设置估计参数的上限和下限。事实上，这样并没有使问题得到解决。解决的只是在上下限之间波动，而并没有提高可靠性。于是，Neuman（1973）在目标函数中引入一项Fp：F=Fh+λFp（9）（10）该等式最初是根据拟合和稳定性提出的。根据Tihonov（1993）的观点，可以将Fp作为正则项（调整项）。在水文学中，Emselem和deMarsily（1971）采用该项来降低摆动。当Ch作为测量误差的协方差矩阵时，可以根据统计平均值得到式（9）。Gavalas等（1976）将模型参数的后验概率密度函数最大化（最大后验法，MAP）来推导式（9）。得到的目标函数是λ＝1时的式（9）。Carrera和Neuman（1986）通过将参数可能性最大化（最大似然估计）推导式（9）。把式（9）放在统计学背景下进行叙述的优点是，得到的方法不仅能估计控制含水层性质的参数，也能估计那些控制不确定性的参数（方差、变异函数和似然性）。事实上，后者和前者同样重要（Zimmerman等，1998）。将目标函数最小化是一项艰巨的工作，这是由于状态变量和参数之间的关系通常是非线性的。这也是反演问题的公式划分为线性和非线性的原因。（二）线性方法不考虑用公式表示，可以将问题线性化为：h(p)=h0+Jhp(p-p*)（11）式中，Jhp是雅可比矩阵行列式（灵敏度矩阵），包括关于p的h的导数。后面的章节会对其计算方式进行讨论（灵敏度方程）。假设这个方程有效，可以将式（2）进行一般化，得到给定h条件下的p的条件期望（Dagan，1985；Rubin和Dagan，1987）。这样，就需要确定数据和参数Qph的交叉协方差以及Qhh的协方差矩阵，假定：

(12)

HYPERLINK"/uploadfile/200731819154341.gif"

(13)则：

(14)通过采用与式（2）相似的方式，扩展矩阵和，就可以将式（14）看作协克里金（Kitanidis和Vomvoris，1983；Hoeksema和Kitanidis，1984），是通过将估计参数的方差最小化得到的。最重要的是式（14）和克里金估计值与几何参数无关，一旦确定Qhh，就可以从理论上估计任一点的p。事实上，Kitanidis和Vomvoris（1983）认为需要估计的唯一参数就是表征log-T场的统计性质的参数。另外，也需要知道测量误差的协方差。（三）非线性方法根据式（11）引出一个问题：在参数先验估计时为什么要进行线性化？事实上，Carrera和Glorioso（1991）证明围绕估计参数本身线性化更好一些，这样更接近于真实值。在实际工作中，可以通过反复迭代进行线性化。Carrera和Glorioso（1991）研究表明，结果与将式（9）最小化的结果相同，式（14）的结果与将式（9）按照高斯－牛顿方法进行一级迭代进行最小化的结果相同。这样就形成一个闭合环。对提及的估计方法的总结见表2：表2估算方法总结方法类型评估对象运算法则（等式编号）研究人员最小二乘法非线性Fh（8）Cooley（1977）、Hill等（1998）最大似然估计非线性-2lnP(p/h)（9）CarreraandNeuman（1986）条件期望线性或非线性E(p/h)（14）Dagan（1985）、Carrera和Glorioso（1991）协克里金法线性或非线性方差估计（14）Kitanidis和Vomvoris（1983）、Carrera等（1993）最大后验法非线性P(p/h)（9）其中λ＝1Gavalas等（1976）各种方法的差异与统计参数的刻画和对不确定性的评价有关。本文没有讨论这些方法如何处理上面两个问题，而是着重分析它们的重要性。五、代码如何工作？许多代码基本上都是按照如下的迭代过程运行的，如：（1）初始化：读取输入数据，设置迭代计数器i=0，初始化参数p0；（2）求解模拟问题h(pi)，计算目标函数Fi及其梯度gi（假定是连续可微分的）和雅可比矩阵行列式Jhp；（3）计算校正矢量d，可能会利用到先前的迭代信息以及gi和Jhp；（4）更新参数，pi+1=pi+d；（5）如果已达到收敛，则停止计算，否则，取i=i+1返回到第（2）步。除了确定校正矢量d以及计算g和Jhp以外，其余步骤都非常简单。最常用方法如下：（一）计算校正矢量dd的计算属于优化范畴。关于数值优化方法的文献很多，但是只有少量应用于求解地下水反演问题（表3）。Cooley（1997）提出采用Marquardt方法（Marquardt，1963），这一方法的基础是把目标函数的二次逼近最小化，同时限制每次迭代的步长大小。这一方法效果非常好，通常经过少量迭代后就会收敛，但是由于需要计算雅可比矩阵，因此比较耗时。Carrera和Neuman（1986）提出将拟牛顿方法和共轭梯度法结合使用。与Marquardt方法相比，这一方法的效果要差一些，但每次迭代的计算量要小得多，仅需要计算目标函数的梯度。方法的选择取决于研究问题。Cooley（1985）对这两种方法进行比较，认为Marquardt方法要好一些。表3地下水反演模拟方法比较方法计算指标收敛级收敛1研究人员高斯－牛顿，MarquardtF2、g3和J42局部Cooley（1977）共轭梯度F和g1局部Carrear和Neuman（1986）模拟退火遗传算法F0全局Rao等（2003），Tsai等（2003）备注：1、全局收敛表示这种方法可以从理论上摆脱局部极小；局部收敛表示这种方法很难摆脱一个局部极小值；2、目标函数；3、梯度矢量；4、雅可比矩阵。大部分代码都利用了这些方法，一般称之为“下降法”，因为这些方法试图在每次迭代中改进目标函数。因此，它们往往易收敛于局部极小。克服局部极小值是许多方法的研究动机，如模拟退火或遗传算法（Tsai等，2003）。已证明在地表水水文学研究中一种非常有效的方法——SCE算法（ShuffledComplexEvolution）（Duan等，1992），但尚未应用在地下水领域中。这些方法在实际工作中尚未得到广泛应用，原因是：当参数数量很大时（大于20），方法的性能会降低，这在含水层模拟工作中经常出现；另外，由于需要对目标函数进行大量评价，因此导致地下水模型的成本过高。比较不同优化方法时，前述问题的复杂性具有一定的关联。例如，如果直接问题是非线性的，则与计算梯度或雅各比矩阵相比，评价目标函数的成本会更高，因为这些计算是无法迭代的。这样在解决非线性问题时，倾向于采用二阶收敛而非一阶或零阶收敛的方法。（二）灵敏度方程灵敏度是状态变量对模型参数的导数。这一方程非常有用，主要有两个原因：一是可以应用于上文提及的一些优化方法中，二是可以得到关于模型和估计参数可靠性的信息。另外，由于根据该方程可以确定哪些数据对哪些参数最灵敏，因此可以用于设计优化网络（见“模型性能如何”一节）。灵敏度一般可以通过以下三种方法获得：直接求导法、伴随状态法和有限差分法。Carrera等（1990）、Carrera和Medina（1994）对这些方法进行了比较。直接求导法是对模拟方程简单求导，在稳定条件下，可以得到：\o"点击图片看全图"（15）需要注意，在式（15）中，可以计算节点水头对参数的导数，但是采用有限元法（FEM）的内插函数可以计算空间上每一点的导数。同时也可以计算与水头有关的任何变量的导数。伴随状态法是根据将F的最小化作为关于h和p的优化问题，同时将状态方程作为等式约束。伴随状态矢量λ是联合优化问题的拉格朗日乘子的集合，可以根据下式获得：\o"点击图片看全图"（16）这样目标函数的梯度就可以写为：HYPERLINK"/uploadfile/2007318192059990.gif"

（17）伴随状态方程也可以用于计算雅可比矩阵，当观测点数量小于参数数量时，这一方法具有一定优势，也可以准确计算Hessian矩阵（目标函数对模型参数的二阶导数）（详见Carrera和Medina，1994；Medina和Carrera，2003的研究成果）。有限差分是根据增量比取近似导数求得的：\o"点击图片看全图"（18）采用有限差分时无需计算雅可比矩阵。直接问题的性质也会影响导数计算方法的选择。如果问题是非线性的（如非饱和流），则计算量会非常大，这样采用有限差分法的成本要比精确方法高（见表4）。然而，有限差分法比较灵活，越来越多的代码使用了这一方法，如UCODE和PEST（Poeter和Hill，1998；Doherty等，2002）。表4导数计算的方法比较方法优点缺点每次迭代成本线性问题非线性问题直接求导法精确难编程10111101伴随状态法精确难编程211n/a有限差分法容易编程不精确1011010备注：1、如果可以利用这样一个事实，即所有参数的系统矩阵是相同的，则可以减小精确求导和伴随状态法的计算成本；2、根据100个参数和20个观测点假设问题的模拟运行时间，估计高斯－牛顿迭代的成本；3、对于直接非线性问题（即非线性控制方程，如非饱和流问题），假定需要进行10次迭代。

六、模型性能如何？根据模拟流程，在对参数进行估计时，需要评价模型概念和估计参数的性能。模型质量受3个因素影响（Beck，1987）：（1）模型结构的不确定性；（2）模型结构中出现的参数值的不确定性；（3）系统预测的不确定性。此处的不确定性不仅指误差的随机波动，也指系统偏差。（一）模型结构的不确定性由于数据不完整和不一致，概念模型具有许多不确定性特征。建模人员被迫做简化假设。在参数化、离散化和选择边界条件等过程中会引入误差。解决这些不确定性经常需要引入几个概念模型。可以根据残差分布、参数相关和参数的置信区间等对不同的模型进行比较。理想情况下，一个好的模型应当与观测结果、不相关残差和合理的参数值相符。现在有几种模型符合这些标准，需要在其中选择最合适的模型。在时间序列分析和地下水应用研究中，已经确定了几种模型选择标准（Akaike，1974，1977；Rissanen，1978；Schwarz，1978；Hannan，1980；Kashyap，1982）。Carrera和Neuman（1986）应用这些标准中的其中4条进行综合测试研究，认为Kashyap标准最佳。Carrera等（1990）以及Median和Carrera（1996）也得到类似的研究结果。然而，研究人员可能会质疑只选择一种模型是否合适。如果其它模型与有效的数据相容但不采用这些模型的话，会不符合逻辑。这种质疑使得许多模型得到认可，并用来确定模型预测的不确定性（Beven和Binley，1992；Beven和Freer，2001）。（二）校准过程中的难点优化问题中最重要的难题是多解性、不可识别性和不（18）采用有限差分时无需计算雅可比矩阵。直接问题的性质也会影响导数计算方法的选择。如果问题是非线性的（如非饱和流），则计算量会非常大，这样采用有限差分法的成本要比精确方法高（见表4）。然而，有限差分法比较灵活，越来越多的代码使用了这一方法，如UCODE和PEST（Poeter和Hill，1998；Doherty等，2002）。表4导数计算的方法比较方法优点缺点每次迭代成本线性问题非线性问题直接求导法精确难编程10111101伴随状态法精确难编程211n/a有限差分法容易编程不精确1011010备注：1、如果可以利用这样一个事实，即所有参数的系统矩阵是相同的，则可以减小精确求导和伴随状态法的计算成本；2、根据100个参数和20个观测点假设问题的模拟运行时间，估计高斯－牛顿迭代的成本；3、对于直接非线性问题（即非线性控制方程，如非饱和流问题），假定需要进行10次迭代。稳定性。当不止一组参数能产生正问题的给定解时，会存在不可识别问题；当不止一组参数可以使目标函数最小化时，会存在多解性问题；当观测值的较小变化会引起估计参数的较大变化时，会存在不稳定性问题，但是通常可以根据求解初始参数确定。另外，优化算法，如Marquardt和共轭梯度法可能会陷入局部极小，不能表现出全局优化。Carrera和Neuman（1986）对这些概念进行了讨论，显示这些概念之间具有较强的相关性。事实上，无论是否是线性问题，都可以用不确定性来刻画。不稳定性和较大的不确定性是不同的概念，但是二者经常与延伸的置信区间有关，可以通过后验协方差矩阵的特征值和特征向量来刻画。该矩阵的特征向量构成了一组正交矢量，每个矢量都与一个特征值相关。与最大特征值相关的矢量表示具有最大不确定性的参数的线性组合，而与最小特征值相关的矢量表示不确定性最小。如果特征值差别极大，则很可能会出现不稳定性。估计参数之间具有某种联系，这样会造成参数的置信区间比假如参数是独立时的置信区间要大很多。解决不稳定性和大的不确定性的方法包括：

正则化（调整参数），这也正是在式（9）的目标函数中引入Fp的原因。Weiss和Smith（1998）认为先验信息是减少不确定性最有效的方法；

减少参数量，这是“需要估计什么”一节讨论的优化方案的动机；

增加数据量和数据类型，已在“概念模型：知识和数据”一节讨论过；

优化观测方案，可以设计观测网和试验来减小模型的不确定性，增加数据的可靠性，来区别对待不同的模型（Knopman和Voss，1989；Usunoff等，1992）。尽管如此，还是不可能得到问题的唯一解。这样就促使一些研究人员不采用估计，而使用随机模拟方法（如，Gómez-Hernández等，1997）。不确定性与所有模拟的总集，而非不确定性的统计特征有关。还有另一种可供选择的方法，Yapo等（1998）采用称为帕累托最优（ParetoOptimum）的概念来研究不稳定性，指出改进目标函数某一组分的参数矢量集合会造成另一个组分的失效。（三）预测系统未来行为的难点所有以上的难点会造成模型参数具有一些不确定性，因此模型预测也就遗传了这种不确定性。多尺度概念模型可能会进一步造成模型预测的不确定性。因此，评价预测的不确定性需要分析参数和模型不确定性的影响。后者通常可以根据模型模拟和评价预测范围来分析（如Medina和Carrera，1996）。当需要进行系统化时，由于所建立的概念模型是非系统的，因此很难进行系统化。这样就需要重点评价参数不确定性的影响。可以采用几种方法来量化预测的不确定性：线性近似、非线性近似和蒙特卡罗方法。线性方法相对简单，根据后验协方差矩阵进行计算：\o"点击图片看全图"（19）式中，Jhp是雅可比矩阵，Ch是测量协方差矩阵，Cp是参数协方差矩阵。如果f是一个模型预测值（f是参数p的一个函数），可以根据下式确定方差的下限：\o"点击图片看全图"（20)

式中，表示独立于参数不确定性的模型误差。该方程只是线性近似，因此实际不确定性可能远大于所得到的结果，这也是需要开发更为复杂方法的原因。然而，根据这一等式可以清楚地量化出不同的因素如何影响预测的不确定性。一般而言，预测不确定性与以下条件有关：（1）预测对模型参数和初始条件的灵敏度，显然，只有当预测对参数响应灵敏时，这些参数才会引起重视，可以根据求得。这也是有时会采用灵敏度分析而非误差分析的原因。（2）模型参数的不确定性（和/或初始状态），可以根据协方差矩阵Σp求得。非线性方法应用难度更大。Vecchia和Cooley（1987）提出计算非线性置信区间的方法，得到的结论如下：（1）线性和非线性置信区间通常是互相偏移或变化的，而且非线性置信区间通常要大一些；（2）非线性置信区间通常比线性置信区间的变化范围大；（3）随着区间的增大，线性和非线性置信区间的差异也随之增大；（4）先验信息可以改变置信区间的大小。Christensen和Cooley（1999）提出了一种量化模型非线性的方法。目前PEST的预测分析包括非线性置信区间。蒙特卡罗方法是计算能力最强的方法，基于不同参数集的正问题评价。这一方法的主要优点是易于理解，可以产生一个概率密度函数，无需进行复杂的假定；主要问题是很难确定所需的模拟量。实例之一是普适似然不确定估计（GLUE）方法（Beven和Freer，2001）。七、实际做什么？（一）应用趋势科技期刊中的文献往往集中在新方法和新解释的发展上，因此不适用于分析应用趋势，只有内部报告或专题研讨会才经常有关于应用的分析，因此，本节分析是基于作者的个人观点，会有一些偏差。可以链接HYPERLINK"/portal.cgi"/portal.cgi的“webofscience”，查询主要期刊上已出版的所有相关文章。根据统计，关于反演模拟的论文数量基本稳定（其中有5%是关于地下水模拟研究），而且在过去的13年里，关于反演模拟研究的论文数量在缓慢而稳定地增加。这些论文覆盖面很广，在前文中已经引用了一些。然而，很难确定应用趋势，唯一可以确定的是地质统计方法反演往往用来分析相对小范围的问题，如解释水力试验，而大尺度模型趋向于分析分区问题。可以将大多数地下水应用软件划分为这两类，将在下文作进一步讨论。（二）地质统计学反演在模拟含水层时，为了确定水力性质（特别是导水系数）的变化，应用到了地质统计学方法。然而，可以从两方面来认识这种需要：（1）约束模型参数；（2）模型预测。尽管这两点并不排斥，但很少放在一起使用。在方法上，地质统计反演遵循Clifton和Neuman（1982）最初提出的步骤，也就是说，首先提出一个随机模型（即不考虑Log-T是否稳定，变异函数和平均值是什么等）；其次，根据方程（2）和（3），采用可以利用的资料来产生一个前验，log-T的最佳估计及其协方差；最后，通过最小化方程（9）或利用方程（14），在前两步分析的基础上获得模型参数的估计。Kitanidis和Vomvoris（1983）、Carrera和Neuman（1986）认为需要对统计参数（误差方差和相关距离等）进行最优估计，因此对这一方法进行了改进。许多研究人员已经认识到这些参数的重要性。事实上，这也是Zimmerman（1998）对不同地质统计学反演技术进行比较之后得出的结论之一。野外数据的成功应用也仅限于相对小尺度的问题，其中包括水力试验解释（Yeh和Liu，2000；Meier等，2001；Vesselinov等，2001）、井捕获区描绘（Vassolo等，1998；Kunstmann等，2002；Harrar等，2003）以及其它应用（Barlebo等，2004）。大尺度含水层很少应用这一方法，这一事实说明，很难用地质统计学方法模拟大尺度含水层，同时也表明了地质统计学方法在应用过程中的两个局限性：不能充分考虑地质信息，不能重现真实的变化性。地质信息通常用术语（如沉积类型、导水特征的方向性等）来表达，这在反演过程中很难说明。这类信息在大尺度范围内可能会很准确，但在试验尺度内，则很少是准确的，在此过程中，唯一可以确定的是渗透性是变化的，这也是地质统计学所公认的。不过，在离散特征得到识别后，可以通过将其作为确定性的特征，而使反演方法得到明显改进，这一点也为Meier等人（2001）所认识，他们采用剪切带应力状态的相关信息来确定裂隙导水系数的各向异性。然而，地质统计学估计的主要问题是不能重现真实的变化。正如在“如何估计模型参数p”一节所讨论的，根据式（9）和（14）可以得到给定数据的log-T的条件期望，同样的，估计不能反映真实的偏差，因为期望值已经将偏差给过滤掉了。当资料充分时（通常需要经过详细的水力试验），可以从一定程度上确定变化类型（如Meier等，2001）；否则，只有估计协方差才能提供关于变异性和导水性质好的区域（溶质运移通道）或导水性质差的区域的信息。解决大尺度地质统计学反演反题的另一个方法是利用地球物理数据（Hubbard和Rubin，2000）。Gómez-Hernández及其合作伙伴（Gómez-Hernández等，1997；Capilla等，1998）使用条件模拟对以上问题进行了研究。理想情况下，所有这些模拟的平均值都等于条件估计值，但是每一个模拟都可以重现假定的变异性。这样，在预测对变异性比较敏感的过程（如污染物运移）时，采用模拟方法就适合得多。（三）基于地质条件的反演含水层尺度模型和试验尺度模型的主要区别是：确定变异性和参数类型的多样性时对地质条件的依赖程度。关于变异性的地质资料一般都不准确，但是不能忽略。不同地层或同一个地层中不同的地质单元，会具有不同的性质。这样，在描绘这些地质单元时，就不能忽略这种信息的价值。然而，很难准确认识地质单元的边界条件。因此，当认识到几何形状非常重要时，就要进行大量的工作对不同的几何形状进行研究。这个过程包括确定几何形状，并通过可利用的资料进行测试，重复进行直到满意为止。这一过程很单调，需要模拟人员和地质人员的互动。事实上，这个工作没有进行很好的备案，因此，如果几年后对模型进行修改，就不能确定当时为什么会选择这种模型结构。如果需要同时确定导水系数和其它类型的参数类型时，问题就会变得更复杂。根据作者的经验，补给（补给量和补给时间）、边界流量、抽水率（抽水率通常都会被低估）以及河流与含水层相互作用等方面经常存在不确定性。Cooley等（1986）和Hill（1998）尝试对这些问题进行研究。可以确定下面的一些趋势：

水力传导系数和导水系数的点值易出现误差。而且，当模拟尺度比抽水试验尺度大得多时，这两个参数用处不大。需要将这些测量值放在相关的地质条件下进行分析；

即使资料不完善，也要在模拟过程中考虑一些主要特征（导水裂隙带和古河道等）；

含水层活动的许多信息包括在离散事件（洪水和强降雨）中。要充分利用这些资料就需要采用瞬时模拟；

模型校准很少是唯一的（即不同的模型结构可能都会与硬数据拟合很好）。在模型预测时，需要承认这种不确定性。为了减小不确定性，就需要补充资料。这种方法具有几个缺点：一方面，没有正确考虑不确定性，模型参数的协方差矩阵不仅取决于硬数据，也取决于模拟人员的主观判断；另一方面，除校准条件外，无法保证模型的有效性。总之，需要对这个过程进行系统化，这一工作可以通过地质统计学来完成，因此，需要寻找能充分利用定性地质资料的地质统计学描述。八、下一步工作通过上述讨论，显而易见，在针对含水层刻画和管理的地下水研究工