CN104644173A - 一种抑郁风险三级预警方法及系统 - Google Patents

Abstract

本发明涉及一种抑郁风险三级预警方法及系统，所述方法包括：采集静息态和任务态的BOLD信号及其对应的时间序列；从BOLD信号中提取静息态BOLD值和任务态BOLD值；根据静息态BOLD值和任务态BOLD值生成静息态BOLD变化率曲线和任务态BOLD变化率曲线；对静息态BOLD变化率曲线和任务态BOLD变化率曲线进行向量分析以得到变化强度值||ΔG||和变化矢量角θ；根据变化强度值||ΔG||和变化矢量角θ与预设的变化强度阈值||ΔG||^阈值和变化矢量角阈值θ^阈值进行预警级别判定。本发明通过对人体与抑郁症有关的生理指标进行定量的分析，实现了对不同程度抑郁高危人群的预警，并且预警稳定性好、正确率高。

Description

一种抑郁风险三级预警方法及系统

技术领域

本发明属于生物信息、疾病预警技术领域，特别涉及一种抑郁风险三级预警方法及系统。

背景技术

随着现代城市生活节奏不断地加快，人们压力越来越大，抑郁症的发病率也在逐年递增。目前抑郁症已成为一种严重危害人类身心健康的精神类疾病。

抑郁高危是介于正常的抑郁情绪与需要治疗的抑郁症之间的状态，这种身体和心理上的不适应对相应群体的生活同样造成影响，更重要的是，抑郁高危的轻、中、重是一个恶性变化的过程。抑郁症高危者处于较恶劣的心理环境中，这种状态会使人更容易进入更为抑郁的心理环境中，所以抑郁高危者极为容易进一步发展为抑郁症。因此，对抑郁高危人群进行风险预警刻不容缓。

近年来，功能影像学技术迅速发展,为精神疾病的研究提供了新思路。与其它功能影像学技术相比，功能性磁共振成像技术(fMRI，functional magnetic resonance imaging)能直接显示激活脑区的部位及激活程度，实现了功能与结构的充分融合，图像的空间、时间分辨率高、重复性好，对人体无离子辐射创伤性。fMRl技术是以脱氧血红蛋白的顺磁性敏感效应为基础的成像技术。人体血液中的血红蛋白是抗磁性物质,脱氧血红蛋白是顺磁性物质。人脑在接受各种感觉刺激、执行运动、进行思维、记忆及情绪体验等精神活动时，脑部特定区域会出现血流动力学变化，如在接受视觉、听觉或触觉等外界特定刺激后其相应的皮质功能区被激活,出现局部脑血流量(region cerebral blood flood，rCBF)增加。这种rCBF增加造成局部氧合血红蛋白增加，但局部氧耗量增加不明显，即脱氧血红蛋白相对减少，使局部磁化率发生改变，在MRI的T₂加权出现局部信号增加，能够很好的反应脑区激活部位和激活程度。

利用fMRI技术能够获得抑郁高危预警时所需的客观生理指标，为生理指标的定量分析提供了可能，从而为抑郁高危人群的预警提供了新的方式，可进一步实现对抑郁症的有效防范。

发明内容

本发明的目的是提供一种抑郁风险三级预警方法及系统，通过对脑区生理指标进行定量的分析，实现对不同程度抑郁高危人群的预警。

为解决上述问题，本发明的一个方面提供了一种抑郁风险三级预警方法，包括：步骤S1，采集静息态和任务态的BOLD信号及其对应的时间序列；步骤S2，从BOLD信号中提取静息态BOLD值和任务态BOLD值；步骤S3，根据静息态BOLD值和任务态BOLD值生成静息态BOLD变化率曲线和任务态BOLD变化率曲线；步骤S4，对静息态BOLD变化率曲线和任务态BOLD变化率曲线进行向量分析以得到变化强度值||ΔG||和变化矢量角θ；步骤S5，根据变化强度值||ΔG||和变化矢量角θ与预设的变化强度阈值||ΔG||^阈值和变化矢量角阈值θ^阈值进行预警级别判定。

本发明的又一个方面提供了一种抑郁风险三级预警系统，包括BOLD信号采集装置、BOLD值提取装置、BOLD变化率曲线生成装置、变化强度值和变化矢量角生成装置和预警级别判断装置；所述BOLD信号采集装置，用于采集静息态和任务态的BOLD信号及其对应的时间序列；所述BOLD值提取装置，用于从所述BOLD信号中提取静息态BOLD值和任务态BOLD值；BOLD变化率曲线生成装置，用于根据静息态BOLD值和任务态BOLD值生成静息态BOLD变化率曲线和任务态BOLD变化率曲线；所述变化强度值和变化矢量角生成装置，用于对静息态BOLD变化率曲线和任务态BOLD变化率曲线进行向量分析以得到变化强度值||ΔG||和变化矢量角θ；和所述预警级别判断装置，用于根据变化强度值||ΔG||和变化矢量角θ与预设的变化强度阈值||ΔG||^阈值和变化矢量角阈值θ^阈值进行预警级别判定。

有益效果：本发明通过fMRI技术，能够获得人体与抑郁症有关的生理指标，这为生理指标的定量分析提供了可能，从而为抑郁高危人群的预警提供了新的方式，可进一步实现对抑郁症的有效防范，并且预警稳定性好、正确率高。

附图说明

图1是根据本发明第一实施方式的流程图；

图2是根据本发明第一实施方式中步骤S3的流程图；

图3是根据本发明第一实施方式中步骤S4的流程图；

图4是根据本发明第二实施方式的系统结构图；

图5是根据本发明第二实施方式的BOLD变化率曲线生成装置的结构图；

图6是根据本发明第二实施方式的变化强度值和变化矢量角生成装置的结构图；

图7是本发明一具体示例的抑郁高危蓝色预警阈值确定曲线图；

