CN103630891A - 利用gpu实现外辐射源雷达中估计目标来波方向的方法 - Google Patents

Abstract

本发明公开了一种利用GPU实现外辐射源雷达中估计目标来波方向的方法。主要解决传统目标来波方向估计方法受限于信号处理平台和算法，难以满足实时性需求的问题。其技术方案是：基于压缩感知的思想，将进行时域干扰相消和距离-多普勒二维相关处理后的单元数据加载到基于GPU的计算系统中，并且构建感知矩阵，根据处理后的单元数据和感知矩阵，在GPU中利用并行的方法实现基于压缩感知的方位向稀疏重构求解，获得求解后的方位向信号矢量，利用该方位向信号矢量重构出目标的来波方向。本发明具有硬件平台搭建简单、稳定性好、运算速度快、扩展性强、易于开发的优点，可用于对外辐射源雷达中的多个目标的来波方向进行估计。

Description

利用GPU实现外辐射源雷达中估计目标来波方向的方法

技术领域

本发明属于雷达技术领域，涉及雷达信号处理方法，具体地说是一种利用GPU实现外辐射源雷达中估计目标来波方向的方法，用于对目标进行精确测向。

技术背景

近年来，利用广播、电视、卫星等民用照射源的无源雷达探测技术越来越受到大家的重视。无源雷达系统本身不需要发射信号，而是利用环境中已有的或者目标本身发射的电磁波信号进行目标探测与定位，因此它具有抗干扰、抗反辐射导弹、抗低空突防和反隐身的综合“四抗”潜力。

依据探测对象的不同，无源雷达主要分为两大类：一类是利用探测目标自身辐射进行探测和测向；另一类是基于外辐射源的无源雷达系统，利用目标反射环境中存在的电磁波信号来进行探测和测向。然而由于外辐射源雷达具有非协作、不可控、不可预知的特性，使得外辐射源雷达的探测性能往往较传统雷达系统差，因此工程中需要采用更复杂的信号处理手段，来有效的完成目标的定位与跟踪，但同时会产生计算量大，处理复杂的问题。

目前，在外辐射源雷达系统中对目标测向主要通过两种算法来实现：一种是多重信号分类MUSIC算法，另一种是旋转不变子空间EPSRIT算法。这两种算法通常需要多快拍积累，处理时间较长难以满足实时需求，同时对多信源和相关信源的估计效果较差以及对信噪比要求高且抗噪声性能差。由于在外辐射源雷达中获取的信号里直达波和多径能量很强，为了抑制干扰和提高目标回波的信噪比，需要对目标来波进行时域干扰相消和距离-多普勒二维相关处理，然而经过时域干扰相消和距离-多普勒二维相关处理以后，能够应用于测向的数据只有单次快拍，此时多重信号分类MUSIC等传统的子空间类算法将失效。同时，在对目标的测向实现过程中主要是通过开发专用信号处理板进行数据处理，但这种硬件平台开发费用相对较高，而且随着数据处理运算量的增加，一般只能依靠成倍增加硬件的规模来实现功能。此外，信号板的开发时间、调试周期较长，可操作性较差，而且随着硬件规模的成倍增加，开发难度变得更大和费用增加，无法满足科研实验与应用的实时性与节约成本等需求。

发明内容

本发明的目的在于克服上述已有技术的不足，提出一种利用GPU实现外辐射源雷达中估计目标来波方向的方法，以简化系统结构，降低开发费用，改进系统可操作性，提高数据处理速度。

为实现上述目的，本发明的实现步骤包括如下：

(1)获取用于估计目标来波方向的外辐射源雷达信号数据y：

(1a)将外辐射源雷达的回波天线各个阵元接收的信号经过数模转换，滤波，抽取变换到数字基带；

(1b)将变换到数字基带后的数据进行时域干扰相消和距离-多普勒二维相关处理，得到单元数据y，再将单元数据y传送到主机端；

(2)在主机端进行数据初始化：

(2a)在主机端为中央处理器CPU分配计算所需内存，为图形处理器GPU分配计算所需显存空间；

