一种基于结构耦合约束的体波和面波三维联合反演方法

技术领域

本发明属于地震信息处理技术领域，尤其涉及一种基于结构耦合约束的体波和面波三维联合反演方法。

背景技术

地震成像是研究地球内部结构的最重要工具之一，通过地震成像获取的三维速度模型为油气矿产资源探测及开发、大地震的深部孕震及发震机理、地质灾害监测及预警、地球内部结构及深部动力学过程等领域提供关键信息。随着探测技术不断发展，野外地震观测也逐渐往多参数、智能化、高精度、超密集台阵探测等方向发展，海量的、高精度的地震数据不断积累，高精度成像技术也不断推进。体波和面波是地震记录上两种重要的信号，充分地利用这两种数据体可以有效的提高反演精度。

基于射线理论的反演方法是体波和面波成像的重要手段，因其数据资料处理简单、正反演技术成熟等优势，广泛应用于地球不同尺度介质非均匀性特征研究。目前体波和面波数据反演主要以单方法反演为主，然而单方法反演往往受限于方法本身的局限，成像分辨率受到限制。体波反演精度受限于地震和台站的分布，在地震发生方位不均匀或台站分布稀疏的情况下，成像分辨率有限。而面波反演方法的纵向分辨率随深度增加而降低、对介质P波速度的灵敏度低于S波速度等方面的约束，仅靠面波单方法反演难以获得高精度的速度结构。

已有的体波和面波联合思路是将体波和面波反演过程中的灵敏核函数构建成一个大型线性稀疏矩阵，通过选取两种地震数据的权重等参数，同步求解出P波和S波的速度结构。该联合反演策略是在同一个反演系统中，设定在同一迭代终止条件，联合反演体波走时和面波频散数据，同步获取地下介质的P波和S波速度结构。然而该方法需要选取合适的权重因子，才能保证在同一的评价体系下中获得最优的P波和S波速度结构，确保其联合反演结果优于其单方法反演结果。权重因子的选择不仅仅包含有不同数据体质量控制权重因子、不同数据体之间权重因子、正则化项权重因子等，如何选取合适的权重因子使得联合反演获取的速度结构优于单方法反演结果给该方法带来了很大的困难。

发明内容

针对现有技术中的上述不足，本发明提供的一种基于结构耦合约束的体波和面波三维联合反演方法克服了传统联合反演两种不同数据体权重因子选择困难的问题，提高了单方法反演成像分辨率。

为了达到上述发明目的，本发明采用的技术方案为：一种基于结构耦合约束的体波和面波三维联合反演方法，包括以下步骤：

S1、获取体波走时数据和面波频散数据；

S2、进行网格剖分，得到剖分结果；

S3、构建第一初始模型和第二初始模型；

S4、根据剖分结果对第一初始模型进行体波走时数据反演，得到体波单方法反演模型；

S5、根据剖分结果对第二初始模型进行面波频散数据反演，得到面波频散单方法反演模型；

S6、将第一初始模型更新为体波单方法反演模型，以及将第二初始模型更新为面波频散单方法反演模型，并判断是否达到初始迭代次数，若是，则进入步骤S7，否则，返回步骤S4；

S7、根据第一初始模型构建面波频散约束反演模型，以及根据第二初始模型构建体波约束反演模型，并基于体波约束反演模型和面波频散约束反演模型完成联合反演。

本发明的有益效果为：本发明充分利用了地震数据，并加入双差走时提高了震源区域速度结构精度；利用交叉梯度函数构建体波约束反演模型和面波频散约束反演模型，并在面波频散单方法反演模型的目标函数的构建中带入体波约束反演模型，在体波单方法反演模型的目标函数的构建中带入面波频散约束反演模型，实现两种方法的约束反演，采用交替约束策略联合反演，最终获取在结构上越来越相似的两种速度模型。交替约束联合反演策略可以在联合反演过程中不断更新约束模型，相比直接用固定模型做约束反演，避免了先验约束模型精度对反演结果的影响。

进一步地，所述步骤S2网格剖分具体为中间目标区域网格间距为0.1°，深度网格在26km内网格间距为2km，边界网格采用不均匀剖分。

