基于块稀疏压缩感知的沙漠地震勘探随机噪声压制方法与流程

本发明属于一种沙漠地震勘探随机噪声压制方法，尤其是指基于块稀疏压缩感知的沙漠地震勘探随机噪声压制方法。

背景技术：

油气资源对于国民经济发展具有重要作用，地震勘探是探查油气资源的重要手段。随着国民经济发展的需要和科学技术水平的提高，对地震勘探技术的要求也由增强有效信号逐步发展到高信噪比、高分辨率、高保真度的“三高”水平。在“三高”要求中，高信噪比是基础。被沙漠覆盖的塔里木盆地作为我国最大的油气后备基地，是目前最为重要的地震勘探地区。沙漠地区采集到的地震资料伴随大量低频随机噪声，由于该噪声与地震子波序列性质相近且与有效信号频带大量重叠，导致在时频域压制低频噪声同时，基本没有可能提取到有效信号。面对复杂的随机噪声干扰，亟需建立沙漠地带低频随机噪声消减新对策，提高地震资料信噪比，满足油气地震勘探“三高”要求。

塔里木盆地沙漠地区地表起伏大，低降速带变化剧烈，受大沙丘与不均匀分布沙丘的影响，采集的地震数据中噪声干扰十分严重，导致地震记录信噪比较低，有效信号淹没在强能量噪声干扰、散射干扰之中。同时，沙漠地区干旱、松散的地表严重影响到激发能量的向下传播和地震检波器耦合，采集获得的有效信号能量很弱。地表起伏的沙丘横向速度变化很大，通过散射、绕射、交混回响等作用产生大量的随机噪声干扰，因此地震资料普遍存在反射波能量弱、信噪比低的问题。

在消减沙漠随机噪声的研究中，选用目前优秀的tfpf技术和小波变换技术，虽然一定程度上提高了信噪比，但效果不理想，存在的问题很多。时频峰值滤波(tfpf)方法是典型的时频分析方法，可以提高有效信号的信噪比和时频分辨率等。但是由于沙漠地震资料采集过程中受到散射、绕射、交混回响等作用产生大量复杂随机噪声干扰，且沙漠噪声幅值往往较大，tfpf在压制噪声的过程对有效信号衰减极为严重。小波变换是一种典型的多尺度分析方法，它有较好的时间频率局部化分析能力，目前小波域去噪技术在地震勘探领域中得到较好应用。但是，该方法在地震数据的消噪方面很大程度上依赖于阈值的选择，阈值选择的正确与否直接影响到噪声压制效果的好坏，并且该变换不具有方向性，在处理复杂的多维信号时具有一定局限性。为了有效提取低品质地震资料中的微弱地震信号，需要建立一套适应性更强的多约束综合信号提取方法。

技术实现要素：

本发明提供一种基于块稀疏压缩感知的沙漠地震勘探随机噪声压制方法，将同相轴定位与滤波相结合的信号提取技术作为沙漠地震信号处理的方法，利用非线性降维方法定位有效信号同相轴，并应用块稀疏压感知方法对定位后信号进行精确重构，提高信噪比。传统的cs理论忽略了源信号的结构特征，因此当处理某些具有特定结构的稀疏信号时重构效率较低，为了改善这一问题，提出了模型压缩感知方法，并实现了块稀疏信号的重构。块稀疏压缩充分利用块结构内部特点和块与块之间相关性，通过更少的观测值对信号进行解释分析，从而提高信号重构效率。充分利用地震信号的结构化块稀疏特性，分析有效信号块结构之间的相关性以及与噪声块结构之间的差异性，实现沙漠地区复杂随机噪声消减与地震有效信号提取。

本发明采取的技术方案是，包括下列步骤：

(1)沙漠地震记录获取：

结合沙漠地区野外实际的条件，进行地震激发，根据沙丘厚度因素调整测线部署及检波器布置，完成沙漠地震记录采集，得到沙漠地震信号记录x；

(2)沙漠地震记录定位处理：

(a)基于非线性降维模型的沙漠地震信号处理

将沙漠地震信号记录x表示为：