(2b)设阵元数为M≥8、方向重构数为N≥100、迭代终止判断值l=0.01；

(2c)设置步长初始值s₁=1、步长的更新因子β=0.5，惩罚因子的初始值t₀=1，惩罚因子的更新系数μ=2；

(2d)在主机端构造一个2M*2N的感知矩阵A与初始方位向信号矢量（x,u）₀ ^T，设（x,u）₀ ^T的值x＝[00...00]_2*N ^T，u＝[11...11]_N ^T，T表示向量的转置；

(2e)将单元数据y加载到中央处理器CPU内存中；

(3)将中央处理器CPU内存中的单元数据y和感知矩阵A以及初始方位向信号矢量（x,u）₀ ^T传输至图形处理器GPU的显存中；

(4)在图形处理器GPU中进行基于压缩感知的方位向稀疏重构求解：

(4a)从图形处理器GPU中获取单元数据y、感知矩阵A与初始方位向信号矢量（x,u）₀ ^T，并设置基于压缩感知的方位向稀疏重构求解迭代次数k的初始值为1；

(4b)计算对偶空隙η：

$\mathrm{\η}={\left|\right|A*x-y\left|\right|}_2^2+{\left|\right|x\left|\right|}_2^2+0.25*{\left|\right|2*(A*x-y)\left|\right|}_2^2+{[2*(A*x-y)]}^T*y,$ 式中，x为方位向信号矢量（x,u）_k-1 ^T中的向量，k为迭代次数；

(4c)判断迭代终止条件是否满足，若满足则跳到步骤(4)，若不满足，则继续下一步操作，终止条件为： $\mathrm{\η}/(-0.25*{\left|\right|2*(A*x-y)\left|\right|}_2^2-{[2*(A*x-y)]}^T*y)\≤l;$

(4d)更新第k次迭代的惩罚因子：t_k＝max(min(4*N*μ/η,μ*t_k-1),t_k-1)；

(4e)根据惩罚因子t_k，构建代价函数J：

$J=t_k{\left|\right|\mathrm{Ax}-y\left|\right|}_2^2+t_k{1}^{T}u+\underset{i=1}{\overset{N}{\mathrm{\Σ}}}\log ({u_i}^2-{u_i}^2-{x_{i+N}}^2),$

式中，u为方位向信号矢量（x,u）_k-1 ^T中的向量，u_i表示向量u的第i个元素，x_i与x_i+N分别表示向量x的第i个和第i+N个元素，1^T表示为1×N维的单位向量；

(4f)计算代价函数J的第k次迭代的共轭梯度 $g_k=\left[\begin{array}{c}{\▿}_{x}J\\ {\▿}_{u}J\end{array}\right]$ 和海森Hessian矩阵H_k＝▽²J，其中▽_xJ表示代价函数J对向量x求一阶偏导，▽_uJ表示代价函数J对向量u求一阶偏导，▽²J表示代价函数J对信号矢量（x,u）_k-1 ^T的二阶导数；

(4g)求出第k次迭代的搜索方向：(Δx,Δu)_k ^T＝-H_k ^-1g_k，其中H_k ^-1为H_k的逆矩阵；

(4h)更新第k次迭代的方位向信号矢量，得到更新后的方位向信号矢量：(x,u)_k ^T＝(x,u)_k-1 ^T+s_k(Δx,Δu)_k ^T；

(4i)把k的值增加1，然后更新回归线性搜索方法的步长：s_k＝s_k-1*β，再返回到(4b)；

(5)将更新后的方位向信号矢量（x,u）_k-1 ^T从图形处理器GPU显存中传送到中央处理器CPU的内存中并保存，再根据该更新后的方位向信号矢量（x,u）_k-1 ^T重构出目标的来波方向；

(6)在主机端释放中央处理器CPU中所分配的内存与图形处理器GPU中所分配的显存。

本发明具有以下优点：

1)本发明由于利用了CPU的逻辑控制能力和GPU的大规模数据并行能力，使得浮点运算能力大，数据扩展性强，能有效提高信号处理速度。

