A plane wave ultrasonic forward wave field construction method based on frequency-wave number domain extrapolation and reverse time migration imaging method

Abstract

The application discloses a plane wave ultrasonic forward wave field construction method and reverse time migration imaging method based on frequency-wavenumber domain extrapolation. A transmission boundary space-time signal is constructed according to the array element position, transmission pulse and deflection angle, a frequency-lateral wavenumber domain boundary spectrum is obtained through preprocessing and two-dimensional Fourier transform; starting from a two-dimensional scalar wave equation, an axial wavenumber is calculated according to the medium sound velocity and dispersion relation, a transmittable spectrum component is screened, and phase extrapolation is applied according to a target depth; through full frequency spectrum recovery, inverse transform, clipping and time alignment, a forward propagation wave field containing lateral position, depth and time dimensions is obtained. Zero time delay cross correlation is performed on the forward propagation wave field and a reverse propagation wave field obtained by time reversal of received data, and multi-angle results are coherently superimposed to obtain a coherent composite reverse time migration imaging result. The application can reduce the numerical calculation and storage cost of the forward wave field at each time step, and keep the propagation phase and imaging time synchronized.

Description

Technical Field

[0001] This invention belongs to the fields of ultrasound imaging, wavefield reconstruction and frequency domain signal processing, and particularly relates to a method for constructing a plane wave ultrasound positive wavefield outside the phase of the frequency-wavenumber domain, and a method for using the constructed wavefield for reverse time migration imaging. Background Technology

[0002] Plane wave ultrasound imaging generates approximate plane waves with preset deflection angles by applying corresponding transmission delays to each element of a linear array transducer. A single transmission can cover a large imaging area, featuring high frame rates and ease of multi-angle coherent recombination. The received array radio frequency data contains propagation time, phase, and amplitude information of the target scatterer, and focusing is typically achieved through delay summation or frequency domain shifting methods.

[0003] Reverse-time migration (RTM) is an imaging method based on a two-way wave equation. For any plane wave emission, traditional plane wave RTM typically requires calculating the forward propagation wavefield of the source in the medium, reversing the received radio frequency data in time to solve for the reverse propagation wavefield, and then using the zero-delay cross-correlation of the two wavefields to form an image. This method can utilize wavefield propagation information, but both forward and reverse propagation require time-step numerical advancement on a discrete spatial grid.

[0004] Under plane wave emission conditions, the emission boundary signal and its propagation direction are known in advance. Current implementations still commonly use finite difference time-domain methods, finite element methods, or k-space pseudospectral time-domain methods to perform complete numerical solutions for the forward wave field. As the imaging grid, number of time samples, and number of deflection angles increase, the computational and storage overhead of the forward wave field increases significantly, limiting the computational efficiency of multi-angle reverse time migration.

[0005] Frequency-wavenumber domain methods can decompose spatiotemporal wavefields into plane wave components with different frequencies and transverse wavenumbers, and perform phase propagation based on dispersion relations. Existing frequency-domain ultrasound methods are mostly used to directly perform migration or beamforming on received radio frequency data, focusing on the direct reconstruction from boundary echoes to the final image; their processing objects, output formats, and time variables differ from the complete forward propagation wavefield required for reverse time migration. If the frequency-domain migration results are directly used, it is impossible to perform zero-delay cross-correlation with the reverse propagation wavefield on a unified time axis.

[0006] Therefore, a frequency-wavenumber domain forward wavefield construction method for plane wave reverse time offset is needed to quickly extrapolate the transmitted boundary signal into a forward propagation wavefield with a complete time dimension within the imaging region, and to solve problems such as propagable spectrum screening, spatial periodic extension, time offset, and synchronization with the reverse propagation wavefield. Summary of the Invention

[0007] Objective: To address the problem of high computational and storage requirements in existing plane wave reverse-time migration methods where the forward wavefield relies on time-step numerical propagation, and to further provide a reverse-time migration imaging method using the aforementioned forward wavefield. This method converts the known transmitted boundary signal to the frequency-transverse wavenumber domain, performs depth phase extrapolation of each spectral component based on the dispersion relation of the wave equation, and then reconstructs the spatiotemporal wavefield corresponding to the grid size of the backpropagation wavefield.

