(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) で表される。
【０００６】
【数１】
0004000343.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) で表される３ｎ×ｍ個の要素をもった行列（リードフィールド行列）である。
【０００９】
【数２】
0004000343.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 に及ぼす磁界〔Ｂ〕との２乗誤差｜〔Ｂｄ〕−〔Ｂ〕｜に、以下に示すペナルティ項を加算した形で表された評価関数ｆを最小にするという条件を付加することにより、線形最小２乗法を用いて電流源〔Ｐ〕を算出する。この評価関数を式(7) に示す。
【００２９】
【数３】
0004000343.tif 000004
【００３０】
評価関数の第１項が磁場の２乗誤差、第２項がペナルティ項である。ペナルティ項中のλはペナルティ項の重み係数、ｗｉは磁気センサから磁場測定面までの距離による影響をキャンセルするための係数である。ペナルティ項は、電流源が固まって存在する程、その値が小さくなるという性質があるので、評価関数中の第１項の磁場の２乗誤差のみを最小にするという条件で電流源を求めた場合に、電流源がバラツキ易いという傾向が抑制されるとともに、離散的に存在するノイズ成分が解として採用されることが軽減される。
【００３１】
上記の評価関数ｆを最小にするという条件を付加することにより、式(4) は
〔Ｐ〕＝Ａ+ 〔Ｂｄ〕                      ………(8)
で表される。
ここで、Ａ+ は次式(9) により求められる。
Ａ+ ＝（Ａt Ａ＋λ・Ｗ）-1  ・Ａt                 ………(9)
(9) 式中のＷは、次式(10)で表されるペナルティ関数行列である。ここでは電流源の要素数が２なのでＷは６４行６４列の行列である。
【００３２】
【数４】
0004000343.tif 000005
【００３３】
ここで、一般に、電流源Ｐｉが磁場測定面に近いほど大きな磁場が測定されるので、行列Ｗを次式(11)のようにすることで、求める電流源Ｐｉについて、磁場測定面との距離による影響（すなわち、電流源Ｐｉが磁場測定面の近くに推定されるという傾向）をキャンセルすることができる。
【００３４】
【数５】
0004000343.tif 000006
【００３５】
また、ペナルティ項の重みλを次式(12)のようにすることで、磁場の２乗誤差の項（Ａt Ａ）とペナルティ項（Ｗ）を同程度の大きさにとることができる。
λ＝｜Ａt Ａ｜／｜Ｗ｜                  ………(12)
【００３６】
上記した(8) 式と(9) 式により電流源〔Ｐ〕が求まると、ステップＳ３に進んで、推定された電流源〔Ｐ〕が磁気センサＳ1 〜Ｓm に及ぼす磁界〔Ｂ〕と測定された磁界〔Ｂｄ〕との２乗誤差が大域的に最小であるか否かを判断する。ここで、磁場の２乗誤差が大域的に最小であるとは、後述するステップＳ４を繰り返す過程で格子点群を複数回移動させた場合に、それぞれの格子点の位置で上述した方法により求めた磁場の２乗誤差の中で、値が最小のものをいう。磁場の２乗誤差が大域的に最小であるか否かの判断は、後述するステップＳ４を繰り返す過程で、再配置された格子点群について、上述したステップＳ２、Ｓ３の処理により、それぞれ求めた磁場の２乗誤差を記憶しておき、それぞれの値を比較して極小となる値を大域的な最小値とすればよい。
【００３７】
最小でないと判断された場合はステップＳ４に進み、格子点数３２個の３次元格子点群の再配置を行う。ステップＳ４では、最小化問題解決に用いられる直接探索法のひとつであるシンプレックス法に基づき、格子点群の再配置が行われることになる。通常、シンプレックス法の場合、ｎ次元空間の（ｎ＋１）個の点の集合が最小化問題対象のシンプレックス（単体）を形成する。ここでは格子点群が３次元であるから、１個のシンプレックスを構成する格子点の数は４個であることになる。したがって、ここでは３２個の格子点を三角錐（４面体）の頂点で８個にグループ分けして８個のシンプレックスＳＰ１〜ＳＰ８を設定する。
【００３８】
そして、設定したシンプレックスＳＰ１〜ＳＰ８の各々について、シンプレックス法のアルゴリズムにおける鏡像・拡張・収縮・縮小の各操作に従って、各グループ毎に推定電流の最も小さな格子点だけを動かして移動後の格子点の電流源の大きさを予測し、その結果から電流源が大きくなると期待される位置を格子点の再配置位置として決定することになる。つまり、４個の格子点のうち推定電流の最も小さな格子点をシンプレックス法のアルゴリズムに従って動かして移動後に推定電流が大きくなるか否かを予測算出し、この予測結果に従って電流源が大きくなると期待される位置を決定するのである。
【００３９】
一方、シンプレックス法のアルゴリズムにおける各操作には目的関数が必要であるが、移動後の電流源Ｐ' ＝（Ｐsa，Ｐsb）の大きさの逆数であるｆｓ＝〔√（Ｐsa2 ，Ｐsb2 ）〕-1が目的関数ｆｓとなる。ここでは電流源は要素２のベクトルであり、電流源が大きくなる位置を求めたいわけであるから、目的関数ｆｓの値が小さくなるほど電流源が大きくなるよう、電流源Ｐ' （Ｐsa，Ｐsb）の大きさの逆数を目的関数ｆｓにしている。
【００４０】
以下では、シンプレックスＳＰ１の４個の格子点の位置，電流源およびその目的関数ｆｓの値は、格子点の位置r1（電流源Ｐs1，目的関数値ｆs1），格子点の位置r2（電流源Ｐs2，目的関数値ｆs2），格子点の位置r3（電流源Ｐs3，目的関数値ｆs3），格子点の位置r4（電流源Ｐs4，目的関数値ｆs4）とし、かつｆs1＞ｆs2＞ｆs3＞ｆs4であるとして説明する。つまり、格子点の位置r1の電流源Ｐs1が最も小さく、格子点の位置r4の電流源Ｐs4が最も大きいものとする。
【００４１】
また、実施例の場合、電流源Ｐ' は次の式(13)により近似算出される。
Ｐ' ＝（Ａ' t Ａ' ＋λＷ')-1Ａ' t （Ｂｄ−Ａ''Ｐ''）    ……(13)
〔但し，Ａ' は移動電流源のリードフィールド行列（１２９行２列），
Ａ''は非移動電流源のリードフィールド行列（１２９行６２列），
Ｐ''は非移動電流源の大きさベクトル（要素数６２）
Ｗ' は移動電流源に対するペナルティ関数行列（２行２列），
Ｂｄは計測磁場ベクトル（要素数１２９），λは重み係数〕
(9) 式の場合は６４行６４列の逆行列を計算する必要があるが、(13)式の場合は２行２列の逆行列を計算であるので計算量が少ない。なお、移動後の位置が、Ａ' とＷ' に因子として入る。他は全て既知である。
【００４２】
続いて、鏡像・拡張・収縮・縮小の各操作に従って各シンプレックスにおける格子点の再配置位置が決定される過程を、図３および図４を参照しながら説明する。図３は実施例での格子点群の再配置の求出経過の際の電流源の移動位置を示す模式図、図４は実施例での格子点群再配置位置の求出経過を示すフローチャートである。
〔ステップＦ１〕最も大きな目的関数値ｆs1（最も小さな電流源Ｐs1）をもつ格子点の位置r1の鏡像の位置raを求めるとともに、鏡像の位置raの目的関数値ｆraを求める。鏡像の位置raは、図３に示すように、位置r1から他の３つの格子点の作る三角形への垂線の上であって、三角形に対する距離の等しい位置である。
【００４３】
〔ステップＦ２〕鏡像の位置raの目的関数値ｆraが、２番目に大きな目的関数値ｆs2以下か否かを調べて、以下であれば次のステップに進み、否であればステップＦ６に進む。
【００４４】
〔ステップＦ３〕鏡像の位置raの目的関数値ｆraが、最小目的関数値ｆs4未満か否かを調べて、未満であれば、鏡像の位置raを格子点の再配置位置として確定する。否であれば、次のステップＦ４に進む。
【００４５】
〔ステップＦ４〕拡張の位置rbを求めるとともに，拡張の位置rbの目的関数値ｆrbを求める。拡張の位置rbは、図３に示すように、位置r1から他の３つの格子点の作る三角形への垂線の上であって、鏡像の位置raより適当な距離だけ先の位置である。
【００４６】
〔ステップＦ５〕拡張の位置rbの目的関数値ｆrbが、鏡像の位置raの目的関数値ｆra未満か否かをチェックする。未満であれば、拡張の位置rbを格子点の再配置位置として確定する。否であれば、鏡像の位置raを格子点の再配置位置として確定する。
【００４７】
〔ステップＦ６〕前記ステップＦ２で、鏡像の位置raの目的関数値ｆraが２番目に大きな目的関数値ｆs2よりも大きいと判断された場合は、鏡像の位置raの目的関数値ｆraが、最大目的関数値ｆs1未満か否かをチェックする。未満であれば、鏡像の位置raを格子点の再配置位置として確定する。否であれば、次のステップＦ７に進む。
【００４８】
〔ステップＦ７〕収縮の位置rcを求めるとともに、収縮の位置rcの目的関数値ｆrcを求める。収縮の位置rcは、図３に示すように、位置r1から他の３つの格子点の作る三角形への垂線の上であって、三角形から適当な距離だけ位置r1側へ近寄った位置である。
【００４９】
〔ステップＦ８〕収縮の位置rcの目的関数値ｆrcが、最大目的関数値ｆs1未満か否かをチェックする。未満であれば、収縮の位置rcを格子点の再配置位置として確定する。否であれば、次のステップＦ９に進む。
【００５０】
〔ステップＦ９〕縮小操作により、位置r1と位置r4の中間位置を位置r1ｍに対応する格子点の再配置位置とし、位置r2と位置r4の中間位置r2ｍを位置r2に対応する格子点の再配置位置とし、位置r3と位置r4の中間位置r3ｍを位置r3に対応する格子点の再配置位置として確定することになる。
以上で、シンプレックスＳＰ１に対する処理が終わることになる。
【００５１】
さらに、残りの７個のシンプレックスＳＰ２〜ＳＰ８についてもシンプレックスＳＰ１と同様の処理が行われる。したがって、ステップＳ４が１回行われると、各シンプレックス毎に１個、合計８個の電流源について処理が行われることになる。
【００５２】
格子点を移動させた後、ステップＳ２に戻って、移動させた格子点群について、上述した評価関数を最小にするという条件を付加して、線形最小２乗法で電流源〔Ｐ〕を新たに求める。電流源Ｐ' は近似的に求めた値であるので、格子点群の再配置後、正確に電流源を算出し直すことになるのである。そして、ステップＳ３で再び２乗誤差が大域的に最小であるか否かの判断を行い、最小でないと判断した場合は上述したステップＳ２〜Ｓ４を繰り返し行う。
【００５３】
ステップＳ３で最小であると判断された場合は、ステップＳ５において大域的に最小となった磁界〔Ｂ〕に対する電流源〔Ｐ〕を真の電流源と推定する。真の電流源として推定された電流源〔Ｐ〕は、図１に示した光磁気ディスク９に記憶されたＸ線ＣＴ像やＭＲ断層像上に重ね合わされて、例えばカラーモニタ１０に表示される。
【００５４】
〔シミュレーション〕
実施例の生体活動電流源推定装置による電流源推定機能を視覚的に確認するためのシミュレーション例を図面を参照しながら説明する。図５はシミュレーションの状況示す図であって、白抜きのバツマークが示す位置に真の電流源（電流双極子）が置かれる。計測磁場の信号対雑音比であるＳ／Ｎ比は１０とした。
図５（ａ）は、格子点群の再配置を行わない最初の電流源の算出状況を示している。第１回目の格子点群の再配置後は、図５（ｂ）に示すように、真の電流源に算出電流源が若干近くなり、第２回目の格子点群の再配置後は、図５（ｃ）に示すように、真の電流源に算出電流源がぐんと近くなり、第３回目の格子点群の再配置後は、図５（ｄ）に示すように、真の電流源に算出電流源が殆ど一致しており、実施例の生体活動電流源推定装置での格子点群の再配置が的確で、良好な推定結果が迅速に得られることが良く分かる。
【００５５】
この発明は上記の実施例に限られるものではなく、以下のように変形して実施することができる。
（１）実施例の場合、１個のシンプレックスに対し、鏡像・拡張・収縮・縮小の操作処理が１回であったが、１個のシンプレックスに対し、鏡像・拡張・収縮・縮小の操作処理が所定の停止条件を満たすまで繰り返し行われる構成の装置が、変形例として挙げられる。
【００５６】
（２）この発明の生体活動電流源推定装置における格子点群の再配置に用いられるシンプレックス法の具体的なアルゴリズムは、実施例に示したアルゴリズムに限られるものではない。
【００５７】
【発明の効果】
以上の説明から明らかなように、この発明の生体活動電流源推定装置によれば、電流源の算出を繰り返す前に行う格子点群の再配置は、いわゆる直接探索法のひとつであるシンプレックス法のアルゴリズムに基づき格子点群がグループ分けされるとともに、シンプレックス法のアルゴリズムに従って格子点群の各グループに対する鏡像・拡張・収縮・縮小の各操作により、推定電流の最も小さな格子点の試験的な移動に伴う電流源の大きさの予測結果に従って各グループ毎に格子点の再配置用の位置が決まるよう構成されているので、特定の推定電流の大きな電流源に他の電流源が画一的に集中する事態が避けられると同時に、従来のように単純なルールによる一律的な移動ではなく試験的な予測に基づいて再配置用の位置が決められることから、真の電流源の状況を十分に反映した的確な格子点群の再配置がなされるようになる。
その結果、診断対象領域における真の電流源が、同程度の大きさの複数個の電流源からなる場合や、広範囲に電流源が分布する非集中的な電流源の場合であっても、格子点群の再配置が的確に行われるので、電流源算出の繰り返し回数が少なくなり、良好な推定結果を迅速に得ることができる。
【図面の簡単な説明】
【図１】この発明に係る生体活動電流源推定装置の一実施例の概略構成を示したブロック図である。
【図２】実施例に係る推定処理全体のフローチャートである。
【図３】実施例での格子点群の再配置の求出経過の際の電流源の移動位置を示す模式図である。
【図４】実施例での格子点群再配置位置の求出経過を示すフローチャートである。
【図５】実施例に係る推定処理のシミュレーション結果を示すグラフである。
【図６】従来例に係る推定処理の説明図である。
【符号の説明】
１…マルチチャンネルＳＱＵＩＤセンサ
２…磁気シールドルーム
４…データ変換ユニット
５…データ収集ユニット
８…データ解析ユニット
１０…カラーモニタ
Ｍ…被検体
Ｓ1 〜Ｓm …磁気センサ
【図１】
0004000343.tif 000007
【図２】
0004000343.tif 000008
【図３】
0004000343.tif 000009
【図４】
0004000343.tif 000010
【図５】
0004000343.tif 000011
【図６】
0004000343.tif 000012