2)本发明由于采用了基于GPU的高性能计算系统，与一些专业的信号处理板相比，开发费用少，开发周期显著降低。

3)本发明采用基于压缩感知的思想，利用稀疏重构求解的方法来获得目标的来波方向，提高了外辐射源雷达中多目标的测向精度。

附图说明

图1是本发明的实现流程图；

图2是本发明的雷达接收系统框图；

图3是用本发明对两个目标进行来波方向估计后的结果图；

图4是用本发明对三个目标进行来波方向估计后的结果图；

图5是用本发明和现有MATLAB处理方法估计两个目标的来波方向所需要的处理时间对比图；

图6是用本发明和现有MATLAB处理方法估计三个目标的来波方向所需要的处理时间对比图。

具体实施方式

参照图1，本发明利用GPU实现外辐射源雷达中估计目标来波方向的方法，其实现步骤如下：

步骤1：获取用于估计目标来波方向的外辐射源雷达信号数据y。

参照图2，本发明应用调频广播FM信号作为外辐射源雷达的机会照射源信号，以一个M阵元的均匀线阵作为雷达回波天线，接收由运动目标反射的回波信号，同时接收FM基站直达波信号、多径等干扰噪声信号。为了检测到目标回波和估计目标的来波方向，需要消除FM基站直达波信号、多径等干扰噪声信号。因此将天线各个阵元接收的信号经过数模转换，滤波，抽取等操作变换到数字基带以后，再通过时域干扰相消和距离-多普勒二维相关处理消除各个阵元接收的直达波和多径等干扰噪声信号，得到相消以后的单元数据y，将该单元数据y传送到主机端。

步骤2：在主机端进行数据初始化。

(2a)在主机端为中央处理器CPU分配计算所需内存，为图形处理器GPU分配计算所需显存空间；

(2b)设阵元数为M≥8、方向重构数为N≥100、迭代终止判断值l=0.01；

(2c)设置步长初始值s₁=1、步长的更新因子β=0.5，惩罚因子的初始值t₀=1，惩罚因子的更新系数μ=2；

(2d)在主机端设置一个2M*2N的感知矩阵A和初始方位向信号矢量（x,u）₀ ^T，其中（x,u）₀ ^T的x＝[00...00]_2*N ^T，表示方位向稀疏重构求解的初始位置，u＝[11...11]_N ^T是一个用于方位向稀疏重构求解的辅助向量，T表示向量的转置；

(2e)将单元数据y加载到中央处理器CPU内存中。

步骤3：选择能支持计算统一设备架构CUDA的图形处理器GPU，并调用计算统一设备架构CUDA中的cudaSetDevice()函数将该图形处理器GPU与方位向稀疏重构算法线程进行连接。

步骤4：在主机端调用计算统一设备架构CUDA中的同步数据拷贝函数cudaMemcpy()把CPU内存中的单元数据y和感知矩阵A以及初始方位向信号矢量（x,u）₀ ^T传输至图形处理器GPU的显存中。

步骤5：在图形处理器GPU中进行基于压缩感知的方位向稀疏重构求解：

(5a)从图形处理器GPU的显存中获取单元数据y、感知矩阵A与初始方位向信号矢量（x,u）₀ ^T，并设置基于压缩感知的方位向稀疏重构求解迭代次数k的初始值为1；

(5b)根据单元数据y、感知矩阵A与方位向信号矢量（x,u）_k-1 ^T，计算对偶空隙η： $\mathrm{\η}={\left|\right|A*x-y\left|\right|}_2^2+{\left|\right|x\left|\right|}_2^2+0.25*{\left|\right|2*(A*x-y)\left|\right|}_2^2+{[2*(A*x-y)]}^T*y,$ 式中，x为方位向信号矢量（x,u）_k-1 ^T中的向量，k为迭代次数，该对偶空隙η用于判断方位向信号矢量（x,u）_k-1 ^T是否为正确的解；

