基于压缩感知的小波域地震数据实时压缩与高精度重构方法

摘要

本发明涉及一种基于压缩感知的小波域地震数据实时压缩与高精度重构方法，包括以下步骤：首先对微震信号在小波域稀疏表示；然后利用Logistic混沌序列构造混沌伯努利测量矩阵(CBMM)，并用测量矩阵对稀疏表示的微震信号压缩观测；最后采用贝叶斯小波树结构压缩感知重构算法(BTSWCS)，恢复出完整原始数据。本发明经实际对比实验结果表明：使用本方法对总采样点为28的数据压缩，压缩时间可缩短至10‑5s，也就是说，若地震仪采样率是1KSPS，CBMM测量矩阵对它采集的0.25s的数据基本可以实现实时压缩。低信噪比情况下，本发明重构算法使PSNR值至少提升5dB，相对贪婪算法，明显提高了峰值信噪比，降低了重构误差。

说明书

技术领域

本发明涉及一种微地震数据的处理方法，特别涉及一种基于压缩感知的小波域地震数据实时压缩与高精度重构方法。

背景技术

数据实时传输是制约无缆存储式地震仪发展的最主要因素。随着地震勘探的不断深入，地震勘探正向多维、多分量和高分辨等方向发展，使得地震勘探数据持续膨胀，给无缆地震仪的采样和无线传输速度、存储器的存储容量以及计算机的处理速度造成极大压力。通过对地震数据进行合理的压缩处理，可以提高实时处理速度，改善地震仪无线通讯数据传输性能。因此，缩减压缩时间、提高重构精度成为迫切需要研究的课题。

CN104749633A公开了一种《高效高精度地震勘探数据压缩方法》对地球物理勘探仪器中采集的信号经24位模数转换后的数据压缩1/3的数据量，把24位3字节数据压缩为16位2字节数据。首先，设24位模数转换器输出的数据为D24；接着，设经压缩后的数据格式为4位阶码E4+16位尾数M16；然后把24位数据在保留符号位的情况下转换为4位阶码E4加16位尾数M16，阶码E4的取值为0～8；最后把4个采样点数据组合成10个字节，即4个样点的阶码占用2个字节，4个尾数占用8个字节。该方法具有压缩误差小，能够保持数据原貌，实现高效高精度压缩的优点。

CN104378118A公开了一种《高效自适应地震数据流无损压缩及解压缩方法》，通过使用编码方式将采样数据原始的24位3字节形式自适应的压缩为1字节或2字节或3字节或者和少量数据转化为4字节；在数据压缩前，首先根据原始数据的数值大小进行所需字节的判定，0～63及-64～-1区间的数据压缩后为1字节，64～8191及-8192～-65区间的数据压缩后为2字节，8192～104875及-104876～-8193区间的数据压缩方法操作后为3字节，占用字节与原来相同，除以上数据范围外，其他24位有符号整型数据可以表示的其他整型数在经压缩算法操作后需要使用4个字节表示。该方法具有大量节省存储空间，显著提高数据传输效率的优点。

2011年石水勇公开了《基于提升小波变换的地震勘探数据有损压缩研究》，以小波变换为基础，研究了不同小波转化为提升方案的原理，用于地震勘探数据提升小波正逆变换。对基于小波变换的三种典型的有损压缩方法一EZW、SPIHT、EBCOT方法，研究了其各自的原理和优缺点，从而选择出SPIHT算法作为地震勘探数据的编码方法。对于地震勘探原始数据和重构数据，设计并使用5种不同的评价标准对重构数据进行分析。2015年徐锋涛公开了《基于EZW算法的地震数据压缩》，对地震信号进行二维小波变换；利用系数之间的相关性对变换后的系数采用EZW方法重新组织；采用自适应算术编码对量化后的系数进行无损熵编码。

虽然上述现有技术可用于微地震信号压缩重构，但由于地球物理数据的特殊性，数据采集时要求尽量保持数据原貌，但这些数据压缩方式均不能用于数据采集阶段，必须对完整的地震数据进行操作，而不能对一个或少量几个采样数据处理，无法对数据流实时压缩，并且解码纯粹是编码的逆运算，压缩倍数是以损失信号精度为代价，无法做到降低压缩比的同时，减少重构误差。

发明内容

本发明目的在于针对上述现有技术的不足，提供一种基于压缩感知的小波域地震数据实时压缩与高精度重构方法。

