基于稀疏表示的地震资料处理方法及系统与流程

本发明涉及地震信号处理领域，具体涉及一种基于稀疏表示的地震资料处理方法和系统。

背景技术：

随着油气勘探开发的逐步深入，勘探目标向着中深层复杂断块油藏、岩性隐蔽性油藏以及地层油藏等方向发展，相应的油气勘探开发难度越来越大。因此深化油气勘探开发对地震成像的精度自然提出了更高的需求。地震成像资料是深化油气勘探开发的重要资料，提高地震成像的精度以满足油气田勘探开发的需求至关重要。地震资料的分辨率问题是地震资料处理需要解决的核心问题，利用现有常规采集资料提高地震分辨率是目前地震资料处理的迫切需求。

现有提高地震资料分辨率的方法主要有反褶积处理方法和反q滤波处理方法。其中，反褶积处理方法存在人为假设，且当拓频要求高时，信号保真性差；反q滤波处理方法需要比较准确的等效q值场，只能处理纵向上的吸收衰减。

因此，期待一种效果更好的地震资料处理方法以获得高分辨率的地震资料图像。

技术实现要素：

本发明的目的是提供一种基于稀疏表示的地震资料处理方法和系统，以提高已知常规地震采集资料图像的分辨率。

本发明一方面提供一种基于稀疏表示的地震资料处理方法，包括：

步骤1：通过k-svd字典训练算法进行训练，建立字典d；

步骤2：从现有地震资料中选择多个亚样本，针对每个亚样本，通过stomp算法分解出解向量集合x，然后通过k-svd字典训练算法更新字典d中的每个原子以及解向量集合x中相应的行；

所述步骤2包括针对字典中的每一列，执行以下子步骤201-204：

子步骤201：根据约束条件确定索引集合εk；

子步骤202：计算误差矩阵ek；

子步骤203：基于误差矩阵ek和索引集合εk计算收缩矩阵针对收缩矩阵进行svd分解，根据分解结果对所述列进行更新；

子步骤204：重复执行子步骤201-203，根据更新的列计算误差矩阵，直到误差矩阵达到最小，将对应于最小误差矩阵的列作为最终列；

步骤3：基于所有最终列获得最终字典，根据最终字典和待处理地震资料得到最终的重建结果。

优选地，所述步骤1包括：

选择亚样本集合，并从s所述亚样本集合中随机抽取多个亚样本；

基于每个亚样本，通过k-svd字典训练算法进行训练，得到每个亚样本对应的子字典d(j),j＝1,……,p，p表示亚样本的数量；

对所有的子字典d(j)进行加权求和，得到字典d：

其中，wj表示第j个字典对应的加权系数，

优选地，所述约束条件为：

其中，n为从现有地震资料中选择的亚样本的数量，表示解向量集合x中与dk相应的第k行，dk表示字典d的第k列。

优选地，所述误差矩阵ek为：

其中，表示n个亚样本的集合。

优选地，所述收缩矩阵为：

其中，ωk＝n×|εk|。

优选地，所述针对收缩矩阵进行svd分解，根据分解结果对所述列进行更新包括：

对矩阵进行svd分解，则有其中矩阵u和v的列矢量均为正交基；

用矩阵u的第一列更新字典d的第k列dk，并用矩阵v的第一列和δ(1,1)的乘积更新其中，表示去掉零输入后的结果。

优选地，从现有地震资料中选择的亚样本数量为500-2000个。

本发明另一方面提供一种基于稀疏表示的地震资料处理系统，该系统包括：

存储器，存储有计算机可执行指令；

处理器，所述处理器运行所述存储器中的计算机可执行指令，执行以下步骤：

步骤1：通过k-svd字典训练算法进行训练，建立字典d；

步骤2：从现有地震资料中选择多个亚样本，针对每个亚样本，通过stomp算法分解出解向量集合x，然后通过k-svd字典训练算法更新字典d中的每个原子以及解向量集合x中相应的行；

所述步骤2包括针对字典中的每一列，执行以下子步骤201-204：

子步骤201：根据约束条件确定索引集合εk；

子步骤202：计算误差矩阵ek；

子步骤203：基于误差矩阵ek和索引集合εk计算收缩矩阵针对收缩矩阵进行svd分解，根据分解结果对所述列进行更新；

子步骤204：重复执行子步骤201-203，根据更新的列计算误差矩阵，直到误差矩阵达到最小，将对应于最小误差矩阵的列作为最终列；

步骤3：基于所有最终列获得最终字典，根据最终字典和待处理地震资料得到最终的重建结果。

优选地，所述步骤1包括：

选择亚样本集合，并从s所述亚样本集合中随机抽取多个亚样本；

基于每个亚样本，通过k-svd字典训练算法进行训练，得到每个亚样本对应的子字典d(j),j＝1,……,p，p表示亚样本的数量；

对所有的子字典d(j)进行加权求和，得到字典d：

其中，wj表示第j个字典对应的加权系数，

优选地，所述约束条件为：

其中，n为从现有地震资料中选择的亚样本的数量，表示解向量集合x中与dk相应的第k行，dk表示字典d的第k列。

