CN116165707A - 一种基于粘弹性波场模拟的隧道工程探测方法和装置 - Google Patents
一种基于粘弹性波场模拟的隧道工程探测方法和装置 Download PDFInfo
- Publication number
- CN116165707A CN116165707A CN202211538943.0A CN202211538943A CN116165707A CN 116165707 A CN116165707 A CN 116165707A CN 202211538943 A CN202211538943 A CN 202211538943A CN 116165707 A CN116165707 A CN 116165707A
- Authority
- CN
- China
- Prior art keywords
- data
- viscoelastic
- wave
- wave field
- seismic
- 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.)
- Pending
Links
- 238000004088 simulation Methods 0.000 title claims abstract description 30
- 238000001514 detection method Methods 0.000 title claims abstract description 22
- 238000000034 method Methods 0.000 claims abstract description 87
- 238000003384 imaging method Methods 0.000 claims abstract description 38
- 238000013508 migration Methods 0.000 claims abstract description 23
- 230000005012 migration Effects 0.000 claims abstract description 23
- 238000010276 construction Methods 0.000 claims abstract description 20
- 238000013500 data storage Methods 0.000 claims abstract description 5
- 238000013480 data collection Methods 0.000 claims abstract description 4
- 238000004364 calculation method Methods 0.000 claims description 37
- 230000005540 biological transmission Effects 0.000 claims description 31
- 239000011159 matrix material Substances 0.000 claims description 28
- 230000006870 function Effects 0.000 claims description 14
- 238000010521 absorption reaction Methods 0.000 claims description 9
- 238000000354 decomposition reaction Methods 0.000 claims description 9
- 230000005641 tunneling Effects 0.000 claims description 9
- 238000007781 pre-processing Methods 0.000 claims description 8
- 238000012545 processing Methods 0.000 claims description 8
- 230000005284 excitation Effects 0.000 claims description 4
- 238000006243 chemical reaction Methods 0.000 claims description 3
- 230000009466 transformation Effects 0.000 claims description 3
- 230000008569 process Effects 0.000 abstract description 12
- 238000013461 design Methods 0.000 description 2
- 238000010586 diagram Methods 0.000 description 2
- 238000012986 modification Methods 0.000 description 2
- 230000004048 modification Effects 0.000 description 2
- 239000011435 rock Substances 0.000 description 2
- 238000005070 sampling Methods 0.000 description 2
- 238000006467 substitution reaction Methods 0.000 description 2
- 230000009172 bursting Effects 0.000 description 1
- 238000005553 drilling Methods 0.000 description 1
- 230000005611 electricity Effects 0.000 description 1
- 238000013277 forecasting method Methods 0.000 description 1
- 238000012544 monitoring process Methods 0.000 description 1
- 230000000149 penetrating effect Effects 0.000 description 1
- 239000002689 soil Substances 0.000 description 1
- XLYOFNOQVPJJNP-UHFFFAOYSA-N water Substances O XLYOFNOQVPJJNP-UHFFFAOYSA-N 0.000 description 1
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
- G01—MEASURING; TESTING
- G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
- G01V1/00—Seismology; Seismic or acoustic prospecting or detecting
- G01V1/28—Processing seismic data, e.g. for interpretation or for event detection
- G01V1/30—Analysis
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F16/00—Information retrieval; Database structures therefor; File system structures therefor
- G06F16/20—Information retrieval; Database structures therefor; File system structures therefor of structured data, e.g. relational data
- G06F16/21—Design, administration or maintenance of databases
- G06F16/215—Improving data quality; Data cleansing, e.g. de-duplication, removing invalid entries or correcting typographical errors
-
- 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/11—Complex mathematical operations for solving equations, e.g. nonlinear equations, general mathematical optimization problems
-
- 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
-
- 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
-
- Y—GENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
- Y02—TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
- Y02T—CLIMATE CHANGE MITIGATION TECHNOLOGIES RELATED TO TRANSPORTATION
- Y02T90/00—Enabling technologies or technologies with a potential or indirect contribution to GHG emissions mitigation
Landscapes
- Engineering & Computer Science (AREA)
- Physics & Mathematics (AREA)
- General Physics & Mathematics (AREA)
- Theoretical Computer Science (AREA)
- Mathematical Physics (AREA)
- Data Mining & Analysis (AREA)
- Mathematical Optimization (AREA)
- Mathematical Analysis (AREA)
- Computational Mathematics (AREA)
- General Engineering & Computer Science (AREA)
- Databases & Information Systems (AREA)
- Pure & Applied Mathematics (AREA)
- Software Systems (AREA)
- Algebra (AREA)
- Evolutionary Computation (AREA)
- Life Sciences & Earth Sciences (AREA)
- Remote Sensing (AREA)
- Computing Systems (AREA)
- Geometry (AREA)
- Operations Research (AREA)
- Computer Vision & Pattern Recognition (AREA)
- Artificial Intelligence (AREA)
- Quality & Reliability (AREA)
- Computer Hardware Design (AREA)
- Medical Informatics (AREA)
- Acoustics & Sound (AREA)
- Environmental & Geological Engineering (AREA)
- Geology (AREA)
- General Life Sciences & Earth Sciences (AREA)
- Geophysics (AREA)
- Geophysics And Detection Of Objects (AREA)
Abstract
本发明隧道工程地质探测技术领域，具体提供一种基于粘弹性波场模拟的隧道工程探测方法和装置，所述方法包括如下步骤：根据隧道内施工环境的具体情况完成震源和传感器设置；每个震源激发后传感器自动进行数据采集并将数据收集存储，将收集的数据进行去噪处理；设计粘弹性波正演方法求解粘弹性地震波方程；使用存储的数据进行全波形反演和逆时偏移成像。能够更好地拟合波传播的真实物理过程从而提高了地质探测中的反演成像精度，通过使用粘弹波作为波在介质中传播的物理模型的正演方法以及基于粘弹波正演方法的全波形反演方法来实现提高地震波反演成像精度的目的。
Description
技术领域
本发明涉及隧道工程地质探测技术领域，具体涉及一种基于粘弹性波场模拟的隧道工程地质探测方法和装置。
背景技术
隧道快速掘进的主要难题是了解前方的地质情况和岩石力学参数，其中隧道轴线的地质界面可能会在施工掘进中发生严重的问题，在隧道施工中，若遇到不良的地质构造很容易造成突水突泥等灾害事故，尤其是当这些灾害交叉发生时，问题会更加严重。严重威胁了施工人员的生命安全同时导致施工仪器的损毁而造成巨大的财产损失。
为尽量减少隧道施工过程中事故的发生，并保障施工建设的安全，隧道超前地质预报就是解决这个难题行之有效的方法。其中隧道超前地质预报方法中的地球物理超前预报方法在隧道内的使用越来越广泛，地球物理超前预报方法就是在隧道内采用地球物理勘探方法对隧道施工掌子面及周围临近区域进行探测，根据围岩与不良地质体的物理特性差异查明不良地质体的性质、位置及规模。目前，隧道超前地质预报的地球物理探测法主要有高密度电、探地雷达法、地震波法、红外探测法和beam(bore-tunnenng elec-trical aheadmonitoring)法等。其中地震波法由于其具有界面成像精度高、探测距离远等优点已经成为隧道超前预报的重要方法之一。
然而传统地震波法是基于经典波动方程的，但很多实际数据证据表明，波在实际介质中的传播会出现波形衰减的现象，而经典的无衰减的波动方程无法对波形衰减的现象做出解释。
发明内容
传统地震波法是基于经典波动方程的，但很多实际数据证据表明，波在实际介质中的传播会出现波形衰减的现象，而经典的无衰减的波动方程无法对波形衰减的现象做出解释，本发明提供一种基于粘弹性波场模拟的隧道工程地质探测方法和装置。
第一方面，本发明技术方案提供一种基于粘弹性波场模拟的隧道工程探测方法，包括如下步骤：
根据隧道内施工环境的具体情况完成震源和传感器设置；
每个震源激发后传感器自动进行数据采集并将数据收集存储，将收集的数据进行去噪处理；
设计粘弹性波正演方法求解粘弹性地震波方程；
使用存储的数据进行全波形反演和逆时偏移成像。
本申请中设计粘弹性波正演方法求解粘弹性地震波方程，解决非局部波动方程计算瓶颈问题。
作为本发明技术方案的一个优选，将收集的数据进行去噪处理的步骤中使用小波阈值去噪方法进行去噪处理，具体包括：
利用小波分解公式进行小波分解；
利用公式进行阈值λ计算；计算公式：
N为分解的层数；
利用软阈值去噪公式进行去噪处理；公式如下：
其中，其中w为小波系数；
设计小波重构公式利用小波重构公式得到去噪后的信号；
小波重构公式：
作为本发明技术方案的一个优选，设计粘弹性波正演方法求解粘弹性地震波方程的步骤包括：
将探测的区域以及震源和传感器的位置在数值模型中完成将探测的区域以及震源和传感器的位置在数值模型中完成粘弹性地震波方程组构建；
将构建出来的粘弹波方程组使用短记忆分裂格式将方程组转化为一个时间局部、空间非局部的等价问题；
使用样条方法对转化出的问题进行数值求解；
使用辅助微分方程完全匹配层方法对边界进行处理得到吸收边界区域方程组。
作为本发明技术方案的一个优选，粘弹性地震波方程组：
以及时间分数阶应力-应变关系：
其中，
u1和u3是形变分量，ε11、ε13、ε33是应变分量，σ11、σ13、σ33是应力分量，f1、f3是源项；
通过将平坦的奇异核t-2γ转化为相对局域化的奇异核y4γ-1，即令
将方程转化为一个时间局部、空间非局部的等价问题。
作为本发明技术方案的一个优选，使用样条方法对转化出的问题进行数值求解的步骤包括：
样条模拟公式为：
令：
得到：
补充边界条件：
求得γν(yj)；
其中，在实际计算中f(x,y)为实际波场值，s(x,y)为模拟波场值，ην,β是常规二维计算的系数，γν(yj)为转换为一维问题后的系数，Bν(x)为x方向样条基函数，Bβ(y)为y方向样条基函数。
作为本发明技术方案的一个优选，使用辅助微分方程完全匹配层方法对边界进行处理得到吸收边界区域方程组的步骤包括：
设定拉伸坐标系的相关计算公式；
根据傅里叶变换的微分定理得到辅助微分方程；
采用积分求解器对辅助微分方程进行求解，最终得到吸收边界区域方程组。
作为本发明技术方案的一个优选，使用存储的数据进行全波形反演和逆时偏移成像步骤中，进行全波形反演的步骤包括：
利用初始模型计算正传波场数据；
记录模拟观测数据和实际观测数据；
将实际观测数据和模拟观测数据作差计算残差数据；
利用残差数据获取逆传波场数据；
将正传波场数据和逆传波场数据进行互相关计算获取模型误差关于波速的梯度矩阵；
将梯度矩阵导入模型更新程序，利用拟牛顿法进行速度模型的更新，更新速度模型后进行误差函数的计算，直至误差函数的计算结果小于或等于准许误差，获得反演速度模型。
作为本发明技术方案的一个优选，将正传波场数据和逆传波场数据进行互相关计算获取模型误差关于波速的梯度矩阵的步骤包括：
将获取的正传波场数据d正和逆传波场数据d逆在每个时间点在二维空间上进行对应相乘运算，再将每个时间节点进行相加运算，即进行互相关运算，导出互相关计算的结果即为全波形反演的梯度矩阵J：
其中，d正是正传波场数据，d逆是逆传波场数据。
作为本发明技术方案的一个优选，使用存储的数据进行全波形反演和逆时偏移成像的步骤中，进行逆时偏移成像的步骤包括：
利用反演所得模型进行逆时偏移成像，在反演速度模型中进行波场正传；
将采集到的数据进行波场的逆传；
将正传波场数据和逆传波场数据进行互相关计算得到逆时偏移成像结果。
本发明技术方案提供一种基于粘弹性波场模拟的隧道工程探测装置，包括数据采集与预处理模块、正演模块、反演成像模块；
震源模块，根据隧道内施工环境的具体情况设置在施工隧道内，用于激发后输出模拟地震信号；
数据采集与预处理模块，用于对每个震源激发后模拟的地震信号自动进行数据采集并将数据收集存储，将收集的数据进行去噪处理；
正演模块，用于设计粘弹性波正演方法求解粘弹性地震波方程；
反演成像模块，用于使用存储的数据进行全波形反演和逆时偏移成像。
基于粘弹性波动方程，相对于经典波动方程在研究诸如波在粘弹性、滞弹性介质中的传播等问题时，能够更好地拟合波传播的真实物理过程从而提高了地质探测中的反演成像精度，通过使用粘弹波作为波在介质中传播的物理模型的正演方法以及基于粘弹波正演方法的全波形反演方法来实现提高地震波反演成像精度的目的。
从以上技术方案可以看出，本发明具有以下优点：1、在波动方程模型方面，本公开利用的粘弹波模型项对于经典波动方程模型而言，能够更加切合波在真实岩土介质中的传播过程，因此在物理模型构建上更为合理；2、在求解分数阶波动方程方面本公开设计的短记忆算子分裂(SMOS)格式，在实际计算中表明，当分数阶很小时(≤0.06)，SMOS可以用32层记忆变量达到原有方法500层记忆变量的精度(L^2误差在2％左右)，并且降低了存储量和计算时间，极大地降低了计算成本；3、在反演成像方面，由于使用的物理模型更加切合实际，所以相对于使用经典波方程作为正演物理模型来说成像的分辨率会更高。
此外，本发明设计原理可靠，结构简单，具有非常广泛的应用前景。
由此可见，本发明与现有技术相比，具有突出的实质性特点和显著地进步，其实施的有益效果也是显而易见的。
附图说明
为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例或现有技术描述中所需要使用的附图作简单地介绍，显而易见地，对于本领域普通技术人员而言，在不付出创造性劳动的前提下，还可以根据这些附图获得其他的附图。
图1是本发明一个实施例的方法的示意性流程图。
图2是本发明一个实施例的方法中进行逆时偏移成像的示意性流程图。
图3是本发明一个实施例的方法中进行全波形反演的示意性流程图。
图4是本发明一个实施例的装置的示意性框图。
图5是本发明实施例的隧道实例示意图。
具体实施方式
传统地震波法是基于经典波动方程的，但很多实际数据证据表明，波在实际介质中的传播会出现波形衰减的现象，而经典的无衰减的波动方程无法对这一现象做出解释。本申请提供基于粘弹性波动方程，相对于经典波动方程在研究诸如波在粘弹性、滞弹性介质中的传播等问题时，能够更好地拟合波传播的真实物理过程从而提高了地质探测中的反演成像精度。为了使本技术领域的人员更好地理解本发明中的技术方案，下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都应当属于本发明保护的范围。
如图1所示，本发明实施例提供一种基于粘弹性波场模拟的隧道工程探测方法，包括如下步骤：
步骤1：根据隧道内施工环境的具体情况完成震源和传感器设置；
步骤2：每个震源激发后传感器自动进行数据采集并将数据收集存储，将收集的数据进行去噪处理；
步骤3：设计粘弹性波正演方法求解粘弹性地震波方程；
步骤4：使用存储的数据进行全波形反演和逆时偏移成像。
本申请中设计粘弹性波正演方法求解粘弹性地震波方程，解决非局部波动方程计算瓶颈问题。
本发明实施例中，将收集的数据进行去噪处理的步骤中使用小波阈值去噪方法进行去噪处理，具体包括：
步骤21：利用小波分解公式进行小波分解；首先用阶梯信号对原始信号进行逼近：
将fj分解为：fj＝ωj-1+fj-1
其中：
最终可得：fj＝ωj-1+ωj-2+…+ω0+f0。
步骤22：利用公式进行阈值λ计算；计算公式：
N为分解的层数；
步骤23：利用软阈值去噪公式进行去噪处理；公式如下：
其中，其中w为小波系数；
步骤24：设计小波重构公式利用小波重构公式得到去噪后的信号；
小波重构公式：
在本发明实例中，设计粘弹性波正演方法求解粘弹性地震波方程的步骤包括：
步骤31：将探测的区域以及震源和传感器的位置在数值模型中完成将探测的区域以及震源和传感器的位置在数值模型中完成粘弹性地震波方程组构建；
步骤32：将构建出来的粘弹波方程组使用短记忆分裂格式将方程组转化为一个时间局部、空间非局部的等价问题；
设计了二维分数阶波动方程的短记忆分裂(Short-Memory Operator Splitting,SMOS)格式；粘弹性地震波方程组：
以及时间分数阶应力-应变关系：
其中，
u1和u3是形变分量，ε11、ε13、ε33是应变分量，σ11、σ13、σ33是应力分量，f1、f3是源项；
通过将平坦的奇异核t-2γ转化为相对局域化的奇异核y4γ-1，即令
将方程转化为一个时间局部、空间非局部的等价问题。
步骤33：使用样条方法对转化出的问题进行数值求解；
使用三次B样条基函数：
该公式可以由Cox-de Boor递归公式获得：
对于一维问题：
其中i＝-1与n+1为边界条件；计算得n＝1,2,…,n-1时
增加如下两个边界条件：
方程组的矩阵计算的形式为：
其中A为：
C和F为：
C＝[c-1,c0,c1,…,cn-1,cn,cn+1]T
F＝[f'(u0),f(u0),f(u1),…,f(un-1),f(un),f'(un+1)]T
对于二维问题，为了加快计算速度，不再使用原来的二维样条模拟方法，改为将二维问题变为一维问题的组合的方式进行求解。
样条模拟公式为：
令：
得到：
补充边界条件：
求得γν(yj)。利用相似的步骤进行第二次求解即可求得所有系数。整个系数计算过程每次只需进行一次求解过程。系数求解完毕之后，将系数与基函数的导数相乘即获得近似空间导数。
步骤34：使用辅助微分方程完全匹配层方法对边界进行处理得到吸收边界区域方程组，具体包括：
步骤341：设定拉伸坐标系的相关计算公式：
其中，
代表拉伸坐标系，此外有：
ki＝1+(kmax-1)m
其中，L为PML层的厚度，N＝2,σ0＝-(N+1)vpmaxlog(R)/2L,vpmax为压力波的速度，R为目标理论反射系数，此处选择为0.1％。取p＝1，αmax＝πf0，其中f0为震源的主频率。kmax通常在1到20之间，我们取m＝2。拉伸坐标下的频率域一阶声波方程：
在频率域有：
令：
得到：
步骤342：根据傅里叶变换的微分定理得到辅助微分方程；辅助微分方程：
整个PML区域的声波方程组为：
步骤343：采用积分求解器对辅助微分方程进行求解，最终得到吸收边界区域方程组：
本发明实施例中，使用存储的数据进行全波形反演和逆时偏移成像步骤中，如图2所示，进行全波形反演的步骤包括：
步骤41：利用初始模型计算正传波场数据，记录模拟观测数据和实际观测数据；并将实际观测数据和模拟观测数据作差计算残差数据；
利用初始模型计算正传波场数据d正，记录模拟观测数据ob模拟，记实际观测数据为ob实际，利用公式计算残差数据：ob残差＝ob实际-ob模拟；
步骤42:利用残差数据ob残差获取逆传波场数据d逆；
步骤43：将正传波场数据和逆传波场数据进行互相关计算获取模型误差关于波速的梯度矩阵；
将获取的正传波场数据d正和逆传波场数据d逆在每个时间点在二维空间上进行对应相乘运算，再将每个时间节点进行相加运算，即进行互相关运算，导出互相关计算的结果即为全波形反演的梯度矩阵J：
其中，d正是正传波场数据，d逆是逆传波场数据。
步骤44：将梯度矩阵导入模型更新程序，利用拟牛顿法进行速度模型的更新，更新速度模型后进行误差函数的计算，直至误差函数的计算结果小于或等于准许误差，获得反演速度模型；具体包括：
将梯度矩阵导入模型更新程序，利用拟牛顿法计算速度模型的更新矩阵dm，首先计算搜索方向b：
然后计算速度更新矩阵dm：λ＝arg minf(m+λb)
dm＝λb
其中D为近似海森矩阵，其更新方式为：dJ＝J(t+1)-J(t)
将原速度模型矩阵和更新矩阵相加获得更新后的速度模型矩阵m1，利用最小二乘方法计算模型误差，将更新后模型参数导入正演程序，记录模拟地震数据，计算更新后模型的地震数据和更新前地震数据的差的二范数，可表达为：Φ＝||f(m1)-f(m0)||2；
其中Φ为计算得的速度模型的误差，将允许模型的准许误差ε设为0.01，若Φ>ε则进行下一步反演计算，将更新后的速度模型再进行一次上述的梯度计算过程以及模型更新过程，直至Φ≤ε，即获得反演速度模型m。
本发明实施例中，使用存储的数据进行全波形反演和逆时偏移成像的步骤中，如图3所示，进行逆时偏移成像的步骤包括：
步骤45：利用反演所得模型进行逆时偏移成像，在反演速度模型中进行波场正传；
步骤46：将采集到的数据进行波场的逆传；
步骤47：将正传波场数据和逆传波场数据进行互相关计算得到逆时偏移成像结果。
如图4所示，本发明实施例提供一种基于粘弹性波场模拟的隧道工程探测装置，包括震源模块、数据采集与预处理模块、正演模块、反演成像模块；
震源模块，根据隧道内施工环境的具体情况设置在施工隧道内，用于激发后输出模拟地震信号；
数据采集与预处理模块，用于对每个震源激发后模拟的地震信号自动进行数据采集并将数据收集存储，将收集的数据进行去噪处理；
正演模块，用于设计粘弹性波正演方法求解粘弹性地震波方程；
反演成像模块，用于使用存储的数据进行全波形反演和逆时偏移成像。隧道实例示意图如图5所示。
第一步通过震源模块，产生主动源地震波信号，震源可使用锤击震源或气动震源等，震源激发频率约在200到400Hz，在本实例中每侧边墙布置两个锤击震源，震源间的间隔为2m。
第二步通过数据采集与预处理模块对地震信号接收、传输、去噪处理。在本实例中于隧道前方边墙距掌子面约5m处进行传感器布置，每隔2m的距离使用电钻进行钻孔，每隔2m，钻一个孔，钻孔直径约为50mm，钻孔的深度约为1.60m，总布置长度为10m，每侧边墙上有5个传感器；在这里传感器也可以是检波器。
每个震源激后数据采集模块自动进行一次数据采集，每次信号接收时长为0.1s，采样频率为8000次，且每次采集到的数据自动储存直至全部震源激发完毕。数据储存完毕后，对收集到的数据进行去噪处理，在预处理模块，使用小波阈值去噪方法进行去噪处理，利用上文中小波分解公式进行小波分解，选用Haar、Daubechies、Coiflets、Symlets四种小波函数将地震数据进行32层小波分解，再进行阈值计算，计算公式为：
其中N为分解的层数，计算得λ＝1.735。
再利用软阈值去噪公式进行去噪处理：
其中，其中w为小波系数；
最后利用上文中的小波重构公式即得去噪后的信号。
第三步是通过正演模型构建模块进行正演模型构建，在本实施例中，将探测的区域以及震源和传感器的位置在数值模型中完成构建，震源设置为雷克子波，公式为：
其中fp为震源主频为300Hz，根据采样频率可得时间间隔dt为：
将构建出来的粘弹波方程组使用短记忆分裂(Short-Memory OperatorSplitting,SMOS)格式进行求解，在数值模拟过程中使用样条方法计算空间导数，使用差分方法进行计算时间导数，吸收边界使用辅助微分方程完美匹配层(ADE-Pml)。
第四步进行反演计算，首先给速度初始模型m0；导入初始速度模型矩阵，在初始速度模型中进行一次正演模拟，模拟时长为0.1s，模拟的时间节点是0s到0.1s，记录模拟波场数据和模拟地震记录，再将模拟地震记录和之前储存的传感器数据进行相减，获得残差数据；利用残差数据在波场中进行逆传波场模拟，模拟时长为0.1s，模拟的时间节点是0.1s到0s，同样记录逆传波场的模拟波场数据；将获取的正传波场数据和逆传波场数据在每个时间点在二维空间上进行对应相乘运算，再将每个时间节点进行相加运算，即进行互相关运算：
导出互相关计算的结果即为全波形反演的梯度矩阵J(目标函数对速度的梯度)，将梯度矩阵导入模型更新程序，利用拟牛顿法计算速度模型的更新矩阵dm，首先计算搜索方向b：
然后计算速度更新矩阵dm：
λ＝arg minf(m+λb)
dm＝λb
其中D为近似海森矩阵，其更新方式为：dJ＝J(t+1)-J(t)
将原速度模型矩阵和更新矩阵相加获得更新后的速度模型矩阵m1，利用最小二乘方法计算模型误差，将更新后模型参数导入正演程序，记录模拟地震数据，计算更新后模型的地震数据和更新前地震数据的差的二范数，可表达为：
Φ＝||f(m1)-f(m0)||2
其中Φ为计算得的速度模型的误差，将允许模型的准许误差设为0.01，若计算所得模型误差大于准许误差，进行下一步反演计算，将更新后的速度模型再进行一次上述的梯度计算过程以及模型更新过程，直至模型误差小于准许误差，得到反演速度模型m，在速度模型m中进行波场正传，再将采集到的数据进行波场的逆传，最终将正传波场数据和逆传波场数据进行互相关计算得到逆时偏移成像结果。
尽管通过参考附图并结合优选实施例的方式对本发明进行了详细描述，但本发明并不限于此。在不脱离本发明的精神和实质的前提下，本领域普通技术人员可以对本发明的实施例进行各种等效的修改或替换，而这些修改或替换都应在本发明的涵盖范围内/任何熟悉本技术领域的技术人员在本发明揭露的技术范围内，可轻易想到变化或替换，都应涵盖在本发明的保护范围之内。因此，本发明的保护范围应以所述权利要求的保护范围为准。
Claims (10)
1.一种基于粘弹性波场模拟的隧道工程探测方法，其特征在于，包括如下步骤：
根据隧道内施工环境的具体情况完成震源和传感器设置；
每个震源激发后传感器自动进行数据采集并将数据收集存储，将收集的数据进行去噪处理；
设计粘弹性波正演方法求解粘弹性地震波方程；
使用存储的数据进行全波形反演和逆时偏移成像。
2.根据权利要求1所述的基于粘弹性波场模拟的隧道工程探测方法，其特征在于，将收集的数据进行去噪处理的步骤中使用小波阈值去噪方法进行去噪处理，具体包括：
利用小波分解公式进行小波分解；
利用公式进行阈值计算；
利用软阈值去噪公式进行去噪处理；
设计小波重构公式利用小波重构公式得到去噪后的信号。
3.根据权利要求2所述的基于粘弹性波场模拟的隧道工程探测方法，其特征在于，设计粘弹性波正演方法求解粘弹性地震波方程的步骤包括：
将探测的区域以及震源和传感器的位置在数值模型中完成将探测的区域以及震源和传感器的位置在数值模型中完成粘弹性地震波方程组构建；
将构建出来的粘弹波方程组使用短记忆分裂格式将方程组转化为一个时间局部、空间非局部的等价问题；
使用样条方法对转化出的问题进行数值求解；
使用辅助微分方程完全匹配层方法对边界进行处理得到吸收边界区域方程组。
4.根据权利要求3所述的基于粘弹性波场模拟的隧道工程探测方法，其特征在于，粘弹性地震波方程组：
以及时间分数阶应力-应变关系：
其中，
u1和u3是形变分量，ε11、ε13、ε33是应变分量，σ11、σ13、σ33是应力分量，f1、f3是源项；
通过将平坦的奇异核t-2γ转化为相对局域化的奇异核y4γ-1，即令
将方程转化为一个时间局部、空间非局部的等价问题。
5.根据权利要求4所述的基于粘弹性波场模拟的隧道工程探测方法，其特征在于，使用样条方法对转化出的问题进行数值求解的步骤包括：
样条模拟公式为：
令：
得到：
补充边界条件：
求得γν(yj)；
其中，在实际计算中f(x,y)为实际波场值，s(x,y)为模拟波场值，ην,β是常规二维计算的系数，γν(yj)为转换为一维问题后的系数，Bν(x)为x方向样条基函数，Bβ(y)为y方向样条基函数。
6.根据权利要求5所述的基于粘弹性波场模拟的隧道工程探测方法，其特征在于，使用辅助微分方程完全匹配层方法对边界进行处理得到吸收边界区域方程组的步骤包括：
设定拉伸坐标系的相关计算公式；
根据傅里叶变换的微分定理得到辅助微分方程；
采用积分求解器对辅助微分方程进行求解，最终得到吸收边界区域方程组。
7.根据权利要求6所述的基于粘弹性波场模拟的隧道工程探测方法，其特征在于，使用存储的数据进行全波形反演和逆时偏移成像步骤中，进行全波形反演的步骤包括：
利用初始模型计算正传波场数据；
记录模拟观测数据和实际观测数据；
将实际观测数据和模拟观测数据作差计算残差数据；
利用残差数据获取逆传波场数据；
将正传波场数据和逆传波场数据进行互相关计算获取模型误差关于波速的梯度矩阵；
将梯度矩阵导入模型更新程序，利用拟牛顿法进行速度模型的更新，更新速度模型后进行误差函数的计算，直至误差函数的计算结果小于或等于准许误差，获得反演速度模型。
8.根据权利要求7所述的基于粘弹性波场模拟的隧道工程探测方法，其特征在于，将正传波场数据和逆传波场数据进行互相关计算获取模型误差关于波速的梯度矩阵的步骤包括：
将获取的正传波场数据d正和逆传波场数据d逆在每个时间点在二维空间上进行对应相乘运算，再将每个时间节点进行相加运算，即进行互相关运算，导出互相关计算的结果即为全波形反演的梯度矩阵J：
其中，d正是正传波场数据，d逆是逆传波场数据。
9.根据权利要求8所述的基于粘弹性波场模拟的隧道工程探测方法，其特征在于，使用存储的数据进行全波形反演和逆时偏移成像的步骤中，进行逆时偏移成像的步骤包括：
利用反演所得模型进行逆时偏移成像，在反演速度模型中进行波场正传；
将采集到的数据进行波场的逆传；
将正传波场数据和逆传波场数据进行互相关计算得到逆时偏移成像结果。
10.一种基于粘弹性波场模拟的隧道工程探测装置，其特征在于，震源模块、数据采集与预处理模块、正演模块、反演成像模块；
震源模块，根据隧道内施工环境的具体情况设置在施工隧道内，用于激发后输出模拟地震信号；
数据采集与预处理模块，用于对每个震源激发后模拟的地震信号自动进行数据采集并将数据收集存储，将收集的数据进行去噪处理；
正演模块，用于设计粘弹性波正演方法求解粘弹性地震波方程；
反演成像模块，用于使用存储的数据进行全波形反演和逆时偏移成像。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202211538943.0A CN116165707A (zh) | 2022-12-02 | 2022-12-02 | 一种基于粘弹性波场模拟的隧道工程探测方法和装置 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202211538943.0A CN116165707A (zh) | 2022-12-02 | 2022-12-02 | 一种基于粘弹性波场模拟的隧道工程探测方法和装置 |
Publications (1)
| Publication Number | Publication Date |
|---|---|
| CN116165707A true CN116165707A (zh) | 2023-05-26 |
Family
ID=86413846
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202211538943.0A Pending CN116165707A (zh) | 2022-12-02 | 2022-12-02 | 一种基于粘弹性波场模拟的隧道工程探测方法和装置 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN116165707A (zh) |
Cited By (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN117406278A (zh) * | 2023-11-01 | 2024-01-16 | 北京美方信科技有限公司 | 一种粘弹性流体因子叠前地震反演方法 |
-
2022
- 2022-12-02 CN CN202211538943.0A patent/CN116165707A/zh active Pending
Cited By (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN117406278A (zh) * | 2023-11-01 | 2024-01-16 | 北京美方信科技有限公司 | 一种粘弹性流体因子叠前地震反演方法 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| EP2572214B1 (en) | Passive monitoring method for seismic events | |
| CN114035228B (zh) | 一种基于深度学习的隧道地震波速反演方法及系统 | |
| Pitarka et al. | Analysis of ground motion from an underground chemical explosion | |
| US11789173B1 (en) | Real-time microseismic magnitude calculation method and device based on deep learning | |
| US9354337B2 (en) | System, machine, and computer-readable storage medium for forming an enhanced seismic trace using a virtual seismic array | |
| NO341717B1 (no) | Stakking av seismiske støydata for å analysere mikroseismiske hendelser | |
| CN116165707A (zh) | 一种基于粘弹性波场模拟的隧道工程探测方法和装置 | |
| CN103913768A (zh) | 基于地震波资料对地表中浅层进行建模的方法及装置 | |
| CN104199088B (zh) | 一种提取入射角道集的方法及系统 | |
| CN117826243A (zh) | 融入混合注意力Unet神经网络隧洞地震数据重建方法 | |
| Liu et al. | An inverted heterogeneous velocity model for microseismic source location in deep buried tunnels | |
| US20240230942A9 (en) | Method and System for Real-Time Calculating a Microseismic Focal Mechanism Based on Deep Learning | |
| Babacan et al. | Evaluation of soil liquefaction potential with a holistic approach: a case study from Araklı (Trabzon, Turkey). | |
| CN115980851A (zh) | 复合震源参数反演的方法、计算机设备及可读存储介质 | |
| CN115453616A (zh) | 基于深度学习的隧道地震勘探数据增强方法及系统 | |
| Chen et al. | Seismic ahead-prospecting based on deep learning of retrieving seismic wavefield | |
| Chen et al. | A frequency-domain full waveform inversion method of elastic waves in quantitative defection investigation | |
| CN113031072A (zh) | 虚同相轴层间的多次波压制方法、装置及设备 | |
| Chakravarty et al. | Hydraulic Fracturing-Driven Infrasound Signals–a New Class of Signal for Subsurface Engineering | |
| Parasyris et al. | Near surface full waveform inversion via deep learning for subsurface imaging | |
| Aditya et al. | Hydraulic Fracturing-driven Infrasound Signals-A New Class of Signal for Subsurface Engineering | |
| Wang et al. | Time-lapse global inversion for surface-waves: A differential approach using a linear approximation of the Rayleigh wave phase velocity | |
| Wang et al. | Seismic Wavelet Signal Noise Reduction Algorithm of Blind Surce Separation Optimization | |
| CN114076990B (zh) | 油页岩反射能量确定方法、系统、存储介质及电子设备 | |
| CN117726182B (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 |