CN107272066B - 一种含噪地震信号初至走时拾取方法及装置 - Google Patents
一种含噪地震信号初至走时拾取方法及装置 Download PDFInfo
- Publication number
- CN107272066B CN107272066B CN201710482427.3A CN201710482427A CN107272066B CN 107272066 B CN107272066 B CN 107272066B CN 201710482427 A CN201710482427 A CN 201710482427A CN 107272066 B CN107272066 B CN 107272066B
- Authority
- CN
- China
- Prior art keywords
- signal
- seismic
- component
- imf
- noise
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Expired - Fee Related
Links
Classifications
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
- G01V1/00—Seismology; Seismic or acoustic prospecting or detecting
- G01V1/28—Processing seismic data, e.g. analysis, for interpretation, for correction
- G01V1/36—Effecting static or dynamic corrections on records, e.g. correcting spread; Correlating seismic signals; Eliminating effects of unwanted energy
- G01V1/362—Effecting static or dynamic corrections; Stacking
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
- G01V2210/00—Details of seismic processing or analysis
- G01V2210/30—Noise handling
- G01V2210/32—Noise reduction
Abstract
本发明公开了一种含噪地震信号初至走时拾取方法及装置，该方法包括以下步骤：采用自适应噪声的完备集合经验模态分解方法对原始含噪地震信号s(t)分解处理得到一系列从高频到低频的有限个IMF分量及残余分量；根据EIF有效指标函数筛选部分高频IMF分量进行小波阈值去噪；将小波阈值法去噪后的高频IMF分量和不做处理的低频IMF分量以及残余分量累加重构，即得到去噪后的地震信号x(t)。采用改进的滑动时窗能量比法初步确定初至波走时t 0；使用AIC信息准则法计算t 0时刻前后时窗范围内的局部极小值，该极小值对应时间即为精确的初至波走时t。本发明中方法可有效改善地震噪声压制的效果和提高初至波走时的拾取精度，为后续地震资料处理、解释的准确性创造条件。
Description
技术领域
本发明涉及地球物理勘探的地震数据处理技术领域，尤其涉及一种含噪地震信号初至走时拾取方法及装置。
背景技术
地震波在有限区域的介质中的传播需要消耗时间，地震信号由震源(激发点)传播到接收(检波点)的初至时刻数据蕴含着介质内部的物性特征(状态、密度等)。地震记录初至波(P波)走时拾取是地震层析正演计算和静校正技术中一项基础而又重要的工作，对震源定位，油藏描述，探明地质缺陷分布等具有重要意义。初至波走时的手工拾取已不符合现今大规模地震勘探数据的需求。一般的计算机自动拾取方法能满足于地震记录背景较平静且初至波形简单的地震资料。但对于具有较强背景噪声干扰，初至波形变化较大，各种波相互干扰的地震记录，往往难以获得准确的初至时间。因此，在运用自动拾取算法前对地震随机噪声进行压制处理，是提高地震初至波走时拾取精度的有效手段。
发明内容
本发明的目的在于提供一种含噪地震信号初至走时拾取方法及装置，以解决自然或人工地震资料中存在随机噪声干扰时难以获得准确的初至时间的问题，该方法可提高地震信号信噪比，提升地震初至波走时拾取精度。
为达到上述目的，一方面，一种含噪地震信号初至走时拾取方法，其步骤：
S1：采用自适应噪声完备集合经验模态分解(Complementary EnsembleEmpirical Mode Decomposition with Adaptive Noise,CEEMDAN)方法对原始含噪地震信号s(t)分解处理得到一系列从高频到低频的有限个本征模态函数(Intrinsic ModeFunction,IMF)分量及残余分量；
S2：随机噪声往往分布在第一个或前几个高频IMF分量中，采用有效指标函数(Effective Index Function,EIF)确定各阶IMF分量中有效信号与噪声的能量分界点，判断含噪较多的高频IMF分量；
S3：依据能量分界点对部分高频IMF分量进行小波阈值重构小波系数，并保持低频IMF分量及残余分量不变；
S4：将小波阈值中去噪后的高频IMF分量和不做处理的低频IMF分量以及残余分量累加重构，即得到去噪后的地震信号x(t)：
式中，n为CEEMDAN分解得到的IMF分量级数；m为需要进行小波阈值去噪的高频IMF分量数目；IMFi′(t)为进行小波阈值去噪后的高频IMF分量；IMFi(t)为不需要处理的低频IMF分量；Rn(t)为CEEMDAN分解后的残余分量；
S5：采用改进的滑动时窗能量比法识别地震事件并初步确定初至波走时为t0；
S6：去噪后地震序列可截断为多个局部平稳部分，其各部分均可表示为自回归(Auto Recursive,AR)模型。当AR模型阶数固定时，采用AIC信息准则法(AkaikeInformation Criterion,AIC)计算AR模型的AIC函数最小值，该最小值位置，可充当两段地震序列的分界点。在对应t0时刻前后取一时窗，使用AIC函数计算在该时窗范围内的局部极小值，即得到精确的初至波走时t；
S7：根据约束准则对拾取的初至波走时进行评价，对含异常初至时间信息的地震道进行约束处理。
所述原始含噪地震信号的CEEMDAN分解步骤包括：
S1_1，定义算子Ek(·)为产生第k个IMF模态分量的操作过程；wi(t)表示零均值单位方差的正负白噪声；i＝1，2，...，I；系数βk可在每一个模态分解阶段选择合适的信噪比(Signal-Noise Ratio，SNR)；k＝1，2，...，K表示分解阶段；
S1_2，向原始地震信号s(t)中添加零均值单位方差的正负白噪声wi(t)，得到处理后信号序列si(t)为：si(t)＝s(t)+β0wi(t)；
S1_3，对si(t)进行I次验模态分解(Empirical Mode Decomposition，EMD)处理，得到第一阶模态分量IMF1(t)为：
S1_4，对于k＝1，计算第一阶模态分量后残余分量R1(t)为：R1(t)＝s(t)-IMF1(t)；
S1_5，对信号R1(t)+β1E1(w1(t))进行EMD分解，得到第二阶模态分量为：
S1_6，对其余每个阶段，即k＝1，2，...，K，与步骤S1_4的计算过程一致，首先计算第k阶残余分量，再计算第k+1阶模态分量，即有：
S1_7，重复以上步骤直至残余分量成为一个单调函数或整个时间序列上极大值点数和极小值点数之和不大于两个，则结束分解过程。
所述EIF有效指标函数定义如下：
式中，x(t)为去噪后的地震信号，N代表信号长度；IMF(t)为分解出的各级模态分量；EIF越小就表示IMF(t)分量越近似于原始信号，EIF的较大值被认为是噪声。
所述小波阈值去噪步骤包括：首先将能量分界点划分出的含噪高频IMF分量信号转换至小波域，求取小波系数wj：
wj＝W(IMF(t))；
其中W(·)表示小波变换，j为小波分解层数。
设置一个阈值，小于该阈值的小波系数被认为是由噪声产生的，将其去除。阈值处理采用软阈值函数
式中，sgn()为符号函数，wj为小波系数，λ为阈值，，其中σ为噪声标准差，N为信号长度。噪声标准差为小波多分辨级分解系数的中值。
最后对小波系数wj进行逆变换，即得到小波阈值去噪后的高频IMF分量信号。
所述改进的滑动时窗能量比法的原理是：地震记录初至波来临之前，大部分为噪声信号，当初至波来临时初至时间前后时窗内的地震能量特征有较大的差异，该时刻前后时窗能量比值A(t)达到最大值，因此，拾取A(t)最大值对应的时间即可获得地震信号的初至走时。在前、后时窗能量中加入稳定因子aw，增强其稳定性，改进的滑动时窗能量比法的计算公式可表示为：
式中，x(t)为地震记录振幅值，α为稳定系数；为地震道相对能量；T0为前时窗起点，T1为前时窗终点或后时窗起点，T2为后时窗终点。
所述去噪后地震信号的AIC值采用下式计算：
式中，k为有效地震信号和随机噪声分量的最佳分界点，σ1和σ2则分别表示随机噪声和有效地震信号分量的方差，M表示AR模型阶数，N为信号长度，C为常数。
在地震波初至来临时波形发生较大变化，随机噪声和地震信号具有较大的数学统计差异，两者拟合度最差，可求得最小AIC值。因此，计算地震波的AIC极小值可用于地震初至走时的确定。
所述初至拾取异常的道采用如下步骤处理：
根据同一炮的相邻两个检波点的初至时间不会有较大突变，相继各道合理的拾取时间应形成折射段或同相轴的先验知识，对初至走时进行评价，然后构造约束准则完成异常道的初至拾取。计算相邻道的初至时间差绝对值和的平均值。式中m为单炮总道数，ti是每道的初至走时，对相邻初至时间差分进行组合分析，如果满足条件(-1≤n≤9,1≤i≤m,1≤j≤m)，则判定该道初至时间的可信度为1，如果不满足，则判定该道初至时间的可信度为0。确定了部分初至时间拾取可信度为1的道，对筛选后被认为可信度为0的道，通过道间样条插值或者线性插值进行重新确定，即得到异常道纠正处理后的初至走时。其中0＜k≤2，针对信噪比高的地震资料，k参数值可取较大值，反之，k取较小值。
一种含噪地震信号初至走时拾取装置，包括：噪声压制单元U1、初至拾取单元U2。
噪声压制单元，用于压制原始地震信号中的随机噪声，是地震初至波走时准确拾取的前提条件；
初至拾取单元，用于拾取地震波初至来临的时间信息。
噪声压制单元与初至拾取单元通过电路连接。
所述噪声压制单元包括：CEEMDAN分解模块M1、IMF分量分析模块M2、小波阈值去噪模块M3、分量重构模块M4；
CEEMDAN分解模块，用于对原始含噪地震信号进行逐级分解得到一系列从高频到低频的有限个IMF分量以及残余分量；
IMF分量分析模块，采用有效的指标函数EIF来确定能量分界点划分分解得到的有效信号或噪声(大部分的高频噪声会出现在前几个模态分量中)：
式中，N代表信号长度；IMF(t)为原信号分解出的各级模态分量；EIF越小就表示IMF(t)分量越近似于原始信号，EIF的较大值被认为是噪声。
小波阈值去噪模块，将IMF分量分析模块中划分出的含噪高频IMF分量再次进行小波阈值去噪；
分量重构模块，用于将去噪后的高频IMF分量和低频IMF分量以及残余量进行累加重构，即得到去噪后的地震信号。
CEEMDAN分解模块、IMF分量分析模块、小波阈值去噪模块和分量重构模块依次按顺序通过电路连接。
所述初至拾取单元包括：改进滑动时窗能量比法预拾取模块M5、AIC法再拾取模块M6和初至评价与异常处理模块M7；
改进滑动时窗能量比法预拾取模块，用于识别地震事件并大致确定波至时刻，地震记录初至波来临之前，大部分为噪声信号，当初至波来临时初至时间前后时窗内的地震能量特征有较大的差异，该时刻前后时窗能量比值A(t)达到最大值，因此，拾取A(t)最大值对应的时间即可获得地震信号的初至走时。在前、后时窗能量中加入稳定因子aw，增强其稳定性，改进的滑动时窗能量比法的计算公式可表示为：
式中，x(t)为地震记录振幅值，α为稳定系数；为地震道相对能量；T0为前时窗起点，T1为前时窗终点或后时窗起点，T2为后时窗终点。
AIC法再拾取模块，用于初至时刻的精确拾取，在改进滑动时窗能量比法预拾取模块的基础上，于预拾取走时对应时刻前后取一时窗，使用AIC信息准则法计算地震序列AR模型在该时窗范围内的局部极小值，该极小值对应时间即为精确的初至波走时，AIC值计算公式如下式：
式中，k为有效地震信号和随机噪声分量的最佳分界点，σ1和σ2则分别表示随机噪声和有效地震信号分量的方差，M表示AR模型阶数，N为信号长度，C为常数。
在地震波初至来临时波形发生较大变化，随机噪声和地震信号具有较大的数学统计差异，两者拟合度最差，可求得最小AIC值。因此，计算地震波的AIC极小值可用于地震初至走时的确定。
初至评价与异常处理模块，用于处理无信号道或干扰道等异常情况出现时，计算机自动初至拾取方法出错的情况。根据同一炮的相邻两个检波点的初至时间不会有较大突变，相继各道合理的拾取时间应形成折射段或同相轴的先验知识，对初至进行评价，然后构造约束准则完成异常道的初至拾取。计算相邻道的初至时间差绝对值和的平均值。式中m为单炮总道数，ti是每道的初至走时，对相邻初至时间差分进行组合分析，如果满足条件，(-1≤n≤9,1≤i≤m,1≤j≤m)，则判定该道初至时间的可信度为1，如果不满足，则判定该道初至时间的可信度为0。确定了部分初至时间拾取可信度为1的道，对筛选后被认为可信度为0的道，通过道间样条插值或者线性插值进行重新确定，即得到异常道纠正处理后的初至走时。其中0＜k≤2，针对信噪比高的地震资料，k参数值可取较大值，反之，k取较小值。
改进滑动时窗能量比法预拾取模块、AIC法再拾取模块和初至评价与异常处理模块依次按顺序通过电路连接。
本发明的优点：依据先验知识对拾取的初至时刻进行评价，将评价可信度高的道作为可信度低道的约束条件采用插值计算进行二次拾取等一系列操作，免人工干扰，抗干扰能力强，自动剔除低信噪比道，保证初至拾取结果的可靠性和提高拾取速度。本发明可运用于不同复杂度的地震资料，静校正效果好，初至拾取质量高，对实际生产、研究应用具有较高意义。
附图说明
为了更清楚地说明本申请实施例或现有技术中的技术方案，下面将对实施例或现有技术描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本申请中记载的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动性的前提下，还可以根据这些附图获得其他的附图。
附图1是一种含噪地震信号初至走时拾取方法实现流程图；
附图2是一种含噪地震信号初至走时拾取装置的模块结构示意图；
附图3是一种含噪地震信号初至走时拾取方法实现示意图；
附图4是本申请一实施例的合成地震道集数据去噪效果对比图；
附图5是本申请一实施例的合成地震道集数据中各道去噪后信噪比变化对比图；
附图6是本申请一实施例的合成地震道集数据中各道去噪后初至拾取效果对比图；
附图7是本发明一实施例的四组实际单炮地震记录；
附图8是本申请一实施例的四组实际单炮地震记录的初至拾取分布示意图；
附图9是本申请一实施例的四组实际单炮地震记录的初至拾取精度ROC曲线示意图；
具体实施方式
下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有作出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。
本实施例是一种含噪地震信号初至走时拾取方法及装置，见附图1至附图9所示。在本发明实施例中，首先，原始含噪地震信号通过CEEMDAN分解得到各级IMF分量，采用EIF有效指标函数确定有效信号和噪声的能量分界点，对部分高频IMF分量再次进行小波阈值去噪，并保持低频IMF分量和残余分量不变；将去噪后的高频IMF分量和低频IMF分量以及残余量进行累加重构，即得到去噪后的地震信号。然后，采用改进滑动时窗能量比法在地震信号上取一滑动时窗计算前后时窗能量比值，当能量比值大于设定阈值时，即可识别地震事件并大致确定波至时刻；在该时刻前后取一时窗，使用AIC方法计算AR模型在该时窗范围内的局部极小值，即得到精确的初至波走时。最后，循环道集判断是否有初至时间拾取异常的道出现，如有，采用异常处理措施进行纠正处理。
如附图1，本请实施例提供一种含噪地震信号初至走时拾取方法。该方法包括：
S1：采用自适应噪声完备集合经验模态分解CEEMDAN方法对原始含噪地震信号s(t)分解处理得到一系列从高频到低频有限个的本征模态函数IMF分量及残余分量。原始含噪地震信号的CEEMDAN分解步骤包括：
S1_1，定义算子Ek(·)为产生第k个IMF模态分量的操作过程；wi(t)表示零均值单位方差的正负白噪声；i＝1，2，...，I；系数βk可在每一个模态分解阶段选择合适的信噪比SNR；k＝1，2，...，K表示分解阶段；
S1_2，向原始地震信号s(t)中添加零均值单位方差的正负白噪声wi(t)，得到处理后信号序列si(t)为：si(t)＝s(t)+β0wi(t)；
S1_3，对si(t)进行I次验模态分解EMD处理，得到第一阶模态分量IMF1(t)为：
S1_4，对于k＝1，计算第一阶模态分量后残余分量R1(t)为：R1(t)＝s(t)-IMF1(t)；
S1_5，对信号R1(t)+β1E1(w1(t))进行EMD分解，得到第二阶模态分量为：
S1_6，对其余每个阶段，即k＝1，2，...，K，与步骤S1_4的计算过程一致，首先计算第k阶残余分量，再计算第k+1阶模态分量，即有：
S1_7，重复以上步骤直至残余分量成为一个单调函数或整个时间序列上极大值点数和极小值点数之和不大于两个，则结束分解过程。
S2：随机噪声往往分布在第一个或前几个高频IMF分量中，采用有效指标函数EIF确定各阶IMF分量中有效信号与噪声的能量分界点，判断含噪较多的高频数IMF分量。EIF有效指标函数定义如下：
式中，x(t)为去噪后的地震信号，N代表信号长度；IMF(t)为分解出的各级模态分量；EIF越小就表示IMF(t)分量越近似于原始信号，EIF的较大值被认为是噪声。
S3：依据能量分界点对部分高频IMF分量进行小波阈值重构小波系数，并保持低频IMF分量及残余分量不变。小波阈值去噪步骤包括：
首先将能量分界点划分出的含噪高频IMF分量信号转换至小波域，求取小波系数wj：
wj＝W(IMF(t))；
其中W(·)表示小波变换，j为小波分解层数；
设置一个阈值，小于该阈值的小波系数被认为是由噪声产生的，将其去除；阈值处理采用软阈值函数
式中，sgn()为符号函数，wj为小波系数，λ为阈值，，其中σ为噪声标准差，N为信号长度；噪声标准差为小波多分辨级分解系数的中值；
最后对小波系数wj进行逆变换，即得到小波阈值去噪后的高频IMF分量信号。
S4：将小波阈值中去噪后的高频IMF分量和不做处理的低频IMF分量以及残余分量累加重构，即得到去噪后的地震信号x(t)：
式中，n为CEEMDAN分解得到的IMF分量级数；m为需要进行小波阈值去噪的高频IMF分量数目；IMFi′(t)为进行小波阈值去噪后的高频IMF分量；IMFi(t)为不需要处理的低频IMF分量；Rn(t)为CEEMDAN分解后的残余分量；
S5：采用改进的滑动时窗能量比法识别地震事件并初步确定初至波走时为t0。改进的滑动时窗能量比法的原理是根据地震记录初至波到达时，该时刻的前后时窗能量比值A(t)达到最大值，拾取A(t)最大值对应的时刻即为初至走时。在前、后时窗能量中加入稳定因子，增强其稳定性，改进公式可表示为：
式中，x(t)为地震记录振幅值，α为稳定系数；为地震道相对能量；T1为第1个时窗起点，T0为第1个时窗终点或第2个时窗起点，T2为第2个时窗终点。
S6：去噪后地震序列可截断为多个局部平稳部分，其各部分均可表示为自回归AR模型；当AR模型阶数固定时，采用AIC信息准则法AIC计算AR模型的AIC函数最小值，该最小值位置，可充当两段地震序列的分界点；在对应t0时刻前后取一时窗，使用AIC函数计算在该时窗范围内的局部极小值，即得到精确的初至波走时t。所述去噪后地震信号的AIC值采用下式计算：
式中，k为有效地震信号和随机噪声分量的最佳分界点，σ1和σ2则分别表示随机噪声和有效地震信号分量的方差，M表示AR模型阶数，N为信号长度，C为常数；
在地震波初至来临时波形发生较大变化，随机噪声和地震信号具有较大的数学统计差异，两者拟合度最差，可求得最小AIC值；因此，计算地震波的AIC极小值可用于地震初至走时的确定。
S7：根据约束准则对拾取的初至波走时进行评价，对含异常初至时间信息的地震道进行约束处理。根据同一炮的相邻两个检波点的初至时间不会有较大突变，相继各道合理的拾取时间应形成折射段或同相轴的先验知识，对初至走时进行评价，然后构造约束准则完成异常道的初至拾取；计算相邻道的初至时间差绝对值和的平均值；式中m为单炮总道数，ti是每道的初至走时，对相邻初至时间差分进行组合分析，如果满足条件，(-1≤n≤9,1≤i≤m,1≤j≤m)，则判定该道初至时间的可信度为1，如果不满足，则判定该道初至时间的可信度为0；确定了部分初至时间拾取可信度为1的道，对筛选后被认为可信度为0的道，通过道间样条插值或者线性插值进行重新确定，即得到异常道纠正处理后的初至走时。其中0＜k≤2，针对信噪比高的地震资料，k参数值可取较大值，反之，k取较小值。
附图3为一种含噪地震信号初至走时拾取方法实现示意图，图3(a)为原始含噪单道地震波形，其信噪比较低(9dB)，初至拾取困难；图3(b)为本申请实施例中CEEMDAN联合小波阈值去噪后地震波形，信噪比提升明显(18dB)，大部分随机噪声被去除，人工可准确拾取地震初至时刻为492ms；图3(c)为本申请实施例中改进滑动时窗能量比法大致确定的地震波初至时间为486ms；图3(d)为本申请实施例中AIC信息准则在486ms前后取一时窗(时窗长度为200ms)计算时窗局部范围内的最小AIC值对应走时为492ms，与人工拾取一致。
附图4为合成地震道集数据去噪效果对比图，图4(a)为纯净的单炮记录(基于Ricker子波人工合成一组单同相轴地震数据，采样点数为1000，道数为40，波速为3000m/s，道间距为50m，采样频率为1000Hz)，图4(b)为加入0db高斯白噪声的含噪记录，图4(c)为EMD去噪结果，图4(d)为本申请实施例中CEEMDAN联合小波阈值去噪结果。从图4(c)中可以发现EMD去噪后的记录可以有效地将淹没在噪声中的同相轴恢复出来，具有一定去噪能力，但一些较小的突变仍存在，模态混叠导致EMD去噪过程不彻底；从图4(d)可以看出本申请实施例的去噪效果要优于EMD方法，有效地平滑了EMD方法中的高频毛刺，使得去噪后的结果与纯净信号更为接近。
附图5为合成单同相轴地震信号采用不同去噪方法处理后各道SNR分布对比图，由图可知，采用本申请实施例的方法去噪后信噪比提升最为明显。表1为各方法去噪后整体道集信噪比SNR及均方根误差RMSE的数值对比，较其余方法SNR整体提升明显，重构去噪后RMSE最小。由以上图表可知本申请实施例中联合去噪方法可有效去除地震信号中的随机噪声干扰，提升地震信号信噪比，降低去噪信号的重构误差。
表1合成地震数据各方法去噪效果对比
其中，信噪比(SNR)和均方根误差(RMSE)采用如下公式计算：
式中，si为原始地震信号，xi为去噪后地震信号，N代表信号长度。
附图6是本申请一实施例的合成地震道集数据中各道去噪后初至拾取效果对比图；表2为对应拾取均方根误差RMSE及其在相应绝对拾取误差范围内的占比情况。由图6和表2数据可知，本申请实施例中CEEMDAN联合小波阈值法去除地震随机噪声后再采用改进滑动时窗能量比联合AIC信息准则对去噪后地震数据进行初至拾取，其拾取误差较小，各道初至时间分布与纯净信号经人工拾取相当，拾取误差明显低于未去噪直接采用其他初至拾取组合模式。本申请实施例中联合去噪方法与本申请实施例中联合初至拾取方法的组合，在50ms内的拾取误差范围内，全部准确拾取；在10ms拾取误差范围内，达95％准确率；在5ms拾取误差范围内，达87％准确率，明显优于其它组合模式，验证了本发明方法的有效性。
表2对应拾取误差参数对比
附图7是本申请实施例中四组实际单炮地震记录，其中(a)和(b)为炸药震源地震波形，(c)和(d)为可控源地震波形。图中M为炮号，N为道数，dt为采样时间间隔，dx为道间距。
附图8为对应四组单炮记录采用不同初至拾取方法的拾取结果分布图。由图8可知本申请实施例中联合去噪方法与本申请实施例中联合初至拾取方法组合的拾取结果与人工拾取相当，相比其他去噪和拾取方法组合具有较低的拾取误差；在(c)和(d)图中，由于无信号道或干扰道等异常情况的出现，导致采用短时窗平均/长时窗平均初至拾取方法(LTA/STA)时相应道出现较大拾取误差的情况。
在本申请实施例中，根据同一炮的相邻两个检波点的初至时间不会有较大突变，相继各道合理的拾取时间应形成折射段或同相轴的先验知识，对初至走时进行评价，然后构造约束准则完成异常道的初至拾取；计算相邻道的初至时间差绝对值和的平均值；式中m为单炮总道数，ti是每道的初至走时，对相邻初至时间差分进行组合分析，如果满足条件，(-1≤n≤9,1≤i≤m,1≤j≤m)，则判定该道初至时间的可信度为1，如果不满足，则判定该道初至时间的可信度为0；确定了部分初至时间拾取可信度为1的道，对筛选后被认为可信度为0的道，通过道间样条插值或者线性插值进行重新确定，即得到异常道纠正处理后的初至走时。其中0＜k≤2，针对信噪比高的地震资料，k参数值可取较大值，反之，k取较小值。采用本申请实施例中联合初至拾取方法可有效避免较大拾取误差问题的发生。
图9为对应拾取精度ROC曲线图，ROC曲线中不需要设定各方法拾取与人工拾取误差阈值即可对算法本身效果做出评价。其中ROC曲线右下方的面积被定义为AUC值，其值越接近1，表示该算法准确性越高。四组实际地震记录中采用本申请实施例中联合去噪方法与本申请实施例中联合初至拾取方法组合的ROC曲线AUC平均值达0.96，明显高于其余组合。综合以上分析可知本申请实施例中联合去噪方法与联合初至拾取方法的组合模式，可取得最好的初至拾取效果。
如附图2，本请实施例还提供一种含噪地震信号初至走时拾取装置，包括：噪声压制单元U1，初至拾取单元U2，噪声压制单元U1与初至拾取单元U2通过电路连接。噪声压制单元U1中包括CEEMDAN分解模块M1、IMF分量分析模块M2、小波阈值去噪模块M3、分量重构模块M4，M1-M4依次按顺序通过电路连接；初至拾取单元U2中包括改进滑动时窗能量比法预拾取模块M5、AIC法再拾取模块M6和初至评价与异常处理模块M7，M5-M7依次按顺序通过电路连接。
噪声压制单元U1，用于压制原始地震信号中的随机噪声，是地震初至波走时准确拾取的前提条件；
初至拾取单元U2，用于拾取地震波初至来临的时间信息。
所述噪声压制单元U1包括：CEEMDAN分解模块M1、IMF分量分析模块M2、小波阈值去噪模块M3、分量重构模块M4。
CEEMDAN分解模块M1，用于对原始含噪地震信号进行逐级分解得到频谱分布由高至低的各级IMF分量以及残余分量；
IMF分量分析模块M2，采用有效的指标函数(EIF)来确定能量分界点划分分解得到的有效信号或噪声。
小波阈值去噪模块M3，将IMF分量分析模块M2中划分出的含噪高频IMF分量再次进行小波阈值去噪；
分量重构模块M4，用于将去噪后的高频IMF分量和低频IMF分量以及残余量进行累加重构，即得到去噪后的地震信号。
所述初至拾取单元U2包括：改进滑动时窗能量比法预拾取模块M5、AIC
法再拾取模块M6和初至评价与异常处理模块M7。
改进滑动时窗能量比法预拾取模块M5，用于识别地震事件并大致确定波至时刻，当初至来临时，该时刻的前后时窗能量比值A(t)将达到最大值，拾取A(t)最大值对应的时间即为初至时间。
AIC法再拾取模块M6，用于初至时刻的精确拾取，在改进滑动时窗能量比法预拾取模块M5的基础上，于预拾取时刻前后取一时窗，使用AIC方法计算地震序列AR模型在该时窗范围内的局部极小值，即得到精确的初至波走时。
初至评价与异常处理模块M7，用于处理无信号道或干扰道等异常情况出现时，计算机自动初至拾取方法出错的情况。
本申请实施例一种含噪地震信号初至走时拾取装置的各组成部分分别用于实现前述实施例方法的各步骤，在前述方法实施例中，已经对各步骤进行了详细说明，在此不再赘述。
在一个或多个示例性的设计中，本申请实施例所描述的上述功能可以在硬件、软件、固件或这三者的任意组合来实现。如果在软件中实现，这些功能可以存储与电脑可读的媒介上，或以一个或多个指令或代码形式传输于电脑可读的媒介上。电脑可读媒介包括电脑存储媒介和便于使得让电脑程序从一个地方转移到其它地方的通信媒介。存储媒介可以是任何通用或特殊电脑可以接入访问的可用媒体。例如，这样的电脑可读媒体可以包括但不限于RAM、ROM、EEPROM、CD-ROM或其它光盘存储、磁盘存储或其它磁性存储装置，或其它任何可以用于承载或存储以指令或数据结构和其它可被通用或特殊电脑、或通用或特殊处理器读取形式的程序代码的媒介。
以上所述的具体实施例，对本申请的目的、技术方案和有益效果进行了进一步详细说明，所应理解的是，以上所述仅为本申请实施例的具体实施例而已，并不用于限定本申请的保护范围，凡在本申请的精神和原则之内，所做的任何修改、等同替换、改进等，均应包含在本申请的保护范围之内。
Claims (5)
1.一种含噪地震信号初至走时拾取方法，其步骤：
S1：采用自适应噪声完备集合经验模态分解CEEMDAN方法对原始含噪地震信号s(t)分解处理得到一系列从高频到低频有限个的本征模态函数IMF分量及残余分量；
S2：随机噪声分布在第一个或前几个高频IMF分量中，采用有效指标函数确定各阶IMF分量中有效信号与噪声的能量分界点，判断含噪较多的高频IMF分量；所述有效指标函数定义如下：
式中，x(t)为去噪后的地震信号，N代表信号长度；IMF(t)为分解出的各级模态分量；EIF表示有效指标函数的值，其值越小就表示IMF(t)分量越近似于原始信号，EIF的较大值被认为是噪声；
S3：依据能量分界点对部分高频IMF分量进行小波阈值重构小波系数，并保持低频IMF分量及残余分量不变；
S4：将小波阈值中去噪后的高频IMF分量和不做处理的低频IMF分量以及残余分量累加重构，即得到去噪后的地震信号x(t)：
式中，n为CEEMDAN分解得到的IMF分量级数；m为需要进行小波阈值去噪的高频IMF分量数目；IMFi′(t)为进行小波阈值去噪后的高频IMF分量；IMFi(t)为不需要处理的低频IMF分量；Rn(t)为CEEMDAN分解后的残余分量；
S5：采用改进的滑动时窗能量比法识别地震事件并初步确定初至波走时为t0；所述改进的滑动时窗能量比法的原理是：地震记录初至波来临之前，大部分为噪声信号，当初至波来临时初至时间前后时窗内的地震能量特征有较大的差异，该时刻前后时窗能量比值A(t)达到最大值，因此，拾取A(t)最大值对应的时间即可获得地震信号的初至走时；在前、后时窗能量中加入稳定因子aw，增强其稳定性，改进的滑动时窗能量比法的计算公式可表示为：
式中，x(t)为地震记录振幅值，α为稳定系数；为地震道相对能量，其中N为地震信号采样点数，T为采样时间；T0为前时窗起点，T1为前时窗终点或后时窗起点，T2为后时窗终点；
S6：去噪后地震序列可截断为多个局部平稳部分，其各部分均可表示为自回归AR模型；当AR模型阶数固定时，采用AIC信息准则法AIC计算AR模型的AIC函数最小值，该最小值位置，可充当两段地震序列的分界点；在对应t0时刻前后取一时窗，使用AIC函数计算在该时窗范围内的局部极小值，即得到精确的初至波走时t；
S7：根据同一炮的相邻两个检波点的初至时间不会有较大突变，相继各道合理的拾取时间应形成折射段或同相轴的先验知识，对初至走时进行评价，然后构造约束准则完成异常道的初至拾取。
2.根据权利要求1所述的一种含噪地震信号初至走时拾取方法，其特征在于：所述原始含噪地震信号的CEEMDAN分解步骤包括：
S1_1，定义算子Ek(·)为产生第k个IMF模态分量的操作过程；wi(t)表示零均值单位方差的正负白噪声；i＝1，2，...，I；系数βk可在每一个模态分解阶段选择合适的信噪比SNR；k＝1，2，...，K表示分解阶段；
S1_2，向原始含噪地震信号s(t)中添加零均值单位方差的正负白噪声wi(t)，得到处理后信号序列si(t)为：si(t)＝s(t)+β0wi(t)，其中β0为第一次添加正负白噪声wi(t)的幅值；
S1_3，对si(t)进行I次经验模态分解EMD处理，得到第一阶模态分量IMF1(t)为：
S1_4，对于k＝1，计算第一阶模态分量后残余分量R1(t)为：R1(t)＝s(t)-IMF1(t)；
S1_5，对信号R1(t)+β1E1(w1(t))进行EMD分解，得到第二阶模态分量为：
S1_6，对其余每个阶段，即k＝2，...，K，与步骤S1_4的计算过程一致，首先计算第k阶残余分量，再计算第k+1阶模态分量，即有：
Rk(t)＝Rk-1(t)-IMFk(t)，
S1_7，重复以上步骤直至残余分量成为一个单调函数或整个时间序列上极大值点数和极小值点数之和不大于两个，则结束分解过程。
3.根据权利要求1所述的一种含噪地震信号初至走时拾取方法，其特征在于：所述小波阈值去噪步骤包括，首先将能量分界点划分出的含噪高频IMF分量信号转换至小波域，求取小波系数wj：
wj＝W(IMF(t))；
其中W(·)表示小波变换，j为小波分解层数；
设置一个阈值，小于该阈值的小波系数被认为是由噪声产生的，将其去除；阈值处理采用软阈值函数
式中，sgn()为符号函数，wj为小波系数，λ为阈值，其中σ为噪声标准差，N为信号长度；噪声标准差Median(|wj|)为小波多分辨级分解系数的中值；
最后对小波系数wj进行逆变换，即得到小波阈值去噪后的高频IMF分量信号。
4.根据权利要求1所述的一种含噪地震信号初至走时拾取方法，其特征在于：所述去噪后地震信号的AIC值采用下式计算：
式中，k为有效地震信号和随机噪声分量的最佳分界点，σ1和σ2则分别表示随机噪声和有效地震信号分量的方差，M表示AR模型阶数，N为信号长度，C为常数；
在地震波初至来临时波形发生较大变化，随机噪声和地震信号具有较大的数学统计差异，两者拟合度最差，可求得最小AIC值；因此，计算地震波的AIC极小值可用于地震初至走时的确定。
5.根据权利要求1所述的一种含噪地震信号初至走时拾取方法，其特征在于：所述初至拾取异常的道采用如下步骤处理：
根据同一炮的相邻两个检波点的初至时间不会有较大突变，相继各道合理的拾取时间应形成折射段或同相轴的先验知识，对初至走时进行评价，然后构造约束准则完成异常道的初至拾取；计算相邻道的初至时间差绝对值和的平均值 式中m为单炮总道数，ti是每道的初至走时，对相邻初至时间差分进行组合分析，如果满足条件当满足-1≦n≦9,1≦i≦m,1≦j≦m，则判定该道初至时间的可信度为1，如果不满足，则判定该道初至时间的可信度为0；确定了部分初至时间拾取可信度为1的道，对筛选后被认为可信度为0的道，通过道间样条插值或者线性插值进行重新确定，即得到异常道纠正处理后的初至走时；其中0＜k≤2，针对信噪比高的地震资料，k参数值可取较大值，反之，k取较小值。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201710482427.3A CN107272066B (zh) | 2017-06-22 | 2017-06-22 | 一种含噪地震信号初至走时拾取方法及装置 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201710482427.3A CN107272066B (zh) | 2017-06-22 | 2017-06-22 | 一种含噪地震信号初至走时拾取方法及装置 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN107272066A CN107272066A (zh) | 2017-10-20 |
| CN107272066B true CN107272066B (zh) | 2019-01-25 |
Family
ID=60069167
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN201710482427.3A Expired - Fee Related CN107272066B (zh) | 2017-06-22 | 2017-06-22 | 一种含噪地震信号初至走时拾取方法及装置 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN107272066B (zh) |
Families Citing this family (31)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN107807387B (zh) * | 2017-10-31 | 2019-08-27 | 中国科学技术大学 | 基于神经网络的地震初至波走时获取方法 |
| CN110031899B (zh) * | 2018-01-11 | 2020-09-08 | 中国石油化工集团有限公司 | 基于压缩感知的弱信号提取算法 |
| CN108267784A (zh) * | 2018-01-18 | 2018-07-10 | 国网江西省电力有限公司电力科学研究院 | 一种地震信号随机噪声压制处理方法 |
| CN108345033B (zh) * | 2018-01-26 | 2019-07-26 | 中国石油大学(华东) | 一种微地震信号时频域初至检测方法 |
| CN108387887A (zh) * | 2018-05-22 | 2018-08-10 | 西安邮电大学 | 一种水声信号的混合降噪方法 |
| CN109031422A (zh) * | 2018-08-09 | 2018-12-18 | 吉林大学 | 一种基于CEEMDAN与Savitzky-Golay滤波的地震信号噪声抑制方法 |
| CN110954941B (zh) * | 2018-09-26 | 2021-08-24 | 中国石油化工股份有限公司 | 一种自动化初至拾取方法及系统 |
| CN109871733B (zh) * | 2018-09-27 | 2023-05-19 | 南京信息工程大学 | 一种自适应海杂波信号去噪方法 |
| CN109239554A (zh) * | 2018-09-28 | 2019-01-18 | 山东康威通信技术股份有限公司 | 一种电力电缆局部放电信号去噪及有效信号提取方法及系统 |
| CN110967739B (zh) * | 2018-09-30 | 2021-11-05 | 中国石油化工股份有限公司 | 基于误差正态分布的微地震识别质量分析方法及系统 |
| CN109730670B (zh) * | 2018-10-12 | 2021-11-30 | 浙江大学宁波理工学院 | 一种基于异构计算的心电信号降噪方法 |
| CN111123356B (zh) * | 2018-10-30 | 2022-07-12 | 中石化石油工程技术服务有限公司 | 基于初至信息的异常道智能判识方法 |
| CN109343118B (zh) * | 2018-11-07 | 2019-12-06 | 西南石油大学 | 一种异常初至时间修正方法 |
| CN109283576B (zh) * | 2018-11-26 | 2020-06-09 | 辽宁工程技术大学 | 一种以振幅为特征函数的自动拾取p波震相方法 |
| CN110007342B (zh) * | 2019-04-03 | 2020-12-15 | 中煤科工集团西安研究院有限公司 | 一种用于低信噪比地震信号的时频域直接拾取初至方法及系统 |
| CN110147766B (zh) * | 2019-05-21 | 2022-06-03 | 东华理工大学 | 基于移不变稀疏编码的低频大地电磁信号去噪方法 |
| CN110297271B (zh) * | 2019-06-26 | 2020-09-11 | 中国矿业大学 | 一种用于矿震报警的单分量探头p波初至到时修正方法 |
| CN110619265A (zh) * | 2019-07-31 | 2019-12-27 | 江西理工大学 | 球磨机筒体振动信号联合去噪方法、装置及存储介质 |
| CN112711074B (zh) * | 2019-10-24 | 2024-03-26 | 中国石油化工股份有限公司 | 一种地震初至波的去噪方法及装置 |
| CN111538082B (zh) * | 2020-06-05 | 2021-12-07 | 吉林大学 | 一种地震波时频域初至自动拾取方法 |
| CN111637938A (zh) * | 2020-06-10 | 2020-09-08 | 南京农业大学 | 一种冲量式谷物质量流量信号降噪方法 |
| CN112200037B (zh) * | 2020-09-29 | 2023-12-26 | 中国科学院上海微系统与信息技术研究所 | 一种微弱信号检测方法、终端和计算机可读存储介质 |
| CN112505781A (zh) * | 2020-10-28 | 2021-03-16 | 中国石油天然气集团有限公司 | 用于地震采集初至拾取的图像处理方法及装置 |
| CN112558159A (zh) * | 2020-12-08 | 2021-03-26 | 中国石油天然气集团有限公司 | 一种声波测井初至拾取方法 |
| CN113297195B (zh) * | 2021-07-28 | 2022-03-15 | 云智慧（北京）科技有限公司 | 一种时间序列异常检测方法、装置及设备 |
| CN113985478A (zh) * | 2021-07-28 | 2022-01-28 | 中石化石油工程技术服务有限公司 | 一种低信噪比地震资料初至自动拾取及修正的方法 |
| CN113723207A (zh) * | 2021-08-03 | 2021-11-30 | 上海海事大学 | 一种基于直方图距离的突变信号检测方法 |
| CN114114416A (zh) * | 2021-11-26 | 2022-03-01 | 同济大学 | 一种基于强噪音弱信号检测的初至拾取方法 |
| CN114486262B (zh) * | 2022-01-28 | 2023-03-21 | 河海大学 | 一种基于cnn-at-lstm的轴承剩余使用寿命预测方法 |
| CN116071881B (zh) * | 2023-03-24 | 2023-06-13 | 光谷技术有限公司 | 基于地波和视频的多模态入侵探测系统 |
| CN117148432B (zh) * | 2023-10-27 | 2024-03-19 | 胜利信科(山东)勘察测绘有限公司 | 基于多尺度分量提取的浅剖数据空间插值方法 |
Citations (5)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN104459779A (zh) * | 2014-12-09 | 2015-03-25 | 中国石油天然气集团公司 | 一种异常地震道的自动识别的方法和装置 |
| CN104636609A (zh) * | 2015-01-30 | 2015-05-20 | 电子科技大学 | 一种基于经验模态分解与小波分析的信号联合去噪方法 |
| WO2017040399A1 (en) * | 2015-09-04 | 2017-03-09 | Saudi Arabian Oil Company | Improvement and automatic quality control of seismic travel time |
| CN106646598A (zh) * | 2016-12-27 | 2017-05-10 | 吉林大学 | 一种fast‑aic法微地震信号拾取方法 |
| CN106798554A (zh) * | 2017-01-12 | 2017-06-06 | 安徽大学 | 一种含噪imf分量及心电信号的去噪方法 |
-
2017
- 2017-06-22 CN CN201710482427.3A patent/CN107272066B/zh not_active Expired - Fee Related
Patent Citations (5)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN104459779A (zh) * | 2014-12-09 | 2015-03-25 | 中国石油天然气集团公司 | 一种异常地震道的自动识别的方法和装置 |
| CN104636609A (zh) * | 2015-01-30 | 2015-05-20 | 电子科技大学 | 一种基于经验模态分解与小波分析的信号联合去噪方法 |
| WO2017040399A1 (en) * | 2015-09-04 | 2017-03-09 | Saudi Arabian Oil Company | Improvement and automatic quality control of seismic travel time |
| CN106646598A (zh) * | 2016-12-27 | 2017-05-10 | 吉林大学 | 一种fast‑aic法微地震信号拾取方法 |
| CN106798554A (zh) * | 2017-01-12 | 2017-06-06 | 安徽大学 | 一种含噪imf分量及心电信号的去噪方法 |
Non-Patent Citations (1)
| Title |
|---|
| 基于经验模态分解的地震P波初至自动识别研究;黄汉明 等;《Proceedings of 2010 The 3rd International Conference on Computational Intelligence and Industrial Application》;20101231;第7卷;第107-108页 |
Also Published As
| Publication number | Publication date |
|---|---|
| CN107272066A (zh) | 2017-10-20 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN107272066B (zh) | 一种含噪地震信号初至走时拾取方法及装置 | |
| CN108549106B (zh) | 混叠噪声压制方法及装置 | |
| WO2008112036A1 (en) | Imaging of multishot seismic data | |
| CN101598812B (zh) | 去除数字检波器单点接收地震记录中的异常噪声方法 | |
| CN108267784A (zh) | 一种地震信号随机噪声压制处理方法 | |
| US10317549B2 (en) | Systems and methods for non-parametric autopicking of seismic wave features from seismic data | |
| CN107144879A (zh) | 一种基于自适应滤波与小波变换结合的地震波降噪方法 | |
| CN105607125A (zh) | 基于块匹配算法和奇异值分解的地震资料噪声压制方法 | |
| CN110780348B (zh) | 基于立体成像条件的一次波和多次波联合成像方法及系统 | |
| CN108345033A (zh) | 一种微地震信号时频域初至检测方法 | |
| CN109557587B (zh) | 一种vsp地震资料井筒波频率域滤波方法及装置 | |
| CN109425897B (zh) | 消除地震资料野值干扰的方法及系统 | |
| CN109507726A (zh) | 时间域弹性波多参数全波形的反演方法及系统 | |
| CN102854530B (zh) | 基于对数时频域双曲平滑的动态反褶积方法 | |
| CN109884691B (zh) | 用于随采地震信号的强单频和随机噪声压制方法及系统 | |
| CN108037533B (zh) | 一种基于块重组期望对数似然的地震勘探噪声压制方法 | |
| CN106950597A (zh) | 基于三边滤波的混合震源数据分离方法 | |
| CN110032968A (zh) | 基于双树复小波和自适应半软阈值法的去噪方法 | |
| CN110007342B (zh) | 一种用于低信噪比地震信号的时频域直接拾取初至方法及系统 | |
| CN111276154B (zh) | 风噪声抑制方法与系统以及炮声检测方法与系统 | |
| CN111538086B (zh) | 一种提高地震数据初至波质量的初至自动拾取方法 | |
| CN110850443B (zh) | 激光雷达温湿度数据阶梯分析处理方法 | |
| CN112327356A (zh) | 基于同相轴迭代追踪提取的混叠记录分离方法 | |
| CN114779334B (zh) | 一种基于统计理论模型的地表一致性振幅处理方法 | |
| CN112200069A (zh) | 时频域谱减法和经验模态分解联合的隧道滤波方法及系统 |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| PB01 | Publication | ||
| PB01 | Publication | ||
| SE01 | Entry into force of request for substantive examination | ||
| SE01 | Entry into force of request for substantive examination | ||
| GR01 | Patent grant | ||
| GR01 | Patent grant | ||
| CF01 | Termination of patent right due to non-payment of annual fee | ||
| CF01 | Termination of patent right due to non-payment of annual fee |
Granted publication date: 20190125 Termination date: 20190622 |