简介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形态会严重变形hrfmodelHRF模型类型canonical可选项包括canonical和fir按分析目的选择deltaHRF采样间隔0.5一般取TR的约数数值越小拟合越精细计算越慢order模型阶数3控制HRF形态复杂度太大会过拟合AR是否建模自回归噪声1建议开启能显著改善噪声模型thr显著性阈值0.05低于阈值的体素可能被判定为无显著HRF响应verbose是否显示迭代日志1排错时打开正式批处理时关掉调参最大的心得是不要为了追求“漂亮结果”而把阶数设得过高。曾经有一次我为了更贴合某个被试的HRF曲线把order提升到8结果单个被试看着拟合得很好但组水平分析时却出现了明显的离群值。后来检查发现是阶数过高导致部分体素把噪声也拟合进去了HRF形态出现不合理的多峰。现在我做组分析时基本固定order3虽然损失一些个体拟合精度但组间统计结果稳定得多。还有个细节是delta参数。它表示HRF估计的时间分辨率默认取TR的约数即可。如果TR2delta0.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加入流程试试也许能解释掉一些一直让你困惑的“噪声”来源。本文还有配套的精品资源点击获取