CN104133241A - 波场分离方法和装置 - Google Patents

Abstract

本发明提供了一种波场分离方法和装置，该方法包括：计算初始模型各向异性参数分布的变异函数值，拟合得到该初始模型的变异函数；利用所述变异函数和参考点搜索策略，从初始模型中选取N个参考模型；将空间域的矢量弹性波场变换至波数域；对选取的N个参考模型中的每个参考模型执行以下操作：截断拟微分算子，在波数域进行波场分离，将纵横波波场反变换回空间域，得到该参考模型下的波场分离结果；利用变异函数计算各个参考模型相对于初始模型的权重系数；在空间域对各个参考模型下的波场分离结果进行加权插值处理得到初始模型的波场分离结果。本发明实施达到了在保证波场分离的分离结果准确度的基础上，减少了计算量的效果。

Description

波场分离方法和装置

技术领域

本发明涉及油气勘探技术领域，特别涉及一种波场分离方法和装置。

背景技术

随着油气勘探程度的加深，对更高精度的成像和反演的需求也越来越迫切。常规的纵波地震勘探虽然应用广泛，但是对于特定的复杂的地质结构的成像精度不够，因此，发展基于弹性波理论的多波多分量地震勘探技术有一定的现实意义。弹性矢量波场逆时偏移技术采用全波波动方程，会获得更加丰富的细节信息，不论在各向同性介质或者是各向异性介质中，弹性矢量波对地下介质的成像效果都要优于单纯的声波成像效果。但是弹性波逆时偏移面临的一个最大的问题是传统的互相关成像条件不分离弹性矢量波场，只是对波场的垂直分量和水平分量分别成像，然而对垂直分量和水平分量分别成像，在物理上毫无意义，因为纵横波混杂在一起会造成串扰和假象，这就需要在成像前，先对矢量波场进行分离。

目前，矢量波场分离主要有两种方法：

1)通过求解Christoffel方程，得到qP，qSV和qSH波的精确频散关系，其中，q是quasi-的缩写，因为在各向异性介质中，只有在特定方向，纵波和横波的偏振方向才是水平和垂直波传播方向的，严格上并不是真正的P波和S波，近似频散关系之后，再进行傅里叶反变换得到相应的波动方程。目前有一种方式是假设横波速度为零，推导出VTI介质qP波频散方程，再经过傅里叶反变换，从而得到VTI介质qP波频率-空间域的波动方程。

2)利用地震波场的偏振特性对矢量波场进行分离，在各向同性介质中，利用Helmholtz定理对矢量波场分别求取散度和旋度得到纵波和横波，其本质是将矢量波场分别向水平波数和垂直波数方向进行投影。根据在各向异性介质中，纵波和横波的偏振方向并不一定是水平或者垂直波数方向的特点，通过求解波数域Christoffel方程，得到纵横波的偏振向量，并分别将矢量波场分别投影到相应的偏振方向上，从而得到qP波和qSV波。后来，将该方法推广到三维各向异性介质中的纵横波分离，并利用拟微分算子，在空间域进行纵横波分离，将Dellinger的理论从均匀介质推广到非均匀介质，为了提高分离效率，还提出了在混合域分离的方法，即，根据模型参数选择一些参考模型，在波数域进行分离，然后在空间域插值得到该模型参数下的纵横波波场。主要的插值方式有两种，一种是构建一个联系各向异性参数和偏振向量的解析函数，计算插值系数，但该函数是基于弱各向异性理论建立的，对于强各向异性介质，存在偏差，将导致插值出现问题，另一种是使用空间插值的方法，采用了IDW(InverseDistance Weights)算法，IDW算法基于相似相近的原则，受非均匀分布的数据点影响比较大，数据点越多，插值越准确，但如果选择的数据点距离一致，其权重也一致，难以有效区分数据点的空间构性。

然而，为了达到很好的精确度，目前对各向异性的处理都是在空间域进行的，这样虽然可以保证精确度，但是计算量很大，计算效率比较低。

发明内容

本发明实施例提供了一种波场分离方法，以减少波场分离的计算量，该方法包括：

计算初始模型各向异性参数分布的变异函数值，拟合得到该初始模型的变异函数；

利用所述变异函数和参考点搜索策略，从初始模型中选取N个参考模型，其中，N为正整数；

将空间域的矢量弹性波场变换至波数域；

对选取的N个参考模型中的每个参考模型执行以下操作：根据参考模型计算得到拟微分算子，并采用自褶积窗函数截断拟微分算子；对该参考模型下的矢量弹性波场在波数域进行波场分离；将在波数域进行分离后的纵横波波场反变换回空间域，得到该参考模型下的波场分离结果；

利用所述变异函数计算各个参考模型相对于初始模型的权重系数；

根据计算得到的各个参考模型的权重系数，在空间域对各个参考模型下的波场分离结果进行加权插值处理；

将插值处理后得到的结果作为初始模型的波场分离结果。

在一个实施例中，所述自褶积窗函数是按照以下方式得到的：

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

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

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

在一个实施例中，按照以下公式计算初始模型各向异性参数分布的变异函数值：

