CN111045076A - 多模式瑞雷波频散曲线并行联合反演方法 - Google Patents

多模式瑞雷波频散曲线并行联合反演方法 Download PDF

Info

Publication number
CN111045076A
CN111045076A CN201911259895.XA CN201911259895A CN111045076A CN 111045076 A CN111045076 A CN 111045076A CN 201911259895 A CN201911259895 A CN 201911259895A CN 111045076 A CN111045076 A CN 111045076A
Authority
CN
China
Prior art keywords
dispersion curve
inversion
order
formula
rayleigh wave
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
Application number
CN201911259895.XA
Other languages
English (en)
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.)
Beijing Research Institute of Uranium Geology
Original Assignee
Beijing Research Institute of Uranium Geology
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 Beijing Research Institute of Uranium Geology filed Critical Beijing Research Institute of Uranium Geology
Priority to CN201911259895.XA priority Critical patent/CN111045076A/zh
Publication of CN111045076A publication Critical patent/CN111045076A/zh
Pending legal-status Critical Current

Links

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
    • 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/62—Physical property of subsurface
    • G01V2210/624—Reservoir parameters

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

本发明属于近地表地震面波勘探技术领域，具体涉及一种瑞雷波多模式频散曲线并行反演方法。本发明包括以下步骤：步骤一、对实际采集的地震记录做二维傅里叶变换，得到实测瑞雷波频散曲线f‑vo Rij；步骤二、证明频散曲线与各阶综合灵敏度关系，频散曲线与各阶综合灵敏度的最高值相对应；步骤三、基于步骤二的结论，对步骤一所求取的实测瑞雷波频散曲线f‑vo Rij，设计反演个体目标函数Φ；步骤四、基于步骤三反演个体目标函数利用阻尼最小二乘法进行迭代优化，更新模型x；步骤五、基于步骤四所产生的新模型，重复步骤三、四，直至达到事先设定的迭代次数T，输出地层模型反演结果x。本发明能够明显提高反演收敛效率和地层模型参数精度。

Description

