CN104181598A - 计算地层不连续性属性值的方法及装置 - Google Patents

Abstract

本发明提供一种计算地层不连续性属性值的方法及装置，所述方法包括：读取三维地震数据体中各个点所对应的地震数据；根据读取的地震数据，获得梯度体结构张量方阵；从梯度体结构张量方阵的每个元素中抽取与三维地震数据体中的一个点的位置相对应的位置处的数据，并将抽取的数据按所对应的元素在梯度体结构张量方阵中的位置排列，以构成所述一个点的结构张量方阵；计算所述一个点的结构张量方阵的特征值；基于所述特征值，计算所述一个点的不连续性属性值。通过所述方法能够快速准确地计算出地层的不连续性属性值，从而能够清晰地反映断层和岩性的变化。

Description

计算地层不连续性属性值的方法及装置

技术领域

本发明属于石油地震勘探数据处理领域，更具体地说，涉及一种计算地层不连续性属性值的方法及装置。

背景技术

据统计，全球石油天然气产量有一半以上分布于裂缝性储层中。在我国，也广泛地发育有这样的裂缝性储层，如四川、华北、长庆、塔里木、克拉玛依、胜利、吉林、辽河、青海、玉门等许多地方都发现了裂缝性油气田。因此，准确地预测作为石油天然气通道的裂缝的发育情况对于理解地下地质储层的真实状态、降低勘探风险、提高钻探成功率有着重大的意义。地层的不连续性属性可反映断层和岩性的变化，被广泛应用于裂缝预测中。然而，目前应用于裂缝预测的方法中，具有计算复杂度高、速度慢、抗噪性能差、不能清晰反映断层和岩性变化等问题。

因此，需要一种计算地层不连续性属性值的方法及装置，以改善上述问题。

发明内容

本发明的目的在于提供一种地层不连续性属性值的计算方法及装置，这种方法计算速度快、抗噪性能强。

本发明的一方面提供一种计算地层不连续性属性值的方法，所述方法包括：a)读取三维地震数据体中各个点所对应的地震数据；b)根据读取的地震数据，获得梯度体结构张量方阵；c)从梯度体结构张量方阵的每个元素中抽取与三维地震数据体中的一个点的位置相对应的位置处的数据，并将抽取的数据按所对应的元素在梯度体结构张量方阵中的位置排列，以构成所述一个点的结构张量方阵；d)计算所述一个点的结构张量方阵的特征值；e)基于所述特征值，计算所述一个点的不连续性属性值。

可选地，所述一个点的不连续性属性值通过如下等式获得：

$c_n=\frac{\underset{i=0}{\overset{2}{\mathrm{\Σ}}}\underset{j=0}{\overset{2}{\mathrm{\Σ}}}{T}_{\mathrm{ns}}(i,j)}{3({\mathrm{\λ}}_{n1}+{\mathrm{\λ}}_{n2}+{\mathrm{\λ}}_{n3})}$

其中，c_n为所述一个点的不连续性属性值，T_ns(i,j)为所述一个点的结构张量方阵中位于第i行第j列的元素，λ_n1、λ_n2、λ_n3为所述一个点的结构张量方阵的三个特征值。

可选地，所述方法还包括：f)分别针对三维地震数据体中的其他点，重复步骤c)-e)来计算所述其他点的不连续性属性值。

可选地，步骤b)包括：b1)分别计算三维地震数据体中各个点所对应的地震数据在线号、道号、时间三个方向的梯度，以得到线号梯度体、道号梯度体、时间梯度体，并形成包括线号梯度体、道号梯度体、时间梯度体的梯度体向量；b2)基于所述梯度体向量，构建梯度体结构张量方阵。

可选地，步骤b)还包括：b3)对构建的梯度体结构张量方阵进行平滑。

可选地，在步骤b1)中，得到线号梯度体的步骤包括：计算三维地震数据体中每个点所对应的地震数据在线号方向的梯度，并将各个点对应的地震数据的梯度按对应点在三维地震数据体中的位置排列，以构成线号梯度体；得到道号梯度体的步骤包括：计算三维地震数据体中每个点所对应的地震数据在道号方向的梯度，并将各个点对应的地震数据的梯度按对应点在三维地震数据体中的位置排列，以构成道号梯度体；得到时间梯度体的步骤包括：计算三维地震数据体中每个点所对应的地震数据在时间方向的梯度，并将各个点对应的地震数据的梯度按对应点在三维地震数据体中的位置排列，以构成时间梯度体；形成梯度体向量的步骤包括：将线号梯度体、道号梯度体、时间梯度体按线号梯度体、道号梯度体、时间梯度体的顺序纵向排列，以构成梯度体向量。

可选地，计算三维地震数据体中每个点所对应的地震数据在线号方向的梯度的步骤包括：将与线号方向对应的卷积核与三维地震数据体中每个点所对应的地震数据沿线号方向进行卷积，以得到三维地震数据体中各个点所对应的地震数据在线号方向的梯度；计算三维地震数据体中每个点所对应的地震数据在道号方向的梯度的步骤包括：将与道号方向对应的卷积核与三维地震数据体中每个点所对应的地震数据沿道号方向进行卷积，以得到三维地震数据体中各个点所对应的地震数据在道号方向的梯度；计算三维地震数据体中每个点所对应的地震数据在时间方向的梯度的步骤包括：将与时间方向对应的卷积核与三维地震数据体中每个点所对应的地震数据沿时间方向进行卷积，以得到三维地震数据体中各个点所对应的地震数据在时间方向的梯度。

可选地，与线号方向对应的卷积核、与道号方向对应的卷积核、与时间方向对应的卷积核分别通过将一维零均值离散高斯核函数的导数在离散变量为相应的取值范围内的各个整数值时的函数值按对应的离散变量的从小到大的顺序排列构成，所述各卷积核中的函数值的计算式为：

$G^{\′}\left(t\right)=-\frac{1}{2\mathrm{\π}{{\mathrm{\σ}}_i}^2}t{e}^{[-t^2/\left(2{{\mathrm{\σ}}_i}^2\right)]}$

其中，G′(t)为一维零均值离散高斯核函数的导数，t为离散变量，t的取值范围为[-R_i，+R_i]，R_i ²＝4²σ_i，R_i为核半径，σ_i为预定尺度因子，i∈{x,y,z}，σ_x为与线号方向对应的预定尺度因子，R_x为与线号方向对应的核半径，σ_y为与道号方向对应的预定尺度因子，R_y为与道号方向对应的核半径，σ_z为与时间方向对应的预定尺度因子，R_z为与时间方向对应的核半径。

可选地，在步骤b2)中通过下面的等式构建梯度体结构张量方阵：

$T={\mathrm{gg}}^T=\left(\begin{array}{ccc}{g}_x{g}_x& g_x{g}_y& g_x{g}_z\\ g_y{g}_x& g_y{g}_y& g_y{g}_z\\ g_z{g}_x& g_z{g}_y& g_z{g}_z\end{array}\right)$

其中，T为梯度体结构张量方阵，g为梯度体向量， $g=\left(\begin{array}{c}{g}_x\\ g_y\\ g_z\end{array}\right),$ g_x为线号梯度体，g_y为道号梯度体，g_z为时间梯度体，g^T为梯度体向量的转置，梯度体结构张量方阵T的任意元素g_ug_v表示g_u中的每个位置的元素与g_v中的相同位置的元素相乘，u∈{x,y,z}，v∈{x,y,z}。

