CN109307865A - A CUDA-based millimeter-wave radar RMA imaging method - Google Patents

Abstract

The invention belongs to the technical field of radar imaging, in particular to a CUDA-based millimeter-wave radar RMA imaging method. The method comprises the steps of: (S2) respectively performing two-dimensional azimuth inverse fast Fourier transform on the horizontal and height dimensions of the three-dimensional scattering data s of the target, and performing parallel translation processing to obtain data s ₁ ; (S3) transforming the data s ₁ in space The non-numerical positions outside the support area of the spectral domain range are set to zero to obtain the data s ₂ ; (S4) the data s ₂ is interpolated in the distance to the wavenumber domain to obtain the data s ₃ ; (S5) the data s ₃ is subjected to three-dimensional inverse fast Fourier Leaf transformation, taking the modulo value of the transformation result, and then taking the maximum scattering value for each dimension of the three-dimensional data for imaging. The invention realizes the innovation of the algorithm acceleration of the array scanning system, solves the disadvantage of slow processing speed under the new system, transfers the dependence of the existing algorithm on the CPU, and realizes the repeated calling and transplantation of the algorithm function.

Description

A CUDA-based millimeter-wave radar RMA imaging method

technical field

The invention belongs to the technical field of radar imaging, and in particular relates to a millimeter-wave radar RMA (Range Migration Algorithm, abbreviation: RMA) imaging method based on CUDA (Unified Computing Device Architecture, abbreviation: CUDA, which is a computing platform launched by graphics card manufacturer NVIDIA). .

Background technique

At present, the GPU (graphics processing unit) radar imaging developed in the world mainly focuses on SAR/ISAR. GPU high-performance computing is a parallel computing technology that has developed rapidly in recent years, because it can provide a fast parallel computing framework for large data processing. , has become a very high data-parallel computing density; supports multi-threading; a general-purpose mathematical coprocessor with multiple computing cores, it not only has a very large computing intensity, but also GPU memory throughput is larger than CPU. In addition, related GPU radar imaging research at home and abroad includes the GPU real-time image processing algorithm using NVIDIA CUDA technology at the Technical University of Madrid for 3-D clothing detection radar in the THz band, corresponding to a scanning area of 50x90cm2, which allows the detection of 6000 The image composed of pixels is refreshed at a rate greater than 8fps. The performance acceleration of traditional time-domain backprojection and Kirchhoff migration algorithms for hidden target imaging is considered. The Bangalore Electronics and Radar Development Agency in Bangalore, India has realized the performance acceleration of traditional time-domain back-projection and Kirchhoff migration algorithms for imaging hidden targets. Since the image intensity is calculated pixel by pixel in these algorithms, the GPU can accelerate other Radar Imaging Algorithms. In addition, when the amount of data is large, Xiong Chao from the National University of Defense Technology transplanted the spectrum analysis and spectrum peak search to the GPU for calculation, which effectively shortened the overall response time of the system.

Existing patent technologies related to radar GPU imaging mainly focus on the realization of subsystems in various special fields. Patent application CN201711246749.4 discloses a GPU-based fast frequency domain backward projection 3D imaging method, which uses 3D frequency domain backward projection. The algorithm uses the graphics processor GPU to use a single-program multi-data instruction mode to schedule its rich hardware resources, and then obtains the three-dimensional SAR image through the three-dimensional inverse Fourier transform, thus greatly improving the backward projection. The efficiency of the algorithm. Patent application CN201710485297.9 discloses a GPU-based synthetic aperture radar image target recognition method, and patent application CN107064930A discloses a GPU-based radar forward-looking imaging method. The observation vector and measurement matrix are obtained from the CPU side. On the GPU side, the CUDA architecture is used to restore the backscattering coefficient vector of the scene target from the observation vector in parallel to complete the radar forward-looking imaging. Patent application CN201711019081.X discloses a single-pixel lidar imaging device and imaging method based on an embedded GPU, mainly describing a single-pixel lidar imaging device based on an embedded GPU. The uniform smooth surface is reconstructed, and the reconstructed three-dimensional smooth surface model is displayed dynamically.

The existing radar GPU imaging algorithm is not very portable, and the program module can only be applied to specific algorithms, and the coverage of the algorithm on the GPU is insufficient, so it is only accelerated in specific steps. For example, Xiong Chao of the National Defense Science and Technology University only performed FFT. The CUDA parallelization of the algorithm, matrix transpose and spectral peak search algorithm does not fully port the signal processing part to the GPU implementation. In the existing near-field array imaging for high-frequency radar, the decrease in wavelength will lead to a decrease in the spacing of array elements and a sharp increase in the number of array elements, resulting in a huge amount of data. The large amount of echo data has higher requirements on the data acquisition system and storage device, and the acquisition time of the signal will be longer, and the efficiency of the original imaging method is low.

