JP6707249B2 - Micro-tomography visualization method and system - Google Patents

Description

The present invention relates to methods and systems that enable microtoming of physical quantities within structures.

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.

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.

JP, 2006-284249, A

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.

It is a figure showing roughly composition of a fault visualization system concerning an example. It is a figure showing the concept of the micro tomography method concerning an example. It is a figure which shows roughly the processing procedure by the recursive cross correlation method. It is a figure which shows roughly the processing procedure by a subpixel analysis. It is a figure which shows a finite element model roughly. It is a figure which shows the processing procedure of reverse analysis schematically. It is a figure which shows the processing procedure of reverse analysis schematically. It is a figure which shows the processing procedure of reverse analysis schematically. It is a figure which shows the processing procedure of reverse analysis schematically. It is a flow chart which shows a flow of physical quantity visualization processing performed by a control operation part. 11 is a flowchart showing in detail the deformation vector calculation process of S12 in FIG. 11 is a flowchart showing in detail the inverse analysis process of S14 in FIG. It is a figure which shows the verification method by a numerical experiment. It is a figure which shows the verification method by a numerical experiment. It is a figure which shows the verification method by a numerical experiment. It is a figure which shows a verification result. It is a figure which shows the error with the correct value by the analysis of a present Example.

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.

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.

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.

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.

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.

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.

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.

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.

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.

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.

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.

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.

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.

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.

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.

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.

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.

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).

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.

By the way, the coherence length l _c, which is the resolution in the optical axis direction (depth direction) of OCT, is determined by the autocorrelation function of the light source. Here, the coherence length l _c can be represented by the following formula (1) with the half width at half maximum of the comprehensive line of the autocorrelation function.
Here, λ _c is the central wavelength of the laser beam, and Δλ is the full width at half maximum of the laser beam.

On the other hand, the resolution in the direction perpendicular to the optical axis (beam scanning direction) is set to 1/2 of the beam spot diameter D based on the focusing performance of the focusing lens. The beam spot diameter D can be expressed by the following equation (2).
Here, d _B is the beam diameter incident on the condenser lens, f is a focal point of the condenser lens. As described above, although the resolution by OCT has a limit, in the present embodiment, the introduction of sub-pixel analysis, which will be described later, makes it possible to display a tomographic display of the temperature distribution on a microscale. The details will be described below.

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.

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.

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).

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.

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.

(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.

The cross-correlation method is a method of evaluating the similarity of local speckle patterns from the correlation value R _ij based on the following equation (3). Therefore, as shown in FIG. 3A, regarding the OCT images before and after the photographing, the inspection region S1 to be the inspection target of the similarity is set in the previous tomographic image (Image1), and the subsequent tomographic image (Image2) is set. Is set to the search area S2 that is the search range of the similarity.
Here, as the spatial coordinates, the Z axis is set in the optical axis direction and the X axis is set in the direction perpendicular to the optical axis. f(X _i ,Z _j ) and g(X _i ,Z _j ) are the inspection area S1 (N _x ×N _z pixels) of the center position (X _i ,Z _j ) set in the OCT image before and after deformation. Shows the speckle pattern.

The correlation value distribution R _i,j (ΔX, ΔZ) in the search region S2 (M _x ×M _z pixels) is calculated, and pattern matching is performed as shown in FIG. Actually, as shown in the following formula (4), the movement amount U _i,j that gives the maximum correlation value is determined as the deformation vector before and after the deformation.
Note that f ⁻ and g ⁻ represent average values of f(X _i , Z _j ) and g(X _i , Z _j ) in the inspection area S1.

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.

In addition to this, the following equation (5) is used to set a threshold value using the average deviation σ of a total of nine deformation vectors including eight surrounding coordinates centered on the coordinate being calculated to eliminate the error vector. Suppress error propagation associated with recursive processing.
Here, Um represents the median of the vector amount, and the coefficient κ serving as the threshold value is set arbitrarily.

(Adjacent cross-correlation multiplication)
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).

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.

