An inverse projection method based on acoustic reflection tomography imaging

Abstract

This invention belongs to the technical field of inverse projection methods, specifically relating to an inverse projection method based on acoustic reflection tomography imaging. It utilizes upsampling and frequency-domain matched filtering techniques to solve the problem of unknown cross-sectional distance indices at various angles in the inverse projection algorithm, aiming to improve the accuracy of cross-sectional distance localization, thereby enhancing the precision and noise resistance of acoustic reflection tomography sonar imaging. This invention fully leverages the high peak value and strong anti-interference and noise capabilities of matched filtering, as well as the one-to-one temporal correspondence between the acoustic reflection tomography cross-sectional distance and the matched filter peak value to estimate the cross-sectional distance index. Furthermore, this invention uses frequency-domain matched filtering to transform the time-domain convolution operation of the received signal and the reference signal into a frequency-domain multiplication, reducing computational complexity. This invention employs frequency-domain zero-padding to upsample the time-domain signal, improving the accuracy of the peak position after matched filtering.

Description

Technical Field

[0001] This invention belongs to the technical field of inverse projection methods, specifically relating to an inverse projection method based on acoustic reflection tomography imaging. Background Technology

[0002] Existing high-resolution underwater target imaging methods utilize sonar movement along a straight line, i.e., synthetic aperture, to improve the resolution of small targets. Strip synthetic aperture sonar (SAS) systems have an illumination angle of only ±45°, primarily limited by the divergence angle of the sonar radiation beam pattern. Focused SAS, on the other hand, illuminates the target area for a longer time than strip SAS by scanning the beam, thus achieving a longer synthetic aperture length and ultimately improving azimuth beam resolution. However, the increase in synthetic aperture length is limited by the observation angle between the sonar carrier and the target. In the case of a straight path, continuously increasing this observation angle is difficult. Moreover, as the sonar moves away from the target, the beam incidence angle continuously increases, significantly weakening target coherence. Simply increasing the synthetic aperture length in a straight direction reduces the effectiveness of improving target resolution and noise immunity.

[0003] X-ray computed tomography-assisted imaging (CT) uses 360° illumination along a cross-section of the target to obtain a projected signal after attenuation. The received signal is a function of the distance and angle of the cross-section. Through inverse projection reconstruction, the spatial distribution of the target's attenuation coefficient can be obtained. In contrast, underwater small-target sonar imaging based on acoustic reflection tomography records projection information provided by the target's back-reflection signal. Accurate estimation of the distance trajectory along the acoustic reflection tomography cross-section is necessary to reconstruct the target's spatial acoustic reflection distribution function using inverse projection algorithms. Summary of the Invention

[0004] To address the technical problem mentioned above regarding underwater small target sonar imaging based on acoustic reflection tomography, which requires accurate estimation of the distance trajectory of the acoustic reflection tomography cross section before the inverse projection algorithm can be used to reconstruct the spatial acoustic reflection distribution function of the target, this invention provides an inverse projection method based on acoustic reflection tomography imaging. This method utilizes upsampling and frequency domain matched filtering techniques to solve the problem of unknown cross section distance indices at various angles in the inverse projection algorithm. The aim is to improve the accuracy of cross section distance positioning, thereby enhancing the precision and noise resistance of acoustic reflection tomography sonar imaging.

[0005] To solve the above-mentioned technical problems, the technical solution adopted by the present invention is as follows:

[0006] A method for inverse projection based on acoustic reflection tomography includes the following steps:

[0007] S1. Set the initial scanning angle θ i , i = 1;

[0008] S2. The acoustic reflection time-domain echo signal s(t,θ) of the target at the initial scanning angle. i Baseband conversion and low-pass filtering are performed, and a two-dimensional function S(k,θ) with respect to wavenumber and angle is obtained through Fourier transform. i );

[0009] S3. Set the reference position as the origin O(0,0) of the target area center, and calculate the acoustic reflection time-domain echo signal s0(t,θ) of the cross-section at the reference position under the same scanning angle. i The signal is then converted to baseband and subjected to Fourier transform to obtain a two-dimensional reference signal S0(k,θ) with respect to wavenumber and angle. i );

[0010] S4. Perform frequency domain matched filtering on the frequency domain data obtained in S2 and S3.

[0011] S5. Update the scanning angle i+1, and repeat S2, S3, and S4 until all angle scans are completed, obtaining the received data matrix S after wavenumber domain matched filtering at each angle. Match The dimension is N×M;

