Matlab地震勘探Marmousi模型：从数据读取到正演模拟全攻略

简介本资源是一套面向数学建模学习者、地球物理专业学生及科研工程师的MATLAB地震波正演模拟实践材料聚焦经典Marmousi速度模型的构建与波动方程数值求解解决地震勘探中地下结构可视化与算法验证的核心问题。压缩包共13个文件含5幅结果图像jpg用于展示不同时刻波场快照与成像效果5个.dat数据文件存储速度模型与波场快照1个核心marmousi.m主程序实现有限差分法求解二维声波方程1个.fig图形文件保留可交互视图1个说明性txt文档整体大小6.81MB结构紧凑、即开即用。已有1311人学习下载资源提供完整可运行代码、多角度可视化结果及典型地层响应特征帮助读者深入理解边界条件设置、时间步进策略、反射/折射物理机制并为后续偏移成像、反演算法开发奠定建模基础。 很多人在平台下载过“数学建模基于Matlab 地震勘探Marmousi模型”这个源码包编号1977期里面通常有一堆.m文件、一个二进制或文本格式的速度模型还有一个写得比较随意的说明文档。第一次打开的人大概率是懵的不知道先运行哪个脚本不知道生成的地震记录有什么用更不知道Marmousi这三个字母背后到底意味着什么。这篇文章就围绕这个源码包展开。我会从Marmousi模型的来历讲起到数据怎么读、正演模拟怎么跑、炮集怎么处理再到我实际调通代码后踩过的几个坑尽量把一条完整的技术链路串起来。内容适合准备数学建模竞赛的同学、刚接触地震勘探仿真或波动方程数值模拟的研究生也适合那些想在Matlab里找一个复杂地球物理模型练手的人。1. Marmousi模型是什么为什么它能成为地震勘探的“标准试卷”1.1 从勘探难题到数值试验田Marmousi模型最早是法国石油研究院IFP等单位在1988年前后发布的一个二维速度模型最初目的是检验当时各种地震偏移成像算法。它的名字来自Marmousi这个地方后来随着版本迭代逐渐成为地球物理界公认的复杂性测试模型。为什么它能成为“标准试卷”因为真实野外采到的地震资料往往不够干净有噪音、有采集脚印、有近地表复杂因素用来对比不同算法的优劣时很难说是算法本身的问题还是数据质量的问题。Marmousi模型的好处是速度场已知地下构造足够复杂你可以在一个“标准环境”里反复测试各种正演、偏移、反演方法。现在很多数学建模题目虽然不会直接说“用Marmousi”但只要涉及到地震波传播、地下速度结构成像、断层识别等话题队伍都会私下用Marmousi模型生成合成数据再拿这套数据去设计算法、评估结果。这比编一个简单水平层状模型要有说服力得多。1.2 速度模型参数看一眼就知道后边要吃什么苦我拿到这个源码包后第一件事是把速度模型文件的结构摸清楚。最常见的一个公开版本规模是这样的参数数值网格数深度384点 × 水平750点网格间距12.5 m横向和纵向一致模型水平长度约9.375 km模型深度约4.8 km速度范围约1500 m/s 到 5500 m/s典型构造断层、盐丘、薄互层、低速异常体速度从浅层的1500 m/s一路升到深层的5500 m/s区别不是渐变这么简单。模型里有一个大盐丘盐丘和围岩之间的速度差会让地震波路径发生强烈扭曲断层面又把地层切断造成起伏反射界面和绕射波。所以这个模型跑出来的单炮记录远不是几根光滑双曲线而是各种波的叠加非常考验后期处理水平。如果压缩包里的数据尺寸和我说的不一样也别慌。不同版本的Marmousi经过重采样网格间距可能是6.25 m、12.5 m或25 m网格数会随之变化。你只需要在读取后打印一下size(vp)再对照dx、dz参数就能确认自己手里是哪一版。2. 解压源码包后第一件事把速度模型读进Matlab并正确画出来2.1 三种常见数据文件的读取方式源码包里通常有几种形式的速度模型文件读取方式差别很大。第一种是二进制文件后缀常为.dat或.bin。这类文件在Matlab里用fread读取但必须先知道存储的数据类型和维度顺序。一个常见组合是以float32格式存储按深度优先列存储尺寸为384×750。读取代码可以写成fid fopen(marmousi_model.bin, rb); if fid -1 error(找不到模型文件请检查路径); end raw fread(fid, [384, 750], float32); fclose(fid); vp raw; % vp(1:nz, 1:nx)行对应深度列对应水平第二种是文本文件后缀为.txt或.asc。速度值一般用空格或逗号分隔数据量大时读取慢一些但便于确认内容。推荐用readmatrixvp readmatrix(marmousi_model.txt);第三种是.mat文件这个最简单load(marmousi_model.mat);不管用哪种方式读完之后建议先检查数据范围fprintf(min vp %.2f m/s, max vp %.2f m/s\n, min(vp(:)), max(vp(:)));如果最大值在5500左右说明速度单位是m/s可以直接用。如果看到0.5、1.5这种数值说明单位是km/s正演计算前要乘1000否则波速全部慢了1000倍模拟时间完全是错的。2.2 显示模型时最容易犯的三个错误正确显示速度模型是后续所有工作的基础。显示代码本身很简单dx 12.5; dz 12.5; nx size(vp, 2); nz size(vp, 1); x (0:nx-1) * dx / 1000; % 转成 km z (0:nz-1) * dz / 1000; figure; imagesc(x, z, vp); axis equal tight; colormap(jet); colorbar; xlabel(水平距离 (km)); ylabel(深度 (km)); set(gca, YDir, reverse); title(Marmousi 速度模型);这里有几个坑我实际踩过也看周围同学反复踩第一个坑是维度顺序搞反。很多人读完raw后习惯性转置结果vp变成750×384再用imagesc画的时候横坐标对不上水平距离纵坐标对不上深度后面设计炮点和检波器位置全部错位。建议统一约定vp的第1维是深度第2维是水平方向后续正演代码也用这个约定。第二个坑是没有设置YDir reverse。如果不加这一句深度方向会从下往上增大模型看起来就像从地下往地表看和地质上“地表在上、深度向下”的习惯相反。看图容易误判构造位置。第三个坑是坐标单位混用。如果x和z用m显示范围是9000和4800横纵比例会拉得很扁转成km显示配合axis equal tight盐丘和断层的相对位置看起来才符合真实比例。这个显示习惯直接关系到后面解释构造特征。3. 声波正演模拟的原理和Matlab实现3.1 波动方程离散化为什么二阶差分就够入门正演模拟的核心是求解声波方程。二维声波方程写成[ \frac{\partial^2 p}{\partial t^2} c^2 \left( \frac{\partial^2 p}{\partial x^2} \frac{\partial^2 p}{\partial z^2} \right) ]其中p是压力波场c是速度。把这个方程离散化最常用的是时间二阶、空间二阶中心差分。时间二阶的意思是用上一时刻和当前时刻的波场推算下一时刻[ p^{n1} 2p^n - p^{n-1} \left(\frac{c \Delta t}{\Delta x}\right)^2 \left(p_{i1,j} p_{i-1,j} p_{i,j1} p_{i,j-1} - 4p_{i,j}\right) ]这个过程可以类比成水面波纹你往水里丢一颗石子水面相邻位置会依次产生扰动。差分格式做的就是把这种“相邻位置扰动”用网格点的数值关系表达出来。空间二阶差分理解起来最直观代码也最容易写对。虽然它比四阶差分更容易出现频散但作为入门和数学建模完全够用。如果之后追求更好的模拟效果可以把空间差分升级到四阶。四阶差分对Laplace项的近似是[ \nabla^2 p \approx \frac{-p_{i2,j} 16p_{i1,j} - 30p_{i,j} 16p_{i-1,j} - p_{i-2,j}}{12\Delta x^2} ]在同一种网格下四阶差分的波前更干净频散更弱代价是计算量变大而且编程时对边界的索引处理更复杂。我给源码包里加过一版四阶差分实测在12.5 m网格、25 Hz主频条件下效果比二阶好不少。3.2 震源子波、稳定性条件与参数计算震源子波我建议用Ricker子波这是地震勘探里最常用的零相位子波。表达式是[ w(t) \left(1 - 2\pi^2 f_0^2 t^2\right) \exp\left(-\pi^2 f_0^2 t^2\right) ]其中f0是主频我取25 Hz。为什么取25 Hz而不是50 Hz要解释一下。网格间距12.5 m时空间上能有效采样的最短波长一般要求网格间距小于波长的十分之一。浅层速度1500 m/s时25 Hz对应波长约60 m比12.5 m的网格大得多深层速度5500 m/s时波长更大。反过来如果主频取80 Hz浅层已经是高速区也可能出现网格频散波前面会出现明显抖动。所以我的经验是12.5 m网格搭配20到30 Hz主频最稳。时间步长也不能乱选。二维声波方程时间显式格式的CFL稳定性条件近似是[ \Delta t \le \frac{\Delta x}{c_{max} \sqrt{2}} ]代入Δx12.5 mc_max5500 m/s算出来约1.6 ms。这只是一个上限实际为了防止数值误差过大我会取0.5 ms或者1 ms。0.5 ms时一个25 Hz子波的周期40 ms覆盖80个时间采样点波形表示足够细腻总时长2 s等于4000步在Matlab里也就是一次循环的事。3.3 吸收边界用“海绵层”先把边界反射压住做正演模拟时模型四边的截断会造成人为反射这些反射会从边界返回区域内部污染真实的反射波。如果不加处理单炮记录里会出现一堆从边框冒出来的平行同相轴跟Marmousi模型本身的复杂波场混在一起。最标准的做法是PML完美匹配层效果最好但代码量大对于数学建模场景有点杀鸡用牛刀。退一步可以用海绵吸收层也就是在模型四周设置一定厚度的衰减带每走一步把边界附近的波场乘一个小数让波碰到外边界前先被“吸收”掉。这个办法虽然不如PML干净但胜在代码短、容易调。海绵层厚度我习惯设为40到60个网格。衰减系数设计成从内到外逐渐增强内边界附近几乎不衰减最外层每步衰减80%左右。这样既不会在内外分界处形成新反射又能把大部分出射波吃掉。3.4 可运行的完整正演循环代码下面这段代码是我基于源码包整理出的最小可运行版本噪声尽量少逻辑尽量直白% 参数设置 dx 12.5; dz 12.5; dt 0.5e-3; nz size(vp, 1); nx size(vp, 2); nt 4000; % 总时间步2 s % 波场初始化 p zeros(nz, nx); pold zeros(nz, nx); pnew zeros(nz, nx); % 雷克子波 f0 25; tw (-2/f0 : dt : 2/f0); src (1 - (pi*f0*tw).^2) .* exp(-(pi*f0*tw).^2); nsrc length(src); % 海绵吸收层 nbound 40; damp ones(nz, nx); for k 1:nbound beta 0.8 * ((nbound - k 1) / nbound)^2; d 1 - beta; damp(k, :) damp(k, :) * d; damp(nz-k1, :) damp(nz-k1, :) * d; damp(:, k) damp(:, k) * d; damp(:, nx-k1) damp(:, nx-k1) * d; end % 震源与检波器位置 sx 375; sz 60; recIdx 10:20:700; nrec length(recIdx); rec zeros(nt, nrec); % 时间递推 for it 1:nt if it nsrc p(sz, sx) p(sz, sx) src(it); end % 空间二阶差分 lap zeros(nz, nx); lap(2:nz-1, 2:nx-1) p(1:nz-2, 2:nx-1) p(3:nz, 2:nx-1) ... p(2:nz-1, 1:nx-2) p(2:nz-1, 3:nx) - ... 4 * p(2:nz-1, 2:nx-1); % 时间递推 pnew 2*p - pold (vp * dt / dx).^2 .* lap; % 边界置零 海绵吸收 pnew(1, :) 0; pnew(nz, :) 0; pnew(:, 1) 0; pnew(:, nx) 0; pnew pnew .* damp; % 检波器记录 rec(it, :) p(sz, recIdx); % 更新波场 pold p; p pnew; end跑完之后rec就是一个二维矩阵横轴是检波器序号纵轴是时间采样点。这个矩阵就是所谓的地震炮集记录。注意我把检波器放在和震源同一深度这样既能看下行波也能看反射波如果你想要贴近地表的记录把recIdx的取值改成p(2, recIdx)并去掉震源附近的强能量干扰即可。这段代码没有做太多性能优化但胜在可直接运行。用这个思路去处理1887行的模型4000步循环在普通电脑上要跑一段时间。我建议先把nx和nz截小一半做个快速验证确认逻辑没问题再全模型跑。4. 炮集记录生成后如何往下做叠加与成像4.1 从单炮记录里能看出哪些波单炮记录生成后先画出来看看figure; imagesc((recIdx-1)*dx/1000, (0:nt-1)*dt*1000, rec); xlabel(偏移距 (km)); ylabel(时间 (s)); set(gca, YDir, reverse); colormap(gray);Marmousi模型的单炮记录非常热闹。你会看到几类波。第一类是直达波出现在记录最上面近似一条斜率固定的斜线由震源沿地表附近直接传播到检波器本文还有配套的精品资源点击获取