JP4617444B2 - 脳内電流源推定方法、脳内電流源推定プログラムおよび脳内電流源推定装置 - Google Patents

Description

本発明は、脳内活動に対応して、頭蓋外部に磁場または電場を生成する脳内電流源の位置や分布を推定する方法、脳内電流源推定プログラムおよび脳内電流源推定装置に関する。

近年の生体計測技術の進歩はめざましく、従来測定が困難で誤差も大きかった脳から発生する微弱電場（脳波）や微弱磁場（脳磁波）の計測精度は、年々向上している。

すなわち、外部からの刺激を受けて、脳内の神経細胞は電流を発生する。この電流により、上述したような微弱な電場や微弱な磁場が発生することになる。このうち「脳波」とは、この神経細胞の電流により、脳から発生する電場のことである。「脳波計（Electroencephalogram:ＥＥＧ）」は、このような脳波を測定する手法である。

一方、「脳磁波」とは、この神経細胞の電流により、脳から発生する磁場のことである。「脳磁計（Magnetoencephalography:ＭＥＧ）」は、このような脳磁波を測定する手法である。脳磁波計測の最大の利点は、磁場が容積導体による影響をほとんど受けないので、頭蓋外からの磁気計測により、脳内電流源位置を３次元的に比較的正確に推定できることが期待されることである。

このような脳磁波の解析においては、発生した磁場を脳の外部から測定することで非侵襲的に脳の活動部位を推定する。

しかしながら、この磁場は非常に微弱であるために、地磁気等の外部の磁場の影響を大きく受けてしまう。そこで、外部の磁場を遮蔽するシールドを作り、その中で超伝導を用いた測定機器である超伝導量子干渉素子（SQUID：Superconducting quantum interference device）により微弱な磁場を測定することが行われる。

このような「脳内電流源位置推定」のアルゴリズムの従来の研究としては、たとえば、「双極子推定法」がある（たとえば、非特許文献１を参照）。しかしながら、この「双極子推定法」は、脳内の電流源が１つか数個の電流双極子で表せると仮定して観測磁場から双極子の位置を推定する方法であり、双極子の数を幾つにすれば良いかを決めるのが難しいという欠点を持っている。

一方、他のアルゴリズムとしては、「空間フィルター法」がある（たとえば、非特許文献２を参照）。「空間フィルター法」は、生理学的な知見から、脳内電流源の場所を限定し、双極子の分布を推定する方法である。しかしながら、この方法の欠点は電流源の深さを正確に推定できない事である。

そこで、本願の発明者により、以上のような問題点を解決するための「脳内電流源位置推定」のアルゴリズムとして、「階層変分ベイズ推定法」が提案されている（たとえば、特許文献１を参照）。

以下、この特許文献１に記載されている「階層変分ベイズ推定法」の基本原理の概念を背景として説明する。

［脳内電流源推定の原理］
（曲面上の電流分布による電磁場の復元）
電流源を金属と磁性体とで作られた理想的な電磁シールド面で囲み込めば、シールド面の外側に電磁場が存在しなくなることは電磁シールドの原理として良く知られている。すなわち、シールド面上を流れる電流と磁性体を構成する微小磁石の集まりが作る電磁場により、電流源がシールド面の外側に作る電磁場が完全に打ち消されることになる。また、微小磁石は仮想的な渦電流で置き換えることができる。これは見方を変えれば、このシールド面に適当な電流を流してやれば電流源が作る電磁場と全く同じ電磁場をシールド面の外側に作れることを示している。逆にシールド面の外側にある電流源が作る電磁場はこのシールド面にどのような電流を流しても完全には復元できない。この事実を以後「電磁場復元の原理」と呼ぶ。

図１４は、電流源が作る電磁場を適当な曲面上で観測している場合の概念図である。

図１４に示すように、観測曲面と電流源の間に仮想曲面を用意すると、上に述べた「電磁場復元の原理」から仮想曲面上に適当な電流を流してやれば電流源が観測面上に作る電磁場を復元できることが分かる。

図１５は、脳内の電流源と複数の仮想曲面との関係を示す概念図である。

図１５に示すとおり、電流が作る電磁場は距離の二乗に逆比例して減衰することを考慮すれば、電流源と等価な仮想曲面（仮想曲面２）上の電流の広がりは仮想曲面が電流源から遠ざかるほど広くなってゆく。したがって、仮想曲面１上の電流の広がりは、仮想曲面２上の電流の広がりよりも大きい。

一方、電流源に対して観測曲面と反対側にある仮想曲面（仮想曲面３）では電流源が観測面に作る電磁場を完全に復元する事が出来ない。

「階層変分ベイズ推定法」では、以上の原理を元に観測面上での電磁場の観測データから電流源を推定する。

（電流源推定の原理）
図１４に示すように、脳内で発生した電流源から作られる磁場（または電場〉を脳の表面近くの観測面で観測していると仮定する。

脳内に仮想曲面を考え、観測された磁場を復元する仮想曲面上の電流分布を求める。この仮想曲面を脳表面から半径を縮めながら脳の内部へ移動させてゆくと、観測磁場を復元する電流分布は、その広がりが小さくなってゆき真の電流源を仮想曲面が含むときに最小になる。また仮想曲面を電流源より深く移動させていくと仮想曲面上での電流分布の広がりは再び広がってゆき、この電流が生ずる磁場と観測磁場との誤差も大きくなってゆく。

すなわち、観測磁場を復元する電流分布の広がりと磁場の復元誤差を見ることにより電流源の深さが特定できる。またこのようにして特定された深さの仮想曲面上の電流分布を求めれば電流源の広がりもわかることになる。以上が本発明の原理である。

（階層変分ベイズ推定法による電流源の深さの特定）
上述した電流源推定の原理を具体的に実施するにあたり、「階層変分ベイズ推定法」では、以下のような確率論に基づく手続きをとる。

ＭＥＧやＥＥＧで観測できるのは、脳表面近くに設置された複数のセンサにより測定される数十カ所から数百カ所の磁場（電場）である。また仮想曲面上の電流分布を近似するために、仮想曲面上に格子点を設定し各格子点に電流双極子（または適当な電流源モデル）を割り当てる。電流分布を高い分解能で推定するためには格子点の数を増やして格子点の密度を大きくする必要がある。

一方、仮想曲面上の格子点の数を観測点の数より増やすと解が一意に決まらない。また観測点より推定点の数の方が多いと、観測データ数よりも推定すべきパラメータの数の方が多くなるので、電流源より深い位置にある仮想曲面でも観測磁場の復元性が良くなる。

このような困難を解決するために、べイズ推定理論を利用して仮想曲面上の電流分布を推定する。ベイズ推定を行う際に、電流源に関する事前知識を表す事前分布を用いる。

具体的には、脳内電流源は複数の場所に局在していると考えられるので、電流分布の広がりができるだけ小さく推定されるような事前分布、すなわち、観測データだけからではその強さが良く特定できないような格子点上の電流双極子は、その強さがゼロに近くなるような事前分布を導入する。これは有効信号路自動決定（Automatic Relevance Determination）事前分布（文献：R.M. Neal, Bayesian Learning for Neural Networks, Springer Verlag., 1996）と呼ばれる階層的事前分布を導入することで実現できる。この事前分布を以後、「局在条件事前分布」と呼ぶことにする。なお、上記特許文献１にも記載されている通り、「階層変分ベイズ推定法」では、このような局在条件以外の事前知識を導入した事前分布を選んで推定を行うことも可能である。

