基于GPU的长轨SAR实时成像算法.pdf

ISSN 10049037，CODEN SCYCE4Journal of Data Acquisition and Processing Vol.38，No.6，Nov.2023，pp.1380-1391DOI：10.16337/j.10049037.2023.06.013 2023 by Journal of Data Acquisition and Processinghttp：/Email：sjcj Tel/Fax：[phone]基于 GPU的长轨 SAR实时成像算法谭运馨1，黄海风2，赖涛2，但琪洪2，欧鹏飞2（1.中山大学系统科学与工程学院，广州 510006；2.中山大学电子与通信工程学院，深圳 518107）摘要：为了满足长轨道超高分辨 W 波段合成孔径雷达（Synthetic aperture radar，SAR）的快速成像需求，本文提出了一种基于图形处理器（Graphics processing unit，GPU）的 K实时成像算法，该算法采用并行架构和双流多线程的处理方式。默认流沿着物理原理的方向进行数据处理，首先对距离补偿、误差校正和补零等操作进行并行化处理，然后采用一层嵌套的插值方式，通过维持上下层的依赖关系和同步管理就能达到约 30的加速比。阻塞流与默认流同时启动，生成默认流所需的参数和函数，并在执行前将其存入显存，极大地缩小了算法的运行时间，同时通过在默认流上设置事件以保持双流的同步并行执行。实验结果表明，算法总的加速比可达 13左右，幅值和相位相对误差趋近 0，不仅具有良好的实时性、聚焦性，还保持了良好的成像效果。关键词：长轨道合成孔径雷达；图形处理器；实时成像；Stolt插值；K算法中图分类号：TN957.52 文献标志码：AGPUBased RealTime Imaging Algorithm for LongTrack SARTAN Yunxin1，HUANG Haifeng2，LAI Tao2，DAN Qihong2，OU Pengfei2(1.School of Systems Science and Engineering,Sun YatSen University,Guangzhou 510006,China;2.School of Electronics and Communcation Engineering,Sun YatSen University,Shenzhen 518107,China)Abstract：To meet the fast imaging requirements of long-orbit ultra-high resolution W-band synthetic aperture radar（SAR），this paper proposes a graphics processing unit（GPU）-based K real-time imaging algorithm which adopts parallel architecture and double stream multithreading processing.The default stream processes data along the direction of the physical principle.Firstly，it parallelizes the rang compensation，error correction，zero filling and other operations，and then adopts one-layer nested interpolation method.By maintaining the upper and lower dependencies and synchronization management，it can achieve a speed ratio of about 30.The blocking stream starts at the same time as the default stream，generates the parameters and functions required by the default stream，and stores them into video memory before execution，which can greatly reduce the running time of the algorithm.Meanwhile，by setting events on the default stream，the two streams can be executed synchronously in parallel.Experimental results show that the total acceleration ratio of the algorithm can reach about 13，and the relative errors of amplitude and phase are close to 0，which not only has good real-time performance and focusing performance，but also maintains good imaging effect.基 金 项 目：广 东 省 重 点 领 域 研 发 计 划 项 目(2019B111101001)；国 家 自 然 科 学 基 金(62071499)；深 圳 市 科 技 计 划 项 目(JCYJ20190807153416984)。收稿日期：20221110；修订日期：20230220谭运馨 等：基于 GPU 的长轨 SAR实时成像算法Key words:long-track synthetic aperture radar(SAR);graphics processing unit(GPU);real-time imaging;Stolt interpolation;-K algorithm引 言长轨道高分辨率成像雷达可用于机场异物检测、形变监测和科学研究等领域12。本文在前期开发完成的 W 波段厘米级分辨率长轨道合成孔径雷达（Synthetic aperture radar，SAR）和子孔径成像算法并行化的基础上，进一步探讨了基于图形处理器（Graphics processing unit，GPU）的实时成像方法。快速实时成像算法是一个热门的研究领域，学者们的研究目标主要是传统机载或者星载体制的雷达。常见的成像算法有距离多普勒算法（RangeDoppler algorithm，RDA）、CSA（Chirp scaling algorithm）算法、KA 等，RDA 和 CSA 具有高效、成熟且易于实现等特点，但这两种算法由于采用了近似泰勒展开，因此难以适应大带宽、大波束角的成像场景。KA 在二维频域通过 Stolt插值来矫正了距离方位耦合与距离向时间和方位向频率的依赖关系，避免了 RDA 和 CSA 的缺陷。因此适用于轨道 SAR 实时成像系统34。K 算法处理雷达图像需要较高的内存和处理时间。Sahay等5提出了一种改进时域参考函数的方法，缓解了 SAR 图像生成处理时间长和高内存的问题。参考函数的生成及运算要求更高的计算量，所以快速傅里叶变换（Fast Fourier transform，FFT）及快速傅里叶逆变换（Inverse fast Fourier transform，IFFT）需要更快的转换时间6。所以仅对参考函数进行改进，对整个系统的实时性来说远远不够。文献 7 在现场可编程逻辑门阵列（Field programmable gate array，FPGA）上实现了对 K 算法的优化，其对 K 算法中的 Stolt插值部分采用数据循环存储，然后进行插值计算。与文献 5 一样，他们仅对 K 算法中的插值部分进行优化，并没有对整个算法进行优化处理。当处理的数据量变大，分辨率变高，不只是插值部分的工作量变大，整个 K 算法的复杂度也会变大，同时采用的 FPGA 平台的计算能力和并行化程度并不高。Zhu 等8在 FPGA 上实现了极坐标格式算法（Polar format algorithm，PFA）实时处理，该方法相比 CPU 处理速度有所提升，但是依然无法满足实时成像的要求，并且 FPGA时序难以规划、程序调试周期长、算法实现复杂。Zou等9在一种数学框架上提出了一种高效准确地基合成孔径雷达（Groundbased synthetic aperture radar，GBSAR）算法，避免了庞大的计算量；但工程实现难度较大，应用于实际工程当中也未必能达到理想的效果。所以从硬件上进行加速方面的研究成了最佳选择。GPU 强大的处理能力、高度并行性及高显存带宽为 K 算法提供了具有发展前景的新型运算平台10。因此众多学者在 GPU 平台上对 SAR 成像算法进行并行化处理1113。本文提出了一种基于GPU 平台的实时成像算法。与其他并行算法不同，本算法依据各个部分的独立性和并行性，或对算法进行改进使各个部分之间体现最大限度的独立性和并行性以进行整体加速。其中对插值部分提出了 3种优化方案并进行对比分析，其中一层嵌套很好地适配了实时成像系统，将大数据量多重循环的插值计算优化成更细粒度的并发运算。利用双流方式将相关参数的生成与相关函数的执行并行进行，避免了设备端与主机端频繁的内存访问，将参数和函数的执行隐藏于默认流的执行当中，极大地缩减了实时成像算法的运行时间。1 K并行化优化处理策略 在速度恒定的条件下，K 算法具有在宽孔径或大斜视角范围内校正沿距离向的距离徙动变化的1381数据采集与处理 Journal of Data Acquisition and Processing Vol.38,No.6,2023能力14。本文提出了在 GPU 端实现的 K 实时成像算法，如图 1所示。图 2为本算法雷达系统工作示意图，表 1为该雷达系统的参数设置。1.1基于非阻塞流的参数生成在 统 一 计 算 设 备 架 构（Compute unified device architecture，CUDA）中，同一流中的操作有严格的执行顺序，而在不同流中的执行顺序不受限制15。如图 1 所示，创建一条与默认流相并行的非阻塞流。在相关函数执行之前，在非阻塞流上提前将所需函数（refer 函数等）、变量（Kr，kx）生成并保存至显存中。图 3所示为采用默认流、阻塞流和非阻塞流执行生成 refer函数的时间坐标。图 3中 3个条形柱的宽度为 refer函数生成的时间，为 18.7 ms左右。由图 3 可知相比于串行的默认流，并行执行的非阻塞流要提前 1 s 运行；即使是阻塞流也要比默认流提前0.94 s；同时非阻塞流比阻塞流节省 0.072 3 s的运行时间。如图 1所示，refer函数、汉明窗 1、汉明窗2 等函数位于与默认流相并行的非阻塞流。两条流是并行执行，同时默认流上设有 4个事件，他们会阻止默认流中的相关操作的执行，直到所需函数和变量已经提前保存在显存当中。图 4为非阻塞流和 CPU 上各个生成函数的运行时间，不难发现：在 CPU 上生成汉明窗的时间成本要更少一些。因为汉明窗生成的工作量很小，以至于内核的启动时间远大于核函数的运行时间。但考虑到 CPU 与 GPU 之间数据传输要花更长的时间，汉明窗的生成不能放到 CPU 端上执行。因为远距离目标的回波比较弱，在模拟域的补偿较少，所以要在数字域进行补偿，补偿参数为距离的 1.5次方。如图 4 所示，距离补偿函数的生成，GPU 可以加速 4.18 倍。采用双流并行的方式，从 GPU 端来说，总时长节省了 61.34 ms；从 CPU 端来说，总时长可以节省 364.395 ms。图 1长轨道 SAR K实时成像算法实现Fig.1Realization of K realtime imaging algorithm for longorbit SAR图 2雷达系统工作示意图Fig.2Radar system working diagram表 1回波数据相关参数Table 1Related parameters of echo data系统参数载波频率/GHz距离向带宽/GHz脉冲重复频率/Hz采样频率 fs/MHz回波信号尺寸（NaNr）平台速度 Va/(ms-1)数值9564002.5 4 6108 0000.302 51382谭运馨 等：基于 GPU 的长轨 SAR实时成像算法1.2方位向、距离向抽取二维频域并行优化策略待处理的数据存在于 CPU 的内存中，需要将其读取并保存到 GPU 的显存当中；由于成像距离比较短，方位向和距离向过采样，所以要进行方位向和距离向抽取以降低采样率。为了进行插值前的数据处理，数据的读取、方位向和距离向抽取后需要进行 FFT。同时 K 算法需要在二维频域对信号进行处理。而在数据量较大的时候 FFT计算量非常大，通常可以利用旋转因子的周期性、对称性等性质进行化简。对 FFT/IFFT 的快速计算进行相关的研究很早就开始进行，Cooley 和Tukey16将 离 散 傅 里 叶 变 换（Discrete Fourier transform，DFT）的算法复杂度从O(NlgN)降低为O(N2)。GPU 因其计算能力强，众多学者尝试在 GPU 上进行 FFT 的快速实现1718。本文将从 CPU端和 GPU 端上分析 FFT的加速方法。Frigo 和 Johnson19开发了一种在 CPU 端使用的快速计算离散傅里叶变换的标准 C 语言程序集FFTW。FFTW 比 CPU 端其他傅里叶变换要更快、更高效，因而被广泛应用于如 Intel的数学库和 Scilib等软件上。GPU 端的 CDUA 库提供了一系列库函数来提高开发人员的开发效率。其中 CUFFT 库提供了一个优化的快速傅里叶变换，其加速效果不亚于 FFTW，并且避免了主机和设备之间频繁相互访问的时间消耗。表 3 和图 5 对 CPU 和 GPU 上两种不同的 FFT 加速库进行了对比。可以发现随着点数的增加，CPU 端 FFT/IFFT 时间也会随之增加，两者之间的加速差距也会明显变大。由图 5右图的加速比曲线可以看出，随着点数增高，CUFFT 的加速效果要优于 FFTW。因此，在 FFT 加速方面，GPU 端的加速图 3默认流、阻塞流和非阻塞流生成 refer函数的时间坐标Fig.3Time coordinate of refer function in default stream,blocking stream and nonblocking stream图 4非阻塞流运行各个参数函数的时间Fig.4Running time of each parameter generating function on nonblocking stream表 2CPU端与 GPU端加速效果对比Table 2Comparison of CPU and GPU acceleration effects点数1e+63e+65e+67e+61e+73e+78e+73e+85e+8CPU FFTW 库耗时/ms30.905 4235.752409.01629.255836.8462 736.727 451.2433 304.8124 228CUDA CUFFT 库耗时/ms17.5832.013 337.379 149.980 455.084101.207249.291224.497540.423加速比1.767.3610.9412.6015.1927.0429.89148.35229.871383数据采集与处理 Journal of Data Acquisition and Processing Vol.38,No.6,2023能力要远超 FFTW。并且若从 CPU 端进行加速，需要将待处理的数据传到 CPU，然后再将 FFT/IFFT后的结果再传回 GPU，导致多出了两次数据传输的时间消耗，不能实现处理过程中的结果暂存，且违背了 GPU 程序编程的 3条通用法则中的第 1和第 3法则20。CUDA 库需要调用设备变量，普通主机函数无法直接访问设备变量，需要先将数据从显存传回CPU 内存中。所以本文通过 CUFFT 库实现了距离向和方位向上的 FFT 和 IFFT，并将其封装成一个宏函数可供主机端调用，因为宏函数的形参可以为设备变量或者主机变量。表 3为文件读取，距离向和方位向抽取、二维频域转换的时间花销对比。可以发现，在 GPU 上进行列 FFT 的时间要比行 FFT 时间要长一些。因为 CUFFT 采用行优先的数组存储方式，这会导致设计列FFT/IFFT 宏函数时，要在 FFT 前后进行两次转置操作。行和列的逆傅里叶变换的结果需要除以点数进行矫正，故其运行时间要比正傅里叶变换的时间长。由表 3可知从 GPU 端进行加速的效果更明显，加速比可达数十倍。因为本算法在离散化之前距离向就已经完成解线性调频，所以距离向已经处在频域。故只要进行方位向傅里叶变换即可进入二维频域。表 3两种平台插值前的数据处理运行时间Table 3Data processing running time before interpolation of two platforms项目数据读取方位向抽取距离向抽取二维频域项目列 FFT列 IFFT列 FFT列 IFFT方位抽取行 FFT行 IFFT距离抽取列 FFT频谱搬移数据量4 6088 0002 3048 0008 0002 3048 0002 3044 0002 3044 0002 3044 0002 3044 0001 1524 0001 7284 0001 728GPU 耗时/ms50.95430.20529.28737.2774.3066.52612.5552.17410.7746.043CPU 耗时/ms1 221.210731.6121 096.0301 099.490113.867229.235249.79757.371204.691139.166加速比23.9724.2237.4229.5026.4435.1319.9026.3919.0023.03图 5FFTW 和 CUFFT 运行时间对比及加速比Fig.5Speed-up ratio and comparison of running time for FFTW and CUFFT1384谭运馨 等：基于 GPU 的长轨 SAR实时成像算法1.3系统误差矫正及距离频域前端补零并行优化策略1.3.1系统误差矫正天线相位中心到目标的距离为Rsys=0.08 m，该距离使得目标回波在空中和射频电缆中的传输存在时间延迟Tsys。根据线性调频中的时频关系，频谱会产生一个Fsys的频移。这是一个与系统有关的固定误差，不会随数据的变化而更新。Tsys=2 Rsys/c（1）Fsys=Kr Tsys（2）式中：Kr为距离向调频率，c为光速。想要补偿频谱中产生的相位差，需要乘以如式（3）所示的补偿参考函数Hrefer=exp(-j2fsys)（3）式中：Hrefer对应图 1中的非阻塞流上生成的 refer函数，为距离向快时间。图 1中设置的事件 1引入了默认流和非阻塞流之间的依赖关系；默认流在运行到事件 1 时，若非阻塞流还未运行到事件 1，则会阻塞默认流的继续执行直到非阻塞流完成 refer函数的生成。因此采用双流并行执行的方式，保证进行系统误差矫正前，refer参考函数已经提前完成生成并保存至 GPU 的显存之中，从而将生成参数的时间隐藏于默认流之中。如图 6所示，refer生成函数（get_t）的运行时间几乎可以完全隐藏于 CPU 到 GPU 的数据传输（Memcpy HtoD）当中。表 4所示为系统误差矫正的时长，由于补偿函数refer函数的生成是相位补偿前在非阻塞流上生成的，所以 GPU 上系统误差矫正总时长仅包含相位补偿过程的时长 6.911 ms。由表 4 可以看出，在GPU 上进行相位补偿的加速比是 212.73。整个系统误差矫正中，GPU 加速比为 246.40。1.3.2距离前端补零、方位向、距离向加窗插值前的数据是非均匀数据，无法直接与回波数据对应，所以无法直接进行逆傅里叶变换。因此需要通过 Stolt插值将非均匀的数据重采样为均匀的数据。进行 Stolt插值后的数据前端会向前突起，突起的部分超出了源数据的范围，如图 7 右图阴影部分。故在Stolt插值前，需要对源数据的前端进行补零，扩大源数据的范围以覆盖突起的部分。同时为了降低旁瓣对主瓣的影响，防止频谱信号能量泄露，提高精度，需要在方位向和距离向进行加窗聚焦。如表 5所示为 GPU 和 CPU 端加窗和补零的运行时间。在 GPU 端上，距离前端补零和距离向加窗图 6refer函数的双流处理Fig.6Operation of refer function in two streams表 4GPU和 CPU上系统误差矫正的运行时间Table 4Running time of system error correction on GPU and CPUms项目CPUGPU加速比参数生成232.71218.595(非阻塞流)12.51相位补偿1 470.186.911(默认流)212.73总时长1 702.8926.911246.40图 7Stolt 映射Fig.7Stolt mapping1385数据采集与处理 Journal of Data Acquisition and Processing Vol.38,No.6,2023同时进行。在 GPU 端进行距离向补零加窗的加速比为 118.61，方位向加窗的加速比为 163.71。由 1.1小节分析可知，在窗函数的生成上，CPU 的时间明显快于 GPU 时间。由于考虑到 CPU 和 GPU 之间数据传输的成本，不能在 CPU 端生成窗函数。同时 GPU 上窗函数生成时间可以隐藏于默认流之中，即GPU 上生成窗函数是几乎没有时间成本。同时默认流中设置的事件 2 和事件 3 确保了窗函数的提前生成。1.4基于动态并行的插值运行Stolt插值采用的插值因子通过距离向频率轴的映射和弯曲来完成21，Stlot插值完成了残余距离徙动矫正、残余二次距离压缩和残余方位压缩。由于 Stolt插值对成像质量的影响较大，所以对于 Stolt插值要求精确计算，这就导致 Stolt插值计算量大22。本算法对比了 3种 Stolt插值的并行化方法。1.4.1逐行进行插值（零层嵌套）插值本身包含两重循环，在加上回波数据的一重循环，共 3重循环。在对子孔径成像算法进行并行化处理中，采用从已知数据或者估计数据中的一个方向进行插值。已知数据为真实时域，是一维数据，而估计数据即虚拟时域为二维数据。如果从虚拟时域方向进行插值，可以同时对整个二维矩阵进行遍历、比较和插值，所以只要虚拟时域点数不超出 GPU 所限制的线程数，运行时间理论上会保持不变。在K 算法 Stolt插值中，已知数据是二维数据矩阵 Na1Nr1（4 0001 728），而估计数据是一维数据数组，且插值的比较范围是已知数据矩阵中的每一行。如果整个已知数据同时进行遍历、比较和插值，则在比较范围以外的其他行，有可能存在更接近于估计数据的临近点，导致该点进行错误的插值，进而导致整个估计数据进行错误的插值。所以 k算法中 Stolt插值部分不能进行二维数据同时插值，但是可逐行进行插值。为了与一层动态嵌套和二层动态嵌套进行对比解释，本文将逐行进行插值称为零层动态嵌套。如图 8所示，零层嵌套对二维频域数据上的单个点进行插值运算，一个估计数据与 Nr1个已知数据同时进行比较，筛选出临近点，进行线性拟合。二维数据总共有 Na1Nr1个点，故总共需要进行 Na1Nr1次点插值运算。1.4.2动态并行进行插值动态并行允许 GPU 端在已有的 GPU 内核中直接创建和同步新的 GPU 内核15。实现核函数中嵌套其他核函数。插值部分的工作量在串行条件下以 3重循环进行点对点的遍历、比较和插值。所以这就限制了其本身最多进行深度为 3的嵌套。表 5距离前端补零、方位向和距离向加窗对比结果Table 5Comparison results of zero padding at the distance front end，azimuth and range windowing平台CPUGPU项目距离前端补零距离向加窗方位向加窗距离向补零加窗加速比方位向加窗加速比结果82.217 2 ms261.874 ms388.3262.901 ms118.612.372 ms163.711386谭运馨 等：基于 GPU 的长轨 SAR实时成像算法（1）一层嵌套动态并行图 9为深度为 2的嵌套示意图，父核函数为 nestinterp1，子核函数为 interp2_GPU。其中子核函数进行每一行的插值运算，父核函数在 Na1个线程上即点处同时启动 Na1个子核函数，同时完成 Na1行插值运算。当运行至处时，子核函数已经完成每一行的插值运算。本实验平台为 NVIDIA GeForce GTX 1060 6 GB，其多流处理器（Streaming multiprocessor，SM）的数量为 10，流处理器（Streaming processor，SP）数量为 128。估计数据的大小为 Na1Nr1。CUDA 程序中要有足够的并行性，块的数量要多于 SM 的数量。父核函数和子核函数均采用一维网格。综上所述，父核函数和子核函数的线程块的中线程应满足(4 000+blocksize1-1)blocksize1 10（4）(1 728+blocksize2-1)blocksize2 10（5）可以求出设备并行性的临界值为blocksize1 444，blocksize2 192。线程块尺寸越小于 444 和192，单位 SM 内的线程块的数量就越多。同时线程块的最大数量也受到硬件资源的限制，所以高占用率并不意味着更高的性能，其他因素同样也限制 GPU 的性能。为了避免线程块太小，每个块至少包含128个线程。为了能极大地提高加载效率，同时要保持每一个线程块中线程数量是线程束大小（32个线程）的倍数。所以从blocksize1 444，blocksize2 192的配置附近寻找最佳执行配置，如表 6所示。图 8零层嵌套插值Fig.8Zero-layer nested nested interpolation图 9一层嵌套工作原理图Fig.9One-layer nested working principle diagram1387数据采集与处理 Journal of Data Acquisition and Processing Vol.38,No.6,2023由表 6可以看出，最佳执行配置是 blocksize1=448，blocksize=128。在 blocksize1=352，blocksize=128 的配置中，线程块数量最多，并行性也是最好的；但其占有率（85.94%和 75%）却不是最高的。同理，最佳配置 blocksize1=448，blocksize=128 并行性不是最好的，占有率（87.50%和 75%）也不是最高的。占有率最高的执行配置 blocksize1=384，blocksize=128却不是最佳配置。所以没有一个单独的指标能直接优化性能，需要从多个相关指标和 GPU 硬件资源之间寻找一个恰当的平衡来实现最佳总体性能。占用率是每一个 SM 中活跃的线程束数量和最大线程束数量之间的比值，即占用率=活跃线程束数量最大线程束数量（6）NVIDIA GeForce GTX 1060 6 GB 上的每一个 SM 上的最大线程数量为 2 048，所以最大线程束数量为 2 048/32=64。由此可以算出最佳配置的父核函数和子核函数活跃线程束数量分别为父核活跃线程束数量=64 0.875=56（7）子核活跃线程束数量=64 0.75=48（8）（2）双层嵌套动态并行由（1）可知，子核函数的最佳线程块大小为 128，父核最佳线程块大小为 448。父网格 interp_lay3即最外层循环的大小为 Na1，中间循环 interp_lay2 即第 1 层嵌套和内层循环 interp_lay1 即第 2 层嵌套的大小为 Nr1。所以第 1、2层嵌套的块大小设置为 128，父网格设置为 448。父网格调用第 1层嵌套，然后第 1层嵌套调用第 2层嵌套，如图 10所示，所以可以进行深度为 3的嵌套。表 6执行配置运行时间Table 6Time of execution configuration point operationms父核子核 128(75%)子核 160(70.31%)子核 192(75%)352（85.94%）1 580.8371 998.4872 357.737384（93.75%）969.2571 588.9572 068.197416（81.25%）1 219.7071 042.0571 785.447448（87.50%）770.417804.967919.517480（93.75%）1 648.5371 778.7471 919.007图 102层嵌套工作原理图Fig.10Twolayer nested working principle diagram1388谭运馨 等：基于 GPU 的长轨 SAR实时成像算法由表 7 可知，0 层嵌套和 1 层嵌套都具有加速效果，而 2层嵌套运行时间反而比 CPU 端的更慢。实际上，实时成像算法需要全天候、全天时对地质体进行连续监测，由此产生的数据量是巨大的。在新的一级网格中启动另一个网格需要额外消耗内存资源和时间花销，同时每一个上层嵌套和下层嵌套之间的同步管理又要耗费大量的设备内存和时间。在本算法的 2层嵌套中，在父网格的所有线程中开辟一个新的网格，新网格的所有线程又启动另一个新的网格，需要额外开辟新的内存资源来启动新的网格和保持新旧网格之间的同步，同时嵌套深度越深，嵌套的各层之间的同步需求就更加频繁，需要维持的依赖关系也会变得更复杂。如图 9、10所示，1层嵌套之间的同步管理只有上、下两层之间，只需维持上层和下层的依赖关系；而 2层嵌套需要保持上中、上下、中下、上中下的同步管理和依赖关系，理论上同步管理的工作量增加了 4倍以上，需要维持的依赖关系也增加到了 4个。由表 7可知 2层以上嵌套的加速效果已经不足以抵消各层之间同步和维持各个层之间的依赖关系所带来的时间消耗。所以对于实时成像系统来说，1层嵌套的加速效果比零层嵌套和双层嵌套的加速效果要好得多。加速比可达 30倍左右。2 结果分析 表 8 为 CPU 端和 GPU 实时成像算法的运行时间，运行大小为 4 6108 000的实测数据，其中CPU 端的运行时间为 30.954 s，GPU 的实时成像算法运行时间为 2.440 s，加速比约 13 倍。图 11是 CPU 与 GPU 处理结果的幅值和相位相对误差帕累托图。不难发现 CPU 端和 GPU 端处理结果的相位和幅值相对误差非常小，其中幅值相对误差有 63.18%的点分布在(10-13，10-12的范围内，幅值相对误差分布在(10-14，10-11的点占 99.28%。同理相位相对误差分布在(10-14，10-12范围内的点数占了 85.99%，而分布在(10-14，10-11的点占了 98.43%。因此 GPU 端的实时成像算法相位和幅值相对误差非常小，几乎为 0。该误差主要来源于 GPU 端和 CPU 端舍入规则和数据截断方式的不同。以上分析说明：为了最大限度地体现 GPU 算法的并行性和独立性所作的并行化处理并没有破坏算法的物理原表 7采用 0、1、2层嵌套的算法总用时Table 7Total time of using zerolayer，onelayer and twolayer nested algorithms项目运行时间/s加速比0层嵌套13.0381.761层嵌套0.77029.782层嵌套319.7460.072CPU 插值22.931表 8CPU端与 GPU端实时成像算法运行时间Table 8CPU and GPU realtime imaging algorithm running timeCPU 端/s30.954GPU 端/s2.440加速比12.69图 11相位和幅值相对误差分布Fig.11Relative error distribution of amplitude and phase1389数据采集与处理 Journal of Data Acquisition and Processing Vol.38,No.6,2023理，同时本算法的成像效果不会因为追求实时性而有所降低。图 12为 CPU 成像算法和 GPU 成像算法的成像效果显示。不难发现，GPU 实时成像算法具有良好的聚焦性、实时性和良好的成像效果。3 结束语 本文提出了一种基于 GPU 的轨道 SAR K 实时成像算法。本文对距离补偿、误差校正和补零等操作进行了优化处理。通过双流的并行执行方式，将系统误差、方位向和距离向加窗、插值和距离补偿所需的参考函数和窗函数放在与默认流相并行的流中执行。将 CPU 端需要 364.395 ms完成的时间花销隐藏于默认流的各个操作和数据传输当中，大大缩小了系统误差矫正等操作的运行时间，同时避免了主机和设备之间的交互。通过 CUFFT 设计了行、列 FFT 和 IFFT 的宏函数，实验证明该宏函数比CPU 端的加速库快数 10倍。本算法提出了 0层、1层和 2层嵌套的插值方案，经过分析对比，发现零层和单层嵌套的插值方案均有加速效果。其中单层的嵌套加速比可达 30倍左右。因为层数的增加，同步管理的需求变多以及依赖关系变得更为复杂，导致双层加速效果低于 CPU 端，所以双层的插值方式并不适合于实时成像系统。基于大小 4 6108 000 的实测数据表明，本算法可以加速 13 倍左右，并且幅值和相位相对误差为百亿分之一，趋近于零。同时本算法具有良好的实时性、聚焦性和成像效果。参考文献：1SERRANO-JUAN A,VZQUEZ-SU E,MONSERRAT O,et al.Gb-SAR interferometry displacement measurements during dewatering in construction works:Case of La Sagrera railway station in Barcelona,SpainJ.Engineering Geology,2016,205:104-115.2LIU Xiangzeng,ZHENG Tian,QIANG Lu,et al.A new affine invariant descriptor framework in shearlets domain for SAR image multiscale registrationJ.AEU-International Journal of Electronics and Communications,2013,67(9):743-753.3ZHANG Lei,SHENG Jialian,XING Mengdao,et al.Wavenumber-domain autofocusing for highly squinted UAV SAR imageryJ.IEEE Sensors Journal,2011,12(5):1574-1588.4XU Gang,XING Mengdao,ZHANG Lei,e