CN104237859A - 利用gpu实现外辐射源雷达多通道时域杂波抑制的方法 - Google Patents

Abstract

本发明涉及一种利用GPU实现外辐射源雷达多通道时域杂波抑制的方法，主要解决雷达探测环境中的杂波抑制问题。主要步骤为：首先获取校准后参考和监测通道数据，并分配计算所需内存、显存；然后加载数据到GPU显存，采用迭代算法计算自相关矩阵R_x，求逆得利用参考和监测通道数据计算互相关矩阵C，计算与C乘积得到矩阵D，计算参考通道数据与矩阵D乘积得到监测通道中直达波与多径分量从监测通道数据中减去得到目标回波分量；最后释放已分配的内存、显存。本发明克服了迭代类杂波抑制算法收敛速度和稳定性等难以确定的困难，降低了ECA及ECA-B算法的时间和空间复杂度，并具有平台搭建简单、稳定性好、运算速度快、扩展性强、易于开发等优点。

Description

利用GPU实现外辐射源雷达多通道时域杂波抑制的方法

技术领域

本发明涉及外辐射源雷达技术领域，尤其涉及一种利用GPU实现外辐射源雷达多通道时域杂波抑制的方法。

背景技术

近年来，利用广播、电视、卫星等民用照射源的外辐射源雷达探测技术越来越受到大家的重视。外辐射源雷达系统本身不需要发射信号，直接利用环境中已有的或者目标本身发射的电磁波信号进行目标探测与定位，因此具有抗干扰、抗反辐射导弹、抗低空突防和反隐身的综合“四抗”潜力。然而由于外辐射源雷达具有非协作、不可控、不可预知的特性，使得外辐射源雷达的探测性能往往较传统雷达系统差，特别是接收到的目标回波幅度弱，易被直达波、多径杂波及干扰淹没，因而在目标检测前需要进行杂波抑制处理，其中杂波抑制的实时处理是外辐射源雷达实时化关键步骤之一。

常用的如LMS，NLMS，BLMS，RLS，LSL等闭环自适应滤波算法属于迭代类算法，处理时数据按一定批次输入自适应滤波器，由于每批次数据量较少，算法时间和空间复杂度相对较低，易于在FPGA、DSP或CPU上实现。但迭代类算法面临收敛速度、稳定性和适应能力等难以确定的问题。F.Colone等人提出基于最小二乘(Least Square,LS)算法的扩展杂波相消(Extensive Cancellation Algorithm,ECA)和扩展相消批处理算法(Extensive Cancellation Algorithm Batches,ECA-B)，算法直接利用输入信号和期待响应的信息计算自适应权值，无需引入输出信息的反馈，避免了迭代类算法的收敛问题，因而不会引起输出误差的发散，是一种稳健的自适应算法。但由于其求解过程中需要计算高维复矩阵乘法和矩阵求逆，耗时较长；且为减少估计误差，通常需要提高样本数目以取得较好的滤波效果，导致空间复杂性增加。原理上，ECA和ECA-B算法针对单通道数据处理，以ECA为例，采用常规算法，对单个通道处理时，时间和空间复杂度分别为O(NK²+K³)和O(NK)，其中N为数据点数，K为杂波子空间自由度。当通道数为M时，时间复杂度和空间复杂度分别为O(MNK²+MK³)和O(MNK)，算法存储空间占用和计算量与通道数目成正比，若采用传统基于FPGA、DSP或ARM等微处理器的实时处理方案，不仅开发周期、费用相对较高，而且存储空间和硬件计算能力需求难以满足。

近年来，随着大规模集成电路技术的发展，中央处理器(CentralProcessing Unit,CPU)和图形处理器(Graphic Processing Unit,GPU)性能与日俱增，为外辐射源雷达的实时信号处理提供了新的解决方案。相对于传统方案，其具有成本低廉、结构简单、开发容易等优点。由于GPU能够并行计算，且具有存储空间充足、数据吞吐量大、浮点运行能力强等优点，ECA和ECA-B算法(以下并称为ECA类算法)更适合利用GPU予以实现。但该类算法处理多通道数据时需要分别对各个单通道进行处理，计算量和空间占用存在大量冗余。

发明内容

本发明主要是解决现有技术所存在的技术问题；提供了一种利用了CPU强大的逻辑控制能力和GPU的大规模数据并行运算能力，使得浮点运算能力强，数据扩展性好，能有效提高信号处理速度的外辐射源雷达多通道时域杂波抑制的方法。

