用Matlab手把手复现MRI并行成像SENSE算法：从k空间欠采样到图像重建全流程-编程实验室

用Matlab手把手复现MRI并行成像SENSE算法：从k空间欠采样到图像重建全流程

在医学影像领域，磁共振成像（MRI）技术的进步始终围绕着两个核心目标：提升图像质量与缩短扫描时间。传统MRI扫描需要完整采集k空间数据，而并行成像技术通过利用多个接收线圈的空间敏感度信息，实现了k空间的欠采样，显著提高了成像效率。本文将带您深入SENSE（SENSitivity Encoding）算法的实现细节，通过Matlab代码逐步演示从k空间欠采样到最终图像重建的全过程。

1. 环境准备与数据加载

1.1 实验数据说明

我们使用的实验数据包含两个关键部分：

brain_coil.mat：5通道全采样大脑MR图像，尺寸120×128×5
coil_sensitivity_map.mat：对应5个线圈的敏感度图，尺寸与MR图像相同

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

1.2 数据可视化

在开始处理前，先观察原始数据特征：

% 显示各线圈全采样MR图像 figure; for i = 1:5 subplot(2,3,i); imagesc(brainCoils(:,:,i)); title(['Coil-',num2str(i),' MR图像']); colormap('gray'); colorbar; end % 显示各线圈敏感度图 figure; for i = 1:5 subplot(2,3,i); imagesc(coilMap(:,:,i)); title(['Coil-',num2str(i),'敏感度图']); colormap('gray'); colorbar; end

关键观察点：

不同线圈的图像呈现区域亮度差异，反映了各线圈的空间敏感度分布
敏感度图显示了各线圈对不同解剖区域的信号接收能力

2. k空间转换与欠采样策略

2.1 全采样k空间转换

将图像域数据转换到k空间是MRI处理的基础步骤：

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

这里需要注意fftshift和ifftshift的正确使用，确保k空间中心对齐。转换后，我们可以观察各线圈的k空间特征：

线圈编号	k空间特征	能量分布
1	中心亮区明显	高频成分较少
2	能量分布均匀	中高频成分丰富
3	边缘信号较强	各频段均衡
4	中心集中	低频主导
5	对角分布	各向异性明显

2.2 欠采样掩模设计

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

downSamplingRate = 5; mask = zeros(size(fullKspace)); for line = 1:size(mask,1) if rem(line, downSamplingRate) == 0 mask(line,:,:) = 1; end end

这种采样模式在相位编码方向（ky）每隔R-1行采集一行，实现了R倍的扫描加速。欠采样会导致k空间信息缺失，进而引起图像域的混叠伪影。

3. 敏感度图计算与验证

3.1 敏感度图生成原理

敏感度图反映了各线圈的空间响应特性，其计算通常通过低分辨率扫描数据获得。在我们的实现中：

BodybrainCoils = rsos(brainCoils); % 全采样组合图像 for i = 1:5 sensitivity_map(:,:,i) = brainCoils(:,:,i) ./ BodybrainCoils; end

敏感度图的关键特性：

在靠近线圈的区域值较大
不同线圈的敏感度分布互补
需进行适当的平滑和归一化处理

3.2 敏感度图验证

通过两种方法生成的敏感度图对比：

直接计算法：简单但易受噪声影响
k空间法：更接近实际MRI系统采集流程

% k空间法生成敏感度图 rawfullKspace = fullKspace; acsBrainCoils = ifftshift(ifft2(fftshift(rawfullKspace))); acsImage = rsos(acsBrainCoils); for i = 1:5 sens_map(:,:,i) = rsos(acsBrainCoils(:,:,i)) ./ acsImage; end

实验表明，k空间法生成的敏感度图与预扫描结果更为接近，噪点更少，对比度更合理。

4. SENSE重建核心算法

4.1 重建数学模型

SENSE算法的核心是求解以下线性方程组：

I = C · ρ

其中：

I是观测到的混叠图像向量（尺寸：线圈数×1）
C是敏感度矩阵（尺寸：线圈数×R）
ρ是待求解的完整图像向量（尺寸：R×1）

对于R=5，Nc=5的情况，这是一个适定方程组，可通过伪逆求解：

cHatPinv = pinv(cHat); vectorRho = cHatPinv * vectorI;

4.2 Matlab实现细节

完整的SENSE重建实现包含以下步骤：

初始化重建矩阵
遍历图像每个像素位置
构建敏感度矩阵
求解线性方程组
填充重建结果

fieldOfView = floor(phaseLength/downSamplingRate); senseRecon = zeros(phaseLength, freqLength); for x = 1:freqLength for y = 1:fieldOfView % 获取当前像素的混叠信号 vectorI = reshape(dsBrainCoils(y, x, :), coilNum, 1); % 构建敏感度矩阵 for k = 1:downSamplingRate cHat(:, k) = reshape(coilMap(y+(k-1)*fieldOfView, x, :), coilNum, 1); end % 伪逆求解 cHatPinv = pinv(cHat); vectorRho = cHatPinv * vectorI; % 填充重建结果 for k = 1:downSamplingRate senseRecon(y+(k-1)*fieldOfView, x) = vectorRho(k); end end end

4.3 结果可视化与评估

重建完成后，我们需要评估重建质量：

senseImage = rsos(senseRecon) * downSamplingRate; original = rsos(brainCoils); delta = original - senseImage; figure; subplot(2,2,1); imagesc(original); title('原始图像'); subplot(2,2,2); imagesc(dsImage); title('欠采样图像'); subplot(2,2,3); imagesc(senseImage); title('SENSE重建'); subplot(2,2,4); imagesc(delta); title('差异图像');

质量评估指标：

结构相似性（SSIM）
峰值信噪比（PSNR）
均方根误差（RMSE）

5. 高级话题与优化方向

5.1 常见问题排查

在实际实现中，可能会遇到以下典型问题：

混叠伪影残留：
- 检查敏感度图准确性
- 验证伪逆计算稳定性
- 考虑添加正则化项
边缘伪影：
- 确保fftshift/ifftshift使用正确
- 检查k空间中心对齐
信噪比下降：
- 优化加速因子选择
- 考虑并行采集校准

5.2 算法优化策略

针对SENSE算法的几种改进方向：

优化方法	实现要点	效果预期
正则化SENSE	在伪逆计算中加入Tikhonov正则化	改善病态问题，抑制噪声
迭代SENSE	使用共轭梯度法等迭代求解	提高重建精度
联合重建	结合压缩感知等稀疏约束	实现更高加速比

% 正则化SENSE示例 lambda = 0.1; % 正则化参数 cHatPinv = (cHat'*cHat + lambda*eye(downSamplingRate)) \ cHat';

5.3 多线圈系统设计考量

线圈数量和布局对SENSE性能有重要影响：

线圈数量：
- 最少需要R个线圈实现R倍加速
- 更多线圈提供冗余，改善重建稳定性
线圈布局：
- 各线圈敏感度分布应有足够差异
- 常用阵列线圈设计：
  - 头部：8-32通道
  - 体部：16-64通道
  - 四肢：4-16通道
几何因子（g-factor）：
- 衡量重建过程中的噪声放大
- 计算公式：g = sqrt(diag(pinv(C'*C)))