CN116842691B - 一种智能提高地下水数值模拟收敛性的松弛方法 - Google Patents

Abstract

本发明公开了一种智能提高地下水数值模拟收敛性的松弛方法，包括以下步骤：构建地下水数值模拟模型；采用全局松弛迭代机制，提高所述地下水数值模拟模型的计算收敛性；对所述地下水数值模拟模型的计算收敛性进行判断，若不满足模拟需求则继续采用逐单元松弛迭代机制，进一步智能提高所述地下水数值模拟模型的计算收敛性。本发明一是引入逐单元松弛迭代机制，可智能确定每个地下水数值模拟网格单元之上的水头松弛程度；二是可灵活组合全局松弛迭代机制和逐单元松弛迭代机制；使用本发明能有效提高复杂地下水系统模拟的成功率和准确性，从而可为后续的地下水资源精细化管理与决策提供支撑。

Description

一种智能提高地下水数值模拟收敛性的松弛方法

技术领域

本发明属于地下水数值模拟技术领域，尤其涉及一种智能提高地下水数值模拟收敛性的松弛方法。

背景技术

近年来全球地下水开采量明显增加，目前全球约有50％的居民生活用水、25％的农业灌溉用水来源于地下水。尽管地下水如此重要，但由于其不可见的特性使得对地下水系统的研究有很大难度，目前对大范围地下水系统的研究主要依赖地下水数值模拟模型，准确的地下水数值模拟是地下水资源精细化管理与决策的重要支撑。

地下水数值模拟模型的本质是对描述地下水运动的方程组进行求解，以反演或预测地下水系统各处的水头。地下水数值模拟模型的计算收敛性是指对地下水运动方程组求解的难易程度，若地下水运动方程组易于求解，说明地下水数值模拟模型的计算收敛性好；反之说明地下水数值模拟模型的计算收敛性差。由于水文地质条件和人类干扰活动的复杂性，使得在许多地下水系统中，描述地下水运动的方程组具有高度非线性，因此难以求解，需采用相关技术提升模型的计算收敛性，也即提升模型模拟成功的可能性。

目前，国际主流的地下水数值模拟模型MODFLOW-2005采用了一种松弛迭代方法来提升模型的计算收敛性。所谓松弛就是指在使用迭代方法求解地下水运动方程组时，用前一次迭代求解的水头值修正本次迭代求解的水头值，从而减小搜索地下水运动方程组真解的步长(速度)，避免由于搜索步长过大而无法搜索到真解的情况，以达到提升模型计算收敛性的目的。该方法的主要缺点在于所有地下水数值模拟网格单元之上的水头松弛程度都必须相同，因此整个模型的松弛程度取决于最难收敛的地下水数值模拟网格单元，方法灵活性较低；同时随着模型松弛程度的加深，模拟结果的水量平衡误差会明显增大，因此当地下水数值模拟模型的计算收敛性较差时，使用该方法会明显降低模拟结果的准确性。实践证明该方法难以适用于各种复杂模拟条件的要求，更具鲁棒性、灵活性的松弛迭代方法对于拓宽地下水数值模拟模型的应用面、提升地下水数值模拟模型的准确性具有重要意义，同时也有助于支撑水资源保护和管理工作。

发明内容

为解决上述技术问题，本发明提出了一种智能提高地下水数值模拟收敛性的松弛方法，该方法具有很好的鲁棒性和灵活性，能提高复杂地下水系统模拟的成功率和准确性。

为实现上述目的，本发明提供了一种智能提高地下水数值模拟收敛性的松弛方法，包括以下步骤：

构建地下水数值模拟模型；

采用全局松弛迭代机制，提高所述地下水数值模拟模型的计算收敛性；

对所述地下水数值模拟模型的计算收敛性进行判断，若不满足模拟需求则继续采用逐单元松弛迭代机制，进一步智能提高所述地下水数值模拟模型的计算收敛性。

可选的，构建地下水数值模拟模型的方法包括：

