Terahertz synthetic aperture radar imaging system and method based on digital signal processing

Abstract

The present invention provides a terahertz synthetic aperture radar imaging system and method based on digital signal processing, belonging to the field of radar imaging technology, the imaging system comprises a preprocessing unit, a plurality of core processors and a shared storage unit; the preprocessing unit is used for performing motion compensation on original echo data; the shared storage unit is used for dividing the motion compensation data into a plurality of parts according to the number of core processors, and sending the parts to each core processor respectively; each core processor comprises a scale transformation module, an azimuth interpolation module and an imaging module; the scale transformation module obtains range transposed data according to motion compensation sub-data; the azimuth interpolation module obtains azimuth transposed data according to range transposed sub-data; the imaging module obtains final imaging data according to azimuth transposed sub-data. The imaging speed of the terahertz synthetic aperture radar is improved by processing data in parallel by a plurality of core processors and splicing and splitting the data of each processing stage in combination with the shared storage unit.

Description

Technical Field

[0001] The present invention relates to the field of radar imaging, and in particular to a terahertz synthetic aperture radar imaging system and method based on digital signal processing. Background Art

[0002] Synthetic Aperture Radar (SAR) is a microwave imaging radar that can achieve modern high-resolution imaging all day and all weather. Terahertz SAR has a narrow beam, good directivity, strong anti-interference ability, and the terahertz band is between millimeter waves and infrared light, with good penetration and high imaging resolution. It is widely used in national defense security, chemical identification, medical testing, bioengineering and other fields. Terahertz SAR imaging requires processing a large amount of data, and the echo needs to be motion compensated after dechirp reception, and the two-dimensional coupling in the frequency domain is serious.

[0003] In addition, terahertz SAR usually requires a large bandwidth to obtain a large amount of information. However, due to the immaturity of terahertz radiation sources, terahertz detectors, and many terahertz functional devices, as well as the low power and low energy conversion rate of terahertz radiation sources, terahertz SAR imaging takes a long time.

[0004] PFA (polar format algorithm) stores data in polar coordinate format. The conversion from polar coordinate format to matrix format is completed by two-dimensional resampling of distance and azimuth, and imaging is completed by two-dimensional fast Fourier transform. However, if the PFA algorithm is not combined with hardware to realize the imaging processing of huge terahertz data, the computational burden will be greatly increased, which is far from meeting the imaging speed requirements.

[0005] Based on the above problems, a new imaging method is urgently needed to improve the imaging speed. Summary of the invention

[0006] The purpose of the present invention is to provide a terahertz synthetic aperture radar imaging system and method based on digital signal processing, which can improve the imaging speed.

[0007] To achieve the above object, the present invention provides the following solutions:

[0008] A terahertz synthetic aperture radar imaging system based on digital signal processing, the terahertz synthetic aperture radar imaging system based on digital signal processing comprising a preprocessing unit, a plurality of core processors and a shared storage unit;

[0009] The preprocessing unit is connected to the synthetic aperture radar, and is used to perform motion compensation on the original echo data to obtain motion compensated data; the original echo data is the echo data collected by the synthetic aperture radar;

[0010] The shared storage unit is connected to the preprocessing unit and each core processor respectively; the shared storage unit is used to divide the motion compensation data into a plurality of motion compensation sub-data according to the number of core processors, and send each motion compensation sub-data to each core processor respectively;

[0011] Each core processor includes: a scale transformation module, an azimuth interpolation module and an imaging module;

[0012] The scale transformation module is connected to the shared storage unit, and is used to perform distance scale transformation on the motion compensation sub-data, and perform transposition processing to obtain distance transposed data, and send the data to the shared storage unit;

[0013] The shared storage unit is also used to perform splicing processing on the distance transposed data sent by each scale transformation module to obtain the distance transposed total data, and divide the distance transposed total data into a plurality of distance transposed sub-data, and send each distance transposed sub-data to each orientation interpolation module respectively;

[0014] The azimuth interpolation module is connected to the shared storage unit, and is used for performing azimuth interpolation processing on the range transposed sub-data using a Sinc interpolation method to obtain azimuth processed data, performing azimuth fast Fourier transform on the azimuth processed data, and performing transposition processing to obtain azimuth transposed data, and sending the data to the shared storage unit;

[0015] The shared storage unit is also used to perform splicing processing on the azimuth transposed data sent by each azimuth interpolation module to obtain azimuth transposed total data, and divide the azimuth transposed total data into a plurality of azimuth transposed sub-data, and send each azimuth transposed sub-data to each imaging module respectively;

[0016] The imaging module is connected to the shared storage unit, and is used to perform one-dimensional Fourier transform on the azimuth transposed sub-data to obtain initial imaging data, perform one-dimensional self-focusing processing on the initial imaging data to obtain final imaging data, and send the final imaging data to the shared storage unit;

[0017] The shared storage unit is also used to perform splicing processing on the imaging data sent by each imaging module to obtain final image information.

[0018] Optionally, the scale transformation module includes:

[0019] A scaling processing sub-module, connected to the shared storage unit, for multiplying the motion compensation sub-data by a scaling function to obtain scaled function processed data; the scaling function is: where, Φ scl (τ) is the scaling function, τ is the range instantaneous fast time, j is the complex number, kr is the frequency modulation slope, δ r is the range scaling factor, Φ ref is the elevation angle of the synthetic aperture radar to the center point of the scene, θ is the instantaneous azimuth angle during the flight of the synthetic aperture radar, Φ is the instantaneous elevation angle of the synthetic aperture radar to the center of the scene during flight, t r is the range fast time, R a is the instantaneous distance between the synthetic aperture radar and the scene, c is the speed of light;

