Spaceborne SAR scene-matched curve imaging multi-gpu time-domain fast processing method

Abstract

The application discloses a spaceborne SAR scene matching curve imaging multi-GPU time domain fast processing method, establishes a terrain matching grid division sequence model based on spaceborne curve SAR echo data designed through a curve imaging planning system, carries out multi-GPU parallel pulse compression and back projection processing on the echo data, obtains coarse images of each sub-aperture, then carries out two-stage spectrum compression-decompression parallel processing on the coarse images obtained through processing through multi-GPU, obtains fine images of each sub-aperture, and effectively reduces the data volume, expands the distributed parallelism, and equivalently realizes the fast imaging effect of large data volume echo data. The application solves the defects that the traditional synthetic aperture radar back projection algorithm has high calculation complexity and leads to long processing time in the scene matching curve imaging mode, and has a wide market application prospect.

Description

Technical Field

[0001] This invention belongs to the field of Synthetic Aperture Radar (SAR) technology, specifically relating to a fast multi-GPU temporal processing method for spaceborne SAR scene matching curve imaging. Background Technology

[0002] Spaceborne SAR scene matching curve imaging significantly improves the observation timeliness of long curve scenes by customizing and generating imaging bands with matching scenes and uniform azimuth resolution, making it an effective method for monitoring maritime routes. However, when using this mode for observation and imaging, the observation geometry is complex, and the slant range and viewing angle vary greatly, resulting in increased echo delay and viewing angle variations.

[0003] Back projection (BP) algorithms are suitable for complex observation geometries and can achieve high-precision spaceborne SAR imaging without approximations. Spaceborne curved SAR has a large swath width / resolution ratio, and orbital curvature and surface curvature make the signal spatially complex. Back projection algorithms can be used to address this problem and obtain high-precision focused SAR images.

[0004] To address the challenges of nonlinear trajectory, large echo data volume, high computational complexity in time-domain imaging, and difficulty in rapid processing of spaceborne curved synthetic aperture radar (SAR), many researchers have improved the BP imaging algorithm based on hierarchical techniques, decomposing SAR images into multiple sub-aperture images, also known as sub-images. However, for spaceborne curved SAR data, these methods still suffer from high computational load and significant time consumption when calculating each level of sub-aperture image. Furthermore, when the number of sub-apertures is large, serial processing of each sub-aperture image is also very time-consuming. These problems limit the practical application of hierarchical BP algorithms. Therefore, it is necessary to develop a fast time-domain imaging processing method based on multi-GPUs for spaceborne curved SAR imaging. This method utilizes multi-GPU parallelism and high-concurrency threads to efficiently and accurately achieve high-resolution SAR imaging, possessing significant engineering value and application prospects. Summary of the Invention

[0005] In view of this, the present invention provides a multi-GPU fast temporal processing method for scene matching curve imaging of spaceborne SAR. It realizes the design of a multi-GPU parallel processing architecture for curve imaging mode under complex observation geometry, solves the problem of long processing time of echo imaging using temporal imaging processing algorithms for curve SAR, and can meet the needs of fast imaging processing of large data volume echoes of spaceborne curve SAR.

[0006] The technical solution of this invention is as follows: A fast temporal processing method for spaceborne SAR scene matching curve imaging using multiple GPUs includes: Step 1: Use the scene matching curve imaging echo data of spaceborne SAR to obtain a fast time-domain imaging processing grid division sequence; Step 2: Perform pulse compression and back projection parallel processing on the echo data of each sub-aperture under each GPU to obtain the coarse image under that sub-aperture; Step 3: Perform two-stage spectral compression-decompression full-process imaging parallel pipeline processing on the coarse image to obtain fine images of each sub-aperture; Step 4: Use multiple GPUs in parallel to perform data sub-block fusion processing and geographic coordinate interpolation projection on each sub-aperture fine image to finally obtain the SAR sub-image.

[0007] Beneficial effects: (1) The present invention adopts a software implementation method of CPU multi-threading and multi-GPU collaboration. While expanding the parallel processing capability of multiple nodes, it reduces data redundancy and memory occupation between multiple GPUs, avoids the problem of large memory resource consumption in traditional multi-GPU solutions, thereby improving system resource utilization and realizing large-scale data processing capability.

