GPU-Parallel Accelerated Method for Space-Based ISAR Autofocusing Fast Motion Compensation

Abstract

The present invention provides a space-based ISAR autofocusing fast motion compensation method accelerated by GPU parallelization. The method utilizes echo data and angle measurement information of space targets measured by space-based radar, and utilizes GPU parallel processing architecture to synchronously calculate and fit radar angle information. Then, by parallelizing polar coordinate resampling, loop statements are reduced, front-to-back coupling correlation is reduced, and efficient two-dimensional motion decoupling is achieved. Then, a polar coordinate format algorithm and a motion error inversion scheme are combined to perform distance downsampling on the signal after two-dimensional motion decoupling, and valid data segments are intercepted to reduce the data volume. Phase gradient autofocusing and image minimum entropy autofocusing are used in stages to estimate the non-space-varying error of the distance. Then, the GPU parallel processing architecture is used to implement parallel processing of phase gradient autofocusing and image minimum entropy autofocusing, and space-varying error is interpolated and estimated using the estimated error, thereby improving estimation accuracy and computational efficiency, and completing robust focused imaging of space targets.

Description

Technical Field

[0001] The present invention relates to the technical field of radar image processing, and in particular to a space-based ISAR autofocus fast motion compensation method accelerated by GPU parallelization. Background Art

[0002] Inverse synthetic aperture radar (ISAR) plays a vital role in the observation of aviation and aerospace targets due to its all-day, all-weather, high-resolution, and long-range capabilities. In space-based ISAR imaging of space targets, the high-speed motion of both the radar and the target complicates the imaging geometry. This leads to strong two-dimensional motion coupling and phase variation in the imaging model, which in turn reduces motion compensation accuracy and causes space-variant defocusing of the imaging results. Furthermore, due to limited satellite payload and low radar power, the received radar echoes suffer from low signal-to-noise ratio when the target is far from the radar. This further reduces the accuracy of space-based ISAR motion compensation and migration correction, resulting in defocusing of the imaging results and a significant reduction in image quality. Therefore, efficiently and accurately implementing two-dimensional motion compensation for space targets has become a key technology that must be addressed to improve imaging quality.

[0003] In the prior art, methods for turntable model migration correction generally include spherical transformation and polar coordinate resampling methods. Traditional spherical transformation achieves migration correction for the distance dimension by means of scale transformation. However, there is strong two-dimensional motion coupling in the ISAR imaging scenario of space targets, and there is still defocusing in the azimuth dimension after the spherical transformation. Therefore, the spherical transformation is not suitable for this scenario. Polar coordinate resampling is usually used to decouple the two-dimensional motion. Polar coordinate resampling is generally implemented by performing one-dimensional interpolation on the distance dimension and the azimuth dimension respectively. Through multiple loop iterations, the distance dimension and the azimuth dimension are interpolated to achieve two-dimensional motion decoupling. This is inefficient, and the algorithm has low parallelism and high time complexity.

[0004] In the prior art, residual translation compensation is achieved by estimating the non-space-variant distance error using image minimum entropy autofocusing after distance downsampling, then estimating the space-variant error through interpolation, and finally achieving residual translation compensation. Directly estimating the non-space-variant error term using image minimum entropy autofocusing after distance downsampling is prone to falling into local optimality, resulting in inaccurate estimations and affecting subsequent interpolation results and imaging quality. Furthermore, using image minimum entropy autofocusing results in a high number of iterations and low efficiency. Summary of the Invention

[0005] To address the above-mentioned problems in the prior art, the present invention provides a GPU-accelerated space-based ISAR autofocus fast motion compensation method. The technical problem to be solved by the present invention is achieved through the following technical solutions:

[0006] The present invention provides a GPU parallel-accelerated space-based ISAR autofocus fast motion compensation method, comprising:

[0007] S100, receiving echo data of space targets and space-based radar angle measurement information in a low signal-to-noise ratio environment, converting all echo data into the wavenumber domain to obtain echo data in the wavenumber domain, and assembling all echo data into an original echo matrix;

[0008] S200, using a GPU parallel processing architecture to concurrently fit the target's rotation angle at each pulse moment from the angle measurement information to form rotation angle information;

[0009] S300 uses the GPU parallel processing architecture to perform resampling processing on the echo matrix of the current iteration in polar coordinates to obtain a two-dimensional interpolated spectrum;

[0010] S400, intercepting the range data in the spectrum after the two-dimensional interpolation and downsampling it, then selecting the valid data segment in the azimuth direction and assigning values ​​thereto to obtain a first echo matrix;

[0011] S500, inputting the first echo matrix into the GPU parallel processing architecture, and estimating the range non-space-varying error of the first stage in the GPU parallel processing architecture;

[0012] S600: Input the non-space-varying distance error of the first stage into the GPU parallel processing architecture, interpolate and solve the two-dimensional space-varying error phase in the GPU parallel processing architecture, and use the two-dimensional space-varying error phase to roughly compensate the echo matrix in S500 to obtain the echo matrix of the current cycle, which is used as the echo matrix input of the next cycle, and repeat the process from S500 to S600 until the number of loop iterations is reached to obtain the echo matrix after the compensation of the first stage;

