具体实施方式

为了使本发明目的、技术方案更加清楚明白，下面结合实施例，对本发明作进一步详细说明。

实施例1：

如图1所示的流程图，本发明所述的地震资料成像域最小二乘逆时偏移成像方法，包括如下步骤：

步骤一：对原始共炮域地震资料进行预处理，预处理包括剔除直达波和压制噪声；

步骤二：利用规则网格有限差分方法计算震源正向传播波场；

步骤三：利用规则网格有限差分方法计算检波点反向传播波场；

步骤四：利用互相关成像条件计算逆时偏移成像剖面；

步骤五：基于偏移速度模型，在空间坐标为(x_i，z_j)的离散点处设置速度扰动，空间坐标(x_i，z_j)的计算公式如下：

其中，i表示x方向离散点的序号，j表示z方向离散点的序号，Δx表示x方向离散点的间距，Δz表示z方向离散点的间距，n_x表示x方向离散点的个数，n_z表示z方向离散点的个数；

步骤六：计算步骤五中离散点对应的点扩散函数；

步骤七：利用步骤六计算的点扩散函数，采用插值方法计算所有离散点的点扩散函数；

步骤八：计算Hessian矩阵，Hessian矩阵的计算公式如下：

H＝[psf₁，…，psf_np]， (2)

其中，psf表示点扩散函数，np表示点扩散函数的总个数；

步骤九：求解最小二乘反演方程，获得最小二乘逆时偏移成像剖面。

其中，所述步骤五中，点扩散函数的计算公式如下：

其中，psf表示点扩散函数，ip表示点扩散函数的序号，W表示剔除与第ip个点扩散函数无关的偏移噪声的算子，is表示震源序号，ns表示震源的总个数，t表示时间，T表示地震资料时间采样长度，P表示偏移速度模型存在速度扰动时的震源正向传播波场，Q表示偏移速度模型存在速度扰动时的检波点反向传播波场。

其中，所述步骤九中，最小二乘反演方程如下：

Hm＝m_mig， (4)

其中，H表示步骤八计算的Hessian矩阵，m表示最小二乘逆时偏移成像剖面，m_mig表示步骤四计算的逆时偏移成像剖面。

实施例2：

下面结合具体实施方式，对于本发明进行解释和说明。

为了进一步说明本方法的实现思路及实现过程并证明方法的有效性，用一个理论速度模型合成资料进行测试。理论速度模型如图2所示。

S1:观测系统：炮间距100m，共200炮。拖缆长6km，道间距10m，共600道。拖缆位于震源的左侧。采用时间二阶、空间十二阶精度的规则网格有限差分方法正演模拟合成地震资料。对合成地震资料进行预处理，剔除直达波。

S2:利用规则网格有限差分方法计算震源正向传播波场。

S3:利用规则网格有限差分方法计算检波点反向传播波场。

S4:利用互相关成像条件计算逆时偏移成像剖面，如图3所示。

S5:利用公式(1)，基于偏移速度模型设置速度扰动。

S6:利用公式(3)计算点扩散函数，如图4所示。

S7:利用步骤六计算的点扩散函数，采用插值方法计算所有离散点的点扩散函数。

S8:利用公式(2)计算Hessian矩阵。

S9：利用公式(4)计算最小二乘逆时偏移成像剖面，如图5所示。

图6是理论反射系数模型图。将图3和图5与图6对比可以看出，相比于逆时偏移成像剖面，利用本发明计算的最小二乘逆时偏移成像剖面更接近理论反射系数模型。本发明计算的最小二乘逆时偏移成像剖面与理论反射系数模型非常相似。这证明了本发明的正确性和有效性。

本发明所述的地震资料成像域最小二乘逆时偏移成像方法，考虑到了常规数据域最小二乘逆时偏移成像方法计算量、极其耗时的缺陷，通过计算点扩散函数来估计Hessian矩阵，解决了常规数据域最小二乘逆时偏移成像方法需要通过多次偏移和反偏移运算估计Hessian矩阵极其耗费时间的问题；本发明计算简单、易于实现、耗时少，具有很强的适应性，且最小二乘逆时偏移结果可靠性高。

以上所述仅为本发明的较佳实施例而己，并不以本发明为限制，凡在本发明的精神和原则之内所作的均等修改、等同替换和改进等，均应包含在本发明的专利涵盖范围内。