JP3889883B2 - 散乱吸収体の内部情報の計測方法及び装置 - Google Patents
散乱吸収体の内部情報の計測方法及び装置 Download PDFInfo
- Publication number
- JP3889883B2 JP3889883B2 JP20548698A JP20548698A JP3889883B2 JP 3889883 B2 JP3889883 B2 JP 3889883B2 JP 20548698 A JP20548698 A JP 20548698A JP 20548698 A JP20548698 A JP 20548698A JP 3889883 B2 JP3889883 B2 JP 3889883B2
- Authority
- JP
- Japan
- Prior art keywords
- light
- coefficient
- calculating
- absorption coefficient
- optical path
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Expired - Lifetime
Links
Images
Description
【０００１】
【発明の属する技術分野】
本発明は、輸送散乱係数や吸収係数、さらには吸収成分の濃度といった散乱吸収体の内部情報を計測するための方法並びに装置に関する。
【０００２】
【従来の技術】
マイクロ・ベア・ランバート則(Microscopic Beer-Lambert Law、以下「ＭＢＬ則」という)に基づいて散乱吸収体の吸収係数あるいは吸収成分の濃度を測定する方法としては、本発明者らによって特開平８−９４５１７号公報、特開平１０−７３４８１号公報、特開平１０−１１１２３８号公報に開示された方法がある。このようなＭＢＬ則に基づく方法は、原理的に、▲１▼媒体形状、▲２▼境界条件及び▲３▼散乱の影響を受けないという大きな特長がある。
【０００３】
しかしながら、上記の方法においては、吸収係数あるいは吸収成分の濃度を定量計測する際に、被計測対象の散乱特性がほぼ等しいかもしくは差があっても無視できる程度に小さいと見做せるような複数の波長の光を使用する必要があるため、これらの波長における散乱特性の差が大きくなると計測誤差が増大し、差が顕著に大きくなると計測不能となる。そのため、多波長分光測定を行う際に散乱特性に波長依存性がある場合には、別の方法で輸送散乱係数を測定あるいは推定する必要がある。ところが、従来の輸送散乱係数の測定方法においては、以下に説明するように演算時間が長いといった問題があり、運動体などについてのリアルタイム計測が困難であるといった点で未だ充分なものではない。
【０００４】
すなわち、一種類の波長の光を用いた一種類又は一回の計測で散乱吸収体（散乱媒体）の輸送散乱係数と吸収係数を定量する方法として、いわゆるフィッティング法が知られている。かかる従来のフィッティング法としては、散乱吸収体に時間幅が十分短いと見なせるパルス光を入射して、散乱吸収体内を散乱伝播した光を検出し、得られた光検出信号の波形を光拡散方程式（あるいはその解）でフィッティングして、誤差が最小になるような輸送散乱係数μ'sと吸収係数μaを定量値として求める方法がある。このような方法に関しては、例えば(1) M. Miwa, Y. Ueda, and B. Chance: Development of time resolved spectroscopy system for quantitative non-invasive tissue measurement, Proc. SPIE, 2389, 142-149 (1995)に記載されている。
【０００５】
この場合、計測される時間波形は媒体のインパルス応答と計測系のインパルス応答（装置関数）とのコンボリューションになる。従って、先ず、計測される時間波形から媒体のインパルス応答（真の波形）を求める必要があり、計測される時間波形を装置関数でデコンボリューションする方法が一般的に用いられている。ところが、この方法は演算時間が長いという欠点がある。さらに、フィッティング演算をするにあたっては輸送散乱係数μ'sと吸収係数μaを更新しながら再計算を繰り返す必要があるため、フィッティング法には長い演算時間が必要であるという問題がある。また、光検出信号の全域（全時間領域）をフィッティング演算に用いると、輸送散乱係数μ'sと吸収係数μaの定量精度が悪くなるという問題もある。このような問題に関しては、例えば(2) R. Cubeddu, A. Pifferi, P. Taroni, A. Torricelli, and G. Valentini: A solid tissue phantom for migration studies, Phys. Med. Biol. 42, 1971-1979 (1997)に記載されている。そこで従来は、フィッティング演算に用いる光検出信号の範囲を経験的に決めているが、輸送散乱係数μ'sと吸収係数μaの定量精度がこのフィッティング範囲に依存するという、精度に関する重大な問題がある。
【０００６】
また、従来の他の方法としては、散乱吸収体に強度変調光を入射し、検出光の中から所定周波数成分の信号を抽出して振幅、位相、変調度などを計測し、強度変調光の変調周波数を変化（掃引又は走査）させた後に再び前記の振幅、位相、変調度などを計測してそれらの変調周波数に対する変化を求める方法がある。この場合、複数の周波数成分を含む強度変調光を入射して、光検出信号から複数の周波数成分の信号を抽出するようにしてもよい。そして、前記で求めた各種パラメータの周波数依存性あるいはスペクトルを光拡散方程式（あるいはその解）でフィッティングして、誤差が最小になるような輸送散乱係数μ'sと吸収係数μaを定量値として求める。しかしながら、この方法においてはかなり広い周波数範囲に対して計測値を求めてフィッティングする必要があるため、複数の周波数について計測して解析するには時間がかかり、またフィッティング演算にも時間がかかるため、実用的でないという問題がある。
【０００７】
【発明が解決しようとする課題】
本発明は、上記従来の課題に鑑みてなされたものであり、上記従来のフィッティング法を用いることなく、散乱吸収体の輸送散乱係数μ'sや吸収係数μaを簡単にかつ高精度に定量することが可能な方法及び装置を提供することを目的とする。
【０００８】
【課題を解決するための手段】
本発明者は、上記目的を達成すべく鋭意研究した結果、任意の波長を有するパルス光を用いてインパルス応答（媒体の時間応答に対して時間幅が十分短いと見做せるパルス光入射に対する媒体の時間応答）を計測し、そのインパルス応答の時間分解波形から演算される平均光路長と光路長の分散（以下、場合により単に「分散」という）あるいは歪度等を用いることによって散乱吸収体の輸送散乱係数μ's及び吸収係数μaを簡単にかつ高精度に定量できることを見出し、本発明に到達した。なお、時間分解波形の時間積分値の吸収係数μaに対する１階の偏微分値が平均光路長、２階の偏微分値が分散、３階の偏微分値が歪度に対応する。
【０００９】
また、本発明者は、任意の周波数成分を有する強度変調光を用いて得られた光検出信号から所定周波数成分の信号を抽出し、その所定周波数成分の信号から演算される遅延値（群遅延又は位相遅延）とその吸収係数に対する１階あるいは２階の偏微分値とを用いることによっても散乱吸収体の輸送散乱係数μ's及び吸収係数μaを簡単にかつ高精度に定量できることを見出し、本発明に到達した。
【００１０】
すなわち、本発明の第１の散乱吸収体の内部情報の計測方法は、
所定波長を有するパルス光を散乱吸収体中に光入射位置から入射する光入射ステップと、
前記散乱吸収体内部を伝播した光を光検出位置で検出して光検出信号を取得する光検出ステップと、
前記光検出信号に基づいて、検出光強度の時間変化を示す波形データを取得する信号処理ステップと、
前記波形データに基づいて、前記検出光を構成する複数光子の平均光路長を演算する平均光路長演算ステップと、
前記波形データに基づいて、前記複数光子の光路長の分散を演算する分散演算ステップと、
前記平均光路長及び前記分散と輸送散乱係数及び吸収係数との間の所定の関係に基づいて、前記散乱吸収体の輸送散乱係数及び吸収係数を算出する輸送散乱係数及び吸収係数算出ステップと、
を含むことを特徴とする方法である。
【００１１】
また、本発明の第１の散乱吸収体の内部情報の計測装置は、
所定波長を有するパルス光を散乱吸収体中に光入射位置から入射する光入射手段と、
前記散乱吸収体内部を伝播した光を光検出位置で検出して光検出信号を取得する光検出手段と、
前記光検出信号に基づいて、検出光強度の時間変化を示す波形データを取得する信号処理手段と、
前記波形データに基づいて、前記検出光を構成する複数光子の平均光路長を演算する平均光路長演算手段と、
前記波形データに基づいて、前記複数光子の光路長の分散を演算する分散演算手段と、
前記平均光路長及び前記分散と輸送散乱係数及び吸収係数との間の所定の関係に基づいて、前記散乱吸収体の輸送散乱係数及び吸収係数を算出する輸送散乱係数及び吸収係数算出手段と、
を備えることを特徴とする装置である。
【００１２】
上記本発明の方法及び装置においては、任意の波長を有するパルス光を用いて得られた光検出信号から検出光強度の時間変化を示す波形データが取得され、その波形データから演算される平均光路長と光路長の分散とを用いて散乱吸収体の輸送散乱係数μ's及び吸収係数μaが直接的に求められる。そのため、本発明の方法及び装置においては、フィッティング法を用いることなく、散乱吸収体の輸送散乱係数μ's及び吸収係数μaが簡単にかつ高精度に定量されることとなり、演算時間が短縮される。
【００１３】
なお、前記輸送散乱係数及び吸収係数算出ステップ（輸送散乱係数及び吸収係数算出手段）において用いられる所定の関係は、散乱吸収体内の光子移動を拡散近似して求められた前記パルス光の入射に対する応答から解析的に求められた関係であることが好ましい。
【００１４】
より具体的には、前記輸送散乱係数及び吸収係数算出ステップ（輸送散乱係数及び吸収係数算出手段）において、以下の式（Ｉ）：
【数９】
及び式（II）：
【数１０】
［式（Ｉ）及び式（II）中、μ'sは輸送散乱係数、Ｌは平均光路長、σ2は光路長の分散、ρは光入射位置から光検出位置までの距離、μaは吸収係数をそれぞれ示す］
の連立関係に基づいて前記輸送散乱係数及び前記吸収係数を算出することが特に好ましい。
【００１５】
また、前記輸送散乱係数及び吸収係数算出ステップ（輸送散乱係数及び吸収係数算出手段）において、以下の式（I'）：
【数１１】
及び式（II'）：
【数１２】
［式（I'）及び式（II'）中、μaは吸収係数、Ｌは平均光路長、σ2は光路長の分散、ρは光入射位置から光検出位置までの距離をそれぞれ示す］
の連立関係に基づいて前記輸送散乱係数及び前記吸収係数を算出することが特に好ましい。
【００１６】
上記本発明の方法及び装置においては、前記光入射ステップ（光入射手段）において２種類以上の波長を有するパルス光を入射し、前記光検出ステップ（光検出手段）において前記２種類以上の波長に関してそれぞれ前記光検出信号を取得し、前記信号処理ステップ（信号処理手段）において前記２種類以上の波長に関してそれぞれ前記波形データを取得し、前記平均光路長演算ステップ（平均光路長演算手段）において前記２種類以上の波長に関してそれぞれ前記平均光路長を演算し、前記分散演算ステップ（分散演算手段）において前記２種類以上の波長に関してそれぞれ前記光路長の分散を演算し、前記輸送散乱係数及び吸収係数算出ステップ（輸送散乱係数及び吸収係数算出手段）において前記２種類以上の波長に関してそれぞれ前記輸送散乱係数を算出してもよい。そしてこの場合、本発明の方法及び装置は、前記輸送散乱係数に基づいて該輸送散乱係数の比率を演算する比率演算ステップ（比率演算手段）と、前記波形データに基づいて前記２種類以上の波長に関してそれぞれ前記検出光の光量を演算する光量演算ステップ（光量演算手段）と、前記輸送散乱係数の比率と前記光量と前記平均光路長と前記２種類以上の波長のパルス光に対する吸収成分の単位濃度当たりの吸収係数の差との間の所定の関係に基づいて前記吸収成分の濃度を演算する濃度演算ステップ（濃度演算手段）と、を更に含むことが好ましい。
【００１７】
本発明の方法及び装置においては、フィッティング法を用いることなく散乱吸収体の輸送散乱係数μ'sが簡単にかつ高精度に定量され、後述するように異なる波長の光に対する輸送散乱係数μ'sの比率が特に高精度で得られる。そのため、ＭＢＬ則から導かれる上記の多波長分光計測を本発明の方法及び装置に適用することによって、輸送散乱係数が未知である散乱吸収体に対しても吸収成分の濃度を高精度で計測することが可能となる。
【００１８】
また、本発明の第２の散乱吸収体の内部情報の計測方法は、
所定周波数成分を有する強度変調光を散乱吸収体中に光入射位置から入射する光入射ステップと、
前記散乱吸収体内部を伝播した光を光検出位置で検出して光検出信号を取得する光検出ステップと、
前記光検出信号から所定周波数成分の信号を抽出する信号処理ステップと、
前記所定周波数成分の信号に基づいて、前記所定周波数成分の信号の群遅延及び位相遅延のうちのいずれかの遅延値を演算する遅延値演算ステップと、
前記所定周波数成分の信号に基づいて、前記群遅延の吸収係数に対する偏微分値及び前記位相遅延の吸収係数に対する偏微分値のうちのいずれかの偏微分値を演算する偏微分値演算ステップと、
前記遅延値及び前記偏微分値と輸送散乱係数及び吸収係数との間の所定の関係に基づいて、前記散乱吸収体の輸送散乱係数及び吸収係数を算出する輸送散乱係数及び吸収係数算出ステップと、
を含むことを特徴とする方法である。
【００１９】
また、本発明の第２の散乱吸収体の内部情報の計測装置は、
所定周波数成分を有する強度変調光を散乱吸収体中に光入射位置から入射する光入射手段と、
前記散乱吸収体内部を伝播した光を光検出位置で検出して光検出信号を取得する光検出手段と、
前記光検出信号から所定周波数成分の信号を抽出する信号処理手段と、
前記所定周波数成分の信号に基づいて、前記所定周波数成分の信号の群遅延及び位相遅延のうちのいずれかの遅延値を演算する遅延値演算手段と、
前記所定周波数成分の信号に基づいて、前記群遅延の吸収係数に対する偏微分値及び前記位相遅延の吸収係数に対する偏微分値のうちのいずれかの偏微分値を演算する偏微分値演算手段と、
前記遅延値及び前記偏微分値と輸送散乱係数及び吸収係数との間の所定の関係に基づいて、前記散乱吸収体の輸送散乱係数及び吸収係数を算出する輸送散乱係数及び吸収係数算出手段と、
を備えることを特徴とする装置である。
【００２０】
上記本発明の方法及び装置においては、任意の周波数成分を有する強度変調光を用いて得られた光検出信号から所定周波数成分の信号が抽出され、その所定周波数成分の信号から演算される遅延値（群遅延又は位相遅延）とその吸収係数に対する偏微分値とを用いて散乱吸収体の輸送散乱係数μ's及び吸収係数μaが直接的に求められる。そのため、本発明の方法及び装置においては、フィッティング法を用いることなく、散乱吸収体の輸送散乱係数μ's及び吸収係数μaが簡単にかつ高精度に定量されることとなり、演算時間が短縮される。
【００２１】
なお、前記輸送散乱係数及び吸収係数算出ステップ（輸送散乱係数及び吸収係数算出手段）において用いられる所定の関係は、散乱吸収体内の光子移動を拡散近似して求められた前記強度変調光の入射に対する応答から解析的に求められた関係であることが好ましい。
【００２２】
【発明の実施の形態】
以下、図面を参照しつつ本発明の原理及び好適な実施形態について詳細に説明する。尚、図面中、同一又は相当部分には同一符号を付することとする。なお、散乱されつつ進行する光に関しては３次元座標を用いて考える必要があるが、以下では説明を簡単にするために場合により２次元座標を用いて説明する。この場合、３次元座標は円柱座標（ρ，θ，ｚ）を考えると便利であり、２次元では（ρ，ｚ）となる。先ず、本発明の原理について説明する。
【００２３】
［本発明の原理］
散乱吸収体内部をジグザグに伝播する光子の生存率は、ジグザグ光路長ｌ（エル）と吸収係数μaとの積の指数関数exp（−μaｌ）になる。このとき、散乱吸収体のインパルス応答ｈ(t)は時間因果関数になり、
【数１３】
と表される。ここで、μsとμaは非等方散乱係数と吸収係数、ｃは媒体中の光速度、ｔは飛行時間、ｌは飛行距離である。飛行時間ｔは時間分解計測によって計測することができる。
【００２４】
インパルス応答の時間積分Ｉ（μs,μa）は、
【数１４】
となる。上記(2b)式のＬ（μs,μa）＝ｃ＜ｔ＞は検出された光子の平均光路長を表す。また＜ｔ＞はインパルス応答波形の重心に相当するもので、インパルス応答の時間分解波形から計算することができる。
【００２５】
つぎに、上記平均光路長Ｌ（μs,μa）の吸収係数依存性を求めると、
【数１５】
となる。ここで、σ2は光路長の分散を表し、平均光路長Ｌ（μs,μa）をμaで偏微分して符号を変えたものに等しい。この分散は、インパルス応答の時間分解波形から計算することができる。同様にして、平均光路長Ｌ（μs,μa）のμaに対する２階偏微分をとれば、光路長分布の歪度に相当する値が得られ、この歪度もインパルス応答の時間分解波形から計算することができる。
【００２６】
なお、通常の計測系では、入射光パルスの時間幅は有限であり、増幅器や計数回路の帯域幅も有限である。この場合、実際の計測で得られる時間波形（検出光強度等の時間変化を示す波形データ）は、散乱吸収体のインパルス応答と計測システムのインパルス応答（装置関数ともよばれる）とのコンボリューションになる。したがって、計測で得られる時間波形から計測装置の特性の影響を取り除いて、真の散乱吸収体のインパルス応答を求める必要がある。この手段としては次のような２種の方法がある。第１の方法は、公知のいわゆるデコンボリューション演算を行い、デコンボリューション演算の結果として得られる波形から真の平均光路長、分散、歪度等を求める方法である。また、第２の方法は、計測系のインパルス応答（装置関数）の平均光路長、分散、歪度等、さらに観測波形から得られる平均光路長、分散、歪度等から、真の波形すなわち媒体のインパルス応答の平均光路長、分散、歪度等を直接求める方法である。
【００２７】
次に、強度変調光（photon density wave）を利用する方法（位相変調法とよばれる）について説明する。媒体の周波数応答を示すシステム関数Ｈ(ω)は、インパルス応答のフーリエ変換で表され、
【数１６】
となる。ここで、Ｒ及びＸはそれぞれ実部及び虚部であり、Ａ及びφはそれぞれ振幅及び位相遅れであり、これらはロックインアンプなどで計測することができる。
【００２８】
そして、
【数１７】
なる関係（コーシー−リーマン(Cauchy-Riemann)の関係式）が成立し、システム関数Ｈ(ω)が正則関数であることが認識される。さらに、(5a)式及び(5b)式から、
【数１８】
を導出することができる。すると、例えば(5c)式から、
【数１９】
が導出される。(6)式中の左辺は観測可能な値である。また、右辺第１項の被積分関数は２つの変調周波数ω1とω2における位相遅れφ1とφ2から近似的に求めることができる。また、右辺第２項は吸収係数がμa＝０のときの値である。そして、群遅延
【数２０】
は前述した平均光路長Ｌ（μs,μa）に相当するもので、ω<<ｃμaのときは
【数２１】
（位相遅延）と近似される。
【００２９】
ここで、前記と同様にして、上記群遅延の吸収係数依存性を求めると、
【数２２】
なる関係が得られる。この場合、(7)式右辺の偏微分項
【数２３】
は、３種の変調周波数を利用して計測することができる。この偏微分項は先に求めた光路長の分散に相当する。また、上記では(5c)式から(6)式を求めたが、同様の関係式を(5a)式、(5b)式及び(5d)式からそれぞれ求めることも可能である。
【００３０】
なお、位相変調法の場合、先に述べた真の波形を求めるためのデコンボリューション演算は、通常不要である。
【００３１】
他方、均一媒体に対する時間依存型の光拡散方程式は次式
【数２４】
で与えられる。
【００３２】
ここで、Ψは光子流動率（フルーエンスレート）[photon/mm2・s]、Ｄは拡散係数、ｒは光入射位置（光源）からの位置ベクトル、μaは吸収係数、μ'sは輸送散乱係数である。なお、Ψ／ｃは光子密度[photon/mm3]に相当する。また、拡散定数Ｄは、Ｄ＝１／３（μa＋μ's）である。また、散乱角の余弦の平均値をｇとしたときμ's＝μs（１−ｇ）である。この拡散方程式の解はグリーン関数（Green's function）であり、(8)式を円柱座標（ρ,ｚ,ｔ）で表したときの解は、
【数２５】
となる。ここで、ｒ2＝ρ2＋ｚ2である。ただし、ここでは２次元モデル化して考えているから、円柱座標（ρ,θ,ｚ）のθは省略してある。
【００３３】
また、媒体の境界における屈折率のミスマッチングによる反射を考慮したときの境界条件は、
【数２６】
となる。ここで、Ψ（ρ,ξ,ｔ）はｚ＝ξにある自由表面
【数２７】
の中の位置（ρ,ξ）における光子流動率、ｒdは屈折率のミスマッチングによる界面の反射係数である。
【００３４】
このとき、光子流動率Ψがゼロになる位置までの距離、つまり外挿距離ｚbは、
【数２８】
である。そして、このときの出力光子流Ｊは、
【数２９】
となる。つまり、出力光子流Ｊは、自由表面でのΨ、あるいはΨの傾き(gradient)を用いた２種の方法で表すことができる。そして、有限の大きさをもつ媒体に対するＪ（ρ,ξ,ｔ）は、パターソンら（Patterson et al.）が提案したミラー光源（正負）による光子流の重ね合わせで表される。この点は、(3) M. S. Patterson, B. Chance, and B. C. Wilson: Time resolved reflectance and transmittance for the noninvasive measurement of tissue optical properties, Applied Optics, 28 2331-2336 (1989)に記載されている。
【００３５】
上記出力光子流Ｊを構成する光子の平均光路長Ｌ（μs,μa）は、拡散方程式を解いて解析的に求めることができる。例えば、半無限空間と見なせる十分大きい散乱体に対して、コリメート光を入射する場合には、前述したパターソンらの方法を用いて次のように解析できる。
【００３６】
つまり、コリメート光を入射する反射型の時間分解計測（ξ＝０）で得られる出力光子流Ｊ（ρ,ｔ）は、
【数３０】
となる。ただし、
【数３１】
なる関係が成立している。すると、Ｊ（ρ,ｔ）の時間積分値Ｉ(ρ)は、
【数３２】
となる。よって、
【数３３】
なる関係が成立する。
【００３７】
また、平均光路長Ｌ（μs,μa）は、単にＬと略記して、
【数３４】
と表される。但し、上式中の第３辺を求めるときはμa<<μ's及びｚ0<<ρを仮定した。また、平均光路長Ｌの吸収係数依存性は、
【数３５】
となる。但し、上式中の第４辺を求めるときもμa<<μ's及びｚ0<<ρを仮定した。
【００３８】
すると、前述の(2b)式、(3)式、(17)式及び(18)式から観測可能値Ｌとσ2に対して、
【数３６】
なる連立方程式が導出される。コンピュータを用いてこの連立方程式を解けば、簡単にμaとμ'sが求められる。この場合、(17)式及び(18)式中の近似式でない精密式を用いることもできる。さらにこの連立方程式は解析的に解くことも可能であり、解は
【数３７】
となる。つまり、計測値から求められる平均光路長Ｌと分散σ2とを用いて、輸送散乱係数μ's及び吸収係数μaを直接的に定量することができる。なお、実用上は(20a)式のμaを省略してもよい。
【００３９】
なお、上記(8)式に示した光拡散方程式においてＤ＝１／３μ'sとしたときは、
【数３８】
なる関係が得られ、この関係から、
【数３９】
なる関係式が得られる。ここで、(20a)及び(20b)式と(22a)及び(22b)式を比較すると、(20a)及び(20b)式を近似したものが(22a)及び(22b)式である。
【００４０】
以上と同様にして、平均光路長Ｌ（μs,μa）のμaに対する２階の偏微分値、つまり歪度に相当する値に対しても、光拡散方程式から所定の関係、すなわち(17)式及び(18)式あるいは(19a)式及び(19b)式と同様な関係式が得られる。従って、上記で求めた３つの関係式｛平均光路長Ｌ（μs,μa）に関する関係式、分散σ2に関する関係式、及び歪度γに関する関係式｝のうちの任意の２つの関係式の連立関係から、上記と同様にして輸送散乱係数μ'sと吸収係数μaを求めることができる。
【００４１】
以上では、半無限空間と見做せる十分大きな散乱体に対してコリメート光を入射する反射型の時間分解計測について説明したが、スラブ状の散乱体に対する透過型の計測に対しても同様にして、関係式を導出することができる。つまり、未知数が２個（μ'sとμa）あるから、２連の連立方程式を解いてこれらの未知数を求めることができる。また、厚さがｄの媒体を透過型で計測する場合は、上記の関係式において、光入射・光検出位置間距離をρ＝ｄ、ｒ＝（ｄ−ｚ0）として近似することができる。通常はｚ0<<ｄであるから、ｒ≒ｄ＝ρとしてよい。
【００４２】
また、上記の平均光路長と分散の場合と同様にして、光子の自由行程１／μ'sに対して十分に大きい散乱吸収体に対して、前記群遅延
【数４０】
や、その偏微分値
【数４１】
を前記パターソンらが提案したモデルを用いて解析的に拡散方程式から求めることができる。したがって、前記の平均光路長、分散及び歪度の場合と同様にして、位相変調法の場合にも、計測値から求められる群遅延とその吸収係数に対する偏微分値とを用いて、例えば下記式：
【数４２】
の連立関係に基づいて輸送散乱係数μ's及び吸収係数μaを直接的に定量することができる。ただし、
【数４３】
である。
【００４３】
さらに、θ≒０の場合には、(23a)式及び(23b)式の連立方程式は
【数４４】
となるので、(25a)式及び(25b)式の連立方程式に基づいて輸送散乱係数μ's及び吸収係数μaを定量することができる。
【００４４】
また、上記群遅延は、ω<<ｃμaのときは位相遅延と近似されるため、計測値から求められる位相遅延とその吸収係数に対する偏微分値とを用いて輸送散乱係数μ's及び吸収係数μaを直接的に定量することができる。
【００４５】
なお、位相変調法の場合には前述したような４種類の関係式（(5a)〜(5d)）があり、これらに対して上記と同様の関係が成立する。
【００４６】
また、光拡散方程式によれば、Ｄ＝１／３（μa＋μ's）あるいはＤ＝１／３μ'sのいずれを用いてもμa<<μ'sである限り、比較的大きい媒体（ｚ0<<ρのとき）に関する光検出信号に対してかなり高い精度で解が得られる。しかし、光拡散方程式の解とモンテカルロシミュレーション（Monte Carlo simulation）の計算結果とを詳しく検討すると、図１に示すように、Ｄ＝１／３（μa＋μ's）あるいはＤ＝１／３μ'sのいずれを用いた場合でも、光検出信号の立ち上がり部分にわずかな不一致が見られる。この点は、例えば(4) M. Bassani, F. Martelli, G. Zaccanti, and D. Contini: Independence of the diffusion coefficient from absorption: experimental and numerical evidence, Opt. Lett. 22, 853-866 (1997)に記載されている。つまり、拡散近似を用いると、媒体内の光速度ｃと直進距離から計算した時間より少し早い時間に光検出信号が現れる。光拡散方程式に現れる光拡散定数Ｄは、最近ではＤ＝１／３μ'sであると考えられているが、本発明者は詳細な検討の結果、一般的には(20a)及び(20b)式を用いた方が精度の良い定量ができることを見出した。これは、光拡散方程式の解に関する前記の不一致点を考慮した場合、(20a)及び(20b)式で表されるμ'sやμaの方が良い精度が得られることを示唆している。
【００４７】
また、異なる二つの波長の光に対する輸送散乱係数μ's1とμ's2の比（ｋ＝μ's1／μ's2）を求める場合には、いずれのＤを用いてもほぼ同じ精度が得られる。従って、媒体形状、境界条件等に関するパラメータを含んでいないＭＢＬ則から導かれる関係式を用いる計測法に上記の輸送散乱係数μ's1とμ's2の比（ｋ＝μ's1／μ's2）を適用することによって、散乱係数に波長依存性がある場合であっても多波長計測の誤差増大を招くことなく、媒体形状、境界条件等に左右されずに精度の高い計測が可能となる。この点は実施形態の説明の際にも詳述する。
【００４８】
なお、以上に求めた関係は、有限の厚さのスラブに対する透過型の計測に対しても成立する。また、上記の関係は、特開平７−２０９１７７号公報に記載の方法のような、等方光を入射する場合にも成立する。
【００４９】
（第１実施形態）
図２〜図４を参照して本発明の好適な一実施形態である第１実施形態について説明する。図２には、被計測媒体である散乱吸収体１の輸送散乱係数μ's及び吸収係数μaを定量する内部情報計測装置２が示してある。先ず、内部情報計測装置２の構成について説明する。また、説明を簡潔にするため、ここでは散乱吸収体１として散乱特性及び吸収特性が一様な媒体を考える。
【００５０】
図２に示す装置は、光入射用の光ガイド３を備えており、光ガイド３の先端が散乱吸収体１の表面に配置されている。そして、光ガイド３には波長選択器４を介して光源５が光学的に接続されており、光源５から発せられたパルス光は波長選択器４で所定波長に波長選択され、光ガイド３を通して位置ｕjから散乱吸収体１中に入射される。このパルス光の時間幅は、光検出信号から平均光路長が導出できる程度に短いものであればよく、普通は１０ｐｓ〜１００ｎｓ程度の範囲で選択される。また、光の波長は計測対象に応じて適宜に選択される。一般に生体では、ヘモグロビンなどの吸収の関係から７００ｎｍ以上の光がよく使用される。光源５には、発光ダイオード、レーザーダイオード、Ｈｅ−Ｎｅレーザー、チタンサファイアレーザー等種々のものが使用できる。本実施形態において使用する光源５は単波長の光を発生するものであってもよいが、２波長以上の光を発生可能なものであってもよい。
【００５１】
また、図２に示す装置は、光検出用の光ガイド６を備えており、光ガイド６の先端が散乱吸収体１の表面に配置されている。そして、光ガイド６には光検出器７が光学的に接続されており、散乱吸収体１内部を散乱されつつ伝播した光は位置ｖkから光ガイド６を介して光検出器７に導かれ、光検出器７で受光信号が光検出信号（電気信号）に変換される。また、光検出器７及び光源５には信号処理部８が電気的に接続されており、信号処理部８で光検出信号に基づいて検出光強度の時間変化を示す波形データが取得される。更に、信号処理部８には演算処理部９が電気的に接続されており、演算処理部９で波形データに基づいて検出光を構成する複数光子の平均光路長並びに分散が演算され、これらの平均光路長及び分散に基づいて輸送散乱係数μ's及び吸収係数μaが定量される。
【００５２】
上記の光ガイド３、波長選択器４及び光源５が本発明にかかる光入射手段、光ガイド６及び光検出器７が本発明にかかる光検出手段、信号処理部８が本発明にかかる信号処理手段、演算処理部９が本発明にかかる平均光路長演算手段、分散演算手段、輸送散乱係数及び吸収係数算出手段をそれぞれ構成する。
【００５３】
なお、散乱吸収体１の表面における光ガイド３による光入射面及び光ガイド６による光検出面以外の場所は、光を吸収あるいは遮光する構造にしておくことが望ましい。また、散乱吸収体１の内部を散乱伝播した光が複数の波長の光を含む場合には、光検出器７と光ガイド６との間に波長選択フィルタ（図示せず）を適宜配置してもよい。
【００５４】
図３は、光検出器７、信号処理部８及び演算処理部９の好適な構成の一例を示す。図３に示す構成は、いわゆる時間相関光電子計数法と呼ばれる方法を用いて高速時間波形計測法を実施するための構成である。光検出器７として光電子増倍管（ＰＭＴ）を用いており、また、信号処理部８がコンスタント・ファラクション・ディスクリミネータ（ＣＦＤ）２１、時間−振幅変換器（ＴＡＣ）２２及びＡＤコンバータ（Ａ／Ｄ）２３で構成されている。そして、ＰＭＴ７の出力信号は、ＣＦＤ２１を介してＴＡＣ２２に導かれて時間に対応したアナログ電圧に変換され、更にＡＤコンバータ２３でデジタル信号に変換される。このデジタル信号は、検出光強度の時間変化を示す波形データに対応するものである。
【００５５】
図３に示す演算処理部９においては、光源５及び信号処理部８にＣＰＵ３０が電気的に接続されており、光入射に同期した光検出のタイミング等がＣＰＵ３０によって制御されると共に、信号処理部８から出力された波形データはＣＰＵ３０に導かれる。また、複数の波長を有する測定光を使用する場合は、入射光の波長もＣＰＵ３０によって制御される。具体的な手法としては、異なる波長の光を時分割で入射させて使用する手法と、異なる波長の光を同時に含む光を使用する手法とがある。具体的な波長選択手段としては、ミラーを用いた光ビーム切り換え器、フィルターを用いた波長切り換え器、光スイッチを用いた光切り換え器等がある。
【００５６】
図３に示す演算処理部９は、更に、オペレーティングシステム（ＯＳ）４１及び後で詳述する内部情報計測プログラム４２が記憶されたプログラムメモリ４０と、各種データファイルが記憶されるたデータファイルメモリ５０と、得られた内部情報を示すデータを記憶するデータメモリ６１と、作業用データを一時的に記憶する作業用メモリ６２と、データの入力を受け付けるキーボード７１及びマウス７２を備える入力装置７０と、得られたデータを出力するディスプレイ８１及びプリンタ８２を備える出力装置８０とを備えており、これらも電気的に接続されているＣＰＵ３０によって制御される。なお、上記のメモリは、コンピュータの内部メモリ（ハードディスク）であっても、フレキシブルディスクであってもよい。
【００５７】
データファイルメモリ５０には、内部情報計測プログラム４２の実行によって得られる波形データ、平均光路長、装置関数（計測系のインパルスレスポンス）、分散、歪度、輸送散乱係数、吸収係数等の諸データが記憶され、また、入力装置７０を用いて予め入力された計測条件や既知値も記憶される。このような入力データとしては、被計測媒体の形状、光入射位置、光検出位置、光入射・光検出位置間距離、計測に用いる光の波長、計測の種類（反射型、透過型など）、計測対象となる吸収成分の所定の波長における消光比などがある。
【００５８】
なお、光検出器７は光電子増倍管のほか、光電管、フォトダイオード、アバランシェフォトダイオード、ＰＩＮフォトダイオード等、あらゆる種類の光検出器を使用することができる。光検出器７の選択に際しては、使用される測定光の波長の光が検出できる分光感度特性をもっていれば良い。更に、光信号が微弱であるときは高感度あるいは高利得の光検出器を使用することが好ましい。また、上記の光入射用光ガイド３や光検出用光ガイド６の代わりに、光ファイバーやレンズなどを利用してもよい。
【００５９】
次に、図４に示す本発明の方法の一実施形態のフローチャート（図３に示す内部情報計測プログラム４２の処理を示すフローチャート）に基づいて以下に詳細に説明する。
【００６０】
図４に示すフローチャートにおいては、先ず、光源５で生成した所定波長のパルス光を光ガイド３を介して散乱吸収体１の光入射位置ｕjに入射し（Ｓ１１０）、散乱吸収体１内部で散乱されつつ伝播した光を光検出位置ｖkに設置した光ガイド６を介して光検出器７で検出する（Ｓ１２０）。
【００６１】
そして、検出された光に対応する光検出信号が光検出器７から発せられ、信号処理部８において検出光強度の時間変化を示す波形データに変換される（Ｓ１３０）。そして、必要に応じて波形データについてデコンボリューション処理を施すことによって計測系の特性の影響が除去された散乱吸収体１のインパルス応答を求める（Ｓ１４０）。このデコンボリューションに際して必要となる計測系のインパルスレスポンス（装置関数）は、予め測定してデータファイルメモリ５０に記憶しておく。この際、計測系のインパルスレスポンス（装置関数）は、図３に示した構成で、散乱吸収体１を取り除いて、光ガイド３の光出力端と光ガイド６の光入力端とを接触結合させて、測定する。従って、計測系のインパルスレスポンス（装置関数）には、光源のパルス幅や検出系の帯域幅などの影響が含まれている。
【００６２】
次に、このインパルス応答の時間分解波形に基づいて、前記検出光を構成する複数光子の平均光路長Ｌ並びにその光路長の分散σ2を演算する（Ｓ１５０及びＳ１５１）。なお、平均光路長は、十分時間幅の短いパルス光入射に対する光検出信号の時間波形の加重平均であるから、例えば前出の(2b)式の第２辺を計算して求められる。また、光路長の分散σ2は、例えば前出の(3)式の第２辺を計算して求められる。
【００６３】
そして、この平均光路長Ｌ及び分散σ2と散乱吸収体の輸送散乱係数μ's及び吸収係数μaとの間の所定の関係、例えば前述した(20a)式及び(20b)式あるいは(22a)式及び(22b)式に示した関係を利用して、散乱吸収体の輸送散乱係数μ's及び吸収係数μaを直接的に算出し（Ｓ１６０）、出力する（Ｓ１７０）。
【００６４】
なお、上記では、検出光強度の時間変化を示す波形データを計測系のインパルスレスポンス（装置関数）でデコンボリューション処理して、散乱吸収体１のインパルスレスポンスを求め（Ｓ１４０）、この波形から平均光路長及び分散を求めている（Ｓ１５０，Ｓ１５１）。しかし、この演算処理は次のような方法で求めることもできる。
【００６５】
すなわち、計測系のインパルスレスポンス（装置関数）から、平均光路長、分散（平均光路長の吸収係数μaに対する１階偏微分に相当）、歪度（平均光路長の吸収係数μaに対する２階偏微分に相当）などを求めておく。そして、検出光強度の時間変化を示す波形データから、平均光路長、分散、歪度などを求める。すると、散乱吸収体１のインパルスレスポンスに対する平均光路長、分散、歪度などは、検出光強度の時間変化を示す波形データから求めた平均光路長、分散、歪度などのそれぞれの値から計測系のインパルスレスポンス（装置関数）に対して求めた平均光路長、分散、歪度などのそれぞれの値を減算した値に等しい。この演算処理は、計測系のインパルスレスポンス（装置関数）及び検出光強度の時間変化を示す波形データを用いて、コンピュータで容易に計算することができ、その計算速度は前述したデコンボリューション法に比べて飛躍的に速くなる。このような演算処理を図４の中にＳ１４１で示した。
【００６６】
なお、以上の実施形態で、前記所定波長の光として複数波長のパルス光を使用して、それぞれの波長における吸収係数μaを求め、これらの値から吸収成分の濃度を定量することもできる。
【００６７】
（第２実施形態）
前記所定波長の光として複数波長の光を使用し、複数波長の光についてそれぞれ第１実施形態と同様の計測を行えば、複数波長の各光に対する輸送散乱係数μ'sや吸収係数μa、さらには複数波長の光に対する輸送散乱係数の比率や吸収係数の比率を計測することができる。本実施形態では、波長λ1とλ2の２つの光で計測する２波長分光計測法を本発明に適用して吸収成分の濃度を計測する方法について説明する。具体的には、１種類の光検出距離（光入射・光検出位置間距離）と２種類の光検出距離を用いるものがある。以下では、これらをまとめて説明する。
【００６８】
すなわち、本実施形態においては、２種類の波長（λ1、λ2）のパルス光を光入射位置から入射し、その光入射位置から距離ρ1およびρ2の光検出位置から光検出する。そして、実施形態１と同様の方法で、２種類の波長（λ1、λ2）及び２種類の光検出距離（ρ1、ρ2）に対応して４個の波形データを取得し、それぞれの波形データに基づいて平均光路長（＜Ｌ1（λ1）＞、＜Ｌ1（λ2）＞、＜Ｌ2（λ1）＞、＜Ｌ2（λ2）＞）及び分散（＜σ2 1（λ1）＞、＜σ2 1（λ2）＞、＜σ2 2（λ1）＞、＜σ2 2（λ2）＞）を演算し、前記２種類の波長（λ1、λ2）に関してそれぞれ輸送散乱係数（μ's1、μ's2）を算出する。
【００６９】
また、本実施形態においては、上記輸送散乱係数の比率（ｋ＝μ's2／μ's1）を演算すると共に、前記４個の波形データに基づいてそれぞれ検出光の光量（＜Ｉ1（λ1）＞、＜Ｉ1（λ2）＞、＜Ｉ2（λ1）＞、＜Ｉ2（λ2）＞）を演算する。
【００７０】
そして、上記の輸送散乱係数の比率、平均光路長及び光量を用いて、下記式：
【数４５】
又は下記式：
【数４６】
［上記両式中、
Ｖは吸収成分の濃度、
ε1は吸収成分の波長λ1の光に対する単位濃度当たりの吸収係数、
ε2は吸収成分の波長λ2の光に対する単位濃度当たりの吸収係数、
＜Ｌ1（λ1）＞は波長λ1の光に対して光検出距離ρ1で得た検出光強度の時間変化を示す波形データから得た平均光路長、
＜Ｌ1（λ2）＞は波長λ2の光に対して光検出距離ρ1で得た検出光強度の時間変化を示す波形データから得た平均光路長、
＜Ｌ2（λ1）＞は波長λ1の光に対して光検出距離ρ2で得た検出光強度の時間変化を示す波形データから得た平均光路長、
＜Ｌ2（λ2）＞は波長λ2の光に対して光検出距離ρ2で得た検出光強度の時間変化を示す波形データから得た平均光路長、
Ｉ1（λ1）は入射光強度Ｂ1、波長λ1の光に対する光検出距離ρ1での検出光量、
Ｉ1（λ2）は入射光強度Ｂ2、波長λ2の光に対する光検出距離ρ1での検出光量、
Ｉ2（λ1）は入射光強度Ｂ1、波長λ1の光に対する光検出距離ρ2での検出光量、
Ｉ2（λ2）は入射光強度Ｂ2、波長λ2の光に対する光検出距離ρ2での検出光量、
ｋは波長λ1の光に対する輸送散乱係数μ's1に対する波長λ2の光に対する輸送散乱係数μ's2の比（μ's2／μ's1）、
ｐは０≦ｐ≦１なる所定の値、
ｑは０≦ｑ≦１なる所定の値、
をそれぞれ示す］
に基づいて、吸収成分の濃度Ｖを演算する。このとき、経験的に定めることができる定数ｐとｑを用いる。また実際の計測では、ｐ＝ｑ＝１／２として、十分な精度が得られる。
【００７１】
（第３実施形態）
本実施形態は本発明を位相変調計測に応用する例を示す。この場合、内部情報計測装置の構成は、前述の図３に示す信号処理部８を、例えばロックインアンプを含む演算装置で置き換えた構成になる。また、光源５は３種の変調周波数成分（ω1，ω2，ω3）を含む変調光を発生する。
【００７２】
図５に、本発明の方法を位相変調計測に応用した実施形態のフローチャートを示す。図５に示すフローチャートにおいては、先ず、光源５で生成した所定波長の強度変調光を光ガイド３を介して散乱吸収体１の光入射位置ｕjに入射し（Ｓ１１０）、散乱吸収体１内部で散乱されつつ伝播した光を光検出位置ｖkに設置した光ガイド６を介して光検出器７で検出する（Ｓ１２０）。そして、検出された光に対応する光検出信号が光検出器７から発せられ、信号処理部８に供給される。
【００７３】
信号処理部８に含まれるロックインアンプは、３種の所定周波数成分の信号を抽出する（Ｓ１３１）と共に、３種の所定周波数成分の信号に対して(4)式に述べた実部Ｒ、虚部Ｘ、振幅Ａ、及び位相遅れφを出力する。次に、この実施形態では、上記の中の３種の所定周波数成分の信号に対する振幅Ａ、位相遅れφ及び３種の変調周波数成分（ω1，ω2，ω3）を用いて、変調周波数成分がω2である検出光を構成する複数光子の群遅延並びに群遅延の吸収係数に対する偏微分値を演算する（Ｓ１５２及びＳ１５３）。
【００７４】
そして、この群遅延及びその吸収係数に対する偏微分値と散乱吸収体の輸送散乱係数μ's及び吸収係数μaとの間の所定の関係、例えば前述した(23a)式及び(23b)式あるいは前述した(25a)式及び(25b)式に示した関係を利用して、散乱吸収体の輸送散乱係数μ's及び吸収係数μaを直接的に算出し（Ｓ１６０）、出力する（Ｓ１７０）。
【００７５】
なお、上記では群遅延及びその吸収係数に対する偏微分値を求めているが、前述のようにω<<ｃμaのときは群遅延は位相遅延と近似されることから、位相遅延並びに位相遅延の吸収係数に対する偏微分値を演算してもよい。
【００７６】
また、位相変調計測する場合にも、前記所定波長の光として複数波長の強度変調光を使用し、複数波長の強度変調光についてそれぞれ第３実施形態と同様の計測を行えば、複数波長の各光に対する輸送散乱係数μ'sや吸収係数μa、さらには複数波長の光に対する輸送散乱係数の比率や吸収係数の比率を計測することができ、第２実施形態のように吸収成分の濃度を計測することも可能である。
【００７７】
以上、本発明の好適な実施形態について説明したが、本発明は勿論上記実施形態に限定されるものではない。
【００７８】
すなわち、上記実施形態においては光入射位置及び光検出位置を固定しているが、光入射位置及び／又は光検出位置を走査させてもよい。また、散乱吸収体の周囲に複数の光入射位置及び／又は光検出位置を配置するようにしてもよい。
【００７９】
【実施例】
以下、実施例に基づいて本発明をより具体的に説明する。
【００８０】
実施例１
本実施例においては、本発明の方法の精度を確認するためにシミュレーションを行なった結果を示す。
【００８１】
すなわち、厚さ40mmのスラブ状の散乱吸収体に対して、以下の諸条件を設定してモンテカルロシミュレーションを行い、得られた輸送散乱係数（μ's1、μ's2）、吸収係数（μa1、μa2）、輸送散乱係数の比率（ｋ＝μ's1／μ's2）、及びその比率の誤差（（ｋ−μ'r）／μ'r）を表１及び表２に示す。なお、光入射位置と光検出位置との間の距離が15mmの反射型計測、並びに光入射位置と光検出位置との間の距離が40mm（媒体の厚さに相当する）の透過型計測の双方についてシミュレーションを行い、前者の結果を表１に、後者の結果を表２にそれぞれ示す。
【００８２】
（設定条件）
散乱係数：μS1＝4.5/mm、μS2＝3.0/mm
散乱角の余弦の平均値：ｇ＝0.85
μ'S1＝0.675、μ'S2＝0.45
境界での反射なし
入射光子数：2×108個
平均光路長Ｌの演算：(2b)式第2辺を利用
光路長の分散σ2の演算：(3)式第2辺を利用
輸送散乱係数μ'sの演算：(20a)式第3辺を利用
吸収係数μaの演算：(20b)式第3辺を利用。
【００８３】
上記条件下における輸送散乱係数の理論値はそれぞれμ's1＝0.675/mm及びμ's2＝0.45/mmである。従って、輸送散乱係数の比率の理論値（μ'r）はμ'r＝μ's1／μ's2＝0.675/0.45＝1.5である。なお、かかる理論値は、(5) H.J.van Staveren, C.J.M.Mooes, J.van Marle, S.S.Prahl, and M.C.J.van Gemert: Light scattering in Intralipid-10% in the wavelength range of 400-1000 nm, Applied Optics, 30, 4507-4514 (1991)に記載されている方法によって計算した値である。また、境界での屈折率差はない（matched-boundary condition）と仮定し、吸収係数がμa＝0.01、0.02及び0.03の３種類（輸送散乱係数とは独立）についてシミュレーションを行った。
【００８４】
【表１】
【表２】
表１及び表２に示した結果から明らかなように、本発明の方法による定量精度は高いことが認識された。また、吸収係数の定量精度、並びに光入射位置と光検出位置との間の距離が40mm程度以上であるときの輸送散乱係数及び輸送散乱係数の比率に対する定量精度が特に優れていることが判明した。
【００８５】
実施例２
本実施例においては、本発明の方法の精度を確認するためにファントムを用いて実験を行なった結果を示す。
【００８６】
すなわち、ファントムとしては、アクリル製の容器（幅120mm、高さ120mm、奥行き40mm）に散乱物質として１％イントラリピッド（Intralipid）溶液（Pharmacia AB社製、商品名：Intralipid）420mlを入れ、そこに吸収物質として グリーニッシュブラウン（greenish brown）インク（中外製薬社製）をインクの総量が0〜0.56mlになるように0.07ml間隔で添加したものを用いた。そして、各ファントムについて、以下の諸条件下で透過型の計測（光入射位置と光検出位置との間の距離：40mm）を行い、得られた輸送散乱係数（μ's1、μ's2）、吸収係数（μa1、μa2）、及び輸送散乱係数の比率（ｋ＝μ's1／μ's2）を表３に示す。なお、実際のファントムの吸収係数は、添加したインクの吸収係数と水の吸収係数との和になる。
【００８７】
（実験条件）
パルス光波長：λ1＝782nm、λ2＝831nm
パルス幅：50ps (FWHM)、繰り返し：5MHz、Peak Power：100mW
検出器：光電子増倍管
光入射ファイバー：200μmφ グレーテッドインデックス(GI)
光検出ファイバー：5mmφ バンドルファイバー
平均光路長Ｌの演算：(2b)式第2辺を利用
光路長の分散σ2の演算：(3)式第2辺を利用
輸送散乱係数μ'sの演算：(20a)式第3辺を利用
吸収係数μaの演算：(20b)式第3辺を利用。
【００８８】
上記条件下における輸送散乱係数の理論値は、λ1＝782nmのときμ's1＝1.0207、λ2＝831nmのときμ's2＝0.9531である。従って、輸送散乱係数の比率の理論値はμ's1／μ's2＝1.071である。なお、かかる理論値は前記文献(5)に記載されている方法によって計算した値である。
【００８９】
【表３】
表３に示した結果から明らかなように、本発明の方法による定量精度は高いことが認識された。また、輸送散乱係数の比率に対する定量精度が特に優れていることが判明した。
【００９０】
【発明の効果】
以上説明したように、本発明によれば、フィッティング法を用いることなく、散乱吸収体の輸送散乱係数μ'sや吸収係数μaを簡単かつ高精度に定量することが可能となる。また、本発明において複数の波長の光を利用することによって、複数波長の光の各々に対する輸送散乱係数μ'sや吸収係数μa、さらには複数波長間における輸送散乱係数の比率や吸収係数の比率を簡単かつ高精度に計測することが可能となる。
【図面の簡単な説明】
【図１】光拡散方程式の解とモンテカルロシミュレーションの計算結果との比較を示すグラフである。
【図２】本発明の散乱吸収体の内部情報計測装置の一実施形態を示す模式図である。
【図３】図２に示す装置の好適な具体的構成の一例を示す模式図である。
【図４】本発明の散乱吸収体の内部情報計測方法の一実施形態を示すフローチャートである。
【図５】本発明の散乱吸収体の内部情報計測方法の他の実施形態を示すフローチャートである。
【符号の説明】
１…散乱吸収体、２…内部情報計測装置、３…光ガイド、４…波長選択器、５…光源、６…光ガイド、７…光検出器、８…信号処理部、９…演算処理部。
【発明の属する技術分野】
本発明は、輸送散乱係数や吸収係数、さらには吸収成分の濃度といった散乱吸収体の内部情報を計測するための方法並びに装置に関する。
【０００２】
【従来の技術】
マイクロ・ベア・ランバート則(Microscopic Beer-Lambert Law、以下「ＭＢＬ則」という)に基づいて散乱吸収体の吸収係数あるいは吸収成分の濃度を測定する方法としては、本発明者らによって特開平８−９４５１７号公報、特開平１０−７３４８１号公報、特開平１０−１１１２３８号公報に開示された方法がある。このようなＭＢＬ則に基づく方法は、原理的に、▲１▼媒体形状、▲２▼境界条件及び▲３▼散乱の影響を受けないという大きな特長がある。
【０００３】
しかしながら、上記の方法においては、吸収係数あるいは吸収成分の濃度を定量計測する際に、被計測対象の散乱特性がほぼ等しいかもしくは差があっても無視できる程度に小さいと見做せるような複数の波長の光を使用する必要があるため、これらの波長における散乱特性の差が大きくなると計測誤差が増大し、差が顕著に大きくなると計測不能となる。そのため、多波長分光測定を行う際に散乱特性に波長依存性がある場合には、別の方法で輸送散乱係数を測定あるいは推定する必要がある。ところが、従来の輸送散乱係数の測定方法においては、以下に説明するように演算時間が長いといった問題があり、運動体などについてのリアルタイム計測が困難であるといった点で未だ充分なものではない。
【０００４】
すなわち、一種類の波長の光を用いた一種類又は一回の計測で散乱吸収体（散乱媒体）の輸送散乱係数と吸収係数を定量する方法として、いわゆるフィッティング法が知られている。かかる従来のフィッティング法としては、散乱吸収体に時間幅が十分短いと見なせるパルス光を入射して、散乱吸収体内を散乱伝播した光を検出し、得られた光検出信号の波形を光拡散方程式（あるいはその解）でフィッティングして、誤差が最小になるような輸送散乱係数μ'sと吸収係数μaを定量値として求める方法がある。このような方法に関しては、例えば(1) M. Miwa, Y. Ueda, and B. Chance: Development of time resolved spectroscopy system for quantitative non-invasive tissue measurement, Proc. SPIE, 2389, 142-149 (1995)に記載されている。
【０００５】
この場合、計測される時間波形は媒体のインパルス応答と計測系のインパルス応答（装置関数）とのコンボリューションになる。従って、先ず、計測される時間波形から媒体のインパルス応答（真の波形）を求める必要があり、計測される時間波形を装置関数でデコンボリューションする方法が一般的に用いられている。ところが、この方法は演算時間が長いという欠点がある。さらに、フィッティング演算をするにあたっては輸送散乱係数μ'sと吸収係数μaを更新しながら再計算を繰り返す必要があるため、フィッティング法には長い演算時間が必要であるという問題がある。また、光検出信号の全域（全時間領域）をフィッティング演算に用いると、輸送散乱係数μ'sと吸収係数μaの定量精度が悪くなるという問題もある。このような問題に関しては、例えば(2) R. Cubeddu, A. Pifferi, P. Taroni, A. Torricelli, and G. Valentini: A solid tissue phantom for migration studies, Phys. Med. Biol. 42, 1971-1979 (1997)に記載されている。そこで従来は、フィッティング演算に用いる光検出信号の範囲を経験的に決めているが、輸送散乱係数μ'sと吸収係数μaの定量精度がこのフィッティング範囲に依存するという、精度に関する重大な問題がある。
【０００６】
また、従来の他の方法としては、散乱吸収体に強度変調光を入射し、検出光の中から所定周波数成分の信号を抽出して振幅、位相、変調度などを計測し、強度変調光の変調周波数を変化（掃引又は走査）させた後に再び前記の振幅、位相、変調度などを計測してそれらの変調周波数に対する変化を求める方法がある。この場合、複数の周波数成分を含む強度変調光を入射して、光検出信号から複数の周波数成分の信号を抽出するようにしてもよい。そして、前記で求めた各種パラメータの周波数依存性あるいはスペクトルを光拡散方程式（あるいはその解）でフィッティングして、誤差が最小になるような輸送散乱係数μ'sと吸収係数μaを定量値として求める。しかしながら、この方法においてはかなり広い周波数範囲に対して計測値を求めてフィッティングする必要があるため、複数の周波数について計測して解析するには時間がかかり、またフィッティング演算にも時間がかかるため、実用的でないという問題がある。
【０００７】
【発明が解決しようとする課題】
本発明は、上記従来の課題に鑑みてなされたものであり、上記従来のフィッティング法を用いることなく、散乱吸収体の輸送散乱係数μ'sや吸収係数μaを簡単にかつ高精度に定量することが可能な方法及び装置を提供することを目的とする。
【０００８】
【課題を解決するための手段】
本発明者は、上記目的を達成すべく鋭意研究した結果、任意の波長を有するパルス光を用いてインパルス応答（媒体の時間応答に対して時間幅が十分短いと見做せるパルス光入射に対する媒体の時間応答）を計測し、そのインパルス応答の時間分解波形から演算される平均光路長と光路長の分散（以下、場合により単に「分散」という）あるいは歪度等を用いることによって散乱吸収体の輸送散乱係数μ's及び吸収係数μaを簡単にかつ高精度に定量できることを見出し、本発明に到達した。なお、時間分解波形の時間積分値の吸収係数μaに対する１階の偏微分値が平均光路長、２階の偏微分値が分散、３階の偏微分値が歪度に対応する。
【０００９】
また、本発明者は、任意の周波数成分を有する強度変調光を用いて得られた光検出信号から所定周波数成分の信号を抽出し、その所定周波数成分の信号から演算される遅延値（群遅延又は位相遅延）とその吸収係数に対する１階あるいは２階の偏微分値とを用いることによっても散乱吸収体の輸送散乱係数μ's及び吸収係数μaを簡単にかつ高精度に定量できることを見出し、本発明に到達した。
【００１０】
すなわち、本発明の第１の散乱吸収体の内部情報の計測方法は、
所定波長を有するパルス光を散乱吸収体中に光入射位置から入射する光入射ステップと、
前記散乱吸収体内部を伝播した光を光検出位置で検出して光検出信号を取得する光検出ステップと、
前記光検出信号に基づいて、検出光強度の時間変化を示す波形データを取得する信号処理ステップと、
前記波形データに基づいて、前記検出光を構成する複数光子の平均光路長を演算する平均光路長演算ステップと、
前記波形データに基づいて、前記複数光子の光路長の分散を演算する分散演算ステップと、
前記平均光路長及び前記分散と輸送散乱係数及び吸収係数との間の所定の関係に基づいて、前記散乱吸収体の輸送散乱係数及び吸収係数を算出する輸送散乱係数及び吸収係数算出ステップと、
を含むことを特徴とする方法である。
【００１１】
また、本発明の第１の散乱吸収体の内部情報の計測装置は、
所定波長を有するパルス光を散乱吸収体中に光入射位置から入射する光入射手段と、
前記散乱吸収体内部を伝播した光を光検出位置で検出して光検出信号を取得する光検出手段と、
前記光検出信号に基づいて、検出光強度の時間変化を示す波形データを取得する信号処理手段と、
前記波形データに基づいて、前記検出光を構成する複数光子の平均光路長を演算する平均光路長演算手段と、
前記波形データに基づいて、前記複数光子の光路長の分散を演算する分散演算手段と、
前記平均光路長及び前記分散と輸送散乱係数及び吸収係数との間の所定の関係に基づいて、前記散乱吸収体の輸送散乱係数及び吸収係数を算出する輸送散乱係数及び吸収係数算出手段と、
を備えることを特徴とする装置である。
【００１２】
上記本発明の方法及び装置においては、任意の波長を有するパルス光を用いて得られた光検出信号から検出光強度の時間変化を示す波形データが取得され、その波形データから演算される平均光路長と光路長の分散とを用いて散乱吸収体の輸送散乱係数μ's及び吸収係数μaが直接的に求められる。そのため、本発明の方法及び装置においては、フィッティング法を用いることなく、散乱吸収体の輸送散乱係数μ's及び吸収係数μaが簡単にかつ高精度に定量されることとなり、演算時間が短縮される。
【００１３】
なお、前記輸送散乱係数及び吸収係数算出ステップ（輸送散乱係数及び吸収係数算出手段）において用いられる所定の関係は、散乱吸収体内の光子移動を拡散近似して求められた前記パルス光の入射に対する応答から解析的に求められた関係であることが好ましい。
【００１４】
より具体的には、前記輸送散乱係数及び吸収係数算出ステップ（輸送散乱係数及び吸収係数算出手段）において、以下の式（Ｉ）：
【数９】
【数１０】
の連立関係に基づいて前記輸送散乱係数及び前記吸収係数を算出することが特に好ましい。
【００１５】
また、前記輸送散乱係数及び吸収係数算出ステップ（輸送散乱係数及び吸収係数算出手段）において、以下の式（I'）：
【数１１】
【数１２】
の連立関係に基づいて前記輸送散乱係数及び前記吸収係数を算出することが特に好ましい。
【００１６】
上記本発明の方法及び装置においては、前記光入射ステップ（光入射手段）において２種類以上の波長を有するパルス光を入射し、前記光検出ステップ（光検出手段）において前記２種類以上の波長に関してそれぞれ前記光検出信号を取得し、前記信号処理ステップ（信号処理手段）において前記２種類以上の波長に関してそれぞれ前記波形データを取得し、前記平均光路長演算ステップ（平均光路長演算手段）において前記２種類以上の波長に関してそれぞれ前記平均光路長を演算し、前記分散演算ステップ（分散演算手段）において前記２種類以上の波長に関してそれぞれ前記光路長の分散を演算し、前記輸送散乱係数及び吸収係数算出ステップ（輸送散乱係数及び吸収係数算出手段）において前記２種類以上の波長に関してそれぞれ前記輸送散乱係数を算出してもよい。そしてこの場合、本発明の方法及び装置は、前記輸送散乱係数に基づいて該輸送散乱係数の比率を演算する比率演算ステップ（比率演算手段）と、前記波形データに基づいて前記２種類以上の波長に関してそれぞれ前記検出光の光量を演算する光量演算ステップ（光量演算手段）と、前記輸送散乱係数の比率と前記光量と前記平均光路長と前記２種類以上の波長のパルス光に対する吸収成分の単位濃度当たりの吸収係数の差との間の所定の関係に基づいて前記吸収成分の濃度を演算する濃度演算ステップ（濃度演算手段）と、を更に含むことが好ましい。
【００１７】
本発明の方法及び装置においては、フィッティング法を用いることなく散乱吸収体の輸送散乱係数μ'sが簡単にかつ高精度に定量され、後述するように異なる波長の光に対する輸送散乱係数μ'sの比率が特に高精度で得られる。そのため、ＭＢＬ則から導かれる上記の多波長分光計測を本発明の方法及び装置に適用することによって、輸送散乱係数が未知である散乱吸収体に対しても吸収成分の濃度を高精度で計測することが可能となる。
【００１８】
また、本発明の第２の散乱吸収体の内部情報の計測方法は、
所定周波数成分を有する強度変調光を散乱吸収体中に光入射位置から入射する光入射ステップと、
前記散乱吸収体内部を伝播した光を光検出位置で検出して光検出信号を取得する光検出ステップと、
前記光検出信号から所定周波数成分の信号を抽出する信号処理ステップと、
前記所定周波数成分の信号に基づいて、前記所定周波数成分の信号の群遅延及び位相遅延のうちのいずれかの遅延値を演算する遅延値演算ステップと、
前記所定周波数成分の信号に基づいて、前記群遅延の吸収係数に対する偏微分値及び前記位相遅延の吸収係数に対する偏微分値のうちのいずれかの偏微分値を演算する偏微分値演算ステップと、
前記遅延値及び前記偏微分値と輸送散乱係数及び吸収係数との間の所定の関係に基づいて、前記散乱吸収体の輸送散乱係数及び吸収係数を算出する輸送散乱係数及び吸収係数算出ステップと、
を含むことを特徴とする方法である。
【００１９】
また、本発明の第２の散乱吸収体の内部情報の計測装置は、
所定周波数成分を有する強度変調光を散乱吸収体中に光入射位置から入射する光入射手段と、
前記散乱吸収体内部を伝播した光を光検出位置で検出して光検出信号を取得する光検出手段と、
前記光検出信号から所定周波数成分の信号を抽出する信号処理手段と、
前記所定周波数成分の信号に基づいて、前記所定周波数成分の信号の群遅延及び位相遅延のうちのいずれかの遅延値を演算する遅延値演算手段と、
前記所定周波数成分の信号に基づいて、前記群遅延の吸収係数に対する偏微分値及び前記位相遅延の吸収係数に対する偏微分値のうちのいずれかの偏微分値を演算する偏微分値演算手段と、
前記遅延値及び前記偏微分値と輸送散乱係数及び吸収係数との間の所定の関係に基づいて、前記散乱吸収体の輸送散乱係数及び吸収係数を算出する輸送散乱係数及び吸収係数算出手段と、
を備えることを特徴とする装置である。
【００２０】
上記本発明の方法及び装置においては、任意の周波数成分を有する強度変調光を用いて得られた光検出信号から所定周波数成分の信号が抽出され、その所定周波数成分の信号から演算される遅延値（群遅延又は位相遅延）とその吸収係数に対する偏微分値とを用いて散乱吸収体の輸送散乱係数μ's及び吸収係数μaが直接的に求められる。そのため、本発明の方法及び装置においては、フィッティング法を用いることなく、散乱吸収体の輸送散乱係数μ's及び吸収係数μaが簡単にかつ高精度に定量されることとなり、演算時間が短縮される。
【００２１】
なお、前記輸送散乱係数及び吸収係数算出ステップ（輸送散乱係数及び吸収係数算出手段）において用いられる所定の関係は、散乱吸収体内の光子移動を拡散近似して求められた前記強度変調光の入射に対する応答から解析的に求められた関係であることが好ましい。
【００２２】
【発明の実施の形態】
以下、図面を参照しつつ本発明の原理及び好適な実施形態について詳細に説明する。尚、図面中、同一又は相当部分には同一符号を付することとする。なお、散乱されつつ進行する光に関しては３次元座標を用いて考える必要があるが、以下では説明を簡単にするために場合により２次元座標を用いて説明する。この場合、３次元座標は円柱座標（ρ，θ，ｚ）を考えると便利であり、２次元では（ρ，ｚ）となる。先ず、本発明の原理について説明する。
【００２３】
［本発明の原理］
散乱吸収体内部をジグザグに伝播する光子の生存率は、ジグザグ光路長ｌ（エル）と吸収係数μaとの積の指数関数exp（−μaｌ）になる。このとき、散乱吸収体のインパルス応答ｈ(t)は時間因果関数になり、
【数１３】
【００２４】
インパルス応答の時間積分Ｉ（μs,μa）は、
【数１４】
【００２５】
つぎに、上記平均光路長Ｌ（μs,μa）の吸収係数依存性を求めると、
【数１５】
【００２６】
なお、通常の計測系では、入射光パルスの時間幅は有限であり、増幅器や計数回路の帯域幅も有限である。この場合、実際の計測で得られる時間波形（検出光強度等の時間変化を示す波形データ）は、散乱吸収体のインパルス応答と計測システムのインパルス応答（装置関数ともよばれる）とのコンボリューションになる。したがって、計測で得られる時間波形から計測装置の特性の影響を取り除いて、真の散乱吸収体のインパルス応答を求める必要がある。この手段としては次のような２種の方法がある。第１の方法は、公知のいわゆるデコンボリューション演算を行い、デコンボリューション演算の結果として得られる波形から真の平均光路長、分散、歪度等を求める方法である。また、第２の方法は、計測系のインパルス応答（装置関数）の平均光路長、分散、歪度等、さらに観測波形から得られる平均光路長、分散、歪度等から、真の波形すなわち媒体のインパルス応答の平均光路長、分散、歪度等を直接求める方法である。
【００２７】
次に、強度変調光（photon density wave）を利用する方法（位相変調法とよばれる）について説明する。媒体の周波数応答を示すシステム関数Ｈ(ω)は、インパルス応答のフーリエ変換で表され、
【数１６】
【００２８】
そして、
【数１７】
【数１８】
【数１９】
【数２０】
【数２１】
【００２９】
ここで、前記と同様にして、上記群遅延の吸収係数依存性を求めると、
【数２２】
【数２３】
【００３０】
なお、位相変調法の場合、先に述べた真の波形を求めるためのデコンボリューション演算は、通常不要である。
【００３１】
他方、均一媒体に対する時間依存型の光拡散方程式は次式
【数２４】
【００３２】
ここで、Ψは光子流動率（フルーエンスレート）[photon/mm2・s]、Ｄは拡散係数、ｒは光入射位置（光源）からの位置ベクトル、μaは吸収係数、μ'sは輸送散乱係数である。なお、Ψ／ｃは光子密度[photon/mm3]に相当する。また、拡散定数Ｄは、Ｄ＝１／３（μa＋μ's）である。また、散乱角の余弦の平均値をｇとしたときμ's＝μs（１−ｇ）である。この拡散方程式の解はグリーン関数（Green's function）であり、(8)式を円柱座標（ρ,ｚ,ｔ）で表したときの解は、
【数２５】
【００３３】
また、媒体の境界における屈折率のミスマッチングによる反射を考慮したときの境界条件は、
【数２６】
【数２７】
【００３４】
このとき、光子流動率Ψがゼロになる位置までの距離、つまり外挿距離ｚbは、
【数２８】
【数２９】
【００３５】
上記出力光子流Ｊを構成する光子の平均光路長Ｌ（μs,μa）は、拡散方程式を解いて解析的に求めることができる。例えば、半無限空間と見なせる十分大きい散乱体に対して、コリメート光を入射する場合には、前述したパターソンらの方法を用いて次のように解析できる。
【００３６】
つまり、コリメート光を入射する反射型の時間分解計測（ξ＝０）で得られる出力光子流Ｊ（ρ,ｔ）は、
【数３０】
【数３１】
【数３２】
【数３３】
【００３７】
また、平均光路長Ｌ（μs,μa）は、単にＬと略記して、
【数３４】
【数３５】
【００３８】
すると、前述の(2b)式、(3)式、(17)式及び(18)式から観測可能値Ｌとσ2に対して、
【数３６】
【数３７】
【００３９】
なお、上記(8)式に示した光拡散方程式においてＤ＝１／３μ'sとしたときは、
【数３８】
【数３９】
【００４０】
以上と同様にして、平均光路長Ｌ（μs,μa）のμaに対する２階の偏微分値、つまり歪度に相当する値に対しても、光拡散方程式から所定の関係、すなわち(17)式及び(18)式あるいは(19a)式及び(19b)式と同様な関係式が得られる。従って、上記で求めた３つの関係式｛平均光路長Ｌ（μs,μa）に関する関係式、分散σ2に関する関係式、及び歪度γに関する関係式｝のうちの任意の２つの関係式の連立関係から、上記と同様にして輸送散乱係数μ'sと吸収係数μaを求めることができる。
【００４１】
以上では、半無限空間と見做せる十分大きな散乱体に対してコリメート光を入射する反射型の時間分解計測について説明したが、スラブ状の散乱体に対する透過型の計測に対しても同様にして、関係式を導出することができる。つまり、未知数が２個（μ'sとμa）あるから、２連の連立方程式を解いてこれらの未知数を求めることができる。また、厚さがｄの媒体を透過型で計測する場合は、上記の関係式において、光入射・光検出位置間距離をρ＝ｄ、ｒ＝（ｄ−ｚ0）として近似することができる。通常はｚ0<<ｄであるから、ｒ≒ｄ＝ρとしてよい。
【００４２】
また、上記の平均光路長と分散の場合と同様にして、光子の自由行程１／μ'sに対して十分に大きい散乱吸収体に対して、前記群遅延
【数４０】
【数４１】
【数４２】
【数４３】
【００４３】
さらに、θ≒０の場合には、(23a)式及び(23b)式の連立方程式は
【数４４】
【００４４】
また、上記群遅延は、ω<<ｃμaのときは位相遅延と近似されるため、計測値から求められる位相遅延とその吸収係数に対する偏微分値とを用いて輸送散乱係数μ's及び吸収係数μaを直接的に定量することができる。
【００４５】
なお、位相変調法の場合には前述したような４種類の関係式（(5a)〜(5d)）があり、これらに対して上記と同様の関係が成立する。
【００４６】
また、光拡散方程式によれば、Ｄ＝１／３（μa＋μ's）あるいはＤ＝１／３μ'sのいずれを用いてもμa<<μ'sである限り、比較的大きい媒体（ｚ0<<ρのとき）に関する光検出信号に対してかなり高い精度で解が得られる。しかし、光拡散方程式の解とモンテカルロシミュレーション（Monte Carlo simulation）の計算結果とを詳しく検討すると、図１に示すように、Ｄ＝１／３（μa＋μ's）あるいはＤ＝１／３μ'sのいずれを用いた場合でも、光検出信号の立ち上がり部分にわずかな不一致が見られる。この点は、例えば(4) M. Bassani, F. Martelli, G. Zaccanti, and D. Contini: Independence of the diffusion coefficient from absorption: experimental and numerical evidence, Opt. Lett. 22, 853-866 (1997)に記載されている。つまり、拡散近似を用いると、媒体内の光速度ｃと直進距離から計算した時間より少し早い時間に光検出信号が現れる。光拡散方程式に現れる光拡散定数Ｄは、最近ではＤ＝１／３μ'sであると考えられているが、本発明者は詳細な検討の結果、一般的には(20a)及び(20b)式を用いた方が精度の良い定量ができることを見出した。これは、光拡散方程式の解に関する前記の不一致点を考慮した場合、(20a)及び(20b)式で表されるμ'sやμaの方が良い精度が得られることを示唆している。
【００４７】
また、異なる二つの波長の光に対する輸送散乱係数μ's1とμ's2の比（ｋ＝μ's1／μ's2）を求める場合には、いずれのＤを用いてもほぼ同じ精度が得られる。従って、媒体形状、境界条件等に関するパラメータを含んでいないＭＢＬ則から導かれる関係式を用いる計測法に上記の輸送散乱係数μ's1とμ's2の比（ｋ＝μ's1／μ's2）を適用することによって、散乱係数に波長依存性がある場合であっても多波長計測の誤差増大を招くことなく、媒体形状、境界条件等に左右されずに精度の高い計測が可能となる。この点は実施形態の説明の際にも詳述する。
【００４８】
なお、以上に求めた関係は、有限の厚さのスラブに対する透過型の計測に対しても成立する。また、上記の関係は、特開平７−２０９１７７号公報に記載の方法のような、等方光を入射する場合にも成立する。
【００４９】
（第１実施形態）
図２〜図４を参照して本発明の好適な一実施形態である第１実施形態について説明する。図２には、被計測媒体である散乱吸収体１の輸送散乱係数μ's及び吸収係数μaを定量する内部情報計測装置２が示してある。先ず、内部情報計測装置２の構成について説明する。また、説明を簡潔にするため、ここでは散乱吸収体１として散乱特性及び吸収特性が一様な媒体を考える。
【００５０】
図２に示す装置は、光入射用の光ガイド３を備えており、光ガイド３の先端が散乱吸収体１の表面に配置されている。そして、光ガイド３には波長選択器４を介して光源５が光学的に接続されており、光源５から発せられたパルス光は波長選択器４で所定波長に波長選択され、光ガイド３を通して位置ｕjから散乱吸収体１中に入射される。このパルス光の時間幅は、光検出信号から平均光路長が導出できる程度に短いものであればよく、普通は１０ｐｓ〜１００ｎｓ程度の範囲で選択される。また、光の波長は計測対象に応じて適宜に選択される。一般に生体では、ヘモグロビンなどの吸収の関係から７００ｎｍ以上の光がよく使用される。光源５には、発光ダイオード、レーザーダイオード、Ｈｅ−Ｎｅレーザー、チタンサファイアレーザー等種々のものが使用できる。本実施形態において使用する光源５は単波長の光を発生するものであってもよいが、２波長以上の光を発生可能なものであってもよい。
【００５１】
また、図２に示す装置は、光検出用の光ガイド６を備えており、光ガイド６の先端が散乱吸収体１の表面に配置されている。そして、光ガイド６には光検出器７が光学的に接続されており、散乱吸収体１内部を散乱されつつ伝播した光は位置ｖkから光ガイド６を介して光検出器７に導かれ、光検出器７で受光信号が光検出信号（電気信号）に変換される。また、光検出器７及び光源５には信号処理部８が電気的に接続されており、信号処理部８で光検出信号に基づいて検出光強度の時間変化を示す波形データが取得される。更に、信号処理部８には演算処理部９が電気的に接続されており、演算処理部９で波形データに基づいて検出光を構成する複数光子の平均光路長並びに分散が演算され、これらの平均光路長及び分散に基づいて輸送散乱係数μ's及び吸収係数μaが定量される。
【００５２】
上記の光ガイド３、波長選択器４及び光源５が本発明にかかる光入射手段、光ガイド６及び光検出器７が本発明にかかる光検出手段、信号処理部８が本発明にかかる信号処理手段、演算処理部９が本発明にかかる平均光路長演算手段、分散演算手段、輸送散乱係数及び吸収係数算出手段をそれぞれ構成する。
【００５３】
なお、散乱吸収体１の表面における光ガイド３による光入射面及び光ガイド６による光検出面以外の場所は、光を吸収あるいは遮光する構造にしておくことが望ましい。また、散乱吸収体１の内部を散乱伝播した光が複数の波長の光を含む場合には、光検出器７と光ガイド６との間に波長選択フィルタ（図示せず）を適宜配置してもよい。
【００５４】
図３は、光検出器７、信号処理部８及び演算処理部９の好適な構成の一例を示す。図３に示す構成は、いわゆる時間相関光電子計数法と呼ばれる方法を用いて高速時間波形計測法を実施するための構成である。光検出器７として光電子増倍管（ＰＭＴ）を用いており、また、信号処理部８がコンスタント・ファラクション・ディスクリミネータ（ＣＦＤ）２１、時間−振幅変換器（ＴＡＣ）２２及びＡＤコンバータ（Ａ／Ｄ）２３で構成されている。そして、ＰＭＴ７の出力信号は、ＣＦＤ２１を介してＴＡＣ２２に導かれて時間に対応したアナログ電圧に変換され、更にＡＤコンバータ２３でデジタル信号に変換される。このデジタル信号は、検出光強度の時間変化を示す波形データに対応するものである。
【００５５】
図３に示す演算処理部９においては、光源５及び信号処理部８にＣＰＵ３０が電気的に接続されており、光入射に同期した光検出のタイミング等がＣＰＵ３０によって制御されると共に、信号処理部８から出力された波形データはＣＰＵ３０に導かれる。また、複数の波長を有する測定光を使用する場合は、入射光の波長もＣＰＵ３０によって制御される。具体的な手法としては、異なる波長の光を時分割で入射させて使用する手法と、異なる波長の光を同時に含む光を使用する手法とがある。具体的な波長選択手段としては、ミラーを用いた光ビーム切り換え器、フィルターを用いた波長切り換え器、光スイッチを用いた光切り換え器等がある。
【００５６】
図３に示す演算処理部９は、更に、オペレーティングシステム（ＯＳ）４１及び後で詳述する内部情報計測プログラム４２が記憶されたプログラムメモリ４０と、各種データファイルが記憶されるたデータファイルメモリ５０と、得られた内部情報を示すデータを記憶するデータメモリ６１と、作業用データを一時的に記憶する作業用メモリ６２と、データの入力を受け付けるキーボード７１及びマウス７２を備える入力装置７０と、得られたデータを出力するディスプレイ８１及びプリンタ８２を備える出力装置８０とを備えており、これらも電気的に接続されているＣＰＵ３０によって制御される。なお、上記のメモリは、コンピュータの内部メモリ（ハードディスク）であっても、フレキシブルディスクであってもよい。
【００５７】
データファイルメモリ５０には、内部情報計測プログラム４２の実行によって得られる波形データ、平均光路長、装置関数（計測系のインパルスレスポンス）、分散、歪度、輸送散乱係数、吸収係数等の諸データが記憶され、また、入力装置７０を用いて予め入力された計測条件や既知値も記憶される。このような入力データとしては、被計測媒体の形状、光入射位置、光検出位置、光入射・光検出位置間距離、計測に用いる光の波長、計測の種類（反射型、透過型など）、計測対象となる吸収成分の所定の波長における消光比などがある。
【００５８】
なお、光検出器７は光電子増倍管のほか、光電管、フォトダイオード、アバランシェフォトダイオード、ＰＩＮフォトダイオード等、あらゆる種類の光検出器を使用することができる。光検出器７の選択に際しては、使用される測定光の波長の光が検出できる分光感度特性をもっていれば良い。更に、光信号が微弱であるときは高感度あるいは高利得の光検出器を使用することが好ましい。また、上記の光入射用光ガイド３や光検出用光ガイド６の代わりに、光ファイバーやレンズなどを利用してもよい。
【００５９】
次に、図４に示す本発明の方法の一実施形態のフローチャート（図３に示す内部情報計測プログラム４２の処理を示すフローチャート）に基づいて以下に詳細に説明する。
【００６０】
図４に示すフローチャートにおいては、先ず、光源５で生成した所定波長のパルス光を光ガイド３を介して散乱吸収体１の光入射位置ｕjに入射し（Ｓ１１０）、散乱吸収体１内部で散乱されつつ伝播した光を光検出位置ｖkに設置した光ガイド６を介して光検出器７で検出する（Ｓ１２０）。
【００６１】
そして、検出された光に対応する光検出信号が光検出器７から発せられ、信号処理部８において検出光強度の時間変化を示す波形データに変換される（Ｓ１３０）。そして、必要に応じて波形データについてデコンボリューション処理を施すことによって計測系の特性の影響が除去された散乱吸収体１のインパルス応答を求める（Ｓ１４０）。このデコンボリューションに際して必要となる計測系のインパルスレスポンス（装置関数）は、予め測定してデータファイルメモリ５０に記憶しておく。この際、計測系のインパルスレスポンス（装置関数）は、図３に示した構成で、散乱吸収体１を取り除いて、光ガイド３の光出力端と光ガイド６の光入力端とを接触結合させて、測定する。従って、計測系のインパルスレスポンス（装置関数）には、光源のパルス幅や検出系の帯域幅などの影響が含まれている。
【００６２】
次に、このインパルス応答の時間分解波形に基づいて、前記検出光を構成する複数光子の平均光路長Ｌ並びにその光路長の分散σ2を演算する（Ｓ１５０及びＳ１５１）。なお、平均光路長は、十分時間幅の短いパルス光入射に対する光検出信号の時間波形の加重平均であるから、例えば前出の(2b)式の第２辺を計算して求められる。また、光路長の分散σ2は、例えば前出の(3)式の第２辺を計算して求められる。
【００６３】
そして、この平均光路長Ｌ及び分散σ2と散乱吸収体の輸送散乱係数μ's及び吸収係数μaとの間の所定の関係、例えば前述した(20a)式及び(20b)式あるいは(22a)式及び(22b)式に示した関係を利用して、散乱吸収体の輸送散乱係数μ's及び吸収係数μaを直接的に算出し（Ｓ１６０）、出力する（Ｓ１７０）。
【００６４】
なお、上記では、検出光強度の時間変化を示す波形データを計測系のインパルスレスポンス（装置関数）でデコンボリューション処理して、散乱吸収体１のインパルスレスポンスを求め（Ｓ１４０）、この波形から平均光路長及び分散を求めている（Ｓ１５０，Ｓ１５１）。しかし、この演算処理は次のような方法で求めることもできる。
【００６５】
すなわち、計測系のインパルスレスポンス（装置関数）から、平均光路長、分散（平均光路長の吸収係数μaに対する１階偏微分に相当）、歪度（平均光路長の吸収係数μaに対する２階偏微分に相当）などを求めておく。そして、検出光強度の時間変化を示す波形データから、平均光路長、分散、歪度などを求める。すると、散乱吸収体１のインパルスレスポンスに対する平均光路長、分散、歪度などは、検出光強度の時間変化を示す波形データから求めた平均光路長、分散、歪度などのそれぞれの値から計測系のインパルスレスポンス（装置関数）に対して求めた平均光路長、分散、歪度などのそれぞれの値を減算した値に等しい。この演算処理は、計測系のインパルスレスポンス（装置関数）及び検出光強度の時間変化を示す波形データを用いて、コンピュータで容易に計算することができ、その計算速度は前述したデコンボリューション法に比べて飛躍的に速くなる。このような演算処理を図４の中にＳ１４１で示した。
【００６６】
なお、以上の実施形態で、前記所定波長の光として複数波長のパルス光を使用して、それぞれの波長における吸収係数μaを求め、これらの値から吸収成分の濃度を定量することもできる。
【００６７】
（第２実施形態）
前記所定波長の光として複数波長の光を使用し、複数波長の光についてそれぞれ第１実施形態と同様の計測を行えば、複数波長の各光に対する輸送散乱係数μ'sや吸収係数μa、さらには複数波長の光に対する輸送散乱係数の比率や吸収係数の比率を計測することができる。本実施形態では、波長λ1とλ2の２つの光で計測する２波長分光計測法を本発明に適用して吸収成分の濃度を計測する方法について説明する。具体的には、１種類の光検出距離（光入射・光検出位置間距離）と２種類の光検出距離を用いるものがある。以下では、これらをまとめて説明する。
【００６８】
すなわち、本実施形態においては、２種類の波長（λ1、λ2）のパルス光を光入射位置から入射し、その光入射位置から距離ρ1およびρ2の光検出位置から光検出する。そして、実施形態１と同様の方法で、２種類の波長（λ1、λ2）及び２種類の光検出距離（ρ1、ρ2）に対応して４個の波形データを取得し、それぞれの波形データに基づいて平均光路長（＜Ｌ1（λ1）＞、＜Ｌ1（λ2）＞、＜Ｌ2（λ1）＞、＜Ｌ2（λ2）＞）及び分散（＜σ2 1（λ1）＞、＜σ2 1（λ2）＞、＜σ2 2（λ1）＞、＜σ2 2（λ2）＞）を演算し、前記２種類の波長（λ1、λ2）に関してそれぞれ輸送散乱係数（μ's1、μ's2）を算出する。
【００６９】
また、本実施形態においては、上記輸送散乱係数の比率（ｋ＝μ's2／μ's1）を演算すると共に、前記４個の波形データに基づいてそれぞれ検出光の光量（＜Ｉ1（λ1）＞、＜Ｉ1（λ2）＞、＜Ｉ2（λ1）＞、＜Ｉ2（λ2）＞）を演算する。
【００７０】
そして、上記の輸送散乱係数の比率、平均光路長及び光量を用いて、下記式：
【数４５】
【数４６】
Ｖは吸収成分の濃度、
ε1は吸収成分の波長λ1の光に対する単位濃度当たりの吸収係数、
ε2は吸収成分の波長λ2の光に対する単位濃度当たりの吸収係数、
＜Ｌ1（λ1）＞は波長λ1の光に対して光検出距離ρ1で得た検出光強度の時間変化を示す波形データから得た平均光路長、
＜Ｌ1（λ2）＞は波長λ2の光に対して光検出距離ρ1で得た検出光強度の時間変化を示す波形データから得た平均光路長、
＜Ｌ2（λ1）＞は波長λ1の光に対して光検出距離ρ2で得た検出光強度の時間変化を示す波形データから得た平均光路長、
＜Ｌ2（λ2）＞は波長λ2の光に対して光検出距離ρ2で得た検出光強度の時間変化を示す波形データから得た平均光路長、
Ｉ1（λ1）は入射光強度Ｂ1、波長λ1の光に対する光検出距離ρ1での検出光量、
Ｉ1（λ2）は入射光強度Ｂ2、波長λ2の光に対する光検出距離ρ1での検出光量、
Ｉ2（λ1）は入射光強度Ｂ1、波長λ1の光に対する光検出距離ρ2での検出光量、
Ｉ2（λ2）は入射光強度Ｂ2、波長λ2の光に対する光検出距離ρ2での検出光量、
ｋは波長λ1の光に対する輸送散乱係数μ's1に対する波長λ2の光に対する輸送散乱係数μ's2の比（μ's2／μ's1）、
ｐは０≦ｐ≦１なる所定の値、
ｑは０≦ｑ≦１なる所定の値、
をそれぞれ示す］
に基づいて、吸収成分の濃度Ｖを演算する。このとき、経験的に定めることができる定数ｐとｑを用いる。また実際の計測では、ｐ＝ｑ＝１／２として、十分な精度が得られる。
【００７１】
（第３実施形態）
本実施形態は本発明を位相変調計測に応用する例を示す。この場合、内部情報計測装置の構成は、前述の図３に示す信号処理部８を、例えばロックインアンプを含む演算装置で置き換えた構成になる。また、光源５は３種の変調周波数成分（ω1，ω2，ω3）を含む変調光を発生する。
【００７２】
図５に、本発明の方法を位相変調計測に応用した実施形態のフローチャートを示す。図５に示すフローチャートにおいては、先ず、光源５で生成した所定波長の強度変調光を光ガイド３を介して散乱吸収体１の光入射位置ｕjに入射し（Ｓ１１０）、散乱吸収体１内部で散乱されつつ伝播した光を光検出位置ｖkに設置した光ガイド６を介して光検出器７で検出する（Ｓ１２０）。そして、検出された光に対応する光検出信号が光検出器７から発せられ、信号処理部８に供給される。
【００７３】
信号処理部８に含まれるロックインアンプは、３種の所定周波数成分の信号を抽出する（Ｓ１３１）と共に、３種の所定周波数成分の信号に対して(4)式に述べた実部Ｒ、虚部Ｘ、振幅Ａ、及び位相遅れφを出力する。次に、この実施形態では、上記の中の３種の所定周波数成分の信号に対する振幅Ａ、位相遅れφ及び３種の変調周波数成分（ω1，ω2，ω3）を用いて、変調周波数成分がω2である検出光を構成する複数光子の群遅延並びに群遅延の吸収係数に対する偏微分値を演算する（Ｓ１５２及びＳ１５３）。
【００７４】
そして、この群遅延及びその吸収係数に対する偏微分値と散乱吸収体の輸送散乱係数μ's及び吸収係数μaとの間の所定の関係、例えば前述した(23a)式及び(23b)式あるいは前述した(25a)式及び(25b)式に示した関係を利用して、散乱吸収体の輸送散乱係数μ's及び吸収係数μaを直接的に算出し（Ｓ１６０）、出力する（Ｓ１７０）。
【００７５】
なお、上記では群遅延及びその吸収係数に対する偏微分値を求めているが、前述のようにω<<ｃμaのときは群遅延は位相遅延と近似されることから、位相遅延並びに位相遅延の吸収係数に対する偏微分値を演算してもよい。
【００７６】
また、位相変調計測する場合にも、前記所定波長の光として複数波長の強度変調光を使用し、複数波長の強度変調光についてそれぞれ第３実施形態と同様の計測を行えば、複数波長の各光に対する輸送散乱係数μ'sや吸収係数μa、さらには複数波長の光に対する輸送散乱係数の比率や吸収係数の比率を計測することができ、第２実施形態のように吸収成分の濃度を計測することも可能である。
【００７７】
以上、本発明の好適な実施形態について説明したが、本発明は勿論上記実施形態に限定されるものではない。
【００７８】
すなわち、上記実施形態においては光入射位置及び光検出位置を固定しているが、光入射位置及び／又は光検出位置を走査させてもよい。また、散乱吸収体の周囲に複数の光入射位置及び／又は光検出位置を配置するようにしてもよい。
【００７９】
【実施例】
以下、実施例に基づいて本発明をより具体的に説明する。
【００８０】
実施例１
本実施例においては、本発明の方法の精度を確認するためにシミュレーションを行なった結果を示す。
【００８１】
すなわち、厚さ40mmのスラブ状の散乱吸収体に対して、以下の諸条件を設定してモンテカルロシミュレーションを行い、得られた輸送散乱係数（μ's1、μ's2）、吸収係数（μa1、μa2）、輸送散乱係数の比率（ｋ＝μ's1／μ's2）、及びその比率の誤差（（ｋ−μ'r）／μ'r）を表１及び表２に示す。なお、光入射位置と光検出位置との間の距離が15mmの反射型計測、並びに光入射位置と光検出位置との間の距離が40mm（媒体の厚さに相当する）の透過型計測の双方についてシミュレーションを行い、前者の結果を表１に、後者の結果を表２にそれぞれ示す。
【００８２】
（設定条件）
散乱係数：μS1＝4.5/mm、μS2＝3.0/mm
散乱角の余弦の平均値：ｇ＝0.85
μ'S1＝0.675、μ'S2＝0.45
境界での反射なし
入射光子数：2×108個
平均光路長Ｌの演算：(2b)式第2辺を利用
光路長の分散σ2の演算：(3)式第2辺を利用
輸送散乱係数μ'sの演算：(20a)式第3辺を利用
吸収係数μaの演算：(20b)式第3辺を利用。
【００８３】
上記条件下における輸送散乱係数の理論値はそれぞれμ's1＝0.675/mm及びμ's2＝0.45/mmである。従って、輸送散乱係数の比率の理論値（μ'r）はμ'r＝μ's1／μ's2＝0.675/0.45＝1.5である。なお、かかる理論値は、(5) H.J.van Staveren, C.J.M.Mooes, J.van Marle, S.S.Prahl, and M.C.J.van Gemert: Light scattering in Intralipid-10% in the wavelength range of 400-1000 nm, Applied Optics, 30, 4507-4514 (1991)に記載されている方法によって計算した値である。また、境界での屈折率差はない（matched-boundary condition）と仮定し、吸収係数がμa＝0.01、0.02及び0.03の３種類（輸送散乱係数とは独立）についてシミュレーションを行った。
【００８４】
【表１】
【００８５】
実施例２
本実施例においては、本発明の方法の精度を確認するためにファントムを用いて実験を行なった結果を示す。
【００８６】
すなわち、ファントムとしては、アクリル製の容器（幅120mm、高さ120mm、奥行き40mm）に散乱物質として１％イントラリピッド（Intralipid）溶液（Pharmacia AB社製、商品名：Intralipid）420mlを入れ、そこに吸収物質として グリーニッシュブラウン（greenish brown）インク（中外製薬社製）をインクの総量が0〜0.56mlになるように0.07ml間隔で添加したものを用いた。そして、各ファントムについて、以下の諸条件下で透過型の計測（光入射位置と光検出位置との間の距離：40mm）を行い、得られた輸送散乱係数（μ's1、μ's2）、吸収係数（μa1、μa2）、及び輸送散乱係数の比率（ｋ＝μ's1／μ's2）を表３に示す。なお、実際のファントムの吸収係数は、添加したインクの吸収係数と水の吸収係数との和になる。
【００８７】
（実験条件）
パルス光波長：λ1＝782nm、λ2＝831nm
パルス幅：50ps (FWHM)、繰り返し：5MHz、Peak Power：100mW
検出器：光電子増倍管
光入射ファイバー：200μmφ グレーテッドインデックス(GI)
光検出ファイバー：5mmφ バンドルファイバー
平均光路長Ｌの演算：(2b)式第2辺を利用
光路長の分散σ2の演算：(3)式第2辺を利用
輸送散乱係数μ'sの演算：(20a)式第3辺を利用
吸収係数μaの演算：(20b)式第3辺を利用。
【００８８】
上記条件下における輸送散乱係数の理論値は、λ1＝782nmのときμ's1＝1.0207、λ2＝831nmのときμ's2＝0.9531である。従って、輸送散乱係数の比率の理論値はμ's1／μ's2＝1.071である。なお、かかる理論値は前記文献(5)に記載されている方法によって計算した値である。
【００８９】
【表３】
【００９０】
【発明の効果】
以上説明したように、本発明によれば、フィッティング法を用いることなく、散乱吸収体の輸送散乱係数μ'sや吸収係数μaを簡単かつ高精度に定量することが可能となる。また、本発明において複数の波長の光を利用することによって、複数波長の光の各々に対する輸送散乱係数μ'sや吸収係数μa、さらには複数波長間における輸送散乱係数の比率や吸収係数の比率を簡単かつ高精度に計測することが可能となる。
【図面の簡単な説明】
【図１】光拡散方程式の解とモンテカルロシミュレーションの計算結果との比較を示すグラフである。
【図２】本発明の散乱吸収体の内部情報計測装置の一実施形態を示す模式図である。
【図３】図２に示す装置の好適な具体的構成の一例を示す模式図である。
【図４】本発明の散乱吸収体の内部情報計測方法の一実施形態を示すフローチャートである。
【図５】本発明の散乱吸収体の内部情報計測方法の他の実施形態を示すフローチャートである。
【符号の説明】
１…散乱吸収体、２…内部情報計測装置、３…光ガイド、４…波長選択器、５…光源、６…光ガイド、７…光検出器、８…信号処理部、９…演算処理部。
Claims (14)
- 所定波長を有するパルス光を散乱吸収体中に光入射位置から入射する光入射ステップと、
前記散乱吸収体内部を伝播した光を光検出位置で検出して光検出信号を取得する光検出ステップと、
前記光検出信号に基づいて、検出光強度の時間変化を示す波形データを取得する信号処理ステップと、
前記波形データに基づいて、前記検出光を構成する複数光子の平均光路長を演算する平均光路長演算ステップと、
前記波形データに基づいて、前記複数光子の光路長の分散を演算する分散演算ステップと、
前記平均光路長及び前記分散と輸送散乱係数及び吸収係数との間の所定の関係に基づいて、前記散乱吸収体の輸送散乱係数及び吸収係数を算出する輸送散乱係数及び吸収係数算出ステップと、
を含むことを特徴とする、散乱吸収体の内部情報の計測方法。 - 前記輸送散乱係数及び吸収係数算出ステップにおいて用いられる所定の関係は、散乱吸収体内の光子移動を拡散近似して求められた前記パルス光の入射に対する応答から解析的に求められた関係であることを特徴とする、請求項１記載の方法。
- 前記輸送散乱係数及び吸収係数算出ステップにおいて、以下の式（Ｉ）：
の連立関係に基づいて前記輸送散乱係数及び前記吸収係数を算出することを特徴とする、請求項１記載の方法。 - 前記輸送散乱係数及び吸収係数算出ステップにおいて、以下の式（I'）：
の連立関係に基づいて前記輸送散乱係数及び前記吸収係数を算出することを特徴とする、請求項１記載の方法。 - 前記光入射ステップにおいて、２種類以上の波長を有するパルス光を入射し、
前記光検出ステップにおいて、前記２種類以上の波長に関してそれぞれ前記光検出信号を取得し、
前記信号処理ステップにおいて、前記２種類以上の波長に関してそれぞれ前記波形データを取得し、
前記平均光路長演算ステップにおいて、前記２種類以上の波長に関してそれぞれ前記平均光路長を演算し、
前記分散演算ステップにおいて、前記２種類以上の波長に関してそれぞれ前記光路長の分散を演算し、
前記輸送散乱係数及び吸収係数算出ステップにおいて、前記２種類以上の波長に関してそれぞれ前記輸送散乱係数を算出すると共に、
前記輸送散乱係数に基づいて、該輸送散乱係数の比率を演算する比率演算ステップと、
前記波形データに基づいて、前記２種類以上の波長に関してそれぞれ前記検出光の光量を演算する光量演算ステップと、
前記輸送散乱係数の比率と、前記光量と、前記平均光路長と、前記２種類以上の波長のパルス光に対する吸収成分の単位濃度当たりの吸収係数の差との間の所定の関係に基づいて、前記吸収成分の濃度を演算する濃度演算ステップと、
を更に含むことを特徴とする、請求項１〜４のうちのいずれか一項記載の方法。 - 所定周波数成分を有する強度変調光を散乱吸収体中に光入射位置から入射する光入射ステップと、
前記散乱吸収体内部を伝播した光を光検出位置で検出して光検出信号を取得する光検出ステップと、
前記光検出信号から所定周波数成分の信号を抽出する信号処理ステップと、
前記所定周波数成分の信号に基づいて、前記所定周波数成分の信号の群遅延及び位相遅延のうちのいずれかの遅延値を演算する遅延値演算ステップと、
前記所定周波数成分の信号に基づいて、前記群遅延の吸収係数に対する偏微分値及び前記位相遅延の吸収係数に対する偏微分値のうちのいずれかの偏微分値を演算する偏微分値演算ステップと、
前記遅延値及び前記偏微分値と輸送散乱係数及び吸収係数との間の所定の関係に基づいて、前記散乱吸収体の輸送散乱係数及び吸収係数を算出する輸送散乱係数及び吸収係数算出ステップと、
を含むことを特徴とする、散乱吸収体の内部情報の計測方法。 - 前記輸送散乱係数及び吸収係数算出ステップにおいて用いられる所定の関係は、散乱吸収体内の光子移動を拡散近似して求められた前記強度変調光の入射に対する応答から解析的に求められた関係であることを特徴とする、請求項６記載の方法。
- 所定波長を有するパルス光を散乱吸収体中に光入射位置から入射する光入射手段と、
前記散乱吸収体内部を伝播した光を光検出位置で検出して光検出信号を取得する光検出手段と、
前記光検出信号に基づいて、検出光強度の時間変化を示す波形データを取得する信号処理手段と、
前記波形データに基づいて、前記検出光を構成する複数光子の平均光路長を演算する平均光路長演算手段と、
前記波形データに基づいて、前記複数光子の光路長の分散を演算する分散演算手段と、
前記平均光路長及び前記分散と輸送散乱係数及び吸収係数との間の所定の関係に基づいて、前記散乱吸収体の輸送散乱係数及び吸収係数を算出する輸送散乱係数及び吸収係数算出手段と、
を備えることを特徴とする、散乱吸収体の内部情報の計測装置。 - 前記輸送散乱係数及び吸収係数算出手段において用いられる所定の関係は、散乱吸収体内の光子移動を拡散近似して求められた前記パルス光の入射に対する応答から解析的に求められた関係であることを特徴とする、請求項８記載の装置。
- 前記輸送散乱係数及び吸収係数算出手段において、以下の式（Ｉ）：
の連立関係に基づいて前記輸送散乱係数及び前記吸収係数を算出することを特徴とする、請求項８記載の装置。 - 前記輸送散乱係数及び吸収係数算出手段において、以下の式（I'）：
の連立関係に基づいて前記輸送散乱係数及び前記吸収係数を算出することを特徴とする、請求項８記載の装置。 - 前記光入射手段において、２種類以上の波長を有するパルス光を入射し、
前記光検出手段において、前記２種類以上の波長に関してそれぞれ前記光検出信号を取得し、
前記信号処理手段において、前記２種類以上の波長に関してそれぞれ前記波形データを取得し、
前記平均光路長演算手段において、前記２種類以上の波長に関してそれぞれ前記平均光路長を演算し、
前記分散演算手段において、前記２種類以上の波長に関してそれぞれ前記光路長の分散を演算し、
前記輸送散乱係数及び吸収係数算出手段において、前記２種類以上の波長に関してそれぞれ前記輸送散乱係数を算出すると共に、
前記輸送散乱係数に基づいて、該輸送散乱係数の比率を演算する比率演算手段と、
前記波形データに基づいて、前記２種類以上の波長に関してそれぞれ前記検出光の光量を演算する光量演算手段と、
前記輸送散乱係数の比率と、前記光量と、前記平均光路長と、前記２種類以上の波長のパルス光に対する吸収成分の単位濃度当たりの吸収係数の差との間の所定の関係に基づいて、前記吸収成分の濃度を演算する濃度演算手段と、
を更に備えることを特徴とする、請求項８〜１１のうちのいずれか一項記載の装置。 - 所定周波数成分を有する強度変調光を散乱吸収体中に光入射位置から入射する光入射手段と、
前記散乱吸収体内部を伝播した光を光検出位置で検出して光検出信号を取得する光検出手段と、
前記光検出信号から所定周波数成分の信号を抽出する信号処理手段と、
前記所定周波数成分の信号に基づいて、前記所定周波数成分の信号の群遅延及び位相遅延のうちのいずれかの遅延値を演算する遅延値演算手段と、
前記所定周波数成分の信号に基づいて、前記群遅延の吸収係数に対する偏微分値及び前記位相遅延の吸収係数に対する偏微分値のうちのいずれかの偏微分値を演算する偏微分値演算手段と、
前記遅延値及び前記偏微分値と輸送散乱係数及び吸収係数との間の所定の関係に基づいて、前記散乱吸収体の輸送散乱係数及び吸収係数を算出する輸送散乱係数及び吸収係数算出手段と、
を備えることを特徴とする、散乱吸収体の内部情報の計測装置。 - 前記輸送散乱係数及び吸収係数算出手段において用いられる所定の関係は、散乱吸収体内の光子移動を拡散近似して求められた前記強度変調光の入射に対する応答から解析的に求められた関係であることを特徴とする、請求項１３記載の装置。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| JP20548698A JP3889883B2 (ja) | 1998-07-21 | 1998-07-21 | 散乱吸収体の内部情報の計測方法及び装置 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| JP20548698A JP3889883B2 (ja) | 1998-07-21 | 1998-07-21 | 散乱吸収体の内部情報の計測方法及び装置 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| JP2000039396A JP2000039396A (ja) | 2000-02-08 |
| JP3889883B2 true JP3889883B2 (ja) | 2007-03-07 |
Family
ID=16507660
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| JP20548698A Expired - Lifetime JP3889883B2 (ja) | 1998-07-21 | 1998-07-21 | 散乱吸収体の内部情報の計測方法及び装置 |
Country Status (1)
| Country | Link |
|---|---|
| JP (1) | JP3889883B2 (ja) |
Families Citing this family (3)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| JP4651186B2 (ja) * | 2000-12-11 | 2011-03-16 | 株式会社日立メディコ | 生体光計測装置 |
| EP2691901B1 (en) | 2011-03-28 | 2016-08-10 | AVL Test Systems, Inc. | Deconvolution method for emissions measurement |
| JP5834704B2 (ja) * | 2011-09-27 | 2015-12-24 | セイコーエプソン株式会社 | 濃度定量装置、光吸収係数算出方法、濃度定量方法、光吸収係数の算出を行うプログラム及び濃度の算出を行うプログラム |
Family Cites Families (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| JP3689496B2 (ja) * | 1995-08-08 | 2005-08-31 | 技術研究組合医療福祉機器研究所 | 計測装置 |
| JP3887668B2 (ja) * | 1996-03-28 | 2007-02-28 | 浜松ホトニクス株式会社 | 光散乱吸光媒質の吸光係数および光散乱係数の測定方法および測定装置 |
-
1998
- 1998-07-21 JP JP20548698A patent/JP3889883B2/ja not_active Expired - Lifetime
Also Published As
| Publication number | Publication date |
|---|---|
| JP2000039396A (ja) | 2000-02-08 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| Bizheva et al. | Path-length-resolved dynamic light scattering in highly scattering random media: The transition to diffusing wave spectroscopy | |
| JP3527600B2 (ja) | 散乱マトリックスの内部に関する分析データを決定するための方法および装置 | |
| US7920252B2 (en) | Method and apparatus for spectrophotometric characterization of turbid materials | |
| EP1081486A1 (en) | Method and device for measuring internal characteristic distribution of scattering/absorbing body | |
| EP0592200A1 (en) | Apparatus for measuring optical information in scattering medium and method therefor | |
| US20060243911A1 (en) | Measuring Technique | |
| JP3950243B2 (ja) | 散乱吸収体の内部情報の計測方法及び装置 | |
| AU769370B2 (en) | Method and apparatus for spectrometric analysis of turbid, pharmaceutical samples | |
| JPH07209177A (ja) | 散乱吸収体計測方法及び散乱吸収体計測装置 | |
| JP6630061B2 (ja) | 拡散スペクトルデータの処理方法及び処理装置 | |
| JP6101678B2 (ja) | 不透明な媒体の吸収係数を決定する方法、システムおよびコンピュータ・プログラム製品 | |
| Hebden et al. | The spatial resolution performance of a time‐resolved optical imaging system using temporal extrapolation | |
| JP4018799B2 (ja) | 散乱吸収体の吸収成分の濃度計測方法及び装置 | |
| EP1119763B1 (en) | Method for measuring locally and superficially the scattering and absorption properties of turbid media | |
| WO2009067043A1 (fr) | Procédé de mesure des dimensions de particules dans un liquide et dispositif de mise en oeuvre | |
| WO2007083691A1 (ja) | 光学分析装置 | |
| JP3889883B2 (ja) | 散乱吸収体の内部情報の計測方法及び装置 | |
| US7937226B2 (en) | Method and device for backscatter spectroscopy | |
| US20030078503A1 (en) | Method of calculating optical path distribution inside scattering absorber | |
| JP2003106979A (ja) | 低コヒーレンス干渉法を用いた動的光散乱測定装置 | |
| Helton et al. | Reconstruction of optical coefficients in turbid media using time-resolved reflectance and calibration-free instrument response functions | |
| Strojnik et al. | Interferometric tissue characterization: II. Experimental | |
| WO2023210521A1 (ja) | 光学特性値計測装置及び光学特性値計測方法 | |
| WO2010061673A1 (ja) | 散乱吸収体計測方法及び散乱吸収体計測装置 | |
| McMaster | Time-gated diffuse optical spectroscopy: experiments on layered media |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| A621 | Written request for application examination |
Free format text: JAPANESE INTERMEDIATE CODE: A621 Effective date: 20050330 |
|
| A977 | Report on retrieval |
Free format text: JAPANESE INTERMEDIATE CODE: A971007 Effective date: 20061027 |
|
| TRDD | Decision of grant or rejection written | ||
| A01 | Written decision to grant a patent or to grant a registration (utility model) |
Free format text: JAPANESE INTERMEDIATE CODE: A01 Effective date: 20061128 |
|
| A61 | First payment of annual fees (during grant procedure) |
Free format text: JAPANESE INTERMEDIATE CODE: A61 Effective date: 20061201 |
|
| R150 | Certificate of patent or registration of utility model |
Free format text: JAPANESE INTERMEDIATE CODE: R150 |
|
| S531 | Written request for registration of change of domicile |
Free format text: JAPANESE INTERMEDIATE CODE: R313532 |
|
| R350 | Written notification of registration of transfer |
Free format text: JAPANESE INTERMEDIATE CODE: R350 |
|
| FPAY | Renewal fee payment (event date is renewal date of database) |
Free format text: PAYMENT UNTIL: 20091208 Year of fee payment: 3 |
|
| FPAY | Renewal fee payment (event date is renewal date of database) |
Free format text: PAYMENT UNTIL: 20101208 Year of fee payment: 4 |
|
| FPAY | Renewal fee payment (event date is renewal date of database) |
Free format text: PAYMENT UNTIL: 20111208 Year of fee payment: 5 |
|
| FPAY | Renewal fee payment (event date is renewal date of database) |
Free format text: PAYMENT UNTIL: 20121208 Year of fee payment: 6 |
|
| FPAY | Renewal fee payment (event date is renewal date of database) |
Free format text: PAYMENT UNTIL: 20131208 Year of fee payment: 7 |
|
| EXPY | Cancellation because of completion of term |