图 1所示的直角坐标系中, SAS平台沿着方位向(垂直方向)运动, 目标位于坐标点$(r,\mathrm{0})$处, 发射阵和$N$个接收子阵沿方位向排列, 发射阵位于接收阵尾端, 其位置为$(\mathrm{0},vt)$。$R_{\mathrm{1}}$表示信号发射路径, $R_{\mathrm{2}}$表示信号实际接收路径, $t^{\mathrm{*}}$表示信号收发期间平台运行时间, $v$表示平台运动速度。为简化处理, 仅考虑第$i$个接收子阵的情况, 此接收子阵与发射阵方位向等效相位中心间隔为$\mathrm{\Delta }{h}_i$。信号传播期间平台移动距离为$v{t}^{\mathrm{*}}$, 第$i$个接收子阵接收到的回波信号$s_i(\tau ,t,r)$可表示为:

图  1  SAS成像模型

Figure 1.  The Mode of SAS Imaging

式中, $A_{\mathrm{0}}$表示回波信号幅度; $c$表示声速; ${\omega }_a\left(t\right)$表示换能器时域加权; $k$表示调频斜率; $R_i\left(t,r\right)=R_{\mathrm{1}}+R_{\mathrm{2}}$表示信号实际传播路径, $R_{\mathrm{1}}=\sqrt[]{{(vt)}^{\mathrm{2}}+r^{\mathrm{2}}}$为去程, $R_{\mathrm{2}}=$$\sqrt[]{(vt+\mathrm{\Delta }{h}_i+v{t}^{\mathrm{*}}{)}^{\mathrm{2}}+r^{\mathrm{2}}}$为回程; $G\left(\cdot \right)$表示发射信号矩形窗包络; $r$表示点目标距离向坐标; $f_{\mathrm{0}}$表示信号中心频率; $\tau$表示距离向快时间; $t$表示方位向慢时间; $T$表示脉冲宽度。

后向投影算法为典型的时域算法, 其思想来源于计算机层析成像, 它首先采用匹配滤波方式实现回波距离向上脉冲压缩, 然后在距离向脉冲压缩-方位时域, 采用同向叠加方式逐点实现每一像素方位向聚焦。它针对成像平面中每一像素进行遍历, 依次计算每一像素目标在各回波中的位置, 取出这些位置回波数据, 补偿相位差异后进行相干叠加, 同向回波得到增强, 进而获得当前点像素值。

距离向脉冲压缩通常在距离频域实现。首先将回波信号在距离向进行傅里叶变换(Fourier transform, FT), 然后将变换后的数据与距离向参考函数$H_r$相乘:

式中, $f_r$表示距离向瞬时频率; $k$表示调频斜率。最后进行傅里叶逆变换回到时域, 此时距离脉冲压缩后的回波可表示为:

式中, $B$表示信号带宽; $\lambda$表示波长。

在进行方位向聚焦之前, 首先需要建立坐标系, 设计合理的成像网格, 一种特殊的网格大小设置是图像网格与回波网格大小相同。图 2为后向投影成像示意图, 建立$(\tau ,t)$直角坐标系, $I(\tau ,t)$和$s(\tau ,t)$分别表示图像空间和距离脉冲压缩后的回波空间。

图  2  后向投影成像示意图

Figure 2.  Diagram of Back Projection Algorithm

以成像空间中的点目标$P$为例介绍其重建过程, 其坐标为$({\tau }_{\mathrm{0}},t_{\mathrm{0}})$, 回波空间中与其对应的点为$s({\tau }_{\mathrm{0}},t_{\mathrm{0}})$。系统工作过程中, $t_{-\mathrm{4}}~t_{\mathrm{4}}$时刻回波中均包含点目标$P$回波信息, 对其重建均有贡献。将$t_{-\mathrm{4}}~t_{\mathrm{4}}$时刻脉冲对$P$点贡献进行相干累积, 得到$P$点的重建结果, 即:

式中, $\mathrm{\Delta }r=c(t^{\mathrm{*}}-{\tau }_{\mathrm{0}})/\mathrm{2}$表示$t_m$时刻点目标$P$所对应的距离徙动量。

那么, 对于任意点目标$(\tau ,t)$的重建表达式为:

式中, $t_s$和$t_e$分别表示目标开始照射和结束照射的时刻。

对于成像区域内的固定点目标$P$, 由于平台运动导致收发阵与固定点目标距离按$R_i\left(t,r\right)=R_{\mathrm{1}}+R_{\mathrm{2}}$发生改变, 以致点目标在不同回波中的位置出现距离徙动现象, 如图 2所示。与单子阵成像方式相比, 多子阵形式下精确距离徙动量$\mathrm{\Delta }r$不再呈现出对称的距离徙动形式, 不能再以理想状态下距离徙动量的特殊形式为基础设计快速成像算法。由于成像空间中每一点$(\tau ,t)$的重建过程完全独立, 因此后向投影成像算法具有天然的并行性。但在早期受限于计算能力, 后向投影算法并未在大规模成像场景下得到应用。