本发明还有一目的是解决现有技术所存在的技术问题；提供了一种开发费用少，开发周期缩短，调试简单的利用GPU实现外辐射源雷达多通道时域杂波抑制的方法。

本发明再有一目的是解决现有技术所存在的技术问题；提供了一种大大减小了计算量和存储空间占用，缩短了信号处理耗时的利用GPU实现外辐射源雷达多通道时域杂波抑制的方法。

本发明的上述技术问题主要是通过下述技术方案得以解决的：

一种利用GPU实现外辐射源雷达多通道时域杂波抑制的方法，其特征在于，包括以下步骤：

步骤1、获取用于杂波抑制处理的外辐射源雷达信号数据：具体是首先将外辐射源雷达的接收阵列各个阵元接收到的回波信号经过模数转换、滤波抽取变换到数字基带，然后将变换到数字基带后的数据传送至主机端，在主机端对各个通道数据乘以相应通道校准值进行校准；

步骤2、在主机端将步骤1接收的数据进行初始化：主机端使用函数malloc为中央处理器CPU分配计算所需内存，使用函数cudaMalloc为图形处理器GPU分配计算所需显存空间；

步骤3、主机端将接收的数据传送至图形处理器GPU显存中，并定义：

监测通道数为M,杂波抑制距离元长度为K，抑制数据长度为N；

参考通道数据为s_ref＝[s_ref(0) s_ref(1) … s_ref(N+K-2)]^T；

监测通道数据为S_sur＝[s_sur,0 s_sur,1 … s_sur,M-1]；

其中，每个监测通道数据为s_sur,m＝[s_sur(K-1) s_sur(K) …s_sur(N+K-2)]^T，m为监测通道编号，T代表转置；

步骤4、在图形处理器GPU中对各通道数据进行并行杂波抑制处理，即杂波抑制结果 ${\tilde{S}}_{\sup }=S_{\mathrm{sur}}-{\tilde{S}}_{\mathrm{sur}}=S_{\mathrm{sur}}-X_{\mathrm{ref}}{\left(X_{\mathrm{ref}}^H{X}_{\mathrm{ref}}\right)}^{-1}{X}_{\mathrm{ref}}^H{S}_{\mathrm{sur}};$ 包括以下子步骤：

步骤4.1、利用参考通道数据采用迭代算法计算自相关矩阵 $R_x=X_{\mathrm{ref}}^H{X}_{\mathrm{ref}},$ 具体是：

步骤4.11、首先计算矩阵R_x第i＝0行值，即R_x(0,j)，其中0≤j≤K-1；

步骤4.12、计算出第i＝i+1行的值；若列序号j<i,则否则通过利用关系式：

$R_x(i,j)-R_x(i-1,j-1)=s_{\mathrm{ref}}^{*}(N+i-1)s_{\mathrm{ref}}(N+j-1)-s_{\mathrm{ref}}^{*}(i-1)s_{\mathrm{ref}}(j-1)$ 及对角线元素R_x(i-1,j-1)计算R_x(i,j)；

步骤4.13、重复步骤4.12，若i＝K-1则计算完毕后停止；

步骤4.2、对自相关矩阵R_x进行原地求逆得到

步骤4.3、利用参考通道数据和监测通道数据进行互相关运算得到 $C=X_{\mathrm{ref}}^H{S}_{\mathrm{sur}};$

步骤4.4、计算及互相关矩阵C的乘积得到矩阵

步骤4.5、计算参考通道数据和D矩阵乘积得到监测通道中直达波和多径杂波分量 ${\tilde{S}}_{\mathrm{sur}}=X_{\mathrm{ref}}D;$

步骤4.6、从原始监测通道数据减去杂波分量得到剩余回波分量，即杂波抑制结果 ${\tilde{S}}_{\sup }=S_{\mathrm{sur}}-{\tilde{S}}_{\mathrm{sur}};$

其中，s(i)＝|s_ref(i) s_ref(i+1) … s_ref(N+i-1)|^T，X_ref＝|s(0) s(1) … s(K-1)|；

步骤5、将杂波抑制结果在图形处理器GPU中继续进行其他处理或传送至中央处理器CPU的内存中并保存；

