基于rjMCMC的一维大地电磁反演：原理、Python实现与避坑指南

简介这份资源面向地球物理反演方向的研究生、科研人员及工程师提供可逆跳跃马尔科夫链蒙特卡洛rjMCMC方法在一维大地电磁反演中的完整实现。它解决传统反演难以处理模型维数变化与多解性的问题通过在不同层数与参数分布间自动跳转获得电性结构的后验分布适合具备一定贝叶斯推断与MATLAB基础的中高级学习者。压缩包共24个文件约7.59MB以14个m脚本和8个mat数据文件为主辅以README说明与LICENSE脚本覆盖正演、似然函数、扰动、收敛指标与结果绘图等模块mat文件保存多组MCMC采样链与模型数据。目前已有253人学习下载。读者可借此理解rjMCMC的完整流程包括先验构建、迭代采样与后验分析并直接运行脚本复现一维反演实验为后续拓展至二维或实际数据应用提供可参考的代码框架与排错思路。1. 从一条测深曲线说起为什么 rjMCMC 值得在地磁反演里占一个位置一维大地电磁MT反演说白了就是拿一条视电阻率和相位随周期变化的曲线去猜地下几十公里内每层的电阻率和厚度。传统做法是固定层数比如三层、四层然后拿最小二乘或 Occam 去拟合。问题在于层数这个事你事先根本不知道。定少了模型欠拟合曲线尾巴压不住定多了参数冗余反演出来的层要么薄得没意义要么电阻率在几个数量级之间乱跳这就是一线常说的“玄学反演”。可逆跳跃马尔科夫链蒙特卡洛rjMCMC也叫 Reversible Jump MCMC解决的就是这个痛点让层数本身变成一个待采样的随机变量。链在“加一层”“删一层”“改电阻率”“改厚度”这几种跳跃之间游走最后你拿到的不是一个最优模型而是一整个模型集合的后验分布。对 MT 这种欠定问题来说后验分布比单点解诚实得多——它直接告诉你哪些深度是数据约束得住的哪些深度是模型自己编的。这个方向适合两类人一类是做深部电性结构、地热或资源勘查需要给出不确定性量化的另一类是已经写过固定层数 MT 反演、想把手里的代码升级成贝叶斯框架的。下面按“先立住原理、再动手复现、最后讲坑”的顺序拆开讲代码用 Python 写核心逻辑不依赖任何商业软件。2. rjMCMC 反演一维 MT 的数学骨架与先验怎么定2.1 从贝叶斯公式到跨维后验一维 MT 的正演常见做法是把地下切成若干水平层每层给电阻率 ρ 和厚度 h最底层半空间只给电阻率。给定模型 m用递推公式算出地表阻抗 Z再转成视电阻率 ρa 和相位 φ。观测数据 d 就是各周期的 ρa 和 φ误差一般假设对数域独立高斯。贝叶斯公式写出来是 p(m|d) ∝ p(d|m) p(m)。固定层数时m 的维度固定MCMC 直接采就行。但层数 k 是变量不同 k 对应不同维度的参数空间普通 Metropolis-Hastings 没法在维度之间跳。rjMCMC 的做法是构造一个跨维的提议分布让链能在 k 和 k1 之间跳同时用一个 Jacobian 项修正维度变化带来的测度改变。接受率公式里那个 Jacobian是整套方法最容易写错的地方后面避坑章会专门讲。似然函数取对数后是import numpy as np def log_likelihood(model, periods, rhoa_obs, phase_obs, rhoa_err, phase_err): model: dict, 含 rho (list, 长度 k), thk (list, 长度 k-1) periods: 周期数组, 秒 rhoa_obs, phase_obs: 观测视电阻率(Ohm·m)与相位(度) rhoa_err, phase_err: 对应的对数域/线性域误差 返回对数似然 rhoa_pred, phase_pred forward_1d(model, periods) # 视电阻率在对数域比较相位在线性域比较 res_rhoa (np.log(rhoa_obs) - np.log(rhoa_pred)) / rhoa_err res_phase (phase_obs - phase_pred) / phase_err return -0.5 * np.sum(res_rhoa**2 res_phase**2)这段的关键点视电阻率跨好几个数量级必须在对数域算残差相位本身是线性量直接减。误差项 rhoa_err 和 phase_err 不是随便填的一般取观测误差的估计值如果数据没给误差常用 0.05对数域和 1.0 度作为默认但这会直接影响后验宽度属于必须交代清楚的参数。2.2 先验分布别让先验替数据说话rjMCMC 里先验分两层一层是层数 k 的先验一层是给定 k 后各参数的先验。层数先验常用泊松分布或均匀分布截断到某个上限比如 k ∈ [1, 15]。我一般用截断泊松λ 取 3 到 5这样链不会一上来就往十几层跑。给定 k 后电阻率通常取对数均匀先验范围比如 0.1 到 10000 Ohm·m厚度取均匀先验范围 10 m 到 100 km且要保证总厚度不超过某个物理上限。先验不是装饰它直接决定后验的边界。如果先验范围设窄了后验会被截断看起来“约束很好”其实是先验在说话。def log_prior(model, k_max15, rho_min0.1, rho_max1e4, h_min10.0, h_max1e5, lam4.0): k len(model[rho]) if k 1 or k k_max: return -np.inf # 层数先验截断泊松 lp k * np.log(lam) - lam - np.math.lgamma(k 1) # 电阻率对数均匀先验 for r in model[rho]: if r rho_min or r rho_max: return -np.inf lp - np.log(np.log(rho_max / rho_min)) # 厚度均匀先验 for h in model[thk]: if h h_min or h h_max: return -np.inf lp - np.log(h_max - h_min) return lp参数说明lam 是泊松均值控制层数偏好rho_min/rho_max 和 h_min/h_max 是先验硬边界超出直接返回负无穷链会被拒绝。这几个值要根据工区实际情况调比如沉积盆地浅层薄层多h_min 可以放到 5 m深部探测 h_max 可以到 200 km。2.3 跨维跳跃的四种提议与 JacobianrjMCMC 的核心是设计跳跃类型。一维 MT 里常用四种更新电阻率、更新厚度、加一层、删一层。加层时从当前模型里选一层把它劈成两层新层的电阻率和厚度由随机数生成删层时反过来把相邻两层合并。接受率里必须包含 Jacobian因为加层是从低维到高维参数空间体积变了。常见做法是加层时用“分裂-合并”配对分裂时新参数由旧参数加扰动生成合并时用逆变换。这样 Jacobian 可以解析算出来。如果偷懒不做配对接受率会系统性偏低链几乎不跳最后退化成固定层数反演。这是血泪经验后面还会展开。3. 用 Python 把 rjMCMC 一维 MT 反演跑起来3.1 正演函数递推阻抗与视电阻率相位正演是反演的地基写错了后面全白搭。一维 MT 用递推公式从最底层往上算def forward_1d(model, periods): 一维 MT 正演返回视电阻率(Ohm·m)和相位(度) model: {rho: [...], thk: [...]} periods: 周期数组, 秒 mu0 4 * np.pi * 1e-7 rho np.array(model[rho], dtypefloat) thk np.array(model[thk], dtypefloat) k len(rho) omega 2 * np.pi / periods rhoa np.zeros_like(periods, dtypefloat) phase np.zeros_like(periods, dtypefloat) for i, w in enumerate(omega): # 最底层半空间的本征阻抗 Z np.sqrt(1j * w * mu0 * rho[-1]) # 从倒数第二层往上递推 for j in range(k - 2, -1, -1): kj np.sqrt(1j * w * mu0 / rho[j]) Zj np.sqrt(1j * w * mu0 * rho[j]) # 层内传播因子 ej np.exp(-2 * kj * thk[j]) R (Zj - Z) / (Zj Z) Z Zj * (1 - R * ej) / (1 R * ej) rhoa[i] np.abs(Z)**2 / (w * mu0) phase[i] np.degrees(np.arctan2(Z.imag, Z.real)) return rhoa, phase逻辑说明从最底层半空间阻抗出发逐层向上用反射系数和传播因子递推。thk 长度是 k-1因为最底层没有厚度。参数单位要统一rho 用 Ohm·mthk 用米periods 用秒。相位返回的是度和观测数据对齐。这个正演对层数不敏感k 变了直接改数组长度就行正好配合 rjMCMC。3.2 四种跳跃的提议与接受率实现下面把加层和删层的核心逻辑写出来更新电阻率和厚度就是普通 Metropolis 提议这里不展开。def propose_add_layer(model, rng): 在当前模型里随机选一层劈成两层返回新模型和 Jacobian 对数 k len(model[rho]) idx rng.integers(0, k) rho_old model[rho][idx] # 新电阻率在对数域加扰动 log_rho_new np.log(rho_old) rng.normal(0, 0.5) rho_new np.exp(log_rho_new) # 厚度分裂旧层厚度按比例分 if idx k - 1: h_old model[thk][idx] u rng.uniform(0.2, 0.8) h1, h2 h_old * u, h_old * (1 - u) new_rho model[rho][:idx] [rho_new, rho_old] model[rho][idx1:] new_thk model[thk][:idx] [h1, h2] model[thk][idx1:] else: # 最底层半空间加层需要给它一个厚度 h_new rng.uniform(10, 5000) new_rho model[rho][:idx] [rho_new, rho_old] new_thk model[thk] [h_new] # Jacobian对数变换 厚度分裂的比例项 log_jac np.log(rho_new) - np.log(rho_old) np.log(h_old) if idx k-1 else np.log(rho_new) - np.log(rho_old) return {rho: new_rho, thk: new_thk}, log_jac参数说明扰动标准差 0.5 是对数域电阻率太大接受率低太小跳不动一般 0.3 到 0.8 之间试。厚度分裂比例 u 限制在 0.2 到 0.8避免产生极薄层导致正演数值问题。Jacobian 这里只写了主要项实际实现里还要乘上提议分布的比值漏掉任何一项接受率都会偏。删层是加层的逆操作选相邻两层合并新电阻率取对数平均新厚度相加Jacobian 取倒数。加层和删层的接受率要成对出现否则细致平衡被破坏后验会有偏。3.3 主循环与收敛判断主循环就是不断提议、算接受率、按概率接受。关键参数总迭代次数、燃烧期、 thinning 间隔。一维 MT 数据量不大一般跑 20 万到 50 万次燃烧前 5 万次之后每 50 步存一个样本。def run_rjmcmc(d, periods, n_iter200000, burn50000, thin50, seed42): rng np.random.default_rng(seed) # 初始模型3 层 model {rho: [100.0, 1000.0, 100.0], thk: [500.0, 5000.0]} samples [] n_accept {add: 0, del: 0, update: 0} for it in range(n_iter): move rng.choice([add, del, update], p[0.3, 0.3, 0.4]) if move add: new_model, log_jac propose_add_layer(model, rng) log_alpha (log_likelihood(new_model, *d) log_prior(new_model) - log_likelihood(model, *d) - log_prior(model) log_jac) if np.log(rng.uniform()) log_alpha: model new_model n_accept[add] 1 elif move del: # 删层逻辑与加层对称此处省略具体实现 pass else: # 更新电阻率或厚度普通 MH pass if it burn and it % thin 0: samples.append(model.copy()) return samples, n_accept逻辑说明move 的概率分配影响混合效率加层和删层各 0.3、更新 0.4 是常用起点。接受率要监控加层接受率低于 5% 说明提议步长太大或 Jacobian 有问题高于 50% 说明步长太小链探索不充分。收敛判断不能只看一条链至少跑 3 到 4 条不同初值的链看层数后验分布是否一致。4. 反演结果怎么读后验分布、层数直方图与不确定性4.1 层数后验直方图告诉你数据支持几层跑完之后第一件事是画层数 k 的直方图。如果数据信噪比高、频带宽直方图会在某个 k 附近出尖峰如果数据差直方图会很平说明层数根本约束不住。这时候你报“最优层数”就是自欺欺人应该报整个后验。常见做法是把层数后验和电阻率-深度剖面叠在一起看层数概率高的地方剖面收得紧层数概率低的地方剖面散得开。这比单点反演的信息量大得多。4.2 电阻率-深度后验剖面的画法后验剖面不是简单把样本叠起来因为每个样本的层边界不一样。标准做法是在深度轴上打网格对每个深度统计所有样本在该深度的电阻率分布取分位数画阴影。def posterior_profile(samples, depth_grid): 把变层数样本插值到统一深度网格返回分位数 rho_at_depth np.zeros((len(samples), len(depth_grid))) for i, m in enumerate(samples): # 把层模型转成深度-电阻率阶梯函数 depths np.concatenate([[0], np.cumsum(m[thk])]) for j, z in enumerate(depth_grid): idx np.searchsorted(depths, z, sideright) - 1 idx min(idx, len(m[rho]) - 1) rho_at_depth[i, j] m[rho][idx] return (np.percentile(rho_at_depth, 5, axis0), np.percentile(rho_at_depth, 50, axis0), np.percentile(rho_at_depth, 95, axis0))参数说明depth_grid 一般取对数等间隔从 10 m 到 100 km。分位数取 5%、50%、95% 是惯例也可以取 16% 和 84% 对应一倍标准差。注意这个剖面在层边界附近会有平滑效应因为不同样本的边界位置不同插值后过渡带变宽这是正常的不是反演失败。4.3 用合成数据验证先跑通再上实测上实测数据之前必须用合成数据验证整套代码。流程是给定一个已知层模型正演生成视电阻率和相位加高斯噪声然后跑 rjMCMC看后验是否覆盖真值。# 合成数据验证 true_model {rho: [50.0, 500.0, 50.0], thk: [300.0, 3000.0]} periods np.logspace(-2, 3, 30) rhoa_true, phase_true forward_1d(true_model, periods) rng np.random.default_rng(0) rhoa_obs np.exp(np.log(rhoa_true) rng.normal(0, 0.05, len(periods))) phase_obs phase_true rng.normal(0, 1.0, len(periods)) d (periods, rhoa_obs, phase_obs, 0.05, 1.0) samples, acc run_rjmcmc(d, periods)如果后验 90% 区间盖不住真值先查正演对不对再查似然函数的误差项有没有写反。合成数据跑通之前不要碰实测数据否则你分不清是代码错还是数据差。5. 避坑与排查rjMCMC 一维 MT 反演最容易翻车的五个地方5.1 加层接受率极低链几乎不跳现象跑完发现层数始终停在初始值附近加层接受率不到 1%。原因通常是 Jacobian 漏项或提议步长过大。解决先检查加层和删层的 Jacobian 是否互为倒数再逐步调小电阻率扰动标准差从 0.5 降到 0.2 试。如果还不行检查先验边界是不是把新层卡死了比如新生成的厚度超出 h_max。5.2 后验层数分布过宽看起来“什么层数都可能”现象层数直方图从 2 层到 10 层都有概率没有明显峰值。原因可能是数据频带太窄或误差估计过大。解决先确认观测误差是不是被高估了如果误差项填得比实际大似然面会变平后验自然散。另外检查周期范围如果只有高频没有低频深部信息本来就少层数约束不住是正常的这时候应该老实报后验而不是硬选一个层数。5.3 电阻率后验剖面在深部发散现象浅部剖面收得很紧深部 5% 和 95% 分位数差好几个数量级。原因是一维 MT 对深部低阻层不敏感这是物理限制不是代码问题。解决在解释时明确标出分辨深度一般用趋肤深度公式估算超过这个深度的后验只能作为定性参考。如果工区确实需要深部约束得加先验信息或联合其他数据。5.4 链的自相关性太强有效样本数不够现象跑了 20 万次算有效样本数只有几百。原因是 thinning 间隔太小或提议步长太小。解决先算自相关函数看滞后多少步降到 0.1 以下thinning 间隔至少取这个滞后步数的两倍。另外可以并行跑多条链每条链不同初值最后合并这样比单链跑很久更有效。5.5 正演数值溢出导致似然变成 NaN现象链跑着跑着突然报 NaN或者接受率骤降。原因是某些模型参数组合下正演里的指数项溢出。解决在正演函数里加保护比如对传播因子 ej 做截断超过 1e10 就当 0 处理或者在似然函数里检查 NaN一旦出现直接返回负无穷让链拒绝这个模型。另外厚度不能取 0先验里 h_min 要设一个正数。6. 进阶技巧用并行链和分层提议把混合效率提上去跑通基础版之后最影响体验的就是混合效率。单链跑 50 万次可能才勉强收敛实际项目里等不起。我一般会做两件事并行多链和分层提议。并行多链不是简单重复跑而是每条链用不同初值和不同随机种子跑完之后用 Gelman-Rubin 统计量判断收敛。如果 R 值接近 1说明链之间一致可以合并如果 R 值大说明还没收敛继续跑。Python 里用 multiprocessing 就能实现注意每个进程独立随机种子避免伪独立。分层提议是针对加层和删层的优化。基础版里加层是随机选一层劈开但实际数据对浅部和深部的敏感度不一样浅部层数应该比深部更容易跳。做法是给每层一个权重浅层权重高深层权重低提议时按权重选。这样浅部结构能更快被探索深部也不会被无意义的加层拖慢。def weighted_choose_layer(model, rng, depth_grid): 按深度权重选一层浅层权重大 depths np.concatenate([[0], np.cumsum(model[thk])]) mid_depths 0.5 * (depths[:-1] depths[1:]) weights 1.0 / (1.0 mid_depths / 1000.0) # 1 km 以内权重高 weights weights / weights.sum() return rng.choice(len(model[rho]), pweights)参数说明权重函数里的 1000 m 是特征深度可以根据工区调整。这个改动不破坏细致平衡因为权重只影响提议分布接受率里会乘上对应的提议概率比值。改完之后加层接受率通常能翻倍有效样本数也明显上升。最后说个习惯每次改完提议分布或先验我都会先用合成数据跑一遍确认后验覆盖真值再上实测。这个后悔药成本很低但能省掉大量排查时间。希望帮到你。本文还有配套的精品资源点击获取