JP5087526B2 - 散乱吸収体計測方法及び散乱吸収体計測装置 - Google Patents

Description

本発明は、生体等の散乱吸収体の内部情報を計測するための散乱吸収体計測方法及び散乱吸収体計測装置に関する。

生体等の散乱吸収体の内部情報を計測する装置として、光マンモグラフィ装置等の時間分解分光法（ＴＲＳ：Time Resolved Spectroscopy）を用いた散乱吸収体計測装置といった、散乱吸収体の内部を伝播したパルス光を用いる装置がある。このような装置で計測された光の時間分解波形には、計測対象とする散乱吸収体の時間応答特性の他、装置自体の時間応答特性（装置関数）が含まれている。従って、データ解析を行う際には、計測波形から装置関数の影響を除去する、即ち、データの校正を行う必要がある。

現状の光マンモグラフィ装置では、実際に計測した装置関数を用いてデータの校正を行っている。装置関数を計測する際、ガントリ部に設置した多数の照射用及び検出用の光ファイバを取り外して装置関数計測用モジュールに付け替え、装置関数の計測終了後に再び光ファイバにガントリ部に設置する必要がある。そのため、この方法では作業時間と計測時間とが長くなる点が問題となる。

この問題点を改善するために、装置関数を導出することなく計測対象とする散乱吸収体の計測波形を校正する、特許文献１に記載された方法が提案された。この方法では、計測対象とする散乱吸収体の計測波形及び光学係数（吸収係数及び散乱係数）が既知であるリファレンス用の散乱吸収体の計測波形とその理論波形とを用いて、計測対象とする散乱吸収体の計測波形の校正を行う。これによって、より簡易に装置関数を考慮した散乱吸収体の内部情報の取得が可能となった。
特開２００８−２９０８号公報

しかしながら、特許文献１に記載された上記の方法を用いると、校正処理を行う際に計測波形にノイズが影響し処理後の波形が大きく歪んでしまうことがあった。処理後の波形の歪みは、その後に算出される内部情報の正確性にも影響を及ぼしえる。このような状況から、計測波形にノイズが含まれる場合であっても、より正確かつ簡易に装置関数を考慮して散乱吸収体の内部情報を取得する方法が求められていた。

本発明は、以上の問題点を解決するためになされたものであり、計測波形にノイズが含まれる場合であっても、より正確かつ簡易に装置関数を考慮して散乱吸収体の内部情報を取得することができる散乱吸収体計測方法及び散乱吸収体計測装置を提供することを目的とする。

上記目的を達成するために、本発明に係る散乱吸収体計測方法は、計測対象の散乱吸収体及びリファレンス用の散乱吸収体に対して、所定波長のパルス光を光入射位置から入射する光入射ステップと、光入射ステップにおいて入射されて各散乱吸収体の内部を伝播したパルス光を光検出位置で検出して光検出信号を取得する光検出ステップと、光検出ステップにおいて検出された光検出信号に基づいて、光強度の時間変化を示す計測波形を取得する信号処理ステップと、信号処理ステップにおいて取得された計測対象の散乱吸収体での計測波形をＰ_ｍ（ｔ）、所定のパラメータを有する予め用意された関数によって示される計測対象の散乱吸収体での理論波形をＰ_ｍ ^{ｔｈｅｏｒｙ}（ｔ）、信号処理ステップにおいて取得されたリファレンス用の散乱吸収体での計測波形をＰ_ｒ（ｔ）、予め用意されたリファレンス用の散乱吸収体での理論波形をＰ_ｒ ^{ｔｈｅｏｒｙ}（ｔ）、コンボリューションを演算子“＊”とそれぞれしたときに、以下の式

に基づいてパラメータを特定する波形処理ステップと、波形処理ステップにおいて特定されたパラメータを有する関数によって示される計測対象の散乱吸収体での理論波形に基づいて、計測対象の散乱吸収体の内部情報を算出する内部情報算出ステップと、を含むことを特徴とする。

計測対象及びリファレンス用の散乱吸収体の計測波形には装置関数の影響が含まれている。従って、上記のような波形の関係に基づいて上記のパラメータを特定することによって、装置関数を算出することなく、計測対象の散乱吸収体の計測波形から装置関数の影響を除去した理論波形を得ることができる。また、計測対象の散乱吸収体の理論波形を所定のパラメータを有する予め用意された関数とすることによって、計測波形のばらつきの影響を防止することができる。即ち、計測波形にノイズが含まれる場合であっても、計測対象の散乱吸収体の理論波形をより正確に求めることができる。従って、本発明によれば、計測波形にノイズが含まれる場合であっても、より正確かつ簡易に装置関数を考慮して散乱吸収体の内部情報を取得することができる。

Ｐ_ｍ ^{ｔｈｅｏｒｙ}（ｔ）を示す関数は、離散的な時刻ｔ_ｎの関数であり、Ａ、α、βをパラメータ、ｃを光速を示す値、γを光入射位置と光検出位置との間の距離を示す値とそれぞれしたときに、以下の式

（ここで、

である）で表されることが望ましい。上記の関数は、散乱吸収体内部の光伝播をより適切に示すものであるので、更に正確に計測対象の散乱吸収体の理論波形を求めることができ、結果として、更に正確に散乱吸収体の内部情報を取得することができる。

ところで、本発明は、上記のように散乱吸収体計測方法の発明として記述できる他に、以下のように散乱吸収体計測装置の発明としても記述することができる。これはカテゴリが異なるだけで、実質的に同一の発明であり、同様の作用及び効果を奏する。

