一种基于张量的自适应秩截断的地震数据压缩方法

摘要

本发明公开一种基于张量自适应秩截断的地震数据压缩方法，通过给定的压缩条件和即时通过高阶奇异值分解获取的不同维度奇异值设置截断秩，同时根据奇异值的分布决定截断秩大小进行张量分解；保证压缩率的同时，提高压缩的峰值信噪比，改善压缩效果；本发明通过设计基于自适应秩截断的张量分解的流程，减少压缩过程中人为调控参数，简化数据压缩工作；本发明基于奇异值分布设定截断秩的大小，针对数据各向异性的特点，引入奇异值的标准差计算不同为维度对应截断秩的权重，提高压缩后的峰值性噪比，改善数据的压缩效果。

说明书

技术领域

本发明属于地震数据处理技术领域，具体涉及一种地震数据压缩方法。

背景技术

随着技术硬件的发展，获取的地震数据规模日益增长。地震数据的特点也从低维度、单属性，变为高纬度、多属性。庞大的地震数据不仅耗费大量的存储空间，也为数据的处理，可视化带来极大的困难。因此高效的数据压缩技术的提出与应用刻不容缓。对于庞大的地震数据通常采用的是高压缩比的有损(lossy)压缩技术。

现有的数据压缩技术主要是基于的域变换和编码的压缩技术。常见的域变换压缩技术包括傅里叶变换、余弦变换及小波变换压缩技术，其基本思想是将数据由多个不同频率的基波信号表出。这类压缩技术的优点是易于实现，且正变换和逆变换可快速实现。但存在变换域受限，在高压缩比的条件下，压缩效果不理想(较低峰值信噪比)的缺陷。基于域变换技术的改进相继出现Curvelet变换以及Dreamlet变换，这类变换技术也被称为类小波变换。因其具有频率-时间-空间域的局域化性质，因此能很好地表示地震数据自身的特性。这些压缩技术在二维地震数据压缩上取得较好的压缩效果，但未能有效的用于在地震体数据的压缩。基于编码的压缩技术主要包括矢量量化。这类压缩技术依据数据本身生成码书对数据进行近似，因此能应对数据的各向异性，但缺点是码书生成时间长，压缩的预处理过程耗时。

现有的压缩技术的核心是利用有限数目的基和对应系数对数据进行近似，即有效地对数据进行特征提取。其中用于数据近似的基可以分为基于预先定义的和基于学习的。预定义的基的优点是耗费较少的预处理时间。由于基是预先定义的，因此与原始数据相关差，不能有效地提取数据中的特征信息。基于学习的基的压缩技术则耗费更长的与处理时间通过数据自身的特征获取基，因此基更能反映数据特征，获得更高的压缩率。目前针对地震数据处理自身的需求，数据重构具有快速、实时的特点，但数据的压缩可以是耗时的。相比于预定义基的压缩技术，基于学习的压缩技术在保留数据特征的同时，更能有效去除数据中的冗余，取得更好的压缩效果。

张量由Tamara提出用以表示不同维度的数据，张量的阶数能直接反映数据的维度。一个三维数据可以用三阶张量表示。在对张量进行分块后，进行张量模型的运算。

张量分解是通过高阶奇异值(Higher-order SVD，HOSVD)分解对原始数据进行特征提取，实现对原始数据的近似。高阶奇异值分解是奇异值分解运算在高维数据的推广，即将高维数据退化为二维矩阵进行奇异值分解运算。高阶奇异值分解进行特征提取的核心思想是基于主成成分分析实现高维数据降维的推广。

主成分分析(Principal component analysis，PCA)作为现代数据分析的标准工具被用于神经科学到计算机图形学等不同领域。因其是基于数据本身的方法，所以与数据本身具有很强的相关性。

