matlab实现鬼波信号压制算法（附鬼波算法压制工具包）代码实践--第一篇频率-空间域自适应鬼波压制

2024-05-30 03:31:34

matlab实现鬼波信号压制算法（附鬼波算法压制工具包）代码实践

涵盖了频率-空间域、频率-波数域、拉东域鬼波压制算法建议实践之前熟练掌握各个域鬼波压制方法的原理，才能对代码有更深入的了解。

本次介绍共分三个篇章，

第一篇介绍基于峭度因子最大化的频率空间域鬼波压制算法

第二篇介绍频率波数域鬼波压制算法

第三篇介绍拉东域鬼波压制算法

具体实现思路代码详见每篇最后附上的代码，另需下载配置鬼波压制算法工具包（https://download.csdn.net/download/zw461450141/16811437?spm=1001.2014.3001.5501），工具包里面内容如下：

配置好工具包后，直接运行文中附上的代码，也可结合自己模拟的数据，修改代码中对应的参数，参数在代码中均有注释

鬼波信号压制算法代码实践

第一篇频率-空间域自适应鬼波压制方法

1.1 频率-空间域自适应鬼波压制思路

1.2 频率-空间域自适应鬼波压制效果对比

1.2.1 原始数据切除直达波对比图

1.2.2 鬼波压制前后单道数据对比图

1.2.3 鬼波压制前后剖面数据对比图

1.2.4 鬼波压制前后数据F-K谱对比图

1.3 代码实现过程

第一篇频率-空间域自适应鬼波压制方法

这里实现的算法是基于峭度最大化的频率-空间域自适应鬼波压制，压制原理可参考相关文献：Deghosting for variable depth streamer based on super-Gaussianity

该方法可用于水平拖缆、斜缆数据的鬼波压制

1.1 频率-空间域自适应鬼波压制思路

1.2 频率-空间域自适应鬼波压制效果对比

1.2.1 原始数据切除直达波对比图

1.2.2 鬼波压制前后单道数据对比图

1.2.3 鬼波压制前后剖面数据对比图

1.2.4 鬼波压制前后数据F-K谱对比图

1.3 代码实现过程

实现过程需下载工具函数包：

结合自己模拟的拖缆数据，进行鬼波压制效果对比，部分代码如下：

%-------------------------------------------------------
% 频率--空间域鬼波压制算法实践
% 结合模拟的海上斜缆采集数据 对算法进行验证
%---------------------------------------------------------------------
%% 加载数据
clc;clear;close all;
addpath('./syn_data');
addpath(genpath('./utility'));
% addpath(genpath('./crewes'));
d = load('rcv_18_101.dat');
%% 直达波切除
d_rawc= derectcut(d,1,200,101,400);
figure; subplot(1,2,1); imagesc(d); colormap('redblue'); title('原始数据'); caxis([-0.05,0.05]);
subplot(1,2,2); imagesc(d_rawc); colormap('redblue'); title('切除直达波数据'); caxis([-0.05,0.05]);
%% 给数据加噪声
d_rawc = addrandnoise(d_rawc,9,10);
%% 数据采集参数设置
dt = 0.0005;  % 采样时间
dx = 4 ;      % 道间距
[nt,ntra] = size(d_rawc);  % nt采样总时长  ntra道数
rcvdepth = 18*ones(1,ntra);  % 检波器深度   水平拖缆  未考虑起伏海面的影响
xx = 0:dx:(ntra-1)*dx;
tt = 0:dt:(nt-1)*dt;
ks = 0.5; % 正则化系数
%% 基于峭度最大化的频率空间域自适应鬼波压制
[d_raw_deghost] = deghost_fx(d_rawc,dt,rcvdepth,ks);
%-------------压制前后  单道地震数据对比
nn = 80;
figure;plot(d_rawc(:,nn),'r');hold on; plot(d_raw_deghost(:,nn),'b');legend('压制前','压制后');
%-------------压制前后 地震剖面对比-----------------
figure; subplot(1,2,1);imagesc(xx,tt,d_rawc);
colormap('jet');caxis([-0.05,0.05]);
title('原始数据');ylabel('时间{\itt} / s');xlabel('距离{\itd} / m');
subplot(1,2,2);imagesc(xx,tt,d_raw_deghost);
colormap('jet');caxis([-0.05,0.05]);title('鬼波压制后的数据'); ylabel('时间{\itt} / s');xlabel('距离{\itd} / m');
%--------------压制前后 地震数据fk谱对比------------------
ff = [0 120];       % 设置频率显示范围
kk = [-0.05,0.05];  % 设置波数显示范围
figure,subplot(1,2,1); imag_fk_spectrum(d_rawc,dt,dx,ff,kk);
title('原始数据');ylabel('频率{\itf} / Hz');xlabel('波数{\itk} / m');
caxis([0,0.5]);
subplot(1,2,2); imag_fk_spectrum(d_raw_deghost,dt,dx,ff,kk);
ylabel('频率{\itf} / Hz');xlabel('波数{\itk} / m');
title('鬼波压制后的数据');caxis([0,0.4]);function dnew = addrandnoise(d,L,snr)
% L 噪声强度
% snr 信噪比
if nargin == 1L = 5;snr = 3;
end
dmax  = max(max(d));
op = hamming(L);
ops = sum(sum(op));
op = op/ops;
Noise = conv2(randn(size(d)),op,'same');
Noisemax  = max(max(Noise));
dnew = d + (dmax/Noisemax)*Noise/snr;function [u_true] = deghost_fx(d,dt,z,ks)
% d   二维地震数据[t，x]
% dt  采样时间
% z   检波器深度   length(z)= size(d,2)
% ks  正则化系数
v_wat = 1500;  % 海水速度
[nt,ntra]=size(d); % nt 时间  ntra 道数
% 将地震数据进行fft变换，从时间域变换到频率域
NFFT = 2^nextpow2(nt); % Next power of 2 from length of N
tmax = (NFFT-1)*dt;
t = 0:dt:tmax;
ntpad = NFFT;
d = [d;zeros(ntpad-nt,ntra)];
[specfx,f] = fftrl(d,t,0,ntpad);  % fft变换  t-x域变换到f-x域
[rp,~] = size(specfx);
%% 对每一道进行搜索 自适应求取最佳鬼波延迟时间  延迟时间和沉放深度有关
%成像条件
R = -1;  % 海表面反射系数
z = ones(rp,1)*z;     %理论沉放深度
kz = f/v_wat;
NZ  = 100;
IR = zeros(1,NZ);
dh_new = zeros(rp,ntra);
dh_true = zeros(rp,ntra);
for j=1:ntratem_dz=linspace(-z(1,j)/3,2,NZ);    %控制深度z的范围  在范围内搜索最佳值for i=1:NZdh_new(:,j)=z(:,j)+tem_dz(i);G1 = 1+R*exp(-2*pi*1i*2*dh_new(:,j).*kz); %计算鬼波算子G12 = abs(G1).^2+ks; u_fx = (specfx(:,j).*conj(G1))./G12;[tem_u,~]=ifftrl(u_fx,f); %------峭度因子-------------ku = kurtosis(tem_u);    % 求鬼波压制后每道的峭度值  yes max%--------峰值因子----------------
%         ma = max(tem_u); %最大值
%         mi = min(tem_u); %最小值
%         pk = ma-mi;  %峰-峰值
%         rm = rms(tem_u); %均方根   yes  min
%         C = pk/rm;   %峰值因子  yes max%------------------------IR(i) = ku;   end[~,index]=max(IR);dh_true(:,j) = z(:,j) + tem_dz(index);
end
%% 采用搜索的自适应鬼波延迟时间  进行鬼波压制
ks1 = 0.05;
kz = f*ones(1,ntra)/v_wat;
G1 = 1+R*exp(-2*pi*1i*2*dh_true.*kz);
G12 = abs(G1).^2+ ks1;
u_fx_true = (specfx.*conj(G1))./G12;
[u_true,~] = ifftrl(u_fx_true,f);