上述进一步方案的有益效果为：网格划分考虑了探测目标体的规模以及数据覆盖范围，平衡了正演和反演的计算精度。

进一步地，所述步骤S4具体为：

S401、根据剖分结果对第一初始模型进行正演计算，得到体波理论模型正演数据：

其中，

S402、将体波理论模型正演数据和体波走时数据作差，得到体波实测走时数据与理论模型走时数据的差：

其中，

S403、根据体波实测走时数与理论模型走时数据的差，获取所有地震事件及其对应的接收地震台站的体波走时差，得到数据走时残差矩阵：

ΔT＝AΔX+CΔU

其中，ΔT为数据走时残差矩阵；A为所有震源参数的扰动组合成的大型矩阵；ΔX为所有震源的空间扰动位置组成的一维矩阵；C为所有地震事件传播路径矩阵；ΔU为所有网格的慢度扰动组成的一维矩阵；

S404、在数据走时残差矩阵中加入双差矩阵和圆滑因子，得到体波单方法反演模型，所述体波单方法反演模型的目标函数为：

其中，L

上述进一步方案的有益效果为：将通过不断修正模型，使得体波理论模型正演数据与体波走时数据越来越逼近，获取满足要求的速度结构，并加入双差走时，提高震源区域速度结构的精度。

进一步地，所述双差矩阵具体为双差走时组成的矩阵；所述双差走时为两个位置相近的震源到同一观测台的实测走时和理论模型走时的时差做差。

上述进一步方案的有益效果为：提高了震源区域速度结构的精度。

进一步地，所述步骤S5具体为：

S501、根据剖分结果对第二初始模型进行正演计算，得到面波理论模型正演数据：

其中，

S502、将面波理论模型正演数据和面波频散数据作差，得到面波实测频散走时数据和理论模型正演走时数据的差：

其中，δt

S503、根据面波实测频散走时数据和理论模型正演走时数据的差，获取所有地震事件及其对应的接收地震台站的面波走时差，得到面波频散走时差矩阵：

Δd

其中，Δd

S504、在面波频散走时差矩阵中加入模型光滑项，得到面波频散单方法反演模型，所述面波频散单方法反演模型的目标函数为：

其中，L

上述进一步方案的有益效果为：将通过不断修正模型，使得面波理论模型正演数据与面波频散数据越来越逼近，获取满足要求的速度结构，在走时差矩阵中加入模型光滑项，提高反演结果的稳定性。

进一步地，所述步骤S7具体为：

S701、构建并根据交叉梯度函数，分别利用第一初始模型和第二初始模型作为约束构建面波频散约束反演目标函数和体波约束反演目标函数；

S702、根据面波频散约束反演目标函数和体波约束反演目标函数分别得到面波频散约束反演模型和体波约束反演模型；

S703、采用交替迭代策略进行联合约束反演，判断是否满足迭代终止条件，若是，则得到P波速度模型和S波速度模型，否则，将第一初始模型更新为体波约束反演模型，以及将第二初始模型更新为面波频散约束反演模型，并返回步骤S701。

上述进一步方案的有益效果为：利用体波约束反演模型和面波频散约束反演模型，在体波单方法反演模型的目标函数中带入面波频散约束反演模型，实现体波约束反演，以及在面波频散单方法反演模型的目标函数中带入体波约束反演模型，实现面波频散约束反演，采用交替约束反演策略联合反演体波走时数据和面波频散数据，交替约束联合反演后获取的两种模型在结构上也越来越相似，相比直接用固定模型做约束反演对约束模型的依赖性小一些，避免了先验约束模型精度对反演结果的影响。

进一步地，所述步骤S701具体为：

S7011、构建交叉梯度函数：

其中，D为交叉梯度函数；W

S7012、根据体波单方法反演模型的目标函数和交叉梯度函数，以第二初始模型为约束构建体波约束反演目标函数：

其中，L

S7013、根据面波频散单方法反演模型和交叉梯度函数，以第一初始模型为约束构建面波频散约束反演目标函数：

其中，L

