当前位置: 首页 > 期刊 > 《生物医学工程学杂志》 > 1999年第4期
编号:10499251
脑磁源的定位计算的最优化
http://www.100md.com 《生物医学工程学杂志》 1999年第4期
     作者:赵双任 蒋大宗

    单位:(西安交通大学,西安 710049)

    关键词:脑磁图;最小模;偶极子;权重函数

    生物医学工程学杂志990414 摘要 脑磁源的定位是脑磁图(Magnetoencephalography,MEG)研究的一个基本问题。该问题属不定问题,即根据探头测量的微小磁场所建立的方程组求得的未知脑磁源有无穷组解。因此需要通过合理的数学模型补充适当信息。目前的两种主要的数学模型为偶极子源搜索法和Lp最小模法。对前者若采用全局搜索费时太多,若采用梯度类搜索法又容易收敛到局部极值点。故一般得用其他方法引导,定出初值,再使用梯度类搜索法。后者算得脑磁源的深度依赖于该方法中权重函数的选择,而该权重函数又有不同的几种,因此不能准确的估计脑磁源的深度。为了克服上述两种方法的困难,采用多目标函数最优化法。该方法的目标函数为,(1)未知源的Lp模,(2)由未知源计算出的磁场与测量的磁场之间的偏差,(3)源的数目。该方法与最小模方法相似也引入权重矩阵。但权重矩阵含多个参数,这些参数同未知源一样也要通过三目标最优化法来确定。适当选用权重函数及其参数,该方法可以准确的估计源的深度
, 百拇医药
    Best Minimization to Locate the Brain Sources of Magnetic Waves

    Zhao Shuangren Jiang Dazong

    (Xi'an Jiaotong University, Xi'an 710049)

    Abstract Locating the sources of brain magnetic fields is a basic question in magnetoencephalography (MEG). This question belongs to the indeterminate problem, that is, infinite solutions can be obtained from the equations of the measured small magnetic fields on the detector. Therefore a model of mathematics is needed to give more information. The major models of mathematics are the search method with many dipoles and the Lp minimum norm method. For the first model, if a global search is used, too much calculation time will be spent; if the gradient search is used, the searches usually converge to a local minimun. So other methods are introduced to find out the initial values, and then the gradient search can be used. For the second model, the depth of the brain sources will dependend on the choice of the weight functions which have several different types, and we do not know in advance which type is correct, so that it is difficult to have a good estimation of the depth of the brain sources. To overcome the difficulties, we use the method with three object functions which are the minimum Lp norm, the deviation between the calculated magnetic fields and the measured magnetic fields, and the number of the sources. This method is similar to the minimum norm method, a weight function is also introduced. But this weight function contains many reference variables which have to be determined through the minimization with three object functions. By choosing suitable weight function and its reference variables, the correct estimation of the depth of the brain sources can be achieved.
, 百拇医药
    Key words Magnetoencephalography Minimum norm Lp dipoles Weight function

    1 引 言

    脑磁图(Magnetoencephalography,MEG)是在脑电图EEG之后,80年代发展起来的又一研究大脑的临床应用设备。由脑神经冲动产生的微小的磁场被头表面超导磁量子(SQUID)[1]探头测量到。脑磁图的逆问题是根据测得的脑表面磁波寻求脑磁的电流源的位置和大小。该问题类似于根据多台天线测得的无线电信号确定电台的位置。对该问题的基本方法可分为两种,一种是多偶极子搜索法[1,2]。另一种是最小模法[3]。最小模法又分L2模法和L1模法。L2模法的解是连续分布的其解可表为广义逆的形式[4~6]。L1模法的解是稀少的[7~9],即只有个别几个网点上的源不为零。
