JPH03277345A - 生体磁場計測における磁場発生源の推定方法 - Google Patents

Abstract

(57)【要約】本公報は電子出願前の出願データであるた
め要約のデータは記録されません。

Description

【発明の詳細な説明】

〔産業上の利用分野〕 本発明は生体表面の磁場分布を測定し、そのデータから
生体内部の磁場発生源の分布を推定し表示する生体磁場
計測装置における磁場発生源の推定方法に関する。 〔従来の技術〕 生体外部の磁場測定値から生体内部の電流分布の推定は
通常以下のように行なわれる。 以下、ここでは生体磁場計測でも代表的な脳磁計測を例
として説明する。 座標系第１図にように定義する。図ではｒｌの位置にベ
クトルｑｎで表わされた電流が存在すると仮定する。脳
磁気計測の分野では脳の境界を球面と仮定できる限りこ
のような孤立した電流ベクトルを仮定でき、これを電流
ダイポールと呼ぶ。 また、第１図でｒｌは測定点の座標である。 さて、脳内にＮ個の電流ダイポールを仮定すると、ｒｌ
Ｉの点での磁束密度ベクトルＢ、はで表わされる。 通常、磁束計で測定できるのはＢ１の法線成分である。 さて、ｒ、の位置で実際に得られる測定値をり、とする
。 また、各ダイポールの位置および電流ベクトルの推定値
を°を付けて表わすと推定値から計算される仮想的な測
定値り、は である。 ここで、 コスト関数を Ｅ（ｒｌｌ　　ｒｚ＋　　−’゛ｒｓ；　　ｑ１ｑｚ＋
　　”’ＱＮ）と定義する。（３）式は推定データと実
際の測定データとの一致度を表わす。 従来、最適推定値は上の（３）式で定義されたコスト関
数を最小にするｒ□、・・・、　ｒＮｌ　ｑｔ＋・・Ｑ
Ｎとして求められた。 〔発明が解決しようとする課題〕 以上の事かられかる通り、従来の推定手法においては、
生体内の磁場発生源、つまりは電流ダイポールの個数Ｎ
は既知でなければならなかった。 もし、実際の電流ダイポールの数と異なる数Ｎを設定し
て上記コスト関数を最小にするｒｉｌ・・生体の電流分
布と異なるものとなってしまう。 そこで、本発明のひとつの目的は磁場発生源の数、つま
り対象とする生体内の電流ダイポールの数が正確に知ら
れていない場合でも正しくその生体内の電流分布を示す
解を得ることができる磁場発生源の推定方法を提供する
にある。 本発明の別の目的は実際に存在する数より多い数の電流
ダイポーは、仮定されて逆問題法による推定値の特定が
成されたとき、実際に存在する数の電流ダイポールを示
す解が自動的に得られる磁場発生源の推定方法を提供す
るにある。 〔課題を解決するための手段〕 °本発明は、対象とする生体の外部の複数位置にてこの
生体から発生する磁場を計測し５一方生体内に磁場発生
源として複数の電流ダイポールを仮定してその座標の推
定値と電流ベクトルの推定値から上記複数位置における
仮想的な磁場の計測値を算出し、実際の磁場の計測値の
分布と上記の算出された仮想的な磁場の計測値の分布と
の差を基本的に示すコスト関数が最小となる上記の推定
値を求めて上記の生体内の電流分布を推定する方法にお
いて、予想される生体内の電流ダイポールの数より多い
数の電流ダイポールを仮定して測定点の位置及び電流ベ
クトルの推定値から仮想的な磁場計測値を算出すること
、及びコスト関数として上述の実際の磁場計測値の分布
と仮想的な磁場計測値の分布との差を示す量に加えて仮
定した測定点の電流総和にかかわる量、 代表的には各電流ベクトルの絶対値の総和、もしくはそ
のα乗（α法線成分の２乗の総和を含むものを用いる点
を特徴とする。 本発明の別の特徴は実施例の説明に示される。 〔作用〕 上記のように、予想されるより多い数の電流ダイポール
の位置と電流ベクトルの推定値から算出した磁場計測値
と実測された磁場計測値との一致度が向上するよう推定
値を変位させて行くと、推定値に過剰な自由度が導入さ
れてしまうので、実際の電流分布と異なる解が得られる
可能性が高い。 これに対し、仮定した測定点の電流ベクトルの絶対値の
総和、もしくはそのα乗の総和は、推定値の過剰な自由
度をおさえる働きを有するので、正しい解が得られる。 つまりは生体内の電流ダイポール（独立した磁場発生源
）の数が正確に分っていなくても正しい電流分布が推定
できる。 る。 ただし、一般的に言えば得られる解は、はぼ同じ位置に
重なる複数の電流ベクトルを含むものとなる。このほぼ
同じ位置に重なる複数の電流ベクトルの合成ベクトルは
、実際に生体内に含まれるひとつの電流ダイポールを示
している。もし０くα〈ｌなるαを用い、コスト関数に
付加される量として仮定された測定点の電流ベクトルの
絶対値のα乗の総和を採用すれば、このコスト関数を最
小とする解は、このようなほぼ同一位置に重なる複数の
電流ベクトルを含まない。つまりこのような複数の電流
ベクトルのうちひとつを除いて残りはベクトル量がゼロ
となるので、自動的に実際に生体に含まれる数の電流ダ
イポーは、推定できる。 〔実施例〕 脳磁計測を例として本発明の詳細な説明する。 脳磁測定は第１図に示されるごとく、磁束計１０を点１
２−１から１２−Ｍのように１点１点動かしながら頭部
のまわりの測定点における磁場の洗練成分を測定して行
く。さて、このようにして測定された頭部周囲の測定点
がら磁場発生源となっている脳内の電流ダイポールの空
間座標と電流ベクトルの推定を本発明を用いて以下の通
りに行なうことができる。 頭部表面の測定点の個数をＭとし、１からＭまで番号付
けを行なう、第ｍ番目の測定点の位置をベクトルｒ、で
表わし、その点における磁場の測定値をり、と表記する
。 電流ダイポーは、脳の内部にＮ個含まれていると仮定す
る１本発明はこのダイ′ボールの正確な個数が既知でな
い場合でも適用できるものであり。 ひとまずＮは脳内部に存在すると予想されるダイポール
の個数より多めに設定しておく。 イ反 このように回定した電流ダイポールの推定位置をｒｌｌ
、推定電流ベクトルをｑｌｌで表わす。ここでｎ＝１．
２．・・・、Ｎである。このように仮定された電流ダイ
ポーは、前記の計測点に作る磁場り。 （ｍ＝１．・・・、Ｍ）を計算により求める。つまりＤ
ヨは推定した電流ダイポールから求めた仮想的な磁場計
測値と呼ぶことができる。頭部表面を球面で近似し、こ
の球の中心を座標原点とすれば。 Ｄ、は次式で計算される。 次に後で詳述する最適化演算により仮定した電流ダイポ
ールの位置推定値’Ｉ、＋電流ベクトル推定値Ｑｎを変
位させて行き、コスト関数が最小となる推定値を求める
。このコスト関数として、本実施例では次式のものを用
いる。 ｍ＝１ （４） 二こで、（４）式右辺第１項は電流ダイポールの推定値
から算出した仮想的な磁場計測値の分布と実際の計測値
データとの一致度を表わす項であり、第２項は推定され
た電流分布の総和に関係した量である。すなわち、本発
明では、ダイポール数を余分に仮定しておき、実際のデ
ータの一致度を最良にするもののうち、さらに電流分布
の総和を最小にするｒ　ｔ＋ｇ　Ｑｎ（ｎ　＝　１　？
・・・、Ｎ）を選ぶ。 では とするのが自然である。ここでＷ、は定数である。 ところで（５）式のごとく関数Ｅｓを定めると場合によ
っては真のダイポール位置よりもより脳表面に近い位置
に誤まった推定が行なわれることがある。これを防ぐ別
の例では、 ｎ＝１ とする。 ここで であり、 さらにｇ′″□は である。すなわちベクトル９Ｉｌの各成分ｇ″″□（ｍ
＝１．・・Ｍ）はｎ番目の電流ダイポーは、ｍ番目の測
定点位置につくる磁場の法線成分（検出コイル面と垂直
な成分）である。したがって（６）式におけるベクトル
９ｒｌの絶対値１ｓｙｒｌＩは、ベクトルの各成分の、
つまりｎ番目の電流ダイポーは、各測定位置につくる磁
場の法線成分の２乗の和の平方根であり、よって１９ｏ
１はｎ番目の電流ダイポールの電流の大きさを各検出器
とダイポールの路数で補正したものに相当する。この補
正された電流の大きさの総和を示すのが（６）式のＥＳ
である。 （５）式あるいは（６）式によって表わされるＥｓ（ｑ
＋・・・？ＱＮ）は実際に含まれるダイポールの個数よ
りも多くのダイポールを仮定したことによって導入され
る過剰な自由度をおさえる働きをする。このことを理論
的に厳密に証明することは、本発明の対象とする脳内の
磁場発生源推定のような非線型逆問題の場合にはむずか
しいが、線型逆問題の場合には数学的に証明することが
できる。 この証明については羽石他ｒＣＴ画像再構成にシミュー
テッドアニーリングを用いる場合のコスト関数について
解析」 アプライドオプティックス印刷中 （Ｈ，Ｈａｎｅｉｓｈｉ　ｅｔ　ａｌ、　”Ａｎ　ａｎ
ａｌｙｓｉｓ　ｏｆ　ｃｏｓｔｆｕｎｃｔｉｏｎ　　ｕ
ｓｅｄ　　ｉｎ　　ｓｉｍｕｌａｔｅｄ　　ａｎｎｅｓ
ｌｉｎｇ　　ｆｏｒ　　ＣＴｉｍａｇｅ　ｒｅｃｏｎｓ
ｔｒｕｃｔｉｏｎ”、Ａｐｐｌｉｅｄ　０ｐｒｉｃｓ　
１ｎＰｒｅｓｓ、） （５）式あるいは（６）式に含まれるＷは（４）式右辺
第１項と第２項の重みづけを行なう定数であり、後はど
述べる最適化演算に際してコスト関数の変化ΔＥを計算
する際に（４）式の第１項あるいは第２項のどちらかか
らだけの寄与に片よらないように、あらかじめ試みにΔ
Ｅを計算してみることによりＷを設定する。 さて、本発明を実施する際に（４）式中のＥｓの項とし
ては（５）式あるいは（６）式を直接用いてもよいが、
これらを拡張した。 あるいは を用いるとさらに好ましい。ここでα法線成分の実数で
ある。 （５）′および（６）２式のαの役割は次の通りである
。 本発明においては脳内部に含まれると予想されるダイポ
ールの個数より多くのダイポールを仮定し推定を行なう
。したがって、複数個のダイポーは、ほぼ同一場所に推
定される場合が起りうる。 （５）２式および（６）２式のαはこのような場合にお
いて得られる解の性質を決めるものである。今、２個の
ダイポールｕ、ｖがほぼ同一場所に推定される場合、す
なわちｍｕ”ｒｖの場合を例として説明すると （１）Ｏ＜αくｌの場合、ｕ、ｖ番目のダイポールの電
推定Ｑｕｔ　Ｑｖのどちらかが非常に小さくななる。す
なわち、空間座標が同一の推定値となる複数個のダイポ
ールはそのうちの一個の電流推定値が本来その位置に存
在するダイポールの電流ベクトル値に近いものとなり、
他のダボールの電流ベクトル値はゼロに近いものとなる
。 （２）α〉１の場合、ｑｕ□ａ１となり、またこれらの
値はｍｕ”ｒｖの位置に本来存在するダイポールの電流
ベクトルの値の約１／２となる。 （３）ａ＝１の場合、ｑｕ＋ｑｖが本来ｒｕ”ｒｖの位
置に存在するダイポールの電流ベクトル値に近いものに
なるがＩｑ、ｌと１ｑ７１の比はどのような値となるか
は予測できない。 さて、この前述の（２）、（３）の場合においては同一
の位置推定結果を得た（複数値の）ダイポールについて
それらのダイポールの電流ベクトル値を加え合わせて、
その位置に存在するダイポールの電流ベクトル値とする
操作が必要となる。 一方、（１）の場合にはほぼ同一位置に推定された複数
の電流ベクトルのうちひとつを除いて残りはほぼゼロと
なるので実際に存在する個数より多い数の電流ダイポー
は、仮定されても自動的に実際に存在する電流ダイポー
は、解として得られる。すなわち０〈α〈１とするのが
最も好ましい。 次に（４）式により、くわしくはさらに（５）式（６）
式、（５）’式、（６）２式のいずれかにより定義され
たコスト関数Ｅ（ｒ工、・・・ｒＮ；（ｎ　＝　１　、
・・・Ｎ）を求める最適化演算について詳述する。コス
ト関数Ｅはその変数に関し非線型であり、局所的な最小
値を持つ。したがって、従来知られている非線型最適化
手法では真の最小値（大域的最小値）を求めることは不
可能である。 本実施例ではこのような場合でも大域的最小値を求める
ことのできる最近提案されたシミュレーティッド・アニ
ーリング（ｓｉｍｕｌａｔｅｄ　ａｎｎｅａｌｊｎｇ）
を採用している。 以下第２図を参照して説明する。 まず、上述した数Ｎの電流ダイポールを仮定し、それぞ
れの座ＩＮｒｎ　（ｎ＝Ｌ　”’ｔ　Ｎ）、および電流
ベクトルロイ（ｎ　＝　１　、・・・、Ｎ）の推定値を
初期値として設定する。なお、ダイポール座標ｒｎ　（
１：＝ｌ、２ｔ　”’ｔ　Ｎ）は変数ｘｎ（ｎ＝１゜２
、・・・、Ｎ）として、また各ダイポールの電流ベクト
ルｑｎ　（ｎ　”　１　＋　２　ｔ　”’＋　Ｎ）は変
数ｘｎ（ｎ＝Ｎ＋１．Ｎ＋２．・・・・・・、２Ｎ）と
して設定される。さらに実測され、記録されている磁場
計測値Ｄ　＊　（ｍ　〜１　ｓ・・・、Ｍ）を用いて、
前に詳述したコスト関数Ｅの値を算出する。（ステップ
１０１）以下、ｘｏの値の微少変化を与え（これを試行
という）、その変位により生じるコスト関数Ｅの変化Δ
Ｅを参照してこの試行を受け入るかどうか判断し、徐々
に各Ｘｎの値を最適値に近付けて行く。 まず、十分大きくパラメータＴを設定する。ＴはΔＥ＞
Ｏでも変化を受け入れる確率を決めるパラメータであり
、温度と呼ばれる。試し計算によりΔＥを計算し、　ｅ
ｘｐ　（−ΔＥ／Ｔ）〜０．８〜０．９となるようにＴ
を設定する。また、変数ｒ工。 ・・・　ｒＮに対して一回の試行における変位量１Δｑ
１を設定する。（ステップ１０２）は予想されるｌａｌ
の最大値の１７１０〜１／１０００程度が適切である。 次に、試行のくり返しの制御のためのパラメータＮＴ、
Ｎ工＋Ｎ２にそれぞれ初期値ゼロを設定する６　（ステ
ップ１０３） 次に試行の対象となる変数Ｘｎをひとつ選び（ステップ
１０６）、変化量ベクトルΔＸｎを決める（ステップ１
０７）。これは変数Ｘｎが座標であるなら を用い、 Ｘｎが電流ベクトルであるなら を用いて、θをＯくθく２７ｃの一様乱数で、ψをＯく
ψ（πの一様乱数を発生して決める。 次にステップ１０８では上記の変数Ｘｎをｘｎ＋ΔＸｎ
でおきかえる。つまり試行を実行する。 さらにステップ１０９では変位Δｘｒｌを与える前と後
のコスト関数Ｅの差ΔＥを ΔＥ　＝Ｅ（ｘｌ、−、ｘｌｌ＋ΔＸｎ＊　　”’ｙ　
　Ｘ２Ｎ）−Ｅ（ＸＬＩ−９Ｘｎｔ−ｔ　　Ｘ２Ｎ）　
　　（１０）から計算する。次にステップ１１１でΔＥ
〈０であれが変位ΔＸｌｌを受け入れる。すなわちｘｎ
＋ΔＸ０を新しいｘｌ、としたままステップ１１２゜１
０５を介して次の変数×１＋１の試行のためにステップ
１０６にもどる。また、ΔＥ＞Ｏである場合はステップ
１１３，１１４および１１６に示されるように温度Ｔに
依存するある確率Ｐ（ΔＥ）＝ｅｘρ（−ΔＥ／Ｔ）に
従がって変位Δｘｒｌを受け入れるか、拒絶するかを決
める。変位を受け入れると決められれば、Ｘｎ＋ΔＸｎ
を新しいｘｎとしたままステップ１１５，１０５を介し
てステップ１０６にもどる。一方変位を拒絶すると決め
られた場合はステップ１１７でｘ、、及びコスト関数Ｅ
を元の値にもどし、ステップ１１８．１０５を介してス
テップ１０６に戻る。 このような一連のステップを各変数について、つまりｎ
＝１〜２Ｎまでひととおり行ない、トータルでこれをＮ
　ＴＩＩ　ａ　２回くり返す。Ｎｔ’″ａｘとしては例
えば１００〜４００程度の値が適当であろう。 ここで、第２図のフローチャートのステップ１１２にお
けるＮ１はコスト関数を減らす方向の変位で受けいれら
れたものの数、ステップ１１５におけるＮ２はコススト
関数を増やす方向の変位で受け入れられたものの数、ス
テップ１１８におけるＮ３は拒絶された変位の数である
。ステップ１０５で、Ｎ丁″″０回の試行が終了したこ
とが判明したときステップ１１９に進み、変位を受け入
れる確率を定める温度Ｔを下げるか否かの判定を行なう
、すなわち、受けいれられた試行のうちコスト関数を増
やした回数Ｎ□がコスト関数を減らした回数Ｎ２に近い
値になるか否かをＩＮ□−Ｎ２１／Ｎの値により判定す
る。この値がε　（εはたとえば０．０２程度）以上で
あれば、温度Ｔを変えずにステップ１０３にもどり、さ
らにＮＴＩＩａｘ回の試行をくる返す。ＩＮよ−Ｎ２１
　／Ｎの値がεを下回った場合にはステップ１２０に進
み、温度Ｔを下げることによりΔＥ＞Ｏとするよう　変
位を受け入れる確率を下げてからステップ１０３にもど
って同様にＮ　Ｔ　”　ａ　１回の試行をくり返えす。 ステップ１２０における温度Ｔの下げ方としては本実施
例ではＴｋをに番目の温度段階として、Ｔ　ｈ　＝ξＴ
ｋ−１とする方法をとる。このときξとしては０．９〜
０．９５程度の値を用いる。別の下げ方Ｔｋ＝Ｔ、／　
（１＋ｋ）あるいはＴｈ＝Ｔｏ／ｌｏｇ（ｅ十ｋ）も文
献 Ｈ，スー「ファースト　シミュレーテッドアニーリング
」フィジックス　レタースＡ　Ｖｏｌ、１２２゜ｐ１５
７．　１９８７（Ｈ，Ｓｚｕ　　ｅｔ　　ａｌ、、　　
”Ｆａｓｔ　　ＳｉｍｕｌａｔｅｄＡｎｎｅａｌｉｎｇ
″、　Ｐｈｙｓｉｃｓ　Ｌｅｔｔｅｒｓ　Ａ、　Ｖｏｌ
、１２２．　Ｐ２Ｓ５゜１９８７） に提案されている。 以上のようにして、温度Ｔを徐々に下げながら各変数の
変位のＮ　Ｔｍ　＆　１回゛の試行を何巡かくり返すに
従がい、各変数は最適値に近づいて行き、受け入れられ
る変位は徐々に少なくなる。ステップ１１８にてＮ１＋
Ｎ、＝Ｏと判定されたとき、つまり試行された変位のう
ち受け入れられたものがなくなったときコスト関数を最
小とする変数、つまり最適な電流ダイポールの位置と電
流ベクトは、求められたので、シミュレーテッドアニー
リングの演算を終了する。 以上、本実施例では最適化演算にｓｉｍｕｌａｔｅｄａ
ｎｎｅａｌｉｎｇを用いる事を例とて説明したが、本発
明はこれに限るものではない。局所的最小値を含む関数
から真の最小値を見出すことのできる他のアルゴリズム
、例えば Ｊ、　Ｈ，Ｈｏｌｌａｎｄ氏提案のＧｅｎｅｔｉｃ　Ａ
ｌｇｏｒｉｔｈｍなども用いることができる。このアル
ゴリズムについては（Ｊ、　Ｈ，Ｉ（ｏｌｌａｎｄ　ｒ
Ａｄａｐｔａｔｉｏｎ　ｉｎ　Ｎａｔｕｒａｌａｎｄ　
Ａｒｆｉｆｉｃｉａｌ　ＳｙｓｔｅｍｓＪ　Ｔｈｅ　Ｕ
ｎｉｖｅｒｓｉｔｙ　ｏｆＭｉｃｈｉｇａｎ　Ｐｒｅｓ
ｓ　１９７８）で提案された議論されている。 またあらかじめダイポールの存在する領域が先験的に知
られている場合には以下のようにしてこの先験情報を最
適化演算の中に組み入れることができる。 すなわち、この場合には（４）式のがわりにで計算され
るＥを最小にする各推定値を求める。 ここでＥＬは次のように定義する。 ’ｉ？・・・＋　　ｒＮがすべてあらかじめ設定された
領域に存在する場合　ＥＬＣｒｌ、−、ｒＮ）　＝＝□
ｒ工、・・・　ｒＮのどれががあらかじめ設定された領
域を出た場合　Ｅｔ、（ｒよ、・・・　ｒＮ）＝Δここ
でＡはそのときの温度Ｔに対して十分大きな値、すなわ
ち　ｅｘｐ　（−Ａ／Ｔ）　＃Ｏであるような値である
。このようにＥＬを定義することにより、各ｒｋがあら
かじめ設定された領域から出るような変位を受け入れる
確率はゼロに近いものとなり、あらかじめ設定された領
域内で Ｅ　（ｒｌｙ　”’　　ｒＮ；　Ｑｔ＋　”’＋　ＱＮ
）を最小とするような推定値を求めることができる。 〔発明の効果〕 以上のごとく本発明によれば脳内に含まれる電流ダイポ
ールの数を正確に知らない場合でもダイポールの座標、
電流ベクトル値を正確に推定できる。 （

