基于图特征学习的海杂波抑制算法

陈鹏， 许震， 曹振新， 王宗新

(东南大学 信息科学与工程学院 毫米波国家重点实验室， 江苏 南京 210096)

收稿日期：2021-10-09

基金项目：江苏省优秀青年科学基金项目(BK20220128)；国家重点研发计划项目(2019YFE0120700)；中国航天科技集团公司第八研究院产学研合作基金资助项目(SAST2021-039)；国家自然科学基金项目(61801112)

0 引言

在海洋探测雷达系统中，雷达照射海面所产生的后向散射回波称为海杂波[1]。由于海洋雷达受到复杂海况所产生的海杂波影响，极大地降低了雷达监视和检测目标的能力[2]。因此，从海洋雷达接收的信号中抑制海杂波分量是亟需解决的问题。

由于海面的复杂运动，导致海杂波具有非线性、非高斯、非平稳特征[3]。传统的随机信号处理方法一般建立在信号是平稳性的假设前提下，单纯采用时域或频域的单一特征无法有效地描述海杂波的复杂特征。通过引入时频变换的理论构造基于时域和频域的联合二维处理方法，将一维的时间信号投射成时间- 频率的二维图像，可以得到描述时间- 频率- 能量三者关系的时频谱，更加直观地观测到信号在时频域上的局部变化规律[4]。目前，针对海杂波时频变换的主要方法有短时傅里叶变换(STFT)、魏格纳- 维尔分布(WVD)和小波变换(WT)等[5]。其中，STFT计算简便，但时频分辨率无法同时达到最优；WVD时频分辨率可以同时达到理想情况，但在分析多分量信号时会产生严重的交叉项干扰；WT的多尺度分辨特性适合于分析非线性、非平稳信号，但受到小波基函数的约束，导致适应性较弱。

随着近年来机器学习领域的迅速发展，人工智能算法也广泛运用于海杂波的处理领域中。针对海杂波的仿真建模，可以通过反向传播神经网络和门控循环神经网络进行幅度预测和参数估计[6-7]。在海面的弱小目标检测中，深度学习和迁移学习被广泛运用[8-10]。在海杂波抑制领域中，与传统的海杂波抑制算法，如基于子空间分解的奇异值分解(SVD)算法[11-12]和基于时域对消的自适应杂波抑制算法[13]比较，运用经典的神经网络模型同样可以产生较好的抑制效果，如基于AlexNet的自适应杂波智能抑制方法[14]。

在人工智能领域，Joan Bruna等在卷积神经网络(CNN)的基础上提出了应用于图结构的图卷积网络(GCN)[15]，使得图结构的特征分析越来越多地用于传统图像的处理中。Aditya Grover等在词向量的数据处理算法word2vec[16]基础上提出了应用于图结构中的数据处理算法node2vec，该算法通过学习图结构的特征构造出高维的节点特征向量，完成了节点分类和链路预测等任务[17]。

本文分析了海杂波的非平稳特性，发现了其时频谱与目标回波信号的时频谱之间具有明显的结构特征区别。因此，提出基于图特征学习的海杂波抑制算法，将雷达接收信号的时频谱构造成图结构的形式，并采用node2vec算法对图进行图嵌入处理，使用K-means聚类算法将图嵌入的节点特征向量分类，实现在时频域中海杂波和目标回波信号的分离，从而抑制了雷达接收信号中的海杂波分量。最后进行算法的仿真以及实测数据对比实验，分析了该算法在不同超参数情况下的节点分类结果以及海杂波抑制效果，并给出了适合工程化应用的研究方向。

1 雷达接收信号模型

如图1所示，机载或舰载雷达对海上目标进行监视或检测时，海杂波会严重恶化雷达探测性能，此时雷达的接收信号r(t)可以表示为三部分：

r(t)=s(t)+x(t)+n(t)

(1)

式中：s(t)为目标回波信号；x(t)为海杂波回波信号；n(t)为复高斯白噪声。海杂波的抑制问题就是在不损失目标回波信号s(t)的情况下，最大程度地抑制杂波信号x(t)，提高信杂比(SCR)。

