JP2018179529A - 応答特性取得装置、応答特性取得方法および応答特性取得プログラム - Google Patents

Abstract

【課題】少ない測定回数で高精度に応答特性行列を生成する。
【解決手段】応答特性取得装置は、光源１００からの光の一部の位相を変調することで互いに位相差を有する第１の光と第２の光を生成する位相変調手段１７０と、第１の光と被検体に照射された第２の光との干渉により生じる信号の測定を行う測定手段２３０と、位相シフト法により上記信号の位相分布を回復し、第１および第２の光の被検体への入射条件が異なる上記信号から回復した複数の位相分布から応答特性行列を生成する処理手段２４０とを有する。処理手段は、応答特性行列の複数の要素としての位相の回復誤差を算出し、該回復誤差に基づいて少なくとも１つの要素としての位相の補正を行う。
【選択図】図２

Description

本発明は、被検体の応答特性を取得する技術に関する。

生体等の被検体である散乱媒質に対して、空間光変調器（ＳＬＭ：Spatial Light Modulator）により様々な波面を成形した光を照射し、媒質を透過した散乱光を撮像素子を通して観察する方法が非特許文献１に開示されている。また、非特許文献１にて開示された方法では、ＳＬＭにより光の位相を変調する変調領域と変調しない非変調領域とを分離し、変調領域に非変調領域に対する位相差を与えて位相シフト法を実施することで、散乱光の位相分布を測定する。様々な入射波面に対して散乱光の位相分布を測定することで、媒質の光に対する入出力応答を表す応答特性行列（透過行列）を取得することができる。この透過行列を用いれば、媒質の背後の任意の位置に光をフォーカスすることが可能となる。

また、非特許文献２には、媒質を透過した散乱光の代わりに、媒質中で発生した光音響信号の応答特性行列を測定する方法が開示されている。応答特性行列を利用して媒質内部の吸収物体の位置に選択的に光を照射し、高ＳＮＲ（Signal-to-Noise Ratio）にて光音響信号を測定することができる。

さらに、特許文献１には、マルチモードファイバの応答特性行列を測定し、該行列を利用してファイバの出射端に生成したフォーカススポットのスキャンを行うことで、ファイバを介して被検体の像を取得する装置が開示されている。加えて、非特許文献３には、被検体を透過または反射した光の位相分布を測定または回復する様々な方法や装置が開示されている。

前述した被検体の応答特性行列を測定するためには、被検体から発せられる信号（透過光、反射光または光音響信号等）の位相を測定する必要がある。この位相測定には、光の干渉計測を行うのが一般的である。例えば、マッハツェンダー干渉計に代表されるように、光源からの光を物体光と参照光とに分離し、前者を被検体に入射させ、その透過光を観察面で参照光と重ね合わせて光強度を測定する。このような測定系において、非特許文献３にて開示されているように位相シフト法を行うことで、被検体を透過した光の位相を測定することができる。

位相シフト法では、物体光または参照光に対して既知の位相遅延量を等間隔に複数回与え、それぞれの位相遅延量で測定した光強度を用いて、非特許文献３に開示された様々なアルゴリズムを用いて位相を回復することができる。また、非特許文献１にて開示されている測定装置のように、物体光と参照光の光路が、被検体への入射から被検体から出射して撮像素子により観察されるまで共通である測定系においても位相シフト法を適用することができる。このようなコモンパスの測定系は、前述の測定系に比べて振動に強く、装置の構成も比較的簡易である。また、非特許文献２にて開示されているように、光音響信号の応答特性行列を測定することも可能である。

米国特許公開第２０１３／０１８２２５３号公報

S. M. Popoff et al., "Measuring the Transmission Matrix in Optics: An Approach to the Study and Control of Light Propagation in Disordered Media", Phys. Rev. Lett., 104(10), 100601 (2010) T. Chaigne et al., "Controlling light in scattering media non-invasively using the photoacoustic transmission matrix", Nat. Photonics, 8, 58 (2014) D. Malacara, OPTICAL SHOP TESTING 3rd Edition, WILEY Y. Choi et al., "Scanner-free and wide-field endoscopic imaging by using a single multimode optical fiber", Phys. Rev. Lett., 109, 203901 (2012)

前述した位相シフト法による位相回復において、与える既知の位相遅延量（位相シフト量）が観測面で観測する際に誤差を生じる場合、既知の位相遅延量を前提として位相を回復（算出）するアルゴリズムを用いると、回復される位相に誤差が生じ精度が劣化する。特に、被検体が散乱媒質であり、測定系が前述したコモンパスである場合は、被検体への光の入射前に与えた位相シフト量が被検体を透過した後にランダムに乱される。この結果、観測面での位相シフト量は、既知として与えた量に対して大きな誤差が生じ、これが位相回復の精度を劣化させる。このため、回復された位相を要素に持つ応答特性行列の精度が劣化する。これにより、前述した媒質の背後でフォーカスされた光の最大強度が低下したり、応答特性行列を用いて再構成される画像の精度が低下したりする等、応答特性行列を利用した様々な機能の精度が低下する。

本発明は、少ない測定回数で高精度に応答特性を取得することが可能な装置等を提供する。

本発明の一側面としての応答特性取得装置は、光源からの光の一部の位相を変調することで互いに位相差を有する第１の光と第２の光を生成する位相変調手段と、第１の光と被検体に照射された第２の光との干渉により生じる信号の測定を行う測定手段と、位相シフト法により上記信号の位相分布を回復し、第１および第２の光の被検体への入射条件が異なる上記信号から回復した複数の位相分布から応答特性行列を生成する処理手段とを有する。処理手段は、応答特性行列の複数の要素としての位相の回復誤差を算出し、該回復誤差に基づいて少なくとも１つの要素としての位相の補正を行うことを特徴とする。

また、本発明の他の一側面としての応答特性取得方法は、光源からの光の一部の位相を変調することで互いに位相差を有する第１の光と第２の光を生成するステップと、第１の光と被検体に照射された第２の光との干渉により生じる信号の測定を行うステップと、位相シフト法により上記信号の位相分布を回復し、第１および第２の光の被検体への入射条件が異なる上記信号から回復した複数の位相分布から応答特性行列を生成する処理ステップとを有する。処理ステップは、応答特性行列の複数の要素としての位相の回復誤差を算出し、該回復誤差に基づいて少なくとも１つの要素としての位相の補正を行うことを特徴とする。