可选地，对构建的梯度体结构张量方阵进行平滑的步骤包括：使用具有预定尺度因子的三维零均值离散高斯核函数确定的卷积核与梯度体结构张量方阵进行卷积，以得到平滑后的梯度体结构张量方阵。

本发明另一方面提供一种计算地层不连续性属性值的装置，所述装置包括：读取数据单元，读取三维地震数据体中各个点所对应的地震数据；梯度体结构张量方阵获得单元，根据读取的地震数据，获得梯度体结构张量方阵；结构张量方阵构成单元，从梯度体结构张量方阵的每个元素中抽取与三维地震数据体中的一个点的位置相对应的位置处的数据，并将抽取的数据按所对应的元素在梯度体结构张量方阵中的位置排列，以构成所述一个点的结构张量方阵；特征值计算单元，计算所述一个点的结构张量方阵的特征值；不连续性属性值计算单元，基于所述特征值，计算所述一个点的不连续性属性值。

可选地，不连续性属性值计算单元通过如下等式获得所述一个点的不连续性属性值：

$c_n=\frac{\underset{i=0}{\overset{2}{\mathrm{\Σ}}}\underset{j=0}{\overset{2}{\mathrm{\Σ}}}{T}_{\mathrm{ns}}(i,j)}{3({\mathrm{\λ}}_{n1}+{\mathrm{\λ}}_{n2}+{\mathrm{\λ}}_{n3})}$

其中，c_n为所述一个点的不连续性属性值，T_ns(i,j)为结构张量方阵构成单元所构成的所述一个点的结构张量方阵中位于第i行第j列的元素，λ_n1、λ_n2、λ_n3为所述一个点的结构张量方阵的三个特征值。

可选地，结构张量方阵构成单元从梯度体结构张量方阵的每个元素中分别抽取与三维地震数据体中的其他各个点的位置相对应的位置处的数据，并分别将抽取的数据按所对应的元素在梯度体结构张量方阵中的位置排列，以分别构成其他各个点的结构张量方阵；特征值计算单元分别计算其他各个点的结构张量方阵的特征值；不连续性属性值计算单元分别基于其他各个点的结构张量方阵的特征值，计算其他各个点的不连续性属性值。

可选地，所述梯度体结构张量方阵获得单元包括：梯度体向量计算单元，分别计算三维地震数据体中各个点所对应的地震数据在线号、道号、时间三个方向的梯度，以得到线号梯度体、道号梯度体、时间梯度体，并形成包括线号梯度体、道号梯度体、时间梯度体的梯度体向量；梯度体结构张量方阵构建单元，基于所述梯度体向量，构建梯度体结构张量方阵。

可选地，所述梯度体结构张量方阵获得单元还包括：平滑单元，对梯度体结构张量方阵构建单元所构建的梯度体结构张量方阵进行平滑。

可选地，梯度体向量计算单元包括：线号梯度体计算单元，计算三维地震数据体中每个点所对应的地震数据在线号方向的梯度，并将各个点对应的地震数据的梯度按对应点在三维地震数据体中的位置排列，以构成线号梯度体；道号梯度体计算单元，计算三维地震数据体中每个点所对应的地震数据在道号方向的梯度，并将各个点对应的地震数据的梯度按对应点在三维地震数据体中的位置排列，以构成道号梯度体；时间梯度体计算单元，计算三维地震数据体中每个点所对应的地震数据在时间方向的梯度，并将各个点对应的地震数据的梯度按对应点在三维地震数据体中的位置排列，以构成时间梯度体；排列单元，将线号梯度体、道号梯度体、时间梯度体按线号梯度体、道号梯度体、时间梯度体的顺序纵向排列，以构成梯度体向量。

可选地，线号梯度体计算单元通过将与线号方向对应的卷积核与三维地震数据体中每个点所对应的地震数据沿线号方向进行卷积，以得到三维地震数据体中各个点所对应的地震数据在线号方向的梯度；道号梯度体计算单元通过将与道号方向对应的卷积核与三维地震数据体中每个点所对应的地震数据沿道号方向进行卷积，以得到三维地震数据体中各个点所对应的地震数据在道号方向的梯度；时间梯度体计算单元通过将与时间方向对应的卷积核与三维地震数据体中每个点所对应的地震数据沿时间方向进行卷积，以得到三维地震数据体中各个点所对应的地震数据在时间方向的梯度。

可选地，与线号方向对应的卷积核、与道号方向对应的卷积核、与时间方向对应的卷积核分别通过将一维零均值离散高斯核函数的导数在离散变量为相应的取值范围内的各个整数值时的函数值按对应的离散变量的从小到大的顺序排列构成，所述各卷积核中的函数值的计算式为：

$G^{\′}\left(t\right)=-\frac{1}{2\mathrm{\π}{{\mathrm{\σ}}_i}^2}t{e}^{[-t^2/\left(2{{\mathrm{\σ}}_i}^2\right)]}$

其中，G′(t)为一维零均值离散高斯核函数的导数，t为离散变量，t的取值范围为[-R_i，+R_i]，R_i ²＝4²σ_i，R_i为核半径，σ_i为预定尺度因子，i∈{x,y,z}，σ_x为与线号方向对应的预定尺度因子，R_x为与线号方向对应的核半径，σ_y为与道号方向对应的预定尺度因子，R_y为与道号方向对应的核半径，σ_z为与时间方向对应的预定尺度因子，R_z为与时间方向对应的核半径。

可选地，梯度体结构张量方阵构建单元通过下面的等式构建梯度体结构张量方阵：

$T={\mathrm{gg}}^T=\left(\begin{array}{ccc}{g}_x{g}_x& g_x{g}_y& g_x{g}_z\\ g_y{g}_x& g_y{g}_y& g_y{g}_z\\ g_z{g}_x& g_z{g}_y& g_z{g}_z\end{array}\right)$

其中，T为梯度体结构张量方阵，g为梯度体向量， $g=\left(\begin{array}{c}{g}_x\\ g_y\\ g_z\end{array}\right),$ g_x为线号梯度体，g_y为道号梯度体，g_z为时间梯度体，g^T为梯度体向量的转置，梯度体结构张量方阵T的任意元素g_ug_v表示g_u中的每个位置的元素与g_v中的相同位置的元素相乘，u∈{x,y,z}，v∈{x,y,z}。

可选地，平滑单元使用具有预定尺度因子的三维零均值离散高斯核函数确定的卷积核与梯度体结构张量方阵进行卷积，以得到平滑后的梯度体结构张量方阵。

根据本发明能够提供一种计算速度快、抗噪性能强的地层不连续性属性值的计算方法及装置，从而通过计算出的地层不连续性属性值能够清晰地反映断层和岩性的变化。

将在接下来的描述中部分阐述本发明另外的方面和/或优点，还有一部分通过描述将是清楚的，或者可以经过本发明的实施而得知。

附图说明

通过下面结合附图进行的对实施例的描述，本发明的上述和/或其它目的和优点将会变得更加清楚，其中：

图1是示出根据本发明示例性实施例的计算地层不连续性属性值的方法的流程图；

图2是示出根据本发明示例性实施例的图1中的根据读取的地震数据获得梯度体结构张量方阵的方法的流程图；

图3是示出根据本发明示例性实施例的图2中的根据读取的地震数据获得包括线号梯度体、道号梯度体、时间梯度体的梯度体向量的方法的流程图；

图4示出采用本发明示例性实施例的计算地层不连续性属性值的方法的计算效率与传统方法的计算效率对比曲线图；