即ち、本発明に係る散乱吸収体計測装置は、計測対象の散乱吸収体及びリファレンス用の散乱吸収体に対して、所定波長のパルス光を光入射位置から入射する光入射手段と、光入射手段によって入射されて各散乱吸収体の内部を伝播したパルス光を光検出位置で検出して光検出信号を取得する光検出手段と、光検出手段によって検出された光検出信号に基づいて、光強度の時間変化を示す計測波形を取得する信号処理手段と、信号処理手段によって取得された計測対象の散乱吸収体での計測波形をＰ_ｍ（ｔ）、所定のパラメータを有する予め用意された関数によって示される計測対象の散乱吸収体での理論波形をＰ_ｍ ^{ｔｈｅｏｒｙ}（ｔ）、信号処理手段によって取得されたリファレンス用の散乱吸収体での計測波形をＰ_ｒ（ｔ）、予め用意されたリファレンス用の散乱吸収体での理論波形をＰ_ｒ ^{ｔｈｅｏｒｙ}（ｔ）、コンボリューションを演算子“＊”とそれぞれしたときに、以下の式

に基づいてパラメータを特定する波形処理手段と、波形処理手段によって特定されたパラメータを有する関数によって示される計測対象の散乱吸収体での理論波形に基づいて、計測対象の散乱吸収体の内部情報を算出する内部情報算出手段と、を備えることを特徴とする。

本発明によれば、装置関数を算出することなく、計測対象の散乱吸収体の計測波形から装置関数の影響を除去した理論波形を得ることができる。また、計測対象の散乱吸収体の理論波形を所定のパラメータを有する予め用意された関数とすることによって、計測波形にノイズが含まれる場合であっても、計測対象の散乱吸収体の理論波形をより正確に求めることができる。従って、本発明によれば、計測波形にノイズが含まれる場合であっても、より正確かつ簡易に装置関数を考慮して散乱吸収体の内部情報を取得することができる。

以下、図面とともに本発明に係る散乱吸収体計測方法及び散乱吸収体計測装置の好適な実施形態について詳細に説明する。なお、図面の説明においては同一要素には同一符号を付し、重複する説明を省略する。また、図面の寸法比率は、説明のものと必ずしも一致していない。

図１に本実施形態に係る散乱吸収体計測装置１０の機能構成を示す。図２に本実施形態に係る散乱吸収体計測装置１０のハードウェア構成を示す。散乱吸収体計測装置１０は、散乱吸収体の内部情報を計測するための装置である。散乱吸収体計測装置１０の計測対象である散乱吸収体（Scattering Medium）は、自身に入射された光を散乱及び吸収するものであり、例えば生体である。また、計測される内部情報は、例えば、散乱吸収体の散乱係数及び吸収係数、並びに散乱吸収体に含まれている各成分の濃度である。本実施形態に係る散乱吸収体計測装置１０は、これらの内部情報を計測することによって、例えば、乳癌の診断に応用される。なお、図２に示すハードウェア構成は、いわゆる時間相関光電子計数法と呼ばれる方法を用いて高速時間波形計測法を実施するための構成である。

図１及び２に示すように、散乱吸収体計測装置１０は、光源２０と、波長選択器２１と、光検出器３０と、光入射用の光ガイド４０ａと、光検出用の光ガイド４０ｂと、信号処理部５０と、演算処理部６０とを備えて構成されている。

光源２０は、散乱吸収体ＳＭに入射するための、所定波長のパルス光を発生させる装置である。光源２０により発生し散乱吸収体ＳＭに入射されるパルス光としては、散乱吸収体ＳＭの内部情報が計測できる程度に短い時間幅のパルス光が用いられ、通常は１ｎｓ以下の範囲の時間幅が選択される。光源２０から発生したパルス光は、波長選択器２１に入力される。光源２０としては、発光ダイオード、レーザダイオード、及び各種のパルスダイオード等、様々なものを使用することができる。

波長選択器２１は、光源２０から入力されたパルス光を、所定波長に波長選択する装置である。光源２０により発生し波長選択器２１により波長選択されるパルス光の波長は、計測対象である散乱吸収体ＳＭに応じて適宜、選択されるが、一般に例えば生体では、生体の透過率と定量すべき吸収成分の分光吸収係数との関係等から、７００〜９００ｎｍ程度の近赤外線領域の波長が用いられる。波長選択器２１により波長選択されたパルス光は、光入射用の光ガイド４０ａの入力端に入力される。なお、複数の成分についての内部情報を計測する場合等、必要があれば光源２０及び波長選択器２１は、複数の波長成分のパルス光を計測光として入射可能に構成される。また、波長選択を行う必要がない場合等には、波長選択器２１については、設置を省略してもよい。

光入射用の光ガイド４０ａは、その入力端で波長選択器２１からパルス光の入力を受け付けて、その出力端から散乱吸収体ＳＭに対して当該パルス光を入射するものである。光入射用の光ガイド４０ａは、その出力端が散乱吸収体ＳＭの表面上の光入射位置Ａに位置するように設置されている。光入射用の光ガイド４０ａは、具体的には、パルス光を伝達することができる光ファイバにより構成される。上記の光源２０、波長選択器２１及び光入射用の光ガイド４０ａは、散乱吸収体ＳＭに対して、所定波長のパルス光を光入射位置Ａから入射する入射手段を構成する。

光検出用の光ガイド４０ｂは、その入力端で散乱吸収体ＳＭの内部を伝播した光を入力させて、その出力端から光検出器３０に当該光を出力するものである。光検出用の光ガイド４０ｂは、その入力端が散乱吸収体ＳＭの表面上の光検出位置Ｂに位置するように設置されている。光検出用の光ガイド４０ｂは、具体的には、光を伝達することができる光ファイバにより構成される。

