FWI 地震数据的认识

news/2024/12/2 9:44:07/

目录

1、数据来源。

1)SEG  系列。

2)OpenFWI 系列。

2、数据量,尺寸。

1) SEG 包含两个数据集:SEGSaltData 和 SimulateData。

2)OpenFWI 包含12个数据集:

3、地震数据对应的观测系统。

1) SEG系列

2)OpenFWI系列

4、显示数据的源码

5、正演的原理及源码

6、我的疑问:


地震数据是非常宝贵的资源,很多真实数据并不是公开的,目前在网上流传的都是合成数据。我将从以下角度来介绍合成数据:

1、数据来源。

1)SEG  系列。

数据源自论文: Yang F, Ma J. Deep-learning inversion: A next-generation seismic velocity model building method DL for velocity model building[J]. Geophysics, 2019, 84(4): R583-R599.

作者基于unet架构,实现了端到端的全波形反演,从地震数据直接获得速度模型,并在盐体数据上测试性能。

论文代码:GitHub - YangFangShu/FCNVMB-Deep-learning-based-seismic-velocity-model-building: Deep-learning inversion: A next-generation seismic velocity model building method

2)OpenFWI 系列。

数据源自论文:Deng C, Feng S, Wang H, et al. OpenFWI: Large-scale Multi-structural Benchmark Datasets for Full Waveform Inversion[J]. Advances in Neural Information Processing Systems, 2022, 35: 6007-6020.

数据和代码网址: https://openfwi-lanl.github.io/

它涵盖了地球物理的多个领域(界面、断层、CO2储层等),涵盖了不同的地质地下构造(平面、曲线等),包含了不同数量的数据样本(2K - 67K),还包括3D FWI数据集。此外,比较了三种二维FWI方法:InversionNet[20]提出了一种全卷积网络来模拟地震反演过程;VelocityGAN[27]采用基于gan的模型求解FWI;UPFWI[31]将正演建模与CNN连接在一个回路中,实现无监督学习,无需地面真值速度图进行训练。还有一种三维FWI方法:InversionNet3D[30]将InversionNet扩展到三维领域。为了减少内存占用和提高计算效率(即3D反演中最具挑战性的两个障碍),该网络在编码器中使用了群卷积,并通过基于加性耦合的可逆层采用了部分可逆架构[46]。

2、数据量,尺寸。

1) SEG 包含两个数据集:SEGSaltData 和 SimulateData。

数据集SEGSaltDataSimulateData
数量(train+test)140(130+10)1700(1600+100)
地震数据尺寸400 x 301201 x 301
速度模型尺寸400 x 301201 x 301

2)OpenFWI 包含12个数据集:

其中第二行的 B系列,比第一行的A系列有更高的难度。从左到右,依次为平面(FlatVel-A、FlatVel-B)、曲面(CurveVel-A 、CurveVel-B)、平面断层(FlatFault-A、 FlatFault-B)、曲面断层(CurveFault-A、CurveFault-B)、复杂的野外合成数据(Style-A、Style-B)、二氧化碳储层(Kimberlina-CO2)、三维地震数据(3D Kimberlina-V1)。

在这里要注意,理解下表的尺寸时,如速度模型70x1x70, 中间为1表明是个二维速度模型。

3、地震数据对应的观测系统。

1) SEG系列

 参考星移论文

2)OpenFWI系列

4、显示数据的源码

我的数据都存在.mat文件中:

1)显示地震地震数据:

2)显示速度模型:

5、正演的原理及源码

原理及波动方程公式后续补上。


