一种LDLT分解的外辐射源雷达杂波抑制并行处理技术

贾 东，温 博，施 健，罗扬静，王海涛

(1.桂林电子科技大学 信息与通信学院，广西 桂林 541002；2.中国电子科技集团公司第五十四研究所，河北 石家庄 050081)

0 引言

外辐射源雷达是一种利用第三方辐射源来探测目标的双多基地雷达,具有隐蔽性好、抗干扰能力强、反隐身能力和防低空突袭能力强等诸多优势,受到了国内外的广泛关注[1]。但外辐射源雷达接收信号中不仅包括目标回波信号,还包括能量很强的直达波和多径杂波,从而导致微弱目标回波信号被杂波副瓣掩盖。因此需要对其进行杂波抑制,消除强杂波对目标检测的影响[2-4]。

目前外辐射源雷达[5]主流的杂波相消算法是扩展相消批处理算法(Extensive Cancellation Algorithm Batches, ECA-B)[6-8],其分段利用参考信号的多普勒和时延矩阵来构建多径和杂波空间,然后将回波信号投影到此空间的正交子空间从而实现杂波消除,该方法不仅可实现静止杂波和慢动目标杂波的消除,同时具有很好的稳健性。但近年来随着可用的照射源信号带宽越来越宽(如传统调频(Frequency Modulation, FM)信号的有效带宽为80 kHz[9],目前第5代移动通信技术(5th Generation Mobile Communication Technology, 5G)信号带宽信号已达到50 MHz[10]),一方面大幅度提升了外辐射源的探测能力,另一方面导致在对消相同距离杂波的情况下,杂波对消算法所需的阶数越来越大和计算复杂度越来越高。因此近年来,在外辐射源雷达实时跟踪系统中,如何提升ECA-B的实时处理性能成为研究热点。

目前高性能并行计算中,中央处理器(Central Processing Unit, CPU)+图形处理器(Graphic Processing Unit, GPU)并行异构由于成本低、开发难度小、体积小、便于携带等特点被广泛应用于外辐射源雷达信号处理系统中[11-14]。

文献[11]利用GPU求解失配滤波因子,其将所有数据分段进行并行处理。文献[12]对拓展的LS算法进行GPU加速实现,其中矩阵相乘采用CUBLAS库函数完成。CUBLAS库中的乘法函数计算行列差距较大的矩阵时,会产生极大的空间冗余。文献[13]采用双GPU对ECA-B进行加速,每个GPU消除一半回拨信号中的杂波。文献[14]在实现过程中不需要构造参考矩阵,节省大量空间,且改进了自相关矩阵的计算,减少了大量计算冗余。上述传统文献都是以高斯消元法为基础原理对矩阵进行求逆,且文献[14]调用统一计算设备架构(Compute Unified Device Architecture, CUDA)流实现并行计算。CUDA流并不能使所有分段中同一计算环节的核函数并行,使得算法计算耗时长。虽然高斯消元法具有计算时间短的特点,但是需要在GPU与CPU之间多次传输数据,在数据量非常大且为双精度的情况下,数据传输时间将会大于计算时间,不利于算法的实时处理[15-16]。因此随着辐射源信号带宽变大而带来的巨大的数据量,且数据精度要求变高,传统算法已经难以满足高精度实时处理需求。为此,本文基于GPU多线程并行处理技术并结合ECA-B各段相同子模块特性,使得ECA-B的各子模块之间并行;然后,针对传统ECA-B求逆过程中数据传输耗时问题,利用自相关矩阵共轭对称特性提出一种基于LDLT的并行迭代求逆方法[17],通过2个CUDA核函数实现求逆处理,解决了双精度数据的矩阵求逆过程中数据传输耗时问题,进一步提升ECA-B的实现效率。实验结果表明,与传统算法相比,本文提出的算法具有更高的时效性和有效性。

1 ECA-B原理

假定参考信道接收信号为Sref(t),目标回波信号为Ssur(t),辐射源信号为S(t)。

Sref(t)=S(t-τ)+nref(t),

(1)

(2)

式中:τ为参考信道回波相对于发射信号的时延,βi、τi分别为M条多径中第i条回波衰减和时延,其中i=0时这条多径信号为直达波;αk、τk、fk分别为K个目标中第k个目标的衰减、时延和多普勒频移,nref(t)、nsurv(t)分别为参考信道和目标回波信道中的噪声,一般情况可以看作高斯白噪声。

