企业官网建设流程全解析

简介：这是一份用于核磁共振（NMR）自由感应衰减信号模拟的MATLAB代码包，围绕Bloch方程数值求解展开，面向需要理解核磁矩演化、T2弛豫及FID信号规律的科研人员与学生。压缩包共6个文件，其中5个为.m脚本/函数，1个为txt结果记录，整体大小仅2KB；脚本内置四阶Runge-Kutta数值积分方法，可稳定求解描述磁矢量变化的微分方程组，另有多个变体脚本对应不同初始条件或脉冲序列设定，txt文件则保存了三种输入情况下的FID数据。已有202人学习下载。通过学习这份代码，读者不仅能动态观察不同参数下FID信号的衰减行为，还可自行修改弛豫时间、射频脉冲等条件，对比自旋锁定、CPMG等不同序列的仿真结果，从而深入理解T1、T2弛豫对信号的影响。这种可运行的仿真模型对NMR实验参数优化、谱图解析以及生物组织物理化学性质研究都具有实用参考价值。

1. FID.zip 和 bloch 方程：先于图像的是一串随时长衰减的信号

把一个叫 FID.zip 的压缩包解开，里面通常不会出现任何一张磁共振图像，排在前面的大多是一段没有注释的微分方程积分程序。做过序列开发的人一眼就能认出这是在跟什么打交道：自由感应衰减（FID），也就是射频脉冲停掉以后，横向磁化矢量在 Bloch 方程约束下边进动边衰减、被接收线圈采成复数点的信号序列。FID.zip 把 Bloch 方程的求解过程从纸面公式搬成了一段能跑的数值代码，让 T1、T2、共振偏置、梯度强度这些参数直接对应到一条可重复计算的曲线上。对 MRI 参数定量、新序列仿真和 NMR 谱仪调试这几个方向来说，它比动不动就上机扫一版参考数据要顺手得多。这篇就围绕标题里的 FID.zip 与 Bloch 方程展开：先说清楚方程怎么被离散求解，再写最小可运行模拟代码，最后给出从 FID 曲线反推弛豫参数的拟合套路。

2. Bloch 方程怎么被 FID.zip 变成 FID：解析替代方案与数值步进

很多人第一次接触 Bloch 方程是在教科书上看旋转坐标系里的三个分量表达式。对自由进动段而言，解析解确实够用，但 FID.zip 这类工具的存在意义恰恰在于：真实脉冲序列里有 B1 不均匀、梯度爬升、共振偏置漂移，这些会让解析解的适用条件逐个失效。要理解代码在算什么，先得把方程的每一项对号入座。

2.1 旋转坐标系下的自由进动：FID 的本质是横向分量的相干衰减

磁化矢量 M 在静态磁场 B0 下以拉莫尔频率围绕 B0 进动。射频激发后，M 偏离纵向，横向分量 Mxy 在接收线圈中感生出电压，这就是 FID 的原始物理来源。设旋转坐标系以拉莫尔频率转动，则零阶近似下，Bloch 方程可写成：

dMx/dt = γ·My·ΔBz − Mx/T2
dMy/dt = −γ·Mx·ΔBz + γ·B1(t)·Mz − My/T2
dMz/dt = −γ·B1(t)·My − (Mz − M0)/T1

其中 ΔBz 是局域场与 B0 的偏差，单位为特斯拉；γ 是旋磁比；B1(t) 是射频场包络。方程右侧第一项描述进动，第二项描述横向弛豫，纵向方程里则是指数趋近 M0 的恢复过程。FID 信号通常取 s(t) = Mx(t) + i·My(t)，因为接收链路是正交解调，实部和虚部一起构成复数谱。

如果整个过程中没有梯度、没有 B1 残余、共振偏置固定，那么这段方程有闭式解，比如 90° 脉冲后的理想 FID 信号就是 s(t) = M0·exp(−t/T2*)·exp(i·2π·f0·t)。但 FID.zip 没有一个完整的理想 RF 脉冲段需要模拟的临床应用。激发期间 B1(t) 不是瞬时完成的，选层梯度也会让不同空间位置的 ΔBz 不同，此时解析积分就需要分段拼接，误差会沿着每个不连续点累积。数值方法通过把时间轴切成小段，在每段内用当前磁场和弛豫状态更新 M，天然绕过了分段拼接的复杂性。

2.2 解析解不够用的三个场景：时变梯度、非均匀 B1 与有限翻转角

第一个失效场景是时变梯度场。梯度波形在脉冲序列里往往是梯形或正弦形，爬升期内的磁场随时间和空间同时变化，空间上不同体素经历的进动频率不同。用解析解处理时，需要按空间位置逐个积分，而且每个时间步的 B 值都不同，闭式解的推导工作量随序列复杂度指数上升。