基于地下水数值模拟网格单元系统结合所述地下水数值模拟网格单元系统的源汇项信息，根据水量平衡关系构建描述目标网格单元处地下水运动的差分方程，联立若干个所述网格单元的地下水运动差分方程构建所述地下水数值模拟模型。

可选的，所述地下水数值模拟网格单元系统由地下水系统空间离散后形成的具有指定行数、列数、层数的若干个网格单元组成。

可选的，所述水量平衡关系包括：在模拟时段内自周围六个相邻网格单元流入所述目标网格单元的地下水量与作用于所述目标网格单元之上的源汇项之和等于所述目标网格单元贮水量的变化。

可选的，根据水量平衡关系构建描述目标网格单元处地下水运动的差分方程的方法为：

其中，下标(i,j,k)为目标网格单元，下标(i,j-1,k)、(i,j+1,k)分别为目标网格单元沿行方向上左、右侧的相邻网格单元；(i+1,j,k)、(i-1,j,k)分别为目标网格单元沿列方向上前、后侧的相邻网格单元；(i,j,k-1)、(i,j,k+1)分别为目标网格单元沿层方向上上、下侧的相邻网格单元；CR、CC、CV分别为沿行、列、层方向上网格单元之间的水力传导系数；h为水头；上标m为本次迭代，上标m-1为上次迭代；P_i,j,k为作用于目标网格单元之上的与水头有关的源汇项相关系数；Q_i,j,k为作用于目标网格单元之上的流量源汇项；SCI_i,j,k为目标网格单元的贮水项；Δt为模拟时段的时间长度；Δh为模拟时段内的水头变化量。

可选的，联立若干个所述网格单元的地下水运动差分方程构建所述地下水数值模拟模型的方法为：

[A]{h}＝{q}

其中，[A]为系数矩阵；{h}为未知水头向量，代表各网格单元处的待求解水头；{q}为右端项向量，表示差分方程组中的常数项和已知项。

可选的，采用全局松弛迭代机制，提高所述地下水数值模拟模型的计算收敛性的方法包括：

其中，为地下水系统中任一目标网格单元直接通过地下水运动差分方程组求解得到的某次迭代的水头值；/>为目标网格单元经过全局松弛后上次迭代计算的松弛水头值；Damp为全局松弛因子；/>为目标网格单元经过全局松弛后本次迭代计算的松弛水头值。

可选的，对所述地下水数值模拟模型的计算收敛性进行判断的方法包括：

采用全局松弛迭代机制后，判断经过全局松弛后任一所述目标网格单元前后两次迭代计算的水头值相差是否都小于设定阈值，若经过若干次数迭代后任一所述目标网格单元前后两次迭代计算的水头值都小于所述设定阈值则所述地下水数值模拟模型的计算收敛，此时若模拟结果的水量平衡误差很小，表明模拟是成功的；若模拟结果的水量平衡误差较大或前后两次迭代计算的水头值相差始终无法在任一所述目标网格单元上都小于所述设定阈值则继续采用逐单元松弛迭代机制进一步智能提高所述地下水数值模拟模型的计算收敛性。

可选的，若不满足模拟需求则继续采用逐单元松弛迭代机制，进一步智能提高所述地下水数值模拟模型的计算收敛性的方法为：

其中，为任一目标网格单元经过全局松弛后本次和前一次迭代计算的水头的相对变化值；/>为逐单元松弛迭代机制下本次迭代计算时目标网格单元的自适应松弛因子；/>为逐单元松弛迭代机制下目标网格单元本次迭代计算的松弛水头值；/>为上次迭代计算时目标网格单元的自适应松弛因子值；α为自适应松弛因子的调减系数；λ为自适应松弛因子的调增系数；β为自适应松弛因子的单次调增值。

本发明技术效果：本发明公开了一种智能提高地下水数值模拟收敛性的松弛方法，一是引入逐单元松弛迭代机制，可智能确定每个地下水数值模拟网格单元之上的水头松弛程度；二是可灵活组合全局松弛迭代机制和逐单元松弛迭代机制；使用本发明能有效提高复杂地下水系统模拟的成功率和准确性，从而可为后续的地下水资源精细化管理与决策提供支撑。

