【干扰抑制】基于稀疏低秩Hankel矩阵分解的FMCW雷达干扰抑制【附MATLAB代码】

文章来源：微信公众号 EW Frontier

摘要-本文研究了采用去调频接收机的调频连续波（FMCW）雷达系统的干扰抑制问题。在去线性调频操作之后，来自目标的散射信号导致拍频信号，即，当干涉导致啁啾样短脉冲时，复指数之和。利用有用信号和干扰之间的这些不同的时间和频率特征，干扰抑制被配制成一个优化问题：通过提升测量值构造的Hankel矩阵的稀疏和低秩分解。然后，提出了一种迭代优化算法来解决这个问题，利用交替方向乘法器（ADMM）计划。与现有方法相比，该方法不需要检测干扰，也提高了分离出的有用信号的估计精度。点目标的数值模拟和分布式目标的实验结果（即，雨滴）来演示和验证其性能。仿真结果表明，该方法对静止和运动目标干扰抑制具有较好的适用性。

索引条款-FMCW雷达，干扰抑制，稀疏性，低秩，汉克尔矩阵。

I.介绍

调频连续波（FMCW）雷达（汽车雷达、生命体征观察、智能灯等）的广泛应用引起了对他们相互干涉的严重关切。结果表明，FMCW雷达中的外部干扰会增加杂波的总体水平，可能导致虚目标的出现，掩盖弱目标，干扰目标回波。已经提出了许多方法来减轻干扰。然而，在给定的信号处理成本的干扰抑制的期望水平尚未实现。在文献中，FMCW雷达系统的干扰减轻（IM）通常通过两种方式来解决：(1)系统/信号（即，天线阵列、波形等）设计方法[1]-[10]和（2）信号处理方法[11]- [24]。与系统/信号设计相关的方法通常要求雷达系统具有敏捷波形调制、新雷达架构、频谱感测模块或协调单元等的能力;因此，雷达系统必须修改或重新设计，系统的复杂性将会增加。另一方面，信号处理方法通过后处理来解决干扰抑制问题，这可以很容易地用于现有的FMCW雷达系统。

本文所介绍的工作属于后一组福尔斯。雷达系统干扰抑制的信号处理方法可以分为三类：滤波方法[11]-[13]，[25]，干扰置零和重构方法[14]-[17]，以及信号分离方法[18]-[24]。滤波方法通常利用特定域（即，时间、频率或空间），并设计适当的滤波器来抑制干扰。它们对于某些固定或缓慢变化的干扰工作良好。对于复杂场景，干扰可能是瞬时的，或者在不同扫描之间快速变化。然后可以利用自适应滤波[13]，其使用参考输入来生成相关干扰分量以用于干扰减轻。然而，在实际应用中，很难找到合适的参考输入信号，因此，它的干扰抑制性能无法保证。干扰归零和重建方法是首先切除受污染的信号样本[14]，然后在切除区域[15]，[17]中重建信号。对于这些方法来说，干扰的精确检测和切出信号样本的精确重构是关键步骤。由于FMCW雷达系统低通滤波器输出的干扰信号和目标信号的幅值通常是未知的，这使得干扰的检测成为一个具有挑战性的问题。

如果干扰不能被精确检测，它们将被部分去除，或者有用的信号被过度消除[16]。此外，即使正确地检测到干扰，也会不可避免地将被干扰污染的目标信号切掉，导致目标信号的功率损失。为了解决这个问题，提出了基于模型的方法，通过使用布尔格方法[15]和矩阵束方法[17]来重建切除区域中的有用信号。然而，随着污染样本在全测量中所占比例的增加，这些方法重建信号的精度会降低.另一方面，信号分离方法通过利用干扰或有用信号在某些基（或域）中的稀疏性来解决干扰减轻，这避免了干扰的显式检测。例如，对于超宽带雷达系统，目标的有用信号可以表现为时间上的稀疏尖峰，而射频干扰（RFI）通常在频域中是稀疏的[19]。

