基于缓变各向异性的井控速度反演方法与流程

1.本发明涉及油气勘探地震资料处理技术领域，特别涉及一种基于缓变各向异性的井控速度反演方法。


背景技术：

2.目前，生产中常用的速度建模方法是基于射线理论的剩余曲率速度分析，主要应用反射波的走时信息，反演的结果往往是以低波数段为主的背景速度场，难以满足高精度地震勘探的需求。如何提高速度模型的精度，实现复杂地区的精确成像是一个急需解决的关键性问题，也是地震成像领域乃至地震勘探领域的一大难点。基于波动方程理论的速度建模方法在理论上可以得到包含高频成分的速度场，分辨率高，能更好地适应速度变化剧烈的地区，但存在许多尚未解决的问题。波动方程偏移速度分析理论和实际应用尚未完善，初始模型问题和对速度模型的敏感度是该方法面临的一大问题，并且计算量庞大，处理分析不灵活。全波形反演虽然在理论上已基本完善，但其计算量太大，严重依赖低频信息和大偏移距地震数据。而对于陆上地震勘探，如何提取地震子波以及模拟复杂波场也是摆在全波形反演面前的难题。复杂地区的地震资料还存在信噪比低的问题，因此，应用全波形反演解决复杂地区的速度建模问题是短期内无法实现的。
3.除上述问题外，许多工区还存在各向异性现象。各向异性作用和速度相互耦合，显著影响深度成像的准确性。


技术实现要素：