[0008] The present invention also provides a plane wave ultrasound reverse time migration imaging method using the above-described plane wave ultrasound forward wave field construction method, in order to accelerate the computational efficiency of plane wave ultrasound reverse time migration imaging.

[0009] Technical solution: This invention provides a method for constructing a plane wave ultrasonic positive wavefield based on frequency-wavenumber domain extrapolation, comprising:

[0010] Step (1): Obtain the transmission boundary spatiotemporal signal of the linear array transducer at a preset plane wave deflection angle. The transmission boundary spatiotemporal signal includes the lateral position of the array element, the transmission pulse, and the transmission delay of the array element corresponding to the plane wave deflection angle.

[0011] Step (2): Perform at least one of the following preprocessing steps on the transmission boundary spatiotemporal signal: DC removal, bandpass filtering, time window truncation, lateral windowing, time dimension zero-padding, and lateral dimension zero-padding to obtain the preprocessed transmission boundary spatiotemporal signal;

[0012] Step (3): Perform Fourier transform on the preprocessed transmission boundary spatiotemporal signal along the time dimension and the array transverse dimension to obtain the frequency-transverse wavenumber domain boundary spectrum; save the conjugate symmetry relationship used to recover the full frequency spectrum before the inverse Fourier transform, while only retaining the positive frequency spectrum;

[0013] Step (4): The sound pressure in the two-dimensional homogeneous medium satisfies the two-dimensional scalar wave equation. Based on the medium sound velocity and sound wave dispersion relationship in the frequency domain-beam domain, the axial wave number is calculated for each frequency-transverse wave number combination in the frequency-transverse wave number domain boundary spectrum. A propagable spectral mask is constructed based on the axial wave number to determine the propagable spectral components and suppress the evanescent spectral components.

[0014] Step (5): According to the depth position of the wave field to be constructed, apply the axial propagation phase extrapolation factor determined by the axial wave number and propagation depth to the propagable spectral components to obtain the frequency-wave number domain wave field at the corresponding depth position.

[0015] Step (6): Recover the full frequency spectrum based on the conjugate symmetry relationship while retaining only the positive frequency spectrum. Perform inverse Fourier transform on the frequency-wavenumber domain wavefield at the corresponding depth position to recover the spatiotemporal domain wavefield. Perform at least one of the following processing on the recovered spatiotemporal domain wavefield: zero-filling region clipping, time offset correction, and time interpolation or resampling, so that the recovered spatiotemporal domain wavefield is consistent with the preset inverse time offset time grid.

[0016] Step (7): Repeat steps (5) and (6) for each depth position of the wave field to be constructed, and combine the spatiotemporal wave fields corresponding to each depth position into a plane wave forward propagation wave field within the imaging area that includes the lateral position, depth and time dimensions.

[0017] Further, in step (1), the spatiotemporal signal of the transmission boundary satisfies:

[0018] ,

[0019] in, In order to transmit boundary spatiotemporal signals, The horizontal position of the array element. To transmit pulses, The plane wave deflection angle. This refers to the aperture weight or element effectiveness weight. The array element transmission delay; the array element transmission delay satisfies:

[0020] ,

[0021] in, The horizontal position of the array element. The plane wave deflection angle. For background sound velocity, Fixed time offset for calculation error

[0022] Furthermore, in step (2), the horizontal windowing uses at least one of the Tukey window, Hann window, or cosine fading window; in step (3), the Fourier transform uses a two-dimensional fast Fourier transform to obtain the boundary spectrum. ,in For transverse wavenumber, ω is the angular frequency.

[0023] Furthermore, in step (4), within a uniform sound velocity medium, the sound pressure in a two-dimensional uniform medium... It satisfies the scalar two-dimensional scalar wave equation:

[0024] ,

[0025] right and After performing a Fourier transform, we obtain information about depth. Ordinary differential equations:

[0026] ,

[0027] The axial wavenumber is calculated according to the following formula:

[0028] ,

[0029] ,

[0030] Among them when When the corresponding frequency-transverse wavenumber combination is determined as the propagable spectral component; when At that time, the corresponding spectral component is set to zero, attenuated, or limited to a preset threshold. It is a propagable spectral mask.

[0031] Furthermore, in step (5), the depth position is... The frequency-wavenumber domain wave field at a given point satisfies:

