基于可调Q因子小波变换的海杂波抑制算法

张俊玲, 董 玫, 陈伯孝

(西安电子科技大学雷达信号处理国家重点实验室, 陕西 西安 710071)

0 引 言

海面检测目标主要存在海面物理散射机制极其复杂、目标回波弱、弱目标模型难以建立、目标与海杂波类别非平衡(目标所占距离单元数远小于海杂波)等难点,因此海面目标检测是一项长期的挑战[1]。早期雷达分辨率较低,海面目标检测通常通过目标信号与海杂波之间的能量竞争完成。而在高分辨率雷达的检测下,传统的基于能量比较的检测方法面临着检测不出微弱目标的问题。针对海杂波背景下的弱目标检测问题,各国学者做了大量的研究工作,提出了许多检测方法,如基于海面分形特征的方法、基于海杂波混沌特征的方法、基于时频分析的方法、基于极化分析的方法、基于多特征联合的人工智能类方法等[1-2]。

海杂波背景下检测目标的传统方法是恒虚警率(constant false alarm rate, CFAR)检测算法[3]。该算法基于杂波均匀的假设,也就是认为每个距离单元的后向散射具有相同的统计分布。然而,实际处理的海杂波通常是非均匀的,这导致该算法的检测性能明显下降。

时频联合处理方法可以有效地描述信号频谱分量随时间的变化。常见的时频分析方法有短时傅里叶变换[4](short-time Fourier transform, STFT)、Wigner-Ville分布[5](Wigner-Ville distribution, WVD)和小波变换等。时频分析方法STFT和WVD由于自身固有的缺点,在抑制海杂波方面仍存在许多难点。小波变换具有良好的时频局部化特性,常常应用于海杂波背景下的微弱目标检测。小波变换中所用小波函数的选择不是唯一的,因此如何选择合适的小波函数是一个难点,而这往往需要通过大量的试验获得[6]。许多学者针对时频信息做了进一步处理,梁壮等[7]利用时频信息进行奇异值分解的海杂波抑制算法,通过提取前Q阶信号分量分离目标信号与海杂波,但是对于阶数的选取会对结果造成影响。左磊等[8]在时频域迭代分解过程中,利用凸包算法融合特征检测慢速弱目标,但仍然包含时频处理过程的难点。

小波变换在频域有一个恒定的相对带宽,称为等Q结构[9],因此小波变换不能根据信号的振荡属性调整Q值。可调Q因子小波变换(tunable Q-factor wavelet transform, TQWT)可以根据信号的振荡属性高低调整Q值,通过Q值的选取使可调Q因子小波对具有特定振荡属性的信号实现最优的匹配。对信号进行TQWT处理之后,在小波子带中呈现出稀疏特性,因此可以将稀疏表示理论[4]引入海杂波抑制算法的处理过程中。Selesnick[9]首先提出TQWT,通过选择品质因子Q、过采样率r和分解层数J匹配具有不同振荡属性的信号,从而实现信号的分离,适用于存在一个目标的情况。Ng等[10]和Rosenberg等[11]采用基追踪去噪(basis pursuit denoising, BPD)的方法来优化信号的稀疏性,并通过Ingara实测数据验证该方法的有效性。Farshchian等[12]应用双可调Q小波变换将橡胶充气船和鸟群两种目标从海杂波中分离,使算法可以适应多个目标存在的情况,此算法得到进一步应用。潘美燕等[13]提出了基于能量选择的改进TQWT海杂波抑制算法,该算法具有更好的检测性能,但是其参数的选取依赖目标信号和海杂波的先验信息。冯云等[14]针对文献[13]进行了改进,运用形态分量分析法分析得到稀疏系数,该方法不依赖目标信号和海杂波的先验信息,但是其选择合适参数的算法复杂度高,在实际应用中并不适用。

文献[9-14]均根据目标信号的振荡特性选取参数,重构目标信号,在目标检测步骤没有介绍具体的目标检测方法。考虑到目标信号能量弱,很难根据其振荡特性选择出最匹配的参数,而海杂波能量强,本文匹配杂波能获得更为准确的重构信息。同时,将原始回波信号与海杂波重构信号的差作为检测样本,采用自适应的阈值检测方法进行目标检测。由于利用了海杂波的振荡特性,本文提出改进的TQWT的海杂波抑制算法,不仅参数不依赖海杂波模型,而且目标检测中阈值具有自适应性。通过对某“雷达对海探测计划数据共享计划”中的实测海杂波数据进行处理,验证了本文所提出算法的有效性和实用性。