对于FMCW雷达，差拍信号在频域中可以是稀疏的。相比之下，在去调频和抗混叠低通滤波操作之后的干扰通常表现为类似啁啾的信号，其在各种场景中在时域、频域或时频域中可能是稀疏的。为了利用频域中的差拍信号和时频域中的干扰的稀疏表示，离散傅里叶变换（DFT）基和短时傅里叶变换（STFT）基是对应域中的规则网格，一般采用分裂增广拉格朗日收缩算法（SALSA）来解决相关的形态成分分析问题[18].虽然这些表示对于通过利用快速傅立叶变换（FFT）来有效实现基于SALSA的干扰减轻方法是有吸引力的，但是固有的所谓“离网”问题，即，真实频率（或时间-频率）分量与离散基之间的失配可能导致较少的稀疏表示，并因此降低干扰减轻性能。最近，瞬态干扰或RFI抑制也已经通过公式化为鲁棒主成分分析（RPCA）问题来解决[26]，[27]，其中由测量构建的矩阵被分解为稀疏矩阵和低秩矩阵的和[20]-[23]。通过使用奇异值阈值（SVT）算法[20]，[22]，重新加权核范数或重新加权Frobenius范数[23]来解决公式化的RPCA问题。在这些方法中，操作直接涉及用信号构造的矩阵的奇异值。因此，在这些方法中需要奇异值分解（SVD），这通常是非常昂贵的计算，特别是对于大型矩阵。

此外，由于用于合成孔径雷达系统的瞬态干扰减轻的多个脉冲重复间隔（PRI）中的测量所形成的信号矩阵通常不是结构化的，所开发的方法几乎不利用矩阵的结构。为了克服传统信号分离方法的“离网”问题和现有基于RPCA的方法计算量大的问题，提出了一种新的基于SPArse和低秩HanKeL矩阵合成的FMCW雷达干扰抑制方法（IM-SPARKLE）.对于所提出的方法，我们制定了FMCW雷达的干扰缓解作为一个RPCA类的问题，通过利用干扰的时间稀疏性和有用信号的频谱稀疏性。受用于指数分量估计的矩阵束方法的启发[17]，[28]，通过最小化用时域中的样本构造的汉克尔矩阵的秩来利用有用信号的频谱稀疏性。

秩最小化问题通常被放松为核为了规避用于核范数最小化的计算昂贵的SVD，我们利用Frobenius范数分解[29]-[31]来放松核范数最小化;因此，与传统的基于RPCA的IM方法相比，开发了更有效的算法。此外，我们的配方直接施加的sparsiy干扰的时域样本，并避免了不均匀的加权效果时，不同的样本使用的干扰分量的汉克尔矩阵。因此，自然不需要像[23]中那样使用重新加权操作。本文的其余部分组织如下。在第二节中，提出了作为信号分离问题的干扰减轻的问题公式，这导致Hankel矩阵的稀疏和低秩分解的优化问题。然后，在第三节中详细地提供了解决这个问题的所提出的算法。之后，第IV节和第V节显示了数值结果和实验测量，以证明所提出的IM方法用于FMCW雷达系统的性能。最后，在第六节中得出结论。

算法流程

MATLAB代码：

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Matrix completion                                    %
% Factorization nuclear norm to sum of Frobenious norm %
% Simulation: Three stationary targets                 %
% Edit by J.Wang, May 13, 2020                         %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
close all;clear; clc;
addpath('./lowRaS_MC/')

path4figure = './figs/';
flag_plot=0;
%%  load data
load('./data/Data4Demo.mat')

t_us = t*1e6;   % change the unit of time to micro-second.
%% ======================Interference Mitigation =========================
%% =============IM-SPARKLE==========================
disp('IM-SPARKLE...')
beta_1 = 0.1;           % 0.5
mu     = 0.02;          % 0.05
tau    = 0.02;          % 0.02
k_beta = 1.6;
k_mu   = 1.2; 
R = 10;

