基于传感器阵列分解和波束成形的源定位方法与流程

本发明是2017年08月08日递交的、申请号为201710672659.5(母案)、发明名称为“基于传感器阵列分解和波束成形的脑磁图源定位方法”的发明专利申请的分案申请。

本发明属于生物医学信息领域中的成像技术和医学神经病学的病灶区定位的交叉技术领域，尤其是涉及一种基于传感器阵列分解和波束成形的脑磁图源定位装置。

背景技术：

脑磁图是一种具有很高的时间分辨率和空间分辨率的功能性神经成像技术，且传感器获取的信号几乎不受头盖骨和头皮干扰的影响。因而利用这种高质量的传感器信号来反推大脑中信号源的位置，这一源定位装置和方法已广泛应用认知神经科学和神经病学等多个领域，其中一个非常重要的应用就是寻找癫痫患者的致痫灶。据现有资料显示，目前脑磁图在新皮层癫痫患者的定位中已取得较好的效果，然而，对于颞叶内侧或岛叶病灶区的癫痫患者来说，脑磁图的源定位方法几乎很难找到准确位置。脑磁图的深度源定位一直以来成为关注的焦点，也是尚未突破的技术难点。

近十多年来，已经有很多源定位方法被提出，其中波束成形方法为解决该问题的一类有效的方法。该类方法之所以有非常好的定位效果，主要应用到协方差矩阵估计和投影后向量的最小方差优化准则，然而利用原有传感器数据进行方差估计存着受噪声干扰等弊端。主要是由于大脑深部或较深部的源发出的信号被外部传感器接收后，这些信号原本就很弱，有些甚至不容易被观察到。因此，得到一种更稳健有效的源定位装置，尤其是深度源定位装置和方法以估计协方差矩阵成为能否定位大脑深部或较深部的源的关键。

技术实现要素：

为解决上述技术问题，本发明提供了一种基于传感器阵列分解和波束成形的脑磁图源定位装置，不仅针对普通源定位，尤其针对大脑深部或较深部的源的定位有非常好的效果。

本发明完整的技术方案包括：

一种基于传感器阵列分解和波束成形的脑磁图源定位装置，尤其是基于传感器阵列分解和波束成形的脑磁图深度源定位装置，包括脑磁图传感器和定位单元，所述的定位单元包括传感器阵列的迭代分解模块、向量波束成形估计器模块和源定位计算模块，所述的传感器阵列的迭代分解模块应用迭代cp(candecomp/parafac)矩阵分解方法重构脑磁图(英文简写：meg)传感器阵列(包含梯度计阵列和磁力计阵列)信号；所述的向量波束成形估计器模块应用重构后的传感器阵列信号和样本的二阶统计量来估计总体的协方差矩阵；所述的源定位计算模块应用投影向量的方差最小优化准则来估计处于深部脑区的源的位置。

验证模块，所述的验证模块在仿真数据和真实的癫痫病例数据上验证定位的准确性和有效性。

所述定位装置利用脑磁图传感器采集得到脑磁图传感器信号，并记为x＝[x1，x2，…，xm]，所述的定位单元根据传感器信号x来反推出源信号的位置，即反问题求解，该脑磁图传感器信号的反问题解模型表示为公式(1)：

x＝ld+ε(1)

具中xi表示在第i时刻点的meg传感器测量值，xi为n×1向量，n为传感器数目，l表示一个n×j的前向解获得矩阵(lead-fieldmatrix)，j表示未知的偶极子矩参数，d是一个在给定时间序列上的j×m偶极子矩矩阵，ε表示n×m的噪声矩阵，反问题求解过程即为根据已知的传感器阵列信号x求出解d；

优选的，所述的传感器阵列的迭代分解模块对meg传感器阵列信号进行如下重构处理：

应用秩1的矩阵分解将信号从一个空间投影到另一个空间，具体如公式(2)和公式(3)描述，将一个原信号矩阵近似分解为r个秩1矩阵之和：

