Cross-correlation correction full-waveform inversion ultrasonic tomography method based on source coding

Abstract

The invention provides a cross-correlation correction full-waveform inversion ultrasonic tomography method based on source coding. The method comprises the following steps of: 1, acquiring an actually measured sound pressure signal, setting an initial sound velocity model and the number of iterations of each stage, and initializing the number of iterations; step 2, calculating travel time difference between a prediction signal dm (xn, t) and a measured signal gm (xn, t) in the first iteration of each stage based on the initial sound velocity model in the step 1 and the number of iterations of each stage, generating an intermediate signal of the stage based on travel time, and not updating the travel time difference in the stage except the first iteration; 3, coding the intermediate signal and the prediction signal based on the intermediate signal and the prediction signal generated in the step 2; 4, calculating an adjoint wave field and gradient based on the intermediate signal and the prediction signal coded in the step 3, determining the step length of sound velocity model iteration, and completing model updating; and 5, obtaining a final model. According to the method, the stability of a traditional FWI-based ultrasonic tomography algorithm is improved while the calculated amount is not remarkably increased.

Description

Technical Field

[0001] The invention belongs to the technical field of ultrasonic tomography, and in particular relates to a source coding-based cross-correlation correction full waveform inversion ultrasonic tomography method. Background Art

[0002] Ultrasound tomography (USCT) is an emerging medical imaging technique, particularly in breast cancer detection. The goal of USCT is to obtain information about the acoustic properties of the interior of a tissue from ultrasound data measured at its boundaries. Empirical studies have shown that the acoustic properties of diseased tissue, such as the speed of sound, are different from those of normal tissue. USCT can quantitatively reconstruct an image of the acoustic properties of tissue, thereby providing valuable diagnostic information. While other medical imaging modalities, such as mammography and MRI, are often good enough in many cases, USCT has several advantages worth exploring, including the absence of ionizing radiation, lower cost, and providing more acoustic performance.

[0003] Due to the similarities between medical USCT and seismic tomography, FWI, which is traditionally used in seismology, has gained increasing attention in recent years in medical USCT imaging. By updating the acoustic model, the predicted signal is matched to the measured signal using the gradient descent technique, thereby reducing the objective function. Since the FWI method takes into account the high-order diffraction effects, the reconstructed image has a higher spatial resolution than the traditional ray-based methods. Therefore, FWI is a promising method for quantitative soft tissue imaging.

[0004] However, in clinical applications, these methods have difficulty converging to the true model and require a high computational burden. FWI uses the adjoint method to calculate the gradient of the objective function, which requires two simulations of the wave field for each transmission in the measured data set to calculate the gradient. In order to increase the speed of the algorithm, the execution time for solving the wave equation can be significantly reduced by using the acceleration provided by the parallel processing of graphics processing units (GPUs). Then, with the implementation of numerical methods, the computational efficiency of FWI has been greatly improved. Moreover, by reducing the number of solutions to the wave equation, the computational complexity can be reduced. In the field of geology, researchers have proposed the use of source coding to reduce the number of solutions to the wave equation, which effectively speeds up the efficiency of FWI. However, source coding also has its limitations. The coded wave field causes the stacked signal to contain multiple arrival events from all emission locations, which increases the risk of convergence to the local minimum.

[0005] On the other hand, in order to adopt FWI in clinical systems, it is necessary to overcome the non-convexity of FWI. The fidelity term of the objective function is usually the least squares loss, which is locally convex when the time shift between the predicted signal and the measured signal is less than half a cycle. This problem is called cycle jump, so an accurate prior model is required to avoid being trapped in the local minimum. When the prior model is too different from the actual model, the stability of the FWI technology will be seriously damaged. In order to improve the stability of the FWI algorithm, the present invention is based on an acoustic model with no mass density change and no absorption, and proposes an ultrasound tomography sound velocity reconstruction algorithm based on source coded cross-correlation corrected FWI, wherein the cross-correlation correction improves the stability of FWI, and the introduction of source coding reduces the computational complexity of the algorithm. Summary of the invention

[0006] The present invention provides a source coding-based cross-correlation correction full waveform inversion ultrasonic tomography method, which improves the stability of a traditional FWI-based ultrasonic tomography algorithm without significantly increasing the amount of calculation.

[0007] The present invention is achieved through the following technical solutions:

[0008] A source-coded cross-correlation corrected full waveform inversion ultrasonic tomography method, the method comprising the following steps:

[0009] Step 1: Collect the measured sound pressure signal, set the initial sound velocity model and the number of iterations in each stage, and initialize the number of iterations;

[0010] Step 2: Based on the initial sound speed model of step 1 and the number of iterations in each stage, the predicted signal d is calculated at the first iteration of each stage. m (x n ,t) and the measured signal g m (x n ,t), and use the travel time difference to generate the intermediate signal of this stage. In this stage, the travel time difference and the intermediate signal are not updated except for the first iteration;

