CN103308941A - 一种基于任意广角波动方程的成像方法及装置 - Google Patents

Abstract

本发明实施例提供一种基于任意广角波动方程的成像方法及装置，所述方法包括：获取建立的速度模型和给定的震源波场；利用正演模拟得到速度模型的合成地震记录；对正演模拟的合成地震记录和给定的震源波场做倾角滤波处理；从频散关系出发，分析任意广角波动方程的成像精度，得到任意广角波动方程的最优化参考速度；采用有限差分方法求解最优化参考速度的任意广角波动方程；采用互相关成像条件进行偏移成像，以提高任意广角波动方程的偏移成像精度。可以有效地压制任意广角波动方程叠前深度偏移的倏逝波成像噪音，提高了成像质量；所提出的参考速度优化方法，提高了任意广角波动方程叠前深度偏移成像精度，因此可以对更陡倾角的地层成像。

Description

一种基于任意广角波动方程的成像方法及装置

技术领域

本发明涉及勘探地球物理领域，具体地说涉及一种基于任意广角波动方程的成像方法及装置。

背景技术

勘探目标的复杂化和勘探精度的提高，要求地震偏移成像方法可以对高陡倾角构造和横向速度变化较剧烈的地质体准确成像，这在很大程度上促进了叠前深度偏移方法的发展。

目前，叠前深度偏移方法主要有基于射线理论的Kirchhoff积分类偏移成像方法和基于波场延拓的波动方程偏移成像方法。积分类方法的几何意义比较直观，具有灵活的目标处理能力和较高的计算效率，目前已成为工业界普遍采用的叠前深度偏移成像方法。但是高频近似假设使其在地下介质比较复杂的情况下得不到满意的成像结果。基于波场延拓的波动方程偏移方法可以很好地保留地震波的运动学和动力学特征，在对复杂构造区域成像时，具有一定优势。

基于波场延拓的波动方程偏移方法主要包含单程波方程和双程波方程两类：基于双程波方程的逆时偏移方法（RTM）完全遵循全波方程，不存在偏移倾角限制，可以处理速度场的剧烈变化，在保幅处理方面也存在一定优势。但是其在计算量和数据读写方面要求很高，因此至今仍未得到工业界大规模广泛应用。与逆时偏移相比，单程波方法的计算效率较高，但是因为存在对平方根算子的近似，所以对成像倾角和速度的横向变化剧烈程度有一定限制。

针对传统偏移方法中存在的一些问题，Guddati（2005）提出基于任意广角波动方程（AWWE）的偏移方法。与逆时偏移方法相比，该偏移方法可以避免低频成像噪音和大存储量问题；与传统的单程波偏移方法相比，该偏移方法更适用于速度横向变化较大的复杂构造区以及陡倾角界面的成像。但是因为AWWE存在对平方根算子的近似，所以它的成像精度有待提高。

发明内容

本发明实施例提供一种基于任意广角波动方程的成像方法及装置，所提出的倾角滤波方法和参考速度优化方法可以提高任意广角波动方程的偏移成像精度。

一方面，本发明实施例提供了一种基于任意广角波动方程的成像方法，所述基于任意广角波动方程的成像的方法，包括：

获取建立的速度模型和给定的震源波场；

利用正演模拟得到速度模型的合成地震记录；

对正演模拟的合成地震记录和给定的震源波场做倾角滤波处理；

从频散关系出发，分析任意广角波动方程的成像精度，得到任意广角波动方程的最优化参考速度；

采用有限差分方法求解最优化参考速度的任意广角波动方程；

采用互相关成像条件进行偏移成像。

可选的，在本发明一实施例中，所述利用正演模拟得到速度模型的合成地震记录，包括：采用交错网格有限差分方法求解时空域各向同性均匀介质的一阶速度—应力声波波动方程，从而得到所述速度模型的合成地震记录。

可选的，在本发明一实施例中，所述对正演模拟的合成地震记录和给定的震源波场做倾角滤波处理，包括：利用频率波数域的扇形滤波器对正演模拟的合成地震记录和给定的震源波场做倾角滤波处理。

