ERT电阻层析成像MATLAB实现：从正问题到CGLS反演全解析

简介：面向电导率/电阻层析成像（ERT）方向的研究者与工程师，提供一套基于MATLAB的ERT成像仿真实现。资源包含主程序main.m及JacobianERT.m、nodeeit.m等核心算法脚本，配合jacobian矩阵、电压实测数据等mat文件，可完成从边界测量到电导率分布图像重建的完整流程；另有COMSOL仿真操作PDF和运行结果JPG，便于对照验证与二次开发。压缩包共18个文件，以m脚本和mat数据为主，附带pdf、txt等辅助资料，整体仅1.55MB，轻量便捷。内容预览中可见空场模型mphbin和多种测量函数，适合具备一定MATLAB基础、需要快速搭建ERT实验或深入理解电阻层析成像算法的读者。目前已有269人学习使用，可直接用作课程设计、算法测试与科研预研的实用参考。

1. 电导率成像的逆问题：ERT电阻层析成像在matlab里怎么落地

一个装满液体的管道内壁等距排布16个电极，从第1对电极注入电流、第3对电极测量电压，轮换一圈能拿到上百个电压读数。电阻层析成像（ERT）要做的事，就是根据这些边界电压，把内部电导率分布重建出来。这套基于MATLAB的ERT源码包，把正问题有限元离散、Jacobian灵敏度矩阵、CGLS反演和结果绘图全部串在main.m里，还附带了COMSOL空模型、参考电压uref.mat和若干测量数据。我在MATLAB 2019b下直接跑通，再换成自己生成的测量数据后，发现电极策略、迭代次数、网格尺寸对重建结果的影响都能直观看到。适合做两相流检测、电化学过程监测，或者想搞懂电导率反问题如何用matlab代码实现的工程师。下面从正问题和Jacobian矩阵开始拆。

2. 从拉普拉斯方程到JacobianERT.m：ERT正问题与灵敏度矩阵的matlab实现

2.1 电流场控制方程与有限元离散的边界条件

在ERT正问题中，激励频率通常在几十kHz以下，感应电场和位移电流能忽略，电流密度与电场强度满足欧姆定律的微分形式，电荷守恒得到 ∇·(σ∇φ)=0。σ是电导率分布，φ是电位。求解域边界上大部分是绝缘面，法向电流密度为0；只有被选中注入电流的一对电极上有非零电流边界条件。这正是标准椭圆型方程，用电位有限元离散很合适。

以最简单的三角形线性单元为例，单元内电位插值基函数是坐标的线性函数，组装成总体刚度矩阵K后，节点电位满足 Kφ = b。这里的b由电流激励位置决定，不同的激励模式只是换掉b的少数非零项。源码包里的nodeeit.m负责节点自由度的排列，Currenteit.m用来生成电流注入向量。如果对多个激励模式循环求解，K只组装一次，用LU分解后在MATLAB里反复回代即可，避免每次都从头解方程组。

需要注意电极模型对正问题精度的影响。点电极模型把电极缩成一个节点，形成边界上的一点电流源，实现简单但会在电极附近产生奇异的电位梯度。完整电极模型则把电极看作等位体，还会引入接触阻抗，更接近真实测量，但代码里要额外处理电极自由度。源码包内的empty-model.mphbin能在COMSOL中建立同样的几何和电极，把导出的电位分布与MATLAB正问题结果对比，可以很快判断用的是哪种模型。

我做仿真时一般会先直接跑通默认参数，然后把某一对电极的电流方向反转，观察测量电压是否变化。理论上线性正问题对换激励和测量位置具有互易性，如果结果连互易性都不满足，多半是节点排序或电极编号出了问题，而不是反演算法的问题。

2.2 Jacobian矩阵的行列含义与JacobianERT.m核心流程

Jacobian矩阵J的每一行对应一次独立电压测量，每一列对应一个有限元单元的电导率变量。J的元素∂V_i/∂σ_j表示第i个测量电压对第j个单元电导率的偏导数。因为ERT反问题是非线性的，J需要在当前电导率分布处计算，迭代过程中还要不断更新。