$\mathrm{\μ}\left(h\right)=\frac{1}{2n}\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{[\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{\μ}({({\mathrm{\ϵ}}_k-\mathrm{\ϵ})}^2+{({\mathrm{\δ}}_k-\mathrm{\δ})}^2+{({\mathrm{\θ}}_k-\mathrm{\θ})}^2)}{\underset{k=1}{\overset{n}{\mathrm{\Σ}}}1/\mathrm{\μ}({({\mathrm{\ϵ}}_k-\mathrm{\ϵ})}^2+{({\mathrm{\δ}}_k-\mathrm{\δ})}^2+{({\mathrm{\θ}}_k-\mathrm{\θ})}^2)}$

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

在一个实施例中，按照以下方式选择插值参考点：

根据所述初始模型计算得到多个分散的变异函数值；

对得到的多个分散的变异函数值进行拟合，得到理论变异函数模型；

根据所述变异函数模型确定变异函数中的各个参数的参数值；

遍历初始模型中不同数值的ε、δ、θ出现的概率，选取在拟合的变异函数临界值之内的，出现概率最大的点作为插值参考点。

本发明实施例还提供了一种波场分离装置，以减少波场分离的计算量，该装置包括：

变异函数求取模块，用于计算初始模型各向异性参数分布的变异函数值，拟合得到该初始模型的变异函数；

参考模型选取模块，用于利用所述变异函数和参考点搜索策略，从初始模型中选取N个参考模型，其中，N为正整数；；

域变换模块，用于将空间域的矢量弹性波场变换至波数域；

第一波场分离模块，用于对选取的N个参考模型中的每个参考模型执行以下操作：根据参考模型计算得到拟微分算子，并采用自褶积窗函数截断拟微分算子；对该参考模型下的矢量弹性波场在波数域进行波场分离；将在波数域进行分离后的纵横波波场反变换回空间域，得到该参考模型下的波场分离结果；

权重系数计算模块，用于利用所述变异函数计算各个参考模型相对于初始模型的权重系数；

插值模块，用于根据计算得到的各个参考模型的权重系数，在空间域对各个参考模型下的波场分离结果进行加权插值处理；

第二波场分离模块，用于将插值处理后得到的结果作为初始模型的波场分离结果。

在一个实施例中，所述第一波场分离模块还用于按照以下方式计算自褶积窗函数：

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

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

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

在一个实施例中，所述变异函数求取模块具体用于按照以下公式计算初始模型各向异性参数分布的变异函数值：

$\mathrm{\μ}\left(h\right)=\frac{1}{2n}\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{[\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{\μ}({({\mathrm{\ϵ}}_k-\mathrm{\ϵ})}^2+{({\mathrm{\δ}}_k-\mathrm{\δ})}^2+{({\mathrm{\θ}}_k-\mathrm{\θ})}^2)}{\underset{k=1}{\overset{n}{\mathrm{\Σ}}}1/\mathrm{\μ}({({\mathrm{\ϵ}}_k-\mathrm{\ϵ})}^2+{({\mathrm{\δ}}_k-\mathrm{\δ})}^2+{({\mathrm{\θ}}_k-\mathrm{\θ})}^2)}$

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

在一个实施例中，所述权重系数计算模块按照以下方式选择插值参考点：

根据所述初始模型计算得到多个分散的变异函数值；

对得到的多个分散的变异函数值进行拟合，得到理论变异函数模型；

根据所述变异函数确定变异函数中的各个参数的参数值；

遍历初始模型中不同数值的ε、δ、θ出现的概率，选取在拟合的变异函数临界值之内的，出现概率最大的点作为插值参考点。

在本发明实施例中，通过将矢量弹性波场由空间域变换至波数域进行波场分离，有效降低了计算复杂度，同时在波数域进行波场分离的时候采用自褶积窗函数进行截断，可以保证截断得到的微分算子的准确度更高，进一步的，采用变异函数和参考点搜索策略，求取插值处理时各个参考模型相对于初始模型的权重系数，使得计算得到的权重系数较为合理，达到了在保证波场分离的分离结果准确度的基础上，减少了计算量的效果。

附图说明

此处所说明的附图用来提供对本发明的进一步理解，构成本申请的一部分，并不构成对本发明的限定。在附图中：

图1是本发明实施例的波场分离方法的流程图；

图2是本发明实施例的波场分离的详细流程图；

图3是本发明实施例的波场分离装置的结构框图。

具体实施方式

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

在本例中提出了一种波场分离方法，如图1所示，包括以下步骤：

步骤101：计算初始模型各向异性参数分布的变异函数值，拟合得到该初始模型的变异函数；

步骤102：利用所述变异函数和参考点搜索策略，从初始模型中选取N个参考模型，其中，N为正整数；

步骤103：将空间域的矢量弹性波场变换至波数域；

步骤104：对选取的N个参考模型中的每个参考模型执行以下操作：根据参考模型计算得到拟微分算子，并采用自褶积窗函数截断拟微分算子；对该参考模型下的矢量弹性波场在波数域进行波场分离；将在波数域进行分离后的纵横波波场反变换回空间域，得到该参考模型下的波场分离结果；

步骤105：利用所述变异函数计算各个参考模型相对于初始模型的权重系数；

步骤106：根据计算得到的各个参考模型的权重系数，在空间域对各个参考模型下的波场分离结果进行加权插值处理；

步骤107：将插值处理后得到的结果作为初始模型的波场分离结果。

在上述实施例中，通过将矢量弹性波场由空间域变换至波数域进行波场分离，有效降低了计算复杂度，同时在波数域进行波场分离的时候采用自褶积窗函数进行截断，可以保证截断得到的微分算子的准确度更高，进一步的，采用变异函数和参考点搜索策略，求取加权插值处理时各个参考模型相对于初始模型的权重系数，使得计算得到的权重系数较为合理，达到了在保证波场分离的分离结果准确度的基础上，减少了计算量的效果。

考虑到因为窗函数幅度响应的主瓣和旁瓣性能直接影响着差分逼近拟微分算子的精度问题，因此设计合适的窗函数的显得很重要。首先，需要明白主瓣和旁瓣的性能是如何影响逼近精度的，其次需要明白如何去控制主瓣和旁瓣，在本例中，提供了一种窗函数的获取方法，即，用自褶积窗函数，具体可以按照以下方式得到自褶积窗函数，包括：

步骤1：选择主瓣和旁瓣性能高于预定阈值的窗函数作为原始窗函数，例如，比较多种已有窗函数幅度响应的主瓣和旁瓣性能，选择最优的窗函数；

步骤2：因此窗函数自褶积会导致旁瓣衰减性能提升，相应的会引起旁瓣变宽，性能变弱，在本例中，对原始窗函数(即，选择的窗函数)做L次自褶积运算得到自褶积后的窗函数，即得到一个新的窗函数，其中，L为正整数；

步骤3：对自褶积后的窗函数与原始窗函数进行加权运算，得到所述自褶积窗函数，选择的加权系数的取值为0到1，可以通过对加权系数的选择确定主瓣性能优先还是旁瓣性能优先。

在按照上述方式获取自褶积窗函数之后就可以用该窗函数去截断得到优化的有限差分算子，还可以引入逼近误差函数，计算并绘制出逼近误差曲线，重点观察期频谱覆盖范围以及逼近精度的稳定性，如果效果不好，可以重新选取L的值或者加权系数的值来生成自褶积窗函数，直至得到满意的窗函数为止。

对于非均匀各向异性而言，理论上每个采样点都需要根据相应的各向异性参数计算偏振向量，才能获得更好的分离结果，但是这样处理计算量太大，因此，可以从初始模型中选择出N个参考模型，分别在波数域进行矢量波场分离，然后反变换回空间域，根据事先计算好的权重系数，重构初始模型的分离结果。因此，考虑到引入描述区域化变量的空间结构性变化和随机性变化的变异函数，代替两点之间的距离计算权重系数。具体的，可以按照以下公式计算初始模型各向异性参数分布的变异函数值：

$\mathrm{\μ}\left(h\right)=\frac{1}{2n}\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{[\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{\μ}({({\mathrm{\ϵ}}_k-\mathrm{\ϵ})}^2+{({\mathrm{\δ}}_k-\mathrm{\δ})}^2+{({\mathrm{\θ}}_k-\mathrm{\θ})}^2)}{\underset{k=1}{\overset{n}{\mathrm{\Σ}}}1/\mathrm{\μ}({({\mathrm{\ϵ}}_k-\mathrm{\ϵ})}^2+{({\mathrm{\δ}}_k-\mathrm{\δ})}^2+{({\mathrm{\θ}}_k-\mathrm{\θ})}^2)}$

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

上述表示位置增量，对于每一个都可以计算出一个以不同的的模为横坐标，为纵坐标，绘制点图，由绘制的点图可以发现：空间距离越小，变异性越小，越大，变异性越大，当增大到一个临界值R时，两个数值相互之间基本不会有太大的影响了，这也就表明在选取参考点，或者需要拟合的点的时候，应该在R范围内选取，否则选取的点没有太多的实际意义。例如：可以根据初始的参考模型计算得到多个分散的变异函数值，对得到的多个分散的变异函数值进行拟合得到理论变异函数模型，再根据拟合得到的变异函数模型确定变异函数中各个参数的参数值，遍历初始模型中不同数值的ε、δ、θ出现的概率，选取在拟合的变异函数临界值之内的，出现概率最大的点作为插值参考点。

下面结合一个具体的实施例进行说明，然而值得注意的是，该具体实施例仅是为了更好地说明本发明，并不构成对本发明的不当限定。

基于偏振特性对矢量波场进行分离的过程中所存在的一个不可避免的问题是如何在保证计算效率的前提下，让波场更准确地投影到相应的偏振方向上。在波数域，称之为投影，而在空间域，称之为空间滤波。在各向同性介质中，矢量波场的分离效果仅与微分算子的精度有关，精度越高，分离效果越好，同时计算量也增大，然而，在各向异性介质中，矢量波场的分离效果主要与各向异性参数和拟微分算子的精度这两个因素有关，考虑到在TI介质中，地震波动力学特征可以表示为各向同性部分与各向异性部分的叠加，可以通过Christoffel方程计算得到偏振向量，然后将其分解为波数向量与偏转向量之和，其中，偏转向量是以波数为自变量的函数，具有空间分布特性，在波数域利用常规二项式窗函数和Gauss窗函数截断逼近拟微分算子，并在空间域利用IDW算法去插值偏振向量的各向异性部分，以降低波场分离的计算量，从而提高波场分离的效率与精度。

具体而言，影响基于偏振特性分离矢量波场的技术的因素主要有两个：其一是截断窗函数的幅频响应特性，窗函数的幅频响应的主瓣形状控制着过渡带的范围，也就是频谱覆盖范围，旁瓣的形状决定了差分算子逼近微分算子的偏差程度，主瓣和旁瓣性能直接影响到了逼近的精度；其二是插值算法的插值效果，对于一种插值算法，既要求其有较高的插值精度，又需要较少的计算量，IDW算法执行效率较高，但是仅考虑到两点之间的距离计算权重，因此插值精度并不是很高，因此需要需求一种在保证执行效率的情况下，又可以提高插值精度的插值算法。

在本例中，提出了一种矢量波场分离方法，主要是为了解决一下两个问题：一是截断窗函数的选择，二是插值算法的选择。

首先，因为窗函数幅度响应的主瓣和旁瓣性能直接影响着差分逼近拟微分算子的精度问题，如果想设计合适的窗函数，首先要明白主瓣和旁瓣的性能如何影响逼近精度，其次是要研究如何去控制主瓣和旁瓣。然后，是插值算法的选择，将空间插值算法的优化引入矢量波场分离技术，因为对于非均匀各向异性介质，如果每一个点都计算一次拟微分算子，分离矢量波场有更好的效果，但是，这将耗费巨大的计算量，假设拟微分算子的大小是n²，模型的大小是N²，那么对模型进行矢量波场分离的计算量是2n²N²，这远远大于m阶精度的有限差分法2mN²的计算量，n<m。因此需要利用一种插值算法，在空间域通过插值的形式重构初始模型的分离结果。

在本例中，所依据的基本原理是：在各向同性介质中，应用Helmholtz定理，分别对波场求取散度和旋度，以分离纵横波，有如下公式：

$P=\&dtri;\&CenterDot;U,$

$S=\&dtri;\&times;U.$

在波数域可以表示为：

$P=\&dtri;\&CenterDot;U={\mathrm{ik}}_x{U}_x+{\mathrm{ik}}_z{U}_z;$

$S=\&dtri;\&times;U={\mathrm{ik}}_z{U}_x-{\mathrm{ik}}_x{U}_z.$

由上述公式可以看出，在各向同性介质中，P波是矢量波场在波数方向的投影，S波是矢量波场在垂直波数方向的投影，P波的偏振向量为(k_x,k_z)，S波的偏振向量与其正交，因而为(k_z,-k_x)。

对于各向异性介质，通过各向异性介质相对应的Christoffel方程，也可以得到各向异性介质中纵横波的偏振向量，以二维TTI介质为例，TTI介质对应的Christoffel方程为：

$\left[\begin{array}{cc}{\mathrm{\&Gamma;}}_{11}-\mathrm{\&rho;}{v}^2& {\mathrm{\&Gamma;}}_{12}\\ {\mathrm{\&Gamma;}}_{12}& {\mathrm{\&Gamma;}}_{22}-\mathrm{\&rho;}{v}^2\end{array}\right]\left[\begin{array}{c}{P}_x\\ P_z\end{array}\right]=0$

其中：

${\mathrm{\&Gamma;}}_{11}=c_{11}{k}_x^2+2c_{15}{k}_x{k}_z+c_{55}{k}_z^2$

${\mathrm{\&Gamma;}}_{12}=c_{15}{k}_x^2+(c_{13}+c_{55})k_x{k}_z+c_{35}{k}_z^2$

${\mathrm{\&Gamma;}}_{22}=c_{55}{k}_x^2+2c_{35}{k}_x{k}_z+c_{33}{k}_z^2$

其中，c₁₁,c₁₅,c₃₃,c₃₅,c₅₅表示弹性系数张量，k_x,k_z表示归一化波数，Γ表示Christoffel矩阵，上述公式是典型的特征值和特征向量的问题，为使公式有非零解，就要使系数行列式为零，因此可以求得偏振向量P＝(P_x,P_z)，这是k_x,k_z的函数，将iP反变换回空间域，便可以得到拟微分算子L，L也称作空间滤波算子。

现在引入表征TI介质(横向各向同性介质)各向异性的Thomsen系数V_P0,V_S0,ε,δ,γ，其中，V_P0,V_S0,ε,δ与qP波和qSV波有关，V_S0,γ与qSH波有关，因为qSH波是解耦的，因此矢量波场的分离仅需要V_P0,V_S0,ε,δ这四个参数。以二维TTI介质(具有倾斜对称轴的横向各向同性介质)为例，还需引入TI介质对称轴的倾角θ。在TI介质中，地震波动力学特征由两部分组成，第一个部分是各向同性的部分，第二个部分为各向异性的部分，可以将其近似表示为：

K≈K^iso+L(ε,δ)+Q(ε,δ,V_S0)

其中，K表示地震波的动力学特征，K^iso表示各向同性部分，ε＝δ＝0；L(ε,δ)+Q(ε,δ,V_S0)表示各向异性部分，L(ε,δ)表示线性部分，Q(ε,δ,V_S0)表示非线性部分。对于偏振向量，同样适用于上述近似公式，即，TI介质中的qP波和qSV波的偏振向量也由两部分组成：各向同性部分和各向异性部分。

对于各向同性部分，可以采用常规的有限差分算子的优化方法，即：使用窗函数截断方法或者最优化方法，这两种方法的本质都是类似的，都是希望在更高的波数范围达到一个较好的逼近精度。

对于各向异性部分，理论上每一个采样点都需要根据相应的各向异性参数计算偏振向量，以获得更好的分离效果，但这必然会带来海量的计算量，因此，考虑到采用插值的方式，即：从初始模型中选择N个参考模型，分别在波数域进行矢量波场分离，然后再反变换回空间域，再根据事先计算好的权重系数，重构初始模型的分离结果。既然是插值算法，插值的精度和插值算法的执行效率都是需要考虑的因素，IDW插值仅考虑到两点之间的距离计算权重，插值精度较差，且仅对均匀模型适用性较好，因此可以考虑引入描述区域化变量的空间结构性变化和随机性变化的变异函数，代替两点之间的距离，计算权重系数.

对上述内容具体阐述如下：

1、各向同性部分：

一个带限的连续信号f(x)在x＝0处的导数，可以被以一个均匀采样的信号f_n表示为：

$\frac{\&PartialD;f}{\&PartialD;x}{|}_{x=0}=\frac{1}{\mathrm{\&Delta;x}}\underset{n=-\&infin;}{\overset{\&infin;}{\mathrm{\&Sigma;}}}[-\frac{1}{n}\cos \left(\mathrm{n\&pi;}\right)]f\left(n\right)$

$\frac{{\&PartialD;}^{2}f}{\&PartialD;x^2}{|}_{x=0}=\frac{1}{\mathrm{\&Delta;}{x}^2}\underset{n=-\&infin;}{\overset{\&infin;}{\mathrm{\&Sigma;}}}[-\frac{2}{n^2}\cos \left(\mathrm{n\&pi;}\right)]f\left(n\right)$

存在一个长度为N+1点的窗函数，N为偶数，去截断上述两个公式，并经过一些简单的处理，可以得到有限差分算子：

$\frac{\&PartialD;f}{\&PartialD;x}{|}_{x=0}=\frac{1}{\mathrm{\&Delta;x}}\underset{n=1}{\overset{N/2}{\mathrm{\&Sigma;}}}{b}_n(f_n-f_{-n})$

$\frac{{\&PartialD;}^{2}f}{\&PartialD;x^2}{|}_{x=0}=\frac{1}{\mathrm{\&Delta;}{x}^2}(c_0{f}_0+\underset{n=1}{\overset{N/2}{\mathrm{\&Sigma;}}}{c}_n(f_n+f_{-n}))$

其中：

$b_n=d_n^{1}w\left(n\right)=-\frac{1}{n}\cos \left(\mathrm{n\&pi;}\right)w\left(n\right),n=\mathrm{1,2}...,N/2$

$c_n=d_n^{2}w\left(n\right)=-\frac{2}{n^2}\cos \left(\mathrm{n\&pi;}\right)w\left(n\right),n=\mathrm{1,2}...,N/2,C_0=-2\underset{n}{\overset{N/2}{\mathrm{\&Sigma;}}}(c_n+c_{-n}).$

为了优化有限差分算子逼近微分算子的精度，可以选择Chebyshev窗函数去截断：

$w_C\left(n\right)=\left\{\begin{array}{c}\frac{1}{N+1}\{\frac{1}{r}+2\underset{i=1}{\overset{N/2}{\mathrm{\&Sigma;}}}{C}_N\left(x_0\cos \left(\frac{\mathrm{i\&pi;}}{N+1}\right)\right)\cos \left(\frac{2\mathrm{ni\&pi;}}{N+1}\right)\},\left|n\right|\&le;N/2\\ 0,\left|n\right|>N/2\end{array}\right.$

$C_N\left(x\right)=\left\{\begin{array}{c}\cos \left[N{\cos }^{-1}\left(x\right)\right],\left|x\right|\&le;1\\ \cosh \left[N{\cosh }^{-1}\left(x\right)\right],\left|x\right|>1\end{array}\right.$

$x_0=\cosh \left[\frac{1}{N}{\cosh }^{-1}\left(\frac{1}{r}\right)\right]$

其中，r表示纹波率，代表旁瓣的衰减程度，N+1表示窗的长度，N为偶数。

不同窗函数幅频响应的主瓣和旁瓣性能影响着差分的逼近精度，具体影响方式包括：

1)主瓣大小和过渡带宽有关：主瓣窄，则过渡带窄，使用该窗函数截断逼近的有限差分算子的精度误差的谱覆盖范围大，可以用低阶算子达到高阶的精度，主瓣宽，则过渡带宽，使用该窗函数截断逼近的有限差分算子的精度误差的谱覆盖范围大。

2)旁瓣衰减和逼近精度稳定性的关系：旁瓣的衰减直接影响到了窗函数截断逼近的有限差分算子的精度误差的稳定性，旁瓣衰减越大，逼近精度误差波动越小，稳定性高，旁瓣衰减越小，逼近精度误差波动越大，稳定性低。

2、各向异性部分：

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

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

其中，α表示相角，传播方向与Z轴的夹角，代表的是各向同性部分，θ表示TTI介质对称轴的倾角， $f=1/2(1-V_{S0}^2/V_{P0}^2).$

去除上述近似公式中的各向同性部分，只保留线性各向异性部分，可以得到：

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

其中，Φ和ε、δ成近似线性关系，和sin2(α-θ)成近似比例关系。

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

$\mathrm{\&Phi;}=\underset{k=1}{\overset{N}{\mathrm{\&Sigma;}}}{w}^k{\mathrm{\&Phi;}}^k$

其中，w^k表示权重系数。

IDW插值算法在计算权重时，仅考虑插值点与参考点的空间位置关系，也就是当不同参考点与插值点空间位置一致时，就有相同的权重系数，这种插值算法虽然执行效率较高，但是没有考虑参考点之间的差异。

在本例中，对权重系数的确定进行了改进，具体包括：假设Φ满足二阶平稳假设或者本征假设，使用变异函数μ(h)来描述Φ随着空间位置不同而变化的特性，按照以下公式计算得到变异函数值：

$\mathrm{\&mu;}\left(h\right)=\frac{1}{2n}\underset{i-1}{\overset{n}{\mathrm{\&Sigma;}}}{[\mathrm{\&Phi;}(\mathrm{\&epsiv;},\mathrm{\&delta;},\mathrm{\&theta;})-\mathrm{\&Phi;}(\mathrm{\&epsiv;}+h_{\mathrm{\&epsiv;}},\mathrm{\&delta;}+h_{\mathrm{\&delta;}},\mathrm{\&theta;}+h_{\mathrm{\&theta;}})]}^2$

其中，表示位置增量，对于每一个都可以计算出一个以不同的的模为横坐标，为纵坐标，绘制点图，由绘制的点图可以发现：空间距离越小，变异性越小，越大，变异性越大，当增大到一个临界值R时，两个数值相互之间基本不会有太大的影响了，这也就说明在选取参考点，或者需要拟合的点的时候，应该在R范围内选取，否则选取的点没有太多的实际意义。这里描述得到的是零散的点，可以使用一些理论变异函数模型去拟合这些点，得到相关的参数。

引入变异函数R后可以得到权重系数w^k：

$w^k=\frac{1/\mathrm{\&mu;}({({\mathrm{\&epsiv;}}_k-\mathrm{\&epsiv;})}^2+{({\mathrm{\&delta;}}_k-\mathrm{\&delta;})}^2+{({\mathrm{\&theta;}}_k-\mathrm{\&theta;})}^2)}{\underset{k=1}{\overset{n}{\mathrm{\&Sigma;}}}1/\mathrm{\&mu;}({({\mathrm{\&epsiv;}}_k-\mathrm{\&epsiv;})}^2+{({\mathrm{\&delta;}}_k-\mathrm{\&delta;})}^2+{({\mathrm{\&theta;}}_k-\mathrm{\&theta;})}^2)}$

因为μ(h)是距离的函数，因此将μ(h)引入IDW插值算法的权重计算中，可以获得更为可靠的插值结果。

对于插值参考点(ε_k,δ_k,θ_k)的选择，可以遍历初始模型中不同(ε,δ,θ)出现的概率，并使用前面提到的变异函数临界值R进行约束，选取在R范围内的，出现概率最大的参考点。

在本例中，进行波场分离的总的流程图如图2所示，主要包括以下流程：

1)各向同性部分：

第一步：比较多种已有窗函数幅度响应的主瓣和旁瓣性能，选择最优的窗函数。

第二步：考虑到窗函数自褶积会导致旁瓣衰减性能提升，响应的会引起主瓣变宽，性能变弱，将第一步选择的最优窗函数做L次自褶积，产生一个新的窗函数。

第三步：将第二步产生的新的窗函数与第一步选择的最优窗函数进行加权组合，得到自褶积组合窗函数。选择加权系数为λ，λ取值为[0,1]，可以通过选择λ，确定是主瓣性能优先还是旁瓣性能优先。

第四步：应用第三步产生的窗函数去截断得到优化的有限差分算子，引入逼近误差函数，计算并绘出逼近误差曲线，重点观察其频谱覆盖范围以及逼近精度的稳定性，如果效果不好，返回第二步，重新开始循环，直到得到满足要求的窗函数。

第五步：得到优选的窗函数。

2)各向异性部分：

第一步：根据初始模型计算变异函数值。

第二步：得到分散的变异函数值后，拟合理论变异函数模型，得到变异函数的各个参数值。

第三步：遍历初始模型的三个参数(ε,δ,θ)，并以变异函数的距离参数临界值R进行约束，选取在R范围内的，出现概率最大的点作为插值参考点。

第四步：计算对应各个参考模型的权重系数w^k。

第五步：确定得到的权重系数是否满足要求，如果不满足，则重新执行第二步到第四步，重新选择理论变异函数模型和参考模型，计算权重系数。

3)综合部分：

