JP4006543B2 - Bioactive current source estimation device - Google Patents

Description

【０００７】
式(1) において、VP_j ＝（Ｐ_jx，Ｐ_jy，Ｐ_jz）
α_ij＝（α_ijx,α_ijy,α_ijz ）
で表される。なお、α_ijは、格子点上にＸ，Ｙ，Ｚ方向の単位大きさの電流源を置いた場合に磁気センサＳ₁ 〜Ｓ_m の各位置で検出される磁界の強さを表す既知の係数である。
【０００８】
ここで、〔Ｂ〕＝（Ｂ₁ ，Ｂ₂ ，…，Ｂ_m ）
〔Ｐ〕＝（Ｐ_1x，Ｐ_1y，Ｐ_1z，Ｐ_2x，Ｐ_2y，Ｐ_2z，…，Ｐ_nx，Ｐ_ny，Ｐ_nz）
のように表すと、(1) 式は(2) 式のような線形の関係式に書き換えられる。
〔Ｂ〕＝Ａ〔Ｐ〕 ………(2)
(2) 式において、Ａは次式(3) で表される３ｎ×ｍ個の要素をもった行列（リードフィールド行列）である。
【０００９】
【数２】

【００１０】
ここで、Ａの逆行列をＡ^- で表すと、〔Ｐ〕は次式(4) で表される。
〔Ｐ〕＝Ａ^- 〔Ｂ〕 ………(4)
ここで、最小ノルム法は、式の個数ｍ（磁気センサＳ₁ 〜Ｓ_m の個数）よりも、未知数の個数３ｎ（各格子点に仮定される電流源のＸ，Ｙ，Ｚ方向の大きさを考慮した場合の未知数）が多い場合を前提として、電流源〔Ｐ〕のノルム｜〔Ｐ〕｜を最小にするという条件を付加することで電流源〔Ｐ〕の解を求めるものである。なお、上述した式の個数ｍと未知数の個数３ｎとを等しくとることで、解は一意的に求めることができるが、かかる場合には、解が非常に不安定となることからこの最小ノルム法が用いられている。
【００１１】
電流源〔Ｐ〕のノルム｜〔Ｐ〕｜を最小にするという条件を付加することで、上式(4) は次式(5) のように表される。
〔Ｐ〕＝Ａ⁺ 〔Ｂ〕 ………(5)
ここで、Ａ⁺ は次式(6) で表される一般逆行列である。
Ａ⁺ ＝Ａ^t （ＡＡ^t ）^-1 ………(6)
ただし、Ａ^t はＡの転置行列である。
【００１２】
上式(5) を解いて各格子点上の電流源VP_j の方向，大きさを推定し、その中で値の最も大きなものを真の電流源に近いものとしている。これが、最小ノルム法による電流源推定方法の原理である。
【００１３】
さらに、最小ノルム法の位置分解能を向上させるために格子点分割を細分しながら最小ノルム解を繰り返し求め、より真の電流源に近い電流源を推定する方法も提案されている。しかし、この場合は(5) 式のベクトル〔Ｐ〕の要素が多くなり、最小ノルム解の計算精度が低下するという難点がある。
そこで、本出願人は、先に特開平６−３４３６１３号及び特開平６−３４３６１４号によって、最小ノルム法によって推定された電流源の内、値の大きな電流源が存在する格子点付近に他の格子点群を移動させて、格子点の数を前回と同じくして、格子点間隔のみを狭くした状態で電流源を推定することで、最小ノルム解の精度を維持しながら、電流源を精度よく推定する格子点移動最小ノルム法を提案している。しかし、この場合も、格子点間隔に基づく収束判定値の適切な値の設定困難等により、安定した推定結果が得られないという難点がある。
【００１４】
そこで、更に本出願人は、特願平６−４７２２０号によって、格子点上の電流源が及ぼす磁界と実際の計測磁界の２乗誤差を最小にするという条件を付加することにより、格子点間隔に基づく収束判定値を不要として、安定した推定結果を得るようにする手法を提案している。
しかし、この手法の場合も、推定された各電流源のモーメントにバラツキが発生し易く、計測磁場の離散的なノイズを解（電流源）として推定してしまう難点があるので、更に出願人は、特開平７−３２７９４３号によって、設定した各格子点上の未知の電流源が及ぼす磁界と実際の計測磁界（磁界データ）の２乗誤差（磁界２乗誤差項）と、電流源の重み係数λ付き２乗和（ペナルティ関数項）との和を最小にするという条件を付加する手法を提案している。
【００１５】
【発明が解決しようとする課題】
しかしながら、上記の特開平７−３２７９４３号に係る従来手法の場合、ペナルティ関数項における重み係数λを適切な値にセットすることが必ずしも容易でないという問題がある。従来、重み係数λは、磁界２乗誤差項とペナルティ関数項が同程度の比重となるような値にセットするか、格子間隔に応じた値にセットしているのであるが、格子点の設定位置に誤差がある場合や、計測磁場に含まれるノイズの大きさによっては、重み係数λが不適切な値となり、精度の十分な推定結果が得られない場合がある。
【００１６】
この発明は、このような事情に鑑みてなされたものであって、ペナルティ関数項における重み係数λが適切な値にセットされて、生体活動電流源の物理量について精度の十分な推定結果を得ることができる生体活動電流源推定装置を提供することを目的としている。
【００１７】
【課題を解決するための手段】
上記の課題を解決するために、請求項１の発明は、生体活動電流によって生じる磁界を複数個の磁気センサで検出することによって、生体活動電流源の位置，大きさ，方向等の物理量を推定する生体活動電流源推定装置であって、（ａ）被検体の診断対象領域に近接配備され、前記診断対象領域内の生体活動電流源による微小磁界を計測する複数個の磁気センサと、（ｂ）前記各磁気センサによって計測された磁界データをデジタルデータに変換するデータ変換手段と、（ｃ）前記デジタルデータに変換された磁界データを収集して記憶するデータ収集手段と、（ｄ）前記診断対象領域に複数個の格子点を設定する格子点設定手段と、（ｅ）前記各格子点上の未知の電流源が及ぼす磁界と前記データ収集手段に記憶された磁界データの２乗誤差と、電流源の重み係数λ付き２乗和との和を最小にするという条件を付加することにより未知の電流源を求める電流源算出手段と、（ｆ）前記電流源算出手段で推定された電流源を、前記被検体の診断対象領域の断層像に重ね合わせて表示する表示手段とを備えているとともに、（ｇ）格子点上の未知の電流源の位置と真の電流源の位置との間の位置誤差が未知の電流源の発生磁界へ与える影響の度合い（磁界影響度）を求出する磁界影響度求出手段と、（ｈ）前記磁界影響度求出手段による磁界影響度を加味して前記重み係数λを求出する重み係数求出手段とを備えており、かつ前記電流源算出手段は重み係数求出手段によって求出された重み係数λを用いて電流源を求めるよう構成されている。
【００１８】
また、請求項２の発明は、請求項１に記載の生体活動電流源推定装置において、重み係数求出手段は、重み係数λを求出する際に磁気センサによって計測された磁界データをも加味し、重み係数λが変数として含まれている最適化基準関数ＧＣＶ（λ）＝（磁界データの２乗誤差）／（磁界データの誤差ベクトルの自由度）の値を最小とする時のλの値を、求出すべき重み係数λとするよう構成されている。
【００１９】
〔作用〕
この発明の生体活動電流源推定装置の作用は次のとおりである。
請求項１の発明の装置では、各磁気センサで計測された微小磁界はデータ変換手段でデジタルデータに変換された後、データ収集手段に記憶される。そして、格子点設定手段により、診断対象領域に複数個の格子点が仮想的に設定される。普通、この格子点の個数は、磁気セサンの個数よりも少なく設定される。一方、磁界影響度求出手段により、各格子点上の未知の電流源の位置と真の電流源の位置との間の位置誤差が未知の電流源の発生磁界へ与える影響の度合い（磁界影響度）が求出されるとともに、重み係数求出手段により磁界影響度を加味した重み係数λが求出される。
続いて、電流源算出手段は、求出した重み係数λを用い、各格子点上の未知の電流源が及ぼす磁界とデータ収集手段に記憶された磁界データとの２乗誤差と、電流源の重み係数λ付き２乗和との和を最小にするという条件を付加して、各格子点上の電流源を求める。その後、推定結果としての電流源は、表示手段により被検体の診断対象領域の断層像に重ね合わせて表示される。
【００２０】
請求項１の発明の装置の場合、電流源算出過程でのペナルティ関数項における重み係数λは、従来の場合ように、磁界２乗誤差項とペナルティ関数項とが同程度の比重となるか、単に格子間隔に対応するといった単純な値ではなく、重み係数λは、未知の電流源と真の電流源との間の位置誤差に起因する磁界影響度が加味された従来より適切な値であるので、生体活動電流源についての推定精度が向上する。
【００２１】
請求項２の発明の生体活動電流源推定装置では、重み係数求出手段によって、磁界影響度だけでなく磁界データも加味され、最適化基準関数ＧＣＶ（λ）＝（磁界データの２乗誤差）／（磁界データの誤差ベクトルの自由度）の値を最小とする時のλの値が、求出すべき重み係数λであるとして決定される。
つまり、決定される重み係数λは、実測される磁界データの中のノイズの影響も考慮されることになり、しかも、最適化基準関数ＧＣＶ（λ）の値が最小であるλは、磁界データの誤差の大小だけでなく、磁界データの誤差ベクトルの自由度（誤差の方向）も考慮されている磁界データのノイズの状況に良くマッチしたものなので、重み係数λ付きのペナルティ関数項が磁界データのノイズも十分に反映した的確なものになる結果、実測された磁界データに磁界ノイズを含んでいても、十分な推定精度が得られる。
【００２２】
【実施例】
以下、図面を参照してこの発明の実施例を説明する。
図１はこの発明に係る生体活動電流源推定装置の一実施例の概略構成を示したブロック図である。図中、符号２は磁気シールドルームであり、この磁気シールドルーム２内に被検体Ｍが仰臥されるベッド３と、被検体Ｍの例えば脳に近接配備され、脳内に生じた生体活動電流源による微小磁界を無侵襲に計測するためのマルチチャンネルＳＱＵＩＤセンサ１とが設けられている。上述したように、マルチチャンネルＳＱＵＩＤセンサ１は、デュアー内に多数の磁気センサが冷媒に浸漬して収納されている。本実施例において、各磁気センサは被検体Ｍの脳を球体と仮定した場合に、その半径方向の磁界成分を検出する一対のコイルでそれぞれ構成されている。
【００２３】
マルチチャンネルＳＱＵＩＤセンサ１で検出された磁界データはデータ変換ユニット４に与えられてデジタルデータに変換された後、データ収集ユニット５に集められる。刺激装置６は、被検体Ｍに電気的刺激（あるいは音、光刺激など）を与えるためのものである。ポジショニングユニット７は、マルチチャンネルＳＱＵＩＤセンサ１を基準とした３次元座標系に対する被検体Ｍの位置関係を把握するための装置である。例えば、被検体Ｍの複数個所に小コイルを取り付け、これらの小コイルにポジショニングユニット７から給電する。そして、各コイルから発生した磁界をマルチチャンネルＳＱＵＩＤセンサ１で検出することにより、マルチチャンネルＳＱＵＩＤセンサ１に対する被検体Ｍの位置関係を把握する。なお、ＳＱＵＩＤセンサ１に対する被検体Ｍの位置関係を把握するための手法は、これ以外に、デュワーに投光器を取り付けて光ビームを被検体Ｍに照射して両者の位置関係を把握するものや、あるいは、特開平５─２３７０６５号、特開平６─７８８９２５号などに開示された種々の手法が用いられる。
【００２４】
データ解析ユニット８は、データ収集ユニット５に集められた磁界データに基づいて、被検体Ｍの診断対象領域内の電流源を推定するためのものである。このデータ解析ユニット８は、後述する説明から明らかになるように、この発明における格子点設定手段、電流源算出手段、磁界影響度求出手段、重み係数求出手段および他の各手段としての機能を備える。データ解析ユニット８に関連して設けられた光磁気ディスク９には、例えばＸ線ＣＴ装置やＭＲＩ装置で得られた断層画像が記憶されており、データ解析ユニット８で推定された電流源が、これらの断層像上に重ね合わされてカラーモニタ１０に表示されたり、あるいはカラープリンタ１１に印字出力されるようになっている。なお、Ｘ線ＣＴ装置やＭＲＩ装置で得られた断層画像は、図１に示した通信回線１２を介してデータ解析ユニット８に直接伝送するように構成してもよい。
【００２５】
上述したように、マルチチャンネルＳＱＵＩＤセンサ１を基準とした３次元座標系に対する被検体Ｍの位置関係を測定して記憶するとともに、マルチチャンネルＳＱＵＩＤセンサ１で被検体Ｍの診断対象領域である例えば脳内の生体活動電流源からの微小磁界を計測して、その磁界データをデータ収集ユニット５に集めた後、データ解析ユニット８で電流源の推定処理が、以下のようなプロセスで実行される。
【００２６】
〔格子点設定プロセス〕
図６に示した従来例と同様に、診断対象領域である例えば脳内に３次元の格子点群Ｎを想定する。ここで、格子点群Ｎについての未知数（各格子点について想定される電流源の個数）が、マルチチャンネルＳＱＵＩＤセンサ１を構成する各磁気センサＳ₁ 〜Ｓ_m の個数ｍよりも少なくなるように、格子点群Ｎの格子点の数が設定される。ここでは、格子点の数は３２個、磁気センサの数は１２９個とする。格子点の設定の仕方は特に限定されないが、臨床的知見に基づく設定アルゴリズムやシンプレックス法などの数学的手法に基づき設定アルゴリズムが挙げられる。
【００２７】
〔リードフィールド行列求出プロセス〕
そして、式(3) で表される行列Ａの各係数をビオ・サバールの法則を使って算出（行列Ａの各係数は、後述する格子点の移動ごとに算出される）する。ここでは、診断対象領域を球体と仮定し、各センサで球体の半径方向成分の磁界を検出しているので、要素数２の電流源が３２個であり、発生磁場が１２９個である。したがって、リードフィールド行列Ａは１２９行６４列となる。
行列Ａが求まると、各格子点上の電流源をペナルティ項付きの線形最小２乗法で求めることになる。
【００２８】
〔電流源算出プロセス〕
【００２９】
上述した式(4) を使って、磁気センサＳ₁ 〜Ｓ_m によって計測され、データ収集ユニット５に記憶された磁界データ〔Ｂｄ〕から各格子点上の電流源〔Ｐ〕を求める。式(4) を改めて以下に示す。
〔Ｐ〕＝Ａ^- 〔Ｂｄ〕 ………(4)
【００３０】
行列Ａは上述した(3) 式で表され、今は１２９行×６４列の個数の要素をもった行列である。(4) 式において、式の個数ｍ（磁気センサＳ₁ 〜Ｓ_m の個数）が、未知数の個数よりも多いために解が求まらない。そこで、測定された磁界〔Ｂｄ〕と、各格子点上に仮定した電流源〔Ｐ〕が磁気センサＳ₁ 〜Ｓ_m に及ぼす磁界（計算によって求められる磁界）〔Ｂ〕との２乗誤差｜〔Ｂｄ〕−〔Ｂ〕｜² に、以下に示す重み係数λ付きペナルティ項を加算した形で表された評価関数ｆを最小にするという条件を付加することにより、線形最小２乗法を用いて電流源〔Ｐ〕を算出する。この評価関数を式(7) に示す。
【００３１】
【数３】