可选的，在本发明一实施例中，所述从频散关系出发，分析任意广角波动方程的成像精度，得到任意广角波动方程的最优化参考速度，包括：从频散关系出发，分析任意广角波动方程的成像精度，利用非线性牛顿法得到任意广角波动方程的最优化参考速度。

可选的，在本发明一实施例中，所述采用有限差分方法求解最优化参考速度的任意广角波动方程，包括：采用有限差分方法求解最优化参考速度的任意广角波动方程，并在边界处采用完美匹配层PML吸收边界条件。

另一方面，本发明实施例提供了一种基于任意广角波动方程的成像装置，所述基于任意广角波动方程的成像的装置，包括：

获取单元，用于获取建立的速度模型和给定的震源波场；

正演模拟单元，用于利用正演模拟得到速度模型的合成地震记录；

倾角滤波单元，用于对正演模拟的合成地震记录和给定的震源波场做倾角滤波处理；

分析单元，用于从频散关系出发，分析任意广角波动方程的成像精度，得到任意广角波动方程的最优化参考速度；

方程求解单元，用于采用有限差分方法求解最优化参考速度的任意广角波动方程；

成像单元，用于采用互相关成像条件进行偏移成像。

可选的，在本发明一实施例中，所述正演模拟单元，进一步用于采用交错网格有限差分方法求解时空域各向同性均匀介质的一阶速度—应力声波波动方程，从而得到所述速度模型的合成地震记录。

可选的，在本发明一实施例中，所述倾角滤波单元，进一步用于利用频率波数域扇形滤波器对正演模拟的合成地震记录和给定的震源波场做倾角滤波处理。

可选的，在本发明一实施例中，所述分析单元，进一步用于从频散关系出发，分析任意广角波动方程的成像精度，利用非线性牛顿法得到任意广角波动方程的最优化参考速度。

可选的，在本发明一实施例中，所述方程求解单元，进一步用于采用有限差分方法求解最优化参考速度的任意广角波动方程，并在边界处采用完美匹配层PML吸收边界条件。

上述技术方案具有如下有益效果：因为采用所述基于任意广角波动方程的成像的方法，包括：获取建立的速度模型和给定的震源波场；利用正演模拟得到速度模型的合成地震记录；对正演模拟的合成地震记录和给定的震源波场做倾角滤波处理；从频散关系出发，分析任意广角波动方程的成像精度，得到任意广角波动方程的最优化参考速度；采用有限差分方法求解最优化参考速度的任意广角波动方程；采用互相关成像条件进行偏移成像的技术手段，所以达到了如下的技术效果：(1)提出的倾角滤波方法可以有效地压制任意广角波动方程叠前深度偏移的倏逝波成像噪音，提高成像质量；(2)优化了任意广角波动方程的参考速度参数，提高了任意广角波动方程叠前深度偏移成像精度，因此使基于任意广角波动方程的叠前深度偏移方法可以对更陡倾角的地层成像。

附图说明

为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例或现有技术描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他的附图。

图1为本发明实施例提供了一种基于任意广角波动方程的成像方法流程图；

图2为本发明实施例一种基于任意广角波动方程的成像装置结构示意图；

图3(a)为本发明应用实例AWWE的频散关系曲线实部之间的关系示意图；

图3(b)为本发明应用实例AWWE的频散关系曲线虚部之间的关系示意图；

图4为本发明应用实例频率波数域扇形滤波器的频率波数响应示意图；

图5(a)为本发明应用实例优化二参数的任意广角波动方程的频散关系曲线示意图；

图5(b)为本发明应用实例优化三参数的任意广角波动方程的频散关系曲线示意图；

图6为本发明应用实例复杂介质速度模型示意图；

图7为本发明应用实例复杂介质速度模型的多炮偏移成像结果示意图。

具体实施方式

下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。

如图1所示，为本发明实施例提供了一种基于任意广角波动方程的成像方法流程图，所述基于任意广角波动方程的成像的方法，包括：

101、获取建立的速度模型和给定的震源波场；