多模式瑞雷波频散曲线并行联合反演方法
技术领域
本发明属于近地表地震面波勘探技术领域，具体涉及一种瑞雷波多模式频散曲线并行反演方法。
背景技术
多道瞬态面波法利用瑞雷波在弹性介质(土壤、岩石、路基)中传播的频散特性，通过提取、反演频散曲线构建近地表S波速度结构模型，具有快速、便捷、非侵入性及较高的浅层分辨率等优点，被广泛应用于地下水、工程地质及环境等探查研究中。瑞雷波在介质界面附近大约0.5个波长深度范围内从左至右以逆进椭圆形式向前传播，其能量不均匀地分布于各个模式，不同模式的频散在各频带内均可能占据主导地位。研究表明，高模式的瑞雷波具有更大的探测深度和反演敏感度，而常规单模式反演受其有限的敏感频带限制，模型各参数并不能有效地向真实参数同步收敛，为了对地层参数进一步约束，瑞雷波联合体波、电法、重力甚至勒夫波等众多地球物理探测手段进行联合反演的方法技术层出不穷，而多模式频散联合反演不失为一种更有效、更经济的探测方法。
传统从基阶到高阶逐步递进的联合反演不仅收敛速度缓慢，而且在处理局部变化梯度较大的实测频散时，各阶频散反演难以相互制约以达到一种平衡状态，使得真实地球物理模型和反演模型并不能完全契合。
发明内容
本发明解决的技术问题：
本发明提供一种瑞雷波多模式频散曲线并行反演方法，定义了频散曲线综合灵敏度函数，进一步证明了其与频谱图中视频散曲线的对应关系，并基于此定义了多模式频散曲线并行反演目标函数，对反演地层模型实行多模式频散的交叉梯度约束可以明显提高反演收敛效率和地层模型参数精度。
为解决上述技术问题，本发明一种瑞雷波多模式频散曲线并行反演方法，包括以下步骤：
步骤一、对实际采集的地震记录做二维傅里叶变换，将地震记录从x-t域转换到f-k域；然后依据k＝f/v生成f-v频谱能量图，对其能量峰值进行识别，得到实测瑞雷波频散曲线f-vo Rij，vo Rij是实测频散曲线，其中i＝1,2,…N，N是实测相速度的个数，j＝1,2,…M，M是拾取的瑞雷波频散曲线的阶数；
步骤二、证明频散曲线与各阶综合灵敏度关系，频散曲线与各阶综合灵敏度的最高值相对应；
步骤三、基于步骤二的结论，对步骤一所求取的实测瑞雷波频散曲线f-vo Rij，设计反演个体目标函数Φ；
步骤四、基于步骤三反演个体目标函数利用阻尼最小二乘法进行迭代优化，并记迭代次数T＝1，其参数修正量求取公式为：
Figure BDA0002311336280000021
其中，I是单位矩阵，x0是反演初始模型参数，μ是阻尼因子，本实施例设其初始值为1，A是雅格比矩阵，通过步骤二灵敏度公式求取；
基于上述求得的参数修正量Δx，利用最速下降法求取参数最优修正步长λ；求取最优修正步长λ后，按照这个公式更新模型x，
则当前所选个体参数为x＝x0+λΔx，并令T＝T+1；
步骤五、基于步骤四所产生的新模型，重复步骤三、四，直至达到事先设定的迭代次数T，输出地层模型反演结果x。
所述步骤一中，从x-t域转换到f-k域的求取公式为：
Figure BDA0002311336280000031
其中，f是频率，u(x，t)是x-t域地震记录；
Figure BDA0002311336280000032
其中，k是波数，U(f，k)是f-k域地震记录。
所述步骤二的具体步骤为，
计算地层模型各层横波速度vSk对于某j阶频散曲线的敏感度；
定义综合灵敏度函数，即对相同阶SvSk(f，j)进行求和；
根据上式计算的各阶频散综合灵敏度曲线与频谱能量图的对比，证明频散曲线与各阶综合灵敏度的最高值相对应。
所述步骤二中，敏感度求取公式为：
Figure BDA0002311336280000033
其中SvSk(f，j)为第k层横波速度vSk在频率f时计算的第j阶频散灵敏度，c(f，j)是频率f时计算的第j阶瑞雷波相速度，α为差商因子。
所述相同阶SvSk(f，j)进行求和公式为
Figure BDA0002311336280000034
其中，k＝1,2,…L，L是地层模型的层数。
所述步骤三中的求取公式为：
Figure BDA0002311336280000035
其中，vo Rij是实测瑞雷波相速度，vc Rij是反演的理论相速度值。
所述步骤四找那个参数最优修正步长λ的求取公式为：
Figure BDA0002311336280000041
式中，gi＝(ATA+μI)Δxi。
本发明的有益技术效果在于：
(1)本发明提供的一种瑞雷波多模式频散曲线并行反演方法，加入的高阶瑞雷波频散曲线信息增强了对地层参数的约束作用，基于多模频散交叉梯度进行的修正使得反演各参数在整个频带范围内有效地向真实模型同步收敛；
(2)本发明提供的一种瑞雷波多模式频散曲线并行反演方法，在一定程度上降低了阻尼最小二乘法对反演初始模型的依赖程度，并在反演效率和反演精度上大大提高，
(3)本发明提供的一种瑞雷波多模式频散曲线并行反演方法，对于多道瞬态面波近地表精细结构探测具有现实意义。
附图说明
图1为实际采集的瑞雷波地震记录；
图2为地震记录频散能量图及提取的多阶瑞雷波频散曲线；
图3为频散能量图及各阶频散综合灵敏度值；
图4为本方法地层模型迭代反演最终结果。
具体实施方式
下面结合附图和实施实例对本发明提供的一种瑞雷波多模式频散曲线并行反演方法作进一步详细说明。
本发明为一种瑞雷波多模式频散曲线并行反演方法，具体包括以下步骤：
步骤一、对实际采集的如图1所示的地震记录做二维傅里叶变换，将地震记录从x-t域转换到f-k域，其求取公式为：
Figure BDA0002311336280000051
其中，f是频率，u(x，t)是x-t域地震记录；
Figure BDA0002311336280000052
其中，k是波数，U(f，k)是f-k域地震记录。
然后依据k＝f/v生成如图2所示的f-v频谱能量图，对其能量峰值进行识别，得到实测瑞雷波频散曲线f-vo Rij，vo Rij是实测频散曲线，其中i＝1,2,…N，N是实测相速度的个数，j＝1,2,…M，M是拾取的瑞雷波频散曲线的阶数。
步骤二、计算地层模型各层横波速度vSk对于某j阶频散曲线的敏感度，其求取公式为：
Figure BDA0002311336280000053
其中SvSk(f，j)为第k层横波速度vSk在频率f时计算的第j阶频散灵敏度，c(f，j)是频率f时计算的第j阶瑞雷波相速度，α为差商因子，本实施例取1.2。
为了综合分析反演参数vS对各阶频散的灵敏度值SvSj，本实施例定义了综合灵敏度函数，即对相同阶SvSk(f，j)进行求和：
Figure BDA0002311336280000054
其中，k＝1,2,…L，L是地层模型的层数，图3是根据上式计算的各阶频散综合灵敏度曲线(虚线，与右侧纵坐标对应)与频谱能量图(能量峰值即为各阶频散曲线，与左侧纵坐标对应)的对比，可以看出频谱图中出现的频散曲线(即称为视频散曲线)与各阶综合灵敏度的最高值相对应。
步骤三、基于步骤二的结论，对步骤一所求取的实测瑞雷波频散曲线f-vo Rij，设计反演个体目标函数Φ，其求取公式为：
Figure BDA0002311336280000061
其中，vo Rij是实测瑞雷波相速度，vc Rij是反演的理论相速度值；
步骤四、基于步骤三反演个体目标函数利用阻尼最小二乘法进行迭代优化，并记迭代次数T＝1，其参数修正量求取公式为：
Figure BDA0002311336280000062
其中，I是单位矩阵，x0是反演初始模型参数，μ是阻尼因子，本实施例设其初始值为1，A是雅格比矩阵，通过步骤二灵敏度公式求取；
基于上述求得的参数修正量Δx，利用最速下降法求取参数最优修正步长λ，其求取公式为：
Figure BDA0002311336280000063
式中，gi＝(ATA+μI)Δxi。求取最优修正步长λ后，按照这个公式更新模型x，则当前所选个体参数为x＝x0+λΔx，并令T＝T+1；
步骤五、基于步骤四所产生的新模型，重复步骤三、四，直至达到事先设定的迭代次数T，输出地层模型反演结果x如图4。
图4所示是本方法耗时38.735s，12次迭代优化后的地层模型，目标函数拟合误差为4.48e-14km/s，可以看出虽然初始模型误差较大且没有反映出地层的结构信息，在各模式交叉梯度的约束下，各参数迅速向真实地层收敛，反演结果和真实模型完全重合，说明本方法的有效性及准确性。

