利用梯度算法估计脑磁源参数
作者：李军 刘正东
单位：李军(浙江大学信息学院光与电磁波研究中心 杭州 310027)；刘正东(浙江大学理学院物理系 杭州 310027)
关键词：脑磁图；参数估计；梯度算法
北京生物医学工程000304 摘 要 脑磁技术是生物电磁研究的一个重要领域。通过对脑磁图的反演可以获得脑神经的功能信息，进而为认知研究及脑神经疾病的诊断提供一种无损探测手段。这里利用模拟计算的方法对梯度算法在推断磁源参数中的应用进行了讨论。在假设已知磁源数目的条件下，对具有两种不同噪声水平(SNR＝20200)的磁场分别进行了计算。结果表明在反演一个磁源的磁场时，两种信噪比下均可以得到理想的结果，而反演两个磁源的磁场时必须有较高的信噪比。在事先不能确定磁源数目时，对反演模型中选择的源数目与实际源数目不一致时的情况亦进行了计算，并根据反演结果中出现的不同情况给出了相应的判别准则。
, 百拇医药
MEG's Source Parameter Estimation by Gradient Algorithm
Li Jun
(Center of Electromagnetic Wave and Light Research, College of Information
Zhejiang University, Hangzhou 310027)
Liu Zhengdong
(Physics Department, College of Science, Zhejiang University, 310027)
Abstract
, 百拇医药
Magnetoencephalography study is an important area of biomagnetic research.By solving the neuromagnetic inverse problem, one can obtain functional information in living human brain, which is a noninvasive detecting technique for the study of understanding and diagnosis of several neural diseases. In this paper, the application of gradient algorithm for magnetic source parameter estimation is discussed.In our computer simulation, Gaussian white noise is added to the field data with different signal-to-noise level(SNR=20 or 200).If the magnetic source number is known to be one, the estimated result is quite well with either of the two SNR levels. If the source number is two, reliable estimation can only be done with high SNR level(SNR≥200).In many cases, when the source number is not known in advance, the selected number in model may not be equal to that of the actual one. In respect of this problem, a calculation is also carried out, and a discrimination criteria is proposed based on results of inverse calculation.
, 百拇医药
Key words:Megnetoencephalogram, Parameter estimation, Gradient algorithm
0 前 言
人体内的生物电流可以激发出生物磁场，这种磁场很弱，目前利用超导量子干涉仪(superconducting quantum interference device, SQUID)可以拾取这一磁场信号[1]。脑磁为其中的重要一类。通过记录和分析这种由于脑内神经元的电活动而产生的脑磁图(magnetoencephalogram, MEG)，可以获得活动神经元的位置及强度等信息，为认知研究与脑神经疾病的诊断提供帮助[2，3]，目前结合核磁共振的脑磁源成像技术作为一种高时间、空间分辨率的无损伤功能成像技术正在迅速发展。
在脑磁的基础研究中，磁源参数的推断是一个重要的问题。目前针对这一问题存在两类解决的方法：一类称为磁源重建[4，5]，这是基于图像重建技术的方法。在这类方法中，假设磁源位于离散化的固定位置处，问题转化为对欠定的线性方程组的求解，其中包括最小范数解及加权最小范数解等。另一类称为偶极子定位方法[6，7]，直接面对磁场方程的非线性，通过非线性优化求解磁源参数。考虑到在目前的脑磁研究中，大多是处理少数几个源的场。因此相对全局优化方法来讲，非线性局域优化方法应是一种经济实用的选择。本文讨论了梯度算法在推断磁源参数中的应用，包括通过对不同信噪比下模拟磁场的反演计算，以确定反演时可以接受的噪声水平。讨论了在反演模型中假设的源数目与实际源数目不一致时出现的情况，及相应分辨的准则。另外对迭代初值的选择方法及收敛到真解的可能性亦进行了一些探讨。
, http://www.100md.com
1 计算模型与步骤
头模型采用电导率呈球对称分布的导体球，相应于脑电问题中的多层球模型。由于在整个脑内磁导率μ近似为常数，所以同样在这种情况下对电磁源参数的推断，利用脑磁比利用脑电要简单得多。磁源以电流偶极子P为模型，脑磁探测仪SQUID探头通常测量的是头表面径向磁感应强度值Br，计算公式为[8]：
(1)
式(1)中，R为测量场点位置，r为磁源位置，直角坐标下(1)式化为
(2)
式(2)中，(X，Y，Z)为磁源位置坐标，(X′，Y′，Z′)为测量场点位置坐标，
, 百拇医药
Pθ，Pφ为偶极子P的两个切向分量。
在对磁源参数估计时，采用极大似然方法。假设各测量点处噪声无关且呈具有相同方差的高斯分布，则磁源参数β(为由源位置及强度参数组成的列向量)应使目标函数F(β)＝
为最小。其中B为多点处磁场的测量值，B′(β)为(1)式代入参数β时的计算值，可见这里参数β的估计即为求最优化问题minβF(β)的解。对这一最优化问题，本文采用局域优化方法中的梯度算法。
梯度算法采用的迭代公式为：β(K+1)＝β(K)+λKPK,其中PK=-JT(β(K))e(β(K))，残差向量e(β(K))＝B-B′(β(K))，J为残差e关于β在β(K)处的Jacobi矩阵。步长因子λK满足：
，可采用一维搜索的方法确定步长因子λK，也可采用其它方法(如下列计算步骤中使用的)，只要沿负梯度方向即可。
, 百拇医药
计算时采用以下几个步骤：
第一步：在解空间内取初始点β(0)，置K=0
第二步：求目标函数F(β)在β＝β(K)的梯度
计算梯度的范数‖gK‖，若‖gK‖<ε或K>K0，则输出β，否则令PK＝-gK。
第三步：取β(K+1)＝β(K)+λKPK
置K＝K+1，返回第二步。
, 百拇医药
在上述第二步中，ε为梯度值接近零的精度，K0为设定的迭代步数。第三步中取步长因子
从以上步骤可以看到，梯度算法搜寻的是极值点，可能是某些局部极值点，也可能是全局最优点，所以具体应用时，尚须考虑对应点的目标函数情况。
2 模拟计算结果与讨论
设头半径为11cm，磁导率系数μ＝1.26×10-6H/m，采用头表面20个点处的磁场数据，分别以单电流偶极子及双电流偶极子为模拟磁场的磁源。磁源参数见表1、表2，其中位置坐标以厘米(cm)、偶极矩以纳安米(nAm)表示，以下情况同。
表1 单电流偶极子磁源参数 X(cm)
Y(cm)
, 百拇医药
Z(cm)
Pθ(nAm)
Pφ(nAm)
4.500
4.500
6.364
0.300
1.000
表2 双电流偶极子磁源参数 源
X(cm)
Y(cm)
Z(cm)
, http://www.100md.com
Pθ(nAm)
Pφ(nAm)
P1
4.500
4.500
6.363
0.300
1.000
P2
-4.950
4.950
0.000
, http://www.100md.com
0.500
0.800
为模拟真实的磁场情况，在上述源产生的径向磁感应强度数值中，加入一定能量的高斯噪声，并假设各测量点处的噪声具有相同的标准差且彼此不相关。这里定义能量为磁场强度值的平方，信噪比SNR为所有测量点处磁场的总能量除以噪声的总能量。
首先假设已知磁源的数目，在两种信噪比情况下，分别对单个源及两个源的磁场进行反演计算。结果见表3。表3 两种信噪比下磁源参数的反演计算结果 源
数
目
信
噪
比
, 百拇医药
源
迭代初值
磁源参数计算值
目标函数
相对值
(目标函数
/磁场能量)
X
Y
Z
Pθ
Pφ
X
, http://www.100md.com
Y
Z
Pθ
Pφ
1
20
P1
1.907
-1.484
1.133
0.784
3.673
4.488
, 百拇医药
4.050
6.147
0.299
1.118
0.034
P1
1.970
-1.075
1.784
0.313
2.591
4.488
, 百拇医药 4.050
6.147
0.299
1.118
0.034
200
P1
-0.085
-2.260
2.170
-0.733
3.745
4.498
, 百拇医药
4.356
6.298
0.301
1.035
0.004
P1
2.509
-5.393
3.067
-0.117
1.642
4.498
, 百拇医药 4.356
6.298
0.301
1.035
0.004
2
20
P1
-0.768
0.911
-0.669
0.716
2.639
, 百拇医药
0.808
4.543
-2.770
-0.834
0.460
P2
0.771
3.080
0.002
1.997
4.049
2.827
, http://www.100md.com 3.482
4.319
1.028
2.034
0.023
P1
-0.300
0.898
-9.690
0.154
-0.180
-1.613
1.935
, http://www.100md.com
-10.669
0.013
-0.314
P2
-2.396
1.860
0.175
1.819
4.204
3.494
3.651
4.983
, 百拇医药 0.740
1.648
0.038
200
P1
4.271
-0.257
4.040
0.326
2.273
4.547
4.419
6.392
, http://www.100md.com
0.327
0.976
P2
-1.425
-2.874
0.171
-0.185
1.448
-4.367
4.580
0.103
0.431
, 百拇医药 0.907
0.0014
P1
-0.080
-0.336
0.023
0.311
1.718
4.547
4.419
6.392
0.327
0.976
, 百拇医药
P2
-0.012
5.358
1.347
-0.197
1.919
-4.367
4.580
0.103
0.431
0.907
0.0014
, 百拇医药 在计算表3时，对于单个源情况，迭代初始值是通过在解空间中随机产生100个点，选取其中目标函数为最小的而产生的。对照表1，可以看到在SNR＝200时的反演结果优于SNR＝20时的。从目标函数相对值中可以看到，这一数值与噪声水平一致，具有相同的数量级，表明此时误差来源于噪声。对于两个源的计算，迭代初始点是从1000个随机点中选取的。在SNR＝20时，尽管目标函数相对值与噪声水平相近，但磁源参数的计算值已有很大的偏差了，说明在这种噪声水平下，对源参数的推断已不可靠。在SNR＝200时，结果尚可，但已不如同样信噪比情况下对单个磁源的反演结果。可见在利用梯度法对磁源参数计算时，对于单个源情况，信噪比要求相对不高，SNR＝20以上即可，但计算两个源时，信噪比必须大于200。在实际处理脑磁数据时，比如对于诱发脑磁实验，一般情况下是通过相同条件下的多次测量取平均以压制噪声，提高信噪比，从而提高对源的反演精度。
在梯度算法使用中，初值的选取方式及相应收敛到真解处的可能性，也是必须考虑的问题，因为这也关系到算法是否实用。对此我们进行了估算，在SNR＝200时，在解空间中随机产生1000个点作为迭代初值，计算最后搜寻到的点。结果显示对于一个源的情况，有756个点是收敛在真解附近的，而对于两个源的情况则仅有86个点收敛在真解附近。可见在对一个源反演时，迭代初值很容易产生，而对于两个源的情况，则须仔细选择，由此进一步推断对于更多源的情况，迭代初值将很难产生。从这一点上看，梯度算法适合于对含有一两个源磁场的反演计算。
, http://www.100md.com
以上假设反演前已知源的数目，其实在实际情况下，有时并不能先验地确定磁源数目，这样在反演模型中则可能出现采用的源数目与真实的情况不一致。以下就这一问题做一些计算与讨论。首先以一个源为假设反演两个源的场，取信噪比为200，结果见表4。在表4中三次计算的初值分别采用在100100010000个随机点中选出一个目标函数最小的为迭代初值。从表4可以看到，从三个不同初值开始的搜寻均收敛在一个点处，表明该点为极值点，显然该点不是真正的解。实际情况下又如何做这样的断定呢?通过观察目标函数相对值情况，可以看到此时该值远高于噪声水平，两者相差一个数量级，由此可以推断这一差异不是源于噪声，而是源于计算模型的不正确，即采用了少于实际数目的源作为反演的源假设，造成了目标函数无法最小优化。
表4 以一个磁源为假设反演两个磁源磁场的结果(SNR＝200) 迭代初值选择方式
迭代初值
磁源参数计算值
, http://www.100md.com
目标函数
相对值
(目标函数
/磁场能量)
X
Y
Z
Pθ
Pφ
X
Y
Z
Pθ
, 百拇医药
Pφ
100∶1
1.733
2.404
1.287
-0.974
4.973
2.148
4.056
4.066
0.704
2.134
0.062
, http://www.100md.com
1000∶1
2.155
-0.312
-0.789
-0.472
5.586
2.148
4.056
4.066
0.704
2.134
0.062
10000∶1
, http://www.100md.com
1.240
2.277
0.759
1.540
6.429
2.148
4.056
4.066
0.704
2.134
0.062
注：在迭代初值选择方式列中，100∶1表示从解空间随机产生的100个点中，选择目标函数最小的为迭代初值。其它类同。
, 百拇医药
另外在反演时对于源数目的选择，还可能出现另一种情况，即模型中采用的数目多于实际的，对此我们用以两个源为假设反演一个源的情况为例，加以讨论说明。设源参数取自表1，信噪比为200，计算结果见表5。迭代初值来自不同的两组1000个点的样本。
表5 以两个磁源为假设反演一个磁源磁场的结果(SNR＝200) 源
迭代初值
磁源参数计算值
目标函数
相对值
(目标函数
/磁场能量)
计算源的
, http://www.100md.com
相对能量
(计算源
能量/磁
场能量)
X
Y
Z
Pθ
Pφ
X
Y
Z
Pθ
, 百拇医药
Pφ
P1
-0.483
-1.845
0.126
0.311
1.781
-0.423
-0.387
1.745
0.249
0.403
, 百拇医药 0.0022
P2
-0.012
5.358
1.347
-0.197
1.919
4.426
4.299
6.228
0.356
1.012
0.0022
, 百拇医药
0.9912
P1
0.689
-0.579
-1.988
-0.802
-1.919
-0.001
-0.001
-0.278
1.353
-2.655
, 百拇医药
0.0017
P2
0.521
2.039
0.894
1.973
5.333
4.418
4.300
6.215
0.350
1.031
, 百拇医药 0.0023
1.0149
观察表5可以发现，两个计算的源中，有一个接近真实的源，而另一个源在两次不同初值的计算结果中是不一致的。同时目标函数均可最小优化，与噪声水平在同一量级上。实际中如何发现并处理这种情况呢?一种方法可以通过几次不同初值的计算，观察解中是否存在相对稳定的点及变化不定的点，若出现这种情况，则表示变化不定的点对应多余的伪源。另一种方法可以通过计算这两个源的能量在整个场能量的比例来进行真伪的区分，从表5可以看到，伪源的相对能量极小，而磁场的能量几乎全部来自于另一个源，所以也可以只通过一次计算，根据两个源对磁场能量的贡献不同而去掉伪源。
在一些实际情况下，若磁源的数目较少，则可以利用梯度法进行源参数的推断，这样在计算时间上远优于一般的全局优化方法。根据上述模拟计算的结果与讨论，在实际测量脑磁数据及反演脑磁源参数时，可以进行如下操作：首先对各脑磁场测量点处的噪声进行测量以获得它的统计特性与能量，然后在外界条件相同的情况下，对脑磁数据进行多次测量并求平均，使数据具有较高的信噪比水平(如200以上)。然后可以进行反演计算，在模型中采用一个源反演时，若出现目标函数无法优化到噪声能量水平的情况，则表示模型中所采用的源数目少于实际的源数目，这时可采用增大源数目的方法。模型中采用两个源反演时，在目标函数优化到噪声能量水平的情况下，若出现其中一个计算源对磁场能量的贡献极小，则表示该源是一个伪源而真实的源只有一个，同时可以断定对磁场贡献大的另一个源接近真实的源，若两个计算出的源对磁场能量均有较大贡献，则表示实际即为两个源且反演是成功的。
, 百拇医药
基金项目：本项研究得到浙江省自然科学基金(No.699025)的资助。
作者简介：李军(1965—)，男，籍贯：辽宁沈阳，博士，现在浙江大学信息学院光与电磁波研究中心做博士后，主要从事脑电磁及生物光学方向研究。
3 参考文献
[1] Hamalainen M, Hari R, et al.Magnetoencephalography theory, instrumentation, and applications to noninvasive studies of the working human brain.Rev Mod Phys,1993,65(2):413
[2] Hari R,Lounasmaa O V.Recording and interpretation of cerebral magnetic fields. Science,1989,244:432
, 百拇医药
[3] Forss N.Magnetoencephalography(MEG)in epilepsy surgery.Acta Neurochur suppl,1997,58:81
[4] Phillips J W,Leahy R M,Mosher J C.MEG-based imaging of focal neuronal current sources.IEEE Trans Biomed Eng,1997,16(3):338-348
[5] Wang J Z,Williamson S J,Kaufman L.Magnetic source images determined by a lead-field analysis:the unique minimum-norm least-squares estimation. IEEE Trans Biomed Eng,1992,39(7)665-675
, 百拇医药
[6] Mosher J C,Lewis P S, Leahy R M.Multiple dipole modeling and localization from spatio-temporal MEG data.IEEE Trans Biomed Eng,1992,39(6):541-557
[7] Sekihara K, Poeppel D,Marantz A, et al.Noise covariance incorporated MEG-MUSIC algorithm:a method for multiple-dipole estimation tolerant of the influnence of background brain activity.IEEE Trans.Biomed Eng,1997,44(9):839-847
[8] Sarvas J. Basic mathematical and electromagnetic concepts of the biomagnetic inverse problem.Phys Med Biol,1987,32(1):11-22
(1999-09-02收稿), 百拇医药
单位：李军(浙江大学信息学院光与电磁波研究中心 杭州 310027)；刘正东(浙江大学理学院物理系 杭州 310027)
关键词：脑磁图；参数估计；梯度算法
北京生物医学工程000304 摘 要 脑磁技术是生物电磁研究的一个重要领域。通过对脑磁图的反演可以获得脑神经的功能信息，进而为认知研究及脑神经疾病的诊断提供一种无损探测手段。这里利用模拟计算的方法对梯度算法在推断磁源参数中的应用进行了讨论。在假设已知磁源数目的条件下，对具有两种不同噪声水平(SNR＝20200)的磁场分别进行了计算。结果表明在反演一个磁源的磁场时，两种信噪比下均可以得到理想的结果，而反演两个磁源的磁场时必须有较高的信噪比。在事先不能确定磁源数目时，对反演模型中选择的源数目与实际源数目不一致时的情况亦进行了计算，并根据反演结果中出现的不同情况给出了相应的判别准则。
, 百拇医药
MEG's Source Parameter Estimation by Gradient Algorithm
Li Jun
(Center of Electromagnetic Wave and Light Research, College of Information
Zhejiang University, Hangzhou 310027)
Liu Zhengdong
(Physics Department, College of Science, Zhejiang University, 310027)
Abstract
, 百拇医药
Magnetoencephalography study is an important area of biomagnetic research.By solving the neuromagnetic inverse problem, one can obtain functional information in living human brain, which is a noninvasive detecting technique for the study of understanding and diagnosis of several neural diseases. In this paper, the application of gradient algorithm for magnetic source parameter estimation is discussed.In our computer simulation, Gaussian white noise is added to the field data with different signal-to-noise level(SNR=20 or 200).If the magnetic source number is known to be one, the estimated result is quite well with either of the two SNR levels. If the source number is two, reliable estimation can only be done with high SNR level(SNR≥200).In many cases, when the source number is not known in advance, the selected number in model may not be equal to that of the actual one. In respect of this problem, a calculation is also carried out, and a discrimination criteria is proposed based on results of inverse calculation.
, 百拇医药
Key words:Megnetoencephalogram, Parameter estimation, Gradient algorithm
0 前 言
人体内的生物电流可以激发出生物磁场，这种磁场很弱，目前利用超导量子干涉仪(superconducting quantum interference device, SQUID)可以拾取这一磁场信号[1]。脑磁为其中的重要一类。通过记录和分析这种由于脑内神经元的电活动而产生的脑磁图(magnetoencephalogram, MEG)，可以获得活动神经元的位置及强度等信息，为认知研究与脑神经疾病的诊断提供帮助[2，3]，目前结合核磁共振的脑磁源成像技术作为一种高时间、空间分辨率的无损伤功能成像技术正在迅速发展。
在脑磁的基础研究中，磁源参数的推断是一个重要的问题。目前针对这一问题存在两类解决的方法：一类称为磁源重建[4，5]，这是基于图像重建技术的方法。在这类方法中，假设磁源位于离散化的固定位置处，问题转化为对欠定的线性方程组的求解，其中包括最小范数解及加权最小范数解等。另一类称为偶极子定位方法[6，7]，直接面对磁场方程的非线性，通过非线性优化求解磁源参数。考虑到在目前的脑磁研究中，大多是处理少数几个源的场。因此相对全局优化方法来讲，非线性局域优化方法应是一种经济实用的选择。本文讨论了梯度算法在推断磁源参数中的应用，包括通过对不同信噪比下模拟磁场的反演计算，以确定反演时可以接受的噪声水平。讨论了在反演模型中假设的源数目与实际源数目不一致时出现的情况，及相应分辨的准则。另外对迭代初值的选择方法及收敛到真解的可能性亦进行了一些探讨。
, http://www.100md.com
1 计算模型与步骤
头模型采用电导率呈球对称分布的导体球，相应于脑电问题中的多层球模型。由于在整个脑内磁导率μ近似为常数，所以同样在这种情况下对电磁源参数的推断，利用脑磁比利用脑电要简单得多。磁源以电流偶极子P为模型，脑磁探测仪SQUID探头通常测量的是头表面径向磁感应强度值Br，计算公式为[8]：
(1)式(1)中，R为测量场点位置，r为磁源位置，直角坐标下(1)式化为
(2)式(2)中，(X，Y，Z)为磁源位置坐标，(X′，Y′，Z′)为测量场点位置坐标，
, 百拇医药
Pθ，Pφ为偶极子P的两个切向分量。
在对磁源参数估计时，采用极大似然方法。假设各测量点处噪声无关且呈具有相同方差的高斯分布，则磁源参数β(为由源位置及强度参数组成的列向量)应使目标函数F(β)＝
为最小。其中B为多点处磁场的测量值，B′(β)为(1)式代入参数β时的计算值，可见这里参数β的估计即为求最优化问题minβF(β)的解。对这一最优化问题，本文采用局域优化方法中的梯度算法。梯度算法采用的迭代公式为：β(K+1)＝β(K)+λKPK,其中PK=-JT(β(K))e(β(K))，残差向量e(β(K))＝B-B′(β(K))，J为残差e关于β在β(K)处的Jacobi矩阵。步长因子λK满足：
，可采用一维搜索的方法确定步长因子λK，也可采用其它方法(如下列计算步骤中使用的)，只要沿负梯度方向即可。, 百拇医药
计算时采用以下几个步骤：
第一步：在解空间内取初始点β(0)，置K=0
第二步：求目标函数F(β)在β＝β(K)的梯度
计算梯度的范数‖gK‖，若‖gK‖<ε或K>K0，则输出β，否则令PK＝-gK。
第三步：取β(K+1)＝β(K)+λKPK
置K＝K+1，返回第二步。, 百拇医药
在上述第二步中，ε为梯度值接近零的精度，K0为设定的迭代步数。第三步中取步长因子
从以上步骤可以看到，梯度算法搜寻的是极值点，可能是某些局部极值点，也可能是全局最优点，所以具体应用时，尚须考虑对应点的目标函数情况。2 模拟计算结果与讨论
设头半径为11cm，磁导率系数μ＝1.26×10-6H/m，采用头表面20个点处的磁场数据，分别以单电流偶极子及双电流偶极子为模拟磁场的磁源。磁源参数见表1、表2，其中位置坐标以厘米(cm)、偶极矩以纳安米(nAm)表示，以下情况同。
表1 单电流偶极子磁源参数 X(cm)
Y(cm)
, 百拇医药
Z(cm)
Pθ(nAm)
Pφ(nAm)
4.500
4.500
6.364
0.300
1.000
表2 双电流偶极子磁源参数 源
X(cm)
Y(cm)
Z(cm)
, http://www.100md.com
Pθ(nAm)
Pφ(nAm)
P1
4.500
4.500
6.363
0.300
1.000
P2
-4.950
4.950
0.000
, http://www.100md.com
0.500
0.800
为模拟真实的磁场情况，在上述源产生的径向磁感应强度数值中，加入一定能量的高斯噪声，并假设各测量点处的噪声具有相同的标准差且彼此不相关。这里定义能量为磁场强度值的平方，信噪比SNR为所有测量点处磁场的总能量除以噪声的总能量。
首先假设已知磁源的数目，在两种信噪比情况下，分别对单个源及两个源的磁场进行反演计算。结果见表3。表3 两种信噪比下磁源参数的反演计算结果 源
数
目
信
噪
比
, 百拇医药
源
迭代初值
磁源参数计算值
目标函数
相对值
(目标函数
/磁场能量)
X
Y
Z
Pθ
Pφ
X
, http://www.100md.com
Y
Z
Pθ
Pφ
1
20
P1
1.907
-1.484
1.133
0.784
3.673
4.488
, 百拇医药
4.050
6.147
0.299
1.118
0.034
P1
1.970
-1.075
1.784
0.313
2.591
4.488
, 百拇医药 4.050
6.147
0.299
1.118
0.034
200
P1
-0.085
-2.260
2.170
-0.733
3.745
4.498
, 百拇医药
4.356
6.298
0.301
1.035
0.004
P1
2.509
-5.393
3.067
-0.117
1.642
4.498
, 百拇医药 4.356
6.298
0.301
1.035
0.004
2
20
P1
-0.768
0.911
-0.669
0.716
2.639
, 百拇医药
0.808
4.543
-2.770
-0.834
0.460
P2
0.771
3.080
0.002
1.997
4.049
2.827
, http://www.100md.com 3.482
4.319
1.028
2.034
0.023
P1
-0.300
0.898
-9.690
0.154
-0.180
-1.613
1.935
, http://www.100md.com
-10.669
0.013
-0.314
P2
-2.396
1.860
0.175
1.819
4.204
3.494
3.651
4.983
, 百拇医药 0.740
1.648
0.038
200
P1
4.271
-0.257
4.040
0.326
2.273
4.547
4.419
6.392
, http://www.100md.com
0.327
0.976
P2
-1.425
-2.874
0.171
-0.185
1.448
-4.367
4.580
0.103
0.431
, 百拇医药 0.907
0.0014
P1
-0.080
-0.336
0.023
0.311
1.718
4.547
4.419
6.392
0.327
0.976
, 百拇医药
P2
-0.012
5.358
1.347
-0.197
1.919
-4.367
4.580
0.103
0.431
0.907
0.0014
, 百拇医药 在计算表3时，对于单个源情况，迭代初始值是通过在解空间中随机产生100个点，选取其中目标函数为最小的而产生的。对照表1，可以看到在SNR＝200时的反演结果优于SNR＝20时的。从目标函数相对值中可以看到，这一数值与噪声水平一致，具有相同的数量级，表明此时误差来源于噪声。对于两个源的计算，迭代初始点是从1000个随机点中选取的。在SNR＝20时，尽管目标函数相对值与噪声水平相近，但磁源参数的计算值已有很大的偏差了，说明在这种噪声水平下，对源参数的推断已不可靠。在SNR＝200时，结果尚可，但已不如同样信噪比情况下对单个磁源的反演结果。可见在利用梯度法对磁源参数计算时，对于单个源情况，信噪比要求相对不高，SNR＝20以上即可，但计算两个源时，信噪比必须大于200。在实际处理脑磁数据时，比如对于诱发脑磁实验，一般情况下是通过相同条件下的多次测量取平均以压制噪声，提高信噪比，从而提高对源的反演精度。
在梯度算法使用中，初值的选取方式及相应收敛到真解处的可能性，也是必须考虑的问题，因为这也关系到算法是否实用。对此我们进行了估算，在SNR＝200时，在解空间中随机产生1000个点作为迭代初值，计算最后搜寻到的点。结果显示对于一个源的情况，有756个点是收敛在真解附近的，而对于两个源的情况则仅有86个点收敛在真解附近。可见在对一个源反演时，迭代初值很容易产生，而对于两个源的情况，则须仔细选择，由此进一步推断对于更多源的情况，迭代初值将很难产生。从这一点上看，梯度算法适合于对含有一两个源磁场的反演计算。
, http://www.100md.com
以上假设反演前已知源的数目，其实在实际情况下，有时并不能先验地确定磁源数目，这样在反演模型中则可能出现采用的源数目与真实的情况不一致。以下就这一问题做一些计算与讨论。首先以一个源为假设反演两个源的场，取信噪比为200，结果见表4。在表4中三次计算的初值分别采用在100100010000个随机点中选出一个目标函数最小的为迭代初值。从表4可以看到，从三个不同初值开始的搜寻均收敛在一个点处，表明该点为极值点，显然该点不是真正的解。实际情况下又如何做这样的断定呢?通过观察目标函数相对值情况，可以看到此时该值远高于噪声水平，两者相差一个数量级，由此可以推断这一差异不是源于噪声，而是源于计算模型的不正确，即采用了少于实际数目的源作为反演的源假设，造成了目标函数无法最小优化。
表4 以一个磁源为假设反演两个磁源磁场的结果(SNR＝200) 迭代初值选择方式
迭代初值
磁源参数计算值
, http://www.100md.com
目标函数
相对值
(目标函数
/磁场能量)
X
Y
Z
Pθ
Pφ
X
Y
Z
Pθ
, 百拇医药
Pφ
100∶1
1.733
2.404
1.287
-0.974
4.973
2.148
4.056
4.066
0.704
2.134
0.062
, http://www.100md.com
1000∶1
2.155
-0.312
-0.789
-0.472
5.586
2.148
4.056
4.066
0.704
2.134
0.062
10000∶1
, http://www.100md.com
1.240
2.277
0.759
1.540
6.429
2.148
4.056
4.066
0.704
2.134
0.062
注：在迭代初值选择方式列中，100∶1表示从解空间随机产生的100个点中，选择目标函数最小的为迭代初值。其它类同。
, 百拇医药
另外在反演时对于源数目的选择，还可能出现另一种情况，即模型中采用的数目多于实际的，对此我们用以两个源为假设反演一个源的情况为例，加以讨论说明。设源参数取自表1，信噪比为200，计算结果见表5。迭代初值来自不同的两组1000个点的样本。
表5 以两个磁源为假设反演一个磁源磁场的结果(SNR＝200) 源
迭代初值
磁源参数计算值
目标函数
相对值
(目标函数
/磁场能量)
计算源的
, http://www.100md.com
相对能量
(计算源
能量/磁
场能量)
X
Y
Z
Pθ
Pφ
X
Y
Z
Pθ
, 百拇医药
Pφ
P1
-0.483
-1.845
0.126
0.311
1.781
-0.423
-0.387
1.745
0.249
0.403
, 百拇医药 0.0022
P2
-0.012
5.358
1.347
-0.197
1.919
4.426
4.299
6.228
0.356
1.012
0.0022
, 百拇医药
0.9912
P1
0.689
-0.579
-1.988
-0.802
-1.919
-0.001
-0.001
-0.278
1.353
-2.655
, 百拇医药
0.0017
P2
0.521
2.039
0.894
1.973
5.333
4.418
4.300
6.215
0.350
1.031
, 百拇医药 0.0023
1.0149
观察表5可以发现，两个计算的源中，有一个接近真实的源，而另一个源在两次不同初值的计算结果中是不一致的。同时目标函数均可最小优化，与噪声水平在同一量级上。实际中如何发现并处理这种情况呢?一种方法可以通过几次不同初值的计算，观察解中是否存在相对稳定的点及变化不定的点，若出现这种情况，则表示变化不定的点对应多余的伪源。另一种方法可以通过计算这两个源的能量在整个场能量的比例来进行真伪的区分，从表5可以看到，伪源的相对能量极小，而磁场的能量几乎全部来自于另一个源，所以也可以只通过一次计算，根据两个源对磁场能量的贡献不同而去掉伪源。
在一些实际情况下，若磁源的数目较少，则可以利用梯度法进行源参数的推断，这样在计算时间上远优于一般的全局优化方法。根据上述模拟计算的结果与讨论，在实际测量脑磁数据及反演脑磁源参数时，可以进行如下操作：首先对各脑磁场测量点处的噪声进行测量以获得它的统计特性与能量，然后在外界条件相同的情况下，对脑磁数据进行多次测量并求平均，使数据具有较高的信噪比水平(如200以上)。然后可以进行反演计算，在模型中采用一个源反演时，若出现目标函数无法优化到噪声能量水平的情况，则表示模型中所采用的源数目少于实际的源数目，这时可采用增大源数目的方法。模型中采用两个源反演时，在目标函数优化到噪声能量水平的情况下，若出现其中一个计算源对磁场能量的贡献极小，则表示该源是一个伪源而真实的源只有一个，同时可以断定对磁场贡献大的另一个源接近真实的源，若两个计算出的源对磁场能量均有较大贡献，则表示实际即为两个源且反演是成功的。
, 百拇医药
基金项目：本项研究得到浙江省自然科学基金(No.699025)的资助。
作者简介：李军(1965—)，男，籍贯：辽宁沈阳，博士，现在浙江大学信息学院光与电磁波研究中心做博士后，主要从事脑电磁及生物光学方向研究。
3 参考文献
[1] Hamalainen M, Hari R, et al.Magnetoencephalography theory, instrumentation, and applications to noninvasive studies of the working human brain.Rev Mod Phys,1993,65(2):413
[2] Hari R,Lounasmaa O V.Recording and interpretation of cerebral magnetic fields. Science,1989,244:432
, 百拇医药
[3] Forss N.Magnetoencephalography(MEG)in epilepsy surgery.Acta Neurochur suppl,1997,58:81
[4] Phillips J W,Leahy R M,Mosher J C.MEG-based imaging of focal neuronal current sources.IEEE Trans Biomed Eng,1997,16(3):338-348
[5] Wang J Z,Williamson S J,Kaufman L.Magnetic source images determined by a lead-field analysis:the unique minimum-norm least-squares estimation. IEEE Trans Biomed Eng,1992,39(7)665-675
, 百拇医药
[6] Mosher J C,Lewis P S, Leahy R M.Multiple dipole modeling and localization from spatio-temporal MEG data.IEEE Trans Biomed Eng,1992,39(6):541-557
[7] Sekihara K, Poeppel D,Marantz A, et al.Noise covariance incorporated MEG-MUSIC algorithm:a method for multiple-dipole estimation tolerant of the influnence of background brain activity.IEEE Trans.Biomed Eng,1997,44(9):839-847
[8] Sarvas J. Basic mathematical and electromagnetic concepts of the biomagnetic inverse problem.Phys Med Biol,1987,32(1):11-22
(1999-09-02收稿), 百拇医药