一种振幅、相位信息可调节的全波形反演方法及装置

1.本发明涉及一种振幅、相位信息可调节的全波形反演方法及装置，属于地震勘探速度建模技术领域。


背景技术：

2.随着我国石油勘探向深层及超深层转移，深部目标层高陡复杂构造的高精度成像问题已成为当前制约我国石油勘探和开发迫在眉睫的瓶颈。速度建模是实现复杂构造成像的灵魂，而诸如速度层析反演等传统速度建模技术虽在过去的数十年中取得了巨大的成功，但已无法满足高分辨率地震成像的需求。因此，亟需发展新一代速度建模技术，以满足复杂构造成像的需求。
3.目前，全波形反演(full waveform inversion，fwi)是公认的最具潜力的新一代速度建模技术，在准确的初始模型下能够恢复出高分辨率的速度模型，可用于提高深层目标层偏移成像的质量，具有广阔的应用前景和潜在的经济效益。全波形反演以层析反演或者速度分析结果为初始模型，可实现高分辨率的速度建模，结合保幅的偏移成像方法，获得高质量的成像剖面，为地质解释提供支撑和依据，因此全波形反演成为近些年来的研究热点之一。然而，计算效率低下是制约全波形反演进行大规模计算以及实际资料处理的主要因素，现有的一种解决方法是使用震源编码策略，其思想是通过某种震源组合方式将很多单炮形成一个超级震源，从而减少波场模拟的次数。其中被广泛应用的一种为基于频率选择的同时源全波形反演方法，它可以获得无串扰的反演结果，但是这种方法如果在初始模型不是很理想时，由于波场振幅信息对初始模型比较敏感，使用全波场信息匹配很容易使得全波形反演陷入局部极小值，所以开展数据子集或者波场属性反演就显得十分需要。由于振幅信息对初始模型的敏感度比相位信息高，所以当初始模型不准确时，直接使用全波形反演，可能导致全波形反演低波数更新不准确，从而导致高波数更新的界面位置不对，降低反演的精度。
4.因此，如何提供一种新的同时源反演方法，能够克服现有技术存在的反演精度低的问题是本领域技术人员亟需解决的技术难题。


技术实现要素：

