简介针对一维大地电磁反演问题的完整MATLAB实现基于可逆跳跃马尔科夫链蒙特卡洛rjMCMC方法面向从事地球物理反演研究的研究生、科研人员及工程师。这种算法能够在不同维数参数空间间跳转自适应地探索层数、厚度与电阻率等模型结构有效应对大地电磁反演中的多解性与模型复杂度不确定问题。程序涵盖从似然函数、先验分布、协方差生成、提议扰动到链收敛评估与后验统计的整套流程并附有测试模型、常用绘图脚本及README说明便于二次开发与结果复现。压缩包共24个文件核心为14个MATLAB源代码.m辅以8个用于样例与中间结果的.mat数据文件以及1个说明文档和1个开源许可证整体大小7.59MB。内容组织紧凑适合作为对比算法或教学示例已有253人学习浏览。借助该代码可在不深入推导复杂采样公式的前提下完成rjMCMC反演实验并可针对不同观测数据调整模型参数快速获得包含模型不确定性的后验反演结果。1. 一维大地电磁反演最头疼的不是深部而是层数rjMCMC 把这个问题交给数据自己回答做过一维 MT 反演的人都有过这种体验给定同一组视电阻率和相位曲线固定 3 层模型和固定 5 层模型反演出来的电阻率结构完全两样而拟合指标却差不了多少。层数成了先验里最主观、最影响结论的参数。这个压缩包标题里的核心词rjMCMCreversible-jump Markov Chain Monte Carlo可逆跳跃马尔科夫链蒙特卡洛解决的就是这件事它允许模型维度——也就是层数——在采样过程中自由变化最终输出的不是一个“最优分层结果”而是一组覆盖不同层数和不同电阻率结构的后验样本。这套方法特别适合 MT 这类非线性强、等值性严重、观测误差往往不可忽略的反演问题。本文按我实际做一维 MT 反演的流程从正演、采样器设计、参数配置到翻车点把这条技术路线完整拆开讲。2. 从固定层数到可变维度rjMCMC 的接受率到底改了什么2.1 为什么固定层数反演会把不确定性藏起来常规的一维 MT 反演无论是最小二乘还是带正则化的 Occam 型反演第一步都是给定层数。常见做法是手工设 3 到 5 层或者用 Occam 从粗糙模型开始逐步细化。问题在于层数一旦固定反演结果的不确定性就只剩下“每层电阻率的扰动范围”层本身的不确定性被完全忽略。真实地下结构在 MT 数据面前往往是等值性的一层薄的高阻层和一层稍厚的中阻层在某个频段上可以给出几乎一样的响应。固定层数时反演器会选一个局部最优解交给你然后你就拿着这个“唯一”结果去做地质解释——但数据本身根本不足以支撑这个唯一性。rjMCMC 的出发点是把反演放到贝叶斯框架里。后验分布写为 p(m|d) ∝ p(d|m)p(m)其中模型 m 不再是一个长度固定的向量而是长度也可变的结构。在变维度模型空间里马尔科夫链的转移核必须满足跨维度的细致平衡条件从层数为 k 的模型跳到层数为 k1 的模型和反向跳回来的概率流必须平衡。这就是可逆跳跃的核心约束也是它和常规 MCMC 的本质差异。2.2 可逆跳跃的核心不同维度之间的接受率怎么守恒常规 Metropolis-Hastings 接受率是α(x, y) min(1, π(y)q(x|y) / (π(x)q(y|x)))这里 x 和 y 在同维度空间里。rjMCMC 把模型空间扩展成不同维度的并集x 和 y 维数不同q(x|y) 和 q(y|x) 的测度基础不一样直接写这个公式没有意义。可逆跳跃的常见做法是引入一个辅助变量 u让新模型的参数由 y g(x, u) 映射出来然后在接受率里补上雅可比行列式α(x, y) min(1, [π(y)q(x|y)] / [π(x)q(y|x)] × |∂g/∂(x,u)|)做过实际实现的人会知道真实落地时很少去显式算雅可比。常见的设计方式是让出生步骤的新参数从先验分布里直接抽取死亡步骤把被删掉的参数直接丢弃这样映射关系里辅助变量和模型参数的对应是单位映射雅可比行列式为 1。我在自己的代码里就是这么做的接受率退化为log α log π(y) − log π(x) log q(x|y) − log q(y|x)其中 q 的比值只包含“选择哪个节点出生/死亡”的组合项以及新参数提议密度与先验密度相消后的剩余。一个重要的实现细节出生时新电阻率必须从先验分布采样而不是从某个以当前模型为中心的窄分布采样。如果从窄分布采样提议密度不会和先验密度抵消接受率公式里就要显式带上 q_prop(ρ_new) 这一项代码里很容易漏。2.3 先验怎么选才不“拉偏”结果一维 MT 的 rjMCMC 实现里先验一般分三块层数先验、电阻率先验、深度网格先验。层数先验我通常取均匀分布U(1, Nmax)。Nmax 设置多少很关键常见做法是取 8 到 12 层。太小会把可能的构造截断太大会让出生步的搜索空间过大、接受率变低。电阻率先验一般在对数域写这是 MT 反演里最基础的习惯电阻率跨好几个数量级线性域采样效率极低。我常用 log10ρ ~ U(−1, 4)对应 0.1 到 10000 Ω·m。深度网格先验则是把层边界的位置放在一个固定网格上网格在对数深度上均匀分布。深度网格范围需要根据最低频估计。最大探测深度可以用趋肤深度公式粗估δ ≈ 503 × sqrt(ρ_avg / f_min) 米。比如平均电阻率约 1000 Ω·m最低频 0.01 Hzδ 约 500 公里量级但一维 MT 实际探测深度远达不到这个理论值我一般取公式值的三分之一到一半作为最大深度再把网格节点按对数均匀分布放进去。网格放好之后每个节点的“是否激活”作为模型的维数指示变量这就是下一章要讲的节点模型。3. 从正演到采样器一维 MT 的 rjMCMC 反演核心实现步骤3.1 一维 MT 正演从底层半空间往上递推阻抗一维 MT 正演本身不复杂复杂的是把阻抗递推、层参数化和维度跳跃组合到一起。先看正演核心。模型由若干层组成每层有电阻率 ρ_k 和厚度 h_k底层是半空间。从最底层开始递推公式是反射系数法先算底层阻抗再依次向上更新每层顶面的阻抗。import numpy as np from scipy.constants import mu_0 def mt_1d_forward(rho, h, freq): rho: 每层电阻率, 单位 Ohm.m, 最后一层为半空间 h: 每层厚度, 单位 m, 长度应为 len(rho)-1 freq: 频率数组, 单位 Hz 返回视电阻率 rho_a 和相位 phase_rad omega 2.0 * np.pi * freq nlayer len(rho) # 底层半空间阻抗 Z sqrt(i*omega*mu_0*rho[-1]) Z np.sqrt(1j * omega * mu_0 * rho[-1]) # 自底向上递推 for k in range(nlayer - 2, -1, -1): k_wave np.sqrt(1j * omega * mu_0 / rho[k]) # 第 k 层波数 # 反射系数下层阻抗与当前层特征阻抗的差异 gamma (Z - 1j * omega * mu_0 / k_wave) / (Z 1j * omega * mu_0 / k_wave) # 当前层顶面的阻抗含层厚相位延迟 Z (1j * omega * mu_0 / k_wave) * \ (1.0 - gamma * np.exp(2j * k_wave * h[k])) / \ (1.0 gamma * np.exp(2j * k_wave * h[k])) rho_a np.abs(Z) ** 2 / (omega * mu_0) phase np.arctan2(Z.imag, Z.real) return rho_a, phase这段代码里两个细节值得说明。第一k_wave是复波数因为电阻率对应的是有损耗介质gamma的相位会随频率变化np.exp(2j * k_wave * h[k])里的双程相位延迟不能省。第二底层半空间的阻抗sqrt(1j * omega * mu_0 * rho[-1])在低频极限下它的辐角接近 45 度这直接决定了最低频的相位响应所以最低频点的相位数据质量对深部约束很关键。反演时正演输出要和观测数据比对。观测数据有两种组织方式一种是视电阻率加相位的形式另一种是复数阻抗形式。我建议在似然函数里直接用复数阻抗或者对数视电阻率加相位分开建模避免视电阻率跨数量级导致拟合失衡。3.2 节点模型与出生/死亡/扰动三种提议的接受率实现维度跳跃最干净的方式是节点模型。先在深度网格上布置 Nmax 个节点每个节点有一个“激活”标记和一个电阻率值。激活的节点把模型分成若干层层电阻率取该区间内节点的值。这样出生步就是激活一个未激活节点死亡步就是取消一个激活节点扰动步就是修改某个激活节点的电阻率值。厚度网格完全不用动雅可比问题也自然绕开。rng np.random.default_rng() def propose_birth(state, grid_n_active, prior_range): 从当前未激活节点中随机选一个电阻率从先验分布采样 返回候选状态和提议比调整项 log_q_ratio inactive_idx np.where(~state.active)[0] if len(inactive_idx) 0: return state, -np.inf node rng.choice(inactive_idx) rho_new rng.uniform(prior_range[0], prior_range[1]) # 对数电阻率 state_new state.copy() state_new.active[node] True state_new.rho_log[node] rho_new # 出生从 (Nmax-k) 个未激活节点选 1 个死亡从 k 个激活节点选 1 个 k_old state.active.sum() log_q_ratio np.log(k_old) - np.log(Nmax - k_old) return state_new, log_q_ratio这里log_q_ratio就是上一章说的组合项。它只和当前层数有关层数越小出生相对越容易被接受层数越大出生组合项越负防止链无限增长。注意出生时的rho_new是从先验均匀分布采样的所以它的先验密度和提议密度在验前后验比中正好抵消log_q_ratio里不需要额外加log prob(rho_new)。死亡步是出生的逆def propose_death(state): 随机删除一个激活节点的电阻率值恢复为未激活 active_idx np.where(state.active)[0] if len(active_idx) 1: return state, -np.inf node rng.choice(active_idx) state_new state.copy() state_new.active[node] False k_old state.active.sum() log_q_ratio np.log(Nmax - k_old 1) - np.log(k_old) return state_new, log_q_ratio扰动步就简单了在某个激活节点的对数电阻率上加一个高斯扰动def propose_perturb(state, sigma_rho0.3): active_idx np.where(state.active)[0] if len(active_idx) 0: return state, -np.inf node rng.choice(active_idx) state_new state.copy() state_new.rho_log[node] rng.normal(0, sigma_rho) log_q_ratio 0.0 # 对称随机游走提议比抵消 return state_new, log_q_ratio扰动步的提议是随机游走正反向提议密度相同log_q_ratio为零。sigma_rho这个参数值得单独调太小会接受率高但移动慢自相关巨大太大会频繁拒绝链卡在原处。我通常先跑短链看接受率目标扰动接受率维持在 20% 到 40%按需调整sigma_rho。3.3 采样器主循环与并行回火温度链主循环要做的就是把三类提议按固定概率混合计算后验比决定接受或拒绝。核心框架如下def mcmc_step(state, data, forward_func, prior): # 按概率选择提议类型 kind rng.choice([perturb, birth, death], p[0.4, 0.3, 0.3]) if kind perturb: cand, log_q_ratio propose_perturb(state) elif kind birth: cand, log_q_ratio propose_birth(state) else: cand, log_q_ratio propose_death(state) if np.isinf(log_q_ratio): return state, False # 对数似然差 对数先验差 提议比 log_like_diff log_likelihood(cand, data) - log_likelihood(state, data) log_prior_diff log_prior(cand) - log_prior(state) log_alpha log_like_diff log_prior_diff log_q_ratio if np.log(rng.uniform()) log_alpha: return cand, True return state, False这里log_prior要包含两部分层数先验 log p(k) 和每层电阻率先验 log p(ρ)。因为出生步的 ρ 从先验抽取ρ 的先验项在log_prior_diff和log_q_ratio里天然匹配实际影响接受率的先验项主要是层数先验。如果把层数先验取均匀分布这一项更简单。但只有这个主循环的一维 MT rjMCMC收敛往往很慢。非线性反演的后验分布常有多个局部峰链容易陷在某个峰里出不来。常见做法是并行回火同时跑多条链每条链的接受率里把对数似然乘以温度系数 β1/T温度从 1 到 Tmax 几何分布冷链跑正式采样热链在平滑后的似然面上探索相邻温度链之间周期性交换模型。交换接受率的实现也很常规def temper_exchange(state_i, state_j, beta_i, beta_j): # 只交换两个模型状态 log_E_i -log_likelihood(state_i) - log_prior(state_i) log_E_j -log_likelihood(state_j) - log_prior(state_j) log_alpha (beta_j - beta_i) * (log_E_i - log_E_j) if np.log(rng.uniform()) log_alpha: return state_j, state_i # 交换成功 return state_i, state_j温度链数我一般用 6 到 8 条Tmax 取 3 到 5。温度太高热链完全随机游走和冷链的交换接受率会掉到接近零交换机制形同虚设温度太低热链的探索范围不够起不到回火作用。判断温度设置是否合适的标准只有一个相邻链的交换接受率应该在 5% 到 30% 之间。4. 跑通最小反演样例模型参数、先验范围与收敛判断4.1 最小反演样例的模型设置与数据准备假设手里有一套某台站的实测数据频率从 320 Hz 到 0.0625 Hz共约 20 个频点包含视电阻率和相位。实测数据要先做预处理人工剔除明显飞点估计每个频点的观测误差。观测误差的估计直接影响后验宽度不要用固定百分比一刀切。低频段的 MT 数据往往受噪声影响大误差应该按频点单独估计。我一般把数据整理成一个文本文件每行四个数频率、视电阻率、视电阻率误差、相位误差。然后按下面的方式启动反演python run_rjmcmc_1d.py \ --data obs_data.txt \ --freqs freqs.txt \ --zmin 20 \ --zmax 50000 \ --nmax 10 \ --prior-rho-min 0.1 \ --prior-rho-max 10000 \ --iter 300000 \ --burn 0.4 \ --thin 10 \ --nchain 6 \ --tmax 4这些参数是我实际跑一维 MT 反演的常用起点。zmin取 20 米是因为浅部高频点的趋肤深度大约在这个量级zmax取 50000 米对应最低频 0.0625 Hz 且平均电阻率约 1000 Ω·m 的探测深度再深就是“看着像深部实际靠先验撑起来”的区域了。nmax10在一维 MT 反演里是一个平衡点既能表达复杂地层又不至于让出生步的搜索空间爆炸。4.2 迭代次数、燃烧期与保存间隔怎么定rjMCMC 的收敛速度比固定维 MCMC 慢迭代次数不能省。我的习惯是先跑一个 5 万步的短链用来看层数直方图是否合理、出生/死亡接受率是否太低、扰动是否卡壳然后再决定要不要把迭代加到 30 万或更多。燃烧期取 20% 到 40% 是常见做法但更好的方式是直接看轨迹对数后验的轨迹在某个区间附近稳定波动不再有持续上升或下降趋势就把前面那段弃掉。保存间隔thinning的作用是降低自相关。rjMCMC 的出生/死亡步本身就会带来混合动力自相关通常不像纯随机游走那样可怕thin10通常够了。有一点要注意如果保存间隔设置后样本量太少后验分布的分位数会不稳。一个判断标准是最终保存的有效样本数要能画出平滑的层数直方图我一般要求每层至少有 200 到 500 个有效样本。4.3 收敛判断轨迹、R-hat 与层数直方图判断收敛不能只看负对数后验的轨迹还要看具体统计量。rjMCMC 里每步的模型维度不同不能直接把所有参数拼在一起算 R-hat。正确做法是分层统计把“层数”这个标量当作统计量对每条链的层数序列算 R-hat再把固定深度点上的电阻率值取出计算 R-hat。def rhat_from_chains(chains): chains: 形状 (n_chain, n_sample) 的统计量序列 例如每个深度点上的电阻率后验样本 n_chain, n_sample chains.shape between_var n_sample * np.var(chains.mean(axis1), ddof1) within_var np.mean(np.var(chains, axis1, ddof1)) var_hat (n_sample - 1) / n_sample * within_var between_var / n_sample return np.sqrt(var_hat / within_var)R-hat 小于 1.05 是我能接受的底线小于 1.1 只能算“勉强能用”。另一个更直观的判断方法是层数直方图分段对比把后验样本按迭代顺序分成前后两半分别画层数直方图如果两个直方图的峰位置和形状基本一致说明链在维度空间里已经充分混合层数分布稳定下来了。这个方法虽然简单但比任何统计量都更能暴露问题——如果后一半的层数集中在 Nmax 附近而前一半集中在中间层数说明链根本没有遍历。5. rjMCMC 一维 MT 反演避坑清单5 种常见翻车场景5.1 现象层数直方图在 Nmax 处堆成一根柱子链跳不出来这是一个很典型的翻车。现象是后验层数分布的最高峰顶着 Nmax结果解释成“地下真有十层”实际上链卡在了高层数状态里出不来。原因有两个一是 Nmax 设得过大出生步从先验抽取的电阻率和当前模型完全不搭接受率极低二是并行回火温度不够链在局部峰之间无法跳跃。解决方法是先压 Nmax 到 6 或 8同时提高 Tmax并检查出生接受率。如果出生接受率只有 1% 到 2%说明出生步的提议分布和数据的支撑区间严重不匹配。这个时候不要急着调接受率公式而是先看先验范围是否合理——先验电阻率范围太宽会让出生的新层大概率落在似然很低的位置。5.2 现象结果对先验范围极度敏感换一个边界结果变一个样现象是电阻率先验从 [0.1, 10000] 改成 [1, 10000] 之后深部后验的结构肉眼可见地变了。这不一定是你实现错了。一维 MT 数据对深部电阻率的分辨能力本来就弱等值性使得“薄低阻层厚高阻层”和“厚中阻层”可以对应几乎相同的响应这时深部信息的约束力很差后验分布主要由先验在撑。这类问题没有算法层面的根治办法。我能给的建议是把先验敏感性分析当作流程的一部分至少先验上下边界各跑一次短链对比结果差异。如果差异很大在写解释报告时就要如实说明“深部结构主要受先验约束数据本身不能区分”。强行上固定层数反演会让结果看起来更确定性但那是在给不确定性“化妆”。5.3 现象误差标准差设得不对后验宽度完全失真观测误差标准差是似然函数里的关键参数。把误差设得过小后验会非常尖锐链在似然面上像走钢丝一样难受混合极差误差设得过大后验会被摊平地层的边界全部模糊掉看起来“什么也分不出来”。一个隐蔽的坑是相位和视电阻率的误差单位不一致相位误差 1 度和视电阻率误差 5% 在数值上完全不是一个量级如果直接相加反演会默认视电阻率的比值权重更大相位几乎不参与拟合。解决方法是把两类数据统一到同一个相对误差框架或者在似然函数里分别归一化。我通常的做法是视电阻率用对数残差建模相位用弧度残差建模各自除以各自的误差标准差再相加。跑完后检查每个频点的标准化残差如果某几个频点系统性偏正或偏负说明数据里还有没剔除的飞点或误差估计偏了。5.4 现象并行回火温度太高交换接受率直接归零现象是日志里相邻温度链的交换接受率常年低于 1%冷链和热链各自独立运行温度链形同虚设。原因通常是 Tmax 取得太大冷链和热链的能量分布重叠区域太小交换概率指数级衰减。解决方法是先看交换接受率如果低于 5%把 Tmax 降下来或者增加中间温度链的条数。温度链的间隔在几何尺度上均匀分布也可以用自适应调温算法让链自己去匹配目标接受率。注意一个细节交换操作应该在温度链之间周期性执行而不是每步都做否则交换本身会引入额外自相关。我一般每 50 步做一次温度交换。5.5 现象低频段残差超大链反复在深部电阻率上大幅抖动低频点是最容易出问题的。现象是正演预测在低频段系统性偏离观测值且后验深部电阻率的分布宽到没有解释价值。常见原因有两个一是低频观测数据本身的噪声被低估二是深度网格最大深度不够模型表达不了真实深部结构低频响应只能靠浅层结构硬凑。排查顺序是先看低频点的误差估计是否合理MT 低频段经常受自然场源信号弱、噪声占比高的影响误差放大 2 到 3 倍是常有的事再看正演响应里最低频点的相位是否在 45 度基线区如果相位明显偏离且数据质量很好那大概率是模型深度截断了。深部截断问题可以改 zmax 加一个数量级重跑对比但要注意计算量也会跟着涨。6. 结果不要只看均值后验概率密度、MAP 与预测检查的四步读法6.1 层数直方图和 MAP 模型先看层数直方图确认后验主峰不是被 Nmax 截断的。然后从后验样本里找出对数后验最大的那组模型这就是 MAP 模型它比后验均值模型更接近真实数据拟合的峰值。做地质解释时我习惯把 MAP 模型的电阻率分层结构和层数直方图放在同一张图上一个代表“数据最支持的模型”一个代表“数据对层数的分辨能力”。6.2 后验概率密度剖面与不确定性带状图把样本按深度-电阻率二维直方图统计就能画出一张后验概率密度剖面。这张图的横轴是对数电阻率纵轴是深度颜色代表后验概率密度。这张图最大的价值是能直观看到等值性如果某个深度区间出现两个分离的高密度峰而不是一条细窄的带说明数据在那个区间只能给出“两个可能电阻率值”这不是反演器坏了而是数据本身的信息极限。不确定性的呈现我一般用中位数剖面配 95% 可信区间带状图。这里有个实践教训不要把后验均值模型直接当成解释模型因为均值会把多峰分布抹平成过渡带看起来像有一个缓慢变化的梯度层实际可能来自两个截然不同的构造模型。6.3 后验预测检查用样本重新正演看预测包络是否盖住观测从后验样本里取出若干模型重新算正演响应把所有预测曲线的分布画出来。如果观测点落在预测包络之外而且不是一两个频点说明似然函数或误差估计有问题。这一步的成本很低却是整个流程里最有效的验证手段比任何收敛诊断都直观。我会把预测包络图直接附在处理报告中让读者一眼看出反演结果到底是在拟合数据还是只是在拟合先验。最后说一个我自己的习惯每次跑 rjMCMC先只上 5 万步短链只画两张图——层数直方图和温度的交换接受率。层数直方图能看到维度混合情况交换接受率能看到回火是否有效。这两张图过了关再加大迭代量到正式发布用的精度。这个习惯帮我避掉了至少两次大规模返工省下的时间足够重跑十几轮完整反演。希望帮到你。本文还有配套的精品资源点击获取