[0032] ,

[0033] in, For a propagable spectral mask, The imaginary unit, The symbol parameters are determined by the Fourier transform symbol convention and the direction of wave field propagation.

[0034] Furthermore, step (6) includes: when only processing the positive frequency spectrum, recovering the full frequency spectrum using the conjugate symmetry relationship of the real signal spectrum; removing the zero-filled regions of the transverse and time dimensions from the recovered spatiotemporal wavefield; performing time alignment based on the transmission start time, fixed time offset, angular delay caused by the array aperture, and sampling start offset; when the target sampling time does not fall on the original sampling point, using linear interpolation, spline interpolation, or frequency domain resampling, resampling the spatiotemporal wavefield to the sampling time of the backpropagation wavefield.

[0035] Beneficial effects:

[0036] 1. This invention addresses the characteristic of known plane wave emission boundaries by constructing a forward propagation wave field through frequency-wavenumber domain decomposition and depth phase extrapolation, eliminating the need for time-step numerical advancement of the forward wave field and reducing the computation time of forward propagation in multi-angle reverse time migration.

[0037] 2. Starting from the two-dimensional scalar wave equation, the present invention determines the axial wave number by the dispersion relation in the frequency domain-beam domain. Each frequency-transverse wave number component propagates according to the corresponding phase, which can maintain the evolution relationship of the plane wave propagation direction, time phase and depth.

[0038] 3. The present invention outputs an internal positive wave field containing the complete time dimension, and synchronizes it with the back propagation wave field through time offset correction and resampling. It can be directly applied to zero-delay cross-correlation imaging, which is different from the direct imaging method that only outputs the final frequency domain offset image.

[0039] 4. This invention supports multi-angle plane waves, spatial and temporal zero-padding, propagation spectrum screening, and real signal conjugate symmetry recovery, and can be flexibly configured according to calculation accuracy, imaging range, and operating platform.

[0040] 5. This invention uses Fast Fourier Transform and element-wise complex multiplication as the main computational operations, which is suitable for parallel implementation on CPU, GPU, DSP and FPGA platforms, and can be combined with existing ultrasonic acquisition systems and inverse wave field solvers.

[0041] The plane wave ultrasonic reverse time migration imaging method provided by this invention can adopt the following technical solutions:

[0042] A plane wave ultrasound reverse time migration imaging method using the aforementioned plane wave ultrasound forward wave field construction method includes:

[0043] Step R1: Construct the forward propagation wave field corresponding to the plane wave deflection angle using the plane wave ultrasonic forward wave field construction method described above;

[0044] Step R2: Obtain the array received radio frequency data corresponding to the plane wave deflection angle, perform time reversal on the array received radio frequency data and extrapolate it in reverse based on the acoustic wave equation to obtain the reverse propagation wave field;

[0045] Step R3: Unify the forward propagation wave field and the reverse propagation wave field corresponding to the same plane wave deflection angle to the same spatial grid and time grid, and perform zero-delay cross-correlation to obtain a single-angle reverse time migration image;

[0046] Step R4 involves coherently stacking single-angle counter-time offset images with multiple deflection angles, and then performing envelope detection, normalization, and logarithmic compression on the stacking results to obtain the final display image.

[0047] Furthermore, in step R3, a single-angle reverse-time offset image is obtained:

[0048] ,

[0049] in Deflection angle The reverse time-shifted image, For a forward propagating wave field, For the backpropagating wave field, Total sampling time The time sampling interval is...

[0050] The final image is displayed in step R4:

[0051] ,

[0052] ,

[0053] in, For multi-angle coherent composite images, This represents the number of plane wave deflection angles. For the first A plane wave deflection angle, This is the final display image after logarithmic compression. This represents the Hilbert transform.

[0054] Beneficial effects:

[0055] 1. The constructed forward propagation wavefield and reverse propagation wavefield perform zero-delay cross-correlation on a unified spatial and temporal grid, which can maintain the matching of the propagation time and phase information of the scattered echo, so that the forward propagation wavefield constructed in the frequency-wavenumber domain can be directly used for reverse time migration imaging.