, http://www.100md.com
    由于脑磁图问题的不定性,给出对脑磁源的分布模型是必要的。常见的一种模型假定脑磁源是稀少的,这正同上述L1最小模法的结果一致。但即使对L1最小模法仍有一个权重函数可自由选择。现有多种不同的方法[8,10]用来规定权重函数。不同的方法的解给出源的深度也不同。没有一种方法对不同深度的源都是适用的[11]

    我们在L1模法[2]的基础上发展起来,该方法三个目标函数,第一个目标函数是L1模。第二个目标函数为由偶极子源计算所得的磁场和测得的磁场的均方偏差。为了使这两个互相低逐的目标函数相结合,还得以某种条件控制这两个目标函数的相互比例,这由建立另一个条件目标函数来完成,即使在上述均方偏差不超过某一指标的条件下偶极子源的数目取最小。

    2 多目标函数最优化方法

    设探头各个通道测量到的磁场为H=[H1,H2,…Hi…,HM]t,t为转置。Hi为第i通道测得的沿线圈轴线方向磁场,M为探头的测量通道数目。将脑磁源所在空间用某种网格划分,例如用方格划分,每一网点布置相互垂直的三个电流偶极子。偶极子强度为X=[X1,X2,…Xj…XN]t,其中Xj为第j个偶极子源的强度,N为源的数目。每个源所在的位置在坐标xj=(xj,yj,zj)处。这些偶极子在探头测量通道处的磁场可由正问题,例如有限元法计算出来,这可线性的表示为
, http://www.100md.com
    BX+S=H (1)

    式中:B为N×M阶矩阵,它由正向问题计算。S=[S1,S2…Si,SM]t为噪音矢量;Si为第i通道的噪音。考虑一个权重因子W=diag([W1,W2…Wi…,WN]),记A=BW,Y=W-1X,此处Y=[Y1,Y2…Yi,YN]t。(1)式可改写为

    AY+S=H (2)

    定义权重函数矩阵为W=G或W=GWn,其中G随参数变化的权重函数,定义为G=diag(g(τ,x1),g(τ,x2),……g(τ,xN)),Wn为规一化的权重函数[B.Jeffs],包括这一规一化因子有时可使问题简化,其定义为Wn=diag(),式中‖b1‖,‖b2‖,……‖bN‖为矩阵A的列矢量的L2模。G为含J个参量τ=τ1,τ2,…τJ的空间权重函数,xj=(xj,yj,zj)为未知源的位置。
, 百拇医药
    M×N阶引导矩阵A可分解为A=UDVT,其中U,D,V为矩阵A的SVD(singular value decomposition)分解矩阵。设Ur,Dr,Vr为M×r,r×r,N×r矩阵,它们为矩阵U,D,V的截断矩阵,对应于矩阵A的前r个较大的SV(singular values)值。r在该方法中对应产生磁场的不为零的偶极子源的数目。利用上述定义式(2)可近似表示为

    VrTy=F (3)

    式中:F=Dr-1UrH,在式(3)的条件下对加权的未知源的模求最小 (4)

, http://www.100md.com     式中:minimize为最小;arc为最小时对应的自变量,对应的未知源为X*=WY*,由于所选空间权重函数g还含有未知参数τ,参量r也未知,X*=X*(r,τ)依赖于这些未知参量。下面给出确定这些参量的方法。定义相对计算磁场与测量磁场的相对误差函数E(r,τ)=,它也是参量r,τ的函数。‖.‖为L2模。首先对r给定一个初值r0,一般可选r0=,使上述相对误差函数取最小。 (5)

    这是参量τ的一个估计,记为τ*。在上述取最小的过程中对不同的参量τ值,我们得不断调用前一个取最小的公式(4)和(3)以得到X*
, http://www.100md.com
    对于固定的τ*值,我们可由下式找出r (6)

    并使其满足条件E(r,τ*)0。此处E0是相对误差函数的上限,事先可以给定,例如取E0=0.03,或E0=0.04等。在式(6)的取极植的过程中得调前面两个极值过程。计算表明参量r与τ基本是独立的,因此只须对固定的τ值计算一次r即可。得到r*后,可以由式(5)重新计算τ*,计算中用r*代替r0

    当得到了参数r*,τ*后,公式(3)和(4)给出对应的Y*=Y*(r**)值,进一步定出对应的X*=X*(r**)值。值得指出的是公式(3)和(4)可由FOCUSS算法实现[8]