[0008] (2) By designing the sub-aperture imaging process in parallel, this invention realizes distributed parallel computing of modules such as echo data pulse compression, coarse image back projection, frequency domain transformation and phase compensation, which effectively improves the processing efficiency of spaceborne curve SAR time domain imaging and solves the problem of low processing speed of traditional single GPU.

[0009] (3) This invention addresses the problem of high computational complexity and long processing time of traditional BP algorithm in scene matching curve SAR imaging mode. By constructing a fast time domain imaging processing flow, it achieves fast imaging of large amount of echo data while ensuring imaging processing capabilities. It has good engineering application value and market application prospects. Attached Figure Description

[0010] Figure 1 Geometric model for fast temporal imaging processing of spaceborne curve SAR; Figure 2 A flowchart illustrating the principle of a fast temporal domain imaging processing method for spaceborne curve SAR. Figure 3 Diagram of the parallel processing architecture for multi-GPU curve SAR imaging; Figure 4 Flowchart of spaceborne curve SAR scene matching and mesh generation; Figure 5 The projection process corresponding to each grid point is parallelized, and each grid point is independently imaged. Figure 6 A schematic diagram of a single-GPU imaging pipeline for spaceborne curve SAR. Figure 7A comparison of the echo data processing and imaging effects of the method described in this invention and the traditional time-domain BP time-domain imaging algorithm; Figure 8 The image shows the curve imaging effect of the fast time-domain imaging algorithm for spaceborne curve SAR. Detailed Implementation

[0011] The present invention will now be described in detail with reference to the accompanying drawings and embodiments.

[0012] In the fast temporal imaging mode of spaceborne curve SAR, traditional temporal imaging processing methods do not consider curve maneuver constraints in their design. This leads to problems such as decreased focusing accuracy and difficulty in compensating for phase errors when imaging large squints, curve trajectories, and complex scenes. Furthermore, traditional temporal BP / FBP algorithms require calculating slant range and phase compensation pixel-by-pixel and azimuth-by-time, resulting in an exponential increase in computational load, making it difficult to meet the rapid processing requirements of massive echo data from spaceborne curve SAR. While graphics processing units (GPUs) possess powerful parallel computing capabilities and are widely used in real-time SAR imaging processing, the computing power and storage capacity of a single GPU are limited, failing to meet the demands of ultra-high bandwidth and ultra-large mapping strip curve SAR imaging processing. The method proposed in this invention first performs auxiliary data reading on the spaceborne curved SAR echo data and calculates scene matching grid division parameters to adapt to the spatial variation of the wavenumber spectrum. Then, for the echo data of each sub-aperture, it performs pulse compression processing, coarse image back projection, and two-stage spectral compression-decompression phase compensation parallel processing to achieve high-precision image focusing. Finally, it performs sub-aperture data fusion, geometric offset correction, and geographic information two-dimensional interpolation projection parallel processing on the image slices to form image slice sub-block data. This solves the problems of low imaging processing efficiency and insufficient accuracy of time-domain algorithms in curved scenes in traditional methods. The geometric model diagram of the fast time-domain imaging processing algorithm for spaceborne curved SAR described in this invention is shown below. Figure 1 As shown, the specific algorithm principle and process are as follows: Figure 2 As shown.

[0013] The imaging processing method of this invention adopts a "CPU + multi-GPU" parallel processing architecture, using a CPU multi-threaded + CUDA hybrid programming mode. The CPU is mainly responsible for task scheduling, data calculation and distribution, and result storage, while each GPU is responsible for parallel imaging calculation and sub-block data fusion. The overall multi-GPU curve SAR imaging parallel architecture is shown in the figure below. Figure 3 As shown, it includes the following steps: Step 1: Utilize spaceborne SAR scene matching curve imaging echo data to obtain a fast temporal imaging processing grid partitioning sequence: The process of mesh generation along the curve is as follows: Figure 4 As shown, this includes beam illumination analysis, non-uniform sub-aperture meshing, and two-level mesh offset calculation. The specific steps for meshing along the curve zone are as follows: Step 1.1: Extract the spatial variation characteristics of the target wavenumber spectrum on the curved imaging strait using feature points. Establish a Cartesian coordinate system at the center of the curved imaging strait, where the Z-axis is perpendicular to the Earth's surface, the Y-axis is set as the reference direction, and the X-axis is set as the extension direction, with the X-axis direction perpendicular to the Y-axis. Distribute N sufficiently dense feature points evenly across the curved imaging strait. At each azimuth time instant... The components of the wavenumber spectrum decomposed along the X-axis and the component along the Y-axis It can be represented as (1) in With the center wave number, This indicates that c is the speed of light. This is the operating frequency of the radar. and These are the coordinates of the feature target point along the X-axis and Y-axis, respectively. For the satellite and the target point at time The slope distance.