[0020] A filtering sub-module, connected to the scaling processing sub-module, for performing a one-dimensional fast Fourier transform on the scaled function processed data and multiplying it by a discrete filter to obtain filtered processed data; the discrete filter is where, H2 is the discrete filter, fr is the range frequency, fc is the carrier frequency, and R0 is the operating range of the synthetic aperture radar;

[0021] A transpose sub-module, connected to the filtering sub-module and the shared storage unit respectively, for transforming the filtered processed data into the time domain using the one-dimensional inverse fast Fourier transform method and multiplying it by the inverse scaling function to obtain range transposed data, and sending it to the shared storage unit; the inverse scaling function is: where, Φ ins is the inverse scaling function.

[0022] Optionally, the imaging module includes:

[0023] An initial imaging sub-module, connected to the shared storage unit, for performing a one-dimensional Fourier transform on the azimuth transposed sub-data to obtain initial imaging data;

[0024] A shifting sub-module, connected to the initial imaging sub-module, for the m-th iteration, determining the strongest scatter point in each range cell of the initial imaging data according to the maximum energy criterion, and performing circular shifting on each strongest scatter point so that the strongest scatter point in each range cell moves to the center of the image to obtain center shifted data; 0 < m < M, M is the iteration threshold;

[0025] A windowing sub-module, connected to the shifting sub-module, for performing windowing processing on the center shifted data to obtain Hamming window processed data;

[0026] An inverse transform submodule, connected to the windowing submodule, for performing an inverse fast Fourier transform on the Hamming window processed data to obtain inverse fast Fourier transform data;

[0027] A conjugate submodule, connected to the inverse transform submodule, for conjugating the inverse fast Fourier transform data to obtain conjugate processed data;

[0028] an error gradient determination submodule, connected to the inverse transformation submodule and the conjugate submodule respectively, and used for multiplying the inverse fast Fourier transform data with the conjugate processed data to obtain a phase error gradient;

[0029] an error determination submodule, connected to the error gradient determination submodule, and configured to determine a phase error value according to the phase error gradient;

[0030] a compensation submodule, connected to the shift submodule and the error determination submodule respectively, for multiplying the phase error value by each strongest scattering point respectively to obtain corresponding compensation processing data, and replacing the corresponding strongest scattering point with the compensation processing data;

[0031] A first transformation submodule is connected to the initial imaging submodule and the compensation submodule respectively, and is used to transform each initial imaging data into a phase history domain by using a one-dimensional inverse fast Fourier transform method when the number of iterations is greater than or equal to an iteration threshold, and multiply the initial imaging data by the current phase error value to obtain error imaging data;

[0032] The second transformation submodule is connected to the first transformation submodule and the shared storage unit respectively, and is used to transform the error imaging data into the image domain by adopting the fast Fourier transform method, obtain the final imaging data, and send it to the shared storage unit.

[0033] To achieve the above object, the present invention also provides the following solution:

[0034] A terahertz synthetic aperture radar imaging method based on digital signal processing, comprising:

[0035] Performing motion compensation on original echo data to obtain motion compensated data; the original echo data is echo data collected by synthetic aperture radar;

[0036] According to the number of core processors, the motion compensation processing data is divided into a plurality of motion compensation sub-data;

[0037] For each piece of motion compensation sub-data, the core processor performs distance-direction scale transformation and transposition processing on the motion compensation sub-data to obtain distance-direction transposed data;

[0038] Performing splicing processing on each distance transposed data to obtain distance transposed total data, and dividing the distance transposed total data into a plurality of distance transposed sub-data;

[0039] For each range transposed sub-data, the core processor uses the Sinc interpolation method to perform azimuth interpolation processing on the range transposed sub-data to obtain azimuth processed data; the azimuth processed data is subjected to azimuth fast Fourier transform and transposed processing to obtain azimuth transposed data;

[0040] Performing splicing processing on each azimuth transposed data to obtain azimuth transposed total data, and dividing the azimuth transposed total data into a plurality of azimuth transposed sub-data;

[0041] For each azimuth transposed sub-data, a core processor performs a one-dimensional Fourier transform on the azimuth transposed sub-data to obtain initial imaging data; and performs a one-dimensional self-focusing process on the initial imaging data to obtain final imaging data;

[0042] The imaging data are stitched together to obtain the final image information.

[0043] Optionally, performing motion compensation on the echo data to obtain motion compensated data specifically includes:

[0044] The original echo data is multiplied by the motion compensation phase to obtain motion compensation processing data; the motion compensation phase is: Among them, Φ rmc (n,t r ) is the motion compensation phase, n is the number of pulses, kr is the frequency modulation slope, B r is the bandwidth, T P is the pulse width, c is the speed of light, fc is the carrier frequency, t r is the distance to fast time, R q =-(R a -R0), R a is the instantaneous distance between the synthetic aperture radar and the scene, and R0 is the effective range of the synthetic aperture radar.

[0045] Optionally, performing distance-direction scale transformation and transposition processing on the motion compensation sub-data to obtain distance-direction transposed data specifically includes:

[0046] The motion compensation sub-data is multiplied by a scale function to obtain scale function processed data; the scale function is: Among them, Φ scl (τ) is the scale function, τ is the instantaneous fast time in the range direction, j is a complex number, kr is the frequency modulation slope, δ r is the distance-to-scale transformation factor, Φ ref The pitch angle of the synthetic aperture radar to the center point of the scene, θ is the instantaneous azimuth angle during the flight of the synthetic aperture radar, Φ is the instantaneous pitch angle from the synthetic aperture radar to the center of the scene during flight, and t r is the fast time in the range direction, and R a is the instantaneous distance between the synthetic aperture radar and the scene, and c is the speed of light;

[0047] Perform one-dimensional fast Fourier transform on the data processed by the scale function and multiply it by the discrete filter to obtain the filtered data; the discrete filter is where H2 is the discrete filter, fr is the range direction frequency, fc is the carrier frequency, and R0 is the operating range of the synthetic aperture radar;

