电容层析成像Tikhonov算法

电容层析成像 (Electrical Capacitance Tomography, ECT) 中使用 Tikhonov 正则化算法 进行图像重建的演示实例。

目标： 从一组边界电容测量值中重建管道横截面内介电常数分布（例如，区分空气和油）。

---

演示实例：管道内两相流重建 (仿真)

1. 问题设置

• 成像区域： 一个圆形管道截面 (半径 R = 1)。
• 传感器： 8 个等间距贴附在管道外壁的电极 (N = 8)。
• 测量模式： 相邻激励-相邻测量 (相邻对模式)。总独立测量数 M = N*(N-1)/2 = 28。
• 真实分布 (仿真目标)： 管道底部有一个高介电常数 (ε_r = 3.0) 的矩形物体 (模拟油)，其余区域为低介电常数 (ε_r = 1.0，模拟空气)。
• 正问题模型： 使用线性近似 (灵敏度矩阵法)。
• 逆问题求解： Tikhonov 正则化。

2. MATLAB 代码

%% 步骤 1: 定义模型和网格
R = 1; % 管道半径
Nelec = 8; % 电极数量
M = Nelec*(Nelec-1)/2; % 独立测量数 (28)

% 创建成像区域的有限元网格 (例如使用 pdetool 或 distmesh)
% 这里简化表示：假设我们有一个包含 P 个像素或单元的网格
% pixels: [Npixels x 2] 存储每个像素中心的坐标
% elem: 定义单元连接性 (用于计算灵敏度)

% 定义真实介电常数分布 (ground truth)
epsilon_true = ones(Npixels, 1); % 初始化为空气 (ε_r=1)
% 在某个区域 (例如底部矩形) 设置高介电常数 (ε_r=3)
rect_indices = find(pixels(:,2) < -0.5 & abs(pixels(:,1)) < 0.4); % 找到矩形区域内的像素索引
epsilon_true(rect_indices) = 3.0;

%% 步骤 2: 计算灵敏度矩阵 S (核心，通常离线计算)
% S 是一个 [M x Npixels] 矩阵
% S(i, j) 表示第 i 个电容测量对第 j 个像素介电常数的微小变化的灵敏度

% 方法：
%   a) 使用有限元法 (FEM) 或有限差分法 (FDM) 求解正问题 (泊松方程)。
%   b) 在均匀背景场 (ε_r=1 everywhere) 下，计算所有电极对组合的电场分布 E_i (对应第 i 个激励模式)。
%   c) 对于每个像素 j，计算 S(i, j) = -∫∫_{Pixel_j} (E_i · E_i) dA / (V_i^2)
%      (根据互易定理和线性近似推导得出)
%   d) 对每个测量模式 i 重复步骤 b 和 c。

% 简化处理 (仅用于演示，实际需用FEM/FDM计算):
%   这里我们假设已经有一个计算好的 S 矩阵 (S_simulated.mat)。实际项目中需要用专业软件计算。
load('S_simulated.mat'); % 假设 S 已保存好，尺寸 [28 x Npixels]

%% 步骤 3: 生成模拟电容测量值 C_meas (正问题)
% 使用线性模型： ΔC = S * Δε
% ΔC: 相对于充满空气 (ε_r=1) 时电容的变化量 (向量, M x 1)
% Δε: 相对于空气的介电常数变化量 (向量, Npixels x 1) = epsilon_true - 1

delta_epsilon_true = epsilon_true - 1; % 介电常数变化量
C_meas = S * delta_epsilon_true; % 模拟电容变化量测量值 (线性近似)

% 添加噪声 (更符合实际)
noise_level = 0.01; % 1% 噪声
C_meas = C_meas + noise_level * max(abs(C_meas)) * randn(size(C_meas));

%% 步骤 4: 定义逆问题并应用 Tikhonov 正则化
% 目标：求解 Δε，使得 ||S * Δε - C_meas||² + λ² ||L * Δε||² 最小化
% λ: 正则化参数 (关键！)
% L: 正则化矩阵 (通常取单位矩阵 I 或梯度算子)

% 4.1 选择 L (这里使用单位矩阵 I，即零阶 Tikhonov / 岭回归)
L = eye(Npixels); % [Npixels x Npixels]

% 4.2 选择 λ (正则化参数) - 这是关键且困难的一步！
% 方法：L-曲线法 (L-curve) 或 广义交叉验证 (GCV) 或 试凑法 (trial-and-error)
% 这里为了演示，我们尝试几个不同的 λ 值观察效果
lambda_values = [1e-5, 1e-4, 1e-3, 1e-2, 0.1];

% 4.3 求解 Tikhonov 正则化最小二乘问题
% 解为： Δε_rec = (SᵀS + λ² LᵀL)⁻¹ Sᵀ C_meas
% 更稳定高效的解法是使用 SVD 或 直接构建方程求解