图8是本发明一具体示例的抑郁高危橙色预警阈值确定曲线图；

图9是本发明一具体示例的抑郁高危红色预警阈值确定曲线图；

图10是本发明一具体示例的抑郁高危极限预警阈值确定曲线图；

图11是本发明抑郁高危风险三级预警模型图。

具体实施方式

为使本发明的目的、技术方案和优点更加清楚明了，下面结合具体实施方式并参照附图，对本发明进一步详细说明。应该理解，这些描述只是示例性的，而并非要限制本发明的范围。此外，在以下说明中，省略了对公知结构和技术的描述，以避免不必要地混淆本发明的概念。

如前所示，脱氧血红蛋白(Deoxyhemoglobin)比氧合血红蛋白(Oxyhemoglobin)更具有顺磁性，其本身就有和组织一样的磁敏感性。因此可以采用脱氧血红蛋白作为对比剂。如果影响大脑的状态使氧摄取和血流之间产生不平衡，并采用对磁场不均匀性敏感的MRI成像序列，就可在脑皮层血管周围得到MRI信号的变化。此技术称作血氧合度依赖值(BloodOxygenation Level Dependent，BOLD)的对比。

本发明中的BOLD值表示血氧水平依赖信号值，其是通过对四个脑区的脑核磁共振进行扫描，获取脑核磁共振图像数据(即BOLD信号)，从而得到BOLD值原始数据。

图1是根据本发明第一实施方式的流程图。

步骤S1，采集静息态和任务态的BOLD信号及其对应的时间序列。

首先设计刺激任务，因抑郁高危者存在快感缺失症状，在奖励期待任务中，其对快乐的体验和期待能力下降，奖励期待任务的变化会引起相关脑区BOLD值的变化，并且不同抑郁高危程度的BOLD值变化存在差异。

基于上述理论，本步骤首先设计引起相关脑区BOLD值变化的刺激任务，其共分为3类任务：奖赏任务、惩罚任务和非奖罚任务，其中奖赏任务共有3个水平(例如“赢得1元”、“赢得3元”和“赢得5元”)，惩罚任务共有3个水平(例如“损失1元”、“损失3元”和“损失5元”)，非奖罚任务有1个水平，其中每个水平任务由20个子任务构成，任务在fMRI扫描过程中进行。

其次，屏幕中央会出现刺激任务提示线索，持续时间为250ms，然后出现“+”号，持续时间2000-2500ms，接着呈现目标按键，按键后会经过1250-1850ms的心理期待，最后屏幕出现1650ms的反馈，根据不同的任务，会显示“您赢得了*元！共*元”或“您输了*元！共*元”或“0，共*元”。当前子任务完成后进入下一子任务。

首先，采集静息态下的四个脑区的BOLD信号)，其中，静息态扫描时长为720000ms，体素大小为3.75*3.75*4.0mm³(即体积大小)。四个脑区包括纹状体、岛叶、伏隔核和杏仁体。

其次，采集任务态下的四个脑区的BOLD信号，其中任务态扫描时长t为720000ms，体素大小为3.75*3.75*4.0mm³，四个脑区包括纹状体、岛叶、伏隔核和杏仁体。

在该实施方式中，可以采用fMRI扫描设备，记录静息态和任务态时四个脑区BOLD值原始数据。例如，采用西门子3.0(3.0T Siemens Allegrasystem)。

步骤S2，从BOLD信号中提取静息态BOLD值和任务态BOLD值。

对采集到的BOLD信号进行预处理包括下述步骤，首先进行数据转换，将原始Dicom图像转换为NIfTI图像，其中NIfTI是一种典型的数据分析格式；其次，经过时间校正和头动校正，校正后的图像配准到SPM软件中的标准蒙特利尔(Montreal Neurological Institute，MNI)空间(体素大小为3*3*3mm³)进行重新采样，进行信号的去除线性漂移和滤波，回归掉白质信号、脑脊液信号以及头动对真实生理信号的影响。

最后，从抑郁高危激活脑区，按照时间序列从BOLD信号中分别提取k个体素的静息态BOLD值和任务态BOLD值。

其中，静息态扫描时长t为720000ms，体素大小3.75*3.75*4.0mm³，任务态扫描时长t为720000ms，体素大小3.75*3.75*4.0mm³；i代表脑区，1为纹状体，2为岛叶，3为伏隔核，4为杏仁体。

步骤S3，根据静息态BOLD值和任务态BOLD值生成静息态BOLD变化率曲线和任务态BOLD变化率曲线。

本步骤中，根据静息态BOLD值和任务态BOLD值，基于时间序列计算四个脑区的BOLD值变化率，根据计算结果生成静息态BOLD变化率曲线和任务态BOLD变化率曲线。

步骤S4，对静息态BOLD变化率曲线和任务态BOLD变化率曲线进行向量分析以得到变化强度值||ΔG||和变化矢量角θ。

其中，所述抑郁高危风险预警阈值包括变化强度阈值||ΔG||^阈值和变化矢量角阈值θ^阈值，阈值的确定是将曲线间差异量化得到的，其计算步骤如下：

预先采集n名轻度抑郁高危、n名中度抑郁高危、n名重度抑郁高危和n名轻度抑郁症患者的抑郁高危激活脑区任务态和静息态BOLD信号，按照时间序列和脑区顺序分别提取k个体素的BOLD值，静息态时长t为720000ms，体素大小3.75*3.75*4.0mm³，任务态时长t为720000ms，体素大小为3.75*3.75*4.0mm³，所述激活脑区包括纹状体(脑区1)、岛叶(脑区2)、伏隔核(脑区3)和杏仁体(脑区4)。将其同一抑郁高危程度相关脑区静息态时同一时间序列t时的n个BOLD值的均值作为相应静息态时间序列t时的BOLD值，计算公式为：