しかし、局在条件事前分布と観測データから事後確率分布を解析的に計算することは不可能である。そこで、「階層変分ベイズ推定法」では、後に説明するように、事後確率分布を近似的に計算する手法として変分法的ベイズ推定法（Variational Bayes method）（文献：Attias, H. Proc. 15th Conference on Uncertainty in Artificial Intelligence pp.21-30, 1999およびMasa-aki Sato, Neural Computation, 13, 1649-1681, 2001）を用いる。

この局在条件事前分布を用いたベイズ推定法を行うことにより、観測データを復元し、かつ出来るだけ広がりの少ない仮想曲面上の電流分布が得られる。また深さの異なる仮想曲面を使って推定を行ったときのモデル事後確率を比べることにより電流源の深さが特定できる。

モデル事後確率の対数は対数尤度項とモデル複雑度項の和で表される。対数尤度項は復元誤差が小さいほど大きくなる。

一方、モデル複雑度項は有効な（すなわち適当なしきい値以上の強さを持つ）格子点上の電流双極子の数が少ないほど大きくなる。上述したように、復元誤差と仮想曲面上の電流分布の広がり（有効な電流双極子の数）は仮想曲面の深さが電流源と一致したときに最小になるので、この時モデル事後確率が最大になることが分かる。すなわちモデル事後確率が最大になる深さの仮想曲面上に電流源が存在することが分かる。

すなわち、以上をまとめると、「階層変分ベイズ推定法」は、脳磁波（ＭＥＧ）や脳波(ＥＥＧ）の観測データから脳内電流源の位置を深さ方向まで含めて推定する方法を提供する。

（複数個の電流源の存在する場合）
以上では、電流源が一つの場合について説明したが、この方法は電流源が複数ある場合にも適用できる。

電流により生ずる電磁場は距離の二乗に逆比例して減衰するので、脳表面に最も近い電流源が脳表面の観測磁場に最も影響を与える。そこで脳表面に近い電流源から順番に特定することが出来る。

仮想曲面を脳表面からだんだん深くしてゆくと脳表面に最も近い電流源（これを第１電流源と呼ぶ）の近くでモデル事後確率が極大になる。電流源が２つ以上ある場合にはこの仮想曲面上の電流分布には電流源の数に対応した局所的な電流双極子の集合が複数個出来る。

この局所的な電流双極子の集合を局在電流分布と呼ぶ。個々の局在電流分布を含む局所面を切り出し、これを深さ方向に動かしてモデル事後確率を求める。第１電流源に対応した局所面を動かしたときにはモデル事後確率が第１電流源の深さで極大になる。しかし、これ以外の局所面を動かしてもモデル事後確率は第１電流源の深さで極大になることはない。このことから第１電流源の位置、すなわち最も表面に近い電流源の位置が特定できる。

２番目に深い電流源を求めるために、第１電流源に対応する局所面を固定し、残りの局所面の深さを同じにして深くしてゆく。すると今度は２番目に深い電流源のところでモデル事後確率が極大になる。ここでまた個々の局所面を深さ方向に動かすと第２電流源に対応する局所面を動かしたときにのみ第２電流源の深さでモデル事後確率が極大になる。こうして第２電流源の位置を特定できる。３番目以降の電流源も同様にして順次特定することが出来る。
Mosher.J.C．，Lewis P.S．and Leahy R.M. IEEE Trans．Biomed．Engng．〈1992〉vol．39，pp.541−557 Toyama K, Yoshizawa K, Yoshida Y, Kondo Y, Tomita S, Takanashi Y, Ejima Y, and Yoshikawa S. Neuroscience. 1999;91(2):pp.405-415 国際公開第ＷＯ０３／０５７０３５Ａ１号パンフレット

しかしながら、現在のＭＥＧ装置等では、一般にセンサの数が数百個程度である。一方、ミリメートル単位の高空問分解能で脳内電流分布を逆推定するには、数千から数万個程度の未知パラメータを推定しなければならない。このため、逆推定問題は高度の不良設定問題になっている。

本発明の目的は、高空問分解能で、観測データから脳内電流源の位置を深さ方向まで含めて推定することが可能な脳内電流源推定方法、脳内電流源推定プログラムおよび脳内電流源推定装置を提供することである。

本発明では、そこでセンサーを計測ごとにアクティブに動かすことで、実質的なセンサー数を数万個のレベルにまで引き上げ高空間分解能での推定精度を大幅に向上させる。

係る目的を達成するために本願発明に係る脳内電流源推定方法は、脳内に存在する電磁場の源となる電流源の位置を推定するための方法であって、観測対象の頭蓋外部において、観測曲面上に配置された複数の観測センサを、観測曲面に沿って指定距離だけ移動させ複数のセンサ位置において各々複数回にわたって電磁場の時系列データを測定して準備するステップと、脳内に脳表面からの深さが互いに異なり、かつ互いに交わらない形状を有する複数の仮想曲面と、各仮想曲面上に格子点を設定するステップと、準備された電磁場の時系列データについて時間平均した平均電磁場が観測される確率を脳内の電流分布により表現したときに、各仮想曲面上において、それぞれ平均電磁場を復元するための電流分布を階層変分ベイズ推定により推定するステップと、複数の仮想曲面上において推定された電流分布の広がりと、電流分布に基づいて復元される電磁場と観測された電磁場との誤差とに基づいて、複数の仮想曲面のうちから広がりと誤差とが極小となる仮想曲面を、真の電流源の存在する曲面として特定するステップとを備える。

好ましくは、電流分布を推定するステップは、事前分布と電磁場の観測データとから、ベイズ推定法により事後確率を求めるステップを含み、真の電流源の存在する曲面として特定するステップは、仮想曲面のうち、対応する事後確率が最大となる仮想曲面を特定するステップを含む。

この発明の他の局面にしたがうと、脳内に存在する電磁場の源となる電流源の位置を推定するコンピュータのためのプログラムであって、観測対象の頭蓋外部において、観測曲面上に配置された複数の観測センサを、観測曲面に沿って指定距離だけ移動させ複数のセンサ位置において各々複数回にわたって電磁場の時系列データを測定して記憶装置に格納するステップと、脳内に脳表面からの深さが互いに異なり、かつ互いに交わらない形状を有する複数の仮想曲面と、各仮想曲面上に格子点を設定するステップと、準備された電磁場の時系列データについて時間平均した平均電磁場が観測される確率を脳内の電流分布により表現したときに、各仮想曲面上において、それぞれ平均電磁場を復元するための電流分布を階層変分ベイズ推定により推定するステップと、複数の仮想曲面上において推定された電流分布の広がりと、電流分布に基づいて復元される電磁場と観測された電磁場との誤差とに基づいて、複数の仮想曲面のうちから広がりと誤差とが極小となる仮想曲面を、真の電流源の存在する曲面として特定するステップとをコンピュータに実行させる。

好ましくは、電流分布を推定するステップは、事前分布と電磁場の観測データとから、ベイズ推定法により事後確率を求めるステップを含み、真の電流源の存在する曲面として特定するステップは、仮想曲面のうち、対応する事後確率が最大となる仮想曲面を特定するステップを含む。

この発明のさらに他の局面に従うと、脳内に存在する電磁場の源となる電流源の位置を推定するための脳内電流源推定装置であって、電磁場の観測データを格納するための記憶装置と、観測曲面上に配置され、電磁場を測定するための複数の観測センサと、観測対象の頭蓋外部において、複数の観測センサを、観測曲面に沿って指定距離だけ移動させるセンサ移動手段と、複数の指定距離だけ複数のセンサを移動させ、複数のセンサ位置において各々複数回にわたって電磁場の時系列データを測定して記憶装置に格納する測定制御手段と、脳内に脳表面からの深さが互いに異なり、かつ互いに交わらない形状を有する複数の仮想曲面と、各仮想曲面上に格子点を設定する仮想曲面設定手段と、準備された電磁場の時系列データについて時間平均した平均電磁場が観測される確率を脳内の電流分布により表現したときに、各仮想曲面上において、それぞれ平均電磁場を復元するための電流分布を階層変分ベイズ推定により推定する電流分布推定手段と、複数の仮想曲面上において推定された電流分布の広がりと、電流分布に基づいて復元される電磁場と観測された電磁場との誤差とに基づいて、複数の仮想曲面のうちから広がりと誤差とが極小となる仮想曲面を、真の電流源の存在する曲面として特定する電流源特定手段とを備える。

好ましくは、電流分布推定手段は、事前分布と電磁場の観測データとから、ベイズ推定法により事後確率を求める事後確率導出手段を含み、電流源特定手段は、仮想曲面のうち、対応する事後確率が最大となる仮想曲面を特定する仮想曲面特定手段を含む。

本発明の脳内電流源推定方法に従えば、ＭＥＧやＥＥＧの観測データから脳内電流源の位置を深さ方向まで含めて、高空問分解能で推定することが可能である。しかも、このような深さ方向の推定は、複数の電流源が存在する場合にも適用可能である。また、初期分解能での推定の後に、位置の推定分解能を逐次的に上げる処理を追加して行うことが可能である。

以下、図面を参照して本発明の実施の形態について説明する。

なお、上述のとおり、脳内活動を外部から観測する手法としては、脳磁計（Magnetoencephalography:ＭＥＧ）や脳波計（Electroencephalogram:ＥＥＧ）がある。以下では、これらのうち、主にＭＥＧについて説明することにする、ただし、ＭＥＧにおける磁場を電場に置き換えればＥＥＧについても、本発明を適用することが可能である。

したがって、本明細書においては、本来「電場」と「磁場」とが複合して存在する「電磁場」に対して、「観測対象」となる物理量が「電場」である場合、または「磁場」である場合とを総称して、「電磁場を観測する」と呼ぶことにする。

［実施の形態］
（本発明のシステム構成）
図１は、本発明の脳磁計システム１０００の構成の一例を示す概念図である。

脳磁計１２は、マルチチャンネルのＳＱＵＩＤ磁束計（超高感度磁気測定器）アレイ２０を備えており、被験者１０の頭蓋表面での磁場を測定する。

コンピュターシステム１００は、この脳磁計１２の動作を制御するとともに、この脳磁計１２からの測定結果を受け取って、測定結果を解析し、脳内の電流源の位置の推定処理を行う。

なお、後に説明するように、ＳＱＵＩＤ磁束計アレイ（以下、「センサアレイ」と呼ぶ）２０には、数百個のＳＱＵＩＤセンサ２０．ｉ（ｉ：自然数）が、被験者１０の脳に近接する観測曲面（たとえば、半球）に沿って配置されており、さらに、このセンサアレイ２０は、コンピュータシステム１００の制御により、この観測曲面に沿って指定された距離だけ移動可能な機構を有するものとする。

図２は、コンピュータシステム１００の構成をブロック図形式で示す図である。

図２を参照してこのコンピュータシステム１００は、ＣＤ−ＲＯＭ（Compact Disc Read-Only Memory ）上の情報を読込むためのＣＤ−ＲＯＭドライブ１０８およびフレキシブルディスク（Flexible Disk、以下ＦＤ）１１６に情報を読み書きするためのＦＤドライブ１０６を備えたコンピュータ本体１０２と、コンピュータ本体１０２に接続された表示装置としてのディスプレイ１０３と、同じくコンピュータ本体１０２に接続された入力装置としてのキーボード１１０およびマウス１１２とを含む。

さらに、図２に示されるように、このコンピュータ１００を構成するコンピュータ本体１０２は、ＣＤ−ＲＯＭドライブ１０８およびＦＤドライブ１０６に加えて、それぞれバスＢＳに接続されたＣＰＵ（Central Processing Unit ）１２０と、ＲＯＭ（Read Only Memory) およびＲＡＭ （Random Access Memory）を含むメモリ１２２と、直接アクセスメモリ装置、たとえば、ハードディスク１２４と、被制御系である脳磁計１２とデータの授受を行うための通信インタフェース４０とを含んでいる。ＣＤ−ＲＯＭドライブ１０８にはＣＤ−ＲＯＭ１１８が装着される。ＦＤドライブ１０６にはＦＤ１１６が装着される。

