CN106291551A - 一种基于gpu的并行结构isar距离对准方法 - Google Patents

Abstract

本发明公开了一种基于GPU的并行结构ISAR距离对准方法，主要解决传统的基于CPU单线程的ISAR距离对准算法效率低的问题，该方法主要步骤包括：步骤一：读取原始回波数据；步骤二：采用流水线异步并行机制，将原始回波数据分流从主机内存拷贝到设备显存；步骤三：距离向压缩，包括距离向FFT、点乘匹配滤波以及IFFT；步骤四：全局最小熵距离对准，包括时域插值、求平均距离像、求距离偏移量、构造频移因子、亚距离单元距离对准；步骤五：将距离对准数据分流从GPU拷贝到CPU；步骤六：利用FreeImage库函数输出图像。本发明利用并行化硬件平台GPU对ISAR距离对准算法进行硬件加速，运算速度得到几十倍的提高。

Description

一种基于GPU的并行结构ISAR距离对准方法

技术领域

本发明涉及一种并行结构ISAR距离对准方法，特别是涉及一种在CUDA平台上用GPU进行ISAR距离对准并行加速算法设计。

背景技术

逆合成孔径雷达(Inverse-Synthetic-Aperture-Radar，简称ISAR)距离对准是消除目标相对于雷达平动造成的相邻回波在距离向上的错位，对准后，目标上同一散射点的信号在不同的回波脉冲中都位于同一距离单元中。目前多数距离对准算法是基于目标相邻回波间相似性，其中较具代表性的有包络互相关法、积累互相关法等。这种相邻回波相关法，是以前一个(几个)一维距离像作为基准，认为前一个(几个)一维像已经严格对准，估计出当前一维距离像相对于前一个(几个)距离像的距离走动。这类算法依赖于距离像之间的相关性，在相关性不强的情况下对准效果不明显，并且一维距离像对准是时间先后的串行机制。

全局最小熵是一种基于平均距离像(Average range profile，简称ARP)熵最小原则的距离对准算法。假设p(r,n)是距离脉压以后的回波信号包络的幅度，其中r代表的是雷达视线(radar line-of-sight，简称RLOS)的距离，n为回波的序号(0≤n≤N)，N是总的回波个数。定义平均距离像ARP是所有脉冲包络的迭加：

$p_{ave}\left(r\right)=\underset{n=0}{\overset{N-1}{\Σ}}p(r+{\mathrm{\Δ}}_r(n),n)---\left(1\right)$

上式中，△_r(n)＝[△_r(0),△_r(1),…,△_r(N-1)]，△_r(n)是第n个回波的距离偏移量。

平均距离像的熵为

$E\left({\mathrm{\Δ}}_r\right)=-\underset{r}{\∫}\frac{p_{ave}\left(r\right)}{S}\·ln\frac{p_{ave}\left(r\right)}{S}dr---\left(2\right)$

其中

$S=\underset{r}{\∫}{p}_{ave}\left(r\right)dr---\left(3\right)$

平均距离像熵值最小转化为求取距离偏移量，使其满足：

${\mathrm{\Δ}}_r\left(n\right)=\arg \{max\[p(r,n)\⊗ln{P}_{ave}(-r)\]\},n=0,1,...,N-1---\left(4\right)$

其中，表示在距离向一维卷积运算。

这种基于平均距离像熵值最小的距离对准算法，主要包含脉冲包络迭加、平均距离像与一维距离像互相关运算、以及求取矩阵最大值问题，其中互相关运算可通过快速傅里叶变换频域共轭相乘实现。每个一维像的距离对准是独立进行的，可以实现并行优化。

以往，ISAR成像距离对准运算主要用软件方式实现。随着技术的发展，工作者需要处理的数据量越来越大，对成像系统的实时性要求也越来越高，软件实现ISAR距离对准的速度远达不到要求。近年来，图形处理器(Graphics Processing Unit，简称GPU)的迅速发展为大规模数据的计算与分析提供了一种新的技术手段。在此基础上，各种合成孔径雷达(Synthetic-Aperture-Radar，简称SAR)成像处理算法在GPU上得到加速实现。目前基于GPU的ISAR距离对准，如文献：刘百玲,江海清,倪书爱等.基于GPU的ISAR成像算法实现[J].电子测量技术,2015,(8):76-78,89.，其实现的距离对准算法是积累互相关法，基本思想是将待对齐的脉冲与已对准后的相邻几个脉冲相关处理。这种下一个待对准的脉冲需依赖上几个已对准脉冲的串行机制距离对准方法，并不能很好利用GPU资源。