[0013] S700, inputting the echo matrix after the compensation in the first stage into the GPU parallel processing architecture, and estimating the range non-space-varying error in the second stage in the GPU parallel processing architecture;

[0014] At step S800, the range non-space-varying error of the second stage is input into the GPU parallel processing architecture, and the two-dimensional space-varying error phase is solved by parallel interpolation. The two-dimensional space-varying error phase is used to perform fine compensation on the echo matrix in step S700 to obtain the echo matrix after the compensation of the second stage, which is used as the echo matrix for the next iteration.

[0015] S900, repeat the process from S300 to S800 until the number of loop iterations is reached, and obtain a two-dimensional motion-compensated echo matrix, which is subjected to a two-dimensional Fourier transform to obtain a focused image.

[0016] Beneficial effects:

[0017] 1. This invention makes full use of the characteristics of GPU computing cores to parallelize and improve the polar coordinate resampling algorithm. For two-dimensional sinc interpolation, the distance dimension and the azimuth dimension are interpolated separately, which reduces the number of loop statements in each dimension interpolation and reduces the front-to-back coupling correlation. When each thread interpolates, the interpolation point is searched by binary search, which can effectively improve the efficiency and reduce the time consumption of polar coordinate resampling from seconds to milliseconds.

[0018] 2. The present invention considers that if the two-dimensional space-variant motion compensation directly uses the image minimum entropy self-focusing method to estimate the non-space-variant error and inverts it to compensate for the space-variant error, there will be too many iterations, low efficiency, and the possibility of falling into a local optimum. To address this problem, it is proposed to use different self-focusing methods in stages to estimate the non-space-variant error and invert and compensate for the two-dimensional space-variant error. This can effectively improve the estimation accuracy and take advantage of the multi-core feature of the GPU to achieve further efficiency improvement. The total time consumption of two-dimensional space-variant motion compensation is increased from seconds to milliseconds.

[0019] The present invention will be further described in detail below with reference to the accompanying drawings and embodiments. BRIEF DESCRIPTION OF THE DRAWINGS

[0020] Figure 1 A flowchart of a GPU parallel-accelerated space-based ISAR autofocus fast motion compensation method provided by the present invention;

[0021] Figure 2a A schematic diagram of the compensation process of the GPU parallel-accelerated space-based ISAR autofocus fast motion compensation method provided by the present invention;

[0022] Figure 2b A flowchart for the rapid implementation of two-dimensional sinc interpolation with GPU parallel acceleration provided by the present invention;

[0023] Figure 3 The simulated point target scattering model diagram provided by the present invention;

[0024] Figure 4 This is the RD imaging result diagram before two-dimensional motion decoupling of the space-based ISAR echo provided by the present invention after de-skew parametrization and rough translation compensation;

[0025] Figure 5 The frequency domain result diagram after resampling in polar coordinates on MATLAB and GPU provided by the present invention;

[0026] Figure 6 This is the RD imaging result diagram after two-dimensional motion compensation on the MATLAB side and GPU side provided by the present invention. DETAILED DESCRIPTION

[0027] The present invention will be further described in detail below with reference to specific examples, but the embodiments of the present invention are not limited thereto.

[0028] The idea of ​​achieving the purpose of the present invention is: based on radar principles, an ISAR turntable model expression is derived, and the correlation between echo data can be used to effectively accumulate echo energy, which is particularly suitable for estimating target motion parameters in a low signal-to-noise ratio echo environment. The echoes are transformed into the wavenumber domain and polar coordinate resampling is used. By interpolating the echoes in polar coordinates into a rectangular coordinate system, the range and azimuth are decoupled, achieving the effects of rotation compensation and two-dimensional motion decoupling. Based on the principle of polar coordinate resampling, a highly parallelized algorithm is designed and implemented to achieve efficient two-dimensional motion decoupling on the GPU side. After resampling, the two-dimensional motion decoupled signal is range-downsampled, and the valid data segment is intercepted to reduce the data volume. The polar coordinate format algorithm is combined with the motion error inversion technology. Phase gradient autofocusing and image minimum entropy autofocusing are used in stages to combine with the polar coordinate resampled echoes to estimate the non-space-varying range error of the downsampled data. The space-varying error is interpolated and estimated through the estimated error. The iterative loop achieves precise compensation for the remaining translation. Based on the above motion compensation principle, phase gradient autofocusing and image minimum entropy autofocusing are parallelized to achieve efficient motion compensation on the GPU side. Finally, the focused ISAR imaging result is obtained by performing a two-dimensional Fourier transform on the motion compensated data.

[0029] Combine Figure 1 and Figure 2a The present invention provides a GPU parallel accelerated space-based ISAR self-focusing fast motion compensation method, including:

[0030] S100, receiving echo data of space targets and space-based radar angle measurement information in a low signal-to-noise ratio environment, converting all echo data into the wavenumber domain to obtain echo data in the wavenumber domain, and assembling all echo data into an original echo matrix;