ハードディスク１２４やメモリ１２２には、コンピュータシステム１００の動作を制御するためのプログラムがインストールされる。また、ハードディスク１２４（または、図示しない外部記憶装置）には、測定結果のデータが格納される。

なお、ＣＤ−ＲＯＭ１１８は、コンピュータ本体に対してインストールされるプログラム等の情報を記録可能な媒体であれば、他の媒体、たとえば、ＤＶＤ−ＲＯＭ（Digital Versatile Disc）やメモリカードなどでもよく、その場合は、コンピュータ本体１０２には、これらの媒体を読取ることが可能なドライブ装置が設けられる。

本発明は、このコンピュータシステム１００でのソフトウェア処理に関するものである。ただし、処理の一部ないし全部は、処理時間の高速化のためにハードウェアにより行われる構成であってもよい。また、このソフトウェアは、特に限定されないが、ＣＤ−ＲＯＭ（Compact Disc Read Only Memory）やＤＶＤ−ＲＯＭ(Digital Versatile Disc Read Only Memory)のような記録媒体上に記録されて、コンピュータシステム１００にインストールされてもよいし、あるいは、通信回線を介して配信されてもよい。

［電流源推定の具体的手続き（階層変分ベイズ推定法）］
以下では、本発明の「脳内電流源推定方法」の手続を説明するための前提として、上述した「国際公開第ＷＯ０３／０５７０３５Ａ１号パンフレット」に開示された内容にしたがって、具体的に「脳内電流源」の位置を特定する手続きについて、さらに詳しく説明しておく。

まず、図３は、センサアレイ２０内のＳＱＵＩＤセンサ２０．ｉの配置を示す図である。

以下の説明では、本発明の動作を説明する前提として、図３に示すとおり、複数のセンサ２０．ｉは、観測曲面３０上で固定された位置に存在するものとして、脳内電流源推定方法の具体的手続を説明する。

［（Ｉ）電流源推定のための準備］
以下に説明する電流源推定のための処理手順を行なうに当たって、まず、その準備となる手続きについて説明する。

すなわち、電流分布の推定のために、以下の手続きを予め行っておく必要がある。

（Ｉ−１） 仮想曲面形状と電流モデルの決定
（Ｉ−２） さらに、仮想曲面を移動させてその仮想的な電流分布を推定するに当たり、仮想曲面上のサンプル点と、深さ方向のサンプル地点の決定を行なう必要がある。

（Ｉ−３） また、以下で詳しく説明するように、電流分布推定のための階層事前分布の分布形状を指定するためのメタパラメータを初期値として予め決定しておく必要がある。

以下では、上述した（Ｉ−１）仮想曲面形状と電流モデルの決定の手続きについて、さらに詳しく説明する。

（Ｉ−１−１） 仮想曲面形状の決定
最も単純な仮想曲面形状は、脳を半球とみなし、半径の異なるいろいろな半球面を仮想曲面とすることにより得られるものである。

核磁気共鳴画像（Magnetic Resonance Imaging:MRI）等の他の測定方法により、大脳皮質等の脳内電流が存在する可能性のある場所がわかっている場合には、それらの情報をもとに仮想曲面形状を決定することもできる。