1 TQWT

目标信号具有一定径向速度,多普勒频谱窄,振荡特性强,而海杂波由于海浪的运动导致多普勒频谱展宽,振荡特性弱。基于此物理机理，利用不同Q因子的小波基函数对海杂波和目标信号进行匹配[13]。由于本文选择的Q因子与海杂波匹配,而与目标失配,因此在原始回波与重构海杂波信号的差中目标信号能量占比提高,可以通过判断结果来确定原始回波中是否存在目标信号。

1.1 模型建立

根据目标信号的有无,接收的原始回波信号Y有Y0和Y1两种情况:

(1)

式中:XT为目标信号；XC为海杂波信号；Xn为噪声分量；Y1为存在目标信号时接收的回波信号;Y0为不存在目标信号时接收的回波信号。由于不存在目标信号,模型中不存在XT分量。

1.2 TQWT

TQWT[15]类似于离散小波变换(discrete wavelet transform, DWT),具有由低通和高通滤波器组成的双通道滤波器组,如图1所示。

图1 双通道滤波器组Fig.1 Two-channel filter banks

在图1中，H0(ω)和H1(ω)分别是低通滤波器和高通滤波器的传输函数。LPS和HPS分别代表低通尺度和高通尺度。α和β分别对应低通滤波器和高通滤波器的尺度变换参数。为了确保小波变换过程中小波系数不会过度冗余,滤波器比例因子α和β应满足以下条件:0<β≤1,0<α<1。同时为了得到完整的重构信号,要求滤波器组进行严格采样,因此要求α+β>1[15]。信号通过x(n)分解滤波器组将其分解为多组小波系数v0(ω)和v1(ω),如图1(a)所示。

TQWT是由分解滤波器组和重构滤波器组的级联组成。通过多层双通道滤波器组将信号x(n)进行多层分解,如图2所示为多层小波分解的实现结构图。将低通小波子带C(j)作为第j+1层双通道滤波器组的输入进行分解,而高通小波子带W(j)作为第j层小波系数的输出,其中j为小波分解的层数。

图2 TQWT的多层分解结构图Fig.2 TQWT’s multi-layer decomposition structure diagram

信号重构是分解过程的逆过程,即TQWT反变换是TQWT变换的逆过程。因此,重构信号x(n)需要将经过稀疏优化后的低通滤波器的小波系数C(j)和高通滤波器的小波系数W(j)作为第j层的输入,进行第j层重构。通过多层如图1(b)所示的重构滤波器组,就可得到信号x(n)的重构信号y(n)。

1.3 TQWT参数

TQWT是由一系列小波基函数组成的,小波基函数由过采样率r、品质因子Q和分解层数J这3个参数控制。3个参数的不同组合可以使小波基函数呈现不同的形状和性质,从而影响TQWT小波系数的稀疏性[16]。

(1) 过采样率r:由上述分析可知,双通道滤波器组由α+β因子过采样。如果对双通道滤波器组的低通输出进行无穷迭代以实现小波变换,则由r=β/(1-α)因子实现过采样,即所有子带的过采样率总和为r。由文献[16]分析得到,过采样率r控制小波系数的冗余度,即如果信号x有N个采样点,则小波系数的点数为rN。但是当r≈1时,小波类似于sinc函数,频域过渡带很窄,导致时域响应受到限制。为了避免出现此问题,一般选择r≥3。

(2) 品质因子Q:Q表述了小波变换的可调性,决定了小波的“振荡”性质。信号的振荡特性越低,Q值越小,相反Q值越大。因此较小的Q值可以更好地检测静止目标,而较大的Q值可以更好地检测运动目标。在实际计算中,Q因子太大会导致计算量的增大,因此不完全匹配的Q值也可获得良好的性能。Q值的选取准则可参考文献[17]。一般情况下,Q值在1到8之间选取[10],目标静止时的Q值为1,目标振荡特性最高时Q值为8。

(3) 分解层数J:分解层数过多会造成系数的冗余,带来巨大的计算量。为了避免此种情况,通常应限制分解层数。输入信号的长度决定了分解层数,即需要保证小波长度不低于被分析的信号长度。结合文献[15]的分析,可以得到J的上限为

(2)

因此选择合适的TQWT参数能够从接收的回波数据中分离出重要的组成部分。换句话说,选定一组与海杂波最匹配的参数Q、r、J能够从接收的回波信号中准确地估计出海杂波分量。

