TOP > 国内特許検索 > 応力可視化装置および力学物性値可視化装置 > 明細書

明細書 :応力可視化装置および力学物性値可視化装置

発行国 日本国特許庁(JP)
公報種別 再公表特許(A1)
発行日 平成29年6月1日(2017.6.1)
発明の名称または考案の名称 応力可視化装置および力学物性値可視化装置
国際特許分類 G01L   1/24        (2006.01)
G01B  11/16        (2006.01)
FI G01L 1/24 Z
G01B 11/16 G
国際予備審査の請求
全頁数 32
出願番号 特願2016-544258 (P2016-544258)
国際出願番号 PCT/JP2015/073457
国際公開番号 WO2016/027874
国際出願日 平成27年8月21日(2015.8.21)
国際公開日 平成28年2月25日(2016.2.25)
優先権出願番号 2014168802
優先日 平成26年8月21日(2014.8.21)
優先権主張国 日本国(JP)
指定国 AP(BW , GH , GM , KE , LR , LS , MW , MZ , NA , RW , SD , SL , ST , SZ , TZ , UG , ZM , ZW) , EA(AM , AZ , BY , KG , KZ , RU , TJ , TM) , EP(AL , AT , BE , BG , CH , CY , CZ , DE , DK , EE , ES , FI , FR , GB , GR , HR , HU , IE , IS , IT , LT , LU , LV , MC , MK , MT , NL , NO , PL , PT , RO , RS , SE , SI , SK , SM , TR) , OA(BF , BJ , CF , CG , CI , CM , GA , GN , GQ , GW , KM , ML , MR , NE , SN , TD , TG) , AE , AG , AL , AM , AO , AT , AU , AZ , BA , BB , BG , BH , BN , BR , BW , BY , BZ , CA , CH , CL , CN , CO , CR , CU , CZ , DE , DK , DM , DO , DZ , EC , EE , EG , ES , FI , GB , GD , GE , GH , GM , GT , HN , HR , HU , ID , IL , IN , IR , IS , JP , KE , KG , KN , KP , KR , KZ , LA , LC , LK , LR , LS , LU , LY , MA , MD , ME , MG , MK , MN , MW , MX , MY , MZ , NA , NG , NI , NO , NZ , OM , PA , PE , PG , PH , PL , PT , QA , RO , RS , RU , RW , SA , SC , SD , SE , SG , SK , SL , SM , ST , SV , SY , TH , TJ , TM , TN , TR , TT , TZ , UA , UG , US
発明者または考案者 【氏名】佐伯 壮一
出願人 【識別番号】506122327
【氏名又は名称】公立大学法人大阪市立大学
個別代理人の代理人 【識別番号】110002273、【氏名又は名称】特許業務法人インターブレイン
審査請求 未請求
テーマコード 2F065
Fターム 2F065AA65
2F065BB22
2F065FF49
2F065FF52
2F065FF61
2F065GG07
2F065HH13
2F065JJ01
2F065JJ05
2F065JJ09
2F065JJ18
2F065LL04
2F065LL12
2F065LL53
2F065LL62
2F065MM16
2F065QQ31
2F065QQ41
要約 制御演算部114は、光学変調器110による偏光状態を変化させることにより測定対象Wに照射する偏光波の位相をシフトさせ、その偏光波の位相シフトごとに得られる干渉信号の位相差に基づき、その干渉信号の位相差断層分布の空間勾配を演算し、その空間勾配の断層分布を応力の断層分布と対応づけるようにして表示装置116に可視化表示させる。制御演算部114は、光干渉信号に基づいて測定対象の断層位置での変形速度ベクトル分布を演算し、さらにひずみ速度テンソルの断層分布を演算する。そして、応力の断層分布とひずみ速度テンソルの断層分布とに基づいて力学物性値の断層分布を演算し、可視化表示させる。
特許請求の範囲 【請求項1】
測定対象に作用又は内在する応力の分布を断層可視化する応力可視化装置であって、
低コヒーレンス光又は波長走査した光を出射する光源と、
前記光源から出射された光を直線偏光する偏光子と、
直線偏光された光を前記測定対象を経由するオブジェクトアームと、参照鏡を経由するリファレンスアームとに分けるビームスプリッタと、
前記オブジェクトアームに設けられ、前記直線偏光された光の偏光状態を変化させることが可能な光学変調器と、
前記測定対象にて反射した物体光と前記参照鏡にて反射した参照光とが重畳された干渉光を、その水平偏光成分と垂直偏光成分とに分けて検出する光検出装置と、
前記光学変調器を制御することにより前記偏光状態を調整する一方、前記測定対象の奥行き方向の位置とその直交方向の位置とにより特定される断層位置における応力分布を、前記光検出装置が検出した干渉信号の水平偏光成分と垂直偏光成分との位相差に基づいて演算する制御演算部と、
前記制御演算部の演算結果に基づいて、前記測定対象における応力分布を断層可視化する態様で表示する表示装置と、
を備え、
前記制御演算部は、前記光学変調器による偏光状態を変化させることにより前記測定対象に照射する偏光波の位相をシフトさせ、その偏光波の位相シフトごとに得られる干渉信号の位相差に基づき、その干渉信号の位相差断層分布の空間勾配を演算し、その空間勾配の断層分布を応力の断層分布と対応づけるようにして前記表示装置に可視化表示させることを特徴とする応力可視化装置。
【請求項2】
前記制御演算部は、
前記測定対象となりうる部材の材質ごとに、前記干渉信号として得られるべき位相差の空間勾配と応力値との対応関係が予め定められた演算マップを保持し、
前記干渉信号の位相差断層分布の空間勾配を演算した後、その演算結果に基づいて前記演算マップを参照することにより応力を算出し、その算出結果としての応力分布を断層可視化することを特徴とする請求項1に記載の応力可視化装置。
【請求項3】
前記測定対象に外的負荷をかけることなく、その測定対象に内在する残留応力を可視化表示することを特徴とする請求項1または2に記載の応力可視化装置。
【請求項4】
高分子基材料を前記測定対象とすることを特徴とする請求項1~3のいずれかに記載の応力可視化装置。
【請求項5】
測定対象の力学物性値の断層分布を可視化する力学物性値可視化装置であって、
光コヒーレンストモグラフィーを用いる光学系を含む光学ユニットと、
前記測定対象を変形させるための変形エネルギーを付与可能な負荷装置と、
前記光学ユニットおよび前記負荷装置の駆動を制御し、それらの駆動に基づいて前記光学ユニットから出力された光干渉信号を処理し、所定の演算処理を行う制御演算部と、
前記制御演算部の演算結果を表示可能な表示装置と、
を備え、
前記制御演算部は、
前記光学ユニットの駆動に伴って得られた光干渉信号に基づいて、前記測定対象に作用又は内在する応力の断層分布を演算し、
前記光学ユニットおよび前記負荷装置の駆動に伴って得られた光干渉信号に基づいて、前記変形エネルギーによる前記測定対象の断層位置での変形速度ベクトル分布又は変形ベクトル分布を演算し、その変形速度ベクトル分布又は変形ベクトル分布に基づいてひずみ速度テンソル又はひずみテンソルの断層分布を演算し、
前記応力の断層分布の演算結果と、前記ひずみ速度テンソル又はひずみテンソルの断層分布の演算結果とに基づいて、前記力学物性値の断層分布を演算し、前記表示装置に可視化表示させることを特徴とする力学物性値可視化装置。
発明の詳細な説明 【技術分野】
【0001】
本発明は応力可視化装置に関し、特に測定対象に負荷される応力又は内在する応力の分布を光学的手法により測定して可視化する装置に関する。本発明はまた、測定対象の力学物性値の断層分布を可視化する力学物性値可視化装置に関する。
【背景技術】
【0002】
近年、航空・宇宙分野をはじめ、自動車,電子実装部品,医療器具など様々な分野において複合材料が急速に普及し、その形状やサイズの多様化および微細化が進んでいる。このような複合材料は、母材と強化材の積層面を有することから、界面剥離やトランスバースクラックなどマイクロスケールの内部欠陥が複雑に関係して破壊に至ることがある。このため、材料内部における応力集中などの力学的状態を定量評価可能な非破壊計測手法の開発が望まれている。
【0003】
このような計測手法として、光学的手法により検査対象表面の残留応力を非破壊非接触にて測定する方法が提案されている(例えば特許文献1参照)。この手法では、検査対象表面に加熱用のレーザ光と変形量測定用のレーザ光とを照射し、検査領域における熱応力開放前後の全変形量を光干渉計測により取得し、その全変形量から熱ひずみによる変形量を差し引くことにより残留応力による変形量を算出する。その残留応力による変形量を加熱前後のヤング率の差により除算することにより、残留応力を算出するものである。
【先行技術文献】
【0004】