SUMMARY OF THE INVENTION

In order to solve the above technical problems, the present invention is based on the NVIDIA CUDA computing platform, and integrates the GPU into the radar near-field array imaging, thereby accelerating the processing of large amounts of data, and has the advantages of relatively stability and higher efficiency. The specific technical solutions are as follows:

A CUDA-based millimeter-wave radar RMA imaging method, comprising the following steps:

(S1) receiving radar echo data, and performing preprocessing, the preprocessing includes performing range correction on the echo data, and then using the CUDA platform library function to perform clutter filtering, and generate target three-dimensional scattering data s;

(S2) respectively perform two-dimensional azimuth inverse fast Fourier transform on the horizontal and height dimensions of the three-dimensional scattering data s of the target, and perform parallel translation processing to obtain data s ₁ ;

(S3) Zeroing the position where the NaN value appears in the data s ₁ outside the support area of the spatial spectral domain range to obtain the data s ₂ ; NaN is the abbreviation of Not a Number in English, and is defined in the IEEE floating point arithmetic standard (IEEE 754), Represents some special value (infinity and not-a-number (NaN)), and the support area represents the shape of the spectral domain.

(S4) Interpolate the distance to the wavenumber domain for the data s ₂ to obtain the data s ₃ ;

(S5) Perform a _three -dimensional inverse fast Fourier transform on the data s3, take a modulo value for the transform result, and then take a maximum scattering value for each dimension of the three-dimensional data for imaging.

Preferably, the two-dimensional orientation fast Fourier transform in the step (S2) is composed of the ifftshift kernel function, the cufft library function, and the fftshift kernel function, which respectively realize the inverse FFT movement, the fast Fourier transform, and the zero-frequency component. Move to the center of the spectrum.

Preferably, the specific process of zeroing the position of the NaN value outside the support area of the spatial spectrum range is realized by the ndgrid kernel function, the getP kernel function and the padarray kernel function, and the ndgrid kernel function realizes the rectangular grid in the N-dimensional space. , the getP kernel function obtains the range frequency in parallel from the electromagnetic plane wave discrete relationship, and the padarray kernel function realizes that the part outside the original data is set to zero.

Preferably, the distance-to-wavenumber domain interpolation for the data s ₂ is implemented by the interp1 kernel function.

Preferably, the three-dimensional inverse fast Fourier transform is implemented by the ifftshift kernel function, the cufft kernel function and the fftshift kernel function.

In order to better understand the present invention, the present invention will be described in detail below in combination with the prior art such as CUDA platform library and MATLAB software.

The method of the invention performs accelerated processing on the original data, such as preprocessing, two-dimensional azimuth fast Fourier transform, zeroing of data outside the spatial spectrum domain, range-to-wavenumber domain interpolation, and three-dimensional inverse fast Fourier transform. The algorithm processing step can enable RMA to achieve the 3D imaging result of the target. The method is to construct a kernel function through the CUDA platform language, use MATLAB to generate a mex file, and create a new function to process the data. After entering the CUDA platform, the original three-dimensional data will be converted into one-dimensional data, and the transformed RMA can also pass through cufft. The three-dimensional data needs to undergo the Fourier transform of the last two dimensions. It only needs to reorder the data, and the dimension that needs to be transformed is moved to The first scale, swap the first dimension and the second dimension; adjust the transformation step value and pointer, that is, change the batch size, the scale of the Fourier transform, and the distance of the batches in the function cufftPlanMany, which together change the pointer handle, After the data is processed through a series of steps, it is then transferred to MATLAB or OPENCV is used for imaging. Figure 2 shows the structure diagram of radar near-field imaging processing.

Using MATLAB to generate mex files and creating new functions to process data involves ".cpp" data input and output files and ".h". The ".cpp" data input and output files are extracted from the real and imaginary parts of the original data, and the main function of the ".cu" file is referenced (the ".cu" file is a file unique to the CUDA platform, which is mainly used to execute the device-side code, where the host and device The terminal code is written in the main function of the ".cu" file), the output processed data, and the data export part. The implementation process of the main function of the ".cu" file is:

The main function of the ".cu" file runs on the CUDA platform and needs to include the .h header file. Where parallelism is required, for loops are used to replace sequential iterations with simultaneous ones, and there should be no connection between indexes. The operations in the kernel function on the Device side can use the operations on the host side, and the functions can be changed according to actual needs, such as changing the memset on the host side to cudaMemset. However, operations on the device side cannot be performed on the host. Functions for taking the maximum and minimum values in CUDA: fmaxf, fminf. First, the three-dimensional data will be converted into one-dimensional data by line on the host side of the CUDA platform, and then use the cudaMalloc function to allocate the corresponding memory to the device side variables, and then use the cudaMemcpy function to transfer the data from the host side to the device side, and write the __global__ kernel function. Or use the existing library functions to process the data on the device side, and then use the cudaMemcpy function to transfer the data from the device side to the host side for output after the processing is completed. Because the data on the device side cannot be printed, it must be transmitted to the host side, and a single data does not need to be passed in and out of the cudaMemcpy function, and can be called directly. And after allocating memory to the variable, the initial value of the variable is generally 0. However, when the variable is superimposed ^on itself, the initial value of the variable will be assigned an initial value of the order of 108. Therefore, if the initial value of the variable needs to be used, it is best to use memset or cudaMemset on the host side or the device side respectively. Assign the initial value to the variable. Each type of NVIDIA GPU has a rated number of threads per block. When setting blocksize, i, j, r in dim3blockSize(i, j, r) is the dimension gridDim of dim3 in the x, y, and z directions The value of .x, gridDim.y, and gridDim.z, their product must be less than or equal to the rated number of threads per block, where j, r can be omitted as 1. Generally use dim3 to define the type of gridSize, Determine the dimensions of gridSize in the x, y, z directions of dim3 blockDim.x, blockDim.y, blockDim.z, x, y, z are the sizes of the data that need to be calculated in the three dimensions, that is, the CPU is the iteration of the for loop frequency. If it is only performed in one x, y, z direction, such as the kernel function that implements the dot product function, just set the blockSize to 1, and specify the gridSize as the dimension of one direction, that is, change the scale of the kernel function. When Debugging, due to some hidden errors when compiling MATLAB, it will not be pointed out in detail. You can compile the ".cu" file in VS, and generally only "Error LNK2019: Unresolved external symbol main, which is in the function __tmainCRTStartup, will appear. referenced" and "error LNK1120: 1 unresolved external command" means that the ".cu" file is fine, but sometimes even if it says "MEX has completed successfully", it will reveal "MATLAB has encounteredan internal problem" and needs to close." (for example, one of the reasons for the error is that the host side uses the device side variable), then you need to change the ".cu" file to a file that can be compiled separately, and then debug it on VS, after confirming that it is correct Change back to the original format. When using MATLAB to generate a mex file, I encountered the error "LIBCMT.lib(excptptr.obj): fatalerror LNK1112: The module computer type 'X86' conflicts with the target computer type 'x64'". This error indicates that the lib library conflicts, solve The solution is to modify the mex command: -lcudart LINKFLAGS="$LINKFLAGS/NODEFAULTLIB:LIBCMT.LIB"-lcufft. Ignore the LIBCMT.lib library and continue compiling. The CPU has a larger area for cache and branch prediction. The GPU is a bunch of computing units with a controller, the program is the same, but the input and output data are different. Because the problem that GPU solves is massively parallel computing. The arithmetic unit (sp, stream processor, stream processor, that is, the ALU unit) does not need to be as complicated as the CPU, and the number is several thousand. Therefore, when using the GPU for simple operations, it can be understood as parallel computing like verilog. For example, when performing simple operations of exchange, it is not necessary to allocate a cache for exchange like using a CPU, and the data can be exchanged directly. If the kernel function is called in the kernel function, an error will be reported at the beginning "error: calling a__global__function("")from a__global__function("")is only allowed on the compute_35architecture or above", then just click the project to open the project properties, change The first setting of CUDA tells CUDA to generate relocatable device code by choosing "yes (-rdc=true)", then it has to tell CUDA to use compute 3.5 by changing the code generation to "compute_35,sm_35", then, the notification runs The time library uses a multithreaded library. If in debug mode, select "Multi-threaded debugging (/MTd)" or in release mode "Multi-threaded (/MT)", and finally in the linker properties, add an additional dependency as "cudadevrt.lib".

The ".h header file" only needs to declare the input and output variables in the ".cu" file, and the function prototype extern "C" can be used to indicate the convention to be used.

