CN107664773A - 一种基于时间移位和熵约束的地震信号编码方法 - Google Patents

Abstract

本发明提供一种基于时间移位和熵约束的地震信号编码方法，解决在有限信道下如何有效传输高质量地震信号问题。本发明包括：将编码过程分为离线训练和在线测试两个阶段:（1）在离线训练阶段，结合离线地震信号数据，通过迹对齐和编码参数训练两个步骤生成在线测试编码所需的字典和概率密度函数；（2）在线测试阶段，通过迹对齐和熵约束稀疏编码生成时间移位参数、低频系数、高频系数、低频系数位置和非零系数位置等信息传输到解码端，其中低频系数和高频系数通过熵编码后再进行传输。在解码端，首先通过熵解码算法恢复低频系数和高频系数，并进一步结合相关位置信息进行地震信号重建，最后通过时间移位恢复得到最终的重构地震信号数据。

Description

一种基于时间移位和熵约束的地震信号编码方法

技术领域

本发明属于地震信号数据传输方法，具体涉及一种基于时间移位和熵约束的地震信号编码方法。

背景技术

世界石油工业发展史和中国油气勘探开发的实践证明，地震技术是寻找地底矿产资源的有效方式之一。通过在地表均匀布置震源及传感器，再结合人工放炮方式制造人工“地震”，由传感器接收地底返回信号从而对地底资源分布情况进行判断。为了提高获得的图像质量，采用高密度传感器组网进行探测是未来的发展趋势。随之而来的问题是采集信号的数据量将会急剧增长，而目前的传输带宽是非常有限的，因此有必要研究有效的地震信号编码技术来减少高密度传感网下地震信号传输的数据量。

现有技术中，文献[1](“Transform methods for seismic data compression,”IEEE Transactions on Geoscience and Remote sensing)通过基于离散余弦的地震信号编码方法，取得了接近于3倍的压缩倍数。为了在压缩的同时，能够对重要特征信息进行保持，可以采用二维自适应局部离散余弦变换方法。(“Adaptive local cosine transformfor seismic image compression,”2006International Conference on AdvancedComputing and Communications)。在地震信号编码中，一个核心技术是需要对地震信号进行稀疏表示。除了上述方法中的离散余弦变换或者自适应离散余弦变换之外，其他方法还包括小波变换，Curvelet变换，Sieslet变换等。这些方法往往通过固定基的形式对信号进行表征，近年来，字典学习在数据编码、去噪等领域得到了广泛的应用。例如，在图像编码中(“Compressibility constrained sparse representation with learnt dictionaryfor low bit-rate image compression,”IEEE Transactions on Circuits and Systemsfor Video Technology)，通过字典学习产生适应于不同信号表征的字典，从而使得信号的表征更为稀疏，因此获得更好的编码性能。因此，相关研究结果证明了在地震信号编码中通过字典学习进行信号表征是有效的。

传统字典学习往往是针对信号的表征，要求其表征尽可能的稀疏。而对于编码应用而言，系数稀疏并不等同于系数经过编码后的码流尽可能的少。熵是衡量编码性能的关键要素之一，因此，有必要将熵引入到字典学习过程中。同时针对地震信号不同迹存在时间移位等现象，在字典学习过程中还可以引入时间移位约束。

发明内容

本发明提出一种基于时间移位和熵约束的地震信号编码方法，属于地震信号编码、传输及采集方法，解决在高密度传感网络中有限信道下地震信号编码传输问题，可以应用于基于无线/有线传感网的地底资源探测中。

本发明的一种基于时间移位和熵约束的地震信号编码方法，将编码过程分为离线训练和在线测试两个阶段，步骤为：

S1.在离线训练阶段，结合离线地震信号数据，通过迹对齐和编码参数训练两个步骤生成在线测试编码所需的字典和概率密度函数；

S2.在在线测试阶段，结合测试地震信号数据，通过迹对齐和熵约束稀疏编码生成时间移位参数、低频系数、高频系数、低频系数位置和非零系数位置等信息传输到解码端，其中低频系数和高频系数通过熵编码后再进行传输；

在解码端包括S3：

S3.首先通过熵解码算法恢复低频系数和高频系数，并进一步结合低频系数位置和非零系数位置通过稀疏重建算法进行地震信号重建，最后通过时间移位恢复(采用时间移位参数)得到最终的重构地震信号数据。

