JP2002333400A - 吸収物質濃度空間分布画像化装置 - Google Patents

Abstract

(57)【要約】 【課題】光を用いて散乱体中の吸収物質濃度の空間分布
を画像化する方法を提供する。 【解決手段】計測された検出光強度から計測値行列Ｙを
生成するｓｔｅｐ１、散乱体内部の吸収物質の吸光係数
と、吸収物質の存在しない散乱体を模擬して得る散乱体
内における光が通過する光路データをモンテカルロ法に
より求め、光路データから、平均光路長の空間分布行列
Ａを生成するｓｔｅｐ２、行列Ｙと行列Ａから未知数行
列Ｘを導出するｓｔｅｐ３、求められた未知数行列Ｘか
ら被検体内部の吸収物質濃度を求めこの空間分布を表示
するｓｔｅｐ４からなる画像化する方法である。モンテ
カルロ法による光路の計算が、ｓｔｅｐ２の処理開始以
前に既になされ、光路データが記憶保存されている場合
には、記憶保存されている光路データを使用する。高分
解能に精度良く、高速に吸収物質濃度の３次元分布を画
像化できる。

Description

【発明の詳細な説明】

【０００１】

【発明の属する技術分野】本発明は、散乱物質中に吸収
物質が存在してなる被検体において、吸収物質濃度の空
間分布を求める方法に関する。特に、光を用いて生体内
部あるいは大気中の吸収物質濃度の空間分布を計測する
方法、装置に関する。

【０００２】

【従来の技術】現在、医療用の画像診断装置として、Ｍ
ＲＩ（核磁気共鳴イメージング装置）Ｘ線ＣＴ、超音波
診断装置等がある。それぞれ、臨床医療の現場で多用さ
れ有用性が認められている。現在、臨床医療に使用され
ている画像診断装置は、主に生体内の組織、構造を断層
像として得るものが多く、生体内に分布する化学物質の
代謝活動に起因する生体機能を計測するための簡便で実
用的な画像診断装置は少ない。生体機能の計測により得
られる情報は、脳疾患やガンの予防診断、脳機能の解明
に有用と期待されており、簡便、実用的な生体機能画像
診断装置の開発が医療の分野では待たれている。散乱物
質中に、光の波長により吸収特性が異なる吸収物質が存
在してなる被検体中の吸収物質の濃度の空間分布を画像
化する光計測装置として、生体内の特定の色素濃度の空
間分布を画像化する生体光計測装置がある。生体中に
は、生体機能に関与する吸収物質が多く存在し、吸収物
質の濃度の空間分布を計測できる生体光計測装置は、生
体機能を画像化して診断に有効となる装置の一つとして
開発が急がれている。

【０００３】生体機能の中でも、生体内酸素飽和度は、
生体中に多く含まれる特定色素の濃度に対応しており、
この色素濃度は可視から近赤外領域における光吸収の量
から計測することが可能である。特定色素として、酸化
型あるいは還元型の別により、異なる光吸収の分光特性
を持つことが知られている、血液中に存在するヘモグロ
ビン（以下では、Ｈｂと略記する）、細胞中に存在する
チトクロームａａ３、筋肉中に存在するミオグロビン等
がある。さらに、光による計測は、光ファイバーや小型
光源を用いることにより、簡便、小型な装置として開発
される可能性を持っている。

【０００４】可視から近赤外の光を用いた生体機能を計
測する装置が、例えば、特開昭５７−１１５２３２号公
報、特開昭６３−２７５３２３号公報に記載されてい
る。生体光計測装置で重要な課題の１つとして、生体内
に存在する吸収物質の空間分布を画像化する方法の確立
が挙げられる。この課題は吸収物質の濃度と位置を求め
る問題である。従来技術での画像化方法として、（１）
Ｘ線ＣＴのアルゴリズムを用いる方法、（２）輸送方程
式を用いる方法、（３）モンテカルロシミュレーション
の結果を用いる方法、に大別される。

【０００５】（１）の方法の具体的例として、時間ゲー
ト法を用いて得た計測データから、Ｘ線ＣＴにおいて使
用されてきた画像再構成のアルゴリズムにより画像を得
る方法（特開平０１−２０９３４２号公報に記載）が挙
げられる。（２）の方法の具体的例として、アリッジ(
Aridge )等の方法（" Reconstruction methods for inf
rared absorption imaging ", Proc. SPIE 1431, 204-2
15(1991)に記載）や、シンガー( Singer )等の方法（"
Image reconstruction ofthe interior of bodies that
diffuse radiation ", Science, 248, 990-993(1990)
に記載）が挙げられる。両者の方法共に、散乱と吸収現
象を表わす数学モデルに基づき、照射光強度と散乱を表
わす変数（以下では、散乱変数と略記する）と吸収物質
の吸収を表わす変数（以下では、吸収変数と略記する）
とから、輸送方程式を立て検出光強度を計算する。ここ
で、散乱変数と吸収変数は散乱体内部の位置の関数であ
り、また、吸収変数は吸収物質の濃度を反映する値であ
る。上記の輸送方程式の散乱変数と吸収変数の値を散乱
体内部の任意の位置で独立に変化させ、計算された出光
強度と実際に計測された検出光強度との差を評価し、そ
の差が最小となるまで計算を繰り返す。この繰り返しの
結果得られた、吸収変数の空間分布が吸収物質濃度の空
間分布を表わしており、その分布を画像化して表示す
る。アリッジ( Aridge )等の方法では、散乱変数として
散乱係数を、吸収変数として吸収係数を用いており、シ
ンガー( Singer )等の方法では、散乱変数として散乱体
内に設定された面積要素もしくは体積要素から出て行く
光子の確率を、吸収変数として要素内で吸収される光子
の確率を用いている。

【０００６】（３）の方法の具体的例として、バーバー
ランダル エル等によるモンテカルロ法を利用した画
像化方法がある（特表平３−５０５９２２号公報に記
載）。この方法の概要は以下の通りである。媒体中の全
ての点（位置及び方向）に対し、フラックスによって与
えられる重みを割当て、その点で検出器の吸収体の存在
しない媒体に対する応答に対する光子の期待される寄与
を割り当てることによって、複雑な媒体を画像化する技
術に関する。照射される輻射線の測定された強度が、予
定される強度と比較されることにより、対応する散乱体
積の透過係数が決定される。予定される強度は、予め規
定されたモデル媒体に対するモンテカルロシュミレーシ
ョンを用いた計算により得られ、透過強度係数が重み関
数に適用される。媒体中における微小な吸収の検出器応
答特性に対する影響は、影響を受ける体積要素に対する
吸収散乱断面積と対応する重み関数との積となるものと
仮定される。重み関数は、モンテカルロ法を用いて計算
され、領域内の全ての点の、与えられた光源−検出器の
配置に対する検出器応答性への相対寄与の３次元マップ
を与える。画像化は、透過係数と、光源−検出器の各配
置ごとに、各体積要素に対応する重み関数との積を計算
することにより達成される。これらの積の総和は、光源
−検出器の配置のすべてに対するすべての体積要素の重
ね合わせに対応する。この計算は、全ての光源の位置の
関数として、全ての領域にわたる相対吸収断面積の３次
元マップを与える。

【０００７】

【発明が解決しようとする課題】生体光計測装置は、Ｈ
ｂ（ヘモグロビン）のように、生体内で機能する吸収物
質濃度の空間分布を、画像として表示して、例えば癌の
様な疾病の予防診断への応用が期待されている。しか
し、生体は光に対して強い散乱特性を有しているため
（散乱係数は１．０〔ｍｍ^-１〕程度である）、生体内
部の吸収物質濃度の空間分布の精度良い画像化は非常に
困難である。生体光計測装置を臨床医療で実用化するた
めには、生体の様な散乱体の内部に存在する吸収物質濃
度の空間分布を精度良く画像化する方法が重要である。
現在、このような生体内部の吸収物質濃度の画像化方法
として、実用化レベルに達していると考えられる方法は
まだ存在しない。従来技術における方法では、様々な問
題があり、上記（１）の方法では、３次元的な散乱現象
を２次元空間上で扱うため再構成された画像の空間分解
能が低い。上記（２）の方法では、繰り返し計算をする
ため長時間の演算が必要で、現状の計算機速度では高分
解能の画像を得ることは実質的に困難であり、また局所
最小値（ local minimun ）に落ちこみ真の最適解が得
られにくい場合がある。上記（３）の方法では、不透明
な媒体中にある標的物体を３次元的に画像化することの
記載はあるが、画像化する具体的な手順は数式的に明確
には示されておらず、さらに媒体中の吸収物質の濃度を
計算することの記載はない。

【０００８】本願発明の目的は、３次元的な散乱現象を
考慮して、特に生体内で機能する吸収物質濃度を精度良
く、短時間に空間分解能で求め、吸収物質の濃度の３次
元空間分布を求める方法、及び装置を提供することにあ
る。

【０００９】

【課題を解決するための手段】本願発明は、被検体を通
過して検出された光強度が、被検体中に存在する吸収物
質の濃度と、散乱に基づく光子の平均光路長との線形結
合により表わされると仮定し、平均光路長の空間分布を
使用して、被検体中に存在する吸収物質の濃度の３次元
分布を求めることに特徴を有する。即ち、光吸収物質を
含んでなる光散乱体からなる被検体に、所定の光照射位
置から所定の波長を有するパルス光又は連続光を照射す
る照射ステップと、被検体を通過する光の強度を所定の
光検出位置で検出する検出ステップと、を有し被検体内
に分布する吸収物質濃度の空間分布画像化方法におい
て、光照射位置から光検出位置までの複数の光路を求め
るステップと、複数の光路から平均光路長を求めるステ
ップと、平均光路長と、被検体に照射される光の強度
と、光検出位置で検出される検出光強度と、被検体に照
射される光に対する吸収物質の光学定数とから、各体積
要素毎の吸収物質の濃度を求めることを特徴とする。平
均光路長の空間分布は、モンテカルロ法により模擬して
計算により求める。この計算では、被検体と実質的に同
一の散乱係数を有し、仮想被検体を包含する大きさの仮
想散乱体内の一点から発する複数光子のそれぞれのたど
る複数光路を求める。仮想被検体は被検体と同一又は類
似した形状を有し吸収物質を含まないものとする。

【００１０】より詳細に説明すると、本願発明の吸収物
質濃度の空間分布画像化方法では、被検体と同一又は類
似した形状を有し、吸収物質を含まない仮想被検体を想
定して、仮想被検体において設定される光照射位置と光
検出位置の複数の対のそれぞれで、光子がたどる光照射
位置から光検出位置に至るまでの複数の光路を求めるス
テップと、複数の光路から光照射位置から光検出位置ま
での仮想被検体内における平均光路長を求めるステップ
と、平均光路長と、被検体に照射される照射光強度と、
光検出位置で検出される検出光強度と、被検体に照射さ
れる光に対する吸収物質の光学定数（吸光係数）と、か
ら吸収物質の濃度の空間分布を求めるステップと、吸収
物質の濃度の空間分布を表示するステップと、を有する
ことに特徴を有する。被検体の形状を表す空間は複数の
体積要素に分割され、各体積要素ごとにに平均光路長が
求められ、各体積要素内の吸収物質濃度が一様であると
仮定し、各体積要素毎の吸収物質濃度が求められる。

