简介：一套面向本硕博教研人群的线性地震反演Matlab仿真资源，聚焦高斯混合马尔科夫-蒙特卡洛（GM-MCMC）算法的编程实现与原理验证。资源包共13个文件，压缩后约1.9MB，其中包含9个M脚本/函数、2个MAT数据文件、1个TXT说明文档以及1个AVI操作录像，源码按算法核心流程拆分为Metropolis采样、协方差矩阵计算、后验概率求解、转移矩阵处理等多个可独立调用的模块，方便学习者逐步读懂并复用。配套的操作录像演示了在Matlab2021a及以上版本中的运行步骤，尤其强调通过Runme.m主脚本启动、避免直接运行子函数等关键注意事项，能有效降低上手门槛。目前已有500人学习浏览，适合需要结合完整仿真实例快速掌握GM-MCMC线性地震反演方法的高校师生和科研人员。
1. 高斯混合马尔可夫链蒙特卡洛：线性地震反演不再只给一个解
地下介质参数对叠前地震数据的响应在常规反演里往往被压成一个最优解，但真实生产中最有用的不是那个点，而是它周围的分布。这个仿真工程把正演模型、高斯混合先验和MCMC采样串在一起，用matlab实现了从合成三层油水介质生成地震道，到采样后验弹性参数的完整闭环。对研究生和刚接触反演的工程师来说，它的价值不在于公式推导，而在于你能看到每个环节的矩阵和概率密度到底怎么衔接。运行主脚本Runme.m，配合操作录像，可以快速复现GM-MCMC的完整行为。适合想搞懂算法边界、又不想把反演当黑箱的人。
2. 线性地震反演的正演链路：从弹性参数到合成记录
2.1 线性化近似是GM-MCMC工作的前提
地震反演里最容易翻车的是把非线性问题直接塞进线性MCMC框架。Zoeppritz方程能准确描述反射振幅与入射角、纵横波速度、密度之间的关系，但对每个MCMC提议都计算一遍Zoeppritz，成本高，而且接受率会随着角度道集的增加明显下降。这个工程选择Aki-Richards三参数线性近似，把正演过程写成d = G*m + w，这和GM-MCMC的采样逻辑正好咬合。GM-MCMC的核心是计算后验比，线性正演算子G可以预先建好，整个采样过程只需要做矩阵乘法，不存在每步重新推导雅可比矩阵的问题。
这里容易误解的是“线性反演”不等于模型简单。data_synth_3layers_oil_water.mat包含的是含油、含水和盖层三层模型。三层介质对弹性参数的需求是三个高斯混合分量：泥岩、含水砂岩、含油砂岩。油层的纵波速度往往比水层低，密度也会变化，但横波速度变化不如纵波那么明显，因此高斯混合先验能把这种多峰相关性表达出来。
2.2 正演模型的MATLAB实现
elasticForwardModel.m是这个工程的第一个关键点。它把vp、vs、rho的纵向剖面转换成反射系数，再与子波卷积，最后输出角度道集和线性算子。按常见做法，Aki-Richards近似可以写成如下形式：
function [syn, G, rc] = elasticForwardModel(vp, vs, rho, theta, wavelet) % 把三层模型扩展成界面反射系数序列 % vp, vs, rho 长度为 nLayer 的向量 % theta 是角度道集的角度列表，单位是度 theta_r = theta(:) * pi / 180; % 转弧度 nTheta = length(theta_r); % 计算界面处的上、下层平均 vpA = 0.5 * (vp(1:end-1) + vp(2:end)); vsA = 0.5 * (vs(1:end-1) + vs(2:end)); rhoA = 0.5 * (rho(1:end-1) + rho(2:end)); dvp = vp(2:end) - vp(1:end-1); dvs = vs(2:end) - vs(1:end-1); drho = rho(2:end) - rho(1:end-1); % 反射系数矩阵 nReflector x nTheta rc = zeros(length(vpA), nTheta); for it = 1:nTheta rc(:, it) = 0.5 * (1 + tan(theta_r(it))^2) .* (dvp ./ vpA) ... - 4 * (vsA ./ vpA).^2 * sin(theta_r(it))^2 .* (dvs ./ vsA) ... + 0.5 * (1 - 4 * (vsA ./ vpA).^2 * sin(theta_r(it))^2) .* (drho ./ rhoA); end % 每个角度道分别卷积子波 syn = zeros(size(rc)); for it = 1:nTheta syn(:, it) = conv(rc(:, it), wavelet(:), 'same'); end这段代码的物理含义是：纵波速度的相对变化控制近偏移距振幅，横波速度和密度则在远偏移距部分起作用。用dvp ./ vpA而不是直接用绝对速度，能够避免不同区块背景速度差异带来的数值尺度问题，也让MCMC在采样vp时更容易控制步长。计算G矩阵时，常见做法是将子波构造成Toeplitz卷积矩阵，把卷积操作并入线性算子，从而得到d = G * m。这样GM-MCMC里的后验似然p(d|m)可以直接用G乘当前模型参数得到。
需要注意，Aki-Richards近似在入射角超过30度以后误差会明显增大。对叠前道集，工程里一般保留30度以内的道集，否则后验估计会偏差。这个仿真包只处理线性正演，因此不要试图往theta里填50度这样的角度数据。
2.3 协方差矩阵怎么放进反演
反演不是把反射系数求出来就结束，后验概率密度还需要协方差矩阵来定义距离。covariance_matrix_exp.m构造的就是这类矩阵，常见形式是平方指数核：
function K = covariance_matrix_exp(x, sigma, l) % 平滑型协方差矩阵 % x : 采样点坐标 % sigma : 幅度，控制参数波动范围 % l : 相关长度，控制平滑程度 n = length(x); K = zeros(n, n); for i = 1:n for j = 1:n K(i, j) = sigma^2 * exp(-(x(i) - x(j)).^2 / (2 * l^2)); end end这段实现虽然用了双重循环，但n一般不会超过500，在matlab里直接跑也没有太大性能压力。l是最值得调的一个参数。l太小，协方差矩阵接近对角阵，反演结果会高频抖动；l太大，矩阵几乎变成常值矩阵，约束过强，薄层信息会被抹平。对三层含油水模型，可以把l设置为一个采样间隔的5到10倍，即让上下两层之间的相关性自然衰减。sigma则参照目标参数的标准差设置，例如纵波速度变化在5%左右，sigma就取背景速度的0.05倍。协方差矩阵用在两个地方：一是似然里的噪声协方差，另一个是高斯混合先验中每个分量的协方差。两者混用时，MATLAB里的变量名经常都是C或Sigma，运行前要看清楚是给哪个环节用的。
3. GM-MCMC采样：把多峰先验注入马尔可夫链
3.1 为什么高斯混合先验比单高斯先验合适
如果只有一层干净的砂岩，用单高斯先验p(m)=N(m|mu, Sigma)就够了。但真实剖面是层状介质：泥岩、含水砂岩、含油砂岩的速度和密度存在多个中心，整个模型空间的分布并不满足单峰假设。高斯混合先验把模型参数先验写成：
p(m) = sum_k w_k * N(m | mu_k, Sigma_k)
每一个k对应一种岩相或流体状态，权值w_k不是随意给的，它表示该分量在剖面上出现的先验比例。对一个三层含油水模型，k=1代表泥岩盖层，k=2代表含水砂岩，k=3代表含油砂岩。这个仿真包里的GaussianMixMCMC_metropolis.m就是围绕这个先验搭建的。后验概率写成：
posterior(m) proportional to likelihood(d | m) * prior(m)
由于先验是多峰的，后验也会是多峰。传统最优化方法容易陷入其中一个局部极大值，而MCMC的目标是生成符合后验分布的样本，所以链有机会在几个峰值之间移动。不过MCMC本身并不能自动解决模态跳跃问题，这就要看提议分布和状态转移矩阵设计得怎么样。
3.2 Metropolis-Hastings 采样器构建
Metropolis-Hastings是这个包的核心采样器。给定当前参数m_cur，先从提议分布中随机抽取m_prop，然后按接受概率决定是否接受。常见的简化版本是随机游走提议：
for iter = 1:nIter % 用预先做好的Cholesky分解生成相关随机扰动 z = randn(nParam, 1); m_prop = m_cur + stepScale * (L_prop * z); % 计算后验对数比 logPost_prop = computeLogPosterior(m_prop, d, G, noiseSigma, gmPrior); logPost_cur = computeLogPosterior(m_cur, d, G, noiseSigma, gmPrior); logAlpha = min(0, logPost_prop - logPost_cur); if log(rand()) < logAlpha m_cur = m_prop; end chain(:, iter) = m_cur; end计算logPosterior时需要把似然、先验都写在log域里。高斯混合先验要用log-sum-exp，避免多个小概率分量乘在一起后直接下溢成0。calculate_posterior_probability.m做的就是这个工作，函数内部会调用covariance2correlation.m来检查分量协方差矩阵是否病态，必要时会把协方差转换为相关系数做诊断。参数stepScale对采样影响极大。stepScale太大，提议经常跑到先验的低概率区域，接受率低；stepScale太小，接受率高但链像蜗牛一样爬，需要很长的链才能覆盖整个后验空间。工程实践中先跑500次试验，把接受率控制在0.2到0.4之间，再决定正式链长。
3.3 转移矩阵与隐状态切换
GM-MCMC比普通MCMC多出来的部分就是岩性或流体状态切换。transpose_trasition_matrix.m和simulate_markov_chain.m负责生成状态序列。拿三分量岩相模型举例，状态转移矩阵T是一个3乘3矩阵，T(i,j)表示当前样本处于分量i时，下一步转移到分量j的概率。生成状态序列的常见做法是：
K = 3; T = [0.8 0.15 0.05; 0.1 0.8 0.10; 0.1 0.1 0.80]; % 每行相加为1 state = zeros(1, nStep); state(1) = randsample(1:K, 1, true, [0.4 0.3 0.3]); for t = 2:nStep state(t) = randsample(1:K, 1, true, T(state(t-1), :)); end这一段模拟的是先验标签的游走过程。在有些实现里，状态序列会作为条件先验的标签参与每个参数的采样；在另外一些实现里，状态转移矩阵被转置后作用于MCMC的提议。文件名中的transpose_trasition_matrix.m，拼写保留了原包的trasition，调用时不要把名字写错。设计T时要注意对角元素不能太接近1，否则链永远停留在同一岩相；也不能过小，否则频繁切换会让混合模型退化成独立采样。对三层模型，对角线取0.8到0.9通常是比较稳的起点。
4. 反演工程落地：数据准备、运行顺序与参数调优
4.1 拿到压缩包后的运行顺序
这个仿真包不是把一堆函数堆在一起就完事，正确顺序是先运行根目录下的主脚本Runme.m，再按需查看函数。主脚本会加载data_synth_3layers_oil_water.mat和cmaps.mat，生成正演记录、初始化GM-MCMC参数、调用GaussianMixMCMC_metropolis.m并画出结果图。直接点击子函数文件运行很容易报“函数未定义”或“找不到变量”，因为很多中间变量只在主脚本中存在。
注意：打开MATLAB后必须把当前文件夹窗口切换到工程所在路径，再运行
Runme.m。即使脚本文件已经在编辑器中打开，当前路径不对也会导致cmaps.mat等数据文件加载失败。
工程里几个核心文件的作用可以先用下面这张表理解：
| 文件 | 在整个反演管线中的角色 |
|---|---|
Runme.m | 主入口，控制数据加载、参数设置、输出结果图 |
elasticForwardModel.m | 从弹性参数正向计算合成角度道集 |
covariance_matrix_exp.m | 构造先验/噪声协方差矩阵 |
GaussianMixMCMC_metropolis.m | MH采样核心，更新模型参数 |
simulate_markov_chain_finalDefined.m | 生成采样链，处理马尔可夫状态序列 |
transpose_trasition_matrix.m | 转置状态转移矩阵 |
calculate_posterior_probability.m | 在log域计算后验概率 |
covariance2correlation.m | 协方差转相关，诊断变量耦合 |
操作录像操作录像0023.avi会演示完整的点击过程。第一次跑的时候先看录像，确认主脚本名字和当前路径，能省掉一半的报错时间。
4.2 参数怎么设：从三层模型出发
三层含油水模型不是随便给一个高斯混合参数就能跑。先看岩石物理趋势：泥岩的纵波速度一般比含水砂岩偏高或接近，含油砂岩因为含流体性质不同，纵波速度会相对低一些，密度也可能降低。下面给一张示例性质的初始值表，具体数值要以压缩包内data_synth_3layers_oil_water.mat中的真实变量为准：
| 岩相 | vp范围(m/s) | vs范围(m/s) | 密度范围(g/cm3) | 先验权值 |
|---|---|---|---|---|
| 泥岩 | 3400~3800 | 1700~2000 | 2.30~2.50 | 0.4 |
| 含水砂岩 | 3200~3600 | 1600~1900 | 2.20~2.40 | 0.3 |
| 含油砂岩 | 3000~3400 | 1500~1800 | 2.15~2.35 | 0.3 |
MH采样器不会直接采样整个速度剖面，而是采样模型参数的相对扰动量。这样协方差矩阵的sigma取相对变化更合理。链长方面，三层模型参数维度不高，1万次采样够用，去掉前2000次作为burn-in，再每隔5个样本保留一个，大约得到1600个有效样本。若看到后验均值还不稳定，先别急着加链长，优先调整提议步长。
4.3 后验统计与收敛诊断
采样完成后，要把后验样本转成可用结果，同时检查链是否收敛。下面是常用的一段处理代码：
burnin = 2000; thin = 5; chainAfterThin = chain(:, burnin+1:thin:end)'; postMean = mean(chainAfterThin, 1); postStd = std(chainAfterThin, 0, 1); postP95 = quantile(chainAfterThin, [0.025 0.975], 1); figure; subplot(2,1,1); plot(chain(1, :)); title('第一维参数的采样轨迹'); subplot(2,1,2); autocorr(chain(1, burnin+1:end));观察轨迹图时，如果变量在某一个值附近长时间徘徊，说明链可能被困在某个局部峰。如果轨迹呈现明显的分段跳跃，说明GM混合分量之间的转移矩阵起作用了。covariance2correlation.m把采样协方差转成相关系数矩阵，可以用来检查vp和rho之间是否高度相关。线性反演中纵横波速度与密度存在固有耦合，相关系数超过0.95时，反演结果的可解释性要打折扣，这时需要固定一个参数或增加先验约束，而不是去调MCMC步长。
5. 验证与绕坑：GM-MCMC仿真前先做三件事
5.1 先跑短链，看变量范围
不要一上来就跑10万次迭代。先用1000次迭代跑通流程，打印模型参数的最小、最大和均值。如果参数在迭代50步后直接飞到1e6量级，问题通常不在地质模型，而在协方差矩阵的数值稳定性或提议分布的尺度上。
5.2 协方差矩阵加一个小的jitter
协方差矩阵exp核在采样点间距过小或l过大时，容易变成数值上接近奇异的矩阵。Cholesky分解时会直接报错。常见做法是给对角加一个1e-6倍的单位阵：
K = K + 1e-6 * eye(n); [L, p] = chol(K, 'lower'); if p ~= 0 warning('协方差矩阵非正定，增大jitter或减小相关长度'); endchol分解成功后再进入MCMC循环。这样既避免了重复分解的开销，也能让提议分布保持正确。
5.3 文件名和路径是最大的坑
包里transpose_trasition_matrix.m里的trasition是原始拼写，调用时保持一致即可，改成transition反而会让主脚本找不到文件。操作录像里最值得留意的一步，是启动MATLAB后先设置当前文件夹。只要当前路径是工程根目录，cmaps.mat和data_synth_3layers_oil_water.mat都能被相对路径找到。若在Windows下用脚本自动跑，建议在Runme.m一开始加入：
cd(fileparts(mfilename('fullpath')));这一行会把工作目录强制切到主脚本所在目录，避免双击打开脚本但路径不对的问题。mfilename只有在脚本文件未运行前使用才有效，放到Runme.m第一行没有问题。
最后再强调一个容易忽略的操作：GM-MCMC的接受率计算建议始终在log域完成，不要把概率密度的原始值乘起来再取对数。calculate_posterior_probability.m内部已经处理了log-sum-exp，自己写子函数时不要贪省事直接用prod。对于三层线性地震反演模型，把接受率稳定在0.3附近后再放大链长，所生成的后验区间才值得往下游解释。
本文还有配套的精品资源，点击获取