なお、上記応答特性取得方法に対応する処理をコンピュータに実行させるコンピュータプログラムも、本発明の他の一側面を構成する。

本発明によれば、少ない測定回数で高精度に応答特性を取得することができる。

本発明の実施例に共通する処理を示すフローチャート。 実施例１の装置の構成を示す図。 実施例１の装置の構成の一部を示す図。 実施例１における応答特性行列生成処理を示すフローチャート。 誤差マップ行列の全要素（ａ）と、誤差が大きい順に並び替えて抽出した一部（ｂ）を示す図。 実施例１における再生フォーカス像の設定を示す図。 実施例１における応答特性行列の補正処理による再生フォーカス像のエンハンスメントを比較する図。 実施例２の装置の構成を示す図。 実施例３の装置の構成を示す図。 実施例３の装置の構成におけるイメージングを示す図。

以下、本発明の実施例について図面を参照しながら説明する。本発明の各実施例は、散乱媒質である被検体の応答特性行列を生成する。応答特性行列ｔは、次式（１）のように、被検体に入力される電場をＥ_ａとし、被検体から出力（例えば透過）される電場をＥ_ｂとするときに、被検体の線形的な入出力応答として定義される。

式（１）に基づいて様々な入力（基底）のＥ_ａに対して出力Ｅ_ｂを測定することで、被検体の応答特性行列ｔを得ることができる。この応答特性行列を利用して入力Ｅ_ａを適切に調整すれば、被検体からの出力（散乱光）を自由に制御することができる。例えば、被検体の背後に光のフォーカススポットを生成することも可能である。被検体の例として、生体組織がある。可視域から近赤外の光を用いて生体内部の光学特性を測定およびイメージングする場合、光は強く散乱されるため、解像度や測定深度（侵達長）が低下する。この課題に対して上述した応答特性行列を利用した散乱光の制御を行えば、生体背後や内部のある特定の位置に光を効率的に送り込むことができる。これにより、生体光計測における散乱に起因する課題を解決することができる。

図１のフローチャートは、応答特性行列の測定から該応答特性行列を利用した被検体の測定およびイメージングを行う場合の処理の例を示している。パーソナルコンピュータ（以下、ＰＣという）がコンピュータプログラムに従って本処理を実行する。Ｓはステップを意味する。

まず、ＰＣは、Ｓ３１０およびＳ３２０において、後述する方法により被検体の応答特性行列を生成する。次に、ＰＣは、Ｓ３３０において応答特性行列の生成結果を解析し、Ｓ３４０において応答特性行列を利用して散乱光を制御した被検体への光照射（例えば、フォーカススポットの生成）を行う。そして、Ｓ３５０で被検体からの信号を取得し、必要に応じて該信号から画像を生成する。以下、具体的な実施例について説明する。

本発明の実施例１である応答特性取得装置および応答特性取得方法について説明する。図２は、本実施例の応答特性取得装置の構成を示している。

光源１００からは所望のビームサイズでコリメートされた光１６０が出力される。光源１００は、例えば４００〜１５００ｎｍの可視域から近赤外域の波長を有し、強度が時間的に一定な連続光（ＣＷ光）を放射するレーザ光源である。光源１００からの光は、目的に応じて、上記範囲以外の波長を有してもよいし、任意の周波数で強度変調された光やパルス光であってもよい。また、光源１００からの光は、コヒーレンス長が十分に長い（例えば、数１０ｃｍ以上）光であることが望ましい。

光源１００からの光は、空間フィルタ１１０およびレンズ１２０によって平行光に変換され、可変絞り１３０によりビームサイズが適切に調整される。また、光源１００からの光の強度は、被検体２００に応じて又は目的に応じて適切に調整されている。ビームサイズと光強度が調整された光１６０は、ミラー１４０，１５０で反射されて位相変調手段としてのＳＬＭ１７０に入射する。

ＳＬＭ１７０としては、例えばliquid crystal on silicon（ＬＣＯＳ）を用いることができる。ＳＬＭ１７０は位相変調を行うデバイスであり、反射型のデバイスであってもよいし透過型のデバイスであってもよい。また、位変調手段としてＤＭＤ（Digital Mirror Device）を用いて、バイナリパターンによる回折像を利用して位相変調を行ってもよい。ＳＬＭ１７０による位相変調を受けた光はその波面が成形される。

ＳＬＭ１７０は、処理手段としてのＰＣ２４０に接続されており、ＰＣ２４０からの制御信号に応じて制御される。ＳＬＭ１７０に入力する光１６０の偏光は、ＳＬＭ１７０による位相変調が行われる偏光方向と一致するように調整されている。

ＳＬＭ１７０で反射し、かつ波面が成形された光１９０は、レンズ１８０を介して被検体２００に入射する。このとき、レンズ１８０とＳＬＭ１７０との間隔およびレンズ１８０と被検体２００の入射面との間隔はいずれもレンズ１８０の焦点距離に等しい。したがって、ＳＬＭ１７０と被検体２００の入射面はフーリエ変換の関係にあり、ＳＬＭ１７０で成形された波面をフーリエ変換した分布の光が被検体２００に入射される。

被検体２００を透過して出射した散乱光２１０は、レンズ２２０を介してＣＣＤセンサ等の撮像デバイス（測定手段：以下、単にＣＣＤという）２３０により測定（検出）される。撮像デバイスとして、ＣＭＯＳセンサ、イメージインテンシファイアを有するアレイセンサ、ＥＭＣＣＤ、ｓＣＭＯＳ等を用いることもできる。また、ＣＣＤ２３０は、ＰＣ２４０に接続され、ＰＣ２４０からの制御信号に応じて制御される。ＣＣＤ２３０による測定により得られた測定データは、ＰＣ２４０に転送されて解析される。測定データおよびその解析結果は、表示部であるモニタ２５０に表示される。

