CN109960776A - 一种用于水力走时和水力信号衰减反演计算的改进算法 - Google Patents

一种用于水力走时和水力信号衰减反演计算的改进算法 Download PDF

Info

Publication number
CN109960776A
CN109960776A CN201910023300.4A CN201910023300A CN109960776A CN 109960776 A CN109960776 A CN 109960776A CN 201910023300 A CN201910023300 A CN 201910023300A CN 109960776 A CN109960776 A CN 109960776A
Authority
CN
China
Prior art keywords
walked
iteration
waterpower
hydraulic
algorithm
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Granted
Application number
CN201910023300.4A
Other languages
English (en)
Other versions
CN109960776B (zh
Inventor
胡睿
裘鹏翔
蒋弥
吴双红
刘泉
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Hohai University HHU
Original Assignee
Hohai University HHU
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Hohai University HHU filed Critical Hohai University HHU
Priority to CN201910023300.4A priority Critical patent/CN109960776B/zh
Publication of CN109960776A publication Critical patent/CN109960776A/zh
Application granted granted Critical
Publication of CN109960776B publication Critical patent/CN109960776B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
    • G06F17/10—Complex mathematical operations
    • G06F17/16—Matrix or vector computation, e.g. matrix-matrix or matrix-vector multiplication, matrix factorization

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Mathematical Physics (AREA)
  • Pure & Applied Mathematics (AREA)
  • Mathematical Analysis (AREA)
  • Mathematical Optimization (AREA)
  • Computational Mathematics (AREA)
  • Data Mining & Analysis (AREA)
  • Theoretical Computer Science (AREA)
  • Computing Systems (AREA)
  • Algebra (AREA)
  • Databases & Information Systems (AREA)
  • Software Systems (AREA)
  • General Engineering & Computer Science (AREA)
  • Management, Administration, Business Operations System, And Electronic Commerce (AREA)
  • Complex Calculations (AREA)

Abstract

本发明公开了一种用于水力走时和水力信号衰减反演计算的改进算法，Cimmino迭代通过调整迭代修正量的构造方法，来避免离散化后的D与走时之间的积分关系得到的矩阵方程的系数矩阵A的第j列所有元素为0，无法生成构造迭代修正量∆x^((k))的对角矩阵的问题；同时提出采用前50次迭代中残差收敛子列中残差最小的迭代作为反演结果的迭代次数，来解决问题（2）反演过程中没有具体准则用于确定终止迭代所需的迭代次数。本发明算法通过修改迭代规则，降低了对水力压力信号传播路径变化的敏感度，同时引入松弛因子，提升了算法的稳定性，并提出了新的标准迭代选取规则，提升了反演结果的准确度和展现含水层非均质性分布规律的能力。

Description

一种用于水力走时和水力信号衰减反演计算的改进算法
技术领域
本发明涉及一种用于水力走时和水力信号衰减反演计算的改进算法，属于水文地质与工 程地质、浅层地热及地下水环境技术领域。
背景技术
水力信号走时和衰减的积分都是采用将瞬态地下水流动方程转换为程函方程的形式进行 推导得出压力脉冲信号走时与介质水力扩散系数(D)以及水力信号衰减与介质的贮水系数 (Ss)的相关性积分方程的，而程函方程又可以使用射线追踪技术来求解，同时射线追踪允 许压力传播沿轨迹进行计算。因此，可以采用同步迭代重建算法(Simultaneous Iterative Reconstruction Technique，下称SIRT)将水力走时反演和衰减反演联合求解，从而得到信号 与激发点和观测点之间的含水层水力扩散系数(D)和贮水系数(Ss)在二维甚至是三维上的 非均质分布。
在使用迭代类反演算法的过程中，反演结果对应的迭代次数的选定主要基于残差的大小 以及以往的反演经验，缺乏足够的理论支撑。SIRT在实际运算中会遇到下列问题：
(1)将研究区域Ω进行网格划分如果网格内某一个格子Ωj没有任何一 个信号通过，即离散化后的D与走时之间的积分关系得到的矩阵方程的系数矩阵A的第j列 所有元素为0，无法生成构造迭代修正量Δx^((k))的对角矩阵；
(2)反演过程中没有具体准则用于确定终止迭代所需的迭代次数，即使在残差收敛之后， 反演结果仍然具有很大的差异。所以一个可行可靠，并且与模型无关的迭代次数判定准则对 于保证反演结果的正确性至关重要。
发明内容
本发明所要解决的技术问题是克服现有技术的缺陷，提供一种用于水力走时和水力信号 衰减反演计算的改进算法，基于Cimmino迭代的SIRT算法，以下简称为SIRT-Cimmino算法， Cimmino迭代通过调整迭代修正量的构造方法，来避免问题(1)；同时提出采用前50次迭代 中残差收敛子列中残差最小的迭代作为反演结果的迭代次数，来解决问题(2)。
为达到上述目的，本发明提供一种用于水力走时和水力信号衰减反演计算的改进算法，
包括以下步骤：
第一步，设置容忍值、初始值xinit和迭代终止准则；
第二步，开始第k次迭代，设初始值x(0)＝xinit，等号左侧上标为当前迭代数k，k为自然数， 使用射线追踪技术并遵循费马最小走时原理，从x(k)构造出矩阵A(k)，并求得b(k)＝A(k)x(k)， 此处b(k)表示第k次迭代后模拟路径上水力压力信号走时；A(k)是m行n列矩阵，记录压力信 号传播路径信息，x(k)是n维慢度向量，x(k)代表着水力扩散系数的分布；
第三步，计算残差Δb(k)＝b-b(k)，向量b记录了观测到的水力压力信号的走时如果满 足判定条件：残差小于容忍值或者超过设定的迭代次数，则算法终止并输出结果x(k)，否则 进入第四步；
第四步，在残差的基础上通过公式＃(1)构造迭代修正量Δx(k)，
Δx(k)＝λkA(k)TMΔb(k)＃(1)
此处，矩阵其中ai，i范围为0～m，代表A(k)第i个行 向量，A(k)T表示矩阵A(k)的转置即A(k)T中第i行j列元素等于A(k)中的第j行i列元素，m是A(k)的 行数，λk是松弛因子，变换得到公式＃(2)：
此处，Δb(k)T表示Δb(k)的转置，即Δb(k)T中第i行j列元素等于Δb(k)中的第j行i列元素；
第五步，使用迭代修正量设置第k+1次迭代的初始值：x(k+1)＝x(k)+Δx(k)，转入第二步，开始第k+1次迭代。
优先地，第一步包括，在缺少含水层扩散系数分布信息时，将初始分布设置为均质分布， 即xinit内每个分量大小相同。
优先地，第三步中，选择固定的迭代次数作为迭代终止准则。
优先地，判定条件为：残差超过固定的迭代次数，则算法终止并输出结果x(k)。
优先地，随着迭代次数k的更新，x(k)也随之更新，并导致矩阵A(k)发生相应的改变，所 以将λk设置为与迭代次数有关的活动参数能更大限度的保证算法的收敛性。
优先地，为了评估模拟路径上水力压力信号走时和观测到的水力压力信号走时之间的相 对误差，引入残差R(R代表Residual)，
其中是第k次迭代内第i个模拟路径上水力压力信号走时，ti是第i个观测到的水力压力 信号走时。
优先地，取前50次迭代中残差收敛子列中残差最小的迭代为标准迭代，对应的反演结果 为标准反演结果。
本发明所达到的有益效果：
本发明算法通过修改迭代规则，降低了对水力压力信号传播路径变化的敏感度，同时引 入松弛因子，提升了算法的稳定性，并提出了新的标准迭代选取规则，提升了反演结果的准 确度和展现含水层异质性分布规律的能力。
附图说明
图1是本发明的流程图；
图2是本发明的数值模拟几何模型示意图；
图3是本发明中模型水力扩散系数预设分布的示意图；
图4是现有技术中SIRT算法在8x6分辨率下的示意图；
图5是本发明SIRT-Cimmino算法在8x6分辨率下的示意图；
图6是现有技术中SIRT算法在8x8分辨率下的示意图；
图7是本发明SIRT-Cimmino算法在8x8分辨率下的示意图；
图8是现有技术中SIRT算法在12x12分辨率下的示意图；
图9是本发明SIRT-Cimmino算法在12x12分辨率下的示意图。
具体实施方式
下面结合附图对本发明作进一步描述。以下实施例仅用于更加清楚地说明本发明的技术 方案，而不能以此来限制本发明的保护范围。
考虑二维情况，假设目标区域Ω为介于激发点和接收点之间的矩形区域。通过网格将区 域Ω分解为n个小矩形格子(以下简称格子)，记为Ω1，…，Ωn，格子j内的扩散系数平均值设 为Dj。假设一共有m个压力信号，记第i个信号在格子Ωj的长度为sij，如果信号没有经过Ωj， 则sij＝0，因此走时反演的积分方程可以离散为：
式中：其中ti为压力响应从激发点到接收点的峰值走时，s为压力响应传播路径，D(s)为 在传播路径上水力学扩散系数函数，c为维度系数，当传播介质为2D时，c＝4，当传播介 质为3D时，c＝6。
上式可改写为矩阵形式
b＝Ax＃
此处b是m维向量，包含观测到的压力信号的走时；A是m行n列矩阵，记录压力信号传播 路径信息；x是n维慢度向量，代表着扩散系数的分布，其中：
本发明所提供的算法就是用于上述水力走时反演方法转化为通过解矩阵方程求未知量x。 上式中，b已知，A和x均未知。不过根据费马最小走时原理：信号沿用时最少的路径传播。 水力走时反演中扩散系数越大，慢度越小，用时越少，即信号倾向于通过高扩散系数区域。 因此在扩散系数分布(慢度向量x)已知的基础上可以通过射线追踪技术(Ray-Tracing)寻找 用时最少路径，进而推断出传播路径，从而构造矩阵A。
一种用于水力走时和水力信号衰减反演计算的改进算法，包括以下步骤：
第一步，设置容忍值、初始值xinit和迭代终止准则；
第二步，开始第k次迭代，设初始值x(0)＝xinit，等号左侧上标为当前迭代数k，k为自然数， 使用射线追踪技术并遵循费马最小走时原理，从x(k)构造出矩阵A(k)，并求得b(k)＝A(k)x(k)， 此处b(k)表示第k次迭代后模拟路径上水力压力信号走时；A(k)是m行n列矩阵，根据费马 最小走时原理：信号沿用时最少的路径传播。水力走时反演中扩散系数越大，慢度越小，用 时越少，即信号倾向于通过高扩散系数区域；因此在扩散系数分布(慢度向量x)已知的基础 上可以通过射线追踪技术(Ray-Tracing)寻找用时最少路径，进而推断出传播路径，从而构 造矩阵A，记录压力信号传播路径信息；x(k)是n维慢度向量，代表着水力扩散系数的分布， 在水力走时反演中D为水力扩散系数；
第三步，计算残差Δb(k)＝b-b(k)，向量b记录了观测到的水力压力信号的走时如果满 足判定条件：残差小于容忍值或者超过设定的迭代次数，则算法终止并输出结果x(k)，否则 进入第四步；
第四步，在残差的基础上通过公式＃(1)构造迭代修正量Δx(k)，
Δx(k)＝λkA(k)TMΔb(k)＃(1)
此处，矩阵其中ai，i范围为0～m，代表A(k)第i个行 向量，A(k)T表示矩阵A(k)的转置即A(k)T中第i行j列元素等于A(k)中的第j行i列元素，m是A(k)的 行数，λk是松弛因子，变换得到公式＃(2)：
此处，Δb(k)T表示Δb(k)的转置，即Δb(k)T中第i行j列元素等于Δb(k)中的第j行i列元素；
第五步，使用迭代修正量设置第k+1次迭代的初始值：x(k+1)＝x(k)+Δx(k)，转入第二步，开始第k+1次迭代。
进一步地，第一步包括，在缺少含水层扩散系数分布信息时，将初始分布设置为均质分 布，即xinit内每个分量大小相同。
进一步地，第三步中，选择固定的迭代次数作为迭代终止准则。
进一步地，判定条件为：残差超过固定的迭代次数，则算法终止并输出结果x(k)。
进一步地，随着迭代次数k的更新，x(k)也随之更新，并导致矩阵A(k)发生相应的改变， 所以将λk设置为与迭代次数有关的活动参数能更大限度的保证算法的收敛性。
进一步地，为了评估模拟路径上水力压力信号走时和观测到的水力压力信号走时之间的 相对误差，引入残差R(R代表Residual)，
其中是第k次迭代内第i个模拟路径上水力压力信号走时，ti是第i个观测到的水力压力 信号走时。
进一步地，取前50次迭代中残差收敛子列中残差最小的迭代为标准迭代，对应的反演结 果为标准反演结果。
为验证此算法，先建立一个水力学参数已知的正向模型并模拟抽水实验。然后，使用反演算 法反向计算出水力学参数，来和已知参数进行对比。
正向模拟部分采用有限元数值软件COMSOL Multiphysic在两个二维轴对称模型上模拟 抽水试验，如图2所示模型几何尺寸均为23.75m×32.8m，其中研究区域尺寸为4m×3.2m如 图2中颜色最深区域。如图3所示，研究区域设置一个模型用于模拟含有一条倾斜的高扩散 系数带的河流成因的砂砾混合含水层。
反向计算使用两种算法对模型走时数据进行反演。图4和图5显示了在8x6分辨率下两 种算法a)SIRT反演结果，b)SIRT-Cimmino反演结果的标准反演结果，图6和图7显示了在8x8分辨率下两种算法a)SIRT反演结果，b)SIRT-Cimmino反演结果的标准反演结果， 图8和图9显示了在12x12分辨率下两种算法a)SIRT反演结果，b)SIRT-Cimmino反演结 果的标准反演结果。
为了衡量反演结果和预设分布之间的差别，引入均方误差(Root Mean SquareError， RMSE)和相关系数(Correlation Coefficient)。均方误差是一个非负实数用于评估两个分布 在数值上的精确度，均方误差越接近于0表明两个分布在数值上越为接近。相关系数则是一 个介于[-1，1]之间的实数用于评估两个分布在结构上的相似度，越接近1表明两个分布相似 度越高。均方误差RMSE和相关系数Corr的定义如下：
其中是反演后第i个网格的扩散系数，是反演后的平均扩散系数 Di是第i个网格的预设扩散系数，是预设平均扩散系数
表1模型均方误差与相关系数
结果分析表明，在相同分辨率下，SIRT-Cimmino的结果能更好的体现中部高扩散系 数带的连通性。此外，表1记录了六个反演结果对应的均方误差与相关系数。两种算法的均 方误差接近，但是SIRT-Cimmino的相关系数更高，说明SIRT-Cimmino标准反演结果在结构 上更接近已知分布。
在水力走时反演中，压力信号的传播路径在每次扩散系数分布更新后会被重新设定，导 致矩阵A发生相应的改变，所以此时将λ设置为与迭代次数有关的活动参数能更大限度的保证 算法的收敛性。同时由于水力走时反演中矩阵A维数相对较小，所以相比信号路径判定，构造 松弛因子所需时间可以忽略不计。
以上所述仅是本发明的优选实施方式，应当指出，对于本技术领域的普通技术人员来说， 在不脱离本发明技术原理的前提下，还可以做出若干改进和变形，这些改进和变形也应视为 本发明的保护范围。