计算J的最直接方式是数值差分：对每个单元加一个小扰动δσ，重解正问题，看电压变化量。这个方法代码短但耗时，16电极、800单元时明显卡顿。常见做法是采用伴随法或解析导数，一次正问题加上一次对全部单元的矩阵向量运算就能得到一整列。

源码包中JacobianERT.m的返回结果应是一个大小为测量数×单元数的矩阵。这里给一个等价的核心计算骨架：

function J = JacobianERT(sigma, Node, Edge, Elec) % sigma : 当前电导率分布，列向量，长度为单元总数 % Node : 节点坐标数组 % Edge : 三角形网格的边矩阵，每行是一条边 % Elec : 电极节点编号 % 返回值 J 的大小为 mea_num × elem_num K = assemble_stiffness(Node, Edge, sigma); % 组装总体刚度矩阵 U = solve_all_drive(Node, Edge, K); % 求解所有激励模式下的节点电位 M = measure_operator(Node, Edge, Elec); % 节点电位到测量电压的映射 J = spalloc(mea_num, length(sigma), mea_num*ceil(length(sigma)/8)); for e = 1:length(sigma) dK = element_stiffness_deriv(Node, Edge, e); % 单元e的刚度矩阵对sigma求导 dU = -K \ (dK * U); % 电位对sigma_e的偏导 J(:, e) = M * dU(:); % 换算成测量电压变化 end end

代码里最关键的是K\这一行，它利用同一K矩阵对多个右端项求解，效率远高于循环内重复inv。dU的每一列对应一个激励模式，所以M * dU(:)把所有模式、所有测量对的灵敏度都合到了一列。measure_operator必须和采集电压时的电极顺序一致，否则后边的反演结果会是错的，这在所有ERT程序里都是最容易忽略的一步。

如果不想自己从头写，可以在源码包基础上把JacobianERT.m中电极编号部分改成自己的实验配置。每一轮激励的注入电极编号、测量电极编号都应该以向量形式存在独立的变量里，例如drive = [1,2; 2,3; ...]、measure = [3,4; 4,5; ...]，这样J的行自然按同一顺序排列。很多数据错乱问题，都是因为文本里记录的测量序列和代码中内置的序列不一致。

2.3 电极数、激励模式与Jacobian矩阵规模

在相邻激励模式下，N个电极会得到 N(N-3)/2 个独立测量。原因是从第1电极注入、第2电极流出，测量电极可以从第3、4到N-2、N-1，去掉对称重复；轮换一圈后总数是 N(N-3)/2。下面是几个典型规模：

从表格能直接看出ERT的先天问题：独立测量数远小于网格单元数，J是矮胖矩阵，反问题严重欠定。这也是为什么后面必须靠CGLS迭代和正则化找稳定解。还有一个实用点：J应该用稀疏矩阵存储，否则网络稍大，两次矩阵乘法就能把内存吃满。建议运行后先用size(J)和nnz(J)/prod(size(J))检查一下，前者确认行列数符合预期，后者确认稀疏程度。

3. 从电压差到电导率差：CGLS迭代反演与正则化参数选择

3.1 病态性与正则化的数学取舍

把正问题写成线性化形式 J Δσ = ΔV，其中 ΔV 是测量电压与参考电压之差，Δσ 是电导率增量。由于J的行数远小于列数，加上测量噪声，直接求最小二乘解会产生巨大伪影。标准处理是给目标函数加惩罚项：

min ||JΔσ - ΔV||² + λ² ||L Δσ||²。

L常常是单位矩阵或一阶差分矩阵。λ是正则化参数，过大则图像太平滑，小目标被抹掉；过小则噪声被放大。工程上经常用L曲线法或试算几组λ来选。源码包中的cgls.m走的是另一条路：不显式引入λ，而是通过提前终止迭代来达到类似效果。