[0031] Assuming that the target has completed the de-oblique echo parametrization and horizontal motion compensation, the moving radar echo is received in a low signal-to-noise ratio environment, and the echo data is converted to the wavenumber domain, which is expressed as

[0032] W(k r ,θ)=P(k r )exp(-jk r (xsinθ+ycosθ))exp(-jk r ΔR(θ))(1)

[0033] Among them, P(k r ) represents the signal envelope, k rrepresents the radial wave number, θ represents the rotation angle of the target, (x, y) represents the position of the target scattering point in space, and ΔR(θ) represents the residual translational momentum.

[0034] S200, using a GPU parallel processing architecture to concurrently fit the target's rotation angle at each pulse moment from the angle measurement information to form rotation angle information;

[0035] Assuming that the number of echo points in range is M and the total number of points in azimuth is N, the number of points required for the angle accumulation pulse number of the i-th ISAR imaging is N i , open up four data storage areas in the GPU memory, where the four data storage areas are data storage areas echo_all, M×N i × 8 bytes of data storage area echo_cur, N × 8 bytes of data storage area Cta and N i ×8-byte data storage area Cta_temp, M is the number of points in the range direction of the echo signal, N is the total number of points in the azimuth direction of the echo data, and the number of points of the angle accumulation pulse number required for imaging the i-th ISAR image is N i ;Use the GPU parallel processing architecture to set up (N×8+512-1) / 512 blocks, each block contains at least 512 threads;

[0036] Run multiple threads in parallel to use the angle fitting formula to fit the angle θ at the nth pulse moment from the radar angle measurement information n , and store the calculation results in the data storage area Cta;

[0037] A single block thread can have up to 1024 threads. The cudaMemcpyAsync() function, address offsets, and multi-stream parallel processing are used to assign the device-side Cta to the device-side Cta_temp. GPU operations use double data types to ensure accuracy. echo_all stores echo data collected over a period of time, echo_cur stores the i-th echo data that currently meets the imaging conditions, Cta stores the rotation angle information of spatial targets over a period of time, and Cta_temp stores the rotation angle information required for the i-th ISAR imaging image.

[0038] S300: Using the GPU parallel processing architecture, the echo matrix of the current iteration is resampled in polar coordinates to obtain a two-dimensional interpolated spectrum.

[0039] Resample the echo matrix in polar coordinates, i.e. k r sinθ i →k x ,k r cosθi →k y , the resampling formula is expressed as:

[0040] W(k y ,k x )=W(k r cosθ i ,k r sinθ i )=P(k y )exp(-jk x x-jk y y)exp(-jφ e )

[0041]

[0042] θ i =arctan(k x / k y )

[0043] where φ e represents the error phase caused by the residual translation, φ r0 (k y0 ,k x ) represents the non-space-varying error, φ rd (k y ,k x ) represents the space-variant error, k y0 Indicates the distance to the wave number center, φ r0 and φ rd Only the values ​​of the independent variables are different, so we can r0 Interpolation yields φ rd ,θ i represents the change in rotation angle after resampling of the i-th ISAR image, k x represents the azimuthal wave number, k y represents the distance wave number, and arctan() represents the inverse tangent.

[0044] In a specific embodiment of the present invention, S300 includes:

[0045] S310, construct the matrix size of the range wave number and the azimuth wave number according to the size of the original echo matrix; r and the angle θ at the time of the i-th pulse i Construct the values ​​of the first matrix of distance wavenumbers and the second matrix of azimuth wavenumbers; allocate M×N matrices of size in the GPU global memory respectively i × 8 bytes of distance to the wave number data storage area kx_new and M × N i×8 bytes of azimuth wave number data storage area ky_new; using GPU parallel processing architecture, set (M×N i +512-1) / 512 blocks, each block containing at least 512 threads; using multiple threads to update the values ​​of the first matrix and the second matrix in parallel to obtain the first updated first matrix and the second matrix, that is, to construct the first matrix of range wavenumbers and the second matrix of azimuth wavenumbers;

[0046] Among them, the size of the original echo matrix is ​​M×N i ×8×2 bytes.

[0047] S320, using GPU parallel processing architecture, set ((ceil(log2M×N i )) 2 +512-1) / 512 blocks, where ceil means rounding up, and each block contains at least 512 threads; using the parallel reduction idea, M×N i The data is divided into M×N levels i / 2 groups, each group contains two data sampling points, using the GPU's M×N i / 2 threads parallely calculate the maximum value of each set of data sampling points; for M×N i / 2 maximum values ​​are grouped in pairs, and the solution process for finding the maximum value is iterated cyclically until one group of maximum values ​​remains; the maximum and minimum values ​​of the updated range wavenumber and azimuth wavenumber are solved using the solution for solving the remaining group of maximum values; and the first matrix and the second matrix updated for the first time are updated according to the maximum and minimum values ​​to obtain the first matrix and the second matrix updated for the second time;

[0048] It is worth noting that: setting shared memory kx_temp

[512] , the maximum shared memory can be set to 1024, assigning global memory to shared memory can further improve efficiency. i If it is not a multiple of 2, kx_temp is supplemented. If the maximum value is required, it is supplemented to INT_MIN. Otherwise, it is supplemented to INT_MAX, so as to solve the maximum and minimum values ​​of the distance wave number.

