An ultrasonic cancellous bone imaging method and system based on full-wave inversion

Abstract

The application relates to the technical field of bone imaging, and provides an ultrasonic cancellous bone imaging method based on full-wave inversion, which comprises the following steps: S1: setting an imaged bone phantom model, emitting an ultrasonic pulse signal, recording the ultrasonic pulse signal passing through the model, and sampling the signal as a simulation experiment wave field signal; S2: establishing a simulation wave field, setting an initial sound velocity model as a current sound velocity model; S3: based on the current sound velocity model, carrying out forward modeling on the simulation wave field, carrying out incremental extrapolation of the forward wave field signal along the propagation direction of the ultrasonic pulse signal, recording the simulation wave field signal after completing one round of emission-reception; S4: constructing a loss function, solving the loss function by using a conjugate gradient method, and updating the current sound velocity model; and S5: judging whether the loss function converges or not, and outputting the current sound velocity model if the loss function converges. High-resolution ultrasonic imaging of bone and microstructure is realized, and the demand of bone medical treatment and diagnosis with high precision requirements is met.

Description

TECHNICAL FIELD

[0001] The present application relates to the technical field of bone imaging, and in particular to an ultrasonic cancellous bone imaging method and system based on full-wave inversion. BACKGROUND

[0002] Bone, as the main part of human hard tissue, is difficult for ultrasound to propagate in due to its higher difference in sound speed compared to surrounding soft tissue and multi-layer irregular structure, thereby greatly increasing the difficulty of obtaining high-resolution bone ultrasound images. Traditional bone imaging methods based on traditional pulse echo are not only subject to the requirement of prior sound speed model, but also cannot effectively image the microstructure of cancellous bone, making it difficult to meet the needs of high-precision diagnosis and treatment in clinical applications.

[0003] Full-wave inversion imaging, which utilizes full waveform information to invert the target, can theoretically achieve a spatial resolution of half the wavelength and has the characteristics of high resolution, and is often used in situations where high image accuracy is required. The principle is to establish a new forward simulation wave field based on the measured experimental wave field information, update the forward wave field using the simulation wave field information and the experimental wave field information, and finally make the forward wave field match the experimental wave field to complete the inversion process. Full-wave inversion has been widely used in the field of geophysics, mainly for analyzing underground structures and realizing energy and resource exploration.

[0004] For bone and surrounding soft tissue with a large difference in sound speed, the reflection phenomenon of ultrasound waves will be significantly enhanced when propagating to the interface, and the attenuation of ultrasound waves caused by multiple diffraction, scattering, non-elastic absorption and elastic mode conversion results in a very complex wave field for bone imaging. Most existing ultrasonic pulse echo bone imaging methods are difficult to achieve high-quality imaging of bone without prior sound speed model, and cannot effectively feedback the microstructure information inside the bone.

[0005] The traditional ultrasonic pulse echo bone imaging method has the following disadvantages: the pulse echo-based medical ultrasonic imaging needs to meet the uniform sound speed assumption, the sound speed of the bone hard tissue (2800-4000 m / s) is greatly different from the sound speed of the surrounding soft tissue (about 1540 m / s), the uniform sound speed assumption cannot be directly applied to bone imaging, otherwise the distortion and deformation of bone imaging will be caused, although part of the research has realized imaging by using the prior sound speed model, but in the experiment and clinic, the bone sound speed model is often unknown, and the method is not practical; the resolution of the pulse echo ultrasonic bone imaging method is limited, the axial resolution is equal to half of the spatial pulse length, and the lateral resolution depends on the width of the ultrasonic beam, which can be improved in actual application, but cannot be reconstructed for bone microstructure (bone trabecula, bone micro-hole). Although the traditional time-domain full-wave inversion can realize high-resolution imaging of bone, it directly processes time-domain data, which leads to a large amount of data and consumes a lot of computing resources, and is difficult to apply in real time in actual scenes. SUMMARY

[0006] In view of the above problems, the purpose of the present application is to provide an ultrasonic cancellous bone imaging method and system based on full-wave inversion, which can realize high-resolution ultrasonic imaging of bone and its microstructure, and meet the needs of bone medical treatment and diagnosis with high precision requirements.

[0007] The above invention purpose of the present application is realized by the following technical scheme:

[0008] An ultrasonic cancellous bone imaging method based on full-wave inversion, comprising the following steps:

[0009] S1: setting a to-be-imaged bone phantom model including surrounding soft tissue, emitting an ultrasonic pulse signal to the to-be-imaged bone phantom model, recording the ultrasonic pulse signal representing the sound field passing through the to-be-imaged bone phantom model, and sampling the ultrasonic pulse signal into an analog experimental wave field signal;