图1 海洋雷达探测示意图Fig.1 Schematic diagram of ocean radar detection

经过基本的雷达信号处理后，海上目标的回波信号s(t)可描述为等效复基带形式

(2)

式中：K为目标个数；a(k)为第k个目标回波信号的雷达散射截面(RCS)；fd(k)为第k个目标的多普勒频率。

为了在海杂波干扰下有效地提取目标信息，采用STFT对雷达接收信号做预处理。STFT是典型的线性时频变换，其计算简便且不存在交叉项干扰的问题[18]，定义为

(3)

式中：y(·)为初始信号；STFTy(t,f)表示针对y(·)信号的短时傅里叶变换；g(·)为窗函数，上标*为复共轭。对式(1)两边同时进行STFT，可以得到

STFTr(t,f)=STFTs(t,f)+STFTx(t,f)+STFTn(t,f)

(4)

式中：STFTr(t,f)、STFTs(t,f)、STFTx(t,f)、STFTn(t,f)分别表示针对r(t)信号、s(t)信号、x(t)信号和n(t)信号的短时傅里叶变换。

2 基于图特征学习的海杂波抑制算法

2.1 时频谱的图构造

对于雷达接收信号中的目标回波信号和非平稳的海杂波信号，在经过STFT后均可以在时频谱中表现出频率随时间变化的规律，且两者的时频谱在频率上均是连续的。

通过对雷达接收信号的时频谱分析，可以发现海杂波与目标回波信号时频谱中的能量均主要分布在低频部分，但由于海杂波的非平稳特性与目标回波信号的确定性，导致两者在时频谱中的能量分布情况存在着较大的结构特征差异，具体表现为：目标回波信号的能量分布比海杂波的能量分布更为集中与连续。这里将离散分布的海杂波时频谱描述为结构相似性，将集中分布的目标回波信号时频谱描述为社群相似性。

因此，雷达接收信号的时频谱具有图的结构形式，可以进行图构造处理。将其时频谱STFTr(t,f)进行图构造的步骤如图2所示。

图2 时频谱的图构造流程
Fig.2 Graph construction flow of time spectrum

时频谱的图构造具体过程为：

步骤1提取式(4)中雷达接收信号时频谱STFTr(t,f)图像中的像素点，构造成二维初始矩阵S，大小为r×c。将时频谱转换为矩阵的形式，便于进行后续的算法处理。矩阵S表示为

(5)

步骤2设定阈值λ，记矩阵S中元素的最大值为sm(即第m个元素为最大值)，对于矩阵S中的任意元素si,j，若si,j/sm<λ，则将si,j的值更新为0。记更新后的矩阵S中非0元素的个数为N。将矩阵中较小的元素赋值为0，以有利于降低矩阵的冗余度和算法的时间复杂度。此时矩阵S表示为

(6)

式中：为简化描述，设置初始矩阵S中任意r×c-N个元素为0。

步骤3将矩阵S中为0的元素不变，其余N个非0元素重新赋值。具体的赋值方法为按照元素的下标顺序，进行从1-N的连续递增的整数赋值。图结构中最重要的是节点之间的位置关系，使用连续递增的整数赋值代替原本矩阵中的非0元素，便于进行图结构中的节点关系表征。赋值更新后的矩阵S表示为

(7)

式中：4

步骤4此时矩阵S中有N个非0的连续元素值，将这N个元素看作图结构中的N个不同节点，根据N个不同节点在矩阵S中的位置关系，构造出相对应的图结构的邻接矩阵Ar。

具体构造过程为：对矩阵S进行遍历，每次所遍历的区域大小为m×m，将遍历区域记作Q，其余区域记作Q′，则矩阵S在遍历过程中由Q和Q′组成，将遍历的移动步长设置为1。在每一步遍历中，将区域Q中节点之间的关系视为相关节点，区域Q中节点与区域Q′中节点之间的关系视为无关节点。在邻接矩阵Ar中，将相关节点之间的权重设置为1，无关节点之间的权重设置为0，进而构造出大小为N×N的邻接矩阵Ar。采取遍历矩阵的方式，可以将m×m范围内的节点看作是结构相似的，从而达到通过构造邻接矩阵Ar来表征雷达接收信号时频谱中不同节点之间结构位置关系的目的。

