MATLAB加随机噪声并去噪

MATLAB加随机噪声并去噪

 

1. 从某数据集中获取实际地震数据。时间采样间隔dt假定为1ms。添加指定信噪比水平的噪声，噪声能量/信号能量(1/SNR) = 0.3 + 10/100；
2. 对添加噪声后的数据做频域空间域去噪处理；
3. 画出查验去噪结果图。

clc,clear;%清除工作区以及命令行内容

load SeismicData3D.mat  %加载实际地震数据集

a=SeismicData3D(:,:,11); %从数据集中选择地震测线

b=dtInSeconds;%获取采样时间间隔

[x,y]=size(a);%获取数据矩阵的行和列大小,x为行大小，y为列大小

t(:,1)=((1:x)-1).*b;%计算总共采样时间

z(:,1)=1:1:y;

%---添加噪声---%

NP=0.3+11/100;%设定信噪比

Noise=zeros(x,y);%预分配内存，加快计算

for i=1:y%从1到y循环y次

    Noise(:,i)=rand(x,1)-0.5;%产生-0.5-0.5之间的随机噪声

    NC=sqrt(NP*sum(a(:,i).^2)/sum(Noise(:,i).^2));%计算随机噪声的缩放系数

    Noise(:,i)=NC*Noise(:,i);%产生最终的随机噪声

    a1(:,i)=a(:,i)+Noise(:,i);%对实际地震数据添加噪声

end

%---去除噪声---%

%频率空间域去噪声的参数设定

LenghtOfOperator = 15;%0perator length算子长度

FreqLow = 1; % Min frequency (Hz)to process最小频率

FreqHigh = 120;% Max frequency (Hz)to process最大频率

mu = 0.01;% trade-off parameter去噪声时通常都会涉及到的折中调节算子

DataFiltered = FXDenoising(a1, b, LenghtOfOperator, mu, FreqLow,FreqHigh);%调用频率空间域滤波

NoisePredicted = a1-DataFiltered;%计算出预测的噪声数据

save data%保存数据集到文件夹

%---画图---%

iDataCaxis = [min(a(:)) max(a(:))]*0.5;%设定画图时的数据范围

iFontSize = 12;%设定字体大小

iFontName = 'Helvetica';%设定字体

iFigureNameXYZ = [3 0.05 0];%设定图题标注在图中的位置

iFigureNameColor = 'yellow';%设定图题标注的字体颜色

figure(1);clf;%注意clf：Clear current figure清空当前图片

set(gcf,'color','white','Units','inches','Position',[0.3 1 12 6]);

% 'color','white' 指定图片的背景色为白色，替代默认的灰色，方便截图

% 'Units','inches' 设定窗口大小位置时所采用的单位，英尺

% 'Position',[0.3 1 12 6] 设定图片的窗口位置和大小

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%画输入的地震数据%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

subplot(1,4,[1]);

imagesc(z,t,a);%画输入的地震数据

text(iFigureNameXYZ(1), iFigureNameXYZ(2), iFigureNameXYZ(3),'Input data',...

    'color',iFigureNameColor,'rotation',0,'FontName',iFontName,'FontSize',...

    iFontSize,'FontWeight','normal');

%用text在任意位置输入图题，设定字体颜色、字体、是否加粗、旋转

caxis(iDataCaxis);%设定画图时数据的幅值范围

colormap(gray);%设定colormap

xlim([min(z) max(z)]);%设定x轴的坐标范围

ylim([min(t) max(t)]);%设定y轴的坐标范围

set(gca,'xTick',[0:50:max(z)]);%设定x轴的坐标刻度

set(gca,'yTick',[0:0.1:max(t)]);%设定y轴的坐标刻度

colorbar('Ticks',[round(min(iDataCaxis)):4:round(max(iDataCaxis))],...

    'location','southoutside');%设定colorbar的显示刻度及位置

% 设定图题

% title('Input data','FontName',iFontName,'FontSize',iFontSize,'FontWeight','bold');

% 设定x轴的标签

xlabel('Trace number','FontName',iFontName,'FontSize',iFontSize,'FontWeight','normal');

% 设定y轴的标签

ylabel('Time (s)','FontName',iFontName,'FontSize',iFontSize,'FontWeight','normal');

set(gca,'Box','on','XAxisLocation','top','ydir','reverse','FontName',iFontName,...

    'FontSize',iFontSize,'FontWeight','normal');

% 'XAxisLocation','top'设定x轴在顶部

% 'ydir','reverse'设定y轴逆序

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%画去噪后的地震数据%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

subplot(1,4,[2]);

imagesc(z,t,DataFiltered);%画去噪后的地震数据