[0014] Step 1.2: Extract the two-level grid offsets of the coarse and fine images from the curve grid after sub-aperture division; Step 1.3: Store the two-level grid offsets of the coarse and fine images calculated in Step 1.2 as a grid partitioning sequence parameter structure. The recorded parameters will be used in the subsequent coarse and fine image imaging processing flow.

[0015] Step 2: Perform pulse compression and back projection parallel processing on the echo data of each sub-aperture under each GPU to obtain the coarse image under that sub-aperture: The principle and process of coarse image imaging processing are as follows: Figure 5 As shown, it includes three steps: grid point beam determination, pulse compression, and parallel processing of back projection.

[0016] Step 2.1: Based on the two-dimensional beamwidth of the radar beam, establish the equation of the elliptical cone surface, and determine the target illumination by the beam at the imaging grid points, as follows: If the Z-value of the imaging grid point in the transformed coordinate system is greater than the Z-value of the X-coordinate and Y-coordinate corresponding to the elliptical cone surface, the target is considered to be inside the elliptical cone surface and within the beam illumination range; otherwise, it is outside the beam illumination range. This process obtains information on whether the imaging grid point is illuminated by the beam at each time step, resulting in the calculated ground beam ellipse template. Here, the two-dimensional kernel function length is divided according to the azimuth direction Na and the range direction Nr. The grid division parameters, the positions of the satellite and wave foot in the scene system, and the X-axis / Y-axis / Z-axis direction vector data of the beam are copied to the video memory. The grid point beam judgment kernel function is then calculated, outputting a fast temporal imaging effective grid point matrix. Step 2.2, pulse compression matched filtering, achieves high resolution in the range direction and suppresses range sidelobes; this is a core step in temporal imaging preprocessing. Under a single GPU, pulse compression is achieved using fast frequency-domain convolution. The computation flow is as follows: (2) (3) In the formula: This is the original echo range-to-FFT result. This is a range-oriented matched filter function. This is the inverse Fourier transform of the distance direction.

[0017] Here, the length is divided according to the azimuth direction Na and the range direction Nr of each sub-aperture. The echo data is read, and the memory of size Na*Nr is allocated on the host to store the range-compressed data. The CufftPlanMany in the CUFFT library is used to formulate a batch FFT plan, and the cufftExecZ2Z function is called to complete the one-dimensional FFT of the echo data in batches by segmenting the Na*Nr length data. Step 2.3: Write a kernel function to perform a dot product of the FFT data and the matched filter coefficients, padding with zeros in the middle to complete the Sinc interpolation, and store the result in the allocated video memory.

[0018] Step 2.4, similar to step 2.2, calls the CufftExeZ2Z function again to perform IFFT in parallel batches on the data after video memory interpolation in step 6, thus completing distance compression.

[0019] Step 2.5 parallelizes the projection process corresponding to each grid point, with each grid point imaged independently. First, beam projection is performed on each grid point to determine its corresponding echo, and then the echo projections are accumulated to obtain the image points. After processing all grid points, a complete coarse image is obtained. The specific back projection process is as follows: Figure 3 As shown.

[0020] The specific steps of back projection are as follows: 1. The back projection operation is performed sequentially along the azimuth direction, with azimuth time... The pulse pressure data and necessary parameters are copied from the CPU to the GPU and the backprojection kernel function is executed.

[0021] 2. At the same azimuth time, for an imaging grid of size Na*Nr, Na*Nr parallel threads are allocated to the kernel function. Each thread in the kernel function independently completes the back projection operation of a point in the grid, including the calculation of the slant range between each grid point and the radar satellite, and projects the signal corresponding to the slant range delay of the curve imaging grid point onto the corresponding grid point.

[0022] 3. Simultaneously compensate for the Doppler phase of the imaging grid points corresponding to this azimuth and time. The pulse compression phase compensation coefficient is: (4) In the formula, wavelength , This represents the distance from the target point to the satellite.