图3 遍历矩阵S的范围和所对应的邻接矩阵ArFig.3 The range of the traversal matrix S and the corresponding adjacency matrix Ar

为简化描述，设置m=2，图3为依次遍历矩阵S中前4个2×2的范围后得到相对应的邻接矩阵Ar。对矩阵S以步长大小为1向右继续移动遍历，到达每行的末端时会向下移动1行，再从行首继续以相同的方式遍历，直到遍历完成整个矩阵S。

步骤5通过步骤1～步骤4，可以得到式(4)中雷达接收信号的时频谱STFTr(t,f)对应的邻接矩阵Ar。通过Ar得到该信号的时频谱对应的无向图Gr，表示为

Gr=(Vr,Er,Wr)

(8)

式中：Vr={vr1vr2vr3…vrN}为图Gr中的所有节点集合；每两个节点之间可能会存在边；Er为这些边的集合；Wr为每条边上的权重共同构成的权重矩阵。这里只画出图3中邻接矩阵Ar对应的图Gr，如图4所示。

图4 邻接矩阵Ar对应的图GrFig.4 Graph Gr of adjacency matrix Ar

2.2 基于node2vec的图嵌入算法

图嵌入是用来挖掘图的结构特征的方法。node2vec是一种典型的图嵌入算法，该算法建立在传统随机游走策略[19-20]基础上，结合了广度优先搜索(BFS)和深度优先搜索(DFS)两种图结构的搜索策略，是一种带有偏置的随机游走策略，可以更加有效地检索分散的相邻节点[21]。

图4所示的图Gr中，演示node2vec算法中BFS和DFS的搜索策略，结果如图5所示。

图5 针对节点n2的两种搜索策略(步长k=3)Fig.5 Two search strategies for nodes n2(step k=3)

在图5中，如果从节点n2出发，当步长k=3时，采用BFS进行搜索所得到的节点集合是{1，n1，n3}；采用DFS进行搜索所得到的节点集合是{2，3，4}。从图5中可以发现，海杂波信号的时频谱STFTx(t,f)的像素点呈现出的结构相似性与BFS的搜索策略是一致的，而目标回波信号的时频谱STFTs(t,f)的像素点呈现出的社群相似性与DFS的搜索策略是一致的。因此，node2vec算法可以对由雷达接收信号的时频谱所生成的无向图进行图嵌入处理，使用高维向量表征每个节点，对同一种信号节点之间的相似程度以及不同信号节点之间的差异程度进行描述。

在node2vec算法中，有两个重要的超参数p和q，下面解释这两个超参数的作用。

图4所示的图Gr中，演示node2vec算法的随机游走搜索策略：当从节点1到达节点n2后，评估离开节点n2时对下一个即将到达的节点的搜索策略，具体如图6所示。

图6 随机游走的搜索策略Fig.6 Random walk search strategy

图6中，α为随机游走的搜索偏置，也可以被看作是从节点n2到下一个节点的转移概率，其取值范围为

(9)

式中：l1x表示节点1和节点x之间的最短距离，这里的节点x是指图6中的节点1、节点2、节点n1和节点n3。

可以发现，转移概率α受到超参数p和q的影响，因此超参数p和q的作用是调节随机游走的搜索策略走向，对其具体的理解为：

1)超参数p控制着再次访问已访问节点的概率(返回节点1)，超参数q控制着远离已访问节点的概率(走向节点2或者节点n3)；

2)p的值越小，随机游走的策略倾向于回溯节点1，当p的值越大，随机游走的策略将更接近DFS，倾向于不回溯节点1，降低了2-hop的冗余度；

3)q的值越小，随机游走的策略将更接近DFS，q的值越大，随机游走的策略将更接近BFS[22]。

