CN107966734B - 多分量地震数据的矢量去噪方法 - Google Patents
多分量地震数据的矢量去噪方法 Download PDFInfo
- Publication number
- CN107966734B CN107966734B CN201710867404.4A CN201710867404A CN107966734B CN 107966734 B CN107966734 B CN 107966734B CN 201710867404 A CN201710867404 A CN 201710867404A CN 107966734 B CN107966734 B CN 107966734B
- Authority
- CN
- China
- Prior art keywords
- vector
- wave vector
- wave
- window
- component
- 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.)
- Active
Links
- 239000013598 vector Substances 0.000 title claims abstract description 205
- 238000000034 method Methods 0.000 title claims abstract description 28
- 238000001914 filtration Methods 0.000 claims abstract description 54
- 238000005070 sampling Methods 0.000 claims description 2
- 230000002123 temporal effect Effects 0.000 claims description 2
- 230000010287 polarization Effects 0.000 description 25
- 230000015572 biosynthetic process Effects 0.000 description 23
- 238000010586 diagram Methods 0.000 description 22
- 238000003786 synthesis reaction Methods 0.000 description 20
- 238000001228 spectrum Methods 0.000 description 17
- 238000000926 separation method Methods 0.000 description 11
- 230000001427 coherent effect Effects 0.000 description 9
- 239000011159 matrix material Substances 0.000 description 5
- NAWXUBYGYWOOIX-SFHVURJKSA-N (2s)-2-[[4-[2-(2,4-diaminoquinazolin-6-yl)ethyl]benzoyl]amino]-4-methylidenepentanedioic acid Chemical compound C1=CC2=NC(N)=NC(N)=C2C=C1CCC1=CC=C(C(=O)N[C@@H](CC(=C)C(O)=O)C(O)=O)C=C1 NAWXUBYGYWOOIX-SFHVURJKSA-N 0.000 description 4
- 238000009499 grossing Methods 0.000 description 4
- 239000002245 particle Substances 0.000 description 4
- 238000012545 processing Methods 0.000 description 4
- 238000004458 analytical method Methods 0.000 description 2
- 230000008859 change Effects 0.000 description 2
- 230000000694 effects Effects 0.000 description 2
- 238000005516 engineering process Methods 0.000 description 2
- 238000005259 measurement Methods 0.000 description 2
- 238000012360 testing method Methods 0.000 description 2
- 150000001875 compounds Chemical class 0.000 description 1
- 239000012141 concentrate Substances 0.000 description 1
- 238000001514 detection method Methods 0.000 description 1
- 239000006185 dispersion Substances 0.000 description 1
- 239000012530 fluid Substances 0.000 description 1
- 230000004927 fusion Effects 0.000 description 1
- 230000006872 improvement Effects 0.000 description 1
- 238000012804 iterative process Methods 0.000 description 1
- 239000000203 mixture Substances 0.000 description 1
- 238000012986 modification Methods 0.000 description 1
- 230000004048 modification Effects 0.000 description 1
- 238000003969 polarography Methods 0.000 description 1
- 230000008569 process Effects 0.000 description 1
- 230000000644 propagated effect Effects 0.000 description 1
- 238000000746 purification Methods 0.000 description 1
- 238000011084 recovery Methods 0.000 description 1
- 238000011160 research Methods 0.000 description 1
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. for interpretation or for event detection
- G01V1/36—Effecting static or dynamic corrections on records, e.g. correcting spread; Correlating seismic signals; Eliminating effects of unwanted energy
- G01V1/364—Seismic filtering
-
- 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. for interpretation or for event detection
- G01V1/36—Effecting static or dynamic corrections on records, e.g. correcting spread; Correlating seismic signals; Eliminating effects of unwanted energy
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
- G01V1/00—Seismology; Seismic or acoustic prospecting or detecting
- G01V1/16—Receiving elements for seismic signals; Arrangements or adaptations of receiving elements
- G01V1/18—Receiving elements, e.g. seismometer, geophone or torque detectors, for localised single point measurements
- G01V1/181—Geophones
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
- G01V2210/00—Details of seismic processing or analysis
- G01V2210/20—Trace signal pre-filtering to select, remove or transform specific events or signal components, i.e. trace-in/trace-out
- G01V2210/24—Multi-trace filtering
- G01V2210/242—F-k filtering, e.g. ground roll
-
- 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
-
- 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
- G01V2210/324—Filtering
Landscapes
- Engineering & Computer Science (AREA)
- Remote Sensing (AREA)
- Physics & Mathematics (AREA)
- Life Sciences & Earth Sciences (AREA)
- Acoustics & Sound (AREA)
- Environmental & Geological Engineering (AREA)
- Geology (AREA)
- General Life Sciences & Earth Sciences (AREA)
- General Physics & Mathematics (AREA)
- Geophysics (AREA)
- Geophysics And Detection Of Objects (AREA)
Abstract
本发明公开了一种多分量地震数据的矢量去噪方法，包括：取得多分量地震数据的波矢量，其中多分量地震数据的波矢量包括有效信号、地滚波和随机噪音的波矢量。利用第一移动时窗，对多分量地震数据的波矢量进行计算出第一平均波矢量、并将第一平均波矢量进行中值滤波，以取得地滚波的真实模量，且将多分量地震数据的波矢量与地滚波的波矢量相减，以取得矢量时间序列。利用第二移动时窗，对所述矢量时间序列进行均值滤波，以取得第二平均波矢量，并对所述第二平均波矢量进行中值滤波以及在同一接收线内进行中值滤波，以压制矢量时间序列中的随机噪音的波矢量，进而取得净化后的有效信号的波矢量。
Description
技术领域
本发明涉及地震勘探的技术领域，尤其涉及一种多分量地震数据的矢量去噪方法。
背景技术
多分量地震检波器可以充分记录弹性波场的空间质点运动，充分利用弹性波场的矢量振幅信息有利于进行随后的储层预测和流体检测。因此，在多分量数据处理过程中保持弹性波场的波矢信息则显得尤为重要，为了避免有效地震信号矢量保真度的损失，多分量地震采集时通常采用单点非组合接收模式，这也将导致获得的多分量地震数据出现较低的信噪比(S/N)。
许多现有的去噪方法可以直接应用于多分量地震数据来提高信噪比，但其中大部分是针对单分量的，而且去噪方法的应用受到有效信号的能量和信噪比的影响，在大多数情况下，垂直Z分量具有最高的信噪比，水平分量相对较差。因此，滤波参数的变化以及每个分量上信号和噪声的变化，导致传统的去噪方法在某种程度上不足以维持矢量振幅信息。所以，在去噪过程中保持矢量振幅信息对于多分量地震数据处理来说仍然是一个挑战。
极化滤波是利用三分量记录测量的矢量振幅完全分离有效信号与噪声的方法之一。极化滤波方法已经被研究了很多年，大多数研究是基于某种椭圆率和方向性的带通滤波器，在这类滤波方法中，通常假设地滚波倾向于在垂直平面中出现椭圆极化，体波线性极化，随机噪声则极化很差，基于这样预定义的极化特性(例如椭圆率或极化方向)，将有效信号与噪音背景进行分离。
极化特性通常在时域或频域(或时间-频率)中获得，在时域中，大多数技术是在给定的两分量或三分量数据的移动时窗内，对构建的数据矩阵进行特征分析来实现的。可以根据数据协方差矩阵的特征值表示的偏振轨迹的半长轴和半短轴，或通过数据矩阵上的SVD的本征图来确定极化类型。这些技术的成功取决于主周期，信噪比和极化类型影响的移动窗口的选择。目前，时域极化滤波所面临的问题主要是难以对同一时间到达同一时窗的有效信号和噪音进行分离。
发明内容
本发明的主要目的在于提供一种多分量地震数据的矢量去噪方法，以解决现有技术存在难以对同一时间到达同一时窗的有效信号和噪音进行分离的问题。
为解决上述问题，本发明实施例提供一种多分量地震数据的矢量去噪方法，包括：取得多分量地震数据的波矢量，其中多分量地震数据的波矢量包括有效信号、地滚波和随机噪音的波矢量。利用第一移动时窗，对多分量地震数据的波矢量进行计算出第一平均波矢量、并将第一平均波矢量进行中值滤波，以取得地滚波的真实模量，且将多分量地震数据的波矢量与地滚波的波矢量相减，以取得矢量时间序列。利用第二移动时窗，对所述矢量时间序列进行均值滤波，以取得第二平均波矢量，并对所述第二平均波矢量进行中值滤波以及在同一接收线内进行中值滤波，以压制矢量时间序列中的随机噪音的波矢量，进而取得净化后的有效信号的波矢量。
根据本发明的技术方案，通过利用第一移动时窗，对多分量地震数据的波矢量进行计算出第一平均波矢量、并将第一平均波矢量进行中值滤波，以取得地滚波的真实模量，且将多分量地震数据的波矢量与地滚波的波矢量相减，以取得矢量时间序列以及利用第二移动时窗，对所述矢量时间序列进行均值滤波，以取得第二平均波矢量，并对所述第二平均波矢量进行中值滤波以及在同一接收线内进行中值滤波，以压制矢量时间序列中的随机噪音的波矢量，进而取得净化后的有效信号的波矢量。如此一来，在不依赖于椭圆度和方向性，也不需要设置任何类型的数据协方差矩阵或复杂的道集的情况下，可以极大压制分量地震数据的波矢量中的噪声，而且不会破坏有效信号的高频和低频信息。
附图说明
此处所说明的附图用来提供对本发明的进一步理解，构成本申请的一部分，本发明的示意性实施例及其说明用于解释本发明，并不构成对本发明的不当限定。在附图中：
图1是根据本发明实施例的平均波矢量集的形成的示意图；
图2是根据本发明实施例的三维空间中平均波矢量的示意图；
图3是根据本发明实施例的第二移动时窗长度设置为地震子波旁瓣的近1/4的示意图；
图4是根据本发明实施例的从炮点三分量记录中的中值滤波到波矢量的示意图；
图5是根据本发明实施例的求解矩形的参数的迭代过程的示意图；
图6是根据本发明实施例的多分量地震数据的矢量去噪方法的流程图；
图7a是根据本发明实施例的主频为30hz的地震信号加上50％随机噪声和地滚波合成的两分量记录；
图7b是根据本发明实施例的地滚波的两分量记录；
图7c是根据本发明实施例的主频为30hz的地震信号加上50％随机噪声的两分量记录；
图8a是根据本发明实施例的分离的地滚波的波形图；
图8b是根据本发明实施例的夹带随机噪音的有效信号的波形图；
图9a是根据本发明实施例的原始两分量记录的矢量图；
图9b是根据本发明实施例的不采用任何缩放因子分离的地滚波矢量图；
图9c是根据本发明实施例的使用比例因子分离的地滚波矢量图；
图10a是根据本发明实施例的随机噪声衰减后的地震信号的波形图；
图10b是滤波后的随机噪声的波形图；
图11a是根据本发明实施例的分离的有效信号与随机噪声的矢量图；
图11b是根据本发明实施例的无缩放的随机噪声衰减后的有效信号矢量图；
图11c是根据本发明实施例的随机噪声衰减后有缩放的有效信号矢量图；
图11d是根据本发明实施例的滤波的随机噪声矢量图；
图12a是根据本发明实施例的Z分量的合成炮集记录；
图12b是根据本发明实施例的R分量的合成炮集记录；
图13a是根据本发明实施例的Z分量的合成炮集记录的频谱；
图13b是根据本发明实施例的R分量的合成炮集记录的频谱；
图14是根据本发明实施例的50ms移动时窗合成炮集记录中计算得到的波矢模量；
图15a是根据本发明实施例的Z分量的分离出的地面滚动能量；
图15b是根据本发明实施例的R分量的分离出的地面滚动能量；
图16a是根据本发明实施例的Z分量的地面滚动能量衰减后的合成炮集记录；
图16b是根据本发明实施例的R分量的地面滚动能量衰减后的合成炮集记录；
图17是根据本发明实施例的地面滚动能量衰减后合成炮集记录中计算得到的波矢模量；
图18a是根据本发明实施例的Z分量的滤波得到的最终炮集记录；
图18b是根据本发明实施例的R分量的滤波得到的最终炮集记录；
图19a是根据本发明实施例的Z分量的滤波得到的最终炮集记录的频谱；
图19b是根据本发明实施例的R分量的滤波得到的最终炮集记录的频谱；
图20a是根据本发明实施例的Z分量的分离出的随机噪音以及其它相关噪音的示意图；
图20b是根据本发明实施例的R分量的分离出的随机噪音以及其它相关噪音的示意图；
图21a是根据本发明实施例的从80ms移动时窗的现场炮集记录中计算得到的波矢模量；
图21b是根据本发明实施例的从地面滚动能量衰减后的现场炮集记录中计算得到的波矢模量；
图22a是根据本发明实施例的原始z分量的示意图；
图22b是根据本发明实施例的滤波后的z分量的示意图；
图22c是根据本发明实施例的z分量分离的噪声的示意图；
图23a是根据本发明实施例的原始z分量的频谱；
图23b是根据本发明实施例的滤波后的z分量的频谱；
图23c是根据本发明实施例的z分量分离的噪声的频谱；
图24a是根据本发明实施例的原始X分量的示意图；
图24b是根据本发明实施例的滤波后的X分量的示意图；
图24c是根据本发明实施例的X分量分离的噪声的示意图；
图25a是根据本发明实施例的原始X分量的频谱；
图25b是根据本发明实施例的滤波后的X分量的频谱；
图25c是根据本发明实施例的X分量分离的噪声的频谱；
图26a是根据本发明实施例的原始T分量的示意图；
图26b是根据本发明实施例的滤波后的T分量的示意图；
图26c是根据本发明实施例的T分量分离的噪声的示意图；
图27a是根据本发明实施例的原始T分量的频谱；
图27b是根据本发明实施例的滤波后的T分量的频谱；
图27c是根据本发明实施例的T分量分离的噪声的频谱。
具体实施方式
为使本发明的目的、技术方案和优点更加清楚，以下结合附图及具体实施例，对本发明作进一步地详细说明。
在三维三分量(3D3C)地震勘探中，测量的3C数据由地震检波器的一个垂直(Z-)分量和两个水平(R和T)分量同时记录，R分量与震源到检波器的方位角同向，T分量位于震源到检波器方向的正交方向。三分量振幅aZ、aR和aT在同一时间t形成信号矢量，如公式(1)所示：
A(t)＝aZ(t)eZ+aR(t)eR+aT(t)eT.， (1)
A(t)表示一个单独的波矢量，它不能反映真实的波矢量特征，因为存在不同的噪声，如地面滚动、随机噪声和任何其它类型的噪声，因此，A(t)是一个复合矢量，近似表示为公式(2)，如下所示：
A(t)＝B(t)+G(t)+N(t)。 (2)
其中，B(t)、G(t)和N(t)分别表示在t时刻，有效信号、地面滚动(地滚波)和随机噪音(连同其它噪音)的波矢量，且这些波矢量如公式(3)所示：
其中，eB、eG和eN表示单位波矢量，B(t)、G(t)和N(t)分别表示波矢量的模。
向3C记录中添加一个移动分析窗口，可以推出一个平均(或中值)矢量来表征数据窗口中的极化特征，并且可以为不同的用途和统计数据定制相应的移动窗口长度，例如，使用较大的移动窗口提供近零的随机噪声模值和相对准确的瞬时地面极化，使用较小的移动窗口，效果则正好相反。在本实施例中，将考虑在不同的移动窗口中使用均值(或中值)矢量来降低3C记录中的噪声。
地滚波一般在高振幅和低频椭圆极化的低速地层中靠近地表面传播。为了分离地滚波椭圆质点轨迹上的波矢量，在3C记录中添加了大量的移动时窗，可以从平均矢量集中导出中值矢量，以表征数据窗口中的极化。图1是根据本发明实施例的平均波矢量集的形成的示意图。在图1中，标号110表示地滚波的质点(圆点)运动轨迹。如图1所示，在第一移动时窗T1内，每个第一平均波矢量Kj可以由公式(4)计算，公式(4)如下所示：
其中，△t表示地震数据的时间采样率。这些第一平均波矢量形成1+T1/(2Δt)维矢量集，其中每个元素接近于地面样本在t时刻处的真实波矢量，然后，通过公式(5)计算一个矢量与所有其它矢量的距离的总和，公式(5)如下所示：
中值矢量由中值滤波法推导出，如公式(6)所示：
M(t)＝argmin(D(Kj))。 (6)
图2是根据本发明实施例的三维空间中平均波矢量的示意图。在图2中，标号210表示由公式(6)定义的中值矢量。如图2所示，M(t)位于均值矢量集的中心，与所有其它平均矢量的距离最小。
中值滤波器将导致矢量模量比真实滤波器稍小，这可以通过公式(7)给出的最小二乘法来恢复，公式(7)如下所示：
其中，mZ、mR和mT分别为振幅在Z、R和T分量上的投影，γ1(t)是用于修改中值矢量模量的第一比例因子，为了最小化目标函数Q(t)，可以利用公式(8)，如下所示：
然后，第一比例因子可由公式(9)求得，公式(9)如下所示：
其中，γ1(t)表示第一比例因子。因此，具有真实模量的地滚波矢量表示为公式(10)，如下所示：
G(t)＝γ1(t)M(t)。 (10)
将公式(10)与公式(2)相减(即公式(2)减公式(10))，可以得到由随机噪声和有效信号矢量组成的地滚分离波的波矢量时间序列C(t)，如公式(11)所示：
C(t)＝A(t)-G(t)＝B(t)+N(t)。 (11)
正确的移动时窗长度对于地滚波矢量特征的恢复至关重要，在计算地滚波的中值矢量过程中，有效信号被处理为噪声和随机噪声，只有当第一移动时窗长度大于波长的一半时，有效信号的较大模数波矢量才能通过中值滤波器进行滤波。因此，将第一移动时窗长度的下限设置为有效信号的最大视周期的一半，另外，第一移动时窗长度也应该小于地滚波的最小视周期的一半，因为第一移动时窗长度太大将导致平均波矢量集里出现较大误差。总的来说，中值波矢量的计算对移动时窗长度的上下限不敏感。
在某些情况下，由于吸收和衰减的影响，体波的带宽与地滚波的带宽会出现很多重叠，在地滚波的移动时窗应用时，这些低频地震波可产生大量假中值矢量。因此，利用极化椭圆率来区分椭圆的或线性的极化，并减少伪中值矢量，在整个移动时窗中，最大模数的波矢量可以被认为是半长轴，并且在与半长轴正交的方向上，容易找到具有最大模量的另一个波矢量作为半短轴，如果质点轨迹的半长轴和半短轴的比例大于阈值，将考虑体波的极化，并将当前矢量的模量设为最小值。
另外，与地滚波不同，为了压制随机噪声，使用较小的时窗来防止有效信号的主峰值被压制。如图3所示，第二移动时窗T2长度可以设置为地震子波旁瓣的近1/4，在第二移动时窗内，采用三段平滑方法来净化体波矢量，在第一阶段，使用均值滤波得到第二平均波矢量M1(t)，如公式(12)所示：
其中，T2表示第二移动时窗。均值滤波将所有波矢量在第二移动时窗内进行叠加，并取其中心平均矢量，三个分量矢量的叠加可以突出体波的矢量特征，压制噪声的矢量振幅。
图4是根据本发明实施例的从炮点三分量记录中的中值滤波到波矢量的示意图。在图4中，标号410表示地震采样点，标号420表示移动窗口，用于随时间进行的中值滤波，标号430表示移动窗口，用于在相同接收测线内的道集上进行的中值滤波。如图4所示，在第二阶段，对所有第二平均波矢量的波矢量进行中值滤波，以压制随机噪声矢量，如公式(13)所示：
其中，M2(t)表示压制随机噪声矢量后的波矢量。然后，在第三阶段，在同一接收测线内使用移动道窗(图4中的标号430)进行中值滤波，以压制非地面滚动相关的相干噪声，如公式(14)所示：
其中，M3(t)表示压制其它非地滚波的相干噪声的波矢量，k为中值滤波的输出道，K为中值滤波的移动道窗长度。如图5所示，由于相邻的地震道在同一接收测线上具有相似的波矢量特性，因此道集上的中值滤波应该在同一地震接收测线内进行。然后，可以通过公式(15)缩放中值矢量M3来得到净化后的有效信号矢量，公式(15)如下所示：
B(t)＝γ2(t)M3(t), (15)
并B且t，＝γ2(t)表M3示t第二比例因子且表示为公式(16)，如下所示：
其中，cZ、cR和cT分别表示滤除地滚波以后的波矢量在Z、R和T分量上的投影，m3Z、m3R和m3T分别表示表示压制其它非地滚波的相干噪声的波矢量在Z、R和T分量上的投影，cZeZ+cReR+cTeT＝C，m3ZeZ+m3ReR+m3TeT＝M3。这里的第二移动时窗T2可以设定为近似一个子波视周期。
以上，大略说明了本发明之实施例所需要运用到的相关公式，以下将提供对应的实施例来进行说明。根据本发明的实施例，提供了一种基于横波极化分析预测裂缝方位角的方法。
图6是根据本发明实施例的多分量地震数据的矢量去噪方法的流程图。
步骤S602，取得多分量地震数据的波矢量，其中多分量地震数据的波矢量包括有效信号、地滚波和随机噪音的波矢量。其中，所述多分量地震数据的波矢量如公式(2)所示。
步骤S604，利用第一移动时窗，对多分量地震数据的波矢量进行计算出第一平均波矢量、并将第一平均波矢量进行中值滤波，以取得地滚波的真实模量，且将多分量地震数据的波矢量与地滚波的波矢量相减，以取得矢量时间序列。在本实施例中，所述第一移动时窗的长度的下限为所述有效信号的最大视周期的一半，且所述第一移动时窗的长度小于所述地滚波的最小视周期的一半。
另外，第一平均波矢量如公式(4)所示。将第一平均波矢量进行中值滤波，以取得地滚波的真实模量可以由公式(5)、公式(6)、公式(7)、公式(8)、公式(9)和公式(10)进行。将多分量地震数据的波矢量与地滚波的波矢量相减，以取得矢量时间序列如公式(11)所示。
步骤S606，利用第二移动时窗，对所述矢量时间序列进行均值滤波，以取得第二平均波矢量，并对所述第二平均波矢量进行中值滤波以及在同一接收线内进行中值滤波，以压制矢量时间序列中的随机噪音的波矢量，进而取得净化后的有效信号的波矢量。其中，第二平均波矢量如公式(12)所示。对所述第二平均波矢量进行中值滤波以及在同一接收线内进行中值滤波可以由公式(13)和公式(14)进行。净化后的有效信号的波矢量可以由公式(15)和公式(16)产生。
图7a是根据本发明实施例的主频为30hz的地震信号加上50％随机噪声和地滚波合成的两分量记录。图7b是根据本发明实施例的地滚波的两分量记录。图7c是根据本发明实施例的主频为30hz的地震信号加上50％随机噪声的两分量记录。在图7a、图7b和图7c中，标号710表示Z分量，标号720表示R分量。如图7a、图7b和图7c所示，合成的1D2C测试记录里加入了50％的随机噪音和强烈的地面滚动干扰。明显可见，地震信号和地面滚动的视周期分别为40ms和100ms。
图8a是根据本发明实施例的分离的地滚波的波形图。图8b是根据本发明实施例的夹带随机噪音的有效信号的波形图。在图8a和图8b中，标号810表示Z分量，标号820表示R分量。如图8a和图8b所示，可以用20ms的移动时窗将地滚波(图8a)和带有随机噪声(图8b)的有效信号分开。在图8a和图8b中，开始20ms和结束20ms时的数据不能被滤除，因为它们位于移动时窗的边界位置。
图9a是根据本发明实施例的原始两分量记录的矢量图。图9b是根据本发明实施例的不采用任何缩放因子分离的地滚波矢量图。图9c是根据本发明实施例的使用比例因子分离的地滚波矢量图。在图9a、图9b和图9c中，标号910表示理论的地滚波偏振轨迹，标号920表示实际的地滚波偏振轨迹。由图9a、图9b和图9c中的极化波形图表明，尽管减小了波矢量的模量，但中值滤波仍可以提高地滚波的椭圆率。因此，在不采用任何比例因子(图9b)的情况下，分离的地滚波将呈现一个较小的椭圆偏振轨迹，但是，时变缩放因子可以使分离的地滚波的偏振轨迹(图9c)与理论情况更符合。
然后，进行了三个阶段的平滑处理并分离有效信号(图10a)和随机噪声(图10b)。在图10a和图10b中，标号1010表示Z分量，标号1020表示R分量。图11a是根据本发明实施例的分离的有效信号与随机噪声的矢量图。图11b是根据本发明实施例的无缩放的随机噪声衰减后的有效信号矢量图。图11c是根据本发明实施例的随机噪声衰减后有缩放的有效信号矢量图。图11d是根据本发明实施例的滤波的随机噪声矢量图。在图11a、图11b、图11c和图11d中，标号1110表示理论的地滚波偏振轨迹，标号1120表示实际的地滚波偏振轨迹。图11a、图11b、图11c和图11d表明，噪声衰减法可以很好地区分体波极化和随机噪声极化，使得体波极化更线性。将图11b与11c进行比较，明显可见，对来自三个阶段平滑的波矢量模量选取适当的缩放因子可以使本实施例的矢量去噪方法具有更好的保幅效果。
图12a是根据本发明实施例的Z分量的合成炮集记录。图12b是根据本发明实施例的R分量的合成炮集记录。利用图12a和图12b中所示的合成两分量炮集记录来说明噪音衰减方法的具体步骤，Z分量和R分量的有效信号主要是PP和PS波，它们是用30赫兹的RicKer子波合成的，加入的地滚波存在明显的频散现象，截止频率约为10Hz(图13a和图13b)，两个分量的随机噪声水平为200％，频带为0-500Hz。此外，还在本次测试的两个分量上加入了由地面覆盖的非地面滚动相关的相干噪声。
首先，选择50ms的移动时窗，在30Hz的Ricker子波半视周期与地滚波半视周期之间分离地滚波。图14是根据本发明实施例的50ms移动时窗合成炮集记录中计算得到的波矢模量。如图14所示，具有50ms移动时窗的合成炮集记录计算的波矢量模量主要反映地滚波能量，而有效信号和随机噪声都被大大地压制。图15a是根据本发明实施例的Z分量的分离出的地面滚动能量。图15b是根据本发明实施例的R分量的分离出的地面滚动能量。图16a是根据本发明实施例的Z分量的地面滚动能量衰减后的合成炮集记录。图16a是根据本发明实施例的R分量的地面滚动能量衰减后的合成炮集记录。从图15a、图15b和图16a、图16b中可见，在不损坏有效信号的前提下，地滚波能量从合成炮集记录中有效地分离了出来。此外，在地面滚动能量衰减之后，可以区分出线性非地面滚动相关的相干噪音。
其次，选择一个5ms的移动时窗，对地面滚动能量衰减后的合成炮集记录进行第一和第二阶段的平滑处理，在第三平滑阶段，选择一个3道的移动道窗。图17是根据本发明实施例的地面滚动能量衰减后合成炮集记录中计算得到的波矢模量如。图17所示，滤波参数有助于突出显示随机噪声背景下有效信号的波矢量模数。图18a是根据本发明实施例的Z分量的滤波得到的最终炮集记录。图18b是根据本发明实施例的R分量的滤波得到的最终炮集记录。图19a是根据本发明实施例的Z分量的滤波得到的最终炮集记录的频谱。图19b是根据本发明实施例的R分量的滤波得到的最终炮集记录的频谱。从最终滤波结果(图18a、图18b)和相应频谱(图19a、图19b)可以看出，本实施例的矢量去噪方法是很有效的，并且可以保持有效信号的高频和低频部分。图20a是根据本发明实施例的Z分量的分离出的随机噪音以及其它相关噪音。图20b是根据本发明实施例的R分量的分离出的随机噪音以及其它相关噪音。从图20a和图20b可以看出，随机噪声和其它相干噪声都被大大滤除，而且不会干扰到有效信号。
将本实施例的矢量去噪方法应用于中国西南四川盆地新场气田采集的3D3C数据的处理。图21a是根据本发明实施例的从80ms移动时窗的现场炮集记录中计算得到的波矢模量。图21b是根据本发明实施例的从地面滚动能量衰减后的现场炮集记录中计算得到的波矢模量。如图21a所示，80ms移动时窗可以很好地突出显示地面波的波矢量模量，但是，浅层中的初至波的波矢量模量也被突出显示。因此，只对图21a中标号2110所示的曲线下的区域中的数据进行地滚波能量衰减，为了衰减随机噪音和其它相干噪声，采用8ms移动时窗和一个3道的移动道窗。
本实施例的矢量去噪方法可以大大压制Z分量上的地滚波能量和随机噪声(如图22a、图22b和图22c所示)，滤除的噪声主要是宽频带部分(如图23a、图23b和图23c所示)。然而，Z分量上的初至波能量也被滤除了一些，这导致了Z分量频谱上发生地震主频的变化(如图23c所示)。在X分量(如图24a、图24b、图24c和图25a、图25b、图25c所示)上，除了地面滚动和随机噪声之外，还充分压制了诸如异常值和单频干扰的噪声。当然，也可以从T分量中找到满意的滤波结果(如图26a、图26b、图26c和图27a、图27b、图27c所示)。总的来说，本实施例的矢量去噪方法可以有效压制地面滚动能量和随机噪声，并保持有效信号的振幅。
综上所述，本发明通过利用第一移动时窗，对多分量地震数据的波矢量进行计算出第一平均波矢量、并将第一平均波矢量进行中值滤波，以取得地滚波的真实模量，且将多分量地震数据的波矢量与地滚波的波矢量相减，以取得矢量时间序列以及利用第二移动时窗，对所述矢量时间序列进行均值滤波，以取得第二平均波矢量，并对所述第二平均波矢量进行中值滤波以及在同一接收线内进行中值滤波，以压制矢量时间序列中的随机噪音的波矢量，进而取得净化后的有效信号的波矢量。如此一来，在不依赖于椭圆度和方向性，也不需要设置任何类型的数据协方差矩阵或复杂的道集的情况下，可以极大压制分量地震数据的波矢量中的噪声，而且不会破坏有效信号的高频和低频信息。
以上所述仅为本发明的实施例而已，并不用于限制本发明，对于本领域的技术人员来说，本发明可以有各种更改和变化。凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的权利要求范围之内。
Claims (1)
1.一种多分量地震数据的矢量去噪方法，其特征在于，包括：
取得多分量地震数据的波矢量，其中所述多分量地震数据的波矢量包括有效信号、地滚波和随机噪音的波矢量；
利用第一移动时窗，对多分量地震数据的波矢量进行计算出第一平均波矢量，并将所述第一平均波矢量进行中值滤波，以取得所述地滚波的真实模量，且将多分量地震数据的波矢量与所述地滚波的波矢量相减，以取得矢量时间序列，所述第一移动时窗的长度的下限为所述有效信号的最大视周期的一半，且所述第一移动时窗的长度小于所述地滚波的最小视周期的一半；
利用第二移动时窗，对所述矢量时间序列进行均值滤波，以取得第二平均波矢量，并对所述第二平均波矢量进行中值滤波以及在同一接收线内进行中值滤波，以压制矢量时间序列中的随机噪音的波矢量，进而取得净化后的有效信号的波矢量，所述第二移动时窗的长度设置为地震子波旁瓣的1/4；
所述多分量地震数据的波矢量满足如下公式：
A(t)＝B(t)+G(t)+N(t)，
其中，t表示时间，B(t)、G(t)和N(t)分别表示在t时刻的有效信号、地面滚动和随机噪音的波矢量；
所述第一平均波矢量满足如下公式：
其中，Kj表示第一平均波矢量，T1表示第一移动时窗，△t表示地震数据的时间采样率，
将所述第一平均波矢量进行中值滤波满足如下公式：
M(t)＝argmin(D(Kj))，
其中，M(t)表示中值矢量，D(Kj)表示所述第一平均波矢量中一个矢量与所有其它矢量的距离的总和且满足如下公式：
其中，所述地滚波的真实模量满足如下公式：
G(t)＝γ1(t)M(t)，
其中，γ1(t)表示第一比例因子；
所述第二平均波矢量满足如下公式：
其中，M1(t)表示第二平均波矢量，T2表示第二移动时窗，C(t)表示所述矢量时间序列，
对所述第二平均波矢量进行中值滤波以及在同一接收线内进行中值滤波满足如下公式：
其中，M2(t)表示压制随机噪声矢量后的波矢量，M3(t)表示压制其它非地滚波的相干噪声的波矢量，k为中值滤波的输出道，K为中值滤波的移动道窗长度，
所述净化后的有效信号的波矢量满足如下公式：
B(t)＝γ2(t)M3(t),
其中，γ2(t)表示第二比例因子。
Priority Applications (2)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201710867404.4A CN107966734B (zh) | 2017-09-22 | 2017-09-22 | 多分量地震数据的矢量去噪方法 |
| US15/944,713 US10761232B2 (en) | 2017-09-22 | 2018-04-03 | Vector denoising method for multicomponent seismic data |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201710867404.4A CN107966734B (zh) | 2017-09-22 | 2017-09-22 | 多分量地震数据的矢量去噪方法 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN107966734A CN107966734A (zh) | 2018-04-27 |
| CN107966734B true CN107966734B (zh) | 2019-04-02 |
Family
ID=61997455
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN201710867404.4A Active CN107966734B (zh) | 2017-09-22 | 2017-09-22 | 多分量地震数据的矢量去噪方法 |
Country Status (2)
| Country | Link |
|---|---|
| US (1) | US10761232B2 (zh) |
| CN (1) | CN107966734B (zh) |
Families Citing this family (5)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US11686872B2 (en) | 2019-12-23 | 2023-06-27 | Saudi Arabian Oil Company | Attenuation of guided waves using polarization filtering |
| CN111427090B (zh) * | 2020-04-29 | 2022-03-25 | 王仰华 | 稳健矢量中值滤波方法 |
| CN112925024B (zh) * | 2021-01-26 | 2023-10-20 | 中国石油化工股份有限公司 | 一种地震记录的方波压制方法 |
| US12000972B2 (en) | 2021-10-20 | 2024-06-04 | Saudi Arabian Oil Company | Attenuation of interface waves using single component seismic data |
| CN115932969B (zh) * | 2023-01-28 | 2023-06-16 | 中国地质大学（北京） | 波场分离方法及装置 |
Citations (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN101644782A (zh) * | 2009-08-25 | 2010-02-10 | 中国石化集团胜利石油管理局 | 基于极化滤波的多波多分量地震资料的去噪方法 |
| CN102692646A (zh) * | 2012-06-19 | 2012-09-26 | 北京多分量地震技术研究院 | 一种三维三分量矢量波场分离的方法和*** |
Family Cites Families (8)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US7042967B2 (en) * | 2003-03-03 | 2006-05-09 | Interdigital Technology Corporation | Reduced complexity sliding window based equalizer |
| US20060103892A1 (en) * | 2004-11-18 | 2006-05-18 | Schulze Mark A | System and method for a vector difference mean filter for noise suppression |
| US8050482B2 (en) * | 2006-09-28 | 2011-11-01 | Siemens Medical Solutions Usa, Inc. | System and method for online optimization of guidewire visibility in fluoroscopic systems |
| KR101092668B1 (ko) * | 2009-06-17 | 2011-12-13 | 서울대학교산학협력단 | 파형 역산을 이용한 지하 구조의 영상화 장치와 방법 |
| CN101893720B (zh) * | 2010-07-02 | 2012-09-05 | 中国科学院地质与地球物理研究所 | 一种地震波的矢量波场分离与合成的方法和*** |
| CN101937100B (zh) * | 2010-08-17 | 2012-10-03 | 中国科学院地质与地球物理研究所 | 一种叠前深度偏移方法 |
| EP2781937B1 (en) * | 2013-03-22 | 2022-09-28 | CGG Services SAS | Vector-dip filtering of seismic data in the time-frequency domain |
| CN104597499B (zh) * | 2013-10-31 | 2017-02-08 | 中国石油天然气集团公司 | 可控震源独立同步激发地震数据邻炮干扰压制方法和装置 |
-
2017
- 2017-09-22 CN CN201710867404.4A patent/CN107966734B/zh active Active
-
2018
- 2018-04-03 US US15/944,713 patent/US10761232B2/en active Active
Patent Citations (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN101644782A (zh) * | 2009-08-25 | 2010-02-10 | 中国石化集团胜利石油管理局 | 基于极化滤波的多波多分量地震资料的去噪方法 |
| CN102692646A (zh) * | 2012-06-19 | 2012-09-26 | 北京多分量地震技术研究院 | 一种三维三分量矢量波场分离的方法和*** |
Non-Patent Citations (1)
| Title |
|---|
| M. Landrø.Attenuation of seismic water-column noise,tested on seismic data from the Grane field.《GEOPHYSICS》.2007, |
Also Published As
| Publication number | Publication date |
|---|---|
| US20190094400A1 (en) | 2019-03-28 |
| CN107966734A (zh) | 2018-04-27 |
| US10761232B2 (en) | 2020-09-01 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN107966734B (zh) | 多分量地震数据的矢量去噪方法 | |
| Tian et al. | Variable-eccentricity hyperbolic-trace TFPF for seismic random noise attenuation | |
| CN101598812B (zh) | 去除数字检波器单点接收地震记录中的异常噪声方法 | |
| CN105911585B (zh) | 一种地震记录规则干扰波的提取方法及装置 | |
| CN109164492B (zh) | 一种提取套管井地层声波速度的方法 | |
| Hu et al. | Ground-roll noise extraction and suppression using high-resolution linear Radon transform | |
| Margrave et al. | The Hussar low-frequency experiment | |
| CN105116443A (zh) | 一种低频信号的能量补偿方法及装置 | |
| EP2245484A2 (en) | Method for detecting and/or processing seismic signals | |
| CN108957551B (zh) | 基于重构地面力信号的可控震源谐波压制方法 | |
| CN106324702B (zh) | 一种地震干涉法成像观测***设计的定量评价方法 | |
| CN105093282A (zh) | 基于频率约束的能量置换面波压制方法 | |
| Chen et al. | Robust adaptive polarization analysis method for eliminating ground roll in 3C land seismics | |
| CN109975873A (zh) | 一种逆时偏移成像去除低频噪音的方法及*** | |
| CN115826039A (zh) | 一种时间切片分类模型训练方法、***及应用方法、*** | |
| CN104502968B (zh) | 基于阈值多级中值滤波的可控震源地震数据检测方法 | |
| CN107402406B (zh) | 一种压制地震数据中大钻噪声的方法 | |
| Douglas | Making the most of the recordings from short-period seismometer arrays | |
| CN106707333A (zh) | 数字检波器的叠前去噪方法及装置 | |
| Du et al. | Research and application of Rayleigh wave extraction method based on microtremors signal analysis | |
| CN105425297A (zh) | 一种压制虚反射信号的方法及装置 | |
| CN110764147A (zh) | 一种基于vmd局部f-x谱分解的沙漠勘探弱信号恢复方法 | |
| CN110231651A (zh) | 一种基于振动信号矢量分解原理的叠前地震数据去噪方法 | |
| Lévy et al. | Rayleigh waves in seismic signals of rockfalls | |
| Tan et al. | Combined adaptive multiple subtraction based on event tracing and Wiener filtering |
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 |