Page 70 - 2023年第54卷第6期

上式表明 ｈ 受到前一时刻水头 ｈ 以及当前时段相邻水头
ｉ，ｊ ｉ，ｊ － １
ｈ 、ｈ 、源汇项 ｗ 的共同影响，系数 ｐ 、ｐ 、ｐ 、β ｉ，ｊ 则
ｉ ＋ １，ｊ ｉ － １，ｊ ｉ，ｊ ｉ，ｊ － ｉ ＋ ，ｊ ｉ － ，ｊ
表明三者影响的相对强弱。实际上，总有 ｐ ＋ ｐ ＋ ｐ ≡１。式
ｉ ＋ ，ｊ
ｉ － ，ｊ
ｉ，ｊ －
（ ３）可看作为一个（ｘ，ｔ）二维格子上的随机行走规则，即：一个
假想粒子 （称为行走元） 位于网格点（ ｉ，ｊ）时，接下来有三种可
能选择，即移动至（ ｉ，ｊ － １ ）、（ｉ ＋ １ ，ｊ）、（ｉ － １ ，ｊ），相邻点移动
概率分别是 ｐ 、ｐ 、ｐ ；并且每经过一次（ｉ，ｊ），均受到网
ｉ ＋ ，ｊ
ｉ － ，ｊ
ｉ，ｊ －
图 １ 控制方程的随机
）。
格点（ ｉ，ｊ）上源汇项影响 （即式 （３）最 后一项：加上 ｗ · β ｉ，ｊ
ｉ，ｊ
行走替代示意
由于从点（ ｉ，ｊ）出发，无论到达三点中的哪个点，都要加上 ｗ ·
ｉ，ｊ
；行走元会在当前格点变量值为已知值（例如到达第一类边界上）时停止，否则继续行走。行走元
β ｉ，ｊ
移动到相邻网格点后有可能再次回到（ｉ，ｊ），则再次受到源汇项 ｗ 的影响。值得一提的是，对于非稳
ｉ，ｊ
定流问题，行走元可以在时间维上移动并穿越图 １中网格的不同行，有可能到达 ｔ ＝ ０对应的行（即初
始条件行）；对于稳定流问题，相当于令 Ｓ ≡０，自动得 ｐ ≡０，图 １中网格仅有 １行，行走元只在该
ｉ，ｊ －
时间行移动。采用大量行走元重复行走，统计这些行走元途经和停止的网格点，即得到出发点变量与
边界条件、初始条件、源汇项之间的平均定量关系，形式为：
ｈ ＝ Ｐ·ｈ ｋｎｏｗｎ β ·ｗ （４）
＋
ｉ，ｊ
式中：ｈ 为网格上边界条件、初始条件等已知水头值组成的向量；Ｐ为随机行走终点所确定的权
ｋｎｏｗｎ
重；ｗ为时空源汇项向量；β为源汇项权重向量，由随机行走平均途径次数确定。即每个计算点 ｘ有
ｉ
５
１０ 个｛（Ｈ ，ｗ），ｈ｝样本用于训练式（４）中的权重，（Ｈ ，ｗ）表示一类边界水头值 Ｈ 与补给强度 ｗ
Ｂ ｉｊ Ｂ Ｂ
的随机采样，ｈ为对应的解。Ｐ和 β实际上传递了含水层系统与水头 ｈ 之间的水力联系信息，这些信
ｉｊ
ｉ，ｊ
息由含水层结构、非均质性以及边界类型决定，与具体水头、流量数值无关。由此确定的式（ ４）可作
为地下水流高效替代模型。由于随机行走规则来自于控制方程（ ３）的离散形式，因此所得式（４）在数值
误差范围内满足式（１）。
该行走可以直接扩展到空间二维或三维问题，再加上时间维，使图 １中行走网格变为时空三维或
四维网格。以上步骤可使用规则或不规则网格的有限元等其他离散格式离散，也可直接扩展到其他类
型控制方程。需要注意的是，由于海水与淡水的密度差异，计算海水水头 ｈ 时需要先统一换算成等
（ｓ）
效淡水水头 ｈ ，即 ［２２］ ：
（ｆ）
? － ｚ? ε （５）
（ｆ）
ｈ ＝ｈ ·ρ （ｓ） ρ （ｆ）
（ｓ）
?（ － 为海水密度；ｚ为标高。为简便起见，下文提到的水头
式中：ε ＝ ρ （ｆ） ρ （ｓ） ρ （ｆ） ）；ρ （ｆ） 为淡水密度；ρ （ｓ）
默认已换算为等效淡水水头，且仅考虑咸淡水突变界面情况。
２．２ 随机行走 ＧＰＵ并行化 该随机行走替代模型方法具有显著的并行性优势，通过在不同计算单元
（ＧＰＵ或多核 ＣＰＵ）上同时模拟多个行走元的随机行走，可以极大减少训练时间成本。近年来 ＧＰＵ计算
发展迅猛，为随机行走模拟效率的提高提供了新驱动力。本研究中把 ＧＰＵ高效矩阵计算优势与本算法并
ｏｌｄ
ｎｅｗ
行优势结合，将大量行走元的随机行走处理为位置矩阵与一个随机位移矩阵之和，即Ｘ ＝ Ｘ ＋ Δ Ｘ ，其
Ｎ Ｎ Ｎ
中 Ｎ为行走元数量；Ｘ 为位置矩阵，其每行表示一个行走元当前位置对应的网格点；Δ Ｘ 为随机位
Ｎ
Ｎ
移矩阵，其每一行表示一个行走元从其当前位置移动到某个相邻点的相对位移，由相邻点移动概率
（式（ ３）中的 ｐ 等）随机取样获得。该运算可以在 ＧＰＵ上以极高效率完成，从而大大降低计算时间消耗。
ｉ，ｊ －
２
２．３ 替代模型精度评价指标 采用相对误差 ｅ、标准化均方根误差 ｅ ［２３］ 与确定系数 Ｒ作为替代模
ＮＲＭＳ
型精度评价指标，其计算公式分别为：
ｙ － ＾ｙ
ｉ
ｉ
ｅ ＝ × １０００‰ （６）
ｍａｘ（ｙ ）
ｉ
ｎ
２
∑ （ｙ －＾ｙ ）?ｎ
ｉ
ｉ
ｉ ＝１
ｅ ＝ （７）
槡 ｍａｘ（ ｙ ）
ＮＲＭＳ
ｉ
８
— ６ ９ —