x≈k1+k2+k3+…+kr(2)

其中r是一个正整数，且ki是第i个秩1矩阵，即该式表示为一个原信号矩阵近似分解为r个秩1矩阵之和，同时，每个秩1矩阵可以进一步表示成分数向量ti和加载向量pi的外积，其表示形式见公式(3)：

其中符号表示外积算子，ti和pi为两个列向量，t表示矩阵或向量的转置，具体的迭代过程描述如下：

1)取原阵列信号的某一个列向量xi作为分数向量t的初始值：tstart＝xi；

2)归一化该向量的能量为1：tnew＝tstart/||tstart||；

3)计算一个新的加载向量pold：pold＝x^tt；此处的t为分数向量；

4)归一化该加载向量的能量为1：pnew＝pold/||pold||，此处的pold等于步骤3)中的pold；

5)计算一个新的分数向量told：told＝xpnew，此处的pnew等于步骤4)中的pnew；

6)归一化该新的分数向量的能量为1：t’new＝told/||told||，此处的told等于步骤5)中的told；

7)比较步骤6)计算出的新分数向量t’new与第2)步的分数向量tnew，若收敛，则转到步骤8)，若不收敛，则转到步骤3)继续计算直到收敛；

8)计算残差矩阵e：e＝x-t’newpnew^t，且r的值由残差矩阵e的范数决定，即当||e||＞0.1，继续，当||e||≤0.1，重新计算；

通过上述步骤的计算，由每一步得到的t’new得出原阵列信号的本质特征t，并重构出新的传感器阵列信号

优选的，所述的向量波束成形估计器模块应用重构后的传感器阵列信号和样本的二阶统计量进行如下处理以估计总体的协方差矩阵：

采用样本的二阶中心矩来估计总体协方差矩阵，具体的表达形式见公式(4)和(5)：

其中公式(5)为样本均值，应用迭代的cp分解重构出的传感器阵列矩阵来替换公式(4)中的阵列矩阵χ。

给出向量波束成形估计器的目标函数如下：

且wl(r0)＝1(6)

其中e(·)表示期望算子，w为权值矩阵，r0表示所有体素中的任意一个；

采用方差最小化得到优化后的权值矩阵w，优化后的结果如公式(7)：

其中w(r)表示r位置处的投影矩阵，(·)^-1表示逆算子，表示重构传感器阵列的协方差矩阵，l(r)表示r位置处的前向解；

根据公式(8)计算出每个位置处的能量值，并根据所有位置点处能量值找出源的位置；

其中pow(r)表示位置r处的能量值，tr{·}表示迹算子。

所述的定位单元根据脑磁图传感器采集到的脑磁图传感器信号进行如下处理：

(1)前向解计算：

1)基于一个单壳(singleshell)近似的皮层约束下，构建出容积传导模型v；

2)基于皮层约束下计算引导场(leadfield)矩阵l；

3)在第i个体素点上定位引导场矩阵li；

(2)脑磁图传感器阵列分解：

4)挑选出第k个时间序列xk，，该时间序列包含棘波，且时间窗的长度约为200ms；

5)利用上述迭代的cp分解算法将数据矩阵xk分解为分数矩阵tk和加载矩阵pk；

6)重构传感器阵列矩阵；

(3)反问题求解：

7)利用公式(4)估计协方差矩阵；

8)基于公式(7)获得r位置处的优化权值矩阵wk(r)；

9)利用公式(8)计算r位置处的能量值powk(r)；

(4)定位结果的呈现：

10)选择最大的能量值所对应的位置点作为最终确定的源位置；

11)将源定位的结果显示在个体磁共振(mri)上。

此外，本发明还提供了一种基于传感器阵列分解和波束成形的源定位方法，所述的方法包括如下步骤：

(1)应用迭代cp(candecomp/parafac)矩阵分解技术重构meg传感器阵列信号；