特に高精度で大脳皮質の形状がわかっていて、かつ、その他の領域には脳内電流源が存在していないと考えられる場合は、以下で説明する深さ方向の探索を省略して大脳皮質点の電流分布推定を行うこともできる。また、一部の領域だけ深さ方向の探索をすることも可能である。

この場合は、深さの異なる仮想曲面の形状が一般に異なっていてもよいが、互いに交わらないように形状を決めておく必要がある。

（Ｉ−１−２） 電流モデルの決定
仮想曲面上の電流モデルとして、ここでは電流双極子モデルを考える。ただし、他の電流モデルを考えることも可能である。

仮想曲面上に適当な格子点（サンプル点）｛Ｘ_n｜ｎ＝１，…，Ｎ｝を考える。各格子点上に電流の強さがｊ_nの電流双極子を考える。このとき、格子点Ｘ_nの電流双極子ｊ_nが脳表面上での観測点Ｙ_i（ｉ＝１，…，Ｉ）に作る磁場は、以下のビオサバールの式で表わされる。

ここで、μは透磁率であり、たとえば、ベクトルＸ_ｎに対し、記号｜Ｘ_ｎ｜は、ベクトルＸ_ｎの絶対値を表す。より正確な式として、脳を導電性溶液がつまった球体とみなした時のサルバスの式（文献：Sarvas J. Phys. Med. Biol. 32, 11-22, 1987）、あるいは脳の細かな構造まで考慮に入れて有限要素法や境界要素法等の数値解を用いることもできる。以下では、説明や表記を簡単にするために、上記のビオサバールの式を用いて説明する。

このとき、仮想曲面上の電流双極子｛ｊ_n｜ｎ＝１，…，Ｎ｝が観測点Ｙ_i（ｉ＝１，…，Ｉ）に作る磁場は、次式で表わされる。

観測点Ｙ_iで観測できる磁場の方向をベクトルＳ_i，ｃ（ｃ＝１，…，Ｃ）とすると、この点での磁場のベクトルＳ_i，ｃ方向成分Ｂ_ｉ，ｃは、次式のようになる。

また、観測磁場として磁場の局所的勾配を測定する場合には、上式のＹ_iによる微分を観測することになる。

格子点Ｘ_n（位置ベクトル）での電流双極子の方向に制限をつける場合には、電流双極子の可能な独立方向をｂ_n，ｄ（ｄ＝１，…，Ｄ）として、電流双極子は以下の式で表わされる。この式において、Ｄ＝３の場合は方向に制限がない場合に対応している。

以上まとめると、仮想曲面上の格子点｛Ｘ_n｜ｎ＝１，…，Ｎ｝での電流双極子が観測点Ｙ_iに作る磁場は次式のように書ける。

ここで、ｊ_n，ｄは、格子点Ｘ_nでの電流双極子の独立な成分である。また、関数Ｇ（ｉ，ｃ；ｎ，ｄ）は格子点Ｘ_nでの電流双極子ｊ_n，ｄが観測点Ｙ_iに作る磁場のＳ_ｉ，ｃ方向成分である。

電流分布推定の問題とは、観測された磁場｛Ｂ_i，ｃ｜ｉ＝１，…，Ｉ；ｃ＝１，…，Ｃ｝から仮想曲面上の電流分布｛ｊ_n，ｄ｜ｎ＝１，…，Ｎ；ｄ＝１，…，Ｄ｝を推定する問題であるといえる。

上述した観測点における磁場の測定値を、行列表現を用いて表わすこととすると、以下のように書くことができる。なお、Ｇはリードフィールド行列と呼ばれる。ビオサバールの式ではなく、サルバスの式、あるいは境界要素法等の数値解を用いる場合には、このリードフィールド行列の値が異なるだけであり、以下の説明のその他の部分は同様に適用することができる。

（電流源確率モデル）
このように電流モデルを決定した上で、さらにこのような電流分布の推定のために、以下に説明する「電流源の確率モデル」を導入する。

（Ｉ−１−３） 電流源確率モデル
上述した仮想曲面上の電流モデルに対する確率モデルを以下のように考える。観測される磁場は仮想曲面上の電流分布Ｊから作られる磁場と観測ノイズの和として表わされるとする。また、観測ノイズは各測定点で独立な分散σ²を持つガウスノイズとして表わされるとする。

より一般的には、各測定点のノイズの相関が共分散行列の形で表される多次元正規分布に従うガウスノイズを考えることもできるが、以下では説明の簡単のために等方等分散のノイズモデルを用いて説明する。

すなわち、観測磁場が以下の式のように表わされるものとする。

この電流モデルに対する確率分布は以下のように与えられる。

まず、ある特定の深さに対する仮想曲面、あるいは各局所面の深さが指定された局所面の集合をＭで表わす。この仮想曲面Ｍ上の電流分布Ｊが与えられたときに、観測される磁場がＢである確率Ｐ（Ｂ｜Ｊ，β，Ｍ）は、以下の式（１）で表わされる。ただし、以下の式において、β＝１／σ²である。

（Ｉ−１−４） 階層事前分布
仮想曲面Ｍ上の電流分布Ｊに対する事前分布として、上述したように、局在条件階層事前分布を仮定することにする。

観測を行なう前の電流分布Ｊに対する事前分布（Ｊが実現する確率）は局在条件階層事前分布の仮定の下で以下の式（２）のように表わされる。

ただし、Γ（…）は、ガンマ分布を表わし、次式のように定義される。次式において、Γ（γ₀）は、ガンマ関数であり、その定義式も併せて示す。

また、αとτは、電流分布Ｊとノイズ逆分散βに対する事前分布をモデル化するために導入されたハイパーパラメータであり、このαとτに対する階層事前分布Ｐ₀（α｜Ｍ）とＰ₀（τ｜Ｍ）は以下のように与えられる。

なお、上式において、αとτに対する階層事前分布Ｐ_０（α｜Ｍ）とＰ_０（τ｜Ｍ）の分布形状を決定するメタパラメータ バーα_０、γ_α０、τ_０、バーγ_τ０については、後に詳しく説明する。なお、変数名の前に「バー」との語句をつけたものは、その変数の頭部に記号「−」がついていることをあらわす。

（Ｉ−１−５） ベイズ推定
磁場のデータＢを観測したときに、電流分布がＪである事後確率分布Ｐ（Ｊ｜Ｂ，Ｍ）は、ベイズの定理を使って以下のように求まる。

ここで、事後確率分布Ｐ（Ｊ，β，α，τ｜Ｂ，Ｍ）は、以下の式で表される。

この事後確率分布を用いて、電流分布の期待値Ｅ［Ｊ｜Ｂ，Ｍ］は以下のように与えられる。

また、Ｐ（Ｂ｜Ｍ）は仮想曲面モデルＭの周辺尤度である。電流源の深さを推定するとき、深さの異なるいくつかのモデルを比べる。このとき、深さに関する事前の知識がないものとする。すなわち、以下では、Ｐ（Ｍ）＝一定であるものとする。

観測データＢが与えられたときに、モデルＭが真のモデルである確率、すなわち、モデル事後確率Ｐ（Ｍ｜Ｂ）は、上の仮定のもとでモデル周辺尤度Ｐ（Ｂ｜Ｍ）に比例し、以下の関係が成り立つ。

（Ｉ−１−６） 変分法的ベイズ推定
確率モデルと階層事前分布が与えられたときに、モデル周辺尤度を解析的に求めることは一般的にはできない。

そこで、モデル周辺尤度を近似的に計算する手法として、変分法的ベイズ法を用いる。モンテカルロ法や、ラプラス近似法等、他の近似方法を使うことも可能である。

