一种基于同相轴形态差异的面波提纯方法及装置与流程

本发明涉及地震数据处理技术领域，具体而言，涉及一种基于同相轴形态差异的面波提纯方法及装置。

背景技术：

面波探测凭借衰减小、信噪比高、无损探测的特点，其频散特性反映了浅层横波速度结构的变化，在探测地表地质结构、工程探测等方面得到广泛应用。为精确提取面波频散曲线，要根据面波与干扰波的特征差异、借助信号处理方法分离波场，与反射波勘探中压制面波的方法不同，面波勘探中要求处理前后不能影响面波频散信息。现有的面波提纯方法只是单纯地基于某些数学变换在变换域进行波场分离，但由于面波和干扰波在变换域的表示系数不够稀疏，存在波场分离效果不佳等问题。

技术实现要素：

有鉴于此，本发明实施例的目的在于提供一种基于同相轴形态差异的面波提纯方法及装置，以改善上述的问题。

本发明采用的技术方案如下：

本发明实施例提供了一种基于同相轴形态差异的面波提纯方法，所述方法包括，

s10：对地震反射波数据yr进行初始化；

s20：根据原始地震数据y及所述反射波数据yr，求取面波在频域高分辨率lrt变换域的稀疏表示系数，并重构面波数据yg；

s30：根据所述原始地震数据y及重构后的所述面波数据yg，求取反射波在时域高分辨率hrt变换域的稀疏表示系数，并重构反射波数据yr；

迭代进行步骤s20、步骤s30以获得面波提纯结果。

进一步地，所述根据原始地震数据y及所述反射波数据yr，求取面波在频域高分辨率lrt变换域的稀疏表示系数，并重构面波数据yg的步骤包括：

根据原始地震数据y以及反射波数据yr计算获得第一残差数据δy1，其中，δy1＝y-yr；

对第一残差数据δy1进行频域高分辨率lrt变换获得第一变换系数zg；

对获得的第一变换系数zg利用第一阈值σg进行硬阈值处理得到面波在频域高分辨率lrt变换域的稀疏表示系数zg'；

对面波稀疏表示系数zg'进行频域反lrt变换，重构面波数据yg。

进一步地，所述根据所述原始地震数据y及重构后的面波数据yg，求取反射波在时域高分辨率hrt变换域的稀疏表示系数，并重构反射波数据yr的步骤包括：

根据原始地震数据y以及重构后的面波数据yg计算获得第二残差数据δy2，其中δy2＝y-yg；

对第二残差数据δy2进行时域高分辨率hrt变换得到第二变换系数zr；

对获得的第二变换系数zr利用第二阈值σr进行硬阈值处理得到反射波在时域高分辨率hrt变换域的稀疏表示系数zr'；

对面波稀疏表示系数zr'进行时域反hrt变换，重构反射波数据yr。

进一步地，设定一循环变量k、给定迭代次数k，循环变量k的初始值k＝0，在所述步骤s20之前，所述方法还包括：

令循环变量k＝k+1。

进一步地，在每次迭代过程中，更新第一阈值σg、第二阈值σr的值，具体包括：

其中，k为循环变量，k为迭代次数。

进一步地，所述对地震反射波数据yr进行初始化包括：

令所述反射波数据初始值为yr＝0；

进一步地，所述迭代进行步骤s20、步骤s30以获得面波提纯结果的步骤包括：

当循环变量k小于或等于迭代次数k时，重复执行步骤s20～s30以获得面波提纯结果。

进一步地，所述方法还包括：

当循环变量k大于迭代次数k时，停止迭代，输出面波提纯数据。

进一步地，在对面波数据yg、反射波数据yr进行初始化之前，所述方法还包括：

获取原始地震数据y。

本发明还提供了一种基于同相轴形态差异的面波提纯装置，其特征在于，所述装置用于执行所述的基于同相轴形态差异的面波提纯方法，所述装置包括：

初始化单元，用于对反射波数据yr进行初始化；

第一重构单元，用于根据原始地震数据y及所述反射波数据yr，求取面波在频域高分辨率lrt变换域的稀疏表示系数，并重构面波数据yg；