[0056] 2. After coherently superimposing single-angle reverse time-shifted images with multiple deflection angles, compared with the 0° single-angle imaging result, the main lobe of the point target is more concentrated, the artifacts in the non-target area are suppressed, and the boundary and internal contrast of the low-echo target are clearer, thereby improving the resolution and contrast performance of the multi-angle reverse time-shifted images.

[0057] 3. Under the same conditions of receiving radio frequency data, backpropagation process, and imaging, the forward wavefield construction time for point target models with 21 deflection angles was reduced from 460.0 s to 76.1 s, a speedup of 6.05 times and a time reduction of 83.46%; the time for point target and low-echo target models was reduced from 728.8 s to 77.1 s, a speedup of 9.45 times and a time reduction of 89.42%. These efficiency improvements are limited to the forward wavefield construction stage, which can reduce the computational burden of multi-angle reverse time migration imaging.

[0058] The present invention also provides a computer device, including a memory, a processor, and a computer program stored in the memory and executable on the processor. When the processor executes the computer program, it implements the aforementioned method for constructing a plane wave ultrasonic forward wavefield based on frequency-wavenumber domain extrapolation and the method for reverse time migration imaging.

[0059] The present invention also provides a computer-readable storage medium having a computer program stored thereon, wherein when the computer program is executed by a processor, the computer program implements the aforementioned method for constructing a plane wave ultrasonic forward wavefield based on frequency-wavenumber domain extrapolation and the method for reverse time migration imaging. Attached Figure Description

[0060] Figure 1 This is a flowchart of the overall process of constructing the forward wavefield of a plane wave based on frequency-wavenumber domain extrapolation and reverse time migration imaging in this invention.

[0061] Figure 2 A schematic diagram for constructing the boundary signal of a plane wave transmission;

[0062] Figure 3 A schematic diagram illustrating the transformation of the transmitted boundary signal from the time-space domain to the frequency-transverse wavenumber domain.

[0063] Figure 4 This is a schematic diagram of depth phase extrapolation and spatiotemporal wavefield reconstruction based on axial wavenumber.

[0064] Figure 5 This is a schematic diagram of the forward wavefield reverse time migration imaging results constructed in this invention;

[0065] Figure 6 The graph shows a comparison of computation time and speedup ratio between numerical forward propagation using existing technology and forward wavefield construction using the frequency-wavenumber domain of this invention. Detailed Implementation

[0066] The present invention will be further described below with reference to the accompanying drawings and specific embodiments. The embodiments are used to illustrate the technical solutions of the present invention, and the listed parameters are one feasible configuration. The transducer structure, transmission bandwidth, number of deflection angles, imaging area, zero-padding factor, and computing platform can be adjusted according to the actual system.

[0067] Example 1: Construction of the forward propagation wavefield of a plane wave based on frequency-wavenumber domain extrapolation

[0068] like Figure 1 As shown, this embodiment provides a method for constructing the forward ultrasonic wavefield of a plane wave based on frequency-wavenumber domain extrapolation, including the following steps:

[0069] Step (1) Set the number of array elements, the lateral coordinates of the array elements, the spacing between array elements, the effective emission aperture, the emission pulse, the sampling frequency, the time sampling interval, and the background sound velocity of the linear array transducer. The imaging lateral range, imaging depth grid, and target time axis are specified. This embodiment uses an array element count of N=128 and an element spacing of... The center frequency is Background sound speed and The sampling frequency, the horizontal coordinates of the array elements, the transmission delay of the array elements, and the spatiotemporal signal of the transmission boundary are determined according to the following formulas:

[0070] ,

[0071] ,

[0072] ,

[0073] in, Number the array elements. For the first Aperture weights of individual array elements; Calculate the fixed time offset caused by the calculation error. Arrange the transmitted pulse according to... After translating along the time axis, a spatiotemporal matrix of the emission boundary can be formed, which represents the time sampling point × the position of the array element.

[0074] Figure 2 The results of transmitting signal construction are given when the deflection angles are -10°, 0° and 10°. Figure 2 The left column represents the delay pulses on each transmission channel. Figure 2 The right column represents the spatiotemporal signals of the transmission boundary formed after the array elements are arranged in the lateral position. When the deflection angle is 0°, each array element transmits synchronously, and the boundary signal is approximately horizontally distributed in the time-space plane; when the deflection angle is positive or negative, the transmission time of the array elements changes monotonically with the lateral position, and the boundary signal is inclined in opposite directions, thereby determining the propagation direction of the plane wave.