Claims (7)

1.一种用于水力走时和水力信号衰减反演计算的改进算法，其特征在于，包括以下步骤：
第一步，设置容忍值、初始值xinit和迭代终止准则；
第二步，开始第k次迭代，设初始值x(0)＝xinit，等号左侧上标为当前迭代数k，k为自然数，使用射线追踪技术并遵循费马最小走时原理，从x(k)构造出矩阵A(k)，并求得b(k)＝A(k)x(k)，此处b(k)表示第k次迭代后模拟路径上水力压力信号走时；A(k)是m行n列矩阵，记录压力信号传播路径信息，x(k)是n维慢度向量，x(k)代表着水力扩散系数的分布；
第三步，计算残差Δb(k)＝b-b(k)，向量b记录了观测到的水力压力信号走时，如果满足判定条件：残差小于容忍值或者超过设定的迭代次数，则算法终止并输出结果x(k)，否则进入第四步；
第四步，在残差的基础上通过公式＃(1)构造迭代修正量Δx(k)，
此处，矩阵其中ai，i范围为0～m，代表A(k)第i个行向量，表示矩阵A(k)的转置即中第i行j列元素等于A(k))中的第j行i列元素，m是A(k)的行数，λk是松弛因子，变换得到公式＃(2)：
此处，表示Δb(k)的转置，即中第i行j列元素等于Δb(k)中的第j行i列元素；
第五步，使用迭代修正量设置第k+1次迭代的初始值：x(k+1)＝x(k)+Δx(k)，转入第二步，开始第k+1次迭代。
2.根据权利要求1所述的一种用于水力走时和水力信号衰减反演计算的改进算法，其特征在于，第一步包括，在缺少含水层扩散系数分布信息时，将初始分布设置为均质分布，即xinit内每个分量大小相同。
3.根据权利要求1所述的一种用于水力走时和水力信号衰减反演计算的改进算法，其特征在于，第三步中，选择固定的迭代次数作为迭代终止准则。
4.根据权利要求3所述的一种用于水力走时和水力信号衰减反演计算的改进算法，其特征在于，判定条件为：残差超过固定的迭代次数，则算法终止并输出结果x(k)。
5.根据权利要求1所述的一种用于水力走时和水力信号衰减反演计算的改进算法，其特征在于，随着迭代次数k的更新，x(k)也随之更新，并导致矩阵A(k)发生相应的改变，所以将λk设置为与迭代次数有关的活动参数能更大限度的保证算法的收敛性。
6.根据权利要求1所述的一种用于水力走时和水力信号衰减反演计算的改进算法，其特征在于，为了评估模拟路径上水力压力信号走时和观测到的水力压力信号走时之间的相对误差，引入残差R(R代表Residual)，
其中是第k次迭代内第i个模拟路径上水力压力信号走时，ti是第i个观测到的水力压力信号走时。
7.根据权利要求1所述的一种用于水力走时和水力信号衰减反演计算的改进算法，其特征在于，取前50次迭代中残差收敛子列中残差最小的迭代为标准迭代，对应的反演结果为标准反演结果。
CN201910023300.4A 2019-01-10 2019-01-10 一种用于水力走时和水力信号衰减反演计算的改进方法 Active CN109960776B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201910023300.4A CN109960776B (zh) 2019-01-10 2019-01-10 一种用于水力走时和水力信号衰减反演计算的改进方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201910023300.4A CN109960776B (zh) 2019-01-10 2019-01-10 一种用于水力走时和水力信号衰减反演计算的改进方法