(2)应用重构后的传感器阵列信号和样本的二阶统计量来估计总体的协方差矩阵；

(3)应用投影向量的方差最小优化准则来估计深脑源的位置。

(4)在仿真数据和真实数据上验证算法的可行性。

上述定位方法，对于观察到的传感器信号x＝[x1，x2，...，xm]，根据传感器信号来反推出源信号的位置，该过程称为反问题求解，该传感器信号的反问题解模型表示为公式(1)：

x＝ld+ε(1)

其中xi表示在第i时刻点的meg传感器测量值，xi为n×1向量，n为传感器数目，l表示一个n×j的前向解获得矩阵(lead-fieldmatrix)，j表示未知的偶极子矩参数，d是一个在给定时间序列上的j×m偶极子矩矩阵，ε表示n×m的噪声矩阵，反问题求解过程即为根据已知的传感器阵列信号x求出解d；

所述的应用迭代cp(candecomp/parafac)矩阵分解技术重构meg传感器阵列信号具体操作方法为：

应用秩1的矩阵分解将信号从一个空间投影到另一个空间，具体如公式(2)和公式(3)描述，将一个原信号矩阵近似分解为r个秩1矩阵之和：

x≈k1+k2+k3+...+kr(2)

其中r是一个正整数，且ki是第i个秩1矩阵，即该式表示为一个原信号矩阵近似分解为r个秩1矩阵之和，同时，每个秩1矩阵可以进一步表示成分数向量ti和加载向量pi的外积，其表示形式见公式(3)：

其中符号表示外积算子，ti和pi为两个列向量，t表示矩阵或向量的转置，具体的迭代过程描述如下：

1)取原阵列信号的某一个列向量xi作为分数向量t的初始值：tstart＝xi；

2)归一化该向量的能量为1：tnew＝tstart/||tstart||；

3)计算一个新的加载向量po1d：pold＝x^tt；此处的t为分数向量；

4)归一化该加载向量的能量为1：pnew＝pold/||pold||，此处的pold等于步骤3)中的pold；

5)计算一个新的分数向量told：told＝xpnew，此处的pnew等于步骤4)中的pnew；

6)归一化该新的分数向量的能量为1：t’new＝told/||told||，此处的told等于步骤5)中的told；

7)比较步骤6)计算出的新分数向量t’new与第2)步的分数向量tnew，若收敛，则转到步骤8)，若不收敛，则转到步骤3)继续计算直到收敛；

8)计算残差矩阵e：e＝x-t’newpnew^t，且r的值由残差矩阵e的范数决定，即当||e||＞0.1，继续，当||e||≤0.1，重新计算；

通过上述步骤的计算，由每一步得到的t’new得出原阵列信号的本质特征t，并重构出新的传感器阵列信号

所述的应用重构后的传感器阵列信号和样本的二阶统计量来估计总体的协方差矩阵具体操作方法为：

采用样本的二阶中心矩来估计总体协方差矩阵，具体的表达形式见公式(4)和(5)：

其中公式(5)为样本均值，应用迭代的cp分解重构出的传感器阵列矩阵来替换公式(4)中的阵列矩阵χ。

给出向量波束成形估计器的目标函数如下：

且wl(r0)＝1(6)

其中e(·)表示期望算子，w为权值矩阵，r0表示所有体素中的任意一个；

采用方差最小化得到优化后的权值矩阵w，优化后的结果如公式(7)：

其中w(r)表示r位置处的投影矩阵，(·)^-1表示逆算子，表示重构传感器阵列的协方差矩阵，l(r)表示r位置处的前向解；

根据公式(8)计算出每个位置处的能量值，并根据所有位置点处能量值找出源的位置；

其中pow(r)表示位置r处的能量值，tr{·}表示迹算子。

上述基于传感器阵列分解和波束成形的源定位方法，具体操作方法为：