更进一步地，步骤S1具体为：

S11：对离线地震信号数据进行迹对齐，包括时间移位稀疏编码和字典更新

两个步骤，从而生成时间移位参数和重新对齐数据，其具体算法如下：

输入参数：离线地震信号Y＝[y₁,...,y_K]，最大迭代次数I，初始字典D₀，最大时间偏移参数S，稀疏性L，每组里面段的数目N。

初始化：将输入信号Y划分为多个组，从而可以表示为Y＝[P₁,...P_b,...P_K/N]，其中P_b＝[y_(b-1)N+1,...y_bN],b∈[1,K/N]，循环初始化参数k＝0，b＝0，j＝-S-1，字典初始化D⁽⁰⁾＝D₀

主循环：令k＝k+1，如果k≤I，在第一循环内做以下运算：

令b＝b+1，如果在第二循环内做以下运算：

令j＝j+1，如果j≤S，在第三循环内做以下运算：

S111：对P_b沿着迹的方向进行整体平移，移动量为j，得到平移后的组P_b ^j；

S112：基于字典D^(k-1)，采用ORMP算法(“Comparison of basis selectionmethods,”in Signals,Systems and Computers,1996.Conference Record of theThirtieth Asilomar Conference)对P_b ^j进行稀疏表示，得到稀疏系数W_b ^j；

S113：计算残差e^j，

第三循环结束

S114：计算残差e^j(j∈[-S,S])的最小值，并将最小值对应的j赋值给变量T^(k) _b；

S115：将P_b整体移动T_b，并用移动后的结果更新P^(k) _b；

第二循环结束

S116：基于所有时间偏移后的地震信号数据，通过MOD方法(“Method of OptimalDirections for Frame Design,”in 1999IEEE International Conference onAcoustics,Speech and Signal Processing(ICASSP))训练并更新字典，生成字典D^(k)；

第一循环结束

输出：重新对齐的地震信号数据P_b＝P^(I) _b(b∈[1,K/N])，时间偏移参数T_b＝T^(I) _b

S12：对编码参数进行训练，生成概率密度函数和字典，包括熵约束稀疏编码、字典更新及概率密度函数更新三个步骤，具体流程如下：

输入参数：重新对齐的地震信号最大迭代次数I，初始字典D₀，初始概率密度函数χ₀，阈值τ，稀疏性L，量化步长Δ

初始化：循环初始化参数k＝0，i＝0，j＝0，α＝-E-1，字典初始化D⁽⁰⁾＝D₀，概率密度函数χ⁽⁰⁾＝χ₀

主循环：令k＝k+1，如果k≤I，在第一循环内做以下运算：

令i＝i+1，如果i≤K，在第二循环内做以下运算：

S121:基于字典D^(k-1)，对重新对齐的地震信号采用ORMP方法进行稀疏表示，得到稀疏系数w_i；

S122:对稀疏系数w_i进行量化得到将中的非零值表示为{z₁,...,z_L}；

S123:初始化参数γ_j＝0,j∈[1,L]；

令j＝j+1，如果j≤L，在第三循环内做以下运算：

如果|z_j|≤τ，在第四循环内做以下运算：

S124:令α＝α+1，如果α≤E，计算

其中，En(z_j+α)＝-plog₂p，p是z_j+α在χ^(k-1)中的概率；

第四循环结束

S125:找到TE_α的最小值对应的值更新z_j:＝z_j+α_j

第三循环结束

S126:对对应的非零值{z₁,...,z_L}进行更新；

第二循环结束

S127:采用MOD方法更新字典D^(k)；

S128:统计每个数值出现的概率，从而对概率密度函数χ^(k-1)进行更新；

第一循环结束

输出：概率密度函数χ＝χ^(I)和字典D＝D^(I)。

更一步地，步骤S2具体为：

S21：对测试地震信号数据进行迹对齐，采用与S11相同的步骤，从而生成时间移位参数和重新对齐数据；

S22:基于步骤S12已经生成的概率密度函数χ和字典D，对重新对齐后的地震信号数据进行熵约束稀疏编码，具体步骤包含步骤S12中的步骤S121至S126(无需步骤S127和步骤S128)，生成稀疏系数