${\mathrm{bold}}_i^{\mathrm{rest\; j}}\left(t\right)=\frac{\underset{1}{\overset{n}{\mathrm{\Σ}}}{\mathrm{bold}}_i^{\mathrm{rest\; J\; j}}\left(t\right)}{n},$

其中i＝1,2,3,4；J＝1,2,…,n；j＝1,2,…,k；t∈[0,720000]；

其中为静息态时间序列t时第J个人第i个脑区第j个体素的BOLD值，n为人数，为计算完成的静息态相应时间序列t时BOLD值。

计算相关脑区任务态相应时间序列t时的BOLD值，计算公式为：

${\mathrm{bold}}_i^{\mathrm{rest\; j}}\left(t\right)=\frac{\underset{1}{\overset{n}{\mathrm{\Σ}}}{\mathrm{bold}}_i^{\mathrm{rest\; J\; j}}\left(t\right)}{n},$

其中i＝1,2,3,4；J＝1,2,…,n；j＝1,2,…,k；t∈[0,720000]；

其中为任务态时间序列t时第J个人第i个脑区第j个体素的BOLD值，n为人数，为计算完成的相关脑区任务态相应时间序列t时BOLD值；

按照计算静息态和任务态BOLD变化率基线值的计算方法，计算出不同抑郁高危程度静息态BOLD基线值和任务态BOLD基线值进一步计算出静息态BOLD变化率和任务态BOLD变化率

按照脑区顺序和时间序列，生成轻度抑郁高危静息态BOLD变化率曲线和轻度抑郁高危任务态BOLD变化率曲线，中度抑郁高危静息态BOLD变化率曲线和中度抑郁高危任务态BOLD变化率曲线，重度抑郁高危静息态BOLD变化率曲线和重度抑郁高危任务态BOLD变化率曲线，轻度抑郁症静息态BOLD变化率曲线和轻度抑郁症任务态BOLD变化率曲线，其中横坐标为时间序列，纵坐标为BOLD变化率值。

步骤S5，根据变化强度值||ΔG||和变化矢量角θ与预设的变化强度阈值||ΔG||^阈值和变化矢量角阈值θ^阈值进行预警级别判定。

具体的，依据已有的抑郁高危风险三级预警阈值进行预警级别判定步骤如下：

抑郁高危风险蓝色预警判定：且满足该条件表示抑郁高危风险的预警级别为低。

抑郁高危风险橙色预警判定：且满足该条件表示抑郁高危风险的预警级别为中。

抑郁高危风险红色预警判定：且满足该条件表示抑郁高危风险的预警级别为重。

其预警级别为：低-中-高，对应抑郁高危程度为：轻-中-重。

图2是根据本发明第一实施方式中步骤S3的流程图。

如图2所示，在本发明的优选实施例中，前述步骤S3进一步包括下述步骤：

步骤S31,计算四个脑区静息态BOLD变化率基线值。

本步骤中，根据以下公式分别计算四个脑区静息态BOLD变化率基线值：

${\mathrm{bold}}_i^{\mathrm{rest}}=\frac{\underset{1}{\overset{k}{\mathrm{\Σ}}}{\mathrm{bold}}_i^{\mathrm{rest\; j}}\left(t\right)}{k},$

i＝1,2,3,4；j＝1,2,…,k；t∈[0,720000]；

其中，为静息态时间序列t时第i个脑区第j个体素的BOLD值，为计算完成的静息态时第i个脑区BOLD变化率基线值。

步骤S32,计算四个脑区任务态BOLD变化率基线值。

本步骤中，根据以下公式分别计算四个脑区任务态BOLD变化率基线值：

${\mathrm{bold}}_i^{\mathrm{rest}}=\frac{\underset{1}{\overset{k}{\mathrm{\Σ}}}{\mathrm{bold}}_i^{\mathrm{rest\; j}}\left(t\right)}{k},$

i＝1,2,3,4；j＝1,2,…,k；t∈[0,720000]；

其中，为任务态时间序列t时第i个脑区第j个体素的BOLD值，为计算完成的任务态时第i个脑区BOLD变化率基线值。

步骤S33,计算四个脑区静息态BOLD变化率。

本步骤中，根据以下公式分别计算四个脑区静息态BOLD变化率：

$R_i^{\mathrm{rest\; j}}\left(t\right)=\frac{{\mathrm{bold}}_i^{\mathrm{rest\; j}}\left(t\right)-{\mathrm{bold}}_i^{\mathrm{rdst}}}{{\mathrm{bold}}_i^{\mathrm{rest}}},$

i＝1,2,3,4；j＝1,2,…,k；t∈[0,720000]；

其中为静息态时间序列t时第i个脑区第j个体素的BOLD值，为静息态时第i个脑区BOLD变化率基线值，为计算完成的相应脑区静息态体素BOLD变化率。

步骤S34,计算四个脑区任务态BOLD变化率。

本步骤中，根据以下公式分别计算四个脑区任务态BOLD变化率：

$R_i^{\mathrm{rest\; j}}\left(t\right)=\frac{{\mathrm{bold}}_i^{\mathrm{rest\; j}}\left(t\right)-{\mathrm{bold}}_i^{\mathrm{rdst}}}{{\mathrm{bold}}_i^{\mathrm{rest}}},$

i＝1,2,3,4；j＝1,2,…,k；t∈[0,720000]；

其中为任务态时间序列t时第i个脑区第j个体素的BOLD值，为任务态时第i个脑区BOLD变化率基线值，为计算完成的相应脑区任务态体素BOLD变化率。

步骤S35,生成静息态BOLD变化率曲线和任务态BOLD变化率曲线。

本步骤中，按照脑区顺序和时间序列，生成静息态BOLD变化率曲线和任务态BOLD变化率曲线，横坐标为时间序列，纵坐标为BOLD变化率值。