第一步：进行弹性矢量波场的数值模拟，得到矢量弹性波场。

第二步：利用傅里叶变换，将空间域的矢量波场变换到波数域。

第三步：根据已选择的参考模型，计算拟微分算子，在波数域利用自褶积窗函数对算子进行截断。

第四步：在截断后进行矢量波场分离，将分离后的波场反变换回空间域。

第五步：重复第一步到第三步N次，计算N个参考模型下的波场分离结果。

第六步：利用计算的权重系数，进行插值，得到最终的波场分离结果。

在本例中，将变异函数理论、条件搜索参考点以及自褶积组合窗函数引入了矢量波场分离技术中，变异函数是距离的函数，描述了变量随着空间位置变化而变化的特性，将变异函数与IDW算法相结合，并应用于波场分离，是一个计算效率和插值精度的优选折中，通过构建变异函数，可以计算临界距离，应用于参考模型的选择策略，可以加强模型选择的约束，提高插值效果；设计的自褶积组合窗在保证旁瓣衰减的前提下，拥有较窄的主瓣，这样的幅度响应特性，可以保证逼近精度拥有较大的谱范围，且精度误差波动较小，稳定性较强。各向同性和各向异性两个方向优化矢量波场分离算法，可以在大幅提升计算效率的前提下，保持波场分离的精度。