[0049] S330, using the GPU parallel processing architecture to set the same block and thread number as in S320, and combining the first matrix and the second matrix updated for the second time to calculate the radial wave number matrix k r Perform parallel updates; allocate M×N size in GPU global memory i ×8 bytes of data storage area ind, and solve in parallel the new k x -k yThe index value ind in the coordinate system is stored in the data storage area ind; the size of M×N is allocated in the GPU global memory i ×8×2 bytes of data storage area s1; using the GPU parallel processing architecture and sinc interpolation principle, setting the same number of blocks and threads as in S320, combining the index value ind to perform one-dimensional range interpolation on the original echo matrix to obtain an interpolated spectrum, and storing it in data storage area s1;

[0050] It is worth noting that the interpolated spectrum is of type complex double. The formula used for interpolation is:

[0051]

[0052] Among them, s represents the echo spectrum of the original echo matrix, x interp Represents the spectrum after one-dimensional interpolation, i represents the number of interpolation rounds for each point, round() represents the rounding function, sinc represents the Singer function,

[0053] S340 uses the GPU parallel architecture to set the same number of blocks and threads as in S320, and combines the first matrix and the second matrix updated for the second time to calculate the rotation angle θ. i Perform parallel updates to obtain the updated angle theta_new; construct a __device__ function on the device side to implement binary search, and call the __device__ function in the kernel function, so that each thread in the kernel function can perform binary search in parallel to obtain the updated angle theta_new. i The neighboring points on the distance are used as index values, and the index values ​​are linearly interpolated in parallel to obtain a new index value ind1; the new index value ind1 and the spectrum s after distance interpolation are combined to obtain a two-dimensional interpolated spectrum ss, the size of which is M×N i ×8×2; returns the two-dimensional interpolated spectrum ss.

[0054] It's worth noting that, while the interpolated points change nonlinearly after distance interpolation, the change between adjacent points is small. Therefore, the change between adjacent points can be approximated as linear. This can be achieved by constructing a __device__ function on the device side to implement a binary search. After returning the two-dimensional interpolated spectrum ss, unneeded memory allocated on the GPU is released.

[0055] S400, intercepting the range data in the spectrum after the two-dimensional interpolation and downsampling it, then selecting the valid data segment in the azimuth direction and assigning values ​​thereto to obtain a first echo matrix;

[0056] In a specific embodiment of the present invention, S400 includes:

[0057] S410, intercepting a portion of data close to the distance wavenumber center in the two-dimensional interpolated spectrum, and downsampling the portion to obtain downsampled data;

[0058] S420 , selecting valid data segments in azimuth from the downsampled data, and assigning values ​​to them in parallel by means of kernel function indexing to obtain a first echo matrix.

[0059] It is worth noting that by changing the φ of formula (3) e In k y =k y0 After Taylor expansion, it can be found that the closer the distance is to the wave number center, the larger the two-dimensional space-variant error φ is. rd The smaller it is, the smaller the k y0 Part of the data is downsampled, and then the valid data segment in the azimuth is selected to further reduce the amount of matrix data. The parallel assignment is realized by the kernel function index to obtain a new downsampled echo matrix yPol_new, whose memory space size is Nr×Na×8×2 bytes, where Nr represents the number of range points after distance downsampling, and Na represents the number of azimuth points after the valid data segment in the azimuth is intercepted, which is recorded as the first echo matrix.

[0060] S500, inputting the first echo matrix into the GPU parallel processing architecture, and estimating the range non-space-varying error of the first stage through the GPU parallel processing architecture;

[0061] In a specific embodiment of the present invention, reference Figure 2b , S500 includes:

[0062] S510, using the parallel architecture of the GPU to set a kernel function, then input the first echo matrix into the kernel function, and use the parallel reduction idea to select special points in the first echo matrix through the kernel function, and form the special points into a second echo matrix; at this time, the size of the second echo matrix becomes Nr1×Na, where Nr1 represents the number of special points.

[0063] S520, allocate a data storage area X of Nr1×Na×8×2 bytes and a data storage area Hpha of Na×8×2 bytes in the GPU global memory, where the data storage area Hpha is used to store the distance non-space-varying error; use the set interpolation multiplier k to interpolate the second echo matrix to obtain a third echo matrix, and store it in the data storage area Y; the third echo matrix is ​​expressed as [0(X.*Hpha)0], where 0 is an all-zero matrix with Nr1 rows and (k-1)×Na / 2 columns.

[0064] S530, calling the FFT function of the CUFFT library to perform FFT transformation on the third echo matrix along the azimuth dimension to obtain a fourth echo matrix Tmp, and then performing Fourier transform on the fourth echo matrix Tmp in the azimuth dimension;

[0065] This step calls cufftHandle in the CUFFT library to formulate an FFT plan, sets the parameter nx to the number of azimuth points Na, sets the parameter batch to the number of range points Nr1, and performs Fourier transform on the fourth echo matrix in the azimuth dimension, expressed as Tmp = CUDA_FFT(Y,2).

