A High Signal-to-Noise Ratio Imaging Method for Synthetic Aperture Data

Abstract

The present invention provides a high signal-to-noise ratio imaging method for synthetic aperture data. First, the STA data of the input biological material is preprocessed to filter out the direct wave; then the preprocessed STA data is extrapolated to the target layer to obtain the STA data of the target layer; the singular value decomposition (SVD) is used for the STA data of the target layer to extract the signal subspace, and then the cross-spectrum matrix (CMS) corresponding to each frequency is calculated; finally, the cross-spectrum matrix of each frequency is deconvolved by using the unsteady point spread function (PSF) to obtain the imaging result corresponding to each frequency, and then the imaging result of the biological material is obtained by iterative calculation frequency by frequency. The imaging method of the present invention uses SVD to eliminate most of the noise in the STA data, and uses PSF for deconvolution to effectively suppress the speckle artifacts in the imaging process, greatly improving the signal-to-noise ratio, enhancing the imaging quality, with small computational amount and high imaging efficiency.

Description

Technical Field

[0001] The present invention belongs to the technical field of detection methods, and particularly relates to a high signal-to-noise ratio imaging method for synthetic aperture data. Background Art

[0002] Ultrasound is an important imaging tool in clinical diagnosis because it can provide real-time anatomical and functional images without causing any side effects. With the development of electronic information technology, synthetic transmit aperture (STA) data has become a research hotspot. It contains the responses of each pair of transmit-receive pairs in a phased array. This data set has been widely used in various scenarios, such as tissue structure imaging, flow estimation, and tissue motion compensation. STA data set is usually obtained by sequentially exciting a single element and receiving by all elements. It uses two-way dynamic focusing for imaging, which greatly improves the lateral spatial resolution and contrast because it adopts the real wave propagation path. In addition, STA data set can be recovered from other beamforming strategies. Generally, they use multiple elements to be excited simultaneously according to certain weights, so that the signal-to-noise ratio (SNR) and frame rate can be improved.

[0003] The STA data contains the maximum possible independent information about the measurement area from the phased array probe. Appropriate imaging algorithms are needed to process it to reconstruct the position of the reflector. The delay and sum (DAS) method is a standard method in medical imaging. It estimates the propagation time by treating the beam as a ray, and superimposes the data according to the time delay and the weights obtained by minimum variance optimization to reconstruct the target. This method has certain advantages in terms of robustness and signal-to-noise ratio, but when the measurement object is a multi-layer structure, time-consuming path calculation is required.

[0004] In addition, the rapid development of machine learning has also provided new options for ultrasonic imaging. However, most researchers use machine learning as a means of preprocessing or postprocessing the DAS results, and the neural networks (NNs) used for ultrasonic data set image reconstruction are very limited. Here, Hyun uses the NN framework to process STA data to obtain imaging results, which can suppress speckle artifacts; Youn also preliminarily attempts to directly process ultrasonic signals through a convolutional neural network to detect scatterers. Since the STA data set contains a large amount of signals, a huge network is needed to obtain a correct and somewhat general mapping from input to output, which depends on a large and diverse training data set. Therefore, the huge workload currently hinders the popularization of this technology.

[0005] Full matrix capture (FMC) evolved from the STA method and has been widely used in industrial inspections. Since the FMC and STA datasets are the same, the new imaging methods of FMC can be transformed into STA, including the phase shift migration (PSM) method, reverse time migration (RTM) method, and model-based method. The PSM method regards the received signal as the wave field on the surface of the measurement area, separates each harmonic component using Fourier transform, reconstructs the wave field of the measurement area through the phase shift operator, and then implements the imaging condition to obtain the focused image of the target. When measuring the defects of multi-layer structures, the PSM shows great efficiency advantages. The RTM method uses the full wave equation to reconstruct the wave field of the measurement area and then executes the imaging condition to generate a focused image. The frequency domain RTM can solve different excitations simultaneously, which significantly reduces the computational complexity. This method is not limited by the complexity of the sound speed model and can obtain high-resolution focused images, but it is very sensitive to the accuracy of the sound speed model. The model-based method uses matrix equations to establish the relationship between the FMC dataset and the measurement target, and adopts an inversion algorithm to recover a high-resolution image. Nevertheless, due to the scale of the FMC or STA dataset, the matrix equation is extremely large, making the computational cost of this method very high and difficult to apply.

[0006] The above imaging methods all assume that the sound speed in the measurement area is uniform. However, the sound speed in biological tissues is usually slightly non-uniform, and the random interference of sub-resolution scattering will cause speckle artifacts in the imaging results. Speckles are usually regarded as noise that deteriorates the imaging results because it reduces the resolution and signal-to-noise ratio of the target. Therefore, in order to obtain a target image with a high signal-to-noise ratio, it is necessary to largely suppress the speckle noise. However, there is no imaging method in the prior art that can well suppress the speckle noise, and the signal-to-noise ratio of the obtained target image is relatively low. Summary of the Invention

[0007] To solve the problems existing in the prior art, the present invention provides a high signal-to-noise ratio imaging method for synthetic aperture data, which can not only eliminate most of the noise in the STA data, but also suppress the speckle artifacts in the imaging process, effectively improving the signal-to-noise ratio and resolution of the imaging results.

[0008] A high signal-to-noise ratio imaging method for synthetic aperture data (STA) includes the following steps:

[0009] (1) Input the synthetic aperture data, shift the synthetic aperture data to the upper surface of the target layer to obtain the target layer data;

