A method for fast and robust imaging of body tissues by waveform inversion from ultrasound examination data

Abstract

A method for imaging body tissue that is fast and robust is provided. [Solution] The method comprises the steps of (a) acquiring ultrasonic inspection data, (b) generating one wave field using multiple CPUs, (c) integrating the squared error between the waveform recorded at the receiving point from the wave field and the observed waveform acquired in step (a) and fc in equation (5) from all CPUs, (d) re-dividing the corrected gradient amount and pseudo-Hessian matrix in the x direction into regions with an equal number of CPUs, (e) re-dividing in the y direction, (f) calculating the step length required for model updating of full waveform inversion using the quasi-Newton method, and (g) repeating steps (b) to (f) to obtain imaging of the viscous Q value in addition to the sound propagation speed. TIFF0007607881000021.tif12145 Here, fc is a cost function, which is the sum of squared residuals between the observation data (d_obs) at all receiving points and all times and the synthetic data (d_syn) calculated using the initial model, and Ns, Nr, and Nt are the number of transmitting points, receiving points, and number of time samples.

Description

[Technical field]

[0001] The present invention relates to full waveform inversion of data acquired by ultrasonic medical equipment, which images the inside of the body using the sound propagation speed and attenuation Q value, unlike the conventional image (echo image) of a reflecting surface (or scattering point). The subject matter claimed in the present invention is a method for realizing high-speed calculation of numerical operations, which is essential for practical application of full waveform inversion for ultrasonic examination data. [Background technology]

[0002] Full waveform inversion is an underground visualization technique that applies wave theory developed in the field of seismology (Non-Patent Documents 1 and 2). In principle, full waveform inversion can be applied to data acquired by ultrasonic medical examination equipment, and several previous patents have been published in conjunction with the examination equipment (Patent Documents 1 to 6).

[0003] Imaging the inside of the body with ultrasound is preferable to X-rays, CT scans, and mammograms, which use X-rays, because there is no risk of radiation exposure. An example of ultrasound medical equipment is the echo machine, which can provide images in real time. However, echo examinations require repeated testing by experienced doctors and technicians to estimate problem areas.

[0004] In mammography, the most common method of breast cancer screening, it is difficult to distinguish between dense breast tissue and cancerous tissue in dense breasts. Full waveform inversion can image body tissue using two physical quantities, sound propagation velocity and attenuation Q value, making it possible to numerically distinguish between breast tissue and cancerous tissue, even in dense breasts. Although MRI has a slightly lower resolution than X-ray CT, there is no risk of radiation exposure and the computational cost is low enough that imaging results can be obtained on the same day of the examination.

[0005] As described in Patent Document 5, to be comparable to MRI, full waveform inversion must provide real-time images, and to provide real-time images, the turnaround time must be reduced by an order of magnitude. In addition, when using data from a limited frequency band (especially high frequency bands that require a propagation distance of more than a few wavelengths relative to the size of the target), the results of full waveform inversion are highly dependent on the initial model, making it difficult to ensure robustness. Full waveform inversion in the time domain has a significantly improved robustness compared to that in the frequency domain, but it has not been thoroughly studied due to the problem of computation time. Due to the problems described above, full waveform inversion has not yet been put to practical use as an ultrasonic tissue imaging technique. [Prior art documents] [Patent documents]

[0006] [Patent Document 1] JP 2020-18789 A [Patent Document 2] Patent Publication No. 2021-19839 [Patent Document 3] JP 2018-528829 A [Patent Document 4] JP 2020-501648 A [Patent Document 5] Special Publication No. 2020-501735 [Patent Document 6] International Publication No. 2013-116854 [Non-patent literature]