% 对文件夹内所有文件进行计算并存储
% function [] = forward(vdir,vname,gdir,gname,start_num,end_num)
%     for i = start_num:end_num
%         vfile = [vdir,'/',vname,num2str(i),'.mat']
%         gfile = [gdir,'/',gname,num2str(i),'.mat']
%         calc(vfile,gfile);
%     end
% end% 对单个速度模型文件计算多炮数据
function [] = forward(vfile,outfile,sn) %sn 炮数
load(vfile,'vmodel');  %加载vfile文件中的vmodel变量
[nz,nx] = size(vmodel);%  201*301
Rec = ones(400,nx,sn);%  400*301 是地震数据的尺寸%tic用来保存当前时间,而后使用toc来记录程序完成时间,两者往往结合使用
tic;
for i = 1:snsx = i*fix(nx/sn);  % 这一步的意思是放炮的位置是均匀放置的么?singleRec = calc(vmodel, sx);Rec(:,:,i) = singleRec;
end
toc;save(outfile,'Rec');
end% 对单炮速度模型进行计算
function singleRec = calc(vmodel, sx)[nz,nx] = size(vmodel); % x方向网格数(从左到右 和 z方向网格数(从上到下) 201*301 
npmlz = 10;            	% 顶部和底部的PML层厚度
npmlx = 10;             % 左边和右边的PML层厚度
% sx = 100;                % 震源的x方向网格数
sz = 0;                % 震源的z方向网格数 表明震源在地表
dx = 10;                 % x方向的网格大小 单位:米
dz = 10;                 % y方向的网格大小 单位:米
nt = 2000;               % 计算的总时间步数   在方程的最后下采样5,最后求的400 对应301*400的地震数据尺寸
dt = 1e-3;            	% 单位时间步长度 单位:秒
ampl = 1.0e0;           % 震源子波的震幅
xrcvr = 1:1:nx;         % 地面上x方向的接收器位置
nodr = 3;               % half of the order number for spatial difference.B = [1 zeros(1,nodr - 1)]';
A = NaN*ones(nodr,nodr);
for i = 1:1:nodrA(i,:) = (1:2:2*nodr - 1).^(2*i - 1);
end
C = A\B; Nz = nz + 2*npmlz;
Nx = nx + 2*npmlx;rho = 1000*ones(Nz,Nx);                                                        % 密度; Unit: kg/m^3.  密度对结果有什么影响? 此处密度1000
rho(fix(Nz/3):end,fix(Nx/2):end) = 500;                                        % 通过修改密度,可以模拟介质中存在不同的物理特性,如不均匀的岩层、气体或液体的存在等
vp = padarray(vmodel,[10,10],'replicate','both');                           % 扩展边界% 边界扩展是为了处理波场模拟中的边界效应,确保在模拟过程中声波传播不会超出vmodel的范围。% 通过在vmodel的周围添加额外的边界值,可以避免波场反射和干扰。padarray函数的'replicate'选项表示将vmodel的边界值复制并填充到扩展后的边界上。f0 = 25;                                                                    % 震源频率 单位 Hz
t0 = 1/f0;                                                                  % the time shift of source Ricker wavelet; Unit: s; Suggest: 0.02 if fm = 50, or 0.05 if fm = 20.
t = dt*(1:1:nt);
src = (1 - 2*(pi*f0.*(t - t0)).^2).*exp( - (pi*f0*(t - t0)).^2);           	% the time series of source wavelet.   雷克子波
% The source wavelet formula refers to the equations (18) of Collino and Tsogka, 2001.%% Perfectly matched layer absorbing factor PML层吸收因子设置
dpml0z = 3*max(vp(:))/dz*(8/15 - 3/100*npmlz + 1/1500*npmlz^2);
dpmlz = zeros(Nz,Nx);
dpmlz(1:npmlz,:) = (dpml0z*((npmlz: - 1:1)./npmlz).^2)'*ones(1,Nx);
dpmlz(npmlz + nz + 1:Nz,:) = dpmlz(npmlz: - 1:1,:);
dpml0x = 3*max(vp(:))/dx*(8/15 - 3/100*npmlx + 1/1500*npmlx^2);
dpmlx = zeros(Nz,Nx);
dpmlx(:,1:npmlx) = ones(Nz,1)*(dpml0x*((npmlx: - 1:1)./npmlx).^2);
dpmlx(:,npmlx + nx + 1:Nx) = dpmlx(:,npmlx: - 1:1);
% The PML formula refers to the equations (2) and (3) of Marcinkovich and Olsen, 2003.%% Wavefield calculating  依据广义胡克定律求的弹性系数% rho1 和 rho2 是密度(rho)的副本。它们用于计算波场的系数,与波场更新方程中的密度项有关。
% Coeffi1 和 Coeffi2 是沿 x 轴和 z 轴方向的PML吸收因子的系数。它们根据PML吸收因子(dpmlx 和 dpmlz)和时间步长(dt)计算得出。这些系数在波场更新方程中用于考虑PML的吸收效果。
% Coeffi3 和 Coeffi4 是与密度(rho)和空间步长(dx 和 dz)相关的系数,用于考虑波场更新方程中的空间导数项。
% Coeffi5 和 Coeffi6 是与密度(rho)和速度(vp)的平方以及空间步长(dx 和 dz)相关的系数,用于考虑波场更新方程中的速度和应力项。
% 这些系数是根据弹性介质的性质和PML吸收层的影响,结合波场更新方程中的相关项,计算得出的。它们的计算方式是基于广义胡克定律和对介质参数的离散化模拟。
% 通过使用这些系数,可以准确地更新波场的值,模拟弹性波在介质中的传播和衰减行为。rho1 = rho;             % or = [(rho(:,1:end - 1) + rho(:,2:end))./2 (2*rho(:,end) - rho(:,end - 1))];
rho2 = rho;             % or = [(rho(1:end - 1,:) + rho(2:end,:))./2; (2*rho(end,:) - rho(end - 1,:))];Coeffi1 = (2 - dt.*dpmlx)./(2 + dt.*dpmlx);
Coeffi2 = (2 - dt.*dpmlz)./(2 + dt.*dpmlz);
Coeffi3 = 1./rho1./dx.*(2*dt./(2 + dt.*dpmlx));
Coeffi4 = 1./rho2./dz.*(2*dt./(2 + dt.*dpmlz));
Coeffi5 = rho.*(vp.^2)./dx.*(2*dt./(2 + dt.*dpmlx));
Coeffi6 = rho.*(vp.^2)./dz.*(2*dt./(2 + dt.*dpmlz));% +++++++++++++++++++++++++++++++++++++ approximate coeffient ++++++++++++++++++++++++++++++++++++++
% Coeffi1 = 1 - dt.*dpmlx;
% Coeffi2 = 1 - dt.*dpmlz;
% Coeffi3 = 1./rho./dx.*dt;
% Coeffi4 = 1./rho./dz.*dt;
% Coeffi5 = rho.*(vp.^2)./dx.*dt;
% Coeffi6 = rho.*(vp.^2)./dz.*dt;
% --------------------------------------------------------------------------------------------------NZ = Nz + 2*nodr;                                                           % All values of the outermost some columns are set to zero to be a boundary condition: all of wavefield values beyond the left and right boundary are null.
NX = Nx + 2*nodr;                                                           % All values of the outermost some rows are set to zero to be a boundary condition: all of wavefield values beyond the top and bottom boundary are null.Znodes = nodr + 1:NZ - nodr;
Xnodes = nodr + 1:NX - nodr;
znodes = nodr + npmlz + 1:nodr + npmlz + nz;
xnodes = nodr + npmlx + 1:nodr + npmlx + nx;
nsrcz = nodr + npmlz + sz;
nsrcx = nodr + npmlx + sx;Ut = NaN*ones(NZ,NX);                                                     % the wavefield value preallocation. 存储波场的值,并在时间步进过程中进行更新
Uz = zeros(NZ,NX);                                                        % The initial condition: all of wavefield values are null before source excitation. 波场在z方向上的速度分量初始条件
Ux = zeros(NZ,NX);                                                        % The initial condition: all of wavefield values are null before source excitation. 波场在x方向上的速度分量初始条件
Vz = zeros(NZ,NX);                                                        % The initial condition: all of wavefield values are null before source excitation. 波场在z方向上的位移分量初始条件
Vx = zeros(NZ,NX);                                                        % The initial condition: all of wavefield values are null before source excitation. 波场在x方向上的位移分量初始条件
Psum = NaN*ones(Nz,Nx);                                                   % 存储波场在不同时间步长上的压力分量的累积和U = NaN*ones(nz,nx,nt);                                                   % 用于存储完整的波场数据% tic;
for it = 1:1:ntUx(nsrcz,nsrcx) = Ux(nsrcz,nsrcx) + ampl*src(it)./2;Uz(nsrcz,nsrcx) = Uz(nsrcz,nsrcx) + ampl*src(it)./2;Ut(:,:) = Ux(:,:) + Uz(:,:);U(:,:,it) = Ut(znodes,xnodes);
end
% toc;syngram(:,:) = U(1,xrcvr,:);
singleRec = syngram';
singleRec = downsample(singleRec,5); % 下采样end% 实验

6、我的疑问:

1)评价指标,目前论文SSIM、MAE、RMSE,   地质勘探的角度,有什么评价指标?尤其在自然数据中,工程应用场景下,更希望通过FWI提供什么信息?