102、利用正演模拟得到速度模型的合成地震记录；

103、对正演模拟的合成地震记录和给定的震源波场做倾角滤波处理；

104、从频散关系出发，分析任意广角波动方程的成像精度，得到任意广角波动方程的最优化参考速度；

105、采用有限差分方法求解最优化参考速度的任意广角波动方程；

106、采用互相关成像条件进行偏移成像。

可选的，所述利用正演模拟得到速度模型的合成地震记录，包括：采用交错网格有限差分方法求解时空域各向同性均匀介质的一阶速度—应力声波波动方程，从而得到所述速度模型的合成地震记录。

可选的，所述对正演模拟的合成地震记录和给定的震源波场做倾角滤波处理，包括：利用频率波数域的扇形滤波器对正演模拟的合成地震记录和给定的震源波场做倾角滤波处理。

可选的，所述从频散关系出发，分析任意广角波动方程的成像精度，得到任意广角波动方程的最优化参考速度，包括：从频散关系出发，分析任意广角波动方程的成像精度，利用非线性牛顿法得到任意广角波动方程的最优化参考速度。

可选的，所述采用有限差分方法求解最优化参考速度的任意广角波动方程，包括：采用有限差分方法求解最优化参考速度的任意广角波动方程，并在边界处采用完美匹配层PML（Perfectly Matched Layer）吸收边界条件。

如图2所示，为本发明实施例一种基于任意广角波动方程的成像装置结构示意图，所述基于任意广角波动方程的成像的装置，包括：

获取单元21，用于获取建立的速度模型和给定的震源波场；

正演模拟单元22，用于利用正演模拟得到速度模型的合成地震记录；

倾角滤波单元23，用于对正演模拟的合成地震记录和给定的震源波场做倾角滤波处理；

分析单元24，用于从频散关系出发，分析任意广角波动方程的成像精度，得到任意广角波动方程的最优化参考速度；

方程求解单元25，用于采用有限差分方法求解最优化参考速度的任意广角波动方程；

成像单元26，用于采用互相关成像条件进行偏移成像。

可选的，所述正演模拟单元22，进一步用于采用交错网格有限差分方法求解时空域各向同性均匀介质的一阶速度—应力声波波动方程，从而得到所述速度模型的合成地震记录。

可选的，所述倾角滤波单元23，进一步用于利用频率波数域的扇形滤波器对正演模拟的合成地震记录和给定的震源波场做倾角滤波处理。

可选的，所述分析单元24，进一步用于从频散关系出发，分析任意广角波动方程的成像精度，利用非线性牛顿法得到任意广角波动方程的最优化参考速度。

可选的，所述方程求解单元25，进一步用于采用有限差分方法求解最优化参考速度的任意广角波动方程，并在边界处采用完美匹配层PML吸收边界条件。

本发明实施例上述技术方案具有如下有益效果：因为采用所述基于任意广角波动方程的成像的方法，包括：获取建立的速度模型和给定的震源波场；利用正演模拟得到速度模型的合成地震记录；对正演模拟的合成地震记录和给定的震源波场做倾角滤波处理；从频散关系出发，分析任意广角波动方程的成像精度，得到任意广角波动方程的最优化参考速度；采用有限差分方法求解最优化参考速度的任意广角波动方程；采用互相关成像条件进行偏移成像的技术手段，所以达到了如下的技术效果：(1)提出的倾角滤波方法可以有效地压制任意广角波动方程叠前深度偏移的倏逝波成像噪音，提高成像质量；(2)优化了任意广角波动方程的参考速度参数，提高了任意广角波动方程叠前深度偏移成像精度，因此使基于任意广角波动方程的叠前深度偏移方法可以对更陡倾角的地层成像。

以下举应用实例进行详细说明：

本发明应用实例的目的在于提供一种提高任意广角波动方程成像精度的方法。提出了一种基于倾角滤波的成像噪音压制方法，提高了偏移成像质量。从频散关系出发，用非线性牛顿法求解单程波算子的最优化问题，使任意广角波动方程的频散关系与准确的平方根算子达到最佳匹配，提高了方程的成像倾角，使之适合更陡倾角界面的成像。