[0048] Use the one-dimensional inverse fast Fourier transform method to transform the filtered data into the time domain and multiply it by the inverse scale function to obtain the range direction transposed data; the inverse scale function is:

[0049] where, Φ ins is the inverse scale function.

[0050] Optionally, perform azimuth interpolation processing on the range direction transposed data according to the following formula to obtain azimuth processed data:

[0051]

[0052] where g(x) is the azimuth processed data, g d (i) is the range direction transposed data, and sinc is the convolution kernel.

[0053] Optionally, the one-dimensional autofocus processing of the initial imaging data to obtain the final imaging data specifically includes:

[0054] For the m-th iteration, according to the energy maximum criterion, determine the strongest scattering point in each range cell of the initial imaging data, and perform circular shift on each strongest scattering point to move the strongest scattering point in each range cell to the center of the image to obtain the center-shifted data; 0 < m < M, where M is the iteration threshold;

[0055] Perform windowing processing on the center-shifted data to obtain the Hamming window processed data;

[0056] Perform inverse fast Fourier transform on the Hamming window processed data to obtain the inverse fast Fourier transform data;

[0057] Take the conjugate of the inverse fast Fourier transform data to obtain the conjugate processed data;

[0058] Multiplying the inverse fast Fourier transform data with the conjugate processed data to obtain a phase error gradient;

[0059] Determining a phase error value according to the phase error gradient;

[0060] Multiplying the phase error value by each strongest scattering point to obtain corresponding compensation processing data;

[0061] Use compensation processing data to replace the corresponding strongest scattering point;

[0062] When the number of iterations is greater than or equal to the iteration threshold, the one-dimensional inverse fast Fourier transform method is used to transform each initial imaging data into the phase history domain, and multiply it with the current phase error value to obtain error imaging data;

[0063] The error imaging data is transformed into an image domain by using a fast Fourier transform method to obtain final imaging data.

[0064] Optionally, performing windowing processing on the center shifted data to obtain Hamming window processed data specifically includes:

[0065] The center shifted data is multiplied by a Hamming window function to obtain Hamming window processed data; the Hamming window function is: Where w(n) is the Han window function and w is the window width.

[0066] According to the specific embodiment provided by the present invention, the present invention discloses the following technical effects: by processing data in parallel through multiple core processors, and combining with a shared storage unit to splice and split the data of each processing stage, the difficulty and long processing time of terahertz SAR imaging are eliminated, and the imaging speed is improved. In addition, the range is transposed through the scale transformation module, and the azimuth is transposed through the azimuth interpolation module, which further improves the efficiency of imaging. BRIEF DESCRIPTION OF THE DRAWINGS

[0067] In order to more clearly illustrate the embodiments of the present invention or the technical solutions in the prior art, the drawings required for use in the embodiments will be briefly introduced below. Obviously, the drawings described below are only some embodiments of the present invention. For ordinary technicians in this field, other drawings can be obtained based on these drawings without paying creative labor.

[0068] Figure 1 It is a schematic diagram of the module structure of the terahertz synthetic aperture radar imaging system of the present invention;

[0069] Figure 2 This is a flow chart of the terahertz synthetic aperture radar imaging method of the present invention;

[0070] Figure 3Flowchart for determining a method for range-transposed data.

[0071] Explanation of symbols:

[0072] Synthetic aperture radar-1, preprocessing unit-2, shared storage unit-3, core processor-4, scale transformation module-41, azimuth interpolation module-42, imaging module-43, secondary storage module-44. DETAILED DESCRIPTION

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

[0074] The purpose of the present invention is to provide a terahertz synthetic aperture radar imaging system and method, which processes data in parallel through multiple core processors and combines a shared storage unit to splice and split the data of each processing stage, thereby getting rid of the difficulty and long time of terahertz SAR imaging processing and improving the imaging speed.

[0075] In order to make the above-mentioned objects, features and advantages of the present invention more obvious and easy to understand, the present invention is further described in detail below with reference to the accompanying drawings and specific embodiments.

[0076] like Figure 1 As shown, the terahertz synthetic aperture radar 1 imaging system based on digital signal processing of the present invention comprises a preprocessing unit 2, a plurality of core processors 4 and a shared storage unit 3. Preferably, the number of the core processors 4 is 8.

[0077] The preprocessing unit 2 is connected to the synthetic aperture radar 1 and is used to perform motion compensation on the original echo data to obtain motion compensated data. The original echo data is the echo data collected by the synthetic aperture radar 1.

[0078] The shared storage unit 3 is connected to the preprocessing unit 2 and each core processor 4 respectively; the shared storage unit 3 is used to divide the motion compensation data into multiple motion compensation sub-data according to the number of core processors 4, and send each motion compensation sub-data to each core processor 4 respectively.

[0079] Each core processor 4 includes: a scale transformation module 41 , an orientation interpolation module 42 and an imaging module 43 .

[0080] The scale transformation module 41 is connected to the shared storage unit 3, and is used to perform distance scale transformation (principle of chirp scaling, PCS) on the motion compensation sub-data, and perform transposition processing to obtain distance transposed data, and send it to the shared storage unit 3. In this embodiment, the scale transformation module 41 performs distance PCS processing to complete distance resampling.

[0081] The shared storage unit 3 is also used to splice the distance transposed data sent by each scale transformation module 41 to obtain the total distance transposed data, and divide the total distance transposed data into multiple distance transposed sub-data, and send each distance transposed sub-data to each orientation interpolation module 42 respectively.

[0082] In addition, the core processor 4 further includes a secondary storage module 44, which is respectively connected to the scale transformation module 41, the azimuth interpolation module 42 and the imaging module 43. The secondary storage module 44 is used to store the range transposed data, the azimuth transposed data and the imaging data.