(1)前向解计算：

1)基于一个单壳(singleshell)近似的皮层约束下，构建出容积传导模型v；

2)基于皮层约束下计算引导场(leadfield)矩阵l；

3)在第i个体素点上定位引导场矩阵li；

(2)传感器阵列分解：

4)挑选出第k个时间序列xk，，该时间序列包含棘波，且时间窗的长度约为200ms；

5)利用上述迭代的cp分解算法将数据矩阵xk分解为分数矩阵tk和加载矩阵pk.

6)重构传感器阵列矩阵；

(3)反问题求解：

7)利用公式(4)估计协方差矩阵；

8)基于公式(7)获得r位置处的优化权值矩阵wk(r)；

9)利用公式(8)计算r位置处的能量值powk(r)；

(4)定位结果的呈现：

10)选择最大的能量值所对应的位置点作为最终确定的源位置；

11)将源定位的结果显示。

本发明相对于现有技术的优点在于：

(1)本发明受传感器信号噪声和干扰的影响小，这主要因为本发明应用了改进的波束成形方法，而这一方法的核心应用了迭代的二维张量分解技术，该技术能很好的提取脑磁图传感器阵列信号的本质特征，这些特征已经很大程度上排除了噪声的干扰。(2)本发明对于大脑深部或较深部的源的定位精度很高，特别在模拟源的精度上达到小于4mm，在浅源定位精度达到小于2mm.这是因为采用了迭代的cp分解技术重构传感器阵列，再结合线性约束的最小方差优化准则使得定位精度得到更多的提高。(3)因为真实的脑磁图信号容易受噪声干扰，相比临床上现有的脑磁图的源定位方法就是偶极子拟合模型(dipolefitting)易受噪声干扰且对脑区深部或较深部的源的定位效果很差的问题，本发明能较好地改进寻找颞叶内侧起源的癫痫患者的病灶区。

附图说明

图1为矩阵x的秩1分解即cp分解。

图2为仿真的传感器数据信噪比在12种不同噪声水平下的均值。

图3为通过传感器信号显示的真实信号、噪声信号和仿真信号。(a)真实信号和噪声信号，且时间窗为600ms.(b)由余弦振荡的真实信号和高斯白噪声的噪声信号(信噪比为0.459)叠加生成的真实信号

图4为本发明所采用的装置、lcmv、music和cp+music四种不同的源定位装置分别在六种不同位置的定位结果比较。

每个子图的横坐标都表示不同的噪声水平，纵坐标表示在相同噪声水平下随机100次的空间精度平均值。(a)当真实源在右额叶位置处，基于四种不同方法定位后的空间精度的比较；(b)当真实源在右颞叶外侧位置处，基于四种不同方法定位后的空间精度的比较；(c)当真实源在右顶叶位置处，基于四种不同方法定位后的空间精度的比较；(d)当真实源在右枕叶位置处，基于四种不同方法定位后的空间精度的比较；(e)当真实源在右颞叶内侧(右海马)位置处，基于四种不同方法定位后的空间精度的比较；(f)当真实源在左颞叶内侧(左海马)位置处，基于四种不同方法定位后的空间精度的比较。

图5为本发明所提的装置与临床上常用的dipolefitting(偶极子)装置在定位颞叶内侧结果的比较，横坐标轴为患者，纵坐标轴为每个患者能定位到颞叶内侧的棘波数目与棘波总数目的比值。

图6为第一个颞叶内侧起源的癫痫患者，使用dipolefitting方法的所有棘波定位结果，r表示右侧，l表示左侧，白色圆框区域表示为海马区。

图7为第一个颞叶内侧起源的癫痫患者，应用本发明提出的装置，基于某一个棘波的定位结果。(a)该子图表示一系列矢状位切片图，将大脑从左向右切片得到该图的自上而下且自右向左的顺序，其中r表示大脑右半球，白色圆框区域为定位的病灶区；(b)该子图表示一系列冠状位切片图，将大脑从前向后切片得到该图的自上而下且自左向右的顺序，其中r表示右侧，l表示左侧，白色圆框区域为定位的病灶区。