CGLS属于Krylov子空间方法，核心运算只有J和J^T的矩阵向量乘。它迭代过程中解从零开始逐渐逼近真解，前期恢复大尺度结构，后期开始拟合噪声。因此，迭代次数越少越稳定，图像越模糊；迭代次数越多越锐利，伪影也越明显。理解这一点后，就不难理解源码包里为什么明明有更直接的最小二乘函数，却仍然用CGLS。

3.2 一个可直接嵌入main.m的CGLS实现

下面这段代码是CGLS的简约版，与源码包内cgls.m的算法思想一致，但变量名更直白，方便改参数。

function [x, iter_used] = cgls_demo(J, b, maxiter, tol) % J: m×n 灵敏度矩阵 % b: m×1 电压差向量 % maxiter: 最大迭代次数，常用值 5~20 % tol: 相对残差阈值，例如 1e-6 x = zeros(size(J,2),1); r = b; % 初始残差 p = J'*r; % 搜索方向 z = r'*r; % 残差平方和 gamma = p'*p; % 方向向量范数 for iter = 1:maxiter q = J*p; alpha = z / (q'*q); % 步长 x = x + alpha*p; r = r - alpha*q; new_z = r'*r; if sqrt(new_z / (b'*b)) < tol iter_used = iter; return; end beta = new_z / z; % 方向更新系数 p = J'*r + beta*p; z = new_z; end iter_used = maxiter; end

参数说明：p是共轭方向，alpha由残差与搜索方向的正交关系确定，beta保证新方向与上一方向关于J'J共轭。tol不能设得太小，否则在16电极ERT场景下很容易迭代到噪声拟合区域。源码包中调用时，一般会把cgls返回的x加到背景电导率上，得到最终分布。

需要说明的是，这段代码没有显式处理非负约束。如果重建出负电导率，可以把负值截断为零，但更好的做法是在外层加约束，或使用带边界约束的变体。源码包的默认场景是电导率相差不大的液体两相流，因此不做约束也能得到可读图像。

3.3 迭代次数与重建效果的对应关系

CGLS的迭代次数就是正则化强度，这个特性对ERT特别有利。下表是基于16电极、104个独立测量、约800个网格单元的典型表现：

操作时建议从5次开始跑，观察残差下降速度：如果前3步残差快速下降，后面几乎不动，那就不需要继续迭代。如果发现重建图的边界出现“光环”状高亮，先降迭代次数到2~3，再判断问题是否出在Jacobian矩阵上。这个顺序能帮你快速区分算法问题和数据问题。

4. main.m数据流复现：文件角色、运行步骤与COMSOL联合调整

4.1 源码包文件角色速览

把源码包展开后，真正需要在MATLAB里运行的是main.m，其他m文件都是它调用的函数。初次接到这套代码，最好先浏览一遍文件结构。

运行前的准备比较固定：把上述文件解压到一个英文路径目录下，例如D:\ert_demo，在MATLAB当前文件夹切到该目录，然后双击main.m。为了防止某个函数不在路径里，可以在命令窗口先执行：

cd 'D:\ert_demo'; addpath(genpath(pwd));

addpath(genpath(pwd))会把当前目录下所有子目录加入搜索路径，避免Undefined function报错。如果文件夹里有中文字符名，某些MATLAB版本在读取uref.mat时会出错，所以我一般会改成纯英文目录。

4.2 main.m内部做的事情：加载数据、计算J、反演、绘图

把main.m的执行逻辑画成数据流就是这样：先加载xy.mat里的网格和坐标，加载uref.mat里的参考电位，接着用JacobianERT.m生成灵敏度矩阵，再从measure1.m得到当前测量电压，计算差值ΔV，用cgls.m迭代得到电导率增量，最后把初始电导率加上增量并作图。

下面是一段简化版的主流程骨架，用于理解参数从哪来、结果到哪去：