具体的，在本步骤中分别计算出四个脑区静息态和任务态BOLD变化率基线值，进一步计算出四个脑区静息态和任务态BOLD变化率，根据脑区顺序和时间序列，生成静息态BOLD变化率曲线和任务态BOLD变化率曲线，其中曲线横坐标为时间序列，纵坐标为BOLD变化率值。

图3是根据本发明第一实施方式中步骤S4的流程图。

如图3所示，在本发明的一个优选实施例中，步骤S4包括下述步骤：

步骤S41,根据静息态BOLD变化率曲线和任务态BOLD变化率曲线得到BOLD变化率曲线图。

将轻度抑郁高危静息态BOLD变化率曲线和轻度抑郁高危任务态BOLD变化率曲线转化到同一坐标中，使其在X轴上具有相同的投影(时间序列相同)，得到相应的脑区BOLD变化率曲线Ⅰ，如图7所示。

步骤S42,对BOLD变化率曲线图进行向量分析以得到变化向量。

对BOLD变化率曲线图进行向量分析，定量刻画曲线间的差异，计算出变化强度值||ΔG||和变化矢量角θ。变化率曲线图中任务态BOLD变化率曲线每相邻的两个坐标构成一个向量，四个脑区共得到4*k个坐标，按照时间序列构成4*k-1个向量，即构成向量集静息态BOLD变化率曲线也构成4*k-1个向量，即构成向量集静息态BOLD变化率曲线到任务态BOLD变化率曲线的变化，可由向量的变化表示，计算公式如下：

H为任务态BOLD变化率曲线构成的向量集，G为静息态BOLD变化率曲线构成的向量集，ΔG表示变化向量。

例如，对于轻度抑郁高危静息态BOLD变化率曲线I来说，轻度抑郁高危静息态BOLD变化率曲线到轻度抑郁高危任务态BOLD变化率曲线的变化向量为：

其中，H_Ⅰ为轻度抑郁高危任务态BOLD变化率曲线构成的向量集，G_Ⅰ为轻度抑郁高危静息态BOLD变化率曲线构成的向量集，ΔG_Ⅰ表示变化向量。

步骤S43，根据变化向量计算变化强度值。

根据变化向量ΔG，计算变化向量ΔG的变化强度值||ΔG||。变化向量ΔG的变化强度值为所有变化向量的模||ΔG||的均值，由欧氏距离(Euclideandistance)决定，其计算公式如下：

其中||ΔG||为变化强度值。

对于曲线I来说，变化向量ΔG_Ⅰ的变化强度值为所有变化向量的模||ΔG_Ⅰ||的均值，计算公式如下：

其中为蓝色预警变化强度阈值。

按照曲线Ⅰ的处理方法，对中度抑郁高危任务态BOLD变化率曲线和中度抑郁高危静息态BOLD变化率曲线进行处理，获得BOLD变化率曲线Ⅱ，进而得到向量集则变化向量为

其变化强度为得到橙色预警变化强度阈值

根据曲线Ⅰ的处理方法，对重度抑郁高危任务态BOLD变化率曲线和重度抑郁高危静息态BOLD变化率曲线进行处理，获得BOLD变化率曲线Ⅲ，进而得到向量集和则变化向量

为

变化强度为得到红色预警变化强度阈值

同样，依据曲线Ⅰ的处理方法，将轻度抑郁症任务态BOLD变化率曲线和轻度抑郁症静息态BOLD变化率曲线进行处理，获得BOLD变化率曲线Ⅳ，进而得到向量集和则变化向量为：

变化强度为得到预警极限阈值

步骤S44，根据向量夹角集计算变化矢量角。

变化矢量角阈值θ^阈值的计算步骤如下：计算曲线Ⅰ中轻度抑郁高危任务态BOLD变化率曲线和轻度抑郁高危静息态BOLD变化率曲线在同一时间序列段m时，两向量的夹角，计算公式如下：

其中和为向量集G和H中同一时间序列段m时对应的向量，θ_m为对应的两向量夹角，得到向量夹角集θ＝(θ₁，θ₂，…，θ_m，…θ_4*k-1)。

计算变化矢量角θ，其计算公式如下：

$\mathrm{\θ}=\frac{\underset{1}{\overset{4*k-1}{\mathrm{\Σ}}}{\mathrm{\θ}}_m}{4*k-1},m=\mathrm{1,2},\·\·\·,4*k-1;$

其中θ_m为静息态BOLD变化率曲线和任务态BOLD变化率曲线同一时间序列段m内两向量的夹角值，θ为计算完成的变化矢量角。

例如，对于曲线Ⅰ计算θ_Ⅰ的阈值，其计算公式如下：

其中θ_m为轻度抑郁高危静息态BOLD变化率曲线和轻度抑郁高危任务态BOLD变化率曲线同一时间序列段m内两向量的夹角值，为蓝色预警变化矢量角阈值。

依据曲线Ⅰ向量角的计算方法，得到曲线Ⅱ的向量夹角集θ_Ⅱ＝(θ_Ⅱ1，θ_Ⅱ2，…，θ_Ⅱ4*k-1)，计算θ_Ⅱ的阈值得到橙色预警变化矢量角阈值

根据曲线Ⅰ向量角的计算方法，得到曲线Ⅲ的向量夹角集θ_Ⅲ＝(θ_Ⅲ1，θ_Ⅲ2，…，θ_Ⅲ4*k-1)，计算θ_Ⅲ的阈值得到红色预警变化矢量角阈值

同样，按照曲线Ⅰ向量角的处理方法，得到曲线Ⅳ的向量夹角集θ_Ⅳ＝(θ_Ⅳ1，θ_Ⅳ2，…，θ_Ⅳ4*k-1)，计算θ_Ⅳ的阈值得到预警极限矢量角阈值

