尧图网站建设 
近红外脑功能成像fNIRS是一种非侵入式、可穿戴的脑功能活动监测技术它通过测量大脑皮层血流动力学变化来间接反映神经活动。对于心理学、神经科学、康复医学等领域的研究者和学生而言掌握从原始数据到可发表图表的完整分析流程是开展研究的关键一步。然而这个过程涉及信号预处理、质量评估、统计建模和可视化等多个环节新手极易在软件选择、参数设置和结果解读上陷入误区。本文将围绕近红外数据分析的核心流程系统性地拆解每一步的操作与原理。无论你是刚开始接触 fNIRS还是已经积累了一些经验但在绘图和统计上遇到瓶颈都可以通过本文构建一个清晰、可复现的分析框架。我们将从理解 fNIRS 数据的基本结构开始逐步完成数据导入、预处理、绘制脑区激活图、以及进行组间统计比较的完整实操并重点解释每个步骤背后的“为什么”帮助你避免常见的分析陷阱。1. 理解 fNIRS 数据从光信号到血红蛋白浓度在动手处理数据之前必须清楚你手中的数据文件究竟代表了什么。fNIRS 设备发射特定波长的近红外光穿透头皮和颅骨部分光被脑组织吸收部分光被散射回来被探测器接收。设备记录的是光强随时间变化的原始信号。1.1 数据文件的核心结构一个典型的 fNIRS 数据文件如.nirs、.snirf或.mat格式通常包含以下几个关键矩阵d: 原始光密度数据。这是一个m×n的矩阵其中m是时间点数n是测量通道数。每个通道对应一个特定的“光源-探测器”对。s: 刺激标记。这是一个m×p的矩阵其中p是实验条件或任务的数量。在任务开始时对应列会在相应时间点标记为 1其余为 0。t: 时间轴。一个m×1的向量表示每个数据点对应的时间通常以秒为单位。aux: 辅助信号。如加速度计、心率等同步记录的生理信号。SD (Source-Detector) 结构体: 这是数据的“地图”它定义了SrcPos: 光源的三维坐标。DetPos: 探测器的三维坐标。MeasList: 一个列表定义了每个通道由哪个光源和哪个探测器组成以及使用的光的波长。% 示例在 MATLAB 中加载一个 .nirs 文件并查看其结构 load(subj01_task.nirs); whos % 你会看到类似以下的变量 % Name Size Bytes Class % d 3000x48 1152000 double % s 3000x2 48000 double % t 3000x1 24000 double % SD 1x1 ... struct1.2 血红蛋白浓度变化的计算原理fNIRS 直接测量的是光密度变化ΔOD。我们需要利用修正的朗伯-比尔定律将其转换为有生理意义的指标——氧合血红蛋白HbO和脱氧血红蛋白HbR的浓度变化。其核心公式基于双波长通常为 730nm 和 850nm 附近测量Δ[HbO] (ε_HbR_λ2 * ΔOD_λ1 - ε_HbR_λ1 * ΔOD_λ2) / (ε_HbO_λ1 * ε_HbR_λ2 - ε_HbO_λ2 * ε_HbR_λ1) * DPF * L Δ[HbR] (ε_HbO_λ1 * ΔOD_λ2 - ε_HbO_λ2 * ΔOD_λ1) / (ε_HbO_λ1 * ε_HbR_λ2 - ε_HbO_λ2 * ε_HbR_λ1) * DPF * L其中ε是血红蛋白在特定波长λ下的消光系数这是一个已知的物理常数。DPF是差分路径因子与光在组织中的散射路径有关通常与年龄和波长相关。L是光源与探测器之间的欧氏距离。注意实际操作中我们几乎从不手动计算这个公式。分析软件如 Homer2, NIRS-SPM, MNE-NIRS内置了此转换函数。理解原理的意义在于当转换结果异常时你可以排查是否是波长设置错误、DPF值不合理或SD距离信息有误。2. 搭建分析环境与数据准备一个稳定、可复现的分析环境至关重要。我们推荐使用MATLAB平台配合Homer2工具包因为它在学术界使用广泛教程丰富且预处理流程成熟。2.1 环境配置与工具包安装安装 MATLAB确保你拥有 MATLAB R2014b 或更高版本。Homer2 对版本有一定要求太旧的版本可能不兼容。下载 Homer2从官方 GitHub 仓库或指定网站下载 Homer2 工具包。添加到 MATLAB 路径将 Homer2 文件夹及其所有子文件夹添加到 MATLAB 的搜索路径中。这是最关键的一步路径添加不正确会导致函数无法调用。% 在 MATLAB 命令窗口中操作 % 假设你的 Homer2 解压在了 D:\Research\Tools\Homer2 addpath(genpath(D:\Research\Tools\Homer2)); savepath; % 保存路径设置下次启动 MATLAB 无需重新添加 % 验证安装运行以下命令应能看到 Homer2 的图形用户界面 homer22.2 数据组织规范在分析开始前建立清晰的项目文件夹结构能极大提升效率避免文件混乱。你的项目根目录/ ├── RawData/ # 存放原始设备导出的数据文件 │ ├── sub-01/ │ │ ├── sub-01_task-rest.nirs │ │ └── sub-01_task-motor.nirs │ ├── sub-02/ │ └── ... ├── DerivedData/ # 存放处理过程中产生的各级数据 │ ├── 01_Preprocessed/ │ ├── 02_HRF_Estimated/ │ └── ... ├── Code/ # 存放分析脚本 │ ├── preprocess_pipeline.m │ └── group_analysis.m ├── Figures/ # 存放生成的图表 │ ├── individual/ │ └── group/ └── README.md # 项目说明文档3. 核心分析流程实操从预处理到统计我们将以一个简单的“手指敲击任务 vs. 静息态”的组块设计实验为例展示完整流程。假设你已经有了sub-01_task-motor.nirs文件。3.1 数据预处理流程预处理的目标是去除噪声提取与神经活动相关的血流信号。标准流程通常按以下顺序进行检测并剔除不良通道信号质量极差如完全脱落的通道应先被标记。将原始光强转换为光密度hmrIntensity2OD运动伪迹校正这是最关键的一步。常用方法有样条插值法hmrMotionCorrectSpline。适用于尖峰状运动伪迹。PCA/ICA法hmrMotionCorrectPCA。适用于周期性或复杂运动。小波滤波法hmrMotionCorrectWavelet。带通滤波去除高频生理噪声如心率 ~1Hz和低频漂移如 Mayer 波 ~0.1Hz。典型滤波带宽为 0.01 - 0.5 Hz。hmrBandpassFilt将光密度转换为血红蛋白浓度hmrOD2Conc% preprocess_pipeline.m 脚本示例片段 function dataProc preprocess_pipeline(rawData, t, SD, s) % 1. 转换光密度 dOD hmrIntensity2OD(rawData); % 2. 运动伪迹校正 (以样条插值为例) % tMotion: 运动检测时间窗口 (秒) % tMask: 标记为运动伪迹后需要校正的时间范围 (秒) % STDEVthresh: 判断为运动伪迹的标准差阈值 % AMPthresh: 判断为运动伪迹的幅度阈值 [dOD_corrected, tIncChAuto] hmrMotionCorrectSpline(dOD, t, SD, tMotion0.5, tMask1.0, STDEVthresh50, AMPthresh5); % 3. 带通滤波 % hpf: 高通滤波截止频率 (Hz) % lpf: 低通滤波截止频率 (Hz) dOD_filtered hmrBandpassFilt(dOD_corrected, t, hpf0.01, lpf0.5); % 4. 转换为血红蛋白浓度 % ppf: 每波长对应的 DPF 因子数组如 [6.0, 6.0] [hb, ~] hmrOD2Conc(dOD_filtered, SD, ppf[6.0, 6.0]); % hb 是一个三维矩阵: [时间点 x 通道 x 血红蛋白种类] % 血红蛋白种类顺序通常是: HbO, HbR, HbT (总血红蛋白) dataProc.hb hb; dataProc.t t; dataProc.SD SD; dataProc.s s; dataProc.tIncChAuto tIncChAuto; % 保存自动标记的运动区间用于后续质量评估 end预处理参数选择背后的逻辑运动校正阈值STDEVthresh和AMPthresh设置过低会误将生理波动判为运动设置过高会漏检真实运动。建议先可视化原始信号观察典型运动伪迹的幅度再进行设置。滤波带宽0.01Hz 高通用于去除缓慢基线漂移0.5Hz 低通用于去除心搏~1.1Hz和呼吸~0.3Hz等高频生理噪声。如果你的任务频率很高可能需要提高低通截止频率。3.2 绘制单个被试的脑激活图获得干净的 HbO/HbR 时间序列后下一步是针对每个任务条件计算其引起的血红蛋白浓度平均变化并映射到大脑空间。构建一般线性模型对于组块设计我们可以使用 GLM 来估计每个通道、每种血红蛋白对任务条件的响应强度β值。% 使用 hmrDeconvHRF_DriftSS 函数进行 GLM 拟合 % t: 时间轴 % s: 刺激标记矩阵 % hb: 预处理后的血红蛋白数据 % SD: 源探测器信息 % trange: 分析的时间范围 [-2, 20] 表示任务前2秒到后20秒 % glmSolveMethod: 求解方法OLS普通最小二乘或AR-IRLS自回归迭代加权最小二乘抗噪声更强 [yavg, yavgstd, tHRF, ~, beta, ~] hmrDeconvHRF_DriftSS(dataProc.hb, dataProc.t, dataProc.s, dataProc.SD, trange[-2, 20], glmSolveMethodAR-IRLS);beta矩阵包含了我们关心的效应大小。它的维度是[通道数 x 条件数 x 血红蛋白种类]。将 β 值映射到脑表面我们需要一个标准大脑模板如 Colin27 或 MNI 脑和每个通道在模板上的坐标。首先根据SD.SrcPos和SD.DetPos利用空间配准方法如 3D 数字化仪记录或基于 EEG 10-20 系统的估计计算出每个通道的中点在大脑模板上的 MNI 坐标。然后使用插值方法如 Shepard‘s 方法将离散通道的 β 值平滑地映射到整个脑皮层表面。% 假设我们已经有了 channels.mni_pos (N通道 x 3坐标) 和 channels.beta_hbo (N通道 x 1 HbO的β值) % 使用 Homer2 的 hmrDisplayData 系列函数或自定义绘图 % 自定义简单映射示例需要 BrainNet Viewer 等工具辅助生成精美图片 % 1. 准备节点文件.node node_coords channels.mni_pos; node_values channels.beta_hbo; % 2. 准备表面文件.surf如 BrainMesh_ICBM152.nv % 3. 使用 BrainNet Viewer 加载并绘制更实用的方法对于初学者使用NIRS-SPM或MNE-NIRS进行统计和绘图可能更直观它们与标准脑模板的集成更好。例如在 NIRS-SPM 中你可以直接得到经过多重比较校正的统计参数图。3.3 组水平统计分析单个被试的结果受个体差异和噪声影响很大。我们需要在组水平上进行统计推断以确定哪些脑区的激活在群体中是稳定且显著的。准备二阶数据对每个被试提取其所有通道在目标条件下的 β 值例如任务条件 vs. 基线形成一个被试 × 通道的矩阵。选择统计检验单样本 t 检验检验组内每个通道的 β 值是否显著不等于 0即是否有激活。ttest(secondLevelData)配对样本 t 检验比较同一组被试在两个不同条件下如任务A vs. 任务B的激活差异。ttest(condA_data, condB_data)双样本 t 检验比较两组不同被试如患者组 vs. 对照组在同一条件下的激活差异。ttest2(group1_data, group2_data)多重比较校正由于我们同时检验了数十甚至上百个通道会大大增加犯 I 类错误假阳性的概率。必须进行校正。错误发现率FDR校正相对宽松控制的是在所有被拒绝的假设中错误拒绝的比例。族错误率FWE校正如 Bonferroni更为严格控制的是整个家族所有通道中至少出现一个假阳性的概率。聚类水平校正在 SPM 或 NIRS-SPM 中常用。它先定义空间上相邻且 t 值超过阈值的通道为一个“聚类”然后对聚类的规模如总通道数进行统计推断比基于单通道的校正更有力。% 组分析脚本示例片段 (使用 NIRS-SPM 风格) % 假设我们已经将所有被试的 beta 值整理成了一个结构体数组 S % S(i).beta (通道数 x 条件数) 代表第 i 个被试的数据 % 1. 指定设计矩阵例如单样本t检验 designMatrix ones(nSubjects, 1); % 所有被试属于同一组 contrast [1]; % 对比向量检验这组的均值是否0 % 2. 在通道水平进行统计NIRS-SPM 内部会处理 % 这通常通过调用 NIRS-SPM 的 nirs_spm_glm 等函数完成涉及指定模型、估计、对比和推断。 % 此处为逻辑示意非可运行代码。 [statMap, correctedP] perform_group_glm(S, designMatrix, contrast, correction, FDR);4. 常见问题、排查路径与最佳实践4.1 预处理阶段常见问题问题现象可能原因检查与排查方法解决方案转换后的 HbO/HbR 信号幅度异常大100 μM或异常小。1. 光源-探测器距离SD.SrcPos, SD.DetPos单位错误应是厘米或毫米。2. DPF 值设置错误。3. 原始数据光强值异常。1. 检查SD.SrcPos和SD.DetPos的数值范围正常在几厘米到十厘米之间。2. 检查hmrOD2Conc函数中ppf参数DPF值。3. 绘制原始光强信号d查看是否有通道饱和达到设备最大值或完全无信号。1. 确认坐标单位并与设备手册核对。2. 根据被试年龄和波长查阅文献设置正确的 DPF。3. 标记并排除信号饱和或无信号的通道。运动校正后信号出现严重畸变或“空洞”。1. 运动校正算法参数tMotion,tMask, 阈值过于激进。2. 数据中存在极长时程的运动干扰。1. 对比校正前后的信号图观察畸变发生的时段是否与tIncChAuto标记的运动区间吻合。2. 可视化原始光密度信号观察运动伪迹的真实形态。1. 调高STDEVthresh和AMPthresh阈值。2. 尝试不同的运动校正算法如 PCA 法。3. 考虑手动标记运动区间并使用hmrMotionArtifact进行分段剔除。滤波后任务相关的信号似乎也被削弱了。高通滤波截止频率 (hpf) 设置过高滤除了与任务相关的低频血流响应。检查任务设计的组块长度。血流动力学响应函数HRF是低频信号主峰约在5-6秒。如果组块长度较长如30秒0.01Hz的高通是安全的。确保高通滤波频率 (hpf) 低于任务频率的倒数。例如对于20秒的组块频率为0.05Hz高通应低于0.05Hz如设为0.01Hz。4.2 绘图与统计阶段常见问题问题现象可能原因检查与排查方法解决方案激活图显示在大脑非采样区域如小脑、深部核团。通道坐标MNI坐标配准错误。检查channels.mni_pos的坐标值。正常大脑皮层坐标的 Z 值矢状面通常为正。重新检查从 SD 结构到 MNI 空间的配准流程确保使用的模板如 Colin27和配准方法正确。组统计结果没有任何通道显著但单个被试图看起来有激活。1. 被试间激活模式变异太大。2. 多重比较校正过于严格如用了 Bonferroni。3. 预处理不一致引入了噪声。1. 查看所有被试单个通道的 β 值分布箱线图看是否方向不一致。2. 尝试不使用校正观察未校正的 p 值图。3. 检查每个被试的预处理日志确保参数和步骤一致。1. 考虑使用更宽松的 FDR 校正或聚类校正。2. 检查实验范式和预处理流程确保真正提取到了与任务锁定的信号。3. 增加样本量。绘制的脑图颜色映射混乱或无法覆盖预期区域。1. 插值方法或参数不当。2. 用于绘图的表面网格文件与坐标系统不匹配。1. 检查插值函数如griddata的输入参数通道坐标、值、网格点。2. 确认表面网格文件.nv, .surf的坐标范围与你的通道 MNI 坐标是否在同一空间。1. 使用成熟的绘图工具包如 NIRS-SPM 的nirsview或 MNE-Python 的plot_3d。2. 确保所有空间数据都转换到了同一标准空间如 MNI152。4.3 分析流程最佳实践清单原始数据备份永远保留一份原始的、未修改的数据文件。所有处理都在副本上进行。脚本化与记录使用 MATLAB 脚本或 Python 脚本记录每一步处理而不是在 GUI 中点击。保存关键的中间参数如运动校正阈值、滤波频率。质量评估先行在正式分析前花时间评估每个通道、每个被试的数据质量。计算并记录信噪比、运动伪迹比例、通道排除比例等指标。预处理流程一致对所有被试使用完全相同的预处理步骤和参数。如果某个被试需要特殊处理如更强的滤波需在论文方法部分明确说明理由。理解统计假设清楚你使用的统计检验如 t 检验的前提条件如正态性、方差齐性并在可能的情况下进行检验。校正多重比较只要进行空间上的多点比较就必须报告使用了何种多重比较校正方法及其阈值。结果可视化与解释并重漂亮的脑激活图是结果的一部分但必须在文中用文字准确描述激活的脑区提供解剖学名称或 Brodmann 分区、效应方向HbO 增加/HbR 减少以及统计效力t 值校正后 p 值。数据公开与代码共享如果条件允许将预处理后的数据和关键分析代码在公开仓库如 OSF, GitHub分享这有助于研究的可重复性。5. 扩展方向与工具选择掌握上述基础流程后你可以根据研究需求向更深入的方向探索功能连接分析计算不同脑区时间序列之间的相关性如皮尔逊相关、波谱相干性研究脑网络。工具推荐Homer2的hmrConnectivity或NIRS-KIT。脑区解码与机器学习使用模式分类如 SVM从 fNIRS 信号中解码被试的状态或意图。工具推荐MVPA-Light(MATLAB) 或scikit-learn(Python)。超扫描与人际神经科学分析两个或多个人在交互时脑活动之间的同步性。这需要专门的实验设计和分析工具如Hyperscanning分析流程。从 MATLAB 转向 Python如果你追求更免费、更现代的分析生态MNE-Python是当前增长最快的 fNIRS 分析工具它提供了从预处理、GLM 到群体统计和三维绘图的完整流水线且能与 EEG 分析无缝整合。无论选择哪条路径核心依然是扎实理解数据本质、审慎对待每个处理步骤、合理解读统计结果。从一份可靠的原始数据开始通过严谨、透明、可重复的分析流程最终得到经得起推敲的科学发现这才是数据分析工作的真正价值所在。 fNIRS数据分析全流程：从预处理到统计绘图与常见问题解决
近红外脑功能成像fNIRS是一种非侵入式、可穿戴的脑功能活动监测技术它通过测量大脑皮层血流动力学变化来间接反映神经活动。对于心理学、神经科学、康复医学等领域的研究者和学生而言掌握从原始数据到可发表图表的完整分析流程是开展研究的关键一步。然而这个过程涉及信号预处理、质量评估、统计建模和可视化等多个环节新手极易在软件选择、参数设置和结果解读上陷入误区。本文将围绕近红外数据分析的核心流程系统性地拆解每一步的操作与原理。无论你是刚开始接触 fNIRS还是已经积累了一些经验但在绘图和统计上遇到瓶颈都可以通过本文构建一个清晰、可复现的分析框架。我们将从理解 fNIRS 数据的基本结构开始逐步完成数据导入、预处理、绘制脑区激活图、以及进行组间统计比较的完整实操并重点解释每个步骤背后的“为什么”帮助你避免常见的分析陷阱。1. 理解 fNIRS 数据从光信号到血红蛋白浓度在动手处理数据之前必须清楚你手中的数据文件究竟代表了什么。fNIRS 设备发射特定波长的近红外光穿透头皮和颅骨部分光被脑组织吸收部分光被散射回来被探测器接收。设备记录的是光强随时间变化的原始信号。1.1 数据文件的核心结构一个典型的 fNIRS 数据文件如.nirs、.snirf或.mat格式通常包含以下几个关键矩阵d: 原始光密度数据。这是一个m×n的矩阵其中m是时间点数n是测量通道数。每个通道对应一个特定的“光源-探测器”对。s: 刺激标记。这是一个m×p的矩阵其中p是实验条件或任务的数量。在任务开始时对应列会在相应时间点标记为 1其余为 0。t: 时间轴。一个m×1的向量表示每个数据点对应的时间通常以秒为单位。aux: 辅助信号。如加速度计、心率等同步记录的生理信号。SD (Source-Detector) 结构体: 这是数据的“地图”它定义了SrcPos: 光源的三维坐标。DetPos: 探测器的三维坐标。MeasList: 一个列表定义了每个通道由哪个光源和哪个探测器组成以及使用的光的波长。% 示例在 MATLAB 中加载一个 .nirs 文件并查看其结构 load(subj01_task.nirs); whos % 你会看到类似以下的变量 % Name Size Bytes Class % d 3000x48 1152000 double % s 3000x2 48000 double % t 3000x1 24000 double % SD 1x1 ... struct1.2 血红蛋白浓度变化的计算原理fNIRS 直接测量的是光密度变化ΔOD。我们需要利用修正的朗伯-比尔定律将其转换为有生理意义的指标——氧合血红蛋白HbO和脱氧血红蛋白HbR的浓度变化。其核心公式基于双波长通常为 730nm 和 850nm 附近测量Δ[HbO] (ε_HbR_λ2 * ΔOD_λ1 - ε_HbR_λ1 * ΔOD_λ2) / (ε_HbO_λ1 * ε_HbR_λ2 - ε_HbO_λ2 * ε_HbR_λ1) * DPF * L Δ[HbR] (ε_HbO_λ1 * ΔOD_λ2 - ε_HbO_λ2 * ΔOD_λ1) / (ε_HbO_λ1 * ε_HbR_λ2 - ε_HbO_λ2 * ε_HbR_λ1) * DPF * L其中ε是血红蛋白在特定波长λ下的消光系数这是一个已知的物理常数。DPF是差分路径因子与光在组织中的散射路径有关通常与年龄和波长相关。L是光源与探测器之间的欧氏距离。注意实际操作中我们几乎从不手动计算这个公式。分析软件如 Homer2, NIRS-SPM, MNE-NIRS内置了此转换函数。理解原理的意义在于当转换结果异常时你可以排查是否是波长设置错误、DPF值不合理或SD距离信息有误。2. 搭建分析环境与数据准备一个稳定、可复现的分析环境至关重要。我们推荐使用MATLAB平台配合Homer2工具包因为它在学术界使用广泛教程丰富且预处理流程成熟。2.1 环境配置与工具包安装安装 MATLAB确保你拥有 MATLAB R2014b 或更高版本。Homer2 对版本有一定要求太旧的版本可能不兼容。下载 Homer2从官方 GitHub 仓库或指定网站下载 Homer2 工具包。添加到 MATLAB 路径将 Homer2 文件夹及其所有子文件夹添加到 MATLAB 的搜索路径中。这是最关键的一步路径添加不正确会导致函数无法调用。% 在 MATLAB 命令窗口中操作 % 假设你的 Homer2 解压在了 D:\Research\Tools\Homer2 addpath(genpath(D:\Research\Tools\Homer2)); savepath; % 保存路径设置下次启动 MATLAB 无需重新添加 % 验证安装运行以下命令应能看到 Homer2 的图形用户界面 homer22.2 数据组织规范在分析开始前建立清晰的项目文件夹结构能极大提升效率避免文件混乱。你的项目根目录/ ├── RawData/ # 存放原始设备导出的数据文件 │ ├── sub-01/ │ │ ├── sub-01_task-rest.nirs │ │ └── sub-01_task-motor.nirs │ ├── sub-02/ │ └── ... ├── DerivedData/ # 存放处理过程中产生的各级数据 │ ├── 01_Preprocessed/ │ ├── 02_HRF_Estimated/ │ └── ... ├── Code/ # 存放分析脚本 │ ├── preprocess_pipeline.m │ └── group_analysis.m ├── Figures/ # 存放生成的图表 │ ├── individual/ │ └── group/ └── README.md # 项目说明文档3. 核心分析流程实操从预处理到统计我们将以一个简单的“手指敲击任务 vs. 静息态”的组块设计实验为例展示完整流程。假设你已经有了sub-01_task-motor.nirs文件。3.1 数据预处理流程预处理的目标是去除噪声提取与神经活动相关的血流信号。标准流程通常按以下顺序进行检测并剔除不良通道信号质量极差如完全脱落的通道应先被标记。将原始光强转换为光密度hmrIntensity2OD运动伪迹校正这是最关键的一步。常用方法有样条插值法hmrMotionCorrectSpline。适用于尖峰状运动伪迹。PCA/ICA法hmrMotionCorrectPCA。适用于周期性或复杂运动。小波滤波法hmrMotionCorrectWavelet。带通滤波去除高频生理噪声如心率 ~1Hz和低频漂移如 Mayer 波 ~0.1Hz。典型滤波带宽为 0.01 - 0.5 Hz。hmrBandpassFilt将光密度转换为血红蛋白浓度hmrOD2Conc% preprocess_pipeline.m 脚本示例片段 function dataProc preprocess_pipeline(rawData, t, SD, s) % 1. 转换光密度 dOD hmrIntensity2OD(rawData); % 2. 运动伪迹校正 (以样条插值为例) % tMotion: 运动检测时间窗口 (秒) % tMask: 标记为运动伪迹后需要校正的时间范围 (秒) % STDEVthresh: 判断为运动伪迹的标准差阈值 % AMPthresh: 判断为运动伪迹的幅度阈值 [dOD_corrected, tIncChAuto] hmrMotionCorrectSpline(dOD, t, SD, tMotion0.5, tMask1.0, STDEVthresh50, AMPthresh5); % 3. 带通滤波 % hpf: 高通滤波截止频率 (Hz) % lpf: 低通滤波截止频率 (Hz) dOD_filtered hmrBandpassFilt(dOD_corrected, t, hpf0.01, lpf0.5); % 4. 转换为血红蛋白浓度 % ppf: 每波长对应的 DPF 因子数组如 [6.0, 6.0] [hb, ~] hmrOD2Conc(dOD_filtered, SD, ppf[6.0, 6.0]); % hb 是一个三维矩阵: [时间点 x 通道 x 血红蛋白种类] % 血红蛋白种类顺序通常是: HbO, HbR, HbT (总血红蛋白) dataProc.hb hb; dataProc.t t; dataProc.SD SD; dataProc.s s; dataProc.tIncChAuto tIncChAuto; % 保存自动标记的运动区间用于后续质量评估 end预处理参数选择背后的逻辑运动校正阈值STDEVthresh和AMPthresh设置过低会误将生理波动判为运动设置过高会漏检真实运动。建议先可视化原始信号观察典型运动伪迹的幅度再进行设置。滤波带宽0.01Hz 高通用于去除缓慢基线漂移0.5Hz 低通用于去除心搏~1.1Hz和呼吸~0.3Hz等高频生理噪声。如果你的任务频率很高可能需要提高低通截止频率。3.2 绘制单个被试的脑激活图获得干净的 HbO/HbR 时间序列后下一步是针对每个任务条件计算其引起的血红蛋白浓度平均变化并映射到大脑空间。构建一般线性模型对于组块设计我们可以使用 GLM 来估计每个通道、每种血红蛋白对任务条件的响应强度β值。% 使用 hmrDeconvHRF_DriftSS 函数进行 GLM 拟合 % t: 时间轴 % s: 刺激标记矩阵 % hb: 预处理后的血红蛋白数据 % SD: 源探测器信息 % trange: 分析的时间范围 [-2, 20] 表示任务前2秒到后20秒 % glmSolveMethod: 求解方法OLS普通最小二乘或AR-IRLS自回归迭代加权最小二乘抗噪声更强 [yavg, yavgstd, tHRF, ~, beta, ~] hmrDeconvHRF_DriftSS(dataProc.hb, dataProc.t, dataProc.s, dataProc.SD, trange[-2, 20], glmSolveMethodAR-IRLS);beta矩阵包含了我们关心的效应大小。它的维度是[通道数 x 条件数 x 血红蛋白种类]。将 β 值映射到脑表面我们需要一个标准大脑模板如 Colin27 或 MNI 脑和每个通道在模板上的坐标。首先根据SD.SrcPos和SD.DetPos利用空间配准方法如 3D 数字化仪记录或基于 EEG 10-20 系统的估计计算出每个通道的中点在大脑模板上的 MNI 坐标。然后使用插值方法如 Shepard‘s 方法将离散通道的 β 值平滑地映射到整个脑皮层表面。% 假设我们已经有了 channels.mni_pos (N通道 x 3坐标) 和 channels.beta_hbo (N通道 x 1 HbO的β值) % 使用 Homer2 的 hmrDisplayData 系列函数或自定义绘图 % 自定义简单映射示例需要 BrainNet Viewer 等工具辅助生成精美图片 % 1. 准备节点文件.node node_coords channels.mni_pos; node_values channels.beta_hbo; % 2. 准备表面文件.surf如 BrainMesh_ICBM152.nv % 3. 使用 BrainNet Viewer 加载并绘制更实用的方法对于初学者使用NIRS-SPM或MNE-NIRS进行统计和绘图可能更直观它们与标准脑模板的集成更好。例如在 NIRS-SPM 中你可以直接得到经过多重比较校正的统计参数图。3.3 组水平统计分析单个被试的结果受个体差异和噪声影响很大。我们需要在组水平上进行统计推断以确定哪些脑区的激活在群体中是稳定且显著的。准备二阶数据对每个被试提取其所有通道在目标条件下的 β 值例如任务条件 vs. 基线形成一个被试 × 通道的矩阵。选择统计检验单样本 t 检验检验组内每个通道的 β 值是否显著不等于 0即是否有激活。ttest(secondLevelData)配对样本 t 检验比较同一组被试在两个不同条件下如任务A vs. 任务B的激活差异。ttest(condA_data, condB_data)双样本 t 检验比较两组不同被试如患者组 vs. 对照组在同一条件下的激活差异。ttest2(group1_data, group2_data)多重比较校正由于我们同时检验了数十甚至上百个通道会大大增加犯 I 类错误假阳性的概率。必须进行校正。错误发现率FDR校正相对宽松控制的是在所有被拒绝的假设中错误拒绝的比例。族错误率FWE校正如 Bonferroni更为严格控制的是整个家族所有通道中至少出现一个假阳性的概率。聚类水平校正在 SPM 或 NIRS-SPM 中常用。它先定义空间上相邻且 t 值超过阈值的通道为一个“聚类”然后对聚类的规模如总通道数进行统计推断比基于单通道的校正更有力。% 组分析脚本示例片段 (使用 NIRS-SPM 风格) % 假设我们已经将所有被试的 beta 值整理成了一个结构体数组 S % S(i).beta (通道数 x 条件数) 代表第 i 个被试的数据 % 1. 指定设计矩阵例如单样本t检验 designMatrix ones(nSubjects, 1); % 所有被试属于同一组 contrast [1]; % 对比向量检验这组的均值是否0 % 2. 在通道水平进行统计NIRS-SPM 内部会处理 % 这通常通过调用 NIRS-SPM 的 nirs_spm_glm 等函数完成涉及指定模型、估计、对比和推断。 % 此处为逻辑示意非可运行代码。 [statMap, correctedP] perform_group_glm(S, designMatrix, contrast, correction, FDR);4. 常见问题、排查路径与最佳实践4.1 预处理阶段常见问题问题现象可能原因检查与排查方法解决方案转换后的 HbO/HbR 信号幅度异常大100 μM或异常小。1. 光源-探测器距离SD.SrcPos, SD.DetPos单位错误应是厘米或毫米。2. DPF 值设置错误。3. 原始数据光强值异常。1. 检查SD.SrcPos和SD.DetPos的数值范围正常在几厘米到十厘米之间。2. 检查hmrOD2Conc函数中ppf参数DPF值。3. 绘制原始光强信号d查看是否有通道饱和达到设备最大值或完全无信号。1. 确认坐标单位并与设备手册核对。2. 根据被试年龄和波长查阅文献设置正确的 DPF。3. 标记并排除信号饱和或无信号的通道。运动校正后信号出现严重畸变或“空洞”。1. 运动校正算法参数tMotion,tMask, 阈值过于激进。2. 数据中存在极长时程的运动干扰。1. 对比校正前后的信号图观察畸变发生的时段是否与tIncChAuto标记的运动区间吻合。2. 可视化原始光密度信号观察运动伪迹的真实形态。1. 调高STDEVthresh和AMPthresh阈值。2. 尝试不同的运动校正算法如 PCA 法。3. 考虑手动标记运动区间并使用hmrMotionArtifact进行分段剔除。滤波后任务相关的信号似乎也被削弱了。高通滤波截止频率 (hpf) 设置过高滤除了与任务相关的低频血流响应。检查任务设计的组块长度。血流动力学响应函数HRF是低频信号主峰约在5-6秒。如果组块长度较长如30秒0.01Hz的高通是安全的。确保高通滤波频率 (hpf) 低于任务频率的倒数。例如对于20秒的组块频率为0.05Hz高通应低于0.05Hz如设为0.01Hz。4.2 绘图与统计阶段常见问题问题现象可能原因检查与排查方法解决方案激活图显示在大脑非采样区域如小脑、深部核团。通道坐标MNI坐标配准错误。检查channels.mni_pos的坐标值。正常大脑皮层坐标的 Z 值矢状面通常为正。重新检查从 SD 结构到 MNI 空间的配准流程确保使用的模板如 Colin27和配准方法正确。组统计结果没有任何通道显著但单个被试图看起来有激活。1. 被试间激活模式变异太大。2. 多重比较校正过于严格如用了 Bonferroni。3. 预处理不一致引入了噪声。1. 查看所有被试单个通道的 β 值分布箱线图看是否方向不一致。2. 尝试不使用校正观察未校正的 p 值图。3. 检查每个被试的预处理日志确保参数和步骤一致。1. 考虑使用更宽松的 FDR 校正或聚类校正。2. 检查实验范式和预处理流程确保真正提取到了与任务锁定的信号。3. 增加样本量。绘制的脑图颜色映射混乱或无法覆盖预期区域。1. 插值方法或参数不当。2. 用于绘图的表面网格文件与坐标系统不匹配。1. 检查插值函数如griddata的输入参数通道坐标、值、网格点。2. 确认表面网格文件.nv, .surf的坐标范围与你的通道 MNI 坐标是否在同一空间。1. 使用成熟的绘图工具包如 NIRS-SPM 的nirsview或 MNE-Python 的plot_3d。2. 确保所有空间数据都转换到了同一标准空间如 MNI152。4.3 分析流程最佳实践清单原始数据备份永远保留一份原始的、未修改的数据文件。所有处理都在副本上进行。脚本化与记录使用 MATLAB 脚本或 Python 脚本记录每一步处理而不是在 GUI 中点击。保存关键的中间参数如运动校正阈值、滤波频率。质量评估先行在正式分析前花时间评估每个通道、每个被试的数据质量。计算并记录信噪比、运动伪迹比例、通道排除比例等指标。预处理流程一致对所有被试使用完全相同的预处理步骤和参数。如果某个被试需要特殊处理如更强的滤波需在论文方法部分明确说明理由。理解统计假设清楚你使用的统计检验如 t 检验的前提条件如正态性、方差齐性并在可能的情况下进行检验。校正多重比较只要进行空间上的多点比较就必须报告使用了何种多重比较校正方法及其阈值。结果可视化与解释并重漂亮的脑激活图是结果的一部分但必须在文中用文字准确描述激活的脑区提供解剖学名称或 Brodmann 分区、效应方向HbO 增加/HbR 减少以及统计效力t 值校正后 p 值。数据公开与代码共享如果条件允许将预处理后的数据和关键分析代码在公开仓库如 OSF, GitHub分享这有助于研究的可重复性。5. 扩展方向与工具选择掌握上述基础流程后你可以根据研究需求向更深入的方向探索功能连接分析计算不同脑区时间序列之间的相关性如皮尔逊相关、波谱相干性研究脑网络。工具推荐Homer2的hmrConnectivity或NIRS-KIT。脑区解码与机器学习使用模式分类如 SVM从 fNIRS 信号中解码被试的状态或意图。工具推荐MVPA-Light(MATLAB) 或scikit-learn(Python)。超扫描与人际神经科学分析两个或多个人在交互时脑活动之间的同步性。这需要专门的实验设计和分析工具如Hyperscanning分析流程。从 MATLAB 转向 Python如果你追求更免费、更现代的分析生态MNE-Python是当前增长最快的 fNIRS 分析工具它提供了从预处理、GLM 到群体统计和三维绘图的完整流水线且能与 EEG 分析无缝整合。无论选择哪条路径核心依然是扎实理解数据本质、审慎对待每个处理步骤、合理解读统计结果。从一份可靠的原始数据开始通过严谨、透明、可重复的分析流程最终得到经得起推敲的科学发现这才是数据分析工作的真正价值所在。