(57)【要約】【課題】サンプル内部の吸収分布または散乱分布を正確に求めることができる光学的測定方法を提供する。【解決手段】本発明の光学的測定方法は、(1) サンプルからの反射光を検出し、複数の時間帯または波長帯におけるサンプル内部の反射位置と反射光量との関係を取得するデータ取得ステップと、(2) この関係に基づいて、入射光量に対する反射光量の減衰量に含まれる吸収成分および散乱成分を互いに分離して求めるデータ処理ステップと、を備える。データ処理ステップでは、反射光量の信号品質指標の値に応じて平滑化パラメータの値を設定して、反射光量，吸収成分，散乱成分，ならびに，平滑化パラメータの設定値が乗ぜられた吸収成分の空間二次微分および散乱成分の空間二次微分を含む目的関数が最小となるように、吸収成分および散乱成分を求める。【選択図】図１
2014102087.tif 000002
【特許請求の範囲】
【請求項１】
サンプルに光を照射または伝搬させたときに該サンプルにおいて生じる反射光を時間または波長の関数として検出して、複数の時間帯または波長帯における前記サンプル内部の反射位置と反射光量との関係を取得するデータ取得ステップと、
このデータ取得ステップにおいて取得された複数の時間帯または波長帯における前記サンプル内部の反射位置と反射光量との関係に基づいて、入射光量に対する前記反射光量の減衰量に含まれる吸収成分および散乱成分を互いに分離して求め、前記サンプル内部の吸収分布または散乱分布を求めるデータ処理ステップと、
を備え、
前記データ処理ステップにおいて、
前記反射光量の信号品質指標の値に応じて平滑化パラメータの値を設定して、
前記反射光量，前記吸収成分，前記散乱成分，ならびに，前記平滑化パラメータの設定値が乗ぜられた前記吸収成分の空間二次微分および前記散乱成分の空間二次微分を含む目的関数が最小となるように、前記吸収成分および前記散乱成分を求める、
ことを特徴とする光学的測定方法。
【請求項２】
前記データ処理ステップにおいて、
前記信号品質指標の基準値に対応する前記平滑化パラメータの基準値を予め記憶装置に格納しておき、
前記信号品質指標の値が前記信号品質指標の基準値と等しい場合に前記平滑化パラメータの設定値を前記平滑化パラメータの基準値と等しくし、
前記信号品質指標の値が前記信号品質指標の基準値より低い場合に前記平滑化パラメータの設定値を前記平滑化パラメータの基準値より小さくし、
前記信号品質指標の値が前記信号品質指標の基準値より高い場合に前記平滑化パラメータの設定値を前記平滑化パラメータの基準値より大きくする、
ことを特徴とする請求項１に記載の光学的測定方法。
【請求項３】
前記データ取得ステップにおいて、前記サンプルにおいて生じる反射光と参照光とを干渉させた光を時間または波長の関数として検出し、この関数をフーリエ変換することで前記サンプル内部の反射位置と反射光量との関係を取得する、ことを特徴とする請求項１または２に記載の光学的測定方法。
【請求項４】
前記データ取得ステップにおいて、前記サンプルにおいて生じる反射光と参照光とを干渉させた光を波長の関数として検出し、この関数を複数の波長帯に分割して波長帯毎にフーリエ変換することで波長帯毎の前記サンプル内部の反射位置と反射光量との関係を取得する、ことを特徴とする請求項１または２に記載の光学的測定方法。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、光学的測定方法に関するものである。
【背景技術】
【０００２】
ＯＣＴ（Optical Coherence Tomography），ＯＴＤＲ（Optical Time Domain Reflectometer）およびＤＯＴ（DiffuseOptical Tomography）等の光学的測定方法では、サンプルに光を照射または伝搬させたときに該サンプルにおいて生じる反射光を検出して、その検出結果に基づいてサンプルの情報を得ることができる。ＯＣＴでは、サンプルにおいて生じる反射光と参照光とを干渉させた光を検出することでサンプルの深さ方向の一次元の光断層画像を取得することができ、また、サンプルへの光照射の位置を走査することで二次元または三次元の光断層画像を取得することができる。
【０００３】
このような光学的測定方法において、サンプルに入射した光量に対してサンプルで生じた反射光の光量は吸収および散乱の双方による減衰が含まれる場合がある。非特許文献１に記載された光学的測定方法は、入射光量に対する反射光量の減衰量に含まれる吸収成分および散乱成分を互いに分離して求め、サンプル内部の吸収分布または散乱分布を求めることができる。
【先行技術文献】
【非特許文献】
【０００４】
【非特許文献１】Chenyang Xu, et al., "Separationof absorption and scattering profiles in spectroscopic optical coherencetomography using a least-squares algorithm," Optics Express, Vol.12, No.20,pp.4790-4803 (2004).
【発明の概要】
【発明が解決しようとする課題】
【０００５】
非特許文献１に記載された光学的測定方法により得られる結果は、計算時のパラメータの設定値によってはサンプル内部の吸収分布または散乱分布を正確に求めることができないという問題点を有していることを、本発明者は見出した。本発明は、上記問題点を解消する為になされたものであり、サンプル内部の吸収分布または散乱分布を正確に求めることができる光学的測定方法を提供することを目的とする。
【課題を解決するための手段】
【０００６】
本発明の光学的測定方法は、(1) サンプルに光を照射または伝搬させたときに該サンプルにおいて生じる反射光を時間または波長の関数として検出して、複数の時間帯または波長帯におけるサンプル内部の反射位置と反射光量との関係を取得するデータ取得ステップと、(2) このデータ取得ステップにおいて取得された複数の時間帯または波長帯におけるサンプル内部の反射位置と反射光量との関係に基づいて、入射光量に対する反射光量の減衰量に含まれる吸収成分および散乱成分を互いに分離して求め、サンプル内部の吸収分布または散乱分布を求めるデータ処理ステップと、を備えることを特徴とする。さらに、本発明の光学的測定方法は、データ処理ステップにおいて、反射光量の信号品質指標の値に応じて平滑化パラメータの値を設定して、反射光量，吸収成分，散乱成分，ならびに，平滑化パラメータの設定値が乗ぜられた吸収成分の空間二次微分および散乱成分の空間二次微分を含む目的関数が最小となるように、吸収成分および散乱成分を求めることを特徴とする。
【０００７】
データ処理ステップにおいて、信号品質指標の基準値に対応する平滑化パラメータの基準値を予め記憶装置に格納しておき、反射光量の信号品質指標の値に応じて平滑化パラメータを増減させるのが好適である。反射光量の信号品質指標の値が信号品質指標の基準値と等しい場合に平滑化パラメータの設定値を平滑化パラメータの基準値と等しくし、信号品質指標の値が信号品質指標の基準値より低い場合に平滑化パラメータの設定値を平滑化パラメータの基準値より小さくし、信号品質指標の値が信号品質指標の基準値より高い場合に平滑化パラメータの設定値を平滑化パラメータの基準値より大きくするのが好適である。また、複数の信号品質指標の基準値と平滑化パラメータの基準値の組み合わせをデータベースとして格納しておき、反射光量の信号品質指標の値に対応する平滑化パラメータをデータベースを参照して算出するのが好適である。その算出方法としては、線形補間や多項式補間を用いる。
【０００８】
データ取得ステップにおいて、サンプルにおいて生じる反射光と参照光とを干渉させた光を時間または波長の関数として検出し、この関数をフーリエ変換することでサンプル内部の反射位置と反射光量との関係を取得するのが好適である。また、データ取得ステップにおいて、サンプルにおいて生じる反射光と参照光とを干渉させた光を波長の関数として検出し、この関数を複数の波長帯に分割して波長帯毎にフーリエ変換することで波長帯毎のサンプル内部の反射位置と反射光量との関係を取得するのが好適である。
【発明の効果】
【０００９】
本発明によれは、サンプル内部の吸収分布または散乱分布を正確に求めることができる光学的測定方法を提供することができる。
【図面の簡単な説明】
【００１０】
【図１】本実施形態の光学的測定方法が適用される光学的測定システム１の概略構成を示す図である。
【図２】光学的測定システム１により取得されるデータ等の例を示す図である。
【図３】吸光係数εａ(λ)および散乱係数εｓ(λ)それぞれの波長依存性の例を示すグラフである。
【図４】行列Ｙの一例を模式的に示す図である。
【図５】行列Ａの一例を模式的に示す図である。
【図６】行列Ｌの一例を模式的に示す図である。
【図７】行列Ｘの一例を模式的に示す図である。
【図８】本実施形態におけるＦａの二次元画像の一例を示す図である。
【図９】作成した干渉信号の例を示す図である。
【図１０】病変を模擬した、所定の光吸収スペクトルを有する吸収性散乱体と、正常組織を模擬した、光吸収を有さない非吸収性散乱体を想定したサンプルそれぞれの吸収スコアの例を示す図である。
【図１１】ＲＯＣ曲線を示す図である。
【図１２】フィッティング品質と平滑化パラメータとの関係を示すグラフである。
【発明を実施するための形態】
【００１１】
以下、添付図面を参照して、本発明を実施するための形態を詳細に説明する。なお、図面の説明において同一の要素には同一の符号を付し、重複する説明を省略する。
【００１２】
図１は、本実施形態の光学的測定方法が適用される光学的測定システム１の概略構成を示す図である。光学的測定システム１は、ＯＣＴ装置２およびプローブ３を備え、サンプル９の光断層画像を取得する。ＯＣＴ装置２は、光源から出力された光を２分岐して観察光および参照光として出力する。プローブ３は、ＯＣＴ装置２の光分岐部から出力された観察光を導光して先端からサンプル９へ照射するとともに、その光照射によって該サンプル９で生じる反射光を先端に入力してＯＣＴ装置２へ導光する。ＯＣＴ装置２は、プローブ３により導光されて到達した反射光と参照光とを干渉させて当該干渉光を検出する。そして、ＯＣＴ装置２は、この干渉光に基づいて、サンプルの深さ方向の光断層画像を取得することができる。また、プローブ３によるサンプル９への光照射の位置を走査（スキャン）することで、ＯＣＴ装置２は二次元または三次元の光断層画像を取得することができる。
【００１３】
本実施形態の光学的測定方法は、データ取得ステップおよびデータ処理ステップを備える。データ取得ステップでは、サンプルに光を照射または伝搬させたときに該サンプルにおいて生じる反射光を時間または波長の関数として検出して、複数の時間帯または波長帯におけるサンプル内部の反射位置と反射光量との関係を取得する。データ処理ステップでは、データ取得ステップにおいて取得された複数の時間帯または波長帯におけるサンプル内部の反射位置と反射光量との関係に基づいて、入射光量に対する反射光量の減衰量に含まれる吸収成分および散乱成分を互いに分離して求め、サンプル内部の吸収分布または散乱分布を求める。
【００１４】
図２は、光学的測定システム１により取得されるデータ等の例を示す図である。データ取得ステップにおいて、干渉光を検出することで、波長λに対する干渉光強度（信号強度）を得ることができる。この例では、波長帯を６分割して波長帯毎にフーリエ変換することで、波長帯毎の光断層画像を取得することができる。続く、データ処理ステップでは、この波長帯毎の光断層画像に基づいて以下のような処理を行う。
【００１５】
各波長帯の光断層画像は、該波長帯における散乱および吸収の双方の情報を含んでいる。ここで、サンプル９における波長λ、深さ方向の位置zに対するログスケールでの反射光量をＹ(λ,ｚ)とし、サンプル９における単位濃度当りの吸光係数をεａ(λ)とし、サンプル９における単位濃度当りの散乱係数をεｓ(λ)とし、サンプル９における吸収性散乱体の濃度のｚ＝０からの積分値をＦａ(ｚ)とし、サンプル９における非吸収性散乱体の濃度のｚ＝０からの積分値をＦｓ(ｚ)とし、サンプル９における波長λに依存しない他の項のｚ＝０からの積分値をＣ(ｚ)とする。サンプル９の表面をｚ＝０とする。これらのパラメータの間の関係は下記（１）式で表される。
【００１６】
【数１】
2014102087.tif 000003
【００１７】
吸光係数εａ(λ)および散乱係数εｓ(λ)が既知であれば、未知数がＦａ(ｚ)、Ｆｓ(ｚ)およびＣ(ｚ)の３つであるので、３つ以上の波長帯で光断層像を得て連立方程式を解けば、Ｆａ(ｚ)およびＦｓ(ｚ)を得ることができる。そして、Ｆａ(ｚ)を微分することで、サンプル９における吸収性散乱体の空間分布を得ることができる。
【００１８】
吸光係数εａ(λ)は、ＦＴ-ＩＲによる実測値として得ることができる。散乱係数εｓ(λ)は、下記（２）式で表される波長λを変数とする一次関数のモデルを用いることができる。（２）式のa,bは散乱係数の一次関数モデルにおけるそれぞれ傾きと切片である。（２）式のλ０は、全波長帯の中心波長とする。図３は、吸光係数εａ(λ)および散乱係数εｓ(λ)それぞれの波長依存性の例を示すグラフである。（２）式を用いると、（１）式は下記（３）式のように表される。また、波長λを変数とする（３）式に替えて、波数ｋを変数とする下記（４）式を用いてもよい。（４）式で、波長λ０に対応する波数がｋ０である。
【００１９】
【数２】
2014102087.tif 000004
【００２０】
【数３】
2014102087.tif 000005
【００２１】
【数４】
2014102087.tif 000006
【００２２】
６分割後の各波長帯の波数をｋ１〜ｋ６とする。このとき、解くべき連立方程式は、下記（５）式および（６）式のように行列を用いて表現することができる。６行３列の行列Ａは、既知の行列要素からなる。３行１列の行列Ｘは、データ処理ステップにおいて求められるべき未知の行列要素からなる。６行１列の行列Ｙは、データ取得ステップにおいて得られた反射光量からなる。
【００２３】
【数５】
2014102087.tif 000007
【００２４】
【数６】
2014102087.tif 000008
【００２５】
（５）式の両辺に行列Ａの転置行列ＡＴを左から掛けると、下記（７）式となる。この（７）式から下記（８）式が得られる。このように行列表現を用いると、解Ｘは表される。
【００２６】
【数７】
2014102087.tif 000009
【００２７】
【数８】
2014102087.tif 000010
【００２８】
Ｆａ(ｚ)を求めるには、方程式の数（すなわち、分割した波長帯の数）は３であればよい。しかし、実際にはノイズの影響があるので、方程式の数が３であると、得られる解が不安定になり易い。そこで、方程式の数（すなわち、分割する波長帯の数）は、３より大きい６としている。ただし、波長帯の分割数が多すぎると光断層像の分解能が低下するので、波長帯の分割数を適切に設定することが好ましい。
【００２９】
未知変数の個数に対して方程式の個数が多いので、最小二乗法を用いて、下記（９）式で表される目的関数Ｓが最小（目的関数の偏微分がゼロ）となるような解Ｘを求める。しかし、実際に計算して解Ｘを求めてみると、光断層像の空間的なノイズが重畳された解Ｘが得られる。
【００３０】
【数９】
2014102087.tif 000011
【００３１】
このノイズ問題に対処する為、Ｆａ(ｚ)，Ｆｓ'(ｚ)それぞれがノイズに比べて空間的に滑らかに変化する（空間的な変化が小さい）ことを仮定して、目的関数を修正する。すなわち、下記（１０）式で表されるＦａ(ｚ)，Ｆｓ'(ｚ)それぞれの二次微分値が小さいと仮定する。修正後の目的関数Ｓ' は下記（１１）式で表される。αは、Ｓに対する二次微分値の重みを決めるパラメータである。α２は、Ｆａ(ｚ)，Ｆｓ'(ｚ)それぞれの二次微分値を小さくするためのパラメータであるので、平滑化パラメータと呼ぶことにする。目的関数Ｓ' は、反射光量Ｙ，吸収成分Ｆａ，散乱成分Ｆｓ'，吸収成分Ｆａの空間二次微分および散乱成分Ｆｓ'の空間二次微分を含む関数である。
【００３２】
【数１０】
2014102087.tif 000012
【００３３】
【数１１】
2014102087.tif 000013
【００３４】
行列を用いて表現すると下記（１２）式のようになる。Ｌは二次微分演算子（ラプラシアン演算子）である。この（１２）式から、解Ｘは下記（１３）式で表される。
【００３５】
【数１２】
2014102087.tif 000014
【００３６】
【数１３】
2014102087.tif 000015
【００３７】
Ｆａ(ｚ)，Ｆｓ'(ｚ)は離散的なデータであるので、或る位置での二次微分値は、該位置のデータに加えて前後のデータも用いた計算により求める。計算により求められる二次微分値は、計算手法によって異なり、また、計算に用いるデータの点数によっても異なる。計算手法としては、差分法、Savitzky-Golay filterを用いた数値微分、前進差分近似、後退差分近似、中心差分近似、Stirlingの公式、などが知られている。
【００３８】
図４は、行列Ｙの一例を模式的に示す図である。図５は、行列Ａの一例を模式的に示す図である。図６は、行列Ｌの一例を模式的に示す図である。図７は、行列Ｘの一例を模式的に示す図である。深さ方向のデータ数は５１２であり、波長帯の数は６である。差分法により３点のデータを用いて二次微分を計算する場合の行列Ｌを示している。二次微分演算子である行列Ｌは、Ｆａ(ｚ)，Ｆｓ'(ｚ)に作用するが、Ｄ(ｚ)には作用しない。
【００３９】
図８は、本実施形態におけるＦａの二次元画像の一例を示す図である。ここでは、７点のデータを用いたSavitzky-Golay filterを用いた数値微分による二次微分演算子と、９点のデータを用いたSavitzky-Golay filterを用いた数値微分による二次微分演算子と、の線形結合で表される二次微分演算子を採用した。
【００４０】
ところで、本発明者の知見によれば、上記（１１）式の目的関数Ｓ'における平滑化パラメータα２の値が適切でないと、フィッティング品質が悪くなり、ノイズ対策の効果が小さくなって、サンプル内部の吸収分布または散乱分布を正確に求めることができない。平滑化パラメータα２の値が大きすぎると、ノイズの影響を低減することができるものの、本来の信号が崩れてしまう。逆に、平滑化パラメータα２の値が小さすぎると、本来の信号の崩れを抑制することができるが、ノイズの影響の低減が不充分となる。したがって、フィッティング品質が良くなるように、平滑化パラメータα２の値を適切に設定する必要がある。
【００４１】
さらに本発明者の知見によれば、フィッティング品質が良くなる平滑化パラメータα２の好適値は、反射光量Ｙの信号品質指標との間に相関を有する。すなわち、反射光量Ｙの信号品質指標が高いほど、平滑化パラメータα２の設定値を大きくすれば、フィッティング品質が良くなって、サンプル内部の吸収分布または散乱分布を正確に求めることができる。
【００４２】
そこで、本実施形態の光学的測定方法では、反射光量Ｙの信号品質指標の値に応じて平滑化パラメータα２の値を設定して、上記（１１）式の目的関数Ｓ'が最小となるように吸収成分および散乱成分を求める。好適には、データ処理ステップにおいて、複数の信号品質指標の基準値と平滑化パラメータの基準値の組み合わせをデータベースとして記憶装置に格納しておき、反射光量Ｙの信号品質指標の値に対応する平滑化パラメータをデータベースを参照して算出するのが好適である。その算出方法としては、線形補間や多項式補間を用いる。
【００４３】
反射光量Ｙの信号品質指標の基準値に対応する平滑化パラメータの基準値は、以下のようにして求めることができる。すなわち、反射光量Ｙの信号品質指標の基準値に対して平滑化パラメータを各値に設定し、平滑化パラメータを各値についてフィッティング品質を求める。そして、フィッティング品質が最も良くなるときの平滑化パラメータを基準値とする。これら信号品質指標の基準値および平滑化パラメータの基準値を予め求めておいて記憶装置に格納しておく。
【００４４】
反射光量Ｙの信号品質指標は、下記（１４）式で定義されるＳＮＲ（Signal-to-NoiseRate）であってもよいし、下記（１５）式で定義されるＣＮＲ（Contrast-to-Noise Rate）であってもよいし、下記（１６）式で定義されるＥＮＲ（Enhancement-to-Noise Ratio）であってもよいし、また、下記（１７）式で定義されるピーク値とバックグラウンドの差Ｐであってもよい。これらの指標は同様の傾向を示す。ｘは注目する領域内の強度の最大値であり、μは平均値であり、σは標準偏差である。これらの値は何れもリニアスケールの値を用いる。添え字のｘは注目する領域を示し、添え字のｂは信号のないノイズ領域を示す。
【００４５】
【数１４】
2014102087.tif 000016
【００４６】
【数１５】
2014102087.tif 000017
【００４７】
【数１６】
2014102087.tif 000018
【００４８】
【数１７】
2014102087.tif 000019
【００４９】
フィッティング品質については以下のとおりである。吸収がない非吸収性散乱体と、散乱および一定の吸収がある吸収性散乱体とを想定する。これら非吸収性散乱体および吸収性散乱体それぞれをサンプルとしてＯＣＴ測定を行ったときに得られる干渉信号（図９）を計算し、この干渉信号に様々なレベルのノイズを加える。図２に示されるように、分割した波長帯毎の断層画像を作成し、各断層画像の信号品質指標を求める。そして、平滑化パラメータα２を各値に設定して、上記（１１）式の目的関数が最小となるように吸収成分を求める。
【００５０】
ノイズが無い理想の条件であれば、非吸収性散乱体では吸収がゼロとなり、吸収性散乱体では設定した一定の吸収となる。しかし、ノイズが加えられた条件では、目的関数の最小化で得られる結果において、非吸収性散乱体および吸収性散乱体の何れの吸収スコアも或る幅を有する（図１０）。そこで、非吸収性散乱体および吸収性散乱体それぞれの吸収スコアの統計的な差を考慮して、ＲＯＣ（Receiver Operating Characteristic）曲線の面積（図１１）のハッチング領域の面積）をフィッティング品質と定義する。ＲＯＣ曲線とは、吸収スコアの境界値を変化させたときの偽陽性率（境界値以上のスコア値を持つ非吸収性散乱体の割合）を横軸にし、陽性率（境界値以上のスコア値を持つ吸収性散乱体の割合）を縦軸にして、プロットしたものである。ハッチング領域の面積が大きいほど、２つのスコア分布の重なりが小さいことを表し、フィッティング品質が良いことを表す。
【００５１】
図１２は、以上のようにして求めたフィッティング品質と平滑化パラメータとの関係を示すグラフである。同図には、信号品質指標としてCNR = 11およびCNR = 9それぞれの場合が示されている。同図から以下のことが判る。CNR = 11およびCNR = 9の何れの場合にも、平滑化パラメータα２の値が小さすぎるとフィッティング品質は悪くなり、平滑化パラメータα２の値が大きすぎてもフィッティング品質は悪くなり、フィッティング品質が良くなる平滑化パラメータα２の好適値が存在する。CNR = 11の場合と比べてCNR = 9の場合には、フィッティング品質は低い。また、CNR = 11の場合と比べてCNR = 9の場合には、平滑化パラメータα２の好適値が小さい。
【００５２】
このように、フィッティング品質が良くなる平滑化パラメータα２の好適値は、反射光量Ｙの信号品質指標との間に相関を有する。すなわち、反射光量Ｙの信号品質指標が高いほど、平滑化パラメータα２の設定値を大きくすれば、フィッティング品質が良くなって、サンプル内部の吸収分布または散乱分布を正確に求めることができる。本実施形態の光学的測定方法では、反射光量Ｙの信号品質指標の値に応じて平滑化パラメータα２の値を設定して、上記（１１）式の目的関数Ｓ'が最小となるように吸収成分および散乱成分を求める。これにより、サンプル内部の吸収分布または散乱分布を正確に求めることができる。
【００５３】
本発明は、上記実施形態に限定されるものではなく、種々の変形が可能である。上記の実施形態ではＯＣＴについて説明したが、本発明の光学的測定方法は、ＯＣＴに限られるものではなく、ＯＴＤＲおよびＤＯＴにも適用が可能である。また、本発明の光学的測定方法は、ＯＣＴのうちでもＴＤ-ＯＣＴおよびＦＤ-ＯＣＴの何れにも適用が可能である。
【符号の説明】
【００５４】
１…光学的測定システム、２…ＯＣＴ装置、３…プローブ、９…サンプル。
【図１】
2014102087.tif 000020
【図３】
2014102087.tif 000021
【図４】
2014102087.tif 000022
【図５】
2014102087.tif 000023
【図６】
2014102087.tif 000024
【図７】
2014102087.tif 000025
【図９】
2014102087.tif 000026
【図１０】
2014102087.tif 000027
【図１１】
2014102087.tif 000028
【図１２】
2014102087.tif 000029
【図２】
2014102087.tif 000030
【図８】
2014102087.tif 000031