CN114415234B - 基于主动源面波频散和h/v确定浅地表横波速度的方法 - Google Patents
基于主动源面波频散和h/v确定浅地表横波速度的方法 Download PDFInfo
- Publication number
- CN114415234B CN114415234B CN202210077648.3A CN202210077648A CN114415234B CN 114415234 B CN114415234 B CN 114415234B CN 202210077648 A CN202210077648 A CN 202210077648A CN 114415234 B CN114415234 B CN 114415234B
- Authority
- CN
- China
- Prior art keywords
- wave
- rayleigh
- surface wave
- velocity
- phase
- 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
Links
- 239000006185 dispersion Substances 0.000 title claims abstract description 48
- 238000000034 method Methods 0.000 title claims abstract description 45
- 238000005516 engineering process Methods 0.000 claims abstract description 35
- 238000004458 analytical method Methods 0.000 claims abstract description 16
- 238000005259 measurement Methods 0.000 claims abstract description 14
- 238000003325 tomography Methods 0.000 claims abstract description 13
- 238000004422 calculation algorithm Methods 0.000 claims abstract description 9
- 229910052704 radon Inorganic materials 0.000 claims abstract description 9
- SYUHGPGVQRZVTB-UHFFFAOYSA-N radon atom Chemical compound [Rn] SYUHGPGVQRZVTB-UHFFFAOYSA-N 0.000 claims abstract description 9
- 230000009466 transformation Effects 0.000 claims abstract description 7
- 238000010586 diagram Methods 0.000 claims description 10
- 239000011159 matrix material Substances 0.000 claims description 9
- 230000008569 process Effects 0.000 claims description 9
- 239000002344 surface layer Substances 0.000 claims description 9
- 238000004364 calculation method Methods 0.000 claims description 8
- 238000001914 filtration Methods 0.000 claims description 4
- 238000011835 investigation Methods 0.000 claims description 3
- 238000012360 testing method Methods 0.000 abstract description 10
- 239000010410 layer Substances 0.000 description 15
- 230000035945 sensitivity Effects 0.000 description 9
- 238000000926 separation method Methods 0.000 description 5
- 238000001228 spectrum Methods 0.000 description 4
- 230000003595 spectral effect Effects 0.000 description 3
- 230000008859 change Effects 0.000 description 2
- 238000011161 development Methods 0.000 description 2
- 230000004044 response Effects 0.000 description 2
- 230000003321 amplification Effects 0.000 description 1
- 230000009286 beneficial effect Effects 0.000 description 1
- 230000015572 biosynthetic process Effects 0.000 description 1
- 238000005056 compaction Methods 0.000 description 1
- 239000002131 composite material Substances 0.000 description 1
- 238000005314 correlation function Methods 0.000 description 1
- 230000008878 coupling Effects 0.000 description 1
- 238000010168 coupling process Methods 0.000 description 1
- 238000005859 coupling reaction Methods 0.000 description 1
- 238000013016 damping Methods 0.000 description 1
- 238000009795 derivation Methods 0.000 description 1
- 238000001514 detection method Methods 0.000 description 1
- 230000000694 effects Effects 0.000 description 1
- 230000007613 environmental effect Effects 0.000 description 1
- 230000036039 immunity Effects 0.000 description 1
- 238000000691 measurement method Methods 0.000 description 1
- 238000003199 nucleic acid amplification method Methods 0.000 description 1
- 230000035515 penetration Effects 0.000 description 1
- 230000001902 propagating effect Effects 0.000 description 1
- 230000009467 reduction Effects 0.000 description 1
- 239000011435 rock Substances 0.000 description 1
- 238000004088 simulation Methods 0.000 description 1
- 239000002689 soil Substances 0.000 description 1
- 238000010183 spectrum analysis Methods 0.000 description 1
- 238000006467 substitution reaction Methods 0.000 description 1
- 230000001131 transforming effect Effects 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
-
- 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/301—Analysis for determining seismic cross-sections or geostructures
-
- 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
-
- 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
- Y02A—TECHNOLOGIES FOR ADAPTATION TO CLIMATE CHANGE
- Y02A90/00—Technologies having an indirect contribution to adaptation to climate change
- Y02A90/30—Assessment of water resources
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)
- Measurement Of Mechanical Vibrations Or Ultrasonic Waves (AREA)
Abstract
本发明涉及基于主动源面波频散和H/V确定浅地表横波速度的方法，包括步骤：利用高分辨率线性拉东变换技术从主动源地震记录中分离得到基阶模式瑞雷面波；应用反褶积和频率‑时间分析技术测量每一道记录的瑞雷面波的H/V和瑞雷面波在任意两地震道之间传播的相旅行时；基于任意两道之间的相旅行时测量结果，利用基于直射线理论的层析成像技术计算位于测线上每一个离散网格的瑞雷面波相速度；利用马尔科夫链‑蒙特卡洛算法联合反演瑞雷面波相速度和H/V，获得每一个网格下方的横波速度结构，并构建拟二维横波速度剖面。本发明基于理论和实例测试验证了从主动源地震记录提取并联合反演面波相速度和H/V的技术的可靠性，证明本方案能更准确确定横波速度结构。
Description
技术领域
本发明涉及浅地表地震面波勘探技术领域，特别涉及一种基于主动源面波频散和H/V确定浅地表横波速度的方法。
背景技术
浅地表横波速度是环境和工程勘探中的重要参数，其与土壤和岩石的硬度直接相关。除此之外，浅地表横波速度结构也是强地面运动模拟的重要参数。面波是地震波场的主要成分，而横波速度是决定面波频散特性的最重要的参数。如今，在浅地表地球物理勘探领域，面波分析已经成为确定横波速度结构的最主要方法之一。
在过去几十年，主动源面波勘探技术的发展受到了广泛的关注，浅地表面波分析最早始于面波频谱分析技术(SASW)，其中信号由一个脉冲源激发并由两道构成的一对台站对接收。但是仅用一对台站对测量的面波相速度频散不可避免的会受到噪声和其它类型地震波信号的影响。为了解决此问题，使用面波多道分析技术(MASW)通过同时分析多道地震面波记录以计算频散能量图，并从中提取不同模式的面波相速度频散。基于MASW方法拾取的基阶模式瑞雷面波相速度频散进行反演计算，得到了地震排列中点下方的一维横波速度结构。更近一步，通过增加更多信息（譬如勒夫面波相速度和高阶模式瑞雷面波相速度）进行联合反演，可以进一步的提高横波速度反演结果的精度。
为了准确拾取多模式瑞雷面波信号，采用高分辨率线性拉东变换(HRLRT)技术以提高面波频散能量图的分辨率，并在频率-速度域上对含多模式的瑞雷面波进行模式分离。基于HRLRT技术得到的模式分离后的瑞雷面波，使用基于射线理论的面波层析成像技术计算了位于测线上的每一个离散网格的相速度频散并最终通过反演得到测线下放的拟二维横波速度剖面。相比较于传统的MASW方法，基于层析成像的方法解析了面波相速度的横向变化，继而获得了更高水平分辨率的横波速度剖面。
除了瑞雷面波相速度频散，瑞雷面波的水平垂直振幅比（本文简称为H/V）的相关研究也受到了广泛关注，并被应用于场地响应特征分析，强地面运动预测和浅部横波速度结构的计算等领域。早期通过计算水平和垂直分量原始噪声记录的傅里叶振幅谱的比值测量H/V，但后续的研究表明这种方式测量的H/V不可避免的会受到不同类型的地震波的干扰。为了准确测量高频瑞雷面波H/V，可基于主动源地震记录，利用HRLRT技术进行瑞雷面波的模式分离，并计算分离后的水平分量和垂直分量的基阶模式面波的傅里叶振幅谱的比值以确定H/V。但是，不含任何噪声的准确的模式分离是非常困难的，基于模式分离之后的面波测量的H/V仍然会受到噪声的影响。
虽然主动源面波勘探已经经历了长足的发展，但是仍然存在一些有待解决的问题和制约。因为高频面波的衰减系数很大并且主动源激发的地震波能量有限，通常只能获得几十赫兹以内的瑞雷面波相速度频散，导致瑞雷面波相速度反演难以准确约束最浅部的横波速度结果（譬如，深度小于0.5m）。此外，如果地表存在显著的速度异常（譬如，坚硬的马路），不准确的地表横波速度亦会导致更深部的横波速度反演结果的偏差。目前，很多关于城市或者公路环境下的主动源地震勘探都是直接将检波器与坚硬的马路耦合。譬如，在美国盐湖城市区的公路环境下开展了若干主动源地震勘探工作，对MASW方法得到的横波速度剖面和地震锥渗透测试集(SCPT)的结果进行了比较。结果表明，用拾取的4-35Hz的瑞雷面波相速度频散进行反演，假如不考虑马路的高速度异常而将地表设置为一个低速层(<150m/s)，会导致反演得到的深度2m以浅的平均横波速度被高估超过20%，而反演得到的深度超过5m的平均横波速度被高估10%左右。另外，不准确的横波速度模型亦会影响地震场地响应的估计。软弱的表层会导致地震波能量的“放大”，而坚硬的地表（譬如公路）会导致地震波能量的“减小”。
综上所述，在复杂地表环境下开展面波勘探，急需一种能够准确确定横波速度结构尤其是地表横波速度的可靠方法。
发明内容
本发明的目的在于准确确定地表横波速度结构，提供一种基于主动源面波频散和H/V确定浅地表横波速度的方法。
为了实现上述发明目的，本发明实施例提供了以下技术方案：
基于主动源面波频散和H/V确定浅地表横波速度的方法，包括以下步骤：
步骤S1：利用高分辨率线性拉东变换技术从主动源地震记录中分离得到基阶模式瑞雷面波；
步骤S2：对分离后的基阶模式瑞雷面波应用反褶积和频率-时间分析技术测量每一道记录的瑞雷面波的H/V和瑞雷面波在任意两地震道之间传播的相旅行时；
步骤S3：基于任意两道之间的相旅行时测量结果，利用基于直射线理论的层析成像技术计算位于测线上每一个离散网格的瑞雷面波相速度；
步骤S4：利用马尔科夫链-蒙特卡洛算法联合反演瑞雷面波相速度和H/V，获得每一个网格下方的横波速度结构，并构建拟二维横波速度剖面。
更进一步地，所述步骤S1具体包括以下步骤：
利用高分辨率线性拉东变换技术分别将垂直分量和径向分量的原始地震记录从时间-距离域转换到频率-慢度域，再利用插值的方式将原始地震记录从频率-慢度域转换到频率-速度域，从而获得瑞雷面波频散能量图；
从瑞雷面波的频率-速度域中拾取基阶模式瑞雷面波并将其反变换到时间-距离域，从而获得分离后的基阶模式瑞雷面波信号。
更进一步地，所述步骤S2具体包括以下步骤：
基于高分辨率线性拉东变换技术得到的分离后的基阶模式瑞雷面波记录，将其中每一道的垂直分量的记录分别与其他所有道的垂直分量和径向分量的记录进行反褶积计算，从而获得所有任意两道间的反褶积函数；所述反褶积函数被认为一个地震道作为虚源，而另一个地震道作为接收点观测的信号；
利用频率-时间分析技术测量两道之间反褶积函数中瑞雷面波的振幅以及相旅行时、群旅行时；
对于任意两道，计算第一道的径向分量-第二道的垂直分量和第一道的垂直分量-第二道的垂直分量的反褶积函数中瑞雷面波信号的振幅比，作为第一道的瑞雷面波H/V；或计算第一道的垂直分量-第二道的径向分量和第一道的垂直分量-第二道的垂直分量的反褶积函数中瑞雷面波信号的振幅比，作为第二道的瑞雷面波H/V；对于任意两道，舍弃低信噪比的反褶积函数或不同分量的瑞雷面波相、群旅行时不一致的反褶积函数。
更进一步地，所述利用频率-时间分析技术测量两道之间反褶积函数中瑞雷面波的振幅以及相旅行时、群旅行时的步骤，包括：
利用频率-时间分析技术基于不同的中心频率进行高斯窄带滤波，对每一个高斯窄带滤波后的波形做希尔伯特变换得到信号的包络；通过寻找包络的最大值从而确定瑞雷面波的振幅，通过测量振幅最大值所在时刻确定瑞雷面波的群旅行时；并且通过测量振幅最大值所在时刻的瞬时相位和瞬时频率确定瑞雷面波的相旅行时。
更进一步地，所述步骤S3具体包括以下步骤：
基于由垂直分量-垂直分量的反褶积函数中瑞雷面波信号测量的任意两道间的相旅行时，采用直射线理论的面波层析成像技术计算不同频率的瑞雷面波相速度的横向变化；
在正演时，将一维测线剖分为若干离散的网格；对于一个特定频率f，给定路径的瑞雷面波的相旅行时表达为该射线穿过每一个网格所需的相旅行时的叠加；穿过研究区域的所有射线的面波相旅行时由以下公式表示：
d=Gm
其中G表示关于射线路径的矩阵，d表示相旅行时的矩阵，m表示与频率f相关的慢度矩阵；通过算子G将模型空间m投影到数据空间d；
以此为正演过程，通过反演计算得到不同频率的相速度分布，其中反演过程中使用的目标函数为。
更进一步地，所述步骤S4具体包括以下步骤：
将每一道的H/V测量结果进行插值得到不同网格的H/V，网格的大小与层析成像中所使用的网格大小相同；使用马尔科夫链-蒙特卡洛算法联合反演瑞雷面波的相速度频散和H/V，获取每一个网格下方的一维横波速度结构；
其中，利用瑞雷面波相速度频散确定联合反演所使用的横波速度的初始模型；且要求以初始模型为基准，在允许范围内进行扰动；
具体的扰动范围设置为：表层横波速度的扰动范围从40m/s到初始模型横波速度的1000%，表层以下横波速度的扰动范围是初始模型的横波速度的±50%。
更进一步地，利用瑞雷面波相速度频散确定横波速度的初始模型，以初始模型为基准，在允许范围内进行扰动获得最终的反演模型的步骤，包括：不同频率的瑞雷面波相速度对横波速度的敏感深度不同，根据经验公式利用不同频率的瑞雷面波相速度确定不同深度的横波速度的初始模型。以初始模型为起点，进行满足高斯概率分布的随机扰动得到新模型，接着利用Metropolis准则对新模型的质量进行评价。如果新模型被接受，则以新模型作为基准进行后续的扰动；如果新模型被舍弃，则仍然以之前的模型为基准进行扰动。最后，将所有被接受的模型的平均值作为最终的反演模型，而所有被接受的模型的标准差作为反演结果的不确定度。
与现有技术相比，本发明的有益效果：
本发明通过对比分析瑞雷面波相速度和H/V对于不同深度的横波速度的敏感度，发现瑞雷面波H/V相较于相同频率的瑞雷面波相速度对于横波速度具有明显更浅的敏感深度。因此本发明提出了一种基于理论和实例测试验证了从主动源地震记录提取并联合反演面波相速度和H/V的技术的可靠性，证明本方案能更准确确定横波速度结构；并且最后讨论了所提出的技术对于在复杂地表环境下开展地震探测的潜力。
附图说明
为了更清楚地说明本发明实施例的技术方案，下面将对实施例中所需要使用的附图作简单地介绍， 应当理解，以下附图仅示出了本发明的某些实施例，因此不应被看作是对范围的限定，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他相关的附图。
图1为本发明流程图；
图2为本发明实施例2基于两层理论模型计算的相同频率的瑞雷面波H/V和相速度对不同深度的横波速度敏感度的对比图，其中图2（a）为基于两层理论模型分别计算的频率为20Hz、40Hz、60Hz的瑞雷面波相速度对不同深度的横波速度敏感度，图2（b）为基于两层理论模型分别计算的频率为20Hz、40Hz、60Hz的瑞雷面波H/V对不同深度的横波速度敏感度；
图3为本发明实施例2比较了通过两道间反褶积函数和FTAN技术测量的H/V与通过单道的傅里叶振幅谱比测量的H/V，其中图3（a）展示了通过傅里叶谱比测量的H/V与理论H/V，图3（b）展示了两道间的反褶积函数和FTAN技术测量的H/V与理论H/V；
图4展示了基于地表为高速层的阶梯状模型的合成记录的测试结果，其中图4（a）展示了观测的相速度频散与基于相速度单独反演结果正演得到的瑞雷面波相速度频散，图4（b）展示了观测的瑞雷面波H/V与基于相速度单独反演结果正演得到的瑞雷面波H/V，图4（c）为相速度单独反演得到的深度1m以浅的结果，图4（d）为相速度单独反演得到的深度10m以浅的结果；
图5为基于瑞雷面波相速度和H/V联合反演的测试结果示意图，其中图5（a）展示了观测的相速度频散与基于联合反演结果正演得到的瑞雷面波相速度频散，图5（b）展示了观测的瑞雷面波相速度H/V与基于联合反演结果正演得到的瑞雷面波H/V，图5（c）为瑞雷面波相速度和H/V联合反演得到的深度1m以浅的结果，图5（d）为瑞雷面波相速度和H/V联合反演得到的深度10m以浅的结果。
具体实施方式
下面将结合本发明实施例中附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。通常在此处附图中描述和示出的本发明实施例的组件可以以各种不同的配置来布置和设计。因此，以下对在附图中提供的本发明的实施例的详细描述并非旨在限制要求保护的本发明的范围，而是仅仅表示本发明的选定实施例。基于本发明的实施例，本领域技术人员在没有做出创造性劳动的前提下所获得的所有其他实施例，都属于本发明保护的范围。
应注意到：相似的标号和字母在下面的附图中表示类似项，因此，一旦某一项在一个附图中被定义，则在随后的附图中不需要对其进行进一步定义和解释。同时，在本发明的描述中，术语“第一”、“第二”等仅用于区分描述，而不能理解为指示或暗示相对重要性，或者暗示这些实体或操作之间存在任何这种实际的关系或者顺序。
实施例1：
本发明通过下述技术方案实现，基于主动源面波频散和h/v确定浅地表横波速度的方法，包括以下步骤：
为了不让权利要求书和说明书一样，可以在以下每个步骤中给出技术实施时的实验例子，再配合实验数据等，就会形成一个对应于技术步骤的完整实施例。
步骤S1：利用高分辨率线性拉东变换(HRLRT)技术从主动源地震记录中分离得到基阶模式瑞雷面波。
为了获得瑞雷面波频散能量图，先利用高分辨率线性拉东变换(HRLRT)技术分别将垂直分量和径向分量的原始地震记录从时间-距离域转换到频率-慢度域，再利用插值的方式将面波数据从频率-慢度域转换到频率-速度域。然后从包含瑞雷面波的频散能量图，也就是频率-速度域中拾取基阶模式瑞雷面波，并将其反变换到时间-距离域，即获得了分离后的基阶模式瑞雷面波信号。
步骤S2：对分离后的基阶模式瑞雷面波应用反褶积和频率-时间分析技术测量每一道记录的瑞雷面波的H/V和瑞雷面波在任意两地震道之间传播的相旅行时。
基于高分辨率线性拉东变换(HRLRT)技术得到的分离后的基阶模式瑞雷面波记录，将其中每一道的垂直分量的记录分别与其他所有道的垂直分量和径向分量的记录进行反褶积计算，获得所有任意两道间的反褶积函数。所述反褶积函数可以被认为是一个地震道作为虚源，而另一个地震道作为接收点观测的信号。
利用频率-时间分析(FTAN)技术测量两道之间反褶积函数中瑞雷面波的振幅以及相、群旅行时。频率-时间分析(FTAN)技术基于不同的中心频率进行一系列的高斯窄带滤波，接着对每一个高斯窄带滤波后的波形做希尔伯特变换得到信号的包络。通过寻找包络的最大值就可以确定瑞雷面波的振幅，通过测量振幅最大值所在时刻即为瑞雷面波的群旅行时。再通过测量振幅最大值所在时刻的瞬时相位和瞬时频率，即可确定瑞雷面波的相旅行时。
对于任意两道，第一道的径向(R)分量-第二道的垂直(Z)分量和第一道的垂直(z)分量-第二道的垂直(Z)分量的反褶积函数中瑞雷面波信号的振幅比，作为第一道的瑞雷面波H/V；或计算第一道的垂直(Z)分量-第二道的径向(R)分量和第一道的垂直(Z)分量-第二道的垂直(Z)分量的反褶积函数中瑞雷面波信号的振幅比，作为第二道的瑞雷面波H/V。对于任意两道，舍弃低信噪比的反褶积函数或不同分量的瑞雷面波相、群旅行时不一致的反褶积函数。
最后，对于每一道，将剩余的H/V测量结果的平均值作为最终的H/V测量结果，标准差作为H/V测量结果的不确定度。然而，瑞雷面波的H/V会受到复杂结构或者方位各向异性的影响。而本方案基于一维测线的观测系统难以充分考虑这些因素的影响，会导致测量的H/V的不确定度偏小。考虑到具有良好方位角覆盖的天然地震或者基于二维台阵观测的背景噪声互相关函数测量的H/V的不确定度通常小于0.02，因此在本方案中如果H/V测量结果的不确定度小于0.025，则统一将其设为0.025。
步骤S3：基于任意两道之间的相旅行时测量结果，利用基于直射线理论的层析成像技术计算位于测线上每一个离散网格的瑞雷面波相速度。
基于由垂直(Z)分量-垂直(Z)分量的反褶积函数测量的瑞雷面波信号在任意两道间的相旅行时，采用直射线理论的面波层析成像技术计算不同频率的瑞雷面波相速度的横向变化。在正演时，将一维测线剖分为若干离散的网格；对于一个特定频率f，给定路径的瑞雷面波的相旅行时表达为该射线穿过每一个网格所需的相旅行时的叠加。
因此穿过研究区域的所有射线的面波相旅行时由以下公式表示：d=Gm，其中G表示关于射线路径的矩阵，d表示相旅行时的矩阵，m表示与频率f相关的慢度矩阵；通过算子G将模型空间m投影到数据空间d。
以此为正演过程，可以通过反演计算得到不同频率的相速度分布。将反演过程中使用的目标函数定义为，通过一系列的检测板测试确定离散网格合适的大小，并利用L-曲线寻找最优的阻尼因子。
从瑞雷面波频散能量图中拾取的面波相速度频散的不确定度通常被定义为频散能量的中值的半宽度。因为面波层析成像技术具有较强的抗噪声能力，并且相对于瑞雷面波频散能量图能更好的分辨相速度水平方向的变化，所以本方案将相速度的不确定度同一设置为由频散能量图拾取的相速度的不确定的0.5倍。
步骤S4：利用马尔科夫链-蒙特卡洛算法(MCMC)联合反演瑞雷面波相速度和H/V，获得每一个网格下方的横波速度结构，并构建拟二维横波速度剖面。
马尔科夫链-蒙特卡洛算法(MCMC)在足够的随机搜索次数的前提下，可以搜索到全局极小值并获得反演结果的不确定度。
为了进行联合反演，首先将每一道的H/V测量结果进行插值得到不同网格的H/V，网格的大小与层析成像中所使用的网格大小相同。使用马尔科夫链-蒙特卡洛算法(MCMC)联合反演瑞雷面波的相速度频散和H/V，获取每一个网格下方的一维横波速度结构。
在本步骤中，通过瑞雷面波相速度频散确定横波速度的初始模型，以初始模型为基准，在允许范围内进行扰动得到最终的横波速度的反演模型。模型的扰动范围设置如下：表层横波速度的扰动范围从40m/s到初始模型横波速度的1000%；表层以下横波速度的扰动范围是初始模型横波速度的±50%。
具体的反演过程为：根据经验公式利用不同频率的瑞雷面波相速度计算不同深度的横波速度的初始模型。以初始模型为起点，进行满足高斯概率分布的随机扰动得到新模型，接着利用Metropolis准则对新模型的质量进行评价。如果新模型被接受，则以新模型作为基准进行后续的扰动；如果新模型被舍弃，则仍然以之前的模型为基准进行扰动。最后，将所有被接受的模型的平均值作为最终的反演模型，而所有被接受的模型的标准差作为反演结果的不确定度。
一般而言，由于地层的压实效应，横波速度随着深度的增加而增加。因此本方案增加了在深度1.5m以下横波速度随深度增加而增加的限制，以减小反演结果的不确定度。需要注意的是增加了这个约束如果导致基于反演结果的预测值不能很好的拟合观测数据，则应该取消该约束。因为横波速度是瑞雷面波频散和H/V最敏感的参数，为了使反演过程更稳定，在反演中假设不同深度的纵波速度和密度是已知的，而只反演得到横波速度。
实施例2：
基于实施例1的技术方案，本实施例进行相应的理论测试。
基于有限差分算子计算二维P-SV波动方程得到所有垂直分量和径向分量的理论合成数据。其中，震源子波函数被定义为雷克子波，主频为15Hz。为了保证有限差分求导的稳定性，理论模型被均匀剖分成0.05×0.05m的网格，时间域上的步长为0.1ms。
本实施例使用一个两层理论模型测试瑞雷面波H/V测量方法的鲁棒性。为了模拟野外记录的地震波场，在这个两层理论模型的合成数据中添加了高斯随机噪声。其中，高斯随机噪声的平均值是0，方差定义为第一道垂直分量理论合成数据的最大振幅的2%。基于加噪后的理论数据，利用HRLRT方法计算频散能量图，并从中分离出基阶模式瑞雷面波。基于分离得到的基阶模式瑞雷面波，使用任意两道间反褶积函数和FTAN技术测量每一道的瑞雷面波H/V。
请参见图2，其中图2（a）为基于两层理论模型分别计算的频率为20Hz、40Hz、60Hz的瑞雷面波相速度对不同深度的横波速度的敏感度，图2（b）为基于两层理论模型分别计算的频率为20Hz、40Hz、60Hz的瑞雷面波H/V对不同深度的横波速度的敏感度。
作为一个例子，图3基于上述加噪记录测量的基阶模式瑞雷面波H/V，比较了通过两道间反褶积函数和FTAN技术测量的H/V与通过单道的傅里叶振幅谱比测量的H/V，其中展示了第20道H/V的测量。请参见图3（a），通过傅里叶谱比测量的H/V与理论H/V在不同频段尤其是>40Hz的频段存在明显的偏差，其中虚线代表利用傅里叶谱比测量的第20道的H/V，实线代表基于两层理论模型计算的H/V。与之相反，请参见图3（b），使用两道间反褶积函数和FTAN技术测量的H/V与理论H/V在不同频段都非常接近，所有偏差都在测量结果的误差棒的范围内，其中圆点代表利用反褶积函数和FTAN技术测量的第20道的H/V，误差棒代表对应的不确定度，实线代表基于两层理论模型计算的H/V。
更进一步地，本实施例使用地表为高速层的阶梯状模型模拟城市地表环境下的复杂结构。根据提出的技术方案，计算瑞雷面波的相速度频散和H/V，然后利用MCMC算法分别进行瑞雷面波相速度频散的单独反演，以及瑞雷面波相速度频散和H/V的联合反演，获得地下横波速度结构，请参见图4基于地表为高速层的阶梯状模型的合成记录的测试结果。
结果表明，请参见图4（a）、4（b），误差棒代表观测的某一网格的瑞雷面波相速度频散和H/V，实线代表基于MCMC反演得到的所有接受模型的正演结果。基于单独反演结果正演预测的相速度可以很好的拟合观测相速度，但是正演预测的H/V则不能拟合观测H/V，请参见图4（c）、4（d），表明单独反演得到的横波速度模型在表层和半空间与理论模型存在明显的偏差，其中虚线代表由瑞雷面波相速度反演得到的深度1m以浅的所有接受模型的平均结果，实线代表理论模型，表示了深度10m以浅的横波速度反演结果。譬如，反演得到的表层横波速度（430m/s）明显小于理论模型的横波速度（1300m/s），且具有较大的不确定度。而反演得到的第二层（0.2~0.5m）、第三层（0.5~1.5m）和半空间（>4.5m）的横波速度则稍大于理论模型的横波速度。
总体而言，基于瑞雷面波相速度频散单独反演得到的横波速度模型与理论模型差异较大。与之相反，本方案通过引入H/V进行联合反演，基于反演结果的预测相速度和H/V都能很好的拟合观测值，请参见图5（a）、图5（b）。并且反演得到的横波速度模型与理论横波速度模型在不同深度都非常接近，请参见图5（c）、图5（d），参数表示同图4。
需要注意的是，联合反演得到的地表横波速度和理论横波速度非常接近，且反演结果的不确定度非常小（<100m/s）。理论测试表明通过引入H/V进行联合反演可以显著改善横波速度反演结果的准确度。
以上所述，仅为本发明的具体实施方式，但本发明的保护范围并不局限于此，任何熟悉本技术领域的技术人员在本发明揭露的技术范围内，可轻易想到变化或替换，都应涵盖在本发明的保护范围之内。因此，本发明的保护范围应所述以权利要求的保护范围为准。
Claims (3)
1.基于主动源面波频散和H/V确定浅地表横波速度的方法，其特征在于：包括以下步骤：
步骤S1：利用高分辨率线性拉东变换技术从主动源地震记录中分离得到基阶模式瑞雷面波；
步骤S2：对分离后的基阶模式瑞雷面波应用反褶积和频率-时间分析技术测量每一道记录的瑞雷面波的H/V和瑞雷面波在任意两地震道之间传播的相旅行时；
所述步骤S2具体包括以下步骤：
基于高分辨率线性拉东变换技术得到的分离后的基阶模式瑞雷面波记录，将其中每一道的垂直分量的记录分别与其他所有道的垂直分量和径向分量的记录进行反褶积计算，从而获得所有任意两道间的反褶积函数；所述反褶积函数被认为一个地震道作为虚源，而另一个地震道作为接收点观测的信号；
利用频率-时间分析技术测量两道之间反褶积函数中瑞雷面波的振幅以及相旅行时、群旅行时；
对于任意两道，计算第一道的径向分量-第二道的垂直分量和第一道的垂直分量-第二道的垂直分量的反褶积函数中瑞雷面波信号的振幅比，作为第一道的瑞雷面波H/V；或计算第一道的垂直分量-第二道的径向分量和第一道的垂直分量-第二道的垂直分量的反褶积函数中瑞雷面波信号的振幅比，作为第二道的瑞雷面波H/V；对于任意两道，舍弃低信噪比的反褶积函数或不同分量的瑞雷面波相、群旅行时不一致的反褶积函数；
所述利用频率-时间分析技术测量两道之间反褶积函数中瑞雷面波的振幅以及相旅行时、群旅行时的步骤，包括：
利用频率-时间分析技术基于不同的中心频率进行高斯窄带滤波，对每一个高斯窄带滤波后的波形做希尔伯特变换得到信号的包络；通过寻找包络的最大值从而确定瑞雷面波的振幅，通过测量振幅最大值所在时刻确定瑞雷面波的群旅行时；并且通过测量振幅最大值所在时刻的瞬时相位和瞬时频率确定瑞雷面波的相旅行时；
步骤S3：基于任意两道之间的相旅行时测量结果，利用基于直射线理论的层析成像技术计算位于测线上每一个离散网格的瑞雷面波相速度；
所述步骤S3具体包括以下步骤：
基于由垂直分量-垂直分量的反褶积函数中瑞雷面波信号测量的任意两道间的相旅行时，采用直射线理论的面波层析成像技术计算不同频率的瑞雷面波相速度的横向变化；
在正演时，将一维测线剖分为若干离散的网格；对于一个特定频率f，给定路径的瑞雷面波的相旅行时表达为该射线穿过每一个网格所需的相旅行时的叠加；穿过研究区域的所有射线的面波相旅行时由以下公式表示：
d=Gm
其中G表示关于射线路径的矩阵，d表示相旅行时的矩阵，m表示与频率f相关的慢度矩阵；通过算子G将模型空间m投影到数据空间d；
以此为正演过程，通过反演计算得到不同频率的相速度分布，其中反演过程中使用的目标函数为，为目标函数中的正则化参数；
步骤S4：利用马尔科夫链-蒙特卡洛算法联合反演瑞雷面波相速度和H/V，获得每一个网格下方的横波速度结构，并构建拟二维横波速度剖面。
2.根据权利要求1所述的基于主动源面波频散和H/V确定浅地表横波速度的方法，其特征在于：所述步骤S1具体包括以下步骤：
利用高分辨率线性拉东变换技术分别将垂直分量和径向分量的原始地震记录从时间-距离域转换到频率-慢度域，再利用插值的方式将原始地震记录从频率-慢度域转换到频率-速度域，从而获得瑞雷面波频散能量图；
从瑞雷面波的频率-速度域中拾取基阶模式瑞雷面波并将其反变换到时间-距离域，从而获得分离后的基阶模式瑞雷面波信号。
3.根据权利要求2所述的基于主动源面波频散和H/V确定浅地表横波速度的方法，其特征在于：所述步骤S4具体包括以下步骤：
将每一道的H/V测量结果进行插值得到不同网格的H/V，网格的大小与层析成像中所使用的网格大小相同；使用马尔科夫链-蒙特卡洛算法联合反演瑞雷面波的相速度频散和H/V，获取每一个网格下方的一维横波速度结构；
其中，利用瑞雷面波相速度频散确定联合反演所使用的横波速度的初始模型；且反演过程中要求以初始模型为基准，在允许范围内进行扰动；
具体的扰动范围为：表层横波速度的扰动范围从40m/s到初始模型横波速度的1000%，表层以下横波速度的扰动范围是初始模型横波速度的±50%。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202210077648.3A CN114415234B (zh) | 2022-01-24 | 2022-01-24 | 基于主动源面波频散和h/v确定浅地表横波速度的方法 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202210077648.3A CN114415234B (zh) | 2022-01-24 | 2022-01-24 | 基于主动源面波频散和h/v确定浅地表横波速度的方法 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN114415234A CN114415234A (zh) | 2022-04-29 |
| CN114415234B true CN114415234B (zh) | 2023-05-12 |
Family
ID=81275075
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202210077648.3A Active CN114415234B (zh) | 2022-01-24 | 2022-01-24 | 基于主动源面波频散和h/v确定浅地表横波速度的方法 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN114415234B (zh) |
Families Citing this family (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN117310810A (zh) * | 2023-09-26 | 2023-12-29 | 西南交通大学 | 联合反演面波频散和h/v获取横波速度及径向各向异性 |
Citations (8)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN101419292A (zh) * | 2007-10-25 | 2009-04-29 | 中国石油天然气集团公司 | 采用纵波源多分量地震数据生成横波地震剖面的方法 |
| CN102419454A (zh) * | 2011-06-30 | 2012-04-18 | 中国科学院地质与地球物理研究所 | 隧道掌子面前方远距离含水目标体的瞬变电磁预报方法 |
| CN102749646A (zh) * | 2012-07-06 | 2012-10-24 | 西安石油大学 | 一种瑞雷面波深度-频率分析方法 |
| CN102759749A (zh) * | 2012-07-06 | 2012-10-31 | 西安石油大学 | 一种瑞雷面波速度分析方法 |
| CN103487825A (zh) * | 2013-09-30 | 2014-01-01 | 雷文太 | 一种运营高速公路路基病害瑞雷波自动检测装置 |
| CN104849752A (zh) * | 2015-05-18 | 2015-08-19 | 中国石油集团川庆钻探工程有限公司地球物理勘探公司 | 获取地震记录中的转换横波的高精度反褶积算子的方法 |
| WO2017165341A2 (en) * | 2016-03-25 | 2017-09-28 | Schlumberger Technology Corporation | Method and device for estimating sonic slowness in a subterranean formation |
| CN107561589A (zh) * | 2017-10-25 | 2018-01-09 | 中国石油化工股份有限公司 | 一种近地表横波层速度模型建立方法 |
Family Cites Families (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US11099290B2 (en) * | 2019-06-12 | 2021-08-24 | Emerson Paradigm Holding Llc | Attenuating surface waves in common shot gathers of seismic data collected by a set of geophones |
-
2022
- 2022-01-24 CN CN202210077648.3A patent/CN114415234B/zh active Active
Patent Citations (8)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN101419292A (zh) * | 2007-10-25 | 2009-04-29 | 中国石油天然气集团公司 | 采用纵波源多分量地震数据生成横波地震剖面的方法 |
| CN102419454A (zh) * | 2011-06-30 | 2012-04-18 | 中国科学院地质与地球物理研究所 | 隧道掌子面前方远距离含水目标体的瞬变电磁预报方法 |
| CN102749646A (zh) * | 2012-07-06 | 2012-10-24 | 西安石油大学 | 一种瑞雷面波深度-频率分析方法 |
| CN102759749A (zh) * | 2012-07-06 | 2012-10-31 | 西安石油大学 | 一种瑞雷面波速度分析方法 |
| CN103487825A (zh) * | 2013-09-30 | 2014-01-01 | 雷文太 | 一种运营高速公路路基病害瑞雷波自动检测装置 |
| CN104849752A (zh) * | 2015-05-18 | 2015-08-19 | 中国石油集团川庆钻探工程有限公司地球物理勘探公司 | 获取地震记录中的转换横波的高精度反褶积算子的方法 |
| WO2017165341A2 (en) * | 2016-03-25 | 2017-09-28 | Schlumberger Technology Corporation | Method and device for estimating sonic slowness in a subterranean formation |
| CN107561589A (zh) * | 2017-10-25 | 2018-01-09 | 中国石油化工股份有限公司 | 一种近地表横波层速度模型建立方法 |
Non-Patent Citations (2)
| Title |
|---|
| 主动源和被动源面波浅勘方法综述;刘庆华;鲁来玉;王凯明;;地球物理学进展(06);第2906-2922页 * |
| 消除近地表地震散射噪音的方法;杨旭明,周熙襄,王克斌,石生林,屠世杰;成都理工学院学报(04);第428-432页 * |
Also Published As
| Publication number | Publication date |
|---|---|
| CN114415234A (zh) | 2022-04-29 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| Park et al. | Multimodal analysis of high frequency surface waves | |
| CN109669212B (zh) | 地震数据处理方法、地层品质因子估算方法与装置 | |
| EP0060029B1 (en) | A method of determining the ratio of the velocities of compressional and shear waves in subterranean formations | |
| US20100265797A1 (en) | Interferometric seismic data processing | |
| US20120113751A1 (en) | Statics calculation | |
| CN105388518A (zh) | 一种质心频率与频谱比联合的井中地震品质因子反演方法 | |
| CN109884709B (zh) | 一种基于面波旅行时层析的转换波静校正方法 | |
| US20120053839A1 (en) | Method of detecting or monitoring a subsurface hydrocarbon reservoir-sized structure | |
| CN107884829A (zh) | 一种联合压制浅海obc地震资料多次波的方法 | |
| Hayashi et al. | CMP spatial autocorrelation analysis of multichannel passive surface-wave data | |
| CN104570116A (zh) | 基于地质标志层的时差分析校正方法 | |
| CN101852864B (zh) | 一种利用地表一致性统计频谱分析技术处理海量地震数据的方法 | |
| Muir et al. | Rayleigh‐wave H/V via noise cross correlation in southern California | |
| Zhang et al. | Retrieval of shallow S-wave profiles from seismic reflection surveying and traffic-induced noise | |
| Yin et al. | Improving horizontal resolution of high-frequency surface-wave methods using travel-time tomography | |
| CN114415234B (zh) | 基于主动源面波频散和h/v确定浅地表横波速度的方法 | |
| Meehan et al. | Reconstruction of historical surface mass balance, 1984–2017 from GreenTrACS multi-offset ground-penetrating radar | |
| EA030770B1 (ru) | Система и способ адаптивной сейсмической оптики | |
| CN102323618A (zh) | 基于分数阶傅里叶变换的相干噪声抑制方法 | |
| CN107918152B (zh) | 一种地震相干层析成像方法 | |
| CN110780346A (zh) | 一种隧道超前探测复杂地震波场的分离方法 | |
| CN113552632B (zh) | 基于小波域卡尔曼滤波的地震互相关信号拾取方法和系统 | |
| CN102890288B (zh) | 一种地震波层速度反演方法 | |
| Woo et al. | Processing Ambient Noise Data Using Phase Cross‐Correlation and Application Toward Understanding Spatiotemporal Environmental Effects | |
| CN103984013A (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 | ||
| GR01 | Patent grant | ||
| GR01 | Patent grant |