Method for imaging body tissue rapidly and robustly using waveform inversion from ultrasound examination data

Abstract

In view of the problems and the like of the prior art, the purpose of the present invention is to provide a method for imaging body tissue rapidly and robustly using waveform inversion from ultrasound examination data.　The method for imaging body tissue rapidly and robustly using waveform inversion from ultrasound examination data according to the present invention comprises: (a) a step for acquiring ultrasound examination data; (b) a step for generating one wave field by means of a plurality of computers (CPUs); (c) a step for integrating squared errors (which constitute fc in expression (5)) from all CPUs, each squared error being the error between a waveform recorded at a reception point from the wave field and an observed waveform obtained in step (a); (d) a step for re-dividing a correction gradient and a pseudo Hessian matrix into regions with an equal number of CPUs in the x-direction; (e) a step for re-dividing from the uniform regions in the x-direction into regions with an equal number of CPUs in the y-direction; (f) a step for calculating a step length required for model updating in full-waveform inversion using a quasi-Newton method while maintaining the division of regions in the x-direction; and (g) a step for repeating the steps (b) to (f) to obtain imaging of a viscosity Q value in addition to a sound wave propagation velocity. Expression 5: Here, fc (scalar quantity) is a cost function representing the sum of squared residuals between the observed data (d_obs) and the synthetic data (d_syn) calculated from the initial model at all times of all reception points, and Ns, Nr, Nt are respectively the number of transmission points, the number of reception points, and the number of time samples.

Description

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

[0001] The present invention relates to full waveform inversion of data acquired by an ultrasonic medical device, which, unlike conventional images (echo images) of reflecting surfaces (or scattering points), images the inside of the body using the sound propagation velocity and attenuation Q value. What is claimed in this invention is a method for realizing high-speed numerical calculations, which are essential for practical use of full waveform inversion for ultrasonic examination data.

[0002] Full waveform inversion is a subsurface 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 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 mammography, which use X-rays, as there is no risk of radiation exposure. An example of ultrasound medical equipment is an ultrasound machine, which can provide images in real time. However, ultrasound examinations require repeated testing by skilled doctors and technicians to identify problem areas.

[0004] Mammography, the most common method of breast cancer screening, has difficulty distinguishing between dense breast tissue and cancerous tissue in dense breasts. Full waveform inversion can image body tissue using two physical quantities: sound wave velocity and attenuation Q value, making it possible to numerically distinguish between breast tissue and cancerous tissue, even in dense breasts. While MRI offers 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, full waveform inversion requires the provision of real-time images to rival MRI, and providing real-time images requires an order of magnitude reduction in turnaround time. Furthermore, when using data from a limited frequency band (especially high-frequency bands requiring propagation distances of several wavelengths or more relative to the target size), the results of full waveform inversion are highly dependent on the initial model, making it difficult to ensure robustness. Time-domain full waveform inversion significantly improves robustness compared to frequency-domain inversion, but has not been thoroughly studied due to issues with computational time. Due to these issues, full waveform inversion has not yet been put to practical use as an ultrasound imaging technique for body tissues.

[0006] Japanese Patent Application Publication No. 2020-18789 Japanese Patent Application Publication No. 2021-19839 Japanese Patent Application Publication No. 2018-528829 Japanese Patent Application Publication No. 2020-501648 Japanese Patent Application Publication No. 2020-501735 International Publication No. 2013-116854

[0007] Lailly, P. , “The seismic inverse problem as a sequence of before stack migrations”, Convergence on Inverse Scattering, Theory and application, Society for Industrial and Applied Mathematics, 1983, p206-220 Tarantola, A. , "Inversion of seismic reflection data in the acoustic approximation," Geophysics, 1984, vol. 49, no. 8, pp. 1259-1266; Kato, Masashi, "Study on numerical simulation of ground vibrations caused by trains based on viscoelastic wave theory," Kyoto University Graduate School of Engineering, Doctoral Dissertation, 2008; Virieux, J., "P-SV wave propagation in heterogeneous media: Velocity-stress dint-difference method," Geophysics, 1986, vol. 51, P889-901Yang, 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 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-149Shin C. , 2 others,“Improved amplitude preservation for prestack depth migration by inverse scattering theory.”Geophysical Prospecting, 2001, vol. 49.5, p592-606,

[0008] SUMMARY OF THE INVENTION In view of the above-mentioned problems of the prior art, an object of the present invention is to provide a method for imaging body tissues quickly and robustly by waveform inversion from ultrasound examination data.

