CN113449440A - 一种最小二乘逆时偏移的梯度道集相关加权预处理方法 - Google Patents
一种最小二乘逆时偏移的梯度道集相关加权预处理方法 Download PDFInfo
- Publication number
- CN113449440A CN113449440A CN202110999970.7A CN202110999970A CN113449440A CN 113449440 A CN113449440 A CN 113449440A CN 202110999970 A CN202110999970 A CN 202110999970A CN 113449440 A CN113449440 A CN 113449440A
- Authority
- CN
- China
- Prior art keywords
- gradient
- gather
- reverse time
- track
- profile
- 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
- 238000013508 migration Methods 0.000 title claims abstract description 31
- 230000005012 migration Effects 0.000 title claims abstract description 31
- 238000000034 method Methods 0.000 title claims abstract description 22
- 238000007781 pre-processing Methods 0.000 title claims abstract description 18
- 230000002194 synthesizing effect Effects 0.000 claims description 3
- 238000003384 imaging method Methods 0.000 abstract description 8
- 238000004364 calculation method Methods 0.000 abstract description 7
- 238000005520 cutting process Methods 0.000 abstract description 4
- 238000001514 detection method Methods 0.000 description 3
- 238000002474 experimental method Methods 0.000 description 3
- 230000009286 beneficial effect Effects 0.000 description 2
- 238000005422 blasting Methods 0.000 description 2
- 238000004422 calculation algorithm Methods 0.000 description 2
- 238000010586 diagram Methods 0.000 description 2
- 238000002679 ablation Methods 0.000 description 1
- 230000000694 effects Effects 0.000 description 1
- 238000005516 engineering process Methods 0.000 description 1
- 230000001788 irregular Effects 0.000 description 1
- 238000004519 manufacturing process Methods 0.000 description 1
- 239000003607 modifier Substances 0.000 description 1
- 238000004088 simulation Methods 0.000 description 1
Images
Classifications
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F30/00—Computer-aided design [CAD]
- G06F30/20—Design optimisation, verification or simulation
-
- 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/301—Analysis for determining seismic cross-sections or geostructures
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
- G06F17/10—Complex mathematical operations
- G06F17/11—Complex mathematical operations for solving equations, e.g. nonlinear equations, general mathematical optimization problems
Landscapes
- Engineering & Computer Science (AREA)
- Physics & Mathematics (AREA)
- General Physics & Mathematics (AREA)
- Remote Sensing (AREA)
- Life Sciences & Earth Sciences (AREA)
- Theoretical Computer Science (AREA)
- Mathematical Physics (AREA)
- Geology (AREA)
- Pure & Applied Mathematics (AREA)
- Geophysics (AREA)
- Data Mining & Analysis (AREA)
- General Life Sciences & Earth Sciences (AREA)
- Acoustics & Sound (AREA)
- General Engineering & Computer Science (AREA)
- Environmental & Geological Engineering (AREA)
- Computational Mathematics (AREA)
- Mathematical Analysis (AREA)
- Mathematical Optimization (AREA)
- Geometry (AREA)
- Operations Research (AREA)
- Evolutionary Computation (AREA)
- Algebra (AREA)
- Computer Hardware Design (AREA)
- Databases & Information Systems (AREA)
- Software Systems (AREA)
- Geophysics And Detection Of Objects (AREA)
Abstract
本发明涉及一种最小二乘逆时偏移的梯度道集相关加权预处理方法,属于地震偏移成像领域,包括以下步骤:首先在最小二乘逆时偏移中得到每一炮的梯度剖面,将同一位置处的梯度道分别抽取出来,并按照道与炮之间的距离由小到大进行排列,得到各个位置处的梯度道集,在每个梯度道集中,对远道的噪音部分进行初步切除,并将梯度道集中的所有梯度道与参考道逐一进行相关,选取相关值较大的N道进行叠加形成一个新的梯度道,将各个新的梯度道放于对应位置处,得到最终的梯度剖面。本方法可提高梯度预处理的精度,降低梯度剖面中低频噪声的干扰,提高模型两侧地层的连续性和深层地层结构的成像精度,可减少迭代次数,解决最小二乘逆时偏移大计算量问题。
Description
技术领域
本发明属于地震速度建模领域,具体涉及一种最小二乘逆时偏移的梯度道集相关加权预处理方法。
背景技术
最小二乘逆时偏移作为一种高精度的地震速度建模算法,其每次迭代过程中,首先基于初始反射系数模型和背景速度,采用双程波方程有限差分数值模拟技术进行线性化正演得到模拟地震记录,然后将模拟地震记录与实测地震记录的差值放入L2范数中建立目标泛函,基于梯度类方法进行迭代反演使得目标泛函最小,从而得到模型修改量完成一次迭代。整个最小二乘逆时偏移通常需要几十次甚至上百次的迭代反演。
最小二乘逆时偏移可以将逆时偏移的成像问题作为最小二乘反演问题求解,通过迭代算法求解一个与观测地震数据最佳匹配的成像结果,相对于逆时偏移而言,其具有压制偏移噪音、改善深部成像、提高同相轴的连续性和均衡性以及改善振幅保真性的优势,因此,最小二乘逆时偏移在未来的实际生产中有着巨大的应用潜力。
常规最小二乘逆时偏移算法求取梯度剖面时,总的梯度剖面是由每一炮的梯度剖面直接叠加而成,实际上,当梯度道和炮点距离很远时,当前炮对此道的梯度几乎无贡献,相反,会带来不规则干扰,影响梯度剖面的质量,进而增加迭代次数和计算时间。当前并无很好地选取有效梯度炮从而提升反演收敛效率的方法。
发明内容
本发明要解决的技术问题在于提供一种最小二乘逆时偏移的梯度道集相关加权预处理方法。首先在最小二乘逆时偏移中得到每一炮的梯度剖面,在各炮的梯度剖面中将同一位置处的梯度道分别抽取出来,并按照道与炮之间的距离由小到大进行排列,以此得到各个位置处的梯度道集,然后在每个梯度道集中,对远道的噪音部分进行初步切除,并将梯度道集中的所有梯度道与参考道(即每个处理后的梯度道集所有道叠加的结果)逐一进行相关,选取相关值较大的N道(N为给定的有效叠加道数,一般取该条地震测线满覆盖次数的1/2)进行叠加,从而形成一个新的梯度道,对各个梯度道集进行上述处理,并将各个新的梯度道放于对应位置处,得到最终的梯度剖面。模型实验结果表明,本方法可提高梯度预处理的精度,减少迭代次数,节约计算时间,从而显著提升最小二乘逆时偏移的效果。
本发明采取以下技术方案:
一种最小二乘逆时偏移的梯度道集相关加权预处理方法,所述方法具体包括以下步骤:
(1)获得背景速度V(x, z)和初始反射系数模型M(x, z),其中x、z表示空间位置坐标,x=1, 2, 3, … , Nx, z=1, 2, 3, … , Nz,Nx、Nz分别代表模型横向和纵向网格点总数;给定雷克子波W i (t),i表示炮的序号,i=1, 2, 3, … , S,S表示总炮数,t表示时间,基于背景速度和初始反射系数模型,应用二阶标量声波方程进行线性化正演,得到各炮的正时波场U i (x, z, t)和合成地震记录Cal i (x, z, t);
(2)获取实际地震记录Obs i (x, z, t),将其与合成地震记录的差值作为逆时扰动进行逆时延拓得到各炮的逆推波场R i (x, z, t);
(4)对于这S个梯度剖面,将每一个剖面中的同一个x位置处的道提取出来,形成Nx个梯度道集,记为GG x (i, z),其中x=1, 2, 3, … , Nx, i=1, 2, 3, … , S,每个梯度道集中,各道按照炮点位置与当前x的距离由小到大排列;
(5)在各个梯度道集GG x (i, z)中,对噪音部分进行切除,然后利用公式(2)将所有道进行叠加,形成参考道D x (z);
(6)利用公式(3)将各个梯度道集GG x (i, z)中的各道分别与该梯度道集所对应的参考道D x (z)做相关,得到各道相对于参考道的相关值A x (i);
(7)在各个梯度道集GG x (i, z)中取相关值最大的前N道利用公式(4)进行叠加,形成新的梯度剖面GN(x, z),其中N为给定的有效叠加道数,一般取该条地震测线满覆盖次数的1/2;
最小二乘逆时偏移每次迭代中的梯度预处理均需要重复以上7个步骤。
本发明与现有技术相比的有益效果:
本发明提出的一种最小二乘逆时偏移的梯度道集相关加权预处理方法。首先在最小二乘逆时偏移中得到每一炮的梯度剖面,在各炮的梯度剖面中将同一位置处的梯度道分别抽取出来,并按照道与炮之间的距离由小到大进行排列,以此得到各个位置处的梯度道集,然后在每个梯度道集中,对远道的噪音部分进行初步切除,并将梯度道集中的所有梯度道与参考道(即每个处理后的梯度道集所有道叠加的结果)逐一进行相关,选取相关值较大的N道(N为给定的有效叠加道数,一般取该条地震测线满覆盖次数的1/2)进行叠加,从而形成一个新的梯度道,对各个梯度道集进行上述处理,并将各个新的梯度道放于对应位置处,得到最终的梯度剖面。模型实验结果表明,本方法可提高梯度预处理的精度,降低梯度剖面中低频噪声的干扰,提高模型两侧地层的连续性和深层地层结构的成像精度,将预处理后的梯度道集带入最小二乘逆时偏移的运算中,可减少迭代次数,节约计算时间,解决最小二乘逆时偏移大计算量问题,从而推动最小二乘逆时偏移在工业领域的应用进程。
附图说明
图1为最小二乘逆时偏移常规梯度预处理方法流程图;
图2为一种最小二乘逆时偏移的梯度道集相关加权预处理方法流程图;
图3 平滑速度模型;
图4为Marmousi局部反射系数模型;
图5为第400道的梯度道集图;
图6为对第400道的梯度道集切除后的图片;
图7为第400道梯度道集的参考道;
图8为常规梯度预处理后的梯度剖面;
图9为基于道集相关加权的梯度预处理后的梯度剖面;
图10 常规梯度预处理的第10次反演结果;
图11 基于梯度道集相关加权的梯度预处理的第5次反演结果。
具体实施方式
本发明以Marmousi局部模型为例(如图4所示)进行阐述,该模型横向长度为4000m,纵向长度为2500m, 模型x、z方向的网格步长均为5m。
基于该模型,一共进行200次放炮,炮间隔为20m,设置800个检波点,检波点以5m为间隔均匀的分布在模型的正上方,炮点深度和检波点深度均为0m。
下面结合图2详细阐述本发明具体实施方式:
(1)将横向长度为4000m,纵向长度为2500m的Marmousi局部模型网格化,模型x、z方向的网格步长均为5m,得到网格大小为800×500的模型,将其平滑处理得到背景速度V(x, z)并算出反射系数模型M(x, z),平滑速度模型见图3所示,其中x、z表示网格位置坐标(x=1, 2, 3, … ,800, z=1, 2, 3, … , 500);给定地震子波W i (t),进行400次放炮,i表示炮的序号(i=1, 2, 3, … , 400),t表示时间,基于背景速度和初始反射系数模型,应用二阶标量声波方程进行线性化正演,得到各炮的正时波场U i (x, z, t)和合成地震记录Cal i (x, z, t);
(2)获得实际地震记录Obs i (x, z, t),将其与合成地震记录的差值作为逆时扰动进行逆时延拓得到各炮的逆推波场R i (x, z, t);
(4)对于这400个梯度剖面,将每一个剖面中的同一个x位置处的道提取出来,形成800个梯度道集,记为GG x (i, z),其中x=1, 2, 3, … , 800, i=1, 2, 3, … , 400,每个梯度道集中,各道按照炮点位置与当前x的距离由小到大排列;图5为第400道的梯度道集图;
(5)在各个梯度道集GG x (i, z)中,对噪音部分进行切除,图6为对第400道的梯度道集切除后的图片,然后利用公式(2)将所有道进行叠加,形成参考道D x (z),图7为第400道梯度道集的参考道;
(6)利用公式(3)将各个梯度道集GG x (i, z)中的各道分别与该梯度道集所对应的参考道D x (z)做相关,得到各道相对于参考道的相关值A x (i);
(7)在各个梯度道集GG x (i, z)中取相关值最大的前400道利用公式(4)进行叠加,形成新的梯度剖面GN(x, z);
为了验证本实施例所述方法的有益效果,这里与常规梯度预处理(具体流程图见图1)后的梯度剖面进行对比。如图8为常规梯度预处理后的梯度剖面,图9为基于梯度道集相关加权的梯度预处理梯度剖面。对比两图可以发现,图9低频噪音得到消除,同相轴均衡性变好,深层成像精度更高。图10为常规梯度预处理的第10次反演结果,图11为基于梯度道集相关加权的梯度预处理的第5次反演结果,对比发现,两图精度相当,模型实验结果表明,本方法可提高梯度预处理的精度,降低梯度剖面中低频噪声的干扰,提高模型两侧地层的连续性和深层地层结构的成像精度,将预处理后的梯度道集带入最小二乘逆时偏移的运算中,可减少迭代次数,节约计算时间,解决最小二乘逆时偏移大计算量问题,从而推动最小二乘逆时偏移在工业领域的应用进程。
Claims (1)
1.一种最小二乘逆时偏移的梯度道集相关加权预处理方法,其特征在于所述方法具体包括以下步骤:
(1) 获得背景速度V(x, z)和初始反射系数模型M(x, z),其中x、z表示空间位置坐标,x=1, 2, 3, … , Nx, z=1, 2, 3, … , Nz,Nx、Nz分别代表模型横向和纵向网格点总数;给定雷克子波W i (t),i表示炮的序号,i=1, 2, 3, … , S,S表示总炮数,t表示时间,基于背景速度和初始反射系数模型,应用二阶标量声波方程进行线性化正演,得到各炮的正时波场U i (x, z, t)和合成地震记录Cal i (x, z, t);
(2)获得实际地震记录Obs i (x, z, t),将其与合成地震记录的差值作为逆时扰动进行逆时延拓得到各炮的逆推波场R i (x, z, t);
(4)对于这S个梯度剖面,将每一个剖面中的同一个x位置处的道提取出来,形成Nx个梯度道集,记为GG x (i, z),其中x=1, 2, 3, … , Nx, i=1, 2, 3, … , S,每个梯度道集中,各道按照炮点位置与当前x的距离由小到大排列;
(5)在各个梯度道集GG x (i, z)中,对噪音部分进行切除,然后利用公式(2)将所有道进行叠加,形成参考道D x (z);
(6)利用公式(3)将各个梯度道集GG x (i, z)中的各道分别与该梯度道集所对应的参考道D x (z)做相关,得到各道相对于参考道的相关值A x (i);
(7)在各个梯度道集GG x (i, z)中取相关值最大的前N道利用公式(4)进行叠加,形成新的梯度剖面GN(x, z),其中N为给定的有效叠加道数,一般取该条地震测线满覆盖次数的1/2;
Priority Applications (2)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN202110999970.7A CN113449440A (zh) | 2021-08-30 | 2021-08-30 | 一种最小二乘逆时偏移的梯度道集相关加权预处理方法 |
ZA2022/02669A ZA202202669B (en) | 2021-08-30 | 2022-03-04 | A gradient gather correlation weighted preprocessing method for least square inverse time migration |
Applications Claiming Priority (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN202110999970.7A CN113449440A (zh) | 2021-08-30 | 2021-08-30 | 一种最小二乘逆时偏移的梯度道集相关加权预处理方法 |
Publications (1)
Publication Number | Publication Date |
---|---|
CN113449440A true CN113449440A (zh) | 2021-09-28 |
Family
ID=77818913
Family Applications (1)
Application Number | Title | Priority Date | Filing Date |
---|---|---|---|
CN202110999970.7A Pending CN113449440A (zh) | 2021-08-30 | 2021-08-30 | 一种最小二乘逆时偏移的梯度道集相关加权预处理方法 |
Country Status (2)
Country | Link |
---|---|
CN (1) | CN113449440A (zh) |
ZA (1) | ZA202202669B (zh) |
Citations (3)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN107783190A (zh) * | 2017-10-18 | 2018-03-09 | 中国石油大学(北京) | 一种最小二乘逆时偏移梯度更新方法 |
CN113031063A (zh) * | 2021-04-09 | 2021-06-25 | 中国海洋大学 | 一种基于成像道集相关加权的逆时偏移成像方法 |
CN113296146A (zh) * | 2021-05-19 | 2021-08-24 | 中国海洋大学 | 一种基于梯度道集相关加权的全波形反演梯度预处理方法 |
-
2021
- 2021-08-30 CN CN202110999970.7A patent/CN113449440A/zh active Pending
-
2022
- 2022-03-04 ZA ZA2022/02669A patent/ZA202202669B/en unknown
Patent Citations (3)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN107783190A (zh) * | 2017-10-18 | 2018-03-09 | 中国石油大学(北京) | 一种最小二乘逆时偏移梯度更新方法 |
CN113031063A (zh) * | 2021-04-09 | 2021-06-25 | 中国海洋大学 | 一种基于成像道集相关加权的逆时偏移成像方法 |
CN113296146A (zh) * | 2021-05-19 | 2021-08-24 | 中国海洋大学 | 一种基于梯度道集相关加权的全波形反演梯度预处理方法 |
Also Published As
Publication number | Publication date |
---|---|
ZA202202669B (en) | 2022-06-29 |
Similar Documents
Publication | Publication Date | Title |
---|---|---|
CN111158049B (zh) | 一种基于散射积分法的地震逆时偏移成像方法 | |
CN109669212B (zh) | 地震数据处理方法、地层品质因子估算方法与装置 | |
CN113031063B (zh) | 一种基于成像道集相关加权的逆时偏移成像方法 | |
CN110058307B (zh) | 一种基于快速拟牛顿法的全波形反演方法 | |
CN108983284B (zh) | 一种适用于海上斜缆数据的f-p域鬼波压制方法 | |
CN111045077B (zh) | 一种陆地地震数据的全波形反演方法 | |
CN101021568A (zh) | 基于最大能量旅行时计算的三维积分叠前深度偏移方法 | |
CN112327362B (zh) | 速度域的海底多次波预测与追踪衰减方法 | |
CN111999764B (zh) | 基于时频域目标函数的盐下构造最小二乘逆时偏移方法 | |
CN110389377A (zh) | 基于波形互相关系数相乘的微震偏移成像定位方法 | |
CN113296146B (zh) | 一种基于梯度道集相关加权的全波形反演梯度预处理方法 | |
CN113449440A (zh) | 一种最小二乘逆时偏移的梯度道集相关加权预处理方法 | |
CN109490961B (zh) | 起伏地表无射线追踪回折波层析成像方法 | |
CN117092693A (zh) | 一种基于地震背景噪声的跨尺度面波全波形反演方法 | |
CN114624766B (zh) | 基于行波分离的弹性波最小二乘逆时偏移梯度求取方法 | |
CN112630830A (zh) | 一种基于高斯加权的反射波全波形反演方法及系统 | |
CN110873895A (zh) | 一种变网格微地震逆时干涉定位方法 | |
CN110888158B (zh) | 一种基于rtm约束的全波形反演方法 | |
CN114924313B (zh) | 基于行波分离的声波最小二乘逆时偏移梯度求取方法 | |
CN118330733B (zh) | 一种基于成像道集互相关校正的Kirchhoff积分偏移成像方法 | |
CN111239828A (zh) | 基于最优双曲线积分路径叠加的多次波压制方法 | |
CN118131322B (zh) | 一种基于Marchenko成像的微震震源定位方法 | |
CN112327361B (zh) | 基于线性同相轴迭代追踪衰减的倾斜干扰剔除方法 | |
CN111983685B (zh) | τ-p域地表非一致性静校正方法 | |
CN115993650B (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 | ||
WD01 | Invention patent application deemed withdrawn after publication |
Application publication date: 20210928 |
|
WD01 | Invention patent application deemed withdrawn after publication |