CN110031895B - 一种基于图像缝合的多点地质统计学随机反演方法及装置 - Google Patents
一种基于图像缝合的多点地质统计学随机反演方法及装置 Download PDFInfo
- Publication number
- CN110031895B CN110031895B CN201910180942.5A CN201910180942A CN110031895B CN 110031895 B CN110031895 B CN 110031895B CN 201910180942 A CN201910180942 A CN 201910180942A CN 110031895 B CN110031895 B CN 110031895B
- Authority
- CN
- China
- Prior art keywords
- simulation
- inversion
- area
- data
- simulation result
- 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.)
- Active
Links
- 238000000034 method Methods 0.000 title claims abstract description 44
- 238000004088 simulation Methods 0.000 claims abstract description 110
- 238000000137 annealing Methods 0.000 claims abstract description 23
- 238000005457 optimization Methods 0.000 claims abstract description 18
- 238000012549 training Methods 0.000 claims description 32
- 230000008569 process Effects 0.000 claims description 13
- 230000002194 synthesizing effect Effects 0.000 claims description 7
- 239000002245 particle Substances 0.000 claims description 6
- 230000000704 physical effect Effects 0.000 claims description 6
- 238000005381 potential energy Methods 0.000 claims description 6
- XOFYZVNMUHMLCC-ZPOLXVRWSA-N prednisone Chemical compound O=C1C=C[C@]2(C)[C@H]3C(=O)C[C@](C)([C@@](CC4)(O)C(=O)CO)[C@@H]4[C@@H]3CCC2=C1 XOFYZVNMUHMLCC-ZPOLXVRWSA-N 0.000 claims description 6
- 238000006243 chemical reaction Methods 0.000 claims description 5
- 230000001186 cumulative effect Effects 0.000 claims description 3
- 238000007781 pre-processing Methods 0.000 claims description 3
- 238000004422 calculation algorithm Methods 0.000 abstract description 11
- 238000005516 engineering process Methods 0.000 abstract description 5
- 238000012512 characterization method Methods 0.000 abstract description 4
- 238000011160 research Methods 0.000 abstract description 4
- 238000013076 uncertainty analysis Methods 0.000 abstract description 2
- 238000012360 testing method Methods 0.000 description 12
- 238000004364 calculation method Methods 0.000 description 7
- 238000002948 stochastic simulation Methods 0.000 description 6
- 230000009471 action Effects 0.000 description 2
- 230000007547 defect Effects 0.000 description 2
- 230000000694 effects Effects 0.000 description 2
- 238000001914 filtration Methods 0.000 description 2
- 238000002679 ablation Methods 0.000 description 1
- 230000015572 biosynthetic process Effects 0.000 description 1
- 230000003750 conditioning effect Effects 0.000 description 1
- 230000008021 deposition Effects 0.000 description 1
- 238000011161 development Methods 0.000 description 1
- 238000010586 diagram Methods 0.000 description 1
- 238000000605 extraction Methods 0.000 description 1
- 238000007429 general method Methods 0.000 description 1
- 230000005484 gravity Effects 0.000 description 1
- 230000006872 improvement Effects 0.000 description 1
- 230000000452 restraining effect Effects 0.000 description 1
- 230000000717 retained effect Effects 0.000 description 1
- 238000004062 sedimentation Methods 0.000 description 1
- 239000013589 supplement Substances 0.000 description 1
- 238000003786 synthesis reaction Methods 0.000 description 1
Images
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
-
- 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/306—Analysis for determining physical properties of the subsurface, e.g. impedance, porosity or attenuation profiles
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
- G01V2210/00—Details of seismic processing or analysis
- G01V2210/50—Corrections or adjustments related to wave propagation
- G01V2210/51—Migration
- G01V2210/512—Pre-stack
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
- G01V2210/00—Details of seismic processing or analysis
- G01V2210/50—Corrections or adjustments related to wave propagation
- G01V2210/51—Migration
- G01V2210/514—Post-stack
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
- G01V2210/00—Details of seismic processing or analysis
- G01V2210/60—Analysis
- G01V2210/61—Analysis by combining or comparing a seismic data set with other data
- G01V2210/616—Data from specific type of measurement
- G01V2210/6169—Data from specific type of measurement using well-logging
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
- G01V2210/00—Details of seismic processing or analysis
- G01V2210/60—Analysis
- G01V2210/62—Physical property of subsurface
- G01V2210/624—Reservoir parameters
- G01V2210/6242—Elastic parameters, e.g. Young, Lamé or Poisson
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
- G01V2210/00—Details of seismic processing or analysis
- G01V2210/70—Other details related to processing
- G01V2210/74—Visualisation of seismic data
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)
- Management, Administration, Business Operations System, And Electronic Commerce (AREA)
- Geophysics And Detection Of Objects (AREA)
Abstract
本发明公开了一种基于图像缝合的多点地质统计学随机反演方法及装置,利用多点地质统计学进行随机模拟,能够更好的刻画目标体形态;引入图像缝合技术,能够快速高效地获得连续性较好的多点地质统计学随机模拟结果;通过正演方法将模拟结果合成地震记录与实际地震记录进行对比,先获得一个能与实际地震记录比较匹配的模拟结果;然后,通过量子退火优化算法,快速地将模拟结果优化,在地震记录的约束下获得最终的随机反演结果,提高了井间区域反演结果的精度和分辨率,多次进行随机反演可以获得多个不同的多点地质统计学随机反演结果,为不确定性分析提供可靠的依据,进而为储层表征和油藏描述等研究提供必要的参考信息。
Description
技术领域
本发明涉及随机反演、储层表征和建模技术领域,特别是一种基于图像缝合的多点地质统计学随机反演方法及装置。
背景技术
地震反演技术在储层表征和建模中能够发挥重要的作用。通常,地震反演可以分为确定性反演和随机反演。在实际采集的资料中,地震数据横向连续性较好,但受到由于地震资料的带限性质的影响,分辨率的提高受到限制;测井数据虽然分辨率高,但其仅能显示井孔位置的特征和参数变化,横向分辨率比较低。因此,地震反演的一项关键重点工作是充分地利用地震和测井数据中的有用信息提高横向和纵向分辨率,更精确地反映储层的特征。随机反演与确定性反演相比,随机反演方法可以充分利用井数据来提高分辨率,能够得到多个反演解,在分辨率提高和不确定性降低上具有一定的优势。
地质统计学是随机反演的一种重要原理。随机反演建立在随机模拟技术基础上。在地震反演中,两点地质统计学通过变差函数表征区域化变量的空间相关性,两点地质统计学方法计算效率快,能对离散变量和连续变量进行模拟和反演,但两点方法利用变差函数描述地质变量的相关性和变异性,仅仅通过某个方向上两点之间的地质变量的变化关系来描述空间的变化特性,对目标体形态的再现能力相对不强。多点地质统计学模拟是两点地质统计学的补充和发展。多点地质统计学方法以像元为基本模拟单元,并且结合基于目标的算法的优点,通过“训练图像”描述区域化变量的空间相关性,对目标体几何形态的刻画能力突出。目前,多点地质统计学随机模拟的主要算法包括单元正态方程模拟,模式模拟,滤波模拟,距离模拟等。其中单元正态方程模拟仅适用于离散变量(如岩相),模式模拟,滤波模拟等算法模拟结果连续性较差,且计算量大,耗时严重;除此之外,多数算法要求训练图像平稳,这与实际情况不符;鉴于上述原因,多点地质统计学随机模拟和随机反演的应用受到限制。随机反演通常的做法是先利用多点地质统计学模拟算法,产生大量的等可能的弹性参数模型;然后,正演合成地震记录与实际地震记录做匹配,若不能满足要求,则重新进行模拟,并合成地震记录与实际地震记录对比,需要进行大量的正演模拟,若直接将这些现有的多点地质统计学随机模拟算法用于随机反演中,则计算量及其大,尤其当实际数据量较大时,计算效率低,运行缓慢,阻碍了多点地质统计学随机反演的实际应用。
总之,目前多点地质统计学随机反演方法及研究存在以下问题和不足:1.现有的多点地质统计学随机模拟算法,模拟结果连续性差,需要进一步提高连续性;且现有模拟算法的计算效率低;不适合直接用于随机反演中;2.目前多点地质统计学随机反演方法在反演过程中进行大量的正演,效率低下;3.多点地质统计学随机反演受以上客观不足影响,适用范围局限,实际应用受到限制等问题。
发明内容
针对上述现有技术中存在的问题,本发明提供一种基于图像缝合的多点地质统计学随机反演方法,以及实现该方法的装置。
为了实现上述任务,本发明采用以下技术方案:
一种基于图像缝合的多点地质统计学随机反演方法,包括以下步骤:
步骤1,定义最大模拟次数、中间误差容忍度、终止条件、模版、重叠区大小以及切除大小、权重;
步骤2,对目标区域的测井曲线进行分析,根据实际需要计算井中弹性参数作为条件数据,划分网格并构建训练图像;
步骤3,利用定义的模板扫描训练图像,以测井数据为条件数据,从扫描结果中随机提取一个模式作为初始位置的模拟结果;
步骤4,根据设置的重叠区大小,提取已模拟区域边界的重叠部分,对重叠部分和训练图像进行卷积,以测井数据为条件数据,在训练图像中寻找与已模拟区域边界最接近的位置;
步骤5,访问已模拟区域的相邻区域,从已找到的与已模拟区域边界最接近的位置中,随机提取一个位置,按照定义的模版选定模式,切除多余部分后,缝合作为该区域的模拟结果;
步骤6,按照顺序访问下一个相邻区域,重复步骤4和步骤5,直到所有区域都完成模拟,则实现了一次基于图像缝合的多点地质统计学随机模拟,从而形成弹性参数模型;
步骤7,将所述弹性参数模型进行时深转换得到时间域的弹性参数模型,利用时间域弹性参数模型计算反射系数,然后通过反射系数根据褶积模型正演合成地震记录;
步骤8,依次利用所述合成地震记录与实际地震数据对比,判断是否满足预先设定的中间误差容忍度;若满足,则输出本次模拟结果;若不满足该容忍度,则判断是否达到最大模拟次数,若未达到,则重复执行步骤3至步骤8直到满足中间误差容忍度;若达到,则输出当前最优模拟结果。
进一步地,所述的方法还包括:
步骤9,对输出的最优模拟结果进行量子退火迭代优化,直到满足终止条件,获得最终的随机反演结果。
进一步地,步骤2所述的弹性参数包括波阻抗、弹性阻抗或纵横波速度和密度。
进一步地,步骤5所述切除多余部分时利用误差最小边界线,误差最小边界线的选择过程包括:
对于位置为(i,j)的像元,累计误差为:
Ee(i,j)=e(i,j)+min(Ee(i-1,j-1),Ee(i-1,j),Ee(i-1,j+1)),j=2,…,b-1
Ee(i,j)=e(i,j)+min(Ee(i-1,j),Ee(i-1,j+1)),j=1
Ee(i,j)=e(i,j)+min(Ee(i-1,j-1),Ee(i-1,j)),j=b
上式中,Ee(i-1,j-1)表示第(i-1,j-1)网格点处的累积最小误差;e(i,j)表示当前重叠区域Bd中第(i,j)网格点处的差异;假设与已模拟区域边界的重叠区域Bd误差最小的位置为Be,两者大小均为a×b;
对每一个i,得到累计误差Ee(i,j)最小时对应的j,连接获得的每一个(i,j),得到最小误差边界线。
进一步地,所述的通过反射系数根据褶积模型正演合成地震记录,表示为:
其中,b表示用于合成地震记录的地震子波,从地震记录中提取获得;R表示利用时间域的弹性参数模型计算的反射系数。
进一步地,步骤9所述的量子退火迭代优化的具体过程为:
1)输入弹性参数,定义初始动能CΓ(t)、退火计划以及终止条件;
其中,dj(j=1,2,…,M)为观测到的地震数据中的第j道数据,M表示地震道数,dj(m(k))表示利用弹性参数第k次的迭代更新值m(k)合成的地震数据,m表示输出的最优模拟结果,k表示量子退火中第k次迭代;
3)更新模型参数:
m(k+1)=m(k)+ξα
其中,ξ为0-1之间的随机数,α为步长;
4)计算势能ΔE,若ΔE<0,则用更新后的模型参数m(k+1)代替m(k),若ΔE>0,则根据概率接受准则判断是否用更新后的模型参数m(k+1)代替m(k);
所述的概率接受准则表达式为:
ρ(H)=(e-H/PΓ(t))P
其中,P表示粒子所处状态的粒子数,Γ(t)为横向场,在随机反演中为反演参数个数;H表示量子体系的总能量,表示为:
H=ΔE+CΓ(t)
其中,ΔE=E(m(k+1))-E(m(k))表示体系的势能,C为常数;
5)判断是否满足终止条件,若满足,则输出模型参数m(k+1)作为最终的物性参数反演结果,若不满足则返回步骤3)继续迭代直到满足终止条件。
在上述技术方案的基础上,本发明进一步公开了一种用于实现上述方法的装置,包括:
参数定义模块,用于定义最大模拟次数、中间误差容忍度、终止条件、模版、重叠区大小以及切除大小、权重;
数据预处理模块,用于对目标区域的测井曲线进行分析,根据实际需要计算井中弹性参数作为条件数据,划分网格并构建训练图像;
扫描模块,用于利用定义的模板扫描训练图像,从扫描结果中随机提取一个模式作为初始位置的模拟结果;
寻找最接近位置模块,用于根据设置的重叠区大小,提取已模拟区域边界的重叠部分,对重叠部分和训练图像进行卷积,以测井数据为条件数据,在训练图像中寻找与已模拟区域边界最接近的位置;
缝合模块,用于访问已模拟区域的相邻区域,从已找到的与已模拟区域边界最接近的位置中,随机提取一个位置,按照定义的模版选定模式,利用误差最小边界线切除多余部分后,缝合作为该区域的模拟结果;
随机模拟模块,用于按照顺序访问下一个相邻区域,重复寻找最接近位置模块和区域模拟结果模块中的执行步骤,直到所有区域都完成模拟,则实现了一次基于图像缝合的多点地质统计学随机模拟,从而形成弹性参数模型;
地震正演模块,用于将所述弹性参数模型进行时深转换得到时间域的弹性参数模型,利用时间域弹性参数模型计算反射系数,然后通过反射系数根据褶积模型正演合成地震记录;
对比匹配模块,用于依次利用所述合成地震记录与实际地震数据对比,判断是否满足预先设定的中间误差容忍度;若满足,则输出本次模拟结果;若不满足该容忍度,则判断是否达到最大模拟次数,若未达到,则重复执行图像缝合模拟模块至对比匹配模块之间所有模块中的执行步骤,直到满足中间误差容忍度;若达到,则输出当前最优模拟结果。
进一步地,所述的装置还包括:
优化模块,用于对输出的最优模拟结果进行量子退火迭代优化,直到满足终止条件,获得最终的随机反演结果。
本发明基于图像缝合实现多点地质统计学随机模拟,并结合量子退火优化算法形成了一种有效的多点地质统计学随机反演方法。与现有技术相比,本发明具有以下技术特点:
1.本发明应用多点地质统计学随机模拟,考虑空间中多点之间的相关性,更好的表征目标的形态等特征。
2.本发明引入图像缝合技术,克服了目前多点地质统计学随机模拟方法计算慢、连续性差的不足,也提高了多点地质统计学随机反演的计算效率和实用性。
3.本发明结合量子退火优化算法,在模拟结果的基础上进行优化,能够减小搜索空间,提高优化的速度。
4.本发明提供的随机反演方法和装置可以较为快速精确地获得反演结果,分辨率高于确定性反演,计算效率高于现有的多点地质统计学随机反演方法。
5.在不需要改变反演流程的前提下,利用叠前地震数据和正演算子很容易即可将随机反演从叠后推广到叠前。
附图说明
图1是本发明方法的流程示意图;
图2是本发明实施例用于测试的模型数据的波阻抗剖面;
图3是本发明实施例用于多点地质统计学随机模拟的井中波阻抗曲线,从图2的模型中抽取获得,其分布位置如图2中黑色直线所示;
图4是本发明实施例中测试模型对应的实际地震记录;
图5是本发明实施例中测试模型波阻抗反演结果;
图6是本发明实施例中测试模型波阻抗反演结果与测试模型之间的误差。
具体实施方式
为使本发明实施例的目的、技术方案和优点更加清楚明白,下面结合附图对本发明实施例做进一步详细说明。在此,本发明的示意性实施例及其说明用于解释本发明,但并不作为对本发明的限定。
本发明提出了一种基于图像缝合的多点地质统计学随机反演方法,通过正演方法将模拟结果合成地震记录与实际地震记录进行对比,先获得一个能与实际地震记录比较匹配的模拟结果,然后,通过量子退火优化算法,快速地将模拟结果优化,在地震记录的约束下获得最终的随机反演结果,提高了井间区域反演结果的精度和分辨率,多次进行随机反演可以获得多个不同的多点地质统计学随机反演结果,为不确定性分析提供可靠的依据,进而为储层表征和油藏描述等研究提供必要的参考信息。如图1所示,本发明具体包括以下步骤:
步骤1,参数定义
定义最大模拟次数、中间误差容忍度、模版、重叠区大小、权重。
本发明引入图像缝合技术实现多点地质统计学随机模拟,以地震数据为约束,最终实现随机反演。首先,定义模拟参数,反演参数和图像缝合参数,其中:
模版是多点地质统计学中用来从训练图像中获取模式的工具,由空间中几个点之间的结构组成;定义模板包括定义模板的大小及结构,即模板的空间形态。
中间误差容忍度用于从随机模拟结果中挑选出一个比较接近地震数据的结果,利用该结果约束搜索空间,作为于量子退火优化的初始模型。
最大模拟次数设置最大所允许的模拟次数。
重叠区大小和切除大小是在图像缝合过程中保证连续性的参数。
设置的权重用于在图像缝合中考虑最相似位置和条件数据各自的比重,取值在0~1之间。
步骤2,数据预处理
对目标区域的测井曲线进行分析,根据实际需要计算井中弹性参数作为条件数据,划分网格并构建训练图像;其中,所述的弹性参数包括波阻抗、弹性阻抗或纵横波速度和密度。
随机反演以测井数据作为硬数据,能够忠实于测井数据,因此根据实际反演任务的要求,对获取的测井曲线进行分析,计算获得需要反演的弹性参数,作为随机反演的硬数据,其中弹性参数可以是波阻抗、弹性阻抗或者纵横波速度和密度;针对反演目标,划分网格,形成基本模拟或反演单元;利用工区的前期研究资料、地质认识和测井等数据,构建训练图像;训练图像是多点地质统计学中用于表述空间中多点之间相关性的基本工具,是可以表述实际储层的结构、几何形态及其分布模式的数字化图像,通常被看作是先验的地质概念,主要包含了待模拟区域的主要的特征模式,不要求与精准的测井信息一致或者服从一定的分布。训练图像可以反映地质体空间分布的一般性特征,是特征模式定量化的表示,可以体现出不同特征模式的空间关系,在多点地质统计学随机模拟和随机反演中具有重要的作用。
步骤3,图像扫描
利用定义的模板扫描训练图像,以测井数据为条件数据,从扫描结果中随机提取一个模式作为初始位置的模拟结果。
利用步骤1定义的数据模板对训练图像进行扫描,获得训练图像中所包含的地质模拟,每一个地质模式都可以看作是一幅子图像;从所有子图像中随机抽取一幅作为多点地质统计学随机模拟初始位置的模拟结果。
步骤4,寻找最接近位置
根据设置的重叠区大小,提取已模拟区域边界的重叠部分,对重叠部分和训练图像进行卷积,以测井数据为条件数据,在训练图像中寻找与已模拟区域边界最接近的位置。
在进行位置模拟时,是分块进行的,每次模拟只在整个需要模拟的区域中完成部分区域的模拟,这部分的大小等于模版的大小,因此每次模拟之后会有一个边界。
该步骤中,根据定义的重叠区大小,从已模拟区域的边界提取重叠区域Bd,利用Bd在训练图像T上按顺序滑动,对应位置的值相乘,然后求和,实现边界区域Bd与训练图像T通过在训练图像T中不同位置处求取的卷积结果;具体的卷积操作如下:
其中,n表示重叠区域中模拟单元的位置,Bd中共有n的网格点,Bd(i)表示第i个网格节点中的值,Tm表示重叠区域Bd在T中滑动的某一位置,大小与Bd相同;从而可以得到T中每一位置与当前重叠区域Bd之间的差异Εd大小,则根据Εd的值即可判断是否与已模拟区域边界是否最接近,Εd越小说明越接近,从而找到与模拟区域边界最接近的位置。
然而,有可能存在多个Εd值相等且最小的情况,因此在下一个步骤中要从中随机提取一个位置,若只有一个Εd最小,则不需要随机提取。
步骤5,图像缝合
访问已模拟区域的相邻区域,从已找到的与已模拟区域边界最接近的位置中,随机提取一个位置,按照定义的模版选定模式,利用误差最小边界线切除多余部分后,缝合作为该区域的模拟结果。
按照模版大小访问已模拟区域的相邻区域,若该分布区域中有井分布,考虑井数据的约束作用,为了使模拟结果能够忠实于硬数据,计算上述步骤中扫描获得的每一个模式中在井点位置处与井数据之间的误差,若有多个井分布在该区域,则将绝对误差相加,得到Ec,然后计算总体差异:
E=(1-w)Ed/D+wEc/C (2)
其中,w表示步骤1中定义的权重,D表示重叠区域的大小,C表示当前模拟区域中分布的硬数据的个数,Εd表示训练图像T中每个位置Tm处与Bd的差异值。选择E最小的位置,按照模版选择模式,利用误差最小边界线原理切除选定模式的多余部分,将切除后的模式放置在当前模拟区域;若有多个位置的E最小,则从中随机选取其中一个。
假设与已模拟区域边界的重叠区域Bd误差最小的位置为Be,两者大小均为a×b,e=(Bd-Be)2,则最小误差边界线选择过程如下:对于位置为(i,j)的像元,累计误差为:
Ee(i,j)=e(i,j)+min(Ee(i-1,j-1),Ee(i-1,j),Ee(i-1,j+1)),j=2,…,b-1
Ee(i,j)=e(i,j)+min(Ee(i-1,j),Ee(i-1,j+1)),j=1 (3)
Ee(i,j)=e(i,j)+min(Ee(i-1,j-1),Ee(i-1,j)),j=b
在式(3)中,参数的下标表示网格点的位置,例如Ee(i-1,j-1)表示第(i-1,j-1)网格点处的累积最小误差;e(i,j)表示当前Bd中第(i,j)网格点处的差异。
对每一个i,得到累计误差Ee(i,j)最小时对应的j,连接获得的每一个(i,j),得到最小误差边界线。
步骤6,随机模拟
按照顺序访问下一个相邻区域,重复步骤4和步骤5,直到所有区域都完成模拟,则实现了一次基于图像缝合的多点地质统计学随机模拟,从而形成弹性参数模型;多次模拟则可形成多个模拟的弹性参数模型。
步骤7,地震正演
将所述弹性参数模型进行时深转换得到时间域的弹性参数模型,利用时间域弹性参数模型计算反射系数,然后通过反射系数根据褶积模型正演合成地震记录。
前述模拟生成的数据是深度域的数据,而地震记录是时间域的数据,因此在合成地震数据之间,要将深度域数据转化为时间域数据,得到时间域的弹性参数数据之后,就可以进行正演获得合成地震记录dsyn,便于与实际地震记录dobs对比,从而选择出能与实际地震记录最匹配的岩相分布和物性参数分布,实现岩相和物性参数的随机反演。
对于叠后波阻抗随机反演,利用时间域的弹性参数模型计算反射系数R,然后利用褶积模型进行正演合成地震记录:
其中,b表示用于合成地震记录的地震子波,从地震记录中提取获得。每一个弹性参数模型都可以合成一个地震记录。
步骤8,对比匹配
依次利用所述合成地震记录与实际地震数据对比,判断是否满足预先设定的中间误差容忍度;若满足,则输出本次模拟结果;若不满足该容忍度,则判断是否达到最大模拟次数,若未达到,则重复执行步骤3至步骤8直到满足中间误差容忍度;若达到,则输出当前最优模拟结果。
合成地震记录与实际地震记录的对比匹配,可以通过计算两者之间的误差实现,叠后波阻抗反演使用的目标函数(合成地震记录与实际地震记录的误差)可以表示为:
O=||dobs-dsyn||2 (5)
若某一次模拟结果的合成地震记录小于或等于中间误差容忍度,则终止模拟,输出对应的多点地质统计学模拟数据;若大于中间误差容忍度,则判断是否达到最大模拟次数,若未达到最大模拟次数,返回步骤3,重复执行步骤3至步骤8直到满足中间误差容忍度;若达到最大模拟次数,则终止模拟,输出当前最优多点地质统计学模拟数据(即当前O最小时对应的模拟数据)。
步骤9,结果优化
对输出的最优模拟结果进行量子退火迭代优化,直到满足终止条件,获得最终的随机反演结果。
由于在多点地质统计学模拟过程中,井间区域的约束性比较差,而地震反演通常又具有多解性。因此将该输出的数据进行优化,使其能与地震数据尽可能的一致,从而提高反演结果的精度。量子退火在退火收敛速度和避免陷入局部极小方面有一定优势,对于提高随机反演的精度和效率都有重要价值。
在地球物理反演中,量子退火优化可以表示为:
H=ΔE+CΓ(t) (6)
其中,ΔE=E(m(k+1))-E(m(k))表示体系的势能,m表示上述输出的最优模拟结果,即在模拟过程中产生的能够使O最小时对应的模拟弹性参数模型,k表示量子退火中第k次迭代;第k次迭代的弹性参数为 为目标函数,dj(j=1,2,…,M)为观测到的地震数据中的第j道数据,M表示地震道数,dj(m(k))表示利用弹性参数第k次的迭代更新值m(k)合成的地震数据(地震记录)的第j道数据;CΓ(t)表示体系的动能,C为常数。
在每一次的搜索过程中,依然通过概率接收准则保留概率较大状态对应的参数,量子退火的概率接受准则由体系在横向场Γ(t)作用下所处某状态的机率给出:
ρ(H)=(e-H/PΓ(t))P (7)
其中,P表示粒子所处状态的粒子数,在随机反演中为反演参数个数;H表示量子体系的总能量,通过式6计算。
因此,在获得步骤8输出的最优模拟结果,即弹性参数后,执行量子退火优化的步骤可以概括为:
1)输入弹性参数,定义初始动能CΓ(t)、退火计划以及终止条件;
2)计算目标函数值E(m(k));
3)更新模型参数:
m(k+1)=m(k)+ξα (8)
其中,ξ为0-1之间的随机数,α为步长;
4)计算势能ΔE,若ΔE<0,则用更新后的模型参数m(k+1)代替m(k),若ΔE>0,则根据概率接受准则(式7)判断是否用更新后的模型参数m(k+1)代替m(k);
5)判断是否满足终止条件,若满足,则输出模型参数m(k+1)作为最终的物性参数反演结果,若不满足则返回步骤3)继续迭代直到满足终止条件。所述的终止条件可以为达到最大迭代次数、满足误差要求、误差连续保持不变的最大迭代次数等。
本发明利用多点地质统计学进行随机模拟,考虑空间中多点之间的相关性,能够更好的刻画目标体形态;引入图像缝合技术,将多点地质统计学随机模拟过程看作图像缝合的过程,能够快速高效地获得连续性较好的多点地质统计学随机模拟结果,计算快,模拟效果好,实用性强。
如图2至图4所示,图2是本发明一实施例用于测试的模型数据的波阻抗剖面。从图中可以看出,该模型是一个河道沉积模型,发育多期河道,河道间为泛滥平原沉积,河道和泛滥平原沉积的岩性不同,波阻抗有差异,深色为河道,浅色为泛滥平原。图3是本发明一实施例用于多点地质统计学随机模拟的井中波阻抗曲线,从图2的模型中抽取获得,用做条件数据,其位置如图2中黑色直线所示;图4是本发明一实施例中测试模型对应的实际地震记录,该地震记录由波阻抗计算出的反射系数与设置的地震子波褶积计算得到,从地震记录中可以大致分辨出河道的位置。
如图5和图6所示,图5是本发明一实施例中测试模型波阻抗反演结果,通过本发明方法反演得到,从图中可以看出,反演结果很好的反演出了河道和泛滥平原的波阻抗特征,有效区分河道和泛滥平原的波阻抗差异,河道波阻抗值较小,而泛滥平原波阻抗值较大,与模型一致,说明测试有效,体现了本发明方法的可行性和有效性。
图6是本发明一实施例中测试模型波阻抗反演结果与测试模型之间的误差。从误差图中可以看出,误差集中在0附近,表明反演结果与测试模型比较接近,充分验证了本发明方法的效果。
Claims (2)
1.一种基于图像缝合的多点地质统计学随机反演方法,其特征在于,包括以下步骤:
步骤1,定义最大模拟次数、中间误差容忍度、模版、重叠区大小以及权重;
步骤2,对目标区域的测井曲线进行分析,根据实际需要计算井中弹性参数作为条件数据,划分网格并构建训练图像;
步骤3,利用定义的模板扫描训练图像,以测井数据为条件数据,从扫描结果中随机提取一个模式作为初始位置的模拟结果;
步骤4,根据设置的重叠区大小,提取已模拟区域边界的重叠部分,对重叠部分和训练图像进行卷积,以测井数据为条件数据,在训练图像中寻找与已模拟区域边界最接近的位置;
步骤5,访问已模拟区域的相邻区域,从已找到的与已模拟区域边界最接近的位置中,随机提取一个位置,按照定义的模版选定模式,切除多余部分后,缝合作为该区域的模拟结果;
步骤6,按照顺序访问下一个相邻区域,重复步骤4和步骤5,直到所有区域都完成模拟,则实现了一次基于图像缝合的多点地质统计学随机模拟,从而形成弹性参数模型;
步骤7,将所述弹性参数模型进行时深转换得到时间域的弹性参数模型,利用时间域弹性参数模型计算反射系数,然后通过反射系数根据褶积模型正演合成地震记录;
步骤8,依次利用所述合成地震记录与实际地震数据对比,判断是否满足预先设定的中间误差容忍度;若满足,则输出本次模拟结果;若不满足该容忍度,则判断是否达到最大模拟次数,若未达到,则重复执行步骤3至步骤8直到满足中间误差容忍度;若达到,则输出当前最优模拟结果;
还包括:
步骤9,对输出的最优模拟结果进行量子退火迭代优化,直到满足终止条件,获得最终的随机反演结果;
步骤2所述的弹性参数包括波阻抗、弹性阻抗或纵横波速度和密度;
步骤5所述切除多余部分时利用误差最小边界线,误差最小边界线的选择过程包括:
对于位置为(i,j)的像元,累计误差为:
Ee(i,j)=e(i,j)+min(Ee(i-1,j-1),Ee(i-1,j),Ee(i-1,j+1)),j=2,…,b-1
Ee(i,j)=e(i,j)+min(Ee(i-1,j),Ee(i-1,j+1)),j=1
Ee(i,j)=e(i,j)+min(Ee(i-1,j-1),Ee(i-1,j)),j=b
上式中,Ee(i-1,j-1)表示第(i-1,j-1)网格点处的累积最小误差;e(i,j)表示当前重叠区域Bd中第(i,j)网格点处的差异;假设与已模拟区域边界的重叠区域Bd误差最小的位置为Be,两者大小均为a×b;
对每一个i,得到累计误差Ee(i,j)最小时对应的j,连接获得的每一个(i,j),得到最小误差边界线;
所述的通过反射系数根据褶积模型正演合成地震记录,表示为:
其中,b表示用于合成地震记录的地震子波,从地震记录中提取获得;R表示利用时间域的弹性参数模型计算的反射系数;
步骤9所述的量子退火迭代优化的具体过程为:
1)输入弹性参数,定义初始动能CΓ(t)、退火计划以及终止条件;
其中,dj(j=1,2,…,M)为观测到的地震数据中的第j道数据,M表示地震道数,dj(m(k))表示利用弹性参数第k次的迭代更新值m(k)合成的地震数据,m表示输出的最优模拟结果,k表示量子退火中第k次迭代;
3)更新模型参数:
m(k+1)=m(k)+ξα
其中,ξ为0-1之间的随机数,α为步长;
4)计算势能ΔE,若ΔE<0,则用更新后的模型参数m(k+1)代替m(k),若ΔE>0,则根据概率接受准则判断是否用更新后的模型参数m(k+1)代替m(k);
所述的概率接受准则表达式为:
ρ(H)=(e-H/PΓ(t))P
其中,P表示粒子所处状态的粒子数,在随机反演中为反演参数个数;H表示量子体系的总能量,表示为:
H=ΔE+CΓ(t)
其中,ΔE=E(m(k+1))-E(m(k))表示体系的势能,C为常数;
5)判断是否满足终止条件,若满足,则输出模型参数m(k+1)作为最终的物性参数反演结果,若不满足则返回步骤3)继续迭代直到满足终止条件。
2.一种用于实现如权利要求1所述方法的装置,包括:
参数定义模块,用于定义最大模拟次数、中间误差容忍度、终止条件、模版、重叠区大小以及切除大小、权重;
数据预处理模块,用于对目标区域的测井曲线进行分析,根据实际需要计算井中弹性参数作为条件数据,划分网格并构建训练图像;
扫描模块,用于利用定义的模板扫描训练图像,从扫描结果中随机提取一个模式作为初始位置的模拟结果;
寻找最接近位置模块,用于根据设置的重叠区大小,提取已模拟区域边界的重叠部分,对重叠部分和训练图像进行卷积,以测井数据为条件数据,在训练图像中寻找与已模拟区域边界最接近的位置;
缝合模块,用于访问已模拟区域的相邻区域,从已找到的与已模拟区域边界最接近的位置中,随机提取一个位置,按照定义的模版选定模式,利用误差最小边界线切除多余部分后,缝合作为该区域的模拟结果;
随机模拟模块,用于按照顺序访问下一个相邻区域,重复寻找最接近位置模块和区域模拟结果模块中的执行步骤,直到所有区域都完成模拟,则实现了一次基于图像缝合的多点地质统计学随机模拟,从而形成弹性参数模型;
地震正演模块,用于将所述弹性参数模型进行时深转换得到时间域的弹性参数模型,利用时间域弹性参数模型计算反射系数,然后通过反射系数根据褶积模型正演合成地震记录;
对比匹配模块,用于依次利用所述合成地震记录与实际地震数据对比,判断是否满足预先设定的中间误差容忍度;若满足,则输出本次模拟结果;若不满足该容忍度,则判断是否达到最大模拟次数,若未达到,则重复执行图像缝合模拟模块至对比匹配模块之间所有模块中的执行步骤,直到满足中间误差容忍度;若达到,则输出当前最优模拟结果;
还包括:
优化模块,用于对输出的最优模拟结果进行量子退火迭代优化,直到满足终止条件,获得最终的随机反演结果。
Priority Applications (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN201910180942.5A CN110031895B (zh) | 2019-03-11 | 2019-03-11 | 一种基于图像缝合的多点地质统计学随机反演方法及装置 |
Applications Claiming Priority (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN201910180942.5A CN110031895B (zh) | 2019-03-11 | 2019-03-11 | 一种基于图像缝合的多点地质统计学随机反演方法及装置 |
Publications (2)
Publication Number | Publication Date |
---|---|
CN110031895A CN110031895A (zh) | 2019-07-19 |
CN110031895B true CN110031895B (zh) | 2020-12-15 |
Family
ID=67235177
Family Applications (1)
Application Number | Title | Priority Date | Filing Date |
---|---|---|---|
CN201910180942.5A Active CN110031895B (zh) | 2019-03-11 | 2019-03-11 | 一种基于图像缝合的多点地质统计学随机反演方法及装置 |
Country Status (1)
Country | Link |
---|---|
CN (1) | CN110031895B (zh) |
Families Citing this family (5)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN111505713B (zh) * | 2020-01-21 | 2021-05-07 | 长江大学 | 基于多点地质统计的叠前地震反演方法 |
CN111273348B (zh) | 2020-01-21 | 2021-02-05 | 长江大学 | 基于更新概率比率恒定理论的多点地质统计叠前反演方法 |
CN111709169B (zh) * | 2020-05-29 | 2021-08-24 | 中国地质大学(武汉) | 基于条件传导概率的多点地质统计学随机模拟方法 |
CN113189677B (zh) * | 2021-04-22 | 2022-05-31 | 西南石油大学 | 一种三维油藏物性参数模型自动更新方法 |
CN115060769B (zh) * | 2022-06-07 | 2024-04-02 | 深圳大学 | 一种基于智能反演的隧道围岩裂隙及松动检测方法、系统 |
Citations (4)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
US8116979B2 (en) * | 2005-03-09 | 2012-02-14 | Baker Hughes Incorporated | System and method for determining a more accurate resistivity model of a geological formation using time-lapse well logging data |
CN103533266A (zh) * | 2013-10-01 | 2014-01-22 | 中国人民解放军国防科学技术大学 | 垂直方向宽视域的360度拼接式全景摄像机 |
CN104199092A (zh) * | 2014-08-31 | 2014-12-10 | 电子科技大学 | 基于多层次框架的三维全层位自动追踪方法 |
CN108645994A (zh) * | 2018-04-25 | 2018-10-12 | 中国石油大学(北京) | 一种基于多点地质统计学的地质随机反演方法及装置 |
Family Cites Families (1)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN102338885B (zh) * | 2011-06-20 | 2016-05-25 | 中国海洋石油总公司 | 三分量vsp资料初至时间自动拾取方法 |
-
2019
- 2019-03-11 CN CN201910180942.5A patent/CN110031895B/zh active Active
Patent Citations (4)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
US8116979B2 (en) * | 2005-03-09 | 2012-02-14 | Baker Hughes Incorporated | System and method for determining a more accurate resistivity model of a geological formation using time-lapse well logging data |
CN103533266A (zh) * | 2013-10-01 | 2014-01-22 | 中国人民解放军国防科学技术大学 | 垂直方向宽视域的360度拼接式全景摄像机 |
CN104199092A (zh) * | 2014-08-31 | 2014-12-10 | 电子科技大学 | 基于多层次框架的三维全层位自动追踪方法 |
CN108645994A (zh) * | 2018-04-25 | 2018-10-12 | 中国石油大学(北京) | 一种基于多点地质统计学的地质随机反演方法及装置 |
Non-Patent Citations (1)
Title |
---|
基于改进图像缝合技术的纹理合成研究;许刚,等;《计算机工程与设计》;20131130;第34卷(第11期);3948-3951、3995 * |
Also Published As
Publication number | Publication date |
---|---|
CN110031895A (zh) | 2019-07-19 |
Similar Documents
Publication | Publication Date | Title |
---|---|---|
CN110031895B (zh) | 一种基于图像缝合的多点地质统计学随机反演方法及装置 | |
CN109709603B (zh) | 地震层位识别与追踪方法、系统 | |
CN110031896B (zh) | 基于多点地质统计学先验信息的地震随机反演方法及装置 | |
WO2021147529A1 (zh) | 基于更新概率比率恒定理论的多点地质统计叠前反演方法 | |
CN108645994B (zh) | 一种基于多点地质统计学的地质随机反演方法及装置 | |
de Vries et al. | Application of multiple point geostatistics to non-stationary images | |
Azevedo et al. | Generative adversarial network as a stochastic subsurface model reconstruction | |
WO2017007924A1 (en) | Improved geobody continuity in geological models based on multiple point statistics | |
CA3122488A1 (en) | Subsurface models with uncertainty quantification | |
AU2020397825B2 (en) | Generation of subsurface representations using layer-space | |
CN111080021B (zh) | 一种基于地质信息库的砂体构型cmm神经网络预测方法 | |
US10984590B1 (en) | Generation of subsurface representations using layer-space | |
CN111596978A (zh) | 用人工智能进行岩相分类的网页显示方法、模块和系统 | |
KR20160024232A (ko) | 지반정보에 기초한 3차원 공간 모델링 방법 | |
US10890688B2 (en) | Method for generating secondary data in geostatistics using observed data | |
CN111505713B (zh) | 基于多点地质统计的叠前地震反演方法 | |
Pyrcz | Integration of geologic information into geostatistical models | |
CN104422969A (zh) | 一种减小电磁测深反演结果非唯一性的方法 | |
CN110927793A (zh) | 一种基于序贯随机模糊模拟的储层预测方法及系统 | |
Yanshu et al. | A three-dimensional model of deep-water turbidity channel in Plutonio oilfield, Angola: From training image generation, optimization to multi-point geostatistical modelling | |
Koneshloo et al. | A workflow for static reservoir modeling guided by seismic data in a fluvial system | |
EP3320450B1 (en) | Improved geobody continuity in geological models based on multiple point statistics | |
Jo et al. | Conditioning stratigraphic, rule-Based models with generative adversarial networks: A deepwater lobe, deep learning example | |
Exterkoetter et al. | Petroleum reservoir connectivity patterns reconstruction using deep convolutional generative adversarial networks | |
Ferreira | Rule-based modeling of turbiditic channel environments constrained with 3D and 4D seismic and well data |
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 | ||
GR01 | Patent grant | ||
GR01 | Patent grant |