(57)【特許請求の範囲】
【請求項１】
各受光プローブの光強度に関する情報であり、基準となる情報である基準光強度情報を格納している基準光強度情報格納部と、
観測対象の物体に対して光を送光する送光プローブと、光を受光する受光プローブの組であるプローブセットを１組以上有する近赤外光計測装置を用いて、前記各プローブセットを構成する各送光プローブから、一の光強度の光を前記物体に対して発した場合に、各受光プローブが受光する光強度に関する情報である光強度情報を取得する光強度情報取得部と、
前記基準光強度情報格納部の基準光強度情報と、前記光強度情報取得部が取得した光強度情報とから、プローブセットごとに、光強度の変化に関する情報である光強度変化情報を取得する光強度変化情報取得部と、
前記プローブセットごとの光強度変化情報を用いて、前記物体の内部の吸光度変化の位置に関する三次元の情報である光学特性変化部位情報を取得する推定部と、
前記推定部が取得した光学特性変化部位情報を出力する出力部とを具備する物体観測装置であって、
前記推定部は、
光学特性変化と光強度変化との関係を示す情報である感度情報を格納しており、当該感度情報を用いて、深さによる感度の減衰に応じた補正を行い、かつ空間解像度を上げるスパース性を導入して前記光学特性変化部位情報を取得する物体観測装置。
【請求項２】
前記推定部は、
光学特性変化と光強度変化との関係を示す情報である感度情報を格納し得る感度情報格納手段と、
前記感度情報を用いて、深さによる感度の減衰に応じた補正を行う補正手段と、
空間解像度を上げるスパース性を導入して前記光学特性変化部位情報を取得するスパース性適用手段とを具備する請求項1記載の物体観測装置。
【請求項３】
前記推定部は、
前記光学特性変化部位情報を取得するために使用する演算式であり、逆問題の解を算出する演算式であり、前記感度情報に応じた補正を行なえるコスト関数である第一の演算式を格納し得る演算式格納手段をさらに具備し、
前記補正手段は、
前記感度情報を前記第一の演算式に代入し、第一の光学特性変化部位情報を取得し、
前記スパース性適用手段は、
前記第一の光学特性変化部位情報を用いて、空間解像度を上げるスパース性を用いて最終的な光学特性変化部位情報を取得し、
前記出力部は、
前記スパース性適用手段が取得した光学特性変化部位情報を出力する請求項２記載の物体観測装置。
【請求項４】
前記演算式格納手段は、
初期値に対して階層ベイズ推定を行う第二の演算式をも格納しており、
前記スパース性適用手段は、
前記第一の光学特性変化部位情報を初期値として前記第二の演算式に代入し、当該第二の演算式を実行することにより階層ベイズ推定を行い、最終的な光学特性変化部位情報を取得する請求項３記載の物体観測装置。
【請求項５】
前記観測対象の物体は、
生物の脳であり、
前記基準光強度情報は、
安静状態において測定された各受光プローブの光強度に関する情報である、生物の脳機能を観測する請求項１記載の物体観測装置。
【請求項６】
前記観測対象の物体は、
乳房であり、
前記基準光強度情報は、
正常な細胞の乳房に対する各受光プローブの光強度に関する情報である、乳がんの検診を行える請求項１記載の物体観測装置。
【請求項７】
前記各送受光プローブの間隔は、２０ｍｍ以下である請求項１から請求項６いずれか記載の物体観測装置。
【請求項８】
記録媒体に、
各受光プローブの光強度に関する情報であり、基準となる情報である基準光強度情報を格納しており、
光強度情報取得部、光強度変化情報取得部、推定部、および出力部により実現される物体観測方法であって、
前記光強度情報取得部が、観測対象の物体に対して光を送光する送光プローブと、光を受光する受光プローブの組であるプローブセットを１組以上有する近赤外光計測装置を用いて、前記各プローブセットを構成する各送光プローブから、一の光強度の光を前記物体に対して発した場合に、各受光プローブが受光する光強度に関する情報である光強度情報を取得する光強度情報取得ステップと、
前記光強度変化情報取得部が、前記記録媒体の基準光強度情報と、前記光強度情報取得ステップで取得された光強度情報とから、プローブセットごとに、光強度の変化に関する情報である光強度変化情報を取得する光強度変化情報取得ステップと、
前記推定部が、前記プローブセットごとの光強度変化情報を用いて、前記物体の内部の吸光度変化の位置に関する三次元の情報である光学特性変化部位情報を取得する推定ステップと、
前記出力部が、前記推定ステップで取得された光学特性変化部位情報を出力する出力ステップとを具備する物体観測方法であって、
前記推定ステップにおいて、
光学特性変化と光強度変化との関係を示す情報である感度情報を用いて、深さによる感度の減衰に応じた補正を行い、かつ空間解像度を上げるスパース性を導入して前記光学特性変化部位情報を取得する物体観測方法。
【請求項９】
記録媒体に、
各受光プローブの光強度に関する情報であり、基準となる情報である基準光強度情報を格納しており、
コンピュータを、
観測対象の物体に対して光を送光する送光プローブと、光を受光する受光プローブの組であるプローブセットを１組以上有する近赤外光計測装置を用いて、前記各プローブセットを構成する各送光プローブから、一の光強度の光を前記物体に対して発した場合に、各受光プローブが受光する光強度に関する情報である光強度情報を取得する光強度情報取得部と、
前記記録媒体の基準光強度情報と、前記光強度情報取得部が取得した光強度情報とから、プローブセットごとに、光強度の変化に関する情報である光強度変化情報を取得する光強度変化情報取得部と、
前記プローブセットごとの光強度変化情報を用いて、前記物体の内部の吸光度変化の位置に関する三次元の情報である光学特性変化部位情報を取得する推定部と、
前記推定部が取得した光学特性変化部位情報を出力する出力部として機能させるためのプログラムを記録した記録媒体であって、
前記推定部は、
光学特性変化と光強度変化との関係を示す情報である感度情報を用いて、深さによる感度の減衰に応じた補正を行い、かつ空間解像度を上げるスパース性を導入して前記光学特性変化部位情報を取得するものとして、コンピュータを機能させるためのプログラムを記録した記録媒体。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、物体内部の局所的な光学特性変化部位を示す３次元情報を出力する物体観測装置等に関するものである。
【背景技術】
【０００２】
近赤外光を用いた計測は、簡便で安全であることにより、非侵襲な脳活動計測（非特許文献１参照）や乳がんの検知（非特許文献２参照）に用いられている。また、近赤外光は生体内で強く散乱されるが、その拡散過程を計算することにより、観測データから生体内部の光学特性変化や構造を再構成しようとする試みが拡散光トモグラフィ（ＤＯＴ）である。また、近赤外光計測には、連続光計測、時間分解計測、周波数分解計測があり、連続光計測は得られる情報が少ないが、コストが安く、照射する光量が多いためにノイズにロバストであり、サンプリング周波数が高いという特徴がある。一方、時間分解計測や周波数分解計測は、得られる情報が豊富であるが、ノイズの影響が大きく、コストが高いという特徴がある。【０００３】
近赤外光は生体内で強く散乱されるため、光を透過させることができない場合が多い。そのため例えば脳機能計測では、頭部表面に送光プローブを置き、拡散により返ってきた光を、同じく頭部表面に置いた受光プローブで観測する反射型の計測を行う。この脳機能計測では、活動領域をカバーするのに多くのプローブを用いた同時計測が必要となるため連続光計測が一般的である。【０００４】
また、従来は脳活動計測に良いとされる３ｃｍ程度の間隔でプローブを配置して、観測結果をトポグラフィとして捉えていたが（非特許文献３参照）、近年の脳機能計測においては、単位面積当たりのプローブ数を増やし、高密度にしてＤＯＴを行う試みが盛んに行われている。それにより、視覚野や体性感覚野における複数の皮質活動領域の識別に成功している（非特許文献４、非特許文献５参照）。
【先行技術文献】
【非特許文献】
【０００５】
【非特許文献１】Y. Hoshi and M. Tamura,"Detection of dynamic changes in cerebral oxygenation coupled to neuronal function during mental work in man," Neurosci. Lett. 150 (1993) 5-8.
【非特許文献２】山下  大輔，小田  元樹，鈴木  裕昭，山中  武志，鈴木  俊彦，上田  之雄，太田  和義，三輪  光春，長谷川  寛，久嶋  浩之，山下    豊「多チャンネル時間分解分光計測システムの開発と乳癌検診への応用」J Jpn Coll Angiol, 47 (2007) 11-16.
【非特許文献３】A. Maki, Y. Yamashita, Y. Ito, E. Watanabe, Y. Mayanagi, and H. Koizumi, "Spatial and temporal analysis of human motor activity using noninvasive NIR topography," Med. Phys. 22 (1995) 1997-2005.
【非特許文献４】B. W. Zeff, B. R. White, H. Dehghani, B. L. Schlaggar, and J. P. Culver, "Retinotopic mapping of adult human visual cortex with high-density diffuse optical tomography," Proc. Natl. Acad. Sci. U.S.A. 104 (2007) 12169-12174.
【非特許文献５】C. Habermehl, S. Holtze, J. Steinbrink, S. P. Koch, H. Obrig, J. Mehnert, and C. H. Schmitz,"Somatosensory activation of two fingers can be discriminated with ultrahigh-density diffuse optical tomography," NeuroImage 59 (2011) 3201-3211.
【発明の概要】
【発明が解決しようとする課題】
【０００６】
しかしながら、上述した反射型の連続光計測は、ＤＯＴの中でも得られる情報が最も乏しく、それにより活動源推定にあたっての不良設定性が高く、推定は困難である。反射型連続光ＤＯＴにおける課題は、主に以下の２点が挙げられる。【０００７】
１点目は、反射型連続光ＤＯＴは反射型であるために、表面にしかプローブを置くことができないため、および光が拡散されるにあたって深さとともに感度が指数関数的に減衰してしまうために、活動源の深さを特定することが難しい。【０００８】
また、２点目は、反射型連続光ＤＯＴは光が物体内で拡散されるために空間的にぼやけてしまい、空間分解能を上げにくい。【０００９】
つまり、従来の物体観測システムにおいては、脳等の物体内部の局所的な光学特性変化部位を三次元情報として高精度に推定することができなかった。
【課題を解決するための手段】
【００１０】
本発明の物体観測装置において、上記１点目の課題に対しては、深さによる感度の減衰に応じた補正を行う。また、上記の２点目の課題に対しては、空間解像度を上げる方法としてスパース推定を行う。本発明の物体観測装置において、二つの対策を組み合わせることにより、深さ方向を含めた高精度三次元ＤＯＴを実現する。【００１１】
さらに具体的には、本第一の発明の物体観測装置は、各受光プローブの光強度に関する情報であり、基準となる情報である基準光強度情報を格納している基準光強度情報格納部と、観測対象の物体に対して光を送光する送光プローブと、光を受光する受光プローブの組であるプローブセットを１組以上有する近赤外光計測装置を用いて、各プローブセットを構成する各送光プローブから、一の光強度の光を物体に対して発した場合に、各受光プローブが受光する光強度に関する情報である光強度情報を取得する光強度情報取得部と、基準光強度情報格納部の基準光強度情報と、光強度情報取得部が取得した光強度情報とから、プローブセットごとに、光強度の変化に関する情報である光強度変化情報を取得する光強度変化情報取得部と、プローブセットごとの光強度変化情報を用いて、物体の内部の吸光度変化の位置に関する三次元の情報である光学特性変化部位情報を取得する推定部と、推定部が取得した光学特性変化部位情報を出力する出力部とを具備する物体観測装置であって、推定部は、光学特性変化と光強度変化との関係を示す情報である感度情報を格納し得る感度情報格納手段と、感度情報を用いて、深さによる感度の減衰に応じた補正を行う補正手段と、空間解像度を上げるスパース性を導入して光学特性変化部位情報を取得するスパース性適用手段とを具備する物体観測装置である。なお、推定部は、上述した感度の減衰に応じた補正と、空間解像度を上げるスパース性の適用を同時に行っても良い。【００１２】
かかる構成により、物体内部の局所的な光学特性変化部位を三次元情報として高精度に推定できる。【００１３】
また、本第二の発明の物体観測装置は、第一の発明に対して、推定部は、光学特性変化部位情報を取得するために使用する演算式であり、逆問題の解を算出する演算式であり、感度情報に応じた補正を行なえるコスト関数である第一の演算式を格納し得る演算式格納手段をさらに具備し、補正手段は、感度情報を第一の演算式に代入し、第一の光学特性変化部位情報を取得し、スパース性適用手段は、第一の光学特性変化部位情報を用いて、空間解像度を上げるスパース性を用いて最終的な光学特性変化部位情報を取得し、出力部は、スパース性適用手段が取得した光学特性変化部位情報を出力する物体観測装置である。【００１４】
かかる構成により、物体内部の局所的な光学特性変化部位を三次元情報として高精度に推定できる。【００１５】
また、本第三の発明の物体観測装置は、第二の発明に対して、演算式格納手段は、初期値に対して階層ベイズ推定を行う第二の演算式をも格納しており、スパース性適用手段は、第一の光学特性変化部位情報を初期値として第二の演算式に代入し、第二の演算式を実行することにより階層ベイズ推定を行い、最終的な光学特性変化部位情報を取得する物体観測装置である。【００１６】
かかる構成により、物体内部の局所的な光学特性変化部位を三次元情報として高精度に推定できる。【００１７】
また、本第四の発明の物体観測装置は、第一から第三いずれかの発明に対して、観測対象の物体は、生物の脳であり、基準光強度情報は、安静状態において測定された各受光プローブの光強度に関する情報である、生物の脳機能を観測する物体観測装置である。【００１８】
かかる構成により、脳活動を高精度に容易に推定できる。【００１９】
また、本第五の発明の物体観測装置は、第一から第三いずれかの発明に対して、観測対象の物体は、乳房であり、基準光強度情報は、正常な細胞の乳房に対する各受光プローブの光強度に関する情報である、乳がんの検診を行える物体観測装置である。【００２０】
かかる構成により、乳がんを高精度に容易に検知できることができる。【００２１】
また、本第六の発明の物体観測装置は、第一から第五いずれかの発明に対して、各送受光プローブの間隔は、２０ｍｍ以下である物体観測装置である。【００２２】
かかる構成により、物体内部の局所的な光学特性変化部位を三次元情報としてより高精度に推定できる。
【発明の効果】
【００２３】
本発明による物体観測システムによれば、物体内部の局所的な光学特性変化部位を三次元情報として高精度に推定できる。
【図面の簡単な説明】
【００２４】
【図１】実施の形態１における物体観測システムのブロック図
【図２】同物体観測装置の動作について説明するフローチャート
【図３】同推定処理の動作について説明するフローチャート
【図４】同物体観測システムの実験の模式図
【図５】同ファントム実験の送光プローブ、受光プローブの位置を説明する図
【図６】同実験１の環境を説明する図
【図７】同正則化項を感度で規格化しない場合のＴｉｋｈｏｎｏｖ正則化法の推定結果を示す図
【図８】同正規化項を感度で規格化したときの推定結果を示す図
【図９】同階層ベイズ法に取り入れた場合の推定結果を示す図
【図１０】同プローブ間隔と吸光体の深さ位置を変えながら行った推定結果を示す図
【図１１】同３次元推定図を２次元で表すことを説明する図
【図１２】同実験２の結果を示す図
【図１３】同実験２の結果を示す図
【図１４】同実験３の結果を示す図
【図１５】同実験３の結果を示す図
【図１６】同実験４の結果を示す図
【図１７】同実験４の結果を示す図
【図１８】同コンピュータシステムの概観図
【図１９】同コンピュータシステムのブロック図
【発明を実施するための形態】
【００２５】
以下、物体観測システム等の実施形態について図面を参照して説明する。なお、実施の形態において同じ符号を付した構成要素は同様の動作を行うので、再度の説明を省略する場合がある。
【００２６】
（実施の形態１）
本実施の形態において、物体に対して、近赤外光を送光し、かつ物体内部で散乱された光を受光し、物体内部の局所的な光学特性変化部位を示す三次元情報を出力する物体観測システムについて説明する。ここでの物体観測システムを構成する物体観測装置は、通常、物体を透過した光を利用するのではなく、物体の内部からの後方散乱光を受光して利用する。
【００２７】
図１は、本実施の形態における物体観測システムのブロック図である。物体観測システムは、近赤外光計測装置１、および物体観測装置２を備える。
【００２８】
近赤外光計測装置１は、送光プローブ１１、受光プローブ１２、受付部１３、および制御部１４を備える。近赤外光計測装置１は、送光プローブ１１と受光プローブ１２の組であるプローブセット（チャネルとも言う）を１組以上有する。なお、近赤外光計測装置１は、通常、複数組のプローブセットを有する。
【００２９】
物体観測装置２は、基準光強度情報格納部２１、光強度情報取得部２２、光強度変化情報取得部２３、推定部２４、および出力部２５を備える。また、推定部２４は、感度情報格納手段２４１、演算式格納手段２４２、補正手段２４３、およびスパース性適用手段２４４を備える。
【００３０】
近赤外光計測装置１を構成する送光プローブ１１は、近赤外光を物体（ここでは、例えば、生体の頭部や乳房等）に照射する。
【００３１】
受光プローブ１２は、通常、物体内部で吸収あるいは散乱され、物体内から出てくる光を受光する。なお、プローブの配置密度が２０ｍｍ以下であることはさらに好適である。さらに、プローブの配置密度が１．５ｃｍ以下であることは極めて好適である。プローブの配置密度とは、送光プローブ１１と受光プローブ１２間の間隔である。なお、プローブの配置密度は密であればある程良いのは、物体内部の状態をとらえるためのチャネルがたくさんできること等からである。
【００３２】
なお、送光プローブ１１は、観測対象の物体に対して一の方向から光を送光する。また、受光プローブ１２は、通常、前記一の方向と同一の方向から光を受光する。なお、ここでの同一の方向とは、厳密に一致する方向とは限らず、角度のずれ（例えば、３０度、４５度等の大きなずれ）が存在しても良い。同一の方向とは、通常、送光プローブ１１が光を送光する側と受光プローブ１２が光を受光する側とが、物体の同一の面側であることである。
【００３３】
受付部１３は、外部からの指示を受け付ける。かかる指示は、近赤外光を生体に照射する指示である照射指示などである。
【００３４】
制御部１４は、受付部１３が照射指示を受け付けた場合に、送光プローブ１１に近赤外光の照射を指示する。なお、送光プローブ１１は、通常、制御部１４から照射指示を受け付けた場合に、近赤外光を物体に照射する。
【００３５】
物体観測装置２を構成する基準光強度情報格納部２１は、基準光強度情報を格納している。基準光強度情報は、各受光プローブの光強度に関する情報であり、基準となる情報である。例えば、物体観測装置２が脳活動を観測する場合は、基準光強度情報は、安静状態において測定された各受光プローブ１２の光強度に関する情報である。また、乳がんを検知する場合は、基準光強度情報は、正常な細胞に対する各受光プローブ１２の光強度に関する情報である。なお、基準光強度情報が基準光強度情報格納部２１に格納される過程は問わない。また、基準光強度情報と、後述する光強度情報とが同時期に取得されても良い。
【００３６】
光強度情報取得部２２は、観測対象の物体に対して光を送光する送光プローブ１１と、光を受光する受光プローブ１２の組であるプローブセットを１組以上有する近赤外光計測装置１を用いて、光強度情報を取得する。この光強度情報は、各プローブセットを構成する各送光プローブ１１から、一の光強度の光を物体に対して発した場合に、各受光プローブが受光する光強度に関する情報である。なお、観測対象の物体とは、例えば、生物の脳や、乳房等である。
【００３７】
光強度変化情報取得部２３は、基準光強度情報格納部２１の基準光強度情報と、光強度情報取得部２２が取得した光強度情報とから、プローブセットごとに、光強度の変化に関する情報である光強度変化情報を取得する。
【００３８】
推定部２４は、プローブセットごとの光強度変化情報を用いて、光学特性変化部位情報を取得する。なお、この光強度変化情報は、光強度変化情報取得部２３が取得した情報である。また、光学特性変化部位情報とは、物体の内部の吸光度変化の位置に関する三次元の情報である。また、光学特性変化部位情報は、物体内部の吸光度の変化量に関する情報である、とも言える。
【００３９】
推定部２４を構成する感度情報格納手段２４１は、感度情報を格納し得る。感度情報は、光学特性変化（例えば、後述するＸが対応）と光強度変化（例えば、後述するＹが対応）との関係を示す情報である。感度情報は、例えば、後述する感度行列である。
【００４０】
演算式格納手段２４２は、物体観測装置２が利用する演算式を格納し得る。演算式格納手段２４２は、例えば、第一の演算式または第二の演算式を格納し得る。第一の演算式は、例えば、光学特性変化部位情報を取得するために使用する演算式であり、逆問題の解を算出する演算式であり、感度情報に応じた補正を行なえるコスト関数である。また、第二の演算式は、例えば、初期値に対して階層ベイズ推定を行う演算式である。第二の演算式は、例えば、Ｌ１ノルムの演算式であっても良い。
【００４１】
補正手段２４３は、感度情報格納手段２４１の感度情報を用いて、深さによる感度の減衰に応じた補正を行う。さらに具体的には、補正手段２４３は、感度情報を第一の演算式に代入し、第一の光学特性変化部位情報を取得する。なお、深さによる感度の減衰に応じた補正を行うための演算式は問わない。
【００４２】
スパース性適用手段２４４は、空間解像度を上げるスパース性を導入して光学特性変化部位情報を取得する。具体的には、スパース性適用手段２４４は、第一の光学特性変化部位情報を入力として、空間解像度を上げるスパース性を用いて最終的な光学特性変化部位情報を取得する。さらに具体的には、スパース性適用手段２４４は、第一の光学特性変化部位情報を初期値として第二の演算式に代入し、第二の演算式を実行することにより階層ベイズ推定を行い、最終的な光学特性変化部位情報を取得する。なお、スパース性の導入とは、階層ベイズ推定が好適であるが、Ｌ１ノルムを用いても良い。スパース性を導入する場合のアルゴリズムは問わない。また、アルゴリズムにおいて補正手段２４３とスパース性適用手段２４４とを同時に実行しても良い。
【００４３】
さらに具体的には、推定部２４は、以下のように動作する。つまり、推定部２４は、ベースラインで計測した光量に関する情報である基準光強度情報と、吸光度変化が起きたときに計測した光量に関する情報（光強度情報）とを用いて、二つの光量の変化値である光強度変化情報から、３次元的な吸光度変化マップを再構成する。
【００４４】
推定部２４は、上記の再構成を行うために、順問題を解くステップと逆問題を解くステップの二つのステップを採る。まず、順問題を解くステップにより、推定部２４は、吸光度変化と観測光強度変化の関係性を取得する。次に、逆問題を解くステップにより、推定部２４は、順問題で求めた関係性と、観測光強度変化から、吸光度変化を推定する。なお、推定する値（例えば、Ｘ）の数より少ない観測値（例えば、Ｙ）を用いて、値（例えば、Ｘ）を推定し、取得する場合の問題は逆問題の中でも不良設定問題と呼ばれ、解くためには正則化等の何らかの工夫が必要となる。
【００４５】
ここで、物体観測装置２は、順問題を解くステップではＲｙｔｏｖ近似を用いる。なお、Ｒｙｔｏｖ近似については、公知技術であるので詳細な説明は省略する。また、Ｒｙｔｏｖ近似により、ベースラインにおいて、ソース（送光プローブ１１）のポジションｒｓから出た光がディテクター（受光プローブ１２）のポジションｒｄに届く光量Φ０（ｒｓ，ｒｄ）が、物体内部の吸光度がδμａだけ変化することにより、Φ（ｒｓ，ｒｄ）になるとしたとき、生じる光量変化Φｐｅｒｔ（ｒｓ，ｒｄ）＝ｌｎ（Φ０（ｒｓ，ｒｄ）／Φ（ｒｓ，ｒｄ））は、以下の数式１により記述される。つまり、数式１は、Ｒｙｔｏｖ近似を示す。
【数１】
0005570674.tif 000002
【００４６】
ここで、本実施の形態において、物体観測装置２は、領域を「２．５ｍｍ×２．５ｍｍ×２．５ｍｍ」のボクセルに離散化する。ボクセルのインデックスをｉ、チャネルのインデックスをｊとし、「Ｙｊ＝Φｐｅｒｔ（ｒｓ（ｊ），ｒｄ（ｊ）），Ｘｉ＝δμａ（ｒｉ），Ａｊ，ｉ＝Φ０（ｒｓ（ｊ），ｒｉ）Φ０（ｒｉ，ｒｄ（ｊ））／Φ０（ｒｓ（ｊ），ｒｄ（ｊ）））」とすると、数式１のＲｙｔｏｖ近似式は「Ｙ＝ＡＸ」と線形方程式で記述される。ここで、「Ｙ＝（Ｙ１，・・・，ＹＭ）Ｔ」は全チャネルにおける観測光強度変化、「Ｘ＝（Ｘ１，・・・，ＸＮ）Ｔ」は全ボクセルにおける吸光度変化、ＡはＡｊ，ｉを要素とする感度行列（感度情報の例）である。なお、感度行列は、感度情報格納手段２４１に格納されている。また、Ｒｙｔｏｖ近似式、または「Ｙ＝ＡＸ」は演算式格納手段２４２に格納されていても良い。なお、以下、適宜、ボクセルの個数はＮ、チャネルの個数はＭとする。
【００４７】
Φ０の値は光子モンテカルロシミュレーションで得ることができるが（「Q. Fang and D. A. Boas,"Monte Carlo simulation of photon migration in 3D turbid media accelerated by graphics processing units,"Opt. Express 17 (2009) 20178-20190.」参照）、媒質の形状によっては拡散近似を用いた解析解でもΦ０の良い近似値を得ることができる（「R. C. Haskell, L. O. Svaasand, T. T. Tsay, T. C. Feng, M. S. McAdams, and B. J. Tromberg,
"Boundary conditions for the diffusion equation in radiative transfer,"J. Opt. Soc. Am. A 11 (1994) 2727-2741.参照）。
【００４８】
次に、推定部２４は、逆問題を以下のように解く。つまり、推定部２４を構成する補正手段２４３は、正則化項を感度項により規格化したＴｉｋｈｏｎｏｖ正則化法の推定解を得る。次に、スパース性適用手段２４４は、補正手段２４３が得た値を初期値として階層ベイズ推定を行う。
【００４９】
まず、補正手段２４３の動作について詳細に説明する。演算式格納手段２４２に、以下の数式２に示すコスト関数が格納されている。そして、補正手段２４３は、コスト関数を最小とするものを解「Ｘ?R＝ａｒｇｍｉｎＸＣＲ（Ｘ；λ）」を算出する。なお、「Ｘ?」の「＾」は、Ｘの真上に存在するものとする。
【数２】
0005570674.tif 000003
【００５０】
数式２の第一項が観測誤差項、第二項がオーバーフィッティングを防ぐためのＬ２正則化項である。なお、オーバーフィッティングとは、観測データにパラメータを過度に最適化することである。また、数式２において、λは正則化の程度を決定するパラメータである。最適なλの値は情報量基準の一つであるＡＢＩＣ最小化により決定され、演算式格納手段２４２に格納されている（「H. Akaike, "Likelihood and the Bayes procedure,"In Bayesian Statistics (eds J. M. Bernardo et al., Univ. Press, Valencia) 143-166 (1980).」参照）。また、Σｙは観測値Ｙの観測ノイズ共分散行列である。
【００５１】
ここで、数式２のＤは感度に応じた規格化項である。感度行列Ａの値は、深さとともに指数減衰するという特徴を持つ。そのため、空間的に一様な正則化を行うと推定結果にバイアスが生じるという問題がある。規格化項は、この問題を防ぐために、各場所における正則化項の強さを感度の大きさと釣り合わせる。なお、数式２のＤは、以下の数式３で示される。
【数３】
0005570674.tif 000004
【００５２】
数式３のｄｉａｇは対角成分から成る行列を表し、右肩のＴは行列の転置を表す。Ｉは単位行列である。βは規格化を行う感度の閾値である。本実施の形態では、以下のすべての解析においてβの値は、例えば、プローブ面から２２．５ｍｍ以上の深さにあるボクセルにおける感度の最大値、「β＝ｍａｘｚ＞２２．５ｍｍ［ｄｉａｇ（ＡＴΣ−１ｙＡ）］」と定める。また、正則化項がＬ２ノルムの形をしているときには、数式４に示すように、推定解Ｘ＾R（＾はＸの真上に存在する）は解析的に得ることができる。つまり、数式４は、演算式格納手段２４２に格納されている。そして、補正手段２４３は、数式４を用いて、推定解Ｘ＾Rを得る。
【数４】
0005570674.tif 000005
【００５３】
次に、スパース性適用手段２４４の動作について詳細に説明する。なお、正則化項を感度項により規格化したＴｉｋｈｏｎｏｖ正則化法による推定解Ｘ＾Rはいくつかの問題がある。それらの問題は、解を局在化できていないことにより、感度に応じた補正が不必要な部分に対し強く働き、その結果推定された位置に大きな誤差が生じるといった問題と、空間解像度が低いという問題である。そこで、解をさらに改善するために、スパース性適用手段２４４は、推定解Ｘ＾Rを初期値とした階層ベイズ推定を行う。つまり、推定解がなるべく少ない領域であらわされるとする「スパース性」を導入する。補正手段２４３の動作により、感度を補正して得られた結果から、重要な場所を絞り込むことにより、高空間解像度を達成できる。かかる動作により、深さもより正しく推定される。
【００５４】
階層ベイズ推定は、観測値の生成だけでなく、正則化についても以下のように階層確率モデルとして表す。まず、内部吸光度変化Ｘにおける観測光強度変化Ｙの観測モデルを以下の数式５に示す。
【数５】
0005570674.tif 000006
【００５５】
数式５において、平均μと共分散行列Σの正規分布をＮ（ｘ｜μ，Σ）と表記した。Σｙは観測値Ｙの観測ノイズ共分散行列である。さらに、解Ｘは平均０、分散が調節可能な独立なパラメータΛ＝ｄｉａｇ（λ１，λ２，・・・，λＮ）から成る正規分布で与えられるものとする。
【数６】
0005570674.tif 000007
【００５６】
このような事前分布は関連度自動決定（Ａｕｔｏｍａｔｉｃ  ｒｅｌｅｖａｎｃｅ  ｄｅｔｅｒｍｉｎａｔｉｏｎ，ＡＲＤ）事前分布と言う。さらに、Λは以下の数式７に示す確率分布から生じるものである。
【数７】
0005570674.tif 000008
【００５７】
数式７において、平均μ，シェイプパラメータγのガンマ分布をΓ（ｘ｜μ，γ）とする。λｐｉは事前分布の平均であり、γｐｉはその信頼度パラメータである。ここで、例えば、「γpｉ＝０」とする。そして、「γpｉ＝０」の場合、数式７は無情報事前分布となり、λpｉは確率モデルには影響を与えない。
【００５８】
ベイズ法では観測値Ｙが与えられた下での、Ｘ，Λの事後分布Ｐ（Ｘ，Λ｜Ｙ）を求める（数式８参照）。
【数８】
0005570674.tif 000009
【００５９】
しかし、事後分布は解析的には求まらない。事後分布の計算は後述するように変分ベイズ法を用いて自由エネルギーの最大化に帰着させることができる。自由エネルギーは以下の数式９であらわされる。数式９のf(Λ)は事前分布項をあらわす。自由エネルギーは後述するように反復計算により最大化させることができる。
【数９】
0005570674.tif 000010
【００６０】
Ｔｉｋｈｏｎｏｖ正則化法、数式２では、すべての内部状態Ｘｉに共通の正則化パラメータλを用い、一様なペナルティをかけることによって、スムーズ化された解が得られる。一方、階層ベイズ推定、数式９ではボクセル毎に異なる正則化パラメータΛ＝ｄｉａｇ（λ１，・・・，λＮ）を用い、それをデータに適応させることによりスパース化の効果が得られる（「A. C. Faul and M. E. Tipping, "Analysis of sparse Bayesian learning," Adv. Neural Inf. Process. Syst. 14 (2002) 383-389.」参照）。スパース性適用手段２４４は、後述する反復計算の初期値λ０ｉにＴｉｋｈｏｎｏｖ正則化法により得られた推定解Ｘ＾Rを代入し、推定を行う。また、正確にいえば、スパース性適用手段２４４は、数式５においてノイズ項のスケール補正と、数式６において平滑化フィルタの導入を行っているが、これらは後述する。
【００６１】
以上が階層ベイズ推定法の骨子である。以下、詳細な設定を説明する。まず、計測に当たって時間サンプルデータがＴ個得られる。そのため「Ｙ＝（Ｙ（１），・・・，Ｙ（ｔ），・・・，Ｙ（Ｔ））」と拡張し、それに応じて吸光度変化も以下のように時間に応じて「Ｘ＝（Ｘ（１），・・・，Ｘ（ｔ），・・・，Ｘ（Ｔ））」と拡張する。また、観測ノイズ共分散行列はベースラインで測定されたノイズ共分散行列の定数倍「Σｙ＝σ−１Σｙ０」とし、さらにその係数σに対し事前分布を与え、その値は推定により決定するものとする（数式１０参照）。
【数１０】
0005570674.tif 000011
【００６２】
また、非線形効果が強い場合、例えばプローブ間隔１３ｍｍの場合においては、それをノイズとして組み込み、観測ノイズΣｙの対角項に観測値の二乗値Ｙ２に比例する項を加える。また、解がスパースになりすぎるのを防ぐため、解Ｘは、解のスパースな表現Ｚに平滑化フィルタＷを掛けたものである「Ｘ＝ＷＺ」とし、そのＺに対して関連度自動決定事前分布を適用する。これによりＸの確率分布は以下で与えられる、
【数１１】
0005570674.tif 000012
【００６３】
数式１１は、数式６をより適切な形に拡張したものである。なお、ここでは、例えば、Ｗは全半値幅５ｍｍのガウシアンフィルタを用いる。かかる場合、事後分布は、以下の数式１２、数式１３で与えられる。
【数１２】
0005570674.tif 000013
【数１３】
0005570674.tif 000014
【００６４】
事後分布は、解析的に計算できないために変分ベイズ法を用いる。変分ベイズ法では、事後分布を直接計算する代わりに、試験事後分布Ｑ（Ｘ，Λ，σ）を用意し、真の事後分布Ｐ（Ｘ，Λ，σ｜Ｙ）を近似させる。この近似は、以下の数式１４、数式１５で定義される自由エネルギーＦ（Ｑ）の最大化によって行われる。
【数１４】
0005570674.tif 000015
【数１５】
0005570674.tif 000016
【００６５】
なお、数式１５は、数式１４を書き換えたものである。数式１５におけるKLはＫｕｌｌｂａｃｋ−Ｌｅｉｂｌｅｒ距離をあらわす。数式１５の第一項はＱに依存しないので、自由エネルギー最大化は、第二項にある真の事後分布と試験事後分布とのＫｕｌｌｂａｃｋ−Ｌｅｉｂｌｅｒ距離の最小化と等価である。
【００６６】
ここで、試験事後分布におけるＸとΛ，σの間に独立分解仮定をおくと、計算を簡単に行うことが可能になる（数式１６参照）。
【数１６】
0005570674.tif 000017
【００６７】
この仮定により自由エネルギーの最大化は、ＱＸに関してＦ（Ｑ）を最大化する「Ｘ最大化ステップ」と，ＱΛに関してＦ（Ｑ）を最大化する「Λ，σ最大化ステップ」を交互に繰り返すことにより実現できる。
【００６８】
また、この場合、自由エネルギーは以下の数式１７であらわされる。
【数１７】
0005570674.tif 000018
【００６９】
なお、数式９は、数式１７を簡略化したものである。ただし、数式１７におけるＸ，Λ，σは分布の期待値をあらわす。
【００７０】
次に、スパース性適用手段２４４が行う階層ベイズ推定のアルゴリズムについて説明する。
【００７１】
初期値λ０ｉ＝（１０・Ｘ＾Rｉ）-２，σ０＝１で、スパース性適用手段２４４は、以下の「Ｘ最大化ステップ」および「Λ，σ最大化ステップ」をＫ（例えば、Ｋ＝１０００）回繰り返す。
【００７２】
なお、「Ｘ最大化ステップ」は、以下の数式１８から数式２０で示される。また、「Λ，σ最大化ステップ」は、以下の数式２１から数式２４で示される。
【００７３】
そして、スパース性適用手段２４４は、最後に「Ｘ最大化ステップ」をして、解Ｘ?B（数式２５参照）を得る。
【数１８】
0005570674.tif 000019
【数１９】
0005570674.tif 000020
【数２０】
0005570674.tif 000021
【数２１】
0005570674.tif 000022
【数２２】
0005570674.tif 000023
【数２３】
0005570674.tif 000024
【数２４】
0005570674.tif 000025
【数２５】
0005570674.tif 000026
なお、数式２１，２３のＴは時間サンプルデータの個数である。また、数式２３において、Ｍはチャネルの個数である。
【００７４】
出力部２５は、スパース性適用手段２４４が取得した光学特性変化部位情報を出力する。出力部２５は、通常、図的に光学特性変化部位情報を出力する。ただし、光学特性変化部位情報の出力の態様は問わない。
【００７５】
出力とは、ディスプレイへの表示、プロジェクターを用いた投影、プリンタでの印字、音出力、外部の装置への送信、記録媒体への蓄積、他の処理装置や他のプログラムなどへの処理結果の引渡しなどを含む概念である。
【００７６】
基準光強度情報格納部２１、感度情報格納手段２４１、および演算式格納手段２４２は、不揮発性の記録媒体が好適であるが、揮発性の記録媒体でも実現可能である。
【００７７】
基準光強度情報格納部２１等に基準光強度情報等が記憶される過程は問わない。例えば、記録媒体を介して基準光強度情報等が基準光強度情報格納部２１等で記憶されるようになってもよく、通信回線等を介して送信された基準光強度情報等が基準光強度情報格納部２１等で記憶されるようになってもよく、あるいは、入力デバイスを介して入力された基準光強度情報等が基準光強度情報格納部２１等で記憶されるようになってもよい。
【００７８】
光強度情報取得部２２、光強度変化情報取得部２３、推定部２４、補正手段２４３およびスパース性適用手段２４４は、通常、ＣＰＵやメモリ等から実現され得る。光強度情報取得部２２等の処理手順は、通常、ソフトウェアで実現され、当該ソフトウェアはＲＯＭ等の記録媒体に記録されている。但し、ハードウェア（専用回路）で実現しても良い。
【００７９】
出力部２５は、ディスプレイ等の出力デバイスを含むと考えても含まないと考えても良い。出力部２５は、出力デバイスのドライバーソフトまたは、出力デバイスのドライバーソフトと出力デバイス等で実現され得る。
【００８０】
次に、物体観測装置２の動作について、図２のフローチャートを用いて説明する。
【００８１】
（ステップＳ２０１）光強度情報取得部２２は、各プローブセットを構成する各送光プローブから一の光強度の光を物体に対して発した場合に、各受光プローブが受光する光強度に関する情報である光強度情報を取得する。
【００８２】
（ステップＳ２０２）光強度変化情報取得部２３は、基準光強度情報格納部２１から基準光強度情報を読み出す。
【００８３】
（ステップＳ２０３）光強度変化情報取得部２３は、ステップＳ２０１で取得された光強度情報と、ステップＳ２０２で読み出された基準光強度情報とから、光強度変化情報を算出する。
【００８４】
（ステップＳ２０４）推定部２４は、演算式格納手段２４２から演算式を読み出す。
【００８５】
（ステップＳ２０５）推定部２４は、感度情報格納手段２４１から感度情報を読み出す。なお、例えば、推定部２４は予め感度情報を算出し、感度情報格納手段２４１に格納している、とする。
【００８６】
（ステップＳ２０６）推定部２４は、光強度変化情報を用いて、光学特性変化部位情報を推定する。かかる処理を推定処理と言う。推定処理について、図３のフローチャートを用いて説明する。
【００８７】
（ステップＳ２０７）出力部２５は、スパース性適用手段２４４が取得した光学特性変化部位情報を出力する。ステップＳ２０１に戻る。
【００８８】
なお、図２のフローチャートにおいて、電源オフや処理終了の割り込みにより処理は終了する。
【００８９】
次に、ステップＳ２０６の推定処理について、図３のフローチャートを用いて説明する。  
【００９０】
（ステップＳ３０１）補正手段２４３は、感度情報格納手段２４１の感度情報を用いて、深さによる感度の減衰に応じた補正を行う。つまり、補正手段２４３は、感度情報を第一の演算式（感度に応じた補正を入れた光学特性変化部位情報の推定処理を行う演算式）に代入し、第一の光学特性変化部位情報を取得する。
【００９１】
（ステップＳ３０２）スパース性適用手段２４４は、ステップＳ３０１で取得された第一の光学特性変化部位情報を入力として、空間解像度を上げるスパース性を導入して、最終的な光学特性変化部位情報を取得する。
【００９２】
以下、本実施の形態における物体観測システムの実験について説明する。図４は、物体観測システムの実験（ファントム実験）の模式図である。ここでのファントム実験において、物体に対して、一の方向から近赤外光を送光し、当該一の方向と同一の方向から後方散乱光を受光している。また、本実験では、水槽ファントムを用いて物体観測装置２が実現する推定法を検証するためのデータを得た。
【００９３】
本実験において、ファントム水槽の底面から、送光プローブ１１、および受光プローブ１２によって近赤外光を送光、および受光する。図５は、送光プローブ１１、受光プローブ１２の位置と吸光体を置いた位置を説明する図である。図５の灰色の丸が送光プローブ１１を表し、白丸が受光プローブ１２を表す。プローブは間隔ｌで配置されている。以下、吸光体を置いた位置を説明する。図５の（ａ）はセンター（座標値（０，０））、（ｂ）は中央ポイント（−ｌ／２，０）、（ｃ）は一のソース（送光プローブ１１）の位置（−ｌ／２，−ｌ／２）、（ｄ）は一のディテクター（受光プローブ１２）の位置（−ｌ／２，＋ｌ／２）、（ｅ）は第二の吸光体の位置（−ｌ／２＋ｄ，０）である。ここで、２吸光体の水平距離をｄと表記した。
【００９４】
本実験において、大脳皮質の灰白質と似た光学特性になるように、水槽の中は１．４％のイントラリピッドとインクの混合溶液で満たされており、分光光度計と近赤外光イメージング装置の計測に基づいた、光学パラメータの推定値は以下である。吸収係数μａ＝０．０１９ｍｍ−１、換算散乱係数μ'ｓ＝１．１ｍｍ−１。
【００９５】
吸光度変化計測には近赤外光イメージング装置（ＦＯＩＲＥ−３０００，株式会社島津製作所）を用いた。今回は８０５ｎｍの波長のみを計測に用いた。近赤外光は水槽の底面「ｚ＝０」に正方格子状に配置された送光プローブ１１から時分割で照射され、受光プローブ１２で受光する。サンプリングレートは１３０ｍｓであり、各々２０秒間計測を行った。プローブ間隔は１３ｍｍ、１８．４ｍｍ（＝１３ｍｍ×√２）、２６ｍｍ（＝１３ｍｍ×２）の３条件で実験を行った。
【００９６】
また、実験では、直径５ｍｍの球状の吸光体を１個ないし２個を溶液の中に配置して実験を行った。１吸光体実験では、（ａ）センター、（ｂ）送光プローブ１１と受光プローブ１２の中央ポイント、（ｃ）送光プローブ１１の真下、（ｄ）受光プローブ１２の真下の４か所で計測を行った（図５参照）。２吸光体実験では、（ｂ）と、（ｂ）からｄだけ距離を離した（ｅ）に吸光体を配置して計測を行った。
（実験１）
【００９７】
ここで、実際のファントム実験観測データを用いた３次元再構成法について紹介する。実験１において、データの計測条件は、ｚ＝０の面に置かれた送光プローブ１１と受光プローブ１２によって計測を行う（図６参照）。実験１におけるプローブの間隔は１８．４ｍｍである。そして、直径５ｍｍの球の吸光体を二つ入れることにより、吸光度変化を生じさせる。二つの吸光体の中心座標（ｘ，ｙ，ｚ）は（−９．２ｍｍ，０ｍｍ，１０ｍｍ）、および（３．３ｍｍ，０ｍｍ，１５ｍｍ）であり、水平に１２．５ｍｍ、垂直に５ｍｍ離しているが、これはプローブ間隔よりも短い。ここで、図６の座標系は図４の座標系を上下逆さまにしたものに対応する。これは脳機能計測において通常プローブ平面は上にあるため、それに合わせたものである。
【００９８】
図７がＤを単位行列としたときのＴｉｋｈｏｎｏｖ正則化法の推定結果である。なお、Ｄを単位行列とすることの意味は、感度の減衰に応じた正則化項の補正を行わないことである。また、図７に示すように、感度の深さによる感度の指数的な減衰により、実際の位置よりもプローブ平面へ近づいた解が得られる。また、図７に示すように、二つの吸光度変化の分離もできていない。
【００９９】
図８は、正規化項を感度で規格化したときの推定結果である。図８において、逆に深い方向にバイアスされていると見えるが、これは意味のある結果である。浅くて弱い活動と深くて強い活動は似た観測値を返すが、Ｔｉｋｈｏｎｏｖ正則化法では解を局在化できないために、これらの重ね合わせ（共存）状態が推定結果として得られる。深くて強い活動の方が目立つために、深い方向に推定されたと見えるだけであると考えられる。
【０１００】
図９は、図８の結果を初期値「λ０ｉ＝（１０・Ｘ＾Rｉ）-２」として階層ベイズ法に取り入れて推定したものである。つまり、図９は、物体観測装置２による推定結果（光学特性変化部位情報の取得結果）である。ここでＴｉｋｈｏｎｏｖ正則化法の推定結果は過度にスムーズ化されているため、その推定値を１０倍したものを初期値とした（以下の実験においても同様である）。推定部２４のスパース性適用手段２４４は、階層ベイズアルゴリズムを用いて、図８で得られるさまざまな解の重ね合わせ状態の中から適切なものを選び出す。このことにより、スパース性適用手段２４４は、より正確な推定結果を得る。本実験１において、二つの吸光体位置をそれぞれ識別でき、そのずれは１ボクセル以内であった。
（実験２）
【０１０１】
実験２において、吸光体を、深さを２．５ｍｍ刻みで変えながらセンターの位置（図５の（ａ））に置き、正しく深さが推定できるのかを検証した。プローブ間隔は２６ｍｍ、１８．４ｍｍ、１３ｍｍの三条件で計測を行った。
【０１０２】
図１０は、プローブ間隔と吸光体の深さ位置を変えながら、物体観測装置２が行った推定の結果（出力部２５の光学特性変化部位情報の出力例）を表している。推定結果は、３次元推定図を２次元で表すために、図１１で示すように、ｙ方向において吸光体を含む中心付近の３ボクセルの最大値のみをそれぞれのｘ−ｚ平面に投影している。図１０において、上の行からプローブ間隔がｌ＝２６ｍｍ、１８．４ｍｍ、１３ｍｍでの推定結果であり、左の列から、吸光体の中心の真の深さがｚ＝７．５ｍｍ、１０ｍｍ、１２．５ｍｍ、１５ｍｍ、１７．５ｍｍ、２０ｍｍ、２２．５ｍｍのときの推定結果である。図１０において、白丸は真の吸光体の位置を表している。推定結果の位置誤差が、ｘ、ｙ、ｚの軸でそれぞれ１ボクセル以内（２．５ｍｍ以内）であり、かつ推定最大値が０．０２５ｍｍ−１以上であれば推定できた、とし、太い枠で囲んでいる。
【０１０３】
図１０の上の行（プローブ間隔２６ｍｍでの結果）において、吸光体の真の深さがｚ＝７．５ｍｍからｚ＝１５ｍｍまでは推定結果は同じ深さに現れている（図１２の左側のグラフを参照）。これは深さの推定が実質的にできていないことを表している。また、ｚ＝１７．５ｍｍ以降は推定できなかった。
【０１０４】
図１０の中間の行（プローブ間隔１８．４ｍｍでの結果）において、吸光体の真の深さがｚ＝７．５ｍｍからｚ＝２０ｍｍまで推定が可能であり、真の深さと推定された深さが１ボクセル以内で一致していた（図１２の中央のグラフを参照）。
【０１０５】
図１０の下の行（プローブ間隔１３ｍｍでの結果）において、吸光体の真の深さがｚ＝７．５ｍｍからｚ＝２２．５ｍｍまで推定が可能であり、真の深さと推定された深さが１ボクセル以内で一致していた（図１２の右側のグラフを参照）。このようにプローブ間隔が１８．４ｍｍや１３ｍｍといった高密度計測であれば吸光度変化部位の深さを正確に推定できることがわかる。
【０１０６】
次に、推定結果のプローブ相対位置依存性を調べるために、吸光体位置を図５の（ａ）センター、（ｂ）中央ポイント、（ｃ）送光プローブ１１の位置、（ｄ）受光プローブ１２の位置の４地点のいずれかに置き、深さを２．５ｍｍずつ深くしていったときに、推定結果の位置誤差がｘ、ｙ、ｚの軸でそれぞれ１ボクセル以内であり、推定最大値が０．０２５ｍｍ−１以上となる限界の深さを調べたのが図１３である。ただし、プローブ間隔２６ｍｍの場合は、深さは推定できないので、ｘ、ｙの値のみが合っていれば良いものとした。高密度計測であれば吸光体位置にそれほど依存せずに正確な推定が可能であることがわかる。また、プローブ間隔が短くなるにつれて、推定可能な深さが増すことが分かる。
（実験３）
【０１０７】
実験３において、物体観測装置２が採用する階層ベイズ法の空間解像度を調べるために、吸光体を２個に増やして実験を行った。それぞれの吸光体のｘ，ｙ座標は、（−ｌ／２，０）、（−ｌ／２＋ｄ，０）とし（図５（ｂ）と（ｅ））、中心間の水平距離をｄと表記する。
【０１０８】
実験３において、プローブ間隔については、１８．４ｍｍと１３ｍｍの二条件のみで実験を行った。２６ｍｍの場合は実験２により三次元再構成ができないことが明らかになったので以後用いない。
【０１０９】
図１４はプローブ間隔１８．４ｍｍでの、物体観測装置２の出力結果（推定結果）である。本実験において、水平間距離ｄと２吸光体の深さｚを２．５ｍｍ刻みで変えながら実験を行った。推定結果は、図１０と同様に、３次元推定図を２次元で表すために、図１１で示すように、ｙ方向において吸光体を含む中心付近の３ボクセルの最大値のみをそれぞれのｘ−ｚ平面に投影している。図１４において、上の行から水平距離がｄ＝１７．５ｍｍ、１５ｍｍ、１２．５ｍｍ、１０ｍｍでの推定結果であり、左の列から、吸光体中心の真の深さがｚ＝１２．５ｍｍ、１５ｍｍ、１７．５ｍｍ、２０ｍｍのときの推定結果である。
【０１１０】
図１４に示す「ｄ＝１７．５ｍｍ、ｚ＝１７．５ｍｍ」の場合のように、プローブ間隔と同程度、２吸光体が離れている場合、ｚ＝１７．５ｍｍの深さまで二つの吸光体は正しい位置に推定された。
【０１１１】
また、プローブ間隔よりも２吸光体の距離を短くしていった場合も推定が可能であった。特にプローブ間隔が１８．４ｍｍのときは、中心間距離がｄ＝１０ｍｍ（５ｍｍ球の端間距離５ｍｍ）と、プローブ間隔のおよそ半分の距離にある２吸光体を識別可能であった（図１４「ｄ＝１０ｍｍ、ｚ＝１５ｍｍ」の図を参照）。
【０１１２】
図１５は、プローブ間隔を１８．４ｍｍ、１３ｍｍの二条件、２吸光体水平距離をｄ＝１５ｍｍ、１２．５ｍｍ、１０ｍｍとしたときに、吸光体の深さｚを２．５ｍｍずつ深くしていったときに、二吸光体とも、推定結果の位置誤差がｘ、ｙ、ｚの軸でそれぞれ１ボクセル以内であり、推定吸光係数変化値が０．０２５ｍｍ−１以上となる限界の深さをまとめたものである。
（実験４）
【０１１３】
複雑な条件でも深さを含めた３次元推定が可能かどうかを確認するために、実験４を行った。実験４において、２つの吸光体の深さを５ｍｍだけずらして実験を行った。それぞれの吸光体のｘ，ｙ座標は、実験３と同様に、（−ｌ／２，０）、（−ｌ／２＋ｄ，０）とし（図５（ｂ）と（ｅ））、中心間の水平距離をｄと表記する。また、（ｂ）地点の吸光体は（ｅ）地点の吸光体よりも５ｍｍだけ浅く置くものとする。実験４において、プローブ間隔については、深さの推定が可能な１８．４ｍｍと１３ｍｍの二条件のみで実験を行った。
【０１１４】
図１６は、プローブ間隔１８．４ｍｍでの推定結果である。実験４において、水平間距離ｄと２吸光体の深さｚを２．５ｍｍ刻みで変えながら実験を行った。推定結果は、図１０と同様に、３次元推定図を２次元で表すために、ｙ方向において吸光体を含む中心付近の３ボクセルの最大値のみをそれぞれのｘ−ｚ平面に投影している。
【０１１５】
図１６の「ｄ＝１７．５ｍｍ、ｚ＝１５＆２０ｍｍ」の場合のように、プローブ間隔と同程度、２吸光体が離れている場合、１吸光体での推定限界深さｚ＝２０ｍｍに吸光体があるときでも、もう一つの浅い位置にある吸光体に邪魔されずに、正しい位置に推定されることが、図１６から分かる。
【０１１６】
また、プローブ間隔よりも２吸光体距離を短くしていった場合も推定は可能であった。このように深さの異なる複数の吸光体があっても、物体観測装置２が採用する推定法によれば、高い空間解像度が得られるものと考えられる。
【０１１７】
図１７は、プローブ間隔を１８．４ｍｍ、１３ｍｍの二条件、２吸光体水平距離をｄ＝１５ｍｍ、１２．５ｍｍ、１０ｍｍとしたときに、吸光体の深さｚを２．５ｍｍずつ深くしていったときに、二吸光体とも、推定結果の位置誤差がｘ、ｙ、ｚの軸でそれぞれ１ボクセル以内であり、推定吸光係数変化値が０．０２５ｍｍ−１以上となる限界の深さをまとめたものである。
【０１１８】
以上、本実施の形態によれば、物体内部の局所的な光学特性変化部位を三次元情報として高精度に推定できる。さらに具体的には、本物体観測システムが採用する、深さ方向を含めた高精度３次元再構成アルゴリズムにより、物体内部の局所的な光学特性変化部位を三次元情報として高精度に推定できる。なお、物体観測システムが採用する再構成アルゴリズムは以下の２ステップから成る。第一に、正則化項を感度で規格化したＴｉｋｈｏｎｏｖ正則化法で推定を行うことにより、空間的に広がった解の候補を求める。第二に、その推定結果を初期値とした階層ベイズ推定を行うことで、解の候補からよりふさわしいものを絞り込む。
【０１１９】
なお、本実施の形態における上記のファントム実験の結果、プローブ間隔が２６ｍｍでは観測情報の不足から再構成アルゴリズムはうまく働かなかったが、１８．４ｍｍと１３ｍｍのプローブ間隔においては、深さを正確に推定することができ、プローブ間隔よりも短い距離にある二つの吸光体を識別可能であった。特に、１８．４ｍｍの間隔では、そのおよそ半分の距離である中心間距離１０ｍｍの二つの吸光体を識別可能であった。
【０１２０】
また、本実施の形態によれば、物体観測装置２において、階層ベイズ法によりスパース性を導入した。しかし、本実施の形態において、Ｌ１ノルムを用いることにより、スパース性を導入することも可能である。例えば、数式２の正規化項の２乗を１乗にすると、感度の規格化とスパース性の導入を同時に行うことができる。
【０１２１】
なお、Ｌ１正則化を用いた推定も実際に行った。Ｌ１正則化でも解がスパースになるため、本実験において、階層ベイズ法を用いた場合の条件と同じ条件で平滑化フィルタＷを導入した。実験データに適用した結果、適切な正則化パラメータλを選べば、階層ベイズ法による場合とほぼ同程度の推定精度が得られることがわかった。なお、かかる実験において、上述した感度による補正とスパース性の適用を同時に行うことができることを示せた。つまり、推定部は、光学特性変化と光強度変化との関係を示す情報である感度情報を格納しており、当該感度情報を用いて、深さによる感度の減衰に応じた補正を行い、かつ空間解像度を上げるスパース性を導入して光学特性変化部位情報を取得すれば良く、当該感度による補正とスパース性の適用とを同時に行っても、順次に行っても良い。
【０１２２】
さらに、本実施の形態における処理は、ソフトウェアで実現しても良い。そして、このソフトウェアをソフトウェアダウンロード等により配布しても良い。また、このソフトウェアをＣＤ−ＲＯＭなどの記録媒体に記録して流布しても良い。なお、このことは、本明細書における他の実施の形態においても該当する。なお、本実施の形態における物体観測装置２を実現するソフトウェアは、以下のようなプログラムである。つまり、このプログラムは、記録媒体に、各受光プローブの光強度に関する情報であり、基準となる情報である基準光強度情報を格納しており、コンピュータを、観測対象の物体に対して光を送光する送光プローブと、光を受光する受光プローブの組であるプローブセットを１組以上有する近赤外光計測装置を用いて、前記各プローブセットを構成する各送光プローブから、一の光強度の光を前記物体に対して発した場合に、各受光プローブが受光する光強度に関する情報である光強度情報を取得する光強度情報取得部と、前記記録媒体の基準光強度情報と、前記光強度情報取得部が取得した光強度情報とから、プローブセットごとに、光強度の変化に関する情報である光強度変化情報を取得する光強度変化情報取得部と、前記プローブセットごとの光強度変化情報を用いて、前記物体の内部の吸光度変化の位置に関する三次元の情報である光学特性変化部位情報を取得する推定部と、前記推定部が取得した光学特性変化部位情報を出力する出力部として機能させるプログラムであって、前記推定部は、光学特性変化と光強度変化との関係を示す情報である感度情報を用いて、深さによる感度の減衰に応じた補正を行う補正手段と、空間解像度を上げるスパース性を導入して前記光学特性変化部位情報を取得するスパース性適用手段とを具備するものとして、コンピュータを機能させるプログラムである。
【０１２３】
また、上記プログラムにおいて、前記推定部は、前記光学特性変化部位情報を取得するために使用する演算式であり、逆問題の解を算出する演算式であり、前記感度情報に応じた補正を行なえるコスト関数である第一の演算式を格納し得る演算式格納手段をさらに具備し、前記補正手段は、前記感度情報を前記第一の演算式に代入し、第一の光学特性変化部位情報を取得し、前記スパース性適用手段は、前記第一の光学特性変化部位情報を用いて、空間解像度を上げるスパース性を用いて最終的な光学特性変化部位情報を取得し、前記出力部は、前記スパース性適用手段が取得した光学特性変化部位情報を出力するものとして、コンピュータを機能させることは好適である。
【０１２４】
また、上記プログラムにおいて、前記演算式格納手段は、初期値に対して階層ベイズ推定を行う第二の演算式をも格納しており、前記スパース性適用手段は、前記第一の光学特性変化部位情報を初期値として前記第二の演算式に代入し、当該第二の演算式を実行することにより階層ベイズ推定を行い、最終的な光学特性変化部位情報を取得するものとして、コンピュータを機能させることは好適である。
【０１２５】
また、上記プログラムにおいて、前記観測対象の物体は、生物の脳であり、前記基準光強度情報は、安静状態において測定された各受光プローブの光強度に関する情報であり、生物の脳機能を観測するものとして、コンピュータを機能させることは好適である。
【０１２６】
また、上記プログラムにおいて、前記観測対象の物体は、乳房であり、前記基準光強度情報は、正常な細胞の乳房に対する各受光プローブの光強度に関する情報であり、乳がんの検診を行えるものとして、コンピュータを機能させることは好適である。
【０１２７】
また、上記プログラムにおいて、前記各送受光プローブの間隔は、２０ｍｍ以下であるものとして、コンピュータを機能させることは好適である。
【０１２８】
また、図１８は、本明細書で述べたプログラムを実行して、上述した実施の形態の物体観測装置２等を実現するコンピュータの外観を示す。上述の実施の形態は、コンピュータハードウェア及びその上で実行されるコンピュータプログラムで実現され得る。図１８は、このコンピュータシステムの概観図であり、図１９は、コンピュータシステムのブロック図である。
【０１２９】
図１８において、コンピュータシステム３４０は、ＦＤドライブ、ＣＤ−ＲＯＭドライブを含むコンピュータ３４１と、キーボード３４２と、マウス３４３と、モニタ３４４とを含む。
【０１３０】
図１９において、コンピュータ３４１は、ＦＤドライブ３４１１、ＣＤ−ＲＯＭドライブ３４１２に加えて、ＣＰＵ３４１３と、ＣＤ−ＲＯＭドライブ３４１２及びＦＤドライブ３４１１に接続されたバス３４１４と、ブートアッププログラム等のプログラムを記憶するためのＲＯＭ３４１５とに接続され、アプリケーションプログラムの命令を一時的に記憶するとともに一時記憶空間を提供するためのＲＡＭ３４１６と、アプリケーションプログラム、システムプログラム、及びデータを記憶するためのハードディスク３４１７とを含む。ここでは、図示しないが、コンピュータ３４１は、さらに、ＬＡＮへの接続を提供するネットワークカードを含んでも良い。
【０１３１】
コンピュータシステム３４０に、上述した実施の形態の物体観測装置等の機能を実行させるプログラムは、ＣＤ−ＲＯＭ３５０１、またはＦＤ３５０２に記憶されて、ＣＤ−ＲＯＭドライブ３４１２またはＦＤドライブ３４１１に挿入され、さらにハードディスク３４１７に転送されても良い。これに代えて、プログラムは、図示しないネットワークを介してコンピュータ３４１に送信され、ハードディスク３４１７に記憶されても良い。プログラムは実行の際にＲＡＭ３４１６にロードされる。プログラムは、ＣＤ−ＲＯＭ３５０１、ＦＤ３５０２またはネットワークから直接、ロードされても良い。
【０１３２】
プログラムは、コンピュータ３４１に、上述した実施の形態の物体観測装置等の機能を実行させるオペレーティングシステム（ＯＳ）、またはサードパーティープログラム等は、必ずしも含まなくても良い。プログラムは、制御された態様で適切な機能（モジュール）を呼び出し、所望の結果が得られるようにする命令の部分のみを含んでいれば良い。コンピュータシステム３４０がどのように動作するかは周知であり、詳細な説明は省略する。
【０１３３】
また、上記プログラムを実行するコンピュータは、単数であってもよく、複数であってもよい。すなわち、集中処理を行ってもよく、あるいは分散処理を行ってもよい。
【０１３４】
また、上記各実施の形態において、各処理（各機能）は、単一の装置（システム）によって集中処理されることによって実現されてもよく、あるいは、複数の装置によって分散処理されることによって実現されてもよい。
【０１３５】
本発明は、以上の実施の形態に限定されることなく、種々の変更が可能であり、それらも本発明の範囲内に包含されるものであることは言うまでもない。
【産業上の利用可能性】
【０１３６】
以上のように、本発明にかかる物体観測システムは、物体内部の局所的な光学特性変化部位を三次元情報として高精度に推定できるという効果を有し、脳機能観測システム等として有用である。
【図１】
0005570674.tif 000027
【図２】
0005570674.tif 000028
【図３】
0005570674.tif 000029
【図４】
0005570674.tif 000030
【図１３】
0005570674.tif 000031
【図１５】
0005570674.tif 000032
【図１７】
0005570674.tif 000033
【図１８】
0005570674.tif 000034
【図１９】
0005570674.tif 000035
【図５】
0005570674.tif 000036
【図６】
0005570674.tif 000037
【図７】
0005570674.tif 000038
【図８】
0005570674.tif 000039
【図９】
0005570674.tif 000040
【図１０】
0005570674.tif 000041
【図１１】
0005570674.tif 000042
【図１２】
0005570674.tif 000043
【図１４】
0005570674.tif 000044
【図１６】
0005570674.tif 000045