【特許文献1】特開2008-134186号公報
【発明の概要】
【発明が解決しようとする課題】
【0005】
しかしながら、特許文献1に記載の技術では、検査対象表面の応力分布が分かるものの、その検査対象内部の応力分布を知ることはできない。また、応力測定のために干渉計測装置に加えて加熱用装置が必要となるため、計測システム全体として大がかりとなる。さらに、検査対象に対して外的に熱応力を負荷するため、仮にそれを弾性変形範囲内に抑えるとしても、検査対象の物性に変性をきたす可能性がある点で問題があった。
【0006】
本発明はこのような課題に鑑みてなされたものであり、その目的の一つは、測定対象に負荷される応力又は内在する応力の分布を光学的手法により簡易に測定および可視化可能な技術を提供することにある。
【課題を解決するための手段】
【0007】
本発明のある態様は、測定対象に作用又は内在する応力の分布を断層可視化する応力可視化装置である。この応力可視化装置は、低コヒーレンス光又は波長走査した光を出射する光源と、光源から出射された光を直線偏光する偏光子と、直線偏光された光を測定対象を経由するオブジェクトアームと、参照鏡を経由するリファレンスアームとに分けるビームスプリッタと、オブジェクトアームに設けられ、直線偏光された光の偏光状態を変化させることが可能な光学変調器と、測定対象にて反射した物体光と参照鏡にて反射した参照光とが重畳された干渉光を、その水平偏光成分と垂直偏光成分とに分けて検出する光検出装置と、光学変調器を制御することにより偏光状態を調整する一方、測定対象の奥行き方向の位置とその直交方向の位置とにより特定される断層位置における応力分布を、光検出装置が検出した干渉信号の水平偏光成分と垂直偏光成分との位相差に基づいて演算する制御演算部と、制御演算部の演算結果に基づいて、測定対象における応力分布を断層可視化する態様で表示する表示装置と、を備える。
【0008】
制御演算部は、光学変調器による偏光状態を変化させることにより測定対象に照射する偏光波の位相をシフトさせ、その偏光波の位相シフトごとに得られる干渉信号の位相差に基づき、その干渉信号の位相差断層分布の空間勾配を演算し、その空間勾配の断層分布を応力の断層分布と対応づけるようにして表示装置に可視化表示させる。
【0009】
この態様によると、測定対象の応力分布の計測のために、いわゆる偏光感受型光干渉断層計(Polarization Sensitive Optical Coherence Tomography:以下「PS-OCT」という)が利用される。それにより、測定対象の内部における2偏波位相差分布を断層計測でき、その断層画像に表れるべき干渉縞の空間周波数から応力の分布を算出し、可視化することができる。すなわち、測定対象に対して応力測定のための別途の応力を負荷する必要がなく、応力断層分布を簡易に測定することができる。そして特に、偏光波の位相をシフトさせることにより、測定対象の内部の任意位置での応力計測をOCTの空間分解能にて実施できるようになる。
【発明の効果】
【0010】
本発明によれば、測定対象に負荷される応力又は残留する応力の分布を光学的手法により簡易に測定可能な技術を提供することができる。
【図面の簡単な説明】
【0011】
【図1】第1実施形態に係る応力可視化装置の構成を表す概略図である。
【図2】PS-OCTによる応力断層計測の原理を表す説明図である。
【図3】PS-OCTによる応力断層計測の原理を表す説明図である。
【図4】PS-OCTによる応力断層計測の原理を表す説明図である。
【図5】光学変調器への印加電圧を変化させて位相シフトを行った場合のS3/S0の変化を演算した結果を例示する図である。
【図6】光学変調器への印加電圧を変化させて位相シフトを行った場合の特定の計測位置(x,z)におけるS3/S0の変化を示す図である。
【図7】PS-OCTを用いた応力計測に用いられる演算マップを示す図である。
【図8】実施例に係る応力可視化装置の全体構成を表す図である。
【図9】制御演算部により実行される応力可視化処理の流れを概略的に示す図である。
【図10】実験に使用された測定対象Wおよびその応力負荷方向を示す図である。
【図11】応力可視化状態を示す画面例を示す図である。
【図12】再帰的相互相関法による処理手順を概略的に示す図である。
【図13】サブピクセル解析による処理手順を概略的に示す図である。
【図14】力学物性値の断層可視化処理の手順を表す図である。
【図15】ひずみ速度可視化画面を示す図である。
【図16】制御演算部4により実行される物理量断層可視化処理の流れを示すフローチャートである。
【図17】図16におけるS214のひずみ速度分布の演算可視化処理を詳細に示すフローチャートである。
【発明を実施するための形態】
【0012】
[第1実施形態]
以下、本実施形態に係る応力可視化方法の原理について説明する。
図1は、第1実施形態に係る応力可視化装置の構成を表す概略図である。本実施形態の応力可視化装置は、2偏波を扱うPS-OCTを利用して測定対象Wの応力断層計測を行うものであり、光源100、偏光子102、オブジェクトアーム104、リファレンスアーム106、ビームスプリッタ108、電気光学変調器(Electro-Optic Modulator:以下「EOM」と表記する)110、光検出装置112、制御演算部114および表示装置116を備える。なお、図示の例では、マイケルソン干渉計をベースとした光学系が示されているが、マッハツェンダー干渉計その他の光学系を採用することもできる。