[0011] Step 3: Based on the intermediate signal and the prediction signal generated in step 2, encode the intermediate signal and the prediction signal;

[0012] Step 4: Based on the intermediate signal and predicted signal encoded in step 3, calculate the accompanying wave field and gradient, determine the step size of the sound velocity model iteration, and complete the model update;

[0013] Step 5: Get the final model.

[0014] Furthermore, the step 1 is specifically that the ultrasonic signal is generated by M circular ultrasonic transducers with uniform spacing, so the data acquisition process of USCT includes M ultrasonic pulse transmissions. In the mth transmission, the acoustic pulse transmitted by the mth transducer unit is recorded as s m (x, t), the sound pressure wave field generated in the plane is recorded as p(x, t). Without considering the change and absorption of mass and density, p(x, t) satisfies the sound wave equation:

[0015]

[0016] In formula (1), x∈Ω m represents the imaging area, t∈[0,T] represents the sampling time, and p m (x, t) represents the sound pressure signal in space generated at position x after the mth transducer emits a pulse, and c(x) is the sound velocity distribution model of the tissue to be measured.

[0017] Furthermore, the step 1 specifically includes the following steps:

[0018] Step 1.1: M ultrasonic transducers emit ultrasonic pulse signals in sequence and collect the generated sound pressure signals g m (x n ,t), where m=1,2,3,…,M; n=1,2,3,…,M;

[0019] Step 1.2: Set the initial sound velocity model c0(x) of the object to be tested and the number of iterations k included in each iteration stage max ;

[0020] Step 1.3: The total number of iterations is reset to zero, that is, n=0, and the phased iteration is recorded as 1, that is, k=1.

[0021] Furthermore, the step 2 is specifically as follows:

[0022] The sound pressure field p under the current sound velocity model is solved by formula (1): m (x,t), at the nth ultrasonic transducer position x n The sound pressure at is recorded as the prediction signal d m (x n ,t), which satisfies formula (2):

[0023] d m (x n ,t)=p m (x n ,t)(2)

[0024] In the mth transmission, the measured signal and the predicted signal are at the spatial position x n The travel time difference at is the maximum value of the cross-correlation, which can be expressed as:

[0025]

[0026] Travel time difference Δt s And the ultrasonic emission frequency f should satisfy the relationship (4):

[0027]

[0028] Therefore, the time threshold δ is selected to satisfy The threshold is the travel time difference between the intermediate signal and the predicted signal, which avoids periodic jumps; thus, the time shift factor Δ corresponding to the intermediate signal generated by the measured signal when the ultrasonic wave is emitted for the mth time can be obtained. m (x n ) satisfies the following formula:

[0029]

[0030] According to formula (5), the intermediate signal can be obtained as g m (x n ,t+Δ m (x n )).

[0031] Furthermore, the step 2 specifically includes the following steps:

[0032] Step 2.1: If k≠1, jump to step 3, otherwise continue to perform the following operations;

[0033] Step 2.2: Use equation (1) to solve the sound pressure field p under the current sound velocity model: m (x,t),

[0034] Step 2.3: Calculate the predicted signal d by equation (2) m (x n ,t);

[0035] Step 2.4: Calculate the measured signal and the predicted signal at the spatial position x by using equation (3) n The travel time difference τ m (x n );

[0036] Step 2.5: Select the time threshold δ and calculate the time shift factor Δ of the intermediate signal relative to the measured signal by equation (5): m (x n ), thereby obtaining the intermediate signal g m (x n ,t+Δ m (x n )).

[0037] Furthermore, the step 3 is specifically as follows:

[0038] The source of emission is represented by:

[0039]

[0040] Based on the linearity of the wave equation with respect to the source term, the predicted acoustic pressure field is expressed as:

[0041]

[0042] Therefore, the prediction signal is encoded as:

[0043]

[0044] And, the intermediate signal is encoded as:

[0045]

[0046] Furthermore, the step 3 specifically includes the following steps:

[0047] Step 3.1: Randomly generate the encoding vector w through Rademacher distribution;

[0048] Step 3.2: Calculate the encoded prediction signal d through equation (8) w (x n ,t);

[0049] Step 3.3: Use equation (9) to convert the intermediate signal g m (x n ,t+Δ m (x n )) is encoded to obtain the encoded intermediate signal

[0050] Furthermore, the step 4 is specifically that the FWI nonlinear numerical optimization problem is expressed as:

[0051]

[0052] Where F(c), R(c) and α represent the data fidelity term, regularization term and regularization parameter respectively. The data fidelity and regularization term together constitute the objective function J(c) = F(c) + αR(c);

[0053] The data fidelity term F(c) is defined as the squared L2 norm of the deviation between the measured signal and the predicted signal. After introducing source coding and intermediate signals, it is reformulated as the expectation of a random quantity:

[0054]

[0055] Wherein, E represents the expectation operator of the residual with respect to the randomly encoded intermediate signal and the randomly encoded prediction signal;

[0056] The regularization term R(c) is expressed as:

[0057]

[0058] Among them, c i,j Represents the speed of sound c at pixel (i, j). To prevent the denominator from being zero when taking the derivative, a small positive number ∈ is added.

[0059] Furthermore, the step 4 specifically includes the following steps:

[0060] Step 4.1: Calculate the adjoint wave field q by the following formula: w (x,t):

[0061]

[0062] Step 4.2: Calculate the gradient matrix D(c):

[0063] D(c)=D F (c)+αD R (c)#(14)

[0064] Among them, D F (c) and D R (c) corresponds to the functional gradient of the data fidelity and regularization term in the objective function with respect to the sound velocity model c(x), which is given by the following formula: F (c) and D R (c) is calculated as follows:

[0065]

[0066] Formula (16) is the value of D at pixel point (i, j): R (c) the value of;

[0067] Step 4.3: Use the steepest descent method to determine the step size λ of the model iteration through straight line search i , the model update is realized by the following formula:

[0068] c i+1 =c i -λ i D(c)#(17)

[0069] Wherein, the subscript i refers to the i-th iteration;

[0070] Step 4.4: Total number of iterations n = n + 1;

[0071] Step 4.5: Determine whether the iteration end condition is met. If so, execute step 5. Otherwise, continue to execute the following steps;

[0072] Step 4.6: The number of stage iterations k is less than kmax , then let k = k + 1, otherwise, let k = 1;

[0073] Step 4.7: Return to step 2.

[0074] Furthermore, the step 5 specifically includes taking the latest updated sound field data after the iteration as the final sound field data, and reconstructing the image of the tissue to be tested with the sound field data.

[0075] The beneficial effects of the present invention are:

[0076] The present invention determines the time shift of the intermediate signal relative to the measured signal according to the travel time difference between the measured signal and the predicted signal, so that the travel time difference between the intermediate signal and the predicted signal is limited to a certain range. The intermediate signal is used as the target of inversion, so that the predicted signal is updated to the intermediate signal through iteration, and the travel time difference between the measured signal and the predicted signal is reduced. FWI will update the sound velocity distribution in the correct direction, cleverly avoiding the periodic jump phenomenon of traditional FWI and improving the stability of FWI imaging.

[0077] The present invention introduces an intermittent update strategy to avoid the problem that the prediction signal cannot be encoded to reduce the amount of calculation due to the need to calculate the travel time difference. The intermittent update strategy divides multiple iterations into a stage, and only calculates the travel time difference in the first iteration of each stage to complete the update of the intermediate signal. After the first iteration of this stage, FWI acceleration with source coding is introduced to reduce the amount of calculation. The introduction of the intermittent update strategy improves the stability of FWI without adding a high computational burden. BRIEF DESCRIPTION OF THE DRAWINGS

[0078] Figure 1 It is a flow chart of the method of the present invention.

[0079] Figure 2 It is a schematic diagram of the USCT data acquisition system in the simulation of the present invention.

[0080] Figure 3 Schematic diagram of the numerical phantom of the present invention.

[0081] Figure 4 It is a schematic diagram of the initial sound velocity model of the present invention.

[0082] Figure 5 It is a schematic diagram of the sound velocity distribution of the simulated medium of the present invention.

[0083] Figure 6 These are the images reconstructed after 10, 50, 100, and 300 FWI-SE iterations of the present invention, where (a) is FWI-SE iteration 10 times, (b) is FWI-SE iteration 50 times, (c) is FWI-SE iteration 100 times, and (d) is FWI-SE iteration 300 times.

[0084] Figure 7 These are the images reconstructed after 10, 50, 100, and 300 CCAFWI-SE iterations, where (a) is CCAFWI-SE iteration 10 times, (b) is CCAFWI-SE iteration 50 times, (c) is CCAFWI-SE iteration 100 times, and (d) is CCAFWI-SE iteration 300 times.

[0085] Figure 8 It is the gradient of FWI-SE and CCAFWI-SE in the first iteration of the present invention, where (a) is the gradient distribution of the first iteration of FWI-SE, and (b) is the gradient distribution of the first iteration of CCAFWI-SE.

[0086] Fig. 9 It is the root mean square error RMSE and similarity SSIM curve of 300 iterations in the simulation of the present invention, wherein (a) is the curve of RMSE changing with the number of iterations, and (b) is the curve of SSIM changing with the number of iterations. DETAILED DESCRIPTION

[0087] The technical solutions in the embodiments of the present invention will be described clearly and completely below in conjunction with the drawings in the embodiments of the present invention. Obviously, the described embodiments are only part of the embodiments of the present invention, not all of the embodiments. Based on the embodiments of the present invention, all other embodiments obtained by ordinary technicians in this field without creative work are within the scope of protection of the present invention.