第二个场景是非均匀 B1。射频线圈的 B1 场在成像体积内有空间分布，实际翻转角不是设定值 α，而是 α(r) = γ·B1(r)·τ。理想脉冲模型假设翻转角处处相同，但 Bloch 方程数值求解只需把每个体素各自的 B1 代入驱动项，不需要为“非理想”单独构造理论模型。

第三个场景是有限翻转角下的激发过程。对于小翻转角可以线性近似，但在 90° 或 180° 脉冲里，Mz 和 Mxy 之间的能量交换高度非线性。若把脉冲当作一个瞬时旋转矩阵，就会丢失脉冲形状带来的过渡响应。数值步进则是把射频脉冲切成数百个微小时间片，每片更新一次 M，脉冲形状、幅度、相位调制都天然包含在内。

2.3 FID.zip 类求解器的最小积分步：四阶龙格库塔与参数对照

常见做法是把进动项和弛豫项合在一个 ODE 右侧，再用四阶龙格库塔（RK4）积分。RK4 每步计算四个斜率，对 Bloch 方程这种维度低但可能存在高进动频率的系统，其离散化误差远小于显式欧拉法。一个可直接替换到 FID.zip 里的步进函数如下：

# 单步 Bloch 方程更新，采用 RK4，M 为 [Mx, My, Mz] def bloch_step(M, wx, wy, wz, T1, T2, dt): # wx, wy：有效横向进动角频率，包含梯度与共振偏移 # wz：B1 导致的旋转项，无射频时置为 0 def rhs(m): Mx, My, Mz = m dMx = wz*My - Mx/T2 dMy = -wz*Mx + wx*Mz - My/T2 dMz = -wx*My - (Mz - 1.0)/T1 return np.array([dMx, dMy, dMz]) k1 = rhs(M) k2 = rhs(M + 0.5*dt*k1) k3 = rhs(M + 0.5*dt*k2) k4 = rhs(M + dt*k3) return M + (dt/6.0)*(k1 + 2*k2 + 2*k3 + k4)

这段代码里 wx 对应 B1 场的旋转项，wz 对应磁化矢量绕 z 轴的进动。实际使用时要特别注意单位：M0 归一化为 1，时间单位是秒，频率单位是弧度/秒。dt 的选取准则是小于最高进动频率周期的 1/10。比如共振偏置为 1 kHz 时，周期是 1 ms，dt 应小于 0.1 ms；留有余量时通常取 1/20 周期。步长过大，信号尾部会出现相位抖动甚至幅度发散；步长过小，则会显著拖慢多体素模拟的速度。

下面参数表列出 FID.zip 里最常调整的量以及它们对输出 FID 的影响：

3. 用 FID.zip 的最小代码结构跑通一段 FID：参数装配、脉冲与读出窗

拿到 FID.zip 这类模拟器时，我一般先把它的输入输出结构拆成三块：序列参数文件、积分主循环、输出结果。序列参数文件描述射频脉冲的翻转角、时间点、梯度波形；积分主循环把 Bloch 方程随时间步进；输出结果则是每一时刻的 Mx 和 My。搞清楚这三个块的边界之后，再跑最小示例就非常直接。

3.1 最小序列配置：为什么要用归一化单位和绝对时间轴

常见做法是把所有时间参数都放在同一个时间轴上，射频脉冲用起始时间和持续时间限定，梯度用阶梯数组表示，读出窗单独开一段。用归一化 M0 的好处是计算结果不依赖具体磁体场强，T1、T2 也能以秒为单位直接写。这样一个 90° 脉冲加读出窗的最小配置可以写成：

import numpy as np # 配置段 T1 = 0.8 # 纵向弛豫时间，单位 s T2 = 0.06 # 横向弛豫时间，单位 s f0 = 300.0 # 共振偏置，单位 Hz dt = 1e-6 # 积分步长，1 us t_obs = 0.2 # 读出总时长，单位 s M = np.array([0.0, 0.0, 1.0]) # 理想 90° 脉冲：沿 x 轴旋转 90°，磁化矢量落到 y 轴 M = np.array([0.0, 1.0, 0.0]) t = np.arange(0.0, t_obs, dt) signal = np.zeros(len(t), dtype=np.complex128)

这个配置里先把 90° 脉冲当作瞬时完成，激发完成后直接进入自由进动。对大多数参数定量实验，这种近似在脉冲持续时间远小于 T2 时成立。若要模拟实际有限脉冲，就需要在循环里加入 B1 分段。

3.2 主循环：把 Bloch 方程推进和信号采集合在一起