步骤6、在主机端使用函数free释放中央处理器CPU所分配的内存，使用cudaFree释放图形处理器GPU中所分配的显存。

在上述的利用GPU实现外辐射源雷达多通道时域杂波抑制的方法，步骤4.3与步骤4.5计算C矩阵和矩阵时直接在GPU核函数中利用s_ref进行计算；

若占用显存量的数据格式为单精度浮点复数，由于直接使用占用显存量为8(N+K-1)的向量s_ref，不需要形成占用显存量为8NK的矩阵X_ref，大大降低显存需求。

若占用显存量的数据格式为双精度浮点复数，由于直接使用占用显存量为16(N+K-1)的向量s_ref，不需要形成占用显存量为16NK的矩阵X_ref，大大降低显存需求；

本算法大大降低显存需求，跟数据格式没有关系。

在上述的利用GPU实现外辐射源雷达多通道时域杂波抑制的方法，所述步骤4.11计算矩阵R_x第i＝0行值，是利用s_ref来计算矩阵R_x第i＝0行值，核函数每个线程块计算num个点，num取值可为1,2或4，共分为K/num个线程块，每个线程块为m个线程，每个线程分配m个寄存器rA，1个寄存器rB,m个寄存器rC；每个线程块分配m+num-1个共享存储器sA；其中，计算第i个线程块包括以下子步骤：

步骤4.111、将结果寄存器数组rC清零，序号j赋值为0；

步骤4.112、从向量s_ref偏移i*num+j处加载(m+num-1)个数据到共享存储器sA中；加载sA前num个数据到寄存器rA中；从向量s_ref偏移j处加载m个数据到线程块各自寄存器rB中；

步骤4.113、计算rA数组共轭与rB乘积结果，并累加到数组rC中

步骤4.114、给j增加m；若j<N,重复执行步骤4.112,4.113,4.114，否则进入步骤4.115.

步骤4.115、依次将rC寄存器数组中的数据加载到共享存储器sA中，进行并行规约求和，并将求和结果保存在R_x矩阵第i*num到第(i+1)*num-1处。

在上述的利用GPU实现外辐射源雷达多通道时域杂波抑制的方法，所述步骤4.4中，计算互相关矩阵C是利用s_ref和S_sur进行计算，核函数分成K/2个线程块，每个线程块包含K个线程；每个线程分配长度为1的寄存器rA，rA0，长度为M的寄存器数组rC；每个线程块分配KM长度的共享存储器sC；每个线程块计算C矩阵2行元素，计算第i和第i+1行元素的线程块按如下步骤进行：

4.31、从s_ref偏移i和i+1处加载2个元素u和v至寄存器rA,rA0中，从S_sur中加载第一行共M个元素pv到GPU共享存储器sC中，M为S_sur矩阵列数；

4.32、分别计算u、v与pv对应元素乘积，保存结果1行存放在寄存器数组rC中，另1行保存在GPU共享存储器sC中；

4.33、继续加载s_ref下2个元素继续保存至寄存器rA,rA0中，从S_sur中加载下一行共M个元素pv到GPU共享存储器sC中，重复步骤4.32，并将计算结果第1行对寄存器数组中结果进行乘加运算，计算结果第2行对GPU共享存储器进行乘加运算；

4.34、重复步骤4.33，直至把S_sur全部N行加载计算完毕后，把寄存器数组rC中结果保存至C矩阵第i行，把共享存储器sC中结果保存至C矩阵第i+1行，计算完毕。

在上述的利用GPU实现外辐射源雷达多通道时域杂波抑制的方法，所述步骤4.5中，利用s_ref和D计算直达波和多径杂波分量核函数分成二维网格(N/K,32)规模的线程块，每个线程块划分成二维线程(32，K/32)，用于计算矩阵的K/32行，每行共M个元素，也即监测通道数目；每个线程分配长度为M的寄存器数组rC；每个线程块分配长度为(K+K/32-1)的共享存储器sA；计算第i个线程块包括以下子步骤：

4.51、将结果寄存器数组rC清零，从向量s_ref偏移处i*(32N/K)处加载K+K/32-1个元素到共享存储器sA中；

4.52、分别以sA[0]，sA[1]，…，sA[K/32-1]为起点计算与D矩阵第j＝0列乘积，并累加到rC寄存器数组中；

4.53、重复步骤1，取D矩阵第j＝1,2，…，M-1列乘积，并累加到rC寄存器数组中；