[0083] Specifically, the motion compensation data is divided into 8 parts according to the number of core processors, and evenly distributed to 8 core processors 4, and each core processor 4 processes Na / 8 pulses. Each core processor transfers the motion compensation sub-data of the processed pulse number to the secondary storage module 44 of each core processor through EDMA. After the processing is completed, each core processor sends the result to the shared storage unit 3, and the cycle is repeated Nr times to complete multi-core synchronization. Finally, the data processed by PCS is transposed to obtain the range transposed data. Starting from the first range gate, Nr range transposed data are transferred to the secondary storage module 44 through EDMA and then sent to the shared storage unit 3 through EDMA. The cycle is repeated Na times to obtain the total range transposed data. Wherein Nr is the number of points in the range, and Na is the number of pulses in the azimuth.

[0084] The azimuth interpolation module 42 is connected to the shared storage unit 3, and is used to perform azimuth interpolation processing on the range transposed sub-data using the Sinc interpolation method to obtain azimuth processed data, perform azimuth fast Fourier transform on the azimuth processed data, perform transposition processing, obtain azimuth transposed data, and send it to the shared storage unit 3. In this embodiment, azimuth resampling is completed by the azimuth interpolation module 42.

[0085] Specifically, according to the following formula, the range transposed sub-data is subjected to azimuth interpolation processing to obtain azimuth processed data:

[0086]

[0087] Among them, g(x) is the azimuth processing data, g d (i) is the distance transposed sub-data, and sinc is the convolution kernel.

[0088] In this embodiment, the total distance transposed data is first divided into 8 equal parts according to the number of parallel processing cores (the number of core processors), and each core processes Nr / 8 distance gate data respectively. During the Sinc interpolation process, each core transfers the distance transposed sub-data to the secondary storage module 44 in each core through EDMA (Enhanced Direct Memory Access). In order to realize hardware acceleration processing, a table lookup method can also be used instead of the Sinc interpolation method to perform azimuth interpolation processing on the distance transposed sub-data. The table used in the table lookup method is a pre-set table. Specifically, the processing process of the table lookup method includes: generating a Sinc coefficient table by calculation in advance, and then determining the required value by address lookup. No need for multiple calculations to save time.

[0089] After the interpolation processing is completed, each core processor performs an azimuth FFT (fast Fourier transform) operation on the azimuth processing data to obtain azimuth FFT processing data. After the azimuth FFT operation is completed, the result is moved from the secondary storage module 44 of the core processor 4 to the shared storage unit 3, and the cycle is repeated Na times to complete multi-core synchronization. Finally, the azimuth FFT processing data is transposed to obtain azimuth transposed data. Starting from the first pulse, Na data is moved to the secondary storage module 44 through EDMA and then moved back to the shared storage unit 3 through EDMA. The cycle is repeated Nr times to obtain the total azimuth transposed data.

[0090] The shared storage unit 3 is also used to splice the azimuth transposed data sent by each azimuth interpolation module 42 to obtain the azimuth transposed total data, and divide the azimuth transposed total data into multiple azimuth transposed sub-data, and send each azimuth transposed sub-data to each imaging module 43 respectively.

[0091] The imaging module 43 is connected to the shared storage unit 3, and is used to perform one-dimensional Fourier transform on the azimuth transposed sub-data to obtain initial imaging data, perform one-dimensional autofocus processing on the initial imaging data to obtain final imaging data, and send the final imaging data to the shared storage unit 3. Specifically, the initial imaging data is PFA imaging data.

[0092] Specifically, the total azimuthally transposed data is evenly distributed to each core processor, and each core processor processes Na / 8 pulses respectively. Each core processor transfers the azimuthally transposed sub-data in the shared storage unit 3 to the secondary storage module 44 within the core through EDMA. After completing the FFT operation, the FFT-processed data is sent to the shared storage unit 3. The cycle is repeated Nr times to complete multi-core synchronization.

[0093] The shared storage unit 3 is also used to perform splicing processing on the imaging data sent by each imaging module 43 to obtain final image information.

[0094] Furthermore, the scale transformation module 41 includes: a scale processing submodule, a filtering submodule and a transposition submodule.

[0095] The scale processing submodule is connected to the shared storage unit 3, and is used to multiply the motion compensation subdata with a scale function to obtain scale function processed data. The scale function is: Among them, Φ scl (τ) is the scale function, τ is the instantaneous fast time in the range direction, j is a complex number, kr is the frequency modulation slope, δ r is the distance-to-scale transformation factor, Φ ref is the pitch angle of the synthetic aperture radar to the center of the scene, θ is the instantaneous azimuth angle of the synthetic aperture radar during flight, Φ is the instantaneous pitch angle of the synthetic aperture radar from the center of the scene during flight, t r is the distance to fast time (vector), R a is the instantaneous distance between the synthetic aperture radar and the scene, and c is the speed of light.

[0096] The filtering submodule is connected to the scale processing submodule, and the filtering submodule is used to perform a one-dimensional fast Fourier transform on the scale function processing data and multiply it with the discrete filter to obtain the filtered processing data. The discrete filter is Among them, H2 is a discrete filter, fr is the range frequency, fc is the carrier frequency, and R0 is the effective range of the synthetic aperture radar.

[0097] The transposition submodule is connected to the filtering submodule and the shared storage unit 3 respectively. The transposition submodule is used to transform the filtered data into the time domain by using a one-dimensional inverse fast Fourier transform method (IFFT), and multiply it with an inverse scaling function to obtain the range transposed data, and send it to the shared storage unit 3; the inverse scaling function is: Among them, Φ ins is the inverse scaling function.