[0007] [Non-Patent Document 1] Lailly, P., “The seismic inverse problem as a sequence of before stack migrations”, Converence on Inverse Scattering, Theory and application, Society for Industrial and Applied Mathematics, 1983, p206-220 [Non-Patent Document 2] Tarantola, A., “Inversion of seismic reflection data in the acoustic approximation”, Geophysics, 1984, vol. 49, no. 8, pp. 1259-1266 [Non-Patent Document 3] Masashi Kato, "Study on Numerical Simulation of Ground Vibration Caused by Train Traffic Based on Viscoelastic Wave Theory", Doctoral Dissertation, Graduate School of Engineering, Kyoto University, 2008 [Non-Patent Document 4] Virieux J., “P-SV wave propagation in heterogeneous media: Velocity-stress dinite-difference method” Geophysics, 1986, vol.51, P889-901 [Non-Patent Document 5] Yang, P., and 3 others, “Wavefield reconstruction in attenuating media: A checkpointing-assisted reverse-forward simulation method”, Gephysics, 2016, vol. 81, no. 6, pp. R349-362 [Non-Patent Document 6] Yang, P., and 3 others, “A review on the systematic formulation of 3-D multiparameter full waveform inversion in viscoelastic medium”, Geophysical Journal International, 2017, vol. 207, p129-149 [Non-Patent Document 7] Shin C., and 2 others, “Improved amplitude preservation for prestack depth migration by inverse scattering theory.”Geophysical Prospecting, 2001, vol. 49.5, p592-606 Summary of the Invention [Problem to be solved by the invention]

[0008] In view of the above problems in the prior art, an object of the present invention is to provide a method for imaging body tissues quickly and robustly using waveform inversion from ultrasonic inspection data. [Means for solving the problem]

[0009] First, to improve the robustness of the full-waveform inversion, we adopt a time-domain computational method. In the time-domain full-waveform inversion, the image pixel size needs to be more than 30 pixels per wavelength. This leads to an increase in computation time, but has the advantage of increasing the image resolution. Instead of the acoustic wave equation cited in previous patents, we adopt the acoustic wave equation with viscosity. This makes it possible to provide the viscous Q value as an imaging result in addition to the acoustic wave propagation velocity.

[0010] Next, to speed up the calculations, parallel calculations are performed on multiple computers. Since each computer can be configured with a general-purpose CPU, the cost is low even if multiple computers are used simultaneously. However, full waveform inversion requires a huge amount of memory space, so the memory must be divided. In other words, the means for speeding up full waveform inversion calculations can be summarized as incorporating a memory space division method into full waveform inversion, as follows.

[0011] (1) When solving the wave equation (for both forward and backward propagation), the memory is divided at equal intervals in the x and y directions. In the divided memory space, the wave equation for the initial model is solved using the finite difference method. For each time step of the finite difference method, data is exchanged at the ends of the divided areas. At the same time, a cost function, which is the error with the observed data, is calculated for each area, and then these are collected from the divided memory areas into one and added together. (2) Calculate the correction gradient and pseudo Hessian for the area divided in (1). (3) Apply Gaussian smoothing to the correction gradient and pseudo-Hessian matrix in the y direction, dividing the x direction into equal intervals. (4) Following (3), Gaussian smoothing is applied in the x direction, with equal divisions in the y direction only. (5) While keeping the memory area divided at equal intervals in the y direction in (4), the correction amount gradient is corrected using the pseudo Hessian matrix. (6) While dividing the y-direction in (4) at equal intervals, calculate the norm of the correction gradient, collect the results from the divided memory areas, and add them together. (7) The direction vector of the correction amount is calculated from the norm obtained in (5), and the correction amount is multiplied by the correction amount gradient to obtain the correction amount of the initial model. The obtained correction amount is used as the initial model to obtain an updated model. (8) Return the memory partition to the form of (1) and repeat (1) to (7). (9) When the cost function converges to a certain threshold or the number of iterations reaches a set upper limit, the update process is terminated and the final imaging result is obtained.