发明内容

发明目的：为了克服现有技术中存在的不足，本发明提供提出一种基于GPU的并行结构ISAR距离对准方法，用于解决现有的基于GPU的ISAR算法不能很好地利用GPU资源导致对准的精度和效率均不高的技术问题。

技术方案：为实现上述目的，本发明采用的技术方案为：

一1、一种基于GPU的并行结构ISAR距离对准方法，其特征在于：包括顺序执行的以下步骤：

步骤一、CPU读取ISAR原始回波数据，大小为Na*Nr，Na为方位向采样点数，Nr为距离向采样点数；

步骤二、采用统一计算设备架构流水线异步并行机制进行如下的分流脉冲压缩过程：将步骤一中读取的原始回波数据按方位向分块后分配给不同的工作流处理，每个工作流单独完成将本块数据从CPU拷贝到GPU中并在GPU中进行距离向压缩的工作得到脉压数据；

步骤三、在GPU中将脉压数据并行进行全局最小熵距离对准；

步骤四、将完成距离对准的数据分流从GPU中拷贝到CPU；

步骤五、利用FreeImage库函数输出图像。

进一步的，在本发明中，所述步骤三包括如下过程：

步骤3.1、将脉压数据转置为Nr列Na行，然后变换到距离向频域补零，再逆傅里叶变换到距离向时域，得到距离向时域插值后的TT*Nr列Na行的脉压数据，其中TT是插值倍数；

步骤3.2、对每个线程块的x维度和y维度均分配TILE_DIM个线程，每个线程网格的x维度分配Na/TILE_DIM个线程块、y维度分配TT*Nr/TILE_DIM个线程块，TILE_DIM固定为16，编写取模值核函数求插值一维距离像包络；

步骤3.3、编写转置核函数，将插值一维距离像包络的矩阵转置为Na列TT*Nr行；

步骤3.4、利用共享内存和归约算法编写包络求和核函数对一维距离像包络按方位向求和，分两步计算：第一步，对每个线程网格x维度分配TT*Nr个线程块、y维度分配Na/threadnum个线程块，每个线程块x维度分配threadnum个线程、y维度分配1个线程，threadnum取值为64，求得TT*Nr*Na/threadnum个和值；第二步，对每个线程网格x维度分配TT*Nr个线程块、y维度分配1个线程块，每个线程块的x维度分配Na/threadnum个线程、y维度分配1个线程，对上一步中求得的TT*Nr*Na/threadnum个和值继续求和，得到TT*Nr个和元素，即为平均距离像ARP；

步骤3.5、对每个线程网格x维度分配TT*Nr/TILE_DIM个线程块、y维度分配1个线程块，对每个线程块x维度分配为TILE_DIM个线程、y维度分配1个线程，编写核函数求平均距离像ARP的对数

步骤3.6、对每个线程网格x维度分配TT*Nr/TILE_DIM个线程块、y维度分配Na/TILE_DIM个线程块，对每个线程块x维度和y维度均分配TILE_DIM个线程，编写与一维距离像互相关核函数，通过频域共轭相乘实现互相关运算；

步骤3.7、利用共享内存和归约算法编写核函数求互相关矩阵每行最大值，分两步计算：第一步，对每个线程网格x维度分配TT*Nr/threadnum个线程块、y维度分配Na个线程块，对每个线程块x维度分配threadnum个线程、y维度分配1个线程，求得Na*TT*Nr/threadnum个最大值及其列索引；第二步，对每个线程网格x维度分配Na个线程块、y维度分配1个线程块，对每个线程块x维度分配TT*Nr/threadnum个线程、y维度分配1个线程，对上一步中求得的Na行TT*Nr/threadnum列最大值矩阵求每行最大值及其列索引，从而求得Na个距离偏移量；

步骤3.8、重复3.2～3.7的迭代过程，当达到设定迭代次数后，利用最近一次迭代所得到的距离偏移量构造频移因子exp(-j2πm△_r(n)/(TTN_r))，m是距离单元标号，△_r(n)是第n个回波的距离偏移量；

