用Matlab一步步复现MRI并行成像SENSE算法：从k空间欠采样到图像重建的保姆级教程

从零实现MRI并行成像SENSE算法：Matlab实战指南与深度调优

开篇：为什么选择SENSE算法动手实践？

在医学影像领域，磁共振成像（MRI）的扫描速度一直是制约临床应用的瓶颈。传统序列扫描需要患者保持静止长达数十分钟，而并行成像技术通过线圈空间敏感度的差异，实现了k空间数据的欠采样，将扫描时间缩短至原来的1/2甚至1/8。SENSE（Sensitivity Encoding）作为最具代表性的并行成像算法之一，其核心思想是利用多个接收线圈的空间敏感度信息来解混叠图像。

动手实现SENSE算法的价值：

• 深入理解线圈敏感度与图像重建的数学关系
• 掌握处理复数k空间数据的实际技巧
• 学习调试医学图像重建中的典型问题
• 为后续研究GRAPPA等高级算法打下基础

本文将使用Matlab从零开始，完整实现一个5倍加速的SENSE重建流程，包含以下关键环节：

% 基础环境准备 clear; close all; clc; addpath(genpath('.')); % 添加工具函数路径

1. 数据准备与敏感度图生成

1.1 加载多通道MRI数据

我们使用包含5个线圈的大脑MRI数据，每个线圈采集的图像尺寸为120×128。原始数据包含两个关键变量：

• brain_coil_tmp: 各线圈的全采样MR图像（120×128×5）
• coil_map: 预扫描得到的线圈敏感度图（120×128×5）

% 加载.mat数据文件 brainCoilsData = load('brain_coil.mat'); brainCoils = brainCoilsData.brain_coil_tmp; coilMapData = load('coil_sensitivity_map.mat'); coilMap = coilMapData.coil_map;

数据可视化技巧： 使用subplot同时显示多个线圈的图像，便于比较各通道信号差异：

figure('Name','各线圈全采样MR图像'); for i = 1:5 subplot(2,3,i); imagesc(abs(brainCoils(:,:,i))); title(['线圈 ',num2str(i)]); colormap gray; axis image; colorbar; end

1.2 敏感度图计算原理

线圈敏感度图反映各空间位置对线圈信号的接收效率差异。理想情况下，敏感度图应满足：

$$ C_j(\mathbf{r}) = \frac{S_j(\mathbf{r})}{S_{\text{body}}(\mathbf{r})} $$

其中$S_j$为第j个线圈图像，$S_{\text{body}}$为所有线圈的平方和根（RSOS）重建图像。

手动计算敏感度图的Matlab实现：

BodybrainCoils = sqrt(sum(abs(brainCoils).^2, 3)); % RSOS重建 sensitivity_maps = zeros(size(brainCoils)); for i = 1:5 sensitivity_maps(:,:,i) = brainCoils(:,:,i) ./ BodybrainCoils; end

注意：实际应用中，敏感度图通常通过低分辨率预扫描数据计算，并经过滤波等后处理以提高稳定性。

2. k空间操作与欠采样模拟

2.1 全采样k空间转换

将图像域数据转换到k空间是重建算法的关键步骤。需要注意FFT前后的shift操作：

fullKspace = ifftshift(fft2(fftshift(brainCoils)));

为什么需要fftshift？

• MRI设备的k空间数据采集通常以中心为原点
• fftshift将直流分量移到频谱中心
• ifftshift是fftshift的逆操作

2.2 构造欠采样掩模

设置加速因子R=5，构造等距欠采样模式：

R = 5; % 加速因子 mask = zeros(size(fullKspace)); for line = 1:size(mask,1) if mod(line, R) == 0 mask(line,:,:) = 1; % 采集线标记为1 end end

欠采样k空间通过元素乘获得：

downSampledKspace = fullKspace .* mask;

欠采样模式对比表：

3. SENSE核心重建算法实现

3.1 数学原理剖析

SENSE重建的核心是求解线性方程组：