因此，通过调节node2vec算法中的超参数p和q，可以控制随机游走序列中下一个节点的选择。超参数p和q的引入使得BFS和DFS可以交互存在，进而实现不同特性的搜索策略，得到适合本文中图Gr所对应的图结构搜索方式。

通过使用node2vec算法，将图Gr中的每个节点实现由低维度向高维度的映射，进而得到每个节点的高维特征向量。将node2vec算法的映射关系记作f，映射过程可以表示为

f:Vr→V′rvri=[i]→v′ri=[v′ri1v′ri2…v′rid]

(10)

(11)

式中：Vr和V′r分别表示图Gr映射前后的节点集合；vri和v′ri分别为在图Gr中映射前后的任意一个节点；d为维度为d的高维空间中的全体实向量集合。

可以发现，在进行图嵌入前，节点集合Vr中节点之间的相关性只能通过邻接矩阵Ar来表示，每个节点本身没有携带任何信息；在进行图嵌入后节点集合V′r中的每个节点都使用了d维向量进行表征，此时可以采用词向量数据处理算法word2vec[16]中衡量向量相似度的方法来分析节点之间的相关性，计算方法如下：

(12)

式中：D(v′ri,v′rj)为节点集合V′r中任意两个节点v′ri和v′rj之间的欧式距离；u用来衡量节点之间的相似度，范围为[0,1]，取值越大，表明节点之间的相似度越大。

2.3 基于K-means的节点分类算法

通过式(12)中节点之间的相似度u，可以得到雷达接收信号时频谱中的像素点在经过图嵌入处理后的相关性，依据相关性的大小，可以采用K-means聚类算法，将经过图嵌入处理后图Gr中的节点进行分类。这里将式(10)中的节点集合V′r划分成n类，表示为

(13)

式中：V′s表示目标回波信号的节点集合；V′x表示海杂波的节点集合。在图Gr中，目标回波信号中同一多普勒频率对应的节点是集中连续的，其节点之间的相似度u较大，通过K-means聚类算法会将这些节点划分为一类。根据式(2)，目标回波信号中有K个不同的多普勒频率分量，因此V′s是由K个节点集合构成的，V′x是由剩下的n-K个节点集合构成的。

节点集合V′r中的N个节点是与式(6)中矩阵S中的N个非0元素一一对应的，因此V′s和V′x也各自对应着矩阵S中的部分非0元素，矩阵S可以表示为

S=Ss+SxV′s→SsV′x→Sx

(14)

式中：矩阵Ss表示目标回波信号的时频谱对应的像素点矩阵；矩阵Sx表示海杂波信号的时频谱对应的像素点矩阵，由此可以对式(4)中雷达接收信号的时频谱STFTr(t,f)进行分离，得到其中的STFTs(t,f)部分和STFTx(t,f)部分，再去除STFTx(t,f)部分，就完成了在时频域上对雷达接收信号中海杂波回波信号分量的抑制。

图7 海杂波抑制算法Fig.7 Sea clutter suppression algorithm

因此，本文所提基于图特征学习的海杂波抑制算法的流程图(见图7)主要包括4个步骤：

1)对目标回波信号和海杂波信号数据进行STFT，得到时频谱；

2)对时频谱进行图构造处理，得到对应的邻接矩阵Ar和图Gr；

3)使用node2vec算法，对邻接矩阵Ar进行图嵌入处理；

4)使用K-means聚类算法，依据图Gr中节点的特征向量进行节点分类，实现海杂波的抑制。

3 仿真与实测结果及分析

3.1 参数设置

本文算法的参数主要由3个部分组成：

1)目标回波信号中的目标个数K，多普勒频率分量的取值集合fd，采用基于实测数据的海杂波建模方法[23-24]生成海杂波的仿真数据时需要设置尺度参数和形状参数；

2) node2vec算法中的超参数p和q，节点特征向量的维度d，随机游走的次数r；

3) K-means聚类算法中的分类参数n。

具体的仿真参数设置如表1所示。

表1 仿真参数设置Table 1 Simulation parameter setting

