CN107942377A - 一种地震数据压缩与重构方法 - Google Patents

Abstract

本发明涉及一种地震数据压缩与重构方法，包括以下步骤：首先，对地震数据小波变换，增加地震数据的可压缩性；然后，根据混沌序列构造出硬件可实现的测量矩阵，并用混沌测量矩阵对小波变换后的地震数据压缩观测；最后，对贝叶斯小波树结构压缩感知重构算法进行改进，用改进后的BTSWCS‑vb算法恢复出完整地震数据。本发明实际数据试验结果表明：相比于常用随机测量矩阵，本发明的混沌测量矩阵是便于硬件实现的，几种混沌序列矩阵不仅可以实现地震数据压缩比为0.2～0.55的实时压缩，而且对重构效果影响也很好，特别是Logistic序列矩阵；本发明改进的重构算法BTSWCS‑vb算法，不仅提高了重构精度，而且明显缩短了重构时间。

Description

一种地震数据压缩与重构方法

技术领域

本发明涉及一种地震数据的处理方法，特别涉及一种地震数据压缩与重构方法，基于贝叶斯压缩感知改进算法。

背景技术

地质勘探过程中海量数据的实时回收给现有地震仪带来了巨大挑战，这也是制约地震仪发展的最主要因素。特别是，随着地震勘探的不断深入，使得地震勘探数据持续膨胀，给地震仪的采样和无线传输速度、存储器的存储容量以及计算机的处理速度造成极大压力。通过压缩感知理论的关键技术，在下位机对地震数据实时压缩，数据传输至上位机，当需要地震资料解释的时候上位机又能够高精度重构出完整原始数据，这不但可以提高数据处理速度和减小存储空间，而且可以改善地震仪无线通讯数据传输性能。因此，地震数据的压缩与重构方法成为迫切需要研究的课题。

CN106772567A公开了《一种用于地震勘探仪的数据传输无损压缩算法》，压缩的实现分两个方面：第一方面是根据前两个采样点对本采样点的预测值与本采样点的真实值之差求得预测差值ΔV，用ΔV替代本采样点24位原始数据；第二方面是根据数据特点，采用游程编码和霍夫曼编码两次压缩指示预测差值ΔV数据位数的指示数据位数段，显著减少一帧中指示数据位数段的长度。其采用了多次压缩，压缩能力强，所需上传数据量减少50％以上。该方法具有在算法的复杂度和时效性方面，计算和实现简单，不改变采集站微处理器实时性的优点。

CN103067022A公开了《一种整型数据无损压缩方法、解压缩方法及装置》，提供一种整型数据压缩方法，包括：位重组处理过程、数据块划分及基本信息存储过程和编码处理过程；整型数据解压缩方法包括：基本信息解析过程、解码处理过程和位恢复处理过程；整型数据压缩装置包括：位重组处理模块、数据块划分及基本信息存储模块和编码处理模块；整型数据解压缩装置包括：基本信息解析模块、解码处理模块和位恢复处理模块。该方法具有算法简单、易于实现，编码解码效率高，压缩效果好，可用于对音频数据、地震勘探信号、地球物理测井波形数据和图像数据等多种类型的整型数据进行无损压缩和解压缩的优点。

2011年于文茂公开了《基于SPIHT改进算法的地震数据压缩》，用该算法对地震数据进行有效的压缩和重构，在减少编解码时间的同时，能够得到较高的压缩比和较好的图像重构效果。为了验证算法的有效性，选取了10幅地震剖面图进行试验，结果表明该算法能够很好地保留图像的客观质量，同时也提高了编解码的效率。2015年Carlos公开了《Seismic Data Compression Using 2D Lifting-Wavelet Algorithms》，在研究二维提升小波压缩地震数据过程中，通过实验发现分解层次越多信噪比越低，非均匀量化在大幅值时容易出错，在地震数据应用中用均匀量化来获得最小熵，当目标是实现信噪比高于40dB时，Huffman编码比算术编码更好。验证了压缩效果与滤波器类型和长度、分解层数、量化方式以及编码方式有关。

虽然上述现有技术可用于地震数据压缩与重构，但常见压缩算法中压缩倍数是以损失信号精度为代价，数据压缩程度直接影响到重构效果，实时压缩与高精度重构的需求无法同时满足。并且鉴于压缩过程更注重实时性，重构过程更注重精确性，常见压缩重构算法是无法满足地震勘探***实时压缩，并且高精度重构出原始数据的需求。

发明内容

本发明目的就在于针对上述现有技术的不足，提供一种地震数据压缩与重构方法。