, http://www.100md.com
    3 各种权重函数

    (1)一个权重参量方法[12]。该方法取权重函数矩阵的对角元素为 (7)

    非对角元素为零,式中R为头盔的最大半径。该方法适用于脑磁源主要分布在某一半径的球面附近,而该半径又未知的情况。

    (2)两个参量的方法。W=GWn,Wn由式(15)定出,而 (8)

    此处ρi=‖xi‖,ρλρτ为两个参数,R为头盔的最大半径,k取某一恒定值例如k=0.1。该方法适用于已知脑磁源分布在两个未知球面附近。
, 百拇医药
    (3)分区差值方法。W=GWn,将脑磁源分布的空间分成K个小的区域,每一区域中选出一点,该点函数g可取高值1,低值取0(或0.1等),其他点的值由这些点差值而得。高值代表该区域代表这个区域里含有源,低值区域正好相反。函数g有2K种变化,这些变化用参数λ表示,λ=1,2,……,2K。例如我们可以把半径分成K=10份,得到10个区域。

    (4)微扰方法。L1最小模方法往往不能准确的计算脑磁源的深度。最强的脑磁源是最重要的。为了准确计算最强的脑磁源的深度,可对它的深度位置进行微扰。

    G=diag{g(λ,X1),g(λ,X2)……g(λ,XJ)} (9) (10)
, http://www.100md.com
    其中:λ为其微扰参数;ρ0=‖xmax‖为前次计算最大源的位置xmax点到原点O的距离;ρ=‖x‖为其他源的位置x点到原点O的距离。

    4 讨论和总结

    上面仅仅叙述了多目标优化方法。但对其背景并没有详细交代。或许有人会问,为什么要引入这样三个目标?为什么要引入上面的含参数权重函数g(τ,x)?为什么要采用SVD过程?多目标法究竟比L1最小模方法有什么改进?

    首先我们知道L1最小模方法中权重矩阵W有许多不同的选择,其结果也强烈的依赖这些选择。这意味着权重矩阵包含着许多未定的因素。为了描述这些未定因素,我们在最小模方法的权重矩阵基础上引入了含参数权重函数g(τ,X)。但是这又产生了怎样确定权重函数的参数τ的问题。从偶极子搜索法的经验中我们知道如果所选偶极子源的数目太多,即使所选偶极子源远离搜索最佳位置时,计算磁场和测量磁场之间的相对误差仍然很小,这很不利偶极子源的搜索。所以偶极子源远离搜索最佳位置时适当减少偶极子源的数目,可引起上述相对误差的增加。较大的相对误差有利于偶极子源的搜索。L1最小模方法实际上也选了M个偶极子源,此处M是探头的数目,这些偶极子源只能在坐标网点位置上。恒定不变地选择源的数目与探头的数目一样是没有根据的,特别是当探头测量通道数目较大时,例如BTI公司生产头盔探头有通道148个之多。对L1最小模法,减少源的数量的方法是减少其约束方程组的数目。但我们不能通过减少探头数目来减少约束方程组的数目,因为这样丢失信息太多。不丢失太多信息,又减少约束方程组的数目的方法为SVD方法。在SVD方法中通过控制SV值的数目r,来控制约束方程组的数目,从而进一步控制源的数目。另一方面约束方程者(3)尽管其数目比式(2)减少了,仍然是约束方程(2)的很好近似,给出了对测量系统的较好描述。L1最小模法的相对误差为零。当适当的减少了源的数目后,我们发现该相对误差不再为零,较好的选择权重函数τ可使该相对误差尽可能的小。这样我们可用相对误差搜索参数τ,这正是式(5)。最后一个问题是究竟应该选多少个源好,即数目r的大小?因为L1最小模方法所选的源仅在坐标的网点上,如果实际源不在网点上,一个实际的源往往需要在实际源周围几个在网点上的源来描述。因此r肯定应大于实际源的数目。这仅是一个原则,实际问题中事先我们并不知道多少源存在。因此我们首先初选择r数目为M/2。公式(6)是控制源的数目尽可能少的情况下,又不致使在最佳参数τ时相对误差太大。该方法克服了L1方法中所求得的源的深度依赖权重函数的选择的问题。