4.本发明的目的是针对常规偏移速度分析结果分辨率低、各向异性和速度相互耦合的问题，提供一种基于缓变各向异性的井控速度反演方法，以反演精细的地下速度场，为深度域地震成像提供技术支持。
5.本发明提供一种基于缓变各向异性的井控速度反演方法，该方法包括：
6.步骤一，整理工区内的测井速度数据，去除其中的异常值，并做平滑处理；
7.步骤二，基于时深转换得到的先验速度模型，对预处理的叠前地震数据进行射线类三维叠前深度偏移，获取深度偏移剖面mig(x,y,z)和偏移距共成像点道集cip(x,y,z,h)，其中，x、y、z、h分别为成像点的主测线方向坐标、联络测线方向坐标、深度和偏移距；
8.步骤三，在所述偏移剖面mig(x,y,z)中对所述成像点进行x/y方向局部倾斜叠加扫描，以最大相似谱对应的倾角作为该点倾角，逐点扫描得x方向倾角场和y方向倾角场
9.步骤四，对各所述成像点道集cip(x,y,z,h)进行γ谱扫描，获得最大相似谱对应的γ场γ(x,y,z)；
10.步骤五，基于所述倾角场和γ场进行综合分析和质量监控，筛选出所述深度偏移剖面mig(x,y,z)中的反射点；
11.步骤六，以各所述反射点作为出射点，进行射线追踪到地表，获得射线路径并基于
射线终点之间的偏移距计算剩余深度差；
12.步骤七，将预处理的测井速度沿所述倾角场和所述外推，加权求和得到井速度插值模型；
13.步骤八，基于射线路径及对应的深度差、平滑约束以及从所述井速度插值模型中提取的信息，构建反演方程组；
14.步骤九，求解所述反演方程组，获取所述速度模型更新量，更新所述速度模型。
15.进一步的，在步骤三中，x/y方向的相似谱的计算公式如下：
[0016][0017]
其中，sem
dipx
、sem
dipy
分别为x、y方向的相似谱，z为深度，k为倾角的正切值，trace为地震道，x0、y0为中心道位置，xj、yj为倾斜叠加的地震道位置，l为扫描窗口道数。
[0018]
进一步的，在步骤四中，γ谱扫描所采用的成像深度z
mig
与γ的关系
[0019]
式为：
[0020][0021]
γ相似谱的计算公式如下：
[0022][0023]
其中，sem
rmo
为γ相似谱，δz为加权平均的窗口宽度。
[0024]
进一步的，在步骤六中，从拾取的倾角场中获得反射点处的x/y方向倾角，计算反射点法线方向的射线以及在法线两侧以相同角度出射的射线，每次计算均以不同出射角和方位角进行射线追踪，并以到达地表边界作为该条射线的终止条件。
[0025]
进一步的，在步骤六中，具体采用的三维射线追踪公式如下：
[0026][0027]
其中，p
x
、py、pz为慢度分量，v是模型速度，上述方程组进一步通过runge-kutta法进行数值求解。
[0028]
进一步的，在步骤六中，剩余深度差δz的计算公式如下：
[0029][0030]
其中，z0为反射点在偏移距成像点道集中的零偏移距深度。
[0031]
进一步的，在步骤八中，基于射线路径及对应的深度差构建的方程如下：
[0032]
l0δs
0 cosβ
h-lhδsh＝2sc·
cosβh·
cosα
·
δz，
[0033]
其中，l0为零偏移距射线路径在各网格内的长度，δs0为零偏移距射线路径上的慢度更新量，βh为偏移距h对应的出射角，lh为偏移距h的射线路径在各网格内的长度，δsh为偏移距h的射线路径上的慢度更新量，sc为反射点慢度，α为地层倾角，δz为剩余深度差。
[0034]
进一步的，在步骤八中，从所述井速度插值模型中提取相邻速度的相对变化约束反演过程，具体方程如下：
[0035][0036]
其中，si、δsi、分别为第i点慢度、慢度更新量和井慢度，s
i+1
、δs
i+1
、分别为第i+1点慢度、慢度更新量和井慢度。
[0037]
进一步的，在所述步骤八中，构建的反演方程组形式如下：
[0038][0039]
其中，l、δt为由权利要求7中的方程构成的系数矩阵和右端项，r为由三维拉普拉斯算子构成的平滑矩阵，w、e为由权利要求8中的方程构成的系数矩阵和右端项，δs表示模型更新量，ε1、ε2为权重系数。
[0040]
进一步的，在所述步骤九中，采用并行lsqr算法求解大型稀疏的反演方程组。
[0041]
本发明有益技术效果如下：
[0042]
本发明实施例提供一种基于缓变各向异性的井控速度反演方法，可以构建偏移速度场。与常规偏移速度分析相比，本发明实施例方法应用测井信息、构造信息，并考虑各向异性影响，所得的各向同性速度模型更符合实际情况，分辨率更高，为深度域的正确成像奠定基础，具有广阔的应用前景。
附图说明
[0043]
图1为本发明实施例中基于缓变各向异性的井控速度反演方法的流程图；
[0044]
图2为本发明实施例中工区内某口井预处理后的测井速度曲线显示图；
[0045]
图3为图1中所述方法的先验速度模型及其偏移结果显示图，其中，图3(a)显示为先验速度模型，图3(b)显示为基于先验模型得到的深度偏移剖面，图3(c)显示为基于先验模型得到的某成像点道集；
[0046]
图4为本发明实施例中基于深度偏移剖面倾斜叠加得到的倾角场显示图，图中黑色背景条纹为偏移剖面中的地震波同相轴；
[0047]
图5为本发明实施例中基于深度偏移剖面γ扫描得到的γ场显示图，图中黑色背景条纹为偏移剖面中的地震波同相轴；
[0048]
图6为本发明实施例中三维射线覆盖显示图；
[0049]
图7为本发明实施例中工区井速度模型显示图，图中黑色背景条纹为偏移剖面中的地震波同相轴；
[0050]
图8为本发明实施例中常规偏移速度反演的最终结果显示图，其中图8(a)为反演速度模型，其中黑圈内为对比区域；图8(b)为最终的深度偏移剖面；图8(c)为最终的成像点道集；
[0051]
图9为图1所示方法的最终结果显示图；
[0052]
图9(a)为反演速度模型，其中黑圈内为对比区域；
[0053]
图9(b)为最终的深度偏移剖面；
[0054]
图9(c)为最终的成像点道集。
具体实施方式
[0055]
本部分将详细描述本发明的具体实施例，本发明之较佳实施例在附图中示出，附图的作用在于用图形补充说明书文字部分的描述，使人能够直观地、形象地理解本发明的每个技术特征和整体技术方案，但其不能理解为对本发明保护范围的限制。
[0056]
如图所示，图1为本发明基于缓变各向异性的井控速度反演方法的流程图。该方法包括以下步骤：
[0057]
步骤一，整理工区内的测井速度数据，去除其中的异常值，并做平滑处理，图2显示为其中一口井处理后的井速度曲线；
[0058]
步骤二，基于时深转换得到先验速度模型，如图3(a)所示。对预处理的叠前地震数据进行射线类三维叠前深度偏移，获取深度偏移剖面mig(x,y,z)如图3(b)所示，和偏移距共成像点道集cip(x,y,z,h)，如图3(c)所示，其中，x、y、z、h分别为成像点的主测线方向坐标、联络测线方向坐标、深度和偏移距。
[0059]
步骤三，在所述偏移剖面mig(x,y,z)中对所述成像点进行x/y方向局部倾斜叠加扫描，以最大相似谱对应的倾角作为该点倾角，逐点扫描得x方向倾角场和y方向倾角场如图4所示。
[0060]
优选的，x/y方向的相似谱的计算公式如下：
[0061][0062]
其中，sem
dipx
、sem
dipy
分别为x、y方向的相似谱，z为深度，k为倾角的正切值，trace为地震道，x0、y0为中心道位置，xj、yj为倾斜叠加的地震道位置，l为扫描窗口道数。
[0063]
步骤四，对各所述成像点道集cip(x,y,z,h)进行γ谱扫描，获得最大相似谱对应的γ场γ(x,y,z)，如图5所示。
[0064]
优选的，γ谱扫描所采用的成像深度z
mig
与γ的关系式为：
[0065][0066]
γ相似谱的计算公式如下：
[0067][0068]
其中，sem
rmo
为γ相似谱，δz为加权平均的窗口宽度。
[0069]
步骤五，基于所述倾角场和γ场进行综合分析和质量监控，筛选出所述深度偏移剖面mig(x,y,z)中的反射点。
[0070]
步骤六，以各所述反射点作为出射点，进行射线追踪到地表，获得射线路径并基于射线终点之间的偏移距计算剩余深度差，如图6为三维射线覆盖显示图。
[0071]
具体地，从拾取的倾角场中获得反射点处的x/y方向倾角，计算反射点法线方向的射线以及在法线两侧以相同角度出射的射线，每次计算均以不同出射角和方位角进行射线追踪，并以到达地表边界作为该条射线的终止条件。
[0072]
作为本发明的一个优选实施例，在步骤六中，具体采用的三维射线追踪公式如下：
[0073][0074]
其中，p
x
、py、pz为慢度分量，v是模型速度，上述方程组进一步通过runge-kutta法进行数值求解。
[0075]
进一步的，在步骤六中，剩余深度差δz的计算公式如下：
[0076][0077]
其中，z0为反射点在偏移距成像点道集中的零偏移距深度。
[0078]
步骤七，将预处理的测井速度沿所述倾角场和所述外推，加权求和得到井速度插值模型。
[0079]
具体地，对工区内的测井速度进行处理，如去除异常值、平滑，并沿倾角场外推，加权求和得到井速度模型，如图7所示。
[0080]
步骤八，基于射线路径及对应的深度差、平滑约束以及从所述井速度插值模型中提取的信息，构建反演方程组。
[0081]
进一步的，在步骤八中，基于射线路径及对应的深度差构建的方程如下：
[0082]
l0δs
0 cosβ
h-lhδsh＝2sc·
cosβh·
cosα
·
δz，
[0083]
其中，l0为零偏移距射线路径在各网格内的长度，δs0为零偏移距射线路径上的慢度更新量，βh为偏移距h对应的出射角，lh为偏移距h的射线路径在各网格内的长度，δsh为偏移距h的射线路径上的慢度更新量，sc为反射点慢度，α为地层倾角，δz为剩余深度差。
[0084]
作为本发明的一个具体实施例，在步骤八中，从井速度插值模型中提取相邻速度的相对变化约束反演过程，具体方程如下：
[0085][0086]
其中，si、δsi、分别为第i点慢度、慢度更新量和井慢度，s
i+1
、δs
i+1
、分别为第i+1点慢度、慢度更新量和井慢度。
[0087]
进一步的，在所述步骤八中，构建的反演方程组形式如下：
[0088][0089]
其中，l、δt为由权利要求7中的方程构成的系数矩阵和右端项，r为由三维拉普拉斯算子构成的平滑矩阵，w、e为由权利要求8中的方程构成的系数矩阵和右端项，δs表示模型更新量，ε1、ε2为权重系数。
[0090]
步骤九，求解所述反演方程组，获取速度模型更新量，更新速度模型。
[0091]
优选的，在所述步骤九中，采用并行lsqr算法求解大型稀疏的反演方程组。
[0092]
经两次迭代后，可得到最终的结果。常规偏移速度反演的最终结果如图8所示，其中，图8(a)为反演速度模型，图8(b)为最终的深度偏移剖面，图8(c)为最终的成像点道集。本发明实施例方法的最终结果如图9所示，其中，图9(a)为反演速度模型，图9(b)为最终的深度偏移剖面，图9(c)为最终的成像点道集。对比图8和图9可知，两种方法的结果均可使同相轴拉平，使深度偏移进行有效成像。而对比两个速度模型图8(a)和图9(a)可知，本发明实施例方法的结果具有更高的分辨率，尤其是黑圈内的区域。
[0093]
对比深度偏移剖面可知，本发明实施例方法的速度模型更符合成像剖面中的构造趋势。该具体示例证明了本发明实施例方法是一种有效的成像域层析反演方法。
[0094]
以上所述仅为本发明的较佳实施例而己，并不以本发明为限制，凡在本发明的精神和原则之内所作的均等修改、等同替换和改进等，均应包含在本发明的专利涵盖范围内。