[0075] Step (2) involves truncating the boundary signal using a time window based on the transmit pulse bandwidth and effective duration, and performing at least one of the following processes: DC removal, effective bandwidth limitation, and amplitude normalization. A Tukey window is applied at the edge of the array's effective aperture to reduce lateral spectral leakage caused by aperture truncation; subsequently, in the lateral dimension... and time dimension Zero-padding is performed before and after the sample. The number of samples after zero-padding satisfies:

[0076] ,

[0077] ,

[0078] in, Take 2 to 8, Choose 1 to 2. In this embodiment, the corresponding zero-padding multiple is selected according to the target lateral range and the maximum propagation time, and the original aperture range, effective time start point and index mapping relationship are retained for the trimming and time alignment in step (6).

[0079] Step (3) Perform a two-dimensional fast Fourier transform on the preprocessed emission boundary spatiotemporal matrix along the time dimension and the array transverse dimension to obtain the frequency-transverse wavenumber domain boundary spectrum:

[0080] ,

[0081] according to The discrete angular frequency vector is established using the sampling frequency and the zero-padding time. ,according to The discrete transverse wavenumber vector is established by the element spacing and the zero-padded transverse length. For real-value transmitted signals, this embodiment preserves the complete positive and negative frequency spectra; in the implementation that only processes the positive frequency spectrum, before the inverse transformation in step (6), the following is used: Restore the negative frequency spectrum.

[0082] like Figure 3 As shown, the left column displays the frequency-spatial domain spectrum at different deflection angles, with the main energy concentrated in the effective frequency band of the transmitted pulse; the right column displays the frequency-transverse wavenumber domain spectrum obtained after further transverse position transformation. At a deflection angle of 0°, the main spectral energy is concentrated near the zero transverse wavenumber; at deflection angles of -10° and 10°, the main spectral energy shifts in opposite directions along the transverse wavenumber axis, and its distribution corresponds to the transverse propagation component of the plane wave. The positive and negative frequency spectra maintain a conjugate symmetry relationship about the origin of the spectrum.

[0083] Step (4) In the frequency domain-beam domain, for each angular frequency With transverse wavenumber Combined, the axial wavenumber is calculated based on the dispersion relation of acoustic waves in a homogeneous medium, and a propagable spectral mask is constructed:

[0084] ,

[0085] ,

[0086] when hour, are real numbers and ;when At that time, the combination is identified as a fleeting spectral component and then... At the same time, a frequency mask is set outside the effective frequency band of the transmitted pulse to prevent invalid frequency components from participating in phase extrapolation.

[0087] Step (5) For each target depth within the imaging region A complex exponential phase extrapolation factor is constructed based on the axial wavenumber and propagation depth. The Fourier transform convention used in this embodiment corresponds to the phase factor when propagating in the positive depth direction as follows: Therefore, depth The frequency-wavenumber domain wave field at that point is:

[0088] ,

[0089] The phase extrapolation factor is multiplied element-by-element by the emission boundary spectrum and the propagable spectrum mask, so that the phase of different frequency-transverse wavenumber components is accumulated and propagated according to their respective axial wavenumbers.

[0090] Step (6) When processing only the positive frequency spectrum, first recover the full frequency spectrum based on the conjugate symmetry relationship, and then perform an inverse two-dimensional Fourier transform on the frequency-wavenumber domain wavefield at the target depth:

[0091] ,

[0092] After the inverse transformation, the real part is taken, and the zero-filled regions in the horizontal and time dimensions are removed, so that the forward propagation wave field and the reverse propagation wave field correspond point by point in spatial coordinates, time intervals and starting time.

[0093] Figure 4 The left column gives the depth The frequency-wavenumber domain wavefield after phase extrapolation is shown. The right column presents snapshots of the spatial wavefield selected according to the theoretical arrival times for each deflection angle. The example times corresponding to -10° and 10° are... The example time corresponding to 0° is The recovered wavefront maintains an approximately planar shape within the imaging region. The 0° wavefront is parallel to the array, while the positive and negative deflection wavefronts are tilted in opposite directions, indicating that frequency-wavenumber domain phase extrapolation can recover the internal positive propagation wavefield consistent with the set emission angle.

