三目标最优化结合磁场信号SVD正交分解用于脑磁源的定位计算(MEG)
作者：赵双任 蒋大宗
单位：西安交通大学 (西安 710049)
关键词：脑磁图；最小模；偶极子；权重函数
北京生物医学工程990206 摘 要 脑磁源的定位是脑磁图(magnetoencephalography简称MEG)研究的一个基本问题。该问题属不定问题，即根据探头测量的微小磁场所建立的方程组求得的未知脑磁源有无穷组解。因此需要通过建立合理的数学模型补充适当信息。脑磁场测量信号SVD正交分解且进行单偶极子搜索方法(music)已被广泛采用。该方法假定脑磁源是许多孤立的偶极子源，且偶极子源之间在统计上是互相独立的。若脑磁源不满足上述假定，一般讲该方法失效。三目标最优化法是最近推出的一方法，若偶极子源分布在脑空间内任意深度平行于探头表面的一层附近，该方法可给出正确的解。三目标最优化法与脑磁场测量信号的SVD正交分解方法相结合可得一新方法，它将脑空间划分成几个子区域，每一区域在平行探头表面不同深度的一层附近。新方法仅要求每一区域内的脑磁源与其它区域内的源在统计上互相独立。而区域内的源是允许互相相关的。与前两种方法相比，新方法进一步放松了对脑磁源的要求，因此有更广泛的适用性。本文详细推导并论证该新方法并讨论它对磁场的源的独立性的要求。本文所述新方法也适用对脑电图EEG源的定位计算。
, 百拇医药
Locating the Brain Sources by Applying the Combination the
Method of Minimization with Three Object Functions and the
Method of Singular Value Decomposition to the Signals of
Brain Magnetic Fields
Zhao Shuangren, Jiang Dazong
(Xi′an Jiaotong University, Xi′an, 710049)
Abstract
, 百拇医药
Locating the sources of brain magnetic fields is a basic problem of magnetoencephalography (abbreviated as MEG). The problem is an “indefinite” one, that is: infinite solutions can be obtained from the equations according to the small magnetic fields measured by detectors. Therefore, a mathematical model is necessary to give more information. The method searching with single dipole applying the technique of the singular value decomposition (SVD) to the signals of the measured magnetic fields (MUSIC), has been widely used. It is supposed that the sources of brain magnetic fields are many compact electric current dipoles, which are independent with each other in respect of statistics, and all of them stay far away with each other. If these conditions are not satistied, this method fails in general. The method, using minimizations with three object functions to locate the sources of brain magnetic fields was developed recently. If all sources of brain magnetic fields stay in the neighborhood of one layer, which is underneath the detectors in any depth, correct results can be achieved with this method. In this paper a new method is developed through the combination of the method, using minimizations with three object functions and the method of the decomposition of the measured brain magnetic fields, using SVD. The new method divides the brain space into many subspaces. Each subspace contains the neighborhood of one layer with fixed depth underneath the surface of the detectors. The requirement to the new method is that the brain sources in one subspace is independent with the brain sources in other subspaces in tespect of statistics. The new method further decreases the requirements for the brain sources, so its applicability is increased. This paper careful derivation of the method, proof of the results and discussion of the requirements of independence for the brain sources. This new method is also suitable for EEG.
, 百拇医药
Key words:Magnetoencephalography (MEG); Minimum norm; LP; Dipoles; Weight function
0 引 言
脑磁图(magnetoencephalography简称MEG)是在脑电图EEG之后，80年代发展起来的又一大脑研究和临床应用设备。由脑神经冲动产生的微小的磁场被头表面超导磁量子干涉仪(squid)[1]探头测量到。基于这样测量来重建神经源冲动的位置和数值大小是现代科学的又一挑战。
脑磁图MEG要解决的问题包括：滤波、差值、正计算问题及源的定位问题。滤波是在时域或频域进一步提取必要的数据，滤去噪声[2]。数据差值主要在空域由差值显示在脑皮层外表面探头测得的数据或其梯度或其Laplace算子的强弱[3]。正计算问题有球模型的解析方法[4]，边界元法和有限元法[5]等数值方法。
, 百拇医药
脑磁图的源定位问题是根据测得的脑表面磁场寻求产生脑磁场的电流源的位置和大小。该问题类似于根据多台天线测得的无线电信号确定电台的位置。在此我们简述现有的定位方法。第一种是多偶极子搜索法[1]，该方法以脑磁源电流偶极子计算所得的磁场和测得的磁场的均方误差为目标函数，调整偶极子的位置和极矩使该均方误差取最小。该法又分全局搜索和梯度类搜索法。全局搜索当所取的偶极子数目超过3个时计算量往往太大，难以得到计算结果。梯度搜索法一般只收敛到局部极值处。第二种方法是脑磁场测量信号SVD正交分解且进行单偶极子搜索[6]。该方法假定脑磁源是许多孤立的偶极子源，且偶极子源之间在统计上是互相独立的，若脑磁源不满足上述假定，一般讲该方法失效。第三种是最小模法[7]。最小模法在脑磁源存在的空间里按一定规律布置网格，每个网格的网点上固定两个或三个小的互相垂直电流偶极子源，这些小电流源恰好产生测得的磁场。满足上述条件的源分布有无数多种选择，其中在某种意义下的模取最小的源的分布是该方法的解。根据模的定义又分L2模法和L1模法。L2模法的解是连续分布的[4，8，9]，即每个网格上都有一源存在，该方法的解与探头上各个通道的测量数据间的关系是线性的，可用广义逆的形式表示。L1模法的解是稀疏的[10～12]，即只有个别几个网点上的源不为零。当模界于L2和L1模之间，其解也界于它们的解之间。其解与探头上各个通道的测量数据间的关系是非线性的，一般得用迭代方法求解[11]。当模界于L1和L0之间时，其结果与L1模情况接近，仅在迭代的收敛速度和收敛性方面有所变化[11]。各种最小模法在模的定义中有一权重因子由不同的归一化方法确定[2，11，13]。不同的归一方法对应的解在源的深度方面是不同的。没有一种权重因子对各种不同深度的源都适用[14]。第四种是在L1模方法基础上发展起来的三目标优化方法[15]。在确定权重函数方面与上述各种规一化方法不同。该方法有三个目标函数，权重函数中含一组待定的参量，该参量对应脑磁源的深度。一个目标函数是L1模；一个为偶极子源计算所得的磁场和测得的磁场的均方误差；还有一个是偶极子源的数目。该方法要求偶极子源分布在脑空间内任意深度平行探头表面的一层附近，其深度可任意给定。改变其参量可在脑空间不同深度搜索脑磁源。若脑磁源分布在不同深度，一般地该方法失效。以上四种方法是在不同的要求或前提下工作因此各有局限性。寻找进一步放松要求的方法是脑磁源定位的关键。
, 百拇医药
本文将第二种方法中脑磁场测量信号SVD正交分解法同第四种方法三目标优化方法有机地结合起来构成一新方法。它将脑空间划分成几个子区域，每一区域在平行探头表面的一层附近。新方法仅要求每一区域内的脑磁源与其它区域内的源在统计上是互相独立的。而区域内的源允许是相关的。与上述第二种方法相比，新方法不再要求所有的源互相正交。与第四种方法相比，新方法允许脑磁源分布在不同的深度层内，只要各层的源与其它层的源是不相关的。新方法进一步放松了对脑磁源的要求，因此有更广泛的适用性。
1 新方法
设探头各个通道测量到的磁场为H=[H1，H2…Hi…，HM]T。Hi为第i通道测得的沿线圈轴线方向磁场，M为探头的测量通道数目。考虑测得的磁场也是时间的函数，设磁场沿时间有T个采样点。Hi为T个元素的矢量，H为M×T阶矩阵。将脑磁源所在空间用某种网格划分，例如用方格划分，每一网点布置相互垂直的三个电流偶极子。偶极子强度为X=[X1，X2…Xj…XN]T，其中N是源的数目，Xj是第j个偶极子源的强度，则它也随时间变化，设时间抽样有T个点，它是T个元素的矢量。X为N×T阶矩阵。每个源Xj所在的位置在坐标xj=(xj,yj,zj)处。这些偶极子在探头测量通道线圈处的磁场可由正问题方法，例如有限元法[5]计算出来，脑磁源和测量的磁场间关系可线性地表示为：
, 百拇医药
(1)
此处B为M×N阶矩阵，称为引导场矩阵，该矩阵的矩阵元代表每个网点某一方向上的单位偶极子对探头某个通道测量值的贡献，它由正向问题计算。S=[S1，S2…Si…SM]T为噪声矩阵，Si为第i通道的噪音，为T个元素的矢量。S与H都是M×T阶矩阵。
考虑磁场矩阵H可用SVD方法作正交分解[6]
(2)
UH为M×M阶正交矩阵，VH为T×T阶正交矩阵。DH为对角矩阵，DH=diag(δ1，δ2，……δM)，δi为第i大奇异值。取D1=diag(δ1，0，……0)，D2=diag(0,δ2，……0)，……，DM=diag(0,0,……δM)。故有DH=D1+D2+……+DM，因此，(2)式又可表示为
, http://www.100md.com
(3)
此处Hi=UHDiVTH，Hi。而Hj,i≠j可看成由统计上不相关的源产生的磁场。
考虑X=WY,W=diag(W1,W2,……WN)是权重矩阵[13]，Y是加权后的未知偶极子源的大小。式(1)可化为
(4)
上式只考虑磁场的一个独立磁场分量。S’为同S类似的噪声矩阵。考虑矩阵A=BW可用SVD方法正交分解成
, 百拇医药
(5)
U为M×M阶正交矩阵，V为N×N阶正交矩阵。D为对角矩阵，D=diag(ε1,ε2,……εM)。ε1，ε2，……εM为按大小排列的奇异值。取前r个奇异值不为零，其他的均为零可得
(6)
Ar称为矩阵A的前r个奇异值的近似，rr为U的前r列子矩阵，它是M×r阶正交矩阵，Vr为V的前r列子矩阵，它为N×r阶正交矩阵，D=diag(ε1，ε2，……εr)。考虑式(6)，式(4)可改写为
, http://www.100md.com
(7)
式中，
。式(7)是含r个方程的方程组，而式(4)是M个方程的方程组。因此式(7)比式(4)规模小，更利于计算。式(4)中权重函数的一个例子[13，15]为
(8)
式中R是头盔的最大半径，τ为待定参数。xi为未知源Xi的位置矢量。满足式(7)的最小模解由
(9)
且满足方程组
, http://www.100md.com
(10)
得到。传统的最小模法采用的约束条件不是(10)式而是AY=Hi，为了区别起见，我们称式(9)-(10)所提供的方法为改进的最小模法。
定义为
(11)
得到加权未知源Y*后可进一步得到未知源本身
(12)
由式(8)-(10)知未知源同待定参数τ,r相关即
, 百拇医药
(13)
当p=1时，式(9)-(10)也可用线性规划的理论求解[12]。由线性规划的理论知，上述X*中不为零的分量不超过r。因此r也可被看成未知源的数量。对于p<1的情况也与p=1时的情况一致，r可看成未知源的数量。p≤1时，传统的最小模方法未知源的数量为M，即探头通道数。现在使用的MEG设备例如BTI公司的头盔探头含148个通道。取脑磁源的数目也多达148个，是没有根据的。一个偶极子源得用它周围几个网点上的偶极子模拟，即使如此取脑磁源大数目同探头数目一致也确实有点过于多了些。改进了的最小模方法式(9)-(10)克服了这一困难，我们可取未知源的数目为r个，r是小于探头通道的任意数。定义相对计算磁场与测量磁场的某一正交分量的误差函数
(14)
, http://www.100md.com 它也是参量τ,r的函数。
为L2模。r=M时A≡Ar，式(10)与BX-Hi=0等效，故有Ei(r,τ)=0。模拟计算表明误差函数Ei(r,τ)随未知源的数目r减小而增大。另外经仿真计算我们发现，若未知源在距离探头表面某一深度层附近，当我们选择参数τ使误差函数Ei(r，τ)取最小的同时，式(9)-(10)的解X*(τ,r)给出正确的未知源的深度。由此导出下面确定这些参数τ,r的方法。
首先对r给定一个初值r0，一般可选
，这意味着我们初选未知源的数目为探头通道数目的一半。进一步使上述相对误差函数取最小
, 百拇医药
(15)
由公式(15)可得上述参量τ的一个估计，记为τ*。在上述取最小的过程中对不同的参量τ值，我们得不断调用前一个取最小的步骤即公式(9)-(11)以得到X*，并用公式(14)计算Ei(r,τ)。式(15)只有一个参数τ变化故可采用遍历所有值的搜索，即全局搜索。
对于固定的τ*值，我们可由下式找出r
(16)
并使其满足条件
(17)
, http://www.100md.com
此处γ0是相对误差函数的上限，事先可以给定，例如取γ0=0.03，或γ0=0.04等。在式(16)的取极值的过程中需要调用前面两个极值过程。只须对固定的τ值计算一次r即可。得到r*后由公式(15)重新计算τ*，计算中用r*代替r0。式(16)(17)表明选取尽量少的源但又不至相对误差超出其上限γ0。
当我们得到了参数r*,τ*后，公式(9)-(12)给出对应的X*=X*(r*,τ*)值。式(9)-(17)中含三个取最小值的过程，分别用来确定未知量X，r,τ。模拟计算表明这些未知量，在取最小值的意义下互相独立，因此经式(9)-(17)一次循环即可将X*,r*,τ*定出。
, http://www.100md.com
上述求解是对测量磁场的一个正交分量Hi进行的。我们还需对其他所有具有较大奇异值的分量做相同的工作。
值得指出的是公式(9)-(10)的最小模法可由类似FOCUSS[12]算法的迭代法实现。其迭代步骤为
(18)
(19)
(20)
式中1为N×N阶单位矩阵，符号“+”为广义逆，
是独立的源产生的磁场，它是2维矩阵，其两个指标一个是探头通道数，一个是时间采样点数。作上述迭代时须对固定的时间一个一个进行，对固定的时间，Hi只是一个M元素的矢量。元素设Y(n)收敛到Y*则由式(12)给出未知元的解X*。当q=0.5对应L1模的解。上式中广义逆可由下式[12]计算
, 百拇医药
(21)
此处G为任意矩阵，η为光滑因子，为接近零的常数。
2 分析和总结
测量磁场的SVD正交分解是首先用于单偶极子搜索。该方法要求偶极子源是互相独立的。此处两源Xt,Yt独立是指
(22)
角标t表示时间。这样各个偶极子源的磁场可通过正交分解SVD算法分解开来。分解开来的磁场只是一个偶极子产生的场，因此可采用单偶极子搜索方法找到解。倘若两个以上偶极子源是相关的，式(22)便不满足。正交分解SVD算法不能将它们的磁场分开。我们知道用单偶极子在两个以上偶极子产生的磁场中搜索往往得不到正确的解。因此，此时该法失效。
, 百拇医药
上述新算法是将三目标优化方法，与测量磁场的SVD正交分解结合起来。三目标优化方法也是一种搜索方法，不过不是一个点一个点的搜索。该方法是对不同的参数τ搜索。由式(8)知，参数τ对应于源所在的深度。因此该方法实际上是对不同深度的一层源进行搜索。这里说的一个深度层实际包括了一个深度和该深度附近的一个区域。只要每层源作为一个整体与其它层的源是互相独立的，它们的磁场就可用SVD正交分解分开。因此由三目标优化方法的层搜索可得到正确的解。该方法允许一个深度层内的源是相关的。尚若相关的源在不同的深度存在，一般的该方法失效。尽管如此该新算法对源独立性的要求比起磁场测量信号SVD正交分解且进行单偶极子搜索的方法还是大大放宽了。
本文将三目标脑磁源定位法SVD正交分解方法相结合，推出了一种新的方法。该方法用深度层区域搜索取代单偶极子搜索从而可放宽对脑磁源相互独立的要求。本文详细论述了这一原理。
作者简介：赵双任：男，45岁，1976年毕业于北京工业学院精密机械与电子专业。1987年西北电讯工程学院电磁场与微波技术专业研究生毕业。自1990年起德国Juelich研究中心客座学者。1997年起被西安交大接收为论文博士。主要工作范围：脑磁图(MEG)计算与研究，断层扫描(CT、SPECT、PET、MRT)算法，电磁场微波计算。
, http://www.100md.com
参考文献
1 Hamalainen M S, et al. Magnetoencephalography-theory, instrumentation, and applications to noninvasive studies of the working human brain.Rev Modern Physics,1993，65(2):413
2 Bertrand O, Tallon-Baudry C, Pernier J. Time-frequency analysis of oscillatory gamma-band activity: wavelet approach and phase-locaking estimation. Biomagnetism 1996, Santa Fe, Biomag Conference, 16-21 Feb 96
3 Zhao S R, Heer H, Loanniedes AA, Wagener M, et al. Interpolation of magnetic fields and their gradients for MEG data with 3D spline function. Bild verarbeitung fuer medizin algorithmen Systeme, Anwendungen, Novermber 8, Proceedings des Aachener Workshops am 8. und 9. November 1996, p 297
, 百拇医药
4 Savas J. Basic mathematical and electromagnetic concepts of the bioelectromagnetic inverse problem. Phys Med Biol, 1987, 32(1):11
5 Jens Haueisen, Ceon Ramon, Piotr Czapaki, et al. On the influence of volume currents and extended sources on neuromagnetic fields: a simulation study. Submitted to the Annals of Biomedical Engineering second revisionsubmitted April 25, 1995.
6 Mosher JC, Lewis PS, Leahy RM. Multiple dipole modeling and localization from spatio-temporal MEG data. IEEE Trans Biomedical Eng, 1992, 39:541
, 百拇医药
7 Jeffs B, Leahy R, and Singh M. An evaluation of methods for neuromagnetic image reconstruction.IEEE Trans BME, 1987, 34(9):713
8 Hamalainen MS and Ilmoniemi R J. Interpreting measured magnetic fields of the brain: estimates of current distributions. Helsinki University of Technology. Finland, Technical Report TKK-F-A559, 1984
9 Wang J Z, Williamson S J, and 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
, http://www.100md.com
10 Loanides A, Bolton J, and Clarke C. Continuous probabilistic solutions to the biomagnetic inverse problem. Inverse Problems, 1990, 6:523
11 Corodnitsky I F, George J S, and Rao B D. Neuromagnetic source imaging with FOCUSS: a recursive weighted minimum norm algorithm. EEG and Clinical Neurophysiol, 1995, 95:231
12 Kanta Matsuura and Yoichi okabe. Selective minimum-norm solution of the biomagnetic inverse problem. IEEE Transactions on Biomedical Engineering, 1995, 42(6):608
, 百拇医药
13 Clarke C J S.Inverse Problems,1989, 5:999. Neuroscience and Tech., Lyon, France, Nov., 1992. on Biomedical Engineering, 1995, 42(6):60
14 Leahy R M, Mosher J C and Phillps J W. Proceedings of the 10th International Conference on Biomagnetism, BIOMAG 96, Sanata Fe, New Mexico, February 1996
15 Zhao S R, Horst Halling. Minimum L1 norm MEG reconstruction minimising signal deviation using a reduced lead field. Proceedings of the 18th Annual International Conference IEEE Engineering in Medicine and Biology Society, Part 2(of 5). (or Annual Interational Conference of the IEEE Engineering in Medicine and Biology Proceedings V 2 1996. IEEE, Piscataway, NJ, USA, 96CB36036. p 814-816)
(1998-06-16收稿，1998-12-03修回), http://www.100md.com
单位：西安交通大学 (西安 710049)
关键词：脑磁图；最小模；偶极子；权重函数
北京生物医学工程990206 摘 要 脑磁源的定位是脑磁图(magnetoencephalography简称MEG)研究的一个基本问题。该问题属不定问题，即根据探头测量的微小磁场所建立的方程组求得的未知脑磁源有无穷组解。因此需要通过建立合理的数学模型补充适当信息。脑磁场测量信号SVD正交分解且进行单偶极子搜索方法(music)已被广泛采用。该方法假定脑磁源是许多孤立的偶极子源，且偶极子源之间在统计上是互相独立的。若脑磁源不满足上述假定，一般讲该方法失效。三目标最优化法是最近推出的一方法，若偶极子源分布在脑空间内任意深度平行于探头表面的一层附近，该方法可给出正确的解。三目标最优化法与脑磁场测量信号的SVD正交分解方法相结合可得一新方法，它将脑空间划分成几个子区域，每一区域在平行探头表面不同深度的一层附近。新方法仅要求每一区域内的脑磁源与其它区域内的源在统计上互相独立。而区域内的源是允许互相相关的。与前两种方法相比，新方法进一步放松了对脑磁源的要求，因此有更广泛的适用性。本文详细推导并论证该新方法并讨论它对磁场的源的独立性的要求。本文所述新方法也适用对脑电图EEG源的定位计算。
, 百拇医药
Locating the Brain Sources by Applying the Combination the
Method of Minimization with Three Object Functions and the
Method of Singular Value Decomposition to the Signals of
Brain Magnetic Fields
Zhao Shuangren, Jiang Dazong
(Xi′an Jiaotong University, Xi′an, 710049)
Abstract
, 百拇医药
Locating the sources of brain magnetic fields is a basic problem of magnetoencephalography (abbreviated as MEG). The problem is an “indefinite” one, that is: infinite solutions can be obtained from the equations according to the small magnetic fields measured by detectors. Therefore, a mathematical model is necessary to give more information. The method searching with single dipole applying the technique of the singular value decomposition (SVD) to the signals of the measured magnetic fields (MUSIC), has been widely used. It is supposed that the sources of brain magnetic fields are many compact electric current dipoles, which are independent with each other in respect of statistics, and all of them stay far away with each other. If these conditions are not satistied, this method fails in general. The method, using minimizations with three object functions to locate the sources of brain magnetic fields was developed recently. If all sources of brain magnetic fields stay in the neighborhood of one layer, which is underneath the detectors in any depth, correct results can be achieved with this method. In this paper a new method is developed through the combination of the method, using minimizations with three object functions and the method of the decomposition of the measured brain magnetic fields, using SVD. The new method divides the brain space into many subspaces. Each subspace contains the neighborhood of one layer with fixed depth underneath the surface of the detectors. The requirement to the new method is that the brain sources in one subspace is independent with the brain sources in other subspaces in tespect of statistics. The new method further decreases the requirements for the brain sources, so its applicability is increased. This paper careful derivation of the method, proof of the results and discussion of the requirements of independence for the brain sources. This new method is also suitable for EEG.
, 百拇医药
Key words:Magnetoencephalography (MEG); Minimum norm; LP; Dipoles; Weight function
0 引 言
脑磁图(magnetoencephalography简称MEG)是在脑电图EEG之后，80年代发展起来的又一大脑研究和临床应用设备。由脑神经冲动产生的微小的磁场被头表面超导磁量子干涉仪(squid)[1]探头测量到。基于这样测量来重建神经源冲动的位置和数值大小是现代科学的又一挑战。
脑磁图MEG要解决的问题包括：滤波、差值、正计算问题及源的定位问题。滤波是在时域或频域进一步提取必要的数据，滤去噪声[2]。数据差值主要在空域由差值显示在脑皮层外表面探头测得的数据或其梯度或其Laplace算子的强弱[3]。正计算问题有球模型的解析方法[4]，边界元法和有限元法[5]等数值方法。
, 百拇医药
脑磁图的源定位问题是根据测得的脑表面磁场寻求产生脑磁场的电流源的位置和大小。该问题类似于根据多台天线测得的无线电信号确定电台的位置。在此我们简述现有的定位方法。第一种是多偶极子搜索法[1]，该方法以脑磁源电流偶极子计算所得的磁场和测得的磁场的均方误差为目标函数，调整偶极子的位置和极矩使该均方误差取最小。该法又分全局搜索和梯度类搜索法。全局搜索当所取的偶极子数目超过3个时计算量往往太大，难以得到计算结果。梯度搜索法一般只收敛到局部极值处。第二种方法是脑磁场测量信号SVD正交分解且进行单偶极子搜索[6]。该方法假定脑磁源是许多孤立的偶极子源，且偶极子源之间在统计上是互相独立的，若脑磁源不满足上述假定，一般讲该方法失效。第三种是最小模法[7]。最小模法在脑磁源存在的空间里按一定规律布置网格，每个网格的网点上固定两个或三个小的互相垂直电流偶极子源，这些小电流源恰好产生测得的磁场。满足上述条件的源分布有无数多种选择，其中在某种意义下的模取最小的源的分布是该方法的解。根据模的定义又分L2模法和L1模法。L2模法的解是连续分布的[4，8，9]，即每个网格上都有一源存在，该方法的解与探头上各个通道的测量数据间的关系是线性的，可用广义逆的形式表示。L1模法的解是稀疏的[10～12]，即只有个别几个网点上的源不为零。当模界于L2和L1模之间，其解也界于它们的解之间。其解与探头上各个通道的测量数据间的关系是非线性的，一般得用迭代方法求解[11]。当模界于L1和L0之间时，其结果与L1模情况接近，仅在迭代的收敛速度和收敛性方面有所变化[11]。各种最小模法在模的定义中有一权重因子由不同的归一化方法确定[2，11，13]。不同的归一方法对应的解在源的深度方面是不同的。没有一种权重因子对各种不同深度的源都适用[14]。第四种是在L1模方法基础上发展起来的三目标优化方法[15]。在确定权重函数方面与上述各种规一化方法不同。该方法有三个目标函数，权重函数中含一组待定的参量，该参量对应脑磁源的深度。一个目标函数是L1模；一个为偶极子源计算所得的磁场和测得的磁场的均方误差；还有一个是偶极子源的数目。该方法要求偶极子源分布在脑空间内任意深度平行探头表面的一层附近，其深度可任意给定。改变其参量可在脑空间不同深度搜索脑磁源。若脑磁源分布在不同深度，一般地该方法失效。以上四种方法是在不同的要求或前提下工作因此各有局限性。寻找进一步放松要求的方法是脑磁源定位的关键。
, 百拇医药
本文将第二种方法中脑磁场测量信号SVD正交分解法同第四种方法三目标优化方法有机地结合起来构成一新方法。它将脑空间划分成几个子区域，每一区域在平行探头表面的一层附近。新方法仅要求每一区域内的脑磁源与其它区域内的源在统计上是互相独立的。而区域内的源允许是相关的。与上述第二种方法相比，新方法不再要求所有的源互相正交。与第四种方法相比，新方法允许脑磁源分布在不同的深度层内，只要各层的源与其它层的源是不相关的。新方法进一步放松了对脑磁源的要求，因此有更广泛的适用性。
1 新方法
设探头各个通道测量到的磁场为H=[H1，H2…Hi…，HM]T。Hi为第i通道测得的沿线圈轴线方向磁场，M为探头的测量通道数目。考虑测得的磁场也是时间的函数，设磁场沿时间有T个采样点。Hi为T个元素的矢量，H为M×T阶矩阵。将脑磁源所在空间用某种网格划分，例如用方格划分，每一网点布置相互垂直的三个电流偶极子。偶极子强度为X=[X1，X2…Xj…XN]T，其中N是源的数目，Xj是第j个偶极子源的强度，则它也随时间变化，设时间抽样有T个点，它是T个元素的矢量。X为N×T阶矩阵。每个源Xj所在的位置在坐标xj=(xj,yj,zj)处。这些偶极子在探头测量通道线圈处的磁场可由正问题方法，例如有限元法[5]计算出来，脑磁源和测量的磁场间关系可线性地表示为：
, 百拇医药
(1)
此处B为M×N阶矩阵，称为引导场矩阵，该矩阵的矩阵元代表每个网点某一方向上的单位偶极子对探头某个通道测量值的贡献，它由正向问题计算。S=[S1，S2…Si…SM]T为噪声矩阵，Si为第i通道的噪音，为T个元素的矢量。S与H都是M×T阶矩阵。
考虑磁场矩阵H可用SVD方法作正交分解[6]
(2)
UH为M×M阶正交矩阵，VH为T×T阶正交矩阵。DH为对角矩阵，DH=diag(δ1，δ2，……δM)，δi为第i大奇异值。取D1=diag(δ1，0，……0)，D2=diag(0,δ2，……0)，……，DM=diag(0,0,……δM)。故有DH=D1+D2+……+DM，因此，(2)式又可表示为
, http://www.100md.com
(3)
此处Hi=UHDiVTH，Hi。而Hj,i≠j可看成由统计上不相关的源产生的磁场。
考虑X=WY,W=diag(W1,W2,……WN)是权重矩阵[13]，Y是加权后的未知偶极子源的大小。式(1)可化为
(4)
上式只考虑磁场的一个独立磁场分量。S’为同S类似的噪声矩阵。考虑矩阵A=BW可用SVD方法正交分解成
, 百拇医药
(5)
U为M×M阶正交矩阵，V为N×N阶正交矩阵。D为对角矩阵，D=diag(ε1,ε2,……εM)。ε1，ε2，……εM为按大小排列的奇异值。取前r个奇异值不为零，其他的均为零可得
(6)
Ar称为矩阵A的前r个奇异值的近似，r
, http://www.100md.com
(7)
式中，
。式(7)是含r个方程的方程组，而式(4)是M个方程的方程组。因此式(7)比式(4)规模小，更利于计算。式(4)中权重函数的一个例子[13，15]为(8)
式中R是头盔的最大半径，τ为待定参数。xi为未知源Xi的位置矢量。满足式(7)的最小模解由
(9)
且满足方程组
, http://www.100md.com
(10)
得到。传统的最小模法采用的约束条件不是(10)式而是AY=Hi，为了区别起见，我们称式(9)-(10)所提供的方法为改进的最小模法。
定义为(11)
得到加权未知源Y*后可进一步得到未知源本身
(12)
由式(8)-(10)知未知源同待定参数τ,r相关即
, 百拇医药
(13)
当p=1时，式(9)-(10)也可用线性规划的理论求解[12]。由线性规划的理论知，上述X*中不为零的分量不超过r。因此r也可被看成未知源的数量。对于p<1的情况也与p=1时的情况一致，r可看成未知源的数量。p≤1时，传统的最小模方法未知源的数量为M，即探头通道数。现在使用的MEG设备例如BTI公司的头盔探头含148个通道。取脑磁源的数目也多达148个，是没有根据的。一个偶极子源得用它周围几个网点上的偶极子模拟，即使如此取脑磁源大数目同探头数目一致也确实有点过于多了些。改进了的最小模方法式(9)-(10)克服了这一困难，我们可取未知源的数目为r个，r是小于探头通道的任意数。定义相对计算磁场与测量磁场的某一正交分量的误差函数
(14)
, http://www.100md.com 它也是参量τ,r的函数。
为L2模。r=M时A≡Ar，式(10)与BX-Hi=0等效，故有Ei(r,τ)=0。模拟计算表明误差函数Ei(r,τ)随未知源的数目r减小而增大。另外经仿真计算我们发现，若未知源在距离探头表面某一深度层附近，当我们选择参数τ使误差函数Ei(r，τ)取最小的同时，式(9)-(10)的解X*(τ,r)给出正确的未知源的深度。由此导出下面确定这些参数τ,r的方法。首先对r给定一个初值r0，一般可选
，这意味着我们初选未知源的数目为探头通道数目的一半。进一步使上述相对误差函数取最小, 百拇医药
(15)
由公式(15)可得上述参量τ的一个估计，记为τ*。在上述取最小的过程中对不同的参量τ值，我们得不断调用前一个取最小的步骤即公式(9)-(11)以得到X*，并用公式(14)计算Ei(r,τ)。式(15)只有一个参数τ变化故可采用遍历所有值的搜索，即全局搜索。
对于固定的τ*值，我们可由下式找出r
(16)
并使其满足条件
(17)
, http://www.100md.com
此处γ0是相对误差函数的上限，事先可以给定，例如取γ0=0.03，或γ0=0.04等。在式(16)的取极值的过程中需要调用前面两个极值过程。只须对固定的τ值计算一次r即可。得到r*后由公式(15)重新计算τ*，计算中用r*代替r0。式(16)(17)表明选取尽量少的源但又不至相对误差超出其上限γ0。
当我们得到了参数r*,τ*后，公式(9)-(12)给出对应的X*=X*(r*,τ*)值。式(9)-(17)中含三个取最小值的过程，分别用来确定未知量X，r,τ。模拟计算表明这些未知量，在取最小值的意义下互相独立，因此经式(9)-(17)一次循环即可将X*,r*,τ*定出。
, http://www.100md.com
上述求解是对测量磁场的一个正交分量Hi进行的。我们还需对其他所有具有较大奇异值的分量做相同的工作。
值得指出的是公式(9)-(10)的最小模法可由类似FOCUSS[12]算法的迭代法实现。其迭代步骤为
(18)
(19)
(20)
式中1为N×N阶单位矩阵，符号“+”为广义逆，
是独立的源产生的磁场，它是2维矩阵，其两个指标一个是探头通道数，一个是时间采样点数。作上述迭代时须对固定的时间一个一个进行，对固定的时间，Hi只是一个M元素的矢量。元素设Y(n)收敛到Y*则由式(12)给出未知元的解X*。当q=0.5对应L1模的解。上式中广义逆可由下式[12]计算, 百拇医药
(21)
此处G为任意矩阵，η为光滑因子，为接近零的常数。
2 分析和总结
测量磁场的SVD正交分解是首先用于单偶极子搜索。该方法要求偶极子源是互相独立的。此处两源Xt,Yt独立是指
(22)
角标t表示时间。这样各个偶极子源的磁场可通过正交分解SVD算法分解开来。分解开来的磁场只是一个偶极子产生的场，因此可采用单偶极子搜索方法找到解。倘若两个以上偶极子源是相关的，式(22)便不满足。正交分解SVD算法不能将它们的磁场分开。我们知道用单偶极子在两个以上偶极子产生的磁场中搜索往往得不到正确的解。因此，此时该法失效。
, 百拇医药
上述新算法是将三目标优化方法，与测量磁场的SVD正交分解结合起来。三目标优化方法也是一种搜索方法，不过不是一个点一个点的搜索。该方法是对不同的参数τ搜索。由式(8)知，参数τ对应于源所在的深度。因此该方法实际上是对不同深度的一层源进行搜索。这里说的一个深度层实际包括了一个深度和该深度附近的一个区域。只要每层源作为一个整体与其它层的源是互相独立的，它们的磁场就可用SVD正交分解分开。因此由三目标优化方法的层搜索可得到正确的解。该方法允许一个深度层内的源是相关的。尚若相关的源在不同的深度存在，一般的该方法失效。尽管如此该新算法对源独立性的要求比起磁场测量信号SVD正交分解且进行单偶极子搜索的方法还是大大放宽了。
本文将三目标脑磁源定位法SVD正交分解方法相结合，推出了一种新的方法。该方法用深度层区域搜索取代单偶极子搜索从而可放宽对脑磁源相互独立的要求。本文详细论述了这一原理。
作者简介：赵双任：男，45岁，1976年毕业于北京工业学院精密机械与电子专业。1987年西北电讯工程学院电磁场与微波技术专业研究生毕业。自1990年起德国Juelich研究中心客座学者。1997年起被西安交大接收为论文博士。主要工作范围：脑磁图(MEG)计算与研究，断层扫描(CT、SPECT、PET、MRT)算法，电磁场微波计算。
, http://www.100md.com
参考文献
1 Hamalainen M S, et al. Magnetoencephalography-theory, instrumentation, and applications to noninvasive studies of the working human brain.Rev Modern Physics,1993，65(2):413
2 Bertrand O, Tallon-Baudry C, Pernier J. Time-frequency analysis of oscillatory gamma-band activity: wavelet approach and phase-locaking estimation. Biomagnetism 1996, Santa Fe, Biomag Conference, 16-21 Feb 96
3 Zhao S R, Heer H, Loanniedes AA, Wagener M, et al. Interpolation of magnetic fields and their gradients for MEG data with 3D spline function. Bild verarbeitung fuer medizin algorithmen Systeme, Anwendungen, Novermber 8, Proceedings des Aachener Workshops am 8. und 9. November 1996, p 297
, 百拇医药
4 Savas J. Basic mathematical and electromagnetic concepts of the bioelectromagnetic inverse problem. Phys Med Biol, 1987, 32(1):11
5 Jens Haueisen, Ceon Ramon, Piotr Czapaki, et al. On the influence of volume currents and extended sources on neuromagnetic fields: a simulation study. Submitted to the Annals of Biomedical Engineering second revisionsubmitted April 25, 1995.
6 Mosher JC, Lewis PS, Leahy RM. Multiple dipole modeling and localization from spatio-temporal MEG data. IEEE Trans Biomedical Eng, 1992, 39:541
, 百拇医药
7 Jeffs B, Leahy R, and Singh M. An evaluation of methods for neuromagnetic image reconstruction.IEEE Trans BME, 1987, 34(9):713
8 Hamalainen MS and Ilmoniemi R J. Interpreting measured magnetic fields of the brain: estimates of current distributions. Helsinki University of Technology. Finland, Technical Report TKK-F-A559, 1984
9 Wang J Z, Williamson S J, and 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
, http://www.100md.com
10 Loanides A, Bolton J, and Clarke C. Continuous probabilistic solutions to the biomagnetic inverse problem. Inverse Problems, 1990, 6:523
11 Corodnitsky I F, George J S, and Rao B D. Neuromagnetic source imaging with FOCUSS: a recursive weighted minimum norm algorithm. EEG and Clinical Neurophysiol, 1995, 95:231
12 Kanta Matsuura and Yoichi okabe. Selective minimum-norm solution of the biomagnetic inverse problem. IEEE Transactions on Biomedical Engineering, 1995, 42(6):608
, 百拇医药
13 Clarke C J S.Inverse Problems,1989, 5:999. Neuroscience and Tech., Lyon, France, Nov., 1992. on Biomedical Engineering, 1995, 42(6):60
14 Leahy R M, Mosher J C and Phillps J W. Proceedings of the 10th International Conference on Biomagnetism, BIOMAG 96, Sanata Fe, New Mexico, February 1996
15 Zhao S R, Horst Halling. Minimum L1 norm MEG reconstruction minimising signal deviation using a reduced lead field. Proceedings of the 18th Annual International Conference IEEE Engineering in Medicine and Biology Society, Part 2(of 5). (or Annual Interational Conference of the IEEE Engineering in Medicine and Biology Proceedings V 2 1996. IEEE, Piscataway, NJ, USA, 96CB36036. p 814-816)
(1998-06-16收稿，1998-12-03修回), http://www.100md.com