第h段参考信号构成的滑动矩阵为:

(3)

式中:L为数据长度,H为数据分段数,h代表分段的第h段,a为杂波抑制的对消阶数。

基于此,每段目标回波信号自适应滤波权值系数可以转化为一个求解凸优化问题:

(4)

求式(4)是个最小二乘优化问题,其解析解为:

(5)

那么,第h段的剩余信号可以通过第h段的回波信号减去第h段的滑动矩阵与权值向量的积得到。

(6)

依照上述算法原理,可将ECA-B的数据处理流程总结如下:

① 对参考信号Sref(t)和回波信号Ssur(t)分段;

② 利用每段参考信号构建滑动矩阵Vh;

③ 利用第h段回波信号Ssur,h和滑动矩阵Vh求解权值向量;

④ 进行杂波抑制得到第h段的剩余信号。

算法具体流程如图1所示。

图1 ECA-B计算流程Fig.1 The calculation flowchart for ECA-B

各分段之间杂波相消是相互独立且处理流程相同,因此段间并行性很高。但是在传统计算自相关矩阵过程中需要构造参考矩阵Vh,采用CUDA流实现各分段并行。而CUDA流只能使数据传输和核函数并行,并不能使各分段同一计算环节之间并行计算。基于ECA-B各分段计算流程相同,本文采用段间并行算法使各分段同一计算环节在同一个核函数中并行实现。此外所有分段一起计算导致算法中的自相关矩阵求逆环节数据量非常大,而传统的自相关矩阵求逆方法是高斯消元法,该方法进行N阶矩阵求逆时,需要传输3×N次数据。在采用双精度浮点数数据形式时,这种设备间传输大批量数据大大增加了杂波抑制的处理时间[15],基于此问题,本文采用LDLT分解法只需传输一次数据,降低了数据传输时间。

2 算法的改进和并行实现

2.1 EDA-B的分段并行实现

目前CUDA并行主要用CUDA流实现并行计算,CUDA流的优点在于数据传输的同时调用核函数,不能使核函数之间直接并行计算。此外,随着辐射源数据量不断增大,且高精度要求,使得ECA-B各子模块处理时长变长,从而整体时长也会大大增大,因此对于ECA-B来说,并行度不高[18]。针对此问题,本文将算法中各分段的同一环节一起计算,提高了各分段间的并行度,减小数据量变大对算法实时处理产生影响。具体来说,此方法利用CUDA编程特点,通过线程索引控制每个线程访问全局内存中具体位置的数据,利用这一特点实现线程块分开计算各分段数据[19]。根据上述特点,把ECA-B所有分段的同一变量存储在同一段连续存储器内。例如将数据分为10段计算,在存储自相关矩阵Rh时,R1、R2、…、Rh、…、R10可按行优先顺序存储,而矩阵Ch、Wh、Sh、Vh在计算机中按列优先顺序存储,Rh的公式为:

(7)

在实现ECA-B的计算过程中,用线程块独立计算不同段的计算环节。例如以Wh计算为例,线程网分配10个线程块,每个线程块开辟10个线程。10个线程块对应计算ECA-B 10段中Wh,利用x维度索引矩阵的行数和列数,y维度索引具体某一段中的矩阵,即线程块Block0(线程号0～255)计算矩阵W0,线程块Block0内线程中迭代计算W0(i,i+k),线程块Blockh(线程号256～511)计算矩阵Wh,线程块Blockh内线程中迭代计算Wh(i,i+k),以此类推。以为Wh例,一个核函数同时计算10个块的示意如图2所示。

图2 GPU计算矩阵C示意Fig.2 Schematic diagram of matrix Cfor GPU calculation

同理,图1中的矩阵Rh、Wh、Vh、Sh与图2类似,也是将数据按一定规则存储,用不同的线程块去读取和计算对应的地址数据,这样即可实现块间并行。此外,此并行框架由于访问内存并不完全为顺序访问,所以计算时间并不是每段单独计算的k倍,而是各分段计算时长中最长的部分。算法实现流程如图3所示。图3中Vh为参考矩阵,Rh为自相关矩阵,Ssur,h为第h段回波信号,Wh为式(4)中的权值向量。ECA-B的数据处理流程总结如下:

图3 算法实现流程Fig.3 Algorithm implementation flowchart

① 在各分段中计算自相关矩阵Rh;

② 对自相关矩阵Rh进行求逆;

③ 利用第h段回波信号Ssur,h,Vh和求逆后的自相关矩阵Rh求解权值向量Wh;

④ 进行杂波抑制得到第h段的剩余信号。

2.2 LDLT分解求逆

传统的ECA-B实现过程中对自相关矩阵求逆通常采用高斯消元法。此方法需要多次进行GPU与CPU之间的数据传输,导致数据传输时间大于数据计算时间。而外辐射源雷达的数据量非常大,使得高斯消元法的数据传输时间更大,因此在CPU与GPU传输带宽小的设备中,高斯消元法不适合处理ECA-B。针对这一问题,本文根据厄米特矩阵Rh的对称特性,采用LDLT分解法,通过2个GPU核函数实现自相关矩阵求逆。其中,一个核函数实现LDLT分解,另一个核函数实现矩阵求逆,并将其应用到段间并行实现中,从而减少数据传输耗时。

LDLT分解法由Cholesky分解(Cholesky Factorization)法改进而来,Cholesky分解法通过对三角矩阵L的列向量归一化可以得到改进后的Cholesky分解法为A=LDLH,其中D为对角矩阵,L为对角线上都是1的下三角矩阵,LH为L的复共轭转置矩阵。当对矩阵A求逆,即对角矩阵D-1的对角元素矩阵D对角元素倒数,计算量很小,所以对矩阵A求逆主要就是对下三角矩阵L求逆。

(8)

根据矩阵乘法计算,再化简可以得到:

(9)

式中:gij为过渡矩阵G的元素,lij为下三角矩阵L的元素,dj为对角矩阵D的元素。过渡矩阵G也是下三角矩阵,辅助计算下三角矩阵L和对角矩阵D。

根据逆矩阵定义LB=E,推导可得:

(10)

根据式(9),当j=1时,矩阵G的第一列等于A的第一列,即gi1=ai1。当计算矩阵G的第j列元素时,需要用到L的第j行的所有非零元素,而每行中的单个元素计算与其他行的元素互不影响,因此改进型Cholesky分解法可以并行迭代实现。一个线程迭代分解得到矩阵G的一行数据,线程中每次迭代计算需要用到上一次的结果,所以需要所有线程同步。而线程同步函数__syncthreads()仅保持同一个线程块内所有线程同步,当矩阵维度大于线程块内最大线程数时,迭代放在CPU端控制,即CPU控制核函数调用,核函数迭代最大线程数次。具体GPU实现流程如下:

① 计算矩阵L的第一列,即gi1=ai1,1≤i≤C;

② 利用G和L前j-1行的值计算矩阵G第j列的值,线程同步;

③ 计算矩阵D和L的第j列的值;

④ 重复计算步骤②和③。

LDLT分解求逆实现示意如图4所示。

图4 LDLT分解求逆示意Fig.4 Inversion process based on LDLT decomposition

3 算法性能分析

为了评估本文算法的性能,在以下环境中进行性能测试:CPU为Intel(R)Core(TM)i5-6200F CPU@2.30 GHz,内存4 GB,图形处理器为NVIDIA GeForce GT 940M显卡,2 048 MB显存;操作系统为Windows 10,配置了Matlab R2018a、Microsoft Visual Studio 2019和CUDA 11.0。

实验中数据为实测接收的8路GSM信号,每路信号数据点数为80 000,取第一组作为参考信号,给每路信号加上6组不同的多径干扰,再加上同一目标回波信号,对8路信号进行波束形成处理之后再消除杂波。采样率fs为250 kHz,数据分段数为10,采用双精度浮点数运算,重复实验50次。为了验证算法的性能,矩阵求逆模块的验证,将本文的LDLT算法与文献[14]中的求逆算法进行对比,ECA-B的整体性能与文献[14]进行对比。性能评价所采用的指标为算法加速比。算法加速比为50次重复实验下对比文献[14]算法与本文算法的平均处理时间之比。