[0088] Embodiment 1

[0089] This embodiment provides a source-coded cross-correlation corrected full waveform inversion ultrasound tomography method, in which an ultrasound transducer transmits an ultrasound signal to the tissue to be measured, and after the ultrasound signal propagates in the tissue to be measured, the transducer obtains the sound field data, i.e., the measured signal. The predicted signal is obtained by solving the acoustic wave equation, and the travel time difference between the predicted signal and the measured signal can be calculated using the maximum value of the cross-correlation between the predicted signal and the measured signal. Based on the size of the travel time difference, an intermediate signal after the measured signal is time-shifted is generated, and the travel time difference between the intermediate signal and the predicted signal is limited to a small threshold, and the intermediate signal is used as the target of waveform inversion. In addition, an intermittent stage-by-stage update of the intermediate signal is adopted, that is, a certain number of iterations are regarded as a stage, and each stage only calculates the travel time difference and updates the intermediate signal in the first iteration, and in each iteration after the current stage, the FWI acceleration program based on source coding is used, and then after multiple stage-by-stage iterations, the predicted signal can be close to the measured signal, thereby obtaining the sound velocity distribution of the tissue to be measured, i.e., ultrasound tomography.

[0090] The method comprises the following steps,

[0091] Step 1: Collect the measured sound pressure signal, set the initial sound velocity model and the number of iterations in each stage, and initialize the number of iterations;

[0092] Specifically, the ultrasonic signal is generated by M circular ultrasonic transducers with uniform spacing. Therefore, the data acquisition process of USCT includes M ultrasonic pulse transmissions. In the mth transmission, the acoustic pulse transmitted by the mth transducer unit is denoted as s m (x, t), the sound pressure wave field generated in the plane is recorded as p(x, t). Without considering the change and absorption of mass and density, p(x, t) satisfies the sound wave equation:

[0093]

[0094] In formula (1), x∈Ω m represents the imaging area, t∈[0,T] represents the sampling time, and p m (x, t) represents the acoustic pressure signal in space generated at x after the mth transducer transmits a pulse, and c(x) is the sound velocity distribution model of the tissue to be measured. In the actual data acquisition process, the acoustic pressure signal cannot be obtained at any position in the imaging area, but can only be obtained at a fixed position of the transducer. In the mth transmission, the pressure signal measured by the nth ultrasonic transducer is recorded as g m (x n ,t).

[0095] The step 1 specifically comprises the following steps:

[0096] Step 1.1: M ultrasonic transducers emit ultrasonic pulse signals in sequence and collect the generated sound pressure signals g m (x n ,t), where m=1,2,3,…,M; n=1,2,3,…,M;

[0097] Step 1.2: Set the initial sound velocity model c0(x) of the object to be tested and the number of iterations k included in each iteration stage max ;

[0098] Step 1.3: The total number of iterations is reset to zero, that is, n=0, and the phased iteration is recorded as 1, that is, k=1.

[0099] Step 2: Based on the initial sound speed model of step 1 and the number of iterations in each stage, the predicted signal d is calculated at the first iteration of each stage. m (x n ,t) and the measured signal g m (x n ,t), and use the travel time difference to generate the intermediate signal of this stage. In this stage, the travel time difference and the intermediate signal are not updated except for the first iteration;

[0100] The step 2 is specifically as follows:

[0101] The sound pressure field p under the current sound velocity model is solved by formula (1): m (x,t), at the nth ultrasonic transducer position x n The sound pressure at is recorded as the prediction signal d m (x n ,t), which satisfies formula (2):

[0102] d m (x n ,t)=p m (x n ,t)(2)

[0103] In the mth transmission, the measured signal and the predicted signal are at the spatial position x n The travel time difference at is the maximum value of the cross-correlation, which can be expressed as:

[0104]

[0105] To avoid cycle jump, the travel time difference Δt s And the ultrasonic emission frequency f should satisfy the relationship (4):

[0106]

[0107] Therefore, the time threshold δ is selected to satisfy The threshold is the travel time difference between the intermediate signal and the predicted signal, which avoids periodic jumps; thus, the time shift factor Δ corresponding to the intermediate signal generated by the measured signal when the ultrasonic wave is emitted for the mth time can be obtained. m (x n ) satisfies the following formula:

[0108]

[0109] According to formula (5), the intermediate signal can be obtained as g m (x n ,t+Δ m (x n )).

[0110] The step 2 specifically includes the following steps:

[0111] Step 2.1: If k≠1, jump to step 3, otherwise continue to perform the following operations;

[0112] Step 2.2: Use equation (1) to solve the sound pressure field p under the current sound velocity model: m (x,t),

[0113] Step 2.3: Calculate the predicted signal d by equation (2)m (x n ,t);