[0012] That is, the method of imaging body tissues using waveform inversion from ultrasound examination data is as follows: (a) obtaining ultrasound examination data; (b) for the acoustic wave equation with viscosity, in memory spaces of the multiple computers (CPUs) divided independently at equal intervals in the x and y directions, data exchange of displacement velocity field and sound pressure field is performed between adjacent CPUs in the x direction and between adjacent CPUs in the y direction, and one wave field is generated by the multiple computers (CPUs); (c) a step of integrating, from all CPUs, the squared error between the waveform recorded at the receiving point from the wave field and the observed waveform obtained in the step (a) and fc in the following formula (5); (d) The correction gradient and the pseudo-Hessian matrix calculated in the memory space divided in the step (b) are After dividing the memory space into equal parts in the x direction equal to the number of CPUs, Apply a Gaussian filter in the y direction Use The stage, (e) dividing the memory space divided at equal intervals only in the x direction in the step (d) into a memory space divided at equal intervals in the y direction equal to the number of CPUs, and then applying a Gaussian filter in the x direction; (f) calculating the step length required for updating the model of the full waveform inversion by the quasi-Newton method in the memory space divided at equal intervals in the y direction in the step (e); (g) repeating steps (b) to (f) to obtain imaging of the viscous Q value in addition to the acoustic wave propagation velocity; The viscous acoustic wave equation is characterized in being expressed by simultaneous differential equations of the following formulas (1) to (4). [Formula 1] TIFF0007607881000002.tif11128 [Formula 2] TIFF0007607881000003.tif13149 [Formula 3] TIFF0007607881000004.tif11128 [Formula 4] TIFF0007607881000005.tif11128Here, the superscript · denotes the derivative in the time direction, ρ and V are the density and sound wave propagation speed of the medium, respectively, and σ, v x , v y are the sound pressure, the displacement velocity in the x direction, and the displacement velocity in the y direction, respectively. γl is a memory variable, τl is the relaxation time converted from L frequencies (ωl) that define the viscous Q value, L is the number of relaxation mechanisms (here, 3), [Formula 5] TIFF0007607881000006.tif12145Here, fc (a scalar quantity) is a cost function, which is the sum of squared residuals between the observational data (d_obs) and the synthetic data (d_syn) calculated using the initial model at all times at all receiving points, and Ns, Nr, and Nt are the number of transmitting points, the number of receiving points, and the number of time samples, respectively. Effect of the Invention

[0013] According to the present invention, a method is provided for fast and robust imaging of body tissue by waveform inversion from ultrasound examination data. [Brief description of the drawings]

[0014] [Figure 1] FIG. [Diagram 2] (a) is an explanatory diagram of the numerical model, and (b) is the physical quantities of body tissue based on experiments (relationship between viscosity Q value and sound wave propagation velocity V). [Diagram 3] Diagram of a wavelet: (a) waveform, (b) spectrum. [Figure 4] FIG. 1 is an explanatory diagram of a staggered grid arrangement. [Diagram 5] A snapshot of the wave field (64 μs after emission). [Figure 6] An explanatory diagram of memory space division: (a) no division, (b) 2 divisions, (c) 4 divisions, (d) 8 divisions, (e) 16 divisions. [Figure 7] FIG. 13 is an explanatory diagram of a data exchange diagram between adjacent CPUs in the x direction. [Figure 8]FIG. 13 is an explanatory diagram of a data exchange diagram between CPUs adjacent in the y direction. [Figure 9] FIG. 1 is a diagram of observation data. [Figure 10] Illustration of the initial model for full waveform inversion. (a) V is the acoustic wave velocity, (b) the viscous Q value, and (c) ρ is the medium density. [Figure 11] FIG. 13 is a diagram of synthetic data calculated with the initial model. [Figure 12] 1 is a plot of the residuals between observed and synthetic data. [Figure 13] This is a snapshot of the wave field after backpropagating the residual (64 μs after the transmission time). [Figure 14] Graph of correction amount gradient: (a) ρV2, (b) Q-1. [Figure 15] An explanatory diagram of how the memory space to be calculated is divided: (a) a four-way division of transmission points 1 to N / 4, (b) an x-directional division of the memory space to be calculated, and (c) a y-directional division of the memory space to be calculated. [Figure 16] Image of the results: (a) sound wave propagation velocity V, (b) viscous Q value. [Figure 17] 1 is a graph of a cost function. [Figure 18] FIG. 13 is an explanatory diagram of a model update. [Figure 19] 2 is a flowchart of the present invention. [Figure 20] This is a benchmark graph. [Figure 21] Snapshots of the frequency domain wave field: (a) forward propagation and (b) backward propagation. [Figure 22] An illustration of the correction gradient in the frequency domain: (a) V wave velocity, (b) viscous Q value. [Figure 23] 2 is a flowchart of the present invention. DETAILED DESCRIPTION OF THE PREFERRED EMBODIMENTS