具体实施方式

下面结合附图和具体实施方式对本发明做进一步说明。

该部分将详细地描述本发明的具体实施方案，本发明公开的基于传感器阵列分解和波束成形的脑磁图深度源定位装置包括脑磁图传感器和定位单元，所述的脑磁图传感器采集传感器信号，所述的定位单元主要包括传感器阵列的迭代cp(candecomp/parafac)分解、向量波束成形估计器、源定位算法三个模块。

1.传感器阵列的迭代分解模块

对于观察到的脑磁图传感器信号x＝[x1，x2，...，xm]，源定位的核心就是根据传感器信号来反推出源信号的位置，这个过程被称为反问题求解。为了更清晰地表达该问题，该信号的反问题解模型被表示为公式(1)：

x＝ld+ε(1)

其中xi表示在第i时刻点的meg传感器测量值，假设是一个n×1向量，n为传感器数目，l表示一个n×j的前向解获得矩阵(lead-fieldmatrix)，j表示未知的偶极子矩参数，d是一个在给定时间序列上的j×m偶极子矩矩阵，并且ε表示n×m的噪声矩阵。那么，反问题求解过程的核心就是根据已知的传感器阵列信号x求出解d。

鉴于原传感器阵列信号很容易受噪声或其他因素的干扰，这将直接导致定位结果的不精确和定位失效。特别地，深部源信号经过传输到达传感器后，对噪声更敏感。因此，为了更准确地定位脑区深部或较深部的源，必须要考虑无失真地去除原信号的干扰。接下来，本发明将应用一种二维张量分解技术来提取原信号的本质特征，并将信号从一个空间投影到另一个空间。该张量分解技术为秩1的矩阵分解，其图示表达可见图1。cp张量分解已经成功应用在很多领域中，并成功恢复出更干净的传感器阵列信号，提取出本质特征。为了详细地描述该分解，通过公式(2)和公式(3)来表达，具体描述如下。cp分解将一个原信号矩阵近似分解为r个秩1矩阵之和。

x≈k1+k2+k3+...+kr(2)

其中r是一个正整数，且ki是第i个秩1矩阵，即该式表示为一个原信号矩阵近似分解为r个秩1矩阵之和。同时，每个秩1矩阵可以进一步表示成分数向量ti和加载向量pi的外积，其表示形式见公式(3)：

其中符号表示外积算子，ti和pi为两个列向量，t表示矩阵或向量的转置。下面将采用一种鲁棒的迭代分解过程来求解，具体的迭代过程描述如下：

1)取原阵列信号的某一个列向量xi作为分数向量t的初始值：tstart＝xi；

2)归一化该向量的能量为1：tnew＝tstart/||tstart||；

3)计算一个新的加载向量pold：pold＝x^tt；此处的t为分数向量；

4)归一化该加载向量的能量为1：pnew＝pold/||pold||，此处的pold等于步骤3)中的pold；

5)计算一个新的分数向量told：told＝xpnew，此处的pnew等于步骤4)中的pnew；

6)归一化该新的分数向量的能量为1：t’new＝told/||told||，此处的told等于步骤5)中的told；

7)比较步骤6)计算出的新分数向量t’new与第2)步的分数向量tnew，若收敛，则转到步骤8)，若不收敛，则转到步骤3)继续计算直到收敛；

8)计算残差矩阵e：e＝x-t’newpnew^t，且r的值由残差矩阵e的范数决定，即当||e||＞0.1，继续，当||e||≤0.1，重新计算；

通过上述步骤的计算，由每一步得到的t’new得出原阵列信号的本质特征t，并重构出新的传感器阵列信号

2.向量波束成形估计器模块