1.4 稀疏优化

由于过采样,得到多组小波系数。为了得到最具有稀疏性的一组小波系数,采用BPD算法[18]对小波系数进行稀疏优化[19]。该方法试图以近似值x最小化所观察到的信号y的平方误差。其中,该平方误差受到稀疏性约束。这个问题的数学公式[18-19]为

(3)

影响BPD优化性能的一个重要参数是正则化参数λ。由于λ值的选取与杂波的能量有关,即杂波能量越大,其值也越大,故本文将λ设置为与小波系数的l2-范数成比例的向量,用来调整小波系数的加权和,即λ=‖w‖2。BPD是一个凸优化问题,文献[18]提出可以使用SALSA(state-action-reward-state-action)算法来解决此问题。本文将影响SALSA算法收敛速度的参数μ选择为2,并使用200次迭代进行优化,参数μ的值和迭代次数是通过大量实验得到的[19]。

为了验证海杂波经过TQWT之后具有稀疏性,将实测海杂波数据、仿真目标信号和加入仿真目标信号的海杂波数据分别做TQWT和BPD稀疏优化,观察各自分解到子带中的小波系数的稀疏性。实测海杂波数据选取某网站的雷达对海探测数据集中的某一个脉冲,添加仿真目标的参数如表1所示。将海杂波数据、仿真目标信号和加入仿真目标信号的海杂波数据分别分解到40个小波子带中,其中子带号从小到大对应频段由高到低。图3(a)为仿真目标信号的各子带能量分布,其高频部分子带能量较高,但整体能量分布比较均匀,各子带的能量占比相差不大都比较低。由于选取的TQWT参数与仿真目标信号的振荡特性失配,因此目标信号的各子带不具有稀疏性。图3(b)为实测海杂波数据的各子带能量分布,其最后一个小波子带的能量占比明显高于其他子带,且低频部分子带的能量比较大。由于选取的参数与海杂波的振荡特性相匹配,因此海杂波信号的各子带具有稀疏性的特征。图3(c)为加入仿真目标信号的海杂波数据的各子带能量分布,由于海杂波占主要能量部分,同时选取的TQWT参数与海杂波相匹配,因此加入仿真目标信号的海杂波数据的各子带能量分布与海杂波数据的各子带能量分布几乎相同,都具有稀疏性的特征。

表1 仿真目标参数

图3 稀疏优化后的子带能量分布Fig.3 Subband energy distribution after sparse optimization

1.5 目标检测

目标检测的过程可以写为

H0yγH1

式中:γ为检测阈值;y是检测样本;H0对应接收回波信号为Y0的情况;H1对应接收回波信号为Y1的情况。由于最终的阈值检测依据能量,而λmax与海杂波能量有关[20-21],本文提出一种阈值检测的准则,即γ=|20lg(kλmax)|,λmax由所选的品质因子Q、过采样率r和分解层数J这3个参数构成的小波系数的二范数向量最大值确定,即λmax=(‖w‖2)max。

为了确定k的大小,将仿真目标信号加入实测海杂波数据,经过TQWT算法处理并做阈值检测,观察不同的k值确定的阈值检测结果。实测海杂波数据选取某网站的雷达对海探测数据中的某一个脉冲,添加仿真目标的参数如表2所示。

推荐理由:人的成长从身体开始，经历情绪、感觉上升到心理、认知，然后由精神升华——这是人内在的不同存在层面。这些奇妙的内在部分将协助儿童自己创造出一个独一无二的自我。而我们成人已经失去了或未曾珍惜创造完整自我的机会，我们眼生羡慕地看着儿童的成长，并由此反思自己。儿童的成长令我们啧啧称奇，也让我们从中重获完整成长的机会。

表2 目标参数

图4(a)为虚警概率一定时,检测概率随k值变化的曲线,图4(b)为检测概率一定时,虚警概率随k值变化的曲线,图4(c)为某一脉冲的原始回波信号,图4(d)为经过TQWT海波抑制算法处理后得到的信号。

图4 阈值的确定Fig.4 Determination of threshold