【００１１】より具体的な吸収物質濃度の空間分布画像
化のステップは、上記対における照射光強度と検出光強
度との比に基づいて計測値行列を生成するステップと、
上記対における平均光路長に基づいて平均光路長の空間
分布行列を生成するステップと、計測値行列と平均光路
長の空間分布行列から吸収物質濃度を求めるステップ、
とからなる。予め与えられている所定の形状と所定の散
乱係数を有する被検体に本発明の手順を適用する場合
は、平均光路長の空間分布行列を生成するステップで求
められ記憶手段に予め記憶されている平均光路長の空間
分布行列を使用する。

【００１２】光路を求めるステップは、被検体と実質的
に同一の散乱係数を有し、仮想被検体を包含する大きさ
の仮想散乱体内の一点から発する光子のたどる光路を模
擬するモンテカルロ法により模擬して求められ、モンテ
カルロ法では、光子が順次散乱される複数の散乱点の内
の連続する２つの散乱点のそれぞれの空間座標位置を関
連づける散乱角度特性を表す関数、及び乱数を使用す
る。モンテカルロ法により求められた複数の光路のそれ
ぞれを表す複数の散乱点の空間座標位置が記憶され、複
数の光路データが記憶手段に記憶される。光路データ
は、仮想散乱体内の一点から発する単一光子のたどる有
限光路長を有する光路データを求めるか、あるいは仮想
被検体の２点を結ぶ最大距離以上の半径を有する球を仮
想散乱体として、球の中心から複数光子を発して、複数
光子のそれぞれの光子のたどる光路のうち、球の内部に
存在する光路データを求める。

【００１３】仮想被検体での光照射位置から光検出位置
までの光路は、仮想被検体での光照射位置と仮想散乱体
内の１点とを一致させ、光検出位置とモンテカルロ法に
より求められた光路とが交差する点を、光路を表す散乱
点の空間座標、又は仮想被検体の形状を表す空間座標の
いずれか一方の空間座標を、他方の空間座標に対して座
標変換（回転移動及び／又は並進移動）して求める。求
められた仮想被検体での光照射位置から仮想被検体での
光検出位置までの複数光路のそれぞれのを、被検体での
対応する光照射位置から光検出位置までの光路とする。
なお、仮想被検体の形状を表す空間は、被検体の形状を
表す空間を複数の体積要素に分割したと同様に体積要素
に分割される。また、光照射位置と光検出位置の対の数
は、被検体、仮想被検体を分割した体積要素の総数以上
の値に設定される。上記の複数光路データから仮想被検
体の各体積要素を通過する平均光路長を体積要素ごとに
求め、被検体の各体積要素を通過する平均光路長とす
る。光照射位置（ｉ）と光検出位置（ｊ）の対（ｉ−
ｊ）に対して得られた光路に関するデータに基づいて、
光照射位置と光検出位置の位置関係が入れ替わった光照
射位置（ｊ）と光検出位置（ｉ）の対（ｊ−ｉ）に対す
る光路に関するデータを、対（ｉ−ｊ）に対して得られ
た光路に関するデータから、座標変換により求めること
ができる。

【００１４】照射ステップにおいて被検体に照射される
光は、４００ｎｍから２０００ｎｍの範囲にある複数の
波長を有する光であり、複数の波長の差が１００ｎｍ以
下とする。被検体を通過する光の強度が時間分解計測さ
れる場合には、検出光強度は時間分解計測された光強度
の所定の時間内での光強度の積分値から求められ、所定
の時間内に光照射位置から光検出位置に到達する光子に
関する光路から求められる平均光路長に基づいて平均光
路長の空間分布行列が生成される。

【００１５】本発明では、計測された検出光強度から計
測値行列Ｙを生成し、散乱体内部の吸収物質の吸光係数
と、吸収物質の存在しない散乱体を模擬して得る散乱体
内における光が通過する光路データをモンテカルロ法に
より求め、光路データから、平均光路長の空間分布行列
Ａを生成し、行列Ｙと行列Ａとから、未知数行列Ｘを導
出して、求められた未知数行列Ｘから被検体内部の吸収
物質濃度を求めこの空間分布を表示する。さらに、異な
る複数位置に配置される光照射位置と光検出位置の対の
数は、被検体、仮想被検体を分割した体積要素の総数以
上の値に設定するので、ある体積要素に注目するとき、
複数方向からの散乱の影響が考慮されるので、従来方法
に比較して、約３倍以上の空間分解能を有する３次元的
な吸収物質の空間分布を得ることができ、１０〜１００
倍程度計算時間が短縮され、さらに定量性の改善が可能
である。本方法では、局所の最適解に落ち込むことなく
解を得ることができる。さらに、モンテカルロ法による
光路の計算が、予め計算された光路データが記憶保存さ
れている場合には、記憶保存されている光路データを使
用するので、従来のモンテカルロ法を適用する方法と比
較して、平均光路長の空間分布を約１００倍以上の速度
で求めることができる。

【００１６】

【発明の実施の形態】従来技術として、光を用いて生体
内部に分布する特定色素濃度の空間分布を計測する生体
光計測装置があり、また大気中に存在する吸収物質濃度
の空間分布の計測として、レーザーレーダがある。この
様に、散乱物質中に存在する吸収物質の濃度の空間分布
を求め画像化する方法を必要とする分野は広い。特に以
下では、生体光計測装置を代表的な例として説明をす
る。

【００１７】本願発明は、光計測装置で計測された検出
光強度から、被検体内部に存在する吸収物質の濃度を画
像化する方法に係り、本願発明における手順は要約する
と、計測された検出光強度から計測値行列Ｙを生成する
ｓｔｅｐ１、散乱体内部の吸収物質の吸光係数と、吸収
物質の存在しない散乱体を模擬して得る散乱体内におけ
る光が通過する光路データをモンテカルロ法により求
め、光路データから、平均光路長の空間分布行列Ａを生
成するｓｔｅｐ２、行列Ｙと行列Ａとから未知数行列Ｘ
を導出するｓｔｅｐ３、求められた未知数行列Ｘから被
検体内部の吸収物質濃度を求めこの空間分布を表示する
ｓｔｅｐ４からなる。モンテカルロ法による光路の計算
が、ｓｔｅｐ２の処理開始以前に既になされ、光路デー
タが記憶保存されている場合には、記憶保存されている
光路データを使用する。

【００１８】まず最初に、被検体によって照射光が散乱
されない場合（即ち、被検体が非散乱体である場合）に
関して説明し、次いで、被検体によって照射光が散乱さ
れる場合（即ち、被検体が散乱体である場合）に関して
説明する。図１に示す様な、光源１から被検体２へ所定
の波長の光を照射し、光路３を通過してきた光の強度を
検出器４で検出する光計測系を考える。図１は、光路３
を含む被検体２の断面を示す。被検体２中に１成分の吸
収物質が一様濃度ｃで存在する場合に、光源１からの照
射光強度をＩ₀、検出器４で検出される検出光強度を
Ｉ_d、吸収物質の吸光係数をε、光照射位置と光検出位
置間の幾何学的距離で決まる光路長をＬ’とすると、
（数１）で定義されるベール・ランバート則が成り立
つ。 Ｉ_d／Ｉ₀＝ｅｘｐ（−εｃＬ’） …（数１） 図１と同様の光計測系を考え、図２に示すように、被検
体２が、吸収物質濃度の異なる２個の体積要素５−１、
５−２に分けられる場合に、体積要素５−１中での吸収
物質濃度をｃ₁、光路長をＬ₁’、体積要素５−２中での
吸収物質濃度をｃ₂、光路長をＬ₂’とし、光源１からの
照射光強度をＩ₀、体積要素５−１から５−２に入射す
る光の強度をＩ_C、検出器４で検出される検出光強度Ｉ_d
とすると、（数１）より、体積要素５−１、体積要素５
−２に関し、それぞれ（数２）、（数３）が成り立つ。 Ｉ_C／Ｉ₀＝ｅｘｐ（−εｃ₁Ｌ₁’） …（数２） Ｉ_d／Ｉ_C＝ｅｘｐ（−εｃ₂Ｌ₂’） …（数３） 体積要素５−１から５−２に入射する光の強度Ｉ_Cは計
測されないので、（数２）、（数３）を用いてＩ_Cを消
去すると、（数４）が成り立つ。 Ｉ_d／Ｉ₀＝ｅｘｐ（−ε（ｃ₁Ｌ₁’＋ｃ₂Ｌ₂’）） …（数４） 次に、被検体が散乱体である場合について記述する（以
下の説明では、特に明記しない場合は、被検体は散乱体
とする）。

【００１９】図３に示すような光計測系で、被検体２中
に１成分の吸収物質が一様な濃度ｃで存在する場合に、
光源１からの照射光強度をＩ₀、検出器４で検出される
検出光強度をＩ_d、吸収物質の吸光係数をε、被検体内
部に吸収物質が存在しないときの光路長（以下光路長と
略す）をＬ、散乱によって被検体２の外に出て行く成分
を減衰定数Ｄ_Sとすると、（数５）が成り立つ（ Victor
Twersky,"Absorptionand Multiple Scattering by Bio
logical Suspensions ", J. Opt. Soc. Am.,60, 1089(1
970)に記載）。 Ｉ_d／Ｉ₀＝Ｄ_S・ｅｘｐ（−εｃＬ） …（数５） ここで、光を光子として扱うと、光源１から照射された
光子は被検体２中で散乱を繰り返して検出器４に到達す
る。光子は散乱の度に進行方向が変わるので、その光路
は、図３中に示す光路３の様にジグザグを描く。なお、
図３では、図１、図２と同様に被検体２は立体であり、
簡単のために、２本の光路３が１つの断面内にある場合
を例示しており、一般に、光路３は３次元的にジグザグ
を描く。散乱による進行方向の変化はランダムであり光
路３は各光子によって異なり、（数５）中の光路長Ｌは
複数個存在し一義的に決定できない。そこで、複数個存
在する光路長の平均値（以下では、平均光路長と略記す
る）を近似値として用いる。検出器に到達する光子数を
ｍ本とすると、光路の数はｍ本存在し、ｍ本ある光路の
うちｋ番目の光路長をＬ_Kとして、平均光路長Ｌ_Aを（数
６）で定義し、さらに、平均光路長Ｌ_Aを用いて（数
５）を書き換えると（数７）になる。 Ｌ_A＝（ΣＬ_K）／ｍ …（数６） （数６）において、加算Σはｋ＝１〜ｍについて行な
う。 Ｉ_d／Ｉ₀＝Ｄ_S・ｅｘｐ（−εｃＬ_A） …（数７） 図３に示す光計測系を、図４に示すように、被検体２
を、ｎ個の体積要素に分割して考え、光源から発した光
子は、全ての体積要素を通過して検出器４に至るものと
仮定する。さらに各体積要素における吸収物質濃度は一
般に異なると考える。図４は、簡単のために、図３で例
示したのと同様に、光源１から照射された光子が被検体
２中で散乱を繰り返して検出器４に到達する光路が全
て、厚さが１体積要素からなる１断面内にある場合を示
し、各体積要素の中心部分に示される太い線分は、この
１断面内に存在する複数個の光路を各体積要素に分割し
て考え、各体積要素において複数個の光路のそれぞれ進
行方向での距離を加算し一定値で割算して得た光路長を
示し（得られる光路長の最大長を体積要素の辺長さ以下
とするために一定値で割算した）、光源１と検出器４を
結ぶ方向に平行に示したものである。言いかえると、図
４は、被検体２として２次元物体を想定し、被検体２を
空間（面積）要素に分割して考え、光源１からの複数光
子のそれぞれが、２次元平面内で散乱の度に進行方向を
変えて光路を形成し進む場合に、各面積要素において複
数光子の描くそれぞれの光路を、それぞれの光子の進行
方向での距離を加算し一定値で割算して得た光路長とし
て太い線分により、光源１と検出器４を結ぶ方向に平行
に示したものである。吸収物質が存在せず、かつ一様な
散乱体では、図４に示すように、この太い線分の長さの
分布は光源１と検出器４を結ぶ直線に対して対称とな
り、この直線から距離が離れる程太い線分の長さは短く
なる。しかし、一般に被検体２は一様な散乱体ではなく
しかも、被検体２の各部位は異なる濃度の吸収物質が存
在するため、このような場合には図４に示す太い線分の
長さは、各要素で異なることになる。

