尧图网站开发 
1. 这个项目到底在模拟什么FLASH、k空间与布洛赫方程如何串成一条线如果你搜FLASH这个词大概率会搜出一堆闪存颗粒型号、网页播放器历史、甚至动画软件的老黄历。但做MRI序列仿真的人听到FLASH脑子里只有四个字快速小角度。Fast Low Angle Shot一种被用到极致的基础梯度回波序列。这个项目的核心就是在Matlab里用布洛赫方程把这套序列的物理过程完整跑一遍并且把采集方式从最常见的笛卡尔k空间换成投影式radialk空间采集。我最初想做这件事的动机很简单在真正上机调序列之前先搞清楚翻转角、TR、TE这些参数会对最终图像产生什么影响。MRI序列的调试成本高一次扫描时间长参数改错了很可能要重扫。而布洛赫仿真就是在电脑上先把信号机制推演一遍让你在躺着等扫描的时候心里已经对结果有了预判。这个项目最适合三类人看一是研究生阶段要搭建脉冲序列仿真的同学二是想搞懂FLASH对比度机制而不是只会用序列参数的工程师三是准备用仿真数据测试重建算法的图像处理背景的朋友。先说清楚整件事的逻辑链。MRI的原始数据是k空间数据也就是磁化矢量的空间频域。FLASH序列负责决定在什么时刻、用什么翻转角、沿着什么轨迹去读取这个k空间。布洛赫方程负责决定每个体素里的磁化矢量在射频激励、梯度编码、T1/T2弛豫的影响下到底变成了什么状态。三者缺一不可没有布洛赫方程你只能假设信号幅度是常数没有k空间轨迹你采集到的只是一串没有空间编码意义的FID没有FLASH的参数逻辑布洛赫方程也不知道自己该在什么时候被激发。所以这项目本质上不是在写一个成像算法而是在搭一条完整的物理链路。这条链路里最容易被忽略的是投影采集这个关键词。常规笛卡尔采集是一行一行地把k空间填满投影采集则是每一条k空间线都穿过原点呈放射状分布。两条路线的物理意义完全不同后面模拟时对梯度编码和重建的处理方式也完全不同。1.1 布洛赫方程在这条链路里充当什么角色布洛赫方程描述的是一根磁化矢量在外加磁场下的运动规律标准形式是dM/dt γM×B - (Mx/T2, My/T2, (Mz-M0)/T1)这个方程看着复杂但物理图像其实很直观磁化矢量会绕着磁场方向进动同时纵向分量往M0恢复横向分量不断衰减。射频脉冲相当于把这个矢量掰一个角度梯度场则是让不同位置的矢量获得不同的进动频率从而在k空间形成编码。在仿真里我们不会真的去用数值积分器解这个微分方程。对于MRI序列这种脉冲-等待-采集的结构化过程完全可以把时间演化拆成一个个离散小步骤先算射频激励造成的旋转再算自由进动造成的相位积累最后算弛豫造成的衰减。这种拆法在物理上有根据在工程上也好实现它让布洛赫方程从偏微分方程变成了一串矩阵运算。1.2 二维这个词到底指什么标题里的二维其实承担了双重含义。第一重是空间维度的二维——模拟的对象是一个x-y平面上的二维对象每个像素位置都有一根独立的磁化矢量。第二重是布洛赫方程的横向平面——Mx和My组成的复数平面。做仿真时两个含义都要顾到空间二维决定网格和k空间坐标横向平面决定信号怎么积累。刚开始搭模型时容易在这两个维度上打架。比如你想模拟一个128×128的对象又要对每个像素维护Mx、My、Mz三个分量最直接的想法是建一个128×128×3的三维数组。但实际写下去你会发现把所有像素的磁化矢量都存成矩阵然后整块做旋转和弛豫运算比用循环逐体素处理要快出好几个数量级代码也更接近布洛赫方程的原始数学形式。1.3 整条链路的任务清单我最后把项目拆成了六个模块顺序绝对不能乱建立二维空间网格定义组织参数M0、T1、T2、B0偏移分布。初始化每组isochromat的磁化矢量。生成投影式k空间轨迹确定每条线的角度和采样点。按FLASH序列逐TR循环每个周期依次完成激发、等待、读出采样、扰相。将采集到的数据按角度和采样序号存入k空间矩阵。通过滤波反投影或网格化重建得到图像与原始phantom对比。这套清单其实就是本文的骨架。下面每一章会把其中一到两个模块彻底拆开包括原理推导、代码骨架和我在调参过程中踩过的坑。2. 布洛赫方程的离散化实现旋转算子与弛豫项的Matlab代码骨架布洛赫方程离散化这件事网上资料很多但大多数版本要么只算一个孤立体素要么把弛豫项写得特别抽象拿到实际图像仿真里根本用不了。我在这里给出一个能直接铺到网格上的版本并且解释每一步为什么这么写。2.1 为什么可以把时间演化拆成旋转弛豫考虑一个很短的时间间隔dt。在这个dt里磁化矢量经历了两种物理过程一是绕磁场方向进动也就是旋转二是纵向和横向的弛豫也就是衰减和恢复。在dt足够小的前提下这两个过程可以近似为依次发生先转一个角度再做弛豫衰减。射频激励脉冲也有两种建模方式。一种是真正的时变B1场把脉冲波形离散成几十个时间点每个时间点做一次小旋转。另一种是理想化处理假设脉冲持续时间远小于TR直接看成一次瞬时旋转。对于FLASH这种翻转角通常在5°到20°之间的小角度序列瞬时旋转的近似误差完全可以接受。我在仿真里用的就是瞬时旋转效果跟分步模拟差别极小但代码量少了一大截。旋转算子的具体形式分两种。射频激励绕x轴旋转角即翻转角αMx Mx My My·cosα - Mz·sinα Mz My·sinα Mz·cosα自由进动绕z轴旋转角由B0偏移和梯度场共同决定。这一步在采集编码时尤其重要因为梯度编码的本质就是让不同位置的矢量积累不同的进动相位。弛豫项则单独处理Mxy Mxy·exp(-dt/T2) Mz M0 (Mz - M0)·exp(-dt/T1)这三条公式是整套仿真的地基。所有后续代码包括梯度编码、B0不均匀性、投影采集都是在这三步算子之上叠加。2.2 isochromat怎么组织用矩阵而不是逐体素循环一开始我的代码长得很蠢三层for循环嵌套先遍历y坐标再遍历x坐标每个像素重建一个布洛赫状态去演化算一次64×64的图像要跑好几分钟。后来我意识到Matlab的强项是矩阵运算布洛赫状态的三个分量完全可以各自拆成一个矩阵所有物理过程都整块作用上去。具体做法是Mx、My、Mz各是一个NX×NY的矩阵。T1、T2、M0也是同尺寸矩阵这样每个像素可以有自己独立的地缘物理参数。旋转操作直接对三个矩阵做线性组合。弛豫操作使用矩阵点乘和点除。这套写法的另一个好处是它让你很自然地处理组织模型。比如我要建一个中心是高信号的phantom中心区域的T1600ms、T270ms外围T11000ms、T2120ms只需要先算一个半径距离矩阵然后用逻辑索引把对应的参数矩阵区域改掉就行。不需要为不同组织写不同的演化分支。2.3 最简可运行的演化代码骨架下面是我在这个项目里实际使用的代码骨架。为了方便阅读我把旋转和弛豫各自封装成了函数。% 网格参数 NX 64; NY 64; FOV 256; % mm x linspace(-FOV/2, FOV/2 - FOV/NX, NX); [X, Y] meshgrid(x, x); % 组织参数中心圆区域T1/T2偏短模拟高信号组织 R sqrt(X.^2 Y.^2); T1 1000 * ones(NX, NY); T2 120 * ones(NX, NY); M0 ones(NX, NY); maskCenter R 40; T1(maskCenter) 600; T2(maskCenter) 70; M0(maskCenter) 1.5; % B0偏移场单位Hz单等色模型里先置零 dB0 zeros(NX, NY); % 磁化矢量初始化 Mx zeros(NX, NY); My zeros(NX, NY); Mz M0; % FLASH序列参数 alpha 12 * pi / 180; TR 15; % ms TE 6; % ms % 射频激励绕x轴旋转 function [Mx, My, Mz] applyExcitation(Mx, My, Mz, alpha) MyNew My .* cos(alpha) - Mz .* sin(alpha); MzNew My .* sin(alpha) Mz .* cos(alpha); Mx Mx; My MyNew; Mz MzNew; end % 弛豫与自由进动演化 function [Mx, My, Mz] applyEvolution(Mx, My, Mz, dt, T1, T2, dB0) % dt单位msdB0单位Hz相位中的时间要转换为秒 phase 2 * pi * dB0 * dt / 1000; Mxy complex(Mx, My); Mxy Mxy .* exp(-dt ./ T2) .* exp(-1i * phase); Mx real(Mxy); My imag(Mxy); Mz M0 (Mz - M0) .* exp(-dt ./ T1); end注意applyEvolution函数里M0是通过外部变量传入的我这里为了简洁省略了传参实际编写时记得把M0也传进去。这里有个细节dB0是Hz单位dt是ms单位相位公式里的时间必须除以1000换成秒否则相位会差1000倍。这个单位坑我踩过后面专门有一章讲。3. 投影式k空间轨迹的坐标生成以及它和笛卡尔采样的本质差异项目标题里最显眼的关键词除了FLASH就是投影。我做这个项目之前对投影采集的理解停留在每条线过原点这个层面真正写代码才发现轨迹数学、采样点数、角度数目、编码相位定义每一步都有讲究。3.1 笛卡尔与径向的k空间覆盖形态笛卡尔采集一次TR只填充一行k空间相邻行间隔均匀ky方向通过相位编码步进。它的优点是数据天然落在矩形网格上重建直接做二维FFT就行。缺点是k空间中心区域和边缘区域的采样密度相同中心欠采样带来的问题不大但长的采样时间让整体扫描变慢。径向采集则完全不同每条k空间线从-kmax到kmax扫过原点角度均匀分布在0到π之间。这些线在k空间中心区域密集交叉越往外越稀疏。这个采样形态带来一个天生优势k空间中心的低空间频率被反复测量等效于极强的信号平均对运动伪影的容忍度也高。投影FLASH在极短TR的实时成像里很常见因为它对呼吸和心跳这类整体运动的敏感性远低于笛卡尔序列。但径向采样的代价也很明显边缘区域采样稀疏不满足各向同性的奈奎斯特判据时会产生放射状的星芒伪影。重建时也不能直接FFT必须先做网格化重采样或者走滤波反投影路线。3.2 径向轨迹的数学生成与参数设定投影轨迹的数学可以一句话说清每条线有一个固定角度θnk空间坐标从-kmax均匀扫到kmax线上每一点同时对应x轴和y轴方向的梯度分量。k空间采样点数Nsamp和视野FOV决定了kmax。以FOV256mm、Nsamp128为例kmax Nsamp / (2 × FOV) 128 / 512 0.25 cycles/mm对应的空间分辨率是1/(2×kmax) 2mm正好等于FOV/Nsamp。这个关系必须成立否则重建图像的空间尺度和像素间距会对不上。角度生成方面初期我用均匀分布角度Nproj 180; angles (0:Nproj-1) * pi / Nproj;后来我也试过黄金角增量每次递增约111.25°。黄金角的优势在于任意连续数据块里的投影角度都近似均匀适合动态成像和实时重建。但稳态仿真第一步用均匀角度就够探索径向轨迹的统计学特性时再上黄金角。采样点k坐标生成Nsamp 128; dk 1 / FOV; kvals (-Nsamp/2 : Nsamp/2 - 1) * dk;这里kval的起点是负值、终点是正值k0正好落在数组中间。这个对称性在后面做ifftshift时非常关键。3.3 编码相位的符号约定——最容易错的一步有了k坐标之后每个体素在采集时积累的空间编码相位为phase 2π × k(s) × (x·cosθ y·sinθ)如果用矩阵化的写法在Matlab里是这样的% 第p条投影角度为theta for s 1:Nsamp phase 2 * pi * kvals(s) * (X * cos(theta) Y * sin(theta)); signal(s) sum((Mx 1i*My) .* exp(-1i * phase), all); end这里的指数符号是负号也就是约定信号s(t) ∫Mxy(r,t)·exp(-i2πk(t)·r)dr。这个约定不能随便更改因为后面的B0偏移相位、重建时的共轭操作全部要跟它保持同一符号体系否则重建图像会左右翻转。还有一个我直到现在都会条件反射检查的坑X和Y矩阵的坐标原点。我用的linspace(-FOV/2, FOV/2-FOV/NX, NX)会把坐标原点放在网格边缘的第一个像素和FOV中心之间偏离半个像素的位置。这个半像素偏移在绝大多数情况下不影响图像视觉效果但在做像素级对比时会造成整体平移。如果追求像素级对齐应该用x (-NX/2 : NX/2-1) * (FOV/NX);这个区别在对称phantom的自检里很容易暴露。后面第6章会教你怎么通过对称性检验发现这个问题。4. FLASH序列参数设定翻转角、TR/TE与对比度的权衡逻辑序列参数是布洛赫仿真的灵魂。同样的物理仿真代码参数不同结论完全不同。这一章不讲教科书的全部推导只讲我实际仿真FLASH时怎么定参数以及为什么这样定。4.1 低翻转角为什么是FLASH的灵魂FLASH从一开始就放弃了90°激励。核心原因是TR很短时90°激励后纵向磁化来不及恢复下一次激励能用的信号就很少。低翻转角把一部分纵向磁化保留下来换来的是每个TR都能有相对稳定的信号。对于稳态梯度回波带扰相稳态纵向磁化和信号幅度可以写成Mz_ss M0·(1-E1) / (1-cosα·E1) S ∝ sinα·Mz_ss·exp(-TE/T2*)其中E1 exp(-TR/T1)。让信号对翻转角求导等于零得到Ernst角αE arccos(exp(-TR/T1))这个公式的重要性在于它告诉你翻转角不是越大越好而是由TR和T1共同决定。举个例子T1700ms、TR15ms时E1 exp(-15/700) ≈ 0.979 αE arccos(0.979) ≈ 11.8°所以我在仿真里把翻转角设在12°既能获得最大稳态信号又保留足够的T1对比度。如果你把翻转角改成25°信号幅度不会更高反而会大幅降低纵向稳态磁化图像整体变暗。这是FLASH这个核最核心的机制。4.2 我的初始参数表我实际用的一组仿真参数如下参数数值说明网格尺寸64×64调试快重建也能看出问题FOV256mm对应像素分辨率4mm够用采样点数Nsamp128单条k线采样数投影数Nproj1800~179°间隔1°TR15ms稳态梯度回波的典型短TRTE6ms落在读取窗中心翻转角12°按Ernst角估算中心组织T1/T2600ms / 70ms模拟高信号实质外围组织T1/T21000ms / 120ms长T1对比dB0偏移0Hz单等色严格T2*见6.1节TE的选择逻辑是梯度回波没有180°重聚脉冲TE越短T2*衰减越少信号越强但TE也不能太短否则磁敏差异造成的相位演化还没来得及体现。6ms在大部分组织里是合理的折中。如果你想看明显的磁敏效应把TE拉长到15~20ms图像里就能看到组织交界处的信号流失。4.3 如何确认模拟进入了稳态FLASH是稳定序列但磁化矢量不是从第一个TR就立刻稳定的。初始时刻MzM0第一次12°激励后的信号会显著高于稳态。如果不处理这个瞬态采集到的第一条投影线幅度异常重建后会出现沿径向方向的条带伪影。处理办法有两个。我的做法是在正式采集前先跑一段热身循环warm-up让纵向磁化收敛到稳态warmupTRs 20; for idx 1:warmupTRs [Mx, My, Mz] applyExcitation(Mx, My, Mz, alpha); [Mx, My, Mz] applyEvolution(Mx, My, Mz, TE, T1, T2, dB0); % 理想扰相清除残余横向磁化 Mx(:,:) 0; My(:,:) 0; [Mx, My, Mz] applyEvolution(Mx, My, Mz, TR-TE, T1, T2, dB0); end每一次都清零Mx和My对应FLASH序列末端的扰相梯度或者射频扰相。如果你要严格模拟无扰相版本比如真正的bSSFP序列就不要清零但稳态条件和信号公式完全不同。稳态是否收敛可以看单个体素的Mz变化幅度。我写过一个脚本记录每次TR结束后的Mz平均值当相邻两次TR的相对变化小于1e-3时就认为进入了稳态。通常20个TR足够TR越短收敛越快。5. 重建路径与验证从径向k空间线到图像仿真采集到的数据是一组原始k空间线形状是Nproj×Nsamp的复数矩阵。接下来要做的是把它变成一张能肉眼判断好坏的图像。这一步如果不走对前面的布洛赫物理全部白搭。5.1 方法Aifft加滤波反投影——最稳径向k空间和图像之间的桥梁是中心切片定理每个角度的k空间线做一维逆傅里叶变换得到的正是该角度下物体的辐射投影。因此最直接的路径是先把每条k线变换为投影数据再用滤波反投影iradon重建。Matlab里实现很简洁% 每条k线先做ifftshift再ifft再ifftshift确保k0对齐到逆傅里叶中心 proj ifftshift(Kdata, 2); proj ifft(proj, [], 2); proj ifftshift(proj, 2); % iradon需要角度单位为度数且投影按探测器维度排列 angles_deg angles * 180 / pi; recon iradon(proj., angles_deg, linear, Ram-Lak);iradon默认内置了Ram-Lak斜坡滤波器正好补偿径向k空间中心过采样造成的低频权重过大。如果你自己写网格化重建这步滤波是无论如何也绕不开的。我第一次跑这个流程时重建出来的图像左右翻转且略有一像素偏移。后来发现两个原因一是ifftshift的方向和频域原点没对齐二是iradon内部对探测器原点的定义和我的网格定义差了半个像素。这类问题不用硬记直接用一个中心对称的phantom做测试看到翻转就翻转回来看到平移就调整ifftshift或对k线补零。5.2 方法B笛卡尔重采样加密度补偿另一种做法是把极坐标下的k线插值到笛卡尔网格再对网格做二维ifft。这种方法的优势是更接近MRI重建的主流实现gridding但需要密度补偿函数否则重建图会有明显的模糊和边缘振铃。密度补偿的基本逻辑是k空间中心区域每条线贡献的采样点过多权重必须降低。理想的权重是采样点周围区域的面积倒数即1/|k|越靠近边缘权重越大。我叫它越稀疏越重要。实现时可以用三角插值算出每个笛卡尔网格点周围受哪些极坐标样本影响然后按1/r归一化。但如果你既不想装JSON工具箱又不想手写复杂的网格化器方法A的iradon在这个仿真规模下完全够用。iradon自带滤波效果不需要显式做密度补偿。所以项目主线我推荐方法A方法B适合你想在重建算法上做文章的后续扩展。5.3 投影数量和星芒伪影的压力测试径向采集的伪影很有辨识度当投影数不够时重建图像会从高强组织向外辐射出尖刺状条纹。原因是径向角度间隔过大k空间外围各向异性采样造成的角域混叠。我做了一组对比仿真网格64×64、采样128点/线投影数分别是60、120、180、256。结果和经验公式基本吻合要避免可见星芒投影数需满足Nproj ≥ π × Nsamp / 2128采样点时π×128/2≈201。所以180条线会有轻微条纹256条才几乎看不见。这就解释了为什么不少径向MRI应用用黄金角或者高角度采样不只是为了动态时间分辨率更是为了压制角方向混叠。6. 调这套仿真时踩过的坑与自检方法最后这部分是纯干货清单。每一段都是我真真实实调试过的问题按排查顺序写希望能帮你少走弯路。6.1 单isochromat模型根本没有T2*信号衰减这是整套仿真里最大的一个物理陷阱。许多教程会告诉你在布洛赫模拟中引入一个B0偏移场dB0就能看到信号衰减模拟出T2*效果。这句话严格来说是错的。如果每个体素只有一根isochromat那么这支磁化矢量在B0偏移下只是整体旋转了一个相位模长完全不变信号幅度不会衰减。真正的T2*衰减来自一个体素内成千上万个自旋它们各自处在不同的微观场环境里进动频率不同随时间推移相位散开宏观信号相干性才逐渐消失。解决这个问题有三条路线。第一条严格路线在每个网格点内部再放几个子isochromat给它们各自的随机dB0偏移把所有子isochromat的信号叠加作为该体素输出。第二条工程近似路线直接给弛豫项换成exp(-dt/T2*)人为加入T2*衰减。第三条折中路线给每个体素分配一个均值附近的随机dB0分布比如高斯分布标准差5Hz体素间保留一定空间平滑性。我在FLASH仿真里用的是第二条加第一条的组合先按T2*简化跑通主流程再做一组子isochromat版本的对照实验验证趋势没变。如果你想严谨强烈建议用第一条因为它能自然模拟出梯度回波随TE衰减的曲线。6.2 单位、相位符号与采样时钟的一致性我在代码里处处强调单位因为它实在太容易出错了。dB0单位是Hz相位公式里要乘时间时间为毫秒时必须换算成秒再乘。t6ms、dB05Hz时相位 2π × 5 × 0.006 ≈ 0.188 rad如果忘了除以1000算出来的是188 rad完全没有物理意义信号会乱闪。相位符号也是如此。空间编码用exp(-i2πk·r)B0偏移相位就用同符号exp(-i2πdB0·t)。反向使用会导致k空间数据在频域里左右颠倒重建图像跟着左右颠倒。这一类错误很难肉眼发现只能靠对称性测试。还有一个隐蔽问题采样时钟和k空间坐标的关系。每条k线并非瞬时采完而是在读出窗内随时间匀速移动的。如果TE不落在采样窗正中心k0对应的就不应该是信号最大值位置而是偏在一边。FLASH的回波定义是TE时刻信号峰值所在的采样点所以要确保采样点序号和TE时刻对齐否则重建相位会错乱。6.3 用对称phantom做自检与参数扫描的收获我强烈建议你在布洛赫仿真的早期就准备一个关于原点对称的phantom。它可以是中心圆加四个对称小方块只要左右翻转后和原图一样就行。然后用仿真重建结果检查翻转后的图像是否等于原图像。如果不等于先查相位符号再查ifftshift方向再查网格原点定义。做完功能验证后我还做了一组翻转角扫描α分别取5°、10°、12°、15°、25°其他参数不变。重建两张图的相对信号强度结论落在12°附近表现最好和Ernst角公式预测的一致。这件事本身没什么科研价值但能直观说明仿真确实捕捉到了FLASH的物理本质——它不是在画画而是在复现一种真实的信号机制。我最后想说的是这套仿真里最有用的不是最终图像本身而是你在推导、编码、排错过程中建立起来的参数直觉。以后扔给你一台没有调试面板的MRI扫描仪你也能在脑子里预估改变翻转角会发生什么——这种能力才是做布洛赫仿真的最大回报。 
1. 这个项目到底在模拟什么FLASH、k空间与布洛赫方程如何串成一条线如果你搜FLASH这个词大概率会搜出一堆闪存颗粒型号、网页播放器历史、甚至动画软件的老黄历。但做MRI序列仿真的人听到FLASH脑子里只有四个字快速小角度。Fast Low Angle Shot一种被用到极致的基础梯度回波序列。这个项目的核心就是在Matlab里用布洛赫方程把这套序列的物理过程完整跑一遍并且把采集方式从最常见的笛卡尔k空间换成投影式radialk空间采集。我最初想做这件事的动机很简单在真正上机调序列之前先搞清楚翻转角、TR、TE这些参数会对最终图像产生什么影响。MRI序列的调试成本高一次扫描时间长参数改错了很可能要重扫。而布洛赫仿真就是在电脑上先把信号机制推演一遍让你在躺着等扫描的时候心里已经对结果有了预判。这个项目最适合三类人看一是研究生阶段要搭建脉冲序列仿真的同学二是想搞懂FLASH对比度机制而不是只会用序列参数的工程师三是准备用仿真数据测试重建算法的图像处理背景的朋友。先说清楚整件事的逻辑链。MRI的原始数据是k空间数据也就是磁化矢量的空间频域。FLASH序列负责决定在什么时刻、用什么翻转角、沿着什么轨迹去读取这个k空间。布洛赫方程负责决定每个体素里的磁化矢量在射频激励、梯度编码、T1/T2弛豫的影响下到底变成了什么状态。三者缺一不可没有布洛赫方程你只能假设信号幅度是常数没有k空间轨迹你采集到的只是一串没有空间编码意义的FID没有FLASH的参数逻辑布洛赫方程也不知道自己该在什么时候被激发。所以这项目本质上不是在写一个成像算法而是在搭一条完整的物理链路。这条链路里最容易被忽略的是投影采集这个关键词。常规笛卡尔采集是一行一行地把k空间填满投影采集则是每一条k空间线都穿过原点呈放射状分布。两条路线的物理意义完全不同后面模拟时对梯度编码和重建的处理方式也完全不同。1.1 布洛赫方程在这条链路里充当什么角色布洛赫方程描述的是一根磁化矢量在外加磁场下的运动规律标准形式是dM/dt γM×B - (Mx/T2, My/T2, (Mz-M0)/T1)这个方程看着复杂但物理图像其实很直观磁化矢量会绕着磁场方向进动同时纵向分量往M0恢复横向分量不断衰减。射频脉冲相当于把这个矢量掰一个角度梯度场则是让不同位置的矢量获得不同的进动频率从而在k空间形成编码。在仿真里我们不会真的去用数值积分器解这个微分方程。对于MRI序列这种脉冲-等待-采集的结构化过程完全可以把时间演化拆成一个个离散小步骤先算射频激励造成的旋转再算自由进动造成的相位积累最后算弛豫造成的衰减。这种拆法在物理上有根据在工程上也好实现它让布洛赫方程从偏微分方程变成了一串矩阵运算。1.2 二维这个词到底指什么标题里的二维其实承担了双重含义。第一重是空间维度的二维——模拟的对象是一个x-y平面上的二维对象每个像素位置都有一根独立的磁化矢量。第二重是布洛赫方程的横向平面——Mx和My组成的复数平面。做仿真时两个含义都要顾到空间二维决定网格和k空间坐标横向平面决定信号怎么积累。刚开始搭模型时容易在这两个维度上打架。比如你想模拟一个128×128的对象又要对每个像素维护Mx、My、Mz三个分量最直接的想法是建一个128×128×3的三维数组。但实际写下去你会发现把所有像素的磁化矢量都存成矩阵然后整块做旋转和弛豫运算比用循环逐体素处理要快出好几个数量级代码也更接近布洛赫方程的原始数学形式。1.3 整条链路的任务清单我最后把项目拆成了六个模块顺序绝对不能乱建立二维空间网格定义组织参数M0、T1、T2、B0偏移分布。初始化每组isochromat的磁化矢量。生成投影式k空间轨迹确定每条线的角度和采样点。按FLASH序列逐TR循环每个周期依次完成激发、等待、读出采样、扰相。将采集到的数据按角度和采样序号存入k空间矩阵。通过滤波反投影或网格化重建得到图像与原始phantom对比。这套清单其实就是本文的骨架。下面每一章会把其中一到两个模块彻底拆开包括原理推导、代码骨架和我在调参过程中踩过的坑。2. 布洛赫方程的离散化实现旋转算子与弛豫项的Matlab代码骨架布洛赫方程离散化这件事网上资料很多但大多数版本要么只算一个孤立体素要么把弛豫项写得特别抽象拿到实际图像仿真里根本用不了。我在这里给出一个能直接铺到网格上的版本并且解释每一步为什么这么写。2.1 为什么可以把时间演化拆成旋转弛豫考虑一个很短的时间间隔dt。在这个dt里磁化矢量经历了两种物理过程一是绕磁场方向进动也就是旋转二是纵向和横向的弛豫也就是衰减和恢复。在dt足够小的前提下这两个过程可以近似为依次发生先转一个角度再做弛豫衰减。射频激励脉冲也有两种建模方式。一种是真正的时变B1场把脉冲波形离散成几十个时间点每个时间点做一次小旋转。另一种是理想化处理假设脉冲持续时间远小于TR直接看成一次瞬时旋转。对于FLASH这种翻转角通常在5°到20°之间的小角度序列瞬时旋转的近似误差完全可以接受。我在仿真里用的就是瞬时旋转效果跟分步模拟差别极小但代码量少了一大截。旋转算子的具体形式分两种。射频激励绕x轴旋转角即翻转角αMx Mx My My·cosα - Mz·sinα Mz My·sinα Mz·cosα自由进动绕z轴旋转角由B0偏移和梯度场共同决定。这一步在采集编码时尤其重要因为梯度编码的本质就是让不同位置的矢量积累不同的进动相位。弛豫项则单独处理Mxy Mxy·exp(-dt/T2) Mz M0 (Mz - M0)·exp(-dt/T1)这三条公式是整套仿真的地基。所有后续代码包括梯度编码、B0不均匀性、投影采集都是在这三步算子之上叠加。2.2 isochromat怎么组织用矩阵而不是逐体素循环一开始我的代码长得很蠢三层for循环嵌套先遍历y坐标再遍历x坐标每个像素重建一个布洛赫状态去演化算一次64×64的图像要跑好几分钟。后来我意识到Matlab的强项是矩阵运算布洛赫状态的三个分量完全可以各自拆成一个矩阵所有物理过程都整块作用上去。具体做法是Mx、My、Mz各是一个NX×NY的矩阵。T1、T2、M0也是同尺寸矩阵这样每个像素可以有自己独立的地缘物理参数。旋转操作直接对三个矩阵做线性组合。弛豫操作使用矩阵点乘和点除。这套写法的另一个好处是它让你很自然地处理组织模型。比如我要建一个中心是高信号的phantom中心区域的T1600ms、T270ms外围T11000ms、T2120ms只需要先算一个半径距离矩阵然后用逻辑索引把对应的参数矩阵区域改掉就行。不需要为不同组织写不同的演化分支。2.3 最简可运行的演化代码骨架下面是我在这个项目里实际使用的代码骨架。为了方便阅读我把旋转和弛豫各自封装成了函数。% 网格参数 NX 64; NY 64; FOV 256; % mm x linspace(-FOV/2, FOV/2 - FOV/NX, NX); [X, Y] meshgrid(x, x); % 组织参数中心圆区域T1/T2偏短模拟高信号组织 R sqrt(X.^2 Y.^2); T1 1000 * ones(NX, NY); T2 120 * ones(NX, NY); M0 ones(NX, NY); maskCenter R 40; T1(maskCenter) 600; T2(maskCenter) 70; M0(maskCenter) 1.5; % B0偏移场单位Hz单等色模型里先置零 dB0 zeros(NX, NY); % 磁化矢量初始化 Mx zeros(NX, NY); My zeros(NX, NY); Mz M0; % FLASH序列参数 alpha 12 * pi / 180; TR 15; % ms TE 6; % ms % 射频激励绕x轴旋转 function [Mx, My, Mz] applyExcitation(Mx, My, Mz, alpha) MyNew My .* cos(alpha) - Mz .* sin(alpha); MzNew My .* sin(alpha) Mz .* cos(alpha); Mx Mx; My MyNew; Mz MzNew; end % 弛豫与自由进动演化 function [Mx, My, Mz] applyEvolution(Mx, My, Mz, dt, T1, T2, dB0) % dt单位msdB0单位Hz相位中的时间要转换为秒 phase 2 * pi * dB0 * dt / 1000; Mxy complex(Mx, My); Mxy Mxy .* exp(-dt ./ T2) .* exp(-1i * phase); Mx real(Mxy); My imag(Mxy); Mz M0 (Mz - M0) .* exp(-dt ./ T1); end注意applyEvolution函数里M0是通过外部变量传入的我这里为了简洁省略了传参实际编写时记得把M0也传进去。这里有个细节dB0是Hz单位dt是ms单位相位公式里的时间必须除以1000换成秒否则相位会差1000倍。这个单位坑我踩过后面专门有一章讲。3. 投影式k空间轨迹的坐标生成以及它和笛卡尔采样的本质差异项目标题里最显眼的关键词除了FLASH就是投影。我做这个项目之前对投影采集的理解停留在每条线过原点这个层面真正写代码才发现轨迹数学、采样点数、角度数目、编码相位定义每一步都有讲究。3.1 笛卡尔与径向的k空间覆盖形态笛卡尔采集一次TR只填充一行k空间相邻行间隔均匀ky方向通过相位编码步进。它的优点是数据天然落在矩形网格上重建直接做二维FFT就行。缺点是k空间中心区域和边缘区域的采样密度相同中心欠采样带来的问题不大但长的采样时间让整体扫描变慢。径向采集则完全不同每条k空间线从-kmax到kmax扫过原点角度均匀分布在0到π之间。这些线在k空间中心区域密集交叉越往外越稀疏。这个采样形态带来一个天生优势k空间中心的低空间频率被反复测量等效于极强的信号平均对运动伪影的容忍度也高。投影FLASH在极短TR的实时成像里很常见因为它对呼吸和心跳这类整体运动的敏感性远低于笛卡尔序列。但径向采样的代价也很明显边缘区域采样稀疏不满足各向同性的奈奎斯特判据时会产生放射状的星芒伪影。重建时也不能直接FFT必须先做网格化重采样或者走滤波反投影路线。3.2 径向轨迹的数学生成与参数设定投影轨迹的数学可以一句话说清每条线有一个固定角度θnk空间坐标从-kmax均匀扫到kmax线上每一点同时对应x轴和y轴方向的梯度分量。k空间采样点数Nsamp和视野FOV决定了kmax。以FOV256mm、Nsamp128为例kmax Nsamp / (2 × FOV) 128 / 512 0.25 cycles/mm对应的空间分辨率是1/(2×kmax) 2mm正好等于FOV/Nsamp。这个关系必须成立否则重建图像的空间尺度和像素间距会对不上。角度生成方面初期我用均匀分布角度Nproj 180; angles (0:Nproj-1) * pi / Nproj;后来我也试过黄金角增量每次递增约111.25°。黄金角的优势在于任意连续数据块里的投影角度都近似均匀适合动态成像和实时重建。但稳态仿真第一步用均匀角度就够探索径向轨迹的统计学特性时再上黄金角。采样点k坐标生成Nsamp 128; dk 1 / FOV; kvals (-Nsamp/2 : Nsamp/2 - 1) * dk;这里kval的起点是负值、终点是正值k0正好落在数组中间。这个对称性在后面做ifftshift时非常关键。3.3 编码相位的符号约定——最容易错的一步有了k坐标之后每个体素在采集时积累的空间编码相位为phase 2π × k(s) × (x·cosθ y·sinθ)如果用矩阵化的写法在Matlab里是这样的% 第p条投影角度为theta for s 1:Nsamp phase 2 * pi * kvals(s) * (X * cos(theta) Y * sin(theta)); signal(s) sum((Mx 1i*My) .* exp(-1i * phase), all); end这里的指数符号是负号也就是约定信号s(t) ∫Mxy(r,t)·exp(-i2πk(t)·r)dr。这个约定不能随便更改因为后面的B0偏移相位、重建时的共轭操作全部要跟它保持同一符号体系否则重建图像会左右翻转。还有一个我直到现在都会条件反射检查的坑X和Y矩阵的坐标原点。我用的linspace(-FOV/2, FOV/2-FOV/NX, NX)会把坐标原点放在网格边缘的第一个像素和FOV中心之间偏离半个像素的位置。这个半像素偏移在绝大多数情况下不影响图像视觉效果但在做像素级对比时会造成整体平移。如果追求像素级对齐应该用x (-NX/2 : NX/2-1) * (FOV/NX);这个区别在对称phantom的自检里很容易暴露。后面第6章会教你怎么通过对称性检验发现这个问题。4. FLASH序列参数设定翻转角、TR/TE与对比度的权衡逻辑序列参数是布洛赫仿真的灵魂。同样的物理仿真代码参数不同结论完全不同。这一章不讲教科书的全部推导只讲我实际仿真FLASH时怎么定参数以及为什么这样定。4.1 低翻转角为什么是FLASH的灵魂FLASH从一开始就放弃了90°激励。核心原因是TR很短时90°激励后纵向磁化来不及恢复下一次激励能用的信号就很少。低翻转角把一部分纵向磁化保留下来换来的是每个TR都能有相对稳定的信号。对于稳态梯度回波带扰相稳态纵向磁化和信号幅度可以写成Mz_ss M0·(1-E1) / (1-cosα·E1) S ∝ sinα·Mz_ss·exp(-TE/T2*)其中E1 exp(-TR/T1)。让信号对翻转角求导等于零得到Ernst角αE arccos(exp(-TR/T1))这个公式的重要性在于它告诉你翻转角不是越大越好而是由TR和T1共同决定。举个例子T1700ms、TR15ms时E1 exp(-15/700) ≈ 0.979 αE arccos(0.979) ≈ 11.8°所以我在仿真里把翻转角设在12°既能获得最大稳态信号又保留足够的T1对比度。如果你把翻转角改成25°信号幅度不会更高反而会大幅降低纵向稳态磁化图像整体变暗。这是FLASH这个核最核心的机制。4.2 我的初始参数表我实际用的一组仿真参数如下参数数值说明网格尺寸64×64调试快重建也能看出问题FOV256mm对应像素分辨率4mm够用采样点数Nsamp128单条k线采样数投影数Nproj1800~179°间隔1°TR15ms稳态梯度回波的典型短TRTE6ms落在读取窗中心翻转角12°按Ernst角估算中心组织T1/T2600ms / 70ms模拟高信号实质外围组织T1/T21000ms / 120ms长T1对比dB0偏移0Hz单等色严格T2*见6.1节TE的选择逻辑是梯度回波没有180°重聚脉冲TE越短T2*衰减越少信号越强但TE也不能太短否则磁敏差异造成的相位演化还没来得及体现。6ms在大部分组织里是合理的折中。如果你想看明显的磁敏效应把TE拉长到15~20ms图像里就能看到组织交界处的信号流失。4.3 如何确认模拟进入了稳态FLASH是稳定序列但磁化矢量不是从第一个TR就立刻稳定的。初始时刻MzM0第一次12°激励后的信号会显著高于稳态。如果不处理这个瞬态采集到的第一条投影线幅度异常重建后会出现沿径向方向的条带伪影。处理办法有两个。我的做法是在正式采集前先跑一段热身循环warm-up让纵向磁化收敛到稳态warmupTRs 20; for idx 1:warmupTRs [Mx, My, Mz] applyExcitation(Mx, My, Mz, alpha); [Mx, My, Mz] applyEvolution(Mx, My, Mz, TE, T1, T2, dB0); % 理想扰相清除残余横向磁化 Mx(:,:) 0; My(:,:) 0; [Mx, My, Mz] applyEvolution(Mx, My, Mz, TR-TE, T1, T2, dB0); end每一次都清零Mx和My对应FLASH序列末端的扰相梯度或者射频扰相。如果你要严格模拟无扰相版本比如真正的bSSFP序列就不要清零但稳态条件和信号公式完全不同。稳态是否收敛可以看单个体素的Mz变化幅度。我写过一个脚本记录每次TR结束后的Mz平均值当相邻两次TR的相对变化小于1e-3时就认为进入了稳态。通常20个TR足够TR越短收敛越快。5. 重建路径与验证从径向k空间线到图像仿真采集到的数据是一组原始k空间线形状是Nproj×Nsamp的复数矩阵。接下来要做的是把它变成一张能肉眼判断好坏的图像。这一步如果不走对前面的布洛赫物理全部白搭。5.1 方法Aifft加滤波反投影——最稳径向k空间和图像之间的桥梁是中心切片定理每个角度的k空间线做一维逆傅里叶变换得到的正是该角度下物体的辐射投影。因此最直接的路径是先把每条k线变换为投影数据再用滤波反投影iradon重建。Matlab里实现很简洁% 每条k线先做ifftshift再ifft再ifftshift确保k0对齐到逆傅里叶中心 proj ifftshift(Kdata, 2); proj ifft(proj, [], 2); proj ifftshift(proj, 2); % iradon需要角度单位为度数且投影按探测器维度排列 angles_deg angles * 180 / pi; recon iradon(proj., angles_deg, linear, Ram-Lak);iradon默认内置了Ram-Lak斜坡滤波器正好补偿径向k空间中心过采样造成的低频权重过大。如果你自己写网格化重建这步滤波是无论如何也绕不开的。我第一次跑这个流程时重建出来的图像左右翻转且略有一像素偏移。后来发现两个原因一是ifftshift的方向和频域原点没对齐二是iradon内部对探测器原点的定义和我的网格定义差了半个像素。这类问题不用硬记直接用一个中心对称的phantom做测试看到翻转就翻转回来看到平移就调整ifftshift或对k线补零。5.2 方法B笛卡尔重采样加密度补偿另一种做法是把极坐标下的k线插值到笛卡尔网格再对网格做二维ifft。这种方法的优势是更接近MRI重建的主流实现gridding但需要密度补偿函数否则重建图会有明显的模糊和边缘振铃。密度补偿的基本逻辑是k空间中心区域每条线贡献的采样点过多权重必须降低。理想的权重是采样点周围区域的面积倒数即1/|k|越靠近边缘权重越大。我叫它越稀疏越重要。实现时可以用三角插值算出每个笛卡尔网格点周围受哪些极坐标样本影响然后按1/r归一化。但如果你既不想装JSON工具箱又不想手写复杂的网格化器方法A的iradon在这个仿真规模下完全够用。iradon自带滤波效果不需要显式做密度补偿。所以项目主线我推荐方法A方法B适合你想在重建算法上做文章的后续扩展。5.3 投影数量和星芒伪影的压力测试径向采集的伪影很有辨识度当投影数不够时重建图像会从高强组织向外辐射出尖刺状条纹。原因是径向角度间隔过大k空间外围各向异性采样造成的角域混叠。我做了一组对比仿真网格64×64、采样128点/线投影数分别是60、120、180、256。结果和经验公式基本吻合要避免可见星芒投影数需满足Nproj ≥ π × Nsamp / 2128采样点时π×128/2≈201。所以180条线会有轻微条纹256条才几乎看不见。这就解释了为什么不少径向MRI应用用黄金角或者高角度采样不只是为了动态时间分辨率更是为了压制角方向混叠。6. 调这套仿真时踩过的坑与自检方法最后这部分是纯干货清单。每一段都是我真真实实调试过的问题按排查顺序写希望能帮你少走弯路。6.1 单isochromat模型根本没有T2*信号衰减这是整套仿真里最大的一个物理陷阱。许多教程会告诉你在布洛赫模拟中引入一个B0偏移场dB0就能看到信号衰减模拟出T2*效果。这句话严格来说是错的。如果每个体素只有一根isochromat那么这支磁化矢量在B0偏移下只是整体旋转了一个相位模长完全不变信号幅度不会衰减。真正的T2*衰减来自一个体素内成千上万个自旋它们各自处在不同的微观场环境里进动频率不同随时间推移相位散开宏观信号相干性才逐渐消失。解决这个问题有三条路线。第一条严格路线在每个网格点内部再放几个子isochromat给它们各自的随机dB0偏移把所有子isochromat的信号叠加作为该体素输出。第二条工程近似路线直接给弛豫项换成exp(-dt/T2*)人为加入T2*衰减。第三条折中路线给每个体素分配一个均值附近的随机dB0分布比如高斯分布标准差5Hz体素间保留一定空间平滑性。我在FLASH仿真里用的是第二条加第一条的组合先按T2*简化跑通主流程再做一组子isochromat版本的对照实验验证趋势没变。如果你想严谨强烈建议用第一条因为它能自然模拟出梯度回波随TE衰减的曲线。6.2 单位、相位符号与采样时钟的一致性我在代码里处处强调单位因为它实在太容易出错了。dB0单位是Hz相位公式里要乘时间时间为毫秒时必须换算成秒再乘。t6ms、dB05Hz时相位 2π × 5 × 0.006 ≈ 0.188 rad如果忘了除以1000算出来的是188 rad完全没有物理意义信号会乱闪。相位符号也是如此。空间编码用exp(-i2πk·r)B0偏移相位就用同符号exp(-i2πdB0·t)。反向使用会导致k空间数据在频域里左右颠倒重建图像跟着左右颠倒。这一类错误很难肉眼发现只能靠对称性测试。还有一个隐蔽问题采样时钟和k空间坐标的关系。每条k线并非瞬时采完而是在读出窗内随时间匀速移动的。如果TE不落在采样窗正中心k0对应的就不应该是信号最大值位置而是偏在一边。FLASH的回波定义是TE时刻信号峰值所在的采样点所以要确保采样点序号和TE时刻对齐否则重建相位会错乱。6.3 用对称phantom做自检与参数扫描的收获我强烈建议你在布洛赫仿真的早期就准备一个关于原点对称的phantom。它可以是中心圆加四个对称小方块只要左右翻转后和原图一样就行。然后用仿真重建结果检查翻转后的图像是否等于原图像。如果不等于先查相位符号再查ifftshift方向再查网格原点定义。做完功能验证后我还做了一组翻转角扫描α分别取5°、10°、12°、15°、25°其他参数不变。重建两张图的相对信号强度结论落在12°附近表现最好和Ernst角公式预测的一致。这件事本身没什么科研价值但能直观说明仿真确实捕捉到了FLASH的物理本质——它不是在画画而是在复现一种真实的信号机制。我最后想说的是这套仿真里最有用的不是最终图像本身而是你在推导、编码、排错过程中建立起来的参数直觉。以后扔给你一台没有调试面板的MRI扫描仪你也能在脑子里预估改变翻转角会发生什么——这种能力才是做布洛赫仿真的最大回报。