[0023] 4. Finally, the contributions of the pulses at each azimuth time point are coherently accumulated to obtain the reconstruction results of each pixel in the imaging grid. This back projection calculation is extremely computationally intensive, so a two-stage spectral compression-decompression is introduced to achieve data dimensionality reduction, and GPU parallel computing is combined to improve processing efficiency.

[0024] Step 3: Perform a two-stage spectral compression-decompression full-process imaging parallel pipeline processing on the coarse image to obtain fine images for each sub-aperture: The flowchart of the two-stage spectrum compression-decompression algorithm is as follows: Figure 6 As shown, it specifically includes four core modules: first-level spectral compression, second-level spectral compression, second-level spectral decompression, and first-level spectral decompression. Each module embeds corresponding phase compensation coefficients to achieve high-precision image focusing processing. The operation steps are as follows: Step 3.1, First-stage spectrum compression filter Represented as (5) in, These are the grid coordinates. Is it a satellite? The position at any given moment.

[0025] Then, apply an additional function. To correct frequency distortion, the function is expressed as: (6) in It is the wavenumber along the reference direction Y-axis after the first stage of compression. This is the offset of the mesh relative to the center position of the scene matching mesh on the Y-axis. By multiplying it with an additional function, the orthogonality between different meshes can be equivalently restored, thus converting the spectrum into a compressible state.

[0026] The first stage of spectral compression performs frequency domain downsampling in the range direction, eliminating redundant frequencies, reducing data dimensionality, and compensating for linear motion phase errors in the range. First, a range-direction FFT transform is performed on the backward-projected data to perform linear motion phase compensation: (7) The length can be divided according to the azimuth (Na) and range (Nr) of each sub-aperture. The coarse image data is read by back-projection of the backplane (BP). The CufftPlanMany function in the CUFFT library is used to formulate a batch column-wise FFT plan. The cufftExecZ2Z function is called to perform a batch one-dimensional range FFT of the coarse image data of length Na*Nr, and then compared with the first-stage spectral compression filter. Multiplying in the kernel function completes the first stage of image spectral compression.

[0027] Step 3.2, the second-stage filter considers the two cases expressed in equation (8). If the upsampling direction is the same as the expansion direction, then the expression for the second-stage filter is: Otherwise, the expression for the second-stage filter is:

[0028] (8) in It is the wave number along the expansion direction after the first stage of compression. and These are the coordinates of the center of the imaging grid in the reference direction and the extended direction, respectively.

[0029] The second-level spectral compression performs dimensional reduction in the azimuth direction and, combined with the characteristics of the nonlinear trajectory, compensates for the nonlinear phase error in the azimuth direction. Specifically, it performs an FFT transform on the azimuth direction to convert the time domain to the frequency domain, and then superimposes the nonlinear trajectory phase compensation. (9) The length can be divided according to the azimuth (Na) and range (Nr) of each sub-aperture. The first-level image spectral compression data is read, and the CufftPlanMany function in the CUFFT library is used to formulate a batch, row-by-row FFT plan. The cufftExecZ2Z function is called to perform a one-dimensional azimuth FFT on the Na*Nr length data segments in batches, and then combined with additional functions. The coefficients of the second-stage filter are multiplied in the kernel function to complete the second-stage image spectral compression.

[0030] Step 3.3: After spectral compression, upsampling is performed in the two-dimensional wavenumber domain using zero padding. Finally, the compressed spectrum is recovered by compensating for the conjugate of the spectral compression filter.

[0031] Step 3.4, the second-level spectrum decompression is the inverse operation of the second-level spectrum compression, recovering the azimuth time-domain signal through azimuth-to-IFFT, and simultaneously compensating the phase compensation coefficients of the second-level spectrum decompression: (10) Step 3.5: The first-level spectral decompression is the inverse operation of the first-level compression. The full-resolution range-time domain signal is recovered by IFFT from the range, compensating for the range-space-variable phase error. (11) Using the same parallel processing method, the range Nr and azimuth Na sampling points are mapped to GPU thread blocks and threads. Each thread is responsible for the azimuth row IFFT, the range column IFFT, the phase coefficient calculation, and the multiplication operation of a single sampling point, realizing the second-level and first-level spectral decompression, and finally outputting the focused fine image data with full phase compensation.

