目录
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。
数据集 | SEGSaltData | SimulateData |
数量(train+test) | 140(130+10) | 1700(1600+100) |
地震数据尺寸 | 400 x 301 | 201 x 301 |
速度模型尺寸 | 400 x 301 | 201 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)地震数据采集后,到目前拿到手上的数据,中间是否经过了别的处理?这种处理是否是加入噪声,或者减少信息量的。