5.针对上述现有技术存在的问题，本发明提供一种振幅、相位信息可调节的全波形反演方法及装置，通过一系列单频子波同时激发形成混叠波场；并采用相敏检测法从混叠波场中解耦出每一个编码信号对应的波场变量，实现无串扰的梯度计算，从而有效提高反演精度。
6.为了实现上述目的，本发明采用的技术方案是：一种振幅、相位信息可调节的全波形反演方法，具体步骤为：
7.a、通过多个位于不同位置的单频子波同时激发作为震源，采用多个检波点组成观测系统，通过观测系统获得观测数据；建立速度场模型(该速度场模型从现有模型中选择一
个)，并根据震源情况计算得出模拟数据；
8.b、基于观测数据和模拟数据，构建表征所述观测数据与所述模拟数据之间数据残差的目标函数；
9.c、根据震源与检波点之间的位置，得出虚拟震源的信息，然后根据观测数据获取震源的正传波场；并根据目标函数获得虚拟震源的反传波场；
10.d、确定目标函数关于速度场模型参数的梯度计算表达式，在不同的反演阶段具有两种计算表达式，其中在反演初期即低频反演阶段，采用具有相位信息反演的梯度表达式，在反演的中后阶段，采用兼顾波场相位和振幅信息的全波形反演的梯度表达式；先对震源的正传波场和虚拟震源的反传波场分别进行无串扰解耦，获得解耦的波场变量；然后根据得到的波场变量和梯度计算表达式，计算得到速度场模型参数的梯度；反演初期和反演中后期具体划分为：多尺度全波形反演一般反演3至5个频带，通常前一到两个频带被认为是低波数恢复，也就是低频反演阶段，中间频带的反演恢复中波数信息，后续频带为高频反演。
11.e、采用共轭梯度方法求取更新方向，根据得到的梯度对速度场模型参数进行迭代更新，直至数据残差达到设定阈值或迭代更新次数达到预设值，停止迭代更新确定速度场模型参数，此时采用该参数的速度场模型进行全波形反演。
12.进一步，所述步骤b中的目标函数采用最小化位于检波点x
r
的观测数据d
obs
(x
r
,ω,t)和模拟数据u(x
r
,ω,t；m)的残差，具体为：
[0013][0014]
其中，x表示计算区域，m(x)表示计算区域x处的速度场模型参数；x
r
为检波点的位置信息；ω表示角频率，t表示时间；m表示速度场模型参数，u(x
r
,ω,t；m)表示在当前速度场模型参数m下获得的模拟数据。
[0015]
进一步，所述步骤a中的震源为同时源模拟波场：
[0016]
给定一系列位于x
s
的单频随机信号s(x
s
,ω)，则震源能表达为：
[0017][0018]
其中，x
s
为震源位置信息，n
s
表示一个震源内编码信号的个数；是表示震源编码信号中第i个信号的位置；ω＝(ω1,ω2,
…
,ω
n
)，ω
i
表示编码信号中第i个信号的角频率。
[0019]
进一步，所述步骤c中震源的正传波场由求解模拟方程得到，具体模拟方程为：
[0020][0021]
其中，x表示计算区域，x
s
震源位置信息，l[
·
]是声波正演模拟算子；
[0022]
震源对应的虚拟震源为：
[0023][0024]
根据确定的虚拟震源，获取伴随方程，所述伴随方程为：
[0025][0026]
其中，表示伴随算子，u
′
(x,ω,t；m)表示反传波场；
[0027]
求解伴随方程，获得虚拟震源的反传波场。
[0028]
进一步，所述步骤d中无串扰解耦具体为：
[0029]
采用相敏检测法先获取第一参考信号和第二参考信号，所述第一参考信号的振幅和频率均与所述第二参考信号的振幅和频率相同，所述第一参考信号的相位与所述第二参考信号的相位相差90
°
；然后根据第一参考信号和第二参考信号，采用三角函数正交性对震源的正传波场和虚拟震源的反传波场分别进行无串扰解耦。
[0030]
一种振幅、相位信息可调节的全波形反演装置，包括：
[0031]
目标函数构建模块，用于基于观测数据和模拟数据，构建表征所述观测数据与所述模拟数据之间数据残差的目标函数；
[0032]
正传波场获取模块，用于获取震源的正传波场；
[0033]
反传波场获取模块，用于获取虚拟震源的反传波场；
[0034]
梯度计算模块，用于根据正传波场、反传波场和目标函数关于速度场模型参数的相位反演、全波场反演梯度计算表达式，计算速度场模型参数的梯度；
[0035]
迭代更新模块，用于根据计算获得的梯度对当前速度场模型参数进行迭代更新，直至所述数据残差达到设定阈值或者迭代更新次数达到预设值，确定速度场模型参数。
[0036]
与现有技术相比，本发明采用一系列单频子波同时激发形成混叠波场；并采用了相敏检测法从混叠波场中解耦出每一个编码信号对应的波场变量，实现了无串扰的梯度计算，从而有利于获得准确的速度场模型参数的更新方向。另外在目标函数关于速度场模型参数的梯度计算表达式计算中，采用在不同的反演阶段使用不同的计算表达式，即在反演初期(即低频反演阶段)主要以相位信息反演为主，在反演的中后阶段，采用全波形反演，可以减弱反演对初始模型的依赖性。
[0037]
综上，本发明通过多个单频编码信号同时激发实现同时源模拟；采用一阶速度—应力声波动方程，在一次正演模拟中完成多个震源子波的波场模拟，通过相敏检测方法分别对正传和反传波场实施波场解耦计算，求取无串扰的梯度结果，通过使用相位反演和全波形反演相结合的策略，可以降低目标函数的复杂度，从而减少反演中的多解性，实现高精度的反演。
附图说明
[0038]
图1是本发明实施例的步骤示意图；
[0039]
图2是overthrust模型下真实速度场模型参数和初始速度场模型参数的示意图；
[0040]
其中，图2a为真实的速度模型，图2b为反演所用的初始速度模型；
[0041]
图3是本发明中采用相敏检测方法对反传源波场进行解耦得到的其对应的解耦波场的振幅和相位信息；
[0042]
其中，图3a表示解耦波场的振幅信息，图3b表示解耦波场的相位信息；
[0043]
图4是本发明实施例运用到overthrust模型得到的反演结果与常规方法得到的反演结果的对比示意图；
[0044]
其中，图4a是仅使用波场振幅信息的同时源反演结果；图4b是仅使用波场相位信息的同时源反演结果；图4c是同时使用波场振幅信息的同时源反演结果，即全波形反演方法；图4d是采用本发明实施例得到的结果；
[0045]
图5是本发明实施例运用到overthrust模型中得到的反演的目标函数下降曲线对比示意图。
具体实施方式
[0046]
下面将结合本文实施例中的附图，对本文实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本文一部分实施例，而不是全部的实施例。基于本文中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本文保护的范围。
[0047]
如图1所示，本发明实施例的反演方法，包括：
[0048]
s1：通过多个位于不同位置的单频子波同时激发作为震源，采用多个检波点组成观测系统，通过观测系统获得观测数据；建立已知的速度场模型(采用overthrust模型)，并根据震源情况计算得出模拟数据；
[0049]
s2：建立表征观测数据与模拟数据之间数据残差的目标函数；
[0050]
所述目标函数采用最小化位于检波点x
r
的观测数据d
obs
(x
r
,ω,t)和模拟数据u(x
r
,ω,t；m)的残差，具体为：
[0051][0052]
其中，x表示计算区域，m(x)表示计算区域x处的速度场模型参数；x
r
为检波点的位置信息；ω表示角频率，t表示时间；m表示速度场模型参数，u(x
r
,ω,t；m)表示在当前速度场模型参数m下获得的模拟数据。
[0053]
s3：根据震源与检波点之间的位置，得出虚拟震源的信息，然后根据观测数据获取震源的正传波场；并根据目标函数获得虚拟震源的反传波场；
[0054]
给定一系列位于x
s
的单频随机信号s(x
s
,ω)，则所述震源可表达为：
[0055][0056]
其中，x
s
为震源位置信息，n
s
表示一个震源内编码信号的个数；是表示震源编码信号中第i个信号的位置；ω＝(ω1,ω2,
…
,ω
n
)，ω
i
表示编码信号中第i个信号的角频率。
[0057]
则震源对应的正传波场模拟可以表示为：
[0058][0059]
其中，u(x,ω,t；m)表示正传波场，l[
·
]为声波正演模拟算子。
[0060]
相应地，伴随方程为：
[0061]
[0062]
式(4)中，表示伴随算子，u
′
(x,ω,t；m)表示反传波场；为虚拟震源，则所述虚拟震源可表达为：
[0063][0064]
式(5)中，式(5)中，表示编码信号中第i个信号对应检波点的位置。求解式(4)和式(5)，即可模拟得到震源的正传波场和虚拟震源的反传波场。
[0065]
通常情况下，由于观测数据和模拟数据的观测系统不一致，所以需要分别从观测和模拟数据对应的单个信号数据中提取相同频率的信号，在观测数据的观测系统约束下求取数据的残差，以此解决观测系统不一致的问题。
[0066]
s4：使用正传波场和反传波场进行梯度计算，迭代计算目标函数数据残差以更新速度场模型参数：
[0067]
具体地，包括如下步骤：
[0068]
s41：获得目标函数对速度场模型参数的梯度计算表达式；
[0069]
s42：根据同时源模拟进行正演和反传计算，获得解耦的波场变量；
[0070]
本实施例中，梯度计算表达式为正传波场与反传波场的特定形式的互相关，可简写为：
[0071][0072]
梯度简写形式g(x,m)表达式中，ψ
s
(x,ω)和ψ
r
(x,ω)分别表示正传波场与反传波场。从而，由(6)式的第二个等式可以看到，等式右端的第一项表示每一个编码信号对应波场数据的自相关运算；第二项表示当前编码信号与其它编码信号的波场互相关运算，即串扰。
[0073]
由(6)式可知，只有从混叠的波场中解耦出所有编码信号的波场，然后只计算当前编码信号对应的正反传波场的自相关，才能从根本上消除串扰噪声。
[0074]
即，有效地解耦混叠波场是解决串扰噪声的关键，本实施例引入相敏检测(psd)法来解决这个问题。psd的主要思想是：利用两个与未知信号具有相同频率但彼此相位相差90
°
的参考信号进行时间积分来提取振幅为sin(ω
i
t)、相位为cos(ω
i
t)的未知信号(即目标信号)。
[0075]
本实施例中，取e
ref
sin(ω
i
t+θ
ref
)为第一参考信号，取e
ref
cos(ω
i
t+θ
ref
)为第二参考信号；其中，e
ref
为第一参考信号和第二参考信号的振幅；θ
ref
为所述第一参考信号和第二参考信号的相位。由此可以看出，第一参考信号的振幅和频率均与第二参考信号的振幅和频率相同，第一参考信号的相位与第二参考信号的相位相差90
°
。
[0076]
由于空间某一点的混叠波场，均可以表示为：
[0077]
[0078]
其中，分别为每个单频编码信号的振幅，分别为每个单频编码信号的频率；以及分别为每个单频编码信号的相位。
[0079]
使|ω
i
‑
ω
j
|≥pδω(i,j∈n
s
,p∈z)，且δω＝2π/t。这里的t是地震记录的长度，p是最小的整数，pδω表示编码信号之间的最小频率间隔。当定义最小的积分区间为t
p
＝t/p，对于任何q(p∈[1,p])倍的t
p
区间qt
p
都作为合理的积分区间。通常，需要一定的时间等待波场达到稳态状况，所以p比q要小。
[0080]
在积分区间qt
p
上，将所述第一参考信号与混叠波场的信号相乘并积分，得到：
[0081][0082]
t
p
为最小积分区间；
[0083]
以及，将所述第二参考信号与所述混叠波场的信号相乘并积分，得到：
[0084][0085]
从而，与第一参考信号和第二参考信号频率相同的编码信号的对应波场便能解耦出来，其计算公式为：
[0086][0087]
上述积分过程主要利用了三角函数正交性，只有与第一参考信号和第二参考信号频率相同的信号其得到积分结果不为零，其他与第一参考信号和第二参考信号频率不同的信号其得到的积分结果为零，从而实现对空间中任意一混叠波场中各频率编码信号其对应波场的解耦。
[0088]
s43：利用目标函数对速度场模型参数的梯度计算表达式和解耦的波场变量，获取模型参数梯度；
[0089]
其中速度模型的梯度计算式为：
[0090][0091]
其中，re表示取向量的实数部分；u(x,ω
i
,t；m)为正传波场中，频率为ω
i
的编码信号解耦得到的其对应的波场变量，u
′
(x,ω
i
,t；m)为反传波场中，频率为ω
i
的编码信号经解耦得到其对应的波场变量。式(11)即是对频率均为ω
i
的两个波场变量进行正相关运算，从而实现无串扰的梯度求解。
[0092]
本实施例采用五个频带的多尺度反演，在前两个频带反演(即反演的初始阶段)，采用具有相位信息反演的梯度表达式，在后面三个频带反演(即反演的中后阶段)，采用兼顾波场相位和振幅信息的全波形反演的梯度表达式。
[0093]
基于上述梯度表达式，拓展出具有相位信息反演的梯度表达式为：
[0094]
[0095]
其中，表示取波场的相位信息，x表示计算区域，ω表示角频率，m表示模型参数，t表示时间，n
s
表示一个震源内编码信号的个数，u(x,ω
i
,t；m)表示正传波场，u
′
(x,ω
i
,t；m)表示反传波场；
[0096]
与之对应的兼顾波场相位和振幅信息的全波形反演的梯度表达式为
[0097][0098]
其中提取波场信息的振幅信息。
[0099]
s5：采用共轭梯度方法求取更新方向，对模型参数梯度进行迭代更新，直至残差达到设定值或者迭代更新次数达到预设值，停止迭代更新输出速度场模型参数，此时采用该参数的速度场模型进行全波形反演。
[0100]
采用共轭梯度方法求取更新方向δm
k+1
，
[0101][0102]
式(14)中的下标k表示迭代次数，s
k
是用于计算共轭方向的一个中间变量，上标t表示矩阵转置，更新后的速度场模型为：
[0103]
m
k+1
＝m
k
+t
k
δm
k+1
ꢀꢀ
(15)
[0104]
即将式(14)中计算得到的δm
k+1
带入到式(15)中，式(15)中的t
k
表示选取的更新步长，t
k
δm
k+1
即为此次迭代过程中速度场模型参数的更新量。
[0105]
上述反演方法所需的反演装置，包括：
[0106]
目标函数构建模块，用于基于观测数据和模拟数据，构建表征所述观测数据与所述模拟数据之间数据残差的目标函数；
[0107]
正传波场获取模块，用于获取震源的正传波场；
[0108]
反传波场获取模块，用于获取虚拟震源的反传波场；
[0109]
梯度计算模块，用于根据正传波场、反传波场和目标函数关于速度场模型参数的相位反演、全波场反演梯度计算表达式，计算速度场模型参数的梯度；
[0110]
迭代更新模块，用于根据计算获得的梯度对当前速度场模型参数进行迭代更新，直至所述数据残差达到设定阈值或者迭代更新次数达到预设值，确定速度场模型参数。
[0111]
效果验证：
[0112]
将本发明实施例运用到overthrust模型得到的反演结果与常规方法得到的反演结果进行对比；由图4b与图4a对比，能看出相位反演相对于振幅反演的精度更高，说明相位反演对初始模型的依赖性较小，可以恢复出有效的低波数信息；图4c与图4d相比，可以看出，采用本发明实施例全波形反演方法其结果比常规反演在界面的清晰度，高波数的连续性等方面更好。
[0113]
如图5所示，采用本发明实施例的全波形反演方法与常规方法对应的目标函数下降曲线对比，由图可知，相比于其他几种反演方法，振幅反演的目标函数在低频下降缓慢，
这说明振幅反演过早陷入局部极值；虽然同时使用相位和振幅的全波形反演在低频有着最大的下降率，但是随着迭代次数的增加，其目标函数的下降曲线逐渐与相位反演相一致；对于本发明实施例的方案，在开始阶段与相位反演吻合度很高，随着迭代次数的增加，其目标函数值逐渐低于相位反演，直至最后的频带、组合反演的目标函数值是四种反演算法中最低的，这说明本发明实施例的反演方法可以降低目标函数的复杂度，使得目标函数下降的更多。