CN107518898B - 基于传感器阵列分解和波束成形的脑磁图源定位装置 - Google Patents
基于传感器阵列分解和波束成形的脑磁图源定位装置 Download PDFInfo
- Publication number
- CN107518898B CN107518898B CN201710672659.5A CN201710672659A CN107518898B CN 107518898 B CN107518898 B CN 107518898B CN 201710672659 A CN201710672659 A CN 201710672659A CN 107518898 B CN107518898 B CN 107518898B
- Authority
- CN
- China
- Prior art keywords
- matrix
- vector
- sensor array
- new
- decomposition
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Active
Links
Images
Classifications
-
- A—HUMAN NECESSITIES
- A61—MEDICAL OR VETERINARY SCIENCE; HYGIENE
- A61B—DIAGNOSIS; SURGERY; IDENTIFICATION
- A61B5/00—Measuring for diagnostic purposes; Identification of persons
- A61B5/05—Detecting, measuring or recording for diagnosis by means of electric currents or magnetic fields; Measuring using microwaves or radio waves
- A61B5/0522—Magnetic induction tomography
-
- A—HUMAN NECESSITIES
- A61—MEDICAL OR VETERINARY SCIENCE; HYGIENE
- A61B—DIAGNOSIS; SURGERY; IDENTIFICATION
- A61B5/00—Measuring for diagnostic purposes; Identification of persons
- A61B5/0033—Features or image-related aspects of imaging apparatus classified in A61B5/00, e.g. for MRI, optical tomography or impedance tomography apparatus; arrangements of imaging apparatus in a room
- A61B5/004—Features or image-related aspects of imaging apparatus classified in A61B5/00, e.g. for MRI, optical tomography or impedance tomography apparatus; arrangements of imaging apparatus in a room adapted for image acquisition of a particular organ or body part
- A61B5/0042—Features or image-related aspects of imaging apparatus classified in A61B5/00, e.g. for MRI, optical tomography or impedance tomography apparatus; arrangements of imaging apparatus in a room adapted for image acquisition of a particular organ or body part for the brain
-
- A—HUMAN NECESSITIES
- A61—MEDICAL OR VETERINARY SCIENCE; HYGIENE
- A61B—DIAGNOSIS; SURGERY; IDENTIFICATION
- A61B5/00—Measuring for diagnostic purposes; Identification of persons
- A61B5/05—Detecting, measuring or recording for diagnosis by means of electric currents or magnetic fields; Measuring using microwaves or radio waves
- A61B5/055—Detecting, measuring or recording for diagnosis by means of electric currents or magnetic fields; Measuring using microwaves or radio waves involving electronic [EMR] or nuclear [NMR] magnetic resonance, e.g. magnetic resonance imaging
-
- A—HUMAN NECESSITIES
- A61—MEDICAL OR VETERINARY SCIENCE; HYGIENE
- A61B—DIAGNOSIS; SURGERY; IDENTIFICATION
- A61B5/00—Measuring for diagnostic purposes; Identification of persons
- A61B5/40—Detecting, measuring or recording for evaluating the nervous system
- A61B5/4076—Diagnosing or monitoring particular conditions of the nervous system
- A61B5/4094—Diagnosing or monitoring seizure diseases, e.g. epilepsy
-
- A—HUMAN NECESSITIES
- A61—MEDICAL OR VETERINARY SCIENCE; HYGIENE
- A61B—DIAGNOSIS; SURGERY; IDENTIFICATION
- A61B5/00—Measuring for diagnostic purposes; Identification of persons
- A61B5/72—Signal processing specially adapted for physiological signals or for diagnostic purposes
- A61B5/7271—Specific aspects of physiological measurement analysis
-
- A—HUMAN NECESSITIES
- A61—MEDICAL OR VETERINARY SCIENCE; HYGIENE
- A61B—DIAGNOSIS; SURGERY; IDENTIFICATION
- A61B2576/00—Medical imaging apparatus involving image processing or analysis
- A61B2576/02—Medical imaging apparatus involving image processing or analysis specially adapted for a particular organ or body part
- A61B2576/026—Medical imaging apparatus involving image processing or analysis specially adapted for a particular organ or body part for the brain
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F2218/00—Aspects of pattern recognition specially adapted for signal processing
- G06F2218/22—Source localisation; Inverse modelling
Landscapes
- Health & Medical Sciences (AREA)
- Life Sciences & Earth Sciences (AREA)
- Physics & Mathematics (AREA)
- Engineering & Computer Science (AREA)
- Molecular Biology (AREA)
- General Health & Medical Sciences (AREA)
- Biophysics (AREA)
- Pathology (AREA)
- Nuclear Medicine, Radiotherapy & Molecular Imaging (AREA)
- Biomedical Technology (AREA)
- Heart & Thoracic Surgery (AREA)
- Medical Informatics (AREA)
- Veterinary Medicine (AREA)
- Surgery (AREA)
- Animal Behavior & Ethology (AREA)
- Public Health (AREA)
- Radiology & Medical Imaging (AREA)
- Neurology (AREA)
- Physiology (AREA)
- Neurosurgery (AREA)
- High Energy & Nuclear Physics (AREA)
- Artificial Intelligence (AREA)
- Computer Vision & Pattern Recognition (AREA)
- Psychiatry (AREA)
- Signal Processing (AREA)
- Measurement And Recording Of Electrical Phenomena And Electrical Characteristics Of The Living Body (AREA)
Abstract
本发明公开了一种基于传感器阵列分解和波束成形的脑磁图深度源定位装置,包括脑磁图传感器和定位单元,定位单元包括传感器阵列的迭代分解模块、向量波束成形估计器模块和深度源定位计算模块,传感器阵列的迭代分解模块应用迭代CP(CANDECOMP/PARAFAC)矩阵分解方法重构MEG传感器阵列信号;的向量波束成形估计器模块应用重构后的传感器阵列信号和样本的二阶统计量来估计总体的协方差矩阵;深度源定位计算模块应用投影向量的方差最小优化准则来估计深脑源的位置,本发明的技术方案受传感器信号噪声和干扰的影响小,对于深度源的定位精度很高,能较好地改进寻找颞叶内侧起源的癫痫患者的病灶区。
Description
技术领域
本发明属于生物医学信息领域中的成像技术和医学神经病学的病灶区定位的交叉技术领域,尤其是涉及一种基于传感器阵列分解和波束成形的脑磁图源定位装置。
背景技术
脑磁图是一种具有很高的时间分辨率和空间分辨率的功能性神经成像技术,且传感器获取的信号几乎不受头盖骨和头皮干扰的影响。因而利用这种高质量的传感器信号来反推大脑中信号源的位置,这一源定位装置和方法已广泛应用认知神经科学和神经病学等多个领域,其中一个非常重要的应用就是寻找癫痫患者的致痫灶。据现有资料显示,目前脑磁图在新皮层癫痫患者的定位中已取得较好的效果,然而,对于颞叶内侧或岛叶病灶区的癫痫患者来说,脑磁图的源定位方法几乎很难找到准确位置。脑磁图的深度源定位一直以来成为关注的焦点,也是尚未突破的技术难点。
近十多年来,已经有很多源定位方法被提出,其中波束成形方法为解决该问题的一类有效的方法。该类方法之所以有非常好的定位效果,主要应用到协方差矩阵估计和投影后向量的最小方差优化准则,然而利用原有传感器数据进行方差估计存着受噪声干扰等弊端。主要是由于大脑深部或较深部的源发出的信号被外部传感器接收后,这些信号原本就很弱,有些甚至不容易被观察到。因此,得到一种更稳健有效的源定位装置,尤其是深度源定位装置和方法以估计协方差矩阵成为能否定位大脑深部或较深部的源的关键。
发明内容
为解决上述技术问题,本发明提供了一种基于传感器阵列分解和波束成形的脑磁图源定位装置,不仅针对普通源定位,尤其针对大脑深部或较深部的源的定位有非常好的效果。
本发明完整的技术方案包括:
一种基于传感器阵列分解和波束成形的脑磁图源定位装置,尤其是基于传感器阵列分解和波束成形的脑磁图深度源定位装置,包括脑磁图传感器和定位单元,所述的定位单元包括传感器阵列的迭代分解模块、向量波束成形估计器模块和源定位计算模块,所述的传感器阵列的迭代分解模块应用迭代CP(CANDECOMP/PARAFAC)矩阵分解方法重构脑磁图(英文简写:MEG)传感器阵列(包含梯度计阵列和磁力计阵列)信号;所述的向量波束成形估计器模块应用重构后的传感器阵列信号和样本的二阶统计量来估计总体的协方差矩阵;所述的源定位计算模块应用投影向量的方差最小优化准则来估计处于深部脑区的源的位置。
验证模块,所述的验证模块在仿真数据和真实的癫痫病例数据上验证定位的准确性和有效性。
所述定位装置利用脑磁图传感器采集得到脑磁图传感器信号,并记为X=[x1,x2,...,xM],所述的定位单元根据传感器信号X来反推出源信号的位置,即反问题求解,该脑磁图传感器信号的反问题解模型表示为公式(1):
X=LD+ε (1)
其中xi表示在第i时刻点的MEG传感器测量值,xi为N×1向量,N为传感器数目,L表示一个N×J的前向解获得矩阵(lead-field matrix),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):
1)取原阵列信号的某一个列向量xi作为分数向量t的初始值:tstart=xi;
2)归一化该向量的能量为1:tnew=tstart/‖tstart‖;
3)计算一个新的加载向量pold:pold=XTt;此处的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,重新计算;
优选的,所述的向量波束成形估计器模块应用重构后的传感器阵列信号和样本的二阶统计量进行如下处理以估计总体的协方差矩阵:
采用样本的二阶中心矩来估计总体协方差矩阵,具体的表达形式见公式(4)和(5):
给出向量波束成形估计器的目标函数如下:
其中E(·)表示期望算子,W为权值矩阵,r0表示所有体素中的任意一个;
采用方差最小化得到优化后的权值矩阵W,优化后的结果如公式(7):
根据公式(8)计算出每个位置处的能量值,并根据所有位置点处能量值找出源的位置;
其中Pow(r)表示位置r处的能量值,tr{·}表示迹算子。
所述的定位单元根据脑磁图传感器采集到的脑磁图传感器信号进行如下处理:
(1)前向解计算:
1)基于一个单壳(single shell)近似的皮层约束下,构建出容积传导模型V;
2)基于皮层约束下计算引导场(lead field)矩阵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-field matrix),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):
1)取原阵列信号的某一个列向量xi作为分数向量t的初始值:tstart=xi;
2)归一化该向量的能量为1:tnew=tstart/‖tstart‖;
3)计算一个新的加载向量pold:pold=XTt;此处的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,重新计算;
所述的应用重构后的传感器阵列信号和样本的二阶统计量来估计总体的协方差矩阵具体操作方法为:
采用样本的二阶中心矩来估计总体协方差矩阵,具体的表达形式见公式(4)和(5):
给出向量波束成形估计器的目标函数如下:
其中E(·)表示期望算子,W为权值矩阵,r0表示所有体素中的任意一个;
采用方差最小化得到优化后的权值矩阵W,优化后的结果如公式(7):
根据公式(8)计算出每个位置处的能量值,并根据所有位置点处能量值找出源的位置;
其中Pow(r)表示位置r处的能量值,tr{·}表示迹算子。
上述基于传感器阵列分解和波束成形的源定位方法,具体操作方法为:
(1)前向解计算:
1)基于一个单壳(single shell)近似的皮层约束下,构建出容积传导模型V;
2)基于皮层约束下计算引导场(lead field)矩阵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为本发明所提的装置与临床上常用的dipole fitting(偶极子)装置在定位颞叶内侧结果的比较,横坐标轴为患者,纵坐标轴为每个患者能定位到颞叶内侧的棘波数目与棘波总数目的比值。
图6为第一个颞叶内侧起源的癫痫患者,使用dipole fitting方法的所有棘波定位结果,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-field matrix),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):
1)取原阵列信号的某一个列向量xi作为分数向量t的初始值:tstart=xi;
2)归一化该向量的能量为1:tnew=tstart/‖tstart‖;
3)计算一个新的加载向量pold:pold=XTt;此处的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,重新计算;
2.向量波束成形估计器模块
向量波束成形估计器主要利用一种线性约束的最小方差估计(LCMV)来优化目标函数。理论上,需要已知总体的协方差矩阵C(x),可在实际中经常用样本的二阶中心矩来估计总体协方差矩阵。具体的表达形式见公式(4)和(5):
接下来,给出向量波束成形估计器的目标函数如下:
其中E(·)表示期望算子,W为权值矩阵,r0表示所有体素中的任意一个。下面将采用方差最小化得到优化后的权值矩阵W,该方法也可简单描述成在r0位置允许单位能量通过,同时抑制所有其他位置源的能量贡献,在该优化公式中已经考虑到降低噪声ε对真实信号的影响,为源定位方法中的常规做法,在此不赘述,且优化后的结果如公式(7):
其中W(r)表示r位置处的投影矩阵,(·)-1表示逆算子,表示重构传感器阵列的协方差矩阵,L(r)表示r位置处的前向解。那么,我们将根据公式(8)计算出每个位置处的能量值,并根据所有位置点处能量值找出源的位置。
其中Pow(r)表示位置r处的能量值,tr{·}表示迹算子。
3.源定位算法模块
根据前两个模块的计算,下面将给出针对脑区深部或较深部的源定位的算法详细流程:
(1)前向解计算:
4)基于一个单壳(single shell)近似的皮层约束下,构建出容积传导模型V;
5)基于皮层约束下计算引导场(lead field)矩阵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例由临床综合诊断为颞叶内侧起源的癫痫患者,验证本发明提出的新方法,基于每个棘波计算出一个源,统计出每个患者能定位在颞叶内侧源的棘波数量所占的比例,并与现有临床上有效的偶极子拟合方法(dipole fitting)进行对比。图5-图7显示,新方法在这些患者上的病灶区定位效果上有一定的改善,通过这些有限病例可以看出,本发明可以为颞叶内侧癫痫患者提供更可靠的术前评估结果。
以上所述,仅是本发明的较佳实施例,并非对本发明作任何限制,凡是根据本发明技术实质对以上实施例所作的任何简单修改、变更以及等效结构变化,均仍属于本发明技术方案的保护范围内。
Claims (4)
1.一种基于传感器阵列分解和波束成形的脑磁图源定位装置,其特征在于,包括脑磁图传感器和定位单元,所述的定位单元包括传感器阵列的迭代分解模块、向量波束成形估计器模块和源定位计算模块,所述的传感器阵列的迭代分解模块应用迭代CP(CANDECOMP/PARAFAC)矩阵分解方法重构MEG传感器阵列信号;所述的向量波束成形估计器模块应用重构后的传感器阵列信号和样本的二阶统计量来估计总体的协方差矩阵;所述的源定位计算模块应用投影向量的方差最小优化准则来估计不同脑区的疾病病灶或功能区/功能亚区的源的位置;
所述的基于传感器阵列分解和波束成形的脑磁图源定位装置为基于传感器阵列分解和波束成形的脑磁图源深度源定位装置。
2.权利要求1所述的一种基于传感器阵列分解和波束成形的脑磁图源定位装置,其特征在于,还包括验证模块,所述的验证模块在仿真数据和真实的癫痫数据上验证定位的准确性。
3.权利要求1或2所述的一种基于传感器阵列分解和波束成形的脑磁图源定位装置,其特征在于,
所述定位装置利用脑磁图传感器采集得到脑磁图传感器信号,并记为X=[x1,x2,...,xM],所述的定位单元根据传感器信号X来反推出源信号的位置,即反问题求解,该脑磁图传感器信号的反问题解模型表示为公式(1):
X=LD+ε (1)
其中以xi表示在第i时刻点的MEG传感器测量值,xi为N×1向量,N为传感器数目,L表示N×J的前向解获得矩阵(lead-field matrix),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矩阵之和,以ti表示第i个分数向量,以pi表示第i个加载向量,每个秩1矩阵能够进一步表示为分数向量ti和加载向量pi的外积,其表示形式见公式(3):
1)取原阵列信号的某一个列向量xi作为分数向量t的初始值:tstart=xi;
2)归一化该向量的能量为1:tnew=tstart/‖tstart‖;
3)计算一个新的加载向量pold:pold=XTt;此处的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,重新计算;
所述的向量波束成形估计器模块应用重构后的传感器阵列信号和样本的二阶统计量进行如下处理以估计总体的协方差矩阵:
采用样本的二阶中心矩来估计总体协方差矩阵,具体的表达形式见公式(4)和(5):
给出向量波束成形估计器的目标函数如下:
其中E(·)表示期望算子,W为权值矩阵,r0表示所有体素中的任意一个;
采用方差最小化得到优化后的权值矩阵W,优化后的结果如公式(7):
根据公式(8)计算出每个位置处的能量值,并根据所有位置点处能量值找出源的位置;
其中Pow(r)表示位置r处的能量值,tr{·}表示迹算子。
4.权利要求3所述的一种基于传感器阵列分解和波束成形的脑磁图源定位装置,其特征在于,所述的定位单元根据脑磁图传感器采集到的脑磁图传感器信号进行如下处理:
(1)前向解计算:
1)基于一个单壳(single shell)近似的皮层约束下,构建出容积传导模型V;
2)基于皮层约束下计算引导场(lead field)矩阵L;
3)在第i个体素点上定位引导场矩阵Li;
(2)脑磁图传感器阵列分解:
4)挑选出第k个时间序列Xk,该时间序列包含棘波;
5)利用上述迭代的CP分解算法将数据矩阵Xk分解为分数矩阵Tk和加载矩阵Pk;
6)重构传感器阵列矩阵;
(3)反问题求解:
7)利用公式(4)估计协方差矩阵;
8)基于公式(7)获得r位置处的优化权值矩阵Wk(r);
9)利用公式(8)计算r位置处的能量值Powk(r);
(4)定位结果的呈现:
10)选择最大的能量值所对应的位置点作为最终确定的源位置;
11)将源定位的结果显示在个体磁共振(MRI)上。
Priority Applications (2)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN201710672659.5A CN107518898B (zh) | 2017-08-08 | 2017-08-08 | 基于传感器阵列分解和波束成形的脑磁图源定位装置 |
CN202010252977.8A CN111528794B (zh) | 2017-08-08 | 2017-08-08 | 基于传感器阵列分解和波束成形的源定位方法 |
Applications Claiming Priority (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN201710672659.5A CN107518898B (zh) | 2017-08-08 | 2017-08-08 | 基于传感器阵列分解和波束成形的脑磁图源定位装置 |
Related Child Applications (1)
Application Number | Title | Priority Date | Filing Date |
---|---|---|---|
CN202010252977.8A Division CN111528794B (zh) | 2017-08-08 | 2017-08-08 | 基于传感器阵列分解和波束成形的源定位方法 |
Publications (2)
Publication Number | Publication Date |
---|---|
CN107518898A CN107518898A (zh) | 2017-12-29 |
CN107518898B true CN107518898B (zh) | 2020-04-28 |
Family
ID=60680902
Family Applications (2)
Application Number | Title | Priority Date | Filing Date |
---|---|---|---|
CN202010252977.8A Active CN111528794B (zh) | 2017-08-08 | 2017-08-08 | 基于传感器阵列分解和波束成形的源定位方法 |
CN201710672659.5A Active CN107518898B (zh) | 2017-08-08 | 2017-08-08 | 基于传感器阵列分解和波束成形的脑磁图源定位装置 |
Family Applications Before (1)
Application Number | Title | Priority Date | Filing Date |
---|---|---|---|
CN202010252977.8A Active CN111528794B (zh) | 2017-08-08 | 2017-08-08 | 基于传感器阵列分解和波束成形的源定位方法 |
Country Status (1)
Country | Link |
---|---|
CN (2) | CN111528794B (zh) |
Families Citing this family (10)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN107550493B (zh) * | 2017-07-14 | 2020-09-08 | 北京大学 | 基于功能磁共振成像的时变约束脑电图或脑磁图溯源方法 |
CN108478208A (zh) * | 2018-03-22 | 2018-09-04 | 武汉大学 | 一种基于脑磁图的时空域双子空间投影的去干扰方法 |
CN108416822B (zh) * | 2018-03-22 | 2021-12-03 | 武汉大学 | 一种基于贝叶斯估计的多层次多尺度层析成像方法 |
CN109239207B (zh) * | 2018-07-23 | 2020-07-24 | 中山大学 | 基于电子鼻的气味识别方法、装置和电子鼻系统 |
CN112674773B (zh) * | 2020-12-22 | 2021-12-24 | 北京航空航天大学 | 基于Tucker分解和ripple时间窗的脑磁图源定位方法和装置 |
CN113995422B (zh) * | 2021-10-15 | 2022-11-25 | 苏州大学 | 基于非负块稀疏贝叶斯学习的瞬态脑电源定位方法及系统 |
CN113948189B (zh) * | 2021-12-22 | 2022-03-15 | 北京航空航天大学杭州创新研究院 | 基于gru神经网络的meg源定位方法 |
CN114376522B (zh) * | 2021-12-29 | 2023-09-05 | 四川大学华西医院 | 构建用于识别青少年肌阵挛癫痫的计算机识别模型的方法 |
CN114065825B (zh) * | 2022-01-17 | 2022-04-19 | 北京航空航天大学杭州创新研究院 | 一种基于结构相似性的脑磁meg源定位方法 |
CN115177256B (zh) * | 2022-09-13 | 2022-11-29 | 北京昆迈医疗科技有限公司 | 一种视频脑磁图系统 |
Citations (4)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
WO2009021308A1 (en) * | 2007-08-15 | 2009-02-19 | Hybrid Strategies Corporation | Magnetoelectric generator within double coaxial halbach cylinder |
CN105147288A (zh) * | 2015-07-23 | 2015-12-16 | 中国科学院苏州生物医学工程技术研究所 | 脑磁源强度定位方法 |
WO2016073985A1 (en) * | 2014-11-07 | 2016-05-12 | The General Hospital Corporation | Deep brain source imaging with m/eeg and anatomical mri |
CN105596000A (zh) * | 2016-02-03 | 2016-05-25 | 王小姗 | 一种脑磁图分析装置及方法 |
Family Cites Families (5)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
US20120083647A1 (en) * | 2010-09-30 | 2012-04-05 | Harry Scheinin | Method for changing an individual's state of consciousness |
US10531806B2 (en) * | 2013-12-17 | 2020-01-14 | University Of Florida Research Foundation, Inc. | Brain state advisory system using calibrated metrics and optimal time-series decomposition |
US20180122507A1 (en) * | 2015-04-14 | 2018-05-03 | University Of Utah Research Foundation | Genetic alterations in ovarian cancer |
TW201711634A (zh) * | 2015-09-17 | 2017-04-01 | 國立中央大學 | 腦功能成像辨識的方法及系統 |
CN106539582A (zh) * | 2015-09-17 | 2017-03-29 | 中央大学 | 脑功能成像辨识的方法及系统 |
-
2017
- 2017-08-08 CN CN202010252977.8A patent/CN111528794B/zh active Active
- 2017-08-08 CN CN201710672659.5A patent/CN107518898B/zh active Active
Patent Citations (4)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
WO2009021308A1 (en) * | 2007-08-15 | 2009-02-19 | Hybrid Strategies Corporation | Magnetoelectric generator within double coaxial halbach cylinder |
WO2016073985A1 (en) * | 2014-11-07 | 2016-05-12 | The General Hospital Corporation | Deep brain source imaging with m/eeg and anatomical mri |
CN105147288A (zh) * | 2015-07-23 | 2015-12-16 | 中国科学院苏州生物医学工程技术研究所 | 脑磁源强度定位方法 |
CN105596000A (zh) * | 2016-02-03 | 2016-05-25 | 王小姗 | 一种脑磁图分析装置及方法 |
Also Published As
Publication number | Publication date |
---|---|
CN107518898A (zh) | 2017-12-29 |
CN111528794B (zh) | 2021-03-19 |
CN111528794A (zh) | 2020-08-14 |
Similar Documents
Publication | Publication Date | Title |
---|---|---|
CN107518898B (zh) | 基于传感器阵列分解和波束成形的脑磁图源定位装置 | |
Meyer et al. | Flexible head-casts for high spatial precision MEG | |
Baillet et al. | Electromagnetic brain mapping | |
Heers et al. | Localization accuracy of distributed inverse solutions for electric and magnetic source imaging of interictal epileptic discharges in patients with focal epilepsy | |
Albera et al. | ICA-based EEG denoising: a comparative analysis of fifteen methods | |
US8032209B2 (en) | Localizing neural sources in a brain | |
Cai et al. | Hierarchical multiscale Bayesian algorithm for robust MEG/EEG source reconstruction | |
WO2012162569A2 (en) | Magnetoencephalography source imaging | |
US20160051161A1 (en) | Method for locating a brain activity associated with a task | |
Jonmohamadi et al. | Comparison of beamformers for EEG source signal reconstruction | |
US20120232376A1 (en) | Methods and systems for channel selection | |
CN116887752A (zh) | 脑磁图方法和系统 | |
Chang et al. | Assessing recurrent interactions in cortical networks: modeling EEG response to transcranial magnetic stimulation | |
Korats et al. | A space-time-frequency dictionary for sparse cortical source localization | |
Cetin et al. | Multimodal based classification of schizophrenia patients | |
Ramírez | Source localization | |
PT1005677E (pt) | Sistema e método para a tomografia da corrente eléctrica primária do cérebro e do coração | |
Belaoucha et al. | Structural connectivity to reconstruct brain activation and effective connectivity between brain regions | |
Mideksa et al. | Source analysis of median nerve stimulated somatosensory evoked potentials and fields using simultaneously measured EEG and MEG signals | |
Owen et al. | Estimating the location and orientation of complex, correlated neural activity using MEG | |
Fujiwara et al. | A hierarchical Bayesian method to resolve an inverse problem of MEG contaminated with eye movement artifacts | |
Gómez et al. | Localization accuracy of a common beamformer for the comparison of two conditions | |
Ramírez et al. | Neuroelectromagnetic source imaging of brain dynamics | |
CN117653030A (zh) | 基于tde-hmm模型的癫痫病灶溯源定位系统、方法、终端及介质 | |
Wan et al. | Electromagnetic source imaging: Backus–Gilbert resolution spread function‐constrained and functional MRI‐guided spatial filtering |
Legal Events
Date | Code | Title | Description |
---|---|---|---|
PB01 | Publication | ||
PB01 | Publication | ||
SE01 | Entry into force of request for substantive examination | ||
SE01 | Entry into force of request for substantive examination | ||
GR01 | Patent grant | ||
GR01 | Patent grant |