一种混合域地震偏移海森矩阵估计方法
阅读说明：本技术 一种混合域地震偏移海森矩阵估计方法 (Mixed domain seismic migration hessian matrix estimation method ) 是由 杨继东 张东林 黄建平 李振春 徐洁 秦善源 于 2021-08-23 设计创作，主要内容包括：本发明公开了一种混合域地震偏移海森矩阵估计方法,涉及勘探地球物理技术领域,该方法首先利用常规偏移方法获得初始成像结果,然后将该偏移结果作为反射率模型,进行线性正演和偏移获得二次成像结果,使用这两个偏移结果在空间-波数域估计非稳相滤波器,用于近似汉森矩阵,并利用混合域反褶积计算汉森矩阵的逆,再使用估计的汉森矩阵的逆在最小二乘偏移中作为泛函的梯度,以提高反演的收敛速率,使得最小二乘反演在4-5次迭代内收敛至准确解。本发明的有益效果是,该混合域地震偏移汉森矩阵估计方法既具有较高的计算效率,又可加快最小二乘偏移收敛速度,降低最小二乘反射率反演的计算成本,具有较好的应用前景。(The invention discloses a mixed domain seismic migration hessian matrix estimation method, which relates to the technical field of exploration geophysics, and is characterized in that an initial imaging result is obtained by utilizing a conventional migration method, then the migration result is used as a reflectivity model, linear forward modeling and migration are carried out to obtain a secondary imaging result, the two migration results are used for estimating a non-stationary phase filter in a space-wave number domain to approximate a hessian matrix, the inverse of the hessian matrix is calculated by utilizing mixed domain deconvolution, and the estimated inverse of the hessian matrix is used as a functional gradient in least square migration to improve the convergence rate of inversion, so that the least square inversion converges to an accurate solution in 4-5 iterations. The mixed domain seismic migration Hansen matrix estimation method has the advantages of being high in calculation efficiency, capable of accelerating the least square migration convergence speed, reducing the calculation cost of least square reflectivity inversion and good in application prospect.)
技术领域
本发明涉及勘探地球物理技术领域，尤其涉及一种混合域地震偏移海森矩阵估计方法。
背景技术
目前，高精度地震偏移成像是油气勘探中重要的关键技术之一，常规地震偏移方法，包括：Kirchhoff偏移、高斯束偏移、单程波偏移以及逆时偏移等，在数学意义上均可看成线性正演模拟的伴随算子，而非真正意义的逆算子。因此，当地震采集不规则、震源及检波器频带较窄以及地下照明不均匀时，基于伴随算子的常规偏移方法很难获得高质量成像结果。基于反演框架的最小二乘偏移通过数据拟合可以逐渐逼近正演算子的逆，进而大幅提高偏移结果的振幅均衡性和空间分辨率，是目前地震成像领域研究的热点。
最小二乘偏移成像的核心在于如何高效计算汉森矩阵并求其逆算子。常规数据域最小二乘偏移首先构建观测数据和预测数据的误差泛函，然后使用线性正演计算预测数据，使用伴随偏移计算泛函梯度，并通过迭代的方式逐步逼近汉森矩阵的逆。该方法通常需要数十次至上百次正演模拟和伴随偏移，计算效率较低，往往很难应用于大规模实际生产处理中。通过设计合理的预条件算子，可以加快最小二乘偏移的收敛速率，但目前常用的震源照明补偿只能改善振幅的相对强弱关系，不能够提高空间分辨率。
为使最小二乘偏移走向生产化应用，急需研发精确高效的预条件算子，使其使用较少次数的迭代就可以收敛至真实解附近。
发明内容
为解决上述最小二乘偏移成像中的收敛速度较慢、计算效率较低的技术问题，本发明公开了一种混合域地震偏移海森矩阵估计方法。
为实现上述目的，本发明采用下述技术方案：
一种混合域地震偏移海森矩阵估计方法，具体包括以下步骤：
步骤一，获取输入数据；
步骤二，通过求解波动方程计算震源波场，计算检波点波场，再采用零延迟互相关成像条件计算初步成像结果；
步骤三，将初步成像结果作为反射率模型，使用线性正演模拟计算合成地震数据；
步骤四，使用步骤三合成的地震数据和偏移速度模型，采用步骤二的成像处理过程进行再次偏移，获得二次成像结果；
步骤五，将初步成像结果和二次成像结果应用空间-波数域变换，获得混合域成像结果；
步骤六，计算近似汉森矩阵的非稳相混合域滤波器；
步骤七，使用混合域反褶积计算汉森矩阵的逆；
步骤八，将估计的汉森矩阵的逆作为预条件算子应用于最小二乘偏移的泛函梯度中，以提高反演的收敛速率；
步骤九，在每次最小二乘迭代中都应用步骤八中的预条件算子，使得最小二乘偏移快速收敛于准确解。
进一步地，步骤一中，输入数据包括：偏移速度模型v(x)、震源函数子波f(t)、野外观测数据dobs(xS,xR,t)，混合域滤波器窗口长度l，xR为检波点位置。
进一步地，步骤二包括：
(2.1)使用震源函数子波f(t)和偏移速度模型v(x)，通过求解波动方程计算震源波场US(x,t)，波动方程如下：
其中，x为地下成像空间位置，xS代表了震源位置，t代表时间变量，δ(x-xs)为狄拉克delta函数。
(2.2)使用野外观测数据dobs(xS,xr,t)和偏移速度模型v(x)，计算检波点波场UR(x,t)，具体表达式如下：
(2.3)采用零延迟互相关成像条件计算初步成像结果m0(x)，且互相关成像条件为：
其中，T为地震记录的长度，US(x,t)为震源波场，为UR(x,t)检波点波场。
进一步地，步骤三中，将初步成像结果m0(x)作为反射率模型，计算合成地震数据dsyn(xS,xr,t)所采用的线性正演模拟表达式如下：
其中，U0(x,t)为背景波场，U1(x,t)为扰动波场。
进一步地，步骤五中，对初次成像结果m0(x)和二次成像结果m1(x)应用空间-波数域变换，该空间-波数域变换表达式为：
其中，x＝[x,y,z]T为空间矢量，k＝[kx,ky,kz]T为波数矢量，x'＝[x',y',z']T为空间-波数域变换后的空间矢量，M0(x',k)、M1(x',k)分别为m0(x)、m1(x)的空间-波数域变换结果，σx、σy和σz分别为在x、y和z方向的窗口长度，设置σx、σy和σz的值为输入混合域滤波器窗口长度l。
进一步地，步骤六中，使用M0(x',k)和M1(x',k)，计算近似汉森矩阵的非稳相混合域滤波器F(x',k)，其表达式如下：
其中，ε为避免除以零的小值，取ε＝0.001*max|M1(x',k)|。
进一步地，步骤七中，汉森矩阵的逆Hinv(x',k)为：
进一步地，步骤八中，预条件算子如下：
其中，g(x',k)为泛函梯度的空间-波数域谱，gpre(x)为应用预条件后的梯度。
本发明的有益效果是，首先计算了均匀介质中解析的汉森矩阵，并通过分析发现当把该矩阵的列向量投影至成像空间时所形成的点拓展函数具有局部化特征，基于该分析，本发明专利提出了一种混合域地震偏移汉森矩阵估计方法，即根据地震偏移的结果构建一系列空间-波数域的滤波器，用于近似局部化的点拓展函数，进而可以通过混合域反褶积计算汉森矩阵的逆。再将该汉森矩阵的逆作为最小二乘偏移的预条件算子，可大幅提高反演的收敛速率，使得最小二乘反演在4-5次迭代内收敛至准确解。该混合域地震偏移汉森矩阵估计方法既具有较高的计算效率，又可加快最小二乘偏移收敛速度，降低最小二乘反射率反演的计算成本，具有较好的应用前景。
附图说明
图1为部分Sigsbee速度模型图；
图2为由传统偏移方法获得的成像结果m0(x)；
图3为由m0(x)作为反射率模型，进行线性正演和再次偏移所获得的成像结果m1(x)；
图4为对m0(x)应用汉森矩阵的逆获得的校正后的偏移结果；
图5为常规最小二乘偏移迭代更新5次后，获得的反演成像结果
图6为使用本发明汉森矩阵作为预条件算子，优化的最小二乘偏移迭代更新5次后，获得的反演成像结果。
具体实施方式
下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。
基于线性反演框架，最小二乘偏移可以消除由于采集不规则造成的偏移假象，并可以改善深部振幅均衡性和整体的空间分辨率。，但常规数据域最小二乘偏移通过迭代数据拟合逐步逼近汉森矩阵的逆，在实际应用中，往往需要数十次正演模拟和伴随偏移，计算量较大，在当前计算条件下仍较难走向实用化生产应用。
针对最小二乘偏移成像中的收敛速度较慢、计算效率较低的问题，本发明公开了一种混合域地震偏移海森矩阵估计方法，具体包括以下步骤：
(1)获取输入数据。
输入数据包括：偏移速度模型v(x)、震源函数子波f(t)、野外观测数据dobs(xS,xR,t)，混合域滤波器窗口长度l，xR为检波点位置。
(2)通过求解波动方程计算震源波场，计算检波点波场，再采用零延迟互相关成像条件计算初步成像结果。
具体地：
(2.1)使用震源函数子波f(t)和偏移速度模型v(x)，通过求解波动方程计算震源波场US(x,t)，波动方程如下：
其中，x为地下成像空间位置，xS代表了震源位置，t代表时间变量，δ(x-xs)为狄拉克delta函数。
(2.2)使用野外观测数据dobs(xS,xr,t)和偏移速度模型v(x)，计算检波点波场UR(x,t)，具体表达式如下：
(2.3)采用零延迟互相关成像条件计算初步成像结果m0(x)，且互相关成像条件为：
其中，T为地震记录的长度，US(x,t)为震源波场，为UR(x,t)检波点波场。
(3)将初步成像结果作为反射率模型，使用线性正演模拟计算合成地震数据。
将初步成像结果m0(x)作为反射率模型，计算合成地震数据dsyn(xS,xr,t)所采用的线性正演模拟表达式如下：
其中，U0(x,t)为背景波场，U1(x,t)为扰动波场。
(4)使用步骤(3)合成的地震数据dsyn(xS,xr,t)和偏移速度模型，采用步骤(2)的成像处理过程进行再次偏移，获得二次成像结果m1(x)；
(5)将初步成像结果和二次成像结果应用空间-波数域变换，获得混合域成像结果M0(x',k)和M1(x',k)。
对初次成像结果m0(x)和二次成像结果m1(x)应用空间-波数域变换，该空间-波数域变换表达式为：
其中，x＝[x,y,z]T为空间矢量，k＝[kx,ky,kz]T为波数矢量，x'＝[x',y',z']T为空间-波数域变换后的空间矢量，M0(x',k)、M1(x',k)分别为m0(x)、m1(x)的空间-波数域变换结果，σx、σy和σz分别为在x、y和z方向的窗口长度，本方法中设置σx、σy和σz的值为输入混合域滤波器窗口长度l。
(6)计算近似汉森矩阵的非稳相混合域滤波器。
使用M0(x',k)和M1(x',k)，计算近似汉森矩阵的非稳相混合域滤波器F(x',k)，其表达式如下：
其中，ε为避免除以零的小值，取ε＝0.001*max|M1(x',k)|。
(7)使用混合域反褶积计算汉森矩阵的逆。
汉森矩阵的逆Hinv(x',k)为：
(8)将估计的汉森矩阵的逆作为预条件算子应用于最小二乘偏移的泛函梯度中，以提高反演的收敛速率。
预条件算子如下：
其中，g(x',k)为泛函梯度的空间-波数域谱，gpre(x)为应用预条件后的梯度。
(9)在每次最小二乘迭代中都应用步骤八中的预条件算子，使得最小二乘偏移快速收敛于准确解。
将本方法应用于国际标准Sigsbee模型地震偏移成像中，取得了较好的成像结果。
图1为部分Sigsbee速度模型图，观测数据的震源函数为15Hz的雷克子波，观测数据集为50个炮集，均匀分布在地表，每个炮集包含500个检波点，间隔为10m。图2为使用观测数据获得的偏移成像结果m0(x)，可以看出，常规偏移对于浅部反射层刻画较好，但对于盐下深部沉积层，受盐丘能量屏蔽作用影响，很难获得高质量成像结果。图3为使用图2所示的成像结果m0(x)作为反射率模型，进行线性正演和偏移获得的二次偏移结果m1(x)，受汉森矩阵影响，二次偏移结果相对初始偏移结果深部振幅更弱、分辨率更差。利用图2和图3所示的偏移结果在空间-波数域估计汉森矩阵，并应用于初始结果上，获得了图4的成像结果，可以看出深部能量得到了显著补偿；此外，整体的空间分辨率也得到了很大提高。图5和图6对比了常规最小二乘和使用本发明专利提出的汉森矩阵作为预条件的最小二乘偏移在5次迭代更新后获得的成像结果。可以看出，常规最小二乘由于收敛较慢，在经过5次迭代后成像结果相对图2中的成像结果改善不大。但是，在使用本发明所述的汉森矩阵作为预条件算子后，经过5次迭代，盐丘和盐下沉积层的成像质量均得到了大幅度改善。
本发明专利提出了一种混合域地震偏移汉森矩阵估计方法，主要解决常规最小二乘偏移中收敛速率较慢、计算效率较低问题，通过将估计的汉森矩阵作为最小二乘偏移中泛函梯度的预条件算子，可使得最小二乘偏移在4-5次迭代后收敛至准确解。
本发明首先利用常规偏移方法获得初始成像结果，然后将该偏移结果作为反射率模型，进行线性正演和偏移获得二次成像结果，使用这两个偏移结果在空间-波数域估计非稳相滤波器，用于近似汉森矩阵，并利用混合域反褶积计算汉森矩阵的逆，再使用估计的汉森矩阵的逆在最小二乘偏移中作为泛函的梯度，可大题提高收敛速率，使其在4-5次迭代内收敛至准确解，最终推进最小二乘偏移的实用化进程。
当然，上述说明并非是对本发明的限制，本发明也并不仅限于上述举例，本技术领域的技术人员在本发明的实质范围内所做出的变化、改型、添加或替换，也应属于本发明的保护范围。