The beneficial effects obtained by adopting the invention: 1. The invention realizes the innovation of algorithm acceleration of the array scanning system, and solves the disadvantage of slow processing speed under this new system; 2. The invention can take advantage of the GPU floating-point calculation function, The dependence of the existing algorithm on the CPU is transferred; 3. The present invention generates the mex file by using MATLAB, thereby realizing the repeated invocation and transplantation of the algorithm function for many times.

Description of drawings

Fig. 1 is the flow chart of the method of the present invention;

Figure 2 is a block diagram of radar near-field imaging processing;

3 is a diagram showing the imaging result of the half-body mannequin by the method of the present invention in the embodiment, wherein the figure (a) is the original image of the half-body mannequin, the figure (b) is the front view of the imaging result of the half-body mannequin, and the figure (c) is the half-body mannequin. Side view of mannequin imaging results.

Detailed ways

Hereinafter, the present invention will be further described with reference to the accompanying drawings and embodiments.

Figure 1 is a flowchart of the present invention, a CUDA-based millimeter-wave radar RMA imaging method, comprising the following steps:

(S1) Receive radar echo data, and perform preprocessing. The preprocessing includes performing range correction on the echo data, and then using the inverse transformation parameter CUFFT_INVERSE in the cufft function cufftExecC2C of the CUDA platform library to convert the range to the range beyond the actual target range. Set the data to zero to remove the influence of clutter, and then use the forward transformation parameter CUFFT_FORWARD in the library cufft function cufftExecC2C to obtain the target three-dimensional scattering data s(x, y, f), where x, y, f are the horizontal, height, and frequency dimensions, respectively. The data preprocessing in the example consists of a linspace kernel function, an Equal kernel function, and a multwith kernel function, which are used to generate linearly spaced arrays in parallel, to remove the last value of the array, and to multiply the data and the correction item in parallel. The two kernel functions of linspace and Equal are used to generate the actual distance sequence. The actual distance sequence is used to limit the data of the echo within the actual distance range. Before performing inverse fast Fourier transform and forward transformation on the distance direction, the multwith kernel is used for this limiting process. The data after the function is executed is then Fourier transformed.

(S2) Perform two-dimensional azimuth inverse fast Fourier transform on the horizontal and height dimensions of the three-dimensional scattering data s of the target respectively, and perform parallel translation processing to obtain data s ₁ ; specifically, use the platform library cufft to respectively transform the scattering data processed in the previous step by using the platform library cufft The horizontal and height dimensions of the data s(x, y, f) are subjected to two-dimensional azimuth inverse fast Fourier transform, and the kernel function that implements the fftshift function in matlab is written to perform parallel translation of the horizontal and height dimensions, respectively, to obtain scattering data. Horizontal and height dimension data s ₁ (X,Y,f).

The formula for calculating the zero-IF signal s(X,Y) is as follows:

where f(x,y) is the target scattering characteristic function, D is the target area, Z ₀ indicates that the target position is on the z=Z ₀ plane, k is the wavenumber, k=2π/λ, λ is the signal wavelength, k _x ,k _y , k _z are the wavenumber components of the wavenumber k along the x, y, and z directions in turn, FT2 represents the two-dimensional Fourier transform, and IFT2 represents the two-dimensional inverse Fourier transform.

(S3) zeroing out the position where the NaN value appears in the data s ₁ outside the support area of the spatial spectrum domain, to obtain the data s ₂ (X, Y, f);

(S4) Interpolate the distance to the wavenumber domain for the data s ₂ to obtain the data s ₃ ;

First, use downward rounding to get the output (frequency dimension, for example, the step value has been reduced by 10 times, and the frequency point has been expanded by 10 times) coordinate F _m (F _m represents the mth of the input frequency F after the step value is reduced by ten times value) in the relative position of the input (frequency dimension) coordinate f _m after the above operations, that is, d _m =(int)F _m -f _m , and then according to the relative position in the adjacent input data (scattering data after the above steps) )s ₂ (X, Y, f) distance weight value, i.e. i=blockIdx.x*blockDim.x+threadIdx.x; j=blockIdx.y*blockDim.y+threadIdx.y; r=blockIdx.z *blockDim.z+threadIdx.z; set idx1, idx0, idx22 to represent the coordinate index number, because the three-dimensional data will become one-dimensional data in CUDA, and the coordinate index number is used to specify the index value in the data.

idx22=sy·(sz·i+r)+j,

