CN102906599B - 用于全波场反演和逆时偏移中的地下参数估计的方法 - Google Patents

Abstract

本发明涉及一种转换地震数据以获得例如体积弹性模量或密度的地下模型的方法。利用地震数据（101）和背景地下介质模型（102）计算（103）目标函数的梯度。在背景模型中计算源照明和接收器照明（104）。利用所述背景模型的速度计算地震分辨率体积（105）。通过源照明和接收器照明、地震分辨率体积和背景地下模型，梯度被转换为差异地下模型参数（106）。这些相同的因子可以用于补偿通过逆时偏移而偏移的地震数据，偏移的地震数据然后可以与地下体积弹性模量模型相关。对于迭代反演，将差异地下模型参数（106）用作预处理的梯度（107）。

Description

用于全波场反演和逆时偏移中的地下参数估计的方法

相关申请的交叉参考

本申请要求2010年2月10日提出申请、标题为“Methods for Subsurface Parameter Estimation in Full Wavefield Inversion and Reverse-Time Migration”的美国临时申请61/303,148的权益，其整个内容为全部目的以参考方式包括在本发明中。

技术领域

本发明一般涉及地球物理勘探领域，更具体地，涉及地震数据处理。具体地，本发明是一种用于在全波场反演和逆时偏移中进行地下参数估计的方法。

背景技术

勘探地震处理中的全波场反演（FWI）依赖于目标函数关于地下模型参数[12]的梯度的计算。通常以L₂范数的形式给出目标函数E：

$E=\frac{1}{2}\∫\∫\∫{|p(r_g,r_s;t)-p_b(r_g,r_s;t)|}^2\mathrm{dt}{\mathrm{dS}}_g{\mathrm{dS}}_s,---\left(1\right)$

其中p和p_b是在接收器位置r_g针对位于r_s的***的背景地下模型中的测量压力（即震波振幅）和建模压力。在迭代反演过程中，背景介质一般是先前的反演过程产生的介质。在非迭代反演过程中或偏移中，背景介质一般利用传统的地震处理技术导出，例如偏移速度分析。目标函数对整个时间t、表面S_g和S_s积分，S_g和S_s由接收器和***(shot)的遍布范围（spread）限定。我们定义K_d(r)=K(r)－K_b(r)和ρ_d(r)＝ρ(r)－ρ_b(r)，其中K（r）和ρ（r）是真实的体积弹性模量和密度，而K_b（r）和ρ_b（r）是背景模型在地下位置r的体积弹性模量和密度。我们还将测量压力和建模压力之差定义为p_d(r_g，r_s；t)－p(r_g，r_s；t)－p_b(r_g,r_s；t)。

测量压力p满足波方程 

$\mathrm{\ρ}\▿\·(\frac{1}{\mathrm{\ρ}}\▿p)-\frac{\mathrm{\ρ}}{K}\stackrel{\·\·}{p}=-q\left(t\right)\mathrm{\δ}(r-r_s)---\left(2\right)$

或 $({\mathrm{\ρ}}_b+{\mathrm{\ρ}}_d)\▿\·(\frac{1}{{\mathrm{\ρ}}_b+{\mathrm{\ρ}}_d}\▿(p_b+p_d))-\frac{{\mathrm{\ρ}}_b+{\mathrm{\ρ}}_d}{K_b+K_d}({\stackrel{\·\·}{p}}_b+{\stackrel{\·\·}{p}}_d)=-q\left(t\right)\mathrm{\δ}(r-r_s),---\left(3\right),$

其中q（t）是震源特征（source signature）。通过扩展微扰项和仅仅保留一阶波恩（Born）近似项，可以导出压力pd的波恩散射方程，

${\mathrm{\ρ}}_b\▿\·\left(\frac{1}{{\mathrm{\ρ}}_b}{\▿p}_d\right)-\frac{{\mathrm{\ρ}}_b}{K_b}{\stackrel{\·\·}{p}}_d=-[\frac{{\mathrm{\ρ}}_b{K}_d}{K_b^2}{\stackrel{\·\·}{p}}_b-{\mathrm{\ρ}}_b\▿\·\left(\frac{{\mathrm{\ρ}}_d}{{\mathrm{\ρ}}_b^2}{\▿p}_b\right)],---\left(4\right)$

由此，p_d满足

$p_d(r_g,r_s;t)=\∫[{\mathrm{\ρ}}_b\left(r^{\′}\right)\frac{K_d\left(r^{\′}\right)}{K_b^2\left(r^{\′}\right)}{\stackrel{\·\·}{p}}_b(r^{\′},r_g;t)-{\mathrm{\ρ}}_b\left(r^{\′}\right)\▿\·\left(\frac{{\mathrm{\ρ}}_d\left(r^{\′}\right)}{{\mathrm{\ρ}}_b^2\left(r^{\′}\right)}{\▿p}_b(r^{\′},r_s;t)\right)]*g_b(r_g,r^{\′};t)d{V}^{\′},---\left(5\right)$

其中V’是r’跨越的体积，而g_b是背景介质中的格林函数。

可以利用方程（5）并通过考虑由于随无穷小体积dV的微小(fractional)变化δK_b和δρ_b而产生的微小变化δp_b，可以导出p_b的梯度方程，

$\frac{\∂p_b}{{\∂K}_b\left(r\right)}=\frac{{\mathrm{\ρ}}_b\left(r\right)\mathrm{dV}}{K_b^2\left(r\right)}{g}_b(r_g,r;t)*{\stackrel{\·\·}{p}}_b(r,r_s;t),---\left(6\right)$

和 $\frac{{\∂p}_b}{{\∂\mathrm{\ρ}}_b\left(r\right)}=F^{-1}\{\frac{\mathrm{dV}}{{\mathrm{\ρ}}_b\left(r\right)}\▿G_b(r_g,r;f)\·\▿P_b(r,r_s;f)],---\left(7\right)$