【００２０】本願発明では、光源１からの複数光子のそ
れぞれは、被検体２を構成する全ての体積要素を一般に
通過していくと想定し、各体積要素が図３に示す被検体
であると見做して（数７）を適用し、さらに複数の異な
る濃度の吸収物質が存在する場合に適用される（数４）
を適用する。即ち、光源１からの複数光子のそれぞれ
は、散乱のためジグザグを描いて進行するが、このとき
各光子は、全ての体積要素を通過して検出器４に至るも
のと仮定すると、（数４）、（数７）を考慮して（数
８）が成立するものとする。但し、（数８）において、
各体積要素中での、吸収物質濃度をｃ_j、平均光路長を
Ｌ_Ajとし、添字ｊは、ｊ番目の体積要素内における吸収
物質濃度、平均光路長の値であることを示し、加算Σは
ｊ＝１〜ｎについて行なう。ｎは、体積要素の総数であ
り、図４に示す例では、５×５＝２５である。 Ｉ_d／Ｉ₀＝Ｄ_S・ｅｘｐ（−ε（ｃ₁Ｌ_A1＋…＋ｃ_nＬ_An）） ＝Ｄ_S・ｅｘｐ（−εΣ（ｃ_jＬ_Aj）） …（数８） なお、光源からの光子は、被検体を構成する全ての体積
要素を一般に通過していくと想定するが、光子が通過し
ない体積要素では、Ｌ_Aj＝０であり、体積要素毎の平均
光路長Ｌ_Aj（ｊ＝１〜ｎ）を求める方法に関しては後述
するので、Ｌ_Aj（ｊ＝１〜ｎ）は既知として扱う。（数
８）において、減衰定数Ｄ_Sを直接求めることは困難で
あるため、照射光として波長の異なる２種類の光を用
い、各波長毎に計測を行ないＤ_Sを消去する。波長λ₁の
照射光、波長λ₂の照射光を用いた計測に関して、それ
ぞれ（数９）、（数１０）が成り立つ。ここで、平均光
路長Ｌ_Ajは波長に依存しないものとする。Ｄ_S、εは波
長に依存する。 Ｉ_d（λ₁）／Ｉ₀（λ₁）＝Ｄ_S（λ₁）・ｅｘｐ（−ε（λ₁）Σ（ｃ_jＬ_Aj）） …（数９） Ｉ_d（λ₂）／Ｉ₀（λ₂）＝Ｄ_S（λ₂）・ｅｘｐ（−ε（λ₂）Σ（ｃ_jＬ_Aj）） …（数１０） （数９）、（数１０）において、加算Σはｊ＝１〜ｎに
ついて行なう。

【００２１】（数９）及び（数１０）を（数１１）及び
（数１２）を用いて変数変換をすると、（数１３）及び
（数１４）が得られる。 μ（ｊ，λ₁）＝ε（λ₁）ｃ_j …（数１１） μ（ｊ，λ₂）＝ε（λ₂）ｃ_j …（数１２） Ｉ_d（λ₁）／Ｉ₀（λ₁）＝Ｄ_S（λ₁）・ｅｘｐ（−Σμ（ｊ，λ₁）Ｌ_Aj）） …（数１３） Ｉ_d（λ₂）／Ｉ₀（λ₂）＝Ｄ_S（λ₂）・ｅｘｐ（−Σμ（ｊ，λ₂）Ｌ_Aj）） …（数１４） 例えば、被検体が生体組織の場合には、波長λ₁、λ₂の
波長差の絶対値（│λ ₁−λ₂│）が、少なくとも１００
ｎｍ、特に３０ｎｍ以内であれば、近似的に次の（数１
５）が成り立つ。 Ｄ_S（λ₁）≒Ｄ_S（λ₂） …（数１５） （Ｉ_d（λ₁）／Ｉ_d（λ₂））（Ｉ₀（λ₂））／Ｉ₀（λ₁）） ＝ｅｘｐ（−Σ（μ（ｊ，λ₁）−μ（ｊ，λ₂））Ｌ_Aj））…（数１６） χ（ｊ，λ₁，λ₂）＝μ（ｊ，λ₁）−μ（ｊ，λ₂） ＝（ε（λ₁）−ε（λ₂））ｃ_j …（数１７） と置いて、（数１６）の両辺の自然対数（Ｌｏｇ_e）と
り（数１８）を得て、 −Ｌｏｇ_e（（Ｉ_d（λ₁）／Ｉ_d（λ₂））（Ｉ₀（λ₂））／Ｉ₀（λ₁））） ＝Σ（χ（ｊ，λ₁，λ₂）Ｌ_Aj） …（数１８） さらに、 ｙ＝−Ｌｏｇ_e（（Ｉ_d（λ₁）／Ｉ_d（λ₂））（Ｉ₀（λ₂））／Ｉ₀（λ₁））） …（数１９） と置いて、 ｙ＝Σ（χ（ｊ，λ₁，λ₂）Ｌ_Aj） …（数２０） を得る。（数２０）において、加算Σはｊ＝１〜ｎにつ
いて行なう。

【００２２】（数２０）は、線形方程式であり、少なく
ともｎ組の同様の方程式を用意すればｎ個ある未知数χ
（ｊ，λ₁，λ₂）（ｊ＝１〜ｎ）を求めることができ
る。１つの方程式は被検体に関する１つの計測位置（光
源１と検出器４の位置関係）毎に成立するので、図５に
示す様に、１つの光源１に対して検出器４を被検体２に
複数配置して複数位置で検出するか、あるいは複数の照
射位置から光を照射し、ｎ個の検出位置で光検出を行う
と、ｎ組の（数２０）の方程式が得られ、未知数χ
（ｊ，λ₁，λ₂）（ｊ＝１〜ｎ）と同数の連立方程式
（数２１）を得ることができる。但し、ｉはｉ番目の光
検出位置を示す添え字である。なお図５は、被検体２が
２次元体である例を示す。 ｙ₁＝Σ（χ（ｊ，λ₁，λ₂）Ｌ_A1j） … ｙ_i＝Σ（χ（ｊ，λ₁，λ₂）Ｌ_Aij） … ｙ_n＝Σ（χ（ｊ，λ₁，λ₂）Ｌ_Anj） …（数２１） （数２２）、（数２３）、（数２４）により、行列Ｙ、
Ａ、Ｘを定義すると、（数２１）で表される連立方程式
は（数２５）となる。以下の説明においては、Ｙを計測
値行列、Ａを光路長の空間分布行列、Ｘを未知数行列と
呼ぶ。 Ｙ＝│ｙ₁…ｙ_i…ｙ_n│_T …（数２２） （数２２）において、Ｔは転置を表し、Ｙはｎ行１列の
行列である。 Ａ＝│Ｌ_Aij│ …（数２３） （数２３）において、Ａはｉ行ｊ列の要素の値をＬ_Aij
とするｎ行ｎ列の行列である。 Ｘ＝│χ（１，λ₁，λ₂）…χ（ｊ，λ₁，λ₂）…χ（ｎ，λ₁，λ₂）│_T …（数２４） （数２４）において、Ｔは転置を表し、Ｘはｎ行１列の
行列である。 Ｙ＝ＡＸ …（数２５） （数２５）から未知数行列Ｘを求めるためには、（数２
６）で示すように、Ａの逆行列を両辺にかければよい。 Ａ^-1Ｙ＝Ａ^-1ＡＸ＝Ｘ …（数２６） ここで、求められた未知数行列Ｘの各成分は、（数２
７）で表される。 χ（１，λ₁，λ₂）＝μ（ｊ，λ₁）−μ（ｊ，λ₂） ＝（ε（λ₁）−ε（λ₂））ｃ_j …（数２７） （数２７）に照射光として用いた各波長λ₁、λ₂に対す
る吸光係数を代入して、各体積要素中の吸収物質濃度ｃ
_jを得る。

【００２３】以上の説明では、被検体２をｎ個の体積要
素に分割して、ｎ個の検出位置（光源１（光照射位置）
と検出器４の対）で光検出を行ない、ｎ組の（数２０）
の方程式が得たが、ｎをｎ’≧ｎなるｎ’に置き換え、
即ちｎ’個の検出位置（光源１（光照射位置）と検出器
４の対）で光検出を行ない、ｎ’組の（数２０）の方程
式が得ても、以下のようにして未知数行列Ｘを得ること
ができる。以上の説明において（数２２）、（数２３）
を以下のように置き換えればよい。即ち、（数２２）
は、 Ｙ＝│ｙ₁…ｙ_i…ｙ_n'│_T …（数２８） とする。（数２８）においてＴは転置を表し、Ｙはｎ’
行１列の行列である。

【００２４】（数２３）は、 Ａ＝│Ｌ_Aij│ …（数２９） とする。（数２９）において、Ａはｉ行ｊ列の要素の値
をＬ_Aijとするｎ’行ｎ列の行列である。このような置
き換えを行なった（数２５）の解は一般に、 Ｘ＝（（Ａ_T）ＷＡ）^-1（（Ａ_T）Ｗ）Ｙ …（数３０） として得られる。なお（数３０）において、Ｔは反転を
表し、Ａ_TはＡの反転行列であり、Ｗは重み行列を表す
ｎ’行ｎ’列の行列である。重み行列Ｗをｎ’行ｎ’列
の単位行列とするとき（数３０）は、 Ｘ＝（（Ａ_T）Ａ）^-1（Ａ_T）Ｙ …（数３１） となる。

【００２５】このように、被検体をｎ個の体積要素に分
割して、ｎ’≧ｎなるｎ’個の各種方向での検出位置
（光源１（光照射位置）と検出器４の対）で光検出を行
なうことにより、精度よく未知行列Ｘ、即ち被検体の各
体積要素中の吸収物質濃度ｃ_jを得ることができる。