Claims (7)

1.一种瑞雷波多模式频散曲线并行反演方法，其特征在于：包括以下步骤：
步骤一、对实际采集的地震记录做二维傅里叶变换，将地震记录从x-t域转换到f-k域；然后依据k＝f/v生成f-v频谱能量图，对其能量峰值进行识别，得到实测瑞雷波频散曲线f-vo Rij，vo Rij是实测频散曲线，其中i＝1,2,…N，N是实测相速度的个数，j＝1,2,…M，M是拾取的瑞雷波频散曲线的阶数；
步骤二、证明频散曲线与各阶综合灵敏度关系，频散曲线与各阶综合灵敏度的最高值相对应；
步骤三、基于步骤二的结论，对步骤一所求取的实测瑞雷波频散曲线f-vo Rij，设计反演个体目标函数Φ；
步骤四、基于步骤三反演个体目标函数利用阻尼最小二乘法进行迭代优化，并记迭代次数T＝1，其参数修正量求取公式为：
Figure FDA0002311336270000011
其中，I是单位矩阵，x0是反演初始模型参数，μ是阻尼因子，本实施例设其初始值为1，A是雅格比矩阵，通过步骤二灵敏度公式求取；
基于上述求得的参数修正量Δx，利用最速下降法求取参数最优修正步长λ；求取最优修正步长λ后，按照这个公式更新模型x，
则当前所选个体参数为x＝x0+λΔx，并令T＝T+1；
步骤五、基于步骤四所产生的新模型，重复步骤三、四，直至达到事先设定的迭代次数T，输出地层模型反演结果x。
2.根据权利要求1所述的一种瑞雷波多模式频散曲线并行反演方法，其特征在于：所述步骤一中，从x-t域转换到f-k域的求取公式为：
Figure FDA0002311336270000021
其中，f是频率，u(x，t)是x-t域地震记录；
Figure FDA0002311336270000022
其中，k是波数，U(f，k)是f-k域地震记录。
3.根据权利要求2所述的一种瑞雷波多模式频散曲线并行反演方法，其特征在于：所述步骤二的具体步骤为，
计算地层模型各层横波速度vSk对于某j阶频散曲线的敏感度；
定义综合灵敏度函数，即对相同阶SvSk(f，j)进行求和；
根据上式计算的各阶频散综合灵敏度曲线与频谱能量图的对比，证明频散曲线与各阶综合灵敏度的最高值相对应。
4.根据权利要求3所述的一种瑞雷波多模式频散曲线并行反演方法，其特征在于：所述步骤二中，敏感度求取公式为：
Figure FDA0002311336270000023
其中SvSk(f，j)为第k层横波速度vSk在频率f时计算的第j阶频散灵敏度，c(f，j)是频率f时计算的第j阶瑞雷波相速度，α为差商因子。
5.根据权利要求4所述的一种瑞雷波多模式频散曲线并行反演方法，其特征在于：所述相同阶SvSk(f，j)进行求和公式为
Figure FDA0002311336270000024
其中，k＝1,2,…L，L是地层模型的层数。
6.根据权利要求5所述的一种瑞雷波多模式频散曲线并行反演方法，其特征在于：所述步骤三中的求取公式为：
Figure FDA0002311336270000031
其中，vo Rij是实测瑞雷波相速度，vc Rij是反演的理论相速度值。
7.根据权利要求6所述的一种瑞雷波多模式频散曲线并行反演方法，其特征在于：所述步骤四找那个参数最优修正步长λ的求取公式为：
Figure FDA0002311336270000032
式中，gi＝(ATA+μI)Δxi。
CN201911259895.XA 2019-12-10 2019-12-10 多模式瑞雷波频散曲线并行联合反演方法 Pending CN111045076A (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201911259895.XA CN111045076A (zh) 2019-12-10 2019-12-10 多模式瑞雷波频散曲线并行联合反演方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201911259895.XA CN111045076A (zh) 2019-12-10 2019-12-10 多模式瑞雷波频散曲线并行联合反演方法

Publications (1)

Publication Number Publication Date
CN111045076A true CN111045076A (zh) 2020-04-21

Family

ID=70235735

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201911259895.XA Pending CN111045076A (zh) 2019-12-10 2019-12-10 多模式瑞雷波频散曲线并行联合反演方法

Country Status (1)

Country Link
CN (1) CN111045076A (zh)

Cited By (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN112765774A (zh) * 2020-12-25 2021-05-07 青岛黄海学院 一种铁路震源瑞雷面波力学模型及其数值模拟方法
CN112861721A (zh) * 2021-02-09 2021-05-28 南方科技大学 一种自动提取背景噪声频散曲线的方法及装置
CN113671570A (zh) * 2021-08-23 2021-11-19 湖南工商大学 一种地震面波走时和重力异常联合反演方法与***
CN113945975A (zh) * 2021-10-09 2022-01-18 中国船舶重工集团公司第七六0研究所 一种基于勒夫波和瑞利波联合反演地层分层结构的方法
CN114185093A (zh) * 2021-12-07 2022-03-15 中国石油大学(北京) 一种基于瑞雷面波反演的近地表速度模型建立方法及装置

Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN104216003A (zh) * 2014-09-20 2014-12-17 中国地质大学（北京） 多道瞬态瑞雷波探测方法
CN104678435A (zh) * 2014-10-27 2015-06-03 李欣欣 一种提取Rayleigh面波频散曲线的方法
CN107703545A (zh) * 2017-09-01 2018-02-16 中煤科工集团西安研究院有限公司 一种三分量地震槽波波场分离方法及***
CN109239773A (zh) * 2018-09-12 2019-01-18 西安石油大学 一种高阶模式瑞雷波的重建方法
CN109799530A (zh) * 2018-12-25 2019-05-24 核工业北京地质研究院 用于地震面波勘探的瑞雷波频散曲线反演方法

Patent Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN104216003A (zh) * 2014-09-20 2014-12-17 中国地质大学（北京） 多道瞬态瑞雷波探测方法
CN104678435A (zh) * 2014-10-27 2015-06-03 李欣欣 一种提取Rayleigh面波频散曲线的方法
CN107703545A (zh) * 2017-09-01 2018-02-16 中煤科工集团西安研究院有限公司 一种三分量地震槽波波场分离方法及***
CN109239773A (zh) * 2018-09-12 2019-01-18 西安石油大学 一种高阶模式瑞雷波的重建方法
CN109799530A (zh) * 2018-12-25 2019-05-24 核工业北京地质研究院 用于地震面波勘探的瑞雷波频散曲线反演方法

Non-Patent Citations (2)

* Cited by examiner, † Cited by third party
Title
梁青: "多模式瑞雷波频散曲线反演研究", 《中国优秀硕士论文全文数据库 基础科学辑》 *
雷宇航: "自适应GA与DLS联合反演瑞雷波频散曲线方法研究", 《中国优秀硕士论文全文数据库 基础科学辑》 *

Cited By (8)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN112765774A (zh) * 2020-12-25 2021-05-07 青岛黄海学院 一种铁路震源瑞雷面波力学模型及其数值模拟方法
CN112765774B (zh) * 2020-12-25 2022-07-01 青岛黄海学院 一种铁路震源瑞雷面波力学模型及其数值模拟方法
CN112861721A (zh) * 2021-02-09 2021-05-28 南方科技大学 一种自动提取背景噪声频散曲线的方法及装置
CN112861721B (zh) * 2021-02-09 2024-05-07 南方科技大学 一种自动提取背景噪声频散曲线的方法及装置
CN113671570A (zh) * 2021-08-23 2021-11-19 湖南工商大学 一种地震面波走时和重力异常联合反演方法与***
CN113945975A (zh) * 2021-10-09 2022-01-18 中国船舶重工集团公司第七六0研究所 一种基于勒夫波和瑞利波联合反演地层分层结构的方法
CN114185093A (zh) * 2021-12-07 2022-03-15 中国石油大学(北京) 一种基于瑞雷面波反演的近地表速度模型建立方法及装置
CN114185093B (zh) * 2021-12-07 2023-05-12 中国石油大学(北京) 一种基于瑞雷面波反演的近地表速度模型建立方法及装置

Similar Documents

Publication Publication Date Title
CN111045076A (zh) 多模式瑞雷波频散曲线并行联合反演方法
KR101548976B1 (ko) 지진 표면파들의 파형들을 사용하는 토양 특성들의 추정
Bretaudeau et al. 2D elastic full‐waveform imaging of the near‐surface: application to synthetic and physical modelling data sets
O’Neill et al. Dominant higher surface-wave modes and possible inversion pitfalls
CN108873066B (zh) 弹性介质波动方程反射波旅行时反演方法
Pan et al. Love-wave waveform inversion in time domain for shallow shear-wave velocity
CN110780351B (zh) 纵波和转换波叠前联合反演方法及***
Mosher et al. Common angle imaging conditions for pre-stack depth migration
CN111025387B (zh) 一种页岩储层的叠前地震多参数反演方法
CN109188519B (zh) 一种极坐标下的弹性波纵横波速度反演***及方法
Yao et al. An effective absorbing layer for the boundary condition in acoustic seismic wave simulation
US20150309200A1 (en) A method for processing acoustic waveforms
CN109655918B (zh) 地面浅井微地震监测观测台站位置确定方法及***
CN105093278A (zh) 基于激发主能量优化算法的全波形反演梯度算子的提取方法
CN111025388B (zh) 一种多波联合的叠前波形反演方法
CN110579795A (zh) 基于被动源地震波形及其逆时成像的联合速度反演方法
Nguyen et al. Detection of bridge-deck delamination using full ultrasonic waveform tomography
CN110658558A (zh) 吸收衰减介质叠前深度逆时偏移成像方法及***
Kumar et al. Resolving phase wrapping by using sliding transform for generation of dispersion curves
Zhang et al. A niching particle swarm optimization strategy for the multimodal inversion of surface waves
CN112649848B (zh) 利用波动方程求解地震波阻抗的方法和装置
Chen et al. A frequency-domain full waveform inversion method of elastic waves in quantitative defection investigation
Song et al. Application of full waveform inversion to land seismic data in Sichuan Basin, Southwest China
Dai et al. Study of an Automatic Picking Method for Multimode Dispersion Curves of Surface Waves Based on an Improved U-Net
Sun et al. Amplitude-preserving imaging condition for scattering-based RTM in acoustic VTI media

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: 20200421