ＳＬＭ１７０の位相変調領域の画素を用いて応答特性行列の基底を表現し、式（１）における入力Ｅ_ａとする。基底としてアダマール基底（Hadamard basis）を用いる。アダマール基底は、全ての要素が＋１または−１の要素から構成されることを特徴とする直交基底である。この特徴を利用して、基底の一方の要素（例えば＋１）に対して、他方の要素に位相差πを与えることで、位相変調を行うＳＬＭ１７０を用いてアダマール基底を生成することができる。このアダマール基底をＳＬＭ１７０に順次設定する。

次に、被検体２００から出射する透過光（信号）の位相測定について説明する。図３は、ＳＬＭ１７０の有効領域（位相変調可能な全画素領域）１７１を模式的に示している。ＳＬＭ１７０の有効領域１７１の中心に光源からの光１６０が照射されるようにアライメントがなされている。ここで、アダマール基底を描画して位相変調する位相変調領域１７３を設定する。位相変調領域１７３の中心は、有効領域１７１および入射光ビーム１６０の中心に一致し、入射光ビーム１６０の照射領域よりも小さく設定されている。この位相変調領域１７３の外側で、かつ入射光ビーム１６０が照射されている領域１７２を、位相変調を行わない非位相変調領域１７２とする。

本実施例では、位相変調領域１７３に、非位相変調領域１７２に対する位相差を与えることで、位相シフト法を行う。この際、位相変調領域１７３で位相変調される光を物体光とし、非位相変調領域１７２で位相変調されずに反射される光を参照光とすることで、入射光ビーム１６０を物体光と参照光とに分割する。

図１に示したＳ３１０において、ＰＣ２４０は、物体光の位相を参照光の位相に対して位相ステップ量αでシフトさせながら、各位相シフト量（位相差）に応じてＣＣＤ２３０で透過光２１０の強度分布を測定する。この位相シフト法を用いた光強度の測定を、Ｎ個のアダマール基底に対して順次行う。測定が終了すると、光源１００の出力を停止する。

次に図１に示したＳ３２０において、ＰＣ２４０は、測定された結果を用いて、後述する方法により各基底で測定した透過光の位相分布を回復し、回復した位相分布を列ベクトルとして行列にスタックすることで応答特性行列を生成する。ＣＣＤ２３０で測定に用いられる画素数をＭとするとき、Ｍ行Ｎ列の応答特性行列ｔが得られる。例えば、応答特性行列のｊ列目の成分は、ｊ番目のアダマール基底で測定した位相分布に対応する。

以下、図１のＳ３２０でＰＣ２４０が行う応答特性行列生成処理について、図４のフローチャートを用いて詳細に説明する。ＰＣ２４０は上述したコンピュータプログラムに含まれる応答特性取得プログラムに従って本処理を実行する。

まずＳ３２１において、ＰＣ２４０は、Ｓ３２１〜Ｓ３２３の処理がＮ個の全ての基底（要素）に対して完了したか否かを確認する。処理が完了していない場合は、処理が未完了のある基底ｊ（１≦ｊ≦Ｎ）の測定結果を選択し、Ｓ３２２において位相回復処理を行う。この位相回復処理は、位相ステップ量αで測定した結果から被検体の透過光の位相分布を回復する処理であり、非特許文献３にて説明されている。非特許文献３によると、ｍ回目（ｍ≧１）の位相ステップで測定される光強度Ｉ_ｍ（ｘ，ｙ）は以下の式（２）のように表せる。

Ｉ′（ｘ，ｙ）は強度が一定のバイアス成分であり、Ｉ″（ｘ，ｙ）は物体光と参照光との干渉によって生じる変調成分（変調度）である。（ｘ，ｙ）はＣＣＤ２３０の撮像面上での離散的な座標値である。φ（ｘ，ｙ）が回復したい位相分布である。
本実施例では、Ｓ３１０の測定において、位相ステップ量αをπ／２、測定回数ｍを４として光強度の測定を行い、Ｉ_１（ｘ，ｙ）、Ｉ_２（ｘ，ｙ）、Ｉ_３（ｘ，ｙ）およびＩ_４（ｘ，ｙ）の４つの測定データ（透過像データ）を得る。そして、この測定をＳＬＭ１７０に設定する基底を変えながら、つまりは物体光と参照光が被検体２００に入射するときの入射条件を変えながらＮ個の基底のそれぞれに対して行う。

各基底に対してＩ_１（ｘ，ｙ）〜Ｉ_４（ｘ，ｙ）の光強度の測定結果が得られた場合に、これら４つの測定結果を用いて、次式（３）のように位相分布φ（ｘ，ｙ）を回復する。

また、バイアス成分Ｉ′（ｘ，ｙ）および変調成分Ｉ″（ｘ，ｙ）についても、次式（４ａ）および（４ｂ）の計算により求めることができる。

次にＳ３２２において、ＰＣ２４０は、位相分布φ（ｘ，ｙ）の回復に加えて、バイアス成分Ｉ′（ｘ，ｙ）および変調成分Ｉ″（ｘ，ｙ）についても算出する。

次にＳ３２３において、ＰＣ２４０は、位相分布を回復した結果について以下の式（５）を用いて回復誤差を評価する。具体的には、式（３），（４ａ）および（４ｂ）による算出結果を式（２）に代入し、式（２）で表される推定値としてのコサインカーブと測定結果との差分をとって誤差を算出する。そして、該誤差を変調度Ｉ″（ｘ，ｙ）で重み付けすることで、誤差評価値としての回復誤差ε（ｘ，ｙ）を算出する。

ただし、Ｉ_k（ｘ，ｙ）はｋ番目（１≦ｋ≦ｍ）の測定結果であり、本実施例では、ｍ＝４、α＝π／２である。このようにして各測定位置（ｘ，ｙ）での回復誤差εを評価する。変調度Ｉ″（ｘ，ｙ）で重み付けした回復誤差εを評価することで、単純に誤差が大きいものを評価するのではなく、差分が大きくかつ物体光と参照光との干渉の寄与が大きい回復誤差εから順に優先して評価する。

なお、回復誤差の算出式は、式（５）に限定されず、式（２）による位相シフト法での回復結果と測定結果との差分を評価できるものであればどのような算出式であってもよい。ただし、変調度で重み付けして干渉の寄与が大きいものを優先することが望ましい。