【００３２】
評価関数の第１項が磁場の２乗誤差、第２項がペナルティ項である。ペナルティ項中のλはペナルティ項の重み係数、ｗｉは磁気センサから磁場測定面までの距離による影響をキャンセルするための係数である。ペナルティ項は、電流源が固まって存在する程、その値が小さくなるという性質があり、評価関数中の第１項の磁場の２乗誤差のみを最小にするという条件で電流源を求めた場合に、電流源がバラツキ易いという傾向を抑えるとともに、離散的に存在するノイズ成分が解として採用されることを抑える。
【００３３】
上記の評価関数ｆを最小にするという条件を付加することにより、式(4) は
〔Ｐ〕＝Ａ⁺ 〔Ｂｄ〕 ………(8)
で表される。
ここで、Ａ⁺ は次式(9) により求められる。
Ａ⁺ ＝（Ａ^t Ａ＋λ・Ｗ）^-1 ・Ａ^t ………(9)
(9) 式中のＷは、次式(10)で表されるペナルティ関数行列である。ここでは、電流源の要素数が２なのでＷは６４行６４列の行列である。(9) は正則化一般逆行列法が適用されており、行列Ａ^t Ａに行列λ・Ｗを加え正則化すると逆行列を求められる。
【００３４】
【数４】