第一步，建立速度模型，给定震源波场。

根据设计要求建立一个理论速度模型，设定好地下介质速度。根据震源子波给定震源波场。

第二步，正演模拟得到速度模型的合成地震记录。

采用交错网格有限差分方法求解时空域各向同性均匀介质的一阶速度—应力声波波动方程，从而得到所设计速度模型的合成地震记录。均匀各项同性介质一阶速度—应力声波方程的形式为

$\frac{\∂p}{\∂t}=-K(\frac{\∂v_x}{\∂x}+\frac{\∂v_z}{\∂z})$

$\frac{\∂v_x}{\∂t}=-\frac{1}{\mathrm{\ρ}}\frac{\∂p}{\∂x}---\left(1\right)$

$\frac{\∂v_z}{\∂t}=-\frac{1}{\mathrm{\ρ}}\frac{\∂p}{\∂z}$

其中K，ρ为介质的体积模量和密度，它们与声波传播速度V之间的关系为K＝ρV²，P,v_x,v_z为介质中压力和质点振动速度在两个坐标方向的分量，x,z分别为水平和垂直坐标，t为时间变量。

第三步，对正演模拟的合成地震记录和给定的震源波场做倾角滤波处理。

如图3(a)所示，为本发明应用实例AWWE的频散关系曲线实部之间的关系示意图，如图3(b)所示，为本发明应用实例AWWE的频散关系曲线虚部之间的关系示意图，其中exact表示准确的频散关系，AWWE表示任意广角波动方程的频散关系。通过分析图3(a)、图3(b)的频散关系曲线和一系列的推导论证，得出偏移成像噪音是由于不准确的倏逝波近似产生的。我们根据倏逝波的慢度特征，在频率—波数域设计滤波器，将频率波数域的扇形滤波器作用于接收到的记录波场和给定的震源波场，这样就可以滤除不准确倏逝波。频率波数域的扇形滤波器的表达式为