ＰＣ２４０は、以上説明したＳ３２２〜Ｓ３２３の処理を、Ｓ３２１においてＮ個の全ての基底に対して完了したと判定するまで反復する。Ｓ３２１で完了を判定するとＳ３２４に進む。

Ｓ３２４では、ＰＣ２４０は、それぞれの基底ｊで回復した２次元の位相分布φ_ｊ（ｘ，ｙ）を１次元の列ベクトルにして、順に行列としてスタックしていくことで、前述したＭ行Ｎ列の応答特性行列ｔを生成する。また、ＰＣ２４０は、各基底ｊで算出（評価）した回復誤差ε_ｊ（ｘ，ｙ）についても同様に１次元の列ベクトルにして順に行列としてスタックして、同じＭ行Ｎ列のサイズの誤差マップ行列ｔ_εを生成する。誤差マップ行列ｔ_εは、応答特性行列ｔの各要素と対応する位置の要素に、式（５）で評価した回復誤差εが記録されている。したがって、誤差マップ行列ｔ_εに基づいて応答特性行列ｔのどの要素が位相回復精度が低いかを判断することができる。

Ｓ３２５以降では、ＰＣ２４０は、誤差マップ行列ｔ_εに基づいて、暫定的に求めた応答特性行列ｔの要素、すなわち位相の補正を行う。

誤差マップ行列ｔ_εの例として、被検体２００をガラス拡散板として実測した結果を図５（ａ），（ｂ）に示す。図５（ａ）は、誤差マップ行列ｔ_εの全要素を横軸にとり、式（５）で評価した回復誤差εを縦軸にとって回復誤差εをプロットした図である。なお、アダマール基底の数はＮ＝２５６、ＣＣＤ２３０の観測画素数はＭ＝２５６である。

図５（ａ）から、応答特性行列の要素によって回復誤差εが大きくばらついていることが分かる。ここで、閾値（所定値）を任意に設定し、該閾値以上の回復誤差εを抽出する。

図５（ｂ）は、抽出した回復誤差εの一部を大きい順に並べ替えてプロットした図である。Ｓ３２５では、ＰＣ２４０は、このように誤差マップ行列ｔ_εから閾値を設定して、該閾値以上の回復誤差を有する少なくとも１つの要素を補正すべき要素（以下、補正対象要素という）として抽出して誤差補正リストを作成する。誤差補正リストは、図５（ｂ）に示すように、回復誤差が大きい順に補正対象要素を並び替え、この順で以下に説明する位相補正を行う。

Ｓ３２６では、ＰＣ２４０は、補正の終了条件を満たしているか否かを確認する。終了条件は、Ｓ３２５で作成したリストにある全ての補正対象要素に対する補正を完了したことである。全補正対象要素に対する補正が終了すれば、ＰＣ２４０は本処理を終了する。一方、Ｓ３２６においてまだ全補正対象要素に対する補正が終了していない場合は、ＰＣ２４０は、誤差補正リストのより上位の補正対象要素から順にＳ３２７〜Ｓ３２９に従って補正対象要素の位相補正を行う。

Ｓ３２７では、ＰＣ２４０は、補正対象要素について、位相探索範囲内でＳ３２８〜Ｓ３２９の処理が完了したか否かを確認する。位相探索範囲は０〜２πでもいし、Ｓ３２２で回復した位相分布に対して適切な±δφを加えた範囲を設定してもよい。例えば、δφは位相ステップ量α以下に設定する。

Ｓ３２８では、ＰＣ２４０は、上記位相探索範囲において、応答特性行列ｔの要素である位相を補正すべき位相に設定し、その位相を補正する。Ｓ３２９では、ＰＣ２４０は、補正された位相が正解に近いか否かを評価する。具体的には、ＰＣ２４０は、応答特性行列ｔを利用して、被検体２００の背後のＣＣＤ２３０にフォーカススポットを生成する。フォーカススポットの生成は、まず、図６に模式的に示すように、ＣＣＤ２３０で再生したい電場の振幅分布４５０（ターゲットＥ_{ｔａｒｇｅｔ}）を設定する。図６の中心に示す白い領域がフォーカススポット領域であり、周囲の黒い部分は振幅０の領域である。また、再生するフォーカススポットの光強度を最大限に向上させるために、このフォーカススポット領域はＳ３１０においてＣＣＤ２３０で観測されるスペックル粒サイズ以下に設定することが望ましい。

次に、ＰＣ２４０は、同スポット分布が得られるためにＳＬＭ１７０で設定すべき入力波面（位相分布）Ｅ_ａを以下の式（６）のように算出する。

ここで、†は転置共役を示す。式（６）で表現される入力波面Ｅ_ａは、図６の電場分布がＣＣＤ２３０から射出して、測定とは逆方向に伝搬し、被検体２００を通ってＳＬＭ１７０まで到達した波面の位相共役波に相当する。ＰＣ２４０は、この入力波面Ｅ_ａの位相分布をＳＬＭ１７０に設定する。

次に、ＰＣ２４０は、再度光源１００の出力をオンとし、ＳＬＭ１７０に設定されている位相分布により再生された入力波面を有する光を被検体２００に照射する。被検体２００を透過した散乱光は、ＣＣＤ２３０上において、設定されたフォーカススポット領域にフォーカスされる。この再生フォーカス像をＣＣＤ２３０で測定する。Ｓ３２９における評価は、次式（７）で定義される、フォーカスによる光強度の向上（エンハンスメント）ηを用いて行う。

ここで、Ｉ_ｍａｘはフォーカススポットでの散乱光の最大強度である。＜Ｉ_０＞は、ＳＬＭ１７０にフラットな位相分布（位相共役波面ではない光であり、例えば、ｊ＝１のアダマール基底）を設定したしたときの透過像（バックグランド）の平均強度を表す。バックグランドの測定は、Ｓ３１０で測定した結果を用いてもよいし、Ｓ３２９の評価時に測定してもよい。

このエンハンスメントηは、応答特性行列の精度を定量的に示す評価値として考えることができる。すなわち、応答特性行列の精度が高ければ、算出される位相共役波面の精度も高く、再生されるフォーカス像の最大強度は向上する。したがって、エンハンスメントηが高ければ応答特性行列の精度が高いことを表している。