其中P_b=F{p_b}，P_d=F{p_d}，G_b=F{g_b}，并且F和F^-1是傅里叶变换运算符和逆傅里叶变换运算符。

利用方程6和7并利用相互关系ρ_b(r)G_b(r_g，r)-ρ_b(r_g)G_b(r，r_g)，

$\frac{\∂E}{{\∂K}_b\left(r\right)}=-\∫\∫\∫p_d\frac{{\∂p}_b}{{\∂K}_b\left(r\right)}\mathrm{dt}{\mathrm{dS}}_g{\mathrm{dS}}_s$

$=-\frac{{\mathrm{\ρ}}_b\left(r\right)\mathrm{dV}}{K_b^2\left(r\right)}\∫\∫\∫{\left(i2\mathrm{\πf}\right)}^2{P}_b(r,r_s;f)G_b(r_g,r;f)P_d^{*}(r_g,r_s;f)\mathrm{df}{\mathrm{dS}}_g{\mathrm{dS}}_s---\left(8\right)$

$=-\frac{{\mathrm{\ρ}}_b\left(r\right)\mathrm{dV}}{K_b^2\left(r\right)}\∫\∫{\stackrel{\·}{p}}_b(r,r_s;t)\∫\frac{{\mathrm{\ρ}}_b\left(r_g\right)}{{\mathrm{\ρ}}_b\left(r\right)}{g}_b(r,r_g;-t)*{\stackrel{\·}{p}}_d(r_g,r_s;t){\mathrm{dS}}_g\mathrm{dt}{\mathrm{dS}}_s,$

和

$\frac{\∂E}{{\∂\mathrm{\ρ}}_b\left(r\right)}=-\∫\∫\∫p_d\frac{{\∂p}_b}{{\∂\mathrm{\ρ}}_b\left(r\right)}\mathrm{dt}{\mathrm{dS}}_g{\mathrm{dS}}_s$

$=-\frac{\mathrm{dV}}{{\mathrm{\ρ}}_b\left(r\right)}\∫\∫\∫P_d^{*}(r_g,r_s;f)\▿P_b(r,r_g;f)\·\▿G_b(r_g,r;f)\mathrm{df}{\mathrm{dS}}_{g}d{S}_s---\left(9\right)$

$=-\frac{\mathrm{dV}}{{\mathrm{\ρ}}_b\left(r\right)}{\∫\∫\▿p}_b(r,r_s;t)\·[\∫\▿(g_b(r,r_g;-t)*\frac{{\mathrm{\ρ}}_b\left(r_g\right)}{{\mathrm{\ρ}}_b\left(r\right)}{p}_d(r_g,r_s;t)){\mathrm{dS}}_g]\mathrm{dt}{\mathrm{dS}}_s.,$

则可以使用方程8和9以迭代方式执行全波场反演。

逆时偏移（RTM）基于与FWI中的梯度计算相似的技术，其中前向传播场与逆时接收场交叉相关。通过这种方法，RTM克服了基于光线的偏移技术如基尔霍夫偏移的限制。在RTM中，给出了在地下位置r的偏移图像场M

M(r)＝∫∫p_b(r，r_s；t)∫g_b(r，r_g；t)*p(r_g，r_s；t)dS_g dt dS_s，    (10)

这与FWI的梯度方程8很相似。

尽管方程8和9提供了用于将数据转换到地下模型的架构，但是该反演过程的收敛很缓慢。而且，由于波场的遍布范围，采用方程10的RTM在深部会遭受弱振幅。为了改善FWI的收敛或改善逆时偏移的振幅，已经通过利用目标函数的赫斯（Hessian）[9]（即目标函数的二阶导数）进行了很多尝试。然而，赫斯的计算不仅仅在计算资源方面相当昂贵，针对现实的三维反演问题也需要相当大的存储空间。而且，采用完全赫斯矩阵的FWI可能导致非最优反演[2]。

任何人可以通过将赫斯的非对角项归并为对角项来执行更稳定的反演[2]。然而，这依然需要计算完全赫斯矩阵或赫斯矩阵的至少一些非对角项，这种计算的成本是昂贵的。尽管可以选择只使用赫斯函数的对角项[11]，但是这仅仅在具有无限孔径的高频渐进区域[1,7]中是有效的。Plessix和Mulder试图通过首先计算近似对角赫斯，然后通过 对其按比例缩放来克服这些困难，其中z是深度，v_p是压缩波速[7]。根据数值实验，他们已经确定最佳的缩放参数 然而，该方法并未提供具有修正单位（correct unit）的地下介质参数的定量反演，因为只应用了近似按比例缩放。而且，该方法应用于只考虑压缩波速改变的RTM，因此可能不适用于诸如密度和剪切波速的其他弹性参数在空间变化的FWI。

发明内容

在一个实施例中，本发明是用于根据从地下区域的地震勘测获取的地震数据的反演或者根据来自地震数据的地震图像的逆时偏移确定地下区域中的物理特性的模型的方法，所述方法包含确定针对物理特性的地震分辨率体积，以及将其用作在计算机上执行的用于

转换反演中的数据失配梯度或

补偿逆时偏移的地震图像

的计算中的乘数比例缩放因子，从而获得物理特性的模型或对假设的模型的更新。

在本发明方法的一些实施例中，将数据失配（misfit）梯度或逆时偏移的地震图像乘以除地震分辨率体积之外的额外缩放比例因子，其中额外的缩放比例因子包括源照明因子、接收器照明因子和背景介质特性因子。这将产生物理特性的模型或对具有修正单位的假设模型的更新。

对于从事本领域的技术人员来说，在本发明的任何实际应用中，地震数据的反演或偏移必须在被具体编程为执行上述操作的计算机上执行是显而易见的。

附图说明

通过参考下面的详细说明和附图将更好地理解本发明及其优点，在附图中：

图1是显示本发明方法的一个实施例中的基本步骤的流程图；

