Page 101 - 2023年第54卷第10期

分析，描述全局和局 部 的 统 计 特 性，建 模 方 便、精 度 较 高。Ｙａｎ等 ［３］ 针 对 理 想 均 质 含 水 层，采 用
Ｋｒｉｇｉｎｇ法建立了非稳定流数值模型的替代模型，验证结果表明建立的替代模型具有较高精度。替代模
型与 ＭＣＭＣ结合可提高参数优选的计算效率，如张江江 ［１３］ 采用自适应稀疏格子插值法构造了替代模
型，并与 ＭＣＭＣ结合优选监测井分布、反演模型参数。
另一种参数反演途径是数据同化，通过融合观测数据与模型预测值信息，自动调整模型参数、状态
和运行轨迹，实现对模型参数的反演估计 ［１４］ 。常用的数据同化方法包括集合卡尔曼滤波 ［１５］ （ＥｎＫＦ）、集
合平滑器 ［１３］ （ＥＳ）和多重数据同化集合平滑器 ［１６］ （ＥＳ － ＭＤＡ）等。ＥｎＫＦ使用当前时刻观测值对上一时
刻的参数和状态进行更新的实时算法，需要多次调用地下水模型，计算过程较为繁琐，且 ＥｎＫＦ对模
型参数和状态同时更新，可能导致更新后的参数和状态的不一致性。ＥＳ作为一种批处理算法，一次
性使用所有时刻的观测值更新模型参数和状态，且只更新参数以避免参数和状态的不一致性。ＥＳ －
ＭＤＡ通过多次迭代进行数据同化，可有效地避免了观测数据过拟合，具有调用模型次数少、计算成本
低的优点。陈冲等 ［１６］ 基于 ＭＯＤＦＬＯＷ 构建的黑河流域中游地下水模型，使用 ＥＳ － ＭＤＡ方法反演了模
型中 １３个参数最优取值，对比了 ＥＳＭ－ ＤＡ不同超参数对反演结果的影响。周念清等 ［１７］ 利用 ＥＳ － ＭＤＡ
方法融合高密度电阻率法采集的观测数据，实现了对污染源参数和渗透系数场的联合反演。然而 ＥＳ －
ＭＤＡ也存在不足，如敏感度低的模型参数反演精确度低。
目前常用的地下水数值模拟软件一般包括 ＧＭＳ、ＶｉｓｕａｌＭＯＤＦＬＯＷ 和 ＦＥＦＬＯＷ 等。这些软件具有
较为实用的交互界面，可以直接在图形界面构建模型，由于它们无法与优化算法直接耦合，在过去的
研究中更多倾向于先优化后模拟，即先选出待优化参数的可能取值，再将参数输入地下水模型进行验
证，还有研究采用松散耦合进行参数反演 ［１８］ 。这些方法无法利用计算机内存在优化算法和模型之间快
［１９］
速共享数据，极大地降低了计算速度和求解效率，ＦｌｏＰｙ 是由 ＭａｒｋＢａｋｋｅｒ等基于 Ｐｙｔｈｏｎ编写的第
三方库，可以和其他优化算法紧密耦合反演地下水数值模型参数。
针对地下水数值模型和参数反演各方法的优缺点，本文以复杂三维含水层（潜水、承压水）为案
例，假设某一含水层渗透系数为异质性空间变化场，在各含水层开采量未知情形下，建立地下水数值
模型和基于 Ｋｒｉｇｉｎｇ方法的替代模型，模拟含水层分层地下水水头变化，建立基于替代模型的 ＭＣＭＣ、
替代模型和数值模型相结合的两阶段 ＭＣＭＣ以及 ＥＳ － ＭＤＡ反演方法，对比各方法对参数和随机场反
演精度，为复杂含水层参数反演方法的选择提供借鉴。

２ 研究方法

２．１ 地下水数值模拟及参数反演问题 对于潜水、承压为水多层含水层，地下水流运动的一般表达式为：

 [ Ｋ （Ｈ － ｚ）  Ｈ ] [  Ｈ ] [ Ｋ （Ｈ － ｚ）  Ｈ ] ＋ ｗ１ ＝ μ  Ｈ

＋
＋
Ｋ （Ｈ － ｚ）
 ｘ ｘ  ｘ  ｙ ｙ  ｙ  ｚ ｚ  ｚ  ｔ
 [  Ｈ ] [  Ｈ ] [  Ｈ ] ＋ ｗ２ ＝ μ   Ｈ


＋
＋
ｙ
ｘ
ｚ
 ｘ ＫＭ  ｘ  ｙ ＫＭ  ｙ  ｚ ＫＭ  ｚ  ｔ
Ｈ（ｘ，ｙ，ｚ，ｔ）｜ ｔ ＝ ０ ＝ Ｈ（ｘ，ｙ，ｚ） （ｘ，ｙ，ｚ） ∈Ｄ，ｔ ＝ ０
０
，ｔ ＝ ０ （１）
Γ １ １
Ｈ（ｘ，ｙ，ｚ，ｔ）｜ ＝ Ｈ（ｘ，ｙ，ｚ，ｔ） （ｘ，ｙ，ｚ） ∈Γ １
 Ｈ
Ｋ ｎ   → ＝ Ｑ（ｘ，ｙ，ｚ，ｔ） （ｘ，ｙ，ｚ） ∈Γ ２ ，ｔ ＝ ０
 ｎ Γ ２
 Ｈ
，ｔ ＝ ０
  → ＋ α Ｈ ＝ β （ｘ，ｙ，ｚ） ∈Γ ３
 ｎ Γ ３
式中：Ｋ、Ｋ、Ｋ分别为潜水或承压含水层渗透系数，ｍ?ｄ；Ｈ、ｚ、Ｍ 分别为含水层水头、潜水含水
ｙ
ｘ
ｚ
层底板高和承压含水层厚度，ｍ；ｗ１、ｗ２分别为潜水、承压含水层源汇项（如入渗补给、开采强度），
－ １

０ (
ｄ ；μ为给水度，无量纲；μ 为储水系数，无量纲；Ｄ为研究区；Ｈ ｘ，ｙ，ｚ ) 为研究区初始水位，
２
— １ ３ ７ —