CN105137486A - 各向异性介质中弹性波逆时偏移成像方法及其装置 - Google Patents

Abstract

本发明提供了一种各向异性介质中弹性波逆时偏移成像方法及其装置，所述方法包括：构建地下炮点正传波场；利用变异函数IDW算法对所述炮点正传波场进行分离，并存储分离后的所述炮点正传波场的有效边界；以所述存储的有效边界作为边界，对分离后的炮点正传波场进行逆时反传重构；获取地震数据中的矢量数据，并以所述矢量数据作为边界条件进行检波点逆时反传模拟；利用变异函数IDW算法对所述检波点反传波场进行分离；利用归一化成像条件对所述炮点反传重构波场和分离后的所述检波点波场进行成像，得到地下成像数据。本发明可有效提高成像效率。

Description

各向异性介质中弹性波逆时偏移成像方法及其装置

技术领域

本发明涉及油气勘探技术领域，特别涉及一种各向异性介质中弹性波逆时偏移成像方法及其装置。

背景技术

弹性波逆时偏移和声波逆时偏移有相似的流程，因为逆时偏移的核心是求解波动方程，所以，弹性波逆时偏移就是基于弹性波动方程重建地下炮点和检波点的波场，并以检波点处的矢量地震数据作为边界条件，通过成像条件提取反射系数成像。

弹性波逆时偏移重建的地下波场包含多种波场信息，如果只是简单地将垂直分量当做纵波，水平分量当成横波，然后直接应用成像条件进行成像，会造成成像结果的耦合性串扰和假象，且没有物理上的实际意义。

目前弹性波逆时偏移存在着两种策略，第一种策略基于标量处理理论，即首先对检波点处的地震数据进行波场分离，然后分别应用标量波动方程理论进行逆时偏移处理。第二种策略是通过各种数值模拟方法，比如有限差分法，重建地下某一点的矢量波场，然后应用成像条件，分别从炮点和检波点的正传波场和反传波场处得到反射系数信息，之后在成像区域重复流程，重建地下构造。

对于第二种策略，成像条件的选择非常重要，由于炮点正传波场和检波点反传波场都是矢量场，因此弹性波的成像条件比声波成像条件更加复杂，传统的互相关成像条件对弹性波矢量场各分量做互相关运算，各种波场混杂在一起，造成耦合性串扰和假象，应用弹性势的成像条件，基于Helmholtz定理将波场分为纵波势和横波势的表达，然后应用互相关成像条件对其进行成像，这样就可以显著消弱不同波场间的耦合性串扰。

第二种策略保持了弹性波场的矢量特征，直接矢量输入，矢量处理。虽然，该策略基于弹性势的成像条件，可以避免波场之间的串扰和假象，但是会极大地增加计算量和存储量，特别是波场重建，波场分离部分更是消耗了巨大的计算资源，这极大地制约了该技术在复杂二维和三维各向异性介质中的应用。

发明内容

本发明提供一种各向异性介质中弹性波逆时偏移成像方法及其装置，以解决上述一项或多项缺失。

本发明提供一种各向异性介质中弹性波逆时偏移成像方法，所述方法包括：进行炮点正传模拟，以构建地下炮点正传波场；利用变异函数IDW算法对所述炮点正传波场进行分离，得到分离后的所述炮点正传波场，并存储分离后的所述炮点正传波场的有效边界；以计算换存储，以存储的所述有效边界作为边界，对分离后的所述炮点正传波场进行炮点逆时反传重构，得到炮点反传重构波场；获取地震数据中的矢量数据，并以所述矢量数据作为边界条件进行检波点逆时反传模拟，得到地下检波点反传波场；利用所述变异函数IDW算法对所述检波点反传波场进行分离，得到分离后的所述检波点反传波场；利用归一化成像条件对所述炮点反传重构波场和分离后的所述检波点波场进行成像，得到地下成像数据，生成成像剖面并输出。

一个实施例中，所述的方法还包括：对所述成像数据进行低通滤波以压制其中的低波数噪音并输出所述成像剖面。

一个实施例中，利用所述变异函数IDW算法对所述炮点正传波场进行波场分离，包括：利用傅里叶变换分别将各时刻的所述炮点正传波场由空间域变换至波数域；选取N个参考模型，根据所述参考模型计算分离算子，并在波数域利用自褶积组合窗函数对所述分离算子进行截断；在波数域利用所述分离算子对所述参考模型下的所述炮点正传波场进行分离，得到所述参考模型下的分离后的所述炮点正传波场；将所述参考模型下的分离后的所述炮点正传波场由波数域变换至空间域，N个所述参考模型对应得到N个所述参考模型下的分离后的所述炮点正传波场，N为正整数；在空间域，根据所述参考模型的权重系数对N个所述参考模型下的分离后的所述炮点正传波场进行加权插值处理，得到最终分离后的所述炮点正传波场。

一个实施例中，利用所述变异函数IDW算法对所述检波点反传波场进行分离，包括：利用傅里叶变换分别将各时刻的所述检波点反传波场由空间域变换至波数域；选取N个参考模型，根据所述参考模型计算分离算子，并在波数域利用自褶积组合窗函数对所述分离算子进行截断；在波数域利用所述分离算子对所述参考模型下的所述检波点反传波场进行分离，得到所述参考模型下的分离后的所述检波点反传波场；将所述参考模型下的分离后的所述检波点反传波场由波数域变换至空间域，N个所述参考模型对应得到N个所述参考模型下的分离后的所述检波点反传波场，N为正整数；在空间域，根据所述参考模型的权重系数对N个所述参考模型下的分离后的所述检波点反传波场进行加权插值处理，得到最终分离后的所述检波点反传波场。

一个实施例中，选取所述参考模型的方法，包括：遍历待数值模拟初始模型的多个各向异性参数，并利用变异函数的临界值进行约束，选取在所述临界值范围内出现概率最大的参考点，根据所述参考点搜索得到所述参考模型。

一个实施例中，所述参考模型的权重系数通过所述变异函数计算得到。

一个实施例中，所述变异函数的获取方法包括：通过计算所述初始模型的各向异性参数分布的变异函数值，拟合得到所述初始模型的变异函数。

一个实施例中，获取所述自褶积组合窗函数的方法，包括：选择主瓣和旁瓣性能高于设定阈值的窗函数作为原始窗函数；对原始窗函数做L次自褶积运算得到自褶积后的窗函数，其中，L为正整数；对自褶积后的窗函数与原始窗函数进行加权运算，得到所述自褶积组合窗函数。

一个实施例中，在确定需要主瓣性能优先的情况下，原始窗函数的加权系数高于自褶积后的窗函数的加权系数；在确定需要旁瓣性能优先的情况下，自褶积后的函数的加权系数高于原始窗函数的加权系数。

本发明还提供一种各向异性介质中弹性波逆时偏移成像装置，所述装置包括：炮点正传模拟单元，用于进行炮点正传模拟，以构建地下炮点正传波场；炮点正传波场分离单元，用于利用变异函数IDW算法对所述炮点正传波场进行分离，得到分离后的所述炮点正传波场，并存储分离后的所述炮点正传波场的有效边界；炮点正传波场反传重构单元，用于以存储的所述有效边界作为边界条件，对分离后的所述炮点正传波场进行炮点逆时反传重构，得到炮点反传重构波场；检波点逆时反传模拟单元，用于获取地震数据中的矢量数据，并以所述矢量数据作为边界条件进行检波点逆时反传模拟，得到地下检波点反传波场；检波点反传波场分离单元，用于利用所述变异函数IDW算法对所述检波点反传波场进行分离，得到分离后的所述检波点反传波场；成像单元，用于利用归一化成像条件对所述炮点反传重构波场和分离后的所述检波点波场进行成像，得到地下成像数据，生成成像剖面并输出。

一个实施例中，所述装置还包括：滤波单元，用于对所述成像数据进行低通滤波以压制其中的低波数噪音并输出所述成像剖面。