[0066] S540, using the parallel architecture and parallel reduction concept of the GPU to select a strong distinctive point in the fourth echo matrix; performing a circular shift on the original echo matrix, and moving the peak of the strong distinctive point to the image center corresponding to the original echo matrix to obtain a fifth echo matrix YImg, the matrix size of which is Nr1×Na;

[0067] This step selects strong distinctive points based on the same parallel reduction idea as in S320.

[0068] S550: Input the fifth echo matrix into the GPU parallel architecture and calculate the sixth echo matrix YImgnew, which is expressed as

[0069]

[0070] Where, CUDA_IFFT([],2) indicates inverse Fourier transform along the azimuth direction, Win indicates window function, and conj() indicates conjugation.

[0071] S560: Using the GPU parallel architecture and the parallel reduction concept, set a block size of Nr1, with a thread size of Na under each block; solve the mean and correlation coefficient for the same distance unit data in each block, and solve the weight in parallel based on the results; use the sixth echo matrix to solve the distance non-space-varying error in parallel for the threads in each block;

[0072] This step uses the GPU parallel architecture and follows the parallel reduction principle. The block size is set to Nr1 and the thread size is Na, where Na <= 1024. The mean, correlation coefficient, and other parameters are calculated for the same distance cell data within each block, and the weights are calculated in parallel based on the results. If Na > 1024, the mean and other parameters need to be solved iteratively. The non-null-varying error function for the parallel distance of threads within each block is expressed as:

[0073]

[0074] In the formula, α represents the weight;

[0075] S570, iterate S520-S560 until the set number of cycles is reached, and output the distance non-space-varying error of the first stage.

[0076] S600: Input the non-space-varying distance error of the first stage into the GPU parallel processing architecture, interpolate and solve the two-dimensional space-varying error phase in the GPU parallel processing architecture, and use the two-dimensional space-varying error phase to perform coarse compensation on the echo matrix in S500 to obtain the echo matrix of the current cycle, which is used as the echo matrix input of the next cycle, and repeat the process from S500 to S600 until the number of loop iterations is reached to obtain the echo matrix after the compensation of the first stage;

[0077] In a specific embodiment of the present invention, reference Figure 2b , S600 includes:

[0078] S610, taking advantage of the GPU's multi-stream parallel processing capabilities, transfers the first-stage distance non-empty error Hpha estimated by the device side back to the host side through the cudaMemcpyAsync() function. The use of multi-stream data asynchronous copy can improve efficiency.

[0079] In step S620, the cubic spline interpolation function of the gsl library is called. The non-space-variant distance error Hpha from the first stage and the rotation angle information before and after resampling are combined to interpolate the non-space-variant distance residual translation error and the two-dimensional space-variant residual translation error before interception in step S400 to obtain an interpolation result. The interpolation result is then transmitted back to the GPU end with the multi-stream parallel function in step S610.

[0080] S630, using the GPU parallel processing architecture to write kernel functions, thread blocks, thread sizes, and construct the residual translation error compensation function H CEr ;

[0081] This step constructs the residual translation error compensation function H according to the following formula CEr :

[0082]

[0083] Where ΔR(arctan(k x / k y )) can be expressed by ΔR(arctan(k x / k y0 ))Interpolation is obtained, φ r0_interp (k y0 ,k x ) is the distance non-space-varying error φ r0 (k y0 ,k x ) interpolation to obtain the two-dimensional space-varying error phase before downsampling.

[0084] S640, using the GPU parallel processing architecture to write a kernel function, set the same thread block and thread size as S630, and compensate the first echo matrix in echo_cur in parallel for each block of threads to obtain the echo matrix of the current cycle, which is used as the echo matrix input for the next cycle. Repeat the process from S500 to S600 until the number of loop iterations is reached to obtain the echo matrix after the first stage compensation.

[0085] S700, inputting the echo matrix after the compensation in the first stage into the GPU parallel processing architecture, and estimating the range non-space-varying error in the second stage in the GPU parallel processing architecture;

[0086] In a specific embodiment of the present invention, S700 includes:

[0087] S710, in the current cycle, using the echo matrix, the number of azimuth points, and the number of range points Nr of the current cycle to concurrently solve various parameters required for the Newton method optimization solution;

[0088] During the loop, the present invention utilizes the GPU parallel processing architecture to write a kernel function, set the thread size to 512, the thread block size to (Nr × Na + 512 - 1) / 512, and set up parallel matrix transposition operations. It then calls cufftHandle in the CUFFT library to formulate the FFT plan, setting the parameter nx to the number of azimuth points Na and the parameter batch to the number of range points Nr, and performing a Fourier transform on the azimuth dimension. The kernel function solves the parameters required by the Newton method in parallel, with the summation function implemented using the parallel reduction concept described above.

[0089] S720 uses the GPU parallel processing architecture and Newton method to optimize the required parameters, and uses the threads of each block to solve the Newton method gradient descent direction in parallel to obtain the non-space-varying distance error of the current cycle;

[0090] S730, determining whether the cycle is completed. If not, compensating the first echo matrix using the non-space-varying error of the distance in the current cycle to obtain the first echo matrix in the next cycle;

