(57)【要約】電流源と観測面の間にある仮想曲面上に適当な電流を流すことによって、電流源が発生する電磁場を復元することができ、またこの曲面が真の電流源に近づくほど曲面上の電流分布の広がりが小さくなることに基づいて、観測データを復元する電流源に対するベイズ推定を行う。この推定の際に、モデル事後確率は曲面か電流源を含むときに最大になること、すなわちモデル事後確率を調べることにより電流源の位置を深さ方向まで含めて推定する。また、他の観測手段による観測データが同時に利用できる場合に、それらの情報をベイズ推定を行う際の階層事前分布に組み入れ、より精度の高い推定を行う。
2003057035.tif 000002
【特許請求の範囲】
【請求項１】
頭蓋外部において観測された電磁場に基づいて、脳内に存在する前記電磁場の源となる電流源の位置を推定するための方法であって、
前記脳内に脳表面からの深さが互いに異なり、かつ互いに交わらない形状を有する複数の仮想曲面と、各前記仮想曲面上に格子点を設定するステップと、
各前記仮想曲面上において、それぞれ前記観測された電磁場を復元するための電流分布を推定するステップと、
前記複数の仮想曲面上において推定された電流分布の広がりと、前記電流分布に基づいて復元される電磁場と前記観測された電磁場との誤差とに基づいて、前記複数の仮想曲面のうちから前記広がりと誤差とが極小となる仮想曲面を、真の前記電流源の存在する曲面として特定するステップとを備える、脳内電流源推定方法。
【請求項２】
前記電流分布を推定するステップは、
事前分布と前記電磁場の観測データとから、ベイズ推定法により事後確率を求めるステップを含み、
前記真の電流源の存在する曲面として特定するステップは、
前記仮想曲面のうち、対応する前記事後確率が最大となる仮想曲面を特定するステップを含む、請求項１記載の脳内電流源推定方法。
【請求項３】
前記電流分布を推定するステップは、
前記複数の仮想曲面のうち、脳表面側の仮想曲面から順次深い側の仮想曲面に移行するのに応じて、前記脳表面に最も近く、かつ事後確率が極大になる第１の仮想曲面を特定するステップを含み、
前記真の電流源の存在する曲面として特定するステップは、
前記第１の仮想曲面において、前記電流分布の極大点に対応する局在電流分布を特定するステップと、
前記局在電流分布を各々が含む複数の局所面を切り出すステップと、
前記複数の局所面のうち、特定対象となる局所面以外を固定し、前記特定対象となる局所面を深さ方向に動かして、前記事後確率が極大となる位置を、前記脳表面に近い側から順次特定するステップとを含む、請求項２記載の脳内電流源推定方法。
【請求項４】
前記電流分布を推定するステップにおいては、第１の空間分解能で、前記電流分布を推定し、
前記複数の仮想曲面の深さ方向の分解能を高め、かつ、前記第１の空間分解能よりも高い第２の空間分解能で、前記電流分布の再推定を行うステップをさらに備える、請求項３記載の脳内電流源推定方法。
【請求項５】
前記電流分布を推定するステップは、前記電流分布をベイズ推定する際に、前記電流源の推定のための前記電磁場の観測とは独立した他の観測方法による観測データを用いて、前記ベイズ推定における階層事前分布を設定するステップを含む、請求項１記載の脳内電流源推定方法。
【請求項６】
頭蓋外部において観測された電磁場に基づいて、脳内に存在する前記電磁場の源となる電流源の位置を推定するコンピュータのためのプログラムであって、
前記脳内に脳表面からの深さが互いに異なり、かつ互いに交わらない形状を有する複数の仮想曲面と、各前記仮想曲面上に格子点を設定するステップと、
各前記仮想曲面上において、それぞれ前記観測された電磁場を復元するための電流分布を推定するステップと、
前記複数の仮想曲面上において推定された電流分布の広がりと、前記電流分布に基づいて復元される電磁場と前記観測された電磁場との誤差とに基づいて、前記複数の仮想曲面のうちから前記広がりと誤差とが極小となる仮想曲面を、真の前記電流源の存在する曲面として特定するステップとを実行させるためのプログラム。
【請求項７】
前記電流分布を推定するステップは、
事前分布と前記電磁場の観測データとから、ベイズ推定法により事後確率を求めるステップを含み、
前記真の電流源の存在する曲面として特定するステップは、
前記仮想曲面のうち、対応する前記事後確率が最大となる仮想曲面を特定するステップを含む、請求項６記載のプログラム。
【請求項８】
前記電流分布を推定するステップは、
前記複数の仮想曲面のうち、脳表面側の仮想曲面から順次深い側の仮想曲面に移行するのに応じて、前記脳表面に最も近く、かつ事後確率が極大になる第１の仮想曲面を特定するステップを含み、
前記真の電流源の存在する曲面として特定するステップは、
前記第１の仮想曲面において、前記電流分布の極大点に対応する局在電流分布を特定するステップと、
前記局在電流分布を各々が含む複数の局所面を切り出すステップと、
前記複数の局所面のうち、特定対象となる局所面以外を固定し、前記特定対象となる局所面を深さ方向に動かして、前記事後確率が極大となる位置を、前記脳表面に近い側から順次特定するステップとを含む、請求項７記載のプログラム。
【請求項９】
前記電流分布を推定するステップにおいては、第１の空間分解能で、前記電流分布を推定し、
前記複数の仮想曲面の深さ方向の分解能を高め、かつ、前記第１の空間分解能よりも高い第２の空間分解能で、前記電流分布の再推定を行うステップをさらに備える、請求項８記載のプログラム。
【請求項１０】
前記電流分布を推定するステップは、前記電流分布をベイズ推定する際に、前記電流源の推定のための前記電磁場の観測とは独立した他の観測方法による観測データを用いて、前記ベイズ推定における階層事前分布を設定するステップを含む、請求項６記載のプログラム。
【請求項１１】
頭蓋外部において観測された電磁場に基づいて、脳内に存在する前記電磁場の源となる電流源の位置を推定するための脳内電流源推定装置であって、
前記脳内に脳表面からの深さが互いに異なり、かつ互いに交わらない形状を有する複数の仮想曲面と、各前記仮想曲面上に格子点を設定する仮想曲面設定手段と、
各前記仮想曲面上において、それぞれ前記観測された電磁場を復元するための電流分布を推定する電流分布推定手段と、
前記複数の仮想曲面上において推定された電流分布の広がりと、前記電流分布に基づいて復元される電磁場と前記観測された電磁場との誤差とに基づいて、前記複数の仮想曲面のうちから前記広がりと誤差とが極小となる仮想曲面を、真の前記電流源の存在する曲面として特定する電流源特定手段とを備える、脳内電流源推定装置。
【請求項１２】
前記電流分布推定手段は、
事前分布と前記電磁場の観測データとから、ベイズ推定法により事後確率を求める事後確率導出手段を含み、
前記電流源特定手段は、
前記仮想曲面のうち、対応する前記事後確率が最大となる仮想曲面を特定する仮想曲面特定手段を含む、請求項１１記載の脳内電流源推定装置。
【請求項１３】
前記電流分布推定手段は、
前記複数の仮想曲面のうち、脳表面側の仮想曲面から順次深い側の仮想曲面に移行するのに応じて、前記脳表面に最も近く、かつ事後確率が極大になる第１の仮想曲面を特定する最浅仮想曲面特定手段を含み、
前記電流源特定手段は、
前記第１の仮想曲面において、前記電流分布の極大点に対応する局在電流分布を特定する局在電流分布特定手段と、
前記局在電流分布を各々が含む複数の局所面を切り出す局所面抽出手段と、
前記複数の局所面のうち、特定対象となる局所面以外を固定し、前記特定対象となる局所面を深さ方向に動かして、前記事後確率が極大となる位置を、前記脳表面に近い側から順次特定する局所面位置特定手段とを含む、請求項１２記載の脳内電流源推定装置。
【請求項１４】
前記電流分布推定手段は、初期的に第１の空間分解能で、前記電流分布を推定した後、前記複数の仮想曲面の深さ方向の分解能を高め、かつ、前記第１の空間分解能よりも高い第２の空間分解能で、前記電流分布の再推定を行う、請求項１３記載の脳内電流源推定装置。
【請求項１５】
前記電流分布推定手段は、前記電流分布をベイズ推定する際に、前記電流源の推定のための前記電磁場の観測とは独立した他の観測方法による観測データを用いて、前記ベイズ推定における階層事前分布を設定する手段を含む、請求項１１記載の脳内電流源推定装置。
【請求項１６】
頭蓋外部において観測された電磁場に基づいて、脳内に存在する前記電磁場の源となる電流源の位置を推定するコンピュータのためのプログラムを記録した記録媒体であって、前記プログラムは、
前記脳内に脳表面からの深さが互いに異なり、かつ互いに交わらない形状を有する複数の仮想曲面と、各前記仮想曲面上に格子点を設定するステップと、
各前記仮想曲面上において、それぞれ前記観測された電磁場を復元するための電流分布を推定するステップと、
前記複数の仮想曲面上において推定された電流分布の広がりと、前記電流分布に基づいて復元される電磁場と前記観測された電磁場との誤差とに基づいて、前記複数の仮想曲面のうちから前記広がりと誤差とが極小となる仮想曲面を、真の前記電流源の存在する曲面として特定するステップとを備える記録媒体。
【請求項１７】
前記電流分布を推定するステップは、
事前分布と前記電磁場の観測データとから、ベイズ推定法により事後確率を求めるステップを含み、
前記真の電流源の存在する曲面として特定するステップは、
前記仮想曲面のうち、対応する前記事後確率が最大となる仮想曲面を特定するステップを含む、請求項１６記載の記録媒体。
【請求項１８】
前記電流分布を推定するステップは、
前記複数の仮想曲面のうち、脳表面側の仮想曲面から順次深い側の仮想曲面に移行するのに応じて、前記脳表面に最も近く、かつ事後確率が極大になる第１の仮想曲面を特定するステップを含み、
前記真の電流源の存在する曲面として特定するステップは、
前記第１の仮想曲面において、前記電流分布の極大点に対応する局在電流分布を特定するステップと、
前記局在電流分布を各々が含む複数の局所面を切り出すステップと、
前記複数の局所面のうち、特定対象となる局所面以外を固定し、前記特定対象となる局所面を深さ方向に動かして、前記事後確率が極大となる位置を、前記脳表面に近い側から順次特定するステップとを含む、請求項１７記載の記録媒体。
【請求項１９】
前記電流分布を推定するステップにおいては、第１の空間分解能で、前記電流分布を推定し、
前記複数の仮想曲面の深さ方向の分解能を高め、かつ、前記第１の空間分解能よりも高い第２の空間分解能で、前記電流分布の再推定を行うステップをさらに備える、請求項１８記載の記録媒体。
【請求項２０】
前記電流分布を推定するステップは、前記電流分布をベイズ推定する際に、前記電流源の推定のための前記電磁場の観測とは独立した他の観測方法による観測データを用いて、前記ベイズ推定における階層事前分布を設定するステップを含む、請求項１６記載の記録媒体。
【発明の詳細な説明】
技術分野
本発明は、脳内活動に対応して、頭蓋外部に磁場または電場を生成する脳内電流源の位置や分布を推定する方法、脳内電流源推定プログラム、当該プログラムを記録した記録媒体および脳内電流源推定装置に関する。
背景技術
近年の生体計測技術の進歩はめざましく、従来測定が困難で誤差も大きかった脳から発生する微弱電場（脳波）や微弱磁場（脳磁波）の計測精度は、年々向上している。
すなわち、外部からの刺激を受けて、脳内の神経細胞は電流を発生する。この電流により、上述したような微弱な電場や微弱な磁場が発生することになる。このうち「脳波」とは、この神経細胞の電流により、脳から発生する電場のことである。「脳波計（Ｅｌｅｃｔｒｏｅｎｃｅｐｈａｌｏｇｒａｍ：ＥＥＧ）」は、このような脳波を測定する手法である。
一方、「脳磁波」とは、この神経細胞の電流により、脳から発生する磁場のことである。「脳磁計（Ｍａｇｎｅｔｏｅｎｃｅｐｈａｌｏｇｒａｐｈｙ：ＭＥＧ）」は、このような脳磁波を測定する手法である。脳磁波計測の最大の利点は、磁場が容積導体による影響をほとんど受けないので、頭蓋外からの磁気計測により、脳内電流源位置を３次元的に比較的正確に推定できることが期待されることである。
このような脳磁波の解析においては、発生した磁場を脳の外部から測定することで非侵襲的に脳の活動部位を推定する。
しかしながら、この磁場は非常に微弱であるために、地磁気等の外部の磁場の影響を大きく受けてしまう。そこで、外部の磁場を遮蔽するシールドを作り、その中で超伝導を用いた測定機器である超伝導量子干渉素子（ＳＱＵＩＤ：Ｓｕｐｅｒｃｏｎｄｕｃｔｉｎｇ  ｑｕａｎｔｕｍ  ｉｎｔｅｒｆｅｒｅｎｃｅ  ｄｅｖｉｃｅ）により微弱な磁場を測定することが行われる。
ただし、このような「脳内電流源位置推定」のアルゴリズムの研究においても、初期のモデルから様々な変法が使用されているにも関わらず、決定的な手法は今だ存在しないのが現状である。
たとえば、このような「脳内電流源位置推定」のアルゴリズムの１つである「双極子推定法」が、文献１：Ｍｏｓｈｅｒ．Ｊ．Ｃ．，Ｌｅｗｉｓ  Ｐ．Ｓ．ａｎｄ  Ｌｅａｈｙ  Ｒ．Ｍ．ＩＥＥＥ  Ｔｒａｎｓ．Ｂｉｏｍｅｄ．Ｅｎｇｎｇ．〈１９９２〉ｖｏｌ．３９，ｐｐ．５４１−５５７に開示されている。しかしながら、この「双極子推定法」は、脳内の電流源が１つか数個の電流双極子で表せると仮定して観測磁場から双極子の位置を推定する方法であり、双極子の数を幾つにすれば良いかを決めるのが難しいという欠点を持っている。
一方、他のアルゴリズムとしては、「空間フィルター法」が、文献２：Ｔｏｙａｍａ  Ｋ，Ｙｏｓｈｉｚａｗａ  Ｋ，Ｙｏｓｈｉｄａ  Ｙ，Ｋｏｎｄｏ  Ｙ，Ｔｏｍｉｔａ  Ｓ，Ｔａｋａｎａｓｈｉ  Ｙ，Ｅｊｉｍａ  Ｙ，ａｎｄ  Ｙｏｓｈｉｋａｗａ  Ｓ，Ｎｅｕｒｏｓｃｉｅｎｃｅ．１９９９；９１（２）：ｐｐ．４０５−４１５に開示されている。「空間フィルター法」は、生理学的な知見から、脳内電流源の場所を限定し、双極子の分布を推定する方法である。しかしながら、この方法の欠点は電流源の深さを正確に推定できない事である。
発明の開示
本発明の目的は、観測データから脳内電流源の位置を深さ方向まで含めて推定することが可能な脳内電流源推定方法を提供することである。
この発明の他の目的は、複数の脳内電流源が存在する場合にも、観測データから複数の脳内電流源の位置を推定することが可能な脳内電流源推定方法を提供することである。
この発明のさらに他の目的は、脳内の電流源の推定のための電磁場の観測とは独立した他の観測方法による観測データが同時に利用できる場合に、複数の観測方法によるデータを組み合わせて、より精度の高い推定を行う脳内電流源推定方法を提供することである。
この発明のさらに他の目的は、観測データから脳内電流源の位置を深さ方向まで含めて推定することが可能な脳内電流源推定プログラム、当該プログラムを記録した記録媒体を提供することである。
この発明のさらに他の目的は、複数の脳内電流源が存在する場合にも、観測データから複数の脳内電流源の位置を推定することが可能な脳内電流源推定プログラム、当該プログラムを記録した記録媒体を提供することである。
この発明のさらに他の目的は、脳内の電流源の推定のための電磁場の観測とは独立した他の観測方法による観測データが同時に利用できる場合に、複数の観測方法によるデータを組み合わせて、より精度の高い推定を行う脳内電流源推定プログラム、当該プログラムを記録した記録媒体を提供することである。
この発明のさらに他の目的は、観測データから脳内電流源の位置を深さ方向まで含めて推定することが可能な脳内電流源推定装置を提供することである。
この発明のさらに他の目的は、複数の脳内電流源が存在する場合にも、観測データから複数の脳内電流源の位置を推定することが可能な脳内電流源推定装置を提供することである。
この発明のさらに他の目的は、脳内の電流源の推定のための電磁場の観測とは独立した他の観測方法による観測データが同時に利用できる場合に、複数の観測方法によるデータを組み合わせて、より精度の高い推定を行う脳内電流源推定装置を提供することである。
係る目的を達成するために本願発明に係る脳内電流源推定方法は、頭蓋外部において観測された電磁場に基づいて、脳内に存在する電磁場の源となる電流源の位置を推定するための方法であって、脳内に脳表面からの深さが互いに異なり、かつ互いに交わらない形状を有する複数の仮想曲面と、各仮想曲面上に格子点を設定するステップと、各仮想曲面上において、それぞれ観測された電磁場を復元するための電流分布を推定するステップと、複数の仮想曲面上において推定された電流分布の広がりと、電流分布に基づいて復元される電磁場と観測された電磁場との誤差とに基づいて、複数の仮想曲面のうちから広がりと誤差とが極小となる仮想曲面を、真の電流源の存在する曲面として特定するステップとを備える。
好ましくは、電流分布を推定するステップは、事前分布と電磁場の観測データとから、ベイズ推定法により事後確率を求めるステップを含み、真の電流源の存在する曲面として特定するステップは、仮想曲面のうち、対応する事後確率が最大となる仮想曲面を特定するステップを含む。
好ましくは、電流分布を推定するステップは、複数の仮想曲面のうち、脳表面側の仮想曲面から順次深い側の仮想曲面に移行するのに応じて、脳表面に最も近く、かつ事後確率が極大になる第１の仮想曲面を特定するステップを含み、真の電流源の存在する曲面として特定するステップは、第１の仮想曲面において、電流分布の極大点に対応する局在電流分布を特定するステップと、局在電流分布を各々が含む複数の局所面を切り出すステップと、複数の局所面のうち、特定対象となる局所面以外を固定し、特定対象となる局所面を深さ方向に動かして、事後確率が極大となる位置を、脳表面に近い側から順次特定するステップとを含む。
好ましくは、電流分布を推定するステップにおいては、第１の空間分解能で、電流分布を推定し、複数の仮想曲面の深さ方向の分解能を高め、かつ、第１の空間分解能よりも高い第２の空間分解能で、電流分布の再推定を行うステップをさらに備える。
好ましくは、電流分布を推定するステップは、電流分布をベイズ推定する際に、電流源の推定のための電磁場の観測とは独立した他の観測方法による観測データを用いて、ベイズ推定における階層事前分布を設定するステップを含む。たとえば、他の観測方法による観測データにより脳内電流源が存在する場所が限られた領域内にあることが分かっている場合には、深さ方向の探索を省略して、限定された領域内における電流分布をベイズ推定により求めることができる。
この発明の他の局面に従うと、頭蓋外部において観測された電磁場に基づいて、脳内に存在する電磁場の源となる電流源の位置を推定するコンピュータのためのプログラムであって、脳内に脳表面からの深さが互いに異なり、かつ互いに交わらない形状を有する複数の仮想曲面と、各仮想曲面上に格子点を設定するステップと、各仮想曲面上において、それぞれ観測された電磁場を復元するための電流分布を推定するステップと、複数の仮想曲面上において推定された電流分布の広がりと、電流分布に基づいて復元される電磁場と観測された電磁場との誤差とに基づいて、複数の仮想曲面のうちから広がりと誤差とが極小となる仮想曲面を、真の電流源の存在する曲面として特定するステップとを実行させる。
好ましくは、電流分布を推定するステップは、事前分布と電磁場の観測データとから、ベイズ推定法により事後確率を求めるステップを含み、
真の電流源の存在する曲面として特定するステップは、仮想曲面のうち、対応する事後確率が最大となる仮想曲面を特定するステップを含む。
好ましくは、電流分布を推定するステップは、複数の仮想曲面のうち、脳表面側の仮想曲面から順次深い側の仮想曲面に移行するのに応じて、脳表面に最も近く、かつ事後確率が極大になる第１の仮想曲面を特定するステップを含み、真の電流源の存在する曲面として特定するステップは、第１の仮想曲面において、電流分布の極大点に対応する局在電流分布を特定するステップと、局在電流分布を各々が含む複数の局所面を切り出すステップと、複数の局所面のうち、特定対象となる局所面以外を固定し、特定対象となる局所面を深さ方向に動かして、事後確率が極大となる位置を、脳表面に近い側から順次特定するステップとを含む。
好ましくは、電流分布を推定するステップにおいては、第１の空間分解能で、電流分布を推定し、複数の仮想曲面の深さ方向の分解能を高め、かつ、第１の空間分解能よりも高い第２の空間分解能で、電流分布の再推定を行うステップをさらに備える。
好ましくは、電流分布を推定するステップは、電流分布をベイズ推定する際に、電流源の推定のための電磁場の観測とは独立した他の観測方法による観測データを用いて、ベイズ推定における階層事前分布を設定するステップを含む。たとえば、他の観測方法による観測データにより脳内電流源が存在する場所が限られた領域内にあることが分かっている場合には、深さ方向の探索を省略して、限定された領域内における電流分布をベイズ推定により求めることができる。
この発明のさらに他の局面に従うと、頭蓋外部において観測された電磁場に基づいて、脳内に存在する電磁場の源となる電流源の位置を推定するための脳内電流源推定装置であって、脳内に脳表面からの深さが互いに異なり、かつ互いに交わらない形状を有する複数の仮想曲面と、各仮想曲面上に格子点を設定する仮想曲面設定手段と、各仮想曲面上において、それぞれ観測された電磁場を復元するための電流分布を推定する電流分布推定手段と、複数の仮想曲面上において推定された電流分布の広がりと、電流分布に基づいて復元される電磁場と観測された電磁場との誤差とに基づいて、複数の仮想曲面のうちから広がりと誤差とが極小となる仮想曲面を、真の電流源の存在する曲面として特定する電流源特定手段とを備える。
好ましくは、電流分布推定手段は、事前分布と電磁場の観測データとから、ベイズ推定法により事後確率を求める事後確率導出手段を含み、電流源特定手段は仮想曲面のうち、対応する事後確率が最大となる仮想曲面を特定する仮想曲面特定手段を含む。
好ましくは、電流分布推定手段は、複数の仮想曲面のうち、脳表面側の仮想曲面から順次深い側の仮想曲面に移行するのに応じて、脳表面に最も近く、かつ事後確率が極大になる第１の仮想曲面を特定する最浅仮想曲面特定手段を含み、電流源特定手段は、第１の仮想曲面において、電流分布の極大点に対応する局在電流分布を特定する局在電流分布特定手段と、局在電流分布を各々が含む複数の局所面を切り出す局所面抽出手段と、複数の局所面のうち、特定対象となる局所面以外を固定し、特定対象となる局所面を深さ方向に動かして、事後確率が極大となる位置を、脳表面に近い側から順次特定する局所面位置特定手段とを含む。
好ましくは、電流分布推定手段は、初期的に第１の空間分解能で、電流分布を推定した後、複数の仮想曲面の深さ方向の分解能を高め、かつ、第１の空間分解能よりも高い第２の空間分解能で、電流分布の再推定を行う。
好ましくは、電流分布推定手段は、電流分布をベイズ推定する際に、電流源の推定のための電磁場の観測とは独立した他の観測方法による観測データを用いて、ベイズ推定における階層事前分布を設定する手段を含む。たとえば、他の観測方法による観測データにより脳内電流源が存在する場所が限られた領域内にあることが分かっている場合には、深さ方向の探索を省略して、限定された領域内における電流分布をベイズ推定により求めることができる。
本発明の脳内電流源推定方法に従えば、ＭＥＧやＥＥＧの観測データから脳内電流源の位置を深さ方向まで含めて推定することが可能である。しかも、このような深さ方向の推定は、複数の電流源が存在する場合にも適用可能である。さらに、電流源が電流双極子のように局在している場合にも、広がりを持つ電流源の場合にも適用できるだけでなく、電流源の広がり具合を推定することも可能である。また、初期分解能での推定の後に、位置の推定分解能を逐次的に上げる処理を追加して行うことが可能である。
発明を実施するための最良の形態
以下、図面を参照して本発明の実施の形態について説明する。
なお、上述のとおり、脳内活動を外部から観測する手法としては、脳磁計（Ｍａｇｎｅｔｏｅｎｃｅｐｈａｌｏｇｒａｐｈｙ：ＭＥＧ）や脳波計（Ｅｌｅｃｔｒｏｅｎｃｅｐｈａｌｏｇｒａｍ：ＥＥＧ）がある。以下では、これらのうち、主にＭＥＧについて説明することにする、ただし、ＭＥＧにおける磁場を電場に置き換えればＥＥＧについても、本発明を適用することが可能である。
したがって、本明細書においては、本来「電場」と「磁場」とが複合して存在する「電磁場」に対して、「観測対象」となる物理量が「電場」である場合、または「磁場」である場合とを総称して、「電磁場を観測する」と呼ぶことにする。
［実施の形態１］
図１は、脳磁計システムの構成の一例を示す概念図である。
脳磁計１２には、マルチチャンネルのＳＱＵＩＤ磁束計（超高感度磁気測定器）アレイを備えており、被験者１０の頭蓋表面での磁場を測定する。この脳磁計１２での測定結果を受け取って、コンピュターシステム２０が、測定結果を解析し、脳内の電流源の位置の推定処理を行う。
本発明は、このコンピュータシステム２０でのソフトウェア処理に関するものである。ただし、処理の一部ないし全部は、処理時間の高速化のためにハードウェアにより行われる構成であってもよい。また、このソフトウェアは、特に限定されないが、ＣＤ−ＲＯＭ（Ｃｏｍｐａｃｔ  Ｄｉｓｃ  Ｒｅａｄ  Ｏｎｌｙ  Ｍｅｍｏｒｙ）やＤＶＤ−ＲＯＭ（Ｄｉｇｉｔａｌ  Ｖｅｒｓａｔｉｌｅ  Ｄｉｓｃ  Ｒｅａｄ  Ｏｎｌｙ  Ｍｅｍｏｒｙ）のような記録媒体２２上に記録されて、コンピュータシステム２０にインストールされてもよいし、あるいは、通信回線を介して配信されてもよい。
［脳内電流源推定の原理］
以下、本発明の「脳内電流源推定方法」の詳細について、説明する前に、本発明の推定方法の原理、および概略についてまとめておく。
（曲面上の電流分布による電磁場の復元）
電流源を金属と磁性体とで作られた理想的な電磁シールド面で囲み込めば、シールド面の外側に電磁場が存在しなくなることは電磁シールドの原理として良く知られている．すなわち、シールド面上を流れる電流と磁性体を構成する微小磁石の集まりが作る電磁場により、電流源がシールド面の外側に作る電磁場が完全に打ち消されることになる。また、微小磁石は仮想的な渦電流で置き換えることができる。これは見方を変えれば、このシールド面に適当な電流を流してやれば電流源が作る電磁場と全く同じ電磁場をシールド面の外側に作れることを示している．逆にシールド面の外側にある電流源が作る電磁場はこのシールド面にどのような電流を流しても完全には復元できない．この事実を以後「電磁場復元の原理」と呼ぶ。
図２は、電流源が作る電磁場を適当な曲面上で観測している場合の概念図である。
図２に示すように、観測曲面と電流源の間に仮想曲面を用意すると、上に述べた「電磁場復元の原理」から仮想曲面上に適当な電流を流してやれば電流源が観測面上に作る電磁場を復元できることが分かる。
図３は、脳内の電流源と複数の仮想曲面との関係を示す概念図である。
図３に示すとおり、電流が作る電磁場は距離の二乗に逆比例して減衰することを考慮すれば、電流源と等価な仮想曲面（仮想曲面２）上の電流の広がりは仮想曲面が電流源から遠ざかるほど広くなってゆく。したがって、仮想曲面１上の電流の広がりは、仮想曲面２上の電流の広がりよりも大きい。
一方、電流源に対して観測曲面と反対側にある仮想曲面（仮想曲面３）では電流源が観測面に作る電磁場を完全に復元する事が出来ない。
本発明では、以上の原理を元に観測面上での電磁場の観測データから電流源を推定する。
（電流源推定の原理）
図２に示すように、脳内で発生した電流源から作られる磁場（または電場）を脳の表面近くの観測面で観測していると仮定する。
脳内に仮想曲面を考え、観測された磁場を復元する仮想曲面上の電流分布を求める。この仮想曲面を脳表面から半径を縮めながら脳の内部へ移動させてゆくと、観測磁場を復元する電流分布は、その広がりが小さくなってゆき真の電流源を仮想曲面が含むときに最小になる。また仮想曲面を電流源より深く移動させていくと仮想曲面上での電流分布の広がりは再び広がってゆき、この電流が生ずる磁場と観測磁場との誤差も大きくなってゆく。
すなわち、観測磁場を復元する電流分布の広がりと磁場の復元誤差を見ることにより電流源の深さが特定できる。またこのようにして特定された深さの仮想曲面上の電流分布を求めれば電流源の広がりもわかることになる。以上が本発明の原理である。
（ベイズ推定法による電流源の深さの特定）
上述した電流源推定の原理を具体的に実施するにあたり、本発明では、以下のような確率論に基づく手続きをとる。その手続きの概略を以下にまとめることにする。
ＭＥＧやＥＥＧで観測できるのは脳表面近くの数十カ所から数百カ所の磁場（電場）である。また仮想曲面上の電流分布を近似するために、仮想曲面上に格子点を設定し各格子点に電流双極子（または適当な電流源モデル）を割り当てる。電流分布を高い分解能で推定するためには格子点の数を増やして格子点の密度を大きくする必要がある。
一方、仮想曲面上の格子点の数を観測点の数より増やすと解が一意に決まらない。また観測点より推定点の数の方が多いと、観測データ数よりも推定すべきパラメータの数の方が多くなるので、電流源より深い位置にある仮想曲面でも観測磁場の復元性が良くなる。
このような困難を解決するために、本発明ではベイズ推定理論を利用して仮想曲面上の電流分布を推定する。ベイズ推定を行う際に、電流源に関する事前知識を表す事前分布を用いる。具体的には、脳内電流源は複数の場所に局在していると考えられるので、電流分布の広がりができるだけ小さく推定されるような事前分布、すなわち、観測データだけからではその強さが良く特定できないような格子点上の電流双極子は、その強さがゼロに近くなるような事前分布を導入する．これは有効信号路自動決定（Ａｕｔｏｍａｔｉｃ  Ｒｅｌｅｖａｎｃｅ  Ｄｅｔｅｒｍｉｎａｔｉｏｎ）事前分布（文献：Ｒ．Ｍ．Ｎｅａｌ，Ｂａｙｅｓｉａｎ  Ｌｅａｒｎｉｎｇ  ｆｏｒ  Ｎｅｕｒａｌ  Ｎｅｔｗｏｒｋｓ，Ｓｐｒｉｎｇｅｒ  Ｖｅｒｌａｇ．，１９９６）と呼ばれる階層的事前分布を導入することで実現できる。この事前分布を以後、「局在条件事前分布」と呼ぶことにする。また、局在条件以外の事前知識を導入した事前分布の選び方については実施の形態２において説明をすることにする。
しかし、局在条件事前分布と観測データから事後確率分布を解析的に計算することは不可能である．そこで、本発明では、後に説明するように、事後確率分布を近似的に計算する手法として変分法的ベイズ推定法（Ｖａｒｉａｔｉｏｎａｌ  Ｂａｙｅｓ  ｍｅｔｈｏｄ）（文献：Ａｔｔｉａｓ，Ｈ．Ｐｒｏｃ．１５ｔｈ  Ｃｏｎｆｅｒｅｎｃｅ  ｏｎ  Ｕｎｃｅｒｔａｉｎｔｙ  ｉｎ  Ａｒｔｉｆｉｃｉａｌ  Ｉｎｔｅｌｌｉｇｅｎｃｅ  ｐｐ．２１−３０，１９９９およびＭａｓａ−ａｋｉ  Ｓａｔｏ，Ｎｅｕｒａｌ  Ｃｏｍｐｕｔａｔｉｏｎ，１３，１６４９−１６８１，２００１）を用いる。ただし、モンテカルロ法等の他の近似手法を用いることも可能である。
この局在条件事前分布を用いたベイズ推定法を行うことにより、観測データを復元し、かつ出来るだけ広がりの少ない仮想曲面上の電流分布が得られる。また深さの異なる仮想曲面を使って推定を行ったときのモデル事後確率を比べることにより電流源の深さが特定できる。
モデル事後確率の対数は対数尤度項とモデル複雑度項の和で表される。対数尤度項は復元誤差が小さいほど大きくなる。
一方、モデル複雑度項は有効な（すなわち適当なしきい値以上の強さを持つ）格子点上の電流双極子の数が少ないほど大きくなる．上述したように、復元誤差と仮想曲面上の電流分布の広がり（有効な電流双極子の数）は仮想曲面の深さが電流源と一致したときに最小になるので、この時モデル事後確率が最大になることが分かる。すなわちモデル事後確率が最大になる深さの仮想曲面上に電流源が存在することが分かる。
すなわち、以上をまとめると、本発明は、脳磁波（ＭＥＧ）や脳波（ＥＥＧ）の観測データから脳内電流源の位置を深さ方向まで含めて推定する方法を提供する。
本発明の基本的な考え方は、上述のとおり、電流源と観測面の間にある曲面上に適当な電流を流すことによって、電流源が発生する電磁場を復元することができ、またこの曲面が電流源に近づくほど曲面上の電流分布の広がりが小さくなるという点に基づく。この考え方をベースに、本発明は、観測データを復元する曲面上の電流分布に対するベイズ推定を行う際、モデル事後確率は曲面が電流源を含むときに最大になること、すなわちモデル事後確率を調べることにより電流源の位置を深さ方向まで含めて推定することを１つの特徴とする。
（複数個の電流源の存在する場合）
以上では、電流源が一つの場合について説明したが、この方法は電流源が複数ある場合にも適用できる。
電流により生ずる電磁場は距離の二乗に逆比例して減衰するので、脳表面に最も近い電流源が脳表面の観測磁場に最も影響を与える．そこで脳表面に近い電流源から順番に特定することが出来る。
仮想曲面を脳表面からだんだん深くしてゆくと脳表面に最も近い電流源（これを第１電流源と呼ぶ）の近くでモデル事後確率が極大になる．電流源が２つ以上ある場合にはこの仮想曲面上の電流分布には電流源の数に対応した局所的な電流双極子の集合が複数個出来る。
この局所的な電流双極子の集合を局在電流分布と呼ぶ。個々の局在電流分布を含む局所面を切り出し、これを深さ方向に動かしてモデル事後確率を求める。第１電流源に対応した局所面を動かしたときにはモデル事後確率が第１電流源の深さで極大になる。しかし、これ以外の局所面を動かしてもモデル事後確率は第１電流源の深さで極大になることはない。このことから第１電流源の位置、すなわち最も表面に近い電流源の位置が特定できる。
２番目に深い電流源を求めるために、第１電流源に対応する局所面を固定し、残りの局所面の深さを同じにして深くしてゆく。すると今度は２番目に深い電流源のところでモデル事後確率が極大になる。ここでまた個々の局所面を深さ方向に動かすと第２電流源に対応する局所面を動かしたときにのみ第２電流源の深さでモデル事後確率が極大になる．こうして第２電流源の位置を特定できる。３番目以降の電流源も同様にして順次特定することが出来る。
この方法は個々の局所面の深さを独立に動かして調べる方法に比べて、計算時間がずっと少なくて済むという利点を持つ。
（逐次的に解像度を上げる方法）
以上のようにして、脳表面に近い側から各電流源の位置を特定することが可能であることに加えて、本発明の脳内電流源推定方法では、電流源推定の解像度を逐次的に上げることが出来る。
すなわち、まず低い解像度で（仮想曲面上の格子点の数を少なくし、深さ方向の標本点の数も少なくして）電流源の位置を大雑把に推定する。この段階で各電流源に対応した局所面の位置が大体定まっている。
次により高い分解能で推定を行う。この時、局所面の面積は元の仮想曲面の面積に比べて小さくなっているので、同じ数の格子点を使っても分解能は高くなる。深さ方向の分解能もあげて電流源の位置の推定精度を上げることが出来る。分解能を上げたときに電流分布がより局在する場合には局所面をより小さくして再び分解能を上げて推定できる。
一方、分解能を上げても電流分布の広がりが余り変わらなければ、この電流源が広がっていることを示している。このようにして電流源の広がり具合に応じて分解能を調節することが可能である。
［電流源推定の具体的手続き］
以下では、以上概説した手続きにしたがって、具体的に「脳内電流源」の位置を特定する手続きについて、さらに詳しく説明する。
［（Ｉ）電流源推定のための準備］
以下に説明する電流源推定のための処理手順を行なうに当たって、まず、その準備となる手続きについて説明する。
すなわち、電流分布の推定のために、以下の手続きを予め行っておく必要がある。
（Ｉ−１）  仮想曲面形状と電流モデルの決定
（Ｉ−２）  さらに、仮想曲面を移動させてその仮想的な電流分布を推定するに当たり、仮想曲面上のサンプル点と、深さ方向のサンプル地点の決定を行なう必要がある。
（Ｉ−３）  また、以下で詳しく説明するように、電流分布推定のための階層事前分布の分布形状を指定するためのメタパラメータを初期値として予め決定しておく必要がある。
以下では、上述した（Ｉ−１）仮想曲面形状と電流モデルの決定の手続きについて、さらに詳しく説明する。
（Ｉ−１−１）  仮想曲面形状の決定
最も単純な仮想曲面形状は、脳を半球とみなし、半径の異なるいろいろな半球面を仮想曲面とすることにより得られるものである。
核磁気共鳴画像（Ｍａｇｎｅｔｉｃ  Ｒｅｓｏｎａｎｃｅ  Ｉｍａｇｉｎｇ：ＭＲＩ）等の他の測定方法により、大脳皮質等の脳内電流が存在する可能性のある場所がわかっている場合には、それらの情報をもとに仮想曲面形状を決定することもできる。
特に高精度で大脳皮質の形状がわかっていて、かつ、その他の領域には脳内電流源が存在していないと考えられる場合は、以下で説明する深さ方向の探索を省略して大脳皮質点の電流分布推定を行うこともできる。また、一部の領域だけ深さ方向の探索をすることも可能である。
この場合は、深さの異なる仮想曲面の形状が一般に異なっていてもよいが、互いに交わらないように形状を決めておく必要がある。
（Ｉ−１−２）  電流モデルの決定
仮想曲面上の電流モデルとして、ここでは電流双極子モデルを考える。ただし、他の電流モデルを考えることも可能である。
仮想曲面上に適当な格子点（サンプル点）｛Ｘｎ｜ｎ＝１，…，Ｎ｝を考える。各格子点上に電流の強さがｊｎの電流双極子を考える。このとき、格子点Ｘｎの電流双極子ｊｎが脳表面上での観測点Ｙｉ（ｉ＝１，…，Ｉ）に作る磁場は、以下のビオサバールの式で表わされる。
2003057035.tif 000003
ここで、μは透磁率であり、たとえば、ベクトルＸｎに対し、記号｜Ｘｎ｜は、ベクトルＸｎの絶対値を表す。より正確な式として、脳を導電性溶液がつまった球体とみなした時のサルバスの式（文献：Ｓａｒｖａｓ  Ｊ．Ｐｈｙｓ．Ｍｅｄ．Ｂｉｏｌ．３２，１１−２２，１９８７）、あるいは脳の細かな構造まで考慮に入れて有限要素法や境界要素法等の数値解を用いることもできる。以下では、説明や表記を簡単にするために、上記のビオサバールの式を用いて説明する。
このとき、仮想曲面上の電流双極子｛ｊｎ｜ｎ＝１，…，Ｎ｝が観測点Ｙｉ（ｉ＝１，…，Ｉ）に作る磁場は、次式で表わされる。
2003057035.tif 000004
観測点Ｙｉで観測できる磁場の方向をベクトルＳｉ，ｃ（ｃ＝１，…，Ｃ）とすると、この点での磁場のベクトルＳｉ，ｃ方向成分Ｂｉ，ｃは、次式のようになる。
2003057035.tif 000005
また、観測磁場として磁場の局所的勾配を測定する場合には、上式のＹｉによる微分を観測することになる。
格子点Ｘｎ（位置ベクトル）での電流双極子の方向に制限をつける場合には、電流双極子の可能な独立方向をｂｎ，ｄ（ｄ＝１，…，Ｄ）として、電流双極子は以下の式で表わされる。この式において、Ｄ＝３の場合は方向に制限がない場合に対応している。
2003057035.tif 000006
以上まとめると、仮想曲面上の格子点｛Ｘｎ｜ｎ＝１，…，Ｎ｝での電流双極子が観測点Ｙｉに作る磁場は次式のように書ける。
2003057035.tif 000007
ここで、ｊｎ，ｄは、格子点Ｘｎでの電流双極子の独立な成分である。また、関数Ｇ（ｉ，ｃ；ｎ，ｄ）は格子点Ｘｎでの電流双極子ｊｎ，ｄが観測点Ｙｉに作る磁場のＳｉ，ｃ方向成分である。
電流分布推定の問題とは、観測された磁場｛Ｂｉ，ｃ｜ｉ＝１，…，Ｉ；ｃ＝１，…，Ｃ｝から仮想曲面上の電流分布｛ｊｎ，ｄ｜ｎ＝１，…，Ｎ；ｄ＝１，…，Ｄ｝を推定する問題であるといえる。
上述した観測点における磁場の測定値を、行列表現を用いて表わすこととすると、以下のように書くことができる。なお、Ｇはリードフィールド行列と呼ばれる。ビオサバールの式ではなく、サルバスの式、あるいは境界要素法等の数値解を用いる場合には、このリードフィールド行列の値が異なるだけであり、以下の説明のその他の部分は同様に適用することができる。
2003057035.tif 000008
（電流源確率モデル）
このように電流モデルを決定した上で、さらにこのような電流分布の推定のために、以下に説明する「電流源の確率モデル」を導入する。
（Ｉ−１−３）  電流源確率モデル
上述した仮想曲面上の電流モデルに対する確率モデルを以下のように考える。観測される磁場は仮想曲面上の電流分布Ｊから作られる磁場と観測ノイズの和として表わされるとする。また、観測ノイズは各測定点で独立な分散σ２を持つガウスノイズとして表わされるとする。
より一般的には、各測定点のノイズの相関が共分散行列の形で表される多次元正規分布に従うガウスノイズを考えることもできるが、以下では説明の簡単のために等方等分散のノイズモデルを用いて説明する。
すなわち、観測磁場が以下の式のように表わされるものとする。
2003057035.tif 000009
この電流モデルに対する確率分布は以下のように与えられる。
まず、ある特定の深さに対する仮想曲面、あるいは各局所面の深さが指定された局所面の集合をＭで表わす。この仮想曲面Ｍ上の電流分布Ｊが与えられたときに、観測される磁場がＢである確率Ｐ（Ｂ｜Ｊ，β，Ｍ）は、以下の式で表わされる。ただし、以下の式において、β＝１／σ２である。
2003057035.tif 000010
（Ｉ−１−４）  階層事前分布
仮想曲面Ｍ上の電流分布Ｊに対する事前分布として、上述したように、局在条件階層事前分布を仮定することにする。
観測を行なう前の電流分布Ｊに対する事前分布（Ｊが実現する確率）は局在条件階層事前分布の仮定の下で以下の式のように表わされる。
2003057035.tif 000011
ただし、Γ（…）は、ガンマ分布を表わし、次式のように定義される。次式において、Γ（γ０）は、ガンマ関数であり、その定義式も併せて示す。
2003057035.tif 000012
また、αとτは、電流分布Ｊとノイズ逆分散βに対する事前分布をモデル化するために導入されたハイパーパラメータであり、このαとτに対する階層事前分布Ｐ０（α｜Ｍ）とＰ０（τ｜Ｍ）は以下のように与えられる。
2003057035.tif 000013
なお、上式において、αとτに対する階層事前分布Ｐ０（α｜Ｍ）とＰ０（τ｜Ｍ）の分布形状を決定するメタパラメータ  バーα０、γα０、τ０、バーγτ０については、後に詳しく説明する。なお、変数名の前に「バー」との語句をつけたものは、その変数の頭部に記号「−」がついていることをあらわす。
（Ｉ−１−５）  ベイズ推定
磁場のデータＢを観測したときに、電流分布がＪある事後確率分布Ｐ（Ｊ｜Ｂ，Ｍ）は、ベイズの定理を使って以下のように求まる。
2003057035.tif 000014
ここで、事後確率分布Ｐ（Ｊ，β，α，τ｜Ｂ，Ｍ）は、以下の式で表される。
2003057035.tif 000015
この事後確率分布を用いて、電流分布の期待値Ｅ［Ｊ｜Ｂ，Ｍ］は以下のように与えられる。
2003057035.tif 000016
また、Ｐ（Ｂ｜Ｍ）は仮想曲面モデルＭの周辺尤度である。電流源の深さを推定するとき、深さの異なるいくつかのモデルを比べる。このとき、深さに関する事前の知識がないものとする。すなわち、以下では、Ｐ（Ｍ）＝一定であるものとする。
観測データＢが与えられたときに、モデルＭが真のモデルである確率、すなわち、モデル事後確率Ｐ（Ｍ｜Ｂ）は、上の仮定のもとでモデル周辺尤度Ｐ（Ｂ｜Ｍ）に比例し、以下の関係が成り立つ。
2003057035.tif 000017
（Ｉ−１−６）  変分法的ベイズ推定
確率モデルと階層事前分布が与えられたときに、モデル周辺尤度を解析的に求めることは一般的にはできない。
そこで、モデル周辺尤度を近似的に計算する手法として、変分法的ベイズ法を用いる。モンテカルロ法や、ラプラス近似法等、他の近似方法を使うことも可能である。
「変分法的ベイズ法」について、具体的な手続きは後述するものとして、以下に簡単に説明しておく。
真の事後分布Ｐ（Ｊ，β，α，τ｜Ｂ，Ｍ）を近似的に求めるために、試験事後分布Ｑ（Ｊ，β，α，τ）を考える。
２つの確率分布Ｐ（Ｊ，β，α，τ｜Ｂ，Ｍ）とＱ（Ｊ，β，α，τ）の近さは、次式で表されるようなＫ−Ｌ距離を用いて計算することができる。
2003057035.tif 000018
Ｋ−Ｌ距離は、２つの分布が等しいときに限りゼロになり、それ以外のときは常に正の値をとる。
上式において試験事後分布Ｑに対する自由エネルギーＦ（Ｑ）を定義すると、以下の不等式が得られることになる。
2003057035.tif 000019
すなわち、自由エネルギーＦ（Ｑ）を最大にする分布Ｑ（Ｊ，β，α，τ）は、真の事後分布Ｐ（Ｊ，β，α，τ｜Ｂ，Ｍ）に等しくなり、このとき自由エネルギーは対数周辺尤度に等しくなる。
変分法的ベイズ法では、Ｑ（Ｊ，β，α，τ）の関数形を以下の形に制限して自由エネルギーの最大化を行なう。
2003057035.tif 000020
上式の右辺第２項Ｑαを固定して、右辺第１項ＱＪに関してＦ（Ｑ）を最大化するステップと、第１項ＱＪを固定して、第２項Ｑαに関してＦ（Ｑ）を最大化するステップを交互に繰返すことにより、自由エネルギーＦ（Ｑ）を極大にする分布Ｑ＊が得られることになる。
［（ＩＩ）電流源推定の手順］
以上のような準備の下に、以下電流源を推定する具体的な手続きについて、図にしたがって説明する。
図４は、本発明の脳内電流源の位置の推定方法の手順の全体的な流れを示すフローチャートである。
図４を参照して、まず、脳内電流源の位置の推定処理が開始されると（ステップＳ１００）、上述したように、電流分布推定のための階層事前分布の分布形状を指定するためのメタパラメータの値と推定すべき変数の初期値を決定する（ステップＳ１０２）。
続いて、初期分解能を用いた電流源の初期推定が行われ、電流源の候補の抽出が行われる（ステップＳ１０４）。
さらに、このようにして抽出された電流源の候補のうち、まず、最も脳表面に近い電流源の位置の特定が行われる（ステップＳ１０６）。その後、他の電流源の深さの特定が順次行われる（ステップＳ１０８）。
以上のようにして、初期分解能による電流源の位置の特定が行われた後、空間分解能を上げて、電流源の位置の再推定が行われ（ステップＳ１１０）、処理が終了する（ステップＳ１１２）。
以下、図４で説明した各ステップの処理について、さらに詳しく説明する。
（ＩＩ−１）  初期分解能を用いた電流源の初期推定（電流源の候補の抽出）
図５は、図４で説明した処理のうち、ステップＳ１０４の処理、すなわち、初期分解能を用いた電流源の初期推定の処理を説明するためのフローチャートである。
図５を参照して、まず、初期分解能での深さ方向のサンプル地点を｛Ｒｋ｜ｋ＝１，…，Ｋ｝とする。各深さＲｋにおける仮想曲面に対して電流分布の推定を行なう。さらに、このようにして推定された電流分布に基づいて、各深さＲｋに対する事後確率を求める。この事後確率が後に説明する自由エネルギー値に相当する（ステップＳ２０２）。
以上のようにして求められた事後確率が最大になる深さＲＭでの電流分布に対して、電流の強さの極大点を求める（ステップＳ２０４）。
さらに、各極大点を含む局所面を決定する。局所面の数をＬ個とすると、各局所面が局在化された電流源の候補となる（ステップＳ２０６）。
以下、図４中のステップＳ１０２の処理および、それに続く図５中のステップＳ２０２の処理について、さらに詳しく説明する。
（ＩＩ−１−１）電流分布推定のための階層事前分布の分布形状を指定するためのメタパラメータの値と推定すべき変数の初期値の決定
上述したとおり、電流源の推定は、変分法的ベイズ法を用いて行われる。したがって、図４のステップＳ１０２の処理は、このような変分法的ベイズ法の適用にあたって、以下のように行われる。
階層事前分布の分布形状を指定するメタパラメータを、たとえば、以下のように決定する。
2003057035.tif 000021
さらに、上述した式に基づいて、各変数の初期設定を以下にように行なう。
2003057035.tif 000022
（ＩＩ−１−２）電流分布推定のための変分法的ベイズ法の具体的処理
（１）  パラメータＪ，βの期待値の計算（Ｊ−ステップの処理）
ここでは、パラメータＪ，βの期待値の計算を行なう。この処理により、ＱＪに関して自由エネルギーＦ（Ｑ）を最大化する処理が行なわれる。
対角行列Ａを以下のように定義することで、ＱＪが以下の式に基づいて導出される。なお、以下の式において、変数の期待値には、その変数名の上に、”−”をつけて表現している。
2003057035.tif 000023
（２）  ハイパーパラメータ期待値、すなわち、パラメータα，τの期待値の計算（α−ステップの処理）
Ｊステップに続いて、ハイパーパラメータα，τの期待値の計算を行なう。すなわち、この処理では、Ｑαに関して自由エネルギーＦ（Ｑ）を最大化する処理を行なうことになる。
この手続は以下の式のように表わされる。
2003057035.tif 000024
（３）  自由エネルギーの計算
以上説明したようなＪステップおよびαステップに基づいて計算されたＱＪおよびＱαを用いて、自由エネルギーを以下のようにして計算する。
2003057035.tif 000025
このようにして、自由エネルギーＦの値が収束するまで、上述したＪステップの処理から自由エネルギーの計算までの処理を繰返す。収束後の自由エネルギーＦ（Ｑ）の値が、対数周辺尤度ｌｏｇＰ（Ｂ｜Ｍ）の近似値を与える。
また、対数モデル事後確率は、ｌｏｇＰ（Ｂ｜Ｍ）と定数分だけの違いしかないので、モデル事後確率が最大になるモデルは、上述したような近似の下では自由エネルギーが最大になるモデルに等しいことになる。
以上のような手続により、深さＲｋに対する事後確率を求めることができる。これに応じて、上述した図５中のステップＳ２０４およびＳ２０６の処理を行なうことにより、各局所面に局在化された電流源の候補を求めることができる。
（ＩＩ−２）  最も脳表面に近い電流源の特定
次に、図４のステップＳ１０６の処理、すなわち、以上のようにして求められた電流源の候補から、まず、最も脳表面に近い電流源の特定を行なう処理について説明する。
図６は、最も脳表面に近い電流源の特定を行なう処理を説明するためのフローチャートである。
まず、初期値として変数ｌの値を１とする（ステップＳ３０２）。つづいて、変数ｌと、変数ｌのとり得る最大値Ｌとの比較が行われ（ステップＳ３０４）、変数ｌが最大値Ｌ以下であれば、処理はステップＳ３０６に移行し、変数ｌが最大値Ｌを超えていれば、処理はステップＳ３２４に移行する（ステップＳ３０４）。すなわち、ステップＳ３０６からＳ３２２までの処理を、ｌ＝１からＬまで繰返す。
ステップＳ３０６では、１番目の局所面以外の局所面の深さを図５のステップＳ２０４で求められた深さＲＭに固定する。
変数ｋの値を１とする（ステップＳ３０８）。つづいて、変数ｋと、変数ｋのとり得る最大値Ｋとの比較が行われ（ステップＳ３１０）、変数ｋが最大値Ｋ以下であれば、処理はステップＳ３１２に移行し、変数ｋが最大値Ｋを超えていれば、処理はステップＳ３２０に移行する。
ステップＳ３１２では、１番目の局所面の深さをまずＲｋとする。
次に、Ｌ個の局所面の集合に対して電流分布の推定を行なう（ステップＳ３１４）。
さらに、この局所面配置に対する事後確率（自由エネルギー）を求める（ステップＳ３１６）。変数ｋの値を１だけインクリメントして（ステップＳ３１８）、処理はステップＳ３１０に復帰する。
このようにして、ｋ＝１からｋ＝Ｋまでの深さの各局所面についての処理を行なった後、事後確率が最大になる１番目の局所面の深さＲＭ（ｌ）を求める（ステップＳ３２０）。変数ｌの値を１だけインクリメントして（ステップＳ３２２）、処理はステップＳ３０４に復帰する。
このようにして、各変数ｌについて求めた局所面の深さＲＭ（ｌ）の中で、最も脳表面に近い（浅い）ｌを求め、これをｌ１とする（ステップＳ３２４）。
以上の処理により、ｌ１番目の局所面が最も脳表面に近い電流源に対応し、その深さの初期推定値がＲＭ（ｌ１）として求められ、処理が終了する（ステップＳ３２６）。
（ＩＩ−３）  各局所面に対応する電流源の深さの特定
図７および図８は、図４のステップＳ１０８の処理、すなわち、各局所面に対応する電流源の深さの特定処理を説明するためのフローチャートである。
図７を参照して、初期値として変数ｓの値を１とする（ステップＳ４０２）。つづいて、変数ｓと、変数ｓのとり得る最大値Ｌとの比較が行われ（ステップＳ４０４）、変数ｓが最大値Ｌ未満であれば、処理はステップＳ４０６に移行し、変数ｓが最大値Ｌ以上であれば、処理はステップＳ４３４に移行する（ステップＳ４０４）。すなわち、ステップＳ４０６からＳ４３２までの処理を、ｓ＝１からＬまで繰返す。
ステップＳ４０６では、これまでに特定された局所面の深さを｛ＲＭ（ｌ１），…，ＲＭ（ｌｓ）｝とするとき、これらの局所面の深さを、それぞれ｛Ｒ，Ｍ（ｌ１），…，ＲＭ，（ｌｓ）｝に固定する（ステップＳ４０６）。
まず、ｌが｛ｌ１，…，ｌｓ｝のいずれとも等しくなく、かつｌは｛１，…，Ｌ｝の集合に属しているものとする（ステップＳ４０８）。次に、変数ｌの値として、可能な値を全て尽くしたか否かを判断し（ステップＳ４１０）、全て尽くしている場合は処理はステップＳ４３０に移行する。一方、全てを尽くしていない場合は、処理はステップＳ４１２に移行する。
ステップＳ４１２では、１番目の局所面以外で、かつ｛ｌ１，…，ｌｓ｝でもない局所面の深さをＲＭに固定する。
変数ｋの値を１とする（ステップＳ４１４）。つづいて、変数ｋと、変数ｋのとり得る最大値Ｋとの比較が行われ（ステップＳ４１６）、変数ｋが最大値Ｋ以下であれば、処理はステップＳ４１８に移行し、変数ｋが最大値Ｋを超えていれば、処理はステップＳ４２６に移行する。すなわち、ステップＳ４１８からＳ４２４までの手続を、ｋ＝１からｋ＝Ｋまで繰返す。
ステップＳ４１８では、１番目の局所面の深さをＲｋにする。
次に、Ｌ個の局所面の集合に対して電流分布推定を行なう（ステップＳ４２０）。
さらに、この局所面配置に対する事後確率（自由エネルギー）を求める（ステップＳ４２２）。
図８を参照して、ステップＳ４２６では、以上のようにしてｋ＝Ｋまでの処理が終わった後に、事後確率が最大になる１番目の局所面の深さＲＭ（ｌ）を求める（ステップＳ４２６）。
｛ｌ１，…，ｌｓ｝のいずれとも等しくなく、かつ｛１，…，Ｌ｝の集合に属している、まだ処理を行っていない他の変数ｌを選択し（ステップＳ４２８）、ステップＳ４１０に復帰する。
以上の手続により、処理を行った全ての変数ｌに対応するＲＭ（ｌ）の中で、脳表面に最も近いｌを求め、これをｌｓ＋１とする（ステップＳ４３０）。
さらに、ｓにｓ＋１を代入し（ステップＳ４３２）、処理はステップＳ４０４に復帰する。
全ての可能なｓについて処理が終了することで、各局所面に対応する電流源の深さの特定が終了する（ステップＳ４３４）。
（ＩＩ−４）  分解能を上げて電流源の再推定
図９および図１０は、図４のステップＳ１１０の処理、すなわち、空間分解能を上げて電流源の位置を再推定する処理を説明するためのフローチャートである。
図９および図１０を参照して、まず、初期分解能を用いて推定した電流源と対応する局面の番号を脳表面に近いものから１，２，３，…，Ｌとなるように番号の並べ替えを行なう（ステップＳ５０２）。次に、各局所面での電流分布の拡がりに対応して局所面を小さくする（ステップＳ５０４）。
図４のステップＳ１０８までの処理で求めた局所面の深さＲＭ（ｌ）をＲＭｏｌｄ（ｌ）とする（ステップＳ５０６）。
新しい深さ方向の分解能と探索幅を、それぞれΔＲおよび（ＫＬ・ΔＲ）とする。また局所面上での格子点の分解能も上げる（ステップＳ５０８）。
初期値として変数ｌの値を１とする（ステップＳ５１０）。つづいて、変数ｌと、変数ｌのとり得る最大値Ｌとの比較が行われ（ステップＳ５１２）、変数ｌが最大値Ｌ以下であれば、処理はステップＳ５１４に移行し、変数ｌが最大値Ｌを超えていれば、処理はステップＳ５３４に移行する（ステップＳ５１２）。すなわち、ステップＳ５１４からＳ５３２までの処理を、ｌ＝１からｌ＝Ｌまで繰返す。
ステップＳ５１４では、１番目以外の局所面ｌ′の深さをＲＭｏｌｄ（ｌ′）に固定する。
その上で、さらに、ｋ＝０，±１，…，±ＫＬについて、以下のステップＳ５１８からステップＳ５２６までの処理を行なう。
まず、変数ｋの値を０とする（ステップＳ５１６）。つづいて、変数ｋの絶対値と、変数ｋの絶対値のとり得る最大値ＫＬとの比較が行われ（ステップＳ５１８）、変数ｋの絶対値が最大値ＫＬ以下であれば、処理はステップＳ５２０に移行し、変数ｋの絶対値が最大値ＫＬを超えていれば、処理はステップＳ５２８に移行する（ステップＳ５１８）。
ステップＳ５２０では、１番目の局所面の深さを（ＲＭｏｌｄ（ｌ）＋ｋ・ΔＲ）とする。
次に、Ｌ個の局所面の集合に対して、現在の分解能を用いて、電流分布推定を行なう（ステップＳ５２２）。
さらに、この局所面配置に対する事後確率（自由エネルギー）を求める（ステップＳ５２４）。その上で、ｋの値を｛±１，…，±ＫＬ｝のうちの次の値とし、処理はステップＳ５１８に復帰する。
以上のような処理がｋ＝±ＫＬまで行なわれた後に、ステップＳ５２８では、事後確率が最大となるｋを求め、これをｋＭとする。
その上で、ＲＭｏｌｄ（ｌ）にＲＭｏｌｄ（ｌ）＋ｋＭ・ΔＲを代入する（ステップＳ５３０）。変数ｌの値を１だけインクリメントして（ステップＳ５３２）、処理はステップＳ５１２に復帰する。
以上のような処理が変数ｌの値が最大値Ｌとなるまで行われた時点で、分解能が最終分解能に達していれば（ステップＳ５３４）、処理を終了する（ステップＳ５３８）。
一方、分解能が最終分解能に達していなければ（ステップＳ５３４）、局所面上での格子点の分解能と深さ方向の分解能を上げて（ステップＳ５３６）、処理は、ステップＳ５１０に復帰する。
以上説明したような「脳内電流源推定方法」を用いれば、ＭＥＧ（またはＥＥＧ）の観測データから脳内電流源の位置を深さ方向まで含めて推定することが可能となる。さらに、このような深さ方向の推定は、複数の電流源が存在する場合にも適用可能である。しかも、電流源が電流双極子のように局在している場合にも、広がりを持つ電流源の場合にも適用できるだけでなく、電流源の広がり具合を推定することも可能である。
また、図９および図１０で説明したような処理を、初期分解能での推定の後に追加して行えば、位置の推定分解能を逐次的に上げることが可能となる。このことは、また、生理学的な知見や他の観測方法による観測データなどにより探索範囲を限定して調べることが可能なことも意味する。
［シミュレーション結果］
以下では、以上説明した脳内電流源推定方法に従って、モデル的に仮定した電流源位置を推定するシミュレーションを行なった結果について説明する。
図１１は、人間の脳を、半径１０．０（任意単位）の半球状であると仮定し、この半球の表面上で観測された動径方向の磁場分布を上面から見た上面図である。
図１１においては、後に説明するように、中心から半径７．０の位置に、単一の電流源が置かれている場合の半球表面上での磁場分布を示している。以下のシミュレーションでは磁場の観測データにＳ／Ｎ比が０．１のノイズを加えている。
図１２は、初期分解能を用いた電流源の初期推定結果を示す図である。
図１２では、図４のステップＳ１０２〜ステップＳ１０４の処理において説明したとおり、初期分解能を用いて電流源の初期推定を行なうために、深さ（半径Ｒ）をステップ状に変化させつつ、各深さの半球上において自由エネルギーが最大となるような電流分布を変分法的ベイズ法に基づいて算出した結果を示す。
図１２（ａ）は半径Ｒ５．０の場合を示し、図１２（ｂ）は半径Ｒ＝６．０の場合を示し、図１２（ｃ）は半径Ｒ＝７．０の場合を示す。
図１３は、さらに、半径がより大きい（脳表面により近く）場合の初期分解能を用いた電流源の初期推定結果を示す図である。
図１３（ｄ）は半径Ｒ＝７．５の場合であり、図１３（ｅ）は半径Ｒ＝８．０の場合であり、図１３（ｆ）は半径Ｒ＝９．０の場合である。
図１４は、このように、脳内に想定した様々な深さの仮想半球面に対して、自由エネルギーを最大とするような電流分布を各々求める際に得られる、各仮想半球面上での自由エネルギーの大きさを示す図である。
図１４に示すとおり、半球面全体を仮想曲面として計算した場合、半径Ｒ＝７から半径Ｒ＝８の間に自由エネルギーの最大点が存在していることがわかる。
図１３（ｄ）に示すとおり、半径Ｒ＝７．５の仮想半球面上には、電流分布に１つの極大点が存在する。
このようにして、初期分解能を用いて電流源が初期推定されたので、自由エネルギーが最大となっている仮想曲面において、電流分布の極大点について、その極大点を含む局所面を求め、さらにこの局所面について電流源の特定処理が行なわれる。
図１５および図１６は、このように、電流分布の極大点を含む局所面について、局所面上の格子点の分解能と深さ方向の分解能を上げて、さらに電流源の特定処理を行なった場合の処理を示す図である。
図１５（ａ）は半径Ｒ＝５．０の、図１５（ｂ）は半径Ｒ＝６．０の、図１５（ｃ）は半径Ｒ＝７．０の場合の局所面における電流分布を示している。同様に、図１６（ｄ）は、半径Ｒ＝７．５の、図１６（ｅ）は半径Ｒ＝８．０の、図１６（ｆ）は半径Ｒ＝９．０の電流分布をそれぞれ示している。
図１７は、このようにして局所面の深さを変えつつ自由エネルギーを求めた結果を示す図である。
図１７に従えば、局所面を用いて計算した結果の自由エネルギーは、半径Ｒ＝７付近に極大値を持ち、結果から、半径Ｒ＝７の位置に１つの電流源が存在していることを特定することが可能である。
このときの電流源における電流分布は、図１５（ｃ）のとおりであり、上述したとおり、電流分布も最も拡がりの狭いものとなっていることがわかる。
図１８は、さらに、脳表面を半球と仮定した場合に、電流源が２つ脳内に存在する場合の脳表面における磁場分布を示す上面図である。
図１８に示した例においては半径Ｒ＝６．０およびＲ＝８．０の位置にそれぞれ電流源が置かれているものとする。
図１９および図２０は、初期分解能を用いた電流源の初期推定結果を示す図である。
すなわち、図１９および図２０では、２個の電流源が脳内に存在している場合に、図４に示した初期分解能を用いた電流源の初期推定を行なう処理（ステップＳ１０４）を説明している。
図１９（ａ）は半径Ｒ＝５．０の場合の、図１９（ｂ）は半径Ｒ＝６．０の場合の、図１９（ｃ）は半径Ｒ＝７．０の場合の仮想半球面上の電流分布を示している。
図２０（ｄ）は、半径Ｒ＝７．５の場合の、図２０（ｅ）は半径Ｒ＝８．０の場合の、図２０（ｆ）は半径Ｒ＝９．０の場合の仮想半球面上の電流分布をそれぞれ示している。
図２１は、このようにして仮想半球面全体について、各仮想半球面ごとに自由エネルギーが最大となるように電流分布を求めた場合の、自由エネルギーの半径依存性を示す図である。
図２１に示すとおり、半球面全体にわたって自由エネルギーを計算した場合、半径Ｒ＝７から半径Ｒ＝８の間の半径７．５付近に自由エネルギーの最大値が存在していることがわかる。
図２０（ｄ）に示すとおり、半径Ｒ＝７．５の仮想半球面上には、電流分布に２つの極大点が存在し、２つの局所面が決定される。
次に、図４におけるステップＳ１０６の処理に対応する計算結果について説明する。
図２２および図２３は、このようにして、最も脳表面に近い電流源を特定するために、事後確率が最大になる１番目の局所面の深さＲＭ（ｌ）を求める処理を説明するための図である。
したがって、図２２および図２３は、半球面全体にわたって計算した自由エネルギーが最大となる面上に他方の局所面を固定したままで、一方の電流源に対応する局所面を移動させた場合の電流分布を示す図である。
図２２（ａ）は半径Ｒ＝５．０の場合の、図２２（ｂ）は半径Ｒ＝６．０の場合の、図２２（ｃ）は半径Ｒ＝７．０の場合の局所面上の電流分布を示している。
図２３（ｄ）は、半径Ｒ＝７．５の場合の、図２３（ｅ）は半径Ｒ＝８．０の場合の、図２３（ｆ）は半径Ｒ＝９．０の場合の局所面上の電流分布をそれぞれ示している。
図２４は、図２２および図２３で示したように仮想局所面を移動させた場合の自由エネルギーの局所面位置（半径）依存性を示す図である。
図２４に示すとおり、半球面全体にわたって計算した自由エネルギーが最大となる面上に他方の局所面を固定したままで、一方の電流源に対応する局所面を半径Ｒ＝８の位置にしたときに、自由エネルギーが最大となる。
したがって、図２４に示した結果から、最も脳表面に近い電流源の位置としては、半径Ｒ＝８．０の位置であることがわかる。
続いて、他方の側の電流源の深さの特定を行なう。
この場合、一方の電流源に対応する局所面の深さは、半径Ｒ＝８．０に固定しておき、他方の電流源に対応する局所面の深さを変化させていく。
図２５および図２６は、特定された面上に一方の局所面を固定したままで、他方の電流源に対応する局所面を移動させた場合の電流分布を示す図である。
図２５（ａ）は半径Ｒ＝５．０の場合の、図２５（ｂ）は半径Ｒ＝５．５の場合の、図２５（ｃ）は半径Ｒ＝６．０の場合の局所面上の電流分布を示している。
図２６（ｄ）は、半径Ｒ＝６．５の場合の、図２６（ｅ）は半径Ｒ＝７．０の場合の、図２６（ｆ）は半径Ｒ＝７．５の場合の局所面上の電流分布をそれぞれ示している。
図２７は、図２５および図２６で示したように仮想局所面を移動させた場合の自由エネルギーの局所面位置（半径）依存性を示す図である。
図２７に示すとおり、この場合半径Ｒ＝６の位置で自由エネルギーが極大となることがわかる。
このような処理により、２つの電流源が脳内に存在している場合も、各々の電流源の深さを特定することが可能であることがわかる。
［実施の形態２］
以下では、上述した実施の形態１における「脳内電流源推定方法」、「脳内電流源推定プログラム」および「脳内電流源推定装置」の構成において、ソフトウェア（または、一部ハードウェアにより処理されてもよい）により実現される処理のうち、「（Ｉ）電流源推定のための準備」から「（３）自由エネルギーの計算」までの処理を変更した本発明の実施の形態２について説明する。
まず、実施の形態２における変更点を説明する前提として、実施の形態１の処理の概念を改めて以下にまとめるとともに、その基本概念のうち、実施の形態２において変更される概念を説明しておく。
（実施の形態１における仮想曲面推定の概念）
実施の形態１においては、仮想曲面上の電流モデルとして、格子点における電流双極子を想定し、観測面上の観測磁場Ｂから電流分布Ｊを推定するという逆問題を解くという処理を行った。このとき、このような逆問題は、一般に、（電流双極子の存在する格子点数）が（観測点数）よりも多いために、いわゆる「不良設定問題」となり、厳密に解を求めることが困難である。
実施の形態１では、このような問題を解くために、「ベイズ推定」の手法を用いた。
ベイズ推定においては、仮想曲面Ｍと観測磁場Ｂが与えられた時の事後確率分布が以下のように表現されることを利用していた。
2003057035.tif 000026
ここで、Ｐ（Ｊ｜Ｂ，Ｍ）は、磁場Ｂが観測されたと言う条件での、仮想曲面Ｍ上での電流の確率分布であって、「事後確率分布」を表す。また、Ｐ（Ｂ｜Ｊ，Ｍ）は「データ尤度」であり、Ｐ（Ｊ｜Ｍ）は事前分布であり、Ｐ（Ｂ｜Ｍ）は周辺尤度である。
周辺尤度は、以下の式で表される。
2003057035.tif 000027
事後確率Ｐ（Ｍ｜Ｂ）は、以下の式で表される。
2003057035.tif 000028
深さに関する事前の知識がないものとすると、Ｐ（Ｍ）＝一定となるので、実施の形態１でも説明したとおり、観測データＢが与えられたときに、モデルＭが真のモデルである確率、すなわち、モデル事後確率Ｐ（Ｍ｜Ｂ）は、上の仮定のもとでモデル周辺尤度Ｐ（Ｂ｜Ｍ）に比例する。
すなわち、実施の形態１において、電流モデルの選択にあたっては、事後確率、すなわち、周辺尤度が最大となるモデルを選ぶという処理を行った。
このような処理を行うために、以下の式で表される対数周辺尤度ｌｏｇＰ（Ｂ｜Ｍ）を計算の対象とする。
2003057035.tif 000029
上式に示すとおり、対数周辺尤度を最大にするということは、概念的には、誤差と実効的なパラメータ自由度とを最小にすること、言いかえると、誤差と電流分布の広がりを最小とすることに相当する。
実施の形態１においては、上記のような対数周辺尤度の最大化の問題を、試験事後分布Ｑを事後分布の近似として導入し、自由エネルギーＦ（Ｑ）を最大化するように事後分布を決定するという、いわゆる「変分法的ベイズ法」を用いて解いていた。
（実施の形態１における事前分布の概念）
実施の形態１においては、脳内電流源に対する事前知識として、電流源は脳内の複数の場所に局在しているという局在条件を用い、これを階層的事前分布の形で表わした。すなわち、ハイパーパラメータαｎ（ｎ＝１，…，Ｎ）を導入して、次式の階層事前分布の形を仮定した。
2003057035.tif 000030
上式からハイパーパラメータαｎは電流Ｊｎの逆分散（分散の逆数）に比例していることがわかる。もし、電流Ｊｎの分散とノイズの逆分散βの値が予めわかっている場合には、αｎの値を電流の逆分散とノイズの逆分散の値とから計算される値に固定した事前分布を用いればよい。しかし、電流の分散の値は通常わからないので、実施の形態１においては、ハイパーパラメータαｎに対する階層事前分布Ｐ０（αｎ｜バーα０，γα０）を仮定してαｎに対するベイズ推定を行なっていた。
ここで、バーα０とγα０は、階層事前分布の分布形状を決めるメタパラメータであり、すべてのαｎに対して共通の値を仮定していた。この局所条件階層事前分布を用いてベイズ推定を行なった場合、前述したようにできるだけ少ない数の電流双極子で観測データを復元するようなベイズ推定が行なわれる。
一方、ＭＥＧの電流源推定においてしばしば用いられる最小ノルム推定法（文献：Ｈａｍａｌａｉｎｅｎ  Ｍ．Ｓ．ａｎｄ  Ｉｌｍｏｎｉｅｍ  Ｒ．Ｊ．，Ｍｅｄ．＆  Ｂｉｏｌ．Ｅｎｇ．＆  Ｃｏｍｐｕｔ．，３２，３５−４２，１９９４）は、ベイズ推定の立場からは、事前分布としてαｎの値をすべての点において共通の値にし、しかもこの値を直接適当な値に定めることに相当する。
最小ノルム推定では、すべての点における電流強度の総和を最小にしようとするので、推定解には小さな電流強度を持つ点が多く現れる。このため、これらの小さな電流分布が真の電流分布を与えているのか、それとも推定誤差のためのノイズとみなしてよいのかを判定することが難しいという欠点を持っている。
局在条件階層事前分布は、最小ノルム推定法のこのような欠点を解消してくれる。
（実施の形態２における変更点）
実施の形態２においては局所条件のみならず連続条件も考慮に入れた事前分布を仮定する。空間分解能がミリメートル単位の高精度の電流源分布推定を行なう場合には、電流分布の連続性を考慮する必要がある。実際、大脳皮質の神経細胞は半径１０ｍｍ程度のコラム構造を形成しており、神経活動が起こる場合には半径１０ｍｍ程度の領域内における神経細胞が同時に発火する確率が高い。このような点を考慮して、以下に説明するように実施の形態２においては、局所条件と連続条件を組合せた階層事前分布を仮定する。
連続条件を導入するために、新たに内部変数Ｚを導入し、電流Ｊは内部変数Ｚを平滑化フィルタＷで平滑化したものとして表現する。ここで、このような平滑化フィルタとしては、特に限定されないが、たとえばガウスフィルタ等を用いることができる。このとき、事前分布は以下のように表わせる。
2003057035.tif 000031
上記の事前分布においてＺに対する積分は容易に実行できて、以下のように書ける。
2003057035.tif 000032
但し、ＡとΛは、それぞれ｛αｎ｜ｎ＝１，…，Ｎ｝と｛λｎ｜ｎ＝１，…，Ｎ｝を対角成分に持つ対角行列である。
すなわち、内部変数Ｚを導入して事前分布Ｐ０（Ｊ｜Ｚ，α，β）Ｐ０（Ｚ｜λ）を仮定することは、電流Ｊに対する共分散行列β−１（Ａ−１＋Ｗ・Λ−１・Ｗ′）を持つ事前分布を仮定することに等しい。但し、内部変数Ｚを用いて事前分布を表わした方が変分法的ベイズ法を適用しやすいので、以下では内部変数Ｚを用いて事前分布を表わす。また、ハイパーパラメータαｎと同様に、ハイパーパラメータλｎの値も事前にはわからないので、以下の形の階層事前分布を仮定して、αｎとλｎに対するベイズ推定を行なう。
2003057035.tif 000033
このように拡張された階層的事前分布は、その特別の場合として実施の形態１の局在条件事前分布を含む。すなわち、平滑化フィルタ行列Ｗが恒等的に０の場合、上記の連続条件と局在条件を組合せた階層的事前分布は、実施の形態１の局在条件事前分布に一致する。
また、階層的事前分布の分布形状を決めるメタパラメータ  バーα０ｎ，γα０ｎ，バーλ０ｎ，γλ０ｎは、より一般的に格子点の場所に依存してもよいと仮定している。ＭＥＧ（あるいはＥＥＧ）の観測データしかない場合は、これらの値に対する事前の知識は特にないので、後述するように、実施の形態１と同様に、実質的にはこれらのメタパラメータを格子点の場所に依存せずに共通に決める。しかし、実施の形態３で説明するように、他の観測手段による観測データが利用できる場合には、それらの情報を用いてメタパラメータの値を格子点毎に決めることができる。
（実施の形態２の具体的な計算手続き）
以下で説明する計算手続きは、その大筋において、実施の形態１と同様である。したがって、以下では、実施の形態１からの変更点を主として説明することにする。
まず、実施の形態１と同様にして、仮想曲面上の電流分布Ｊが与えられたときに、観測される磁場がＢであるモデル確率分布Ｐ（Ｂ｜Ｊ，β）は、対数表現を用いると、以下の式で表わされる。
2003057035.tif 000034
［階層事前分布］
実施の形態２の階層事前分布は、以下の式で表される。
2003057035.tif 000035
ここで、βは実施の形態１と同様に、ノイズの逆分散を示す。
［試験事後分布の計算］
このとき、試験事後分布は、Ｑ＝ＱＪ（Ｊ，β）ＱＺ（Ｚ）Ｑα（α，λ，τ）のように、ＱＪとＱＺとＱαとの積の形に表されると仮定する。
そして、自由エネルギーＦ（Ｑ）を順番にＱＪとＱＺとＱαの各々に関して最大化する。まず、第１ステップで、ＱＺとＱαとを固定して、ＱＪに関してＦ（Ｑ）を最大化すると、以下の式が得られる。
2003057035.tif 000036
次に、第２ステップで、ＱＪとＱαとを固定して、ＱＺに関してＦ（Ｑ）を最大化すると、以下の式が得られる。
2003057035.tif 000037
さらに、第３ステップで、ＱＪとＱＺとを固定して、Ｑαに関してＦ（Ｑ）を最大化すると、以下の式が得られる。
2003057035.tif 000038
なお、上式における未知変数の具体的計算方法については、以下に説明する。
実施の形態１と同様にして、実施の形態２においても、繰り返し計算を用いて自由エネルギーＦ（Ｑ）を極大化する試験事後分布を求める。
具体的には、以下に説明する手続きを用いて、上述した試験事後分布Ｑを以下の第１〜第３のステップで計算し、自由エネルギーＦ（Ｑ）の値が収束するまでこの手続きを繰り返す。なお、その際の定数項の決め方については、後に説明する。
（１）第１のステップ（Ｊ−ステップと呼ぶ）の手続き
電流Ｊと逆分散βの期待値および変数γβについて以下の式に従い計算する。
2003057035.tif 000039
（２）第２のステップ（Ｚ−ステップ）の手続き
続いて、内部変数Ｚの期待値は以下のとおり計算される。
2003057035.tif 000040
（３）第３のステップ（α−ステップ）の手続き
続いて、ハイパーパラメータαｎ、γαｎ、λｎ，γλｎ、γτについて、以下の計算を行う。
2003057035.tif 000041
以上の（１）〜（３）の手続きにより、試験事後分布Ｑの計算を行うことができる。
［自由エネルギーの計算］
実施の形態２の自由エネルギーは以下の式で計算される。
2003057035.tif 000042
期待対数尤度ＬＰおよびモデル複雑度項ＨＪ，ＨＺ，Ｈα，Ｈλ，Ｈβ，Ｈτは、以下のように計算される。
2003057035.tif 000043
実施の形態２でも、事後分布Ｐ（Ｊ｜Ｂ）は、Ｆ（Ｑ）を最大化するＱＪに対応し、対数周辺尤度ｌｏｇ（Ｐ（Ｂ））が最大化されたＦ（Ｑ）に対応する。
以上の手続きにより、ある深さの仮想曲面についての対数周辺尤度ｌｏｇ（Ｐ（Ｂ））を求めることができる。
他の手続きは、実施の形態１と同様であるので、その説明は繰り返さない。
［階層事前分布のメタパラメータ］
階層事前分布の計算において、メタパラメータの選択について説明する。
他の観測方法による観測データが利用できないときには、実施の形態１と同様にメタパラメータの値は格子点によらずに共通に選ぶ。
すなわち、バーτ０，γβ０，γτ０は実施の形態１と同様に選ぶことができる。また、メタパラメータ  バーα０ｎ，γα０ｎ，バーλ０ｎ，γλ０ｎは、以下のように表わせる。
2003057035.tif 000044
ここで、バーα０とγα０は実施の形態１と同様に選ぶことができる。
（実施の形態３）
実施の形態１および２においては他の観測方法による観測データが利用できないという前提に立って、脳内電流源の推定方法を説明してきた。
しかし、近年、脳の神経活動を観測するためのさまざまな手法が開発されており、同じ実験条件の下で複数の観測手段を用いた観測データが得られることも可能になってきた。そこで、実施の形態３においては、ＭＲＩやｆＭＲＩやＰＥＴ等の他の観測方法による観測データが利用できる場合において、実施の形態１および２に対する変更点を説明する。
まず、同一人物に対するＭＲＩ観測データが得られる場合には、ＭＲＩ画像から大脳皮質の位置やその他の脳部位の位置に関する情報が得られる。このような場合、実施の形態１および２において仮定された仮想曲面の形を大脳皮質の形や他の脳部位の形に合せることができる。
また、たとえば、ある実験条件下での神経活動が大脳皮質で起こっていることが確かなときには、深さ方向の探索を省略し、大脳皮質上にのみ格子点を置いて電流分布推定を行なうことができる。
ｆＭＲＩやＰＥＴでは、神経活動の場所に関する情報を得ることができる。しかし、これらの情報の取扱には注意が必要である。なぜなら、ｆＭＲＩでは血流量を、ＰＥＴでは注入した放射性同位元素の量を測っており、これらは神経活動が高いところほど多くなると考えられるものの、神経活動による電流強度とは測定しているものが異なるからである。
ｆＭＲＩやＰＥＴでの活動度と神経活動による電流強度はある程度の相関はあると考えられるが、その相関が線形的なものか、あるいは非線形な関係になっているかは明らかでない。また、ＭＥＧやＥＥＧはミリ秒単位の時間分解能を持っているが、ｆＭＲＩは秒単位の時間分解能であり、ＰＥＴは数十秒単位の時間分解能しかない。すなわち、ｆＭＲＩやＰＥＴで得られるのはＭＥＧやＥＥＧでの測定結果を長時間平均したものである。
Ｄａｌｅ等は、文献：Ｄａｌｅ  Ａ．Ｍ．ｅｔ．Ａｌ．Ｎｅｕｒｏｎ．２６，ｐｐ．５５−６７（２０００）において、ＭＲＩやｆＭＲＩやＰＥＴの観測結果をＭＥＧやＥＥＧの電流源推定に利用する方法を提案している。しかし、その方法は、ｆＭＲＩ等の活動情報を電流の分散情報として直接与えるものである。
これは、本発明におけるベイズ推定の立場からは、逆分散ハイパーパラメータαｎの値を直接指定することに対応する。しかし、上述したようにｆＭＲＩ等で得られる情報は電流そのものに対するものではなく、また時間的にも、ＭＥＧやＥＥＧで測っている活動を長時間平均したものである。
これに対して、本発明においては、ｆＭＲＩ等から得られる情報を電流の分散情報として直接指定するのではなく、階層事前分布のメタパラメータのレベルで指定する。メタパラメータγα０ｎとγλ０ｎは、階層事前分布により与える情報の信頼度を制御するメタパラメータである。すなわち、γα０ｎやγλ０ｎの値が小さいほど信頼性が低いことを表わし、γα０ｎやγλ０ｎの値が大きくなるほど信頼性が高いことを表わしている。
また、実施の形態２で説明したように、局在条件と連続条件を組合せた階層事前分布は、事前分布における電流の共分散行列としてβ−１（Ａ−１＋Ｗ・Λ−１・Ｗ′）を仮定していることに等しい。但し、ＡとΛはそれぞれ｛αｎ｜ｎ＝１，…，Ｎ｝と｛λｎ｜ｎ＝１，…，Ｎ｝を対角成分に持つ対角行列である。
すなわち、αｎの値が小さいほど、電流の分散は大きくなり、λｎの値が小さいほど電流の連続条件による共分散成分が大きくなる。メタパラメータ  バーα０ｎとバーλ０ｎは、ＭＥＧ（あるいはＥＥＧ）による観測データがないときのαｎとλｎの期待値を表わしている。
ｆＭＲＩや他の観測方法により電流格子点ｎにおける神経活動が高いことが予想される場合は、メタパラメータ  バーαＵｎとバーλＤｎの値を小さくして、この電流格子点における電流が大きい値をとりやすいようにバイアス情報を入れることができる。
このように、階層事前分布のメタパラメータのレベルで測っている量や時間スケールの異なる他の観測方法による情報を入れることで、曖昧な情報として他の観測手段の情報を取り入れることができる。
（シミュレーション結果）
以下、実施の形態２および３の手続によるシミュレーション結果について示す。
以下では、ある人物のＭＲＩ画像から得られた大脳皮質の後頭部にある視覚野を含む領域を切り出したものを仮想曲面とし、深さ方向の探索は省略している。
シミュレーションデータとして、視覚野のＶ１，Ｖ２，Ｖ５の領域にそれぞれ電流強度が（１．０，１．０，０．５）（単位は任意）の電流を流した場合と、（１．０，０．０，０．５）の電流を流した場合について、これらの電流がＭＥＧセンサ上に作る磁場を計算し、これにＳ／Ｎ比が０．１のガウスノイズを加えて観測データとした。
図２８および図２９は、局在条件のみを用いてベイズ推定した結果を示す。図中の各三角形の頂点は、この仮想曲面上に用意した格子点を表わし、この各々の格子点上に電流双極子が割当てられている。各格子点の間隔は平均約３ｍｍであり空間分解能は高い。また、この仮想曲面は大脳皮質の一部を切り出したもので複雑な三次元形状をしている。図は、この三次元形状を二次元投影したものである。図中で、Ｖ１，Ｖ２，Ｖ５と示されているのは、それぞれ視覚野のＶ１，Ｖ２，Ｖ５領域を示している。図中で太い黒丸で示されている点は、推定電流強度が高い点を表わす。
図３０および図３１は、実施の形態２で説明した局在条件と連続条件を組合せた階層事前分布を用いてベイズ推定を行なった結果を示す。この推定結果は、ほぼ真の値に近い電流分布の位置と広がりを再現している。一方、図２８と図２９に示されているように、局在条件を用いてベイズ推定を行なった場合、電流の位置は正しく推定されているが、電流の広がりが過小評価されていることがわかる。このシミュレーションから空間分解能が高い場合に、連続条件を入れる有効性がわかる。
次に、電流分布はそのままにして、ノイズのＳ／Ｎ比を０．２にしてシミュレーションを行なった。この場合、局在条件と連続条件を組合せた階層事前分布を用いても、強いノイズの影響で電流分布の広がりが正確には再現されなかった。
続けて、このような劣悪の条件下で他の観測手段による観測データの情報が利用できる場合についてシミュレーションを行なった。シミュレーションでは、ｆＭＲＩによる観測データからＶ１とＶ２とＶ５の領域の活動が高いという情報が得られたと仮定した。ｆＭＲＩの活動領域と電流活動が必ずしも一致せず、またｆＭＲＩはＭＥＧ観測データの長時間平均に対応していることを考慮して、ｆＭＲＩで高い活動を示している領域は、Ｖ１とＶ２とＶ５の各々の領域を半径２０ｍｍ程度広げた近傍領域であるとした。また、Ｖ１，Ｖ２，Ｖ５の電流強度がそれぞれ（１．０，１．０，０．５）と（１．０，０．０，０．５）の両方の場合に同じｆＭＲＩ情報を用いた。
以上の設定により、ｆＭＲＩの活動領域と電流強度は必ずしも一致しないが、ある程度の相関はあるという状況をシミュレーションすることになる。このような設定の下で、ｆＭＲＩの活性化領域に対応する格子点のメタパラメータをバーα０ｎ＝バーλ０ｎ＝２×Ｔｒ（Ｇ′・Ｇ）／ティルダＮとし、そうでない格子点のメタパラメータをバーα０ｎ＝バーλ０ｎ＝５００００×Ｔｒ（Ｇ′・Ｇ）／ティルダＮとした。なお、変数の前に「ティルダ」との語句をつけた場合、その変数の頭部に「〜」がつけえられていることを示す。
また  γα０ｎおよびγλ０ｎは全点共通にγα０ｎ＝γλ０ｎ＝０．１とした。
このときの推定結果を示したのが、図３２と図３３である。この図から、ノイズが非常に大きい場合でも他の観測手段の情報を加えることにより、正確な推定を行なえることがわかる。
この発明を詳細に説明し示してきたが、これは例示のためのみであって、限定となってはならず、発明の精神と範囲は添付の請求の範囲によってのみ限定されることが明らかに理解されるであろう。
【図面の簡単な説明】
図１は、脳磁計システムの構成の一例を示す概念図である。
図２は、電流源が作る電磁場を適当な曲面上で観測している場合の概念図である。
図３は、脳内の電流源と複数の仮想曲面との関係を示す概念図である。
図４は、本発明の脳内電流源推定方法の手順の全体的な流れを示すフローチャートである。
図５は、初期分解能を用いた電流源の初期推定の処理を説明するためのフローチャートである。
図６は、最も脳表面に近い電流源の特定を行なう処理を説明するためのフローチャートである。
図７は、各局所面に対応する電流源の深さの特定処理を説明するための第１のフローチャートである。
図８は、各局所面に対応する電流源の深さの特定処理を説明するための第２のフローチャートである。
図９は、空間分解能を上げて電流源の位置を再推定する処理を説明するための第１のフローチャートである。
図１０は、空間分解能を上げて電流源の位置を再推定する処理を説明するための第２のフローチャートである。
図１１は、人間の脳を半球状であると仮定し、この半球の表面上で観測された磁場分布を上面から見た上面図である。
図１２（ａ）、図１２（ｂ）および図１２（ｃ）は、それぞれ、半径Ｒ５．０、半径Ｒ＝６．０および半径Ｒ＝７．０の場合において、初期分解能を用いた電流源の初期推定結果を示す図である。
図１３（ｄ）、図１３（ｅ）および図１３（ｆ）は、それぞれ、半径Ｒ＝７．５、半径Ｒ＝８．０および半径Ｒ＝９．０の場合であって、脳表面により近い場合の初期分解能を用いた電流源の初期推定結果を示す図である。
図１４は、自由エネルギーを最大とするような電流分布を各々求める際に得られる、各仮想半球面上での自由エネルギーの大きさを示す図である。
図１５（ａ）は半径Ｒ＝５．０、図１５（ｂ）は半径Ｒ＝６．０、図１５（ｃ）は半径Ｒ＝７．０の場合において、極大点を含む局所面についてさらに電流源の特定処理を行なった場合の処理をそれぞれ示す図である。
図１６（ｄ）は半径Ｒ＝７．５、図１６（ｅ）は半径Ｒ＝８．０、図１６（ｆ）は半径Ｒ＝９．０の場合において、極大点を含む局所面についてさらに電流源の特定処理を行なった場合の処理をそれぞれ示す図である。
図１７は、局所面の深さを変えつつ自由エネルギーを求めた結果を示す図である。
図１８は、脳表面を半球と仮定した場合に、電流源が２つ脳内に存在する場合の脳表面における磁場分布を示す上面図である。
図１９（ａ）は半径Ｒ＝５．０、図１９（ｂ）は半径Ｒ＝６．０、図１９（ｃ）は半径Ｒ＝７．０の場合において、初期分解能を用いた電流源の初期推定結果を示す図である。
図２０（ｄ）は半径Ｒ＝７．５、図２０（ｅ）は半径Ｒ＝８．０、図２０（ｆ）は半径Ｒ＝９．０の場合において、初期分解能を用いた電流源の初期推定結果を示す図である。
図２１は、各仮想半球面ごとに自由エネルギーが最大となるように電流分布を求めた場合の、自由エネルギーの半径依存性を示す図である。
図２２（ａ）は半径Ｒ＝５．０、図２２（ｂ）は半径Ｒ＝６．０、図２２（ｃ）は半径Ｒ＝７．０の場合の局所面上の電流分布をそれぞれ示し、最も脳表面に近い電流源を特定するために、事後確率が最大になる１番目の局所面の深さを求める処理を説明する図である。
図２３（ｄ）は、半径Ｒ＝７．５、図２３（ｅ）は半径Ｒ＝８．０、図２３（ｆ）は半径Ｒ＝９．０の場合の局所面上の電流分布をそれぞれ示し、最も脳表面に近い電流源を特定するために、事後確率が最大になる１番目の局所面の深さを求める処理を説明する図である。
図２４は、仮想局所面を移動させた場合の自由エネルギーの局所面位置（半径）依存性を示す図である。
図２５（ａ）は半径Ｒ＝５．０、図２５（ｂ）は半径Ｒ＝５．５、図２５（ｃ）は半径Ｒ＝６．０の場合の局所面上の電流分布を示し、特定された面上に一方の局所面を固定したままで、他方の電流源に対応する局所面を移動させた場合の電流分布を示す図である。
図２６（ｄ）は、半径Ｒ＝６．５、図２６（ｅ）は半径Ｒ＝７．０、図２６（ｆ）は半径Ｒ＝７．５の場合の局所面上の電流分布をそれぞれ示し、特定された面上に一方の局所面を固定したままで、他方の電流源に対応する局所面を移動させた場合の電流分布を示す図である。
図２７は、仮想局所面を移動させた場合の自由エネルギーの局所面位置（半径）依存性を示す図である。
図２８は、局在条件のみを考慮したシミュレーション結果を示す第１の図である。
図２９は、局在条件のみを考慮したシミュレーション結果を示す第２の図である。
図３０は、連続条件および局在条件の双方を考慮したシミュレーションを行った結果を示す第１の図である。
図３１は、連続条件および局在条件の双方を考慮したシミュレーションを行った結果を示す第２の図である。
図３２は、連続条件および局在条件の双方を考慮し、かつ電流源が存在する可能性の高い領域内で電流が推定されやすくなる階層事前分布を用いた場合のシミュレーションを行った結果を示す第１の図である。
図３３は、連続条件および局在条件の双方を考慮し、かつ電流源が存在する可能性の高い領域内で電流が推定されやすくなる階層事前分布を用いた場合のシミュレーションを行った結果を示す第２の図である。
【図１】
2003057035.tif 000045
【図２】
2003057035.tif 000046
【図３】
2003057035.tif 000047
【図４】
2003057035.tif 000048
【図５】
2003057035.tif 000049
【図６】
2003057035.tif 000050
【図７】
2003057035.tif 000051
【図８】
2003057035.tif 000052
【図９】
2003057035.tif 000053
【図１０】
2003057035.tif 000054
【図１１】
2003057035.tif 000055
【図１２（ａ）】
2003057035.tif 000056
【図１２（ｂ）】
2003057035.tif 000057
【図１２（ｃ）】
2003057035.tif 000058
【図１３（ｄ）】
2003057035.tif 000059
【図１３（ｅ）】
2003057035.tif 000060
【図１３（ｆ）】
2003057035.tif 000061
【図１４】
2003057035.tif 000062
【図１５（ａ）】
2003057035.tif 000063
【図１５（ｂ）】
2003057035.tif 000064
【図１５（ｃ）】
2003057035.tif 000065
【図１６（ｄ）】
2003057035.tif 000066
【図１６（ｅ）】
2003057035.tif 000067
【図１６（ｆ）】
2003057035.tif 000068
【図１７】
2003057035.tif 000069
【図１８】
2003057035.tif 000070
【図１９（ａ）】
2003057035.tif 000071
【図１９（ｂ）】
2003057035.tif 000072
【図１９（ｃ）】
2003057035.tif 000073
【図２０（ｄ）】
2003057035.tif 000074
【図２０（ｅ）】
2003057035.tif 000075
【図２０（ｆ）】
2003057035.tif 000076
【図２１】
2003057035.tif 000077
【図２２（ａ）】
2003057035.tif 000078
【図２２（ｂ）】
2003057035.tif 000079
【図２２（ｃ）】
2003057035.tif 000080
【図２３（ｄ）】
2003057035.tif 000081
【図２３（ｅ）】
2003057035.tif 000082
【図２３（ｆ）】
2003057035.tif 000083
【図２４】
2003057035.tif 000084
【図２５（ａ）】
2003057035.tif 000085
【図２５（ｂ）】
2003057035.tif 000086
【図２５（ｃ）】
2003057035.tif 000087
【図２６（ｄ）】
2003057035.tif 000088
【図２６（ｅ）】
2003057035.tif 000089
【図２６（ｆ）】
2003057035.tif 000090
【図２７】
2003057035.tif 000091
【図２８】
2003057035.tif 000092
【図２９】
2003057035.tif 000093
【図３０】
2003057035.tif 000094
【図３１】
2003057035.tif 000095
【図３２】
2003057035.tif 000096
【図３３】
2003057035.tif 000097
【国際調査報告】
2003057035.tif 000098
2003057035.tif 000099
2003057035.tif 000100
2003057035.tif 000101