一个实施例中，所述炮点正传波场分离单元，包括：第一空间域至波数域变换模块，用于利用傅里叶变换分别将各时刻的所述炮点正传波场由空间域变换至波数域；第一分离算子生成模块，用于选取N个参考模型，根据所述参考模型计算分离算子，并在波数域利用自褶积组合窗函数对所述分离算子进行截断；第一波场分离模块，用于在波数域利用所述分离算子对所述参考模型下的所述炮点正传波场进行分离，得到所述参考模型下的分离后的所述炮点正传波场；第一波数域至空间域变换模块，用于将所述参考模型下的分离后的所述炮点正传波场由波数域变换至空间域，N个所述参考模型对应得到N个所述参考模型下的分离后的所述炮点正传波场，N为正整数；第一加权插值处理模块，用于在空间域，根据所述参考模型的权重系数对N个所述参考模型下的分离后的所述炮点正传波场进行加权插值处理，得到最终分离后的所述炮点正传波场。

一个实施例中，检波点反传波场分离单元包括：第二空间域至波数域变换模块，用于利用傅里叶变换分别将各时刻的所述检波点反传波场由空间域变换至波数域；第二分离算子生成模块，用于选取N个参考模型，根据所述参考模型计算分离算子，并在波数域利用自褶积组合窗函数对所述分离算子进行截断；第二波场分离模块，用于在波数域利用所述分离算子对所述参考模型下的所述检波点反传波场进行分离，得到所述参考模型下的分离后的所述检波点反传波场；第二波数域至空间域变换模块，用于将所述参考模型下的分离后的所述检波点反传波场由波数域变换至空间域，N个所述参考模型对应得到N个所述参考模型下的分离后的所述检波点反传波场，N为正整数；第二加权插值处理模块，用于在空间域，根据所述参考模型的权重系数对N个所述参考模型下的分离后的所述检波点反传波场进行加权插值处理，得到最终分离后的所述检波点反传波场。

本发明的各向异性介质中弹性波逆时偏移成像方法及其装置，综合运用自褶积组合窗函数，提高弹性波数值模拟的精度和效率。针对弹性波逆时偏移中波场计算量巨大的情况，通过引入变异函数IDW插值算法，在混合域分离波场，极大地提高分离效率。运用有效边界存储策略，不存储炮点分离后的波场，只存储有效边界，在成像时重构炮点分离波场。

附图说明

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

图1是本发明实施例的各向异性介质中弹性波逆时偏移成像方法的流程示意图；

图2是本发明另一实施例的各向异性介质中弹性波逆时偏移成像方法的流程示意图；

图3是本发明一实施例中波场分离方法的流程示意图；

图4是本发明一实施例中波场分离方法的流程示意图；

图5是本发明一实施例中获取自褶积组合窗函数的方法的流程示意图；

图6是本发明一实施例中纵波速度模型的示意图；

图7是本发明一实施例中横波速度模型的示意图；

图8是本发明一实施例中密度模型的示意图；

图9至图12分别是本发明一实施例中PP波、PS波、SP波、SS波的单炮成像结果示意图；

图13是本发明实施例的各向异性介质中弹性波逆时偏移成像装置的结构示意图；

图14是本发明一实施例的各向异性介质中弹性波逆时偏移成像装置的结构示意图；

图15是本发明一实施例中波场分离单元的结构示意图；

图16是本发明又一实施例中波场分离单元的结构示意图。

具体实施方式

为使本发明实施例的目的、技术方案和优点更加清楚明白，下面结合附图对本发明实施例做进一步详细说明。在此，本发明的示意性实施例及其说明用于解释本发明，但并不作为对本发明的限定。

为了避免波场之间的耦合性串扰，在弹性波逆时偏移中采用基于弹性势的成像条件，需要将数值模拟得到的矢量波场分解为纵波波场和横波波场。然而，利用偏振特性分离波场的策略需要极大的计算量。假设用来分离波场的拟微分算子/分离算子的大小是n²，模型的大小是N²，那么在空间域每个网格点对模型进行波场分离，相当于空间滤波，需要n²次乘法运算，n²-1次加法运算，计算量是(2n²-1)N²，这远远大于m阶精度的有限差分法(2m-1)N²的计算量。拟微分算子分离矢量波场只有在n²≥51*51时，才有较好的分离效果。常规有限差分算子，即使采用24阶精度，计算量也远远小于拟微分算子分离矢量波场的计算量。

因此，发明人考虑采用混合域分离的方式，即在波数域分离矢量波场，在空间域采用插值的方式重构分离波场，可以有效降低计算量。

逆时偏移的原理决定了需要保存所有时刻的炮点正传矢量波场，对于弹性波来说，需要同时保存P波和S波的波场，这对计算机的硬件有非常大的需求，主要体现在容量和I/O(输入/输出)。对于小模型勉强可以存储所有时刻的波场，但对于复杂模型，甚至是三维模型，存储所有时刻的波场是不可行的。所以，弹性波逆时偏移另一个需要解决的问题就是存储问题。

因此，发明人考虑，在边界存储策略的基础上，采用有效边界存储策略，以进一步减少边界存储量。在构建炮点正传波场时，不需要存储每个时刻的波场，极大地减少了I/O(输入/输出)。

图1是本发明实施例的各向异性介质中弹性波逆时偏移成像方法的流程示意图。如图1所示，本发明实施例的各向异性介质中弹性波逆时偏移成像方法，包括步骤：

S101：进行炮点正传模拟，以构建地下炮点正传波场；

S102：利用变异函数IDW(InverseDistanceWeighted，反距离权重)算法对所述炮点正传波场进行分离，得到分离后的所述炮点正传波场，并存储分离后的所述炮点正传波场的有效边界；

S103：以存储的所述有效边界作为边界条件，对分离后的所述炮点正传波场进行炮点逆时反传重构，得到炮点反传重构波场；

S104：获取地震数据中的矢量数据，并以所述矢量数据作为边界条件进行检波点逆时反传模拟，得到地下检波点反传波场；

S105：利用所述变异函数IDW算法对所述检波点反传波场进行分离，得到分离后的所述检波点反传波场；

S106：利用归一化成像条件对所述炮点反传重构波场和分离后的所述检波点波场进行成像，得到地下成像数据，生成成像剖面并输出。

上述步骤S102和S105中，该变异函数IDW算法为通过IDW算法优化的变异函数算法。在其他实施例中，还可以选择多种其他插值算法进行波场分离。

对于复杂模型，甚至是三维模型，存储所有时刻的波场是不可行的，计算量巨大，因此不存储炮点正传波场而只是将其有效边界提取出来，并该有效边界作为边界条件进行炮点逆时反传重构。

本发明实施例的各向异性介质中弹性波逆时偏移成像方法，不存储所有时刻的边界，而是存储有效边界，可显著减少边界存储量；本发明实施例，通过运用变异函数IDW算法插值运算分离波场可有效提高分离效率。如此一来，各向异性介质中弹性波的逆时偏移成像的效率得到显著提高。

针对计算量问题，还可以将各向异性的非均匀介质近似为均匀介质，以降低逆时偏移成像的计算量；针对存储量问题，还可以采用随机边界条件等构建炮点反传波场。

图2是本发明另一实施例的各向异性介质中弹性波逆时偏移成像方法的流程示意图。如图2所示，图1所示的各向异性介质中弹性波逆时偏移成像方法，还可包括步骤：

S107：对所述成像数据进行低通滤波以压制其中的低波数噪音并输出所述成像剖面。

本发明实施例中，通过对成像数据进行低通滤波，压制低波数噪声，可以提高逆时偏移成像的品质。

图3是本发明一实施例中波场分离方法的流程示意图。如图3所示，图1所示的各向异性介质中弹性波逆时偏移成像方法，步骤S102中对炮点正传波场进行波场分离，利用上述变异函数IDW算法进行波场分离的方法，可包括步骤：

S301：利用傅里叶变换分别将各时刻的炮点正传波场由空间域变换至波数域；