[0010] S2: establishing a simulation wave field and setting an initial sound speed model based on the simulation wave field, taking the initial sound speed model as a current sound speed model;

[0011] S3: based on the current sound speed model, performing forward modeling on the simulation wave field, performing incremental extrapolation of the forward wave field signal along the propagation direction of the ultrasonic pulse signal, recording all simulation wave field signals after completing a round of emission-reception of the ultrasonic pulse signal;

[0012] S4: constructing a loss function according to the analog experimental wave field signal and the simulation wave field signal, and solving the loss function by using a conjugate gradient method to obtain an optimal solution, and updating the current sound speed model based on the iterative search direction of the conjugate gradient method;

[0013] S5: determining whether the loss function converges, outputting the current sound velocity model as a final sound velocity model if the loss function converges, completing the full-wave inversion process, otherwise repeating steps S3 and S4 until the loss function converges.

[0014] Further, in step S1, the ultrasonic pulse signal is emitted to the bone phantom model to be imaged, and the ultrasonic pulse signal representing the sound field passing through the bone phantom model to be imaged is recorded, specifically:

[0015] Two identical linear array ultrasonic transducers are arranged in parallel, each of the linear array ultrasonic transducers includes N array elements, and the bone phantom model to be imaged is arranged between the two linear array ultrasonic transducers.

[0016] All array elements on one side of the linear array ultrasonic transducer sequentially emit the ultrasonic pulse signal with a fixed center frequency to the bone phantom model to be imaged, while simultaneously recording the ultrasonic pulse signal received on the other side of the linear array ultrasonic transducer.

[0017] Further, in step S1, the ultrasonic pulse signal is sampled into the simulated experimental wave field signal, specifically:

[0018] The ultrasonic pulse signal is sampled into N groups of simulated experimental wave field signals with a size of N*N t by k-space pseudospectral method, where N t is the sampling number.

[0019] Further, in step S1, the simulated experimental wave field signal is further transformed from time domain to frequency domain, specifically:

[0020] The Fourier transform of the discrete frequency points is performed on each group of simulated experimental wave field signals with a size of N*N t to obtain frequency domain sound field signals at different frequency points;

[0021] where the frequency points are equally spaced points within the -6dB bandwidth of the center frequency of the linear array ultrasonic transducer, and the total number is n f .

[0022] Further, in step S2, the simulation wave field is established, and the initial sound velocity model based on the simulation wave field is set, specifically:

[0023] Two linear array ultrasonic transducers with the same position and parameters as in step S1 are set in the simulation wave field.

[0024] The initial sound velocity model is an initial slowness distribution model s(x,z) representing the reciprocal of sound velocity, where x is the horizontal axis of the two-dimensional wave field, and z is the wave field vertical axis, i.e., the propagation direction.

[0025] Further, in step S3, based on the current sound speed model, forward the simulation wave field, extrapolate the incremental forward wave field signal along the propagation direction of the ultrasonic pulse signal, record all the simulation wave field signals after completing a round of emission-reception of the ultrasonic pulse signal, specifically:

[0026] The two linear array ultrasonic transducers are transmitted and received in the same way as in step S1, and the forward process of the simulation wave field is started based on the current sound speed model, and the emission frequency of the excitation source of the forward wave field is the frequency point used for the Fourier transform;

[0027] The forward process of the linear array ultrasonic transducer transceiver each time adopts the angular spectrum method, and the wave field extrapolation is performed along the ultrasonic propagation direction z direction to obtain the frequency domain sound field at the incremental Δz, and the extrapolation formula is:

[0028]

[0029] Where x and z are two coordinates of the two-dimensional wave field, exp is the exponential function, W sim,i (x,z+Δz,f) is the frequency domain sound field with coordinates (x,z+Δz) in the two-dimensional wave field, j is the imaginary unit of the complex number, s(x,z) represents the slowness at the coordinates (x,z) in the current sound field, The average slowness at the propagation depth z, F and F -1 Respectively represent the Fourier transform and inverse transform operation, The average slowness at the propagation depth z, k x is the transverse spatial frequency of the frequency domain wave field, A series of frequency points of the Fourier transform and inverse transform operation, after completing a round of emission-reception process, record all the simulation wave field signals Indicates a group of a row b column complex number array.

[0030] Further, in step S4, the loss function is constructed according to the simulation experiment wave field signal and the simulation wave field signal, and the conjugate gradient method is used to solve the loss function to obtain the optimal solution, and the current sound speed model is updated based on the iterative search direction of the conjugate gradient method, specifically:

[0031] The imaging process of full-wave inversion is the optimization and update process of model parameters, until the parameter value of the real model is approximated, and the loss function of inversion is constructed as follows:

[0032]

[0033] Where, is the frequency domain simulation experiment wave field signal, Wsim,i (s, f) is the simulated wavefield signal, N is the number of array elements;

[0034] solving the loss function by using the conjugate gradient method to obtain an optimal solution wherein the conjugate gradient method can be expressed as:

[0035]

[0036] wherein d k+1 represents the (k+1)th iteration search direction, d k represents the kth iteration search direction, L(s k +1) is a loss function obtained based on the (k+1)th iteration slow model s k +1, L(s k ) is a loss function obtained based on the kth iteration slow model s k , T represents a matrix transposition, is a gradient operator, and β k is a momentum coefficient of the current search direction, based on the iteration search direction of the conjugate gradient method, the current sound speed model is updated.

[0037] A full-wave inversion based ultrasonic cancellous bone imaging system for performing the full-wave inversion based ultrasonic cancellous bone imaging method as described above, comprising:

[0038] an analog experimental wavefield signal acquisition module, configured to set a to-be-imaged bone phantom model surrounding soft tissue, emit an ultrasonic pulse signal to the to-be-imaged bone phantom model, record the ultrasonic pulse signal representing a sound field passing through the to-be-imaged bone phantom model, and sample the ultrasonic pulse signal into an analog experimental wavefield signal;

[0039] a simulated wavefield establishment module, configured to establish a simulated wavefield and set an initial sound speed model based on the simulated wavefield, and take the initial sound speed model as a current sound speed model;

[0040] a simulated wavefield forward modeling module, configured to perform forward modeling on the simulated wavefield based on the current sound speed model, perform incremental extrapolation of a forward wavefield signal along a propagation direction of the ultrasonic pulse signal, and record all simulated wavefield signals after completing one round of emission-reception of the ultrasonic pulse signal;

[0041] a loss function construction module, configured to construct a loss function according to the analog experimental wavefield signal and the simulated wavefield signal, solve the loss function by using a conjugate gradient method to obtain an optimal solution, and update the current sound speed model based on an iteration search direction of the conjugate gradient method;

[0042] A convergence judgment module is configured to judge whether the loss function converges, output the current sound velocity model as a final sound velocity model if the loss function converges, complete the full-wave inversion process, or repeat the simulation wave field forward modeling module and the loss function construction module until the loss function converges.

[0043] A computer device comprises a memory and one or more processors, the memory stores computer code, and the computer code is executed by the one or more processors to make the one or more processors execute the method as described above.

[0044] A computer readable storage medium stores computer code, and the computer code is executed to execute the method as described above.

[0045] Compared with the prior art, the present application has the following beneficial effects:

[0046] By using the frequency domain full-wave inversion method to image the bone, the data quantity is significantly smaller than the time domain signal due to the use of discrete frequency point frequency domain signal, thereby greatly improving the calculation efficiency; the simulation sound field is reconstructed based on the full-wave inversion theory to match the actual wave field, which can display the microstructure of the bone and realize high-resolution imaging of the bone. BRIEF DESCRIPTION OF DRAWINGS

[0047] Figure 1 It is a whole flow chart of the present application based on full-wave inversion ultrasonic cancellous bone imaging method;

[0048] Figure 2 It is a specific example diagram of the bone phantom model to be imaged of the present application;

[0049] Figure 3 It is a position diagram of the linear array ultrasonic transducer and the model of the present application;

[0050] Figure 4 It is an imaging result diagram of the present application;

[0051] Figure 5 It is the sound velocity inversion effect of the imaging result and the phantom model at z=0 depth of the present application;

[0052] Figure 6 It is a whole structure diagram of the present application based on full-wave inversion ultrasonic cancellous bone imaging system. DETAILED DESCRIPTION

[0053] In order to make the purposes, technical solutions and advantages of the embodiments of the present application clearer, the technical solutions in the embodiments of the present application will be described clearly and completely below with reference to the drawings in the embodiments of the present application. Obviously, the described embodiments are some but not all of the embodiments of the present application. Based on the embodiments in the present application, all other embodiments obtained by those of ordinary skill in the art without creative efforts belong to the scope of the present application.

[0054] Those skilled in the art can understand that the singular forms "a", "an" and "the" used herein include plural forms, unless specifically stated otherwise. It should be further understood that the use of the term "comprise" in the specification of the present application means that the features, integers, steps, operations, elements and / or components exist, but does not exclude the presence or addition of one or more other features, integers, steps, operations, elements, components and / or groups thereof.

