JPWO2016027874A1 - 応力可視化装置および力学物性値可視化装置 - Google Patents

Abstract

制御演算部１１４は、光学変調器１１０による偏光状態を変化させることにより測定対象Ｗに照射する偏光波の位相をシフトさせ、その偏光波の位相シフトごとに得られる干渉信号の位相差に基づき、その干渉信号の位相差断層分布の空間勾配を演算し、その空間勾配の断層分布を応力の断層分布と対応づけるようにして表示装置１１６に可視化表示させる。制御演算部１１４は、光干渉信号に基づいて測定対象の断層位置での変形速度ベクトル分布を演算し、さらにひずみ速度テンソルの断層分布を演算する。そして、応力の断層分布とひずみ速度テンソルの断層分布とに基づいて力学物性値の断層分布を演算し、可視化表示させる。

Description

本発明は応力可視化装置に関し、特に測定対象に負荷される応力又は内在する応力の分布を光学的手法により測定して可視化する装置に関する。本発明はまた、測定対象の力学物性値の断層分布を可視化する力学物性値可視化装置に関する。

近年、航空・宇宙分野をはじめ、自動車，電子実装部品，医療器具など様々な分野において複合材料が急速に普及し、その形状やサイズの多様化および微細化が進んでいる。このような複合材料は、母材と強化材の積層面を有することから、界面剥離やトランスバースクラックなどマイクロスケールの内部欠陥が複雑に関係して破壊に至ることがある。このため、材料内部における応力集中などの力学的状態を定量評価可能な非破壊計測手法の開発が望まれている。

このような計測手法として、光学的手法により検査対象表面の残留応力を非破壊非接触にて測定する方法が提案されている（例えば特許文献１参照）。この手法では、検査対象表面に加熱用のレーザ光と変形量測定用のレーザ光とを照射し、検査領域における熱応力開放前後の全変形量を光干渉計測により取得し、その全変形量から熱ひずみによる変形量を差し引くことにより残留応力による変形量を算出する。その残留応力による変形量を加熱前後のヤング率の差により除算することにより、残留応力を算出するものである。

特開２００８−１３４１８６号公報

しかしながら、特許文献１に記載の技術では、検査対象表面の応力分布が分かるものの、その検査対象内部の応力分布を知ることはできない。また、応力測定のために干渉計測装置に加えて加熱用装置が必要となるため、計測システム全体として大がかりとなる。さらに、検査対象に対して外的に熱応力を負荷するため、仮にそれを弾性変形範囲内に抑えるとしても、検査対象の物性に変性をきたす可能性がある点で問題があった。

本発明はこのような課題に鑑みてなされたものであり、その目的の一つは、測定対象に負荷される応力又は内在する応力の分布を光学的手法により簡易に測定および可視化可能な技術を提供することにある。

本発明のある態様は、測定対象に作用又は内在する応力の分布を断層可視化する応力可視化装置である。この応力可視化装置は、低コヒーレンス光又は波長走査した光を出射する光源と、光源から出射された光を直線偏光する偏光子と、直線偏光された光を測定対象を経由するオブジェクトアームと、参照鏡を経由するリファレンスアームとに分けるビームスプリッタと、オブジェクトアームに設けられ、直線偏光された光の偏光状態を変化させることが可能な光学変調器と、測定対象にて反射した物体光と参照鏡にて反射した参照光とが重畳された干渉光を、その水平偏光成分と垂直偏光成分とに分けて検出する光検出装置と、光学変調器を制御することにより偏光状態を調整する一方、測定対象の奥行き方向の位置とその直交方向の位置とにより特定される断層位置における応力分布を、光検出装置が検出した干渉信号の水平偏光成分と垂直偏光成分との位相差に基づいて演算する制御演算部と、制御演算部の演算結果に基づいて、測定対象における応力分布を断層可視化する態様で表示する表示装置と、を備える。

制御演算部は、光学変調器による偏光状態を変化させることにより測定対象に照射する偏光波の位相をシフトさせ、その偏光波の位相シフトごとに得られる干渉信号の位相差に基づき、その干渉信号の位相差断層分布の空間勾配を演算し、その空間勾配の断層分布を応力の断層分布と対応づけるようにして表示装置に可視化表示させる。

この態様によると、測定対象の応力分布の計測のために、いわゆる偏光感受型光干渉断層計（Polarization Sensitive Optical Coherence Tomography：以下「ＰＳ−ＯＣＴ」という）が利用される。それにより、測定対象の内部における２偏波位相差分布を断層計測でき、その断層画像に表れるべき干渉縞の空間周波数から応力の分布を算出し、可視化することができる。すなわち、測定対象に対して応力測定のための別途の応力を負荷する必要がなく、応力断層分布を簡易に測定することができる。そして特に、偏光波の位相をシフトさせることにより、測定対象の内部の任意位置での応力計測をＯＣＴの空間分解能にて実施できるようになる。

本発明によれば、測定対象に負荷される応力又は残留する応力の分布を光学的手法により簡易に測定可能な技術を提供することができる。

第１実施形態に係る応力可視化装置の構成を表す概略図である。 ＰＳ−ＯＣＴによる応力断層計測の原理を表す説明図である。 ＰＳ−ＯＣＴによる応力断層計測の原理を表す説明図である。 ＰＳ−ＯＣＴによる応力断層計測の原理を表す説明図である。 光学変調器への印加電圧を変化させて位相シフトを行った場合のＳ３/Ｓ０の変化を演算した結果を例示する図である。 光学変調器への印加電圧を変化させて位相シフトを行った場合の特定の計測位置(x,z)におけるＳ３/Ｓ０の変化を示す図である。 ＰＳ−ＯＣＴを用いた応力計測に用いられる演算マップを示す図である。 実施例に係る応力可視化装置の全体構成を表す図である。 制御演算部により実行される応力可視化処理の流れを概略的に示す図である。 実験に使用された測定対象Ｗおよびその応力負荷方向を示す図である。 応力可視化状態を示す画面例を示す図である。 再帰的相互相関法による処理手順を概略的に示す図である。 サブピクセル解析による処理手順を概略的に示す図である。 力学物性値の断層可視化処理の手順を表す図である。 ひずみ速度可視化画面を示す図である。 制御演算部４により実行される物理量断層可視化処理の流れを示すフローチャートである。 図１６におけるＳ２１４のひずみ速度分布の演算可視化処理を詳細に示すフローチャートである。

［第１実施形態］
以下、本実施形態に係る応力可視化方法の原理について説明する。
図１は、第１実施形態に係る応力可視化装置の構成を表す概略図である。本実施形態の応力可視化装置は、２偏波を扱うＰＳ−ＯＣＴを利用して測定対象Ｗの応力断層計測を行うものであり、光源１００、偏光子１０２、オブジェクトアーム１０４、リファレンスアーム１０６、ビームスプリッタ１０８、電気光学変調器（Electro-Optic Modulator：以下「ＥＯＭ」と表記する）１１０、光検出装置１１２、制御演算部１１４および表示装置１１６を備える。なお、図示の例では、マイケルソン干渉計をベースとした光学系が示されているが、マッハツェンダー干渉計その他の光学系を採用することもできる。

光源１００には、低コヒーレンス光（時間的に低コヒーレンスな光）を出射可能なものが採用される。なお、図示の例では、光遅延走査方式の分類におけるＴＤ−ＯＣＴ（Time Domain OCT）が示されているが、ＳＤ−ＯＣＴ（Spectral Domain OCT）、ＳＳ−ＯＣＴ（Swept Source OCT）その他のＯＣＴを用いることもできる。