【0013】
光源100には、低コヒーレンス光(時間的に低コヒーレンスな光)を出射可能なものが採用される。なお、図示の例では、光遅延走査方式の分類におけるTD-OCT(Time Domain OCT)が示されているが、SD-OCT(Spectral Domain OCT)、SS-OCT(Swept Source OCT)その他のOCTを用いることもできる。

【0014】
ここで、「TD-OCT」は、リファレンスアーム106において参照鏡118を駆動して光遅延走査を行うことで、オブジェクトアーム104の後方散乱光との干渉光を計測し、測定対象Wの深さ方向の反射光強度を得るものである。「SD-OCT」は、TD-OCTのような機械的な光遅延走査を行うことなく、分光器にて干渉光を分光し、フーリエ変換を行うことにより反射光強度を得るものである。「SS-OCT」は、同様に機械的な光遅延走査を行うことなく、擬似的に低コヒーレンス光(広帯域光源)とするとともに、フーリエ変換により反射光強度を得るものである。

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

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

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

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

【0019】
制御演算部114は、測定対象Wのx-z平面上で特定される断層位置に応じた応力分布を、光検出装置112の検出値に基づいて演算する。また、この演算により計測される断層画像の解像度を高めるために、測定対象Wに照射される光の偏光状態(進相軸と遅相軸との位相差)を適宜シフトさせる位相シフト法が採用される。この位相シフトは、EOM110を制御することにより行われる。

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