第二重构单元，用于根据所述原始地震数据y及重构后的所述面波数据yg，求取反射波在时域高分辨率hrt变换域的稀疏表示系数，并重构反射波数据yr；

所述第一重构单元及所述第二重构单元还用于进行多次迭代重构以获得面波提纯结果。

相对现有技术，本发明具有以下有益效果：

本发明提供的一种基于同相轴形态差异的面波提纯方法及装置，通过对反射波数据yr进行初始化；根据原始地震数据y及所述反射波数据yr，求取面波在频域高分辨率lrt变换域的稀疏表示系数，并重构面波数据yg；根据所述原始地震数据y及重构后的所述面波数据yg，求取反射波在时域高分辨率hrt变换域的稀疏表示系数，并重构反射波数据yr；迭代进行重构计算步骤以求解mca框架下的稀疏优化问题并获得面波提纯结果。通过对面波和反射波进行多次变换、重构，当面波和反射波的频率、视速度接近时，仍能有效提取面波，降低反射波对面波频散能谱的干扰，为频散分析提供了高精度的面波数据。

为使本发明的上述目的、特征和优点能更明显易懂，下文特举较佳实施例，并配合所附附图，作详细说明如下。

附图说明

为了更清楚地说明本发明实施方式的技术方案，下面将对实施方式中所需要使用的附图作简单地介绍，应当理解，以下附图仅示出了本发明的某些实施例，因此不应被看作是对范围的限定，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他相关的附图。

图1示出了本发明所提供的一种地震数据计算设备示意图。

图2示出了基于同相轴形态差异的面波提纯方法的流程图。

图3示出了图2中步骤s20的子步骤流程图。

图4示出了图2中步骤s30的子步骤流程图。

图5示出了基于同相轴形态差异的面波提纯装置的功能模块示意图。

图标：100-地震数据计算设备；101-存储器；102-存储控制器；103-处理器；104-外设接口；105-显示单元；106-输入输出单元；200-基于同相轴形态差异的面波提纯装置；210-获取单元；220-初始化单元；230-第一重构单元；240-第二重构单元；250-判断单元；260-输出单元。

具体实施方式

下面将结合本发明实施例中附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。通常在此处附图中描述和示出的本发明实施例的组件可以以各种不同的配置来布置和设计。因此，以下对在附图中提供的本发明的实施例的详细描述并非旨在限制要求保护的本发明的范围，而是仅仅表示本发明的选定实施例。基于本发明的实施例，本领域技术人员在没有做出创造性劳动的前提下所获得的所有其他实施例，都属于本发明保护的范围。

应注意到：相似的标号和字母在下面的附图中表示类似项，因此，一旦某一项在一个附图中被定义，则在随后的附图中不需要对其进行进一步定义和解释。

在本发明的描述中，还需要说明的是，在本文中，诸如第一和第二等之类的关系术语仅仅用来将一个实体或者操作与另一个实体或操作区分开来，而不一定要求或者暗示这些实体或操作之间存在任何这种实际的关系或者顺序。术语“包括”、“包含”或者其任何其他变体意在涵盖非排他性的包含，从而使得包括一系列要素的过程、方法、物品或者设备不仅包括那些要素，而且还包括没有明确列出的其他要素，或者是还包括为这种过程、方法、物品或者设备所固有的要素。在没有更多限制的情况下，由语句“包括一个……”限定的要素，并不排除在包括所述要素的过程、方法、物品或者设备中还存在另外的相同要素。对于本领域的普通技术人员而言，可以具体情况理解上述术语在本发明中的具体含义。

下面结合附图，对本发明的一些实施方式作详细说明。在不冲突的情况下，下述的实施例及实施例中的特征可以相互组合。

图1示出本发明较佳实施例提供的地震数据计算设备100的方框示意图。所述地震数据计算设备100可以是台式计算机、笔记本电脑、平板电脑、智能手机、个人数字助理(personaldigitalassistant，pda)等。所述地震数据计算设备100包括基于同相轴形态差异的面波提纯装置200、存储器101、存储控制器102、处理器103、外设接口104、显示单元105、输入输出单元106。