[0055] The ultrasonic cancellous bone imaging method based on full-wave inversion of the present application updates the simulation model by establishing a simulation wave field, comparing the received information of the simulation wave field with the experimental information, and then updating the simulation model to approximate the real model. Compared with the time-domain full-wave inversion method, the frequency-domain full-wave inversion method proposed by the present application only uses partial frequency point full-wave field information, which not only significantly reduces the data information storage amount, but also greatly saves the calculation time due to the characteristics of parallel computing of all frequency point wave field information.

[0056] Briefly, the inversion process of the present application is as follows: set a known bone phantom model to be imaged with a sound speed related parameter (such as the sound speed of the cancellous bone region is set to 2800 m / s; the sound speed of the surrounding soft tissue is set to 1525 m / s), collect a group of simulated experimental wave field signals, then set another group of initial sound speed models in the simulation wave field (which can be arbitrary, and does not need to be too large. For example, the initial sound speed model can be set to have a sound speed of 1525 m / s), and the simulation wave field obtains a different group of simulation wave field signals using the initial sound speed model. There is an error between the simulated experimental wave field signals and the simulation wave field signals. Based on this error, the sound speed model in the simulation wave field is updated, and finally the simulated experimental wave field signals and the simulation wave field signals are basically consistent, that is, the sound speed model in the simulation wave field is almost the same as the known sound speed model set before.

[0057] The following is described by specific embodiments:

[0058] First embodiment

[0059] As shown in Figure 1 The present embodiment provides an ultrasonic cancellous bone imaging method based on full-wave inversion, comprising the following steps:

[0060] S1: Set up a bone phantom model to be imaged, including soft tissues around it, emit an ultrasonic pulse signal to the bone phantom model to be imaged, record the ultrasonic pulse signal representing the sound field passing through the bone phantom model to be imaged, and sample the ultrasonic pulse signal as a simulated experimental wave field signal.

[0061] First, a bone phantom model to be imaged, including surrounding soft tissue, needs to be set up. In this embodiment, such as... Figure 2 As shown, a specific example of a bone phantom model to be imaged is provided. Figure 2 In the example, the bone phantom model to be imaged is cancellous bone, surrounded by soft tissue. The relevant parameters of the example bone phantom model are: the sound velocity in the cancellous bone region is set to 2800 m / s; the sound velocity in the surrounding soft tissue is set to 1525 m / s. Simulation experimental wavefield signals are acquired using a bone phantom model with known parameters.

[0062] Secondly, after setting up the bone phantom model to be imaged, a linear array ultrasonic transducer capable of emitting ultrasonic pulse signals needs to be configured. The specific configuration method is as follows: Figure 3 As shown, two identical linear array ultrasound transducers are arranged parallel to each other. Each linear array ultrasound transducer contains N array elements. The bone phantom model to be imaged is positioned between the two linear array ultrasound transducers. All array elements on one side of the linear array ultrasound transducer sequentially transmit ultrasound pulse signals with a fixed center frequency to the bone phantom model to be imaged, while simultaneously recording the ultrasound pulse signals received on the other side of the linear array ultrasound transducer. In this embodiment, each linear array ultrasound transducer has 128 array elements. All array elements of the transmitting end of the linear array ultrasound transducer sequentially transmit ultrasound pulse signals with a center frequency of 0.5 MHz, and the center distance between adjacent array elements is 0.3 mm. The sampling rate of the receiving end of the linear array ultrasound transducer is 25 MHz.

[0063] Next, the ultrasonic pulse signal needs to be sampled into the simulated experimental wavefield signal. Specifically, the ultrasonic pulse signal is sampled into N groups of N*N using the k-space pseudospectral method. t The simulated experimental wave field signal, wherein N tis the number of samples. The process of obtaining the simulated experimental wavefield signal can be achieved through the k-space pseudospectral method. In the k-space pseudospectral method, a discrete grid representing the wavefield needs to be defined first. This grid is uniformly sampled in k-space, including both frequency and wavenumber directions. Frequency represents the frequency components of the wavefield signal, while wavenumber represents the variation of the wavefield signal in space. Next, the wave equation is converted from the time domain to the frequency domain by using the Fourier transform or other related transforms. This will make the wave equation into algebraic equations in the frequency domain. In the frequency domain, these algebraic equations can be solved using discrete numerical methods such as the difference method. By numerically solving these algebraic equations, the wavefield response at different frequencies and wavenumbers can be calculated. By inverse transform (such as the inverse Fourier transform), the wavefield response can be converted from the frequency domain back to the time domain to obtain the simulated experimental wavefield signal. The advantage of the k-space pseudospectral method in simulating the experimental wavefield signal is that it can effectively handle the frequency and wavenumber distribution of the wavefield, especially suitable for the simulation of complex media and complex wave phenomena. It provides an efficient and accurate calculation method, which can provide detailed frequency distribution and spatial distribution information of the wavefield signal.