后向投影成像算法从实现上来说主要包含两步: (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后的回波进行频谱搬移操作, 参考函数初始化时将其频率坐标范围设置为$\left[\mathrm{0},F_s\right]$, 其中$F_s$为距离向采样率。FFT后的回波与参考函数相乘过程如图 3所示, 其中回波数据块大小为$N_r\times N_a$, 线程块大小为$T_r\times T_a$, 当前线程块索引为$(B_i,B_j)$, 当前线程索引为$(T_i,T_j)$。

图  3  距离向回波与参考函数相乘

Figure 3.  Echo Multiplied with Reference in Range

根据回波数据块和线程块大小, 推算出总线程块行列数分别为$F((N_a+T_a-\mathrm{1})/T_a)$和$F((N_r+T_r-\mathrm{1})/T_r)$, 其中$F(\cdot )$表示向下取整操作。为使得计算任务均匀分配, 每一线程完成单点数据乘法处理。根据线程索引与数据索引之间的关系, 令当前线程处理的数据索引为$(C_i,C_j)$, 则$C_i=B_i\cdot T_r+T_i$, $C_j=B_j\cdot T_a+T_j$。为防止核函数所访问的处理数据越界, 需满足$C_i<N_r$且$C_j<N_a$, 对访问数据下标进行合法性检测。

启动核函数, 将变换后的数据在距离向上同时与距离向参考函数$H_r$相乘后, 再次调用cufftExecC2C对整块数据进行IFFT, 完成距离向上的脉冲压缩。

完成回波数据距离向脉冲压缩后, 后向投影成像算法方位向聚焦实质是找准成像空间中每一点在回波数据中的位置, 提取数据, 补偿因距离徙动产生的相位差, 将相干叠加值作为当前像素点值。由于成像空间中每一像素点的重建过程完全相同, 本文仍然采用与脉冲压缩相同的二维线程块来完成方位向同向叠加, 每一线程完成成像空间内单个像素点的处理。当前线程$(T_i,T_j)$完成成像空间中像素点$I(T_i,T_j)$的处理, 核函数具体执行步骤为:

1) 计算当前像素点$I(T_i,T_j)$所对应的成像空间距离向坐标$r_p$。

2) 推算当前位置的合成孔径长度$L$, 确定目标点所在回波数据的起始行数$L_s$和终止行数$L_e$。

3) 根据当前处理数据行数推算接收数据子阵编号, 进而确定收发阵间隔$\mathrm{\Delta }{h}_i$。

4) 根据收发阵间隔$\mathrm{\Delta }{h}_i$和距离向坐标$r_p$计算精确时延$t^{\mathrm{*}}$, 进而得到距离徙动量$\mathrm{\Delta }r$$=$$c{t}^{\mathrm{*}}/\mathrm{2}-r_p$。

5) 根据距离徙动量$\mathrm{\Delta }r$, 取出对应最邻近的数据点, 补偿因距离徙动产生的相位差, 将其累加至当前像素点$I(T_i,T_j)$上, 作为当前回波对成像结果的贡献。

6) 重复步骤1)~5), 直至当前孔径内所有脉冲回波数据完成同向叠加。

异构环境下的多子阵SAS精确后向投影成像算法信号处理流程如图 4所示。从流程图中可看出, CPU主要负责计算前的初始化工作, 以及将数据上传至GPU全局存储器。GPU主要负责FFT、IFFT、点乘等数据密集型计算。核函数执行完毕计算出结果后, CPU再从GPU全局存储器中取出计算结果, 整个成像过程完成。

图  4  异构环境下的成像算法流程

Figure 4.  Flowchart of Imaging Algorithm in Heterogeneous Environment

为验证本文所提异构成像算法的性能, 进行了多子阵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)所示, 其结果与串行成像方式获得的剖面图几乎相同。

图  5  仿真点目标成像结果

Figure 5.  Imaging Results for the Simulated Target

为进一步比较两种成像方式下点目标成像的性能, 对距离向和方位向冲击响应宽度(impulse response width, IRW)、峰值旁瓣比(peak sidelobe ratio, PSLR)和积分旁瓣比(integrated sidelobe ratio, ISLR)进行了测量, 其结果如表 2所示。从表 2中可看出, 两种成像方式下的点目标距离向指标相同, 方位向指标中PSLR存在细微差别, 主要原因在于两种成像方式下所使用的FFT函数不一样, 导致成像结果指标不完全相同, 但对成像结果影响不大, 进一步验证了本文所提并行成像算法的正确性。

实测数据实验主要用于验证并行成像算法的效率。实验数据选用2010年7月海军工程大学在浙江千岛湖进行干涉SAS海试样机实验所获取的两段原始回波数据, 系统工作时的基本参数如表 3所示。小场景原始数据块距离向点数为9 600, 方位向点数为2 880, 共72个脉冲, 数据采集时间为23.04 s。大场景数据块距离向点数为9 600, 方位向点数为10 368, 共216个脉冲, 数据采集时间为69.12 s。

分析实测数据效率时, 本文比较了串行成像算法、共享内存并行成像算法和异构并行成像算法3种方式下的计算时间, 其结果如表 4所示。共享内存并行成像算法是在串行成像算法基础上, 采用OpenMP并行指令, 利用多核计算资源协同加速成像过程。采用共享内存方式成像时, 计算核数与计算平台总的逻辑核数一致, 设置为40。采用异构平台成像时, 线程块大小设置为32×32。

采用异构环境下的成像算法分别对小场景数据和大场景数据进行成像处理, 结果如图 6(a)和图 6(b)所示, 从成像结果中可以看出小场景区域为大场景区域中的局部区域, 成像算法成功实现了两块成像区域的水下地形清晰重构, 图 6中道路和水塘清晰可见, 验证了本文算法的正确性。

图  6  实测数据成像结果

Figure 6.  Image Results for the Real Data

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环境下的并行化方法研究。