【００３５】
ここで、一般に、電流源Ｐｉが磁場測定面に近いほど大きな磁場が測定されるので、行列Ｗを次式(11)のようにすることで、求める電流源Ｐｉについて、磁場測定面との距離による影響（すなわち、電流源Ｐｉが磁場測定面の近くに推定されるという傾向）をキャンセルすることができる。
そして、実施例の装置の場合、ペナルティ項の重み係数λは次のようにして求出する。
【００３６】
【数５】

【００３７】
〔重み係数λの求出プロセス〕
実施例の装置のペナルティ項の重み係数λの求出過程を、図２のフローチャートを参照しながら説明する。
〔ステップＳ１〕先ず格子点上の未知の電流源の位置と真の電流源の位置との間の位置誤差が未知の電流源の発生磁界へ与える影響の度合い（磁界影響度）を示す行列Ａ’を求める。行列Ａ’は１２９行６４列の行列であり、隣接する電流源のリードフィールド行列Ａの要素の差を要素としており、行列Ａ’のｉ番目の列ベクトルα’ⁱ ，行列Ａのｉ番目の列ベクトルαⁱ ，行列Ａのｊ番目の列ベクトルα^j とすると、α’ⁱ ＝｜αⁱ −α^j ｜である。ｉ＝１〜６４，ｊ＝１〜６４であるが、ｉとｊは等しくない。
【００３８】
〔ステップＳ２〕下式(12)の最適化基準関数ＧＣＶ（λ）の算出を開始する。
ＧＣＶ（λ）＝
｜ [I −（Ａ＋Ａ’ ）Ａ⁺ ] Bd｜² ÷｛trace [I−（Ａ＋Ａ’ )]Bd ｝² ……(12)
但し：(12)式中のＩは単位行列，trace は行列の対角要素の和である。
Ａ⁺ は(9) 式で示したようにλを含み、このλの値を変化させてＧＣＶ（λ）の値を次々と算出してゆく。またＧＣＶはGeneralized Cross-Validationの略である。
【００３９】
〔ステップＳ３〕ＧＣＶ（λ）の値の算出結果に基づいてＧＣＶ（λ）を最小とする時のλの値を確定する。
【００４０】
〔ステップＳ４〕ステップＳ３で確定したλの値を、求出すべき重み係数λとして決定する。
【００４１】
最適化基準関数ＧＣＶ（λ）には、行列Ａ’が含まれていることにより、未知の電流源と真の電流源との間の位置誤差（モデル化ノイズ）に起因する磁界影響度が加味され、従来より適切なλが求められる。また、最適化基準関数ＧＣＶ（λ）には、Ｂｄが含まれていて、実測される磁界データの中のノイズ（実際のノイズ）の影響も考慮されることになる。
しかも、最適化基準関数ＧＣＶ（λ）の値が最小であるλは、磁界の２乗誤差に対応する｜ [I −（Ａ＋Ａ’ ）Ａ⁺ ] Bd｜² による磁界データの誤差の大小だけでなく、磁界データの誤差ベクトル、つまり〔（Ｂｄｉ−Ｂｉ）の１２９個の要素（ｉ＝１〜１２９）を持つベクトル〕に対応する｛trace [I−（Ａ＋Ａ’ )]Bd ｝² により、磁界データの誤差ベクトルの自由度（誤差の方向）も考慮された磁界のノイズ状況に良くマッチした的確なものになる。
そして、こうして決定した重み係数λと上記した(8) 式と(9) 式により電流源〔Ｐ〕が求められる。
【００４２】
〔算出電流源の判定プロセス〕
求まった電流源〔Ｐ〕が適切か否かを次のようにチェックする。すなわち、推定された電流源〔Ｐ〕が磁気センサＳ₁ 〜Ｓ_m に及ぼす磁界〔Ｂ〕と測定された磁界〔Ｂｄ〕との２乗誤差と、上述した重み係数λ付きペナルティ項との和が、大域的に最小であるか否かを判断する。ここで、磁場の２乗誤差とペナルティ項との和が大域的に最小であるとは、後述する反復プロセスを繰り返す過程で格子点群を複数回移動させた場合に、それぞれの格子点の位置で上述した方法により求めた磁場の２乗誤差とペナルティ項との和の中で、値が最小のものをいう。磁場の２乗誤差と重み係数λ付きペナルティ項との和が大域的に最小であるか否かの判断は、後述する反復プロセスを繰り返す過程で、再配置された格子点群について、上記の過程中の必要プロセスに従って求めた磁場の２乗誤差と重み係数λ付きペナルティ項との和を記憶しておき、それぞれの値を比較して極小となる値を大域的な最小値とすればよい。
【００４３】
〔反復プロセス〕
最小でないと判断された場合は、次のようにして反復プロセスを行う。前回求めた電流源〔Ｐ〕の中、値の大きな電流源が存在する格子点の周りに他の格子点上の電流源を移動させる。これにより、格子点数は元の格子点群Ｎと同じで、間隔が密になった新たな格子点群が得られる。なお、値の大きな電流源は必ずしも一つではなく、通常は複数個ある。
【００４４】
値の大きな電流源が存在する格子点に、他の格子点群を近づけるための手法は特に限定しないが、例えば、次のような手法が例示される。前回求められた各格子点の電流源の大きさを質量と考え、重力によって各格子点間に引力が働くと仮定する。そうすると、各格子点は、質量の大きな格子点に近づいて行き、質量の大きな格子点に近い程、密度の高い格子点群が得られる。各格子点の移動距離は適宜に設定される。さらに、別の格子点移動法としては、例えば、推定電流の大きい格子点を中心とする立方体の頂点に推定電流の小さい各格子点８個を移動する。このとき、立方体の大きさは、例えば、推定電流の大きい格子点と、これに最も近い格子点との距離の半分とする。
【００４５】
格子点を移動させた後、移動させた格子点群について、再び、電流源〔Ｐ〕を新たに求めるとともに、２乗誤差と重み係数λ付きペナルティ項との和が大域的に最小であるか否かの判断を行い、最小でないと判断した場合は、また反復プロセスを繰り返し行う。
【００４６】
〔電流源の推定・表示プロセス〕
２乗誤差と重み係数λ付きペナルティ項との和が大域的に最小であると判断された場合は、判定プロセスにおいて大域的に最小となった磁界〔Ｂ〕に対する電流源〔Ｐ〕を真の電流源と推定する。真の電流源として推定された電流源〔Ｐ〕は、図１に示した光磁気ディスク９に記憶されたＸ線ＣＴ像やＭＲ断層像上に重ね合わされて、例えばカラーモニタ１０に表示される。
【００４７】
〔シミュレーション〕
実施例の生体活動電流源推定装置による電流源推定の精度向上を視覚的に確認するためのシミュレーション例を図面を参照しながら説明する。図３はシミュレーション結果を示す図であって、図３に白抜きのバツマークＱＲが示す位置に真の電流源（電流双極子）が置かれる。
第１回目のシミュレーションでは、計測磁場の信号対雑音比であるＳ／Ｎ比は５とし、反復プロセスは行わず、最初の算出結果を推定結果として、推定結果の各電流源を図３に実線矢印で示した。矢印が長いほど大きな電流源であることを示している。第１回目のシミュレーションの場合、重み係数λは９．１８２８×１０^-5となった。
【００４８】
また、Ｓ／Ｎ比を１０とした他は、第１回目と同様にして、第２回目のシュミレーションを行い、推定結果の各電流源を図４に実線矢印で示した。なお、重み係数λはλ＝９．１３５８９×１０^-5となった。
さらに、Ｓ／Ｎ比を２０とした他は、第１回目と同様にして、第３回目のシュミレーションを行い、推定結果の各電流源を図５に実線矢印で示した。なお、重み係数λは９．０８７４３×１０^-5となった。
【００４９】
図３〜図５に示された電流源の状況を観察すれば、３回のシミュレーションは全て、推定結果の各電流源のうちバツマークＱＲの示す真の電流源の近傍に位置する電流源が大きくて、推定精度が良いことが分かる。また、ＳＮ比が５〜２０と変化しても、推定結果は変わらないので、磁界ノイズの多い診断対象領域の場合にも、十分な推定精度が得られることが分かる。
【００５０】
この発明は上記の実施例に限られるものではなく、以下のように変形して実施することができる。
（１）実施例装置の場合、格子点を移動させて電流源の算出を繰り返す反復プロセスを実行できる構成であったが、反復プロセスを実行しない構成の装置が、変形例として挙げられる。診断対象領域内の真の電流源の位置が予め把握可能な場合、変形例の装置によって最初の電流源の算出で良好な推定結果が反復プロセスを実行せずに出せる。
【００５１】
（２）実施例装置の場合、重み係数求出手段は、磁界影響度と実測磁界データの両方を加味して重み係数λを求出する構成であったが、重み係数求出手段が磁界影響度だけを加味して重み係数λを求出する構成の装置が、変形例として挙げられる。
【００５２】
【発明の効果】
以上の説明から明らかなように、請求項１の発明の生体活動電流源推定装置によれば、電流源算出過程で用いられる重み係数λ付きペナルティ関数項の重み係数λが、未知の電流源と真の電流源との間の位置誤差に起因する磁界影響度が加味されて決定されるので、生体活動電流源についての推定精度が向上する。
【００５３】
また、請求項２の発明の生体活動電流源推定装置によれば、ペナルティ関数項の重み係数λは、実測磁界データ中のノイズの影響も考慮されていて、しかも、磁界データの誤差の大小だけでなく、磁界データの誤差ベクトルの自由度（誤差の方向）も考慮され、磁界データのノイズ状況に良くマッチしたものであるので、磁界データのノイズも十分に反映したより的確なものになり、計測された磁界データに磁界ノイズが含まれていても、十分な推定精度が得られる。
【図面の簡単な説明】
【図１】実施例に係る生体活動電流源推定装置の概略構成を示したブロック図である。
【図２】実施例装置における重み係数求出プロセスを示すフローチャートである。
【図３】実施例に係る第１回目の推定シミュレーション結果を示すグラフである。
【図４】実施例に係る第２回目の推定シミュレーション結果を示すグラフである。
【図５】実施例に係る第３回目の推定シミュレーション結果を示すグラフである。
【図６】従来例に係る推定処理の説明図である。
【符号の説明】
１…マルチチャンネルＳＱＵＩＤセンサ
２…磁気シールドルーム
４…データ変換ユニット
５…データ収集ユニット
８…データ解析ユニット
１０…カラーモニタ
Ｍ…被検体
Ｓ₁ 〜Ｓ_m …磁気センサ[0001]
[Industrial application fields]
The present invention relates to a life activity current source estimation device that estimates the position, orientation, and size of a life activity current source, and relates to a technique for improving the estimation accuracy of a physical quantity of a life activity current source.
[0002]
[Prior art]
When a living body is stimulated, the polarization formed across the cell membrane breaks and a biological activity current flows. This bioactive current appears in the brain and heart and is recorded as an electroencephalogram or an electrocardiogram. The magnetic field generated by the bioactive current is recorded as a magnetoencephalogram or magnetocardiogram.
[0003]
In recent years, a sensor using a SQUID (Superconducting Quantum Interface Device) has been developed as a sensor for measuring a minute magnetic field in a living body. This sensor can be placed outside the head, and a minute magnetic field generated by a current dipole (hereinafter also simply referred to as a current source) generated in the brain can be measured non-invasively with the sensor. The position, orientation, and size of the current source related to the lesion are estimated from the measured magnetic field data, and the estimated current source is displayed on a tomographic image obtained by an X-ray CT apparatus or an MRI apparatus, and the physics of the affected part, etc. It is used to specify the target position.
[0004]
Conventionally, there is a method using a minimum norm method as one of current source estimation methods (for example, WHKullmann, KDJandt, K. Rehm, HASchlitt, WJDallas and WESmith, Advances in Biomagnetism, pp.571-574, Plenum Pless , New York, 1989). Hereinafter, a conventional current source estimation method using the minimum norm method will be described with reference to FIG.
As shown in FIG. 6, the multi-channel SQUID sensor 1 is provided in the vicinity of the subject M. The multi-channel SQUID sensor 1 stores a large number of magnetic sensors (pickup coils) S _{1 to} S _m immersed in a refrigerant such as liquid nitrogen in a container called a dewar.
[0005]
On the other hand, a large number of grid points (1) to (n) are set in the diagnosis target region of the subject M, for example, the brain, and an unknown current source (current dipole) is assumed at each grid point. Is represented by a three-dimensional vehicle VP _j (j = 1 to n). Then, the magnetic fields B _{i to} B _m detected by the magnetic sensors S _{1 to} S _m of the SQUID sensor 1 are expressed by the following equation (1).
[0006]
[Expression 1]