[0064] Finally, the time-domain to frequency-domain transformation of the simulated experimental wavefield signal needs to be completed, specifically: performing Fourier transform of the discrete frequency points on each group of simulated experimental wavefield signals with a size of N*N t to obtain the frequency-domain sound field signals at different frequency points; wherein the frequency points are equally spaced points within the -6dB bandwidth of the center frequency of the linear array ultrasonic transducer, and the total number is n f .

[0065] S2: Establish a simulation wavefield and set an initial sound speed model based on the simulation wavefield, taking the initial sound speed model as the current sound speed model.

[0066] Specifically, two linear array ultrasonic transducers with the same position and parameters as in step S1 are set in the simulation wavefield; the initial sound speed model is a distribution model s(x,z) of the initial slowness representing the reciprocal of sound speed, where x is the horizontal axis of the two-dimensional wavefield and z is the wavefield vertical axis, i.e. the propagation direction.

[0067] For example, the initial sound speed model can be set with a sound speed of 1525m / s. The part of the bone (2800m / s) will gradually add to the sound speed as the loss function converges.

[0068] S3: Based on the current sound speed model, forward the simulation wavefield, and extrapolate the incremental wavefield signal along the propagation direction of the ultrasonic pulse signal. After one round of emission-reception of the ultrasonic pulse signal, record all the simulation wavefield signals.

[0069] Specifically, in the present embodiment, the two linear array ultrasonic transducers are transmitted and received in the same way as in step S1, and the forward process of the simulated wave field is started based on the current sound speed model, and the transmission frequency of the excitation source of the forward wave field is each frequency point used in the Fourier transform;

[0070] The forward process of the linear array ultrasonic transducer transceiver each time uses the angular spectrum method, and the wave field extrapolation is carried out along the ultrasonic propagation direction z direction to obtain the frequency domain sound field at the increment Δz, and the extrapolation formula is:

[0071]

[0072] Where x and z are the two coordinates of the two-dimensional wave field, exp is the exponential function, W sim,i (x,z+Δz,f) is the frequency domain sound field with coordinates (x,z+Δz) in the two-dimensional wave field, j is the imaginary unit of the complex number, s(x,z) represents the slowness at the coordinates (x,z) in the current sound field, The average slowness at the propagation depth z, F and F -1 respectively represent the Fourier transform and inverse transform operation, is the average slowness at the propagation depth z, k x is the transverse spatial frequency of the frequency domain wave field, is a series of frequency points of the Fourier transform and inverse transform operation. After completing a round of transmission-reception process, all the simulated wave field signals represent a complex number array of a group of a rows and b columns.

[0073] Where, the angular spectrum method (Angular Spectrum Method) is a commonly used sound wave propagation simulation and reconstruction method, which is suitable for analyzing and simulating the propagation process of sound waves in space. It is based on the principle of Fourier transform, which can handle complex sound wave scenes and propagation media. The basic idea of angular spectrum method is to simulate the propagation and interference effect of sound waves by representing the sound wave field as a superposition of a series of plane waves. It uses the frequency and wave vector relationship (angular spectrum relationship) of sound wave propagation to convert the propagation process of sound wave field in space into phase and amplitude calculation in frequency domain.

[0074] The following are the basic steps of the angular spectrum method:

[0075] (1) Determine the initial wave field: define the initial wave field, including the position, waveform and frequency of the sound source and other parameters.

[0076] (2) Fourier transform: Fourier transform the initial wave field to get the frequency spectrum representation in the wave number domain.

[0077] (3) Apply the propagation function: In the wave number domain, apply the propagation function to calculate the propagation phase for each frequency based on the frequency and wave vector relationship.

[0078] (4) Inverse Fourier transform: Apply the inverse Fourier transform to the propagation function to convert the propagation phase in the wave number domain back to the time domain.

[0079] (5) Repeat steps (3) and (4): Repeat steps (3) and (4) to simulate the propagation process of sound waves step by step according to the required propagation distance and resolution.

[0080] (6) Reconstruct the wave field: By superimposing the wave field at each propagation distance, the final sound field distribution is obtained.