ここで、「ＴＤ−ＯＣＴ」は、リファレンスアーム１０６において参照鏡１１８を駆動して光遅延走査を行うことで、オブジェクトアーム１０４の後方散乱光との干渉光を計測し、測定対象Ｗの深さ方向の反射光強度を得るものである。「ＳＤ−ＯＣＴ」は、ＴＤ−ＯＣＴのような機械的な光遅延走査を行うことなく、分光器にて干渉光を分光し、フーリエ変換を行うことにより反射光強度を得るものである。「ＳＳ−ＯＣＴ」は、同様に機械的な光遅延走査を行うことなく、擬似的に低コヒーレンス光（広帯域光源）とするとともに、フーリエ変換により反射光強度を得るものである。

測定対象Ｗとしては、例えば高分子樹脂材料等の高分子基材料が挙げられる。プラスチックなどの結晶性材料でもよい。本装置は、特に非破壊検査をマイクロスケールにて行うものに有効である。内部に反射する構造（界面など）がある材料であればよく、ナノ粒子などの散乱体が樹脂に内在する複合材料など、透明でない材料にも有効である。

偏光子１０２は、光源１００から出射された光を整えて所望の直線偏光を得るために設置される。その直線偏光は、ビームスプリッタ１０８にて分けられ、その一方がオブジェクトアーム１０４に導かれ、他方がリファレンスアーム１０６に導かれる。

オブジェクトアーム１０４に導かれた直線偏光は、ＥＯＭ１１０を経ることでその偏光状態が調整され、円偏光や楕円偏光となって測定対象に照射される。この照射光は、測定対象Ｗにおいてリファレンスアーム１０６の光学距離と一致した任意断面にて後方散乱光として反射され、物体光としてビームスプリッタ１０８に戻る。一方、リファレンスアーム１０６に導かれた直線偏光は、参照鏡１１８にて反射され、参照光としてビームスプリッタ１０８に戻る。これら物体光と参照光とがビームスプリッタ１０８にて合波（重畳）され、その干渉光が光検出装置１１２により検出される。なお、マッハツェンダー干渉計を採用する場合、物体光と参照光とはビームスプリッタに戻ることなく合波されることになる。

光検出装置１１２は、干渉光を水平偏光成分Ｅ_Ｈと垂直偏光成分Ｅ_Ｖの２偏波に分波し、干渉光強度の水平偏光成分Ｉ_Ｈ(x,z)と垂直偏光成分Ｉ_Ｖ(x,z)とを検出する。ここで、測定対象Ｗに照射される光の光軸方向をｚ軸、その光軸と直交するようにｘ軸およびｙ軸として空間座標を定義する。ｚ軸は測定対象Ｗの奥行き方向を示し、ｘ軸は測定対象Ｗに対する照射光の走査方向を示す。本実施形態では、測定対象Ｗのx-ｚ平面に沿った断層画像を可視化するため、干渉光強度もx,z成分に依存した形で検出される。

制御演算部１１４は、測定対象Ｗのx-ｚ平面上で特定される断層位置に応じた応力分布を、光検出装置１１２の検出値に基づいて演算する。また、この演算により計測される断層画像の解像度を高めるために、測定対象Ｗに照射される光の偏光状態（進相軸と遅相軸との位相差）を適宜シフトさせる位相シフト法が採用される。この位相シフトは、ＥＯＭ１１０を制御することにより行われる。

表示装置１１６は、例えば液晶ディスプレイからなり、制御演算部１１４にて演算された応力分布を測定対象Ｗの断層画像の態様で画面に表示させる。

以下、このような演算処理の原理について順に説明する。図２〜図４は、ＰＳ−ＯＣＴによる応力断層計測の原理を表す説明図である。本実施形態の断層計測は、３次元的に応力分布を計測する光弾性散乱光法にＰＳ−ＯＣＴを適用することにより行われる。光弾性散乱光法は、測定対象Ｗに偏光波を照射したときの後方散乱光を用いてその測定対象Ｗの内部の光波の位相差を計測することにより、３次元的に応力分布を計測する手法である。

図２に示すように、この光弾性散乱光法では、３次元応力場である測定対象Ｗに対して空間座標(x,y,z)が設定される。偏光波の入射軸をｚ軸とし、これに垂直となるようｘ軸，ｙ軸を設定する。測定対象Ｗが複屈折性を有する高分子材料であるため、その測定対象Ｗに対して円偏光を入射すると、その偏光波の偏光状態が変化する。この偏光状態の変化は、主にｚ軸に垂直なｘ−ｙ断面の屈折率変化によるものと考えられる。

このため、測定対象Ｗに作用する主応力が変化した場合、あるいは局所的に主応力の変化が内在している場合には、それが干渉光の水平偏光成分Ｅ_Ｈと垂直偏光成分Ｅ_Ｖとの位相差Φ(x,y,z)として表れる。すなわち、その２偏波の位相差Φ(x,y,z)にしたがってその２偏波に対応する干渉信号が得られる。この位相差Φ(x,y,z)は、下記式（１）として表すことができる。

ここで、λは測定対象Ｗへの入射光の波長であり、σ'(x,y,z)は、測定対象Ｗのｘ−ｙ断面に作用する主応力σ_１とσ_２との差（主応力差）である。Ｃは光弾性定数である。すなわち、測定対象Ｗの内部を進行する光波の位相差Φ(x,y,z)は、力学複屈折効果に基づいて主応力差σ'(x,y,z)の干渉信号を得る任意断面までの往復積分に比例する。

一方、本実施形態のＰＳ−ＯＣＴにおいては、干渉光強度の水平偏光成分Ｉ_Ｈ(x,z)と垂直偏光成分Ｉ_Ｖ(x,z)とについて、下記式（２）および（３）が成立する。

ここで、Ｒ(x,z)は反射率分布、Φ(x,z)は位相差Φ(x,y,z)のx,y成分、ｌ_ｃは出射光のコヒーレンス長である。Δｌ_０は、オブジェクトアーム１０４とリファレンスアーム１０６との光路長差である。

このとき、偏光状態を定式化したストークスパラメータの成分のうち、Ｓ_０とＳ_３を下記式（４）および（５）に基づいて算出することができる。ここで、Ｓ_０は計測される光全体の強度を示し、Ｓ_３は直交成分の直線偏光の強度差を示すが、このことは周知であるため、詳細な説明は省略する。

ここで、Ｓ_３(x,z)をＳ_０(x,z)で除算することにより、下記式（６）に示すように、反射強度の情報を除去し、力学複屈折効果に依存する偏光波の位相差Φ(x,z)のみに依存する余弦関数を抽出することができる。

すなわち、上記式（１）により、位相差分布Φ(x,z)は、光軸に垂直なｘ−ｙ断面の主応力差分布σ'(x,z)の積分値に比例することが分かる。一方、上記式（６）により、ストークスパラメータマップＳ_３/Ｓ_０(x,z)の強度分布が、位相差分布Φ(x,z)の余弦関数の縞模様として表れることが理解できる。このため、図３に示すように、Ｓ_３/Ｓ_０(x,z)における干渉縞の空間周波数分布を検出（空間勾配を算出）することにより、主応力差分布σ'(x,z)を断層可視化することができる。なお、位相差分布Φ(x,z)は、ストークスパラメータマップＳ_３/Ｓ_０(x,z)を表す余弦関数の位相であると捉えることもできる。

図４には、測定対象Ｗに対してｘ方向の引張応力を作用させた場合に得られる断層画像が例示されている。図示のように、測定対象Ｗへの負荷応力が大きくなるほど、干渉縞の空間周波数は高くなる。このため、その断層画像により測定対象Ｗの応力分布が可視化されることになる。なお、図示の例では、測定対象Ｗに対して外的負荷がかかった場合が示されているが、測定対象Ｗに外的負荷がかからず、内部に残留応力が作用している場合にも、その残留応力を断層計測して可視化することが可能である。