[0015] To explain the morphology, a numerical simulation is performed for the numerical model and piezoelectric element arrangement shown in Figure 1. Figure 2 shows the physical quantities of body tissue based on experiments. The sound propagation speeds of water and tumors are similar, but there is a large difference in the viscous Q value. This suggests that when a tumor is present in an area of ​​body tissue with a high water content (blood vessels and mammary glands), it is difficult to detect the tumor using only imaging of the sound propagation speed, and that the imaging results of the viscous Q value play a major role in tumor detection.

[0016] The ultrasonic signal shown in FIG. 3 is an ultrasonic signal with a narrow band of about two octaves or less (in the case of FIG. 2, the range is approximately 200 kHz to 800 kHz, but the band can be selected at any frequency) that can be generated by a general-purpose piezoelectric element.

[0017] In Fig. 1, the geometry is such that 256 piezoelectric elements are grouped into one segment, and the body tissue is surrounded by eight segments (a total of 2048 piezoelectric elements in the clockwise direction). From a practical standpoint, there are some gaps between the segments, and data cannot be acquired in these gap sections (radius 11.7 cm, receiver interval 0.336 mm, gap between segments 5.71 mm to 6.38 mm).

[0018] To understand full waveform inversion, consider the case where an ultrasonic signal having the waveform shown in FIG. 5 is emitted from a piezoelectric element at the position indicated by the arrow (TR1) in FIG.

[0019] The pressure wave emitted from the piezoelectric element is expressed by four simultaneous differential equations: Equation (1) Hooke's law, Equation (2) memory variable (γ), and Equations (3) and (4) equations of motion. The acoustic wave equation with viscosity is expressed by the finite difference method, and the sound pressure (σ) and displacement velocity (v x , v y ) is calculated. [Formula 1] TIFF0007607881000007.tif11128 [Formula 2] TIFF0007607881000008.tif13149 [Formula 3] TIFF0007607881000009.tif11128 [Formula 4] TIFF0007607881000010.tif11128

