TOP > 国内特許検索 > マイクロ断層可視化方法およびシステム > 明細書

明細書 :マイクロ断層可視化方法およびシステム

発行国 日本国特許庁(JP)
公報種別 特許公報(B2)
特許番号 特許第6707249号 (P6707249)
公開番号 特開2016-170080 (P2016-170080A)
登録日 令和2年5月22日(2020.5.22)
発行日 令和2年6月10日(2020.6.10)
公開日 平成28年9月23日(2016.9.23)
発明の名称または考案の名称 マイクロ断層可視化方法およびシステム
国際特許分類 G01N  21/17        (2006.01)
G01N  23/046       (2018.01)
FI G01N 21/17 630
G01N 23/046
請求項の数または発明の数 9
全頁数 29
出願番号 特願2015-050553 (P2015-050553)
出願日 平成27年3月13日(2015.3.13)
審査請求日 平成30年2月13日(2018.2.13)
特許権者または実用新案権者 【識別番号】599002043
【氏名又は名称】学校法人 名城大学
発明者または考案者 【氏名】佐伯 壮一
個別代理人の代理人 【識別番号】110002273、【氏名又は名称】特許業務法人インターブレイン
審査官 【審査官】越柴 洋哉
参考文献・文献 特開2008-136880(JP,A)
特開2007-244533(JP,A)
特開2004-057652(JP,A)
特開2004-057653(JP,A)
特開2011-177535(JP,A)
国際公開第2015/033058(WO,A1)
調査した分野 G01N 21/00-21/61
G01N 23/00-23/2276
A61B 1/00- 1/32
JSTPlus/JMEDPlus/JST7580(JDreamIII)
特許請求の範囲 【請求項1】
構造を対象とする順解析の過程で変位関連ベクトルが演算される所定の数値解析において解析条件となり得る物理量に関し、前記構造の内部における前記物理量の分布を断層可視化するマイクロ断層可視化方法であって、
前記構造の変形前後の断層画像を取得する工程と、
前記断層画像に基づき変位関連ベクトルの断層分布を演算する工程と、
前記断層画像に対して関心領域を設定する工程と、
前記関心領域を要素分割する工程と、
演算された変位関連ベクトルを、分割された要素の節点に補間するとともに境界条件として与える工程と、
補間後の変位関連ベクトルの断層分布を用いて前記数値解析の逆解析を境界条件の設定を別途行わずに実行することにより、各要素における前記物理量の分布を前記物理量の断層分布として演算する工程と、
演算された前記物理量の分布を表示装置に断層可視化する工程と、
を備えることを特徴とするマイクロ断層可視化方法。
【請求項2】
前記関心領域が、前記断層画像の内部に対して設定されることを特徴とする請求項1に記載のマイクロ断層可視化方法。
【請求項3】
前記数値解析が、状態変化又は境界条件変化による要因に伴う前記構造の変形を考慮した応力解析であることを特徴とする請求項1または2に記載のマイクロ断層可視化方法。
【請求項4】
前記物理量が温度であり、
前記数値解析が有限要素法による熱変形を考慮した熱応力解析であることを特徴とする請求項3に記載のマイクロ断層可視化方法。
【請求項5】
前記物理量が前記構造の物理定数であり、
前記数値解析が有限要素法による応力解析であることを特徴とする請求項3に記載のマイクロ断層可視化方法。
【請求項6】
光コヒーレンストモグラフィーを用いて前記構造の断層画像を取得することを特徴とする請求項1~5のいずれかに記載のマイクロ断層可視化方法。
【請求項7】
前記構造として、少なくとも高分子基材料又は生体組織を対象として含むことを特徴とする請求項1~6のいずれかに記載のマイクロ断層可視化方法。
【請求項8】
前記構造における所定の領域において前記物理量の真値を与える工程と、
前記逆解析により得られた物理量のうち、前記領域に対応する物理量である推定値と、前記領域において与えられた真値との関係に基づいて、前記逆解析により得られた物理量の値を補正する工程と、
をさらに備えることを特徴とする請求項1~7のいずれかに記載のマイクロ断層可視化方法。
【請求項9】
測定対象となる構造の内部における所定の物理量の分布を断層可視化するマイクロ断層可視化システムであって、
前記構造の変形前後の断層画像を取得するための断層画像取得装置と、
前記断層画像取得装置から前記構造の変形前後の断層画像データを取得し、その断層画像データに基づいて前記物理量の断層分布を演算する制御演算部と、
前記制御演算部の演算結果に基づいて、前記構造における前記物理量の分布を断層可視化する態様で表示する表示装置と、
を備え、
前記物理量は、前記構造を対象とする順解析の過程で変位関連ベクトルが演算される所定の数値解析において解析条件となり得るものであり、
前記制御演算部は、
取得した断層画像データに基づき変位関連ベクトルの断層分布を演算し、
前記断層画像に対して関心領域を設定し、
前記関心領域を要素分割し、
演算された変位関連ベクトルを、分割された要素の節点に補間するとともに境界条件として与え、
補間後の変位関連ベクトルの断層分布を用いて前記数値解析の逆解析を境界条件の設定を別途行わずに実行することにより、各要素における前記物理量の分布を前記物理量の断層分布として演算し、
演算された前記物理量の分布を前記表示装置に断層可視化させることを特徴とするマイクロ断層可視化システム。
発明の詳細な説明 【技術分野】
【0001】
本発明は、構造内部における物理量のマイクロ断層可視化を可能にする方法およびシステムに関する。
【背景技術】
【0002】
パーソナルコンピュータや携帯端末などの電子機器には、いわゆるMEMS(Micro Electro Mechanical Systems)と呼ばれる電子実装部品が組み込まれている。このような電子実装部品の小型化,高密度化に伴うジュール発熱密度の増大により、熱ひずみによる材料内部の微視的破壊が問題となっている。特に、導電性フィラーを含有する繊維強化プラスチックなど(CFRP,GFRP等)、温度依存導電性を有する高分子基複合材料からなる部品について、その問題が顕著となっている。また、例えばリチウムイオン二次電池用キャパシタは、薄板電極間における温度上昇がドライアップ故障を引き起こすといった問題を有する。このため、これらの工業製品にはマイクロスケールでの最適熱設計が必要とされる。
【0003】
このような工業製品の熱設計においては、構造内部の温度評価が重要となるが、内部の温度情報が所望の位置でマイクロスケールにて得られる実計測システムは未だ存在せず、有限要素法等の数値解析に頼るところが大きい(例えば特許文献1参照)。また、設計という観点からは、このような熱設計だけでなく、強度(剛性)など力学特性を考慮した設計も重要となるが、これも同様の数値解析に頼るところが大きい。さらに、工業製品を離れると、例えば生体組織ではその力学物理量を測定することが困難であるため、これも同様の数値解析に依存する。
【先行技術文献】
【0004】