[0098] Further, the imaging module 43 includes: an initial imaging sub-module, a shifting sub-module, a windowing sub-module, an inverse transform sub-module, a conjugate sub-module, an error gradient determination sub-module, an error determination sub-module, a compensation sub-module, a first transform sub-module, and a second transform sub-module.

[0099] Among them, the initial imaging sub-module is connected to the shared storage unit 3, and the initial imaging sub-module is used to perform a one-dimensional Fourier transform on the azimuth transposed data to obtain initial imaging data.

[0100] The shifting sub-module is connected to the initial imaging sub-module. For the nth iteration, the shifting sub-module is used to determine the strongest scattering point in each range cell of the initial imaging data according to the maximum energy criterion, and perform a cyclic shift on each strongest scattering point so that the strongest scattering point in each range cell moves to the center of the image to obtain center-shifted data; 0 < n < N, where N is the iteration threshold. In this embodiment, the iteration threshold is 6.

[0101] Specifically, according to the maximum energy criterion, the strongest scattering point in each range cell is found, generally about 5% of the total number of range gates. The selected data is transferred to the secondary storage module 44, and then the selected data is cyclically shifted so that the strongest scattering point in each range cell moves to the center of the image to obtain center-shifted processed data.

[0102] The windowing sub-module is connected to the shifting sub-module, and the windowing sub-module is used to perform windowing processing on the center-shifted data to obtain Hamming window processed data. By performing windowing processing on the center-shifted data, the signal-to-noise ratio is improved, and the influence caused by interfering targets and noise is reduced.

[0103] The construction method of the Hamming window is as follows: where w is the window width.

[0104] The inverse transform sub-module is connected to the windowing sub-module, and the inverse transform sub-module is used to perform an inverse fast Fourier transform on the Hamming window processed data to obtain inverse fast Fourier transform data.

[0105] The conjugate sub-module is connected to the inverse transform sub-module, and the conjugate sub-module is used to find the conjugate of the inverse fast Fourier transform data to obtain conjugate processed data.

[0106] The error gradient determination sub-module is respectively connected to the inverse transform sub-module and the conjugate sub-module, and the error gradient determination sub-module is used to multiply the inverse fast Fourier transform data by the conjugate processed data to obtain a phase error gradient.

[0107] The error determination submodule is connected to the error gradient determination submodule, and the error determination submodule is used to determine the phase error value according to the phase error gradient. Specifically, the phase error gradient result is matrix-transposed, summed according to the range gate, and the inverse tangent of the sum is calculated to obtain the desired phase error value.

[0108] The compensation submodule is connected to the shift submodule and the error determination submodule respectively, and is used to multiply the phase error value by each strongest scattering point respectively to obtain corresponding compensation processing data, and use the compensation processing data to replace the corresponding strongest scattering point.

[0109] The first transformation submodule is connected to the initial imaging submodule and the compensation submodule respectively. The first transformation submodule is used to transform each initial imaging data into a phase history domain by using a one-dimensional inverse fast Fourier transform method when the number of iterations is greater than or equal to an iteration threshold, and multiply it with the current phase error value to obtain error imaging data.

[0110] The second transformation submodule is connected to the first transformation submodule and the shared storage unit 3 respectively, and the second transformation submodule is used to transform the error imaging data into the image domain by adopting the fast Fourier transform method, obtain the final imaging data, and send it to the shared storage unit 3.

[0111] During the processing, the total azimuth transposed data is first divided equally among the core processors. Each core processor processes Nr / 8 range gate data. According to the maximum energy criterion, about 5% of the range units in the total azimuth transposed data are selected and stored in the shared storage unit to complete multi-core synchronization. Then the selected data is divided into 8 parts according to the number of cores and distributed to each core processor. Each core processor moves the selected data to the secondary storage module in the core through EDMA, performs center shift, windowing, phase error gradient estimation and 6 rounds of correction to obtain phase error data. Finally, the PFA imaging data is multiplied by the phase error after the azimuth IFFT operation. The result is then subjected to the azimuth FFT operation to obtain the final imaging data.

[0112] The preprocessing unit is mainly used to perform motion compensation on the original echo data. The instantaneous distance R between the synthetic aperture radar and the scene a There is a deviation between the effective range R0 of the synthetic aperture radar. In fact, R0 is the R used in real-time and real-time compensation. a Quantitative estimation of the SAR, usually R0 is used as the demodulation signal, and the difference between the two is usually expressed as the motion compensation phase Φ rmc (n,t r ) represents the azimuth frequency of the scene center. aThe changes in R0 over the azimuth aperture range can seriously affect the azimuth signal. The residual range frequency and azimuth phase of the signal received from the scatterer in the center of the scene can be eliminated through fine motion compensation. After motion compensation, the residual range frequency and azimuth phase are eliminated, and the imaging details are better.

[0113] like Figure 2 As shown, the terahertz synthetic aperture radar imaging method based on digital signal processing of the present invention includes:

[0114] S1: performing motion compensation on original echo data to obtain motion compensated data; the original echo data is echo data collected by synthetic aperture radar.

[0115] S2: Divide the motion compensation data into a plurality of motion compensation sub-data according to the number of core processors.

[0116] S3: for each piece of motion compensation sub-data, the core processor performs distance-direction scale transformation and transposition processing on the motion compensation sub-data to obtain distance-direction transposed data.

[0117] Terahertz data processing in the beamforming mode is stored in polar coordinate format, and direct two-dimensional fast Fourier transform imaging will reduce the image focusing quality. Therefore, the present invention first uses distance PCS (principle of chirpscaling) processing to complete distance resampling, thereby improving the quality of image focusing.

[0118] S4: performing concatenation processing on each distance transposed data to obtain distance transposed total data, and dividing the distance transposed total data into a plurality of distance transposed sub-data.

