反演成像速度与品质因子粘声最小二乘逆时偏移成像方法

1.本发明涉及石油勘探技术领域，尤其涉及一种反演成像速度与品质因子粘声最小二乘逆时偏移成像方法。


背景技术：

2.地下介质具有严重的粘弹性，对地震波有很强的吸收和衰减作用。q衰减导致反射层振幅衰减和成像错位。声波介质下的最小二乘逆时偏移方法只能反演反射系数模型，常规粘声最小二乘逆时偏移方法大多基于常q品质因子推导而来，对于q品质因子的初始精度要求较高，但现实情况下无法得到精度较高的q品质因子。


技术实现要素：

3.为解决上述技术问题，本发明公开了一种反演成像速度与品质因子粘声最小二乘逆时偏移成像方法，设计合理，克服了现有技术的不足，具有良好的效果。
4.为实现上述目的，本发明采用下述技术方案：
5.一种反演成像速度与品质因子粘声最小二乘逆时偏移成像方法，包括如下步骤：
6.s1、输入速度场及品质因子参数场，建立观测系统；
7.s2、推导q衰减线性born算子；
8.s3、推导q补偿伴随算子；
9.s4、推导q补偿梯度算子；
10.s5、共轭梯度法计算步长；
11.s6、更新成像速度与品质因子。
12.可选地，步骤s2中，基于二阶粘声波动方程，如下式：
[0013][0014]
其中，p表示应力，x和z分别表示空间水平分量和垂直分量的空间坐标，t表示时间，表示拉普拉斯算子，v表示速度，σ表示正则化算子，τ表示无量纲变量，由下式计算得到：
[0015]
τ＝τ
ε
/τ
σ-1
ꢀꢀꢀ
(2)
[0016][0017][0018]
其中，τ
ε
表示应变松弛时间，τ
σ
表示应力松弛时间，q表示品质因子；
[0019]
在式(1)中，ε表示q衰减和q补偿转换参数，由下式给出：
[0020][0021]
在式(1)中，如果参考速度v0被v代替，粘声方程简化为有损方程；
[0022]
v0和v之间的关系由下式表达：
[0023][0024]
其中，ω0表示参考角频率；
[0025]
基于式(6)的maclaurin级数展开，得到：
[0026][0027]
考虑到：
[0028][0029]
式(7)简化为：
[0030][0031]
考虑v和τ(δv,δτ)的扰动，v和τ的反偏移方程由下式得到：
[0032][0033]
其中，vs和τs分别表示v和τ的背景场，表示的背景场；
[0034]
将ε＝-1带入式(1)，计算得到
[0035][0036]
其中，f表示震源函数；
[0037]
在式(10)中，m(v)和m(τ)分别表示v和τ的扰动，由下式表达：
[0038][0039]
式(10)、(11)、(12)共同构成了q衰减线性born算子。
[0040]
可选地，步骤s2中，基于born近似，参数扰动δm(v,τ)＝m-ms导致扰动波场δp；
[0041]
ps和p＝ps δp由下式计算得到：
[0042][0043]
[0044]
实施泰勒展开后，得到：
[0045][0046]
其中，o(δv)表示v的高阶项。
[0047]
可选地，步骤s3中，令式(1)中ε＝1，应用伴随状态理论，得到q补偿伴随方程：
[0048][0049]
其中，表示的伴随形式。
[0050]
可选地，步骤s4中，q补偿形式的波动方程如下所示：
[0051][0052]
其中，表示q补偿正演算子，表示q补偿伴随波场，表示重构的震源项，由下式给出：
[0053][0054][0055]
其中，由下式给出：
[0056][0057][0058]
基于伴随状态理论：
[0059][0060]
其中，表示粘声方程中的逆时延拓算子，表示的伴随算子，r
*
表示将炮记录扩展到整个空间的运算符，表示q衰减伴随波场；
[0061]
式(22)中，q衰减合成记录通过线性born正演算子计算得到：
[0062][0063]
为了补偿q衰减，将q补偿梯度方程从式(22)更改为：
[0064][0065]
因此，关于v和τ的梯度方向从式(24)得出：
[0066]
[0067][0068]
可选地，步骤s5中，第k次迭代的共轭梯度方程由下式得到：
[0069][0070]
其中，c表示预处理算子，β
(k)
表示最速下降法的第k次迭代的步长，由下式计算得到：
[0071][0072]
可选地，步骤s6中，共轭梯度法通过下式来更新模型：
[0073][0074]
其中，α表示共轭梯度法的步长，通过下式获得：
[0075][0076]
本发明的有益效果是，
[0077]
本发明提出了一种同时反演成像速度与品质因子的粘声最小二乘逆时偏移成像方法，解决了以往针对地层的衰减作用提出的粘介质偏移成像方法无法同时反演成像速度与品质因子的问题。该方法所需的初始速度模型与q品质因子模型精度较低，反演得到的成像结果精度高，信噪比高，振幅更加均衡，震源效应更弱。
附图说明
[0078]
图1为本发明一实施例中带有气块的三层模型；
[0079]
图2为本发明一实施例中的输入模型，其中，(a)为v，(b)为τ；
[0080]
图3为本发明一实施例中的模型扰动参数，其中，(a)为v，(b)为τ；
[0081]
图4为本发明一实施例中的600ms时刻的正演波场快照；
[0082]
图5为本发明一实施例中的正演模拟炮记录，其中，(a)为粘声介质，(b)为声介质，(c)为取(a)左半部分，取(b)右半部分；
[0083]
图6为本发明一实施例中的born正演模拟炮记录，(a)为粘声介质；(b)为声介质；(c)为取(a)左半部分，取(b)右半部分；
[0084]
图7为本发明一实施例中的粘声介质下的反传波场炮记录，(a)为补偿，(b)为未补偿；
[0085]
图8为本发明一实施例中的第1次迭代的补偿偏移结果，(a)为v，(b)为τ；
[0086]
图9为本发明一实施例中的第50次迭代的补偿偏移结果，(a)为v，(b)为τ；
[0087]
图10为本发明一实施例中的第50次迭代的未补偿偏移结果，(a)为v，(b)为τ；
[0088]
图11为本发明一实施例中的补偿lsrtm与未补偿lsrtm的残差收敛曲线；
[0089]
图12为本发明一实施例中的实际工区模型，(a)为v，(b)为τ；
[0090]
图13为本发明一实施例中的输入炮记录，(a)为粘声介质，(b)为声介质；
[0091]
图14为本发明一实施例中的输入模型，(a)为v，(b)为τ；
[0092]
图15为本发明一实施例中的拉普拉斯滤波后第1次迭代的补偿偏移结果，(a)为v，(b)为τ；
[0093]
图16为本发明一实施例中的拉普拉斯滤波后第50次迭代的补偿偏移结果，(a)为v，(b)为τ；
[0094]
图17为本发明一实施例中的反偏移炮记录，(a)为第1次迭代，(b)为第50次迭代，(c)为输入；
[0095]
图18为本发明一实施例中的拉普拉斯滤波后未补偿偏移结果，(a)为第1次迭代，(b)为第50次迭代；
[0096]
图19为本发明一实施例中的1200m处未补偿lsrtm与本方法对比图，(a)为波形图，(b)为频谱图。
具体实施方式
[0097]
下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。
[0098]
一种反演成像速度与品质因子粘声最小二乘逆时偏移成像方法，包括如下步骤：
[0099]
s1、输入速度场及品质因子参数场，建立观测系统；
[0100]
s2、推导q衰减线性born算子；
[0101]
可选地，步骤s2中，基于二阶粘声波动方程，如下式：
[0102][0103]
其中，p表示应力，x和z分别表示空间水平分量和垂直分量的空间坐标，t表示时间，表示拉普拉斯算子，v表示速度，σ表示正则化算子，τ表示无量纲变量，由下式计算得到：
[0104]
τ＝τ
ε
/τ
σ-1
ꢀꢀꢀ
(2)
[0105][0106][0107]
其中，τ
ε
表示应变松弛时间，τ
σ
表示应力松弛时间，q表示品质因子；
[0108]
在式(1)中，ε表示q衰减和q补偿转换参数，由下式给出：
[0109][0110]
在式(1)中，如果参考速度v0被v代替，粘声方程简化为有损方程；
[0111]
v0和v之间的关系由下式表达：
[0112]
[0113]
其中，ω0表示参考角频率；
[0114]
基于式(6)的maclaurin级数展开，得到：
[0115][0116]
考虑到：
[0117][0118]
式(7)简化为：
[0119][0120]
考虑v和τ(δv,δτ)的扰动，v和τ的反偏移方程由下式得到：
[0121][0122]
其中，vs和τs分别表示v和τ的背景场，表示的背景场；
[0123]
将ε＝-1带入式(1)，计算得到
[0124][0125]
其中，f表示震源函数；
[0126]
在式(10)中，m(v)和m(τ)分别表示v和τ的扰动，由下式表达：
[0127][0128]
式(10)、(11)、(12)共同构成了q衰减线性born算子。
[0129]
可选地，步骤s2中，基于born近似，参数扰动δm(v,τ)＝m-ms导致扰动波场δp；
[0130]
ps和p＝ps δp由下式计算得到：
[0131][0132][0133]
实施泰勒展开后，得到：
[0134][0135]
其中，o(δv)表示v的高阶项。
[0136]
s3、推导q补偿伴随算子；
[0137]
可选地，步骤s3中，令式(1)中ε＝1，应用伴随状态理论，得到q补偿伴随方程：
[0138][0139]
其中，表示的伴随形式。
[0140]
s4、推导q补偿梯度算子；
[0141]
可选地，步骤s4中，q补偿形式的波动方程如下所示：
[0142][0143]
其中，表示q补偿正演算子，表示q补偿伴随波场，表示重构的震源项，由下式给出：
[0144][0145][0146]
其中，由下式给出：
[0147][0148][0149]
基于伴随状态理论：
[0150][0151]
其中，表示粘声方程中的逆时延拓算子，表示的伴随算子，r
*
表示将炮记录扩展到整个空间的运算符，表示q衰减伴随波场；
[0152]
式(22)中，q衰减合成记录通过线性born正演算子计算得到：
[0153][0154]
为了补偿q衰减，将q补偿梯度方程从式(22)更改为：
[0155][0156]
因此，关于v和τ的梯度方向从式(24)得出：
[0157][0158][0159]
s5、共轭梯度法计算步长；
[0160]
可选地，步骤s5中，第k次迭代的共轭梯度方程由下式得到：
[0161][0162]
其中，c表示预处理算子，β
(k)
表示最速下降法的第k次迭代的步长，由下式计算得到：
[0163][0164]
s6、更新成像速度与品质因子。
[0165]
可选地，步骤s6中，共轭梯度法通过下式来更新模型：
[0166][0167]
其中，α表示共轭梯度法的步长，通过下式获得：
[0168][0169]
本发明一种同时反演成像速度与品质因子的粘声最小二乘逆时偏移成像方法，应用于粘声介质模型，取得了理想的计算效果。图1展示了带有气块的三层模型；图2展示了输入模型，其中，图2(a)为速度参数，图2(b)为q品质因子参数；图3展示了模型扰动参数，其中，图3(a)为速度参数，图3(b)为q品质因子参数；图4展示了600ms时刻的正演波场快照，图4(a)为未补偿情况下的，图4(b)为补偿情况下的；图5展示了正演单炮炮记录，图5(a)为未补偿情况下的，图5(b)为补偿情况下的，图5(c)展示了图5(a)取左半部分，图5(b)取右半边部分；图6展示了born正演单炮炮记录，其中，图6(a)为未补偿情况下的，图6(b)为补偿情况下的，图6(c)展示了图6(a)取左半部分，图6(b)取右半边部分；图7展示了粘声介质下的补偿的反传波场炮记录，其中，图7(a)为补偿情况下的，图7(b)为未补偿情况下的；图8展示了第1次迭代的补偿偏移结果，其中，图8(a)为速度参数，图8(b)为q品质因子参数；图9展示了第50次迭代的补偿偏移结果，其中，图9(a)为速度参数，图9(b)为q品质因子参数；图10展示了第50次迭代的未补偿偏移结果，其中，图10(a)为速度参数，图10(b)为q品质因子参数；图11展示了补偿lsrtm与未补偿lsrtm的残差收敛曲线。
[0170]
图12所示的实际工区模型，其中，图12(a)为速度参数，图12(b)为q品质因子参数。图13所示的输入炮记录，其中，图13(a)为粘声介质，图13(b)为声波介质。图14所示的输入模型参数，其中，图14(a)为速度参数，图14(b)为q品质因子参数。图16(a)为拉普拉斯滤波后第50次迭代补偿的偏移结果，图18(b)为拉普拉斯滤波后第50次迭代未补偿的偏移结果。通过对比可以发现本发明的偏移结果比未补偿的偏移结果具有更高的分辨率、更高的信噪比、更均衡的成像幅度和更弱的震源效应。
[0171]
当然，上述说明并非是对本发明的限制，本发明也并不仅限于上述举例，本技术领域的技术人员在本发明的实质范围内所做出的变化、改型、添加或替换，也应属于本发明的保护范围。