【特許文献1】特開2006-284249号公報
【発明の概要】
【発明が解決しようとする課題】
【0005】
しかしながら、このような熱解析や力学解析などの数値解析をマイクロスケールにて行う場合、それらにおける境界条件の与え方が難しい。また、局所的なデータの取得を目的とする場合にも部品全体を解析対象とする必要があり、演算効率上の問題も生じる。
【0006】
本発明はこのような課題に鑑みてなされたものであり、その目的の一つは、構造内部における物理量を簡易な手法にてマイクロ断層可視化可能な技術を提供することにある。
【課題を解決するための手段】
【0007】
本発明のある態様は、構造を対象とする順解析の過程で変位関連ベクトルが演算される所定の数値解析において解析条件又は解析結果となり得る物理量に関し、構造の内部における物理量の分布を断層可視化するマイクロ断層可視化方法である。この方法は、構造の変形前後の断層画像を取得する工程と、断層画像に基づき変位関連ベクトルの断層分布を演算する工程と、演算された変位関連ベクトルの断層分布を用いて数値解析の逆解析を実行することにより、物理量の断層分布を演算する工程と、演算された物理量の分布を表示装置に断層可視化する工程と、を備える。なお、ここでいう「変位関連ベクトル」は、変形ベクトルであってもよいし、変形速度ベクトルであってもよい。
【0008】
この態様によると、実測により得られる断層画像から変位関連ベクトルの分布が演算され、その変位関連ベクトル分布に基づき所定の数値解析の逆解析を行うことで物理量の分布が得られる。すなわち、実測と数値解析の逆解析とのハイブリッドにより物理量を断層可視化することが可能となる。実測により内部の変位関連ベクトルが得られるため、数値解析の順解析のように境界部に境界条件を設定する必要もない。このため、所望の断層位置について物理量を簡易に算出することができる。すなわち、この態様によれば、構造内部における物理量を簡易な手法にてマイクロ断層可視化することができる。
【0009】
本発明の別の態様は、測定対象となる構造の内部における所定の物理量の分布を断層可視化するマイクロ断層可視化システムである。このシステムは、構造の変形前後の断層画像を取得するための断層画像取得装置と、断層画像取得装置から構造の変形前後の断層画像データを取得し、その断層画像データに基づいて物理量の断層分布を演算する制御演算部と、制御演算部の演算結果に基づいて、構造における物理量の分布を断層可視化する態様で表示する表示装置と、を備える。
【0010】
物理量は、構造を対象とする順解析の過程で変位関連ベクトルが演算される所定の数値解析において解析条件又は解析結果となり得るものである。なお、ここでいう「変位関連ベクトル」は、変形ベクトルであってもよいし、変形速度ベクトルであってもよい。制御演算部は、取得した断層画像データに基づき変位関連ベクトルの断層分布を演算し、演算された変位関連ベクトルの断層分布を用いて数値解析の逆解析を実行することにより、物理量の断層分布を演算し、演算された物理量の分布を表示装置に断層可視化させる。
【0011】
この態様によると、断層画像取得装置にて取得される断層画像から変位関連ベクトルの分布が演算され、その変位関連ベクトル分布に基づき所定の数値解析の逆解析を行うことで物理量の分布が得られる。すなわち、断層画像取得装置による実測と、制御演算部による数値解析の逆解析とのハイブリッドにより、物理量を断層可視化することが可能となる。実測により内部の変位関連ベクトルが得られるため、数値解析の順解析のように境界部に境界条件を設定する必要もない。このため、所望の断層位置について物理量を簡易に算出できる。すなわち、この態様によれば、構造内部における物理量を簡易な手法にてマイクロ断層可視化することができる。
【発明の効果】
【0012】
本発明によれば、構造内部における物理量を簡易な手法にてマイクロ断層可視化可能な技術を提供することができる。
【図面の簡単な説明】
【0013】
【図1】実施例に係る断層可視化システムの構成を概略的に表す図である。
【図2】実施例に係るマイクロ断層可視化方法の概念を表す図である。
【図3】再帰的相互相関法による処理手順を概略的に示す図である。
【図4】サブピクセル解析による処理手順を概略的に示す図である。
【図5】有限要素モデルを概略的に示す図である。
【図6】逆解析の処理手順を概略的に示す図である。
【図7】逆解析の処理手順を概略的に示す図である。
【図8】逆解析の処理手順を概略的に示す図である。
【図9】逆解析の処理手順を概略的に示す図である。
【図10】制御演算部により実行される物理量可視化処理の流れを示すフローチャートである。
【図11】図10におけるS12の変形ベクトル演算処理を詳細に示すフローチャートである。
【図12】図10におけるS14の逆解析処理を詳細に示すフローチャートである。
【図13】数値実験による検証方法を示す図である。
【図14】数値実験による検証方法を示す図である。
【図15】数値実験による検証方法を示す図である。
【図16】検証結果を示す図である。
【図17】本実施例の解析による正解値との誤差を示す図である。
【発明を実施するための形態】
【0014】
本発明の一実施形態は、構造内部の物理量の分布をマイクロスケールにて断層可視化する方法である。なお、ここでいう「構造」は、材料構造や生体構造(生体組織)を含み得る。「断層可視化」は、少なくともマイクロスケールにて可能であるが、それ以外のスケールを排除するものではない。断層画像取得装置の解像度に応じてナノスケールにて可能なものも含み得る。「物理量」は、所定の数値解析において解析条件又は解析結果の一つとなり得るものであればよく、後述する温度分布など構造の状態変化に係るものでもよいし、物理定数など構造の物性値に係るものでもよい。その場合、「物理定数」は、「弾性係数」,「粘性係数」,「透水係数」等を含み得る。「断層画像」は、光コヒーレンストモグラフィー(Optical Coherence Tomography:以下「OCT」という)によるものでもよい。あるいは、X線CTや超音波断層像によるものでもよい。

【0015】
この方法は、数値解析の順解析の過程で変位関連ベクトルの分布が演算される点と、断層画像取得装置により実測される断層画像から変位関連ベクトル分布を演算可能である点に着目し、両者の変位関連ベクトルを結び付ける。ここで、「変位関連ベクトル」は、構造内部の変位や変形に関連するベクトルであり、変形ベクトルであってもよいし、変形速度ベクトルであってもよい。「変形ベクトル」は「変位ベクトル」の概念を含む。すなわち、実測と数値解析とを後者の「逆解析」という形で結び付け、両者のハイブリッドによる演算処理を実現している。ここでいう「逆解析」は、順解析の一部を逆方向にたどることにより実現されてよい。