本发明首先给出了基于压缩感知理论的地震数据压缩方法；然后，给出了地震数据重构方法。根据压缩感知理论可知：若信号是稀疏的或者是在某个变换域内稀疏，那么信号可以由测量矩阵将原始信号变换为低维信号，再通过重构算法优化求解得到原始信号。由此可以看出：压缩方法即测量矩阵的构造，重构方法即重构优化算法的设计。本发明与常见地震数据压缩重构算法的区别在于：可对地震数据实现边采集边压缩，压缩所用的测量矩阵易于硬件实现，并且压缩过程和重构过程是两个独立的算法实现，重构精度不完全受压缩比大小的制约，地震数据可以在下位机压缩，回收压缩后的数据，当后续需要地震资料解释时，再在上位机高精度重构出完整原始数据。

本发明的目的是通过以下技术方案实现的：

一种地震数据压缩与重构方法，包括以下步骤：

A、输入地震监测实际数据x；

B、对输入数据用小波基Ψ稀疏表示，得到稀疏系数θ＝Ψx；

C、构造混沌序列测量矩阵Φ；

D、利用测量矩阵对稀疏系数观测，得到压缩后的数据y＝Φθ+n，n表示噪声；

E、设计改进的贝叶斯小波树结构压缩感知重构算法(BTSWCS-vb)；

F、用重构算法求解完整数据的稀疏系数

G、对求得的稀疏系数反变换，得到完整的微震数据

H、与其他算法作比较，分析试验结果，对测量矩阵和重构算法进行评价，评价标准是压缩时间(Time)、峰值信噪比(PSNR)和均方根误差(RMSE)。

与现有技术相比，本发明的有益效果在于：

本发明基于压缩感知理论，构造混沌序列测量矩阵对地震数据的变换系数压缩，满足地震数据实时压缩的需求，并且算法简单，易于硬件实现。本发明运用改进的贝叶斯小波树结构压缩感知重构算法高精度重构出完整原始数据，压缩过程与重构过程之间联系较少，有效地解决了常规压缩方式不能同时满足低压缩比压缩和高精度重构的问题，并且算法相对于改进之前收敛速度更快，计算效率更高。按照本发明对地震数据进行压缩时，其压缩比在0.2～0.55之间，不仅极大地缓解了地震数据传输速度和存储压力，而且减小了地震数据传输和监测成本等方面的负担。

附图说明

图1为地震数据压缩与重构方法流程图；

图2为小波系数空间方向树父子关系示意图；

图3为混沌序列测量矩阵构造流程图；

图4a为几种混沌序列测量矩阵性能比较——PSNR折线图；

图4b为几种混沌序列测量矩阵性能比较——RMSE折线图；

图5为改进的贝叶斯小波树结构压缩感知重构算法先验模型；

图6a为BTSWCS-vb算法与BTSWCS-mcmc性能比较——Time折线图；

图6b为BTSWCS-vb、BTSWCS-mcmc和CoSaMP算法比较——PSNR折线图；

图6c为BTSWCS-vb、BTSWCS-mcmc和CoSaMP算法比较——RMSE折线图；

图7a为试验所用的实际地震数据图；

图7b为用Logistic矩阵做0.25的观测后，再用BTSWCS-vb算法恢复出的地震数据。

具体实施方式

下面结合附图对本发明进行详细说明。

如图1所示，一种地震数据压缩与重构方法，包括以下步骤：

A、输入地震监测实际数据x；

B、对输入数据用小波基Ψ稀疏表示，得到稀疏系数θ＝Ψx，

稀疏表示过程：θ＝Ψx表示，实际地震数据x经过小波变换后得到稀疏系数θ；

C、构造混沌序列测量矩阵Φ；

矩阵构造：混沌序列测量矩阵的构造是地震数据压缩的关键，根据本发明测量矩阵的构造步骤，选择合适的混沌序列和设置合理的初始值，来满足测量矩阵所需的RIP条件。

D、利用测量矩阵对稀疏系数观测，得到压缩后的数据y＝Φθ+n，n表示噪声，

压缩过程：y＝Φθ+n表示，在成功构造矩阵的基础上，用列数远大于行数的测量矩阵与稀疏系数乘积，从而达到降维的效果，得到压缩后的值y；

E、设计改进的贝叶斯小波树结构压缩感知重构算法(BTSWCS-vb)；

算法改进：原有的贝叶斯小波树结构算法重构效果远优于常见贪婪类算法，不足之处在于重构所用时间相对要长；本发明改进后的算法克服了这一缺点，明显的缩短了重构时间。