[0094] Step (7) Repeat steps (5) and (6) according to the imaging depth grid, and then... The wavefields at each depth location are recovered and time-aligned, and then written into a 3D wavefield array:

[0095] ,

[0096] This results in a size of Complete forward propagation wave field The calculation process can be performed layer by layer according to depth, or the phase extrapolation factors of multiple depths can be combined into a three-dimensional array for batch complex multiplication and inverse transformation. When storage resources are limited, the wavefield can be recovered in blocks according to depth or time, and after recovery, it can be directly cross-correlated and accumulated with the corresponding backpropagation wavefield.

[0097] Example 2: Plane Wave Reverse Time Migration Imaging and Simulation Verification

[0098] This embodiment provides a plane wave ultrasound reverse time migration imaging method based on the plane wave ultrasound forward wave field construction method described in Embodiment 1.

[0099] This plane-wave ultrasonic reverse time migration imaging method uses Field II to establish a linear array transducer model and generate array-received radio frequency data. (Number of elements in the linear array transducer is missing.) The center frequency is The spacing between array elements is The gap between array elements is The array element height is The sampling frequency is The background sound velocity is The transducer pulse response uses a Gaussian pulse with a relative bandwidth of 0.67, and the transmission excitation duration is one center frequency cycle. Both transmission and reception use full aperture, and different deflection angles are controlled only by the element transmission delay.

[0100] The multi-angle simulation uses 21 plane wave deflection angles within the range of -10° to 10°, with adjacent angles spaced 1° apart; Figures 2 to 4 -10°, 0°, and 10° were selected as typical angles to illustrate the boundary signal, spectral distribution, and depth extrapolation results. The k-Wave spatial grid was set with approximately 3 grid points per wavelength, and a fully matched layer with a thickness of 15 grid points was set around the perimeter of the computational region to absorb sound waves propagating to the boundary and reduce artificial boundary reflections.

[0101] The point target model contains 3 scattering points, all of which are located laterally. The depths are respectively , and The scattering amplitude is set to 1 for all values. Point target and low-echo target models are set up in a random scattering background to simultaneously evaluate target localization, spatial resolution, sidelobe suppression, and low-echo region contrast. Field II output data is arranged by time sampling points and receiver elements, and DC removal, effective bandwidth limitation, and start-time correction are performed on the RF data at each angle.

[0102] The model for point targets and low-echo targets includes three strong echo scattering points and two low-echo targets. The lateral positions of the strong echo scattering points are as follows: , and The depth position is The lateral center position of the low-echo target is The depth positions are respectively and Their radii are respectively and .

[0103] In traditional comparison methods, both the forward and backward propagation wavefields are numerically calculated by k-Wave based on the same sound velocity model and spatial grid. In the method of this invention, the backward propagation wavefield is still calculated by k-Wave, while the forward propagation wavefield is constructed using the frequency-wavenumber domain extrapolation method described in Example 1, and mapped to an imaging grid consistent with k-Wave through two-dimensional interpolation. Both methods use the same imaging area, received radio frequency data, backward propagation process, and imaging conditions, ensuring that the comparison results primarily reflect the differences in the construction method of the forward wavefield.

[0104] The forward propagation wave field corresponding to the same deflection angle θ and the backpropagating wave field By multiplying the results at each imaging grid point at different time intervals and accumulating them along the time dimension, a single-angle zero-delay cross-correlation image is obtained:

[0105] ,

[0106] After constructing the forward wavefield, performing backpropagation, and cross-correlation imaging at each of the 21 deflection angles, the images from each single angle are coherently superimposed.

[0107] ,

[0108] After superimposing the results, the envelope is obtained, followed by amplitude normalization and logarithmic compression:

[0109] ,

[0110] in, This represents the number of plane wave deflection angles. The time sampling interval is... This represents the Hilbert transform. This yields the final FK-RTM display image.

