(57)【要約】【課題】周期性ノイズを低減することができる光学的測定方法を提供する。【解決手段】本発明の光学的測定方法は、(1) サンプルからの反射光を検出し、複数の時間帯または波長帯におけるサンプル内部の反射位置と反射光量との関係を取得するデータ取得ステップと、(2) この関係に基づいて、入射光量に対する反射光量の減衰量に含まれる吸収成分および散乱成分を互いに分離して求めるデータ処理ステップと、を備える。データ処理ステップでは、(a) 反射光量，吸収成分，散乱成分，吸収成分の空間二次微分および散乱成分の空間二次微分を含む目的関数が最小となるように、吸収成分および散乱成分を求め、(b) 目的関数における吸収成分の空間二次微分および散乱成分の空間二次微分を求める際に、互いに異なる空間周波数特性を有する複数の二次微分演算子を吸収成分および散乱成分に作用させる。【選択図】図１
2014092438.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) このデータ取得ステップにおいて取得された複数の時間帯または波長帯におけるサンプル内部の反射位置と反射光量との関係に基づいて、入射光量に対する反射光量の減衰量に含まれる吸収成分および散乱成分を互いに分離して求め、サンプル内部の吸収分布または散乱分布を求めるデータ処理ステップと、を備えることを特徴とする。さらに、本発明の光学的測定方法は、データ処理ステップにおいて、(a) 反射光量，吸収成分，散乱成分，吸収成分の空間二次微分および散乱成分の空間二次微分を含む目的関数が最小となるように、吸収成分および散乱成分を求め、(b) 目的関数における吸収成分の空間二次微分および散乱成分の空間二次微分を求める際に、互いに異なる空間周波数特性を有する複数の二次微分演算子を吸収成分および散乱成分に作用させる、ことを特徴とする。
【０００７】
データ処理ステップにおいて、目的関数における吸収成分の空間二次微分および散乱成分の空間二次微分を求める際に、互いに異なる空間周波数特性を有する複数の二次微分演算子の線形結合で表される二次微分演算子を吸収成分および散乱成分に作用させるのが好適である。また、データ処理ステップにおいて、目的関数における吸収成分の空間二次微分および散乱成分の空間二次微分を求める際に、空間の位置によって空間周波数特性が異なる二次微分演算子を吸収成分および散乱成分に作用させるのも好適である。
【０００８】
吸収成分の空間二次微分または散乱成分の空間二次微分は、サンプル内部の深さ方向についての空間分解能当りに含まれるデータ数が１００以下であるのが好適である。データ取得ステップにおいて、サンプルにおいて生じる反射光と参照光とを干渉させた光を時間または波長の関数として検出し、この関数をフーリエ変換することでサンプル内部の反射位置と反射光量との関係を取得するのが好適である。また、データ取得ステップにおいて、サンプルにおいて生じる反射光と参照光とを干渉させた光を波長の関数として検出し、この関数を複数の波長帯に分割して波長帯毎にフーリエ変換することで波長帯毎のサンプル内部の反射位置と反射光量との関係を取得するのが好適である。
【発明の効果】
【０００９】
本発明によれは、周期性ノイズを低減することができる光学的測定方法を提供することができる。
【図面の簡単な説明】
【００１０】
【図１】本実施形態の光学的測定方法が適用される光学的測定システム１の概略構成を示す図である。
【図２】光学的測定システム１により取得されるデータ等の例を示す図である。
【図３】吸光係数εａ(λ)および散乱係数εｓ(λ)それぞれの波長依存性の例を示すグラフである。
【図４】行列Ｙの一例を模式的に示す図である。
【図５】行列Ａの一例を模式的に示す図である。
【図６】行列Ｌの一例を模式的に示す図である。
【図７】行列Ｘの一例を模式的に示す図である。
【図８】Savitzky-Golay filterを用いた数値微分により７点のデータを用いて二次微分を計算する際の重み係数を示すグラフである。
【図９】Savitzky-Golay filterを用いた数値微分により７点のデータを用いて二次微分を計算する際の重み係数をフーリエ変換して絶対値をとって求めた空間周波数特性を示すグラフである。
【図１０】Ｆａの二次元画像の一例を示す図である。
【図１１】Savitzky-Golay filterを用いた数値微分により７点のデータを用いて二次微分を計算する際の重み係数をフーリエ変換して絶対値をとって求めた空間周波数特性と、Ｆａの二次元画像のうち深さ方向に沿った或る直線上のデータをフーリエ変換して絶対値をとって求めた空間周波数特性と、を示すグラフである。
【図１２】本実施形態におけるＦａの二次元画像の一例を示す図である。
【図１３】Savitzky-Golay filterを用いた数値微分により７点のデータを用いて二次微分を計算する際の重み係数をフーリエ変換して絶対値をとって求めた空間周波数特性と、Savitzky-Golay filterを用いた数値微分により９点のデータを用いて二次微分を計算する際の重み係数をフーリエ変換して絶対値をとって求めた空間周波数特性と、を示すグラフである。
【発明を実施するための形態】
【００１１】
以下、添付図面を参照して、本発明を実施するための形態を詳細に説明する。なお、図面の説明において同一の要素には同一の符号を付し、重複する説明を省略する。
【００１２】
図１は、本実施形態の光学的測定方法が適用される光学的測定システム１の概略構成を示す図である。光学的測定システム１は、ＯＣＴ装置２およびプローブ３を備え、サンプル９の光断層画像を取得する。ＯＣＴ装置２は、光源から出力された光を２分岐して観察光および参照光として出力する。プローブ３は、ＯＣＴ装置２の光分岐部から出力された観察光を導光して先端からサンプル９へ照射するとともに、その光照射によって該サンプル９で生じる反射光を先端に入力してＯＣＴ装置２へ導光する。ＯＣＴ装置２は、プローブ３により導光されて到達した反射光と参照光とを干渉させて当該干渉光を検出する。そして、ＯＣＴ装置２は、この干渉光に基づいて、サンプルの深さ方向の光断層画像を取得することができる。また、プローブ３によるサンプル９への光照射の位置を走査（スキャン）することで、ＯＣＴ装置２は二次元または三次元の光断層画像を取得することができる。
【００１３】
本実施形態の光学的測定方法は、データ取得ステップおよびデータ処理ステップを備える。データ取得ステップでは、サンプルに光を照射または伝搬させたときに該サンプルにおいて生じる反射光を時間または波長の関数として検出して、複数の時間帯または波長帯におけるサンプル内部の反射位置と反射光量との関係を取得する。データ処理ステップでは、データ取得ステップにおいて取得された複数の時間帯または波長帯におけるサンプル内部の反射位置と反射光量との関係に基づいて、入射光量に対する反射光量の減衰量に含まれる吸収成分および散乱成分を互いに分離して求め、サンプル内部の吸収分布または散乱分布を求める。
【００１４】
図２は、光学的測定システム１により取得されるデータ等の例を示す図である。データ取得ステップにおいて、干渉光を検出することで、波長λに対する干渉光強度（信号強度）を得ることができる。この例では、波長帯を６分割して波長帯毎にフーリエ変換することで、波長帯毎の光断層画像を取得することができる。続く、データ処理ステップでは、この波長帯毎の光断層画像に基づいて以下のような処理を行う。
【００１５】
各波長帯の光断層画像は、該波長帯における散乱および吸収の双方の情報を含んでいる。ここで、ログスケールでのＯＣＴ断層像信号をＹ(λ,ｚ)とし、サンプル９における単位濃度当りの吸光係数をεａ(λ)とし、サンプル９における単位濃度当りの散乱係数をεｓ(λ)とし、サンプル９における吸収体の濃度のｚ＝０からの積分値をＦａ(ｚ)とし、サンプル９における散乱体の濃度のｚ＝０からの積分値をＦｓ(ｚ)とし、サンプル９における波長λに依存しない他の項のｚ＝０からの積分値をＣ(ｚ)とする。ｚは、サンプル９における深さ方向の位置を表す変数である。サンプル９の表面をｚ＝０とする。これらのパラメータの間の関係は下記（１）式で表される。
【００１６】
【数１】
2014092438.tif 000003
【００１７】
吸光係数εａ(λ)および散乱係数εｓ(λ)が既知であれば、未知数がＦａ(ｚ)、Ｆｓ(ｚ)およびＣ(ｚ)の３つであるので、３つ以上の波長帯で光断層像を得て連立方程式を解けば、Ｆａ(ｚ)およびＦｓ(ｚ)を得ることができる。そして、Ｆａ(ｚ)を微分することで、サンプル９における吸収体の空間分布を得ることができる。
【００１８】
吸光係数εａ(λ)は、ＦＴ-ＩＲによる実測値として得ることができる。散乱係数εｓ(λ)は、下記（２）式で表される波長λを変数とする一次関数のモデルを用いることができる。（２）式のλ０は、全波長帯の中心波長とする。図３は、吸光係数εａ(λ)および散乱係数εｓ(λ)それぞれの波長依存性の例を示すグラフである。（２）式を用いると、（１）式は下記（３）式のように表される。また、波長λを変数とする（３）式に替えて、波数ｋを変数とする下記（４）式を用いてもよい。（４）式で、波長λ０に対応する波数がｋ０である。
【００１９】
【数２】
2014092438.tif 000004
【００２０】
【数３】
2014092438.tif 000005
【００２１】
【数４】
2014092438.tif 000006
【００２２】
６分割後の各波長帯の波数をｋ１〜ｋ６とする。このとき、解くべき連立方程式は、下記（５）式および（６）式のように行列を用いて表現することができる。６行３列の行列Ａは、既知の行列要素からなる。３行１列の行列Ｘは、データ処理ステップにおいて求められるべき未知の行列要素からなる。６行１列の行列Ｙは、データ取得ステップにおいて得られたＯＣＴ断層像信号からなる。
【００２３】
【数５】
2014092438.tif 000007
【００２４】
【数６】
2014092438.tif 000008
【００２５】
（５）式の両辺に行列Ａの転置行列ＡＴを左から掛けると、下記（７）式となる。この（７）式から下記（８）式が得られる。このように行列表現を用いると、解Ｘは表される。
【００２６】
【数７】
2014092438.tif 000009
【００２７】
【数８】
2014092438.tif 000010
【００２８】
Ｆａ(ｚ)を求めるには、方程式の数（すなわち、分割した波長帯の数）は３であればよい。しかし、実際にはノイズの影響があるので、方程式の数が３であると、得られる解が不安定になり易い。そこで、方程式の数（すなわち、分割する波長帯の数）は、３より大きい６としている。ただし、波長帯の分割数が多すぎると光断層像の分解能が低下するので、波長帯の分割数を適切に設定することが好ましい。
【００２９】
未知変数の個数に対して方程式の個数が多いので、最小二乗法を用いて、下記（９）式で表される目的関数Ｓが最小（目的関数の偏微分がゼロ）となるような解Ｘを求める。しかし、実際に計算して解Ｘを求めてみると、光断層像の空間的なノイズが重畳された解Ｘが得られる。
【００３０】
【数９】
2014092438.tif 000011
【００３１】
このノイズ問題に対処する為、Ｆａ(ｚ)，Ｆｓ'(ｚ)それぞれがノイズに比べて空間的に滑らかに変化する（空間的な変化が小さい）ことを仮定して、目的関数を修正する。すなわち、下記（１０）式で表されるＦａ(ｚ)，Ｆｓ'(ｚ)それぞれの二次微分値が小さいと仮定する。修正後の目的関数Ｓ' は下記（１１）式で表される。αは、Ｓに対する二次微分値の重みを決めるパラメータである。目的関数Ｓ' は、反射光量Ｙ，吸収成分Ｆａ，散乱成分Ｆｓ'，吸収成分Ｆａの空間二次微分および散乱成分Ｆｓ'の空間二次微分を含む関数である。
【００３２】
【数１０】
2014092438.tif 000012
【００３３】
【数１１】
2014092438.tif 000013
【００３４】
行列を用いて表現すると下記（１２）式のようになる。Ｌは二次微分演算子（ラプラシアン演算子）である。この（１２）式から、解Ｘは下記（１３）式で表される。
【００３５】
【数１２】
2014092438.tif 000014
【００３６】
【数１３】
2014092438.tif 000015
【００３７】
Ｆａ(ｚ)，Ｆｓ'(ｚ)は離散的なデータであるので、或る位置での二次微分値は、該位置のデータに加えて前後のデータも用いた計算により求める。計算により求められる二次微分値は、計算手法によって異なり、また、計算に用いるデータの点数によっても異なる。計算手法としては差分法およびSavitzky-Golay filterを用いた数値微分が知られている。
【００３８】
図４は、行列Ｙの一例を模式的に示す図である。図５は、行列Ａの一例を模式的に示す図である。図６は、行列Ｌの一例を模式的に示す図である。図７は、行列Ｘの一例を模式的に示す図である。深さ方向のデータ数は５１２であり、波長帯の数は６である。差分法により３点のデータを用いて二次微分を計算する場合の行列Ｌを示している。二次微分演算子である行列Ｌは、Ｆａ(ｚ)，Ｆｓ'(ｚ)に作用するが、Ｄ(ｚ)には作用しない。
【００３９】
図８は、Savitzky-Golay filterを用いた数値微分により７点のデータを用いて二次微分を計算する際の重み係数を示すグラフである。この場合、７点のデータそれぞれに重み係数を乗じて足し合わせた結果の値が二次微分値となる。図９は、Savitzky-Golay filterを用いた数値微分により７点のデータを用いて二次微分を計算する際の重み係数をフーリエ変換して絶対値をとって求めた空間周波数特性を示すグラフである。
【００４０】
図１０は、Ｆａの二次元画像の一例を示す図である。ここでは、７点のデータを用いたザSavitzky-Golay filterを用いた数値微分による二次微分演算子を採用した。同図に示されるように、Ｆａの二次元画像は、深さ方向に周期的なノイズが重畳されている。図１１は、Savitzky-Golay filterを用いた数値微分により７点のデータを用いて二次微分を計算する際の重み係数をフーリエ変換して絶対値をとって求めた空間周波数特性と、Ｆａの二次元画像のうち深さ方向に沿った或る直線上のデータをフーリエ変換して絶対値をとって求めた空間周波数特性と、を示すグラフである。
【００４１】
同図に示されるように、Ｆａの二次元画像（図１０）において重畳されている周期的なノイズは、二次微分の空間周波数特性の零点に相当する周波数を有していることが判る。
【００４２】
そこで、本実施形態では、目的関数Ｓ' における吸収成分の空間二次微分および散乱成分の空間二次微分を求める際に、互いに異なる空間周波数特性を有する複数の二次微分演算子を吸収成分および散乱成分に作用させる。好適には、互いに異なる空間周波数特性を有する複数の二次微分演算子の線形結合で表される二次微分演算子を吸収成分および散乱成分に作用させる。線形結合には平均だけでなく加重平均も含まれる。または、空間の位置によって空間周波数特性が異なる二次微分演算子を吸収成分および散乱成分に作用させるのも好適である。
【００４３】
図１２は、本実施形態におけるＦａの二次元画像の一例を示す図である。ここでは、７点のデータを用いたSavitzky-Golay filterを用いた数値微分による二次微分演算子Ｌ７と、９点のデータを用いたSavitzky-Golay filterを用いた数値微分による二次微分演算子Ｌ９と、の線形結合で表される二次微分演算子を採用した。図１０と比べると、図１２に示されるＦａの二次元画像では、深さ方向に周期的なノイズの発生が抑制されている。
【００４４】
Savitzky-Golay filterを用いた数値微分および差分法の何れの場合にも、計算に用いるデータ点数に依存した零点が空間周波数特性に現れる。Savitzky-Golay filterを用いた数値微分および差分法の何れの場合にも、データ点数を増やすと零点の数が増えて、零点と零点との間の平均的間隔が狭まる。また、データ点数によって零点の位置は異なる。
【００４５】
図１３は、Savitzky-Golay filterを用いた数値微分により７点のデータを用いて二次微分を計算する際の重み係数をフーリエ変換して絶対値をとって求めた空間周波数特性と、Savitzky-Golay filterを用いた数値微分により９点のデータを用いて二次微分を計算する際の重み係数をフーリエ変換して絶対値をとって求めた空間周波数特性と、を示すグラフである。２つの二次微分演算子Ｌ７，Ｌ９それぞれの空間周波数特性は、互いに異なっており、互いに異なる周波数で零点を有している。
【００４６】
このような空間周波数特性の異なる複数の二次微分演算子を用いることにより、深さ方向に周期的なノイズの発生が抑制されたＦａ，Ｆｓ'を得ることができる。そして、Ｆａ，Ｆｓ'を微分することで、サンプルにおける吸収分布、散乱分布を得ることができる。
【００４７】
本発明の方法である空間周波数特性の異なる複数（ｎ個）の二次微分演算子を組み合わせる方法では、複数の二次微分演算子をＬ１（i = 0,1,2,・・・,n）、線形結合または空間位置ごとのそれぞれの二次微分演算子の重み係数を決める行列をＡｉ（i = 0,1,2,・・・,n）とすると、Ａｉはｎ次の対角行列によって表され、前記（１２）式のＬとして、（１４）式で表される行列を用いれば良い。また、複数の二次微分演算子を組み合わせる方法として、(１２)式のＬＴＬとして、(１５)式で表される行列を用いても良い。
【００４８】
【数１４】
2014092438.tif 000016
【００４９】
【数１５】
2014092438.tif 000017
【００５０】
上記実施形態では、Savitzky-Golay filterを用いた数値微分の場合について記載しているが、本発明の方法は、他の数値微分手法を用いる場合にも適用できることは、当業者には明白である。そのような他の手法としては、前進差分近似、後退差分近似、中心差分近似、Stirlingの公式、などが知られている。
【００５１】
本発明は、上記実施形態に限定されるものではなく、種々の変形が可能である。上記の実施形態ではＯＣＴについて説明したが、本発明の光学的測定方法は、ＯＣＴに限られるものではなく、ＯＴＤＲおよびＤＯＴにも適用が可能である。また、本発明の光学的測定方法は、ＯＣＴのうちでもＴＤ-ＯＣＴおよびＦＤ-ＯＣＴの何れにも適用が可能である。
【符号の説明】
【００５２】
１…光学的測定システム、２…ＯＣＴ装置、３…プローブ、９…サンプル。
【図１】
2014092438.tif 000018
【図３】
2014092438.tif 000019
【図４】
2014092438.tif 000020
【図５】
2014092438.tif 000021
【図６】
2014092438.tif 000022
【図７】
2014092438.tif 000023
【図８】
2014092438.tif 000024
【図９】
2014092438.tif 000025
【図１１】
2014092438.tif 000026
【図１３】
2014092438.tif 000027
【図２】
2014092438.tif 000028
【図１０】
2014092438.tif 000029
【図１２】
2014092438.tif 000030