なお、本実施形態ではビーム走査をｘ方向に１次元的に行うため、２次元断層計測に基づく２次元応力分布を可視化することになるが、ビームをｘ−ｙ平面にて２次元的に走査させれば３次元断層計測を行うことができ、３次元応力分布を可視化することができる。

以上の手法において、詳細な応力分布を得るためには、ストークスパラメータマップＳ_３/Ｓ_０(x,z)における干渉縞の空間周波数を高空間分解能にて算出する必要がある。そこで本実施形態では、以下に述べる位相シフト法を組み入れる。この位相シフト法は、干渉縞の位相を変化させた画像を用い、その位相の変化を参照しつつ非線形パラメータ同定を行うことで干渉縞の空間周波数を解析する手法である。

この位相シフト処理においては、ストークスパラメータマップＳ_３/Ｓ_０(x,z)における干渉縞の位相を変化させた断層像を撮影する。その際、ＥＯＭ１１０を用いて円偏光から位相差δ_EOM(V)を変化させた光を入射させる。このとき、取得されるストークスパラメータマップＳ_３/Ｓ_０(x,z)は、上記式（６）を変形して下記式（７）のようになる。

ここで、比例係数Ａ(x,z)と定数Ｂ(x,z)を定義しているが、これは位相シフト法を適用した際に位相差分布Φ(x,z)を求め、その位相差分布Φ(x,z)の空間勾配を求め、その空間勾配から主応力差σ'(x,z)をより正確に得るためのものである。これらＡ(x,z)およびＢ(x,z)については、例えば最小二乗法を用いて求めることができる。なお、位相差分布Φ(x,z)は、ストークスパラメータマップＳ_３/Ｓ_０(x,z)における位相断層分布Φ(x,z)として算出される。

図５は、ＥＯＭ１１０への印加電圧を変化させて位相シフトを行った場合のＳ_３/Ｓ_０の変化を演算した結果を例示する図である。（Ａ）は印加電圧を５０Ｖとした場合、（Ｂ）は１００Ｖとした場合、（Ｃ）は１５０Ｖとした場合を示す。各図の横軸は測定対象Ｗにおけるｚ軸方向の深さを示し、縦軸はＳ_３/Ｓ_０(x,z)の値を示す。図中の実線は各印加電圧を示す。図中の波線は参照電圧として０Ｖ、つまりＥＯＭ１１０を動作させなかった場合の演算結果を示している。この演算においては、ストークスパラメータマップＳ_３/Ｓ_０(x,z)をｘ軸方向に平均化し、ｚ軸方向に関する位相差変化が算出されている。

ここでは、測定対象Ｗに対してｘ軸方向に４．５Mpaの引張応力を外的に負荷している。また、印加電圧が０Ｖのときに、測定対象Ｗへの入射光の偏光状態が円偏光になるように偏光状態を設定した。図示の結果より、ＥＯＭ１１０への印加電圧を変化させることにより、ストークスパラメータマップＳ_３/Ｓ_０(x,z)の概形を維持したまま、各計測位置(x,z)での偏光波の位相をシフトできることが分かる。

図６は、ＥＯＭ１１０への印加電圧を変化させて位相シフトを行った場合の特定の計測位置(x,z)におけるＳ_３/Ｓ_０の変化を示す図である。同図の横軸はＥＯＭ１１０への印加電圧を示し、縦軸はＳ_３/Ｓ_０(x,z)の値を示す。この演算結果は、上記式（７）において、座標(x,z)を特定座標に固定したときの結果を示す。

この演算結果より、ＥＯＭ１１０への印加電圧に対して特定の座標におけるストークスパラメータマップＳ_３/Ｓ_０(x,z)が正弦関数にしたがって変化することが分かる。すなわち、ＰＳ−ＯＣＴによる計測結果に対して位相シフト法が適用可能であり、偏光波の位相をシフトさせることにより任意の断層位置について応力を算出できることが分かる。それにより、測定対象Ｗの内部の任意位置での応力計測をＯＣＴの空間分解能にて実施でき、応力分布の可視化に際して極めて高い解像度を得ることが可能となる。

図７は、ＰＳ−ＯＣＴを用いた応力計測に用いられる演算マップを示す図である。（Ａ）はある材料に適用される演算マップを示す。同図の横軸はＳ_３/Ｓ_０(x,z)の計測から得られる空間周波数（1/μm）を示し、縦軸はその計測位置における主応力差σ'(x,z)（MPa）を示す。（Ｂ）は同材料について公知の応力−ひずみ線図である。同図の横軸はひずみ（％）を示し、縦軸は応力（MPa）を示す。

図７（Ａ）に示すように、本実施形態では、Ｓ_３/Ｓ_０(x,z)の計測から得られる空間周波数（後述する位相断層分布Φ(x,z)の空間勾配）と、その計測位置における主応力差σ'(x,z)とを対応づけた演算マップを、測定対象Ｗとなりうる材質ごとに用意する。この演算マップは、測定対象Ｗへの負荷応力と、ＰＳ−ＯＣＴにより得られる干渉縞の空間周波数との相関を示すキャリブレーション曲線を実験的に求めることにより得られる。図７（Ａ）および（Ｂ）を対比すると、両者の曲線が高精度に一致していることが分かる。

以下、図面を参照しつつ本実施形態を具体化した実施例について詳細に説明する。
図８は、実施例に係る応力可視化装置の全体構成を表す図である。本実施例の応力可視化装置１は、マイケルソン干渉光学系をベースとしたＴＤ−ＰＳ−ＯＣＴを利用して測定対象の応力分布を演算するシステムとして構成されている。応力可視化装置１は、ＯＣＴによる光学的処理を実行する光学系２と、その処理結果に基づいて測定対象Ｗの応力分布を演算する制御演算部４とを備える。光学系２は、光源６、偏光子８、カプラー１０、リファレンスアーム１２、オブジェクトアーム１４および光検出装置１６等の光学要素を含む。各光学要素は、ファイバーにて互いに接続されている。

光源６は、スーパールミネッセントダイオード（Super Luminessent Diode：以下「ＳＬＤ」と表記する）からなる広帯域光源である。ＳＬＤは、高輝度な低コヒーレンス近赤外光を発する。なお、変形例においては、低コヒーレンス光を出射可能な他の光源を採用してもよい。

偏光子８は、光源６からの光を直線偏光にする直線偏光子（Linear Polarizer：以下「ＬＰ」と表記する）であり、光源６とカプラー１０との間に配置されている。カプラー１０は、光源６からの光をリファレンスアーム１２とオブジェクトアーム１４とに向けて分波するビームスプリッタとしての機能と、リファレンスアーム１２およびオブジェクトアーム１４のそれぞれから戻ってきた光を合波して干渉させる機能とを有する。

リファレンスアーム１２には、偏波コントローラ（Polarization Controller：以下「ＰＣ」と表記する）１８、ＬＰ２６、ＥＯＭ２８、集光レンズ２２が順次接続されている。集光レンズ２２は、参照鏡２４に対向配置されている。ＬＰ２６は、カプラー１０を経てリファレンスアーム１２に導かれた光の偏光状態を整える。ＥＯＭ２８は、印加電圧を線形変化させることで、単一の高周波なキャリア周波数を得ることを可能にするものであり、ＥＯＰＭ（Electro‐Optic Phase Modulator）とも呼ばれる（以下「ＥＯＰＭ２８」とも表記する）。これにより、ＯＣＴによる干渉信号を高感度にてヘテロダイン検出することができる。

ＬＰ２６は、ＥＯＰＭ２８により偏光状態を変調しても、リファレンスアーム１２から戻ってくる光の強度が変化しないようＰＣ１８と協働する。具体的には、ＰＣ１８は、参照鏡２４からの参照光が再びカプラー１０に入射する際の偏光状態を４５度の直線偏光となるように調整する。

