JP6707249B2 - Micro-tomography visualization method and system - Google Patents
Micro-tomography visualization method and system Download PDFInfo
- Publication number
- JP6707249B2 JP6707249B2 JP2015050553A JP2015050553A JP6707249B2 JP 6707249 B2 JP6707249 B2 JP 6707249B2 JP 2015050553 A JP2015050553 A JP 2015050553A JP 2015050553 A JP2015050553 A JP 2015050553A JP 6707249 B2 JP6707249 B2 JP 6707249B2
- Authority
- JP
- Japan
- Prior art keywords
- distribution
- analysis
- physical quantity
- deformation
- tomographic
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Active
Links
- 238000010603 microCT Methods 0.000 title claims description 15
- 238000007794 visualization technique Methods 0.000 title claims description 6
- 238000000034 method Methods 0.000 claims description 154
- 238000009826 distribution Methods 0.000 claims description 147
- 238000004458 analytical method Methods 0.000 claims description 135
- 239000013598 vector Substances 0.000 claims description 135
- 238000012014 optical coherence tomography Methods 0.000 claims description 68
- 238000004364 calculation method Methods 0.000 claims description 52
- 238000006073 displacement reaction Methods 0.000 claims description 51
- 230000008569 process Effects 0.000 claims description 32
- 239000000463 material Substances 0.000 claims description 26
- 230000008859 change Effects 0.000 claims description 18
- 238000012800 visualization Methods 0.000 claims description 13
- 230000002441 reversible effect Effects 0.000 claims description 8
- 229920000642 polymer Polymers 0.000 claims description 5
- 230000008685 targeting Effects 0.000 claims description 3
- 238000003325 tomography Methods 0.000 claims description 3
- 230000035882 stress Effects 0.000 claims 2
- 230000008646 thermal stress Effects 0.000 claims 1
- 238000005259 measurement Methods 0.000 description 46
- 230000003287 optical effect Effects 0.000 description 40
- 238000007689 inspection Methods 0.000 description 39
- 239000011159 matrix material Substances 0.000 description 29
- 238000012545 processing Methods 0.000 description 26
- 230000007246 mechanism Effects 0.000 description 14
- 238000012795 verification Methods 0.000 description 10
- 238000010586 diagram Methods 0.000 description 9
- 238000002474 experimental method Methods 0.000 description 8
- XUIMIQQOPSSXEZ-UHFFFAOYSA-N Silicon Chemical compound [Si] XUIMIQQOPSSXEZ-UHFFFAOYSA-N 0.000 description 7
- 239000010703 silicon Substances 0.000 description 7
- 229910052710 silicon Inorganic materials 0.000 description 7
- 239000002131 composite material Substances 0.000 description 6
- 238000012986 modification Methods 0.000 description 6
- 230000004048 modification Effects 0.000 description 6
- 230000008602 contraction Effects 0.000 description 5
- 238000013461 design Methods 0.000 description 5
- 239000000758 substrate Substances 0.000 description 5
- 238000011426 transformation method Methods 0.000 description 5
- 239000000470 constituent Substances 0.000 description 4
- RYGMFSIKBFXOCR-UHFFFAOYSA-N Copper Chemical compound [Cu] RYGMFSIKBFXOCR-UHFFFAOYSA-N 0.000 description 3
- 229910052802 copper Inorganic materials 0.000 description 3
- 239000010949 copper Substances 0.000 description 3
- 238000001514 detection method Methods 0.000 description 3
- 230000009467 reduction Effects 0.000 description 3
- 210000001519 tissue Anatomy 0.000 description 3
- 238000012546 transfer Methods 0.000 description 3
- 230000009466 transformation Effects 0.000 description 3
- 238000005311 autocorrelation function Methods 0.000 description 2
- 210000000845 cartilage Anatomy 0.000 description 2
- 230000001419 dependent effect Effects 0.000 description 2
- 230000000694 effects Effects 0.000 description 2
- 238000004100 electronic packaging Methods 0.000 description 2
- 238000013178 mathematical model Methods 0.000 description 2
- 238000002715 modification method Methods 0.000 description 2
- 239000013307 optical fiber Substances 0.000 description 2
- NAWXUBYGYWOOIX-SFHVURJKSA-N (2s)-2-[[4-[2-(2,4-diaminoquinazolin-6-yl)ethyl]benzoyl]amino]-4-methylidenepentanedioic acid Chemical compound C1=CC2=NC(N)=NC(N)=C2C=C1CCC1=CC=C(C(=O)N[C@@H](CC(=C)C(O)=O)C(O)=O)C=C1 NAWXUBYGYWOOIX-SFHVURJKSA-N 0.000 description 1
- 206010003210 Arteriosclerosis Diseases 0.000 description 1
- 239000004593 Epoxy Substances 0.000 description 1
- 229920002430 Fibre-reinforced plastic Polymers 0.000 description 1
- 102100040287 GTP cyclohydrolase 1 feedback regulatory protein Human genes 0.000 description 1
- 101710185324 GTP cyclohydrolase 1 feedback regulatory protein Proteins 0.000 description 1
- HBBGRARXTFLTSG-UHFFFAOYSA-N Lithium ion Chemical compound [Li+] HBBGRARXTFLTSG-UHFFFAOYSA-N 0.000 description 1
- 238000000342 Monte Carlo simulation Methods 0.000 description 1
- 206010028980 Neoplasm Diseases 0.000 description 1
- 238000002940 Newton-Raphson method Methods 0.000 description 1
- 210000003484 anatomy Anatomy 0.000 description 1
- 208000011775 arteriosclerosis disease Diseases 0.000 description 1
- 210000000988 bone and bone Anatomy 0.000 description 1
- 238000004422 calculation algorithm Methods 0.000 description 1
- 201000011510 cancer Diseases 0.000 description 1
- 239000003990 capacitor Substances 0.000 description 1
- 239000004918 carbon fiber reinforced polymer Substances 0.000 description 1
- 208000019425 cirrhosis of liver Diseases 0.000 description 1
- 239000011231 conductive filler Substances 0.000 description 1
- 239000002178 crystalline material Substances 0.000 description 1
- 238000012888 cubic function Methods 0.000 description 1
- 230000006378 damage Effects 0.000 description 1
- 230000007423 decrease Effects 0.000 description 1
- 230000001066 destructive effect Effects 0.000 description 1
- ALKZAGKDWUSJED-UHFFFAOYSA-N dinuclear copper ion Chemical compound [Cu].[Cu] ALKZAGKDWUSJED-UHFFFAOYSA-N 0.000 description 1
- 239000003822 epoxy resin Substances 0.000 description 1
- 239000011151 fibre-reinforced plastic Substances 0.000 description 1
- 238000001914 filtration Methods 0.000 description 1
- 230000020169 heat generation Effects 0.000 description 1
- 238000003384 imaging method Methods 0.000 description 1
- 238000012886 linear function Methods 0.000 description 1
- 239000004973 liquid crystal related substance Substances 0.000 description 1
- 229910001416 lithium ion Inorganic materials 0.000 description 1
- 239000002184 metal Substances 0.000 description 1
- 229910052751 metal Inorganic materials 0.000 description 1
- 239000000203 mixture Substances 0.000 description 1
- 239000002105 nanoparticle Substances 0.000 description 1
- 230000035699 permeability Effects 0.000 description 1
- 230000000704 physical effect Effects 0.000 description 1
- 239000004033 plastic Substances 0.000 description 1
- 229920003023 plastic Polymers 0.000 description 1
- 229920000647 polyepoxide Polymers 0.000 description 1
- 239000002952 polymeric resin Substances 0.000 description 1
- 238000003672 processing method Methods 0.000 description 1
- 230000003252 repetitive effect Effects 0.000 description 1
- 229920005989 resin Polymers 0.000 description 1
- 239000011347 resin Substances 0.000 description 1
- SBEQWOXEGHQIMW-UHFFFAOYSA-N silicon Chemical compound [Si].[Si] SBEQWOXEGHQIMW-UHFFFAOYSA-N 0.000 description 1
- 210000003491 skin Anatomy 0.000 description 1
- 230000003595 spectral effect Effects 0.000 description 1
- 230000003068 static effect Effects 0.000 description 1
- 229920003002 synthetic resin Polymers 0.000 description 1
- 238000002076 thermal analysis method Methods 0.000 description 1
- 238000001931 thermography Methods 0.000 description 1
- 239000012780 transparent material Substances 0.000 description 1
Landscapes
- Analysing Materials By The Use Of Radiation (AREA)
- Investigating Or Analysing Materials By Optical Means (AREA)
Description
本発明は、構造内部における物理量のマイクロ断層可視化を可能にする方法およびシステムに関する。 The present invention relates to methods and systems that enable microtoming of physical quantities within structures.
パーソナルコンピュータや携帯端末などの電子機器には、いわゆるMEMS(Micro Electro Mechanical Systems)と呼ばれる電子実装部品が組み込まれている。このような電子実装部品の小型化,高密度化に伴うジュール発熱密度の増大により、熱ひずみによる材料内部の微視的破壊が問題となっている。特に、導電性フィラーを含有する繊維強化プラスチックなど(CFRP,GFRP等)、温度依存導電性を有する高分子基複合材料からなる部品について、その問題が顕著となっている。また、例えばリチウムイオン二次電池用キャパシタは、薄板電極間における温度上昇がドライアップ故障を引き起こすといった問題を有する。このため、これらの工業製品にはマイクロスケールでの最適熱設計が必要とされる。 Electronic equipment such as a so-called MEMS (Micro Electro Mechanical Systems) is incorporated in electronic devices such as personal computers and mobile terminals. Due to the increase in Joule's heat generation density due to the miniaturization and high density of electronic packaging components, microscopic destruction inside the material due to thermal strain has become a problem. In particular, the problem is remarkable in parts made of a polymer-based composite material having temperature-dependent conductivity, such as fiber reinforced plastics containing a conductive filler (CFRP, GFRP, etc.). Further, for example, a lithium ion secondary battery capacitor has a problem that a temperature rise between thin plate electrodes causes a dry-up failure. Therefore, these industrial products require optimal thermal design at the microscale.
このような工業製品の熱設計においては、構造内部の温度評価が重要となるが、内部の温度情報が所望の位置でマイクロスケールにて得られる実計測システムは未だ存在せず、有限要素法等の数値解析に頼るところが大きい(例えば特許文献1参照)。また、設計という観点からは、このような熱設計だけでなく、強度(剛性)など力学特性を考慮した設計も重要となるが、これも同様の数値解析に頼るところが大きい。さらに、工業製品を離れると、例えば生体組織ではその力学物理量を測定することが困難であるため、これも同様の数値解析に依存する。 In the thermal design of such industrial products, it is important to evaluate the temperature inside the structure, but there is no actual measurement system that can obtain the internal temperature information at the desired position on the microscale, and the finite element method, etc. There is a large reliance on the numerical analysis of (see, for example, Patent Document 1). Further, from the viewpoint of design, not only such thermal design but also design considering mechanical properties such as strength (rigidity) is important, but this also relies heavily on similar numerical analysis. Furthermore, when the industrial product is left, it is difficult to measure the mechanical and physical quantity of, for example, a living tissue, and this also depends on the same numerical analysis.
しかしながら、このような熱解析や力学解析などの数値解析をマイクロスケールにて行う場合、それらにおける境界条件の与え方が難しい。また、局所的なデータの取得を目的とする場合にも部品全体を解析対象とする必要があり、演算効率上の問題も生じる。 However, when performing numerical analysis such as thermal analysis or mechanical analysis on a microscale, it is difficult to give boundary conditions in them. Further, even when the purpose is to acquire local data, it is necessary to analyze the entire component, which causes a problem in calculation efficiency.
本発明はこのような課題に鑑みてなされたものであり、その目的の一つは、構造内部における物理量を簡易な手法にてマイクロ断層可視化可能な技術を提供することにある。 The present invention has been made in view of such problems, and one of the objects thereof is to provide a technique capable of visualizing a micro-tomography by a simple method for a physical quantity inside a structure.
本発明のある態様は、構造を対象とする順解析の過程で変位関連ベクトルが演算される所定の数値解析において解析条件又は解析結果となり得る物理量に関し、構造の内部における物理量の分布を断層可視化するマイクロ断層可視化方法である。この方法は、構造の変形前後の断層画像を取得する工程と、断層画像に基づき変位関連ベクトルの断層分布を演算する工程と、演算された変位関連ベクトルの断層分布を用いて数値解析の逆解析を実行することにより、物理量の断層分布を演算する工程と、演算された物理量の分布を表示装置に断層可視化する工程と、を備える。なお、ここでいう「変位関連ベクトル」は、変形ベクトルであってもよいし、変形速度ベクトルであってもよい。 An aspect of the present invention relates to a physical quantity that can be an analysis condition or an analysis result in a predetermined numerical analysis in which a displacement-related vector is calculated in the process of a forward analysis targeting a structure, and visualizes the distribution of the physical quantity inside the structure as a tomographic image. It is a micro-tomography visualization method. This method includes the steps of obtaining tomographic images before and after structural deformation, calculating a tomographic distribution of displacement-related vectors based on the tomographic images, and performing an inverse analysis of numerical analysis using the calculated tomographic distribution of displacement-related vectors. By executing the step of calculating the tomographic distribution of the physical quantity, and the step of visualizing the calculated distribution of the physical quantity on the display device. Note that the “displacement-related vector” here may be a deformation vector or a deformation velocity vector.
この態様によると、実測により得られる断層画像から変位関連ベクトルの分布が演算され、その変位関連ベクトル分布に基づき所定の数値解析の逆解析を行うことで物理量の分布が得られる。すなわち、実測と数値解析の逆解析とのハイブリッドにより物理量を断層可視化することが可能となる。実測により内部の変位関連ベクトルが得られるため、数値解析の順解析のように境界部に境界条件を設定する必要もない。このため、所望の断層位置について物理量を簡易に算出することができる。すなわち、この態様によれば、構造内部における物理量を簡易な手法にてマイクロ断層可視化することができる。 According to this aspect, the distribution of the displacement-related vector is calculated from the tomographic image obtained by the actual measurement, and the physical quantity distribution is obtained by performing the inverse analysis of the predetermined numerical analysis based on the displacement-related vector distribution. That is, the physical quantity can be visualized by a hybrid of the actual measurement and the inverse analysis of the numerical analysis. Since the internal displacement-related vector is obtained by the actual measurement, it is not necessary to set the boundary condition at the boundary portion unlike the numerical analysis forward analysis. Therefore, it is possible to easily calculate the physical quantity for the desired tomographic position. That is, according to this aspect, it is possible to visualize the physical quantity inside the structure by a simple method.
本発明の別の態様は、測定対象となる構造の内部における所定の物理量の分布を断層可視化するマイクロ断層可視化システムである。このシステムは、構造の変形前後の断層画像を取得するための断層画像取得装置と、断層画像取得装置から構造の変形前後の断層画像データを取得し、その断層画像データに基づいて物理量の断層分布を演算する制御演算部と、制御演算部の演算結果に基づいて、構造における物理量の分布を断層可視化する態様で表示する表示装置と、を備える。 Another aspect of the present invention is a micro tomography system that visualizes a distribution of a predetermined physical quantity inside a structure to be measured. This system acquires a tomographic image acquisition device for acquiring tomographic images before and after structural deformation, and acquires tomographic image data before and after structural deformation from the tomographic image acquisition device, and based on the tomographic image data, a tomographic distribution of physical quantities. And a display device that displays the distribution of the physical quantity in the structure in a form of visualizing the tomography based on the calculation result of the control calculation unit.
物理量は、構造を対象とする順解析の過程で変位関連ベクトルが演算される所定の数値解析において解析条件又は解析結果となり得るものである。なお、ここでいう「変位関連ベクトル」は、変形ベクトルであってもよいし、変形速度ベクトルであってもよい。制御演算部は、取得した断層画像データに基づき変位関連ベクトルの断層分布を演算し、演算された変位関連ベクトルの断層分布を用いて数値解析の逆解析を実行することにより、物理量の断層分布を演算し、演算された物理量の分布を表示装置に断層可視化させる。 The physical quantity can be an analysis condition or an analysis result in a predetermined numerical analysis in which a displacement-related vector is calculated in the process of forward analysis for a structure. Note that the “displacement-related vector” here may be a deformation vector or a deformation velocity vector. The control calculation unit calculates the tomographic distribution of the displacement-related vector based on the acquired tomographic image data, and performs the inverse analysis of the numerical analysis using the calculated tomographic distribution of the displacement-related vector to obtain the tomographic distribution of the physical quantity. The calculation is performed, and the distribution of the calculated physical quantity is visualized on the display device.
この態様によると、断層画像取得装置にて取得される断層画像から変位関連ベクトルの分布が演算され、その変位関連ベクトル分布に基づき所定の数値解析の逆解析を行うことで物理量の分布が得られる。すなわち、断層画像取得装置による実測と、制御演算部による数値解析の逆解析とのハイブリッドにより、物理量を断層可視化することが可能となる。実測により内部の変位関連ベクトルが得られるため、数値解析の順解析のように境界部に境界条件を設定する必要もない。このため、所望の断層位置について物理量を簡易に算出できる。すなわち、この態様によれば、構造内部における物理量を簡易な手法にてマイクロ断層可視化することができる。 According to this aspect, the distribution of the displacement related vector is calculated from the tomographic image acquired by the tomographic image acquisition device, and the physical quantity distribution is obtained by performing the inverse analysis of the predetermined numerical analysis based on the displacement related vector distribution. .. That is, the physical quantity can be visualized by a hybrid of the actual measurement by the tomographic image acquisition device and the inverse analysis of the numerical analysis by the control calculation unit. Since the internal displacement-related vector is obtained by the actual measurement, it is not necessary to set the boundary condition at the boundary portion unlike the numerical analysis forward analysis. Therefore, the physical quantity can be easily calculated for the desired tomographic position. That is, according to this aspect, it is possible to visualize the physical quantity inside the structure by a simple method.
本発明によれば、構造内部における物理量を簡易な手法にてマイクロ断層可視化可能な技術を提供することができる。 According to the present invention, it is possible to provide a technique capable of visualizing a micro tom by a simple method for a physical quantity inside a structure.
本発明の一実施形態は、構造内部の物理量の分布をマイクロスケールにて断層可視化する方法である。なお、ここでいう「構造」は、材料構造や生体構造(生体組織)を含み得る。「断層可視化」は、少なくともマイクロスケールにて可能であるが、それ以外のスケールを排除するものではない。断層画像取得装置の解像度に応じてナノスケールにて可能なものも含み得る。「物理量」は、所定の数値解析において解析条件又は解析結果の一つとなり得るものであればよく、後述する温度分布など構造の状態変化に係るものでもよいし、物理定数など構造の物性値に係るものでもよい。その場合、「物理定数」は、「弾性係数」,「粘性係数」,「透水係数」等を含み得る。「断層画像」は、光コヒーレンストモグラフィー(Optical Coherence Tomography:以下「OCT」という)によるものでもよい。あるいは、X線CTや超音波断層像によるものでもよい。 One embodiment of the present invention is a method of visualizing the distribution of physical quantities inside a structure on a microscale. The “structure” here may include a material structure and a biological structure (living tissue). "Tomographic visualization" is possible at least on a microscale, but does not exclude other scales. Those that can be performed at the nanoscale may be included depending on the resolution of the tomographic image acquisition apparatus. The “physical quantity” may be one that can be one of the analysis conditions or the analysis results in a predetermined numerical analysis, and may be related to a change in the state of the structure such as a temperature distribution described later, or a physical property value of the structure such as a physical constant. It may be related. In that case, the "physical constant" may include "elastic coefficient", "viscosity coefficient", "permeability coefficient", and the like. The “tomographic image” may be based on optical coherence tomography (hereinafter referred to as “OCT”). Alternatively, X-ray CT or an ultrasonic tomographic image may be used.
この方法は、数値解析の順解析の過程で変位関連ベクトルの分布が演算される点と、断層画像取得装置により実測される断層画像から変位関連ベクトル分布を演算可能である点に着目し、両者の変位関連ベクトルを結び付ける。ここで、「変位関連ベクトル」は、構造内部の変位や変形に関連するベクトルであり、変形ベクトルであってもよいし、変形速度ベクトルであってもよい。「変形ベクトル」は「変位ベクトル」の概念を含む。すなわち、実測と数値解析とを後者の「逆解析」という形で結び付け、両者のハイブリッドによる演算処理を実現している。ここでいう「逆解析」は、順解析の一部を逆方向にたどることにより実現されてよい。 This method pays attention to the point that the distribution of displacement-related vectors is calculated in the process of forward analysis of numerical analysis, and that the displacement-related vector distribution can be calculated from the tomographic image actually measured by the tomographic image acquisition device. Connect the displacement-related vectors of. Here, the “displacement related vector” is a vector related to the displacement or deformation inside the structure, and may be a deformation vector or a deformation velocity vector. “Deformation vector” includes the concept of “displacement vector”. That is, the actual measurement and the numerical analysis are linked in the latter form of "inverse analysis" to realize the arithmetic processing by the hybrid of both. The “inverse analysis” here may be realized by tracing a part of the forward analysis in the reverse direction.
すなわち、順解析においては、例えば上記物理量が所定の数理モデルの解析条件として入力され、その数理モデルに沿った演算処理を経ることで解が算出される。その演算過程で変位関連ベクトルが演算される。「逆解析」は、逆にその変位関連ベクトルの分布を既知として順解析を逆方向にたどり、解析条件の一つを構成する「物理量」の分布を算出する。その変位関連ベクトル分布として、実測された断層画像に基づいて算出されたものが用いられる。なお、「変形ベクトル」の時間微分が「変形速度ベクトル」であるが、「変形速度ベクトル」を変形ベクトルの時間変化とみて「変形ベクトル」の概念で捉えてもよい。いずれにしても、演算過程において変形ベクトルが演算される。本実施形態によれば、OCTやCT等の断層可視化手法から逆解析に結び付ける斬新な手法が提供される。 That is, in the forward analysis, for example, the above-mentioned physical quantity is input as an analysis condition of a predetermined mathematical model, and a solution is calculated by performing arithmetic processing according to the mathematical model. The displacement-related vector is calculated in the calculation process. In the “inverse analysis”, on the contrary, the distribution of the displacement-related vector is known, and the forward analysis is followed in the reverse direction to calculate the distribution of the “physical quantity” that constitutes one of the analysis conditions. As the displacement-related vector distribution, one calculated based on the actually measured tomographic image is used. Note that the time derivative of the “deformation vector” is the “deformation speed vector”, but the “deformation speed vector” may be regarded as the time change of the deformation vector and may be grasped by the concept of the “deformation vector”. In any case, the deformation vector is calculated in the calculation process. According to this embodiment, a novel method for linking back analysis with a tomographic visualization method such as OCT or CT is provided.
この方法では、断層画像取得装置を介して構造の変形前後の断層画像が取得される。なお、「変形前後の断層画像」については、構造の変形前と変形後のそれぞれについて変化が停止している状態で取得してもよい。あるいは、変形中の前後の時点で取得してもよい。すなわち、本手法は、後者のように時間的に変化する現象にも適用可能である。具体的には、OCT断層画像の動画中の2枚の断層画像に対して本手法を適用することにより、それらの断層画像に対応する温度分布を断層可視化することが可能である。前後の断層画像を連続的に取得すれば、温度分布をリアルタイムに断層可視化することもできる。 In this method, tomographic images before and after structural deformation are acquired via a tomographic image acquisition device. The “tomographic images before and after deformation” may be acquired in a state where the change is stopped before and after the deformation of the structure. Alternatively, it may be acquired before and after the deformation. That is, this method can be applied to the latter phenomenon that changes with time. Specifically, by applying this method to two tomographic images in the moving image of the OCT tomographic image, it is possible to visualize the temperature distribution corresponding to these tomographic images. If the front and back tomographic images are continuously acquired, the temperature distribution can be visualized in real time.
そして、その断層画像に基づいて断層位置に対応した変位関連ベクトルの分布が演算される一方、断層画像の内部に対して関心領域(Region of Interest:以下「ROI」と表記する)が設定される。また、上記数値解析の逆解析を利用するためにROIに対して要素分割がなされる。これらROIおよび要素分割の設定は、入力装置を介したユーザの指示入力によりなされる。 Then, the distribution of the displacement-related vector corresponding to the tomographic position is calculated based on the tomographic image, while the region of interest (hereinafter referred to as “ROI”) is set inside the tomographic image. .. Further, element division is performed on the ROI in order to utilize the inverse analysis of the above numerical analysis. The setting of these ROI and element division is made by the user's instruction input through the input device.
断層画像に基づいて演算された変位関連ベクトルは、分割された要素の節点に補間される。この補間後の変形ベクトルの分布に対して上記逆解析を施すことにより、各要素における物理量の分布が演算される。そして、このとき演算された物理量の分布が、表示装置に断層可視化される。 The displacement-related vector calculated based on the tomographic image is interpolated to the nodes of the divided elements. By performing the inverse analysis on the distribution of the deformation vector after the interpolation, the distribution of the physical quantity in each element is calculated. Then, the distribution of the physical quantity calculated at this time is visualized as a tomographic image on the display device.
上記数値解析は、状態変化又は境界条件変化による要因に伴う構造の変形を考慮した応力解析であってもよい。例えば、有限要素法(Finite Element Method:以下「FEM」と表記する)、有限差分法(Finite Difference Method:以下「FDM」と表記する)、あるいは境界要素法(Boundary Element Method:以下「BEM」と表記する)によるものでもよい。 The above-mentioned numerical analysis may be a stress analysis in which the deformation of the structure due to the factors due to the state change or the boundary condition change is taken into consideration. For example, a finite element method (Finite Element Method: hereinafter referred to as “FEM”), a finite difference method (Finite Difference Method: hereinafter referred to as “FDM”), or a boundary element method (Boundary Element Method: hereinafter referred to as “BEM”) Notation).
上記手法は、構造における所定の領域において物理量の真値を与える工程と、逆解析により得られた物理量のうち上記領域に対応する物理量である推定値と、上記領域において与えられた真値との関係に基づき、逆解析により得られた物理量の値を補正する工程と、をさらに備えるものでもよい。 The method is a step of giving a true value of a physical quantity in a predetermined area in the structure, an estimated value which is a physical quantity corresponding to the area in the physical quantity obtained by the inverse analysis, and a true value given in the area. A step of correcting the value of the physical quantity obtained by the inverse analysis based on the relationship may be further provided.
「所定の領域」は、複数設定されてもよく、構造の表面又はその近傍に設定されてもよいし、構造内部の特定の部材の周囲に沿って設定されてもよい。「真値」は、実測値であってもよいし、解析結果(解析値)であってもよく、正しいと推定される値であればよい。例えば、「物理量」が温度であり、赤外線サーモグラフィーや熱電対によりROI近傍の温度値を計測できる場合、その実測値である温度値を真値としてもよい。あるいは、銅配線の電流と電圧によって算出されるジュール熱などに基づいてROI近傍の計測データを算出できる場合、その演算結果を真値としてもよい。この「真値」は、ROI内部の値であることが望ましいが、ROI近傍の温度値(実験計測値)であってもよい。また、熱伝導率などを仮定してモデリングし、ROI内部の任意領域の温度値を推定し、その推定値を「真値」としてもよい。 A plurality of “predetermined regions” may be set, may be set on the surface of the structure or in the vicinity thereof, or may be set along the periphery of a specific member inside the structure. The “true value” may be an actual measurement value, an analysis result (analysis value), or any value that is estimated to be correct. For example, when the “physical quantity” is temperature and the temperature value near the ROI can be measured by infrared thermography or a thermocouple, the temperature value that is the measured value may be the true value. Alternatively, when the measurement data in the vicinity of the ROI can be calculated based on the Joule heat calculated by the current and voltage of the copper wiring, the calculation result may be the true value. The “true value” is preferably a value inside the ROI, but may be a temperature value (experimental measurement value) near the ROI. Further, modeling may be performed assuming thermal conductivity and the like, and a temperature value in an arbitrary region inside the ROI may be estimated, and the estimated value may be set as a “true value”.
上記数値解析は、状態変化又は境界条件変化による要因に伴う構造の変形を考慮した応力解析であってもよい。ここでいう「要因」は、構造への通電であってもよいし、構造の周囲温度の変化であってもよい。あるいは、構造に外力が付与されることでもよい。 The above-mentioned numerical analysis may be a stress analysis in which the deformation of the structure due to the factors due to the state change or the boundary condition change is taken into consideration. The "factor" referred to here may be a current supplied to the structure or a change in ambient temperature of the structure. Alternatively, an external force may be applied to the structure.
例えば、力学的物理量のマイクロ断層可視化に際して下記式(a)を利用する。
[K]{U}={F} ・・・(a)
ここで、[K]は全体剛性マトリクスであり、{U}は要素分割されたROIの各節点における変形ベクトルであり、{F}は外力ベクトルである。上記式(a)において、{F}に関し、ROI内部の節点の荷重はゼロ(合力ゼロ)であり、ROI境界の節点の荷重に対しては変位境界条件(変位拘束条件)を与える。{U}については、変形前後の断層画像に基づいて演算される変位関連ベクトル分布(変形ベクトルの断層分布等)から得られるため、それを代入することができる。上記式(a)を逆解析することにより、物理量としての[K]を推定することができる。
For example, the following equation (a) is used when visualizing a micro-slice of a mechanical physical quantity.
[K]{U}={F} (a)
Here, [K] is the overall stiffness matrix, {U} is the deformation vector at each node of the element-divided ROI, and {F} is the external force vector. In the above equation (a), regarding {F}, the load at the node inside the ROI is zero (zero resultant force), and a displacement boundary condition (displacement constraint condition) is given to the load at the node at the ROI boundary. Since {U} is obtained from the displacement-related vector distribution (such as the deformation vector tomographic distribution) calculated based on the tomographic images before and after deformation, it can be substituted. By inversely analyzing the above equation (a), [K] as a physical quantity can be estimated.
以下、図面を参照しつつ本実施形態を具体化した実施例について詳細に説明する。
[実施例]
図1は、実施例に係る断層可視化システムの構成を概略的に表す図である。本実施例の断層可視化システムは、電子実装部品の内部の特定の物理量をマイクロスケールにて断層可視化するものである。この断層可視化にOCTを利用する。
Hereinafter, examples embodying the present embodiment will be described in detail with reference to the drawings.
[Example]
FIG. 1 is a diagram schematically illustrating a configuration of a tomographic visualization system according to an embodiment. The tomographic visualization system of the present embodiment visualizes a specific physical quantity inside an electronic packaging component on a microscale. OCT is used for this tomographic visualization.
図1に示すように、断層可視化システム1は、OCTを用いる光学系を含む光学ユニット2、光学ユニット2に接続される光学機構4、OCTにより得られた光干渉データに基づいて物理量を演算する制御演算部6、ユーザの指示入力を受け付ける入力装置7、および物理量を断層可視化する態様で表示する表示装置8を備える。光学ユニット2および光学機構4は、「断層画像取得装置」を構成する。図示の例では、光学ユニット2としてマッハツェンダー干渉計をベースとした光学系が示されているが、マイケルソン干渉計その他の光学系を採用することもできる。 As shown in FIG. 1, the tomographic visualization system 1 calculates a physical quantity based on an optical unit 2 including an optical system using OCT, an optical mechanism 4 connected to the optical unit 2, and optical interference data obtained by OCT. A control calculation unit 6, an input device 7 that receives a user's instruction input, and a display device 8 that displays a physical quantity in a manner of visualizing a tomographic image are provided. The optical unit 2 and the optical mechanism 4 form a “tomographic image acquisition device”. In the illustrated example, an optical system based on a Mach-Zehnder interferometer is shown as the optical unit 2, but a Michelson interferometer and other optical systems can also be adopted.
本実施例では、OCTとして波長走査型レーザを光源とするSS−OCT(Swept Source OCT)を用いるが、TD−OCT(Time Domain OCT)、SD−OCT(Spectral Domain OCT)その他のOCTを用いることもできる。SS−OCTは、TD−OCTとは異なり、参照鏡走査などの機械的な光遅延走査を必要としないため、高い時間分解能と高い位置検出精度が得られる点で好ましい。 In this embodiment, SS-OCT (Swept Source OCT) using a wavelength scanning laser as a light source is used as the OCT, but TD-OCT (Time Domain OCT), SD-OCT (Spectral Domain OCT) or other OCT is used. You can also Unlike TD-OCT, SS-OCT does not require mechanical optical delay scanning such as reference mirror scanning, and is therefore preferable in that high time resolution and high position detection accuracy can be obtained.
光学ユニット2は、光源10、オブジェクトアーム12、リファレンスアーム14、および光検出器16を備える。各光学要素は、光ファイバーにて互いに接続されている。光源10から出射された光は、カプラ18(ビームスプリッタ)にて分けられ、その一方がオブジェクトアーム12に導かれる物体光となり、他方がリファレンスアーム14に導かれる参照光となる。オブジェクトアーム12に導かれた物体光は、サーキュレータ20を介して光学機構4に導かれ、測定対象Wに照射される。この物体光は、測定対象Wの表面および断面にて後方散乱光として反射されてサーキュレータ20に戻り、カプラ22に導かれる。 The optical unit 2 includes a light source 10, an object arm 12, a reference arm 14, and a photodetector 16. Each optical element is connected to each other by an optical fiber. The light emitted from the light source 10 is split by the coupler 18 (beam splitter), one of which serves as the object light guided to the object arm 12, and the other serves as the reference light guided to the reference arm 14. The object light guided to the object arm 12 is guided to the optical mechanism 4 via the circulator 20, and is irradiated onto the measurement target W. This object light is reflected as backscattered light on the surface and cross section of the measurement target W, returns to the circulator 20, and is guided to the coupler 22.
一方、リファレンスアーム14に導かれた参照光は、サーキュレータ24を介して反射鏡26に導かれる。このとき、参照光は、コリメートレンズ28を経て集光レンズ30にて反射鏡26に集光される。この参照光は、反射鏡26にて反射されてサーキュレータ24に戻り、カプラ22に導かれる。すなわち、物体光と参照光とがカプラ22にて合波(重畳)され、その干渉光が光検出器16により検出される。 On the other hand, the reference light guided to the reference arm 14 is guided to the reflecting mirror 26 via the circulator 24. At this time, the reference light is condensed on the reflecting mirror 26 by the condenser lens 30 via the collimator lens 28. This reference light is reflected by the reflecting mirror 26, returns to the circulator 24, and is guided to the coupler 22. That is, the object light and the reference light are combined (superposed) by the coupler 22, and the interference light is detected by the photodetector 16.
光学機構4は、光学ユニット2のオブジェクトアーム12を構成し、サーキュレータ20から延びる光ファイバー34が接続されている。光学機構4は、光学ユニット2からの光を測定対象Wに導いて走査させる機構と、その機構を駆動するための駆動部(アクチュエータ)を備える。光学機構4は、例えばガルバノ鏡を含むガルバノ装置であってもよい。光学機構4は、測定対象Wに向けて物体光を照射する。 The optical mechanism 4 constitutes the object arm 12 of the optical unit 2 and is connected to the optical fiber 34 extending from the circulator 20. The optical mechanism 4 includes a mechanism that guides the light from the optical unit 2 to the measurement target W and scans it, and a drive unit (actuator) for driving the mechanism. The optical mechanism 4 may be, for example, a galvanometer device including a galvanometer mirror. The optical mechanism 4 irradiates the object of measurement W with the object light.
測定対象Wは、MEMSに対応した電子実装部品であり、ICチップを基板に埋め込んで構成される。その基板は、高分子樹脂材料等の高分子基材料やプラスチック等の結晶性材料からなるものでもよい。その基板は、内部に反射する構造(界面など)がある材料からなるものであればよく、ナノ粒子などの散乱体が樹脂に内在する複合材料など、透明でない材料からなるものでもよい。 The measurement target W is an electronic mounting component compatible with MEMS, and is configured by embedding an IC chip in a substrate. The substrate may be made of a polymer-based material such as a polymer resin material or a crystalline material such as plastic. The substrate may be made of a material having a structure (interface or the like) that reflects inside, or may be made of a non-transparent material such as a composite material in which a scattering body such as nanoparticles is contained in a resin.
カプラ22にて合波された干渉光は、光検出器16に入力される。光検出器16は、これを光干渉信号(干渉光強度を示す信号)として検出する。この光干渉信号は、A/D変換器40を介して制御演算部6に入力される。A/D変換器40は、光検出器16から出力されたアナログ信号をデジタル信号に変換して制御演算部6へ出力する。 The interference light multiplexed by the coupler 22 is input to the photodetector 16. The photodetector 16 detects this as an optical interference signal (a signal indicating the intensity of interference light). This optical interference signal is input to the control calculation unit 6 via the A/D converter 40. The A/D converter 40 converts the analog signal output from the photodetector 16 into a digital signal and outputs the digital signal to the control calculation unit 6.
制御演算部6は、CPU、ROM、RAM、ハードディスクなどを有し、これらのハードウェア、およびソフトウェアによって、光学ユニット2の光学系全体の制御と、光学機構4の駆動制御と、OCTによる画像出力のための演算処理を行う。制御演算部6の指令信号は、D/A変換器42を介して光学機構4等に入力される。D/A変換器42は、制御演算部6から出力されたデジタル信号をアナログ信号に変換して光学機構4等へ出力する。制御演算部6は、光学機構4の駆動に基づいて光学ユニット2から出力された光干渉信号を処理し、OCTによる測定対象Wの断層画像を取得する。そして、その断層画像データに基づき、後述の手法により測定対象Wの内部における特定の物理量の断層分布を演算する。 The control calculation unit 6 has a CPU, a ROM, a RAM, a hard disk, etc., and controls the entire optical system of the optical unit 2, drive control of the optical mechanism 4, and image output by OCT by these hardware and software. Calculation processing for. The command signal of the control calculation unit 6 is input to the optical mechanism 4 and the like via the D/A converter 42. The D/A converter 42 converts the digital signal output from the control calculation unit 6 into an analog signal and outputs the analog signal to the optical mechanism 4 and the like. The control calculation unit 6 processes the optical interference signal output from the optical unit 2 based on the driving of the optical mechanism 4, and acquires a tomographic image of the measurement target W by OCT. Then, based on the tomographic image data, a tomographic distribution of a specific physical quantity inside the measuring object W is calculated by a method described later.
入力装置7は、キーボードやタッチパネル等からなり、ユーザの指示入力を受け付け、制御演算部6に入力する。表示装置8は、例えば液晶ディスプレイからなり、制御演算部6にて演算された測定対象Wの内部の物理量を断層可視化する態様で画面に表示させる。 The input device 7 includes a keyboard, a touch panel, etc., receives an instruction input by the user, and inputs it to the control calculation unit 6. The display device 8 is composed of, for example, a liquid crystal display, and displays the physical quantity inside the measurement target W calculated by the control calculation unit 6 on the screen in a form of visualizing a tomographic image.
以下、上記物理量の断層分布を可視表示するまでの演算処理方法について説明する。
図2は、実施例に係るマイクロ断層可視化方法の概念を表す図である。本実施例では、測定対象Wの内部の温度を上記物理量としてその分布を断層可視化する。この手法は、測定対象Wの熱変形に基づいて内部の温度分布を演算するものであり、OCTにて取得した断層画像から変形ベクトル分布を演算する工程(変形ベクトル演算工程)と、その変形ベクトル分布に対して有限要素法(以下「FEM」と表記する)の逆解析を施すことにより、測定対象Wの内部の温度分布を演算する工程(逆解析工程)とを有する。
Hereinafter, a calculation processing method for visually displaying the tomographic distribution of the physical quantity will be described.
FIG. 2 is a diagram showing the concept of the micro tomography visualization method according to the embodiment. In the present embodiment, the temperature inside the measuring object W is used as the physical quantity to visualize the distribution of the tomographic image. This method calculates the internal temperature distribution based on the thermal deformation of the measurement target W, and the step of calculating the deformation vector distribution from the tomographic image acquired by OCT (deformation vector calculation step) and the deformation vector thereof. By performing an inverse analysis of a finite element method (hereinafter referred to as “FEM”) on the distribution, there is a step of calculating a temperature distribution inside the measurement object W (inverse analysis step).
まず、変形ベクトル演算工程について説明する。上述のように、OCTにおいて、オブジェクトアーム12を経た物体光(測定対象Wからの反射光)と、リファレンスアーム14を経た参照光とが合波され、光検出器16により光干渉信号として検出される。制御演算部6は、この光干渉信号を干渉光強度に基づく測定対象Wの断層画像として取得することができる。 First, the deformation vector calculation step will be described. As described above, in OCT, the object light that has passed through the object arm 12 (reflected light from the measurement target W) and the reference light that has passed through the reference arm 14 are combined and detected by the photodetector 16 as an optical interference signal. It The control calculation unit 6 can acquire this optical interference signal as a tomographic image of the measurement target W based on the interference light intensity.
ところで、OCTの光軸方向(奥行き方向)の分解能であるコヒーレンス長lcは、光源の自己相関関数によって決定される。ここでは、コヒーレンス長lcを自己相関関数の包括線の半値半幅とし、下記式(1)にて表すことができる。
一方、光軸垂直方向(ビーム走査方向)の分解能は、集光レンズによる集光性能に基づき、ビームスポット径Dの1/2とされる。そのビームスポット径Dは、下記式(2)にて表すことができる。
まず、OCTを利用したマイクロスケールの力学特性検出法について説明する。この検出法は、測定対象Wの変形前後の2枚のOCT断層画像にデジタル相互相関法を適用することで変形ベクトル分布を算出する。この手法によれば、マイクロスケールにて測定対象Wのひずみ断層分布を検出することができる。 First, a microscale mechanical property detection method using OCT will be described. In this detection method, the deformation vector distribution is calculated by applying the digital cross-correlation method to the two OCT tomographic images before and after the deformation of the measurement target W. According to this method, the strain fault distribution of the measurement target W can be detected on a microscale.
変形ベクトル分布を算出する際には、繰り返し相互相関処理を実施する再帰的相互相関法(Recursive Cross-correlation method)を適用する。これは、低解像度において算出された変形ベクトルを参照し、探査領域を限定するとともに階層的に検査領域を縮小して相互相関法を適用する手法である。これにより、高解像度な変形ベクトルを取得することができる。さらに、スペックル・ノイズ低減法として、隣接検査領域の相関値分布との乗算を行う隣接相互相関乗法(Adjacent Cross-correlation Multiplication)を用いる。そして、乗算されることによって高SN化した相関値分布から最大相関値を探索する。 When calculating the deformation vector distribution, a recursive cross-correlation method that performs repetitive cross-correlation processing is applied. This is a method of applying the cross-correlation method by referring to the deformation vector calculated at low resolution, limiting the search area, and hierarchically reducing the inspection area. Thereby, a high-resolution deformation vector can be acquired. Furthermore, as a speckle noise reduction method, Adjacent Cross-correlation Multiplication is used which multiplies the correlation value distribution of the adjacent inspection area. Then, the maximum correlation value is searched for from the correlation value distribution that has been increased in SN by multiplication.
また、マイクロスケールにおける微小変形解析では、変形ベクトルのサブピクセル精度が重要となる。このため、輝度勾配を利用する風上勾配法(Up-stream Gradientmethod)と、伸縮および剪断を考慮した画像変形法(Image Deformation method)の両サブピクセル解析法を併用し、変形ベクトルの高精度検出を実現する。なお、ここでいう「風上勾配法」は、勾配法(オプティカルフロー法)の一種である。 In addition, sub-pixel accuracy of the deformation vector is important in micro deformation analysis on the micro scale. Therefore, the up-stream gradient method that uses the brightness gradient and the image deformation method that considers expansion and contraction (Image Deformation method) are used together to detect the deformation vector with high accuracy. To achieve. The “upwind gradient method” here is a kind of gradient method (optical flow method).
このようにして得られた変形ベクトル分布を空間微分することにより、ひずみ断層分布を演算することができる。なお、本実施例では後述のように、変形ベクトル分布に対してFEMの逆解析を施すことにより、測定対象Wの内部の温度分布を演算する。このため、ひずみ断層分布まで演算する必要はなく、その演算過程である変形ベクトル分布の演算までが利用されることとなる。 The strain fault distribution can be calculated by spatially differentiating the deformation vector distribution thus obtained. In this embodiment, the temperature distribution inside the measurement object W is calculated by performing the FEM inverse analysis on the deformation vector distribution as described later. Therefore, it is not necessary to calculate the strain fault distribution, and the calculation of the deformation vector distribution, which is the calculation process, is used.
以下、各手法について詳細に説明する。図3は、再帰的相互相関法による処理手順を概略的に示す図である。図4は、サブピクセル解析による処理手順を概略的に示す図である。 Hereinafter, each method will be described in detail. FIG. 3 is a diagram schematically showing a processing procedure by the recursive cross-correlation method. FIG. 4 is a diagram schematically showing a processing procedure by subpixel analysis.
(再帰的相互相関法)
図3(A)〜(C)は、再帰的相互相関法による処理過程を示している。各図にはOCTにより撮影される前後の断層画像が示されている。左側には先の断層画像(Image1)が示され、右側には後の断層画像(Image2)が示されている。
(Recursive cross-correlation method)
3(A) to 3(C) show a processing process by the recursive cross-correlation method. Each figure shows tomographic images before and after being photographed by OCT. The previous tomographic image (Image1) is shown on the left side, and the subsequent tomographic image (Image2) is shown on the right side.
相互相関法とは、局所的なスペックル・パターンの類似度を下記式(3)に基づく相関値Rijより評価する方法である。そのため、図3(A)に示すように、撮影された前後のOCT画像について、先の断層画像(Image1)に類似度の検査対象となる検査領域S1が設定され、後の断層画像(Image2)に類似度の探査範囲となる探査領域S2が設定される。
探査領域S2(Mx×Mzピクセル)内における相関値分布Ri,j(ΔX,ΔZ)を算出し、図3(B)に示すようにパターンマッチングを行う。実際には、下記式(4)に示すように、最大相関値を与える移動量Ui,jを変形前後の変形ベクトルとして決定する。
本手法では、検査領域S1を縮小しながら相互相関処理を繰り返して空間解像度を高める再帰的相互相関法を採用している。なお、本実施例では解像度を上げる際に空間解像度が倍になるようにしている。図3(C)に示すように、検査領域S1を1/4に分割し、前階層にて算出された変形ベクトルを参照し、探査領域S2を縮小している。ここで探査領域S2も1/4に分割している。再帰的相互相関法を用いることで、高解像度において多発する過誤ベクトルの抑制を可能にしている。このような再帰的相互相関処理を施すことにより、変形ベクトルの解像度を高めることができる。 This method employs a recursive cross-correlation method in which the cross-correlation process is repeated while reducing the inspection area S1 to increase the spatial resolution. In this embodiment, the spatial resolution is doubled when the resolution is increased. As shown in FIG. 3C, the inspection area S1 is divided into quarters, the deformation vector calculated in the previous layer is referred to, and the search area S2 is reduced. Here, the search area S2 is also divided into 1/4. By using the recursive cross-correlation method, it is possible to suppress error vectors that frequently occur at high resolution. By performing such recursive cross-correlation processing, the resolution of the deformation vector can be increased.
また、これに加え、下記式(5)により、演算中の座標を中心とする周囲8つの座標を含む合計9つの変形ベクトルの平均偏差σを用いた閾値を設定し、過誤ベクトルの除去を行い再帰処理に伴う誤差伝播を抑制する。
(隣接相互相関乗法)
本実施例では、スペックルノイズの影響を受けたランダム性の強い相関値分布から正確な最大相関値を決定する方法として、隣接相互相関乗法を導入している。下記式(6)により、隣接相互相関乗法では検査領域S1における相関値分布Ri,j(Δx,Δz)と、その検査領域S1にオーバーラップする隣接検査領域に対する相関値分布Ri+Δi,j(Δx,Δz)とRi,j+Δj(Δx,Δz)の乗算を行い、新たな相関値分布R'i,j(Δx,Δz)を用いて最大相関値を検索する。
In the present embodiment, the adjacent cross-correlation multiplication method is introduced as a method of accurately determining the maximum correlation value from the correlation value distribution having strong randomness affected by speckle noise. According to the following equation (6), the correlation value distribution R i,j (Δx, Δz) in the inspection area S1 and the correlation value distribution R i +Δi,j (for the adjacent inspection area overlapping the inspection area S1 in the adjacent cross correlation multiplication method Δx, Δz) is multiplied by R i,j+Δj (Δx, Δz), and the maximum correlation value is searched using the new correlation value distribution R′ i,j (Δx, Δz).
これにより、相関値同士の乗算によってランダム性を低減させることが可能になる。上述した検査領域S1の縮小と共に干渉強度分布の情報量も減少するため、スペックル・ノイズを原因とする複数相関ピークの出現が計測精度の悪化を招いていると考えられる。一方、隣接境界同士の移動量には相関があるため、最大相関値座標付近では強い相関値が残存する。この隣接相互相関乗法の導入によって最大相関値ピークが明瞭化され計測精度が向上し、正確な移動座標を抽出することが可能となる。また、この隣接相互相関乗法をOCTの各ステージに導入することで、誤差伝播が抑制され、スペックル・ノイズに対するロバスト性が向上する。それにより、高空間解像度においても高精度な変形ベクトル分布(移動量分布)の算出が可能となる。 This makes it possible to reduce the randomness by multiplying the correlation values. Since the amount of information of the interference intensity distribution decreases with the reduction of the inspection area S1 described above, it is considered that the appearance of a plurality of correlation peaks due to speckle noise deteriorates the measurement accuracy. On the other hand, since the movement amounts of the adjacent boundaries have a correlation, a strong correlation value remains near the maximum correlation value coordinate. By introducing this adjacent cross-correlation multiplication method, the maximum correlation value peak is clarified, the measurement accuracy is improved, and accurate moving coordinates can be extracted. Also, by introducing this adjacent cross-correlation multiplication method in each stage of OCT, error propagation is suppressed and robustness against speckle noise is improved. Thereby, it is possible to calculate the deformation vector distribution (movement amount distribution) with high accuracy even at high spatial resolution.
(風上勾配法)
図4(A)〜(C)は、サブピクセル解析による処理過程を示している。各図にはOCTにより連続的に撮影される前後の断層画像が示されている。左側には先の断層画像(Image1)が示され、右側には後の断層画像(Image2)が示されている。
(Upwind gradient method)
FIGS. 4A to 4C show the process of subpixel analysis. Each figure shows tomographic images before and after being continuously photographed by OCT. The previous tomographic image (Image1) is shown on the left side, and the subsequent tomographic image (Image2) is shown on the right side.
本実施例では、サブピクセル解析のために風上勾配法と画像変形法を採用する。最終的な移動量の算出は後述の画像変形法によるが、計算の収束性の問題から、画像変形法より先に風上勾配法を適用する。検査領域サイズが小さく高空間解像度の条件において、サブピクセル移動量を高精度検出する画像変形法及び風上勾配法を適用している。画像変形法におけるサブピクセル移動量の検出が困難な場合において、風上勾配法によりサブピクセル移動量を算出する。 In this embodiment, the upwind gradient method and the image transformation method are adopted for subpixel analysis. The final amount of movement is calculated by the image modification method described later, but due to the problem of convergence of calculation, the upwind gradient method is applied before the image modification method. Under the condition that the inspection area size is small and the spatial resolution is high, the image transformation method and the upwind gradient method that detect the sub-pixel movement amount with high accuracy are applied. When it is difficult to detect the subpixel movement amount in the image transformation method, the subpixel movement amount is calculated by the upwind gradient method.
サブピクセル解析では、注目点における変形前後の輝度差が各成分の輝度勾配と移動量によって表される。このため、検査領域S1内の輝度勾配データより最小二乗法を用いてサブピクセル移動量を決定することができる。本実施例では、輝度勾配を求める際に、サブピクセル変形前の風上側の輝度勾配を与える風上差分法を採用している。すなわち、サブピクセル解析は様々な手法が存在するが、本実施例では検査領域サイズが小さく高空間解像度の条件においても、サブピクセル移動量を高精度検出する勾配法を採用している。 In the sub-pixel analysis, the brightness difference before and after the deformation at the point of interest is represented by the brightness gradient and the movement amount of each component. Therefore, the sub-pixel movement amount can be determined from the brightness gradient data in the inspection area S1 by using the least square method. In the present embodiment, when obtaining the brightness gradient, the upwind difference method which gives the brightness gradient on the upwind side before the sub-pixel deformation is adopted. That is, although there are various methods for subpixel analysis, the present embodiment employs the gradient method for detecting the subpixel movement amount with high accuracy even under the condition that the inspection region size is small and the spatial resolution is high.
風上勾配法は、検査領域S1内の注目点の移動を、図4(A)に示すピクセル精度に留まらず、図4(B)に示すサブピクセル精度にて算出するものである。なお、図中の各格子は1ピクセルを表している。実際には図示の断層画像と比較して相当小さいが、説明の便宜上、大きく表記している。この風上勾配法は、微小変形前後における輝度分布の変化を輝度勾配と移動量によって定式化する手法であり、fを輝度とすると、微小変形f(x+Δx,z+Δz)をテイラー展開する下記式(7)として表される。
上記式(7)は、注目点の変形前後の輝度差が変形前の輝度勾配と移動量によって表されることを示している。なお、移動量(Δx,Δz)については上記式(7)のみでは決定できないため、検査領域S1内で移動量が一定と考え、最小二乗法を適用して算出している。 The above formula (7) indicates that the difference in brightness before and after the deformation of the target point is represented by the brightness gradient and the amount of movement before the deformation. Since the movement amount (Δx, Δz) cannot be determined only by the above equation (7), the movement amount is considered to be constant in the inspection region S1 and is calculated by applying the least squares method.
上記式(7)を用いて移動量を算出する際には、右辺の各注目点における移動前後の輝度差は一意にしか求まらない。そのため輝度勾配をどれだけ正確に算出するかが移動量の精度に直結する。輝度勾配の差分化では、一次精度風上差分を用いている。差分化において高次差分を適用すると、多くのデータが必要になり、ノイズが含まれていた際に影響を大きく受けてしまうためである。また、検査領域S1内の各点を基準とした高次差分では、検査領域S1外のデータを多く使用することとなり、検査領域S1そのものの移動量ではなくなってしまうという問題点も存在するからである。 When the movement amount is calculated using the above equation (7), the brightness difference before and after the movement at each attention point on the right side can be uniquely obtained. Therefore, how accurately the brightness gradient is calculated is directly related to the accuracy of the movement amount. The first-order accuracy windward difference is used for the difference in the brightness gradient. This is because if a high-order difference is applied in the subtraction, a large amount of data is needed and it will be greatly affected when noise is included. In addition, in the high-order difference using each point in the inspection area S1 as a reference, a large amount of data outside the inspection area S1 is used, and there is a problem that the movement amount of the inspection area S1 itself is lost. is there.
輝度勾配を求める際に変形前の風上側の輝度勾配が移動することによって注目点の輝度差が生まれると考えることができるので、変形前は風上側の差分を適用する。ここでいう風上は、実際の移動方向ではなく、ピクセル移動量に対するサブピクセル移動量の向きのことであり、最大相関値ピークに放物線近似を施すことによって風上側を決定する。逆に、変形後の風下側の輝度勾配が逆に移動することによって注目点の輝度差が生じると考えることができるため、変形後は風下側の差分を適用する。 When obtaining the brightness gradient, it can be considered that the brightness difference on the windward side before deformation moves to cause a brightness difference at the point of interest, and therefore the difference on the windward side is applied before the deformation. The upwind here is not the actual moving direction but the direction of the subpixel moving amount with respect to the pixel moving amount, and the upwind side is determined by performing parabolic approximation on the maximum correlation value peak. On the contrary, since it can be considered that the difference in brightness at the point of interest occurs due to the reverse movement of the brightness gradient on the leeward side after the deformation, the difference on the leeward side is applied after the deformation.
変形前の風上差分と変形後の風下差分を用いて2通りの解を求め、それらの平均をとった。さらに、実際には移動量が軸方向に沿わない場合には、変形前や変形後の輝度勾配が注目点と同一軸上に無く、ずれた位置の勾配を求める必要がある。つまり、X方向の勾配を求めたい場合には、Z方向の移動も考慮して勾配を求めなければならない。そのため、輝度の内挿による輝度勾配の推定をすることで、精度向上を図っている。基本的には変形前(もしくは変形後)の位置を予測し、その位置での勾配を内挿により求める。 Two solutions were obtained using the upwind difference before the deformation and the downwind difference after the deformation, and the average of them was taken. Furthermore, when the amount of movement does not actually follow the axial direction, there is no brightness gradient before and after deformation on the same axis as the point of interest, and it is necessary to obtain a gradient at a shifted position. That is, in order to obtain the gradient in the X direction, the gradient must be obtained in consideration of the movement in the Z direction. Therefore, the accuracy is improved by estimating the brightness gradient by the brightness interpolation. Basically, the position before (or after) deformation is predicted, and the gradient at that position is obtained by interpolation.
変形前(後)の注目点の位置は放物線近似を施した際のサブピクセル移動量により求める。その注目点位置が囲まれる8つの座標を用い、それらの比によって輝度勾配を算出する。具体的には、下記式(8)を用いる。そのようにして算出された輝度勾配と、輝度変化を用いて最小二乗法を適用し移動量を決定した。
(画像変形法)
上述した風上勾配法までは検査領域S1の形状は変更せず、正方形を保ったまま変形ベクトルの算出を行っている。しかし、現実には計測対象の変形に合わせて検査領域S1も変形していると考えられるため、検査領域S1の微小変形を考慮したアルゴリズムを導入し、変形ベクトル算出を高精度にて算出する必要がある。このため、本実施例ではサブピクセル精度での変形ベクトルの算出に画像変形法を導入している。すなわち、材料の変形前の検査領域S1と変形後の伸縮及びせん断変形を考慮した検査領域S1とで相互相関を実施し、相関値ベースの反復計算によってサブピクセル変形量を決定している。なお、検査領域S1の伸縮及びせん断変形は線形で近似している。
(Image transformation method)
The shape of the inspection region S1 is not changed until the upwind gradient method described above, and the deformation vector is calculated while keeping the square shape. However, in reality, it is considered that the inspection region S1 is also deformed in accordance with the deformation of the measurement target. Therefore, it is necessary to introduce an algorithm that takes into account minute deformation of the inspection region S1 and calculate the deformation vector with high accuracy. There is. For this reason, in this embodiment, the image transformation method is introduced to calculate the transformation vector with sub-pixel accuracy. That is, the cross-correlation is performed between the inspection region S1 before the material is deformed and the inspection region S1 in which the expansion and contraction and the shear deformation after the deformation are considered, and the sub-pixel deformation amount is determined by the iterative calculation based on the correlation value. The expansion and contraction and the shear deformation of the inspection area S1 are linearly approximated.
画像変形法は一般的に、材料の表面ひずみ計測法に用いられ、ランダムパターンを塗布した材料表面を高空間分解能なカメラで撮影した画像に対して適用される。一方、OCT断層像はスペックルノイズを多く含み、また、複合材料での局所的な変形が発生するため、検査領域S1に変形が伴う場合には解の収束解が著しく低下する。本手法における検査領域S1の縮小は、局所的な力学特性を検出するために必要不可欠である。そこで画像変形法では、風上勾配法で得られた変形量を収束計算の初期値として採用し、さらに、輝度分布の双三次関数補間によって検査領域S1を縮小した際でも低ロバスト性を実現している。なお、変形例においては、双三次関数補間以外の補間関数を用いてもよい。 The image deformation method is generally used for measuring the surface strain of a material, and is applied to an image of a material surface coated with a random pattern taken by a camera with high spatial resolution. On the other hand, the OCT tomographic image contains a large amount of speckle noise, and since local deformation occurs in the composite material, when the inspection region S1 is deformed, the convergent solution is significantly reduced. The reduction of the inspection area S1 in this method is indispensable for detecting local mechanical characteristics. Therefore, in the image deformation method, the deformation amount obtained by the upwind gradient method is adopted as the initial value of the convergence calculation, and further, the low robustness is realized even when the inspection area S1 is reduced by the bicubic function interpolation of the luminance distribution. ing. In the modification, an interpolation function other than the bicubic function interpolation may be used.
より詳細には、以下の手順にて演算を実行する。まず、材料変形前のOCT断層像の輝度分布に双三次関数補間法を適用し、輝度分布の連続化を実施する。双三次関数補間法とは、sinc関数を区分的に三次関数近似した畳み込み関数を用い、輝度情報の空間連続性を再現する手法である。本来は連続的な輝度分布を画像計測する際には光学系に依存した点広がり関数が畳み込まれるため、sinc関数を用いた逆畳み込みを行うことにより、本来の連続的な輝度分布が復元される。離散的な一軸信号f(x)の補間をする場合、畳み込み関数h(x)は下記式(9)にて表される。
なお、OCT計測条件の違いによって輝度補間関数h(x)の形状も変更する必要がある。そこで輝度補間関数h(x)のx=1での微係数aを可変とし、aの値を変更することで輝度補間関数h(x)の形状を変更可能なアルゴリズムとした。本実施例では、擬似OCT断層像を用いた数値実験による検証結果を元にし、aの値を決定した。以上のように画像補間をすることで、伸縮及びせん断変形を考慮した検査領域S1の各点にて、OCT輝度値を求めることが可能となる。 The shape of the brightness interpolation function h(x) needs to be changed depending on the difference in OCT measurement conditions. Therefore, the differential coefficient a at x=1 of the luminance interpolation function h(x) is made variable, and the shape of the luminance interpolation function h(x) can be changed by changing the value of a. In this example, the value of a was determined based on the verification result by the numerical experiment using the pseudo OCT tomographic image. By performing the image interpolation as described above, it becomes possible to obtain the OCT luminance value at each point of the inspection region S1 in consideration of expansion and contraction and shear deformation.
伸縮及びせん断変形を考慮して算出した検査領域S1は、図4(C)に示すように、移動とともに変形を伴う。材料変形前のOCT断層像におけるある検査領域S1内の整数ピクセル位置での座標(x,z)が変形後に座標(x*,z*)に移動すると考えると、x*,z*の値は下記式(10)にて表される。
ここで、Δ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)にて表される。
収束の判定には、反復計算で随時得られる漸近解が収束解の近傍で十分小さくなることを用いる。しかし、スペックルパターンの変化が激しい領域においては、線形変形では追従できないために正しい収束解が得られない場合がある。その場合、本実施例では風上勾配法によって求めたサブピクセル移動量を採用している。以上のようにして、サブピクセル精度の変形ベクトル分布が得られる。 For the determination of convergence, it is used that the asymptotic solution obtained by iterative calculation is sufficiently small near the convergent solution. However, in a region where the speckle pattern changes drastically, it may not be possible to obtain a correct convergent solution because linear transformation cannot follow it. In this case, in this embodiment, the subpixel movement amount obtained by the upwind gradient method is used. As described above, the deformation vector distribution with sub-pixel accuracy is obtained.
次に、逆解析工程について説明する。
図5は、有限要素モデルを概略的に示す図である。(A)は三角形要素を示し、(B)は簡易モデルを示す。図6〜図9は、逆解析の処理手順を概略的に示す図である。各図の(A)〜(C)は、その処理過程を示す。図8は、図7(C)の部分拡大図に対応する。
Next, the inverse analysis step will be described.
FIG. 5 is a diagram schematically showing a finite element model. (A) shows a triangular element, (B) shows a simple model. 6 to 9 are diagrams schematically showing the procedure of inverse analysis processing. (A) to (C) of each figure show the processing steps. FIG. 8 corresponds to a partially enlarged view of FIG.
本実施例では、OCT断層像から求められる変形ベクトル断層分布とFEMとを組み合わせることで温度の断層分布を逆解析推定する。図6(A)には、OCTにより得られる電子実装部品(測定対象W)の断層画像が示されている。横軸がx座標を示し、縦軸がz座標を示す。同図には、高分子基材料からなる基板にICチップが埋設された電子実装部品の一部が示されている。 In this embodiment, the deformation vector tomographic distribution obtained from the OCT tomographic image and the FEM are combined to inversely estimate the temperature tomographic distribution. FIG. 6A shows a tomographic image of the electronic mounting component (measurement target W) obtained by OCT. The horizontal axis represents the x coordinate and the vertical axis represents the z coordinate. The figure shows a part of an electronic mounting component in which an IC chip is embedded in a substrate made of a polymer-based material.
図6(B)には、OCTにより得られた前後の断層画像を用いて上記演算処理を行うことにより得られる変形ベクトル分布が矢印にて示されている。このような変形ベクトル分布に対してFEMの逆解析を施すことにより、測定対象内部の温度分布を算出し、断層表示することができる。 In FIG. 6(B), the deformation vector distribution obtained by performing the above-described arithmetic processing using the tomographic images before and after obtained by OCT is indicated by arrows. By performing an FEM inverse analysis on such a deformation vector distribution, the temperature distribution inside the measurement target can be calculated and displayed in a tomographic manner.
ここで、「逆解析」の意義を明確にするために、FEMの順解析について理解しておく必要がある。この順解析では、解析対象となる領域(解析領域)を要素分割した後、材料パラメータを設定するとともに荷重等の境界条件を設定し、各要素における温度を与えた後、多元連立方程式となる全体剛性方程式を解く。それにより、各節点における変位(未知項)を算出することができる。なお、本実施例では温度分布と変形分布とが連成するため、順解析は次のようになる。すなわち、解析領域を要素分割した後、熱伝導係数分布を設定するとともに熱(温度)の境界条件を設定し、方程式を解いて温度分布を算出する過程を経る。また、解析領域を要素分割した後、材料パラメータ(剛性、線膨張係数など),力学的境界条件,温度分布を設定し、方程式を解いて変形分布を算出する過程を経る。 Here, in order to clarify the significance of “inverse analysis”, it is necessary to understand the forward analysis of FEM. In this forward analysis, after the area to be analyzed (analysis area) is divided into elements, material parameters are set, boundary conditions such as loads are set, and the temperature at each element is given. Solve the stiffness equation. Thereby, the displacement (unknown term) at each node can be calculated. In this example, since the temperature distribution and the deformation distribution are coupled, the forward analysis is as follows. That is, after the analysis region is divided into elements, the process of calculating the temperature distribution by setting the heat conduction coefficient distribution, setting the heat (temperature) boundary condition, and solving the equation. Further, after the analysis region is divided into elements, a process of setting material parameters (rigidity, coefficient of linear expansion, etc.), mechanical boundary conditions, temperature distribution, solving equations and calculating deformation distribution is performed.
これに関し、本実施例では上述のように、OCTの断層画像に基づいて変形ベクトル分布が算出されるため、解析領域に設定された節点に対してそれらの変形ベクトルを補間することにより、境界条件を設定しなくとも節点上の変位を算出することができる。すなわち、全体剛性方程式における変位ベクトル(変形ベクトル)の項が既知となる。この点に着目し、この既知となった変位ベクトルを用いてFEMの逆解析を実行することにより温度項を算出する。このとき、全節点における変形ベクトルが既知であるため、全要素について温度項を算出することができる。つまり、全要素に対応した温度分布を算出することができる。以下、具体的に説明する。 In this regard, in the present embodiment, as described above, the deformation vector distribution is calculated based on the tomographic image of the OCT, and therefore, by interpolating the deformation vectors with respect to the nodes set in the analysis region, the boundary condition is calculated. The displacement on the node can be calculated without setting. That is, the term of the displacement vector (deformation vector) in the overall stiffness equation becomes known. Focusing on this point, the temperature term is calculated by performing the inverse analysis of the FEM using the known displacement vector. At this time, since the deformation vectors at all nodes are known, the temperature terms can be calculated for all the elements. That is, the temperature distribution corresponding to all the elements can be calculated. The details will be described below.
(FEMによる温度項を含む全体剛性方程式)
二次元平面応力問題に関し、熱膨張を考慮した応力-ひずみ関係式として下記式(12)が得られる。なお、ここではOCTの断層画像に合わせてx−z平面を考える。
Regarding the two-dimensional plane stress problem, the following equation (12) is obtained as a stress-strain relational equation considering thermal expansion. Note that, here, the xz plane is considered according to the tomographic image of OCT.
そして、上記式(12)に仮想仕事の原理を適用して変形し、要素剛性マトリクス[k]6×6を用いて表すと、下記式(13)が得られる。
本実施例ではFEMを利用するに際し、要素分割される各要素を図5(A)に示す三角形要素にて表し、また、隣接する要素間の関係を図5(B)に示す簡易モデルにて表す。上記式(13)において、P,Qは等価節点力であり、u,vはそれぞれ節点位置のx方向,z方向の変形ベクトル成分である。なお、P,Q,u,vの添え字1〜3は、図5(A)に示す要素分割に係る三角形要素の局所節点番号である。a1〜b3は三角形要素において物理量を近似する際の定数であり、それぞれ下記式(14)により表される。
本実施例では、図7(A)に示すように、測定対象の内部にROIを設定し、図7(B)に示すように、そのROIに対してFEMに基づく要素分割を実行する。OCTにより得られた変形ベクトルは、図7(C)に示すようにROIに分布する。 In this embodiment, an ROI is set inside the measurement target as shown in FIG. 7A, and element division based on FEM is executed for the ROI as shown in FIG. 7B. The deformation vector obtained by OCT is distributed in the ROI as shown in FIG.
上記式(13)を温度の項について展開し、図5(B)に示す簡易モデルを用いて全体方程式を作成すると、下記式(15)が成立する。
上記式(15)を簡略化および一般化すると下記式(16)となり、線形化された観測方程式となる。
ただし、minはROI境界上を除くROI内の節点数であり、2minはそれらの節点上に求まる変形ベクトルの成分数である。mbはROI境界上の節点数であり、2mbはそれらの節点上に求まる変形ベクトルの成分数である。したがって、2min+2mbはROI以内の全節点に求まる変形ベクトルの成分数となる。nはROI内の要素数である。上記式(16)の連立方程式に境界変位条件を導入することで修正した全体剛性方程式に対し、種々の逆解析手法を用いることによって、各要素における温度変化ΔTを推定することができる。 However, m in is the number of nodes in the ROI except on the ROI boundary, and 2 m in is the number of components of the deformation vector obtained on those nodes. m b is the number of nodes on the ROI boundary, and 2 m b is the number of components of the deformation vector obtained on those nodes. Therefore, 2m in +2m b is the number of components of the deformation vector obtained at all nodes within the ROI. n is the number of elements in the ROI. The temperature change ΔT in each element can be estimated by using various inverse analysis methods for the overall stiffness equation modified by introducing the boundary displacement condition into the simultaneous equations of the above equation (16).
(境界変位条件の導入)
ROI内の有限要素モデルに上記式(16)を適用した後、各節点における変形ベクトル成分{U}(2min+2mb)×1を節点補間により与える。このとき、ROI境界上の節点に得られる変形ベクトル成分u,vを境界強制変位(拘束節点)とする必要がある。そのため、上記式(16)における全体剛性マトリクス[K](2min+2mb)×(2min+2mb)、および観測マトリクス[A](2min+2mb)×nを修正してマトリクスの縮小を行う。
(Introduction of boundary displacement condition)
After applying the above equation (16) to the finite element model in the ROI, the deformation vector component {U} (2min+2mb)×1 at each node is given by the node interpolation. At this time, it is necessary to use the deformation vector components u and v obtained at the nodes on the ROI boundary as boundary forced displacements (restraint nodes). Therefore, the overall rigidity matrix [K] (2min+2mb)×(2min+2mb) and the observation matrix [A] (2min+2mb)×n in the above equation (16) are corrected to reduce the matrix.
すなわち、マトリクス成分の添え字の意味を改変した上で上記式(16)を改めて書き下すと下記式(17)となる。
この全体剛性方程式を解く過程で境界変位条件を満足するように修正する。ここで、ROI境界上の節点に得られた変形ベクトル成分をul,vlとする。このとき、{U}(2min+2mb)×1において境界上の節点に求まる変形ベクトル成分を表す行はulが2l−1行目、vlが2l行目となる。ROI境界の節点では境界変位条件(変位拘束境界)として、ROI境界上節点の既知の変位量ul,vlについて左右両辺が満たされるように、全体剛性マトリクス[K](2min+2mb)×(2min+2mb)、および観測マトリクス[A](2min+2mb)×nを修正する。具体的には、ROI境界上節点における変位量ul,vlが,左辺と右辺においてul=ul,vl=vlとなるように、[K]の成分と、[A]における同じ行の成分を下記式(18)に従って修正する。
また、ベクトル{F}の境界上の節点に対応する2l−1,2l行目も、それぞれul,vlと修正して左右両辺が満たされるようにする。さらに、ROIが物体内部に設定され、かつ静的変形により体積力が作用しないのであれば、{F}に関してROI境界上を除くROI内の等価節点力はゼロとなる。また、ROI境界を自由表面に設定した場合は自由変形となるため、表面力は{F}=0となる。このため、上記式(17)は下記式(19)のようになる。
以上のようにして、全体剛性方程式に境界変位条件を満足させる。なお、上記式(19)において、両辺の2l−1,2l行目は拘束節点に帰属し、境界変位条件を満たすため、それらを省略することで成分を並び替え、マトリクスの縮小を行う。このようにして得られる観測マトリクス[A]を[A']と表し、全体剛性マトリクス[K]と変形ベクトル{U}の積を{K'U'}と表すと、下記式(20)が得られる。
(節点補間)
この逆解析のためには、上記式(16)における測定項{U}を与えなければならない。一方、図8(A)に示すように、OCTにより得られた変形ベクトル分布は、ROI内に分割された要素の節点に対応していないため、これをそのまま測定項{U}とすることはできない。そこで、図8(B)に示すようにROIの節点上に変形ベクトルを設定するための補間処理を実行する。
(Nodal interpolation)
For this inverse analysis, the measurement term {U} in the above equation (16) must be given. On the other hand, as shown in FIG. 8(A), since the deformation vector distribution obtained by OCT does not correspond to the nodes of the elements divided in the ROI, this may not be used as it is as the measurement term {U}. Can not. Therefore, as shown in FIG. 8B, an interpolation process for setting a deformation vector on the node of the ROI is executed.
すなわち、OCTにより取得した変形ベクトルの断層分布を用いて各節点上に変形ベクトルを補間する。ここでは、各節点を中心とする所定半径(本実施例では300μm)を「影響半径」とし、その影響半径に収まる変形ベクトルをその節点に補間した。この補間に際しては、最小二乗法により下記式(21)に示す近似関数f(x,z)を求め、その近似関数f(x,z)に則って各節点位置の変形ベクトルを算出した。
この近似関数f(x,z)について下記式(22)を適用する。
具体的には、Qを算出した後、Qを近似関数の各未知係数について偏微分することにより、下記式(23)に示す正規式を得る。
この正規式を解き、未知係数を算出することで近似関数を求める。求まった近似関数に節点座標を参照させ、節点位置での変形量成分を求めることで節点補間がなされる。ただし、表面近傍のROI境界や深部のROI境界では輝度情報にノイズが多く含まれるため、変形ベクトルの算出精度が低い。さらに、ROI境界近傍、特にROIの頂点部分(角部分)では使用できる変形ベクトルの数が少ないため、ROI境界外に変形ベクトルを外挿補間している。このような手法により、ROI内の節点上の各変形量成分を求めて節点補間をする。なお、上記式(16)の{U}は、OCTによる測定データとして得られる。 An approximate function is obtained by solving this normal expression and calculating an unknown coefficient. Nodal interpolation is performed by referring to the nodal coordinates in the found approximate function and finding the deformation amount component at the nodal position. However, since a lot of noise is included in the luminance information at the ROI boundary near the surface and the ROI boundary at the deep part, the calculation accuracy of the deformation vector is low. Furthermore, since the number of deformation vectors that can be used is small in the vicinity of the ROI boundary, particularly in the vertex portion (corner portion) of the ROI, the deformation vector is extrapolated outside the ROI boundary. By such a method, each deformation amount component on the node in the ROI is obtained and the node interpolation is performed. Note that {U} in the above equation (16) is obtained as measurement data by OCT.
(逆解析)
上記式(20)における観測マトリクス[A']は一般に2min×nの行列である。minはROI境界を除いたROI内部の節点x,z座標の数であり、nはROI内の要素数である。min<nとなるため、[A']は正則ではない。そのため、本実施例では、逆解析手法として擬似逆行列であるムーアペンローズ一般逆行列を求め、逆解析を行った。ここで、観測マトリクス[A']がフルランクでランクがnの場合、[A']の一般逆行列は下記式(24)にて表される。
The observation matrix [A′] in the above equation (20) is generally a matrix of 2m in ×n. m in is the number of node x and z coordinates inside the ROI excluding the ROI boundary, and n is the number of elements in the ROI. [A′] is not regular because m in <n. Therefore, in this embodiment, a Moore-Penrose generalized inverse matrix, which is a pseudo-inverse matrix, is obtained as an inverse analysis method, and the inverse analysis is performed. Here, when the observation matrix [A′] has a full rank and the rank is n, the general inverse matrix of [A′] is represented by the following equation (24).
上記式(24)に基づき、温度の逆解析を下記式(25)を用いて実行する。
ただし、ムーアペンローズ一般逆行列により算出された推定値は規格化されており、物理量として単位をもたない値となる。そこで、この推定値を有用な物理量の温度値に較正する。すなわち、図9(A)に点線枠にて示すように、推定値の一部に対応する温度の真値を取得し、その真値に対応する推定値(「対応推定値」ともいう)とその真値とに基づいて、全推定値を線形補間することにより較正を行う。この真値は実際に計測した値(実計測値)であってよい。 However, the estimated value calculated by the Moore-Penrose generalized inverse matrix is standardized and has a unitless value as a physical quantity. Therefore, this estimated value is calibrated to a temperature value of a useful physical quantity. That is, as shown by a dotted line frame in FIG. 9A, a true value of temperature corresponding to a part of the estimated value is acquired, and an estimated value (also referred to as “corresponding estimated value”) corresponding to the true value is obtained. Calibration is performed by linearly interpolating all estimates based on their true values. This true value may be a value actually measured (actual measured value).
すなわち、逆解析にて得た温度分布を、下記式(26)に基づいて較正する。
制御演算部6は、以上のような演算処理を実行し、算出された温度値T* realizedを表示装置8に断層可視化する。それにより、例えば図9(C)に示すような色分けされた温度分布が表示可能となる。 The control calculation unit 6 executes the above-described calculation process to visualize the calculated temperature value T * realized on the display device 8. Thereby, for example, the color-coded temperature distribution as shown in FIG. 9C can be displayed.
次に、制御演算部6が実行する具体的処理の流れについて説明する。
図10は、制御演算部6により実行されるマイクロ断層可視化処理の流れを示すフローチャートである。なお、この断層可視化処理は、上述したOCTによる断層撮影がなされることを前提に実行される。ここでは、測定対象Wが熱変形したときの変形前後の断層画像を取得し、物理量として測定対象Wの内部の温度分布を可視化する場合を例に説明する。
Next, the flow of specific processing executed by the control calculation unit 6 will be described.
FIG. 10 is a flowchart showing the flow of the micro tomography visualization process executed by the control calculator 6. It should be noted that this tomographic visualization processing is executed on the premise that the tomographic imaging by OCT described above is performed. Here, a case will be described as an example in which a tomographic image before and after the deformation of the measurement target W is thermally deformed and the temperature distribution inside the measurement target W is visualized as a physical quantity.
制御演算部6は、OCTにより撮影された前後2枚のOCT断層画像を読み込む(S10)。続いて、それらの断層画像に基づいて変形ベクトル分布を算出するための変形ベクトル演算処理を実行する(S12)。続いて、その変形ベクトル分布を用いて逆解析処理を実行する(S14)。そして、その逆解析により算出された物理量の断層分布(温度分布)を表示装置8に表示させる(S16)。 The control calculation unit 6 reads two front and rear OCT tomographic images captured by OCT (S10). Subsequently, a deformation vector calculation process for calculating the deformation vector distribution based on those tomographic images is executed (S12). Then, an inverse analysis process is executed using the deformation vector distribution (S14). Then, the fault distribution (temperature distribution) of the physical quantity calculated by the inverse analysis is displayed on the display device 8 (S16).
図11は、図10におけるS12の変形ベクトル演算処理を詳細に示すフローチャートである。制御演算部6は、変形前後の断層画像に基づき、再帰的相互相関法による処理を実行する。ここではまず、最小解像度(最大サイズの検査領域)での相互相関処理を実行し、相関係数分布を求める(S20)。続いて、隣接相互相関乗法により、隣接する相関係数分布の積を演算する(S22)。このとき、標準偏差フィルタ等の空間フィルタにより過誤ベクトルの除去をし(S24)、最小二乗法等により除去ベクトルの内挿補間を実行する(S26)。続いて、検査領域を小さくすることによって解像度を上げて相互相関処理を継続する(S28)。すなわち、低解像度での参照ベクトルを基に相互相関処理を実行する。このときの解像度が予め定める最高解像度でなければ(S30のN)、S22に戻る。 FIG. 11 is a flowchart showing in detail the deformation vector calculation process of S12 in FIG. The control calculation unit 6 executes processing by the recursive cross-correlation method based on the tomographic images before and after the deformation. Here, first, cross-correlation processing is performed at the minimum resolution (inspection area having the maximum size) to obtain the correlation coefficient distribution (S20). Then, the product of the adjacent correlation coefficient distributions is calculated by the adjacent cross correlation multiplication method (S22). At this time, the error vector is removed by a spatial filter such as a standard deviation filter (S24), and the removal vector is interpolated by the least square method or the like (S26). Subsequently, the inspection area is reduced to increase the resolution and continue the cross-correlation processing (S28). That is, the cross-correlation processing is executed based on the reference vector with low resolution. If the resolution at this time is not the predetermined maximum resolution (N in S30), the process returns to S22.
そして、S22〜S28の処理を繰り返し、最高解像度での相互相関処理が完了すると(S30のY)、サブピクセル解析を実行する。すなわち、最高解像度(最小サイズの検査領域)での変形ベクトルの分布に基づき、風上勾配法によるサブピクセル移動量を演算する(S32)。そして、このとき算出されたサブピクセル移動量に基づき、画像変形法によるサブピクセル変形量を演算する(S34)。続いて、最大相互相関値によるフィルタ処理により過誤ベクトルの除去をし(S36)、最小二乗法等により除去ベクトルの内挿補間を実行する(S38)。このようにして変形ベクトル分布が得られる。 Then, when the processes of S22 to S28 are repeated and the cross-correlation process at the highest resolution is completed (Y of S30), the sub-pixel analysis is executed. That is, the sub-pixel movement amount by the upwind gradient method is calculated based on the distribution of the deformation vector at the highest resolution (minimum size inspection area) (S32). Then, the sub-pixel deformation amount by the image deformation method is calculated based on the sub-pixel movement amount calculated at this time (S34). Subsequently, the error vector is removed by the filtering process using the maximum cross-correlation value (S36), and the removal vector is interpolated by the least square method or the like (S38). In this way, the deformation vector distribution is obtained.
図12は、図10におけるS14の逆解析処理を詳細に示すフローチャートである。制御演算部6は、演算対象となった断層画像に対応するROIを設定する(S40:図7(A)参照)。この設定は、当該解析を行うユーザの入力指示にしたがって行われる。ユーザは、測定対象Wにおいて温度分布を表示させたい箇所をROIとして設定することができ、入力装置7を介してそれを指示することができる。 FIG. 12 is a flowchart showing in detail the inverse analysis process of S14 in FIG. The control calculation unit 6 sets the ROI corresponding to the tomographic image that is the calculation target (S40: see FIG. 7A). This setting is performed according to the input instruction of the user who performs the analysis. The user can set a portion of the measurement object W where the temperature distribution is to be displayed as the ROI, and can instruct it via the input device 7.
続いて、制御演算部6は、設定されたROIに対して要素分割を設定する(S42:図7(B)参照)。この要素分割は、上述したFEMに対応するものであり、ユーザの入力指示にしたがって行われる。ユーザは、求める解析精度に応じた要素分割を入力装置7を介して指示することができる。 Then, the control calculation part 6 sets element division with respect to the set ROI (S42: refer FIG. 7(B)). This element division corresponds to the FEM described above, and is performed according to the input instruction from the user. The user can instruct the element division according to the desired analysis accuracy via the input device 7.
続いて、制御演算部6は、S38を経て得られた変形ベクトル分布を、分割された要素の節点に補間する(S44:図8参照)。また、ROIに対して弾性率および線膨張係数を設定する(S46)。本実施例では、これらの値を測定対象Wの材料に固有の値として設定する。 Subsequently, the control calculation unit 6 interpolates the deformation vector distribution obtained through S38 into the nodes of the divided elements (S44: see FIG. 8). Further, an elastic modulus and a linear expansion coefficient are set for the ROI (S46). In this embodiment, these values are set as values unique to the material of the measurement object W.
続いて、制御演算部6は、仮想仕事の原理に基づいて上記式(16)の全体剛性マトリクス[K]を演算する一方(S48)、形状関数,力学パラメータ(弾性率など),温度依存パラメータ(線膨張係数)から温度変化項の全体マトリクス観測行列[A]を決定する(S50)。そして、上述したマトリクスの縮小などを行い、逆解析(例えば逆行列演算)により温度分布を演算する(S52)。以上のようにして算出された温度分布が、図10のS16にて断層可視化される。 Subsequently, the control calculation unit 6 calculates the overall stiffness matrix [K] of the above formula (16) based on the principle of virtual work (S48), while the shape function, the mechanical parameter (elastic modulus etc.), the temperature dependent parameter are calculated. The entire matrix observation matrix [A] of the temperature change term is determined from (linear expansion coefficient) (S50). Then, the matrix is reduced as described above, and the temperature distribution is calculated by inverse analysis (for example, inverse matrix calculation) (S52). The temperature distribution calculated as described above is visualized as a slice in S16 of FIG.
次に、本実施例による効果について説明する。発明者は、本実施例による手法、つまりOCTによる測定とFEMの逆解析とによるハイブリッド手法の有効性を確認するために、数値実験による精度検証を行った。以下、その検証結果について説明する。図13〜図15は、数値実験による検証方法を示す図である。各図の(A)〜(C)はその検証過程を示す。図16は、検証結果を示す図である。(A)は本実施例のハイブリッド手法による解析結果を示し、(B)はFEMによる正解値を示す。図17は、本実施例の解析による正解値との誤差を示す図である。(A)は絶対誤差を示し、(B)は相対誤差を示す。 Next, the effect of this embodiment will be described. The inventor conducted accuracy verification by numerical experiments in order to confirm the effectiveness of the method according to the present embodiment, that is, the hybrid method based on the measurement by OCT and the inverse analysis of FEM. The verification result will be described below. 13 to 15 are diagrams showing a verification method by a numerical experiment. (A) to (C) of each figure show the verification process. FIG. 16 is a diagram showing a verification result. (A) shows the analysis result by the hybrid method of the present embodiment, and (B) shows the correct value by FEM. FIG. 17 is a diagram showing an error from the correct answer value by the analysis of the present embodiment. (A) shows an absolute error, (B) shows a relative error.
上記ハイブリッド手法の精度を評価するためには、その比較対象となる正解値が必要となる。本来であれば、この比較対象を温度センサ等による実測値(真値)として得るのが望ましいが、物体内部の断層に沿ってそれを行うのは不可能に等しい。そこで、ここでは測定対象のモデルに対してFEMによる順解析を実行し、それにより得られる温度分布を正解値とする。そして、その正解値からOCT断層画像を疑似的に作成し、その断層画像(以下「疑似OCT断層画像」ともいう)をOCTによる実際の断層画像と仮定して上記逆解析を実行する。その逆解析によって得られた温度分布と、正解値としての温度分布(以下「正解温度分布」ともいう)とを比較することで、上記ハイブリッド手法の精度を評価した。 In order to evaluate the accuracy of the hybrid method, the correct answer value to be compared is required. Originally, it is desirable to obtain this comparison target as a measured value (true value) by a temperature sensor or the like, but it is almost impossible to do it along the fault inside the object. Therefore, here, the forward analysis by FEM is executed on the model to be measured, and the temperature distribution obtained by this is taken as the correct value. Then, an OCT tomographic image is pseudo-created from the correct answer value, and the inverse analysis is executed assuming that the tomographic image (hereinafter also referred to as “pseudo-OCT tomographic image”) is an actual tomographic image by OCT. The accuracy of the hybrid method was evaluated by comparing the temperature distribution obtained by the inverse analysis with the temperature distribution as the correct value (hereinafter also referred to as “correct temperature distribution”).
(数値実験)
具体的には、モンテカルロシミュレーションによって散乱体からのランダムスッペクルパターンを演算することで、図14(B)に示す疑似OCT断層画像を作成し、FEMを用いて散乱体の熱変形移動量の算出を行うこととした。
(Numerical experiment)
Specifically, a pseudo-OCT tomographic image shown in FIG. 14B is created by calculating a random speckle pattern from a scatterer by Monte Carlo simulation, and the thermal deformation movement amount of the scatterer is calculated using FEM. Decided to do.
まず、数値モデルに基づいたFEMによる伝熱解析を行った。図13(A)に示すように、測定対象として基板にICチップ(シリコンチップ)が埋め込まれた電子実装部品を模擬したモデルを想定している。このモデルは、エポキシ樹脂(Epoxy)を基材とし、シリコンチップのシリコン(Silicon)、配線の銅(Copper)を内包し、表面にはレジスト材(Resist)が存在する複合材料とされている。 First, heat transfer analysis by FEM based on a numerical model was performed. As shown in FIG. 13A, it is assumed that a model that simulates an electronic mounting component in which an IC chip (silicon chip) is embedded in a substrate as a measurement target. This model is a composite material in which epoxy resin (Epoxy) is used as a base material, silicon (Silicon) of a silicon chip and copper (Copper) of wiring are included, and a resist material (Resist) is present on the surface.
そして、図13(B)に示すように、このモデルに対してFEMの順解析による要素分割を施した。なお、モデルを構成する各材料の材料定数は、下記の表1に示すとおりである。
温度の境界条件については、シリコンチップおよび配線が発熱していると仮定した。すなわち、シリコンおよび銅の温度を250[℃]とし、表面(上下面)であるz=0,2000[μm]の温度を外気温に等しい25[℃]とした。また、熱伝達係数を10[W/(mK)]とした。ただし、モデルの側面(x=0,5000[μm])では熱の出入りはないと仮定した。このような設定のもとFEMの伝熱解析を実行することで、図13(C)に示す温度分布を得た。この温度分布が数値実験における「正解温度分布」とされる。 Regarding the boundary condition of temperature, it is assumed that the silicon chip and the wiring are generating heat. That is, the temperature of silicon and copper was set to 250 [° C.], and the temperature of the surface (upper and lower surfaces) of z=0,2000 [μm] was set to 25 [° C.] equal to the outside air temperature. Further, the heat transfer coefficient was set to 10 [W/(mK)]. However, it was assumed that heat did not flow in and out on the side surface of the model (x=0,5000 [μm]). The temperature distribution shown in FIG. 13(C) was obtained by executing the FEM heat transfer analysis under such settings. This temperature distribution is called "correct temperature distribution" in the numerical experiment.
この温度分布データを用いてFEMによる力学解析を実行し、モデルを熱変形させることとした。その際、モデルの境界拘束条件として、z=0[μm]の位置を全方向に固定することとした。ICチップを内包する電子実装部品は、マザー基板へ実装される際にその底面がはんだ付けされることを考慮したものである。 It was decided to carry out a mechanical analysis by FEM using this temperature distribution data to thermally deform the model. At that time, as a boundary constraint condition of the model, the position of z=0 [μm] was fixed in all directions. The electronic mounting component including the IC chip is designed considering that the bottom surface of the electronic mounting component is soldered when mounted on the mother board.
図14(A)には、このような連成解析により算出された熱変形ベクトル分布が示されている。この熱変形量データを図14(B)の各スペックルに対応させ、スペックルの移動を決定した。各スペックルの熱変形移動は、スペックルを内包する有限要素の節点変形ベクトルおよび形状関数を用いて算出することができる。それにより、図14(C)に示す熱変形後の疑似OCT断層画像を作成した。 FIG. 14A shows the thermal deformation vector distribution calculated by such a coupled analysis. The heat deformation amount data was made to correspond to each speckle in FIG. 14B, and the movement of the speckle was determined. The thermal deformation movement of each speckle can be calculated using a nodal deformation vector and a shape function of a finite element including the speckle. Thereby, the pseudo OCT tomographic image after thermal deformation shown in FIG. 14C was created.
このようにして、図14(B)および(C)に示すように、熱変形前後の疑似OCT断層画像を得た。なお、ここでいう「熱変形前」とは、シリコンチップおよび配線が発熱する前の状態、つまり両者の温度が実質的に外気温である25[℃]に等しい状態を意味する。「熱変形後」とは、シリコンチップおよび配線が上述のように250[℃]に発熱した状態を意味する。熱変形前後の疑似OCT断層画像は、それぞれの状態を仮定して得たものである。この2つの疑似OCT断層画像の取得を出発点として上記ハイブリッド手法を適用した。すなわち、変形前後の疑似OCT断層画像をOCTによる実際の断層画像と仮定し、図10〜図12に示した処理を実行した。そして、図10のS16にて得られた温度分布を、図13(C)に示した正解温度分布と比較することで精度を評価した。 In this way, pseudo OCT tomographic images before and after thermal deformation were obtained as shown in FIGS. The term "before thermal deformation" as used herein means a state before the silicon chip and the wiring generate heat, that is, a state in which the temperature of both is substantially equal to the ambient temperature of 25[°C]. “After thermal deformation” means a state in which the silicon chip and the wiring generate heat at 250[° C.] as described above. The pseudo-OCT tomographic images before and after thermal deformation are obtained assuming each state. The above hybrid method was applied with the acquisition of these two pseudo OCT tomographic images as a starting point. That is, assuming that the pseudo OCT tomographic images before and after the deformation are actual tomographic images by OCT, the processes shown in FIGS. 10 to 12 were executed. Then, the accuracy was evaluated by comparing the temperature distribution obtained in S16 of FIG. 10 with the correct temperature distribution shown in FIG.
具体的には、図15(A)に示す疑似OCT断層画像に対し、図15(B)に示すようにROIを設定し、そのROIについて要素分割を行った。なお、このROI内の弾性率分布,線膨張係数分布を定義する必要があるが、弾性率および線膨張係数は材料に固有の値であるため、ここでは上記表1に示す材料定数を設定した。 Specifically, with respect to the pseudo OCT tomographic image shown in FIG. 15(A), an ROI is set as shown in FIG. 15(B), and the ROI is divided into elements. Although it is necessary to define the elastic modulus distribution and the linear expansion coefficient distribution within this ROI, since the elastic modulus and the linear expansion coefficient are values unique to the material, the material constants shown in Table 1 above are set here. ..
図15(C)には、前後の疑似OCT断層画像に基づいて算出された変形ベクトル分布が示されている。この変形ベクトル分布に基づいて上述した逆解析処理を実行し、図16(A)に示すように温度分布を断層可視化した。この温度分布を図16(B)に示す正解温度分布と比較すると、それらが定性的に一致していることが分かる。すなわち、いずれも基材中心部分の温度がその周囲より低くなっているのに対し、シリコンチップおよび銅の付近で周囲より高い温度分布が得られていることが確認できる。 In FIG. 15C, the deformation vector distribution calculated based on the front and rear pseudo OCT tomographic images is shown. The above-mentioned inverse analysis processing was executed based on this deformation vector distribution, and the temperature distribution was visualized as a cross section as shown in FIG. 16(A). When this temperature distribution is compared with the correct temperature distribution shown in FIG. 16B, it is found that they are qualitatively consistent. That is, it can be confirmed that the temperature of the central portion of the base material is lower than that of its surroundings, while the temperature distribution near the silicon chip and copper is higher than that of the surroundings.
図17(A)には、本実施例のハイブリッド手法により得られた温度分布(「推定温度分布」ともいう)と、正解温度分布との絶対誤差が示されている。図17(B)には両者の相対誤差が示されている。なお、「絶対誤差」とは両者の温度差を意味し、「相対誤差」とは絶対誤差を正解値により除算した値を意味する。 FIG. 17A shows the absolute error between the temperature distribution (also referred to as “estimated temperature distribution”) obtained by the hybrid method of this embodiment and the correct temperature distribution. FIG. 17(B) shows the relative error between them. The "absolute error" means the temperature difference between the two, and the "relative error" means a value obtained by dividing the absolute error by the correct value.
これらの結果に基づき、両者の誤差平均(RMS誤差)は14%程度と比較的小さく収まることが確認できた。なお、ROIの境界近傍でエラーが大きくなる傾向にあるが、これは境界でのノイズの影響によるものと考えられる。なお、このような温度推定の誤差は、温度昇降を動的に検出してノイズ低減を施す方法や、ノイズを考慮した逆推定手法であるカルマンフィルタなどの導入により改善可能である。 Based on these results, it was confirmed that the average error (RMS error) between them was relatively small, about 14%. The error tends to increase near the ROI boundary, which is considered to be due to the effect of noise at the boundary. Note that such an error in temperature estimation can be improved by introducing a method of dynamically detecting temperature rise and fall to reduce noise, a Kalman filter that is an inverse estimation method considering noise, or the like.
以上に説明したように、本実施例によれば、OCTによる断層計測(実計測データ)を用いることにより、材料内部の温度断層分布を逆解析することが可能となる。すなわち、実測と解析とのハイブリッドにより物理量を断層可視化できる。物体内部における任意の領域の温度分布を、境界条件を考慮することなく断層可視化できることから、特に境界条件が不明確な状況における内部温度の検出に好適となる。また、このような温度断層分布がマイクロスケールで非侵襲・非破壊にて検出可能であるため、MEMS等の電子実装部品のように境界条件の付与が難しい工業製品への応用が期待される。 As described above, according to the present embodiment, it is possible to inversely analyze the temperature tomographic distribution inside the material by using the tomographic measurement (actual measurement data) by OCT. That is, a physical quantity can be visualized by a hybrid of actual measurement and analysis. Since it is possible to visualize the temperature distribution of an arbitrary region inside the object without considering the boundary conditions, it is suitable for detecting the internal temperature particularly in a situation where the boundary conditions are unclear. Further, since such a temperature fault distribution can be detected on a microscale in a non-invasive and non-destructive manner, it is expected to be applied to an industrial product such as an electronic mounting component such as MEMS in which boundary conditions are difficult to apply.
以上、本発明の好適な実施例について説明したが、本発明はその特定の実施例に限定されるものではなく、本発明の技術思想の範囲内で種々の変形が可能であることはいうまでもない。 The preferred embodiments of the present invention have been described above, but it goes without saying that the present invention is not limited to the specific embodiments and various modifications can be made within the scope of the technical idea of the present invention. Nor.
上記実施例では、OCT断層画像として、構造の変形前後の静止画を用いる例を示した。変形例においては、構造の変形中を撮影したOCT断層画像の動画に対して上記手法を適用してもよい。すなわち、動画中の異なる時刻の2枚のOCT断層画像を取得し、それらに対して上記逆解析を適用して物理量を演算し、断層可視化してもよい。上記ハイブリッド手法は、このように時間変化する現象に対しても適用可能である。 In the above-described embodiment, an example in which still images before and after structural deformation are used as the OCT tomographic image has been shown. In a modified example, the above method may be applied to a moving image of an OCT tomographic image of the structure being deformed. That is, two OCT tomographic images at different times in the moving image may be acquired, and the inverse analysis may be applied to them to calculate a physical quantity and visualize the tomographic image. The hybrid method can be applied to such a time-varying phenomenon.
上記実施例では述べなかったが、測定対象の温度場が変化しているような非定常な状況では、変形している場(分布)が時々刻々と変化する。そのような状況下においてOCT断層画像を連続撮影して本手法を適用することにより、時々刻々変化する非定常温度場(分布)をマイクロ断層可視化することが可能となる。この演算過程において、連続するOCT断層画像と時空間最小二乗法を用いた節点補間を行うなどが考えられる。 Although not described in the above embodiment, in an unsteady situation where the temperature field of the measurement target is changing, the deforming field (distribution) changes momentarily. Under such a circumstance, by continuously photographing OCT tomographic images and applying the present method, it becomes possible to visualize a micro-tomographical view of an unsteady temperature field (distribution) that changes from moment to moment. In this calculation process, continuous OCT tomographic images and nodal interpolation using the spatiotemporal least squares method may be performed.
上記実施例では、変形ベクトルの断層分布に基づいて逆解析を行い、それにより温度分布を算出する例を示した。変形例においては、変形ベクトルの時間微分である変形速度ベクトルの断層分布に基づいて逆解析を行い、それにより物理量の分布を算出してもよい。 In the above-described embodiment, an example has been shown in which the inverse analysis is performed based on the fault distribution of the deformation vector, and the temperature distribution is calculated accordingly. In the modified example, an inverse analysis may be performed based on the tomographic distribution of the deformation velocity vector, which is the time derivative of the deformation vector, and the distribution of the physical quantity may be calculated accordingly.
上記実施例では、OCTによる断層計測と数値解析の逆解析とにより部品内部の温度断層分布を可視化する例を示した。変形例においては、同様の手法により部品内部の弾性係数や粘弾性係数の分布を可視表示することも可能である。例えば、上記式(20)に関して温度変化{ΔT}がゼロであるなど既知の場合、測定項{U}の取得により全体剛性マトリクス[K]を算出し、弾性係数の分布を可視表示することができる。それにより、部品内部に存在する材料を特定することも可能となる。 In the above-described embodiment, the example in which the temperature fault distribution inside the component is visualized by the fault measurement by OCT and the inverse analysis of the numerical analysis is shown. In the modified example, the distribution of the elastic coefficient and the viscoelastic coefficient inside the component can be visually displayed by the same method. For example, when it is known that the temperature change {ΔT} is zero with respect to the above equation (20), the overall stiffness matrix [K] can be calculated by acquiring the measurement term {U}, and the elastic modulus distribution can be visually displayed. it can. This also makes it possible to identify the material that exists inside the component.
上記実施例では、OCTによる断層画像から得た変形ベクトル分布を逆解析することにより物理量を算出する例を示したが、OCTに代えてX線CT、超音波、MRIなどを採用することもできる。例えば、電子実装部品などの高分子基複合材や金属との複合材に適用する際にはマイクロX線CTによる断層画像などが十分適用できると考えられる。 In the above-mentioned embodiment, the example in which the physical quantity is calculated by inversely analyzing the deformation vector distribution obtained from the tomographic image by OCT has been shown, but X-ray CT, ultrasonic waves, MRI or the like may be adopted instead of OCT. .. For example, when applied to a polymer-based composite material such as an electronic mounting component or a composite material with a metal, it is considered that a tomographic image by micro X-ray CT can be sufficiently applied.
上記実施例では、OCTによる断層画像を二次元で取得する例を示したが、三次元で取得してもよい。すなわち、変位関連ベクトルとしての変形ベクトルあるいは変形速度ベクトルを三次元データとして処理してもよいことは言うまでもない。 In the above-described embodiment, the example in which the tomographic image by OCT is acquired in two dimensions is shown, but it may be acquired in three dimensions. That is, it goes without saying that the deformation vector or the deformation velocity vector as the displacement-related vector may be processed as three-dimensional data.
上記実施例では、物理量を断層可視化する対象として部品構造を例示した。変形例においては、その他の材料構造を対象としてもよい。また、体内の軟骨やしこり(硬化箇所)など、生体構造(生体組織)を対象としてもよい。生体構造の力学特性の変化を断層可視化することで、例えば動脈硬化,癌,肝硬変,皮膚,軟骨,骨などの多くの症例に対する診断材料とすることが可能となる。 In the above-described embodiment, the component structure is exemplified as the target of which the physical quantity is tomographically visualized. Other material structures may be targeted in the modification. In addition, biological structures (living tissues) such as cartilage and lumps (hardened areas) in the body may be targeted. By visualizing the changes in the mechanical properties of the anatomy by tomography, it becomes possible to use it as a diagnostic material for many cases such as arteriosclerosis, cancer, liver cirrhosis, skin, cartilage, and bone.
なお、本発明は上記実施例や変形例に限定されるものではなく、要旨を逸脱しない範囲で構成要素を変形して具体化することができる。上記実施例や変形例に開示されている複数の構成要素を適宜組み合わせることにより種々の発明を形成してもよい。また、上記実施例や変形例に示される全構成要素からいくつかの構成要素を削除してもよい。 It should be noted that the present invention is not limited to the above-described embodiments and modifications, and the constituent elements can be modified and embodied without departing from the scope of the invention. Various inventions may be formed by appropriately combining a plurality of constituent elements disclosed in the above-described embodiments and modifications. Further, some constituent elements may be deleted from all the constituent elements shown in the above-described embodiments and modifications.
1 断層可視化システム、2 光学ユニット、4 光学機構、6 制御演算部、7 入力装置、8 表示装置、10 光源、12 オブジェクトアーム、14 リファレンスアーム、16 光検出器、26 反射鏡、ROI 関心領域、S1 検査領域、S2 探査領域。 1 tomographic visualization system, 2 optical unit, 4 optical mechanism, 6 control operation part, 7 input device, 8 display device, 10 light source, 12 object arm, 14 reference arm, 16 photodetector, 26 reflector, ROI region of interest, S1 inspection area, S2 exploration area.
Claims (9)
前記構造の変形前後の断層画像を取得する工程と、
前記断層画像に基づき変位関連ベクトルの断層分布を演算する工程と、
前記断層画像に対して関心領域を設定する工程と、
前記関心領域を要素分割する工程と、
演算された変位関連ベクトルを、分割された要素の節点に補間するとともに境界条件として与える工程と、
補間後の変位関連ベクトルの断層分布を用いて前記数値解析の逆解析を境界条件の設定を別途行わずに実行することにより、各要素における前記物理量の分布を前記物理量の断層分布として演算する工程と、
演算された前記物理量の分布を表示装置に断層可視化する工程と、
を備えることを特徴とするマイクロ断層可視化方法。 Regarding a physical quantity that can be an analysis condition in a predetermined numerical analysis in which a displacement-related vector is calculated in the process of a forward analysis targeting a structure, a micro-tomographic visualization method for visualizing the distribution of the physical quantity inside the structure by tomography,
Acquiring tomographic images before and after deformation of the structure,
Calculating a tomographic distribution of displacement-related vectors based on the tomographic image,
Setting a region of interest for the tomographic image,
Dividing the region of interest into elements,
Interpolating the calculated displacement-related vector to the nodes of the divided elements and giving it as a boundary condition ,
A step of calculating a distribution of the physical quantity in each element as a fault distribution of the physical quantity by performing an inverse analysis of the numerical analysis using the fault distribution of the displacement-related vector after interpolation without separately setting boundary conditions. When,
A step of visualizing the computed distribution of the physical quantity on a display device;
A method for visualizing micro-tomography, comprising:
前記数値解析が有限要素法による熱変形を考慮した熱応力解析であることを特徴とする請求項3に記載のマイクロ断層可視化方法。 The physical quantity is temperature,
The method for visualizing micro-tomography according to claim 3, wherein the numerical analysis is a thermal stress analysis considering thermal deformation by a finite element method.
前記数値解析が有限要素法による応力解析であることを特徴とする請求項3に記載のマイクロ断層可視化方法。 The physical quantity is a physical constant of the structure,
The method for visualizing a micro tomography according to claim 3, wherein the numerical analysis is a stress analysis by a finite element method.
前記逆解析により得られた物理量のうち、前記領域に対応する物理量である推定値と、前記領域において与えられた真値との関係に基づいて、前記逆解析により得られた物理量の値を補正する工程と、
をさらに備えることを特徴とする請求項1〜7のいずれかに記載のマイクロ断層可視化方法。 Providing a true value of the physical quantity in a predetermined region of the structure,
Of the physical quantity obtained by the inverse analysis, the value of the physical quantity obtained by the inverse analysis is corrected based on the relationship between the estimated value, which is the physical quantity corresponding to the area, and the true value given in the area. The process of
The method for visualizing micro tomography according to any one of claims 1 to 7, further comprising:
前記構造の変形前後の断層画像を取得するための断層画像取得装置と、
前記断層画像取得装置から前記構造の変形前後の断層画像データを取得し、その断層画像データに基づいて前記物理量の断層分布を演算する制御演算部と、
前記制御演算部の演算結果に基づいて、前記構造における前記物理量の分布を断層可視化する態様で表示する表示装置と、
を備え、
前記物理量は、前記構造を対象とする順解析の過程で変位関連ベクトルが演算される所定の数値解析において解析条件となり得るものであり、
前記制御演算部は、
取得した断層画像データに基づき変位関連ベクトルの断層分布を演算し、
前記断層画像に対して関心領域を設定し、
前記関心領域を要素分割し、
演算された変位関連ベクトルを、分割された要素の節点に補間するとともに境界条件として与え、
補間後の変位関連ベクトルの断層分布を用いて前記数値解析の逆解析を境界条件の設定を別途行わずに実行することにより、各要素における前記物理量の分布を前記物理量の断層分布として演算し、
演算された前記物理量の分布を前記表示装置に断層可視化させることを特徴とするマイクロ断層可視化システム。 A micro-tomography system that visualizes the distribution of a given physical quantity inside the structure to be measured,
A tomographic image acquisition device for acquiring a tomographic image before and after the deformation of the structure,
A control calculation unit that acquires tomographic image data before and after the deformation of the structure from the tomographic image acquisition device, and calculates a tomographic distribution of the physical quantity based on the tomographic image data,
A display device that displays the distribution of the physical quantity in the structure in a form of a tomographic visualization based on the calculation result of the control calculation unit;
Equipped with
The physical quantity can be an analysis condition in a predetermined numerical analysis in which a displacement-related vector is calculated in the process of a forward analysis targeting the structure,
The control calculation unit,
Calculate the tomographic distribution of displacement-related vectors based on the acquired tomographic image data,
Set a region of interest for the tomographic image,
Dividing the region of interest into elements,
The calculated displacement related vector is interpolated into the nodes of the divided elements and given as a boundary condition,
By performing the reverse analysis of the numerical analysis using the fault distribution of the displacement-related vector after interpolation without separately setting the boundary condition, the distribution of the physical quantity in each element is calculated as the fault distribution of the physical quantity,
A micro-tomography visualization system characterized by causing the display device to visualize the distribution of the calculated physical quantity.
Priority Applications (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
JP2015050553A JP6707249B2 (en) | 2015-03-13 | 2015-03-13 | Micro-tomography visualization method and system |
Applications Claiming Priority (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
JP2015050553A JP6707249B2 (en) | 2015-03-13 | 2015-03-13 | Micro-tomography visualization method and system |
Publications (2)
Publication Number | Publication Date |
---|---|
JP2016170080A JP2016170080A (en) | 2016-09-23 |
JP6707249B2 true JP6707249B2 (en) | 2020-06-10 |
Family
ID=56983546
Family Applications (1)
Application Number | Title | Priority Date | Filing Date |
---|---|---|---|
JP2015050553A Active JP6707249B2 (en) | 2015-03-13 | 2015-03-13 | Micro-tomography visualization method and system |
Country Status (1)
Country | Link |
---|---|
JP (1) | JP6707249B2 (en) |
Families Citing this family (3)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
JP6498796B1 (en) * | 2018-01-11 | 2019-04-10 | エムテックスマツムラ株式会社 | Semiconductor device inspection method, semiconductor device manufacturing method, semiconductor device inspection device |
JP7423180B2 (en) * | 2018-06-26 | 2024-01-29 | 公益財団法人鉄道総合技術研究所 | High-precision position correction method and system for waveform data |
KR102303643B1 (en) | 2021-06-04 | 2021-09-17 | 주식회사 테크디알 | Inspection device and method for detecting air void at the lead film of battery |
Family Cites Families (6)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
JP4258015B2 (en) * | 2002-07-31 | 2009-04-30 | 毅 椎名 | Ultrasonic diagnostic system, strain distribution display method, and elastic modulus distribution display method |
JP4221555B2 (en) * | 2002-07-31 | 2009-02-12 | 毅 椎名 | Ultrasonic diagnostic system, strain distribution display method, and elastic modulus distribution display method |
ES2379468T3 (en) * | 2004-08-24 | 2012-04-26 | The General Hospital Corporation | Procedure, system and software configuration to determine the modulus of elasticity |
JP4257982B2 (en) * | 2006-03-14 | 2009-04-30 | 国立大学法人山口大学 | Strain distribution measurement system, elastic modulus distribution measurement system and methods thereof |
JP4389091B2 (en) * | 2008-01-29 | 2009-12-24 | 毅 椎名 | Ultrasonic diagnostic system, strain distribution display method, and elastic modulus distribution display method |
FR3010218B1 (en) * | 2013-09-03 | 2016-12-30 | Univ Joseph Fourier - Grenoble 1 | METHOD OF IMAGE PROCESSING BASED ON THE TECHNIQUE OF FINISHED ELEMENTS FOR THE DIRECT RESOLUTION OF REVERSE PROBLEMS IN MECHANICAL STRUCTURES. |
-
2015
- 2015-03-13 JP JP2015050553A patent/JP6707249B2/en active Active
Also Published As
Publication number | Publication date |
---|---|
JP2016170080A (en) | 2016-09-23 |
Similar Documents
Publication | Publication Date | Title |
---|---|---|
Reu et al. | DIC challenge: developing images and guidelines for evaluating accuracy and resolution of 2D analyses | |
JP6707258B2 (en) | Stress visualization system and mechanical property visualization system | |
Yoneyama | Basic principle of digital image correlation for in-plane displacement and strain measurement | |
Romeo | Two-dimensional digital image correlation for asphalt mixture characterisation: interest and limitations | |
Thompson et al. | Tissue perfusion measurements: multiple-exposure laser speckle analysis generates laser Doppler–like spectra | |
WO2015129336A1 (en) | Ultrasonic image pickup device and method | |
US20110262015A1 (en) | Image processing apparatus, image processing method, and storage medium | |
Sun et al. | Digital image correlation–based optical coherence elastography | |
Mudassar et al. | Improved digital image correlation method | |
Swillens et al. | Comparison of non-invasive methods for measurement of local pulse wave velocity using FSI-simulations and in vivo data | |
Dai et al. | Strain field estimation based on digital image correlation and radial basis function | |
JP7051111B2 (en) | Estimating electromechanical quantities using filter technology based on digital images and models | |
Wijesinghe et al. | Computational optical palpation: a finite-element approach to micro-scale tactile imaging using a compliant sensor | |
JP6707249B2 (en) | Micro-tomography visualization method and system | |
Clements et al. | Evaluation of model-based deformation correction in image-guided liver surgery via tracked intraoperative ultrasound | |
Chen et al. | Noninvasive, three-dimensional full-field body sensor for surface deformation monitoring of human body in vivo | |
WO2018212115A1 (en) | Device and method for tomographic visualization of tissue viscoelasticity | |
Mcgarry et al. | An inverse approach to determining spatially varying arterial compliance using ultrasound imaging | |
Xu et al. | Numerical simulations of flow patterns in the human left ventricle model with a novel dynamic mesh morphing approach based on radial basis function | |
Riemer et al. | Determining haemodynamic wall shear stress in the rabbit aorta in vivo using contrast-enhanced ultrasound image velocimetry | |
Tan et al. | Influence of scan resolution, thresholding, and reconstruction algorithm on computed tomography-based kinematic measurements | |
Claus et al. | Large-field-of-view optical elastography using digital image correlation for biological soft tissue investigation | |
Campo et al. | Digital image correlation for full-field time-resolved assessment of arterial stiffness | |
Solav et al. | Duodic: 3d digital image correlation in Matlab | |
Kleinendorst et al. | Mechanical Shape Correlation: A novel integrated digital image correlation approach |
Legal Events
Date | Code | Title | Description |
---|---|---|---|
RD02 | Notification of acceptance of power of attorney |
Free format text: JAPANESE INTERMEDIATE CODE: A7422 Effective date: 20150406 |
|
A521 | Request for written amendment filed |
Free format text: JAPANESE INTERMEDIATE CODE: A821 Effective date: 20150407 |
|
RD02 | Notification of acceptance of power of attorney |
Free format text: JAPANESE INTERMEDIATE CODE: A7422 Effective date: 20160613 |
|
A521 | Request for written amendment filed |
Free format text: JAPANESE INTERMEDIATE CODE: A821 Effective date: 20160614 |
|
A621 | Written request for application examination |
Free format text: JAPANESE INTERMEDIATE CODE: A621 Effective date: 20180213 |
|
A977 | Report on retrieval |
Free format text: JAPANESE INTERMEDIATE CODE: A971007 Effective date: 20181213 |
|
A131 | Notification of reasons for refusal |
Free format text: JAPANESE INTERMEDIATE CODE: A131 Effective date: 20190108 |
|
A601 | Written request for extension of time |
Free format text: JAPANESE INTERMEDIATE CODE: A601 Effective date: 20190218 |
|
A521 | Request for written amendment filed |
Free format text: JAPANESE INTERMEDIATE CODE: A523 Effective date: 20190507 |
|
A711 | Notification of change in applicant |
Free format text: JAPANESE INTERMEDIATE CODE: A712 Effective date: 20190724 |
|
A02 | Decision of refusal |
Free format text: JAPANESE INTERMEDIATE CODE: A02 Effective date: 20191001 |
|
A521 | Request for written amendment filed |
Free format text: JAPANESE INTERMEDIATE CODE: A523 Effective date: 20191226 |
|
A711 | Notification of change in applicant |
Free format text: JAPANESE INTERMEDIATE CODE: A711 Effective date: 20191226 |
|
A521 | Request for written amendment filed |
Free format text: JAPANESE INTERMEDIATE CODE: A821 Effective date: 20191226 |
|
A521 | Request for written amendment filed |
Free format text: JAPANESE INTERMEDIATE CODE: A523 Effective date: 20200220 |
|
A911 | Transfer to examiner for re-examination before appeal (zenchi) |
Free format text: JAPANESE INTERMEDIATE CODE: A911 Effective date: 20200311 |
|
TRDD | Decision of grant or rejection written | ||
A01 | Written decision to grant a patent or to grant a registration (utility model) |
Free format text: JAPANESE INTERMEDIATE CODE: A01 Effective date: 20200407 |
|
A61 | First payment of annual fees (during grant procedure) |
Free format text: JAPANESE INTERMEDIATE CODE: A61 Effective date: 20200501 |
|
R150 | Certificate of patent or registration of utility model |
Ref document number: 6707249 Country of ref document: JP Free format text: JAPANESE INTERMEDIATE CODE: R150 |
|
R250 | Receipt of annual fees |
Free format text: JAPANESE INTERMEDIATE CODE: R250 |
|
R250 | Receipt of annual fees |
Free format text: JAPANESE INTERMEDIATE CODE: R250 |