[0020] The superscript · denotes the derivative in the time direction. ρ and V are the density of the medium and the sound wave propagation speed, respectively. σ, v x , v y are the sound pressure, the displacement velocity in the x direction, and the displacement velocity in the y direction, respectively. L is the number of relaxation mechanisms (assumed to be 3). A rigorous derivation method and the L frequencies (ω l ) to the relaxation time (τ l The method for converting this to the above is shown in the inventor's (Non-Patent Document 3).

[0021] The finite difference method is a technique for solving the above differential equations (1) to (4) by directly discretizing them (∂x→Δx; ∂t→Δt). In order to improve the calculation accuracy of the finite difference method, a staggered grid arrangement shown in FIG. 4 is adopted. The staggered grid arrangement was devised by the inventor's supervisor, Virieux, J. (Non-Patent Document 4). The discretization accuracy adopted is fourth order in the spatial direction and second order in the time direction.

[0022] Figure 5 shows a snapshot of the sound pressure field 64 μs after the transmission time of the forward propagation calculation.

[0023] To speed up the calculation, Memory space is divided among multiple computers at the dividing points in Fig. 6 (Fig. 6 shows the case of no division and the case of division into 16 computers, but the number of divisions can be freely set according to the capabilities of the CPU and network). Fig. 7 and Fig. 8 show diagrams of data exchange between CPUs.

[0024] Figure 9 shows the time series of sound pressure received by 2048 piezoelectric elements.

[0025] Next, we will explain how to perform imaging by full waveform inversion from simulated observed waveforms (Figure 8) observed at an unknown target (for the sake of explanation, the models in Figures 1 and 2 will be used as the target). Since full waveform inversion requires an initial model, an initial model of the sound velocity is obtained from travel time tomography, a common method for obtaining initial models. Since the viscous Q value cannot be obtained from travel time tomography, the initial model is made up of water and fat from Figure 2, excluding the tumor and skin. The density is assumed to be constant at the value of water, and is not solved by full waveform inversion.

[0026] The initial model obtained by the travel-time tomography method is shown in Fig. 10. By calculating this "forward propagating wave field in the initial model", the synthetic observation data shown in Fig. 9 is created at the observation point (Fig. 11). Figure 12 shows the residuals between the observation data obtained with the true model (Figure 9) and the synthetic observation data obtained with the initial model (Figure 11). With the initial model, sound waves can be synthesized in the same way as the observation data when there is no obstruction by body tissue between the transmitting and receiving points. On the other hand, since the initial model does not have clear boundaries of body tissue, it can be seen that almost no reflected waves are synthesized.

[0027] The residual is calculated by calculating the sum of squared residuals between the observed data (d_obs) and the synthetic data (d_syn) calculated using the initial model at all times for all receiving points using the following formula 5. [Formula 5] TIFF0007607881000011.tif12145Ns, Nr, and Nt are the number of transmission points, reception points, and time samples, respectively. fc is a scalar quantity, and is also called the cost function because it is the objective function to be minimized by the subsequent quasi-Newton method. The cost function is the sum of the divided calculations. Memory space Since it is necessary to keep the values ​​of all the computers together, the values ​​of all the computers are integrated so that all the computers have the same value.

[0028] This observed waveform residual is back-propagated (Figure 13 shows a snapshot of the wave field after back-propagating the residual (64 μs after the transmission time)). The cross-correlation is calculated in the time direction at the imaging point for the "forward propagating wave field in the initial model" and the "wave field generated by back-propagating the observed waveform residual". The zero-lag value of this cross-correlation is the gradient of the model correction amount in full waveform inversion. The forward propagating wave field is not stored in computer memory, but rather the total correction amount gradient is calculated computationally efficiently by performing a back-calculation of the forward propagating wave field as shown in the inventor's non-patent document 5. At the same time, the pseudo-Hessian matrix is ​​also calculated.

[0029] In the present invention, imaging of the viscous Q value is performed in addition to the acoustic propagation velocity. To do this, it is necessary to decompose the total correction amount gradient into the correction amount gradients of the acoustic propagation velocity and the viscous Q value. This decomposition method is not shown in Non-Patent Documents 1 and 2. Therefore, the inventor has shown this method in Non-Patent Document 6.

[0030] The above-mentioned forward propagation, back propagation, and calculation of the modified gradient magnitude are illustrated with the following comparative examples to illustrate the technical features of performing full waveform inversion in the time domain rather than the frequency domain (as is used in most prior art).

[0031] Once the correction gradients for the sound wave propagation velocity and the viscous Q value are obtained, a Gaussian filter (smoothing) is applied to each correction gradient and the pseudo-Hessian matrix. Memory Space As it is, it is not possible to calculate it on individual computers. Memory space Divide in the x direction and perform Gaussian smoothing in the y direction. Next, calculate Memory space Divide in the y direction and perform Gaussian smoothing in the x direction (Figure 15). Due to the nature of Gaussian smoothing, the result will be the same even if the order (x, y) is reversed.

[0032] After Gaussian smoothing, the inverses of the two pseudo-Hessian matrices for the sound velocity and the viscous Q value (H -1) to correct the correction gradient. This correction operation using a pseudo-Hessian matrix can be said to have the effect of increasing the convergence of the correction gradient to a solution by using only the diagonal elements of the Hessian matrix, which have a dominant effect on the correction gradient, without explicitly solving the inverse matrix of the Hessian matrix, which requires a huge amount of calculation (Non-Patent Document 7). However, it only makes the calculation possible numerically, and is theoretically incomplete, but works well empirically. In addition, it is expected to have the effect of suppressing the generation of locally erroneous images caused by the heterogeneity of the geometry.

[0033] [Formula 6] Δm=-H -1 G (6) A corrected model can be obtained by adding the correction amount (Δm) of the above formula 6 obtained in this way to the initial model. However, since this correction method is a linear approximation, the calculation is repeated several tens to several hundreds of times with the corrected model to obtain the final imaging. When correcting a model, the correction amount gradient is generally multiplied by a step length (α) to speed up the convergence. k ) to the k+1th model (m k+1 ) is expressed by the following Equation 7. [Formula 7] m k+1 = m k + αΔm k (7) There is an idea to obtain the step length from a cost function (evaluated only at the receiving point position) (as adopted in Patent Document 5), but it is better to evaluate it on the whole system, so the product (scalar) of the correction amount gradient and the correction amount after correction (8) below is used. TIFF0007607881000012.tif9128 is the correction amount gradient of (9) below. 2 Power(scalar) TIFF0007607881000013.tif9128 is used. i is the grid number and N is the total number of grids. The cost function is used only to determine whether convergence has occurred. Up to this point, the method has been steepest descent, and α leads the objective function to a minimum solution in the tangential direction. However, the correction gradient calculated in the time domain has little local error, so it is possible to further speed up convergence by adopting a quasi-Newton method called the limited-memory Broyden Fletcher Goldfarb Shanno (LBFGS) method and correcting α in the tangential direction of the quadratic function approximation of the objective function.

