CN109557582A - 一种二维多分量地震资料偏移成像方法及系统 - Google Patents
一种二维多分量地震资料偏移成像方法及系统 Download PDFInfo
- Publication number
- CN109557582A CN109557582A CN201811539943.6A CN201811539943A CN109557582A CN 109557582 A CN109557582 A CN 109557582A CN 201811539943 A CN201811539943 A CN 201811539943A CN 109557582 A CN109557582 A CN 109557582A
- Authority
- CN
- China
- Prior art keywords
- wave
- iteration
- wave field
- components
- particle velocity
- 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
- 238000003384 imaging method Methods 0.000 title claims abstract description 59
- 238000013508 migration Methods 0.000 claims abstract description 112
- 230000005012 migration Effects 0.000 claims abstract description 112
- 238000000034 method Methods 0.000 claims abstract description 83
- 239000013598 vector Substances 0.000 claims abstract description 82
- 238000000354 decomposition reaction Methods 0.000 claims abstract description 11
- 239000002245 particle Substances 0.000 claims description 115
- 238000006243 chemical reaction Methods 0.000 claims description 14
- 239000011159 matrix material Substances 0.000 claims description 8
- 238000005457 optimization Methods 0.000 claims description 8
- 230000015572 biosynthetic process Effects 0.000 claims description 4
- 238000003786 synthesis reaction Methods 0.000 claims description 4
- 230000017105 transposition Effects 0.000 claims description 3
- 101100391172 Dictyostelium discoideum forA gene Proteins 0.000 claims description 2
- 238000012546 transfer Methods 0.000 claims description 2
- XLYOFNOQVPJJNP-UHFFFAOYSA-N water Substances O XLYOFNOQVPJJNP-UHFFFAOYSA-N 0.000 claims 1
- 230000008569 process Effects 0.000 abstract description 4
- 239000002131 composite material Substances 0.000 abstract description 2
- 230000008901 benefit Effects 0.000 description 6
- 238000010586 diagram Methods 0.000 description 5
- 230000008859 change Effects 0.000 description 3
- 230000000694 effects Effects 0.000 description 3
- 238000010276 construction Methods 0.000 description 2
- 230000007423 decrease Effects 0.000 description 2
- 238000011161 development Methods 0.000 description 2
- 238000005516 engineering process Methods 0.000 description 2
- 238000004880 explosion Methods 0.000 description 2
- 239000012530 fluid Substances 0.000 description 2
- 230000006872 improvement Effects 0.000 description 2
- 238000011160 research Methods 0.000 description 2
- 241000446313 Lamella Species 0.000 description 1
- 238000004458 analytical method Methods 0.000 description 1
- 238000004364 calculation method Methods 0.000 description 1
- 238000007796 conventional method Methods 0.000 description 1
- 230000008878 coupling Effects 0.000 description 1
- 238000010168 coupling process Methods 0.000 description 1
- 238000005859 coupling reaction Methods 0.000 description 1
- 238000001514 detection method Methods 0.000 description 1
- 235000013399 edible fruits Nutrition 0.000 description 1
- 238000013213 extrapolation Methods 0.000 description 1
- 238000004519 manufacturing process Methods 0.000 description 1
- 230000010287 polarization Effects 0.000 description 1
- 238000002360 preparation method Methods 0.000 description 1
- 238000012545 processing Methods 0.000 description 1
- 230000000750 progressive effect Effects 0.000 description 1
- 230000002035 prolonged effect Effects 0.000 description 1
- 230000005477 standard model Effects 0.000 description 1
- 239000004575 stone Substances 0.000 description 1
- 238000003325 tomography Methods 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
-
- 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
-
- 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)
- Geophysics And Detection Of Objects (AREA)
Abstract
本发明公开一种二维多分量地震资料偏移成像方法及系统。该方法包括：构建矢量波场延拓算子；基于矢量波场分解算子构建标量成像条件，建立梯度计算方程；基于矢量及标量波场延拓算子构建预测反射波场延拓方程，利用波场合成算子获得预测多分量波场；建立目标函数和反演流程。本发明将最小二乘反演引入多分量资料叠前深度逆时偏移成像方法，通过建立准确的多分量波场叠前深度逆时偏移及逆时反偏移算子，可获得高精度、高分辨率、高信噪比、振幅保真的多分量叠前深度逆时偏移剖面，该偏移剖面与多分量资料叠前深度逆时偏移物理意义一致，直接反映了地下介质的纵波及转换横波的反射系数信息。该方法成像精度高，偏移剖面振幅保真，物理意义明确。
Description
技术领域
本发明涉及勘探地震领域，特别是涉及一种二维多分量地震资料偏移成像方法及系统。
背景技术
地震勘探方法是利用地震波在地下介质中的传播规律研究地球内部构造、探测石油和天然气等地下资源的重要手段。地震波为弹性波，包括纵波和横波，两种波型包含不同的地下介质属性信息。因受制于技术水平和经济效益等方面因素的影响，长期以来地震勘探主要是仅利用纵波信息的纵波勘探方法，几乎忽略了包含重要应用价值的横波信息。随着油气勘探难度日益加大以及勘探技术水平的不断提高，多波多分量地震勘探技术越来越受到重视，其技术研究和实际应用也取得了长足的进步。与纵波勘探相比，多波多分量地震勘探技术可以同时利用纵波和横波信息，纵横波之间的耦合更好地保持了地震波场的运动学(走时、路径等)和动力学(波形、振幅、相位、频率、偏振特性等)特征，因此能够提供更多的地下介质信息。与纵波相比，横波(转换波)对于小断层、小幅构造、缝洞等具有更高的分辨率，对于气云区、复杂薄层构造可以获得更清晰的成像，在一些地区横波资料能能够提供更详细的构造形态、内部形变等特点，具有更好的储层刻画能力，此外，在流体识别与描述、裂缝分布估计、岩性估计、各向异性分析、亮点反射设别等方面，横波资料都更加有效，因此充分利用多波多分量地震资料可以有效提高地震勘探的精度和分辨率，减少地震勘探的多解性。
随着多波多分量地震勘探技术水平的不断提高以及高性能计算技术的迅速发展，多波多分量地震偏移成像方法逐渐在由叠前时间偏移转向叠前深度偏移、由基于射线理论的弹性波Kirchhoff偏移转向基于波动方程的弹性波逆时偏移。随着油气勘探开发程度的逐渐提高，如何对多波多分量地震资料进行准确地偏移成像已成为地震勘探亟待解决的科学难题。与基于射线理论的弹性波Kirchhoff叠前深度偏移相比，基于波动方程的弹性波叠前深度逆时偏移能够更加准确地描述地震波在地下介质中的传播规律，更完整地保持了纵波和横波的动力学和运动学信息，对地下任意复杂构造均具有成像能力，该方法能够同时获得地下介质的纵波和转换横波反射系数信息，是最具发展潜力的一类多波多分量地震偏移成像方法。尽管此类方法具有明显的理论优势，但由于实际地震数据并不完整且存在噪音，采集孔径有限，介质速度变化剧烈，偏移算子近似误差大等因素，使得弹性波叠前深度逆时偏移剖面存在明显采集脚印、分辨率低、假象严重、振幅失衡等问题，这对后续的数据处理和解释造成很大影响，使其很难直接用于实际生产，也使得多分量勘探的优势难以发挥。为此，必须建立一套新的、可获得高质量的、能够直接表征地下介质的纵波和转换横波反射系数信息的多分量地震资料叠前深度偏移方法或系统。
发明内容
本发明的目的是提供一种二维多分量地震资料偏移成像方法及系统，以克服现有技术存在的上述缺陷，可获得高质量、能够直接表征地下介质的纵波和转换横波反射系数信息。
为实现上述目的，本发明提供了如下方案：
一种二维多分量地震资料偏移成像方法，所述二维多分量地震资料偏移成像方法包括：
获取待偏移成像的观测多炮多分量观测地震记录；
获取地震工区的观测系统参数、纵波偏移速度、横波偏移速度和偏移密度模型和偏移参数；
根据所述观测系统参数、所述纵波偏移速度、所述横波偏移速度和所述偏移密度模型，获得每一炮对应的多分量正传质点速度波场；
对于第i次迭代，获取前一次迭代更新得到的预测多分量地震记录；
根据所述待偏移成像的观测多炮多分量观测地震记录和前一次迭代更新得到的预测多分量地震记录，获得当前迭代对应的多分量反传质点速度波场；
根据所述当前迭代对应的多分量反传质点速度波场和所述每一炮对应的多分量正传质点速度波场，获得当前第i次迭代对应的梯度剖面；
根据所述当前第i次迭代对应的梯度剖面，采用最优化反演算法，获得当前第i次迭代对应的下降方向剖面；
确定最优化步长；
根据所述最优化步长和所述当前第i次迭代对应的下降方向剖面，更新当前第i次迭代的偏移剖面；
判断当前第i次迭代是否满足收敛标准，得到第一判断结果；
当所述第一判断结果表示当前第i次迭代满足收敛标准，将所述当前第i次迭代的偏移剖面确定为最终偏移剖面；
当所述第一判断结果表示当前第i次迭代不满足收敛标准，迭代次数加1，返回“获取前一次迭代更新得到的预测多分量地震记录”步骤，进入第i+1次迭代。
可选的，所述根据所述观测系统参数、所述纵波偏移速度、所述横波偏移速度和所述偏移密度模型，获得每一炮对应的多分量正传质点速度波场，具体包括：
根据所述观测系统参数确定每一炮的炮点坐标；
在所述炮点坐标位置设置震源子波；
根据所述纵波偏移速度、横波偏移速度、偏移密度模型和观测系统参数，基于数值方法求解二维矢量波动方程，实现对所述炮点波场的顺时延拓，获得单炮每一个时刻的多分量正传波场，进而获得所述单炮每一个时刻的多分量正传质点速度波场，从而获得每一炮对应的多分量正传质点速度波场。
可选的，所述根据所述纵波偏移速度、横波偏移速度、偏移密度模型和观测系统参数，基于数值方法求解二维矢量波动方程，实现对所述炮点波场的顺时延拓，获得单炮每一个时刻的多分量正传波场，进而获得所述单炮每一个时刻的多分量正传质点速度波场，具体包括：
利用偏微分方程数值方法求解：
所述方程(1)中，x＝(x,z)炮点坐标矢量，t表示波传播时间，xs＝(xs,zs＝0)表示震源位置矢量；QF＝(VFx,VFz,σFxx,σFzz,σFxz)T表示多分量正传波场矢量；VF＝(VFx,VFz)T表示多分量正传质点速度波场矢量；f＝(fx,fz)T表示震源矢量；上标T表示矩阵转置运算；x和z分别表示水平方向和垂直方向；参数矩阵A1和A2具体为：
ρ为介质密度，vp为纵波偏移速度，vs为横波偏移速度。
可选的，所述根据所述待偏移成像的观测多炮多分量观测地震记录和前一次迭代更新得到的预测多分量地震记录，获得当前迭代对应的多分量反传质点速度波场，具体包括：
针对每一炮，根据所述待偏移成像的观测多炮多分量观测地震记录D和前一次迭代更新得到的预测多分量地震记录di，利用公式Δdi(xr,t；xs)＝di-D计算当前迭代对应的单炮的多分量地震记录残差Δdi(xr,t；xs)，其中Δdi(xr,t；xs)＝(Δdx,Δdz)，xs＝(xs,zs＝0)表示震源位置矢量，xr＝(xr,zr＝0)表示检波点位置矢量，t表示波传播时间；
基于数值方法对所述单炮的多分量地震记录残差进行逆时延拓，获得所述单炮的每一时刻的多分量反传波场，进而得到所述单炮的多分量反传质点速度波场；
根据所有单炮的多分量反传质点速度波场，得到当前迭代对应的多分量反传质点速度波场。
可选的，所述根据所述当前迭代对应的多分量反传质点速度波场和所述每一炮对应的多分量正传质点速度波场，获得当前第i次迭代对应的梯度剖面，具体包括：
针对每一炮，在相同时刻对单炮对应的多分量正传质点速度波场和多分量反传质点速度波场进行纵横波场分解，获得相同时刻的矢量纵波正传质点速度波场、矢量横波正传质点速度波场、矢量纵波反传质点速度波场和矢量横波反传质点速度波场；
根据所述单炮的矢量纵波正传质点速度波场、矢量纵波反传质点速度波场和矢量横波反传质点速度波场，应用梯度计算方程，获得所述单炮对应的单炮梯度剖面，所述单炮梯度剖面包括纵波剖面和转换横波剖面；
依次获得每一炮对应的单炮梯度剖面；
将所有单炮梯度剖面依据观测系统位置信息进行叠加，获得当前第i次迭代对应的梯度剖面。
可选的，所述针对每一炮，在相同时刻对单炮对应的多分量正传质点速度波场和多分量反传质点速度波场进行纵横波场分解，获得相同时刻的矢量纵波正传质点速度波场、矢量横波正传质点速度波场、矢量纵波反传质点速度波场和矢量横波反传质点速度波场，具体包括：
利用波场矢量分解方程进行纵横波场分解，具体为：
所述方程(6)中，和散度和旋度算子；PF(x,t；xs)和SF(x,t；xs)分别表示正传波场纵波势和横波势，PB(x,t；xs)和SB(x,t；xs)分别表示反传波场纵波势和横波势，为矢量纵波正传质点速度波场，为矢量横波正传质点速度波场，为矢量纵波反传质点速度波场，为矢量横波反传质点速度波场；
所述根据所述单炮的矢量纵波正传质点速度波场、矢量纵波反传质点速度波场和矢量横波反传质点速度波场，应用梯度计算方程，获得所述单炮对应的单炮梯度剖面，所述单炮梯度剖面包括纵波下降方向剖面和转换横波下降方向剖面，具体采用下列公式：
所述方程(7)中，“·”表示矢量内积运算，gPP(x)为纵波下降方向剖面，gPS(x)为转换横波下降方向剖面。
可选的，所述确定最优化步长，之后还包括：
针对每一炮，利用波场标量分解方程对所述多分量正传质点速度波场进行纵横波分解，获得该炮的纵波正传波场和转换横波正传波场；
根据该炮对应的纵波下降方向剖面和纵波正传波场，利用所述纵波偏移速度、偏移密度模型和偏移参数，基于数值方法求解二维标量波动方程，实现对该炮点预测纵波场的顺时延拓，获得该炮的每一个时刻的预测纵波场；
根据该炮对应的转换横波下降方向剖面和转换横波正传波场，利用所述横波偏移速度、偏移密度模型及偏移参数，基于数值方法求解二维标量波动方程，实现对该炮点预测转换横波场的顺时延拓，获得该炮的每一个时刻的预测转换横波场；
在相同时刻，利用波场合成方程对所述预测纵波场和所述预测转换横波场进行纵横波场合成，获得该炮的每一个时刻的预测多分量纵波场和预测多分量转换横波场，进而获得该炮的每一个时刻的预测多分量波场；
对所述预测多分量波场进行采样，获得该炮的在第i次迭代的多分量地震记录增量δdi(xr,t；xs)，δdi(xr,t；xs)＝(δdx,δdz)；其中，δdx和δdz分别表示多分量地震记录增量δd在水平方向和垂直方向的分量；
获得每一炮在当前第i次迭代的多分量地震记录增量；
利用公式di＝di-1+αiδdi更新当前第i次迭代的预测多分量地震记录di，其中，d0＝0，di-1为前一次迭代的预测多分量地震记录，αi为最优化步长，δdi为当前第i次迭代的多分量地震记录增量。
可选的，所述根据所述最优化步长和所述当前第i次迭代对应的下降方向剖面，更新当前第i次迭代的偏移剖面，具体包括：
利用公式mi＝mi-1+αidgi更新当前第i次迭代的偏移剖面mi，其中mi-1为前一次迭代的偏移剖面，αi为最优化步长，dgi为当前第i次迭代对应的下降方向剖面。
可选的，所述判断当前第i次迭代是否满足收敛标准，具体包括：
基于L-2范数，计算第i次迭代的目标函数值misfiti；
判断当前第i次迭代后的目标函数值是否满足其中，Relerr表示迭代停止的阈值标准；
所述第i次迭代后的目标函数值满足即为满足收敛标准；所述第i次迭代后的目标函数值不满足即为不满足收敛标准。
一种二维多分量地震资料偏移成像系统，所述二维多分量地震资料偏移成像系统包括：
观测多炮多分量观测地震记录获取模块，用于获取待偏移成像的观测多炮多分量观测地震记录；
参数获取模块，用于获取地震工区的观测系统参数、纵波偏移速度、横波偏移速度和偏移密度模型；
多分量正传质点速度波场获取模块，用于根据所述观测系统参数、所述纵波偏移速度、所述横波偏移速度和所述偏移密度模型，获得每一炮对应的多分量正传质点速度波场；
预测多分量地震记录获取模块，用于对于第i次迭代，获取前一次迭代更新得到的预测多分量地震记录；
多分量反传质点速度波场获取模块，用于根据所述待偏移成像的观测多炮多分量观测地震记录和前一次迭代更新得到的预测多分量地震记录，获得当前迭代对应的多分量反传质点速度波场；
梯度剖面获取模块，用于根据所述当前迭代对应的多分量反传质点速度波场和所述每一炮对应的多分量正传质点速度波场，获得当前第i次迭代对应的梯度剖面；
下降方向剖面获取模块，用于根据所述当前第i次迭代对应的梯度剖面，采用最优化反演算法，获得当前第i次迭代对应的下降方向剖面；
最优化步长确定模块，用于确定最优化步长；
偏移剖面更新模块，用于根据所述最优化步长和所述当前第i次迭代对应的下降方向剖面，更新当前第i次迭代的偏移剖面；
第一判断模块，用于判断当前第i次迭代是否满足收敛标准，得到第一判断结果；
最终偏移剖面确定模块，当所述第一判断结果表示当前第i次迭代满足收敛标准，用于将所述当前第i次迭代的偏移剖面确定为最终偏移剖面；
返回模块，用于当所述第一判断结果表示当前第i次迭代不满足收敛标准，迭代次数加1，返回“获取前一次迭代更新得到的预测多分量地震记录”步骤，进入第i+1次迭代。
根据本发明提供的具体实施例，本发明公开了以下技术效果：
1)本发明方法为一种二维多分量地震资料反演偏移方法及系统，与常规叠前深度偏移相比，本发明可以获得高精度、高分辨率、高信噪比、振幅保真的偏移剖面；2)本发明以地下介质的纵波和转换横波的反射系数信息为反演目标，通过弹性波逆时偏移和弹性波逆时反偏移，可以获得两个高质量的偏移剖面，该偏移剖面直接反映了地下介质的纵波和转换横波反射系数信息，可以直接用于后续的资料解释，极大提高了储层预测及流体识别的精度。
附图说明
为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动性的前提下，还可以根据这些附图获得其他的附图。
图1是本发明提供的二维多分量地震资料偏移成像方法的流程示意图；
图2是本发明实施例1提供的纵波速度模型；
图3是本发明实施例1提供的横波速度模型；
图4是本发明实施例1利用逆时偏移方法所得的PP偏移剖面；
图5是本发明实施例1利用逆时偏移方法所得的PS偏移剖面；
图6是本发明实施例1利用本发明方法所得的PP偏移剖面；
图7是本发明实施例1利用本发明方法所得的PS偏移剖面；
图8是本发明实施例1对应的多分量地震记录观测的水平分量；
图9是本发明实施例1对应的多分量地震记录观测的垂直分量；
图10是本发明实施例1利用逆时偏移方法所得偏移剖面预测的水平分量；
图11是本发明实施例1利用逆时偏移方法所得偏移剖面预测的垂直分量；
图12是本发明实施例1利用本发明方法所得偏移剖面预测的水平分量；
图13是本发明实施例1利用本发明方法所得偏移剖面预测的垂直分量；
图14是本发明实施例2提供的纵波速度模型；
图15是本发明实施例2提供的横波速度模型；
图16是本发明实施例2提供的密度模型；
图17是本发明实施例2利用逆时偏移方法所得的PP偏移剖面；
图18是本发明实施例2利用逆时偏移方法所得的PS偏移剖面；
图19是本发明实施例2利用本发明方法所得的PP偏移剖面；
图20是本发明实施例2利用本发明方法所得的PS偏移剖面；
图21是本发明实施例2利用逆时偏移方法所得的PP偏移剖面的局部放大图；
图22是本发明实施例2利用逆时偏移方法所得的PS偏移剖面的局部放大图；
图23是本发明实施例2利用本发明方法所得的PP偏移剖面的局部放大图；
图24是本发明实施例2利用本发明方法所得的PS偏移剖面的局部放大图；
图25是本发明实施例2利用本发明方法所得的目标函数收敛曲线示意图；
图26是本发明二维多分量地震资料偏移成像系统的结构示意图。
具体实施方式
下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。
为使本发明的上述目的、特征和优点能够更加明显易懂，下面结合附图和具体实施方式对本发明作进一步详细的说明。
图1是本发明提供的二维多分量地震资料偏移成像方法的流程示意图。如图1所示，所述二维多分量地震资料偏移成像方法具体包括：
步骤101：获取待偏移成像的观测多炮多分量观测地震记录。依据勘探目标及工区的地质地球物理条件，确定待偏移成像的观测多炮多分量观测地震记录D(xr,t；xs)＝(Dx,Dz)，其中，Dx和Dz分别表示观测的多分量观测地震记录D在水平方向和垂直方向的分量，xs＝(xs,zs＝0)表示震源位置矢量，xr＝(xr,zr＝0)表示检波点位置矢量，t表示波传播时间。
步骤102：获取地震工区的观测系统参数、纵波偏移速度、横波偏移速度和偏移密度模型和偏移参数。
步骤103：根据观测系统参数、纵波偏移速度、横波偏移速度和偏移密度模型，获得每一炮对应的多分量正传质点速度波场。具体过程如下：
针对每一炮，基于确定的观测系统参数获得该炮的炮点坐标，在该炮对应炮点位置设置震源子波，利用上述获取的纵波偏移速度、横波偏移速度、偏移密度模型及观测系统参数，基于数值方法求解二维矢量波动方程，实现对该炮点波场的顺时延拓，获得该炮的每一个时刻的多分量正传波场QF(x,t；xs)，进而获得该炮的每一个时刻的多分量正传质点速度波场VF(x,t；xs)，从而获得每一炮对应的多分量正传质点速度波场；
其中所述的基于数值方法求解二维矢量波动方程，实现对该炮点波场的顺时延拓，获得该炮的每一个时刻的多分量正传波场QF(x,t；xs)，进而获得该炮的每一个时刻的多分量正传质点速度波场VF(x,t；xs)，具体为利用偏微分方程数值方法求解：
所述方程(1)中，x＝(x,z)表示地下成像点位置坐标矢量；QF＝(VFx,VFz,σFxx,σFzz,σFxz)T表示多分量正传波场矢量；VF＝(VFx,VFz)T表示正传质点速度波场矢量；f＝(fx,fz)T表示震源矢量；上标T表示矩阵转置运算；x和z分别表示水平方向和垂直方向；参数矩阵A1和A2具体为：
ρ为介质密度，vp为纵波偏移速度，vs为横波偏移速度。
步骤104：获取前一次迭代更新得到的预测多分量地震记录。当前为第i次迭代，前一次迭代即第i-1次迭代更新所得的预测多分量地震记录为di(xr,t；xs)＝(dx,dz)。
步骤105：根据待偏移成像的观测多炮多分量观测地震记录和前一次迭代更新得到的预测多分量地震记录，获得当前迭代对应的多分量反传质点速度波场。具体如下：
当前的迭代次数为i，针对每一炮，利用该炮的观测多分量地震记录D及第i-1次迭代更新所得的预测多分量地震记录di(xr,t；xs)＝(dx,dz)，计算多分量地震记录残差Δdi(xr,t；xs)＝(Δdx,Δdz)，其计算方法为Δdi(xr,t；xs)＝di-D；基于数值方法对该炮的多分量地震记录残差Δdi进行波场逆时延拓，获得该炮的每一个时刻的多分量反传波场QB(x,t；xs)，进而可以获得多分量反传质点速度波场VB(x,t；xs)。具体为利用偏微分方程数值方法求解：
所述方程(5)中，QB＝(VBx,VBz,σBxx,σBzz,σBxz)T表示多分量反传波场矢量；VB＝(VBx,VBz)T表示反传质点速度波场矢量。
步骤106：根据当前迭代对应的多分量反传质点速度波场和每一炮对应的多分量正传质点速度波场，获得当前第i次迭代对应的梯度剖面。具体过程如下：
首先，针对每一炮，在相同时刻对单炮对应的多分量正传质点速度波场和多分量反传质点速度波场进行纵横波场分解，获得相同时刻的矢量纵波正传质点速度波场矢量横波正传质点速度波场矢量纵波反传质点速度波场和矢量横波反传质点速度波场具体的利用波场矢量分解方程进行纵横波场分解，分解方程如下：
所述方程(6)中，和散度和旋度算子；PF(x,t；xs)和SF(x,t；xs)分别表示正传波场纵波势和横波势，PB(x,t；xs)和SB(x,t；xs)分别表示反传波场纵波势和横波势。
然后，根据所述单炮的矢量纵波正传质点速度波场、矢量纵波反传质点速度波场和矢量横波反传质点速度波场，应用梯度计算方程，获得所述单炮对应的单炮梯度剖面，所述单炮梯度剖面包括纵波下降方向剖面gPP(x)和转换横波下降方向剖面gPS(x)；进而获得每一炮对应的单炮梯度剖面；所有单炮剖面依据观测系统位置信息进行叠加构成了本次迭代的梯度剖面，也即第i次迭代的梯度剖面其中，dx和dz分别表示预测的多分量地震记录d在水平方向和垂直方向的分量；Δdx和Δdz分别表示多分量地震记录残差Δd在水平方向和垂直方向的分量；表示第i次迭代的纵波梯度剖面，表示第i次迭代的转换横波梯度剖面。纵波下降方向剖面gPP(x)和转换横波下降方向剖面gPS(x)，具体计算公式为：
所述方程(7)中，“·”表示矢量内积运算。
步骤107：根据当前第i次迭代对应的梯度剖面，采用最优化反演算法，获得当前第i次迭代对应的下降方向剖面。基于第i次迭代的梯度剖面gi(x)，利用最优化反演算法，获得第i次迭代的下降方向剖面其中，表示第i次迭代的纵波下降方向剖面，表示第i次迭代的转换横波下降方向剖面。
步骤108：确定最优化步长。
步骤109：根据最优化步长和当前第i次迭代对应的下降方向剖面，更新当前第i次迭代的偏移剖面。具体的，利用公式mi＝mi-1+αidgi更新当前第i次迭代的偏移剖面mi，其中mi-1为前一次迭代的偏移剖面，αi为最优化步长，dgi为当前第i次迭代对应的下降方向剖面。偏移剖面包含PP分量剖面和PS分量剖面
步骤1010：判断当前第i次迭代是否满足收敛标准；如果是，执行步骤1011；如果否，执行步骤1012。具体的，基于L-2范数，计算第i次迭代的目标函数值misfiti，所述的收敛标准具体为：
所述方程(17)中，Relerr表示迭代停止的阈值标准，通常选取1.0e-5。
步骤1011：将当前第i次迭代的偏移剖面确定为最终偏移剖面。
步骤1012：迭代次数加1，返回步骤104，进入第i+1次迭代。
其中，本发明上述方法还包括：
读取上述计算的该炮的多分量正传质点速度波场，利用波场标量分解方程对该炮的多分量正传质点速度波场进行纵横波分解，获得该炮的纵波正传波场和转换横波正传波场基于该炮的纵波下降方向剖面dgPP和该炮的纵波正传波场利用上述确定的纵波偏移速度、偏移密度模型及偏移参数，基于数值方法求解二维标量波动方程，实现对该炮点预测纵波场的顺时延拓，获得该炮的每一个时刻的预测纵波场基于该炮的转换横波下降方向剖面dgPS和该炮的转换横波正传波场利用上述确定的横波偏移速度、偏移密度模型及偏移参数，基于数值方法求解二维标量波动方程，实现对该炮点预测转换横波场的顺时延拓，获得该炮的每一个时刻的预测转换横波场在相同的时刻，利用波场合成方程对所得的预测纵波场和预测转换横波场进行纵横波场合成，获得该炮的每一个时刻的预测多分量纵波场和预测多分量转换横波场进而获得该炮的每一个时刻的预测多分量波场VR(x,t；xs)，具体为；对该炮的预测多分量波场进行采样，获得第该炮的在第i次迭代的多分量地震记录增量δdi(xr,t；xs)＝(δdx,δdz)；进而获得每一炮在第i次迭代的多分量地震记录增量；其中，δdx和δdz分别表示多分量地震记录增量δd在水平方向和垂直方向的分量。利用最优化步长αi及多分量地震记录增量δdi，更新第i次迭代的预测多分量地震记录di＝di-1+αiδdi，其中，d0＝0。
具体的，利用波场标量分解方程对该炮的多分量正传质点速度波场进行纵横波分解，获得该炮的纵波正传波场和转换横波正传波场公式如下：
基于该炮的纵波下降方向剖面dgPP和该炮的纵波正传波场利用上述确定的纵波偏移速度模型、偏移密度模型及偏移参数，基于数值方法求解二维标量波动方程，实现对该炮点预测纵波场的顺时延拓，获得该炮的每一个时刻的预测纵波场公式如下：
所述方程(9)中，预测纵波场参数矩阵和具体为：
基于该炮的转换横波下降方向剖面dgPS和该炮的转换横波正传波场利用上述确定的横波偏移速度模型、偏移密度模型及偏移参数，基于数值方法求解二维标量波动方程，实现对该炮点预测转换横波场的顺时延拓，获得该炮的每一个时刻的预测转换横波场公式如下：
所述方程(12)中，预测转换横波场参数矩阵和具体为：
在相同的时刻，利用波场合成方程对所得的预测纵波场和预测转换横波场进行纵横波场合成，获得该炮的每一个时刻的预测多分量纵波场和预测多分量转换横波场公式如下：
所述方程(15)中，预测多分量纵波场其中和分别表示预测纵波场的水平和垂直质点速度场分量；预测多分量转换横波场其中和分别表示预测转换横波场的水平和垂直质点速度场分量；
对该炮的预测多分量波场进行采样，获得第该炮的在第i次迭代的多分量地震记录增量δdi(xr,t；xs)＝(δdx,δdz)，公式如下：
本发明涉及一种多分量地震资料偏移成像方法，本发明包括：构建矢量波场延拓算子；基于矢量波场分解算子构建标量成像条件，建立梯度计算方程；基于矢量及标量波场延拓算子构建预测反射波场延拓方程，利用波场合成算子获得预测多分量波场；建立目标函数和反演流程。本发明将最小二乘反演引入多分量资料叠前深度逆时偏移成像方法，通过建立准确的多分量波场叠前深度逆时偏移及逆时反偏移算子，可获得高精度、高分辨率、高信噪比、振幅保真的多分量叠前深度逆时偏移剖面，该偏移剖面与多分量资料叠前深度逆时偏移物理意义一致，直接反映了地下介质的纵波及转换横波的反射系数信息。该方法成像精度高，偏移剖面振幅保真，物理意义明确，可应用到各种复杂地质条件中。
为进一步说明本发明的可行性和有效性，下面举两个实例：
实施例1：采用二维层状介质模型偏移速度模型
图2-图3是二维层状介质模型偏移速度模型图，图2为纵波速度模型，图3为横波速度模型。在此模型上设置49个爆炸震源，震源子波设定为雷克子波，主频为20赫兹，起始震源点的位于(150m，100m)处，炮间隔为100m。采用中间放炮两边接收观测系统，单边最大偏移距2000m，最小偏移距为150m，道间距为10m。图4-7是二维层状介质模型的多炮叠加偏移剖面：其中，图4是利用逆时偏移方法所得的PP偏移剖面，图5是利用逆时偏移方法所得的PS偏移剖面，图6是利用本发明方法所得的PP偏移剖面，图7是利用本发明方法所得的PS偏移剖面。从图4和图5中可以看出，剖面存在比较明显的噪声，剖面的分辨率较低，振幅不均衡。从图6和图7中可以看出，本发明方法所得的偏移剖面具有很高的精度、分辨率和信噪比，而且振幅是均衡性很好，也证明了本发明方法的可行性及有效性。图8-图13是二维层状介质模型的多分量地震记录剖面：其中，图8为观测的水平分量，图9为观测的垂直分量，图10为利用逆时偏移方法所得偏移剖面预测的水平分量，图11为利用逆时偏移方法所得偏移剖面预测的垂直分量，图12为利用本发明方法所得偏移剖面预测的水平分量，图13为利用本发明方法所得偏移剖面预测的垂直分量。由图8-图13可以看出，利用逆时偏移方法所得偏移剖面进行多分量波场预测时，所得的预测多分量记录存在明显的虚假同相轴，而且，同相轴的振幅和相位均与观测的多分量记录不一致；与之相比，利用本发明方法所得偏移剖面进行多分量波场预测所得的预测多分量记录与观测的多分量地震记录一致性很好，这间接证明了本发明方法的有效性。
实施例2：采用Marmousi-2偏移模型
图14-图16是本发明提供的Marmousi-2偏移模型，图14为纵波速度模型，图15为横波速度模型，图16为密度模型。该模型是验证各种偏移方法成像效果的国际标准模型之一。在此模型上设置109个爆炸震源，震源子波设定为雷克子波，主频为15赫兹，起始震源点的位于(2550m，150m)处，炮间隔为200m。采用中间放炮两边接收观测系统，单边最大偏移距2500m，道间距为10m。图17-图20是图Marmousi-2模型的多炮叠加偏移剖面：其中，图17是利用逆时偏移方法所得的PP偏移剖面，图18是利用逆时偏移方法所得的PS偏移剖面，图19是利用本发明方法所得的PP偏移剖面，图20是利用本发明方法所得的PS偏移剖面。由图17-图20可以看出，传统方法的成像剖面振幅严重失衡，深部剖面相对振幅并不保真。而本发明方法所得成像剖面，效果更好，分辨率和精度更高，信噪比也更好，振幅均衡性更佳。图21-图24是图17-图20所示偏移剖面的局部放大图：其中，图21是利用逆时偏移方法所得的PP偏移剖面，图22是利用逆时偏移方法所得的PS偏移剖面，图23是利用本发明方法所得的PP偏移剖面，图24是利用本发明方法所得的PS偏移剖面。从图21-图24可以清晰看出，本发明方法所得的偏移剖面在薄层成像中具有更加明显的优势，此外断层清晰，断点明确。图25是利用本发明方法偏移Marmousi-2模型的目标函数收敛曲线，从中可以看出，本发明方法在复杂模型中具有很稳定的收敛性。综上所述，本发明方在复杂的地质地球物理模型中具有很好的可行性及实用性。
图26是本发明二维多分量地震资料偏移成像系统的结构示意图，如图26所示，所述二维多分量地震资料偏移成像系统包括：
观测多炮多分量观测地震记录获取模块201，用于获取待偏移成像的观测多炮多分量观测地震记录；
参数获取模块202，用于获取地震工区的观测系统参数、纵波偏移速度、横波偏移速度和偏移密度模型；
多分量正传质点速度波场获取模块203，用于根据所述观测系统参数、所述纵波偏移速度、所述横波偏移速度和所述偏移密度模型，获得每一炮对应的多分量正传质点速度波场；
预测多分量地震记录获取模块204，用于对于第i次迭代，获取前一次迭代更新得到的预测多分量地震记录；
多分量反传质点速度波场获取模块205，用于根据所述待偏移成像的观测多炮多分量观测地震记录和前一次迭代更新得到的预测多分量地震记录，获得当前迭代对应的多分量反传质点速度波场；
梯度剖面获取模块206，用于根据所述当前迭代对应的多分量反传质点速度波场和所述每一炮对应的多分量正传质点速度波场，获得当前第i次迭代对应的梯度剖面；
下降方向剖面获取模块207，用于根据所述当前第i次迭代对应的梯度剖面，采用最优化反演算法，获得当前第i次迭代对应的下降方向剖面；
最优化步长确定模块208，用于确定最优化步长；
偏移剖面更新模块209，用于根据所述最优化步长和所述当前第i次迭代对应的下降方向剖面，更新当前第i次迭代的偏移剖面；
第一判断模块2010，用于判断当前第i次迭代是否满足收敛标准，得到第一判断结果；
最终偏移剖面确定模块2011，当所述第一判断结果表示当前第i次迭代满足收敛标准，用于将所述当前第i次迭代的偏移剖面确定为最终偏移剖面；
返回模块2012，用于当所述第一判断结果表示当前第i次迭代不满足收敛标准，迭代次数加1，返回“获取前一次迭代更新得到的预测多分量地震记录”步骤，进入第i+1次迭代。
本说明书中各个实施例采用递进的方式描述，每个实施例重点说明的都是与其他实施例的不同之处，各个实施例之间相同相似部分互相参见即可。对于实施例公开的系统而言，由于其与实施例公开的方法相对应，所以描述的比较简单，相关之处参见方法部分说明即可。
本文中应用了具体个例对本发明的原理及实施方式进行了阐述，以上实施例的说明只是用于帮助理解本发明的方法及其核心思想；同时，对于本领域的一般技术人员，依据本发明的思想，在具体实施方式及应用范围上均会有改变之处。综上所述，本说明书内容不应理解为对本发明的限制。
Claims (10)
1.一种二维多分量地震资料偏移成像方法，其特征在于，所述二维多分量地震资料偏移成像方法包括：
获取待偏移成像的观测多炮多分量观测地震记录；
获取地震工区的观测系统参数、纵波偏移速度、横波偏移速度和偏移密度模型和偏移参数；
根据所述观测系统参数、所述纵波偏移速度、所述横波偏移速度和所述偏移密度模型，获得每一炮对应的多分量正传质点速度波场；
对于第i次迭代，获取前一次迭代更新得到的预测多分量地震记录；
根据所述待偏移成像的观测多炮多分量观测地震记录和前一次迭代更新得到的预测多分量地震记录，获得当前迭代对应的多分量反传质点速度波场；
根据所述当前迭代对应的多分量反传质点速度波场和所述每一炮对应的多分量正传质点速度波场，获得当前第i次迭代对应的梯度剖面；
根据所述当前第i次迭代对应的梯度剖面，采用最优化反演算法，获得当前第i次迭代对应的下降方向剖面；
确定最优化步长；
根据所述最优化步长和所述当前第i次迭代对应的下降方向剖面，更新当前第i次迭代的偏移剖面；
判断当前第i次迭代是否满足收敛标准，得到第一判断结果；
当所述第一判断结果表示当前第i次迭代满足收敛标准，将所述当前第i次迭代的偏移剖面确定为最终偏移剖面；
当所述第一判断结果表示当前第i次迭代不满足收敛标准，迭代次数加1，返回“获取前一次迭代更新得到的预测多分量地震记录”步骤，进入第i+1次迭代。
2.根据权利要求1所述的二维多分量地震资料偏移成像方法，其特征在于，所述根据所述观测系统参数、所述纵波偏移速度、所述横波偏移速度和所述偏移密度模型，获得每一炮对应的多分量正传质点速度波场，具体包括：
根据所述观测系统参数确定每一炮的炮点坐标；
在所述炮点坐标位置设置震源子波；
根据所述纵波偏移速度、横波偏移速度、偏移密度模型和观测系统参数，基于数值方法求解二维矢量波动方程，实现对所述炮点波场的顺时延拓，获得单炮每一个时刻的多分量正传波场，进而获得所述单炮每一个时刻的多分量正传质点速度波场，从而获得每一炮对应的多分量正传质点速度波场。
3.根据权利要求2所述的二维多分量地震资料偏移成像方法，其特征在于，所述根据所述纵波偏移速度、横波偏移速度、偏移密度模型和观测系统参数，基于数值方法求解二维矢量波动方程，实现对所述炮点波场的顺时延拓，获得单炮每一个时刻的多分量正传波场，进而获得所述单炮每一个时刻的多分量正传质点速度波场，具体包括：
利用偏微分方程数值方法求解：
所述方程(1)中，x＝(x,z)炮点坐标矢量，t表示波传播时间，xs＝(xs,zs＝0)表示震源位置矢量；QF＝(VFx,VFz,σFxx,σFzz,σFxz)T表示多分量正传波场矢量；VF＝(VFx,VFz)T表示多分量正传质点速度波场矢量；f＝(fx,fz)T表示震源矢量；上标T表示矩阵转置运算；x和z分别表示水平方向和垂直方向；参数矩阵A1和A2具体为：
ρ为介质密度，vp为纵波偏移速度，vs为横波偏移速度。
4.根据权利要求1所述的二维多分量地震资料偏移成像方法，其特征在于，所述根据所述待偏移成像的观测多炮多分量观测地震记录和前一次迭代更新得到的预测多分量地震记录，获得当前迭代对应的多分量反传质点速度波场，具体包括：
针对每一炮，根据所述待偏移成像的观测多炮多分量观测地震记录D和前一次迭代更新得到的预测多分量地震记录di，利用公式Δdi(xr,t；xs)＝di-D计算当前迭代对应的单炮的多分量地震记录残差Δdi(xr,t；xs)，其中Δdi(xr,t；xs)＝(Δdx,Δdz)，xs＝(xs,zs＝0)表示震源位置矢量，xr＝(xr,zr＝0)表示检波点位置矢量，t表示波传播时间；
基于数值方法对所述单炮的多分量地震记录残差进行逆时延拓，获得所述单炮的每一时刻的多分量反传波场，进而得到所述单炮的多分量反传质点速度波场；
根据所有单炮的多分量反传质点速度波场，得到当前迭代对应的多分量反传质点速度波场。
5.根据权利要求1所述的二维多分量地震资料偏移成像方法，其特征在于，所述根据所述当前迭代对应的多分量反传质点速度波场和所述每一炮对应的多分量正传质点速度波场，获得当前第i次迭代对应的梯度剖面，具体包括：
针对每一炮，在相同时刻对单炮对应的多分量正传质点速度波场和多分量反传质点速度波场进行纵横波场分解，获得相同时刻的矢量纵波正传质点速度波场、矢量横波正传质点速度波场、矢量纵波反传质点速度波场和矢量横波反传质点速度波场；
根据所述单炮的矢量纵波正传质点速度波场、矢量纵波反传质点速度波场和矢量横波反传质点速度波场，应用梯度计算方程，获得所述单炮对应的单炮梯度剖面，所述单炮梯度剖面包括纵波剖面和转换横波剖面；
依次获得每一炮对应的单炮梯度剖面；
将所有单炮梯度剖面依据观测系统位置信息进行叠加，获得当前第i次迭代对应的梯度剖面。
6.根据权利要求5所述的二维多分量地震资料偏移成像方法，其特征在于，所述针对每一炮，在相同时刻对单炮对应的多分量正传质点速度波场和多分量反传质点速度波场进行纵横波场分解，获得相同时刻的矢量纵波正传质点速度波场、矢量横波正传质点速度波场、矢量纵波反传质点速度波场和矢量横波反传质点速度波场，具体包括：
利用波场矢量分解方程进行纵横波场分解，具体为：
所述方程(6)中，▽、▽·和▽×分别表示梯度、散度和旋度算子；PF(x,t；xs)和SF(x,t；xs)分别表示正传波场纵波势和横波势，PB(x,t；xs)和SB(x,t；xs)分别表示反传波场纵波势和横波势，为矢量纵波正传质点速度波场，为矢量横波正传质点速度波场，为矢量纵波反传质点速度波场，为矢量横波反传质点速度波场；
所述根据所述单炮的矢量纵波正传质点速度波场、矢量纵波反传质点速度波场和矢量横波反传质点速度波场，应用梯度计算方程，获得所述单炮对应的单炮梯度剖面，所述单炮梯度剖面包括纵波下降方向剖面和转换横波下降方向剖面，具体采用下列公式：
所述方程(7)中，“·”表示矢量内积运算，gPP(x)为纵波下降方向剖面，gPS(x)为转换横波下降方向剖面。
7.根据权利要求6所述的二维多分量地震资料偏移成像方法，其特征在于，所述确定最优化步长，之后还包括：
针对每一炮，利用波场标量分解方程对所述多分量正传质点速度波场进行纵横波分解，获得该炮的纵波正传波场和转换横波正传波场；
根据该炮对应的纵波下降方向剖面和纵波正传波场，利用所述纵波偏移速度、偏移密度模型和偏移参数，基于数值方法求解二维标量波动方程，实现对该炮点预测纵波场的顺时延拓，获得该炮的每一个时刻的预测纵波场；
根据该炮对应的转换横波下降方向剖面和转换横波正传波场，利用所述横波偏移速度、偏移密度模型及偏移参数，基于数值方法求解二维标量波动方程，实现对该炮点预测转换横波场的顺时延拓，获得该炮的每一个时刻的预测转换横波场；
在相同时刻，利用波场合成方程对所述预测纵波场和所述预测转换横波场进行纵横波场合成，获得该炮的每一个时刻的预测多分量纵波场和预测多分量转换横波场，进而获得该炮的每一个时刻的预测多分量波场；
对所述预测多分量波场进行采样，获得该炮的在第i次迭代的多分量地震记录增量δdi(xr,t；xs)，δdi(xr,t；xs)＝(δdx,δdz)；其中，δdx和δdz分别表示多分量地震记录增量δd在水平方向和垂直方向的分量；
获得每一炮在当前第i次迭代的多分量地震记录增量；
利用公式di＝di-1+αiδdi更新当前第i次迭代的预测多分量地震记录di，其中，d0＝0，di-1为前一次迭代的预测多分量地震记录，αi为最优化步长，δdi为当前第i次迭代的多分量地震记录增量。
8.根据权利要求1所述的二维多分量地震资料偏移成像方法，其特征在于，所述根据所述最优化步长和所述当前第i次迭代对应的下降方向剖面，更新当前第i次迭代的偏移剖面，具体包括：
利用公式mi＝mi-1+αidgi更新当前第i次迭代的偏移剖面mi，其中mi-1为前一次迭代的偏移剖面，αi为最优化步长，dgi为当前第i次迭代对应的下降方向剖面。
9.根据权利要求1所述的二维多分量地震资料偏移成像方法，其特征在于，所述判断当前第i次迭代是否满足收敛标准，具体包括：
基于L-2范数，计算第i次迭代的目标函数值misfiti；
判断当前第i次迭代后的目标函数值是否满足其中，Relerr表示迭代停止的阈值标准；
所述第i次迭代后的目标函数值满足即为满足收敛标准；所述第i次迭代后的目标函数值不满足即为不满足收敛标准。
10.一种二维多分量地震资料偏移成像系统，其特征在于，所述二维多分量地震资料偏移成像系统包括：
观测多炮多分量观测地震记录获取模块，用于获取待偏移成像的观测多炮多分量观测地震记录；
参数获取模块，用于获取地震工区的观测系统参数、纵波偏移速度、横波偏移速度和偏移密度模型；
多分量正传质点速度波场获取模块，用于根据所述观测系统参数、所述纵波偏移速度、所述横波偏移速度和所述偏移密度模型，获得每一炮对应的多分量正传质点速度波场；
预测多分量地震记录获取模块，用于对于第i次迭代，获取前一次迭代更新得到的预测多分量地震记录；
多分量反传质点速度波场获取模块，用于根据所述待偏移成像的观测多炮多分量观测地震记录和前一次迭代更新得到的预测多分量地震记录，获得当前迭代对应的多分量反传质点速度波场；
梯度剖面获取模块，用于根据所述当前迭代对应的多分量反传质点速度波场和所述每一炮对应的多分量正传质点速度波场，获得当前第i次迭代对应的梯度剖面；
下降方向剖面获取模块，用于根据所述当前第i次迭代对应的梯度剖面，采用最优化反演算法，获得当前第i次迭代对应的下降方向剖面；
最优化步长确定模块，用于确定最优化步长；
偏移剖面更新模块，用于根据所述最优化步长和所述当前第i次迭代对应的下降方向剖面，更新当前第i次迭代的偏移剖面；
第一判断模块，用于判断当前第i次迭代是否满足收敛标准，得到第一判断结果；
最终偏移剖面确定模块，当所述第一判断结果表示当前第i次迭代满足收敛标准，用于将所述当前第i次迭代的偏移剖面确定为最终偏移剖面；
返回模块，用于当所述第一判断结果表示当前第i次迭代不满足收敛标准，迭代次数加1，返回“获取前一次迭代更新得到的预测多分量地震记录”步骤，进入第i+1次迭代。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201811539943.6A CN109557582B (zh) | 2018-12-17 | 2018-12-17 | 一种二维多分量地震资料偏移成像方法及系统 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201811539943.6A CN109557582B (zh) | 2018-12-17 | 2018-12-17 | 一种二维多分量地震资料偏移成像方法及系统 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN109557582A true CN109557582A (zh) | 2019-04-02 |
| CN109557582B CN109557582B (zh) | 2019-11-29 |
Family
ID=65870101
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN201811539943.6A Active CN109557582B (zh) | 2018-12-17 | 2018-12-17 | 一种二维多分量地震资料偏移成像方法及系统 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN109557582B (zh) |
Cited By (6)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN110568487A (zh) * | 2019-09-19 | 2019-12-13 | 中国科学技术大学 | 基于天然地震波形的活跃断层结构成像方法 |
| CN112305614A (zh) * | 2020-10-20 | 2021-02-02 | 中海石油(中国)有限公司 | 一种气云区空间展布范围刻画方法和系统 |
| CN112327359A (zh) * | 2020-10-14 | 2021-02-05 | 山东省科学院海洋仪器仪表研究所 | 一种基于成像能流矢量的弹性逆时偏移方法 |
| CN112462428A (zh) * | 2020-11-13 | 2021-03-09 | 中国地质科学院 | 一种多分量地震资料偏移成像方法及系统 |
| CN112904426A (zh) * | 2021-03-27 | 2021-06-04 | 中国石油大学(华东) | 一种解耦弹性波逆时偏移方法、系统及应用 |
| CN114814944A (zh) * | 2022-03-15 | 2022-07-29 | 长江大学 | 基于散度和旋度的弹性纵横波场分离和逆时偏移成像方法 |
Citations (4)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US6785612B1 (en) * | 2003-05-29 | 2004-08-31 | Pgs Americas, Inc. | Seismic velocity update for anisotropic depth migration |
| CN105974470B (zh) * | 2016-07-04 | 2017-06-16 | 中国石油大学(华东) | 一种多分量地震资料最小二乘逆时偏移成像方法及系统 |
| CN107589443A (zh) * | 2017-08-16 | 2018-01-16 | 东北石油大学 | 基于弹性波最小二乘逆时偏移成像的方法及系统 |
| CN108802813A (zh) * | 2018-06-13 | 2018-11-13 | 中国石油大学(华东) | 一种多分量地震资料偏移成像方法及系统 |
-
2018
- 2018-12-17 CN CN201811539943.6A patent/CN109557582B/zh active Active
Patent Citations (4)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US6785612B1 (en) * | 2003-05-29 | 2004-08-31 | Pgs Americas, Inc. | Seismic velocity update for anisotropic depth migration |
| CN105974470B (zh) * | 2016-07-04 | 2017-06-16 | 中国石油大学(华东) | 一种多分量地震资料最小二乘逆时偏移成像方法及系统 |
| CN107589443A (zh) * | 2017-08-16 | 2018-01-16 | 东北石油大学 | 基于弹性波最小二乘逆时偏移成像的方法及系统 |
| CN108802813A (zh) * | 2018-06-13 | 2018-11-13 | 中国石油大学(华东) | 一种多分量地震资料偏移成像方法及系统 |
Non-Patent Citations (2)
| Title |
|---|
| B. GU* ET AL.: "A Wavefield-separation-based Elastic Least Squares Reverse Time Migration", 《80TH EAGE CONFERENCE & EXHIBITION 2018》 * |
| ZONGCAI FENG ET AL.: "Elastic least-squares reverse time migration", 《ELASTIC LEAST-SQUARES REVERSE TIME MIGRATION》 * |
Cited By (11)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN110568487A (zh) * | 2019-09-19 | 2019-12-13 | 中国科学技术大学 | 基于天然地震波形的活跃断层结构成像方法 |
| CN112327359A (zh) * | 2020-10-14 | 2021-02-05 | 山东省科学院海洋仪器仪表研究所 | 一种基于成像能流矢量的弹性逆时偏移方法 |
| CN112327359B (zh) * | 2020-10-14 | 2022-06-14 | 山东省科学院海洋仪器仪表研究所 | 一种基于成像能流矢量的弹性逆时偏移方法 |
| CN112305614A (zh) * | 2020-10-20 | 2021-02-02 | 中海石油(中国)有限公司 | 一种气云区空间展布范围刻画方法和系统 |
| CN112305614B (zh) * | 2020-10-20 | 2024-03-29 | 中海石油(中国)有限公司 | 一种气云区空间展布范围刻画方法和系统 |
| CN112462428A (zh) * | 2020-11-13 | 2021-03-09 | 中国地质科学院 | 一种多分量地震资料偏移成像方法及系统 |
| CN112462428B (zh) * | 2020-11-13 | 2024-02-20 | 中国地质科学院 | 一种多分量地震资料偏移成像方法及系统 |
| CN112904426A (zh) * | 2021-03-27 | 2021-06-04 | 中国石油大学(华东) | 一种解耦弹性波逆时偏移方法、系统及应用 |
| CN112904426B (zh) * | 2021-03-27 | 2022-09-30 | 中国石油大学(华东) | 一种解耦弹性波逆时偏移方法、系统及应用 |
| CN114814944A (zh) * | 2022-03-15 | 2022-07-29 | 长江大学 | 基于散度和旋度的弹性纵横波场分离和逆时偏移成像方法 |
| CN114814944B (zh) * | 2022-03-15 | 2022-10-21 | 长江大学 | 基于散度和旋度的弹性纵横波场分离和逆时偏移成像方法 |
Also Published As
| Publication number | Publication date |
|---|---|
| CN109557582B (zh) | 2019-11-29 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN109557582B (zh) | 一种二维多分量地震资料偏移成像方法及系统 | |
| CN108802813B (zh) | 一种多分量地震资料偏移成像方法及系统 | |
| JP6982103B2 (ja) | 地下構造の検出 | |
| Kent et al. | Distribution of magma beneath the East Pacific Rise between the Clipperton transform and the 9 17′ N Deval from forward modeling of common depth point data | |
| US8406081B2 (en) | Seismic imaging systems and methods employing tomographic migration-velocity analysis using common angle image gathers | |
| CN105974470B (zh) | 一种多分量地震资料最小二乘逆时偏移成像方法及系统 | |
| CN100354654C (zh) | 用于从单色波场计算有限频率地震迁移传播时间的方法 | |
| US6188964B1 (en) | Method for using global optimization to the estimation of surface-consistent residual statics | |
| US10768323B2 (en) | Methods and systems for seismic data analysis using a tilted transversely isotropic (TTI) model | |
| Cheng et al. | Two-dimensional pre-stack reverse time imaging based on tunnel space | |
| US11609349B2 (en) | Determining properties of a subterranean formation using an acoustic wave equation with a reflectivity parameterization | |
| Hobiger et al. | Site characterization of Swiss strong‐motion stations: The benefit of advanced processing algorithms | |
| CN108594299A (zh) | 高铁智能预警方法、装置及系统 | |
| CN104991268A (zh) | 一种真振幅偏移成像方法 | |
| Kuo et al. | Comparison of three different methods in investigating shallow shear-wave velocity structures in Ilan, Taiwan | |
| CN112462427B (zh) | 多分量地震资料保幅角度域共成像点道集提取方法及系统 | |
| CN109085644A (zh) | 基于双射线走时的真地表成像方法 | |
| CN108375794A (zh) | 基于对称观测的vsp缝洞绕射成像技术方法 | |
| Fang et al. | A single station approach for subsurface anomaly detection | |
| RU2645790C1 (ru) | Способ определения границ субвертикальных протяженных объектов в геологической среде | |
| CN112462428B (zh) | 一种多分量地震资料偏移成像方法及系统 | |
| Ming‐Zhou et al. | THE THREE‐DIMENSIONAL IMAGING TECHNOLOGY AND ITS IMPROVEMENT FOR GEOLOGICAL STRUCTURE BASED ON TRT SYSTEM | |
| Pei et al. | Research and application of 5D seismic prediction technology | |
| US20230305176A1 (en) | Determining properties of a subterranean formation using an acoustic wave equation with a reflectivity parameterization | |
| Davis et al. | Multicomponent Seismic Monitoring |
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 |