S302：选取N个参考模型，根据所述参考模型计算分离算子，并在波数域利用自褶积组合窗函数对所述分离算子进行截断；

S303：在波数域利用所述分离算子对所述参考模型下的所述炮点正传波场进行分离，得到所述参考模型下的分离后的所述炮点正传波场；

S304：将所述参考模型下的分离后的所述炮点正传波场由波数域变换至空间域，N个所述参考模型对应得到N个所述参考模型下的分离后的所述炮点正传波场，N为正整数；

S305：在空间域，根据所述参考模型的权重系数对N个所述参考模型下的分离后的所述炮点正传波场进行加权插值处理，得到最终分离后的所述炮点正传波场。

图4是本发明一实施例中波场分离方法的流程示意图。如图4所示，图1所示的各向异性介质中弹性波逆时偏移成像方法，上述步骤S105中对检波点反传波场进行波场分离，利用上述变异函数IDW算法进行波场分离的方法，可包括步骤：

S401：利用傅里叶变换分别将各时刻的检波点反传波场由空间域变换至波数域；

S402：选取N个参考模型，根据所述参考模型计算分离算子，并在波数域利用自褶积组合窗函数对所述分离算子进行截断；

S403：在波数域利用所述分离算子对所述参考模型下的所述检波点反传波场进行分离，得到所述参考模型下的分离后的所述检波点反传波场；

S404：将所述参考模型下的分离后的所述检波点反传波场由波数域变换至空间域，N个所述参考模型对应得到N个所述参考模型下的分离后的所述检波点反传波场，N为正整数；

S405：在空间域，根据所述参考模型的权重系数对N个所述参考模型下的分离后的所述检波点反传波场进行加权插值处理，得到最终分离后的所述检波点反传波场。

本发明实施例中，在波数域利用分离算子分离波场，在空间域运用变异函数IDW算法插值运算得到最终分离后的波场，可以有效降低波场分离的计算复杂度。本发明实施例中，在波数域采用自褶积组合窗函数截断分离算子，再分离波场，可以提高波场分离的效率，可以保证截断得到的分离算子的准确度更高。

一个实施例中，可以通过遍历待数值模拟初始模型的多个各向异性参数，并利用变异函数的临界值R进行约束，选取在所述临界值范围内出现概率最大的参考点，根据所述参考点搜索得到上述参考模型。

上述各项异性参数可以包括ε、δ、θ三个参数，其中，参数ε、δ表示横向各向同性TI介质各向异性的系数，参数θ表示具有倾斜对称轴的横向各向同性TTI介质对称轴的倾角。

本发明实施例中，选取参考模型时，不仅考虑到两点之间的距离也考虑了空间位置关系，使得进行波场分离时插值精度和效率都较高。

一个实施例中，上述参考模型的权重系数可以通过上述变异函数计算得到。

本发明实施例中，采用变异函数和参考点搜索策略选取参考模型，再利用变异函数计算得到各个参考模型相对于初始模型的权重系数，使得计算得到的权重系数较为合理，达到了在保证波场分离的分离结果准确度的基础上，减少了计算量的效果。

一个实施例中，可以通过计算上述初始模型的各向异性参数分布的变异函数值，然后拟合得到上述初始模型的变异函数。

例如，按照以下公式计算初始模型各向异性参数分布的变异函数值：

$\mathrm{\μ}\left(\stackrel{-}{h}\right)=\frac{1}{2n}\underset{i-1}{\overset{n}{\Σ}}{\[\mathrm{\Φ}(\mathrm{\ϵ},\mathrm{\δ},\mathrm{\θ})-\mathrm{\Φ}(\mathrm{\ϵ}+h_{\mathrm{\ϵ}},\mathrm{\δ}+h_{\mathrm{\δ}},\mathrm{\θ}+h_{\mathrm{\θ}})\]}^2,$

其中，μ(h)表示变异函数值，Φ＝f[δ+2(ε-δ)sin₂(α-θ)]sin2(α-θ)，δ表示横向各向同性TI介质各向异性的系数，θ表示具有倾斜对称轴的横向各向同性TTI介质对称轴的倾角，表示插值参考点的位置增量，i＝1,2...n，其中，n表示插值参考点的个数。

例如，按照以下公式计算各个参考模型相对于初始模型的权重系数：

$w^k=\frac{1/\mathrm{\μ}({\left({\mathrm{\ϵ}}_k-\mathrm{\ϵ}\right)}^2+{\left({\mathrm{\δ}}_k-\mathrm{\δ}\right)}^2+{\left({\mathrm{\θ}}_k-\mathrm{\θ}\right)}^2)}{\underset{k=1}{\overset{n}{\Σ}}1/\mathrm{\μ}({\left({\mathrm{\ϵ}}_k-\mathrm{\ϵ}\right)}^2+{\left({\mathrm{\δ}}_k-\mathrm{\δ}\right)}^2+{\left({\mathrm{\θ}}_k-\mathrm{\θ}\right)}^2)},$

其中，w^k表示权重系数，(ε_k,δ_k,θ_k)表示插值参考点。

例如，插值参考点的选择方法可包括：根据所述初始模型计算得到多个分散的变异函数值；对得到的多个分散的变异函数值进行拟合，得到理论变异函数模型；根据所述变异函数模型确定变异函数中的各个参数的参数值；遍历待数值模拟初始模型中不同数值的ε、δ、θ出现的概率，选取在拟合的变异函数临界值之内的，出现概率最大的点作为插值参考点。

图5是本发明一实施例中获取自褶积组合窗函数的方法的流程示意图。如图5所示，图3及图4所示的波场分离的方法中，获取自褶积组合窗函数的方法，可包括步骤：

S501：选择主瓣和旁瓣性能高于设定阈值的窗函数作为原始窗函数；

S502：对原始窗函数做L次自褶积运算得到自褶积后的窗函数，其中，L为正整数；

S503：对自褶积后的窗函数与原始窗函数进行加权运算，得到所述自褶积组合窗函数。

在上述实施例中，对选定的原始窗函数进行自褶积运算，自褶积运算使得原始窗函数的主瓣性能变差，旁瓣性能变好，然后根据当前的需要对原始窗函数和自褶积后的窗函数进行加权运算，从而得到一个满足要求的窗函数。利用该种自褶积组合窗函数的到的分离算子可以提高波场分离的精度和稳定性。

一个实施例中，在上述步骤S502中，进行加权运算的加权系数可以具备特点：在确定需要主瓣性能优先的情况下，原始窗函数的加权系数高于自褶积后的窗函数的加权系数；在确定需要旁瓣性能优先的情况下，自褶积后的函数的加权系数高于原始窗函数的加权系数。以此在维持适当的主瓣宽度的前提下，可以尽可能增大旁瓣衰减。

为使本领域的技术人员更容易理解本发明的技术方案，下面将以一具体实施例说明本发明的构思、技术方案及其优点，但并非限定本发明的保护范围。

从本构方程、运动微分方程及几何方程出发采用规则网格有限差分算法实现弹性波各向异性数值模拟。在线性弹性理论的假设下，可以建立应变张量和位移之间的联系，可得几何方程：

${\mathrm{\ϵ}}_{kl}=\frac{1}{2}\[{\∂}_k{u}_l+{\∂}_l{u}_k\],k,l=1,2,3,---\left(1\right)$

在公式(1)中，ε_kl是应变张量，u_l和u_k代表位移分量，是偏导数算子。

在线性弹性体内，应变张量和应力张量有如下的本构关系：

σ_ij＝c_ijklε_kl，(2)

在公式(2)中，σ_ij为应力张量，c_ijkl为弹性系数张量。

从应力表示的动力学平衡方程出发，不考虑体积力的影响，可得到弹性动力学方程，即运动微分方程：

$\mathrm{\ρ}{\∂}_{tt}^2{u}_i={\∂}_j{\mathrm{\σ}}_{ij}.---\left(3\right)$