[0032] Step 4: Utilize multi-GPU parallel processing to perform data sub-block fusion and geographic coordinate interpolation projection on the fine images of each sub-aperture, ultimately obtaining the SAR sub-image: The specific steps are as follows: Step 4.1: Calculate the current image segment information based on the current number of GPUs, store it as an image segment parameter structure, and record the parameters for use in the subsequent sub-block data fusion and geographic information interpolation projection algorithm process.

[0033] Step 4.2: Start each GPU thread, load the sub-aperture data required for the current image segment, put the data into the sub-block image according to the distance offset of the fine mesh, and generate the corresponding image file after the sub-block data is fused.

[0034] Step 4.3: Load the sub-block data fusion data, perform geometric offset correction, bilinear interpolation, and geographic coordinate projection on the image tile data, and finally generate a geographic information interpolated projection image file.

[0035] Experimental Verification: The echo data parameters of "Luojia-2 01" satellite, which were verified by the multi-GPU satellite-borne curve SAR fast time-domain imaging processing experiment, are shown in Table 1.

[0036] Table 1. List of experimental parameters for fast temporal imaging processing of multi-GPU satellite-borne curve SAR

[0037] First, to verify the advantages of the multi-GPU satellite-borne curve SAR fast temporal imaging processing method in solving the problem of high-resolution SAR imaging along curve scenes with beam maneuver constraints, a set of echo data was processed temporally using the traditional back projection algorithm (i.e., the method of directly performing pulse compression and back projection on the echo data) under the parameters in Table 1. The imaging results were compared as follows: Figure 7 As shown, the method proposed in this paper has higher focusing accuracy. Next, we use the Luojia echo data to perform fast time-domain imaging processing using pure CPU, single GPU, and multi-GPU methods, and statistically analyze the time required for different methods. Figure 8(a) shows the results of pure CPU-based fast temporal imaging. Figure 8 (b) shows the fast temporal imaging results using a single GPU. Figure 8 (c) shows the results of fast temporal imaging using multiple GPUs. The results demonstrate that all three methods can achieve sub-aperture imaging processing, sub-block data fusion, and geographic information interpolation projection, effectively ensuring imaging accuracy. However, this algorithm is more efficient in terms of processing speed, achieving both imaging accuracy and significantly reduced processing time. Furthermore, to verify the sub-aperture imaging processing time and speedup ratio under different numbers of GPUs, specific statistics are shown in Table 2. As the number of GPUs increases, the processing time decreases significantly, with a speedup ratio of 15.125 for 4 GPUs, approaching linear acceleration. This indicates that the proposed multi-GPU parallel architecture has reasonable task partitioning, low communication overhead, and excellent parallel efficiency.

[0038] Table 1. List of experimental parameters for fast temporal imaging processing of multi-GPU satellite-borne curve SAR

[0039] In summary, the above are merely preferred embodiments of the present invention and are not intended to limit the scope of protection of the present invention. Any modifications, equivalent substitutions, improvements, etc., made within the spirit and principles of the present invention should be included within the scope of protection of the present invention.

Claims

1. A fast temporal processing method for spaceborne SAR scene matching curve imaging using multiple GPUs, characterized in that, Includes the following steps: Step 1: Divide the grid along the curved zone to generate a curved SAR terrain matching gridded sequence. This includes extracting the spatial variation features of the target wavenumber spectrum on the curved imaging zone using feature points, extracting the two-level grid offsets of the coarse and fine images, and storing them as a grid division sequence parameter structure. Step 2: Perform grid point beam determination, pulse compression and back projection in parallel to form a coarse image. This includes establishing the elliptical cone equation based on the two-dimensional beamwidth of the radar beam, determining the beam illumination target for the imaging grid points, using fast frequency domain convolution to achieve pulse compression, and parallelizing the projection process corresponding to the grid points so that each grid point is imaged independently. Step 3: Perform two-stage spectral compression-decompression phase compensation to achieve high-precision image focusing and form a fine image. This includes four core modules: first-stage spectral compression, second-stage spectral compression, second-stage spectral decompression, and first-stage spectral decompression. Each module is embedded with corresponding phase compensation coefficients. Step 4: Perform data sub-block fusion processing and geographic coordinate interpolation projection to generate sub-block image data with geographic coordinate information registration, thus obtaining SAR sub-images.