, 百拇医药
    也许还有人要问,本文方法与偶极子源搜索法相比,实际上把搜索偶极子源改为搜索权重函数的参数,这样做又有什么优点呢?本文方法确实是偶极子源搜索法与L1最小模方法相结合的产物。但选择少量几个参数,例如仅一个就可给出对该MEG问题一个很好的描述[12]。采用偶极子方法即使最粗糙的模型选择3个偶极子,每个偶极子含5个非线性参量,三个位置参量两个方向参量,和一个线性参量,即偶极子的大小,总共18个参量。18个参量要全局优化是相当费时的,因此一般得采用梯度类算法,但这类算法往往收敛到局部极值处。最优化方法的参数可少得多。

    我们详细的推导了用于脑磁图定位问题的三目标最优化方法。比较了三目标最优化方法相对L1最小模方法的主要优越性。讨论了三目标最优化方法进一步发展的几种可能性。

    参考文献

    1 Matti S. Hamalainen.Magnetoencephalography-theory, instrumentation, and applications to noninvasive studies of the working human brain. Rev Modern Physics, 1993;65(2)∶413
, 百拇医药
    2 Mosher JC, Lewis PS, Leahy RM. Multiple dipole modeling and localization from spatio-temporal MEG Data. IEEE Trans Biomedical Eng, 1992; 39∶541

    3 Jeffs B, Leahy R, Singh M. An evaluation f methods for neuromagnetic image reconstruction. IEEE Trans BME, 1987; 34(9)∶713

    4 Matti S, Hamalainen, Risto J. Ilmoniemi. Interpreting measured magnetic fields of the brain: Estimates of current distributions. Finland:Helsinki University of Technology, Technical Report TKK-F-A559, 1984
, http://www.100md.com
    5 Savas J. Basic mathematical and electromagnetic concepts of the bioelectromagnetic inverse problem. Phys Med Biol, 1987; 32(1)∶11

    6 Wang JZ, Williamson SJ, Kaufman L. Magnetic Source Images Determined by Lead-Field Analysis: The Unique Minimum-Norm Least-Squares Estimation. IEEE Trans Bio Eng, 1992; 39(7)∶665

    7 Ioanides A, Bolton J, Clarke C. Continuous probabilistic solutions to the biomagnetic inverse problem. Inverse Problems, 1990; 6∶523
, 百拇医药
    8 Corodnitsky IF, George JS, Rao BD. Neuromagnetic source imaging with FOCUSS:a recursive weighted minimum norm algorithm. EEG and Clinical Neurophysiol, 1995; 95∶231

    9 Kanta Matsuura, Yoichi Okabe. Selective minimum-norm solution of the biomagnetic inverse problem. IEEE Transactions on Biomedical Engineering, 1995; 42(6)∶608

    10 Clarke CJS. Inverse Problems, 1989;(5)∶999

    11 Leahy RM, Mosher JC, Phillps JW. Proceedings of the 10th International Conference on Biomagnetism, BIOMA-G 96, Sanata Fe, New Mexico, February 1996
, 百拇医药
    12 Zhao Shuangren,Horst Halling. Minimum L1 Norm MEG Reconstruction Minimising Signal Deviation Using a Reduced Lead field, 18th Annual International Conference IEEE Engineering in Medicine and Biology Society. October 31-Novermber 3,1996, Amsterdam-The Neth-

    erlands, in CD-ROM

    (收稿:1998-07-17 修回:1998-11-23), 百拇医药