CN116699678A - 一种多信息融合的速度建模方法及装置 - Google Patents
一种多信息融合的速度建模方法及装置 Download PDFInfo
- Publication number
- CN116699678A CN116699678A CN202210167182.6A CN202210167182A CN116699678A CN 116699678 A CN116699678 A CN 116699678A CN 202210167182 A CN202210167182 A CN 202210167182A CN 116699678 A CN116699678 A CN 116699678A
- Authority
- CN
- China
- Prior art keywords
- data
- seismic
- logging
- interpolation
- speed
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Pending
Links
- 238000000034 method Methods 0.000 title claims abstract description 103
- 230000004927 fusion Effects 0.000 title claims abstract description 69
- 230000005012 migration Effects 0.000 claims abstract description 43
- 238000013508 migration Methods 0.000 claims abstract description 43
- 238000003384 imaging method Methods 0.000 claims abstract description 40
- 238000012545 processing Methods 0.000 claims abstract description 36
- 238000013213 extrapolation Methods 0.000 claims abstract description 32
- 238000010276 construction Methods 0.000 claims abstract description 16
- 238000003708 edge detection Methods 0.000 claims abstract description 8
- 238000005516 engineering process Methods 0.000 claims abstract description 8
- 238000001914 filtration Methods 0.000 claims description 51
- 238000009499 grossing Methods 0.000 claims description 15
- 238000013528 artificial neural network Methods 0.000 claims description 13
- 230000008569 process Effects 0.000 claims description 8
- 208000035126 Facies Diseases 0.000 claims description 6
- 239000011435 rock Substances 0.000 claims description 5
- 230000010365 information processing Effects 0.000 claims description 3
- 230000006870 function Effects 0.000 description 16
- 230000008901 benefit Effects 0.000 description 4
- 238000004364 calculation method Methods 0.000 description 4
- 230000008859 change Effects 0.000 description 4
- 238000010586 diagram Methods 0.000 description 4
- 230000000694 effects Effects 0.000 description 3
- 238000007499 fusion processing Methods 0.000 description 3
- 238000013507 mapping Methods 0.000 description 3
- 238000007500 overflow downdraw method Methods 0.000 description 3
- 230000004044 response Effects 0.000 description 3
- 230000015572 biosynthetic process Effects 0.000 description 2
- 238000002939 conjugate gradient method Methods 0.000 description 2
- 230000001419 dependent effect Effects 0.000 description 2
- 230000001747 exhibiting effect Effects 0.000 description 2
- 238000012986 modification Methods 0.000 description 2
- 230000004048 modification Effects 0.000 description 2
- 238000005070 sampling Methods 0.000 description 2
- NAWXUBYGYWOOIX-SFHVURJKSA-N (2s)-2-[[4-[2-(2,4-diaminoquinazolin-6-yl)ethyl]benzoyl]amino]-4-methylidenepentanedioic acid Chemical compound C1=CC2=NC(N)=NC(N)=C2C=C1CCC1=CC=C(C(=O)N[C@@H](CC(=C)C(O)=O)C(O)=O)C=C1 NAWXUBYGYWOOIX-SFHVURJKSA-N 0.000 description 1
- 101100269850 Caenorhabditis elegans mask-1 gene Proteins 0.000 description 1
- 238000010521 absorption reaction Methods 0.000 description 1
- 230000009471 action Effects 0.000 description 1
- 230000007547 defect Effects 0.000 description 1
- 238000013461 design Methods 0.000 description 1
- 238000001514 detection method Methods 0.000 description 1
- 238000011161 development Methods 0.000 description 1
- 238000005553 drilling Methods 0.000 description 1
- 230000005484 gravity Effects 0.000 description 1
- 230000006872 improvement Effects 0.000 description 1
- 230000000877 morphologic effect Effects 0.000 description 1
- 230000035939 shock Effects 0.000 description 1
- 238000012876 topography Methods 0.000 description 1
- 230000007704 transition Effects 0.000 description 1
Classifications
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
- G01V1/00—Seismology; Seismic or acoustic prospecting or detecting
- G01V1/28—Processing seismic data, e.g. for interpretation or for event detection
- G01V1/282—Application of seismic models, synthetic seismograms
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
- G01V1/00—Seismology; Seismic or acoustic prospecting or detecting
- G01V1/28—Processing seismic data, e.g. for interpretation or for event detection
- G01V1/30—Analysis
- G01V1/303—Analysis for determining velocity profiles or travel times
Landscapes
- Engineering & Computer Science (AREA)
- Remote Sensing (AREA)
- Physics & Mathematics (AREA)
- Life Sciences & Earth Sciences (AREA)
- Acoustics & Sound (AREA)
- Environmental & Geological Engineering (AREA)
- Geology (AREA)
- General Life Sciences & Earth Sciences (AREA)
- General Physics & Mathematics (AREA)
- Geophysics (AREA)
- Geophysics And Detection Of Objects (AREA)
Abstract
本发明公开了一种多信息融合的速度建模方法与装置，所述方法包括，获取地震偏移处理后的成像数据A、数据相控类型蒙板B、地震构造解释层位和断层数据C、测井速度曲线数据D和地震速度场数据E；通过边缘检测技术提取成像数据A的边缘特征，得到地震图像的构造骨架数据F；将相控类型蒙板B与构造骨架数据F进行融合，获得相控构造导向数据G；在相控构造导向数据G中加入地震构造解释数据C，形成空间插值引导数据H；通过两步插值外推，获得与空间插值引导数据H相对应的空间插值距离权场Q和测井速度插值填充模型P；为测井速度插值填充模型P建立分类蒙板M；通过对类蒙板M进行两步填充，获得多信息融合模型R。
Description
技术领域
本发明涉及地震资料处理技术领域，尤其涉及一种多信息融合的速度建模方法及装置。
背景技术
实现复杂构造和储层成像的核心方法是叠前深度偏移技术。叠前深度偏移成像的准确性极大的依赖于深度域速度建模的可靠性。中国陆上复杂探区存在地形起伏剧烈，地表出露岩性变化快，近地表厚度、速度横向变化大，地下构造破碎，断层发育，地层倾角大，产状变化多，速度横向变化剧烈，局部垂向速度反转等复杂地震地质特征，迫切需要发挥不同勘探资料(地质、地震、非地震)优势，通过多种信息融合建模的技术手段有效应对复杂性。
地质认识(构造、岩性)、地震速度场和测井速度的融合建模是一项技术挑战。地质认识由层位、断层和岩性体边界特征组成，提供地下大尺度地质框架结构信息；地震速度场具有中低频特征，横向连续性好，大套速度结构与形态特征明显，但纵向分辨率低，主要提供大于菲涅尔带尺度的速度信息；测井速度具有高频特征，纵向分辨率高，但横向探测范围小，主要提供局部范围内小尺度速度高频细节信息。可见，在地质认识(构造、岩性)、地震速度场和测井速度融合建模过程中，不同尺度的信息融合是需要解决的关键技术问题。
常规实现地质构造、地震速度场和测井速度融合的方式有两种：一种是将测井速度平滑到与地震速度相当的尺度，再通过适当的插值方法将平滑后的测井速度沿层外推填充模型，在填充区域内，得到与平滑的测井速度细节完全一致的速度模型；另一种是用神经网络的方式学习地震速度和测井速度之间的映射关系，然后将这种局部映射关系应用到整个模型，获得带有细节特征的速度模型。这两种技术思路存在三个共性问题；一是，没有充分发挥地震速度场的作用；二是，试图用井位处测井速度的“一孔之见”表达地下不同位置的复杂情况，致使无论是直接插值填充形成的模型细节，还是局部映射关系全局应用形成的细节都与实际情况存在较大差异；三是，两种思路都不考虑岩性相带控制。
发明内容
本发明提出一种多信息融合的速度建模方法及装置，该方法充分考虑地质认识(构造、岩性)、地震速度场和测井速度资料特点，一方面综合利用地质构造和岩性相带特征对井位处测井速度外推范围和形态加以控制，另一方面充分考虑地震速度场的作用，实现地质认识(包括构造和岩性相带)、测井速度和地震速度场的自适应加权融合。为实现上述目的，本发明提供如下技术方案：
一种多信息融合的速度建模方法，所述方法包括以下步骤，
获取地震偏移处理后的成像数据A；
获取岩相分类解释数据B，称为相控类型蒙板；
获取地震构造解释层位和断层数据C；
获取测井速度曲线数据D，并进行平滑处理；
获取地震速度场数据E；
通过边缘检测技术提取所述成像数据A的边缘特征，得到地震图像的构造骨架数据F；
将所述相控类型蒙板B与所述构造骨架数据F进行融合，获得相控构造导向数据G；
在所述相控构造导向数据G中加入地震构造解释数据C，形成空间插值引导数据H；
通过两步插值外推，获得与空间插值引导数据H相对应的空间插值距离权场Q和测井速度插值填充模型P；
为测井速度插值填充模型P建立分类蒙板M；
通过对类蒙板M进行两步填充，获得多信息融合模型R。
优选的，所述成像数据A是通过用户提供的叠前地震道集数据体和用户提供的偏移速度场获取的。
优选的，所述岩相分类解释数据B是根据成像数据A和用户提供的工区岩心录井柱状图/测井岩相图数据，通过测井岩相内插外推法、测井岩相神经网络智能预测法或地震数据地震相智能预测法获得的。
优选的，所述地震构造解释层位和断层数据C是通过人工解释法或神经网络智能解释法对所述成像数据A进行解释获得的。
优选的，所述测井速度曲线数据D是通过用户提供的测井速度数据进行平滑处理得到的。
优选的，所述平滑处理包括，对测井速度数据进行中值滤波；对中值滤波所得数据做高斯递归滤波平滑。
优选的，所述中值滤波包括对一个含奇数个点的滑动窗口的数值进行排序，将中值赋值给中心点；所述中值滤波的滑动窗口形状包括线状、方形、圆形、十字形或用户提供的特殊几何形状。
优选的，所述高斯递归滤波步骤包括：对中值滤波数据进行前向滤波；对前向滤波数据进行后向滤波；其中，
所述前向滤波公式为：
GSForward(n)＝a1*Input1(n)+a2*Input1(n-1)-a3*GSForward(n-1)-a4*GSForward(n-2)
式中，GSForward是前向滤波生成的数据，Input1是前向滤波输入数据，这里为中值滤波数据，n是数据样点位置，a1、a2、a3、a4是滤波器系数，根据不同工区实际需要进行设定；
所述后向滤波公式为：
GSBackward(n)＝b1*Input2(n)+b2*Input2(n+1)-b3*GSBackward(n+1)-b4*GSBackward(n+2)
式中，GSBackward是后向滤波生成的数据，Input2是后向滤波输入数据，这里为前向滤波数据，n是数据样点位置，b1、b2、b3、b4是滤波器系数，根据不同工区实际需要进行设定。
优选的，所述地震速度场数据E是与所述地震偏移处理后的成像数据A对应的偏移速度场数据，通过用户提供的叠前地震道集数据体Data处理获得。
优选的，所述通过两步插值外推获得与空间插值引导数据H相对应的空间插值距离权场Q和测井速度插值填充模型P包括，
在空间插值引导数据H中加入测井速度曲线数据D，判断测井数据空间和时间范围及与相控类型蒙板B最近邻相控区域的位置关系。
优选的，所述两步插值外推具体方法如下：按照空间插值引导数据H提供的局部构造骨架特征，
通过径向基函数方法对落入相控区域的测井数据进行外推处理；
通过克里金插值法建立所有测井数据的外推关系，实现测井速度沿构造和相关相带的外推填充。
优选的，所述径向基函数方法公式如下：
其中，ε表示径向基函数的尺度因子，这里取1；r表示插值点和井位置的欧式距离。
优选的，所述分类蒙板M的区域划分为0值指示区和1值指示区域；
所述0值指示区域表示测井速度插值填充模型P的未填充区域；
所述1值指示区域表示通过加权最小二乘法对地震速度模型和测井速度插值填充模型进行加权融合后的速度值的填充区域。
一种多信息融合的速度建模方法的装置，所述装置包括，
信息收集装置用于获取地震偏移处理后的成像数据A、岩相分类解释数据B、地震构造解释层位和断层数据C、测井速度曲线数据D、获取地震速度场数据E；
信息处理装置用于通过边缘检测技术提取所述成像数据A的边缘特征，得到地震图像的构造骨架数据F，将所述相控类型蒙板B与所述构造骨架数据F进行融合，获得相控构造导向数据G，在所述相控构造导向数据G中加入地震构造解释数据C，形成空间插值引导数据H；
信息建模装置用于通过两步插值外推，获得与空间插值引导数据H相对应的空间插值距离权场Q和测井速度插值填充模型P，为测井速度插值填充模型P建立分类蒙板M，通过对分类蒙板M进行两步填充，获得多信息融合模型R。
本发明的技术效果和优点：
本发明充分考虑地质认识(构造、岩性)、地震速度场和测井速度资料特点，一方面综合利用地质构造和岩性相带特征对井位处测井速度外推范围和形态加以控制，另一方面充分考虑地震速度场的作用，实现地质认识(包括构造和岩性相带)、测井速度和地震速度场的自适应加权融合。
本发明的其它特征和优点将在随后的说明书中阐述，并且，部分地从说明书中变得显而易见，或者通过实施本发明而了解。本发明的目的和其他优点可通过在说明书、权利要求书以及附图中所指出的结构来实现和获得。
附图说明
图1是本发明提供的多信息融合速度建模方法流程图；
图2是本实施例中所采用的地震成像数据图；
图3a为本实施例中用户提供的表示岩相的地震属性示意图；
图3b为本实施例中人工智能地质体雕刻形成的表示岩相的封闭区域示意图；
图4是本实施例中所采用的地震构造解释数据图；
图5是本实施例中所采用的四口测井速度曲线数据与地震成像数据的位置关系图；
图6是本实施例中所采用的地震速度场数据图；
图7a为本实施例中1#测井位置处地震速度、原始测井速度和平滑后的测井速度的叠合显示图；
图7b为本实施例中2#测井位置处地震速度、原始测井速度和平滑后的测井速度的叠合显示图；
图7c为本实施例中3#测井位置处地震速度、原始测井速度和平滑后的测井速度的叠合显示图；
图7d为本实施例中4#测井位置处地震速度、原始测井速度和平滑后的测井速度的叠合显示图；
图8是本发明所获得的空间插值距离权场图；
图9是本发明中获得的测井速度插值填充模型图；
图10是本发明中所得的相对中频的融合速度模型图；
图11是本发明中所得的相对高频的融合速度模型图；
图12是本发明中测井位置处的地震速度、测井插值速度和相对中频融合速度的比较图；
图13是本发明中测井位置处的地震速度、测井插值速度和相对高频融合速度的比较图。
具体实施方式
下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。
为解决现有技术的不足，本发明公开了一种多信息速度融合方法及装置，所述装置包括，信息收集装置用于获取地震偏移处理后的成像数据A、岩相分类解释数据B、地震构造解释层位和断层数据C、测井速度曲线数据D、获取地震速度场数据E；信息处理装置用于通过边缘检测技术提取所述成像数据A的边缘特征，得到地震图像的构造骨架数据F，将所述相控类型蒙板B与所述构造骨架数据F进行融合，获得相控构造导向数据G，在所述相控构造导向数据G中加入地震构造解释数据C，形成空间插值引导数据H；信息建模装置用于通过两步插值外推，获得与空间插值引导数据H相对应的空间插值距离权场Q和测井速度插值填充模型P，为测井速度插值填充模型P建立分类蒙板M，通过对分类蒙板M进行两步填充，获得多信息融合模型R。
如图1所示，本发明的方法包括以下步骤，获取地震偏移处理后的成像数据A；获取岩相分类解释数据B，称为相控类型蒙板；获取地震构造解释层位和断层数据C；获取测井速度曲线数据D，并进行平滑处理；获取地震速度场数据E；通过边缘检测技术提取所述成像数据A的边缘特征，得到地震图像的构造骨架数据F；将所述相控类型蒙板B与所述构造骨架数据F进行融合，获得相控构造导向数据G；在所述相控构造导向数据G中加入地震构造解释数据C，形成空间插值引导数据H；通过两步插值外推，获得与空间插值引导数据H相对应的空间插值距离权场Q和测井速度插值填充模型P；为测井速度插值填充模型P建立分类蒙板M；通过对分类蒙板M进行两步填充，获得多信息融合模型R。
进一步地，获取地震偏移处理后的成像数据A，如图2所示，图中，Xline(trace)表示地震道，time表示时间，Amplitude表示振幅，图像中存在明暗变化的位置指示了地下反射层位置及响应振幅，图中不同位置反射层的响应振幅大小与反射层对应的地质岩性相关，不同地质岩性对地震波的吸收衰减不同，因此成像响应振幅不相同。常用的地震叠前深度偏移方法有射线类偏移和波动方程偏移两大类，射线类方法主要包括Kirchhoff偏移、高斯束偏移、控制束偏移等；波动方程类方法主要包括单程波偏移和逆时偏移。射线类偏移和波动方程偏移两类方法都以波动方程为理论基础，不同之处在于射线类偏移利用几何射线理论来计算波场的振幅以及相位信息，从而实现波场的延拓成像，而波动方程类偏移则是基于波动方程的数值解法。一般来说，波动方程类偏移具有更高的成像精度，而射线类偏移具有更高的计算效率和灵活性。主流地震资料处理商业软件，如CGG，Omega，AGT，Paradigm等都提供面向两类偏移技术的处理模块，同时，开源软件Seismic Unix和Madagascar也提供面向两类偏移技术的处理模块，用户提供叠前地震道集数据体Data和偏移速度场Velocity的条件下，任意选择处理软件的偏移算法模块，即可获得所述地震偏移处理后的成像数据A，方法如下：
A＝ImageModule(Data,Velocity) (1)
其中ImageModule是商业软件或开源软件的成像处理模块，Data是用户提供的叠前地震道集数据体，Velocity是用户提供的偏移速度场，A是所述经过地震偏移处理后获得的成像数据体。
进一步地，获取岩相分类解释数据B，称为相控类型蒙板；获得岩相分类解释数据的方法有测井岩相内插外推法、测井岩相神经网络智能预测法、地震数据地震相智能预测法三大类，可根据不同工区实际资料情况，选择合适的方法模块产生岩相分类解释数据B。工作流程描述为：首先，将工区钻/测井位置的岩心录井柱状图/测井岩相图Well_Lithology_Model与所述地震偏移处理后的成像数据A进行标定匹配，然后，通过测井岩相内插外推法或测井岩相神经网络智能预测法或地震数据地震相智能预测法获得岩相分类解释数据B，方法如下：
B＝LithologyModule(A,Well_Lithology_Model) (2)
式中，LithologyModule是测井岩相内插外推法、测井岩相神经网络智能预测法或地震数据地震相智能预测法处理模块，A是所述经过地震偏移处理后获得的成像数据体，Well_Lithology_Model是用户提供的工区岩心录井柱状图/测井岩相图数据，B是所述岩相分类解释获得的岩相分类解释数据体。
相控类型蒙板可以是如图3a所示表示岩相的地震属性特征，相控类型蒙板也可以是如图3b所示人工智能地质体雕刻形成的表示岩相的封闭区域，还可以是非地震方法(重力、磁法、电法、地球化学法)提供的各种类型岩性相带特征或属性特征，图中Xline(trace)表示地震道，time表示时间。本例提到的相控类型蒙板只用于提供关键地质岩相边界特征，对实现方案不做限定。
获取地震构造解释层位和断层数据C；获得地震构造解释层位和断层数据的方法有人工解释法、神经网络智能解释法两大类，可根据不同工区实际资料情况，选择合适的方法模块产生地震构造解释层位和断层数据C。工作流程描述为：通过人工解释法或神经网络智能解释法对所述经过地震偏移处理后获得的成像数据A进行解释，获得所述地震构造解释层位和断层数据C，方法如下：
C＝StructureModule(A) (3)
式中，StructureModule是人工解释法或神经网络智能解释法处理模块，A是所述经过地震偏移处理后获得的成像数据体，C是所述地震构造解释获得的地震构造解释层位和断层数据。如图4所示本实施例共包含2个断层，10个层位，图中，Xline(trace)表示地震道，time表示时间，Amplitude表示振幅。
进一步地，获取测井速度曲线数据D，进行适当平滑处理；如图5所示，图中，Xline(trace)表示地震道，time表示时间，Amplitude表示振幅。图中显示区域内4口测井(从左到右依次为1#测井，2#测井，3#测井和4#测井)速度曲线数据与地震成像数据标定后的位置关系，横向上，可见4口测井位于不同位置，纵向上，通过井震标定解释，确定4口测井速度曲线的起始位置高低不同，同时，测井速度曲线用垂直条状图表示，其上明暗关系反映速度值随深度的变化，用户获得测井速度数据WellVelocity后，可以通过主流地震资料处理商业软件CGG，Omega，AGT，Paradigm等或开源软件Seismic Unix和Madagascar等的滤波器处理模块获得平滑后的测井速度数据，平滑程度应根据不同工区的实际要求进行选择，这里不做限定，实施方法如下：
D＝FilterModule(WellVelocity) (4)
式中，FilterModule是商业软件或开源软件的滤波器处理模块，WellVelocity是用户提供的测井速度数据，D是所述经过滤波后获得的平滑的测井速度数据。
通常测井速度数据离散分布在空间中，没有任何的几何结构，甚至可能不位于均匀采样的地震网格上。对图5中的4口测井(从左到右依次为1#，2#，3#和4#)速度曲线数据进行平滑处理，测井速度平滑前后与地震速度的关系如图7a-7d所示，图中，横标velocity表示速度，Time表示时间，，根据图7a-7d可以看出其中波动较大的曲线表示原始测井速度数据、平滑线表示井点对应位置处的地震速度数据、波动较小线表示平滑后的测井速度数据，其中波动较小的线在波动较大的线中间。三条曲线对比可见，原始测井速度数据与地震速度数据(比较波动较大线和平滑线)两者在频率成分相差较大。地震速度数据呈现低频特征，而原始测井速度数据呈现极高频特征。平滑后的测井速度数据(波动较小线)既包含了一些速度细节信息，频率成分又不至于过高而与地震速度偏差太大。
这里所述平滑处理包括两步，首先对测井速度数据做中值滤波，然后对中值滤波所得数据做高斯递归滤波做平滑，其中，中值滤波通常采用一个含奇数个点的滑动窗口，用窗口的中的数值的中值来代替中心点的数值，也就是对这个窗口中的数值进行排序，然后将中值赋值给中心点即可。常用的中值滤波窗口形状有线状、方形、圆形以及十字形或用户提供的特殊几何形状等，工作方法描述为：
中值滤波数据＝MedianFilter(窗口形状，窗口大小) (5)
其中MedianFilter是商业软件或开源软件提供的中值滤波器。
所述高斯递归滤波器具体步骤分为两步：第一步，对MedianFilter获得的中值滤波数据做一次前向滤波，第二步，对第一步得到的前向滤波数据再做一次后向滤波。所述前向滤波公式为：
GSForward(n)＝a1*Input1(n)+a2*Input1(n-1)-a3*GSForward(n-1)-a4*GSForward(n-2) (6)
式中，GSForward是前向滤波生成的数据，Input1是前向滤波输入数据，这里为公式5得到的中值滤波数据，n是数据样点位置，a1,a2,a3,a4是滤波器系数，根据不同工区实际需要进行设定；
所述后向滤波公式为：
GSBackward(n)＝b1*Input2(n)+b2*Input2(n+1)-b3*GSBackward(n+1)-b4*GSBackward(n+2) (7)
式中，GSBackward是后向滤波生成的数据，Input2是后向滤波输入数据，这里为公式6得到的前向滤波数据，n是数据样点位置，b1,b2,b3,b4是滤波器系数，根据不同工区实际需要进行设定。
进一步地，获取地震速度场数据E，如图6所示，图6展示了地震速度模型，该速度模型是一种层速度模型，该模型浅部呈现低速特征，深部呈现高速特征，在断层左侧范围内呈现中高速特征，在断层右侧范围内呈现中低速特征。所述地震速度场数据是与所述地震偏移处理后的成像数据A对应的偏移速度场数据。用户提供叠前地震道集数据体Data的条件下，可以通过主流地震资料处理商业软件，如CGG，Omega，AGT，Paradigm等，或开源软件Seismic Unix和Madagascar等的偏移速度建模处理模块直接获得，方法如下：
E＝VelocityModule(Data) (8)
式中，VelocityModule是商业软件或开源软件的偏移速度建模处理模块，Data是用户提供的叠前地震道集数据体，E是所述经过偏移速度建模处理获得的地震速度场数据体。
进一步地，通过边缘检测技术提取成像数据A的边缘特征，得到地震图像的构造骨架数据F；将相控类型蒙板B与构造骨架数据F进行融合，获得相控构造导向数据G；在相控构造导向数据G中加入地震构造解释数据C，形成空间插值引导数据H。
进一步地，为了获得空间插值距离权场如图8所示和测井速度插值填充模型如图9所示，图中Xline(trace)表示地震道，Distance field表示距离场，interpolatedvelocity表示测井速度插值，time表示时间，图8所示空间插值距离权场示意图中的明暗变化表示距离井点不同位置的插值权重大小，通常在井位附近设定为最大权重，距离井点越远，插值权重越小；图9所示测井速度插值填充模型是井点处速度曲线数据根据图8所示空间插值距离权场提供的权重关系外推插值获得。在判断测井数据空间和时间范围及与相控类型蒙板主要相控区域的位置关系后，首先对落入主要相控区域的测井数据进行外推处理，即在相控区域内，按照空间插值引导数据H提供的局部构造骨架特征，通过径向基函数方法实现井位处测井速度局部外推填充，外推过程延伸到相控区域边界终止；其次，按照空间插值引导数据H提供的整体构造导向特征，通过克里金插值法建立所有测井数据的外推关系，实现测井速度沿构造和相关相带的外推填充。通过如上两步插值外推，获得与空间插值引导数据H相对应的空间插值距离权场Q和测井速度插值填充模型P；
径向基函数方法公式如下：
其中，ε表示径向基函数的尺度因子，这里取1；r表示插值点和井位置的欧式距离。通过该方法插值的速度具有构造信息。
进一步地，这里的径向基函数可以有很多选择，比如，可以选择常规径向基函数，也可以选择径向基函数神经网络，也可以根据具体需要进行新型径向基函数设计。在本实施例中，我们设计了一种距离相关的径向基函数，表达式如下：
公式(10)中，我们假设研究区共有N口测井，其中第i口测井落入主要相控区域，设第i口测井的速度曲线段上共有Ki个采样点，(xk，yk，zk)表示速度曲线段上每个采样点坐标，(x，y，z)表示第i口测井落入主要相控区域的待插值点坐标，wi(x，y，z)表示第i口测井待插值点(x，y，z)上由径向基函数(10)式计算生成的局部空间插值距离权场。∈是径向参数，用来控制计算精度和平衡不同加权方案的关系。
其次，按照空间插值引导数据提供的整体构造导向特征，我们通过克里金插值法建立所有测井数据之间的外推关系，同时生成主要相控区域以外构造和相带相关待插值点(x，y，z)的空间插值距离权场，由径向基函数和克里金插值函数共同计算生成的空间插值距离权场如图8所示，用w(x，y，z表示。空间插值距离权场与井的位置有关联，主要用来考虑了井距离产生的影响。这是因为在做井外推时，距离井较远处的速度是不可靠的。使用空间插值距离权场能更好地控制距离井较远处的速度对模型融合的影响，从图8中可以看出该距离权场的分布与井的分布一致。
假设第i口测井速度数据用vi(x，y，z)表示，则全区N口测井速度插值填充模型按照如下公式计算生成：
图9为公式(11)计算获得测井速度插值填充模型，图中可见，测井速度插值填充模型能够很好地反映层位、断层、地震相和岩相的综合特征，垂向分辨率高，模型细节丰富，为后续多信息速度融合提供良好的数据基础。
进一步地，为测井速度插值填充模型P建立分类蒙板M，在M中，测井速度填充的区域用1指示，测井速度未填充的区域用0指示；按照测井速度插值填充模型分类蒙板M的区域划分，在测井速度插值填充模型P基础上，用地震速度场对M蒙板0值指示区域进行填充，同时，在M蒙板1值指示区域，通过加权最小二乘法对地震速度模型和测井速度插值填充模型进行加权融合，并用融合后的速度值对区域进行填充。通过如上两步填充，最终获得本发明提出的多信息融合模型R。
进一步地，获取的地震速度场，将测井速度插值填充模型和地震速度模型进行融合，可以获得相对中频的融合速度模型如图10所示和相对高频的融合速度模型如图11所示，图中Xline(trace)表示地震道，time表示时间，Merged velocity表示融合速度。图10所示的融合速度模型呈现中频，既反映了背景速度场的信息，还体现了测井速度的细节构造。图11所示融合速度模型呈现高频，主要体现了测井速度插值填充模型的特征，可以很好地反映细节。图10和图11所示的两个速度模型在空间的变化趋势与构造特征一致。
我们将多个速度模型融合问题表达成估计一个统一的速度模型，使得该速度模型同时拟合已知的多个速度模型的过程。本发明提供的融合方法不限制模型数量和维数，可以是任意多个速度模型，可以是一维、二维、三维或多维模型，融合方法可以根据参加融合的速度模型数量和维数进行相似扩展。同时，本发明提出用加权最小二乘法实施速度融合，可以是常规加权最小二乘法，可以是改进的加权最小二乘法，也可以是用神经网络实现的加权最小二乘法，三者技术效果相同，这里不做限定。本实施例中，我们设计了一个改进的加权最小二乘法，由两个一维速度模型v1(x)和v2(x)生成融合模型v(x)的过程用如下公式表达：
公式(12)中，假设v1(x)表示地震速度模型(图6)，v2(x)表示测井速度插值填充模型(图9)。w1(x)和w2(x)为空间变化的质量图，值域均为[0,1]，分别表征速度模型v1(x)和v2(x)在空间上的置信度。比如在v2(x)速度模型中，远离测井或者测井数据缺失的区域，速度的置信度是相对较低的，在该区域中w2(x)应设置为接近于0的值，在多个速度模型融合中给予较小的权重。因此w2(x)可定义为空间中某点x到最近测井距离的函数，距离越大，w2(x)越小。此外，在多个速度模型融合过程中，为了消除融合过程中产生的人为假象或者不连续性，我们引入一个正则化项D。在已知的两个速度模型v1和v2冲突(或者差异性)较大的区域，该正则化项有助于获得一个平滑过渡的最终的融合速度模型v(x)。最后，λ1,λ2,和∈均为0至1范围内的常数，用于平衡各方程。
令方程(12)可以简化为：
为了求解拟合以上公式(13)中多个方程的最小二乘解，我们求解其对应的如下法方程：
(LTL+∈2DTD)v＝LTd (14)
该方程的解可以写成如下形式：
v＝(LTL+∈2DTD)-1LTd (15)
为了在速度融合过程中更好的引入已有的地质认识，获得我们期待的融合速度模型，我们在以上求解过程中引入正则化项S：
S＝(I+∈2DTD)-1 (16)
该正则化项可以实现为一个通用的平滑算子，该算子用于有效引入我们的已有认识。比如，我们要估计的速度模型在断层两边应是不连续的，沿地层构造方向应具备一定的空间连续性。这时，我们可以将S实现为一个空变的平滑算子，使其用于增强和保持断层位置处速度模型的不连续性，并保证速度模型在沿地层构造方向的空间连续性。以上正则化项可重写为
∈2DTD＝S-1-I (17)
将其带入公式(15)中，可以得到以下基于正则化算子的求解形式：
v＝(LTL+S-1-I)-1LTd＝[I+S(LTL-I)]-1SLTd (18)
我们发现假设当S＝I时，v＝(LTL)-1LTd，此时融合速度为不加正则化的最小二乘解。当(LTL)＝I时，v＝SLTd，此时不需要反演。当S＝λI时，在λ→0的情况下，v≈λLTd。由于在实际问题中，L可能有物理单位，需要乘尺度因子1/λ，此时公式(18)变为：
v＝[λ2I+S(LTL-λ2I)]-1SLTd (19)
由于用共轭梯度法迭代反演时，通常需要正定且对称的算子，令S＝HHT，此时
v＝H[λ2I+HT(LTL-λ21)]-1HTLTd (20)
我们可以通过共轭梯度法对方程(20)进行反演，从而获得速度融合结果。
图6展示了地震速度模型，图中，Xline(trace)表示地震道，time表示时间，vint表示时间域的层速度。该速度模型是一种层速度模型，该模型浅部呈现低速特征，深部呈现高速特征，在断层左侧范围内呈现中高速特征，在断层右侧范围内呈现中低速特征。断层附近的速度差异较为明显。我们采取两种权重将图6所示地震速度模型和图9所示测井速度插值填充模型进行融合，得到了相对中频的融合速度模型(图10)和相对高频的融合速度模型(图11)。图10所示的融合速度模型呈现中频，既反映了背景速度场的信息，还体现了测井速度的细节构造。图11所示融合速度模型呈现高频，主要体现了测井速度插值填充模型的特征，可以很好地反映细节。图10和图11所示的两个速度模型在空间的变化趋势与构造特征一致。
我们还比较了两个融合速度模型与地震速度模型和测井速度插值填充模型在井位置处的差异。图12是本发明中测井位置处的地震速度、测井插值速度和相对中频融合速度的比较，图中横标Velocity表示速度，纵标Time表示时间，图13是本发明中测井位置处的地震速度、测井插值速度和相对高频融合速度的比较，图中横标Velocity表示速度，纵标Time表示时间，在图12和图13中我们可以发现：在井的上下两侧范围内，融合速度(最长的线)与地震速度(平滑线)保持吻合，这可以有效保证融合速度的稳定性。在测井曲线区域内，融合速度(最长的线)介于测井插值速度(波动最大的线)和地震速度(平滑线)之间。并且融合速度(最长的线)的整体趋势变化和地震速度(平滑线)保持一致。在细节部分，融合速度(最长的线)的波形变化和测井插值速度(波动最大的线)吻合较好。图12和图13提供的测井位置处的速度波形比较更好地印证了我们提出的融合方法的可靠性。
最后应说明的是：以上所述仅为本发明的优选实施例而已，并不用于限制本发明，尽管参照前述实施例对本发明进行了详细的说明，对于本领域的技术人员来说，其依然可以对前述各实施例所记载的技术方案进行修改，或者对其中部分技术特征进行等同替换，凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。
Claims (14)
1.一种多信息融合的速度建模方法，其特征在于，所述方法包括以下步骤，
获取地震偏移处理后的成像数据A；
获取岩相分类解释数据B，称为相控类型蒙板；
获取地震构造解释层位和断层数据C；
获取测井速度曲线数据D，并进行平滑处理；
获取地震速度场数据E；
通过边缘检测技术提取所述成像数据A的边缘特征，得到地震图像的构造骨架数据F；
将所述相控类型蒙板B与所述构造骨架数据F进行融合，获得相控构造导向数据G；
在所述相控构造导向数据G中加入地震构造解释数据C，形成空间插值引导数据H；
通过两步插值外推，获得与空间插值引导数据H相对应的空间插值距离权场Q和测井速度插值填充模型P；
为测井速度插值填充模型P建立分类蒙板M；
通过对类蒙板M进行两步填充，获得多信息融合模型R。
2.根据权利要求1所述的一种多信息融合的速度建模方法，其特征在于，所述成像数据A通过用户提供的叠前地震道集数据体和用户提供的偏移速度场获取。
3.根据权利要求1所述的一种多信息融合的速度建模方法，其特征在于，
所述岩相分类解释数据B根据成像数据A和用户提供的工区岩心录井柱状图/测井岩相图数据，通过测井岩相内插外推法、测井岩相神经网络智能预测法或地震数据地震相智能预测法获得。
4.根据权利要求1所述的一种多信息融合的速度建模方法，其特征在于，
所述地震构造解释层位和断层数据C是通过人工解释法或神经网络智能解释法对所述成像数据A进行解释获得的。
5.根据权利要求1所述的一种多信息融合的速度建模方法，其特征在于，
所述测井速度曲线数据D是通过用户提供的测井速度数据进行平滑处理得到的。
6.根据权利要求5所述的一种多信息融合的速度建模方法，其特征在于，
所述平滑处理包括：对测井速度数据进行中值滤波；对中值滤波所得数据做高斯递归滤波平滑。
7.根据权利要求6所述的一种多信息融合的速度建模方法，其特征在于，
所述中值滤波包括对一个含奇数个点的滑动窗口的数值进行排序，将中值赋值给中心点；所述中值滤波的滑动窗口形状包括线状、方形、圆形、十字形或用户提供的特殊几何形状。
8.根据权利要求6所述的一种多信息融合的速度建模方法，其特征在于，所述高斯递归滤波步骤包括：对中值滤波数据进行前向滤波；对前向滤波数据进行后向滤波；其中，
所述前向滤波公式为：
GSForward(n)＝a1*Input1(n)+a2*Input1(n-1)-a3*GSForward(n-1)-a4*GSForward(n-2)
式中，GSForward是前向滤波生成的数据，Input1是前向滤波输入数据，这里为中值滤波数据，n是数据样点位置，a1、a2、a3、a4是滤波器系数，根据不同工区实际需要进行设定；
所述后向滤波公式为：
GSBackward(n)＝b1*Input2(n)+b2*Input2(n+1)-b3*GSBackward(n+1)-b4*GSBackward(n+2)
式中，GSBackward是后向滤波生成的数据，Input2是后向滤波输入数据，这里为前向滤波数据，n是数据样点位置，b1、b2、b3、b4是滤波器系数，根据不同工区实际需要进行设定。
9.根据权利要求1所述的一种多信息融合的速度建模方法，其特征在于，所述地震速度场数据E是与所述地震偏移处理后的成像数据A对应的偏移速度场数据，通过用户提供的叠前地震道集数据体Data处理获得。
10.根据权利要求1所述的一种多信息融合的速度建模方法，其特征在于，所述通过两步插值外推获得与空间插值引导数据H相对应的空间插值距离权场Q和测井速度插值填充模型P包括，
在空间插值引导数据H中加入测井速度曲线数据D，判断测井数据空间和时间范围及与相控类型蒙板B最近邻相控区域的位置关系。
11.根据权利要求1或10所述的一种多信息融合的速度建模方法，其特征在于，
所述两步插值外推具体方法如下：
按照空间插值引导数据H提供的局部构造骨架特征，
通过径向基函数方法对落入相控区域的测井数据进行外推处理；
通过克里金插值法建立所有测井数据的外推关系，实现测井速度沿构造和相关相带的外推填充。
12.根据权利要求11所述的一种多信息融合的速度建模方法，其特征在于，所述径向基函数方法公式如下：
其中，ε表示径向基函数的尺度因子，这里取1；r表示插值点和井位置的欧式距离。
13.根据权利要求1所述的一种多信息融合的速度建模方法，其特征在于，所述分类蒙板M的区域划分为0值指示区和1值指示区域；
所述0值指示区域表示测井速度插值填充模型P的未填充区域；
所述1值指示区域表示通过加权最小二乘法对地震速度模型和测井速度插值填充模型进行加权融合后的速度值的填充区域。
14.一种多信息融合的速度建模方法的装置，其特征在于，所述装置包括，
信息收集装置用于获取地震偏移处理后的成像数据A、岩相分类解释数据B、地震构造解释层位和断层数据C、测井速度曲线数据D、获取地震速度场数据E；
信息处理装置用于通过边缘检测技术提取所述成像数据A的边缘特征，得到地震图像的构造骨架数据F，将岩相分类解释数据B与所述构造骨架数据F进行融合，获得相控构造导向数据G，在所述相控构造导向数据G中加入地震构造解释数据C，形成空间插值引导数据H；
信息建模装置用于通过两步插值外推，获得与空间插值引导数据H相对应的空间插值距离权场Q和测井速度插值填充模型P，为测井速度插值填充模型P建立分类蒙板M，通过对分类蒙板M进行两步填充，获得多信息融合模型R。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202210167182.6A CN116699678A (zh) | 2022-02-23 | 2022-02-23 | 一种多信息融合的速度建模方法及装置 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202210167182.6A CN116699678A (zh) | 2022-02-23 | 2022-02-23 | 一种多信息融合的速度建模方法及装置 |
Publications (1)
| Publication Number | Publication Date |
|---|---|
| CN116699678A true CN116699678A (zh) | 2023-09-05 |
Family
ID=87839754
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202210167182.6A Pending CN116699678A (zh) | 2022-02-23 | 2022-02-23 | 一种多信息融合的速度建模方法及装置 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN116699678A (zh) |
-
2022
- 2022-02-23 CN CN202210167182.6A patent/CN116699678A/zh active Pending
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| US11693139B2 (en) | Automated seismic interpretation-guided inversion | |
| EP0891562B1 (en) | 3-d geologic modelling | |
| US8363509B2 (en) | Method for building velocity models for pre-stack depth migration via the simultaneous joint inversion of seismic, gravity and magnetotelluric data | |
| US20150066460A1 (en) | Stratigraphic function | |
| US20100149917A1 (en) | Method For Geophysical and Geological Interpretation of Seismic Volumes In The Domains of Depth, Time, and Age | |
| CN108139499A (zh) | Q-补偿的全波场反演 | |
| Hamid et al. | Structurally constrained impedance inversion | |
| EP3639062A1 (en) | A method for validating geological model data over corresponding original seismic data | |
| EA011519B1 (ru) | Способы моделирования геологической среды и построения сейсмических изображений с использованием итерационного и избирательного обновления | |
| CN107462924B (zh) | 一种不依赖于测井资料的绝对波阻抗反演方法 | |
| Thiel et al. | 2D reservoir imaging using deep directional resistivity measurements | |
| CN113552625B (zh) | 一种用于常规陆域地震数据的多尺度全波形反演方法 | |
| WO2005026776A1 (en) | Wide-offset-range pre-stack depth migration method for seismic exploration | |
| CN109884710A (zh) | 针对激发井深设计的微测井层析成像方法 | |
| WO2022272058A1 (en) | Method and system for seismic imaging using s-wave velocity models and machine learning | |
| Pereira et al. | Iterative geostatistical seismic inversion incorporating local anisotropies | |
| EA030770B1 (ru) | Система и способ адаптивной сейсмической оптики | |
| Guo et al. | Becoming effective velocity-model builders and depth imagers, Part 2—The basics of velocity-model building, examples and discussions | |
| Cheng et al. | Multiscale fracture prediction technique via deep learning, seismic gradient disorder, and aberrance: Applied to tight sandstone reservoirs in the Hutubi block, southern Junggar Basin | |
| Li et al. | An efficient deep learning method for VSP wavefield separation: A DAS-VSP case | |
| Qu et al. | Using a synthetic data trained convolutional neural network for predicting subresolution thin layers from seismic data | |
| CN115880455A (zh) | 基于深度学习的三维智能插值方法 | |
| CN116699678A (zh) | 一种多信息融合的速度建模方法及装置 | |
| Li et al. | Automatic extraction of seismic data horizon across faults | |
| Fam et al. | Multi-Focusing stacking using the Very Fast Simulated Annealing global optimization algorithm |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| PB01 | Publication | ||
| PB01 | Publication | ||
| SE01 | Entry into force of request for substantive examination | ||
| SE01 | Entry into force of request for substantive examination |