光検出器３０は、光検出用の光ガイド４０ｂから入力された光を検出するための装置である。光検出器３０は、検出した光の光強度等を示す光検出信号Ｓ１を生成する。生成された光検出信号Ｓ１は信号処理部５０に入力される。光検出器３０としては、図３に示す光電子増倍管（ＰＭＴ：Photomultiplier）の他、フォトダイオード、アバランシェフォトダイオード、ＰＩＮフォトダイオード等、様々なものを使用することができる。光検出器３０の選択については、散乱吸収体計測に使用されるパルス光の波長の光が十分に検出できる分光感度特性を有していればよい。また、光信号が微弱であるときは、高感度あるいは高利得の光検出器を使用することが好ましい。上記の光検出器３０及び光検出用の光ガイド４０ｂは、光入射手段により入射されて散乱吸収体ＳＭの内部を伝播したパルス光を光検出位置Ｂで検出して光検出信号を取得する光検出手段を構成する。

上述した光入射用の光ガイド４０ａ及び光検出用の光ガイド４０ｂは、それぞれ必ずしも一つでなくてもよい。例えば、乳癌の診断のために用いられる散乱吸収体計測装置１０では、計測対象の散乱吸収体ＳＭである乳房を入れることができるカップ状のガントリに複数のフォルダが設けられており、当該フォルダそれぞれに光ガイド４０ａ，４０ｂが差し込まれる。なお、各フォルダに差し込まれるのは、入射用の光ガイド４０ａと光検出用の光ガイド４０ｂとが１本になった光ガイドとして、演算処理部６０からの制御等により何れのフォルダ位置からでも入射及び検出の何れかができるようにしておくのがよい。

信号処理部５０は、光検出器３０及び光源２０と電気的に接続されており、光検出器３０により検出された光検出信号及び光源２０からのパルス光出射のトリガ信号Ｓ２等に基づいて、光強度の時間変化を示す計測波形を取得する信号処理手段である。計測波形の情報は電子データＤ１として取得され、演算処理部６０に入力される。

信号処理部５０は、図２に示すように、コンスタント・フラクション・ディスクリミネータ（ＣＦＤ）５１、時間−振幅変換器（ＴＡＣ）５２及びＡＤコンバータ（Ａ／Ｄ）５３で構成されている。光検出器（ＰＭＴ）３０からの出力信号である光検出信号は、ＣＦＤ５１を介してＴＡＣ５２に導かれて時間に対応したアナログ電圧に変換され、更にＡＤコンバータ５３でデジタル信号に変換される。このデジタル信号は、検出光強度の時間変化を示す計測波形のデータに対応するものである。

なお、本実施形態に係る散乱吸収体計測装置１０において、パルス光が入射されて当該パルス光の光強度の時間変化を示す計測波形が計測されるのは、後述するように、計測対象の散乱吸収体ＳＭ及びリファレンス用の散乱吸収体ＳＭである。リファレンス用の散乱吸収体ＳＭは、計測対象の散乱吸収体ＳＭの計測波形から装置関数の影響を除去するために計測されるものである。リファレンス用の散乱吸収体ＳＭとしては、計測対象の散乱吸収体ＳＭと屈折率がほぼ等しく、また、吸光係数や散乱係数等の理論波形を算出するのに必要な光学パラメータが既知であるものを用いる。

演算処理部６０は、信号処理部５０と電気的に接続されており、信号処理部５０から計測波形の情報の入力を受け付けて、当該計測波形の情報を用いて散乱吸収体ＳＭの内部情報を算出する演算処理を実行する。演算処理部６０は、上記の演算処理を行うために、図１に示すように、波形処理部６１と、計測波形格納部６２と、理論波形格納部６３と、内部情報算出部６４とを備えて構成されている。また、演算処理部６０は、それ以外にも、光源２０や光検出器３０等、上記の各構成要素を制御する機能を更に備えているとよい。

波形処理部６１は、信号処理部５０により取得された計測波形に基づき、装置関数の影響を除去した計測対象の散乱吸収体ＳＭの波形を取得するための波形処理手段である。即ち、波形処理部６１は、計測対象の散乱吸収体ＳＭの理論波形を取得する手段である。波形処理部６１は、計測対象の散乱吸収体ＳＭでの理論波形Ｐ_ｍ ^{ｔｈｅｏｒｙ}（ｔ）を、所定のパラメータを有する予め用意された関数で示されるものとする。波形処理部６１は、この関数を予め記憶している。ここで、上記のｔはパルス光が散乱吸収体ＳＭに入射された時刻（≒光源２０からパルス光が出射した時刻）を原点（ｔ＝０）として経過時間であり、関数値Ｐ_ｍ ^{ｔｈｅｏｒｙ}（ｔ）は各時刻における光強度を示す。

波形処理部６１は、装置関数の影響を除去した計測対象の散乱吸収体ＳＭの波形を取得するために、計測対象の散乱吸収体ＳＭでの理論波形とリファレンス用の散乱吸収体ＳＭでの計測波形とのコンボリューション演算の結果が、リファレンス用の散乱吸収体ＳＭでの理論波形と計測対象の散乱吸収体ＳＭでの計測波形とのコンボリューション演算の結果と等しくなるように演算処理を行う。即ち、波形処理部６１は、計測対象の散乱吸収体ＳＭでの計測波形をＰ_ｍ（ｔ）、リファレンス用の散乱吸収体ＳＭでの計測波形をＰ_ｒ（ｔ）、予め用意されたリファレンス用の散乱吸収体ＳＭでの理論波形をＰ_ｒ ^{ｔｈｅｏｒｙ}（ｔ）、コンボリューションを演算子“＊”とそれぞれしたときに、以下の式