[0091] S740, repeat S710-S730 until the loop ends to obtain the first echo matrix of the end loop;

[0092] S750: Input the first echo matrix after the loop is completed into the GPU parallel processing architecture, and estimate the range non-space-varying error of the second stage in the GPU parallel processing architecture.

[0093] At step S800, the range non-space-varying error of the second stage is input into the GPU parallel processing architecture, and a two-dimensional space-varying error phase is solved by parallel interpolation. The two-dimensional space-varying error phase is then used to precisely compensate the echo matrix in step S700 to obtain the echo matrix after compensation in the second stage.

[0094] This step is similar to S600 in solving the two-dimensional space-varying error phase φ rd (k y ,k x ) process is the same, and the processing process will not be repeated here.

[0095] S900, repeat the process from S300 to S800 until the number of loop iterations is reached, and obtain a two-dimensional motion-compensated echo matrix, which is subjected to a two-dimensional Fourier transform to obtain a focused image.

[0096] The present invention fully utilizes the multi-core computing capability of GPUs and improves the parallel design of the algorithm. It parallelizes and improves the two-dimensional sinc interpolation method algorithm for polar coordinate resampling in two-dimensional space-variable translation compensation for ISAR. By performing one-dimensional interpolation on the range and azimuth dimensions, loop statements are reduced, and the front-to-back coupling correlation is lowered. At the same time, a binary search method is used to query the interpolation points, thereby improving efficiency and reducing the time consumption of polar coordinate resampling from seconds to nearly milliseconds.

[0097] Considering the problem that two-dimensional space-varying motion compensation directly uses image minimum entropy self-focusing to estimate non-space-varying errors, it is easy to fall into local optimality, resulting in inaccurate compensation and low efficiency. The present invention can improve the estimation efficiency and estimation accuracy by downsampling the distance direction, and at the same time adopts a staged method to estimate the non-space-varying error. After the first resampling, the phase gradient self-focusing method is used to roughly estimate and roughly compensate the non-space-varying error. Subsequently, the image minimum entropy self-focusing method is used for estimation, which can improve the estimation accuracy and efficiency. At the same time, the effective data segment in the azimuth direction is intercepted to reduce the amount of data and improve the estimation efficiency. The above algorithm is also improved by parallelization, and the characteristics of multiple GPU computing cores are utilized to achieve further improvement in efficiency. The total time consumption of two-dimensional space-varying motion compensation is increased from seconds to milliseconds. The present invention adopts a multi-stream parallel processing method to perform multi-stream parallel design on parameter generation and host-side and device-side data transmission, thereby improving operation and processing efficiency.

[0098] The effectiveness of the present invention can be further illustrated by the following simulation.

[0099] Simulation content and result analysis:

[0100] Emulated device: CPU model is Intel i9-10885H, GPU model is NVIDIA RTX 5000 with Max-Q Design.

[0101] Simulation experiment: The simulated radar parameters are shown in Table 1, and the simulation results are shown in Table 2, Table 3, Figures 4 to 6 .

[0102] Table 1 Simulated radar parameters

[0103]

[0104] Table 2 Comparison of entropy values ​​of imaging results after two-dimensional motion compensation

[0105] Polar coordinate resampling interpolation method 2D sinc interpolation MATLAB image entropy 7.5396 GPU image entropy 7.5398

[0106] Table 3 Comparison of two-dimensional autofocus motion compensation time between MATLAB and GPU

[0107]

[0108]

[0109] Figure 3 This is a diagram of the simulated point target scattering model; Figure 4 This is the RD imaging result of the space-based ISAR echo after de-skew parametrization and coarse translation compensation, but before two-dimensional motion decoupling. It can be seen that the general outline of the observed target has been revealed, but the solar panels on both sides of the target are significantly widened and rotated, and there is still significant defocus at various points. Figure 5 Figure (a) shows the frequency domain result after polar coordinate resampling on the MATLAB side, and Figure (b) shows the frequency domain result after polar coordinate resampling on the GPU side. It can be found that the two are basically consistent. After subtracting the two, the maximum difference is 0.0027, and the average relative error is about 1.6729e-05%, which is basically not much different. Figure 6 Figure (a) shows the RD imaging result after two-dimensional motion compensation on the MATLAB side, and Figure (b) shows the RD imaging result after two-dimensional motion compensation on the GPU side. It can be seen that the imaging results are basically consistent, the rotation of the solar panels on both sides of the target is greatly alleviated, and each scattering point is relatively clear.

[0110] Table 2 compares the image entropy values ​​after two-dimensional motion compensation on the MATLAB and GPU sides, showing essentially identical results. Table 3 compares the efficiency of two-dimensional autofocus motion compensation. The input echo matrices are all complex double. In the simulation, the input matrix size is 512*1024. The two-dimensional autofocus motion compensation process consists of two polar coordinate resampling passes, four phase gradient autofocus passes, and two image minimum entropy autofocus passes. The GPU side's polar coordinate resampling efficiency increased by approximately 70 times, the phase gradient autofocus efficiency by approximately 7 times, and the image minimum entropy autofocus efficiency by approximately 10 times. The overall process efficiency increased by approximately 16 times, from seconds to milliseconds.