2)盐体数据中,地层很薄,对这类的识别是否很重要?

3)地震数据的噪声来源:采集设备、地理环境,哪些噪声对地震数据的影响特别大?

4)地震数据采集后,到目前拿到手上的数据,中间是否经过了别的处理?这种处理是否是加入噪声,或者减少信息量的。


http://www.ppmy.cn/news/498795.html

相关文章

CSPM项目管理专业人员能力评价证书报考条件与考试形势

2021年10月,中共中央、国务院发布的《国家标准化发展纲要》明确提出构建多层次从业人员培养培训体系,开展专业人才培养培训和国家质量基础设施综合教育。建立健全人才的职业能力评价和激励机制。由中国标准化协会(CAS)组织开展的项…

如何选显卡

1. NVIDIA Tesla(k40/k80) 专门用来做深度学习 2.Quadro (专门用来做图像设计的显卡) 3. GeForce (专业游戏显卡) 那么问题来了,GeForce 能做深度学习吗? 答案是你自己用可以,但是不能在机房用&#xff0c…

matlab外接显卡,联想发布首款外接显卡坞,可让笔记本性能暴增

描述 (文章来源:驱动之家) 1月6日消息,CES 2020展会期间,联想发布了首款外接显卡坞Legion BoostStation,进军eGPU市场。 此前,Razer已经发布Core X系列eGPU扩展坞,可搭配轻薄笔记本使用,增强后者…