が成り立つように上記のパラメータを特定することによって計測対象の散乱吸収体の理論波形Ｐ_ｍ ^{ｔｈｅｏｒｙ}（ｔ）を算出する。計測対象の散乱吸収体ＳＭでの理論波形を表すための関数、及び波形処理部６１による演算処理について、より具体的には後述する。波形処理部６１により処理された波形の情報は、内部情報算出部６４に出力される。

計測波形格納部６２は、信号処理部５０で所得された計測対象の散乱吸収体ＳＭの計測波形及びリファレンス用の散乱吸収体ＳＭの計測波形の情報（信号処理部５０から演算処理部６０に入力された上記の電子データＤ１）を格納するものである。計測波形格納部６２に格納された情報は、波形処理部６１による処理において適宜参照される。

理論波形格納部６３は、予め用意されたリファレンス用の散乱吸収体ＳＭの理論波形の情報を格納するものである。リファレンス用の散乱吸収体ＳＭの理論波形については、例えば、リファレンス用の散乱吸収体ＳＭの吸光係数や散乱係数等の光学パラメータから予め算出しておく。この算出は、例えば、有限要素法で行われる。それ以外でも、モンテカルロ法、差分法及びその他の解析解法等が用いられてもよい。理論波形格納部６３に格納された情報は、波形処理部６１による処理において適宜参照される。

内部情報算出部６４は、波形処理部６１により処理された計測対象の散乱吸収体ＳＭの計測波形に基づいて、計測対象の散乱吸収体ＳＭの内部情報を算出する内部情報算出手段である。内部情報の算出は、例えば、検出光の時間分解波形を利用する時間分解波形を利用する時間分解計測法（ＴＲＳ法：Time Resolved Spectroscopy）、あるいは、変調光を利用する位相変調計測法（ＰＭＳ法：Phase Modulation Spectroscopy）等による解析演算を適用することにより行われる。

上記の各構成要素を備える演算処理部６０は、例えば、図２に示すようなハードウェア構成により実現される。演算処理部６０においては、光源２０及び信号処理部５０にＣＰＵ（Central Processing Unit）７０が電気的に接続されている。これによって、光入射に同期した光検出のタイミング等がＣＰＵ７０によって制御されると共に、信号処理部５０から出力された計測波形の情報は、ＣＰＵ７０に導かれて所定の演算処理が行われる。また、入射パルス光の波長等、計測光の入射条件についても、このＣＰＵ７０によって制御あるいは選択される。

図２に示す演算処理部６０は、更に、オペレーティングシステム（ＯＳ）７１ａ及び所定の演算処理を行うための演算プログラム７１ｂが記憶されたプログラムメモリ７１と、各種データファイルが記憶されるデータファイルメモリ７２と、得られた散乱吸収体ＳＭの内部情報を示す情報を記憶するデータメモリ７３と、作業用データを一時的に記憶する作業用メモリ７４とを備えている。

プログラムメモリ７１の演算プログラム７１ｂには、上述した波形処理の演算や、内部情報を算出するための解析演算をそれぞれ実行するためのプログラム等が含まれている。また、データファイルメモリ７２は、計測波形の情報を格納する計測波形格納部６２、及び予め用意された理論波形の情報を格納する理論波形格納部６３が含まれている。また、必要な物理量等の諸データや、予め入力された計測条件や既知値等のデータ等も記憶される。

また、演算処理部６０は、データの入力を受け付けるキーボード７５ａ及びマウス７５ｂを備える入力装置７５と、得られたデータを出力するディスプレイ７６ａ及びプリンタ７６ｂを備える出力装置７６とを備えている。これらの演算処理部６０の各部の動作は、何れもＣＰＵ７０によって制御される。なお、上記の各メモリについては、コンピュータ内部に設置されているハードディスク等であってもよい。また、散乱吸収体計測装置１０の具体的なハードウェア構成については図２に示すものに限られるものではなく、必要に応じて変形又は拡張を行ってもよい。

ここで、計測波形と、散乱吸収体に対する理論波形と、装置関数との関係、及び波形処理部６１における演算処理について説明する。図５は、計測波形Ｐ（ｔ）、理論波形Ｐ^{ｔｈｅｏｒｙ}（ｔ）、及び装置関数ＩＦ（ｔ）の関係の一例を模式的に示すグラフである。このグラフにおいて、横軸は、パルス光が散乱吸収体ＳＭに入射された時刻を原点（ｔ＝０）とした経過時間ｔを示し、縦軸は、各時刻における光強度を示している。この光強度の経過時間ｔに依存した時間変化が、図示したそれぞれの時間波形となっている。

計測対象の散乱吸収体ＳＭに対してパルス光を照射し、散乱吸収体ＳＭ内部を伝播して出射された拡散光は、理想的には理論波形Ｐ^{ｔｈｅｏｒｙ}（ｔ）で表される。散乱吸収体の内部情報は、この理論波形Ｐ^{ｔｈｅｏｒｙ}（ｔ）に基づいて算出されるべきである。しかしながら、計測装置を用いた散乱吸収体計測では、計測対象である散乱吸収体による時間応答とは別に、計測装置の各部、各回路に起因する装置関数ＩＦ（ｔ）として、時間遅れや時間分散等の装置自体の時間応答を生じる。例えば、図１に示した散乱吸収体計測装置１０では、光源２０、波長選択器２１、光検出器３０、光ガイド４０ａ，４０ｂ、信号処理部５０、及びそれらの間の配線や回路等において、それぞれ時間応答を生じる。これらの時間応答を全て合わせたものを、散乱吸収体計測装置１０自体の時間応答を表す装置関数とする。従って、信号処理部５０で取得されて演算処理部６０へと入力される計測波形Ｐ（ｔ）は、散乱吸収体での理想的な計測波形（理論波形）Ｐ^{ｔｈｅｏｒｙ}（ｔ）に対して、上記の装置関数ＩＦ（ｔ）が重畳（コンボリューション）されたのものとなる。