向量波束成形估计器主要利用一种线性约束的最小方差估计(lcmv)来优化目标函数。理论上，需要已知总体的协方差矩阵c(x)，可在实际中经常用样本的二阶中心矩来估计总体协方差矩阵。具体的表达形式见公式(4)和(5)：

其中公式(5)为样本均值，并应用迭代的cp分解重构出的传感器阵列矩阵来替换公式(4)中的阵列矩阵χ。

接下来，给出向量波束成形估计器的目标函数如下：

且wl(r0)＝1(6)

其中e(·)表示期望算子，w为权值矩阵，r0表示所有体素中的任意一个。下面将采用方差最小化得到优化后的权值矩阵w，该方法也可简单描述成在r0位置允许单位能量通过，同时抑制所有其他位置源的能量贡献，在该优化公式中已经考虑到降低噪声ε对真实信号的影响，为源定位方法中的常规做法，在此不赘述，且优化后的结果如公式(7)：

其中w(r)表示r位置处的投影矩阵，(·)^-1表示逆算子，表示重构传感器阵列的协方差矩阵，l(r)表示r位置处的前向解。那么，我们将根据公式(8)计算出每个位置处的能量值，并根据所有位置点处能量值找出源的位置。

其中pow(r)表示位置r处的能量值，tr{·}表示迹算子。

3.源定位算法模块

根据前两个模块的计算，下面将给出针对脑区深部或较深部的源定位的算法详细流程：

(1)前向解计算：

4)基于一个单壳(singleshell)近似的皮层约束下，构建出容积传导模型v；

5)基于皮层约束下计算引导场(leadfield)矩阵l；

6)在第i个体素点上定位引导场矩阵li；

(2)脑磁图传感器阵列分解：

4)挑选出第k个时间序列xk，若在癫痫患者中，该时间序列包含棘波，且时间窗的长度大约为200ms；

5)利用上述迭代的cp分解算法将数据矩阵xk分解为分数矩阵tk和加载矩阵pk.

6)重构传感器阵列矩阵；

(3)反问题求解：

7)利用公式(4)估计协方差矩阵；

8)基于公式(7)获得r位置处的优化权值矩阵wk(r)；

9)利用公式(8)计算r位置处的能量值powk(r)，并以所求到的能量值powk(r)表达d；

(4)定位结果的呈现：

10)选择最大的能量值所对应的位置点作为最终确定的源位置；

11)将源定位的结果显示在个体磁共振(mri)上。

4.本发明验证模块

根据第3个模块提出的源定位算法，本模块将分别在模拟数据源和真实的癫痫病例数据中验证该算法的定位有效性和定位精度。

基于fieldtrip工具包生成的模拟源，如图2和图3所示，分别由颞叶外侧、额叶、顶叶、枕叶、颞叶内侧位置的源生成的模拟信号，且该模拟信号含有不同信噪比的噪声干扰。将本发明提出的方法与效果很好的线性约束的最小方差法(lcmv)、多重信号分类(music)、cp分解结合多重信号分类(cp+music)三种源定位方法比较，图5表明新方法的空间精度最高，浅源的空间定位精度小于2mm，深源的空间定位精度小于4mm。

此外，针对10例由临床综合诊断为颞叶内侧起源的癫痫患者，验证本发明提出的新方法，基于每个棘波计算出一个源，统计出每个患者能定位在颞叶内侧源的棘波数量所占的比例，并与现有临床上有效的偶极子拟合方法(dipolefitting)进行对比。图5-图7显示，新方法在这些患者上的病灶区定位效果上有一定的改善，通过这些有限病例可以看出，本发明可以为颞叶内侧癫痫患者提供更可靠的术前评估结果。

以上所述，仅是本发明的较佳实施例，并非对本发明作任何限制，凡是根据本发明技术实质对以上实施例所作的任何简单修改、变更以及等效结构变化，均仍属于本发明技术方案的保护范围内。