【００２６】さらに以上の説明において、計測対象であ
る被検体が多成分の吸収物質を含む場合には、照射光に
用いる波長の種類を増やして、多成分の吸収物質の濃度
を求めることができる。例えば、被検体がαとβの２種
類の吸収物質を含む場合、（数５）を拡張して（数３
２）が成り立つ。但し、吸収物質α、βの分子吸光係数
と濃度をそれぞれ、ε_α、ε_β、ｃ_α、ｃ_βで表す。 Ｉ_d／Ｉ₀＝Ｄ_S・ｅｘｐ（−（ε_αｃ_α＋ε_βｃ_β）Ｌ_A） …（数３２） 図４に示すように、被検体２をｎ個の体積要素に分割す
る場合には、吸収物質が１種類ある場合と同様にして、
（数３２）から（数３３）が得られる。 Ｉ_d／Ｉ₀＝Ｄ_S・ｅｘｐ｛−ε_α〔ｃ_α1Ｌ_A1＋…＋ｃ_αnＬ_An〕 ＋ε_β〔ｃ_β1Ｌ_A1＋…＋ｃ_βnＬ_An〕｝ ＝Ｄ_S・ｅｘｐ｛−Σ（ε_αｃ_αj＋ε_βｃ_βj）Ｌ_Aj｝ …（数３３） （数３３）において、加算Σはｊ＝１…ｎについて行な
う。

【００２７】減衰定数Ｄ_Sを２種類の波長λ₁、λ₂の照
射光を用いて消去すると（数３４）が得られ、（数３
５）の置換を行なうと、（数３４）は（数３６）で表さ
れる。 （Ｉ_d（λ₁）／Ｉ_d（λ₂））（Ｉ₀（λ₂））／Ｉ₀（λ₁）） ＝ｅｘｐ｛−Σ〔（ε_α（λ₁）−ε_α（λ₂））ｃ_αj ＋（ε_β（λ₁）−ε_β（λ₂））ｃ_βj）〕Ｌ_Aj｝ …（数３４） χ（ｊ，λ₁，λ₂）＝（ε_α（λ₁）−ε_α（λ₂））ｃ_αj ＋（ε_β（λ₁）−ε_β（λ₂））ｃ_βj） …（数３５） （Ｉ_d（λ₁）／Ｉ_d（λ₂））（Ｉ₀（λ₂））／Ｉ₀（λ₁）） ＝ｅｘｐ｛−Σχ（ｊ，λ₁，λ₂）Ｌ_Aj｝ …（数３６） 同様にして２種類の波長λ₁、λ₃の照射光を用いて、
（数３６）と同様の式を得る。

【００２８】（数３４）、（数３６）において、加算Σ
はｊ＝１…ｎについて行なう。（数３６）の両辺の自然
対数（Ｌｏｇ_e）をとると線形な式となるので、以下、
先に説明した１種類の吸収物質がある場合と同様に、
ｎ’≧ｎなるｎ’個の検出位置（光源１（光照射位置）
と検出器４の対）で光検出を行ない、（数３６）を解く
ことができ、各体積要素内における変数χ（ｊ，λ₁，
λ₂）が求まる。ここで、波長λ₁と波長λ₃から得られ
る各体積要素内における変数χ（ｊ，λ₁，λ₃）と前記
の変数χ（ｊ，λ₁，λ₂）により連立方程式（数３７）
が成り立つ。 χ（ｊ，λ₁，λ₂）＝（ε_α（λ₁）−ε_α（λ₂））ｃ_αj ＋（ε_β（λ₁）−ε_β（λ₂））ｃ_βj） χ（ｊ，λ₁，λ₃）＝（ε_α（λ₁）−ε_α（λ₃））ｃ_αj＋ ＋（ε_β（λ₁）−ε_β（λ₃））ｃ_βj） …（数３７） ここで、吸収物質α、βの吸光係数が既知のとき、連立
方程式（数３７）に吸光係数の値を代入して、各体積要
素内の２種類の吸収物質の濃度ｃ_α、ｃ_βを求めること
ができる。従って、被検体内部にγ種類の吸収物質が存
在する場合、（γ＋１）種類の波長の照射光を用いて、
各吸収物質の濃度分布を求めることができる。

【００２９】以上の説明した被検体を複数の体積要素に
分割した場合において、特に被検体を分割しない場合
（即ち、分割の数をｎ＝１とする）にも、以上の方法は
適用できることは言うまでもない。この場合、（数２
９）はｎ’行１列の行列となり、異なる各種の方向位置
に配置される光照射位置と光検出位置の対により計測さ
れるｎ’個の計測値からなる計測値行列と、上記の対に
おけるそれぞれの平均光路長とから吸収物質の濃度が求
められる。もちろん、被検体が複数種類の吸収物質を含
む場合にも、被検体を分割せず分割の数をｎ＝１とし
て、（数３７）を適用できる。このようにして求められ
る吸収物質の濃度は、被検体の２次元的空間（被検体の
特定の断面の周囲に２次元的に上記の対が配置されると
き）、あるいは３次元的空間（被検体の外周に３次元的
に上記の対が配置されるとき）での平均的な値が求めら
れることになる。

【００３０】以上の説明において、行列Ｘを求めるため
には、行列Ａが既知でなければならない。行列Ａを構成
する成分Ｌ_Aijは、ｉ番目の計測位置（光源と検出器の
位置関係）におけるｊ番目の体積要素内の平均光路長を
表している。即ち、Ｌ_Aijは、各計測位置毎に検出器に
到達する光子の光路が既知であれば求めることができ
る。図６は、光源１から被検体２を通過し、検出器４に
到達する光子の光路３を、被検体２が２次元体である例
をとって模式的に表した図である。

【００３１】本発明では、光源から検出器に至る複数の
光路を求める方法として、被検体内の光散乱を模擬する
モンテカルロ法シミュレーションを用いる。被検体内の
光散乱を模擬するモンテカルロ法シミュレーションの例
として、Ｂ．Ｃ．Ｗｉｌｓｏｎ等の方法（" A Monte Ca
rlo model for the absorption and flux distribution
s of light in tissue ", Med. Phys., 10(6), 824ー830
(1983)に記載）が知られている。本願発明で使用するモ
ンテカルロ法は、よく知られた方法であるが、この方法
の概要を以下に説明する。光源から発した１光子のたど
る経路は各散乱点を表す、極座標（ｒ，θ，φ）で示さ
れる。各散乱点では、１光子が次に達する散乱点を与え
る、ｒ、θ、φの各座標値が確率的に定められ、散乱点
間の距離は散乱距離Ｌ_Sから、光子の進行方向を定める
ための角度方向は、θ方向の散乱角度θ_Sとφ方向の散
乱角度φ_Sとから決定される。なお、光子がたどる物体
中には吸収物質は含まれないものとする。各散乱毎の散
乱距離Ｌ_Sは散乱係数μ_Sを用いて（数３８）で表わさ
れ、 Ｌ_S＝−（ｌｏｇ_e（Ｒ））／μ_S …（数３８） である。Ｒは一様乱数（０＜Ｒ＜１）である。θ方向の
散乱角度θ_Sは、例えばＭｉｅ散乱理論から近似して得
られる位相関数の近似式である、Ｈｅｎｅｙ−Ｇｒｅｅ
ｎｓｔｅｉｎの式を用いて、一様乱数Ｒ（０＜Ｒ＜１）
と位相関数( phase function )（即ち１回散乱における
散乱強度の角度分布を示す関数）のコサイン平均ｇとを
使用して（数３９）で表わされ、 θ_S＝ｃｏｓ^-1〔｛（１＋ｇ²）−（１＋ｇ²）²／（１−ｇ＋２ｇＲ）²｝ ／（２ｇ）〕 …（数３９） である。なお、ｇは、位相関数をＰ_fとし、積分∫を立
体角ω＝０〜４πの範囲について行なうものとして、ｇ
＝（∫Ｐ_f・ｃｏｓθｄω）／（∫Ｐ_fｄω）により与え
れる。一般の生体関連物質では、ｇ≒０．９であること
が知られている。特に、ｇ＝０の場合には等方的散乱を
表し、（数３９）は、 θ_S＝ｃｏｓ^-1（２Ｒ−１） …（数４０） となる。φ方向の散乱角度φ_Sは一様な確率分布とみな
すことができ、一様乱数Ｒを用いた（数４１）で表わさ
れる。 φ_S＝２πＲ …（数４１） 光子のたどる各散乱点は、このように、一様乱数（０＜
Ｒ＜１）Ｒにより、確率的に与えられる（Ｌ_S，θ_S，φ
_S）により、順次与えられる。以上が、本願発明で使用
するモンテカルロ法の概要である。

【００３２】このモンテカルロ法による光子のたどる経
路シミュレーションを用いて、ｉ番目の計測位置（光照
射位置（光源）と光検出位置の組み合わせ）における１
光路を求めた時に、ｊ番目の体積要素内での光路長をｐ
（ｉ，ｊ）とする。実際の被検体内での光の伝播を実質
的に再現するに十分な光路数ｍ本（被検体形状に依存す
るが、少なくとも１０００本程度とする）を、上記で説
明したモンテカルロ法シミュレーションを用いて計算し
記憶しておき、各光路毎にｐ（ｉ，ｊ）を求める。ｋ
（ｋ＝１〜ｍ）番目の光路におけるｐ（ｉ，ｊ）をｐ
（ｉ．ｊ．ｋ）とすると、Ｌ_Aijは（数４２）で計算さ
れる。 Ｌ_Aij＝（Σｐ（ｉ，ｊ，ｋ））／ｍ …（数４２） （数４２）において、加算Σはｋ＝１〜ｍについて行な
う。全ての計測位置（ｉ）と、全ての体積要素（ｊ）に
関してＬ_Aijを計算すれば、平均光路長の空間分布行列
Ａが求められる。

【００３３】以上説明した本発明の手順は、照射光とし
て連続光を用い被検体を通過する光強度を計測信号とす
る場合、あるいは、照射光としてパルス光を用い、被検
体を通過するパルス光に対応する全光強度を計測信号と
する場合であるが、被検体を通過する光強度を時間分解
計測し、計測された光強度の１部を計測信号とする場合
にも、同様の手順で吸収物質の濃度分布を求めることが
できる。

【００３４】図７に示す様な時間スペクトルを有するδ
関数状のパルス光を、被検体、例えば生体の所定部分に
照射し、照射位置と対向する側の位置で通過光強度を、
ストリークカメラを用いて時間分解計測すると、図８に
示す様に時間的に拡がった通過光強度の時間スペクトル
Ｉ_d（t）が得られる。ただし、横軸はパルス光が照射さ
れてから通過光が検出されるまでの時間（検出時間）を
表す。検出光強度の時間スペクトルＩ_d（t）が、時間的
に広がる理由は、被検体が散乱体である場合には、各光
子の光路長に応じて検出される時間が変化するためであ
る。図６の様な系で、ｉ番目の各計測位置毎に検出光強
度の時間スペクトルＩ_di（t）を計測し、検出時間ｔ₁か
らｔ₂までＩ_di（t）を積分した値を検出光強度Ｉ_diとし
て用い計測値行列Ｙを求め、検出時間ｔ₁からｔ₂までの
平均光路長の空間分布行列Ａを求めれば、同様の手順で
未知数行列Ｘを求めることができる。検出時間ｔ₁から
ｔ₂までの平均光路長の空間分布行列Ａは、検出時間ｔ₁
からｔ₂までの間に検出器に到達する光子の光路のみに
ついて（数４２）を用いて計算すれば求めることができ
る。このように検出時間ｔ₁からｔ₂まで、時間ゲートを
設定してＩ_di（t）を積分した値を求める際に、時間ゲ
ート幅、（ｔ₂−ｔ₁）の大きさを、小さくすればする
程、最終的に再構成される被検体中の吸収体の３次元空
間分布の空間分解能、及び定量性を向上させることが可
能となる。