(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.

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.

The windward gradient method calculates the movement of the point of interest in the inspection area S1 not only with the pixel precision shown in FIG. 4A but also with the subpixel precision shown in FIG. 4B. Each grid in the figure represents one pixel. Actually, it is considerably smaller than the tomographic image shown in the figure, but it is shown in a large size for convenience of explanation. The upwind gradient method is a method for formulating the change in the luminance distribution before and after the minute deformation by the luminance gradient and the movement amount. When f is the luminance, the following equation (6) that Taylor-expands the minute deformation f(x+Δx, z+Δz) ( Represented as 7).

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.

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.

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.

The position of the point of interest before (after) the transformation is determined by the amount of subpixel movement when the parabolic approximation is performed. Eight coordinates that surround the position of the point of interest are used, and the brightness gradient is calculated by the ratio thereof. Specifically, the following formula (8) is used. The amount of movement was determined by applying the least squares method using the luminance gradient calculated in this way and the luminance change.

(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.

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.

More specifically, the calculation is executed in the following procedure. First, the bicubic function interpolation method is applied to the luminance distribution of the OCT tomographic image before material deformation to make the luminance distribution continuous. The bicubic function interpolation method is a method of reproducing spatial continuity of luminance information by using a convolution function obtained by piecewise approximating a sinc function by a cubic function. Originally, when measuring a continuous luminance distribution, the point spread function depending on the optical system is convolved. Therefore, the original continuous luminance distribution is restored by performing deconvolution using the sinc function. It When the discrete uniaxial signal f(x) is interpolated, the convolution function h(x) is expressed by the following equation (9).

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.

The inspection region S1 calculated in consideration of the expansion and contraction and the shear deformation is accompanied by the deformation along with the movement, as shown in FIG. Considering that the coordinates (x,z) at the integer pixel position in a certain inspection area S1 in the OCT tomographic image before material deformation moves to the coordinates (x ^* ,z ^* ) after deformation, the values of x ^* ,z ^* are It is represented by the following formula (10).

Here, Δx and Δz are the distances from the center of the inspection area S1 to the coordinates x and z, u and v are the deformation amounts in the x and z directions, respectively, and ∂u/∂x and ∂v/∂z are respectively x and Deformation amounts of the inspection region S1 in the z direction in the vertical direction, and ∂u/∂z and ∂v/∂x are deformation amounts of the inspection region S1 in the shear direction in the x and z directions, respectively. The Newton-Raphson method is used for the numerical solution, and the correlation value differential coefficient for 6 variables (u, v, ∂u/∂x, ∂u/∂z, ∂v/∂x, ∂v/∂z) is 0. That is, iterative calculation is performed to obtain the maximum correlation value. In order to improve the convergence of iterative calculation, the sub-pixel moving amount obtained by the upwind gradient method is used as the initial moving amount in the x and z directions. When the Hessian matrix for the correlation value R is H and the Jacobian vector for the correlation value is ∇R, the update amount ΔPi obtained in one iteration is represented by the following equation (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.

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.

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.

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.

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.

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.

(Overall stiffness equation including temperature terms by FEM)
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.
Here, σ _x is the x component of stress, σ _z is the z component of stress, and τ _xz is the shear stress. ε _x is strain in the x direction, ε _z is strain in the z direction, and γ _xz is shear strain. E is Young's modulus, ν is Poisson's ratio, α is linear expansion coefficient, and ΔT is temperature change.

Then, when the principle of virtual work is applied to the above equation (12) to transform it and expressed by using the element stiffness matrix [k] _6×6 , the following equation (13) is obtained.

In the present embodiment, when the FEM is used, each element divided into elements is represented by a triangular element shown in FIG. 5(A), and the relationship between adjacent elements is represented by a simple model shown in FIG. 5(B). Represent In the above equation (13), P and Q are equivalent nodal forces, and u and v are deformation vector components in the x direction and z direction of the node position, respectively. The subscripts 1 to 3 of P, Q, u, and v are the local node numbers of the triangular elements related to the element division shown in FIG. 5(A). a _{1 to} b ₃ are constants for approximating the physical quantity in the triangular element, and each is represented by the following formula (14).
Here, subscripts 1 to 3 of x and z are local node numbers of triangular elements related to element division.

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.

When the above equation (13) is expanded with respect to the term of temperature and a general equation is created using the simple model shown in FIG. 5(B), the following equation (15) is established.
Here, the superscript and the subscript represent the element number and the node number, respectively.

The above equation (15) is simplified and generalized into the following equation (16), which is a linearized observation equation.
Here, [A] _{(2min+2mb)×n} represents an observation matrix of the temperature change term, which depends on the [B] matrix determined by the shape function of the FEM, and those components are represented by the nodal coordinates. {ΔT} _n×1 is the temperature change in each element, [K] _{(2min+2mb)×(2min+2mb)} is the overall stiffness matrix, and {U} _{(2min+2mb)×1} is the deformation vector component at each node. Since this deformation vector component is obtained based on the OCT tomographic image, the right side is the observation vector.

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).

(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.

That is, if the meaning of the subscript of the matrix component is modified and the above equation (16) is rewritten, the following equation (17) is obtained.

In the process of solving this overall stiffness equation, it is modified to satisfy the boundary displacement condition. Here, the deformation vector components obtained at the nodes on the ROI boundary are u _l and v _l . At this time, in {U} _{(2min+2mb)×1} , u ₁ is the 21-th row and v ₁ is the _21-th row that represents the deformation vector component found at the node on the boundary. At the node of the ROI boundary, as a boundary displacement condition (displacement constraint boundary), the entire stiffness matrix [K] _{(2min+2mb)×(2min+2mb)} so that both the left and right sides of the known displacement amounts u _l , v _l of the node on the ROI boundary are satisfied. ₎ , and the observation matrix [A] _{(2min+2mb)×n} . Specifically, the components of [K] and [A] such that the displacements u _l and v _{l at the} nodes on the ROI boundary are u _l =u _l and v _l =v _{l on} the left and right sides, respectively. The components in the same row are modified according to the following equation (18).

Also, the 2l-1 and 2l lines corresponding to the nodes on the boundary of the vector {F} are also modified to u _l and v _l so that both the left and right sides are satisfied. Further, if the ROI is set inside the object and the volume force does not act due to the static deformation, the equivalent nodal force within the ROI except for the ROI boundary with respect to {F} becomes zero. In addition, when the ROI boundary is set on the free surface, free deformation occurs, so that the surface force is {F}=0. Therefore, the above equation (17) becomes the following equation (19).

As described above, the boundary displacement condition is satisfied in the overall stiffness equation. In the above equation (19), the 2l-1 and 2l lines on both sides belong to the constraint nodes and satisfy the boundary displacement condition. Therefore, by omitting them, the components are rearranged and the matrix is reduced. When the observation matrix [A] thus obtained is represented as [A'] and the product of the overall stiffness matrix [K] and the deformation vector {U} is represented as {K'U'}, the following equation (20) is obtained. can get.

(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.

That is, the deformation vector is interpolated on each node using the tomographic distribution of the deformation vector acquired by OCT. Here, a predetermined radius (300 μm in this embodiment) centered on each node is set as an “influence radius”, and a deformation vector that falls within the influence radius is interpolated at that node. In this interpolation, the approximation function f(x,z) shown in the following equation (21) was obtained by the least square method, and the deformation vector at each node position was calculated according to the approximation function f(x,z).
Here, f(x,z) is an approximate function that is a two-dimensional quadratic polynomial. a to g are the unknown coefficients, and x and z are the starting point coordinate values of the deformation vector obtained by OCT.

The following expression (22) is applied to this approximation function f(x,z).
Here, Q is the sum of squares of the residual between the approximate function f(x,z) and the measured value y _i . The measured value y _i is a deformation amount component obtained by OCT, and N is the number of data of deformation vectors within the influence radius. An approximate function f(x,z) when Q becomes the minimum in the above equation (22) is obtained, and the deformation amount component on the node in the ROI is calculated according to the f(x,z).

Specifically, after Q is calculated, Q is partially differentiated with respect to each unknown coefficient of the approximation function to obtain the normal expression shown in the following Expression (23).

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.

(Inverse analysis)
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).

Inverse analysis of the temperature is executed using the following equation (25) based on the above equation (24).
In this example, the temperature fault distribution {ΔT} was calculated as an estimated value by multiplying both sides of the above equation (20) by the Moore-Penrose generalized inverse matrix of [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).

That is, the temperature distribution obtained by the inverse analysis is calibrated based on the following formula (26).
Here, T ^* is an estimated value calculated by inverse analysis, and T ^* _realized is a temperature value (actual value) after calibration. As conceptually shown in FIG. 9B, T _sample1 and T _sample2 are true values at arbitrary two points used for calibration, and T ^* _sample1 and T ^* _sample2 are corresponding estimated values corresponding to the true values. Is. It should be noted that although the linear (linear function) calibration is illustrated in the present embodiment, the calibration may be performed using a polynomial (polynomial approximation).

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.

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.

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).

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.

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.

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.

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.

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.

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.

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.

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”).

(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.

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.

Then, as shown in FIG. 13B, element division was performed by FEM forward analysis on this model. The material constants of each material forming the model are as shown in Table 1 below.

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.

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.

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.

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.

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. ..

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.

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.

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.

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.

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.

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.

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.

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.

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 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: The method of claim 1, wherein the region of interest is set inside the tomographic image. The method for visualizing micro-tomography according to claim 1 or 2, wherein the numerical analysis is a stress analysis considering deformation of the structure due to a factor due to a state change or a boundary condition change. 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. 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. The micro tomographic visualization method according to claim 1, wherein a tomographic image of the structure is obtained by using optical coherence tomography. 7. The method for visualization of micro-tomography according to claim 1, wherein the structure includes at least a polymer-based material or a biological tissue. 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.