$$ \mathbf{I} = \mathbf{C} \mathbf{\rho} $$

其中：

• $\mathbf{I}$是各线圈测量的混叠图像向量（尺寸：Nc×1）
• $\mathbf{C}$是敏感度矩阵（尺寸：Nc×R）
• $\mathbf{\rho}$是待求解的完整图像像素（尺寸：R×1）

通过伪逆（pseudo-inverse）求解：

$$ \mathbf{\rho} = (\mathbf{C}^H \mathbf{C})^{-1} \mathbf{C}^H \mathbf{I} $$

3.2 Matlab逐像素实现

[phaseLength, freqLength, coilNum] = size(brainCoils); fieldOfView = floor(phaseLength/R); senseRecon = zeros(phaseLength, freqLength); for x = 1:freqLength for y = 1:fieldOfView % 获取当前像素各线圈测量值 vectorI = reshape(dsBrainCoils(y,x,:), coilNum, 1); % 构建敏感度矩阵 cHat = zeros(coilNum, R); for k = 1:R pos = y + (k-1)*fieldOfView; cHat(:,k) = reshape(coilMap(pos,x,:), coilNum, 1); end % 伪逆求解 cHatPinv = pinv(cHat); vectorRho = cHatPinv * vectorI; % 填充重建结果 for k = 1:R senseRecon(y+(k-1)*fieldOfView, x) = vectorRho(k); end end end

关键调试经验：

1. 矩阵维度必须严格匹配（特别是复数数据）
2. 伪逆求解对敏感度图精度敏感
3. 最终结果需要乘以R补偿信号强度

4. 结果评估与优化技巧

4.1 重建质量可视化

original = sqrt(sum(abs(brainCoils).^2, 3)); dsImage = sqrt(sum(abs(dsBrainCoils).^2, 3)); senseImage = abs(senseRecon) * R; figure; subplot(1,3,1); imshow(original,[]); title('原始图像'); subplot(1,3,2); imshow(dsImage,[]); title('欠采样混叠图像'); subplot(1,3,3); imshow(senseImage,[]); title('SENSE重建');

4.2 常见问题解决方案

问题1：边缘伪影

• 原因：敏感度图在边缘区域信噪比低
• 解决：对敏感度图进行低通滤波

问题2：噪声放大

• 原因：几何因子（g-factor）在敏感度变化剧烈区域较大
• 解决：加入正则化项，修改伪逆计算：

lambda = 0.1; % 正则化参数 cHatPinv = (cHat'*cHat + lambda*eye(R)) \ cHat';

问题3：重建时间过长

• 优化：将逐像素循环改为矩阵运算
• 利用Matlab的并行计算工具箱

5. 高级扩展与实战建议

5.1 不同采样模式对比

除等距采样外，还可实现：

1. 随机采样：减少相干伪影

   mask = rand(phaseLength,1) < 1/R; mask = repmat(mask, [1 freqLength coilNum]);

2. 中心密集采样：保留更多低频信息

5.2 自动敏感度图校准

实际应用中，可通过以下方式改进敏感度图：

1. 低分辨率全采样校准扫描
2. 自适应校准（ACS lines）
3. 迭代联合估计敏感度图与图像

5.3 与现代深度学习结合

传统SENSE的局限催生了以下研究方向：

1. DL-SENSE：用神经网络学习重建映射

   # 伪代码示例 model = UNet() output = model(torch.cat([under_sampled, sensitivity_maps], dim=1))

2. 端到端优化：联合优化采样模式与重建算法

调试备忘录：那些年踩过的坑

1. 复数数据处理：忘记取模导致图像显示异常

   % 错误做法 imagesc(senseRecon); % 正确做法 imagesc(abs(senseRecon));

2. 矩阵维度不匹配：reshape操作遗漏线圈维度

3. fftshift遗漏：导致k空间中心错位

4. 敏感度图归一化：避免除以零风险

   BodybrainCoils = BodybrainCoils + eps; % 添加极小值

5. 加速因子选择：R超过线圈数量会导致重建失败