[0081] The angular spectrum method is suitable for simulating the propagation and interference process of sound waves, especially in complex acoustic scenes and media. It can handle phenomena such as sound wave diffraction, scattering, transmission, etc., and has important significance for applications such as acoustic imaging, acoustic sensing, and sound wave diffraction. It should be noted that when simulating sound wave propagation, the requirements of spatial sampling and computing resources, as well as the selection of frequency sampling and propagation distance, should be considered to obtain accurate and reliable simulation results.

[0082] S4: Construct a loss function based on the simulated experimental wave field signal and the simulated wave field signal, and solve the loss function using the conjugate gradient method to obtain the optimal solution, and update the current sound speed model based on the iterative search direction of the conjugate gradient method.

[0083] The imaging process of full-wave inversion is the optimization and update process of model parameters, until the parameter value of the true model is approximated, and the loss function of inversion is constructed as follows:

[0084]

[0085] where, is the simulated experimental wave field signal in the frequency domain, W sim,i (s, f) is the simulated wave field signal, and N is the number of elements.

[0086] Solve the loss function using the conjugate gradient method to obtain the optimal solution where the conjugate gradient method can be represented as:

[0087]

[0088] where d k+1 represents the (k+1)th iteration search direction, d k represents the kth iteration search direction, L(s k +1) is the loss function based on the (k+1)th iteration slow model s k +1, and L(sk ) is a slow model based on the kth iteration s k The resulting loss function, T represents the matrix transposition, is a gradient operator, β k is the momentum coefficient of the current search direction, based on the iterative search direction of the conjugate gradient method, updating the current sound speed model.

[0089] S5: Determine whether the loss function converges, if it converges, output the current sound speed model as the final sound speed model, complete the full-wave inversion process, otherwise repeat steps S3 and S4 until the loss function converges.

[0090] Figure 4 is the imaging result obtained by the model simulation based on the embodiment. As Figure 4 shown, the outer boundary of cancellous bone can be clearly inverted, and the overall size and shape are highly close to the real model; the internal microstructure is also visualized, and the bone microstructure such as trabecula and bone micro-pore inside the cancellous bone can be clearly observed. For Figure 5 the imaging area shown, the method can achieve an imaging result with a root mean square error of 202.91 m / s and a relative error of 7.37% compared with the real model, which shows that the algorithm can clearly and accurately reflect the bone microstructure and its distribution.

[0091] In the embodiment, the full-wave inversion based ultrasonic cancellous bone imaging method of the application can achieve high-resolution and accurate imaging of complex hard tissue models such as bone and its microstructure. Not only does it overcome the difficulty of imaging caused by the large difference in sound speed between soft and hard tissues, but the high-resolution feature of the full-wave inversion technology itself also enables the visualization of bone microstructure. In addition, the use of frequency wave field for full-wave inversion update can significantly reduce the data volume and greatly improve the computational efficiency, speeding up the imaging process.

[0092] The bone imaging method based on frequency domain full-wave inversion can achieve high-precision bone imaging, and this method has high imaging resolution, smaller data storage and processing volume, and can complete the inversion calculation process faster. The use of angular spectrum method for frequency domain multi-frequency point wave field extrapolation greatly improves the operation efficiency and speeds up the imaging process. Without the need for a priori sound speed model of cancellous bone, not only can the external contour be imaged, but also the internal microstructure details can be displayed.

[0093] Second embodiment

[0094] As Figure 6 shown, the embodiment provides a full-wave inversion based ultrasonic cancellous bone imaging system for performing the full-wave inversion based ultrasonic cancellous bone imaging method as in the first embodiment, comprising:

[0095] The simulation experiment wave field signal acquisition module 1 is used for setting a bone phantom model to be imaged including surrounding soft tissues, emitting an ultrasonic pulse signal to the bone phantom model to be imaged, recording the ultrasonic pulse signal representing a sound field passing through the bone phantom model to be imaged, and sampling the ultrasonic pulse signal into a simulation experiment wave field signal.

[0096] The simulation wave field establishment module 2 is used for establishing a simulation wave field, and setting an initial sound velocity model based on the simulation wave field, taking the initial sound velocity model as a current sound velocity model.

[0097] The simulation wave field forward modeling module 3 is used for performing forward modeling on the simulation wave field based on the current sound velocity model, performing incremental extrapolation of the forward wave field signal along the propagation direction of the ultrasonic pulse signal, recording all simulation wave field signals after completing a round of emission-reception of the ultrasonic pulse signal.

[0098] The loss function construction module 4 is used for constructing a loss function according to the simulation experiment wave field signal and the simulation wave field signal, and solving the loss function by using a conjugate gradient method to obtain an optimal solution, and updating the current sound velocity model based on the iterative search direction of the conjugate gradient method.

