做磁共振成像MRI序列开发的同事应该都有体会一个新想法直接搬到扫描仪上做实验调试周期长、机时成本高信号链里混入的各种系统误差还会干扰你对序列本身的判断。所以业界通行的做法是先在Matlab里做布洛赫模拟Bloch simulation把脉冲序列在“理想世界”里的行为跑透再决定要不要上机。这个项目做的就是这件事用二维布洛赫模拟实现FLASH序列的投影k空间radial k-space采集把序列参数、体模定义、k空间填充轨迹、重建图像全部放进一个Matlab工程里跑通还附带可直接复用的源码工程。先说清楚这个工程适合谁。如果你是刚接触MRI物理的硕士生想搞明白FLASH序列为什么用低翻转角、k空间为什么可以沿径向填充这份代码能让你一步步看到从自旋磁化矢量到k空间信号再到图像的完整链路。如果你是在做序列开发或图像重建的研究人员这个模拟工程可以用作新序列的验证基线——在仿真里确认正确性再动手改真机序列。核心价值一句话把“FLASH 径向k空间采样 图像重建”这条链路以参数完全可控、结果可反复检查的方式完整复现在电脑上。1. 项目背景与序列设计思路1.1 FLASH序列为什么能“快”FLASH全称Fast Low Angle Shot是梯度回波GRE家族的经典序列。它的设计目标用一句话概括在尽量短的重复时间TR内拿到可用的图像对比度。为此它做了两个关键取舍一是用梯度反转代替180°重聚脉冲来形成回波省掉了自旋回波序列里耗时的重聚脉冲二是把激发翻转角压得很低通常5°~30°让纵向磁化在每次激发后不会彻底耗尽稳态下仍有充裕的纵向磁化供下次激发使用。这两个特点叠加FLASH才能在几十毫秒甚至几毫秒的TR下连续成像成为临床3D扰相梯度回波如腹部T1增强扫描、颅脑高分辨T1的基本盘。你可能会问翻转角低不会让信号变弱吗这里的关键是稳态。短TR下若用90°翻转角一次激发后纵向磁化几乎清零TR余量不够它恢复下一拍的信号反而很弱倒是小角度激发让每次“消耗”的纵向磁化不多稳态时能持续产生稳定的横向磁化。这就是恩斯特角Ernst angle的理论来源后面我会专门算一笔账。1.2 投影k空间采集与笛卡尔采样的取舍常规MRI的k空间是逐行采的每次TR采一条平行线效率高、重建简单一维FFT即可。但笛卡尔采集有个明显短板对运动特别敏感扫描中一旦发生移动对应的相位编码线直接错乱伪影沿着相位编码方向蔓延整幅图几乎报废。径向投影采集则换了个思路每次TR沿一条穿过k空间中心、角度旋转的直线采样很像CT里的Radon变换测量。好处有三个一是k空间中心被反复过采信号能量最集中的区域天然得到强化对运动造成的相位误差有更强的容忍度二是运动伪影从“整图错位”变成相对柔和的放射状条纹信息损失更局部三是配合黄金角角度步进可以随时截断数据做动态重建时间分辨率灵活。代价也同样明显重建不能直接FFT需要网格化gridding或滤波反投影spoke数不足时图像上会出现明显的放射状条纹伪影。这个项目把radial读取的每个细节——角度步进、读出梯度持续时间、k空间填充坐标——全部纳入布洛赫模拟目的就是把这些权衡看得清清楚楚。2. 布洛赫方程与FLASH信号的定量关系要模拟MRI信号绕不开布洛赫方程。MRI信号本质上是宏观磁化矢量M(Mx, My, Mz)绕主磁场B0进动时在接收线圈中感应出的电动势布洛赫方程就是描述这个矢量时间演化的动力学方程dMx/dt γ(M × B)x − Mx/T2 dMy/dt γ(M × B)y − My/T2 dMz/dt γ(M × B)z − (Mz − M0)/T1其中γ是旋磁比B是总磁场主磁场射频场梯度场T1和T2分别是纵向和横向弛豫时间。这个方程组看着不复杂但射频脉冲、梯度、弛豫同时作用时很少有人手解数值仿真成了标配。2.1 把布洛赫方程拆开看布洛赫方程可以拆成三项理解。进动项γ(M×B)描述自旋绕磁场转动这是信号产生的来源——横向磁化绕主磁场进动在线圈中感应出可测量的复数信号梯度场能实现空间编码正是因为它改变不同位置的进动频率让磁化在读出时间内累积出与位置成正比的相位。弛豫项−Mx/T2、−My/T2描述横向磁化散失T2越短信号衰得越快纵向项−(Mz−M0)/T1描述纵向磁化向平衡值恢复T1越短恢复越快。在FLASH的一个TR周期里这三类过程按时间顺序重复上演射频脉冲把纵向磁化“扳倒”成横向磁化TE等待期间横向磁化一边进动一边按T2*衰减TR余量内纵向磁化按T1恢复。仿真要做的事就是把这三种算子按序列时间轴依次作用到每个体素的磁化矢量上。2.2 稳态信号公式与恩斯特角对扰相FLASHTR末把横向残留磁化spoiler掉来说稳态下的横向磁化有很漂亮的解析表达式Mxy M0 · sin(α) · (1 − E1) / (1 − cos(α)·E1) · E2其中α是翻转角E1 exp(−TR/T1)E2 exp(−TE/T2*)。当TR远小于T1时使信号最大化所需翻转角近似为α(Ernst) ≈ sqrt(2·TR/T1)假设T1800 ms、TR10 ms估算出来翻转角约5.8°考虑到T2*加权和组织差异实际常取8°~15°。模拟时我习惯先用这个公式算一个理论最优角再用布洛赫模拟扫一组翻转角验证二者趋势一致说明模拟实现没有原则性错误。这里多说一句解析公式假设了理想扰相、无偏共振、读出瞬间完成真实序列行为与它有偏差是正常的模拟的价值恰恰在于把这些“非理想”因素逐个加回去。2.3 为什么有了解析公式还要数值模拟既然稳态有解析解为什么还要费劲做布洛赫模拟因为解析解做了太多理想化假设忽略偏共振失相位化学位移、主磁场不均、忽略读出采样点之间的横向磁化衰减、忽略序列启动阶段的瞬态过程。真要设计新序列——比如改扰相梯度、加脂肪抑制准备脉冲、调读出方向——这些假设全都站不住必须回到数值解。这个项目的做法是逐体素更新磁化矢量在各体素上按时间顺序施加射频激发、梯度相位累积、T1/T2弛豫算子完整还原磁化矢量的“旅程”。虽然比套公式慢但能看到每一个中间状态这才是理解序列行为最直观的路径。3. Matlab仿真框架设计与实现这一节是实操核心。我把仿真架构分成五层参数配置层、体模定义层、序列事件编排层、布洛赫求解层、k空间填充与重建层。参数配置层集中存放TR、TE、翻转角、读出点数、spoke数、FOV、网格分辨率体模定义层生成一个有真实T1/T2/质子密度分布的二维目标序列事件编排层把“激发—等待TE—读出采样—恢复等待”循环按时间轴排成事件清单布洛赫求解层在每个体素上按事件清单推进磁化矢量最后把每个读出采样点的信号按坐标放入k空间。3.1 总体架构与数据流数据流一句话概括体模每个网格点一个磁化矢量→ 序列事件驱动磁化矢量演化 → 读出窗口内每隔一个采样间隔Δt采一个复数信号样本 → 按投影角度填入k空间 → 所有spoke采完后重建成图像。为了让代码既清晰又算得动两个细节我建议这样处理。第一每个体素看成一个独立的自旋等色群isochromat体素内的T1/T2和偏共振频率取该格子的中心值。加密体素数量可以减少磁化矢量“离散化”带来的误差但计算量按网格尺寸平方增长64×64和128×128的差别非常明显起步阶段用小网格调通逻辑再放大。第二梯度编码不通过逐时刻积分整个梯度波形实现而是直接按“k空间位置 读出采样点 × 对应方向梯度矩”计算相位。这与MR物理中k空间的定义完全一致也大幅节约计算量——毕竟真正的布洛赫积分要细到微秒级时间步全网格跑一遍代价很高。3.2 体模选择与参数设置仿真实验体模我建议先用Shepp-Logan数字体模它包含多个不同T1/T2/质子密度的区域足以观察对比度和伪影行为。Matlab自带phantom函数生成几何形状自己做一层T1/T2映射表即可。要注意体模网格分辨率必须与k空间读出点数匹配图像矩阵若为128×128体模也定义在128×128上读出点数至少取256沿直径方向过采样一倍spoke数按奈奎斯特准则至少需要约π/2倍矩阵点数也就是大约201条。实践里我常取的组合是128×128体模、256读出点、256条spoke既有明显的径向采集特征又不至于让计算量失控。我常用的基准参数如下表参数数值说明主磁场场强1.5 T旋磁比按γ42.58 MHz/TTR / TE10 ms / 3 ms短TR快速成像典型值翻转角10°接近恩斯特角附近采样带宽±31.25 kHz决定采样间隔Δt读出采样点数256每条spoke样本数投影spoke数256按π/2倍矩阵尺寸过采样图像矩阵128×128重建网格预热TR数50让磁化矢量进入稳态3.3 序列事件怎么编排序列的时间轴是布洛赫模拟的“剧本”。一个TR循环包含三段射频激发段、回波形成与读出段、驰豫恢复段。每段对应不同的更新算子。激发段射频脉冲把磁化矢量绕某个轴旋转旋转角就是翻转角。矩形脉冲的射频持续约0.5~1 ms仿真里可以简化为瞬时旋转相比TR内的弛豫时间很短误差可接受。旋转矩阵的具体形式是绕xy平面内相位为φ的轴转α角用标准三维旋转矩阵即可实现。读出段这一段的操作最密集。横向磁化绕z轴按“偏共振频率梯度场贡献”的总频率旋转同时按T2*衰减纵向磁化按T1恢复。读出窗口内每个采样间隔Δt采一次信号信号值就是当前横向磁化矢量的复数表示Mx i·My。梯度场贡献在投影方向上的体现是体素位置沿投影方向的坐标乘以梯度矩得到累积相位这正是k空间编码的本质。恢复段没有射频和梯度只有T1/T2弛豫。对扰相FLASH这一段的末尾还要模拟扰相操作把残余横向磁化清零确保下一个TR开始时横向磁化从零开始。3.4 核心Matlab代码实现下面给出几个最关键的函数片段。第一个是磁化矢量的RF旋转function M applyRF(M, alpha, phi) % alpha: 翻转角单位rad % phi: RF脉冲的初始相位单位rad % 标准三维旋转绕xy平面内方位角为phi的轴旋转alpha c cos(phi); s sin(phi); ca cos(alpha); sa sin(alpha); R [ca c^2*(1-ca), s*(1-ca), s*sa; s*(1-ca), ca s^2*(1-ca), -c*sa; -s*sa, c*sa, ca]; M R * M; end第二个是包含弛豫与偏共振的推进函数自由进动和弛豫可以合并成闭式解function M relax(M, dt, T1, T2star, dfreq) % dt: 持续时间 % T1,T2star: 纵向/有效横向弛豫常数 % dfreq: 偏共振频率(Hz) E1 exp(-dt/T1); E2 exp(-dt/T2star); phi 2*pi*dfreq*dt; % 偏共振导致的相位累积 Mc E2 * (M(1)1i*M(2)) * exp(1i*phi); M(1) real(Mc); M(2) imag(Mc); M(3) M0 - E1*(M0 - M(3)); % M0为平衡磁化 end第三个是主循环展示每个TR里体素磁化矢量的演进和信号采集for iSpoke 1:Nspokes theta (iSpoke-1) * pi / Nspokes; % 投影角度均匀覆盖180° Gx G0 * cos(theta); Gy G0 * sin(theta); % 5.1 射频激发 M applyRF(M, alpha, 0); % 5.2 TE等待自由进动 弛豫 M relax(M, TE, T1, T2star, dfreq); % 5.3 读出采样每个采样点对应一个k空间坐标 for iS 1:Nread kpos (iS - Nread/2 - 0.5) * dk; % 梯度编码相位位置在投影方向的贡献 phase_grad 2*pi * kpos * (Gx*X Gy*Y); sig (M(1) 1i*M(2)) .* exp(-1i*phase_grad); signal(iSpoke, iS) sum(sig(:)); % 所有体素相干叠加 end % 5.4 恢复段TR剩余时间 理想扰相 M relax(M, TR - TE - readTime, T1, T2, 0); M(1) 0; M(2) 0; % 理想扰相清空横向磁化 end注意上面这个写法是教学化的双层循环变量名里X、Y是meshgrid之后的位置坐标数组relax函数要支持对整数组体素同时运算。实际工程一定要矢量化把体素展开成一维数组所有体素同时推进信号通过矩阵乘法累加。我实测128×128体模、256条spoke、50个预热TRMatlab向量化后单次完整模拟大约1~2分钟完全可以接受。3.5 预热TR与稳态收敛这是新手最容易忽略的一步直接从M0开始模拟并记录信号得到的是瞬态信号不是FLASH真正工作的稳态信号。磁化矢量要经过多个TR才能进入动态平衡。处理办法是在正式采数前先跑20~50个“预热TR”只推进磁化矢量不记录信号。预热TR数量取决于T1/TR比值。以T1800 ms、TR10 ms为例20个TR后纵向磁化恢复到平衡值的1−exp(−20×10/800)≈22%还差得远想要98%的恢复大约需要800/10×ln(50)≈31个TR。所以我统一取50个预热TR代价只是增加约1/4的总计算量换来的稳态准确性对仿真可信度很重要。4. 仿真结果怎么看从k空间到图像仿真做完不能只在Matlab工作区里存几个数组要把过程“看见”。我通常做三件事画k空间填充轨迹、看每条spoke的原始信号幅度、重建成图像并做参数对比。4.1 k空间填充轨迹验证早期最常见的bug之一就是k空间坐标算错导致重建图像旋转或折叠。验证方法很简单把每条spoke的(kx, ky)坐标画成散点图正常应是以原点为中心、辐射状均匀覆盖半径范围的多段直线。角度步长为均匀π/Nspokes时相邻spoke的角度间隔应完全相等。如果看到坐标点乱窜、或某一方向的“扇区”明显稀疏问题多半出在梯度矩计算或采样索引偏移上。这一关必须在重建前过掉否则后面所有图像结论都不可靠。4.2 径向k空间重建从数据到图像径向k空间不能直接做笛卡尔FFT。常用两条路一是网格化重建gridding把非笛卡尔坐标上的样本用卷积核凯撒窗等插值到笛卡尔网格再做FFT二是滤波反投影FBP把每条spoke数据先做一维傅里叶变换再沿对应角度反投影到图像域本质是逆Radon变换。Matlab里做FBP可以直接用radon和iradon函数配合处理做gridding则可以用较新版本内置的nufft% ktraj为每条spoke每个采样点的归一化k空间坐标Nx×Ny为重建网格 ktraj buildKtraj(Nspokes, Nread); % 按投影角度和读出位置生成 data signal(:); % 采集信号展平 img nufft2(data, ktraj, [Nx Ny]); % 2D NUFFT重建用NUFFT重建时要注意坐标归一化范围不同版本的nufft函数对输入坐标域的定义不完全一致建议先用一个已知解析信号比如点扩散函数仿真验证重建正确性再正式使用。4.3 参数对比翻转角、spoke数、T2*的直观效果这部分是该模拟最有教学价值的地方。固定其他参数把翻转角从5°、10°、20°、35°扫一遍重建图的平均信号强度变化应当符合恩斯特角曲线把spoke数从256降到64再降到32图像会明显出现放射状条纹伪影把T2从30 ms改到10 ms读出时间内信号衰减加剧图像会变模糊、边缘变钝这就是典型的T2模糊效应。这三组对比都以子图形式保存既方便写实验报告也方便后续调参时做对比参照。5. 调试中踩过的坑与实战经验仿真不难难在让仿真结果可信。以下是我在这个项目里实际踩过、也常见同行掉进去的坑。5.1 单位制混乱导致所有参数失效梯度矩、k空间坐标、频率这三个量之间单位极易搞混。一个稳妥的口令是全部用国际单位制m、s、T、Hz旋磁比按γ42.58×10^6 Hz/T注意是Hz/T若要用rad单位需要乘2π。读出采样间隔Δt1/带宽读出梯度幅值Gread带宽/(γ·FOV)k空间步长Δk1/FOV。所有推导都建立在kγ∫Gdt这个基本关系上。一旦出现GHz·m级别的奇葩数量级先查单位不要着急查逻辑。我在一个早期版本里把旋磁比写成42.58×10^6 rad/s/T结果所有梯度矩偏小2π倍图像整体错位查了整整半天才定位至今记忆犹新。5.2 稳态没收敛就采数图像出现“渐变”伪影预热TR不足时前面若干条spoke的横向磁化比后面小k空间中央的信号权重异常重建图像呈现亮度不均或奇怪的渐变伪影。排查方法是画出每条spoke的中心信号幅度k空间原点样本随spoke序号的变化正常应是一条先快速下降后趋平的曲线如果看到持续波动把预热数加倍再看。5.3 扰相不彻底导致T1/T2混合对比FLASH是扰相GRE序列前提是每个TR末把横向磁化残留清掉。如果扰相梯度矩不足、或代码里压根没实现扰相步骤稳态信号里会混入相干横向磁化的贡献对比度就会偏向bSSFP如FIESTA/TrueFISP那种T1/T2混合。判断方法很直接对多个偏共振频率做模拟看信号是否随频率振荡——扰相不彻底时信号会出现明显振荡。仿真里直接把Mx、My置零是“理想扰相”真机需要用梯度扰相加射频相位循环RF spoiling来逼近这个落差本身也是很好的研究对象。5.4 读出点数与FOV不匹配导致混叠折返radial采集的每条spoke本质上是沿某方向的1D频率编码。若读出点数乘以Δk小于需要覆盖的k空间半径约1/(2×体素尺寸)高频信息就缺失图像边缘会出现混叠折返。经验公式Nread至少等于图像矩阵尺寸才能完整覆盖直径方向想要过采样余量就用1.5~2倍。我在工程注释里特意写了这条经验公式因为这是改参数时最容易忽略的连锁反应。5.5 偏共振与T2*衰减缺失结果过于“理想”不少初版模拟只做RF激发和梯度相位忽略偏共振和T2*。这样仿出来的图像特别干净漂亮但一上真机就对不上。我的建议是从一开始就在模拟里加入两个“现实因子”一是每个体素的偏共振频率dfreq可由主磁场不均匀、化学位移、脂肪成分等引起在弛豫算子里体现为额外相位旋转二是读出期间的T2衰减不是T2而是T2因为存在不均匀性每个采样点时横向磁化额外乘exp(−t/T2*)。这两个因子加入后仿真结果才具备预测真机行为的能力。6. 常见问题速查表症状可能原因解决办法重建图像整体旋转投影角度坐标旋转方向反了检查theta与cos/sin配对确认坐标系定义一致图像中心有亮斑/十字k空间原点过采后加权不均重建时做密度补偿用ramp filter或Voronoi权重放射状条纹过多spoke数不足增加到π/2倍矩阵尺寸以上图像模糊、边缘不锐T2*衰减过短或读出窗过长缩短TE或读出时长提高读出带宽整体亮度偏低翻转角远小于恩斯特角用α≈sqrt(2·TR/T1)重新估算翻转角周期带状伪影扰相不彻底检查扰相梯度矩加入RF相位循环对比度与解析公式不一致预热TR不足dummy TR加到50个以上再采数7. 后续扩展方向这个工程做完之后我建议往三个方向扩展性价比都不错。第一是三维化把体素网格从2D扩到3D梯度编码增加一个kz维度径向采集改成3D径向或锥形采样stack-of-stars重建用3D NUFFT这套组合在动态对比增强和超短回波UTE成像里非常主流。第二是加运动模型在逐TR循环里让体模做刚性平移或旋转直接观察径向采集对运动的鲁棒性这是radial序列相对笛卡尔的核心卖点仿真里能定量展示。第三是接深度学习用本工程批量生成不同参数下的仿真数据作为重建网络的训练配对或数据增强来源这类数据的物理一致性比纯GAN生成的数据强得多。我个人实操中的体会是布洛赫模拟这个方向最忌讳把代码写得“只能跑一次”。参数配置层、序列编排层、求解层分离得好后续要加bSSFP、EPI、Dixon这类新序列时只需要替换序列编排层的几个函数求解层和重建层完全复用。仿真代码的复用性和可读性往往比单次运行速度更值得花时间去打磨。最后分享一个小习惯每个版本的模拟配置和结果图我都会按日期归档参数变化用git记录这样复查问题时能立刻回到产生异常结果的原始环境省下的时间远比维护成本多。
做磁共振成像MRI序列开发的同事应该都有体会一个新想法直接搬到扫描仪上做实验调试周期长、机时成本高信号链里混入的各种系统误差还会干扰你对序列本身的判断。所以业界通行的做法是先在Matlab里做布洛赫模拟Bloch simulation把脉冲序列在“理想世界”里的行为跑透再决定要不要上机。这个项目做的就是这件事用二维布洛赫模拟实现FLASH序列的投影k空间radial k-space采集把序列参数、体模定义、k空间填充轨迹、重建图像全部放进一个Matlab工程里跑通还附带可直接复用的源码工程。先说清楚这个工程适合谁。如果你是刚接触MRI物理的硕士生想搞明白FLASH序列为什么用低翻转角、k空间为什么可以沿径向填充这份代码能让你一步步看到从自旋磁化矢量到k空间信号再到图像的完整链路。如果你是在做序列开发或图像重建的研究人员这个模拟工程可以用作新序列的验证基线——在仿真里确认正确性再动手改真机序列。核心价值一句话把“FLASH 径向k空间采样 图像重建”这条链路以参数完全可控、结果可反复检查的方式完整复现在电脑上。1. 项目背景与序列设计思路1.1 FLASH序列为什么能“快”FLASH全称Fast Low Angle Shot是梯度回波GRE家族的经典序列。它的设计目标用一句话概括在尽量短的重复时间TR内拿到可用的图像对比度。为此它做了两个关键取舍一是用梯度反转代替180°重聚脉冲来形成回波省掉了自旋回波序列里耗时的重聚脉冲二是把激发翻转角压得很低通常5°~30°让纵向磁化在每次激发后不会彻底耗尽稳态下仍有充裕的纵向磁化供下次激发使用。这两个特点叠加FLASH才能在几十毫秒甚至几毫秒的TR下连续成像成为临床3D扰相梯度回波如腹部T1增强扫描、颅脑高分辨T1的基本盘。你可能会问翻转角低不会让信号变弱吗这里的关键是稳态。短TR下若用90°翻转角一次激发后纵向磁化几乎清零TR余量不够它恢复下一拍的信号反而很弱倒是小角度激发让每次“消耗”的纵向磁化不多稳态时能持续产生稳定的横向磁化。这就是恩斯特角Ernst angle的理论来源后面我会专门算一笔账。1.2 投影k空间采集与笛卡尔采样的取舍常规MRI的k空间是逐行采的每次TR采一条平行线效率高、重建简单一维FFT即可。但笛卡尔采集有个明显短板对运动特别敏感扫描中一旦发生移动对应的相位编码线直接错乱伪影沿着相位编码方向蔓延整幅图几乎报废。径向投影采集则换了个思路每次TR沿一条穿过k空间中心、角度旋转的直线采样很像CT里的Radon变换测量。好处有三个一是k空间中心被反复过采信号能量最集中的区域天然得到强化对运动造成的相位误差有更强的容忍度二是运动伪影从“整图错位”变成相对柔和的放射状条纹信息损失更局部三是配合黄金角角度步进可以随时截断数据做动态重建时间分辨率灵活。代价也同样明显重建不能直接FFT需要网格化gridding或滤波反投影spoke数不足时图像上会出现明显的放射状条纹伪影。这个项目把radial读取的每个细节——角度步进、读出梯度持续时间、k空间填充坐标——全部纳入布洛赫模拟目的就是把这些权衡看得清清楚楚。2. 布洛赫方程与FLASH信号的定量关系要模拟MRI信号绕不开布洛赫方程。MRI信号本质上是宏观磁化矢量M(Mx, My, Mz)绕主磁场B0进动时在接收线圈中感应出的电动势布洛赫方程就是描述这个矢量时间演化的动力学方程dMx/dt γ(M × B)x − Mx/T2 dMy/dt γ(M × B)y − My/T2 dMz/dt γ(M × B)z − (Mz − M0)/T1其中γ是旋磁比B是总磁场主磁场射频场梯度场T1和T2分别是纵向和横向弛豫时间。这个方程组看着不复杂但射频脉冲、梯度、弛豫同时作用时很少有人手解数值仿真成了标配。2.1 把布洛赫方程拆开看布洛赫方程可以拆成三项理解。进动项γ(M×B)描述自旋绕磁场转动这是信号产生的来源——横向磁化绕主磁场进动在线圈中感应出可测量的复数信号梯度场能实现空间编码正是因为它改变不同位置的进动频率让磁化在读出时间内累积出与位置成正比的相位。弛豫项−Mx/T2、−My/T2描述横向磁化散失T2越短信号衰得越快纵向项−(Mz−M0)/T1描述纵向磁化向平衡值恢复T1越短恢复越快。在FLASH的一个TR周期里这三类过程按时间顺序重复上演射频脉冲把纵向磁化“扳倒”成横向磁化TE等待期间横向磁化一边进动一边按T2*衰减TR余量内纵向磁化按T1恢复。仿真要做的事就是把这三种算子按序列时间轴依次作用到每个体素的磁化矢量上。2.2 稳态信号公式与恩斯特角对扰相FLASHTR末把横向残留磁化spoiler掉来说稳态下的横向磁化有很漂亮的解析表达式Mxy M0 · sin(α) · (1 − E1) / (1 − cos(α)·E1) · E2其中α是翻转角E1 exp(−TR/T1)E2 exp(−TE/T2*)。当TR远小于T1时使信号最大化所需翻转角近似为α(Ernst) ≈ sqrt(2·TR/T1)假设T1800 ms、TR10 ms估算出来翻转角约5.8°考虑到T2*加权和组织差异实际常取8°~15°。模拟时我习惯先用这个公式算一个理论最优角再用布洛赫模拟扫一组翻转角验证二者趋势一致说明模拟实现没有原则性错误。这里多说一句解析公式假设了理想扰相、无偏共振、读出瞬间完成真实序列行为与它有偏差是正常的模拟的价值恰恰在于把这些“非理想”因素逐个加回去。2.3 为什么有了解析公式还要数值模拟既然稳态有解析解为什么还要费劲做布洛赫模拟因为解析解做了太多理想化假设忽略偏共振失相位化学位移、主磁场不均、忽略读出采样点之间的横向磁化衰减、忽略序列启动阶段的瞬态过程。真要设计新序列——比如改扰相梯度、加脂肪抑制准备脉冲、调读出方向——这些假设全都站不住必须回到数值解。这个项目的做法是逐体素更新磁化矢量在各体素上按时间顺序施加射频激发、梯度相位累积、T1/T2弛豫算子完整还原磁化矢量的“旅程”。虽然比套公式慢但能看到每一个中间状态这才是理解序列行为最直观的路径。3. Matlab仿真框架设计与实现这一节是实操核心。我把仿真架构分成五层参数配置层、体模定义层、序列事件编排层、布洛赫求解层、k空间填充与重建层。参数配置层集中存放TR、TE、翻转角、读出点数、spoke数、FOV、网格分辨率体模定义层生成一个有真实T1/T2/质子密度分布的二维目标序列事件编排层把“激发—等待TE—读出采样—恢复等待”循环按时间轴排成事件清单布洛赫求解层在每个体素上按事件清单推进磁化矢量最后把每个读出采样点的信号按坐标放入k空间。3.1 总体架构与数据流数据流一句话概括体模每个网格点一个磁化矢量→ 序列事件驱动磁化矢量演化 → 读出窗口内每隔一个采样间隔Δt采一个复数信号样本 → 按投影角度填入k空间 → 所有spoke采完后重建成图像。为了让代码既清晰又算得动两个细节我建议这样处理。第一每个体素看成一个独立的自旋等色群isochromat体素内的T1/T2和偏共振频率取该格子的中心值。加密体素数量可以减少磁化矢量“离散化”带来的误差但计算量按网格尺寸平方增长64×64和128×128的差别非常明显起步阶段用小网格调通逻辑再放大。第二梯度编码不通过逐时刻积分整个梯度波形实现而是直接按“k空间位置 读出采样点 × 对应方向梯度矩”计算相位。这与MR物理中k空间的定义完全一致也大幅节约计算量——毕竟真正的布洛赫积分要细到微秒级时间步全网格跑一遍代价很高。3.2 体模选择与参数设置仿真实验体模我建议先用Shepp-Logan数字体模它包含多个不同T1/T2/质子密度的区域足以观察对比度和伪影行为。Matlab自带phantom函数生成几何形状自己做一层T1/T2映射表即可。要注意体模网格分辨率必须与k空间读出点数匹配图像矩阵若为128×128体模也定义在128×128上读出点数至少取256沿直径方向过采样一倍spoke数按奈奎斯特准则至少需要约π/2倍矩阵点数也就是大约201条。实践里我常取的组合是128×128体模、256读出点、256条spoke既有明显的径向采集特征又不至于让计算量失控。我常用的基准参数如下表参数数值说明主磁场场强1.5 T旋磁比按γ42.58 MHz/TTR / TE10 ms / 3 ms短TR快速成像典型值翻转角10°接近恩斯特角附近采样带宽±31.25 kHz决定采样间隔Δt读出采样点数256每条spoke样本数投影spoke数256按π/2倍矩阵尺寸过采样图像矩阵128×128重建网格预热TR数50让磁化矢量进入稳态3.3 序列事件怎么编排序列的时间轴是布洛赫模拟的“剧本”。一个TR循环包含三段射频激发段、回波形成与读出段、驰豫恢复段。每段对应不同的更新算子。激发段射频脉冲把磁化矢量绕某个轴旋转旋转角就是翻转角。矩形脉冲的射频持续约0.5~1 ms仿真里可以简化为瞬时旋转相比TR内的弛豫时间很短误差可接受。旋转矩阵的具体形式是绕xy平面内相位为φ的轴转α角用标准三维旋转矩阵即可实现。读出段这一段的操作最密集。横向磁化绕z轴按“偏共振频率梯度场贡献”的总频率旋转同时按T2*衰减纵向磁化按T1恢复。读出窗口内每个采样间隔Δt采一次信号信号值就是当前横向磁化矢量的复数表示Mx i·My。梯度场贡献在投影方向上的体现是体素位置沿投影方向的坐标乘以梯度矩得到累积相位这正是k空间编码的本质。恢复段没有射频和梯度只有T1/T2弛豫。对扰相FLASH这一段的末尾还要模拟扰相操作把残余横向磁化清零确保下一个TR开始时横向磁化从零开始。3.4 核心Matlab代码实现下面给出几个最关键的函数片段。第一个是磁化矢量的RF旋转function M applyRF(M, alpha, phi) % alpha: 翻转角单位rad % phi: RF脉冲的初始相位单位rad % 标准三维旋转绕xy平面内方位角为phi的轴旋转alpha c cos(phi); s sin(phi); ca cos(alpha); sa sin(alpha); R [ca c^2*(1-ca), s*(1-ca), s*sa; s*(1-ca), ca s^2*(1-ca), -c*sa; -s*sa, c*sa, ca]; M R * M; end第二个是包含弛豫与偏共振的推进函数自由进动和弛豫可以合并成闭式解function M relax(M, dt, T1, T2star, dfreq) % dt: 持续时间 % T1,T2star: 纵向/有效横向弛豫常数 % dfreq: 偏共振频率(Hz) E1 exp(-dt/T1); E2 exp(-dt/T2star); phi 2*pi*dfreq*dt; % 偏共振导致的相位累积 Mc E2 * (M(1)1i*M(2)) * exp(1i*phi); M(1) real(Mc); M(2) imag(Mc); M(3) M0 - E1*(M0 - M(3)); % M0为平衡磁化 end第三个是主循环展示每个TR里体素磁化矢量的演进和信号采集for iSpoke 1:Nspokes theta (iSpoke-1) * pi / Nspokes; % 投影角度均匀覆盖180° Gx G0 * cos(theta); Gy G0 * sin(theta); % 5.1 射频激发 M applyRF(M, alpha, 0); % 5.2 TE等待自由进动 弛豫 M relax(M, TE, T1, T2star, dfreq); % 5.3 读出采样每个采样点对应一个k空间坐标 for iS 1:Nread kpos (iS - Nread/2 - 0.5) * dk; % 梯度编码相位位置在投影方向的贡献 phase_grad 2*pi * kpos * (Gx*X Gy*Y); sig (M(1) 1i*M(2)) .* exp(-1i*phase_grad); signal(iSpoke, iS) sum(sig(:)); % 所有体素相干叠加 end % 5.4 恢复段TR剩余时间 理想扰相 M relax(M, TR - TE - readTime, T1, T2, 0); M(1) 0; M(2) 0; % 理想扰相清空横向磁化 end注意上面这个写法是教学化的双层循环变量名里X、Y是meshgrid之后的位置坐标数组relax函数要支持对整数组体素同时运算。实际工程一定要矢量化把体素展开成一维数组所有体素同时推进信号通过矩阵乘法累加。我实测128×128体模、256条spoke、50个预热TRMatlab向量化后单次完整模拟大约1~2分钟完全可以接受。3.5 预热TR与稳态收敛这是新手最容易忽略的一步直接从M0开始模拟并记录信号得到的是瞬态信号不是FLASH真正工作的稳态信号。磁化矢量要经过多个TR才能进入动态平衡。处理办法是在正式采数前先跑20~50个“预热TR”只推进磁化矢量不记录信号。预热TR数量取决于T1/TR比值。以T1800 ms、TR10 ms为例20个TR后纵向磁化恢复到平衡值的1−exp(−20×10/800)≈22%还差得远想要98%的恢复大约需要800/10×ln(50)≈31个TR。所以我统一取50个预热TR代价只是增加约1/4的总计算量换来的稳态准确性对仿真可信度很重要。4. 仿真结果怎么看从k空间到图像仿真做完不能只在Matlab工作区里存几个数组要把过程“看见”。我通常做三件事画k空间填充轨迹、看每条spoke的原始信号幅度、重建成图像并做参数对比。4.1 k空间填充轨迹验证早期最常见的bug之一就是k空间坐标算错导致重建图像旋转或折叠。验证方法很简单把每条spoke的(kx, ky)坐标画成散点图正常应是以原点为中心、辐射状均匀覆盖半径范围的多段直线。角度步长为均匀π/Nspokes时相邻spoke的角度间隔应完全相等。如果看到坐标点乱窜、或某一方向的“扇区”明显稀疏问题多半出在梯度矩计算或采样索引偏移上。这一关必须在重建前过掉否则后面所有图像结论都不可靠。4.2 径向k空间重建从数据到图像径向k空间不能直接做笛卡尔FFT。常用两条路一是网格化重建gridding把非笛卡尔坐标上的样本用卷积核凯撒窗等插值到笛卡尔网格再做FFT二是滤波反投影FBP把每条spoke数据先做一维傅里叶变换再沿对应角度反投影到图像域本质是逆Radon变换。Matlab里做FBP可以直接用radon和iradon函数配合处理做gridding则可以用较新版本内置的nufft% ktraj为每条spoke每个采样点的归一化k空间坐标Nx×Ny为重建网格 ktraj buildKtraj(Nspokes, Nread); % 按投影角度和读出位置生成 data signal(:); % 采集信号展平 img nufft2(data, ktraj, [Nx Ny]); % 2D NUFFT重建用NUFFT重建时要注意坐标归一化范围不同版本的nufft函数对输入坐标域的定义不完全一致建议先用一个已知解析信号比如点扩散函数仿真验证重建正确性再正式使用。4.3 参数对比翻转角、spoke数、T2*的直观效果这部分是该模拟最有教学价值的地方。固定其他参数把翻转角从5°、10°、20°、35°扫一遍重建图的平均信号强度变化应当符合恩斯特角曲线把spoke数从256降到64再降到32图像会明显出现放射状条纹伪影把T2从30 ms改到10 ms读出时间内信号衰减加剧图像会变模糊、边缘变钝这就是典型的T2模糊效应。这三组对比都以子图形式保存既方便写实验报告也方便后续调参时做对比参照。5. 调试中踩过的坑与实战经验仿真不难难在让仿真结果可信。以下是我在这个项目里实际踩过、也常见同行掉进去的坑。5.1 单位制混乱导致所有参数失效梯度矩、k空间坐标、频率这三个量之间单位极易搞混。一个稳妥的口令是全部用国际单位制m、s、T、Hz旋磁比按γ42.58×10^6 Hz/T注意是Hz/T若要用rad单位需要乘2π。读出采样间隔Δt1/带宽读出梯度幅值Gread带宽/(γ·FOV)k空间步长Δk1/FOV。所有推导都建立在kγ∫Gdt这个基本关系上。一旦出现GHz·m级别的奇葩数量级先查单位不要着急查逻辑。我在一个早期版本里把旋磁比写成42.58×10^6 rad/s/T结果所有梯度矩偏小2π倍图像整体错位查了整整半天才定位至今记忆犹新。5.2 稳态没收敛就采数图像出现“渐变”伪影预热TR不足时前面若干条spoke的横向磁化比后面小k空间中央的信号权重异常重建图像呈现亮度不均或奇怪的渐变伪影。排查方法是画出每条spoke的中心信号幅度k空间原点样本随spoke序号的变化正常应是一条先快速下降后趋平的曲线如果看到持续波动把预热数加倍再看。5.3 扰相不彻底导致T1/T2混合对比FLASH是扰相GRE序列前提是每个TR末把横向磁化残留清掉。如果扰相梯度矩不足、或代码里压根没实现扰相步骤稳态信号里会混入相干横向磁化的贡献对比度就会偏向bSSFP如FIESTA/TrueFISP那种T1/T2混合。判断方法很直接对多个偏共振频率做模拟看信号是否随频率振荡——扰相不彻底时信号会出现明显振荡。仿真里直接把Mx、My置零是“理想扰相”真机需要用梯度扰相加射频相位循环RF spoiling来逼近这个落差本身也是很好的研究对象。5.4 读出点数与FOV不匹配导致混叠折返radial采集的每条spoke本质上是沿某方向的1D频率编码。若读出点数乘以Δk小于需要覆盖的k空间半径约1/(2×体素尺寸)高频信息就缺失图像边缘会出现混叠折返。经验公式Nread至少等于图像矩阵尺寸才能完整覆盖直径方向想要过采样余量就用1.5~2倍。我在工程注释里特意写了这条经验公式因为这是改参数时最容易忽略的连锁反应。5.5 偏共振与T2*衰减缺失结果过于“理想”不少初版模拟只做RF激发和梯度相位忽略偏共振和T2*。这样仿出来的图像特别干净漂亮但一上真机就对不上。我的建议是从一开始就在模拟里加入两个“现实因子”一是每个体素的偏共振频率dfreq可由主磁场不均匀、化学位移、脂肪成分等引起在弛豫算子里体现为额外相位旋转二是读出期间的T2衰减不是T2而是T2因为存在不均匀性每个采样点时横向磁化额外乘exp(−t/T2*)。这两个因子加入后仿真结果才具备预测真机行为的能力。6. 常见问题速查表症状可能原因解决办法重建图像整体旋转投影角度坐标旋转方向反了检查theta与cos/sin配对确认坐标系定义一致图像中心有亮斑/十字k空间原点过采后加权不均重建时做密度补偿用ramp filter或Voronoi权重放射状条纹过多spoke数不足增加到π/2倍矩阵尺寸以上图像模糊、边缘不锐T2*衰减过短或读出窗过长缩短TE或读出时长提高读出带宽整体亮度偏低翻转角远小于恩斯特角用α≈sqrt(2·TR/T1)重新估算翻转角周期带状伪影扰相不彻底检查扰相梯度矩加入RF相位循环对比度与解析公式不一致预热TR不足dummy TR加到50个以上再采数7. 后续扩展方向这个工程做完之后我建议往三个方向扩展性价比都不错。第一是三维化把体素网格从2D扩到3D梯度编码增加一个kz维度径向采集改成3D径向或锥形采样stack-of-stars重建用3D NUFFT这套组合在动态对比增强和超短回波UTE成像里非常主流。第二是加运动模型在逐TR循环里让体模做刚性平移或旋转直接观察径向采集对运动的鲁棒性这是radial序列相对笛卡尔的核心卖点仿真里能定量展示。第三是接深度学习用本工程批量生成不同参数下的仿真数据作为重建网络的训练配对或数据增强来源这类数据的物理一致性比纯GAN生成的数据强得多。我个人实操中的体会是布洛赫模拟这个方向最忌讳把代码写得“只能跑一次”。参数配置层、序列编排层、求解层分离得好后续要加bSSFP、EPI、Dixon这类新序列时只需要替换序列编排层的几个函数求解层和重建层完全复用。仿真代码的复用性和可读性往往比单次运行速度更值得花时间去打磨。最后分享一个小习惯每个版本的模拟配置和结果图我都会按日期归档参数变化用git记录这样复查问题时能立刻回到产生异常结果的原始环境省下的时间远比维护成本多。
