本发明是关于一种基于结构张量导引的高维s变换方法，涉及地震勘探领域，特别涉及高维、冗余地震数据的特征稀疏表达领域。
背景技术
时频分析作为新兴的非平稳信号处理方法，近年来受到越来越多的关注，时频分析或称时频分布，是描述信号频率随时间变化的信号处理方法。采用时间-频率联合分布表示信号，将一维的时间域信号映射到一个二维的时频平面，在时频域内对信号进行分析，全面反映观测信号的时间-频率联合特征，使研究人员同时掌握信号的时域及频域信息，而且还可以了解信号的频率是如何随时间变化的。
传统的时频变换方法有短时傅里叶变换、gabor变换、小波变换、s变换等。其中，s变换结合了短时傅立叶变换和小波变换的优点而避免了它们的不足：频率的倒数决定了s变换中的高斯窗的尺度大小，因此具有小波变换的多分辨率分析；而且s变换含有相位因子，保留了每个频率的绝对相位特征，这是小波变换所不具有的特性。另外，s变换对于多分量信号具有可加性，是一种线性变换，具有无损可逆性。但是，一维的时频变换缺乏空间分辨率。此外，传统时频变换只是在基函数种类固定的局限下，在处理非平稳信号时，其窗口调节的灵活性略显不足，物理意义也不明确。因此，有必要将一维s变换推广到高维情况。
技术实现要素：
针对上述问题，本发明的目的是提供一种能够突显地震数据稀疏特征，增加地震信息的利用率，有助于地震数据属性提取的基于结构张量导引的高维s变换方法。
为实现上述目的，本发明采取以下技术方案：一种基于结构张量导引的高维s变换方法，包括以下内容：
步骤1)：构建结构张量矩阵，对地震数据通过结构张量矩阵的特征值与特征向量进行各向异性参数提取；
步骤2)：构建基于结构张量导引的高维s变换算子；
步骤3)：将地震数据采用基于结构张量导引的高维s变换算子进行高维s变换，得到符合地震数据特征的稀疏表征数据。
进一步地，所述步骤1)的具体过程为；假设由地震振幅构成的标量场为u，定义结构张量矩阵：
其中，x，z为地震数据的空间坐标；
结构张量矩阵s的特征矩阵为：
(s-λe)x＝0
其中，λ为矩阵s对应的特征值，e为单位矩阵，x则为矩阵s的特征向量；结构张量的特征值为：
根据特征值得到特征向量：
通过结构张量的特征值与特征向量，得到地震数据的各向异性参数：
其中，p,α,β为各向异性参数，ε1和ε2为平衡系数，取值在0～1之间。进一步地，高维s变换算子为：
(x-x′)p＝(x-x′)·cosp+(z-z′)·sinp
(z-z′)p＝-(x-x′)·sinp+(z-z′)·cosp
式中，f为频率，h(x,z)为地震数据，(x,z)与(x′,z′)分别为s变换前后的地震数据空间坐标。
本发明由于采取以上技术方案，其具有以下优点：1、本发明引入了结构张量的概念，增加基函数灵活性的同时，更加稀疏的突显高维、冗余地震数据的特征，从而将物理意义明确、稀疏性更强的地震数据用于地震资料处理解释中，增加地震信息的利用率，提高处理解释的效果。2、本发明将结构张量矩阵与高维s变换相结合，凸显地震数据特征，得到更加符合地震数据特征的稀疏表征数据。3、本发明可以根据结构张量分析出的地震数据的特征信息，对地震数据采用改进的高维s变换，从而使地震数据特征更加突显、细化以及稀疏，提高了地震数据在处理解释时的信息利用率，有助于地震数据属性(方向、尺度、走时等)的提取，具有重要的理论和应用价值，对油气、矿产资源勘探的发展有一定的推动作用。
附图说明
图1是本发明的基于结构张量导引的高维s变换方法流程图；
图2是本发明实施例包含3个地震同相轴的地震剖面效果图；
图3是本发明实施例基于结构张量导引的高维s变换稀疏特征表达的地震数据效果图，图(a)、(b)和(c)分别代表同一频率下，不同方向的特征表达的地震数据效果图。
具体实施方式
以下结合附图来对本发明进行详细的描绘。然而应当理解，附图的提供仅为了更好地理解本发明，它们不应该理解成对本发明的限制。
如图1所示，本发明的基于结构张量导引的高维s变换方法，包括以下内容：
1、构建结构张量矩阵，对地震数据通过结构张量矩阵的特征值与特征向量进行各向异性参数提取。
最初，由dizenzo和forstner先后于1986年和1987年提出结构张量的定义，并将其应用于特征检测。假设由地震振幅构成的标量场为u，定义结构张量矩阵(梯度平方矩阵)：
其中，x，z为地震数据的空间坐标。
结构张量矩阵s的特征矩阵为：
(s-λe)x＝0(2)
其中，λ为结构张量矩阵s对应的特征值，e为单位矩阵，x则为结构张量矩阵s的特征向量。
有非零解的充分必要条件是(2)对应的特征多项式为零，即
利用求根公式：
结构张量的特征值描述了图像在一定邻域内在特征向量方向上的数值变化平均量。当λ1＞＞λ2≈0时，表示该点附近的图像呈线性结构；当λ1＞＞0且λ2＞＞0时，表示该点附近至少存在两个线性结构；当λ1≈0且λ2≈0时，表示该点附近图像缓变，没有线性结构。
根据特征值(4)可以得到特征向量：
较大特征值对应的特征向量x1的方向为图像中梯度变化最大的方向(orthogonaltolinearfeatures)，较小特征值对应的特征向量x2与x1垂直(paralleltolinearfeatures)。
通过结构张量给出的特征值与特征向量，便能根据公式(6)得到地震数据的各向异性参数：
其中，p,α,β为各向异性参数，ε1和ε2为平衡系数，调节各向异性参数的幅值，一般在0～1之间。
2、构建物理意义更加明确的基于结构张量导引的高维s变换算子。
传统的高维s变换(2d为例)是stockwell等人1997年在一维s变换的基础上提出的，二维s变换继承了一维s变换的优点，在空域与频域都有很高的分辨率。这使得它在二维图像处理领域具有广阔的应用前景。
传统二维s变换可写成：
其中，h(x,z)为地震数据，(x,z)与(x′,z′)分别为s变换前后的地震数据空间坐标，kx和kz为x方向和z方向上的波数，公式(7)可以看作是两次一维s变换。
在传统二维s变换基础上，将传统的时频分析扩展到时间-空间频率-各向异性空间，力图通过扩展空间，进一步凸显细化数据的某些特征。然后，引入了结构张量矩阵，在高维、冗余地震信号中稀疏提取特征明确的数据。
改进后的高维s变换算子为：
(x-x′)p＝(x-x′)·cosp+(z-z′)·sinp(9)
(z-z′)p＝-(x-x′)·sinp+(z-z′)·cosp(10)
其中，f表示频率，与公式(7)对比，在传统s变换的基础上，引进了多个物理意义明确的参数，可以反映地震数据的空间方向，尺度信息。
通过公式(8)，便可以扩展时频空间到更高维更冗余的数据空间，从而将地震数据的特征细化突显。
3、采用基于结构张量导引的高维s变换对地震数据进行稀疏特征表达。
改进后的高维s变换公式(8)中的各向异性参数(p,α,β)可以通过公式(6)得到。将地震数据采用基于结构张量导引的高维s变换算子进行高维s变换，便能得到更加符合地震数据特征的稀疏表征数据。
如图2～3所示，下面通过具体实施例详细说明本发明的基于结构张量导引的高维s变换方法的具体过程。本实施例对包含3个地震同相轴的地震剖面进行测试，具体过程：
(1)输入地震剖面数据，根据公式(4)和(5)计算得到对应地震剖面的结构张量矩阵的特征值与特征向量，从而得到相应的地震数据的各向异性参数(p,α,β)；
(2)根据给出频率值f，一般满足香农采样定理即可。然后，结合地震剖面结构张量矩阵提供的各向异性参数与公式(8)，便能够对地震剖面进行高维s变换，提取特征更加细化的、稀疏的、不同频率、不同方向、不同尺度的地震数据信息，有助于进一步的地震资料处理解释工作。
综上所述，与传统的二维s变换方法相比，本发明的基于结构张量导引的高维s变换方法，不仅含有地震数据的时频-空间信息，还突显了地震数据的方向、尺度信息。这样，使得地震数据得到了细化，针对地震资料处理解释的需求，可以更加稀疏的特征表达地震数据，提高了地震数据信息的利用率。
上述各实施例仅用于说明本发明，其中方法的各实施步骤等都是可以有所变化的，凡是在本发明技术方案的基础上进行的等同变换和改进，均不应排除在本发明的保护范围之外。