一种用于常规陆域地震数据的多尺度全波形反演方法

1.本发明涉及地震勘探技术领域，具体涉及一种用于常规陆域地震数据的多尺度全波形反演方法。


背景技术：

2.目前，全波形反演已成功应用于海域地震数据中，但是在常规陆域地震数据的应用上还存在诸多难点。与海域地震数据相比，常规陆域地震数据偏移距小且缺少低频信息，由于常规工业检波器的频率为10hz，而陆域地震数据的频率往往高于5hz，所以常规陆域地震数据对于全波形反演方法而言频率太高，采用全波形反演方法处理常规陆域地震数据，无法得到较好的效果。
3.同时，全波形反演方法的应用受限于梯度类反演方法普遍存在的“周跳”现象。“周跳”现象是指反演时，当初始记录与实际记录之间差值小于1/2周期时，反演收敛得到全局最优值，而当初始记录与实际记录之间差值大于1/2周期时，反演收敛得到局部最优值，无法得到正确的结果。因此，亟需提出一种用于常规陆域地震数据的多尺度全波形反演方法，提高全波形反演结果的准确性。


技术实现要素：

4.针对现有技术中全波形反演方法受限于梯度类反演方法普遍存在的“周跳”难题，本发明提出了一种用于常规陆域地震数据的多尺度全波形反演方法，以传统频率域多尺度为基础，通过增加层剥离多尺度和偏移距多尺度，实现了“由近到远，由浅到深，由低频到高频”的多尺度反演，提高了全波形反演的精度，有利于常规陆域的地质解释。
5.为了实现上述目的，本发明采用如下技术方案：
6.一种用于常规陆域地震数据的多尺度全波形反演方法，具体包括如下步骤：
7.步骤1，输入原始地震数据，获取原始地震数据的初至旅行时，基于二维射线追踪地震层析成像方法，确定层析成像结果，对层析成像结果进行周跳风险评价，确定全波形反演的速度模型，根据全波形反演模型的速度模型，确定全波形反演的密度模型，建立全波形反演模型；
8.步骤2，对原始地震数据进行预处理，获得校正后的地震数据，利用全波形反演的速度模型进行正演模拟得到正演模拟地震数据，基于均方根振幅最小二乘拟合方法，对正演模拟地震数据与校正后的地震数据进行振幅匹配，并将校正后的地震数据由时间域变换至频率域，作为多尺度全波形反演数据；
9.步骤3，输入多尺度全波形反演数据和全波形反演模型，设置频率组、拉普拉斯常数组、偏移距组和迭代次数，基于层剥离的laplace
‑
fourier域多尺度全波形反演方法，通过设置由近及远遍历偏移距组内偏移距的内部循环、由小到大遍历拉普拉斯常数组内拉普拉斯常数的中间循环以及由低到高的遍历频率组内频率的外部循环，利用全波形反演模型进行全波形反演，获得多尺度全波形反演结果，根据多尺度全波形反演结果进行地质解释。
10.优选地，所述步骤1中，具体包括以下步骤：
11.步骤1.1，输入原始地震数据，拾取原始地震数据的初至旅行时获得实际旅行时，结合原始地震数据选取初至数据，设置层析成像参数；
12.步骤1.2，根据初至数据，基于二维射线追踪地震层析成像方法，得到层析成像结果，确定层析成像速度模型；
13.步骤1.3，对层析成像速度模型进行正演模拟，计算各炮点激发检波器位置处的初至旅行时，获得各炮点的正演旅行时，结合各炮点的实际旅行时，计算相同炮点和检波器位置处正演旅行时与实际旅行时之间的差值；
14.步骤1.4，设置周跳风险评价频率，计算评价周期，进行周跳风险评价，判断层析成像速度模型是否会产生周跳现象，若差值小于1/2个评价周期，则层析成像速度模型不会产生周跳现象，执行步骤1.6；若差值不小于1/2个评价周期，则层析成像速度模型会产生周跳现象，执行步骤1.5；
15.步骤1.5，选取层析成像速度模型中差值不小于1/2个评价周期的区域，增加炮点数据，继续进行正演模拟，返回步骤1.3；
16.步骤1.6，输出层析成像速度模型作为全波形反演模型中的速度模型，基于gardner经验关系，确定全波形反演模型中的密度模型，设置品质因子模型，建立全波形反演模型。
17.优选地，所述步骤1.3中，基于频率域2d粘弹声波方程进行正演模拟，如式(1)所示：
[0018][0019]
式中，ω为角频率；κ(x,z)为体积模量；p(x,z,ω)为频率域地震波场；ρ(x,z)为密度；s(x,z,ω)为震源项。
[0020]
优选地，所述步骤1.6中，基于gardner经验关系，利用全波形反演模型的速度模型确定全波形反演模型中的密度模型，计算公式为：
[0021][0022]
式中，ρ为密度，单位为g/cm3；v
p
为纵波速度，单位为m/s。
[0023]
优选地，所述步骤2中，具体包括以下步骤：
[0024]
步骤2.1，原始地震数据预处理，预处理包括静校正、去噪、去面波、振幅补偿、一次性校正和初至切除，得到校正后的地震数据；
[0025]
步骤2.2，利用全波形反演的速度模型对校正后的地震数据进行正演模拟，得到正演模拟地震数据，基于均方根振幅最小二乘拟合方法，将正演模拟地震数据的振幅与校正后地震数据的振幅逐炮匹配，使得正演模拟地震数据与校正后地震数据的尺度相同；
[0026]
正演模拟地震数据与校正后地震数据的振幅匹配公式为：
[0027]
y＝cx
ꢀꢀ
(3)
[0028]
式中，x为正演模拟地震数据的均方根振幅；y为校正后地震数据的均方根振幅；c为缩放系数；
[0029]
步骤2.3，将校正后的地震数据由时间域变换至频率域，得到多尺度全波形反演数
据；
[0030]
校正后地震数据的变换公式如下：
[0031][0032]
式中，ω为角频率；t为时间；f(t)为校正后地震数据的时间域信号；f(ω)为校正后地震数据的频率域信号。
[0033]
优选地，所述步骤3中，具体包括以下步骤：
[0034]
步骤3.1，输入多尺度全波形反演数据和全波形反演模型，分别设置频率组、拉普拉斯常数组、偏移距组和迭代次数，其中，频率组内按照频率由低到高的顺序设置有a个频率，拉普拉斯常数组内按照数值由小到大的顺序设置有b个拉普拉斯常数，偏移距组内按照数值由小到大的顺序设置有c个偏移距，分别对频率组、拉普拉斯常数组和偏移距组内各参数进行编号；
[0035]
步骤3.2，分别设置频率、拉普拉斯常数和偏移距的初始值，其中，将频率组内的最低频率设置为频率的初始值，将拉普拉斯常数设置为拉普拉斯常数组内的最小值，将偏移距设置为偏移距组内的最小值；
[0036]
步骤3.3，设置梯度吸收因子ε和高斯平滑因子ζ；
[0037]
步骤3.4，根据当前设置的偏移距和频率从多尺度全波形反演数据中提取地震数据，利用全波形反演模型进行全波形反演，得到全波形反演结果，并将本次迭代计算得到的全波形反演结果设置为全波形反演模型；
[0038]
步骤3.5，判断是否完成设置的迭代次数；
[0039]
若未完成设置的迭代次数，则返回步骤3.4继续进行全波形反演；若已完成设置的迭代次数，则执行步骤3.6；
[0040]
步骤3.6，判断是否遍历偏移距组内全部的偏移距；
[0041]
若未遍历完，则根据偏移距组对当前偏移距o
i
进行更新，i＝1，
…
，c
‑
1，i为当前偏移距在偏移距组内的序号，更新后的偏移距为o
i+1
，返回步骤3.4继续进行全波形反演；若已遍历完，则对偏移距进行初始化，设置偏移距为o1，重新遍历偏移距组内的偏移距，执行步骤3.7；
[0042]
步骤3.7，判断是否遍历完拉普拉斯常数组内全部的拉普拉斯常数；
[0043]
若未遍历完，则根据拉普拉斯常数组对当前拉普拉斯常数l
j
进行更新，j＝1，
…
，b
‑
1，j为当前拉普拉斯常数在拉普拉斯常数组内的序号，更新后的拉普拉斯常数为l
j+1
，返回步骤3.4继续进行全波形反演；若已遍历完，则对拉普拉斯常数进行初始化，设置拉普拉斯常数为l1，重新遍历拉普拉斯常数组内的拉普拉斯常数，执行步骤3.8；
[0044]
步骤3.8，判断是否遍历完频率组内全部的频率；
[0045]
若未遍历完，则根据频率组对当前频率f
n
进行更新，n＝1，
…
，a
‑
1，n为当前频率在频率组内的序号，更新后的频率为f
n+1
，返回步骤3.4继续进行全波形反演；若已遍历完，则将当前全波形反演模型作为反演速度模型，执行步骤3.9；
[0046]
步骤3.9，输出多尺度全波形反演结果，并根据多尺度全波形反演结果进行地质解释。
[0047]
优选地，所述步骤3.4中，反演过程中全波形反演模型的参数迭代更新公式为：
[0048][0049]
式中，d
obs
为实测地震数据，d
cal
为正演模拟地震数据；c(m)为加权最小二乘目标函数，为转置操作，w
d
为数据预处理算子；m
k+1
为第k+1次迭代得到的属性模型，m
k
为第k次迭代得到的属性模型，δm
k
为第k次模型修正量；为目标函数的最速下降方向，α
k
为沿最速下降方向的步长，ε为梯度吸收因子，i为单位矩阵，ζ为高斯平滑因子，为高斯
‑
牛顿海森矩阵的近似对角矩阵。
[0050]
本发明所带来的有益技术效果：
[0051]
本发明为了克服现全波形反演方法的周跳难题，利用初始模型控制初始数据与实测数据之间的差值，充分利用实测数据中的低频信息，基于层剥离方法在laplace
‑
fourier域进行全波形反演，通过依次优化小偏移距浅层地震数据、大偏移距浅层地震数据、小偏移距深层地震数据和大偏移距深层地震数据，采用“由近到远，由浅到深，由低频到高频”的多尺度反演策略，有效提高了速度反演的精度，为常规陆域的地质解释提供了依据。
附图说明
[0052]
图1为一种用于常规陆域地震数据的多尺度全波形反演方法流程图。
[0053]
图2为全波形反演模型的构建流程图。
[0054]
图3为基于层剥离的laplace
‑
fourier多尺度全波形反演方法流程图。
[0055]
图4为反演速度与地温曲线对比图；图(a)为地温曲线，图(b)为反演速度曲线。
具体实施方式
[0056]
下面结合附图以及具体实施方式对本发明作进一步详细说明：
[0057]
实施例1
[0058]
本发明提出了一种用于常规陆域地震数据的多尺度全波形反演方法，如图1所示，具体包括如下步骤：
[0059]
步骤1，建立全波形反演模型，如图2所示，具体包括以下步骤：
[0060]
步骤1.1，输入原始地震数据，拾取原始地震数据的初至旅行时获得实际旅行时，结合原始地震数据选取初至数据，设置层析成像参数。
[0061]
步骤1.2，根据初至数据，基于二维射线追踪地震层析成像方法，得到层析成像结果，确定层析成像速度模型。
[0062]
步骤1.3，对层析成像速度模型进行正演模拟，计算各炮点激发检波器位置处的初至旅行时，获得各炮点的正演旅行时，结合各炮点的实际旅行时，计算相同炮点和检波器位置处正演旅行时与实际旅行时之间的差值。
[0063]
步骤1.4，设置周跳风险评价频率，计算评价周期，进行周跳风险评价，判断层析成
像速度模型是否会产生周跳现象，若差值小于1/2个评价周期，则层析成像速度模型不会产生周跳现象，执行步骤1.6；若差值不小于1/2个评价周期，则层析成像速度模型会产生周跳现象，执行步骤1.5。
[0064]
步骤1.5，选取层析成像速度模型中差值不小于1/2个评价周期的区域，增加炮点数据，继续进行正演模拟，返回步骤1.3。
[0065]
步骤1.6，输出层析成像速度模型作为全波形反演模型中的速度模型，基于gardner经验关系，确定全波形反演模型中的密度模型，设置品质因子模型，建立全波形反演模型。
[0066]
步骤2，通过对原始地震数据进行预处理，获取多尺度全波形反演数据，具体包括以下步骤：
[0067]
步骤2.1，原始地震数据预处理，预处理包括静校正、去噪、去面波、振幅补偿、一次性校正和初至切除，得到校正后的地震数据。
[0068]
步骤2.2，利用全波形反演的速度模型对校正后的地震数据进行正演模拟，得到正演模拟地震数据，基于均方根振幅最小二乘拟合方法，将正演模拟地震数据的振幅与校正后地震数据的振幅逐炮匹配，使得正演模拟地震数据与校正后地震数据的尺度相同。
[0069]
步骤2.3，将校正后的地震数据由时间域变换至频率域，获得多尺度全波形反演数据。
[0070]
步骤3，多尺度全波形反演，如图3所示，具体包括以下步骤：
[0071]
步骤3.1，输入多尺度全波形反演数据和全波形反演模型，分别设置频率组f＝[f1，
…
，f
a
]、拉普拉斯常数组l＝[l1，
…
，l
b
]、偏移距组o＝[o1，
…
，o
c
]和迭代次数，其中，频率组内的频率按照由低频到高频的顺序排列，拉普拉斯常数组内的拉普拉斯常数按照由小到大的顺序排列，偏移距组内的偏移距按照由小到大的顺序排列，分别对频率组、拉普拉斯常数组和偏移距组内的参数进行编号。
[0072]
步骤3.2，分别设置反演频率、拉普拉斯常数和偏移距的初始值，其中，将频率组内的最低频率设置为频率的初始值，将拉普拉斯常数设置为拉普拉斯常数组内的最小值，将偏移距设置为偏移距组内的最小值，即设置频率的初始值为f1、拉普拉斯常数的初始值为l1、偏移距的初始值为o1。
[0073]
步骤3.3，设置梯度吸收因子ε和高斯平滑因子ζ。
[0074]
步骤3.4，根据当前设置的偏移距和频率从多尺度全波形反演数据中提取地震数据，利用全波形反演模型进行全波形反演，得到全波形反演结果，并将本次迭代计算得到的全波形反演结果设置为全波形反演模型。
[0075]
步骤3.5，判断是否完成设置的迭代次数。
[0076]
若未完成设置的迭代次数，则返回步骤3.4继续进行全波形反演；若已完成设置的迭代次数，则执行步骤3.6。
[0077]
步骤3.6，判断是否遍历偏移距组内全部的偏移距。
[0078]
若未遍历完，则根据偏移距组对当前偏移距o
i
进行更新，i＝1，
…
，c
‑
1，i为当前偏移距在偏移距组内的序号，更新后的偏移距为o
i+1
，返回步骤3.4继续进行全波形反演；若已遍历完，则对偏移距进行初始化，设置偏移距为o1，重新遍历偏移距组内的偏移距，执行步骤3.7。
[0079]
步骤3.7，判断是否遍历完拉普拉斯常数组内全部的拉普拉斯常数。
[0080]
若未遍历完，则根据拉普拉斯常数组对当前拉普拉斯常数l
j
进行更新，j＝1，
…
，b
‑
1，j为当前拉普拉斯常数在拉普拉斯常数组内的序号，更新后的拉普拉斯常数为l
j+1
，返回步骤3.4继续进行全波形反演；若已遍历完，则对拉普拉斯常数进行初始化，设置拉普拉斯常数为l1，重新遍历拉普拉斯常数组内的拉普拉斯常数，执行步骤3.8。
[0081]
步骤3.8，判断是否遍历完频率组内全部的频率。
[0082]
若未遍历完，则根据频率组对当前频率f
n
进行更新，n＝1，
…
，a
‑
1，k为当前频率在频率组内的序号，更新后的频率为f
n+1
，返回步骤3.4继续进行全波形反演；若已遍历完，则将当前全波形反演模型作为反演速度模型，执行步骤3.9。
[0083]
步骤3.9，输出多尺度全波形反演结果，并根据多尺度全波形反演结果进行地质解释。
[0084]
实施例2
[0085]
将本发明提出的一种用于常规陆域地震数据的多尺度全波形反演方法应用于理想模型中进行测试，取得了理想的反演效果。
[0086]
根据某地区的地震剖面和岩石物性资料，构建二维纵波速度模型，基于gardner经验关系，将纵波速度转化为密度，构建密度模型，设置品质因子模型为常数，进行多尺度全波形反演，多尺度全波形反演过程中所采用的参数如表1所示。
[0087]
多尺度全波形反演过程中，根据该地区实际地震数据的观测系统，设置正演地震数据的观测系统最大偏移距为14.4km，正演地震数据观测系统的最大偏移距与实际地震数据观测系统相同。设置偏移距组内偏移距分别为5km和14.4km，设置laplace常数组内laplace常数分别为0.5和1.0；梯度吸收因子分为两组，分别为1e
‑
2和1e
‑
6；高斯平滑因子分为两组，水平方向分别为1.2km和0.6km，垂直方向分别为0.6km和0.3km；频率采用渐进式由低频到高频多尺度策略，从2hz到20hz分为16个频组。对实际模型进行高斯平滑，作为反演的初始模型。
[0088]
表1理论模型多尺度全波形反演所采用的参数
[0089][0090]
由反演结果可得，本发明方法可以很好的反演出模型的构造和速度值，在水平方向速度变化剧烈的地方，反演效果也较理想。通过抽取模型上不同位置处的一维速度进行对比，可以看出，在深度2.5km以上地层中，反演结果能够很好的与实际数据吻合，反演速度相对于初始模型的速度有了很大的提高。
[0091]
实施例3
[0092]
将本发明提出的一种用于常规陆域地震数据的多尺度全波形反演方法应用于实际数据处理中，取得了理想的效果。
[0093]
选取某地区陆域深反射地震数据40km，据测线900米位置处有一口地热井，井深2500米，作为反演结果的验证井。该地震数据激发震源为炸药震源，检波器为常规工业生产中使用的10hz模拟检波器，道距20m，炮距100m，接收道数1440道，最小偏移距10m，最大偏移距14.4km。
[0094]
对实际数据进行多尺度全波形反演，具体包括以下步骤：
[0095]
步骤1，建立初始模型；
[0096]
输入原始地震数据，获取原始地震数据的初至旅行时，基于二维射线追踪地震层析成像方法，确定层析成像结果，对层析成像结果进行周跳风险评价，确定全波形反演的速度模型，根据全波形反演模型的速度模型，确定全波形反演的密度模型，建立全波形反演模型。
[0097]
步骤2，地震数据预处理；
[0098]
首先，对原始地震进行常规地震数据处理，包括：静校正、去噪、球面振幅补偿、一致性校正、陷波、面波去除、初至以上切除，获得校正后的地震数据；再利用全波形反演的速度模型进行正演模拟得到正演模拟地震数据；由于正演模拟地震数据与校正后的地震数据两者的振幅不在同一尺度内，因此，需要对正演模拟地震数据与校正后的地震数据进行振幅匹配，匹配完成后，为了为laplace
‑
fourier域多尺度全波形反演做准备，还需要将校正后的地震数据由时间域变换至频率域，获得多尺度全波形反演数据。
[0099]
步骤3，多尺度全波形反演；
[0100]
根据偏移距和频率从多尺度全波形反演数据内抽取数据，作为多尺度全波形反演的输入数据，进行多尺度全波形反演，多尺度全波形反演过程中所使用的参数如表2所示。
[0101]
表2实际数据多尺度全波形反演所采用的参数
[0102][0103]
通过分别对三组不同高斯平滑因子进行测试获得反演结果，得到采用数值大的平滑参数进行平滑，反演结果更平滑但是局部缺少细节，采用数值小的平滑参数进行平滑，反演结果细节较多但是深层也加入了更多扰动，通过合理设置平滑参数，可以得到高分辨率的反演结果。
[0104]
利用正演地震记录对比、叠前深度偏移分析和测井数据验证对多尺度全波形反演结果进行评价。
[0105]
正演地震记录对比：通过正演炮点位置x＝18km和20km处的地震记录，得到相比于初始模型地震记录，反演结果的地震记录在中浅层加入了丰富的反射波场信息，初至旅行时和中浅层反射信息与实际数据较一致。
[0106]
叠前深度偏移分析：分别利用初始速度模型和全波形反演速度模型进行地震炮集的叠前深度偏移，对比后得到，反演速度模型的叠前深度偏移结果在深部成像效果明显好于初始速度模型的成像效果，并且根据偏移道集上可以得到，利用反演模型偏移的道集同相轴较利用初始模型偏移的道集更平，显示了反演的速度更准确，更接近于实际地层速度。
[0107]
测井数据验证：利用d13地热井的测井数据验证反演模型，通过对比发现反演结果
与测井曲线吻合的很好，尤其是在1100
‑
1200m位置处，即测井解释的地热储层位置，同时，该位置对应电阻率曲线上的低阻异常位置。再对比d13井地温曲线与反演的速度曲线，得到在区域基岩界面以上(约900米)，地温梯度为3.4℃/100m，符合该区域正常地温梯度，而在地热储层位置处(1100
‑
1200米)，地温梯度降到约1.0℃/100m，对应于速度的低阻异常区，符合该区地热储层物性规律。
[0108]
因此，通过结合正演地震记录对比、叠前深度偏移分析和测井数据验证验证了本发明方法的应用效果，验证了本技术方法充分利用地震数据里的低频信息，有利于常规陆域地震数据的处理，为常规陆域的地质解释提供了依据。
[0109]
当然，上述说明并非是对本发明的限制，本发明也并不仅限于上述举例，本技术领域的技术人员在本发明的实质范围内所做出的变化、改型、添加或替换，也应属于本发明的保护范围。