S23:根据阈值τ，将稀疏系数中大于阈值的划分为低频系数，小于或者等于阈值且不为0的划分为高频系数；

S24:将稀疏系数中非零系数位置，低频系数位置进行传输，也可以对其进行算术编码后传输；

S25:将时间偏移参数进行传输；

S26:对低频系数进行差分编码，对高频系数基于已有的概率密度函数进行自适应算术编码；

更一步地，步骤S3具体为：

S31：对低频和高频编码码流分别进行差分解码和自适应算术解码(基于已有的概率密度函数χ)，生成低频系数和高频系数；

S32：基于非零系数位置和低频系数位置，重建稀疏系数

S33：基于已有的字典D和重建稀疏系数生成重建信号

S34：根据所传输的时间偏移参数，对重建信号进行时间移位恢复(逆向偏移)，得到最终的重建地震信号数据。

与现有技术相比，本发明具有以下优点和有益效果：在相同压缩倍数下，本发明可以获得更佳的重建地震信号质量，从而适用于各种在有线/无线信道下对高质量地震信号实时传输的应用场合，例如基于高密度传感器组网的地底矿产资源探测系统。

附图说明

图1是本发明的流程图；

图2是仅仅通过时间移位进行字典学习所生成的字典；

图3是时间移位字典与MOD算法生成字典间的差异；

图4是时间移位字典学习算法和MOD+ORMP字典学习算法在不同稀疏度下重建地震信号质量的对比情况；

图5是不同方法的重建地震信号质量对比结果：图5(a)为原始信号，图5(b)为MOD+ORMP算法的重建结果，图5(c)为基于熵约束算法的重建结果，图5(d)为基于时间移位和熵约束算法的重建结果；

图6是不同方法的重建信号与原始信号之间的差异。

具体实施方式

下面结合附图所示的实施例对本发明作进一步说明。

实施例1：