これに対して、本実施形態の散乱吸収体計測装置１０では、（波形処理部６１が）計測対象の散乱吸収体ＳＭの理論波形とリファレンス用の散乱吸収体ＳＭでの計測波形とのコンボリューション演算の結果が、リファレンス用の散乱吸収体ＳＭでの理論波形と計測対象の散乱吸収体ＳＭでの計測波形とのコンボリューション演算の結果と等しくなるように演算処理を行うことにより、計測対象の散乱吸収体ＳＭの計測波形から装置関数の影響を除去する。具体的には、以下のように行う。

上述したように、計測対象の散乱吸収体ＳＭでの計測波形をＰ_ｍ（ｔ）、リファレンス用の散乱吸収体ＳＭでの計測波形をＰ_ｒ（ｔ）とし、それぞれの理論波形をＰ_ｍ ^{ｔｈｅｏｒｙ}（ｔ）、Ｐ_ｒ ^{ｔｈｅｏｒｙ}（ｔ）とする。また、装置関数をＩＦ（ｔ）とする。ここで、リファレンス用の散乱吸収体ＳＭでの理論波形Ｐ_ｒ ^{ｔｈｅｏｒｙ}（ｔ）は、散乱吸収体ＳＭの光学係数が既知であるため、予め得られている。計測波形は、理論波形に装置関数がコンボリューションされたものであるため、各々の波形の間には、次の式（１），式（２）に示す関係が成り立つ。

これらの式をフーリエ変換すると、以下の式（３）と式（４）とが得られる。但し、フーリエ変換をＦ［・］と表記する。

式（３），式（４）から次の式（５）に示す関係が得られる。

式（５）を変形すると以下の式（６）になる。

式（６）の逆フーリエ変換を行うと、次の式（７）が得られる。

式（７）は、計測対象の散乱吸収体ＳＭの理論波形とリファレンス用の散乱吸収体ＳＭでの計測波形とのコンボリューション演算の結果が、リファレンス用の散乱吸収体ＳＭでの理論波形と計測対象の散乱吸収体ＳＭでの計測波形とのコンボリューション演算の結果と等しくなることを表している。従って、式（７）を満足するようにして求められたＰ_ｍ ^{ｔｈｅｏｒｙ}（ｔ）が、計測対象の散乱吸収体ＳＭでの計測波形から装置関数の影響が除去されたものになる。

ここで、Ｐ_ｍ ^{ｔｈｅｏｒｙ}（ｔ）は、変数をｔとし所定のパラメータを有する予め用意された関数として扱われる。波形処理部６１による、式（７）を満足させるような関数の決定は、上記のパラメータを決定することで行われる。なお、Ｐ_ｍ ^{ｔｈｅｏｒｙ}（ｔ）を示す関数に含まれる決定対象となるパラメータは必ずしも１つでなくてもよく、後述するように複数のパラメータが含まれていてもよい。具体的には、例えば、最小二乗法を適用して、以下に示す式（８）を最小にするようなパラメータを決定する。ここでは、ｔを離散的な値ｔ_ｎ（ｎ＝０，１，…）（ｔ_０＝０）とし、上記の決定対象のパラメータ（の集合）をｐとする。

式（８）において、

である。また、式（８）におけるＮは、最小二乗法の解析で使用する時間範囲の最大値を表しており、予め設定された値である。例えば、Ｎ＝１０２４であれば、ｔ_０≦ｔ_ｎ≦ｔ_１０２３の範囲の値を用いて解析を行う。この時間範囲には、上記の各理論波形及び計測波形が含まれるようにされる。なお、ｔを離散的な値としているのは、通常、信号処理部５０で得られる波形が、デジタル値である時刻に対応する光強度の値であるためであり、連続的な値として考慮されてもよい。

Ｐ_ｍ ^{ｔｈｅｏｒｙ}（ｔ_ｎ）を示す関数としては、上記のパラメータを与えれば波形が一意に決まり、散乱吸収体ＳＭ内部の光伝播を近似できる（表すことができる）関数を用いる。具体的には例えば、パラメータをＡ，α，βとして、
（即ち、列ベクトルの形で表したパラメータｐは、

である（Ｔは転置行列を意味する））、Ｐ_ｍ ^{ｔｈｅｏｒｙ}（ｔ_ｎ）を示す関数ｆ（ｐ；ｔ_ｎ）として以下の式（９）を用いることができる。

ここで、ｃは光速を示す値、γは光入射位置Ａと光検出位置Ｂとの間の距離を示す値であり、設定値として予め波形処理部６１に記憶されている情報である。また、

である。散乱吸収体ＳＭ内部を伝播した光の時間的な強度変化（波形）は、およそ式（９）のような関数で表すことができる。従って、Ｐ_ｍ ^{ｔｈｅｏｒｙ}（ｔ_ｎ）を近似するには式（９）のような形の関数が適している。

式（９）は、半無限媒質からの反射光として光拡散方程式から解析的に得られる式、

において、

としたものである。ここでＲ（ｐ；ｔ_ｎ）は時刻ｔ_ｎのときの光強度（ｐは、上記の光拡散方程式から解析的に得られる式におけるパラメータの集合を示す）、μ_ａは散乱吸収体の吸収係数、Ｄは散乱吸収体の拡散係数（＝１／３μ´_ｓ）、μ´_ｓは散乱吸収体の散乱係数、ρは光入射位置と光検出位置との間の距離、ｃは光速をそれぞれ示す値である。なお、本実施形態においては、半無限媒質でない形状の散乱吸収体ＳＭに対して得られた波形を解析の対象とする場合にも式（９）を使用するため、式（９）における各パラメータに物理的な意味はない。