主循环里每个时间步调用一次 bloch_step，然后取 Mx + i·My 作为该时刻的 FID 信号。读出窗的频率响应等效于接收带宽，采样步长 dt 决定 Nyquist 频率：1/(2·dt)。f0 超过这个上限就会在频谱上折叠回低频率，所以 dt 的选取同时受进动频率和期望谱宽约束。

for i, ti in enumerate(t): # 自由进动过程：无 B1 驱动，只有进动和弛豫 wx = 0.0 wy = 0.0 wz = 2.0 * np.pi * f0 M = bloch_step(M, wx, wy, wz, T1, T2, dt) signal[i] = M[0] + 1j * M[1]

这一步完成后，signal 就是一个复数数组。对 signal 做 FFT 可得到频域谱线：峰值出现在 f0 处，谱线宽度约等于 1/(π·T2*)。若要模拟磁场不均匀导致的附加衰减，把 wz 改成包含 ΔB0 的形式即可，不需要改动步进函数本体。

3.3 多体素模拟：把梯度效应装进每个自旋的进动频率

做梯度实验时，空间位置 r 处感受到的磁场是 B0 + G(t)·r。因为 Bloch 方程是对单个自旋积分，多体素模拟就是对每个位置分别调用同一套步进逻辑。假设一个一维样本跨 64 个体素，梯度强度 G 在读出窗期间恒定，那么每个体素的进动频率是 γ·G·r。信号累加时把各体素贡献的复数信号求和。

n_voxel = 64 gamma = 42.58e6 # Hz/T G = 0.005 # mT/m 转换为 T/m 后使用 positions = np.linspace(-0.05, 0.05, n_voxel) # 单位 m freq_offset = gamma * G * positions # 每个体素对应的频率偏移 total_signal = np.zeros(len(t), dtype=np.complex128) for r_idx, df in enumerate(freq_offset): Mv = np.array([0.0, 1.0, 0.0]) for i, ti in enumerate(t): wz = 2.0 * np.pi * (f0 + df) Mv = bloch_step(Mv, 0.0, 0.0, wz, T1, T2, dt) total_signal[i] += Mv[0] + 1j * Mv[1]

这段循环里对每个体素重复积分，精度有保证，但代价是计算量随体素数线性增加。更高效的做法是用矩阵或 GPU 批量更新，不过那属于优化阶段。初学者先跑通这个逐体素版本，再考虑向量化。

4. 从 FID 输出反推 T2* 与 T1：拟合代码、窗宽选择和 T2/T2* 的边界

模拟本身不是终点。FID.zip 的典型用途是把模拟曲线当成“已知答案的合成数据”，用来验证拟合算法。这里最容易踩的坑是把 T2* 当成 T2 直接解读，二者相差一个与 B0 不均匀性相关的附加项。

4.1 包络拟合 T2*：先取模值还是直接用复数

T2* 的估计通常基于包络衰减。取模值的好处是不受共振偏置和初始相位影响，缺点是当噪声接近信号幅度时，模值平方会引入 Rician 偏置，尾部拟合出虚高的 T2*。简单做法是取信号的实部或复数拟合，把模型写为：

s(t) = S0·exp(−t/T2*)·exp(i·(2π·f0·t + φ))

先去掉激发刚结束的一小段数据，因为那一段受有限脉冲和数字滤波过渡影响严重。最小二乘拟合中固定 f0 和 φ 可以让 T2* 估计更稳。用 Python 直接线性化拟合包络：

from numpy.polynomial import polynomial as P cut = int(0.01 / dt) # 去掉前 10 ms，约等于 RF 脉冲过渡 t_fit = t[cut:] amp = np.abs(signal[cut:]) amp = amp[amp > 1e-3] # 剔除接近零的尾部 t_fit = t_fit[:len(amp)] # 线性拟合 ln(amp) = ln(S0) - t/T2* coeffs = np.polyfit(t_fit[:len(amp)], np.log(amp), 1) T2_star = -1.0 / coeffs[0]

用 np.polyfit 拟合对数线性模型计算快，但会放大尾部噪声。想更稳妥就改用 scipy 的 curve_fit 直接拟合指数模型，并设置初始值。线性拟合版本适合快速预检，最终报告建议用非线性拟合。

4.2 为什么 FID 直接量到的是 T2*，而不是 T2

没有梯度、没有体素内频散时 FID 衰减由 T2 决定；实际样品总存在局部磁场差异，体素内不同自旋的进动频率略有不同，横向磁化逐渐失去相位一致性，信号更快衰减。用 1/T2* = 1/T2 + γ·ΔB0 近似描述，其中 ΔB0 是体素内磁场分布的标准差。因此 FID 拟合得到的指数时间常数天然小于等于真 T2。

