做MRI序列仿真的人，搜索“FLASH Matlab”时大概率会翻车。搜出来前几页全是nand flash、spi flash、flash download failed之类的嵌入式内容，很难找到真正想找的MRI里的FLASH序列。这个项目标题里的FLASH，是Fast Low Angle Shot的缩写，也就是快速小角度激发梯度回波序列。标题想做的事情很明确：用Matlab做一次“投影k空间采集”的二维布洛赫模拟，把FLASH序列从RF激发、梯度散相到k空间信号产生的完整物理过程在代码里跑一遍。
这个方向适合谁？两类人最需要。一类是刚进入磁共振成像物理方向的研究生，课程里学了布洛赫方程和k空间理论，但一直没亲手把方程变成可运行的程序，对“中心切片定理”“梯度回波形成”这些概念停留在公式层面。另一类是做序列开发或图像重建算法的工程师，需要快速验证某个序列想法，或者需要一个干净的仿真环境来给重建算法喂数据。这篇博文会给出完整可复现的Matlab实现思路，并解释每一步背后的物理逻辑，以及我在这类仿真里踩过的坑。
1. 先理清概念：FLASH、布洛赫方程与投影k空间
1.1 FLASH核到底是什么
FLASH序列是磁共振成像里最基础的梯度回波序列之一。它的核心特点是小翻转角激发加梯度读出，TR可以做到非常短，所以适合快速成像。相比自旋回波序列用180度重聚焦脉冲来消除主磁场不均匀性，FLASH不重聚焦，而是利用梯度反向来产生回波，因此对T2*敏感，对主磁场不均匀也会比较敏感。
FLASH的信号强度可以写成经典的稳态小翻转角公式：
S = M0 * sin(alpha) * (1 - exp(-TR/T1)) / (1 - cos(alpha) * exp(-TR/T1)) * exp(-TE/T2*)
这里面有几个关键点。首先，翻转角alpha的选择会影响图像对比度和信噪比。当翻转角等于Ernst角时，信号强度达到最大值，Ernst角的计算公式是：
alpha_E = arccos(exp(-TR/T1))
在短TR情况下，Ernst角一般远小于90度。比如TR=12ms、T1=1s时，alpha_E大约在10度到15度之间。这也是为什么FLASH序列通常使用小翻转角，而不是自旋回波那种90度加180度的组合。
标题里的“核”，我理解指的是FLASH序列的核心环节。在仿真里，这个核心就是一条完整的采集链：RF激发把纵向磁化翻转到横向平面，然后读出梯度对自旋进行空间编码，接收线圈采集到的信号就对应k空间的一条线。布洛赫模拟要做的，就是把这个物理过程里每个自旋等色体的磁化矢量演化都算出来，最后叠加得到观测信号。
1.2 布洛赫方程：把磁化矢量运动写成可计算的算符
布洛赫方程是描述宏观磁化矢量在磁场中运动的偏微分方程。在旋转坐标系下，忽略主磁场B0的拉莫尔进动后，方程可以写成：
dMx/dt = gamma * (M x B)_x - Mx/T2 dMy/dt = gamma * (M x B)_y - My/T2 dMz/dt = gamma * (M x B)_z - (Mz - M0)/T1
其中B是实际感受到的磁场，包括梯度场、射频场等。在序列仿真里，事情可以大大简化。梯度场引起的效应是让不同空间位置的自旋感受到不同的磁场，从而产生相位发散；射频场的作用是把磁化矢量绕某个轴旋转一个角度；弛豫则是让横向磁化按T2衰减、纵向磁化按T1恢复。
这种拆分思路在数值实现里非常实用。射频激发可以用旋转矩阵来描述，弛豫可以用指数衰减来描述，梯度散相可以等效成每个体素累积了不同的相位。于是，布洛赫方程就变成了三个独立的算符：RF旋转算符、弛豫算符、进动相位算符。在每个采样间隔内依次施加这三个算符，就可以模拟磁化矢量的完整演化。
二维布洛赫模拟的“二维”，指的是在二维图像空间网格上求解布洛赫方程。假设重建矩阵是128x128，那就在这个二维网格上放置128x128个磁化矢量，每个矢量有三个分量Mx、My、Mz。每个体素有自己的T1、T2和初始磁化强度M0。这种逐体素模拟的好处是，可以直接利用图像域的质子密度图、弛豫参数图来生成具有真实对比度的仿真数据，后续就算不做完整序列模拟，也能拿来做重建算法验证。
1.3 投影k空间采集为什么值得做
投影k空间采集，也叫径向采集。传统的笛卡尔采样是在k空间一条一条水平线扫描，而径向采集是让每条读出线从k空间中心出发，指向某个角度方向。关键的理论基础是中心切片定理：物体二维傅里叶变换沿某个方向的中心切片，等于物体在该方向的投影的一维傅里叶变换。
这个定理意味着，如果我们沿某个角度方向采集k空间数据，那么对这组数据做一维逆傅里叶变换，就得到该方向的投影，也就是Radon变换域的一条数据。把所有角度的投影收集齐，就可以用滤波反投影或者迭代重建算法还原图像。
径向采集在实际扫描里有很多优势。第一，k空间中心被过采样，而中心区域决定了图像的低频信息，所以径向采集对运动伪影不那么敏感，而且中心区域的高信噪比对图像整体质量有帮助。第二，如果投影数不够，径向采样的伪影表现为放射状条纹，而不是笛卡尔采样那种明显的混叠重影，视觉上更容易接受。第三，径向采集天然支持欠采样的加速扫描，是很多压缩感知和低秩重建方法的基础。
所以在仿真里实现投影k空间采集，除了能加深对梯度编码的理解，后续还可以直接在这个框架上加运动和欠采样伪影，用于伪影校正算法的研究。
2. 仿真方案设计与参数选取
2.1 整体架构：从物理过程到程序模块
我在设计这个仿真时，刻意把它拆成四个独立模块：参数初始化、布洛赫算符、采集循环、图像重建。这四个模块各自独立，后续想改序列参数或者加物理效应，只需要动其中一个模块就行。
参数初始化模块负责定义所有仿真的物理量，包括梯度强度、采样点数、TR、TE、翻转角、图像矩阵大小、视野等。布洛赫算符模块实现射频激发、弛豫和进动三个基本操作。采集循环是主程序，按投影角度循环执行：先对磁化矢量施加RF激发，然后在每个采样时刻计算当前k空间坐标，对各体素的相位累积，最后把整个图像域的横向磁化矢量叠加，得到一个复数值，这就是该采样点的k空间信号。重建模块负责把采集到的径向k空间数据变成一张图像，这里可以直接用Matlab的iradon函数。
这种模块化设计还有一个好处，就是方便查错。如果重建图像不对，可以单独检查采集循环输出的k空间数据是否符合中心切片定理，或者单独检查布洛赫算符在无梯度情况下的信号衰减曲线是否和理论公式匹配。
2.2 参数表与计算过程
下面给出我实际使用的参考参数，这些参数大致模拟一个常规的头部FLASH扫描，FOV是256mm，矩阵尺寸128x128。
| 参数 | 符号 | 数值 | 说明 |
|---|---|---|---|
| 视野 | FOV | 0.256 m | 成像区域尺寸 |
| 图像矩阵 | N | 128 | 重建矩阵大小 |
| 翻转角 | alpha | 15° | 接近Ernst角 |
| 重复时间 | TR | 12 ms | 相邻RF激发间隔 |
| 回波时间 | TE | 5 ms | 从RF到回波中心 |
| 投影数 | nProj | 180 | 180度范围内采集条数 |
| 每投影采样点 | samplesPerProj | 128 | 一条径向上的采样点数 |
| T1（灰质） | T1 | 1.0 s | 纵向弛豫时间 |
| T2（灰质） | T2 | 80 ms | 横向弛豫时间 |
k空间的分辨率由FOV和矩阵大小共同决定。采样间隔delta_k = 1/FOV，最大k值kmax = N / (2 * FOV)。代入数值，delta_k约等于3.9 1/m，kmax约等于250 1/m。这个kmax对应的空间分辨率正好是FOV/N，也就是2mm。
关于投影数，极坐标采样的Nyquist准则要求总投影数至少是pi/2 * N的量级，128矩阵至少需要201条投影，180条算是一个兼顾时间与质量的折中值。如果投影数少于这个值，重建图像会出现可见的放射状条纹伪影，后面调试部分会演示这个问题。
Ernst角这里算一下。把TR=0.012s、T1=1.0s代入alpha_E = arccos(exp(-TR/T1))，得到约14.06度。所以15度的翻转角是一个很合理的选择，信噪比高而且T1权重适中。
2.3 仿真对象的设定
仿真对象我直接用Matlab自带的Shepp-Logan幻影，这个幻影在图像重建领域是标准测试对象，内部有多个椭圆结构，可以用来验证空间分辨率和对比度。把幻影当作质子密度图M0，再给不同区域赋予不同的T1和T2值，就能得到一个有真实对比度的仿真对象。
注意一个容易困惑的点：M0是空间变量，表示每个体素的自旋密度，也就是平衡状态下的纵向磁化强度。在稳态FLASH中，如果只模拟单个TR，我们可以直接从M0出发，施加alpha翻转后，纵向磁化变成M0cos(alpha)，横向磁化变成M0sin(alpha)。如果想更精确地模拟连续多TR的稳态，还需要在多TR循环里加入T1恢复和再次激发，这部分扩展在第五节单独讨论。
3. Matlab核心实现：从零构建布洛赫模拟
3.1 初始化与幻影准备
第一步先把参数和仿真对象准备好。我用的是伪代码风格，实际运行要放进一个脚本里。
%% 参数初始化 N = 128; % 矩阵尺寸 FOV = 0.256; % 视野，单位m dx = FOV / N; x = linspace(-FOV/2, FOV/2 - dx, N); [xg, yg] = meshgrid(x, x); T1 = 1.0; % 灰质T1，单位s T2 = 0.08; % 灰质T2，单位s M0 = double(phantom(N)); % 质子密度图，Shepp-Logan幻影 % 序列参数 TR = 0.012; % 重复时间 12ms TE = 0.005; % 回波时间 5ms alpha = deg2rad(15); % 翻转角 15度 nProj = 180; % 投影数 samplesPerProj = N; % 每条投影采样点数 readoutTime = 2 * TE; % 读出总时长，这里简化取10ms dt = readoutTime / (samplesPerProj - 1); % 采样间隔这里读出总时长我取了2*TE，也就是说回波中心落在读出窗口中间，对应k空间原点。这样每条投影线从-kmax扫到+kmax，正好是完整的k空间直径，重建时用普通的1D逆傅里叶变换就能得到正确的投影剖面。
3.2 布洛赫算符的代码实现
布洛赫算符是核心中的核心。RF激发我用旋转矩阵实现，弛豫用指数衰减实现，进动用相位旋转实现。三个函数可以单独写，也可以直接内联在采集循环里。
射频激发绕y轴翻转alpha角度的代码：
function M = apply_rf(M, alpha) % M是三维数组(N,N,3)，分别对应Mx,My,Mz My = M(:, :, 2); Mz = M(:, :, 3); M(:, :, 2) = My .* cos(alpha) - Mz .* sin(alpha); M(:, :, 3) = My .* sin(alpha) + Mz .* cos(alpha); end弛豫算符在每个时间步内更新横向和纵向磁化：
function M = apply_relaxation(M, M0, T1, T2, dt) M(:, :, 1) = M(:, :, 1) .* exp(-dt / T2); M(:, :, 2) = M(:, :, 2) .* exp(-dt / T2); M(:, :, 3) = M0 - (M0 - M(:, :, 3)) .* exp(-dt / T1); end进动算符把每个体素的横向磁化矢量绕z轴旋转一个由梯度累积决定的相位。在读出梯度存在时，不同位置的体素相位变化率不同：
function M = apply_precession(M, xg, yg, Gx, Gy, dt) gamma_bar = 42.58e6; % 旋磁比，单位Hz/T phase = 2 * pi * gamma_bar * (Gx * xg + Gy * yg) * dt; Mxy = M(:, :, 1) + 1j * M(:, :, 2); Mxy = Mxy .* exp(-1j * phase); M(:, :, 1) = real(Mxy); M(:, :, 2) = imag(Mxy); end这里有一个经验性选择：我一个时间步内先做进动再做弛豫。对于比T1和T2都小很多的时间步长，这种分裂算符方法的误差可以忽略。如果采样点数是128，读出时长10ms，那么时间步长大约78微秒，远小于80ms的T2，精度没问题。
3.3 投影数据采集主循环
采集循环是仿真的主干。整体逻辑是：对每个投影角度，先初始化磁化矢量到平衡态，施加RF激发，然后沿径向量化方向逐步模拟读出过程。
这里我采用了完整的布洛赫演化方式，也就是每个采样时刻都对整个图像域的磁化矢量做一次进动和弛豫更新，然后把横向磁化矢量在图像域内求和，作为该时刻的k空间信号。这种做法的好处是物理过程透明，后续要加扩散、流动等效应也方便。
%% 初始化径向k空间数据矩阵 kSpaceRadial = zeros(nProj, samplesPerProj); % 旋磁比 gamma_bar = 42.58e6; % 读出梯度幅度，由kmax和读出时长决定 kmax = N / (2 * FOV); G0 = kmax / (gamma_bar * (readoutTime / 2)); % 主循环 for p = 1:nProj theta = (p - 1) * pi / nProj; Gx = G0 * cos(theta); Gy = G0 * sin(theta); % 初始化磁化矢量：纵向为M0，横向为0 M = zeros(N, N, 3); M(:, :, 3) = M0; % RF激发：绕y轴翻转alpha M = apply_rf(M, alpha); % 读出采样 for s = 1:samplesPerProj % 记录当前采样点的k空间信号（横向磁化矢量的总和） Mxy = M(:, :, 1) + 1j * M(:, :, 2); kSpaceRadial(p, s) = sum(Mxy, 'all'); % 施加进动和弛豫，进入下一个采样时刻 if s < samplesPerProj M = apply_precession(M, xg, yg, Gx, Gy, dt); M = apply_relaxation(M, M0, T1, T2, dt); end end end这段代码实际上是在逐个时刻模拟FID信号。在读出梯度的作用下，不同位置的体素相位发散，横向磁化矢量互相抵消，信号逐步衰减。但有意思的地方在于，如果后续对Mxy求和前不施加任何逆相位，那么得到的kSpaceRadial实际上是以梯度累积为k坐标的k空间数据。这正是中心切片定理的直接体现：每个投影角度的信号序列就是该角度方向的k空间剖面。
有一点需要注意：信号求和得到的是一个复数，里面包含了该方向k空间数据的幅度和相位。在后续重建中，这些数据要当作复数处理，不能只取实部或模值。
3.4 径向k空间重建与效果验证
采集到径向k空间数据后，重建思路分两步。第一步是把这个极坐标采样的k空间数据转换成sinogram，也就是Radon变换域的投影数据。对每条k空间直径线做一维逆傅里叶变换，得到的是该方向的一维投影剖面。第二步是用iradon函数做滤波反投影重建。
%% 将径向k空间数据转换为sinogram % 每条投影线从-kmax到+kmax，需要ifftshift调整FFT原点 kSpaceShifted = ifftshift(kSpaceRadial, 2); proj = real(ifft(kSpaceShifted, [], 2)); proj = proj * sqrt(2) * kmax; % 幅度校正因子 %% iradon重建 theta_deg = linspace(0, 180 - 180/nProj, nProj); img = iradon(proj, theta_deg, 'linear', 'Ram-Lak', 1, N); %% 显示结果 figure; subplot(1,2,1); imshow(M0, []); title('原始幻影'); subplot(1,2,2); imshow(img, []); title('径向重建结果');重建结果应该和原始Shepp-Logan幻影基本一致，只在投影数不足时出现轻微的放射状条纹。如果出现图像上下颠倒或者左右翻转，需要检查theta_deg的方向定义是否和采集循环里的角度定义一致。我在采集循环里用theta = (p-1)*pi/nProj，对应逆时针从x轴正方向开始，而iradon默认也是逆时针从x轴正方向开始，两者是匹配的。
重建结果的验证可以从两个维度看。第一是空间分辨率，128矩阵重建后应该能清晰分辨Shepp-Logan幻影里的几个小椭圆结构。第二是对比度，因为所有体素T1和T2都一样，重建图像对比度应该和M0的对比度一致。如果想看T1加权效果，可以把不同组织区域设置不同的T1值，然后观察纵向磁化恢复速度对信号的影响。
4. 调试记录与避坑指南
4.1 重建图像出现条纹和拖影
最常见的重建问题就是放射状条纹，尤其在图像边缘特别明显。这个伪影几乎都是投影数不足造成的。径向采集的欠采样伪影是放射状的，条纹数量和投影数呈反比。我试过64条投影重建128矩阵，图像中间还可以，边缘全是放射状亮斑，完全没法看。
解决方法是增加投影数，至少满足nProj >= pi/2 * N这个经验公式。128矩阵用201条投影，如果还嫌条纹明显就加到256条。另外，重建时一定要用滤波函数，我默认用'Ram-Lak'，也就是斜坡滤波，这是标准的滤波反投影选择。如果不用滤波，图像会非常模糊，而且本质上是低通滤波的效果。
还有一种隐藏的条纹来源是k空间数据没有做ifftshift。Matlab的fft假设空间原点在数组第一个元素，而径向采样的k空间原点在数组中心。如果不做ifftshift直接逆FFT，得到的投影剖面是错位的，重建图像会有严重的混叠伪影。
4.2 信号幅度异常或全为零
如果发现采集到的kSpaceRadial全是零或者小到接近机器精度，大概率是相位计算出了问题。最常见的原因是角度变量进了三角函数但没有转弧度，比如把theta_deg直接丢给cos和sin。Matlab的cos和sin默认用弧度，如果传入角度值，相位就会完全错误，所有体素的相位互相干扰，信号被抵消成接近零。
另一个坑是旋磁比的单位搞混。我用的是gamma_bar = 42.58e6 Hz/T，也就是以Hz为单位的旋磁比，这样和k空间的公式能匹配上。如果用rad/s单位，公式里2pi的处理就要相应调整。我在写初版代码时吃过这个亏，导致信号频率差了一个2pi的因子。
还有一个物理上的坑：单个TR内，如果RF激发后没有立刻开始采集，而是等到TE时刻才开始，那在RF和TE之间有一段自由进动时间，横向磁化会散相。我在上面的简化代码里是RF后立刻开始读出，回波中心落在读出窗口中间，这有点类似自旋回波的时间组织方式。更标准的FLASH模拟应该是在RF和读出之间插入prephasing梯度的时间，让k空间轨迹从-kmax开始，到中心时信号达到最大值。如果这段prephase时间被省略，信号衰减特征会不一样，重建结果虽然还能出来，但对比度和标准FLASH会有偏差。
4.3 仿真慢到怀疑人生
二维布洛赫模拟最直接的性能瓶颈是双层循环里反复处理128x128x3的数组。180条投影乘128个采样点，就是23040次阵列更新，每次都要做矩阵乘法和指数运算，跑起来确实磨人。
实际优化可以从三个方向入手。第一个是向量化，把内层采样点的循环尽量用矩阵运算代替，比如先计算整个采样时刻的相位数组，再一次应用到所有体素，而不是每个体素单独更新。第二个是降采样验证，如果只是验证思路，先用32x32矩阵、60条投影把流程跑通，确认无误后再上完整参数。第三个是用提前计算好的查表来代替重复的三角函数计算。比如每条投影的cos(theta)和sin(theta)可以提前算好放在数组里，而不是在每次循环里重新调用cos和sin。
我实测下来，32x32矩阵加60条投影，整个流程几秒就能跑完，非常适合调试。128x128加180条投影，在普通笔记本上大约需要几分钟，这个速度在做科研复现时是完全可以接受的。
4.4 快速验证正确性的三个小技巧
在跑完整仿真之前，我建议先做三个快速自检，能省下大量排查时间。
第一个技巧是无梯度验证。把读出梯度设成零，那么所有体素同相位进动，横向磁化矢量的叠加信号应该是一条按T2衰减的自由感应衰减曲线。把时间信号画出来和理论曲线exp(-t/T2)对比，如果吻合说明弛豫算符和信号求和的量纲都没问题。
第二个技巧是单像素验证。把M0设成一个只在图像中心有值、其余全为零的矩阵，这时采集到的信号应该是单一频率的复指数，幅度按T2衰减。经过重建后应该是一个点扩散函数，可以通过这个点的形状看重建链条是否完整。
第三个技巧是直接验证中心切片定理。把原始幻影做二维傅里叶变换，然后沿着某个角度截取中心线，与采集循环里对应角度输出的k空间数据对比。两者应该完全一致（除了缩放因子和数值精度误差）。如果这里对不上，说明采集循环或k空间采样方式有问题。
5. 从单次采集到多TR稳态：扩展方向与个人体会
5.1 多TR稳态布洛赫模拟的改造思路
上面的实现默认每个TR开始时纵向磁化都恢复到M0，这其实是一个很强的简化假设。真实FLASH序列TR很短，纵向磁化在一个TR内不可能完全恢复，经过几个TR后会进入稳态。稳态时的纵向磁化强度不再是M0，而是：
Mz_ss = M0 * (1 - exp(-TR/T1)) / (1 - cos(alpha) * exp(-TR/T1))
要模拟多TR的稳态过程，只需要在外层加一个TR循环，每次RF激发前把当前纵向磁化代入，RF激发后让磁化经历T1恢复和T2衰减，然后进入下一个TR。这个改造的重点是处理好状态变量在不同TR之间的传递，避免每次重新初始化副作用导致稳态无法建立。
通常模拟连续几十个TR就能收敛到稳态，收敛速度由T1和TR的相对大小决定。这个多TR框架做好之后，就能模拟更真实的FLASH序列行为，比如不同翻转角的信号差异、进入稳态之前的瞬态振荡效应，以及反转恢复准备脉冲对对比度的影响。
5.2 进一步加料：伪影、流动与RF非理想性
当前的模型里，所有体素的T1、T2是均匀的，RF是理想的硬脉冲，梯度是完美的梯形脉冲。实际场景中这些假设都太理想。想做伪影校正研究的话，可以在仿真里加入运动项，让M0在读取过程中发生位移，就能模拟出运动伪影。可以在体素维度加随机频率偏移，模拟主磁场不均匀性，这时的径向采集对场不均匀性的表现会和笛卡尔采样很不一样，是很好的研究方向。
RF非理想性也有实际意义。真实序列的RF脉冲有幅值误差和相位误差，接收通道也有灵敏度图。把这些因素加进布洛赫算符里，就能研究翻转角误差对不同组织对比度的影响，对定量成像研究很有帮助。
5.3 给入门者的一句话
我做磁共振物理仿真这几年，最大的体会是：代码本身不难，难的是把每个公式和代码行对应起来。布洛赫方程、k空间、中心切片定理，这些概念刚学的时候觉得抽象，但只要亲手写一遍仿真，把信号和图像跑出来，很多模糊的理解自然就通了。如果这篇里的代码能帮你在某个深夜突然想通一个困扰很久的物理细节，那这篇博文就没白写。从最简单的单TR投影模拟开始，一点点把磁共振序列的世界打开，这个过程本身就很有意思。