式（９）のように、所定のパラメータで波形が決まる関数を用いてＰ_ｍ ^{ｔｈｅｏｒｙ}（ｔ）を近似することにより、計測波形に含まれるノイズに影響されずに理論波形を求めることができる。即ち、Ｐ_ｍ ^{ｔｈｅｏｒｙ}（ｔ_ｎ）を式（９）のような関数で近似せずに、直接求めた場合（各時刻におけるＰ_ｍ ^{ｔｈｅｏｒｙ}（ｔ_ｎ）の値そのものを求めた場合）には、各時刻における計測データのばらつき（ノイズ）が影響を与える。一方、Ｐ_ｍ ^{ｔｈｅｏｒｙ}（ｔ_ｎ）を式（９）で近似した場合、パラメータＡ，α，βを与えれば、各時刻におけるＰ_ｍ ^{ｔｈｅｏｒｙ}（ｔ_ｎ）の値が一意に求まる。従って、各時刻における計測データのばらつきには影響されない。

上記のようにＰ_ｍ ^{ｔｈｅｏｒｙ}（ｔ_ｎ）を近似する関数は式（９）を用いることが望ましいが、式（９）に限定されるものではない。所定のいくつかのパラメータを与えれば波形が一意に決まり、散乱吸収体ＳＭ内部の光伝播を近似できる関数であればどのようなものでも構わない。また、モンテカルロ法、有限差分法及び有限要素法等の数値計算によって得られた波形を示す関数を用いることも可能である。

引き続いて、図４のフローチャートを用いて、散乱吸収体計測装置１０を用いて実行される本実施形態に係る散乱吸収体計測方法を説明する。散乱吸収体計測方法は、計測対象の散乱吸収体ＳＭの内部情報を計測する方法である。

まず、本実施形態に係る散乱吸収体計測方法では、リファレンス用の散乱吸収体ＳＭでの計測波形Ｐ_ｒ（ｔ）の取得が行われる（Ｓ０１）。計測波形の取得は、具体的には、以下のように行われる。まず、散乱吸収体計測装置１０の所定位置にリファレンス用の散乱吸収体ＳＭを配置する（Ｓ０１ａ）。この配置は、散乱吸収体計測装置１０により散乱吸収体ＳＭへのパルス光の入射及び散乱吸収体ＳＭからの光検出ができるように行われる。例えば、散乱吸収体計測装置１０が備えるガントリに光入射位置及び光検出位置が設けられる場合は、ガントリにリファレンス用の散乱吸収体ＳＭを満たす。なお、上述したようにリファレンス用の散乱吸収体ＳＭでの理論波形Ｐ_ｒ ^{ｔｈｅｏｒｙ}（ｔ）は、予め算出されており理論波形格納部６３に格納されている。

続いて、散乱吸収体計測装置１０では、光源２０により発生し波長選択器２１により波長選択されたパルス光が、光入射用の光ガイド４０ａを通って所定の光入射位置からリファレンス用の散乱吸収体ＳＭに入射される（Ｓ０１ｂ、光入射ステップ）。入射されたパルス光は、散乱吸収体内部を散乱し、吸収されつつ伝播する。散乱吸収体計測装置１０では、所定の光検出位置において、光検出用の光ガイド４０ｂを通って光検出器３０により検出される（Ｓ０１ｃ、光検出ステップ）。続いて、信号処理部５０により、検出された光検出信号に基づいて、リファレンス用の散乱吸収体ＳＭでの光強度の時間変化を示す計測波形Ｐ_ｒ（ｔ）が取得される（Ｓ０１ｄ、信号処理ステップ）。取得された計測波形Ｐ_ｒ（ｔ）の情報は、演算処理部６０に出力され、計測波形格納部６２に格納される。

ここで、光入射位置と光検出位置との組み合わせにより、散乱吸収体ＳＭの内部情報を十分に取得できるように、パルス光を入射させる位置を何個所か移動させて上記の処理（Ｓ０１ｂ〜Ｓ０１ｄ）を繰り返し行い、複数の計測波形を取得するようにしてもよい。また、１個所からのパルス光の入射に対して、複数の光検出位置での光検出を行ってもよい。例えば、ガントリにおける１０箇所の光入射位置でのパルス光の入射を行うようにし、１つの光入射位置に対して、３０箇所の光検出位置での光検出を行うこととしてもよい。その場合、１０×３０の計測波形が取得される。

続いて、計測対象の散乱吸収体ＳＭでの計測波形計測波形をＰ_ｍ（ｔ）の取得が行われる（Ｓ０２）。計測対象の散乱吸収体ＳＭでの計測波形計測波形をＰ_ｍ（ｔ）の取得は、以下のように行われる。まず、リファレンス用の散乱吸収体ＳＭでの計測と同様に散乱吸収体計測装置１０の所定位置に、計測対象の散乱吸収体ＳＭを配置する（Ｓ０２ａ）。例えば、散乱吸収体計測装置１０が備えるガントリに計測対象の散乱吸収体ＳＭである乳房を配置する。続いて、リファレンス用の散乱吸収体ＳＭでの計測と同じ条件で、上記と同様に光入射、光検出及び計測波形Ｐ_ｍ（ｔ）の取得が行われる（Ｓ０２ｂ〜Ｓ０２ｄ、光入射ステップ、光検出ステップ、信号処理ステップ）。取得された計測波形Ｐ_ｍ（ｔ）の情報は、演算処理部６０に出力され、計測波形格納部６２に格納される。なお、リファレンス用の散乱吸収体ＳＭでの計測波形Ｐ_ｒ（ｔ）の取得（Ｓ０１）と、計測対象の散乱吸収体ＳＭでの計測波形Ｐ_ｍ（ｔ）の取得（Ｓ０２）とは順序が逆になってもよい。