本文将超参数p的取值集合设置为{1，2，4}，将超参数q的取值集合设置为{0.25，0.50}，主要有两方面原因：

1)为减少在随机游走的过程中对同一个节点重复回溯的搜索行为，降低node2vec算法的冗余度，因此超参数p需要较大的取值；

2)通过分析本文需要进行图嵌入操作的图结构特征，可以发现目标回波信号节点之间的有序性要远大于海杂波节点之间的有序性，因此目标回波信号节点之间的结构关系更容易被探索，且目标回波信号时频谱的结构表现为社群相似性，故本文的node2vec算法需要更倾向于DFS的搜索策略，超参数p需要较大的取值，超参数q需要较小的取值。

3.2 评价指标

本文主要对node2vec算法、K-means聚类算法以及雷达接收信号中海杂波分量的抑制效果这3个方面进行评价。

1)对node2vec算法的有效性进行评价

计算雷达接收信号中目标回波信号节点之间的相似度平均值u1和目标回波信号节点与海杂波节点之间的相似度平均值u2，通过比较u1和u2的大小，可以表明node2vec算法所构造的特征向量的有效性。

2)对K-means聚类算法的有效性进行评价

将目标回波信号对应的节点作为正类节点，将海杂波回波信号对应的节点作为负类节点，则本文的节点分类问题本质上是一个二分类问题。针对本文的二分类模型，定义F1分数对其进行性能衡量，计算公式为

(15)

式中：P为分类模型的精确率；R为分类模型的召回率；TP表示将实际类别是正类的节点预测为正类的节点数目；FP表示将实际类别是负类的节点预测为正类的节点数目；FN表示将实际类别是正类的节点预测为负类的节点数目。

F1分数可以具体分为宏观F1分数F1Macro,s和微观F1分数F1Micro,s，计算公式分别为

(16)

式中：F1T和F1N分别表示正类和负类的F1分数。从二者的计算方式上可以看出，宏观F1分数是平等地看待各个类别，它的值会受到稀有类别的影响；微观F1分数则更容易受到常见类别的影响。宏观F1分数和微观F1分数的取值范围均为[0,1]，取值越大，则表明分类的效果越好。

因此，本文采用F1Macro,s和F1Micro,s对K-means聚类算法的有效性进行评价。

3)对雷达接收信号中海杂波分量的抑制效果进行评价

根据2.3节，在时频域中可以分别得到经过分类模型处理后的目标回波信号功率和海杂波回波信号功率。为了衡量本文算法对雷达接收信号中海杂波分量的抑制效果，提出了信号杂波比(SCR)的改善因子ISCR(ISCR)[25]，定义为

ISCR=SCRo-SCRi

(17)

式中：SCRi和SCRo分别表示使用本文的海杂波抑制算法处理雷达接收信号前后的信杂比。

3.3 仿真与实测结果分析

用于验证本文海杂波抑制算法的海杂波数据有以下两个来源：

1)仿真数据：采用基于实测数据的海杂波建模方法，所生成的100组海杂波仿真数据进行仿真验证；

2)实测数据：采用海杂波实测数据[24]，取第1 000个距离单元至1 100个距离单元的100组海杂波实测数据进行仿真验证。

将集合fd中的多普勒频率分量进行组合，得到表1中K在不同取值时的目标回波信号，分别添加至仿真数据和实测数据中。在K取值为1、2和3时，雷达接收信号的平均SCR分别为-15 dB、-12 dB和-9 dB。

本节将通过对超参数p和q的优化，时频谱的直观结果对比以及本文算法与现有算法的对比这3个方面进行仿真与实测结果的分析。

从3.1节中算法参数的设置可以看出，超参数p和q的取值是一个集合，因此首先需要通过仿真结果的评价指标选取本文算法中的最优超参数p和q。

首先，对K=1时的单目标环境进行仿真分析。使用node2vec算法进行图嵌入处理后所得节点向量之间的相似度u1和u2对比情况如表2和表3所示。从表2和表3中可以发现：当多普勒频率相同时，在不同的超参数情况下，相似度u1的值均远大于u2，因此使用node2vec算法所构造的特征向量可以将两种信号的节点进行有效区分；当超参数{p,q}={4，0.25}时，可以同时取到相似度u1的最大值和相似度u2的最小值，表明此时两种信号节点之间的区分度最大。