采用时间2阶，空间2N阶的规则有限差分格式，时间二阶导数写和空间2N阶格式分别表示为：

$\begin{array}{c}\frac{{\∂}^2{u}_x}{\∂t^2}=\frac{u_x^{n+1}-2u_x^n+u_x^{n-1}}{{\mathrm{\Δt}}^2},\\ \frac{{\∂}^2{u}_z}{\∂t^2}=\frac{u_z^{n+1}-2u_z^n+u_z^{n-1}}{{\mathrm{\Δt}}^2},\end{array}---\left(4\right)$

$\begin{array}{c}\frac{\∂u_x}{\∂x}{|}_{x=0}=\frac{1}{\mathrm{\Δ}x}\underset{m=1}{\overset{N/2}{\Σ}}{a}_m(u_x^{i+m,k,n}-u_x^{i-m,k,n}),\\ \frac{\∂u_z}{\∂z}{|}_{z=0}=\frac{1}{\mathrm{\Δ}z}\underset{m=1}{\overset{N/2}{\Σ}}{a}_m(u_z^{i,k+m,n}-u_z^{i,k-m,n}),\end{array}---\left(5\right)$

其中，u_x和u_z是离散的时间域波场，为时间域的优化系数，w(m)为窗函数。Δx,Δz为空间步长，i,k表示空间位置，n表示时间步。从公式(4)、(5)可以看出，每求解一次导数，需要用到周围N+1个网格点的波场值。

首先应用几何方程公式(1)，由位移计算应变，再根据本构方程公式(2)，计算应力，最后将应力代入运动微分方程公式(3)，通过时间步迭代，计算下一时刻的位移值，直到完成所有时间的计算。最后的弹性波数值模拟公式为：

$\begin{array}{c}{u}_x^{i,k,n+1}=2u_x^{i,k,n}-u_x^{i,k,n-1}+{\mathrm{\Δt}}^2\left[\begin{array}{c}\frac{1}{\mathrm{\Δ}x}\underset{m=1}{\overset{N/2}{\Σ}}{a}_m({\mathrm{\σ}}_{xx}^{i+m,k,n}-{\mathrm{\σ}}_{xx}^{i-m,k,n})+\mathrm{...}\\ \frac{1}{\mathrm{\Δ}z}\underset{m=1}{\overset{N/2}{\Σ}}{a}_m({\mathrm{\σ}}_{xz}^{i,k+m,n}-{\mathrm{\σ}}_{xz}^{i,k-m,n})\end{array}\right],\\ u_z^{i,k,n+1}=2u_z^{i,k,n}-u_z^{i,k,n-1}+{\mathrm{\Δt}}^2\left[\begin{array}{c}\frac{1}{\mathrm{\Δ}x}\underset{m=1}{\overset{N/2}{\Σ}}{a}_m({\mathrm{\σ}}_{zx}^{i+m,k,n}-{\mathrm{\σ}}_{zx}^{i-m,k,n})+\mathrm{...}\\ \frac{1}{\mathrm{\Δ}z}\underset{m=1}{\overset{N/2}{\Σ}}{a}_m({\mathrm{\σ}}_{zz}^{i,k+m,n}-{\mathrm{\σ}}_{zz}^{i,k-m,n})\end{array}\right],\end{array}---\left(6\right)$

为了便于重复使用，可以将弹性系数矩阵的独立变量，存储在GPU的寄存器(register)中，加快读取速度。因为每个炮之间是天然解耦的，所以每个GPU节点处理一炮数据，一炮数据的数值模拟，涉及多个GPU模块，包括空间高阶有限差分模块，吸收边界GPU模块，输出记录GPU模块等。

同理，三维情况下的弹性波数值模拟公式为：

$\begin{array}{c}{u}_x^{i,j,k,n+1}=2u_x^{i,j,k,n}-u_x^{i,j,k,n-1}+{\mathrm{\Δt}}^2\left[\begin{array}{c}\frac{1}{\mathrm{\Δ}x}\underset{m=1}{\overset{N/2}{\Σ}}{a}_m({\mathrm{\σ}}_{xx}^{i+m,j,k,n}-{\mathrm{\σ}}_{xx}^{i-m,j,k,n})+\mathrm{...}\\ \frac{1}{\mathrm{\Δ}y}\underset{m=1}{\overset{N/2}{\Σ}}{a}_m({\mathrm{\σ}}_{xy}^{i,j+m,k,n}-{\mathrm{\σ}}_{xy}^{i,j-m,k,n})+\mathrm{...}\\ \frac{1}{\mathrm{\Δ}z}\underset{m=1}{\overset{N/2}{\Σ}}{a}_m({\mathrm{\σ}}_{xz}^{i,j,k+m,n}-{\mathrm{\σ}}_{xz}^{i,j,k-m,n})\end{array}\right],\\ u_y^{i,j,k,n+1}=2u_y^{i,j,k,n}-u_y^{i,j,k,n-1}+{\mathrm{\Δt}}^2\left[\begin{array}{c}\frac{1}{\mathrm{\Δ}y}\underset{m=1}{\overset{N/2}{\Σ}}{a}_m({\mathrm{\σ}}_{yx}^{i+m,j,k,n}-{\mathrm{\σ}}_{yx}^{i-m,j,k,n})+\mathrm{...}\\ \frac{1}{\mathrm{\Δ}y}\underset{m=1}{\overset{N/2}{\Σ}}{a}_m({\mathrm{\σ}}_{yy}^{i,j+m,k,n}-{\mathrm{\σ}}_{yy}^{i,j-m,k,n})+\mathrm{...}\\ \frac{1}{\mathrm{\Δ}z}\underset{m=1}{\overset{N/2}{\Σ}}{a}_m({\mathrm{\σ}}_{yz}^{i,j,k+m,n}-{\mathrm{\σ}}_{yz}^{i,j,k-m,n})\end{array}\right],\\ u_z^{i,j,k,n+1}=2u_z^{i,j,k,n}-u_z^{i,j,k,n-1}+{\mathrm{\Δt}}^2\left[\begin{array}{c}\frac{1}{\mathrm{\Δ}z}\underset{m=1}{\overset{N/2}{\Σ}}{a}_m({\mathrm{\σ}}_{zx}^{i+m,j,k,n}-{\mathrm{\σ}}_{zx}^{i-m,j,k,n})+\mathrm{...}\\ \frac{1}{\mathrm{\Δ}y}\underset{m=1}{\overset{N/2}{\Σ}}{a}_m({\mathrm{\σ}}_{zy}^{i,j+m,k,n}-{\mathrm{\σ}}_{zy}^{i,j-m,k,n})+\mathrm{...}\\ \frac{1}{\mathrm{\Δ}z}\underset{m=1}{\overset{N/2}{\Σ}}{a}_m({\mathrm{\σ}}_{zz}^{i,j,k+m,n}-{\mathrm{\σ}}_{zz}^{i,j,k-m,n})\end{array}\right],\end{array}---\left(7\right)$

假设采用12阶的有限差分格式，可以看到，每计算一个波场的值，需要读取周围36个网格点的应力值，再写入，读写冗余度达到38，另外，需要计算波场的三个分量，造成计算消耗非常大。为此，考虑在这一个环节使用共享内存(sharedmemory)，提高IO效率。但是，考虑到共享内存的容量有限，无法把整个三维数据体存储进共享内存，可以将x，y维的数据存储到共享内存中，最慢维z方向的网格点存储在寄存器中(register)中，然后在z方向循环，直到遍历所有的网格点。如果计算区域为M×N，可以设置共享内存大小为(M+12)×(N+12)。

在波数域采用自褶积组合窗函数截断逼近拟微分算子/分离算子，分离波场，然后反变换回空间域，运用变异函数IDW算法插值实现最终波场的分离。计算权重时，不仅考虑到两点之间的距离也考虑了空间位置关系，插值精度和效率都较高。