続いて、演算処理部６０の波形処理部６１により計測波形格納部６２から各計測波形Ｐ_ｍ（ｔ）及びＰ_ｒ（ｔ）とリファレンス用の散乱吸収体ＳＭでの理論波形Ｐ_ｒ ^{ｔｈｅｏｒｙ}（ｔ）を示す情報が取得される。また、波形処理部６１により、決定対象のパラメータを有する（自身により記憶されている）理論波形Ｐ_ｍ ^{ｔｈｅｏｒｙ}（ｔ）を示す関数が読み出される。続いて、上述したように、波形処理部６１により、計測対象の散乱吸収体ＳＭでの理論波形Ｐ_ｍ ^{ｔｈｅｏｒｙ}（ｔ）とリファレンス用の散乱吸収体ＳＭでの計測波形Ｐ_ｒ（ｔ）とのコンボリューション演算の結果が、リファレンス用の散乱吸収体ＳＭでの理論波形Ｐ_ｒ ^{ｔｈｅｏｒｙ}（ｔ）と計測対象の散乱吸収体ＳＭでの計測波形Ｐ_ｍ（ｔ）とのコンボリューション演算の結果と等しくなるように演算処理が行われて、計測対象の散乱吸収体ＳＭでの理論波形Ｐ_ｍ ^{ｔｈｅｏｒｙ}（ｔ）のパラメータが算出され、特定される。このようにパラメータが特定された理論波形Ｐ_ｍ ^{ｔｈｅｏｒｙ}（ｔ）は、計測波形Ｐ_ｍ（ｔ）から装置関数ＩＦ（ｔ）の影響が除去されたものである（Ｓ０３、波形処理ステップ）。当該処理により得られた（装置関数ＩＦ（ｔ）の影響が除去された）計測対象の散乱吸収体ＳＭでの時間応答波形Ｐ_ｍ ^{ｔｈｅｏｒｙ}（ｔ）は、内部情報算出部６４に出力される。

続いて、内部情報算出部６４により、上記のように特定されたパラメータを含む関数により示される、計測対象の散乱吸収体ＳＭでの理論波形Ｐ_ｍ ^{ｔｈｅｏｒｙ}（ｔ）に基づいて、計測対象の散乱吸収体ＳＭの内部情報が算出される（Ｓ０４、内部情報算出ステップ）。算出される内部情報は、散乱吸収体ＳＭの吸光係数や散乱係数等である。算出された内部情報は、例えば、ユーザにより参照可能なように、分布化されて演算処理部６０のディスプレイ７６ａに出力される。

上述したように本実施形態に係る散乱吸収体計測方法及び散乱吸収体計測装置１０では、上記の式（７）のような波形の関係に基づいて、計測対象の散乱吸収体ＳＭでの理論波形Ｐ_ｍ ^{ｔｈｅｏｒｙ}（ｔ）を示す関数のパラメータを特定することによって、装置関数を算出することなく、計測対象の散乱吸収体の計測波形から装置関数の影響を除去した理論波形を得ることができる。即ち、装置関数ＩＦ（ｔ）を算出することなく、装置関数ＩＦ（ｔ）の影響が除去された波形に基づいて、計測対象の散乱吸収体ＳＭの内部情報を算出することができる。

また、計測対象の散乱吸収体の理論波形Ｐ_ｍ ^{ｔｈｅｏｒｙ}（ｔ）を上記のような所定のパラメータを有する予め用意された関数とすることによって、計測波形のばらつきの影響を防止することができる。即ち、計測波形Ｐ_ｍ（ｔ）にノイズが含まれる場合であっても、Ｐ_ｍ（ｔ）に対応する計測対象の散乱吸収体ＳＭの理論波形Ｐ_ｍ ^{ｔｈｅｏｒｙ}（ｔ）をより正確に求めることができる。従って、本発明によれば、計測波形Ｐ_ｍ（ｔ）にノイズが含まれる場合であっても、より正確かつ簡易に装置関数ＩＦ（ｔ）を考慮して散乱吸収体の内部情報を取得することができる。

また、本実施形態のように計測対象の散乱吸収体ＳＭでの理論波形Ｐ_ｍ ^{ｔｈｅｏｒｙ}（ｔ）を示す関数を式（９）のような、散乱吸収体内部の光伝播をより適切に示すものとすることが望ましい。このような構成にすれば、更に正確に計測対象の散乱吸収体の理論波形を求めることができ、結果として、更に正確に散乱吸収体の内部情報を取得することができる。

本実施形態に係る方法によって、散乱吸収体の理論波形Ｐ_ｍ ^{ｔｈｅｏｒｙ}（ｔ）を算出した例を図５〜図７を用いて説明する。図５に、計測対象の散乱吸収体ＳＭでの計測波形Ｐ_ｍ（ｔ）、リファレンス用の散乱吸収体ＳＭでの計測波形Ｐ_ｒ（ｔ）、及びリファレンス用の散乱吸収体ＳＭでの理論波形Ｐ_ｒ ^{ｔｈｅｏｒｙ}（ｔ）のグラフを示す。図５のグラフにおいて、横軸はチャネル（経過時間、時刻に相当）であり、縦軸は光強度（ａ．ｕ．）である（図６及び図７でも同様）。