现今张量近似的主要方法基于Tucker模型的张量分解和张量重构。张量的Tucker分解是将张量沿不同的方向展开成矩阵，再对矩阵做奇异值分解，将张量表示为因子矩阵和核心张量。每一个因子矩阵反映数据在该方向的结构特征。相比早先的基于CP模型的张量分解，Tucker分解产生具有数据特征的基，更便于数据的特征提取，因此Tucker分解根据实用性。张量重构基于Tucker分解进行秩截断，将截断后因子矩阵和重新计算的核心张量通过张量的乘积运算还原为张量实现对原始张量的近似。整个过程产生的截断因子矩阵和核心张量的数据规模小于原始数据，整个过程可以看作是对原始数据进行压缩与重构。

在多分辨率数据可视化中，基于Tucker模型的张量近似相比于离散小波变换，在系数数量相同的条件下，张量近似更能反映数据在不同维度的特征。在基于张量近似的多分辨率体绘制应用中，K.Suter提出基于全局因子矩阵的张量近似，在保证较高压缩率的同时，提供整体提取结构特征的方法，使得张量产生的基具有全局性。通过实现GPU加速的张量重构，大大提高了体绘制的实时性。相比于基于域变换的压缩技术，张量近似具有基于原始数据学习的特点，因此与原始数据较强的相关性，在压缩过程中更能有效去除数据中的冗余，保留其结构特征。

目前基于张量分解的压缩在实现流程如图1所示，通常采用预定义截断秩大小，即在进行张量近似前预先设置截断秩的大小。截断秩的大小决定了压缩率和压缩后的数据失真度。截断秩较大虽然能最大程度保留原始数据的特征，但压缩率较小；若截断秩较小，虽然能有较高的压缩率，但数据细节有较大的损失。

对三维地震数据，其数据特点具有各向异性的特点，即数据不同维度的结构特征不同。实际的数据分析显示，对统一地震数据的不同维度进行主成成分分析，其特征大小，分布具有差异。对地震数据的张量分解，对应维度采用不同大小的截断秩能取得更好的压缩效果。然而实际压缩过程中无法事先预知数据的特征，因此基于预定义截断秩大小的张量分解无法产生较好的压缩效果，且需要人为调整压缩参数(截断秩大小)。

发明内容

本发明为解决上述技术问题，提出了一种基于张量的自适应秩截断的地震数据压缩方法，通过奇异值分布决定截断秩对地震数据进行压缩，在给定的压缩率条件下，提高数据压缩效果。

本发明采用的技术方案是：一种基于张量的自适应秩截断的地震数据压缩方法，包括：

S1、对原始三维地震数据进行分块处理，将分块后的数据作为输入的原始张量；

S2、步骤S1输入的每一个原始张量沿三个维度方向展开为二维矩阵，由奇异值分解求得三组不同的奇异值和对应的奇异向量组成的因子矩阵；

S3、根据步骤S2求得的奇异值确定对应截断秩的大小；

S4、根据步骤S3确定的截断秩对因子矩阵进行列选取截断的因子矩阵；

S5、根据截断的因子矩阵求得核心张量，并得到高阶奇异值分解后的三个因子矩阵；

S6、根据高阶正交迭代算法得到优化后的因子矩阵和核心张量。

进一步地，步骤S1所述分块处理选取边长为2的整数幂的立方体。

进一步地，步骤S3还包括选取不同维度的奇异值取标准差量化奇异值的分布，并根据奇异值分布决定截断秩大小。

更进一步地，所述根据奇异值分布决定截断秩大小，具体包括以下分步骤：

A1、输入数据的压缩率C；

A2、计算截断秩总和R；

A3、奇异值取对数，并计算每组奇异值的标准差；

A4、对每组奇异值的标准差求倒数，并计算每组奇异值对应的权重w_i，i表示第i个维度，i＝1，2，3；

A5、根据每组奇异值对应的权重w_i，计算得到每个维度截断秩的大小。

进一步地，步骤S6所述高阶正交迭代算法具体包括以下分步骤：

B1、输入一个原始张量M；

B2、对原始张量M进行高阶奇异值分解得到初始的因子矩阵；

B3、计算原始张量M与步骤S5得到的核心张量的范数差；

B4、将原始张量M、步骤B3得到的范数差以及最大迭代次数作为收敛条件；

