Ultrasonic tomography method

Abstract

The present application provides an ultrasonic tomography method, comprising: S1, using a ring array to perform full-matrix data acquisition on a target to be measured to determine observation data of N frequency points; S2, selecting observation data of the i-th frequency point; S3, encoding the observation data of two adjacent sound sources at the i-th frequency point to obtain first data; S4, encoding the two adjacent sound sources in an inversion process at the i-th frequency point based on the Green's function to obtain second data; S5, determining a target function based on the first data and the second data; the target function is used to quantify the difference between the observation data of the i-th frequency point and the synthetic data predicted based on the initial sound velocity parameter image; S6, iterating the initial sound velocity parameter image m0 to determine the sound velocity parameter image of the i-th frequency point; S7, using the sound velocity parameter image as the initial sound velocity parameter image at the i+1-th frequency point, executing the S2-S6 process, traversing the observation data of the N frequency points, and obtaining an ultrasonic tomography image.

Description

Technical Field

[0001] The present application relates to the field of medical ultrasound technology, and in particular to an ultrasonic tomography method. Background Art

[0002] Ultrasound computed tomography (USCT) is a non-invasive, radiation-free, and low-cost detection technology that uses a surrounding ultrasound transducer array to reveal the internal structure and acoustic properties of an object. USCT is of great value in the diagnosis of early-stage lesions. This technology can be broadly divided into two categories: imaging methods based on ray assumptions, such as delay and sum (DAS) and time of flight tomography (TOFT), which are suitable for rapid imaging; and full waveform inversion (FWI) algorithms based on wave theory, which provide high-resolution quantitative images. However, these methods are computationally intensive, have slow imaging speeds, and have not yet widely met clinical needs.

[0003] To improve imaging efficiency, FWI can be inverted in the frequency domain, a technique called Frequency-Domain Full Waveform Inversion (FD-FWI). This employs a multi-scale inversion strategy, starting with low-frequency data and gradually transitioning to high-frequency data. This effectively reduces computational resource consumption and makes image reconstruction possible on a personal computer. Although this strategy has significantly reduced computational requirements, for annular arrays consisting of hundreds or even thousands of transducers, the algorithm still faces challenges in meeting the image reconstruction time requirements of current clinical applications. The primary challenge lies in the large number of wavefield numerical solutions required to be calculated due to the numerous transmitting transducers. Summary of the Invention

[0004] To address the specific needs of data acquisition and FWI, this application proposes an ultrasound tomography method from a frequency domain perspective. This method employs static source coding for an annular array and implements multi-source coding throughout the annular array's transmission and final inversion processes, comprehensively improving the imaging efficiency of frequency-domain FWI-based USCT. This method optimizes the imaging process and enhances the quality and reliability of medical imaging.

[0005] In a first aspect, an embodiment of the present application provides an ultrasonic tomography method, the method comprising: S1, in an immersion environment, using a multi-element ultrasonic annular array to perform full-matrix data acquisition on a target to be measured therein, and determining the frequency domain observation data of the position of each array element, the frequency domain observation data including the observation data of N frequency points; S2, selecting the observation data of the i-th frequency point from the observation data of the N frequency points; the observation data of the N frequency points are arranged in order from low frequency to high frequency, the frequency of the i-th frequency point is less than the frequency of the i+1-th frequency point; i is a natural number, i+1≤N; S3, determining a source coding strategy, the source coding strategy includes determining the frequency of any two adjacent sound sources emitted together time interval and spatial distribution position; encoding the observation data of two adjacent sound sources at the i-th frequency point according to the source coding strategy to obtain the first data; S4, encoding and predicting the two adjacent sound sources in the inversion process at the i-th frequency point based on the Green function and the source coding strategy to obtain the second data; the Green function is solved based on the initial sound speed parameter image m0; S5, determining the objective function based on the first data and the second data; the objective function is used to quantify the difference between the observation data of the i-th frequency point and the synthetic data predicted based on the initial sound speed parameter image m0; S6, using the full waveform inversion algorithm to iteratively update the initial sound speed parameter image m0 to determine the sound speed parameter image m0 of the i-th frequency point i ; S7, the speed of sound parameter image m i As the initial sound speed parameter image m0 at the i+1th frequency point, let i=i+1, execute S2-S6 process, traverse the observation data of N frequency points, and calculate the sound speed parameter image m0 according to the sound speed parameter image m0. n The ultrasonic tomographic image of the target to be measured is obtained by inversion.

[0006] In some embodiments, a multi-element ultrasonic annular array is used to perform full matrix data acquisition on the target to be measured, and the frequency domain observation data at the location of each array element is determined; the method includes: one array element in the annular array transmits an ultrasonic signal to the target to be measured, and all array elements receive the echo signal, and all array elements are traversed in sequence, and there is no overlap between the ultrasonic signals transmitted by each array element; and obtaining the observation sound pressure signal D at the location of each array element. obs (t,x r ,x s ), the amount of data is N t ×N r ×N s ; Use discrete Fourier transform to obtain the frequency domain observation data set d obs (f,x r ,x s ); where N t is the number of time samples, N r is the number of receiving array elements and N s is the number of transmitting array elements; x sRepresents the position coordinates of the emission point and x r is the position coordinate of the receiving point; t represents the signal sequence corresponding to the time domain {t i , i=1,2,…N t}, the unit is seconds; f represents the signal sequence corresponding to the frequency domain {f i , i=1,2,…N t}, unit is Hertz.