为什么amd显卡便宜却买的人少_为什么这三张显卡没人用? 性能高居榜首, 却无人问津? 网友: 我都没见人提过...

最近峰哥网上冲浪时,看到网友提问,图中AMD的三款显卡怎么没有人用甚至都很少有人提起。 其实根据AMD这3款显卡在显卡天梯图中的位置,就可以说明它们是旗舰类型的显卡,它们的性能自然不必多说,AMD显卡的旗舰巅峰性能。峰…

esxi能直通的显卡型号_最便宜能高画质“吃鸡”的显卡是这个型号:RX470矿卡!...

随着时间的推移,PC端的“吃鸡”游戏玩家在不断减少,但是依然有一群庞大的游戏迷在玩这款游戏,对于预算有限的玩家来说,想要吃鸡最重要的就是更换电脑显卡配件,如果把成本降到最低,只能买矿卡了,…

为什么amd显卡便宜却买的人少_为什么不推荐人选择AMD?

最近真的是被AMD坑的不行,很影响视频的制作时间了,如果可以的话我等10代酷睿出了之后ryzen卖掉上Intel平台。 内存兼容性 这个有一说一,我遇到了很多次,我自己手上有两组不同批次的金士顿骇客神条 8G 2666的内存,一组是…

C盘文件恢复怎么做?数据恢复,就看这4招!

我一般比较重要的文件都会保存到c盘中。最近电脑有点卡顿,想清理一下不需要的文件,但不小心删除了一个很重要的文件,c盘删除的文件还能恢复吗?谁可以帮我想想c盘中的文件如何恢复呢? C盘对于电脑来说是个很重要的磁盘&…

为什么amd显卡便宜却买的人少_买完3080都喊亏!AMD的新显卡用价格砸懵了所有人。。。...

AMD 今天的显卡发布会,逆天了! 可能有小伙伴不太清楚这句话是什么意思,托尼先给大家来点儿前情提要: 作为一家同时研发 CPU 和显卡的科技企业, AMD 多年以来一直被行业里的两大巨头打压: CPU 的性能不如专攻此道的英特尔,显卡性能也不如行业老大英伟达。 虽然从账面上看…