JPH0515504A - 磁場発生源の推定方法 - Google Patents
磁場発生源の推定方法Info
- Publication number
- JPH0515504A JPH0515504A JP3167883A JP16788391A JPH0515504A JP H0515504 A JPH0515504 A JP H0515504A JP 3167883 A JP3167883 A JP 3167883A JP 16788391 A JP16788391 A JP 16788391A JP H0515504 A JPH0515504 A JP H0515504A
- Authority
- JP
- Japan
- Prior art keywords
- magnetic field
- estimating
- covariance matrix
- noise
- measurement
- 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.)
- Pending
Links
Landscapes
- Measuring Magnetic Variables (AREA)
- Measurement And Recording Of Electrical Phenomena And Electrical Characteristics Of The Living Body (AREA)
Abstract
(57)【要約】
【目的】 複数の検出コイルを持ち、複数の測定点を同
時に測定して、得られた生体表面の磁場データから生体
内の磁場源を推定して表示する装置において、電流分布
の推定精度を向上させる。 【構成】 生体磁場計測において、生体外の磁場源から
の磁場はノイズとして働き、生体内の磁場源推定精度を
低下させる。本発明では、複数の検出コイルにより複数
点の磁場を同時計測する場合、上記ノイズはチャネル間
で相関を持つことに着目し、ノイズの共分散行列Ｃを予
め推定しておき（１０２，１０３）、式（３）に示す対
数尤度Ｅを最大とする磁場源パラメータを求める（１０
４）。
時に測定して、得られた生体表面の磁場データから生体
内の磁場源を推定して表示する装置において、電流分布
の推定精度を向上させる。 【構成】 生体磁場計測において、生体外の磁場源から
の磁場はノイズとして働き、生体内の磁場源推定精度を
低下させる。本発明では、複数の検出コイルにより複数
点の磁場を同時計測する場合、上記ノイズはチャネル間
で相関を持つことに着目し、ノイズの共分散行列Ｃを予
め推定しておき（１０２，１０３）、式（３）に示す対
数尤度Ｅを最大とする磁場源パラメータを求める（１０
４）。
Description
【０００１】
【産業上の利用分野】本発明は、生体表面の磁場を観測
し、そのデータから生体内の電流分布を推定して表示す
る装置における電流分布の推定方法に関し、特に多チャ
ネルを用いた生体磁場計測時におけるチャネル間のノイ
ズの相関を考慮することによって、電流分布の推定精度
を向上させるのに好適な磁場発生源の推定方法に関す
る。
し、そのデータから生体内の電流分布を推定して表示す
る装置における電流分布の推定方法に関し、特に多チャ
ネルを用いた生体磁場計測時におけるチャネル間のノイ
ズの相関を考慮することによって、電流分布の推定精度
を向上させるのに好適な磁場発生源の推定方法に関す
る。
【０００２】
【従来の技術】従来、生体外部の磁場測定値から生体内
部の電流分布の推定は通常以下のように行われる。ここ
では、生体磁場計測でも代表的な脳磁計測を例にとって
述べる。例えば、座標系を図２のように定義する。図２
において、２１は測定位置を示し、ベクトルＲnの位置
にベクトルＱnで表わされた電流が存在すると仮定す
る。生体磁気計測の分野では、このような孤立した電流
ベクトルを仮定し、これを電流ダイポールと呼ぶ。ｎは
ダイポールに付けられた番号であり、ｎ＝１，…，Ｎと
する。また、Ｒmは測定点を表わす位置ベクトルであ
り、ｍは測定点に付けられた番号であって、ｍ＝１，
…，Ｍである。さて、脳内にＮ個の電流ダイポールを仮
定すると、Ｒmの点での磁束密度ベクトルの法線成分Ｂm
は
部の電流分布の推定は通常以下のように行われる。ここ
では、生体磁場計測でも代表的な脳磁計測を例にとって
述べる。例えば、座標系を図２のように定義する。図２
において、２１は測定位置を示し、ベクトルＲnの位置
にベクトルＱnで表わされた電流が存在すると仮定す
る。生体磁気計測の分野では、このような孤立した電流
ベクトルを仮定し、これを電流ダイポールと呼ぶ。ｎは
ダイポールに付けられた番号であり、ｎ＝１，…，Ｎと
する。また、Ｒmは測定点を表わす位置ベクトルであ
り、ｍは測定点に付けられた番号であって、ｍ＝１，
…，Ｍである。さて、脳内にＮ個の電流ダイポールを仮
定すると、Ｒmの点での磁束密度ベクトルの法線成分Ｂm
は
【数１】
で表わされる。Ｒmの位置で実際に測定される磁場の値
をＤmとし、第ｎ番目のダイポールの位置および電流ベ
クトルの推定値をＳＲn，ＳＱnとし、これらから式
（１）を用いて計算された位置Ｒmにおける磁場の法線
成分をＳＢmとすると、ＳＲn，ＳＱnの最適推定値は関
数
をＤmとし、第ｎ番目のダイポールの位置および電流ベ
クトルの推定値をＳＲn，ＳＱnとし、これらから式
（１）を用いて計算された位置Ｒmにおける磁場の法線
成分をＳＢmとすると、ＳＲn，ＳＱnの最適推定値は関
数
【数２】
を最小とするＳＲn，ＳＱn（ｎ＝１，２，…，Ｎ）とし
て求める。なお、式（２）は推定データと測定データと
の一致度を表わし、いわゆる最小２乗法を示すものであ
る。
て求める。なお、式（２）は推定データと測定データと
の一致度を表わし、いわゆる最小２乗法を示すものであ
る。
【０００３】
【発明が解決しようとする課題】上記従来技術における
最小２乗法は、統計学で最尤推定法と呼ばれる方法の特
別な場合、すなわち測定データＤmに混入しているノイ
ズが互いに無相関で、正規ノイズの場合に導かれるもの
である。一方、複数個の検出コイルにより複数の測定点
を同時に計測する磁束計（以下、多チャネル磁束計と呼
ぶ）においては、全チャネルで同時に測定するため、測
定領域外の磁気ノイズは全チャネルに入力し、各チャネ
ルで相関のあるノイズとなる。従って、最小２乗法をそ
のまま多チャネル磁束計で測定したデータに適用して
も、各チャネル間でのノイズの相関が考慮されていない
ため、推定結果に誤差を含む。本発明の目的は、ノイズ
のチャネル間での相関を考慮することによって、このよ
うな問題点を改善し、精度の高い電流分布推定が可能な
生体磁気計測装置における磁場源の推定方法を提供する
ことにある。
最小２乗法は、統計学で最尤推定法と呼ばれる方法の特
別な場合、すなわち測定データＤmに混入しているノイ
ズが互いに無相関で、正規ノイズの場合に導かれるもの
である。一方、複数個の検出コイルにより複数の測定点
を同時に計測する磁束計（以下、多チャネル磁束計と呼
ぶ）においては、全チャネルで同時に測定するため、測
定領域外の磁気ノイズは全チャネルに入力し、各チャネ
ルで相関のあるノイズとなる。従って、最小２乗法をそ
のまま多チャネル磁束計で測定したデータに適用して
も、各チャネル間でのノイズの相関が考慮されていない
ため、推定結果に誤差を含む。本発明の目的は、ノイズ
のチャネル間での相関を考慮することによって、このよ
うな問題点を改善し、精度の高い電流分布推定が可能な
生体磁気計測装置における磁場源の推定方法を提供する
ことにある。
【０００４】
【課題を解決するための手段】上記目的を達成するた
め、本発明の磁場発生源の推定方法は、対象とする生体
外部の複数の測定点にて、この生体から発生する磁場を
複数のチャネルを持つ磁束計により測定し、生体内の電
流分布を推定する方法において、各チャネルで測定した
データに混入しているノイズのチャネル間での相関を考
慮し、例えば、共分散行列の（ｍ，ｍ’）要素を推定す
る際、磁束計のｍ番目の検出コイル出力とｍ’番目の検
出コイル出力の積を時間平均した量を用いて、検出コイ
ル出力間でのノイズの共分散行列の（ｍ，ｍ’）要素を
推定し、その共分散行列の逆行列のｍｍ’の要素を示す
第１の量（Ａmm’）と、ｍ番目の検出器出力の測定値と
生体内の電流分布推定値から算出されたｍ番目の検出器
出力の差を示す第２の量（Ｄm−ＳＢm）と、ｍ’番目の
検出器出力の測定値と生体内の電流分布推定値から算出
されたｍ’番目の検出コイル出力の差を示す第３の量
（Ｄm’−ＳＢm’）との積を、ｍｍ’の全組み合わせに
ついて積分した量を含む対数尤度Ｅを求めて、そのＥを
最大とする磁場源パラメータ（ＳＲn,ＳＱｎ）により、
生体内の電流分布を推定することに特徴がある。また、
上記ノイズの共分散行列を推定する場合、心磁波形や自
発脳磁等の計測では、測定対象の生体磁場を検出しない
位置に磁束計を設置して計測を行い、誘発脳磁計測で
は、誘発刺激の付与以前から連続してデータの取り込み
を行うことに特徴があり、さらに、心磁波形では、心拍
の各周期の間にある信号の極めて微弱な部分からノイズ
の共分散行列を推定することに特徴がある。
め、本発明の磁場発生源の推定方法は、対象とする生体
外部の複数の測定点にて、この生体から発生する磁場を
複数のチャネルを持つ磁束計により測定し、生体内の電
流分布を推定する方法において、各チャネルで測定した
データに混入しているノイズのチャネル間での相関を考
慮し、例えば、共分散行列の（ｍ，ｍ’）要素を推定す
る際、磁束計のｍ番目の検出コイル出力とｍ’番目の検
出コイル出力の積を時間平均した量を用いて、検出コイ
ル出力間でのノイズの共分散行列の（ｍ，ｍ’）要素を
推定し、その共分散行列の逆行列のｍｍ’の要素を示す
第１の量（Ａmm’）と、ｍ番目の検出器出力の測定値と
生体内の電流分布推定値から算出されたｍ番目の検出器
出力の差を示す第２の量（Ｄm−ＳＢm）と、ｍ’番目の
検出器出力の測定値と生体内の電流分布推定値から算出
されたｍ’番目の検出コイル出力の差を示す第３の量
（Ｄm’−ＳＢm’）との積を、ｍｍ’の全組み合わせに
ついて積分した量を含む対数尤度Ｅを求めて、そのＥを
最大とする磁場源パラメータ（ＳＲn,ＳＱｎ）により、
生体内の電流分布を推定することに特徴がある。また、
上記ノイズの共分散行列を推定する場合、心磁波形や自
発脳磁等の計測では、測定対象の生体磁場を検出しない
位置に磁束計を設置して計測を行い、誘発脳磁計測で
は、誘発刺激の付与以前から連続してデータの取り込み
を行うことに特徴があり、さらに、心磁波形では、心拍
の各周期の間にある信号の極めて微弱な部分からノイズ
の共分散行列を推定することに特徴がある。
【０００５】
【作用】本発明においては、最尤推定法で、対数尤度と
呼ばれる量を最大とすることにより、最適推定値を求め
る。一般に、ノイズのチャネル間での相関を考慮した場
合の対数尤度Ｅは、
呼ばれる量を最大とすることにより、最適推定値を求め
る。一般に、ノイズのチャネル間での相関を考慮した場
合の対数尤度Ｅは、
【数３】
で与えられることが知られている。ここで、Ａmm'は共
分散行列Ｃの逆行列の第（ｍ，ｍ’）要素である。ま
た、最適推定値はＥを最大とする推定値として求められ
る。この際、チャネル間の相関が無視できる場合には、 Ａmm'＝Ｖ（ｍ）（ｍ＝ｍ’），Ａmm'＝０（ｍ≠ｍ’） となる。但し、Ｖ（ｍ）は第ｍチャネル出力に混入して
いるノイズの分散である。さらに、Ｖ（ｍ）が全てのチ
ャネルで等しい値Ｖである場合、式（３）は、結局、
分散行列Ｃの逆行列の第（ｍ，ｍ’）要素である。ま
た、最適推定値はＥを最大とする推定値として求められ
る。この際、チャネル間の相関が無視できる場合には、 Ａmm'＝Ｖ（ｍ）（ｍ＝ｍ’），Ａmm'＝０（ｍ≠ｍ’） となる。但し、Ｖ（ｍ）は第ｍチャネル出力に混入して
いるノイズの分散である。さらに、Ｖ（ｍ）が全てのチ
ャネルで等しい値Ｖである場合、式（３）は、結局、
【数４】
となり、式（２）に示した最小２乗法と等しくなる。従
って、複数のチャネルを持つ検出器を用いた測定データ
から脳内の磁場源を推定する場合、ノイズの共分散行列
Ｃを測定し、式（３）を用いて、Ｅを最大とする磁場源
パラメータを求めることにより、ノイズの相関を考慮し
た磁場源の推定を行うことができ、精度の高い電流分布
推定を可能にする。
って、複数のチャネルを持つ検出器を用いた測定データ
から脳内の磁場源を推定する場合、ノイズの共分散行列
Ｃを測定し、式（３）を用いて、Ｅを最大とする磁場源
パラメータを求めることにより、ノイズの相関を考慮し
た磁場源の推定を行うことができ、精度の高い電流分布
推定を可能にする。
【０００６】
【実施例】以下、本発明の一実施例を図面により説明す
る。本実施例では、生体磁気計測で代表的な脳磁計測に
ついて述べる。図３は、本発明の一実施例における多チ
ャネル生体磁気計測システムの構成図である。図３にお
いて、３１は被験者、３２は液体ヘリウムの満たされた
デュワー、３３はＳＱＵＩＤ(Superconducting Quautum
Interference Deuise)素子に接続された検出コイル、
３４はＦＬＬ回路、３５は磁気共鳴像、３６はＣＲＴデ
ィスプレイ、３７はＣＰＵ、メモリ等から構成された計
算機である。特に、計算機３７には、図１のステップ１
０２〜１０４に示す処理機能を設ける。このような構成
によって、被検者２１からの信号は複数個の検出コイル
２３で同時に検出され、ＳＱＵＩＤ素子およびＦＬＬ回
路２４を介して、計算機２７に入力される。ここで、磁
場源を推定した後、ＣＲＴディスプレイ２６に推定磁場
源を表示する。なお、表示はＸ線像あるいは磁気共鳴像
２５と重ね合わせて行われる。
る。本実施例では、生体磁気計測で代表的な脳磁計測に
ついて述べる。図３は、本発明の一実施例における多チ
ャネル生体磁気計測システムの構成図である。図３にお
いて、３１は被験者、３２は液体ヘリウムの満たされた
デュワー、３３はＳＱＵＩＤ(Superconducting Quautum
Interference Deuise)素子に接続された検出コイル、
３４はＦＬＬ回路、３５は磁気共鳴像、３６はＣＲＴデ
ィスプレイ、３７はＣＰＵ、メモリ等から構成された計
算機である。特に、計算機３７には、図１のステップ１
０２〜１０４に示す処理機能を設ける。このような構成
によって、被検者２１からの信号は複数個の検出コイル
２３で同時に検出され、ＳＱＵＩＤ素子およびＦＬＬ回
路２４を介して、計算機２７に入力される。ここで、磁
場源を推定した後、ＣＲＴディスプレイ２６に推定磁場
源を表示する。なお、表示はＸ線像あるいは磁気共鳴像
２５と重ね合わせて行われる。
【０００７】次に、計算機３７による磁場源の推定方法
について述べる。図１は、本発明の一実施例における磁
場源の推定方法を示すフローチャート、図４は誘発脳磁
の波形を模式的に示した説明図、図５は心磁波形を模式
的に示した説明図である。本実施例では、図１のよう
に、誘発脳磁界において、多チャネル検出器による測定
により、聴視覚等の感覚刺激付与以前から複数回の測定
データ（外部磁場源からのノイズのみを含む測定デー
タ、および計測対象の信号を含む測定データ）を得て計
算機３７に入力する（１０１）。計算機３７では、ま
ず、共分散行列Ｃを推定する（１０２）。ここで、外部
磁場源によるノイズを含むｉ番目のチャネルからの出力
をＮi、ｊ番目のチャネルからの出力をＮjとすると、共
分散行列Ｃの（ｉ，ｊ）成分ＣijはＣij＝〈ＮiＮj〉で
与えられる。但し、記号〈 〉は集合平均を意味する。
従って、Ｃijの推定は、多チャネル検出器を用いてノイ
ズのみの測定を多数回行い、刺激付与以前のデータから
外部磁場源からのノイズを含む出力の積ＮiＮjを計算し
て、その結果を加算平均することにより行う。このよう
に、誘発脳磁界の刺激付与以前の測定データは信号を含
まないので、ｉチャネルの時刻ｔにおける測定値をＮi
（ｔ）とし、刺激付与以前の測定点の時刻をｔ1，ｔ2，
…，ｔLとすると、Ｃijは、
について述べる。図１は、本発明の一実施例における磁
場源の推定方法を示すフローチャート、図４は誘発脳磁
の波形を模式的に示した説明図、図５は心磁波形を模式
的に示した説明図である。本実施例では、図１のよう
に、誘発脳磁界において、多チャネル検出器による測定
により、聴視覚等の感覚刺激付与以前から複数回の測定
データ（外部磁場源からのノイズのみを含む測定デー
タ、および計測対象の信号を含む測定データ）を得て計
算機３７に入力する（１０１）。計算機３７では、ま
ず、共分散行列Ｃを推定する（１０２）。ここで、外部
磁場源によるノイズを含むｉ番目のチャネルからの出力
をＮi、ｊ番目のチャネルからの出力をＮjとすると、共
分散行列Ｃの（ｉ，ｊ）成分ＣijはＣij＝〈ＮiＮj〉で
与えられる。但し、記号〈 〉は集合平均を意味する。
従って、Ｃijの推定は、多チャネル検出器を用いてノイ
ズのみの測定を多数回行い、刺激付与以前のデータから
外部磁場源からのノイズを含む出力の積ＮiＮjを計算し
て、その結果を加算平均することにより行う。このよう
に、誘発脳磁界の刺激付与以前の測定データは信号を含
まないので、ｉチャネルの時刻ｔにおける測定値をＮi
（ｔ）とし、刺激付与以前の測定点の時刻をｔ1，ｔ2，
…，ｔLとすると、Ｃijは、
【数５】
から近似的に計算できる。また、通常の脳磁界は非常に
微弱であるため、磁場源推定に用いるデータとして測定
データを多数回加算平均したものを用いることが一般に
行われる。この場合は、式（５）におけるＮi（ｔ）Ｎj
（ｔ）としては、磁場源推定に用いられるデータと同じ
回数程度、加算平均を行ったものを用いる。一方、心磁
波形あるいは自発脳磁等、刺激を用いない信号の場合に
は、誘発脳磁と同じ方法を用いることはできない。この
場合、例えば、磁束計を信号源である生体から離して、
まず、共分散行列Ｃを推定するためのデータを測定し、
次に生体に十分近づけて、磁場源推定用のデータを測定
する方法を用いる。なお、共分散行列の推定には、式
（５）と同様の演算を行う。また、心磁波形等では、図
５のＡで示すように、心拍の一同期と次の同期の間に信
号の極めて弱い部分が存在し、この部分からノイズの共
分散行列を推定できる。このＡで示した部分における測
定データから、上記Ｎi（ｔ）Ｎj（ｔ）を計算し、この
部分におけるｔに関しての加算平均をとることにより、
ノイズの共分散行列を推定できる。次に、こうして推定
された共分散行列から、例えば、掃き出し法を用いてそ
の逆行列の要素Ａijを求める（１０３）。この掃き出し
法については、例えば「伊東由文著，線形代数学，共立
出版」において詳述されている。なお、共分散行列の固
有値に他と比べて非常に小さなものを含む場合は、その
固有値の僅かな大きさの違いが逆行列の計算結果に大き
く影響することがある。この場合、他と比べて小さな固
有値をゼロとして偽似的な逆行列を求め、その（ｉ，
ｊ）要素をＡijとして式（３）を計算する。さて、信号
を含む測定値を得た後、これら測定値から磁場発生源の
位置および電流ベクトル値等のパラメータを推定する
（１０４）。すなわち、図２に示したＲmの位置で実際
に測定される測定値をＤmとし、第ｎ番目のダイポール
の位置および電流ベクトルの推定値をＳＲn，ＳＱnとし
て、これらから式（１）を用いて計算された位置Ｒmに
おける磁場の法線成分をＳＢmとすると、ＳＲn，ＳＱn
の最適推定値は、式（３）によって表わされたＥを最小
とするＳＲn，ＳＱn（ｎ＝１，２，…，Ｎ）を求めるこ
とにより得られる。なお、式（３）によって表わされた
Ｅを最大とするパラメータＳＲn，ＳＱnの値を求めるに
は、例えば、シミュレーテッド・アニーリングを用い
る。このアルゴリズムについては、例えば「ピー・ジェ
ー・エム・バン ラホーベン、イー・エイチ・エル・アーツ
著、シミュレーテッドアニーリング：セオリー アンド
アプリケーション、デー・レーデル パブリッシング
カンパニー(P.J.M.van Laarhoven,E.H.L.Aarts、Simu
lated Annealing:Thory and Applications、D.Reidel P
ublishing Company)」や「関原謙介、大山永昭著、生体
磁場逆問題のためのシミュレーテッド アニーリング
アルゴリズム、電気学会研究会資料、ＭＡＧ−９０−１
３８」等に詳述されており、既に多くの提案がなされて
いる。ここで、本実施例（処理１０４）に用いるシミュ
レーテッド・アニーリングの処理の流れを示す。図６は
本発明の一実施例におけるシミュレーテッド・アニーリ
ングを示すフローチャートである。図６において、ｆは
ベクトルであって推定値の組を表わし、ｆ’はｆと僅か
に異なる推定値である。また、Ｔは仮想的な温度であっ
て関数Ｅの変化ΔＥが０より大きい場合でも推定値の変
化を受け入れる確率を決めるためのパラメータである。
また、式（２）で用いたＳＲn、ＳＱn（ｎ＝１，２，
…，Ｎ）はｆの各要素に対応する。本実施例のアルゴリ
ズムにおいては、まずｆに対して初期値をセットし（６
０１）、温度Ｔに対して十分大きい初期値をセットし
（６０２）、一連の処理（６０６〜６０９）の試行回数
のカウントＮを０としておく（６０３）。ここで、現在
の推定値ｆと僅かに異なる推定値をｆ’とし（６０
６）、ｆのｆ’の変化に対する関数Ｅの変化ΔＥを計算
する（６０７）。そして、確率Ｐ（ΔＥ）に従い、変化
ｆ’を受け入れるか否かを決める（６０８，６０９）。
なお、受け入れるとは、推定値ｆをｆ’で置き換えるこ
とを言う。また、Ｐ（ΔＥ）は次の式で計算する。
微弱であるため、磁場源推定に用いるデータとして測定
データを多数回加算平均したものを用いることが一般に
行われる。この場合は、式（５）におけるＮi（ｔ）Ｎj
（ｔ）としては、磁場源推定に用いられるデータと同じ
回数程度、加算平均を行ったものを用いる。一方、心磁
波形あるいは自発脳磁等、刺激を用いない信号の場合に
は、誘発脳磁と同じ方法を用いることはできない。この
場合、例えば、磁束計を信号源である生体から離して、
まず、共分散行列Ｃを推定するためのデータを測定し、
次に生体に十分近づけて、磁場源推定用のデータを測定
する方法を用いる。なお、共分散行列の推定には、式
（５）と同様の演算を行う。また、心磁波形等では、図
５のＡで示すように、心拍の一同期と次の同期の間に信
号の極めて弱い部分が存在し、この部分からノイズの共
分散行列を推定できる。このＡで示した部分における測
定データから、上記Ｎi（ｔ）Ｎj（ｔ）を計算し、この
部分におけるｔに関しての加算平均をとることにより、
ノイズの共分散行列を推定できる。次に、こうして推定
された共分散行列から、例えば、掃き出し法を用いてそ
の逆行列の要素Ａijを求める（１０３）。この掃き出し
法については、例えば「伊東由文著，線形代数学，共立
出版」において詳述されている。なお、共分散行列の固
有値に他と比べて非常に小さなものを含む場合は、その
固有値の僅かな大きさの違いが逆行列の計算結果に大き
く影響することがある。この場合、他と比べて小さな固
有値をゼロとして偽似的な逆行列を求め、その（ｉ，
ｊ）要素をＡijとして式（３）を計算する。さて、信号
を含む測定値を得た後、これら測定値から磁場発生源の
位置および電流ベクトル値等のパラメータを推定する
（１０４）。すなわち、図２に示したＲmの位置で実際
に測定される測定値をＤmとし、第ｎ番目のダイポール
の位置および電流ベクトルの推定値をＳＲn，ＳＱnとし
て、これらから式（１）を用いて計算された位置Ｒmに
おける磁場の法線成分をＳＢmとすると、ＳＲn，ＳＱn
の最適推定値は、式（３）によって表わされたＥを最小
とするＳＲn，ＳＱn（ｎ＝１，２，…，Ｎ）を求めるこ
とにより得られる。なお、式（３）によって表わされた
Ｅを最大とするパラメータＳＲn，ＳＱnの値を求めるに
は、例えば、シミュレーテッド・アニーリングを用い
る。このアルゴリズムについては、例えば「ピー・ジェ
ー・エム・バン ラホーベン、イー・エイチ・エル・アーツ
著、シミュレーテッドアニーリング：セオリー アンド
アプリケーション、デー・レーデル パブリッシング
カンパニー(P.J.M.van Laarhoven,E.H.L.Aarts、Simu
lated Annealing:Thory and Applications、D.Reidel P
ublishing Company)」や「関原謙介、大山永昭著、生体
磁場逆問題のためのシミュレーテッド アニーリング
アルゴリズム、電気学会研究会資料、ＭＡＧ−９０−１
３８」等に詳述されており、既に多くの提案がなされて
いる。ここで、本実施例（処理１０４）に用いるシミュ
レーテッド・アニーリングの処理の流れを示す。図６は
本発明の一実施例におけるシミュレーテッド・アニーリ
ングを示すフローチャートである。図６において、ｆは
ベクトルであって推定値の組を表わし、ｆ’はｆと僅か
に異なる推定値である。また、Ｔは仮想的な温度であっ
て関数Ｅの変化ΔＥが０より大きい場合でも推定値の変
化を受け入れる確率を決めるためのパラメータである。
また、式（２）で用いたＳＲn、ＳＱn（ｎ＝１，２，
…，Ｎ）はｆの各要素に対応する。本実施例のアルゴリ
ズムにおいては、まずｆに対して初期値をセットし（６
０１）、温度Ｔに対して十分大きい初期値をセットし
（６０２）、一連の処理（６０６〜６０９）の試行回数
のカウントＮを０としておく（６０３）。ここで、現在
の推定値ｆと僅かに異なる推定値をｆ’とし（６０
６）、ｆのｆ’の変化に対する関数Ｅの変化ΔＥを計算
する（６０７）。そして、確率Ｐ（ΔＥ）に従い、変化
ｆ’を受け入れるか否かを決める（６０８，６０９）。
なお、受け入れるとは、推定値ｆをｆ’で置き換えるこ
とを言う。また、Ｐ（ΔＥ）は次の式で計算する。
【数６】
上記のステップ６０６〜６０９をＮmax回繰返し、平行
状態に達しているか否かを判断する（６１０）。つま
り、ΔＥ＞０である試行の数をＮ1、ΔＥ＜０で受け入
れられる試行の数をＮ2として、｜Ｎ1−Ｎ2｜／｜Ｎ1＋
Ｎ2｜がある値εより小さいか否かで判断する。この場
合、例えばＮmax＝２００、ε＝０.０１の値にしてお
く。さらに、ステップ６１０で、平衡状態にあると判定
され、Ｎ1＋Ｎ2≠０ならば（６１１）、温度を僅かに下
げて（６１２）、同様の手順を繰り返す。ここで、ステ
ップ６１１は、アルゴリズムの終了の判断を示し、アル
ゴリズムの終了は受け入れられる試行の数（Ｎ1＋Ｎ2）
がゼロとなる場合とする。また、温度の下げ方として
は、
状態に達しているか否かを判断する（６１０）。つま
り、ΔＥ＞０である試行の数をＮ1、ΔＥ＜０で受け入
れられる試行の数をＮ2として、｜Ｎ1−Ｎ2｜／｜Ｎ1＋
Ｎ2｜がある値εより小さいか否かで判断する。この場
合、例えばＮmax＝２００、ε＝０.０１の値にしてお
く。さらに、ステップ６１０で、平衡状態にあると判定
され、Ｎ1＋Ｎ2≠０ならば（６１１）、温度を僅かに下
げて（６１２）、同様の手順を繰り返す。ここで、ステ
ップ６１１は、アルゴリズムの終了の判断を示し、アル
ゴリズムの終了は受け入れられる試行の数（Ｎ1＋Ｎ2）
がゼロとなる場合とする。また、温度の下げ方として
は、
【数７】
あるいは
【数８】
等が知られている。但し、ａは０.９〜０.９９程度の値
である。また、Ｔkはｋ番目の温度を意味する。なお、
本実施例では、式（３）によって表わされたＥを最大と
する電流ダイポールパラメータの値を求めるのに、シミ
ュレーテッド・アニーリングを用いたが、これに限られ
るものではなく、例えば、ジェネティック・アルゴリズ
ム等を用いることもできる。これについては、例えば
「ゴールドバーク著、ジェネティック アルゴリズム
イン サーチ，オプチミゼーション，アンド マシン
ラーニング(David E.Goldberg,Genetic Algorithms in
Search,Optimization,and Machine Learning)」におい
て述べられている。また、生体磁場発生に関与している
電流ダイポールが１個である場合には、関数Ｅに局所的
な最大値を含まないので、Ｅの最大値を求めるのにシミ
ュレーテッド アニーリングやジェネティック アルゴ
リズム等を用いる必要はなく、通常よく知られているア
ルゴリズム、例えば最急上昇法、共役勾配法、ニュート
ン法等を用いることができる。これらについては、例え
ば「エイチ・エム・ワグナー著、プリンシプル オブ オ
ペレーション リサーチ、プレンティス−ホール イン
ク イングルウッド クリフ ニュージャージ(H.M.Wag
ner,Principles of Operations Research,Prentice-Hal
l,Inc.,Englewood Cliffs,New Jersey)」において述べ
られている。さらに、本実施例では、生体内の関心領域
に含まれる磁場源の数は予め知られているものとしてい
るが、磁場源の正確な個数が分からない場合に拡張する
こともできる。すなわち、特願平２−７６８８３号およ
び特願平２−１１００１３号で提案されている方法を利
用する。これは、磁場源の正確な個数がわからない場合
でも磁場源パラメータを推定できる方法であって、推定
値と測定値との一致度を表わす２乗誤差項に付加的な項
を加え、その和を最小とする推定値を最適推定値として
求めるものである。この２乗誤差項を式（３）で表わし
たＥの符号を逆にした−Ｅで置き換えることにより、本
実施例を磁場源の正確な個数がわからない場合に拡張で
きる。こうして、式（３）によって表わされたＥを最大
とする磁場源パラメータＳＲn，ＳＱnの値を求めると、
図３に示したようにＣＲＴディスプレイ３６に表示する
（１０５）。
である。また、Ｔkはｋ番目の温度を意味する。なお、
本実施例では、式（３）によって表わされたＥを最大と
する電流ダイポールパラメータの値を求めるのに、シミ
ュレーテッド・アニーリングを用いたが、これに限られ
るものではなく、例えば、ジェネティック・アルゴリズ
ム等を用いることもできる。これについては、例えば
「ゴールドバーク著、ジェネティック アルゴリズム
イン サーチ，オプチミゼーション，アンド マシン
ラーニング(David E.Goldberg,Genetic Algorithms in
Search,Optimization,and Machine Learning)」におい
て述べられている。また、生体磁場発生に関与している
電流ダイポールが１個である場合には、関数Ｅに局所的
な最大値を含まないので、Ｅの最大値を求めるのにシミ
ュレーテッド アニーリングやジェネティック アルゴ
リズム等を用いる必要はなく、通常よく知られているア
ルゴリズム、例えば最急上昇法、共役勾配法、ニュート
ン法等を用いることができる。これらについては、例え
ば「エイチ・エム・ワグナー著、プリンシプル オブ オ
ペレーション リサーチ、プレンティス−ホール イン
ク イングルウッド クリフ ニュージャージ(H.M.Wag
ner,Principles of Operations Research,Prentice-Hal
l,Inc.,Englewood Cliffs,New Jersey)」において述べ
られている。さらに、本実施例では、生体内の関心領域
に含まれる磁場源の数は予め知られているものとしてい
るが、磁場源の正確な個数が分からない場合に拡張する
こともできる。すなわち、特願平２−７６８８３号およ
び特願平２−１１００１３号で提案されている方法を利
用する。これは、磁場源の正確な個数がわからない場合
でも磁場源パラメータを推定できる方法であって、推定
値と測定値との一致度を表わす２乗誤差項に付加的な項
を加え、その和を最小とする推定値を最適推定値として
求めるものである。この２乗誤差項を式（３）で表わし
たＥの符号を逆にした−Ｅで置き換えることにより、本
実施例を磁場源の正確な個数がわからない場合に拡張で
きる。こうして、式（３）によって表わされたＥを最大
とする磁場源パラメータＳＲn，ＳＱnの値を求めると、
図３に示したようにＣＲＴディスプレイ３６に表示する
（１０５）。
【０００８】
【発明の効果】本発明によれば、複数の検出コイルを持
つ磁束計で生体磁場を計測し、その結果から生体内の磁
場発生源を推定する場合、検出コイル間でのノイズの相
関を測定し、これを磁場源推定の演算に組み入れること
により、推定誤差を低減することができる。
つ磁束計で生体磁場を計測し、その結果から生体内の磁
場発生源を推定する場合、検出コイル間でのノイズの相
関を測定し、これを磁場源推定の演算に組み入れること
により、推定誤差を低減することができる。
【０００９】
【図１】本発明の一実施例における磁場源の推定方法を
示すフローチャートである。
示すフローチャートである。
【図２】従来の脳磁計測における磁場発生源の推定方法
の説明図である。
の説明図である。
【図３】本発明の一実施例における多チャネル生体磁気
計測システムの構成図である。
計測システムの構成図である。
【図４】誘発脳磁の波形を模式的に示した説明図であ
る。
る。
【図５】心磁波形を模式的に示した説明図である。
【図６】本発明の一実施例におけるシミュレーテッド・
アニーリングを示すフローチャートである。
アニーリングを示すフローチャートである。
２１ 測定位置
３１ 被験者
３２ デュワー
３３ 検出コイル
３４ ＦＬＬ回路
３５ 磁気共鳴像
３６ ＣＲＴディスプレイ
３７ 計算機
Claims (7)
- 【請求項１】 生体の活動に伴って発生する磁場分布
を、複数の磁束計により生体外部の複数個の点で同時に
測定し、測定結果から生体内部の電流分布を推定する方
法において、上記磁束計の検出コイル出力間でのノイズ
の共分散行列を推定し、該共分散行列から求まる量を含
む測定磁場分布値の２次形式で表わされる量を含む関数
を求めて、該関数を最大か最小とする磁場源パラメータ
により、上記電流分布を推定することを特徴とする磁場
発生源の推定方法。 - 【請求項２】 上記検出コイル出力間でのノイズの共分
散行列から求まる量として、共分散行列の逆行列の要素
を用いることを特徴とする請求項１記載の磁場発生源の
推定方法。 - 【請求項３】 上記検出コイルの番号をｍ、ｍ’とした
場合、上記２次形式で表わされる量として、共分散行列
の逆行列のｍｍ’の要素を示す第１の量と、ｍ番目の検
出器出力の測定値と生体内の電流分布推定値から算出さ
れたｍ番目の検出器出力の差を示す第２の量と、ｍ’番
目の検出器出力の測定値と生体内の電流分布推定値から
算出されたｍ’番目の検出コイル出力の差を示す第３の
量との積を、ｍｍ’の全組み合わせについて積分した量
を用いることを特徴とする請求項２記載の磁場発生源の
推定方法。 - 【請求項４】 上記ノイズの共分散行列を、誘発刺激を
用いない信号による測定結果から推定する場合、測定対
象の生体磁場を検出しない位置に磁束計を設置して計測
を行い、該計測の結果からノイズの共分散行列を推定す
ることを特徴とする請求項１記載の磁場発生源の推定方
法。 - 【請求項５】 上記ノイズの共分散行列を、誘発脳磁計
測による測定結果から推定する場合、誘発刺激の付与以
前から連続してデータの取り込みを行い、付与以前のデ
ータからノイズの共分散行列を推定することを特徴とす
る請求項１記載の磁場発生源の推定方法。 - 【請求項６】 上記ノイズの共分散行列を、心磁波形を
含む周期的な生体磁気信号の測定結果から推定する場
合、該信号の微弱部分からノイズの共分散行列を推定す
ることを特徴とする請求項１記載の磁場発生源の推定方
法。 - 【請求項７】 上記共分散行列の（ｍ，ｍ’）要素を推
定する際、ｍ番目の検出コイル出力とｍ’番目の検出コ
イル出力の積を時間平均した量を用いることを特徴とす
る請求項１〜６記載の磁場発生源の推定方法。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| JP3167883A JPH0515504A (ja) | 1991-07-09 | 1991-07-09 | 磁場発生源の推定方法 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| JP3167883A JPH0515504A (ja) | 1991-07-09 | 1991-07-09 | 磁場発生源の推定方法 |
Publications (1)
| Publication Number | Publication Date |
|---|---|
| JPH0515504A true JPH0515504A (ja) | 1993-01-26 |
Family
ID=15857849
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| JP3167883A Pending JPH0515504A (ja) | 1991-07-09 | 1991-07-09 | 磁場発生源の推定方法 |
Country Status (1)
| Country | Link |
|---|---|
| JP (1) | JPH0515504A (ja) |
Cited By (3)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US6370414B1 (en) | 1998-01-23 | 2002-04-09 | Ctf Systems, Inc. | System and method for measuring, estimating and displaying RMS current density maps |
| US6697660B1 (en) * | 1998-01-23 | 2004-02-24 | Ctf Systems, Inc. | Method for functional brain imaging from magnetoencephalographic data by estimation of source signal-to-noise ratio |
| US10575738B2 (en) | 2016-07-13 | 2020-03-03 | Advantest Coporation | Magnetic field measurement apparatus and method for noise environment |
-
1991
- 1991-07-09 JP JP3167883A patent/JPH0515504A/ja active Pending
Cited By (3)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US6370414B1 (en) | 1998-01-23 | 2002-04-09 | Ctf Systems, Inc. | System and method for measuring, estimating and displaying RMS current density maps |
| US6697660B1 (en) * | 1998-01-23 | 2004-02-24 | Ctf Systems, Inc. | Method for functional brain imaging from magnetoencephalographic data by estimation of source signal-to-noise ratio |
| US10575738B2 (en) | 2016-07-13 | 2020-03-03 | Advantest Coporation | Magnetic field measurement apparatus and method for noise environment |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| US5269325A (en) | Analysis of biological signals using data from arrays of sensors | |
| US6230037B1 (en) | Biomagnetic field measuring method and apparatus | |
| US8838225B2 (en) | Analysis of multi-channel measurement data using orthogonal virtual channels | |
| US5526811A (en) | Apparatus and process for determining the sources of biomagnetic activity | |
| US4977896A (en) | Analysis of biological signals using data from arrays of sensors | |
| US6370414B1 (en) | System and method for measuring, estimating and displaying RMS current density maps | |
| Auranen et al. | Bayesian analysis of the neuromagnetic inverse problem with ℓp-norm priors | |
| US5887074A (en) | Local principal component based method for detecting activation signals in functional MR images | |
| Long et al. | State-space solutions to the dynamic magnetoencephalography inverse problem using high performance computing | |
| US5687724A (en) | Apparatus and method for determining variations in measured physical parameters of signal-generators | |
| JP6890484B2 (ja) | 磁界計測装置および計測磁界表示方法 | |
| Mohsen et al. | AI aided noise processing of spintronic based IoT sensor for magnetocardiography application | |
| Sekihara et al. | Maximum-likelihood estimation of current-dipole parameters for data obtained using multichannel magnetometer | |
| Ioannides | Estimates of Brain Activity. Using Magnetic Field Tomography and | |
| US10429479B2 (en) | Rapid measurement of perfusion using optimized magnetic resonance fingerprinting | |
| JPH0515504A (ja) | 磁場発生源の推定方法 | |
| Singh et al. | Neuromagnetic localization using magnetic resonance images | |
| JP3519141B2 (ja) | 生体内電流源推定方法 | |
| Nagarajan et al. | Magnetoencephalographic imaging | |
| JP2752884B2 (ja) | 生体活動電流源推定方法 | |
| Ogata et al. | Study of spatial filter for magnetocardiography measurements without a magnetically shielded room | |
| JP3227958B2 (ja) | 生体活動電流源推定方法 | |
| JP4006543B2 (ja) | 生体活動電流源推定装置 | |
| JP7002416B2 (ja) | 磁界計測装置 | |
| JP3409551B2 (ja) | 生体磁気計測装置 |