[0034] However, when the number of iterations is k=1, there is no history of the cost function, so the model is updated using the steepest descent method. Furthermore, when the number of iterations is k=2, the model is updated using the conjugate gradient method. As a result, from k=3 onwards, the cost function and correction amount required for the LBFGS method can be obtained at least twice.

[0035] When performing repeated inversion, the divided memory space for each computer can be divided in the y direction after the Gaussian filter. In the calculation of model update using this inversion technique, amount slope Distribution 2 Sum of powers, correction amount slope Distribution This is done by comparing the sum of the products of the correction amounts after correction with their history. Therefore, by integrating the L2-norm from the memory space divided in the y direction,

[0036] The imaging results obtained as described above are shown in FIG.

[0037] The cost function is plotted against the number of iterations in Figure 17. We can see that the cost function decreases with the number of iterations.

[0038] The updated models after 10th, 20th, and 30th iterations are arranged in Fig. 18. The tumor focus can be captured within 10 iterations, but 20 to 30 iterations are required for the skin to come into focus.

[0039] The sound wave propagation speed is imaged accurately. There is a certain error in the viscous Q value. Part of this error is due to the constant density model. But the majority of the problem is that the image is vibrating, as can be seen around the skin, so the wavelength of the sound wave is too short for the size of the target. However, the relative magnitude of the viscous Q value image preserves the true model, and is accurate enough for diagnostic purposes.

[0040] The flowchart is summarized in Figure 19. Memory space The decomposition method is illustrated along with the full waveform inversion flowchart.