(5c)将单元数据y、感知矩阵A、方位向信号矢量（x,u）_k-1 ^T以及对偶空隙η代入迭代终止条件： $\mathrm{\η}/(-0.25*{\left|\right|2*(A*x-y)\left|\right|}_2^2-{[2*(A*x-y)]}^T*y)\≤l,$ 判断是否满足迭代终止条件，若满足则跳到步骤6，若不满足则继续下一步操作；

(5d)根据对偶空隙η更新第k次迭代的惩罚因子：

t_k＝max(min(4*N*μ/η,μ*t_k-1),t_k-1)；

(5e)根据单元数据y、感知矩阵A、方位向信号矢量（x,u）_k-1 ^T与更新后的惩罚因子t_k构建代价函数： $J=t_k{\left|\right|\mathrm{Ax}-y\left|\right|}_2^2+t_k{1}^{T}u+\underset{i=1}{\overset{N}{\mathrm{\Σ}}}\log ({u_i}^2-{u_i}^2-{x_{i+N}}^2),$

式中，u为方位向信号矢量（x,u）_k-1 ^T中的向量，u_i表示向量u的第i个元素，x_i与x_i+N分别表示向量x的第i个和第i+N个元素，1^T表示为1×N维的单位向量；

(5f)计算代价函数J的第k次迭代的共轭梯度 $g_k=\left[\begin{array}{c}{\▿}_{x}J\\ {\▿}_{u}J\end{array}\right]$ 和海森Hessian矩阵H_k＝▽²J，其中▽_xJ表示代价函数J对向量x求一阶偏导，▽_uJ表示代价函数J对向量u求一阶偏导，▽²J表示代价函数J对信号矢量（x,u）_k-1 ^T的二阶导数；

(5g)根据共轭梯度g_k和海森Hessian矩阵H_k求出第k次迭代的搜索方向：(Δx,Δu)_k ^T＝-H_k ^-1g_k，其中H_k ^-1为H_k的逆矩阵，其实现步骤如下：

(5g1)在图形处理器GPU中构造一个3N*3N的单位矩阵I，将海森Hessian矩阵H_k和单位矩阵I的所有元素加载到显存中，用于对海森Hessian矩阵H_k求逆，其中N为方向重构数；

(5g2)将用于对海森Hessian矩阵H_k求逆的线程格grid的大小设为((3N+15)/16，(3N+15)/16)，将线程格grid中的线程块block大小设为(16,16)，通过公式H_k*H_k ^-1＝I对海森Hessian矩阵H_k求逆，得到H_k的逆矩阵

(5g3)根据所得到海森Hessian矩阵H_k的逆矩阵

求出第k次迭代的搜索方向：(Δx,Δu)_k ^T＝-H_k ^-1g_k，其中g_k表示代价函数J的第k次迭代的共轭梯度；

(5h)根据搜索方向(Δx,Δu)_k ^T更新第k次迭代的方位向信号矢量，得到更新后的方位向信号矢量：(x,u)_k ^T＝(x,u)_k-1 ^T+s_k(Δx,Δu)_k ^T；

(5i)把k的值增加1，然后更新回归线性搜索方法的步长：s_k＝s_k-1*β，再返回到(5b)。

步骤6：调用统一计算设备架构CUDA中的同步数据拷贝函数cudaMemcpy()，将基于压缩感知的方位向稀疏重构求解得到更新后的方位向信号矢量(x,u)_k-1 ^T从图形处理器GPU显存中传送到中央处理器CPU的内存中并保存，再根据更新后的方位向信号矢量(x,u)_k-1 ^T重构出目标的来波方向。

步骤7：在主机端上调用Windows应用编程接口API中的free()函数和统一设备计算架构CUDA中的cudaFree()函数，释放为CPU所分配的内存和为GPU所分配的显存资源。

本发明的效果可通过以下实验结果进一步说明：

1）实验条件：

在实验中，硬件平台选用HP Z820工作站，图形处理器GPU型号为NVIDA TelsaC2075，中央处理器CPU型号为Intel Xeon，Win7系统，软件平台为Visual Stdio2008+CUDA4.0和MATLAB2009b，实验采用的数据是外辐射源雷达仿真数据中的一个距离-多普勒单元数据，该数据为经过时域干扰相消和距离-多普勒二维相关处理后的单元数据。