[0114] Step 2.4: Calculate the measured signal and the predicted signal at the spatial position x by using equation (3) n The travel time difference τ m (x n );

[0115] Step 2.5: Select the time threshold δ and calculate the time shift factor Δ of the intermediate signal relative to the measured signal by equation (5): m (x n ), thereby obtaining the intermediate signal g m (x n ,t+Δ m (x n )).

[0116] Step 3: Based on the intermediate signal and the prediction signal generated in step 2, encode the intermediate signal and the prediction signal;

[0117] Specifically, step 3 introduces source coding to reduce the computational cost of FWI. In each iteration of source coding FWI, the ultrasonic signals emitted by M ultrasonic units are encoded by a random coding vector w at the same time, which has a zero mean and a unit covariance matrix. In this case, the source of the emission is expressed as:

[0118]

[0119] Based on the linearity of the wave equation with respect to the source term, the predicted acoustic pressure field is expressed as:

[0120]

[0121] Therefore, the prediction signal is encoded as:

[0122]

[0123] And, the intermediate signal is encoded as:

[0124]

[0125] The step 3 specifically includes the following steps:

[0126] Step 3.1: Randomly generate the encoding vector w through Rademacher distribution;

[0127] Step 3.2: Calculate the encoded prediction signal d through equation (8) w (x n ,t);

[0128] Step 3.3: Use equation (9) to convert the intermediate signal gm (x n ,t+Δ m (x n )) is encoded to obtain the encoded intermediate signal

[0129] Step 4: Based on the intermediate signal and predicted signal encoded in step 3, calculate the accompanying wave field and gradient, determine the step size of the sound velocity model iteration, and complete the model update;

[0130] Specifically, step 4 is that FWI expresses image reconstruction as a large-scale non-convex optimization problem, which is essentially fitting the waveform of the predicted signal in the forward model of sound propagation with the waveform of the measured signal. The FWI nonlinear numerical optimization problem is expressed as:

[0131]

[0132] Where F(c), R(c) and α represent the data fidelity term, regularization term and regularization parameter respectively. The data fidelity and regularization term together constitute the objective function J(c) = F(c) + αR(c);

[0133] The data fidelity term F(c) is defined as the squared L2 norm of the deviation between the measured signal and the predicted signal. After introducing source coding and intermediate signals, it is reformulated as the expectation of a random quantity:

[0134]

[0135] Where E represents the expectation operator of the residual with respect to the randomly encoded intermediate signal and the randomly encoded prediction signal; its optimization is a stochastic optimization problem that can be solved by various stochastic gradient-based algorithms;

[0136] The regularization term R(c) is expressed as:

[0137]

[0138] Among them, c i,j Represents the speed of sound c at pixel (i, j). To prevent the denominator from being zero when taking the derivative, a small positive number ∈ is added.

[0139] The step 4 specifically comprises the following steps:

[0140] Step 4.1: Calculate the adjoint wave field q by the following formula: w (x,t):

[0141]

[0142] Step 4.2: Calculate the gradient matrix D(c):

[0143] D(c)=DF (c)+αD R (c)#(14)

[0144] Among them, D F (c) and D R (c) corresponds to the functional gradient of the data fidelity and regularization term in the objective function with respect to the sound velocity model c(x), which is given by the following formula: F (c) and D R (c) is calculated as follows:

[0145]

[0146] Formula (16) is the value of D at pixel point (i, j): R (c) the value of;

[0147] Step 4.3: Use the steepest descent method to determine the step size λ of the model iteration through straight line search i , the model update is realized by the following formula:

[0148] c i+1 =c i -λ i D(c)#(17)

[0149] Wherein, the subscript i refers to the i-th iteration;

[0150] Step 4.4: Total number of iterations n = n + 1;

[0151] Step 4.5: Determine whether the iteration end condition is met. If so, execute step 5. Otherwise, continue to execute the following steps;

[0152] Step 4.6: The number of stage iterations k is less than k max , then let k = k + 1, otherwise, let k = 1;

[0153] Step 4.7: Return to step 2.

[0154] Step 5: Get the final model;

[0155] Specifically, the step 5 is to use the latest updated sound field data after the iteration as the final sound field data, and to reconstruct the image of the tissue to be tested with the sound field data.

[0156] Embodiment 2

[0157] This embodiment performs simulation according to the technical solution provided in the first embodiment. Specifically, the example of the present invention completes the reconstruction of the sound velocity by means of simulation. The numerical phantom is acoustically reconstructed through two-dimensional simulation, and the annular distribution shape of the measurement array element is used to simulate the actual USCT system. In the simulation, a ring transducer with 128 uniformly distributed units is used to generate the source signal. These elements emit acoustic pulses in sequence and scan the numerical phantom in sequence. The data acquisition consists of 128 transmissions. In each transmission, all the remaining transducers should record the acoustic wave field signal. The schematic diagram is shown in the figure. Figure 2 shown.