1、测试地震信号数据来源于2004BP地震信号数据库(http://software.seg.org/datasets/2D/2004_BP_Vel Benchmark/)；

2、将地震信号分为多个组分别进行测试，每组大小为960×20。对每组数据，将其分为多个段进行字典学习，每个段的大小为16×1，此时对应每组数据对应的离线地震信号可以表示为Y＝[y₁,...,y₁₂₀₀]；首先验证时间移位字典学习过程对整体性能的影响：

输入参数：Y＝[y₁,...,y₁₂₀₀]，最大迭代次数I＝10，初始字典D₀大小为16×60(基于离散余弦变换函数构造)，最大时间移位参数值S＝5，稀疏性L＝8，每组里面段的数目N＝2；

初始化：将输入信号Y划分为600组，从而可以表示为Y＝[P₁,...P₆₀₀]，其中P_b＝[y_2(b-1)+1,...y_2b],b∈[1,600]，循环初始化参数k＝0，b＝0，j＝-S-1，字典初始化D⁽⁰⁾＝D₀

主循环：令k＝k+1，如果k≤10，在第一循环内做以下运算：

令b＝b+1，如果b≤600，在第二循环内做以下运算：

令j＝j+1，如果j≤5，在第三循环内做以下运算：

S111：对P_b沿着迹的方向进行整体平移，移动量为j，得到平移后的组

S112：基于字典D^(k-1)，采用ORMP算法(“Comparison of basis selectionmethods,”in Signals,Systems and Computers,1996.Conference Record of theThirtieth Asilomar Conference)对P_b ^j进行稀疏表示，得到稀疏系数W_b ^j；

S113：计算残差e^j，

第三循环结束

S114：计算残差e^j(j∈[-S,S])的最小值，并将最小值对应的j赋值给变量T^(k) _b；

S115：将P_b整体移动T_b，并用移动后的结果更新P^(k) _b；

第二循环结束

S116：基于所有时间偏移后的地震信号数据，通过MOD方法(“Method of OptimalDirections for Frame Design,”in 1999IEEE International Conference onAcoustics,Speech and Signal Processing(ICASSP))训练并更新字典，生成字典D^(k)；

第一循环结束

输出：重新对齐的地震信号数据P_b＝P⁽¹⁰⁾ _b(b∈[1,600])，时间偏移参数T_b＝T⁽¹⁰⁾ _b

3、在此情况下生成的字典如图2所示。为了证明与传统方法有所差异，本实施例将仅仅通过MOD算法进行字典学习，所生成的字典与图2字典的差异在图3中进行显示；

4、进一步地改变上述实验中的稀疏度L(取值包括4、6、8)，生成相应的字典D和稀疏系数W并重建地震信号，图4给出了上述算法和MOD+ORMP字典学习算法在不同稀疏度下重建地震信号质量的对比情况。从图4可以看出，引入时间移位到字典学习中可以提高重建地震信号质量。

实施例2：

1、在本实施例中，选用一部分数据作为离线训练数据，已通过训练生成在线测试所需的字典D(字典大小为16×100)和概率密度函数χ。正则化参数λ＝8；

2、下面介绍如何在在线测试阶段，通过迹对齐和熵约束稀疏编码生成时间移位参数、低频系数、高频系数、低频系数位置和非零系数位置等信息。

输入参数：Y＝[y₁,...,y₁₂₀₀]，最大迭代次数I＝10，字典为D，最大时间移位参数值S＝5，稀疏性L＝8，每组里面段的数目N＝2；

初始化：将输入信号Y划分为600组，从而可以表示为Y＝[P₁,...P₆₀₀]，其中P_b＝[y_2(b-1)+1,...y_2b],b∈[1,600]，循环初始化参数k＝0，b＝0，j＝-S-1

主循环：令k＝k+1，如果k≤10，在第一循环内做以下运算：

令b＝b+1，如果b≤600，在第二循环内做以下运算：

令j＝j+1，如果j≤5，在第三循环内做以下运算：

S111：对P_b沿着迹的方向进行整体平移，移动量为j，得到平移后的组P_b ^j；

S112：基于字典D，采用ORMP算法(“Comparison of basis selection methods,”in Signals,Systems and Computers,1996.Conference Record of the ThirtiethAsilomar Conference)对P_b ^j进行稀疏表示，得到稀疏系数W_b ^j；

S113：计算残差e^j，

第三循环结束

S114：计算残差e^j(j∈[-S,S])的最小值，并将最小值对应的j赋值给变量T^(k) _b；

S115：将P_b整体移动T_b，并用移动后的结果更新P^(k) _b；

第二循环结束

第一循环结束

输出：重新对齐的地震信号数据P_b＝P⁽¹⁰⁾ _b(b∈[1,600])，时间偏移参数T_b＝T⁽¹⁰⁾ _b

将P_b按照段的方式排列成为了后续量化设计方便，所有数据均同时乘以2¹⁶；

输入参数：重新对齐的地震信号数据最大迭代次数30，字典为D，初始概率密度函数χ，阈值τ＝2¹²，稀疏性L＝8，量化步长Δ＝32

初始化：循环初始化参数k＝0，i＝0，j＝0，α＝-17，字典D，概率密度函χ

主循环：令k＝k+1，如果k≤30，在第一循环内做以下运算：

令i＝i+1，如果i≤1200，在第二循环内做以下运算：

S121:基于字典D，对重新对齐的地震信号采用ORMP方法进行稀疏表示，得到稀疏系数w_i；

S122:对稀疏系数w_i进行量化得到将中的非零值表示为{z₁,...,z₈}；

S123:初始化参数γ_j＝0,j∈[1,8]；

令j＝j+1，如果j≤8，在第三循环内做以下运算：

如果|z_j|≤τ，在第四循环内做以下运算：

S124:令α＝α+1，如果α≤16，计算

其中，En(z_j+α)＝-plog₂p，p是z_j+α在χ^(k-1)中的概率；

第四循环结束

S125:找到TE_α的最小值对应的值更新z_j:＝z_j+α_j

第三循环结束

S126:对对应的非零值{z₁,...,z_L}进行更新；

第二循环结束

第一循环结束

输出：中的非零值(大于τ的为低频系数，小于τ的为高频系数)，及对应的位置信息。

3、将稀疏系数中非零系数位置，低频系数位置进行传输，也可以对其进行算术编码后传输；

4、将时间偏移参数T_b进行传输；

5、对低频系数进行差分编码，对高频系数基于已有的概率密度函数进行自适应算术编码；

6、图5给出了不同方法的重建地震信号质量对比结果；

7、图6给出了不同方法的重建信号与原始信号之间的差异。可以看出，本发明从整体上而言具有最小的重建差异。

Claims (6)

1.一种基于时间移位和熵约束的地震信号编码方法，其特征在于，将编码过程分为离线训练和在线测试两个阶段，包括S1和S2:

S1.在离线训练阶段，结合离线地震信号数据，通过迹对齐和编码参数训练两个步骤生成在线测试编码所需的字典和概率密度函数；

S2.在在线测试阶段，结合测试地震信号数据，通过迹对齐和熵约束稀疏编码生成时间移位参数、低频系数、高频系数、低频系数位置和非零系数位置信息传输到解码端，其中低频系数和高频系数通过熵编码后再进行传输；

在解码端包括S3：

S3.首先通过熵解码算法恢复低频系数和高频系数，并进一步结合低频系数位置和非零系数位置通过稀疏重建算法进行地震信号重建，最后通过时间移位恢复得到最终的重构地震信号数据。

2.根据权利要求1所述的一种基于时间移位和熵约束的地震信号编码方法，其特征在于，所述步骤S1具体包括以下内容:

S11：对离线地震信号数据进行迹对齐，包括时间移位稀疏编码和字典更新两个步骤，从而生成时间移位参数和重新对齐数据；

S12：对编码参数进行训练，生成概率密度函数和字典，包括熵约束稀疏编码、字典更新及概率密度函数更新三个步骤。

3.根据权利要求2所述的一种基于时间移位和熵约束的地震信号编码方法，其特征在于，所述步骤S11具体算法如下：

输入参数：离线地震信号Y＝[y₁,...,y_K]，最大迭代次数I，初始字典D₀，最大时间偏移参数S，稀疏性L，每组里面段的数目N；

初始化：将输入信号Y划分为多个组，从而可以表示为Y＝[P₁,...,P_b,...,P_K/N]，其中P_b＝[y_(b-1)N+1,...,y_bN],b∈[1,K/N]，k和j为用于控制循环次数的变量，循环初始化k＝0，b＝0，j＝-S-1，字典初始化D⁽⁰⁾＝D₀；

主循环：令k＝k+1，如果k≤I，在第一循环内做以下运算：

令b＝b+1，如果在第二循环内做以下运算：

令j＝j+1，如果j≤S，在第三循环内做以下运算：

S111：对P_b沿着迹的方向进行整体平移，移动量为j，得到平移后的组P_b ^j；

S112：基于字典D^(k-1)，采用ORMP算法对P_b ^j进行稀疏表示，得到稀疏系数W_b ^j；

S113：计算残差ej，其中||·||_F代表Frobenius范数，即矩阵元素绝对值的平方和再开平方；

第三循环结束

S114：计算残差e^j(j∈[-S,S])的最小值，并将最小值对应的j赋值给变量T^(k) _b，T^(k) _b代表在循环第k次时计算出的时间偏移参数；

S115：将P_b整体移动T^(k) _b，并用移动后的结果更新P_b；

第二循环结束

S116：基于所有时间偏移后的地震信号数据，通过MOD方法训练并更新字典，生成字典D^(k)；

第一循环结束

输出：重新对齐的地震信号数据P_b(b∈[1,K/N])，时间偏移参数

4.根据权利要求3所述的一种基于时间移位和熵约束的地震信号编码方法，其特征在于，所述步骤S12具体流程如下：

输入参数：重新对齐的地震信号最大迭代次数I，初始字典D₀，初始概率密度函数χ₀，阈值τ，稀疏性L，量化步长Δ；

初始化：k，i，j和α为用于控制循环次数的变量，循环初始化参数k＝0，i＝0，j＝0，α＝-E-1，E为常数，字典初始化D⁽⁰⁾＝D₀，概率密度函数χ⁽⁰⁾＝χ₀；

主循环：令k＝k+1，如果k≤I，在第一循环内做以下运算：

令i＝i+1，如果i≤K，在第二循环内做以下运算：

S121:基于字典D^(k-1)，对重新对齐的地震信号采用ORMP方法进行稀疏表示，得到稀疏系数w_i；

S122:对稀疏系数w_i进行量化得到将中的非零值表示为{z₁,...,z_L}；

S123:初始化参数γ_j＝0,j∈[1,L]；

令j＝j+1，如果j≤L，在第三循环内做以下运算：

如果|z_j|≤τ(z_j为中的第j个非零值)，在第四循环内做以下运算：

S124:令α＝α+1，如果α≤E，计算

<mrow> <msub> <mi>TE</mi> <mi>&amp;alpha;</mi> </msub> <mo>=</mo> <mo>|</mo> <mo>|</mo> <msub> <mover> <mi>y</mi> <mo>^</mo> </mover> <mi>i</mi> </msub> <mo>-</mo> <munder> <mi>&amp;Sigma;</mi> <mrow> <mi>q</mi> <mo>&amp;NotEqual;</mo> <mi>j</mi> </mrow> </munder> <msub> <mi>d</mi> <mi>q</mi> </msub> <mrow> <mo>(</mo> <msub> <mi>z</mi> <mi>q</mi> </msub> <mo>+</mo> <msub> <mi>&amp;gamma;</mi> <mi>q</mi> </msub> <mo>)</mo> </mrow> <mi>&amp;Delta;</mi> <mo>-</mo> <msub> <mi>d</mi> <mi>j</mi> </msub> <mrow> <mo>(</mo> <msub> <mi>z</mi> <mi>j</mi> </msub> <mo>+</mo> <mi>&amp;alpha;</mi> <mo>)</mo> </mrow> <mi>&amp;Delta;</mi> <mo>|</mo> <msubsup> <mo>|</mo> <mn>2</mn> <mn>2</mn> </msubsup> <mo>+</mo> <mi>&amp;lambda;</mi> <mi>E</mi> <mi>n</mi> <mrow> <mo>(</mo> <msub> <mi>z</mi> <mi>j</mi> </msub> <mo>+</mo> <mi>&amp;alpha;</mi> <mo>)</mo> </mrow> </mrow>

其中，TE_α为所需要求解的代价函数，d_q代表字典D^(k-1)中的第q个字典单元，z_q为d_q对应的非零值，γ_q为d_q对应非零值的增量，d_j代表字典D^(k-1)中的第j个字典单元，λ是正则化参数，En(z_j+α)＝-plog₂p，p是z_j+α在χ^(k-1)中的概率；

第四循环结束

S125:找到TE_α的最小值对应的值更新z_j:＝z_j+α_j，α_j为代价函数TE_α最小时对应的α值；

第三循环结束

S126:对对应的非零值{z₁,...,z_L}进行更新；

第二循环结束

S127:采用MOD方法更新字典D^(k)；

S128:统计每个数值出现的概率，从而对概率密度函数χ^(k-1)进行更新；

第一循环结束

输出：概率密度函数χ＝χ^(I)和字典D＝D^(I)，χ^(I)和D^(I)分别表示经过I次循环后更新得到的概率密度函数和字典。

5.如权利要求4所述的一种基于时间移位和熵约束的地震信号编码方法，其特征在于，所述步骤S2具体包括以下步骤：

S21：对测试地震信号数据进行迹对齐，采用与S11相同的步骤，从而生成时间移位参数和重新对齐数据；

S22:基于步骤S12已经生成的概率密度函数χ和字典D，对重新对齐后的地震信号数据进行熵约束稀疏编码，具体步骤包含步骤S12中的步骤S121至S126，生成稀疏系数

S23:根据阈值τ，将稀疏系数中大于阈值的划分为低频系数，小于或者等于阈值且不为0的划分为高频系数；

S24:将稀疏系数中非零系数位置，低频系数位置进行传输，也可以对其进行算术编码后传输；

S25:将时间偏移参数进行传输；

S26:对低频系数进行差分编码，对高频系数基于已有的概率密度函数进行自适应算术编码。

6.权利要求5所述的一种基于时间移位和熵约束的地震信号编码方法，其特征在于，所述步骤S3具体包括以下步骤：

S31：对低频和高频编码码流分别进行差分解码和自适应算术解码，生成低频系数和高频系数；

S32：基于非零系数位置和低频系数位置，重建稀疏系数

S33：基于已有的字典D和重建稀疏系数生成重建信号

S34：根据所传输的时间偏移参数，对重建信号进行时间移位恢复，得到最终的重建地震信号数据。