(57)【特許請求の範囲】
【請求項１】
構造を対象とする順解析の過程で変位関連ベクトルが演算される所定の数値解析において解析条件となり得る物理量に関し、前記構造の内部における前記物理量の分布を断層可視化するマイクロ断層可視化方法であって、
前記構造の変形前後の断層画像を取得する工程と、
前記断層画像に基づき変位関連ベクトルの断層分布を演算する工程と、
前記断層画像に対して関心領域を設定する工程と、
前記関心領域を要素分割する工程と、
演算された変位関連ベクトルを、分割された要素の節点に補間するとともに境界条件として与える工程と、
補間後の変位関連ベクトルの断層分布を用いて前記数値解析の逆解析を境界条件の設定を別途行わずに実行することにより、各要素における前記物理量の分布を前記物理量の断層分布として演算する工程と、
演算された前記物理量の分布を表示装置に断層可視化する工程と、
を備えることを特徴とするマイクロ断層可視化方法。
【請求項２】
前記関心領域が、前記断層画像の内部に対して設定されることを特徴とする請求項１に記載のマイクロ断層可視化方法。
【請求項３】
前記数値解析が、状態変化又は境界条件変化による要因に伴う前記構造の変形を考慮した応力解析であることを特徴とする請求項１または２に記載のマイクロ断層可視化方法。
【請求項４】
前記物理量が温度であり、
前記数値解析が有限要素法による熱変形を考慮した熱応力解析であることを特徴とする請求項３に記載のマイクロ断層可視化方法。
【請求項５】
前記物理量が前記構造の物理定数であり、
前記数値解析が有限要素法による応力解析であることを特徴とする請求項３に記載のマイクロ断層可視化方法。
【請求項６】
光コヒーレンストモグラフィーを用いて前記構造の断層画像を取得することを特徴とする請求項１〜５のいずれかに記載のマイクロ断層可視化方法。
【請求項７】
前記構造として、少なくとも高分子基材料又は生体組織を対象として含むことを特徴とする請求項１〜６のいずれかに記載のマイクロ断層可視化方法。
【請求項８】
前記構造における所定の領域において前記物理量の真値を与える工程と、
前記逆解析により得られた物理量のうち、前記領域に対応する物理量である推定値と、前記領域において与えられた真値との関係に基づいて、前記逆解析により得られた物理量の値を補正する工程と、
をさらに備えることを特徴とする請求項１〜７のいずれかに記載のマイクロ断層可視化方法。
【請求項９】
測定対象となる構造の内部における所定の物理量の分布を断層可視化するマイクロ断層可視化システムであって、
前記構造の変形前後の断層画像を取得するための断層画像取得装置と、
前記断層画像取得装置から前記構造の変形前後の断層画像データを取得し、その断層画像データに基づいて前記物理量の断層分布を演算する制御演算部と、
前記制御演算部の演算結果に基づいて、前記構造における前記物理量の分布を断層可視化する態様で表示する表示装置と、
を備え、
前記物理量は、前記構造を対象とする順解析の過程で変位関連ベクトルが演算される所定の数値解析において解析条件となり得るものであり、
前記制御演算部は、
取得した断層画像データに基づき変位関連ベクトルの断層分布を演算し、
前記断層画像に対して関心領域を設定し、
前記関心領域を要素分割し、
演算された変位関連ベクトルを、分割された要素の節点に補間するとともに境界条件として与え、
補間後の変位関連ベクトルの断層分布を用いて前記数値解析の逆解析を境界条件の設定を別途行わずに実行することにより、各要素における前記物理量の分布を前記物理量の断層分布として演算し、
演算された前記物理量の分布を前記表示装置に断層可視化させることを特徴とするマイクロ断層可視化システム。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、構造内部における物理量のマイクロ断層可視化を可能にする方法およびシステムに関する。
【背景技術】
【０００２】
パーソナルコンピュータや携帯端末などの電子機器には、いわゆるＭＥＭＳ(Micro Electro Mechanical Systems)と呼ばれる電子実装部品が組み込まれている。このような電子実装部品の小型化，高密度化に伴うジュール発熱密度の増大により、熱ひずみによる材料内部の微視的破壊が問題となっている。特に、導電性フィラーを含有する繊維強化プラスチックなど（CFRP，GFRP等）、温度依存導電性を有する高分子基複合材料からなる部品について、その問題が顕著となっている。また、例えばリチウムイオン二次電池用キャパシタは、薄板電極間における温度上昇がドライアップ故障を引き起こすといった問題を有する。このため、これらの工業製品にはマイクロスケールでの最適熱設計が必要とされる。
【０００３】
このような工業製品の熱設計においては、構造内部の温度評価が重要となるが、内部の温度情報が所望の位置でマイクロスケールにて得られる実計測システムは未だ存在せず、有限要素法等の数値解析に頼るところが大きい（例えば特許文献１参照）。また、設計という観点からは、このような熱設計だけでなく、強度（剛性）など力学特性を考慮した設計も重要となるが、これも同様の数値解析に頼るところが大きい。さらに、工業製品を離れると、例えば生体組織ではその力学物理量を測定することが困難であるため、これも同様の数値解析に依存する。
【先行技術文献】
【特許文献】
【０００４】
【特許文献１】特開２００６−２８４２４９号公報
【発明の概要】
【発明が解決しようとする課題】
【０００５】
しかしながら、このような熱解析や力学解析などの数値解析をマイクロスケールにて行う場合、それらにおける境界条件の与え方が難しい。また、局所的なデータの取得を目的とする場合にも部品全体を解析対象とする必要があり、演算効率上の問題も生じる。
【０００６】
本発明はこのような課題に鑑みてなされたものであり、その目的の一つは、構造内部における物理量を簡易な手法にてマイクロ断層可視化可能な技術を提供することにある。
【課題を解決するための手段】
【０００７】
本発明のある態様は、構造を対象とする順解析の過程で変位関連ベクトルが演算される所定の数値解析において解析条件又は解析結果となり得る物理量に関し、構造の内部における物理量の分布を断層可視化するマイクロ断層可視化方法である。この方法は、構造の変形前後の断層画像を取得する工程と、断層画像に基づき変位関連ベクトルの断層分布を演算する工程と、演算された変位関連ベクトルの断層分布を用いて数値解析の逆解析を実行することにより、物理量の断層分布を演算する工程と、演算された物理量の分布を表示装置に断層可視化する工程と、を備える。なお、ここでいう「変位関連ベクトル」は、変形ベクトルであってもよいし、変形速度ベクトルであってもよい。
【０００８】
この態様によると、実測により得られる断層画像から変位関連ベクトルの分布が演算され、その変位関連ベクトル分布に基づき所定の数値解析の逆解析を行うことで物理量の分布が得られる。すなわち、実測と数値解析の逆解析とのハイブリッドにより物理量を断層可視化することが可能となる。実測により内部の変位関連ベクトルが得られるため、数値解析の順解析のように境界部に境界条件を設定する必要もない。このため、所望の断層位置について物理量を簡易に算出することができる。すなわち、この態様によれば、構造内部における物理量を簡易な手法にてマイクロ断層可視化することができる。
【０００９】
本発明の別の態様は、測定対象となる構造の内部における所定の物理量の分布を断層可視化するマイクロ断層可視化システムである。このシステムは、構造の変形前後の断層画像を取得するための断層画像取得装置と、断層画像取得装置から構造の変形前後の断層画像データを取得し、その断層画像データに基づいて物理量の断層分布を演算する制御演算部と、制御演算部の演算結果に基づいて、構造における物理量の分布を断層可視化する態様で表示する表示装置と、を備える。
【００１０】
物理量は、構造を対象とする順解析の過程で変位関連ベクトルが演算される所定の数値解析において解析条件又は解析結果となり得るものである。なお、ここでいう「変位関連ベクトル」は、変形ベクトルであってもよいし、変形速度ベクトルであってもよい。制御演算部は、取得した断層画像データに基づき変位関連ベクトルの断層分布を演算し、演算された変位関連ベクトルの断層分布を用いて数値解析の逆解析を実行することにより、物理量の断層分布を演算し、演算された物理量の分布を表示装置に断層可視化させる。
【００１１】
この態様によると、断層画像取得装置にて取得される断層画像から変位関連ベクトルの分布が演算され、その変位関連ベクトル分布に基づき所定の数値解析の逆解析を行うことで物理量の分布が得られる。すなわち、断層画像取得装置による実測と、制御演算部による数値解析の逆解析とのハイブリッドにより、物理量を断層可視化することが可能となる。実測により内部の変位関連ベクトルが得られるため、数値解析の順解析のように境界部に境界条件を設定する必要もない。このため、所望の断層位置について物理量を簡易に算出できる。すなわち、この態様によれば、構造内部における物理量を簡易な手法にてマイクロ断層可視化することができる。
【発明の効果】
【００１２】
本発明によれば、構造内部における物理量を簡易な手法にてマイクロ断層可視化可能な技術を提供することができる。
【図面の簡単な説明】
【００１３】
【図１】実施例に係る断層可視化システムの構成を概略的に表す図である。
【図２】実施例に係るマイクロ断層可視化方法の概念を表す図である。
【図３】再帰的相互相関法による処理手順を概略的に示す図である。
【図４】サブピクセル解析による処理手順を概略的に示す図である。
【図５】有限要素モデルを概略的に示す図である。
【図６】逆解析の処理手順を概略的に示す図である。
【図７】逆解析の処理手順を概略的に示す図である。
【図８】逆解析の処理手順を概略的に示す図である。
【図９】逆解析の処理手順を概略的に示す図である。
【図１０】制御演算部により実行される物理量可視化処理の流れを示すフローチャートである。
【図１１】図１０におけるＳ１２の変形ベクトル演算処理を詳細に示すフローチャートである。
【図１２】図１０におけるＳ１４の逆解析処理を詳細に示すフローチャートである。
【図１３】数値実験による検証方法を示す図である。
【図１４】数値実験による検証方法を示す図である。
【図１５】数値実験による検証方法を示す図である。
【図１６】検証結果を示す図である。
【図１７】本実施例の解析による正解値との誤差を示す図である。
【発明を実施するための形態】
【００１４】
本発明の一実施形態は、構造内部の物理量の分布をマイクロスケールにて断層可視化する方法である。なお、ここでいう「構造」は、材料構造や生体構造（生体組織）を含み得る。「断層可視化」は、少なくともマイクロスケールにて可能であるが、それ以外のスケールを排除するものではない。断層画像取得装置の解像度に応じてナノスケールにて可能なものも含み得る。「物理量」は、所定の数値解析において解析条件又は解析結果の一つとなり得るものであればよく、後述する温度分布など構造の状態変化に係るものでもよいし、物理定数など構造の物性値に係るものでもよい。その場合、「物理定数」は、「弾性係数」，「粘性係数」，「透水係数」等を含み得る。「断層画像」は、光コヒーレンストモグラフィー（Optical Coherence Tomography：以下「ＯＣＴ」という）によるものでもよい。あるいは、Ｘ線ＣＴや超音波断層像によるものでもよい。
【００１５】
この方法は、数値解析の順解析の過程で変位関連ベクトルの分布が演算される点と、断層画像取得装置により実測される断層画像から変位関連ベクトル分布を演算可能である点に着目し、両者の変位関連ベクトルを結び付ける。ここで、「変位関連ベクトル」は、構造内部の変位や変形に関連するベクトルであり、変形ベクトルであってもよいし、変形速度ベクトルであってもよい。「変形ベクトル」は「変位ベクトル」の概念を含む。すなわち、実測と数値解析とを後者の「逆解析」という形で結び付け、両者のハイブリッドによる演算処理を実現している。ここでいう「逆解析」は、順解析の一部を逆方向にたどることにより実現されてよい。
【００１６】
すなわち、順解析においては、例えば上記物理量が所定の数理モデルの解析条件として入力され、その数理モデルに沿った演算処理を経ることで解が算出される。その演算過程で変位関連ベクトルが演算される。「逆解析」は、逆にその変位関連ベクトルの分布を既知として順解析を逆方向にたどり、解析条件の一つを構成する「物理量」の分布を算出する。その変位関連ベクトル分布として、実測された断層画像に基づいて算出されたものが用いられる。なお、「変形ベクトル」の時間微分が「変形速度ベクトル」であるが、「変形速度ベクトル」を変形ベクトルの時間変化とみて「変形ベクトル」の概念で捉えてもよい。いずれにしても、演算過程において変形ベクトルが演算される。本実施形態によれば、ＯＣＴやＣＴ等の断層可視化手法から逆解析に結び付ける斬新な手法が提供される。
【００１７】
この方法では、断層画像取得装置を介して構造の変形前後の断層画像が取得される。なお、「変形前後の断層画像」については、構造の変形前と変形後のそれぞれについて変化が停止している状態で取得してもよい。あるいは、変形中の前後の時点で取得してもよい。すなわち、本手法は、後者のように時間的に変化する現象にも適用可能である。具体的には、ＯＣＴ断層画像の動画中の２枚の断層画像に対して本手法を適用することにより、それらの断層画像に対応する温度分布を断層可視化することが可能である。前後の断層画像を連続的に取得すれば、温度分布をリアルタイムに断層可視化することもできる。
【００１８】
そして、その断層画像に基づいて断層位置に対応した変位関連ベクトルの分布が演算される一方、断層画像の内部に対して関心領域（Region of Interest：以下「ＲＯＩ」と表記する）が設定される。また、上記数値解析の逆解析を利用するためにＲＯＩに対して要素分割がなされる。これらＲＯＩおよび要素分割の設定は、入力装置を介したユーザの指示入力によりなされる。
【００１９】
断層画像に基づいて演算された変位関連ベクトルは、分割された要素の節点に補間される。この補間後の変形ベクトルの分布に対して上記逆解析を施すことにより、各要素における物理量の分布が演算される。そして、このとき演算された物理量の分布が、表示装置に断層可視化される。
【００２０】
上記数値解析は、状態変化又は境界条件変化による要因に伴う構造の変形を考慮した応力解析であってもよい。例えば、有限要素法（Finite Element Method：以下「ＦＥＭ」と表記する）、有限差分法（Finite Difference Method：以下「ＦＤＭ」と表記する）、あるいは境界要素法（Boundary Element Method：以下「ＢＥＭ」と表記する）によるものでもよい。
【００２１】
上記手法は、構造における所定の領域において物理量の真値を与える工程と、逆解析により得られた物理量のうち上記領域に対応する物理量である推定値と、上記領域において与えられた真値との関係に基づき、逆解析により得られた物理量の値を補正する工程と、をさらに備えるものでもよい。
【００２２】
「所定の領域」は、複数設定されてもよく、構造の表面又はその近傍に設定されてもよいし、構造内部の特定の部材の周囲に沿って設定されてもよい。「真値」は、実測値であってもよいし、解析結果（解析値）であってもよく、正しいと推定される値であればよい。例えば、「物理量」が温度であり、赤外線サーモグラフィーや熱電対によりＲＯＩ近傍の温度値を計測できる場合、その実測値である温度値を真値としてもよい。あるいは、銅配線の電流と電圧によって算出されるジュール熱などに基づいてＲＯＩ近傍の計測データを算出できる場合、その演算結果を真値としてもよい。この「真値」は、ＲＯＩ内部の値であることが望ましいが、ＲＯＩ近傍の温度値（実験計測値）であってもよい。また、熱伝導率などを仮定してモデリングし、ＲＯＩ内部の任意領域の温度値を推定し、その推定値を「真値」としてもよい。
【００２３】
上記数値解析は、状態変化又は境界条件変化による要因に伴う構造の変形を考慮した応力解析であってもよい。ここでいう「要因」は、構造への通電であってもよいし、構造の周囲温度の変化であってもよい。あるいは、構造に外力が付与されることでもよい。
【００２４】
例えば、力学的物理量のマイクロ断層可視化に際して下記式（ａ）を利用する。
[Ｋ]{Ｕ}＝{Ｆ}  ・・・（ａ）
ここで、[Ｋ]は全体剛性マトリクスであり、{Ｕ}は要素分割されたＲＯＩの各節点における変形ベクトルであり、{Ｆ}は外力ベクトルである。上記式（ａ）において、{Ｆ}に関し、ＲＯＩ内部の節点の荷重はゼロ（合力ゼロ）であり、ＲＯＩ境界の節点の荷重に対しては変位境界条件（変位拘束条件）を与える。{Ｕ}については、変形前後の断層画像に基づいて演算される変位関連ベクトル分布（変形ベクトルの断層分布等）から得られるため、それを代入することができる。上記式（ａ）を逆解析することにより、物理量としての[Ｋ]を推定することができる。
【００２５】
以下、図面を参照しつつ本実施形態を具体化した実施例について詳細に説明する。
［実施例］
図１は、実施例に係る断層可視化システムの構成を概略的に表す図である。本実施例の断層可視化システムは、電子実装部品の内部の特定の物理量をマイクロスケールにて断層可視化するものである。この断層可視化にＯＣＴを利用する。
【００２６】
図１に示すように、断層可視化システム１は、ＯＣＴを用いる光学系を含む光学ユニット２、光学ユニット２に接続される光学機構４、ＯＣＴにより得られた光干渉データに基づいて物理量を演算する制御演算部６、ユーザの指示入力を受け付ける入力装置７、および物理量を断層可視化する態様で表示する表示装置８を備える。光学ユニット２および光学機構４は、「断層画像取得装置」を構成する。図示の例では、光学ユニット２としてマッハツェンダー干渉計をベースとした光学系が示されているが、マイケルソン干渉計その他の光学系を採用することもできる。
【００２７】
本実施例では、ＯＣＴとして波長走査型レーザを光源とするＳＳ−ＯＣＴ（Swept Source OCT）を用いるが、ＴＤ−ＯＣＴ（Time Domain OCT）、ＳＤ−ＯＣＴ（Spectral Domain OCT）その他のＯＣＴを用いることもできる。ＳＳ−ＯＣＴは、ＴＤ−ＯＣＴとは異なり、参照鏡走査などの機械的な光遅延走査を必要としないため、高い時間分解能と高い位置検出精度が得られる点で好ましい。
【００２８】
光学ユニット２は、光源１０、オブジェクトアーム１２、リファレンスアーム１４、および光検出器１６を備える。各光学要素は、光ファイバーにて互いに接続されている。光源１０から出射された光は、カプラ１８（ビームスプリッタ）にて分けられ、その一方がオブジェクトアーム１２に導かれる物体光となり、他方がリファレンスアーム１４に導かれる参照光となる。オブジェクトアーム１２に導かれた物体光は、サーキュレータ２０を介して光学機構４に導かれ、測定対象Ｗに照射される。この物体光は、測定対象Ｗの表面および断面にて後方散乱光として反射されてサーキュレータ２０に戻り、カプラ２２に導かれる。
【００２９】
一方、リファレンスアーム１４に導かれた参照光は、サーキュレータ２４を介して反射鏡２６に導かれる。このとき、参照光は、コリメートレンズ２８を経て集光レンズ３０にて反射鏡２６に集光される。この参照光は、反射鏡２６にて反射されてサーキュレータ２４に戻り、カプラ２２に導かれる。すなわち、物体光と参照光とがカプラ２２にて合波（重畳）され、その干渉光が光検出器１６により検出される。
【００３０】
光学機構４は、光学ユニット２のオブジェクトアーム１２を構成し、サーキュレータ２０から延びる光ファイバー３４が接続されている。光学機構４は、光学ユニット２からの光を測定対象Ｗに導いて走査させる機構と、その機構を駆動するための駆動部（アクチュエータ）を備える。光学機構４は、例えばガルバノ鏡を含むガルバノ装置であってもよい。光学機構４は、測定対象Ｗに向けて物体光を照射する。
【００３１】
測定対象Ｗは、ＭＥＭＳに対応した電子実装部品であり、ＩＣチップを基板に埋め込んで構成される。その基板は、高分子樹脂材料等の高分子基材料やプラスチック等の結晶性材料からなるものでもよい。その基板は、内部に反射する構造（界面など）がある材料からなるものであればよく、ナノ粒子などの散乱体が樹脂に内在する複合材料など、透明でない材料からなるものでもよい。
【００３２】
カプラ２２にて合波された干渉光は、光検出器１６に入力される。光検出器１６は、これを光干渉信号（干渉光強度を示す信号）として検出する。この光干渉信号は、Ａ／Ｄ変換器４０を介して制御演算部６に入力される。Ａ／Ｄ変換器４０は、光検出器１６から出力されたアナログ信号をデジタル信号に変換して制御演算部６へ出力する。
【００３３】
制御演算部６は、ＣＰＵ、ＲＯＭ、ＲＡＭ、ハードディスクなどを有し、これらのハードウェア、およびソフトウェアによって、光学ユニット２の光学系全体の制御と、光学機構４の駆動制御と、ＯＣＴによる画像出力のための演算処理を行う。制御演算部６の指令信号は、Ｄ／Ａ変換器４２を介して光学機構４等に入力される。Ｄ／Ａ変換器４２は、制御演算部６から出力されたデジタル信号をアナログ信号に変換して光学機構４等へ出力する。制御演算部６は、光学機構４の駆動に基づいて光学ユニット２から出力された光干渉信号を処理し、ＯＣＴによる測定対象Ｗの断層画像を取得する。そして、その断層画像データに基づき、後述の手法により測定対象Ｗの内部における特定の物理量の断層分布を演算する。
【００３４】
入力装置７は、キーボードやタッチパネル等からなり、ユーザの指示入力を受け付け、制御演算部６に入力する。表示装置８は、例えば液晶ディスプレイからなり、制御演算部６にて演算された測定対象Ｗの内部の物理量を断層可視化する態様で画面に表示させる。
【００３５】
以下、上記物理量の断層分布を可視表示するまでの演算処理方法について説明する。  
図２は、実施例に係るマイクロ断層可視化方法の概念を表す図である。本実施例では、測定対象Ｗの内部の温度を上記物理量としてその分布を断層可視化する。この手法は、測定対象Ｗの熱変形に基づいて内部の温度分布を演算するものであり、ＯＣＴにて取得した断層画像から変形ベクトル分布を演算する工程（変形ベクトル演算工程）と、その変形ベクトル分布に対して有限要素法（以下「ＦＥＭ」と表記する）の逆解析を施すことにより、測定対象Ｗの内部の温度分布を演算する工程（逆解析工程）とを有する。
【００３６】
まず、変形ベクトル演算工程について説明する。上述のように、ＯＣＴにおいて、オブジェクトアーム１２を経た物体光（測定対象Ｗからの反射光）と、リファレンスアーム１４を経た参照光とが合波され、光検出器１６により光干渉信号として検出される。制御演算部６は、この光干渉信号を干渉光強度に基づく測定対象Ｗの断層画像として取得することができる。
【００３７】
ところで、ＯＣＴの光軸方向（奥行き方向）の分解能であるコヒーレンス長ｌｃは、光源の自己相関関数によって決定される。ここでは、コヒーレンス長ｌｃを自己相関関数の包括線の半値半幅とし、下記式（１）にて表すことができる。
【数１】
0006707249.tif 000002
ここで、λｃはレーザビームの中心波長であり、Δλはレーザビームの半値全幅である。
【００３８】
一方、光軸垂直方向（ビーム走査方向）の分解能は、集光レンズによる集光性能に基づき、ビームスポット径Ｄの１／２とされる。そのビームスポット径Ｄは、下記式（２）にて表すことができる。
【数２】
0006707249.tif 000003
ここで、ｄＢは集光レンズに入射するビーム径、ｆは集光レンズの焦点である。このようにＯＣＴによる分解能には限界があるところ、本実施例では後述するサブピクセル解析の導入などにより、温度分布の断層表示をマイクロスケールにて行うことを可能にしている。以下、その詳細について説明する。
【００３９】
まず、ＯＣＴを利用したマイクロスケールの力学特性検出法について説明する。この検出法は、測定対象Ｗの変形前後の２枚のＯＣＴ断層画像にデジタル相互相関法を適用することで変形ベクトル分布を算出する。この手法によれば、マイクロスケールにて測定対象Ｗのひずみ断層分布を検出することができる。
【００４０】
変形ベクトル分布を算出する際には、繰り返し相互相関処理を実施する再帰的相互相関法(Recursive Cross-correlation method)を適用する。これは、低解像度において算出された変形ベクトルを参照し、探査領域を限定するとともに階層的に検査領域を縮小して相互相関法を適用する手法である。これにより、高解像度な変形ベクトルを取得することができる。さらに、スペックル・ノイズ低減法として、隣接検査領域の相関値分布との乗算を行う隣接相互相関乗法(Adjacent Cross-correlation Multiplication)を用いる。そして、乗算されることによって高ＳＮ化した相関値分布から最大相関値を探索する。
【００４１】
また、マイクロスケールにおける微小変形解析では、変形ベクトルのサブピクセル精度が重要となる。このため、輝度勾配を利用する風上勾配法(Up-stream Gradientmethod)と、伸縮および剪断を考慮した画像変形法(Image Deformation method)の両サブピクセル解析法を併用し、変形ベクトルの高精度検出を実現する。なお、ここでいう「風上勾配法」は、勾配法（オプティカルフロー法）の一種である。
【００４２】
このようにして得られた変形ベクトル分布を空間微分することにより、ひずみ断層分布を演算することができる。なお、本実施例では後述のように、変形ベクトル分布に対してＦＥＭの逆解析を施すことにより、測定対象Ｗの内部の温度分布を演算する。このため、ひずみ断層分布まで演算する必要はなく、その演算過程である変形ベクトル分布の演算までが利用されることとなる。
【００４３】
以下、各手法について詳細に説明する。図３は、再帰的相互相関法による処理手順を概略的に示す図である。図４は、サブピクセル解析による処理手順を概略的に示す図である。
【００４４】
（再帰的相互相関法）
図３（Ａ）〜（Ｃ）は、再帰的相互相関法による処理過程を示している。各図にはＯＣＴにより撮影される前後の断層画像が示されている。左側には先の断層画像（Image1）が示され、右側には後の断層画像（Image2）が示されている。
【００４５】
相互相関法とは、局所的なスペックル・パターンの類似度を下記式（３）に基づく相関値Ｒｉｊより評価する方法である。そのため、図３（Ａ）に示すように、撮影された前後のＯＣＴ画像について、先の断層画像（Image1）に類似度の検査対象となる検査領域Ｓ１が設定され、後の断層画像（Image2）に類似度の探査範囲となる探査領域Ｓ２が設定される。
【数３】
0006707249.tif 000004
ここで、空間座標として、光軸方向にＺ軸、光軸と垂直方向にＸ軸を設定している。ｆ(Ｘｉ,Ｚｊ)とｇ(Ｘｉ,Ｚｊ)は、変形前後のＯＣＴ画像において設置された中心位置(Ｘｉ,Ｚｊ)の検査領域Ｓ１(Ｎｘ×Ｎｚピクセル)のスペックル・パターンを表している。
【００４６】
探査領域Ｓ２(Ｍｘ×Ｍｚピクセル)内における相関値分布Ｒi,j(ΔＸ,ΔＺ)を算出し、図３（Ｂ）に示すようにパターンマッチングを行う。実際には、下記式（４）に示すように、最大相関値を与える移動量Ｕｉ,jを変形前後の変形ベクトルとして決定する。
【数４】
0006707249.tif 000005
なお、ｆ−とｇ−はｆ(Ｘｉ,Ｚｊ)とｇ(Ｘｉ,Ｚｊ)の検査領域Ｓ１内の平均値を表している。
【００４７】
本手法では、検査領域Ｓ１を縮小しながら相互相関処理を繰り返して空間解像度を高める再帰的相互相関法を採用している。なお、本実施例では解像度を上げる際に空間解像度が倍になるようにしている。図３（Ｃ）に示すように、検査領域Ｓ１を１／４に分割し、前階層にて算出された変形ベクトルを参照し、探査領域Ｓ２を縮小している。ここで探査領域Ｓ２も１／４に分割している。再帰的相互相関法を用いることで、高解像度において多発する過誤ベクトルの抑制を可能にしている。このような再帰的相互相関処理を施すことにより、変形ベクトルの解像度を高めることができる。
【００４８】
また、これに加え、下記式（５）により、演算中の座標を中心とする周囲８つの座標を含む合計９つの変形ベクトルの平均偏差σを用いた閾値を設定し、過誤ベクトルの除去を行い再帰処理に伴う誤差伝播を抑制する。
【数５】
0006707249.tif 000006
ここで、Ｕmはベクトル量の中央値を表し、閾値となる係数κは任意に設定した。
【００４９】
（隣接相互相関乗法）
本実施例では、スペックルノイズの影響を受けたランダム性の強い相関値分布から正確な最大相関値を決定する方法として、隣接相互相関乗法を導入している。下記式（６）により、隣接相互相関乗法では検査領域Ｓ１における相関値分布Ｒi,j(Δx,Δz)と、その検査領域Ｓ１にオーバーラップする隣接検査領域に対する相関値分布Ｒi＋Δi,j(Δx,Δz)とＲi,j＋Δj(Δx,Δz)の乗算を行い、新たな相関値分布Ｒ'i,j(Δx,Δz)を用いて最大相関値を検索する。
【数６】
0006707249.tif 000007
【００５０】
これにより、相関値同士の乗算によってランダム性を低減させることが可能になる。上述した検査領域Ｓ１の縮小と共に干渉強度分布の情報量も減少するため、スペックル・ノイズを原因とする複数相関ピークの出現が計測精度の悪化を招いていると考えられる。一方、隣接境界同士の移動量には相関があるため、最大相関値座標付近では強い相関値が残存する。この隣接相互相関乗法の導入によって最大相関値ピークが明瞭化され計測精度が向上し、正確な移動座標を抽出することが可能となる。また、この隣接相互相関乗法をＯＣＴの各ステージに導入することで、誤差伝播が抑制され、スペックル・ノイズに対するロバスト性が向上する。それにより、高空間解像度においても高精度な変形ベクトル分布（移動量分布）の算出が可能となる。
【００５１】
（風上勾配法）
図４（Ａ）〜（Ｃ）は、サブピクセル解析による処理過程を示している。各図にはＯＣＴにより連続的に撮影される前後の断層画像が示されている。左側には先の断層画像（Image1）が示され、右側には後の断層画像（Image2）が示されている。
【００５２】
本実施例では、サブピクセル解析のために風上勾配法と画像変形法を採用する。最終的な移動量の算出は後述の画像変形法によるが、計算の収束性の問題から、画像変形法より先に風上勾配法を適用する。検査領域サイズが小さく高空間解像度の条件において、サブピクセル移動量を高精度検出する画像変形法及び風上勾配法を適用している。画像変形法におけるサブピクセル移動量の検出が困難な場合において、風上勾配法によりサブピクセル移動量を算出する。
【００５３】
サブピクセル解析では、注目点における変形前後の輝度差が各成分の輝度勾配と移動量によって表される。このため、検査領域Ｓ１内の輝度勾配データより最小二乗法を用いてサブピクセル移動量を決定することができる。本実施例では、輝度勾配を求める際に、サブピクセル変形前の風上側の輝度勾配を与える風上差分法を採用している。すなわち、サブピクセル解析は様々な手法が存在するが、本実施例では検査領域サイズが小さく高空間解像度の条件においても、サブピクセル移動量を高精度検出する勾配法を採用している。
【００５４】
風上勾配法は、検査領域Ｓ１内の注目点の移動を、図４（Ａ）に示すピクセル精度に留まらず、図４（Ｂ）に示すサブピクセル精度にて算出するものである。なお、図中の各格子は１ピクセルを表している。実際には図示の断層画像と比較して相当小さいが、説明の便宜上、大きく表記している。この風上勾配法は、微小変形前後における輝度分布の変化を輝度勾配と移動量によって定式化する手法であり、ｆを輝度とすると、微小変形ｆ(x＋Δx,z＋Δz)をテイラー展開する下記式（７）として表される。
【数７】
0006707249.tif 000008
【００５５】
上記式（７）は、注目点の変形前後の輝度差が変形前の輝度勾配と移動量によって表されることを示している。なお、移動量(Δx,Δz)については上記式（７）のみでは決定できないため、検査領域Ｓ１内で移動量が一定と考え、最小二乗法を適用して算出している。
【００５６】
上記式（７）を用いて移動量を算出する際には、右辺の各注目点における移動前後の輝度差は一意にしか求まらない。そのため輝度勾配をどれだけ正確に算出するかが移動量の精度に直結する。輝度勾配の差分化では、一次精度風上差分を用いている。差分化において高次差分を適用すると、多くのデータが必要になり、ノイズが含まれていた際に影響を大きく受けてしまうためである。また、検査領域Ｓ１内の各点を基準とした高次差分では、検査領域Ｓ１外のデータを多く使用することとなり、検査領域Ｓ１そのものの移動量ではなくなってしまうという問題点も存在するからである。
【００５７】
輝度勾配を求める際に変形前の風上側の輝度勾配が移動することによって注目点の輝度差が生まれると考えることができるので、変形前は風上側の差分を適用する。ここでいう風上は、実際の移動方向ではなく、ピクセル移動量に対するサブピクセル移動量の向きのことであり、最大相関値ピークに放物線近似を施すことによって風上側を決定する。逆に、変形後の風下側の輝度勾配が逆に移動することによって注目点の輝度差が生じると考えることができるため、変形後は風下側の差分を適用する。
【００５８】
変形前の風上差分と変形後の風下差分を用いて２通りの解を求め、それらの平均をとった。さらに、実際には移動量が軸方向に沿わない場合には、変形前や変形後の輝度勾配が注目点と同一軸上に無く、ずれた位置の勾配を求める必要がある。つまり、Ｘ方向の勾配を求めたい場合には、Ｚ方向の移動も考慮して勾配を求めなければならない。そのため、輝度の内挿による輝度勾配の推定をすることで、精度向上を図っている。基本的には変形前（もしくは変形後）の位置を予測し、その位置での勾配を内挿により求める。
【００５９】
変形前（後）の注目点の位置は放物線近似を施した際のサブピクセル移動量により求める。その注目点位置が囲まれる８つの座標を用い、それらの比によって輝度勾配を算出する。具体的には、下記式（８）を用いる。そのようにして算出された輝度勾配と、輝度変化を用いて最小二乗法を適用し移動量を決定した。
【数８】
0006707249.tif 000009
【００６０】
（画像変形法）
上述した風上勾配法までは検査領域Ｓ１の形状は変更せず、正方形を保ったまま変形ベクトルの算出を行っている。しかし、現実には計測対象の変形に合わせて検査領域Ｓ１も変形していると考えられるため、検査領域Ｓ１の微小変形を考慮したアルゴリズムを導入し、変形ベクトル算出を高精度にて算出する必要がある。このため、本実施例ではサブピクセル精度での変形ベクトルの算出に画像変形法を導入している。すなわち、材料の変形前の検査領域Ｓ１と変形後の伸縮及びせん断変形を考慮した検査領域Ｓ１とで相互相関を実施し、相関値ベースの反復計算によってサブピクセル変形量を決定している。なお、検査領域Ｓ１の伸縮及びせん断変形は線形で近似している。
【００６１】
画像変形法は一般的に、材料の表面ひずみ計測法に用いられ、ランダムパターンを塗布した材料表面を高空間分解能なカメラで撮影した画像に対して適用される。一方、ＯＣＴ断層像はスペックルノイズを多く含み、また、複合材料での局所的な変形が発生するため、検査領域Ｓ１に変形が伴う場合には解の収束解が著しく低下する。本手法における検査領域Ｓ１の縮小は、局所的な力学特性を検出するために必要不可欠である。そこで画像変形法では、風上勾配法で得られた変形量を収束計算の初期値として採用し、さらに、輝度分布の双三次関数補間によって検査領域Ｓ１を縮小した際でも低ロバスト性を実現している。なお、変形例においては、双三次関数補間以外の補間関数を用いてもよい。
【００６２】
より詳細には、以下の手順にて演算を実行する。まず、材料変形前のＯＣＴ断層像の輝度分布に双三次関数補間法を適用し、輝度分布の連続化を実施する。双三次関数補間法とは、ｓｉｎｃ関数を区分的に三次関数近似した畳み込み関数を用い、輝度情報の空間連続性を再現する手法である。本来は連続的な輝度分布を画像計測する際には光学系に依存した点広がり関数が畳み込まれるため、ｓｉｎｃ関数を用いた逆畳み込みを行うことにより、本来の連続的な輝度分布が復元される。離散的な一軸信号ｆ(x)の補間をする場合、畳み込み関数ｈ(x)は下記式（９）にて表される。
【数９】
0006707249.tif 000010
【００６３】
なお、ＯＣＴ計測条件の違いによって輝度補間関数ｈ(x)の形状も変更する必要がある。そこで輝度補間関数ｈ(x)のｘ＝１での微係数ａを可変とし、ａの値を変更することで輝度補間関数ｈ(x)の形状を変更可能なアルゴリズムとした。本実施例では、擬似ＯＣＴ断層像を用いた数値実験による検証結果を元にし、ａの値を決定した。以上のように画像補間をすることで、伸縮及びせん断変形を考慮した検査領域Ｓ１の各点にて、ＯＣＴ輝度値を求めることが可能となる。
【００６４】
伸縮及びせん断変形を考慮して算出した検査領域Ｓ１は、図４（Ｃ）に示すように、移動とともに変形を伴う。材料変形前のＯＣＴ断層像におけるある検査領域Ｓ１内の整数ピクセル位置での座標(x,z)が変形後に座標(x＊,z＊)に移動すると考えると、x＊,z＊の値は下記式（１０）にて表される。
【数１０】
0006707249.tif 000011
【００６５】
ここで、Δ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は下記式（１１）にて表される。
【数１１】
0006707249.tif 000012
【００６６】
収束の判定には、反復計算で随時得られる漸近解が収束解の近傍で十分小さくなることを用いる。しかし、スペックルパターンの変化が激しい領域においては、線形変形では追従できないために正しい収束解が得られない場合がある。その場合、本実施例では風上勾配法によって求めたサブピクセル移動量を採用している。以上のようにして、サブピクセル精度の変形ベクトル分布が得られる。
【００６７】
次に、逆解析工程について説明する。
図５は、有限要素モデルを概略的に示す図である。（Ａ）は三角形要素を示し、（Ｂ）は簡易モデルを示す。図６〜図９は、逆解析の処理手順を概略的に示す図である。各図の（Ａ）〜（Ｃ）は、その処理過程を示す。図８は、図７（Ｃ）の部分拡大図に対応する。
【００６８】
本実施例では、ＯＣＴ断層像から求められる変形ベクトル断層分布とＦＥＭとを組み合わせることで温度の断層分布を逆解析推定する。図６（Ａ）には、ＯＣＴにより得られる電子実装部品（測定対象Ｗ）の断層画像が示されている。横軸がｘ座標を示し、縦軸がｚ座標を示す。同図には、高分子基材料からなる基板にＩＣチップが埋設された電子実装部品の一部が示されている。
【００６９】
図６（Ｂ）には、ＯＣＴにより得られた前後の断層画像を用いて上記演算処理を行うことにより得られる変形ベクトル分布が矢印にて示されている。このような変形ベクトル分布に対してＦＥＭの逆解析を施すことにより、測定対象内部の温度分布を算出し、断層表示することができる。
【００７０】
ここで、「逆解析」の意義を明確にするために、ＦＥＭの順解析について理解しておく必要がある。この順解析では、解析対象となる領域（解析領域）を要素分割した後、材料パラメータを設定するとともに荷重等の境界条件を設定し、各要素における温度を与えた後、多元連立方程式となる全体剛性方程式を解く。それにより、各節点における変位（未知項）を算出することができる。なお、本実施例では温度分布と変形分布とが連成するため、順解析は次のようになる。すなわち、解析領域を要素分割した後、熱伝導係数分布を設定するとともに熱（温度）の境界条件を設定し、方程式を解いて温度分布を算出する過程を経る。また、解析領域を要素分割した後、材料パラメータ（剛性、線膨張係数など），力学的境界条件，温度分布を設定し、方程式を解いて変形分布を算出する過程を経る。
【００７１】
これに関し、本実施例では上述のように、ＯＣＴの断層画像に基づいて変形ベクトル分布が算出されるため、解析領域に設定された節点に対してそれらの変形ベクトルを補間することにより、境界条件を設定しなくとも節点上の変位を算出することができる。すなわち、全体剛性方程式における変位ベクトル（変形ベクトル）の項が既知となる。この点に着目し、この既知となった変位ベクトルを用いてＦＥＭの逆解析を実行することにより温度項を算出する。このとき、全節点における変形ベクトルが既知であるため、全要素について温度項を算出することができる。つまり、全要素に対応した温度分布を算出することができる。以下、具体的に説明する。
【００７２】
（ＦＥＭによる温度項を含む全体剛性方程式）
二次元平面応力問題に関し、熱膨張を考慮した応力-ひずみ関係式として下記式（１２）が得られる。なお、ここではＯＣＴの断層画像に合わせてｘ−ｚ平面を考える。
【数１２】
0006707249.tif 000013
ここで、σｘは応力のｘ成分、σｚは応力のｚ成分、τｘｚはせん断応力である。εｘはｘ方向のひずみ、εｚはｚ方向のひずみ、γｘｚはせん断ひずみである。Ｅはヤング率、νはポアソン比、αは線膨張係数、ΔＴは温度変化である。
【００７３】
そして、上記式（１２）に仮想仕事の原理を適用して変形し、要素剛性マトリクス［ｋ］６×６を用いて表すと、下記式（１３）が得られる。
【数１３】
0006707249.tif 000014
【００７４】
本実施例ではＦＥＭを利用するに際し、要素分割される各要素を図５（Ａ）に示す三角形要素にて表し、また、隣接する要素間の関係を図５（Ｂ）に示す簡易モデルにて表す。上記式（１３）において、Ｐ，Ｑは等価節点力であり、ｕ，ｖはそれぞれ節点位置のｘ方向，ｚ方向の変形ベクトル成分である。なお、Ｐ，Ｑ，ｕ，ｖの添え字1〜3は、図５（Ａ）に示す要素分割に係る三角形要素の局所節点番号である。ａ１〜ｂ３は三角形要素において物理量を近似する際の定数であり、それぞれ下記式（１４）により表される。
【数１４】
0006707249.tif 000015
ここで、ｘ，ｚの添え字1〜3は要素分割に係る三角形要素の局所節点番号である。
【００７５】
本実施例では、図７（Ａ）に示すように、測定対象の内部にＲＯＩを設定し、図７（Ｂ）に示すように、そのＲＯＩに対してＦＥＭに基づく要素分割を実行する。ＯＣＴにより得られた変形ベクトルは、図７（Ｃ）に示すようにＲＯＩに分布する。
【００７６】
上記式（１３）を温度の項について展開し、図５（Ｂ）に示す簡易モデルを用いて全体方程式を作成すると、下記式（１５）が成立する。
【数１５】
0006707249.tif 000016
ここで、上付き文字，下付文字はそれぞれ要素番号と節点番号を表している。
【００７７】
上記式（１５）を簡略化および一般化すると下記式（１６）となり、線形化された観測方程式となる。
【数１６】
0006707249.tif 000017
ここで、[Ａ](２ｍin＋２ｍb)×ｎは、温度変化項の観測マトリクスを表し、ＦＥＭの形状関数によって決定される[Ｂ]マトリクスに依存し、それらの成分は節点座標により表される。{ΔＴ}ｎ×１は各要素における温度変化、[Ｋ](２ｍin＋２ｍb)×(２ｍin＋２ｍb)は全体剛性マトリクス、{Ｕ}(２ｍin＋２ｍb)×１は各節点における変形ベクトル成分となる。この変形ベクトル成分はＯＣＴ断層像に基づいて得られるため、右辺は観測ベクトルとなる。
【００７８】
ただし、ｍinはＲＯＩ境界上を除くＲＯＩ内の節点数であり、２ｍinはそれらの節点上に求まる変形ベクトルの成分数である。ｍbはＲＯＩ境界上の節点数であり、２ｍbはそれらの節点上に求まる変形ベクトルの成分数である。したがって、２ｍin＋２ｍbはＲＯＩ以内の全節点に求まる変形ベクトルの成分数となる。ｎはＲＯＩ内の要素数である。上記式（１６）の連立方程式に境界変位条件を導入することで修正した全体剛性方程式に対し、種々の逆解析手法を用いることによって、各要素における温度変化ΔＴを推定することができる。
【００７９】
（境界変位条件の導入）
ＲＯＩ内の有限要素モデルに上記式（１６）を適用した後、各節点における変形ベクトル成分{Ｕ}(２ｍin＋２ｍb)×１を節点補間により与える。このとき、ＲＯＩ境界上の節点に得られる変形ベクトル成分ｕ，ｖを境界強制変位（拘束節点）とする必要がある。そのため、上記式（１６）における全体剛性マトリクス[Ｋ](２ｍin＋２ｍb)×(２ｍin＋２ｍb)、および観測マトリクス[Ａ](２ｍin＋２ｍb)×ｎを修正してマトリクスの縮小を行う。
【００８０】
すなわち、マトリクス成分の添え字の意味を改変した上で上記式（１６）を改めて書き下すと下記式（１７）となる。
【数１７】
0006707249.tif 000018
【００８１】
この全体剛性方程式を解く過程で境界変位条件を満足するように修正する。ここで、ＲＯＩ境界上の節点に得られた変形ベクトル成分をｕｌ，ｖｌとする。このとき、{Ｕ}(２ｍin＋２ｍb)×１において境界上の節点に求まる変形ベクトル成分を表す行はｕｌが２ｌ−１行目、ｖｌが２ｌ行目となる。ＲＯＩ境界の節点では境界変位条件（変位拘束境界）として、ＲＯＩ境界上節点の既知の変位量ｕｌ，ｖｌについて左右両辺が満たされるように、全体剛性マトリクス[Ｋ](２ｍin＋２ｍb)×(２ｍin＋２ｍb)、および観測マトリクス[Ａ](２ｍin＋２ｍb)×ｎを修正する。具体的には、ＲＯＩ境界上節点における変位量ｕｌ，ｖｌが，左辺と右辺においてｕｌ＝ｕｌ，ｖｌ＝ｖｌとなるように、[Ｋ]の成分と、[Ａ]における同じ行の成分を下記式（１８）に従って修正する。
【数１８】
0006707249.tif 000019
【００８２】
また、ベクトル{Ｆ}の境界上の節点に対応する２ｌ−１，２ｌ行目も、それぞれｕｌ，ｖｌと修正して左右両辺が満たされるようにする。さらに、ＲＯＩが物体内部に設定され、かつ静的変形により体積力が作用しないのであれば、{Ｆ}に関してＲＯＩ境界上を除くＲＯＩ内の等価節点力はゼロとなる。また、ＲＯＩ境界を自由表面に設定した場合は自由変形となるため、表面力は{Ｆ}＝０となる。このため、上記式（１７）は下記式（１９）のようになる。
【数１９】
0006707249.tif 000020
【００８３】
以上のようにして、全体剛性方程式に境界変位条件を満足させる。なお、上記式（１９）において、両辺の２ｌ−１，２ｌ行目は拘束節点に帰属し、境界変位条件を満たすため、それらを省略することで成分を並び替え、マトリクスの縮小を行う。このようにして得られる観測マトリクス[Ａ]を[Ａ']と表し、全体剛性マトリクス[Ｋ]と変形ベクトル{Ｕ}の積を{Ｋ'Ｕ'}と表すと、下記式（２０）が得られる。
【数２０】
0006707249.tif 000021
【００８４】
（節点補間）
この逆解析のためには、上記式（１６）における測定項｛Ｕ｝を与えなければならない。一方、図８（Ａ）に示すように、ＯＣＴにより得られた変形ベクトル分布は、ＲＯＩ内に分割された要素の節点に対応していないため、これをそのまま測定項｛Ｕ｝とすることはできない。そこで、図８（Ｂ）に示すようにＲＯＩの節点上に変形ベクトルを設定するための補間処理を実行する。
【００８５】
すなわち、ＯＣＴにより取得した変形ベクトルの断層分布を用いて各節点上に変形ベクトルを補間する。ここでは、各節点を中心とする所定半径（本実施例では３００μｍ）を「影響半径」とし、その影響半径に収まる変形ベクトルをその節点に補間した。この補間に際しては、最小二乗法により下記式（２１）に示す近似関数ｆ(x,z)を求め、その近似関数ｆ(x,z)に則って各節点位置の変形ベクトルを算出した。
【数２１】
0006707249.tif 000022
ここで、ｆ(x,z)は二次元二次多項式とした近似関数である。ａ〜ｇはその未知係数であり、ｘ，ｚはＯＣＴにより得られる変形ベクトルの始点座標値である。
【００８６】
この近似関数ｆ(x,z)について下記式（２２）を適用する。
【数２２】
0006707249.tif 000023
ここで、Ｑは、近似関数ｆ(x,z)と測定値ｙｉとの残差の二乗和である。測定値ｙｉはＯＣＴにより得られる変形量成分であり、Ｎは影響半径内にある変形ベクトルのデータ数である。上記式（２２）においてＱが最小となるときの近似関数ｆ(x,z)を求め、そのｆ(x,z)に則ってＲＯＩにある節点上の変形量成分を算出する。
【００８７】
具体的には、Ｑを算出した後、Ｑを近似関数の各未知係数について偏微分することにより、下記式（２３）に示す正規式を得る。
【数２３】
0006707249.tif 000024
【００８８】
この正規式を解き、未知係数を算出することで近似関数を求める。求まった近似関数に節点座標を参照させ、節点位置での変形量成分を求めることで節点補間がなされる。ただし、表面近傍のＲＯＩ境界や深部のＲＯＩ境界では輝度情報にノイズが多く含まれるため、変形ベクトルの算出精度が低い。さらに、ＲＯＩ境界近傍、特にＲＯＩの頂点部分（角部分）では使用できる変形ベクトルの数が少ないため、ＲＯＩ境界外に変形ベクトルを外挿補間している。このような手法により、ＲＯＩ内の節点上の各変形量成分を求めて節点補間をする。なお、上記式（１６）の{Ｕ}は、ＯＣＴによる測定データとして得られる。
【００８９】
（逆解析）
上記式（２０）における観測マトリクス[Ａ']は一般に２ｍin×ｎの行列である。ｍinはＲＯＩ境界を除いたＲＯＩ内部の節点ｘ，ｚ座標の数であり、ｎはＲＯＩ内の要素数である。ｍin＜ｎとなるため、[Ａ']は正則ではない。そのため、本実施例では、逆解析手法として擬似逆行列であるムーアペンローズ一般逆行列を求め、逆解析を行った。ここで、観測マトリクス[Ａ']がフルランクでランクがｎの場合、[Ａ']の一般逆行列は下記式（２４）にて表される。
【数２４】
0006707249.tif 000025
【００９０】
上記式（２４）に基づき、温度の逆解析を下記式（２５）を用いて実行する。
【数２５】
0006707249.tif 000026
本実施例では、この[Ａ']のムーアペンローズ一般逆行列を上記式（２０）の両辺に乗算することにより、温度断層分布｛ΔＴ｝を推定値として算出した。
【００９１】
ただし、ムーアペンローズ一般逆行列により算出された推定値は規格化されており、物理量として単位をもたない値となる。そこで、この推定値を有用な物理量の温度値に較正する。すなわち、図９（Ａ）に点線枠にて示すように、推定値の一部に対応する温度の真値を取得し、その真値に対応する推定値（「対応推定値」ともいう）とその真値とに基づいて、全推定値を線形補間することにより較正を行う。この真値は実際に計測した値（実計測値）であってよい。
【００９２】
すなわち、逆解析にて得た温度分布を、下記式（２６）に基づいて較正する。
【数２６】
0006707249.tif 000027
ここで、Ｔ*は逆解析により算出される推定値であり、Ｔ*realizedは較正後の温度値（実値）である。図９（Ｂ）に概念的に示すように、Ｔsample1，Ｔsample2は、較正に用いる任意の２点における真値であり、Ｔ*sample1，Ｔ*sample2はその真値に対応する対応推定値である。なお、本実施例では線形（一次関数）による較正を例示しているが、多項式にて較正（多項式近似）してもよい。
【００９３】
制御演算部６は、以上のような演算処理を実行し、算出された温度値Ｔ*realizedを表示装置８に断層可視化する。それにより、例えば図９（Ｃ）に示すような色分けされた温度分布が表示可能となる。
【００９４】
次に、制御演算部６が実行する具体的処理の流れについて説明する。
図１０は、制御演算部６により実行されるマイクロ断層可視化処理の流れを示すフローチャートである。なお、この断層可視化処理は、上述したＯＣＴによる断層撮影がなされることを前提に実行される。ここでは、測定対象Ｗが熱変形したときの変形前後の断層画像を取得し、物理量として測定対象Ｗの内部の温度分布を可視化する場合を例に説明する。
【００９５】
制御演算部６は、ＯＣＴにより撮影された前後２枚のＯＣＴ断層画像を読み込む（Ｓ１０）。続いて、それらの断層画像に基づいて変形ベクトル分布を算出するための変形ベクトル演算処理を実行する（Ｓ１２）。続いて、その変形ベクトル分布を用いて逆解析処理を実行する（Ｓ１４）。そして、その逆解析により算出された物理量の断層分布（温度分布）を表示装置８に表示させる（Ｓ１６）。
【００９６】
図１１は、図１０におけるＳ１２の変形ベクトル演算処理を詳細に示すフローチャートである。制御演算部６は、変形前後の断層画像に基づき、再帰的相互相関法による処理を実行する。ここではまず、最小解像度（最大サイズの検査領域）での相互相関処理を実行し、相関係数分布を求める（Ｓ２０）。続いて、隣接相互相関乗法により、隣接する相関係数分布の積を演算する（Ｓ２２）。このとき、標準偏差フィルタ等の空間フィルタにより過誤ベクトルの除去をし（Ｓ２４）、最小二乗法等により除去ベクトルの内挿補間を実行する（Ｓ２６）。続いて、検査領域を小さくすることによって解像度を上げて相互相関処理を継続する（Ｓ２８）。すなわち、低解像度での参照ベクトルを基に相互相関処理を実行する。このときの解像度が予め定める最高解像度でなければ（Ｓ３０のＮ）、Ｓ２２に戻る。
【００９７】
そして、Ｓ２２〜Ｓ２８の処理を繰り返し、最高解像度での相互相関処理が完了すると（Ｓ３０のＹ）、サブピクセル解析を実行する。すなわち、最高解像度（最小サイズの検査領域）での変形ベクトルの分布に基づき、風上勾配法によるサブピクセル移動量を演算する（Ｓ３２）。そして、このとき算出されたサブピクセル移動量に基づき、画像変形法によるサブピクセル変形量を演算する（Ｓ３４）。続いて、最大相互相関値によるフィルタ処理により過誤ベクトルの除去をし（Ｓ３６）、最小二乗法等により除去ベクトルの内挿補間を実行する（Ｓ３８）。このようにして変形ベクトル分布が得られる。
【００９８】
図１２は、図１０におけるＳ１４の逆解析処理を詳細に示すフローチャートである。制御演算部６は、演算対象となった断層画像に対応するＲＯＩを設定する（Ｓ４０：図７（Ａ）参照）。この設定は、当該解析を行うユーザの入力指示にしたがって行われる。ユーザは、測定対象Ｗにおいて温度分布を表示させたい箇所をＲＯＩとして設定することができ、入力装置７を介してそれを指示することができる。
【００９９】
続いて、制御演算部６は、設定されたＲＯＩに対して要素分割を設定する（Ｓ４２：図７（Ｂ）参照）。この要素分割は、上述したＦＥＭに対応するものであり、ユーザの入力指示にしたがって行われる。ユーザは、求める解析精度に応じた要素分割を入力装置７を介して指示することができる。
【０１００】
続いて、制御演算部６は、Ｓ３８を経て得られた変形ベクトル分布を、分割された要素の節点に補間する（Ｓ４４：図８参照）。また、ＲＯＩに対して弾性率および線膨張係数を設定する（Ｓ４６）。本実施例では、これらの値を測定対象Ｗの材料に固有の値として設定する。
【０１０１】
続いて、制御演算部６は、仮想仕事の原理に基づいて上記式（１６）の全体剛性マトリクス［Ｋ］を演算する一方（Ｓ４８）、形状関数，力学パラメータ（弾性率など），温度依存パラメータ（線膨張係数）から温度変化項の全体マトリクス観測行列［Ａ］を決定する（Ｓ５０）。そして、上述したマトリクスの縮小などを行い、逆解析（例えば逆行列演算）により温度分布を演算する（Ｓ５２）。以上のようにして算出された温度分布が、図１０のＳ１６にて断層可視化される。
【０１０２】
次に、本実施例による効果について説明する。発明者は、本実施例による手法、つまりＯＣＴによる測定とＦＥＭの逆解析とによるハイブリッド手法の有効性を確認するために、数値実験による精度検証を行った。以下、その検証結果について説明する。図１３〜図１５は、数値実験による検証方法を示す図である。各図の（Ａ）〜（Ｃ）はその検証過程を示す。図１６は、検証結果を示す図である。（Ａ）は本実施例のハイブリッド手法による解析結果を示し、（Ｂ）はＦＥＭによる正解値を示す。図１７は、本実施例の解析による正解値との誤差を示す図である。（Ａ）は絶対誤差を示し、（Ｂ）は相対誤差を示す。
【０１０３】
上記ハイブリッド手法の精度を評価するためには、その比較対象となる正解値が必要となる。本来であれば、この比較対象を温度センサ等による実測値（真値）として得るのが望ましいが、物体内部の断層に沿ってそれを行うのは不可能に等しい。そこで、ここでは測定対象のモデルに対してＦＥＭによる順解析を実行し、それにより得られる温度分布を正解値とする。そして、その正解値からＯＣＴ断層画像を疑似的に作成し、その断層画像（以下「疑似ＯＣＴ断層画像」ともいう）をＯＣＴによる実際の断層画像と仮定して上記逆解析を実行する。その逆解析によって得られた温度分布と、正解値としての温度分布（以下「正解温度分布」ともいう）とを比較することで、上記ハイブリッド手法の精度を評価した。
【０１０４】
（数値実験）
具体的には、モンテカルロシミュレーションによって散乱体からのランダムスッペクルパターンを演算することで、図１４（Ｂ）に示す疑似ＯＣＴ断層画像を作成し、ＦＥＭを用いて散乱体の熱変形移動量の算出を行うこととした。
【０１０５】
まず、数値モデルに基づいたＦＥＭによる伝熱解析を行った。図１３（Ａ）に示すように、測定対象として基板にＩＣチップ（シリコンチップ）が埋め込まれた電子実装部品を模擬したモデルを想定している。このモデルは、エポキシ樹脂（Epoxy）を基材とし、シリコンチップのシリコン（Silicon）、配線の銅（Copper）を内包し、表面にはレジスト材（Resist）が存在する複合材料とされている。
【０１０６】
そして、図１３（Ｂ）に示すように、このモデルに対してＦＥＭの順解析による要素分割を施した。なお、モデルを構成する各材料の材料定数は、下記の表１に示すとおりである。
【表１】
0006707249.tif 000028
【０１０７】
温度の境界条件については、シリコンチップおよび配線が発熱していると仮定した。すなわち、シリコンおよび銅の温度を２５０[℃]とし、表面（上下面）であるｚ＝０，２０００[μｍ]の温度を外気温に等しい２５[℃]とした。また、熱伝達係数を１０[Ｗ/(ｍＫ)]とした。ただし、モデルの側面（ｘ＝０，５０００[μｍ]）では熱の出入りはないと仮定した。このような設定のもとＦＥＭの伝熱解析を実行することで、図１３（Ｃ）に示す温度分布を得た。この温度分布が数値実験における「正解温度分布」とされる。
【０１０８】
この温度分布データを用いてＦＥＭによる力学解析を実行し、モデルを熱変形させることとした。その際、モデルの境界拘束条件として、ｚ＝０[μｍ]の位置を全方向に固定することとした。ＩＣチップを内包する電子実装部品は、マザー基板へ実装される際にその底面がはんだ付けされることを考慮したものである。
【０１０９】
図１４（Ａ）には、このような連成解析により算出された熱変形ベクトル分布が示されている。この熱変形量データを図１４（Ｂ）の各スペックルに対応させ、スペックルの移動を決定した。各スペックルの熱変形移動は、スペックルを内包する有限要素の節点変形ベクトルおよび形状関数を用いて算出することができる。それにより、図１４（Ｃ）に示す熱変形後の疑似ＯＣＴ断層画像を作成した。
【０１１０】
このようにして、図１４（Ｂ）および（Ｃ）に示すように、熱変形前後の疑似ＯＣＴ断層画像を得た。なお、ここでいう「熱変形前」とは、シリコンチップおよび配線が発熱する前の状態、つまり両者の温度が実質的に外気温である２５[℃]に等しい状態を意味する。「熱変形後」とは、シリコンチップおよび配線が上述のように２５０[℃]に発熱した状態を意味する。熱変形前後の疑似ＯＣＴ断層画像は、それぞれの状態を仮定して得たものである。この２つの疑似ＯＣＴ断層画像の取得を出発点として上記ハイブリッド手法を適用した。すなわち、変形前後の疑似ＯＣＴ断層画像をＯＣＴによる実際の断層画像と仮定し、図１０〜図１２に示した処理を実行した。そして、図１０のＳ１６にて得られた温度分布を、図１３（Ｃ）に示した正解温度分布と比較することで精度を評価した。
【０１１１】
具体的には、図１５（Ａ）に示す疑似ＯＣＴ断層画像に対し、図１５（Ｂ）に示すようにＲＯＩを設定し、そのＲＯＩについて要素分割を行った。なお、このＲＯＩ内の弾性率分布，線膨張係数分布を定義する必要があるが、弾性率および線膨張係数は材料に固有の値であるため、ここでは上記表１に示す材料定数を設定した。
【０１１２】
図１５（Ｃ）には、前後の疑似ＯＣＴ断層画像に基づいて算出された変形ベクトル分布が示されている。この変形ベクトル分布に基づいて上述した逆解析処理を実行し、図１６（Ａ）に示すように温度分布を断層可視化した。この温度分布を図１６（Ｂ）に示す正解温度分布と比較すると、それらが定性的に一致していることが分かる。すなわち、いずれも基材中心部分の温度がその周囲より低くなっているのに対し、シリコンチップおよび銅の付近で周囲より高い温度分布が得られていることが確認できる。
【０１１３】
図１７（Ａ）には、本実施例のハイブリッド手法により得られた温度分布（「推定温度分布」ともいう）と、正解温度分布との絶対誤差が示されている。図１７（Ｂ）には両者の相対誤差が示されている。なお、「絶対誤差」とは両者の温度差を意味し、「相対誤差」とは絶対誤差を正解値により除算した値を意味する。
【０１１４】
これらの結果に基づき、両者の誤差平均（ＲＭＳ誤差）は１４％程度と比較的小さく収まることが確認できた。なお、ＲＯＩの境界近傍でエラーが大きくなる傾向にあるが、これは境界でのノイズの影響によるものと考えられる。なお、このような温度推定の誤差は、温度昇降を動的に検出してノイズ低減を施す方法や、ノイズを考慮した逆推定手法であるカルマンフィルタなどの導入により改善可能である。
【０１１５】
以上に説明したように、本実施例によれば、ＯＣＴによる断層計測（実計測データ）を用いることにより、材料内部の温度断層分布を逆解析することが可能となる。すなわち、実測と解析とのハイブリッドにより物理量を断層可視化できる。物体内部における任意の領域の温度分布を、境界条件を考慮することなく断層可視化できることから、特に境界条件が不明確な状況における内部温度の検出に好適となる。また、このような温度断層分布がマイクロスケールで非侵襲・非破壊にて検出可能であるため、ＭＥＭＳ等の電子実装部品のように境界条件の付与が難しい工業製品への応用が期待される。
【０１１６】
以上、本発明の好適な実施例について説明したが、本発明はその特定の実施例に限定されるものではなく、本発明の技術思想の範囲内で種々の変形が可能であることはいうまでもない。
【０１１７】
上記実施例では、ＯＣＴ断層画像として、構造の変形前後の静止画を用いる例を示した。変形例においては、構造の変形中を撮影したＯＣＴ断層画像の動画に対して上記手法を適用してもよい。すなわち、動画中の異なる時刻の２枚のＯＣＴ断層画像を取得し、それらに対して上記逆解析を適用して物理量を演算し、断層可視化してもよい。上記ハイブリッド手法は、このように時間変化する現象に対しても適用可能である。
【０１１８】
上記実施例では述べなかったが、測定対象の温度場が変化しているような非定常な状況では、変形している場（分布）が時々刻々と変化する。そのような状況下においてＯＣＴ断層画像を連続撮影して本手法を適用することにより、時々刻々変化する非定常温度場（分布）をマイクロ断層可視化することが可能となる。この演算過程において、連続するＯＣＴ断層画像と時空間最小二乗法を用いた節点補間を行うなどが考えられる。
【０１１９】
上記実施例では、変形ベクトルの断層分布に基づいて逆解析を行い、それにより温度分布を算出する例を示した。変形例においては、変形ベクトルの時間微分である変形速度ベクトルの断層分布に基づいて逆解析を行い、それにより物理量の分布を算出してもよい。
【０１２０】
上記実施例では、ＯＣＴによる断層計測と数値解析の逆解析とにより部品内部の温度断層分布を可視化する例を示した。変形例においては、同様の手法により部品内部の弾性係数や粘弾性係数の分布を可視表示することも可能である。例えば、上記式（２０）に関して温度変化｛ΔＴ｝がゼロであるなど既知の場合、測定項｛Ｕ｝の取得により全体剛性マトリクス［Ｋ］を算出し、弾性係数の分布を可視表示することができる。それにより、部品内部に存在する材料を特定することも可能となる。
【０１２１】
上記実施例では、ＯＣＴによる断層画像から得た変形ベクトル分布を逆解析することにより物理量を算出する例を示したが、ＯＣＴに代えてＸ線ＣＴ、超音波、ＭＲＩなどを採用することもできる。例えば、電子実装部品などの高分子基複合材や金属との複合材に適用する際にはマイクロＸ線ＣＴによる断層画像などが十分適用できると考えられる。
【０１２２】
上記実施例では、ＯＣＴによる断層画像を二次元で取得する例を示したが、三次元で取得してもよい。すなわち、変位関連ベクトルとしての変形ベクトルあるいは変形速度ベクトルを三次元データとして処理してもよいことは言うまでもない。
【０１２３】
上記実施例では、物理量を断層可視化する対象として部品構造を例示した。変形例においては、その他の材料構造を対象としてもよい。また、体内の軟骨やしこり（硬化箇所）など、生体構造（生体組織）を対象としてもよい。生体構造の力学特性の変化を断層可視化することで、例えば動脈硬化，癌，肝硬変，皮膚，軟骨，骨などの多くの症例に対する診断材料とすることが可能となる。
【０１２４】
なお、本発明は上記実施例や変形例に限定されるものではなく、要旨を逸脱しない範囲で構成要素を変形して具体化することができる。上記実施例や変形例に開示されている複数の構成要素を適宜組み合わせることにより種々の発明を形成してもよい。また、上記実施例や変形例に示される全構成要素からいくつかの構成要素を削除してもよい。
【符号の説明】
【０１２５】
１  断層可視化システム、２  光学ユニット、４  光学機構、６  制御演算部、７  入力装置、８  表示装置、１０  光源、１２  オブジェクトアーム、１４  リファレンスアーム、１６  光検出器、２６  反射鏡、ＲＯＩ  関心領域、Ｓ１  検査領域、Ｓ２  探査領域。
【図１】
0006707249.tif 000029
【図２】
0006707249.tif 000030
【図３】
0006707249.tif 000031
【図４】
0006707249.tif 000032
【図５】
0006707249.tif 000033
【図６】
0006707249.tif 000034
【図７】
0006707249.tif 000035
【図８】
0006707249.tif 000036
【図９】
0006707249.tif 000037
【図１０】
0006707249.tif 000038
【図１１】
0006707249.tif 000039
【図１２】
0006707249.tif 000040
【図１３】
0006707249.tif 000041
【図１４】
0006707249.tif 000042
【図１５】
0006707249.tif 000043
【図１６】
0006707249.tif 000044
【図１７】
0006707249.tif 000045