このことを利用して、ＰＣ２４０は、Ｓ３２７〜Ｓ３２９において補正対象要素の位相を順次補正して（Ｓ３２８）、再生フォーカス像を評価し（Ｓ３２９）、エンハンスメントηが位相探索範囲内で最大となる位相を探索する。Ｓ３２７で位相探索範囲内の探索が完了すると、ＰＣ２４０は、Ｓ３２７ｂにおいて、上述のエンハンスメントηが最大となる位相を補正前の応答特性行列の要素の位相と置換する。図７は、誤差補正リストの１番上の要素について、上述した補正を行う前（図７（ａ））と補正を行った後（図７（ｂ））にフォーカス像を再生した結果の１次元断面をプロットした結果である。図中のηは上述したエンハンスメントの値であり、エンハンスメントηがＳ３２７〜Ｓ３２９までの処理によって向上していることが確認できる。

Ｓ３２６で終了条件が満足すると、ＰＣ２４０は本処理を終了するとともに、要素が補正された行列を補正後の応答特性行列ｔとして出力する。

このように、本実施例では、エンハンスメントηを応答特性行列の測定精度の評価値として利用し、評価値の値が向上するように応答特性行列の要素（位相）を補正する。これは、上述したように、応答特性行列を利用した再生結果（フォーカススポット）が、応答特性行列の精度に直結するためである。再生像のフォーカススポットを評価する評価値として、エンハンスメントη以外にも、フォーカススポットの最大強度、形状およびサイズ等を評価してもよい。応答特性行列の精度が向上すれば、スポット形状は、歪でランダムな形状からガウス分布に近づき、スポットサイズはスペックル粒サイズに近づく。また、Ｓ３２７〜Ｓ３２９での補正すべき位相の探索において、上述したような逐次探索する方法以外にも、最急降下法等の勾配法を用いた最適化や、遺伝的アルゴリズム等の最適化アルゴリズムを用いてもよい。

応答特性行列の要素は、上述の位相シフト法による位相分布を回復した結果である。したがって、行列の要素の精度は、位相分布の回復精度である。前述したように、被検体が散乱媒質であり、かつ本実施例のように測定系がコモンパスてある場合は、位相シフト法で既知として与える位相シフト量と、観測位置で測定される位相シフト量とに誤差が生じる。この条件下で、与えた位相シフト量に基づいて位相分布を回復すると位相分布の回復精度が低下する。位相分布の回復精度は、通常は、測定回数を増やせば向上する。例えば、Ｎ個の基底の全てにおいて位相シフト法の測定回数を増やすことでも応答特性行列の測定精度を向上させることができる。ただし、この場合は測定回数が増える分、測定時間もより多く要する。

これに対して本実施例では、できるだけ少ない測定回数（例えば、ｍ＝３〜５）の位相シフト法で波面を回復する。そして、得られた結果から位相分布の回復精度を評価し、回復誤差が大きい応答特性行列の要素に注目して、効率的にその要素（位相）を補正する。これにより、Ｎ個の基底全てで測定回数を増やして応答特性行列を測定する場合に比べて、少ない回数で効率的に精度良く応答特性行列を生成することが可能である。測定回数を減らすことで、応答特性行列を取得する時間を短縮することができる。

被検体２００の散乱特性が時間によって変化する場合、応答特性行列の測定時と再生時とで散乱特性が変化するため、上述した再生フォーカス等、応答特性行列を利用した効果が得られなくなる。このような観点から、できるだけ少ない測定回数で高精度に応答特性行列を生成することが重要となる。

本実施例では、応答特性行列の基底としてアダマール基底を用いて説明したが、用いる基底はこれ以外の基底であってもよい。なお、使用する基底は、直交していることが望ましい。また、Ｓ３２２の位相回復処理も、測定回数ｍに応じて、非特許文献３にて開示された様々なアルゴリズムを用いて行うことができる。例えば、任意の測定回数ｍ（ｍ≧３）に対して最小二乗法を用いて位相分布を回復してもよい。

また、本実施例では、物体光と参照光が共通の光路を伝搬するコモンパスの測定系を用い、位相シフト法による光強度の測定（Ｓ３１０）と再生フォーカス像の評価（Ｓ３２９）とを同じ測定系で容易に実施することができる。

また、本発明の実施例は、被検体２００の透過光２１０についての位相分布の測定に限定されず、被検体２００から反射した反射光（または後方散乱光）についても同様に位相分布を回復したり、さらには応答特性行列を取得したりすることもできる。

上述のようにして図１のＳ３２０で応答特性行列を取得したＰＣ２４０は、補正後の応答特性行列を利用して散乱光を制御した光照射と測定を行うことができる。例えば、Ｓ３３０で応答特性行列を特異値分解し、透過率（または反射率）が高い上位の特異ベクトルの入射波面の光を被検体に照射し（Ｓ３４０）、その応答をＳ３５０で測定してもよい。また、補正後の応答特性行列を用いて任意の位置にフォーカススポットを再生し、同フォーカススポットをスキャンすることで画像を生成してもよい。

本発明の実施例２である応答特性取得装置および応答特性行列生成定方法について説明する。図８は、本実施例の応答特性取得装置としての光音響イメージング装置の構成を示している。

被検体５２０は、生体組織であり、散乱粒子５２１を含む。この被検体５２０は、可視域から近赤外域の光に対して散乱媒質である。光源５００からは可視域から近赤外域の波長を有する数ｎｓのパルス光が放射される。また、光源５００は、被検体５２０の主要な構成成分である水、脂肪、タンパク質、酸化ヘモグロビンおよび還元ヘモグロビン等の吸収スペクトルに応じた複数の波長を選択することができる。 光源５００から発せられた光は、ビームスプリッタ５０１を透過してＳＬＭ５０２に入射する。ＳＬＭ５０２は、ＰＣ等により構成される制御部（処理手段）５０５によって制御される。ＳＬＭ５０２に入射した光は、実施例１（図３）と同様に、その一部が位相変調を受けて残りが位相変調されずに反射される。ＳＬＭ５０２で反射した光５１０は、ビームスプリッタ５０１で反射し、光学系５０３を介して被検体５２０に入射する。ＳＬＭ５０２は、光学系５０３の瞳面に配置されている。