図６に、図５に示すデータを用いて、特許文献１に記載の方法を用いて算出した計測対象の散乱吸収体ＳＭでの理論波形、及び実際の理論波形（真値）を示す。特許文献１に記載された方法では、フーリエ変換した波形の演算が計測波形に含まれるノイズに敏感であるため、計測波形の取り扱いが煩雑であった。また、そのままの計測波形を使用して演算を行った場合、ノイズの影響により、図６に示すように処理後の波形（算出された理論波形）が真値から大きくずれて歪んでしまう場合があった。

図７に、図５に示すデータを用いて、本実施形態に係る方法を用いて算出した計測対象の散乱吸収体ＳＭでの理論波形Ｐ_ｍ ^{ｔｈｅｏｒｙ}（ｔ）、及び実際の理論波形（真値）を示す。本実施形態に係る方法では、処理後の波形（算出された理論波形）が真値とほぼ一致しており、この結果からも本実施形態による方法によれば、ノイズが含まれる場合であっても正確な理論波形を導出することができ、その結果、正確な散乱吸収体の内部情報を取得することができることがわかる。

本発明の実施形態に係る散乱吸収体計測装置の機能構成を示す図である。 本発明の実施形態に係る散乱吸収体計測装置のハードウェア構成を示す図である。 計測波形、理論波形、及び装置関数の一例を模式的に示すグラフである。 本発明の実施形態に係る散乱吸収体計測方法を示すフローチャートである。 計測対象の散乱吸収体での計測波形、リファレンス用の散乱吸収体での計測波形、及びリファレンス用の散乱吸収体での理論波形の例を示すグラフである。 計測対象の散乱吸収体での理論波形（真値）、及び従来の方法で校正された計測対象の散乱吸収体での計測波形（算出された理論波形）の例を示すグラフである。 計測対象の散乱吸収体での理論波形（真値）、及び本発明に係る方法で校正された計測対象の散乱吸収体での計測波形（算出された理論波形）の例を示すグラフである。

符号の説明

１０…散乱吸収体計測装置、２０…光源、２１…波長選択器、３０…光検出器、４０ａ，４０ｂ…光ガイド、５０…信号処理部、６０…演算処理部、６１…波形処理部、６２…計測波形格納部、６３…理論波形格納部、６４…内部情報算出部。

Claims (3)

1. 計測対象の散乱吸収体及びリファレンス用の散乱吸収体に対して、所定波長のパルス光を光入射位置から入射する光入射ステップと、
   前記光入射ステップにおいて入射されて前記各散乱吸収体の内部を伝播したパルス光を光検出位置で検出して光検出信号を取得する光検出ステップと、
   前記光検出ステップにおいて検出された光検出信号に基づいて、光強度の時間変化を示す計測波形を取得する信号処理ステップと、
   前記信号処理ステップにおいて取得された前記計測対象の散乱吸収体での計測波形をＰ_ｍ（ｔ）、所定のパラメータを有する予め用意された関数によって示される前記計測対象の散乱吸収体での理論波形をＰ_ｍ ^{ｔｈｅｏｒｙ}（ｔ）、前記信号処理ステップにおいて取得された前記リファレンス用の散乱吸収体での計測波形をＰ_ｒ（ｔ）、予め用意された前記リファレンス用の散乱吸収体での理論波形をＰ_ｒ ^{ｔｈｅｏｒｙ}（ｔ）、コンボリューションを演算子“＊”とそれぞれしたときに、以下の式
   

   

   
   に基づいて前記パラメータを特定する波形処理ステップと、
   前記波形処理ステップにおいて特定されたパラメータを有する前記関数によって示される前記計測対象の散乱吸収体での理論波形に基づいて、前記計測対象の散乱吸収体の内部情報を算出する内部情報算出ステップと、
   を含む散乱吸収体計測方法。
2. Ｐ_ｍ ^{ｔｈｅｏｒｙ}（ｔ）を示す前記関数は、離散的な時刻ｔ_ｎの関数であり、Ａ、α、βを前記パラメータ、ｃを光速を示す値、γを前記光入射位置と前記光検出位置との間の距離を示す値とそれぞれしたときに、以下の式
   

   

   
   （ここで、
   

   

   
   である）で表されることを特徴とする請求項１に記載の散乱吸収体計測方法。
3. 計測対象の散乱吸収体及びリファレンス用の散乱吸収体に対して、所定波長のパルス光を光入射位置から入射する光入射手段と、
   前記光入射手段によって入射されて前記各散乱吸収体の内部を伝播したパルス光を光検出位置で検出して光検出信号を取得する光検出手段と、
   前記光検出手段によって検出された光検出信号に基づいて、光強度の時間変化を示す計測波形を取得する信号処理手段と、
   前記信号処理手段によって取得された前記計測対象の散乱吸収体での計測波形をＰ_ｍ（ｔ）、所定のパラメータを有する予め用意された関数によって示される前記計測対象の散乱吸収体での理論波形をＰ_ｍ ^{ｔｈｅｏｒｙ}（ｔ）、前記信号処理手段によって取得された前記リファレンス用の散乱吸収体での計測波形をＰ_ｒ（ｔ）、予め用意された前記リファレンス用の散乱吸収体での理論波形をＰ_ｒ ^{ｔｈｅｏｒｙ}（ｔ）、コンボリューションを演算子“＊”とそれぞれしたときに、以下の式
   

   

   
   に基づいて前記パラメータを特定する波形処理手段と、
   前記波形処理手段によって特定されたパラメータを有する前記関数によって示される前記計測対象の散乱吸収体での理論波形に基づいて、前記計測対象の散乱吸収体の内部情報を算出する内部情報算出手段と、
   を備える散乱吸収体計測装置。
   

Images