观察图4(a)和图4(b)发现,k<1时,检测概率和虚警概率都在0附近;1≤k≤2时,检测概率为1,但是其虚警概率较高;k>2时,虽然虚警概率变化不大,但是检测概率急速下降。只有k=1时检测概率接近1,同时虚警概率保持在10-4以下。根据奈曼-皮尔逊准则的思想,k值应该取1,即阈值为γ=|20lgλmax|。观察图4(c)和图4(d),经过TQWT海杂波抑制算法处理,第200和第1 500距离单元的目标信号清晰可见。选择不同的阈值进行阈值检测,其中阈值1为γ=|20lg(0.5λmax)|,阈值2为γ=|20lgλmax|,阈值3为λ=|20lg(2λmax)|。当阈值为阈值1和阈值3时,由于阈值过高,不能将第200和第1 500距离单元的目标信号全部检测到,造成漏检。当选择阈值2时,可以检测到在第200和第1 500距离单元的目标信号,同时不会将噪声判定为目标信号,进一步验证了k=1的正确性。

1.6 算法实现及框图

TQWT算法的具体实现步骤如下所示，流程框图如图5所示。

步骤 1根据海杂波的多普勒运动特性,选择TQWT参数,即品质因子Q、过采样率r和分解层数J,进行TQWT,将回波信号分解到小波域[17];

步骤 4将雷达接收回波信号与海杂波重构信号的差R作为检测样本,与本文提出的自适应阈值λ=|20lgλmax|进行比较,判断是否存在目标信号。

图5 基于TQWT的海杂波抑制算法Fig.5 Sea clutter suppression algorithm based on TQWT

2 实验及分析

2.1 实验数据

2.2 实验结果及性能分析

将仿真的运动目标信号分别添加到20210106155330_01_staring数据包中的amplitude_complex_T2矩阵和20210105160634_01_staring数据包中的amplitude_complex_T2矩阵中,分别验证在高海况条件和低海况条件下的海杂波抑制算法的性能。添加仿真目标的参数如表3所示。

表3 实验目标参数

2.2.1 高海况条件下运动目标的检测

在20210106155330_01_staring数据包中的amplitude_complex_T2矩阵中添加目标编号1和目标编号2的运动目标,由于高海况条件下海杂波的多普勒频率高,振荡特性强,根据文献[17],选取TQWT的参数为Q=8、r=3、J=17,实验结果如图6～图9所示。

图6 实测数据杂波+仿真目标信号时域分析(Q=8、r=3、J=17)Fig.6 Time domain analysis of measured data clutter+ simulation target signal (Q=8、r=3、J=17)

图7 实测数据杂波+仿真目标信号时频域分析(Q=8、r=3、J=17)Fig.7 Time-frequency domain analysis of measured data clutter+simulation target signal (Q=8、r=3、J=17)

图8 实测数据杂波+仿真目标信号TQWT抑制时域分析(Q=8、r=3、J=17)Fig.8 Time domain analysis by TQWT suppression of measured data clutter+simulation target signal (Q=8、r=3、J=17)

图9 实测数据杂波+仿真目标信号TQWT抑制 时频域分析(Q=8、r=3、J=17)Fig.9 Time-frequency domain analysis by TQWT suppression measured data clutter+simulation target signal (Q=8、r=3、J=17)

图6为实测海杂波数据的第1 000个脉冲的时域图,可以看出在时域中两目标均被海杂波淹没。图7(a)和图7(b)分别为第200和第1 500距离单元的时频图。从图7中可以看出,目标信号的中心频率被海杂波的频率覆盖,目标信号被海杂波淹没,传统基于频域滤波的算法不能将目标信号检测出。由于海况为高海况,选取与其振荡特性相匹配的参数,提高在原始回波信号与重构海杂波信号的差中失配的目标信号能量占比。图8为经过TQWT算法的杂波抑制处理后第1 000个脉冲的时域图,可以看出在时域中第200和第1 500距离单元的目标均被检测出来,几乎没有虚警存在,只是目标的幅度被抑制。图9(a)和图9(b)分别为经过TQWT算法的杂波抑制处理后,第200和第1 500距离单元的时频图。可以清晰地观察到,两目标的多普勒频率分别为62 Hz和186 Hz,也就是分别以1 m/s和3 m/s的速度靠近雷达。因此,在高海况条件下的低信杂比情况中,TQWT算法能有效抑制海杂波,提高信杂比,提高雷达检测弱目标的能力。

2.2.2 低海况条件下运动目标的检测

在20210105160634_01_staring数据包中的amplitude_complex_T2矩阵中添加目标编号1和目标编号2的运动目标,由于低海况条件下海杂波的多普勒频率低,振荡特性弱,根据文献[17],选取TQWT的参数为Q=2、r=3、J=5,实验结果如图10～图13所示。

