本发明公开一种地震资料偏移成像方法及系统。该方法包括：获取待进行偏移成像的地震资料数据、观测系统参数、偏移速度模型参数、偏移密度模型参数和偏移参数；获取每炮对应的震源波场；获取当前每炮对应的模拟地震记录数据，获取每炮的每一个记录时刻的伴随波场；获得当前每炮对应的单炮梯度剖面；将所有炮的第一梯度剖面叠加，获得第二梯度剖面；获得下降方向剖面、最优化步长、偏移剖面、目标函数值；判断当前目标函数值是否满足收敛标准，如果是，将当前偏移剖面确定为最终的偏移剖面；否则，再次获取每炮对应的模拟地震记录数据。本发明方法及系统，可以获得高精度、高分辨率、高信噪比、振幅保真的叠前深度偏移剖面，且成像精度和分辨率高。
1.一种地震资料偏移成像方法，其特征在于，所述方法包括：获取待进行偏移成像的地震资料数据、观测系统参数、偏移速度模型参数、偏移密度模型参数和偏移参数；所述地震资料数据为多炮的观测地震记录数据；
获取每炮对应的震源波场；所述每炮对应的震源波场包括每炮的所有记录时刻的震源波场；
获取当前每炮对应的模拟地震记录数据，所述模拟地震记录数据根据前一次每炮对应的模拟地震记录数据利用非线性反偏移算子获得；
根据当前每炮对应的观测地震记录数据和对应的模拟地震记录数据，获取每炮的每一个记录时刻的伴随波场；
根据每炮对应的震源波场和伴随波场，获得当前每炮对应的第一梯度剖面；所述第一梯度剖面为每炮对应的单炮梯度剖面；
根据所述观测系统参数和所述偏移参数，将所有炮的第一梯度剖面叠加，获得当前第二梯度剖面；
利用最优化反演算法获得当前下降方向剖面；具体包括：利用δgi＝gi-gi-1获得第i次的第二梯度剖面的增量δgi，其中gi为第i次的第二梯度剖面，gi-1为第i-1次的第二梯度剖面；
利用 获得特征参数Ai；其中αi-1为第i-1次的最优化步长，ri-1为第i-1次的下降方向剖面，gi为第i次的剖面梯度；利用 获得特征参数Bi；利用δmi＝ξiαi-1ri-1获得第i次的偏移剖面增量δmi；其中 利用获得特征参数Ci，τ为大于1的正数；利用
根据所述当前最优化步长和所述当前下降方向剖面获得当前偏移剖面；
根据 获得当前目标函数值E；其中δd表示当前地震记录数据的残差；T表示转置运算；
判断所述当前目标函数值是否满足收敛标准，得到第一判断结果；
当所述第一判断结果表示所述当前目标函数值满足收敛标准时，将所述当前偏移剖面确定为偏移成像最终的偏移剖面；
当所述第一判断结果表示所述当前目标函数值不满足收敛标准时，返回获取每炮对应的模拟地震记录数据步骤，再次获取每炮对应的模拟地震记录数据。
2.根据权利要求1所述的方法，其特征在于，所述获取每炮对应的震源波场，具体包括：根据观测系统参数获取每炮的炮点位置坐标；
根据所述偏移速度模型参数、偏移密度模型参数和偏移参数，利用交错网格有限差分法数值求解各向同性介质声波方程，实现对该炮点波场的正向延拓，获得每炮的所有记录时刻的震源波场。
3.根据权利要求1所述的方法，其特征在于，所述获取当前每炮对应的模拟地震记录数据，具体包括：利用非线性反偏移算子获得前一次的地震记录数据增量；
据 其中第i次为当前次， 为第i-1次的每炮对应的模拟地震记录数据，αi-1为第i-
1次的最优化步长，Δdi-1为第i-1次的地震记录数据增量。
4.根据权利要求3所述的方法，其特征在于，所述利用非线性反偏移算子获得前一次的地震记录数据增量，具体包括：对于第i次的第k炮，获得所述第k炮的震源波场；
根据第i次的下降方向剖面、震源波场和偏移速度模型参数构建所述第k炮的虚拟震源，建立反偏移算子方程；
根据所述反偏移算子方程，利用所述偏移速度模型参数、偏移密度模型参数对所述第k炮进行波场正向延拓，获得所述第k炮的反偏移模拟波场；
根据所述第k炮的震源波场和反偏移模拟波场，利用时间域非线性Born近似方法，确定所述第k炮的非线性模拟波场；
对所述第k炮的非线性模拟波场进行记录采样，获得第k炮的地震记录数据增量；
5.根据权利要求1所述的方法，其特征在于，所述根据当前每炮对应的观测地震记录数据和对应的模拟地震记录数据，获取每炮的每一个记录时刻的伴随波场，具体包括：根据第i次每炮对应的观测地震记录数据和对应的模拟地震记录数据，利用获得第i次第k炮对应的地震记录数据的残差δdi,k；其中 为第k炮对应的观测地震记录数据， 为第i次第k炮对应的模拟地震记录数据；
基于伴随波动方程，对所述第k炮对应的地震记录数据的残差δd进行波场逆时延拓，获得所述第k炮的每一个记录时刻的伴随波场；
6.根据权利要求5所述的方法，其特征在于，所述根据每炮对应的震源波场和伴随波场，获得当前每炮对应的第一梯度剖面，具体包括：对于第k炮，获取所述第k炮的震源波场；
根据所述第k炮的震源波场与第i次伴随波场中记录时刻的对应关系，所述第i次为当前次，利用 获得第i次所述第k炮的第一梯度剖面gi,k；ρ表示偏移密度模型参数，v表示偏移速度模型参数，i表示当前为第i次，其中Qi,k表示第i次第k炮的伴随波场，t表示时间；V＝(Vx,Vy,Vz)表示质点速度场矢量；x，y和z分别表示笛卡尔坐标的x，y和z三个方向；tmax为记录的最大时刻；
7.根据权利要求1所述的方法，其特征在于，所述利用最优化步长计算方法获取当前最优化步长，具体包括：利用 获得第i次的最优化步长αi，其中Δdi为第i次的地震记录数据增量，L为非线性反偏移算子；ri为第i次的下降方向剖面，T表示转置运算。
8.根据权利要求1所述的方法，其特征在于，所述判断所述当前目标函数值是否满足收敛标准，得到第一判断结果，具体包括：判断所述当前目标函数值是否满足 其中Ei为所述当前目标函数值，Ei-1为第i-1次的目标函数值；Theshold为设定的阈值。
9.一种地震资料偏移成像系统，其特征在于，所述系统包括：初始化模块，用于获取待进行偏移成像的地震资料数据、观测系统参数、偏移速度模型参数、偏移密度模型参数和偏移参数；所述地震资料数据为多炮的观测地震记录数据；
震源波场获取模块，用于获取每炮对应的震源波场；所述每炮对应的震源波场包括每炮的所有记录时刻的震源波场；
模拟地震记录数据获取模块，用于获取当前每炮对应的模拟地震记录数据，所述模拟地震记录数据根据前一次每炮对应的模拟地震记录数据利用非线性反偏移算子获得；
伴随波场获取模块，用于根据当前每炮对应的观测地震记录数据和对应的模拟地震记录数据，获取每炮的每一个记录时刻的伴随波场；
第一梯度剖面获取模块，用于根据每炮对应的震源波场和伴随波场，获得当前每炮对应的第一梯度剖面；所述第一梯度剖面为每炮对应的单炮梯度剖面；
第二梯度剖面获取模块，用于根据所述观测系统参数和所述偏移参数，将所有炮的第一梯度剖面叠加，获得当前第二梯度剖面；
下降方向剖面获取模块，用于利用最优化反演算法获得当前下降方向剖面；具体包括：利用δgi＝gi-gi-1获得第i次的第二梯度剖面的增量δgi，其中gi为第i次的第二梯度剖面，gi-1为第i-1次的第二梯度剖面；利用 获得特征参数Ai；其中αi-1为第i-1次的最优化步长，ri-1为第i-1次的下降方向剖面，gi为第i次的剖面梯度；利用获得特征参数Bi；利用δmi＝ξiαi-1ri-1获得第i次的偏移剖面增量δmi；其中 利用 获得特征参数Ci，τ为大于1的正数；利用 获得第i
最优化步长获取模块，用于利用最优化步长计算方法获取当前最优化步长；
偏移剖面获取模块，用于根据所述当前最优化步长和所述当前下降方向剖面获得当前偏移剖面；
目标函数值获取模块，用于根据 获得当前目标函数值E；其中δd表示当前地震记录数据的残差；T表示转置运算；
收敛判断模块，用于判断所述当前目标函数值是否满足收敛标准，得到第一判断结果；
最终的偏移剖面确定模块，用于当所述第一判断结果表示所述当前目标函数值满足收敛标准时，将所述当前偏移剖面确定为偏移成像最终的偏移剖面；
返回模块，用于当所述第一判断结果表示所述当前目标函数值不满足收敛标准时，返回获取每炮对应的模拟地震记录数据步骤，再次获取每炮对应的模拟地震记录数据。
一种地震资料偏移成像方法及系统
技术领域
[0001] 本发明涉及勘探地震领域，特别是涉及一种地震资料偏移成像方法及系统。
背景技术
[0002] 地震资料偏移成像是利用波场正向传播算子的伴随算子来实现的，该做法本质上是利用波场正向传播算子的伴随算子来代替它的逆。尽管这种近似在实际中非常有用，但它是不精确的。实际中，地震数据不是完整的，采集孔径也是有限的，而且地震数据中存在假频、噪音，数据频带有限，介质速度变化剧烈，对算子的求解也存在误差。正是由于这些原因，利用伴随算子代替逆这种做法才存在较大的误差，导致偏移剖面存在采集脚印、分辨率低、假象严重。
[0003] 现有的最小二乘偏移成像方法是针对上述问题所提出的一种地震资料叠前深度偏移成像方法，该方法利用反演的思想求解地震资料偏移问题，一定程度上常规偏移成像方法所存在的偏移剖面存在采集脚印、分辨率低、假象严重等问题。但是，现有的最小二乘偏移方法基于线性Born近似理论，属于一种线性反演偏移方法，故对偏移模型的依赖性较高，即偏移模型需与真实模型非常接近才可获得较好的偏移效果。实际中，偏移模型往往包含较大误差，因此对于实际资料，现有的最小二乘偏移方法并不能获得理想的偏移剖面。为此，必须建立一套新的、以“高精度、高分辨率、高信噪比、保幅、低速度依赖性”为特点的地震资料偏移成像方法。
发明内容
[0004] 本发明的目的是提供一种地震资料偏移成像方法及系统，以提高地震资料偏移成像的精度。
[0005] 为实现上述目的，本发明提供了如下方案：
[0006] 一种地震资料偏移成像方法，所述方法包括：
[0007] 获取待进行偏移成像的地震资料数据、观测系统参数、偏移速度模型参数、偏移密度模型参数和偏移参数；所述地震资料数据为多炮的观测地震记录数据；
[0008] 获取每炮对应的震源波场；所述每炮对应的震源波场包括每炮的所有记录时刻的震源波场；
[0009] 获取当前每炮对应的模拟地震记录数据，所述模拟地震记录数据根据前一次每炮对应的模拟地震记录数据利用非线性反偏移算子获得；
[0010] 根据当前每炮对应的观测地震记录数据和对应的模拟地震记录数据，获取每炮的每一个记录时刻的伴随波场；
[0011] 根据每炮对应的震源波场和伴随波场，获得当前每炮对应的第一梯度剖面；所述第一梯度剖面为每炮对应的单炮梯度剖面；
[0012] 根据所述观测系统参数和所述偏移参数，将所有炮的第一梯度剖面叠加，获得当前第二梯度剖面；
[0013] 利用最优化反演算法获得当前下降方向剖面；
[0014] 利用最优化步长计算方法获取当前最优化步长；
[0015] 根据所述当前最优化步长和所述当前下降方向剖面获得当前偏移剖面；
[0016] 根据 获得当前目标函数值E；其中δd表示当前地震记录数据的残差；T表示转置运算；
[0017] 判断所述当前目标函数值是否满足收敛标准，得到第一判断结果；
[0018] 当所述第一判断结果表示所述当前目标函数值满足收敛标准时，将所述当前偏移剖面确定为偏移成像最终的偏移剖面；
[0019] 当所述第一判断结果表示所述当前目标函数值不满足收敛标准时，返回获取每炮对应的模拟地震记录数据步骤，再次获取每炮对应的模拟地震记录数据。
[0020] 可选的，所述获取每炮对应的震源波场，具体包括：
[0021] 根据观测系统参数获取每炮的炮点位置坐标；
[0022] 在每个炮点位置坐标处分别设置震源子波；
[0023] 根据所述偏移速度模型参数、偏移密度模型参数和偏移参数，利用交错网格有限差分法数值求解各向同性介质声波方程，实现对该炮点波场的正向延拓，获得每炮的所有记录时刻的震源波场。
[0024] 可选的，所述获取当前每炮对应的模拟地震记录数据，具体包括：
[0025] 利用非线性反偏移算子获得前一次的地震记录数据增量；
[0026] 利用公式 获得第i次的每炮对应的模拟地震记录数据 其中第i次为当前次， 为第i-1次的每炮对应的模拟地震记录数据，αi-1为第i-1次的最优化步长，Δdi-1为第i-1次的地震记录数据增量。
[0027] 可选的，所述利用非线性反偏移算子获得前一次的地震记录数据增量，具体包括：
[0028] 对于第k炮，获得所述第k炮的震源波场；
[0029] 根据第i次的下降方向剖面、震源波场和偏移速度模型参数构建所述第k炮的虚拟震源，建立反偏移算子方程；
[0030] 根据所述反偏移算子方程，利用所述偏移速度模型参数、偏移密度模型参数对所述第k炮进行波场正向延拓，获得所述第k炮的反偏移模拟波场；
[0031] 根据所述第k炮的震源波场和反偏移模拟波场，利用时间域非线性Born近似方法，确定所述第k炮的非线性模拟波场；
[0032] 对所述第k炮的非线性模拟波场进行记录采样，获得第k炮的地震记录数据增量；
[0033] 依次获得所有炮的第i次的地震记录数据增量。
[0034] 可选的，所述根据当前每炮对应的观测地震记录数据和对应的模拟地震记录数据，获取每炮的每一个记录时刻的伴随波场，具体包括：
[0035] 根据第i次每炮对应的观测地震记录数据和对应的模拟地震记录数据，利用获得第i次第k炮对应的地震记录数据的残差δdi,k；其中d 为第k炮对应的观测地震记录数据， 为第i次第k炮对应的模拟地震记录数据；
[0036] 基于伴随波动方程，对所述第k炮对应的地震记录数据的残差δd进行波场逆时延拓，获得所述第k炮的每一个记录时刻的伴随波场；
[0037] 依次获得每炮的每一个记录时刻的伴随波场。
[0038] 可选的，所述根据每炮对应的震源波场和伴随波场，获得当前每炮对应的第一梯度剖面，具体包括：
[0039] 对于第k炮，获取所述第k炮的震源波场；
[0040] 根据所述第k炮的震源波场与第i次伴随波场中记录时刻的对应关系，所述第i次为当前次，利用 获得第i次所述第k炮的第一梯度剖面gi,k；ρ表示偏移密度模型参数，v表示偏移速度模型参数，i表示当前为第i次，其中Qi,k表示第i次第k炮的伴随波场，t表示时间；V＝(Vx,Vy,Vz)表示质点速度场矢量；x，y和z分别表示笛卡尔坐标的x，y和z三个方向；tmax为记录的最大时刻；
[0041] 依次获得第i次每炮对应的第一梯度剖面。
[0042] 可选的，所述利用最优化反演算法获得当前下降方向剖面，具体包括：
[0043] 利用δgi＝gi-gi-1获得第i次的第二梯度剖面的增量δgi，其中gi为第i次的第二梯度剖面，gi-1为第i-1次的第二梯度剖面；
[0044] 利用 获得特征参数Ai；其中αi-1为第i-1次的最优化步长，ri-1为第i-1次的下降方向剖面，gi为第i次的剖面梯度；
[0045] 利用 获得特征参数Bi；
[0046] 利 用δmi ＝ξiαi-1ri- 1获得第i次的偏 移剖面增量δmi ；其中[0047] 利用 获得特征参数Ci，τ为大于1的正数；
[0048] 利用 获得第i次的下降方向剖面ri。
[0049] 可选的，所述利用最优化步长计算方法获取当前最优化步长，具体包括：
[0050] 利用 获得第i次的最优化步长αi，其中Δdi为第i次的地震记录数据增量，L为非线性反偏移算子；ri为第i次的下降方向剖面，T表示转置运算。
[0051] 可选的，所述判断所述当前目标函数值是否满足收敛标准，得到第一判断结果，具体包括：
[0052] 判断所述当前目标函数值是否满足 其中Ei为所述当前目标函数值，Ei-1为第i-1次的目标函数值；Theshold为设定的阈值。
[0053] 本发明还提供一种地震资料偏移成像系统，所述系统应用于上述的方法，所述系统包括：
[0054] 初始化模块，用于获取待进行偏移成像的地震资料数据、观测系统参数、偏移速度模型参数、偏移密度模型参数和偏移参数；所述地震资料数据为多炮的观测地震记录数据；
[0055] 震源波场获取模块，用于获取每炮对应的震源波场；所述每炮对应的震源波场包括每炮的所有记录时刻的震源波场；
[0056] 模拟地震记录数据获取模块，用于获取当前每炮对应的模拟地震记录数据，所述模拟地震记录数据根据前一次每炮对应的模拟地震记录数据利用非线性反偏移算子获得；
[0057] 伴随波场获取模块，用于根据当前每炮对应的观测地震记录数据和对应的模拟地震记录数据，获取每炮的每一个记录时刻的伴随波场；
[0058] 第一梯度剖面获取模块，用于根据每炮对应的震源波场和伴随波场，获得当前每炮对应的第一梯度剖面；所述第一梯度剖面为每炮对应的单炮梯度剖面；
[0059] 第二梯度剖面获取模块，用于根据所述观测系统参数和所述偏移参数，将所有炮的第一梯度剖面叠加，获得当前第二梯度剖面；
[0060] 下降方向剖面获取模块，用于利用最优化反演算法获得当前下降方向剖面；
[0061] 最优化步长获取模块，用于利用最优化步长计算方法获取当前最优化步长；
[0062] 偏移剖面获取模块，用于根据所述当前最优化步长和所述当前下降方向剖面获得当前偏移剖面；
[0063] 目标函数值获取模块，用于根据 获得当前目标函数值E；其中δd表示当前地震记录数据的残差；T表示转置运算；
[0064] 收敛判断模块，用于判断所述当前目标函数值是否满足收敛标准，得到第一判断结果；
[0065] 最终的偏移剖面确定模块，用于当所述第一判断结果表示所述当前目标函数值满足收敛标准时，将所述当前偏移剖面确定为偏移成像最终的偏移剖面；
[0066] 返回模块，用于当所述第一判断结果表示所述当前目标函数值不满足收敛标准时，返回获取每炮对应的模拟地震记录数据步骤，再次获取每炮对应的模拟地震记录数据。
[0067] 根据本发明提供的具体实施例，本发明公开了以下技术效果：
[0068] 1)本发明将反演的思想引入波逆时偏移中，与常规逆时偏移相比，本发明可以获得高精度、高分辨率、高信噪比、振幅保真的偏移剖面；2)本发明所用时间域非线性Born近似理论，较好地弥补了常规Born近似理论因忽略多次散射波所造成的成像误差，降低偏移方法对偏移模型的依赖程度；3)本发明所用的最优化反演算法为自适应共轭梯度算法，能够自主确定偏移算法在每次迭代中的下降方向剖面，与常规梯度类反演算法相比具有更高的收敛速率和稳定性；4)本发明所得成像剖面物理意义非常明确，便于后期地质解译；5)本发明可以广泛用于油气勘探领域中，特别是对于复杂构造深部的成像效果更加明显。
附图说明
[0069] 为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动性的前提下，还可以根据这些附图获得其他的附图。
[0070] 图1为本发明地震资料偏移成像方法的流程示意图；
[0071] 图2为本发明地震资料偏移成像系统的结构示意图；
[0072] 图3是本发明提供的二维洼陷偏移速度模型图；
[0073] 图4是图3所示二维洼陷模型的多炮叠加偏移剖面：其中，图4(a)是利用逆时偏移方法所得的偏移剖面，图4(b)是利用本发明方法所得的偏移剖面；
[0074] 图5是本发明提供的Marmousi-2偏移速度模型；
[0075] 图6是图5所示Marmousi-2模型的多炮叠加偏移剖面：其中，图6(a)是利用逆时偏移方法所得的偏移剖面；图6(b)是利用本发明方法所得的偏移剖面。
具体实施方式
[0076] 下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。
[0077] 为使本发明的上述目的、特征和优点能够更加明显易懂，下面结合附图和具体实施方式对本发明作进一步详细的说明。
[0078] 图1为本发明地震资料偏移成像方法的流程示意图。如图1所示，所述方法包括：
[0079] 步骤101：获取待进行偏移成像的地震资料数据、观测系统参数、偏移速度模型参数、偏移密度模型参数和偏移参数；所述地震资料数据为多炮的观测地震记录数据。此步骤为初始化过程，首先读取地震工区观测系统参数、速度模型及观测的多炮地震资料，由于读取的是全部的资料，需要根据偏移成像的要求，从中获取待进行偏移成像的地震资料数据，针对此步骤获取的数据，进行后续的处理。
[0080] 步骤102：获取每炮对应的震源波场。所述每炮对应的震源波场包括每炮的所有记录时刻的震源波场。针对每一炮，基于读取的观测系统参数获得该炮的炮点坐标，在该炮对应炮点位置设置震源子波，根据偏移速度模型参数、偏移密度模型参数和偏移参数，利用交错网格有限差分法数值求解各向同性介质声波方程，实现对该炮点波场的正向延拓，获得每炮的所有记录时刻的震源波场。其中所述的各向同性介质声波方程为：
[0081]
[0082] 所述方程(1)中，P表示波场，V＝(Vx,Vy,Vz)表示质点速度场矢量；t表示时间；x，y和z分别表示笛卡尔坐标的x，y和z三个方向。
[0083] 步骤103：获取每炮的每一个记录时刻的伴随波场。从此步骤开始进行迭代，设置当前的迭代次数为i，针对每一炮，利用该炮的观测地震记录dobs及第i-1次迭代后所得的模拟地震记录 计算地震记录残差δdi， 基于伴随波动方程，对该炮的地震记录残差δdi进行波场逆时延拓，获得该炮的每一个记录时刻的伴随波场，依次获得每一炮的每一个记录时刻的伴随波场。
[0084] 其中具体过程：
[0085] (1)利用公式 获得第i-1次迭代后所得的第i次的每炮对应的模拟地震记录数据 其中第i次为当前次， 为第i-1次的每炮对应的模拟地震记录数据，αi-1为第i-1次的最优化步长，Δdi-1为第i-1次的地震记录数据增量。
[0086] (2)数据增量的获取方式(以Δdi为例)为：
[0087] 对于第k炮，获得所述第k炮的震源波场；
[0088] 根据第i次的下降方向剖面、震源波场和偏移速度模型参数构建所述第k炮的虚拟震源，建立反偏移算子方程；
[0089] 根据所述反偏移算子方程，利用所述偏移速度模型参数、偏移密度模型参数对所述第k炮进行波场正向延拓，获得所述第k炮的反偏移模拟波场；
[0090] 根据所述第k炮的震源波场和反偏移模拟波场，利用时间域非线性Born近似方法，确定所述第k炮的非线性模拟波场；
[0091] 对所述第k炮的非线性模拟波场进行记录采样，获得第k炮的地震记录数据增量；
[0092] 依次获得所有炮的第i次的地震记录数据增量。
[0093] 其中，根据第i次的下降方向剖面、震源波场和偏移速度模型参数构建所述第k炮的虚拟震源，建立反偏移算子方程，具体为：
[0094]
[0095] 所述方程(2)中，δP表示反偏移波场，δV＝(δVx,δVy,δVz)表示反偏移质点速度场矢量；
[0096] 所述虚拟震源，具体为：
[0097]
[0098] 其中根据所述第k炮的震源波场和反偏移模拟波场，利用时间域非线性Born近似方法，确定所述第k炮的非线性模拟波场；对所述第k炮的非线性模拟波场进行记录采样，获得第k炮的地震记录数据增量，具体为：
[0099]
[0100] 式(4)中，f(·)表示傅里叶正变换，f-1(·)表示傅里叶逆变换，δ(·)为δ函数，x为地下介质空间位置矢量，xr为检波点空间位置矢量。
[0101] (3)伴随波动方程为：
[0102]
[0103] 所述方程(5)中，Q表示伴随波场，U＝(Ux,Uy,Uz)表示伴随质点速度场矢量；t表示时间；x，y和z分别表示笛卡尔坐标的x，y和z三个方向。
[0104] 步骤104：根据每炮对应的震源波场和伴随波场，获得当前每炮对应的第一梯度剖面。所述第一梯度剖面为每炮对应的单炮梯度剖面。
[0105] 对于第k炮，获取所述第k炮的震源波场；
[0106] 根据所述第k炮的震源波场与第i次伴随波场中记录时刻的对应关系，所述第i次为当前次，也就是按照相同时刻数据对应的关系 ，利用梯度计算公式获得第i次所述第k炮的第一梯度剖面gi,k；
ρ表示偏移密度模型参数，v表示偏移速度模型参数，i表示当前为第i次，其中Qi,k表示第i次第k炮的伴随波场，t表示时间，tmax为记录的最大时刻；V＝(Vx,Vy,Vz)表示质点速度场矢量；
x，y和z分别表示笛卡尔坐标的x，y和z三个方向；
[0107] 依次获得第i次每炮对应的第一梯度剖面，也就是当前每炮对应的弹炮梯度剖面。
[0108] 步骤105：根据所述观测系统参数和所述偏移参数，将所有炮的第一梯度剖面叠加，获得当前第二梯度剖面。第二梯度剖面为叠加后的梯度剖面，也即为本次迭代获得的梯度剖面。
[0109] 步骤106：利用最优化反演算法获得当前下降方向剖面。此步骤应用的最优化反演算法为一种自适应共轭梯度法，具体过程为：
[0110] 利用δgi＝gi-gi-1获得第i次的第二梯度剖面的增量δgi，其中gi为第i次的第二梯度剖面，gi-1为第i-1次的第二梯度剖面；
[0111] 利用 获得特征参数Ai；其中αi-1为第i-1次的最优化步长，ri-1为第i-1次的下降方向剖面，gi为第i次的剖面梯度；
[0112] 利用 获得特征参数Bi；
[0113] 利 用δmi ＝ξiαi-1ri- 1获得第i次的偏 移剖面增量δmi ；其中[0114] 利用 获得特征参数Ci，τ为大于1的正数；
[0115] 利用 获得第i次的下降方向剖面ri。
[0116] 步骤107：利用最优化步长计算方法获取当前最优化步长。具体采用获得第i次的最优化步长αi，其中Δdi为第i次的地震记录数据增量，L为非线性反偏移算子；ri为第i次的下降方向剖面，T表示转置运算。
[0117] 步骤108：根据所述当前最优化步长和所述当前下降方向剖面获得当前偏移剖面。利用mi＝mi-1+αiri获得当前的偏移剖面mi。
[0118] 步骤109：根据 获得当前目标函数值E；其中δd表示当前地震记录数据的残差；T表示转置运算；确定目标函数值是为了后续判断是否满足收敛标准。也可以采用其他的目标函数值，只要能实现判断是否满足收敛标准就可以。
[0119] 步骤1010：判断所述当前目标函数值是否满足收敛标准，如果是，执行步骤1011；否则，返回步骤103，再次获取每炮对应的模拟地震记录数据，开始下一次迭代。具体需要判断所述当前目标函数值Ei是否满足 其中Ei为所述当前目标函数值，
Ei-1为第i-1次的目标函数值；Theshold为设定的阈值，通常选取1.0e-5。
[0120] 步骤1011：将当前偏移剖面确定为偏移成像最终的偏移剖面。
[0121] 图2为本发明地震资料偏移成像系统的结构示意图。如图2所示，所述系统包括：
[0122] 初始化模块201，用于获取待进行偏移成像的地震资料数据、观测系统参数、偏移速度模型参数、偏移密度模型参数和偏移参数；所述地震资料数据为多炮的观测地震记录数据；
[0123] 震源波场获取模块202，用于获取每炮对应的震源波场；所述每炮对应的震源波场包括每炮的所有记录时刻的震源波场；
[0124] 模拟地震记录数据获取模块203，用于获取当前每炮对应的模拟地震记录数据，所述模拟地震记录数据根据前一次每炮对应的模拟地震记录数据利用非线性反偏移算子获得；
[0125] 伴随波场获取模块204，用于根据当前每炮对应的观测地震记录数据和对应的模拟地震记录数据，获取每炮的每一个记录时刻的伴随波场；
[0126] 第一梯度剖面获取模块205，用于根据每炮对应的震源波场和伴随波场，获得当前每炮对应的第一梯度剖面；所述第一梯度剖面为每炮对应的单炮梯度剖面；
[0127] 第二梯度剖面获取模块206，用于根据所述观测系统参数和所述偏移参数，将所有炮的第一梯度剖面叠加，获得当前第二梯度剖面；
[0128] 下降方向剖面获取模块207，用于利用最优化反演算法获得当前下降方向剖面；
[0129] 最优化步长获取模块208，用于利用最优化步长计算方法获取当前最优化步长；
[0130] 偏移剖面获取模块209，用于根据所述当前最优化步长和所述当前下降方向剖面获得当前偏移剖面；
[0131] 目标函数值获取模块2010，用于根据 获得当前目标函数值E；其中δd表示当前地震记录数据的残差；T表示转置运算；
[0132] 收敛判断模块2011，用于判断所述当前目标函数值是否满足收敛标准，得到第一判断结果；
[0133] 最终的偏移剖面确定模块2012，用于当所述第一判断结果表示所述当前目标函数值满足收敛标准时，将所述当前偏移剖面确定为偏移成像最终的偏移剖面；
[0134] 返回模块2013，用于当所述第一判断结果表示所述当前目标函数值不满足收敛标准时，返回获取每炮对应的模拟地震记录数据步骤，再次获取每炮对应的模拟地震记录数据。
[0135] 具体实施方式一：
[0136] 图3是本实施方式提供的二维洼陷偏移速度模型图。在此模型上设置49个爆炸震源，震源子波设定为雷克子波，主频为15赫兹，起始震源点的位于(150m，100m)处，炮间隔为100m。采用中间放炮两边接收观测系统，单边最大偏移距2300m，最小偏移距为150m，道间距为10m。
[0137] 图4是图3所示二维洼陷模型的多炮叠加偏移剖面，其中，图4(a)是利用逆时偏移方法所得的偏移剖面，图4(b)是利用本发明方法所得的偏移剖面。如图4所示，从图4(a)中可以看出，剖面存在比较明显的噪声，剖面的分辨率较低，振幅不均衡。从图4(b)中可以看出，本发明方法所得的偏移剖面具有很高的精度、分辨率和信噪比，而且振幅是均衡性很好，也证明了本发明方法的可行性及有效性。
[0138] 具体实施方式二：
[0139] 图5是本实施方式提供的Marmousi-2偏移速度模型；该模型是验证各种偏移方法成像效果的国际标准模型之一。在此模型上设置107个爆炸震源，震源子波设定为雷克子波，主频为20赫兹，起始震源点的位于(20m，20m)处，炮间隔为250m。采用中间放炮两边接收观测系统，单边最大偏移距2680m，道间距为10m。
[0140] 图6是图5所示Marmousi-2模型的多炮叠加偏移剖面：其中，图6(a)是利用逆时偏移方法所得的偏移剖面；图6(b)是利用本发明方法所得的偏移剖面。由图6也可以看出，传统方法的成像剖面振幅严重失衡，深部剖面相对振幅并不保真。而本发明方法所得成像剖面，效果更好，分辨率和精度更高，信噪比也更好，振幅均衡性更佳。
[0141] 本发明包括：确定时间域非线性反偏移算子；确定时间域伴随方程；构建梯度公式；建立目标函数和反演流程；确定最优化步长计算方法。本发明将时间域非线性Born近似理论引入最小二乘偏移成像中，与常规最小二乘偏移成像相比，可有效降低常规最小二乘偏移成像方法对偏移模型的依赖程度，获得高精度、高分辨率、高信噪比、振幅保真的叠前深度偏移剖面。该方法可应用到各种复杂地质条件的地震资料偏移中，模型依赖性弱，成像精度和分辨率高，剖面意义明确，便于后期地质解译。
[0142] 本说明书中各个实施例采用递进的方式描述，每个实施例重点说明的都是与其他实施例的不同之处，各个实施例之间相同相似部分互相参见即可。对于实施例公开的系统而言，由于其与实施例公开的方法相对应，所以描述的比较简单，相关之处参见方法部分说明即可。
[0143] 本文中应用了具体个例对本发明的原理及实施方式进行了阐述，以上实施例的说明只是用于帮助理解本发明的方法及其核心思想；同时，对于本领域的一般技术人员，依据本发明的思想，在具体实施方式及应用范围上均会有改变之处。综上所述，本说明书内容不应理解为对本发明的限制。