图4是根据本发明第二实施方式的系统结构图。

如图4所示，本发明第二实施方式提供的抑郁风险三级预警系统，包括BOLD信号采集装置1、BOLD值提取装置2、BOLD变化率曲线生成装置3、变化强度值和变化矢量角生成装置4和预警级别判断装置5。

BOLD信号采集装置1用于采集静息态和任务态下四个脑区的BOLD信号及其对应的时间序列。

具体的，采集静息态四个脑区(纹状体、岛叶、伏隔核和杏仁体)体素扫描时的BOLD值 和任务态四个脑区(纹状体、岛叶、伏隔核和杏仁体)体素扫描时的BOLD值 及其对应的时间序列t₁，t₂，…，t_4*k。

在本发明的实施方式中，BOLD信号采集装置1可采用西门子3.0(3.0TSiemens Allegra system)fMRI扫描设备。

BOLD值提取装置2与所述BOLD信号采集装置1连接，用于对采集到的BOLD信号进行预处理，并从中提取四个脑区体素的BOLD值。

BOLD变化率曲线生成装置3与所述BOLD信号采集装置1和所述BOLD值提取装置2连接，用于根据时间序列计算四个脑区的BOLD值变化率，根据计算结果生成静息态BOLD变化率曲线和任务态BOLD变化率曲线。

变化强度值和变化矢量角生成装置4与所述BOLD变化率曲线生成装置3连接，用于对静息态BOLD变化率曲线和任务态BOLD变化率曲线进行向量分析，定量刻画曲线间的差异，计算出变化强度值||ΔG||和变化矢量角θ。

预警级别判断装置5，与所述变化强度值和变化矢量角生成装置4连接，用于根据计算出的变化强度值||ΔG||和变化矢量角θ与预设的抑郁高危风险预警模型中的变化强度阈值||ΔG||^阈值和变化矢量角阈值θ^阈值进行预警级别判定。

具体的，所述预警级别判断装置5包括蓝色预警模块51、橙色预警模块52和红色预警模块53，用于依据抑郁高危风险三级模型中已有的变化强度阈值||ΔG||^阈值和变化矢量角阈值θ^阈值进行预警级别的判断。

所述蓝色预警模块51，用于当且时，则判定抑郁高危风险预警级别为蓝色预警。

所述橙色预警模块52，用于当且时，则判定抑郁高危风险预警级别为橙色预警。

所述红色预警模块53，用于当且时，则判定抑郁高危风险预警级别为红色预警。

图5是根据本发明第二实施方式的BOLD变化率曲线生成装置3的结构图。

如图5所示，BOLD变化率曲线生成装置3包括BOLD值提取模块31、静息态BOLD变化率基线值计算模块32、任务态BOLD变化率基线值计算模块33、静息态BOLD变化率计算模块34、任务态BOLD变化率计算模块35和BOLD变化率曲线生成模块36。

BOLD值提取模块31用于从BOLD信号中提取静息态BOLD值和任务态BOLD值。具体的，BOLD值提取模块31根据时间序列，从四个脑区的脑核磁共振图像中分别提取k个体素的BOLD值，包括静息态时间序列t时第i个脑区第j个体素的BOLD值以及任务态时间序列t时第i个脑区第j个体素的BOLD值

在一个具体实施例中，静息态扫描时长t为720000ms，体素大小3.75*3.75*4.0mm³，任务态扫描时长t为720000ms，体素大小3.75*3.75*4.0mm³；i为脑区，1为纹状体，2为岛叶，3为伏隔核，4为杏仁体。

静息态BOLD变化率基线值计算模块32与所述BOLD值提取模块31连接，用于根据静息态BOLD值获取静息态BOLD变化率基线值。具体来说，静息态BOLD变化率基线值计算模块32根据以下公式：

${\mathrm{bold}}_i^{\mathrm{rest}}=\frac{\underset{1}{\overset{k}{\mathrm{\Σ}}}{\mathrm{bold}}_i^{\mathrm{rest\; j}}\left(t\right)}{k},i=\mathrm{1,2,3,4};j=\mathrm{1,2},\·\·\·,k;t\∈\left[\mathrm{0,720000}\right];$

分别计算四个脑区静息态BOLD变化率基线值，其中，为静息态时间序列t时第i个脑区第j个体素的BOLD值，为计算完成的静息态时第i个脑区BOLD变化率基线值。

任务态BOLD变化率基线值计算模块33与所述BOLD值提取模块31连接，用于根据任务态BOLD值获取任务态BOLD变化率基线值。

具体来说，任务态BOLD变化率基线值计算模块33，根据以下公式：

${\mathrm{bold}}_i^{\mathrm{rest}}=\frac{\underset{1}{\overset{k}{\mathrm{\Σ}}}{\mathrm{bold}}_i^{\mathrm{rest\; j}}\left(t\right)}{k},i=\mathrm{1,2,3,4};j=\mathrm{1,2},\·\·\·,k;t\∈\left[\mathrm{0,720000}\right];$

分别计算四个脑区任务态BOLD变化率基线值，其中为任务态时间序列t时第i个脑区第j个体素的BOLD值，为计算完成的任务态时第i个脑区BOLD变化率基线值。

静息态BOLD变化率计算模块34与所述BOLD值提取模块31和所述静息态BOLD变化率基线值计算模块32连接，用于根据静息态BOLD值和静息态BOLD变化率基线值获取静息态BOLD变化率。具体来说，静息态BOLD变化率计算模块34，根据以下公式：

$R_i^{\mathrm{rest\; j}}\left(t\right)=\frac{{\mathrm{bold}}_i^{\mathrm{rest\; j}}\left(t\right)-{\mathrm{bold}}_i^{\mathrm{rest}}}{{\mathrm{bold}}_i^{\mathrm{rest}}},i=\mathrm{1,2,3,4};j=\mathrm{1,2},\·\·\·k;t\∈\left[\mathrm{0,720000}\right];$