[0012] S6. Calculate the wavenumber domain interval and determine the grid point position of the target area based on the transmitted signal bandwidth and scanning angle range. Determine the reference time after matched filtering based on the reference position. Determine the new sampling interval based on the upsampling rate U.

[0013] S7. Set the initial grid point position (x p ,y j ), p=1, j=1, set the initial value of the target spatial distribution function to be reconstructed at this location SAS(p,j) to 0;

[0014] S8. Set the initial scanning angle θ i , i = 1;

[0015] S9. Take the frequency domain data at this angle calculated by S5, and pad it with 0 at the beginning and end of the data. The length is the upsampling rate U minus 1 and then multiplied by half of the original data length.

[0016] S10. Perform an inverse Fourier transform on the frequency domain data padded with zeros in S10 to obtain the upsampled time domain signal s. u (t,θ i );

[0017] S11. Calculate the time delay of the received signal at the current grid position and the current scanning angle, and subtract it from the reference time after matched filtering determined in S6. Determine the index corresponding to the cross-sectional distance based on the new sampling interval.

[0018] S12. Accumulate the value at this index of the upsampled matched filter time-domain signal obtained in S10 into the variable set in S7;

[0019] S13. Increment the scanning angle and repeat S9 to S12 until the cross-sectional distance corresponding index is estimated for all scanning angles and the values ​​at the corresponding indexes of the upsampled matched filter time domain signal for all scanning angles are accumulated.

[0020] S14. Update the grid point positions, and repeat S8 to S13 until the reconstruction of the target spatial distribution function at all locations is completed.

[0021] The wavenumber domain data obtained in S4 is subjected to wavenumber domain matched filtering.

[0022] In step S6, the wavenumber domain interval is calculated based on the transmitted signal bandwidth B and the scanning angle range Θ, and the grid point spacing Δx, Δy and the number of points N in the target area are determined using the following formula. x N y The reference time after matched filtering is determined based on the reference cross-section position. c is the speed of sound in water, and the new sampling interval is determined based on the upsampling rate U.

[0023] k x =2kcosθ,k y =2ksinθ

[0024]

[0025] Δx=2π / (max(k x )-min(k x ),Δy=2π / (max(k y )-min(k y ))

[0026]

[0027] R0 is the center distance of the target area, θ is the scanning angle, Θ∈[-π,π], Δθ is the scanning interval, N is the number of frequency sampling points, and k x The matrix representing the component matrix of wavenumber k along the x-axis at various angles, wherein k y The Δk represents the component matrix of wavenumber k along the y-axis at various angles. x and Δk y These represent the uniform wavenumber sampling intervals of the reconstructed target in rectangular coordinates.

[0028] In S9, the wavenumber domain data at this angle calculated by S5 is taken as S(k,θ). iThe data is padded with zeros at the beginning and end, with a length of [length missing].

[0029] In step S11, the time delay τ of the received signal at the current grid position and current scanning angle is calculated and subtracted from the matched-filtered reference time determined in step S6. The index corresponding to the cross-sectional distance is then determined based on the new sampling interval. i ;

[0030]

[0031]

[0032] The x p Let be the coordinates of the p-th spatial grid point in the x-direction of a Cartesian coordinate system, and let y be the coordinates of the p-th spatial grid point in the x-direction. j Let τ be the coordinate of the j-th spatial grid point in the y-direction. ref The reference time after matched filtering, determined by the reference cross-section position, is... The Δt is the time-domain sampling interval. This is the new sampling interval after increasing the sampling rate U.

[0033] In S12, the value at that index of the upsampled matched filter time-domain signal obtained in S10 is accumulated into the variable set in S7;

[0034] SAS(p,j)=SAS(p,j)+s u (Index i ,θ i )

[0035] The s u (Index i ,θ i ) represents the value at this index of the upsampled matched-filtered time-domain signal.

[0036] In S13, the scanning angle θ is updated. i , i+1, repeat S9 to S12 until the estimation of the cross-sectional distance corresponding to all scanning angles i=1…M is completed, and the accumulation of the values ​​at the corresponding indexes of the matched filter time domain signals upsampled under all scanning angles is completed.

[0037] In S14, the grid point position (x) is updated. p ,y j ), p=1…N x ,j=1…N y Repeat steps S8 to S13 until the spatial distribution function of the target at all locations is reconstructed.

[0038] Compared with the prior art, the beneficial effects of this invention are:

[0039] This invention fully utilizes the high peak value and strong anti-interference and noise capabilities of matched filters, as well as the one-to-one temporal correspondence between the cross-sectional distance in acoustic reflection tomography and the peak value of the matched filter to estimate the cross-sectional distance index. Furthermore, this invention employs frequency-domain matched filtering, transforming the time-domain convolution operation between the received signal and the reference signal into frequency-domain multiplication, thus reducing computational complexity. This invention also uses frequency-domain zero-padding to upsample the time-domain signal, improving the accuracy of the peak position after matched filtering. Attached Figure Description

[0040] To more clearly illustrate the embodiments of the present invention or the technical solutions in the prior art, the drawings used in the description of the embodiments or the prior art will be briefly introduced below. Obviously, the drawings in the following description are merely exemplary, and those skilled in the art can derive other embodiments based on the provided drawings without creative effort.

[0041] The structures, proportions, sizes, etc. illustrated in this specification are only for the purpose of assisting those skilled in the art in understanding and reading the content disclosed herein, and are not intended to limit the conditions under which the present invention can be implemented. Therefore, they have no substantial technical significance. Any modifications to the structure, changes in the proportions, or adjustments to the size, without affecting the effects and objectives that the present invention can produce, should still fall within the scope of the technical content disclosed in the present invention.

[0042] Figure 1 This is a working model diagram of the acoustic reflection tomography imaging sonar of the present invention;

[0043] Figure 2 This is a reconstruction error curve after upsampling when the signal-to-noise ratio is -35dB according to the present invention;

[0044] Figure 3 This is a graph showing the change in the average reconstruction error as a function of the signal-to-noise ratio in this invention.

[0045] Figure 4 This is a flowchart of the method of the present invention. Detailed Implementation

[0046] To make the objectives, technical solutions, and advantages of the embodiments of the present invention clearer, the technical solutions in the embodiments of the present invention will be clearly and completely described below. Obviously, the described embodiments are only some embodiments of this application, and not all embodiments. These descriptions are only for further illustrating the features and advantages of the present invention, and not for limiting the claims of the present invention. Based on the embodiments in this application, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of this application.

[0047] The specific embodiments of the present invention will be described in further detail below with reference to the accompanying drawings and examples. The following examples are for illustrative purposes only and are not intended to limit the scope of the invention.

[0048] In this embodiment, as Figure 1 The diagram illustrates the application scenario of the algorithm in this embodiment. An underwater transceiver platform moves in a circular orbit around a target, with a distance R0 relative to the target area's center. The scanning angle is θ, ranging from Θ∈[-π,π], the scanning interval is Δθ, the number of angles is M, the transmitted signal is a linear frequency modulated pulse with pulse width T, a sampling time Δt, a bandwidth of B, and a center frequency f. c The number of frequency sampling points is N. Based on this working method, the inverse projection method based on acoustic reflection tomography imaging in this embodiment includes the following steps:

[0049] Step 1, set the initial scanning angle θ i , i = 1.

[0050] Step 2: Calculate the acoustic reflection time-domain echo signal s(t,θ) of the target at this angle. i Baseband conversion and low-pass filtering are performed, and a two-dimensional function S(k,θ) with respect to wavenumber and angle is obtained through Fourier transform. i ).

[0051] Step 3: Set the reference position as the origin O(0,0) of the target area center, and calculate the acoustic reflection time-domain echo signal s0(t,θ) of the cross-section at the reference position under the same scanning angle. i The signal is then converted to baseband and subjected to Fourier transform to obtain a two-dimensional reference signal S0(k,θ) with respect to wavenumber and angle. i ).

[0052] Step 4: Perform wavenumber domain matched filtering on the wavenumber domain data obtained in steps 2 and 3.

[0053] Step 5: Update the scanning angle i+1, and repeat steps 2, 3, and 4 until all angle scans are completed, obtaining the received data matrix S after wavenumber domain matched filtering at each angle. Match The dimension is N×M.

[0054] Step 6: Calculate the wavenumber domain interval based on the transmitted signal bandwidth B and the scanning angle range Θ, and determine the grid point spacing Δx, Δy and the number of points N in the target area using the following formula. x N y The reference time after matched filtering is determined based on the reference cross-section position. c is the speed of sound in water, and the new sampling interval is determined based on the upsampling rate U.

[0055] kx =2kcosθ,k y =2ksinθ

[0056]