「変分法的ベイズ法」について、具体的な手続きは後述するものとして、以下に簡単に説明しておく。

真の事後分布Ｐ（Ｊ，β，α，τ｜Ｂ，Ｍ）を近似的に求めるために、試験事後分布Ｑ（Ｊ，β，α，τ）を考える。

２つの確率分布Ｐ（Ｊ，β，α，τ｜Ｂ，Ｍ）とＱ（Ｊ，β，α，τ）の近さは、次式で表されるようなＫ−Ｌ距離を用いて計算することができる。

Ｋ−Ｌ距離は、２つの分布が等しいときに限りゼロになり、それ以外のときは常に正の値をとる。

上式において試験事後分布Ｑに対する自由エネルギーＦ（Ｑ）を定義すると、以下の不等式が得られることになる。

すなわち、自由エネルギーＦ（Ｑ）を最大にする分布Ｑ（Ｊ，β，α，τ）は、真の事後分布Ｐ（Ｊ，β，α，τ｜Ｂ，Ｍ）に等しくなり、このとき自由エネルギーは対数周辺尤度に等しくなる。

変分法的ベイズ法では、Ｑ（Ｊ，β，α，τ）の関数形を以下の形に制限して自由エネルギーの最大化を行なう。

上式の右辺第２項Ｑαを固定して、右辺第１項Ｑ_Jに関してＦ（Ｑ）を最大化するステップと、第１項Ｑ_Jを固定して、第２項Ｑαに関してＦ（Ｑ）を最大化するステップを交互に繰返すことにより、自由エネルギーＦ（Ｑ）を極大にする分布Ｑ^*が得られることになる。

［（ＩＩ）電流源推定の手順］
以上のような準備の下に、以下電流源を推定する具体的な手続きについて、図にしたがって説明する。

図４は、本発明の脳内電流源の位置の推定方法の手順の全体的な流れを示すフローチャートである。

図４を参照して、まず、脳内電流源の位置の推定処理が開始されると（ステップＳ１００）、上述したように、電流分布推定のための階層事前分布の分布形状を指定するためのメタパラメータの値と推定すべき変数の初期値を決定する（ステップＳ１０２）。

続いて、初期分解能を用いた電流源の初期推定が行われ、電流源の候補の抽出が行われる（ステップＳ１０４）。

さらに、このようにして抽出された電流源の候補のうち、まず、最も脳表面に近い電流源の位置の特定が行われる（ステップＳ１０６）。その後、他の電流源の深さの特定が順次行われる（ステップＳ１０８）。

以上のようにして、初期分解能による電流源の位置の特定が行われた後、空間分解能を上げて、電流源の位置の再推定が行われ（ステップＳ１１０）、処理が終了する（ステップＳ１１２）。

以下、図４で説明した各ステップの処理について、さらに詳しく説明する。

（ＩＩ−１） 初期分解能を用いた電流源の初期推定（電流源の候補の抽出）
図５は、図４で説明した処理のうち、ステップＳ１０４の処理、すなわち、初期分解能を用いた電流源の初期推定の処理を説明するためのフローチャートである。

図５を参照して、まず、初期分解能での深さ方向のサンプル地点を｛Ｒ_k｜ｋ＝１，…，Ｋ｝とする。各深さＲ_kにおける仮想曲面に対して電流分布の推定を行なう。さらに、このようにして推定された電流分布に基づいて、各深さＲ_kに対する事後確率を求める。この事後確率が後に説明する自由エネルギー値に相当する（ステップＳ２０２）。

以上のようにして求められた事後確率が最大になる深さＲ_Mでの電流分布に対して、電流の強さの極大点を求める（ステップＳ２０４）。

さらに、各極大点を含む局所面を決定する。局所面の数をＬ個とすると、各局所面が局在化された電流源の候補となる（ステップＳ２０６）。

以下、図４中のステップＳ１０２の処理および、それに続く図５中のステップＳ２０２の処理について、さらに詳しく説明する。

（ＩＩ−１−１）電流分布推定のための階層事前分布の分布形状を指定するためのメタパラメータの値と推定すべき変数の初期値の決定
上述したとおり、電流源の推定は、変分法的ベイズ法を用いて行われる。したがって、図４のステップＳ１０２の処理は、このような変分法的ベイズ法の適用にあたって、以下のように行われる。

階層事前分布の分布形状を指定するメタパラメータを、たとえば、以下のように決定する。

さらに、上述した式に基づいて、各変数の初期設定を以下にように行なう。

（ＩＩ−１−２）電流分布推定のための変分法的ベイズ法の具体的処理
（１） パラメータＪ，βの期待値の計算（Ｊ−ステップの処理）
ここでは、パラメータＪ，βの期待値の計算を行なう。この処理により、Ｑ_Jに関して自由エネルギーＦ（Ｑ）を最大化する処理が行なわれる。

対角行列Ａを以下のように定義することで、Ｑ_Jが以下の式に基づいて導出される。なお、以下の式において、変数の期待値には、その変数名の上に、”−”をつけて表現している。

（２） ハイパーパラメータ期待値、すなわち、パラメータα，τの期待値の計算（α−ステップの処理）
Ｊステップに続いて、ハイパーパラメータα，τの期待値の計算を行なう。すなわち、この処理では、Ｑαに関して自由エネルギーＦ（Ｑ）を最大化する処理を行なうことになる。

この手続は以下の式のように表わされる。

（３） 自由エネルギーの計算
以上説明したようなＪステップおよびαステップに基づいて計算されたＱ_JおよびＱαを用いて、自由エネルギーを以下のようにして計算する。

このようにして、自由エネルギーＦの値が収束するまで、上述したＪステップの処理から自由エネルギーの計算までの処理を繰返す。収束後の自由エネルギーＦ（Ｑ）の値が、対数周辺尤度ｌｏｇＰ（Ｂ｜Ｍ）の近似値を与える。

また、対数モデル事後確率は、ｌｏｇＰ（Ｂ｜Ｍ）と定数分だけの違いしかないので、モデル事後確率が最大になるモデルは、上述したような近似の下では自由エネルギーが最大になるモデルに等しいことになる。

以上のような手続により、深さＲ_kに対する事後確率を求めることができる。これに応じて、上述した図５中のステップＳ２０４およびＳ２０６の処理を行なうことにより、各局所面に局在化された電流源の候補を求めることができる。

（ＩＩ−２） 最も脳表面に近い電流源の特定
次に、図４のステップＳ１０６の処理、すなわち、以上のようにして求められた電流源の候補から、まず、最も脳表面に近い電流源の特定を行なう処理について説明する。

図６は、最も脳表面に近い電流源の特定を行なう処理を説明するためのフローチャートである。

まず、初期値として変数ｌの値を１とする（ステップＳ３０２）。つづいて、変数ｌと、変数ｌのとり得る最大値Ｌとの比較が行われ（ステップＳ３０４）、変数ｌが最大値Ｌ以下であれば、処理はステップＳ３０６に移行し、変数ｌが最大値Ｌを超えていれば、処理はステップＳ３２４に移行する（ステップＳ３０４）。すなわち、ステップＳ３０６からＳ３２２までの処理を、ｌ＝１からＬまで繰返す。

ステップＳ３０６では、ｌ番目の局所面以外の局所面の深さを図５のステップＳ２０４で求められた深さＲ-_Mに固定する。

変数ｋの値を１とする（ステップＳ３０８）。つづいて、変数ｋと、変数ｋのとり得る最大値Ｋとの比較が行われ（ステップＳ３１０）、変数ｋが最大値Ｋ以下であれば、処理はステップＳ３１２に移行し、変数ｋが最大値Ｋを超えていれば、処理はステップＳ３２０に移行する。