分别计算四个脑区静息态BOLD变化率，其中为静息态时间序列t时第i个脑区第j个体素的BOLD值，为静息态时第i个脑区BOLD变化率基线值，为计算完成的相应脑区静息态体素BOLD变化率。

任务态BOLD变化率计算模块35与所述BOLD值提取模块31和所述任务态BOLD变化率基线值计算模块33连接，用于根据任务态BOLD值和任务态BOLD变化率基线值获取任务态BOLD变化率。具体来说，任务态BOLD变化率计算模块35根据以下公式：

$R_i^{\mathrm{rest\; j}}\left(t\right)=\frac{{\mathrm{bold}}_i^{\mathrm{rest\; j}}\left(t\right)-{\mathrm{bold}}_i^{\mathrm{rest}}}{{\mathrm{bold}}_i^{\mathrm{rest}}},i=\mathrm{1,2,3,4};j=\mathrm{1,2},\·\·\·k;t\∈\left[\mathrm{0,720000}\right];$

分别计算四个脑区任务态BOLD变化率，其中为任务态时间序列t时第i个脑区第j个体素的BOLD值，为任务态时第i个脑区BOLD变化率基线值，为计算完成的相应脑区任务态体素BOLD变化率。

BOLD变化率曲线生成模块36与所述静息态BOLD变化率计算模块34和所述任务态BOLD变化率计算模块35连接，用于根据静息态BOLD变化率和任务态BOLD变化率生成静息态BOLD变化率曲线和任务态BOLD变化率曲线。BOLD变化率曲线生成模块36按照脑区顺序和时间序列，生成静息态BOLD变化率曲线和任务态BOLD变化率曲线，横坐标为时间序列，纵坐标为BOLD变化率值。

图6是根据本发明第二实施方式的变化强度值和变化矢量角生成装置的结构图。

如图6所示，变化强度值和变化矢量角生成装置4包括BOLD变化率曲线图生成模块41、变化向量计算模块42、变化强度值计算模块43和变化矢量角计算模块44。

BOLD变化率曲线图生成模块41用于根据静息态BOLD变化率曲线和任务态BOLD变化率曲线得到BOLD变化率曲线图。具体的，BOLD变化率曲线图生成模块41将静息态BOLD变化率曲线和任务态BOLD变化率曲线转化到同一坐标中，使其在X轴上具有相同的投影，得到相应的BOLD变化率曲线图。

变化向量计算模块42与所述变化向量计算模块42连接，用于对BOLD变化率曲线图进行向量分析以得到变化向量。具体来说，变化向量计算模块42对BOLD变化率曲线图进行向量分析，定量刻画曲线间的差异，计算出变化强度值||ΔG||和变化矢量角θ；BOLD变化率曲线中，任务态BOLD变化率曲线每相邻的两个坐标构成一个向量，四个脑区共得到4*k个坐标，根据时间序列构成4*k-1个向量，即构成向量集静息态BOLD变化率曲线也构成4*k-1个向量，即向量集静息态BOLD变化率曲线到任务态BOLD变化率曲线的变化，可由向量的变化表示，计算公式如下：

H为任务态BOLD变化率曲线构成的向量集，G为静息态BOLD变化率曲线构成的向量集，ΔG表示变化向量。

变化强度值计算模块43与所述变化强度值计算模块43连接，用于根据变化向量ΔG计算变化强度值||ΔG||，由欧氏距离决定，计算公式如下：

其中||ΔG||为变化强度值。

变化矢量角计算模块44与所述BOLD变化率曲线图生成模块41连接，用于根据向量夹角集计算变化矢量角。其中，所述向量夹角集是静息态BOLD变化率曲线和任务态BOLD变化率曲线同一时间序列段内两向量的夹角值构成的集合。

具体的，变化矢量角计算模块44计算任务态BOLD变化率曲线和静息态BOLD变化率曲线在同一时间序列段m时，两向量的夹角，计算公式如下：

其中和为向量集G和H中同一时间序列段m时对应的向量，θ_m为对应的两向量夹角；

根据上述计算公式，得到向量夹角集θ＝(θ₁，θ₂，…，θ_m，…θ_4*k-1)；

根据向量夹角集，计算变化矢量角θ，计算公式如下：

$\mathrm{\θ}=\frac{\underset{1}{\overset{4*k-1}{\mathrm{\Σ}}}{\mathrm{\θ}}_m}{4*k-1},m=\mathrm{1,2},\·\·\·,4*k-1;$

其中θ_m为静息态BOLD变化率曲线和任务态BOLD变化率曲线同一时间序列段m内两向量的夹角值，θ为计算完成的变化矢量角。

图7为一具体示例的抑郁高危蓝色预警阈值确定曲线图，图8为一具体示例的抑郁高危橙色预警阈值确定曲线图，图9为一具体示例的抑郁高危红色预警阈值确定曲线图，图10为一具体示例的抑郁高危极限预警阈值确定曲线图。

根据图7、图8、图9和图10中不同抑郁高危程度变化向量模的均值，确定出蓝色预警变化强度阈值橙色预警变化强度阈值红色预警变化强度阈值和预警极限变化强度阈值根据静息态和任务态同一时间序列段内两向量夹角的均值，确定出蓝色预警变化矢量角阈值橙色预警变化矢量角阈值红色预警变化矢量角阈值以及预警极限变化矢量角阈值

图11是本发明抑郁高危风险三级预警模型图。

其中，抑郁高危风险三级预警模型提供了抑郁高危预警所需的预警阈值及其分布情况，其抑郁高危判定方法与前述步骤S5的描述类似，在此不再赘述。

如上所述，本发明提供了一种抑郁风险三级预警方法及系统，通过本发明，可以对人体与抑郁症有关的生理指标进行定量的分析，实现了对不同程度抑郁高危人群的预警，并且预警稳定性好、正确率高。