基于同一发明构思，本发明实施例中还提供了一种波场分离装置，如下面的实施例所述。由于波场分离装置解决问题的原理与波场分离方法相似，因此波场分离装置的实施可以参见波场分离方法的实施，重复之处不再赘述。以下所使用的，术语“单元”或者“模块”可以实现预定功能的软件和/或硬件的组合。尽管以下实施例所描述的装置较佳地以软件来实现，但是硬件，或者软件和硬件的组合的实现也是可能并被构想的。图3是本发明实施例的波场分离装置的一种结构框图，如图3所示，包括：变异函数求取模块301、参考模型选取模块302、域变换模块303、第一波场分离模块304、权重系数计算模块305、、插值模块306、第二波场分离模块307，下面对该结构进行说明。

变异函数求取模块301，用于计算初始模型各向异性参数分布的变异函数值，拟合得到该模型的变异函数；

参考模型选取模块302，用于利用变异函数以及参考点搜索策略，从初始模型中选取N个参考模型，其中，N为正整数；

域变换模块303，用于将空间域的矢量弹性波场变换至波数域；

第一波场分离模块304，用于对选取的N个参考模型中的每个参考模型执行以下操作：根据参考模型计算得到拟微分算子，并采用自褶积窗函数截断拟微分算子；对该参考模型下的矢量弹性波场在波数域进行波场分离；将在波数域进行分离后的纵横波波场反变换回空间域，得到该参考模型下的波场分离结果；