其中，x表示沙漠地震信号记录x的数据点，原始记录大小为m×n，表示为m个行向量，每个行向量由n个数据点描述，则这个数据就记为一个n维的记录；

计算x矩阵中每个数据点的欧式距离ai,j与局部切空间夹角bi,j，i,j为数据点的下角标，其范围为1≤i≤m，1≤j≤n，得到相似性矩阵w＝[wi,j]，其中，

矩阵l＝w-d，对l进行普分解，取最小2个特征值对应得到的特征向量，即得到沙漠地震信号降维后信号表示：

其中，u表示降维后信号的数据点，

(b)基于k-means聚类的沙漠地震信号定位

对于降维后的沙漠地震信号记录u，随机选择两个数据点作为质心μ1，μ2，计算每个数据点分别到两个质心的距离，并按照最近邻原则将其划分到距离最近的质心所属类别中，利用均值方法更新质心，不断迭代，直到质心不再发生变化，选出最佳质心，划分出两个质心所属的数据点，完成地震信号与噪声之间的分类，保留地震信号部分，重构得到定位后的沙漠地震信号x0，实现沙漠地震信号定位的目的；

(3)基于块稀疏压缩感知的沙漠随机噪声压制

(a)定位后沙漠信号块稀疏表示：

针对定位后的沙漠地震信号x0，选择过完备稀疏原子字典对其进行稀疏表示，稀疏原子字典表示为大小为m×m维的矩阵ψ＝[ψ1ψ2...ψm]∈r^m×m，其中m维向量ψi∈r^m(i＝1,2...m)为字典里的原子；

则信号x0可以在由ψ组成的稀疏域里唯一线性表示：

x0＝ψ·s，s为对应的稀疏系数表示，是信号在过完备稀疏原子字典上的接近最佳的稀疏表示；

块稀疏信号指的是在普通稀疏信号的基础上限定了分块稀疏的特性，即信号的零值和非零值的分布呈现簇类特性，非零值只在某些特定的分块位置上出现；

考虑地震信号的空间结构特点，充分利用地震信号的结构化块稀疏特性，对地震信号进行块稀疏化表示，得到块稀疏化表示为θ；

地震信号的块稀疏表示可以有效的降低地震信号稀疏度，使得地震信号的稀疏性增强，提高算法效率的同时，有效提高地震信号的重构精度；

(b)感知矩阵φ的设计

对于标准的稀疏信号，其感知矩阵需满足rip条件才能从测量向量中精确重构原始信号，针对块稀疏信号，要设定块rip条件；针对感知矩阵φ∈r^n×m,若对任意块稀疏信号x0∈r^m×n，存在常数δr满足条件：

其中，x0为信号x0中数据点，则称矩阵φ为满足块rip条件的感知矩阵；

(c)沙漠地震信号重构记录获取

针对得到的满足块rip条件的感知矩阵φ，根据压缩感知经典公式y＝φ·x0＝φ·ψ·s＝φ·ψ·θ,对信号进行测量，得到具有更少测量值的测量信号y；针对得到的测量信号y及感知矩阵φ，结合统计学中基于概率模型的块稀疏贝叶斯学习算法，定义s的先验分布为p(s)～n(0,∑0)，其中σ0＝diag{γ1b1,γ2b2...γq-1bq-1,γqbq}是一个主对角矩阵，q为下角标，其主对角线元素为每个块内元素的协方差。根据贝叶斯公式求解s的后验分布，构建如下对数似然函数：

其中，i为单位矩阵，在贝叶斯学习过程中利用估计出的参数λ，γi和bi，可获得上述参数的后验概率函数，从而准确地估计出s的值，得到根据估算出的值，得到噪声压制后的沙漠地震信号完成对地震信号的精确重构。

所述步骤(3)中：

其中θ表示块稀疏系数，表示第一列的第h1块，此块的长度为s为系数点，g1为第一列某一系数点的横坐标，n为系数点纵坐标，