B5、根据步骤B4的收敛条件判断当前迭代是否收敛；若是则停止迭代，并输出优化的因子矩阵以及核心张量；否则转至步骤B6；

B6、判断当前迭代是否为第一次迭代，若是，则对步骤B2得到的初始因子矩阵进行转置，并与原始张量M相乘，得到当前迭代的优化后的张量；否则对上一次迭代的因子矩阵进行转置，并与原始张量M相乘，得到当前迭代的优化后的张量；

B7、对步骤B6得到的优化后的张量进行高阶奇异值分解，得到优化后的因子矩阵和核心张量；

B8、计算原始张量M与步骤B7得到的优化后的核心张量的范数差，并转至步骤B5。

本发明的有益效果：本发明通过给定的压缩条件和即时通过高阶奇异值分解获取的不同维度奇异值设置截断秩，同时根据奇异值的分布决定截断秩大小进行张量分解；保证压缩率的同时，提高压缩的峰值信噪比，改善压缩效果；本发明通过设计基于自适应秩截断的张量分解的流程，减少压缩过程中人为调控参数，简化数据压缩工作；本发明基于奇异值分布设定截断秩的大小，针对数据各向异性的特点，引入奇异值的标准差计算不同为维度对应截断秩的权重，提高压缩后的峰值性噪比，改善数据的压缩效果。

附图说明

图1为现有的基于张量分解的压缩流程图。

图2为本发明提供的基于自适应秩截断的张量分解流程图。

图3为本发明提供的基于奇异值分布的截断秩大小流程图。

图4为本发明提供的高阶正交迭代算法流程图。

具体实施方式

为便于本领域技术人员理解本发明的技术内容，下面结合附图对本发明内容进一步阐释。

本发明通过如图2所示的基于自适应秩截断的张量分解方案，具体包括：

S1、对原始三维地震数据进行分块处理，将分块后的数据作为输入的原始张量；分块通常选取边长为2的整数幂的立方体。三维数据的不同维度实际指的就是长宽高，维度为数学描述的名词。

S2、步骤S1输入的每一个原始张量沿三个维度方向展开为二维矩阵，对展开的三个二维矩阵由奇异值分解得到三组不同的奇异值和对应的奇异向量组成的因子矩阵。依据主成分分析的原理，在对数据进行高阶奇异值分解后得到因子矩阵和奇异值。因子矩阵的奇异向量对应的奇异值代表该奇异向量在所有奇异向量中的权重，在做秩截断时，选取奇异值大对应的向量组成截断后的因子矩阵。

S3、根据步骤S2求得的奇异值确定对应截断秩的大小；在对不同维度的因子矩阵进行秩截断时，依据奇异值的分布确定截断秩的大小。奇异值的分布反映了数据在不同维度的特征分布：奇异值分布分散，数据在该维度的特征分布明显，需要提取的特征较少；反之数据在该维度特征不明显，需要提取特征较多。

本发明选取不同维度的奇异值取标准差量化奇异值的分布，基于奇异值分布决定截断秩大小的流程如图3所示。首先对于奇异值做10为底数的对数运算。对不同维度(mode-n)的奇异值分别计算标准差s，对每组(共三个维度即三组)奇异值的标准差求倒数，记为λ；根据式2-1计算每个维度奇异值对应的权重w_i：

数据的压缩率C由式2-2求得，其中，I_i表示原始张量(数据)在第i个维度的大小，i＝1，2，3，R_i表示不同维度的截断秩大小。因此当不同维度的截断秩的和R为定值时，若不同维度的截断秩相等时，压缩率C有最小值。

由公式2-2可知：分子部分为定值，因为I_i的大小本身是确定的，分母中因为假设截断秩R_i的和R为定值，所以分母左半部分为定值，则右半部分是决定C大小的变量。由正数极值定理(2个或以上推广)：R_i的和R为定值时，取R_i相等时，i＝1，2，3，取得极大值。则本申请式2-2中分母极大，C取极小，这正是本申请的期望。在现有的方法中都是将R_i提前取一样的大小，本申请的方法是先求得一组合适的R_i，并且R_i在大多数情况下是互不相同的，然后根据这组R_i之和，确定一个相同的R_i大小，这样不仅能保证压缩率C，而且可以使得对应的压缩效果(峰值信噪比)也得到提升。