上述进一步方案的有益效果为：利用交叉梯度函数构建体波约束反演目标函数和面波频散约束反演目标函数，并在体波单方法反演模型的目标函数中带入面波频散约束反演目标函数，在面波频散单方法反演模型的目标函数中带入体波约束反演目标函数，采用交替约束策略联合反演体波走时数据和面波频散数据，经交替约束联合反演后获取的两种模型在结构上也越来越相似，相比直接用固定模型做约束反演对约束模型的依赖性小一些，避免了先验约束模型精度对反演结果的影响。

进一步地，所述步骤S704中迭代终止条件为满足数据拟合差的最小值，或相邻两次迭代的第一初始模型和第二初始模型之间的值差满足设定的阈值，或满足设定的最大迭代次数。

上述进一步方案的有益效果为：设定迭代终止条件，为联合反演提供结束出口，保证最终P波速度模型和S波速度模型的获取。

附图说明

图1为本发明的方法流程图。

图2为本发明提供的三维正、反演网格剖分示意图。

图3为本发明实施例提供的三维理论模型示意图。

图4为本发明实施例提供的联合反演后P波速度和S波速度交叉梯度分布图。

具体实施方式

下面对本发明的具体实施方式进行描述，以便于本技术领域的技术人员理解本发明，但应该清楚，本发明不限于具体实施方式的范围，对本技术领域的普通技术人员来讲，只要各种变化在所附的权利要求限定和确定的本发明的精神和范围内，这些变化是显而易见的，一切利用本发明构思的发明创造均在保护之列。

如图1所示，在本发明的一个实施例中，一种基于结构耦合约束的体波和面波三维联合反演方法，包括以下步骤：

S1、获取体波走时数据和面波频散数据；

S2、进行网格剖分，得到剖分结果；

S3、构建第一初始模型和第二初始模型；

S4、根据剖分结果对第一初始模型进行体波走时数据反演，得到体波单方法反演模型；

S5、根据剖分结果对第二初始模型进行面波频散数据反演，得到面波频散单方法反演模型；

S6、将第一初始模型更新为体波单方法反演模型，以及将第二初始模型更新为面波频散单方法反演模型，并判断是否达到初始迭代次数，若是，则进入步骤S7，否则，返回步骤S4；

S7、根据第一初始模型构建面波频散约束反演模型，以及根据第二初始模型构建体波约束反演模型，并基于体波约束反演模型和面波频散约束反演模型完成联合反演。

本实施例中，体波走时数据可以从区域和地方地震台网、流动密集台阵的连续地震波形数据中获取。这两种数据来源对体波走时的提取方式不一样。对于区域和地方地震台网观测数据，选取观测报告中震级大等于1.0且走时残差小等于0.5s的震相数据。对于密集流动地震台阵观测数据，采用统软件拾取波形初至。将所有走时数据画时距曲线，剔除走时曲线中离散较大的震相，确保每个地震至少有5个台站记录到P波初至到时。

本实施例中，面波处理所需的连续地震波形资料来源与体波走时资料相同，但处理手段不一样。具体为：对连续数据重采样、去仪器响应、去均值、去趋势、滤波；对时间序列进行谱白化和时间归一化处理；将连续信号按天或者小时分割，对相同时间段数据进行互相关计算，获取该时段台站对之间的互相关函数。然后将同一台站对不同时间段的互相关函数进行叠加，获取该台站对的经验格林函数；对经验格林函数进行滤波、归一化、时间窗选取等处理。采用多窗口测量法获取理论格林函数和经验格林函数不同频段的走时残差。

所述步骤S2网格剖分具体为中间目标区域网格间距为0.1°，深度网格在26km内网格间距为2km，边界网格采用不均匀剖分。

本实施例中，网格剖分是正、反演计算精度的关键因素之一。密集网格可以提高正演精度，但会增加反演的多解性。网格稀疏会导致正演精度降低，从而影响反演精度。同时网格划分还需要重点考虑探测目标体的规模以及数据覆盖范围。图2为本发明提供的三维正、反演网格剖分示意图，同时考虑了台站、地震分布、分辨率等需求进行网格剖分，在网格剖分方面中间目标区网格间距为0.1°，深度网格在26km内网格间距为2km，边界网格采用不均匀剖分。