权重系数计算模块305，用于利用变异函数计算各个参考模型相对于初始模型的权重系数；

插值模块306，用于根据计算得到的各个参考模型的权重系数，在空间域对各个参考模型下的波场分离结果进行加权插值处理；

第二波场分离模块307，用于将插值处理后得到的结果作为初始模型的波场分离结果。

在一个实施例中，第一波场分离模块304还用于按照以下方式计算自褶积窗函数：

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

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

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

在一个实施例中，变异函数求取模块301具体用于按照以下公式计算计算初始模型各向异性参数分布的变异函数值：

$\mathrm{\&mu;}\left(h\right)=\frac{1}{2n}\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{[\mathrm{\&Phi;}(\mathrm{\&epsiv;},\mathrm{\&delta;},\mathrm{\&theta;})-\mathrm{\&Phi;}(\mathrm{\&epsiv;}+h_{\mathrm{\&epsiv;}},\mathrm{\&delta;}+h_{\mathrm{\&delta;}},\mathrm{\&theta;}+h_{\mathrm{\&theta;}})]}^2$

其中，μ(h)表示变异函数值，Φ＝f[δ+2(ε-δ)sin²(α-θ)]sin2(α-θ)， ε、δ表示横向各向同性TI介质各向异性的系数，θ表示具有倾斜对称轴的横向各向同性TTI介质对称轴的倾角，表示插值参考点的位置增量，i＝1,2...n，其中，n表示插值参考点的个数，在计算得到多个分散的变异函数值后，根据这多个变异函数值进行拟合，得到理论变异函数模型。