$H(\mathrm{\&omega;},k)=\left\{\begin{array}{cc}1,& \left|k\right|\&le;k_1\\ g\left(k\right),& k_1<\left|k\right|<k_2\\ 0,& \left|k\right|\&GreaterEqual;k_2\end{array}\right.---\left(2\right)$

其中，k₁＝ω/c-10dk,k₂＝ω/c+10dk，H(ω,k)为频率波数域的扇形滤波器的频率—波数响应函数，k为圆波数，c为检波器所在位置处的地层速度，g(k)为余弦镶边函数，k₁,k₂为镶边函数的左右边界，镶边函数的长度取为二十倍采样波数间隔dk。如图4所示，为本发明应用实例频率波数域的扇形滤波器的频率波数响应示意图。为了避免吉普斯效应，我们对频率波数域的扇形滤波器进行镶边，采用的余弦镶边函数为

$g(\mathrm{\&omega;},k)=\frac{1}{2}(1+\cos \frac{\left|k\right|-k_1}{k_2-k_1}\mathrm{\&pi;})---\left(3\right)$

第四步，从频散关系出发，分析任意广角波动方程的成像精度，利用非线性牛顿法求解任意广角波动方程参考速度的最优化问题，得到任意广角波动方程的最优化参考速度。

时空域上行任意广角波动方程为

$\frac{{\&PartialD;}^{2}u}{\&PartialD;z\&PartialD;t}-\frac{1}{c}({\mathrm{\&Lambda;}}_1+{\mathrm{\&Lambda;}}_2)\frac{{\&PartialD;}^{2}u}{\&PartialD;t^2}+c{\mathrm{\&Lambda;}}_2\frac{{\&PartialD;}^{2}u}{\&PartialD;x^2}=0---\left(4\right)$

其中，d＝(1 0 0 … 0)^T，u＝(u u₁ u₂ … u_n-1)^T，T表示转置，u表示波场值，u₁,u₂,…,u_n-1是推导(4)式时引入的辅助变量；n-1表示辅助变量的个数，x,z分别表示空间水平和深度坐标，t表示时间变量，c是地震波在地下介质中的传播速度，

c₁,c₂,…,c_n-1,c_n是推导(4)式引入的参考速度。

波动方程的理论频散关系为

$k_x^2+k_z^2=\frac{{\mathrm{\&omega;}}^2}{c^2}---\left(7\right)$

其中k_x是水平波数，k_z是垂直波数，ω是角频率。令X＝ck_x/ω＝sinθ，R＝ck_z/ω，则方程(7)可以简化为

$R=\sqrt{1-X^2}=\cos \mathrm{\&theta;}---\left(8\right)$

其中θ表示波的传播方向与界面法线之间的夹角。对方程(4)中的x、z、t做傅里叶变换，得到频率—波数域中的任意广角波动方程

$k_z\mathrm{du}-\frac{\mathrm{\&omega;}}{c}\{{\mathrm{\&Lambda;}}_1+[1-{\left(\frac{c{k}_x}{\mathrm{\&omega;}}\right)}^2\left]{\mathrm{\&Lambda;}}_2\right\}u=0---\left(9\right)$

整理方程(9)可以得到任意广角波动方程对平方根算子的近似函数F_n(θ,c₁,…,c_n)。当选取两个参考速度，并归一化速度后得到

$F_2(\mathrm{\&theta;},c_1,c_2)=(\frac{1}{2c_1}+\frac{c_1}{2})-\frac{c_1}{2}{\sin }^2\mathrm{\&theta;}-\frac{{[(\frac{1}{c_1}-c_1)+c_1{\sin }^2\mathrm{\&theta;}]}^2}{2[(\frac{1}{c_1}+\frac{1}{c_2})+c_1+c_2]-2(c_1+c_2){\sin }^2\mathrm{\&theta;}}.---\left(10\right)$

从平方根算子表达的频散关系出发，对参考速度进行优化。优化的准则是在最小平方意义下，使任意广角波动方程对平方根算子的近似与准确的平方根算子达到最佳的匹配。下面在两个参考速度的情况下，详细地说明参数的优化过程。(8)式与(10)式之间的相对近似误差为

$E_r=F_2-R=(\frac{1}{2c_1}+\frac{c_1}{2})-\frac{c_1}{2}{\sin }^2\mathrm{\&theta;}-\frac{{[(\frac{1}{c_1}-c_1)+c_1{\sin }^2\mathrm{\&theta;}]}^2}{2[(\frac{1}{c_1}+\frac{1}{c_2})+c_1+c_2]-2(c_1+c_2){\sin }^2\mathrm{\&theta;}}-\cos \mathrm{\&theta;}.---\left(11\right)$

定义目标函数

$\mathrm{Ir}(c_1,c_2)={\&Integral;}_{{\mathrm{\&theta;}}_1}^{{\mathrm{\&theta;}}_2}{E}_r^2\left(\mathrm{\&theta;}\right)\mathrm{d\&theta;}---\left(12\right)$

其中θ₁和θ₂是所选择的角度范围，通常选取θ₁＝0°，θ₂＝90°。通过求取目标函数Ir的最小值，得到参考速度c₁和c₂的值。在最小平方意义约束下，得到

$\left\{\begin{array}{c}\frac{\&PartialD;\mathrm{Ir}(c_1,c_2)}{\&PartialD;c_1}=0\\ \frac{\&PartialD;\mathrm{Ir}(c_1,c_2)}{\&PartialD;c_2}=0\end{array}\right.---\left(13\right)$

采用牛顿法求解非线性方程组(13)。首先给定c₁和c₂的初值

和

然后采用(14)式更新参考速度。

$\left\{\begin{array}{c}{c}_1^{(k+1)}=c_1^{\left(k\right)}+\mathrm{\&Delta;}{c_1}^{\left(k\right)}\\ c_2^{(k+1)}=c_2^{\left(k\right)}+\mathrm{\&Delta;}{c_2}^{\left(k\right)}\end{array}\right.,k=\mathrm{0,1,2}...---\left(14\right)$

其中

表示第k次的参考速度迭代值，Δc₁ ^(k)，Δc₂ ^(k)表示第k次的参考速度更正值，是第k+1次的参考速度迭代值，参考速度的更正值Δc₁ ^(k)，Δc₂ ^(k)由(15)式求出。

其中Δc^(k)＝[Δc₁ ^(k),Δc₂ ^(k)]是参考速度更正向量，g(c^(k))＝[g(c₁),g(c₂)]^(k)是目标函数的梯度向量， $H^{\left(k\right)}={\left[\begin{array}{cc}{h}_{11}& h_{12}\\ h_{21}& h_{22}\end{array}\right]}^{\left(k\right)}$ 是目标函数的Hessian矩阵，梯度向量和Hessian矩阵中各元素由(16)—(20)式定义，即

$g\left(c_1\right)=\frac{\&PartialD;\mathrm{Ir}}{\&PartialD;c_1}={\&Integral;}_{{\mathrm{\&theta;}}_1}^{{\mathrm{\&theta;}}_2}2E_r\left(\mathrm{\&theta;}\right)\frac{\&PartialD;E_r}{\&PartialD;c_1}\mathrm{d\&theta;}---\left(16\right)$

$g\left(c_2\right)=\frac{\&PartialD;\mathrm{Ir}}{\&PartialD;c_2}={\&Integral;}_{{\mathrm{\&theta;}}_1}^{{\mathrm{\&theta;}}_2}2E_r\left(\mathrm{\&theta;}\right)\frac{\&PartialD;E_r}{\&PartialD;c_2}\mathrm{d\&theta;}---\left(17\right)$

$h_{11}=\frac{{\&PartialD;}^2\mathrm{Ir}}{\&PartialD;{c_1}^2}=\frac{\&PartialD;g_1}{\&PartialD;c_1}=2{\&Integral;}_{{\mathrm{\&theta;}}_1}^{{\mathrm{\&theta;}}_2}[{\left(\frac{\&PartialD;E_r}{\&PartialD;c_1}\right)}^2+E_r\left(\mathrm{\&theta;}\right)\frac{\&PartialD;E_r^2}{\&PartialD;{c_1}^2}]\mathrm{d\&theta;},---\left(18\right)$

$h_{12}=h_{21}=\frac{{\&PartialD;}^2\mathrm{Ir}}{\&PartialD;c_1\&PartialD;c_2}=\frac{\&PartialD;g_1}{\&PartialD;c_2}=2{\&Integral;}_{{\mathrm{\&theta;}}_1}^{{\mathrm{\&theta;}}_2}[\left(\frac{\&PartialD;E_r}{\&PartialD;c_1}\right)\left(\frac{\&PartialD;E_r}{\&PartialD;c_2}\right)+E_r\left(\mathrm{\&theta;}\right)\frac{\&PartialD;E_r^2}{\&PartialD;c_1\&PartialD;c_2}]\mathrm{d\&theta;},---\left(19\right)$

$h_{22}=\frac{{\&PartialD;}^2\mathrm{Ir}}{\&PartialD;{c_2}^2}=\frac{\&PartialD;g_2}{\&PartialD;c_2}=2{\&Integral;}_{{\mathrm{\&theta;}}_1}^{{\mathrm{\&theta;}}_2}[{\left(\frac{\&PartialD;E_r}{\&PartialD;c_2}\right)}^2+E_r\left(\mathrm{\&theta;}\right)\frac{\&PartialD;E_r^2}{\&PartialD;{c_2}^2}]\mathrm{d\&theta;}.---\left(20\right)$

把(11)式代入到(16)—(20)式中，再根据(15)式就可以得到每一次的参考速度更正值。

当参考速度的更正值小于给定的值（本发明选取10^-6）时，终止迭代，并输出此时的参考速度。输出的参考速度可以使任意广角波动方程的精度在所选择的角度范围内达到全局最优。如图5(a)所示，为本发明应用实例优化二参数的任意广角波动方程的频散关系曲线示意图，如图5(b)所示，为本发明应用实例优化三参数的任意广角波动方程的频散关系曲线示意图，可以看出方程的成像倾角得到了大大地提高。

第五步，采用有限差分方法求解最优化参考速度的任意广角波动方程。

将第四步得到的参考速度带入到方程(4)中，然后采用有限差分方法求解方程(4)，边界处采用PML吸收边界条件，同时将震源波场和记录波场沿着深度方向做波场延拓。

第六步，采用互相关成像条件进行偏移成像。如图6所示，为本发明应用实例复杂介质速度模型示意图，如图7所示，为本发明应用实例复杂介质速度模型的多炮偏移成像结果示意图，可以看出主要的速度界面都得到了较好的成像。

本发明应用实例与现有技术相比，具有如下的有益效果：(1)提出的倾角滤波方法可以有效地压制任意广角波动方程叠前深度偏移的倏逝波成像噪音，提高成像质量；(2)优化了任意广角波动方程的参考速度参数，提高了任意广角波动方程叠前深度偏移成像精度，因此使基于任意广角波动方程的叠前深度偏移方法可以对更陡倾角的地层成像。

本领域技术人员还可以了解到本发明实施例列出的各种说明性逻辑块（illustrativelogical block），单元，和步骤可以通过电子硬件、电脑软件，或两者的结合进行实现。为清楚展示硬件和软件的可替换性（interchangeability），上述的各种说明性部件（illustrativecomponents），单元和步骤已经通用地描述了它们的功能。这样的功能是通过硬件还是软件来实现取决于特定的应用和整个系统的设计要求。本领域技术人员可以对于每种特定的应用，可以使用各种方法实现所述的功能，但这种实现不应被理解为超出本发明实施例保护的范围。

本发明实施例中所描述的各种说明性的逻辑块，或单元都可以通过通用处理器，数字信号处理器，专用集成电路（ASIC），现场可编程门阵列或其它可编程逻辑装置，离散门或晶体管逻辑，离散硬件部件，或上述任何组合的设计来实现或操作所描述的功能。通用处理器可以为微处理器，可选地，该通用处理器也可以为任何传统的处理器、控制器、微控制器或状态机。处理器也可以通过计算装置的组合来实现，例如数字信号处理器和微处理器，多个微处理器，一个或多个微处理器联合一个数字信号处理器核，或任何其它类似的配置来实现。

本发明实施例中所描述的方法或算法的步骤可以直接嵌入硬件、处理器执行的软件模块、或者这两者的结合。软件模块可以存储于RAM存储器、闪存、ROM存储器、EPROM存储器、EEPROM存储器、寄存器、硬盘、可移动磁盘、CD-ROM或本领域中其它任意形式的存储媒介中。示例性地，存储媒介可以与处理器连接，以使得处理器可以从存储媒介中读取信息，并可以向存储媒介存写信息。可选地，存储媒介还可以集成到处理器中。处理器和存储媒介可以设置于ASIC中，ASIC可以设置于用户终端中。可选地，处理器和存储媒介也可以设置于用户终端中的不同的部件中。

在一个或多个示例性的设计中，本发明实施例所描述的上述功能可以在硬件、软件、固件或这三者的任意组合来实现。如果在软件中实现，这些功能可以存储与电脑可读的媒介上，或以一个或多个指令或代码形式传输于电脑可读的媒介上。电脑可读媒介包括电脑存储媒介和便于使得让电脑程序从一个地方转移到其它地方的通信媒介。存储媒介可以是任何通用或特殊电脑可以接入访问的可用媒体。例如，这样的电脑可读媒体可以包括但不限于RAM、ROM、EEPROM、CD-ROM或其它光盘存储、磁盘存储或其它磁性存储装置，或其它任何可以用于承载或存储以指令或数据结构和其它可被通用或特殊电脑、或通用或特殊处理器读取形式的程序代码的媒介。此外，任何连接都可以被适当地定义为电脑可读媒介，例如，如果软件是从一个网站站点、服务器或其它远程资源通过一个同轴电缆、光纤电缆、双绞线、数字用户线（DSL）或以例如红外、无线和微波等无线方式传输的也被包含在所定义的电脑可读媒介中。所述的碟片（disk）和磁盘（disc）包括压缩磁盘、镭射盘、光盘、DVD、软盘和蓝光光盘，磁盘通常以磁性复制数据，而碟片通常以激光进行光学复制数据。上述的组合也可以包含在电脑可读媒介中。

以上所述的具体实施方式，对本发明的目的、技术方案和有益效果进行了进一步详细说明，所应理解的是，以上所述仅为本发明的具体实施方式而已，并不用于限定本发明的保护范围，凡在本发明的精神和原则之内，所做的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。

Claims (10)

1.一种基于任意广角波动方程的成像方法，其特征在于，所述基于任意广角波动方程的成像的方法，包括：

获取建立的速度模型和给定的震源波场；

利用正演模拟得到速度模型的合成地震记录；

对正演模拟的合成地震记录和给定的震源波场做倾角滤波处理；

从频散关系出发，分析任意广角波动方程的成像精度，得到任意广角波动方程的最优化参考速度；

采用有限差分方法求解最优化参考速度的任意广角波动方程；

采用互相关成像条件进行偏移成像。

2.如权利要求1所述基于任意广角波动方程的成像方法，其特征在于，所述利用正演模拟得到速度模型的合成地震记录，包括：

采用交错网格有限差分方法求解时空域各向同性均匀介质的一阶速度—应力声波波动方程，从而得到所述速度模型的合成地震记录。

3.如权利要求1所述基于任意广角波动方程的成像方法，其特征在于，所述对正演模拟的合成地震记录和给定的震源波场做倾角滤波处理，包括：

利用频率波数域的扇形滤波器对正演模拟的合成地震记录和给定的震源波场做倾角滤波处理。

4.如权利要求1所述基于任意广角波动方程的成像方法，其特征在于，所述从频散关系出发，分析任意广角波动方程的成像精度，得到任意广角波动方程的最优化参考速度，包括：

从频散关系出发，分析任意广角波动方程的成像精度，利用非线性牛顿法得到任意广角波动方程的最优化参考速度。

5.如权利要求1所述基于任意广角波动方程的成像方法，其特征在于，所述采用有限差分方法求解最优化参考速度的任意广角波动方程，包括：

采用有限差分方法求解最优化参考速度的任意广角波动方程，并在边界处采用完美匹配层PML吸收边界条件。

6.一种基于任意广角波动方程的成像装置，其特征在于，所述基于任意广角波动方程的成像的装置，包括：

获取单元，用于获取建立的速度模型和给定的震源波场；

正演模拟单元，用于利用正演模拟得到速度模型的合成地震记录；

倾角滤波单元，用于对正演模拟的合成地震记录和给定的震源波场做倾角滤波处理；

分析单元，用于从频散关系出发，分析任意广角波动方程的成像精度，得到任意广角波动方程的最优化参考速度；

方程求解单元，用于采用有限差分方法求解最优化参考速度的任意广角波动方程；

成像单元，用于采用互相关成像条件进行偏移成像。

7.如权利要求6所述基于任意广角波动方程的成像装置，其特征在于，

所述正演模拟单元，进一步用于采用交错网格有限差分方法求解时空域各向同性均匀介质的一阶速度—应力声波波动方程，从而得到所述速度模型的合成地震记录。

8.如权利要求6所述基于任意广角波动方程的成像装置，其特征在于，

所述倾角滤波单元，进一步用于利用频率波数域的扇形滤波器对正演模拟的合成地震记录和给定的震源波场做倾角滤波处理。

9.如权利要求6所述基于任意广角波动方程的成像装置，其特征在于，

所述分析单元，进一步用于从频散关系出发，分析任意广角波动方程的成像精度，利用非线性牛顿法得到任意广角波动方程的最优化参考速度。

10.如权利要求6所述基于任意广角波动方程的成像装置，其特征在于，

所述方程求解单元，进一步用于采用有限差分方法求解最优化参考速度的任意广角波动方程，并在边界处采用完美匹配层PML吸收边界条件。

Images