步骤3.9、对每个线程网格x维度分配Na/TILE_DIM个线程块、y维度分配Nr/TILE_DIM个线程块，对每个线程块x维度和y维度均分配TILE_DIM个线程，编写未插值一维距离像亚距离单元对准核函数，每个线程索引独立完成未插值一维距离像第n个脉冲第m个距离门处元素值与频移因子exp(-j2πm△_r(n)/(TTN_r))的相乘运算，完成未插值一维距离像的亚距离单元对准；

步骤3.10、将亚距离单元对准后的未插值一维距离像距离向逆傅里叶变换，再转置为Na列Nr行的距离对准数据。

有益效果：

本发明将CUDA平台下CPU+GPU的异构模式引入ISAR距离对准，实现多线程并行结构的亚距离单元距离对准，充分利用GPU并行运算能力，提高了ISAR距离对准的精度和效率，使得基于全局最小熵的亚距离单元距离对准算法能够达到数十倍的加速比的提升，其处理速度可以达到实时成像的要求。具体体现在以下几点：

(1)将原始回波数据分子块，利用CUDA多流异步执行技术，进行距离匹配滤波压缩，并行执行数据拷贝和核函数运行，有效隐藏了数据在主机和设备之间的传输时间。

(2)脉压数据插值处理后，利用插值矩阵求取距离偏移量，构造未插值矩阵的频移因子，实现亚距离单元的距离对准。

(3)每个一维距离像的对准都是与平均距离像相关，相关向量最大值即为距离偏移量。通过连续执行多个内核函数，并行实现全部脉冲的距离对准。

附图说明

图1为本发明流程图；

图2是传统的基于CPU的亚距离单元对准的全局最小熵距离对准算法流程图；

图3是卫星仿真散点分布图；

图4是CPU距离对准图像；

图5是GPU距离对准图像；

图6是CPU距离对准后相位梯度自聚焦补偿后最终图像；

图7是GPU距离对准后相位梯度自聚焦补偿后最终图像。

具体实施方式

下面结合附图对本发明作更进一步的说明。

本发明提出的一种基于GPU的并行结构ISAR距离对准方法的流程图如图1所示。主要流程包括原始回波数据传输、距离向压缩、距离向插值、求距离偏移量、构造平移因子、亚距离单元距离对准、对准数据传回以及图像输出。

首先，原始回波数据传输和距离向压缩属于分流脉冲压缩过程，具体包括如下过程：

CPU读取ISAR原始回波数据，大小为Na×Nr，Na为方位向采样点数，Nr为距离向采样点数；然后采用统一计算设备架构流水线即CUDA的异步并行机制，将读取的原始回波数据分流从CPU中拷贝到GPU中；接着在GPU中对原始回波数据进行距离向压缩，包括距离向FFT、点乘匹配滤波器滤波以及距离向IFFT三步。

为使说明简单明确，我们创建4个流。将ISAR原始回波数据按方位向分块，并分别交由流0、流1、流2、流3处理，每个流单独完成本块数据的脉冲压缩操作，包括CPU到GPU数据拷贝、距离向FFT、点乘匹配滤波器滤波、距离向IFFT共4个过程。经过实验测算，单独CPU处理脉冲压缩耗时2321.0001ms，基于GPU的不分流脉冲压缩耗时69.682ms，利用本发明优化后的分流脉冲压缩耗时48.541ms且相比于CPU耗时获得了近48倍的加速比。

接下来，对传统基于CPU的亚距离单元对准的全局最小熵距离对准算法流程进行分析。如图2所示，设置每个脉冲回波距离偏移量△_r(n)的初始值为零，初始的平均距离像ARP是所有未移位的输入距离像包络求和，对初始平均距离像求其对数将和每个脉冲包络FFT共轭相乘实现互相关运算，使互相关函数绝对值最大的△_r(n)就是满足该次迭代要求的距离偏移量。将一维距离像循环移位△_r(n)，更新平均距离像，若平均距离像的熵值有所减少，重复上述迭代过程求得距离偏移量；当平均距离像熵值不再减少时停止迭代，此轮迭代得到距离偏移量△_r(n)的最终估计值。