本发明的优点在于：将基于块稀疏结构的压缩感知理论应用到地震信号处理领域，并将统计学方法融合到块稀疏重构算法中，在精确重构有效信号的同时有效压制了复杂强随机噪声。传统的压缩感知算法没有充分考虑信号的时空域结构特性，使重构的有效信号幅值与相位保持不理想，甚至产生波形畸变，导致同相轴恢复不连续。本专利采用基于块稀疏结构的压缩感知理论，在充分利用地震信号的块稀疏空间结构特性的基础上，根据有效信号块结构之间的相关性以及与噪声块结构之间的差异性，结合统计学的贝叶斯重构算法，实现复杂随机噪声的消减与有效信号重构相结合的信号提取技术。

附图说明

图1是截取的塔里木盆地沙漠地带随机噪声记录，其中采样频率为500hz，截取道数为24道，采样点数为800点；

图2是模拟合成的含有不同事件的纯净地震信号，其中采样频率为500hz，采样点数为800点，道数为24道；

图3是将实际的沙漠地带随机噪声记录与模拟的纯净地震信号合成的沙漠地带地震信号含噪记录；

图4是基于小波变换算法压制沙漠随机噪声的效果图；

图5是基于tfpf算法压制沙漠随机噪声的效果图；

图6(a)是基于非线性降维算法的沙漠地震有效信号定位图；

图6(b)是基于传统压缩感知算法压制沙漠随机噪声的效果图；

图6(c)是基于块稀疏压缩感知算法压制沙漠随机噪声的效果图。

具体实施方式

包括下列步骤：

(1)沙漠地震记录获取：

塔里木盆地沙漠地区，地表大多为沙漠覆盖，地表沙丘厚度大，起伏大，对地震波的吸收衰减非常严重，对塔里木盆地沙漠地区进行地震勘探，首先进行地表结构调查，通过潜水面测量准确获取潜水面的埋深情况，综合微测井、小折射成果，合理准确设计井深，并通过沙丘曲线调查建立好沙丘曲线数据库，解决大沙丘的深井激发设计问题，然后，结合野外实际的条件，进行地震激发，最后，根据沙丘厚度等因素调整测线部署及检波器布置，完成沙漠地震记录采集，得到沙漠地震信号记录x；

(2)沙漠地震记录定位处理：

(a)基于非线性降维模型的沙漠地震信号处理

将沙漠地震信号记录x表示为：

其中，x表示沙漠地震信号记录x的数据点，原始记录大小为m×n，表示为m个行向量，每个行向量由n个数据点描述，则这个数据就记为一个n维的记录；

计算x矩阵中每个数据点的欧式距离ai,j与局部切空间夹角bi,j，i,j为数据点的下角标，其范围为1≤i≤m，1≤j≤n，得到相似性矩阵w＝[wi,j]，其中，

矩阵l＝w-d，对l进行普分解，取最小2个特征值对应得到的特征向量，即得到沙漠地震信号降维后信号表示：

其中，u表示降维后信号的数据点，

(b)基于k-means聚类的沙漠地震信号定位

对于降维后的沙漠地震信号记录u，随机选择两个数据点作为质心μ1，μ2，计算每个数据点分别到两个质心的距离为：

其中,j为迭代次数，t为降维后信号u中的数据点下标标号，μi表示质心(i＝1,2)；

并按照最近邻原则将其划分到距离最近的质心所属类别中，利用均值方法更新质心，不断迭代，直到质心不再发生变化，选出最佳质心，划分出两个质心所属的数据点，完成地震信号与噪声之间的分类，保留地震信号部分，重构得到定位后的沙漠地震信号x0，实现沙漠地震信号定位的目的；

(3)基于块稀疏压缩感知的沙漠随机噪声压制

(a)定位后沙漠信号块稀疏表示：

针对定位后的沙漠地震信号x0，选择过完备稀疏原子字典对其进行稀疏表示，稀疏原子字典表示为大小为m×m维的矩阵ψ＝[ψ1ψ2...ψm]∈r^m×m，其中m维向量ψi∈r^m(i＝1,2...m)为字典里的原子；

则信号x0可以在由ψ组成的稀疏域里唯一线性表示：

x0＝ψ·s，s为对应的稀疏系数表示，是信号在过完备稀疏原子字典上的接近最佳的稀疏表示；

块稀疏信号指的是在普通稀疏信号的基础上限定了分块稀疏的特性，即信号的零值和非零值的分布呈现簇类特性，非零值只在某些特定的分块位置上出现；