应当理解的是，本发明的上述具体实施方式仅仅用于示例性说明或解释本发明的原理，而不构成对本发明的限制。因此，在不偏离本发明的精神和范围的情况下所做的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。此外，本发明所附权利要求旨在涵盖落入所附权利要求范围和边界、或者这种范围和边界的等同形式内的全部变化和修改例。

Claims (24)

1.一种抑郁风险三级预警方法，其特征在于，包括：

步骤S1，采集静息态和任务态的BOLD信号及其对应的时间序列；

步骤S2，从BOLD信号中提取静息态BOLD值和任务态BOLD值；

步骤S3，根据静息态BOLD值和任务态BOLD值生成静息态BOLD变化率曲线和任务态BOLD变化率曲线；

步骤S4，对静息态BOLD变化率曲线和任务态BOLD变化率曲线进行向量分析以得到变化强度值||ΔG||和变化矢量角θ；

步骤S5，根据变化强度值||ΔG||和变化矢量角θ与预设的变化强度阈值||ΔG||^阈值和变化矢量角阈值θ^阈值进行预警级别判定。

2.根据权利要求1所述的抑郁风险三级预警方法，其特征在于，所述步骤S3包括：

步骤S31,计算四个脑区静息态BOLD变化率基线值；

步骤S32,计算四个脑区任务态BOLD变化率基线值；

步骤S33,计算四个脑区静息态BOLD变化率；

步骤S34,计算四个脑区任务态BOLD变化率；以及

步骤S35,生成静息态BOLD变化率曲线和任务态BOLD变化率曲线。

3.根据权利要求2所述的抑郁风险三级预警方法，其特征在于，所述步骤S31根据以下公式分别计算四个脑区静息态BOLD变化率基线值：

i＝1,2,3,4；j＝1,2,…,k；t∈[0,720000]；

其中，为静息态时间序列t时第i个脑区第j个体素的BOLD 值，为计算完成的静息态时第i个脑区BOLD变化率基线值。

4.根据权利要求2所述的抑郁风险三级预警方法，其特征在于，所述步骤S32根据以下公式分别计算四个脑区任务态BOLD变化率基线值：

i＝1,2,3,4；j＝1,2,…,k；t∈[0,720000]；

其中，为任务态时间序列t时第i个脑区第j个体素的BOLD值，为计算完成的任务态时第i个脑区BOLD变化率基线值。

5.根据权利要求3所述的抑郁风险三级预警方法，其特征在于，所述步骤S33根据以下公式分别计算四个脑区静息态BOLD变化率：

i＝1,2,3,4；j＝1,2,…,k；t∈[0,720000]；

其中为静息态时间序列t时第i个脑区第j个体素的BOLD值，为静息态时第i个脑区BOLD变化率基线值，为计算完成的相应脑区静息态体素BOLD变化率。

6.根据权利要求4所述的抑郁风险三级预警方法，其特征在于，所述步骤S34根据以下公式分别计算四个脑区任务态BOLD变化率：

i＝1,2,3,4；j＝1,2,…,k；t∈[0,720000]；

其中为任务态时间序列t时第i个脑区第j个体素的BOLD值，为任务态时第i个脑区BOLD变化率基线值，为计算完成的相应脑区任务态体素BOLD变化率。

7.根据权利要求1所述的抑郁风险三级预警方法，其特征在于，所述步骤S4包括：

步骤S41,根据静息态BOLD变化率曲线和任务态BOLD变化率曲线得到BOLD变化率曲线图；

步骤S42,对BOLD变化率曲线图进行向量分析以得到变化向量；

步骤S43，根据变化向量计算变化强度值；以及

步骤S44，根据向量夹角集计算变化矢量角。

8.根据权利要求7所述的抑郁风险三级预警方法，所述步骤S42通过下式计算所述变化向量：

h为任务态BOLD变化率曲线每相邻的两个坐标构成的向量，g为静息态BOLD变化率曲线每相邻的两个坐标构成的向量，H为任务态BOLD变化率曲线构成的向量集，G为静息态BOLD变化率曲线构成的向量集，ΔG表示变化向量。

9.根据权利要求7所述的抑郁风险三级预警方法，所述步骤S43通过下式计算所述变化强度值：

其中，h为任务态BOLD变化率曲线每相邻的两个坐标构成的向量，g为静息态BOLD变化率曲线每相邻的两个坐标构成的向量，||ΔG||为变化强度值。

10.根据权利要求7所述的抑郁风险三级预警方法，所述步骤S44通过下式计算所述变化矢量角：

m＝1，2，…，4*k-1；

其中，θ_m为静息态BOLD变化率曲线和任务态BOLD变化率曲线同一时间序列段m内两向量的夹角值，θ为计算完成的变化矢量角。

11.根据权利要求1所述的抑郁风险三级预警方法，所述步骤S5根据预设的变化强度阈值||ΔG||^阈值和变化矢量角阈值θ^阈值进行预警级别的判断，包括：

若且则判定抑郁高危风险预警级别为蓝色预警；

若且则判定抑郁高危风险预警级别为橙色预警；

若且则判定抑郁高危风险预警级别为红色预警。

12.根据权利要求1-11中任一项所述的抑郁风险三级预警方法，所述BOLD信号为四个脑区的BOLD信号，所述四个脑区包括纹状体、岛叶、伏隔核和杏仁体。

13.一种抑郁风险三级预警系统，其特征在于包括：

BOLD信号采集装置(1)，用于采集静息态和任务态的BOLD信号及其对应的时间序列；

BOLD值提取装置(2)，用于从所述BOLD信号中提取静息态BOLD值和任务态BOLD值；