2）实验内容及效果：

仿真1，采用本发明提出的方法，在Visual Stdio2008+CUDA4.0的软件平台上，采用外辐射源雷达仿真数据中经过时域干扰相消和距离-多普勒二维相关处理后的单元数据，该数据中含有两个目标的来波方向信号，两个目标的来波方向分别为-0.1°和-5.8°。该回波天线阵列的阵元数M=10，设方向重构数N=180，运用本发明的方法得到两个目标的来波方向如图3所示。

从图3中可以看出，本发明的方法能够对两个目标进行有效分辨，并且测量得到的两个目标的方向分别为0°和-6°，与真实目标来波方向基本符合。

仿真2，在与仿真1相同条件下，对含有三个目标的来波方向信号的数据进行来波方向估计，这三个目标的来波方向分别为-6.8°、-0.1°和5.6°。该回波天线阵列的阵元数M=15，设方向重构数N=180，运用本发明的方法得到两个目标的来波方向如图4所示。

从图4中可以看出，本发明的方法能够对三个目标进行有效分辨，并且测量得到的三个目标的方向分别为-7°、0°和6°。需要指出是本试验中我们设定正确估计出目标的来波方向是指能够将目标的进行有效分辨，并且所求出的目标来波方向与目标真实的来波方向之间的差值小于1°。

仿真3，采用本发明提出的方法，在Visual Stdio2008+CUDA4.0的软件平台上，对含有两个目标的来波方向信号的数据进行来波方向估计，该回波天线阵列的阵元数M=10，设方向重构数N=180，在相同条件下，分别进行15次重复试验，记录每次试验的运行时间，并与现有MATLAB处理方法估计两个目标的来波方向的时间进行对比，得到了如图5所示的时间对比图。

仿真4，在与仿真3相同的条件下，对含有三个目标的来波方向信号的数据进行来波方向估计，该回波天线阵列的阵元数M=15，设方向重构数N=180，在相同条件下，分别进行15次重复试验，记录每次试验的运行时间，并与现有MATLAB处理方法估计三个目标的来波方向的时间进行对比，得到了如图6所示的时间对比图。

由图5、图6可知，使用本发明提出的方法处理时间均显著小于现有MATLAB处理估计目标来波方向方法的时间，平均加速比为21倍，而且处理时间稳定，验证了本发明提出的方法具有良好的稳定性和加速性能。

Claims (2)

1.一种利用GPU实现外辐射源雷达中估计目标来波方向的方法，包括如下步骤：

(1)获取用于估计目标来波方向的外辐射源雷达信号数据y：

(1a)将外辐射源雷达的回波天线各个阵元接收的信号经过数模转换，滤波，抽取变换到数字基带；

(1b)将变换到数字基带后的数据进行时域干扰相消和距离-多普勒二维相关处理，得到单元数据y，再将单元数据y传送到主机端；

(2)在主机端进行数据初始化：

(2a)在主机端为中央处理器CPU分配计算所需内存，为图形处理器GPU分配计算所需显存空间；

(2b)设阵元数为M≥8、方向重构数为N≥100、迭代终止判断值l=0.01；

(2c)设置步长初始值s₁=1、步长的更新因子β=0.5，惩罚因子的初始值t₀=1，惩罚因子的更新系数μ=2；

(2d)在主机端构造一个2M*2N的感知矩阵A与初始方位向信号矢量（x,u）₀ ^T，设（x,u）₀ ^T的值x＝[00...00]_2*N ^T，u＝[11...11]_N ^T，T表示向量的转置；

(2e)将单元数据y加载到中央处理器CPU内存中；

(3)将中央处理器CPU内存中的单元数据y和感知矩阵A以及初始方位向信号矢量（x,u）₀ ^T传输至图形处理器GPU的显存中；

(4)在图形处理器GPU中进行基于压缩感知的方位向稀疏重构求解：