由压缩感知理论可知：若信号是稀疏的或者是在某个变换域内稀疏，那么信号可以由测量矩阵将原始信号变换为低维信号，再通过重构算法优化求解得到原始信号。本发明首先，在小波域对微震数据稀疏表示；然后，构造测量矩阵对稀疏系数边采集边压缩；最后，设计重构算法高精度恢复出完整数据。本发明与常规压缩编解码算法的区别在于：不是先采集后压缩，而是边采集变压缩；解码不受编码的制约，而是独立的重构算法恢复出采集的完整数据。

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

一种基于压缩感知的小波域地震数据实时压缩与高精度重构方法，包括以下步骤：

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

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

C、构造混沌伯努利测量矩阵(CBMM)Φ对数据进行压缩；

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

E、设计贝叶斯小波树结构压缩感知重构算法(BTSWCS)恢复信号；

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

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

H、与其他算法作比较，分析试验结果，对测量矩阵和重构算法进行评价。

优选地，步骤B，所述小波基Ψ为bior4.4小波基。

与现有技术相比，本发明的有益效果在于：本发明基于压缩感知理论，采用构造测量矩阵的方式，对稀疏表示后的微震信号边采样边压缩，有效地解决了常规压缩方式需要先采集再压缩的问题。同时，通过设计重构算法对观测数据高精度恢复重构，且具有良好的抗噪性能。通过对比实验可以看出，本发明实现实时压缩，对低信噪比数据也能高精度重构，可以广泛地应用于地球物理勘探仪器和其他场合，有利于改善地震仪无线通讯数据传输性能、降低地球物理仪器野外部件的数据存储空间、提高数据传输效率的优点。

附图说明

图1为基于压缩感知的小波域地震数据实时压缩与高精度重构流程图；

图2为中国西北部某区的地震勘探实际数据图；

图3a为当数据长度为2¹⁴时，分别用SPIHT编码算法和CBMM压缩所有时间柱状图；

图3b为当数据长度为2¹²时，分别用SPIHT编码算法和CBMM压缩所有时间柱状图；

图3c为当数据长度为2¹⁰时，分别用SPIHT编码算法和CBMM测量矩阵压缩所有时间柱状图；

图3d为当数据长度为2⁸时，分别用SPIHT编码算法和CBMM测量矩阵压缩所有时间柱状图；

图4为BTSWCS算法精确重构所需的MCMC迭代次数曲线图；

图5a为压缩比0.2时，贪婪算法重构地震记录图；

图5b为压缩比0.2时，BTSWCS算法重构地震记录图；

图6a为实际数据小波变换产生的小波系数图；

图6b为BTSWCS算法重构出的小波系数图；

图7a1为实际数据第2道记录的时域波形图；

图7a2为BTSWCS算重构的实际数据第2道记录的时域波形图；

图7a3为实际数据第2道记录的频谱图；

图7a4为BTSWCS算重构的实际数据第2道记录的频谱图；

图7b1为实际数据第15道记录的时域波形图；

图7b2为BTSWCS算重构的实际数据第15道记录的时域波形图；

图7b3为实际数据第15道记录的频谱图；

图7b4为BTSWCS算重构的实际数据第15道记录的频谱图；

图8a为贪婪算法、SPIHT解码算法和BTSWCS算法重构PSNR折线图；

图8b为贪婪算法、SPIHT解码算法和BTSWCS算法重构RMSE折线图。

具体实施方式

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

如图1所示，一种基于压缩感知的小波域地震数据实时压缩与高精度重构方法，包括以下步骤：

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

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

C、构造混沌伯努利测量矩阵(CBMM)Φ；

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

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

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

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

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

上述步骤中有三个关键技术：稀疏表示、构造测量矩阵、设计重构算法，下面对它们的实现原理介绍。

1、本发明微震数据用小波基稀疏表示原理：

微震数据压缩方面小波分析是其理论基础，小波理论提供的自适应时频分析法能自动调节时-频窗，以适应实际分析的需要。由于地震信号通常表现出分频带性强、幅度变化大、在二维空间上相关性高等特性，且小波基与地震信号具有相似性，所以对地震信号进行小波变换后，得到的小波系数具有很好的稀疏性(可压缩性)。对一维信号(单道微震信号)的离散小波多分辨率分析，Mallet算法描述公式为：

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

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

2、本发明混沌伯努利测量矩阵(CBMM)的构造原理：

已知Logistic混沌系统如式(3)所示：

该Logistic系统在参数λ∈[1.872,2.0]区间的条件下，当初值u₀＝0.23,0.37,0.7时，其Lyapunov指数均大于0，此时的Logistic系统是混沌系统。

由该Logistic系统产生的随机序列{u_n}经过式(4)符号函数映射成序列{a_n}，其中：

