一种控制全声波方程反演过程的地震勘探数据处理方法

摘要

本发明涉及一种控制全声波方程反演过程的地震勘探数据处理方法，采集理论模型地震数据，获得深度域体积模量和密度初始模型；用伪谱法模拟声学介质中地震波场；计算理论模型与初始模型模拟数据的差；给定δ＝1.0e－6，计算误差能量E，当E≤δ，停止反演，输出反演结果，如果E＞δ，继续下列步骤；计算残差数据；计算体积模量和密度模型的共轭修改量；计算修改步长，获得梯度；对初始模型进行修改；将修改后的模型数据作为新的初始模型，通过对反演过程进行自动控制，使声波方程反演稳定、快速收敛，反演收敛速度比常规方法提高了近5倍。

主权项

1.一种控制全声波方程反演过程的地震勘探数据处理方法，其特征在于：具体步骤包括：(1).野外激发、接收获得地震波数据，用常规处理方法对该数据进行静校正、地表一致性振幅补偿和叠前去除噪音，获得炮集数据；(2).从步骤(1)获得的炮集数据中抽取共中心点道集数据，进行速度分析，时间域层速度计算和时间-深度转换，获得深度域体积模量和密度初始模型；①步骤(2)中速度分析为常规地震速度分析；②步骤(2)中时间域层速度用速度分析获得的均方根速度和Dix公式进行计算，$v_i^2=\frac{v_{r,i}^2{t}_i-v_{r,i-1}^2{t}_{i-1}}{t_i-t_{i-1}}---\left(1\right)$其中t_i为第i层双程旅行时，v_i为第i层的速度；③步骤(2)中时间-深度转换是按下式将时间域的层速度转换成深度域层速度，$h_i=\frac{1}{2}(t_i-t_{i-1})v_i---\left(2\right)$其中t_i为第i层双程旅行时，h_i为第i层的厚度，v_i为第i层的速度；④步骤(2)中密度初始模型用Gardner公式计算，ρ＝0.31v^0.25                  (3)其中ρ为密度，v为纵波速度。⑤步骤(2)中体积模量初始模型用下列方法计算K＝ρv²                        (4)其中K为体积模量，ρ为密度，v为纵波速度。(3).利用步骤2建立的初始体积模量和密度模型，进行地震波场正演模拟；①步骤(3)中地震波场正演模拟利用的是声学介质中的波动方程；②步骤(3)中用常规的伪谱法进行波场模拟，即对空间变量进行付氏变换，在付氏域中计算声波方程中波场对空间变量的偏导数；声波方程中波场对时间变量的偏导数用二阶中心差分法计算；③步骤(3)中震源函数采用零相位雷克子波；(4).按下式计算步骤(3)模拟的炮集数据与步骤1获得测量数据的差(残差数据)：Δd＝d_obs-d_cal                (5)其中Δd为残差数据，d_obs为测量地震数据，d_cal为计算地震数据。(5).按下式计算误差能量，$E=\frac{1}{2}{\mathrm{\Delta d}}^t{C}_D^{-1}\mathrm{\Delta d}---\left(6\right)$式中E为误差能量，Δd^t为Δd的转置，C_D为数据协方差矩阵；当E＜δ时，输出反演结果，并停止反演，否则，继续步骤6～步骤9，直到E＜δ。δ为任意给的一个非常小的数，一般取δ＝1.0e-4～1.0e-6；步骤(5)中数据协方差矩阵用下列方法计算$C_D=\frac{1}{X}·\frac{{\sigma }_d^2}{t^{2p}}---\left(7\right)$其中X为到测量点的距离；t为地震波到达的时间；σ_d为数据方差；指数p，对于二维问题，取p≥0.5，对于三维情况，取p≥1；(6).将步骤(4)获得的残差数据作为源数据，用正演模拟类似的方法进行逆时传播，获得残差数据的剩余地震波场；①步骤(6)中的逆时传播与步骤(3)的波场正演过程相反，即解波动方程，按反时间方向进行波场传播；②与骤(3)正演模拟所用的震源相对应，步骤(6)中进行波场模拟时，震源用的是步骤(4)计算的残差数据；(7).将步骤3获得的正演地震波场和步骤(6)获得剩余地震波场分别对时间变量进行求导，进行零延迟互相关，然后对时间变量进行积分，对激发炮数进行累加，获得初始密度参数的共轭修改量；将步骤(3)获得的正演地震波场和步骤(6)获得剩余地震波场分别对空间变量求导，并进行零延迟互相关，将互相关结果对时间变量进行积分，对激发炮数进行累加，获得初始体积参数的共轭修改量，即按下式计算模型参数的共轭修改量；$\delta \hat{K}\left(x\right)=\frac{1}{K^2\left(x\right)}\underset{s}{\Sigma }{\int }_0^T\mathrm{dt}\stackrel{·}{P}(x,t;x_s)\stackrel{·}{\psi }(x,t;x_s)---\left(8\right)$$\delta \hat{\rho }\left(x\right)=\frac{1}{{\rho }^2\left(x\right)}\underset{s}{\Sigma }{\int }_0^T\mathrm{dt}\mathrm{grad}\psi (x,t;x_s)·\mathrm{grad}P(x,t;x_s)$其中P为正演波场，ψ为回传的剩余波场，或表示波场对时间的一阶偏导数，grad(·)为梯度运算，T为地震记录长度，K和ρ为模型参数，和为模型参数共轭量；(8).利用步骤(7)获得的模型参数的共轭修改量，按下列方法的迭代修改模型参数，m_n+1＝m_n-α_ng_n               (9)其中，m＝(K，ρ)为模型，m_n+1和m_n分别为第n次修改后的模型和当前模型；α_n为第n次迭代修改步长；g_n为第n次迭代的梯度；或①步骤(8)中梯度g_n用下列方法获得，$g_n=(\delta {\hat{K}}_n,\delta {\hat{\rho }}_n)---\left(10\right)$其中和是第n步迭代的模型参数的共轭修改量，用骤(7)公式计算。②步骤(8)中修改步长α_n为反演过程的控制，用下列方法计算其中，为映射函数，k为反馈增益系数，J为地震数据对模型参数的导数，E为误差能量；③步骤(8)②中映射函数用下列方法计算当n＞1时，函数当n＝1时，④步骤(8)②中反馈增益系数取满足k＞(10/t_s)的任一数，t_s为系统固有时间，t_s取5～20范围任一值；(9).将修改后的模型数据作为新的初始模型，重复步骤(3)～步骤(5)。