[0007] In some embodiments, determining the source coding strategy includes: determining the time interval between any two sound sources to be: τ = (2n + 1) / (4f i ); where n is an integer; determine the spatial distribution position information of the two sound sources are adjacent and The encoding function is determined as:

[0008]

[0009] where ω i is the corresponding frequency f i The angular frequency below.

[0010] In some embodiments, the observation data of two adjacent sound sources at the i-th frequency point are encoded according to the source coding strategy to obtain the first data; comprising: encoding the observation data of two adjacent sound sources at the i-th frequency point obtained in the full matrix acquisition mode using the source coding strategy, and superimposing the encoded data to obtain the first data as

[0011]

[0012] where s j The sequence {s j , j=1,2…N s}, s j ′ represents the sequence of sound sources after coding and superposition {s j ′, j′=1,2…N s ′},N s′ Indicates the total number of multi-source coding transmissions. When two sound sources are superimposed for coding

[0013] In some embodiments, encoding and predicting two adjacent sound sources in the inversion process at the i-th frequency point based on the Green function and the source coding strategy to obtain the second data includes: and The second data is encoded using the source coding strategy and transmitted sequentially at the same frequency with a time interval τ. as follows:

[0014]

[0015] Where G represents the Green's function, which is solved by the Helmholtz equation. The corresponding expression is as follows:

[0016]

[0017] Where ω represents the angular frequency, m represents the sound velocity parameter image of the medium, Δ represents the Laplace operator, and x is the coordinate of the calculation domain. is the Dirac delta function.

[0018] In some embodiments, determining the target function based on the first data and the second data includes: determining a single frequency point f i The objective function C se (m) is calculated as follows:

[0019]

[0020] Among them, the L2 norm is used to characterize the first data With the second data difference.

[0021] In some embodiments, iteratively updating the sound speed parameter image m in the objective function includes iteratively updating the initial sound speed parameter image m0 using a local gradient optimization algorithm as follows:

[0022]

[0023] Where k represents the number of iterations; α k is the iteration step size; is the descending direction;

[0024] Gradient optimized using previous iteration information Get the current descent direction

[0025] The gradient is calculated using the adjoint state method as follows:

[0026]

[0027] where u se (x) and They represent the forward wavefield after source coding and the residual reverse wavefield at the receiving point in the inversion process respectively; R represents the real part; * represents the conjugate.

[0028] In a second aspect, an embodiment of the present application provides an electronic device comprising: at least one memory for storing programs; and at least one processor for executing the programs stored in the memory. When the programs stored in the memory are executed, the processor is used to execute the method provided in the first aspect.

[0029] In a third aspect, an embodiment of the present application provides a computer storage medium, in which instructions are stored. When the instructions are executed on a computer, the computer executes the method provided in the first aspect. BRIEF DESCRIPTION OF THE DRAWINGS

[0030] In order to more clearly illustrate the technical solutions of the multiple embodiments disclosed in this specification, the following will briefly introduce the drawings required for the description of the embodiments. Obviously, the drawings described below are only the multiple embodiments disclosed in this specification. For ordinary technicians in this field, other drawings can be obtained based on these drawings without any creative work.

[0031] The following is a brief introduction to the drawings required for describing the embodiments or prior art.

[0032] Figure 1 A flowchart of an ultrasonic tomography method proposed in an embodiment of the present application;

[0033] Figure 2 A schematic diagram of two-source coded short-interval transmission provided in an embodiment of the present application;

[0034] Figure 3 A top view of the annular array observation system and pork tissue provided in Example 1 of the present application;

[0035] Figure 4 This is a flowchart of the application of the ultrasonic tomography method provided in Example 1 of the present application;

[0036] Figure 5 Schematic diagram of pork tissue sound velocity reconstruction using two adjacent sources, four sources, and eight adjacent sources simultaneously transmitting under traditional FWI and no sound source coding;

[0037] Figure 6 Schematic diagram of pork tissue sound velocity reconstruction using traditional FWI and multi-source coding with two adjacent sources, four sources, and eight adjacent sources with short intervals.

[0038] Figure 7 Schematic diagram of pork tissue sound velocity reconstruction using traditional FWI and fixed-time interval multi-source coding with two adjacent sources, four sources, and eight sources with short intervals. DETAILED DESCRIPTION

[0039] In order to make the purpose, 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 below with reference to the accompanying drawings.

