基于主动源面波频散和H/V确定浅地表横波速度的方法

基于主动源面波频散和h/v确定浅地表横波速度的方法
技术领域
1.本发明涉及浅地表地震面波勘探技术领域，特别涉及一种基于主动源面波频散和h/v确定浅地表横波速度的方法。


背景技术：

2.浅地表横波速度是环境和工程勘探中的重要参数，其与土壤和岩石的硬度直接相关。除此之外，浅地表横波速度结构也是强地面运动模拟的重要参数。面波是地震波场的主要成分，而横波速度是决定面波频散特性的最重要的参数。如今，在浅地表地球物理勘探领域，面波分析已经成为确定横波速度结构的最主要方法之一。
3.在过去几十年，主动源面波勘探技术的发展受到了广泛的关注，浅地表面波分析最早始于面波频谱分析技术(sasw)，其中信号由一个脉冲源激发并由两道构成的一对台站对接收。但是仅用一对台站对测量的面波相速度频散不可避免的会受到噪声和其它类型地震波信号的影响。为了解决此问题，使用面波多道分析技术(masw)通过同时分析多道地震面波记录以计算频散能量图，并从中提取不同模式的面波相速度频散。基于masw方法拾取的基阶模式瑞雷面波相速度频散进行反演计算，得到了地震排列中点下方的一维横波速度结构。更近一步，通过增加更多信息（譬如勒夫面波相速度和高阶模式瑞雷面波相速度）进行联合反演，可以进一步的提高横波速度反演结果的精度。
4.为了准确拾取多模式瑞雷面波信号，采用高分辨率线性拉东变换(hrlrt)技术以提高面波频散能量图的分辨率，并在频率-速度域上对含多模式的瑞雷面波进行模式分离。基于hrlrt技术得到的模式分离后的瑞雷面波，使用基于射线理论的面波层析成像技术计算了位于测线上的每一个离散网格的相速度频散并最终通过反演得到测线下放的拟二维横波速度剖面。相比较于传统的masw方法，基于层析成像的方法解析了面波相速度的横向变化，继而获得了更高水平分辨率的横波速度剖面。
5.除了瑞雷面波相速度频散，瑞雷面波的水平垂直振幅比（本文简称为h/v）的相关研究也受到了广泛关注，并被应用于场地响应特征分析，强地面运动预测和浅部横波速度结构的计算等领域。早期通过计算水平和垂直分量原始噪声记录的傅里叶振幅谱的比值测量h/v，但后续的研究表明这种方式测量的h/v不可避免的会受到不同类型的地震波的干扰。为了准确测量高频瑞雷面波h/v，可基于主动源地震记录，利用hrlrt技术进行瑞雷面波的模式分离，并计算分离后的水平分量和垂直分量的基阶模式面波的傅里叶振幅谱的比值以确定h/v。但是，不含任何噪声的准确的模式分离是非常困难的，基于模式分离之后的面波测量的h/v仍然会受到噪声的影响。
6.虽然主动源面波勘探已经经历了长足的发展，但是仍然存在一些有待解决的问题和制约。因为高频面波的衰减系数很大并且主动源激发的地震波能量有限，通常只能获得几十赫兹以内的瑞雷面波相速度频散，导致瑞雷面波相速度反演难以准确约束最浅部的横波速度结果（譬如，深度小于0.5m）。此外，如果地表存在显著的速度异常（譬如，坚硬的马路），不准确的地表横波速度亦会导致更深部的横波速度反演结果的偏差。目前，很多关于
城市或者公路环境下的主动源地震勘探都是直接将检波器与坚硬的马路耦合。譬如，在美国盐湖城市区的公路环境下开展了若干主动源地震勘探工作，对masw方法得到的横波速度剖面和地震锥渗透测试集(scpt)的结果进行了比较。结果表明，用拾取的4-35hz的瑞雷面波相速度频散进行反演，假如不考虑马路的高速度异常而将地表设置为一个低速层(《150m/s)，会导致反演得到的深度2m以浅的平均横波速度被高估超过20%，而反演得到的深度超过5m的平均横波速度被高估10%左右。另外，不准确的横波速度模型亦会影响地震场地响应的估计。软弱的表层会导致地震波能量的“放大”，而坚硬的地表（譬如公路）会导致地震波能量的“减小”。
7.综上所述，在复杂地表环境下开展面波勘探，急需一种能够准确确定横波速度结构尤其是地表横波速度的可靠方法。


