-
合成孔径声呐(synthetic aperture sonar, SAS)是一种高分辨率水下成像声呐, 其原理是利用小尺寸基阵沿航迹方向匀速直线运动来虚拟大孔径基阵, 借助信号处理手段实现回波信号相干处理, 获取与成像距离和声波波长无关的高分辨率二维声呐图像[1-3]。因其方位向分辨率比常规成像声呐高1~2个数量级, SAS在水底地形地貌测绘[4]、水下小目标探测[5](如猎雷、考古等)、水底地质鉴别与分类、海洋石油管线巡查等方面具有广阔应用前景。
SAS技术的关键在于如何根据多接收阵回波数据快速完成成像场景的二维图像重构。目前SAS成像算法主要分为时域算法[6-7]和频域算法[8-10]两大类。在合成孔径声呐发展早期, 因信号处理机性能不足, 时域成像算法难以满足计算速度需求, 导致了相关频域算法的出现。但频域成像算法都是以各种近似为基础, 通过降低成像质量来达到快速成像的目的。常用频域成像算法有距离多普勒算法[9]、线频调变标算法[10]、波数域算法[11]和线性调频Z变换算法[12], 它们的共同缺点是难以与运动补偿相结合, 成为限制成像质量进一步提升的瓶颈。时域成像算法中的后向投影算法是一种精确逐点成像算法, 其特点是具有精确聚焦能力, 适合宽波束大场景成像, 并且易于进行运动补偿。但早期SAS后向投影成像快速算法主要利用了距离徙动曲线的对称性, 借助集群计算达到快速计算目的, 存在信号处理设备体积大, 难以进行运动补偿的缺点[13]。近年来, 随着计算机技术的发展, 特别是图形处理器(graphics processing unit, GPU)技术的出现, 为数值型计算提供了强大的计算能力, 并且其性价比较高, 在众多领域得到广泛应用[14-15]。时域成像算法中的后向投影算法又恰好是一种纯数值型计算, 导致其又重新引起了科研人员的重视。
为解决多子阵SAS精确后向投影成像算法效率低的问题, 本文提出了一种异构环境下的多子阵SAS精确后向投影成像快速算法, 满足高精度SAS实时成像需求。首先, 给出了多子阵SAS成像模型和后向投影成像算法; 然后, 针对后向投影成像算法中的距离向脉冲压缩和方位向聚焦过程进行了并行化设计。仿真和实测数据成像实验验证了所提算法的正确性和高效性。
-
图 1所示的直角坐标系中, SAS平台沿着方位向(垂直方向)运动, 目标位于坐标点
处, 发射阵和 个接收子阵沿方位向排列, 发射阵位于接收阵尾端, 其位置为 。 表示信号发射路径, 表示信号实际接收路径, 表示信号收发期间平台运行时间, 表示平台运动速度。为简化处理, 仅考虑第 个接收子阵的情况, 此接收子阵与发射阵方位向等效相位中心间隔为 。信号传播期间平台移动距离为 , 第 个接收子阵接收到的回波信号 可表示为: 式中,
表示回波信号幅度; 表示声速; 表示换能器时域加权; 表示调频斜率; 表示信号实际传播路径, 为去程, 为回程; 表示发射信号矩形窗包络; 表示点目标距离向坐标; 表示信号中心频率; 表示距离向快时间; 表示方位向慢时间; 表示脉冲宽度。 -
后向投影算法为典型的时域算法, 其思想来源于计算机层析成像, 它首先采用匹配滤波方式实现回波距离向上脉冲压缩, 然后在距离向脉冲压缩-方位时域, 采用同向叠加方式逐点实现每一像素方位向聚焦。它针对成像平面中每一像素进行遍历, 依次计算每一像素目标在各回波中的位置, 取出这些位置回波数据, 补偿相位差异后进行相干叠加, 同向回波得到增强, 进而获得当前点像素值。
距离向脉冲压缩通常在距离频域实现。首先将回波信号在距离向进行傅里叶变换(Fourier transform, FT), 然后将变换后的数据与距离向参考函数
相乘: 式中,
表示距离向瞬时频率; 表示调频斜率。最后进行傅里叶逆变换回到时域, 此时距离脉冲压缩后的回波可表示为: 式中,
表示信号带宽; 表示波长。 在进行方位向聚焦之前, 首先需要建立坐标系, 设计合理的成像网格, 一种特殊的网格大小设置是图像网格与回波网格大小相同。图 2为后向投影成像示意图, 建立
直角坐标系, 和 分别表示图像空间和距离脉冲压缩后的回波空间。 以成像空间中的点目标
为例介绍其重建过程, 其坐标为 , 回波空间中与其对应的点为 。系统工作过程中, 时刻回波中均包含点目标 回波信息, 对其重建均有贡献。将 时刻脉冲对 点贡献进行相干累积, 得到 点的重建结果, 即: 式中,
表示 时刻点目标 所对应的距离徙动量。 那么, 对于任意点目标
的重建表达式为: 式中,
和 分别表示目标开始照射和结束照射的时刻。 对于成像区域内的固定点目标
, 由于平台运动导致收发阵与固定点目标距离按 发生改变, 以致点目标在不同回波中的位置出现距离徙动现象, 如图 2所示。与单子阵成像方式相比, 多子阵形式下精确距离徙动量 不再呈现出对称的距离徙动形式, 不能再以理想状态下距离徙动量的特殊形式为基础设计快速成像算法。由于成像空间中每一点 的重建过程完全独立, 因此后向投影成像算法具有天然的并行性。但在早期受限于计算能力, 后向投影算法并未在大规模成像场景下得到应用。 -
后向投影成像算法从实现上来说主要包含两步: (1)距离向脉冲压缩; (2)方位向聚焦。
-
距离向脉冲压缩是对每一接收阵元数据依次进行快速傅里叶变换(fast Fourier transform, FFT), 与距离参考函数相乘和逆快速傅里叶变换(inverse fast Fourier transform, IFFT)的过程。为利用计算统一设备架构(compute unified device architecture, CUDA)中高效FFT函数库, 先调用cufftExecC2C对整块数据进行FFT, 以充分利用GPU的多核计算资源。距离向参考函数为一维数组, 为避免重复计算, 直接采用主机中央处理器(central processing unit, CPU)计算完毕后上传至显存。为避免处理过程中多次对FFT后的回波进行频谱搬移操作, 参考函数初始化时将其频率坐标范围设置为
, 其中 为距离向采样率。FFT后的回波与参考函数相乘过程如图 3所示, 其中回波数据块大小为 , 线程块大小为 , 当前线程块索引为 , 当前线程索引为 。 根据回波数据块和线程块大小, 推算出总线程块行列数分别为
和 , 其中 表示向下取整操作。为使得计算任务均匀分配, 每一线程完成单点数据乘法处理。根据线程索引与数据索引之间的关系, 令当前线程处理的数据索引为 , 则 , 。为防止核函数所访问的处理数据越界, 需满足 且 , 对访问数据下标进行合法性检测。 启动核函数, 将变换后的数据在距离向上同时与距离向参考函数
相乘后, 再次调用cufftExecC2C对整块数据进行IFFT, 完成距离向上的脉冲压缩。 -
完成回波数据距离向脉冲压缩后, 后向投影成像算法方位向聚焦实质是找准成像空间中每一点在回波数据中的位置, 提取数据, 补偿因距离徙动产生的相位差, 将相干叠加值作为当前像素点值。由于成像空间中每一像素点的重建过程完全相同, 本文仍然采用与脉冲压缩相同的二维线程块来完成方位向同向叠加, 每一线程完成成像空间内单个像素点的处理。当前线程
完成成像空间中像素点 的处理, 核函数具体执行步骤为: 1) 计算当前像素点
所对应的成像空间距离向坐标 。 2) 推算当前位置的合成孔径长度
, 确定目标点所在回波数据的起始行数 和终止行数 。 3) 根据当前处理数据行数推算接收数据子阵编号, 进而确定收发阵间隔
。 4) 根据收发阵间隔
和距离向坐标 计算精确时延 , 进而得到距离徙动量 。 5) 根据距离徙动量
, 取出对应最邻近的数据点, 补偿因距离徙动产生的相位差, 将其累加至当前像素点 上, 作为当前回波对成像结果的贡献。 6) 重复步骤1)~5), 直至当前孔径内所有脉冲回波数据完成同向叠加。
异构环境下的多子阵SAS精确后向投影成像算法信号处理流程如图 4所示。从流程图中可看出, CPU主要负责计算前的初始化工作, 以及将数据上传至GPU全局存储器。GPU主要负责FFT、IFFT、点乘等数据密集型计算。核函数执行完毕计算出结果后, CPU再从GPU全局存储器中取出计算结果, 整个成像过程完成。
-
为验证本文所提异构成像算法的性能, 进行了多子阵SAS回波数据成像实验, 实验环境为: 处理器为Intel(R) Xeon(R) CPU E5-2650 V3@2.3 GB (2个CPU, 20个内核, 40个逻辑处理器); 内存64 GB; 显卡为NVIDIA GeForce GTX 1660; 操作系统为Windows 10专业工作站版; 软件环境为VS2015+CUDA10。成像算法执行时间取10次执行时间的平均值。
-
仿真实验主要用于验证本文所提成像算法的正确性。仿真系统参数如表 1所示, 点目标位于距离向115 m、方位向38 m处。分别采用CPU串行成像方式和GPU并行成像方式对仿真点目标回波进行成像, 结果如图 5所示。图 5(a)为串行方式获得的点目标成像结果, 其距离向和方位向剖面图分别如图 5(b)和图 5(c)所示, 与理想剖面图一致。采用GPU并行成像方式获得的点目标成像结果如图 5(d)所示, 与之相对应的距离向和方位向剖面图分别如图 5(e)和图 5(f)所示, 其结果与串行成像方式获得的剖面图几乎相同。
表 1 仿真系统参数
Table 1. System Parameter for Simulation
参数 数值 参数 数值 带宽/kHz 20 脉冲间隔/ms 200 载频/kHz 150 子阵长/m 0.04 脉宽/ms 20 子阵个数 25 声速/(m∙s-1) 1 500 速度/(m∙s-1) 2.5 为进一步比较两种成像方式下点目标成像的性能, 对距离向和方位向冲击响应宽度(impulse response width, IRW)、峰值旁瓣比(peak sidelobe ratio, PSLR)和积分旁瓣比(integrated sidelobe ratio, ISLR)进行了测量, 其结果如表 2所示。从表 2中可看出, 两种成像方式下的点目标距离向指标相同, 方位向指标中PSLR存在细微差别, 主要原因在于两种成像方式下所使用的FFT函数不一样, 导致成像结果指标不完全相同, 但对成像结果影响不大, 进一步验证了本文所提并行成像算法的正确性。
表 2 仿真点目标成像性能比较
Table 2. Imaging Quality Comparison of Simulated Target
方法 距离向 方位向 IRW
/cmPSLR
/dBISLR
/dBIRW
/cmPSLR
/dBISLR
/dB串行 3.4 -13.7 -10.8 2.2 -13.0 -10.2 并行 3.4 -13.7 -10.8 2.2 -13.1 -10.2 -
实测数据实验主要用于验证并行成像算法的效率。实验数据选用2010年7月海军工程大学在浙江千岛湖进行干涉SAS海试样机实验所获取的两段原始回波数据, 系统工作时的基本参数如表 3所示。小场景原始数据块距离向点数为9 600, 方位向点数为2 880, 共72个脉冲, 数据采集时间为23.04 s。大场景数据块距离向点数为9 600, 方位向点数为10 368, 共216个脉冲, 数据采集时间为69.12 s。
表 3 实验系统基本参数
Table 3. System Parameter of the Trial
参数 数值 参数 数值 带宽/kHz 20 脉冲间隔/ms 320 载频/kHz 150 子阵长/m 0.04 脉宽/ms 20 子阵个数 40 声速/(m∙s-1) 1446 速度/(m∙s-1) 2.5 采样率/kHz 40 采样距离/m 51 231 分析实测数据效率时, 本文比较了串行成像算法、共享内存并行成像算法和异构并行成像算法3种方式下的计算时间, 其结果如表 4所示。共享内存并行成像算法是在串行成像算法基础上, 采用OpenMP并行指令, 利用多核计算资源协同加速成像过程。采用共享内存方式成像时, 计算核数与计算平台总的逻辑核数一致, 设置为40。采用异构平台成像时, 线程块大小设置为32×32。
表 4 不同成像方法效率比较
Table 4. Efficiency Comparison of the Different Imaging Algorithms
成像场景 成像方式 数据上传时间/ms 脉冲压缩时间/ms 方位叠加时间/ms 数据下载时间/ms 总成像时间/ms 加速比 小场景 串行 3 541 3 153 097 3 156 638 1 共享内存并行 298 148 049 148 347 21.3 异构并行 50 33 9 553 38 9 674 326.3 大场景 串行 10 025 10 555 120 10 565 145 1 共享内存并行 641 456 303 456 944 23.1 异构并行 151 90 31 311 116 31 668 333.6 采用异构环境下的成像算法分别对小场景数据和大场景数据进行成像处理, 结果如图 6(a)和图 6(b)所示, 从成像结果中可以看出小场景区域为大场景区域中的局部区域, 成像算法成功实现了两块成像区域的水下地形清晰重构, 图 6中道路和水塘清晰可见, 验证了本文算法的正确性。
3种成像方法效率比较见表 4, 串行成像方式所需时间最长, 成像总时间分别为3 156 638 ms和10 565 145 ms。成像算法并行前后, 成像所需时间主要消耗在方位向叠加处理步骤。以串行算法成像时间为基准, 采用共享内存方式并行后, 总成像时间显著减少, 总成像时间分别降为148 347 ms和456 944 ms, 加速比分别达到21.3和23.1。但与原始回波数据获取时间相比, 仍然不能满足实时数据处理需求。采用异构方式并行计算后, 与串行和共享内存并行成像方式相比, 虽然增加了原始数据的上传和计算结果的下载操作, 但因传输带宽高, 所消耗的时间分别仅为88 ms和267 ms。方位向叠加时间也分别进一步降低到9 553 ms和31 311 ms, 总成像时间分别仅为9 674 ms和31 668 ms。与串行方式下的成像时间相比, 异构并行成像方法的加速比分别高达326.3和333.6, 远小于各自原始数据获取时间, 满足实时合成孔径成像需求。
-
本文提出了一种异构环境下的多子阵SAS精确逐点成像快速算法, 在分析精确逐点成像算法中后向投影算法的基础上, 将成像过程中脉冲压缩和方位向叠加过程采用GPU的众核计算资源计算以达到快速计算目的, 并通过点目标仿真实验验证了所提并行成像方法的正确性。实测数据的成像实验结果表明, 与传统串行成像方式相比, 并行成像加速比分别高达326.3和333.6, 满足实时SAS成像需要。为进一步提升算法适用性, 后续主要开展逐点成像算法与运动补偿集成方法, 以及多GPU环境下的并行化方法研究。
A Fast Accurate Back-Projection Algorithm for Multi-receiver Synthetic Aperture Sonar in Heterogeneous Environment
-
摘要: 针对多子阵合成孔径声呐精确后向投影成像算法效率低的问题, 提出了一种异构环境下的精确多子阵合成孔径声呐后向投影成像快速方法。在分析精确逐点后向投影成像算法原理的基础上, 将脉冲压缩和方位向聚焦过程改造为单指令多线程模式, 借助图形处理器(graphics processor unit, GPU)强大的多核计算能力加速成像过程。通过仿真和实测数据的成像实验验证了所提快速成像算法的正确性和高效性, 与串行成像算法相比, 其加速比分别高达326.3和333.6。对于大规模数据成像处理, 所提方法体现出优异的加速性能, 满足实时信号处理需求, 同时为后续开展运动补偿奠定了基础。Abstract:
Objectives Synthetic aperture sonar (SAS) is one kind of high resolution underwater imaging sonar. Its principle is to use small array to simulate large aperture array by uniforming linear motion along the direction of flight track, and realize coherent processing of echo signal to obtain high resolution two-dimensional sonar image, which is independent of imaging distance and acoustic wavelength. The back projection algorithm is an accurate point-by-point imaging algorithm, which is characterized by accurate focusing ability, suitable for wide beam and large scene imaging, and easy to realize motion compensation, but it has the disadvantages of low computational efficiency. In recent years, with the development of computer technology, especially the appearance of graphics processor unit, it has powerful computing power and provides a new way to accelerate the back projection algorithm. Methods We propose a fast accurate back projection algorithm for multi-receiver SAS in heterogeneous environment. On the basis of analyzing the principle of accurate back-projection imaging algorithm, the pulse compression and time-consuming azimuth accumulation are transformed into a single instruction multi-threaded mode, and the imaging process is accelerated by the powerful multi-core computing ability of graphics processor unit. Results The validity and efficiency of the proposed fast imaging algorithm are verified by the imaging experiments performed on simulated and real SAS data. Compared with the serial imaging algorithm, its acceleration ratio is as high as 326.3 and 333.6. For large-scale data imaging processing, it shows excellent acceleration performance. Conclusions The proposed back projection imaging algorithm greatly improves the imaging efficiency and meets the needs of real-time SAS imaging, which lays the foundation for subsequent motion compensation. -
表 1 仿真系统参数
Table 1. System Parameter for Simulation
参数 数值 参数 数值 带宽/kHz 20 脉冲间隔/ms 200 载频/kHz 150 子阵长/m 0.04 脉宽/ms 20 子阵个数 25 声速/(m∙s-1) 1 500 速度/(m∙s-1) 2.5 表 2 仿真点目标成像性能比较
Table 2. Imaging Quality Comparison of Simulated Target
方法 距离向 方位向 IRW
/cmPSLR
/dBISLR
/dBIRW
/cmPSLR
/dBISLR
/dB串行 3.4 -13.7 -10.8 2.2 -13.0 -10.2 并行 3.4 -13.7 -10.8 2.2 -13.1 -10.2 表 3 实验系统基本参数
Table 3. System Parameter of the Trial
参数 数值 参数 数值 带宽/kHz 20 脉冲间隔/ms 320 载频/kHz 150 子阵长/m 0.04 脉宽/ms 20 子阵个数 40 声速/(m∙s-1) 1446 速度/(m∙s-1) 2.5 采样率/kHz 40 采样距离/m 51 231 表 4 不同成像方法效率比较
Table 4. Efficiency Comparison of the Different Imaging Algorithms
成像场景 成像方式 数据上传时间/ms 脉冲压缩时间/ms 方位叠加时间/ms 数据下载时间/ms 总成像时间/ms 加速比 小场景 串行 3 541 3 153 097 3 156 638 1 共享内存并行 298 148 049 148 347 21.3 异构并行 50 33 9 553 38 9 674 326.3 大场景 串行 10 025 10 555 120 10 565 145 1 共享内存并行 641 456 303 456 944 23.1 异构并行 151 90 31 311 116 31 668 333.6 -
[1] Hayes M P, Gough P T. Synthetic Aperture Sonar: A Review of Current Status[J]. IEEE Journal of Oceanic Engineering, 2009, 34(3): 207-224 doi: 10.1109/JOE.2009.2020853 [2] 刘经南, 阳凡林, 赵建虎. 浅析合成孔径声纳与干涉合成孔径声纳[J]. 海洋测绘, 2003, 23(2): 1-4 doi: 10.3969/j.issn.1671-3044.2003.02.001 Liu Jingnan, Yang Fanlin, Zhao Jianhu. Elementary Introduction to Synthetic Aperture Sonar and Interferometric Synthetic Aperture Sonar[J]. Hydrographic Surveying and Charting, 2003, 23(2): 1-4 doi: 10.3969/j.issn.1671-3044.2003.02.001 [3] Wang V T, Hayes M P. Synthetic Aperture Sonar Track Registration Using SIFT Image Correspondences[J]. IEEE Journal of Oceanic Engineering, 2017, 42(4): 901-913 doi: 10.1109/JOE.2016.2634078 [4] Myers V, Quidu I, Zerr B, et al. Synthetic Aperture Sonar Track Registration with Motion Compensation for Coherent Change Detection[J]. IEEE Journal of Oceanic Engineering, 2019, 45(3): 1045-1062 [5] Ødegård Ø, Hansen R E, Singh H, et al. Archaeological Use of Synthetic Aperture Sonar on Deepwater Wreck Sites in Skagerrak[J]. Journal of Archaeological Science, 2018, 89: 1-13 doi: 10.1016/j.jas.2017.10.005 [6] Ulander L M H, Hellsten H, Stenstrom G. SyntheticAperture Radar Processing Using Fast Factorized Back-Projection[J]. IEEE Transactions on Aerospace and Electronic Systems, 2003, 39 (3) : 760-776 doi: 10.1109/TAES.2003.1238734 [7] Yegulalp A F. Fast Back Projection Algorithm for Synthetic Aperture Radar[C]//IEEE Radar Conference, Waltham, MA, USA, 1999 [8] Wu H R, Tang J S, Zhong H P. Moderate Squint Imaging Algorithm for the Multiple-Hydrophone SAS with Receiving Hydrophone Dependence[J]. IET Radar, Sonar and Navigation, 2019, 13(1): 139-147 doi: 10.1049/iet-rsn.2018.5055 [9] Tian Z, Tang J S, Zhong H P, et al. Extended Range Doppler Algorithm for Multiple-Receiver Synthetic Aperture Sonar Based on Exact Analytical Two-Dimensional Spectrum[J]. IEEE Journal of Oceanic Engineering, 2016, 41(1): 164-174 doi: 10.1109/JOE.2015.2402053 [10] Zhang X B, Tang J S, Zhong H P. Multireceiver Correction for the Chirp Scaling Algorithm in Synthetic Aperture Sonar[J]. IEEE Journal of Oceanic Engineering, 2014, 39(3): 472-481 doi: 10.1109/JOE.2013.2251809 [11] 王金波, 唐劲松, 张森, 等. 一种宽带大斜视STOLT插值及距离变标补偿方法[J]. 电子与信息学报, 2018, 40(7): 1575-1582 https://www.cnki.com.cn/Article/CJFDTOTAL-DZYX201807008.htm Wang Jinbo, Tang Jinsong, Zhang Sen, et al. Range Scaling Compensation Method Based on STOLT Interpolation in Broadband Squint SAS Imaging[J]. Journal of Electronics and Information Technology, 2018, 40(7): 1575-1582 https://www.cnki.com.cn/Article/CJFDTOTAL-DZYX201807008.htm [12] Li W J, Liao G S, Zhu S Q, et al. A Novel Helicopter-Borne RoSAR Imaging Algorithm Based on the Azimuth Chirp Z Transform[J]. IEEE Geoscience and Remote Sensing Letters, 2019, 16(2): 226-230 doi: 10.1109/LGRS.2018.2871379 [13] Li B L, Liu W, Liu J Y, et al. Real-Time Implementation of Synthetic Aperture Sonar Imaging on High Performance Clusters[C]//The11th ACIS International Conference on Software Engineering, Artificial Intelligence, Networking and Parallel/Distributed Computing, London, UK, 2010 [14] 方留杨, 王密, 潘俊. CPU和GPU协同的多光谱影像快速波段配准方法[J]. 武汉大学学报·信息科学版, 2018, 43(7): 1000-1007 doi: 10.13203/j.whugis20160218 Fang Liuyang, Wang Mi, Pan Jun. CPU/GPU Cooperative Fast Band Registration Method for Multispectral Imagery[J]. Geomatics and Information Science of Wuhan University, 2018, 43 (7) : 1000-1007 doi: 10.13203/j.whugis20160218 [15] 王鸿琰, 关雪峰, 吴华意. 一种面向CPU/GPU异构环境的协同并行空间插值算法[J]. 武汉大学学报·信息科学版, 2017, 42(12): 1688-1695 doi: 10.13203/j.whugis20150361 Wang Hongyan, Guan Xuefeng, Wu Huayi. A Collaborative Parallel Spatial Interpolation Algorithmon Oriented Towards the Heterogeneous CPU/GPU System[J]. Geomatics and Information Science of Wuhan University, 2017, 42(12): 1688-1695 doi: 10.13203/j.whugis20150361 -