4.54、保存rC中结果到矩阵第i*K/32到第(i+1)*K/32-1行。

因此，本发明具有如下优点：1、本发明利用了CPU强大的逻辑控制能力和GPU的大规模数据并行运算能力，使得浮点运算能力强，数据扩展性好，能有效提高信号处理速度；2、本发明由于采用了GPU的高性能计算系统，与传统基于嵌入式信号处理开发板的方案相比，开发费用少，开发周期缩短，调试简单；3、本发明利用了向量列扩展的思想，对原ECA类算法进行扩展，使其能够同时对多通道或多波束数据进行处理，且提出迭代计算自相关矩阵的方法，大大减小了计算量和存储空间占用，缩短了信号处理耗时。

附图说明

图1为本发明方法中杂波抑制处理流程图。

图2为本发明方法计算流程示意图。

图3为本发明方法计算自相关矩阵迭代算法框图。

图4为本发明方法计算矩阵原地求逆算法框图。

图5为本发明方法GPU端kernel实现总框图。

图6为本发明方法中涉及计算互相关矩阵C计算框图。

图7为本发明实施仿真例程中杂波抑制前通道1距离多普勒谱。

图8为本发明实施仿真例程中杂波抑制后通道1距离多普勒谱。

图9为本发明实施仿真例程中杂波抑制CPU端和GPU处理绝对误差图。

图10为本发明实施实测例程中杂波抑制前通道1距离多普勒谱。

图11为本发明实施实测例程中杂波抑制后通道1距离多普勒谱。

具体实施方式

下面通过实施例，并结合附图，对本发明的技术方案作进一步具体的说明。

实施例：

下面以具体实施例结合附图对本发明作进一步说明。

图1为本发明利用GPU实现外辐射源雷达多通道时域杂波抑制的方法流程图。具体步骤如下：

步骤1：获取用于杂波抑制处理的外辐射源雷达信号数据：(1a)将外辐射源雷达的接收阵列各个阵元接收到得回波信号经过模数转换、滤波抽取变换到数字基带；(1b)将变换到数字基带后的数据传送至主机端，然后进行校准。

步骤2：在主机端进行数据初始化：(2a)为中央处理器CPU分配计算所需内存，为图形处理器GPU分配计算所需显存空间，初始化统一设备架构CUDA线性计算代数库CUBLAS；(2b)设监测通道数为M,杂波抑制距离元长度为K，抑制数据长度为N。

步骤3：将主机端参考和监测通道数据传送至图形处理器GPU显存中，其中参考通道s_ref＝[s_ref(0) s_ref(1) …s_ref(N+K-2)]^T，监测通道S_sur＝[s_sur,0 s_sur,1 … s_sur,M-1]，其中每个监测通道数据为s_sur,m＝[s_sur(K-1) s_sur(K-1) … s_sur(N+K-2)]^T，m为监测通道编号，T代表转置。

步骤4：在图形处理器GPU中对各通道数据进行并行杂波抑制处理，得到杂波抑制结果 ${\tilde{S}}_{\sup }=S_{\mathrm{sur}}-{\tilde{S}}_{\mathrm{sur}}=S_{\mathrm{sur}}-X_{\mathrm{ref}}{\left(X_{\mathrm{ref}}^H{X}_{\mathrm{ref}}\right)}^{-1}{X}_{\mathrm{ref}}^H{S}_{\mathrm{sur}},$ 其中X_ref＝|s(0) s(1) … s(K-1)|，s(i)＝|s_ref(i) s_ref(i+1) … s_ref(N+i-1)|^T。参照图2，按以下步骤进行计算：(4a)按照图3所示，利用参考通道数据采用迭代算法计算自相关矩阵(4b)对自相关矩阵进行原地求逆得到如图4所示；(4c)利用参考通道数据和监测通道数据进行互相关运算得到(4d)计算及互相关矩阵C的乘积得到矩阵(4e)计算参考通道数据和D矩阵乘积得到监测通道中直达波和多径杂波分量(4f)从原始监测通道数据减去杂波分量得到剩余回波分量，即杂波抑制结果 ${\tilde{S}}_{\sup }=S_{\mathrm{sur}}-{\tilde{S}}_{\mathrm{sur}}.$

步骤5：将杂波抑制结果在图形处理器GPU中继续进行其他处理或传送至中央处理器CPU的内存中并保存。

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