idx0=sx·(sy·r+j)+d _m ,

idx1=sx·(sy·r+j)+ _dm +1;

Then, s ₃ (X,Y,f) _idx22 =wd·s ₂ (X,Y,f) _idx0 +wu·s ₂ (X,Y,f) _idx1 , where sx,sy,sz are the three-dimensional data Scale, blockDim and blockIdx are the dimensions and coordinates of the thread block operation in the x, y, and z directions, threadIdx is the x, y, and z direction coordinates of the thread operation in the thread block, and wd and wu are the distance weights. Using the kernel function to perform operations at the same time, you can Corrects for distance curvature.

(S5) respectively perform cufft on the data s ₃ (X, Y, f) after the above processing in the three dimensions of height, level and frequency, and when executing the corresponding dimension, the data needs to be moved to the first dimension for calculation, For example, the height dimension is in the first dimension, and the horizontal dimension is in the second dimension. When performing horizontal cufft, you need to exchange the height of the horizontal dimension and the first dimension, and then set the pointer and step value to get the scattering data s ₂ (x,y,t) . Take the absolute value of the scattering data s ₂ (x, y, t), and then take the maximum scattering value for each dimension of the three-dimensional data for imaging. Because the imaging function imagesc in matlab only accepts the input value as a real number, and the transformation result s ₂ (x, y, t) is a complex number, the modulo value is taken. Because the image formed by imagesc is a two-dimensional image of two-dimensional data, the maximum scattering value of each dimension of the three-dimensional data s ₂ (x, y, t) is taken for imaging separately.

As shown in FIG. 3 , it is the image of the imaging result of the half-body mannequin using the method of the present invention in the embodiment (coordinates indicate the position, unit: m), the original data is the complex scattering data of the frequency, height and horizontal three-dimensional scanning, and then the data is processed RMA processing, this processing achieves the function of GPU-accelerated imaging algorithm by building a kernel function on the CUDA platform and processing the data in parallel. The average imaging original time of 5 times is MATLAB2017b on CPU E5-2683 takes 17.839 seconds, and TITAN XP GPU takes 6.134 seconds, and the speed is increased by 2.9 times. The method of the present invention has verified the measured data through the above experiments, and the effect is ideal.

Claims (5)

1. a CUDA-based millimeter wave radar RMA imaging method, is characterized in that comprising the following steps: (S1) receiving radar echo data, and performing preprocessing, the preprocessing includes performing range correction on the echo data, and then using the CUDA platform library function to perform clutter filtering, and generate target three-dimensional scattering data s; (S2) respectively perform two-dimensional azimuth inverse fast Fourier transform on the horizontal and height dimensions of the three-dimensional scattering data s of the target, and perform parallel translation processing to obtain data s ₁ ; (S3) zeroing the position where the data s ₁ appears non-numerical outside the support region of the spatial spectral domain, to obtain the data s ₂ ; (S4) Interpolate the distance to the wavenumber domain for the data s ₂ to obtain the data s ₃ ; (S5) Perform a _three -dimensional inverse fast Fourier transform on the data s3, take a modulo value for the transform result, and then take a maximum scattering value for each dimension of the three-dimensional data for imaging. 2. a kind of CUDA-based millimeter-wave radar RMA imaging method as claimed in claim 1, is characterized in that: the two-dimensional azimuth fast Fourier transform in described step (S2), by ifftshift kernel function, cufft library function , fftshift kernel function, respectively realize inverse FFT shift, fast Fourier transform, and move the zero-frequency component to the center of the spectrum. 3. a kind of CUDA-based millimeter wave radar RMA imaging method as claimed in claim 1, is characterized in that: the concrete process that the position of non-numerical value appears outside described space spectral domain support area is zeroed by ndgrid kernel function, getP The kernel function and the padarray kernel function are realized, the ndgrid kernel function realizes a rectangular grid in the N-dimensional space, the getP kernel function obtains the range frequency in parallel by the electromagnetic plane wave discrete relationship, and the padarray kernel function realizes the non-numerical value outside the original data. position is set to zero. 4 . The CUDA-based millimeter-wave radar RMA imaging method according to claim 1 , wherein the interpolation of the data s ₂ in the range to the wavenumber domain is realized by the interp1 kernel function. 5 . 5 . The CUDA-based millimeter-wave radar RMA imaging method according to claim 1 , wherein the three-dimensional inverse fast Fourier transform is realized by ifftshift kernel function, cufft kernel function and fftshift kernel function. 6 .