% main.m 核心流程（示意） load('xy.mat'); % 节点坐标 load('uref.mat'); % 参考场电位/电压 sigma0 = ones(nElem,1); % 初始电导率分布 drv = [1 2; 2 3; 3 4]; % 电流注入电极对 msr = [3 4; 4 5; 5 6]; % 电压测量电极对 J = JacobianERT(sigma0, Node, Edge, Elec); % 灵敏度矩阵 b = measure1(drv, msr) - uref; % 电压差 dsigma = cgls(J, b, 8, 1e-6); % 8次迭代 sigma_recon = sigma0 + dsigma; % 更新电导率 show_image(xy, sigma_recon); % 显示重建图像

这里drv和msr必须与采集数据时的激励、测量轮换顺序严格一致。如果实际实验是“1-2注入、3-4测量”开始的，那么J第一行对应的就应该是这个组合。measure1函数名来自源码包，实际内容可能是从文本读取，也可能是生成仿真数据，不影响这条数据流。

运行完成后会出现类似“运行结果.jpg”的效果图，通常是两个并排图：左边是设定的真实电导率分布，右边是重建结果。如果只有右图，可以自己写colorbar标注电导率单位。电导率本身的单位是S/m，但仿真中常常用相对值，图像关注的是空间分布和对比度。

4.3 与COMSOL联合仿真时的数据对齐和问题排查

源码包里出现的empty-model.mphbin可以用COMSOL打开。常见做法是在COMSOL中建好几何、电极和网格，导出节点电位或测量电压，再用MATLAB里的text.m读取。但这一步最容易出问题的不是计算，而是单位。

COMSOL默认长度单位是m，而MATLAB网格坐标可能直接来自millimeter或centimeter的网格文件。如果单位不一致，Jacobian矩阵敏感度会整体偏移，导致重建图像形状失真。建议先把COMSOL导出的坐标和xy.mat中的坐标画在同一张图上，确认两者重合再用。

运行中的常见报错和处理方式如下：

最后一行的“重建结果全为背景色”在仿真数据里最经常出现，原因是参考电位uref.mat是用某个特定电导率场算出来的，而measure1得到的测量数据又来自另一个场，两者相减后如果很小，CGLS第一步就认为已经收敛。这时候可以先把ΔV的范数打印出来，看是不是数量级过小。

5. 重建质量验证：残差曲线、网格叠加与J矩阵顺序探针

5.1 用相对残差决定该不该停止迭代

CGLS没有显式目标函数值，最直接的验证指标是相对残差：

V_rec = J * dsigma; rel_res = norm(measure_voltage - uref - V_rec) / norm(measure_voltage - uref); fprintf('相对残差: %.4f\n', rel_res);

把每步迭代的rel_res画出来：前几步快速下降说明J和ΔV构造正确；如果第一步残差就很小，那大概率是数据顺序或参考电位搞错了。如果残差一直高位徘徊，多半是电极处网格太粗。

5.2 把重建图像叠加到网格上看边界伪影

用trisurf把网格和电导率重建结果画在一起：

trisurf(mesh.tri, Node(:,1), Node(:,2), sigma_recon, 'EdgeColor', 'none'); view(2); colorbar; axis equal;

如果重建云图沿模型外边界出现一圈高亮亮带，说明电极边界条件和网格剖分不匹配。常见原因是点电极模型在电极附近产生过大灵敏度，而真实接触阻抗没有被建模。这时可以把电极周围的网格加密，或者改用完整电极模型。

5.3 一个能省几小时排查的J矩阵顺序探针

这套代码里最隐蔽的问题是J的行顺序与测量电压向量的顺序不一致。我的验证做法是：在某个固定单元k上把电导率提高1%，重新做一次正问题，得到新的电压向量V_new；然后把(V_new - V_ref)与J的第k列做线性相关检查。相关系数接近1说明J列与数据方向一致，若出现负相关或错位，就把drv/msr列表重新排列。

这个技巧不挑版本，无论源码包里的JacobianERT.m是伴随法还是差分法，都能用。运行一次正问题的时间通常比调半天错要短得多。建议每一步改完激励或测量顺序后都跑一次，直到图像不再出现棋盘格状错乱。

本文还有配套的精品资源，点击获取