matlab实现鬼波信号压制算法（附鬼波算法压制工具包）代码实践--第一篇频率-空间域自适应鬼波压制相关推荐

matlab实现鬼波信号压制算法代码实践--第二篇频率-波数域鬼波压制
第二篇 matlab实现频率-波数域鬼波压制方法本篇用来介绍频率-波数域鬼波压制的实现思路和压制效果算法实现思路见2.3节,除了文中代码外,需配置鬼波压制算法工具包(https://downlo ...
python 处理锯齿波信号
python 处理锯齿波信号预期学习目标(ILO): 您应该了解傅里叶分析在时域中周期性信号的基础,即如何将周期性时域信号分解为基频. 能够从头编程周期波形,并将其与"scipy&quo ...
自适应陷波滤波器（Notch滤波器）
自适应滤波器以其自身调节参数的能力以及对信号和噪声的先验知识要求少的优势得到了广泛的应用. 自适应陷波滤波器相比于普通的数字滤波器具有较窄的阻带和通带,带外衰减较快.具有滤波输出(Y输出)和陷波输出( ...
推荐算法(四)——经典模型 DeepFM 模型详解及代码实践
目录 1 介绍 2 模型结构 3 实验结果 4 总结 5 代码实践 1 介绍 DeepFM 是华为诺亚方舟实验室在 2017 年提出的模型. 论文传送门: A Factorization-Machin ...
matlab实现鬼波信号压制算法代码实践--第三篇拉东域鬼波压制
第三篇 matlab实现拉东域鬼波压制方法本篇用来介绍拉东域鬼波压制的实现思路和压制效果该方法可用于水平拖缆.斜缆数据的鬼波压制算法实现思路见3.3节,除了文中代码外,需配置鬼波压制算法工具包 ...
【信号去噪】基于蚁群算法优化小波阈值实现信号去噪附matlab代码
1 简介基于硬阈函数和软阈函数的小波去噪算法处理的信号分别存在着偏差和方差过大的缺点,为有效解决这一问题,提出基于蚁群算法优化小波变换去噪算法.并采用常用的信号用matlab对去噪效果进行了仿真.仿 ...
matlab decomposition filters,MATLAB小波去噪求助（附算法和显示图片）！不知自己哪个地方出了问题，求指点！ - 信息科学 - 小木虫 - 学术科研互动社区...
star2883 不是非常明白楼主的意思.我的理解是: 问题1:为什么中间两幅图,显示的都是空白? 答:你是指Figure2和Figure3的显示中,全是白色吗?这个是因为你用了colormap(ma ...
matlab复杂噪声产生实验报告,基于小波信号的噪声消除matlab实验报告.docx
基于小波信号的噪声消除matlab实验报告.docx 南京师范大学物理科学与技术学院医用电子学论文论文名称:基于小波变换的心电信号噪声消除院系:物科院专业:电路与系统姓名:聂梦雅学号:12100204 ...
基于ITD实现的轴承故障信号分解并附Matlab源码
基于ITD实现的轴承故障信号分解并附Matlab源码轴承是旋转机械中重要的支撑部件,其状态的健康程度对设备运行的性能和寿命有着至关重要的影响.因此,轴承状态监测和故障诊断成为了研究热点.本文基于IT ...

最新文章

热门文章