JP3139414B2 - 生体内活動部位推定方法及び生体内活動部位推定装置 - Google Patents
生体内活動部位推定方法及び生体内活動部位推定装置Info
- Publication number
- JP3139414B2 JP3139414B2 JP09140362A JP14036297A JP3139414B2 JP 3139414 B2 JP3139414 B2 JP 3139414B2 JP 09140362 A JP09140362 A JP 09140362A JP 14036297 A JP14036297 A JP 14036297A JP 3139414 B2 JP3139414 B2 JP 3139414B2
- Authority
- JP
- Japan
- Prior art keywords
- estimating
- electromagnetic field
- grid point
- field distribution
- loss function
- 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.)
- Expired - Fee Related
Links
Landscapes
- Measurement And Recording Of Electrical Phenomena And Electrical Characteristics Of The Living Body (AREA)
Description
【０００１】
【発明の属する技術分野】本発明は、生体体表面頭皮上
で測定された電磁場分布に基づいて、生体内の電流ダイ
ポールを推定する方法及び装置に関する。
で測定された電磁場分布に基づいて、生体内の電流ダイ
ポールを推定する方法及び装置に関する。
【０００２】
【従来の技術】従来より、生体内電流ダイポール推定に
おいては、特に、人間の頭皮上で測定された電磁場分布
から脳内の活動部位を推定する手法として用いられてお
り、この活動部位を知ることによって脳の高次機能や脳
内疾患部位に関する知見を得ることができる。
おいては、特に、人間の頭皮上で測定された電磁場分布
から脳内の活動部位を推定する手法として用いられてお
り、この活動部位を知ることによって脳の高次機能や脳
内疾患部位に関する知見を得ることができる。
【０００３】その一般的な方法は以下のようなものであ
る。
る。
【０００４】脳内の神経細胞（ニューロン）が外部から
の刺激によって興奮すると、ニューロン同士を結ぶ結線
（軸索）にパルス状の電流が流れるが、その電流によっ
て頭部の周囲に微弱な電磁場が生じる。その電磁場の源
泉を電流ダイポール（以下、ダイポールと称する）と呼
ばれる電流素片でモデル化し、ダイポールの６個のパラ
メータを推定することによって脳内活動部位を推定す
る。ここで、６個のパラメータとは、位置を指定する３
成分（ｘ’）、方向を指定する２成分（ｅ）、およびダ
イポールの大きさ（ｑ）である。ダイポールの方向と大
きさとをまとめてｑｅと書き、モーメント・ベクトルと
呼ぶ場合もある。以下では、これらの６個のパラメータ
をまとめてθで表すものとする。
の刺激によって興奮すると、ニューロン同士を結ぶ結線
（軸索）にパルス状の電流が流れるが、その電流によっ
て頭部の周囲に微弱な電磁場が生じる。その電磁場の源
泉を電流ダイポール（以下、ダイポールと称する）と呼
ばれる電流素片でモデル化し、ダイポールの６個のパラ
メータを推定することによって脳内活動部位を推定す
る。ここで、６個のパラメータとは、位置を指定する３
成分（ｘ’）、方向を指定する２成分（ｅ）、およびダ
イポールの大きさ（ｑ）である。ダイポールの方向と大
きさとをまとめてｑｅと書き、モーメント・ベクトルと
呼ぶ場合もある。以下では、これらの６個のパラメータ
をまとめてθで表すものとする。
【０００５】ダイポールのパラメータを推定する従来の
方法は、以下のように行われるのが一般的である。以下
では、頭皮上で測定した磁場分布に基づいて脳内の活動
部位を推定する方法を例に挙げて説明するが、脳以外の
生体部位においても同様な方法で活動部位の推定が行わ
れる。また、磁場分布の代わりに電位分布を用いる場合
でも本質的な違いはない。
方法は、以下のように行われるのが一般的である。以下
では、頭皮上で測定した磁場分布に基づいて脳内の活動
部位を推定する方法を例に挙げて説明するが、脳以外の
生体部位においても同様な方法で活動部位の推定が行わ
れる。また、磁場分布の代わりに電位分布を用いる場合
でも本質的な違いはない。
【０００６】従来の手法においては、以下に定義される
二乗誤差Ｅ
二乗誤差Ｅ
【０００７】
【数１】 を最小にするパラメータを推定値とする。ここに、
【０００８】
【外１】 は磁場の測定値、ｘiはｉ番目の測定点の座標、ｎは測
定点の数、ｍはダイポールの数である。ｆ（ｘi，θj）
は、パラメータθjで指定されるｊ番目のダイポールに
よってｉ番目の測定点に生じる磁場の理論値であり、以
下のように計算することができる。
定点の数、ｍはダイポールの数である。ｆ（ｘi，θj）
は、パラメータθjで指定されるｊ番目のダイポールに
よってｉ番目の測定点に生じる磁場の理論値であり、以
下のように計算することができる。
【０００９】
【数２】
【００１０】ここで、ｎiはｉ番目の測定点における磁
場センサの法線成分であり、‖・‖はベクトルのノルム
を表す。これらのダイポール・パラメータを推定する方
法としては、例えば、Levenberg-Marquardt法やPowell
法などの非線形最適化アルゴリズムが用いられている。
これら非線形最適化アルゴリズムにおいては、例えば今
野浩、山下浩著、「非線形計画法」、日科技連、１９８
７年（以下、文献１と称する）に詳しく記載されてい
る。
場センサの法線成分であり、‖・‖はベクトルのノルム
を表す。これらのダイポール・パラメータを推定する方
法としては、例えば、Levenberg-Marquardt法やPowell
法などの非線形最適化アルゴリズムが用いられている。
これら非線形最適化アルゴリズムにおいては、例えば今
野浩、山下浩著、「非線形計画法」、日科技連、１９８
７年（以下、文献１と称する）に詳しく記載されてい
る。
【００１１】磁場を計算するその他の方法として、マト
リクスを用いた方法がある。この場合、脳の形状を近似
した多面体を用い、多面体の各格子点に大きさが１のダ
イポールを向きを固定して配置する。ダイポールの向き
は通常多面体の面に垂直に固定する。このとき、ダイポ
ールによって生成される磁場分布ｆは、
リクスを用いた方法がある。この場合、脳の形状を近似
した多面体を用い、多面体の各格子点に大きさが１のダ
イポールを向きを固定して配置する。ダイポールの向き
は通常多面体の面に垂直に固定する。このとき、ダイポ
ールによって生成される磁場分布ｆは、
【００１２】
【数３】 となる。ここでｆはｎ個の測定点における磁場の理論値
を並べたｎ次元実ベクトル、ｑは、ｍ個のダイポールの
大きさを表す変数を並べたｍ次元実ベクトルである。ｎ
行ｍ列のマトリクスＦの成分は、固定されたダイポール
の位置および方向から、以下のように表される。
を並べたｎ次元実ベクトル、ｑは、ｍ個のダイポールの
大きさを表す変数を並べたｍ次元実ベクトルである。ｎ
行ｍ列のマトリクスＦの成分は、固定されたダイポール
の位置および方向から、以下のように表される。
【００１３】
【数４】
【００１４】この定式化によって、推定すべきパラメー
タは、各格子点にあるダイポールの大きさｑのみとな
る。この場合、二乗誤差‖Ｂ0−Ｆq‖2を最小にするダ
イポール分布
タは、各格子点にあるダイポールの大きさｑのみとな
る。この場合、二乗誤差‖Ｂ0−Ｆq‖2を最小にするダ
イポール分布
【００１５】
【外２】 は、線形方程式（正規方程式）
【００１６】
【数５】 を解くことによって得られる。しかしながら、生体内活
動部位推定では通常Ｆが特異であるため、一意な解が得
られない場合が多い。そこで、二乗誤差を最小にする解
の中で、その大きさ−ｑ−の最小な解が最適であるとい
う制約条件をつけて一意な解を得る方法が知られてい
る。この方法においては、Ｆの擬逆行列（または一般逆
行列と呼ばれる）Ｆ+を用いて、最小二乗解が
動部位推定では通常Ｆが特異であるため、一意な解が得
られない場合が多い。そこで、二乗誤差を最小にする解
の中で、その大きさ−ｑ−の最小な解が最適であるとい
う制約条件をつけて一意な解を得る方法が知られてい
る。この方法においては、Ｆの擬逆行列（または一般逆
行列と呼ばれる）Ｆ+を用いて、最小二乗解が
【００１７】
【数６】 で与えられる。
【００１８】擬逆行列を用いたパラメータ推定として
は、例えばG.H.Golub and C.F.Van Loan, "Matrix Comp
utation (Third Edition)," Johns Hopkins University
Press, 1996（以下、文献２と称する）に詳しく記載さ
れている。
は、例えばG.H.Golub and C.F.Van Loan, "Matrix Comp
utation (Third Edition)," Johns Hopkins University
Press, 1996（以下、文献２と称する）に詳しく記載さ
れている。
【００１９】上述したような従来の手法においては、二
乗誤差Ｅを最小にするパラメータが最適なパラメータで
あるとされていた。しかしながら、統計学の教えるとこ
ろによると、測定したデータ（ここでは磁場分布）に対
して設定したモデル（ここではダイポール）で用いられ
るパラメータの数が多くなれば多くなるほど、二乗誤差
が小さく出ることが知られている。これは、余分なパラ
メータのために、真の測定データとは無関係な雑音にま
でモデルが過剰適合してしまうからである。したがって
この場合、二乗誤差を最小にするパラメータが得られた
としても本来のデータが持つ構造をとらえているとは限
らない。すなわち、ダイポール数が未知の場合に推定の
評価関数として二乗誤差を用いると、正しいダイポール
数を推定することができない。
乗誤差Ｅを最小にするパラメータが最適なパラメータで
あるとされていた。しかしながら、統計学の教えるとこ
ろによると、測定したデータ（ここでは磁場分布）に対
して設定したモデル（ここではダイポール）で用いられ
るパラメータの数が多くなれば多くなるほど、二乗誤差
が小さく出ることが知られている。これは、余分なパラ
メータのために、真の測定データとは無関係な雑音にま
でモデルが過剰適合してしまうからである。したがって
この場合、二乗誤差を最小にするパラメータが得られた
としても本来のデータが持つ構造をとらえているとは限
らない。すなわち、ダイポール数が未知の場合に推定の
評価関数として二乗誤差を用いると、正しいダイポール
数を推定することができない。
【００２０】脳内活動部位推定の例で言えば、余分なダ
イポールを用いることによって、本来活動していない部
位まで活動部位として推定してしまう可能性がある。こ
れを医療に応用した場合は、正常な部位を切除してしま
うような重大な過失を招く虞れがある。
イポールを用いることによって、本来活動していない部
位まで活動部位として推定してしまう可能性がある。こ
れを医療に応用した場合は、正常な部位を切除してしま
うような重大な過失を招く虞れがある。
【００２１】推定の評価基準として二乗誤差以外の量を
用いる方法としては、例えば、特開平３−２７７３４５
号公報（以下、文献３と称する）に開示されているもの
がある。
用いる方法としては、例えば、特開平３−２７７３４５
号公報（以下、文献３と称する）に開示されているもの
がある。
【００２２】図５は、特開平３−２７７３４５号公報に
開示されているダイポール推定方法を示すフローチャー
トである。
開示されているダイポール推定方法を示すフローチャー
トである。
【００２３】なお、ここではその概略だけを述べるにと
どめる。
どめる。
【００２４】この方法においては、磁場の測定値と理論
値との二乗誤差と、ダイポール・モーメントの総和をコ
スト関数Ｅとして計算する（ステップＳ１０１）。
値との二乗誤差と、ダイポール・モーメントの総和をコ
スト関数Ｅとして計算する（ステップＳ１０１）。
【００２５】コスト関数Ｅを式で表せば、
【００２６】
【数７】 となる。ここでωは定数、αは正の実数である。乱数を
発生させてダイポール・パラメータ（位置またはモーメ
ント・ベクトル）をΔＸだけ変化させる（ステップＳ１
０８）。
発生させてダイポール・パラメータ（位置またはモーメ
ント・ベクトル）をΔＸだけ変化させる（ステップＳ１
０８）。
【００２７】ΔＸの大きさとしては、あらかじめ予想さ
れる変数Ｘの最大値を基準として、その１０分の１から
１０００分の１とする。その後、このパラメータの変位
によるコスト関数の変化量ΔＥを計算し（ステップＳ１
１０）、ΔＥの値が負ならばその変位を受け入れ（ステ
ップＳ１１１）、ΔＥが正の場合はΔＥから計算される
確率分布に従ってその変位を受け入れる（ステップＳ１
１６）。
れる変数Ｘの最大値を基準として、その１０分の１から
１０００分の１とする。その後、このパラメータの変位
によるコスト関数の変化量ΔＥを計算し（ステップＳ１
１０）、ΔＥの値が負ならばその変位を受け入れ（ステ
ップＳ１１１）、ΔＥが正の場合はΔＥから計算される
確率分布に従ってその変位を受け入れる（ステップＳ１
１６）。
【００２８】この方法の特徴としては、二乗誤差の代わ
りに、二乗誤差とダイポール・モーメントの大きさとの
総和を用いていることと、ダイポール・パラメータを離
散化し、シミュレーテッド・アニーリングと呼ばれる方
法で最適なパラメータを探索することである。シミュレ
ーテッド・アニーリングに関しては、S.Kirkpatricketa
l., "Optimization by Simulated Annealing," Scienc
e, Vol.220, pp.671-680, 1983（以下、文献と称する
４）に詳しく記載されている。
りに、二乗誤差とダイポール・モーメントの大きさとの
総和を用いていることと、ダイポール・パラメータを離
散化し、シミュレーテッド・アニーリングと呼ばれる方
法で最適なパラメータを探索することである。シミュレ
ーテッド・アニーリングに関しては、S.Kirkpatricketa
l., "Optimization by Simulated Annealing," Scienc
e, Vol.220, pp.671-680, 1983（以下、文献と称する
４）に詳しく記載されている。
【００２９】
【発明が解決しようとする課題】しかしながら、上述し
たような従来のものにおいては以下に記載するような問
題点がある。
たような従来のものにおいては以下に記載するような問
題点がある。
【００３０】例えば、文献３においては、ダイポール・
モーメントの大きさを評価関数に取り込むことによって
自由度が制限されるため、ダイポールの数が未知である
場合においても最適なダイポール数が得られるとされて
いるが、これには理論的根拠が無い。このことは、式
（７）における二つの定数ω、αの値の決め方が明確に
与えられていないことからもわかる。
モーメントの大きさを評価関数に取り込むことによって
自由度が制限されるため、ダイポールの数が未知である
場合においても最適なダイポール数が得られるとされて
いるが、これには理論的根拠が無い。このことは、式
（７）における二つの定数ω、αの値の決め方が明確に
与えられていないことからもわかる。
【００３１】また、ダイポール・モーメントの総和がな
るべく小さくなるような解を得ようとする制約条件は、
擬逆行列を用いた方法と同じであるが、この方法では得
られた活動部位が広がってしまい、はっきりとした活動
部位が得られないことが知られている。
るべく小さくなるような解を得ようとする制約条件は、
擬逆行列を用いた方法と同じであるが、この方法では得
られた活動部位が広がってしまい、はっきりとした活動
部位が得られないことが知られている。
【００３２】さらに、パラメータの変更方法の効率が悪
いため、計算時間がかかるという問題点もある。例え
ば、あるダイポールの大きさを０にするためには、最悪
の場合１０００回の繰り返しが必要となる。
いため、計算時間がかかるという問題点もある。例え
ば、あるダイポールの大きさを０にするためには、最悪
の場合１０００回の繰り返しが必要となる。
【００３３】本発明は、上述したような従来の技術が有
する問題点に鑑みてなされたものであって、ダイポール
数が未知の場合でも生体内活動部位の位置と個数とを正
確にかつ高速に推定できる生体内活動部位推定方法及び
生体内活動推定装置を提供することを目的とする。
する問題点に鑑みてなされたものであって、ダイポール
数が未知の場合でも生体内活動部位の位置と個数とを正
確にかつ高速に推定できる生体内活動部位推定方法及び
生体内活動推定装置を提供することを目的とする。
【００３４】
【課題を解決するための手段】上記目的を達成するため
に本発明は、生体体表面上で測定された電磁場分布及び
形状データを入力とし、前記電磁場分布の源泉として１
個または複数の電流ダイポールを生体内部に仮定し、前
記電流ダイポールの方向、大きさ及び個数を推定する生
体内活動部位推定方法において、前記電流ダイポールに
よって生じる電磁場分布の理論値と前記電磁場分布との
誤差Ｅを最小にするように各電流ダイポールの方向及び
大きさを推定する第１の工程と、前記電流ダイポールの
位置する格子点に前記電流ダイポールの削除を示す有効
格子点変数を割り当て、前記電流ダイポールの数と、前
記電磁場分布の測定点の数及び前記誤差Ｅとから損失関
数Ｌを計算する第２の工程と、前記損失関数Ｌが減少す
るように前記有効格子点変数を更新する第３の工程とを
有することを特徴とする。
に本発明は、生体体表面上で測定された電磁場分布及び
形状データを入力とし、前記電磁場分布の源泉として１
個または複数の電流ダイポールを生体内部に仮定し、前
記電流ダイポールの方向、大きさ及び個数を推定する生
体内活動部位推定方法において、前記電流ダイポールに
よって生じる電磁場分布の理論値と前記電磁場分布との
誤差Ｅを最小にするように各電流ダイポールの方向及び
大きさを推定する第１の工程と、前記電流ダイポールの
位置する格子点に前記電流ダイポールの削除を示す有効
格子点変数を割り当て、前記電流ダイポールの数と、前
記電磁場分布の測定点の数及び前記誤差Ｅとから損失関
数Ｌを計算する第２の工程と、前記損失関数Ｌが減少す
るように前記有効格子点変数を更新する第３の工程とを
有することを特徴とする。
【００３５】また、前記第１の工程にて、各格子点に配
置された電流ダイポールが各電磁場分布の測定点に生成
する電磁場の大きさを指定するマトリクスＦを用いて電
磁場分布の理論値を算出し、前記第１の工程、前記第２
の工程及び前記第３の工程を実行し、前記第３の工程で
更新された有効格子点変数に基づいて、前記マトリクス
Ｆから列を削除する第４の工程を実行し、前記第１の工
程、前記第２の工程、前記第３の工程及び前記第４の工
程を所定の終了条件が満たされるまで繰り返し実行する
ことを特徴とする。
置された電流ダイポールが各電磁場分布の測定点に生成
する電磁場の大きさを指定するマトリクスＦを用いて電
磁場分布の理論値を算出し、前記第１の工程、前記第２
の工程及び前記第３の工程を実行し、前記第３の工程で
更新された有効格子点変数に基づいて、前記マトリクス
Ｆから列を削除する第４の工程を実行し、前記第１の工
程、前記第２の工程、前記第３の工程及び前記第４の工
程を所定の終了条件が満たされるまで繰り返し実行する
ことを特徴とする。
【００３６】また、前記第３の工程にて、前記有効格子
点変数の更新時に、更新の対象となる格子点をランダム
に選択し、所定の条件を満足した場合に前記選択された
格子点に位置する電流ダイポールを削除することを特徴
とする。
点変数の更新時に、更新の対象となる格子点をランダム
に選択し、所定の条件を満足した場合に前記選択された
格子点に位置する電流ダイポールを削除することを特徴
とする。
【００３７】また、前記第３の工程にて、各格子点に位
置する電流ダイポールの大きさが小さいものほど選択さ
れる確率が高くなる確率分布を前記格子点に割り当て、
前記確率分布に基づいて更新の対象となる格子点を選択
することを特徴とする。
置する電流ダイポールの大きさが小さいものほど選択さ
れる確率が高くなる確率分布を前記格子点に割り当て、
前記確率分布に基づいて更新の対象となる格子点を選択
することを特徴とする。
【００３８】また、前記第３の工程にて、各格子点に位
置する１個の電流ダイポールが、単独で生成する電磁場
分布と前記電磁場分布との単独誤差Ｅ1を計算し、前記
単独誤差Ｅ1が大きいほど選択される確率が高くなる確
率分布を前記格子点に割り当て、該確率分布に基づいて
更新の対象となる格子点を選択することを特徴とする。
置する１個の電流ダイポールが、単独で生成する電磁場
分布と前記電磁場分布との単独誤差Ｅ1を計算し、前記
単独誤差Ｅ1が大きいほど選択される確率が高くなる確
率分布を前記格子点に割り当て、該確率分布に基づいて
更新の対象となる格子点を選択することを特徴とする。
【００３９】また、前記第３の工程にて、前記有効格子
点変数の更新時に、更新の対象となる格子点をランダム
に選択し、前記選択された格子点を更新することによる
前記損失関数Ｌの変化量ΔＬを計算し、前記損失関数の
変化量ΔＬから更新受理確率を算出し、前記更新受理確
率に基づいて、前記格子点の更新を受理することを特徴
とする。
点変数の更新時に、更新の対象となる格子点をランダム
に選択し、前記選択された格子点を更新することによる
前記損失関数Ｌの変化量ΔＬを計算し、前記損失関数の
変化量ΔＬから更新受理確率を算出し、前記更新受理確
率に基づいて、前記格子点の更新を受理することを特徴
とする。
【００４０】また、前記第３の工程にて、各格子点に位
置する電流ダイポールの大きさが小さいものほど選択さ
れる確率が高くなる確率分布を前記格子点に割り当て、
前記確率分布に基づいて更新の対象となる格子点を選択
し、前記選択された格子点を更新することによる前記損
失関数Ｌの変化量ΔＬを計算し、前記損失関数の変化量
ΔＬから更新受理確率を算出し、前記更新受理確率に基
づいて、前記格子点の更新を受理することを特徴とす
る。
置する電流ダイポールの大きさが小さいものほど選択さ
れる確率が高くなる確率分布を前記格子点に割り当て、
前記確率分布に基づいて更新の対象となる格子点を選択
し、前記選択された格子点を更新することによる前記損
失関数Ｌの変化量ΔＬを計算し、前記損失関数の変化量
ΔＬから更新受理確率を算出し、前記更新受理確率に基
づいて、前記格子点の更新を受理することを特徴とす
る。
【００４１】また、前記第３の工程にて、各格子点に位
置する１個の電流ダイポールが、単独で生成する電磁場
分布と前記電磁場分布との単独誤差Ｅ1を計算し、前記
単独誤差Ｅ1が大きいほど選択される確率が高くなる確
率分布を前記格子点に割り当て、該確率分布に基づいて
更新の対象となる格子点を選択し、前記選択された格子
点を更新することによる前記損失関数Ｌの変化量ΔＬを
計算し、前記損失関数の変化量ΔＬから更新受理確率を
算出し、前記更新受理確率に基づいて、前記格子点の更
新を受理することを特徴とする。
置する１個の電流ダイポールが、単独で生成する電磁場
分布と前記電磁場分布との単独誤差Ｅ1を計算し、前記
単独誤差Ｅ1が大きいほど選択される確率が高くなる確
率分布を前記格子点に割り当て、該確率分布に基づいて
更新の対象となる格子点を選択し、前記選択された格子
点を更新することによる前記損失関数Ｌの変化量ΔＬを
計算し、前記損失関数の変化量ΔＬから更新受理確率を
算出し、前記更新受理確率に基づいて、前記格子点の更
新を受理することを特徴とする。
【００４２】また、前記損失関数Ｌとして、前記有効格
子点の数と前記電磁場分布の測定点の数および前記誤差
Ｅから計算される記述長を用いることを特徴とする。
子点の数と前記電磁場分布の測定点の数および前記誤差
Ｅから計算される記述長を用いることを特徴とする。
【００４３】また、前記損失関数Ｌとして、前記有効格
子点の数と前記電磁場分布の測定点の数及び前記誤差Ｅ
から計算される情報量基準を用いることを特徴とする。
子点の数と前記電磁場分布の測定点の数及び前記誤差Ｅ
から計算される情報量基準を用いることを特徴とする。
【００４４】また、前記損失関数Ｌとして、前記有効格
子点の数と前記電磁場分布の測定点の数及び前記誤差Ｅ
とから計算されるストラクチュラル・リスク（Structur
ul Risk）を用いることを特徴とする。
子点の数と前記電磁場分布の測定点の数及び前記誤差Ｅ
とから計算されるストラクチュラル・リスク（Structur
ul Risk）を用いることを特徴とする。
【００４５】また、前記損失関数Ｌとして、磁場分布の
測定点の数及び前記誤差Ｅとの和を用いることを特徴と
する。
測定点の数及び前記誤差Ｅとの和を用いることを特徴と
する。
【００４６】また、被検体の周囲に設けられ、前記被検
体の電磁場分布及び形状データを検出する検出手段と、
該検出手段にて検出されたデータを測定する測定手段
と、該測定手段において測定されたデータに基づいて演
算を行う演算手段とを有し、生体体表面上で測定された
電磁場分布及び形状データを入力とし、前記電磁場分布
の源泉として１個または複数の電流ダイポールを生体内
部に仮定し、前記電流ダイポールの方向、大きさ及び個
数を推定する生体内活動部位推定装置において、前記演
算手段は、前記電流ダイポールによって生じる電磁場分
布の理論値と前記電磁場分布との誤差Ｅを最小にするよ
うに各電流ダイポールの方向及び大きさを推定し、前記
電流ダイポールの位置する格子点に前記電流ダイポール
の削除を示す有効格子点変数を割り当て、前記電流ダイ
ポールの数と、前記電磁場分布の測定点の数及び前記誤
差Ｅとから損失関数Ｌを計算し、前記損失関数Ｌが減少
するように前記有効格子点変数を更新することを特徴と
する。
体の電磁場分布及び形状データを検出する検出手段と、
該検出手段にて検出されたデータを測定する測定手段
と、該測定手段において測定されたデータに基づいて演
算を行う演算手段とを有し、生体体表面上で測定された
電磁場分布及び形状データを入力とし、前記電磁場分布
の源泉として１個または複数の電流ダイポールを生体内
部に仮定し、前記電流ダイポールの方向、大きさ及び個
数を推定する生体内活動部位推定装置において、前記演
算手段は、前記電流ダイポールによって生じる電磁場分
布の理論値と前記電磁場分布との誤差Ｅを最小にするよ
うに各電流ダイポールの方向及び大きさを推定し、前記
電流ダイポールの位置する格子点に前記電流ダイポール
の削除を示す有効格子点変数を割り当て、前記電流ダイ
ポールの数と、前記電磁場分布の測定点の数及び前記誤
差Ｅとから損失関数Ｌを計算し、前記損失関数Ｌが減少
するように前記有効格子点変数を更新することを特徴と
する。
【００４７】（作用）ダイポールの個数が未知の場合
に、推定の良さの指標として、ダイポールの個数ｍを考
慮した次の損失関数ＬMDLを用いる。
に、推定の良さの指標として、ダイポールの個数ｍを考
慮した次の損失関数ＬMDLを用いる。
【００４８】
【数８】 ここで、Ｅは二乗誤差、ｎは測定点の数、ｍはダイポー
ルの数、ｌｎは自然対数である。この量は、情報理論に
おいて知られている記述長最小（Minimum Description
Length, MDL ）原理における損失関数であり、この量Ｌ
MDLを最小にするダイポール・パラメータおよびダイポ
ール数を最適であるとする。
ルの数、ｌｎは自然対数である。この量は、情報理論に
おいて知られている記述長最小（Minimum Description
Length, MDL ）原理における損失関数であり、この量Ｌ
MDLを最小にするダイポール・パラメータおよびダイポ
ール数を最適であるとする。
【００４９】上述したＬMDLの定義式（８）の第１項
は、二乗誤差に比例する項で、データとモデルとのフィ
ッティングの良さを表す。多数のダイポールを用いると
小さな二乗誤差を得ることができるが、ＬMDLの第２項
目はパラメータ数（ダイポールの数）に比例して大きく
なるため、ＬMDLの最小化という意味では不利になる。
すなわち、第２項目は過剰なパラメータ数に対するペナ
ルティ項の役割を担っている。この機構によって、過剰
なダイポールの使用を抑制することができ、適切なダイ
ポール数を推定することが可能となる。
は、二乗誤差に比例する項で、データとモデルとのフィ
ッティングの良さを表す。多数のダイポールを用いると
小さな二乗誤差を得ることができるが、ＬMDLの第２項
目はパラメータ数（ダイポールの数）に比例して大きく
なるため、ＬMDLの最小化という意味では不利になる。
すなわち、第２項目は過剰なパラメータ数に対するペナ
ルティ項の役割を担っている。この機構によって、過剰
なダイポールの使用を抑制することができ、適切なダイ
ポール数を推定することが可能となる。
【００５０】この記述長を用いたパラメータ数の推定に
関する理論的基盤は、例えば、Rissanen, Modeling by
Shortest Data Description, Automatica, Vol.14, pp.
465-471, 1978 （以下、文献５と称する）に詳しく記載
されている。
関する理論的基盤は、例えば、Rissanen, Modeling by
Shortest Data Description, Automatica, Vol.14, pp.
465-471, 1978 （以下、文献５と称する）に詳しく記載
されている。
【００５１】しかしながら、ダイポール数が未知の場
合、ＬMDLを最小化する過程においては、ダイポールの
個数を変化させながら全ての組み合わせを試みて、その
中からＬMDLを最小にするダイポールの組を見いださな
くてはならない。例えば、活動部位の候補となりうるダ
イポールの数が１００個ある場合、その組み合わせの数
Ｎは
合、ＬMDLを最小化する過程においては、ダイポールの
個数を変化させながら全ての組み合わせを試みて、その
中からＬMDLを最小にするダイポールの組を見いださな
くてはならない。例えば、活動部位の候補となりうるダ
イポールの数が１００個ある場合、その組み合わせの数
Ｎは
【００５２】
【数９】 となり、莫大な数の組み合わせを考慮しなければならな
い（
い（
【００５３】
【外３】 はｎ個の中からｒ個を選ぶ組み合わせの数）。
【００５４】本発明では、ＬMDLの最小化を効率的に行
うために、以下のような方法を用いる。まず、活動部位
のモデルとなるダイポールを多数（ダイポール数をｓと
する）用意し、磁場の測定値と理論値との二乗誤差を最
小にするダイポール・パラメータ
うために、以下のような方法を用いる。まず、活動部位
のモデルとなるダイポールを多数（ダイポール数をｓと
する）用意し、磁場の測定値と理論値との二乗誤差を最
小にするダイポール・パラメータ
【００５５】
【外４】 を擬逆行列を用いて求める。この場合、ダイポールの位
置は脳の形状を近似した格子点にあり、方向も固定され
ているので、未知パラメータは各ダイポールの大きさｑ
だけである。
置は脳の形状を近似した格子点にあり、方向も固定され
ているので、未知パラメータは各ダイポールの大きさｑ
だけである。
【００５６】次に、各格子点にあるダイポールの有効性
を指定する有効格子点変数σ＝（σ 1，・・・，σs）T
を導入する（Ｔは転置を表す）。各σi（ｉ＝１，・・
・，ｓ）は、ｉ番目の格子点に存在するダイポールを使
用するかしないかを表し、使用する場合は１、使用しな
い場合は０となる。この有効格子点変数を考慮すると、
二乗誤差Ｅは、
を指定する有効格子点変数σ＝（σ 1，・・・，σs）T
を導入する（Ｔは転置を表す）。各σi（ｉ＝１，・・
・，ｓ）は、ｉ番目の格子点に存在するダイポールを使
用するかしないかを表し、使用する場合は１、使用しな
い場合は０となる。この有効格子点変数を考慮すると、
二乗誤差Ｅは、
【００５７】
【数１０】 となる。また、損失関数ＬMDLは以下のように表すこと
ができる。
ができる。
【００５８】
【数１１】 ここで、
【００５９】
【数１２】 である。式（１３）を最小化するには、εMDL＝Ｅexp
(ＣMDLｍ）を最小化すればよい。εMDLを有効格子点変
数σを用いて表すと、
(ＣMDLｍ）を最小化すればよい。εMDLを有効格子点変
数σを用いて表すと、
【００６０】
【数１３】 ここで、
【００６１】
【数１４】 とする。
【００６２】この損失関数εMDLが減少するように有効
格子変数を変更するには、以下のようにする。
格子変数を変更するには、以下のようにする。
【００６３】まず、変更の対象となる更新点を選択す
る。これを例えばｋ番目の格子点とする。
る。これを例えばｋ番目の格子点とする。
【００６４】次に、ｋ番目の格子点にあるダイポールを
有効とするかどうかを示す有効格子点変数σkをσ'kに
変更する。この変更による損失関数の変化ΔεMDLは、
有効とするかどうかを示す有効格子点変数σkをσ'kに
変更する。この変更による損失関数の変化ΔεMDLは、
【００６５】
【数１５】 となる。ここで、
【００６６】
【数１６】 である。Ｅk＞０かつｙk＞０のときにσ'k＝０となるよ
うに有効格子点変数を更新すれば、
うに有効格子点変数を更新すれば、
【００６７】
【数１７】 となり、変更前の値σkにかかわらず、εMDLは減少する
か、もしくは変化しない。
か、もしくは変化しない。
【００６８】このようにして本発明では、不要なダイポ
ールを削除し、損失関数εMDLすなわちＬMDLを最小にす
るダイポール分布を効率よく推定することができる。
ールを削除し、損失関数εMDLすなわちＬMDLを最小にす
るダイポール分布を効率よく推定することができる。
【００６９】ここでは損失関数としてＬMDLを用いて説
明したが、後述する実施の形態から明らかになるよう
に、本発明では、式（８）に示す損失関数の他に、例え
ば、赤池の情報量基準ＬAICに基づく評価関数も使用す
ることが可能であり、また、ストラクチュラル・リスク
最小化原理に基づく損失関数ＬSRMも使用することが可
能である。
明したが、後述する実施の形態から明らかになるよう
に、本発明では、式（８）に示す損失関数の他に、例え
ば、赤池の情報量基準ＬAICに基づく評価関数も使用す
ることが可能であり、また、ストラクチュラル・リスク
最小化原理に基づく損失関数ＬSRMも使用することが可
能である。
【００７０】
【発明の実施の形態】以下に、本発明の実施の形態につ
いて図面を参照して説明する。なお、以下に用いる記号
は特に言及しない限り今まで用いてきた記号の説明に従
うものとする。特に、損失関数Ｌ、あるいはεについて
は、陽にＬMDL、ＬAICなどと記していない場合には、ど
の損失関数を用いてもよいことを意味する。また、以下
の実施の形態では、生体（以下、被験者と呼ぶ）の頭皮
上で測定された磁場分布から、脳内の活動部位を推定す
る方法について説明するが、例えば、被験者の胸部で電
磁場分布を測定し、心臓内の活動部位を推定する場合で
も同様な方法で行うことができる。さらには、消化器官
や筋肉中での活動部位の推定にも用いることができる。
いて図面を参照して説明する。なお、以下に用いる記号
は特に言及しない限り今まで用いてきた記号の説明に従
うものとする。特に、損失関数Ｌ、あるいはεについて
は、陽にＬMDL、ＬAICなどと記していない場合には、ど
の損失関数を用いてもよいことを意味する。また、以下
の実施の形態では、生体（以下、被験者と呼ぶ）の頭皮
上で測定された磁場分布から、脳内の活動部位を推定す
る方法について説明するが、例えば、被験者の胸部で電
磁場分布を測定し、心臓内の活動部位を推定する場合で
も同様な方法で行うことができる。さらには、消化器官
や筋肉中での活動部位の推定にも用いることができる。
【００７１】図１は、本発明の生体内活動部位推定方法
の実施の一形態を示すフローチャートである。
の実施の一形態を示すフローチャートである。
【００７２】まず、ステップＳ１において、被験者の頭
皮上のｎ箇所で測定された電磁場分布データおよび被験
者の頭部形状、脳の形状などに関するデータを読み込
む。
皮上のｎ箇所で測定された電磁場分布データおよび被験
者の頭部形状、脳の形状などに関するデータを読み込
む。
【００７３】電磁場分布データの測定においては、電場
分布（すなわち脳波）を測定する場合は脳波計が、磁場
分布を測定する場合は超伝導量子干渉素子デバイス（Su
perconducting Quantum Interference Device, SQUID）
がそれぞれ用いられる。通常、測定点は被験者の頭皮上
に２０箇所から１５０箇所に設定される。測定の際に
は、測定点の座標や頭部の基準点（例えば鼻根、左右の
耳など）の位置も同時に記録しておく。
分布（すなわち脳波）を測定する場合は脳波計が、磁場
分布を測定する場合は超伝導量子干渉素子デバイス（Su
perconducting Quantum Interference Device, SQUID）
がそれぞれ用いられる。通常、測定点は被験者の頭皮上
に２０箇所から１５０箇所に設定される。測定の際に
は、測定点の座標や頭部の基準点（例えば鼻根、左右の
耳など）の位置も同時に記録しておく。
【００７４】また、頭部形状は、例えば、ＭＲＩ装置や
Ｘ線ＣＴスキャンなどを用いて測定する。頭部形状の測
定には、上述したような装置が必要であることや時間が
かかるという短所があるため、実際の頭部形状の代わり
として、例えば、頭部を球で近似したモデルを用いても
よい。球モデルを用いた場合には、推定精度が多少劣る
が、計算負荷を軽減する効果がある。
Ｘ線ＣＴスキャンなどを用いて測定する。頭部形状の測
定には、上述したような装置が必要であることや時間が
かかるという短所があるため、実際の頭部形状の代わり
として、例えば、頭部を球で近似したモデルを用いても
よい。球モデルを用いた場合には、推定精度が多少劣る
が、計算負荷を軽減する効果がある。
【００７５】次に、ステップＳ２において、ステップＳ
１にて読み込んだ頭部形状データに基づいて、被験者の
頭部を多面体、例えば三角形の集合で近似し、格子点を
生成する。
１にて読み込んだ頭部形状データに基づいて、被験者の
頭部を多面体、例えば三角形の集合で近似し、格子点を
生成する。
【００７６】ここで、頭部モデルは、ＭＲＩを用いて撮
影された頭部断層画像から構成されたものが考えられ
る。
影された頭部断層画像から構成されたものが考えられ
る。
【００７７】さらに、脳内活動部位のモデルとして、各
三角形の面内に大きさ１の単位ダイポールを配置する。
このときのモデルダイポールの数をｍとする。ダイポー
ルの方向は未知数として求めることもできるが、例え
ば、その方向を三角形の面に対して垂直に固定して推定
することも可能である。ダイポールの方向をこのように
設定するのは、皮質にある神経細胞の軸索が、皮質の表
面に対してほぼ垂直になっているという生理学的知見に
基づいている。以下、本形態においては、ダイポールの
方向を上述のように三角形の面に対して垂直に固定し、
大きさのみを求める場合について説明するが、方向を未
知として推定する場合も全く同様な方法で推定すること
が可能である。
三角形の面内に大きさ１の単位ダイポールを配置する。
このときのモデルダイポールの数をｍとする。ダイポー
ルの方向は未知数として求めることもできるが、例え
ば、その方向を三角形の面に対して垂直に固定して推定
することも可能である。ダイポールの方向をこのように
設定するのは、皮質にある神経細胞の軸索が、皮質の表
面に対してほぼ垂直になっているという生理学的知見に
基づいている。以下、本形態においては、ダイポールの
方向を上述のように三角形の面に対して垂直に固定し、
大きさのみを求める場合について説明するが、方向を未
知として推定する場合も全く同様な方法で推定すること
が可能である。
【００７８】次に、ステップＳ３において、格子点上に
配置されたダイポールを推定に使用するか否かを設定す
る有効格子点変数を初期化する。
配置されたダイポールを推定に使用するか否かを設定す
る有効格子点変数を初期化する。
【００７９】推定の開始時点では全てのダイポールを用
いるので、全ての有効格子点変数を１にする。σt＝
（１，…，１）T。
いるので、全ての有効格子点変数を１にする。σt＝
（１，…，１）T。
【００８０】それ以降は、変数の上付添字ｔによってｔ
回目のループにおける変数の値を表すことにする。ま
た、損失関数の最小値を記録する変数Ｌminをなるべく
大きな数に設定する。Ｌminの初期値としては、例え
ば、倍精度変数を用いる場合Ｌmin＝２64−１などとす
る。
回目のループにおける変数の値を表すことにする。ま
た、損失関数の最小値を記録する変数Ｌminをなるべく
大きな数に設定する。Ｌminの初期値としては、例え
ば、倍精度変数を用いる場合Ｌmin＝２64−１などとす
る。
【００８１】次に、ステップＳ４において、カウンタｔ
を０に初期化する。
を０に初期化する。
【００８２】次に、ステップＳ５において、擬逆行列を
用いて、二乗誤差Ｅ＝‖Ｂ0−Ｆｑ‖2を最小にするダイ
ポール分布を推定する。
用いて、二乗誤差Ｅ＝‖Ｂ0−Ｆｑ‖2を最小にするダイ
ポール分布を推定する。
【００８３】ダイポール分布の推定は以下のように行
う。
う。
【００８４】推定の開始時点では、全てのダイポールを
用いるので、擬逆行列Ｆ+を測定データＢ0に掛けるだけ
でよい。擬逆行列は、マトリクスＦの特異値分解（Sing
ular Value Decomposition, SVD）
用いるので、擬逆行列Ｆ+を測定データＢ0に掛けるだけ
でよい。擬逆行列は、マトリクスＦの特異値分解（Sing
ular Value Decomposition, SVD）
【００８５】
【数１８】 を利用する。ここで、Ｕ，Λ，Ｖはそれぞれｎ×ｍ，ｍ
×ｍ，ｍ×ｍのマトリクスであり、ｕi，ｖiは、それぞ
れマトリクスＵ，Ｖのｉ列を表す列ベクトルである。ま
た、diag（λ1，λ2，…）は、Ｆの特異値λ1，λ2，…
を対角成分とする対角行列を表し、ｒはＦのランクを表
す。ＳＶＤの結果を用いて、Ｆの一般逆行列Ｆ+は次の
ように書ける。
×ｍ，ｍ×ｍのマトリクスであり、ｕi，ｖiは、それぞ
れマトリクスＵ，Ｖのｉ列を表す列ベクトルである。ま
た、diag（λ1，λ2，…）は、Ｆの特異値λ1，λ2，…
を対角成分とする対角行列を表し、ｒはＦのランクを表
す。ＳＶＤの結果を用いて、Ｆの一般逆行列Ｆ+は次の
ように書ける。
【００８６】
【数１９】 さらに、この一般逆行列を用いれば二乗誤差Ｅを最小に
するパラメータ
するパラメータ
【００８７】
【外５】 が
【００８８】
【数２０】 で与えられる。式（２９）における和は、必ずしもｒま
で取る必要はなく、例えば、最大特異値λmaxとｉ番目
の特異値との比が、ある正数ｅよりも大きい間だけ和を
取る方法や、二乗誤差があらかじめ推定される雑音分布
の大きさよりも大きい間だけ取る方法なども可能であ
る。これによって雑音の影響を小さくする効果がある。
で取る必要はなく、例えば、最大特異値λmaxとｉ番目
の特異値との比が、ある正数ｅよりも大きい間だけ和を
取る方法や、二乗誤差があらかじめ推定される雑音分布
の大きさよりも大きい間だけ取る方法なども可能であ
る。これによって雑音の影響を小さくする効果がある。
【００８９】次に、ステップＳ６において、ステップＳ
５で推定したダイポールが生成する磁場と磁場の測定値
との二乗誤差Ｅ、測定点の数ｎおよびダイポールの数
５で推定したダイポールが生成する磁場と磁場の測定値
との二乗誤差Ｅ、測定点の数ｎおよびダイポールの数
【００９０】
【数２１】 を用いて、損失関数を計算する。
【００９１】この損失関数としては、例えば、次式で定
義される記述長
義される記述長
【００９２】
【数２２】 または赤池の情報量基準（Akaike Information Criteri
on, AIC ）
on, AIC ）
【００９３】
【数２３】 などを用いることができる。この赤池の情報量基準を用
いたパラメータ数の推定に関する理論的基盤は、例えば
文献、石黒ほか、情報量統計学、共立出版、１９８３
（以下、文献６と称する）に詳しく記載されている。
いたパラメータ数の推定に関する理論的基盤は、例えば
文献、石黒ほか、情報量統計学、共立出版、１９８３
（以下、文献６と称する）に詳しく記載されている。
【００９４】また、次式で定義されるストラクチュラル
・リスク最小化（Structural RiskMinimization, SRM
）原理の損失関数を用いても良い。
・リスク最小化（Structural RiskMinimization, SRM
）原理の損失関数を用いても良い。
【００９５】
【数２４】 ここで、ｄVC(m)はパラメータ数ｍ個のモデルの“表現
力”の指標となる量で、Vapnik-Chervonenkis（ＶＣ）
次元と呼ばれるものである。ηはＥ0を観測データに依
らない真の二乗誤差として、不等式
力”の指標となる量で、Vapnik-Chervonenkis（ＶＣ）
次元と呼ばれるものである。ηはＥ0を観測データに依
らない真の二乗誤差として、不等式
【００９６】
【数２５】 が確率１−ηで成立することを表し、［ｘ］+＝Ｍａｘ
（ｘ，０）である。定数ｃ＝ａ（ｐ）τは、以下の量か
ら計算される。
（ｘ，０）である。定数ｃ＝ａ（ｐ）τは、以下の量か
ら計算される。
【００９７】
【数２６】 ここで、ｄμ（Ｂ0）はデータＢ0の発生する確率分布で
ある。分布ｄμ（Ｂ0）は一般に未知であるため、τ
は、
ある。分布ｄμ（Ｂ0）は一般に未知であるため、τ
は、
【００９８】
【数２７】 で計算した値を用いてもよい。評価関数ＬSRMが個数推
定に有効であることの理論的基盤は、例えば、V.N.Vapn
ik, The Nature of Statistical Learning Theory, Spr
inger, NY, 1995 （以下、文献７と称する）に詳しく記
載されている。
定に有効であることの理論的基盤は、例えば、V.N.Vapn
ik, The Nature of Statistical Learning Theory, Spr
inger, NY, 1995 （以下、文献７と称する）に詳しく記
載されている。
【００９９】損失関数としてＬSRMを用いた場合、作用
で説明した方法をそのまま適用することができないた
め、以下のようにする。
で説明した方法をそのまま適用することができないた
め、以下のようにする。
【０１００】例えば、損失関数を二乗誤差とダイポール
の数の和
の数の和
【０１０１】
【数２８】 とし、αは常に不等式ＬSRM≦Ｌαを満たすように定め
る。このようにすれば、Ｌαを最小にするパラメータを
探索することによって、ＬSRMを最小化することが可能
となる。
る。このようにすれば、Ｌαを最小にするパラメータを
探索することによって、ＬSRMを最小化することが可能
となる。
【０１０２】次に、損失関数が減少するように有効格子
点変数を更新する（ステップＳ７）。
点変数を更新する（ステップＳ７）。
【０１０３】更新の方法は、作用の項で説明した方法を
用いる。また、ここで用いられる記号は、作用の項の説
明に準ずる。
用いる。また、ここで用いられる記号は、作用の項の説
明に準ずる。
【０１０４】図２は、図１に示した生体内活動部位推定
方法における有効格子点変数の更新方法の実施の一形態
を示すフローチャートである。
方法における有効格子点変数の更新方法の実施の一形態
を示すフローチャートである。
【０１０５】ステップＳ１３においては、ｓ個の格子点
の中から更新の対象となる格子点（例えばｋ番目の格子
点）を一様乱数を用いて選択する。格子点の選択方法は
この限りでなく、例えば、格子点に存在するダイポール
のモーメントが小さいほど高い確率で選択されるような
確率分布を割り当てることも可能である。これによっ
て、あまり重要でないダイポールを優先的に削除するこ
とができ、効率的である。また、各格子点に存在するダ
イポールが単独で生成する磁場と測定磁場との二乗誤差
（これを単独二乗誤差とする）を計算し、この単独二乗
誤差が大きいものほど高い確率で選択されるような確率
分布を割り当てることも可能である。これによって、二
乗誤差の減少に貢献しないダイポールを優先的に削除す
ることができる。
の中から更新の対象となる格子点（例えばｋ番目の格子
点）を一様乱数を用いて選択する。格子点の選択方法は
この限りでなく、例えば、格子点に存在するダイポール
のモーメントが小さいほど高い確率で選択されるような
確率分布を割り当てることも可能である。これによっ
て、あまり重要でないダイポールを優先的に削除するこ
とができ、効率的である。また、各格子点に存在するダ
イポールが単独で生成する磁場と測定磁場との二乗誤差
（これを単独二乗誤差とする）を計算し、この単独二乗
誤差が大きいものほど高い確率で選択されるような確率
分布を割り当てることも可能である。これによって、二
乗誤差の減少に貢献しないダイポールを優先的に削除す
ることができる。
【０１０６】次に、ステップＳ１４において、
【０１０７】
【数２９】 を計算し、ステップＳ１５へ進む。
【０１０８】ステップＳ１５においては、ｙk＞０なら
ばステップＳ１６に進み、ｙk≦０ならばステップＳ８
へ進む。
ばステップＳ１６に進み、ｙk≦０ならばステップＳ８
へ進む。
【０１０９】ステップＳ１６においては、格子点ｋにあ
るダイポールの二乗誤差Ｅへの寄与Ｅk＝Ｅ−ｙkを計算
し、ステップＳ１７へ進む。
るダイポールの二乗誤差Ｅへの寄与Ｅk＝Ｅ−ｙkを計算
し、ステップＳ１７へ進む。
【０１１０】ステップＳ１７においては、Ｅk＞０なら
ばステップＳ１８ヘ進み、Ｅk≦０ならばステップＳ８
へ進む。
ばステップＳ１８ヘ進み、Ｅk≦０ならばステップＳ８
へ進む。
【０１１１】ステップＳ１８においては、
【０１１２】
【外６】 ＝０とし、ｋ番目の格子点にあるダイポールを削除す
る。
る。
【０１１３】上述したステップＳ１３からステップＳ１
８までの工程は、以下に述べるように確率的に行うこと
も可能である。
８までの工程は、以下に述べるように確率的に行うこと
も可能である。
【０１１４】図３は、図１に示した生体内活動部位推定
方法における有効格子点変数の更新方法の実施の他の形
態を示すフローチャートであり、確率的な有効格子点変
数の変更方法を示している。
方法における有効格子点変数の更新方法の実施の他の形
態を示すフローチャートであり、確率的な有効格子点変
数の変更方法を示している。
【０１１５】ステップＳ１９においては、ステップＳ１
３における方法と同様な方法で変更対象となる格子点を
選択し、ステップＳ２０に進む。
３における方法と同様な方法で変更対象となる格子点を
選択し、ステップＳ２０に進む。
【０１１６】ステップＳ２０においては、選択されたｋ
番目の有効格子点変数の値を反転させる。
番目の有効格子点変数の値を反転させる。
【０１１７】なお、
【０１１８】
【外７】
【０１１９】
【外８】 値の反転とは、現在の値が１ならば０に、１ならば０に
する事を言う。
する事を言う。
【０１２０】ステップＳ２１においては、ステップＳ２
０における変更による損失関数の変化Δεを計算し、ス
テップＳ２２に進む。
０における変更による損失関数の変化Δεを計算し、ス
テップＳ２２に進む。
【０１２１】ステップＳ２２においては、変更受理確率
Ｐ＝ｅｘｐ［−Δε／Ｔt］／Ｚを計算する。ここで、
Ｚは正規化定数、Ｔtはループのカウントがｔの時の温
度パラメータである。この温度パラメータの値が大きけ
れば、損失関数の変化（増加または減少）の方向がラン
ダムになり、小さければ減少する方向への変化が生じや
すくなる。温度パラメータＴtは、ループのカウンタの
増加とともに減少するように設定する。その減少のさせ
方としては、例えば、ｃ／ｌｏｇ（１＋ｔ）≧Ｔtを満
たすように減少させる方法や、ｃ／（１＋ｔ）＝Ｔtと
する方法などがある。ここでｃは正数である。
Ｐ＝ｅｘｐ［−Δε／Ｔt］／Ｚを計算する。ここで、
Ｚは正規化定数、Ｔtはループのカウントがｔの時の温
度パラメータである。この温度パラメータの値が大きけ
れば、損失関数の変化（増加または減少）の方向がラン
ダムになり、小さければ減少する方向への変化が生じや
すくなる。温度パラメータＴtは、ループのカウンタの
増加とともに減少するように設定する。その減少のさせ
方としては、例えば、ｃ／ｌｏｇ（１＋ｔ）≧Ｔtを満
たすように減少させる方法や、ｃ／（１＋ｔ）＝Ｔtと
する方法などがある。ここでｃは正数である。
【０１２２】このように、損失関数が増加する方向への
変化も確率的に許容することにより、局所最小解に陥る
ことを防止し、確実に最小解を探索することが可能とな
る。この温度パラメータの減少規則に関しては、例え
ば、S.Geman and D.Geman, Stochastic Relaxation, Gi
bbs Distributions, and the Bayesian Restoration of
Images, IEEE Trans. Patt. Ana. Mach. Int., PAMI-6,
P.721, 1984 （以下、文献８と称する）などに詳しく
記載されている。
変化も確率的に許容することにより、局所最小解に陥る
ことを防止し、確実に最小解を探索することが可能とな
る。この温度パラメータの減少規則に関しては、例え
ば、S.Geman and D.Geman, Stochastic Relaxation, Gi
bbs Distributions, and the Bayesian Restoration of
Images, IEEE Trans. Patt. Ana. Mach. Int., PAMI-6,
P.721, 1984 （以下、文献８と称する）などに詳しく
記載されている。
【０１２３】ステップＳ２３においては、一様乱数ｒを
発生させ、ステップＳ２４に進む。
発生させ、ステップＳ２４に進む。
【０１２４】ステップＳ２４においては、ｒとＰとを比
較し、ｒ＞ＰならばステップＳ２５に進み、そうでなけ
ればステップＳ２６に進む。
較し、ｒ＞ＰならばステップＳ２５に進み、そうでなけ
ればステップＳ２６に進む。
【０１２５】ステップＳ２５においては、反転された格
子点変数を受理し、
子点変数を受理し、
【０１２６】
【外９】 に変更を保存してステップＳ８へ進む。
【０１２７】ステップＳ２６においては、反転された格
子点変数は受理せず、
子点変数は受理せず、
【０１２８】
【外１０】 を変更せずにステップＳ８へ進む。
【０１２９】ここで、図１に戻り、ステップＳ８以降の
処理について説明する。
処理について説明する。
【０１３０】ステップＳ８においては、現時点での損失
関数の値ＬtとＬminとを比較し、Ｌ t≦Ｌminならばステ
ップＳ９へ進み、Ｌt＞ＬminならばステップＳ１１に進
む。
関数の値ＬtとＬminとを比較し、Ｌ t≦Ｌminならばステ
ップＳ９へ進み、Ｌt＞ＬminならばステップＳ１１に進
む。
【０１３１】ステップＳ９では、現時点での損失関数の
値ＬtをＬminに、有効格子点σtの値をσminにそれぞれ
保存し、ステップＳ１０に進む。
値ＬtをＬminに、有効格子点σtの値をσminにそれぞれ
保存し、ステップＳ１０に進む。
【０１３２】ステップＳ１０においては、ループのカウ
ンタｔの値を一つ増やし、ステップＳ７へ戻る。
ンタｔの値を一つ増やし、ステップＳ７へ戻る。
【０１３３】ステップＳ１１においては、所定の終了条
件が満たされているかどうかを調べ、終了条件が満たさ
れているならば処理を終了し、満たされていなければス
テップＳ１２へ進む。
件が満たされているかどうかを調べ、終了条件が満たさ
れているならば処理を終了し、満たされていなければス
テップＳ１２へ進む。
【０１３４】終了条件としては、ループのカウンタｔが
所定のある一定値ｔmaxに達すること、あるいは損失関
数がある一定期間減少しないこと、あるいは温度パラメ
ータＴが所定の最小値Ｔminに到達すること、などとす
る。
所定のある一定値ｔmaxに達すること、あるいは損失関
数がある一定期間減少しないこと、あるいは温度パラメ
ータＴが所定の最小値Ｔminに到達すること、などとす
る。
【０１３５】ステップＳ１２においては、有効格子点変
数σtが０となっている格子点に対応する列をマトリク
スＦから削除し、ステップＳ５へ進み、更新されたマト
リクスＦを用いて、再度ダイポール分布を推定する。ダ
イポールを削除することによって損失関数Ｌは減少する
が、削除後のダイポール分布は、二乗誤差をなるべく小
さくするという意味ではまだ改善の余地がある。そこ
で、削除されていないダイポールのパラメータを再度推
定し直すことにより、さらに二乗誤差を小さくし、その
結果として損失関数を小さくすることができる。
数σtが０となっている格子点に対応する列をマトリク
スＦから削除し、ステップＳ５へ進み、更新されたマト
リクスＦを用いて、再度ダイポール分布を推定する。ダ
イポールを削除することによって損失関数Ｌは減少する
が、削除後のダイポール分布は、二乗誤差をなるべく小
さくするという意味ではまだ改善の余地がある。そこ
で、削除されていないダイポールのパラメータを再度推
定し直すことにより、さらに二乗誤差を小さくし、その
結果として損失関数を小さくすることができる。
【０１３６】以下に、上述した生体内活動部位推定方法
が適用された生体内部位推定装置の一構成例について簡
単に説明する。
が適用された生体内部位推定装置の一構成例について簡
単に説明する。
【０１３７】図４は、図１に示した生体内活動部位推定
方法が適用された生体内部位推定装置の実施の一形態を
示す図である。
方法が適用された生体内部位推定装置の実施の一形態を
示す図である。
【０１３８】図４に示すように、被検体の頭部４の周囲
に設けられ、被検体頭部４の電磁場分布及び形状データ
を検出する検出手段であるピックアップコイル３と、ピ
ックアップコイル３において検出されたデータを測定す
る測定手段２と、測定手段２において測定されたデータ
に基づいて演算を行う演算手段１とから構成されてお
り、演算手段１において、上述した一連の演算動作が行
われる。
に設けられ、被検体頭部４の電磁場分布及び形状データ
を検出する検出手段であるピックアップコイル３と、ピ
ックアップコイル３において検出されたデータを測定す
る測定手段２と、測定手段２において測定されたデータ
に基づいて演算を行う演算手段１とから構成されてお
り、演算手段１において、上述した一連の演算動作が行
われる。
【０１３９】
【発明の効果】本発明は、以上説明したように構成され
ているので、二乗誤差を評価関数として用いる従来方法
では不可能であったダイポールの個数推定が可能とな
る。
ているので、二乗誤差を評価関数として用いる従来方法
では不可能であったダイポールの個数推定が可能とな
る。
【０１４０】また、本発明においては、まず擬逆行列を
用いて推定を行い、その後に不要なダイポールを変更す
るという方法を用いているため、繰り返し個数の変更と
評価関数の比較をする必要がある従来手法よりも高速な
個数推定が可能となる。
用いて推定を行い、その後に不要なダイポールを変更す
るという方法を用いているため、繰り返し個数の変更と
評価関数の比較をする必要がある従来手法よりも高速な
個数推定が可能となる。
【図１】本発明の生体内活動部位推定方法の実施の一形
態を示すフローチャートである。
態を示すフローチャートである。
【図２】図１に示した生体内活動部位推定方法における
有効格子点変数の更新方法の実施の一形態を示すフロー
チャートである。
有効格子点変数の更新方法の実施の一形態を示すフロー
チャートである。
【図３】図１に示した生体内活動部位推定方法における
有効格子点変数の更新方法の実施の他の形態を示すフロ
ーチャートである。
有効格子点変数の更新方法の実施の他の形態を示すフロ
ーチャートである。
【図４】図１に示した生体内活動部位推定方法が適用さ
れた生体内部位推定装置の実施の一形態を示す図であ
る。
れた生体内部位推定装置の実施の一形態を示す図であ
る。
【図５】特開平３−２７７３４５号公報に開示されてい
るダイポール推定方法を示すフローチャートである。
るダイポール推定方法を示すフローチャートである。
１ 演算手段 ２ 測定手段 ３ ピックアップコイル ４ 被検体頭部
───────────────────────────────────────────────────── フロントページの続き (56)参考文献 特開 平９−66038（ＪＰ，Ａ) 特開 平９−56688（ＪＰ，Ａ) 特開 平４−8347（ＪＰ，Ａ) 特開 平３−277345（ＪＰ，Ａ) 特開 平６−121776（ＪＰ，Ａ) 特開 平８−299295（ＪＰ，Ａ) 特開 平６−311975（ＪＰ，Ａ) (58)調査した分野(Int.Cl.7，ＤＢ名) A61B 5/05
Claims (13)
- 【請求項１】 生体体表面上で測定された電磁場分布及
び形状データを入力とし、前記電磁場分布の源泉として
１個または複数の電流ダイポールを生体内部に仮定し、
前記電流ダイポールの方向、大きさ及び個数を推定する
生体内活動部位推定方法において、 前記電流ダイポールによって生じる電磁場分布の理論値
と前記電磁場分布との誤差Ｅを最小にするように各電流
ダイポールの方向及び大きさを推定する第１の工程と、 前記電流ダイポールの位置する格子点に前記電流ダイポ
ールの削除を示す有効格子点変数を割り当て、前記電流
ダイポールの数と、前記電磁場分布の測定点の数及び前
記誤差Ｅとから損失関数Ｌを計算する第２の工程と、 前記損失関数Ｌが減少するように前記有効格子点変数を
更新する第３の工程とを有することを特徴とする生体内
活動部位推定方法。 - 【請求項２】 請求項１に記載の生体内部位測定方法に
おいて、 前記第１の工程にて、各格子点に配置された電流ダイポ
ールが各電磁場分布の測定点に生成する電磁場の大きさ
を指定するマトリクスＦを用いて電磁場分布の理論値を
算出し、 前記第１の工程、前記第２の工程及び前記第３の工程を
実行し、 前記第３の工程で更新された有効格子点変数に基づい
て、前記マトリクスＦから列を削除する第４の工程を実
行し、 前記第１の工程、前記第２の工程、前記第３の工程及び
前記第４の工程を所定の終了条件が満たされるまで繰り
返し実行することを特徴とする生体内活動部位推定方
法。 - 【請求項３】 請求項１または請求項２に記載の生体内
活動部位推定方法において、 前記第３の工程にて、前記有効格子点変数の更新時に、
更新の対象となる格子点をランダムに選択し、所定の条
件を満足した場合に前記選択された格子点に位置する電
流ダイポールを削除することを特徴とする生体内活動部
位推定方法。 - 【請求項４】 請求項１または請求項２に記載の生体内
活動部位推定方法において、 前記第３の工程にて、各格子点に位置する電流ダイポー
ルの大きさが小さいものほど選択される確率が高くなる
確率分布を前記格子点に割り当て、前記確率分布に基づ
いて更新の対象となる格子点を選択することを特徴とす
る生体内活動部位推定方法。 - 【請求項５】 請求項１または請求項２に記載の生体内
活動部位推定方法において、 前記第３の工程にて、各格子点に位置する１個の電流ダ
イポールが、単独で生成する電磁場分布と前記電磁場分
布との単独誤差Ｅ1を計算し、前記単独誤差Ｅ1が大きい
ほど選択される確率が高くなる確率分布を前記格子点に
割り当て、該確率分布に基づいて更新の対象となる格子
点を選択することを特徴とする生体内活動部位推定方
法。 - 【請求項６】 請求項１または請求項２に記載の生体内
活動部位推定方法において、 前記第３の工程にて、前記有効格子点変数の更新時に、
更新の対象となる格子点をランダムに選択し、 前記選択された格子点を更新することによる前記損失関
数Ｌの変化量ΔＬを計算し、 前記損失関数の変化量ΔＬから更新受理確率を算出し、 前記更新受理確率に基づいて、前記格子点の更新を受理
することを特徴とする生体内活動部位推定方法。 - 【請求項７】 請求項１または請求項２に記載の生体内
活動部位推定方法において、 前記第３の工程にて、各格子点に位置する電流ダイポー
ルの大きさが小さいものほど選択される確率が高くなる
確率分布を前記格子点に割り当て、前記確率分布に基づ
いて更新の対象となる格子点を選択し、 前記選択された格子点を更新することによる前記損失関
数Ｌの変化量ΔＬを計算し、 前記損失関数の変化量ΔＬから更新受理確率を算出し、 前記更新受理確率に基づいて、前記格子点の更新を受理
することを特徴とする生体内活動部位推定方法。 - 【請求項８】 請求項１または請求項２に記載の生体内
活動部位推定方法において、 前記第３の工程にて、各格子点に位置する１個の電流ダ
イポールが、単独で生成する電磁場分布と前記電磁場分
布との単独誤差Ｅ1を計算し、前記単独誤差Ｅ1が大きい
ほど選択される確率が高くなる確率分布を前記格子点に
割り当て、該確率分布に基づいて更新の対象となる格子
点を選択し、 前記選択された格子点を更新することによる前記損失関
数Ｌの変化量ΔＬを計算し、 前記損失関数の変化量ΔＬから更新受理確率を算出し、 前記更新受理確率に基づいて、前記格子点の更新を受理
することを特徴とする生体内活動部位推定方法。 - 【請求項９】 請求項１乃至３のいずれか１項に記載の
生体内活動部位推定方法において、 前記損失関数Ｌとして、前記有効格子点の数と前記電磁
場分布の測定点の数および前記誤差Ｅから計算される記
述長を用いることを特徴とする生体内活動部位推定方
法。 - 【請求項１０】 請求項１乃至３のいずれか１項に記載
の生体内活動部位推定方法において、 前記損失関数Ｌとして、前記有効格子点の数と前記電磁
場分布の測定点の数及び前記誤差Ｅから計算される情報
量基準を用いることを特徴とする生体内活動部位推定方
法。 - 【請求項１１】 請求項１乃至３のいずれか１項に記載
の生体内活動部位推定方法において、 前記損失関数Ｌとして、前記有効格子点の数と前記電磁
場分布の測定点の数及び前記誤差Ｅとから計算されるス
トラクチュラル・リスク（Structurul Risk）を用いる
ことを特徴とする生体内活動部位推定方法。 - 【請求項１２】 請求項１乃至３のいずれか１項に記載
の生体内活動部位推定方法において、 前記損失関数Ｌとして、磁場分布の測定点の数及び前記
誤差Ｅとの和を用いることを特徴とする生体内活動部位
推定方法。 - 【請求項１３】 被検体の周囲に設けられ、前記被検体
の電磁場分布及び形状データを検出する検出手段と、 該検出手段にて検出されたデータを測定する測定手段
と、 該測定手段において測定されたデータに基づいて演算を
行う演算手段とを有し、 生体体表面上で測定された電磁場分布及び形状データを
入力とし、前記電磁場分布の源泉として１個または複数
の電流ダイポールを生体内部に仮定し、前記電流ダイポ
ールの方向、大きさ及び個数を推定する生体内活動部位
推定装置において、 前記演算手段は、 前記電流ダイポールによって生じる電磁場分布の理論値
と前記電磁場分布との誤差Ｅを最小にするように各電流
ダイポールの方向及び大きさを推定し、 前記電流ダイポールの位置する格子点に前記電流ダイポ
ールの削除を示す有効格子点変数を割り当て、前記電流
ダイポールの数と、前記電磁場分布の測定点の数及び前
記誤差Ｅとから損失関数Ｌを計算し、 前記損失関数Ｌが減少するように前記有効格子点変数を
更新することを特徴とする生体内活動部位推定装置。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| JP09140362A JP3139414B2 (ja) | 1997-05-29 | 1997-05-29 | 生体内活動部位推定方法及び生体内活動部位推定装置 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| JP09140362A JP3139414B2 (ja) | 1997-05-29 | 1997-05-29 | 生体内活動部位推定方法及び生体内活動部位推定装置 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| JPH10328155A JPH10328155A (ja) | 1998-12-15 |
| JP3139414B2 true JP3139414B2 (ja) | 2001-02-26 |
Family
ID=15267069
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| JP09140362A Expired - Fee Related JP3139414B2 (ja) | 1997-05-29 | 1997-05-29 | 生体内活動部位推定方法及び生体内活動部位推定装置 |
Country Status (1)
| Country | Link |
|---|---|
| JP (1) | JP3139414B2 (ja) |
Cited By (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN108498089A (zh) * | 2018-05-08 | 2018-09-07 | 北京邮电大学 | 一种基于深度神经网络的无创连续血压测量方法 |
Families Citing this family (4)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US7403265B2 (en) * | 2005-03-30 | 2008-07-22 | Asml Netherlands B.V. | Lithographic apparatus and device manufacturing method utilizing data filtering |
| US8406848B2 (en) * | 2009-10-06 | 2013-03-26 | Seiko Epson Corporation | Reconstructing three-dimensional current sources from magnetic sensor data |
| JP6566513B2 (ja) * | 2015-04-21 | 2019-08-28 | 学校法人 関西大学 | 心容積及び心拍出量の推定装置 |
| JP7181454B2 (ja) | 2018-11-12 | 2022-12-01 | 富士通株式会社 | 最適化装置、最適化装置の制御方法および最適化装置の制御プログラム |
-
1997
- 1997-05-29 JP JP09140362A patent/JP3139414B2/ja not_active Expired - Fee Related
Cited By (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN108498089A (zh) * | 2018-05-08 | 2018-09-07 | 北京邮电大学 | 一种基于深度神经网络的无创连续血压测量方法 |
| CN108498089B (zh) * | 2018-05-08 | 2022-03-25 | 北京邮电大学 | 一种基于深度神经网络的无创连续血压测量方法 |
Also Published As
| Publication number | Publication date |
|---|---|
| JPH10328155A (ja) | 1998-12-15 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| US8406848B2 (en) | Reconstructing three-dimensional current sources from magnetic sensor data | |
| US20230301542A1 (en) | Brain atlas individualization method and system based on magnetic resonance and twin graph neural network | |
| EP2499585B1 (en) | Methods and systems for channel selection | |
| JP3033508B2 (ja) | 生体内活動部位推定方法 | |
| Jun et al. | Spatiotemporal Bayesian inference dipole analysis for MEG neuroimaging data | |
| CN114052668B (zh) | 一种基于脑磁图数据的脑功能分析方法 | |
| Wismüller et al. | The deformable feature map-a novel neurocomputing algorithm for adaptive plasticity in pattern analysis | |
| Anninou et al. | Modeling of Parkinson's disease using fuzzy cognitive maps and non-linear Hebbian learning | |
| Park et al. | Bayesian active learning of neural firing rate maps with transformed gaussian process priors | |
| Freestone et al. | Patient-specific neural mass modeling-stochastic and deterministic methods | |
| US11497429B2 (en) | Iterative process for calibrating a direct neural interface | |
| JP3139414B2 (ja) | 生体内活動部位推定方法及び生体内活動部位推定装置 | |
| CN114463493A (zh) | 基于编解码结构的经颅磁刺激电场快速成像方法和模型 | |
| CN115500841A (zh) | 一种融合时域与频域特征深度学习的室性早搏定位方法 | |
| Haneishi et al. | Multiple current dipole estimation using simulated annealing | |
| Wang et al. | Vector-entropy optimization-based neural-network approach to image reconstruction from projections | |
| JP3156772B2 (ja) | 生体内部活動領域推定方法、装置及びその記録媒体 | |
| CN107909653B (zh) | 一种基于稀疏主成分分析的心脏软组织三维重建方法 | |
| Humber et al. | Nonsmooth formulation of the support vector machine for a neural decoding problem | |
| CN114881930A (zh) | 基于降维定位的3d目标检测方法、装置、设备和存储介质 | |
| US6226544B1 (en) | Living body internal active source estimation apparatus | |
| Irimia | Electric field and potential calculation for a bioelectric current dipole in an ellipsoid | |
| Antigoni et al. | Non linear Hebbian learning techniques and fuzzy cognitive maps in modeling the Parkinson's disease | |
| Wang et al. | Multiobjective neural network for image reconstruction | |
| US20210128002A1 (en) | Dipole group quantification method and dipole group display system |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| FPAY | Renewal fee payment (prs date is renewal date of database) |
Free format text: PAYMENT UNTIL: 20071215 Year of fee payment: 7 |
|
| FPAY | Renewal fee payment (prs date is renewal date of database) |
Free format text: PAYMENT UNTIL: 20081215 Year of fee payment: 8 |
|
| FPAY | Renewal fee payment (prs date is renewal date of database) |
Free format text: PAYMENT UNTIL: 20091215 Year of fee payment: 9 |
|
| FPAY | Renewal fee payment (prs date is renewal date of database) |
Free format text: PAYMENT UNTIL: 20091215 Year of fee payment: 9 |
|
| FPAY | Renewal fee payment (prs date is renewal date of database) |
Free format text: PAYMENT UNTIL: 20101215 Year of fee payment: 10 |
|
| LAPS | Cancellation because of no payment of annual fees |