在一个实施例中，权重系数计算模块305具体用于按照以下公式计算各个参考模型相对于初始模型的权重系数：

$w^k=\frac{1/\mathrm{\&mu;}({({\mathrm{\&epsiv;}}_k-\mathrm{\&epsiv;})}^2+{({\mathrm{\&delta;}}_k-\mathrm{\&delta;})}^2+{({\mathrm{\&theta;}}_k-\mathrm{\&theta;})}^2)}{\underset{k=1}{\overset{n}{\mathrm{\&Sigma;}}}1/\mathrm{\&mu;}({({\mathrm{\&epsiv;}}_k-\mathrm{\&epsiv;})}^2+{({\mathrm{\&delta;}}_k-\mathrm{\&delta;})}^2+{({\mathrm{\&theta;}}_k-\mathrm{\&theta;})}^2)}$

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

在一个实施例中，权重系数计算模块305按照以下方式选择插值参考点：

根据所述初始模型计算得到多个分散的变异函数值；

对得到的多个分散的变异函数值进行拟合，得到理论变异函数模型；

根据所述变异函数模型确定变异函数中的各个参数的参数值；

遍历初始模型中不同数值的ε、δ、θ出现的概率，选取在拟合的变异函数临界值之内的，出现概率最大的点作为插值参考点。