ステップＳ３１２では、ｌ番目の局所面の深さをまずＲ_kとする。

次に、Ｌ個の局所面の集合に対して電流分布の推定を行なう（ステップＳ３１４）。

さらに、この局所面配置に対する事後確率（自由エネルギー）を求める（ステップＳ３１６）。変数ｋの値を１だけインクリメントして（ステップＳ３１８）、処理はステップＳ３１０に復帰する。

このようにして、ｋ＝１からｋ＝Ｋまでの深さの各局所面についての処理を行なった後、事後確率が最大になるｌ番目の局所面の深さＲ_M（ｌ）を求める（ステップＳ３２０）。変数ｌの値を１だけインクリメントして（ステップＳ３２２）、処理はステップＳ３０４に復帰する。

このようにして、各変数ｌについて求めた局所面の深さＲ_M（ｌ）の中で、最も脳表面に近い（浅い）ｌを求め、これをｌ₁とする（ステップＳ３２４）。

以上の処理により、ｌ₁番目の局所面が最も脳表面に近い電流源に対応し、その深さの初期推定値がＲ_M（ｌ_１）として求められ、処理が終了する（ステップＳ３２６）。

（ＩＩ−３） 各局所面に対応する電流源の深さの特定
図７および図８は、図４のステップＳ１０８の処理、すなわち、各局所面に対応する電流源の深さの特定処理を説明するためのフローチャートである。

図７を参照して、初期値として変数ｓの値を１とする（ステップＳ４０２）。つづいて、変数ｓと、変数ｓのとり得る最大値Ｌとの比較が行われ（ステップＳ４０４）、変数ｓが最大値Ｌ未満であれば、処理はステップＳ４０６に移行し、変数ｓが最大値Ｌ以上であれば、処理はステップＳ４３４に移行する（ステップＳ４０４）。すなわち、ステップＳ４０６からＳ４３２までの処理を、ｓ＝１からＬまで繰返す。

ステップＳ４０６では、これまでに特定された局所面の深さを｛Ｒ_M（ｌ₁），…，Ｒ_M（ｌ_s）｝とするとき、これらの局所面の深さを、それぞれ｛Ｒ_M（ｌ₁），…，Ｒ_M（ｌ_s）｝に固定する（ステップＳ４０６）。

まず、ｌが｛ｌ₁，…，ｌ_s｝のいずれとも等しくなく、かつｌは｛１，…，Ｌ｝の集合に属しているものとする（ステップＳ４０８）。次に、変数ｌの値として、可能な値を全て尽くしたか否かを判断し（ステップＳ４１０）、全て尽くしている場合は処理はステップＳ４３０に移行する。一方、全てを尽くしていない場合は、処理はステップＳ４１２に移行する。

ステップＳ４１２では、ｌ番目の局所面以外で、かつ｛ｌ₁，…，ｌ_s｝でもない局所面の深さをＲ_Mに固定する。

変数ｋの値を１とする（ステップＳ４１４）。つづいて、変数ｋと、変数ｋのとり得る最大値Ｋとの比較が行われ（ステップＳ４１６）、変数ｋが最大値Ｋ以下であれば、処理はステップＳ４１８に移行し、変数ｋが最大値Ｋを超えていれば、処理はステップＳ４２６に移行する。すなわち、ステップＳ４１８からＳ４２４までの手続を、ｋ＝１からｋ＝Ｋまで繰返す。

ステップＳ４１８では、ｌ番目の局所面の深さをＲ_kにする。

次に、Ｌ個の局所面の集合に対して電流分布推定を行なう（ステップＳ４２０）。

さらに、この局所面配置に対する事後確率（自由エネルギー）を求める（ステップＳ４２２）。

図８を参照して、ステップＳ４２６では、以上のようにしてｋ＝Ｋまでの処理が終わった後に、事後確率が最大になるｌ番目の局所面の深さＲ_M（ｌ）を求める（ステップＳ４２６）。

｛ｌ₁，…，ｌ_s｝のいずれとも等しくなく、かつ｛１，…，Ｌ｝の集合に属している、まだ処理を行っていない他の変数ｌを選択し（ステップＳ４２８）、ステップＳ４１０に復帰する。

以上の手続により、処理を行った全ての変数ｌに対応するＲ_M（ｌ）の中で、脳表面に最も近いｌを求め、これをｌ_s+1とする（ステップＳ４３０）。

さらに、ｓにｓ＋１を代入し（ステップＳ４３２）、処理はステップＳ４０４に復帰する。

全ての可能なｓについて処理が終了することで、各局所面に対応する電流源の深さの特定が終了する（ステップＳ４３４）。

（ＩＩ−４） 分解能を上げて電流源の再推定
図９および図１０は、図４のステップＳ１１０の処理、すなわち、空間分解能を上げて電流源の位置を再推定する処理を説明するためのフローチャートである。

図９および図１０を参照して、まず、初期分解能を用いて推定した電流源と対応する局面の番号を脳表面に近いものから１，２，３，…，Ｌとなるように番号の並べ替えを行なう（ステップＳ５０２）。次に、各局所面での電流分布の拡がりに対応して局所面を小さくする（ステップＳ５０４）。

図４のステップＳ１０８までの処理で求めた局所面の深さＲ_M（ｌ）をＲ_M ^old（ｌ）とする（ステップＳ５０６）。

新しい深さ方向の分解能と探索幅を、それぞれΔＲおよび（Ｋ_L・ΔＲ）とする。また局所面上での格子点の分解能も上げる（ステップＳ５０８）。

初期値として変数ｌの値を１とする（ステップＳ５１０）。つづいて、変数ｌと、変数ｌのとり得る最大値Ｌとの比較が行われ（ステップＳ５１２）、変数ｌが最大値Ｌ以下であれば、処理はステップＳ５１４に移行し、変数ｌが最大値Ｌを超えていれば、処理はステップＳ５３４に移行する（ステップＳ５１２）。すなわち、ステップＳ５１４からＳ５３２までの処理を、ｌ＝１からｌ＝Ｌまで繰返す。

ステップＳ５１４では、ｌ番目以外の局所面ｌ′の深さをＲ_M ^old（ｌ′）に固定する。

その上で、さらに、ｋ＝０，±１，…，±Ｋ_Lについて、以下のステップＳ５１８からステップＳ５２６までの処理を行なう。

まず、変数ｋの値を０とする（ステップＳ５１６）。つづいて、変数ｋの絶対値と、変数ｋの絶対値のとり得る最大値Ｋ_Lとの比較が行われ（ステップＳ５１８）、変数ｋの絶対値が最大値Ｋ_L以下であれば、処理はステップＳ５２０に移行し、変数ｋの絶対値が最大値Ｋ_Lを超えていれば、処理はステップＳ５２８に移行する（ステップＳ５１８）。

ステップＳ５２０では、ｌ番目の局所面の深さを（Ｒ_M ^old（ｌ）＋ｋ・ΔＲ）とする。

次に、Ｌ個の局所面の集合に対して、現在の分解能を用いて、電流分布推定を行なう（ステップＳ５２２）。

さらに、この局所面配置に対する事後確率（自由エネルギー）を求める（ステップＳ５２４）。その上で、ｋの値を｛±１，…，±Ｋ_L｝のうちの次の値とし、処理はステップＳ５１８に復帰する。

以上のような処理がｋ＝±Ｋ_Lまで行なわれた後に、ステップＳ５２８では、事後確率が最大となるｋを求め、これをｋ_Mとする。