为了验证本文算法的杂波抑制效果,分别用Matlab和GPU进行杂波消除。图5是Matlab对信号进行杂波抑制的距离-多普勒频谱。图6是GPU对信号做杂波抑制的距离-多普勒频谱图,杂波被消除,目标清晰可见。可以看出二者几乎一样,都得到了目标的距离-多普勒频谱。图7给出了Matlab和GPU计算结果的绝对误差,误差低于7×10-15。实验结果说明了本文算法的有效性。

图5 Matlab杂波抑制距离-多普勒频谱Fig.5 Clutter suppression range-Doppler spectrum in Matlab

图6 GPU杂波抑制距离-多普勒频谱Fig.6 Clutter suppression range-Doppler spectrum by GPU

图7 GPU和Matlab计算的绝对误差Fig.7 Absolute error calculated by GPU and Matlab

随着辐射源信号的不断增高,外辐射源处理的数据量也越来越大,且数据精度要求越来越高。针对此,文献[14]中的矩阵求逆算法传输耗时问题使其不再适用于ECA-B中的自相关矩阵求逆。本文将LDLT算法与文献[14]中的矩阵求逆算法进行了对比,实验数据由Matlab生成的双精度Hermitian矩阵。实验结果如图8所示,由图8可以看出,矩阵阶数越大,LDLT算法耗时比文献[14]中的矩阵求逆算法越少。矩阵阶数从64阶增长至512阶,LDLT算法时间增大了162倍,而文献[14]中的矩阵求逆算法增大了424倍,相比之下,LDLT算法更适用于处理大批高精度数据。

图8 求逆实验结果对比Fig.8 Comparison of inversion experiment results

本文段间并行实现将所有分段一起计算,而文献[14]采用CUDA流的方法实现了ECA-B的段间串行算法。所有分段一起计算也导致了矩阵求逆环节的数据量增大了h倍,同时文献[14]中的矩阵求逆算法的数据传输时间增大。基于此,本文在段间并行实现中将2种求逆算法再次进行对比,当对消阶数为128时,段间并行算法中的矩阵求逆采用文献[14]中的矩阵求逆算法,ECA-B处理时间为133 ms,而当求逆方法采用LDLT分解法时,ECA-B处理时间则为45 ms。很明显,在段间并行算法中,LDLT分解法优于献[14]中的矩阵求逆算法。本文算法与文献[14]仿真结果对比,结果如图9所示,最大加速比为3.5倍,大大缩短了算法运算时间。从图9中可以看出,当对消阶数在128阶以内时,段间并行算法的处理时间维持在45 ms以内,而当对消阶数增大到256时,段间并行算法也仅仅只有296 ms,可以看出对消阶数的增大对段间并行算法的影响较小,而随着对消阶数的不断变大,对文献[14]的段间串行算法影响较大。由此可以看出,本文的段间并行算法和LDLT求逆算法有着非常好的时效性能。

图9 杂波抑制实验结果对比Fig.9 Comparison of clutter suppression experiment results

为了验证ECA-B消除杂波性能,本文将ECA-B与传统算法中的SMI算法性能进行对比,将目标回波信号经过杂波抑制处理之后进行距离多普勒处理。图10为ECA-B的距离-多普勒处理中多普勒维,图11为SMI算法的距离-多普勒处理中多普勒维。从图中可以看出,二者都可以清晰看出目标信号,并且ECA-B在0多普勒通道产生明显零陷,直达波和多径杂波被完全抑制,杂波抑制效果较好。ECA-B和SMI算法的区别在于ECA-B将信号进行分段处理,因此适合GPU将所有分段并行处理。

图10 ECA-B的多普勒维Fig.10 Doppler dimension of ECA-B

图11 SMI算法的多普勒维Fig.11 Doppler dimension of SMI algorithm

4 结束语

随着辐射源带宽的不断增加,ECA-B对消阶数也在不断增加,使得算法所需处理的数据量不断变大,对数据精度的要求也在不断变高。ECA-B的分块处理过程需要极大的时间和空间复杂度,本文将ECA-B的所有分段在GPU上统一进行段间并行计算,并在GPU上实现矩阵的LDLT分解求逆,矩阵的LDLT分解求逆相对于Gauss-Jordan顺序消去法有着很好的加速效果。ECA-B很好地满足了外辐射源雷达对实时性的要求。