【0021】
以下、このような演算処理の原理について順に説明する。図2~図4は、PS-OCTによる応力断層計測の原理を表す説明図である。本実施形態の断層計測は、3次元的に応力分布を計測する光弾性散乱光法にPS-OCTを適用することにより行われる。光弾性散乱光法は、測定対象Wに偏光波を照射したときの後方散乱光を用いてその測定対象Wの内部の光波の位相差を計測することにより、3次元的に応力分布を計測する手法である。

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

【0023】
このため、測定対象Wに作用する主応力が変化した場合、あるいは局所的に主応力の変化が内在している場合には、それが干渉光の水平偏光成分Eと垂直偏光成分Eとの位相差Φ(x,y,z)として表れる。すなわち、その2偏波の位相差Φ(x,y,z)にしたがってその2偏波に対応する干渉信号が得られる。この位相差Φ(x,y,z)は、下記式(1)として表すことができる。
【数1】
JP2016027874A1_000003t.gif

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

【0025】
一方、本実施形態のPS-OCTにおいては、干渉光強度の水平偏光成分I(x,z)と垂直偏光成分I(x,z)とについて、下記式(2)および(3)が成立する。
【数2】
JP2016027874A1_000004t.gif

【0026】
ここで、R(x,z)は反射率分布、Φ(x,z)は位相差Φ(x,y,z)のx,y成分、lは出射光のコヒーレンス長である。Δlは、オブジェクトアーム104とリファレンスアーム106との光路長差である。

【0027】
このとき、偏光状態を定式化したストークスパラメータの成分のうち、SとSを下記式(4)および(5)に基づいて算出することができる。ここで、Sは計測される光全体の強度を示し、Sは直交成分の直線偏光の強度差を示すが、このことは周知であるため、詳細な説明は省略する。
【数3】
JP2016027874A1_000005t.gif

【0028】
ここで、S(x,z)をS(x,z)で除算することにより、下記式(6)に示すように、反射強度の情報を除去し、力学複屈折効果に依存する偏光波の位相差Φ(x,z)のみに依存する余弦関数を抽出することができる。
【数4】
JP2016027874A1_000006t.gif

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

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

【0031】
なお、本実施形態ではビーム走査をx方向に1次元的に行うため、2次元断層計測に基づく2次元応力分布を可視化することになるが、ビームをx-y平面にて2次元的に走査させれば3次元断層計測を行うことができ、3次元応力分布を可視化することができる。

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

【0033】
この位相シフト処理においては、ストークスパラメータマップS/S(x,z)における干渉縞の位相を変化させた断層像を撮影する。その際、EOM110を用いて円偏光から位相差δEOM(V)を変化させた光を入射させる。このとき、取得されるストークスパラメータマップS/S(x,z)は、上記式(6)を変形して下記式(7)のようになる。
【数5】
JP2016027874A1_000007t.gif

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