图5A示出未采用本发明示例性实施例的计算地层不连续性属性值的方法的不连续性属性效果图；

图5B示出采用本发明示例性实施例的计算地层不连续性属性值的方法的不连续性属性效果图；

图6是示出根据本发明示例性实施例的计算地层不连续性属性值的装置的方框图；

图7是示出根据本发明示例性实施例的图6中的梯度体结构张量方阵获得单元的方框图；

图8是示出根据本发明示例性实施例的图7中的梯度体向量计算单元的方框图；

图9示出三维地震数据体的空间模型。

具体实施方式

现将详细描述本发明的示例性实施例，所述实施例的示例在附图中示出，其中，相同的标号指示相同的部分。以下将通过参照附图来说明所述实施例，以便解释本发明。

图1是示出根据本发明示例性实施例的计算地层不连续性属性值的方法的流程图。

如图1所示，在步骤101，读取三维地震数据体中各个点所对应的地震数据。图9示出三维地震数据体的空间模型。在三维地震勘探中，得到的三维地震数据体的空间模型如图9所示，图9中的黑点代表地震数据，可将地震数据表示为D(x,y,z)。这里的x表示线号，y表示道号，z表示时间。地震数据D可理解为由线号(Inline)、道号(Xline)、时间(Time)所确定的点(即，采样点)所对应的数据。在这里，地震数据可采用叠后地震数据，具体地可以为振幅、瞬时相位等。可以理解，x、y、z为取整数的离散变量。

在步骤102，根据步骤101读取的地震数据，获得梯度体结构张量方阵。这里，获得梯度体结构张量方阵的方法将在后面结合图2和图3进行说明。

在步骤103，从步骤102得到的梯度体结构张量方阵的每个元素中抽取与三维地震数据体中的一个点的位置相对应的位置处的数据，并将抽取的数据按所对应的元素在梯度体结构张量方阵中的位置排列，以构成所述一个点的结构张量方阵。这里，所说的一个点可以是三维地震数据体中的任意点。根据后述内容可知梯度体结构张量方阵为3×3的方阵，所以所述一个点的结构张量方阵也为3×3的方阵。

在步骤104，计算步骤103的所述一个点的结构张量方阵的特征值。这里，可通过现有的各种方法来计算结构张量方阵的特征值。由于结构张量方阵为3×3的方阵，所以将得到三个特征值。

在步骤105，基于步骤104得到的特征值，计算所述一个点的不连续性属性值。这里，所述一个点的不连续性属性值可通过下面的等式(1)获得：

$c_n=\frac{\underset{i=0}{\overset{2}{\mathrm{\Σ}}}\underset{j=0}{\overset{2}{\mathrm{\Σ}}}{T}_{\mathrm{ns}}(i,j)}{3({\mathrm{\λ}}_{n1}+{\mathrm{\λ}}_{n2}+{\mathrm{\λ}}_{n3})}---\left(1\right)$

这里，c_n为所述一个点的不连续性属性值，T_ns(i,j)为所述一个点的结构张量方阵中位于第i行第j列的元素，λ_n1、λ_n2、λ_n3为所述一个点的结构张量方阵的三个特征值。

在图1所示的根据本发明示例性实施例的计算地层不连续性属性值的方法中还可包括分别针对三维地震数据体中的其他点，重复步骤103-105来计算所述其他点的不连续性属性值(即，将步骤103-105中的所述一个点替换为所述其他点中的每个点)，由此可得到三维地震数据体中所有点的地震数据的不连续性属性值。

下面结合图2和图3对前述步骤102中获得梯度体结构张量方阵的方法进行详细说明。

图2是示出根据本发明示例性实施例的图1中的根据读取的地震数据获得梯度体结构张量方阵的方法的流程图。如图2所示，在步骤201，利用地震数据得到梯度体向量。具体地说，分别计算读取的三维地震数据体中各个点所对应的地震数据在线号、道号、时间三个方向的梯度，以得到线号梯度体、道号梯度体、时间梯度体，并形成包括线号梯度体、道号梯度体、时间梯度体的梯度体向量。

在本发明的一个实施例中，可根据图3中所示的方法来计算并获得梯度体向量。

图3是示出根据本发明示例性实施例的获得梯度体向量的方法的流程图。

如图3所示，在步骤301，计算三维地震数据体中每个点所对应的地震数据在线号方向的梯度，并将各个点对应的地震数据的梯度按对应点在三维地震数据体中的位置排列，以构成线号梯度体。这里，可通过各种梯度计算方法来计算三维地震数据体中各个点所对应的地震数据在线号方向的梯度，从而得到线号梯度体。

优选地，利用本发明提出的方式来计算所述各个点的地震数据在线号方向的梯度。具体地说，将与线号方向对应的卷积核与三维地震数据体中每个点所对应的地震数据沿线号方向进行卷积，以得到三维地震数据体中各个点所对应的地震数据在线号方向的梯度，计算式如下式(2)：

h_x(x,y,z)＝f_x*D(x,y,z)                (2)

这里，h_x(x,y,z)为三维地震数据体中的任意点的地震数据在线号方向的梯度，f_x为与线号方向对应的卷积核，*为卷积符号，D(x,y,z)为三维地震数据体中由x、y、z的值确定的任意点的地震数据。

这里，与线号方向对应的卷积核f_x通过将一维零均值离散高斯核函数的导数在离散变量为与线号方向对应的取值范围内的各个整数值时的函数值按对应的离散变量的从小到大的顺序排列构成。具体地，对一维零均值离散高斯核函数求导数，得到等式(3)：

$G^{\′}\left(t\right)=-\frac{1}{2\mathrm{\π}{{\mathrm{\σ}}_i}^2}t{e}^{[-t^2/\left(2{{\mathrm{\σ}}_i}^2\right)]}---\left(3\right)$

这里，G′(t)为一维零均值离散高斯核函数的导数，t为离散变量，t的取值范围为[-R_x，+R_x]，R_x ²＝4²σ_x，R_x为与线号方向对应的核半径，σ_x为与线号方向对应的预定尺度因子。

接下来，将一维零均值离散高斯核函数的导数G′(t)在离散变量t为与线号方向对应的取值范围[-R_x，+R_x]内的各个整数值时的函数值按对应的离散变量t的从小到大的顺序排列得到卷积核。例如，假设与线号方向对应的预定尺度因子σ_x为1，则相应的核半径R_x等于4。那么，离散变量t的取值范围为[-4，+4]，t的从小到大的取值为{-4，-3，-2，-1，0，+1，+2，+3，+4}。将离散变量t的各个取值代入到等式(3)中则得到相应的一维零均值离散高斯核函数的导数的函数值分别为{0.0002，0.0053，0.043，0.0965，0，-0.0965，-0.043，-0.0053，-0.0002}。最后，按对应的离散变量t的从小到大的顺序将所述函数值排列得到与线号方向对应的卷积核f_x为(0.0002，0.0053，0.043，0.0965，0，-0.0965，-0.043，-0.0053，-0.0002)。

优选地，σ_x可以设为道号相邻的两个道之间的间距的三倍。

具体地，等式(2)可进一步表示为等式(4)：

$\begin{array}{c}{h}_x(x,y,z)\\ =\underset{n=1}{\overset{N}{\mathrm{\Σ}}}{f}_x\left(n\right)D(x+(\frac{2n-1}{2}-\frac{N}{2}),y,z)\end{array}---\left(4\right)$