[0041] Figure 20 shows the results of a benchmark test on calculation time and memory. In the benchmark test, the number of pixels in the x direction was 3005, the number of pixels in the y direction was 3005, the number of samples in the time direction was 5625, and the number of transmission points was 24. The computer used was a PC cluster consisting of two Intel XEON (E5-2690 v4 clock frequency 2.6 GHz, 14 core) CPUs and 256 GByte of memory, with 24 servers connected by LAN cables. The parallel numbers were 24, 48, 96, 192, and 384, and normalized values ​​for a parallel number of 24 are shown. Each is for finite difference method calculations. Memory Space The division numbers are 1, 2, 4, 8, and 16. Memory Space When the number of divisions is 2, the calculation time is reduced by more than half, and the improvement in efficiency within the computer due to the reduction in memory exceeds the loss due to communication. Memory Space When the number of divisions is 2 or 4, calculations require communication Memory space The length of the edges of is uniform across all computers, but for 8 and 16, the length of the communication edges becomes longer for computers that calculate the central part. This shows that the loss due to the heterogeneity of the communication volume becomes large. However, both the calculation time and memory can be reduced almost in proportion to the number of parallel nodes.

[0042] Comparative Example [Full waveform inversion in the frequency domain] We explain the difference between performing full waveform inversion in the time domain and the frequency domain. For the full-waveform inversion in the frequency domain, five frequency samples (251770, 375748, 501633, 625610, and 751495 Hz) shown in Fig. 2(b) were used. amount slope Distribution The calculations were done.

[0043] Figure 21 shows the wave field in the frequency (cosine transform) domain at a frequency of 501633 Hz. The wave field in the frequency domain has a striped pattern because it is a simple harmonic response. It is noteworthy that strong responses appear outside the body tissue and outside the circle of the transmitting and receiving point arrangement. However, we decided to perform backpropagation on the records of receiving points within an aperture width of 90° from the transmitting point. This is because if the records of all receiving points are backpropagated with an aperture width of 180°, the corrected amount slope Distribution Because the results were worse.

[0044] Corrected to Figure 22 amount slope Distribution In the frequency domain, the modified amount slope Distribution , which is the real part of the complex product of the forward propagation wave field and the backward propagation wave field, results in large values ​​remaining outside the circle of the transmitting and receiving points arrangement outside the body tissue. This is because the correction amount slope Distribution This is the cycle skip problem, which has long been a problem in full waveform inversion. Various smoothing and regularization methods have been devised to reduce this problem, but the inventors have not found a solution to this problem in their trials.

[0045] The tissues in the body are also corrected at regular intervals. amount slope AllocationIt becomes larger and shows a pattern like concentric circles. This means that the five frequency samples used are in phase during the same period. Unfortunately, even if the frequency samples are changed within the above range, the wavelength is very short (about 3 mm) compared to the size of the target (16 cm), so there will always be periods where the phases are in phase (strengthening or weakening each other). To avoid this, all frequency samples can be calculated, but attempting this is the same as a full waveform inversion in the time domain, which is much heavier than time domain calculations. The cycle skip problem can be solved by lowering the frequency band of the signal, but this requires generating sound waves that penetrate the target at 1 to 10 wavelengths. In order to mechanically generate low frequencies, the piezoelectric element must be enlarged, but this goes against the merit of ultrasonic medical equipment, which is that it is a low-cost inspection equipment.

[0046] There is also a secondary wave field generated by the gaps between the segments. amount slope Distribution It is clear that the reaction is large. Therefore, although the patent documents state that the method can be implemented in either the frequency domain or the time domain, the frequency domain is insufficient in that it is possible if all frequencies can be calculated. For example, in Patent Document 4, full waveform inversion is referred to as characterization of the boundary region, and is limited to reinforcing X-ray CT or travel-time tomography images (initial models), and anything that is not imaged to some extent in the initial model cannot be imaged by full waveform inversion. This implies that the robustness of full waveform inversion in the frequency domain is low. Therefore, the applicant claims that it is possible to obtain robust imaging results only in the time domain (Figure 23). Furthermore, it is believed that the probability of tumor detection can be improved by obtaining an image of the viscous Q value in addition to the sound wave propagation speed. [Industrial Applicability]