BOLD变化率曲线生成装置(3)，用于根据静息态BOLD值和任务态BOLD 值生成静息态BOLD变化率曲线和任务态BOLD变化率曲线；

变化强度值和变化矢量角生成装置(4)，用于对静息态BOLD变化率曲线和任务态BOLD变化率曲线进行向量分析以得到变化强度值||ΔG||和变化矢量角θ；和

预警级别判断装置(5)，用于根据变化强度值||ΔG||和变化矢量角θ与预设的变化强度阈值||ΔG||^阈值和变化矢量角阈值θ^阈值进行预警级别判定。

14.根据权利要求13所述的抑郁风险三级预警系统，其特征在于，所述BOLD变化率曲线生成装置(3)包括：

BOLD值提取模块(31)，用于从BOLD信号中提取静息态BOLD值和任务态BOLD值；

静息态BOLD变化率基线值计算模块(32)，用于根据静息态BOLD值获取静息态BOLD变化率基线值；

任务态BOLD变化率基线值计算模块(33)，用于根据任务态BOLD值获取任务态BOLD变化率基线值；

静息态BOLD变化率计算模块(34)，用于根据静息态BOLD值和静息态BOLD变化率基线值获取静息态BOLD变化率；

任务态BOLD变化率计算模块(35)，用于根据任务态BOLD值和任务态BOLD变化率基线值获取任务态BOLD变化率；和

BOLD变化率曲线生成模块(36)，用于根据静息态BOLD变化率和任务态BOLD变化率生成静息态BOLD变化率曲线和任务态BOLD变化率曲线。

15.根据权利要求14所述的抑郁风险三级预警系统，其特征在于，所述静息态BOLD变化率基线值计算模块(32)根据以下公式分别计算四 个脑区静息态BOLD变化率基线值：

i＝1,2,3,4；j＝1,2,…,k；t∈[0,720000]；

其中，为静息态时间序列t时第i个脑区第j个体素的BOLD值，为计算完成的静息态时第i个脑区BOLD变化率基线值。

16.根据权利要求14所述的抑郁风险三级预警系统，其特征在于，所述任务态BOLD变化率基线值计算模块(33)根据以下公式分别计算四个脑区任务态BOLD变化率基线值：

i＝1,2,3,4；j＝1,2,…,k；t∈[0,720000]；

其中，为任务态时间序列t时第i个脑区第j个体素的BOLD值，为计算完成的任务态时第i个脑区BOLD变化率基线值。

17.根据权利要求15所述的抑郁风险三级预警系统，其特征在于，所述静息态BOLD变化率计算模块(34)根据以下公式分别计算四个脑区静息态BOLD变化率：

i＝1,2,3,4；j＝1,2,…,k；t∈[0,720000]；

其中为静息态时间序列t时第i个脑区第j个体素的BOLD值，为静息态时第i个脑区BOLD变化率基线值，为计算完成的相应脑区静息态体素BOLD变化率。

18.根据权利要求15所述的抑郁风险三级预警系统，其特征在于，所述任务态BOLD变化率计算模块(35)根据以下公式分别计算四个脑区任务态BOLD变化率：

i＝1,2,3,4；j＝1,2,…,k；t∈[0,720000]；

其中为任务态时间序列t时第i个脑区第j个体素的BOLD值，为任务态时第i个脑区BOLD变化率基线值，为计算完成的相应脑区任务态体素BOLD变化率。

19.根据权利要求13所述的抑郁风险三级预警系统，其特征在于，所述变化强度值和变化矢量角生成装置(4)包括:

BOLD变化率曲线图生成模块(41)，用于根据静息态BOLD变化率曲线和任务态BOLD变化率曲线得到BOLD变化率曲线图；

变化向量计算模块(42)，用于对BOLD变化率曲线图进行向量分析以得到变化向量；

变化强度值计算模块(43)，用于根据变化向量ΔG计算变化强度值||ΔG||；和

变化矢量角计算模块(44)，用于根据向量夹角集计算变化矢量角。

20.根据权利要求19所述的抑郁风险三级预警系统，其特征在于，所述变化向量计算模块(42)通过下式计算所述变化向量：

h为任务态BOLD变化率曲线每相邻的两个坐标构成的向量，g为静息态BOLD变化率曲线每相邻的两个坐标构成的向量，H为任务态BOLD变化率曲线构成的向量集，G为静息态BOLD变化率曲线构成的向量集，ΔG表示变化向量。

21.根据权利要求19所述的抑郁风险三级预警系统，其特征在于，所述变化强度值计算模块(43)，通过下式计算所述变化强度值：

其中，h为任务态BOLD变化率曲线每相邻的两个坐标构成的向量，g为静息态BOLD变化率曲线每相邻的两个坐标构成的向量，||ΔG||为变化强度值。

22.根据权利要求19所述的抑郁风险三级预警系统，其特征在于，所述变化矢量角计算模块(44)，通过下式计算所述变化矢量角：

m＝1，2，…，4*k-1；

其中，θ_m为静息态BOLD变化率曲线和任务态BOLD变化率曲线同一时间序列段m内两向量的夹角值，θ为计算完成的变化矢量角。

23.根据权利要求13所述的抑郁风险三级预警系统，所述预警级别判断装置(5)根据预设的变化强度阈值||ΔG||^阈值和变化矢量角阈值θ^阈值进行预警级别的判断，其包括：

蓝色预警模块(51)，用于当且则判定抑郁高危风险预警级别为蓝色预警；

橙色预警模块(52)，用于当且则判定抑郁高危风险预警级别为橙色预警；和

红色预警模块(53)，用于当且 则判定抑郁高危风险预警级别为红色预警。

24.根据权利要求13-23中任一项所述的抑郁风险三级预警系统， 其特征在于，所述BOLD信号为四个脑区的BOLD信号，所述四个脑区包括纹状体、岛叶、伏隔核和杏仁体。