要精确获取 T2，需要在 180° 脉冲构成的回波峰值处采样，即自旋回波或 CPMG 序列。FID.zip 模拟中也可以通过加一段 180° 脉冲和回波窗来验证：主循环里在时刻 τ 处执行翻转 180°，再继续推进到 2τ。回波峰值处的幅度拟合出的才是纯 T2。参数定量项目里，这个区别如果不在读数据前理清，实验结果会系统性偏高。

4.3 拟合窗宽、基线和数据裁剪的 3 个边界

第一个边界是拟合起点必须避开 RF 脉冲脱落和接收器相位过渡段，常见做法是弃置 5～10 倍脉冲持续时间。第二个边界是终点不能拉太长，信号衰减到噪声底附近后，对数线性拟合的残差会被噪声主导，此时多出的数据点不会提高精度，反而拉偏斜率。经验值是把拟合窗设为 2～3 倍预估 T2*。第三个边界是基线，实际信号会有直流偏置，取模或者对数前先减掉基线均值：

baseline = np.mean(amp[-int(0.01/dt):]) # 用最后 10 ms 估计噪声底 amp_corrected = np.maximum(amp - baseline, 0)

基线估计区间要选在信号已经完全衰减之后，若 T2* 为 50 ms，则观察窗至少要 250 ms 以上。

下表列出拟合中常见现象与对应处理方向：

T1 拟合则需要改变重复时间或采用饱和恢复序列。最小实现是在每次激发前把 Mz 重置为 0，然后在不同恢复时间读取信号。FID.zip 里改的是主循环之外的恢复等待段，将 Mz 指数趋近 M0 的过程与 FID 读出一段分离，能显著降低模拟代码复杂度。

5. 把 FID.zip 当作序列预验证器：梯度项、有限脉冲与三个自检方法

原理和基础拟合都跑通之后，FID.zip 的真正价值在于缩短序列迭代周期。与其在扫描仪上反复测试梯度爬升和翻转角误差，不如先在模拟工具里把时序跑一遍，把明显的问题暴露出来。

5.1 按时间轴分段施加 B1 和梯度，序列逻辑更接近真实扫描

把脉冲序列抽象成时间轴上的事件列表。每个事件有类型和持续时长。射频脉冲事件内 wx 不为 0，梯度事件内 wz 增加一处随位置变化的项。B1 不均匀模拟则在射频段乘一个空间调制系数。这样写出的模拟器已经不是单纯的 FID 计算器，而是可以验证梯度回波、自旋回波等序列的通用平台。

# 事件表驱动模拟伪结构，单个事件内更新参数 for event in sequence: if event.kind == 'rf': wx = gamma * event.b1 # b1 为包络幅度 wy = 0.0 elif event.kind == 'gradient': wx = 0.0 # wz 额外叠加 gamma * G * r elif event.kind == 'readout': signal_buffer.append(M[0] + 1j*M[1]) # 每个 dt 调用 bloch_step 推进

把事件和物理参数分离之后，后期替换脉冲形状或梯度波形都不需要改动积分函数。

5.2 三个必做的自检：解析对照、翻转角扫描和回波时间扫描

第一个自检是让模拟回到解析解。关掉 B1 不均匀和梯度，取 f0 = 0，拟合模拟 FID 的幅度，得到的时间常数应当接近设定的 T2。偏差若超过 1%–2%，基本可以断定是 dt 过大或拟合窗处理不当。

第二个自检是翻转角扫描。固定序列，把 90° 脉冲换成不同翻转角 α，重复多次模拟后测量初始信号幅度。理想情况下幅度应正比于 sin(α)。若曲线偏离正弦，说明 B1 不均匀或脉冲持续时间较长时间内弛豫不可忽略。这项检查直接暴露射频校准偏差。

第三个自检是回波时间扫描。在自旋回波模拟中改变 TE，记录回波峰值幅度，绘制 ln(S) 对 TE 的散点。如果 CPMG 模拟数据落在了直线上，说明 T2 拟合前提成立；若呈上凸或下凹，代表存在扩散效应或射频误差。

5.3 步长、存储与体素数的取舍留给后续优化

多体素多帧模拟的数据量增长比想象快。一次 256×256 体素、T2* 为 30 ms、采样步长 1 μs 的 FID 模拟，输出矩阵将超过几 GB。常见做法是在循环内直接以固定间隔取点存储，而不是保存全部时间步；代表性地每 100 μs 存一个复数点，已经能够保留完整的谱信息。若需要频谱分析，再在离线阶段补零插值。FID.zip 的定位从一开始就是物理验证工具，跟成像重建的高吞吐数据流并非同一类任务，所以优先保证物理量正确，其次再优化存储和批量计算。

本文还有配套的精品资源，点击获取