[0007]
In equation (1), VP _j = (P _jx , P _jy , P _jz )
α _ij = (α _ijx, α _ijy, α _ijz )
It is represented by Α _ij is a known value representing the strength of the magnetic field detected at each position of the magnetic sensors S _{1 to} S _m when a current source of unit size in the X, Y, and Z directions is placed on the lattice point. It is a coefficient.
[0008]
[B] = (B ₁ , B ₂ ,..., B _m )
[P] = (P _1x , P _1y , P _1z , P _2x , P _2y , P _2z ,..., P _nx , P _ny , P _nz )
In this way, equation (1) can be rewritten as a linear relational equation like equation (2).
[B] = A [P] (2)
In the equation (2), A is a matrix (lead field matrix) having 3n × m elements expressed by the following equation (3).
[0009]
[Expression 2]

[0010]
Here, when the inverse matrix of A is represented by A ⁻ , [P] is represented by the following equation (4).
(P) = A ^- (B) ......... (4)
Here, in the minimum norm method, the number of unknowns 3n (the magnitude of the current source assumed in each lattice point in the X, Y, and Z directions) is larger than the number m (the number of magnetic sensors S _{1 to} S _m ) in the equation. The solution of the current source [P] is obtained by adding a condition that the norm | [P] | It should be noted that the solution can be obtained uniquely by making the number m of the above-mentioned equations equal to the number of unknowns 3n. However, in such a case, since the solution becomes very unstable, this minimum norm method is used. Is used.
[0011]
By adding the condition that the norm | [P] | of the current source [P] is minimized, the above equation (4) is expressed as the following equation (5).
[P] = A ⁺ [B] (5)
Here, A ⁺ is a general inverse matrix expressed by the following equation (6).
A ⁺ = A ^t (AA ^t ) ^-1 (6)
However, A ^t is the transpose matrix of A.
[0012]
Equation (5) is solved to estimate the direction and size of the current source VP _j on each lattice point, and the largest value among them is assumed to be close to the true current source. This is the principle of the current source estimation method by the minimum norm method.
[0013]
Further, in order to improve the position resolution of the minimum norm method, a method of repeatedly obtaining a minimum norm solution while subdividing the grid point division and estimating a current source closer to a true current source has been proposed. However, in this case, the number of elements of the vector [P] in equation (5) increases, and there is a problem that the calculation accuracy of the minimum norm solution is lowered.
Therefore, the applicant of the present invention previously disclosed another method in the vicinity of a lattice point where a current source having a large value exists among the current sources estimated by the minimum norm method according to JP-A-6-343613 and JP-A-6-343614. By moving the grid point group and estimating the current source with the same number of grid points as the previous time and narrowing only the grid point interval, the current source is accurate while maintaining the accuracy of the minimum norm solution. A well-estimated lattice point moving minimum norm method is proposed. However, even in this case, there is a difficulty that a stable estimation result cannot be obtained due to difficulty in setting an appropriate value for the convergence determination value based on the lattice point interval.
[0014]
Therefore, according to Japanese Patent Application No. Hei 6-47220, the present applicant adds a condition that the square error between the magnetic field exerted by the current source on the lattice points and the actual measurement magnetic field is minimized, thereby We propose a method that eliminates the need for a convergence judgment value based on, and obtains stable estimation results.
However, even in this method, the estimated moment of each current source is likely to vary, and there is a difficulty in estimating the discrete noise of the measured magnetic field as a solution (current source). According to Japanese Patent Laid-Open No. 7-327943, a square error (magnetic field square error term) between a magnetic field exerted by an unknown current source on each set grid point and an actual measurement magnetic field (magnetic field data), and a weighting factor of the current source A method of adding a condition of minimizing the sum with the sum of squares with λ (penalty function term) is proposed.
[0015]
[Problems to be solved by the invention]
However, in the case of the conventional method according to the above Japanese Patent Laid-Open No. 7-327943, there is a problem that it is not always easy to set the weighting factor λ in the penalty function term to an appropriate value. Conventionally, the weighting factor λ is set to a value such that the magnetic field square error term and the penalty function term have the same specific gravity or a value corresponding to the lattice spacing. When there is an error in the position or depending on the magnitude of noise included in the measurement magnetic field, the weighting factor λ may be an inappropriate value, and an accurate estimation result may not be obtained.
[0016]
The present invention has been made in view of such circumstances, and the weighting factor λ in the penalty function term is set to an appropriate value to obtain a sufficiently accurate estimation result for the physical quantity of the life activity current source. It is an object of the present invention to provide a bioactive current source estimation device capable of performing the above.
[0017]
[Means for Solving the Problems]
In order to solve the above problems, the invention of claim 1 estimates physical quantities such as the position, size, direction, etc. of a biological activity current source by detecting a magnetic field generated by the biological activity current with a plurality of magnetic sensors. (A) a plurality of magnetic sensors that are arranged close to a diagnosis target region of a subject and measure a micro magnetic field generated by the life active current source in the diagnosis target region; Data conversion means for converting the magnetic field data measured by each magnetic sensor into digital data, (c) data collection means for collecting and storing the magnetic field data converted into the digital data, and (d) the diagnosis Grid point setting means for setting a plurality of grid points in the target area; and (e) a square error between the magnetic field exerted by the unknown current source on each grid point and the magnetic field data stored in the data collection means. And a current source calculation means for obtaining an unknown current source by adding a condition that minimizes the sum of the sum of squares with a weight coefficient λ of the current source, and (f) estimated by the current source calculation means Display means for superimposing and displaying the current source on the tomographic image of the diagnosis target region of the subject, and (g) the position of the unknown current source on the lattice point and the position of the true current source, A magnetic field influence degree obtaining means for obtaining a degree of influence (magnetic field influence degree) of a position error between the current source and the magnetic field generated by an unknown current source; In addition, weighting factor finding means for finding the weighting factor λ is included, and the current source calculating means finds the current source using the weighting factor λ found by the weighting factor finding means. It is configured.
[0018]
The invention according to claim 2 is the life activity current source estimation device according to claim 1, wherein the weighting factor obtaining means also takes into account the magnetic field data measured by the magnetic sensor when obtaining the weighting factor λ. Then, the optimization criterion function GCV (λ) = (magnetic field data square error) / (magnetic field data error vector degree of freedom) in which the weighting coefficient λ is included as a variable is minimized. The value is configured to be a weighting factor λ to be obtained.
[0019]
[Action]
The operation of the bioactive current source estimation apparatus of the present invention is as follows.
In the apparatus according to the first aspect of the invention, the minute magnetic field measured by each magnetic sensor is converted into digital data by the data conversion means and then stored in the data collection means. Then, a plurality of grid points are virtually set in the diagnosis target area by the grid point setting means. Usually, the number of lattice points is set to be smaller than the number of magnetic sensins. On the other hand, the degree of influence of the position error between the position of the unknown current source on each grid point and the position of the true current source on the magnetic field generated by the unknown current source (field effect Degree) and a weighting factor λ taking into account the magnetic field influence degree is obtained by the weighting factor obtaining means.
Subsequently, the current source calculation means uses the obtained weighting factor λ, the square error between the magnetic field exerted by the unknown current source on each lattice point and the magnetic field data stored in the data collection means, and the current source A current source on each lattice point is obtained by adding a condition that the sum of the square sum with the weighting factor λ is minimized. Thereafter, the current source as the estimation result is displayed by being superimposed on the tomographic image of the diagnosis target region of the subject by the display means.
[0020]
In the case of the device of the invention of claim 1, the weighting factor λ in the penalty function term in the current source calculation process is, as in the conventional case, the magnetic field square error term and the penalty function term have the same specific gravity, The weighting factor λ is not a simple value that simply corresponds to the lattice spacing, but is a more appropriate value than the conventional value in consideration of the magnetic field influence due to the position error between the unknown current source and the true current source. As a result, the estimation accuracy for the bioactive current source is improved.
[0021]
In the biological activity current source estimation device according to the second aspect of the present invention, not only the magnetic field influence degree but also the magnetic field data are taken into account by the weighting factor obtaining means, and the optimization criterion function GCV (λ) = (square error of the magnetic field data). The value of λ when the value of / (degree of freedom of error vector of magnetic field data) is minimized is determined as the weighting factor λ to be obtained.
That is, the determined weighting factor λ takes into account the influence of noise in the actually measured magnetic field data, and λ having the minimum value of the optimization reference function GCV (λ) is the magnetic field data. The penalty function term with a weighting factor λ is the magnetic field data because it matches well with the noise situation of the magnetic field data taking into account not only the magnitude of the error but also the degree of freedom (error direction) of the error vector of the magnetic field data. As a result, it is possible to obtain a sufficient estimation accuracy even if magnetic field noise is included in actually measured magnetic field data.
[0022]
【Example】
Embodiments of the present invention will be described below with reference to the drawings.
FIG. 1 is a block diagram showing a schematic configuration of an embodiment of a life activity current source estimation apparatus according to the present invention. In the figure, reference numeral 2 denotes a magnetic shield room, a bed 3 on which the subject M is supine in the magnetic shield room 2, and a biologically active current source generated in the brain in the proximity of the subject M, for example, the brain. And a multi-channel SQUID sensor 1 for non-invasively measuring a minute magnetic field due to. As described above, the multi-channel SQUID sensor 1 has a large number of magnetic sensors immersed in a refrigerant and housed in a dewar. In this embodiment, each magnetic sensor is composed of a pair of coils that detect a magnetic field component in the radial direction when the brain of the subject M is assumed to be a sphere.
[0023]
Magnetic field data detected by the multichannel SQUID sensor 1 is applied to the data conversion unit 4 and converted into digital data, and then collected in the data collection unit 5. The stimulation device 6 is for applying electrical stimulation (or sound, light stimulation, etc.) to the subject M. The positioning unit 7 is a device for grasping the positional relationship of the subject M with respect to the three-dimensional coordinate system with the multi-channel SQUID sensor 1 as a reference. For example, small coils are attached to a plurality of locations of the subject M, and power is supplied from the positioning unit 7 to these small coils. Then, the magnetic field generated from each coil is detected by the multichannel SQUID sensor 1, thereby grasping the positional relationship of the subject M with respect to the multichannel SQUID sensor 1. In addition to this, a method for grasping the positional relationship of the subject M with respect to the SQUID sensor 1 is to attach a projector to the Dewar and irradiate the subject M with a light beam to grasp the positional relationship between the two, Alternatively, various techniques disclosed in JP-A-5-237065, JP-A-6-788925, and the like are used.
[0024]
The data analysis unit 8 is for estimating a current source in the diagnosis target region of the subject M based on the magnetic field data collected by the data collection unit 5. The data analysis unit 8 functions as a grid point setting unit, a current source calculation unit, a magnetic field influence degree calculation unit, a weighting factor calculation unit, and other units as will be apparent from the following description. Is provided. A magneto-optical disk 9 provided in association with the data analysis unit 8 stores, for example, a tomographic image obtained by an X-ray CT apparatus or an MRI apparatus, and a current source estimated by the data analysis unit 8 is These tomographic images are superimposed on each other and displayed on the color monitor 10 or printed out on the color printer 11. Note that the tomographic image obtained by the X-ray CT apparatus or the MRI apparatus may be directly transmitted to the data analysis unit 8 via the communication line 12 shown in FIG.
[0025]
As described above, the positional relationship of the subject M with respect to the three-dimensional coordinate system based on the multichannel SQUID sensor 1 is measured and stored, and the diagnosis target region of the subject M is detected by the multichannel SQUID sensor 1, for example, the brain After measuring a minute magnetic field from the internal biological activity current source and collecting the magnetic field data in the data collection unit 5, the data analysis unit 8 performs the current source estimation process in the following process.
[0026]
[Lattice point setting process]
Similar to the conventional example shown in FIG. 6, a three-dimensional lattice point group N is assumed in the diagnosis target region, for example, the brain. Here, the unknown (the number of current sources assumed for each grid point) for the grid point group N is smaller than the number m of the magnetic sensors S _{1 to} S _m constituting the multichannel SQUID sensor 1. The number of grid points of the grid point group N is set. Here, the number of lattice points is 32, and the number of magnetic sensors is 129. The method of setting the lattice points is not particularly limited, but examples include setting algorithms based on clinical knowledge and mathematical algorithms such as the simplex method.
[0027]
[Lead field matrix search process]
Then, each coefficient of the matrix A expressed by the equation (3) is calculated using Bio-Savart's law (each coefficient of the matrix A is calculated every time a lattice point is moved later). Here, the diagnosis target region is assumed to be a sphere, and the magnetic field of the radial component of the sphere is detected by each sensor. Therefore, there are 32 current sources with 2 elements and 129 generated magnetic fields. Therefore, the lead field matrix A is 129 rows and 64 columns.
When the matrix A is obtained, the current source on each lattice point is obtained by the linear least square method with a penalty term.
[0028]
[Current source calculation process]
[0029]
The current source [P] on each lattice point is obtained from the magnetic field data [Bd] measured by the magnetic sensors S _{1 to} S _m and stored in the data acquisition unit 5 using the above-described equation (4). Equation (4) is shown below again.
(P) = A ^- [Bd] ......... (4)
[0030]
The matrix A is expressed by the above-described equation (3) and is now a matrix having 129 rows × 64 columns of elements. In equation (4), since the number m of equations (the number of magnetic sensors S _{1 to} S _m ) is greater than the number of unknowns, no solution is obtained. Therefore, the square error between the measured magnetic field [Bd] and the magnetic field (magnetic field obtained by calculation) [B] exerted on the magnetic sensors S _{1 to} S _m by the assumed current source [P] on each lattice point | By adding a condition that minimizes the evaluation function f expressed by adding a penalty term with a weighting coefficient λ shown below to [Bd] − [B] | ² , the linear least square method is used. The current source [P] is calculated. This evaluation function is shown in Equation (7).
[0031]
[Equation 3]