F、用重构算法求解完整数据的稀疏系数

重构过程：用改进的贝叶斯小波树结构重构算法求解出完整实际地震数据的稀疏系数；

G、对求得的稀疏系数反变换，得到完整的微震数据

H、与其他算法作比较，分析试验结果，对测量矩阵和重构算法进行评价，评价标准是压缩时间(Time)、峰值信噪比(PSNR)和均方根误差(RMSE)。

上述步骤中有三个关键技术：小波变换、构造测量矩阵、改进重构算法，下面对它们的实现原理介绍。

1、本发明地震信号小波变换原理：

由于地震信号在时间域上是不稀疏的，根据压缩感知理论，需要对它进行稀疏表示。根据地震信号分频带性强、幅度变化大、在二维空间上相关性高等特性，且小波基与地震信号具有相似性，所以对地震信号进行小波变换，得到的小波系数具有很好的可压缩性。对一维信号(单道微震信号)的离散小波多分辨率分析，Mallet算法描述公式为：

式中，C_j,k、D_j,k分别是原始信号在尺度2^-j下分解得到的低频分量和高频分量，h(m-2k)、g(m-2k)分别是低通和高通滤波器的系数。相应的重建公式如下：

对二维地震信号进行分解与重构可以分别对行和列进行一维小波变换得到，其中行号和列号分别对应地震剖面的时间和道号。本发明选择“db1”小波基对地震信号稀疏表示。

小波变换除了增加地震信号的可压缩性外，如图2小波系数空间方向树父子关系所示，小波系数还有结构统计特性。地震信号经过离散小波变换后生成四叉树结构的小波系数，即每个小波系数作为“父节点”通常有四个“子节点”，基于贝叶斯压缩感知的重构算法就利用了小波系数这一结构特性。我们用箭头方向来表明树中不同尺度上系数之间的关系，箭头是由父节点指向子节点。尺度S＝1对应的系数叫父节点，最大尺度S＝L(L＝3)上的系数叫子节点。在1≤S≤L-1范围类的小波系数都有尺度S+1上有四个“儿子”；只有低频区域(S＝0)中每一组(2×2)矩阵中左上方那个系数和高频区域(S＝3)中的所有系数没有“儿子”。

综上，地震信号小波变换后的系数能量很集中，多数数据是趋近于0的，信号中可以忽略的小波系数往往聚集在一起；另一方面，父子节点有一定的关系，比如，如果一个系数在一个尺度上是可以忽略的，那他的子节点通常也是可以忽略的。

2、本发明混沌序列测量矩阵的构造原理：

压缩感知理论是通过极少的数据来获取尽可能多的信息以精确恢复原始信号，因此测量矩阵的作用至关重要。虽然目前所普遍采用的高斯、伯努利等随机测量矩阵的重构结果比较好，但在仿真试验中，测量矩阵的不确定性直接对实验结果的鲁棒性产生影响；而在实际工程应用中，还会造成计算复杂度高、存储空间大、硬件不易实现等问题。

混沌被认为是具有随机性行为的确定性动力***的现象。由于混沌序列自身良好的伪随机性特点，并且构造混沌序列测量矩阵的复杂度远远低于随机测量矩阵，所以其工程实现意义重大。如图3所示，本发明构造混沌序列测量矩阵步骤为：

步骤1，根据相关混沌映射的变换核，选取合适的序列初始值和***参数使混沌***达到理想混沌状态来产生混沌序列{u_n}，序列长度为M×N-1。

步骤2，将步骤1生成的混沌序列{u_n}通过符号函数映射成序列{a_n}。

步骤3，将步骤2生成的序列{a_n}取N长截断形成M×N维测量矩阵Φ。

本发明使用如下所述的几种混沌序列构造测量矩阵：

第一种：Logistic(罗切斯特)映射变换核：x_n+1＝ax_n(1-x_n) (3)

初始值设定：x₀＝0.37，a＝2；符号函数：

第二种：Henon(埃农)映射变换核：

初始值设定：x₀＝y₀＝0，a＝1.4，b＝0.314；符号函数：

第三种：Tentmap(帐篷)映射变换核：x_n+1＝a-(1+a)|x_n|， (5)

初始值设定：x₀＝0.01，a＝0.99；符号函数：

第四种：Kent(肯特)映射变换核：

初始值设定：x₀＝0.36，a＝0.8；符号函数：