与现有技术相比，本发明的有益效果在于：本发明的处理方法使得处理后的地震资料图像具有良好的清晰度；通过随机抽样可以降低单个k-svd训练算法的误差，增强了字典训练的准确性和稳定性；通过设定稀疏字典的大小，在保证计算效率的同时，达到良好的计算效果。

本发明的方法具有其它的特性和优点，这些特性和优点从并入本文中的附图和随后的具体实施例中将是显而易见的，或者将在并入本文中的附图和随后的具体实施例中进行详细陈述，这些附图和具体实施例共同用于解释本发明的特定原理。

附图说明

通过结合附图对本发明示例性实施例进行更详细的描述，本发明的上述以及其它目的、特征和优势将变得更加明显。

图1显示根据本发明示例性实施例的基于稀疏表示的地震资料处理方法的流程图。

具体实施方式

下面将参照附图更详细地描述本发明。虽然附图中显示了本发明的优选实施例，然而应该理解，可以以各种形式实现本发明而不应被这里阐述的实施例所限制。相反，提供这些实施例是为了使本发明更加透彻和完整，并且能够将本发明的范围完整地传达给本领域的技术人员。

在地震图像采集过程中，由于仪器性能不完善，以及受到图像成像、传输、存储等处理过程中不同类型噪声的干扰，从而导致图像分辨率的下降，影响了地震图像质量。地震图像复原问题可以表示为如下数学模型：

y＝φhx+σ

其中，x表示原始的高分辨率图像；φ表示图像下采样矩阵；h表示退化函数；σ表示噪声；y表示低分辨率图像。

在忽略噪声影响的情况下地震图像复原问题可以表示为：

y＝φx+σ

由于自然界中的图像包含很多重复的部分，即同一个图像块在同一尺度和不同尺度都可以找到相同或相似的图像块，因此超分辨率复原大多是基于图像块进行的。以这种思想为基础，本发明提出了一种基于图像自相似性与稀疏表示模型的超分辨率地震图像处理方法，在本发明的方法中，利用输入的低分辨率图像y的s多级下采样图像y-1,y-2,...y-s作为训练集，训练出复原字典d，再结合稀疏表示理论，将超分辨率复原模型表示为：

x≈dα

其中，λ为正则化参数；x为待估计的超分辨率图像。

考虑到图像的结构自相似性，利用非局部均值滤波的思想，认为图像x中的所有图像块都可以由邻域内非局部自相似块的加权平均近似表示，即：

其中，s为图像的非局部自相似矩阵，表示邻域内非局部自相似块的加权平均近似。

由上式得到基于图像自相似性与稀疏表示的超分辨率复原模型：

x≈dα

y＝φx

为进一步提高图像复原质量，对x和α进行迭代求解可以得到最终的复原图像。

图1显示根据本发明示例性实施例的基于稀疏表示的地震资料处理方法的流程图，其包括以下步骤：

步骤1：通过k-svd字典训练算法进行训练，建立字典d。

具体地，步骤1包括以下子步骤：

子步骤101：选择亚样本集合，并从亚样本集合中随机抽取多个亚样本；

子步骤102：基于每个亚样本，通过k-svd字典训练算法进行训练，得到每个亚样本对应的子字典d(j),j＝1,……,p，p表示亚样本的数量；

本发明实施例的字典训练方法基于aharon和elad提出的k-svd字典训练算法。该算法的目的是找到一组数据能量分布的共同方向，作为字典基底，得到信号的稀疏表示。该算法首先要求初始化字典，然后用该字典表示测试数据，得到满足稀疏度但存在一定误差的稀疏系数矩阵。然后，逐行剥离字典中基底对数据表示的贡献，计算当前误差矩阵，用方法分解收缩后的误差矩阵，取能量最大的行列更新字典与系数。对未达到要求的字典，在开始更新前应重新生成系数矩阵，对于字典中贡献小、有重复部分的列，用数据中被表示的最少的列来填充。

本发明实施例在字典训练过程中使用装袋算法思想，通过引入随机抽样的多个字典进行训练。根据均匀分布，从数据训练样本中抽取n个样本得到n个亚样本，每个亚样本集合和原始数据的分布一样，抽样过程使用一种有放回的策略。

子步骤103：对所有的子字典d(j)进行加权求和，得到字典d：

其中wj表示第j个字典对应的加权系数，

基于每个亚样本进行训练分别得到一个字典，各字典之间相互独立，互相不干扰。为了加快运行速度，可以考虑使用并行计算的方式，多个字典的训练并行执行，各个字典之间的参数互相不干扰，只是在最终的结果输出阶段，把所有字典进行加权，得到最终的加权字典。

步骤2：从现有地震资料中选择多个亚样本，针对每个亚样本，通过stomp(分段正交匹配追踪)算法分解出解向量集合x，然后通过k-svd字典训练算法更新字典d中的每个原子以及解向量集合x中相应的行；

在实施例中，以表示n个亚样本的集合，表示y的解向量集合，其是稀疏矩阵，n表示亚样本的数量。