【図面の簡単な説明】

第１図は本発明の実施例の脳磁計測の概念図と。 リングのフローチャートを示す。

Claims (1)

1. 【特許請求の範囲】 １、生体の活動にともなって発生する磁場を複数の測定
   点で計測し、得られた複数の磁場計測値から前記生体の
   内部の電流分布を推定する推定方法において、前記生体
   内に存在すると予想されるより多くの数の電流ダイポー
   ルを仮定してその位置推定値と電流ベクトル推定値を設
   定し、前記位置推定値と電流ベクトル推定値から前記各
   測定点での仮想的な磁場計測値を算出し、前記の実際に
   計測された磁場計測値の分布と算出された仮想的な磁場
   計測値の分布との差を示す第１の量と前記電流ベクトル
   推定値の各々から算出される量を示す第２の量を少なく
   とも含むコスト関数が最小となる各推定値を求めて前記
   電流分布を推定することを特徴とする生体磁場計測にお
   ける磁場発生源の推定方法。 ２、前記第２の量は前記各電流ベクトルの絶対値の総和
   であること特徴とする請求項１に記載の推定方法。 ３、前記第３の量は前記各電流ベクトルの絶対値のα乗
   （αは正の実数）の総和であることを特徴とする。 ４、α≠１とすることを特徴とする請求項３の推定方法
   。 ５、０＜α＜１とすることを特徴とする請求項３の推定
   方法。 ６、前記第２の量は、前記の仮定された各電流ダイポー
   ルが前記各測定点に作る磁場の法線成分の２乗の計測点
   に関する和の平方根の各電流ダイポールに関する総和で
   あることを特徴とする請求項１の推定方法。 ７、前記第２の量は、前記の仮定された各電流ダイポー
   ルが前記各測定点に作る磁場の法線成分の２乗の測定点
   に関する和の平方根のα乗（αは正の実数）の各電流ダ
   イポールに関する総和であることを特徴とする請求項１
   の推定方法。 ８、α≠１である請求項７の推定方法。 ９、０＜α＜１である請求項７の推定方法。 １０、前記コスト関数は前記第１の量と前記第２の量の
   線形和である請求項１の推定方法。 １１、生体の活動にともなって発生する磁場を複数の測
   定点で計測し、得られた複数の磁場計測値から前記生体
   の内部の電流分布を推定する推定方法において、前記生
   体内に存在すると予想されるより多くの数の電流ダイポ
   ールを仮定してその位置指定値と電流ベクトル推定値を
   設定し、前記位置推定値と電流ベクトル推定値から前記
   各計測点での仮想的な磁場計測値を算出し、前記の実際
   に計測された磁場計測値の分布と算出された仮想的な分
   布との差を示す第１の量と前記電流ダイポールの電流値
   の総和に関する第２の量を少なくとも含むコスト関数の
   変化を参照しながら前記推定値の各々の変位の試行をく
   り返し、変位が前記コスト関数を減するもののときこの
   変位を受け入れ、一方変位が前コスト関数を増加させる
   もののとき温度Ｔに関するある確率でこの変位を受け入
   れ、上記温度Ｔを徐々に低下させながら上記変位の試行
   をくり返して前記コスト関数を最小とする各推定値に特
   定することを特徴とする生体磁場計測における磁場発生
   源の推定方法。