则每个维度截断秩R_i大小可由式2-3求得，其中，[]为取整运算：

R_i＝[w_iR]>

S4、根据步骤S3确定的截断秩对步骤S2得到的因子矩阵进行列选取截断的因子矩阵；例如一个分块32*32*32，得到不截断的三个因子矩阵都是32*32(矩阵理论奇异值分解就中是左奇异矩阵)，现在求得的截断秩为7,15,14。那么截断后的三个因子矩阵大小分别为32*7,32*15,32*14。

S5、根据截断的因子矩阵求得核心张量，并得到高阶奇异值分解后的三个因子矩阵；根据截断后的三个因子矩阵大小分别为32*7,32*15,32*14，求得对核心张量是一个大小为7*15*14的长方体张量。

S6、根据高阶正交迭代算法得到优化后的因子矩阵和核心张量。高阶正交迭代算法(Higher-order Orthogonal Iteration，HOOI)利用交替最小二乘原理对因子矩阵和核心张量进行优化，其流程如图4所示。迭代算法以原始张量(块数据)与核心张量的F-范数差、迭代次数作为收敛条件，原始张量与因子矩阵的转置做张量乘积(Tensor Times Matrix Multiplication，TTM积)得到优化后的张量P，再对P求高阶奇异值分解得到优化后的因子矩阵和核心张量，因子矩阵和核心张量在中张量重构中做张量乘积运算，得到近似原始张量(数据)。

步骤S6所述高阶正交迭代算法采用交替最小二乘原理对因子矩阵和核心张量进行优化，以其中一个原始张量的处理为例进行说明，具体包括以下分步骤：

B1、输入一个原始张量M；

B2、对原始张量M进行高阶奇异值分解得到初始的因子矩阵；

B3、计算原始张量M与步骤S5得到的核心张量的范数差；F-范数差是原始张量与重构张量的范数差。但因为重构张量由因子矩阵和核心张量计算得到，因子矩阵本身是正交性，由矩阵理论中范数不变定理，转化为原始张量与核心张量的范数差。

B4、将原始张量M、步骤B3得到的范数差以及最大迭代次数作为收敛条件；

B5、根据步骤B4的收敛条件判断当前迭代是否收敛；若是则停止迭代，并输出优化的因子矩阵以及核心张量；否则转至步骤B6；

B6、判断当前迭代是否为第一次迭代，若是，则对步骤B2得到的初始因子矩阵进行转置，并与原始张量M相乘，得到当前迭代的优化后的张量；否则对上一次迭代的因子矩阵进行转置，并与原始张量M相乘，得到当前迭代的优化后的张量。

B7、对步骤B6得到的优化后的张量进行高阶奇异值分解，得到优化后的因子矩阵和核心张量；

B8、计算原始张量M与步骤B7得到的优化后的核心张量的范数差，并转至步骤B5。

从上述实施例可以看出，本发明提出了基于张量的自适应秩截断的地震数据压缩方法，通过奇异值分布决定截断秩对地震数据进行压缩，在给定的压缩率条件下，提高数据压缩效果，可以达到以下有益效果：

通过设计基于自适应秩截断的张量分解的流程，减少压缩过程中人为调控参数，简化数据压缩工作。

本发明基于奇异值分布设定截断秩的大小，即利用不同维度方向的奇异值求得对应标准差，以标准差的倒数作为权重值，在给定压缩率的条件下计算不同维度截断秩的大小，针对数据各向异性的特点，引入奇异值的标准差计算不同维度对应截断秩的权重，提高压缩后的峰值性噪比，改善数据的压缩效果。

本领域的普通技术人员将会意识到，这里所述的实施例是为了帮助读者理解本发明的原理，应被理解为本发明的保护范围并不局限于这样的特别陈述和实施例。对于本领域的技术人员来说，本发明可以有各种更改和变化。凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的权利要求范围之内。