附图说明

构成本申请的一部分的附图用来提供对本申请的进一步理解，本申请的示意性实施例及其说明用于解释本申请，并不构成对本申请的不当限定。在附图中：

图1为本发明实施例智能提高地下水数值模拟收敛性的松弛方法的流程示意图；

图2为本发明实施例基于网格单元的空间离散技术和源汇项示意图，其中(a)为本发明实施例的地下水系统沿行方向、列方向的网格化离散技术与模拟源汇项示意图，(b)为本发明实施例的地下水系统沿行方向、层方向的网格化离散技术示意图；

图3为本发明实施例任意一个地下水数值模拟网格单元与其周围六个相邻地下水数值模拟网格单元的示意图；

图4为本发明实施例的水头空间分布和等水头线模拟结果示意图。

具体实施方式

需要说明的是，在不冲突的情况下，本申请中的实施例及实施例中的特征可以相互组合。下面将参考附图并结合实施例来详细说明本申请。

需要说明的是，在附图的流程图示出的步骤可以在诸如一组计算机可执行指令的计算机系统中执行，并且，虽然在流程图中示出了逻辑顺序，但是在某些情况下，可以以不同于此处的顺序执行所示出或描述的步骤。

如图1所示，本实施例中提供一种智能提高地下水数值模拟收敛性的松弛方法，包括以下步骤：

构建地下水数值模拟模型；

采用全局松弛迭代机制，提高所述地下水数值模拟模型的计算收敛性；

对所述地下水数值模拟模型的计算收敛性进行判断，若不满足模拟需求则继续采用逐单元松弛迭代机制，进一步智能提高所述地下水数值模拟模型的计算收敛性。

S1、构建地下水数值模拟模型；

步骤S1将地下水系统空间上离散为具有一定行数、列数、层数的地下水数值模拟网格单元系统如图2所示；结合地下水系统的源汇项信息如图2(a)所示，对任意一个处于第i行、第j列、第k层的网格单元(i,j,k)，按照“在模拟时段内自周围六个相邻网格单元流入网格单元(i,j,k)地下水量(可正可负)与作用于网格单元(i,j,k)之上的源汇项之和等于网格单元(i,j,k)贮水量的变化”这一水量平衡关系如图3所示，建立如下描述网格单元(i,j,k)处地下水运动的差分方程，

其中，下标(i,j-1,k)、(i,j+1,k)分别代表网格单元(i,j,k)沿行方向上左、右侧的相邻网格单元；(i+1,j,k)、(i-1,j,k)分别代表网格单元(i,j,k)沿列方向上前、后侧的相邻网格单元；(i,j,k-1)、(i,j,k+1)分别代表网格单元(i,j,k)沿层方向上上、下侧的相邻网格单元；CR、CC、CV分别为沿行、列、层方向上网格单元之间的水力传导系数(L²T^-1)；h为水头(L)；上标m代表本次迭代，上标m-1代表上次迭代；P_i,j,k代表作用于网格单元(i,j,k)之上的与水头有关的源汇项相关系数(L²T^-1)；Q_i,j,k代表作用于网格单元(i,j,k)之上的流量源汇项(L³T^-1)；SCI_i,j,k代表网格单元(i,j,k)的贮水项(L²)；Δt为模拟时段的时间长度(T)；Δh为模拟时段内的水头变化量(L)。

对地下水系统所涉及的所有网格单元逐个列出以上差分方程，形成描述整个地下水系统水流运动的差分方程组(矩阵方程形式)，也即建立了地下水数值模拟模型，

[A]{h}＝{q} (2)

其中，[A]为系数矩阵；{h}为未知水头向量，代表各网格单元处的待求解水头，也即地下水系统各处的待求解水头；{q}为右端项向量，表示差分方程组中的常数项和已知项。

