CN117607951A - 一种适用于三维隧道模型的多尺度全波形反演方法及系统 - Google Patents
一种适用于三维隧道模型的多尺度全波形反演方法及系统 Download PDFInfo
- Publication number
- CN117607951A CN117607951A CN202311330357.1A CN202311330357A CN117607951A CN 117607951 A CN117607951 A CN 117607951A CN 202311330357 A CN202311330357 A CN 202311330357A CN 117607951 A CN117607951 A CN 117607951A
- Authority
- CN
- China
- Prior art keywords
- wave velocity
- inversion
- dimensional
- objective function
- velocity distribution
- 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
- 238000000034 method Methods 0.000 title claims abstract description 98
- 230000009466 transformation Effects 0.000 claims abstract description 68
- 238000001914 filtration Methods 0.000 claims abstract description 29
- 230000002146 bilateral effect Effects 0.000 claims abstract description 23
- 230000008569 process Effects 0.000 claims description 38
- 238000012937 correction Methods 0.000 claims description 23
- 238000001514 detection method Methods 0.000 claims description 12
- 238000004364 calculation method Methods 0.000 claims description 11
- 239000011435 rock Substances 0.000 claims description 8
- 238000010276 construction Methods 0.000 claims description 6
- 238000012545 processing Methods 0.000 claims description 6
- 230000000694 effects Effects 0.000 claims description 5
- 238000004458 analytical method Methods 0.000 claims description 3
- 238000003825 pressing Methods 0.000 claims description 2
- 230000002159 abnormal effect Effects 0.000 abstract description 2
- 230000006870 function Effects 0.000 description 64
- 239000011159 matrix material Substances 0.000 description 4
- 238000004422 calculation algorithm Methods 0.000 description 3
- 238000010586 diagram Methods 0.000 description 2
- 239000011229 interlayer Substances 0.000 description 2
- 230000035772 mutation Effects 0.000 description 2
- 238000011426 transformation method Methods 0.000 description 2
- 230000009286 beneficial effect Effects 0.000 description 1
- 230000007547 defect Effects 0.000 description 1
- 230000002950 deficient Effects 0.000 description 1
- 238000009795 derivation Methods 0.000 description 1
- 238000013461 design Methods 0.000 description 1
- 238000002474 experimental method Methods 0.000 description 1
- 238000009499 grossing Methods 0.000 description 1
- 230000010354 integration Effects 0.000 description 1
- 238000012986 modification Methods 0.000 description 1
- 230000004048 modification Effects 0.000 description 1
Classifications
-
- 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/282—Application of seismic models, synthetic seismograms
-
- 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
- G01V1/303—Analysis for determining velocity profiles or travel times
-
- 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
- G01V1/306—Analysis for determining physical properties of the subsurface, e.g. impedance, porosity or attenuation profiles
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
- G01V2210/00—Details of seismic processing or analysis
- G01V2210/60—Analysis
- G01V2210/61—Analysis by combining or comparing a seismic data set with other data
- G01V2210/616—Data from specific type of measurement
- G01V2210/6161—Seismic or acoustic, e.g. land or sea measurements
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
- G01V2210/00—Details of seismic processing or analysis
- G01V2210/60—Analysis
- G01V2210/62—Physical property of subsurface
- G01V2210/622—Velocity, density or impedance
- G01V2210/6222—Velocity; travel time
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
- G01V2210/00—Details of seismic processing or analysis
- G01V2210/60—Analysis
- G01V2210/66—Subsurface modeling
Landscapes
- Engineering & Computer Science (AREA)
- Remote Sensing (AREA)
- Physics & Mathematics (AREA)
- Life Sciences & Earth Sciences (AREA)
- Acoustics & Sound (AREA)
- Environmental & Geological Engineering (AREA)
- Geology (AREA)
- General Life Sciences & Earth Sciences (AREA)
- General Physics & Mathematics (AREA)
- Geophysics (AREA)
- Radar Systems Or Details Thereof (AREA)
Abstract
本发明公开了一种适用于三维隧道模型的多尺度全波形反演方法及系统,本发明利用积分变换计算一阶和二阶积分变换观测数据与源子波;利用各阶积分变换观测数据与各阶积分变换源子波进行全波形反演,构造目标函数,计算更新梯度和步长,更新波速模型,计算目标函数,所述目标函数包括数据拟合项与模型拟合项,数据拟合项通过积分变换目标函数计算,模型拟合项包括模糊聚类、全变分和双边滤波的正则化项。本发明能够解决隧道全波形反演的多解性问题,提高隧道预测区域异常地质体的预测准确度和稳定性。
Description
技术领域
本发明涉及隧道地震波超前探测数据处理反演技术领域,特别是涉及一种适用于三维隧道模型的地震波多尺度全波形反演方法及系统。
背景技术
本部分的陈述仅仅是提供了与本发明相关的背景技术信息,不必然构成在先技术。
对于隧道三维地震波单炮主动源探测而言,常规全波形反演的结果是有缺陷的,仅能更新界面位置波速,而难以更新层间波速。常规全波形反演结果对应的正演波形与观测波形十分相近,说明反演结果的缺陷是隧道单炮探测特殊观测模式数据量偏少、偏移距偏小的特殊条件造成的,属于全波形反演的多解性问题,很难从数据角度进一步改善,需要研究专用的反演方法。
发明内容
为了解决上述问题,本发明提出了一种适用于三维隧道模型的多尺度全波形反演方法及系统,本发明能够解决隧道全波形反演的多解性问题,提高隧道预测区域异常地质体的预测准确度和稳定性。
为了实现上述目的,本发明采用如下技术方案:
第一方面,本发明提供一种适用于三维隧道模型的多尺度全波形反演方法,包括以下步骤:
(1)根据隧道周边围岩波速构造均匀初始波速模型,获取隧道单炮主动源探测数据,提取源子波函数;
(2)利用积分变换计算一阶和二阶积分变换观测数据与积分变换源子波;
(3)利用二阶积分变换观测数据与二阶积分变换源子波进行全波形反演,构造目标函数,计算更新梯度和步长,更新波速模型,计算目标函数,所述目标函数包括数据拟合项与模型拟合项,数据拟合项通过积分变换目标函数计算,模型拟合项包括模糊聚类、全变分和双边滤波的正则化项;
(4)迭代进行步骤(3),每迭代到第一设定次数,对反演所得波速分布进行结构校正后,继续迭代,直到达到第二设定次数;
(5)将观测数据和震源子波替换为一阶积分变换观测数据与一阶积分变换震源子波,重复步骤(3)和步骤(4)的反演过程,直到达到第三设定次数;
(6)将观测数据和震源子波替换为原始观测数据和震源子波,重复步骤(3)和步骤(4)的反演过程,直到达到第四设定次数;
(7)将目标函数中的数据拟合项替换为常规目标函数,重复步骤(3)和步骤(4)的反演过程,直到达到第五设定次数,得到的波速分布即为最终反演结果。
作为可选择的实施方式,所述积分变换为三维归一化积分,用于提取观测数据中的低频信息,以实现从低频到高频的多尺度反演,积分变换方法相较于传统低频滤波方法更为稳定,反演效果更优。
作为可选择的实施方式,目标函数的数据拟合项部分通过积分变换目标函数计算,代替了常规的最小二乘目标函数,使用新的标准衡量正演数据和观测数据间的差距,提升了反演过程的稳定性,反演不易陷入局部极小值。
作为可选择的实施方式,所述模糊聚类正则化项的构造过程包括:采用模糊C均值聚类,对反演区域的样本进行聚类,得到聚类中心和每个网格相对于聚类中心的隶属度,基于聚类结果对波速分布进行校正,得到新的波速分布,将聚类分析和正则化方法结合,形成模糊聚类正则化方法,优化反演迭代过程,将目标函数对波速分布求偏导,得到模糊聚类正则化方法的梯度。
作为可选择的实施方式,所述全变分正则化项的构造过程包括:在每次反演迭代过程中,利用全变分去噪处理波速分布,得到新的参考波速分布,并以此构造目标函数中的正则化项。
作为可选择的实施方式,所述双边滤波的正则化项的处理过程包括,将双边滤波的正则化项的目标函数两侧对波速分布求偏导,得到相应的双边滤波正则化梯度。
作为可选择的实施方式,所述目标函数包括数据拟合项与模型拟合项,所述数据拟合项通过积分变换目标函数计算,所述模型拟合项乘以正则化参数,以平衡数据拟合项与模型拟合项的影响。
作为可选择的实施方式,对反演所得波速分布进行结构校正的具体过程包括:将三维波速分布进行z方向切片,得到二维波速分布;
对二维波速分布情况,选取一个纵轴坐标,沿隧道轴线提取一维波速剖面;
将提取出的一维波速剖面自变量由距离转换为时间;
计算一维波速分布相对于时间的导数;
将低于预设下限的波速导数进行压制,变为原有导数数值的设定倍,得到校正后的导数值;
根据校正后的波速导数值,计算以时间为自变量的校正后一维波速分布;
计算以距离为自变量的校正后一维波速分布;
取下一个纵轴网格,重复上述步骤校正水平一维波速分布,直至完成所有网格一维波速剖面的校正。
作为进一步的,进行波速结构校正时,每个网格的一维波速剖面均需依次校正
作为进一步的,将自变量转换为时间后再进行波速校正,使波速突变位置在时间上保持恒定,以保证正演数据中反射波同相轴到时基本不变。
第二方面,本发明提供一种适用于三维隧道模型的多尺度全波形反演系统,包括:
初始模型构建模块,被配置为根据隧道周边围岩波速构造均匀初始波速模型,获取隧道单炮主动源探测数据并提取源子波函数;
积分变换计算模块,被配置为利用积分变换计算一阶和二阶积分变换观测数据与积分变换源子波;
反演计算模块,被配置为利用二阶积分变换观测数据与二阶积分变换源子波进行全波形反演,构造目标函数,计算更新梯度和步长,更新波速模型,计算目标函数,所述目标函数包括数据拟合项与模型拟合项,数据拟合项通过积分变换目标函数计算,模型拟合项包括模糊聚类、全变分和双边滤波的正则化项;
迭代模块,被配置为迭代进行反演计算模块,每迭代到第一设定次数,对反演所得波速分布进行结构校正后,继续迭代,直到达到第二设定次数;
一阶反演模块,被配置为将观测数据和震源子波替换为一阶积分变换观测数据与一阶积分变换震源子波,重复反演过程,直到达到第三设定次数;
观测反演模块,被配置为将观测数据和震源子波替换为原始观测数据和震源子波,重复反演过程,直到达到第四设定次数;
常规反演模块,被配置为将目标函数中的数据拟合项替换为常规目标函数,重复反演过程,直到达到第五设定次数,得到的波速分布即为最终反演结果。
与现有技术相比,本发明的有益效果为:
提出基于二阶和一阶积分变换的多尺度全波形反演方法,利用二阶和一阶积分变换提取三维地震波观测数据中的低频信息,并依照从低频到高频的顺序分步反演,有效提升了全波形反演的稳定性,降低了反演结果的多解性,成功实现了层间波速信息的有效更新,提升了隧道三维地震波探测的反演效果。
提出了基于积分变换目标函数的数据拟合项计算方法,代替了常规的最小二乘目标函数,使用新的标准衡量正演数据和观测数据间的差距,提升了反演过程的稳定性,反演不易陷入局部极小值。
提出波速结构校正方法和以模糊聚类为主的正则化方法,其中反演波速校正针对反演一定次数后的波速分布,直接进行结构校正,使其在正演波形基本不变的情况下波速分布更为接近实际情况。而正则化作用于反演过程之中,可使反演结果更加准确。二者结合,使得最终反演结果更为接近实际波速分布情况。
本发明附加方面的优点将在下面的描述中部分给出,部分将从下面的描述中变得明显,或通过本发明的实践了解到。
附图说明
构成本发明的一部分的说明书附图用来提供对本发明的进一步理解,本发明的示意性实施例及其说明用于解释本发明,并不构成对本发明的不当限定。
图1为本发明实施例1提供的更新波速的流程示意图;
图2为本发明实施例1提供的目标函数随时间偏移量的变换曲线;
图3为本发明实施例1提供的三维波速结构校正方法流程图;
图4为本发明实施例1提供的隧道轴线反演波速剖面;
图5为本发明实施例2提供的三维速度模型;
图6为本发明实施例2提供的三维隧道速度模型与反演速度模型对比图;
图7为本发明实施例2提供的三维隧道速度模型与反演速度模型切片对比图;
图8为本发明实施例2提供的校正后波速正演波形与观测波形对比图;
图9为本发明实施例2提供的波速收敛变化曲线。
具体实施方式
下面结合附图与实施例对本发明做进一步说明。
应该指出,以下详细说明都是示例性的,旨在对本发明提供进一步的说明。除非另有指明,本文使用的所有技术和科学术语具有与本发明所属技术领域的普通技术人员通常理解的相同含义。
需要注意的是,这里所使用的术语仅是为了描述具体实施方式,而非意图限制根据本发明的示例性实施方式。如在这里所使用的,除非上下文另外明确指出,否则单数形式也意图包括复数形式,此外,还应当理解的是,术语“包括”和“具有”以及他们的任何变形,意图在于覆盖不排他的包含,例如,包含了一系列步骤或单元的过程、方法、系统、产品或设备不必限于清楚地列出的那些步骤或单元,而是可包括没有清楚地列出的或对于这些过程、方法、产品或设备固有的其它步骤或单元。
在不冲突的情况下,本发明中的实施例及实施例中的特征可以相互组合。
本发明首先在波速结构校正和正则化反演方法中引入三维网格的概念,然后对波速结构校正以及波速校正、模糊聚类、双边滤波和全变分正则化反演过程中的公式进行推导,改进使其适用于三维模型。
本发明为了解决隧道地震探测模式下全三维波形反演多解性强的问题,对三维波速分布做出了两条假设:假设掌子面前方岩层、异常体等结构内部波速较为均匀;假设结构之间波速变化较为剧烈。
下面以具体实施例进行描述。
实施例1
一种适用于三维隧道模型的多尺度全波形反演方法,本文的多尺度指从低频数据到高频数据分步反演,包括以下步骤:
1、假设已知隧道周边围岩波速、隧道单炮主动源探测数据、源子波函数。根据隧道周边围岩波速构造均匀初始波速模型。
2、利用积分变换计算2阶和1阶的积分变换观测数据与积分变换源子波。
3、利用2阶积分变换观测数据与2阶积分变换源子波进行全波形反演,构造目标函数,计算更新梯度和步长,更新波速分布,如图1所示。目标函数包含数据拟合项与模型拟合项(式(2)),数据拟合项通过归一化积分目标函数计算(式(24)),模型拟合项为模糊聚类(式(7))、全变分(式(14))和双边滤波(式(17))3个正则化项。
下面详细介绍目标函数的相关内容。
式中,F为目标函数,m为模型参数分布,Source为震源,Receive为检波器,Time为接收时间。r为模型参数为m时由震源Source发射,检波器Receive接收,到时为Time的时间域正演波形数据,d为相同位置相同时间的观测数据。
全波形反演正则化方法通常是在原本的关于观测数据的目标函数(式(1))基础上,增加与波速分布相关的目标函数项:
F(v)=Fd(v)+CFinv(v) (2)
式中F为总目标函数,Fd为与观测数据相关的目标函数,Finv为与反演波速分布相关的目标函数,β为平衡二者影响的正则化参数。CFinv(v)为相对于常规全波形反演添加的正则化项。
针对隧道单炮主动源探测模式,根据所做两条波速结构假设,提出了3种正则化方法,分别为模糊聚类(Fuzzy C-Means Cluster ing)正则化方法,全变分(Total Variation)正则化方法和双边滤波(Bi lateral Fi lters)正则化方法:
F(v)=Fd(v)+C1Ffcm(v)+C2Ftv(v)+C3Fbf(v) (3)
式中,Ffcm为模糊聚类正则化项,Ftv为全变分正则化项,Fbf为双边滤波正则化项,C1、C2、C3分别为三种正则化项的控制参数。
模糊C均值聚类正则化
模糊C均值聚类是一种应用广泛的模糊聚类方法,可以处理一组样本,选取数个聚类中心,并计算每个样本对于各聚类中心的隶属度。由于模糊C均值聚类适用于处理一维矩阵,所以采用三维速度模型转换成一维速度矩阵的方式将模糊C均值聚类应用到三维速度模型的处理:
x((i*Nx+j)*Ny+k)=v(i,j,k) (4)
式中x为一维速度矩阵,v为三维速度矩阵,Nx为v的第一维度大小,Ny为v的第二维度大小。
要对Nx*Ny*Nz个样本xi进行聚类,模糊C均值聚类会计算聚类中心cj和样本xi对于聚类中心的隶属度rij,使得下述目标函数最小:
式中C为聚类中心个数,m>1表示模糊水平。将其应用于三维模型中
模糊C均值聚类通过迭代算法不断优化聚类中心位置和隶属度:
当聚类中心和隶属度随迭代变换较小时,认为达到局部最优状态,停止迭代,给出聚类结果。
设置聚类中心个数C=7,模糊水平m=2,使用上式及逆行聚类。得到聚类中心和每个网格相对于聚类中心的隶属度后,基于聚类结果对波速分布进行校正,得到新的波速分布vfcm,i,l:
式中为第j个聚类中心的波速值。基于模糊C均值聚类的波速分布具有明确的地质意义。这里的7个聚类中心相当于7种不同的基础介质。
将聚类分析和正则化方法结合,可以形成模糊聚类正则化方法,优化反演迭代过程:
将目标函数对波速分布求偏导,可得模糊聚类正则化方法的梯度:
全变分正则化
全变分正则化与模糊聚类正则化方法类似,在每次反演迭代过程中,利用全变分去噪处理波速分布,得到新的参考波速分布,并以此构造目标函数中的正则化项。全变分去噪是一种可以保留图像边缘信息的去噪算法,可在保持图形边界信息的基础上降低随机噪声。推导得到三维全变分如下:
式中v是波速分布,dx表示横向微分,dy表示纵向微分,dz表示竖向微分
根据与原波速分布相差不大且全变分降低的原则,可建立全变分去噪目标函数:
式中vtv是速度分布,v0是初始速度,λ是一个调节平滑度程度与和噪声图像差异的值,进一步推导可得:
式中ux、ut、uz表示x、y、z三个方向的波场,vtv是速度分布,v0是初始速度,通过两侧对x、y、z分别求解偏导可得到:
进一步求解可得:
式中vx、vy、vz表示x、y、z方向的速度,vxx、vyy、vzz表示速度对x、y、z方向的二阶偏导,利用迭代法求解全变分去噪后的波速分布vtv:
式中Δh为迭代更新步长。本文取和λ=0并进行1000次迭代。
将全变分去噪结果(式13)代入全波形反演目标函数:
将上式的目标函数两侧对波速分布v求偏导,可以得到相应的全变分正则化梯度为:
双边滤波法正则化
全变分正则化与模糊聚类正则化方法类似,可用双边滤波方法处理每次反演得到的波速分布v,得到滤波后的波速分布vbf。双边滤波是空间域高斯滤波和图像域数值滤波相结合的图像滤波方法。借由二维双边滤波的原理推导出三维双边滤波的核函数可表示为:
式中i,j,k为中心点x、y、z三维坐标,l,m,n为周围点的x、y、z坐标。σd和σv分别为空间位置和波速数值上的平滑因子,本实施例取核函数规模为21*21*21,σd=100且σv=500。利用核函数做卷积处理波速分布v,可得到滤波后的波速分布vbf。
基于双边滤波的全波形反演正则化目标函数为:
将上式的目标函数两侧对波速分布v求偏导,可以得到相应的双边滤波正则化梯度为:
将上述三种正则化方法结合,综合三者的优点,可以得到更优的反演结果。将三种正则化项式(7)、式(14)、式(17)代入到式(3)中可以得到最终正则化目标函数:
并对波速分布vi,j,k求偏导,可得包含三种正则化方法的梯度:
其中,式(7)、式(14)和式(17)中三个正则化参数C1、C2和C3控制着三种正则化方法对反演过程的影响程度,模糊聚类正则化参数Ci的取值公式如下:
其中k1为一常数,可控制正则化强度,通过对比实验选定k1=k2=k3=0.1。
传统的N-阶时间积分波场方法中的N阶积分变换波场在频率域的计算公式为:
式中N为积分阶数,fft为快速傅里叶变换,ifft为快速傅里叶反变换,i为虚数单位,ω=2πf为角频率,UntN为N阶积分算子。
本发明采用的三维归一化积分目标函数是将传统的最小二乘目标函数式(1)替换为一种更加光滑稳定的目标函数,使用一个新的标准衡量正演数据和观测数据间的差距,对于时间域声波全波形反演而言,归一化积分目标函数可表示为:
其中Q为归一化积分算子,其具体表达式如下:
式中Unt1指1阶积分变换(式(22)),max指取最大值。P是一种非负化算子,其作用是将正负交替变化的波形转化为不存在负数的波形。
为了准确衡量正演数据和观测数据的距离大小,交替使用两种softplus变换
P(u)=log(e-r+1) (26)
和
P(u)=log(er+1) (27)
作为非负化算子。
4、每反演迭代一定次数(设为20次),对反演所得波速分布进行结构校正(如图3所示),继续迭代。
具体过程包括:
(1)选取某一个竖轴坐标z0,将三维波速分布沿隧道径向方向提取出二维波速剖面其中x为横轴坐标,y为纵轴坐标。
(2)对二维波速分布情况,选取一个纵轴坐标y0,沿隧道轴线提取一维波速剖面x为横轴坐标。图4即为y0=10m时的一维波速剖面。进行波速结构校正时,每个网格的一维波速剖面均需依次校正。
(3)将提取出的一维波速剖面自变量由距离转换为时间。这里的时间t相当于在最左端发射一个信号,一维传播到距离x所使用的时间。将自变量转换为时间后再进行波速校正,可以使波速突变位置在时间上保持恒定,以保证正演数据中反射波同相轴到时基本不变。对于任意两横纵坐标x1和x2,与其对应的一维波场到达时间存在如下关系:
根据式(27),可将随距离x变化的单道波速剖面转换为随时间t变化。应当注意这种转化仅在一维传播中成立。隧道单道主动源探测主要关心掌子面前方地质情况,将其一维波速剖面通过该变换方法进行处理,不会带来大的误差,可以满足后续的全波形反演迭代要求。
(4)计算一维波速分布相对于时间的导数/>
(5)将低于预设下限的波速导数vt′进行压制,变为原有导数数值的1/5,得到校正后的导数值
(6)根据校正后的波速导数值计算以时间为自变量的校正后一维波速分布
(7)根据式(27),计算以距离为自变量的校正后一维波速分布
(8)取下一个纵轴网格y1,重复第(2)到第(6)步校正水平一维波速分布,直至完成所有网格一维波速剖面的校正。
5、2阶积分变换波场反演达到指定迭代次数后(设为40次),将观测数据和震源子波替换为1阶积分变换观测数据与1阶积分变换震源子波,重复第(3)和第(4)步进行反演。
6、1阶积分变换波场反演达到指定迭代次数后(设为40次),将观测数据和震源子波替换为原始观测数据和震源子波,重复第(3)和第(4)步进行反演。
7、原始波场反演达到指定迭代次数后(设为40次),将目标函数中的数据拟合项替换为常规目标函数(式(1)),重复第(3)和第(4)步进行反演。
8、常规目标函数达到指定迭代次数后(设为40次)所得波速分布即为最终反演结果。
需要说明的是,上述实施例中各个参数的设定值,可以根据具体情况调整。
实施例2
一种适用于三维隧道模型的多尺度全波形反演系统,包括:
初始模型构建模块,被配置为根据隧道周边围岩波速构造均匀初始波速模型,获取隧道单炮主动源探测数据并提取源子波函数;
积分变换计算模块,被配置为利用积分变换计算一阶和二阶积分变换观测数据与源子波;
反演计算模块,被配置为利用二阶积分变换观测数据与二阶积分变换源子波进行全波形反演,构造目标函数,计算更新梯度和步长,更新波速模型,计算目标函数,所述目标函数包括数据拟合项与模型拟合项,数据拟合项通过积分变换目标函数计算,模型拟合项包括模糊聚类、全变分和双边滤波的正则化项;
迭代模块,被配置为迭代进行反演计算模块,每迭代到第一设定次数,对反演所得波速分布进行结构校正后,继续迭代,直到达到第二设定次数;
一阶反演模块,被配置为将观测数据和震源子波替换为一阶积分变换观测数据与一阶积分变换震源子波,重复反演过程,直到达到第三设定次数;
观测反演模块,被配置为将观测数据和震源子波替换为原始观测数据和震源子波,重复反演过程,直到达到第四设定次数;
常规反演模块,被配置为将目标函数中的数据拟合项替换为常规目标函数,重复反演过程,直到达到第五设定次数,得到的波速分布即为最终反演结果。
如图5-图9,为本系统执行结果示意图。
此处需要说明的是,上述模块对应于实施例1中所述的步骤,上述模块与对应的步骤所实现的示例和应用场景相同,但不限于上述实施例1所公开的内容。需要说明的是,上述模块作为系统的一部分可以在诸如一组计算机可执行指令的计算机系统中执行。
在更多实施例中,还提供:
一种电子设备,包括存储器和处理器以及存储在存储器上并在处理器上运行的计算机指令,所述计算机指令被处理器运行时,完成实施例1中所述的方法。为了简洁,在此不再赘述。
应理解,本实施例中,处理器可以是中央处理单元CPU,处理器还可以是其他通用处理器、数字信号处理器DSP、专用集成电路ASIC,现成可编程门阵列FPGA或者其他可编程逻辑器件、分立门或者晶体管逻辑器件、分立硬件组件等。通用处理器可以是微处理器或者该处理器也可以是任何常规的处理器等。
存储器可以包括只读存储器和随机存取存储器,并向处理器提供指令和数据、存储器的一部分还可以包括非易失性随机存储器。例如,存储器还可以存储设备类型的信息。
一种计算机可读存储介质,用于存储计算机指令,所述计算机指令被处理器执行时,完成实施例1中所述的方法。
实施例1中的方法可以直接体现为硬件处理器执行完成,或者用处理器中的硬件及软件模块组合执行完成。软件模块可以位于随机存储器、闪存、只读存储器、可编程只读存储器或者电可擦写可编程存储器、寄存器等本领域成熟的存储介质中。该存储介质位于存储器,处理器读取存储器中的信息,结合其硬件完成上述方法的步骤。为避免重复,这里不再详细描述。
本领域普通技术人员可以意识到,结合本实施例描述的各示例的单元即算法步骤,能够以电子硬件或者计算机软件和电子硬件的结合来实现。这些功能究竟以硬件还是软件方式来执行,取决于技术方案的特定应用和设计约束条件。专业技术人员可以对每个特定的应用来使用不同方法来实现所描述的功能,但是这种实现不应认为超出本申请的范围。
上述虽然结合附图对本发明的具体实施方式进行了描述,但并非对本发明保护范围的限制,所属领域技术人员应该明白,在本发明的技术方案的基础上,本领域技术人员不需要付出创造性劳动即可做出的各种修改或变形仍在本发明的保护范围以内。
Claims (10)
1.一种适用于三维隧道模型的多尺度全波形反演方法,其特征是,包括以下步骤:
(1)根据隧道周边围岩波速构造均匀初始波速模型,获取隧道单炮主动源探测数据并提取源子波函数;
(2)利用积分变换计算一阶和二阶积分变换观测数据与积分变换源子波;
(3)利用二阶积分变换观测数据与二阶积分变换源子波进行全波形反演,构造目标函数,计算更新梯度和步长,更新波速模型,计算目标函数,所述目标函数包括数据拟合项与模型拟合项,数据拟合项通过积分变换目标函数计算,模型拟合项包括模糊聚类、全变分和双边滤波的正则化项;
(4)迭代进行步骤(3),每迭代到第一设定次数,对反演所得波速分布进行结构校正,将三维波速分布进行z方向切片,得到二维波速分布,提取一维波速剖面,并对一维波速剖面进行时间和距离的校正后,继续迭代,直到达到第二设定次数;
(5)将观测数据和震源子波替换为一阶积分变换观测数据与一阶积分变换震源子波,重复步骤(3)和步骤(4)的反演过程,直到达到第三设定次数;
(6)将观测数据和震源子波替换为原始观测数据和震源子波,重复步骤(3)和步骤(4)的反演过程,直到达到第四设定次数;
(7)将目标函数中的数据拟合项替换为常规目标函数,重复步骤(3)和步骤(4)的反演过程,直到达到第五设定次数,得到的波速分布即为最终反演结果。
2.如权利要求1所述的一种适用于三维隧道模型的多尺度全波形反演方法,其特征是,所述积分变换为三维归一化积分,用于提取观测数据中的低频信息,以实现从低频到高频的多尺度反演;所述目标函数数据拟合项通过积分变换目标函数计算,使用新的标准衡量正演数据和观测数据间的差距,提升反演过程的稳定性;
所述目标函数为:
式中,F为目标函数,m为模型参数分布,Source为震源,Receive为检波器,Time为接收时间;r为模型参数为m时由震源Source发射,检波器Receive接收,到时为Time的时间域正演波形数据,d为相同位置相同时间的观测数据。
3.如权利要求1所述的一种适用于三维隧道模型的多尺度全波形反演方法,其特征是,所述模糊聚类正则化项的构造过程包括:采用模糊C均值聚类,对反演区域的样本进行聚类,得到聚类中心和每个网格相对于聚类中心的隶属度,基于聚类结果对波速分布进行校正,得到新的波速分布,将聚类分析和正则化方法结合,形成模糊聚类正则化方法,优化反演迭代过程,将目标函数对波速分布求偏导,得到模糊聚类正则化方法的梯度。
4.如权利要求1所述的一种适用于三维隧道模型的多尺度全波形反演方法,其特征是,所述全变分正则化项的构造过程包括:在每次反演迭代过程中,利用全变分去噪处理波速分布,得到新的参考波速分布,并以此构造目标函数中的正则化项;
三维全变分为
式中v是波速分布,dx表示横向微分,dy表示纵向微分,dz表示竖向微分。
5.如权利要求1所述的一种适用于三维隧道模型的多尺度全波形反演方法,其特征是,所述双边滤波的正则化项的处理过程包括,将双边滤波的正则化项的目标函数两侧对波速分布求偏导,得到相应的双边滤波正则化梯度;
vbf为滤波后的波速分布,波速分布v,i,j,k为中心点x、y、z三维坐标。
6.如权利要求1所述的一种适用于三维隧道模型的多尺度全波形反演方法,其特征是,所述目标函数包括数据拟合项与模型拟合项,所述模型拟合项乘以正则化参数,以平衡数据拟合项与模型拟合项的影响。
7.如权利要求1所述的一种适用于三维隧道模型的多尺度全波形反演方法,其特征是,对反演所得波速分布进行结构校正的具体过程包括:将三维波速分布进行z方向切片,得到二维波速分布;
对二维波速分布情况,选取一个纵轴坐标,沿隧道轴线提取一维波速剖面;
将提取出的一维波速剖面自变量由距离转换为时间;
计算一维波速分布相对于时间的导数;
将低于预设下限的波速导数进行压制,变为原有导数数值的设定倍,得到校正后的导数值;
根据校正后的波速导数值,计算以时间为自变量的校正后一维波速分布;
计算以距离为自变量的校正后一维波速分布;
取下一个纵轴网格,重复上述步骤校正水平一维波速分布,直至完成所有网格一维波速剖面的校正。
8.如权利要求1所述的一种适用于三维隧道模型的多尺度全波形反演方法,其特征是,进行波速结构校正时,每个网格的一维波速剖面均需依次校正。
9.如权利要求1所述的一种适用于三维隧道模型的多尺度全波形反演方法,其特征是,将自变量转换为时间后再进行波速校正,使波速突变位置在时间上保持恒定,以保证正演数据中反射波同相轴到时基本不变。
10.一种适用于三维隧道模型的多尺度全波形反演系统,其特征是,包括:
初始模型构建模块,被配置为根据隧道周边围岩波速构造均匀初始波速模型,获取隧道单炮主动源探测数据并提取源子波函数;
积分变换计算模块,被配置为利用积分变换计算一阶和二阶积分变换观测数据与源子波;
反演计算模块,被配置为利用二阶积分变换观测数据与二阶积分变换源子波进行全波形反演,构造目标函数,计算更新梯度和步长,更新波速模型,计算目标函数,所述目标函数包括数据拟合项与模型拟合项,数据拟合项通过积分变换目标函数计算,模型拟合项包括模糊聚类、全变分和双边滤波的正则化项;
迭代模块,被配置为迭代进行反演计算模块,每迭代到第一设定次数,对反演所得波速分布进行结构校正,将三维波速分布进行z方向切片,得到二维波速分布,提取一维波速剖面,并对一维波速剖面进行时间和距离的校正后,继续迭代,直到达到第二设定次数;
一阶反演模块,被配置为将观测数据和震源子波替换为一阶积分变换观测数据与一阶积分变换震源子波,重复反演过程,直到达到第三设定次数;
观测反演模块,被配置为将观测数据和震源子波替换为原始观测数据和震源子波,重复反演过程,直到达到第四设定次数;
常规反演模块,被配置为将目标函数中的数据拟合项替换为常规目标函数,重复反演过程,直到达到第五设定次数,得到的波速分布即为最终反演结果。
Priority Applications (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN202311330357.1A CN117607951A (zh) | 2023-10-13 | 2023-10-13 | 一种适用于三维隧道模型的多尺度全波形反演方法及系统 |
Applications Claiming Priority (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN202311330357.1A CN117607951A (zh) | 2023-10-13 | 2023-10-13 | 一种适用于三维隧道模型的多尺度全波形反演方法及系统 |
Publications (1)
Publication Number | Publication Date |
---|---|
CN117607951A true CN117607951A (zh) | 2024-02-27 |
Family
ID=89956786
Family Applications (1)
Application Number | Title | Priority Date | Filing Date |
---|---|---|---|
CN202311330357.1A Pending CN117607951A (zh) | 2023-10-13 | 2023-10-13 | 一种适用于三维隧道模型的多尺度全波形反演方法及系统 |
Country Status (1)
Country | Link |
---|---|
CN (1) | CN117607951A (zh) |
-
2023
- 2023-10-13 CN CN202311330357.1A patent/CN117607951A/zh active Pending
Similar Documents
Publication | Publication Date | Title |
---|---|---|
CN108921170B (zh) | 一种有效的图像噪声检测和去噪方法及系统 | |
CN116774292B (zh) | 一种地震波走时确定方法、系统、电子设备及存储介质 | |
CN114663373A (zh) | 一种用于零件表面质量检测的点云配准方法及装置 | |
CN111079893B (zh) | 用于干涉条纹图滤波的生成器网络的获取方法和装置 | |
CN110807428B (zh) | 煤类样品的识别方法、装置、服务器及存储介质 | |
US11662492B2 (en) | Seismic random noise attenuation | |
CN111768349A (zh) | 一种基于深度学习的espi图像降噪方法及系统 | |
CN113917540B (zh) | 基于稀疏约束的抗假频射线束进行地震数据去噪的方法 | |
CN109490954B (zh) | 波场正演模拟方法及装置 | |
CN117607951A (zh) | 一种适用于三维隧道模型的多尺度全波形反演方法及系统 | |
CN116385892A (zh) | 基于目标上下文卷积神经网络的数字高程模型提取方法 | |
CN109856673B (zh) | 一种基于优势频率迭代加权的高分辨Radon变换数据分离技术 | |
CN116068644A (zh) | 一种利用生成对抗网络提升地震数据分辨率和降噪的方法 | |
CN115857018A (zh) | 基于密度聚类的偶极子声波测井频散校正方法及相关装置 | |
CN112578458B (zh) | 叠前弹性阻抗随机反演方法、装置、存储介质及处理器 | |
CN115629426A (zh) | 一种提取局部重力异常的方法、系统、装置及介质 | |
CN113866827A (zh) | 一种解释性速度建模地震成像方法、系统、介质和设备 | |
CN114428343A (zh) | 基于归一化互相关的Marchenko成像方法及系统 | |
CN113219534B (zh) | 一种叠前深度偏移速度质控方法、装置、介质及电子设备 | |
CN117388921B (zh) | 弹性参数的叠前反演方法、装置和电子设备 | |
CN117132478B (zh) | 一种基于法向量二范数特征参数的轨道点云去噪方法 | |
CN111812708B (zh) | 地震波成像方法及装置 | |
CN113050187B (zh) | 滤波方法及装置、计算机设备及计算机可读存储介质 | |
CN116840890A (zh) | 基于随机梯度下降回归算法的全波形反演方法 | |
CN107783184B (zh) | 一种基于多流优化的gpu逆时偏移方法及系统 |
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 |