所述存储器101、存储控制器102、处理器103、外设接口104、显示单元105、输入输出单元106各元件相互之间直接或间接地电性连接，以实现数据的传输或交互。例如，这些元件相互之间可通过一条或多条通讯总线或信号线实现电性连接。所述基于同相轴形态差异的面波提纯装置200包括至少一个可以软件或固件(firmware)的形式存储于所述存储器101中或固化在所述地震数据计算设备100的操作系统(operatingsystem，os)中的软件功能模块。所述处理器103用于执行存储器101中存储的可执行模块，例如所述基于同相轴形态差异的面波提纯装置200包括的软件功能模块或计算机程序。

其中，存储器101可以是，但不限于，随机存取存储器(randomaccessmemory，ram)，只读存储器(readonlymemory，rom)，可编程只读存储器(programmableread-onlymemory，prom)，可擦除只读存储器(erasableprogrammableread-onlymemory，eprom)，电可擦除只读存储器(electricerasableprogrammableread-onlymemory，eeprom)等。其中，存储器101用于存储程序，所述处理器103在接收到执行指令后，执行所述程序，本发明任一实施例揭示的由过程定义的服务器所执行的方法可以应用于处理器103中，或者由处理器103实现。

处理器103可能是一种集成电路芯片，具有信号的处理能力。上述的处理器103可以是通用处理器，包括中央处理器(centralprocessingunit，简称cpu)、网络处理器(networkprocessor，简称np)等；还可以是数字信号处理器(dsp)、专用集成电路(asic)、现场可编程门阵列(fpga)或者其他可编程逻辑器件、分立门或者晶体管逻辑器件、分立硬件组件。可以实现或者执行本发明实施例中的公开的各方法、步骤及逻辑框图。通用处理器可以是微处理器或者该处理器103也可以是任何常规的处理器103等。

所述外设接口104将各种输入/输出装置耦合至处理器103以及存储器101。在一些实施例中，外设接口104，处理器103以及存储控制器102可以在单个芯片中实现。在其他一些实例中，他们可以分别由独立的芯片实现。

显示单元105在所述地震数据计算设备100与用户之间提供一个交互界面(例如用户操作界面)或用于显示图像数据给用户参考。在本实施例中，所述显示单元105可以是液晶显示器或触控显示器。若为触控显示器，其可为支持单点和多点触控操作的电容式触控屏或电阻式触控屏等。支持单点和多点触控操作是指触控显示器能感应到来自该触控显示器上一个或多个位置处同时产生的触控操作，并将该感应到的触控操作交由处理器103进行计算和处理。

输入输出单元106用于提供给用户输入数据实现用户与所述地震数据计算设备100的交互。例如，可以用于设定迭代次数，设定阈值数值等等。

面波具有横向传播衰减小、信噪比高等特点，且频散特性反映了浅层横波速度结构的变化，在壳幔速度结构探测、工程探测等方面得到广泛应用。石油地震中，也有研究利用z分量上的基阶面波频散曲线反演浅层横波速度结构，用于计算ps波静校正量。而近些年基于高、基阶面波对地层参数的不同响应特征，利用基阶和高阶面波频散曲线进行联合反演，增加反演的稳定性，提高反演精度成为领域的研究热点之一。

多分量地震采集中，基阶面波在z分量上具有较强的能量，而高阶面波常常在x分量上发育；由于x分量上面波一般受反射波干扰，频散能谱连续性差，很难提取到准确的频散曲线，导致利用高阶面波反演精细的横波速度结构存在不确定性。为精确提取面波频散曲线，需要根据面波与干扰波的特征差异，提取较纯净的面波。由于ps反射波的频率、视速度与高阶面波接近，仅依靠频率、视速度差异难以有效提取面波。

因此，为利用x分量的面波，高精度的面波提取方法是关键。

第一实施例

请参阅图2，本实施例提供了一种基于同相轴形态差异的面波提纯方法，所述方法包括步骤s01～步骤s50。

步骤s01：获取原始地震数据。

获取原始地震数据y，于本实施例中，原始地震数据y指记录的单炮地震数据y。一般地，单炮地震数据y中包括反射波数据及面波数据。

步骤s10：对地震反射波数据yr进行初始化。