地下水运动差分方程组(矩阵方程形式)的系数矩阵[A]与当前待求解的水头值有关,因此地下水运动差分方程组无法直接求解，而是采用Picard迭代法求解。Picard迭代法具体步骤如下，在进行模拟时段的首次迭代计算时，根据给定的初始水头值来确定系数矩阵[A]中的元素值，并求解地下水运动差分方程组后得到首次迭代计算的水头值h¹，然后再根据h¹来更新系数矩阵[A]，再次求解地下水运动差分方程组得到第二次迭代计算的水头值h²。以上过程反复进行可求得h³…h^m，正常情况下，相邻两次迭代计算的水头值之差的绝对值(|h^m-h^m-1|)将逐渐变小，最终近乎相等。当在所有网格单元上本次迭代和上次迭代计算的水头值相差(|h^m-h^m-1|)小于事先规定的阈值(hclose)时，认为模拟时段的迭代计算收敛(模型收敛)，即地下水运动差分方程组求解完成，h^m即为模拟时段末的模拟水头值。

在仅启用全局松弛迭代机制时，判断收敛与否依据的是经过全局松弛后前后两次迭代计算的水头值相差是否小于hclose，而不是根据地下水运动差分方程前后两次迭代直接求解的水头值相差是否小于hclose判断；同理在同时启用全局松弛迭代机制和逐单元松弛迭代机制时，判断收敛与否依据的是经过全局松弛和逐单元松弛后前后两次迭代计算的水头值相差是否小于hclose。

S2、采用全局松弛迭代机制，提高地下水数值模拟模型的计算收敛性。

步骤S2为所有地下水数值模拟网格单元指定统一的全局松弛因子，对S1中由地下水运动差分方程组直接求解得到的某次迭代的水头值，根据全局松弛因子和上次迭代求得的水头值对其进行如下松弛，

其中：为地下水系统中任一网格单元(i,j,k)直接通过地下水运动差分方程组求解得到的某次迭代的水头值(L)；/>为该网格单元经过全局松弛后上次迭代计算的松弛水头值(L)；Damp为介于0.0-1.0之间的全局松弛因子，Damp＝1.0时表示不进行全局松弛，Damp值越小全局松弛程度越深；/>为该单元经过全局松弛后本次迭代计算的松弛水头值(L)。

S3、在步骤S2启用全局松弛迭代机制后，若模型仍然难以收敛或模拟结果的水量平衡误差较大，则继续采用逐单元松弛迭代机制，进一步智能提高地下水数值模拟模型的计算收敛性。采用逐单元松弛迭代机制，对每个网格单元的水头值都独立进行松弛，计算公式为，

其中：为任一网格单元(i,j,k)经过全局松弛后本次和前一次迭代计算的水头的相对变化值(L)；/>为逐单元松弛迭代机制下本次迭代计算时该单元的自适应松弛因子，为小于1的值，越小松弛程度越深；/>为逐单元松弛迭代机制下该单元本次迭代计算的松弛水头值(L)；

自适应松弛因子的计算如下：

其中：为上次迭代计算时网格单元(i,j,k)的自适应松弛因子值；α为自适应松弛因子的调减系数，一般取值为0.35-0.95；λ为自适应松弛因子的调增系数，一般取值为1-2；β为自适应松弛因子的单次调增值，一般取值为0.0-0.2。上式中自适应松弛因子/>将根据模拟是否出现数值震荡智能调减或调增，数值震荡(Oscilation)是指经过全局松弛后本次迭代计算的水头值变化方向与经过全局松弛后上次迭代计算的水头值变化方向相反，即：

其中：为网格单元(i,j,k)经过全局松弛后上一次迭代计算的水头变化值(L)。数值震荡一旦发生，表明地下水数值模拟模型的计算收敛性变差，因为水头误差无法按照一定的方向趋势性减少，而是反复时正时负，难以趋于稳定。

抑制数值模拟中的数值震荡，可以有效提高模型的计算收敛性。以上公式(5)表明，对任意一个网格单元(i,j,k)，每当发生数值震荡时，该单元的自适应松弛因子ω都将以α值为比例调减，相当于减缓迭代计算的速度从而抑制数值震荡；而每当连续无数值震荡的迭代次数累计超过NIT次时，该单元的自适应松弛因子ω将根据指定的λ和β值调增，相当于增加迭代计算的速度。NIT的作用也十分重要，该参数能够控制自适应松弛因子ω调增的频次，即每次调增了自适应松弛因子后，让模型先保持当前自适应松弛因子迭代计算几次，没有问题的话下次迭代可以继续调增，这样可以显著提高模型的计算收敛性。