所述步骤S4具体为：

S401、根据剖分结果对第一初始模型进行正演计算，得到体波理论模型正演数据：

其中，

S402、将体波理论模型正演数据和体波走时数据作差，得到体波实测走时数据与理论模型走时数据的差：

其中，

S403、根据体波实测走时数与理论模型走时数据的差，获取所有地震事件及其对应的接收地震台站的体波走时差，得到数据走时残差矩阵：

ΔT＝AΔX+CΔU

其中，ΔT为数据走时残差矩阵；A为所有震源参数的扰动组合成的大型矩阵；ΔX为所有震源的空间扰动位置组成的一维矩阵；C为所有地震事件传播路径矩阵；ΔU为所有网格的慢度扰动组成的一维矩阵；

S404、在数据走时残差矩阵中加入双差矩阵和圆滑因子，得到体波单方法反演模型，所述体波单方法反演模型的目标函数为：

其中，L

本实施例中，正演计算是本项目开发的算法的重要过程之一。在理论模型测试中，需要通过正演计算获取模型的正演数据，并将理论模型正演数据与观测数据做差，带入反演计算流程。在反演计算中需要通过计算每次反演模型的正演数据，计算理论数据和观测数据的走时差，通过走时差用于评估反演效果。以下简要介绍体波和面波正演计算过程。基于射线理论的走时成像正演计算理论可简化成射线从源(虚拟源)到地震台站的旅行事件。体波和面波走时的计算区别在于慢度的计算。体波慢度的获取是介质的速度的倒数，通过直接获取介质的速度物性即可。而面波慢度的计算需要对面波进行分频段滤波，分别对每个频段计算地表网格面波传播速度，对地表网格速度求倒数即可获得面波慢度。

本实施例中，采用的地震体波主要针对的是近震体波，即震源在研究区域之内。对于震源在研究区外的事件，可以通过射线追踪到研究区域内的边界网格上，将边界网格交点作为源，然后当成区域内的震源进行处理。

本实施例中，基于射线理论的体波反演问题跟正演问题中的已知量和未知量刚好相反，在反演过程中，采用一阶平滑、二阶平滑或高斯平滑等方式进行模型光滑可以提高反演的稳定性和反演结果的可靠性；对于每个地震信息的ΔX(所有震源的空间扰动位置组成的一维矩阵)是相互独立的，不能进行光滑处理，此处处理方式为补0。

所述双差矩阵具体为双差走时组成的矩阵；所述双差走时为两个位置相近的震源到同一观测台的实测走时和理论模型走时的时差做差。

所述步骤S5具体为：

S501、根据剖分结果对第二初始模型进行正演计算，得到面波理论模型正演数据：

其中，

S502、将面波理论模型正演数据和面波频散数据作差，得到面波实测频散走时数据和理论模型正演走时数据的差：

其中，δt

S503、根据面波实测频散走时数据和理论模型正演走时数据的差，获取所有地震事件及其对应的接收地震台站的面波走时差，得到面波频散走时差矩阵：

Δd

其中，Δd

S504、在面波频散走时差矩阵中加入模型光滑项，得到面波频散单方法反演模型，所述面波频散单方法反演模型的目标函数为：

其中，L

本实施例中，模型光滑项可以采用一阶、二阶、三阶、高斯等平滑公式，模型光滑项可以提高反演的稳定性和反演结果的可靠性。

所述步骤S7具体为：

S701、构建并根据交叉梯度函数，分别利用第一初始模型和第二初始模型作为约束构建面波频散约束反演目标函数和体波约束反演目标函数；

S702、根据面波频散约束反演目标函数和体波约束反演目标函数分别得到面波频散约束反演模型和体波约束反演模型；

S703、采用交替迭代策略进行联合约束反演，判断是否满足迭代终止条件，若是，则得到P波速度模型和S波速度模型，否则，将第一初始模型更新为体波约束反演模型，以及将第二初始模型更新为面波频散约束反演模型，并返回步骤S701。

所述步骤S701具体为：