考虑地震信号的空间结构特点，充分利用地震信号的结构化块稀疏特性，对地震信号进行块稀疏化表示：

其中θ表示块稀疏系数，表示第一列的第h1块，此块的长度为s为系数点，g1为第一列某一系数点的横坐标，n为系数点纵坐标，以此类推，得到块稀疏信号的块稀疏表示形式θ；

地震信号的块稀疏表示可以有效的降低地震信号稀疏度，使得地震信号的稀疏性增强，提高算法效率的同时，有效提高地震信号的重构精度；

(b)感知矩阵φ的设计

对于标准的稀疏信号，其感知矩阵需满足rip条件才能从测量向量中精确重构原始信号，针对块稀疏信号，要设定块rip条件；针对感知矩阵φ∈r^n×m,若对任意块稀疏信号x0∈r^m×n，存在常数δr满足条件：

其中，x0为信号x0中数据点，则称矩阵φ为满足块rip条件的感知矩阵；

(c)沙漠地震信号重构记录获取

针对得到的满足块rip条件的感知矩阵φ，根据压缩感知经典公式y＝φ·x0＝φ·ψ·s＝φ·ψ·θ,对信号进行测量，得到具有更少测量值的测量信号y；针对得到的测量信号y及感知矩阵φ，结合统计学中基于概率模型的块稀疏贝叶斯学习算法，定义s的先验分布为p(s)～n(0,∑0)，其中σ0＝diag{γ1b1,γ2b2...γq-1bq-1,γqbq}是一个主对角矩阵，q为下角标，其主对角线元素为每个块内元素的协方差。根据贝叶斯公式求解s的后验分布，构建如下对数似然函数：

其中，i为单位矩阵，在贝叶斯学习过程中利用估计出的参数λ，γi和bi，可获得上述参数的后验概率函数，从而准确地估计出s的值，得到根据估算出的值，得到噪声压制后的沙漠地震信号完成对地震信号的精确重构。

应用举例：

模拟合成纯净地震信号，其中采样频率为500hz，采样点数为800点，主频为35hz，道间距为25m，幅值为1米，道数为24道，如图2所示。将纯净地震信号与截取的塔里木盆地沙漠地区的实际随机噪声记录相结合，合成沙漠地带地震信号含噪记录，如图3所示。从图中可以看出沙漠噪声频率相较于有效信号频率低，幅值大，在噪声去除方面存在着一定的困难。

本发明提出的基于块稀疏压缩感知的沙漠地震勘探随机噪声压制方法，具体步骤如下：

(1)沙漠地震记录定位处理

将图3合成的沙漠地震信号含噪记录进行非线性降维处理，得到降维后数据。然后利用基于k-means的聚类算法将有效信号与噪声进行分类，保留有效信号部分，实现信号的无损定位，如6(a)所示；

(2)基于块稀疏压缩感知的沙漠随机噪声压制

由于定位后信号同相轴上依旧含有部分噪声，所以本发明应用基于块稀疏压缩感知的算法对定位后信号进行重构，压制同相轴上的沙漠随机噪声，得到去噪后的有效信号，如图6(c)所示。

将本发明中基于块稀疏压缩感知的沙漠随机噪声压制方法与基于传统压缩感知的沙漠随机噪声压制方法进行对比，从图6(c)与图6(b)中可以看出，本发明的有效信号恢复较完整，幅度保持更好，同相轴更加连续清晰。说明了块稀疏压缩感知算法相较于传统压缩感知算法的优越性。

利用小波变换算法，tfpf算法处理沙漠地震信号含噪记录，从图4小波变换算法处理的效果图中可以看出，噪声还是被较多的保留，并没有达到压制噪声的目的。从图5tfpf算法处理的效果图中看出，相较于小波去噪效果，tfpf算法噪声压制效果相对较好一些，但是信号保幅较差，有效信号能量损失严重，有畸变产生。

综上所述，本发明方法在压制沙漠地区随机噪声，提高地震资料的信噪比，分辩率方面上具有明显效果，更加有利于发现地震资料中断层的存在，为沙漠地区复杂环境下的石油勘探、天然气勘探提供压制噪声新对策。