一种成像域井震联合多尺度层析反演方法与流程

1.本发明涉及油气勘探地震资料处理技术领域，特别涉及一种成像域井震联合多尺度层析反演方法。


背景技术：

2.目前，生产中常用的速度建模方法是基于射线理论的剩余曲率速度分析，主要应用反射波的走时信息，反演的结果往往是以低波数段为主的背景速度场，难以满足高精度地震勘探的需求。如何提高速度模型的精度，实现复杂地区的精确成像是一个急需解决的关键性问题，也是地震成像领域乃至地震勘探领域的一大难点。
3.基于波动方程理论的速度建模方法在理论上可以得到包含高频成分的速度场，分辨率高，能更好地适应速度变化剧烈的地区，但存在许多尚未解决的问题。波动方程偏移速度分析理论和实际应用尚未完善，初始模型问题和对速度模型的敏感度是该方法面临的一大问题，并且计算量庞大，处理分析不灵活。全波形反演虽然在理论上已基本完善，但其计算量太大，严重依赖低频信息和大偏移距地震数据。而对于陆上地震勘探，如何提取地震子波以及模拟复杂波场也是摆在全波形反演面前的难题。复杂地区的地震资料还存在信噪比低的问题，因此，应用全波形反演解决复杂地区的速度建模问题是短期内无法实现的。
4.另一种提高速度模型精度的方式是井震联合。井震联合在反演方面研究和应用十分广泛，但主要是反演储层参数、提高成像分辨率等。在层析反演方面的应用，主要是应用井位处的层位信息建立各向异性模型。将测井速度应用于偏移速度建模，应考虑测井速度与偏移速度的差异。测井速度信息只在有限位置处存在，如何约束整个三维偏移速度体，是一个需要解决的问题。现有技术中提出了一种基于强模型的反演方法，称为传播4d，其主要目的是将井内信息传播到数据集中，它不是作为4d反演和解释的第一步而发展起来的，而是一种整合来自油井的稳健先验信息以获得更详细和更高频解的方法。该方法并未通过引入先验地统计学关系来强制连续性，而是让数据驱动连续性：首先使用井反演所在位置的记录道；然后选择其邻域，并以已有结果为约束进行反演，依次类推，直到全部区域反演完成。


技术实现要素：