图10 实测数据杂波+仿真目标信号时域分析(Q=2、r=3、J=5)Fig.10 Time domain analysis of measured data clutter+ simulation target signal (Q=2、r=3、J=5)

图11 实测数据杂波+仿真目标信号时频域分析(Q=2、r=3、J=5)Fig.11 Time-frequency domain analysis of measured data clutter+simulation target signal (Q=2、r=3、J=5)

图12 实测数据杂波+仿真目标信号TQWT抑制时域分析(Q=2、r=3、J=5)Fig.12 Time domain analysis by TQWT suppression of measured data clutter+simulation target signal (Q=2、r=3、J=5)

图13 实测数据杂波+仿真目标信号TQWT抑制 时频域分析(Q=2、r=3、J=5)Fig.13 Time-frequency domain analysis by TQWT suppression of measured data clutter+simulation target signal (Q=2、r=3、J=5)

图10为实测海杂波数据的第1 000个脉冲的时域图,可以看出在时域中两目标均被海杂波淹没。图11(a)和图11(b)分别为第200和第1 500距离单元的时频图,目标的中心频率略高于海杂波的频率,但是直接将目标提取出来并不容易。由于海况为低海况,选取与其振荡特性相匹配的参数,提高在原始回波信号与重构海杂波信号的差中失配的目标信号能量占比。图12为经过TQWT算法的杂波抑制处理后第1 000个脉冲的时域图,同样可以看出在时域中第200和第1 500距离单元的目标均被检测出来,几乎没有虚警存在,只是目标的幅度被抑制。经过TQWT算法处理,图13(a)和图13(b)同样可以清晰地观察到,两目标的多普勒频率分别为62 Hz和186 Hz,也就是分别以1 m/s和3 m/s的速度靠近雷达。因此,在低海况条件下本文提出的算法也可检测出具有不同多普勒频率的弱运动目标,提高了雷达检测弱目标的能力。

2.2.3 算法对比分析

为了进一步说明本文TQWT算法的有效性,对实测海杂波数据分别进行CFAR检测[3]、自适应归一化匹配滤波器(adaptive normalized matched filter, ANMF)算法[24-25]处理、基于能量选择的TQWT算法[13]和本文基于杂波重构的TQWT算法处理,并将以上4种算法的处理结果进行对比。在虚警概率为10-4的条件下,对实测强海杂波数据分别做1 000次蒙特卡罗实验,实验结果如图14所示。

图14 算法检测概率对比Fig.14 Algorithms’ detection probability comparison

图14为4种检测算法性能的对比图,即强海杂波下信杂比与检测概率的关系。可知本文基于杂波重构TQWT算法的检测性能优于CFAR检测算法、ANMF处理算法和基于能量选择的TQWT算法。在信杂比为-10 dB的条件下,基于杂波重构的TQWT算法的检测概率可达到1,而此时基于能量选择的TQWT算法检测概率在0.3左右,CFAR检测和ANMF处理算法的检测概率接近于0。若要求0.9的检测概率,本文基于杂波重构的TQWT算法所需的信杂比为-19 dB,基于能量选择的TQWT算法所需信杂比为-5 dB,ANMF算法所需信杂比为5 dB,CFAR算法所需信杂比为12 dB。本文添加的目标为仿真目标,与真实目标相比,其与海杂波的相关性更弱,在重构海杂波时,能量损失更少。因此在处理真实目标时本文提出的改进的TQWT算法性能略有下降,需要进一步验证。

3 结 论

本文针对以海杂波为背景的弱目标检测中目标回波弱(即信杂比低)的问题,改进了基于TQWT的海杂波抑制算法。通过选取与海杂波振荡特性相匹配的参数,使经过TQWT后的小波子带呈现稀疏性,经过BPD做稀疏优化处理和TQWT反变换进行海杂波重构,最后将回波信号与海杂波重构信号的差作为检测样本,通过本文提出的自适应阈值检测准则判断目标信号是否存在。该算法不依赖于具体的海杂波模型,并且可以对接收的回波信号进行实时处理。对高海况和低海况下的实测海杂波数据进行处理的结果表明,本文算法的检测性能优于CFAR检测算法、ANMF算法和基于能量选择的TQWT算法。由于本文算法时间复杂度稍大,在硬件条件允许的情况下,其性能将会优于传统算法。