オブジェクトアーム１４には、ＰＣ３０、１／４波長板（Quarter Wave Plate：以下「ＱＷＰ」と表記する）４４、ＥＯＭ４６、コリメートレンズ３４、ガルバノ装置３６、対物レンズ３８が順次接続されている。ガルバノ装置３６は、固定鏡４０およびガルバノ鏡４２を含む。対物レンズ３８は、測定対象Ｗに対向配置される。カプラー１０を介してオブジェクトアーム１４に入射された光は、２軸のガルバノ鏡４２によりｘ軸方向に走査されて測定対象Ｗに照射される。測定対象Ｗからの反射光は、物体光として再びカプラー１０に戻り、参照光と重畳されて干渉光として光検出装置１６に送られる。

ＰＣ３０は、カプラー１０からオブジェクトアーム１４に導かれた光の偏光状態を円偏光や楕円偏光となり得るように調整する。ＱＷＰ４４は、４５度の直線偏光を通過させることでその進相軸と遅相軸との位相差を１／４波長ずらし、円偏光に変化させる。ＥＯＭ４６は、印加電圧に応じて透過光の偏光状態を任意に変化させるものであり、ＥＯＡＭ（Electro‐Optic Amplitude Modulator）とも呼ばれる。これにより、ＱＷＰ４４を通過した円偏光を楕円偏光に変化させることができるとともに、印加電圧に応じてその楕円偏光の進相軸と遅相軸との位相差を任意に変化させることができる。

光検出装置１６は、偏光ビームスプリッタ（Polarization Beam Splitter：以下「ＰＢＳ」と表記する）５０、光検出器５２（「第１光検出器」として機能する）、光検出器５４（「第２光検出器」として機能する）を含む。ＰＢＳ５０は、カプラー１０にて合波された干渉光を水平偏光成分Ｅ_Ｈと垂直偏光成分Ｅ_Ｖとに分波し、それぞれを光検出器５２，５４へ出力する。そして、光検出器５２が干渉光強度の水平偏光成分Ｉ_Ｈ(x,z)を検出し、光検出器５４が干渉光強度の垂直偏光成分Ｉ_Ｖ(x,z)を検出する。これらの検出は同時に行われる。その水平偏光成分Ｉ_Ｈ(x,z)は、バンドパスフィルタ（Band-Pass Filter:以下「ＢＰＦ」と表記する）５６、アンプ５８（Lock-in Amp）およびＡ／Ｄ変換器６０を介して制御演算部４に入力される。一方、その垂直偏光成分Ｉ_Ｖ(x,z)は、ＢＰＦ６２、アンプ６４（Lock-in Amp）およびＡ／Ｄ変換器６０を介して制御演算部４に入力される。

制御演算部４は、本実施例ではパーソナルコンピュータからなり、ユーザによる各種設定入力を受け付けるための入力装置６８、画像処理用の演算プログラムにしたがった演算処理を実行する演算処理部７０、および演算結果を表示させる表示装置７２を含む。演算処理部７０は、ＣＰＵ、ＲＯＭ、ＲＡＭ、ハードディスクなどを有し、これらのハードウェア、およびソフトウェアによって、光学系２全体の制御と光学処理結果に基づく画像出力のための演算処理を行うことができる。

制御演算部４は、上述した原理に基づき、ＰＳ−ＯＣＴによる計測結果に対して位相シフト演算を行うことによりストークスパラメータマップＳ_３/Ｓ_０(x,z)の空間周波数を算出する。そして、算出された空間周波数を用いて図７に示した演算マップを参照することにより、測定対象Ｗの任意の位置に作用する応力を算出し、表示装置７２に表示させる。

図９は、制御演算部４により実行される応力可視化処理の流れを概略的に示す図である。この応力可視化処理において、制御演算部４は、まず、ＥＯＭ４６への印加電圧を予め定める初期値（例えば０Ｖの非通電状態）に設定したうえで、ＰＳ−ＯＣＴによる断層画像を取得する。すなわち、干渉信号（干渉光強度）の水平偏光成分Ｉ_Ｈ(x,z)と垂直偏光成分Ｉ_Ｖ(x,z)とを取得する（Ｓ１０）。続いて、ＥＯＭ４６の印加電圧を変化させて偏光波の位相をシフトさせ、断層画像を取得する（Ｓ１２）。この位相シフトによる断層画像の取得を設定回数実行する（Ｓ１４のＮ）。

それにより設定数の断層画像が取得されると（Ｓ１４のＹ）、各断層画像について前処理としてアンサンブル平均を施した後（Ｓ１６）、移動最小二乗法によるノイズ除去処理と平滑化処理を実行する（Ｓ１８）。そして、各断層画像についてストークスパラメータマップＳ_３/Ｓ_０(x,z)を算出する（Ｓ２０）。

そして、各断層画像のストークスパラメータマップＳ_３/Ｓ_０(x,z)について、以下の処理を実行する。すなわち、Ｓ_３/Ｓ_０(x,z)の１軸断層（ＯＣＴにより走査されたｘ座標の一つの断層）に対してフーリエ解析もしくは最小二乗法を適用し、そのｚ方向（奥行き方向）の位相Φ(z)を算出する（Ｓ２２）。この処理をｘ方向（ＯＣＴによる走査方向）の各ｘ座標について行うことにより、位相断層分布Φ(x,z)を算出する（Ｓ２４）。このような処理により、上記式（７）におけるＳ_３/Ｓ_０(x,z)の各係数が求まる。このとき、位相接続処理を行っておく（Ｓ２６）。すなわち、位相が０〜２πで変化を繰り返すため、光軸方向に向かって位相が正に変化していくように位相接続をする。

このようにして得られた位相断層分布Φ(x,z)に対して標準偏差フィルタ等の空間フィルタによるノイズ除去処理を施し（Ｓ２８）、その除去された部分について位相断層分布Φ(x,z)の空間補間を行う（Ｓ３０）。例えば、ある空間領域（関心領域）の位相差分布に多項式近似の最小二乗法を施し、その関心領域における補間関数を求める。