[0040] In the description of the embodiments of the present application, words such as "exemplary," "for example," or "for example" are used to indicate examples, illustrations, or descriptions. Any embodiment or design described as "exemplary," "for example," or "for example" in the embodiments of the present application should not be construed as being preferred or advantageous over other embodiments or designs. Rather, the use of words such as "exemplary," "for example," or "for example" is intended to present the relevant concepts in a concrete manner.

[0041] In the description of the embodiments of this application, the term "and / or" is simply a description of the association relationship between associated objects, indicating that three relationships can exist. For example, A and / or B can represent the following three situations: A exists alone, B exists alone, and A and B exist at the same time. In addition, unless otherwise specified, the term "plurality" means two or more. For example, "multiple systems" refers to two or more systems, and "multiple terminals" refers to two or more terminals.

[0042] Furthermore, the terms "first" and "second" are used for descriptive purposes only and should not be construed as indicating or implying relative importance or implicitly identifying the technical features being referred to. Thus, features specified as "first" or "second" may explicitly or implicitly include one or more of such features. The terms "include," "comprising," "having," and their variations all mean "including but not limited to," unless otherwise specifically emphasized.

[0043] In the description of the embodiments of the present application, reference is made to “some embodiments”, which describe a subset of all possible embodiments, but it can be understood that “some embodiments” may be the same subset or different subsets of all possible embodiments, and may be combined with each other without conflict.

[0044] In the description of the embodiments of the present application, the terms "first\second\third, etc." or module A, module B, module C, etc. are only used to distinguish similar objects and do not represent a specific ordering of the objects. It can be understood that the specific order or sequence can be interchanged where permitted so that the embodiments of the present application described here can be implemented in an order other than that illustrated or described here.

[0045] In the description of the embodiments of the present application, the numbers representing the steps, such as S110, S120, etc., do not necessarily mean that the steps must be executed in this manner. If permitted, the order of the previous and next steps can be interchanged, or they can be executed simultaneously.

[0046] Unless otherwise defined, all technical and scientific terms used herein have the same meaning as commonly understood by those skilled in the art to which this application pertains. The terms used herein are for the purpose of describing the embodiments of this application only and are not intended to limit this application.

[0047] In order to make the purpose, 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 below with reference to the accompanying drawings.

[0048] In the field of medical imaging, the use of source coding technology for signal acquisition and inversion imaging is crucial for improving the efficiency of practical applications. To address the computational efficiency issues associated with FWI, the present invention proposes an ultrasonic tomography method that utilizes a source coding algorithm to address the high computational complexity and resource requirements of existing technologies, thereby improving data acquisition efficiency. By implementing a source coding scheme in the frequency domain, this method promotes the application of annular array-based FWI imaging technology in clinical settings, providing a faster and more accurate diagnostic tool.

[0049] In the ultrasonic tomography method proposed in the embodiment of the present application, a source coding algorithm is implemented in the frequency domain. Considering the large memory requirement when solving the frequency domain Helmholtz equation, a method based on the convergent Bern series (CBS) is used to numerically calculate the frequency domain wave field. In order to reduce the number of solutions of the forward wave field during the FWI calculation process, a multi-source excitation algorithm based on source coding is derived through theoretical approximation to improve the computational efficiency of imaging.

[0050] In this application, the embodiment specifically focuses on reconstructing a sound velocity parameter image of the medium, using sound velocity as a key indicator for identifying different tissues. Because different media have different sound velocities, each sound velocity value on the sound velocity parameter image corresponds to a specific tissue type within the medium, allowing us to clearly distinguish and identify them.

[0051] Figure 1 This is a flow chart of an ultrasonic tomography method proposed in an embodiment of the present application. Figure 1As shown, the method includes: S1, in an immersion environment, using a ring array of multi-element ultrasonic transducers to perform full-matrix data acquisition on a target to be measured therein, and determining frequency domain observation data of the position of each array element, where the frequency domain observation data includes observation data of N frequency points; S2, selecting observation data of the i-th frequency point from the observation data of the N frequency points; the observation data of the N frequency points are arranged in order from low frequency to high frequency, and the frequency of the i-th frequency point is less than the frequency of the i+1-th frequency point; i is a natural number, i+1≤n; S3, determining a source coding strategy, where the source coding strategy includes determining the time interval and spatial distribution position of the simultaneous emission of any two adjacent sound sources; According to the source coding strategy, the observation data of two adjacent sound sources corresponding to the i-th frequency point are encoded to obtain the first data; S4, the two adjacent sound sources emitted at a time interval in the inversion process are encoded and synthesized based on the Green function at the i-th frequency point to obtain the second data; the Green function is solved based on the initial sound speed parameter image m0; S5, the objective function is determined based on the first data and the second data; the objective function is used to quantify the difference between the observation data of the i-th frequency point and the synthetic data predicted based on the initial sound speed parameter image m0; S6, the initial sound speed parameter image m0 is iteratively updated using the full waveform inversion algorithm to determine the sound speed parameter image m0 of the i-th frequency point i ; S7, the speed of sound parameter image m i As the initial sound velocity parameter image m0 of the next frequency point, let i=i+1, execute steps S2-S6, traverse the observation data of the selected N frequency points, and calculate the initial sound velocity parameter image m0 according to the sound velocity parameter image m0. n Invert the ultrasonic tomographic image of the target to be measured.

