CN111665552A - 用于山体滑坡危险性评价的声学参数获取方法 - Google Patents
用于山体滑坡危险性评价的声学参数获取方法 Download PDFInfo
- Publication number
- CN111665552A CN111665552A CN201910172714.3A CN201910172714A CN111665552A CN 111665552 A CN111665552 A CN 111665552A CN 201910172714 A CN201910172714 A CN 201910172714A CN 111665552 A CN111665552 A CN 111665552A
- Authority
- CN
- China
- Prior art keywords
- gradient
- data
- parameter
- wave
- error
- 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 48
- 238000011156 evaluation Methods 0.000 title claims abstract description 9
- 238000004088 simulation Methods 0.000 claims abstract description 9
- 238000001914 filtration Methods 0.000 claims description 7
- 239000011159 matrix material Substances 0.000 claims description 7
- 238000012545 processing Methods 0.000 claims description 6
- 239000006185 dispersion Substances 0.000 claims description 5
- 230000001629 suppression Effects 0.000 claims description 5
- 238000001514 detection method Methods 0.000 claims description 4
- 230000005484 gravity Effects 0.000 claims description 4
- 238000012502 risk assessment Methods 0.000 claims description 4
- 238000012886 linear function Methods 0.000 claims description 3
- 230000006870 function Effects 0.000 description 7
- 238000005070 sampling Methods 0.000 description 5
- 230000009466 transformation Effects 0.000 description 4
- 238000004364 calculation method Methods 0.000 description 3
- 238000012937 correction Methods 0.000 description 3
- 238000011161 development Methods 0.000 description 3
- 230000008569 process Effects 0.000 description 3
- 238000004422 calculation algorithm Methods 0.000 description 2
- 238000010586 diagram Methods 0.000 description 2
- 230000001902 propagating effect Effects 0.000 description 2
- 239000000523 sample Substances 0.000 description 2
- 230000002159 abnormal effect Effects 0.000 description 1
- 238000013459 approach Methods 0.000 description 1
- 238000003491 array Methods 0.000 description 1
- 230000008859 change Effects 0.000 description 1
- 230000001427 coherent effect Effects 0.000 description 1
- 238000002939 conjugate gradient method Methods 0.000 description 1
- 238000010276 construction Methods 0.000 description 1
- 230000007547 defect Effects 0.000 description 1
- 238000013461 design Methods 0.000 description 1
- 238000006073 displacement reaction Methods 0.000 description 1
- 238000005516 engineering process Methods 0.000 description 1
- 230000005284 excitation Effects 0.000 description 1
- 238000003384 imaging method Methods 0.000 description 1
- 229910052500 inorganic mineral Inorganic materials 0.000 description 1
- 239000011707 mineral Substances 0.000 description 1
- 238000012986 modification Methods 0.000 description 1
- 230000004048 modification Effects 0.000 description 1
- 230000000704 physical effect Effects 0.000 description 1
- 238000007781 pre-processing Methods 0.000 description 1
- 239000011435 rock Substances 0.000 description 1
- 238000006467 substitution reaction Methods 0.000 description 1
- 230000002123 temporal effect Effects 0.000 description 1
- XLYOFNOQVPJJNP-UHFFFAOYSA-N water Substances O XLYOFNOQVPJJNP-UHFFFAOYSA-N 0.000 description 1
Images
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
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)
- Geophysics And Detection Of Objects (AREA)
Abstract
本发明涉及一种用于山体滑坡危险性评价的声学参数获取方法，包括：获取山体探测区域内震源震动的多个单炮数据；利用时窗对多个单炮数据中的直达波、浅层反射波、浅层折射波进行截取，得到观测数据；获取声学参数初始模型，基于声学参数初始模型对震源波形进行正演，得到正演模拟数据；根据观测数据和正演模拟数据计算波场残差，根据波场残差构造误差泛函；利用波场残差反传播到声学参数初始模型空间，得到残差反传播数据；利用正演模拟数据和残差反传播数据分别计算误差泛函的速度梯度、密度梯度、波阻抗梯度、衰减因子梯度、拉梅参数梯度；分别利用各个参数梯度更新声学参数初始模型，得到每个声学参数对应的精确模型。
Description
技术领域
本发明涉及地球物理勘探技术领域，尤其涉及一种用于山体滑坡危险性评价的声学参数获取方法。
背景技术
矿产和围岩具有明显的声学参数差异,不同的物体具有不同的声学参数特征，声学参数值各自具有一定的范围，可以根据不同物体的声学参数差异判断各自的物性。
目前较为常用的是基于射线类的初至走时层析反演的方法。这类方法计算简单，不要求准确的背景场，且应用复杂地表情况，因此一直以来是解决近地表及浅层声学参数建模的常用方法。但这种方法是基于高频假设的最短路径原理，在低速体发育的情况下，该技术存在建模的‘盲区’。此外，射线类方法对射线密度敏感，当速度变化剧烈，即使是高速异常体，由于全发射现象的发生，将严重影响射线的密度，继而降低反演的精度。更何况，在实际资料应用中，初至的拾取工作量巨大，人工拾取存在误差，且当地表复杂的情况下初至很难拾取准确。
另外一种解决这类问题的新方法就是全波形反演的方法，这种方法基于波动方程，可以较为真实的模拟波的传播波场，因此，不受射线密度影响。但这一新方法的实现还有很多的限制，如要求大偏移距的观测系统、低频信息缺失等。全波形反演的实用化还有很长的路程，尤其陆上资料的实用化。
对山体及其附近的状况，需要定时或者不定时进行监测，根据监测到的信息判断是否可能发生山体滑坡。例如，在地震/洪水以后，通过判断裂隙或者含水性/量，结合山体形状，可以监测危险程度，评估山体滑坡的可能性和进行预警。目前缺少一种获取用于山体滑坡危险性评价的声学参数的解决方案。
发明内容
本发明的目的是针对现有技术中存在的缺陷，提供一种用于山体滑坡危险性评价的声学参数获取方法。
为实现上述目的，第一方面,本发明提供了一种用于山体滑坡危险性评价的声学参数获取方法，包括：
获取山体探测区域内震源震动的多个单炮数据，所述单炮数据包括单炮单道数据和单炮多道数据；
利用时窗对所述多个单炮数据中的直达波、浅层反射波、浅层折射波进行截取，得到初始观测数据；
对所述初始观测数据进行滤波处理和多次波压制处理，得到观测数据；
获取声学参数初始模型，并基于所述声学参数初始模型对震源波形进行正演，得到正演模拟数据，所述声学参数包括声波传播速度、介质密度、波阻抗、衰减因子和拉梅参数；
根据所述观测数据和所述正演模拟数据计算波场残差，并根据所述波场残差构造误差泛函；
利用所述波场残差反传播到所述声学参数初始模型空间，得到残差反传播数据；
利用所述正演模拟数据和所述残差反传播数据分别计算所述误差泛函的速度梯度、密度梯度、波阻抗梯度、衰减因子梯度、拉梅参数梯度；
分别利用所述速度梯度、密度梯度、波阻抗梯度、衰减因子梯度、拉梅参数梯度更新所述声学参数初始模型，得到每个声学参数对应的精确模型；
根据每个声学参数对应的精确模型生成所述山体探测区域的山体滑坡危险性评价信息。
进一步的，所述基于所述声学参数初始模型对震源波形进行正演，得到正演模拟数据具体包括:
利用交错网格有限差分法对所述震源波形的波动方程进行时间域离散,得到离散后的波动方程；
根据所述离散后的波动方程和所述声学参数初始模型确定交错网格各个时刻空间分布的波场值。
进一步的，所述根据所述波场残差构造误差泛函具体包括：
根据公式
计算误差泛函，其中，E(m)为误差泛函，b(m)为线性函数表示正演模拟的结果数据，dobs为观测数据，b(m)-dobs为波场残差，CD为数据协方差矩阵，CM为模型的协方差矩阵，m为声学参数初始模型的模型参数，mprior为先验信息模型，λ为先验信息比重参数。
进一步的，所述利用所述正演模拟数据和所述残差反传播数据分别计算所述误差泛函的速度梯度、密度梯度、波阻抗梯度、衰减因子梯度、拉梅参数梯度具体包括：
根据公式
计算所述误差泛函的速度梯度；其中，
Pf为正演模拟数据，Pb残差反传播数据，ω为频率，VP为速度，k、ρ为初始模型参数，E为误差泛函。
进一步的，所述利用所述正演模拟数据和所述残差反传播数据分别计算所述误差泛函的速度梯度、密度梯度、波阻抗梯度、衰减因子梯度、拉梅参数梯度具体包括：
根据公式
计算所述误差泛函的密度梯度；其中，
Pf为正演模拟数据，Pb残差反传播数据，ω为频率，ρ为密度，k为模量，E为误差泛函，Vp为速度。
进一步的，所述利用所述正演模拟数据和所述残差反传播数据分别计算所述误差泛函的速度梯度、密度梯度、波阻抗梯度、衰减因子梯度、拉梅参数梯度具体包括：
根据公式
计算所述误差泛函的波阻抗梯度；其中，
Pf为正演模拟数据，Pb残差反传播数据，ω为频率，IP为波阻抗，k、ρ为初始模型参数，E为误差泛函。
进一步的，所述利用所述正演模拟数据和所述残差反传播数据分别计算所述误差泛函的速度梯度、密度梯度、波阻抗梯度、衰减因子梯度、拉梅参数梯度具体包括：
根据公式
计算所述误差泛函的衰减因子梯度；其中，
E为误差泛函，Qj为衰减因子，ω为频率，ωr为共振频率，ρ为密度，vj为速度,Pf为正演模拟数据，Pb残差反传播数据。
进一步的，所述利用所述正演模拟数据和所述残差反传播数据分别计算所述误差泛函的速度梯度、密度梯度、波阻抗梯度、衰减因子梯度、拉梅参数梯度具体包括：
根据公式
计算所述误差泛函的拉梅参数梯度，其中，m为初始模型参数，u为正演模拟播场，B为正演算子，Δd为波场残差，B-1t为残差反传播数据。
进一步的，所述分别利用所述速度梯度、密度梯度、波阻抗梯度、衰减因子梯度、拉梅参数梯度更新所述声学参数初始模型，得到每个声学参数对应的精确模型具体包括：
确定最优迭代步长和迭代终止条件；
分别根据所述速度梯度、密度梯度、波阻抗梯度、衰减因子梯度、拉梅参数梯度和所述最优迭代步长对所述声学参数初始模型进行更新；
当所述声学参数初始模型满足所述迭代终止条件时，得到每个声学参数对应的精确模型。
本发明提供的用于山体滑坡危险性评价的声学参数获取方法，获取山体探测区域内震源震动的多个单炮数据，单炮数据包括单炮单道数据和单炮多道数据；利用时窗对多个单炮数据中的直达波、浅层反射波、浅层折射波进行截取，得到观测数据；获取声学参数初始模型，并基于声学参数初始模型对震源波形进行正演，得到正演模拟数据，声学参数包括声波传播速度、介质密度、波阻抗、衰减因子和拉梅参数；根据观测数据和正演模拟数据计算波场残差，并根据波场残差构造误差泛函；利用波场残差反传播到声学参数初始模型空间，得到残差反传播数据；利用正演模拟数据和残差反传播数据分别计算误差泛函的速度梯度、密度梯度、波阻抗梯度、衰减因子梯度、拉梅参数梯度；分别利用速度梯度、密度梯度、波阻抗梯度、衰减因子梯度、拉梅参数梯度更新声学参数初始模型，得到每个声学参数对应的精确模型。本发明提供的方法，根据地质任务需求利用初至以后一段时间内波场的运动学和动力学特征完成波形反演，实现浅层的高精度声学参数建模。
附图说明
图1为本发明实施例提供的用于山体滑坡危险性评价的声学参数获取方法流程图；
图2为本发明实施例提供的高斯函数图；
图3为本发明实施例提供的二维声波各物理量和声学参数的网格配置示意图；
图4为本发明实施例提供的三维声波各物理量和声学参数的网格配置示意图。
具体实施方式
为了使本发明的目的、技术方案和优点更加清楚，下面将结合附图对本发明作进一步地详细描述，显然，所描述的实施例仅仅是本发明一部份实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其它实施例，都属于本发明保护的范围。
本发明针对浅层建模的难点，尤其在速度变化剧烈，低速绕动体发育的情况下，避开基于射线理论的弱点，综合利用初至以后一段时间内的波场信息实现浅层的高精度建模，为中深层建模及偏移成像提供可靠的支撑。
图1为本发明实施例一提供的用于山体滑坡危险性评价的声学参数获取方法流程图。如图1所示，具体包括如下步骤：
步骤101，获取山体探测区域内震源震动的多个单炮数据；
将地震勘探仪器布设于山体探测区域内，每放一炮采集一炮的单炮数据，单炮数据包括单炮单道数据和单炮多道数据。根据具体情况和施工需要，发射多炮以采集多组单炮数据。其中，地震勘探仪器具体为408ULS电缆仪器等。
其中，单炮数据为浅剖数据，利用高斯函数将采集到的单炮多道数据处理成单炮单道数据，取高斯函数的一半覆盖整条测线，然后将对应的离散的点，作为权重系数，将权重系数做归一化，使其和为1。整理成一个超道集，如公式(1)所示，从而将单炮多道数据处理成单炮单道数据。
其中，μ为炮点，i为检波点，σ为高斯窗口因子。
选取高斯函数图(如图2所示)中的x轴正半轴部分，横坐标对应检波点，每个点的系数为高斯系数M，如公式(2)所示。
步骤102，利用时窗对多个单炮数据中的直达波、浅层反射波、浅层折射波进行截取，得到初始观测数据；
具体的，在声波传播过程中加以时窗控制，得到近偏移距在近地表以及中浅层传播的波形信息。利用固定时长或非固定时长的移动时窗截取单炮数据中的直达波、浅层反射波、浅层折射波，得到初始观测数据。
所截取的波场并不特指某一类型的波，它包含了很多波的信息，如直达波、首播、透射波、回折波、折射波等。没有经过反射界面改造，其波形信息中含有该区域内速度异常体的丰富信息，因此本发明利用动态时窗所截取的波场对低速体发育的区域进行高精度速度建模。
在利用时窗对单炮数据中的直达波、浅层反射波、浅层折射波进行截取，得到初始观测数据之后，还包括对初始观测数据进行预处理的步骤，具体包括：利用小波变换对初始观测数据进行滤波处理，得到处理后的初始观测数据；对处理后的初始观测数据进行多次波压制处理，得到观测数据。
具体的，利用小波变换对初始观测数据进行分频去噪，小波变化可无限细分，彼此正交，对含有相干干扰的初始观测数据进行分频，去噪处理可仅限于很窄的频带中进行，使去噪后对有效波的损失最大限度的减少，不存在傅里叶变换的频泄现象。
采用共中心点叠加法对处理后的观测数据进行多次波压制处理。共中心点叠加法是依据动校正后一次波和多次波之间剩余时差的差异，将不同接收点收到的来自地下同一反射点的不同激发点的信号，经动校正后叠加起来，进行多次波压制。用一次波的速度作动校正，一次波同相轴被校平而多次波仍有剩余时差，通过叠加使一次波得到增强而多次波得到削弱。
另外，还可以采用倾角滤波、速度滤波、扇形滤波等二维滤波法进行多次波压制，滤除多次波保留一次波。
步骤103，获取声学参数初始模型，并基于声学参数初始模型对震源波形进行正演，得到正演模拟数据；
其中，声学参数包括声波传播速度、介质密度、波阻抗、衰减因子和拉梅参数。
具体的，对地下进行网格化剖分，模型的大小是m*n，m表示水平方向的网格点数，n表示纵向的网格点数，网格间距为h，声学参数初始模型的尺寸为水平m*h米，纵向是n*h米。声学参数初始模型的矩阵形式如公式(3)所示:
声波的波动方程如公式(4)所示:
其中,p为压力场,vx和vz分别为横向和纵向速度场；k＝ρv2；ts为声波传播时间。
在地震勘探领域，声波方程交错网格的思想是把不同的波场分量和地下声学参数配置在不同的网格节点上，时间步进也是采用时间交错步进的。本发明采用的二维声波各物理量和声学参数的网格配置如图3所示，三维声波各物理量和声学参数的网格配置如图4所示。
在声波方程离散前，首先需要采用泰勒展开法推导声波场的空间导数在规则网格和交错网格的高精度近似。
设u(x)有2N+1阶导数，则u(x)在x0+△x,与x0-△x处的2N+1阶泰勒展开式为：
(5)、(6)两式做差得到公式(7)：
同理有：
则任意2N阶精度中心有限差分系数的计算公式为：
令
则有：
化简得：
同理可推导交错网格任意2N阶精度有限差分差分格式和差分系数计算式。
根据泰勒展开式有：
化简得：
其中，an为差分系数
因此采用交错网格有限差分，三维一阶速度-应力声波方程可以离散为：
其中，Δx,Δy,Δz,Δt分别为空间和时间采样间隔，p为每个时刻的应力波场，v表示每个时刻的位移波场，x,y,z分别表示不同的方向，ρ表示密度，f表示震源函数。
当给定初值和边值条件时，用上述差分格式可以递推得到各个时刻波场的空间分布。
在有限差分波动方程正演模拟过程中，为了避免数值噪声以及不稳定，对于一个给定的频带宽度，有限差分的网格大小和时间步长需要分别满足频散关系以及稳定性条件。本发明技术方案采用有限差分频散关系需要满足每个最小波长至少需要5个网格点，即为了避免网格色散空间采样间隔需要满足公式(29):
其中Δx为空间网格大小，λmin是最短波长，
为最小纵波速度，fmax是最大频率。
空间采样网格大小确定之后，还需要选择合适的时间采样大小来满足有限差分数值稳定性条件：
其中Δt为时间采样间隔，
为最大纵波速度。
步骤104，根据观测数据和正演模拟数据计算波场残差，并根据波场残差构造误差泛函；
具体的，将观测数据与正演模拟数据做差，得到波场残差。观测数据与正演模拟数据均为二维的数组，将二维数组内对应的点上的数据相减，得到波场残差。
根据公式(31)构造误差泛函：
其中，E(m)为误差泛函，b(m)为线性函数表示正演模拟的结果数据，dobs为观测数据，b(m)-dobs为波场残差，CD为数据协方差矩阵，CM为模型的协方差矩阵，m为声学参数初始模型的模型参数，mprior为先验信息模型，λ为先验信息比重参数，用来调节模型项和先验信息项的比重。
步骤105，利用波场残差反传播到声学参数初始模型空间，得到残差反传播数据；
将反传算子作用于波场残差,得到初始速度模型空间的残差反传播数据。
步骤104得到检波点所在位置处的波场残差，将这些点处的波场残差作为震源加载到时间域正演模拟过程中，进行时间反方向传播,得到残差反传播数据。
步骤106，利用正演模拟数据和残差反传播数据分别计算误差泛函的速度梯度、密度梯度、波阻抗梯度、衰减因子梯度、拉梅参数梯度；
梯度计算是参数反演的关键部分,代表着模型的更新方向，梯度导引类反演方法通过目标泛函对模型参数的导数,来寻找迭代更新方向,以实现对模型的更新。本发明基于伴随状态法进行梯度求取，将正传波场和反传波场的数据残差作为新的震源进行正演,以求取目标函数对模型的梯度。
具体的，根据公式(32)计算误差泛函的速度梯度；
其中，
Pf为正演模拟数据，Pb残差反传播数据，ω为频率，VP为速度，k、ρ为初始模型参数，E为误差泛函。
根据公式(36)计算误差泛函的密度梯度；
其中，
Pf为正演模拟数据，Pb残差反传播数据，ω为频率，ρ为密度，k为模量，E为误差泛函，Vp为速度。
根据公式(40)计算误差泛函的波阻抗梯度；
其中，
Pf为正演模拟数据，Pb残差反传播数据，ω为频率，IP为波阻抗，k、ρ为初始模型参数，E为误差泛函。
根据公式(44)计算误差泛函的衰减因子梯度；
其中，
E为误差泛函，Qj为衰减因子，ω为频率，ωr为共振频率，ρ为密度，vj为速度,Pf为正演模拟数据，Pb残差反传播数据。
根据公式(46)计算误差泛函的拉梅参数梯度；
其中，m为初始模型参数，u为正演模拟播场，B为正演算子，Δd为波场残差，B-1t为残差反传播数据。
步骤107，分别利用速度梯度、密度梯度、波阻抗梯度、衰减因子梯度、拉梅参数梯度更新声学参数初始模型，得到每个声学参数对应的精确模型。
具体的，令梯度函数
得到扰动模型，最终的精确模型为初始模型和扰动模型之和。
根据最快下降法、共轭梯度法、拟牛顿法LBFGS确定最优迭代步长。分别根据速度梯度、密度梯度、波阻抗梯度、衰减因子梯度、拉梅参数梯度和最优迭代步长对声学参数初始模型进行更新；当声学参数初始模型满足迭代终止条件时，得到每个声学参数对应的精确模型。
利用公式(48)进行模型更新：
mupdate＝mbefore-αkdk (48)
其中，αk为第k步的最优迭代步长，dk为第k步模型的梯度。
迭代终止条件如公式(49)所示：
公式(9)规定了声学参数m迭代收敛规则，第n次更新的量小于n-1次迭代的一定比例，例如，ε＝0.001，则迭代收敛，反演终止；否则，将更新结果作为输入，进行下一次迭代。
将每个声学参数对应的精确模型分别与对应的标准声学参数进行比对，根据比对结果确定山体滑坡的危险程度，根据危险程度生成山体探测区域的山体滑坡危险性评价信息，从而评估山体滑坡的可能性和进行预警。
本发明提供的用于山体滑坡危险性评价的声学参数获取方法，获取山体探测区域内震源震动的多个单炮数据，单炮数据包括单炮单道数据和单炮多道数据；利用时窗对多个单炮数据中的直达波、浅层反射波、浅层折射波进行截取，得到观测数据；获取声学参数初始模型，并基于声学参数初始模型对震源波形进行正演，得到正演模拟数据，声学参数包括声波传播速度、介质密度、波阻抗、衰减因子和拉梅参数；根据观测数据和正演模拟数据计算波场残差，并根据波场残差构造误差泛函；利用波场残差反传播到声学参数初始模型空间，得到残差反传播数据；利用正演模拟数据和残差反传播数据分别计算误差泛函的速度梯度、密度梯度、波阻抗梯度、衰减因子梯度、拉梅参数梯度；分别利用速度梯度、密度梯度、波阻抗梯度、衰减因子梯度、拉梅参数梯度更新声学参数初始模型，得到每个声学参数对应的精确模型。本发明提供的方法，根据地质任务需求利用初至以后一段时间内波场的运动学和动力学特征完成波形反演，实现浅层的高精度声学参数建模。
专业人员应该还可以进一步意识到，结合本文中所公开的实施例描述的各示例的单元及算法步骤，能够以电子硬件、计算机软件或者二者的结合来实现，为了清楚地说明硬件和软件的可互换性，在上述说明中已经按照功能一般性地描述了各示例的组成及步骤。这些功能究竟以硬件还是软件方式来执行，取决于技术方案的特定应用和设计约束条件。专业技术人员可以对每个特定的应用来使用不同方法来实现所描述的功能，但是这种实现不应认为超出本发明的范围。
结合本文中所公开的实施例描述的方法或算法的步骤可以用硬件、处理器执行的软件模块，或者二者的结合来实施。软件模块可以置于随机存储器(RAM)、内存、只读存储器(ROM)、电可编程ROM、电可擦除可编程ROM、寄存器、硬盘、可移动磁盘、CD-ROM、或技术领域内所公知的任意其它形式的存储介质中。
以上所述的具体实施方式，对本发明的目的、技术方案和有益效果进行了进一步详细说明，所应理解的是，以上所述仅为本发明的具体实施方式而已，并不用于限定本发明的保护范围，凡在本发明的精神和原则之内，所做的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。
Claims (9)
1.一种用于山体滑坡危险性评价的声学参数获取方法，其特征在于，所述方法包括：
获取山体探测区域内震源震动的多个单炮数据，所述单炮数据包括单炮单道数据和单炮多道数据；
利用时窗对所述多个单炮数据中的直达波、浅层反射波、浅层折射波进行截取，得到初始观测数据；
对所述初始观测数据进行滤波处理和多次波压制处理，得到观测数据；
获取声学参数初始模型，并基于所述声学参数初始模型对震源波形进行正演，得到正演模拟数据，所述声学参数包括声波传播速度、介质密度、波阻抗、衰减因子和拉梅参数；
根据所述观测数据和所述正演模拟数据计算波场残差，并根据所述波场残差构造误差泛函；
利用所述波场残差反传播到所述声学参数初始模型空间，得到残差反传播数据；
利用所述正演模拟数据和所述残差反传播数据分别计算所述误差泛函的速度梯度、密度梯度、波阻抗梯度、衰减因子梯度、拉梅参数梯度；
分别利用所述速度梯度、密度梯度、波阻抗梯度、衰减因子梯度、拉梅参数梯度更新所述声学参数初始模型，得到每个声学参数对应的精确模型；
根据每个声学参数对应的精确模型生成所述山体探测区域的山体滑坡危险性评价信息。
2.根据权利要求1所述的方法，其特征在于，所述基于所述声学参数初始模型对震源波形进行正演，得到正演模拟数据具体包括:
利用交错网格有限差分法对所述震源波形的波动方程进行时间域离散,得到离散后的波动方程；
根据所述离散后的波动方程和所述声学参数初始模型确定交错网格各个时刻空间分布的波场值。
3.根据权利要求1所述的方法，其特征在于，所述根据所述波场残差构造误差泛函具体包括：
根据公式
计算误差泛函，其中，E(m)为误差泛函，b(m)为线性函数表示正演模拟的结果数据，dobs为观测数据，b(m)-dobs为波场残差，CD为数据协方差矩阵，CM为模型的协方差矩阵，m为声学参数初始模型的模型参数，mprior为先验信息模型，λ为先验信息比重参数。
4.根据权利要求1所述的方法，其特征在于，所述利用所述正演模拟数据和所述残差反传播数据分别计算所述误差泛函的速度梯度、密度梯度、波阻抗梯度、衰减因子梯度、拉梅参数梯度具体包括：
根据公式
计算所述误差泛函的速度梯度；其中，
k＝ρVp 2，Pf为正演模拟数据，Pb残差反传播数据，ω为频率，VP为速度，k、ρ为初始模型参数，E为误差泛函。
5.根据权利要求1所述的方法，其特征在于，所述利用所述正演模拟数据和所述残差反传播数据分别计算所述误差泛函的速度梯度、密度梯度、波阻抗梯度、衰减因子梯度、拉梅参数梯度具体包括：
根据公式
计算所述误差泛函的密度梯度；其中，
k＝ρVp 2；Pf为正演模拟数据，Pb残差反传播数据，ω为频率，ρ为密度，k为模量，E为误差泛函，Vp为速度。
6.根据权利要求1所述的方法，其特征在于，所述利用所述正演模拟数据和所述残差反传播数据分别计算所述误差泛函的速度梯度、密度梯度、波阻抗梯度、衰减因子梯度、拉梅参数梯度具体包括：
根据公式
计算所述误差泛函的波阻抗梯度；其中，
Pf为正演模拟数据，Pb残差反传播数据，ω为频率，IP为波阻抗，k、ρ为初始模型参数，E为误差泛函。
7.根据权利要求1所述的方法，其特征在于，所述利用所述正演模拟数据和所述残差反传播数据分别计算所述误差泛函的速度梯度、密度梯度、波阻抗梯度、衰减因子梯度、拉梅参数梯度具体包括：
根据公式
计算所述误差泛函的衰减因子梯度；其中，
E为误差泛函，Qj为衰减因子，ω为频率，ωr为共振频率，ρ为密度，vj为速度,Pf为正演模拟数据，Pb残差反传播数据。
8.根据权利要求1所述的方法，其特征在于，所述利用所述正演模拟数据和所述残差反传播数据分别计算所述误差泛函的速度梯度、密度梯度、波阻抗梯度、衰减因子梯度、拉梅参数梯度具体包括：
根据公式
计算所述误差泛函的拉梅参数梯度，其中，m为初始模型参数，u为正演模拟播场，B为正演算子，Δd为波场残差，B-1t为残差反传播数据。
9.根据权利要求1所述的方法，其特征在于，所述分别利用所述速度梯度、密度梯度、波阻抗梯度、衰减因子梯度、拉梅参数梯度更新所述声学参数初始模型，得到每个声学参数对应的精确模型具体包括：
确定最优迭代步长和迭代终止条件；
分别根据所述速度梯度、密度梯度、波阻抗梯度、衰减因子梯度、拉梅参数梯度和所述最优迭代步长对所述声学参数初始模型进行更新；
当所述声学参数初始模型满足所述迭代终止条件时，得到每个声学参数对应的精确模型。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201910172714.3A CN111665552A (zh) | 2019-03-07 | 2019-03-07 | 用于山体滑坡危险性评价的声学参数获取方法 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201910172714.3A CN111665552A (zh) | 2019-03-07 | 2019-03-07 | 用于山体滑坡危险性评价的声学参数获取方法 |
Publications (1)
| Publication Number | Publication Date |
|---|---|
| CN111665552A true CN111665552A (zh) | 2020-09-15 |
Family
ID=72382219
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN201910172714.3A Pending CN111665552A (zh) | 2019-03-07 | 2019-03-07 | 用于山体滑坡危险性评价的声学参数获取方法 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN111665552A (zh) |
Citations (9)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN101598809A (zh) * | 2008-06-04 | 2009-12-09 | 中国石油天然气集团公司 | 一种自适应消除线性规则噪声以及多次波干扰的方法 |
| CN103091711A (zh) * | 2013-01-24 | 2013-05-08 | 中国石油天然气集团公司 | 全波形反演方法及装置 |
| CN103630933A (zh) * | 2013-12-09 | 2014-03-12 | 中国石油天然气集团公司 | 基于非线性优化的时空域交错网格有限差分方法和装置 |
| WO2017034433A1 (en) * | 2015-08-25 | 2017-03-02 | Saudi Arabian Oil Company | Three-dimensional elastic frequency-domain iterative solver for full waveform inversion |
| CN106662664A (zh) * | 2014-06-17 | 2017-05-10 | 埃克森美孚上游研究公司 | 快速粘声波和粘弹性全波场反演 |
| CN106842295A (zh) * | 2015-12-04 | 2017-06-13 | 中国石油化工股份有限公司 | 测井信息约束的波形反演方法 |
| CN108181652A (zh) * | 2017-12-29 | 2018-06-19 | 中国海洋大学 | 一种海底节点地震资料上下行波场数值模拟方法 |
| CN108680964A (zh) * | 2018-03-30 | 2018-10-19 | 吉林大学 | 一种基于结构约束的归一化重磁电震联合反演方法 |
| US20180356548A1 (en) * | 2017-06-12 | 2018-12-13 | Institute Of Geology And Geophysics Chinese Academy Of Sciences | Inversion velocity model, method for establishing the same and method for acquiring images of underground structure |
-
2019
- 2019-03-07 CN CN201910172714.3A patent/CN111665552A/zh active Pending
Patent Citations (9)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN101598809A (zh) * | 2008-06-04 | 2009-12-09 | 中国石油天然气集团公司 | 一种自适应消除线性规则噪声以及多次波干扰的方法 |
| CN103091711A (zh) * | 2013-01-24 | 2013-05-08 | 中国石油天然气集团公司 | 全波形反演方法及装置 |
| CN103630933A (zh) * | 2013-12-09 | 2014-03-12 | 中国石油天然气集团公司 | 基于非线性优化的时空域交错网格有限差分方法和装置 |
| CN106662664A (zh) * | 2014-06-17 | 2017-05-10 | 埃克森美孚上游研究公司 | 快速粘声波和粘弹性全波场反演 |
| WO2017034433A1 (en) * | 2015-08-25 | 2017-03-02 | Saudi Arabian Oil Company | Three-dimensional elastic frequency-domain iterative solver for full waveform inversion |
| CN106842295A (zh) * | 2015-12-04 | 2017-06-13 | 中国石油化工股份有限公司 | 测井信息约束的波形反演方法 |
| US20180356548A1 (en) * | 2017-06-12 | 2018-12-13 | Institute Of Geology And Geophysics Chinese Academy Of Sciences | Inversion velocity model, method for establishing the same and method for acquiring images of underground structure |
| CN108181652A (zh) * | 2017-12-29 | 2018-06-19 | 中国海洋大学 | 一种海底节点地震资料上下行波场数值模拟方法 |
| CN108680964A (zh) * | 2018-03-30 | 2018-10-19 | 吉林大学 | 一种基于结构约束的归一化重磁电震联合反演方法 |
Non-Patent Citations (16)
| Title |
|---|
| ALBERT TARANTOLA: "Inversion of seismic reflection data in the acoustic approximation", 《GEOPHYSICS》 * |
| ALBERT TARANTOLA: "Inversion of seismic reflection data in the acoustic approximation", 《GEOPHYSICS》, vol. 49, no. 8, 31 August 1984 (1984-08-31), pages 1259 - 1266, XP009173614, DOI: 10.1190/1.1441754 * |
| J. VIRIEUX ET AL.: "An overview of full-waveform inversion in exploration geophysics", 《GEOPHYSICS》 * |
| J. VIRIEUX ET AL.: "An overview of full-waveform inversion in exploration geophysics", 《GEOPHYSICS》, vol. 74, no. 6, 31 December 2009 (2009-12-31), pages 1 * |
| 丰超: "波形反演应用初探", 中国优秀硕士学位论文全文数据库 基础科学辑, pages 7 - 12 * |
| 于长澎: "地层吸收系数的全波形反演研究", 《中国优秀硕士学位论文全文数据库 基础科学辑》 * |
| 于长澎: "地层吸收系数的全波形反演研究", 《中国优秀硕士学位论文全文数据库 基础科学辑》, no. 1, 15 January 2012 (2012-01-15), pages 45 - 55 * |
| 冯波等: "《数据域波动方程层析速度反演方法研究》", 31 August 2017, 上海：同济大学出版社, pages: 58 * |
| 刘玉柱等: "有效利用初至信息的偏移距加权地震层析成像方法", 《石油物探》 * |
| 刘玉柱等: "有效利用初至信息的偏移距加权地震层析成像方法", 《石油物探》, vol. 53, no. 1, 31 January 2014 (2014-01-31), pages 99 - 105 * |
| 吴从辉等: "任意偶数阶精度交错网格声波方程数值模拟", 东华理工大学学报(自然科学版), vol. 38, no. 1, pages 103 - 104 * |
| 智敏等: "一阶声波方程高阶交错网格有限差分数值模拟方法", 科技论坛, pages 84 - 85 * |
| 朱瑞: "基于逆时偏移的地震反射波成像研究", 中国优秀硕士学位论文全文数据库 基础科学辑, pages 12 - 14 * |
| 杨积忠等: "变密度声波方程多参数全波形反演策略", 《地球物理学报》 * |
| 杨积忠等: "变密度声波方程多参数全波形反演策略", 《地球物理学报》, vol. 57, no. 2, 28 February 2014 (2014-02-28), pages 628 - 643 * |
| 郝奇: "三维TTI介质地震波场正演模拟", 中国优秀硕士学位论文全文数据库 基础科学辑, no. 04, pages 17 - 21 * |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| KR101548976B1 (ko) | 지진 표면파들의 파형들을 사용하는 토양 특성들의 추정 | |
| CA2690373C (en) | Method for velocity analysis using waveform inversion in laplace domain for geophysical imaging | |
| RU2518577C2 (ru) | Непрерывный адаптивный анализ поверхностных волн в случае трехмерных сейсмических данных | |
| EP0799429B1 (en) | Source signature determination and multiple reflection reduction | |
| MX2011006036A (es) | Uso de inversion de forma de onda para determinar las propiedades de un medio en el subsuelo. | |
| CN111665556B (zh) | 地层声波传播速度模型构建方法 | |
| KR101820850B1 (ko) | 직접 파형 역산의 반복 적용을 이용한 탄성파 영상화 장치 및 방법 | |
| Poppeliers et al. | Three‐dimensional seismic‐wave gradiometry for scalar waves | |
| CN113466933B (zh) | 基于深度加权的地震斜率层析成像方法 | |
| Park et al. | Refraction traveltime tomography based on damped wave equation for irregular topographic model | |
| Adinolfi et al. | A functional tool to explore the reliability of micro-earthquake focal mechanism solutions for seismotectonic purposes | |
| CN111665546B (zh) | 用于可燃冰探测的声学参数获取方法 | |
| Shen et al. | Characteristics of the far-field source signature and its ghost effect of a marine single air gun | |
| Zheng et al. | Double-beam stacking to infer seismic properties of fractured reservoirs | |
| CN114721044B (zh) | 一种多频率接收函数和振幅比联合反演地壳结构的方法及系统 | |
| CN111665550A (zh) | 地下介质密度信息反演方法 | |
| CN111665549A (zh) | 地层声波衰减因子反演方法 | |
| NO346771B1 (en) | Methods and systems to separate wavefields using pressure wavefield data | |
| CN111665551B (zh) | 用于桥梁基底探测的声学参数获取方法 | |
| CN111665552A (zh) | 用于山体滑坡危险性评价的声学参数获取方法 | |
| CN111665544A (zh) | 用于地下采空区探测的声学参数获取方法 | |
| CN111665548A (zh) | 用于海底探测的声学参数获取方法 | |
| CN111665545A (zh) | 用于稀有金属探测的声学参数获取方法 | |
| CN111665543A (zh) | 用于地铁危险性检测的声学参数获取方法 | |
| CN111665554A (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 | ||
| RJ01 | Rejection of invention patent application after publication | ||
| RJ01 | Rejection of invention patent application after publication |
Application publication date: 20200915 |