[0032]
The first term of the evaluation function is the square error of the magnetic field, and the second term is the penalty term. Λ in the penalty term is a weighting factor of the penalty term, and wi is a coefficient for canceling the influence due to the distance from the magnetic sensor to the magnetic field measurement surface. The penalty term has a property that the value becomes smaller as the current source is solidified, and the current source is obtained under the condition that only the square error of the magnetic field of the first term in the evaluation function is minimized. In addition, the tendency that the current sources are likely to vary is suppressed, and the noise components that exist discretely are suppressed from being adopted as a solution.
[0033]
By adding the condition of minimizing the above evaluation function f, the equation (4) becomes [P] = A ⁺ [Bd] (8)
It is represented by
Here, A ⁺ is obtained by the following equation (9).
A ⁺ = (A ^t A + λ · W) ⁻¹ · A ^t (9)
W in equation (9) is a penalty function matrix expressed by the following equation (10). Here, since the number of elements of the current source is 2, W is a 64 × 64 matrix. In (9), the regularized general inverse matrix method is applied. When the matrix λ · W is added to the matrix A ^t A and regularized, the inverse matrix can be obtained.
[0034]
[Expression 4]

[0035]
Here, in general, the closer the current source Pi is to the magnetic field measurement surface, the larger the magnetic field is measured. Therefore, by setting the matrix W to the following equation (11), the distance from the magnetic field measurement surface to the current source Pi to be obtained is determined. (That is, the tendency that the current source Pi is estimated near the magnetic field measurement surface) can be canceled.
In the case of the apparatus according to the embodiment, the weighting factor λ of the penalty term is obtained as follows.
[0036]
[Equation 5]