接下来，重点分析全局最小熵距离对准算法的可并行化。根据上述分析，全局最小熵距离对准算法主要涉及对脉压数据频域补零、距离向FFT/IFFT、一维距离像包络求和、互相关运算、最大值运算、乘以频移因子运算，其中频域补零是为实现时域插值运算，利用傅氏变换的时频移位对称性质，通过频域乘以线性相位因子即频移因子exp(-j2πm△_r(n)/(TTNr))(m是距离单元标号，△_r(n)是第n个回波的距离偏移量)；实现亚距离单元距离对准。其中，迭代运算可在CPU上运行，其它大矩阵的简单数值计算非常适合GPU并行计算，因此利用CUDA对其进行加速。

最后，详细阐述GPU如何加速全局最小熵距离对准算法。

首先将脉压数据传输到GPU内存，调用核函数转置(Nr列Na行，Na为方位向点数，Nr为距离向点数)，按距离向进行频域补零，再IFFT变换回到距离时域，实现离向在时域的插值运算。

接着，按脉压数据的大小分配线程块和线程，调用核函数完成一维距离像包络求和，得到平均距离像。再分配TT×Na×Nr个线程，每个线程独立访问平均距离像和Na个一维距离像相应位置元素，利用频域共轭相乘实现互相关运算。接下来，合理分配程块和线程，调用核函数寻找互相关矩阵每行(每个一维距离像)最大值，对应的列索引即为距离偏移量，将补零矩阵循环移位。迭代几次后，得到最终估计的距离偏移量。最后，分配Na×Nr个线程，每个线程索引独立完成与距离偏移量构成的频移因子的相乘。

在基于全局最小熵的亚距离单元距离对准求取平均距离像时，每个距离单元需要对所有脉冲求和，可以使用共享内存和并行归约实现，分配TT×Nr个线程块，每个块分配Na个线程，将每个距离单元数据传入共享内存中，使用并行归约计算块内元素和。考虑到实际情况Na可能超过线程块最大线程数，因此用两步核函数并行归约计算脉冲包络和。先求得Na×TT×Nr/threadnum个元素和值，再求出Na个一维距离像包络和(threadnum为第一步核函数线程块内分配的线程数)。距离对准中求取互相关矩阵最大值时，采用了类似的并行归约算法。大规模并行(128个线程以及256个线程)运算时，有效地优化程序性能。

下面通过卫星仿真数据进一步说明本发明的有效性。图3是基于IDL软件的卫星仿真散点分布图，利用STK软件获得运动卫星与地面雷达瞬时斜距，从而在IDL中实现回波模拟。雷达参数如表1所示：

表1

GPU处理型号是Tesla C2075(计算能力2.0,1.15GHz时钟频率，6GB GDDR5显存，448个CUDA计算核心，1.03Tflops单精度峰值性能)。图4和图5是CPU和GPU处理的亚距离单元全局最小熵距离对准结果对比。图6是将GPU距离对准数据输出，采用改进的ROPE算法进行相位补偿的成像结果图，并与基于CPU端的成像结果对比，肉眼难以区分差异，利用图像对比度准则来判断图像质量。图像对比度公式如下：

$C=\frac{\mathrm{\σ}\[{\left|I(q,k)\right|}^2\]}{E\[{\left|I(q,k)\right|}^2\]}=\frac{\sqrt{E\{\[{\left|I(q,k)\right|}^2\]-E{\[\left|I(q,k)\right|\]}^2\]\}}}{E\[{\left|I(q,k)\right|}^2\]}---\left(5\right)$

式中，σ(·)表示方差，E(·)表示均值，|I(q,k)|²表示图像各店像素强度。计算得到基于CPU图像对比度为5.1914，基于GPU图像对比度为5.1964。从图6和图7可见，GPU与CPU成像质量相当。

GPU和CPU处理时间对比(不包含矩阵分配和其它初始化)，如表2所示：

表2

可以看出，基于CUDA平台的ISAR距离对准获得了理想的加速比和实时处理效率。

以上所述仅是本发明的优选实施方式，应当指出：对于本技术领域的普通技术人员来说，在不脱离本发明原理的前提下，还可以做出若干改进和润饰，这些改进和润饰也应视为本发明的保护范围。

Claims (2)

1.一种基于GPU的并行结构ISAR距离对准方法，其特征在于：包括顺序执行的以下步骤：