【0016】
すなわち、順解析においては、例えば上記物理量が所定の数理モデルの解析条件として入力され、その数理モデルに沿った演算処理を経ることで解が算出される。その演算過程で変位関連ベクトルが演算される。「逆解析」は、逆にその変位関連ベクトルの分布を既知として順解析を逆方向にたどり、解析条件の一つを構成する「物理量」の分布を算出する。その変位関連ベクトル分布として、実測された断層画像に基づいて算出されたものが用いられる。なお、「変形ベクトル」の時間微分が「変形速度ベクトル」であるが、「変形速度ベクトル」を変形ベクトルの時間変化とみて「変形ベクトル」の概念で捉えてもよい。いずれにしても、演算過程において変形ベクトルが演算される。本実施形態によれば、OCTやCT等の断層可視化手法から逆解析に結び付ける斬新な手法が提供される。

【0017】
この方法では、断層画像取得装置を介して構造の変形前後の断層画像が取得される。なお、「変形前後の断層画像」については、構造の変形前と変形後のそれぞれについて変化が停止している状態で取得してもよい。あるいは、変形中の前後の時点で取得してもよい。すなわち、本手法は、後者のように時間的に変化する現象にも適用可能である。具体的には、OCT断層画像の動画中の2枚の断層画像に対して本手法を適用することにより、それらの断層画像に対応する温度分布を断層可視化することが可能である。前後の断層画像を連続的に取得すれば、温度分布をリアルタイムに断層可視化することもできる。

【0018】
そして、その断層画像に基づいて断層位置に対応した変位関連ベクトルの分布が演算される一方、断層画像の内部に対して関心領域(Region of Interest:以下「ROI」と表記する)が設定される。また、上記数値解析の逆解析を利用するためにROIに対して要素分割がなされる。これらROIおよび要素分割の設定は、入力装置を介したユーザの指示入力によりなされる。

【0019】
断層画像に基づいて演算された変位関連ベクトルは、分割された要素の節点に補間される。この補間後の変形ベクトルの分布に対して上記逆解析を施すことにより、各要素における物理量の分布が演算される。そして、このとき演算された物理量の分布が、表示装置に断層可視化される。

【0020】
上記数値解析は、状態変化又は境界条件変化による要因に伴う構造の変形を考慮した応力解析であってもよい。例えば、有限要素法(Finite Element Method:以下「FEM」と表記する)、有限差分法(Finite Difference Method:以下「FDM」と表記する)、あるいは境界要素法(Boundary Element Method:以下「BEM」と表記する)によるものでもよい。

【0021】
上記手法は、構造における所定の領域において物理量の真値を与える工程と、逆解析により得られた物理量のうち上記領域に対応する物理量である推定値と、上記領域において与えられた真値との関係に基づき、逆解析により得られた物理量の値を補正する工程と、をさらに備えるものでもよい。

【0022】
「所定の領域」は、複数設定されてもよく、構造の表面又はその近傍に設定されてもよいし、構造内部の特定の部材の周囲に沿って設定されてもよい。「真値」は、実測値であってもよいし、解析結果(解析値)であってもよく、正しいと推定される値であればよい。例えば、「物理量」が温度であり、赤外線サーモグラフィーや熱電対によりROI近傍の温度値を計測できる場合、その実測値である温度値を真値としてもよい。あるいは、銅配線の電流と電圧によって算出されるジュール熱などに基づいてROI近傍の計測データを算出できる場合、その演算結果を真値としてもよい。この「真値」は、ROI内部の値であることが望ましいが、ROI近傍の温度値(実験計測値)であってもよい。また、熱伝導率などを仮定してモデリングし、ROI内部の任意領域の温度値を推定し、その推定値を「真値」としてもよい。

【0023】
上記数値解析は、状態変化又は境界条件変化による要因に伴う構造の変形を考慮した応力解析であってもよい。ここでいう「要因」は、構造への通電であってもよいし、構造の周囲温度の変化であってもよい。あるいは、構造に外力が付与されることでもよい。

【0024】
例えば、力学的物理量のマイクロ断層可視化に際して下記式(a)を利用する。
[K]{U}={F} ・・・(a)
ここで、[K]は全体剛性マトリクスであり、{U}は要素分割されたROIの各節点における変形ベクトルであり、{F}は外力ベクトルである。上記式(a)において、{F}に関し、ROI内部の節点の荷重はゼロ(合力ゼロ)であり、ROI境界の節点の荷重に対しては変位境界条件(変位拘束条件)を与える。{U}については、変形前後の断層画像に基づいて演算される変位関連ベクトル分布(変形ベクトルの断層分布等)から得られるため、それを代入することができる。上記式(a)を逆解析することにより、物理量としての[K]を推定することができる。

【0025】
以下、図面を参照しつつ本実施形態を具体化した実施例について詳細に説明する。
[実施例]
図1は、実施例に係る断層可視化システムの構成を概略的に表す図である。本実施例の断層可視化システムは、電子実装部品の内部の特定の物理量をマイクロスケールにて断層可視化するものである。この断層可視化にOCTを利用する。

【0026】
図1に示すように、断層可視化システム1は、OCTを用いる光学系を含む光学ユニット2、光学ユニット2に接続される光学機構4、OCTにより得られた光干渉データに基づいて物理量を演算する制御演算部6、ユーザの指示入力を受け付ける入力装置7、および物理量を断層可視化する態様で表示する表示装置8を備える。光学ユニット2および光学機構4は、「断層画像取得装置」を構成する。図示の例では、光学ユニット2としてマッハツェンダー干渉計をベースとした光学系が示されているが、マイケルソン干渉計その他の光学系を採用することもできる。

【0027】
本実施例では、OCTとして波長走査型レーザを光源とするSS-OCT(Swept Source OCT)を用いるが、TD-OCT(Time Domain OCT)、SD-OCT(Spectral Domain OCT)その他のOCTを用いることもできる。SS-OCTは、TD-OCTとは異なり、参照鏡走査などの機械的な光遅延走査を必要としないため、高い時間分解能と高い位置検出精度が得られる点で好ましい。

【0028】
光学ユニット2は、光源10、オブジェクトアーム12、リファレンスアーム14、および光検出器16を備える。各光学要素は、光ファイバーにて互いに接続されている。光源10から出射された光は、カプラ18(ビームスプリッタ)にて分けられ、その一方がオブジェクトアーム12に導かれる物体光となり、他方がリファレンスアーム14に導かれる参照光となる。オブジェクトアーム12に導かれた物体光は、サーキュレータ20を介して光学機構4に導かれ、測定対象Wに照射される。この物体光は、測定対象Wの表面および断面にて後方散乱光として反射されてサーキュレータ20に戻り、カプラ22に導かれる。