[0158] In the simulation, the acoustic excitation pulse emitted by the working component is a sine function with a center frequency of 1.5MHz, modulated by a Gaussian kernel, and a standard deviation of 0.05μs. The time domain expression of the excitation pulse is:

[0159]

[0160] In the formula, f c =1.5MHz,σ=0.05μs,t c =3μs.

[0161] The parameters of the simulation process are shown in Table 1. The numerical phantom used in the simulation is a contrast-enhanced cone-beam breast CT image with a tumor near the center of the breast. The horizontal and vertical dimensions of the numerical phantom are 46 mm and 43 mm, respectively. Acoustic properties such as sound attenuation are not considered in the numerical phantom. Figure 3 shown.

[0162] Table 1 Simulation parameters

[0163]

[0164] The above is the end of the simulation setting.

[0165] Execute step 1: Use the sequence emission method to collect the sound pressure signal of the above simulation environment to obtain the measured signal. Initialize the number of iterations n to 0, and the number of staged iterations k to 1. Figure 4is the initial sound velocity model in the simulation, which is a uniform model with a sound velocity of 1500m / s, the same as the background. The black dot indicates the position of the transducer element. In the process of signal acquisition, data beyond the range of the array element emission angle is generally discarded. On the one hand, the actual signal acquisition equipment has a limited focusing angle of the transmitting array element, and the received signal contains a lot of noise; on the other hand, the signal received outside the angle does not contain the phantom information and is redundant data for data fitting. Therefore, in order to be closer to the real physical experiment, for the measured signal, only the signal received at the 64 receivers opposite the transmitter is used. However, when using source coded FWI, it is required that the signal cannot be discarded. Therefore, when performing data fitting, the signal measured in a pure water environment is used to supplement the discarded signal. This will not affect the final reconstructed image, because when using UCT to image breast tissue, it will be immersed in water, and the discarded signal does not contain the phantom information. It is in line with the actual situation to use the sound pressure data in a pure water environment to supplement the discarded signal.