その上で、Ｒ_M ^old（ｌ）にＲ_M ^old（ｌ）＋ｋ_M・ΔＲを代入する（ステップＳ５３０）。変数ｌの値を１だけインクリメントして（ステップＳ５３２）、処理はステップＳ５１２に復帰する。

以上のような処理が変数ｌの値が最大値Ｌとなるまで行われた時点で、分解能が最終分解能に達していれば（ステップＳ５３４）、処理を終了する（ステップＳ５３８）。

一方、分解能が最終分解能に達していなければ（ステップＳ５３４）、局所面上での格子点の分解能と深さ方向の分解能を上げて（ステップＳ５３６）、処理は、ステップＳ５１０に復帰する。

以上説明したような「脳内電流源推定方法」を用いれば、ＭＥＧ（またはＥＥＧ）の観測データから脳内電流源の位置を深さ方向まで含めて推定することが可能となる。さらに、このような深さ方向の推定は、複数の電流源が存在する場合にも適用可能である。しかも、電流源が電流双極子のように局在している場合にも、広がりを持つ電流源の場合にも適用できるだけでなく、電流源の広がり具合を推定することも可能である。

また、図９および図１０で説明したような処理を、初期分解能での推定の後に追加して行えば、位置の推定分解能を逐次的に上げることが可能となる。このことは、また、生理学的な知見や他の観測方法による観測データなどにより探索範囲を限定して調べることが可能なことも意味する。

［本発明の処理：センサ位置を観測曲面上で複数回移動させて測定する場合］
図１１は、本発明の脳磁計システム１０００において、センサアレイ２０を観測曲面３０に沿って指定距離だけ動かした場合のＳＱＩＤセンサ２０．ｉの配置を説明するための概念図である。

図１１においては、初期位置でのＳＱＩＤセンサ２０．ｉの配置を実線で示し、コンピュータシステム１００からの指示により、第１の指定距離だけ動かした場合のＳＱＩＤセンサ２０．ｉの配置を点線で示し、初期位置から第２の指定距離だけ動かした場合のＳＱＩＤセンサ２０．ｉの配置を破線で示している。

本発明では、測定をするにあたり、センサアレイ２０をある位置に設定して、たとえば、Ｎ回分だけＭＥＧ時系列信号の変化パターンを複数回測定し、各測定データを記憶装置に保存する。そののち、センサアレイ２０を観測曲面３０上で指定距離だけ移動させ、再び、Ｎ回分だけＭＥＧ時系列信号の変化パターンを測定して保存する。このようなセンサアレイの移動をＫ回行なって、各センサアレイ位置でのＭＥＧ時系列信号の複数回の測定を繰り返し、最終的に、Ｎ×Ｋ回分のＭＥＧ時系列信号の変化パターンを測定して保存することで、観測データを得る、という手続を行なう。

なお、各移動ごと繰り返す測定回数を、説明の簡単のために同じ回数のＮ回であるとして説明するが、必ずしも移動ごとに測定回数が異なったとしても、本発明は同様に適用可能なものである。

図１２は、本発明において、脳磁計システム１０００により観測されるデータの処理の概念を示す図である。

すなわち、階層ベイズ法に最適なように脳磁計１２において、センサアレイ２０が観測曲面に沿って移動するごとにＭＥＧデータを複数回測定する手続を繰り返し測定手続１、繰り返し測定手続２、繰り返し測定手続３と表現している。これら、繰り返し測定手続１〜３で得られたすべての測定データを、最終的に統合して推定することとしている。

上述したとおり、従来の全脳型脳磁計（ＭＥＧ装置）はセンサの数が数百個程度である。一方、ミリメートル単位の高空間分解能で脳内電流分布を逆推定するには、たとえば、２００００個程度の未知パラメータを推定しなければならず、推定問題は高度の不良設定問題になっている。上述したようなセンサアレイ２０を固定して行なう「階層変分ベイズ推定法」は、ＭＥＧ時系列データから平均電流強度を推定して、電流強度が強いところに逆フィルタをオートフォーカスさせることで実質的なモデルの自由度を下げている。

しかし、非常に多くの脳内領域が同時に活動する時にはどうしても推定精度が落ちてしまう。一方、センサ数を増やすためにセンサを超小型化しようとするとセンサの感度が落ちてしまう。このため、超伝導ＳＱＩＤセンサー技術ではセンサ数を数百個以上に上げることは困難である。

そこで、本発明では、センサアレイ２０を計測ごとにアクティブに動かすことが出来る構成とし、実質的なセンサ数を２００００個のレベルにまで引き上げ高空間分解能での推定精度を大幅に向上させることが可能となる。

従来のＭＥＧ推定では同じ測定を１００回程度繰り返して測定し、１００回のデータの加算平均を取ることでＳＮ比を上げている。Ｗｉｅｎｅｒフィルタ法などの線形フィルタ法を使う場合、各時刻で独立に推定を行うため、センサ位置を毎回変えて１００回測定すると個々のデータのＳＮ比が下がってしまう。このためセンサ位置を変えて測定しても推定精度の向上は余り望めない。

しかし、階層変分ベイズ推定法を使う場合、たとえば、１００ミリ秒間のＭＥＧ時系列データをもとに平均電流強度を推定するので、ＭＥＧデータの加算平均をしなくても平均電流強度の推定精度が確保できる。また、階層変分ベイズ推定法は、センサ位置の異なるデータをセンサ空間ではなく電流空間で平均化できるため、毎回センサ位置を変えて測定したデータを用いて電流推定精度を上げることが可能である。

図１３は、観測曲面上のセンサ配置と推定対象の電流格子点との位置関係を示す図である。

つまり、図１３に示すように、半球状の観測曲面３０上にセンサ２０．ｉが１８８個配置されており、推定の対象となる電流格子点が、センサ個数よりも多い３１１個のような場合でも、良好に電流推定を行うことが出来る。

また、階層変分ベイズ推定法とセンサ位置のアクティブ制御を組み合わせることで、従来のＭＥＧ装置等では実現できなかった高精度の脳活動ムービーを作成することも可能になる。

（センサ位置を観測曲面上で複数回移動させて測定する場合の推定手続）
以下では、繰り返し行なう測定の各測定時刻ｔをあらわに計算式に持ち込んで表現することにする。

上述したように、センサ位置のアクティブ制御を組み合わせることで、繰り返し測定したデータが得られれば、基本的には、センサ位置を固定した場合と同様に観測データから電流源を推定する手続は同様なものとなる。

まず、ある実験条件での脳内電源分布の時系列を以下の式で表す。

また、Ｋ回センサアレイ２０の位置を変え、各センサ位置ごとにＮ回測定した場合の各観測磁場を以下の式で表す。

このとき、上述した式（１）に相当して、電流分布Ｊ（１：Ｔ）が与えられたときに、Ｋ回×Ｎ回の観測で観測される磁場がＢ（１：Ｔ，１：Ｋ，１：Ｎ）である確率Ｐは、以下の式で表される。

なお、（…）´は転置を表し、式（３）の各記号は以下の意味を表す。

ここで、Ｇ（ｉ）（１≦ｉ≦Ｋ）は、センサを移動させた各位置ごとのリードフィールド行列であり、Σ_G(k)は、センサを移動させた各位置ごとの電流Ｊに対する共分散行列である。

上述したように、等方等分散の仮定の下では、この共分散行列は、以下のようになる。

ここで、Ａは、｛α_n｝を対角成分に持つ対角行列である。

したがって、式（３）と、上記式（１）および（２）とを比べると、センサアレイの位置を動かしつつ、測定を行なった場合の式（３）は、時間平均をとれば、上記式（１）および（２）と全く同一のものとなることが分かる。

