做MRI序列仿真这件事我断断续续折腾了快两年。这次的任务是用Matlab实现一个基于FLASH序列的二维布洛赫模拟采集方式是投影k空间也就是常说的radial采样。乍一听有点绕说白了就是逐个体素地计算磁化矢量在不同时刻的状态把射频激励、梯度编码、弛豫衰减全部按时间顺序算出来最后得到k空间原始数据并重建出一张图。这一个模拟跑通之后你对“信号从哪来、ghost伪影和条纹伪影为什么会出现、翻转角和TR怎么影响图像对比度”这些问题的理解会立刻上一个台阶。它非常适合两类人一类是刚接触MRI物理、需要把布洛赫方程和实际序列对上的学生另一类是做序列开发或成像算法验证、想快速确认某个想法是否靠谱的工程师。我属于后者所以这篇文章里的很多取舍都是从“先说结论、再扣细节”的角度出发的。1. 项目思路与总体设计1.1 FLASH序列到底在模拟什么FLASH的全称是Fast Low Angle Shot快速小角度激发。它属于梯度回波序列的一种核心是“小翻转角激励 读出梯度来回切换产生梯度回波”。和自旋回波序列不同FLASH不做180度重聚焦脉冲所以图像对比度由T1加权、T2*加权或质子密度加权决定具体取哪种权重取决于翻转角、TR和TE怎么搭配。布洛赫模拟做的事情就是在每个空间位置初始化一个磁化矢量M然后严格按照序列的时序去更新它。更新内容包括三个方面射频脉冲对M的旋转作用、梯度磁场引起的空间相位累积、以及纵向和横向弛豫导致的指数恢复与衰减。为什么我不直接用解析公式推算信号强度而要在每个像素点上做完整模拟因为解析公式只能给出理想均匀场、无流动、无离共振情况下的稳态信号而一旦加入空间变化的T1、T2、离共振频率或者非笛卡尔k空间轨迹解析路径就断掉了。数值模拟可以保留这些空间信息也能让你亲眼看到每个中间步骤的磁化矢量变化。1.2 为什么选择投影k空间采集传统自旋回波或FLASH序列最常用的是笛卡尔k空间扫描每次TR只填充一行k空间。这种方式的优点是重建简单直接做二维傅里叶逆变换就行。但缺点也明显对运动敏感尤其是相位编码方向的运动会产生经典的重影伪影同时k空间中心区域的采样密度和边缘一致采样效率并不算最优。投影采集则不同每次TR读取的是过中心的一条直径线旋转一个固定角度后继续扫下一条。单看一条线它更像CT的投影数据但把所有角度的线叠在一起中心区域会被反复覆盖天然带有过采样对外围k空间的采样相对稀疏。这种分布在运动伪影和欠采样表现上有天然优势代价是需要处理非笛卡尔k空间数据的重建问题。在模拟层面投影采集还有另一个好处便于验证角度间隔、投影线数、读出点数这三者之间的关系。把线性数从64改成16伪影肉眼可见地变严重对教学和观察来说非常直观。1.3 整体程序框架与设计取舍我把整个程序拆成了四个模块参数定义、体模生成、主循环采集、重建与显示。参数定义模块负责TR、TE、翻转角、FOV、矩阵大小、读出点数、投影线数等所有标量。体模生成模块在一个二维网格上填入不同组织的T1、T2参数模拟一个简单的“水膜背景”结构。主循环模块是整个模拟的核心按每条投影线循环依次执行激励、等待TE、读出采样、填充k空间。最后的重建模块把极坐标分布的k空间数据重新网格化到笛卡尔坐标再做傅里叶逆变换。在设计时我有几个取舍。第一初始磁化矢量设为平衡态[0, 0, 1]不做额外的预稳态处理而是在循环起始额外跑若干次“ dummy”脉冲让系统进入稳态这个后面会细说。第二读出方向采用解析相位累积而非小步长积分减少计算量。第三没有用复杂的高阶网格化算法而是在笛卡尔网格上直接做邻域插值保证核心逻辑清楚、方便比对结果。2. 布洛赫方程数值求解与序列时序设计2.1 布洛赫方程的形式与拆分策略磁化矢量的宏观演化方程很简洁[ \frac{d\mathbf{M}}{dt} \gamma \mathbf{M} \times \mathbf{B} \begin{pmatrix} -M_x/T_2 \ -M_y/T_2 \ (M_0 - M_z)/T_1 \end{pmatrix} ]方程里第一项描述的是磁化矢量在外磁场中的进动第二项是纵向和横向弛豫。数值求解可以直接用小时间步长积分但那样很慢而且难以保证长时间演化后矢量长度不漂移。更实用的做法是把演化过程拆成“旋转”和“弛豫”两个独立步骤。旋转部分在忽略弛豫时有一个解析解激励脉冲相当于绕某个转轴旋转一个角度读出梯度期间的进动相当于绕z轴旋转一个随时间累积的相位。弛豫部分也有解析解横向分量按exp(-t/T2)衰减纵向分量按1-exp(-t/T1)恢复并趋近M0。这种“拆开算”的方式和自旋回波、梯度回波序列的物理过程是对应的因为RF脉冲持续时间通常远小于T1和T2在脉冲期间忽略弛豫不会造成明显误差梯度持续时间也较短弛豫的影响可以通过在每个离散时间点末尾统一做一次衰减来补偿。2.2 离散时间步长怎么定模拟中最重要的两个时间点是TE和读出窗口的持续时间。TE决定了梯度回波信号峰值出现的时刻读出窗口则决定了k空间一条线覆盖的范围。在FLASH序列里时序是先发射一个射频脉冲磁化矢量从纵向翻转到横向平面然后经过一个等待时间达到TE时刻在TE附近打开读出梯度产生梯度回波并在读出窗口内连续采样。在模拟中我把时间离散成三个关键区间激励瞬间用一个旋转矩阵处理持续时间为0。TE等待期持续时间为TE期间同时考虑进动和弛豫。读出期持续时间为读出窗口总长度每隔一个采样间隔记录一次横向磁化分量。摆在面前的问题是读出期内怎样处理自旋的连续进动。我采用的是分点相位累积法不是小步长积分。每一采样间隔Δt内磁化矢量绕z轴转过的相位是γGxΔtx是该体素的空间坐标G是读出梯度的幅度。对每个采样点分别算出累积相位再统一乘上对应的横向弛豫衰减exp(-t/T2*)就能得到该点的信号贡献。步长选择上我建议采样间隔不要超过读出窗口的1/128否则k空间高频部分的相位误差会明显增大。实际实现中采样点数直接对应k空间读出点数比如每线128个点那么Δt就是读出窗口时间除以128精度足够。2.3 空间位置与k空间坐标的映射布洛赫模拟的输入是空间域输出是k空间两者通过相位因子的傅里叶关系关联。一个位于坐标x的自旋在读出梯度Gx作用下经过时间t后累积的相位是γGx∫G(t)dt这个累积量等于2π乘以x和k空间坐标的乘积。在投影采集中读出方向沿着角度θ方向所以第n个采样点的k空间位置是[ k_{x,n} \gamma G_{max} \tau_n \cos\theta ] [ k_{y,n} \gamma G_{max} \tau_n \sin\theta ]其中τ_n是该采样点相对于读出中心的时刻正负对称。这样每条投影线在k空间是一根过中心的直线所有角度的线合起来就是放射状的k空间覆盖。理解了这层关系模拟程序的结构就变得清楚每个空间像素的磁化矢量按相位因子exp(i2πk·r)贡献到对应的k空间采样点所有像素叠加就得到该采样点的复信号。和我平时用的CT模拟不同这里不需要单独做系统矩阵更接近按像素累加的简单方式。3. Matlab实现路径与核心代码3.1 参数定义与体模生成参数模块我习惯集中写在一个结构体里方便后续批量扫描。下面是核心参数的一段代码% sequence parameters seq.FOV 0.2; % field of view, 0.2 m seq.N 64; % reconstruction matrix size seq.TR 20e-3; % repetition time, 20 ms seq.TE 5e-3; % echo time, 5 ms seq.flip 20; % flip angle in degrees seq.Tread 4e-3; % readout duration, 4 ms seq.Nsamples 128; % samples per projection seq.Nproj 128; % number of projections seq.dummy 20; % dummy pulses to reach steady state体模我用了一个简单的同心圆结构外层是T1/T2较短的背景内部放进一个高信号圆盘还有一个离中心的圆点当作可辨认的特征。生成体模时注意T1和T2的单位都是秒和序列参数保持一致。% phantom definition [xg, yg] meshgrid((0:seq.N-1)/seq.N*seq.FOV - seq.FOV/2); disc1 sqrt((xg).^2 yg.^2) 0.04; disc2 sqrt((xg - 0.02).^2 (yg 0.01).^2) 0.01; T1map ones(seq.N)*1.2 0.2*disc1; T2map ones(seq.N)*0.08 0.25*disc1; T2map T2map - 0.06*disc2;这里体模的背景T1设成1.2秒模拟类似脑脊液或缓慢弛豫的组织圆盘区域的T2应变大模拟液体类的高信号结构。实际跑的时候你会发现T2对最终图像对比度的影响极其直接。3.2 主循环逐条投影线更新磁化矢量每条投影线的主循环逻辑如下先计算当前投影角度然后沿着读出方向确定每个采样点的k空间坐标。紧接着对每个空间体素做激励、TE演化、读出采样。kspace zeros(seq.N, seq.N); angles linspace(0, pi, seq.Nproj 1); angles angles(1:end-1); for i 1:seq.Nproj theta angles(i); kdir [cos(theta), sin(theta)]; % RF excitation: rotate Mz into transverse plane Mz_after Mz * cosd(seq.flip); Mxy_after Mz * sind(seq.flip); % complex transverse magnetization % TE evolution: decay and phase accumulation Mxy_after Mxy_after * exp(-seq.TE / T2map); Mz_after Mz_after * exp(-seq.TE / T1map) M0 * (1 - exp(-seq.TE / T1map)); % readout sampling for n 1:seq.Nsamples % time relative to echo center tau (n - seq.Nsamples/2 - 0.5) * seq.Tread / seq.Nsamples; k gamma * Gmax * tau * kdir; % k-space coordinate % phase from gradient phase 2 * pi * (k(1) * xg k(2) * yg); S(n) sum(sum(Mxy_after .* exp(-1i * phase))); % additional T2* decay during readout S(n) S(n) * exp(-abs(tau) / 0.03); end kspace kspace scatter_into_cartesian(kspace, S, theta); % end of TR: longitudinal recovery Mz M0 (Mz_after - M0) * exp(-(seq.TR - seq.TE) / T1map); end这段代码有几个地方需要解释。TE演化期间横向磁化幅度乘以exp(-TE/T2map)纵向磁化向M0恢复用的是指数恢复公式。读出采样阶段每个采样点算一次全域相位叠加相当于对所有像素做一次“离散傅里叶投影”。这里省去了一个细节初始时刻Mxy为0Mz为M0但在第一轮循环前我额外跑了20个dummy脉冲让序列达到稳态。第一次跑的时候不处理采样结果只更新磁化矢量。dummy脉冲数量太少的话前几条投影线的信号会偏弱重建图上能看到一个模糊的中央亮斑。3.3 径向k空间数据的网格化重建投影采集得到的k空间数据分布在一个径向网格上不能直接做二维FFT必须先把它重采样到笛卡尔网格。最简单的做法是找到每个笛卡尔网格点邻近的极坐标采样点用距离加权插值然后填值。function kcart gridding(kspace_radial, sample_points, nx, ny) kcart zeros(nx, ny); for ix 1:nx for iy 1:ny kx (ix - nx/2 - 1) / FOV; ky (iy - ny/2 - 1) / FOV; % find nearest sample point dist sqrt((sample_points.kx - kx).^2 (sample_points.ky - ky).^2); [~, idx] min(dist); kcart(ix, iy) kspace_radial(idx); end end end这种最朴素的近邻插值在小矩阵上表现还可以但会产生一定的网格伪影。实际我推荐稍微改进一下用距离反比加权的多个近邻点或者直接用MATLAB自带的scatteredInterpolant函数做插值代码量更少且结果更平滑。插值完成之后图像重建就一行img fftshift(ifft2(ifftshift(kcart)));这里要注意k空间中心在fftshift前必须位于矩阵中心否则图像会出现半个像素的偏移。我一开始没注意后面发现重建图像的边缘有一圈相位翻转调整数组索引后就好了。3.4 结果对比参数变化的影响模拟的最大价值在于能让你系统性地观察参数变化。作为验证我做了三组对比实验。第一组是把翻转角从10度增加到30度TR保持20ms。结果非常明显大翻转角情况下短T1组织的信号增强图像的T1权重更显著。这和教科书上FLASH序列的T1加权规律完全一致。第二组是改变投影线数。用128条线重建图像清晰锐利降到32条线时图像边缘出现放射性条纹伪影。原因是k空间外围区域采样密度严重不足高频信息缺失这是radial采集欠采样最典型的表现。第三组是把读出点数从128降到64。此时k空间覆盖范围不变但每条线的分辨率下降图像细节变模糊边缘出现振铃。这说明读出采样率必须与k空间最大频率匹配否则高频区域被截断。这三组对比让我对序列参数和图像质量之间的关系有了直观把握比单纯读公式有效得多。4. 常见问题、调试技巧与使用心得4.1 磁化矢量数值发散或越界布洛赫模拟最常见的错误是磁化矢量长度越界表现为Mz大于1或者Mxy的模在数轮迭代后超过1。解决办法首先是检查弛豫公式。纵向弛豫要用M0 (M_initial - M0)exp(-t/T1)不能简化成Mzexp(-t/T1)否则长时间序列的稳态值会错误图像亮度和期望对不上。其次激励旋转矩阵必须保证正交性如果手写旋转矩阵时把cos和sin的位置写错矢量长度会持续扩大最终发散。我建议在每次射频激励后检查一下sqrt(Mx^2 My^2 Mz^2)是否约为1。如果在调试阶段发现偏差超过0.01就说明数值积分或旋转矩阵计算有误。4.2 重建图像中心过暗或过亮径向k空间中心被多条投影线重复覆盖幅度天然偏高。如果重建时不做任何补偿图像中心会出现一个异常的亮斑或者整体对比度被中心低频分量主导。解决方法是施加密度补偿权重。每条线中心的采样点密度高权重应该低外侧采样点密度低权重应该高。最简单的权重计算方式是按半径的倒数设置即w(k) ∝ 1/|k|。在我的程序里我在网格化之前给每个采样点乘上了半径倒数权重重建后的图像均匀性明显改善。4.3 重建图像的条纹伪影和振铃条纹伪影主要来自欠采样这是radial采集的固有现象增加投影线数可以缓解但会增加扫描时间。振铃伪影则多来自读出采样率不足或重建时使用了矩形窗截断。如果你只是想快速验证序列参数而不追求最高图像质量一条实用经验是投影线数至少是矩阵大小的两倍。64×64矩阵对应128条线基本看不到明显条纹32条线下伪影就很扎眼了。这个说法不一定严谨但作为初筛足够实用。4.4 计算效率优化建议二维布洛赫模拟的计算瓶颈在逐条投影线累加信号。按128×128像素、128条投影线、每线128个采样点计算纯Matlab循环需要执行约两千万次复数累加在我日常用的笔记本上跑了大概两三秒体模和数据量再翻倍就明显变慢。优化思路有几个。第一是向量化把空间像素的相位叠加写成矩阵乘法避免for循环内每个像素单独计算。第二是用parfor并行投影线的循环每条线之间没有依赖关系并行起来非常自然。第三是适当减小模拟矩阵或投影线数做初步调试用64×64和64条线跑通后再加大规模。我实际测试下来把内层采样循环向量化之后计算时间几乎缩短了一个数量级。后续如果你要做3D布洛赫模拟这个优化技巧几乎是必须的。4.5 我个人操作中的一些心得这次模拟做完后我最大的体会是序列仿真不是越多代码越好而是要对每一步想清楚为什么。比如我在初始版本里没做dummy脉冲直接进入正式采集结果前几条线信号偏弱重建图像中央区域有不自然的条纹。后来翻资料才意识到FLASH序列的稳态磁化不是一上来就建立好的开头几个TR内纵向磁化会有一个明显的过渡过程。这一点如果你只读信号公式很难注意到但一旦亲手模拟就能看到磁化矢量从初始状态向稳态摆动的整个过程。另外如果你想把这个模拟扩展成更接近实际临床序列的工具还有几个方向可以试。加入离共振效应在TE演化里给不同空间位置设置不同的进动频率可以很好地模拟场不均匀导致的图像变形加入速度项在读出梯度期间让像素相位随时间线性变化可以模拟流动增强或流动伪影改成3D体模则可以把投影采集扩展到球面视角但计算量会成倍提升。最后再分享一个容易踩坑的细节Matlab里角度函数默认使用弧度而实际序列参数一般用度。我在编写射频激励时用cosd/sind在计算k空间坐标时用角度的弧度制转换如果哪个地方混用了sin和sind重建出的图像会变成一团没有规律的噪声。这个坑我踩过一次之后直接把所有角度单位在代码开头固定成弧度并在变量名里标注清楚后续再也没有出过问题。