所述对地震反射波数据yr进行初始化，即是对反射波数据yr赋予初始值，于本实施例中，令反射波数据yr＝0。

在所述步骤s20之前，所述方法还包括，步骤s101。

步骤s101：令循环变量k＝k+1。

由于需要进行多次迭代循环计算，于本实施例中，设定一循环变量k、给定迭代次数k，循环变量k的初始值k＝0，令循环变量k＝k+1。即是说，在第一次迭代过程中，k＝k+1＝0+1＝1，所述迭代次数可以任意设置，于本实施例中，迭代次数k＝3。

分别采用频域高分辨率线性radon变换(linearradontransform，lrt)和时域高分辨率双曲radon变换(hyperbolicradontransform，hrt)稀疏表示面波和反射波，在形态成分分析(morphologicalcomponentanalysis，mca)框架下构建稀疏优化问题；采用stark等改进的分块坐标松弛算法求解。

步骤s20：根据原始地震数据y及所述反射波数据yr，求取面波在频域高分辨率lrt变换域的稀疏表示系数，并重构面波数据yg。

请参阅图3，于本实施例中，步骤s20包括子步骤s201～步骤s204。

步骤s201：根据原始地震数据y以及反射波数据yr计算获得第一残差数据δy1，其中，δy1＝y-yr。

原始地震数据y包括反射波数据yr及面波数据yg，利用y-yr，去除原始地震数据中的反射波数据的影响，余下的即是面波数据。

步骤s202：对第一残差数据δy1进行频域高分辨率lrt变换获得第一变换系数zg。

对第一残差数据δy1进行频域lrt，得到第一变换系数zg，即求取频散能谱。

步骤s203：对获得的第一变换系数zg利用第一阈值σg进行硬阈值处理得到面波在频域高分辨率lrt变换域的稀疏表示系数zg'。

对获得的第一变换系数zg利用第一阈值σg进行硬阈值处理得到面波在频域高分辨率lrt变换域的稀疏表示系数zg'。一般地，第一阈值σg的初值设置为0.4～0.5，在每一次的迭代过程中，其中，k为循环变量，k为迭代次数。对获得的第一变换系数zg利用第一阈值σg进行硬阈值处理，将第一变换系数zg中小于第一阈值的数值舍去，保留大于第一阈值σg的值，得到面波在频域高分辨率lrt变换域的稀疏表示系数zg'。

步骤s204：对面波稀疏表示系数zg'进行频域反lrt变换，重构面波数据yg。

对阈值处理后的第一变换系数zg'进行频域反lrt变换，对面波数据yg进行重构，频域反lrt公式为

d(f)＝l(f)m(f)；

式中：d(f)为指定频率f的时空域地震道的傅立叶系数向量，长度为nx×1；m(f)为指定频率f的radon域地震道的傅立叶系数向量，长度为np×1；l(f)为nx×np阶复矩阵

其中vi(i＝1,2,...,np)为视速度，xj(j＝1,2,...,nx)为偏移距，nx、np分别是地震数据的道数、速度参数的个数。频域高分辨率lrt满足方程：

其中，对角阵wd为数据加权矩阵diag(wd)i＝|(d-lm)i|^-1/2，反映数据的标准差；对角阵wm为radon系数加权矩阵diag(wm)i＝|mi|^-1/2，决定变换系数的稀疏性；i为单位矩阵；参数λ平衡时空域数据误差和radon变换系数稀疏性的比重；l^h是l的共轭转置矩阵。此方程采用迭代重加权最小二乘算法求解。

步骤s30：根据所述原始地震数据y及重构后的面波数据yg，求取反射波在时域高分辨率hrt变换域的稀疏表示系数，并重构反射波数据yr。

请参阅图4，于本实施例中，步骤s30包括子步骤s301～步骤s304。重构反射波数据yr的步骤与重构面波数据的步骤相似，不再进行详细说明。

步骤s301：根据原始地震数据y以及重构后的面波数据yg计算获得第二残差数据δy2，其中δy2＝y-yg。

步骤s302：对第二残差数据δy2进行时域高分辨率hrt变换得到第二变换系数zr。