TI各向异性介质中，地震波的偏振向量的各向异性部分由线性部分和非线性部分构成，对于TTI介质，去掉非线性部分，也就是在弱各向异性的条件下，quasi-P波的偏振角可以被近似表示为：

v_p＝α+f[δ+2(ε-δ)sin²(α-θ)]sin2(α-θ)，(8)

其中，α是相位角，指的是波数向量的方向与垂直方向的夹角，代表地震波的传播方向，θ指的是TTI介质对称轴与垂直方向的夹角，即倾角，TTI介质偏振角α-θ指的是波数向量的方向与TTI介质对称轴的夹角。

去除公式(8)的各向同性部分，只保留线性各向异性部分，可以得到：

Φ＝f[δ+2(ε-δ)sin²(α-θ)]sin2(α-θ)，(9)

如此一来，就可以假设初始模型为m＝{ε，δ，θ},在初始模型条件下偏振角的各向异性部分值为Φ。可以条件选择N个参考模型m^k，并计算得到相应的Φ^k，通过如下公式去插值得到Φ：

$\mathrm{\Φ}=\underset{k=1}{\overset{N}{\Σ}}{w}^k{\mathrm{\Φ}}^k,---\left(10\right)$

为了得到偏振向量的空间差异性，可使用弱各向异性下的偏振向量值近似各种强度各向异性下的偏振向量值。

假设Φ满足二阶平稳假设或者本征假设，运用变异函数来描述Φ随着空间位置不同而变化的特性，计算变异函数值：

$\mathrm{\μ}\left(\stackrel{-}{h}\right)=\frac{1}{2n}\underset{i-1}{\overset{n}{\Σ}}{\[\mathrm{\Φ}(\mathrm{\ϵ},\mathrm{\δ},\mathrm{\θ})-\mathrm{\Φ}(\mathrm{\ϵ}+h_{\mathrm{\ϵ}},\mathrm{\δ}+h_{\mathrm{\δ}},\mathrm{\θ}+h_{\mathrm{\θ}})\]}^2,---\left(11\right)$

基于相近相似的原则，使用空间位置增量替代各向异性参数增量，计算不同空间间隔下，对应的变异函数值，基于这个思想，可以将公式(11)改写为：

$\mathrm{\μ}\left(h\right)=\frac{1}{2n}\underset{i-1}{\overset{n}{\Σ}}{\[\mathrm{\Φ}(x,z)-\mathrm{\Φ}(x+h_x,z+h_x)\]}^2---\left(12\right)$

可以使用一些理论变异函数模型去拟合这些点，得到相关的参数。

较多采用的理论变异函数模型是球状模型以及指数模型，本发明实施例使用球状模型去拟合离散变异函数值，球状模型公式如下所示：