【0029】
一方、リファレンスアーム14に導かれた参照光は、サーキュレータ24を介して反射鏡26に導かれる。このとき、参照光は、コリメートレンズ28を経て集光レンズ30にて反射鏡26に集光される。この参照光は、反射鏡26にて反射されてサーキュレータ24に戻り、カプラ22に導かれる。すなわち、物体光と参照光とがカプラ22にて合波(重畳)され、その干渉光が光検出器16により検出される。

【0030】
光学機構4は、光学ユニット2のオブジェクトアーム12を構成し、サーキュレータ20から延びる光ファイバー34が接続されている。光学機構4は、光学ユニット2からの光を測定対象Wに導いて走査させる機構と、その機構を駆動するための駆動部(アクチュエータ)を備える。光学機構4は、例えばガルバノ鏡を含むガルバノ装置であってもよい。光学機構4は、測定対象Wに向けて物体光を照射する。

【0031】
測定対象Wは、MEMSに対応した電子実装部品であり、ICチップを基板に埋め込んで構成される。その基板は、高分子樹脂材料等の高分子基材料やプラスチック等の結晶性材料からなるものでもよい。その基板は、内部に反射する構造(界面など)がある材料からなるものであればよく、ナノ粒子などの散乱体が樹脂に内在する複合材料など、透明でない材料からなるものでもよい。

【0032】
カプラ22にて合波された干渉光は、光検出器16に入力される。光検出器16は、これを光干渉信号(干渉光強度を示す信号)として検出する。この光干渉信号は、A/D変換器40を介して制御演算部6に入力される。A/D変換器40は、光検出器16から出力されたアナログ信号をデジタル信号に変換して制御演算部6へ出力する。

【0033】
制御演算部6は、CPU、ROM、RAM、ハードディスクなどを有し、これらのハードウェア、およびソフトウェアによって、光学ユニット2の光学系全体の制御と、光学機構4の駆動制御と、OCTによる画像出力のための演算処理を行う。制御演算部6の指令信号は、D/A変換器42を介して光学機構4等に入力される。D/A変換器42は、制御演算部6から出力されたデジタル信号をアナログ信号に変換して光学機構4等へ出力する。制御演算部6は、光学機構4の駆動に基づいて光学ユニット2から出力された光干渉信号を処理し、OCTによる測定対象Wの断層画像を取得する。そして、その断層画像データに基づき、後述の手法により測定対象Wの内部における特定の物理量の断層分布を演算する。

【0034】
入力装置7は、キーボードやタッチパネル等からなり、ユーザの指示入力を受け付け、制御演算部6に入力する。表示装置8は、例えば液晶ディスプレイからなり、制御演算部6にて演算された測定対象Wの内部の物理量を断層可視化する態様で画面に表示させる。

【0035】
以下、上記物理量の断層分布を可視表示するまでの演算処理方法について説明する。
図2は、実施例に係るマイクロ断層可視化方法の概念を表す図である。本実施例では、測定対象Wの内部の温度を上記物理量としてその分布を断層可視化する。この手法は、測定対象Wの熱変形に基づいて内部の温度分布を演算するものであり、OCTにて取得した断層画像から変形ベクトル分布を演算する工程(変形ベクトル演算工程)と、その変形ベクトル分布に対して有限要素法(以下「FEM」と表記する)の逆解析を施すことにより、測定対象Wの内部の温度分布を演算する工程(逆解析工程)とを有する。

【0036】
まず、変形ベクトル演算工程について説明する。上述のように、OCTにおいて、オブジェクトアーム12を経た物体光(測定対象Wからの反射光)と、リファレンスアーム14を経た参照光とが合波され、光検出器16により光干渉信号として検出される。制御演算部6は、この光干渉信号を干渉光強度に基づく測定対象Wの断層画像として取得することができる。

【0037】
ところで、OCTの光軸方向(奥行き方向)の分解能であるコヒーレンス長lは、光源の自己相関関数によって決定される。ここでは、コヒーレンス長lを自己相関関数の包括線の半値半幅とし、下記式(1)にて表すことができる。
【数1】
JP0006707249B2_000002t.gif
ここで、λはレーザビームの中心波長であり、Δλはレーザビームの半値全幅である。

【0038】
一方、光軸垂直方向(ビーム走査方向)の分解能は、集光レンズによる集光性能に基づき、ビームスポット径Dの1/2とされる。そのビームスポット径Dは、下記式(2)にて表すことができる。
【数2】
JP0006707249B2_000003t.gif
ここで、dは集光レンズに入射するビーム径、fは集光レンズの焦点である。このようにOCTによる分解能には限界があるところ、本実施例では後述するサブピクセル解析の導入などにより、温度分布の断層表示をマイクロスケールにて行うことを可能にしている。以下、その詳細について説明する。

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

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

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

【0042】
このようにして得られた変形ベクトル分布を空間微分することにより、ひずみ断層分布を演算することができる。なお、本実施例では後述のように、変形ベクトル分布に対してFEMの逆解析を施すことにより、測定対象Wの内部の温度分布を演算する。このため、ひずみ断層分布まで演算する必要はなく、その演算過程である変形ベクトル分布の演算までが利用されることとなる。

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

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

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

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

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

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

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

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

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

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

【0053】
サブピクセル解析では、注目点における変形前後の輝度差が各成分の輝度勾配と移動量によって表される。このため、検査領域S1内の輝度勾配データより最小二乗法を用いてサブピクセル移動量を決定することができる。本実施例では、輝度勾配を求める際に、サブピクセル変形前の風上側の輝度勾配を与える風上差分法を採用している。すなわち、サブピクセル解析は様々な手法が存在するが、本実施例では検査領域サイズが小さく高空間解像度の条件においても、サブピクセル移動量を高精度検出する勾配法を採用している。

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

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

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

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

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

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

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

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

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

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

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

【0065】
ここで、Δ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は下記式(11)にて表される。
【数11】
JP0006707249B2_000012t.gif

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

【0067】
次に、逆解析工程について説明する。
図5は、有限要素モデルを概略的に示す図である。(A)は三角形要素を示し、(B)は簡易モデルを示す。図6~図9は、逆解析の処理手順を概略的に示す図である。各図の(A)~(C)は、その処理過程を示す。図8は、図7(C)の部分拡大図に対応する。