步骤s303：对获得的第二变换系数zr利用第二阈值σr进行硬阈值处理得到反射波在时域高分辨率hrt变换域的稀疏表示系数zr'。

对获得的第二变换系数zr利用第二阈值σr进行硬阈值处理得到反射波在时域高分辨率hrt变换域的稀疏表示系数zr'。一般地，第二阈值σr的初值设置为0.1～0.2，在每一次的迭代过程中，其中，k为循环变量，k为迭代次数。对获得的第二变换系数zr利用第二阈值σr进行硬阈值处理，将第二变换系数zr中小于第二阈值σr的数值舍去，保留大于第二阈值σr的值，得到反射波在时域高分辨率hrt变换域的稀疏表示系数zr'。

步骤s304：对反射波稀疏表示系数zr'进行时域反hrt变换，重构反射波数据yr。

时域反hrt公式为

其中d(t,x)为时空域地震数据，x是偏移距，t是地震波双程旅行时，m(τ,v)为radon变换系数，v是均方根速度，τ为时间截距，写成矩阵-向量形式

d＝lm，

d是地震数据按道重排成的向量，长度(nx×nt)×1；m是双曲radon系数按速度重排成的向量，长度(nv×nτ)×1，nx、nt、nv、nτ分别是数据的道数、每道的采样点数、速度参数的个数、截距参数的个数；算子l不是按矩阵形式存取，仅代表d＝lm式所示的radon系数域中振幅叠加算法。在sabbione和sacchi提出的受限空间中求解时域hrt满足的方程：

式中地震数据d整体归一化，受限空间定义为

其中t是阈值，取值范围0<t<1。此方程采用左预条件正则方程共轭梯度算法求解。

步骤s40：确认循环变量k是否小于或等于迭代次数k。

确认迭代次数是否达到了设定的迭代次数k。由于在每一次迭代计算过程中，循环变量的值可以表征已经迭代的次数，因此将循环变量的值与迭代次数k进行比较，若当循环变量k小于或等于迭代次数k时，即还未达到预设的迭代次数，则执行步骤s10。若已经达到预设的迭代次数，即当循环变量k大于迭代次数k时，执行步骤s50。

步骤s50：停止迭代运算，输出面波提纯数据。

第二实施例

本实施例提供了一种基于同相轴形态差异的面波提纯装置200，所述基于同相轴形态差异的面波提纯装置200用于执行第一实施例提供的基于同相轴形态差异的面波提纯方法。

需要说明的是，本实施例提供的基于同相轴形态差异的面波提纯装置200，其基本原理与技术效果与第一实施例提供的基于同相轴形态差异的面波提纯方法基本相同，为简要描述，本实施例不再进行详细说明，本实施例未介绍详尽之处，请参阅第一实施例中的相关内容。

请参阅图5，图5示出了本实施例提供的基于同相轴形态差异的面波提纯装置200的功能单元示意图。

所述基于同相轴形态差异的面波提纯装置200包括获取单元210、初始化单元220、第一重构单元230、第二重构单元240、判断单元250，以及输出单元260。

其中，获取单元210，用于获取原始地震数据。

可以理解地，在一种优选实施例中，获取单元210可以用于执行步骤s01。

初始化单元220，用于对反射波数据yr进行初始化。

可以理解地，在一种优选实施例中，初始化单元220可以用于执行步骤s10。

第一重构单元230，用于根据原始地震数据y及所述反射波数据yr，求取面波在频域高分辨率lrt变换域的稀疏表示系数，并重构面波数据yg。

可以理解地，在一种优选实施例中，第一重构单元230可以用于执行步骤s20。

第二重构单元240，用于根据所述原始地震数据y及重构后的所述面波数据yg，求取反射波在时域高分辨率hrt变换域的稀疏表示系数，并重构反射波数据yr。

可以理解地，在一种优选实施例中，第二重构单元240可以用于执行步骤s30。

判断单元250，用于判断确认循环变量是否小于或等于迭代次数。

可以理解地，在一种优选实施例中，判断单元250可以用于执行步骤s40。

所述第一重构单元230及所述第二重构单元240还用于当循环变量小于或等于迭代次数时进行多次迭代重构以获得面波提纯结果。