[0009] First, to improve the robustness of full waveform inversion, a time-domain calculation method is adopted. In time-domain full waveform inversion, an image pixel size of 30 pixels or more per wavelength is required. This increases the calculation time, but it has the advantage of improving image resolution. Instead of the acoustic wave equation cited in previous patents, an acoustic wave equation with viscosity is adopted. This makes it possible to provide the viscosity Q value in addition to the acoustic wave propagation velocity as an imaging result.

[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 when 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 (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. Data is exchanged at the ends of the divided area for each time step of the finite difference method. At the same time, a cost function, which is the error with the observed data, is calculated for each area, and then collected from the divided memory area and added together. (2) The correction gradient and pseudo Hessian matrix are calculated using the area divided in (1). (3) Gaussian smoothing is applied to the correction gradient and pseudo Hessian matrix in the y direction. In this case, the area is divided at equal intervals only in the x direction. (4) Following (3), Gaussian smoothing is applied in the x direction. In this case, the area is divided at equal intervals only in the y direction. (5) While the memory area remains divided at equal intervals in the y direction in (4), the correction amount gradient is corrected using a pseudo-Hessian matrix. (6) While the memory area remains divided at equal intervals in the y direction in (4), the norm of the correction amount gradient is calculated, and the values ​​are collected from the divided memory area and added together. (7) The correction amount directional vector is calculated from the norm obtained in (5), and multiplied by the correction amount gradient to obtain the correction amount of the initial model. An updated model is obtained using the obtained correction amount as the initial model. (8) The divided memory area is returned to the form of (1), and steps (1) to (7) are repeated. (9) When the cost function converges below 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 for imaging body tissues by waveform inversion from ultrasound examination data is characterized by comprising the steps of: (a) acquiring ultrasound examination data; (b) generating one wave field using a 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 step (a), fc in the following formula (5); (d) re-dividing the corrected gradient and pseudo-Hessian matrix in the x direction into regions with an equal number of CPUs; (e) re-dividing the regions with an equal number of CPUs in the y direction from the regions with an equal number of CPUs; (f) calculating the step length required for model updating of full waveform inversion by the quasi-Newton method while the regions are still divided in the x direction; and (g) repeating steps (b) to (f) to obtain imaging of the viscous Q value in addition to the sound wave propagation velocity [Formula 5]. Here, fc (a scalar quantity) is a cost function, which is 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 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.

[0013] The present invention provides a method for fast and robust imaging of body tissue by waveform inversion from ultrasound examination data.

[0014] 1 is a diagram of the arrangement of transmitting and receiving points. (a) An explanatory diagram of a numerical model and (b) physical quantities of body tissues based on experiments (relationship between viscosity Q value and sound wave propagation velocity V). Wavelet diagram (a) waveform, (b) spectrum. An explanatory diagram of a staggered grid arrangement. A diagram of a snapshot of the wave field (64 μs after the transmission time). An explanatory diagram of domain division: (a) no division, (b) 2 division, (c) 4 division, (d) 8 division, (e) 16 division. An explanatory diagram of a data exchange diagram between adjacent CPUs in the x direction. An explanatory diagram of a data exchange diagram between adjacent CPUs in the y direction. A diagram of observed data. An explanatory diagram of the initial model of full waveform inversion: (a) V sound wave propagation velocity, (b) viscosity Q value, (c) ρ medium density. A diagram of synthetic data calculated with the initial model. A diagram of the residual between the observed data and synthetic data. A diagram of a snapshot of the wave field (64 μs after the transmission time) after the residual has been back-propagated. Diagram of correction gradient (a) ρV ２ , (b) Q －１ . An explanatory diagram of the division method of the calculation domain (a) 4 division of transmission points 1 to N, (b) division of the calculation domain in the x direction, (c) division of the calculation domain in the y direction. Diagram of the result image (a) V sound wave propagation speed, (b) viscosity Q value. Graph of the cost function. An explanatory diagram of model update. Flowchart of the present invention. Graph of benchmark. Diagram of snapshot of frequency domain wave field (a) forward propagation, (b) backward propagation. Diagram of correction gradient amount in the frequency domain (a) V sound wave propagation speed, (b) viscosity Q value. Flowchart of the present invention.

[0015] To explain the morphology, we performed a numerical simulation using 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 viscosity Q value. This suggests that when a tumor is present in a location with a high water content in body tissue (blood vessels or mammary glands), it is difficult to detect the tumor using only imaging of the sound propagation speed, and that the results of imaging the viscosity 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 any frequency band can be selected) that can be generated by a general-purpose piezoelectric element.

[0017] 1, the geometry is such that body tissue is surrounded by eight segments (a total of 2048 piezoelectric elements in a clockwise direction), with 256 piezoelectric elements in each segment group. 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 spacing 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 solved by the finite difference method to obtain the sound pressure (σ) and displacement velocity (v ｘ , v ｙ ) is calculated. [Formula 2] [Formula 3] [Formula 4]

[0020] The superscript "・" indicates the time derivative. ρ and V are the density of the medium and the sound wave propagation velocity, respectively. σ, v ｘ , v ｙ are the acoustic 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 three). A rigorous derivation method and the L frequencies (ω ｌ ) to the relaxation time (τ ｌ The method for converting the value of the saturation energy into the value of the saturation energy 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). To improve the calculation accuracy of the finite difference method, a staggered grid arrangement as shown in Figure 4 is adopted. The staggered grid arrangement was devised by Virieux, J. (Non-Patent Document 4), the inventor's supervisor. The discretization accuracy is fourth order in the spatial direction and second order in the time direction.