被検体５２０に入射した光５１１は、散乱されながら被検体５２０の内部を伝搬し、その一部のエネルギーは被検体５２０の内部のある位置としての局所領域５１２の吸収体で吸収される。これにより、局所領域５１２の温度が上昇し、その局所領域５１２の体積が膨張して音響波（光音響信号）５１３が発生する。超音波トランスデューサを含む超音波装置（測定手段）５０４は、この光音響信号５１３を測定する。このとき、制御部５０５は、超音波装置５０４を制御して、被検体５２０内の局所領域５１２からの光音響信号５１３を含む信号が検出されるように、超音波トランスデューサのフォーカスを制御する。

超音波トランスデューサは、例えばリニアアレイ探触子により構成され、アレイ探触子を用いた電子フォーカスによって被検体５２０の内部の任意の位置に超音波フォーカス領域を生成することが可能である。トランスデューサとしては、圧電現象を用いたトランスデューサ、光の共振を用いたトランスデューサ、容量変化を用いたトランスデューサ等を用いることができる。また、超音波トランスデューサは、被検体５２０と音響的に整合されている。

被検体５２０の内部における光５１１の入射位置からの深さがｚの位置（ｘ，ｙ，ｚ）における光音響信号Ｐ(ｘ，ｙ，ｚ)は、以下の式（１１）ように表わされる。式（１１）において、Ｉ（ｘ，ｙ，ｚ）は位置（ｘ，ｙ，ｚ）における光強度であり、μ_ａは該位置（ｘ，ｙ，ｚ）にある吸収体の吸収係数である。また、Γは熱から音響波への変換効率を表すグリュナイゼン係数である。

式（８）から分かるように、光音響信号Ｐ(ｘ，ｙ，ｚ)は、位置（ｘ，ｙ，ｚ）における光強度Ｉ（ｘ，ｙ，ｚ）に比例する。したがって、位置（ｘ，ｙ，ｚ）での光強度の変調に応じて光音響信号が変調する。このことを利用して、非特許文献２では、入力をアダマール基底で変調された光とし、出力を光音響信号として応答特性行列が測定できることが開示されている。非特許文献２では、ある基底の入射光を被検体に照射し、照射したタイミングから超音波装置５０４で測定する時間差（遅延時間）を測ることで、超音波装置５０４の伝搬軸に沿った複数（Ｍ個）の位置の光音響信号を同時に測定する。つまり、１回の光照射でＭ個の位置（画素）の光音響信号の応答を測定する。これをＮ個の基底を用いて、順次、光音響信号を測定すれば、Ｎ×Ｍの応答特性行列を算出することができる。この光音響信号に基づいた応答特性行列の算出についても、位相シフト法を利用する。

測定される光音響信号の変調は、光音響信号発生位置(ｘ，ｙ，ｚ)での光の強度分布に比例する。したがって、位相シフト法を用いて光音響信号の強度の変調を計測することで、位置（ｘ，ｙ，ｚ）での光の位相分布を回復する。光音響信号の強度は、非特許文献２で開示されているように、例えば音圧信号の最大値と最小値の差（peak-to-peak）を測定すればよい。

ＳＬＭ５０２は、図３を用いて説明したのと同様にして、入射した光を位相変調する変調領域（物体光）と位相変調しない非変調領域（参照光）とに分割し、物体光に参照光に対する位相差（位相シフト）を与えることで位相シフト法を行う。

制御部５０５は、実施例１で図１を用いて説明した処理と同様の処理により応答特性行列を生成する。すなわち、制御部５０５は、Ｓ３１０において、位相ステップ量αをπ／２とし、測定回数ｍを５として超音波装置５０４に光音響信号の測定を行わせ、Ｐ_１（ｘ，ｙ，ｚ）〜Ｐ_５（ｘ，ｙ，ｚ）の測定データをそれぞれ異なる遅延時間に相当する複数の位置で得る。そして、Ｓ３２０において、各基底についての測定により得られた５つの測定データを用いて応答特性行列ｔを生成する。

Ｓ３２０の詳細は図４と同様である。すなわち、制御部５０５は、Ｓ３２２で非特許文献３に開示されている方法を用いて位相分布φ（ｘ，ｙ，ｚ）を回復し、Ｓ３２３で回復した位相分布と、バイアス成分および変調成分の結果と、Ｓ３１０での測定結果とを比較して回復誤差を算出する。この処理を全ての基底に対して行った後、Ｓ３２４で応答特性行列ｔと誤差マップ行列ｔ_εを生成する。

次にＳ３２５において、制御部５０５は、誤差マップ行列ｔ_εの全要素に対して回復誤差εが大きい順に補正対象要素を並び替えて誤差補正リストを作成する。さらに、Ｓ３２６では、制御部５０５は、補正の終了条件として、実施例１のように閾値以上の要素ではなく、Ｓ３２９における評価値が有意に変化するまでとする。例えば、評価値の変化が、測定ノイズレベル以下であればＳ３２６で補正の終了条件を満足すると判定する。Ｓ３２６にて評価値の変化が測定ノイズレベルより大きければ処理を反復する。測定ノイズレベルは予め測定しておくことが望ましい。

次にＳ３２７では、制御部５０５は、実施例１と同様に、補正する位相の位相探索範囲内を確認した上で、Ｓ３２８において補正すべき応答特性行列ｔの要素（位相）を補正する。そして、Ｓ３０９において、制御部５０５は、補正後の応答特性行列ｔを用いて、被検体５２０内部の局所領域５１２に光がフォーカスされるように、ＳＬＭ５０２に設定すべき入力波面（位相共役波面）を式（６）を用いて算出してＳＬＭ５０２に設定する。このとき、局所領域５１２は、Ｓ３１０で測定した結果を考慮して、測定された光音響信号の強度が高い位置を選択するとよい。

このようにして成形された入力波面を有する光（位相共役波）５１１は被検体５２０に照射され、被検体５２０の内部において散乱されながら局所領域５１２にフォーカスされる。式（８）に示したように、光音響信号Ｐ（ｘ，ｙ，ｚ）は、位置（ｘ，ｙ，ｚ）での光強度Ｉ（ｘ，ｙ，ｚ）に比例する。したがって、位相共役波の被検体５２０への照射により発生した光音響信号Ｐ（ｘ，ｙ，ｚ）は、ＳＬＭ５０２に位相共役波面以外の位相（例えばフラットな位相分布）を設定したときの光照射により発生した光音響信号（バックグランド）と比べて信号強度が向上する。この信号強度の向上（エンハンスメント）をηとしてＳ３２９で評価する。