【００３５】以下、以上説明した内容を実際に適用した
実施例について詳細に説明する。まず簡単のため、本発
明を、生体の散乱特性を模擬した立方体の被検体に適用
した例を示す。なお、本願発明は被検体が立方体の形状
に限定されるものではなく、任意の形状に適用できるこ
とは言うまでもない。任意形状の被検体では、被検体の
形状を小さい立方体の集合として近似して表してもよい
し、被検体の内部を小さい立方体の集合とし、被検体の
外面を含む部分を、一部を欠損する小さい立方体で近似
して表してもよい。さらに、検査対象物体の特定部分の
みの空間、例えば、特定組織部分を含む空間部分、生体
表面もしくは表面近傍にある腫瘍を含む空間部分等を、
被検体部位として上記と同様の方法で表すことも可能で
ある。

【００３６】図９は、本実施例に用いた被検体の形状を
示す図である。被検体２は１辺２５ｍｍの立方体で、被
検体内部の位置を表す座標系（Ｘ０，Ｙ０，Ｚ０）は図
９に示す通りである。被検体２の各面を図９に示す様
に、Ｐ１〜Ｐ６とする。ここで、Ｐ１はＺ０＝０（ｍ
ｍ）の面、Ｐ２はＸ０＝２５（ｍｍ）の面、Ｐ３はＺ０
＝２５（ｍｍ）の面、Ｐ４はＸ０＝０（ｍｍ）の面、Ｐ
５はＹ０＝２５（ｍｍ）の面、Ｐ６はＹ０＝０（ｍｍ）
の面である。被検体２の散乱特性として、（数３８）で
使用する散乱係数を生体に近い値であるμ_S＝１．０
〔ｍｍ^-1〕とし、（数３９）で使用する先に説明した位
相関数のコサイン平均は、ｇ＝０とした。即ち、１回散
乱における散乱光強度の角度分布は一様であり等方的で
あるとした。

【００３７】また、図１０に示す様に、被検体２内部を
５×５×５の仮想の体積要素に分割し、それぞれの体積
要素の位置を（Ｖ_x，Ｖ_y，Ｖ_z）で表し、各体積要素番
号ｊは（数４３）で定義する。例えば、図１０の体積要
素５の位置は（５，５，１）であり、体積要素番号ｊ＝
２５となる。 ｊ＝Ｖ_x＋５（Ｖ_y−１）＋２５（Ｖ_z−１） …（数４３） 図１１は、照射光７を図９のＰ１の面上にある光照射位
置６−１から照射していることを示している。本実施例
では、全ての光照射位置、及び光検出位置は、被検体表
面に接する任意の体積要素の面の重心位置に設定されて
いる。図９の様に、光を照射するｉ番目の面Ｐｉに固有
の座標（Ｘ，Ｙ）_Piを定義すると、光照射位置６−１は
（２，２）_P1、光照射位置６−２は（２，２）_P2、光照
射位置６−３は（３，３）_P5で表される。また、本実施
例では、光検出位置を光照射位置が対向する側の被検体
外面にある全ての体積要素面の重心位置に設けた。例え
ば、図１１において、Ｐ１面上の光照射位置６−１から
照射光７を照射した時には、対向側のＰ３面上に２５点
の光検出位置を設けており、それぞれの光検出位置で検
出光強度が計測される（面Ｐｉの定義については図９を
参照）。本実施例では、被検体２を合計１２５個の体積
要素で分割しているので、それぞれの体積要素内で一様
な吸収物質濃度を未知数とすると、求めるべき１２５個
の未知数があることになる。未知数を決定するために
は、光計測位置（光源と光検出器の位置関係）も１２５
組だけ必要となる。そこで、本実施例では、照射点を
（２，２） _P1、（２，２）_P2、（４，４）_P3、（４，
２）_P4、（３，３）_P5の５点とし、光検出位置は各光照
射位置の対向側の面上に２５点設定し、合計１２５（５
×２５）組の光計測位置を用意した（面Ｐｉの定義につ
いては図９を参照）。

【００３８】本実施例では、被検体として、吸収物質の
位置の異なる２通りを用意した。これら被検体内部の吸
収物質の位置を図１２と図１３に示す。図１２と図１３
における吸収物質８の位置を（Ｖ_x，Ｖ_y，Ｖ_z）で表す
と、図１２の被検体では、吸収物質が（３，３，３）の
体積要素内に均一濃度で存在し、図１３の被検体では、
吸収物質が（３，３，３）と（１，３，１）の体積要素
内に均一濃度でそれぞれ存在している。照射光として２
波長λ₁、λ₂のパルス光を用い、各検出位置において通
過光強度を時間分解計測した。波長λ₁、λ₂のそれぞれ
に対する吸収物質の吸光係数は、ε（λ₁）＝１．０
〔(ｍｍ・ｍＭ)^-1〕、ε（λ₂）＝０．０〔(ｍｍ・ｍ
Ｍ)^-1)〕とし、吸収物質濃度は全て０．０１〔ｍＭ〕と
した。図１４に、先に説明したモンテカルロシュミレー
ションにより得られた、透過光強度（検出光強度）の時
間スペクトルＩ_di（λ₁，ｔ）、Ｉ_di（λ₂，ｔ）の１例
を示す。横軸は検出時間を表しており、光子が散乱され
ずに直接検出位置に到達する時刻を原点としている。な
お、図７、図８と同様に、図１４の横軸は、それぞれの
横軸の値に、光子の進行する媒体での光速度を乗算し
て、光子の進行距離として示すこともできる。本実施例
では、ｉ番目の計測位置における検出光強度Ｉ
_di（λ₁）、Ｉ_di（λ₂）として、各計測位置（ｉ）毎に
得られる検出光強度の時間スペクトルを検出時間ｔ₁＝
０〔ｐｓ〕から検出時間ｔ₂＝１５００〔ｐｓ〕まで積
分した値を用いた。なお、ｔ₂を１５００〔ｐｓ〕よ
り、小さい値に設定することにより、最終的に得られる
吸収体の３次元分布の空間分解能、及び定量精度は、更
に向上する。

【００３９】計測された検出光強度から、被検体内部に
存在する吸収物質の濃度を画像化する手順の概要を図１
５に示す。この画像化の手順は大きく分けると、計測値
行列Ｙを生成するｓｔｅｐ１、平均光路長の空間分布行
列Ａを生成するｓｔｅｐ２、行列Ｙと行列Ａ逆行列から
未知数行列Ｘを導出するｓｔｅｐ３、求められた未知数
行列Ｘから被検体内部の吸収物質濃度を求めこの空間分
布を表示するｓｔｅｐ４からなる。なお、モンテカルロ
法による光路の計算が、ｓｔｅｐ２の処理開始以前に既
になされ、光路データが記憶保存されている場合には、
記憶保存されている光路データを使用することは言うま
でもない。以下、各ｓｔｅｐにおける詳細な手順を、図
１６と図１７にフローチャートにより説明する。

【００４０】〔ｓｔｅｐ１〕： 計測値行列Ｙの生成。

【００４１】ｓｔｅｐ１−１： 各計測位置（ｉ）毎に
通過光強度を計測して、それぞれ検出光強度Ｉ
_di（λ₁）、Ｉ_di（λ₂）を求める。

【００４２】ｓｔｅｐ１−２： 各計測位置（ｉ）毎に
照射光に用いた２波長の照射光強度Ｉ_d0（λ₁）、Ｉ_d0
（λ₂）を計測する。

【００４３】ｓｔｅｐ１−３： 全計測位置（ｉ）に関
して計測終了後、ｓｔｅｐ１−１、及びｓｔｅｐ１−２
で求められた検出光強度Ｉ_di（λ₁）、Ｉ_di（λ₂）、及
びＩ _d0（λ₁）、Ｉ_d0（λ₂）を（数１９）に代入してｙ
_iを得る。

【００４４】ｓｔｅｐ１−４： ｓｔｅｐ１−３で得た
ｙ_i（i＝１〜１２５）から計測値行列Ｙ（数２２）を生
成する。

【００４５】〔ｓｔｅｐ２〕： 平均光路長の空間分布
行列Ａの生成。

【００４６】ｓｔｅｐ２−１： 計測対象である被検体
の３次元的な外形形状を計測して、外形形状を表す座標
データを演算処理装置に入力し、仮想被検体データを生
成する。外形形状を計測する方法としては、光切断法、
ＭＲＩ、Ｘ線ＣＴ等その他方法を用いる。本実施例では
被検体は１辺２５ｍｍの立方体であり、この外形形状デ
ータを入力する。

【００４７】ｓｔｅｐ２−２： 計測時に設定される計
測位置（光照射位置と光検出位置）の座標を入力する。
本実施例では先に説明した１２５組の計測位置を演算処
理装置に入力する。

【００４８】ｓｔｅｐ２−３： 計測対象となる被検体
の散乱特性として、被検体の散乱係数μ_S、及び位相関
数のコサイン平均ｇを演算処理装置に入力する。本実施
例では、散乱係数μ_S＝１．０〔ｍｍ^-1〕、指向性を表
わす変数としてｇ＝０を用いる。

【００４９】ｓｔｅｐ２−４： ｓｔｅｐ２−１〜ｓｔ
ｅｐ２−３で入力したパラメータを基に、モンテカルロ
法を用いて吸収物質の存在しない被検体内の光散乱を模
擬する計算を行なう。設定された各計測位置（ｉ）毎に
複数の光路を計算し、各計測位置（ｉ）毎の各体積要素
（ｊ）毎における平均光路長Ｌ_Aij（平均光路長の空間
分布）をメモリに記憶しておく。この時、通過光強度を
時間分解計測する場合には、計測信号として用いる検出
時間ｔ₁からｔ₂の間に到達する光子の光路のみから、平
均光路長の空間分布Ｌ_Aijを求めて、メモリに記憶して
おく。本実施例では、各計測位置毎に３００００本の光
路を計算し、検出時間は０〔ｐｓｅｃ〕から１５００
〔ｐｓｅｃ〕とした。本ステップで用いるモンテカルロ
法に関しては先に概要を説明した。なお、モンテカルロ
法による光路の計算が、ｓｔｅｐ２−４の処理開始以前
に既になされ、光路データが記憶保存されている場合に
は、記憶保存されている光路データを使用することは言
うまでもない。光路データが記憶保存されている場合
の、各計測位置（ｉ）毎の各体積要素（ｊ）毎における
平均光路長Ｌ_Aijを求めるための処理方法に関しては、
後で詳細に説明する。

【００５０】ｓｔｅｐ２−５： ｓｔｅｐ２−４で全計
測位置（ｉ）と全体積要素内（ｊ）に関して平均光路長
Ｌ_Aij（ｉ＝１〜１２５、ｊ＝１〜１２５）を求めた
後、平均光路長の空間分布行列Ａ（数２３）を生成す
る。ただし、計測対象である被検体の形状が予想される
場合、あるいは被検体の形状と類似した形状を用いる場
合には、あらかじめ平均光路長の空間分布行列Ａを計算
しておく。例えば、生体頭部や乳房が被検体である場合
には、複数種類の頭部形状や乳房形状について、平均光
路長の空間分布行列Ａを計算しておき、実際に計測対象
となる生体頭部や乳房と最も近い形状に関して計算した
平均光路長の空間分布行列Ａを選択する。また、生体頭
部等は、球体あるいは楕円体に近似して計算する。