[0022] FIG. 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, the calculation domain is divided among multiple computers at the division points shown in Figure 6. (In Figure 6, the calculation domain is divided into 16 computers, starting from no division, but the number of divisions can be freely set depending on the CPU and network capabilities.) Figures 7 and 8 show diagrams of data exchange between CPUs.

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

[0025] Next, we will explain a method for performing imaging using 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 sound velocity is obtained from travel time tomography, a common method for obtaining an initial model. Since the viscosity Q value cannot be obtained from travel time tomography, the initial model is made up of water and fat, excluding the tumor and skin from Figure 2. The density is assumed to be constant at the value of water, and is not solved by full waveform inversion.

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

[0027] The residual sum of squares of the observed data (d_obs) and the synthesized data (d_syn) calculated using the initial model is calculated at all times at all reception points using the following formula 5: [Formula 5] Ns, Nr, and Nt are the number of transmission points, the number of reception points, and the number of 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. Since the cost function requires that the values ​​of the divided calculation domains be combined, the values ​​of all computers are integrated so that all computers have the same value.

[0028] This observed waveform residual is backpropagated (Figure 13 shows a snapshot (64 μs after the transmission time) of the wavefield after backpropagating the residual). The cross-correlation is calculated in the time direction at the imaging point for both the "forward propagating wavefield in the initial model" and the "wavefield generated by backpropagating 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 wavefield 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 wavefield 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 velocity. To do this, it is necessary to decompose the total correction gradient into the correction gradients of the acoustic velocity and the viscous Q value. This decomposition method is not disclosed in Non-Patent Documents 1 and 2. Therefore, the inventor has disclosed this method in Non-Patent Document 6.

[0030] The calculation of the forward propagation, backward propagation, and modified gradient magnitudes described above are presented 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 much of the prior art).

[0031] Once the correction gradients for the sound wave propagation velocity and viscosity Q value have been determined, a Gaussian filter (smoothing) is applied to each correction gradient and pseudo-Hessian matrix. Since Gaussian smoothing cannot be calculated on individual computers using the divided domains, the calculation domain is divided in the x direction and Gaussian smoothing is performed in the y direction. Next, the calculation domain is divided in the y direction and Gaussian smoothing is performed in the x direction (Figure 15). Due to the characteristics of Gaussian smoothing, the results are the same even if the order (x, y) is reversed.

