Page 72 - 2023年第54卷第6期

图 ２ 案例 １人工增储海岛概念模型及两种计算方法结果比较

的含水系统动态变化。为了验证本文替代建模方法在潮汐动态及含水层非均质影响下的有效性和精
度，本案例针对海岛水平分层越流含水层在潮汐作用下的动态行为进行替代建模，并与 Ｚｈａｏ等 ［２５］ 的
海岛分层越流含水层非稳定流模型进行对比。不失一般性，取含水层层数为 ２，模型结构如图 ３所示。
２
参考 Ｚｈａｏ等 ［２５］ 的设定，取上下含水层导水系数均为 Ｔ ＝ ２０００ｍ ?ｄ，上层潜水含水层给水度 Ｓ ＝ ０．１ ，
１
下层含水层贮水系数 Ｓ＝ ０．００１ ，所夹弱透水层的越流系数 ｑ ＝ ０．１?ｄ ，海岛两侧潮汐波动表示为 Ｈ（０，ｔ） ＝
２
Ｈ（Ｌ，ｔ） ＝Ａ·ｃｏｓ（２ π ｔ），Ａ ＝ ０．６５ｍ为潮汐振幅，Ｌ ＝ ２００ｍ为海岛水平长度。
Ｚｈａｏ等求解了该模型上下两层含水层中各处水头随时
间变化的解析表达式（详见文献［２５］中式（１１））。建立替代
５
模型时，取行走元数量 Ｎ ＝ １０，水平方向离散步长 Δ ｘ ＝
５
１ｍ ，时间步长 Δ ｔ ＝ ０．０２ｄ ，即每个计算点 ｘ有 １０ 个
ｉ
｛（ Ｈ（０，ｔ），Ｈ（Ｌ，ｔ）），ｈ｝样本用于训练式（４），从而
ｉ
构建随机行走替代模型。图 ４（ａ）（ｂ）展示了解析解与替
代模型获得的上下两层含水层中不同位置处水头随时间
图 ３ 案例 ２海岛分层越流含水层模型示意
的变化（实线表示解析解，空心圆表示替代模型解），图
４（ｃ）（ｄ）展示了替代模型相对于解析解的相对误差。ｘ ＝ ０ 代表海岛边缘处，ｘ ＝ Ｌ?２ ＝ １００ｍ代表海岛中
心处；由于模型空间对称性和时间周期性，出于简洁性考虑，仅展示海岛左半边［ ０，Ｌ?２］上不同点的
水头一个周期内的变化。结果表明：（１）相比于下层含水层中潮汐波动几乎同步、无损传播，上层含
水层受到更强储水能力的缓冲作用，越向岛内水头振幅和相位受到的阻滞越明显，而随机行走替代模
型很好地捕获了含水层系统时空变化特征。（ ２）在上下两层所有计算点上，替代模型与解析解都吻合得
２
－ ３
很好，相对误差介于±６‰之间，确定系数 Ｒ≈０．９９９９８，ｅ ＝ １．６６ × １０ ，均表明替代模型具有良好的精
ＮＲＭＳ
度。值得一提的是，本例中水头参照解均为随时间变化的周期函数（周期为 １ｄ），而同样周期变化的边
界水头导致替代模型式（ ４）（其中 ｗ ＝ ０）也是严格的同周期函数，导致两种解的差值即替代误差也是周期
函数，但振幅和相位各不相同（图 ４（ｃ）（ｄ））。在计算耗时方面，Ｚｈａｏ等 ［２５］ 解析解的计算由于涉及双曲
－ ３
正弦函数及循环，耗时约为 ２．９７ｓ。而替代模型计算耗时约为 ２．５３ × １０ ｓ，但该替代模型训练耗时约
４４．０ｓ，因此当调用模型次数大于 １５次时，替代模型效率更高。值得注意的是，在精度允许的情况下
４
使用更少的行走元（例如 Ｎ ＝ １０）训练替代模型时，可显著降低训练耗时（详见下文 ３．４节）。
３．３ 案例 ３：海岛随机非均质含水层 实际海岛含水层往往具有比上述分区?分层均质模型更为复杂的
空间变异性。为此本案例在前例海岛分层含水层基础上进一步考虑了随机非均质特征，即含水层渗透
系数参考场取自地质统计学的空间随机场模拟实现，具有更为复杂的非均质结构。假设某海岛含水层
长度 Ｌ ＝ ２０００ｍ ，厚度 Ｌ＝ ３０ｍ ，离散为 ２００ × ３的二维网格，每个单元格大小为 １０ｍ × １０ｍ；含水层对
ｚ
数渗透系数 ｌｎＫ由 ＧＳＬＩＢ工具箱 ［２６］ 根据指数型变异函数生成，均值为 ４．９０，标准差 ３．６８，水平相关长
度为 ５００ｍ，竖直相关长度为 ３０ｍ，得到的渗透系数 Ｋ随位置在 １４～４６４ｍ?ｄ之间变化，分布如图 ５
所示。边界条件方面，海岛两侧潮汐波动表示为 Ｈ（０，ｔ） ＝ Ｈ（Ｌ，ｔ） ＝ Ａ·ｓｉｎ（２ π ｔ），潮汐振幅 Ａ ＝ ０．６５ｍ；

０
— ７ ０ —