【００５１】〔ｓｔｅｐ３〕： 未知数行列Ｘの導出。

【００５２】ｓｔｅｐ３−１： ｓｔｅｐ１−４で生成
された計測値行列Ｙとｓｔｅｐ２−５で生成された平均
光路長の空間分布行列Ａの逆行列Ａ^-1を求め（数２６）
に代入して、未知数行列Ｘを求める。もちろん（数３
０）により未知数行列Ｘを求めてもよい。

【００５３】ｓｔｅｐ３−２： 未知数行列Ｘの各成分
（数２７）より、照射光の各波長に対する吸収物質の吸
光係数を用いて、各体積要素毎の吸収物質濃度を求め
る。

【００５４】〔ｓｔｅｐ４〕： 画像表示。

【００５５】ｓｔｅｐ３で得られた吸収物質濃度の空間
分布から、任意断層面の吸収物質濃度分布を画像表示す
る。

【００５６】図１８から図２３に、上記で説明したフロ
ーチャートに従って求められた吸収物質濃度の空間分布
の結果を表示する。図１８から図２０は、図１２に示し
た被検体における吸収物質濃度の空間分布を求めた結果
である。図１８はＶ_y＝２の面、図１９はＶ_y＝３の面、
図２０はＶ_y＝４の面を表示している。図２１から図２
３は、図１３に示した被検体における吸収物質濃度の空
間分布を求めた結果である。図２１はＶ_y＝２の面、図
２２はＶ_y＝３の面、図２３はＶ_y＝４の面を表示してい
る。図１８から図２３において、ｘ＝、ｙ＝、ｚ＝のそ
れぞれに続いて記される数は対応する座標の値を示す。
以上の計算では、被検体の散乱係数として、生体に近い
散乱係数、μ＝１．０〔ｍｍ^-1〕を設定し、（数３９）
中の位相関数のコサイン平均をｇ＝０（等方散乱）と仮
定している。一般の生体関連物質では、ｇ≒０．９であ
ることが知られているが、上記の計算ではｇ＝０の仮定
にもかかわらず、得られた結果は、明確に予め設定され
た被検体における吸収物質濃度の空間分布を再現してお
り、吸収物質濃度が異なる５×５×５のサイズの空間が
明確に識別されており、が空間分解能も良好である。

【００５７】図２４に、Ｘ線ＣＴのアルゴリズム（ Fil
tered Back Projection 法）を用いた結果例を示す。図
２４は、外径２５ｍｍ、高さ１００ｍｍを有する円柱
（散乱係数１．０ｍｍ^-1である）の中心軸位置に、外径
５．ｍｍ、高さ１００ｍｍを有する円柱（空間分解能の
評価において広く仮定される完全吸収体とし、吸収系数
は無限大とする）の吸収体を配置した被検体に関する再
構成画像である。画像の再構成に使用した計測データ
は、被検体に照射する光をδ関数と見做して、先に概要
を説明したモンテカルロ法を使用したシミュレーション
により、図１４に示すような検出光強度の時間スペクト
ルを生成した。画像の空間分解能を上げるために、時間
ゲート幅を小さく設定して、検出光強度の時間スペクト
ルの０〔ｐｓｅｃ〕から３００〔ｐｓｅｃ〕までの積分
信号を計測信号とした。断層像は、計算により円柱の高
さ方向の中心位置で Filtered Back Projection 法で求
めた。

【００５８】図２４では、再構成画像の結果を立体図で
示しており、高さ方向の値の最大値を１．０、高さ方向
の表示の最小間隔を０．０５として示している。図１８
〜図２３と図２４との比較から、本発明による画像化方
法は従来の方法と比較して明らかに空間分解能が向上し
ていることがわかる。 Filtered Back Projection 法で
は、散乱の影響を十分考慮することができず、検出光強
度の時間スペクトルの積分時間間隔を仮に短くした計測
信号を使用しても、Ｓ／Ｎ比が劣化してしまい高い分解
能を得ることは困難である。図２４の結果では、外径５
ｍｍの吸収体は、空間分解能は半値幅で評価すると１５
ｍｍであり、５ｍｍ×５ｍｍ×５ｍｍの大きさの吸収体
を再構成した、図１８〜図２３に示す結果では、吸収体
は明確に鋭く（シャープ）に検出されている。このよう
に、本発明での空間分解能は、 Filtered Back Project
ion 法に比較して約３倍以上改善されている。

【００５９】図２５は、本発明の方法が適用される装置
構成を説明する図である。光源９は少なくとも２波長の
光を発する光源で、例えばチタンサファイアレーザ、も
しくは複数の半導体レーザで構成される。波長の切り替
えは、計算機２４で制御する。まず、光源９から１波長
（λ₁）の光が選択され、光ファイバー、レンズ等の光
学系で構成される光源用導波路１０(以下、特に断りの
無い場合、導波路は光ファイバーあるいは光学系で構成
される)を通り、光は分波器１１で参照光用導波路１２
と光源用光ファイバー１３の２系統に分岐される。スト
リークカメラ等の光検出器２０により、参照光用導波路
１２に分岐した光は照射光強度Ｉ₀をもつ参照光として
計測される。光源用光ファイバー１３に分岐した光は光
スイッチ１４に入射する。光スイッチ１４は、被検体２
の表面上に配置された複数の光照射用光ファイバー１５
の１本を選択し、光源用光ファイバー１３と接続され、
光源用光ファイバー１３からの光は被検体２の任意位置
から被検体２内部に照射される。

【００６０】ここで、光照射用光ファイバー１５は、光
ファイバー固定ホルダー１６で位置を固定維持されてい
る。被検体２を通過した光は、被検体２の表面上に配置
された複数の光検出用光ファイバー１７に入射する。光
スイッチ１８は、複数の光検出用光ファイバー１７の中
から必要な光ファイバーのみを選択し、検出器用光ファ
イバー１９と接続する。選択された光検出用光ファイバ
ー１７を通過する検出光は光検出器２０に導かれ、各検
出器用光ファイバー１９毎に検出光強度の時間分解信号
が計測される。光源の強度信号と、検出光強度の時間分
解信号は信号伝送ケーブル２１でＡ／Ｄ変換器２２に転
送されデジタル信号に変換される。デジタル信号に変換
された光源の強度信号と検出光強度の時間分解信号は、
信号伝送ケーブル２３で計算機２４に転送され、計算機
内の記憶装置に保存される。

【００６１】引続き、上記で用いた波長（λ₁）と異な
る波長（λ₂）の光を選択し、上記と同様の計測を行な
う。１照射位置に関して計測終了後、前回の照射位置あ
るいは検出位置と異なる位置における光計測をするた
め、光スイッチ１４及び１８を光スイッチ制御ケーブル
２６及び２７を介して計算機２４により制御して、光源
用光ファイバー１３と光照射用光ファイバー１５の接続
位置及び検出器用光ファイバー１９と光検出用光ファイ
バー１７の接続位置を切り替える。必要な計測が終了し
た後、先に説明した図１５、図１６、図１７のフローチ
ャートの手順に従い、吸収物質濃度の空間分布を求め表
示装置２５に表示する。

【００６２】図２６に、光照射用光ファイバー１５と光
検出用光ファイバー１７の配置例を示す。被検体の面の
番号Ｐｉを図９と同様に定義し、照射位置と検出位置を
図９と同様にする場合には、照射用光ファイバー１３を
（２，２）_P1、（２，２）_P2、（４，４）_P3、（４，
２）_P4、（３，３）_P5の５点に配置し、検出用光ファイ
バー１７はＰ１〜Ｐ５の各面の各体積要素面の重心位置
に配置する。ただし、前記の５点の光照射位置には、光
照射用光ファイバー１５と光検出用光ファイバー１７の
２本の光ファイバーが配置されるので、それぞれの光フ
ァイバーの先端位置を多少ずらして配置するものとす
る。光照射用光ファイバー１５、及び光検出用光ファイ
バー１７には、例えば径１ｍｍの石英ファイバーを用い
る。

【００６３】図２７に、被検体２を生体頭部とする場合
の、光照射用光ファイバー１５と光検出用光ファイバー
１７の配置例を示す。被検体が生体頭部の場合には、ヘ
ルメット状の光ファイバー固定具２８を用いて光照射用
光ファイバー１５、あるいは光検出用光ファイバー１７
を固定維持する。光ファイバー固定具２８には、光照射
用光ファイバー１５、あるいは光検出用光ファイバー１
７を挿入する穴が設けられ、各光ファイバーを任意の穴
位置に挿入する。被検体が生体の様な場合にも、仮想の
体積要素５を立方体とするが、図２７に示す様に、表面
に接する体積要素５は立方体とはならない。

【００６４】ここで、被検体を１０００個の体積要素に
分割する場合には、計測位置（照射位置と検出位置の関
係）を少なくとも１０００組用意する必要があり、計測
位置を１０００組設定するために、図２５に示す光照射
用ファイバー１５と光検出用光ファイバー１７の配置方
法としては、様々の方法がある。例えば、光照射用光フ
ァイバー１５を少なくとも２０本、光検出用光ファイバ
ー１７を少なくとも１０００本、それぞれ配置し、少な
くとも２０本の光照射用光ファイバー１５から光を照射
した時に、少なくとも１０００本の光検出用光ファイバ
ー１７の中から少なくとも５０本を選択して、検出光強
度を計測すればよい。この計測は、図２５と同様な装置
構成で実現できる。

【００６５】生体に照射光として近赤外光を用いると、
吸収物質としては酸化ヘモグロビン（以下ＨｂＯ₂と略
す）と脱酸素化ヘモグロビン（以下Ｈｂと略す）が支配
的であり、これら２成分の吸収物質濃度の空間分布を求
めることができる。この場合は、照射光として３種類の
波長（λ₁，λ₂，λ₃）の光を用いて、各波長毎に光計
測を行い、（λ₁，λ₂）の波長を照射光として用いた計
測から図１５、図１６、図１７に記載の、ｓｔｅｐ１か
らｓｔｅｐ３までと同様のフローチャートの手順に従
い、未知行列Ｘ、即ちχ（ｊ，λ₁，λ₂）を求め、次に
（λ₁，λ₃）の波長を照射光として用いた計測から同様
にしてχ（ｊ，λ₁，λ₃）を求め、得られたχ（ｊ，λ
₁，λ₂）、χ（ｊ，λ₁，λ₃）から連立方程式（数３
７）を立て、連立方程式（数３７）に各波長の光に対す
る吸収物質の吸光係数εを代入して解き、各体積要素毎
のＨｂの濃度ｃ（Ｈｂ，ｊ）と、ＨｂＯ₂の濃度ｃ（Ｈ
ｂＯ₂，ｊ）を求める。

【００６６】図２８に、連立方程式を解く際に用いるＨ
ｂとＨｂＯ₂の吸光係数の波長スペクトルを示す（ S.Wr
ay 等, " Characterization of the near infrared abs
orption spectra of cytochrome aa3 and hemoglobin f
or the non-invasive monitoring of cerebral oxygena
tion "; Biochim. Biophys. Acta 933, 184-192(1988)
より数値を抜粋してグラフ化した）。図２８中で、実線
がＨｂの吸光係数の波長スペクトルを、点線がＨｂＯ₂
の吸光係数の波長スペクトルを、それぞれ示す。