图5为该发明方法全部GPU实现流程图。图6为计算步骤(4c)GPU实现框图。参考通道数据s_ref和多个监测通道数据S_sur均保存在GPU全局显存中。核函数每个线程计算C矩阵2行元素，以计算第i和第i+1行为例说明，计算按如下步骤进行：首先从s_ref偏移i和i+1处加载2个元素u和v至寄存器中，从S_sur中加载第一行共M个元素pv到GPU共享存储器中，然后分别计算u、v与pv对应元素乘积，保存结果1行存放在寄存器数组中，另1行保存在GPU共享存储器中。然后继续加载s_ref下2个元素继续保存至寄存器中，从S_sur中加载下一行共M个元素pv到GPU共享存储器中，再计算u、v与pv对应元素乘积，计算结果第1行对寄存器数组中结果进行乘加运算，计算结果第2行对GPU共享存储器进行乘加运算。依次类推，直至把S_sur全部N行加载计算完毕后，把寄存器数组中结果保存至C矩阵第i行，把共享存储器中结果保存至C矩阵第i+1行，计算完毕。

步骤(4a)中迭代计算自相关矩阵的第一行或第一列元素及步骤(4e)计算所采用方法与此类似。此方法能够最大化利用GPU中带宽较大的寄存器资源和共享存储器资源，重用了数据，减少了全局显存的访问次数，从而降低算法运行时间。

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

1)实验软硬件配置如表1所示

表1软硬件配置情况

外辐射源雷达仿真数据为采样率为312.5KHz的8路FM信号，1路为参考通道，其余7路为监测通道，数据保存为单精度浮点类型。仿真时各个监测通道中多径和目标参数设置相同，数据详细参数如表2所示：

表2仿真信号详细参数

实验采用C语言实现本发明提出的多通道时域杂波抑制算法，分为CPU串行版本和GPU并行版本。为了对比算法的时效性，同时测试了未采用迭代算法而直接利用矩阵相乘计算自相关矩阵的常规计算方案。

2)实验内容及效果：

图7为仿真信号杂波抑制前通道1距离多普勒谱，目标完全被直达波及多径杂波旁瓣淹没。图8为仿真信号杂波抑制后通道1距离多普勒谱，仿真加入的两个目标清晰可见，准确的得到目标距离和多普勒频移信息。图9为本发明所提串行算法和并行算法计算绝对误差对比图。误差很小，满足杂波抑制算法的精度要求。

表3和表4分别给出了GPU常规算法和本发明提出的改进算法各步骤精确耗时及算法加速比。其中CPU端和GPU端程序运行时间均在release模式下测得。CPU端编译器设置为“使速度最大化(/O2)”，并利用系统高精度测时函数QueryPerformanceFrequency和QueryPerformanceCounter多次运行取平均求得。GPU端则使用NVIDIA官方提供的NVIDIA VisualProfiler直接测出结果，再多次取平均，具有更高的准确性。从结果来看，本发明提出的算法同时对7个通道数据(每个通道点数为256K点，抑制距离元数为128)进行杂波抑制耗时仅17.215ms。若加上数据从CPU到GPU的传输耗时(5.4ms)，总耗时为22.6ms(GPU端抑制结果由于要继续进行匹配滤波处理，因此理论上无需拷回CPU端)，最终加速比为152.0X，显然本发明所提改进算法大大降低了算法计算量，减小了算法耗时。

表3GPU常规和改进算法各步骤(序号参见图5)耗时(us)

表4算法总耗时与加速比

表5给出了GPU常规算法和本发明提出的改进算法各部分及总显存占用量。改进算法显存占用近似为常规算法的(2M-1)/(K+2M-1)。以仿真实验为例，各点为单精度浮点复数(占用8字节)，N＝256*1024，K＝128，M＝8。则常规算法和改进算法总显存占用量近似为286MB和30MB，改进算法显存占用约为常规算法的1/9.5，显然本发明提出的改进算法大大降低了显存需求。

表5算法显存占用量(单位：点)

本发明实施例的效果可通过空中目标探测的外场实验结果进一步说明。