这里，f_x(n)表示与线号方向对应的卷积核f_x中的第n个数值，N为所述卷积核f_x中包含的数值的个数，即离散变量t在取值范围[-R_x，+R_x]内可以取得的整数值的个数，可以理解，N为奇数。

在计算三维地震数据体中线号(即，x)小于x_min+(N-1)/2的点处的地震数据在线号方向的梯度时(这里，x_min为三维地震数据体中的地震数据的线号的最小值)，根据式(4)可知计算所述点处的地震数据在线号方向的梯度需要使用三维地震数据体外的点处的地震数据(例如，在计算D(x_min,y,z)在线号方向的梯度时，需要使用的(x_min-(N-1)/2,y,z)～(x_min-1,y,z)的地震数据)，由于所述点处无地震数据，所以将所述点处的地震数据的值均取为D(x_min,y,z)。

同时，在计算三维地震数据体中线号大于x_max-(N-1)/2的点处的地震数据在线号方向的梯度时(这里，x_max为三维地震数据体中的地震数据的线号的最大值)，根据式(4)可知计算所述点处的地震数据在线号方向的梯度需要使用三维地震数据体外的点处的地震数据(例如，在计算D(x_max,y,z)在线号方向的梯度时，需要使用的(x_max+1,y,z)～(x_max+(N-1)/2,y,z)的地震数据)，由于所述点处无地震数据，所以将所述点处的地震数据的值均取为D(x_max,y,z)。

可以理解，对三维地震数据体中线号小于x_min+(N-1)/2或大于x_max-(N-1)/2的点处的地震数据计算梯度时所使用到的三维地震数据体外的点处的地震数据可根据需要进行取值，例如，使所述点处的地震数据的取值均为0。

在步骤302，计算三维地震数据体中每个点所对应的地震数据在道号方向的梯度，并将各个点对应的地震数据的梯度按对应点在三维地震数据体中的位置排列，以构成道号梯度体。这里，可通过各种梯度计算方法来计算三维地震数据体中各个点所对应的地震数据在道号方向的梯度，从而得到道号梯度体。

优选地，利用本发明提出的方式来计算所述各个点的地震数据在道号方向的梯度。具体地说，将与道号方向对应的卷积核与三维地震数据体中每个点所对应的地震数据沿道号方向进行卷积，以得到三维地震数据体中各个点所对应的地震数据在道号方向的梯度，计算式如下式(5)：

h_y(x,y,z)＝f_y*D(x,y,z)                   (5)

这里，h_y(x,y,z)为三维地震数据体中的任意点的地震数据在道号方向的梯度，f_y为与道号方向对应的卷积核。

这里，与道号方向对应的卷积核f_y通过将一维零均值离散高斯核函数的导数在离散变量为与道号方向对应的取值范围内的各个整数值时的函数值按对应的离散变量的从小到大的顺序排列构成，所述卷积核的获得方法与步骤301中的与线号方向对应的卷积核的获得方法相同，区别在于，与线号方向对应的预定尺度因子σ_x变更为与道号方向对应的预定尺度因子σ_y，相应的，与线号方向对应的核半径R_x变更为与道号方向对应的核半径R_y。

优选地，σ_y可以设为线号相邻的两条线之间的间距的三倍。

具体地，等式(5)可进一步表示为等式(6)：

$\begin{array}{c}{h}_y(x,y,z)\\ =\underset{m=1}{\overset{M}{\mathrm{\Σ}}}{f}_y\left(m\right)D(x,y+(\frac{2m-1}{2}-\frac{M}{2}),z)\end{array}---\left(6\right)$

这里，f_y(m)表示与道号方向对应的卷积核f_y中的第m个数值，M为所述卷积核f_y中包含的数值的个数，即离散变量t在取值范围[-R_y，+R_y]内可以取得的整数值的个数，可以理解，M为奇数。

在计算三维地震数据体中道号(即，y)小于y_min+(M-1)/2的点处的地震数据在道号方向的梯度时(这里，y_min为三维地震数据体中的地震数据的道号的最小值)，根据式(6)可知计算所述点处的地震数据在道号方向的梯度需要使用三维地震数据体外的点处的地震数据(例如，在计算D(x,y_min,z)在道号方向的梯度时，需要使用的(x,y_min-(M-1)/2,z)～(x,y_min-1,z)的地震数据)，由于所述点处无地震数据，所以将所述点处的地震数据的值均取为D(x,y_min,z)。

同时，在计算三维地震数据体中道号大于y_max-(M-1)/2的点处的地震数据在道号方向的梯度时(这里，y_max为三维地震数据体中的地震数据的道号的最大值)，根据式(6)可知计算所述点处的地震数据在道号方向的梯度需要使用三维地震数据体外的点处的地震数据(例如，在计算D(x,y_max,z)在道号方向的梯度时，需要使用的(x,y_max+1,z)～(x,y_max+(M-1)/2,z)的地震数据)，由于所述点处无地震数据，所以将所述点处的地震数据的值均取为D(x,y_max,z)。

可以理解，对三维地震数据体中道号小于y_min+(M-1)/2或大于y_max-(M-1)/2的点处的地震数据计算梯度时所使用到的三维地震数据体外的点处的地震数据可根据需要进行取值，例如，使所述点处的地震数据的取值均为0。

在步骤303，计算三维地震数据体中每个点所对应的地震数据在时间方向的梯度，并将各个点对应的地震数据的梯度按对应点在三维地震数据体中的位置排列，以构成时间梯度体。这里，可通过各种梯度计算方法来计算三维地震数据体中各个点所对应的地震数据在时间方向的梯度，从而得到时间梯度体。

优选地，利用本发明提出的方式来计算所述各个点的地震数据在时间方向的梯度。具体地说，将与时间方向对应的卷积核与三维地震数据体中每个点所对应的地震数据沿时间方向进行卷积，以得到三维地震数据体中各个点所对应的地震数据在时间方向的梯度，计算式如下式(7)：

h_z(x,y,z)＝f_z*D(x,y,z)                           (7)

这里，h_z(x,y,z)为三维地震数据体中的任意点的地震数据在时间方向的梯度，f_z为与时间方向对应的卷积核。

这里，与时间方向对应的卷积核f_z通过将一维零均值离散高斯核函数的导数在离散变量为与时间方向对应的取值范围内的各个整数值时的函数值按对应的离散变量的从小到大的顺序排列构成，所述卷积核的获得方法与步骤301中的与线号方向对应的卷积核的获得方法相同，区别在于，与线号方向对应的预定尺度因子σ_x变更为与时间方向对应的预定尺度因子σ_z，相应的，与线号方向对应的核半径R_x变更为与道号方向对应的核半径R_z。

优选地，σ_z可以设为时间轴(即，z轴)采样间隔的三倍。

具体地，等式(7)可进一步表示为等式(8)：

$\begin{array}{c}{h}_z(x,y,z)\\ =\underset{r=1}{\overset{R}{\mathrm{\Σ}}}{f}_z\left(r\right)D(x,y,z+(\frac{2r-1}{2}-\frac{R}{2}))\end{array}---\left(8\right)$

这里，f_z(r)表示与时间方向对应的卷积核f_z中的第r个数值，R为所述卷积核f_z中包含的数值的个数，即离散变量t在取值范围[-R_z，+R_z]内可以取得的整数值的个数，可以理解，R为奇数。