[0099] The convergence judgment module 5 is used for judging whether the loss function converges or not, outputting the current sound velocity model as a final sound velocity model if the loss function converges, completing the full-wave inversion process, otherwise repeating the simulation wave field forward modeling module and the loss function construction module until the loss function converges.

[0100] A computer readable storage medium stores computer code, when the computer code is executed, the above method is executed. Those skilled in the art can understand that all or part of the steps of the above method can be completed by a program instructing related hardware, and the program can be stored in a computer readable storage medium, and the storage medium can include a read only memory (ROM), a random access memory (RAM), a magnetic disk or an optical disk.

[0101] The above only describes the preferred embodiments of the present application, and the protection scope of the present application is not limited to the above embodiments. Any technical solution falling within the concept of the present application belongs to the protection scope of the present application. It should be noted that, for ordinary skilled in the art, some improvements and refinements without departing from the principles of the present application are also considered as the protection scope of the present application.

[0102] Any technical features in the above-described embodiments can be combined, and for the sake of brevity, not all possible combinations are described, however, any combination of the technical features is considered to be within the scope of the present disclosure.

[0103] It should be noted that the above-described embodiments can be freely combined as needed. The above-described only is the preferred embodiments of the present application, it should be noted that for ordinary skilled in the art, without departing from the principles of the present application, can make a number of improvements and refinements, these improvements and refinements should be considered as the scope of protection of the present application.

Claims

1. A method for ultrasound imaging of cancellous bone based on full-wave inversion, characterized in that, Includes the following steps: S1: Set up a bone phantom model to be imaged, surrounded by soft tissue. Use two linear array ultrasonic transducers arranged opposite each other to emit ultrasonic pulse signals to the bone phantom model to be imaged. Record the ultrasonic pulse signals that pass through the bone phantom model to be imaged, representing the sound field. Sample the ultrasonic pulse signals as simulated experimental wave field signals. Use Fourier transform to complete the transformation of the experimental wave field signals from the time domain to the frequency domain. S2: Establish a simulated wave field and set an initial sound velocity model based on the simulated wave field. Use the initial sound velocity model as the current sound velocity model. Specifically, set two linear array ultrasonic transducers with the same position and parameters as in step S1 in the simulated wave field. The initial sound velocity model is a distribution model s(x,z) representing the initial slowness of the reciprocal of the sound velocity, where x is the horizontal axis of the two-dimensional wave field and z is the vertical axis of the wave field, i.e., the propagation direction. S3: Based on the current sound speed model, perform forward modeling on the simulated wavefield, extrapolate the incremental forward modeled wavefield signal along the propagation direction of the ultrasonic pulse signal, and after completing one round of ultrasonic pulse signal transmission-reception, record all simulated wavefield signals as follows: The two linear array ultrasonic transducers transmit and receive in the same manner as in step S1, and start the forward modeling process of the simulated wave field based on the current sound velocity model. The transmission frequency of the excitation source of the forward modeling wave field is the frequency point used when performing the Fourier transform. Each forward modeling process of the linear array ultrasonic transducer is performed using the angular spectrum method. Wavefield extrapolation is performed along the ultrasonic propagation direction z to obtain the frequency domain sound field at the increment Δz. The extrapolation formula is as follows: Where x and z are the two coordinates of the two-dimensional wave field, exp is the exponential function, and W sim,i (x,z+Δz,f) represents the frequency domain sound field with coordinates (x,z+Δz) in a two-dimensional wave field, where j is the imaginary unit of the complex number, and s(x,z) refers to the slowness at coordinates (x,z) in the current sound field. The average slowness at propagation depth z, F and F -1 The Fourier transform and inverse transform operations are respectively represented by k. x It is the transverse spatial frequency of the frequency domain wave field. After completing one round of transmit-receive process for a series of frequency points of the Fourier transform and inverse transform operations, all the simulated wavefield signals are recorded. Let N represent a complex array of a rows and b columns, where N is a number of N*N integers. t The number of groups of the simulated experimental wavefield signals, N t Let n be the number of samples. f The total number of frequency points with equal intervals within a -6dB bandwidth that serves as the transmission center frequency of the linear ultrasonic transducer. S4: Construct a loss function based on the simulated experimental wave field signal and the simulated wave field signal, and solve the loss function using the conjugate gradient method to obtain the optimal solution. Update the current sound speed model based on the iterative search direction of the conjugate gradient method. S5: Determine whether the loss function has converged. If it has converged, output the current sound speed model as the final sound speed model to complete the full-wave inversion process. Otherwise, repeat steps S3 and S4 until the loss function converges.