图10为武汉大学外辐射源雷达实验站某段FM实测信号杂波抑制前距离多普勒谱，其中信号中心频率为103.8MHz，带宽为500KHz。图11为该段数据杂波抑制后通道1距离多普勒谱,抑制后多个目标凸显，证实本发明所提方法的有效性。采用本发明所提方法对7个通道同时进行杂波抑制(每个通道数据点数为256K点，约0.524s)，抑制距离元数为256，采用双精度浮点计算，总耗时仅116.5ms。

本发明提出的杂波抑制处理方法算法耗时及显存占用显著小于现有算法，且处理时间稳定，验证了本发明方法具有良好的稳定性和加速性能，能够满足外辐射源雷达杂波抑制处理的实时性要求。

本文中所描述的具体实施例仅仅是对本发明精神作举例说明。本发明所属技术领域的技术人员可以对所描述的具体实施例做各种各样的修改或补充或采用类似的方式替代，但并不会偏离本发明的精神或者超越所附权利要求书所定义的范围。

Claims (5)

1.一种利用GPU实现外辐射源雷达多通道时域杂波抑制的方法，其特征在于，包括以下步骤：

步骤1、获取用于杂波抑制处理的外辐射源雷达信号数据：具体是首先将外辐射源雷达的接收阵列各个阵元接收到的回波信号经过模数转换、滤波抽取变换到数字基带，然后将变换到数字基带后的数据传送至主机端，在主机端对各个通道数据乘以相应通道校准值进行校准；

步骤2、在主机端将步骤1接收的数据进行初始化：主机端使用函数malloc为中央处理器CPU分配计算所需内存，使用函数cudaMalloc为图形处理器GPU分配计算所需显存空间；

步骤3、主机端将接收的数据传送至图形处理器GPU显存中，并定义：

监测通道数为M,杂波抑制距离元长度为K，抑制数据长度为N；

参考通道数据为s_ref＝[s_ref(0) s_ref(1) … s_ref(N+K-2)]^T；

监测通道数据为S_sur＝[s_sur,0 s_sur,1 … s_sur,M-1]；

其中，每个监测通道数据为s_sur,m＝[s_sur(K-1) s_sur(K) …s_sur(N+K-2)]^T，m为监测通道编号，T代表转置；

步骤4、在图形处理器GPU中对各通道数据进行并行杂波抑制处理，即杂波抑制结果 ${\tilde{S}}_{\sup }=S_{\mathrm{sur}}-{\tilde{S}}_{\mathrm{sur}}=S_{\mathrm{sur}}-X_{\mathrm{ref}}{\left(X_{\mathrm{ref}}^H{X}_{\mathrm{ref}}\right)}^{-1}{X}_{\mathrm{ref}}^H{S}_{\mathrm{sur}};$ 包括以下子步骤：

步骤4.1、利用参考通道数据采用迭代算法计算自相关矩阵 $R_x=X_{\mathrm{ref}}^H{X}_{\mathrm{ref}},$ 具体是：

步骤4.11、首先计算矩阵R_x第i＝0行值，即R_x(0,j)，其中0≤j≤K-1；

步骤4.12、计算出第i＝i+1行的值；若列序号j<i,则否则通过利用关系式：

$R_x(i,j)-R_x(i-1,j-1)=s_{\mathrm{ref}}^{*}(N+i-1)s_{\mathrm{ref}}(N+j-1)-s_{\mathrm{ref}}^{*}(i-1)s_{\mathrm{ref}}(j-1)$ 及对角线元素R_x(i-1,j-1)计算R_x(i,j)；

步骤4.13、重复步骤4.12，若i＝K-1则计算完毕后停止；

步骤4.2、对自相关矩阵R_x进行原地求逆得到

步骤4.3、利用参考通道数据和监测通道数据进行互相关运算得到 $C=X_{\mathrm{ref}}^H{S}_{\mathrm{sur}};$

步骤4.4、计算及互相关矩阵C的乘积得到矩阵

步骤4.5、计算参考通道数据和D矩阵乘积得到监测通道中直达波和多径杂波分量 ${\tilde{S}}_{\mathrm{sur}}=X_{\mathrm{ref}}D;$

步骤4.6、从原始监测通道数据减去杂波分量得到剩余回波分量，即杂波抑制结果 ${\tilde{S}}_{\sup }=S_{\mathrm{sur}}-{\tilde{S}}_{\mathrm{sur}};$

其中，s(i)＝|s_ref(i) s_ref(i+1) … s_ref(N+i-1)|^T，X_ref＝|s(0) s(1) … s(K-1)|；