[0111] Figure 5 The top row shows the 0° single-angle FK-RTM image and the 21° coherent superimposed image of the point target model, while the bottom row shows the corresponding results of the point target and the low-echo target model. After multi-angle coherent superposition, the main lobe of the point target is more concentrated, artifacts in the non-target area are suppressed, and the boundary and internal contrast of the low-echo area are clearer, indicating that the constructed forward propagation wavefield can directly participate in reverse time migration cross-correlation imaging.

[0112] Figure 6 A comparison of the forward wavefield construction time for two simulation models is presented. The timing range only includes the forward wavefield construction process; the time for Field II RF data generation, backpropagation, cross-correlation imaging, and post-processing is not included. The traditional numerical forward modeling time is the k-Wave forward propagation calculation time; the time of the method in this invention includes the source matrix construction, frequency-wavenumber domain core extrapolation, and interpolation mapping to a unified k-Wave grid.

[0113] For the point target model, the traditional numerical forward modeling at 21 angles takes a cumulative time of 460.0 s, while the method of this invention takes a cumulative time of 76.1 s, resulting in an overall speedup of 6.05 times and a reduction in forward wavefield construction time of approximately 83.46%. For the point target and low-echo target models, the traditional numerical forward modeling takes a cumulative time of 728.8 s, while the method of this invention takes a cumulative time of 77.1 s, resulting in an overall speedup of 9.45 times and a reduction in forward wavefield construction time of approximately 89.42%. Under different deflection angles, the time consumption of the method of this invention remains at a relatively stable low level, verifying the computational efficiency advantage of this method in multi-angle plane wave reverse-time migration.

[0114] In summary, the method of the present invention, while maintaining the synchronization relationship of plane wave propagation direction, frequency components, time phase and reverse time offset, can effectively replace the time-step numerical propagation of the forward wave field by using fast Fourier transform, element-wise complex multiplication and inverse transform as the main operations, and can be implemented in combination with existing Field II, k-Wave or other ultrasonic simulation and imaging platforms.

Claims

1. A method for constructing a plane wave ultrasonic positive wavefield based on frequency-wavenumber domain extrapolation, characterized in that, include: Step (1): Obtain the transmission boundary spatiotemporal signal of the linear array transducer at a preset plane wave deflection angle. The transmission boundary spatiotemporal signal includes the lateral position of the array element, the transmission pulse, and the array element transmission delay corresponding to the plane wave deflection angle. Step (2) Perform at least one of the following preprocessing steps on the transmission boundary spatiotemporal signal: DC removal, bandpass filtering, time window truncation, lateral windowing, time dimension zero-padding, and lateral dimension zero-padding, to obtain the preprocessed transmission boundary spatiotemporal signal. Step (3): Perform Fourier transform on the preprocessed transmission boundary spatiotemporal signal along the time dimension and the array transverse dimension to obtain the frequency-transverse wavenumber domain boundary spectrum; save the conjugate symmetry relationship used to recover the full frequency spectrum before the inverse Fourier transform, while only retaining the positive frequency spectrum; Step (4): In a two-dimensional homogeneous medium, the sound pressure satisfies the two-dimensional scalar wave equation. Based on the medium sound velocity and sound wave dispersion relationship in the frequency domain-beam domain, the axial wave number is calculated for each frequency-transverse wave number combination in the frequency-transverse wave number domain boundary spectrum. A propagable spectral mask is constructed based on the axial wave number to determine the propagable spectral components and suppress the evanescent spectral components. Step (5): According to the depth position of the wave field to be constructed, apply the axial propagation phase extrapolation factor determined by the axial wave number and propagation depth to the propagable spectral components to obtain the frequency-wave number domain wave field at the corresponding depth position. Step (6): Recover the full frequency spectrum based on the conjugate symmetry relationship while retaining only the positive frequency spectrum. Perform inverse Fourier transform on the frequency-wavenumber domain wavefield at the corresponding depth position to recover the spatiotemporal domain wavefield. Perform at least one of the following processing on the recovered spatiotemporal domain wavefield: zero-filling region clipping, time offset correction, and time interpolation or resampling, so that the recovered spatiotemporal domain wavefield is consistent with the preset inverse time offset time grid. Step (7): Repeat steps (5) and (6) for each depth position of the wave field to be constructed, and combine the spatiotemporal wave fields corresponding to each depth position into a plane wave forward propagation wave field within the imaging area that includes the lateral position, depth and time dimensions.