2. The ultrasound cancellous bone imaging method based on full-wave inversion according to claim 1, characterized in that, In step S1, the ultrasonic pulse signal is emitted to the bone phantom model to be imaged, and the ultrasonic pulse signal representing the sound field passing through the bone phantom model to be imaged is recorded, specifically: Two identical linear array ultrasound transducers are arranged in parallel relative to each other. Each linear array ultrasound transducer contains N array elements. The bone phantom model to be imaged is placed between the two linear array ultrasound transducers. All the array elements on one side of the linear array ultrasound transducer sequentially transmit ultrasound pulse signals with a fixed center frequency to the bone phantom model to be imaged, while simultaneously recording the ultrasound pulse signals received on the other side of the linear array ultrasound transducer.

3. The ultrasound cancellous bone imaging method based on full-wave inversion according to claim 2, characterized in that, In step S1, the ultrasonic pulse signal is sampled as the simulated experimental wave field signal, specifically as follows: The ultrasonic pulse signal was sampled into N groups of size N*N using the k-space pseudospectral method. t The simulated experimental wave field signal, wherein N t The number of samples.

4. The ultrasound cancellous bone imaging method based on full-wave inversion according to claim 3, characterized in that, Step S1 further includes: transforming the simulated experimental wavefield signal from the time domain to the frequency domain, specifically: The size of each group is N*N t The simulated experimental wave field signal was subjected to discrete frequency Fourier transform to obtain frequency domain sound field signals at different frequency points; Wherein, the frequency points are the equally spaced points within the -6dB bandwidth of the transmission center frequency of the linear array ultrasonic transducer, and the total number is n. f .

5. The ultrasound cancellous bone imaging method based on full-wave inversion according to claim 1, characterized in that, In step S4, the loss function is constructed based on the simulated experimental wave field signal and the simulated wave field signal, and the loss function is solved using the conjugate gradient method to obtain the optimal solution. The current sound speed model is updated based on the iterative search direction of the conjugate gradient method, specifically as follows: The full-wave inversion imaging process is a process of optimizing and updating the model parameters until they approximate the parameter values ​​of the true model. The loss function for inversion is constructed as follows: in, For the simulated experimental wavefield signal in the frequency domain, W sim,i (s,f) represents the simulated wave field signal, and N is the number of array elements; The optimal solution is obtained by solving the loss function using the conjugate gradient method. The conjugate gradient method can be expressed as: Where, d k+1 d represents the search direction in the (k+1)th iteration. k L(s) represents the search direction in the k-th iteration. k+1 ) is based on the slow model s of the (k+1)th iteration. k+1 The resulting loss function, L(s) k ) is based on the slow model s of the k-th iteration. k The resulting loss function, where T represents the matrix rank transformation. For the gradient operator, β k It is the momentum coefficient of the current search direction. Based on the iterative search direction of the conjugate gradient method, the current sound speed model is updated.

6. A full-wave inversion-based ultrasound cancellous bone imaging system for performing the full-wave inversion-based ultrasound cancellous bone imaging method as described in any one of claims 1-5, characterized in that, include: The simulated experimental wave field signal acquisition module is used to set up a bone phantom model to be imaged, including soft tissues, emit ultrasonic pulse signals to the bone phantom model to be imaged, record the ultrasonic pulse signals representing the sound field passing through the bone phantom model to be imaged, and sample the ultrasonic pulse signals as a simulated experimental wave field signal. The simulation wave field establishment module is used to establish a simulation wave field and set an initial sound speed model based on the simulation wave field, and use the initial sound speed model as the current sound speed model. The simulated wave field forward modeling module is used to perform forward modeling of the simulated wave field based on the current sound speed model, extrapolate the incremental forward modeled wave field signal along the propagation direction of the ultrasonic pulse signal, and record all simulated wave field signals after completing one round of ultrasonic pulse signal transmission-reception. The loss function construction module is used to construct a loss function based on the simulated experimental wave field signal and the simulated wave field signal, solve the loss function using the conjugate gradient method to obtain the optimal solution, and update the current sound speed model based on the iterative search direction of the conjugate gradient method. The convergence judgment module is used to determine whether the loss function has converged. If it has converged, the current sound velocity model is output as the final sound velocity model to complete the full wave inversion process. Otherwise, the simulation wave field forward modeling module and the loss function construction module are repeated until the loss function converges.

7. A computer device comprising a memory and one or more processors, the memory storing computer code that, when executed by the one or more processors, causes the one or more processors to perform the method as described in any one of claims 1 to 5.

8. A computer-readable storage medium storing computer code, wherein when the computer code is executed, the method of any one of claims 1 to 5 is performed.

Images