本发明具有明显的智能化特点，使用本发明进行地下水数值模拟时每个网格单元上的水头松弛程度是不一样的，并且可以根据地下水运动差分方程组求解难易程度的变化，动态智能调整每个网格单元上的水头松弛程度，当数值震荡发生时(即方程迭代求解困难时)，能够立即响应调减自适应松弛因子ω，以增加收敛性；而当累计一定的次数没有发生数值震荡时(即方程迭代求解比较顺利时)，又能够合理地逐渐调增自适应松弛因子ω，以达到减少水量平衡误差的目的。该方法灵活高效，各技术参数物理意义明确。实例测试表明，该技术不仅能够大幅提高地下水数值计算的计算收敛性，同时也可以显著减少地下水的数值模拟误差。

全局松弛迭代机制和逐单元松弛迭代机制可以单独使用，也可配合使用，只要模型收敛且水量平衡误差很小，不同的松弛迭代参数组合产生的模拟结果只有微小差别。

本实施例地下水系统总面积23.7km²，地下水系统与湖泊、河流之间存在水力联系，地下水自东向西流出该地下水系统如图2(a)所示，地下水系统厚度介于32.7-88.8m之间如图2(b)所示。将该地下水系统空间上离散为72行、84列、8层网格单元，考虑降水入渗补给、定水头、湖泊、河流源汇项建立描述该地下水系统中水流运动的差分方程组(图2)，指定水头收敛阈值hclose＝0.00001m，要求模拟稳定状态下该地下水系统的水头分布情况。

如表1和表2所示，本实施例网格单元较多，在只启用全局松弛迭代机制的情况下，当Damp＝0.1时，地下水数值模拟模型虽然能够勉强收敛，但模拟结果的相对水量平衡误差超过20％，模拟结果的准确性很差，模拟并不成功，因此需要再启用逐单元松弛迭代机制以进一步智能提升模型的收敛性。经反复尝试，当Damp＝0.8,α＝0.7,λ＝1.5,β＝0.001,NIT＝5时模型能够收敛，并且相对水量平衡误差很小，表明模拟结果的准确性高，模拟是成功的，模拟水头如图4所示。

表1

表2

上述参数组合并不是唯一使模型模拟成功的参数组合，当Damp＝1，α＝0.7，λ＝1.5，β＝0.001，NIT＝5时，模型只启用了逐单元松弛迭代机制，未启用全局松弛迭代机制，但模型仍然是收敛的，并且模拟结果的相对水量平衡误差很小。由表2可知，在水量平衡误差很小的情况下，不同松弛迭代方案的模拟结果几乎完全一致，因此应用本发明时只需要找到一组合理的松弛迭代参数即可，方法灵活性高，表2中正值代表地下水系统的补给量，负值代表排泄量。

以上，仅为本申请较佳的具体实施方式，但本申请的保护范围并不局限于此，任何熟悉本技术领域的技术人员在本申请揭露的技术范围内，可轻易想到的变化或替换，都应涵盖在本申请的保护范围之内。因此，本申请的保护范围应该以权利要求的保护范围为准。

Claims (7)

1.一种智能提高地下水数值模拟收敛性的松弛方法，其特征在于，包括以下步骤：

构建地下水数值模拟模型；

采用全局松弛迭代机制，提高所述地下水数值模拟模型的计算收敛性，具体包括：

其中，为地下水系统中任一目标网格单元直接通过地下水运动差分方程组求解得到的某次迭代的水头值；/>为目标网格单元经过全局松弛后上次迭代计算的松弛水头值；Damp为全局松弛因子；/>为目标网格单元经过全局松弛后本次迭代计算的松弛水头值；