图2到图5关于本发明的第一示例应用，其中图2显示利用方程8计算的目标函数关于以Pa m⁴s为单位的体积弹性模量的梯度；

图3显示利用方程18和图2中的梯度计算的以Pa为单位的体积弹性模量更新<K_d(r)>；

图4显示利用方程24和图2中的梯度计算的以Pa为单位的体积弹性模量更新<K_d(r)>；

图5显示利用方程9计算的目标函数关于以Pa²m⁷s/kg为单位的密 度的梯度；

图6和图7关于本发明的第二示例应用，其中图6显示利用方程28和图5中的梯度计算的以kg/m³为单位的密度更新<ρ_d(r)>；和

图7显示利用方程34和图5中的梯度计算的以kg/m³为单位的密度更新<ρ_d(r)>。

将结合示例实施例描述本发明。然而，在这个意义上，下面的详细说明书具体到本发明的特定实施例或特定使用，这仅仅是为了说明的目的，而并不应理解为限制本发明的范围。相反，这是为了覆盖可以包含在本发明范围内的所有替代物、改进、和等价物，如所附权利要求所限定的。

具体实施方式

在本发明中，我们利用方程8和方程9导出K_d和ρ_d的反演方程。这是通过首先利用以下事实完成的：可以利用方程5中的波恩近似扩展方程8和9中的p_d。忽视K_d和ρ_d之间的串扰（crosstalk）分量，方程8和9可以近似为：