[0052] The above steps of an ultrasonic tomography method proposed in an embodiment of the present application are discussed in detail below.

[0053] In some implementations, step S1 may be implemented by the following steps:

[0054] S11, an array element in the annular array transmits an ultrasonic signal to the target to be measured, and all array elements receive the echo signal. All array elements are traversed in turn, and there is no overlap between the ultrasonic signals transmitted by each array element.

[0055] For example, in an immersion environment, the first array element in the circular array can transmit an ultrasonic signal to the target to be measured. After all other array elements receive the echo signal, the next array element transmits an ultrasonic signal to the target to be measured. All other array elements receive the echo signal, and the process traverses all array elements in turn. There is no overlap between the ultrasonic signals transmitted by each array element.

[0056] S12, obtain the observed sound pressure signal D at each array element location obs (t,x r ,x s ), the amount of data is Nt ×N r ×N s Where N t ,N r ,N s They represent the number of time samples, the number of receiving array elements, and the number of transmitting array elements respectively.

[0057] S13, use discrete Fourier transform to obtain the corresponding frequency domain observation data set d for these observation sound pressure signals. obs (f,x r ,x s ), where x s , x r Represent the spatial position coordinates of the transmitting array element and the receiving array element respectively; t represents the signal sequence corresponding to the time domain {t i , i=1,2,…N t}, the unit is seconds; f represents the signal sequence corresponding to the frequency domain {f i , i=1,2,…N t}, unit is Hertz.

[0058] In some implementations, step S2 may be implemented by the following steps:

[0059] S21, observe the data set d in the frequency domain according to the frequency band information of the transmitted signal obs (f,x r ,x s ) to select observation data of multiple frequency points.

[0060] S22, arranging the observation data of the multiple frequency points in order from low frequency to high frequency.

[0061] S23, first select the observation data d of the i-th frequency point from the observation data of multiple frequency points obs (f i ,x r ,x s ); the frequency of the i-th frequency point < the frequency of the i+1-th frequency point; i is a natural number.

[0062] For example, observation data of each frequency point with a frequency band interval of 0.3 to 1.0 MHz and a frequency interval of Δf=0.05 MHz can be selected. Based on the multi-scale inversion strategy, observation data can be selected starting from the observation data of the lowest frequency point 0.3 MHz and gradually transitioning to the observation data of the highest frequency point 1.0 MHz. The frequency band interval of each step is Δf=0.05 MHz.

[0063] After obtaining the observation data of the i-th frequency point, in order to obtain a better inversion imaging effect, it is necessary to ensure that the crosstalk term is an imaginary number to reduce the proportion of the crosstalk term in the imaging. The method provided in the embodiment of the application uses a source coding strategy in step S3 to encode the observation data of the i-th frequency point and the inverted sound source separately.

[0064] like Figure 2 As shown, for the ultrasonic transducer ring array composed of multiple array elements, theoretical approximate deduction and experimental test show that the best position of two sound sources when transmitting in a group of coded signals is adjacent distribution. Therefore, the spatial position information of the two sound sources can be determined to be adjacent x j and x j+1 The time interval between two adjacent sound sources at the same frequency is τ=(2n+1) / (4f i ), where n is an integer. This source coding strategy can avoid interference effects between sound sources and crosstalk that affect imaging.

[0065] In some embodiments, based on the above Figure 2 Given the theoretical approximate derivation and experimental test results, step S3 can be achieved by the following steps:

[0066] S31, determining a source coding strategy, where the source coding strategy includes determining the emission time interval between any two adjacent sound sources and the positions of the two sound sources in a distribution space.

[0067] Among them, the delay interval τ between any two adjacent sound sources is: τ=(2n+1) / (4f i ); where n is an integer, and the spatial distribution position information of two adjacent sound sources are respectively j and x j+1 ;

[0068] S32, determine the encoding function ε(f i )as follows:

[0069]

[0070] S33, the sound sources collected by the full matrix are grouped into two adjacent pairs, and the emission time interval τ of the two adjacent sound sources and the spatial distribution position x of the sound sources are calculated. j and x j+1 , using the encoding function ε(f i ) for frequency domain observation data d obs (f i ,x r ,x s )and Encoding, a new set of frequency domain observation data after encoding is

[0071]

[0072] where s j The sequence {s j , j=1,2…N s}, s j ′ represents the sequence of sound sources after coding combination {s j ′, j′=1,2…N s ′}, when encoding the frequency domain observation data of two adjacent sound source combinations The new observation data can be recorded as Recorded as the first data.

[0073] For example, taking the first two adjacent array elements of the full matrix acquisition as a group of sound sources, the spatial position of the distribution is The new observation data obtained by encoding the observation data using the encoding function of formula (1) is:

[0074]

[0075] Since it is impossible to accurately obtain the sound source information during the processing of actual observation data, a sound source estimation strategy can be adopted, and the corresponding source coding function can be applied to the sound source estimation process. Step S4 can be implemented by the following steps:

[0076] S41, in the inversion process, any two consecutive adjacent sound source information and The coding function and Green's function are used to jointly encode and obtain the synthesized data after the sound source is encoded. The corresponding expressions are as follows:

[0077]

[0078] The synthesized data after the sound source is encoded Denoted as the second data, where G represents the Green's function, which is solved by the Helmholtz equation, and the corresponding expression is as follows:

[0079]

[0080] Where ω represents the angular frequency, m represents the sound velocity parameter image of the target to be measured, Δ represents the Laplace operator, x is the coordinate of the calculation domain, δ(x s -x) is the Dirac delta function.

[0081] In a water environment, the water sound speed parameter image can be imaged as a uniform sound speed of 1540 m / s. The water sound speed parameter image can be used as the initial sound speed parameter image and iteratively updated in the subsequent inversion process.

[0082] S42, the convergent Bern series CBS method is used to simulate the frequency domain wave field of the Helmholtz equation (5) above, and the Green function G(f i ,x,x s ).

[0083] S43, N s The array elements are divided into N s / 2 groups of adjacent array element pairs, each array element is used as a transmitting unit, and the inverted synthetic data is transmitted at the same frequency with a delay interval τ

[0084] For example, let the first two adjacent sound sources g(f i ,x1) and g(f i ,x2) adopts the coding function of formula (1) and the Green function obtained by solving formula (5) to jointly encode, and transmits at a time interval of τ, and the synthetic data of the first two adjacent sound sources obtained by inversion is The calculation formula is as follows:

[0085]

[0086] in is the Green function, G is obtained by solving the Helmholtz equation, and the corresponding expression is as follows:

[0087]

[0088] So, for an N s The ultrasonic transducer ring array has N array elements. In the inversion process, s The ultrasonic transducer ring array composed of array elements only needs to transmit N s / 2 times, N s The array elements are divided into N s / 2 groups of adjacent array element pairs, each pair of array elements is used as a transmitting unit, which not only greatly reduces the number of wave field numerical solutions that need to be calculated, but also avoids the impact of crosstalk and interference effects between two adjacent sound sources on imaging by adding a time interval τ for delay coding.

[0089] Next, step S5 is executed to determine the objective function according to the first data and the second data.

[0090] In some embodiments, the first data may be characterized using the L2 norm. The second data in the inversion process The difference between the two frequencies is defined as a single frequency point f i The objective function C of the down-inverted sound velocity parameter image m se The formula for (m) is shown in (7):

[0091]

[0092] Objective function C se (m) is used to quantify the difference between the observed data at the i-th frequency point and the synthetic data predicted based on the initial sound speed parameter image m0, N s′ Indicates the total number of multi-source coded transmissions.

[0093] In step S6, the full waveform inversion algorithm is used to iteratively update the initial sound velocity parameter image m0, including the following steps:

[0094] S61, use the local gradient optimization algorithm to iteratively update the sound speed parameter image m0, as shown in formula (8):

[0095]

[0096] Among them, k is the number of iterations, α k The line search method can be used to obtain the iteration step size; is the descending direction;

[0097] S62, by selecting the conjugate gradient method or L-BFGS method to use the previous iteration information to optimize the current gradient To get a better descent direction

[0098] S63, the adjoint state method is used to obtain the gradient, as shown in formula (9):

[0099]

[0100] where u se (f i ,x,x s′ )and They represent the forward wavefield after source coding and the residual reverse wavefield at the receiving point in the inversion process respectively; R represents the real part; * represents the conjugate.

[0101] S64, performing multiple iterative calculations under the convergence judgment condition until the iteration stop condition of the current frequency point is met, and determining the sound speed parameter image m.

[0102] In some embodiments, the condition for stopping the iteration may be defined as reaching a set number of iterations.

[0103] For step S7, the ultrasonic tomography method provided in the embodiment of the present application adopts a multi-scale inversion strategy, starting from low-frequency point observation data and gradually transitioning to high-frequency point observation data for inversion, and using the sound velocity parameter image obtained by low-frequency inversion as the next higher-frequency initial sound velocity parameter image, repeating the process of some or all of the above S2-S6 implementation methods, and obtaining the ultrasonic tomography of the target to be measured based on the final sound velocity parameter image inversion result.

[0104] The ultrasonic tomography method provided in the embodiment of the present application addresses the problem that the computational efficiency of the existing FWI-based USCT imaging method is difficult to meet current clinical needs, and a solution is derived through theoretical approximation. This solution effectively reduces the number of forward wave field solutions during the FWI calculation process by implementing a multi-source coded emission strategy. Furthermore, a simple source coding scheme is designed so that adjacent sound sources are emitted at short intervals according to a set delay at a specific frequency, thereby effectively reducing the crosstalk effect during the inversion process. This strategy is not only comparable to the traditional FWI method (full matrix acquisition and full matrix inversion) in imaging accuracy, but also speeds up the entire imaging process and thus improves overall work efficiency. This helps to promote the acceleration of FWI imaging based on annular arrays from data acquisition to inversion.

