一种最小二乘逆时偏移梯度更新方法技术

本发明专利技术涉及一种最小二乘逆时偏移梯度更新方法，其步骤：根据初始速度模型建立二范数意义下波场残差目标泛函；正演模拟声波方程，存储震源边界波场和最后时刻的全部区域波场；逆时重建震源波场，反传残差波场获得伴随状态变量，利用逆散射条件更新单炮目标泛函相对于模型参数的梯度，判断是否满足多炮循环条件，如果是，重复执行并累加单炮目标泛函梯度，如果否，利用目标泛函梯度计算共轭梯度方向；正演模拟震源波场和共轭梯度波场数据，判断是否满足多炮循环条件，如果是，重复并累加每炮共轭梯度波场数据，如果否，计算迭代步长；更新偏移剖面，迭代更新；判断是否满足迭代条件，如果否，终止迭代并输出迭代剖面，如果是则返回正演模拟存储边界波场。

【技术实现步骤摘要】
一种最小二乘逆时偏移梯度更新方法
本专利技术涉及一种石油地震勘探数据高效率、低存储、高精度成像领域，特别是关于一种最小二乘逆时偏移梯度更新方法。
技术介绍
随着勘探目标越来越复杂，对于成像结果的精确性、分辨率和保幅性也提出了更高的要求。常规的一些偏移方法，例如逆时偏移、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)根据真实速度模型得到初始速度模型，建立二范数意义下波场残差目标泛函；2)采用基于CUDA计算平台，CPU/GPU协同并行加速方法正演模拟声波方程，存储震源边界波场和最后时刻的全部区域波场；3)利用存储的边界波场和最后时刻的全部区域波场逆时重建震源波场，利用伴随状态法逆时反传残差波场获得伴随状态变量，利用逆散射条件更新单炮目标泛函相对于模型参数的梯度，判断是否满足多炮循环条件，如果是，重复执行步骤3)并累加单炮目标泛函梯度，如果否，利用目标泛函梯度计算共轭梯度方向；4)基于线性Born近似正演模拟震源波场和共轭梯度波场数据，判断是否满足多炮循环条件，如果是，重复步骤4)并累加每炮共轭梯度波场数据，如果否，计算迭代步长；5)更新偏移剖面和波场残差，进入下一步迭代更新；6)判断是否满足迭代条件，如果否，终止迭代并输出迭代剖面，如果是，返回步骤2)迭代更新。

【技术特征摘要】
1.一种最小二乘逆时偏移梯度更新方法，其特征在于包括以下步骤：1)根据真实速度模型得到初始速度模型，建立二范数意义下波场残差目标泛函；2)采用基于CUDA计算平台，CPU/GPU协同并行加速方法正演模拟声波方程，存储震源边界波场和最后时刻的全部区域波场；3)利用存储的边界波场和最后时刻的全部区域波场逆时重建震源波场，利用伴随状态法逆时反传残差波场获得伴随状态变量，利用逆散射条件更新单炮目标泛函相对于模型参数的梯度，判断是否满足多炮循环条件，如果是，重复执行步骤3)并累加单炮目标泛函梯度，如果否，利用目标泛函梯度计算共轭梯度方向；4)基于线性Born近似正演模拟震源波场和共轭梯度波场数据，判断是否满足多炮循环条件，如果是，重复步骤4)并累加每炮共轭梯度波场数据，如果否，计算迭代步长；5)更新偏移剖面和波场残差，进入下一步迭代更新；6)判断是否满足迭代条件，如果否，终止迭代并输出迭代剖面，如果是，返回步骤2)迭代更新。2.如权利要求1所述的一种最小二乘逆时偏移梯度更新方法，其特征在于：所述步骤1)中，波场残差目标泛函建立过程如下：采用高斯函数对真实速度模型进行平滑处理，得到初始平滑速度模型；波场残差为反偏移方程数值模拟的地震记录和野外观测的地震记录振幅差值，二范数意义下波场残差目标泛函如下：其中，dsyn为反偏移的理论数据，dobs为观测的去除直达波的地震记录，ng为炮点的数目，ns为检波点的数目，xr为检波点，xs为震源，w角频率。3.如权利要求1所述的一种最小二乘逆时偏移梯度更新方法，其特征在于：所述步骤2)中，采用基于CUDA计算平台，CPU/GPU协同并行加速技术，把数值模拟计算区域分割成多个16×16网格节点的子区域GPU显存，每个子区域开存储空间为28×28网格节点的共享存储器模拟整个区域波场，每个时刻把子区域数据拷贝到共享存储器中，同时拷贝相邻的子区域6个网格节点数据到该共享存储器相应的位置，把该时刻共享存储器模拟出来的波场拷贝到GPU显存，同时保存该时刻的计算区域四周边界向内部计算区域连续6个网格节点波场到计算机内存中，每一个时刻重复着上述的步骤直达最大时刻，保存最大时刻的全部计算区域波场值到GPU显存中。4.如权利要求1所述的一种最小二乘逆时偏移梯度更新方法，其特征在于：所述步骤3)中，具体处理过程如下：3.1)采用保存的边界波场以及最大时刻全部波场逆时重建震源波场，把每一个时刻的CPU端的边界波场拷贝到GPU显存相应的子区域边界波场处，再用步骤3)相同的方式把子区域数值拷贝到共享存储器上，利用波动方程时间可逆性，重建最大时刻到零时刻每一个历史时刻的内部计算区域的震源波场；3.2)利用伴随状态方法，反传残差波场获得每一个时刻伴随状态变量，利用逆散射成像条件求取单炮目标泛函相对于模型参数的梯度；3.3)根据目标泛函相对于模型参数的梯度，获得第k次迭代共轭梯度dk。5.如权利要求4所述的一种最小二乘逆时偏移梯度更新方法，其特征在于：所述步骤3.2)中，第k次迭代梯度(逆散射成像...

【专利技术属性】
技术研发人员：方修政，吴迪，钮凤林，
申请(专利权)人：中国石油大学北京，
类型：发明
国别省市：北京,11