企业官网建设流程全解析

做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Δt，x是该体素的空间坐标，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为0，Mz为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，重建出的图像会变成一团没有规律的噪声。这个坑我踩过一次之后，直接把所有角度单位在代码开头固定成弧度，并在变量名里标注清楚，后续再也没有出过问题。