其中，di表示字典d的第i列，xi表示x的第i行。

所述步骤2包括针对字典d中的每一列，执行以下子步骤201-204：

子步骤201：根据约束条件确定索引集合εk；

根据约束条件确定索引集合εk，其是用到dk图像特征集合{yi}的索引集合；

子步骤202：计算误差矩阵ek；

通过下式计算误差矩阵ek：

其中，误差矩阵ek表示的是剥离字典d中原子成分dk的成本在所有n个样本中造成的误差矩阵，表示x中与dk相应的第k行；

子步骤203：基于误差矩阵ek和索引集合εk计算收缩矩阵

通过获得矩阵ek去掉零输入后的收缩结果

通过获得矩阵ek去掉零输入后的收缩结果

其中，ωk＝n×|εk|；

然后，对收缩矩阵进行svd分解，根据分解结果对字典d的第k列dk进行更新；

对矩阵进行svd分解，则有其中矩阵u和v的列矢量均为正交基；

用矩阵u的第一列更新字典d的第k列，并用矩阵v的第一列和δ(1,1)的乘积更新其中，表示去掉零输入后的结果，△代表奇异值；

子步骤204：重复执行子步骤201-203，根据更新的列计算误差矩阵，直到误差矩阵达到最小，将对应于最小误差矩阵的列作为最终列；

步骤3：基于所有最终列获得最终字典，根据最终字典和待处理地震资料得到最终的重建结果。

本发明实施例还提供一种基于稀疏表示的地震资料处理系统，该系统包括：

存储器，存储有计算机可执行指令；

处理器，所述处理器运行所述存储器中的计算机可执行指令，执行以下步骤：

步骤1：通过k-svd字典训练算法进行训练，建立字典d；

步骤2：从现有地震资料中选择多个亚样本，针对每个亚样本，通过stomp算法分解出解向量集合x，然后通过k-svd字典训练算法更新字典d中的每个原子以及解向量集合x中相应的行；

所述步骤2包括针对字典中的每一列，执行以下子步骤201-204：

子步骤201：根据约束条件确定索引集合εk；

子步骤202：计算误差矩阵ek；

子步骤203：基于误差矩阵ek和索引集合εk计算收缩矩阵针对收缩矩阵进行svd分解，根据分解结果对所述列进行更新；

子步骤204：重复执行子步骤201-203，根据更新的列计算误差矩阵，直到误差矩阵达到最小，将对应于最小误差矩阵的列作为最终列；

步骤3：基于所有最终列获得最终字典，根据最终字典和待处理地震资料得到最终的重建结果。

作为优选方案，步骤1包括：

选择亚样本集合，并从s所述亚样本集合中随机抽取多个亚样本；

基于每个亚样本，通过k-svd字典训练算法进行训练，得到每个亚样本对应的子字典d(j),j＝1,……,p，p表示亚样本的数量；

对所有的子字典d(j)进行加权求和，得到字典d：

其中，wj表示第j个字典对应的加权系数，

作为优选方案，约束条件为：

其中，n为从现有地震资料中选择的亚样本的数量，表示解向量集合x中与dk相应的第k行，dk表示字典d的第k列。

作为优选方案，所述误差矩阵ek为：

其中，表示n个亚样本的集合。

作为优选方案，所述收缩矩阵为：

其中，ωk＝n×|εk|。

作为优选方案，所述针对收缩矩阵进行svd分解，根据分解结果对所述列进行更新包括：

对矩阵进行svd分解，则有其中矩阵u和v的列矢量均为正交基；

用矩阵u的第一列更新字典d的第k列dk，并用矩阵v的第一列和δ(1,1)的乘积更新其中，表示去掉零输入后的结果。

作为优选方案，从现有地震资料中选择的亚样本数量为500-2000个。

实施例1：样本数量和字典大小的试验：

采用地震图像进行试验，随机选择500、1000、2000、3000、4000、5000个样本进行训练学习，结果如表1所示。

表1不同训练样本数峰值信噪比psnr结果对比

当训练样本个数大于2000的时候，结果变化不大。稀疏字典增大时计算量也相应增加。为了保证计算效果及计算效率，稀疏字典的大小选为500至2000，本发明中稀疏字典的大小选为1500。

实施例2：地震图像处理质量试验：

从训练图像集中随机选取50000对不同的高、低分辨率图像对，并且求取特征向量进行训练。为了有效对比处理后图像和输入图像，并平衡计算效率和结果质量，将稀疏字典大小固定为1500，采样率取0.5。由表2可以看出，本发明方法处理后的地震图像峰值信噪比明显提高。

表2本发明与常规方法的结果对比

以上已经描述了本发明的实施例，上述说明是示例性的，并非穷尽性的，并且也不限于所披露的各实施例。在不偏离所说明的各实施例的范围和精神的情况下，对于本技术领域的普通技术人员来说许多修改和变更都是显而易见的。本文中所用术语的选择，旨在最好地解释各实施例的原理、实际应用或对市场中的技术的改进，或者使本技术领域的其它普通技术人员能理解本文披露的各实施例。