[0119] S5: For each range transposed sub-data, the core processor uses the Sinc interpolation method to perform azimuth interpolation processing on the range transposed sub-data to obtain azimuth processed data; the azimuth processed data is subjected to azimuth fast Fourier transform and transposed processing to obtain azimuth transposed data.

[0120] S6: performing splicing processing on each azimuthally transposed data to obtain azimuthally transposed total data, and dividing the azimuthally transposed total data into a plurality of azimuthally transposed sub-data.

[0121] S7: for each azimuthally transposed sub-data, perform one-dimensional Fourier transform on the azimuthally transposed sub-data through a core processor to obtain initial imaging data; perform one-dimensional self-focusing processing on the initial imaging data to obtain final imaging data.

[0122] S8: performing stitching processing on each imaging data to obtain final image information.

[0123] Furthermore, step S1 specifically includes: multiplying the original echo data with the motion compensation phase to obtain motion compensation processing data; the motion compensation phase is: Among them, Φ rmc (n,t r ) is the motion compensation phase, n is the number of pulses, kr is the frequency modulation slope, B r is the bandwidth, T P is the pulse width, c is the speed of light, fc is the carrier frequency, t r is the distance to fast time, R q =-(R a -R0), R a is the instantaneous distance between the synthetic aperture radar and the scene, and R0 is the effective range of the synthetic aperture radar.

[0124] Furthermore, if Figure 3 As shown, step S3 specifically includes:

[0125] S31: Multiply the motion compensation sub-data S(τ,t) by a scaling function to obtain scaling function processed data. The scaling function is: Among them, Φ scl (τ) is the scale function, τ is the instantaneous fast time in the range direction, j is a complex number, kr is the frequency modulation slope, δ r is the distance-to-scale transformation factor, Φ ref is the pitch angle of the synthetic aperture radar to the center of the scene, θ is the instantaneous azimuth angle of the synthetic aperture radar during flight, Φ is the instantaneous pitch angle of the synthetic aperture radar from the center of the scene during flight, t r is the distance to fast time, R a is the instantaneous distance between the synthetic aperture radar and the scene, and c is the speed of light.

[0126] S32: Perform a one-dimensional fast Fourier transform on the scale function processed data and multiply it with a discrete filter to obtain filtered processed data. The discrete filter is Among them, H2 is a discrete filter, fr is the range frequency, fc is the carrier frequency, and R0 is the effective range of the synthetic aperture radar.

[0127] S33: transform the filtered data into the time domain using a one-dimensional inverse fast Fourier transform method, and multiply it with an inverse scaling function to obtain range transposed data. The inverse scaling function is:

[0128] Among them, Φ ins is the inverse scaling function.

[0129] Specifically, according to the following formula, azimuth interpolation processing is performed on the range transposed sub-data to obtain azimuth processed data:

[0130]

[0131] where g(x) is the azimuth processed data, g d (i) is the range transposed sub-data, and sinc is the convolution kernel.

[0132] In step S7, one-dimensional autofocus processing is performed on the initial imaging data to obtain the final imaging data, which specifically includes:

[0133] S701: For the m-th iteration, according to the energy maximum criterion, determine the strongest scattering point in each range cell of the initial imaging data, and perform circular shifting on each strongest scattering point so that the strongest scattering point in each range cell moves to the image center to obtain center-shifted data; 0 < m < M, where M is the iteration threshold.

[0134] S702: Perform windowing processing on the center-shifted data to obtain Hamming window processed data. Specifically, multiply the center-shifted data by the Hamming window function to obtain the Hamming window processed data; the Hamming window function is: where w(n) is the Hamming window function and w is the window width.

[0135] S703: Perform inverse fast Fourier transform on the Hamming window processed data to obtain inverse fast Fourier transform data.

[0136] S704: Take the conjugate of the inverse fast Fourier transform data to obtain conjugate processed data.

[0137] S705: Multiply the inverse fast Fourier transform data by the conjugate processed data to obtain the phase error gradient.

[0138] S706: Determine the phase error value according to the phase error gradient.

[0139] S707: Multiply the phase error value by each strongest scattering point respectively to obtain the corresponding compensation processed data.

[0140] S708: Replace the corresponding strongest scattering point with the compensation processed data.

[0141] S709: When the number of iterations is greater than or equal to the iteration threshold, use the one-dimensional inverse fast Fourier transform method to transform each initial imaging data into the phase history domain and multiply it by the current phase error value to obtain error imaging data.

[0142] S710: using a fast Fourier transform method to transform the error imaging data into an image domain to obtain final imaging data.

[0143] After resampling in range and azimuth, the imaging results are improved, and the image focusing effect is good. Therefore, the present invention is suitable for terahertz SAR imaging. By utilizing the high-speed signal processing capability of multi-core DSP to divide the huge amount of terahertz data into memory, data transfer and multi-core parallelism to achieve hardware acceleration processing, it has practical significance such as fast imaging speed, good imaging quality and engineering research. Table 1 shows the processing time of the terahertz imaging process based on multi-core DSP, from which it can be seen that the acceleration effect is obvious.

[0144] Table 1 Statistics of the measured processing time of each module

[0145]

[0146] Compared with the prior art, the terahertz synthetic aperture radar imaging method of the present invention has the same beneficial effects as the terahertz synthetic aperture radar imaging system, which will not be described in detail here.

[0147] In this specification, each embodiment is described in a progressive manner, and each embodiment focuses on the differences from other embodiments. The same or similar parts between the embodiments can be referred to each other. For the system disclosed in the embodiment, since it corresponds to the method disclosed in the embodiment, the description is relatively simple, and the relevant parts can be referred to the method part.

[0148] The principles and implementation methods of the present invention are described in this article using specific examples. The description of the above embodiments is only used to help understand the method and core idea of ​​the present invention. At the same time, for those skilled in the art, according to the idea of ​​the present invention, there will be changes in the specific implementation methods and application scope. In summary, the content of this specification should not be understood as limiting the present invention.

