(57)【特許請求の範囲】
【請求項１】
生体活動電流によって生じる磁界を複数個の磁気センサで検出することによって、生体活動電流源の位置，大きさ，方向等の物理量を推定する生体活動電流源推定装置であって、（ａ）被検体の診断対象領域に近接配備され、前記診断対象領域内の生体活動電流源による微小磁界を計測する複数個の磁気センサと、（ｂ）前記各磁気センサによって計測された磁界データをデジタルデータに変換するデータ変換手段と、（ｃ）前記デジタルデータに変換された磁界データを収集して記憶するデータ収集手段と、（ｄ）前記診断対象領域に複数個の格子点を設定する格子点設定手段と、（ｅ）前記各格子点上の未知の電流源が及ぼす磁界と前記データ収集手段に記憶された磁界データの２乗誤差と、電流源の重み係数λ付き２乗和との和を最小にするという条件を付加することにより未知の電流源を求める電流源算出手段と、（ｆ）前記電流源算出手段で推定された電流源を、前記被検体の診断対象領域の断層像に重ね合わせて表示する表示手段とを備えているとともに、（ｇ）格子点上の未知の電流源の位置と真の電流源の位置との間の位置誤差が未知の電流源の発生磁界へ与える影響の度合い（磁界影響度）を求出する磁界影響度求出手段と、（ｈ）前記磁界影響度求出手段による磁界影響度を加味して前記重み係数λを求出する重み係数求出手段とを備えており、かつ前記電流源算出手段は重み係数求出手段によって求出された重み係数λを用いて電流源を求めるよう構成されていることを特徴とする生体活動電流源推定装置。
【請求項２】
請求項１に記載の生体活動電流源推定装置において、重み係数求出手段は、重み係数λを求出する際に磁気センサによって計測された磁界データをも加味し、重み係数λが変数として含まれている最適化基準関数ＧＣＶ（λ）＝（磁界データの２乗誤差）／（磁界データの誤差ベクトルの自由度）の値を最小とする時のλの値を、求出すべき重み係数λとするよう構成されている生体活動電流源推定装置。
【発明の詳細な説明】
【０００１】
【産業上の利用分野】
この発明は、生体活動電流源の位置，向き，大きさを推定する生体活動電流源推定装置に係り、生体活動電流源の物理量についての推定精度を向上させるための技術に関する。
【０００２】
【従来の技術】
生体に刺激を与えると、細胞膜を挟んで形成されている分極が壊れて生体活動電流が流れる。この生体活動電流は、脳や心臓において現れ、脳波，心電図として記録される。また、生体活動電流によって生じる磁界は、脳磁図，心磁図として記録される。
【０００３】
近年、生体内の微小な磁界を計測するセンサとして、ＳＱＵＩＤ（Superconduc-ting Quantum Interface Device ：超電導量子干渉計）を用いたセンサが開発されている。このセンサを頭部の外側に置き、脳内に生じた生体活動電流源である電流双極子（以下、単に電流源とも称する）による微小磁界をそのセンサで無侵襲に計測することができる。計測された磁界データから病巣に関連した電流源の位置, 向き, 大きさを推定し、推定した電流源をＸ線ＣＴ装置やＭＲＩ装置で得られた断層像上に表示させて患部等の物理的位置の特定などに用いている。
【０００４】
従来、電流源の推定方法の一つとして、最小ノルム法を用いた手法がある（例えば、W.H.Kullmann, K.D.Jandt, K.Rehm, H.A.Schlitt, W.J.Dallas and W.E.Smith, Advances in Biomagnetism, pp.571-574, Plenum Pless, New York, 1989) 。以下、図６を参照して、最小ノルム法を用いた従来の電流源推定方法を説明する。
図６に示すように、被検体Ｍに近接してマルチチャンネルＳＱＵＩＤセンサ１が配備される。マルチチャンネルＳＱＵＩＤセンサ１は、デュアーと呼ばれる容器内に多数の磁気センサ（ピックアップコイル）Ｓ1 〜Ｓm を液体窒素などの冷媒に浸漬して収納している。
【０００５】
一方、被検体Ｍの診断対象領域である例えば脳に、多数の格子点(1) 〜(n) を設定し、各格子点に未知の電流源（電流双極子）を仮定し、各電流源を３次元ベトクルVPj （ｊ＝１〜ｎ）で表す。そうすると、ＳＱＵＩＤセンサ１の各磁気センＳ1 〜Ｓm で検出される磁界Ｂi 〜Ｂm は、次式(1) で表される。
【０００６】
【数１】
0004006543.tif 000002
【０００７】
式(1) において、VPj ＝（Ｐjx，Ｐjy，Ｐjz）
αij＝（αijx,αijy,αijz ）
で表される。なお、αijは、格子点上にＸ，Ｙ，Ｚ方向の単位大きさの電流源を置いた場合に磁気センサＳ1 〜Ｓm の各位置で検出される磁界の強さを表す既知の係数である。
【０００８】
ここで、〔Ｂ〕＝（Ｂ1 ，Ｂ2 ，…，Ｂm ）
〔Ｐ〕＝（Ｐ1x，Ｐ1y，Ｐ1z，Ｐ2x，Ｐ2y，Ｐ2z，…，Ｐnx，Ｐny，Ｐnz）
のように表すと、(1) 式は(2) 式のような線形の関係式に書き換えられる。
〔Ｂ〕＝Ａ〔Ｐ〕                    ………(2)
(2) 式において、Ａは次式(3) で表される３ｎ×ｍ個の要素をもった行列（リードフィールド行列）である。
【０００９】
【数２】
0004006543.tif 000003
【００１０】
ここで、Ａの逆行列をＡ- で表すと、〔Ｐ〕は次式(4) で表される。
〔Ｐ〕＝Ａ- 〔Ｂ〕                  ………(4)
ここで、最小ノルム法は、式の個数ｍ（磁気センサＳ1 〜Ｓm の個数）よりも、未知数の個数３ｎ（各格子点に仮定される電流源のＸ，Ｙ，Ｚ方向の大きさを考慮した場合の未知数）が多い場合を前提として、電流源〔Ｐ〕のノルム｜〔Ｐ〕｜を最小にするという条件を付加することで電流源〔Ｐ〕の解を求めるものである。なお、上述した式の個数ｍと未知数の個数３ｎとを等しくとることで、解は一意的に求めることができるが、かかる場合には、解が非常に不安定となることからこの最小ノルム法が用いられている。
【００１１】
電流源〔Ｐ〕のノルム｜〔Ｐ〕｜を最小にするという条件を付加することで、上式(4) は次式(5) のように表される。
〔Ｐ〕＝Ａ+ 〔Ｂ〕                  ………(5)
ここで、Ａ+ は次式(6) で表される一般逆行列である。
Ａ+ ＝Ａt （ＡＡt ）-1            ………(6)
ただし、Ａt はＡの転置行列である。
【００１２】
上式(5) を解いて各格子点上の電流源VPj の方向，大きさを推定し、その中で値の最も大きなものを真の電流源に近いものとしている。これが、最小ノルム法による電流源推定方法の原理である。
【００１３】
さらに、最小ノルム法の位置分解能を向上させるために格子点分割を細分しながら最小ノルム解を繰り返し求め、より真の電流源に近い電流源を推定する方法も提案されている。しかし、この場合は(5) 式のベクトル〔Ｐ〕の要素が多くなり、最小ノルム解の計算精度が低下するという難点がある。
そこで、本出願人は、先に特開平６−３４３６１３号及び特開平６−３４３６１４号によって、最小ノルム法によって推定された電流源の内、値の大きな電流源が存在する格子点付近に他の格子点群を移動させて、格子点の数を前回と同じくして、格子点間隔のみを狭くした状態で電流源を推定することで、最小ノルム解の精度を維持しながら、電流源を精度よく推定する格子点移動最小ノルム法を提案している。しかし、この場合も、格子点間隔に基づく収束判定値の適切な値の設定困難等により、安定した推定結果が得られないという難点がある。
【００１４】
そこで、更に本出願人は、特願平６−４７２２０号によって、格子点上の電流源が及ぼす磁界と実際の計測磁界の２乗誤差を最小にするという条件を付加することにより、格子点間隔に基づく収束判定値を不要として、安定した推定結果を得るようにする手法を提案している。
しかし、この手法の場合も、推定された各電流源のモーメントにバラツキが発生し易く、計測磁場の離散的なノイズを解（電流源）として推定してしまう難点があるので、更に出願人は、特開平７−３２７９４３号によって、設定した各格子点上の未知の電流源が及ぼす磁界と実際の計測磁界（磁界データ）の２乗誤差（磁界２乗誤差項）と、電流源の重み係数λ付き２乗和（ペナルティ関数項）との和を最小にするという条件を付加する手法を提案している。
【００１５】
【発明が解決しようとする課題】
しかしながら、上記の特開平７−３２７９４３号に係る従来手法の場合、ペナルティ関数項における重み係数λを適切な値にセットすることが必ずしも容易でないという問題がある。従来、重み係数λは、磁界２乗誤差項とペナルティ関数項が同程度の比重となるような値にセットするか、格子間隔に応じた値にセットしているのであるが、格子点の設定位置に誤差がある場合や、計測磁場に含まれるノイズの大きさによっては、重み係数λが不適切な値となり、精度の十分な推定結果が得られない場合がある。
【００１６】
この発明は、このような事情に鑑みてなされたものであって、ペナルティ関数項における重み係数λが適切な値にセットされて、生体活動電流源の物理量について精度の十分な推定結果を得ることができる生体活動電流源推定装置を提供することを目的としている。
【００１７】
【課題を解決するための手段】
上記の課題を解決するために、請求項１の発明は、生体活動電流によって生じる磁界を複数個の磁気センサで検出することによって、生体活動電流源の位置，大きさ，方向等の物理量を推定する生体活動電流源推定装置であって、（ａ）被検体の診断対象領域に近接配備され、前記診断対象領域内の生体活動電流源による微小磁界を計測する複数個の磁気センサと、（ｂ）前記各磁気センサによって計測された磁界データをデジタルデータに変換するデータ変換手段と、（ｃ）前記デジタルデータに変換された磁界データを収集して記憶するデータ収集手段と、（ｄ）前記診断対象領域に複数個の格子点を設定する格子点設定手段と、（ｅ）前記各格子点上の未知の電流源が及ぼす磁界と前記データ収集手段に記憶された磁界データの２乗誤差と、電流源の重み係数λ付き２乗和との和を最小にするという条件を付加することにより未知の電流源を求める電流源算出手段と、（ｆ）前記電流源算出手段で推定された電流源を、前記被検体の診断対象領域の断層像に重ね合わせて表示する表示手段とを備えているとともに、（ｇ）格子点上の未知の電流源の位置と真の電流源の位置との間の位置誤差が未知の電流源の発生磁界へ与える影響の度合い（磁界影響度）を求出する磁界影響度求出手段と、（ｈ）前記磁界影響度求出手段による磁界影響度を加味して前記重み係数λを求出する重み係数求出手段とを備えており、かつ前記電流源算出手段は重み係数求出手段によって求出された重み係数λを用いて電流源を求めるよう構成されている。
【００１８】
また、請求項２の発明は、請求項１に記載の生体活動電流源推定装置において、重み係数求出手段は、重み係数λを求出する際に磁気センサによって計測された磁界データをも加味し、重み係数λが変数として含まれている最適化基準関数ＧＣＶ（λ）＝（磁界データの２乗誤差）／（磁界データの誤差ベクトルの自由度）の値を最小とする時のλの値を、求出すべき重み係数λとするよう構成されている。
【００１９】
〔作用〕
この発明の生体活動電流源推定装置の作用は次のとおりである。
請求項１の発明の装置では、各磁気センサで計測された微小磁界はデータ変換手段でデジタルデータに変換された後、データ収集手段に記憶される。そして、格子点設定手段により、診断対象領域に複数個の格子点が仮想的に設定される。普通、この格子点の個数は、磁気セサンの個数よりも少なく設定される。一方、磁界影響度求出手段により、各格子点上の未知の電流源の位置と真の電流源の位置との間の位置誤差が未知の電流源の発生磁界へ与える影響の度合い（磁界影響度）が求出されるとともに、重み係数求出手段により磁界影響度を加味した重み係数λが求出される。
続いて、電流源算出手段は、求出した重み係数λを用い、各格子点上の未知の電流源が及ぼす磁界とデータ収集手段に記憶された磁界データとの２乗誤差と、電流源の重み係数λ付き２乗和との和を最小にするという条件を付加して、各格子点上の電流源を求める。その後、推定結果としての電流源は、表示手段により被検体の診断対象領域の断層像に重ね合わせて表示される。
【００２０】
請求項１の発明の装置の場合、電流源算出過程でのペナルティ関数項における重み係数λは、従来の場合ように、磁界２乗誤差項とペナルティ関数項とが同程度の比重となるか、単に格子間隔に対応するといった単純な値ではなく、重み係数λは、未知の電流源と真の電流源との間の位置誤差に起因する磁界影響度が加味された従来より適切な値であるので、生体活動電流源についての推定精度が向上する。
【００２１】
請求項２の発明の生体活動電流源推定装置では、重み係数求出手段によって、磁界影響度だけでなく磁界データも加味され、最適化基準関数ＧＣＶ（λ）＝（磁界データの２乗誤差）／（磁界データの誤差ベクトルの自由度）の値を最小とする時のλの値が、求出すべき重み係数λであるとして決定される。
つまり、決定される重み係数λは、実測される磁界データの中のノイズの影響も考慮されることになり、しかも、最適化基準関数ＧＣＶ（λ）の値が最小であるλは、磁界データの誤差の大小だけでなく、磁界データの誤差ベクトルの自由度（誤差の方向）も考慮されている磁界データのノイズの状況に良くマッチしたものなので、重み係数λ付きのペナルティ関数項が磁界データのノイズも十分に反映した的確なものになる結果、実測された磁界データに磁界ノイズを含んでいても、十分な推定精度が得られる。
【００２２】
【実施例】
以下、図面を参照してこの発明の実施例を説明する。
図１はこの発明に係る生体活動電流源推定装置の一実施例の概略構成を示したブロック図である。図中、符号２は磁気シールドルームであり、この磁気シールドルーム２内に被検体Ｍが仰臥されるベッド３と、被検体Ｍの例えば脳に近接配備され、脳内に生じた生体活動電流源による微小磁界を無侵襲に計測するためのマルチチャンネルＳＱＵＩＤセンサ１とが設けられている。上述したように、マルチチャンネルＳＱＵＩＤセンサ１は、デュアー内に多数の磁気センサが冷媒に浸漬して収納されている。本実施例において、各磁気センサは被検体Ｍの脳を球体と仮定した場合に、その半径方向の磁界成分を検出する一対のコイルでそれぞれ構成されている。
【００２３】
マルチチャンネルＳＱＵＩＤセンサ１で検出された磁界データはデータ変換ユニット４に与えられてデジタルデータに変換された後、データ収集ユニット５に集められる。刺激装置６は、被検体Ｍに電気的刺激（あるいは音、光刺激など）を与えるためのものである。ポジショニングユニット７は、マルチチャンネルＳＱＵＩＤセンサ１を基準とした３次元座標系に対する被検体Ｍの位置関係を把握するための装置である。例えば、被検体Ｍの複数個所に小コイルを取り付け、これらの小コイルにポジショニングユニット７から給電する。そして、各コイルから発生した磁界をマルチチャンネルＳＱＵＩＤセンサ１で検出することにより、マルチチャンネルＳＱＵＩＤセンサ１に対する被検体Ｍの位置関係を把握する。なお、ＳＱＵＩＤセンサ１に対する被検体Ｍの位置関係を把握するための手法は、これ以外に、デュワーに投光器を取り付けて光ビームを被検体Ｍに照射して両者の位置関係を把握するものや、あるいは、特開平５─２３７０６５号、特開平６─７８８９２５号などに開示された種々の手法が用いられる。
【００２４】
データ解析ユニット８は、データ収集ユニット５に集められた磁界データに基づいて、被検体Ｍの診断対象領域内の電流源を推定するためのものである。このデータ解析ユニット８は、後述する説明から明らかになるように、この発明における格子点設定手段、電流源算出手段、磁界影響度求出手段、重み係数求出手段および他の各手段としての機能を備える。データ解析ユニット８に関連して設けられた光磁気ディスク９には、例えばＸ線ＣＴ装置やＭＲＩ装置で得られた断層画像が記憶されており、データ解析ユニット８で推定された電流源が、これらの断層像上に重ね合わされてカラーモニタ１０に表示されたり、あるいはカラープリンタ１１に印字出力されるようになっている。なお、Ｘ線ＣＴ装置やＭＲＩ装置で得られた断層画像は、図１に示した通信回線１２を介してデータ解析ユニット８に直接伝送するように構成してもよい。
【００２５】
上述したように、マルチチャンネルＳＱＵＩＤセンサ１を基準とした３次元座標系に対する被検体Ｍの位置関係を測定して記憶するとともに、マルチチャンネルＳＱＵＩＤセンサ１で被検体Ｍの診断対象領域である例えば脳内の生体活動電流源からの微小磁界を計測して、その磁界データをデータ収集ユニット５に集めた後、データ解析ユニット８で電流源の推定処理が、以下のようなプロセスで実行される。
【００２６】
〔格子点設定プロセス〕
図６に示した従来例と同様に、診断対象領域である例えば脳内に３次元の格子点群Ｎを想定する。ここで、格子点群Ｎについての未知数（各格子点について想定される電流源の個数）が、マルチチャンネルＳＱＵＩＤセンサ１を構成する各磁気センサＳ1 〜Ｓm の個数ｍよりも少なくなるように、格子点群Ｎの格子点の数が設定される。ここでは、格子点の数は３２個、磁気センサの数は１２９個とする。格子点の設定の仕方は特に限定されないが、臨床的知見に基づく設定アルゴリズムやシンプレックス法などの数学的手法に基づき設定アルゴリズムが挙げられる。
【００２７】
〔リードフィールド行列求出プロセス〕
そして、式(3) で表される行列Ａの各係数をビオ・サバールの法則を使って算出（行列Ａの各係数は、後述する格子点の移動ごとに算出される）する。ここでは、診断対象領域を球体と仮定し、各センサで球体の半径方向成分の磁界を検出しているので、要素数２の電流源が３２個であり、発生磁場が１２９個である。したがって、リードフィールド行列Ａは１２９行６４列となる。
行列Ａが求まると、各格子点上の電流源をペナルティ項付きの線形最小２乗法で求めることになる。
【００２８】
〔電流源算出プロセス〕
【００２９】
上述した式(4) を使って、磁気センサＳ1 〜Ｓm によって計測され、データ収集ユニット５に記憶された磁界データ〔Ｂｄ〕から各格子点上の電流源〔Ｐ〕を求める。式(4) を改めて以下に示す。
〔Ｐ〕＝Ａ- 〔Ｂｄ〕                  ………(4)
【００３０】
行列Ａは上述した(3) 式で表され、今は１２９行×６４列の個数の要素をもった行列である。(4) 式において、式の個数ｍ（磁気センサＳ1 〜Ｓm の個数）が、未知数の個数よりも多いために解が求まらない。そこで、測定された磁界〔Ｂｄ〕と、各格子点上に仮定した電流源〔Ｐ〕が磁気センサＳ1 〜Ｓm に及ぼす磁界（計算によって求められる磁界）〔Ｂ〕との２乗誤差｜〔Ｂｄ〕−〔Ｂ〕｜2 に、以下に示す重み係数λ付きペナルティ項を加算した形で表された評価関数ｆを最小にするという条件を付加することにより、線形最小２乗法を用いて電流源〔Ｐ〕を算出する。この評価関数を式(7) に示す。
【００３１】
【数３】
0004006543.tif 000004
【００３２】
評価関数の第１項が磁場の２乗誤差、第２項がペナルティ項である。ペナルティ項中のλはペナルティ項の重み係数、ｗｉは磁気センサから磁場測定面までの距離による影響をキャンセルするための係数である。ペナルティ項は、電流源が固まって存在する程、その値が小さくなるという性質があり、評価関数中の第１項の磁場の２乗誤差のみを最小にするという条件で電流源を求めた場合に、電流源がバラツキ易いという傾向を抑えるとともに、離散的に存在するノイズ成分が解として採用されることを抑える。
【００３３】
上記の評価関数ｆを最小にするという条件を付加することにより、式(4) は
〔Ｐ〕＝Ａ+ 〔Ｂｄ〕                      ………(8)
で表される。
ここで、Ａ+ は次式(9) により求められる。
Ａ+ ＝（Ａt Ａ＋λ・Ｗ）-1  ・Ａt                 ………(9)
(9) 式中のＷは、次式(10)で表されるペナルティ関数行列である。ここでは、電流源の要素数が２なのでＷは６４行６４列の行列である。(9) は正則化一般逆行列法が適用されており、行列Ａt Ａに行列λ・Ｗを加え正則化すると逆行列を求められる。
【００３４】
【数４】
0004006543.tif 000005
【００３５】
ここで、一般に、電流源Ｐｉが磁場測定面に近いほど大きな磁場が測定されるので、行列Ｗを次式(11)のようにすることで、求める電流源Ｐｉについて、磁場測定面との距離による影響（すなわち、電流源Ｐｉが磁場測定面の近くに推定されるという傾向）をキャンセルすることができる。
そして、実施例の装置の場合、ペナルティ項の重み係数λは次のようにして求出する。
【００３６】
【数５】
0004006543.tif 000006
【００３７】
〔重み係数λの求出プロセス〕
実施例の装置のペナルティ項の重み係数λの求出過程を、図２のフローチャートを参照しながら説明する。
〔ステップＳ１〕先ず格子点上の未知の電流源の位置と真の電流源の位置との間の位置誤差が未知の電流源の発生磁界へ与える影響の度合い（磁界影響度）を示す行列Ａ’を求める。行列Ａ’は１２９行６４列の行列であり、隣接する電流源のリードフィールド行列Ａの要素の差を要素としており、行列Ａ’のｉ番目の列ベクトルα’i ，行列Ａのｉ番目の列ベクトルαi ，行列Ａのｊ番目の列ベクトルαj とすると、α’i ＝｜αi −αj ｜である。ｉ＝１〜６４，ｊ＝１〜６４であるが、ｉとｊは等しくない。
【００３８】
〔ステップＳ２〕下式(12)の最適化基準関数ＧＣＶ（λ）の算出を開始する。
ＧＣＶ（λ）＝
｜ [I −（Ａ＋Ａ’ ）Ａ+ ] Bd｜2 ÷｛trace [I−（Ａ＋Ａ’ )]Bd ｝2 ……(12)
但し：(12)式中のＩは単位行列，trace は行列の対角要素の和である。
Ａ+ は(9) 式で示したようにλを含み、このλの値を変化させてＧＣＶ（λ）の値を次々と算出してゆく。またＧＣＶはGeneralized Cross-Validationの略である。
【００３９】
〔ステップＳ３〕ＧＣＶ（λ）の値の算出結果に基づいてＧＣＶ（λ）を最小とする時のλの値を確定する。
【００４０】
〔ステップＳ４〕ステップＳ３で確定したλの値を、求出すべき重み係数λとして決定する。
【００４１】
最適化基準関数ＧＣＶ（λ）には、行列Ａ’が含まれていることにより、未知の電流源と真の電流源との間の位置誤差（モデル化ノイズ）に起因する磁界影響度が加味され、従来より適切なλが求められる。また、最適化基準関数ＧＣＶ（λ）には、Ｂｄが含まれていて、実測される磁界データの中のノイズ（実際のノイズ）の影響も考慮されることになる。
しかも、最適化基準関数ＧＣＶ（λ）の値が最小であるλは、磁界の２乗誤差に対応する｜ [I −（Ａ＋Ａ’ ）Ａ+ ] Bd｜2 による磁界データの誤差の大小だけでなく、磁界データの誤差ベクトル、つまり〔（Ｂｄｉ−Ｂｉ）の１２９個の要素（ｉ＝１〜１２９）を持つベクトル〕に対応する｛trace [I−（Ａ＋Ａ’ )]Bd ｝2 により、磁界データの誤差ベクトルの自由度（誤差の方向）も考慮された磁界のノイズ状況に良くマッチした的確なものになる。
そして、こうして決定した重み係数λと上記した(8) 式と(9) 式により電流源〔Ｐ〕が求められる。
【００４２】
〔算出電流源の判定プロセス〕
求まった電流源〔Ｐ〕が適切か否かを次のようにチェックする。すなわち、推定された電流源〔Ｐ〕が磁気センサＳ1 〜Ｓm に及ぼす磁界〔Ｂ〕と測定された磁界〔Ｂｄ〕との２乗誤差と、上述した重み係数λ付きペナルティ項との和が、大域的に最小であるか否かを判断する。ここで、磁場の２乗誤差とペナルティ項との和が大域的に最小であるとは、後述する反復プロセスを繰り返す過程で格子点群を複数回移動させた場合に、それぞれの格子点の位置で上述した方法により求めた磁場の２乗誤差とペナルティ項との和の中で、値が最小のものをいう。磁場の２乗誤差と重み係数λ付きペナルティ項との和が大域的に最小であるか否かの判断は、後述する反復プロセスを繰り返す過程で、再配置された格子点群について、上記の過程中の必要プロセスに従って求めた磁場の２乗誤差と重み係数λ付きペナルティ項との和を記憶しておき、それぞれの値を比較して極小となる値を大域的な最小値とすればよい。
【００４３】
〔反復プロセス〕
最小でないと判断された場合は、次のようにして反復プロセスを行う。前回求めた電流源〔Ｐ〕の中、値の大きな電流源が存在する格子点の周りに他の格子点上の電流源を移動させる。これにより、格子点数は元の格子点群Ｎと同じで、間隔が密になった新たな格子点群が得られる。なお、値の大きな電流源は必ずしも一つではなく、通常は複数個ある。
【００４４】
値の大きな電流源が存在する格子点に、他の格子点群を近づけるための手法は特に限定しないが、例えば、次のような手法が例示される。前回求められた各格子点の電流源の大きさを質量と考え、重力によって各格子点間に引力が働くと仮定する。そうすると、各格子点は、質量の大きな格子点に近づいて行き、質量の大きな格子点に近い程、密度の高い格子点群が得られる。各格子点の移動距離は適宜に設定される。さらに、別の格子点移動法としては、例えば、推定電流の大きい格子点を中心とする立方体の頂点に推定電流の小さい各格子点８個を移動する。このとき、立方体の大きさは、例えば、推定電流の大きい格子点と、これに最も近い格子点との距離の半分とする。
【００４５】
格子点を移動させた後、移動させた格子点群について、再び、電流源〔Ｐ〕を新たに求めるとともに、２乗誤差と重み係数λ付きペナルティ項との和が大域的に最小であるか否かの判断を行い、最小でないと判断した場合は、また反復プロセスを繰り返し行う。
【００４６】
〔電流源の推定・表示プロセス〕
２乗誤差と重み係数λ付きペナルティ項との和が大域的に最小であると判断された場合は、判定プロセスにおいて大域的に最小となった磁界〔Ｂ〕に対する電流源〔Ｐ〕を真の電流源と推定する。真の電流源として推定された電流源〔Ｐ〕は、図１に示した光磁気ディスク９に記憶されたＸ線ＣＴ像やＭＲ断層像上に重ね合わされて、例えばカラーモニタ１０に表示される。
【００４７】
〔シミュレーション〕
実施例の生体活動電流源推定装置による電流源推定の精度向上を視覚的に確認するためのシミュレーション例を図面を参照しながら説明する。図３はシミュレーション結果を示す図であって、図３に白抜きのバツマークＱＲが示す位置に真の電流源（電流双極子）が置かれる。
第１回目のシミュレーションでは、計測磁場の信号対雑音比であるＳ／Ｎ比は５とし、反復プロセスは行わず、最初の算出結果を推定結果として、推定結果の各電流源を図３に実線矢印で示した。矢印が長いほど大きな電流源であることを示している。第１回目のシミュレーションの場合、重み係数λは９．１８２８×１０-5となった。
【００４８】
また、Ｓ／Ｎ比を１０とした他は、第１回目と同様にして、第２回目のシュミレーションを行い、推定結果の各電流源を図４に実線矢印で示した。なお、重み係数λはλ＝９．１３５８９×１０-5となった。
さらに、Ｓ／Ｎ比を２０とした他は、第１回目と同様にして、第３回目のシュミレーションを行い、推定結果の各電流源を図５に実線矢印で示した。なお、重み係数λは９．０８７４３×１０-5となった。
【００４９】
図３〜図５に示された電流源の状況を観察すれば、３回のシミュレーションは全て、推定結果の各電流源のうちバツマークＱＲの示す真の電流源の近傍に位置する電流源が大きくて、推定精度が良いことが分かる。また、ＳＮ比が５〜２０と変化しても、推定結果は変わらないので、磁界ノイズの多い診断対象領域の場合にも、十分な推定精度が得られることが分かる。
【００５０】
この発明は上記の実施例に限られるものではなく、以下のように変形して実施することができる。
（１）実施例装置の場合、格子点を移動させて電流源の算出を繰り返す反復プロセスを実行できる構成であったが、反復プロセスを実行しない構成の装置が、変形例として挙げられる。診断対象領域内の真の電流源の位置が予め把握可能な場合、変形例の装置によって最初の電流源の算出で良好な推定結果が反復プロセスを実行せずに出せる。
【００５１】
（２）実施例装置の場合、重み係数求出手段は、磁界影響度と実測磁界データの両方を加味して重み係数λを求出する構成であったが、重み係数求出手段が磁界影響度だけを加味して重み係数λを求出する構成の装置が、変形例として挙げられる。
【００５２】
【発明の効果】
以上の説明から明らかなように、請求項１の発明の生体活動電流源推定装置によれば、電流源算出過程で用いられる重み係数λ付きペナルティ関数項の重み係数λが、未知の電流源と真の電流源との間の位置誤差に起因する磁界影響度が加味されて決定されるので、生体活動電流源についての推定精度が向上する。
【００５３】
また、請求項２の発明の生体活動電流源推定装置によれば、ペナルティ関数項の重み係数λは、実測磁界データ中のノイズの影響も考慮されていて、しかも、磁界データの誤差の大小だけでなく、磁界データの誤差ベクトルの自由度（誤差の方向）も考慮され、磁界データのノイズ状況に良くマッチしたものであるので、磁界データのノイズも十分に反映したより的確なものになり、計測された磁界データに磁界ノイズが含まれていても、十分な推定精度が得られる。
【図面の簡単な説明】
【図１】実施例に係る生体活動電流源推定装置の概略構成を示したブロック図である。
【図２】実施例装置における重み係数求出プロセスを示すフローチャートである。
【図３】実施例に係る第１回目の推定シミュレーション結果を示すグラフである。
【図４】実施例に係る第２回目の推定シミュレーション結果を示すグラフである。
【図５】実施例に係る第３回目の推定シミュレーション結果を示すグラフである。
【図６】従来例に係る推定処理の説明図である。
【符号の説明】
１…マルチチャンネルＳＱＵＩＤセンサ
２…磁気シールドルーム
４…データ変換ユニット
５…データ収集ユニット
８…データ解析ユニット
１０…カラーモニタ
Ｍ…被検体
Ｓ1 〜Ｓm …磁気センサ
【図１】
0004006543.tif 000007
【図２】
0004006543.tif 000008
【図３】
0004006543.tif 000009
【図４】
0004006543.tif 000010
【図５】
0004006543.tif 000011
【図６】
0004006543.tif 000012