本実施例では、実施例１における透過光の代わりに光音響信号で応答特性行列を構成し、その要素の補正も光音響信号のエンハンスメントηを用いて評価する。実施形態と本実施形態とでは、測定する物理量は異なっているが、エンハンスメントηで、応答特性行列の精度を評価できる点ではどちらも同じである。これを利用して、上述の終了条件を満たすまでＳ３２６以降の反復処理を実行することで、応答特性行列ｔを補正する。

Ｓ３２０で応答特性行列を生成した後は、Ｓ３３０で任意の位置で光音響信号の強度が向上（光がフォーカス）するように位相分布を式（６）に基づいて算出し、これをＳＬＭ５０２に設定する。

続いてＳ３４０において、制御部５０５は、ＳＬＭ５０２に設定された位相分布を有する光を被検体５２０に照射して、Ｓ３５０で光音響信号を測定する。任意の位置でのフォーカススポットの作成を繰り返し、該スポットをスキャンしながら、測定した光音響信号の３次元分布を画像化してモニタ５０６に表示してもよい。被検体５２０の内部の任意の位置に効率的に光をフォーカスさせて光音響信号を測定することができるので、高ＳＮＲに光音響信号を測定して画像化することができる。

また、制御部５０５は、Ｓ３３０において、非特許文献２で開示されているように、応答特性行列を特異値分解する。このうち上位の特異値は、被検体５２０内の局所位置の吸収物体から発生した光音響信号に対応する。これを利用して、入力の特異ベクトルを上位から順次ＳＬＭ５０２に設定して光を被検体５２０に照射することで、各吸収物体を選択的に照射して光音響信号を測定することができる。

また、任意の複数の波長を用いて上述した処理を行い、被検体５２０の内部の吸収分光特性を測定し、酸化ヘモグロビン、還元ヘモグロビンおよび水等の成分比率や、酸素飽和度等の代謝情報を求め、これを３次元画像としてイメージングすることも可能である。

さらに、光音響信号を用いた測定により生成した応答特性行列を、他のイメージングに利用することも可能である。例えば、得られた応答特性行列を利用して、被検体の内部の任意の位置に光をフォーカスさせ、該位置から発生する蛍光、透過光または反射光等の光信号を、別途用意した光検出器で測定して可視化することもできる。

本発明の実施例３である応答特性取得装置および応答特性取得方法について説明する。図９は、本実施例の応答特性取得装置であって、マルチモードファイバを利用したイメージング装置の構成を示している。

光源６００は可視域から近赤外域の波長を有するＣＷ光を発するレーザ光源である。光源６００からは平行光が発せられ、ＳＬＭ６０１の有効領域に応じてそのビーム径が調整される。ＳＬＭ６０１は透過型の位相変調デバイスであり、図３にて説明したように、光源６００からＳＬＭ６０１に入射するビームの一部を物体光とし、残りを参照光とする。また、実施例１で説明したように、物体光の位相を参照光に対してシフトさせることで位相シフト法を行う。逆に物体光の位相を固定して、参照光の位相をシフトさせてもよい。ＳＬＭ６０１は、ＰＣ等により構成される制御部（処理手段）６１２によって制御される。

また、ＳＬＭ６０１において物体光を生成する領域である位相変調領域は、応答特性行列を生成するための基底が表示される。ＳＬＭ６０１を透過した物体光および参照光６１０は、ビームスプリッタ６０６および光学系６０２を透過してマルチモードファイバ６０３にカップリングされる。

マルチモードファイバ６０３は、例えばコア径Ｄが５０〜１０００μｍで、長さが１〜２ｍ、ＮＡ（開口数）が０．１〜０．５のマルチモードファイバである。マルチモードファイバ６０３の内部では、λを入射光６１０の波長とするとき、ＮＡ×Ｄ／λに比例した数の伝搬モードで光が伝搬する。マルチモードファイバ６０３から出射した光６１１は、ＮＡ×Ｄ／λに比例した数の伝搬モードが重ね合わされた光である。制御部６１２は、マルチモードファイバ６０３の出射端の光強度分布を、光学系６０４を介してアレイセンサ６０５に測定させる。

このように構成された装置において、マルチモードファイバ６０３に対して、実施例１および実施例２で説明した方法により、応答特性行列を生成する。応答特性行列は、図１に示したＳ３１０および図１と図４に示したＳ３２０の処理によって生成される。これにより、より短時間で高精度にマルチモードファイバ６０３の応答特性行列を生成することが可能である。

応答特性行列を生成した制御部６１２は、次のステップ（図１に示すＳ３３０以降）でこの応答特性行列を利用する。これにより、本実施例のイメージング装置を、特許文献１にて開示されているように、マルチモードファイバ６０３を用いた内視鏡装置として使用することができる。医療用途では人体・動物等の生体を被検体として、工業用途では検査対象物を被検体として、マルチモードファイバ６０３の出射端を被検体の内部に挿入する。

図１０は、本実施例の装置を用いて被検体の内部の観察対象６２０の位置までマルチモードファイバ６０３を挿入した状態を示している。マルチモードファイバ６０３の出射端のある１点に光がフォーカスされるような入射光の位相分布を応答特性行列から算出し、ＳＬＭ６０１に設定する（Ｓ３３０）。次に、その位相分布において入射光６１０をマルチモードファイバ６０３に入射させ、再生されるフォーカス光を観察対象６２０に照射する（Ｓ３４０）。そして、観察対象６２０からの反射光をマルチモードファイバ６０３を介してビームスプリッタ６０６に導き、ビームスプリッタ６０６で反射した光を光学系６０７を介してアレイセンサ６０８に導く。制御部６１２は、アレイセンサ６０８に該光を測定させる（Ｓ３５０）。アレイセンサ６０８は、図９に示したアレイセンサ６０５と同じものであってもよい。