技术实现要素：

8.本发明的目的在于准确确定地表横波速度结构，提供一种基于主动源面波频散和h/v确定浅地表横波速度的方法。
9.为了实现上述发明目的，本发明实施例提供了以下技术方案：基于主动源面波频散和h/v确定浅地表横波速度的方法，包括以下步骤：步骤s1：利用高分辨率线性拉东变换技术从主动源地震记录中分离得到基阶模式瑞雷面波；步骤s2：对分离后的基阶模式瑞雷面波应用反褶积和频率-时间分析技术测量每一道记录的瑞雷面波的h/v和瑞雷面波在任意两地震道之间传播的相旅行时；步骤s3：基于任意两道之间的相旅行时测量结果，利用基于直射线理论的层析成像技术计算位于测线上每一个离散网格的瑞雷面波相速度；步骤s4：利用马尔科夫链-蒙特卡洛算法联合反演瑞雷面波相速度和h/v，获得每一个网格下方的横波速度结构，并构建拟二维横波速度剖面。
10.更进一步地，所述步骤s1具体包括以下步骤：利用高分辨率线性拉东变换技术分别将垂直分量和径向分量的原始地震记录从时间-距离域转换到频率-慢度域，再利用插值的方式将原始地震记录从频率-慢度域转换到频率-速度域，从而获得瑞雷面波频散能量图；从瑞雷面波的频率-速度域中拾取基阶模式瑞雷面波并将其反变换到时间-距离域，从而获得分离后的基阶模式瑞雷面波信号。
11.更进一步地，所述步骤s2具体包括以下步骤：基于高分辨率线性拉东变换技术得到的分离后的基阶模式瑞雷面波记录，将其中每一道的垂直分量的记录分别与其他所有道的垂直分量和径向分量的记录进行反褶积计算，从而获得所有任意两道间的反褶积函数；所述反褶积函数被认为一个地震道作为虚源，而另一个地震道作为接收点观测的信号；利用频率-时间分析技术测量两道之间反褶积函数中瑞雷面波的振幅以及相旅行时、群旅行时；对于任意两道，计算第一道的径向分量-第二道的垂直分量和第一道的垂直分量-第二道的垂直分量的反褶积函数中瑞雷面波信号的振幅比，作为第一道的瑞雷面波h/v；或
计算第一道的垂直分量-第二道的径向分量和第一道的垂直分量-第二道的垂直分量的反褶积函数中瑞雷面波信号的振幅比，作为第二道的瑞雷面波h/v；对于任意两道，舍弃低信噪比的反褶积函数或不同分量的瑞雷面波相、群旅行时不一致的反褶积函数。
12.更进一步地，所述利用频率-时间分析技术测量两道之间反褶积函数中瑞雷面波的振幅以及相旅行时、群旅行时的步骤，包括：利用频率-时间分析技术基于不同的中心频率进行高斯窄带滤波，对每一个高斯窄带滤波后的波形做希尔伯特变换得到信号的包络；通过寻找包络的最大值从而确定瑞雷面波的振幅，通过测量振幅最大值所在时刻确定瑞雷面波的群旅行时；并且通过测量振幅最大值所在时刻的瞬时相位和瞬时频率确定瑞雷面波的相旅行时。
13.更进一步地，所述步骤s3具体包括以下步骤：基于由垂直分量-垂直分量的反褶积函数中瑞雷面波信号测量的任意两道间的相旅行时，采用直射线理论的面波层析成像技术计算不同频率的瑞雷面波相速度的横向变化；在正演时，将一维测线剖分为若干离散的网格；对于一个特定频率f，给定路径的瑞雷面波的相旅行时表达为该射线穿过每一个网格所需的相旅行时的叠加；穿过研究区域的所有射线的面波相旅行时由以下公式表示：d=gm其中g表示关于射线路径的矩阵，d表示相旅行时的矩阵，m表示与频率f相关的慢度矩阵；通过算子g将模型空间m投影到数据空间d；以此为正演过程，通过反演计算得到不同频率的相速度分布，其中反演过程中使用的目标函数为。
14.更进一步地，所述步骤s4具体包括以下步骤：将每一道的h/v测量结果进行插值得到不同网格的h/v，网格的大小与层析成像中所使用的网格大小相同；使用马尔科夫链-蒙特卡洛算法联合反演瑞雷面波的相速度频散和h/v，获取每一个网格下方的一维横波速度结构；其中，利用瑞雷面波相速度频散确定联合反演所使用的横波速度的初始模型；且要求以初始模型为基准，在允许范围内进行扰动；具体的扰动范围设置为：表层横波速度的扰动范围从40m/s到初始模型横波速度的1000%，表层以下横波速度的扰动范围是初始模型的横波速度的
±
50%。
15.更进一步地，利用瑞雷面波相速度频散确定横波速度的初始模型，以初始模型为基准，在允许范围内进行扰动获得最终的反演模型的步骤，包括：不同频率的瑞雷面波相速度对横波速度的敏感深度不同，根据经验公式利用不同频率的瑞雷面波相速度确定不同深度的横波速度的初始模型。以初始模型为起点，进行满足高斯概率分布的随机扰动得到新模型，接着利用metropolis准则对新模型的质量进行评价。如果新模型被接受，则以新模型作为基准进行后续的扰动；如果新模型被舍弃，则仍然以之前的模型为基准进行扰动。最后，将所有被接受的模型的平均值作为最终的反演模型，而所有被接受的模型的标准差作为反演结果的不确定度。
16.与现有技术相比，本发明的有益效果：本发明通过对比分析瑞雷面波相速度和h/v对于不同深度的横波速度的敏感度，
发现瑞雷面波h/v相较于相同频率的瑞雷面波相速度对于横波速度具有明显更浅的敏感深度。因此本发明提出了一种基于理论和实例测试验证了从主动源地震记录提取并联合反演面波相速度和h/v的技术的可靠性，证明本方案能更准确确定横波速度结构；并且最后讨论了所提出的技术对于在复杂地表环境下开展地震探测的潜力。
附图说明
17.为了更清楚地说明本发明实施例的技术方案，下面将对实施例中所需要使用的附图作简单地介绍， 应当理解，以下附图仅示出了本发明的某些实施例，因此不应被看作是对范围的限定，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他相关的附图。
18.图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以浅的结果。
具体实施方式
19.下面将结合本发明实施例中附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。通常在此处附图中描述和示出的本发明实施例的组件可以以各种不同的配置来布置和设计。因此，以下对在附图中提供的本发明的实施例的详细描述并非旨在限制要求保护的本发明的范围，而是仅仅表示本发明的选定实施例。基于本发明的实施例，本领域技术人员在没有做出创造性劳动的前提下所获得的所有其他实施例，都属于本发明保护的范围。
20.应注意到：相似的标号和字母在下面的附图中表示类似项，因此，一旦某一项在一个附图中被定义，则在随后的附图中不需要对其进行进一步定义和解释。同时，在本发明的描述中，术语“第一”、“第二”等仅用于区分描述，而不能理解为指示或暗示相对重要性，或
者暗示这些实体或操作之间存在任何这种实际的关系或者顺序。
21.实施例1：本发明通过下述技术方案实现，基于主动源面波频散和h/v确定浅地表横波速度的方法，包括以下步骤：为了不让权利要求书和说明书一样，可以在以下每个步骤中给出技术实施时的实验例子，再配合实验数据等，就会形成一个对应于技术步骤的完整实施例。
22.步骤s1：利用高分辨率线性拉东变换(hrlrt)技术从主动源地震记录中分离得到基阶模式瑞雷面波。
23.为了获得瑞雷面波频散能量图，先利用高分辨率线性拉东变换(hrlrt)技术分别将垂直分量和径向分量的原始地震记录从时间-距离域转换到频率-慢度域，再利用插值的方式将面波数据从频率-慢度域转换到频率-速度域。然后从包含瑞雷面波的频散能量图，也就是频率-速度域中拾取基阶模式瑞雷面波，并将其反变换到时间-距离域，即获得了分离后的基阶模式瑞雷面波信号。
24.步骤s2：对分离后的基阶模式瑞雷面波应用反褶积和频率-时间分析技术测量每一道记录的瑞雷面波的h/v和瑞雷面波在任意两地震道之间传播的相旅行时。
25.基于高分辨率线性拉东变换(hrlrt)技术得到的分离后的基阶模式瑞雷面波记录，将其中每一道的垂直分量的记录分别与其他所有道的垂直分量和径向分量的记录进行反褶积计算，获得所有任意两道间的反褶积函数。所述反褶积函数可以被认为是一个地震道作为虚源，而另一个地震道作为接收点观测的信号。
26.利用频率-时间分析(ftan)技术测量两道之间反褶积函数中瑞雷面波的振幅以及相、群旅行时。频率-时间分析(ftan)技术基于不同的中心频率进行一系列的高斯窄带滤波，接着对每一个高斯窄带滤波后的波形做希尔伯特变换得到信号的包络。通过寻找包络的最大值就可以确定瑞雷面波的振幅，通过测量振幅最大值所在时刻即为瑞雷面波的群旅行时。再通过测量振幅最大值所在时刻的瞬时相位和瞬时频率，即可确定瑞雷面波的相旅行时。
27.对于任意两道，第一道的径向(r)分量-第二道的垂直(z)分量和第一道的垂直(z)分量-第二道的垂直(z)分量的反褶积函数中瑞雷面波信号的振幅比，作为第一道的瑞雷面波h/v；或计算第一道的垂直(z)分量-第二道的径向(r)分量和第一道的垂直(z)分量-第二道的垂直(z)分量的反褶积函数中瑞雷面波信号的振幅比，作为第二道的瑞雷面波h/v。对于任意两道，舍弃低信噪比的反褶积函数或不同分量的瑞雷面波相、群旅行时不一致的反褶积函数。
28.最后，对于每一道，将剩余的h/v测量结果的平均值作为最终的h/v测量结果，标准差作为h/v测量结果的不确定度。然而，瑞雷面波的h/v会受到复杂结构或者方位各向异性的影响。而本方案基于一维测线的观测系统难以充分考虑这些因素的影响，会导致测量的h/v的不确定度偏小。考虑到具有良好方位角覆盖的天然地震或者基于二维台阵观测的背景噪声互相关函数测量的h/v的不确定度通常小于0.02，因此在本方案中如果h/v测量结果的不确定度小于0.025，则统一将其设为0.025。
29.步骤s3：基于任意两道之间的相旅行时测量结果，利用基于直射线理论的层析成像技术计算位于测线上每一个离散网格的瑞雷面波相速度。
30.基于由垂直(z)分量-垂直(z)分量的反褶积函数测量的瑞雷面波信号在任意两道间的相旅行时，采用直射线理论的面波层析成像技术计算不同频率的瑞雷面波相速度的横向变化。在正演时，将一维测线剖分为若干离散的网格；对于一个特定频率f，给定路径的瑞雷面波的相旅行时表达为该射线穿过每一个网格所需的相旅行时的叠加。
31.因此穿过研究区域的所有射线的面波相旅行时由以下公式表示：d=gm，其中g表示关于射线路径的矩阵，d表示相旅行时的矩阵，m表示与频率f相关的慢度矩阵；通过算子g将模型空间m投影到数据空间d。
32.以此为正演过程，可以通过反演计算得到不同频率的相速度分布。将反演过程中使用的目标函数定义为，通过一系列的检测板测试确定离散网格合适的大小，并利用l-曲线寻找最优的阻尼因子。
33.从瑞雷面波频散能量图中拾取的面波相速度频散的不确定度通常被定义为频散能量的中值的半宽度。因为面波层析成像技术具有较强的抗噪声能力，并且相对于瑞雷面波频散能量图能更好的分辨相速度水平方向的变化，所以本方案将相速度的不确定度同一设置为由频散能量图拾取的相速度的不确定的0.5倍。
34.步骤s4：利用马尔科夫链-蒙特卡洛算法(mcmc)联合反演瑞雷面波相速度和h/v，获得每一个网格下方的横波速度结构，并构建拟二维横波速度剖面。
35.马尔科夫链-蒙特卡洛算法(mcmc)在足够的随机搜索次数的前提下，可以搜索到全局极小值并获得反演结果的不确定度。
36.为了进行联合反演，首先将每一道的h/v测量结果进行插值得到不同网格的h/v，网格的大小与层析成像中所使用的网格大小相同。使用马尔科夫链-蒙特卡洛算法(mcmc)联合反演瑞雷面波的相速度频散和h/v，获取每一个网格下方的一维横波速度结构。
37.在本步骤中，通过瑞雷面波相速度频散确定横波速度的初始模型，以初始模型为基准，在允许范围内进行扰动得到最终的横波速度的反演模型。模型的扰动范围设置如下：表层横波速度的扰动范围从40m/s到初始模型横波速度的1000%；表层以下横波速度的扰动范围是初始模型横波速度的
±
50%。
38.具体的反演过程为：根据经验公式利用不同频率的瑞雷面波相速度计算不同深度的横波速度的初始模型。以初始模型为起点，进行满足高斯概率分布的随机扰动得到新模型，接着利用metropolis准则对新模型的质量进行评价。如果新模型被接受，则以新模型作为基准进行后续的扰动；如果新模型被舍弃，则仍然以之前的模型为基准进行扰动。最后，将所有被接受的模型的平均值作为最终的反演模型，而所有被接受的模型的标准差作为反演结果的不确定度。
39.一般而言，由于地层的压实效应，横波速度随着深度的增加而增加。因此本方案增加了在深度1.5m以下横波速度随深度增加而增加的限制，以减小反演结果的不确定度。需要注意的是增加了这个约束如果导致基于反演结果的预测值不能很好的拟合观测数据，则应该取消该约束。因为横波速度是瑞雷面波频散和h/v最敏感的参数，为了使反演过程更稳定，在反演中假设不同深度的纵波速度和密度是已知的，而只反演得到横波速度。
40.实施例2：基于实施例1的技术方案，本实施例进行相应的理论测试。
41.基于有限差分算子计算二维p-sv波动方程得到所有垂直分量和径向分量的理论
合成数据。其中，震源子波函数被定义为雷克子波，主频为15hz。为了保证有限差分求导的稳定性，理论模型被均匀剖分成0.05
×
0.05m的网格，时间域上的步长为0.1ms。
42.本实施例使用一个两层理论模型测试瑞雷面波h/v测量方法的鲁棒性。为了模拟野外记录的地震波场，在这个两层理论模型的合成数据中添加了高斯随机噪声。其中，高斯随机噪声的平均值是0，方差定义为第一道垂直分量理论合成数据的最大振幅的2%。基于加噪后的理论数据，利用hrlrt方法计算频散能量图，并从中分离出基阶模式瑞雷面波。基于分离得到的基阶模式瑞雷面波，使用任意两道间反褶积函数和ftan技术测量每一道的瑞雷面波h/v。
43.请参见图2，其中图2（a）为基于两层理论模型分别计算的频率为20hz、40hz、60hz的瑞雷面波相速度对不同深度的横波速度的敏感度，图2（b）为基于两层理论模型分别计算的频率为20hz、40hz、60hz的瑞雷面波h/v对不同深度的横波速度的敏感度。
44.作为一个例子，图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。
45.更进一步地，本实施例使用地表为高速层的阶梯状模型模拟城市地表环境下的复杂结构。根据提出的技术方案，计算瑞雷面波的相速度频散和h/v，然后利用mcmc算法分别进行瑞雷面波相速度频散的单独反演，以及瑞雷面波相速度频散和h/v的联合反演，获得地下横波速度结构，请参见图4基于地表为高速层的阶梯状模型的合成记录的测试结果。
46.结果表明，请参见图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）的横波速度则稍大于理论模型的横波速度。
47.总体而言，基于瑞雷面波相速度频散单独反演得到的横波速度模型与理论模型差异较大。与之相反，本方案通过引入h/v进行联合反演，基于反演结果的预测相速度和h/v都能很好的拟合观测值，请参见图5（a）、图5（b）。并且反演得到的横波速度模型与理论横波速度模型在不同深度都非常接近，请参见图5（c）、图5（d），参数表示同图4。
48.需要注意的是，联合反演得到的地表横波速度和理论横波速度非常接近，且反演结果的不确定度非常小（《100m/s）。理论测试表明通过引入h/v进行联合反演可以显著改善横波速度反演结果的准确度。
49.以上所述，仅为本发明的具体实施方式，但本发明的保护范围并不局限于此，任何熟悉本技术领域的技术人员在本发明揭露的技术范围内，可轻易想到变化或替换，都应涵盖在本发明的保护范围之内。因此，本发明的保护范围应所述以权利要求的保护范围为准。