[0105] Example 1

[0106] Figure 3 The annular array observation system and pork tissue top view provided in Example 1 of this application. Figure 3 As shown, in a water immersion environment, the block of pork belly is wrapped into a nearly hollow cylindrical shape, and a 512-element circular array with an inner diameter of 22 cm is used for data acquisition to perform ultrasonic tomography imaging of the pork tissue.

[0107] Figure 4 This is an application flow chart of the ultrasonic tomography method provided in Example 1 of the present application. Figure 4 As shown, the following steps are included:

[0108] S81, for Figure 3 The 512-element annular array with a center frequency of 0.75 MHz is used to collect full matrix observation signals of pork tissue in a water immersion environment. The block of pork belly is wrapped into a hollow cylindrical shape for data collection, and the observation sound pressure signal D at the position of the array element is obtained. obs (t,x r ,x s ), the data volume is 4955×512×512. The corresponding frequency domain observation data set d is obtained by discrete Fourier transform. obs (f,x r ,x s). Set the initial sound speed parameter image m0 to the sound speed parameter image of water, with a sound speed of 1540 m / s.

[0109] S82, execute LOOP1: select the observation data of each frequency point with a frequency interval of Δf=0.05MHz from 0.3 to 1.0MHz, and based on the multi-scale inversion strategy, gradually transition from the observation data of the lowest frequency point 0.3MHz to the observation data of the highest frequency point 1.0MHz. The frequency band interval of each observation data is Δf=0.05MHz, f=f1:Δf:f n =0.3:0.05:1.0. For the specific implementation, please refer to steps S21-S23.

[0110] S83, iteratively implementing LOOP2 at the i-th frequency point selected by LOOP1, including the following steps:

[0111] S831, using a determined source coding function, encodes the observation data of two adjacent sound sources corresponding to the frequency point, so that the observation data of one sound source is not processed, and the observation data of the other sound source is multiplied by the coding function ε(f i )=exp(-iω i / (4f i )), we get a new set of observation data as The remaining observation data are processed in the same way in sequence. During the inversion process, 512 array elements only need to be calculated 256 times. For a specific implementation, please refer to steps S31-S33.

[0112] S832, at the same frequency, two adjacent sound sources transmitted at short intervals are encoded at the corresponding frequencies so that one of the sound sources is not processed and the other sound source is multiplied by the encoding function ε(f i )=exp(-iω i / (4f i )), the convergent Bern series CBS method is used to simulate the frequency domain wave field of two adjacent source coded sound sources with interval transmission, solve the Green function, and obtain the synthetic data The remaining sound sources and observation data are processed in the same way in sequence. During the inversion process, 512 array elements only need to transmit 256 times. For a specific implementation, please refer to steps S41-S43.

[0113] S833: Determine the objective function based on the encoded observation data and the synthetic data from the forward modeling. For a specific implementation, refer to step S5. S834: Calculate the parameter gradient using the adjoint state method for the sound velocity parameter image from the forward modeling of the frequency domain wavefield based on the CBS.

[0114] S835: Use the local optimal algorithm (conjugate gradient) and line search method to obtain the update direction and step size of the sound velocity parameter image during the iteration process, and update the sound velocity parameter image until the iteration stop condition is met. For detailed implementation, please refer to steps S61-S64.

[0115] S836, perform convergence judgment, and determine whether the conditions for stopping LOOP2 iteration are met based on the results. If the conditions for stopping LOOP2 iteration are met, execute step S84; if the conditions for stopping LOOP2 iteration are not met, return to step S83.

[0116] S84, determine the next frequency point f n+1 =f n +Δf, the sound velocity parameter image obtained at the current frequency point (such as 0.3MHz) is used as the initial model of the next frequency point (such as 0.35MHz).

[0117] S85, judging whether the conditions for ending LOOP1 are met according to the current frequency point, if the conditions for stopping LOOP1 are met, proceeding to step S86, if the conditions for stopping LOOP1 are not met, returning to executing step S82.

[0118] S86, after traversing all selected frequency points, the inversion is terminated to obtain the ultrasonic tomographic imaging of the pork tissue.

[0119] The above coding strategy is also extended to the case of coded transmission of 4 and 8 adjacent sources. Under the same multi-scale inversion strategy, the reconstruction of sound velocity parameter images using full matrix transmission and adjacent short-interval transmission of two, four, and eight sources with and without delay coding is tested.

[0120] Figure 5 Schematic diagram of the reconstructed pork tissue sound velocity using two adjacent sources, four sources, and eight sources simultaneously transmitting under traditional FWI and no sound source delay coding. The title 1in1 indicates that one sound source is transmitted at a time; 2in1 indicates that two sound sources are transmitted at a time, and the others are similar. SE represents encoding, and 1 and 0 indicate whether the encoding function ε(f i )=exp(-iω i / (4f i )).