【００６７】図２８に示すような吸収物質に関する吸光
係数の波長スペクトルをあらかじめ計算機内の記憶装置
に記憶しておき、照射光として用いる波長に対する吸光
係数を求める。頭部以外の生体組織、例えば乳房等の様
な部位を被検体とする場合も、専用の光ファイバー固定
具を使用し、照射光として近赤外光を用いて、被検体内
部のＨｂとＨｂＯ₂の空間分布を求めることができる。

【００６８】また、被検体中にγ種類の吸収物質が存在
する場合には、照射光として（γ＋１）種類の波長の光
を用いて、被検体中のγ種類の吸収物質濃度の空間分布
を得ることができる。例えば、近赤外光を用いて筋肉中
のＨｂとＨｂＯ₂とミオグロビンの空間分布を求める場
合には、４種類の波長の照射光を用いればよい。

【００６９】被検体の一部の領域についての吸収物質濃
度の空間分布を求める方法を、図２９を用いて説明す
る。図２９では、生体頭部を被検体２とし、被検体の一
部である注目領域３１内の吸収物質濃度の空間分布を求
めることを示している。図２９において、注目領域３１
の形状は２５ｍｍ角の立方体で、５×５×５の合計１２
５個の体積要素に分割し、各体積要素内の吸収物質濃度
を求める。注目領域３１の近傍に配置した光照射用光フ
ァイバー１５から照射光２９を被検体２に照射し、被検
体２内で散乱された光が、注目領域３１の近傍に配置さ
れた光検出用光ファイバー１７に到達し検出光３０とし
て検出される（この時検出される光は、光照射位置と同
一方向に配置された検出器により検出される、反射光で
ある）。但し、光照射用光ファイバー１５から光検出用
光ファイバー１７までの光路３は、注目領域３１外の領
域を通過する光路もあるので、光が到達する領域３２
（以下では、光到達領域と略記する）の内部全体につい
ての吸収物質濃度の空間分布を求めなければならない。
光到達領域３２の形状は、あらかじめモンテカルロ法に
よるシミュレーションや、拡散方程式を用いたシミュレ
ーション等で決定する。

【００７０】例えば、光照射用光ファイバー１５から光
検出用光ファイバー１７に到達する光の９５％が通過す
る領域を包絡するように、光到達領域３２を定める。

【００７１】光到達領域３２内の吸収物質濃度の空間分
布を求めるには、光到達領域３２を注目領域３１と同様
に体積要素に分割して求める方法と、注目領域３１外の
光到達領域３２に含まれる領域を１体積要素と扱って求
める方法がある。両方法とも体積要素数が異なるだけで
あり、本質的には同じであるので、本実施例では注目領
域３１外の領域を１体積要素として説明する。従って、
光吸収物質濃度の空間分布を求めるべき光到達領域３２
内には、体積要素数が１２６個存在することになり、異
なる計測位置（照射位置と検出位置の位置関係）を少な
くとも１２６組用意しなければならない。計測位置の配
置方法は無数にあるが、図２９において、例えば光照射
用光ファイバー１５を注目領域３１付近に少なくとも９
本配置し、光検出用光ファイバー１７を注目領域付近に
少なくとも１４本配置すれば、少なくとも１２６（９×
１４）組の計測位置を用意できる。図２５と同様な装置
構成で、各計測位置毎に被検体２を通過してくる検出光
３０の強度を計測し、図１５、図１６、図１７と同様の
フローチャートの手順で、光到達領域３２内の吸収物質
濃度の空間分布を求めることができる。ただし、照射光
として近赤外光を用い、各体積要素中のＨｂＯ₂とＨｂ
の濃度を求める場合には、照射光として３種類の波長を
用いる。また、平均光路長の空間分布行列Ａを計算する
際には、光到達領域３２の外部形状が入力される。光到
達領域３２の外部形状は、あらかじめ求めることができ
るので、平均光路長の空間分布行列Ａも、あらかじめ計
算してメモリに記憶しておいても構わない。

【００７２】以下、被検体の任意の光照射位置から任意
の光検出位置までの、複数の光路を計算して一旦メモリ
に光路データとして記憶しておき、記憶された光路デー
タに基づいて、平均光路長の空間分布行列Ａをもとめる
ための手順（以下では、簡単のため回転移動型モンテカ
ルロ法と略記する）について図を用いて、詳細に説明す
る。任意の光照射位置から任意の光検出位置までの複数
光路は、被検体内の光散乱を被検体の散乱特性を用いて
模擬するモンテカルロ法で計算して求めた。従来の光散
乱を模擬し光路を計算する従来のモンテカルロ法は、ウ
ィルソン等の行った方法が一般的であり、先に説明し
た。被検体の任意の光照射位置から任意の光検出位置ま
での光路を、従来のモンテカルロ法を用いて、必要に応
じてその都度計算してもよいが、その計算のために膨大
な計算時間を要するため実用的ではない。本発明では、
光路データを記憶しておき、これを使用することより、
その都度光路を計算する従来方法より、約１００倍以上
高速化できる。

【００７３】以下では被検体形状として、図１２あるい
は図１３と同様な立方体を用いて説明する。図３０は、
回転移動型モンテカルロ法を説明するために想定した被
検体形状と、照射光２９を照射する光照射位置３３と検
出光３０を検出する光検出位置３４の相対位置関係を示
す図である。図３０において光照射位置３３から光検出
位置３４に至るまでの複数の光路を回転移動型モンテカ
ルロ法で求める手順を、代表例として以下説明する。図
３０において、立方体の被検体２は、仮想的に５×５×
５の体積要素に分割され、各体積要素５は立方体であ
る。回転移動型モンテカルロ法の概念を図３１に示す。
図３１は、簡単のため、光照射位置３３と光検出位置３
４の位置を含む断面を示しており、この断面内に存在す
る光路３７が示される。最初に、半径が被検体２の任意
の２点をとりその間の距離の最大値（図３０の被検体形
状では、頂点３５−１と頂点３５−２の距離）と等し
い、球形の仮想散乱体３６を計算機内で発生する。球形
の仮想散乱体３６（図３０では図示せず）の中心Ｏから
発射される１光子は、被検体２と等価な散乱特性で多重
散乱を繰り返し、仮想散乱体３６の境界面まで到達す
る。ここで与えている散乱特性とは、先に説明した被検
体の散乱係数と位相関数で与えられる。

【００７４】仮想散乱体３６の境界面に到達するまでの
光子の散乱点を、光路データとして記憶保存しておくこ
とにより、１光子の仮想散乱体３６中における光路３７
を再現することができる。発射された１光子が境界面に
到達した後、同様の過程を繰り返して複数の光子に関す
る光路を得る。統計的な精度を実質的に補償するに十分
な光路数が得られた後、図３１に示す様に、被検体２の
外形形状を有する仮想被検体３８を発生する。被検体２
上の光照射位置３３と光検出位置３４の位置に対応し
て、仮想被検体３８上に光照射位置３３と光検出位置３
４を設定する。続いて、仮想被検体３８上の光照射位置
３３を、球形の仮想散乱体３６の中心位置に配置し、球
形の仮想散乱体３６の中心Ｏを回転中心として仮想被検
体３８を回転させて、仮想被検体３８と光路の相対位置
を変化させ、仮想被検体３８上の光検出位置３４と、上
記の過程で記憶保存されている光路とが交差する点を求
め、中心Ｏからこの交差点に至るまでの光路を再現する
散乱点の座標を仮想被検体３８に設定されている座標系
に変換する。即ち、球形の仮想散乱体３６の中心Ｏを回
転中心として、仮想被検体３８上の光照射位置３３と光
検出位置３４を結ぶ線分を回転させて、この線分の先端
点である光検出位置３４と記憶保存されている光路とが
交差する点を求める。上記の保存されている光路データ
のそれぞれに関して、同様の操作を繰り返し、仮想被検
体３８上の光照射位置３３から光検出位置３４に至るま
での複数の光路が求められる。例えば、図３１中に示す
様に、１光路に関して処理が終了した後、仮想散乱体３
８を中心Ｏの周りに矢印の方向に回転し、次の光路デー
タに関する処理を行う。求められた仮想被検体３６上の
光照射位置３３から光検出位置３４に至るまでの光路デ
ータ（各光路に属する散乱点のデータ）から、各体積要
素毎の光路長を計算することができ、各体積要素毎の光
路長を記憶しておき、被検体２上の光照射位置３３と光
検出位置３４の計測位置における、平均光路長の空間分
布を（数４２）から計算することができる。以上が、回
転移動型モンテカルロ法の２次元での説明であるが、３
次元においても同様の処理ができることは言うまでもな
い。被検体の形状が一般の場合にも、被検体を複数の空
間要素に分割して、被検体の外形形状を有する仮想被検
体を処理装置において生成して、以上の説明と同様の処
理を行なって、任意形状の被検体に関しての、光照射位
置と光検出位置とが設定されたとき、この光照射位置と
光検出位置とを通過する平均光路長の空間分布を（数４
２）から計算することができる。

【００７５】以上説明した処理の手順をフローチャート
として図３２、及び図３３に示し、以下にその詳細を各
ステップ毎に説明する。ただし、ここでは図１７に示さ
れている様に、既に被検体の外形形状を表すデータが入
力されており、さらに、被検体の散乱特性も入力されて
いるものとする。

【００７６】ｓｔｅｐ１： 被検体の表面上に任意の２
点を取り、その距離の最大値が半径となるような球の仮
想散乱体３６を発生する。仮想散乱体３６の中心から光
子を出射し、仮想散乱体境界面から光子が出ていくまで
の光路データを保存する。統計精度を実質的に補償する
ために、少なくとも１０００本の光路データを保存す
る。仮想散乱体３６の散乱特性は、被検体の散乱特性を
反映する。

【００７７】ｓｔｅｐ２： 体積要素に分割された仮想
被検体を発生し、被検体上の光照射位置に対応して、光
照射位置が仮想被検体に設定される。仮想被検体上の光
照射位置を球形の仮想散乱体中心に配置する。

【００７８】ｓｔｅｐ３： ｓｔｅｐ１で保存されてい
る任意の光路データから、光路の１本を選択し、被検体
上の光照射位置と同位置に設定される仮想被検体上の光
検出位置と、光路が交差するように被検体を回転する。
この時、回転中心は仮想散乱体の中心（即ち光照射位
置）とする。本ステップで任意光路を選ぶ際に、既に選
択された光路は重複して選択しない。

【００７９】ｓｔｅｐ４： ｓｔｅｐ３で光照射位置と
光検出位置を通過した光路が、光路中で仮想被検体の境
界面と交差するか否かを判断する。交差する場合はｓｔ
ｅｐ３に戻り、次の光路に関して処理を行ない、交差し
ない場合はｓｔｅｐ５に進む。

【００８０】ｓｔｅｐ５： 光照射位置から光検出位置
への到達時間ｔを総光路長から計算し、検出時間条件ｔ
₁、ｔ₂に対して、ｔ₁≦ｔ≦ｔ₂の条件を満たしているか
否かを判断する。条件を満たしていない場合は、ｓｔｅ
ｐ３に戻り次の光路データに関して処理を行ない、条件
を満たしている場合は、ｓｔｅｐ６に進む。 ｓｔｅｐ６： 光照射位置から光検出位置に至るまでの
光路の散乱点の座標を、仮想散乱体の座標系から仮想被
検体の座標系に変換した後、各体積要素毎の光路長ｐ
（ｉ，ｊ，ｋ）を計算し、記憶保存する。 ｓｔｅｐ７： 保存されている全光路データに関する以
上の各ｓｔｅｐの処理が終了したか否かを判断する。計
算終了でない場合はｓｔｅｐ３に戻り次の光路に関して
処理を行ない、処理終了の場合はｓｔｅｐ９へ進む。