[0111] The above is a further detailed description of the present invention in conjunction with specific preferred embodiments, and the specific implementation of the present invention should not be considered to be limited to these descriptions. For those skilled in the art of the present invention, without departing from the concept of the present invention, several simple deductions or substitutions can be made, which should be considered to fall within the scope of protection of the present invention.

Claims

1. A GPU-accelerated space-based ISAR autofocus fast motion compensation method, characterized by: include: S100, receiving echo data of space targets and space-based radar angle measurement information in a low signal-to-noise ratio environment, converting all echo data into the wavenumber domain to obtain echo data in the wavenumber domain, and assembling all echo data into an original echo matrix; S200, using a GPU parallel processing architecture to concurrently fit the target's rotation angle at each pulse moment from the angle measurement information to form rotation angle information; S300 uses the GPU parallel processing architecture to perform resampling processing on the echo matrix of the current iteration in polar coordinates to obtain a two-dimensional interpolated spectrum; S400, intercepting the range data in the spectrum after the two-dimensional interpolation and downsampling it, then selecting the valid data segment in the azimuth direction and assigning values ​​thereto to obtain a first echo matrix; S500, inputting the first echo matrix into the GPU parallel processing architecture, and estimating the range non-space-varying error of the first stage in the GPU parallel processing architecture; S600: Input the non-space-varying distance error of the first stage into the GPU parallel processing architecture, interpolate and solve the two-dimensional space-varying error phase in the GPU parallel processing architecture, and use the two-dimensional space-varying error phase to roughly compensate the echo matrix in S500 to obtain the echo matrix of the current cycle, which is used as the echo matrix input of the next cycle, and repeat the process from S500 to S600 until the number of loop iterations is reached to obtain the echo matrix after the compensation of the first stage; S700, inputting the echo matrix after the compensation in the first stage into the GPU parallel processing architecture, and estimating the range non-space-varying error in the second stage in the GPU parallel processing architecture; At step S800, the range non-space-varying error of the second stage is input into the GPU parallel processing architecture, and the two-dimensional space-varying error phase is solved by parallel interpolation. The two-dimensional space-varying error phase is used to perform fine compensation on the echo matrix in step S700 to obtain the echo matrix after the compensation of the second stage, which is used as the echo matrix for the next iteration. S900, repeat the process from S300 to S800 until the number of loop iterations is reached, and obtain a two-dimensional motion-compensated echo matrix, which is subjected to a two-dimensional Fourier transform to obtain a focused image.

2. The GPU parallel accelerated space-based ISAR autofocus fast motion compensation method according to claim 1, characterized in that: S200 includes: Four data storage areas are opened in the GPU memory, where the four data storage areas are of size Bytes of data storage area 、 Bytes of data storage area 、 Bytes of data storage area and Bytes of data storage area , is the number of distance points of the echo signal, is the total number of points in the azimuth direction of the echo data, The number of points of angle accumulation pulses required for imaging an ISAR image is ; Utilize GPU parallel processing architecture settings blocks, each block contains at least 512 threads; Run multiple threads in parallel to use the angle fitting formula to fit the first angle from the radar angle information. The angle of the pulse moment and store the calculation results in the data storage area middle; in, Store the echo signals collected within a period of time, Stores the first Echo data, Store the target's turning angle information over a period of time. Storage The rotation angle information required for ISAR imaging.

3. The GPU parallel accelerated space-based ISAR autofocus fast motion compensation method according to claim 2, characterized in that: S300 includes: S310, constructing the matrix size of the range wave number and the azimuth wave number according to the size of the original echo matrix; Hedi i Angle measurement at pulse time Construct the values ​​of the first matrix of distance wavenumbers and the second matrix of azimuth wavenumbers; allocate the size of Bytes of distance to the wave number data storage area and Bytes of azimuth wave number data storage area ; Using GPU parallel processing architecture, set blocks, each block containing at least 512 threads; using the multiple threads to update the values ​​of the first matrix and the second matrix in parallel to obtain the first updated first matrix and the second matrix; S320, using GPU parallel processing architecture, set blocks, among which, Indicates rounding up, each block contains at least 512 threads; using the parallel reduction idea, The data is divided into groups, each group contains two data sampling points, using GPU threads to parallelly calculate the maximum value of each set of data sampling points; The extreme values ​​are grouped in pairs, and the solution process for finding the extreme values ​​is iterated cyclically until one group of extreme values ​​remains; the maximum and minimum values ​​of the updated range wavenumber and azimuth wavenumber are solved using the solution for solving the remaining group of extreme values; and the first matrix and the second matrix updated for the first time are updated according to the maximum and minimum values ​​to obtain the first matrix and the second matrix updated for the second time; S330, using the GPU parallel processing architecture to set the same block and thread number as in S320, and combining the first matrix and the second matrix updated for the second time to calculate the radial wave number matrix Perform parallel updates; allocate a size of 1 in the GPU global memory Bytes of data storage area , and solve in parallel in the new - Index value in the coordinate system and store it in the data store ; Open up a global memory in GPU with a size of Bytes of data storage area ; Using the GPU parallel processing architecture and sinc interpolation principle, set the same block and thread number as in S320, and combine the index value Perform one-dimensional range interpolation on the original echo matrix to obtain the interpolated spectrum , and store it in the data storage area ; S340 uses the GPU parallel architecture to set the same number of blocks and threads as in S320, and combines the first and second matrices updated for the second time to rotate the angle. Perform parallel updates to get the updated corners ; Construct a __device__ function on the device side to implement binary search, and call the __device__ function in the kernel function, so that each thread in the kernel function can perform binary search in parallel. The neighboring points on the grid are used as index values, and the index values ​​are linearly interpolated in parallel to obtain new index values. ; Combine the new index value and the spectrum after distance interpolation , get the spectrum after two-dimensional interpolation , the size is ; Returns the two-dimensional interpolated spectrum .