text(iFigureNameXYZ(1), iFigureNameXYZ(2), iFigureNameXYZ(3),...

    'Data after FX decon','color',iFigureNameColor,'rotation',0,'FontName',...

    iFontName,'FontSize',iFontSize,'FontWeight','normal');

caxis(iDataCaxis);%设定画图时数据的幅值范围

colormap(gray);%设定colormap

xlim([min(z) max(z)]);%设定x轴的坐标范围

ylim([min(t) max(t)]);%设定y轴的坐标范围

set(gca,'xTick',[0:50:max(z)]);%设定x轴的坐标刻度

set(gca,'yTick',[0:0.1:max(t)]);%设定y轴的坐标刻度

colorbar('Ticks',[round(min(iDataCaxis)):4:round(max(iDataCaxis))],...

    'location','southoutside');%设定colorbar的显示刻度及位置

% 设定图题

% title('Input data','FontName',iFontName,'FontSize',iFontSize,'FontWeight','bold');

% 设定x轴的标签

xlabel('Trace number','FontName',iFontName,'FontSize',iFontSize,'FontWeight','normal');

% 设定y轴的标签

ylabel('Time (s)','FontName',iFontName,'FontSize',iFontSize,'FontWeight','normal');

set(gca,'Box','on','XAxisLocation','top','ydir','reverse','FontName',iFontName,...

    'FontSize',iFontSize,'FontWeight','normal');

% 'XAxisLocation','top'设定x轴在顶部

% 'ydir','reverse'设定y轴逆序

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%画预测的噪声%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

subplot(1,4,[3]);

imagesc(z,t,NoisePredicted);%画预测的噪声

text(iFigureNameXYZ(1), iFigureNameXYZ(2), iFigureNameXYZ(3),'Estimated noise',...

    'color',iFigureNameColor,'rotation',0,'FontName',iFontName,'FontSize',...

    iFontSize,'FontWeight','normal');

caxis(iDataCaxis);%设定画图时数据的幅值范围

colormap(gray);%设定colormap

xlim([min(z) max(z)]);%设定x轴的坐标范围

ylim([min(t) max(t)]);%设定y轴的坐标范围

set(gca,'xTick',[0:50:max(z)]);%设定x轴的坐标刻度

set(gca,'yTick',[0:0.1:max(t)]);%设定y轴的坐标刻度

colorbar('Ticks',[round(min(iDataCaxis)):4:round(max(iDataCaxis))],...

    'location','southoutside');%设定colorbar的显示刻度及位置

% 设定图题

% title('Input data','FontName',iFontName,'FontSize',iFontSize,'FontWeight','bold');

% 设定x轴的标签

xlabel('Trace number','FontName',iFontName,'FontSize',iFontSize,'FontWeight','normal');

% 设定y轴的标签

ylabel('Time (s)','FontName',iFontName,'FontSize',iFontSize,'FontWeight','normal');

set(gca,'Box','on','XAxisLocation','top','ydir','reverse','FontName',iFontName,...

    'FontSize',iFontSize,'FontWeight','normal');

% 'XAxisLocation','top'设定x轴在顶部

% 'ydir','reverse'设定y轴逆序

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%画加入噪声后的数据%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

subplot(1,4,[4]);

imagesc(z,t,a1);%画加入噪声后的数据

text(iFigureNameXYZ(1), iFigureNameXYZ(2), iFigureNameXYZ(3),'Orignal segy data',...

    'color',iFigureNameColor,'rotation',0,'FontName',iFontName,'FontSize',...

    iFontSize,'FontWeight','normal');

caxis(iDataCaxis);%设定画图时数据的幅值范围

colormap(gray);%设定colormap

xlim([min(z) max(z)]);%设定x轴的坐标范围

ylim([min(t) max(t)]);%设定y轴的坐标范围

set(gca,'xTick',[0:50:max(z)]);%设定x轴的坐标刻度

set(gca,'yTick',[0:0.1:max(t)]);%设定y轴的坐标刻度

colorbar('Ticks',[round(min(iDataCaxis)):4:round(max(iDataCaxis))],...

    'location','southoutside');%设定colorbar的显示刻度及位置

% 设定图题

% title('Input data','FontName',iFontName,'FontSize',iFontSize,'FontWeight','bold');

% 设定x轴的标签

xlabel('Trace number','FontName',iFontName,'FontSize',iFontSize,'FontWeight','normal');

% 设定y轴的标签

ylabel('Time (s)','FontName',iFontName,'FontSize',iFontSize,'FontWeight','normal');

set(gca,'Box','on','XAxisLocation','top','ydir','reverse','FontName',iFontName,...

    'FontSize',iFontSize,'FontWeight','normal');

% 'XAxisLocation','top'设定x轴在顶部

% 'ydir','reverse'设定y轴逆序