[0037]
[Calculation process of weighting factor λ]
A process for obtaining the weighting factor λ of the penalty term in the apparatus of the embodiment will be described with reference to the flowchart of FIG.
[Step S1] First, a matrix A indicating the degree of influence (magnetic field influence degree) that the position error between the position of the unknown current source on the lattice point and the position of the true current source has on the magnetic field generated by the unknown current source. Ask for '. The matrix A ′ is a matrix of 129 rows and 64 columns, and the difference between the elements of the lead field matrix A of the adjacent current source is an element, and the i-th column vector α ′ ^{i of} the matrix A ′ and the i-th column of the matrix A Assuming that the column vector α ^{i is} the j-th column vector α ^j of the matrix A, α ′ ⁱ = | α ⁱ −α ^j |. i = 1 to 64 and j = 1 to 64, but i and j are not equal.
[0038]
[Step S2] Calculation of the optimization reference function GCV (λ) of the following equation (12) is started.
GCV (λ) =
｜ [I − (A + A ′ ) A ⁺ ] Bd | ² ÷ {trace [I− (A + A ′)] Bd} ² …… (12)
However, in equation (12), I is the unit matrix, and trace is the sum of the diagonal elements of the matrix.
A ⁺ includes λ as shown in equation (9), and the value of λ is changed to calculate the value of GCV (λ) one after another. GCV is an abbreviation for Generalized Cross-Validation.
[0039]
[Step S3] The value of λ when GCV (λ) is minimized is determined based on the calculation result of the value of GCV (λ).
[0040]
[Step S4] The value of λ determined in Step S3 is determined as the weighting factor λ to be obtained.
[0041]
Since the optimization reference function GCV (λ) includes the matrix A ′, the influence of the magnetic field due to the position error (modeling noise) between the unknown current source and the true current source is taken into consideration. Thus, a more appropriate λ is required than before. Further, the optimization reference function GCV (λ) includes Bd, and the influence of noise (actual noise) in the actually measured magnetic field data is also taken into consideration.
Moreover, λ having the minimum value of the optimization criterion function GCV (λ) corresponds to the square error of the magnetic field | [I− (A + A ′ ) A ⁺ ] Bd | Not only corresponds to the magnitude of the error of the magnetic field data due to ² , but also corresponds to the error vector of the magnetic field data, that is, a vector having 129 elements (i = 1 to 129) of (Bdi−Bi). By {trace [I− (A + A ′)] Bd} ² , the accuracy becomes well matched to the noise situation of the magnetic field in consideration of the degree of freedom (error direction) of the error vector of the magnetic field data.
Then, the current source [P] is obtained by the weighting factor λ determined in this way and the above equations (8) and (9).
[0042]
[Calculation process of calculated current source]
Whether the obtained current source [P] is appropriate is checked as follows. That is, the sum of the square error between the magnetic field [B] and the measured magnetic field [Bd] exerted on the magnetic sensors S _{1 to} S _m by the estimated current source [P] and the penalty term with the weighting factor λ described above. Is determined to be globally minimal. Here, the sum of the square error of the magnetic field and the penalty term is the smallest globally when the lattice point group is moved a plurality of times in the process of repeating the iterative process described later. In the sum of the square error of the magnetic field and the penalty term obtained by the method described above, the one having the smallest value. Whether or not the sum of the square error of the magnetic field and the penalty term with the weighting factor λ is the global minimum is determined by repeating the above-described process for the relocated grid point group in the process of repeating the iterative process described later. The sum of the square error of the magnetic field obtained according to the necessary process and the penalty term with the weighting factor λ may be stored, and the minimum value obtained by comparing the respective values may be set as the global minimum value.
[0043]
[Iterative process]
If it is determined that it is not the minimum, the iterative process is performed as follows. The current source on another lattice point is moved around the lattice point where the current source having a large value exists in the current source [P] obtained last time. As a result, a new lattice point group having the same number of lattice points as the original lattice point group N and a close interval is obtained. Note that the current source having a large value is not necessarily one, and there are usually a plurality of current sources.
[0044]
A method for bringing another lattice point group closer to a lattice point where a current source having a large value exists is not particularly limited. For example, the following method is exemplified. The current source size of each lattice point obtained last time is considered as mass, and it is assumed that gravitational force acts between each lattice point by gravity. Then, each lattice point approaches a lattice point having a large mass, and a lattice point group having a higher density is obtained as the lattice point becomes closer to a large mass. The moving distance of each lattice point is set appropriately. Further, as another lattice point moving method, for example, 8 lattice points each having a small estimated current are moved to the vertices of a cube centering on the lattice point having a large estimated current. At this time, the size of the cube is, for example, half of the distance between the lattice point having a large estimated current and the lattice point closest thereto.
[0045]
After moving the grid point, a new current source [P] is again obtained for the moved grid point group, and whether the sum of the square error and the penalty term with the weighting factor λ is globally minimum. If it is determined that it is not minimum, the iterative process is repeated.
[0046]
[Current source estimation and display process]
If it is determined that the sum of the square error and the penalty term with the weighting factor λ is globally minimal, the current source [P] for the magnetic field [B] that is globally minimized in the determination process is Estimated current source. The current source [P] estimated as a true current source is superimposed on the X-ray CT image or MR tomographic image stored in the magneto-optical disk 9 shown in FIG. .
[0047]
〔simulation〕
A simulation example for visually confirming improvement in accuracy of current source estimation by the life activity current source estimation apparatus of the embodiment will be described with reference to the drawings. FIG. 3 is a diagram showing a simulation result, and a true current source (current dipole) is placed at the position indicated by the white cross mark QR in FIG.
In the first simulation, the S / N ratio, which is the signal-to-noise ratio of the measured magnetic field, is set to 5, the iterative process is not performed, the first calculation result is used as the estimation result, and the current sources of the estimation result are shown as solid lines in FIG. Shown with an arrow. A longer arrow indicates a larger current source. In the case of the first simulation, the weighting factor λ becomes 9.1828 × 10 ⁻⁵ .
[0048]
Further, the second simulation was performed in the same manner as the first time except that the S / N ratio was set to 10, and each current source of the estimation result is indicated by a solid line arrow in FIG. The weighting factor λ is λ = 9.135989 × 10 ⁻⁵ .
Further, the third simulation was performed in the same manner as the first time except that the S / N ratio was set to 20, and each current source of the estimation result is indicated by a solid arrow in FIG. The weighting factor λ was 9.08743 × 10 ⁻⁵ .
[0049]
Observing the state of the current source shown in FIGS. 3 to 5, all three simulations have a large current source located in the vicinity of the true current source indicated by the cross mark QR among the current sources of the estimation results. Thus, it can be seen that the estimation accuracy is good. Further, even if the SN ratio changes from 5 to 20, the estimation result does not change, so that it can be understood that sufficient estimation accuracy can be obtained even in the diagnosis target region with a lot of magnetic field noise.
[0050]
The present invention is not limited to the above-described embodiments, and can be carried out with the following modifications.
(1) In the case of the embodiment apparatus, the configuration is such that an iterative process in which the calculation of the current source is repeated by moving the grid point, but an apparatus having a configuration in which the iterative process is not performed can be given as a modification. When the position of the true current source in the diagnosis target region can be grasped in advance, a favorable estimation result can be obtained without performing an iterative process by calculating the first current source by the modified apparatus.
[0051]
(2) In the case of the embodiment apparatus, the weighting factor obtaining means is configured to obtain the weighting factor λ in consideration of both the magnetic field influence degree and the measured magnetic field data. An apparatus having a configuration in which the weighting factor λ is obtained with only the degree taken into account is given as a modification.
[0052]
【The invention's effect】
As is clear from the above description, according to the life activity current source estimation device of the invention of claim 1, the weighting factor λ of the penalty function term with the weighting factor λ used in the current source calculation process is the same as the unknown current source. Since the magnetic field influence degree resulting from the position error with the true current source is taken into consideration, the estimation accuracy of the life activity current source is improved.
[0053]
Further, according to the life activity current source estimation device of the invention of claim 2, the weighting factor λ of the penalty function term takes into consideration the influence of noise in the measured magnetic field data, and only the magnitude of the error in the magnetic field data is considered. In addition, the degree of freedom of the error vector of the magnetic field data (the direction of the error) is also taken into account and matches well with the noise situation of the magnetic field data, so it becomes more accurate that sufficiently reflects the noise of the magnetic field data, Even if magnetic field noise is included in the measured magnetic field data, sufficient estimation accuracy can be obtained.
[Brief description of the drawings]
FIG. 1 is a block diagram illustrating a schematic configuration of a life activity current source estimation apparatus according to an embodiment.
FIG. 2 is a flowchart illustrating a weighting factor obtaining process in the embodiment apparatus.
FIG. 3 is a graph showing a first estimation simulation result according to the example.
FIG. 4 is a graph showing a second estimation simulation result according to the example.
FIG. 5 is a graph showing a third estimation simulation result according to the example.
FIG. 6 is an explanatory diagram of an estimation process according to a conventional example.
[Explanation of symbols]
1 ... multichannel SQUID sensor 2 ... magnetic shield room 4 ... data conversion unit 5 ... data collection unit 8 ... data analysis unit 10 ... color monitor M ... subject S ₁ to S _m ... magnetic sensor