【００８１】ｓｔｅｐ８： 本ステップまでに処理され
た光路の累積数ｍは、（数３８）の右辺の分母に相当す
るｍであり、検出器への到達時間〔ｔ₁〜ｔ₂〕の時間区
間に到達する光子数ｍと等しく、ｓｔｅｐ６で保存され
ている各体積要素毎の光路長ｐ（ｉ，ｊ，ｋ）から、各
体積要素毎の累積光路長（（数３８）の右辺の分子）を
計算し、光路の累積数ｍで除算する。本ｓｔｅｐでｉ番
目の光計測位置（光照射位置と光検出位置の組）におけ
る、ｊ番目の体積要素毎の平均光路長Ｌ_Aij（平均光路
長の空間分布）を求めたことになる。異なる計測位置に
関して平均光路長の空間分布を計算する場合には、ｓｔ
ｅｐ９に進む。

【００８２】ｓｔｅｐ９： 全計測位置（光照射位置と
光検出位置の組）に関して、処理を終了したか否かを判
断する。処理終了していない場合は、ｓｔｅｐ２に戻り
次の計測位置に関して処理を行なう。処理終了の場合に
は、ｓｔｅｐ１０に進む。

【００８３】ｓｔｅｐ１０： 全計測位置に関する平均
光路長の空間分布を記憶保存する。

【００８４】以上のフローチャートに示す手順では、全
計測位置に関して必ずしも処理する必要はない。幾何学
的に等価な関係にある異なる計測位置（光照射位置と光
検出位置の組）が複数ある場合に関しては、その中から
１つの計測位置を選択して、光路を求め、得られた光路
が、他の幾何学的に等価な計測位置を通過するように座
標変換して、処理時間を短縮できる。この処理は、前記
のｓｔｅｐ６で行えばよい。例えば、図３０において、
光照射位置３３が（３，３）_P1であり、光検出位置３４
が（４，３）_P2である計測位置は、光照射位置が（３，
３）_P1であり、光検出位置が面Ｐ２と平行な面Ｐ４上の
点である（４，３）_P4、等の計測位置と幾何学的に等価
と考える（面Ｐｉの定義は図９参照）。

【００８５】以上が、回転移動型モンテカルロ法の説明
であるが、同様の考え方により、図３４に示すように、
仮想散乱体３６（十分大きいサイズを有する空間であれ
ば任意の形状でよい）の１点（Ｏ₀）から、発する１光
子の複数の散乱点を、光路データとして記憶保存してお
き、光路３７を再現することができる（この方法を、移
動回転型モンテカルロ法と呼ぶ）。なお、図３４では、
図３１と同様に、簡単のため、図３０に示される光照射
位置３３と光検出位置３４の位置を含む断面を示してお
り、この断面内に存在する光路３７が示される。仮想散
乱体３６中における仮想被検体３８上の光照射位置３３
を、仮想散乱体３６の任意の位置Ｏに配置し、仮想散乱
体３６の中心Ｏを回転中心として仮想被検体３８を回転
させて、仮想被検体３８と光路３７の相対位置を変化さ
せ、仮想被検体３８上の光検出位置３４と、保存された
光路３７とが交差する点を求め、点Ｏからこの交差点に
至るまでの光路を再現する散乱点の座標を仮想被検体３
８に設定されている座標系に変換し、各体積要素毎の光
路長ｐ（ｉ，ｊ，ｋ）を計算し、記憶保存する。次い
で、光照射位置３３の位置を点Ｏから、光路の沿って点
Ｏ’に移動させて同様の処理を行なう。点Ｏと点Ｏ’と
の間の距離は、仮想被検体の光照射位置３３と光検出位
置３４との間の距離より大きくすることが望ましい。そ
の他の処理に関しては、先に説明した回転移動型モンテ
カルロ法の場合と同様である。以上が、移動回転型モン
テカルロ法の２次元での説明であるが、３次元において
も同様の処理ができることは言うまでもない。同様にし
て、被検体の形状が一般の場合にも、被検体に光照射位
置と光検出位置とが設定されたとき、この光照射位置と
光検出位置とを通過する平均光路長の空間分布を（数３
８）から計算することができる。

【００８６】

【発明の効果】従来の画像化方法では、空間分解能、定
量性、計算時間の点で問題点があったが、本発明の方法
では、従来方法に比較して、空間分解能で約３倍以上の
改善がなされた３次元的な吸収物質の空間分布を得るこ
とができ、１０〜１００倍程度計算時間が短縮され、局
所の最適解に落ち込むことなく解を得ることができ、さ
らに定量性の改善が可能である。また、回転移動型、も
しくは移動回転型モンテカルロ法を用いることにより、
従来のモンテカルロ法を適用する方法と比較して、平均
光路長の空間分布を約１００倍以上の速度で求めること
ができる。

【図面の簡単な説明】

【図１】吸収物質を一様に含む被検体中での光路を含む
断面図。

【図２】異なる層状の吸収物質からなる被検体中での光
路を含む断面図。

【図３】散乱体に吸収物質を含む被検体中での光路を含
む断面図。

【図４】散乱体に吸収物質を含む被検体中を複数体積要
素に分割したときの各体積要素での光路を模式的に示す
図。

【図５】被検体に光源と検出器を配置する例を示す図。

【図６】検出器に到達する光子の光路例を示す図。

【図７】照射光の時間スペクトルの例を示す図。

【図８】被検体を通過し検出される光の時間スペクトル
の例を示す図。

【図９】実施例に用いた被検体の形状を示す図。

【図１０】実施例に用いた被検体を体積要素に分割する
例を示す図。

【図１１】実施例に用いた光照射位置の配置例を示す
図。

【図１２】実施例に用いた被検体内での吸収体の位置を
示す図。

【図１３】実施例に用いた被検体内での吸収体の位置を
示す図。

【図１４】実施例に用いた計算により求めた検出光の時
間スペクトル例を示す図。

【図１５】被検体内の吸収物質濃度の空間分布を求める
処理を示すフローチャート。

【図１６】被検体内の吸収物質濃度の空間分布を求める
処理を示すフローチャート。

【図１７】被検体内の吸収物質濃度の空間分布を求める
処理を示すフローチャート。

【図１８】本発明により得られた吸収物質濃度の空間分
布の例を示す図。

【図１９】本発明により得られた吸収物質濃度の空間分
布の例を示す図。

【図２０】本発明により得らえた吸収物質濃度の空間分
布の例を示す図。

【図２１】本発明により得らえた吸収物質濃度の空間分
布の例を示す図。

【図２２】本発明により得られた吸収物質濃度の空間分
布の例を示す図。

【図２３】本発明により得られた吸収物質濃度の空間分
布の例を示す図。

【図２４】従来方法により得られた吸収物質濃度の空間
分布の例を示す図。

【図２５】本発明が適用される装置の構成例を示す図。

【図２６】被検体に光照射用、及び光検出用光ファイバ
ーを配置する例を示す図。

【図２７】生体頭部に光照射用光、及び光検出用ファイ
バーを配置する例を示す図。

【図２８】酸素化ヘモグロビンと脱酸素化ヘモグロビン
の吸光係数の波長スペクトル。

【図２９】本発明を生体の１部の領域に適用する例を示
す図。

【図３０】本発明による回転移動型モンテカルロ法を適
用する被検体の例を示す図。

【図３１】本発明による回転移動型モンテカルロ法の概
念を説明するための図。

【図３２】本発明による回転移動型モンテカルロ法によ
り平均光路長の空間分布を求める手順を示すフローチャ
ート。

【図３３】本発明による回転移動型モンテカルロ法によ
り平均光路長の空間分布を求める手順を示すフローチャ
ート。

【図３４】本発明による移動回転型モンテカルロ法の概
念を説明するための図。

【符号の説明】

Ｐ１…被検体の第１面、Ｐ２…被検体の第２面、Ｐ３…
被検体の第３面、Ｐ４…被検体の第４面、Ｐ５…被検体
の第５面、Ｐ６…被検体の第６面、１…光源、２…被検
体、３…光路、４…検出器、５…体積要素、６…光照射
位置、７…照射光、８…吸収物質、９…光源、１０…光
源用導波路、１１…分波器、１２…参照光用導波路、１
３…光源用光ファイバー、１４…光スイッチ、１５…光
照射用光ファイバー、１６…光ファイバー固定具、１７
…光検出用光ファイバー、１８…光スイッチ、１９…検
出器用光ファイバー、２０…光検出器、２１…信号伝送
ケーブル、２２…Ａ／Ｄ変換器、２３…信号伝送ケーブ
ル、２４…計算機、２５…表示装置、２６…光スイッチ
制御ケーブル、２７…光スイッチ制御ケーブル、２８…
光ファイバー固定具、２９…照射光、３０…検出光、３
１…注目領域、３２…注目領域外の領域、３３…光照射
位置、３４…光検出位置、３５…被検体の頂点、３６…
仮想散乱体、３７…光路、３８…仮想被検体。

フロントページの続き (72)発明者 伊藤 嘉敏 東京都国分寺市東恋ケ窪１丁目280番地 株式会社日立製作所中央研究所内 (72)発明者 山下 優一 東京都国分寺市東恋ケ窪１丁目280番地 株式会社日立製作所中央研究所内 (72)発明者 平林 由紀子 東京都国分寺市東恋ケ窪１丁目280番地 株式会社日立製作所中央研究所内 (72)発明者 小泉 英明 東京都国分寺市東恋ケ窪１丁目280番地 株式会社日立製作所中央研究所内 (72)発明者 川口 文男 東京都国分寺市東恋ケ窪１丁目280番地 株式会社日立製作所中央研究所内 (72)発明者 藤井 秀司 東京都国分寺市東恋ケ窪１丁目280番地 株式会社日立製作所中央研究所内 Ｆターム(参考） 2G059 AA01 AA05 BB02 BB12 CC16 EE02 EE11 FF01 FF02 FF08 GG01 GG08 JJ11 JJ17 JJ22 KK04 MM01 MM03 MM09 MM10 PP04

Claims (1)

【特許請求の範囲】
1. 【請求項１】光吸収物質を含んでなる光散乱体からなる
   被検体に、所定の光照射位置から所定の波長を有するパ
   ルス光又は連続光を照射する第１の光ファイバを含む照
   射手段と、前記被検体を通過する光の強度を所定の光検
   出位置で検出する第２の光ファイバを含む検出手段と、
   前記第１と第２の光ファイバを固定する固定具と、仮想
   被検体を想定して、前記光照射位置から前記光検出位置
   に至るまでの散乱光の光路を求め前記光路の座標から光
   路長を特定し記憶する記憶手段と、予め求められた前記
   光吸収物質の光学定数とから前記光吸収物質の濃度の空
   間分布を求める演算手段と、前記光吸収物質の濃度の空
   間分布を表示する表示手段とを有することを特徴とする
   吸収物質濃度空間分布画像化装置。