Claims

1. A terahertz synthetic aperture radar imaging system based on digital signal processing, characterized in that: The terahertz synthetic aperture radar imaging system based on digital signal processing includes a preprocessing unit, multiple core processors, and a shared storage unit; The preprocessing unit is connected to the synthetic aperture radar. The preprocessing unit is used to perform motion compensation on the original echo data to obtain motion-compensated data; the original echo data is the echo data collected by the synthetic aperture radar; The shared storage unit is respectively connected to the preprocessing unit and each core processor; the shared storage unit is used to divide the motion-compensated data into multiple pieces of motion-compensated sub-data according to the number of core processors, and send each piece of motion-compensated sub-data to each core processor respectively; Each core processor includes: a scaling transformation module, an azimuth interpolation module, and an imaging module; The scaling transformation module is connected to the shared storage unit. The scaling transformation module is used to perform range-direction scaling transformation on the motion-compensated sub-data, and perform a transpose process to obtain range-direction transposed data, and send it to the shared storage unit; The shared storage unit is further used to perform splicing processing on the range-direction transposed data sent by each scaling transformation module to obtain range-direction transposed total data, and divide the range-direction transposed total data into multiple pieces of range-direction transposed sub-data, and send each piece of range-direction transposed sub-data to each azimuth interpolation module respectively; The azimuth interpolation module is connected to the shared storage unit. The azimuth interpolation module is used to perform azimuth-direction interpolation processing on the range-direction transposed sub-data by using the Sinc interpolation method to obtain azimuth-direction processed data, perform azimuth-direction fast Fourier transform on the azimuth-direction processed data, and perform a transpose process to obtain azimuth-direction transposed data, and send it to the shared storage unit; The shared storage unit is further used to perform splicing processing on the azimuth-direction transposed data sent by each azimuth interpolation module to obtain azimuth-direction transposed total data, and divide the azimuth-direction transposed total data into multiple pieces of azimuth-direction transposed sub-data, and send each piece of azimuth-direction transposed sub-data to each imaging module respectively; The imaging module is connected to the shared storage unit. The imaging module is used to perform one-dimensional Fourier transform on the azimuth-direction transposed sub-data to obtain initial imaging data, perform one-dimensional autofocus processing on the initial imaging data to obtain the final imaging data, and send it to the shared storage unit; The imaging module includes: An initial imaging sub-module, connected to the shared storage unit, for performing one-dimensional Fourier transform on the azimuth-direction transposed sub-data to obtain initial imaging data; A shifting sub-module, connected to the initial imaging sub-module, for the m-th iteration, according to the energy maximum criterion, determining the strongest scattering point in each range cell of the initial imaging data, and performing cyclic shifting on each strongest scattering point to move the strongest scattering point in each range cell to the image center to obtain center-shifted data; 0 < m < M, M is the iteration threshold; A windowing sub-module, connected to the shifting sub-module, for performing windowing processing on the center-shifted data to obtain Hamming window processed data; An inverse transform submodule, connected to the windowing submodule, for performing an inverse fast Fourier transform on the Hamming window processed data to obtain inverse fast Fourier transform data; A conjugate submodule, connected to the inverse transform submodule, for conjugating the inverse fast Fourier transform data to obtain conjugate processed data; an error gradient determination submodule, connected to the inverse transformation submodule and the conjugate submodule respectively, and used for multiplying the inverse fast Fourier transform data with the conjugate processed data to obtain a phase error gradient; an error determination submodule, connected to the error gradient determination submodule, and configured to determine a phase error value according to the phase error gradient; a compensation submodule, connected to the shift submodule and the error determination submodule respectively, for multiplying the phase error value by each strongest scattering point respectively to obtain corresponding compensation processing data, and replacing the corresponding strongest scattering point with the compensation processing data; A first transformation submodule is connected to the initial imaging submodule and the compensation submodule respectively, and is used to transform each initial imaging data into a phase history domain by using a one-dimensional inverse fast Fourier transform method when the number of iterations is greater than or equal to an iteration threshold, and multiply the initial imaging data by the current phase error value to obtain error imaging data; a second transformation submodule, connected to the first transformation submodule and the shared storage unit respectively, and configured to transform the error imaging data into an image domain by adopting a fast Fourier transform method, obtain final imaging data, and send the final imaging data to the shared storage unit; The shared storage unit is also used to perform splicing processing on the imaging data sent by each imaging module to obtain final image information.

2. The terahertz synthetic aperture radar imaging system based on digital signal processing according to claim 1, characterized in that: The scale transformation module comprises: The scale processing submodule is connected to the shared storage unit and is used to multiply the motion compensation subdata with a scale function to obtain scale function processed data; the scale function is: Among them, Φ scl (τ) is the scale function, τ is the instantaneous fast time in the range direction, j is a complex number, kr is the frequency modulation slope, δ r is the distance-to-scale transformation factor, Φ ref is the pitch angle of the synthetic aperture radar to the center of the scene, θ is the instantaneous azimuth angle of the synthetic aperture radar during flight, Φ is the instantaneous pitch angle of the synthetic aperture radar from the center of the scene during flight, t r is the distance to fast time, R a is the instantaneous distance between the synthetic aperture radar and the scene, c is the speed of light; The filtering submodule is connected to the scale processing submodule and is used to perform a one-dimensional fast Fourier transform on the scale function processing data and multiply it with a discrete filter to obtain filtered processing data; the discrete filter is Among them, H2 is a discrete filter, fr is the range frequency, fc is the carrier frequency, and R0 is the effective range of the synthetic aperture radar; The transposition submodule is connected to the filtering submodule and the shared storage unit respectively, and is used to transform the filtered data into the time domain by using a one-dimensional inverse fast Fourier transform method, and multiply it with an inverse scaling function to obtain the range transposed data, and send it to the shared storage unit; the inverse scaling function is: Among them, Φ ins is the inverse scaling function.