【0068】
本実施例では、OCT断層像から求められる変形ベクトル断層分布とFEMとを組み合わせることで温度の断層分布を逆解析推定する。図6(A)には、OCTにより得られる電子実装部品(測定対象W)の断層画像が示されている。横軸がx座標を示し、縦軸がz座標を示す。同図には、高分子基材料からなる基板にICチップが埋設された電子実装部品の一部が示されている。

【0069】
図6(B)には、OCTにより得られた前後の断層画像を用いて上記演算処理を行うことにより得られる変形ベクトル分布が矢印にて示されている。このような変形ベクトル分布に対してFEMの逆解析を施すことにより、測定対象内部の温度分布を算出し、断層表示することができる。

【0070】
ここで、「逆解析」の意義を明確にするために、FEMの順解析について理解しておく必要がある。この順解析では、解析対象となる領域(解析領域)を要素分割した後、材料パラメータを設定するとともに荷重等の境界条件を設定し、各要素における温度を与えた後、多元連立方程式となる全体剛性方程式を解く。それにより、各節点における変位(未知項)を算出することができる。なお、本実施例では温度分布と変形分布とが連成するため、順解析は次のようになる。すなわち、解析領域を要素分割した後、熱伝導係数分布を設定するとともに熱(温度)の境界条件を設定し、方程式を解いて温度分布を算出する過程を経る。また、解析領域を要素分割した後、材料パラメータ(剛性、線膨張係数など),力学的境界条件,温度分布を設定し、方程式を解いて変形分布を算出する過程を経る。

【0071】
これに関し、本実施例では上述のように、OCTの断層画像に基づいて変形ベクトル分布が算出されるため、解析領域に設定された節点に対してそれらの変形ベクトルを補間することにより、境界条件を設定しなくとも節点上の変位を算出することができる。すなわち、全体剛性方程式における変位ベクトル(変形ベクトル)の項が既知となる。この点に着目し、この既知となった変位ベクトルを用いてFEMの逆解析を実行することにより温度項を算出する。このとき、全節点における変形ベクトルが既知であるため、全要素について温度項を算出することができる。つまり、全要素に対応した温度分布を算出することができる。以下、具体的に説明する。

【0072】
(FEMによる温度項を含む全体剛性方程式)
二次元平面応力問題に関し、熱膨張を考慮した応力-ひずみ関係式として下記式(12)が得られる。なお、ここではOCTの断層画像に合わせてx-z平面を考える。
【数12】
JP0006707249B2_000013t.gif
ここで、σは応力のx成分、σは応力のz成分、τxzはせん断応力である。εはx方向のひずみ、εはz方向のひずみ、γxzはせん断ひずみである。Eはヤング率、νはポアソン比、αは線膨張係数、ΔTは温度変化である。

【0073】
そして、上記式(12)に仮想仕事の原理を適用して変形し、要素剛性マトリクス[k]6×6を用いて表すと、下記式(13)が得られる。
【数13】
JP0006707249B2_000014t.gif

【0074】
本実施例ではFEMを利用するに際し、要素分割される各要素を図5(A)に示す三角形要素にて表し、また、隣接する要素間の関係を図5(B)に示す簡易モデルにて表す。上記式(13)において、P,Qは等価節点力であり、u,vはそれぞれ節点位置のx方向,z方向の変形ベクトル成分である。なお、P,Q,u,vの添え字1~3は、図5(A)に示す要素分割に係る三角形要素の局所節点番号である。a~bは三角形要素において物理量を近似する際の定数であり、それぞれ下記式(14)により表される。
【数14】
JP0006707249B2_000015t.gif
ここで、x,zの添え字1~3は要素分割に係る三角形要素の局所節点番号である。

【0075】
本実施例では、図7(A)に示すように、測定対象の内部にROIを設定し、図7(B)に示すように、そのROIに対してFEMに基づく要素分割を実行する。OCTにより得られた変形ベクトルは、図7(C)に示すようにROIに分布する。

【0076】
上記式(13)を温度の項について展開し、図5(B)に示す簡易モデルを用いて全体方程式を作成すると、下記式(15)が成立する。
【数15】
JP0006707249B2_000016t.gif
ここで、上付き文字,下付文字はそれぞれ要素番号と節点番号を表している。

【0077】
上記式(15)を簡略化および一般化すると下記式(16)となり、線形化された観測方程式となる。
【数16】
JP0006707249B2_000017t.gif
ここで、[A](2min+2mb)×nは、温度変化項の観測マトリクスを表し、FEMの形状関数によって決定される[B]マトリクスに依存し、それらの成分は節点座標により表される。{ΔT}n×1は各要素における温度変化、[K](2min+2mb)×(2min+2mb)は全体剛性マトリクス、{U}(2min+2mb)×1は各節点における変形ベクトル成分となる。この変形ベクトル成分はOCT断層像に基づいて得られるため、右辺は観測ベクトルとなる。

【0078】
ただし、minはROI境界上を除くROI内の節点数であり、2minはそれらの節点上に求まる変形ベクトルの成分数である。mbはROI境界上の節点数であり、2mbはそれらの節点上に求まる変形ベクトルの成分数である。したがって、2min+2mbはROI以内の全節点に求まる変形ベクトルの成分数となる。nはROI内の要素数である。上記式(16)の連立方程式に境界変位条件を導入することで修正した全体剛性方程式に対し、種々の逆解析手法を用いることによって、各要素における温度変化ΔTを推定することができる。

【0079】
(境界変位条件の導入)
ROI内の有限要素モデルに上記式(16)を適用した後、各節点における変形ベクトル成分{U}(2min+2mb)×1を節点補間により与える。このとき、ROI境界上の節点に得られる変形ベクトル成分u,vを境界強制変位(拘束節点)とする必要がある。そのため、上記式(16)における全体剛性マトリクス[K](2min+2mb)×(2min+2mb)、および観測マトリクス[A](2min+2mb)×nを修正してマトリクスの縮小を行う。

【0080】
すなわち、マトリクス成分の添え字の意味を改変した上で上記式(16)を改めて書き下すと下記式(17)となる。
【数17】
JP0006707249B2_000018t.gif

【0081】
この全体剛性方程式を解く過程で境界変位条件を満足するように修正する。ここで、ROI境界上の節点に得られた変形ベクトル成分をu,vとする。このとき、{U}(2min+2mb)×1において境界上の節点に求まる変形ベクトル成分を表す行はuが2l-1行目、vが2l行目となる。ROI境界の節点では境界変位条件(変位拘束境界)として、ROI境界上節点の既知の変位量u,vについて左右両辺が満たされるように、全体剛性マトリクス[K](2min+2mb)×(2min+2mb)、および観測マトリクス[A](2min+2mb)×nを修正する。具体的には、ROI境界上節点における変位量u,vが,左辺と右辺においてu=u,v=vとなるように、[K]の成分と、[A]における同じ行の成分を下記式(18)に従って修正する。
【数18】
JP0006707249B2_000019t.gif