从构造过程中符号函数的使用，可以看出：相比于常用随机测量矩阵，本发明的混沌测量矩阵是便于硬件实现的，因为矩阵的元素只有0、1和-1。为评价几种混沌测量矩阵的性能，用不同的测量矩阵对地震数据进行不同程度的压缩，然后统一用CoSaMP算法重构，结果如图4a所示，几种混沌测量矩阵对重构结果PSNR值的影响，图4b是对重构结果RMSE值的影响。可以看出，按本发明的矩阵构造方法和初始值、参数的选取，相对来说Logistic矩阵效果最好。

3、本发明改进的贝叶斯小波树结构压缩感知重构算法(BTSWCS-vb)设计原理：

为克服低信噪比下，常用贪婪算法重构性能差的问题，本发明考虑的是贝叶斯框架下的压缩感知重构算法，贝叶斯小波树结构压缩感知算法(BTSWCS)。该算法首先假设贝叶斯先验模型稀疏系数θ中每个元素都服从钉铆先验分布π_i表示权重，δ₀部分表示为0的系数，然后用马尔可夫链蒙特卡罗(Markov Chain Monte Carlo，MCMC)推理对先验模型参数后验估计，但是，采用MCMC推理所需的重构时间长，收敛速度慢，计算效率低。

变分贝叶斯(Variational Bayesian，VB)推理近年来在贝叶斯分析中被广泛应用，因为它不仅可以提供一个精确的后验，而且VB推理收敛速度快，计算效率高。本发明所用的改进BTSWCS算法后验估计采用的就是VB推理，与MCMC推理相比迭代相同的次数，重构效果明显提高，重构时间明显减少。

本发明为了便于使用VB推理，提出了一个与钉铆先验等价的先验模型，模型中没有δ₀组件部分。先验模型假设为 符号表示Hadamard乘积，ω表示不为零的稀疏系数，z是指示参数，为0或1，z_i～Bernoulli(π_i)；稀疏系数不为0时，组件z＝1；反之，组件z＝0。模型可以归纳为：

α₀～Gamma(a₀,b₀) (9a)

α_s～Gamma(c₀,d₀),0≤s≤L (9b)

下标(s,i)表示小波变换尺度s上的第i个元素，pa(s,i)表示相关的父节点。模型中π_s,i应用了小波树结构特性，不同的父节点值pa(s,i)，就有不同的指示参数z，先验参数π也就不同。当系数的父节点值为零时，π取π^s0；当系数的父节点值不为零时，π取π^s1。模型的图形形式参考图5。VB推理使用一组分布q(θ)来近似估计真实的后验分布p(θv)，运用下限F以接近模型的真实对数似然函数logp(vθ)。算法通过迭代更新q(θ)，来使F趋近于logp(vθ)，直到收敛。

为评价改进重构算法性能，首先对实际数据统一用Tentmap矩阵进行不同压缩比的观测，再分别用三种不同算法重构，实验结果取30次实验的平均值，对重构结果进行比较。图6a比较的是BTSWCS-vb算法和BTSWCS-mcmc算法重构所以时间，很容易看出：BTSWCS-vb算法所用时间更短，计算速度更快；图6b是BTSWCS-vb算法、BTSWCS-mcmc算法和CoSaMP算法重构结果PSNR值比较，图6c是三种算法重构结果RMSE值比较，很直观的看出：BTSWCS-vb算法重构峰值信噪比最高，重构误差最小，重构效果最好。本发明试验所用实际数据如图7a所示，用Logistic矩阵做压缩比为0.25的观测后，再用BTSWCS-vb算法恢复出完整数据，结果如图7b所示，重构精度很高，误差只有0.0058。综上所述，从重构所用时间、重构结果PSNR值和重构结果RMSE值来看，改进算法BTSWCS-vb优势是显而易见的。

Claims (1)

1.一种地震数据压缩与重构方法，其特征在于，包括以下步骤：

A、输入地震监测实际数据x；

B、对输入数据用小波基Ψ稀疏表示，得到稀疏系数θ＝Ψx；

C、构造混沌序列测量矩阵Φ；

D、利用测量矩阵对稀疏系数观测，得到压缩后的数据y＝Φθ+n，n表示噪声；

E、设计改进的贝叶斯小波树结构压缩感知重构算法(BTSWCS-vb)；

F、用重构算法求解完整数据的稀疏系数

G、对求得的稀疏系数反变换，得到完整的微震数据

H、与其他算法作比较，分析试验结果，对测量矩阵和重构算法进行评价，评价标准是压缩时间(Time)、峰值信噪比(PSNR)和均方根误差(RMSE)。