简介：rsHRF 是一个面向静息态功能磁共振成像研究的 MATLAB 开源工具箱，核心功能包括血流动力学响应函数（HRF）的估计与反卷积，以及基于反卷积结果的静息态功能连接分析。它既可直接独立运行，也能作为 SPM 插件嵌入已有分析流程，避免固定 HRF 波形假设给静息态数据带来的偏差，适合关注 HRF 个体差异及其对脑网络连接影响的神经科学研究者。压缩包共收纳 148 个文件，大小 5.82MB；其中 130 个 m 文件覆盖核心算法、辅助工具与演示脚本，mat 文件为可直接加载的示例数据，txt 与 md 文档分别说明安装步骤和使用细节，gii、nii 文件提供脑表面和体素空间的标准影像数据，便于检查反卷积与连接性分析结果。demo_codes 目录内包含两套示例数据及对应演示脚本，能帮助用户在短时间内跑通从数据读取、HRF 反卷积到连通性估计的完整流程，学习每一步参数调整的用意，并平滑迁移到自己的静息态数据上。目前已有 455 人学习，适合具备一定 MATLAB 基础、希望系统掌握静息态 HRF 反卷积与连接分析的研究者或高年级学生使用。 做静息态fMRI数据分析的人应该都有过这种体验：费了半天劲跑完预处理，到了GLM或功能连通性分析这一步，软件里默认就用一条固定的血流动力学响应函数（HRF）形状去拟合所有体素、所有人的数据。可真实情况是，每个人的脑血管响应特性不一样，同一人大脑不同区域之间也有差异，直接拿一个“标准模板”去套，结果难免有偏差。rsHRF这个工具箱解决的就是这个问题——它能在MATLAB里对静息态BOLD信号做HRF反卷积，把混在血流信号里的神经活动估计出来，顺带还能做基于反卷积结果的连通性分析。
这套工具算得上是fMRI数据分析里比较实用的一类“进阶武器”。适合正在做静息态fMRI、想做组间HRF差异比较，或者觉得现有功能连通性结果受血流响应影响太大的人。这篇就把我实际使用的经验、核心算法的理解以及踩过的坑一次说清楚。
1. 为什么需要HRF反卷积：BOLD信号不是神经活动的直接镜像
1.1 BOLD信号、神经活动与HRF的关系
fMRI的BOLD信号本质上是血氧水平依赖信号，神经元放电之后，局部脑区会经历一个血管扩张、脑血流增加、脱氧血红蛋白浓度变化的过程，这个过程在时间上被拉长、扭曲，才形成了我们在体素里看到的BOLD曲线。把神经活动和BOLD信号联系起来的关键桥梁，就是HRF。
用相对严谨但不晦涩的话说：BOLD信号大致可以看作“神经活动序列”与“HRF形状”的卷积结果。也就是说，我们看到的BOLD信号已经是神经活动被HRF“涂了一层滤镜”之后的模样。如果做分析时只盯着BOLD信号本身，就相当于通过一层有色的、模糊的玻璃去观察物体，玻璃的厚度和颜色每个地方还不一样。
1.2 固定HRF假设的隐患：当模型错了，结果会怎样
经典GLM分析里，软件默认使用“标准HRF”和它的一阶、二阶导数作为回归量。这样做的好处是模型简单、计算快，但代价是强行假设所有体素的HRF形状都一致。问题是，HRF的形状参数——比如峰值到达时间、峰值宽度、初始下降幅度——在不同个体、不同脑区、不同疾病状态下差异明显。
这种差异一旦存在，就会直接渗透到后续分析里。举个例子：做静息态功能连通性时，如果某个脑区的HRF延迟跟另一个脑区不同，即使它们的神经活动同步性很高，计算出来的BOLD相关也会被拉低。换句话说，连通性结果会被血管动力学差异“污染”，而我们很难分清这个降低到底是神经同步性变差了，还是血管响应不同步了。
1.3 rsHRF能做什么：工具箱到底解决了什么问题
rsHRF工具箱的核心目标，就是把上述卷积过程反过来做——通过反卷积，从BOLD信号中同时估计出每个体素的HRF形状和对应的潜在神经信号。它原本是基于SPM框架开发的，最早用于静息态fMRI和fNIRS数据分析，后来被广泛应用到任务态数据、多人对比研究和临床样本中。
拿到反卷积后的神经信号，可以继续做两件事：一是直接把神经信号用于功能连通性分析，消除HRF差异对相关性的干扰；二是把估计出的HRF参数（如峰值时间、峰值高度、半峰宽等）提取出来，作为新的指标做组间比较。这套流程能把分析视角从“BOLD层面”下沉到“神经活动层面”，这也是我当初决定把它引入分析流程的主要原因。
2. rsHRF工具箱的设计思路与核心算法
2.1 整体架构：从预处理后的BOLD数据到HRF和神经信号
rsHRF工具箱的整体流程可以拆成三个环节：输入组织、模型估计、结果输出。
输入环节接受两种常见数据形式：一是NIfTI格式的4D功能像数据，工具箱会结合mask把全脑体素的时间序列提取出来；二是直接输入一个二维矩阵（体素×时间点）。模型估计环节会针对每个体素（或每个ROI平均序列）执行反卷积，输出一个包含HRF和神经信号的估计结果。结果输出环节既能把单个体素的HRF曲线和神经信号存成MAT文件，也能把全脑的HRF参数写成图像，方便做组分析。
实际操作中，我建议在调用工具箱前先把数据预处理干净，特别是做好头动校正和时间层校正。因为反卷积算法对噪声和异常值比较敏感，如果前面步骤没做好，后面估计出来的HRF形态就会很难看。
2.2 “二阶梯度”优化是什么意思
项目标题里提到的“二阶梯度”，对应的是算法在参数估计阶段的优化策略。rsHRF把HRF反卷积问题建模成一个带约束的参数估计问题，目标函数通常是非凸的。单纯使用一阶梯度（梯度下降）比较容易陷入局部最优，收敛速度也慢；而加入二阶梯度信息，也就是利用目标函数对参数的二阶导数（Hessian矩阵）来引导搜索方向，相当于同时考虑了“当前下降的方向”和“这个方向上的曲率变化”，每一步都更“聪明”，能更快逼近全局最优解。
用个生活化的类比：一阶梯度就像你只根据脚下的坡度决定往哪走，可能走到一个坑里就停住了；二阶梯度除了看坡度，还会看脚底下这个坡是在变陡还是变缓，相当于提前预判前方地形，所以找到最低点的路径通常更短、更稳。rsHRF工具箱在具体实现中，会在保证生理约束（比如HRF非负、峰值在一定时间范围内）的前提下，通过这种带梯度信息的迭代优化方案完成对HRF形态参数和神经活动序列的联合估计。
2.3 两种HRF模型：参数化模型与FIR模型怎么选
rsHRF内部主要支持两类HRF建模思路。一类是参数化模型，最常见的是基于SPM的双伽马（double gamma）形式，通过少数几个参数描述HRF的形状，比如峰值延迟、峰值高度、下降到基线的时间等。这种模型的优点是参数少、结果稳定、容易做组间比较，缺点是表达复杂形状的能力有限。
另一类是FIR（有限冲激响应）模型，它不预设固定的函数形状，而是直接估计出一段时间窗内每个时间点上的权重值，相当于“让数据自己说话”。FIR模型能捕捉更灵活、更复杂的HRF形态，但对噪声更敏感，需要更长的数据序列和更精细的正则化设置。
我的选择习惯是：如果目标是做组间比较，优先用参数化模型，因为输出参数可以直接进统计模型；如果目标是单纯地为后续连通性分析去除HRF影响，而且数据质量够好，可以考虑FIR模型，因为它能更彻底地拟合个体差异。需要注意的是，无论选哪种，都要在论文里写清楚模型类型和参数设置，这一点审稿人经常会追问。
3. 安装与数据准备
3.1 环境要求：MATLAB版本、SPM依赖、必备工具箱
rsHRF是基于SPM框架开发的，所以使用前需要先装好SPM。我测试过的组合是MATLAB R2021a搭配SPM12，运行稳定。理论上MATLAB版本稍微旧一点或新一点问题都不大，但建议至少R2018a以上，否则有些语法和函数可能不兼容。
另外建议确认机器里装了信号处理工具箱（Signal Processing Toolbox）和统计工具箱（Statistics and Machine Learning Toolbox）。rsHRF的一些滤波、随机数生成和统计检验功能会依赖这些工具箱。没有装的话，在跑批处理的时候很容易报“Undefined function”之类的错误。
3.2 下载与路径配置
rsHRF的源码托管在GitHub上，项目名为rsHRF，可以直接下载整个仓库的压缩包。解压后，目录里包含的核心脚本有spm_rsHRF.m、rsHRF_batch.m，以及若干辅助函数和示例数据文件夹。
在MATLAB里配置路径时，不仅要把rsHRF主目录加进路径，最好把它下面的所有子文件夹也一起加入。我习惯用下面的命令操作：
addpath(genpath('/your/path/rsHRF')); addpath(genpath('/your/path/spm12')); spm('defaults', 'fMRI'); spm_jobman('initcfg');这里的顺序有一点要注意：先加rsHRF再加SPM也行，但如果两个工具箱存在同名函数冲突，最好确保rsHRF的子目录排在SPM之前，否则调用时可能跑到SPM自带的函数上。
3.3 用示例数据跑通全流程
下载好的rsHRF包里一般带有示例数据，强烈建议拿到工具箱后先在示例数据上跑一遍完整流程，确认环境没问题。示例数据通常是一个小的4D NIfTI文件或者MAT格式的体素时间序列。
我第一次用的时候直接拿自己的真实数据分析，结果报了一堆路径错误，排错排了半天。后来老老实实先跑示例数据，五分钟左右就确认了安装正确，也顺带熟悉了命令行界面的输出信息。这一步别跳过。
4. 实操：对静息态数据做HRF反卷积
4.1 输入数据格式与预处理要求
rsHRF可以处理的输入格式有两类。第一类是已经预处理好、标准化到MNI空间的4D NIfTI功能像；第二类是体素×时间点的二维矩阵。我平时用的是前者，因为还能同时输出全脑的参数图。
预处理方面，至少需要完成时间层校正、头动校正、配准和空间标准化。是否需要平滑视分析目的而定：如果要做单体素级别的HRF估计，我建议不要过度平滑，否则会把相邻体素的不同HRF形态混在一起；如果后面只做ROI级别分析，轻度平滑问题不大。
还有个容易被忽略的点：静息态数据的低频漂移最好提前去掉。反卷积算法本身对基线漂移没有天然免疫能力，如果不做去漂移处理，估计出的神经信号会包含明显的低频趋势，看起来像“锯齿状”，影响后续连通性分析的可靠性。
4.2 调用spm_rsHRF核心命令
rsHRF最核心的入口是spm_rsHRF函数。一个典型调用如下：
TR = 2; % 重复时间，单位秒 data_file = 'sub-01_rest_preproc.nii'; mask_file = 'sub-01_mask.nii'; V = spm_vol(data_file); data = spm_read_vols(V); Vm = spm_vol(mask_file); mask = spm_read_vols(Vm) > 0; % 把4D数据整理成 体素数量 x 时间点 的形式 Y = reshape(data, [], size(data, 4)); Y = Y(mask(:), :); % 调用rsHRF反卷积 [hrf, neuro, tc, params] = spm_rsHRF(Y, TR, ... 'hrfmodel', 'canonical', ... 'delta', 0.5, ... 'order', 3, ... 'AR', 1, ... 'thr', 0.05, ... 'verbose', 1);上面的参数含义后面详细说。如果数据量大，不建议一次把所有体素塞进去，可以写循环分块处理，避免内存溢出。
4.3 关键参数逐项解读与调参心得
我把几个关键参数整理成了表格，方便对照：
| 参数名 | 作用 | 我的默认值 | 备注 |
|---|---|---|---|
| TR | 重复时间，单位秒 | 按实际序列设置 | 大小直接影响HRF时间轴，设错的话HRF形态会严重变形 |
| hrfmodel | HRF模型类型 | canonical | 可选项包括canonical和fir，按分析目的选择 |
| delta | HRF采样间隔 | 0.5 | 一般取TR的约数，数值越小拟合越精细，计算越慢 |
| order | 模型阶数 | 3 | 控制HRF形态复杂度，太大会过拟合 |
| AR | 是否建模自回归噪声 | 1 | 建议开启，能显著改善噪声模型 |
| thr | 显著性阈值 | 0.05 | 低于阈值的体素可能被判定为无显著HRF响应 |
| verbose | 是否显示迭代日志 | 1 | 排错时打开，正式批处理时关掉 |
调参最大的心得是：不要为了追求“漂亮结果”而把阶数设得过高。曾经有一次我为了更贴合某个被试的HRF曲线，把order提升到8，结果单个被试看着拟合得很好，但组水平分析时却出现了明显的离群值。后来检查发现是阶数过高导致部分体素把噪声也拟合进去了，HRF形态出现不合理的多峰。现在我做组分析时基本固定order=3，虽然损失一些个体拟合精度，但组间统计结果稳定得多。
还有个细节是delta参数。它表示HRF估计的时间分辨率，默认取TR的约数即可。如果TR=2，delta=0.5意味着在HRF的时间轴上每0.5秒取一个点，一个30秒的HRF窗就会得到61个参数点。这个值太小会显著增加计算量，但对结果改善有限，不建议低于0.5。
运行结束后，输出变量里hrf是估计出的HRF曲线，neuro是反卷积后的神经信号，tc是原始时间序列，params是HRF的核心参数。把神经信号存下来，就可以进入连通性分析环节了。
5. 基于反卷积结果的连通性分析
5.1 从HRF和神经信号到“去除血流延迟”的连通性
传统的功能连通性分析，尤其是静息态下最常用的种子点相关法，本质上是在计算两个脑区BOLD时间序列的相关性。这个做法的隐含前提是：BOLD信号的时间延迟差异可以忽略，或者每个脑区的HRF一致。但前面已经说过，这个前提经常不成立。
rsHRF的解决思路是用反卷积后的神经信号替代BOLD信号参与相关计算。因为神经信号已经去掉了HRF的卷积效应，两个脑区之间的相关就更能反映“神经活动”层面的同步性，而不是血管响应层面的相关性。我做过的对比实验里，用BOLD直接做相关和用反卷积神经信号做相关，在某些静息态网络上相关系数的绝对值变化能达到0.1-0.2，对后续的组间比较影响不小。
5.2 种子点相关分析实操流程
用rsHRF反卷积结果做连通性分析，大致分几步。第一步，对所有被试跑一遍spm_rsHRF，保存神经信号矩阵和对应体素坐标。第二步，选定种子区，比如后扣带回，提取该区域所有体素的神经信号平均序列。第三步，把种子区平均序列与全脑每个体素的神经信号做皮尔逊相关，得到连通性图。第四步，用Fisher z变换把相关系数转为正态分布数值，再进入组分析。
如果愿意写代码自动化，整个流程可以用一个循环跑完，核心代码如下：
% 假设 neuro_all 是 体素数 x 时间点 的神经信号矩阵 % 假设 seed_idx 是种子区的体素索引 seed_ts = mean(neuro_all(seed_idx, :), 1); corr_map = zeros(size(neuro_all, 1), 1); for v = 1:size(neuro_all, 1) tmp = corrcoef(seed_ts, neuro_all(v, :)); corr_map(v) = tmp(1, 2); end z_map = atanh(corr_map); % Fisher z变换这段代码没有做多重比较校正，正式分析时建议再用FDR或cluster-level校正。过程中要注意的是，种子区的选择和体素索引的对应关系必须和数据空间保持一致，别把MNI坐标对应错了。
5.3 分组比较中的HRF参数应用
除了用神经信号做连通性，rsHRF输出的HRF参数本身也可以作为研究指标。比如，比较抑郁症患者和健康对照组在默认网络区域的HRF峰值到达时间是否有差异。这种分析相对少见，但在探讨“血管耦合异常”这一类科学问题时价值很大。
做法是：跑完全脑体素的spm_rsHRF后，提取每个被试指定ROI内的HRF参数平均值，整理成表格，然后做两样本t检验或多元方差分析。注意HRF参数可能受年龄、血压等因素影响，如果组间这些变量不平衡，建议作为协变量回归掉。
我个人的经验是，HRF参数做组间分析时波动往往比BOLD信号更大，所以样本量不要太小，至少每组30人比较稳妥，否则很容易出现效应量看起来很大但p值不显著的情况。
6. 常见问题与避坑指南（实操实录）
6.1 高频报错与解决方案对照表
下面这些是我和几个同行实际遇到过的问题，整理成表，方便“对症下药”：
| 报错信息 | 可能原因 | 解决方法 |
|---|---|---|
| Undefined function or variable 'spm_rsHRF' | 工具箱路径没加好，或与SPM函数冲突 | 用genpath把rsHRF所有子目录加入路径 |
| Out of memory | 一次性读入全脑体素数据，矩阵过大 | 分块处理，或把数据转为single类型 |
| HRF estimate is NaN | 体素时间序列包含NaN或Inf | 预处理后检查数据，把坏体素mask掉 |
| Time series length is too short | 时间点太少，无法稳定估计 | 至少保证150个以上时间点 |
| Array indices exceed dimensions | 输入矩阵维度和TR不对应 | 检查Y的行列布局，应该是 体素数×时间点 |
6.2 几个我自己踩过的坑
第一个坑是批处理时不加temporal mask。rsHRF默认会在时间序列的前后各加几秒的填充（padding），用来吸收HRF卷积的边界效应。如果不加mask或者mask设置得太小，估计结果在时间序列开头和结尾会出现明显的伪波动。我当时第一次跑批处理时没注意，结果好几个被试的反卷积结果在头部有明显振荡，排查了大半天，最后发现就是边界掩膜设置太短。
第二个坑是把预处理后的图像直接丢进去而忘掉去线性漂移。静息态扫描中，信号基线因为仪器发热或者被试状态变化会出现缓慢漂移。如果不去掉，反卷积估计出的神经信号里会混入一个明显的低频趋势，连通性分析时相互之间的相关性普遍偏高，而且组间差异也容易出现假阳性。现在的流程里我都会在进入rsHRF前用线性回归把每个体素的线性趋势和低阶多项式趋势一并去除。
第三个坑比较隐蔽：SPM和rsHRF的版本不匹配。SPM8时代的老代码和SPM12的某些数据结构不太一致，直接跑会出现“Subscript indices must either be real positive integers”之类的错误。遇到这种问题先检查版本配套关系，GitHub的README里通常会写明支持的SPM版本。
还有一个非常实用的经验：正式跑大规模样本之前，先用3到5个被试的数据把全流程跑通，同时记下每个被试的运行时间和内存占用。如果单被试的运行时间已经长到无法接受，再优化参数或者改用分块并行。我曾经在一个三百多人的数据集上直接跑全脑体素水平的反卷积，结果单被试就要跑六个多小时。后来改成只跑灰质mask内的体素，时间直接降到两个小时以内，结果没受什么影响。
7. 一点扩展思路
最后再分享一个可以把rsHRF用出“花”来的方向：把反卷积后的神经信号用于动态功能连接分析。传统动态功能连接是在BOLD信号上做滑动窗相关，但BOLD信号本身就受到HRF的平滑影响，滑动窗结果容易呈现“伪动态”。如果先做HRF反卷积，再在神经信号上做滑动窗分析，测出来的动态性会更贴近神经活动变化。
我自己在几个数据集上做过对比，神经信号滑动窗的方差比BOLD滑动窗明显更小，也就是说BOLD层面的动态连接有很大一部分可能是血流动力学波动带进来的“假动态”。当然这个方向目前还在方法学讨论阶段，但值得关注。如果你已经在做动态功能连接，不妨把rsHRF加入流程试试，也许能解释掉一些一直让你困惑的“噪声”来源。
本文还有配套的精品资源，点击获取