在计算三维地震数据体中时间(即，z)小于z_min+(R-1)/2的点处的地震数据在时间方向的梯度时(这里，z_min为三维地震数据体中的地震数据的时间的最小值)，根据式(8)可知计算所述点处的地震数据在时间方向的梯度需要使用三维地震数据体外的点处的地震数据(例如，在计算D(x,y,z_min)在时间方向的梯度时，需要使用的(x,y,z_min-(R-1)/2)～(x,y,z_min-1)的地震数据)，由于所述点处无地震数据，所以将所述点处的地震数据的值均取为D(x,y,z_min)。

同时，在计算三维地震数据体中时间大于z_max-(R-1)/2的点处的地震数据在时间方向的梯度时(这里，z_max为三维地震数据体中的地震数据的时间的最大值)，根据式(8)可知计算所述点处的地震数据在时间方向的梯度需要使用三维地震数据体外的点处的地震数据(例如，在计算D(x,y,z_max)在时间方向的梯度时，需要使用的(x,y,z_max+1)～(x,y,z_max+(R-1)/2)的地震数据)，由于所述点处无地震数据，所以将所述点处的地震数据的值均取为D(x,y,z_max)。

可以理解，对三维地震数据体中时间小于z_min+(R-1)/2或大于z_max-(R-1)/2的点处的地震数据计算梯度时所使用到的三维地震数据体外的点处的地震数据可根据需要进行取值，例如，使所述点处的地震数据的取值均为0。

在本发明中，对步骤301、步骤302、步骤303的执行顺序不做限制，可以相同或不同。

在步骤304，将步骤301、步骤302、步骤303中得到的线号梯度体、道号梯度体、时间梯度体按线号梯度体、道号梯度体、时间梯度体的顺序纵向排列，以构成梯度体向量。例如，在使用g_x表示线号梯度体，g_y表示道号梯度体，g_z表示时间梯度体时，按线号梯度体、道号梯度体、时间梯度体的顺序纵向排列所构成的梯度体向量g表示为等式(9)：

$g=\left(\begin{array}{c}{g}_x\\ g_y\\ g_z\end{array}\right)---\left(9\right)$

返回图2，在步骤202，基于在步骤201得到的梯度体向量，构建梯度体结构张量方阵。这里，可通过下面的等式(10)构建梯度体结构张量方阵：

$T={\mathrm{gg}}^T=\left(\begin{array}{ccc}{g}_x{g}_x& g_x{g}_y& g_x{g}_z\\ g_y{g}_x& g_y{g}_y& g_y{g}_z\\ g_z{g}_x& g_z{g}_y& g_z{g}_z\end{array}\right)---\left(10\right)$

这里，T为梯度体结构张量方阵，g^T为梯度体向量g的转置。梯度体结构张量方阵T的任意元素g_ug_v表示g_u中的每个位置的元素与g_v中的相同位置的元素相乘，u∈{x,y,z}，v∈{x,y,z}。

本发明中将梯度体g_u与梯度体g_v相乘定义为两个梯度体中的相同位置处的元素相乘。

在图2所示的根据本发明示例性实施例的图1中的根据读取的地震数据获得梯度体结构张量方阵的方法中还可包括对构建的梯度体结构张量方阵进行平滑。可以采用各种平滑方法对构建的梯度体结构张量方阵进行平滑。优选地，采用具有预定尺度因子的三维零均值离散高斯核函数确定的卷积核与梯度体结构张量方阵进行卷积，以得到平滑后的梯度体结构张量方阵。这里，三维零均值离散高斯核函数的表达式为下面的等式(11)：

$G(o,p,q)=e^{[-\frac{{(o/{\mathrm{\σ}}_{3x})}^2+{(p/{\mathrm{\σ}}_{3y})}^2+{(q/{\mathrm{\σ}}_{3z})}^2}{2}]}---\left(11\right)$

这里，G(o,p,q)为三维零均值离散高斯核函数，o为与线号方向对应的离散变量，p为与道号方向对应的离散变量，q为与时间方向对应的离散变量，σ_3x为与线号方向对应的预定尺度因子，σ_3y为与道号方向对应的预定尺度因子，σ_3z为与时间方向对应的预定尺度因子。

三维零均值离散高斯核函数确定的卷积核为将离散变量o、p、q的所有取值的组合(即，由o、p、q表示的所有点，离散变量o、p、q的每个取值组合可表示空间中一个点，点的三维坐标可表示为(o,p,q))代入三维零均值离散高斯核函数得到的函数值的三维矩阵，其中，o、p、q为整数，离散变量o的取值范围为[-R_o，+R_o]，R_o ²＝4²σ_3x；离散变量p的取值范围为[-R_p，+R_p]，R_p ²＝4²σ_3y；离散变量q的取值范围为[-R_q，+R_q]，R_q ²＝4²σ_3z。

应该理解，任意函数值在三维矩阵中的相对位置与该函数值所对应的点(o,p,q)在取值空间中的相对位置相同。点(o,p,q)的取值空间为离散变量o、p、q的取值范围形成的三维空间。

例如，假设σ_3x＝σ_3y＝σ_3z＝1，则R_o＝R_p＝R_q＝4。那么，离散变量o的取值范围为[-4，+4]，o的从小到大的取值为{-4，-3，-2，-1，0，+1，+2，+3，+4}；离散变量p的取值范围为[-4，+4]，p的从小到大的取值为{-4，-3，-2，-1，0，+1，+2，+3，+4}；离散变量q的取值范围为[-4，+4]，q的从小到大的取值为{-4，-3，-2，-1，0，+1，+2，+3，+4}。可将离散变量o、p、q的所有取值的组合代入等式(11)从而得到一个三维矩阵作为卷积核。例如，将点(-4，-4，-4)所对应的o＝-4，p＝-4，q＝-4代入到等式(11)得到函数值G(-4，-4，-4)为3.78E-11。3.78E-11在三维矩阵中的相对位置与点(-4，-4，-4)在由o、p、q的上述取值范围确定的三维空间中的相对位置相同。

优选地，σ_3x可以设为道号相邻的两个道之间的间距的三倍。σ_3y可以设为线号相邻的两条线之间的间距的三倍。σ_3z可以设为时间轴(即，z轴)采样间隔的三倍。

对构建的梯度体结构张量方阵进行平滑可通过下面的等式(12)进行：

$T_s=G_3*T=\left(\begin{array}{ccc}{G}_3*g_x{g}_x& G_3*g_x{g}_y& G_3*g_x{g}_z\\ G_3*g_y{g}_x& G_3*g_y{g}_y& G_3*g_y{g}_z\\ G_3*g_z{g}_x& G_3*g_z{g}_y& G_3*g_z{g}_z\end{array}\right)---\left(12\right)$

这里，T_s为平滑后的梯度体结构张量方阵，G₃为前述的具有预定尺度因子的三维零均值离散高斯核函数确定的卷积核。

图4示出采用本发明示例性实施例的计算地层不连续性属性值的方法的计算效率与传统方法的计算效率对比曲线图。该对比曲线是在同一台计算机设备上测得的，由曲线可看出，计算不连续性属性值的点越多，传统方法所需时间就越大于采用本发明示例性实施例的方法所需时间。图5A示出未采用本发明示例性实施例的计算地层不连续性属性值的方法的不连续性属性效果图。图5B示出采用本发明示例性实施例的计算地层不连续性属性值的方法的不连续性属性效果图。对比可知，本发明示例性实施例的计算地层不连续性属性值的方法的抗噪性能强，且通过计算出的地层不连续性属性值能够清晰地反映断层和岩性的变化。