5.本发明的目的是针对常规偏移速度分析结果分辨率低的问题，提供一种成像域井震联合多尺度层析反演方法，以反演精细的地下速度场，为深度域地震成像提供技术支持。
6.本发明提供的一种成像域井震联合多尺度层析反演方法，具体步骤如下：
7.步骤一，基于先验速度，进行叠前深度偏移，获取深度偏移剖面和偏移距共成像点道集；
8.步骤二，在所述偏移剖面中做局部倾斜叠加，得x向倾角场和y向倾角场；
9.步骤三，对各所述成像点道集进行γ扫描，获得γ场；
10.步骤四，基于所述倾角场和所述γ场进行综合分析和质量监控，筛选出所述深度
偏移剖面中的反射点；
11.步骤五，以各所述反射点作为出射点，进行射线追踪到地表，获得射线路径并基于射线终点之间的偏移距计算剩余深度差；
12.步骤六，以各井为中心，划分反演区域；
13.步骤七，确定所述反演区域，基于射线路径及剩余深度差以及测井速度信息，构建局部区域的方程组，求解所述方程组，更新局部速度。
14.步骤八，判断所述区域是否全部完成反演，若是则更新走时残差，否则迭代所述步骤七。
15.步骤九，确定反演网格尺度，重整射线路径，并结合平滑约束以及测井约束，构建反演方程组并求解，获取模型更新量。
16.步骤十，判断是否更新所述尺度，若是，更新走时残差并迭代计算所述步骤九，若否，求和各所述尺度更新量，更新模型并输出。
17.进一步的，在所述步骤六，以各井为中心，划分反演区域以各井为中心，具体包括：以给定的宽度划分整个工区区域，以距离井最近的区域作为首先反演的区域，以距离井次近的区域作为下一次反演的区域，以此类推，由近及远地划分整个工区范围。
18.进一步的，在步骤七中，求解所述方程组采用的目标函数o的具体形式如下：
19.o＝cd||δt-gδs||2 ε1||s δs-s
well
||2 ε2||tr(s δs)||2ꢀꢀꢀ
(1)
20.其中，g为由射线路径构成的敏感度矩阵，δs为慢度更新量，δt为走时残差，cd为数据协方差矩阵，s为模型慢度，s
well
为测井慢度，t为构造方向旋转矩阵，r为正则化算子，ε1、ε2为权重系数，上式中，单项走时残差的计算公式为：
21.δt＝2sc·
δz
·
cosβ
·
cosα
ꢀꢀꢀ
(2)
22.sc为反射点慢度，β为出射角，α为地层倾角，δz为剩余深度差；
23.所述目标函数对应的所述方程组如下：
[0024][0025]
进一步的，在所述步骤七中，进行的局部速度反演为一维反演，反演前将各个网格的射线长度统计并转化为按水平层划分的射线长度。
[0026]
进一步的，在所述步骤八至步骤十中，更新走时残差的公式如下：
[0027]
δt
new
＝δt
old-gδs
ꢀꢀꢀ
(4)，
[0028]
其中，δt
new
为更新后的走时残差，δt
old
为更新前的走时残差，g为敏感度矩阵，δs为慢度更新量。
[0029]
进一步的，在所述步骤九中，初次反演前设定初始尺度、迭代次数、尺度减小倍数，下次迭代中，尺度以设定的倍数减小，直到达到迭代次数。
[0030]
进一步的，在步骤九中，每次迭代中，将射线路径投影到当前尺度的网格中，并统计各网格内的射线长度，以构建符合当前尺度的方程组。
[0031]
进一步的，在步骤九中，构建的反演方程组如式(3)所示，反演中采用并行lsqr算法求解。
[0032]
进一步的，在所述步骤七和步骤九中，反演使用的是经过预处理的测井速度。
[0033]
本发明实施例具有以下技术效果：
[0034]
本发明实施例为一种成像域井震联合多尺度层析反演方法，可以构建精细的偏移速度场。与常规偏移速度分析相比，本方法应用测井信息、构造信息，所得的速度模型更符合认识，为深度域精确成像奠定基础，具有广阔的应用前景。
附图说明
[0035]
图1为本发明实施例基于成像域井震联合多尺度层析反演方法流程图；
[0036]
图2为先验速度模型及其偏移结果显示图；
[0037]
图2(a)显示为先验速度模型，图2(b)显示为基于先验模型得到的深度偏移剖面；
[0038]
图3为基于深度偏移剖面倾斜叠加得到的倾角场显示图，图中黑色背景条纹为偏移剖面中的地震波同相轴；
[0039]
图4为基于深度偏移剖面γ扫描得到的γ场显示图，其中，白色背景条纹为偏移剖面中的地震波同相轴；
[0040]
图5为三维射线覆盖显示图；
[0041]
图6为工区测井速度曲线显示图；
[0042]
图7为以井为中心开展的区域划分图；
[0043]
图8为反演结果的对比显示图，其中，图8(a)为常规方法的结果，图8(b)为本文方法的结果，图中黑圈内为对比区域。
具体实施方式
[0044]
本部分将详细描述本发明的具体实施例，本发明之较佳实施例在附图中示出，附图的作用在于用图形补充说明书文字部分的描述，使人能够直观地、形象地理解本发明的每个技术特征和整体技术方案，但其不能理解为对本发明保护范围的限制。
[0045]
如图1所示，图1为本发明基于成像域井震联合多尺度层析反演方法流程图。该方法包括以下步骤：
[0046]
具体步骤如下：
[0047]
步骤一，基于先验速度，进行叠前深度偏移，获取深度偏移剖面和偏移距共成像点道集.
[0048]
本发明实施例中，如图2(a)所示为先验速度模型，基于该先验速度模型，进行叠前深度偏移，获取深度偏移剖面，如图2(b)所示。
[0049]
步骤二，在偏移剖面中做局部倾斜叠加，得x向倾角场和y向倾角场。
[0050]
本发明实施例中，在偏移剖面中做局部倾斜叠加，得x向倾角场和y向倾角场显示图，图中黑色背景条纹为偏移剖面中的地震波同相轴。
[0051]
步骤三，对各所述成像点道集进行γ扫描，获得γ场；
[0052]
本发明实施例中，如图4所示为基于深度偏移剖面γ扫描得到的γ场显示图，图中白色背景条纹为偏移剖面中的地震波同相轴。
[0053]
可以理解的是，本发明实施例中步骤二和步骤三可同时进行也可以按先后顺序进行。
[0054]
步骤四，基于所述倾角场和所述γ场进行综合分析和质量监控，筛选出所述深度偏移剖面中的反射点。
[0055]
步骤五，以各所述反射点作为出射点，进行射线追踪到地表，获得射线路径并基于射线终点之间的偏移距计算剩余深度差；
[0056]
本发明实施例中，如图5所示为三维射线覆盖显示图。
[0057]
步骤六，以各井为中心，划分反演区域；
[0058]
具体的，在所述步骤六，以各井为中心，划分反演区域以各井为中心，具体包括：以给定的宽度划分整个工区区域，以距离井最近的区域作为首先反演的区域，以距离井次近的区域作为下一次反演的区域，以此类推，由近及远地划分整个工区范围。
[0059]
在本发明的一个具体示例中，以各井为中心，划分反演区域，其中，测井数据共56口，图6展示的是测井速度曲线，图7展示的是基于井的区域划分，确定了局部反演的先后顺序，由近及远，覆盖整个工区；
[0060]
步骤七，确定所述反演区域，基于射线路径及剩余深度差以及测井速度信息，构建局部区域的方程组，求解所述方程组，更新局部速度；
[0061]
具体的，求解所述方程组采用的目标函数o的具体形式如下：
[0062]
o＝cd||δt-gδs||2 ε1||s δs-s
well
||2 ε2||tr(s δs)||2ꢀꢀꢀ
(1)
[0063]
其中，g为由射线路径构成的敏感度矩阵，δs为慢度更新量，δt为走时残差，cd为数据协方差矩阵，s为模型慢度，s
well
为测井慢度，t为构造方向旋转矩阵，r为正则化算子，ε1、ε2为权重系数，上式中，单项走时残差的计算公式为：
[0064]
δt＝2sc·
δz
·
cosβ
·
cosα
ꢀꢀꢀ
(2)
[0065]
sc为反射点慢度，β为出射角，α为地层倾角，δz为剩余深度差；
[0066]
所述目标函数对应的所述方程组如下：
[0067][0068]
在本发明的优选实施例中，在所述步骤七中，进行的局部速度反演为一维反演，反演前将各个网格的射线长度统计并转化为按水平层划分的射线长度。
[0069]
步骤八，判断所述区域是否全部完成反演，若是则更新走时残差，否则迭代步骤七。
[0070]
步骤九，确定反演网格尺度，重整射线路径，并结合平滑约束以及测井约束，构建反演方程组并求解，获取模型更新量。
[0071]
步骤十，判断是否更新所述尺度，若是，更新走时残差并迭代计算所述步骤九，若否，求和各所述尺度更新量，更新模型并输出。
[0072]
具体的，在所述步骤八至步骤十中，更新走时残差的公式如下：
[0073]
δt
new
＝δt
old-gδs
ꢀꢀꢀ
(4)，
[0074]
其中，δt
new
为更新后的走时残差，δt
old
为更新前的走时残差，g为敏感度矩阵，δs为慢度更新量。
[0075]
本发明实施例中，在所述步骤九中，初次反演前设定初始尺度、迭代次数、尺度减
小倍数，下次迭代中，尺度以设定的倍数减小，直到达到迭代次数。
[0076]
优选的，在步骤九中，每次迭代中，将射线路径投影到当前尺度的网格中，并统计各网格内的射线长度，以构建符合当前尺度的方程组。
[0077]
优选的，在步骤九中，构建的反演方程组如式(3)所示，反演中采用并行lsqr算法求解。
[0078]
具体的，在所述步骤七和步骤九中，反演使用的是经过预处理的测井速度。
[0079]
可以理解的是，预处理可以包括去除异常值、平滑等。
[0080]
图8为反演结果的对比显示图，其中，图8(a)为常规方法的反演速度模型，图8(b)为本发明实施例基于成像域井震联合多尺度层析反演方法的反演速度模型。对比两个速度模型图8(a)和图8(b)可知，本发明实施例方法的结果具有更高的分辨率，且更符合地下构造，能将黑圈内的速度反转情况很好的呈现出来。该例证明本发明实施例方法是一种有效的成像域层析反演方法。
[0081]
以上所述仅为本发明的较佳实施例而己，并不以本发明为限制，凡在本发明的精神和原则之内所作的均等修改、等同替换和改进等，均应包含在本发明的专利涵盖范围内。