$\begin{array}{c}\frac{\&PartialD;E}{{\&PartialD;K}_b\left(r\right)}\&ap;-\frac{{\mathrm{\&rho;}}_b\left(r\right)\mathrm{dV}}{K_b^2\left(r\right)}\&Integral;\&Integral;\&Integral;\&Integral;\frac{{\mathrm{\&rho;}}_b\left(r^{\&prime;}\right)K_d\left(r^{\&prime;}\right)}{K_b^2\left(r^{\&prime;}\right)}[g_b(r_g,r^{\&prime;};t)*{\stackrel{\&CenterDot;\&CenterDot;}{p}}_b(r^{\&prime;},r_s;t)]\\ \&times;[g_b(r_g,r;t)*{\stackrel{\&CenterDot;\&CenterDot;}{p}}_b(r,r_s;t)\mathrm{dt}{\mathrm{dS}}_g{\mathrm{dS}}_{s}d{V}^{\&prime;},\end{array}---\left(11\right)$

和

$\begin{array}{c}\frac{\&PartialD;E}{{\&PartialD;\mathrm{\&rho;}}_b\left(r\right)}\&ap;-\frac{\mathrm{dV}}{{\mathrm{\&rho;}}_b\left(r\right)}\&Integral;\&Integral;\&Integral;\&Integral;\frac{{\mathrm{\&rho;}}_d\left(r^{\&prime;}\right)}{{\mathrm{\&rho;}}_b\left(r^{\&prime;}\right)}[{\&dtri;}^{\&prime;}{P}_b^{*}(r^{\&prime;},r_s;f)\&CenterDot;{\&dtri;}^{\&prime;}{G}_b^{*}(r_g,r^{\&prime;};f)]\\ \&times;[{\&dtri;P}_b(r,r_s;f)\&CenterDot;\&dtri;G_b(r_g,r;f)]\mathrm{df}{\mathrm{dS}}_g{\mathrm{dS}}_{s}d{V}^{\&prime;}.\end{array}---\left(12\right)$

通过改变积分阶，在频域中可以将方程11重新记为

$\frac{\&PartialD;E}{{\&PartialD;K}_b\left(r\right)}\&ap;-\frac{{\mathrm{\&rho;}}_b\left(r\right)\mathrm{dV}}{K_b^2\left(r\right)}\&Integral;\&Integral;{\left(2\mathrm{\&pi;f}\right)}^4\frac{{\mathrm{\&rho;}}_b\left(r^{\&prime;}\right)\mathrm{dV}{K}_d\left(r^{\&prime;}\right)}{K_b^2\left(r^{\&prime;}\right)}$

$\&times;\{\&Integral;\frac{{\mathrm{\&rho;}}_b\left(r_g\right)}{{\mathrm{\&rho;}}_b\left(r\right)}{G}_b(r,r_g;f)G_b^{*}(r_g,r^{\&prime;};f){\mathrm{dS}}_s\}\left\{{\&Integral;P}_b(r,r_s;f)P_b^{*}(r^{\&prime;},r_s;f){\mathrm{dS}}_s\right\}{\mathrm{dV}}^{\&prime;}\mathrm{df}.$

（13）

方程13中的第一积分项

$\&Integral;\frac{{\mathrm{\&rho;}}_b\left(r_g\right)}{{\mathrm{\&rho;}}_b\left(r\right)}{G}_b(r,r_g;f)G_b^{*}(r_g,r^{\&prime;};f){\mathrm{dS}}_g---\left(14\right)$

是在r’的脉冲源产生的、在表面S_g之上测量的、之后反向传播至r的场的逆时反向传播的近似（例如，参见参考文献[8,3]）。归因于该项的波场朝着位于r’的脉冲源反向传播，并且如果积分表面S_g包含点r’，则波场的表现与t=0时的空间△（delta）函数δ(r-r’)相似。该波场与归因于第二项∫P_b(r，r_s；f)P’_b(r’，r_s；f)dS_s的波场相关，从而在r=r’附近形成梯度。然后，在r=r’附近，第一项与第二项的相关性迅速衰减。本发明认识到，相关项的振幅不可忽略的区域可通过勘测的地震分辨率来确定。在本发明中，可以做出近似

$\&Integral;{\left(2\mathrm{\&pi;f}\right)}^4\{\&Integral;\frac{{\mathrm{\&rho;}}_b\left(r_g\right)}{{\mathrm{\&rho;}}_b\left(r\right)}{G}_b(r,r_g;f)G_b^{*}(r_g,r^{\&prime;};f){\mathrm{dS}}_g\}\{\&Integral;P_b(r,r_s;f)P_b^{*}(r^{\&prime;},r_s;f){\mathrm{dS}}_s\}\mathrm{df}---\left(15\right)$

$\&ap;I_K\left(r\right)V_K\left(r\right)\mathrm{\&delta;}(r-r^{\&prime;}),$

其中

$I_K\left(r\right)={\&Integral;\left(2\mathrm{\&pi;f}\right)}^4\{\&Integral;\frac{{\mathrm{\&rho;}}_b\left(r_g\right)}{{\mathrm{\&rho;}}_b\left(r\right)}{G}_b(r,r_g;f)G_b^{*}(r_g,r;f){\mathrm{dS}}_g\}$ (16)

$\&times;\{\&Integral;P_b(r,r_s;f)P_b^{*}(r,r_s;f){\mathrm{dS}}_s\}\mathrm{df},$

V_K（r）是在地下位置r的地震分辨率。假设当非对角分量处于对角分量的地震分辨率体积内时非对角分量等于对角分量，则方程15相当于高斯-牛顿赫斯矩阵的质量集中，并且在分辨率体积以外的非对角分量等于零。换句话说，方程15相当于通过利用勘测的地震分辨率体积隐含地计算振幅有效的每个第i行中高斯-牛顿赫斯矩阵的非对角分量的数量N_i，然后将第i行的对角分量与N_i相乘。

可以将地震分辨率体积看作是r处在给定地震数据获取参数下地震成像***可以分解的最小体积。处于彼此所属的一个地震分辨率体积内的两个小目标通常不被分解，并且在地震成像***中呈现为一个目标。由于辐射方向图的差异，不同介质参数的分辨率体积是不同的。例如，由于体积弹性模量波动引起的目标产生单极辐射方向图，而由于密度波动引起的那些目标产生偶极子辐射方向图。例如利用相对廉价的射线近似法[6,4]可以计算体积弹性模量的地震分辨率体积V_K（r）。从事本领域的技术人员可能知道估计分辨率体积的其他方式。例如，一个本领域技术人员能够通过在背景介质中分配点目标和通过研究目标在地震图像中的遍布凭经验估计分辨率体积。如果背景介质包含由于反演的迭代特性 引起的速度不连续，则背景介质可能需要光滑化以进行光线追踪。技术人员也可以进行简单的假设，即波数覆盖是均匀的。在这种情况下地震分辨率体积是半径为σ≈(5/18π)^0.5v_p(r)/f_p的球体，其中f_p是峰值频率[6]。本领域技术人员也可以使用符合雷达分辨率方程[5]的近似σ≈v_p(r)T/4-v_p(r)/4B，其中T和B是源波形的有效持续时间和有效带宽。

方程11可以通过方程15简化为

$\frac{\&PartialD;E}{{\&PartialD;K}_b\left(r\right)}\&ap;-\frac{{\mathrm{\&rho;}}_b^2\left(r\right)\mathrm{dV}}{K_b^4\left(r\right)}<K_d\left(r\right)>I_K\left(r\right)V_K\left(r\right),---\left(17\right)$

由此

$<K_d\left(r\right)>\&ap;-\frac{K_b^4\left(r\right)}{{\mathrm{\&rho;}}_b^2\left(r\right)\mathrm{dV}}\frac{1}{I_K\left(r\right)V_K\left(r\right)}\frac{\&PartialD;E}{{\&PartialD;K}_b\left(r\right)},---\left(18\right)$

其中<K_d(r)>是K_d对在空间位置r的地震分辨率的空间平均值。

如果使用自由空间的格林函数（19）并假设S_g对向半个立体角，则方程16可以被进一步简化。

$G(r_g,r;f)=\frac{1}{4\mathrm{\&pi;}|r-r_g|}{e}^{\mathrm{ik}|r-r_g|},---\left(19\right)$

之后，方程16简化为

I_K(r)≈I_K，s(r)I_K，b(r)，    (20)

其中

$I_{K,s}\left(r\right)={\&Integral;\left|{\stackrel{\&CenterDot;\&CenterDot;}{p}}_b(r,r_s;t)\right|}^2\mathrm{dt},---\left(21\right)$

并且

$I_{K,g}\left(r\right)=\&Integral;\frac{{\mathrm{\&rho;}}_b\left(r_g\right)}{{\mathrm{\&rho;}}_b\left(r\right)}{G}_b(r,r_g;f)G_b^{*}(r_g,r;f){\mathrm{dS}}_g\&ap;\frac{1}{8\mathrm{\&pi;}}\frac{{\mathrm{\&rho;}}_b\left(r_g\right)}{{\mathrm{\&rho;}}_b\left(r\right)}.---\left(22\right)$

可以认为项I_K，s(r)是背景模型中的源照明，以及可以将I_K，g(r)理解为接收器照明。本领域技术人员还能够根据勘测几何改变在每个地下位置r的积分的立体角。于是，方程11变为

$\frac{\&PartialD;E}{{\&PartialD;K}_b\left(r\right)}\&ap;-\frac{{\mathrm{\&rho;}}_b^2\left(r\right)\mathrm{dV}}{K_b^4\left(r\right)}<K_d\left(r\right)>I_{K,s}\left(r\right)I_{K,g}\left(r\right)V_K\left(r\right)---\left(23\right)$

和

$<K_d\left(r\right)>\&ap;-\frac{K_b^4\left(r\right)}{{\mathrm{\&rho;}}_b^2\left(r\right)\mathrm{dV}}\frac{1}{I_{K,s}\left(r\right)I_{K,g}\left(r\right)V_K\left(r\right)}\frac{\&PartialD;E}{{\&PartialD;K}_b\left(r\right)}.---\left(24\right)$

方程18和24示出，本领域技术人员可通过按比例缩放源照明和接 收器照明、分辨率体积和背景介质特性的梯度将梯度 转换为介质参数<K_d(r)>。如果反演过程不是迭代的，那么技术人员应该能够使用方程24进行参数反演。如果反演过程是迭代的，技术人员可以使用方程24中的<K_d(r)>作为优化技术的预先处理的梯度，优化技术如最速下降法、共轭梯度（CG）法或牛顿CG法。重要的是，注意方程（18）和（24）产生具有修正单位的体积弹性模量，即，修正单位在尺寸上是修正的，因为已经考虑所有项，而并不像某些公开的方法那样为简化计算而忽略任何一项。所公开的将源照明、接收器照明、背景介质特性和地震分辨率体积中的一个或多个项忽略的方法将不会产生修正单位，因此在它们可用于迭代反演或非迭代反演之前将需要进行某种临时修复（ad hoc fix up）。

对于密度梯度，我们做出与方程15中相似的假设，

$\&Integral;\&Integral;[{\&dtri;}^{\&prime;}{P}_b^{*}(r^{\&prime;},r_s;f)\&CenterDot;{\&dtri;}^{\&prime;}{G}_b^{*}(r_g,r^{\&prime;};f)]$

(25)

$\&times;[\&dtri;P_b(r,r_s;f)\&CenterDot;\&dtri;G_b(r_g,r;f)]{\mathrm{dS}}_g{\mathrm{dS}}_s\mathrm{df}\&ap;I_{\mathrm{\&rho;}}\left(r\right)V_{\mathrm{\&rho;}}\left(r\right)\mathrm{\&delta;}(r-r^{\&prime;}),$

其中

$I_{\mathrm{\&rho;}}\left(r\right)={\&Integral;\&Integral;\&Integral;|\&dtri;P_b(r,r_s;f)\&CenterDot;{\&dtri;G}_b(r_g,r;f)|}^2{\mathrm{dS}}_g{\mathrm{dS}}_s\mathrm{df},---\left(26\right)$

并且V_ρ（r）是在地下位置r的密度ρ的地震分辨率。分辨率体积V _ρ（r）与V_k（r）不同，因为当入射场和散射场几乎彼此垂直时波数会消失。这归因于先前讨论的密度波动的偶极子辐射方向图。像对于V_K（r）一样，技术人员可以使用光线追踪计算分辨率体积V_ρ（r），同时考虑这些消失的接近垂直的波数。可替代地，技术人员能够通过忽略波数覆盖的差异来假设v_K(r)≈v_p(r)。

然后，可以将梯度方程12重写为

$\frac{\&PartialD;E}{{\&PartialD;\mathrm{\&rho;}}_b\left(r\right)}\&ap;-\frac{\mathrm{dV}}{{\mathrm{\&rho;}}_b^2\left(r\right)}<{\mathrm{\&rho;}}_d\left(r\right)>I_{\mathrm{\&rho;}}\left(r\right)V_{\mathrm{\&rho;}}\left(r\right),---\left(27\right)$

因此

$<{\mathrm{\&rho;}}_d\left(r\right)>\&ap;-\frac{{\mathrm{\&rho;}}_b^2\left(r\right)}{\mathrm{dV}}\frac{1}{I_{\mathrm{\&rho;}}\left(r\right)V_{\mathrm{\&rho;}}\left(r\right)}\frac{\&PartialD;E}{{\&PartialD;\mathrm{\&rho;}}_b\left(r\right)},---\left(28\right)$

其中<ρ_d(r)>是ρ_d（r）在地震分辨率体积V_ρ（r）上的空间平均。

我们可以通过利用矢量恒等式(a·b)(c·d)=(a·d)(b·c)+(a×c)·(b×d)进一步 简化方程25从而获得

$[{\&dtri;}^{\&prime;}{P}_b^{*}(r^{\&prime;},r_s;f)\&CenterDot;{\&dtri;}^{\&prime;}{G}_b^{*}(r_g,r^{\&prime;};f)][\&dtri;P_b(r,r_s;f)-\&dtri;G_b(r_g,r;f)]$

$\text{=}[{\&dtri;}^{\&prime;}{P}_b^{*}(r^{\&prime;},r_s;f)\&CenterDot;\&dtri;P_b(r,r_s;f)][{\&dtri;}^{\&prime;}{G}_b^{*}(r_g,r^{\&prime;};f)\&dtri;\&dtri;G_b(r_g,r;f)]---\left(29\right)$

$+[{\&dtri;}^{\&prime;}{P}_b^{*}(r^{\&prime;},r_s;f)\&times;\&dtri;G_b(r_g,r^{};f)]\&CenterDot;[{\&dtri;}^{\&prime;}{G}_b^{*}(r_g,r^{\&prime;};f)\&times;\&dtri;P_b(r,r_s;f)].$

方程29右侧的第二项是散射场的偶极子辐射方向图的修正项，因此当 和 彼此垂直时方程29右侧的第二项达到最大值。忽略该修正项，

$[{\&dtri;}^{\&prime;}{P}_b^{*}(r^{\&prime;},r_s;f)\&CenterDot;{\&dtri;}^{\&prime;}{G}_b^{*}(r_g,r^{\&prime;};f)][\&dtri;P_b(r,r_s;f)\&CenterDot;\&dtri;G_b(r_g,r;f)]$

(30)

$\&ap;[{\&dtri;}^{\&prime;}{P}_b^{*}(r^{\&prime;},r_s;f)\&CenterDot;\&dtri;P_b(r,r_s;f)][{\&dtri;}^{\&prime;}{G}_b^{*}(r_g,r^{\&prime;};f)\&CenterDot;\&dtri;G_b(r_g,r;f)].$

于是，方程12可以近似为

$\frac{\&PartialD;E}{{\&PartialD;\mathrm{\&rho;}}_b\left(r\right)}\&ap;-\frac{\mathrm{dV}}{{\mathrm{\&rho;}}_b\left(r\right)}\&Integral;\&Integral;\frac{{\mathrm{\&rho;}}_d\left(r^{\&prime;}\right)}{{\mathrm{\&rho;}}_b\left(r^{\&prime;}\right)}\{\&Integral;[{\&dtri;}^{\&prime;}{G}_b^{*}(r_g,r^{\&prime;};f)\&CenterDot;\&dtri;G_b(r_g,r;f)\left]{\mathrm{dS}}_g\right\}$

(31)

$\&times;\{\&Integral;[{\&dtri;}^{\&prime;}{P}_b^{*}(r^{\&prime;},r_s;f)\&CenterDot;\&dtri;P_b(r,r_s;f)\left]{\mathrm{dS}}_s\right\}\mathrm{df}{\mathrm{dV}}^{\&prime;}$

我们可以通过利用自由空间的格林函数来近似方程26中在S_g上的积分，

${\&Integral;\&dtri;G}_b^{*}(r_g,r;f)\&CenterDot;\&dtri;G_b(r_g,r;f){\mathrm{dS}}_g\&ap;\frac{{\mathrm{\&rho;}}_b\left(r_g\right)}{{\mathrm{\&rho;}}_b\left(r\right)}\&Integral;\&dtri;G_b^{*}(r_g,r;f)\&CenterDot;\&dtri;G_b(r,r_g;f){\mathrm{dS}}_g\&ap;\frac{{\mathrm{\&rho;}}_b\left(r_g\right)}{{\mathrm{\&rho;}}_b\left(r\right)}\frac{{\left(2\mathrm{\&pi;f}\right)}^2}{8\mathrm{\&pi;}{v}_p^2\left(r\right)}$

（32）

在假设 和ρ_b(r_g)是沿着S_g的常量的情况下，在半个立体角上进行积分。

之后，以上梯度方程可以改写为

$\frac{\&PartialD;E}{{\&PartialD;\mathrm{\&rho;}}_b\left(r\right)}\&ap;-\frac{\mathrm{dV}}{{\mathrm{\&rho;}}_b^2\left(r\right)}<{\mathrm{\&rho;}}_d\left(r\right)>I_{\mathrm{\&rho;},s}\left(r\right)I_{\mathrm{\&rho;},g}\left(r\right)V_{\mathrm{\&rho;}}\left(r\right)---\left(33\right)$

因此

$<{\mathrm{\&rho;}}_d\left(r\right)>\&ap;-\frac{{\mathrm{\&rho;}}_b^2\left(r\right)}{\mathrm{dV}}\frac{1}{I_{\mathrm{\&rho;},s}\left(r\right)I_{\mathrm{\&rho;},g}\left(r\right)V_{\mathrm{\&rho;}}\left(r\right)}\frac{\&PartialD;E}{{\&PartialD;\mathrm{\&rho;}}_b\left(r\right)},---\left(34\right)$

其中

$I_{\mathrm{\&rho;},s}\left(r\right)={|\&dtri;\stackrel{\&CenterDot;}{p}(r,r_s;t)|}^2\mathrm{dt}{\mathrm{dS}}_s,---\left(35\right)$

和

$I_{\mathrm{\&rho;},g}\left(r\right)=\frac{{\mathrm{\&rho;}}_b\left(r_g\right)}{{\mathrm{\&rho;}}_b\left(r\right)}\frac{1}{8\mathrm{\&pi;}{v}_p^2\left(r\right)}.---\left(36\right)$

类似<K_d(r)>的情形，方程28或34可以用作非迭代反演的反演公式，或用作迭代反演的预处理的梯度方程。重要的是，要注意到这些方程产生 具有修正单位的密度，即，修正单位在尺寸上是修正的，因为已经考虑所有项，而没有像某些公开的方案那样为了简化计算而忽视任何一项。这对于体积弹性模量的方程18和24也成立。所公开的将源照明、接收器照明、背景介质特性和地震分辨率体积中一个或多个项忽略的方案不会产生修正单位，因此，在用于迭代反演或非迭代反演之前需要进行某种临时修复。

由于FWI的第一迭代类似于RTM，因此可应用此处提供的方法以较少的改动分析RTM中的振幅项。包括RTM的地震偏移通常用于成像地下的结构，因此，通常丢弃偏移图像中的振幅信息。我们示出当利用此处提供的方法适当地按比例缩放时，RTM振幅表示地下的真实压缩波速和背景模型速度之差。

我们注意到，RTM方程10缺少方程8中入射场的双重微分。该双重微分表示在经典的瑞利（Rayleigh）散射区域[Refs.[10,13]]中高频分量比低频分量更加有效地散射。因此，我们可以将方程10视为方程8中的梯度计算操作，部分地忽略散射场的频率依赖，

$\frac{\&PartialD;E}{{\&PartialD;K}_b\left(r\right)}\&ap;-\frac{{\mathrm{\&rho;}}_b\mathrm{dV}{\left(i2\mathrm{\&pi;}{f}_c\right)}^2}{K_b^2\left(r\right)}\&Integral;\&Integral;\&Integral;P_b(r,r_s;f)G_b(r,r_g;f)P^{*}(r_g,r_s;f)\mathrm{df}{\mathrm{dS}}_g{\mathrm{dS}}_s$

(37)

$\&ap;-\frac{{\mathrm{\&rho;}}_b\mathrm{dV}{\left(i2\mathrm{\&pi;}{f}_c\right)}^2}{K_b^2\left(r\right)}M\left(r\right),$

其中f_c是源波形的中心频率。频率依赖被部分忽略，因为虽然已经忽略前向场p_b中的频率依赖，但是无法忽略接收场p_s中暗含的频率依赖。在RTM中一般不考虑密度的空间变化，所以在方程37中假设ρ_b是常量。

现在可以将用于导出方程17和23的相同近似应用于方程10，

$\frac{\&PartialD;E}{{\&PartialD;K}_b\left(r\right)}\&ap;-\frac{{\mathrm{\&rho;}}_b^2{I}_{K,g}\left(r\right)V_K\left(r\right)\mathrm{dV}<K_d\left(r\right)>{\left(i2\mathrm{\&pi;}{f}_c\right)}^2}{K_b^4\left(r\right)}{\&Integral;\&Integral;\left(i2\mathrm{\&pi;f}\right)}^2{P}_b(r,r_s;f)P_b^{*}(r,r_s;f)\mathrm{df}{\mathrm{dS}}_s$

（38）

该方程连同方程37产生

$<K_d\left(r\right)>\&ap;\frac{K_b^2\left(r\right)}{{\mathrm{\&rho;}}_b{I}_{K,g}\left(r\right)V_K\left(r\right)}\frac{M\left(r\right)}{\&Integral;\&Integral;{\left|{\stackrel{\&CenterDot;}{p}}_b(r,r_s;t)\right|}^2\mathrm{dt}{\mathrm{dS}}_s}.---\left(39\right)$

方程39能够定量分析逆时偏移的图像中的振幅。更具体地，它能够将振幅反演为地下的差异体积弹性模量。

图1是示出了本发明的方法的一个实施例的基本步骤的流程图。

在步骤103中，利用输入的地震记录（101）和关于背景地下介质（102）的信息计算目标函数的梯度。在步骤104中，计算背景模型中的源照明和接收器照明。在步骤105中，利用背景模型的速度计算地震分辨率体积。在步骤106中，利用步骤104的源照明和接收器照明、步骤105的地震分辨率体积和背景地下模型（102）将步骤103的梯度转换为差异地下模型参数。如果要执行迭代反演过程，则在步骤107中将步骤106中的差异地下模型参数用作迭代反演过程的预处理的梯度。

示例

我们考虑以下情形：K_b=9MPa和ρ_b=1000kg/m³的均质介质中的30m×30m×30m“完美”波恩散射体。目标集中在（x,y,z）=（0,0,250m），其中x和y是两个水平坐标，z是深度。在图2到图7中可以看到目标是位于每个图形中心的小的方形3×3阵列。我们假设源和接收器共同位于x和y方向均具有10m间距的-500m≤x≤500m和-500m≤y≤500m区间内。我们假设源小波在1Hz到51Hz频带内在1m处具有1Pz/Hz的均匀振幅。

在第一示例中，我们假设目标具有给定的体积弹性模量波动K_d=900kPa。图2显示采用方程8时沿着y=0平面的梯度 已经利用方程5计算方程8中的散射场p_d。图2中的梯度具有Pa m⁴s的单位，因此不能与K_d直接相关。如在以上的“背景技术”部分中所述的，这是在一些公开的根据目标函数的梯度计算模型更新的尝试中遇到的难题。

图3显示采用本发明的方程18时的<K_d(r)>。已经假设，地震分辨率体积V_K（r）是半径为σ＝v_p(r)／4B=15m的球体。可以看到，图3是目标的模糊图像，因为<K_d(r)>是地震分辨率体积上的平均特性。图4是使用不太严谨的方程24时的<K_d(r)>。可以看到，图3和图4中的<K_d(r)>彼此良好地吻合。在图3中的目标中心的<K_d(r)>值是752kPa，在图4中是735kPa，两者均在900kPa的真实值的20%内。

第二示例是目标具有ρ_d=100kg/m³的密度波动的情况。图5显示采用方程9时沿着y=0平面的梯度 已经利用图5计算方程9中 的散射场。图5中的梯度具有Pa²m⁷s/kg的单位。如同第一实例，不同的单位阻止梯度与密度更新直接相关，再次图示了公开的方法中遇到的难题。

图6示出采用方程28时的<ρ_d(r)>。已经将地震分辨率体积V_ρ（r）假设为与V_k（r）相同。图7是使用方程34时的<ρ_d(r)>。采用方程34得到的<ρ_d(r)>的估计将导致反演的准确性比利用方程28更低，因为方程31中忽略了偶极子照明项。

前述的专利申请针对本发明的特定实施例，目的在于对本发明进行图解说明。然而，对于本领域的技术人员来说，可以对本文中描述的实施例做出多种修改和变更是显而易见的。所有的修改和变更均在本发明的范围内，如所附权利要求所限定的。本领域技术人员将容易地认识到，在本发明的实际应用中，可以在计算机上或借助计算机执行本发明方法的至少某些步骤，即，本发明是计算机可实施的。

参考书目

[1]G.Beylkin,"Imaging of discontinuities in the inverse scattring problem by inversion of a causal generalized Radon transform,"J.Math.Phys.26,99-108(1985)。

[2]G.Chavent和R.-E.Plessix,"An optimal true-amplitude least-squares prestack depth-migration operator,"Geophysics 64(2),508-515(1999)。

[3]D.R.Jackson和D.R.Dowling,"Phase conjugation in underwater acoustics,"J.Acoust.Soc.Am.89(1),171-181(1991)。

[4]I.Lecomte,"Resolution and illumination analyses in PSDM:A ray-based approach,"The Leading Edge,pages 650—663(May 2008)。

[5]N.Levanon,Radar Principles,chapter 1,pages 1-18,John Wiley&Sons,New York 25(1988)。

[6]M.A.Meier和P.J.Lee,"Converted wave resolution,"Geophysics74(2),Q1-Q16(2009)。

[7]R.E.Plessix和W.A.Mulder,"Frequency-domain finite-difference amplitude-preserving migration,"Geophys.J.Int.157,975-987(2004)。

[8]R.P.Porter,"Generalized holography with application to inverse  scattering and inverse source problems,"Progress in Optics XXVII,E.Wolf，editor,pages 317-397,Elsevier(1989)。

[9]R.G.Pratt,C.Shin,和G.J.Hicks,"Gauss-Newton and full Newton methods in frequency-space seismic waveform inversion,"Geophys.J.Int.133,341-362(1998)。

[10]J.W.S.Rayleigh,"On the transmission of light through an atmosphere containing small particles in suspension,and on the origin of the blue of the sky,"Phil.Mag.47,375-384(1899)。

[11]C.Shin,S.Jang,和D.-J.Min,"Waveform inversion using a logarithmic wavefield,"Geophysics 49,592-606(2001)。

[12]A.Tarantola,"Inversion of seismic reflection data in the acoustic approximation,"Geophysics 49，1259-1266(1984)。

[13]R.J.Urick,Principles of Underwater Sound，chapter 9,pages291-327,McGraw-Hill,New York,3rd edition(1983)。

Claims (21)

1.一种根据从地下区域的地震勘测获得的地震数据的反演或者根据来自所述地震数据的地震图像的逆时偏移确定地下区域的物理特性的模型的方法，所述方法包括：

确定针对所述物理特性的地震分辨率体积，并将其用作在计算机上执行的用于

转换反演中的数据失配梯度或

补偿逆时偏移的地震图像

的计算中的乘数比例缩放因子，从而获得所述物理特性的模型或假设的模型的更新。

2.根据权利要求1所述的方法，进一步包括用额外的比例缩放因子与数据失配梯度或逆时偏移的地震图像相乘以获得所述物理特性的模型或假设的模型的更新，其中所述额外的比例缩放因子包括源照明因子、接收器照明因子和背景介质特性因子。

3.根据权利要求1所述的方法，其中所述地震分辨率体积是通过光线追踪利用背景介质模型的速度和利用假设的地震小波的频率的函数确定的。

4.根据权利要求2所述的方法，其中所述模型是根据地震数据的反演确定的，所述方法进一步包括：

假设地下区域的初始模型，其指定所述地下区域中的离散单元位置处的模型参数；

形成数学的目标函数以测量所测量的地震数据和根据模型计算的地震数据之间的失配；

选择提供对初始模型的调整即更新的数学关系，其中所述调整会降低失配，所述数学关系使所述调整与所述目标函数的按比例缩放的梯度相关，所述梯度是关于所述模型参数的，所述比例缩放包含四个比例缩放因子，即

地震分辨率体积因子，

源照明因子，

接收器照明因子，和

背景介质特性因子，

所有这些因子在所述数学关系中表现为乘数因子，其对所述目标函数的梯度进行比例缩放，从而产生所述模型参数的调整；以及

利用计算机根据所述数学关系计算所述调整，然后，利用计算的调整更新所述初始模型。

5.根据权利要求4所述的方法，其中所述物理特性即所述模型参数是体积弹性模量或密度，或体积弹性模量和密度的组合。

6.根据权利要求4所述的方法，其中所述数学关系依赖于所述物理特性。

7.根据权利要求4所述的方法，其中当所述物理特性是体积弹性模量时，所述背景介质特性因子包括除以密度平方的体积弹性模量的四次方，以及当所述物理特性是密度时，所述背景介质特性因子包括密度的平方。

8.根据权利要求4所述的方法，其中当所述物理特性是体积弹性模量时，所述接收器照明因子近似为(1/8π)(ρ_b(r_g)/ρ_b(r))，其中ρ_b(r)是在位置r的背景密度，r_g是接收器的位置；以及当所述物理特性是密度时，接收器照明因子I_ρ,g(r)近似为：

$I_{\mathrm{\&rho;},g}\left(r\right)=\frac{{\mathrm{\&rho;}}_b\left(r_g\right)}{{\mathrm{\&rho;}}_b\left(r\right)}\frac{1}{8\mathrm{\&pi;}{v}_p^2\left(r\right)},$

其中v_p(r)是在位置r的速度。

9.根据权利要求4所述的方法，进一步包括重复所述方法进行至少一次迭代，其中来自先前迭代的更新的模型取代所述初始模型。

10.根据权利要求9所述的方法，其中所述目标函数的函数形式和所述数学关系的函数形式从一次迭代到下一次迭代不改变。

11.根据权利要求4所述的方法，其中对所述初始模型的调整是通过利用目标函数的赫斯最小化所述目标函数来计算的，其中所述目标函数的赫斯可产生赫斯矩阵，其中当处于地震分辨率体积外部时，所述赫斯矩阵的非对角元素被忽略。

12.根据权利要求11所述的方法，其中假设地震分辨率体积内的所述赫斯矩阵的非对角元素等于对应的对角元素，导致只计算对角元素。

13.根据权利要求2所述的方法，其中当所述物理特性是体积弹性模量时，通过自由空间中格林函数在所述地震勘测的接收器遍布范围所限定的表面上的积分近似所述接收器照明因子。

14.根据权利要求4所述的方法，其中当所述物理特性是密度时，通过自由空间中格林函数的梯度在所述地震勘测的接收器遍布范围所限定的表面上的积分近似所述接收器照明因子。

15.根据权利要求1所述的方法，其中基于均匀的波数覆盖的假设，将所述地震分辨率体积近似为球体。

16.根据权利要求4所述的方法，进一步包含在迭代优化技术中利用所述更新的模型预处理所述梯度。

17.根据权利要求2所述的方法，其中所述模型是根据来自所述地震数据的地震图像的逆时偏移确定的。

18.根据权利要求17所述的方法，其中所述物理特性是体积弹性模量。

19.根据权利要求17所述的方法，其中所述背景介质特性因子包含除以密度的体积弹性模量平方。

20.根据权利要求17所述的方法，其中通过(1/8π)(ρ_b(r_g)/ρ_b(r))近似所述接收器照明因子，其中ρ_b(r)是在位置r的背景密度，r_g是接收器的位置。

21.根据权利要求17所述的方法，其中偏移图像的地震振幅被转换为差异体积弹性模量或差异压缩波速。