[0010] (2) Perform singular value decomposition (SVD) on the target layer data, extract the signal subspace, and calculate the cross-spectrum matrix (CMS) for each frequency;

[0011] (3) Deconvolve the cross-spectrum matrix for each frequency to obtain the imaging result corresponding to each frequency, and superimpose the imaging results corresponding to each frequency to obtain the final imaging result.

[0012] In the above steps, the signal subspace and the noise subspace are divided by analyzing the eigenvalues, and any division method can be used. For a homogeneous medium, in actual operation, it is not necessary to offset the synthetic aperture data, and the subsequent processing can be directly performed on the measured synthetic aperture data.

[0013] Preferably, in step (2), the form of performing singular value decomposition on the target layer data is:

[0014]

[0015] where D(ω) is the target layer data at frequency ω; the subscripts sig and n represent the signal subspace and the noise subspace respectively; ω is the frequency; U sig (ω), V sig (ω) represent two groups of eigenvectors of the signal subspace respectively; ∑ sig (ω) is the eigenvalue of the signal subspace; U n (ω), V n (ω) are two groups of eigenvectors of the noise subspace respectively; ∑ n (ω) is the eigenvalue of the noise subspace; the superscript H represents transpose and complex conjugate.

[0016] The synthetic transmit aperture (STA) uses each element in the phased array to sequentially transmit spherical waves to the target area, and at the same time, all elements in the phased array receive the reflected waves from the medium. For a linear array of N elements, the acquired signal obtained in this acquisition method contains three dimensions, which are represented here as where N t is the number of sampling points in the time domain. Therefore, the STA data is a response matrix that contains the response relationships between all elements. As shown in (a) of Figure 2 , the wave is emitted from the element located at R s =(x s ,0), propagates in the medium, reflects after encountering a reflector, and the reflected wave is received by the element located at R r =(x r ,0).

[0017] The propagation of sound waves between the transmitting element R s and the receiving element R r can be modeled in the frequency domain as:

[0018]

[0019] where M is the number of reflectors, and the reflector m is located at R m =(x m , z m ), and the reflection coefficient is ε m ; F(ω) is the spectrum of the excitation signal; D s (θ s , ω) and D r (θ r , ω) are the directivity coefficients of the transmitting element and the receiving element respectively; θ s and θ r are the azimuth angles of the transmitting element and the receiving element respectively; ω is the frequency.

[0020] Among them, the directivity coefficient is defined as:

[0021]

[0022] where a is the element width and c is the speed of sound.

[0023] G is the Green's equation, obtained from the following formula:

[0024]

[0025] where H0 is the Hankel equation of the first kind of zero order; i is the imaginary unit; the Green's equation G(R m , R r , ω) is the response to the point source at R m at the point R r .

[0026] In the frequency domain, the N×N response matrix of all elements can be written as:

[0027]

[0028] where:

[0029] g m =[D1(θ1,ω)g(R1,R m , ω), D2(θ2,ω)G(R2,R m , ω), …, D N (θ1,ω)G(R N , R m , ω)] T (5)

[0031] where H is the transpose and complex conjugate; T is the transpose.

[0032] G is an N×M matrix, expressed as:

[0033] G = [g1, g2, …, g M(6)

[0034] M is the number of reflectors; N is the number of vibration elements.

[0035] Since the directivity coefficients in different azimuths corresponding to the imaging area do not differ much, the directivity coefficient of the vibration element can be ignored, so G is an orthogonal matrix.

[0036] Q is an M×M matrix, written as:

[0037]

[0038] Therefore, the measured response matrix D(ω) can be modeled as:

[0039] D(ω) = [D(x s ,x r ,ω)] = GQG H +N(ω) (8)

[0040] where D(ω) is an N×N matrix, which is the set of all receive-transmit vibration element pairs in D(x s ,x r ,ω) at frequency ω, and D(x s ,x r ,ω) can be obtained by performing a one-dimensional Fourier transform on D(x s ,x r ,t) in the time dimension.

[0041] Here, it is assumed that the number of reflection points is less than the number of vibration elements, i.e., M < N. Therefore, P(ω) is an ideal response matrix of rank M, but this ideal response matrix is contaminated by noise and becomes D(ω), which is full rank. Theoretically, the first M eigenvalues of D(ω) are much larger than the remaining N - M eigenvalues.

[0042] The measured response matrix can be written in the following eigenvalue decomposition form:

[0043]

[0044] where the subscripts sig and n represent the signal subspace and the noise subspace respectively, and the signal subspace and the noise subspace are distinguished according to the magnitudes of the eigenvalues; ∑ sig (ω), ∑ n (ω) are the eigenvalues of the signal subspace and the noise subspace respectively, both of which are M×M diagonal matrices corresponding to M reflection points; U sig (ω) and V sig (ω) are two sets of eigenvectors of the signal subspace, which are N×M matrices. Theoretically, they should be the same, but the noise in actual measurement will cause differences between them. Here, V is ignored.sig (ω), take U sig (ω) as the signal subspace for subsequent calculation of the cross-spectral matrix; U n (ω), V n (ω) are respectively two eigenvectors of the noise subspace, and are also N×M dimensional matrices.

[0045] Therefore, SVD can separate the scattered signal from the noise, and the following formula can be obtained:

[0046]

[0047] Comparing the three terms in formula (10), V sig (ω), ∑ sig (ω) and U sig (ω) can be respectively interpreted as the signal of the incident wave emitted from the excitation element, the reflection coefficient of the reflector, and the received signal of the reflected wave.

[0048] Therefore, U sig (ω) = [u sig1 , …, u sigM can be regarded as the signal of the wave actively emitted from the scattering point received by the sensor, as shown in Figure 2 (b). Among them, u sigm is the signal vector in the signal subspace U sig (ω), m ∈ [1, M], and M is the number of reflectors.

[0049] Therefore, extracting U sig (ω) through SVD not only filters out most of the noise, but also transforms the passive heterogeneous detection problem into an active sound source localization problem, and then the cross-spectral matrix (CMS) can be used for imaging.

[0050] Preferably, in step (2), the calculation formula of the cross-spectral matrix C M (ω) is:

[0051]

[0052] In the formula, M is the number of reflection points; ω is the frequency; u sigm is the signal vector in the signal subspace; the superscript H is the transpose and complex conjugate.

[0053] Beamforming is an important method for active sound source localization. The wave equation of sound waves in a two-dimensional homogeneous medium can be written as:

[0054]

[0055] Among them, p(x, z, t) is the sound pressure field at point (x, z), q(t) is the point source located at (x0, z0), and the wave field in the time domain is as follows:

[0056]

[0057] Among them, the distance is The time delay is

[0058] Therefore, the signal received by the phased array element is the wave field p(x r , 0, t) on the surface, and the sound source distribution can be reconstructed by inverting the time delay.

[0059] The possible sound sources at point (x, z) can be estimated by the following formula:

[0060]

[0061] Among them, (x j , 0) is the position of the nth element; t j = r j / c.

[0062] For After performing the Fourier transform, the expression of beamforming in the frequency domain is as follows:

[0063]

[0064] Among them, h(x, z, ω) and Y(x, z, ω) are N×1 vectors, written as:

[0065] h(x, z, ω) = [2πr1exp(iωt1), …, 2πr N exp(iωt N )] T (15)

[0066] Y(x, z, ω) = [P(x1, 0, ω), …, P(x N , 0, ω)] T (16)

[0067] is a complex number, and its energy can be calculated by the following formula:

[0068]

[0069] Among them, the superscript * represents the conjugate; C(ω) is an N×N matrix containing the measured signals, C(ω) is the cross-spectral matrix (CMS), and any element C ij can also be obtained by the following formula:

[0070] C ij = P(x i , 0, ω)P * (x j , 0, ω) (18)

[0071] For STA data, U sig (ω) is obtained by the SVD method and is regarded as the signal received by the phased array after the wave actively emitted by the heterogeneous body (reflector) at time t = 0. The imaging condition is defined as follows:

[0072] I(x, z) = ∫I(x, z, ω)dω = ∫h H (x, z, ω)C M (ω)h(x, z, ω)dω (19)

[0073] Where:

[0074]

[0075] Where, u sigm is the signal vector in the signal subspace U sig (ω).

[0076] According to the extracted U sig (ω), the cross-spectral matrix C M (ω) is calculated first, and then the imaging is performed using the imaging condition and C M (ω) to obtain the imaging result corresponding to the frequency ω.

[0077] Preferably, in step (3), the cross-spectral matrix of each frequency obtained is deconvolved using the non-stationary point spread function (PSF).

[0078] The imaging of a point source by the phased array is called the point spread function (PSF), which is related to the hardware parameters, such as element size, element spacing, and element reception performance. The imaging result of the phased array is the convolution result of the PSF and the real image.

[0079] For a point source λ = (x λ , z λ ), a CSM introduced by the sound source is:

[0080] C λ (ω) = Y(x λ , z λ , ω)Y H (x λ , z λ , ω) (21)

[0081] Therefore, the point spread function PSF can be written as:

[0082] I λ (x, z, ω) = h H (x, z, ω)C λ (ω)h(x, z, ω) (22)

[0083] Among them, the PSF is non-steady because it is a function of the light source position. The PSF blurs the image, making adjacent scattering points indistinguishable. Therefore, it is necessary to use deconvolution to eliminate the influence of the PSF and improve the resolution.

[0084] Preferably, in step (3), when performing deconvolution on the cross-spectrum matrix of each frequency, the iteration is terminated when the following conditions are met:

[0085]

[0086] Among them, is the information contained in the cross-spectrum matrix obtained in the (j + 1)-th iteration.

[0087] The sparse deconvolution method can effectively suppress speckle artifacts and improve the signal-to-noise ratio. The CLEAN algorithm is a two-dimensional orthogonal matching pursuit (OMP), and its function is similar to the sparse deconvolution of the l1 norm. The CLEAN algorithm assumes that the source consists of scattering points. It finds the source with the largest amplitude in each iteration and uses PSF deconvolution to remove the contribution of this point to the imaging result.

[0088] As Figure 3 shown, the beamformed image generated by the convolution between the real image (obtained by imaging the cross-spectrum matrix) and the PSF can be regarded as the superposition of the PSFs at the corresponding source positions. The beamformed image is defined as a "dirty" image, in which three sources are almost inseparable. Then, the point with the largest amplitude is extracted into the "clean" image, and the corresponding PSF of this point is subtracted from the "dirty" image. The above steps are repeated until the real image (i.e., the "clean" image) is reconstructed convergently.

[0089] During implementation, the effect of the PSF is removed by updating the CMS:

[0090]

[0091] Among them, is the CMS at the j-th iteration; γ is a control parameter, which determines the magnitude of the removed power; is the maximum value of the pixel intensity in the j-th "dirty" image (the dirty image corresponding to the j-th iteration), and are the horizontal and vertical coordinates of the peak point (the point with the largest amplitude), respectively.

[0092] Then, the "dirty" image I at the (j + 1)-th iteration can be obtained j+1 (x, z, ω):

[0093]

[0094] Furthermore, for the "clean" image I clean The peak point information is collected as follows:

[0095]

[0096] When the CSM obtained in the current iteration contains more information than the CSM obtained in the previous iteration, the iteration terminates. This termination iteration condition is written as the following formula:

[0097]

[0098] In the formula, is the information contained in the CSM obtained in the j-th iteration.

[0099] Preferably, in step (3), the calculation formula for superimposing the imaging results corresponding to each frequency is as follows:

[0100] I(x, z) = ∫I clean (x, z, ω)dω

[0101] where I(x, z) is the final imaging result; I clean (x, z, ω) is the imaging result corresponding to each frequency; ω is the frequency; dω is the discrete step size in the frequency domain.

[0102] The above imaging method uses a frequency-by-frequency iterative method to calculate the sound pressure fields corresponding to each discrete frequency within the selected frequency band range [ω min , ω max , where the frequency band range is determined by the full matrix data, and the discrete step size dω in the frequency domain is determined by the sampling frequency.

[0103] Preferably, in step (1), the direct wave is filtered from the synthetic aperture data first, and then migration is performed. The method of filtering the direct wave can be any filtering method. Here, a feasible method is given, and a window function can be used to select the effective signal section.

[0104] Preferably, in step (1), when migrating the synthetic aperture data to the upper surface of the target layer, the synthetic aperture data is first converted to the frequency-wavenumber domain, extrapolated, and then the inverse Fourier transform is performed to obtain the target layer data in the time-space domain;

[0105] In the frequency-wavenumber domain, the synthetic aperture data is migrated according to the following formula:

[0106]

[0107] Among them, is the target layer data; is the synthetic aperture data in the frequency-wavenumber domain; z Tl-1 is the depth of the upper surface of the target layer; ω is the frequency; is the wavenumber of the excitation element in the horizontal direction; is the wavenumber of the receiving element in the horizontal direction; d j is the thickness of the j-th layer; is the wavenumber of the j-th layer in the vertical direction; i is the imaginary unit.

[0108] The layered structure in inhomogeneous media can cause severe acoustic refraction and reverberation, and deteriorate the image quality. Therefore, a wavefield extrapolation method for STA data has been developed to improve the imaging quality of imaging methods in layered structure media.

[0109] The measured data D(x s , x r , t) can be regarded as the wavefield on the surface. The excitation element and the receiving element each occupy two dimensions in space. Therefore, the wavefield of STA data has five dimensions, denoted as P(x s , z s , x r , z r , t), where x s and z s are the horizontal and vertical coordinates of the excitation element respectively, and x r and z r are the horizontal and vertical coordinates of the receiving element respectively.

[0110] The wavefield on the surface is expressed as:

[0111] P(x s , 0, x r , 0, t) = D(x s , x r , t) (27)

[0112] The wavefield on the surface is obtained through a three-dimensional Fourier transform to get the wavefield in the frequency-wavenumber domain as:

[0113]

[0114] Among them, and are the horizontal wavenumbers of the excitation element and the receiving element respectively; t is the time.

[0115] The downward migration process of the surface wavefield is as Figure 4As shown in (a) and (b), the imaging target is located in the second medium. To perform wavenumber shaping, the wave field at the upper interface z1 of the second layer must be obtained.

[0116] First, the received wave field is migrated according to the following formula:

[0117]

[0118] where d1 is the thickness of the first layer; is the vertical wavenumber of the receiving element, defined as:

[0119]

[0120] In the formula, c1 is the sound speed of the first layer; ω is the frequency.

[0121] After migration, the interface diagrams of the receiving and transmitting elements are as shown in Figure 4 (c) and (d), where the receiving element is migrated to the interface z1, while the excitation element remains on the measurement surface.

[0122] Furthermore, the excitation element is migrated to the interface z1 according to the following formula:

[0123]

[0124] where is the vertical wavenumber of the excitation element, defined as:

[0125]

[0126] At this time, as shown in the cross-sectional diagrams in Figure 4 (e) and (f), both the excitation and receiving elements are located at the interface.

[0127] Generally speaking, the wave field on the surface can be directly migrated to the depth z1:

[0128]

[0129] where k z is the total wavenumber for downward migration, defined as:

[0130]

[0131] Then, the wave field at the depth z1 is converted back to the spatial domain to obtain the equivalent STA data measured at this depth:

[0132]

[0133] which can be used to image the inhomogeneities in the second layer.

[0134] Compared with the prior art, the beneficial effects of the present invention are as follows:

[0135] The high signal-to-noise ratio imaging method for synthetic aperture data of the present invention uses singular value decomposition (SVD) to filter out most of the noise in the STA data to extract the signal subspace, and transforms the passive heterogeneity detection problem into an active source localization problem. The cross-spectral matrix is calculated based on the extracted signal subspace and used for beamforming. At the same time, the cross-spectral matrix is deconvolved using the unsteady point spread function (PSF) to further suppress speckle artifacts, greatly improving the signal-to-noise ratio and imaging quality. In addition, phase shift is introduced to enable the imaging method to meet the imaging requirements of layered media and has wide applicability. BRIEF DESCRIPTION OF THE DRAWINGS

[0136] Figure 1 is a flowchart of the imaging method according to an embodiment of the present invention;

[0137] Figure 2 where: (a) is a schematic diagram of wave propagation during synthetic aperture data acquisition; (b) is a schematic diagram of wave propagation when the signal subspace is equivalent to the active sound source localization problem;

[0138] Figure 3 is to perform PSF deconvolution on the beamformed image to reconstruct the real image;

[0139] Figure 4 where: (a) is a cross-sectional schematic diagram of the STA data corresponding to the same excitation element; (b) is a cross-sectional schematic diagram of the STA data corresponding to the same receiving element; (c) is a cross-sectional schematic diagram of the STA data corresponding to the same excitation element after shifting the receiving element downward to depth z1; (d) is a cross-sectional schematic diagram of the STA data corresponding to the same receiving element after shifting the receiving element downward to depth z1; (e) is a cross-sectional schematic diagram of the STA data corresponding to the same excitation element after shifting both the receiving and transmitting elements downward to depth z1; (f) is a cross-sectional schematic diagram of the STA data corresponding to the same receiving element after shifting both the receiving and transmitting elements downward to depth z1;

[0140] Figure 5 where: (a) is a schematic diagram of direct measurement of pork; (b) is the STA data directly measured in (a); (c) is a schematic diagram of assisted measurement of pork with a wedge; (d) is the STA data obtained by wedge-assisted measurement in (b); Results of image reconstruction of the direct measurement data in (a) by different methods: (e) is DAS, (f) is PSM, (g) is RTM, (h) is the method in the embodiment of the present invention; Results of image reconstruction of the wedge-assisted measurement data in (b) by different methods: (i) is DAS, (j) is PSM, (k) is RTM, (l) is the method in the embodiment of the present invention. DETAILED DESCRIPTION OF THE INVENTION

[0141] As Figure 1 shown, a high signal-to-noise ratio imaging method for synthetic aperture data (STA) includes the following steps:

[0142] (1) Preprocessing: Input the synthetic aperture data D(x s , x r , t), and use a window function to select the effective signal section and filter out the direct wave to obtain

[0143] (2) First, perform a three-dimensional Fourier transform on the preprocessed synthetic aperture data to the frequency-wavenumber domain to obtain Shift through the following formula to obtain the data of the target layer (Tl layer) in the frequency-wavenumber domain The derivation process is as follows:

[0144] The layered structure in the inhomogeneous medium will cause serious acoustic refraction and reverberation, and deteriorate the image quality. Therefore, a wavefield extrapolation method for STA data has been developed to improve the imaging quality of the imaging method in the layered structure medium.

[0145] The measured data D(x s , x r , t) can be regarded as the wavefield on the surface. The excitation element and the receiving element each occupy two dimensions in space. Therefore, the wavefield of STA data has five dimensions, denoted as P(x s , z s , x r , z r , t), where x s and z s are the horizontal and vertical coordinates of the excitation element respectively, and x r and z r are the horizontal and vertical coordinates of the receiving element respectively.

[0146] The wavefield on the surface is expressed as:

[0147] P(x s , 0, x r , 0, t) = D(x s , x r , t) (27)

[0148] The wavefield on the surface is transformed through a three-dimensional Fourier transform to obtain the wavefield in the frequency-wavenumber domain as:

[0149]

[0150] Among them, and are the horizontal wavenumbers of the excitation element and the receiving element, respectively; t is the time.

[0151] The downward migration process of the surface wave field is as Figure 4 shown in (a) and (b). It is assumed that the imaging target is located in the second layer of the medium. To perform wavenumber shaping, the wave field at the upper interface z1 of the second layer must be obtained.

[0152] First, the received wave field is migrated according to the following formula:

[0153]

[0154] where d1 is the thickness of the first layer; is the vertical wavenumber of the receiving element, defined as:

[0155]

[0156] In the formula, c1 is the sound speed of the first layer.

[0157] After migration, the interface diagrams of the receiving element and the transmitting element are as Figure 4 shown in (c) and (d), where the receiving element is migrated to the interface z1, while the excitation element remains on the measurement surface.

[0158] Furthermore, the excitation element is migrated to the interface z1 according to the following formula:

[0159]

[0160] where is the vertical wavenumber of the excitation element, defined as:

[0161]

[0162] At this time, as Figure 4 shown in the cross-sectional diagrams in (e) and (f), both the excitation and receiving elements are located at the interface.

[0163] Generally speaking, the wave field on the surface can be directly migrated to the depth z1:

[0164]

[0165] where k z is the total wavenumber of downward migration, defined as:

[0166]

[0167] Then, the wave field at the depth z1 is converted back to the spatial domain to obtain the equivalent STA data measured at this depth:

[0168]

[0169] It can be used for imaging heterogeneities in the second layer.

[0170] From the above derivation, the extrapolation formula in the frequency-wavenumber domain is as follows:

[0171]

[0172] In the formula, and are the horizontal wavenumbers of the excitation element and the receiving element respectively; z Tl-1 is the depth of the upper surface of the target layer; ω is the frequency; d j is the thickness of the j-th layer; is the wavenumber of the j-th layer in the vertical direction; i is the imaginary unit.

[0173] After that, is subjected to an inverse Fourier transform to obtain the synthetic aperture data (target layer data) of the target layer

[0174]

[0175] (3) Perform a singular value decomposition (SVD) on to extract the signal subspace.

[0176] The synthetic transmit aperture (STA) uses each element in the phased array to sequentially transmit spherical waves to the target area, and at the same time all elements in the phased array receive the reflected waves from the medium. For a linear array of N elements, the acquired signal obtained in this acquisition mode contains three dimensions, denoted here as where N t is the number of sampling points in the time domain. Therefore, the STA data is a response matrix that contains the response relationships between all elements. As shown in Figure 2 (a), the wave is emitted from the element located at R s =(x s ,0), propagates in the medium, reflects after encountering a reflector, and the reflected wave is received by the element located at R r =(x r ,0).

[0177] The propagation of the acoustic wave between the transmit element R s and the receive element R r can be modeled in the frequency domain as:

[0178]

[0179] where M is the number of reflectors, and the reflector is located at R m =(x m ,zm ) and the reflection coefficient is ε m ; F(ω) is the spectrum of the excitation signal; D s (θ s , ω) and D t (θ r , ω) are the directivity coefficients of the transmitting element and the receiving element respectively; θ s and θ r are the azimuth angles of the transmitting element and the receiving element respectively; ω is the frequency.

[0180] Among them, the directivity coefficient is defined as:

[0181]

[0182] Among them, a is the element width and c is the speed of sound.

[0183] G is the Green's equation, obtained from the following formula:

[0184]

[0185] Among them, H0 is the zero-order Hankel equation of the first kind; i is the imaginary unit; the Green's equation G(R m , R r , ω) represents the response to a point source at R m at point R r .

[0186] In the frequency domain, the N×N response matrix of all elements can be written as:

[0187]

[0188] Among them:

[0189] g m = [D1(θ1, ω)G(R1, R m , ω), D2(θ2, ω)G(R2, R m , ω), …, D N (θ1, ω)G(R N , R m , ω)] T (5)

[0191] Among them, H is the transpose and complex conjugate; T is the transpose.

[0192] G is an N×M matrix, expressed as:

[0193] G = [g1, g2, …, g M (6)

[0194] M is the number of reflectors; N is the number of elements.

[0195] Since the direction coefficients in different orientations corresponding to the imaging region do not differ much, the direction coefficient of the vibration element can be ignored, and thus G is an orthogonal matrix.

[0196] Q is an M×M dimensional matrix, written as:

[0197]

[0198] Therefore, the measured response matrix D(ω) can be modeled as:

[0199] D(ω) = [D(x s ,x r ,ω)] = GQG H +N(ω) (8)

[0200] where D(ω) is an N×N matrix, which is the set of all receive-transmit vibration element pairs in D(x s ,x r ,ω) at frequency ω, and D(x s ,x r ,ω) can be obtained by performing a one-dimensional Fourier transform on D(x s ,x r ,t) in the time dimension.

[0201] Here it is assumed that the number of reflection points is less than the number of vibration elements, that is, M < N. Therefore, P(ω) is an ideal response matrix of rank M, but this ideal response matrix becomes D(ω) after being contaminated by noise, and the latter is full rank. In theory, the first M eigenvalues of D(ω) are much larger than the remaining N - M ones.

[0202] The measured response matrix can be written in the following eigenvalue decomposition form:

[0203]

[0204] where the subscripts sig and n represent the signal subspace and the noise subspace respectively, and the signal subspace and the noise subspace are distinguished according to the size of the eigenvalues; ∑ sig (ω), ∑ n (ω) are the eigenvalues of the signal subspace and the noise subspace respectively, both of which are M×M dimensional diagonal matrices corresponding to M reflection points; U sig (ω) and V sig (ω) are two sets of eigenvectors of the signal subspace, which are N×M dimensional matrices. In theory, they should be the same, however, the noise in actual measurement will cause differences between them. Here, V sig (ω) is ignored, and U sig (ω) is used as the signal subspace for subsequent cross-spectrum matrix calculations; Un (ω) and V n (ω) are two groups of eigenvectors of the noise subspace, and are also N×M-dimensional matrices.

[0205] Therefore, SVD can separate the scattered signal from the noise, and the following formula can be obtained:

[0206]

[0207] Comparing the three terms in formula (10), V sig (ω), ∑ sig (ω) and U sig (ω) can be interpreted as the signal emitted by the incident wave from the excitation element, the reflection coefficient of the reflector, and the received signal of the reflected wave, respectively.

[0208] Therefore, U sig (ω) = [u sig1 , …, u sigM can be regarded as the signal of the wave actively emitted from the scattering point received by the sensor, as shown in Figure 2 (b). Among them, u sigm is the signal vector in the signal subspace U sig (ω), m ∈ [1, M], and M is the number of reflectors.

[0209] Therefore, extracting U sig (ω) by SVD not only filters out most of the noise, but also transforms the passive heterogeneous detection problem into an active sound source localization problem.

[0210] (4) Calculate the cross-spectral matrix at each frequency according to the extracted signal subspace.

[0211] Beamforming is an important method for active sound source localization. The wave equation of sound waves in a two-dimensional homogeneous medium can be written as:

[0212]

[0213] Among them, p(x, z, t) is the sound pressure field at the point (x, z), q(t) is the point source located at (x0, z0), and the wave field in the time domain is:

[0214]

[0215] Among them, the distance is The time delay is

[0216] Therefore, the signal received by the phased array element is the wave field p(x r, 0, t), the sound source distribution can be reconstructed by time reversal. The possible sound source at the point (x, z) can be estimated by the following formula:

[0217]

[0218] where (x j , 0) is the position of the nth transducer element; t j = r j / c.

[0219] After performing the Fourier transform on , the expression of beamforming in the frequency domain is obtained as follows:

[0220]

[0221] where h(x, z, ω) and Y(x, z, ω) are N×1 vectors, written as:

[0222] h(x, z, ω) = [2πr1exp(iωt1), …, 2πr N exp(iωt N )] T (15)

[0223] Y(x, z, ω) = [P(x1, 0, ω), …, P(x N , 0, ω)] T (16)

[0224] is a complex number, and its energy can be calculated by the following formula:

[0225]

[0226] where the superscript * represents the conjugate; C(ω) is an N×N matrix containing the measured signals, C(ω) is the cross-spectral matrix (CMS), and any element C ij at the coordinate (i, j) can also be obtained by the following formula:

[0227] C ij = P(x i , 0, ω)P * (x j , 0, ω) (18)

[0228] For STA data, U sig (ω) is extracted by the SVD method and is regarded as the signal received by the phased array after the wave actively emitted from the heterogeneous body (reflector) at time t = 0. Here, the imaging condition is defined as follows:

[0229] I(x,z) = ∫I(x,z,ω)dω = ∫h H (x,z,ω)C M (ω)h(x,z,ω)dω (19)

[0230] Where:

[0231]

[0232] Among them, C M (ω) is the cross-power matrix at frequency ω; u sigm is the signal vector in the signal subspace U sig (ω).

[0233] According to the extracted U sig (ω), first calculate the cross-power matrix C M (ω), and then use the imaging condition and C M (ω) for imaging to obtain the imaging result corresponding to frequency ω.

[0234] (5) Use the non-stationary point spread function (PSF) to deconvolve the cross-spectrum matrix obtained at each frequency to obtain the high signal-to-noise ratio imaging result corresponding to each frequency.

[0235] The imaging of a phased array on a point source is called the point spread function (PSF), which is related to hardware parameters such as element size, element spacing, and element reception performance. The imaging result of the phased array is the convolution result of the PSF and the real image.

[0236] For a point source λ = (x λ , z λ ), the sound source introduces a CSM as:

[0237] C λ (ω) = Y(x λ , z λ , ω)Y H (x λ , z λ , ω)

[0238] Therefore, the point spread function PSF can be written as:

[0239] I λ (x,z,ω) = h H (x,z,ω)C λ (ω)h(x,z,ω)

[0240] Among them, the PSF is non-stationary because it is a function of the light source position. The PSF blurs the image and makes adjacent scattering points indistinguishable. Therefore, it is necessary to use deconvolution to eliminate the influence of the PSF and improve the resolution.

[0241] The sparse deconvolution method can effectively suppress speckle artifacts and improve the signal-to-noise ratio. The CLEAN algorithm is a two-dimensional orthogonal matching pursuit (OMP), and its function is similar to the sparse deconvolution of the l1 norm. The CLEAN algorithm assumes that the source consists of scattering points. It finds the source with the largest amplitude at each iteration and uses PSF deconvolution to remove the contribution of this point to the imaging result.

[0242] As Figure 3 shown, the beamformed image generated by the convolution between the real image (obtained by cross-spectrum matrix imaging) and the PSF can be regarded as the superposition of the PSFs at the corresponding source positions. The beamformed image is defined as a "dirty" image, in which the three sources are almost inseparable. Then, the point with the largest amplitude is extracted into the "clean" image, and the corresponding PSF of this point is subtracted from the "dirty" image. The above steps are repeated until the real image (i.e., the "clean" image) is reconstructed convergently.

[0243] During implementation, the effect of the PSF is removed by updating the CMS:

[0244]

[0245] where is the CMS at the j-th iteration; γ is a control parameter that determines the magnitude of the removed power; is the maximum value of the pixel intensity in the j-th "dirty" image (the dirty image corresponding to the j-th iteration), and are the horizontal and vertical coordinates of the peak point (the point with the largest amplitude), respectively.

[0246] Then, the "dirty" image I j+1 (x, z, ω) at the (j + 1)-th iteration can be obtained:

[0247]

[0248] Furthermore, the "clean" image I clean collects the peak point information:

[0249]

[0250] When the CSM obtained in the current iteration contains more information than the CSM obtained in the previous iteration, the iteration terminates. This iteration termination condition is written as the following formula:

[0251]

[0252] where is the information contained in the CSM obtained at the j-th iteration.

[0253] (6)Superimpose the imaging results corresponding to each frequency to obtain the final imaging result.

[0254] The superimposing calculation formula is as follows:

[0255] I(x,z) = ∫I clean (x,z,ω)dω

[0256] where I(x,z) is the final imaging result; I clean (x,z,ω) is the imaging result corresponding to each frequency; dω is the discrete step size in the frequency domain.

[0257] Detection experiment

[0258] In vitro experiments were conducted to evaluate the proposed method. Pork loin was used as the basic tissue, and three carbon rods (Φ = 1mm) were inserted into the tissue as targets (reflectors). Pork has weak heterogeneity, and its overall sound speed was calibrated to 1650 m / s. First, as Figure 5 shown in (a) below, a 64-element array transducer (Shantou Ultrasonic) with a center frequency of 2.25 MHz and an element interval of 0.75 mm was directly attached to the surface of the pork, and 64 / 64 OEM-PA (AOS.Ltd, America) was used to collect Figure 5 the STA dataset in (b) below, where the sampling frequency was 50 MHz and the time range was 32 μs. In addition, a wedge (SC63-NL-Z20) with a sound speed of 2337 m / s was installed on the phased array to assist in the measurement, as Figure 5 shown in (c) below. Then, a second STA dataset was collected, and the time range was increased to 42 μs, as Figure 5 shown in (d) below.

[0259] To quantitatively analyze the imaging results, the signal-to-noise ratio index SNR was first introduced. It reflects the contrast between the measured target area (D f ) and the background area (N - D f ) in the image, and it can quantitatively characterize the signal-to-noise ratio. The calculation formula is as follows:

[0260]

[0261] where, and are the pixel amplitudes of the target area and the background area respectively, and mean and std are the mean value and the standard deviation respectively.

[0262] These two STA datasets were processed by four imaging methods, namely DAS, PSM, RTM, and the method of this embodiment. The imaging results are as Figure 5 shown in (e)-(l) below. In the directly measured imaging results, onlyFigure 5 In (e), the target on the far left cannot be recognized in the image reconstructed by DAS. The background artifacts in the imaging results of PSM and RTM are very serious, showing as long strips and spots in Figure 5 (f) and (g). In Figure 5 (h), the method of this embodiment can clearly reconstruct three targets and suppress most of the artifacts.

[0263] From the imaging results of the wedge-assisted measurement, Figure 5 the target on the far left in (i) still cannot be recognized by DAS. As shown in Figure 5 (j) and (k), the spot artifacts in the imaging results of PSM and RTM are suppressed to a certain extent, but the stripe artifacts still contaminate the image. Figure 5 (l), the method of this embodiment reconstructs a clear image of the target after suppressing all artifacts.

[0264] Table 1 Figure 5 Signal-to-noise ratio and calculation time corresponding to the imaging results in

[0265]

[0266] Table 1 lists the signal-to-noise ratio and calculation time of different imaging methods in the experiment. Except for RTM, the signal-to-noise ratio of the results of all methods in the wedge-assisted measurement is higher than that in the direct measurement. Therefore, the wedge is helpful in this case. The signal-to-noise ratio of the imaging results of the method of this embodiment is much higher than that of other methods. Although the signal-to-noise ratio of the image reconstructed by PSM is the lowest, its calculation time is much lower than that of other methods, while RTM requires too much calculation time. In addition, the method of this embodiment requires much less time than DAS. Therefore, the imaging method of this embodiment not only improves the signal-to-noise ratio of the imaging results, but also is very competitive in terms of efficiency.

Claims

1. A high signal-to-noise ratio imaging method for synthetic aperture data, characterized in that, It includes the following steps: (1) Input the synthetic aperture data, shift the synthetic aperture data to the upper surface of the target layer to obtain the target layer data; (2) Perform singular value decomposition on the target layer data, extract the signal subspace, and calculate the cross-spectrum matrix for each frequency; (3) Perform deconvolution on the cross-spectrum matrix for each frequency to obtain the imaging result corresponding to each frequency, and superimpose the imaging results corresponding to each frequency to obtain the final imaging result; In step (1), when shifting the synthetic aperture data to the upper surface of the target layer, first convert the synthetic aperture data to the frequency-wavenumber domain, perform extrapolation, and then perform inverse Fourier transform to obtain the target layer data in the time-space domain; In step (2), the cross-spectrum matrix C M (ω) is calculated as follows: where M is the number of reflection points; ω is the frequency; u sigm is the signal vector in the signal subspace; the superscript H represents transpose and complex conjugate; In step (3), use the non-steady state point spread function to perform deconvolution on the cross-spectrum matrix obtained for each frequency; The calculation formula for superimposing the imaging results corresponding to each frequency is as follows: I(x,z) = ∫I clean (x,z,ω)dω Among them, I(x,z) is the final imaging result; I clean (x,z,ω) is the imaging result corresponding to each frequency; ω is the frequency.

2. The high signal-to-noise ratio imaging method for synthetic aperture data according to claim 1, characterized in that In step (2), the form of performing singular value decomposition on the target layer data is: where \(D(\omega)\) is the data of the target layer at frequency \(\omega\); the subscripts sig and n represent the signal subspace and the noise subspace respectively; \(\omega\) is the frequency; \(U\) sig \((\omega)\), \(V\) sig \((\omega)\) are two groups of eigenvectors of the signal subspace respectively; \(\sum\) sig \((\omega)\) is the eigenvalue of the signal subspace; \(U\) n \((\omega)\), \(V\) n \((\omega)\) are two groups of eigenvectors of the noise subspace respectively; \(\sum\) n \((\omega)\) is the eigenvalue of the noise subspace; the superscript \(H\) represents the transpose and complex conjugate.

3. The high signal-to-noise ratio imaging method for synthetic aperture data according to claim 1, wherein In step (3), when performing deconvolution on the cross-spectrum matrix for each frequency, terminate the iteration when the following conditions are met: Among them, is the information contained in the cross-spectrum matrix obtained in the (j + 1)-th iteration; is the cross-spectrum matrix obtained in the (j + 1)-th iteration; is the cross-spectrum matrix obtained in the j-th iteration; M is the number of reflection points; ω is the frequency.

4. The high signal-to-noise ratio imaging method for synthetic aperture data according to claim 1, wherein In step (1), first filter the direct wave from the synthetic aperture data and then perform the shift.

5. The high signal-to-noise ratio imaging method for synthetic aperture data according to claim 1, wherein In step (1), in the frequency-wavenumber domain, shift the synthetic aperture data according to the following formula: Among them, is the data of the target layer; is the synthetic aperture data in the frequency-wavenumber domain; z Tl-1 is the depth of the upper surface of the target layer; ω is the frequency; is the wavenumber of the excitation element in the horizontal direction; is the wavenumber of the receiving element in the horizontal direction; d j is the thickness of the j-th layer; is the wavenumber of the j-th layer in the vertical direction; i is the imaginary unit; ∏ represents the product; Tl is the target layer; e is the natural constant.

Images