4. The GPU parallel accelerated space-based ISAR autofocus fast motion compensation method according to claim 1, characterized in that: S400 includes: S410, intercepting a portion of data close to the distance wavenumber center in the two-dimensional interpolated spectrum, and downsampling the portion to obtain downsampled data; S420 , selecting valid data segments in azimuth from the downsampled data, and assigning values ​​to them in parallel by means of kernel function indexing, to obtain a new downsampled first echo matrix.

5. The GPU parallel accelerated space-based ISAR autofocus fast motion compensation method according to claim 4, characterized in that: S500 includes: S510, using the parallel architecture of the GPU to set a kernel function, then input the first echo matrix into the kernel function, and using the parallel reduction concept to select distinctive points in the first echo matrix through the kernel function, and compose the distinctive points into a second echo matrix; S520, the size of the GPU global memory is Bytes of data storage area ,and Bytes of data storage area , where the data storage area Used to store the non-space-varying error of distance; use the set interpolation multiple The second echo matrix is ​​interpolated to obtain a third echo matrix and stored in the data storage area Among them, It represents the number of azimuth points after the effective azimuth data segment is intercepted, recorded as the first echo matrix, Indicates the number of highlighted points; S530, calling the FFT function of the CUFFT library to perform FFT transformation on the third echo matrix along the azimuth direction to obtain a fourth echo matrix , and then the fourth echo matrix Perform Fourier transform in the azimuthal dimension; S540, using the parallel architecture and parallel reduction concept of the GPU to select a strong distinctive point in the fourth echo matrix; performing a circular displacement on the original echo matrix, and moving the peak of the strong distinctive point to the image center corresponding to the original echo matrix to obtain a fifth echo matrix , the matrix size is ; S550, inputting the fifth echo matrix into the GPU parallel architecture to calculate a sixth echo matrix ; S560, using GPU parallel architecture and following the parallel reduction idea, sets the size to The thread size under each block is ; Solve the mean and correlation coefficient of the same distance unit data in each block, and solve the weight in parallel based on the results; Use the sixth echo matrix to solve the distance non-space-varying error in parallel for the threads in each block; S570, loop iterate S520-S560 until the set number of loops is reached, and output the distance non-empty variation error of the first stage .

6. The GPU parallel accelerated space-based ISAR autofocus fast motion compensation method according to claim 1, characterized in that: S600 includes: S610, using the characteristics of GPU multi-stream parallel processing, converts the first stage distance non-space-variant error estimated by the device side into Pass it back to the host side through the cudaMemcpyAsync() function; S620, calling the cubic spline interpolation function of the gsl library, combined with the distance non-space-varying error of the first stage and the angle information before and after resampling, interpolating the distance non-space-variant residual translation error and the two-dimensional space-variant residual translation error before interception in S400 to obtain an interpolation result, and then transmitting the interpolation result back to the GPU end; S630, using GPU parallel processing architecture to write kernel functions, thread blocks, thread sizes, and construct residual translation error compensation functions : S640 uses GPU parallel processing architecture to write kernel functions, sets the same thread block and thread size as S630, and parallelizes threads in each block. The first echo matrix in the current cycle is compensated to obtain the echo matrix of the current cycle, which is used as the echo matrix input of the next cycle, and the process from S500 to S600 is repeated until the number of cycle iterations is reached to obtain the echo matrix after the first stage compensation.

7. The GPU parallel accelerated space-based ISAR autofocus fast motion compensation method according to claim 1, characterized in that: The S700 includes: S710, in the current cycle, using the echo matrix, azimuth points and range points of the current cycle Parallel solution of various parameters required for Newton's method optimization solution; S720 uses the GPU parallel processing architecture and Newton method to optimize the required parameters, and uses the threads of each block to solve the Newton method gradient descent direction in parallel to obtain the non-space-varying distance error of the current cycle; S730, determining whether the cycle is completed. If not, compensating the first echo matrix using the non-space-varying error of the distance in the current cycle to obtain the first echo matrix in the next cycle; S740, repeat S710-S730 until the loop ends to obtain the first echo matrix of the end loop; S750: Input the first echo matrix after the loop is completed into the GPU parallel processing architecture, and estimate the range non-space-varying error of the second stage in the GPU parallel processing architecture.

Images