本发明涉及一种石油地震勘探数据高效率、低存储、高精度成像领域，特别是关于一种最小二乘逆时偏移梯度更新方法。
背景技术：
随着勘探目标越来越复杂，对于成像结果的精确性、分辨率和保幅性也提出了更高的要求。常规的一些偏移方法，例如逆时偏移、kirchhoff偏移等，其精确性通常会受到地震数据噪音、不规则采样、有限偏移距等影响。
最小二乘偏移(leastsquaremigration，lsm)作为一个反演问题不断的迭代拟合正演模拟的地震记录和实际观测的地震记录，隐含考虑了hessian矩阵对振幅的影响，可以减少常规偏移过程中产生的假象，具有更高的成像精度以及保幅性等优势，可以获得更加精确的反演目标。tarantola(1984)理论推导了二范数意义下波场残差最小二乘目标函数，通过目标函数的下降，不断的迭代反演成像剖面，证明了常规的偏移结果只是最小二乘偏移的第一步迭代。nemethetal.(1999)对不完整数据进行了最小二乘kirchhoff偏移并取得了很好的效果，进一步展示了最小二乘偏移的优势。目前，最小二乘逆时偏移(leastsquarereversetimemigration，lsrtm)的研究也取得了很大进展。但是，当速度体存在强反射界面或者地震记录含有大量折射波时，基于伴随状态法求取的目标函数梯度会产生很强的低频噪音，在迭代过程中会干扰成像精度、降低收敛速度。
技术实现要素：
针对伴随状态法求取目标函数梯度，收敛速度和计算效率受低频噪音影响的问题，本发明的目的是提供一种最小二乘逆时偏移梯度更新方法，该方法可以加快目标泛函梯度收敛速度，提高偏移剖面的更新效率。
为实现上述目的，本发明采取以下技术方案：一种最小二乘逆时偏移梯度更新方法，其特征在于包括以下步骤：1)根据真实速度模型得到初始速度模型，建立二范数意义下波场残差目标泛函；2)采用基于cuda计算平台，cpu/gpu协同并行加速方法正演模拟声波方程，存储震源边界波场和最后时刻的全部区域波场；3)利用存储的边界波场和最后时刻的全部区域波场逆时重建震源波场，利用伴随状态法逆时反传残差波场获得伴随状态变量，利用逆散射条件更新单炮目标泛函相对于模型参数的梯度，判断是否满足多炮循环条件，如果是，重复执行步骤3)并累加单炮目标泛函梯度，如果否，利用目标泛函梯度计算共轭梯度方向；4)基于线性born近似正演模拟震源波场和共轭梯度波场数据，判断是否满足多炮循环条件，如果是，重复步骤4)并累加每炮共轭梯度波场数据，如果否，计算迭代步长；5)更新偏移剖面和波场残差，进入下一步迭代更新；6)判断是否满足迭代条件，如果否，终止迭代并输出迭代剖面，如果是，返回步骤2)迭代更新。
进一步，所述步骤1)中，波场残差目标泛函建立过程如下：采用高斯函数对真实速度模型进行平滑处理，得到初始平滑速度模型；波场残差为反偏移方程数值模拟的地震记录和野外观测的地震记录振幅差值，二范数意义下波场残差目标泛函如下：
其中，dsyn为反偏移的理论数据，dobs为观测的去除直达波的地震记录，ng为炮点的数目，ns为检波点的数目，xr为检波点，xs为震源，w角频率。
进一步，所述步骤2)中，采用基于cuda计算平台，cpu/gpu协同并行加速技术，把数值模拟计算区域分割成多个16×16网格节点的子区域gpu显存，每个子区域开存储空间为28×28网格节点的共享存储器模拟整个区域波场，每个时刻把子区域数据拷贝到共享存储器中，同时拷贝相邻的子区域6个网格节点数据到该共享存储器相应的位置，把该时刻共享存储器模拟出来的波场拷贝到gpu显存，同时保存该时刻的计算区域四周边界向内部计算区域连续6个网格节点波场到计算机内存中，每一个时刻重复着上述的步骤直达最大时刻，保存最大时刻的全部计算区域波场值到gpu显存中。
进一步，所述步骤3)中，具体处理过程如下：3.1)采用保存的边界波场以及最大时刻全部波场逆时重建震源波场，把每一个时刻的cpu端的边界波场拷贝到gpu显存相应的子区域边界波场处，再用步骤3)相同的方式把子区域数值拷贝到共享存储器上，利用波动方程时间可逆性，重建最大时刻到零时刻每一个历史时刻的内部计算区域的震源波场；3.2)利用伴随状态方法，反传残差波场获得每一个时刻伴随状态变量，利用逆散射成像条件求取单炮目标泛函相对于模型参数的梯度；3.3)根据目标泛函相对于模型参数的梯度，获得第k次迭代共轭梯度dk。
进一步，所述步骤3.2)中，第k次迭代梯度(逆散射成像条件)更新表达式如下：
其中，λk为第k次迭代伴随状态变量，v0(x)为背景参考速度模型，ω是圆周频率，s(ω)为震源函数，g0(xr；x,ω)为散射点x到检波点xr的格林函数，g0(x；xs,ω)为从震源xs到散射点x的格林函数，▽为关于空间步长的梯度算子，ng,ns分别代表炮点和检波点的数目。
进一步，所述步骤3.3)中，第k次迭代共轭梯度dk的计算公式为：
其中，k为迭代次数；βk为共轭梯度系数，▽j(mk)为第k次迭代目标泛函梯度。
进一步，所述步骤4)中，迭代步长具体计算过程如下：基于声波方程线性born近似模拟计算反射地震数据，第k次迭代单炮反射波数据lmk如下：
lmk＝ω2∫g0(xr；x,ω)mk(x)g0(x；xs,ω)s(ω)dx
同理计算单炮共轭梯度波场数据ldk：
ldk＝ω2∫g0(xr；x,ω)dk(x)g0(x；xs,ω)s(ω)dx
式中，xr为检波点位置、xs为震源位置，mk为第k次迭代反射率模型，dk第k次迭代共轭梯度方向，s(ω)为震源子波，g0(xr；x,ω)为散射点x到检波点xr的格林函数，g0(x；xs,ω)为从震源xs到散射点x的格林函数，l为线性born正演算子如下：
l＝ω2∫g0(x；x',ω)g0(x'；xs,ω)s(ω)dx'
累加迭代步长ak表达式分子和分母，累加所有炮共轭梯度波场数据和残差数据的乘积和累加所有炮共轭梯度波场数据和共轭梯度波场数据乘积，得到第k次迭代步长ak。
进一步，所述第k次迭代步长ak为：
进一步，所述步骤5)中，根据获得的第k次迭代步长ak，利用第k次迭代共轭梯度方向dk和当前迭代反射率模型mk，计算更新偏移剖面，偏移剖面的更新表达式为：
mk+1＝mk+akdk。
进一步，所述步骤6)中，给定最大迭代次数终止条件n，当迭代次数k小于迭代终止条件n，继续迭代，当迭代次数k大于等于迭代终止条件n，输出最终的迭代剖面成像结果。
本发明由于采取以上技术方案，其具有以下优点：与常规互相关成像条件最小二乘逆时偏移相比，本发明高效的最小二乘逆时偏移梯度更新方法，考虑伴随状态法目标函数梯度会受到低频噪音干扰收敛速度和计算效率，引入线性born逆散射成像条件应用于最小二乘逆时偏移方法中，可以有效的压制低频噪音对最小二乘逆时偏移成像质量的干扰，加快收敛速度，提高剖面的更新效率。同时cpu/gpu协同并行加速技术以及共享存储器的引入可以提高计算效率。
附图说明
图1是本发明的整体流程示意图；
图2是本发明的salt2d反演剖面结果示意图，其中，图a是salt2d速度模型，图b是第一次迭代结果(rtm),图c是本发明20次迭代结果(lsrtm)。
具体实施方式
下面结合附图和实施例对本发明进行详细的描述。
如图1所示，本发明提供一种最小二乘逆时偏移梯度更新方法，其包括以下步骤：
1)根据真实速度模型得到较为准确的初始速度模型，建立二范数意义下波场残差目标泛函；
建立过程如下：
采用高斯函数对真实速度模型进行平滑处理，得到较为准确的初始平滑速度模型。波场残差为反偏移方程数值模拟的地震记录和野外观测的地震记录振幅差值，二范数意义下波场残差目标泛函如下：
其中，dsyn为反偏移的理论数据，dobs为观测的去除直达波的地震记录，ng为炮点的数目，ns为检波点的数目，xr为检波点，xs为震源，w为角频率。
2)正演存储边界波场：采用基于cuda计算平台，cpu/gpu协同并行加速方法正演模拟声波方程，存储震源边界波场和最后时刻的全部区域波场；
具体如下：
采用空间12阶精度，时间2阶精度，完全匹配层(pml)吸收边界条件数值模拟声波方程，则频域常密度声波方程为：
其中，v(x)为速度模型，ω是圆周频率，p(x,ω)为频率压力波场，s(ω)为震源函数，▽为关于空间步长的梯度算子，δ(x-xs)为狄拉克δ函数。
采用基于cuda计算平台，cpu/gpu协同并行加速技术，把数值模拟计算区域分割成多个16×16网格节点的子区域(block)gpu显存，每个子区域开存储空间为28×28网格节点的共享存储器模拟整个区域波场，每个时刻把子区域(block)数据拷贝到共享存储器(share-memory)中，同时拷贝相邻的子区域(block)6个网格节点数据到该共享存储器相应的位置，保证该共享存储器中间内部区域16×16网格节点数据可以利用有限差分数值模拟出来，把该时刻共享存储器模拟出来的波场拷贝到gpu显存，同时保存该时刻的计算区域四周边界向内部计算区域连续6个网格节点波场到计算机内存中，每一个时刻重复着上述的步骤直达最大时刻，保存最大时刻的全部计算区域波场值到gpu显存中。
3)利用存储的边界波场和最后时刻的全部区域波场逆时重建震源波场，利用伴随状态法逆时反传残差波场获得伴随状态变量，利用逆散射条件更新单炮目标泛函相对于模型参数的梯度，判断炮数is(is为循环的炮数)是否满足多炮循环终止条件ns(ns为循环的总炮数)，如果满足循环条件，重复执行步骤3)并累加单炮目标泛函梯度；如果不满足循环条件，则炮循环终止；利用目标泛函梯度计算共轭梯度方向；
具体如下：
3.1)采用步骤2)保存的边界波场以及最大时刻全部波场逆时重建震源波场，把每一个时刻的cpu端的边界波场拷贝到gpu显存相应的子区域(block)边界波场处，再用步骤2)相同的方式把子区域数值拷贝到共享存储器上，利用波动方程时间可逆性，重建最大时刻到零时刻每一个历史时刻的内部计算区域的震源波场。
3.2)利用伴随状态方法，反传残差波场获得每一个时刻伴随状态变量，利用逆散射成像条件求取单炮目标泛函相对于模型参数的梯度，第k次迭代梯度更新表达式如下：
其中，λk为第k次迭代伴随状态变量，v0(x)为背景参考速度模型，ω是圆周频率，s(ω)为震源函数，g0(xr；x,ω)为散射点x到检波点xr的格林函数，g0(x；xs,ω)为从震源xs到散射点x的格林函数，▽为关于空间步长的梯度算子，ng,ns分别代表炮点和检波点的数目。
3.3)根据目标泛函相对于模型参数的梯度，利用以下公式获得第k次迭代共轭梯度dk：
其中，k为迭代次数；βk为共轭梯度系数，▽j(mk)为第k次迭代目标泛函梯度。
4)基于线性born近似正演模拟震源波场和共轭梯度波场数据，判断炮数is(is为循环的炮数)是否满足多炮循环终止条件ns(ns为循环的总炮数)，如果满足循环条件，重复步骤4)并累加每炮共轭梯度波场数据；如果不满足循环条件，则终止炮循环，进行迭代步长的计算；
迭代步长具体计算过程如下：
基于声波方程线性born近似模拟计算反射地震数据(也叫反偏移方程)，第k次迭代单炮反射波数据lmk如下：
lmk＝ω2∫g0(xr；x,ω)mk(x)g0(x；xs,ω)s(ω)dx
同理计算单炮共轭梯度波场数据ldk表达式如下：
ldk＝ω2∫g0(xr；x,ω)dk(x)g0(x；xs,ω)s(ω)dx
式中：xr为检波点位置、xs为震源位置，mk为第k次迭代反射率模型，dk第k次迭代共轭梯度方向，s(ω)为震源子波，g0(xr；x,ω)为散射点x到检波点xr的格林函数，g0(x；xs,ω)为从震源xs到散射点x的格林函数，l为线性born正演算子(反偏移算子)如下：
l＝ω2∫g0(x；x',ω)g0(x'；xs,ω)s(ω)dx'
累加迭代步长ak表达式分子和分母，累加所有炮共轭梯度波场数据和残差数据的乘积(步长表达式分子)和累加所有炮共轭梯度波场数据和共轭梯度波场数据乘积(步长表达式分母)，第k次迭代步长ak表达式如下：
5)更新偏移剖面，进入步骤6)迭代更新。
具体如下：
根据步骤4)获得的第k次迭代步长ak，利用第k次迭代共轭梯度方向dk和当前迭代反射率模型mk，计算更新偏移剖面。偏移剖面的更新表达式为：
mk+1＝mk+akdk。
6)判断是否满足迭代终止条件，如果是，终止迭代并输出迭代剖面，如果否，返回步骤2)迭代更新；
具体如下：
给定最大迭代次数n终止条件，当迭代次数k小于迭代终止条件n，继续迭代，当迭代次数k大于等于迭代次数n，输出最终的迭代剖面成像结果。
如图2所示，为本发明对salt2d速度模型采用最小二乘逆时偏移得到的偏移剖面，图a为salt2d速度模型。图b为逆散射条件逆时偏移剖面，可以发现低频被很好地压制。图c为逆散射条件，最大迭代次数n＝20的最小二乘逆时偏移剖面。
综上所述，目前广泛使用的伴随状态法求取目标泛函梯度会产生低频噪音，干扰迭代收敛速度和计算效率。本发明利用线性born逆散射成像条件迭代更新最小二乘逆时偏移梯度，可以有效的压制低频噪音对最小二乘逆时偏移成像质量的干扰，加快收敛速度，提高剖面的更新效率。同时cpu/gpu协同并行加速技术以及共享存储器的引入可以提高计算效率。
上述各实施例仅用于说明本发明，各个步骤都是可以有所变化的，在本发明技术方案的基础上，凡根据本发明原理对个别步骤进行的改进和等同变换，均不应排除在本发明的保护范围之外。