[0032] After applying Gaussian smoothing, the inverse matrices of the two pseudo-Hessian matrices for the sound propagation velocity and the viscous Q value (H －１ The correction gradient (G) is corrected by dividing the correction gradient (G) by . This pseudo-Hessian matrix correction operation uses only the diagonal components of the Hessian matrix, which have a dominant effect on the correction gradient, and can be said to have the effect of improving the convergence of the correction gradient to a solution without explicitly solving the inverse matrix of the Hessian matrix, which requires a huge amount of calculation (Non-Patent Document 7). However, this only enables numerical calculation, 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 geometric inhomogeneity.

[0033] [Formula 6] Δm=-H －１ 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 the model, the correction amount gradient is generally multiplied by the step length (α) to speed up the convergence. After the kth iteration, the model (m ｋ ) to the k+1th model (m ｋ＋１ The calculation for updating to m is expressed by the following formula 7. ｋ＋１ = mｋ　 + αΔm ｋ (7) There is also an idea (adopted in Patent Document 5) that the step length is calculated from a cost function (evaluated only at the receiving point position), but it is better to evaluate it in the entire system, so the product (scalar) of the correction amount gradient and the correction amount (8) below after correction is used. The square of the correction gradient (scalar) of (9) below The value obtained by dividing by 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, we have used the steepest descent method, and α guides the objective function tangentially to a minimal solution. 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 known as 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 cost function history, 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 functions and correction amounts required for the LBFGS method can be obtained at least twice.

[0035] When performing repeated inversion, the divided domains for each computer can remain divided in the y direction after the Gaussian filter. The calculation of model update using this inversion technique is performed by comparing the sum of the squares of the correction gradients and the sum of the products of the correction gradients and the corrected corrections with their history. Therefore, it is necessary to integrate the L2-norm from the divided domains in the y direction so that it has the same value on all computers.

[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. It can be seen that the cost function decreases with the number of iterations.

[0038] Figure 18 shows the updated models after 10, 20, and 30 iterations. 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 velocity is imaged accurately. There is a certain error in the viscous Q value. Part of this error is due to the constant density model. However, the majority of the problem is that the image vibrates, as can be seen around the skin, and the sound wave wavelength is too short for the size of the target. However, the relative magnitude relationship 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. The decomposition method of the computational domain of the present invention is illustrated along with the flowchart of full waveform inversion.

[0041] Figure 20 shows the results of benchmark tests for calculation time and memory. The benchmark tests used 3,005 pixels in the x direction, 3,005 pixels in the y direction, 5,625 samples in the time direction, and 24 transmission points. The computers used were a PC cluster consisting of two Intel XEON (E5-2690 v4, 2.6 GHz, 14-core) CPUs with 256 GB of memory, forming one server. 24 servers were connected via LAN cables. The parallelism numbers were 24, 48, 96, 192, and 384, and the values ​​shown are normalized to 24. These represent the domain division numbers for finite difference calculations of 1, 2, 4, 8, and 16, respectively. When the domain division number was 2, the calculation time was reduced by more than half, demonstrating that the improvement in computer efficiency due to reduced memory requirements outweighed the loss due to communication. When the number of domain divisions is 2 or 4, the length of the sides of the calculation domain that require communication is uniform across all computers, but when the number of domain divisions is 8 or 16, the length of the communication sides becomes longer for computers that calculate the central part. This shows that the loss due to the non-uniformity of communication volume becomes large. However, both the calculation time and memory can be reduced almost in proportion to the number of parallel processes.

[0042] [Comparative Example] [Frequency Domain Full Waveform Inversion] The difference between full waveform inversion in the time domain and the frequency domain will be explained. The correction gradient was calculated for full waveform inversion in the frequency domain using the five frequency samples (251,770, 375,748, 501,633, 625,610, and 751,495 Hz) shown in Figure 2(b).

[0043] Figure 21 shows the wave field in the frequency (cosine transform) domain at a frequency of 501,633 Hz. The wave field in the frequency domain is a simple harmonic response, resulting in a striped pattern. It is noteworthy that strong responses appear outside the body tissue and outside the circle where the transmitting and receiving points are arranged. However, we decided to perform backpropagation on records from receiving points within an aperture width of 90° from the transmitting point. This is because backpropagating records from all receiving points with an aperture width of 180° would have resulted in poorer results for the correction gradient amount.

[0044] The correction gradient is shown in Figure 22. In the frequency domain, the correction gradient is the real part of the complex product of the forward propagation wave field and the backward propagation wave field, but large values ​​remain outside the body tissue, even outside the circle of the transmitting and receiving point arrangement. This is a completely different situation from the correction gradient calculated in the time domain in Figure 14. This is the cycle skip problem, which has long been considered an issue in full waveform inversion. Various smoothing and regularization calculations have been devised to alleviate this problem, but the inventors have not been able to find a solution in their trials.

[0045] The correction gradient also increases periodically within the tissue, forming a concentric pattern. This indicates the period in which the five frequency samples used are in phase. Unfortunately, even when the frequency samples are changed within the above range, the wavelength is very short (approximately 3 mm) compared to the target size (16 cm), so periods of phase alignment (constructive or destructive) are inevitable. To avoid this, all frequency samples could be calculated, but attempting this is equivalent to full waveform inversion in the time domain, which is significantly more time-consuming than time-domain calculations. The cycle skip problem can be solved by lowering the signal frequency band, but this requires generating sound waves that penetrate within 1 to 10 wavelengths of the target size. Mechanically generating low frequencies requires larger piezoelectric elements, which contradicts the purpose of ultrasound medical devices as low-cost examination equipment.

[0046] Furthermore, there is a secondary wave field generated by the gaps between segments. It can also be seen that the correction gradient responds significantly to this artificial noise. Therefore, although patent documents state that this method can be performed in either the frequency domain or the time domain, the frequency domain method is insufficient in that it is only possible if all frequencies can be calculated. For example, in Patent Document 4, full waveform inversion is referred to as characterizing the boundary region and is only used to reinforce X-ray CT or traveltime tomography images (initial model). Therefore, full waveform inversion cannot image anything that is not sufficiently imaged in the initial model. This implies the low robustness of frequency domain full waveform inversion. Therefore, the applicant claims that robust imaging results can only be obtained in the time domain (Figure 23). Furthermore, it is believed that obtaining an image of the viscous Q value in addition to the acoustic wave propagation velocity can improve the probability of tumor detection.

[0047] According to the present invention, a method for fast and robust imaging of body tissues by waveform inversion from ultrasound examination data can be provided, which enables safer, lower cost, and higher resolution imaging of body tissues than conventional examination systems such as X-ray CT and MRI. In particular, the present invention has great potential for industrial application, as it can solve the problem of the computation time required for waveform inversion in the time domain by effectively using a large number of CPUs for parallel computation, thereby obtaining analysis results in a computation time comparable to that of conventional examination systems.

[0048] T / R1: The first transmitting and receiving point of the snapshot (Fig. 5) T / R: Transmitting and receiving point (A) Tumor (B) Fat (C) Skin (D) Water V: Sound wave propagation speed Q: Viscosity value ρ: Density of the medium

Claims

1. A method for imaging body tissues by waveform inversion from ultrasonic examination data, comprising: (a) obtaining ultrasonic examination data; (b) generating one wave field with a plurality of computers (CPUs); (c) integrating, for each CPU, the squared error between the waveform recorded at the receiving point from the wave field and the observed waveform obtained in step (a); (d) re-dividing the correction gradient amount and the pseudo-Hessian matrix in the x direction into regions with an equal number of CPUs; (e) re-dividing the regions with an equal number of CPUs in the x direction into regions with an equal number of CPUs in the y direction; (f) calculating, by the quasi-Newton method, the step length required for model updating of full waveform inversion while keeping the regions divided in the x direction; and (g) repeating steps (b) to (f) to obtain imaging of the viscous Q value in addition to the sound wave propagation velocity, characterized by a fast and robust method for imaging body tissues. [Equation 5] Here, fc (scalar quantity) is a cost function, which is the sum of the squared residuals of the observed data (d_obs) and the synthetic data (d_syn) calculated with the initial model at all times for all receiving points. Ns, Nr, and Nt are the number of transmitting points, the number of receiving points, and the number of time samples, respectively.

2. In the step (b), data exchange is performed at each time step, the number of divisions of the divided region is set according to the computer environment (number of CPUs, memory), the divided region is applied to all three wave field calculations of forward propagation, backward propagation, and backward calculation of forward propagation of full waveform inversion, and the correction amount gradient and the pseudo-Hessian matrix are obtained while remaining in the divided region. The method for imaging a body tissue that is fast and robust according to claim 1.

3. After the step (d), a Gaussian filter is applied in the y direction, and after further dividing the x direction from a uniform region into a uniform region of the number of CPUs in the y direction and then applying a Gaussian filter in the y direction. The method for imaging a body tissue that is fast and robust according to claim 1.

4. The sum of squares shown in the following (9) of the amount of the corrected gradient calculated in the step (f), the sum of products of the amount of the corrected gradient and the pseudo-Hessian matrix shown in the following (8) are integrated from all CPUs, and the method for imaging a fast and robust body tissue according to claim 3, wherein: [8] [9] 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 corrected gradient.

5. The ultrasonic examination data is ultrasonic breast cancer examination data, and the body tissue includes skin, mammary gland, and cancerous tissue. The method for imaging a body tissue that is fast and robust according to claim 4.

Images