图6是示出根据本发明示例性实施例的计算地层不连续性属性值的装置的方框图。

如图6所示，根据本发明示例性实施例的计算地层不连续性属性值的装置包括：读取数据单元601、梯度体结构张量方阵获得单元602、结构张量方阵构成单元603、特征值计算单元604、不连续性属性值计算单元605。

读取数据单元601，读取三维地震数据体中各个点所对应的地震数据。图9示出三维地震数据体的空间模型。在三维地震勘探中，得到的三维地震数据体的空间模型如图9所示，图9中的黑点代表地震数据，可将地震数据表示为D(x,y,z)。这里的x表示线号，y表示道号，z表示时间。地震数据D可理解为由线号(Inline)、道号(Xline)、时间(Time)所确定的点(即，采样点)所对应的数据。在这里，读取数据单元601所读取的地震数据可为叠后地震数据，具体地可以为振幅、瞬时相位等。可以理解，x、y、z为取整数的离散变量。

梯度体结构张量方阵获得单元602，根据读取数据单元601所读取的地震数据，获得梯度体结构张量方阵。这里，梯度体结构张量方阵获得单元602将在后面结合图7和图8进行说明。

结构张量方阵构成单元603，从梯度体结构张量方阵获得单元602得到的梯度体结构张量方阵的每个元素中抽取与三维地震数据体中的一个点的位置相对应的位置处的数据，并将抽取的数据按所对应的元素在梯度体结构张量方阵中的位置排列，以构成所述一个点的结构张量方阵。这里，所说的一个点可以是三维地震数据体中的任意点。根据后述内容可知梯度体结构张量方阵为3×3的方阵，所以所述一个点的结构张量方阵也为3×3的方阵。

特征值计算单元604，计算结构张量方阵构成单元603得到的所述一个点的结构张量方阵的特征值。这里，特征值计算单元604可通过现有的各种方法来计算结构张量方阵的特征值。由于结构张量方阵为3×3的方阵，所以将得到三个特征值。

不连续性属性值计算单元605，基于特征值计算单元604得到的特征值，计算所述一个点的不连续性属性值。这里，所述一个点的不连续性属性值可通过前述等式(1)计算。

此外，图6所示的根据本发明示例性实施例的计算地层不连续性属性值的装置中的结构张量方阵构成单元603从梯度体结构张量方阵的每个元素中分别抽取与三维地震数据体中的其他各个点的位置相对应的位置处的数据，并分别将抽取的数据按所对应的元素在梯度体结构张量方阵中的位置排列，以分别构成其他各个点的结构张量方阵；特征值计算单元604分别计算其他各个点的结构张量方阵的特征值；不连续性属性值计算单元605分别基于其他各个点的结构张量方阵的特征值，计算其他各个点的不连续性属性值。由此可得到三维地震数据体中所有点的地震数据的不连续性属性值。

下面结合图7和图8对前述梯度体结构张量方阵获得单元602进行详细说明。

图7是示出根据本发明示例性实施例的图6中的梯度体结构张量方阵获得单元的方框图。如图7所示，梯度体结构张量方阵获得单元602包括梯度体向量计算单元701和梯度体结构张量方阵构建单元702。

梯度体向量计算单元701，利用地震数据得到梯度体向量。具体地说，梯度体向量计算单元701分别计算由读取数据单元601读取的三维地震数据体中各个点所对应的地震数据在线号、道号、时间三个方向的梯度，以得到线号梯度体、道号梯度体、时间梯度体，并形成包括线号梯度体、道号梯度体、时间梯度体的梯度体向量。

在本发明的一个实施例中，可根据图8中所示的梯度体向量计算单元701来计算并获得梯度体向量。

图8是示出根据本发明示例性实施例的图7中的梯度体向量计算单元的方框图。

如图8所示，梯度体向量计算单元701包括线号梯度体计算单元801、道号梯度体计算单元802、时间梯度体计算单元803、排列单元804。

线号梯度体计算单元801，计算三维地震数据体中每个点所对应的地震数据在线号方向的梯度，并将各个点对应的地震数据的梯度按对应点在三维地震数据体中的位置排列，以构成线号梯度体。这里，线号梯度体计算单元801可通过各种梯度计算方法来计算三维地震数据体中各个点所对应的地震数据在线号方向的梯度，从而得到线号梯度体。

优选地，线号梯度体计算单元801利用本发明提出的方式来计算所述各个点的地震数据在线号方向的梯度。具体地说，线号梯度体计算单元801通过将与线号方向对应的卷积核与三维地震数据体中每个点所对应的地震数据沿线号方向进行卷积，以得到三维地震数据体中各个点所对应的地震数据在线号方向的梯度。

道号梯度体计算单元802，计算三维地震数据体中每个点所对应的地震数据在道号方向的梯度，并将各个点对应的地震数据的梯度按对应点在三维地震数据体中的位置排列，以构成道号梯度体。这里，道号梯度体计算单元802可通过各种梯度计算方法来计算三维地震数据体中各个点所对应的地震数据在道号方向的梯度，从而得到道号梯度体。

优选地，道号梯度体计算单元802利用本发明提出的方式来计算所述各个点的地震数据在道号方向的梯度。具体地说，道号梯度体计算单元802通过将与道号方向对应的卷积核与三维地震数据体中每个点所对应的地震数据沿道号方向进行卷积，以得到三维地震数据体中各个点所对应的地震数据在道号方向的梯度。

时间梯度体计算单元803，计算三维地震数据体中每个点所对应的地震数据在时间方向的梯度，并将各个点对应的地震数据的梯度按对应点在三维地震数据体中的位置排列，以构成时间梯度体。这里，时间梯度体计算单元803可通过各种梯度计算方法来计算三维地震数据体中各个点所对应的地震数据在时间方向的梯度，从而得到时间梯度体。

优选地，时间梯度体计算单元803利用本发明提出的方式来计算所述各个点的地震数据在时间方向的梯度。具体地说，时间梯度体计算单元803通过将与时间方向对应的卷积核与三维地震数据体中每个点所对应的地震数据沿时间方向进行卷积，以得到三维地震数据体中各个点所对应的地震数据在时间方向的梯度。

排列单元804，将线号梯度体计算单元801、道号梯度体计算单元802、时间梯度体计算单元803得到的线号梯度体、道号梯度体、时间梯度体按线号梯度体、道号梯度体、时间梯度体的顺序纵向排列，以构成梯度体向量。例如，在使用g_x表示线号梯度体，g_y表示道号梯度体，g_z表示时间梯度体时，按线号梯度体、道号梯度体、时间梯度体的顺序纵向排列所构成的梯度体向量g表示为前述等式(9)。

返回图7，梯度体结构张量方阵构建单元702，基于梯度体向量计算单元701得到的梯度体向量，构建梯度体结构张量方阵。这里，可通过前述等式(10)构建梯度体结构张量方阵。

在图7所示的根据本发明示例性实施例的图6中的梯度体结构张量方阵获得单元602还可包括平滑单元，用于对梯度体结构张量方阵构建单元702所构建的梯度体结构张量方阵进行平滑。平滑单元可以采用各种平滑方法对构建的梯度体结构张量方阵进行平滑。优选地，采用具有预定尺度因子的三维零均值离散高斯核函数确定的卷积核与梯度体结构张量方阵进行卷积，以得到平滑后的梯度体结构张量方阵。这里，三维零均值离散高斯核函数的表达式为前述等式(11)。