2. The method for constructing a plane wave ultrasonic positive wave field according to claim 1, characterized in that, In step (1), the spatiotemporal signal at the transmission boundary satisfies: ， in, In order to transmit boundary spatiotemporal signals, The horizontal position of the array element. To transmit pulses, The plane wave deflection angle. This refers to the aperture weight or element effectiveness weight. The array element transmission delay; the array element transmission delay satisfies: ， in, The horizontal position of the array element. The plane wave deflection angle. For background sound velocity, A fixed time offset is used to calculate the error.

3. The method for constructing a plane wave ultrasonic positive wave field according to claim 1, characterized in that, In step (2), the horizontal windowing uses at least one of the Tukey window, Hann window, or cosine fading window; in step (3), the Fourier transform uses a two-dimensional fast Fourier transform to obtain the boundary spectrum. ,in For transverse wavenumber, ω is the angular frequency.

4. The method for constructing a plane wave ultrasonic positive wave field according to claim 1, characterized in that, In step (4), within a uniform sound velocity medium, the sound pressure in a two-dimensional uniform medium... It satisfies the scalar two-dimensional scalar wave equation: ， right and After performing a Fourier transform, we obtain information about depth. Ordinary differential equations: ， The axial wavenumber is calculated according to the following formula: ， ， Among them when When the corresponding frequency-transverse wavenumber combination is determined as the propagable spectral component; when At that time, the corresponding spectral component is set to zero, attenuated, or limited to a preset threshold. It is a propagable spectral mask.

5. The method for constructing a plane wave ultrasonic positive wave field according to claim 1, characterized in that, In step (5), the depth position is The frequency-wavenumber domain wave field at a given point satisfies: ， in, For a propagable spectral mask, The imaginary unit, The symbol parameters are determined by the Fourier transform symbol convention and the direction of wave field propagation.

6. The method for constructing a plane wave ultrasonic positive wave field according to claim 1, characterized in that, Step (6) includes: when only processing the positive frequency spectrum, recovering the full frequency spectrum using the conjugate symmetry relationship of the real signal spectrum; removing the zero-filled regions of the transverse and time dimensions from the recovered spatiotemporal wavefield; performing time alignment based on the transmission start time, fixed time offset, angular delay caused by the array aperture, and sampling start offset; when the target sampling time does not fall on the original sampling point, using linear interpolation, spline interpolation, or frequency domain resampling, resampling the spatiotemporal wavefield to the sampling time of the backpropagation wavefield.

7. A plane wave ultrasound reverse time migration imaging method using the plane wave ultrasound forward wavefield construction method according to any one of claims 1 to 6, characterized in that, include: Step R1: Construct the forward propagation wave field corresponding to the plane wave deflection angle using the plane wave ultrasonic forward wave field construction method described above; Step R2: Obtain the array received radio frequency data corresponding to the plane wave deflection angle, perform time reversal on the array received radio frequency data and extrapolate it in reverse based on the acoustic wave equation to obtain the reverse propagation wave field; Step R3: Unify the forward propagation wave field and the reverse propagation wave field corresponding to the same plane wave deflection angle to the same spatial grid and time grid, and perform zero-delay cross-correlation to obtain a single-angle reverse time migration image; Step R4: Coherently overlay single-angle counter-time offset images with multiple deflection angles, and perform envelope detection, normalization, and logarithmic compression on the overlay results to obtain the final display image.

8. The plane wave ultrasonic reverse time migration imaging method according to claim 7, characterized in that, In step R3, a single-angle reverse-time offset image is obtained: ， in Deflection angle The reverse time-shifted image, For a forward propagating wave field, For the backpropagating wave field, Total sampling time The time sampling interval is... The final image is displayed in step R4: ， ， in, For multi-angle coherent composite images, This represents the number of plane wave deflection angles. For the first A plane wave deflection angle, This is the final display image after logarithmic compression. This represents the Hilbert transform.

9. A computer device comprising a memory, a processor, and a computer program stored in the memory and executable on the processor, characterized in that, When the processor executes the computer program, it implements the method according to any one of claims 1 to 8.

10. A computer-readable storage medium having a computer program stored thereon, characterized in that, When the computer program is executed by a processor, it implements the method described in any one of claims 1 to 8.

Images