[0047] According to the present invention, a method for imaging body tissues quickly and robustly by waveform inversion from ultrasonic inspection data can be provided, which allows imaging of body tissues safely, at low cost, and with higher resolution than conventional inspection systems such as X-ray CT and MRI. In particular, the problem of the calculation time required for waveform inversion in the time domain can be solved by effectively using a large number of CPUs for parallel calculation, which has a great potential for industrial application since analysis results can be obtained in the same calculation time as conventional inspection systems. [Explanation of symbols]

[0048] T / R1 Snapshot (Fig. 5) Transmission point, first reception point T / R transmitting / receiving point (A) Tumor (B) Fat (C) Skin (D) Water V sound wave velocity Q Viscosity value ρ Density of the medium

Claims

1. 1. A method for imaging body tissue using waveform inversion from ultrasound examination data, comprising: (a) acquiring ultrasound examination data; (b) for the viscous acoustic wave equation, in a memory space divided independently at equal intervals in the x direction and the y direction of the plurality of computers (CPUs), data exchange of displacement velocity field and sound pressure field is performed between adjacent CPUs in the x direction and between adjacent CPUs in the y direction, and one wave field is generated by the plurality of computers (CPUs); (c) integrating, from all CPUs, the squared error between the waveform recorded at the receiving point from the wave field and the observed waveform acquired in the step (a) and f c in the following formula (5): (d) dividing the memory space into equal intervals in the x direction equal to the number of CPUs, and then applying a Gaussian filter in the y direction to the correction gradient and the pseudo-Hessian matrix calculated in the divided memory space in (b); (e) dividing the memory space divided at equal intervals only in the x direction in the step (d) into memory spaces at equal intervals in the y direction equal to the number of CPUs, and then applying a Gaussian filter in the x direction; (f) calculating a step length required for updating the model of the full waveform inversion by the quasi-Newton method in the memory space equally divided in the y direction in the step (e); and (g) repeating steps (b) to (f) to obtain imaging of the viscous Q value in addition to the acoustic wave propagation velocity; A method for high-speed and robust imaging of body tissue, characterized in that the viscous acoustic wave equation is expressed by the simultaneous differential equations of the following equations (1) to (4). [Formula 1] [Formula 2] [Formula 3] [Formula 4] Here, the superscript ・ denotes the derivative in the time direction, ρ and V are the density and sound wave propagation speed of the medium, respectively, and σ, v ｘ , v ｙ are the sound pressure, the displacement velocity in the x direction, and the displacement velocity in the y direction, respectively; γl is a memory variable; τl is a relaxation time converted from L frequencies (ωl) that define the viscous Q value; and L is the number of relaxation mechanisms (here, three). [Formula 5] Here, fc (a scalar quantity) is a cost function, which is the sum of squared residuals between the observation data (d_obs) and the synthetic data (d_syn) calculated using the initial model at all times at all receiving points, and Ns, Nr, and Nt are the number of transmission points, the number of receiving points, and the number of time samples, respectively.

2. The method for fast and robust imaging of body tissues according to claim 1, characterized in that in step (b), data is exchanged at each time step, the number of divisions of the divided memory space is set according to the computer environment consisting of the number of CPUs and memory space, the divided memory space is applied to all three wave field calculations, namely forward propagation and backward propagation of full waveform inversion and forward propagation unwinding calculation, and the correction amount gradient and pseudo-Hessian matrix are calculated while keeping the divided memory space as it is.

3. The method for fast and robust imaging of body tissues according to claim 2, characterized in that the sum of squares of the correction gradient calculated in step (f) shown in the following (9) and the sum of multiplications of the correction gradient and the pseudo-Hessian matrix shown in the following (8) are integrated from all CPUs. 〔８〕 〔９〕 Here, i is the grid number, N is the total number of grids, and H －１ is the inverse matrix of the pseudo-Hessian matrix and G is the correction gradient.

4. The ultrasound examination data is ultrasound breast cancer examination data, The method of claim 3, wherein the body tissue includes skin, breast, and cancerous tissue.

Images