Claims (2)

A biological activity current source estimation device that estimates physical quantities such as position, size, direction, and the like of a biological activity current source by detecting a magnetic field generated by the biological activity current with a plurality of magnetic sensors, and (a) a subject A plurality of magnetic sensors that are disposed in proximity to the diagnosis target area and measure a minute magnetic field generated by a bioactive current source in the diagnosis target area; and (b) convert the magnetic field data measured by each magnetic sensor into digital data. Data conversion means, (c) data collection means for collecting and storing magnetic field data converted into the digital data, and (d) grid point setting means for setting a plurality of grid points in the diagnosis target region (E) Minimize the sum of the magnetic field exerted by the unknown current source on each grid point, the square error of the magnetic field data stored in the data collection means, and the sum of squares with the weighting factor λ of the current source And (f) a current source estimated by the current source calculation unit is superimposed on a tomographic image of the diagnosis target region of the subject. (G) degree of influence of the position error between the position of the unknown current source on the lattice point and the position of the true current source on the generated magnetic field of the unknown current source. (H) a weighting factor obtaining unit for obtaining the weighting factor λ in consideration of the magnetic field influence by the magnetic field influence degree finding unit. A life activity current source estimation device comprising: the current source calculation unit configured to obtain a current source using the weighting factor λ obtained by the weighting factor obtaining unit. 2. The life activity current source estimation device according to claim 1, wherein the weighting factor obtaining means takes into account the magnetic field data measured by the magnetic sensor when obtaining the weighting factor λ, and includes the weighting factor λ as a variable. Weighting factor λ to be obtained as the value of λ when the value of the optimization criterion function GCV (λ) = (square error of magnetic field data) / (degree of freedom of error vector of magnetic field data) is minimized A biologically active current source estimation device configured as described above.

Images