したがって、以後の計算手続は、センサアレイ２０の位置を固定した場合に行なう、計算手続と全く同様に行なうことが可能であるので、以後の手続の説明は省略する。

このような手続により、ＭＥＧやＥＥＧの観測データから脳内電流源の位置を深さ方向まで含めて、高空問分解能で推定することが可能である。しかも、このような深さ方向の推定は、複数の電流源が存在する場合にも適用可能である。

なお、上述した特許文献１では、上述したような「局在条件事前分布」以外にも、「局所条件と連続条件を組合せた階層事前分布」に対しても、「階層的変分ベイズ推定法」を適用する場合が開示されている。本願発明のようにセンサアレイ位置をアクティブに変化させて観測データを得て、推定を行う手法は、このような「局所条件と連続条件を組合せた階層事前分布」を用いた場合にも適応可能である。したがって、この場合は、電流源が電流双極子のように局在している場合にも、広がりを持つ電流源の場合にも適用できるだけでなく、電流源の広がり具合を推定することも可能である。

今回開示された実施の形態はすべての点で例示であって制限的なものではないと考えられるべきである。本発明の範囲は上記した説明ではなくて特許請求の範囲によって示され、特許請求の範囲と均等の意味および範囲内でのすべての変更が含まれることが意図される。

脳磁計システムの構成の一例を示す概念図である。 コンピュータシステム１００の構成をブロック図形式で示す図である。 センサアレイ２０内のＳＱＵＩＤセンサ２０．ｉの配置を示す図である。 脳内電流源推定方法の手順の全体的な流れを示すフローチャートである。 初期分解能を用いた電流源の初期推定の処理を説明するためのフローチャートである。 最も脳表面に近い電流源の特定を行なう処理を説明するためのフローチャートである。 各局所面に対応する電流源の深さの特定処理を説明するための第１のフローチャートである。 各局所面に対応する電流源の深さの特定処理を説明するための第２のフローチャートである。 空間分解能を上げて電流源の位置を再推定する処理を説明するための第１のフローチャートである。 空間分解能を上げて電流源の位置を再推定する処理を説明するための第２のフローチャートである。 本発明の脳磁計システム１０００において、センサアレイ２０を観測曲面３０に沿って指定距離だけ動かした場合のＳＱＩＤセンサ２０．ｉの配置を説明するための概念図である。 本発明において、脳磁計システム１０００により観測されるデータの処理の概念を示す図である。 観測曲面上のセンサ配置と推定対象の電流格子点との位置関係を示す図である。 電流源が作る電磁場を適当な曲面上で観測している場合の概念図である。 脳内の電流源と複数の仮想曲面との関係を示す概念図である。

符号の説明

１０ 被験者、１２ 脳磁計、２０ センサアレイ、１００ コンピュータシステム。

Claims (6)

1. 脳内に存在する電磁場の源となる電流源の位置を推定するための方法であって、
   観測対象の頭蓋外部において、観測曲面上に配置された複数の観測センサを、前記観測曲面に沿って指定距離だけ移動させ複数のセンサ位置において各々複数回にわたって前記電磁場の時系列データを測定して準備するステップと、
   前記脳内に脳表面からの深さが互いに異なり、かつ互いに交わらない形状を有する複数の仮想曲面と、各前記仮想曲面上に格子点を設定するステップと、
   準備された前記電磁場の前記時系列データについて時間平均した平均電磁場が観測される確率を前記脳内の電流分布により表現したときに、各前記仮想曲面上において、それぞれ前記平均電磁場を復元するための電流分布を階層変分ベイズ推定により推定するステップと、
   前記複数の仮想曲面上において推定された前記電流分布の広がりと、前記電流分布に基づいて復元される電磁場と前記観測された電磁場との誤差とに基づいて、前記複数の仮想曲面のうちから前記広がりと誤差とが極小となる仮想曲面を、真の前記電流源の存在する曲面として特定するステップとを備える、脳内電流源推定方法。
2. 前記電流分布を推定するステップは、
   事前分布と前記電磁場の観測データとから、ベイズ推定法により事後確率を求めるステップを含み、
   前記真の電流源の存在する曲面として特定するステップは、
   前記仮想曲面のうち、対応する前記事後確率が最大となる仮想曲面を特定するステップを含む、請求項１記載の脳内電流源推定方法。
3. 脳内に存在する電磁場の源となる電流源の位置を推定するコンピュータのためのプログラムであって、
   観測対象の頭蓋外部において、観測曲面上に配置された複数の観測センサを、前記観測曲面に沿って指定距離だけ移動させ複数のセンサ位置において各々複数回にわたって前記電磁場の時系列データを測定して記憶装置に格納するステップと、
   前記脳内に脳表面からの深さが互いに異なり、かつ互いに交わらない形状を有する複数の仮想曲面と、各前記仮想曲面上に格子点を設定するステップと、
   準備された前記電磁場の前記時系列データについて時間平均した平均電磁場が観測される確率を前記脳内の電流分布により表現したときに、各前記仮想曲面上において、それぞれ前記平均電磁場を復元するための電流分布を階層変分ベイズ推定により推定するステップと、
   前記複数の仮想曲面上において推定された電流分布の広がりと、前記電流分布に基づいて復元される電磁場と前記観測された電磁場との誤差とに基づいて、前記複数の仮想曲面のうちから前記広がりと誤差とが極小となる仮想曲面を、真の前記電流源の存在する曲面として特定するステップとをコンピュータに実行させるためのプログラム。
4. 前記電流分布を推定するステップは、
   事前分布と前記電磁場の観測データとから、ベイズ推定法により事後確率を求めるステップを含み、
   前記真の電流源の存在する曲面として特定するステップは、
   前記仮想曲面のうち、対応する前記事後確率が最大となる仮想曲面を特定するステップを含む、請求項３記載のプログラム。
5. 脳内に存在する電磁場の源となる電流源の位置を推定するための脳内電流源推定装置であって、
   電磁場の観測データを格納するための記憶装置と、
   観測曲面上に配置され、前記電磁場を測定するための複数の観測センサと、
   観測対象の頭蓋外部において、前記複数の観測センサを、前記観測曲面に沿って指定距離だけ移動させるセンサ移動手段と、
   複数の前記指定距離だけ前記複数のセンサを移動させ、複数のセンサ位置において各々複数回にわたって前記電磁場の時系列データを測定して前記記憶装置に格納する測定制御手段と、
   前記脳内に脳表面からの深さが互いに異なり、かつ互いに交わらない形状を有する複数の仮想曲面と、各前記仮想曲面上に格子点を設定する仮想曲面設定手段と、
   準備された前記電磁場の前記時系列データについて時間平均した平均電磁場が観測される確率を前記脳内の電流分布により表現したときに、各前記仮想曲面上において、それぞれ前記平均電磁場を復元するための電流分布を階層変分ベイズ推定により推定する電流分布推定手段と、
   前記複数の仮想曲面上において推定された電流分布の広がりと、前記電流分布に基づいて復元される電磁場と前記観測された電磁場との誤差とに基づいて、前記複数の仮想曲面のうちから前記広がりと誤差とが極小となる仮想曲面を、真の前記電流源の存在する曲面として特定する電流源特定手段とを備える、脳内電流源推定装置。
6. 前記電流分布推定手段は、
   事前分布と前記電磁場の観測データとから、ベイズ推定法により事後確率を求める事後確率導出手段を含み、
   前記電流源特定手段は、
   前記仮想曲面のうち、対応する前記事後確率が最大となる仮想曲面を特定する仮想曲面特定手段を含む、請求項５記載の脳内電流源推定装置。

Images