基于Tau-p变换的矢量中值滤波方法
【技术领域】
[0001] 本发明是一种基于Tau-p变换的矢量中值滤波方法，实现对地震数据的高保真滤 波去噪。
【背景技术】
[0002] 在野外地震数据中，为了消除直达波和折射波对浅层的反射波造成的干扰。常用 处理方法为直接切除，但直接切除的同时，很可能会切除远偏移距的有效波。目前，人们越 来越重视远偏移距的有效波信息，直接切除已经无法达到人们的期望。为了在消除直达波 和折射波干扰的同时，能够保护远偏移距信息不被破坏甚至完全切除，使用滤波法达到去 噪目的的次数越来越多。
[0003] 矢量中值滤波是一种由中值滤波发展而来的基于排序的多参数滤波，由于其在消 除噪声，保护边界有良好的稳健性，而且实现结构相对简单，因此在地震数据处理中引起了 一定的关注。但是由于矢量中值滤波没有考虑数据的速度信息，在对含有低速度波的地震 数据进行处理时，由于各道检波器接收到的地震波到达时间差距较大，造成了有效波衰减 的问题。
[0004] 本算法的目的在于矢量中值滤波之前，使用Tau-p变换对地震数据进行波场分 离，这样既能分别得到不同波的地震数据，又能确定对不同波的瞬时斜率，进而确定矢量滤 波窗口的方向，使得矢量滤波窗口内地震波的相关性最强，提高矢量中值滤波对地震数据 的滤波处理效果。
【发明内容】
[0005] 本方法可以解决传统矢量中值滤波算法中固定窗口移动方向所导致的速度较低 的有效波易被衰减的问题，能够在去除噪声的同时保留有效波，实现对地震数据的高保真 滤波去噪。
[0006] 基于Tau-p变换的矢量中值滤波方法的基本实现步骤如下：
[0007] Tau-p变换部分包括：
[0008] 步骤一：将时间-空间域域中的数据d经过FFT变换到频率-空间域，得到数据 D；
[0009] 步骤二：将频率-空间域中的数据D通过Radon变换到频率-参数域，得到数据M；
[0010] 步骤三：将频率-参数域中的数据M通过IFFT变换到Tau-p域，得到数据DD;
[0011] 步骤四：分析DD在Tau-p域中的分布形式，得到不同波的瞬时斜率，并在Tau-p域 中加窗口，将不同波分离开来；
[0012] 步骤五：通过Tau-p反变换将不同波反变换到时间-空间域，得到不同波分离后对 应的地震数据其中j为不同波对应的地震数据的编号，并记录波类型的总数m;
[0013] 矢量中值滤波部分包括：
[0014] 步骤六：从分离的波场信息中，选择第j个地震数据子集&，根据其瞬时斜率，设 计矢量滤波窗口，其中j对应于步骤五中地震数据子集的编号，其初始值为1;
[0015] 步骤七：根据步骤六设计的矢量滤波窗口从Xj中取出数据，并求出窗口内中地震 数据的平均值；
[0016] 步骤八：计算窗口中每个数据到平均值的距离，距离最小的数据为该窗口中数据 的中值，并代替窗口中原中心值输出，移动窗口，直至处理完Xj中所有数据；
[0017] 步骤九：对j进行加1，重复步骤六至步骤八，直至j为波类型的总数，即对所有波 都完成滤波，之后将所有滤波后的地震数据子集进行叠加，得到最终滤波处理后的地震数 据。
[0018] 进一步地，所述步骤一中，FFT变换即快速傅里叶变换的具体公式为： D=ffd 其中d为时间-空间域的地震数据组成的数据矩阵，列向量表示每道地震数据的时域 采样序列，行向量表示每个采样时刻不同道的采样序列，矩阵维数为NXQ，N为采样点数，Q 为地震道数；D为FFT变换后的频率-空间域数据矩阵，列向量表示每道地震数据的频域序 列，行向量表示每个频率点上不同道的采样序列，矩阵维数为NXQ;W为傅里叶变换矩阵， 维数为NXN该矩阵中的各元素为
【主权项】
1. 一种基于Tau-p变换的矢量中值滤波方法，用于地震数据处理，实现对地震数据的 高保真滤波去噪，其特征在于，包括Tau-p变换部分和矢量中值滤波两个部分，其中所述 Tau-p变换部分包括； 步骤一；将时间-空间域域中的数据d经过FFT变换到频率-空间域，得到数据D;步骤二：将频率-空间域中的数据D通过Radon变换到频率-参数域，得到数据M; 步骤S;将频率-参数域中的数据M通过IFFT变换到Tau-p域，得到数据孤； 步骤四：分析DD在Tau-p域中的分布形式，得到不同波的瞬时斜率，并在Tau-p域中加 窗口，将不同波分离开来； 步骤五：通过Tau-p反变换将不同波反变换到时间-空间域，得到不同波分离后对应的 地震数据子集Xj.，其中j为不同波对应的地震数据子集的编号，并记录波类型的总数m; 所述矢量中值滤波部分包括： 步骤六；从分离的波场信息中，选择第j个地震数据子集Xj.，根据其瞬时斜率，设计矢 量滤波窗口，其中j对应于步骤五中地震数据子集的编号，其初始值为1 ; 步骤走；根据步骤六设计的矢量滤波窗口从Xj.中取出数据，并求出窗口内中地震数据 的平均值； 步骤八：计算窗口中每个数据到平均值的距离，距离最小的数据为该窗口中数据的中 值，并代替窗口中原中屯、值输出，移动窗口，直至处理完Xj.中所有数据； 步骤九：对j进行加1，重复步骤六至步骤八，直至j等于步骤五中得到的总数m，即对 所有波都完成滤波，之后将所有滤波后的时距曲线进行叠加，得到最终滤波处理后的地震 数据。
2. 如权利要求1所述的方法，其特征在于，所述步骤一中，FFT变换即快速傅里叶变换 的具体公式为： D=Wd 其中d为时间-空间域的地震数据组成的数据矩阵，列向量表示每道地震数据的时域 采样序列，行向量表示每个采样时刻不同道的采样序列，矩阵维数为NXQ，N为采样点数，Q 为地震道数；D为FFT变换后的频率-空间域数据矩阵，列向量表示每道地震数据的频域序 列，行向量表示每个频率点上不同道的采样序列，矩阵维数为NXQ;W为傅里叶变换矩阵， 维数为NXN该矩阵中的各元素为二e-Zmik/A/，其中n为正整数，表示时域点数且n =1, 2,--N，k为正整数，表示频域点数且k= 1，2,--N。
3. 如权利要求1所述的方法，其特征在于，所述步骤二中，使用阻巧最小二乘法，将频 率-空间域中的数据矩阵D通过Radon变换到频率-参数域得到数据矩阵M，具体形式为：
其中阻巧因子，I为PXP的对角矩阵，〇T是D的转置，4是Lk的转置，L巧第 k个频率点对应的转换系数矩阵、其维数为PXQ、矩阵中的元素为fk为f-x域 第k个频率点对应的频率值，k= 1，2,…，N，Pi为第i个射线参数，i= 1，2,…，P，P为射 线参数的总个数，X。为f-x域中第q道的空间偏移量，q= 1，2,…，Q。
4. 如权利要求1所述的方法，其特征在于，所述步骤S中，所述的IFFT变换为FFT变换 的反变换，其具体公式为：
其中妒是W的共辆转置，IFFT的点数与FFT-致，反变换后得到的数据矩阵DD即为地 震数据的Tau-p域表示。
5. 如权利要求1所述的方法，其特征在于，所述步骤四中，Tau-p域中横坐标的几何意 义为时间-空间域中地震波前表现出的时间-距离曲线上对应点的瞬时斜率，纵坐标的几 何意义为该曲线对应的截距；如果时间-空间域中地震波波前表现出的时间-距离曲线为 直线，变换后在化U-P域中表现为一个点，直线的斜率不同，点对应的横坐标位置也不同； 如果时间-空间域中地震波前表现出的时间-距离曲线为双曲线，变换后在Tau-p域中表 现为一个楠圆；由于波前在时间-空间域中的时间-距离曲线不同，表现出来的瞬时斜率也 不同，在Tau-p域中加窗口后，能将各种波分离开。
6. 如权利要求1所述的方法，其特征在于，所述步骤五中，利用化U-P域中所加的窗口， 将各种波所在窗口内的数据进行Tau-p反变换到时间-空间域，其反变换的具体步骤为， FFT变换一Radon反变换一IFFT变换，经过反变换后，得到不同类型波分离后各类波的地震 数据子集Xj.，其中j为不同类型波对应的地震数据的编号，并记录波类型的总数m。
7. 如权利要求1所述的方法，其特征在于，所述步骤六中，地震数据子集Xj.的矢量滤波 窗口的方向性由该数据在Tau-p域中的瞬时斜率决定。
8. 如权利要求1所述的方法，其特征在于，所述步骤走中，在矢量滤波窗口内对数据沿 自左向右、自上而下的顺序取取数，并求出窗口内的地震数据子集Xj.的平均值。
9. 如权利要求1所述的方法，其特征在于，所述步骤八中，在滤波窗口中求每个数据到 平均值的距离，距离最小的数据值为该滤波窗口内地震数据中值，代替该滤波窗口内的原 中屯、值输出，然后将该滤波窗口按斜率方向移动一个矢量距离，直至完成对Xj.内所有数据 的滤波处理。
10. 如权利要求1所述的方法，其特征在于，所述步骤九中，对j加1，重复步骤六至步 骤八，直至j等于步骤五中的总数m，将所有滤波后的地震数据子集进行叠加，得到最终滤 波处理后的地震数据。
【专利摘要】本发明公开一种基于Tau-p变换的矢量中值滤波方法。本算法包括两个部分：第一部分是通过Tau-p变换，将地震数据从时间-空间域变换到Tau-p域，在Tau-p域中，得到不同波的瞬时斜率，并通过加窗的方法，将斜率不同的各种波分离开来，分离后的数据通过Tau-p反变换回时间-空间域，分别得到其对应的时间-空间域数据子集；第二部分是利用第一部分所得的瞬时斜率，确定滤波时窗的移动方向，对相应的时间-空间域数据子集进行矢量中值滤波，最后将滤波后的数据子集叠加，得到处理后的地震数据。本方法可以解决传统矢量中值滤波算法中固定窗口移动方向所导致的速度较低的有效波易被衰减的问题，能够在去除噪声的同时保留有效波，实现对地震数据的高保真滤波去噪。
【IPC分类】G01V1-36
【公开号】CN104849760
【申请号】CN201510294266
【发明人】郑晶, 杜文凤, 朱逸馨, 姜天琪
【申请人】中国矿业大学(北京)
【公开日】2015年8月19日
【申请日】2015年6月2日