【0082】
また、ベクトル{F}の境界上の節点に対応する2l-1,2l行目も、それぞれu,vと修正して左右両辺が満たされるようにする。さらに、ROIが物体内部に設定され、かつ静的変形により体積力が作用しないのであれば、{F}に関してROI境界上を除くROI内の等価節点力はゼロとなる。また、ROI境界を自由表面に設定した場合は自由変形となるため、表面力は{F}=0となる。このため、上記式(17)は下記式(19)のようになる。
【数19】
JP0006707249B2_000020t.gif

【0083】
以上のようにして、全体剛性方程式に境界変位条件を満足させる。なお、上記式(19)において、両辺の2l-1,2l行目は拘束節点に帰属し、境界変位条件を満たすため、それらを省略することで成分を並び替え、マトリクスの縮小を行う。このようにして得られる観測マトリクス[A]を[A']と表し、全体剛性マトリクス[K]と変形ベクトル{U}の積を{K'U'}と表すと、下記式(20)が得られる。
【数20】
JP0006707249B2_000021t.gif

【0084】
(節点補間)
この逆解析のためには、上記式(16)における測定項{U}を与えなければならない。一方、図8(A)に示すように、OCTにより得られた変形ベクトル分布は、ROI内に分割された要素の節点に対応していないため、これをそのまま測定項{U}とすることはできない。そこで、図8(B)に示すようにROIの節点上に変形ベクトルを設定するための補間処理を実行する。

【0085】
すなわち、OCTにより取得した変形ベクトルの断層分布を用いて各節点上に変形ベクトルを補間する。ここでは、各節点を中心とする所定半径(本実施例では300μm)を「影響半径」とし、その影響半径に収まる変形ベクトルをその節点に補間した。この補間に際しては、最小二乗法により下記式(21)に示す近似関数f(x,z)を求め、その近似関数f(x,z)に則って各節点位置の変形ベクトルを算出した。
【数21】
JP0006707249B2_000022t.gif
ここで、f(x,z)は二次元二次多項式とした近似関数である。a~gはその未知係数であり、x,zはOCTにより得られる変形ベクトルの始点座標値である。

【0086】
この近似関数f(x,z)について下記式(22)を適用する。
【数22】
JP0006707249B2_000023t.gif
ここで、Qは、近似関数f(x,z)と測定値yとの残差の二乗和である。測定値yはOCTにより得られる変形量成分であり、Nは影響半径内にある変形ベクトルのデータ数である。上記式(22)においてQが最小となるときの近似関数f(x,z)を求め、そのf(x,z)に則ってROIにある節点上の変形量成分を算出する。

【0087】
具体的には、Qを算出した後、Qを近似関数の各未知係数について偏微分することにより、下記式(23)に示す正規式を得る。
【数23】
JP0006707249B2_000024t.gif

【0088】
この正規式を解き、未知係数を算出することで近似関数を求める。求まった近似関数に節点座標を参照させ、節点位置での変形量成分を求めることで節点補間がなされる。ただし、表面近傍のROI境界や深部のROI境界では輝度情報にノイズが多く含まれるため、変形ベクトルの算出精度が低い。さらに、ROI境界近傍、特にROIの頂点部分(角部分)では使用できる変形ベクトルの数が少ないため、ROI境界外に変形ベクトルを外挿補間している。このような手法により、ROI内の節点上の各変形量成分を求めて節点補間をする。なお、上記式(16)の{U}は、OCTによる測定データとして得られる。

【0089】
(逆解析)
上記式(20)における観測マトリクス[A']は一般に2min×nの行列である。minはROI境界を除いたROI内部の節点x,z座標の数であり、nはROI内の要素数である。min<nとなるため、[A']は正則ではない。そのため、本実施例では、逆解析手法として擬似逆行列であるムーアペンローズ一般逆行列を求め、逆解析を行った。ここで、観測マトリクス[A']がフルランクでランクがnの場合、[A']の一般逆行列は下記式(24)にて表される。
【数24】
JP0006707249B2_000025t.gif

【0090】
上記式(24)に基づき、温度の逆解析を下記式(25)を用いて実行する。
【数25】
JP0006707249B2_000026t.gif
本実施例では、この[A']のムーアペンローズ一般逆行列を上記式(20)の両辺に乗算することにより、温度断層分布{ΔT}を推定値として算出した。

【0091】
ただし、ムーアペンローズ一般逆行列により算出された推定値は規格化されており、物理量として単位をもたない値となる。そこで、この推定値を有用な物理量の温度値に較正する。すなわち、図9(A)に点線枠にて示すように、推定値の一部に対応する温度の真値を取得し、その真値に対応する推定値(「対応推定値」ともいう)とその真値とに基づいて、全推定値を線形補間することにより較正を行う。この真値は実際に計測した値(実計測値)であってよい。

【0092】
すなわち、逆解析にて得た温度分布を、下記式(26)に基づいて較正する。
【数26】
JP0006707249B2_000027t.gif
ここで、T*は逆解析により算出される推定値であり、T*realizedは較正後の温度値(実値)である。図9(B)に概念的に示すように、Tsample1,Tsample2は、較正に用いる任意の2点における真値であり、T*sample1,T*sample2はその真値に対応する対応推定値である。なお、本実施例では線形(一次関数)による較正を例示しているが、多項式にて較正(多項式近似)してもよい。

【0093】
制御演算部6は、以上のような演算処理を実行し、算出された温度値T*realizedを表示装置8に断層可視化する。それにより、例えば図9(C)に示すような色分けされた温度分布が表示可能となる。

【0094】
次に、制御演算部6が実行する具体的処理の流れについて説明する。
図10は、制御演算部6により実行されるマイクロ断層可視化処理の流れを示すフローチャートである。なお、この断層可視化処理は、上述したOCTによる断層撮影がなされることを前提に実行される。ここでは、測定対象Wが熱変形したときの変形前後の断層画像を取得し、物理量として測定対象Wの内部の温度分布を可視化する場合を例に説明する。

【0095】
制御演算部6は、OCTにより撮影された前後2枚のOCT断層画像を読み込む(S10)。続いて、それらの断層画像に基づいて変形ベクトル分布を算出するための変形ベクトル演算処理を実行する(S12)。続いて、その変形ベクトル分布を用いて逆解析処理を実行する(S14)。そして、その逆解析により算出された物理量の断層分布(温度分布)を表示装置8に表示させる(S16)。

【0096】
図11は、図10におけるS12の変形ベクトル演算処理を詳細に示すフローチャートである。制御演算部6は、変形前後の断層画像に基づき、再帰的相互相関法による処理を実行する。ここではまず、最小解像度(最大サイズの検査領域)での相互相関処理を実行し、相関係数分布を求める(S20)。続いて、隣接相互相関乗法により、隣接する相関係数分布の積を演算する(S22)。このとき、標準偏差フィルタ等の空間フィルタにより過誤ベクトルの除去をし(S24)、最小二乗法等により除去ベクトルの内挿補間を実行する(S26)。続いて、検査領域を小さくすることによって解像度を上げて相互相関処理を継続する(S28)。すなわち、低解像度での参照ベクトルを基に相互相関処理を実行する。このときの解像度が予め定める最高解像度でなければ(S30のN)、S22に戻る。

【0097】
そして、S22~S28の処理を繰り返し、最高解像度での相互相関処理が完了すると(S30のY)、サブピクセル解析を実行する。すなわち、最高解像度(最小サイズの検査領域)での変形ベクトルの分布に基づき、風上勾配法によるサブピクセル移動量を演算する(S32)。そして、このとき算出されたサブピクセル移動量に基づき、画像変形法によるサブピクセル変形量を演算する(S34)。続いて、最大相互相関値によるフィルタ処理により過誤ベクトルの除去をし(S36)、最小二乗法等により除去ベクトルの内挿補間を実行する(S38)。このようにして変形ベクトル分布が得られる。

【0098】
図12は、図10におけるS14の逆解析処理を詳細に示すフローチャートである。制御演算部6は、演算対象となった断層画像に対応するROIを設定する(S40:図7(A)参照)。この設定は、当該解析を行うユーザの入力指示にしたがって行われる。ユーザは、測定対象Wにおいて温度分布を表示させたい箇所をROIとして設定することができ、入力装置7を介してそれを指示することができる。

【0099】
続いて、制御演算部6は、設定されたROIに対して要素分割を設定する(S42:図7(B)参照)。この要素分割は、上述したFEMに対応するものであり、ユーザの入力指示にしたがって行われる。ユーザは、求める解析精度に応じた要素分割を入力装置7を介して指示することができる。

【0100】
続いて、制御演算部6は、S38を経て得られた変形ベクトル分布を、分割された要素の節点に補間する(S44:図8参照)。また、ROIに対して弾性率および線膨張係数を設定する(S46)。本実施例では、これらの値を測定対象Wの材料に固有の値として設定する。

【0101】
続いて、制御演算部6は、仮想仕事の原理に基づいて上記式(16)の全体剛性マトリクス[K]を演算する一方(S48)、形状関数,力学パラメータ(弾性率など),温度依存パラメータ(線膨張係数)から温度変化項の全体マトリクス観測行列[A]を決定する(S50)。そして、上述したマトリクスの縮小などを行い、逆解析(例えば逆行列演算)により温度分布を演算する(S52)。以上のようにして算出された温度分布が、図10のS16にて断層可視化される。

【0102】
次に、本実施例による効果について説明する。発明者は、本実施例による手法、つまりOCTによる測定とFEMの逆解析とによるハイブリッド手法の有効性を確認するために、数値実験による精度検証を行った。以下、その検証結果について説明する。図13~図15は、数値実験による検証方法を示す図である。各図の(A)~(C)はその検証過程を示す。図16は、検証結果を示す図である。(A)は本実施例のハイブリッド手法による解析結果を示し、(B)はFEMによる正解値を示す。図17は、本実施例の解析による正解値との誤差を示す図である。(A)は絶対誤差を示し、(B)は相対誤差を示す。

【0103】
上記ハイブリッド手法の精度を評価するためには、その比較対象となる正解値が必要となる。本来であれば、この比較対象を温度センサ等による実測値(真値)として得るのが望ましいが、物体内部の断層に沿ってそれを行うのは不可能に等しい。そこで、ここでは測定対象のモデルに対してFEMによる順解析を実行し、それにより得られる温度分布を正解値とする。そして、その正解値からOCT断層画像を疑似的に作成し、その断層画像(以下「疑似OCT断層画像」ともいう)をOCTによる実際の断層画像と仮定して上記逆解析を実行する。その逆解析によって得られた温度分布と、正解値としての温度分布(以下「正解温度分布」ともいう)とを比較することで、上記ハイブリッド手法の精度を評価した。

【0104】
(数値実験)
具体的には、モンテカルロシミュレーションによって散乱体からのランダムスッペクルパターンを演算することで、図14(B)に示す疑似OCT断層画像を作成し、FEMを用いて散乱体の熱変形移動量の算出を行うこととした。

【0105】
まず、数値モデルに基づいたFEMによる伝熱解析を行った。図13(A)に示すように、測定対象として基板にICチップ(シリコンチップ)が埋め込まれた電子実装部品を模擬したモデルを想定している。このモデルは、エポキシ樹脂(Epoxy)を基材とし、シリコンチップのシリコン(Silicon)、配線の銅(Copper)を内包し、表面にはレジスト材(Resist)が存在する複合材料とされている。

【0106】
そして、図13(B)に示すように、このモデルに対してFEMの順解析による要素分割を施した。なお、モデルを構成する各材料の材料定数は、下記の表1に示すとおりである。
【表1】
JP0006707249B2_000028t.gif

【0107】
温度の境界条件については、シリコンチップおよび配線が発熱していると仮定した。すなわち、シリコンおよび銅の温度を250[℃]とし、表面(上下面)であるz=0,2000[μm]の温度を外気温に等しい25[℃]とした。また、熱伝達係数を10[W/(mK)]とした。ただし、モデルの側面(x=0,5000[μm])では熱の出入りはないと仮定した。このような設定のもとFEMの伝熱解析を実行することで、図13(C)に示す温度分布を得た。この温度分布が数値実験における「正解温度分布」とされる。

【0108】
この温度分布データを用いてFEMによる力学解析を実行し、モデルを熱変形させることとした。その際、モデルの境界拘束条件として、z=0[μm]の位置を全方向に固定することとした。ICチップを内包する電子実装部品は、マザー基板へ実装される際にその底面がはんだ付けされることを考慮したものである。

【0109】
図14(A)には、このような連成解析により算出された熱変形ベクトル分布が示されている。この熱変形量データを図14(B)の各スペックルに対応させ、スペックルの移動を決定した。各スペックルの熱変形移動は、スペックルを内包する有限要素の節点変形ベクトルおよび形状関数を用いて算出することができる。それにより、図14(C)に示す熱変形後の疑似OCT断層画像を作成した。

【0110】
このようにして、図14(B)および(C)に示すように、熱変形前後の疑似OCT断層画像を得た。なお、ここでいう「熱変形前」とは、シリコンチップおよび配線が発熱する前の状態、つまり両者の温度が実質的に外気温である25[℃]に等しい状態を意味する。「熱変形後」とは、シリコンチップおよび配線が上述のように250[℃]に発熱した状態を意味する。熱変形前後の疑似OCT断層画像は、それぞれの状態を仮定して得たものである。この2つの疑似OCT断層画像の取得を出発点として上記ハイブリッド手法を適用した。すなわち、変形前後の疑似OCT断層画像をOCTによる実際の断層画像と仮定し、図10~図12に示した処理を実行した。そして、図10のS16にて得られた温度分布を、図13(C)に示した正解温度分布と比較することで精度を評価した。

【0111】
具体的には、図15(A)に示す疑似OCT断層画像に対し、図15(B)に示すようにROIを設定し、そのROIについて要素分割を行った。なお、このROI内の弾性率分布,線膨張係数分布を定義する必要があるが、弾性率および線膨張係数は材料に固有の値であるため、ここでは上記表1に示す材料定数を設定した。

【0112】
図15(C)には、前後の疑似OCT断層画像に基づいて算出された変形ベクトル分布が示されている。この変形ベクトル分布に基づいて上述した逆解析処理を実行し、図16(A)に示すように温度分布を断層可視化した。この温度分布を図16(B)に示す正解温度分布と比較すると、それらが定性的に一致していることが分かる。すなわち、いずれも基材中心部分の温度がその周囲より低くなっているのに対し、シリコンチップおよび銅の付近で周囲より高い温度分布が得られていることが確認できる。

【0113】
図17(A)には、本実施例のハイブリッド手法により得られた温度分布(「推定温度分布」ともいう)と、正解温度分布との絶対誤差が示されている。図17(B)には両者の相対誤差が示されている。なお、「絶対誤差」とは両者の温度差を意味し、「相対誤差」とは絶対誤差を正解値により除算した値を意味する。

【0114】
これらの結果に基づき、両者の誤差平均(RMS誤差)は14%程度と比較的小さく収まることが確認できた。なお、ROIの境界近傍でエラーが大きくなる傾向にあるが、これは境界でのノイズの影響によるものと考えられる。なお、このような温度推定の誤差は、温度昇降を動的に検出してノイズ低減を施す方法や、ノイズを考慮した逆推定手法であるカルマンフィルタなどの導入により改善可能である。

【0115】
以上に説明したように、本実施例によれば、OCTによる断層計測(実計測データ)を用いることにより、材料内部の温度断層分布を逆解析することが可能となる。すなわち、実測と解析とのハイブリッドにより物理量を断層可視化できる。物体内部における任意の領域の温度分布を、境界条件を考慮することなく断層可視化できることから、特に境界条件が不明確な状況における内部温度の検出に好適となる。また、このような温度断層分布がマイクロスケールで非侵襲・非破壊にて検出可能であるため、MEMS等の電子実装部品のように境界条件の付与が難しい工業製品への応用が期待される。

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

【0117】
上記実施例では、OCT断層画像として、構造の変形前後の静止画を用いる例を示した。変形例においては、構造の変形中を撮影したOCT断層画像の動画に対して上記手法を適用してもよい。すなわち、動画中の異なる時刻の2枚のOCT断層画像を取得し、それらに対して上記逆解析を適用して物理量を演算し、断層可視化してもよい。上記ハイブリッド手法は、このように時間変化する現象に対しても適用可能である。

【0118】
上記実施例では述べなかったが、測定対象の温度場が変化しているような非定常な状況では、変形している場(分布)が時々刻々と変化する。そのような状況下においてOCT断層画像を連続撮影して本手法を適用することにより、時々刻々変化する非定常温度場(分布)をマイクロ断層可視化することが可能となる。この演算過程において、連続するOCT断層画像と時空間最小二乗法を用いた節点補間を行うなどが考えられる。

【0119】
上記実施例では、変形ベクトルの断層分布に基づいて逆解析を行い、それにより温度分布を算出する例を示した。変形例においては、変形ベクトルの時間微分である変形速度ベクトルの断層分布に基づいて逆解析を行い、それにより物理量の分布を算出してもよい。

【0120】
上記実施例では、OCTによる断層計測と数値解析の逆解析とにより部品内部の温度断層分布を可視化する例を示した。変形例においては、同様の手法により部品内部の弾性係数や粘弾性係数の分布を可視表示することも可能である。例えば、上記式(20)に関して温度変化{ΔT}がゼロであるなど既知の場合、測定項{U}の取得により全体剛性マトリクス[K]を算出し、弾性係数の分布を可視表示することができる。それにより、部品内部に存在する材料を特定することも可能となる。

【0121】
上記実施例では、OCTによる断層画像から得た変形ベクトル分布を逆解析することにより物理量を算出する例を示したが、OCTに代えてX線CT、超音波、MRIなどを採用することもできる。例えば、電子実装部品などの高分子基複合材や金属との複合材に適用する際にはマイクロX線CTによる断層画像などが十分適用できると考えられる。

【0122】
上記実施例では、OCTによる断層画像を二次元で取得する例を示したが、三次元で取得してもよい。すなわち、変位関連ベクトルとしての変形ベクトルあるいは変形速度ベクトルを三次元データとして処理してもよいことは言うまでもない。

【0123】
上記実施例では、物理量を断層可視化する対象として部品構造を例示した。変形例においては、その他の材料構造を対象としてもよい。また、体内の軟骨やしこり(硬化箇所)など、生体構造(生体組織)を対象としてもよい。生体構造の力学特性の変化を断層可視化することで、例えば動脈硬化,癌,肝硬変,皮膚,軟骨,骨などの多くの症例に対する診断材料とすることが可能となる。

【0124】
なお、本発明は上記実施例や変形例に限定されるものではなく、要旨を逸脱しない範囲で構成要素を変形して具体化することができる。上記実施例や変形例に開示されている複数の構成要素を適宜組み合わせることにより種々の発明を形成してもよい。また、上記実施例や変形例に示される全構成要素からいくつかの構成要素を削除してもよい。
【符号の説明】
【0125】
1 断層可視化システム、2 光学ユニット、4 光学機構、6 制御演算部、7 入力装置、8 表示装置、10 光源、12 オブジェクトアーム、14 リファレンスアーム、16 光検出器、26 反射鏡、ROI 関心領域、S1 検査領域、S2 探査領域。
図面
【図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