特許文献１にて開示されているように、マルチモードファイバ６０３の出射端のフォーカス位置のスキャンを行い、順次、観察対象６２０からの反射光を測定することで、被検体６２０をイメージングすることが可能である。さらに、被検体からの反射光以外にも、照射光によって励起される蛍光信号を用いてイメージングしてもよい。この場合、必要に応じて、光源６００からの光として、ＣＷ光ではなく、パルス光を用いてもよい。また、非特許文献４にて開示されているように、任意の入射光で被検体を照射し、アレイセンサ６０８を通じて測定した結果に対して、応答特性行列の逆行列（または疑似逆行列）を適用して被検体の像を回復することでイメージングを行ってもよい。

このように、実施例１，２で説明した応答特性取得方法は、マルチモードファイバの応答特性行列の生成にも適用可能であり、生成した応答特性行列を用いて１本のマルチモードファイバを用いたイメージング装置としての内視鏡装置を実現することができる。

以上説明したように、上記各実施例は、特に測定対象の被検体が散乱体であり、位相シフト法の測定において大きな誤差が生じる場合でも、比較的少ない測定回数で高精度に応答特性行列を生成することができる。これは応答特性行列の精度を定量的に評価できる評価値を利用することが可能であるからであり、その評価値が向上するように応答特性行列を補正するためである。（その他の実施例）
本発明は、上述の実施形態の１以上の機能を実現するプログラムを、ネットワーク又は記憶媒体を介してシステム又は装置に供給し、そのシステム又は装置のコンピュータにおける１つ以上のプロセッサーがプログラムを読出し実行する処理でも実現可能である。また、１以上の機能を実現する回路（例えば、ＡＳＩＣ）によっても実現可能である。

以上説明した各実施例は代表的な例にすぎず、本発明の実施に際しては、各実施例に対して種々の変形や変更が可能である。

１００ 光源
１７０ 空間光変調器
１８０ レンズ
２３０ 撮像デバイス（ＣＣＤ）
２４０ パーソナルコンピュータ（ＰＣ）
２００ 被検体

Claims (13)

1. 光源からの光の一部の位相を変調することで互いに位相差を有する第１の光と第２の光を生成する位相変調手段と、
   前記第１の光と前記被検体に照射された前記第２の光との干渉により生じる信号の測定を行う測定手段と、
   位相シフト法により前記信号の位相分布を回復し、前記第１および第２の光の前記被検体への入射条件が異なる前記信号から回復した複数の前記位相分布から応答特性行列を生成する処理手段とを有し、
   前記処理手段は、
   前記応答特性行列の複数の要素に対して前記位相の回復誤差を算出し、
   該回復誤差に基づいて少なくとも１つの前記要素に対する前記位相の補正を行うことを特徴とする応答特性取得装置。
2. 前記測定手段は、前記応答特性行列を用いて算出された位相共役波面を前記被検体に照射して前記信号を測定し、
   前記処理手段は、前記信号の測定結果から評価値を求め、該評価値に基づいて前記補正を行うことを特徴とする請求項１に記載の応答特性取得装置。
3. 前記評価値は、前記信号の最大強度であること、または前記位相共役波面ではない光の前記被検体への照射により測定される前記信号の強度に対する前記最大強度の比であることを特徴とする請求項２に記載の応答特性取得装置。
4. 前記処理手段は、前記回復誤差が大きい前記要素から順に前記補正を行うことを特徴とする請求項１から３のいずれか一項に記載の応答特性取得装置。
5. 前記処理手段は、所定値より大きい前記回復誤差を有する前記要素に対して前記補正を行うことを特徴とする請求項１から３のいずれか一項に記載の応答特性取得装置。
6. 前記処理手段は、前記回復誤差の変化が前記測定におけるノイズレベル以下になるまで前記補正を行うことを特徴とする請求項１から５のいずれか一項に記載の応答特性取得装置。
7. 前記処理手段は、前記回復誤差を、前記測定により得られた前記信号と前記位相シフト法により回復された前記位相に基づいて算出された前記信号の推定値との差分を用いて算出することを特徴とする請求項１から６のいずれか一項に記載の応答特性取得装置。
8. 前記処理手段は、前記回復誤差を、前記差分に対して前記推定値の変調度で重み付けして算出することを特徴とする請求項７に記載の応答特性取得装置。
9. 前記信号は、前記被検体からの透過光、反射光および光音響信号のいずれかであることを特徴とする請求項１から８のいずれか一項に記載の応答特性取得装置。
10. 前記信号は、前記透過光または前記反射光であり、
    前記測定手段は、前記透過光または前記反射光を測定するためのアレイセンサと、該アレイセンサに前記透過光または前記反射光を導く光学系とを含むことを特徴とする請求項９に記載の応答特性取得装置。
11. 前記信号は、前記光音響信号であり、
    前記測定手段は、前記光音響信号を測定するための超音波トランスデューサを含むことを特徴とする請求項９に記載の応答特性取得装置。
12. 光源からの光の一部の位相を変調することで互いに位相差を有する第１の光と第２の光を生成するステップと、
    前記第１の光と前記被検体に照射された前記第２の光との干渉により生じる信号の測定を行うステップと、
    位相シフト法により前記信号の位相分布を回復し、前記第１および第２の光の前記被検体への入射条件が異なる前記信号から回復した複数の前記位相分布から応答特性行列を生成する処理ステップとを有し、
    前記処理ステップは、
    前記応答特性行列の複数の要素としての前記位相の回復誤差を算出し、
    該回復誤差に基づいて少なくとも１つの前記要素としての前記位相の補正を行うことを特徴とする応答特性取得方法。
13. 光源からの光の一部の位相を変調することで互いに位相差を有する第１の光と第２の光を生成する位相変調手段と、前記第１の光と前記被検体に照射された前記第２の光との干渉により生じる信号の測定を行う測定手段とともに用いられるコンピュータに処理を実行させるコンピュータプログラムであって、
    前記処理は、
    位相シフト法により前記信号の位相分布を回復し、前記第１および第２の光の前記被検体への入射条件が異なる前記信号から回復した複数の前記位相分布から応答特性行列を生成する処理と、
    前記応答特性行列の複数の要素としての前記位相の回復誤差を算出し、該回復誤差に基づいて少なくとも１つの前記要素としての前記位相の補正を行う処理とを含むことを特徴とする応答特性取得プログラム。