(19)【発行国】日本国特許庁(JP)
(12)【公報種別】特許公報(B2)
(11)【特許番号】
(24)【登録日】2022-10-07
(45)【発行日】2022-10-18
(54)【発明の名称】組織の粘弾性を断層可視化する装置および方法
(51)【国際特許分類】
G01N 21/17 20060101AFI20221011BHJP
【ＦＩ】
G01N21/17 630
【請求項の数】
7
(21)【出願番号】P 2019518761
(86)(22)【出願日】2018-05-14
(86)【国際出願番号】 JP2018018455
(87)【国際公開番号】W WO2018212115
(87)【国際公開日】2018-11-22
【審査請求日】2021-05-13
(31)【優先権主張番号】P 2017096134
(32)【優先日】2017-05-15
(33)【優先権主張国・地域又は機関】JP
(73)【特許権者】
【識別番号】599002043
【氏名又は名称】学校法人 名城大学
(74)【代理人】
【識別番号】110002273
【氏名又は名称】特許業務法人インターブレイン
(72)【発明者】
【氏名】佐伯 壮一
【審査官】横尾 雅一
(56)【参考文献】
【文献】国際公開第２０１６／０３１６９７（ＷＯ，Ａ１）
【文献】特開２００８－１７０３６３（ＪＰ，Ａ）
【文献】米国特許出願公開第２０１５／０１４８６５４（ＵＳ，Ａ１）
【文献】QU Y Q , et al.，Acoustic Radiation Force Optical Coherence Elastography of Corneal Tissue，IEEE JOURNAL OF SELECTED TOPICS IN QUANTUM ELECTRONICS，2016年02月08日，VOL. 22, NO. 3，6803507，Digital Object Identifier 10.1109/JSTQE.2016.2524618
【文献】ZHU J , et al.，Imaging and characterizing shear wave and shear modulus under orthogonal acoustic radiation force excitation using OCT Doppler variance method，OPTICS LETTERS，2015年04月30日，Vol. 40, No. 9，p.2099-2102
(58)【調査した分野】(Int.Cl.，ＤＢ名)
Ｇ０１Ｎ ２１／００ － Ｇ０１Ｎ ２１／６１
Ａ６１Ｂ １／００
Ａ６１Ｂ ８／００ － Ａ６１Ｂ ８／１５
ＪＳＴＰｌｕｓ／ＪＭＥＤＰｌｕｓ／ＪＳＴ７５８０（ＪＤｒｅａｍＩＩＩ）
(57)【特許請求の範囲】
【請求項1】
光コヒーレンストモグラフィーを用いる光学系を含み、測定対象における組織の粘弾性を断層可視化する粘弾性可視化装置であって、前記組織に対して変形エネルギーを付与するための負荷装置と、
前記測定対象を経由するオブジェクトアームに設けられ、光源からの光を前記測定対象に導いて走査させる第１光学機構と、
前記測定対象を経由しないリファレンスアームに設けられ、前記光源からの光を参照鏡に導いて反射させる第２光学機構と、
前記測定対象にて反射した物体光と前記参照鏡にて反射した参照光とが重畳された干渉光を検出する光検出装置と、
前記負荷装置および各光学機構の駆動を制御し、それらの駆動に基づいて前記光検出装置から入力された光干渉信号を処理することにより、前記組織の粘弾性の断層分布を演算する制御演算部と、
前記組織の粘弾性を断層可視化する態様で表示する表示装置と、
を備え、
前記負荷装置が、前記測定対象の内部に設定された焦点に向けて超音波を出力することにより、前記組織に対して音響放射圧を負荷する非接触負荷装置であり、
前記制御演算部は、前記光検出装置から入力された光干渉信号に対して周波数解析を実行し、解析から得られるドップラー角周波数シフト量に基づいて前記組織の変形速度を算出し、その変形速度に基づいて得られるひずみ情報と、前記負荷装置による音響放射圧情報とに基づいて前記組織の粘弾性を演算することを特徴とする粘弾性可視化装置。
【請求項2】
前記測定対象に向けて照射される光の軸が、前記測定対象における前記音響放射圧の焦点を通るように設定されることを特徴とする請求項１に記載の粘弾性可視化装置。
【請求項3】
前記制御演算部は、前記組織の変形によって生じる所定の力学特徴量、又は前記組織の変形に伴う所定の力学特徴量の変化を、測定対象の断層位置に対応づけて演算し、その演算結果に基づいて前記組織の粘弾性の断層分布を演算することを特徴とする請求項１または２に記載の粘弾性可視化装置。
【請求項4】
前記光学機構は、前記測定対象の表面側から光を照射し、前記負荷装置は、前記測定対象の裏面側から超音波を出力することを特徴とする請求項１～３のいずれかに記載の粘弾性可視化装置。
【請求項5】
前記測定対象として再生組織を収容するための透光性および音波透過性を有する容器を支持する支持台を備え、前記支持台は、前記容器の裏面側を露出させるための孔を有することを特徴とする請求項４に記載の粘弾性可視化装置。
【請求項6】
前記制御演算部は、前記音響放射圧による加振力の位相と、前記ひずみの位相との位相差に基づき、前記組織の粘弾性を演算することを特徴とする請求項１～５のいずれかに記載の粘弾性可視化装置。
【請求項7】
測定対象の組織の粘弾性を断層可視化する粘弾性可視化方法であって、前記組織に対して超音波の音響放射圧にて荷重を負荷する負荷工程と、
光コヒーレンストモグラフィーを用いることにより、前記荷重の負荷により変形する前記組織からの後方散乱光に基づく光干渉信号を取得する干渉信号取得工程と、
前記光干渉信号を処理し、前記組織の変形に伴う所定の力学特徴量の変化に基づいて前記組織の粘弾性の断層分布を演算する演算工程と、
前記組織の粘弾性を断層可視化する表示工程と、
を備え、
前記演算工程は、前記光干渉信号に対して周波数解析を実行し、解析から得られるドップラー角周波数シフト量に基づいて前記組織の変形速度を算出し、その変形速度に基づいて得られるひずみ情報と、前記音響放射圧の情報とに基づいて前記組織の粘弾性を演算することを特徴とする粘弾性可視化方法。
【発明の詳細な説明】
【技術分野】
【0001】
本発明は、測定対象の組織の粘弾性をマイクロスケールにて断層表示可能な装置および方法に関する。
【背景技術】
【0002】
近年、医療診断技術の更なる発展に向けて光干渉断層画像法（Optical Coherence Tomography：以下「ＯＣＴ」という）の臨床応用が進められている。ＯＣＴは、低コヒーレンス光干渉を利用した断層画像法であり、マイクロスケールの高空間分解能にて生体組織の形態分布を可視化できる。また、２次元ＯＣＴの取得レートはビデオレート以上であり、高時間分解能も有している。
【0003】
そこで、発明者らは、ＯＣＴを用いて生体組織の力学特性を断層可視化する手法を提案している（特許文献１）。軟骨や皮膚等の生体組織では基質量の変化よりも組織の力学特性の変化が有意に現れるため、組織内部のマイクロバイオメカニクスの定量評価が診断技術の向上に寄与すると考えたものである。本手法は、生体組織の特性（粘弾性）を評価可能なＯＣＴシステムの開発を目指したものである。
【先行技術文献】
【特許文献】
【0004】
【文献】国際公開第２０１６／０３１６９７号
【発明の概要】
【発明が解決しようとする課題】
【0005】
ところで、このような力学特性評価を行う場合、ＯＣＴ計測に際して生体組織を変形させるアルゴリズムが組み込まれる。そのため、特許文献１の技術では、圧電素子を測定対象の表面に接触させて負荷する方法を採用している。しかし、例えば再生組織など、厳格な品質保証を要する測定対象を扱う場合、このような圧電素子の接触が組織の汚染を招き、システムの実用化を妨げるおそれがある。発明者は、このような知見のもと、上記ＯＣＴシステムの改善が必要であるとの認識に到った。
【0006】
本発明はこのような課題に鑑みてなされたものであり、その目的の一つは、厳格な品質保証を要する測定対象に対しても、ＯＣＴ計測による力学特性評価を可能とする技術を提供することにある。
【課題を解決するための手段】
【0007】
本発明のある態様は、ＯＣＴを用いる光学系を含み、測定対象における組織の粘弾性を断層可視化する粘弾性可視化装置に関する。この装置は、光源からの光を組織に導いて走査させるための光学機構と、組織に対して変形エネルギーを付与するための負荷装置と、負荷装置および光学機構の駆動を制御し、それらの駆動に基づいて光学系による光干渉信号を処理することにより、組織の粘弾性の断層分布を演算する制御演算部と、組織の粘弾性を断層可視化する態様で表示する表示装置と、を備える。負荷装置は、測定対象の内部に設定された焦点に向けて超音波を出力することにより、組織に対して音響放射圧による加振力を負荷する非接触負荷装置である。
【0008】
本発明の別の態様は、測定対象の組織の粘弾性を断層可視化する粘弾性可視化方法に関する。この方法は、組織に対して超音波の音響放射圧にて荷重を負荷する負荷工程と、ＯＣＴを用いることにより、荷重の負荷により変形する組織からの後方散乱光に基づく光干渉信号を取得する干渉信号取得工程と、光干渉信号を処理し、組織の変形に伴う所定の力学特徴量の変化に基づいて組織の粘弾性の断層分布を演算する演算工程と、組織の粘弾性を断層可視化する表示工程と、を備える。
【発明の効果】
【0009】
本発明によれば、品質保証を要する測定対象に対しても、ＯＣＴ計測による力学特性評価を可能とする技術を提供できる。
【図面の簡単な説明】
【0010】
【図1】実施例に係る粘弾性可視化装置の構成を概略的に表す図である。
【図2】対象に非接触にて荷重負荷を行う方法を表す説明図である。
【図3】再帰的相互相関法による処理手順を概略的に示す図である。
【図4】サブピクセル解析による処理手順を概略的に示す図である。
【図5】制御演算部により実行される粘弾性断層可視化処理の流れを示すフローチャートである。
【図9】生体組織を模した試料についての実験結果を表す図である。
【図10】生体組織を模した試料についての実験結果を表す図である。
【図11】剛性を層状に変化させた層状試料についての実験結果を表す図である。
【図12】剛性を層状に変化させた層状試料についての実験結果を表す図である。
【図13】粘性評価に関する実験結果を表す図である。
【図14】変形例に係る荷重負荷方法を表す図である。
【発明を実施するための形態】
【0011】
本発明の一実施形態は、ＯＣＴを用いる粘弾性可視化装置である。この装置は、光源、光学機構、負荷装置、制御演算部および表示装置を備え、測定対象における組織の粘弾性を断層可視化する。「測定対象」は、皮膚や軟骨等の生体組織であってもよいし、再生皮膚や再生軟骨等の再生組織（培養細胞からなる組織）であってもよい。負荷装置は、測定対象の内部に設定される焦点に向けて超音波を出力し、組織に対して音響放射圧を負荷する。焦点は、測定対象の内部において音響放射圧による加振力が集中するポイントであり、制御演算部により設定される。
【0012】
この態様では、組織の粘弾性を断層可視化するために、音響放射圧による負荷とＯＣＴによる検出が行われる。音響放射圧による負荷は、動的荷重（加振力）であってよい。音響放射圧により測定対象の所望の断層位置に非接触にて荷重を付与し、組織変形を誘発できる。また、ＯＣＴにより非接触・非侵襲にて断層計測が可能となる。すなわち、この態様によれば、負荷装置および光学機構の双方を測定対象に対して非接触とすることができる。このため、測定対象が厳格な品質保証を要するものであったとしても、その汚染を防止でき、力学特性（粘弾性）を評価できる。
【0013】
測定対象に向けて照射される光の軸が、測定対象における音響放射圧の焦点を通るように設定されてもよい。それにより、その焦点位置にある組織の変形をダイレクトかつリアルタイムに検出できる。あるいは、その光軸を焦点から所定距離の位置に設定し、その距離に応じた粘弾性の評価基準を設定してもよい。
【0014】
制御演算部は、組織の変形によって生じる所定の力学特徴量、又は組織の変形に伴う力学特徴量の変化を測定対象の断層位置に対応づけて演算し、その演算結果に基づいて組織の粘弾性の断層分布を演算してもよい。ここでいう「力学特徴量」は、組織の変形量や変形ベクトルに基づいて得られるものでよい。例えば、変形ベクトルを空間微分したひずみテンソルであってもよい。ひずみの振幅や位相であってもよい。あるいは、変形ベクトルを時間微分した変形速度ベクトルであってもよいし、変形速度ベクトルをさらに空間微分したひずみ速度テンソルであってもよい。それらの振幅や位相であってもよい。音響放射圧による加振力の振幅と、ＯＣＴにより検出されるひずみの振幅と、それらの位相差とに基づき、粘弾性を算出してもよい。音響放射圧による加振力が負荷装置の駆動により生成されるものであるため、その加振力の振幅および位相が既知であることを利用できる。
【0015】
負荷変動に対する、変形速度（ひずみ速度）の変動の大きさ（振幅値）と、加振力の変動に対する変形速度（ひずみ速度）の変動の追従性（応答性）が、組織の粘弾性に対応して変化する傾向がある。具体的には、粘性が低下するほど、変形速度（ひずみ速度）の振幅値は大きくなり、時間遅れ（位相遅れ）は小さくなる傾向にある。このような傾向に基づいて粘弾性を評価することもできる。
【0016】
あるいは、測定対象への超音波の印可によって発生する剪断波を断層検出し、その伝播速度から組織の粘弾性を演算してもよい。音響放射圧による非接触エラストグラフィーを利用し、組織の粘弾性をＯＣＴにより検出するものである。
【0017】
より具体的には、光学機構が測定対象の表面側から光を照射し、負荷装置が測定対象の裏面側から超音波を出力する構成としてもよい。それにより、光学機構と負荷装置との物理的な干渉を防止し易くなり、装置の設計が容易となる。
【0018】
このような構成において再生組織を測定対象とする場合、その再生組織を収容する容器を支持する支持台を設けてもよい。この容器として透光性および音波透過性を有するものを採用し、支持台にはその容器の裏面側を露出させるための孔を設けてもよい。それにより、測定対象の表面側から光を受け入れ、裏面側から超音波を受け入れることができる。再生組織の品質保証のため、この容器は密閉可能であることが好ましい。支持台をクリーンベンチ内に配置してもよい。
【0019】
この装置は、第１光学機構、第２光学機構および光検出装置を備えてもよい。第１光学機構は、測定対象を経由するオブジェクトアームに設けられ、光源からの光を測定対象に導いて走査させる。第２光学機構は、測定対象を経由しないリファレンスアームに設けられ、光源からの光を参照鏡に導いて反射させる。光検出装置は、測定対象にて反射した物体光と参照鏡にて反射した参照光とが重畳された干渉光を検出する。制御演算部は、光検出装置から入力された光干渉信号に対して周波数解析を実行し、ドップラー周波数を用いて組織の変形速度を算出する。そして、その変形速度に基づいて得られるひずみ情報と、前記負荷装置による音響放射圧情報とに基づいて粘弾性を演算してもよい。
【0020】
以下、図面を参照しつつ、本実施形態を具体化した実施例について詳細に説明する。
［実施例］
図１は、実施例に係る粘弾性可視化装置の構成を概略的に表す図である。本実施例の装置は、再生組織の粘弾性をマイクロスケールにて断層計測し、可視化する。この断層計測のために、音響放射圧による荷重負荷とＯＣＴによる検出を利用する。
［実施例］
図１は、実施例に係る粘弾性可視化装置の構成を概略的に表す図である。本実施例の装置は、再生組織の粘弾性をマイクロスケールにて断層計測し、可視化する。この断層計測のために、音響放射圧による荷重負荷とＯＣＴによる検出を利用する。
【0021】
図１に示すように、ＯＣＴ装置１は、光源２、オブジェクトアーム４、リファレンスアーム６、光学機構８，１０、光検出装置１２、負荷装置１３、制御演算部１４および表示装置１６を備える。各光学要素は、光ファイバにて互いに接続されている。なお、図示の例では、マッハツェンダー干渉計をベースとした光学系が示されているが、マイケルソン干渉計その他の光学系を採用することもできる。また、本実施例では、ＯＣＴとしてＴＤ－ＯＣＴ（Time Domain OCT）を用いるが、ＳＳ－ＯＣＴ（Swept Source OCT）、ＳＤ－ＯＣＴ（Spectral Domain OCT）その他のＯＣＴを用いてもよい。
【0022】
光源２から出射された光は、カプラー１８（ビームスプリッタ）にて分けられ、その一方がオブジェクトアーム４に導かれて物体光となり、他方がリファレンスアーム６に導かれて参照光となる。オブジェクトアーム４の物体光は、サーキュレータ２０を介して光学機構８に導かれ、測定対象（以下、「対象Ｓ」という）に照射される。対象Ｓは、本実施例では再生組織（培養組織）である。この物体光は、対象Ｓの表面および断面にて後方散乱光として反射されてサーキュレータ２０に戻り、カプラー２２に導かれる。
【0023】
一方、リファレンスアーム６の参照光は、サーキュレータ２４を介して光学機構１０に導かれる。この参照光は、光学機構１０の参照鏡２６にて反射されてサーキュレータ２４に戻り、カプラー２２に導かれる。すなわち、物体光と参照光とがカプラー２２にて合波（重畳）され、その干渉光が光検出装置１２により検出される。光検出装置１２は、これを光干渉信号（干渉光強度を示す信号）として検出する。この光干渉信号は、Ａ／Ｄ変換器２８を介して制御演算部１４に入力される。
【0024】
制御演算部１４は、光学系全体の制御と、負荷装置（超音波加振装置）の制御と、ＯＣＴによる画像出力のための演算処理を行う。制御演算部１４の指令信号は、図示略のＤ／Ａ変換器を介して光源２、光学機構８，１０、負荷装置１３等に入力される。制御演算部１４は、光検出装置１２に入力された光干渉信号を処理し、ＯＣＴによる対象Ｓの断層画像を取得する。そして、その断層画像データに基づき、後述の手法により対象Ｓにおける粘弾性の断層分布を演算する。
【0025】
より詳細には以下のとおりである。
光源２は、スーパールミネッセントダイオード（Super Luminessent Diode：以下「ＳＬＤ」と表記する）からなる広帯域光源である。光源２から出射された光は、カプラー１８にて物体光と参照光に分けられ、それぞれオブジェクトアーム４、リファレンスアーム６に導かれる。オブジェクトアーム４およびリファレンスアーム６には光ファイバとして、偏光による影響を抑えることが可能な偏波保持光ファイバ（Polarization Maintaining Fiber）が用いられている。
光源２は、スーパールミネッセントダイオード（Super Luminessent Diode：以下「ＳＬＤ」と表記する）からなる広帯域光源である。光源２から出射された光は、カプラー１８にて物体光と参照光に分けられ、それぞれオブジェクトアーム４、リファレンスアーム６に導かれる。オブジェクトアーム４およびリファレンスアーム６には光ファイバとして、偏光による影響を抑えることが可能な偏波保持光ファイバ（Polarization Maintaining Fiber）が用いられている。
【0026】
光学機構８は、オブジェクトアーム４を構成し、光源２からの光を対象Ｓに導いて走査させる機構と、その機構を駆動するための駆動部を備える。光学機構８は、コリメータレンズ４０、２軸のガルバノミラー４２、および対物レンズ４４を含む。対物レンズ４４は、対象Ｓに対向配置される。サーキュレータ２０を経た光は、コリメータレンズ４０を介してガルバノミラー４２に導かれ、ｘ軸方向やｙ軸方向に走査されて対象Ｓに照射される。対象Ｓからの反射光は、物体光としてサーキュレータ２０に戻り、カプラー２２に導かれる。
【0027】
光学機構１０は、ＲＳＯＤ方式（Rapid Scanning Optical Delay Line）の機構であり、リファレンスアーム６を構成する。このＲＳＯＤは、光が後述の回折格子５２に往復１回ずつ照射されるシングルパス型でもよいし、往復２回ずつ照射されるダブルパス型でもよい。光学機構１０は、コリメータレンズ５０、回折格子５２、レンズ５４および参照鏡２６を含む。サーキュレータ２４を経た光は、コリメータレンズ５０を介して回折格子５２に導かれる。この光は、回折格子５２によって波長ごとに分光され、それぞれレンズ５４によって参照鏡２６上の異なる位置に集光される。参照鏡２６を微小角にて回転させることで、高速光路走査および周波数変調が可能となる。参照鏡２６からの反射光は、参照光としてサーキュレータ２４に戻り、カプラー２２に導かれる。そして、物体光と重畳されて干渉光として光検出装置１２に送られる。なお、変形例においては、回折格子５２を経た光をレンズ５４に代えて湾曲ミラーによって参照鏡２６に集光してもよい。
【0028】
光検出装置１２は、光検出器６０、フィルタ６２およびアンプ６４を含む。カプラー２２を経ることで得られた干渉光は、光検出器６０にて光干渉信号として検出される。この光干渉信号は、フィルタ６２によりノイズが除去又は低減され、アンプ６４およびＡ／Ｄ変換器２８を経て制御演算部１４に入力される。
【0029】
負荷装置１３は、ＯＣＴの光軸上に対象Ｓを保持するための支持台７０、対象Ｓに向けて超音波を出力する超音波デバイス７２、および超音波デバイス７２を駆動する駆動回路７４を含む。超音波デバイス７２の駆動（発振）により、対象Ｓの内部に向けて超音波ビームが出力される。駆動回路７４は、制御演算部１４からの指令に基づき、ＯＣＴの検出と同期させるようにその超音波ビームを縦横に走査させる。それにより、対象Ｓの組織の任意位置に非接触にて音響放射圧による動的荷重（振動荷重、加振力）を負荷できる。この発振による音響放射圧負荷信号は、振幅変調させておく必要がある。ここでは、音響放射圧の音圧量（圧縮荷重量）に直接依存する振幅を正弦波で変調して送信する。この振幅変調にて組織の変形が促されるため、この変調による変形速度をＯＣＴにより検出する。
【0030】
制御演算部１４は、ＣＰＵ、ＲＯＭ、ＲＡＭ、ハードディスクなどを有する。制御演算部１４は、これらのハードウェアおよびソフトウェアによって、光学系全体の制御と、負荷装置１３の駆動制御と、ＯＣＴによる画像出力のための演算処理を行う。制御演算部１４は、負荷装置１３および光学機構８，１０の駆動を制御し、それらの駆動に基づいて光検出装置１２から出力された光干渉信号を処理し、ＯＣＴによる対象Ｓの断層画像を取得する。そして、その断層画像データに基づき、後述の手法により対象Ｓにおける組織の粘弾性を演算する。
【0031】
表示装置１６は、例えば液晶ディスプレイからなり、制御演算部１４にて演算された対象Ｓの組織の粘弾性を断層可視化する態様で画面に表示する。
【0032】
図２は、対象Ｓに非接触にて荷重負荷を行う方法を表す説明図である。（Ａ）は、負荷装置１３およびその周辺の構成を模式的に示す。（Ｂ）は、音響放射圧の負荷状態を示す。
図２（Ａ）に示すように、対象Ｓは容器８０に収容され、その容器８０が支持台７０にセットされる。本実施形態では、対象Ｓが再生組織であるため、荷重負荷やＯＣＴ計測に際して汚染を生じさせないようにするものである。容器８０は、対象Ｓが播種された基材８２を培養液Ｃに浸した状態で収容し、蓋８４により閉止（密閉）されている。容器８０，基材８２および蓋８４は、透光性および音波透過性を有する材質からなる。なお、光学機構８および負荷装置１３の少なくとも一部をクリーンベンチ内に収容し、容器８０内への塵埃や雑菌の混入を確実に防止できるのが好ましい。
図２（Ａ）に示すように、対象Ｓは容器８０に収容され、その容器８０が支持台７０にセットされる。本実施形態では、対象Ｓが再生組織であるため、荷重負荷やＯＣＴ計測に際して汚染を生じさせないようにするものである。容器８０は、対象Ｓが播種された基材８２を培養液Ｃに浸した状態で収容し、蓋８４により閉止（密閉）されている。容器８０，基材８２および蓋８４は、透光性および音波透過性を有する材質からなる。なお、光学機構８および負荷装置１３の少なくとも一部をクリーンベンチ内に収容し、容器８０内への塵埃や雑菌の混入を確実に防止できるのが好ましい。
【0033】
支持台７０には、上下方向の貫通孔８６が設けられている。容器８０は、貫通孔８６を跨ぐように支持台７０に載置され、その貫通孔８６を介して裏面側を露出させる。対象Ｓは、少なくともその計測範囲が貫通孔８６の軸線方向の投影面に含まれるように配置される。超音波デバイス７２は、貫通孔８６を介して容器８０の裏面に対向し、対象Ｓに向けて超音波を出力する。
【0034】
超音波デバイス７２は、その上面に曲面状の振動子アレイ９０を有する。振動子アレイ９０は、曲面状に配置された複数の素子９２を有する。素子９２は、加振用の超音波を発生させる超音波素子（ピエゾ素子）からなる。図２（Ｂ）にも示すように、各素子９２からの超音波の出力を調整することにより、対象Ｓ内の所望位置を焦点Ｆとして音響放射圧を集束させ、その焦点Ｆにて組織を加振する（変位させる）ことができる（白抜き矢印参照）。この音響放射圧の生成によって、焦点Ｆを起点として対象Ｓの表面と平行な方向にせん断波が発生する（二点鎖線矢印参照）。本実施例では、ＯＣＴによる光の軸が、対象Ｓの焦点Ｆを通るように設定される。焦点Ｆの位置は、複数の素子９２のうち駆動対象（通電対象）を選択する（変化させる）ことで調整（走査）できる。なお、変形例においては、超音波素子を平面状等に配置し、音響レンズを用いることで音響放射圧を焦点Ｆに集束させてもよい。
【0035】
制御演算部１４は、対象Ｓにおける音響放射圧の焦点Ｆ、駆動する超音波素子、音響放射圧による加振周波数等を設定し、これらの設定に基づいて駆動回路７４へ制御信号（パルス）を出力する。それにより、駆動回路７４から超音波デバイス７２へ駆動信号（パルス）が出力され、焦点Ｆに向けて超音波（加振パルス）が出力される。
【0036】
以下、対象Ｓの粘弾性の演算処理方法について詳細に説明する。
上述のように、ＯＣＴにおいて、オブジェクトアーム４を経た物体光（測定対象からの反射光）と、リファレンスアーム６を経た参照光とが合波され、光検出装置１２により光干渉信号として検出される。制御演算部１４は、この光干渉信号を干渉光強度に基づく対象Ｓの断層画像として取得することができる。この断層分布は三次元で演算することもできるが、ここでは一次元又は二次元での演算について説明する。
上述のように、ＯＣＴにおいて、オブジェクトアーム４を経た物体光（測定対象からの反射光）と、リファレンスアーム６を経た参照光とが合波され、光検出装置１２により光干渉信号として検出される。制御演算部１４は、この光干渉信号を干渉光強度に基づく対象Ｓの断層画像として取得することができる。この断層分布は三次元で演算することもできるが、ここでは一次元又は二次元での演算について説明する。
【0037】
ＯＣＴの光軸方向（奥行方向）の分解能であるコヒーレンス長ｌｃは、光源の自己相関関数によって決定される。ここでは、コヒーレンス長ｌｃを自己相関関数の包括線の半値半幅とし、下記式（１）にて表すことができる。
ここで、λｃはビームの中心波長であり、Δλはビームの半値全幅である。
【数1】
【0038】
一方、光軸垂直方向（ビーム走査方向）の分解能は、集光レンズによる集光性能に基づき、ビームスポット径Ｄの１／２とされる。そのビームスポット径ΔΩは、下記式（２）にて表すことができる。
ここで、ｄは集光レンズに入射するビーム径、ｆは集光レンズの焦点距離である。
【数2】
【0039】
ＯＣＴによりドップラー変調信号を検出することにより、組織の変形速度分布を算出できる。すなわち、参照光の電場Ｅ'r(t)は下記式（３）にて表される。
【数3】
【0040】
ここで、Ａ(t)は振幅、ωｃは光源の中心角周波数である。ωｒはＲＳＯＤの参照鏡２６の回転により生じるドップラー角周波数である。
【0041】
また、組織の変形速度により生じるドップラー角周波数シフト量ωｄを考慮し、物体光の電場Ｅ'ｏ(t)は下記式（４）にて表される。
【数4】
【0042】
光強度は電場強度の二乗の時間平均であるため、検出される干渉光強度Ｉ'ｄ(t)は下記式（５）で表される。
【数5】
【0043】
また、検出される断層干渉信号Ｉｄ(x,y,z)は下記式（６）で表される。
【数6】
【0044】
一方、干渉信号の角周波数はωｒ－ωｄとなる。この干渉信号のキャリア角周波数ωｒを用いることにより、変形速度により生じたドップラー角周波数シフト量ωｄを検出する。そして、検出されたドップラー角周波数シフト量ωｄを用いて下記式（７）を演算することにより、変形速度ｖを得ることができる。
ここで、λｃは光源の中心波長、θ(x,y,z)は座標(x,y,z)における変形速度の方向とビームの入射方向とがなす角度である。ｎは組織内部の平均屈折率である。
【数7】
【0045】
（隣接自己相関法）
本実施例では上述のように、ＲＳＯＤを用いた奥行方向（Ｚ軸方向）走査手法を採用する。その際の変形速度検出能の劣化を防止するために、ヒルベルト変換および隣接自己相関法を適用する。すなわち、空間的に隣接する干渉信号にヒルベルト変換を適用して得られる解析信号（複素信号）Γ(t)に対し、隣接自己相関法を適用する。それにより、任意座標における所定時間間隔（時間差ΔＴ）での位相差Δφを求め、変形速度によるドップラー角周波数シフト量ωｄを検出する。
本実施例では上述のように、ＲＳＯＤを用いた奥行方向（Ｚ軸方向）走査手法を採用する。その際の変形速度検出能の劣化を防止するために、ヒルベルト変換および隣接自己相関法を適用する。すなわち、空間的に隣接する干渉信号にヒルベルト変換を適用して得られる解析信号（複素信号）Γ(t)に対し、隣接自己相関法を適用する。それにより、任意座標における所定時間間隔（時間差ΔＴ）での位相差Δφを求め、変形速度によるドップラー角周波数シフト量ωｄを検出する。
【0046】
すなわち、ヒルベルト変換を適用して得られたj，j+1番目の解析信号をそれぞれΓj，Γj+1とすると、それぞれの干渉信号は下記式（８）として表される。
ここで、ｓ(t)は解析信号Γ(t)の実部を表し、ｓ^(t)は虚部を表す。ΔＴはj,j+1番目の干渉信号の取得時間間隔を表し、Ａは干渉信号の振幅（つまり後方散乱強度）を表す。
【数8】
【0047】
物体光の電場において位相変調が生じていない場合（ωｄ＝０）、それぞれの干渉信号の位相差はゼロとなる。ΓjとΓj+1の干渉信号における位相差は、変形速度によって生じるドップラー変調による位相変化量ωｄΔＴに相当する。すなわち、変形速度によって生じるドップラー角周波数シフト量ωｄは下記式（９）として表される。
【数9】
【0048】
なお，偏角Δφ＝ωdΔＴであるため、位相の変化量は－π～πの範囲にて検出可能である。また、下記式（１０）に従い、ｎ本の干渉信号についてアンサンブル平均処理を施すことにより、ドップラー角周波数シフト量ωｄの検出能を向上させることができる。
【数10】
【0049】
以上のようにして得られたドップラー角周波数シフト量ωｄを上記式（７）に代入することにより、変形速度ｖを算出できる。この変形速度ｖに基づいて組織の粘弾性を演算できる。
【0050】
すなわち、粘弾性体の複素弾性率が下記式（１１）のＥ＊で表されることは周知である。複素弾性率Ｅ＊の実数部分Ｅ'は応力とひずみが同位相の成分に相当し、虚数部分Ｅ"は応力とひずみが９０度位相ずれている成分に相当する。前者のＥ'は弾性成分に相当し、「貯蔵弾性率」と呼ばれる（以下「貯蔵弾性率Ｅ'」という）。後者のＥ"は応力とひずみ速度が同位相となる粘性成分に相当し、「損失弾性率」と呼ばれる（以下「損失弾性率Ｅ"」という）。したがって、貯蔵弾性率Ｅ'および損失弾性率Ｅ"の分布を、粘弾性の分布として演算できる。
【0051】
ここで、組織の変形速度ｖの分布（Ｚ方向の分布）から最小二乗法により各時刻のひずみ量（ひずみ速度）を算出でき、そのひずみ量の振幅（ひずみ振幅ε0）は断層検出できる。また、超音波の音響放射圧の大きさ（応力振幅σ0）は、制御演算部１４が超音波デバイス７２に対して設定する制御量であるため、把握できる。一方、制御演算部１４にて超音波の発振タイミングとＯＣＴの計測タイミングとを把握しているため、応力とひずみの位相差（損失角δ）を算出可能である。したがって、粘弾性パラメータである貯蔵弾性率Ｅ'と損失弾性率Ｅ"を、それぞれ下記式（１１）より算出できる。
【数11】
【0052】
このようにして得られた貯蔵弾性率Ｅ'と損失弾性率Ｅ"の断層分布が、組織の粘弾性の断層分布を実質的に示すことになる。制御演算部１４は、その組織の粘弾性を断層可視化する。
【0053】
一方、上記処理とは異なり、二次元断層画像を予め取得したうえでのポスト処理にはなるが、組織の粘弾性を二次元又は三次元的に断層可視化することができる。この手法は、計測対象の変形前後の２枚のＯＣＴ断層画像にデジタル相互相関法を適用することで変形ベクトル分布を算出し、マイクロスケールにて組織のひずみ速度テンソル分布を断層検出する手法である。
【0054】
変形ベクトル分布を算出する際には、繰り返し相互相関処理を実施する再帰的相互相関法(Recursive Cross-correlation method)を適用する。これは、低解像度において算出された変形ベクトルを参照し、探査領域を限定するとともに階層的に検査領域を縮小して相互相関法を適用する手法である。これにより、高解像度な変形ベクトルを取得することができる。さらに、スペックル・ノイズ低減法として、隣接検査領域の相関値分布との乗算を行う隣接相互相関乗法(Adjacent Cross-correlation Multiplication)を用いる。そして、乗算されることによって高ＳＮ化した相関値分布から最大相関値を探索する。
【0055】
また、マイクロスケールにおける微小変形解析では、変形ベクトルのサブピクセル精度が重要となる。このため、輝度勾配を利用する風上勾配法(Up-stream Gradientmethod)と、伸縮および剪断を考慮した画像変形法(Image Deformation method)の両サブピクセル解析法を併用し、変形ベクトルの高精度検出を実現する。なお、ここでいう「風上勾配法」は、勾配法（オプティカルフロー法）の一種である。
【0056】
【0057】
（再帰的相互相関法）
図３（Ａ）～（Ｃ）は、再帰的相互相関法による処理過程を示している。各図にはＯＣＴにより連続的に撮影される前後の断層画像が示されている。左側には先の断層画像（Image1）が示され、右側には後の断層画像（Image2）が示されている。
図３（Ａ）～（Ｃ）は、再帰的相互相関法による処理過程を示している。各図にはＯＣＴにより連続的に撮影される前後の断層画像が示されている。左側には先の断層画像（Image1）が示され、右側には後の断層画像（Image2）が示されている。
【0058】
相互相関法とは、局所的なスペックル・パターンの類似度を下記式（１２）に基づく相関値Ｒｉｊより評価する方法である。そのため、図３（Ａ）に示すように、連続的に撮影された前後のＯＣＴ画像について、先の断層画像（Image1）に類似度の検査対象となる検査領域Ｓ１が設定され、後の断層画像（Image2）に類似度の探査範囲となる探査領域Ｓ２が設定される。
ここで、空間座標として、光軸方向にＺ軸、光軸と垂直方向にＸ軸を設定している。ｆ(Ｘｉ,Ｚｊ)とｇ(Ｘｉ,Ｚｊ)は、変形前後のＯＣＴ画像において設置された中心位置(Ｘｉ,Ｚｊ)の検査領域Ｓ１(Ｎｘ×Ｎｚピクセル)のスペックル・パターンを表している。
なお、ｆ－とｇ－はｆ(Ｘｉ,Ｚｊ)とｇ(Ｘｉ,Ｚｊ)の検査領域Ｓ１内の平均値を表している。
【数12】
なお、ｆ－とｇ－はｆ(Ｘｉ,Ｚｊ)とｇ(Ｘｉ,Ｚｊ)の検査領域Ｓ１内の平均値を表している。
【0059】
探査領域Ｓ２(Ｍｘ×Ｍｚピクセル)内における相関値分布Ｒi,j(ΔＸ,ΔＺ)を算出し、図３（Ｂ）に示すようにパターンマッチングを行う。実際には、下記式（１３）に示すように、最大相関値を与える移動量Ｕｉ,jを変形前後の変形ベクトルとして決定する。
【数13】
【0060】
本手法では、検査領域Ｓ１を縮小しながら相互相関処理を繰り返して空間解像度を高める再帰的相互相関法を採用している。なお、本実施例では解像度を上げる際には空間解像度が倍になるようにしている。図３（Ｃ）に示すように、検査領域Ｓ１を１／４に分割し、前階層にて算出された変形ベクトルを参照し、探査領域Ｓ２を縮小している。ここで探査領域Ｓ２も１／４に分割している。再帰的相互相関法を用いることで、高解像度において多発する過誤ベクトルの抑制を可能にしている。このような再帰的相互相関処理を施すことにより、変形ベクトルの解像度を高めることができる。
【0061】
また、これに加え、下記式（１４）により、演算中の座標を中心とする周囲８つの座標を含む合計９つの変形ベクトルの平均偏差σを用いた閾値を設定し、過誤ベクトルの除去を行い再帰処理に伴う誤差伝播を抑制する。
ここで、Ｕmはベクトル量の中央値を表し、閾値となる係数κは任意に設定した。
【数14】
【0062】
（隣接相互相関乗法）
本実施例では、スペックルノイズの影響を受けたランダム性の強い相関値分布から正確な最大相関値を決定する方法として、隣接相互相関乗法を導入している。下記式（１５）により、隣接相互相関乗法では検査領域Ｓ１における相関値分布Ｒi,j(Δx,Δz)と、その検査領域Ｓ１にオーバーラップする隣接検査領域に対するＲi＋Δi,j(Δx,Δz)とＲi,j＋Δj(Δx,Δz)の乗算を行い、新たな相関値分布Ｒ'i,j(Δx,Δz)を用いて最大相関値を検索する。
本実施例では、スペックルノイズの影響を受けたランダム性の強い相関値分布から正確な最大相関値を決定する方法として、隣接相互相関乗法を導入している。下記式（１５）により、隣接相互相関乗法では検査領域Ｓ１における相関値分布Ｒi,j(Δx,Δz)と、その検査領域Ｓ１にオーバーラップする隣接検査領域に対するＲi＋Δi,j(Δx,Δz)とＲi,j＋Δj(Δx,Δz)の乗算を行い、新たな相関値分布Ｒ'i,j(Δx,Δz)を用いて最大相関値を検索する。
【数15】
【0063】
これにより、相関値同士の乗算によってランダム性を低減させることが可能になる。上述した検査領域Ｓ１の縮小と共に干渉強度分布の情報量も減少するため、スペックル・ノイズを原因とする複数相関ピークの出現が計測精度の悪化を招いていると考えられる。一方、隣接境界同士の移動量には相関があるため、最大相関値座標付近では強い相関値が残存する。この隣接相互相関乗法の導入によって最大相関値ピークが明瞭化され計測精度が向上し、正確な移動座標を抽出することが可能となる。また、この隣接相互相関乗法をＯＣＴの各ステージに導入することで、誤差伝播が抑制され、スペックル・ノイズに対するロバスト性が向上する。それにより、高空間解像度においても高精度な変形ベクトルとひずみ分布の算出が可能となる。
【0064】
（風上勾配法）
図４（Ａ）～（Ｃ）は、サブピクセル解析による処理過程を示している。各図にはＯＣＴにより連続的に撮影される前後の断層画像が示されている。左側には先の断層画像（Image1）が示され、右側には後の断層画像（Image2）が示されている。
図４（Ａ）～（Ｃ）は、サブピクセル解析による処理過程を示している。各図にはＯＣＴにより連続的に撮影される前後の断層画像が示されている。左側には先の断層画像（Image1）が示され、右側には後の断層画像（Image2）が示されている。
【0065】
本実施例では、サブピクセル解析のために風上勾配法と画像変形法を採用する。最終的な移動量の算出は後述の画像変形法によるが、計算の収束性の問題から、画像変形法より先に風上勾配法を適用する。検査領域サイズが小さく高空間解像度の条件において、サブピクセル移動量を高精度検出する画像変形法及び風上勾配法を適用している。画像変形法におけるサブピクセル移動量の検出が困難な場合において、風上勾配法によりサブピクセル移動量を算出する。
【0066】
サブピクセル解析では、注目点における変形前後の輝度差が各成分の輝度勾配と移動量によって表される。このため、検査領域Ｓ１内の輝度勾配データより最小二乗法を用いてサブピクセル移動量を決定することができる。本実施例では、輝度勾配を求める際に、サブピクセル変形前の風上側の輝度勾配を与える風上差分法を採用している。
【0067】
サブピクセル解析は計測誤差からだけではなく、組織の散乱効果や偏光特性に複雑に関係したスペックル・ノイズに強く依存している。さらに、ひずみテンソル分布の算出には移動量の空間微係数が必要となるため、サブピクセルエラーが微係数算出の数値不安定を増大させ、ひずみテンソル分布の不安定化を招いてしまう。サブピクセル解析は様々な手法が存在するが、本実施例では検査領域サイズが小さく高空間解像度の条件においても、サブピクセル移動量を高精度検出する勾配法を採用している。
【0068】
風上勾配法は、検査領域Ｓ１内の注目点の移動を、図４（Ａ）に示すピクセル精度に留まらず、図４（Ｂ）に示すサブピクセル精度にて算出するものである。なお、図中の各格子は１ピクセルを表している。実際には図示の断層画像と比較して相当小さいが、説明の便宜上、大きく表記している。この風上勾配法は、微小変形前後における輝度分布の変化を輝度勾配と移動量によって定式化する手法であり、ｆを輝度とすると、微小変形ｆ(x＋Δx,z＋Δz)をテイラー展開する下記式（１６）として表される。
【数16】
【0069】
上記式（１６）は、注目点の変形前後の輝度差が変形前の輝度勾配と移動量によって表されることを示している。なお、移動量(Δx,Δz)については上記式（１６）のみでは決定できないため、検査領域Ｓ１内で移動量が一定と考え、最小二乗法を適用して算出している。
【0070】
上記式（１６）を用いて移動量を算出する際には、右辺の各注目点における移動前後の輝度差は一意にしか求まらない。そのため輝度勾配をどれだけ正確に算出するかが移動量の精度に直結する。輝度勾配の差分化では、一次精度風上差分を用いている。差分化において高次差分を適用すると、多くのデータが必要になり、ノイズが含まれていた際に影響を大きく受けてしまうためである。また、検査領域Ｓ１内の各点を基準とした高次差分では、検査領域Ｓ１外のデータを多く使用することとなり、検査領域Ｓ１そのものの移動量ではなくなってしまうという問題点も存在するからである。
【0071】
輝度勾配を求める際に変形前の風上側の輝度勾配が移動することによって注目点の輝度差が生まれると考えることができるので、変形前は風上側の差分を適用する。ここでいう風上は、実際の移動方向ではなく、ピクセル移動量に対するサブピクセル移動量の向きのことであり、最大相関値ピークに放物線近似を施すことによって風上側を決定する。逆に、変形後の風下側の輝度勾配が逆に移動することによって注目点の輝度差が生じると考えることができるため、変形後は風下側の差分を適用する。
【0072】
変形前の風上差分と変形後の風下差分を用いて２通りの解を求め、それらの平均をとった。さらに、実際には移動量が軸方向に沿わない場合には、変形前や変形後の輝度勾配が注目点と同一軸上に無く、ずれた位置の勾配を求める必要がある。つまり、Ｘ方向の勾配を求めたい場合には、Ｚ方向の移動も考慮して勾配を求めなければならない。そのため、輝度の内挿による輝度勾配の推定をすることで、精度向上を図っている。基本的には変形前（もしくは変形後）の位置を予測し、その位置での勾配を内挿により求める。
【0073】
変形前（後）の注目点の位置は放物線近似を施した際のサブピクセル移動量により求める。その注目点位置が囲まれる８つの座標を用い、それらの比によって輝度勾配を算出する。具体的には、下記式（１７）を用いる。そのようにして算出された輝度勾配と、輝度変化を用いて最小二乗法を適用し移動量を決定した。
【数17】
【0074】
（画像変形法）
上述した風上勾配法までは検査領域Ｓ１の形状は変更せず、正方形を保ったまま変形ベクトルの算出を行っている。しかし、現実には計測対象の変形に合わせて検査領域Ｓ１も変形していると考えられるため、検査領域Ｓ１の微小変形を考慮したアルゴリズムを導入し、変形ベクトル算出を高精度にて算出する必要がある。このため、本実施例ではサブピクセル精度での変形ベクトルの算出に画像変形法を導入している。すなわち、組織変形前の検査領域Ｓ１と組織変形後の伸縮及びせん断変形を考慮した検査領域Ｓ１とで相互相関を実施し、相関値ベースの反復計算によってサブピクセル変形量を決定している。なお、検査領域Ｓ１の伸縮及びせん断変形は線形で近似している。
上述した風上勾配法までは検査領域Ｓ１の形状は変更せず、正方形を保ったまま変形ベクトルの算出を行っている。しかし、現実には計測対象の変形に合わせて検査領域Ｓ１も変形していると考えられるため、検査領域Ｓ１の微小変形を考慮したアルゴリズムを導入し、変形ベクトル算出を高精度にて算出する必要がある。このため、本実施例ではサブピクセル精度での変形ベクトルの算出に画像変形法を導入している。すなわち、組織変形前の検査領域Ｓ１と組織変形後の伸縮及びせん断変形を考慮した検査領域Ｓ１とで相互相関を実施し、相関値ベースの反復計算によってサブピクセル変形量を決定している。なお、検査領域Ｓ１の伸縮及びせん断変形は線形で近似している。
【0075】
画像変形法は一般的に、材料の表面ひずみ計測法に用いられ、ランダムパターンを塗布した材料表面を高空間分解能なカメラで撮影した画像に対して適用される。一方、ＯＣＴ断層像はスペックルノイズを多く含むだけでなく、特に生体組織では組織内の基質や水分の流動も伴い屈折率が変化するため、スペックルパターンに対する変形が大きい。このため、複雑組織にて局所的な変形が発生し、検査領域Ｓ１に膨張、収縮、剪断等の変形が伴う場合には解の収束解が著しく低下する。本手法における検査領域Ｓ１の縮小は、局所的な組織力学特性を検出するために必要不可欠である。そこで画像変形法では、風上勾配法で得られた変形量を収束計算の初期値として採用し、さらに、輝度分布の双三次関数補間によって検査領域Ｓ１を縮小した際でも低ロバスト性を実現している。なお、変形例においては、双三次関数補間以外の補間関数を用いてもよい。
【0076】
より詳細には、以下の手順にて演算を実行する。まず、組織変形前のＯＣＴ断層像の輝度分布に双三次関数補間法を適用し、輝度分布の連続化を実施する。双三次関数補間法とは、ｓｉｎｃ関数を区分的に三次関数近似した畳み込み関数を用い、輝度情報の空間連続性を再現する手法である。本来は連続的な輝度分布を画像計測する際には光学系に依存した点広がり関数が畳み込まれるため、ｓｉｎｃ関数を用いた逆畳み込みを行うことにより、本来の連続的な輝度分布が復元される。離散的な一軸信号ｆ(x)の補間をする場合、畳み込み関数ｈ(x)は下記式（１８）にて表される。
【数18】
【0077】
なお、ＯＣＴ計測条件の違いによって輝度補間関数ｈ(x)の形状も変更する必要がある。そこで輝度補間関数ｈ(x)のｘ＝１での微係数ａを可変とし、ａの値を変更することで輝度補間関数ｈ(x)の形状を変更可能なアルゴリズムとした。本実施例では、擬似ＯＣＴ断層像を用いた数値実験による検証結果を元にし、ａの値を決定した。以上のように画像補間をすることで、伸縮及びせん断変形を考慮した検査領域Ｓ１の各点にて、ＯＣＴ輝度値を求めることが可能となる。
【0078】
伸縮及びせん断変形を考慮して算出した検査領域Ｓ１は、図４（Ｃ）に示すように、移動とともに変形を伴う。組織変形前のＯＣＴ断層像におけるある検査領域Ｓ１内の整数ピクセル位置での座標(x,z)が組織変形後に座標(x＊,z＊)に移動すると考えると、x＊,z＊の値は下記式（１９）にて表される。
【数19】
【0079】
ここで、Δx，Δzはそれぞれ検査領域Ｓ１中心から座標x，zまでの距離、u，vはそれぞれx，z方向への変形量、∂u／∂x，∂v／∂zはそれぞれx，z方向における検査領域Ｓ１の垂直方向への変形量、∂u／∂z，∂v／∂xはそれぞれx，ｚ方向における検査領域Ｓ１のせん断方向への変形量である。数値解法にはNewton-Raphson法を用い、６変数(u，v，∂u／∂x，∂u／∂z，∂v／∂x，∂v／∂z)での相関値微係数が０となるように、すなわち最大相関値を得るように反復計算を行う。なお、反復計算の収束性を高めるため、x，z方向の移動量初期値には風上勾配法で得られたサブピクセル移動量を用いる。相関値Ｒに対するヘッセ行列をＨ、相関値対するヤコビベクトルを▽Ｒとすると、１回の反復で得られる更新量ΔＰiは下記式（２０）にて表される。
【数20】
【0080】
収束の判定には、反復計算で随時得られる漸近解が収束解の近傍で十分小さくなることを用いる。しかし、スペックルパターンの変化が激しい領域においては、線形変形では追従できないために正しい収束解が得られない場合がある。その場合、本実施例では風上勾配法によって求めたサブピクセル移動量を採用する。このようにして得られるサブピクセル精度の変形ベクトルを時間微分することにより、変形速度ベクトル分布を算出できる。
【0081】
（時空間移動最小二乗法）
ひずみ速度テンソルの算出には移動最小二乗法（Moving Least Square Method:以下「ＭＬＳＭ」という）を用いる。ＭＬＳＭは、移動量分布を平滑化すると共に、微係数の安定算出を可能とする手法である。ＭＬＳＭにおいて用いる二乗誤差式は下記式（２１）にて表される。
ひずみ速度テンソルの算出には移動最小二乗法（Moving Least Square Method:以下「ＭＬＳＭ」という）を用いる。ＭＬＳＭは、移動量分布を平滑化すると共に、微係数の安定算出を可能とする手法である。ＭＬＳＭにおいて用いる二乗誤差式は下記式（２１）にて表される。
【数21】
【0082】
上記式（２１）において、Ｓ(x,z,t)を最小とするａ～ｋのパラメータを求める。すなわち、近似関数として水平方向ｘ，奥行き方向ｚ，時間方向ｔの３変数２次多項式として下記式（２２）を採用する。そして、最小二乗近似に基づき、下記式（２３）から最適な微係数を算出して平滑化する。
【数22】
【0083】
この微係数を用いることにより、下記式（２４）に示すひずみ速度テンソルを算出することができる。ｆx，ｆzは各軸のひずみ増分を示し、その時間変化量からひずみ速度を算出している。
【数23】
【0084】
上記式（２４）にて算出されたひずみ速度を時間積分することにより、ひずみ量を算出でき、そのひずみ量の振幅（ひずみ振幅ε0）は断層検出できる。また、上記式（１１）に関連して説明したように、超音波の音響放射圧の大きさ（応力振幅σ0）は、制御演算部１４や駆動回路７４にて変調振幅として与えることができるため既知であり、応力とひずみの位相差（損失角δ）も算出可能である。よって、粘弾性パラメータである貯蔵弾性率Ｅ'と損失弾性率Ｅ"を算出できる。このようにして得られた貯蔵弾性率Ｅ'と損失弾性率Ｅ"の断層分布が、組織の粘弾性の断層分布を実質的に示すことになる。制御演算部１４は、その組織の粘弾性を断層可視化する。なお、ひずみ速度のままでも定式化は可能である。つまり、粘弾性パラメータとの対応をとることができる。
【0085】
次に、制御演算部１４が実行する具体的処理の流れについて説明する。
図５は、制御演算部１４により実行される粘弾性断層可視化処理の流れを示すフローチャートである。本処理は、所定の演算周期で繰り返し実行される。制御演算部１４は、光源２、光学機構８，１０および負荷装置１３を駆動制御する。それにより、適正に振幅変調された音響放射圧負荷信号を組織に対して入力するとともに（Ｓ１０）、ＯＣＴによる光干渉信号を取得する（Ｓ１１）。
図５は、制御演算部１４により実行される粘弾性断層可視化処理の流れを示すフローチャートである。本処理は、所定の演算周期で繰り返し実行される。制御演算部１４は、光源２、光学機構８，１０および負荷装置１３を駆動制御する。それにより、適正に振幅変調された音響放射圧負荷信号を組織に対して入力するとともに（Ｓ１０）、ＯＣＴによる光干渉信号を取得する（Ｓ１１）。
【0086】
制御演算部１４は、周波数空間での演算のためにその光干渉信号にフーリエ変換を施して周波数解析する（Ｓ１２）。続いて、光学機構１０での変調周波数に対応づけるようにバンドパスフィルタリング処理を実行して信号ＳＮ比を向上させた後（Ｓ１４）、ヒルベルト変換を実行する（Ｓ１６）。このヒルベルト変換によって得られた解析信号を用いて、自己相関型リアルタイム粘弾性演算処理を実行する（Ｓ１８）。また、これと並行して粘弾性断層分布演算処理を実行する（Ｓ２０）。リアルタイム粘弾性演算処理は、ＯＣＴ断層画像を取得する過程において、対象Ｓの組織における一軸方向（Ｚ方向）の粘弾性をほぼリアルタイムに演算し、可視化する処理である。粘弾性断層分布演算処理は、ＯＣＴ断層画像（二次元画像）を複数取得し、それらに基づいてポスト処理的に組織の粘弾性を断層可視化する処理である。
【0087】
図６は、図５におけるＳ１８のリアルタイム粘弾性演算処理を詳細に示すフローチャートである。制御演算部１４は、Ｓ１６で取得した解析信号を用いて隣接自己相関処理を施して任意座標における位相差を求め（Ｓ３０）、ドップラー周波数（ドップラー角周波数シフト量）を得る（Ｓ３２）。続いて、空間平滑化や空間平均化処理により、組織における変形速度の断層分布を算出する（Ｓ３４）。この変形速度分布に基づき、最小二乗法を用いてひずみ分布（ひずみ速度分布）を算出し（Ｓ３６）、上記式（１１）に基づいて貯蔵弾性率Ｅ'および損失弾性率Ｅ"のＺ方向分布を算出し（Ｓ３８）、一軸方向粘弾性を断層可視化する（Ｓ４０）。
【0088】
図７は、図５におけるＳ２０の粘弾性断層分布演算処理を詳細に示すフローチャートである。制御演算部１４は、解析信号の包絡線（振幅）に基づき、干渉光強度を算出する（Ｓ５０）。そして、ＯＣＴ断層画像が所定数取得されると（Ｓ５２のＹ）、連続撮影された時刻の異なる前後２枚の断層画像Ｉ(x,z,t)とＩ(x,z,t+Δt)を読み込む（Ｓ５４）。続いて、再帰的相互相関法による処理を実行する。ここではまず、最小解像度（最大サイズの検査領域）での相互相関処理を実行し、相関係数分布を求める（Ｓ５６）。続いて、隣接相互相関乗法により、隣接する相関係数分布の積を演算する（Ｓ５８）。このとき、標準偏差フィルタ等の空間フィルタにより過誤ベクトルの除去をし（Ｓ６０）、最小二乗法等により除去ベクトルの内挿補間を実行する（Ｓ６２）。続いて、検査領域を小さくすることによって解像度を上げて相互相関処理を継続する（Ｓ６４）。すなわち、低解像度での参照ベクトルを基に相互相関処理を実行する。このときの解像度が予め定める最高解像度でなければ（Ｓ６６のＮ）、Ｓ５８に戻る。
【0089】
そして、Ｓ５８～Ｓ６６の処理を繰り返し、最高解像度での相互相関処理が完了すると（Ｓ６６のＹ）、サブピクセル解析を実行する。すなわち、最高解像度（最小サイズの検査領域）での変形ベクトルの分布に基づき、風上勾配法によるサブピクセル移動量を演算する（Ｓ６８）。そして、このとき算出されたサブピクセル移動量に基づき、画像変形法によるサブピクセル変形量を演算する（Ｓ７０）。続いて、最大相互相関値によるフィルタ処理により過誤ベクトルの除去をし（Ｓ７２）、最小二乗法等により除去ベクトルの内挿補間を実行する（Ｓ７４）。そして、このようにして得られた変形ベクトルの時間微分を行い、その断層画像について変形速度ベクトルの断層分布Ｕ(x,z,t)，Ｗ(x,z,t)を算出する（Ｓ７６）。ここで、Ｕ(x,z,t)は変形速度ベクトル分布のｚ方向成分であり、Ｗ(x,z,t)は変形速度ベクトル分布のｘ方向成分である。以上の処理を後続の断層画像についても実行する（Ｓ７８のＮ）。そして、設定数の断層画像について変形ベクトルの断層分布が求まると（Ｓ７８のＹ）、粘弾性演算提示処理を実行する（Ｓ８０）。ＯＣＴ断層画像が所定数取得されていない場合（Ｓ５２のＮ）、本処理を一旦終了する。
【0090】
図８は、図７におけるＳ８０の粘弾性演算提示処理を詳細に示すフローチャートである。制御演算部６は、時空間移動最小二乗法により、Ｓ７６にて算出した変形速度ベクトル分布の平滑化を実行する（Ｓ９０）。そして、その平滑化された変形速度ベクトルにおいて、音響放射圧による加振周波数ω（振幅変調周波数）に同期する成分のみを抽出する（Ｓ９２）。そして、変形速度ベクトルに対して空間微分をすることにより、ひずみ速度テンソルを算出する（Ｓ９４）。
【0091】
続いて、Ｓ９４で得られたひずみ速度を時間積分することにより、ひずみ量を算出し、そのひずみ量の振幅（ひずみ振幅ε0）の断層分布を演算する。そして上述のように、超音波の音響放射圧の大きさ（応力振幅σ0）と、応力とひずみの位相差（損失角δ）とを用いて貯蔵弾性率Ｅ'および損失弾性率Ｅ"の断層分布を演算し（Ｓ９６）、組織の粘弾性を断層可視化する（Ｓ９８）。なお、ひずみ速度のままでも定式化は可能である。つまり、粘弾性パラメータとの対応をとることができる。
【0092】
次に、本実施例の効果を評価するために行った検証実験について説明する。
図９および図１０は、生体組織を模した試料についての実験結果を表す図である。図９は硬質ゲル試料を用いたときの結果を示し、図１０は軟質ゲル試料を用いたときの結果を示す。各図において（Ａ）はＯＣＴ断層像を示し、（Ｂ）は変形速度（上記式（７）参照）の断層分布を示す。
図９および図１０は、生体組織を模した試料についての実験結果を表す図である。図９は硬質ゲル試料を用いたときの結果を示し、図１０は軟質ゲル試料を用いたときの結果を示す。各図において（Ａ）はＯＣＴ断層像を示し、（Ｂ）は変形速度（上記式（７）参照）の断層分布を示す。
【0093】
本実験では、光源２として中心波長λｃ＝１３１７ｎｍ、半値全幅Δλ＝１００ｎｍのＳＬＤ広帯域光源を用いた。参照鏡２６の振れ角を±１．９ｄｅｇ、走査周波数を４ｋＨｚに設定し、試料における奥行き走査範囲（Ｚ方向検出範囲）を９００μｍとした。サンプリング周波数を１０ＭＨｚとし、Δｘ＝２．５μｍ間隔で奥行き１軸干渉信号を取得し、１５００×９００（μｍ）の断層像を得た。一方、超音波デバイス７２を１ＭＨｚの基本波にて駆動し、印可音圧をハイドロフォンにて実験前に実測した。
【0094】
対象Ｓとして、ヒト皮膚の光学特性を模擬したゲル試料（１５ｍｍ×１５ｍｍ×５ｍｍ）を用いた。なお、ゲル試料として、ゲル硬化剤の重量パーセントが０．５％の軟質ゲル試料と、ゲル硬化剤の重量パーセントをその４倍にした硬質ゲル試料の２種類を用意し、剛性依存性について断層可視化評価を行った。
【0095】
図９および図１０に示す実験結果によれば、高剛性の硬質ゲル試料では変形速度が小さく、低剛性の軟質ゲル試料では変形速度が大きいことが分かる。このことは、加振力に対して生体組織の低剛性部分は大きく振動し（振幅大）、高剛性部分は小さく振動する（振幅小）という実態に沿ったものと言える。なお、軟質ゲル試料において変形速度に分布があるのは、音響放射圧に分布があるためと考えられる。
【0096】
図１１および図１２は、剛性を層状に変化させた層状試料についての実験結果を表す図である。図１１は音響放射圧を付与しない状態での実験結果を示し、図１２は音響放射圧を付与した状態での実験結果を示す。各図において（Ａ）はＯＣＴ断層像を示し、（Ｂ）は変形速度の断層分布を示す。層状試料として、表面から２００μｍ程度までを軟質ゲル層、それより下１０００μｍ程度までを硬質ゲル層としたものを用いた。
【0097】
図１１の実験結果によれば、音響放射圧が無負荷の状態では、試料の両層間で違いが見られないことが確認できる。一方、図１２の実験結果によれば、軟質ゲル層と硬質ゲル層とで変形速度が異なる（軟質ゲル層のほうが変形速度（振幅）が大きい）ことが分かる。このことは、逆に変形速度の断層可視化によって測定対象における剛性の断層分布（弾性の分布，すなわち上記式（１１）におけるＥ'の断層分布）を評価できることを意味する。
【0098】
図１３は、粘性評価に関する実験結果を表す図である。（Ａ）はＯＣＴ断層像を示し、（Ｂ）は変形速度の断層分布を示す。（Ｃ）は、（Ｂ）において奥行方向（Ｚ方向）に変形速度の位相遅れがあることを二点鎖線にて示す。
本実験では、図１２に示したものと同様の試料を用いた。ＯＣＴにおいてＸ方向走査を行いつつ、超音波による音響放射圧の振幅を周期的に変化させ、得られる変形速度を断層可視化した。ここでは、その振幅変調周波数を２．４Ｈｚとした。
本実験では、図１２に示したものと同様の試料を用いた。ＯＣＴにおいてＸ方向走査を行いつつ、超音波による音響放射圧の振幅を周期的に変化させ、得られる変形速度を断層可視化した。ここでは、その振幅変調周波数を２．４Ｈｚとした。
【0099】
この実験結果により、奥行き方向に向けて変形速度の位相が遅れる（上層と下層との間で位相がずれる）ことが分かった（二点鎖線の位置が同位相）。これは、音響放射圧により試料に負荷される応力と、それにより発生する変形（ひずみ）との間に位相差（損失角δ）があることを示し、試料に粘性があることを示している。すなわち、以上の実験結果から、粘性の断層分布（粘性係数の分布、すなわち上記式（１１）におけるＥ"の断層分布）を評価できることを意味する。上記実施例によって生体組織の粘弾性の断層検出が可能であることが示唆された。
【0100】
以上説明したように、本実施例によれば、音響放射圧による非接触の荷重負荷と、ＯＣＴによる非接触・非侵襲の断層計測により、対象Ｓの組織の粘弾性をマイクロ断層可視化する。このため、対象Ｓが厳格な品質保証を要するものであったとしても、その汚染を防止でき、高精度に力学特性評価を行うことができる。このことが、再生組織のバイオメカニクスの解明にも寄与すると考えられる。
【0101】
以上、本発明の好適な実施例について説明したが、本発明はその特定の実施例に限定されるものではなく、本発明の技術思想の範囲内で種々の変形が可能であることはいうまでもない。
【0102】
図１４は、変形例に係る荷重負荷方法を表す図である。（Ａ）は第１変形例を示し、（Ｂ）は第２変形例を示し、（Ｃ）は第３変形例を示す。
第１変形例では、光学機構８と超音波デバイス１７２とが対象Ｓに対して同じ側に配置されつつ、ＯＣＴによる光軸が焦点Ｆを通るように設定される。この構成では、振動子アレイ１９０がドーナツ状に形成され、その中央に孔１９２が形成される。ＯＣＴの光は、この孔１９２を通って対象Ｓに照射される。このような構成によれば、負荷装置の支持台に上記実施例のような貫通孔を設ける必要がない。
第１変形例では、光学機構８と超音波デバイス１７２とが対象Ｓに対して同じ側に配置されつつ、ＯＣＴによる光軸が焦点Ｆを通るように設定される。この構成では、振動子アレイ１９０がドーナツ状に形成され、その中央に孔１９２が形成される。ＯＣＴの光は、この孔１９２を通って対象Ｓに照射される。このような構成によれば、負荷装置の支持台に上記実施例のような貫通孔を設ける必要がない。
【0103】
第２変形例では、光学機構８と超音波デバイス７２とが対象Ｓに対して反対側に配置され、ＯＣＴによる断層計測の過程で光軸が焦点Ｆを通らない場合が存在する。例えば、超音波デバイス７２の焦点Ｆを固定し、ＯＣＴのＸ方向走査を行う場合が該当する。本変形例では、焦点ＦからＸ方向に所定距離ずれた位置の組織の変形が検出される。一方、第３変形例では、光学機構８と超音波デバイス７２とが対象Ｓに対して同じ側に配置されつつ、ＯＣＴによる断層計測の過程で光軸が焦点Ｆを通らないように設定される。
【0104】
上記実施例や第１変形例のように、ＯＣＴによる光軸が焦点Ｆを通るように設定すると、ＯＣＴによる検出箇所について、音響放射圧の振幅および位相が把握し易いため、演算処理が容易となり、その精度も高くし易い。ただし、機構やデバイスの設置スペースや作動スペースを考慮すると、第２，第３変形例のほうがシステムを構築し易いといったメリットがある。
【0105】
なお、第２，第３変形例の構成を用いることにより、音響放射圧による組織の変形ではなく、せん断波の伝播速度を検出してもよい。そして、その伝播速度に基づいて組織の粘弾性を判断してもよい。すなわち、音響放射圧による非接触エラストグラフィーを利用し、組織の粘弾性をＯＣＴにより検出してもよい。
【0106】
上記実施例および変形例では、支持台７０を固定し、ＯＣＴの光軸および超音波の焦点を移動させる例を示した。変形例においては、支持台７０をＸ方向やＹ方向に移動制御し、対象Ｓにおける超音波の焦点やＯＣＴの検出位置を変化させてもよい。
【0107】
上記実施例では、測定対象を再生組織（培養細胞）とする例を示した。変形例においては、上記装置を生体の再生治療に適用し、組織の生着性をモニタリングしてもよい。上記装置は、軟骨診断や皮膚診断等の医療分野だけでなく、美容整形分野やコズメティック産業分野など、様々な用途および分野に応用できる。
【0108】
上記実施例では、非接触負荷装置として超音波デバイスを採用し、音響放射圧により組織変形を誘発する例を示した。変形例においては、電気伝導性や磁性を有する粒子（分子標的薬など）を生体組織に取り込ませ、電磁界を照射することにより組織変形を誘発してもよい。
【0109】
上記実施例では、貯蔵弾性率Ｅ'および損失弾性率Ｅ"から粘弾性評価を行う例を示したが、これら貯蔵弾性率Ｅ'および損失弾性率Ｅ"から粘弾性物理量を表現する値を新たに定義し、それに基づいて粘弾性評価を行ってもよい。
【0110】
上記実施例では、ＯＣＴによる断層画像を二次元で取得する例を示したが、三次元で取得してもよい。すなわち、奥行方向（Ｚ方向）とＸ方向のみならず、Ｙ方向に走査し、組織の粘弾性を断層可視化してもよい。
【0111】
なお、本発明は上記実施例や変形例に限定されるものではなく、要旨を逸脱しない範囲で構成要素を変形して具体化することができる。上記実施例や変形例に開示されている複数の構成要素を適宜組み合わせることにより種々の発明を形成してもよい。また、上記実施例や変形例に示される全構成要素からいくつかの構成要素を削除してもよい。
【図1】
【図2】
【図3】
【図4】
【図5】
【図6】
【図7】
【図8】
【図9】
【図10】
【図11】
【図12】
【図13】
【図14】