Publications (2)

Publication Number Publication Date
CN109960776A true CN109960776A (zh) 2019-07-02
CN109960776B CN109960776B (zh) 2023-06-16

Family

ID=67023485

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201910023300.4A Active CN109960776B (zh) 2019-01-10 2019-01-10 一种用于水力走时和水力信号衰减反演计算的改进方法

Country Status (1)

Country Link
CN (1) CN109960776B (zh)

Cited By (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN116842691A (zh) * 2023-05-24 2023-10-03 中国水利水电科学研究院 一种智能提高地下水数值模拟收敛性的松弛方法

Citations (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN103698810A (zh) * 2013-12-10 2014-04-02 山东科技大学 混合网最小走时射线追踪层析成像方法
US20140342370A1 (en) * 2012-01-13 2014-11-20 Advanced Genomic Technology, Llc Transgenic non-human animal model for accelerated aging and/or age-related symptom, and use thereof

Patent Citations (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US20140342370A1 (en) * 2012-01-13 2014-11-20 Advanced Genomic Technology, Llc Transgenic non-human animal model for accelerated aging and/or age-related symptom, and use thereof
CN103698810A (zh) * 2013-12-10 2014-04-02 山东科技大学 混合网最小走时射线追踪层析成像方法

Non-Patent Citations (1)

* Cited by examiner, † Cited by third party
Title
何林等: "水汽层析代数重构算法", 《测绘学报》 *

Cited By (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN116842691A (zh) * 2023-05-24 2023-10-03 中国水利水电科学研究院 一种智能提高地下水数值模拟收敛性的松弛方法
CN116842691B (zh) * 2023-05-24 2024-03-08 中国水利水电科学研究院 一种智能提高地下水数值模拟收敛性的松弛方法

Also Published As

Publication number Publication date
CN109960776B (zh) 2023-06-16

Similar Documents

Publication Publication Date Title
US6230101B1 (en) Simulation method and apparatus
Isebor et al. Generalized field-development optimization with derivative-free procedures
AU2011283192B2 (en) Methods and systems for machine-learning based simulation of flow
Chen et al. Streamline tracing and applications in embedded discrete fracture models
CN107590853A (zh) 一种城市建筑群震害高真实度展示方法
AU2011283193A1 (en) Methods and systems for machine-learning based simulation of flow
GB2458761A (en) Reservoir simulation using dynamic and static region values
AU2011283190A1 (en) Methods and systems for machine-learning based simulation of flow
CN104504754B (zh) 一种油气储层多点统计建模的方法及装置
CA2803067A1 (en) Methods and systems for machine-learning based simulation of flow
CN105093278B (zh) 基于激发主能量优化算法的全波形反演梯度算子提取方法
US9589081B2 (en) Karstification simulation
Chi et al. Back analysis of the permeability coefficient of a high core rockfill dam based on a RBF neural network optimized using the PSO algorithm
AU2017382546A1 (en) Subsurface modeler workflow and tool
CN107894618A (zh) 一种基于模型平滑算法的全波形反演梯度预处理方法
CN111123374A (zh) 一种基于匹配滤波的探地雷达全波形反演方法
CN106503407A (zh) 存在部分连通断层的线性水侵油藏的试井分析方法及装置
CN109655890A (zh) 一种深度域浅中深层联合层析反演速度建模方法及系统
CN109960776A (zh) 一种用于水力走时和水力信号衰减反演计算的改进算法
CN108460838A (zh) 三维可视化技术与数值模拟技术融合的实现方法与系统
Gong et al. Application of multi-level and high-resolution fracture modeling in field-scale reservoir simulation study
CN114547958A (zh) 基于深度神经网络的井震结合裂缝预测方法及装置
CN111751886A (zh) 一种基于微地震监测数据的页岩气藏裂缝建模方法
Mohajerani et al. An efficient computational model for simulating stress-dependent flow in three-dimensional discrete fracture networks
CN108765570A (zh) 一种三维可视化与数值模拟融合的实现方法与系统

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant