一种地下结构弹性波正演声波反演成像方法及设备

1.本发明涉及地球物理勘探技术领域，特别涉及一种地下结构弹性波正演声波反演成像方法及设备。


背景技术：

2.地震勘探属于地球物理勘探领域，其探测结果具有较高的分辨率和可靠性，被广泛应用于石油勘探、矿产勘察、工程勘探和自然灾害防护等领域。地震波全波形反演技术是地震勘探领域中一种高精度反演方法，能有效刻画地下地质构造的细节，具有反演精度高和计算成本大的特点。近年来，随着计算计算力的整体提升，计算成本下降，全波形反演技术被逐渐应用于各类实际勘探工作中。
3.全波形反演方法的实现流程大致可分为以下几个步骤：
4.1)、给出初始模型，计算模拟信号，利用模拟信号和输入信号计算伴随源；
5.2)、将伴随源作为震源信号反传，得到反传波场，反传波场与正传波场做互相关，得到用于更新模型的梯度项；
6.3)、选择梯度类最优化方案(如最速下降法、共轭梯度法等)，更新模型直至收敛。
7.在实际地震勘探中，输入信号一般由位于地表的检波器观测得到。常见的单分量检波器接收的信号属于单分量弹性波信号(垂向速度分量)。多分量检波器由于成本问题，目前尚未被广泛应用。弹性波信号中含有较强的面波信息，这部分信息在基于体波的全波形反演中被视为一种强干扰信号，不利于反演。
8.目前，应用于实际勘探领域的全波形反演方法大多都基于声波方程。基于弹性波方程的方法由于计算成本太大，而达不到实际生产的需求。在得到观测信号后，一般需要将信号进行一系列预处理，例如归一化、去噪、面波压制等，提升信号的信噪比的同时使其更接近于声波信号，方便后续利用声波方程进行反演。目前常见的声波全波形反演方案还存在以下不足：
9.1)、实际观测信号属于弹性波信号，利用声波进行近似存在一定误差；
10.2)、观测信号含有面波信号，声波信号则不含面波信息；这部分缺失的面波信息会导致反演时两个信号间波形(面波、直达波、反射波等)的能量对比关系存在差异，难以正确恢复信号的幅值信息。
11.因此，现有技术还有待于改进和发展。


技术实现要素：

12.鉴于上述现有技术的不足，本发明的目的在于提供一种地下结构弹性波正演声波反演成像方法及设备，旨在解决现有地下结构成像方法存在计算成本大以及反演结果的准确性和精度不足的问题。
13.本发明的技术方案如下：
14.一种地下结构弹性波正演声波反演成像方法，其中，包括步骤：
15.给定初始模型，通过弹性波方程正演计算得到模拟信号，所述初始模型包括纵波速度模型，横波速度模型和密度模型；
16.根据所述模拟信号和输入信号计算得到伴随源；
17.通过声波方程计算反传波场和正传波场，并将所述反传波场和所述正传波场进行互相关，得到用于模型更新的梯度项；
18.根据所述梯度项对所述初始模型进行更新，直至最小二乘形式的目标函数收敛，输出反演成像结果。
19.所述的地下结构弹性波正演声波反演成像方法，其中，通过弹性波方程正演计算得到模拟信号的步骤包括：
20.在弹性波方程的应力项上加载子波震源信号，对波场传播进行数值模拟，记录检波器位置处弹性波场的垂向速度分量，所述垂向速度分量等价于模拟信号。
21.所述的地下结构弹性波正演声波反演成像方法，其中，所述弹性波方程为：
[0022][0023]
其中，v
x
(x,t)表示t时刻位于空间位置x处的粒子在x方向上的振动速度；ρ(x)为介质密度；τ
xx
(x,t)和τ
zz
(x,t)分别表示水平和垂向的正应力项；τ
xz
(x,t)为切应力项；λ和μ为拉梅常数，可通过模型的纵波速度(v
p
)、横波速度(vs)和密度ρ，按公式μ＝ρv
s2
和计算得到；f(xs,t)为外加的震源项；vz(x,t)为弹性波的垂向速度分量。
[0024]
所述的地下结构弹性波正演声波反演成像方法，其中，根据所述模拟信号和输入信号计算得到伴随源的步骤为：
[0025]
将所述模拟信号与输入信号进行相减，得到伴随源，所述输入信号是指实际工作中由检波器观测到的地震信号，或是在给定一个假设的真实模型条件下，通过正演计算所得到的信号。
[0026]
所述的地下结构弹性波正演声波反演成像方法，其中，通过声波方程计算反传波场和正传波场的步骤包括：
[0027]
假设弹性波方程中横波速度为零vs＝0，且两个正应力项相等τ
xx
(x,t)＝τ
zz
(x,t)＝p(x,t)，推出声波方程为：
[0028][0029]
将子波震源作为正传波场的震源，由0时刻计算到t时刻，按时间正传计算得到正传波场；
[0030]
将所述伴随源作为反传波场的震源，由t时刻计算到0时刻，按时间反传计算得到反传波场。
[0031]
所述的地下结构弹性波正演声波反演成像方法，其中，将所述反传波场和所述正传波场进行互相关，得到用于模型更新的梯度项的步骤包括：
[0032]
将同一时刻t的正传波场和反传波场相乘，累加0到t所有时刻波场相乘的结果，得到用于模型更新的梯度项。
[0033]
所述的地下结构弹性波正演声波反演成像方法，其中，根据所述梯度项对初始模型进行更新，直至最小二乘形式的目标函数收敛，输出反演成像结果的步骤包括：
[0034]
基于所述梯度项，采用梯度类最优化方法对初始模型进行更新，每次更新后所述最小二乘形式的目标函数的值会下降，若迭代若干次后，所述目标函数值相对前一次迭代不再下降，则判定此时迭代收敛，终止迭代程序，得到更新模型；
[0035]
根据所述更新模型输出反演成像结果。
[0036]
所述的地下结构弹性波正演声波反演成像方法，其中，所述最小二乘形式的目标函数：其中，j为目标函数；u为模拟信号；d为输入信号；xr为第r个检波器的空间位置，r＝1,2,3
…
n；m为参数模型；t为信号某一时刻的时间，t∈[0,t]，目标函数越小，模拟信号越接近于观测信号，此时生成模拟信号的模型接近于实际地质模型。
[0037]
一种计算机可读存储介质，其中，所述计算机可读存储介质存储有一个或者多个程序，所述一个或者多个程序可被一个或者多个处理器执行，以实现本发明所述地下结构弹性波正演声波反演成像方法中的步骤。
[0038]
一种地下结构弹性波正演声波反演成像设备，其中，包括：处理器、存储器及通信总线；所述存储器上存储有可被所述处理器执行的计算机可读程序；
[0039]
所述通信总线实现处理器和存储器之间的连接通信；
[0040]
所述处理器执行所述计算机可读程序时实现如本发明所述地下结构弹性波正演声波反演成像方法中的步骤。
[0041]
有益效果：本发明提供了一种地下结构弹性波正演声波反演成像方法，其核心思想是在计算伴随源时采用弹性波方程，使得模拟信号能尽可能地接近输入信号，这样信号之间的各类波形信息能有一个清晰的对应关系，从而提升反演的准确性和精度；另一方面，在计算模型梯度时固定采用声波方程，从而保证整体计算成本并符合实际生产需求。
附图说明
[0042]
图1为一种地下结构弹性波正演声波反演成像方法的流程图；
[0043]
图2a为弹性marmousi纵波速度模型(vp)；
[0044]
图2b为弹性marmousi横波速度模型(vs)；
[0045]
图2c为弹性marmousi密度模型(ρ)；
[0046]
图3为海洋模型输入信号图；
[0047]
图4为初始模型(vp)；
[0048]
图5a为真实模型vs/vp比值图；
[0049]
图5b为真实模型ρ/vp比值图；
[0050]
图6a为声波模拟信号与输入信号对比图；
[0051]
图6b为弹性波模拟信号与输入信号对比图；
[0052]
图7a为海洋模型全声波方案反演结果；
[0053]
图7b为海洋模型声波-声波方案反演结果；
[0054]
图7c为海洋模型弹性-声波方案反演结果；
[0055]
图8为陆地模型输入信号图；
[0056]
图9为水平平滑后的陆地模型输入信号图；
[0057]
图10a为陆地模型全声波方案反演结果；
[0058]
图10b为陆地模型声波-声波方案反演结果；
[0059]
图10c为陆地模型弹性-声波方案反演结果；
[0060]
图11为一种地下结构弹性波正演声波反演成像设备的原理图。
具体实施方式
[0061]
本发明提供一种地下结构弹性波正演声波反演成像方法及设备，为使本发明的目的、技术方案及效果更加清楚、明确，以下对本发明进一步详细说明。应当理解，此处所描述的具体实施例仅仅用以解释本发明，并不用于限定本发明。
[0062]
在实际地震勘探中，输入信号一般由位于地表的检波器观测得，所述输入信号属于弹性波信号，因此采弹性波信号作为模拟信号拟合伴随源，应当比声波更为准确，同时也更能保证信号间各波形信息的对应关系。但基于弹性波方程的全波形反演方法计算量太大，因此在计算模型梯度项时，可以考虑采用声波方程以减少计算成本。
[0063]
基于此，本发明提供了一种地下结构弹性波正演声波反演成像方法，如图1所示，其包括步骤：
[0064]
s10、给定初始模型，通过弹性波方程正演计算得到模拟信号，所述初始模型包括纵波速度模型，横波速度模型和密度模型；
[0065]
s20、根据所述模拟信号和输入信号计算得到伴随源；
[0066]
s30、通过声波方程计算反传波场和正传波场，并将所述反传波场和所述正传波场进行互相关，得到用于模型更新的梯度项；
[0067]
s40、根据所述梯度项对所述初始模型进行更新，直至最小二乘形式的目标函数收敛，输出反演成像结果。
[0068]
具体来讲，本发明在计算伴随源时，将传统声波全波形反演方法所采用的声波方
程替换为了弹性波方程，这样设计的优点有：采用弹性波方程计算模拟信号将更接近于实际输入信号，拟合伴随源时更加准确；通过弹性波方程计算的模拟信号和输入信号一样都存在面波信息，信号间各波形信息的对应关系比采用声波更准确；本发明在计算梯度项时依然采用了声波方程，因此整体的计算成本不大，符合实际生产需求。
[0069]
具体来讲，在实际勘探中，全波形反演方法的纵波速度模型可以通过偏移速度分析，测井数据或其他历史勘探等前期工作得到。数值测试中通过对真实模型进行模糊化得到相应的初始模型，所述初始模型包括纵波速度模型，横波速度模型和密度模型。在本实施例中，所述横波速度模型和密度模型可根据纵波速度模型按一定比例给定，计算公式为：[vs]＝[v
p
]
·
[β1],[ρ]＝[v
p
]
·
[β2]，式中的β1＝vs/v
p
和β2＝ρ/v
p
为比例参数，大小与纵波速度模型相等。实际工作中，这两个参数矩阵可简单按通过地表附近的数据给定，或依靠测井数据给定更为准确的值。
[0070]
在一些实施方式中，通过弹性波方程正演计算得到模拟信号的步骤包括：在弹性波方程的应力项上加载子波震源信号，对波场传播进行数值模拟，记录检波器位置处弹性波场的垂向速度分量，所述垂向速度分量等价于模拟信号。在本实施例中，所述弹性波方程为二维一阶速度-应力形式的弹性波动方程：
[0071][0072]
其中，v
x
(x,t)表示t时刻位于空间位置x处的粒子在x方向上的振动速度；ρ(x)为介质密度；τ
xx
(x,t)和τ
zz
(x,t)分别表示水平和垂向的正应力项；τ
xz
(x,t)为切应力项；λ和μ为拉梅常数，可通过模型的纵波速度(v
p
)、横波速度(vs)和密度ρ，按公式μ＝ρv
s2
和计算得到；f(xs,t)为外加的震源项；vz(x,t)为弹性波的垂向速度分量。
[0073]
在本实施例中，将所述模拟信号u(xr,t,m)与输入信号d(xr,t)进行相减，得到伴随源，所述输入信号是指实际工作中由检波器观测到的地震信号，或是在给定一个假设的真实模型条件下，通过正演计算所得到的信号。所述模拟信号u(xr,t,m)等价于检波器位置处弹性波场的垂向速度分量vz(xr,t,m)。
[0074]
在一些实施方式中，通过声波方程计算反传波场和正传波场的步骤包括：假设弹性波方程中横波速度为零vs＝0，且两个正应力项相等τ
xx
(x,t)＝τ
zz
(x,t)＝p(x,t)，推出声波方程为：
[0075][0076]
将子波震源作为正传波场的震源，由0时刻计算到t时刻，按时间正传计算得到正传波场；
[0077]
将所述伴随源作为反传波场的震源，由t时刻计算到0时刻，按时间反传计算得到反传波场。
[0078]
在本实施例中，将同一时刻t的正传波场和反传波场相乘，累加0到t所有时刻波场相乘的结果，得到用于模型更新的梯度项。具体的，全波形方法中，利用伴随状态法(adjoint-state method)和波恩近似(bonn approximate)可以推导出梯度项的求解公式为：
[0079][0080]
，其中[u(xr,t,m)-d(xr,t)]为伴随源；g为格林函数，*为卷积符号；即梯度项等于正传波场u(x,t,m)与反传波场g*[u(xr,t,m)-d(xr,t)]之间的互相关结果。
[0081]
在一些实施方式中，基于所述梯度项，采用梯度类最优化方法对初始模型进行更新，每次更新后所述最小二乘形式的目标函数的值会下降，若迭代若干次后，所述目标函数值相对前一次迭代不再下降，则判定此时迭代收敛，终止迭代程序，得到更新模型；根据所述更新模型输出反演成像结果。
[0082]
具体来讲，可采用共轭梯度法或拟牛顿法按一定步长更新初始模型的纵波速度模型，对应横波速度模型和密度模型按上文所述的基于纵波速度模型乘以相应的比例参数给出。以最速下降法为例：出。以最速下降法为例：其中，mn(x)为第n次迭代的模型；α为迭代步长；为fr
é
chet导数，表示目标函数对于模型的梯度。
[0083]
在本实施例中，所述最小二乘形式的目标函数：其中，j为目标函数；u为模拟信号；d为输入信号；xr为第r个检波器的空间位置，r＝1,2,3
…
n；m为参数模型；t为信号某一时刻的时间，t∈[0,t]，目标函数越小，模拟信号越接近于观测信号，此时生成模拟信号的模型接近于实际地质模型。
[0084]
在本实施例中，当目标函数收敛后能输出纵波速度模型，若反演效果好，此时的纵波速度模型将接近于实际地下介质的纵波速度分布，可以依据图像和数值结合地质信息辅助判断地下介质和构造分布。实际工作中，程序的输入数据为检波器接收和记录到的信号
以及一个大致正确的初始纵波速度模型(可以通过偏移速度分析，测试数据，历史勘探记录等方法得到)，本发明则是根据地震波传播规律，矫正和优化这一初始模型，使其更准确和清晰，输出优化后的模型，从而可获得更精准的反演结果。
[0085]
下面通过具体实施例对本发明做进一步的解释说明：
[0086]
实施例1
[0087]
共测试和对比了三种全波形反演方案，分别是(1)声波信号作为输入信号，采用声波方程计算伴随源和梯度项，下文中简称为全声波方案；(2)弹性波信号作为输入信号，采用声波方程计算伴随源和梯度项，下文中简称为声波-声波方案；(3)弹性波信号作为输入信号，采用弹性波方程计算伴随源，声波方程计算梯度项，下文中简称为弹性-声波方案。其中，实际观测信号属于弹性波信号，全声波方案的输入信号是声波信号，因此作为理想情况下的对照组；声波-声波方案代表常规方法；弹性-声波方案代表改进方法。
[0088]
1、基础设置：
[0089]
1)、采用了seg公开的弹性marmousi模型(图2a-2c)作为假设的真实模型。原模型的网格大小为13601*2801，我们将其稀疏为681*141，设置网格间距为30m
[0090]
2)、震源和检波器都设置在模型的表层(水面或地表)，炮间距为600m，道间距为30m，共33个炮点，681个检波点；
[0091]
3)、采用4hz的雷克子波作为震源信号，加载于声波或弹性波方程的应力项上。检波器记录的输入信号为弹性波的垂向速度分量(vz)，记录的时间间隔为0.002s，总记录时长为4s。图3和图8分别展示了海洋和陆地的输入信号(第15炮)，分别代表有面波和无面波的两种情况；
[0092]
4)、将真实模型小幅平滑，得到用于反演的初始模型，如图4所示；
[0093]
5)、针对弹性-声波方案，利用初始纵波模型(vp)按一定比例生成横波(vs)和密度(ρ)模型。经测试，这一比例的可选择区间较大，但与真实模型偏离太多时依然会使反演结果变差。为了在实际工作中提供一个可行的参考方案，我们选择真实模型在地表附近的大致模型比例作为弹性-声波方案的生成vs和ρ模型的比例(这里采用vp:vs:ρ＝1:0.3:1.2)。图5a和图5b分别为真实模型vs/vp的值和ρ/vp的值。需要注意，真实模型每个点的比例(vp:vs:ρ)都不同，而弹性-声波方案生成的模型每个点的比例都是固定的，这增大了弹性-声波方案在实际工作中的可行性；
[0094]
2、模拟海洋勘探(无面波影响)
[0095]
1)、利用初始模型生成模拟信号，分别将声波和弹性波模拟信号与输入信号进行对比，对比结果如图6所示；
[0096]
2)、各方案的反演结果如图7所示。
[0097]
3、模拟陆地勘探(有面波影响)
[0098]
1)、将marmousi模型上部的水层去除，其他参数设置不变。陆地模型得到的输入信号(含面波)如图8所示
[0099]
2)、采用水平平滑压制面波，压制后效果如图9所示；
[0100]
3)、各方案反演结果的如图10所示。
[0101]
从模拟信号和输入信号的对比图(图6)和反演结果(图7)可以得出，弹性波比声波信号更接近于输入信号，反演结果也更清晰。
[0102]
陆地模型中弹性-声波方案的反演结果(图10c)明显优于声波-声波方案(图10b)。其主要原因是，声波没有面波信息，在输入信号中面波存在的区域，声波与弹性波的波形幅值信息难以匹配，导致反演结果较差。而弹性-声波方案则可以避免这一问题，压制面波后，其反演结果大体接近于全声波反演结果(图10a)。
[0103]
在一些实施方式中，还提供一种地下结构弹性波正演声波反演成像系统，其包括:模拟信号计算模块，用于给定初始模型，通过弹性波方程正演计算得到模拟信号；伴随源计算模块，用于根据所述模拟信号和输入信号计算得到伴随源；梯度项计算模块，用于通过声波方程计算反传波场和正传波场，并将所述反传波场和所述正传波场进行互相关，得到用于模型更新的梯度项；模型更新模块，用于根据所述梯度项对速度模型进行更新，直至最小二乘形式的目标函数收敛，输出反演成像结果。
[0104]
在一些实施方式中，还提供一种计算机可读存储介质，其中，所述计算机可读存储介质存储有一个或者多个程序，所述一个或者多个程序可被一个或者多个处理器执行，以实现本发明所述地下结构弹性波正演声波反演成像方法中的步骤。
[0105]
在一些实施方式中，本技术还提供了一种地下结构弹性波正演声波反演成像设备，如图11所示，其包括至少一个处理器(processor)20；显示屏21；以及存储器(memory)22，还可以包括通信接口(communications interface)23和总线24。其中，处理器20、显示屏21、存储器22和通信接口23可以通过总线24完成相互间的通信。显示屏21设置为显示初始设置模式中预设的用户引导界面。通信接口23可以传输信息。处理器20可以调用存储器22中的逻辑指令，以执行上述实施例中的方法。
[0106]
此外，上述的存储器22中的逻辑指令可以通过软件功能单元的形式实现并作为独立的产品销售或使用时，可以存储在一个计算机可读取存储介质中。
[0107]
存储器22作为一种计算机可读存储介质，可设置为存储软件程序、计算机可执行程序，如本公开实施例中的方法对应的程序指令或模块。处理器20通过运行存储在存储器22中的软件程序、指令或模块，从而执行功能应用以及数据处理，即实现上述实施例中的方法。
[0108]
存储器22可包括存储程序区和存储数据区，其中，存储程序区可存储操作系统、至少一个功能所需的应用程序；存储数据区可存储根据终端设备的使用所创建的数据等。此外，存储器22可以包括高速随机存取存储器，还可以包括非易失性存储器。例如，u盘、移动硬盘、只读存储器(read-only memory，rom)、随机存取存储器(random access memory，ram)、磁碟或者光盘等多种可以存储程序代码的介质，也可以是暂态存储介质。
[0109]
此外，存储介质以及终端设备中的多条指令处理器加载并执行的具体过程在上述方法中已经详细说明，在这里就不再一一陈述。
[0110]
应当理解的是，本发明的应用不限于上述的举例，对本领域普通技术人员来说，可以根据上述说明加以改进或变换，所有这些改进和变换都应属于本发明所附权利要求的保护范围。