【0035】
図5は、EOM110への印加電圧を変化させて位相シフトを行った場合のS/Sの変化を演算した結果を例示する図である。(A)は印加電圧を50Vとした場合、(B)は100Vとした場合、(C)は150Vとした場合を示す。各図の横軸は測定対象Wにおけるz軸方向の深さを示し、縦軸はS/S(x,z)の値を示す。図中の実線は各印加電圧を示す。図中の波線は参照電圧として0V、つまりEOM110を動作させなかった場合の演算結果を示している。この演算においては、ストークスパラメータマップS/S(x,z)をx軸方向に平均化し、z軸方向に関する位相差変化が算出されている。

【0036】
ここでは、測定対象Wに対してx軸方向に4.5Mpaの引張応力を外的に負荷している。また、印加電圧が0Vのときに、測定対象Wへの入射光の偏光状態が円偏光になるように偏光状態を設定した。図示の結果より、EOM110への印加電圧を変化させることにより、ストークスパラメータマップS/S(x,z)の概形を維持したまま、各計測位置(x,z)での偏光波の位相をシフトできることが分かる。

【0037】
図6は、EOM110への印加電圧を変化させて位相シフトを行った場合の特定の計測位置(x,z)におけるS/Sの変化を示す図である。同図の横軸はEOM110への印加電圧を示し、縦軸はS/S(x,z)の値を示す。この演算結果は、上記式(7)において、座標(x,z)を特定座標に固定したときの結果を示す。

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

【0039】
図7は、PS-OCTを用いた応力計測に用いられる演算マップを示す図である。(A)はある材料に適用される演算マップを示す。同図の横軸はS/S(x,z)の計測から得られる空間周波数(1/μm)を示し、縦軸はその計測位置における主応力差σ'(x,z)(MPa)を示す。(B)は同材料について公知の応力-ひずみ線図である。同図の横軸はひずみ(%)を示し、縦軸は応力(MPa)を示す。

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

【0041】
以下、図面を参照しつつ本実施形態を具体化した実施例について詳細に説明する。
図8は、実施例に係る応力可視化装置の全体構成を表す図である。本実施例の応力可視化装置1は、マイケルソン干渉光学系をベースとしたTD-PS-OCTを利用して測定対象の応力分布を演算するシステムとして構成されている。応力可視化装置1は、OCTによる光学的処理を実行する光学系2と、その処理結果に基づいて測定対象Wの応力分布を演算する制御演算部4とを備える。光学系2は、光源6、偏光子8、カプラー10、リファレンスアーム12、オブジェクトアーム14および光検出装置16等の光学要素を含む。各光学要素は、ファイバーにて互いに接続されている。

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

【0043】
偏光子8は、光源6からの光を直線偏光にする直線偏光子(Linear Polarizer:以下「LP」と表記する)であり、光源6とカプラー10との間に配置されている。カプラー10は、光源6からの光をリファレンスアーム12とオブジェクトアーム14とに向けて分波するビームスプリッタとしての機能と、リファレンスアーム12およびオブジェクトアーム14のそれぞれから戻ってきた光を合波して干渉させる機能とを有する。

【0044】
リファレンスアーム12には、偏波コントローラ(Polarization Controller:以下「PC」と表記する)18、LP26、EOM28、集光レンズ22が順次接続されている。集光レンズ22は、参照鏡24に対向配置されている。LP26は、カプラー10を経てリファレンスアーム12に導かれた光の偏光状態を整える。EOM28は、印加電圧を線形変化させることで、単一の高周波なキャリア周波数を得ることを可能にするものであり、EOPM(Electro‐Optic Phase Modulator)とも呼ばれる(以下「EOPM28」とも表記する)。これにより、OCTによる干渉信号を高感度にてヘテロダイン検出することができる。

【0045】
LP26は、EOPM28により偏光状態を変調しても、リファレンスアーム12から戻ってくる光の強度が変化しないようPC18と協働する。具体的には、PC18は、参照鏡24からの参照光が再びカプラー10に入射する際の偏光状態を45度の直線偏光となるように調整する。

【0046】
オブジェクトアーム14には、PC30、1/4波長板(Quarter Wave Plate:以下「QWP」と表記する)44、EOM46、コリメートレンズ34、ガルバノ装置36、対物レンズ38が順次接続されている。ガルバノ装置36は、固定鏡40およびガルバノ鏡42を含む。対物レンズ38は、測定対象Wに対向配置される。カプラー10を介してオブジェクトアーム14に入射された光は、2軸のガルバノ鏡42によりx軸方向に走査されて測定対象Wに照射される。測定対象Wからの反射光は、物体光として再びカプラー10に戻り、参照光と重畳されて干渉光として光検出装置16に送られる。