for i = 1:length(lambda_values)
    lambda = lambda_values(i);
    
    % 方法 1: 直接求解法 (对小规模问题可行)
    % A = S' * S + lambda^2 * (L' * L); % [Npixels x Npixels]
    % b = S' * C_meas; % [Npixels x 1]
    % delta_epsilon_rec = A \ b; % 求解线性方程组
    
    % 方法 2: 使用 SVD (更稳定，揭示问题本质)
    [U, Sigma, V] = svd(S, 'econ'); % S = U * Sigma * V'
    % Sigma 是 [M x Npixels] 的对角阵 (奇异值)
    diagS = diag(Sigma);
    % Tikhonov 滤波因子
    filter_factors = diagS.^2 ./ (diagS.^2 + lambda^2); % 当 L=I 时
    % 重建解
    delta_epsilon_rec = V * diag(1./diagS .* filter_factors) * U' * C_meas; 
    
    % 4.4 计算重建的介电常数分布
    epsilon_rec = 1 + delta_epsilon_rec; % 加上背景
    
    % 4.5 可视化重建结果 (对比真实分布)
    figure(i);
    subplot(1, 2, 1);
    scatter(pixels(:,1), pixels(:,2), 50, epsilon_true, 'filled');
    axis equal; axis([-R R -R R]); title('真实介电常数分布'); colorbar;
    
    subplot(1, 2, 2);
    scatter(pixels(:,1), pixels(:,2), 50, epsilon_rec, 'filled');
    axis equal; axis([-R R -R R]); title(['Tikhonov重建, \lambda = ', num2str(lambda)]); colorbar;
end

3. 预期结果与讨论

• λ 过小 (如 1e-5):
  • 解对噪声非常敏感。
  • 重建图像会出现严重的振荡、伪影和噪声放大，几乎无法识别目标物体。
  • 这是因为病态问题中小的奇异值被过度放大了噪声分量。
• λ 适中 (如 1e-4, 1e-3):
  • 噪声和伪影得到有效抑制。
  • 管道底部高介电常数的矩形区域能够被大致重建出来，位置和形状相对清晰。
  • 边界可能不够锐利，物体内部的介电常数分布可能不够均匀 (趋向于背景值)。
  • 这是通常希望达到的“最佳”或“可接受”的重建效果。L-曲线法通常能帮助找到这个范围内的 λ。
• λ 过大 (如 1e-2, 0.1):
  • 正则化过度，解过于平滑。
  • 重建图像变得非常模糊，目标物体的对比度显著降低，尺寸可能变小，边界完全消失，几乎与背景融为一体。
  • 丢失了图像细节和分辨率。

参考代码 电容层析成像Tikhonov算法的演示实例 www.youwenfan.com/contentcnc/54998.html

4. 关键点总结

1. 灵敏度矩阵 (S): ECT 图像重建的核心和计算瓶颈。通常需要精确的数值方法 (FEM/FDM) 离线计算。其病态性是逆问题困难的根源。
2. Tikhonov 正则化: 通过引入正则项 λ² ||L Δε||² 来稳定病态问题的求解。本质上是过滤掉小的奇异值对应的噪声主导的分量。
3. 正则化参数 λ: 控制着解的平滑度与对测量数据的拟合程度之间的平衡。选择 λ 至关重要：
   • 太小： 解不稳定，噪声放大 -> 欠正则化 (Under-regularization)。
   • 太大： 解过度平滑，丢失细节 -> 过正则化 (Over-regularization)。
   • 方法： L-曲线法 (寻找曲率最大的拐点)、广义交叉验证 (GCV)、试凑法、基于先验知识的启发式规则是最常用的方法。
4. 正则化矩阵 (L):
   • L = I (零阶): 惩罚解本身的范数 (大小)，倾向于小解。最简单常用。
   • L = 梯度算子 (一阶): 惩罚解的梯度 (变化)，倾向于平滑解。能更好地保留边缘 (理论上)，但计算更复杂，有时对噪声更敏感。
5. 线性模型局限性: 本例使用线性近似 ΔC = S * Δε。实际 ECT 是非线性的 (电容与介电常数关系非线性，灵敏度依赖于当前分布)。对于高对比度或复杂分布，非线性迭代算法 (如 Landweber, Newton-Raphson) 结合 Tikhonov 或其他正则化效果更好。
6. 网格与像素化: 重建结果是离散网格上的介电常数分布。网格的选择 (类型、大小) 会影响计算精度、速度和重建图像的表观质量。
7. 噪声: 实际测量必然包含噪声。Tikhonov 正则化的主要作用之一就是抑制噪声影响。噪声水平会影响最佳 λ 的选择。

5. 可视化示例 (概念图)

真实分布:      λ 过小 (e.g., 1e-5):       λ 适中 (e.g., 1e-3):       λ 过大 (e.g., 0.1):
+----------+    +----------+            +----------+            +----------+
|   Air    |    |**Noisy***|            |  Smooth  |            |   Very   |
|(ε_r=1)  |    |*Artifact*|            |  Object  |            |  Smooth  |
|  +----+  |    |***Oscill.|            |  at Bot  |            |   Blob   |
|  |Oil |  |    |***ation**|            | (ε_r≈2-3)|            |(ε_r≈1.5)|
|  |(ε=3)|  |    |*******__|            |  +----+  |            |__________|
|  +----+  |    |__________|            |__________|            |          |
+----------+                            |          |            |          |
                                       +----------+            +----------+
(清晰矩形在底部) (布满噪声和伪影)      (较平滑矩形在底部)         (几乎均匀的模糊团块)

---