在3种多普勒频率下，图8和图9为使用K-means聚类算法后得到的仿真结果F1Macro,s和F1Micro,s。可以发现，二者均是超参数{p，q}={4，0.25}时取得最大值，表明此时得到了最好的分类效果。图10为使用本文海杂波抑制算法前后的ISCR对比图。由图10可以发现，当超参数{p,q}={4,0.25}时，雷达接收信号的信杂比得到了最大提升，其中仿真数据的ISCR和实测数据的ISCR最高可分别达到11.59 dB和10.13 dB。

表2 海杂波数据来源于仿真数据Table 2 Sea clutter data comes from simulation data

表3 海杂波数据来源于实测数据Table 3 Sea clutter data comes from measured data

图8 Macro-F1 scoreFig.8 Macro-F1 score

图9 Micro-F1 scoreFig.9 Micro-F1 score

通过对上述3种评价指标的分析，均是在超参数为{4，0.25}时取得了最好的分类效果，因此本文算法的最优超参数选取为{p，q}={4，0.25}。

图10 ISCRFig.10 ISCR

在3种不同的多普勒频率下，取最优超参数时，使用本文的海杂波抑制算法处理实测海杂波数据的时频谱如图11所示，其中第1行图片是使用抑制算法前雷达接收信号时频谱，第2行图片是使用抑制算法后从雷达接收信号中所分离出的目标回波信号时频谱。由图11可以发现，在单目标环境下，本文的海杂波抑制算法可以实现在时频域上进行雷达接收信号中海杂波分量的抑制。

设置本文算法的超参数为最优超参数的情况下，与SVD算法、自适应杂波抑制算法和基于AlexNet的自适应杂波智能抑制方法进行仿真比较，得到的ISCR对比结果如图12所示。

由图12可以看出，针对上述两种不同来源的海杂波数据，本文算法的ISCR均优于其他3种算法，进一步验证了该算法的有效性。

对K=2和K=3时的多目标环境进行上述相同的仿真分析，得到的最优超参数均为{p，q}={4，0.25}。再使用本文的海杂波抑制算法处理实测海杂波数据的时频谱如图13所示，其中第1行图片是雷达接收信号的时频谱，第2行图片是从雷达接收信号中分离出的目标回波信号的时频谱。可以发现，在多目标环境下，本文的海杂波抑制算法同样可以实现在时频域上进行雷达接收信号中海杂波分量的抑制。

图11 单目标环境下使用海杂波抑制算法前后的时频谱对比Fig.11 Time spectrum comparison before and after sea clutter suppression algorithm in single target environment

图12 4种算法进行海杂波抑制后的ISCR对比Fig.12 Four algorithms are used to compare ISCR after sea clutter suppression

图13 多目标环境下使用海杂波抑制算法前后的时频谱对比Fig.13 Time spectrum comparison before and after sea clutter suppression algorithm in multiple target environment

本节在单目标和多目标的环境下，分别对两种不同来源的海杂波数据下的雷达接收信号模型进行算法的仿真分析。首先对比在不同超参数情况下的分类效果，得出了最优超参数的取值集合；然后在最优超参数下，使用本文的海杂波抑制算法处理实测海杂波数据，从时频谱的角度验证了本文算法的可行性；最后将本文算法与其他3种算法进行分析对比，进一步验证了本文算法对抑制雷达接收信号中海杂波回波信号分量的有效性。

4 结论

通过对雷达接收信号时频谱的结构分析，将图特征学习的理论引入到海杂波的抑制领域中。本文通过对时频谱的图嵌入和分类算法，实现了雷达接收信号中海杂波分量的有效抑制，同时验证了图论方法应用于海杂波处理中的可行性。新算法流程直观简单，但需要进行时间复杂度较高的海杂波数据的图构造处理，因此在后续的研究中需要进一步优化数据预处理环节。