続いて、最小二乗法などにより位相断層分布Φ(x,z)の変化率、つまり位相断層分布Φ(x,z)の空間勾配断層分布(∂Φ/∂x,∂Φ/∂z）を演算する（Ｓ３２）。このようにして得られた空間勾配断層分布(∂Φ/∂x,∂Φ/∂z）を空間周波数とみなし、図７に示した該当する材質の演算マップを参照して応力分布を算出する。そして、その応力分布を断層可視化する（Ｓ３４）。

次に、本実施形態の応力断層可視化処理に基づく実験結果を示す。図１０は、実験に使用された測定対象Ｗおよびその応力負荷方向を示す。（Ａ）は測定対象Ｗとしての試験片を示し、（Ｂ）は荷重の負荷方法を示す。図１１は、応力可視化画面を示す図である。（Ａ）は本実施形態による応力可視化画面を示す。（Ｂ）は比較例として、ＦＥＭにより本実施形態と同じ条件にて応力解析をした結果を示す。（Ｃ）は、（Ａ）の画像を得る前に取得したＰＳ−ＯＣＴによる断層画像を示す。各図の横軸はｘ方向の位置を、縦軸はｚ方向の位置をそれぞれ示している。

本実験においては、試験片Ｗ１として図１０（Ａ）に示すようなダンベル片を用いた。この試験片Ｗ１は、全長１４９mm、最大幅２０mm、厚み１．５mmを有する低密度ポリエチレンからなる。図１０（Ｂ）に示すように、ＰＳ−ＯＣＴによる偏光波の入射方向にｚ軸をとり、それと直交するようにｘ軸およびｙ軸を設定した。試験片Ｗ１の厚み方向（奥行き方向）の端面がｘ−ｙ平面に延在するように設置される。また、偏光波の光軸と直交する走査方向をｘ軸方向とし、試験片Ｗ１の長手方向に一致させている。応力集中が発生するよう、試験片Ｗ１の長手方向中央の下面に深さ２００μmの切欠きＣを形成した。切欠きＣは、ｙ軸方向に延びる直線状に形成されている。この試験片Ｗ１に対して長手方向の引張荷重Ｐを与え、３０分間応力緩和させた後、上述したＰＳ−ＯＣＴによる計測を行い、応力断層可視化を行った。

その結果、図１１（Ａ）に示すように、応力分布の断層画像を可視化することができた。図中右欄に応力値が示されているが、図１１（Ｂ）に示すＦＥＭ解析の結果と精度良く対応していることが確認できる。なお、図１１（Ａ）では応力値そのものを可視表示しているが、図１１（Ｃ）に示すように、Ｓ_３/Ｓ_０(x,z)の断層画像を表示させることにより、断層位置における干渉縞の空間周波数を把握することができ、それにより、応力の高低の分布を把握することができる。

以上に説明したように、本実施形態によれば、ＰＳ−ＯＣＴに位相シフト法を適用することにより、測定対象Ｗの任意の断層位置について応力を算出し、可視化することができる。つまり、測定対象Ｗの内部の任意位置での応力計測をＯＣＴの空間分解能にて実施できるようになる。このため、応力分布の可視化に際して極めて高い解像度を得ることが可能となる。

［第２実施形態］
次に、第２実施形態について説明する。
本実施形態では、上述した応力断層分布に加え、ひずみ速度断層分布を演算し、これらを用いることにより、測定対象Ｗの内部の任意位置での力学物性値の断層分布を算出および可視化する。ここでいう「力学物性値」は、応力およびひずみ速度（又はひずみ）から得られる係数であって、弾性係数や粘性係数を含む。

まず、ひずみ速度断層分布の演算に関し、ＯＣＴを利用したマイクロスケールの力学特性検出法について説明する。この検出法は、測定対象Ｗの変形前後の２枚のＯＣＴ断層画像にデジタル相互相関法を適用することで変形ベクトル分布を算出し、マイクロスケールにて測定対象内部のひずみ速度テンソル分布を断層検出する手法である。

変形ベクトル分布を算出する際には、繰り返し相互相関処理を実施する再帰的相互相関法(Recursive Cross-correlation method)を適用する。これは、低解像度において算出された変形ベクトルを参照し、探査領域を限定するとともに階層的に検査領域を縮小して相互相関法を適用する手法である。これにより、高解像度な変形ベクトルを取得することができる。さらに、スペックル・ノイズ低減法として、隣接検査領域の相関値分布との乗算を行う隣接相互相関乗法(Adjacent Cross-correlation Multiplication)を用いる。そして、乗算されることによって高ＳＮ化した相関値分布から最大相関値を探索する。

また、マイクロスケールにおける微小変形解析では、変形ベクトルのサブピクセル精度が重要となる。このため、輝度勾配を利用する風上勾配法(Up-stream Gradientmethod)と、伸縮および剪断を考慮した画像変形法(Image Deformation method)の両サブピクセル解析法を併用し、変形ベクトルの高精度検出を実現する。なお、ここでいう「風上勾配法」は、勾配法（オプティカルフロー法）の一種である。

このようにして得られた変形ベクトルを時間微分することにより変形速度ベクトル分布を演算し、それをさらに空間微分することによりひずみ速度テンソル分布を演算する。

以下、各手法について詳細に説明する。図１２は、再帰的相互相関法による処理手順を概略的に示す図である。図１３は、サブピクセル解析による処理手順を概略的に示す図である。

（再帰的相互相関法）
図１２（Ａ）〜（Ｃ）は、再帰的相互相関法による処理過程を示している。各図にはＯＣＴにより連続的に撮影される前後の断層画像が示されている。左側には先の断層画像（Image1）が示され、右側には後の断層画像（Image2）が示されている。

相互相関法とは、局所的なスペックル・パターンの類似度を下記式（８）に基づく相関値Ｒ_ｉｊより評価する方法である。そのため、図１２（Ａ）に示すように、連続的に撮影された前後のＯＣＴ画像について、先の断層画像（Image1）に類似度の検査対象となる検査領域Ｓ１が設定され、後の断層画像（Image2）に類似度の探査範囲となる探査領域Ｓ２が設定される。

ここで、空間座標として、光軸方向にＺ軸、光軸と垂直方向にＸ軸を設定している。ｆ(Ｘ_ｉ,Ｚ_ｊ)とｇ(Ｘ_ｉ,Ｚ_ｊ)は、変形前後のＯＣＴ画像において設置された中心位置(Ｘ_ｉ,Ｚ_ｊ)の検査領域Ｓ１(Ｎ_ｘ×Ｎ_ｚピクセル)のスペックル・パターンを表している。

探査領域Ｓ２(Ｍ_ｘ×Ｍ_ｚピクセル)内における相関値分布Ｒ_i,j(ΔＸ,ΔＺ)を算出し、図１２（Ｂ）に示すようにパターンマッチングを行う。実際には、下記式（９）に示すように、最大相関値を与える移動量Ｕ_ｉ,jを変形前後の変形ベクトルとして決定する。

なお、ｆ⁻とｇ⁻はｆ(Ｘ_ｉ,Ｚ_ｊ)とｇ(Ｘ_ｉ,Ｚ_ｊ)の検査領域Ｓ１内の平均値を表している。

本手法では、検査領域Ｓ１を縮小しながら相互相関処理を繰り返して空間解像度を高める再帰的相互相関法を採用している。なお、本実施例では解像度を上げる際には空間解像度が倍になるようにしている。図１２（Ｃ）に示すように、検査領域Ｓ１を１／４に分割し、前階層にて算出された変形ベクトルを参照し、探査領域Ｓ２を縮小している。ここで探査領域Ｓ２も１／４に分割している。再帰的相互相関法を用いることで、高解像度において多発する過誤ベクトルの抑制を可能にしている。このような再帰的相互相関処理を施すことにより、変形ベクトルの解像度を高めることができる。

また、これに加え、下記式（１０）により、演算中の座標を中心とする周囲８つの座標を含む合計９つの変形ベクトルの平均偏差σを用いた閾値を設定し、過誤ベクトルの除去を行い再帰処理に伴う誤差伝播を抑制する。

ここで、Ｕmはベクトル量の中央値を表し、閾値となる係数κは任意に設定した。

（隣接相互相関乗法）
本実施例では、スペックルノイズの影響を受けたランダム性の強い相関値分布から正確な最大相関値を決定する方法として、隣接相互相関乗法を導入している。下記式（１１）により、隣接相互相関乗法では検査領域Ｓ１における相関値分布Ｒi,j(Δx,Δz)と、その検査領域Ｓ１にオーバーラップする隣接検査領域に対するＲi＋Δi,j(Δx,Δz)とＲi,j＋Δj(Δx,Δz)の乗算を行い、新たな相関値分布Ｒ'i,j(Δx,Δz)を用いて最大相関値を検索する。

これにより、相関値同士の乗算によってランダム性を低減させることが可能になる。上述した検査領域Ｓ１の縮小と共に干渉強度分布の情報量も減少するため、スペックル・ノイズを原因とする複数相関ピークの出現が計測精度の悪化を招いていると考えられる。一方、隣接境界同士の移動量には相関があるため、最大相関値座標付近では強い相関値が残存する。この隣接相互相関乗法の導入によって最大相関値ピークが明瞭化され計測精度が向上し、正確な移動座標を抽出することが可能となる。また、この隣接相互相関乗法をＯＣＴの各ステージに導入することで、誤差伝播が抑制され、スペックル・ノイズに対するロバスト性が向上する。それにより、高空間解像度においても高精度な変形ベクトルとひずみ分布の算出が可能となる。

（風上勾配法）
図１３（Ａ）〜（Ｃ）は、サブピクセル解析による処理過程を示している。各図にはＯＣＴにより連続的に撮影される前後の断層画像が示されている。左側には先の断層画像（Image1）が示され、右側には後の断層画像（Image2）が示されている。

本実施例では、サブピクセル解析のために風上勾配法と画像変形法を採用する。最終的な移動量の算出は後述の画像変形法によるが、計算の収束性の問題から、画像変形法より先に風上勾配法を適用する。検査領域サイズが小さく高空間解像度の条件において、サブピクセル移動量を高精度検出する画像変形法及び風上勾配法を適用している。画像変形法におけるサブピクセル移動量の検出が困難な場合において、風上勾配法によりサブピクセル移動量を算出する。

サブピクセル解析では、注目点における変形前後の輝度差が各成分の輝度勾配と移動量によって表される。このため、検査領域Ｓ１内の輝度勾配データより最小二乗法を用いてサブピクセル移動量を決定することができる。本実施例では、輝度勾配を求める際に、サブピクセル変形前の風上側の輝度勾配を与える風上差分法を採用している。

すなわち、サブピクセル解析は様々な手法が存在するが、本実施例では検査領域サイズが小さく高空間解像度の条件においても、サブピクセル移動量を高精度検出する勾配法を採用している。

風上勾配法は、検査領域Ｓ１内の注目点の移動を、図１３（Ａ）に示すピクセル精度に留まらず、図１３（Ｂ）に示すサブピクセル精度にて算出するものである。なお、図中の各格子は１ピクセルを表している。実際には図示の断層画像と比較して相当小さいが、説明の便宜上、大きく表記している。この風上勾配法は、微小変形前後における輝度分布の変化を輝度勾配と移動量によって定式化する手法であり、ｆを輝度とすると、微小変形ｆ(x＋Δx,z＋Δz)をテイラー展開する下記式（１２）として表される。

上記式（１２）は、注目点の変形前後の輝度差が変形前の輝度勾配と移動量によって表されることを示している。なお、移動量(Δx,Δz)については上記式（１２）のみでは決定できないため、検査領域Ｓ１内で移動量が一定と考え、最小二乗法を適用して算出している。

上記式（１２）を用いて移動量を算出する際には、右辺の各注目点における移動前後の輝度差は一意にしか求まらない。そのため輝度勾配をどれだけ正確に算出するかが移動量の精度に直結する。輝度勾配の差分化では、一次精度風上差分を用いている。差分化において高次差分を適用すると、多くのデータが必要になり、ノイズが含まれていた際に影響を大きく受けてしまうためである。また、検査領域Ｓ１内の各点を基準とした高次差分では、検査領域Ｓ１外のデータを多く使用することとなり、検査領域Ｓ１そのものの移動量ではなくなってしまうという問題点も存在するからである。

輝度勾配を求める際に変形前の風上側の輝度勾配が移動することによって注目点の輝度差が生まれると考えることができるので、変形前は風上側の差分を適用する。ここでいう風上は、実際の移動方向ではなく、ピクセル移動量に対するサブピクセル移動量の向きのことであり、最大相関値ピークに放物線近似を施すことによって風上側を決定する。逆に、変形後の風下側の輝度勾配が逆に移動することによって注目点の輝度差が生じると考えることができるため、変形後は風下側の差分を適用する。

変形前の風上差分と変形後の風下差分を用いて２通りの解を求め、それらの平均をとった。さらに、実際には移動量が軸方向に沿わない場合には、変形前や変形後の輝度勾配が注目点と同一軸上に無く、ずれた位置の勾配を求める必要がある。つまり、Ｘ方向の勾配を求めたい場合には、Ｚ方向の移動も考慮して勾配を求めなければならない。そのため、輝度の内挿による輝度勾配の推定をすることで、精度向上を図っている。基本的には変形前（もしくは変形後）の位置を予測し、その位置での勾配を内挿により求める。

変形前（後）の注目点の位置は放物線近似を施した際のサブピクセル移動量により求める。その注目点位置が囲まれる８つの座標を用い、それらの比によって輝度勾配を算出する。具体的には、下記式（１３）を用いる。そのようにして算出された輝度勾配と、輝度変化を用いて最小二乗法を適用し移動量を決定した。

（画像変形法）
上述した風上勾配法までは検査領域Ｓ１の形状は変更せず、正方形を保ったまま変形ベクトルの算出を行っている。しかし、現実には計測対象の変形に合わせて検査領域Ｓ１も変形していると考えられるため、検査領域Ｓ１の微小変形を考慮したアルゴリズムを導入し、変形ベクトル算出を高精度にて算出する必要がある。このため、本実施例ではサブピクセル精度での変形ベクトルの算出に画像変形法を導入している。すなわち、材料の変形前の検査領域Ｓ１と変形後の伸縮及びせん断変形を考慮した検査領域Ｓ１とで相互相関を実施し、相関値ベースの反復計算によってサブピクセル変形量を決定している。なお、検査領域Ｓ１の伸縮及びせん断変形は線形で近似している。

画像変形法は一般的に、材料の表面ひずみ計測法に用いられ、ランダムパターンを塗布した材料表面を高空間分解能なカメラで撮影した画像に対して適用される。一方、ＯＣＴ断層像はスペックルノイズを多く含み、また、複合材料での局所的な変形が発生するため、検査領域Ｓ１に変形が伴う場合には解の収束解が著しく低下する。本手法における検査領域Ｓ１の縮小は、局所的な力学特性を検出するために必要不可欠である。そこで画像変形法では、風上勾配法で得られた変形量を収束計算の初期値として採用し、さらに、輝度分布の双三次関数補間によって検査領域Ｓ１を縮小した際でも低ロバスト性を実現している。なお、変形例においては、双三次関数補間以外の補間関数を用いてもよい。

より詳細には、以下の手順にて演算を実行する。まず、材料変形前のＯＣＴ断層像の輝度分布に双三次関数補間法を適用し、輝度分布の連続化を実施する。双三次関数補間法とは、ｓｉｎｃ関数を区分的に三次関数近似した畳み込み関数を用い、輝度情報の空間連続性を再現する手法である。本来は連続的な輝度分布を画像計測する際には光学系に依存した点広がり関数が畳み込まれるため、ｓｉｎｃ関数を用いた逆畳み込みを行うことにより、本来の連続的な輝度分布が復元される。離散的な一軸信号ｆ(x)の補間をする場合、畳み込み関数ｈ(x)は下記式（１４）にて表される。

なお、ＯＣＴ計測条件の違いによって輝度補間関数ｈ(x)の形状も変更する必要がある。そこで輝度補間関数ｈ(x)のｘ＝１での微係数ａを可変とし、ａの値を変更することで輝度補間関数ｈ(x)の形状を変更可能なアルゴリズムとした。本実施例では、擬似ＯＣＴ断層像を用いた数値実験による検証結果を元にし、ａの値を決定した。以上のように画像補間をすることで、伸縮及びせん断変形を考慮した検査領域Ｓ１の各点にて、ＯＣＴ輝度値を求めることが可能となる。

伸縮及びせん断変形を考慮して算出した検査領域Ｓ１は、図１３（Ｃ）に示すように、移動とともに変形を伴う。材料変形前のＯＣＴ断層像におけるある検査領域Ｓ１内の整数ピクセル位置での座標(x,z)が変形後に座標(x^＊,z^＊)に移動すると考えると、x^＊,z^＊の値は下記式（１５）にて表される。

ここで、Δ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は下記式（１６）にて表される。

収束の判定には、反復計算で随時得られる漸近解が収束解の近傍で十分小さくなることを用いる。しかし、スペックルパターンの変化が激しい領域においては、線形変形では追従できないために正しい収束解が得られない場合がある。その場合、本実施例では風上勾配法によって求めたサブピクセル移動量を採用している。

以上のようにして得られるサブピクセル精度の変形ベクトルを時間微分することにより、変形速度ベクトル分布を算出することができる。そして、その変形速度ベクトル分布を空間微分することにより、ひずみ速度テンソルを算出することが可能となる。

（時空間移動最小二乗法）
ひずみ速度テンソルの算出には移動最小二乗法（Moving Least Square Method:以下「ＭＬＳＭ」という）を用いる。ＭＬＳＭは、移動量分布を平滑化すると共に、微係数の安定算出を可能とする手法である。ＭＬＳＭにおいて用いる二乗誤差式は下記式（１７）にて表される。

上記式（１７）において、Ｓ(x,z,t)を最小とするａ〜ｋのパラメータを求める。すなわち、近似関数として水平方向ｘ，奥行き方向ｚ，時間方向ｔの３変数２次多項式として下記式（１８）を採用する。そして、最小二乗近似に基づき、下記式（１９）から最適な微係数を算出して平滑化する。

この微係数を用いることにより、下記式（２０）に示すひずみ速度テンソルを算出することができる。ｆx，ｆzは各軸のひずみ増分を示し、その時間変化量からひずみ速度を算出している。

以上のようにして得られるひずみ速度分布と、第１実施形態にて得られた応力分布とを用いることにより、力学物性値を断層可視化する処理を実行する。以下、本実施形態を具体化した実施例について説明する。図１４は、力学物性値の断層可視化処理の手順を表す図である。なお、図示の断層画像は、測定対象Ｗとして図１０および図１１に示した試験片Ｗ１を用いた場合の実験結果を示す。

制御演算部４は、測定対象Ｗの断層画像を連続的に取得し（図１４（Ａ），（Ｂ））、上述のようにして演算された変形速度ベクトル分布を順次記憶するとともに（図１４（Ｂ））、各変形速度ベクトル分布を空間微分して得られたひずみ速度テンソル分布を順次記憶する（図１４（Ｃ））。なお、ひずみ速度テンソル分布は、図示のように、所定範囲ごとに識別可能に断層画像の態様で分布表示することができる。

図１５は、ひずみ速度可視化画面を示す図である。（Ａ）は本実施形態によるひずみ速度可視化画面を示す。（Ｂ）は比較例として、ＦＥＭにより本実施形態と同じ条件にてひずみ速度解析をした結果を示す。各図の横軸はｘ方向の位置を、縦軸はｚ方向の位置をそれぞれ示している。

図１５（Ａ）および（Ｂ）から、ひずみ速度分布の演算および可視化処理に関し、本実施形態によるＯＣＴによる演算結果が、ＦＥＭ解析の結果と精度良く対応していることが確認できる。応力分布について高精度な結果が得られることは第１実施形態で述べたとおりであるから、これら応力分布とひずみ速度分布とに基づいて、力学物性値の断層分布を精度良く得ることが可能となる。

次に、制御演算部４が実行する具体的処理の流れについて説明する。
図１６は、制御演算部４により実行される力学物性値可視化処理の流れを示すフローチャートである。制御演算部４は、ＯＣＴによる測定対象Ｗの断層画像を連続的に取得する（Ｓ２１０）。そして、その断層画像を用いて応力の断層分布を演算し、画面に表示する（Ｓ２１２）。また、その断層画像を用いひずみ速度の断層分布を演算し、画面に表示する（Ｓ２１４）。そして、それら応力断層分布およびひずみ速度断層分布の演算結果から力学物性値の断層分布を演算し、画面に表示する（Ｓ２１６）。ここで、Ｓ２１２の応力断層分布の演算処理等については、第１実施形態の処理を適用することができる。このため、その説明については省略する。

図１７は、図１６におけるＳ２１４のひずみ速度分布の演算可視化処理を詳細に示すフローチャートである。制御演算部４は、連続撮影された時刻の異なる前後２枚のＯＣＴ断層画像Ｉ(x,z,t)とＩ(x,z,t+Δt)を取得し、再帰的相互相関法による処理を実行する（Ｓ２２０）。各ＯＣＴ断層画像は、図９のＳ１０で取得された水平偏光成分Ｉ_Ｈ(x,z)および垂直偏光成分Ｉ_Ｖ(x,z)の二乗和として得られる。すなわち、上記式（４）のストークスパラメータＳ_０が、そのＯＣＴ断層画像に対応する。

ここではまず、最小解像度（最大サイズの検査領域）での相互相関処理を実行し、相関係数分布を求める。続いて、隣接相互相関乗法により、隣接する相関係数分布の積を演算する（Ｓ２２２）。このとき、標準偏差フィルタ等の空間フィルタにより過誤ベクトルの除去をし（Ｓ２２４）、最小二乗法等により除去ベクトルの内挿補間を実行する（Ｓ２２６）。続いて、検査領域を小さくすることによって解像度を上げて相互相関処理を継続する（Ｓ２２８）。すなわち、低解像度での参照ベクトルを基に相互相関処理を実行する。このときの解像度が予め定める最高解像度でなければ（Ｓ２３０のＮ）、Ｓ２２２に戻る。

そして、Ｓ２２２〜Ｓ２２８の処理を繰り返し、最高解像度での相互相関処理が完了すると（Ｓ２３０のＹ）、サブピクセル解析を実行する。すなわち、最高解像度（最小サイズの検査領域）での変形ベクトルの分布に基づき、風上勾配法によるサブピクセル移動量を演算する（Ｓ２３２）。そして、このとき算出されたサブピクセル移動量に基づき、画像変形法によるサブピクセル変形量を演算する（Ｓ２３４）。続いて、最大相互相関値によるフィルタ処理により過誤ベクトルの除去をし（Ｓ２３６）、最小二乗法等により除去ベクトルの内挿補間を実行する（Ｓ２３８）。そして、このようにして得られた変形ベクトルの時間微分を行い、その断層画像について変形速度ベクトルの断層分布を算出する。そして、その変形速度ベクトルに対して空間微分をすることにより、ひずみ速度分布を算出し（Ｓ２４０）、表示装置に断層可視化する（Ｓ２４２）。

そして、図１６のＳ２１６において、応力断層分布およびひずみ速度断層分布の演算結果から力学物性値の断層分布を演算し、画面に表示する。すなわち、応力断層分布のデータを、ひずみ速度断層分布データ（もしくはひずみ断層分布データ）で除算することにより、弾性係数等の力学物性値の断層分布を算出し、可視表示させることができる。

ここでは、Ｓ２１２で得られる応力分布と、Ｓ２１４で得られるひずみ速度分布とが、いずれもＳ２１０で取得された同一の断層画像に基づくものとして対応づけられている。このため、Ｓ２１２の応力断層分布と、Ｓ２１４のひずみ速度断層分布とはほぼ同時に得られることとなる。そして、両断層分布の対応する断層位置について上述した除算処理がなされることにより、力学物性値の断層分布を算出することができる。

なお、本実施形態では、連続的に取得される前後の断層画像による演算を行ったが、変形例においては、所定の時間間隔で取得された前後の断層画像による演算処理を行ってもよい。例えば、時刻ｔ１で取得した断層画像と、時刻ｔ２で取得した断層画像とについて平均近似処理を行うなどして力学物性値の断層分布を算出してもよい。具体的には、時刻ｔ１，ｔ２のそれぞれ断層画像から各時刻の応力分布を演算し、平均近似処理などにより時刻ｔ＝（ｔ１＋ｔ２）／２の応力分布を取得してもよい。その一方、時刻ｔ１，ｔ２のそれぞれ断層画像から時刻ｔ＝（ｔ１＋ｔ２）／２のひずみ速度（ひずみ）分布を取得してもよい。そして、このように同時計測されたひずみと応力の両者の空間情報から、粘性係数や弾性係数等の物性パラメータを力学物性値の断層分布として算出し、可視表示してもよい。

以上に説明したように、本実施形態によれば、応力断層分布、ひずみ速度断層分布をマイクロ断層可視化でき、さらに両断層分布から力学物性値の断層分布をも断層可視化することができる。すなわち、高分子基材料の内部における応力・ひずみ分布を同時に断層検出し、さらに動的弾性率の断層可視化をも可能にするシステムを構築することができる。本実施形態によれば、測定対象の物性に変性をきたすことなく、その測定対象に負荷される応力の分布、ひずみの分布、さらには力学物性値の分布を光学的手法により簡易に測定および可視化可能な技術を提供できる。

以上、本発明の好適な実施形態について説明したが、本発明はその特定の実施形態に限定されるものではなく、本発明の技術思想の範囲内で種々の変形が可能であることはいうまでもない。

上記実施形態では、図９のＳ１０において、２軸の断層画像Ｉ_Ｈ(x,z)，Ｉ_Ｖ(x,z)を取得して処理を進め、Ｓ３４にて応力分布をｘ−ｚ平面に二次元的に可視化する例を示した。変形例においては、特定のｘ座標についての１軸断層画像Ｉ_Ｈ(z)を取得し、ｚ方向の応力分布として可視化してもよい。

上記第２実施形態の技術思想は、例えば以下のように表現できる。
測定対象の力学物性値の断層分布を可視化する力学物性値可視化装置であって、
光コヒーレンストモグラフィー（ＯＣＴ）を用いる光学系を含む光学ユニットと、
前記測定対象を変形させるための変形エネルギー（負荷）を付与可能な負荷装置と、
前記光学ユニットおよび前記負荷装置の駆動を制御し、それらの駆動に基づいて前記光学ユニットから出力された光干渉信号を処理し、所定の演算処理を行う制御演算部と、
前記制御演算部の演算結果を表示可能な表示装置と、
を備え、
前記制御演算部は、
前記光学ユニットの駆動に伴って得られた光干渉信号に基づいて、前記測定対象に作用又は内在する応力の断層分布を演算し、
前記光学ユニットおよび前記負荷装置の駆動に伴って得られた光干渉信号に基づいて、前記変形エネルギーによる前記測定対象の断層位置での変形速度ベクトル分布（又は変形ベクトル分布）を演算し、その変形速度ベクトル分布（又は変形ベクトル分布）に基づいてひずみ速度テンソル（又はひずみテンソル）の断層分布を演算し、
前記応力の断層分布の演算結果と、前記ひずみ速度テンソル（又はひずみテンソル）の断層分布の演算結果とに基づいて、前記力学物性値の断層分布を演算し、前記表示装置に可視化表示させることを特徴とする力学物性値可視化装置。

ここで、「負荷装置」は、測定対象に所定の荷重を負荷可能な荷重機構であってもよいし、超音波、光音響波、電磁波等による負荷（加振力）を測定対象に付与する装置であってもよい。制御演算部は、第１実施形態に示した位相シフト法を用いて応力断層分布を演算してもよいし、位相シフト法を用いずにそれを演算してもよい。後者の場合、応力分布の演算に際して前記負荷装置の駆動を必須としてもよい。ＯＣＴは、ＰＳ−ＯＣＴであってもよいし、他のＯＣＴであってもよい。また、測定対象は高分子基材料に限らず、その他の材料構造を対象としてもよい。また、体内の軟骨やしこり（硬化箇所）など、生体構造（生体組織）を対象としてもよい。生体構造の力学特性の変化を断層可視化することで、例えば動脈硬化，癌，肝硬変，皮膚，軟骨，骨などの多くの症例に対する診断材料とすることが可能となる。

なお、本発明は上記実施形態や変形例に限定されるものではなく、要旨を逸脱しない範囲で構成要素を変形して具体化することができる。上記実施形態や変形例に開示されている複数の構成要素を適宜組み合わせることにより種々の発明を形成してもよい。また、上記実施形態や変形例に示される全構成要素からいくつかの構成要素を削除してもよい。

Claims (5)

1. 測定対象に作用又は内在する応力の分布を断層可視化する応力可視化装置であって、
   低コヒーレンス光又は波長走査した光を出射する光源と、
   前記光源から出射された光を直線偏光する偏光子と、
   直線偏光された光を前記測定対象を経由するオブジェクトアームと、参照鏡を経由するリファレンスアームとに分けるビームスプリッタと、
   前記オブジェクトアームに設けられ、前記直線偏光された光の偏光状態を変化させることが可能な光学変調器と、
   前記測定対象にて反射した物体光と前記参照鏡にて反射した参照光とが重畳された干渉光を、その水平偏光成分と垂直偏光成分とに分けて検出する光検出装置と、
   前記光学変調器を制御することにより前記偏光状態を調整する一方、前記測定対象の奥行き方向の位置とその直交方向の位置とにより特定される断層位置における応力分布を、前記光検出装置が検出した干渉信号の水平偏光成分と垂直偏光成分との位相差に基づいて演算する制御演算部と、
   前記制御演算部の演算結果に基づいて、前記測定対象における応力分布を断層可視化する態様で表示する表示装置と、
   を備え、
   前記制御演算部は、前記光学変調器による偏光状態を変化させることにより前記測定対象に照射する偏光波の位相をシフトさせ、その偏光波の位相シフトごとに得られる干渉信号の位相差に基づき、その干渉信号の位相差断層分布の空間勾配を演算し、その空間勾配の断層分布を応力の断層分布と対応づけるようにして前記表示装置に可視化表示させることを特徴とする応力可視化装置。
2. 前記制御演算部は、
   前記測定対象となりうる部材の材質ごとに、前記干渉信号として得られるべき位相差の空間勾配と応力値との対応関係が予め定められた演算マップを保持し、
   前記干渉信号の位相差断層分布の空間勾配を演算した後、その演算結果に基づいて前記演算マップを参照することにより応力を算出し、その算出結果としての応力分布を断層可視化することを特徴とする請求項１に記載の応力可視化装置。
3. 前記測定対象に外的負荷をかけることなく、その測定対象に内在する残留応力を可視化表示することを特徴とする請求項１または２に記載の応力可視化装置。
4. 高分子基材料を前記測定対象とすることを特徴とする請求項１〜３のいずれかに記載の応力可視化装置。
5. 測定対象の力学物性値の断層分布を可視化する力学物性値可視化装置であって、
   光コヒーレンストモグラフィーを用いる光学系を含む光学ユニットと、
   前記測定対象を変形させるための変形エネルギーを付与可能な負荷装置と、
   前記光学ユニットおよび前記負荷装置の駆動を制御し、それらの駆動に基づいて前記光学ユニットから出力された光干渉信号を処理し、所定の演算処理を行う制御演算部と、
   前記制御演算部の演算結果を表示可能な表示装置と、
   を備え、
   前記制御演算部は、
   前記光学ユニットの駆動に伴って得られた光干渉信号に基づいて、前記測定対象に作用又は内在する応力の断層分布を演算し、
   前記光学ユニットおよび前記負荷装置の駆動に伴って得られた光干渉信号に基づいて、前記変形エネルギーによる前記測定対象の断層位置での変形速度ベクトル分布又は変形ベクトル分布を演算し、その変形速度ベクトル分布又は変形ベクトル分布に基づいてひずみ速度テンソル又はひずみテンソルの断層分布を演算し、
   前記応力の断層分布の演算結果と、前記ひずみ速度テンソル又はひずみテンソルの断層分布の演算結果とに基づいて、前記力学物性値の断層分布を演算し、前記表示装置に可視化表示させることを特徴とする力学物性値可視化装置。

Images