CN114839675A - 一种建立三维速度模型的方法 - Google Patents
一种建立三维速度模型的方法 Download PDFInfo
- Publication number
- CN114839675A CN114839675A CN202110132180.9A CN202110132180A CN114839675A CN 114839675 A CN114839675 A CN 114839675A CN 202110132180 A CN202110132180 A CN 202110132180A CN 114839675 A CN114839675 A CN 114839675A
- Authority
- CN
- China
- Prior art keywords
- shot
- point
- ray
- model
- travel time
- 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.)
- Granted
Links
- 238000000034 method Methods 0.000 title claims abstract description 34
- 238000001514 detection method Methods 0.000 claims abstract description 90
- 238000012937 correction Methods 0.000 claims abstract description 10
- 238000007781 pre-processing Methods 0.000 claims abstract description 4
- 239000000523 sample Substances 0.000 claims description 29
- 239000013598 vector Substances 0.000 claims description 13
- 238000004364 calculation method Methods 0.000 claims description 8
- 230000008569 process Effects 0.000 claims description 7
- 238000012360 testing method Methods 0.000 claims description 5
- 239000006185 dispersion Substances 0.000 claims description 4
- 230000009466 transformation Effects 0.000 claims description 4
- 230000001902 propagating effect Effects 0.000 claims description 3
- 238000012545 processing Methods 0.000 abstract description 2
- 230000000644 propagated effect Effects 0.000 abstract 1
- 238000003384 imaging method Methods 0.000 description 10
- 238000003325 tomography Methods 0.000 description 8
- 238000013508 migration Methods 0.000 description 5
- 230000005012 migration Effects 0.000 description 5
- 230000009286 beneficial effect Effects 0.000 description 2
- 238000010586 diagram Methods 0.000 description 2
- 230000000694 effects Effects 0.000 description 2
- 238000005516 engineering process Methods 0.000 description 2
- 238000004458 analytical method Methods 0.000 description 1
- 230000008859 change Effects 0.000 description 1
- 230000001427 coherent effect Effects 0.000 description 1
- 238000010276 construction Methods 0.000 description 1
- 230000001419 dependent effect Effects 0.000 description 1
- 235000013399 edible fruits Nutrition 0.000 description 1
- 238000007689 inspection Methods 0.000 description 1
- 239000011159 matrix material Substances 0.000 description 1
- 238000012986 modification Methods 0.000 description 1
- 230000004048 modification Effects 0.000 description 1
- 230000003068 static effect Effects 0.000 description 1
- 230000001629 suppression Effects 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
- G01V1/30—Analysis
- G01V1/303—Analysis for determining velocity profiles or travel times
-
- 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
- G01V1/305—Travel times
-
- 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/622—Velocity, density or impedance
- G01V2210/6222—Velocity; travel time
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
本发明属于地震资料数据处理技术领域,提供了一种建立三维速度模型的方法,包括:对采集的原始三维地震资料进行预处理;在共炮点道集和共检波点道集拾取局部相关反射波走时和斜率数据,形成观测数据空间;建立离散速度模型,形成模型空间,并给定初始离散速度模型;在当前速度模型中,用拾取的斜率数据分别计算从炮点和检波点向地下传播的初始射线方向,通过求解三维程函方程获得分别从炮点与检波点的射线路径、走时和对应的射线参数;建立反演目标函数及反演方程,求解获得对速度模型的修正量,更新速度模型;判断当前速度模型是否符合要求,若是,则输出当前速度模型,本发明效率高、对复杂地质构造和低信噪比资料有很强适应性。
Description
技术领域
本发明属于地震资料数据处理技术领域,特别涉及一种建立三维速度模型的方法。
背景技术
地震勘探资料叠前深度偏移成像效果极大地依赖于所用的地下速度模型的准确度。地震走时层析反演是一种利用地震波走时资料建立地下低频速度模型的有效技术。常规反射波走时层析反演需要连续追踪各个同相轴,并拾取各个同相轴的走时数据,特别是需要事先建立这些同相轴与地下地层界面之间的对应关系,这对于复杂地质构造地区或低信噪比地震资料,往往难以实现。此外,反射波走时层析反演中,需要通过射线追踪计算当前速度模型的地震波走时以及反演方程中的核函数矩阵。确定一条从炮点到检波点的射线路径,除了要知道速度模型、炮点和检波点位置外,还需要知道射线的出射角和入射角。常规反射波走时层析反演通过若干次试验确定出射角和入射角,或者采用先求波前面再由波前面确定射线的方式,使射线追踪的计算量很大,导致走时层析反演的效率很低。
发明的内容
本发明的目的在于至少解决现有技术中存在的技术问题之一,提供一种使用地震波走时和斜率数据建立三维速度模型的方法,该方法易于计算机自动化实现,为复杂构造地区的地震成像提供了一种有效的三维速度建模方法。
本发明一种实施例提供的一种建立三维速度模型的方法,包括如下步骤:
步骤一,对原始三维地震资料进行预处理;
步骤二,拾取炮点和检波点的地震走时和斜率数据,形成观测数据空间,具体包括:利用基于曲波变换的局部相关同相轴的斜率和走时拾取算法,分别在共炮点道集和共检波点道集,追踪局部相关同相轴,并拾取相应的走时和斜率数据,对三维地震资料分别沿主测线(x方向)和交叉测线(y方向)进行操作后,获得如下的观测数据空间:
其中,TSR_obs为观测的从炮点S到检波点R的地震波反射走时;PSX_obs和PSY_obs分别为观测的在炮点处沿x和y方向的地震斜率;PRX_obs,、PRY_obs分别为观测的在检波点处沿x和y方向的地震斜率,N为拾取数据的炮点检波点对个数,即炮检对个数。
步骤三,建立离散速度模型,形成模型空间,并给定初始离散速度模型,具体包括:根据观测***资料的高程数据,确定起伏的观测面,再根据研究区的范围,确定三维速度模型在x、y和z方向的尺寸;根据地下构造的复杂程度和拟建立的速度模型分辨率要求,确定离散网格大小;对所述离散速度模型进行离散化,设置各个网格节点的速度值,所述离散速度模型中任意一点的速度值用所述网格节点速度值的三次样条函数表示,各个网格节点的速度值为模型参量,即未知量,所述模型空间为:
其中,Vj为第j个网格节点的速度值,M为网格节点数,给各网格点的速度赋值,形成初值速度模型;
步骤四,根据拾取的炮点斜率数据确定射线在炮点和在检波点的出射方向,即获得射线参数的初值,然后通过求解三维程函方程,获得当前速度模型的射线路径、走时和对应的射线参数数据,具体包括:
使用所述拾取的炮点斜率数据,用下式计算射线在所述炮点的出射方向角,以及所述炮点处的射线参数垂直分量:
其中,为炮点射线方向与过x轴平面的夹角;θS为射线方向在过x轴平面的投影与z轴的夹角;V0S为炮点处的速度;PSX_obs和PSY_obs分别为拾取的炮点处沿x和y方向的地震斜率;PSZ0为炮点处垂直慢度分量;sin-1表示反正弦函数;cos表示余弦函数,炮点处的初始射线参数PS0为
PS0=[PSX_obs,PSY_obs,PSZ0]
使用拾取的检波点斜率数据,用下式计算射线在检波点的出射方向角,以及检波点处的射线参数垂直分量:
其中,为检波点射线方向与过x轴平面的夹角;θR为射线方向在过x轴平面的投影与z轴的夹角;V0R为检波点处的速度;PRX_obs和PRY_obs分别为拾取的检波点处沿x和y方向的地震斜率;PRZ0为检波点处垂直慢度分量;sin-1表示反正弦函数;cos表示余弦函数。这时,检波点处的初始射线参数PR0为
PR0=[PRX_obs,PRY_obs,PRZ0]
分别从炮点S和检波点R处开始,使用在所述步骤S3和S4求得的炮点和检波点初始射线参数,采用Runge-Kutta法求解下列三维程函方程,
其中,x=(X,Y,Z)为空间位置向量;p(x)=[PX(x),PY(x),PZ(x)]为由三个坐标轴方向的射线参数或慢度分量组成的向量;V(x)为x处的速度;t为从炮点或检波点的地震波走时;▽表示梯度算子。
求解上述微分方程,不断求得从炮点S向地下传播的射线路径和单程走时tS,以及从检波点R向地下传播的射线路径和单程走时tR,当tS与tR之和等于该炮检对的观测双程走时TSR_obs时,停止计算,这时炮点S的射线端点位置记为xS,检波点R的射线端点位置记为xR,存储炮点和检波点的射线在相应的各个网格单元的射线长度、射线参数;
步骤五,计算分别从所述炮点和检波点确定的反射点误差及观测反射走时与计算反射走时的误差,判断当前速度模型是否满足要求,确定反演过程是否终止。
进一步,所述步骤五计算分别从所述炮点和检波点确定的反射点误差及观测反射走时与计算反射走时的误差,判断当前速度模型是否满足要求,确定反演过程是否终止,具体包括以下步骤:
S1,如果xS与xR重合,且分别沿炮点和检波点的射线的单程走时之和与观测的双程反射走时一致时,说明所用的所述离散速度模型正确,用炮检对的炮点和检波点射线末端位置的距离表示从炮点和检波点确定的反射点误差,用下列式子分别计算反射点误差及反射走时误差:
其中,(XSi,YSi,ZSi)和(XRi,YRi,ZRi)分别为从第i个炮检对的炮点和检波点的射线末端的坐标;tSi和tRi分别为从第i个炮检对的炮点和检波点到各自射线末端的走时;TSRi obs为观测的第i个炮检对的反射波走时;
S2,给定误差限,当计算的errD和errT小于给定的误差限时,输出当前速度模型,停止运行;否则,执行下列步骤;
S3,建立特定的目标函数和反演方程,求解得到对当前速度模型的修正量,更新速度模型;
S4,建立下列目标函数:
目
DSRi=[(XSi-XRi)2+(YSi-YRi)2+(ZSi-ZRi)2]1/2
TSRi=tSi+tRi-TSRi_obs
其中,N为炮检对的个数;tSi和tRi分别为从第i个炮检对的炮点和检波点到各自射线末端的走时;TSRi_obs为观测的第i个炮检对的反射波双程走时;TSRi为第i个炮检对的两个单程走时之和与观测的双程走时之差;xSi=(XSi,YSi,ZSi)和xRi=(XRi,YRi,ZRi)分别为第i个炮检对的炮点和检波点的射线末端的坐标向量;DSRi为射线末端点xSi与xRi之间的距离;v为由网格节点速度组成的模型参数向量;λ为权系数;
S5,采用高斯-牛顿法求解上述非线性目标函数的极小值解,将所述目标函数进行泰勒级数展开取至一阶项,得到线性反演方程如下:
G·Δv=Δd
其中,
Bij=λ(LSij+LRij),i=1,2,...,N,j=1,2,...,M;
Δd=[DSR1 DSR2 … DSRN TSR1 TSR2 … TSRN]T;
Δv=[V1 V2 … VM]T;
DsRi=[(XSi-XRi)2+(YSi-YRi)2+(ZSi-ZRi)2]1/2,i=1,2,...,N;
TSRi=tSi+tRi-TSRi_obs,i=1,2,...,N;
tSi和tRi分别为从第i个炮检对的炮点和检波点到各自射线末端的走时;TSRi_obs为观测的第i个炮检对的反射波走时;(XSi,YSi,ZSi)和(XRi,YRi,ZRi)分别为第i个炮检对的炮点和检波点的射线末端的坐标;Vi是第i个离散网格的速度值;VSi和VRi分别为第i个炮检对的炮点和检波点射线末端所在单元的速度;PSXi,PSYi和PSZi分别为第i个炮检对的炮点射线末端沿x,y和z方向的射线参数;PRXi,PRYi和PRZi分别为第i个炮检对的检波点射线末端沿x,y和z方向的射线参数;LSij和LRij分别为第i个炮检对的炮点射线和检波点射线在第j个单元的射线段长度;LSXij,LSYij和LSZij分别为第i个炮检对的炮点射线在第j个单元的射线段在x,y和z方向的投影长度;LRXij,LRYij和LRZij分别为第i个炮检对的检波点点射线在第j个单元的射线段在x,y和z方向的投影长度;N为炮-检对总数;M为离散模型的网格节点总数;
S6,使用LSQR算法求解上述线性反演方程,得到离散速度模型更新量Δv,用该模型更新量修改当前离散速度模型,就得到修正后的模型,计算公式如下:
v(k+1)=v(k)+Δv
其中,v(k)是当前离散速度模型,即上次(第k次)迭代的离散速度模型,Δv是本次迭代反演得到的离散速度模型修正量;v(k+1)是本次迭代得到的修正后的离散速度模型,再把修正后的离散速度模型作为新的初始模型,转到步骤四。
本发明的有益效果是:采用本发明所述的使用地震反射波走时和斜率数据建立三维速度模型的方法,用局部相关反射波同相轴的斜率数据计算从炮点和检波点向地下传播方向和初始射线参数,提高了计算射线路径和走时的效率;使用的局部相关反射波同相轴易于自动化追踪,其走时和斜率数据易于自动化拾取,对复杂地下构造和低信噪比地震资料有很强的适应性。本发明实施例易于计算机自动化实现,为复杂构造地区的地震成像提供一种有效的三维速度建模方法。
附图说明
下面结合附图和实施例对本发明进一步地说明;
图1为本发明实施例一种建立三维速度模型的方法流程图;
图2为图1所示实施例利用走时和斜率数据建立的三维速度模型切片图;
图3为常规利用走时数据的层析反演方法建立的三维速度模型切片图;
图4为使用本发明实施例建立的三维速度模型进行偏移后的共成像点道集图;
图5为使用常规走时层析反演建立的三维速度模型进行偏移后的共成像点道集图。
具体实施方式
为了使本发明目的、技术方案更加清楚明白,下面结合附图,对本发明作进一步详细说明。
在地震数据中,共炮道集和共检波点道集容易追踪局部相关同相轴,局部相关同相轴的走时随偏移距的变化率被称为地震斜率。根据射线理论,地震斜率是射线参数,也即慢度矢量的水平分量。如果已知炮点和检波点处的速度值,就可用斜率参数计算出反射射线在炮点的出射角和在检波点的入射角。而且,局部相干同相轴易于追踪,其走时和斜率也易于拾取,对低信噪比地震资料有很强的适应性。若采用光滑速度模型,通过地震斜率对射线进行约束,只需要局部相关的同相轴走时和斜率数据,而不需要事先建立各个同相轴与地层界面之间的对应关系,从而提高反射波走时层析反演的自动化程度和实用性。因此,如果把地震斜率数据引入到反射波走时层析反演中,用地震斜率确定射线出射角和入射角,就可以提高射线追踪的效率,进而形成一种实用化的建立地下低频速度模型的方法,为叠前深度偏移速度分析和波形反演提供低频或初始速度模型。
本发明实施例把地震斜率数据引入到三维地震反射走时层析反演中,用地震斜率确定射线出射角和入射角,根据当速度模型准确时,分别从炮点和检波点向地下传播的射线将相遇在反射点,且两者的走时之和必等于实测的双程反射走时的情况,提出了一种使用地震反射波走时和斜率数据建立三维速度模型的方法,为复杂构造地区地震勘探或低信噪比地震资料的偏移成像,提供一种建立地下三维速度模型的技术。
参照如图1所示,为本发明实施例一种建立三维速度模型的方法,包括如下步骤:
步骤一,对原始三维地震资料进行预处理;
步骤二,拾取炮点和检波点的地震走时和斜率数据,形成观测数据空间,具体包括:利用基于曲波变换的局部相关同相轴的斜率和走时拾取算法,分别在共炮点道集和共检波点道集,追踪局部相关同相轴,并拾取相应的走时和斜率数据,对三维地震资料分别沿主测线(x方向)和交叉测线(y方向)进行操作后,获得如下的观测数据空间:
其中,TSR_obs为观测的从炮点S到检波点R的地震波反射走时;PSX_obs和PSY_obs分别为观测的在炮点处沿x和y方向的地震斜率;PRX_obs,、PRY_obs分别为观测的在检波点处沿x和y方向的地震斜率,N为拾取数据的炮点检波点对个数,即炮检对个数。
步骤三,建立离散速度模型,形成模型空间,并给定初始离散速度模型,具体包括:根据观测***资料的高程数据,确定起伏的观测面,再根据研究区的范围,确定三维速度模型在x、y和z方向的尺寸;根据地下构造的复杂程度和拟建立的速度模型分辨率要求,确定离散网格大小;对所述离散速度模型进行离散化,设置各个网格节点的速度值,所述离散速度模型中任意一点的速度值用所述网格节点速度值的三次样条函数表示,各个网格节点的速度值为模型参量,即未知量,所述模型空间为:
其中,Vj为第j个网格节点的速度值,M为网格节点数,给各网格点的速度赋值,形成初值速度模型;
步骤四,根据拾取的炮点斜率数据确定射线在炮点和在检波点的出射方向,即获得射线参数的初值,然后通过求解三维程函方程,获得当前速度模型的射线路径、走时和对应的射线参数数据,具体包括:
使用所述拾取的炮点斜率数据,用下式计算射线在所述炮点的出射方向角,以及所述炮点处的射线参数垂直分量:
其中,为炮点射线方向与过x轴平面的夹角;θS为射线方向在过x轴平面的投影与z轴的夹角;V0S为炮点处的速度;PSX_obs和PSY_obs分别为拾取的炮点处沿x和y方向的地震斜率;PSZ0为炮点处垂直慢度分量;sin-1表示反正弦函数;cos表示余弦函数,炮点处的初始射线参数PS0为
PS0=[PSX_obs,PSY_obs,PSZ0]
使用拾取的检波点斜率数据,用下式计算射线在检波点的出射方向角,以及检波点处的射线参数垂直分量:
其中,为检波点射线方向与过x轴平面的夹角;θR为射线方向在过x轴平面的投影与z轴的夹角;V0R为检波点处的速度;PRX_obs和PRY_obs分别为拾取的检波点处沿x和y方向的地震斜率;PRZ0为检波点处垂直慢度分量;sin-1表示反正弦函数;cos表示余弦函数。这时,检波点处的初始射线参数PR0为
PR0=[PRX_obs,PRY_obs,PRZ0]
分别从炮点S和检波点R处开始,使用在所述步骤S3和S4求得的炮点和检波点初始射线参数,采用Runge-Kutta法求解下列三维程函方程,
其中,x=(X,Y,Z)为空间位置向量;p(x)=[PX(x),PY(x),PZ(x)]为由三个坐标轴方向的射线参数或慢度分量组成的向量;V(x)为x处的速度;t为从炮点或检波点的地震波走时;表示梯度算子。
求解上述微分方程,不断求得从炮点S向地下传播的射线路径和单程走时tS,以及从检波点R向地下传播的射线路径和单程走时tR,当tS与tR之和等于该炮检对的观测双程走时TSR_obs时,停止计算,这时炮点S的射线端点位置记为xS,检波点R的射线端点位置记为xR,存储炮点和检波点的射线在相应的各个网格单元的射线长度、射线参数;
步骤五,计算分别从所述炮点和检波点确定的反射点误差及观测反射走时与计算反射走时的误差,判断当前速度模型是否满足要求,确定反演过程是否终止。
进一步,所述步骤五计算分别从所述炮点和检波点确定的反射点误差及观测反射走时与计算反射走时的误差,判断当前速度模型是否满足要求,确定反演过程是否终止,具体包括以下步骤:
S1,如果xS与xR重合,且分别沿炮点和检波点的射线的单程走时之和与观测的双程反射走时一致时,说明所用的所述离散速度模型正确,用炮检对的炮点和检波点射线末端位置的距离表示从炮点和检波点确定的反射点误差,用下列式子分别计算反射点误差及反射走时误差:
其中,(XSi,YSi,ZSi)和(XRi,YRi,ZRi)分别为从第i个炮检对的炮点和检波点的射线末端的坐标;tSi和tRi分别为从第i个炮检对的炮点和检波点到各自射线末端的走时;TSRi_obs为观测的第i个炮检对的反射波走时;
S2,给定误差限,当计算的errD和errT小于给定的误差限时,输出当前速度模型,停止运行;否则,执行下列步骤;
S3,建立特定的目标函数和反演方程,求解得到对当前速度模型的修正量,更新速度模型;
S4,建立下列目标函数:
且
DSRi=[(XSi-XRi)2+(YSi-YRi)2+(ZSi-ZRi)2]1/2
TSRi=tSi+tRi-TSRi_obs
其中,N为炮检对的个数;tSi和tRi分别为从第i个炮检对的炮点和检波点到各自射线末端的走时;TSRi_obs为观测的第i个炮检对的反射波双程走时;TSRi为第i个炮检对的两个单程走时之和与观测的双程走时之差;xSi=(XSi,YSi,ZSi)和xRi=(XRi,YRi,ZRi)分别为第i个炮检对的炮点和检波点的射线末端的坐标向量;DSRi为射线末端点xSi与xRi之间的距离;v为由网格节点速度组成的模型参数向量;λ为权系数;
S5,采用高斯-牛顿法求解上述非线性目标函数的极小值解,将所述目标函数进行泰勒级数展开取至一阶项,得到线性反演方程如下:
G·Δv=Δd
其中,
Bij=λ(LSij+LRij),i=1,2,...,N,j=1,2,...,M;
Δd=[DSR1 DSR2 … DSRN TSR1 TSR2 … TSRN]T;
Δv=[V1 V2 … VM]T;
DSRi=[(XSi-XRi)2+(YSi-YRi)2+(ZSi-ZRi)2]1/2,i=1,2,...,N;
TSRi=tSi+tRi-TSRi_obs,i=1,2,...,N;
tSi和tRi分别为从第i个炮检对的炮点和检波点到各自射线末端的走时;TSRi_obs为观测的第i个炮检对的反射波走时;(XSi,YSi,ZSi)和(XRi,YRi,ZRi)分别为第i个炮检对的炮点和检波点的射线末端的坐标;Vi是第i个离散网格的速度值;VSi和VRi分别为第i个炮检对的炮点和检波点射线末端所在单元的速度;PSXi,PSYi和PSZi分别为第i个炮检对的炮点射线末端沿x,y和z方向的射线参数;PRXi,PRYi和PRZi分别为第i个炮检对的检波点射线末端沿x,y和z方向的射线参数;LSij和LRij分别为第i个炮检对的炮点射线和检波点射线在第j个单元的射线段长度;LSXij,LsYij和LSZij分别为第i个炮检对的炮点射线在第j个单元的射线段在x,y和z方向的投影长度;LRXij,LRYij和LRZij分别为第i个炮检对的检波点点射线在第j个单元的射线段在x,y和z方向的投影长度;N为炮-检对总数;M为离散模型的网格节点总数;
S6,使用LSQR算法求解上述线性反演方程,得到离散速度模型更新量Δv,用该模型更新量修改当前离散速度模型,就得到修正后的模型,计算公式如下:
v(k+1)=v(k)+Δv
其中,v(k)是当前离散速度模型,即上次(第k次)迭代的离散速度模型,Δv是本次迭代反演得到的离散速度模型修正量;v(k+1)是本次迭代得到的修正后的离散速度模型,再把修正后的离散速度模型作为新的初始模型,转到步骤四。
作为本发明实施例的一个具体示例,用我国西南某地区的三维油气地震勘探资料对本发明实施例方法进行了测试,并与仅用反射走时数据的反射层析成像的结果进行比较。
步骤100,本工区长16km,宽8km。均匀布设25条炮线和59条检波点线,炮线间距300m,检波线间距150m,炮间距100m,检波点间距50m。每炮激发时,8条检波线且每条检波线240个检波点接收。对原始地震资料进行了压制噪声、振幅补偿、反褶积、静校正等预处理。
步骤110,利用基于曲波变换的局部相关同相轴斜率和走时拾取方法,分别沿检波线和炮线方向,在共炮点域和共检波点域拾取走时和斜率数据,获得236152个炮-检对的反射走时和斜率数据,形成数据空间。
步骤120,建立大小为16km×8km×4km的速度模型,并用250m×250m×250m的网格离散化该三维速度模型,形成模型空间。模型表面用对应的高程数据构建起伏地形。模型初始速度设置成随深度线性增加,表面最高处速度为1.0km/s,模型底界处速度4.5km/s,形成初始速度模型。
步骤130,利用初始速度和拾取的炮点、检波点的斜率数据,获得分别从炮点和检波点向地下传播的地震射线方向,给出射线出射的初始条件。
步骤140,根据初始速度模型和射线初始条件,采用Runge-Kutta法求解三维程函方程,得到各个炮-检对的射线末端位置、走时以及对应的各个网格单元的射线路径和射线参数。
步骤150,用各个炮-检对的炮点和检波点射线末端位置的距离计算反射点误差,用各个炮-检对的炮点和检波点射线末端走时及实测的反射波走时计算反射走时误差。
步骤160,把计算的反射点误差和反射走时误差分别与事先设置的误差限对比,判断当前速度模型是否满足要求。若计算的反射点误差和反射走时误差都小于设置的误差限,则输出当前速度模型,终止计算;否则,进行下列步骤。
步骤170,利用观测的反射走时和斜率数据、在S5步计算的各个炮-检对的射线末端位置、走时以及对应的各个网格单元的射线路径和射线参数,建立反演方程。求解反演方程,得到对当前速度模型的修正量。
步骤180,用求得的模型修正量对当前模型进行修正,得到更新的模型。把更新模型作为新的初始模型,转到步骤S4。
图2是本发明在迭代30次后输出的地下速度模型沿中间测线的垂直切片,图3是仅用反射走时数据建立的地下速度模型沿中间测线的垂直切片。对比这两个模型发现,本发明使用走时和斜率数据建立的速度模型与常规的仅用走时数据建立的速度模型,在横向上和在垂直方向上都有较大的差别。
为了进一步说明本发明的有益效果,分别利用本发明使用走时与斜率数据建立的速度模型和常规的仅用走时数据建立的速度模型,对地震数据进行叠前深度偏移成像。图4所示是用本发明建立的速速模型进行偏移,得到沿中间测线的共成像点道集,图5所示是用常规走时层析反演建立的速度模型进行偏移,得到的沿中间测线的共成像点道集。图5中的同相轴下拉与上翘现象较为明显,同相轴基本都呈现非水平形态,而图4中的同相轴则相对较为水平和光滑。在共成像点道集中,常常通过同相轴的水平程度来判断偏移时所用速度模型的准确性。同相轴越趋于水平,说明偏移所使用的速度模型越准确。比较图4与图5可见,本发明建立的速度模型更准确,表明了本发明的良好应用效果。
以上所述仅为本发明的较佳实施例而己,并不以本发明为限制,凡在本发明的精神和原则之内所作的均等修改、等同替换和改进等,均应包含在本发明的专利涵盖范围内。
Claims (2)
1.一种利用地震数据建立三维速度模型的方法,其特征在于,包括以下步骤:
步骤一,对原始三维地震资料进行预处理;
步骤二,拾取炮点和检波点的地震走时和斜率数据,形成观测数据空间,具体包括:利用基于曲波变换的局部相关同相轴的斜率和走时拾取算法,分别在共炮点道集和共检波点道集,追踪局部相关同相轴,并拾取相应的走时和斜率数据,对三维地震资料分别沿主测线(x方向)和交叉测线(y方向)进行操作后,获得如下的观测数据空间:
其中,TSR_obs为观测的从炮点S到检波点R的地震波反射走时;PSX_obs和PSY_obs分别为观测的在炮点处沿x和y方向的地震斜率;PRX_obs,、PRY_obs分别为观测的在检波点处沿x和y方向的地震斜率,N为拾取数据的炮点检波点对个数,即炮检对个数。
步骤三,建立离散速度模型,形成模型空间,并给定初始离散速度模型,具体包括:根据观测***资料的高程数据,确定起伏的观测面,再根据研究区的范围,确定三维速度模型在x、y和z方向的尺寸;根据地下构造的复杂程度和拟建立的速度模型分辨率要求,确定离散网格大小;对所述离散速度模型进行离散化,设置各个网格节点的速度值,所述离散速度模型中任意一点的速度值用所述网格节点速度值的三次样条函数表示,各个网格节点的速度值为模型参量,即未知量,所述模型空间为:
其中,Vj为第j个网格节点的速度值,M为网格节点数,给各网格点的速度赋值,形成初值速度模型;
步骤四,根据拾取的炮点斜率数据确定射线在炮点和在检波点的出射方向,即获得射线参数的初值,然后通过求解三维程函方程,获得当前速度模型的射线路径、走时和对应的射线参数数据,具体包括:
使用所述拾取的炮点斜率数据,用下式计算射线在所述炮点的出射方向角,以及所述炮点处的射线参数垂直分量:
其中,为炮点射线方向与过x轴平面的夹角;θS为射线方向在过x轴平面的投影与z轴的夹角;V0S为炮点处的速度;PSX_obs和PSY_obs分别为拾取的炮点处沿x和y方向的地震斜率;PSZ0为炮点处垂直慢度分量;sin-1表示反正弦函数;cos表示余弦函数,炮点处的初始射线参数PS0为
PS0=[PSX_obs,PSY_obs,PSZ0]
使用拾取的检波点斜率数据,用下式计算射线在检波点的出射方向角,以及检波点处的射线参数垂直分量:
其中,为检波点射线方向与过x轴平面的夹角;θR为射线方向在过x轴平面的投影与z轴的夹角;V0R为检波点处的速度;PRX_obs和PRY_obs分别为拾取的检波点处沿x和y方向的地震斜率;PRZ0为检波点处垂直慢度分量;sin-1表示反正弦函数;cos表示余弦函数。这时,检波点处的初始射线参数PR0为
PR0=[PRX_obs,PRY_obs,PRZ0]
分别从炮点S和检波点R处开始,使用在所述步骤S3和S4求得的炮点和检波点初始射线参数,采用Runge-Kutta法求解下列三维程函方程,
其中,x=(X,Y,Z)为空间位置向量;p(x)=[PX(x),PY(x),PZ(x)]为由三个坐标轴方向的射线参数或慢度分量组成的向量;V(x)为x处的速度;t为从炮点或检波点的地震波走时;▽表示梯度算子。
求解上述微分方程,不断求得从炮点S向地下传播的射线路径和单程走时tS,以及从检波点R向地下传播的射线路径和单程走时tR,当tS与tR之和等于该炮检对的观测双程走时TSR_obs时,停止计算,这时炮点S的射线端点位置记为xS,检波点R的射线端点位置记为xR,存储炮点和检波点的射线在相应的各个网格单元的射线长度、射线参数;
步骤五,计算分别从所述炮点和检波点确定的反射点误差及观测反射走时与计算反射走时的误差,判断当前速度模型是否满足要求,确定反演过程是否终止。
2.根据权利要求1所述的一种使用地震波走时和斜率数据建立三维速度模型的方法,其特征在于,所述步骤五计算分别从所述炮点和检波点确定的反射点误差及观测反射走时与计算反射走时的误差,判断当前速度模型是否满足要求,确定反演过程是否终止,具体包括以下步骤:
S1,如果xS与xR重合,且分别沿炮点和检波点的射线的单程走时之和与观测的双程反射走时一致时,说明所用的所述离散速度模型正确,用炮检对的炮点和检波点射线末端位置的距离表示从炮点和检波点确定的反射点误差,用下列式子分别计算反射点误差及反射走时误差:
其中,(XSi,YSi,ZSi)和(XRi,YRi,ZRi)分别为从第i个炮检对的炮点和检波点的射线末端的坐标;tSi和tRi分别为从第i个炮检对的炮点和检波点到各自射线末端的走时;TSRi_obs为观测的第i个炮检对的反射波走时;
S2,给定误差限,当计算的errD和errT小于给定的误差限时,输出当前速度模型,停止运行;否则,执行下列步骤;
S3,建立特定的目标函数和反演方程,求解得到对当前速度模型的修正量,更新速度模型;
S4,建立下列目标函数:
且
DSRi=[(XSi-XRi)2+(YSi-YRi)2+(ZSi-ZRi)2]1/2
TSRi=tSi+tRi-TSRi_obs
其中,N为炮检对的个数;tSi和tRi分别为从第i个炮检对的炮点和检波点到各自射线末端的走时;TSRi_obs为观测的第i个炮检对的反射波双程走时;TSRi为第i个炮检对的两个单程走时之和与观测的双程走时之差;xSi=(XSi,YSi,ZSi)和xRi=(XRi,YRi,ZRi)分别为第i个炮检对的炮点和检波点的射线末端的坐标向量;DSRi为射线末端点xSi与xRi之间的距离;v为由网格节点速度组成的模型参数向量;λ为权系数;
S5,采用高斯-牛顿法求解上述非线性目标函数的极小值解,将所述目标函数进行泰勒级数展开取至一阶项,得到线性反演方程如下:
G·Δv=Δd
其中,
Bij=λ(LSij+LRij),i=1,2,…,N,j=1,2,…,M;
Δd=[DSR1 DSR2 … DSRN TSR1 TSR2 … TSRN]T;
Δv=[V1 V2 … VM]T;
DSRi=[(XSi-XRi)2+(YSi-YRi)2+(ZSi-ZRi)2]1/2,i=1,2,…,N;
TSRt=tSt+tRt-TSRt_obs,i=1,2,…,N;
tSi和tRi分别为从第i个炮检对的炮点和检波点到各自射线末端的走时;TSRi_obs为观测的第i个炮检对的反射波走时;(XSi,YSi,ZSi)和(XRi,YRi,ZRi)分别为第i个炮检对的炮点和检波点的射线末端的坐标;Vi是第i个离散网格的速度值;VSi和VRi分别为第i个炮检对的炮点和检波点射线末端所在单元的速度;PSXi,PSYi和PSZi分别为第i个炮检对的炮点射线末端沿x,y和z方向的射线参数;PRXi,PRYi和PRZi分别为第i个炮检对的检波点射线末端沿x,y和z方向的射线参数;LSij和LRij分别为第i个炮检对的炮点射线和检波点射线在第j个单元的射线段长度;LSXij,LSYij和LSZij分别为第i个炮检对的炮点射线在第j个单元的射线段在x,y和z方向的投影长度;LRXij,LRYij和LRZij分别为第i个炮检对的检波点点射线在第j个单元的射线段在x,y和z方向的投影长度;N为炮-检对总数;M为离散模型的网格节点总数;
S6,使用LSQR算法求解上述线性反演方程,得到离散速度模型更新量Δv,用该模型更新量修改当前离散速度模型,就得到修正后的模型,计算公式如下:
v(k+1)=v(k)+Δv
其中,v(k)是当前离散速度模型,即上次(第k次)迭代的离散速度模型,Δv是本次迭代反演得到的离散速度模型修正量;v(k+1)是本次迭代得到的修正后的离散速度模型,再把修正后的离散速度模型作为新的初始模型,转到步骤四。
Priority Applications (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN202110132180.9A CN114839675B (zh) | 2021-01-31 | 2021-01-31 | 一种建立三维速度模型的方法 |
Applications Claiming Priority (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN202110132180.9A CN114839675B (zh) | 2021-01-31 | 2021-01-31 | 一种建立三维速度模型的方法 |
Publications (2)
Publication Number | Publication Date |
---|---|
CN114839675A true CN114839675A (zh) | 2022-08-02 |
CN114839675B CN114839675B (zh) | 2023-09-05 |
Family
ID=82560816
Family Applications (1)
Application Number | Title | Priority Date | Filing Date |
---|---|---|---|
CN202110132180.9A Active CN114839675B (zh) | 2021-01-31 | 2021-01-31 | 一种建立三维速度模型的方法 |
Country Status (1)
Country | Link |
---|---|
CN (1) | CN114839675B (zh) |
Citations (9)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
US20140200814A1 (en) * | 2013-01-11 | 2014-07-17 | Cgg Services Sa | Dip tomography for estimating depth velocity models by inverting pre-stack dip information present in migrated/un-migrated pre-/post-stack seismic data |
CN105093279A (zh) * | 2014-05-16 | 2015-11-25 | 中国石油化工股份有限公司 | 针对山前带的三维地震初至波菲涅尔体层析反演方法 |
CN105319589A (zh) * | 2014-07-25 | 2016-02-10 | 中国石油化工股份有限公司 | 一种利用局部同相轴斜率的全自动立体层析反演方法 |
EP3067718A1 (en) * | 2015-03-12 | 2016-09-14 | CGG Services SA | Boundary layer tomography method and device |
CN107505651A (zh) * | 2017-06-26 | 2017-12-22 | 中国海洋大学 | 地震初至波和反射波联合斜率层析成像方法 |
WO2018021630A1 (ko) * | 2016-07-29 | 2018-02-01 | 가톨릭대학교 산학협력단 | 맞춤형 인공치아 제조 방법 및 맞춤형 인공치아 |
CN109387868A (zh) * | 2018-09-28 | 2019-02-26 | 中国海洋石油集团有限公司 | 一种基于地震波同相轴斜率信息的三维层析成像方法 |
CN109444956A (zh) * | 2019-01-09 | 2019-03-08 | 中国海洋大学 | 三维起伏观测面地震斜率层析成像方法 |
WO2019071504A1 (zh) * | 2017-10-12 | 2019-04-18 | 南方科技大学 | 一种基于两点射线追踪的地震走时层析反演方法 |
-
2021
- 2021-01-31 CN CN202110132180.9A patent/CN114839675B/zh active Active
Patent Citations (9)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
US20140200814A1 (en) * | 2013-01-11 | 2014-07-17 | Cgg Services Sa | Dip tomography for estimating depth velocity models by inverting pre-stack dip information present in migrated/un-migrated pre-/post-stack seismic data |
CN105093279A (zh) * | 2014-05-16 | 2015-11-25 | 中国石油化工股份有限公司 | 针对山前带的三维地震初至波菲涅尔体层析反演方法 |
CN105319589A (zh) * | 2014-07-25 | 2016-02-10 | 中国石油化工股份有限公司 | 一种利用局部同相轴斜率的全自动立体层析反演方法 |
EP3067718A1 (en) * | 2015-03-12 | 2016-09-14 | CGG Services SA | Boundary layer tomography method and device |
WO2018021630A1 (ko) * | 2016-07-29 | 2018-02-01 | 가톨릭대학교 산학협력단 | 맞춤형 인공치아 제조 방법 및 맞춤형 인공치아 |
CN107505651A (zh) * | 2017-06-26 | 2017-12-22 | 中国海洋大学 | 地震初至波和反射波联合斜率层析成像方法 |
WO2019071504A1 (zh) * | 2017-10-12 | 2019-04-18 | 南方科技大学 | 一种基于两点射线追踪的地震走时层析反演方法 |
CN109387868A (zh) * | 2018-09-28 | 2019-02-26 | 中国海洋石油集团有限公司 | 一种基于地震波同相轴斜率信息的三维层析成像方法 |
CN109444956A (zh) * | 2019-01-09 | 2019-03-08 | 中国海洋大学 | 三维起伏观测面地震斜率层析成像方法 |
Non-Patent Citations (2)
Title |
---|
金昌昆等: "地震立体层析成像的实现方法及效果分析", 《CT理论与应用研究》 * |
金昌昆等: "地震立体层析成像的实现方法及效果分析", 《CT理论与应用研究》, vol. 23, no. 06, 30 November 2014 (2014-11-30), pages 939 - 950 * |
Also Published As
Publication number | Publication date |
---|---|
CN114839675B (zh) | 2023-09-05 |
Similar Documents
Publication | Publication Date | Title |
---|---|---|
US7523003B2 (en) | Time lapse marine seismic surveying | |
CN108802813B (zh) | 一种多分量地震资料偏移成像方法及*** | |
US8792299B2 (en) | Method and device for processing seismic data | |
US8406081B2 (en) | Seismic imaging systems and methods employing tomographic migration-velocity analysis using common angle image gathers | |
US20140244177A1 (en) | Time lapse marine seismic surveying employing interpolated multicomponent streamer pressure data | |
US9869783B2 (en) | Structure tensor constrained tomographic velocity analysis | |
MX2011006036A (es) | Uso de inversion de forma de onda para determinar las propiedades de un medio en el subsuelo. | |
CN109444956B (zh) | 三维起伏观测面地震斜率层析成像方法 | |
EP1879052A2 (en) | Time lapse marine seismic surveying employing interpolated multicomponent streamer pressure data | |
CN109557582B (zh) | 一种二维多分量地震资料偏移成像方法及*** | |
CN110187382B (zh) | 一种回折波和反射波波动方程旅行时反演方法 | |
EA032186B1 (ru) | Сейсмическая адаптивная фокусировка | |
CN111665556B (zh) | 地层声波传播速度模型构建方法 | |
US9658354B2 (en) | Seismic imaging systems and methods employing correlation-based stacking | |
CN113466933B (zh) | 基于深度加权的地震斜率层析成像方法 | |
Zheng et al. | Double-beam stacking to infer seismic properties of fractured reservoirs | |
CN114839675B (zh) | 一种建立三维速度模型的方法 | |
CN115598704A (zh) | 一种基于最小二乘逆时偏移生成保幅角道集的方法、设备及可读存储介质 | |
CN111665546B (zh) | 用于可燃冰探测的声学参数获取方法 | |
AU2015201786A1 (en) | Methods and systems to separate wavefields using pressure wavefield data | |
CN107589446B (zh) | 利用高斯束计算波路径的层析成像速度建模方法 | |
CN115061197A (zh) | 二维海面鬼波水体成像测量方法、***、终端及测流设备 | |
CN111665550A (zh) | 地下介质密度信息反演方法 | |
CN111665549A (zh) | 地层声波衰减因子反演方法 | |
CN111665551B (zh) | 用于桥梁基底探测的声学参数获取方法 |
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 |