$\mathrm{\μ}\left(h\right)=\left\{\begin{array}{cc}{C}_0+C(\frac{3h}{2R}-\frac{1}{2}\frac{h^3}{R^3})& 0\≤h\≤R\\ C_0+C& h>R\end{array}\right.,---\left(13\right)$

在公式(13)中，C₀代表块金效应值，C₀+C为基台值，R为变程。

首先将公式(13)改写为：

y＝a₀+a₁x₁+a₂x₂,(14)在公式(14)中，y＝μ(h)，a₀＝C₀，x₁＝h，x₂＝-h³，且a₀≥0,a₁≥0,a₂≥0。

按照最小二乘拟合原理，误差最小，代入公式(14)，应满足：

$\underset{i=0}{\overset{m-1}{\Σ}}{(y_i-(a_0+a_1{x}_{1i}+a_2{x}_{2i}\left)\right)}^2=min,---\left(15\right)$

求出a₀,a₁,a₂后，就可按照球状模型的公式计算出C₀,C,R，表示如下：

C₀＝a₀,

$R=\sqrt{\frac{a_1}{3a_2}},$

$C=\frac{2a_1}{3}\sqrt{\frac{a_1}{3a_2}},$

将变异函数μ(h)引入公式(10)中的权重系数w^k，计算可以得到：

$w^k=\frac{1/\mathrm{\μ}({\left({\mathrm{\ϵ}}_k-\mathrm{\ϵ}\right)}^2+{\left({\mathrm{\δ}}_k-\mathrm{\δ}\right)}^2+{\left({\mathrm{\θ}}_k-\mathrm{\θ}\right)}^2)}{\underset{k=1}{\overset{n}{\Σ}}1/\mathrm{\μ}({\left({\mathrm{\ϵ}}_k-\mathrm{\ϵ}\right)}^2+{\left({\mathrm{\δ}}_k-\mathrm{\δ}\right)}^2+{\left({\mathrm{\θ}}_k-\mathrm{\θ}\right)}^2)},---\left(16\right)$

公式(16)可以计算得到权重系数。

将分离后的波场采用以下成像条件进行归一化互相关运算：

$\begin{array}{c}{I}_{pp}(x,z)=\frac{\underset{time}{\Σ}{S}_p\left(x,z,t\right)R_p\left(x,z,t\right)}{\underset{time}{\Σ}{S_p}^2\left(x,z,t\right)},\\ I_{ps}(x,z)=\frac{\underset{time}{\Σ}{S}_p\left(x,z,t\right)R_s\left(x,z,t\right)}{\underset{time}{\Σ}{S_p}^2\left(x,z,t\right)},\\ I_{sp}(x,z)=\frac{\underset{time}{\Σ}{S}_s\left(x,z,t\right)R_p\left(x,z,t\right)}{\underset{time}{\Σ}{S_s}^2\left(x,z,t\right)},\\ I_{ss}(x,z)=\frac{\underset{time}{\Σ}{S}_s\left(x,z,t\right)R_s\left(x,z,t\right)}{\underset{time}{\Σ}{S_s}^2\left(x,z,t\right)},\end{array}---\left(17\right)$

以二维各向同性介质弹性波波动方程离散差分格式为例说明有效边界存储策略。采用高阶常规网格进行时间二阶，空间N阶离散，公式如下：

$\begin{array}{c}{u}_x^{i,k,n+1}=2u_x^{i,k,n}-u_x^{i,k,n-1}+{\mathrm{\Δt}}^2\left[\begin{array}{c}\frac{1}{\mathrm{\Δ}x}\underset{m=1}{\overset{N/2}{\Σ}}{a}_m({\mathrm{\σ}}_{xx}^{i+m,k,n}-{\mathrm{\σ}}_{xx}^{i-m,k,n})+\mathrm{...}\\ \frac{1}{\mathrm{\Δ}z}\underset{m=1}{\overset{N/2}{\Σ}}{a}_m({\mathrm{\σ}}_{xz}^{i,k+m,n}-{\mathrm{\σ}}_{xz}^{i,k-m,n})\end{array}\right],\\ u_z^{i,k,n+1}=2u_z^{i,k,n}-u_z^{i,k,n-1}+{\mathrm{\Δt}}^2\left[\begin{array}{c}\frac{1}{\mathrm{\Δ}x}\underset{m=1}{\overset{N/2}{\Σ}}{a}_m({\mathrm{\σ}}_{zx}^{i+m,k,n}-{\mathrm{\σ}}_{zx}^{i-m,k,n})+\mathrm{...}\\ \frac{1}{\mathrm{\Δ}z}\underset{m=1}{\overset{N/2}{\Σ}}{a}_m({\mathrm{\σ}}_{zz}^{i,k+m,n}-{\mathrm{\σ}}_{zz}^{i,k-m,n})\end{array}\right],\end{array}---\left(18\right)$

$\left[\begin{array}{c}{\mathrm{\σ}}_{xx}\\ {\mathrm{\σ}}_{zz}\\ {\mathrm{\σ}}_{zx}\end{array}\right]=c\left[\begin{array}{c}{\mathrm{\ϵ}}_{xx}\\ {\mathrm{\ϵ}}_{zz}\\ {\mathrm{\ϵ}}_{zx}\end{array}\right]---\left(19\right)$

$\begin{array}{c}{\mathrm{\ϵ}}_{xx}=\frac{\∂u_x}{\∂x},\\ {\mathrm{\ϵ}}_{zz}=\frac{\∂u_z}{\∂z},\\ {\mathrm{\ϵ}}_{zx}=\frac{1}{2}\left(\frac{\∂u_z}{\∂x}+\frac{\∂u_x}{\∂z}\right),\end{array}---\left(20\right)$

其中，u_x和u_z是离散的时间域波场， $a_m=-\frac{1}{m}\cos \left(\mathrm{m\π}\right)w\left(m\right),m=\mathrm{1,2}...,N/2,$ 为时间域的优化系数，σ＝(σ_xx,σ_zz,σ_zx)为应力张量，ε＝(ε_xx,ε_zz,ε_zx)为应变张量，c为弹性系数张量，Δx,Δz为空间步长，i,k表示空间位置，n表示时间步。

公式(18)，(19)，(20)是正传模拟时的公式，对于反传模拟公式，只要将公式(18)中n+1时刻和n-1时刻对换位置：

$\begin{array}{c}{u}_x^{i,k,n-1}=2u_x^{i,k,n}-u_x^{i,k,n+1}+{\mathrm{\Δt}}^2\left[\begin{array}{c}\frac{1}{\mathrm{\Δ}x}\underset{m=1}{\overset{N/2}{\Σ}}{a}_m({\mathrm{\σ}}_{xx}^{i+m,k,n}-{\mathrm{\σ}}_{xx}^{i-m,k,n})+\mathrm{...}\\ \frac{1}{\mathrm{\Δ}z}\underset{m=1}{\overset{N/2}{\Σ}}{a}_m({\mathrm{\σ}}_{xz}^{i,k+m,n}-{\mathrm{\σ}}_{xz}^{i,k-m,n})\end{array}\right],\\ u_z^{i,k,n-1}=2u_z^{i,k,n}-u_z^{i,k,n+1}+{\mathrm{\Δt}}^2\left[\begin{array}{c}\frac{1}{\mathrm{\Δ}x}\underset{m=1}{\overset{N/2}{\Σ}}{a}_m({\mathrm{\σ}}_{zx}^{i+m,k,n}-{\mathrm{\σ}}_{zx}^{i-m,k,n})+\mathrm{...}\\ \frac{1}{\mathrm{\Δ}z}\underset{m=1}{\overset{N/2}{\Σ}}{a}_m({\mathrm{\σ}}_{zz}^{i,k+m,n}-{\mathrm{\σ}}_{zz}^{i,k-m,n})\end{array}\right],\end{array}---\left(21\right)$

以为例，说明有限边界存储策略。计算对x的一阶空间导数假设采用N阶空间差分算子，需要左右沿着X方向各N个点的值，如果是边界上的点，则需要边界区域和内部区域内至少N/2个点才能保证对x的一阶空间导数有相应的精度。逆时偏移仅对内部网格成像，只要内部网格的波场值正确，就有较好的成像结果。因此，在正传时只保存临近边界的内部区域里的N/2个值，每反传一个时间步长，就用保存的边界点替换边界，就可以重建反传波场。可以应用本构方程和几何方程，通过波场位移值和计算得到。因此，可以直接保存波场的离散位移值和只需要保存和上下左右临近边界的内部区域里的N个点，就可以重建反传波场。

以弹性波逆时偏移Marmousi2模型为例，模型大小为488×1000，横向采样间隔为6.245m，纵向采样间隔为6.245m。采用主频为15Hz的Ricker子波，纵波速度模型如图6所示，横波速度模型如图7所示，密度模型如图8所示，定义各向异性参数为ε＝0.33,δ＝-0.21,θ＝0。

采用分离步骤的单炮TI各向异性逆时偏移成像结果，应用弹性势的成像条件，运用窗函数和变异函数IDW插值算法将波场分为纵波势和横波势，使用有效边界条件存储部分边界，重构炮点分离波场，然后与检波点的分离波场进行归一化互相关运算，提取反射系数，采用爆炸震源，可高效率地得到如图7至图10所示的单炮成像结果。

图9至图12分别是本发明一实施例中PP波、PS波、SP波、SS波的单炮成像结果。如图9至图12所示，四种成像结果PP,PS,SP,SS是不一致的，因为纵波和横波有不同的照明角度，不同的组合有不同的结果，表明本发明实施例的逆时偏移成像方法具有较佳的成像效果。

弹性波数值模拟、分离以及波场存储严重制约着弹性波数值模拟的常规化应用，各向异性介质弹性波数值模拟涉及到海量的波场传播计算，且其是细粒度并行的算法，适合GPU多核并行的架构。

因此，在上述各实施例中，一个或多个步骤可以在图形处理器GPU中实现，例如，图1所示的各向异性介质中弹性波逆时偏移成像方法中，步骤S102、S103、S105、S106中的一个或多个步骤可以在图形处理器GPU中实现，以提高逆时偏移成像的效率。

本发明实施例的各向异性介质中弹性波逆时偏移成像方法，综合运用自褶积组合窗函数，GPU/CPU并行加速策略，提高弹性波数值模拟的精度和效率。针对弹性波逆时偏移中矢量波场计算量巨大的情况，通过引入变异函数IDW插值算法，在混合域分离矢量波场，极大地提高分离效率。运用有效边界存储策略，不存储炮点分离后的波场，只存储有效边界，在成像时重构炮点分离波场。

基于与图1所示的各向异性介质中弹性波逆时偏移成像方法相同的发明构思，本申请实施例还提供了一种各向异性介质中弹性波逆时偏移成像装置，如下面实施例所述。由于该各向异性介质中弹性波逆时偏移成像装置解决问题的原理与各向异性介质中弹性波逆时偏移成像方法相似，因此该各向异性介质中弹性波逆时偏移成像装置的实施可以参见各向异性介质中弹性波逆时偏移成像方法的实施，重复之处不再赘述。

图13是本发明实施例的各向异性介质中弹性波逆时偏移成像装置的结构示意图。如图13所示，本发明实施例的逆时偏移成像装置包括炮点正传模拟单元1210、炮点正传波场分离单元1220、炮点正传波场反传重构单元1230、检波点逆时反传模拟单元1240、检波点反传波场分离单元1250及成像单元1260。上述单元依次顺序连接。

炮点正传模拟单元1210用于进行炮点正传模拟，以构建地下炮点正传波场。

炮点正传波场分离单元1220用于利用变异函数IDW算法对所述炮点正传波场进行分离，得到分离后的所述炮点正传波场，并存储分离后的所述炮点正传波场的有效边界。

炮点正传波场反传重构单元1230用于以存储的所述有效边界作为边界条件，对分离后的所述炮点正传波场进行炮点逆时反传重构，得到炮点反传重构波场。

检波点逆时反传模拟单元1240用于获取地震数据中的矢量数据，并以所述矢量数据作为边界条件进行检波点逆时反传模拟，得到地下检波点反传波场。

检波点反传波场分离单元1250用于利用所述变异函数IDW算法对所述检波点反传波场进行分离，得到分离后的所述检波点反传波场。

成像单元1260用于利用归一化成像条件对所述炮点反传重构波场和分离后的所述检波点波场进行成像，得到地下成像数据，生成成像剖面并输出。

本发明实施例的各向异性介质中弹性波逆时偏移成像装置，炮点正传波场分离单元存储有效边界，可以减少边界存储量，炮点正传波场分离单元及检波点反传波场分离单元运用变异函数IDW算法插值运算分离波场可有效提高分离效率。

图14是本发明一实施例的各向异性介质中弹性波逆时偏移成像装置的结构示意图。如图14所示，如图13所示的各向异性介质中弹性波逆时偏移成像装置，还可包括滤波单元1270。

滤波单元1270用于对所述成像数据进行低通滤波以压制其中的低波数噪音并输出所述成像剖面，以提高成像品质。

图15是本发明一实施例中波场分离单元的结构示意图。如图15所示，上述的炮点正传波场分离单元1220，可包括第一空间域至波数域变换模块1221、第一分离算子生成模块1222、第一波场分离模块1223、第一波数域至空间域变换模块1224及第一加权插值处理模块1225。上述各模块依次顺序连接。

第一空间域至波数域变换模块1221用于利用傅里叶变换分别将各时刻的波场由空间域变换至波数域。

第一分离算子生成模块1222用于选取N个参考模型，根据所述参考模型计算分离算子，并在波数域利用自褶积组合窗函数对所述分离算子进行截断。

第一波场分离模块1223用于在波数域利用所述分离算子对所述参考模型下的所述波场进行分离，得到所述参考模型下的分离后的所述波场。

第一波数域至空间域变换模块1224用于将所述参考模型下的分离后的所述波场由波数域变换至空间域，N个所述参考模型对应得到N个所述参考模型下的分离后的所述波场，N为正整数。

第一加权插值处理模块1225用于在空间域，根据所述参考模型的权重系数对N个所述参考模型下的分离后的所述波场进行加权插值处理，得到最终分离后的所述波场。

图16是本发明又一实施例中波场分离单元的结构示意图。如图16所示，上述的检波点反传波场分离单元1250，可包括第二空间域至波数域变换模块1251、第二分离算子生成模块1252、第二波场分离模块1253、第二波数域至空间域变换模块1254及第二加权插值处理模块1255。上述各模块依次顺序连接。

第二空间域至波数域变换模块1251用于利用傅里叶变换分别将各时刻的波场由空间域变换至波数域。

第二分离算子生成模块1252用于选取N个参考模型，根据所述参考模型计算分离算子，并在波数域利用自褶积组合窗函数对所述分离算子进行截断。

第二波场分离模块1253用于在波数域利用所述分离算子对所述参考模型下的所述波场进行分离，得到所述参考模型下的分离后的所述波场。

第二波数域至空间域变换模块1254用于将所述参考模型下的分离后的所述波场由波数域变换至空间域，N个所述参考模型对应得到N个所述参考模型下的分离后的所述波场，N为正整数。

第二加权插值处理模块1255用于在空间域，根据所述参考模型的权重系数对N个所述参考模型下的分离后的所述波场进行加权插值处理，得到最终分离后的所述波场。

本发明实施例中，波场分离模块在波数域利用分离算子分离波场，加权插值处理模块在空间域运用变异函数IDW算法插值运算得到最终分离后的波场，可以有效降低波场分离的计算复杂度。分离算子生成模块在波数域采用自褶积组合窗函数截断分离算子，再分离波场，可以提高波场分离的效率，可以保证截断得到的分离算子的准确度更高。

一个实施例中，上述装置中，炮点正传波场分离单元1220可包括第一参考模型选取模块，第一参考模型选取模块与其中的第一分离算子生成模块1222连接；以及检波点反传波场分离单元1250可包括第二参考模型选取模块，上述第二参考模型选取模块与上述第二分离算子生成模块1252连接。

上述第一参考模型选取模块及第二参考模型选取模块，均可用于遍历待数值模拟初始模型的多个各向异性参数，并利用变异函数的临界值进行约束，选取在所述临界值范围内出现概率最大的参考点，根据所述参考点搜索得到所述参考模型。

本发明实施例中，参考模型选取模块选取参考模型时，不仅考虑到两点之间的距离也考虑了空间位置关系，使得进行波场分离时插值精度和效率都较高。

又一实施例中，上述装置中，炮点正传波场分离单元1220还可包括第一权重系数生成模，其与上述第一加权插值处理模块1225连接；检波点反传波场分离单元1250还可包括第二权重系数生成模，其与上述第二加权插值处理模块1255连接。

第一权重系数生成模块及第二权重系数生成模块，均可用于通过上述变异函数计算得到所述参考模型的权重系数。

本发明实施例中，权重系数生成模块利用变异函数计算得到各个参考模型相对于初始模型的权重系数，使得计算得到的权重系数较为合理，达到了在保证波场分离的分离结果准确度的基础上，减少了计算量的效果。

再一实施例中，所述装置中，炮点正传波场分离单元1220还可包括第一变异函数生成模块，其与上述第一权重系数生成模块连接；检波点反传波场分离单元1250还可包括第二变异函数生成模块，其与上述第二权重系数生成模块连接。

第一变异函数生成模块及第二变异函数生成模块，均可用于通过计算所述初始模型的各向异性参数分布的变异函数值，拟合得到所述初始模型的变异函数。

另一实施例中，上述炮点正传波场分离单元1220还可包括第一自褶积组合窗函数生成模块，其与上述第一分离算子生成模块1222连接；上述检波点反传波场分离单元1250还可包括第二自褶积组合窗函数生成模块，其与上述第二分离算子生成模块1252连接。

上述第一自褶积组合窗函数生成模块及第二自褶积组合窗函数生成模块，均可用于：选择主瓣和旁瓣性能高于设定阈值的窗函数作为原始窗函数；对原始窗函数做L次自褶积运算得到自褶积后的窗函数，其中，L为正整数；对自褶积后的窗函数与原始窗函数进行加权运算，得到所述自褶积组合窗函数。

本发明实施例中，自褶积组合窗函数生成模块使得原始窗函数的主瓣性能变差，旁瓣性能变好，然后根据当前的需要对原始窗函数和自褶积后的窗函数进行加权运算，从而得到一个满足要求的窗函数。

一个实施例中，上述第一自褶积组合窗函数生成模块及第二自褶积组合窗函数生成模块，还均可用于：在确定需要主瓣性能优先的情况下，原始窗函数的加权系数高于自褶积后的窗函数的加权系数；在确定需要旁瓣性能优先的情况下，自褶积后的函数的加权系数高于原始窗函数的加权系数。

在上述各实施例中，一个或多个模块可以在图形处理器GPU中实现，以提高逆时偏移成像的效率。

本发明提出了新的策略以提高数值模拟、波场分离的效率以及减小波场存储。本发明的策略采用CPU/GPU协同并行的加速方式，提高弹性波正演的计算效率，并采用自褶积组合窗函数截断优化的方法，提高差分算子逼近微分算子的精度，压制数值频散。本发明首次将变异函数IDW插值算法分离波场的方法引入弹性波逆时偏移技术中，运用变异函数和参考点搜索策略选取参考模型并计算权重，变异函数IDW插值算法在计算权重时不仅考虑到了插值点和参考点之间的距离，也考虑到了同一方向不同尺度的多层次的空间变异性，可以有更可靠的插值效果。针对弹性波逆时偏移存储量大的问题，采用有效随机边界理论，进一步减少边界存储量，反传重建炮点波场，充分利用GPU能进行快速数值模拟计算的优势。

本领域内的技术人员应明白，本发明的实施例可提供为方法、系统、或计算机程序产品。因此，本发明可采用完全硬件实施例、完全软件实施例、或结合软件和硬件方面的实施例的形式。而且，本发明可采用在一个或多个其中包含有计算机可用程序代码的计算机可用存储介质(包括但不限于磁盘存储器、CD-ROM、光学存储器等)上实施的计算机程序产品的形式。

本发明是参照根据本发明实施例的方法、设备(系统)、和计算机程序产品的流程图和/或方框图来描述的。应理解可由计算机程序指令实现流程图和/或方框图中的每一流程和/或方框、以及流程图和/或方框图中的流程和/或方框的结合。可提供这些计算机程序指令到通用计算机、专用计算机、嵌入式处理机或其他可编程数据处理设备的处理器以产生一个机器，使得通过计算机或其他可编程数据处理设备的处理器执行的指令产生用于实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能的装置。

这些计算机程序指令也可存储在能引导计算机或其他可编程数据处理设备以特定方式工作的计算机可读存储器中，使得存储在该计算机可读存储器中的指令产生包括指令装置的制造品，该指令装置实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能。

这些计算机程序指令也可装载到计算机或其他可编程数据处理设备上，使得在计算机或其他可编程设备上执行一系列操作步骤以产生计算机实现的处理，从而在计算机或其他可编程设备上执行的指令提供用于实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能的步骤。

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

Claims (13)

1.一种各向异性介质中弹性波逆时偏移成像方法，其特征在于，所述方法包括：

进行炮点正传模拟，以构建地下炮点正传波场；

利用变异函数IDW算法对所述炮点正传波场进行分离，得到分离后的所述炮点正传波场，并存储分离后的所述炮点正传波场的有效边界；

以存储的所述有效边界作为边界条件，对分离后的所述炮点正传波场进行炮点逆时反传重构，得到炮点反传重构波场；

获取地震数据中的矢量数据，并以所述矢量数据作为边界条件进行检波点逆时反传模拟，得到地下检波点反传波场；

利用所述变异函数IDW算法对所述检波点反传波场进行分离，得到分离后的所述检波点反传波场；

利用归一化成像条件对所述炮点反传重构波场和分离后的所述检波点反传波场进行成像，得到地下成像数据，生成成像剖面并输出。

2.如权利要求1所述的各向异性介质中弹性波逆时偏移成像方法，其特征在于，所述的方法还包括：

对所述成像数据进行低通滤波以压制其中的低波数噪音并输出所述成像剖面。

3.如权利要求1所述的各向异性介质中弹性波逆时偏移成像方法，其特征在于，利用变异函数IDW算法对所述炮点正传波场进行分离，包括：

利用傅里叶变换分别将各时刻的所述炮点正传波场由空间域变换至波数域；

选取N个参考模型，根据所述参考模型计算分离算子，并在波数域利用自褶积组合窗函数对所述分离算子进行截断；

在波数域利用所述分离算子对所述参考模型下的所述炮点正传波场进行分离，得到所述参考模型下的分离后的所述炮点正传波场；

将所述参考模型下的分离后的所述炮点正传波场由波数域变换至空间域，N个所述参考模型对应得到N个所述参考模型下的分离后的所述炮点正传波场，N为正整数；

在空间域，根据所述参考模型的权重系数对N个所述参考模型下的分离后的所述炮点正传波场进行加权插值处理，得到最终分离后的所述炮点正传波场。

4.如权利要求1所述的各向异性介质中弹性波逆时偏移成像方法，其特征在于，利用所述变异函数IDW算法对所述检波点反传波场进行分离，包括：

利用傅里叶变换分别将各时刻的所述检波点反传波场由空间域变换至波数域；

选取N个参考模型，根据所述参考模型计算分离算子，并在波数域利用自褶积组合窗函数对所述分离算子进行截断；

在波数域利用所述分离算子对所述参考模型下的所述检波点反传波场进行分离，得到所述参考模型下的分离后的所述检波点反传波场；

将所述参考模型下的分离后的所述检波点反传波场由波数域变换至空间域，N个所述参考模型对应得到N个所述参考模型下的分离后的所述检波点反传波场，N为正整数；

在空间域，根据所述参考模型的权重系数对N个所述参考模型下的分离后的所述检波点反传波场进行加权插值处理，得到最终分离后的所述检波点反传波场。

5.如权利要求3或4所述的各向异性介质中弹性波逆时偏移成像方法，其特征在于，选取所述参考模型的方法，包括：

遍历待数值模拟初始模型的多个各向异性参数，并利用变异函数的临界值进行约束，选取在所述临界值范围内出现概率最大的参考点，根据所述参考点搜索得到所述参考模型。

6.如权利要求5所述的各向异性介质中弹性波逆时偏移成像方法，其特征在于，所述参考模型的权重系数通过所述变异函数计算得到。

7.如权利要求6所述的各向异性介质中弹性波逆时偏移成像方法，其特征在于，所述变异函数的获取方法包括：

通过计算所述初始模型的各向异性参数分布的变异函数值，拟合得到所述初始模型的变异函数。

8.如权利要求3或4所述的各向异性介质中弹性波逆时偏移成像方法，其特征在于，获取所述自褶积组合窗函数的方法，包括：

选择主瓣和旁瓣性能高于设定阈值的窗函数作为原始窗函数；

对原始窗函数做L次自褶积运算得到自褶积后的窗函数，其中，L为正整数；

对自褶积后的窗函数与原始窗函数进行加权运算，得到所述自褶积组合窗函数。

9.如权利要求8所述的各向异性介质中弹性波逆时偏移成像方法，其特征在于，

在确定需要主瓣性能优先的情况下，原始窗函数的加权系数高于自褶积后的窗函数的加权系数；

在确定需要旁瓣性能优先的情况下，自褶积后的函数的加权系数高于原始窗函数的加权系数。

10.一种各向异性介质中弹性波逆时偏移成像装置，其特征在于，所述装置包括：

炮点正传模拟单元，用于进行炮点正传模拟，以构建地下炮点正传波场；

炮点正传波场分离单元，用于利用变异函数IDW算法对所述炮点正传波场进行分离，得到分离后的所述炮点正传波场，并存储分离后的所述炮点正传波场的有效边界；

炮点正传波场反传重构单元，用于以存储的所述有效边界作为边界条件，对分离后的所述炮点正传波场进行炮点逆时反传重构，得到炮点反传重构波场；

检波点逆时反传模拟单元，用于获取地震数据中的矢量数据，并以所述矢量数据作为边界条件进行检波点逆时反传模拟，得到地下检波点反传波场；

检波点反传波场分离单元，用于利用所述变异函数IDW算法对所述检波点反传波场进行分离，得到分离后的所述检波点反传波场；

成像单元，利用归一化成像条件对所述炮点反传重构波场和分离后的所述检波点波场进行成像，得到地下成像数据，生成成像剖面并输出。

11.如权利要求10所述的各向异性介质中弹性波逆时偏移成像装置，其特征在于，所述装置还包括：

滤波单元，用于对所述成像数据进行低通滤波以压制其中的低波数噪音并输出所述成像剖面。

12.如权利要求10所述的各向异性介质中弹性波逆时偏移成像装置，其特征在于，所述炮点正传波场分离单元包括：

第一空间域至波数域变换模块，用于利用傅里叶变换分别将各时刻的所述炮点正传波场由空间域变换至波数域；

第一分离算子生成模块，用于选取N个参考模型，根据所述参考模型计算分离算子，并在波数域利用自褶积组合窗函数对所述分离算子进行截断；

第一波场分离模块，用于在波数域利用所述分离算子对所述参考模型下的所述炮点正传波场进行分离，得到所述参考模型下的分离后的所述炮点正传波场；

第一波数域至空间域变换模块，用于将所述参考模型下的分离后的所述炮点正传波场由波数域变换至空间域，N个所述参考模型对应得到N个所述参考模型下的分离后的所述炮点正传波场，N为正整数；

第一加权插值处理模块，用于在空间域，根据所述参考模型的权重系数对N个所述参考模型下的分离后的所述炮点正传波场进行加权插值处理，得到最终分离后的所述炮点正传波场。

13.如权利要求10所述的各向异性介质中弹性波逆时偏移成像装置，其特征在于，检波点反传波场分离单元包括：

第二空间域至波数域变换模块，用于利用傅里叶变换分别将各时刻的所述检波点反传波场由空间域变换至波数域；

第二分离算子生成模块，用于选取N个参考模型，根据所述参考模型计算分离算子，并在波数域利用自褶积组合窗函数对所述分离算子进行截断；

第二波场分离模块，用于在波数域利用所述分离算子对所述参考模型下的所述检波点反传波场进行分离，得到所述参考模型下的分离后的所述检波点反传波场；

第二波数域至空间域变换模块，用于将所述参考模型下的分离后的所述检波点反传波场由波数域变换至空间域，N个所述参考模型对应得到N个所述参考模型下的分离后的所述检波点反传波场，N为正整数；

第二加权插值处理模块，用于在空间域，根据所述参考模型的权重系数对N个所述参考模型下的分离后的所述检波点反传波场进行加权插值处理，得到最终分离后的所述检波点反传波场。