对所述地下水数值模拟模型的计算收敛性进行判断，若不满足模拟需求则继续采用逐单元松弛迭代机制，进一步智能提高所述地下水数值模拟模型的计算收敛性；对所述地下水数值模拟模型的计算收敛性进行判断的方法包括：

采用所述全局松弛迭代机制后，判断经过全局松弛后任一所述目标网格单元前后两次迭代计算的水头值相差是否都小于设定阈值，若经过若干次数迭代后任一所述目标网格单元前后两次迭代计算的水头值都小于所述设定阈值则所述地下水数值模拟模型的计算收敛，此时若模拟结果的水量平衡误差很小，表明模拟是成功的；若模拟结果的水量平衡误差较大或前后两次迭代计算的水头值相差始终无法在任一所述目标网格单元上都小于所述设定阈值则继续采用逐单元松弛迭代机制进一步智能提高所述地下水数值模拟模型的计算收敛性。

2.如权利要求1所述的智能提高地下水数值模拟收敛性的松弛方法，其特征在于，构建地下水数值模拟模型的方法包括：

基于地下水数值模拟网格单元系统结合所述地下水数值模拟网格单元系统的源汇项信息，根据水量平衡关系构建描述目标网格单元处地下水运动的差分方程，联立若干个所述网格单元的地下水运动差分方程构建所述地下水数值模拟模型。

3.如权利要求2所述的智能提高地下水数值模拟收敛性的松弛方法，其特征在于，所述地下水数值模拟网格单元系统由地下水系统空间离散后形成的具有指定行数、列数、层数的若干个网格单元组成。

4.如权利要求2所述的智能提高地下水数值模拟收敛性的松弛方法，其特征在于，所述水量平衡关系包括：在模拟时段内自周围六个相邻网格单元流入所述目标网格单元的地下水量与作用于所述目标网格单元之上的源汇项之和等于所述目标网格单元贮水量的变化。

5.如权利要求2所述的智能提高地下水数值模拟收敛性的松弛方法，其特征在于，根据水量平衡关系构建描述目标网格单元处地下水运动的差分方程的方法为：

其中，下标(i,j,k)为目标网格单元，下标(i,j-1,k)、(i,j+1,k)分别为目标网格单元沿行方向上左、右侧的相邻网格单元；(i+1,j,k)、(i-1,j,k)分别为目标网格单元沿列方向上前、后侧的相邻网格单元；(i,j,k-1)、(i,j,k+1)分别为目标网格单元沿层方向上上、下侧的相邻网格单元；CR、CC、CV分别为沿行、列、层方向上网格单元之间的水力传导系数；h为水头；上标m为本次迭代，上标m-1为上次迭代；P_i,j,k为作用于目标网格单元之上的与水头有关的源汇项相关系数；Q_i,j,k为作用于目标网格单元之上的流量源汇项；SCI_i,j,k为目标网格单元的贮水项；Δt为模拟时段的时间长度；Δh为模拟时段内的水头变化量。

6.如权利要求2所述的智能提高地下水数值模拟收敛性的松弛方法，其特征在于，联立若干个所述网格单元的地下水运动差分方程构建所述地下水数值模拟模型的方法为：

[A]{h}＝{q}

其中，[A]为系数矩阵；{h}为未知水头向量，代表各网格单元处的待求解水头；{q}为右端项向量，表示差分方程组中的常数项和已知项。

7.如权利要求1所述的智能提高地下水数值模拟收敛性的松弛方法，其特征在于，若不满足模拟需求则继续采用逐单元松弛迭代机制，进一步智能提高所述地下水数值模拟模型的计算收敛性的方法为：

其中，为任一目标网格单元经过全局松弛后本次和前一次迭代计算的水头的相对变化值；/>为逐单元松弛迭代机制下本次迭代计算时目标网格单元的自适应松弛因子；/>为逐单元松弛迭代机制下目标网格单元本次迭代计算的松弛水头值；/>为上次迭代计算时目标网格单元的自适应松弛因子值；α为自适应松弛因子的调减系数；λ为自适应松弛因子的调增系数；β为自适应松弛因子的单次调增值。