在另外一个实施例中，还提供了一种软件，该软件用于执行上述实施例及优选实施方式中描述的技术方案。

在另外一个实施例中，还提供了一种存储介质，该存储介质中存储有上述软件，该存储介质包括但不限于：光盘、软盘、硬盘、可擦写存储器等。

从以上的描述中，可以看出，本发明实施例实现了如下技术效果：通过将矢量弹性波场由空间域变换至波数域进行波场分离，有效降低了计算复杂度，同时在波数域进行波场分离的时候采用自褶积窗函数进行截断，可以保证截断得到的微分算子的准确度更高，进一步的，采用变异函数和参考点搜索策略，求取加权插值处理时各个参考模型相对于初始模型的权重系数，使得计算得到的权重系数较为合理，达到了在保证波场分离的分量结果准确度的基础上，减少了计算量的效果。

显然，本领域的技术人员应该明白，上述的本发明实施例的各模块或各步骤可以用通用的计算装置来实现，它们可以集中在单个的计算装置上，或者分布在多个计算装置所组成的网络上，可选地，它们可以用计算装置可执行的程序代码来实现，从而，可以将它们存储在存储装置中由计算装置来执行，并且在某些情况下，可以以不同于此处的顺序执行所示出或描述的步骤，或者将它们分别制作成各个集成电路模块，或者将它们中的多个模块或步骤制作成单个集成电路模块来实现。这样，本发明实施例不限制于任何特定的硬件和软件结合。

以上所述仅为本发明的优选实施例而已，并不用于限制本发明，对于本领域的技术人员来说，本发明实施例可以有各种更改和变化。凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。

Claims (10)

1.一种波场分离方法，其特征在于，包括：

计算初始模型各向异性参数分布的变异函数值，拟合得到该初始模型的变异函数；

利用所述变异函数和参考点搜索策略，从初始模型中选取N个参考模型，其中，N为正整数；

将空间域的矢量弹性波场变换至波数域；

对选取的N个参考模型中的每个参考模型执行以下操作：根据参考模型计算得到拟微分算子，并采用自褶积窗函数截断拟微分算子；对该参考模型下的矢量弹性波场在波数域进行波场分离；将在波数域进行分离后的纵横波波场反变换回空间域，得到该参考模型下的波场分离结果；