三维零均值离散高斯核函数确定的卷积核为将离散变量o、p、q的所有取值的组合(即，由o、p、q表示的所有点，离散变量o、p、q的每个取值组合可表示空间中一个点，点的三维坐标可表示为(o,p,q))代入三维零均值离散高斯核函数得到的函数值的三维矩阵，其中，o、p、q为整数，离散变量o的取值范围为[-R_o，+R_o]，R_o ²＝4²σ_3x；离散变量p的取值范围为[-R_p，+R_p]，R_p ²＝4²σ_3y；离散变量q的取值范围为[-R_q，+R_q]，R_q ²＝4²σ_3z。

应该理解，任意函数值在三维矩阵中的相对位置与该函数值所对应的点(o,p,q)在取值空间中的相对位置相同。点(o,p,q)的取值空间为离散变量o、p、q的取值范围形成的三维空间。

优选地，σ_3x可以设为道号相邻的两个道之间的间距的三倍。σ_3y可以设为线号相邻的两条线之间的间距的三倍。σ_3z可以设为时间轴(即，z轴)采样间隔的三倍。

对构建的梯度体结构张量方阵进行平滑可通过前述等式(12)进行。

根据本发明的计算地层不连续性属性值的装置，能够利用前述计算地层不连续性属性值的方法快速准确的计算出地层的不连续性属性值，从而通过计算出的地层不连续性属性值能够清晰地反映断层和岩性的变化。

应该理解，根据本发明的示例性实施例的根据地震数据确定方向性的装置中的各个单元可被实现为硬件组件。本领域技术人员根据限定的各个单元所执行的处理，可以使用例如现场可编程门阵列(FPGA)或专用集成电路(ASIC)来实现各个单元。

此外，根据本发明的上述方法可以被实现为计算机可读记录介质中的计算机代码。本领域技术人员可以根据对上述方法的描述来实现所述计算机代码。当所述计算机代码在计算机中被执行时实现本发明的上述方法。

本发明的以上实施例仅仅是示例性的，而本发明并不受限于此。本领域技术人员应该理解：在不脱离本发明的原理和精神的情况下，可对这些实施例进行改变，其中，本发明的范围在权利要求及其等同物中限定。

Claims (20)

1.一种计算地层不连续性属性值的方法，所述方法包括：

a)读取三维地震数据体中各个点所对应的地震数据；

b)根据读取的地震数据，获得梯度体结构张量方阵；

c)从梯度体结构张量方阵的每个元素中抽取与三维地震数据体中的一个点的位置相对应的位置处的数据，并将抽取的数据按所对应的元素在梯度体结构张量方阵中的位置排列，以构成所述一个点的结构张量方阵；

d)计算所述一个点的结构张量方阵的特征值；

e)基于所述特征值，计算所述一个点的不连续性属性值。

2.如权利要求1所述的方法，其特征在于，所述一个点的不连续性属性值通过如下等式获得：

$c_n=\frac{\underset{i=0}{\overset{2}{\mathrm{\Σ}}}\underset{j=0}{\overset{2}{\mathrm{\Σ}}}{T}_{\mathrm{ns}}(i,j)}{3({\mathrm{\λ}}_{n1}+{\mathrm{\λ}}_{n2}+{\mathrm{\λ}}_{n3})}$

其中，c_n为所述一个点的不连续性属性值，T_ns(i,j)为所述一个点的结构张量方阵中位于第i行第j列的元素，λ_n1、λ_n2、λ_n3为所述一个点的结构张量方阵的三个特征值。

3.如权利要求1所述的方法，其特征在于，所述方法还包括：

f)分别针对三维地震数据体中的其他点，重复步骤c)-e)来计算所述其他点的不连续性属性值。

4.如权利要求1所述的方法，其特征在于，步骤b)包括：

b1)分别计算三维地震数据体中各个点所对应的地震数据在线号、道号、时间三个方向的梯度，以得到线号梯度体、道号梯度体、时间梯度体，并形成包括线号梯度体、道号梯度体、时间梯度体的梯度体向量；

b2)基于所述梯度体向量，构建梯度体结构张量方阵。

5.如权利要求4所述的方法，其特征在于，步骤b)还包括：

b3)对构建的梯度体结构张量方阵进行平滑。

6.如权利要求4所述的方法，其特征在于，在步骤b1)中，

得到线号梯度体的步骤包括：计算三维地震数据体中每个点所对应的地震数据在线号方向的梯度，并将各个点对应的地震数据的梯度按对应点在三维地震数据体中的位置排列，以构成线号梯度体；

得到道号梯度体的步骤包括：计算三维地震数据体中每个点所对应的地震数据在道号方向的梯度，并将各个点对应的地震数据的梯度按对应点在三维地震数据体中的位置排列，以构成道号梯度体；

得到时间梯度体的步骤包括：计算三维地震数据体中每个点所对应的地震数据在时间方向的梯度，并将各个点对应的地震数据的梯度按对应点在三维地震数据体中的位置排列，以构成时间梯度体；

形成梯度体向量的步骤包括：将线号梯度体、道号梯度体、时间梯度体按线号梯度体、道号梯度体、时间梯度体的顺序纵向排列，以构成梯度体向量。

7.如权利要求6所述的方法，其特征在于，

计算三维地震数据体中每个点所对应的地震数据在线号方向的梯度的步骤包括：将与线号方向对应的卷积核与三维地震数据体中每个点所对应的地震数据沿线号方向进行卷积，以得到三维地震数据体中各个点所对应的地震数据在线号方向的梯度；

计算三维地震数据体中每个点所对应的地震数据在道号方向的梯度的步骤包括：将与道号方向对应的卷积核与三维地震数据体中每个点所对应的地震数据沿道号方向进行卷积，以得到三维地震数据体中各个点所对应的地震数据在道号方向的梯度；

计算三维地震数据体中每个点所对应的地震数据在时间方向的梯度的步骤包括：将与时间方向对应的卷积核与三维地震数据体中每个点所对应的地震数据沿时间方向进行卷积，以得到三维地震数据体中各个点所对应的地震数据在时间方向的梯度。

8.如权利要求7所述的方法，其特征在于，与线号方向对应的卷积核、与道号方向对应的卷积核、与时间方向对应的卷积核分别通过将一维零均值离散高斯核函数的导数在离散变量为相应的取值范围内的各个整数值时的函数值按对应的离散变量的从小到大的顺序排列构成，所述各卷积核中的函数值的计算式为：

$G^{\′}\left(t\right)=-\frac{1}{2\mathrm{\π}{{\mathrm{\σ}}_i}^2}t{e}^{[-t^2/\left(2{{\mathrm{\σ}}_i}^2\right)]}$

其中，G′(t)为一维零均值离散高斯核函数的导数，t为离散变量，t的取值范围为[-R_i，+R_i]，R_i ²＝4²σ_i，R_i为核半径，σ_i为预定尺度因子，i∈{x,y,z}，σ_x为与线号方向对应的预定尺度因子，R_x为与线号方向对应的核半径，σ_y为与道号方向对应的预定尺度因子，R_y为与道号方向对应的核半径，σ_z为与时间方向对应的预定尺度因子，R_z为与时间方向对应的核半径。

9.如权利要求4所述的方法，其特征在于，在步骤b2)中通过下面的等式构建梯度体结构张量方阵：

$T={\mathrm{gg}}^T=\left(\begin{array}{ccc}{g}_x{g}_x& g_x{g}_y& g_x{g}_z\\ g_y{g}_x& g_y{g}_y& g_y{g}_z\\ g_z{g}_x& g_z{g}_y& g_z{g}_z\end{array}\right)$