tic;
[x_lowRaS, i_lowRaS, rerr] = lowRaS_Hankel(sig_full_trc, R, beta_1, mu, tau, k_beta, k_mu);
toc;
hf = figure
semilogy(1:length(rerr), rerr,'r-','linewidth',1);
xlabel('Iteration','FontSize', 11)
ylabel('Relative error','fontsize',11)
grid on
set(gcf,'PaperUnits','inches','PaperPosition',[0 0 4 3],'PaperSize',[4 3])
if flag_plot==1
print(hf,'-dpng', '-r300', [path4figure 'simu_SNR_' num2str(SNR)  '_converg.png']);
print(hf,'-dpdf', '-r300', [path4figure 'simu_SNR_' num2str(SNR) '_converg.pdf'],'-opengl');
saveas(hf, [path4figure 'simu_SNR_' num2str(SNR)  '_converg.fig'])
end
%% Evalutation metric

SINR_0 = 20*log10(norm(sig_Rx_trc)/norm(sig_Rx_trc- sig_full_trc) )

SINR_lowRaS = 20*log10(norm(sig_Rx_trc)/norm(sig_Rx_trc.'- x_lowRaS) )
corr_lowRaS = (x_lowRaS)'*(sig_Rx_trc.')/(norm(x_lowRaS) * norm(sig_Rx_trc))


%% Image display
hfig_1 =figure;
Ftsz = 11;
set(gcf,'PaperUnits','inches','PaperPosition',[0 0 4 3],'PaperSize',[4 3])
plot(t_us, real(sig_full_trc), 'b-')%t_us, real(sig_Rx_trc), 'b-',...
hold on
%x_rect1 = [70 104 104 70 70];       y_rect1 = [-7 -7 7 7 -7];
x_rect1 = [101 152 152 101 101];       y_rect1 = [-7 -7 7 7 -7];
x_rect2 = [185 245 245 185 185];    y_rect2 = [-7.5 -7.5 7.5 7.5 -7.5];
x_rect3 = [280 320 320 280 280];    y_rect3 = [-6 -6 6 6 -6];
plot(x_rect1, y_rect1,'r--',...
     x_rect2, y_rect2,'r--',...
     x_rect3, y_rect3,'r--')
text(x_rect2(3),y_rect2(3)-0.2, '\leftarrow Interference','Color','red','Fontsize', Ftsz)
grid on
axis tight
ylim([-8 8])
xlabel('Time [\mus]', 'FontSize', Ftsz)
ylabel('Amplitude', 'FontSize', Ftsz)
title('Real part', 'FontSize', Ftsz)
%legend('reference', 'Contaminated Sig')
if flag_plot==1
print(hfig_1,'-dpng', '-r300', [path4figure 'simu_sig_full_SNR_' num2str(SNR)  '.png']);
print(hfig_1,'-dpdf', '-r300', [path4figure 'simu_sig_full_SNR_' num2str(SNR) '.pdf'],'-opengl');
saveas(hfig_1, [path4figure 'simu_sig_full_SNR_' num2str(SNR)  '.fig'])
end

hfig_2 =figure;
set(gcf,'PaperUnits','inches','PaperPosition',[0 0 6 4],'PaperSize',[6 4])
subplot(211)
plot(t_us, real(sig_Rx_trc),'-r',...
     t_us, real(x_lowRaS),'g-.')
% % hold on
% % x_rect_inset = [100 110 110 100 100]; y_rect_inset = [-2.7 -2.7 2.7 2.7 -2.7];
% % plot(x_rect_inset,y_rect_inset,'r--')
grid on
% axis tight
ylim([-1.8, 1.8])
legend('ref sig','Proposed'); %,'Location','south','Orientation','horizontal'
xlabel('Time [\mus]', 'FontSize', Ftsz)
ylabel('Amplitude', 'FontSize', Ftsz)
title('Recovered signal', 'FontSize', Ftsz)

%axes('Position',[0.45 0.78 0.45 0.12])
subplot(212)
II = t_us >=200 & t_us<=206;
% II = t_us >=250 & t_us<=255;
plot(t_us(II), real(sig_Rx_trc(II)),'r-',...
     t_us(II), real(x_lowRaS(II)), 'g-.')
grid on
axis tight
xlabel('Time [\mus]', 'Fontsize', Ftsz)
ylabel('Amplitude', 'Fontsize', Ftsz)
title('Close-up of the recovered signal')
if flag_plot==1
print(hfig_2,'-dpng', '-r300', [path4figure 'simu_extrac_usable_sig_SNR_' num2str(SNR) '.png']);
print(hfig_2,'-dpdf', '-r300', [path4figure 'simu_extrac_usable_sig_SNR_' num2str(SNR) '.pdf'],'-opengl');
saveas(hfig_2, [path4figure, 'simu_extrac_usable_sig_SNR_' num2str(SNR) '.fig'])
end
%%
hfig_3 =figure;
set(gcf,'PaperUnits','inches','PaperPosition',[0 0 4 3],'PaperSize',[4 3])
plot(t_us,real(x_lowRaS),'b--')
grid on
xlabel('Time [\mus]', 'FontSize', Ftsz)
ylabel('Amplitude', 'FontSize', Ftsz)
title('Extracted usable signal', 'FontSize', Ftsz)
if flag_plot==1
print(hfig_3,'-dpng', '-r300', [path4figure 'simu_usable_sig.png']);
print(hfig_3,'-dpdf', '-r300', [path4figure 'simu_usable_sig.pdf'],'-opengl');
end

%% Range profile
sig_len = length(sig_Rx_trc);
Num_fft_RP = 2^(nextpow2(sig_len)+1);

d = (0:Num_fft_RP-1)/Num_fft_RP * f_s /sweep_slope * c/2; 

RP_ref = ifft(sig_Rx_trc, Num_fft_RP);
RP_sig_full = ifft(sig_full_trc, Num_fft_RP);
RP_lowRaS = ifft(x_lowRaS, Num_fft_RP);



RP_max = max(abs([RP_ref, RP_sig_full, RP_lowRaS.']));

RP_ref_nor = db(abs(RP_ref)/RP_max);
RP_sig_full_nor = db(abs(RP_sig_full)/RP_max);
RP_lowRaS_nor = db(abs(RP_lowRaS)/RP_max);

Num_fft_disp = floor(Num_fft_RP/2);

d_1 = d(1:Num_fft_disp)/1e3;
hfig_5 = figure
plot(d_1, RP_ref_nor(1:Num_fft_disp), 'r-',...
     d_1, RP_sig_full_nor(1:Num_fft_disp), 'b--',...
     d_1, RP_lowRaS_nor(1:Num_fft_disp), 'g--')
grid on
axis([0 8 -60 2])
xlabel('Range [km]', 'Fontsize', Ftsz)
ylabel('Amplitude [dB]', 'Fontsize', Ftsz);
title('Range profiles', 'Fontsize', Ftsz);
legend('ref', 'sig\_Interf', 'Proposed', 'Location','southeast')

axes('Position', [0.68 0.72 0.2 0.18])
Ind = d_1>3.494 & d_1<3.506;
plot(d_1(Ind), RP_ref_nor(Ind), 'r-',...
     d_1(Ind), RP_sig_full_nor(Ind), 'b--',...
     d_1(Ind), RP_lowRaS_nor(Ind), 'g--')
grid on


axes('Position', [0.38 0.72 0.2 0.18])
Ind = d_1>1.994 & d_1<2.006;
plot(d_1(Ind), RP_ref_nor(Ind), 'r-',...
     d_1(Ind), RP_sig_full_nor(Ind), 'b--',...
     d_1(Ind), RP_lowRaS_nor(Ind), 'g--')
grid on

set(gcf, 'PaperUnits','inches', 'PaperSize', [4 3], 'PaperPosition', [0 0 4 3])
if flag_plot==1
    print(hfig_5, '-dpng', '-r300', [path4figure 'simu_PtTar_RP_SNR_' num2str(SNR) '.png']);
    print(hfig_5, '-dpdf', '-r300', [path4figure 'simu_PtTar_RP_SNR_' num2str(SNR) '.pdf']);
    saveas(hfig_5, [path4figure 'simu_PtTar_RP_SNR_' num2str(SNR) '.fig'])
end

%%

仿真结果

主要作者：Jianping Wang and Alexander Yarovoy are with the Faculty of ElectricalEngineering, Mathematics and Computer Science, Delft University of Technology

Digital Object Identifier 10.1109/TSP.2022.3147863