利用所述变异函数计算各个参考模型相对于初始模型的权重系数；

根据计算得到的各个参考模型的权重系数，在空间域对各个参考模型下的波场分离结果进行加权插值处理；

将插值处理后得到的结果作为初始模型的波场分离结果。

2.如权利要求1所述的方法，其特征在于，所述自褶积窗函数是按照以下方式得到的：

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

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

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

3.如权利要求1所述的方法，其特征在于，按照以下公式计算初始模型各向异性参数分布的变异函数值：

$\mathrm{\&mu;}\left(h\right)=\frac{1}{2n}\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{[\mathrm{\&Phi;}(\mathrm{\&epsiv;},\mathrm{\&delta;},\mathrm{\&theta;})-\mathrm{\&Phi;}(\mathrm{\&epsiv;}+h_{\mathrm{\&epsiv;}},\mathrm{\&delta;}+h_{\mathrm{\&delta;}},\mathrm{\&theta;}+h_{\mathrm{\&theta;}})]}^2$

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

4.如权利要求3所述的方法，其特征在于，按照以下公式计算各个参考模型相对于初始模型的权重系数：

$w^k=\frac{1/\mathrm{\&mu;}({({\mathrm{\&epsiv;}}_k-\mathrm{\&epsiv;})}^2+{({\mathrm{\&delta;}}_k-\mathrm{\&delta;})}^2+{({\mathrm{\&theta;}}_k-\mathrm{\&theta;})}^2)}{\underset{k=1}{\overset{n}{\mathrm{\&Sigma;}}}1/\mathrm{\&mu;}({({\mathrm{\&epsiv;}}_k-\mathrm{\&epsiv;})}^2+{({\mathrm{\&delta;}}_k-\mathrm{\&delta;})}^2+{({\mathrm{\&theta;}}_k-\mathrm{\&theta;})}^2)}$

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

5.如权利要求4所述的方法，其特征在于，按照以下方式选择插值参考点：

根据所述初始模型计算得到多个分散的变异函数值；

对得到的多个分散的变异函数值进行拟合，得到理论变异函数模型；

根据所述变异函数模型确定变异函数中的各个参数的参数值；

遍历初始模型中不同数值的ε、δ、θ出现的概率，选取在拟合的变异函数临界值之内的，出现概率最大的点作为插值参考点。

6.一种波场分离装置，其特征在于，包括：

变异函数求取模块，用于计算初始模型各向异性参数分布的变异函数值，拟合得到该初始模型的变异函数；

参考模型选取模块，用于利用所述变异函数和参考点搜索策略，从初始模型中选取N个参考模型，其中，N为正整数；

域变换模块，用于将空间域的矢量弹性波场变换至波数域；

第一波场分离模块，用于对选取的N个参考模型中的每个参考模型执行以下操作：根据参考模型计算得到拟微分算子，并采用自褶积窗函数截断拟微分算子；对该参考模型下的矢量弹性波场在波数域进行波场分离；将在波数域进行分离后的纵横波波场反变换回空间域，得到该参考模型下的波场分离结果；

权重系数计算模块，用于利用所述变异函数计算各个参考模型相对于初始模型的权重系数；

插值模块，用于根据计算得到的各个参考模型的权重系数，在空间域对各个参考模型下的波场分离结果进行加权插值处理；

第二波场分离模块，用于将插值处理后得到的结果作为初始模型的波场分离结果。

7.如权利要求6所述的装置，其特征在于，所述第一波场分离模块还用于按照以下方式计算自褶积窗函数：

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

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

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

8.如权利要求6所述的装置，其特征在于，所述变异函数求取模块具体用于按照以下公式计算初始模型各向异性参数分布的变异函数值；

$\mathrm{\&mu;}\left(h\right)=\frac{1}{2n}\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{[\mathrm{\&Phi;}(\mathrm{\&epsiv;},\mathrm{\&delta;},\mathrm{\&theta;})-\mathrm{\&Phi;}(\mathrm{\&epsiv;}+h_{\mathrm{\&epsiv;}},\mathrm{\&delta;}+h_{\mathrm{\&delta;}},\mathrm{\&theta;}+h_{\mathrm{\&theta;}})]}^2$

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

9.如权利要求8所述的装置，其特征在于，所述权重系数计算模块具体用于按照以下公式计算各个参考模型相对于初始模型的权重系数：

$w^k=\frac{1/\mathrm{\&mu;}({({\mathrm{\&epsiv;}}_k-\mathrm{\&epsiv;})}^2+{({\mathrm{\&delta;}}_k-\mathrm{\&delta;})}^2+{({\mathrm{\&theta;}}_k-\mathrm{\&theta;})}^2)}{\underset{k=1}{\overset{n}{\mathrm{\&Sigma;}}}1/\mathrm{\&mu;}({({\mathrm{\&epsiv;}}_k-\mathrm{\&epsiv;})}^2+{({\mathrm{\&delta;}}_k-\mathrm{\&delta;})}^2+{({\mathrm{\&theta;}}_k-\mathrm{\&theta;})}^2)}$

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

10.如权利要求9所述的装置，其特征在于，所述权重系数计算模块按照以下方式选择插值参考点：

根据所述初始模型计算得到多个分散的变异函数值；

对得到的多个分散的变异函数值进行拟合，得到理论变异函数模型；

根据所述变异函数模型确定变异函数中的各个参数的参数值；

遍历初始模型中不同数值的ε、δ、θ出现的概率，选取在拟合的变异函数临界值之内的，出现概率最大的点作为插值参考点。