简介面向电阻抗断层成像EIT逆问题研究的一套完整MATLAB源码包适用于生物医学成像、工业过程层析成像等领域的科研人员与研究生。压缩包共21个文件涵盖15个m脚本、3个mat数据文件、2个COMSOL模型文件mph及1个mphbin文件总计约41.67MB可直接在MATLAB与COMSOL环境下运行调试。核心算法覆盖吉洪诺夫正则化、Landweber迭代、L1稀疏重构及共轭梯度法CGLS并配有激发模式、雅可比矩阵计算、有限元节点处理等关键函数满足从数据加载、正问题构建到逆问题求解的完整流程。已有313人浏览学习适合需要快速对比多种正则化策略或验证EIT重建算法的使用者。资源内部包括Landweber、Tikhonov、L1等方法对应主程序以及fasta稀疏最小二乘工具便于扩展和改进。1. 从 COMSOL 场解到重建图像EIT 逆问题的完整路径说实话纯靠边界电压把内部的电导率分布还原出来这个问题困扰了 EIT电阻抗断层成像研究者三十多年。逆问题本身是非线性的加上灵敏度矩阵严重病态直接把测量方程组求逆的结果根本不能用。这套 EITtext_EIT 源码把问题拆成了几条可运行的求解路线Tikhonov 正则化、Landweber 迭代、CGLS 共轭梯度最小二乘以及用 FASTA 求解的 L1 稀疏重构每一条都有对应的入口脚本和配套 COMSOL 正问题模型。它的典型使用场景是你手里有一组电极电压测量数据或者想在仿真里验证某种正则化策略对重建分辨率的影响。适合正在做 EIT 成像仿真、刚接管实验台需要快速搭起重建管线的工程师和研究生。下面按「先立正问题再攻逆问题最后盯数据质量」的顺序把每个文件的实际作用讲透。2. 正问题与灵敏度矩阵Jacobian_ERT_2D 的组装逻辑EIT 重建不是直接从边界电压画彩图而是先通过有限元或 COMSOL 算出正问题场分布再用灵敏度矩阵把「电导率变化 → 边界电压变化」的线性化关系显式写出来。源码包里model0.mph、model1.mph 是 COMSOL 模型文件model0.m、model1.m 是配套的 LiveLink 脚本负责在 MATLAB 环境里调用模型并导出电极电位而 excitation_mode.m、node2D.m、uel.mat 则是把模型的物理输出组织成逆问题可用的矩阵和向量。2.1 激励模式与电极索引excitation_mode.m 里的测量安排激励模式决定了哪些电极对注入电流、哪些电极对测量电压。这套源码采用的是相邻激励-相邻测量模式也是 EIT 实验中最常见的设置每次只在相邻两个电极上注入恒定电流其余电极按相邻电极对依次测量电位差然后切换到下一对注入电极继续。16 电极系统一次完整扫描可以得到 16 组注入对和 16×(16−3) 个独立测量值。excitation_mode.m 干的就是生成这些索引对的工作% excitation_mode.m —— 生成相邻激励-相邻测量的电极对索引 function [stimPairs, measPairs] excitation_mode(nElec) stimPairs zeros(nElec, 2); for k 1:nElec stimPairs(k, :) [k, mod(k, nElec) 1]; % 电流注入电极对 end measPairs cell(nElec, 1); for k 1:nElec used [k, mod(k, nElec) 1]; used [used, mod(k-2, nElec) 1, mod(k2, nElec) 1]; idx setdiff(1:nElec, unique(used)); if numel(idx) 2 measPairs{k} [idx(1:end-1); idx(2:end)]; end end endstimPairs的每一行是一对注入电极measPairs是元胞数组第 k 个元素对应第 k 次激励时应该取哪些相邻电极对的电位差。这里用setdiff排除与注入电极共享节点的测量对是为了避免信号里混入接触阻抗带来的误差。实际使用中你可以在 8、16、32 电极之间切换索引会自动生长但要注意data_cirs.mat里测量向量的行顺序必须和measPairs的行顺序严格一致否则重建时会错位到错误的灵敏度行上。激励模式的选择直接改变灵敏度分布。相邻激励在边界附近的灵敏度最高中心区域显著衰减如果目标靠近成像区域中心可以改成相对激励或交叉激励代价是等效测量数变化、矩阵结构不再保持简单的循环对称。激励模式电流注入路径灵敏度分布特征相对噪声容限相邻激励相邻电极之间边界高、中心低中相对激励正对电极之间中心区域增强较高交叉激励间隔多个电极分布更均匀低2.2 节点与单元node2D.m 和 uel.mat 的映射关系灵敏度矩阵是按网格单元编号的。node2D.m 负责生成二维圆形成像域的节点坐标uel.mat 保存每个三角形单元由哪三个全局节点构成。假设你要把半径 0.15 m 的圆形区域剖成三角形网格常见做法是沿半径分层布置节点再连接成三角形单元。节点数决定了解空间的维度单元数决定了灵敏度矩阵的列数两者的比值直接影响逆问题的病态程度% 读入网格并查看维度 node node2D(0.15, 8); % 半径 0.15 m径向 8 层节点 load(uel.mat); % 每个单元由 3 个节点编号构成 fprintf(节点数 %d单元数 %d\n, size(node,1), size(uel,1));node2D第一个参数是成像域半径第二个参数是径向分层数。分层越多网格越密解空间维数越高灵敏度矩阵条件数也随之增大。uel.mat每一行是三角形单元的 3 个全局节点编号通常按逆时针排列这样面积和梯度计算都不会出现符号问题。你可以执行hist(node(:,2))看看节点在 y 方向的分布是否均匀——如果中心区域过稀而边界过密灵敏度矩阵各列的范数会差出几个数量级直接导致后期正则化参数在各区域失衡。2.3 Geselowitz 灵敏度定理与 Jacobian_ERT_2D 的实现灵敏度矩阵 J 的物理含义是第 m 个边界电压测量值对第 e 个单元内电导率变化的偏导数。二维 EIT 中最常用的计算依据是 Geselowitz 灵敏度定理当一个单元电导率发生微小扰动时边界电压的变化近似等于该单元内正向电场与伴随场梯度向量点积的积分。写成 MATLAB 循环就是 Jacobian_ERT_2D.m 的核心% Jacobian_ERT_2D.m —— 二维灵敏度矩阵的核心计算循环 function J Jacobian_ERT_2D(node, elem, phiA, phiB) % phiA: 正向激励下的节点电位场 % phiB: 伴随对偶激励下的节点电位场 nelem size(elem, 1); nmeas size(phiA, 1); J zeros(nmeas, nelem); for m 1:nmeas for e 1:nelem gradA shape_gradient(phiA(m,:), e, node, elem); gradB shape_gradient(phiB(m,:), e, node, elem); J(m, e) -gradA * gradB * elem_area(node, elem(e,:)); end end endshape_gradient根据三角形单元的三个顶点电位计算单元内的近似恒定梯度elem_area用叉积计算单元面积。这个双重循环的时间复杂度是 O(nmeas × nelem)网格加密后增长很快所以工程上通常把与测量无关的形函数梯度矩阵和面积提前算好缓存在uel.mat里避免每个测量向量进来都重新组装。拿到的 J 矩阵行数等于测量数列数等于单元数后面所有逆问题求解器都以它为基础。3. Tikhonov 与 Landweber经典正则化的两条实现路线有了灵敏度矩阵 J 和边界电压差 V核心问题变成解 Jx V。直接求伪逆会因条件数过高而失败边界上 1% 的测量噪声可能放大成内部 100% 甚至更高的电导率波动。正则化的本质是在「拟合测量值」和「约束解的结构」之间取折中。这一章拆解源码包里最基础的两条 l2 路线。3.1 吉洪诺夫正则化最小二乘与 L2 罚项Tikhonov 方法的思路非常直接目标函数写作min_x ||Jx - V||_2² α·||Lx||_2²其中 L 是正则化矩阵。L 取单位阵时只限制解的能量取离散差分矩阵时则限制解的平滑度。对应到 Tikhonov.m核心只有一行% Tikhonov.m —— 标准吉洪诺夫正则化 function sigma Tikhonov(J, V, alpha, L) if nargin 4 || isempty(L) L speye(size(J, 2)); % 默认零阶正则化 end sigma (J*J alpha*(L*L)) \ (J*V); endalpha是正则化参数它越大解越偏向零或越平滑它越小解越接近原始病态最小二乘解。需要注意L*L如果是稠密矩阵会破坏稀疏性所以工程上一般用稀疏差分矩阵构造 L让括号内保持稀疏结构。\运算在 MATLAB 中会自动选 Cholesky 分解当单元数超过一万个时建议先对J*J alpha*(L*L)做一次 LDL 分解并保留因子后续每组测量向量进来只需回代耗时从秒级降到毫秒级。提示alpha 建议按 10 倍步长做对数扫描从max(abs(J*V))附近往下衰减通常跨 3 到 4 个数量级就能看到解从过平滑过渡到噪声主导。3.2 半收敛特性Landweber 迭代的步长与停止条件Landweber 方法不直接解方程组而用梯度下降反复修正当前解迭代公式是x_{k1} x_k ω·J·(V - J·x_k)每次迭代只有矩阵乘法和向量加内存占用远低于 Tikhonov 的分解过程因此非常适合嵌入式或实时 EIT 系统。源码包里的 landweber.m 与 main_inverse_landweber.m 是这样配合的% landweber.m —— 步长固定的 Landweber 迭代 function sigma landweber(J, V, omega, iter, sigma0) sigma sigma0; r V - J * sigma; for k 1:iter sigma sigma omega * (J * r); r V - J * sigma; if norm(r) / norm(V) 1e-5 fprintf(第 %d 次迭代收敛\n, k); break; end end endomega是松弛因子必须满足 0 ω 2/||J||₂²一般取1/norm(J, 2)^2保证收敛但偏保守。iter是最大迭代次数。Landweber 有典型半收敛特性迭代前期误差快速下降后期开始拟合噪声所以不能等残差降到极低才停。我一般把 iter 限制在 20 到 200 之间先把残差曲线画出来选残差下降开始变缓的那个点作为停止位置比固定轮数可靠得多。3.3 两条路线的定位对比对比项TikhonovLandweber求解方式一次线性方程组求解迭代逼近核心参数alpha、Lomega、iter、初值内存占用存储分解因子仅 J 和向量边缘保持弱边缘被平滑适度受半收敛限制适用规模单元数 1 万以下较合适单元数几万以上更合适选择依据很实际Tikhonov 结果对 alpha 极敏感但算法确定性高同样参数两次运行结果完全一致Landweber 的初值影响收敛速度但不改变最终收敛的极小位置。快速验证数据是否正常先用 Tikhonov 跑一遍如果要实时迭代或后续接分辨率增强模块Landweber 更顺手。4. CGLS 与 L1 稀疏重构从共轭梯度到 FASTATikhonov 和 Landweber 都属于 l2 正则化默认假设解的能量分布均匀所以重建结果边缘偏模糊。当被测对象只包含少数电导率异常比如一块混凝土里的空洞、一个反应釜里的局部相变区l1 正则化在分辨率和伪影抑制上明显更优。源码包里的 main_inverse.m 和 main_inverse_L1.m 分别对应这两条路线。4.1 共轭梯度最小二乘cgls.m 的迭代逻辑CGLS 求解的仍然是 min||Jx - V||₂但通过 Krylov 子空间构造共轭方向收敛速度远快于 Landweber。它同样有迭代正则化特性早期迭代重建低频大尺度结构后续迭代逐步加入高频细节和噪声。cgls.m 的最小实现如下% cgls.m —— 共轭梯度最小二乘 function x cgls(J, V, iter, tol) x zeros(size(J, 2), 1); r V; % 残差初始为测量向量 p J * r; rho p * p; for k 1:iter q J * p; alpha rho / (q * q); x x alpha * p; r r - alpha * q; s J * r; rhoNew s * s; if sqrt(rhoNew) tol * norm(V) break; end p s (rhoNew / rho) * p; rho rhoNew; end end这里的p是残差在 J 空间的共轭方向alpha是最优步长。相比 Landweber 的反复震荡CGLS 通常几十步内就能达到同等残差水平。但 CGLS 对测量异常值非常敏感如果 V 里某一通道有坏点迭代早期重建图就会出现沿场方向的条状伪影。工程上建议先对测量向量做中值滤波或离群通道剔除再把 CGLS 当精细求解器用而不是唯一依赖。4.2 L1 最小化与近端算子fasta_sparseLeastSquares.m 的包装L1 回归的目标函数是min_x 0.5·||Jx - V||₂² λ·||x||₁l1 罚项利于产生稀疏解绝大多数单元的电导率扰动被压到接近零的位置只有少数单元扰动较大这正好匹配「个别异常目标」的物理直觉。fasta.m 是通用的近端梯度求解器fasta_sparseLeastSquares.m 在外面做了一层 EIT 接口% fasta_sparseLeastSquares.m —— L1 正则化 EIT 逆问题包装 function sigma fasta_sparseLeastSquares(J, V, lambda, maxIter) gradf (x) J * (J * x - V); % 平滑部分梯度 prox (x, t) sign(x) .* max(abs(x) - t * lambda, 0); % 近端算子不同版本 FASTA 参数签名有差异以下载包内 fasta.m 为准 sigma fasta(gradf, prox, maxIter, ... backtracking, true); endgradf是数据保真项梯度prox是 l1 近端算子软阈值backtracking开启后每轮迭代会自动调整步长省去手动调松弛因子的步骤。lambda 的物理含义是稀疏度阈值lambda 越接近max(abs(J*V))解越稀疏降到其千分之一以下解基本退化成 CGLS 的样子。实际操作中我会先算corr J*V看它绝对值分布把 lambda 放在[0.05*max(abs(corr)), 0.5*max(abs(corr))]区间这样能兼顾稀疏性和目标不漏检。4.3 LBP 线性反投影最快速的粗定位手段lbp.m 走的是极简路线把灵敏度矩阵转置直接乘测量向量% lbp.m —— 线性反投影快速重建 function sigma lbp(J, V) Jn J ./ sum(abs(J), 2); % 行归一化削弱灵敏度不均影响 sigma Jn * V; end没有正则化没有迭代重建结果在数值上不保真但异常区域的质心位置一般不会跑偏。它的价值体现在两处一是作为迭代法初值能让 CGLS 和 Landweber 少迭代一半以上二是快速检查测量通道是否存在接反或断路——LBP 图出现扇形状伪影时多半是某个电极的相位反了。这个技巧在调实验数据时比任何高级求解器都管用。方法目标函数输出特性关键参数CGLSminJx-VL1-FASTAmin 0.5Jx-VLBPxJV只定性、假幅值行归一化策略5. 数据流串起来之后datasave.m 与重建质量验证前几章拆开了每个文件现在把它们串成一个完整流程。main_inverse.m 是典型入口加载 data_cirs.mat 得到仿真测量电压加载 e_xyz.mat 得到坐标和参考分布调用 excitation_mode.m 拿到激励索引再按路线进入 Tikhonov、CGLS 或 Landweber 求解。datasave.m 负责把解出的电导率分布写成可被 COMSOL 重新导入或直接用 imagesc 绘图的格式。5.1 从加载数据到输出图像脚本依赖顺序% 串起完整流程的最小示例 load(data_cirs.mat); % 包含 v_meas: 测量电压差向量 load(e_xyz.mat); % 包含 x y 坐标与电导率真值 J Jacobian_ERT_2D(node, elem, phiA, phiB); sigma_l1 fasta_sparseLeastSquares(J, v_meas, 1e-2, 300); datasave(result_l1.mat, sigma_l1);注意顺序必须先组装 J再载入测量向量最后调用求解器。很多新手把data_cirs.mat里的数组当成电压值直接用忽略了它可能还包含激励编号和测量编号导致与 J 的行数对不上建议第一步就检查size(v_meas,1) size(J,1)。5.2 用 LBP 初值加速迭代类求解器迭代类求解器对初值的依赖比参数更隐蔽。把lbp(J, V)的输出作为 Landweber 或 CGLS 的sigma0可以在相同迭代轮数内明显提高目标区域灰度sigma0 lbp(J, v_meas); sigma_cgls cgls(J, v_meas, 30, 1e-6); % 对比两种初值下的重建均方误差确认是否值得引入 LBP 初值注意LBP 初值只用于加速不要把它当作定量结果直接出报告。5.3 三点快速判断重建结果是否可信判断重建质量不要只看图像像不像。第一计算重建目标质心与真实目标质心的距离偏差超过一个相邻电极间距说明要么激励模式选错要么测量通道顺序错位。第二统计目标区域外重建值的最大振幅这是伪影强度正常应低于目标区域均值的 30%超过就需要加大正则化参数。第三对比目标区域与背景区域的灰度均值比这个比值低于 1.5 时说明正则化过强目标被拉平过高则说明噪声主导了重建。这三个指标互相制衡比肉眼看伪彩色图可靠得多。本文还有配套的精品资源点击获取
EIT逆问题求解全解析：从COMSOL建模到正则化图像重建
简介面向电阻抗断层成像EIT逆问题研究的一套完整MATLAB源码包适用于生物医学成像、工业过程层析成像等领域的科研人员与研究生。压缩包共21个文件涵盖15个m脚本、3个mat数据文件、2个COMSOL模型文件mph及1个mphbin文件总计约41.67MB可直接在MATLAB与COMSOL环境下运行调试。核心算法覆盖吉洪诺夫正则化、Landweber迭代、L1稀疏重构及共轭梯度法CGLS并配有激发模式、雅可比矩阵计算、有限元节点处理等关键函数满足从数据加载、正问题构建到逆问题求解的完整流程。已有313人浏览学习适合需要快速对比多种正则化策略或验证EIT重建算法的使用者。资源内部包括Landweber、Tikhonov、L1等方法对应主程序以及fasta稀疏最小二乘工具便于扩展和改进。1. 从 COMSOL 场解到重建图像EIT 逆问题的完整路径说实话纯靠边界电压把内部的电导率分布还原出来这个问题困扰了 EIT电阻抗断层成像研究者三十多年。逆问题本身是非线性的加上灵敏度矩阵严重病态直接把测量方程组求逆的结果根本不能用。这套 EITtext_EIT 源码把问题拆成了几条可运行的求解路线Tikhonov 正则化、Landweber 迭代、CGLS 共轭梯度最小二乘以及用 FASTA 求解的 L1 稀疏重构每一条都有对应的入口脚本和配套 COMSOL 正问题模型。它的典型使用场景是你手里有一组电极电压测量数据或者想在仿真里验证某种正则化策略对重建分辨率的影响。适合正在做 EIT 成像仿真、刚接管实验台需要快速搭起重建管线的工程师和研究生。下面按「先立正问题再攻逆问题最后盯数据质量」的顺序把每个文件的实际作用讲透。2. 正问题与灵敏度矩阵Jacobian_ERT_2D 的组装逻辑EIT 重建不是直接从边界电压画彩图而是先通过有限元或 COMSOL 算出正问题场分布再用灵敏度矩阵把「电导率变化 → 边界电压变化」的线性化关系显式写出来。源码包里model0.mph、model1.mph 是 COMSOL 模型文件model0.m、model1.m 是配套的 LiveLink 脚本负责在 MATLAB 环境里调用模型并导出电极电位而 excitation_mode.m、node2D.m、uel.mat 则是把模型的物理输出组织成逆问题可用的矩阵和向量。2.1 激励模式与电极索引excitation_mode.m 里的测量安排激励模式决定了哪些电极对注入电流、哪些电极对测量电压。这套源码采用的是相邻激励-相邻测量模式也是 EIT 实验中最常见的设置每次只在相邻两个电极上注入恒定电流其余电极按相邻电极对依次测量电位差然后切换到下一对注入电极继续。16 电极系统一次完整扫描可以得到 16 组注入对和 16×(16−3) 个独立测量值。excitation_mode.m 干的就是生成这些索引对的工作% excitation_mode.m —— 生成相邻激励-相邻测量的电极对索引 function [stimPairs, measPairs] excitation_mode(nElec) stimPairs zeros(nElec, 2); for k 1:nElec stimPairs(k, :) [k, mod(k, nElec) 1]; % 电流注入电极对 end measPairs cell(nElec, 1); for k 1:nElec used [k, mod(k, nElec) 1]; used [used, mod(k-2, nElec) 1, mod(k2, nElec) 1]; idx setdiff(1:nElec, unique(used)); if numel(idx) 2 measPairs{k} [idx(1:end-1); idx(2:end)]; end end endstimPairs的每一行是一对注入电极measPairs是元胞数组第 k 个元素对应第 k 次激励时应该取哪些相邻电极对的电位差。这里用setdiff排除与注入电极共享节点的测量对是为了避免信号里混入接触阻抗带来的误差。实际使用中你可以在 8、16、32 电极之间切换索引会自动生长但要注意data_cirs.mat里测量向量的行顺序必须和measPairs的行顺序严格一致否则重建时会错位到错误的灵敏度行上。激励模式的选择直接改变灵敏度分布。相邻激励在边界附近的灵敏度最高中心区域显著衰减如果目标靠近成像区域中心可以改成相对激励或交叉激励代价是等效测量数变化、矩阵结构不再保持简单的循环对称。激励模式电流注入路径灵敏度分布特征相对噪声容限相邻激励相邻电极之间边界高、中心低中相对激励正对电极之间中心区域增强较高交叉激励间隔多个电极分布更均匀低2.2 节点与单元node2D.m 和 uel.mat 的映射关系灵敏度矩阵是按网格单元编号的。node2D.m 负责生成二维圆形成像域的节点坐标uel.mat 保存每个三角形单元由哪三个全局节点构成。假设你要把半径 0.15 m 的圆形区域剖成三角形网格常见做法是沿半径分层布置节点再连接成三角形单元。节点数决定了解空间的维度单元数决定了灵敏度矩阵的列数两者的比值直接影响逆问题的病态程度% 读入网格并查看维度 node node2D(0.15, 8); % 半径 0.15 m径向 8 层节点 load(uel.mat); % 每个单元由 3 个节点编号构成 fprintf(节点数 %d单元数 %d\n, size(node,1), size(uel,1));node2D第一个参数是成像域半径第二个参数是径向分层数。分层越多网格越密解空间维数越高灵敏度矩阵条件数也随之增大。uel.mat每一行是三角形单元的 3 个全局节点编号通常按逆时针排列这样面积和梯度计算都不会出现符号问题。你可以执行hist(node(:,2))看看节点在 y 方向的分布是否均匀——如果中心区域过稀而边界过密灵敏度矩阵各列的范数会差出几个数量级直接导致后期正则化参数在各区域失衡。2.3 Geselowitz 灵敏度定理与 Jacobian_ERT_2D 的实现灵敏度矩阵 J 的物理含义是第 m 个边界电压测量值对第 e 个单元内电导率变化的偏导数。二维 EIT 中最常用的计算依据是 Geselowitz 灵敏度定理当一个单元电导率发生微小扰动时边界电压的变化近似等于该单元内正向电场与伴随场梯度向量点积的积分。写成 MATLAB 循环就是 Jacobian_ERT_2D.m 的核心% Jacobian_ERT_2D.m —— 二维灵敏度矩阵的核心计算循环 function J Jacobian_ERT_2D(node, elem, phiA, phiB) % phiA: 正向激励下的节点电位场 % phiB: 伴随对偶激励下的节点电位场 nelem size(elem, 1); nmeas size(phiA, 1); J zeros(nmeas, nelem); for m 1:nmeas for e 1:nelem gradA shape_gradient(phiA(m,:), e, node, elem); gradB shape_gradient(phiB(m,:), e, node, elem); J(m, e) -gradA * gradB * elem_area(node, elem(e,:)); end end endshape_gradient根据三角形单元的三个顶点电位计算单元内的近似恒定梯度elem_area用叉积计算单元面积。这个双重循环的时间复杂度是 O(nmeas × nelem)网格加密后增长很快所以工程上通常把与测量无关的形函数梯度矩阵和面积提前算好缓存在uel.mat里避免每个测量向量进来都重新组装。拿到的 J 矩阵行数等于测量数列数等于单元数后面所有逆问题求解器都以它为基础。3. Tikhonov 与 Landweber经典正则化的两条实现路线有了灵敏度矩阵 J 和边界电压差 V核心问题变成解 Jx V。直接求伪逆会因条件数过高而失败边界上 1% 的测量噪声可能放大成内部 100% 甚至更高的电导率波动。正则化的本质是在「拟合测量值」和「约束解的结构」之间取折中。这一章拆解源码包里最基础的两条 l2 路线。3.1 吉洪诺夫正则化最小二乘与 L2 罚项Tikhonov 方法的思路非常直接目标函数写作min_x ||Jx - V||_2² α·||Lx||_2²其中 L 是正则化矩阵。L 取单位阵时只限制解的能量取离散差分矩阵时则限制解的平滑度。对应到 Tikhonov.m核心只有一行% Tikhonov.m —— 标准吉洪诺夫正则化 function sigma Tikhonov(J, V, alpha, L) if nargin 4 || isempty(L) L speye(size(J, 2)); % 默认零阶正则化 end sigma (J*J alpha*(L*L)) \ (J*V); endalpha是正则化参数它越大解越偏向零或越平滑它越小解越接近原始病态最小二乘解。需要注意L*L如果是稠密矩阵会破坏稀疏性所以工程上一般用稀疏差分矩阵构造 L让括号内保持稀疏结构。\运算在 MATLAB 中会自动选 Cholesky 分解当单元数超过一万个时建议先对J*J alpha*(L*L)做一次 LDL 分解并保留因子后续每组测量向量进来只需回代耗时从秒级降到毫秒级。提示alpha 建议按 10 倍步长做对数扫描从max(abs(J*V))附近往下衰减通常跨 3 到 4 个数量级就能看到解从过平滑过渡到噪声主导。3.2 半收敛特性Landweber 迭代的步长与停止条件Landweber 方法不直接解方程组而用梯度下降反复修正当前解迭代公式是x_{k1} x_k ω·J·(V - J·x_k)每次迭代只有矩阵乘法和向量加内存占用远低于 Tikhonov 的分解过程因此非常适合嵌入式或实时 EIT 系统。源码包里的 landweber.m 与 main_inverse_landweber.m 是这样配合的% landweber.m —— 步长固定的 Landweber 迭代 function sigma landweber(J, V, omega, iter, sigma0) sigma sigma0; r V - J * sigma; for k 1:iter sigma sigma omega * (J * r); r V - J * sigma; if norm(r) / norm(V) 1e-5 fprintf(第 %d 次迭代收敛\n, k); break; end end endomega是松弛因子必须满足 0 ω 2/||J||₂²一般取1/norm(J, 2)^2保证收敛但偏保守。iter是最大迭代次数。Landweber 有典型半收敛特性迭代前期误差快速下降后期开始拟合噪声所以不能等残差降到极低才停。我一般把 iter 限制在 20 到 200 之间先把残差曲线画出来选残差下降开始变缓的那个点作为停止位置比固定轮数可靠得多。3.3 两条路线的定位对比对比项TikhonovLandweber求解方式一次线性方程组求解迭代逼近核心参数alpha、Lomega、iter、初值内存占用存储分解因子仅 J 和向量边缘保持弱边缘被平滑适度受半收敛限制适用规模单元数 1 万以下较合适单元数几万以上更合适选择依据很实际Tikhonov 结果对 alpha 极敏感但算法确定性高同样参数两次运行结果完全一致Landweber 的初值影响收敛速度但不改变最终收敛的极小位置。快速验证数据是否正常先用 Tikhonov 跑一遍如果要实时迭代或后续接分辨率增强模块Landweber 更顺手。4. CGLS 与 L1 稀疏重构从共轭梯度到 FASTATikhonov 和 Landweber 都属于 l2 正则化默认假设解的能量分布均匀所以重建结果边缘偏模糊。当被测对象只包含少数电导率异常比如一块混凝土里的空洞、一个反应釜里的局部相变区l1 正则化在分辨率和伪影抑制上明显更优。源码包里的 main_inverse.m 和 main_inverse_L1.m 分别对应这两条路线。4.1 共轭梯度最小二乘cgls.m 的迭代逻辑CGLS 求解的仍然是 min||Jx - V||₂但通过 Krylov 子空间构造共轭方向收敛速度远快于 Landweber。它同样有迭代正则化特性早期迭代重建低频大尺度结构后续迭代逐步加入高频细节和噪声。cgls.m 的最小实现如下% cgls.m —— 共轭梯度最小二乘 function x cgls(J, V, iter, tol) x zeros(size(J, 2), 1); r V; % 残差初始为测量向量 p J * r; rho p * p; for k 1:iter q J * p; alpha rho / (q * q); x x alpha * p; r r - alpha * q; s J * r; rhoNew s * s; if sqrt(rhoNew) tol * norm(V) break; end p s (rhoNew / rho) * p; rho rhoNew; end end这里的p是残差在 J 空间的共轭方向alpha是最优步长。相比 Landweber 的反复震荡CGLS 通常几十步内就能达到同等残差水平。但 CGLS 对测量异常值非常敏感如果 V 里某一通道有坏点迭代早期重建图就会出现沿场方向的条状伪影。工程上建议先对测量向量做中值滤波或离群通道剔除再把 CGLS 当精细求解器用而不是唯一依赖。4.2 L1 最小化与近端算子fasta_sparseLeastSquares.m 的包装L1 回归的目标函数是min_x 0.5·||Jx - V||₂² λ·||x||₁l1 罚项利于产生稀疏解绝大多数单元的电导率扰动被压到接近零的位置只有少数单元扰动较大这正好匹配「个别异常目标」的物理直觉。fasta.m 是通用的近端梯度求解器fasta_sparseLeastSquares.m 在外面做了一层 EIT 接口% fasta_sparseLeastSquares.m —— L1 正则化 EIT 逆问题包装 function sigma fasta_sparseLeastSquares(J, V, lambda, maxIter) gradf (x) J * (J * x - V); % 平滑部分梯度 prox (x, t) sign(x) .* max(abs(x) - t * lambda, 0); % 近端算子不同版本 FASTA 参数签名有差异以下载包内 fasta.m 为准 sigma fasta(gradf, prox, maxIter, ... backtracking, true); endgradf是数据保真项梯度prox是 l1 近端算子软阈值backtracking开启后每轮迭代会自动调整步长省去手动调松弛因子的步骤。lambda 的物理含义是稀疏度阈值lambda 越接近max(abs(J*V))解越稀疏降到其千分之一以下解基本退化成 CGLS 的样子。实际操作中我会先算corr J*V看它绝对值分布把 lambda 放在[0.05*max(abs(corr)), 0.5*max(abs(corr))]区间这样能兼顾稀疏性和目标不漏检。4.3 LBP 线性反投影最快速的粗定位手段lbp.m 走的是极简路线把灵敏度矩阵转置直接乘测量向量% lbp.m —— 线性反投影快速重建 function sigma lbp(J, V) Jn J ./ sum(abs(J), 2); % 行归一化削弱灵敏度不均影响 sigma Jn * V; end没有正则化没有迭代重建结果在数值上不保真但异常区域的质心位置一般不会跑偏。它的价值体现在两处一是作为迭代法初值能让 CGLS 和 Landweber 少迭代一半以上二是快速检查测量通道是否存在接反或断路——LBP 图出现扇形状伪影时多半是某个电极的相位反了。这个技巧在调实验数据时比任何高级求解器都管用。方法目标函数输出特性关键参数CGLSminJx-VL1-FASTAmin 0.5Jx-VLBPxJV只定性、假幅值行归一化策略5. 数据流串起来之后datasave.m 与重建质量验证前几章拆开了每个文件现在把它们串成一个完整流程。main_inverse.m 是典型入口加载 data_cirs.mat 得到仿真测量电压加载 e_xyz.mat 得到坐标和参考分布调用 excitation_mode.m 拿到激励索引再按路线进入 Tikhonov、CGLS 或 Landweber 求解。datasave.m 负责把解出的电导率分布写成可被 COMSOL 重新导入或直接用 imagesc 绘图的格式。5.1 从加载数据到输出图像脚本依赖顺序% 串起完整流程的最小示例 load(data_cirs.mat); % 包含 v_meas: 测量电压差向量 load(e_xyz.mat); % 包含 x y 坐标与电导率真值 J Jacobian_ERT_2D(node, elem, phiA, phiB); sigma_l1 fasta_sparseLeastSquares(J, v_meas, 1e-2, 300); datasave(result_l1.mat, sigma_l1);注意顺序必须先组装 J再载入测量向量最后调用求解器。很多新手把data_cirs.mat里的数组当成电压值直接用忽略了它可能还包含激励编号和测量编号导致与 J 的行数对不上建议第一步就检查size(v_meas,1) size(J,1)。5.2 用 LBP 初值加速迭代类求解器迭代类求解器对初值的依赖比参数更隐蔽。把lbp(J, V)的输出作为 Landweber 或 CGLS 的sigma0可以在相同迭代轮数内明显提高目标区域灰度sigma0 lbp(J, v_meas); sigma_cgls cgls(J, v_meas, 30, 1e-6); % 对比两种初值下的重建均方误差确认是否值得引入 LBP 初值注意LBP 初值只用于加速不要把它当作定量结果直接出报告。5.3 三点快速判断重建结果是否可信判断重建质量不要只看图像像不像。第一计算重建目标质心与真实目标质心的距离偏差超过一个相邻电极间距说明要么激励模式选错要么测量通道顺序错位。第二统计目标区域外重建值的最大振幅这是伪影强度正常应低于目标区域均值的 30%超过就需要加大正则化参数。第三对比目标区域与背景区域的灰度均值比这个比值低于 1.5 时说明正则化过强目标被拉平过高则说明噪声主导了重建。这三个指标互相制衡比肉眼看伪彩色图可靠得多。本文还有配套的精品资源点击获取