步骤5、将杂波抑制结果在图形处理器GPU中继续进行其他处理或传送至中央处理器CPU的内存中并保存；

步骤6、在主机端使用函数free释放中央处理器CPU所分配的内存，使用cudaFree释放图形处理器GPU中所分配的显存。

2.根据权利要求1所述的利用GPU实现外辐射源雷达多通道时域杂波抑制的方法，其特征在于，步骤4.3与步骤4.5计算C矩阵和矩阵时直接在GPU核函数中利用s_ref进行计算。

3.根据权利要求2所述的利用GPU实现外辐射源雷达多通道时域杂波抑制的方法，其特征在于，所述步骤4.11计算矩阵R_x第i＝0行值，是利用s_ref来计算矩阵R_x第i＝0行值，核函数每个线程块计算num个点，num取值可为1,2或4，共分为K/num个线程块，每个线程块为m个线程，每个线程分配m个寄存器rA，1个寄存器rB,m个寄存器rC；每个线程块分配m+num-1个共享存储器sA；其中，计算第i个线程块包括以下子步骤：

步骤4.111、将结果寄存器数组rC清零，序号j赋值为0；

步骤4.112、从向量s_ref偏移i*num+j处加载(m+num-1)个数据到共享存储器sA中；加载sA前num个数据到寄存器rA中；从向量s_ref偏移j处加载m个数据到线程块各自寄存器rB中；

步骤4.113、计算rA数组共轭与rB乘积结果，并累加到数组rC中

步骤4.114、给j增加m；若j<N,重复执行步骤4.112,4.113,4.114，否则进入步骤4.115.

步骤4.115、依次将rC寄存器数组中的数据加载到共享存储器sA中，进行并行规约求和，并将求和结果保存在R_x矩阵第i*num到第(i+1)*num-1处。

4.根据权利要求1所述的利用GPU实现外辐射源雷达多通道时域杂波抑制的方法，其特征在于，所述步骤4.4中，计算互相关矩阵C是利用s_ref和S_sur进行计算，核函数分成K/2个线程块，每个线程块包含K个线程；每个线程分配长度为1的寄存器rA，rA0，长度为M的寄存器数组rC；每个线程块分配KM长度的共享存储器sC；每个线程块计算C矩阵2行元素，计算第i和第i+1行元素的线程块按如下步骤进行：

4.31、从s_ref偏移i和i+1处加载2个元素u和v至寄存器rA,rA0中，从S_sur中加载第一行共M个元素pv到GPU共享存储器sC中，M为S_sur矩阵列数；

4.32、分别计算u、v与pv对应元素乘积，保存结果1行存放在寄存器数组rC中，另1行保存在GPU共享存储器sC中；

4.33、继续加载s_ref下2个元素继续保存至寄存器rA,rA0中，从S_sur中加载下一行共M个元素pv到GPU共享存储器sC中，重复步骤4.32，并将计算结果第1行对寄存器数组中结果进行乘加运算，计算结果第2行对GPU共享存储器进行乘加运算；

4.34、重复步骤4.33，直至把S_sur全部N行加载计算完毕后，把寄存器数组rC中结果保存至C矩阵第i行，把共享存储器sC中结果保存至C矩阵第i+1行，计算完毕。

5.根据权利要求1所述的利用GPU实现外辐射源雷达多通道时域杂波抑制的方法，其特征在于，所述步骤4.5中，利用s_ref和D计算直达波和多径杂波分量核函数分成二维网格(N/K,32)规模的线程块，每个线程块划分成二维线程(32，K/32)，用于计算矩阵的K/32行，每行共M个元素，也即监测通道数目；每个线程分配长度为M的寄存器数组rC；每个线程块分配长度为(K+K/32-1)的共享存储器sA；计算第i个线程块包括以下子步骤：

4.51、将结果寄存器数组rC清零，从向量s_ref偏移处i*(32N/K)处加载K+K/32-1个元素到共享存储器sA中；

4.52、分别以sA[0]，sA[1]，…，sA[K/32-1]为起点计算与D矩阵第j＝0列乘积，并累加到rC寄存器数组中；

4.53、重复步骤1，取D矩阵第j＝1,2，…，M-1列乘积，并累加到rC寄存器数组中；

4.54、保存rC中结果到矩阵第i*K/32到第(i+1)*K/32-1行。