2. The method according to claim 1, characterized in that, Step 1 specifically includes: (1) Establish a Cartesian coordinate system at the center of the curved imaging zone, with the Z-axis perpendicular to the Earth's surface, the Y-axis set as the reference direction, and the X-axis set as the extension direction and perpendicular to the Y-axis. Distribute N feature points evenly on the curved imaging zone. (2) Use feature points to extract the spatial variation features of the target wavenumber spectrum on the curved imaging zone, the components of the wavenumber spectrum center along the X-axis and the components along the Y-axis; (3) Extract the two-level grid offsets of the coarse image and the fine image, and store them as a grid partitioning sequence parameter structure.

3. The method according to claim 1, characterized in that, The pulse compression is achieved using frequency domain fast convolution, including: performing a fast Fourier transform on the echo data using FFT; performing a dot product of the FFT data and the matched filter coefficients using a kernel function and padding with zeros in the middle to complete Sinc interpolation; and calling the IFFT function to complete range compression.

4. The method according to claim 1, characterized in that, The parallel processing of the back projection includes: back projection operations are performed sequentially along the azimuth direction; at the same azimuth time, Na·Nr parallel threads are allocated to the imaging grid, and each thread independently completes the slant range calculation and signal projection of the grid points; simultaneously, the Doppler phase of the imaging grid points is compensated, and the pulse compression phase compensation coefficient is: ； In the formula, wavelength , The distance from the target point to the satellite is given; finally, the pulse contribution at each azimuth time is coherently accumulated to obtain the reconstruction result of each pixel in the imaging grid.

5. The method according to claim 1, characterized in that, The first-level spectrum compression includes: (1) First-stage spectrum compression filter Represented as ； in, These are the grid coordinates. Is it a satellite? The position at that moment; (2) Using additional functions To correct frequency distortion, the function is expressed as: ； in It is the wavenumber along the reference direction Y-axis after the first stage of compression. It is the offset of the mesh relative to the center position of the scene matching mesh on the Y-axis; (3) Perform range FFT transformation on the back projection data and perform linear dynamic phase compensation; the first-level spectrum compression performs frequency domain downsampling on the range direction, eliminates redundant frequency points, and reduces the data dimension.

6. The method according to claim 1, characterized in that, The second-level spectrum compression includes: (1) The expression of the second-stage filter is determined based on the relationship between the upsampling direction and the expansion direction, including two-dimensional processing of the spectrum; (2) Dimension reduction is performed on the azimuth direction. Combined with the nonlinear trajectory characteristics of the curve, the nonlinear phase error of the azimuth direction is compensated. The specific method is to perform FFT transformation on the azimuth direction to realize the conversion from the time domain to the frequency domain of the azimuth direction, and then superimpose the nonlinear trajectory phase compensation. (3) After spectral compression, upsampling is performed in the two-dimensional wavenumber domain by zero filling.

7. The method according to claim 1, characterized in that, The two-stage spectral decompression is the inverse operation of the two-stage spectral compression: (1) The second-stage spectral decompression recovers the azimuth time-domain signal through azimuth-to-IFFT, and simultaneously compensates for the phase compensation coefficients of the second-stage spectral decompression: ； (2) The first-stage spectral decompression recovers the full-resolution range-time domain signal via range-to-IFFT, compensating for range-space-variable phase errors: ； (3) The compressed spectrum is recovered by compensating the conjugate of the spectrum compression filter, and the focused fine image data with full phase compensation is output.

8. The method according to claim 1, characterized in that, The data sub-block fusion process includes: calculating image fragment information based on the current number of GPUs and storing it as an image fragment parameter structure; starting each GPU thread, loading sub-aperture data, and placing the data into the sub-block image according to the distance offset of the fine mesh; and generating the corresponding image file after the sub-block data is fused.

9. The method according to claim 1, characterized in that, The geographic coordinate interpolation projection includes: geometric correction based on grid division sequence parameters; geocoding of SAR images; and image interpolation mapping based on geographic coordinate information to generate sub-block image data registered with geographic coordinate information.

10. A software system for implementing the method according to any one of claims 1 to 9, characterized in that, It adopts a "CPU + multi-GPU" parallel processing architecture, including: a CPU main control module for multi-threaded parallel task scheduling and data computation and distribution; multiple GPU processing modules, each GPU corresponding to a sub-aperture imaging processing; and a CPU multi-threaded + CUDA hybrid programming mode, with each GPU responsible for parallel imaging computation and sub-block data fusion.

Images