[0166] Execute step 2: Determine whether the current iteration is the first iteration of this iteration phase, that is, whether k is 1. If so, it is necessary to complete the prediction signal d by solving the M-th acoustic wave equation. m (x n ,t) to calculate the predicted signal d m (x n ,t) and the measured signal g m (x n The time difference of the intermediate signal relative to the measured signal is determined by the time difference. m (x n ), and then generate an intermediate signal. If it is not the first iteration of this iteration stage, directly execute step 3.

[0167] Execute step 3: randomly generate a coding vector w through Rademacher distribution, encode the prediction signal and the intermediate signal, and obtain the encoded prediction signal d w (x n ,t) and intermediate signal

[0168] Execute step 4: After introducing source coding, calculate the accompanying wave field and gradient, use the steepest descent method, determine the step size of the sound velocity model iteration through straight line search, complete the model update, and increase the number of iterations n by one. Determine whether the total number of iterations reaches the preset maximum value of 300. If it reaches the maximum value, jump to step 5; if it does not reach the maximum value, update k according to the size of the staged iteration number k, and jump to step 2.

[0169] Execute step five: use the latest updated sound field data after the iteration as the final sound field data, and reconstruct the image of the simulation phantom with this sound field data.

[0170] By iteratively executing the above steps, the distribution of the sound velocity model to be measured can be obtained by inverting the acoustic wave equation and the measured sound pressure signal, thereby obtaining ultrasonic tomography.

[0171] Figure 5 To simulate the medium sound velocity, the black dots indicate the locations of the transducer elements. To eliminate the influence of artifacts, the imaging area is selected as a white rectangular area. Figure 6 These are the reconstruction results of the traditional source coded full waveform inversion (FWI-SE) after 10, 50, 100, and 300 iterations. Figure 7 The reconstruction results after 10, 50, 100 and 300 iterations using source-coded cross-correlation corrected full waveform inversion (CCAFWI-SE) are shown. Obviously, the reconstruction results of FWI-SE are very different from the original model, while CCAFWI-SE reconstructs the sound velocity model better.

[0172] The significant difference between the initial model and the true model causes the distance between the predicted signal and the measured signal to exceed half a cycle. The traditional FWI-SE will produce a cycle jump, which can be observed from the gradient at the first iteration, such as Figure 8 shown. Figure 8 In the figure, the red rectangle corresponds to the location of the tumor, and a positive gradient in the area enclosed by the rectangle indicates that the sound velocity model will decrease, while a negative gradient indicates that the sound velocity model will increase. Figure 8 In (a), the area near the tumor is mainly positive, indicating that the sound velocity in this area should be reduced. However, in the simulation setting, the sound velocity of the tumor is greater than the sound velocity of the background, and the sound velocity in this area should be increased to iterate towards the correct model. This shows that due to the cycle jump, the iteration direction of FWI-SE is opposite to the required iteration direction, and it falls into a local minimum, and cannot reconstruct the tumor area.

[0173] In contrast, the reconstruction results of CCAFWI-SE are not distorted. The gradient of the first iteration is as follows Figure 8 (b). It can be seen that the gradient is significantly different from the gradient of FWI-SE. The gradient calculated using the intermediate signal is mainly negative near the tumor, which means that there is an increasing trend in the sound velocity model. The intermediate signal is regenerated at the beginning of each stage. Then the entire stage uses source encoding to reduce the computational cost. After completing 10 stage-by-stage iterations, the sound velocity model is close enough to the true model, and there is no need to use cross-correlation to adjust the gradient. Therefore, CCAFWI-SE can correctly update the sound velocity model at a lower computational cost.

[0174] The curves of the root mean square error (RMSE) and similarity (SSIM) of 300 iterations in the simulation with the number of iterations are as follows: Fig. 9 (a) and (b). Fig. 9In (a), the RMSE using FWI-SE does not decrease due to periodic jumps. Finally, the RMSE is 28.47 m / s, which is larger than the RMSE of the initial model. This shows that the gradient direction is not toward the true model, causing the reconstruction result to reach a local minimum. CCAFWI-SE uses cross-correlation calculations to adjust the gradient. Therefore, the RMSE using CCAFWI-SE can be continuously reduced and finally reaches 4.79 m / s. The results show that the sound speed model of CCAFWI-SE can be updated toward the true model.

[0175] The SSIM curves also differ significantly: after 300 iterations, the SSIM using FWI-SE is 0.73, and the SSIM using CCAFWI-SE is 0.91. Similarly, at the beginning of the inversion, both SSIM curves show a downward trend due to source coding. Subsequently, the SSIM curve of FWI-SE shows a trend of first rising and then falling, while the SSIM curve of CCAFWI-SE continues to rise.

[0176] It can be seen that the cross-correlation corrected FWI imaging method based on source coding in the present invention not only introduces the intermediate signal generated by the cross-correlation between the measured signal and the predicted signal, and takes the intermediate signal as the inversion target to ensure the correctness of the iteration direction, but also creatively designs a way to update the intermediate signal in stages, introduces source coding, reduces the calculation cost, and improves the real-time performance of the traditional FWI algorithm, which is of great significance in ensuring the stability of the USCT system for tissue imaging.

Claims

1. A source-coded cross-correlation corrected full waveform inversion ultrasonic tomography method, characterized in that: The method comprises the following steps, Step 1: Collect the measured sound pressure signal, set the initial sound velocity model and the number of iterations in each stage, and initialize the number of iterations; Step 2: Based on the initial sound speed model of step 1 and the number of iterations in each stage, the predicted signal d is calculated at the first iteration of each stage. m (x n ,t) and the measured signal g m (x n ,t), and use the travel time difference to generate the intermediate signal of this stage. In this stage, the travel time difference and the intermediate signal are not updated except for the first iteration; Step 3: Based on the intermediate signal and the prediction signal generated in step 2, encode the intermediate signal and the prediction signal; Step 4: Based on the intermediate signal and predicted signal encoded in step 3, calculate the accompanying wave field and gradient, determine the step size of the sound velocity model iteration, and complete the model update; Step 5: Get the final model.

2. The method of full waveform inversion ultrasonic tomography based on source coding cross-correlation correction according to claim 1, characterized in that: Specifically, the ultrasonic signal is generated by M circular ultrasonic transducers with uniform spacing. Therefore, the data acquisition process of USCT includes M ultrasonic pulse transmissions. In the mth transmission, the acoustic pulse transmitted by the mth transducer unit is denoted as s m (x, t), the sound pressure wave field generated in the plane is recorded as p(x, t). Without considering the change and absorption of mass and density, p(x, t) satisfies the sound wave equation: In formula (1), x∈Ω m represents the imaging area, t∈[0,T] represents the sampling time, and p m (x, t) represents the sound pressure signal in space generated at position x after the mth transducer emits a pulse, and c(x) is the sound velocity distribution model of the tissue to be measured.

3. The method of full waveform inversion ultrasonic tomography based on source coding cross-correlation correction according to claim 2, characterized in that: The step 1 specifically comprises the following steps: Step 1.1: M ultrasonic transducers emit ultrasonic pulse signals in sequence and collect the generated sound pressure signals g m (x n ,t), where m=1,2,3,…,M; n=1,2,3,…,M; Step 1.2: Set the initial sound velocity model c0(x) of the object to be tested and the number of iterations k included in each iteration stage max ; Step 1.3: The total number of iterations is reset to zero, that is, n=0, and the phased iteration is recorded as 1, that is, k=1.

4. The method of full waveform inversion ultrasonic tomography based on source coding cross-correlation correction according to claim 2, characterized in that: The step 2 is specifically as follows: The sound pressure field p under the current sound velocity model is solved by formula (1): m (x,t), at the nth ultrasonic transducer position x n The sound pressure at is recorded as the prediction signal d m (x n ,t), which satisfies formula (2): d m (x n ,t)=p m (x n ,t) (2) In the mth transmission, the measured signal and the predicted signal are at the spatial position x n The travel time difference at is the maximum value of the cross-correlation, which can be expressed as: Travel time difference Δt s And the ultrasonic emission frequency f should satisfy the relationship (4): Therefore, the time threshold δ is selected to satisfy The threshold is the travel time difference between the intermediate signal and the predicted signal, which avoids periodic jumps. Thus, the time shift factor Δ corresponding to the intermediate signal generated by the measured signal when the ultrasonic wave is emitted for the mth time can be obtained. m (x n ) satisfies the following formula: According to formula (5), the intermediate signal can be obtained as g m (x n ,t+Δ m (x n )).

5. The method of full waveform inversion ultrasonic tomography based on source coding cross-correlation correction according to claim 4, characterized in that: The step 2 specifically includes the following steps: Step 2.1: If k≠1, jump to step 3, otherwise continue to perform the following operations; Step 2.2: Use equation (1) to solve the sound pressure field p under the current sound velocity model: m (x,t), Step 2.3: Calculate the predicted signal d by equation (2) m (x n ,t); Step 2.4: Calculate the measured signal and the predicted signal at the spatial position x by using equation (3) n The travel time difference τ m (x n ); Step 2.5: Select the time threshold δ and calculate the time shift factor Δ of the intermediate signal relative to the measured signal by equation (5): m (x n ), thereby obtaining the intermediate signal g m (x n ,t+Δ m (x n )).

6. The method of full waveform inversion ultrasonic tomography based on source coding cross-correlation correction according to claim 4, characterized in that: The step 3 is specifically as follows: The source of emission is represented by: Based on the linearity of the wave equation with respect to the source term, the predicted acoustic pressure field is expressed as: Therefore, the prediction signal is encoded as: And, the intermediate signal is encoded as:

7. The method of full waveform inversion ultrasonic tomography based on source coding cross-correlation correction according to claim 6, characterized in that: The step 3 specifically comprises the following steps: Step 3.1: Randomly generate the encoding vector w through Rademacher distribution; Step 3.2: Calculate the encoded prediction signal d through equation (8) w (x n ,t); Step 3.3: Use equation (9) to convert the intermediate signal g m (x n ,t+Δ m (x n )) is encoded to obtain the encoded intermediate signal 8. The method of full waveform inversion ultrasonic tomography based on source coding cross-correlation correction according to claim 6, characterized in that: Specifically, the step 4 is that the FWI nonlinear numerical optimization problem is expressed as: Where F(c), R(c) and α represent the data fidelity term, regularization term and regularization parameter respectively. The data fidelity and regularization term together constitute the objective function J(c) = F(c) + αR(c); The data fidelity term F(c) is defined as the squared L2 norm of the deviation between the measured signal and the predicted signal. After introducing source coding and intermediate signals, it is reformulated as the expectation of a random quantity: Wherein, E represents the expectation operator of the residual with respect to the randomly encoded intermediate signal and the randomly encoded prediction signal; The regularization term R(c) is expressed as: Among them, c i,j Represents the speed of sound c at pixel (i, j). To prevent the denominator from being zero when taking the derivative, a small positive number ∈ is added.

9. The method of full waveform inversion ultrasonic tomography based on source coding cross-correlation correction according to claim 8, characterized in that: The step 4 specifically comprises the following steps: Step 4.1: Calculate the adjoint wave field q by the following formula: w (x,t): Step 4.2: Calculate the gradient matrix D(c): D(c)=D F (c)+αD R (c)#(14) Among them, D F (c) and D R (c) corresponds to the functional gradient of the data fidelity and regularization term in the objective function with respect to the sound velocity model c(x), which is given by the following formula: F (c) and D R (c) is calculated as follows: Formula (16) is the value of D at pixel point (i, j): R (c) the value of; Step 4.3: Use the steepest descent method to determine the step size λ of the model iteration through straight line search i , the model update is realized by the following formula: c i+1 =c i -λ i D(c)#(17) Wherein, the subscript i refers to the i-th iteration; Step 4.4: Total number of iterations n = n + 1; Step 4.5: Determine whether the iteration end condition is met. If so, execute step 5. Otherwise, continue to execute the following steps; Step 4.6: The number of stage iterations k is less than k max , then let k = k + 1, otherwise, let k = 1; Step 4.7: Return to step 2.

10. The method of full waveform inversion ultrasonic tomography based on source coding and cross-correlation correction according to claim 6, characterized in that: Specifically, the step 5 is to use the latest updated sound field data after the iteration as the final sound field data, and to reconstruct the image of the tissue to be tested with the sound field data.

Images