S7011、构建交叉梯度函数：

其中，D为交叉梯度函数；W

S7012、根据体波单方法反演模型的目标函数和交叉梯度函数，以第二初始模型为约束构建体波约束反演目标函数：

其中，L

S7013、根据面波频散单方法反演模型和交叉梯度函数，以第一初始模型为约束构建面波频散约束反演目标函数：

其中，L

本实施例中，从Gallardo和Meju(2003)提出的不同物性参数的结构耦合关系式出发，所述结构耦合关系式为：

其中，m

公式(1)中m可以是m

其中，W

公式(3)即为交叉梯度函数离散化的线性方程组。本发明将公式(3)加入体波和面波单方法反演目标函数中，即可实现约束反演的目标函数构建.

所述步骤S704中迭代终止条件为满足数据拟合差的最小值，或相邻两次迭代的第一初始模型和第二初始模型之间的值差满足设定的阈值，或满足设定的最大迭代次数。

本实施例中，本发明采用交替迭代的策略实现体波走时和面波频散的联合反演。该策略是以约束反演目标函数为主线，在两种方法迭代过程中，将两种方法反演迭代过程中获取的更新模型作为约束模型去约束另一种方法，直至满足迭代终止条件。本发明迭代终止条件主要采用三种判断：(1)满足数据拟合差的最小值；(2)相邻两次迭代的模型值差满足满足设定的阈值；(3)满足设定的最大迭代次数。在迭代过程中，随着每种方法获取的模型越来越接近于真实模型，带入另一种方法的约束模型也越来越可靠，交替约束联合反演后获取的两种模型在结构上也越来越相似，相比直接用固定模型做约束反演对约束模型的依赖性小一些，避免先验约束模型精度对反演结果的影响。

本实施例中，为了考察本发明的算法对不同规模、埋深、异常强度的异常体的分辨率能力，本项目设计多种理论物性模型(包括规模不同、异常强度不同、埋深不同的单个或多个异常体)，通过理论模型合成数据进行反演试算，理论模型即为第一初始模型和第二初始模型的原始状态，所述理论模型极为地震波的模拟模型，第一初始模型和第二初始模型未开始反演前均为理论模型。图3为本发明实施例提供的三维理论模型示意图。该模型设计思路是为了充分发挥体波和面波的优势，浅部同一深度设计了4个规模一样的低速异常体(异常值低于背景速度3％)，用于考察两种方法横向分辨率的差异。在四个异常体下方设计了一个规模较大的高速异常体(异常值高于背景速度3％)，用于考察两种方法的垂向分辨率的差异。

本实施例中，由于反演迭代开始时，P波速度和S波速度模型跟初始模型相近，交叉梯度约束的作用小，因此刚开始迭代的时候，先进行地震体波走时和面波频散的单方法反演，当经过一定次数迭代后，启动两种方法的约束反演模块，将单方法反演的结果作为约束，约束另一种方法，在迭代更新过程中，两种方法的约束模型交替更新，实现两种方法的约束联合反演。

本实施例中，图4为联合反演后P波速度模型和S波速度模型的交叉梯度分布图，图中显示了交叉梯度值主要集中在这两种物性边界区域，这也间接说明在约束反演过程中，边界信息会通过交叉梯度项影响目标函数，从而影响反演结果，获得两种物性结构上更相似的模型。在一些结构不相似或物性没有变换的区域交叉梯度值几乎为零，在反演过程中，就不会影响目标函数的求解，从而保留了单方法的优势。这也是基于交叉梯度实现联合反演的目的，在联合反演过程中寻求两种物性相似的区域，使其结构更为耦合，最终获取形态上更为相似的多种物理模型。

本发明通过提供了一种地震体波走时和面波频散的联合反演方法，该方法不仅可以充分利用已有的地震数据，发挥不同地震数据体的优势，获取更高精度的速度结构，推进多种地震数据融合的高精度反演技术。同时其获取的高精度P波和S波速度模型还可以为油气矿产资源探测及开发、大地震的深部孕震及发震机理、地质灾害监测及预警、地球内部结构及深部动力学过程等领域提供关键信息。