(4a)从图形处理器GPU中获取单元数据y、感知矩阵A与初始方位向信号矢量（x,u）₀ ^T，并设置基于压缩感知的方位向稀疏重构求解迭代次数k的初始值为1；

(4b)计算对偶空隙η：

$\mathrm{\η}={\left|\right|A*x-y\left|\right|}_2^2+{\left|\right|x\left|\right|}_2^2+0.25*{\left|\right|2*(A*x-y)\left|\right|}_2^2+{[2*(A*x-y)]}^T*y,$

式中，x为方位向信号矢量（x,u）_k-1 ^T中的向量，k为迭代次数；

(4c)判断迭代终止条件是否满足，若满足则跳到步骤(4)，若不满足，则继续下一步操作，终止条件为： $\mathrm{\η}/(-0.25*{\left|\right|2*(A*x-y)\left|\right|}_2^2-{[2*(A*x-y)]}^T*y)\≤l;$

(4d)更新第k次迭代的惩罚因子：t_k＝max(min(4*N*μ/η,μ*t_k-1),t_k-1)；

(4e)根据惩罚因子t_k，构建代价函数J：

$J=t_k{\left|\right|\mathrm{Ax}-y\left|\right|}_2^2+t_k{1}^{T}u+\underset{i=1}{\overset{N}{\mathrm{\Σ}}}\log ({u_i}^2-{u_i}^2-{x_{i+N}}^2),$

式中，u为方位向信号矢量（x,u）_k-1 ^T中的向量，u_i表示向量u的第i个元素，x_i与x_i+N分别表示向量x的第i个和第i+N个元素，1^T表示为1×N维的单位向量；

(4f)计算代价函数J的第k次迭代的共轭梯度 $g_k=\left[\begin{array}{c}{\▿}_{x}J\\ {\▿}_{u}J\end{array}\right]$ 和海森Hessian矩阵H_k＝▽²J，其中▽_xJ表示代价函数J对向量x求一阶偏导，▽_uJ表示代价函数J对向量u求一阶偏导，▽²J表示代价函数J对信号矢量（x,u）_k-1 ^T的二阶导数；

(4g)求出第k次迭代的搜索方向：(Δx,Δu)_k ^T＝-H_k ^-1g_k，其中H_k ^-1为H_k的逆矩阵；

(4h)更新第k次迭代的方位向信号矢量，得到更新后的方位向信号矢量：(x,u)_k ^T＝(x,u)_k-1 ^T+s_k(Δx,Δu)_k ^T；

(4i)把k的值增加1，然后更新回归线性搜索方法的步长：s_k＝s_k-1*β，再返回到(4b)；

(5)将更新后的方位向信号矢量（x,u）_k-1 ^T从图形处理器GPU显存中传送到中央处理器CPU的内存中并保存，再根据该更新后的方位向信号矢量（x,u）_k-1 ^T重构出目标的来波方向；

(6)在主机端释放中央处理器CPU中所分配的内存与图形处理器GPU中所分配的显存。

2.根据权利要求1所述的利用GPU实现外辐射源雷达中估计目标来波方向的方法，其特征在于所述步骤(4g)中求出第k次迭代的搜索方向，按如下步骤进行：

(4g1)在图形处理器GPU中构造一个3N*3N的单位矩阵I，将海森Hessian矩阵H_k和单位矩阵I的所有元素加载到显存中，用于对海森Hessian矩阵H_k求逆，其中N为方向重构数；

(4g2)将用于对海森Hessian矩阵H_k求逆的线程格grid的大小设为((3N+15)/16，(3N+15)/16)，将线程格grid中的线程块block大小设为(16,16)，通过公式H_k*H_k ^-1＝I对海森Hessian矩阵H_k求逆，得到H_k的逆矩阵

(4g3)根据所得到海森Hessian矩阵H_k的逆矩阵

求出第k次迭代的搜索方向：(Δx,Δu)_k ^T＝-H_k ^-1g_k，其中g_k表示代价函数J的第k次迭代的共轭梯度。

Images