3. A terahertz synthetic aperture radar imaging method based on digital signal processing, applied to the terahertz synthetic aperture radar imaging system based on digital signal processing according to any one of claims 1-2, characterized in that: The terahertz synthetic aperture radar imaging method based on digital signal processing comprises: Performing motion compensation on original echo data to obtain motion compensated data; the original echo data is echo data collected by synthetic aperture radar; According to the number of core processors, the motion compensation data is divided into a plurality of motion compensation sub-data; For each piece of motion compensation sub-data, the core processor performs distance-direction scale transformation and transposition processing on the motion compensation sub-data to obtain distance-direction transposed data; Performing splicing processing on each distance transposed data to obtain distance transposed total data, and dividing the distance transposed total data into a plurality of distance transposed sub-data; For each range transposed sub-data, the core processor uses the Sinc interpolation method to perform azimuth interpolation processing on the range transposed sub-data to obtain azimuth processed data; the azimuth processed data is subjected to azimuth fast Fourier transform and transposed processing to obtain azimuth transposed data; Performing splicing processing on each azimuth transposed data to obtain azimuth transposed total data, and dividing the azimuth transposed total data into a plurality of azimuth transposed sub-data; For each azimuth transposed sub-data, a core processor performs a one-dimensional Fourier transform on the azimuth transposed sub-data to obtain initial imaging data; and performs a one-dimensional self-focusing process on the initial imaging data to obtain final imaging data; Stitch the imaging data to obtain the final image information.

4. The terahertz synthetic aperture radar imaging method based on digital signal processing according to claim 3, characterized in that: Motion compensation is performed on the original echo data to obtain motion-compensated data, which specifically includes: The original echo data is multiplied by the motion compensation phase to obtain the motion compensation data; the motion compensation phase is: Among them, Φ rmc (n,t r ) is the motion compensation phase, n is the number of pulses, kr is the frequency modulation slope, B r is the bandwidth, T P is the pulse width, c is the speed of light, fc is the carrier frequency, t r is the distance to fast time, R q =-(R a -R0), R a is the instantaneous distance between the synthetic aperture radar and the scene, and R0 is the effective range of the synthetic aperture radar.

5. The terahertz synthetic aperture radar imaging method based on digital signal processing according to claim 3, characterized in that: Range scaling transformation and transposition processing are performed on the motion-compensated sub-data to obtain range-transposed data, which specifically includes: The motion compensation sub-data is multiplied by a scale function to obtain scale function processed data; the scale function is: Among them, Φ scl (τ) is the scale function, τ is the instantaneous fast time in the range direction, j is a complex number, kr is the frequency modulation slope, δ r is the distance-to-scale transformation factor, Φ ref is the pitch angle of the synthetic aperture radar to the center of the scene, θ is the instantaneous azimuth angle of the synthetic aperture radar during flight, Φ is the instantaneous pitch angle of the synthetic aperture radar from the center of the scene during flight, t r is the distance to fast time, R a is the instantaneous distance between the synthetic aperture radar and the scene, c is the speed of light; The scale function processed data is subjected to a one-dimensional fast Fourier transform and multiplied by a discrete filter to obtain filtered processed data; the discrete filter is Among them, H2 is a discrete filter, fr is the range frequency, fc is the carrier frequency, and R0 is the effective range of the synthetic aperture radar; The filtered data is transformed into the time domain by using the one-dimensional inverse fast Fourier transform method and multiplied by the inverse scaling function to obtain range-transposed data; the inverse scaling function is: Among them, Φ ins is the inverse scaling function.

6. The terahertz synthetic aperture radar imaging method based on digital signal processing according to claim 3, characterized in that: According to the following formula, azimuth interpolation processing is performed on the range-transposed sub-data to obtain azimuth-processed data: Among them, g(x) is the azimuth processing data, g d (i) is the distance transposed sub-data, and sinc is the convolution kernel.

7. The terahertz synthetic aperture radar imaging method based on digital signal processing according to claim 3, characterized in that: One-dimensional autofocus processing is performed on the initial imaging data to obtain the final imaging data, which specifically includes: For the m-th iteration, according to the maximum energy criterion, the strongest scattering point in each range cell of the initial imaging data is determined, and circular shifting is performed on each strongest scattering point so that the strongest scattering point in each range cell moves to the image center to obtain center-shifted data; 0 < m < M, where M is the iteration threshold; Windowing processing is performed on the center-shifted data to obtain Hamming window-processed data; Inverse fast Fourier transform is performed on the Hamming window-processed data to obtain inverse fast Fourier transform data; The conjugate of the inverse fast Fourier transform data is obtained to obtain conjugate-processed data; The inverse fast Fourier transform data is multiplied by the conjugate-processed data to obtain a phase error gradient; According to the phase error gradient, the phase error value is determined; The phase error value is multiplied by each strongest scattering point to obtain corresponding compensation-processed data; The corresponding strongest scattering point is replaced with the compensation-processed data; When the number of iterations is greater than or equal to the iteration threshold, the initial imaging data is transformed into the phase history domain by using the one-dimensional inverse fast Fourier transform method and multiplied by the current phase error value to obtain error imaging data; The error imaging data is transformed into the image domain by using the fast Fourier transform method to obtain the final imaging data.

8. The terahertz synthetic aperture radar imaging method based on digital signal processing according to claim 7, characterized in that: Windowing processing is performed on the center-shifted data to obtain Hamming window-processed data, which specifically includes: The center shifted data is multiplied by a Hamming window function to obtain Hamming window processed data; the Hamming window function is: Where w(n) is the Han window function and w is the window width.

Images