CN114861519B - 复杂地质条件下初始地应力场加速优化反演方法 - Google Patents

复杂地质条件下初始地应力场加速优化反演方法 Download PDF

Info

Publication number
CN114861519B
CN114861519B CN202210216125.2A CN202210216125A CN114861519B CN 114861519 B CN114861519 B CN 114861519B CN 202210216125 A CN202210216125 A CN 202210216125A CN 114861519 B CN114861519 B CN 114861519B
Authority
CN
China
Prior art keywords
stress
value
boundary condition
calculation
inversion
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.)
Active
Application number
CN202210216125.2A
Other languages
English (en)
Other versions
CN114861519A (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.)
Chengdu Univeristy of Technology
Original Assignee
Chengdu Univeristy of Technology
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 Chengdu Univeristy of Technology filed Critical Chengdu Univeristy of Technology
Priority to CN202210216125.2A priority Critical patent/CN114861519B/zh
Publication of CN114861519A publication Critical patent/CN114861519A/zh
Application granted granted Critical
Publication of CN114861519B publication Critical patent/CN114861519B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F30/00—Computer-aided design [CAD]
    • G06F30/20—Design optimisation, verification or simulation
    • G06F30/27—Design optimisation, verification or simulation using machine learning, e.g. artificial intelligence, neural networks, support vector machines [SVM] or training a model
    • 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/18—Complex mathematical operations for evaluating statistical data, e.g. average values, frequency distributions, probability functions, regression analysis
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F18/00—Pattern recognition
    • G06F18/20—Analysing
    • G06F18/24—Classification techniques
    • G06F18/241—Classification techniques relating to the classification model, e.g. parametric or non-parametric approaches
    • G06F18/2413—Classification techniques relating to the classification model, e.g. parametric or non-parametric approaches based on distances to training or reference patterns
    • G06F18/24133—Distances to prototypes
    • G06F18/24137—Distances to cluster centroïds
    • G06F18/2414—Smoothing the distance, e.g. radial basis function networks [RBFN]
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F30/00—Computer-aided design [CAD]
    • G06F30/20—Design optimisation, verification or simulation
    • G06F30/23—Design optimisation, verification or simulation using finite element methods [FEM] or finite difference methods [FDM]
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06N—COMPUTING ARRANGEMENTS BASED ON SPECIFIC COMPUTATIONAL MODELS
    • G06N3/00—Computing arrangements based on biological models
    • G06N3/02—Neural networks
    • G06N3/04—Architecture, e.g. interconnection topology
    • G06N3/048—Activation functions
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06N—COMPUTING ARRANGEMENTS BASED ON SPECIFIC COMPUTATIONAL MODELS
    • G06N3/00—Computing arrangements based on biological models
    • G06N3/02—Neural networks
    • G06N3/08—Learning methods
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F2111/00—Details relating to CAD techniques
    • G06F2111/04—Constraint-based CAD
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F2111/00—Details relating to CAD techniques
    • G06F2111/06—Multi-objective optimisation, e.g. Pareto optimisation using simulated annealing [SA], ant colony algorithms or genetic algorithms [GA]
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F2111/00—Details relating to CAD techniques
    • G06F2111/10—Numerical modelling
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F2119/00—Details relating to the type or aim of the analysis or the optimisation
    • G06F2119/14—Force analysis or force optimisation, e.g. static or dynamic forces

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • General Physics & Mathematics (AREA)
  • Data Mining & Analysis (AREA)
  • General Engineering & Computer Science (AREA)
  • Evolutionary Computation (AREA)
  • Mathematical Physics (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Artificial Intelligence (AREA)
  • Software Systems (AREA)
  • Health & Medical Sciences (AREA)
  • Pure & Applied Mathematics (AREA)
  • Bioinformatics & Cheminformatics (AREA)
  • Bioinformatics & Computational Biology (AREA)
  • Computational Mathematics (AREA)
  • Computer Vision & Pattern Recognition (AREA)
  • Evolutionary Biology (AREA)
  • Mathematical Analysis (AREA)
  • Mathematical Optimization (AREA)
  • Computing Systems (AREA)
  • Molecular Biology (AREA)
  • General Health & Medical Sciences (AREA)
  • Computational Linguistics (AREA)
  • Computer Hardware Design (AREA)
  • Biophysics (AREA)
  • Geometry (AREA)
  • Biomedical Technology (AREA)
  • Databases & Information Systems (AREA)
  • Algebra (AREA)
  • Medical Informatics (AREA)
  • Probability & Statistics with Applications (AREA)
  • Operations Research (AREA)
  • Management, Administration, Business Operations System, And Electronic Commerce (AREA)

Abstract

本发明公开了复杂地质条件下初始地应力场加速优化反演方法，包括步骤：建立用于应力场反演计算的三维地质模型；确定构造应力边界的施加方式和范围值，再通过数值计算求得各应力边界条件下的测点应力计算值；将应力边界条件作为网络输入，以测点应力计算值作为网络输出，建立代理模型的学习样本，并将该学习样本带入代理模型中进行训练；基于训练好的代理模型，使用遗传算法在训练好的样本范围内进行快速寻优，得到最优应力边界条件；以最优应力边界条件反演计算得到区域范围内的初始地应力场。通过本发明提供的反演方法，可快速高效的实现地应力场的反演，还能保证计算结果和实测值之间较高的精度。

Description

复杂地质条件下初始地应力场加速优化反演方法
技术领域
本发明涉及地质工程技术领域，尤其涉及一种复杂地质条件下初始地应力场加速优化反演方法。
背景技术
初始地应力场是地下工程设计、地质稳定性评价及岩土施工等方面的重要指标，实测地应力是提供工程区地应力状态的一种直接方法，但由于测点数量限制，实测值只能代表测点附近一定范围内的应力情况，很难反映区域整体的应力场。针对此问题，通常是在一定实测地应力数据的基础上，结合工程区地质条件，利用数值分析方法开展区域应力场反演，以获得准确且范围较广的地应力分布特征。通过三维数值模拟软件进行含构造地应力反演是其中一种有效的途径，但目前现场实测地应力普遍测点较少，同时利用三维数值模拟软件对地应力进行反演时，存在反演计算量较大，计算精度不足等问题。此外，若为了平衡计算速度与计算精度，所建立的三维模型往往尺寸较小，难以达到工程大尺度，实际计算效果与精度会受到局限。
目前，用于应力场反演的数据分析方法还包括边界荷载调整法、应力函数法、多元线性回归分析法、神经网络反分析法、遗传算法等。然而，目前的分析方法均不同程度的存在反演计算速度较慢或者精度不足等问题，难以进行较好的推广应用。
发明内容
本发明的目的在于解决目前的初始应力场反演方法存在反演效果不理想的问题，提供一种可应用于长大深埋隧道地应力场和复杂地质地应力场，可提升复杂地质条件下初始应力场反演计算的精度和效率的复杂地质条件下初始地应力场加速优化反演方法。
本发明的目的是这样实现的，提供一种复杂地质条件下初始地应力场加速优化反演方法，包括步骤：
S1、建立用于工程区应力场反演计算的三维地质模型；
S2、结合工程区的区域构造特征，确定应力边界条件参数的数值范围，并在该数值范围内进行参数组合得到多组应力边界条件，施加应力边界条件，通过数值计算求得各应力边界条件下的测点地应力计算值；
S3、将应力边界条件作为网络输入，以测点应力计算值作为网络输出，建立代理模型的学习样本，并将该学习样本带入代理模型中进行训练；
S4、基于步骤S3训练好的代理模型，使用遗传算法在样本范围内进行快速寻优，得到最优应力边界条件；
S5、以最优应力边界条件反演计算得到工程区内初始地应力场。
优选的，步骤S2中，确定应力边界条件的范围包括：
S21、将水平构造应力设置在以下范围：
Figure BDA0003534782470000021
Figure BDA0003534782470000022
其中бT、бt分别为水平最大构造应力及水平最小构造应力，бHбh分别为实测水平最大主应力及水平最小主应力，γ为岩石容重，H为埋深，μ为泊松比；
S22、设定实测水平最大主应力、水平最小主应力随埋深为以下线性分布：
σH＝A*H+B，σh＝C*H+D；
其中，A、B、C、D为线性回归常数，将上述水平应力回归式代入水平构造应力不等式中，得到水平构造应力取值范围为：
Figure BDA0003534782470000031
Figure BDA0003534782470000032
S23、设定水平构造应力与埋深的对应的线性函数为：
σT＝a*H+b，σt＝c*H+d；
结合水平构造应力的取值范围，确定a、b、c、d的取值范围为：
Figure BDA0003534782470000033
B-1≤b≤B+1
Figure BDA0003534782470000034
D-1≤d≤D+1
以a、b、c、d作为应力边界条件参数，获得应力边界条件参数的数值范围。
优选的，步骤S3中，采用径向基神经网络作为代理模型，将应力边界条件和测点高程作为网络输入，测点应力计算值作为网络输出，每个测点的应力计算值包括最大主应力值、最小地应力值及垂直主应力值。
优选的，步骤S4包括：
S41、确定遗传算法的参数，并在应力边界条件范围内随机产生多组应力边界条件作为父代群体；
S42、利用训练好的代理模型对每组应力边界条件进行预测，得到每组应力边界条件所对应的应力值；
S43、基于步骤S42中得到应力值对每组应力边界条件进行适应度评价，适应度均不合理，则迭代产生新的N组应力边界条件并进行适应度评价，直到产生适应度合理的应力边界条件；
S44、在适应度合理的应力边界条件中选择适应度最高的应力边界条件，并通过数值计算得到该应力边界条件下的应力值，将所计算出的应力值与实际测量应力值进行对比，判断两者之间的误差是否满足设定的要求，满足则该应力边界条件即为最优应力边界条件。
优选的，当步骤S44中计算出的应力边界条件下的应力值与实际测量应力值之间的误差不满足设定的要求时，将该应力边界条件及对应的应力值作为一组样本加入到代理模型的训练样本中，重新训练代理模型，同时以该应力边界条件迭代产生新的N组应力边界条件，重复步骤S42到S44。
优选的，步骤S43中，当适应度均不合理时，从所评价的应力边界条件中选择适应度最高的应力边界条件进行迭代。
优选的，应力边界条件的适应度的评价方法为：使用代理模型预测一组应力边界条件下不同高程对应的水平最大主应力、水平最小主应力、垂直应力，并将这些应力值与实际测量值作差，以差值绝对值之和的负数作为该组应力边界条件的适应度。
本发明的显著进步性至少体现在：
本发明提供的复杂地质条件下初始地应力场加速优化反演方法，在整个数值模拟及计算方面的加速效果明显，提高了计算效率。方法实施中所涉及的模型不需非常精细化建模的前提下，还能保证计算结果和实测值之间较高的精度，可较好的应用于长大深埋隧道地应力场和复杂地质地应力场，有效提升复杂地质条件下初始应力场反演计算的精度和效率。
附图说明
图1为本发明一种实施例的初始应力场加速优化反演方法的流程图；
图2为本发明另一实施例的初始应力场加速优化反演方法的流程图；
图3为本发明实施例中第一构造应力计算式下的水平最大构造应力与埋深的线性拟合图；
图4为本发明实施例中第一构造应力计算式下的水平最小构造应力与埋深的线性拟合图；
图5为本发明实施例中第二构造应力计算式下的水平最大构造应力与埋深的线性拟合图；
图6为本发明实施例中第二构造应力计算式下的水平最小构造应力与埋深的线性拟合图；
图7为实施例反演方法获得的最大主应力值与实测值得对比图；
图8为实施例反演方法获得的最小主应力值与实测值得对比图；
图9为实施例反演方法获得的垂直主应力值与实测值得对比图。
具体实施方式
下面结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整的描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部实施例。基于本发明的实施例，本领域技术人员在没有做出创造性劳动的前提下所获得的所有其它实施方式，都属于本发明保护的范围。
本发明实施例提供的复杂地质条件下初始地应力场加速优化反演方法，包括步骤：
S1、建立用于工程区应力场反演计算的三维地质模型；
S2、结合工程区的区域构造特征，确定应力边界条件参数的数值范围，并在该数值范围内进行参数组合得到多组应力边界条件，施加应力边界条件，通过数值计算求得各应力边界条件下的测点地应力计算值；
S3、将应力边界条件作为网络输入，以测点应力计算值作为网络输出，建立代理模型的学习样本，并将该学习样本带入代理模型中进行训练；
S4、基于步骤S3训练好的代理模型，使用遗传算法在样本范围内进行快速寻优，得到最优应力边界条件；
S5、以最优应力边界条件反演计算得到工程区内初始地应力场。
应该说明的是，上述方案中，将应力边界条件确定在合适的范围内，可避免样本构造的盲目性，较为明显的提升样本训练的有效性，加快最优应力边界条件的确定。本发明实施例所建立的代理模型以应力边界条件作为网络输入，以测点应力计算值作为网络输出，通过学习样本带入代理模型中进行训练，从而建立输入的造应力边界条件与输出的应力值之间的非线性映射关系，其所建立的是对原有数值模型的一种局部代理模型，直接利用局部代理模型来替代模型的数值计算，则不再对数值模型(三维地质模型)本身求解。不仅能够大幅度地减少计算负荷，还能节省大量时间。数值计算中，可利用有限差分FLAC3D软件进行数值计算，得到不同应力边界条件下的地应力值。进一步的，使用遗传算法在样本范围内进行寻优，可快速得到最优应力边界条件，提高反演效率。综上，本实施例中，数值模拟及计算方面的时间效率高，体现了该方法的加速计算优势；相比相关现有技术，该方法的显著进步还体现在对于三维地质模型不需精细化建模(施加自重应力的基础上，只需训练和构建边界条件即可)的前提下，还能保证计算结果和实测值之间较高的精度。
进一步应该说明的是，在建立三维地质模型后，可先对三维地质模型进行初步的运算和检验，具体的，将模型通过单元及网格的划分后，导入基于有限差分原理的FLAC3D软件中进行数值计算。编写自重应力命令流实现初始条件设置和边界条件约束，在X与Y面边界施加固定约束，在Z轴方向施加自重应力，来检验复杂地质体自重应力状态及其变化规律，也便于检查模型内部的好坏。
由于影响岩体初始地应力场的因素极多，导致地应力场的分布规律非常复杂，本发明实施例中以隧址区地应力场反演为例，可以从考虑地形地貌和构造主应力的方向和量值方面来对三维地质模型施加应力边界条件，根据对实测地应力的描述，隧址区地应力场以水平构造应力为主。本发明实施例中通过对X向与Y向施加构造作用来对地应力场进行模拟，在Z轴竖直施加自重应力条件时，模型的其余面均施加位移约束；在沿着X轴或Y轴方向施加构造应力条件时，对竖直垂直轴的侧面施加构造应力梯度，其余侧面施加位移约束。
可以理解的是，在明确主应力的构成和自重应力容易确定的情况下，如何确定合理的构造应力尤其关键。作为一种优选的实施例方案，步骤S2中，确定应力边界条件的范围包括：
S21、将水平构造应力确定在合适的范围
本实施例中，通过第一构造应力计算式确定水平构造应力的下限，即水平构造应力中的水平最大构造应力的下限及水平最小构造应力的下限的分别为：
σT1＝σH-γ*H，σt1＝σh-γ*H；
通过第二构造应力计算式确定水平构造应力的上限，即水平最大构造应力的上限及水平最小构造应力的上限分别为：
Figure BDA0003534782470000081
由此，将本实施例中的水平构造应力确定在以下范围
Figure BDA0003534782470000082
Figure BDA0003534782470000083
其中：бT、бt分别为水平最大构造应力及水平最小构造应力，бH、бh分别为实测水平最大主应力及水平最小主应力，γ为岩石容重，H为埋深，μ为泊松比；
S22、设定实测水平最大主应力、水平最小主应力随埋深为以下线性分布：
σH＝A*H+B，σh＝C*H+D；
其中，A、B、C、D为线性回归常数，将上述水平应力回归式代入水平构造应力不等式中，得到水平构造应力取值范围为：
Figure BDA0003534782470000091
Figure BDA0003534782470000092
S23、设定水平构造应力与埋深的对应的线性函数为：
Figure BDA0003534782470000093
Figure BDA0003534782470000094
根据上式计算确定a、b、c、d的取值范围，并将该四个数值作为应力边界条件的参数，确定应力边界条件的范围。
以上实施例中，将应力边界条件设置在特定的范围内，充分考虑了实际构造应力的复杂性，能够较为准确的反应出复杂地质体的地质构造应力，从而有助于通过该构造应力条件获得较为准确合理的初始地应力分布特征。此外，通过实施例中确定的应力边界条件的范围可在确保准确性的情况下进一步提升确定最优应力边界条件的效率。
作为一种优选的实施例方案，为更好的实现代理模型对数值计算的替代，步骤S3中，本实施例选择采用径向基神经网络作为代理模型。进一步的，将应力边界条件和测点高程作为网络输入，测点应力计算值作为网络输出，每个测点的应力计算值包括最大主应力值、最小地应力值及垂直主应力值。
作为一种优选的实施例方案，步骤S4包括：
S41、确定遗传算法的参数，并在应力边界条件范围内随机产生多组应力边界条件作为父代群体；
S42、利用训练好的代理模型对每组应力边界条件进行预测，得到每组应力边界条件所对应的应力值；
S43、基于步骤S42中得到应力值对每组应力边界条件进行适应度评价，适应度均不合理，则迭代产生新的N组应力边界条件并进行适应度评价，直到产生适应度合理的应力边界条件；
S44、在适应度合理的应力边界条件中选择适应度最高的应力边界条件，并通过数值计算得到该应力边界条件下的应力值，将所计算出的应力值与实际测量应力值进行对比，判断两者之间的误差是否满足设定的要求，满足则该应力边界条件即为最优应力边界条件；不满足，则将该应力边界条件及对应的应力值作为一组样本加入到代理模型的训练样本中，重新训练代理模型，同时以该应力边界条件迭代产生新的N组应力边界条件，重复步骤S42到S44。
参阅图1所示为基于本实施例的初始应力场加速优化反演方法的流程图，可以理解的是，其中的建立三维模型即为步骤S1中的建立三维地质模型的过程，其中的构造边界条件的组合以及地应力值的计算对应上述步骤S2，其中的建立替代模型(代理模型)学习样本和训练代理模型对应为步骤S3，在此不作一一赘述。
作为一种优选的实施例方案，所述步骤S2还包括：通过数值计算求得各应力边界条件下的测点应力计算值后，寻找出相对最优应力边界条件，并将该相对最优应力边界条件对应的应力值与实测应力进行对比，判断两者之间的误差是否满足设定的要求，满足则执行步骤S5，不满足则执行步骤S3。进一步优选的，在代理模型(替代模型)完成样本训练后，进行代理模型是否达到收敛的判断，收敛，则样本训练学习结束，找出其中最优的应力边界条件并进行反演计算。参阅图2所示为基于本实施例的初始应力场加速优化反演方法的流程图。可以理解的是，本实施例方案中，步骤S2中通过数值计算求得各应力边界条件下的测点应力计算值，该步骤相当于传统的应力试算，根据该试算寻找出的相对最优应力边界条件，其对应的应力值与实测应力值之间的误差在实际测试中大概率是无法否满足设定要求的。
作为优选的，步骤S43中，当适应度均不合理时，从所评价的应力边界条件中选择适应度最高的应力边界条件进行迭代。迭代产生的新个体(应力边界条件)数量N可依据实际需求进行设置，综合考虑应力边界条件寻优的效率及合理性，本实施例中建议N取值为50-100之间的整数。
应该说明的是，以上方案中，基于训练好的代理模型，并对应力边界条件进行适应度的评价，可快速有效的筛选出合理的边界条件，而当应力边界条件均不合理时，则选择适应度最高的应力边界条件进行迭代，有助于快速得到最优的应力边界条件。
作为一种优选的实施例方案，应力边界条件的适应度的评价方法为：使用代理模型预测一组应力边界条件下不同高程对应的水平最大主应力、水平最小主应力、垂直应力，并将这些应力值与实际测量值作差，以各差值绝对值之和的负数作为该组应力边界条件的适应度，即各差值绝对值之和越大，适应度越低。可以理解的是，适应度是否合理，可按照实际需求设置适应度判定为合理的标准范围，适应度在标准范围内，则认定为合理的适应，在标准范围内则判定为不合理的适应度。采用本实施例中的适应度评价方法，可实现对应力边界条件快速有效的评价。
下面，针对隧址区情况给出一种具体的实际测试应用例。
(1)建立用于地应力场计算的隧址区三维地质模型，并对模型进行初步的运算和检验；
(2)确定应力边界条件的范围，并在该范围内分为多组边界条件组合，通过数值计算求得各应力边界条件组合下的测点应力计算值；
具体的，将隧址区地应力实测数据(包括各测点的埋深以及实测水平最大主应力、水平最小主应力)代入第一构造应力计算式中，计算得到各测点的水平构造应力计算值，计算结果如表1所示：
表1
Figure BDA0003534782470000121
Figure BDA0003534782470000131
表1中σv为根据上覆岩石埋深计算的垂向主应力，计算中的岩石容重取26.5kN/m3。通过上述方法得到的水平最大构造应力及水平最小构造应力随埋深的变化趋势分别见图3和图4。
经过整理可得到如下回归关系式(1)：
σT1＝0.007H-1.959，σt1＝-0.002H-0.844。
将隧址区地应力实测数据代入第二构造应力计算式中，计算得到各测点的水平构造应力计算值，计算结果如表2所示：
表2
Figure BDA0003534782470000132
Figure BDA0003534782470000141
表2中σv为根据上覆岩石埋深计算的垂向主应力，计算中的岩石容重取26.5kN/m3。通过上述方法得到的水平最大构造应力及水平最小构造应力随埋深的变化趋势见图5和图6。
通过整理可得到如下回归关系式(2)：
σT2＝0.024H-2.034，σt2＝0.015H-0.918
根据以上回归关系式(1)和(2)，确定a、b、c、d的取值范围为：
0.007≤a≤0.024，-2.034≤b≤-1.959
-0.002≤c≤0.015，-0.918≤d≤-0.844。
为了避免参数范围的局限性，可将上述得到的参数范围适当扩大，a的下限值增至0.005、上限值增至0.025，b的上限值增至-1，d的下限值增至-1、上限值增至0。最终确定a、b、c、d的取值范围分别为[0.005，0.025]、[-2，-1]、[0，0.015]、[-1，0]，将确定的取值范围进行均匀设计，将反演参数划分为50组，均匀分组见表3，作为本测试例反演方法的应力边界条件范围，进一步的，计算50组应力边界条件对应的应力值。
表3
样本序号 a b c d
1 0.005 -2 0 -1
2 0.005 -1.75 0 -0.75
3 0.01 -2 0.004 -1
4 0.01 -1.75 0.004 -0.75
5 0.015 -2 0.008 -1
6 0.015 -1.75 0.008 -0.75
7 0.02 -2 0.012 -1
8 0.02 -1.75 0.012 -0.75
9 0.025 -1.75 0.015 -0.75
10 0.025 -1.5 0.015 -0.5
11 0.005 -2 0.004 -1
12 0.005 -1.75 0.004 -0.75
13 0.01 -2 0.008 -1
14 0.01 -1.75 0.008 -0.75
15 0.015 -2 0.012 -1
16 0.015 -1.75 0.012 -0.75
17 0.02 -2 0.015 -0.75
18 0.02 -1.75 0.015 -0.5
19 0.025 -1.75 0 -1
20 0.025 -1.5 0 -0.75
21 0.005 -2 0.008 -1
(3)将应力边界条件作为网络输入，以测点应力计算值作为网络输出，建立代理模型的学习样本，并将该学习样本带入代理模型中进行训练；
具体的，包括(a)构建训练样本：将表3中的50组应力边界条件组合分别作为网络输入，以通过数值计算得到的50组应力边界条件下不同埋深对应的水平最大主应力、水平最小主应力、垂直应力作为径向基神经网络的输出，建立测点应力与边界条件非线性映射关系的网络模型。(b)模型训练：使用径向基函数网络作为替代模型，首先对输入进行归一化，使用每一类输入的最大值和最小值将输入归一化到0至1之间，对目标也归一化到0至1之间。整个训练过程都一样，输入层节点数为4，输出层节点数为3，使用高斯函数作为径向基函数(激活函数)，高斯函数扩展参数使用固定值为0.1，使用K均值聚类进行无监督学习确定聚类中心，聚类中心个数为15，训练时计算各数据点的激活值，构建激活值矩阵；通过最小二乘法来优化径向基网络训练，进而转化为矩阵方程求解得到网络权重，由此就可达到误差的极小值，再由权重向量和径向基函数值的点积，即可得到网络的输出。
(4)基于训练好的代理模型，使用遗传算法在样本范围内进行快速寻优，得到最优应力边界条件；
具体的，使用遗传算法进行寻优，适应性评估函数为负的拟合误差，基因个数为4，即4个应力边界条件，取值范围为(0，1)，种群大小为100，变异率为0.1，选择率为0.6，代数为50，选择轮盘赌策略。遗传算法每次迭代产生50组应力边界条件(个体)，使用替代模型估计每一组边界条件下26个埋深对应的3种应力值，并将这些应力与实际测量值作差，差值绝对值之和的负数为该边界条件的适应度，遗传算法每次基于适应度高的个体进行迭代，经过50次遗传算法迭代，最高适应度的边界条件即为遗传算法的输出。
(5)通过数值计算得到高适应度边界条件的应力值，与实际测量值对比，判断是否满足要求，满足则该边界条件即为目标边界条件，否则将该边界条件及对应的应力值作为一组数据加入到替代模型的训练数据中，重复(3)、(4)、(5)步。最终得到的最优应力边界条件为：a＝0.0186，b＝-1.581，c＝0.013，d＝-0.076。
对以上隧址区初始地应力场反演结果分析，将26个实测点的地应力实测值与计算值对比见表4：
表4
Figure BDA0003534782470000171
Figure BDA0003534782470000181
由上表可得：最小误差为1％，当相对误差最大，为16％时，可见实测值与计算值的量值小于2MPa，为1.50MPa，此处体现在最小主应力；当误差量值最大时为3.79MPa，相对误差为-12％，此处体现在最大主应力。再将数据以直观的方式展示并比较相对误差见图7-9。
由图7-9可看出：最大主应力值的最大误差为12％，但只有少数4个点的误差值超过了10％，其余测点的误差值均在10％的误差范围内；最小主应力值的相对误差最大达到16％，但只有少数4个点的误差超过了10％，误差值总体在10％左右；垂直主应力值得相对误差最大达到11％，但只有少数2个点的误差值超过了10％，其余误差值都均控制在10％以内。因此可见，反演计算获得的三个主应力值与实测值较为接近，说明本发明实施例的反演法的反演结果较为合理。
地应力反演精度和效率的评价
在评价该反演方法的合理性、有效性的时候，主要是要对比所得的反演计算结果是否与实测值相接近，以上的反演结果已表明了本发明实施例方法的反演精确度较高，反演结果较为合理。此外，反演的高效性也是实际应用中较为关注的评价指标，下面，将本实施例提出的反演方法与现有的较好的智能反演方法(即以仅基于RBF模型的反演方法作为对比方法)在实现过程的时间效率和精确度上作同等对比。采用典型的评价函数再对两种方法得出的结果进行对比。
用多种误差函数分别进行计算和对比，结果见表5，再将它们进行时间效率上的对比，结果见表6：
表5
Figure BDA0003534782470000191
表6
Figure BDA0003534782470000192
由以上对两种反演方法在精度与时间效率上的对比可以看出，在平方误差、绝对误差和相对误差的对比下，实施例方法反演出各主应力的各种误差均明显小于对比方法的反演误差值，可见，在反演精度方面，本发明实施例的反演方法高于对比的反演方法。在时间效率的对比下，实施例方法训练一次的时间(即训练完选定组数的样本所用的时间，表6中为训练50组样本的时间)为0.057s，而对比方法训练一次时间是0.092s，数值计算时间(得到最终计算结果的所有计算总耗时)：实施例方法为90000s，对比方法为180000s，本发明实施例方法的耗时仅为对比方法的一半。由上述分析可得出：本发明实施例提出的方法在拟合精度和时间效率上均明显优于对比方法。
值得说明的是，对比方法(RBF反演法)作为目前较为流行和高效的反演方法，基于RBF神经网络进行，需要大量的训练样本，而训练样本需通过大量的数值计算而得，因此，依然存在数值计算量较大的问题。进一步的，在本发明实施例方法中，只需选取少量样本训练代理模型，使得代理模型形成应力边界条件与地应力值之间的映射关系即可，从而使得通过数值计算产生训练样本的总时间大幅减少；本发明实施例方法的应力边界条件寻优并不依赖于训练样本的数量，而是通过遗传算法迭代若干组应力边界条件进行适应度评价，可快速筛选掉适应度不合理的应力边界条件，每次基于筛选后获得的最高适应度的应力边界条件进行最后的数值计算，数值计算的耗时较少，从而实现通过少量训练样本就能够寻到最优的应力边界条件。还应说明的是，在以上对比测试中，虽然对比方法单次训练样本的时间和本发明实施例方法接近，但是现有RBF方法的神经网络输入(应力值)输出(边界条件)与本发明实施例方法的输入(边界条件)输出(应力值)是相反的，其是通过不断训练大量的样本，当神经网络训练收敛时进行边界条件的寻优，其所需的训练样本数量远大于本发明实施例方法所需的样本数量，因此其数值计算的总时间较长。此外，因为无法避免RBF神经网络陷入局部极小值的固有弊端，RBF反演后的计算结果不能保证为全局最优解。本发明实施例中使用代理模型结合遗传算法既能解决神经网络算法在单独使用中存在的缺陷，又能保证计算结果有较高的精度和效率，能用更少的样本数量，在更短的时间内寻找到最优的边界组合条件。
最后还应该说明的是，尽管已经示出和描述了本发明的实施例，对于本领域的普通技术人员而言，可以理解在不脱离本发明的原理和精神的情况下可以对这些实施例进行多种变化、修改、替换和变型，本发明的范围由所附权利要求及其等同物限定。

Claims (8)

1.一种复杂地质条件下初始地应力场加速优化反演方法，其特征在于，包括步骤：
S1、建立用于工程区应力场反演计算的三维地质模型；
S2、结合工程区的区域构造特征，确定应力边界条件参数的数值范围，并在该数值范围内进行参数组合得到多组应力边界条件，施加应力边界条件，通过数值计算求得各应力边界条件下的测点地应力计算值；
S3、将应力边界条件作为网络输入，以测点应力计算值作为网络输出，建立代理模型的学习样本，并将该学习样本带入代理模型中进行训练；
S4、基于步骤S3训练好的代理模型，使用遗传算法在样本范围内进行快速寻优，得到最优应力边界条件；
S5、以最优应力边界条件反演计算得到工程区内初始地应力场；
步骤S2中，通过以下方式确定应力边界条件参数的数值范围：
S21、将水平构造应力设置在以下范围：
Figure FDA0004228153710000011
Figure FDA0004228153710000012
其中σT、σt分别为水平最大构造应力及水平最小构造应力，σH、σh分别为实测水平最大主应力及水平最小主应力，γ为岩石容重，H为埋深，μ为泊松比；
S22、设定实测水平最大主应力、水平最小主应力随埋深为以下线性分布：
σH＝A*H+B，σh＝C*H+D；
其中，A、B、C、D为线性回归常数，将上述水平构造应力回归式代入水平构造应力不等式中，得到水平构造应力取值范围为：
Figure FDA0004228153710000021
Figure FDA0004228153710000022
S23、设定水平构造应力与埋深的对应的线性函数为：
σT＝a*H+b，σt＝c*H+d；
结合水平构造应力的取值范围，确定a、b、c、d的取值范围为：
Figure FDA0004228153710000023
Figure FDA0004228153710000024
以a、b、c、d作为应力边界条件参数，获得应力边界条件参数的数值范围；
步骤S4包括：
S41、确定遗传算法的参数，并在应力边界条件范围内随机产生多组应力边界条件作为父代群体；
S42、利用训练好的代理模型对每组应力边界条件进行预测，得到每组应力边界条件所对应的应力值；
S43、基于步骤S42中得到应力值对每组应力边界条件进行适应度评价，适应度均不合理，则迭代产生新的N组应力边界条件并进行适应度评价，直到产生适应度合理的应力边界条件；
S44、在适应度合理的应力边界条件中选择适应度最高的应力边界条件，并通过数值计算得到该应力边界条件下的应力值，将所计算出的应力值与实际测量应力值进行对比，判断两者之间的误差是否满足设定的要求，满足则该应力边界条件即为最优应力边界条件。
2.根据权利要求1所述的复杂地质条件下初始地应力场加速优化反演方法，其特征在于，
步骤S3中，采用径向基神经网络作为代理模型，将应力边界条件和测点高程作为网络输入，测点应力计算值作为网络输出，每个测点的应力计算值包括最大主应力值、最小地应力值及垂直主应力值。
3.根据权利要求1所述的复杂地质条件下初始地应力场加速优化反演方法，其特征在于，
当步骤S44中计算出的应力边界条件下的应力值与实际测量应力值之间的误差不满足设定的要求时，将该应力边界条件及对应的应力值作为一组样本加入到代理模型的训练样本中，重新训练代理模型，同时以该应力边界条件迭代产生新的N组应力边界条件，重复步骤S42到S44。
4.根据权利要求1所述的复杂地质条件下初始地应力场加速优化反演方法，其特征在于，
步骤S43中，当适应度均不合理时，从所评价的应力边界条件中选择适应度最高的应力边界条件进行迭代。
5.根据权利要求1所述的复杂地质条件下初始地应力场加速优化反演方法，其特征在于，
应力边界条件的适应度的评价方法为：使用代理模型预测一组应力边界条件下不同高程对应的水平最大主应力、水平最小主应力、垂直应力，并将这些应力值与实际测量值作差，以差值绝对值之和的负数作为该组应力边界条件的适应度。
6.根据权利要求1所述的复杂地质条件下初始地应力场加速优化反演方法，其特征在于，
使用遗传算法进行寻优，适应性评估函数为负的拟合误差，基因个数为4，即4个应力边界条件，取值范围为(0，1)，种群大小为100，变异率为0.1，选择率为0.6，代数为50，选择轮盘赌策略。
7.根据权利要求1所述的复杂地质条件下初始地应力场加速优化反演方法，其特征在于，
所述步骤S2还包括：通过数值计算求得各应力边界条件下的测点应力计算值后，寻找出相对最优应力边界条件，并将该相对最优应力边界条件对应的应力值与实测应力进行对比，判断两者之间的误差是否满足设定的要求，满足则执行步骤S5，不满足则执行步骤S3。
8.根据权利要求2所述的复杂地质条件下初始地应力场加速优化反演方法，其特征在于，
代理模型使用高斯函数作为径向基函数，高斯函数扩展参数使用固定值为0.1，使用训练数据训练径向基函数，使用K均值聚类进行无监督学习，确定15个聚类中心作为径向基函数的中心，然后计算径向基函数的激活值，随后使用最小二乘法求解径向基函数网络的权重。
CN202210216125.2A 2022-03-07 2022-03-07 复杂地质条件下初始地应力场加速优化反演方法 Active CN114861519B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202210216125.2A CN114861519B (zh) 2022-03-07 2022-03-07 复杂地质条件下初始地应力场加速优化反演方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202210216125.2A CN114861519B (zh) 2022-03-07 2022-03-07 复杂地质条件下初始地应力场加速优化反演方法

Publications (2)

Publication Number Publication Date
CN114861519A CN114861519A (zh) 2022-08-05
CN114861519B true CN114861519B (zh) 2023-06-30

Family

ID=82627416

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202210216125.2A Active CN114861519B (zh) 2022-03-07 2022-03-07 复杂地质条件下初始地应力场加速优化反演方法

Country Status (1)

Country Link
CN (1) CN114861519B (zh)

Families Citing this family (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN115859484B (zh) * 2023-02-23 2023-07-14 西安航天动力研究所 一种发动机力学环境适应性确定方法、装置及电子设备
CN116609828B (zh) * 2023-03-07 2024-03-12 中南大学 深部岩体的应力场计算方法

Citations (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN103605900A (zh) * 2013-11-28 2014-02-26 金川集团股份有限公司 跨尺度复杂地质体地应力场识别方法及装置
CN104965969A (zh) * 2015-04-20 2015-10-07 广西大学 一种大型洞室群围岩力学参数反演方法

Family Cites Families (8)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US7882745B2 (en) * 2006-09-20 2011-02-08 Schlumberger Technology Corporation Method and system to invert tectonic boundary or rock mass field in in-situ stress computation
FR3036518B1 (fr) * 2015-05-20 2018-07-06 Services Petroliers Schlumberger Inversion pour contrainte tectonique
CN106908322B (zh) * 2017-02-23 2019-04-23 成都理工大学 一种基于全应力-应变曲线的岩石脆性指数评价方法
CN106709219B (zh) * 2017-03-06 2019-07-05 中国科学院武汉岩土力学研究所 复杂地质条件下区域初始地应力场反演方法及装置
CN108693572B (zh) * 2018-03-12 2020-07-03 太原理工大学 一种基于三维建模的地应力场反演方法
CN110244354A (zh) * 2019-07-11 2019-09-17 东北大学 一种金属矿山地下开采扰动应力场定量动态反演方法
CN112949000B (zh) * 2021-02-26 2022-08-02 北京理工大学 基于卷积神经网络模型的构件残余应力反演方法
CN113868923A (zh) * 2021-10-13 2021-12-31 西南石油大学 油气储层压前三维地质评价方法

Patent Citations (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN103605900A (zh) * 2013-11-28 2014-02-26 金川集团股份有限公司 跨尺度复杂地质体地应力场识别方法及装置
CN104965969A (zh) * 2015-04-20 2015-10-07 广西大学 一种大型洞室群围岩力学参数反演方法

Also Published As

Publication number Publication date
CN114861519A (zh) 2022-08-05

Similar Documents

Publication Publication Date Title
CN114861519B (zh) 复杂地质条件下初始地应力场加速优化反演方法
CN111563706A (zh) 一种基于lstm网络的多变量物流货运量预测方法
US20040167721A1 (en) Optimal fitting parameter determining method and device, and optimal fitting parameter determining program
CN108563837B (zh) 一种冲积河流水沙模型的模型参数实时校正方法和系统
CN109885916B (zh) 一种基于lssvm的混合试验在线模型更新方法
CN105183928A (zh) 铸铝构件中的残余应力和变形的快速分析
CN112926265A (zh) 基于遗传算法优化神经网络的大气多孔探针测量校准方法
CN112016253A (zh) 一种适用于cfd不确定度量化的高保真度混沌多项式修正方法
CN112329349B (zh) 一种边坡可靠度评估方法
CN112801281A (zh) 基于量子化生成模型和神经网络的对抗生成网络构建方法
CN113032902A (zh) 一种基于机器学习优化的高速列车气动头部外形设计方法
CN112214817B (zh) 考虑层间条件以及横观各向同性的多层位移响应确定方法
CN111914943B (zh) 倾倒式岩溶危岩稳定综合判别的信息向量机方法及装置
CN106568647A (zh) 一种基于神经网络的混凝土强度预测方法
CN115640888A (zh) 一种递减函数嵌入式门限序列网络的产量预测方法
CN113486591B (zh) 一种卷积神经网络结果的重力多参量数据密度加权反演方法
CN112862063A (zh) 一种基于深度信念网络的复杂管网泄漏定位方法
CN116303786B (zh) 一种基于多维数据融合算法的区块链金融大数据管理系统
CN103983332A (zh) 一种基于hgsa-bp算法的传感器误差补偿方法
CN117076887A (zh) 一种泵站机组运行状态预测和健康评估方法及系统
CN113821863B (zh) 一种桩基竖向极限承载力预测方法
CN112016956B (zh) 基于bp神经网络的矿石品位估值方法及装置
CN115100233A (zh) 基于生成对抗网络重采样粒子滤波的雷达目标跟踪方法
CN113673015B (zh) 梁柱端板连接节点优化设计智能系统构建及参数识别方法
Ahmadi Prediction of Orthogonal Cutting Forces Based on Multi-fidelity Modeling

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