【0047】
PC30は、カプラー10からオブジェクトアーム14に導かれた光の偏光状態を円偏光や楕円偏光となり得るように調整する。QWP44は、45度の直線偏光を通過させることでその進相軸と遅相軸との位相差を1/4波長ずらし、円偏光に変化させる。EOM46は、印加電圧に応じて透過光の偏光状態を任意に変化させるものであり、EOAM(Electro‐Optic Amplitude Modulator)とも呼ばれる。これにより、QWP44を通過した円偏光を楕円偏光に変化させることができるとともに、印加電圧に応じてその楕円偏光の進相軸と遅相軸との位相差を任意に変化させることができる。

【0048】
光検出装置16は、偏光ビームスプリッタ(Polarization Beam Splitter:以下「PBS」と表記する)50、光検出器52(「第1光検出器」として機能する)、光検出器54(「第2光検出器」として機能する)を含む。PBS50は、カプラー10にて合波された干渉光を水平偏光成分Eと垂直偏光成分Eとに分波し、それぞれを光検出器52,54へ出力する。そして、光検出器52が干渉光強度の水平偏光成分I(x,z)を検出し、光検出器54が干渉光強度の垂直偏光成分I(x,z)を検出する。これらの検出は同時に行われる。その水平偏光成分I(x,z)は、バンドパスフィルタ(Band-Pass Filter:以下「BPF」と表記する)56、アンプ58(Lock-in Amp)およびA/D変換器60を介して制御演算部4に入力される。一方、その垂直偏光成分I(x,z)は、BPF62、アンプ64(Lock-in Amp)およびA/D変換器60を介して制御演算部4に入力される。

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

【0050】
制御演算部4は、上述した原理に基づき、PS-OCTによる計測結果に対して位相シフト演算を行うことによりストークスパラメータマップS/S(x,z)の空間周波数を算出する。そして、算出された空間周波数を用いて図7に示した演算マップを参照することにより、測定対象Wの任意の位置に作用する応力を算出し、表示装置72に表示させる。

【0051】
図9は、制御演算部4により実行される応力可視化処理の流れを概略的に示す図である。この応力可視化処理において、制御演算部4は、まず、EOM46への印加電圧を予め定める初期値(例えば0Vの非通電状態)に設定したうえで、PS-OCTによる断層画像を取得する。すなわち、干渉信号(干渉光強度)の水平偏光成分I(x,z)と垂直偏光成分I(x,z)とを取得する(S10)。続いて、EOM46の印加電圧を変化させて偏光波の位相をシフトさせ、断層画像を取得する(S12)。この位相シフトによる断層画像の取得を設定回数実行する(S14のN)。

【0052】
それにより設定数の断層画像が取得されると(S14のY)、各断層画像について前処理としてアンサンブル平均を施した後(S16)、移動最小二乗法によるノイズ除去処理と平滑化処理を実行する(S18)。そして、各断層画像についてストークスパラメータマップS/S(x,z)を算出する(S20)。

【0053】
そして、各断層画像のストークスパラメータマップS/S(x,z)について、以下の処理を実行する。すなわち、S/S(x,z)の1軸断層(OCTにより走査されたx座標の一つの断層)に対してフーリエ解析もしくは最小二乗法を適用し、そのz方向(奥行き方向)の位相Φ(z)を算出する(S22)。この処理をx方向(OCTによる走査方向)の各x座標について行うことにより、位相断層分布Φ(x,z)を算出する(S24)。このような処理により、上記式(7)におけるS/S(x,z)の各係数が求まる。このとき、位相接続処理を行っておく(S26)。すなわち、位相が0~2πで変化を繰り返すため、光軸方向に向かって位相が正に変化していくように位相接続をする。

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

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

【0056】
次に、本実施形態の応力断層可視化処理に基づく実験結果を示す。図10は、実験に使用された測定対象Wおよびその応力負荷方向を示す。(A)は測定対象Wとしての試験片を示し、(B)は荷重の負荷方法を示す。図11は、応力可視化画面を示す図である。(A)は本実施形態による応力可視化画面を示す。(B)は比較例として、FEMにより本実施形態と同じ条件にて応力解析をした結果を示す。(C)は、(A)の画像を得る前に取得したPS-OCTによる断層画像を示す。各図の横軸はx方向の位置を、縦軸はz方向の位置をそれぞれ示している。