步骤一、CPU读取ISAR原始回波数据，大小为Na*Nr，Na为方位向采样点数，Nr为距离向采样点数；

步骤二、采用统一计算设备架构流水线异步并行机制进行如下的分流脉冲压缩过程：将步骤一中读取的原始回波数据按方位向分块后分配给不同的工作流处理，每个工作流单独完成将本块数据从CPU拷贝到GPU中并在GPU中进行距离向压缩的工作得到脉压数据；

步骤三、在GPU中将脉压数据并行进行全局最小熵距离对准；

步骤四、将完成距离对准的数据分流从GPU中拷贝到CPU；

步骤五、利用FreeImage库函数输出图像。

2.根据权利要求1所述的一种基于GPU的并行结构ISAR距离对准方法，其特征在于：所述步骤三包括如下过程：

步骤3.1、将脉压数据转置为Nr列Na行，然后变换到距离向频域补零，再逆傅里叶变换到距离向时域，得到距离向时域插值后的TT*Nr列Na行的脉压数据，其中TT是插值倍数；

步骤3.2、对每个线程块的x维度和y维度均分配TILE_DIM个线程，每个线程网格的x维度分配Na/TILE_DIM个线程块、y维度分配TT*Nr/TILE_DIM个线程块，TILE_DIM固定为16，编写取模值核函数求插值一维距离像包络；

步骤3.3、编写转置核函数，将插值一维距离像包络的矩阵转置为Na列TT*Nr行；

步骤3.4、利用共享内存和归约算法编写包络求和核函数对一维距离像包络按方位向求和，分两步计算：第一步，对每个线程网格x维度分配TT*Nr个线程块、y维度分配Na/threadnum个线程块，每个线程块x维度分配threadnum个线程、y维度分配1个线程，threadnum取值为64，求得TT*Nr*Na/threadnum个和值；第二步，对每个线程网格x维度分配TT*Nr个线程块、y维度分配1个线程块，每个线程块的x维度分配Na/threadnum个线程、y维度分配1个线程，对上一步中求得的TT*Nr*Na/threadnum个和值继续求和，得到TT*Nr个和元素，即为平均距离像ARP；

步骤3.5、对每个线程网格x维度分配TT*Nr/TILE_DIM个线程块、y维度分配1个线程块，对每个线程块x维度分配为TILE_DIM个线程、y维度分配1个线程，编写核函数求平均距离像ARP的对数

步骤3.6、对每个线程网格x维度分配TT*Nr/TILE_DIM个线程块、y维度分配Na/TILE_DIM个线程块，对每个线程块x维度和y维度均分配TILE_DIM个线程，编写与一维距离像互相关核函数，通过频域共轭相乘实现互相关运算；

步骤3.7、利用共享内存和归约算法编写核函数求互相关矩阵每行最大值，分两步计算：第一步，对每个线程网格x维度分配TT*Nr/threadnum个线程块、y维度分配Na个线程块，对每个线程块x维度分配threadnum个线程、y维度分配1个线程，求得Na*TT*Nr/threadnum个最大值及其列索引；第二步，对每个线程网格x维度分配Na个线程块、y维度分配1个线程块，对每个线程块x维度分配TT*Nr/threadnum个线程、y维度分配1个线程，对上一步中求得的Na行TT*Nr/threadnum列最大值矩阵求每行最大值及其列索引，从而求得Na个距离偏移量；

步骤3.8、重复3.2～3.7的迭代过程，当达到设定迭代次数后，利用最近一次迭代所得到的距离偏移量构造频移因子exp(-j2πm△_r(n)/(TTN_r))，m是距离单元标号，△_r(n)是第n个回波的距离偏移量；

步骤3.9、对每个线程网格x维度分配Na/TILE_DIM个线程块、y维度分配Nr/TILE_DIM个线程块，对每个线程块x维度和y维度均分配TILE_DIM个线程，编写未插值一维距离像亚距离单元对准核函数，每个线程索引独立完成未插值一维距离像第n个脉冲第m个距离门处元素值与频移因子exp(-j2πm△_r(n)/(TTN_r))的相乘运算，完成未插值一维距离像的亚距离单元对准；

步骤3.10、将亚距离单元对准后的未插值一维距离像距离向逆傅里叶变换，再转置为Na列Nr行的距离对准数据。