输出单元260，用于当循环变量大于迭代次数时，输出面波提纯结果。

可以理解地，在一种优选实施例中，输出单元260可以用于执行步骤s50。

综上所述，本发明提供了一种基于同相轴形态差异的面波提纯方法及装置。通过对反射波数据yr进行初始化；根据原始地震数据y及所述反射波数据yr，求取面波在频域高分辨率lrt变换域的稀疏表示系数，并重构面波数据yg；根据所述原始地震数据y及重构后的所述面波数据yg，求取反射波在时域高分辨率hrt变换域的稀疏表示系数，并重构反射波数据yr；迭代进行重构计算步骤以求解mca框架下的稀疏优化问题并获得面波提纯结果。通过对面波和反射波进行多次变换、重构，当面波和反射波的频率、视速度接近时，仍能有效提取面波，降低反射波对面波频散能谱的干扰，为频散分析提供了高精度的面波数据。

在本申请所提供的几个实施例中，应该理解到，所揭露的装置和方法，也可以通过其它的方式实现。以上所描述的装置实施例仅仅是示意性的，例如，附图中的流程图和框图显示了根据本发明的多个实施例的装置、方法和计算机程序产品的可能实现的体系架构、功能和操作。在这点上，流程图或框图中的每个方框可以代表一个模块、程序段或代码的一部分，所述模块、程序段或代码的一部分包含一个或多个用于实现规定的逻辑功能的可执行指令。也应当注意，在有些作为替换的实现方式中，方框中所标注的功能也可以以不同于附图中所标注的顺序发生。例如，两个连续的方框实际上可以基本并行地执行，它们有时也可以按相反的顺序执行，这依所涉及的功能而定。也要注意的是，框图和/或流程图中的每个方框、以及框图和/或流程图中的方框的组合，可以用执行规定的功能或动作的专用的基于硬件的系统来实现，或者可以用专用硬件与计算机指令的组合来实现。

另外，在本发明各个实施例中的各功能模块可以集成在一起形成一个独立的部分，也可以是各个模块单独存在，也可以两个或两个以上模块集成形成一个独立的部分。

所述功能如果以软件功能模块的形式实现并作为独立的产品销售或使用时，可以存储在一个计算机可读取存储介质中。基于这样的理解，本发明的技术方案本质上或者说对现有技术做出贡献的部分或者该技术方案的部分可以以软件产品的形式体现出来，该计算机软件产品存储在一个存储介质中，包括若干指令用以使得一台计算机设备(可以是个人计算机，服务器，或者网络设备等)执行本发明各个实施例所述方法的全部或部分步骤。而前述的存储介质包括：u盘、移动硬盘、只读存储器(rom，read-onlymemory)、随机存取存储器(ram，randomaccessmemory)、磁碟或者光盘等各种可以存储程序代码的介质。

需要说明的是，在本文中，诸如第一和第二等之类的关系术语仅仅用来将一个实体或者操作与另一个实体或操作区分开来，而不一定要求或者暗示这些实体或操作之间存在任何这种实际的关系或者顺序。而且，术语“包括”、“包含”或者其任何其他变体意在涵盖非排他性的包含，从而使得包括一系列要素的过程、方法、物品或者设备不仅包括那些要素，而且还包括没有明确列出的其他要素，或者是还包括为这种过程、方法、物品或者设备所固有的要素。在没有更多限制的情况下，由语句“包括一个……”限定的要素，并不排除在包括所述要素的过程、方法、物品或者设备中还存在另外的相同要素。

以上所述仅为本发明的优选实施例而已，并不用于限制本发明，对于本领域的技术人员来说，本发明可以有各种更改和变化。凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。应注意到：相似的标号和字母在下面的附图中表示类似项，因此，一旦某一项在一个附图中被定义，则在随后的附图中不需要对其进行进一步定义和解释。

以上所述，仅为本发明的具体实施方式，但本发明的保护范围并不局限于此，任何熟悉本技术领域的技术人员在本发明揭露的技术范围内，可轻易想到变化或替换，都应涵盖在本发明的保护范围之内。因此，本发明的保护范围应所述以权利要求的保护范围为准。