[0121] like Figure 5 The reconstruction results of pork tissue are shown (1 and 0 in the figure represent whether the encoding function ε(f i )=exp(-iω i / (4f iIn the absence of source delay coding, as the number of simultaneously transmitting sound sources increases, the reconstruction results deteriorate accordingly. Crosstalk artifacts cause the inversion to fall into a local minimum, which in turn causes the final result of the iteration to differ significantly from the actual situation.

[0122] Figure 6 Schematic diagram of pork tissue sound velocity reconstruction using traditional FWI and adjacent two-source, four-source, and eight-source short-interval transmissions with sound source delay coding.

[0123] like Figure 6 As shown in the figure, after using the encoding strategy, the phenomenon of crosstalk artifacts causing the inversion to fall into local minima is improved, the reconstruction results are closer to the actual situation, such as the sound velocity at the skin of pork tissue, and some artifacts are suppressed.

[0124] To be more practical, a multi-source coded transmission strategy is adopted during the annular array acquisition, and the delay interval of all frequencies is set to a fixed 1 / (4f0), where f0 is selected as the center frequency of the sound source.

[0125] Figure 7 Schematic diagram of pork tissue sound velocity reconstruction using short-interval transmission of two adjacent sources, four sources, and eight sources under traditional FWI and fixed delay interval τ = 1 / (4f0) sound source coding.

[0126] Figure 7 The results show that the method using a fixed delay interval is consistent with the 1 / (4f i ) method, its performance is not significantly inferior and is more in line with the needs and conditions of practical applications.

[0127] This application presents an ultrasonic tomography method based on source-coded FD-FWI (Frequency Wiring) USCT image reconstruction technology, suitable for scenarios with multi-source, short-interval coded transmissions. Compared to the sound velocity results reconstructed using traditional FWI, under the same conditions, multi-source, short-interval coded transmissions can achieve similar inversion accuracy to traditional FWI while reducing the number of forward wavefield acquisitions or calculations, thereby reducing computational complexity and accelerating the entire imaging process.

[0128] It is understood that the processor in the embodiments of the present application may be a central processing unit (CPU), or may be other general-purpose processors, digital signal processors (DSP), application-specific integrated circuits (ASIC), field programmable gate arrays (FPGA), or other programmable logic devices, transistor logic devices, hardware components, or any combination thereof. The general-purpose processor may be a microprocessor or any conventional processor.

[0129] The method steps in the embodiments of the present application can be implemented by hardware or by a processor executing software instructions. The software instructions can be composed of corresponding software modules, which can be stored in random access memory (RAM), flash memory, read-only memory (ROM), programmable read-only memory (PROM), erasable programmable read-only memory (EPROM), electrically erasable programmable read-only memory (EEPROM), registers, hard disks, mobile hard disks, CD-ROMs or any other form of storage medium known in the art. An exemplary storage medium is coupled to the processor so that the processor can read information from the storage medium and write information to the storage medium. Of course, the storage medium can also be a component of the processor. The processor and the storage medium can be located in an ASIC.

[0130] In the above embodiments, it can be implemented in whole or in part by software, hardware, firmware or any combination thereof. When implemented using software, it can be implemented in whole or in part in the form of a computer program product. The computer program product includes one or more computer instructions. When the computer program instructions are loaded and executed on a computer, the process or function described in the embodiment of the present application is generated in whole or in part. The computer can be a general-purpose computer, a special-purpose computer, a computer network, or other programmable device. The computer instructions can be stored in a computer-readable storage medium or transmitted via the computer-readable storage medium. The computer instructions can be transmitted from one website, computer, server or data center to another website, computer, server or data center via a wired (e.g., coaxial cable, optical fiber, digital subscriber line (DSL)) or wireless (e.g., infrared, wireless, microwave, etc.) method. The computer-readable storage medium can be any available medium that a computer can access or a data storage device such as a server or data center that includes one or more available media integrated. The available medium can be a magnetic medium (e.g., a floppy disk, a hard disk, a tape), an optical medium (e.g., a DVD), or a semiconductor medium (e.g., a solid state drive (SSD)).

[0131] It will be understood that the various numerical numbers involved in the embodiments of the present application are merely distinctions for the convenience of description and are not intended to limit the scope of the embodiments of the present application.

Claims

1. An ultrasonic tomography method, characterized in that: The method comprises: S1, in an immersion environment, using a multi-element ultrasonic annular array to perform full-matrix data acquisition on a target to be measured, and determining frequency domain observation data at the location of each array element, wherein the frequency domain observation data includes observation data of N frequency points; S2, selecting the observation data of the i-th frequency point from the observation data of the N frequency points; the observation data of the N frequency points are arranged in order from low frequency to high frequency, and the frequency of the i-th frequency point is less than the frequency of the i+1-th frequency point; i is a natural number, i+1≤N; S3, determining a source coding strategy, wherein the source coding strategy includes determining a time interval and a spatial distribution position of simultaneous transmissions of any two adjacent sound sources; encoding the observation data of the two adjacent sound sources at the i-th frequency point according to the source coding strategy to obtain first data; S4, performing coding prediction on the two adjacent sound sources in the inversion process at the i-th frequency point based on the Green's function and the source coding strategy to obtain second data; the Green's function is solved based on the initial sound speed parameter image m0; S5, determining an objective function based on the first data and the second data; the objective function is used to quantify the difference between the observed data of the i-th frequency point and the synthetic data predicted based on the initial sound speed parameter image m0; S6, using the full waveform inversion algorithm to iteratively update the initial sound velocity parameter image m0 to determine the i-th frequency point sound velocity parameter image m i ; S7, the sound speed parameter image m i As the initial sound speed parameter image m0 at the i+1th frequency point, let i=i+1, execute S2-S6 process, traverse the observation data of the N frequency points, and calculate the initial sound speed parameter image m0 according to the sound speed parameter image m0. n An ultrasonic tomographic image of the target to be measured is obtained by inversion.

2. The method according to claim 1, characterized in that The method uses a multi-element ultrasonic annular array to perform full matrix data acquisition on the target to be measured, and determines the frequency domain observation data of the position of each array element; including: One array element in the annular array transmits an ultrasonic signal to the target to be measured, and all array elements receive echo signals, and all array elements are traversed in sequence, and there is no overlap between the ultrasonic signals transmitted by each array element; Get the observed sound pressure signal D at each array element location obs (t,x r ,x s ), the amount of data is N t ×N r ×N s ; The observed sound pressure signal is subjected to discrete Fourier transform to obtain a frequency domain observation data set d obs (f,x r ,x s ); Among them, N t is the number of time samples, N r is the number of receiving array elements and N s is the number of transmitting array elements; x s Represents the position coordinates of the emission point and x r is the position coordinate of the receiving point; t represents the signal sequence corresponding to the time domain {t i , i=1,2,…N t }, the unit is seconds; f represents the signal sequence corresponding to the frequency domain {f i , i=1,2,…N t }, unit is Hertz.

3. The method according to claim 1, characterized in that Determining the source coding strategy includes: The time interval between the emission of any two sound sources is determined as: τ = (2n + 1) / (4f i ); where n is an integer; Determine that the spatial distribution position information of the two sound sources are adjacent and The encoding function is determined as: where ω i is the corresponding frequency f i The angular frequency below.

4. The method according to claim 1 or 3, characterized in that The step of encoding the observation data of the two adjacent sound sources at the i-th frequency point to obtain first data according to the source coding strategy includes: The source coding strategy is used to encode the observation data of the two adjacent sound sources of the i-th frequency point obtained in the full matrix acquisition mode, and the encoded data are superimposed to obtain the first data: where s j The sequence {s j , j=1,2…N s }, s j ′ represents the sequence of sound sources after coding and superposition {s j ′, j′=1,2…N′ s }, N s′ Indicates the total number of multi-source coding transmissions. When two sound sources are superimposed for coding 5. The method according to claim 1 or 3, characterized in that The step of performing coding prediction on the two adjacent sound sources in the inversion process at the i-th frequency point based on the Green's function and the source coding strategy to obtain the second data includes: For two adjacent sound source signals and The source coding strategy is used to encode and transmit the second data sequentially at the same frequency and at the time interval τ to obtain the second data as follows: Where G represents the Green's function, which is solved by the Helmholtz equation. The corresponding expression is as follows: Where ω represents the angular frequency, m represents the sound velocity parameter image of the medium, Δ represents the Laplace operator, and x is the coordinate of the calculation domain. is the Dirac delta function.

6. The method according to claim 1 or 3, characterized in that The determining of the objective function according to the first data and the second data includes: Determine the single frequency point f i The objective function C is described below se (m) is calculated as follows: Among them, the L2 norm is used to characterize the first data With the second data difference.

7. The method according to claim 1 or 3, characterized in that The iterative updating of the sound speed parameter image m in the objective function includes: The initial sound velocity parameter image m0 is iteratively updated using a local gradient optimization algorithm as follows: Where k represents the number of iterations; α k is the iteration step size; is the descending direction; Gradient optimized using previous iteration information Get the current descent direction The gradient is calculated using the adjoint state method as follows: where u se (x) and They represent the forward wavefield after source coding and the residual reverse wavefield at the receiving point in the inversion process respectively; R represents the real part; * represents the conjugate.

8. An electronic device comprising: at least one memory for storing a program; At least one processor is configured to execute a program stored in a memory. When the program stored in the memory is executed, the processor is configured to execute the method according to any one of claims 1 to 7.

9. A computer storage medium storing instructions, wherein when the instructions are executed on a computer, the computer is caused to execute the method according to any one of claims 1 to 7.

Images