【0057】
本実験においては、試験片W1として図10(A)に示すようなダンベル片を用いた。この試験片W1は、全長149mm、最大幅20mm、厚み1.5mmを有する低密度ポリエチレンからなる。図10(B)に示すように、PS-OCTによる偏光波の入射方向にz軸をとり、それと直交するようにx軸およびy軸を設定した。試験片W1の厚み方向(奥行き方向)の端面がx-y平面に延在するように設置される。また、偏光波の光軸と直交する走査方向をx軸方向とし、試験片W1の長手方向に一致させている。応力集中が発生するよう、試験片W1の長手方向中央の下面に深さ200μmの切欠きCを形成した。切欠きCは、y軸方向に延びる直線状に形成されている。この試験片W1に対して長手方向の引張荷重Pを与え、30分間応力緩和させた後、上述したPS-OCTによる計測を行い、応力断層可視化を行った。

【0058】
その結果、図11(A)に示すように、応力分布の断層画像を可視化することができた。図中右欄に応力値が示されているが、図11(B)に示すFEM解析の結果と精度良く対応していることが確認できる。なお、図11(A)では応力値そのものを可視表示しているが、図11(C)に示すように、S/S(x,z)の断層画像を表示させることにより、断層位置における干渉縞の空間周波数を把握することができ、それにより、応力の高低の分布を把握することができる。

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

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

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

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

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

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

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

【0066】
(再帰的相互相関法)
図12(A)~(C)は、再帰的相互相関法による処理過程を示している。各図にはOCTにより連続的に撮影される前後の断層画像が示されている。左側には先の断層画像(Image1)が示され、右側には後の断層画像(Image2)が示されている。

【0067】
相互相関法とは、局所的なスペックル・パターンの類似度を下記式(8)に基づく相関値Rijより評価する方法である。そのため、図12(A)に示すように、連続的に撮影された前後のOCT画像について、先の断層画像(Image1)に類似度の検査対象となる検査領域S1が設定され、後の断層画像(Image2)に類似度の探査範囲となる探査領域S2が設定される。
【数6】
JP2016027874A1_000008t.gif
ここで、空間座標として、光軸方向にZ軸、光軸と垂直方向にX軸を設定している。f(X,Z)とg(X,Z)は、変形前後のOCT画像において設置された中心位置(X,Z)の検査領域S1(N×Nピクセル)のスペックル・パターンを表している。

【0068】
探査領域S2(M×Mピクセル)内における相関値分布Ri,j(ΔX,ΔZ)を算出し、図12(B)に示すようにパターンマッチングを行う。実際には、下記式(9)に示すように、最大相関値を与える移動量Ui,jを変形前後の変形ベクトルとして決定する。
【数7】
JP2016027874A1_000009t.gif
なお、fとgはf(X,Z)とg(X,Z)の検査領域S1内の平均値を表している。

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

【0070】
また、これに加え、下記式(10)により、演算中の座標を中心とする周囲8つの座標を含む合計9つの変形ベクトルの平均偏差σを用いた閾値を設定し、過誤ベクトルの除去を行い再帰処理に伴う誤差伝播を抑制する。
【数8】
JP2016027874A1_000010t.gif
ここで、Umはベクトル量の中央値を表し、閾値となる係数κは任意に設定した。

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

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

【0073】
(風上勾配法)
図13(A)~(C)は、サブピクセル解析による処理過程を示している。各図にはOCTにより連続的に撮影される前後の断層画像が示されている。左側には先の断層画像(Image1)が示され、右側には後の断層画像(Image2)が示されている。

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

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

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

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

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

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

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

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

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

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

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

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

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

【0087】
伸縮及びせん断変形を考慮して算出した検査領域S1は、図13(C)に示すように、移動とともに変形を伴う。材料変形前のOCT断層像におけるある検査領域S1内の整数ピクセル位置での座標(x,z)が変形後に座標(x,z)に移動すると考えると、x,zの値は下記式(15)にて表される。
【数13】
JP2016027874A1_000015t.gif

【0088】
ここで、Δx,Δzはそれぞれ検査領域S1中心から座標x,zまでの距離、u,vはそれぞれx,z方向への変形量、∂u/∂x,∂v/∂zはそれぞれx,z方向における検査領域S1の垂直方向への変形量、∂u/∂z,∂v/∂xはそれぞれx,z方向における検査領域S1のせん断方向への変形量である。数値解法にはNewton-Raphson法を用い、6変数(u,v,∂u/∂x,∂u/∂z,∂v/∂x,∂v/∂z)での相関値微係数が0となるように、すなわち最大相関値を得るように反復計算を行う。なお、反復計算の収束性を高めるため、x,z方向の移動量初期値には風上勾配法で得られたサブピクセル移動量を用いる。相関値Rに対するヘッセ行列をH、相関値対するヤコビベクトルを▽Rとすると、1回の反復で得られる更新量ΔPiは下記式(16)にて表される。
【数14】
JP2016027874A1_000016t.gif

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

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