[0057] Δx=2π / (max(k x )-min(k x ),Δy=2π / (max(k y )-min(k y ))

[0058]

[0059] Where: R0 is the center distance of the target area, θ is the scanning angle, Θ∈[-π,π], Δθ is the scanning interval, N is the number of frequency sampling points, and k x This represents the component matrix of wavenumber k along the x-axis at various angles, k y Δk represents the component matrix of wavenumber k along the y-axis at various angles. x and Δk y These represent the uniform wavenumber sampling intervals of the reconstructed target in rectangular coordinates.

[0060] Step 7, Set the initial grid point positions (x p ,y j ), p=1, j=1.

[0061] Step 8: Set the initial value of the target spatial distribution function to be reconstructed at this location, SAS(p,j), to 0.

[0062] Step 9, Set the initial scanning angle θ i , i = 1.

[0063] Step 10: Take the wavenumber domain data S(k,θ) at this angle obtained in Step 5. i The data is padded with zeros at the beginning and end, with a length of [length missing].

[0064] Step 11: Perform an inverse Fourier transform on the wavenumber domain data padded with zeros in Step 10 to obtain the upsampled time-domain signal s. u (t,θ i ).

[0065] Step 12: Calculate the received signal delay τ at the current grid position and current scanning angle according to the following formula, and subtract it from the matched-filtered reference time determined in Step 6. Determine the index corresponding to the cross-sectional distance based on the new sampling interval. i .

[0066]

[0067]

[0068] Where: x p Let be the coordinates of the p-th spatial grid point in the x-direction of a Cartesian coordinate system, and y be the coordinates of the grid point in the x-direction. j Let τ be the coordinate of the j-th spatial grid point in the y-direction. ref The reference time after matched filtering is determined by the reference cross-section position. Δt is the time-domain sampling interval. This is the new sampling interval after increasing the sampling rate U.

[0069] Step 13: Add the value at that index of the upsampled matched filter time-domain signal obtained in Step 11 to the variable set in Step 8.

[0070] SAS(p,j)=SAS(p,j)+s u (Index i ,θ i )

[0071] Step 14, update the scanning angle θ i , i+1, repeat steps 10 to 13 until the estimation of the cross-sectional distance corresponding to all scanning angles i=1…M is completed, and the accumulation of the values ​​at the corresponding indices of the matched filter time domain signals upsampled at all scanning angles is completed.

[0072] Step 15, update grid point positions (x p ,y j ), p=1…N x ,j=1…N y Repeat steps 8 to 14 until the spatial distribution function of the target at all locations is reconstructed.

[0073] In this embodiment, as Figure 2 As shown, when the target is a single point source target, the cross-sectional position is estimated using two processing methods: upsampling and frequency domain matched filtering, and no upsampling. The spatial distribution positioning error curve is reconstructed using the inverse projection algorithm at a signal-to-noise ratio of -35dB. Figure 3 As shown, when the signal-to-noise ratio (SNR) varies from -40dB to -10dB, 30 Monte Carlo simulations are performed at each SNR. Two methods are used to estimate the cross-sectional position: upsampling with frequency-domain matched filtering and no upsampling. The average error curve of the spatial distribution positioning is reconstructed using the inverse projection algorithm. It can be seen that the proposed method significantly improves reconstruction accuracy, noise immunity, and robustness compared to the method without upsampling. Figure 4 This is a flowchart of the inverse projection algorithm based on acoustic reflection tomography.

[0074] The above description only illustrates the preferred embodiments of the present invention. However, the present invention is not limited to the above embodiments. Within the scope of knowledge possessed by those skilled in the art, various changes can be made without departing from the spirit of the present invention, and all such changes should be included within the protection scope of the present invention.

Claims

1. A method for inverse projection based on acoustic reflection tomography imaging, characterized in that: Includes the following steps: S1. Set the initial scanning angle , ; S2. Calculate the acoustic reflection time-domain echo signal of the target at the initial scanning angle. Baseband conversion and low-pass filtering are performed, and a two-dimensional function of wavenumber and angle is obtained through Fourier transform. ; S3. Set the reference position as the origin of the target area's center. Calculate the acoustic reflection time-domain echo signal of the cross section at the reference position under the same scanning angle. Then, baseband conversion and Fourier transform are performed to obtain a two-dimensional reference signal with respect to wavenumber and angle. ; S4. Perform frequency domain matched filtering on the frequency domain data obtained in S2 and S3. S5, Update Scan Angle Repeat steps S2, S3, and S4 until all angle scans are completed, obtaining the received data matrix after wavenumber domain matched filtering at each angle. Dimension is N is the number of frequency sampling points. Number of scanning angles; S6. Calculate the wavenumber domain interval and determine the grid point position of the target area based on the transmitted signal bandwidth and scanning angle range. Determine the reference time after matched filtering based on the reference position. Determine the new sampling interval based on the upsampling rate U. S7. Set the initial grid point position , Set the initial value of the target spatial distribution function to be reconstructed at this location. =0; S8. Set the initial scanning angle , ; S9. Take the frequency domain data at this angle calculated by S5, and pad it with 0 at the beginning and end of the data. The length is the upsampling rate U minus 1 and then multiplied by half of the original data length. S10. Perform an inverse Fourier transform on the frequency domain data padded with zeros in S9 to obtain the upsampled time domain signal. ; S11. Calculate the time delay of the received signal at the current grid position and the current scanning angle, and subtract it from the reference time after matched filtering determined in S6. Determine the index corresponding to the cross-sectional distance based on the new sampling interval. S12. Accumulate the value at this index of the upsampled matched filter time-domain signal obtained in S10 into the variable set in S7; S13. Increment the scanning angle and repeat S9 to S12 until the cross-sectional distance corresponding index is estimated for all scanning angles and the values ​​at the corresponding indexes of the upsampled matched filter time domain signal for all scanning angles are accumulated. S14. Update the grid point positions, and repeat S8 to S13 until the reconstruction of the target spatial distribution function at all locations is completed.

2. The inverse projection method based on acoustic reflection tomography imaging according to claim 1, characterized in that: The wavenumber domain data obtained in S4 is subjected to wavenumber domain matched filtering. .

3. The inverse projection method based on acoustic reflection tomography imaging according to claim 1, characterized in that: In S6, the transmission signal bandwidth is used as a reference. and scanning angle range Calculate the wavenumber domain interval and determine the grid point spacing of the target region using the following formula. and points The reference time after matched filtering is determined based on the reference cross-section position. , To determine the speed of sound in water, a new sampling interval is determined based on the upsampling rate U. , The time-domain sampling interval; The The distance to the center of the target area, For the scanning angle, the The The scanning interval is N, where N is the number of frequency sampling points, and k is... x The matrix representing the component matrix of wavenumber k along the x-axis at various angles, wherein k y The Δk represents the component matrix of wavenumber k along the y-axis at various angles. x and Δk y These represent the uniform wavenumber sampling intervals of the reconstructed target in rectangular coordinates.

4. The inverse projection method based on acoustic reflection tomography imaging according to claim 1, characterized in that: The wavenumber domain data at this angle obtained from S5 in S9 is as follows: Pad the data with zeros at the beginning and end, with a length of [length missing]. .

5. The inverse projection method based on acoustic reflection tomography imaging according to claim 1, characterized in that: In step S11, the time delay of the received signal at the current grid position and the current scanning angle is calculated. The difference between this time and the reference time determined by the matched filter in S6 is used to determine the index corresponding to the cross-sectional distance based on the new sampling interval. ; The x p Let be the coordinates of the p-th spatial grid point in the x-direction of a Cartesian coordinate system, and let y be the coordinates of the p-th spatial grid point in the x-direction. j Let τ be the coordinate of the j-th spatial grid point in the y-direction. ref The reference time after matched filtering, determined by the reference cross-section position, is... , Distance from the center of the target area Where is the speed of sound in water, and N is the number of frequency sampling points. The time-domain sampling interval, the This is the new sampling interval after increasing the sampling rate U.

6. The inverse projection method based on acoustic reflection tomography imaging according to claim 1, characterized in that: In S12, the value at that index of the upsampled matched filter time-domain signal obtained in S10 is accumulated into the variable set in S7; The s u (Index i ,θ i ) represents the value at this index of the upsampled matched-filtered time-domain signal.

7. The inverse projection method based on acoustic reflection tomography imaging according to claim 1, characterized in that: The scanning angle is updated in S13. , Repeat steps S9 to S12 until all scanning angles are completed. Estimate the index corresponding to the lower cross section distance, and accumulate the values ​​at the corresponding indexes of the matched filtered time-domain signals sampled down from all scanning angles.

8. The inverse projection method based on acoustic reflection tomography imaging according to claim 1, characterized in that: In S14, the grid point positions are updated. , Repeat steps S8 to S13 until the spatial distribution function of the target at all locations is reconstructed.

Images