其中，T为梯度体结构张量方阵，g为梯度体向量， $g=\left(\begin{array}{c}{g}_x\\ g_y\\ g_z\end{array}\right),$ g_x为线号梯度体，g_y为道号梯度体，g_z为时间梯度体，g^T为梯度体向量的转置，梯度体结构张量方阵T的任意元素g_ug_v表示g_u中的每个位置的元素与g_v中的相同位置的元素相乘，u∈{x,y,z}，v∈{x,y,z}。

10.如权利要求5所述的方法，其特征在于，对构建的梯度体结构张量方阵进行平滑的步骤包括：使用具有预定尺度因子的三维零均值离散高斯核函数确定的卷积核与梯度体结构张量方阵进行卷积，以得到平滑后的梯度体结构张量方阵。

11.一种计算地层不连续性属性值的装置，所述装置包括：

读取数据单元，读取三维地震数据体中各个点所对应的地震数据；

梯度体结构张量方阵获得单元，根据读取的地震数据，获得梯度体结构张量方阵；

结构张量方阵构成单元，从梯度体结构张量方阵的每个元素中抽取与三维地震数据体中的一个点的位置相对应的位置处的数据，并将抽取的数据按所对应的元素在梯度体结构张量方阵中的位置排列，以构成所述一个点的结构张量方阵；

特征值计算单元，计算所述一个点的结构张量方阵的特征值；

不连续性属性值计算单元，基于所述特征值，计算所述一个点的不连续性属性值。

12.如权力要求11所述的装置，其特征在于，不连续性属性值计算单元通过如下等式获得所述一个点的不连续性属性值：

$c_n=\frac{\underset{i=0}{\overset{2}{\mathrm{\Σ}}}\underset{j=0}{\overset{2}{\mathrm{\Σ}}}{T}_{\mathrm{ns}}(i,j)}{3({\mathrm{\λ}}_{n1}+{\mathrm{\λ}}_{n2}+{\mathrm{\λ}}_{n3})}$

其中，c_n为所述一个点的不连续性属性值，T_ns(i,j)为结构张量方阵构成单元所构成的所述一个点的结构张量方阵中位于第i行第j列的元素，λ_n1、λ_n2、λ_n3为所述一个点的结构张量方阵的三个特征值。

13.如权利要求11所述的装置，其特征在于，结构张量方阵构成单元从梯度体结构张量方阵的每个元素中分别抽取与三维地震数据体中的其他各个点的位置相对应的位置处的数据，并分别将抽取的数据按所对应的元素在梯度体结构张量方阵中的位置排列，以分别构成其他各个点的结构张量方阵；特征值计算单元分别计算其他各个点的结构张量方阵的特征值；不连续性属性值计算单元分别基于其他各个点的结构张量方阵的特征值，计算其他各个点的不连续性属性值。

14.如权利要求11所述的装置，其特征在于，所述梯度体结构张量方阵获得单元包括：

梯度体向量计算单元，分别计算三维地震数据体中各个点所对应的地震数据在线号、道号、时间三个方向的梯度，以得到线号梯度体、道号梯度体、时间梯度体，并形成包括线号梯度体、道号梯度体、时间梯度体的梯度体向量；

梯度体结构张量方阵构建单元，基于所述梯度体向量，构建梯度体结构张量方阵。

15.如权利要求14所述的装置，其特征在于，所述梯度体结构张量方阵获得单元还包括：

平滑单元，对梯度体结构张量方阵构建单元所构建的梯度体结构张量方阵进行平滑。

16.如权利要求14所述的装置，其特征在于，梯度体向量计算单元包括：

线号梯度体计算单元，计算三维地震数据体中每个点所对应的地震数据在线号方向的梯度，并将各个点对应的地震数据的梯度按对应点在三维地震数据体中的位置排列，以构成线号梯度体；

道号梯度体计算单元，计算三维地震数据体中每个点所对应的地震数据在道号方向的梯度，并将各个点对应的地震数据的梯度按对应点在三维地震数据体中的位置排列，以构成道号梯度体；

时间梯度体计算单元，计算三维地震数据体中每个点所对应的地震数据在时间方向的梯度，并将各个点对应的地震数据的梯度按对应点在三维地震数据体中的位置排列，以构成时间梯度体；

排列单元，将线号梯度体、道号梯度体、时间梯度体按线号梯度体、道号梯度体、时间梯度体的顺序纵向排列，以构成梯度体向量。

17.如权利要求16所述的装置，其特征在于，

线号梯度体计算单元通过将与线号方向对应的卷积核与三维地震数据体中每个点所对应的地震数据沿线号方向进行卷积，以得到三维地震数据体中各个点所对应的地震数据在线号方向的梯度；

道号梯度体计算单元通过将与道号方向对应的卷积核与三维地震数据体中每个点所对应的地震数据沿道号方向进行卷积，以得到三维地震数据体中各个点所对应的地震数据在道号方向的梯度；

时间梯度体计算单元通过将与时间方向对应的卷积核与三维地震数据体中每个点所对应的地震数据沿时间方向进行卷积，以得到三维地震数据体中各个点所对应的地震数据在时间方向的梯度。

18.如权利要求17所述的装置，其特征在于，与线号方向对应的卷积核、与道号方向对应的卷积核、与时间方向对应的卷积核分别通过将一维零均值离散高斯核函数的导数在离散变量为相应的取值范围内的各个整数值时的函数值按对应的离散变量的从小到大的顺序排列构成，所述各卷积核中的函数值的计算式为：

$G^{\′}\left(t\right)=-\frac{1}{2\mathrm{\π}{{\mathrm{\σ}}_i}^2}t{e}^{[-t^2/\left(2{{\mathrm{\σ}}_i}^2\right)]}$

其中，G′(t)为一维零均值离散高斯核函数的导数，t为离散变量，t的取值范围为[-R_i，+R_i]，R_i ²＝4²σ_i，R_i为核半径，σ_i为预定尺度因子，i∈{x,y,z}，σ_x为与线号方向对应的预定尺度因子，R_x为与线号方向对应的核半径，σ_y为与道号方向对应的预定尺度因子，R_y为与道号方向对应的核半径，σ_z为与时间方向对应的预定尺度因子，R_z为与时间方向对应的核半径。

19.如权利要求14所述的装置，其特征在于，梯度体结构张量方阵构建单元通过下面的等式构建梯度体结构张量方阵：

$T={\mathrm{gg}}^T=\left(\begin{array}{ccc}{g}_x{g}_x& g_x{g}_y& g_x{g}_z\\ g_y{g}_x& g_y{g}_y& g_y{g}_z\\ g_z{g}_x& g_z{g}_y& g_z{g}_z\end{array}\right)$

其中，T为梯度体结构张量方阵，g为梯度体向量， $g=\left(\begin{array}{c}{g}_x\\ g_y\\ g_z\end{array}\right),$ g_x为线号梯度体，g_y为道号梯度体，g_z为时间梯度体，g^T为梯度体向量的转置，梯度体结构张量方阵T的任意元素g_ug_v表示g_u中的每个位置的元素与g_v中的相同位置的元素相乘，u∈{x,y,z}，v∈{x,y,z}。

20.如权利要求15所述的装置，其特征在于，平滑单元使用具有预定尺度因子的三维零均值离散高斯核函数确定的卷积核与梯度体结构张量方阵进行卷积，以得到平滑后的梯度体结构张量方阵。