【0091】
(時空間移動最小二乗法)
ひずみ速度テンソルの算出には移動最小二乗法(Moving Least Square Method:以下「MLSM」という)を用いる。MLSMは、移動量分布を平滑化すると共に、微係数の安定算出を可能とする手法である。MLSMにおいて用いる二乗誤差式は下記式(17)にて表される。
【数15】
JP2016027874A1_000017t.gif

【0092】
上記式(17)において、S(x,z,t)を最小とするa~kのパラメータを求める。すなわち、近似関数として水平方向x,奥行き方向z,時間方向tの3変数2次多項式として下記式(18)を採用する。そして、最小二乗近似に基づき、下記式(19)から最適な微係数を算出して平滑化する。
【数16】
JP2016027874A1_000018t.gif

【0093】
この微係数を用いることにより、下記式(20)に示すひずみ速度テンソルを算出することができる。fx,fzは各軸のひずみ増分を示し、その時間変化量からひずみ速度を算出している。
【数17】
JP2016027874A1_000019t.gif

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

【0095】
制御演算部4は、測定対象Wの断層画像を連続的に取得し(図14(A),(B))、上述のようにして演算された変形速度ベクトル分布を順次記憶するとともに(図14(B))、各変形速度ベクトル分布を空間微分して得られたひずみ速度テンソル分布を順次記憶する(図14(C))。なお、ひずみ速度テンソル分布は、図示のように、所定範囲ごとに識別可能に断層画像の態様で分布表示することができる。

【0096】
図15は、ひずみ速度可視化画面を示す図である。(A)は本実施形態によるひずみ速度可視化画面を示す。(B)は比較例として、FEMにより本実施形態と同じ条件にてひずみ速度解析をした結果を示す。各図の横軸はx方向の位置を、縦軸はz方向の位置をそれぞれ示している。

【0097】
図15(A)および(B)から、ひずみ速度分布の演算および可視化処理に関し、本実施形態によるOCTによる演算結果が、FEM解析の結果と精度良く対応していることが確認できる。応力分布について高精度な結果が得られることは第1実施形態で述べたとおりであるから、これら応力分布とひずみ速度分布とに基づいて、力学物性値の断層分布を精度良く得ることが可能となる。

【0098】
次に、制御演算部4が実行する具体的処理の流れについて説明する。
図16は、制御演算部4により実行される力学物性値可視化処理の流れを示すフローチャートである。制御演算部4は、OCTによる測定対象Wの断層画像を連続的に取得する(S210)。そして、その断層画像を用いて応力の断層分布を演算し、画面に表示する(S212)。また、その断層画像を用いひずみ速度の断層分布を演算し、画面に表示する(S214)。そして、それら応力断層分布およびひずみ速度断層分布の演算結果から力学物性値の断層分布を演算し、画面に表示する(S216)。ここで、S212の応力断層分布の演算処理等については、第1実施形態の処理を適用することができる。このため、その説明については省略する。

【0099】
図17は、図16におけるS214のひずみ速度分布の演算可視化処理を詳細に示すフローチャートである。制御演算部4は、連続撮影された時刻の異なる前後2枚のOCT断層画像I(x,z,t)とI(x,z,t+Δt)を取得し、再帰的相互相関法による処理を実行する(S220)。各OCT断層画像は、図9のS10で取得された水平偏光成分I(x,z)および垂直偏光成分I(x,z)の二乗和として得られる。すなわち、上記式(4)のストークスパラメータSが、そのOCT断層画像に対応する。

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

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

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

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

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

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

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

【0107】
上記実施形態では、図9のS10において、2軸の断層画像I(x,z),I(x,z)を取得して処理を進め、S34にて応力分布をx-z平面に二次元的に可視化する例を示した。変形例においては、特定のx座標についての1軸断層画像I(z)を取得し、z方向の応力分布として可視化してもよい。

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

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

【0110】
なお、本発明は上記実施形態や変形例に限定されるものではなく、要旨を逸脱しない範囲で構成要素を変形して具体化することができる。上記実施形態や変形例に開示されている複数の構成要素を適宜組み合わせることにより種々の発明を形成してもよい。また、上記実施形態や変形例に示される全構成要素からいくつかの構成要素を削除してもよい。
図面
【図1】
0
【図2】
1
【図3】
2
【図4】
3
【図5】
4
【図6】
5
【図7】
6
【図8】
7
【図9】
8
【図10】
9
【図11】
10
【図12】
11
【図13】
12
【図14】
13
【図15】
14
【図16】
15
【図17】
16