当λ＝2.0时，由该Logistic系统产生的序列{a_n}满足Bernoulli分布，同时也满足RIP性质，则由Logistic系统产生的序列{a_n}可作为CS的测量矩阵。

本发明构造测量矩阵算法步骤为：

步骤1，根据Logistic映射扩频序列的相关分布，取λ＝2.0，u₀＝0.37，通过该Logistic混沌系统来产生混沌序列{u_n}，序列长度为M×N-1。

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

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

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

为克服低信噪比下，常用贪婪算法重构性能差的问题，本发明考虑的是贝叶斯框架下的压缩感知恢复算法，利用参数的先验概率和样本信息所满足的后验概率，继而获得总体的概率密度函数。

BTSWCS算法步骤为：

步骤1，假设贝叶斯先验模型，设定稀疏系数θ中每个元素都服从钉铆先验分布

步骤2，对先验模型参数后验估计，本发明用MCMC推理，该方法是通过足够的样本数来估计后验分布，通过迭代每个随机变量(模型参数和中间变量)来收集这些样本，迭代的随机变量是从它的条件后验分布给出的所有其他随机变量的最新值。

p(α_n|-)＝Gamma(a₀+N/2,b₀+(y-Φθ)^T(y-Φθ)/2)(14)

本发明实验所用实际数据如图2所示，为证明用CBMM算法压缩要比常用SPIHT编码算法更具时效性，选用长度分别为2¹⁴、2¹²、2¹⁰和2⁸的数据做压缩比为0.2、0.25、0.3、0.35、0.4、0.45、0.5和0.55的处理，记录CBMM算法和SPIHT编码算法完成压缩所用时间。

因为测量矩阵是由随机数组成的，所以结果取50次试验的平均值，如图3a、3b、3c和3d所示。从柱状对比图可以清晰看出：压缩相同数据，达到相同压缩比情况下，CBMM算法要比SPIHT编码算法所用时间更短。特别是，随着待压缩数据量的减少，CBMM压缩的优势更加突出，对长度为2⁸的数据用CBMM矩阵观测，所用压缩时间更是缩短至10^-5s级；也就是说，若地震仪采样率是1K，CBMM测量矩阵对它采集的0.25s的数据基本可以实现实时压缩。

下位机通过CBMM对地震数据边采集变压缩，上位机通过BTSWCS重构算法恢复出完整记录。本发明实验选用的实际数据(图2)总长2¹⁴，在压缩比分别为0.2、0.25、0.3、0.35、0.4、0.45、0.5和0.55的情况下，BTSWCS算法精确重构所需的MCMC迭代次数曲线如图4所示。为了验证BTSWCS重构算法在低信噪比下的优越性，首先对图2所示的实际数据用CBMM进行压缩比为0.2的压缩，然后用贪婪算法恢复(OMP)，结果如图5a所示。实际数据在压缩比为0.2的情况下，根据图4，BTSWCS算法200次迭代用重构出完整数据，结果如图5b所示。可以很直观的看出，用BTSWCS算法重构效果明显比贪婪算法好，婪算法对本身含有噪声和含有测量噪声的数据进行重构处理时，重构性能明显下降，不能达到重构精度要求。图6和图7是对图5b重构结果的细节进行评价。图6a是实际数据小波系数，图6b是重构出的小波系数；图7a1是实际数据第2道记录的时域波形，图7a2是重构的第2道记录的时域波形，图7a3是实际数据第2道记录的频谱图，图7a4是重构的第2道记录的频谱图；图7b1是实际数据第15道记录的时域波形，图7b2是重构的第15道记录的时域波形，图7b3是实际数据第15道记录的频谱图，图7b4是重构的第15道记录的频谱图。

用贪婪算法、SPIHT解码和BTSWCS算法对压缩的实际数据恢复重构，结果取30次实验平均值，三种算法重构性能比较如图8所示。图8a是三类算法重构实际数据时的PSNR对比曲线，对实际数据的处理，贪婪算法整体峰值信噪比小于65dB，SPIHT解码和BTSWCS算法重构的PSNR值均大于贪婪算法重构。尤其是BTSWCS算法，它的PSNR至少比贪婪算法高出5dB，随着压缩比或迭代次数的增加，PSNR值会更高。图8b是三类算法重构实际数据时的RMSE对比曲线，贪婪算法在重构实际数据时，重构误差均大于0.1，重构误差很大，SPIHT解码和BTSWCS算法要比贪婪算法重构误差小很多，特别是BTSWCS，随着压缩比或迭代次数的增加，误差基本趋近于0。