Damaged Radar Signal Reconstruction Method Based on Variational Mode Decomposition and Compressed Sensing

1.   引言

对雷达信号的准确分析处理是雷达对抗的核心任务之一。现有侦察采集设备特别是升空设备在信号实时采集过程中，受一系列弱观测因素(如低信噪比、信号传输和存储短猝故障、平台姿态断续对准目标等)的影响，导致部分时间内采集数据发生连续丢失或者随机丢失，同时伴有噪声干扰，严重影响了采集信号的完整性和准确性。因此，在进行雷达信号的分析与处理之前，必须进行合理的降噪与重构处理。如果能够去除实际采集数据中的噪声干扰，修复采集过程中丢失的数据，则将有效提高雷达对抗侦察信号处理结果的正确性和可靠性。

针对雷达信号的降噪问题，学者先后提出了短时傅里叶变换(Short-Time Fourier Transform, STFT)、卡尔曼滤波(Kalman Filtering, KF)、小波域变换(Wavelet Transform, WT)^[1]、经验模态分解(Empirical Mode Decomposition, EMD)^[2]等方法。雷达信号一般为非平稳信号，但STFT不能完全表征非平稳信号的时间频率变换特性，故这类方法存在一定的局限性。KF虽然可以处理非平稳信号，但在处理缺失信号或快速变化信号时，该类方法误差较大^[3]。WT具有良好的时频局部特性，但由于其基函数难以自适应信号变化的特点^[4]，这类方法亦有不足。EMD具有良好的自适应性，可以极好地反映信号局部频率特征，但存在端点效应、模态混叠等问题^[5]。2014年，文献[6,7]提出的变分模态分解(Variational Mode Decomposition, VMD)方法自适应效果好，可处理非线性、非平稳信号。该类方法可克服EMD的模态混叠问题，运算效率高，降噪能力强。目前，VMD已应用于激光雷达信号、监测信号、地震勘探信号、语音信号等领域的降噪处理^[8-11]。

对于弱观测条件下雷达信号重构问题，即采集过程中丢失数据修复问题的传统解决方法主要有插值法和曲线拟合法等。该类方法首先基于已有数据信息提取数据的分布，然后根据分布估计缺失的数据，然而该类方法需要求解高维方程组，随着数据量的增大，求解过程会越来越困难。此外，该类方法并未结合雷达信号的数据特点，因此其重构信号与原始信号相差较大，难以满足信号精细分析处理需求，影响雷达信号脉冲宽度、重复频率、脉内调制类型、码元宽度及编码规律等重要特征的判断和参数估计。压缩感知理论(Compressed Sensing, CS)将原始采集数据投影到只有少量测量值的稀疏空间，通过采用适当的重构算法得到完整信号。CS理论目前已经在诸多领域得到研究和应用，如图像处理领域，龚忠毅^[12]研究了基于子带自适应稀疏表示的图像压缩感知重构方法；医疗信号监测领域，齐林等人^[13]提出了基于CS的心电信号稀疏采样和重构方法；雷达领域，雷励^[14]基于CS理论设计了MIMO雷达波形，蒋莹等人^[15]提出的分布式压缩感知算法有效地解决了基于稀疏阵列的宽带欠定信号到达角估计问题；数据修复领域，张新鹏等人^[16]提出一种基于压缩感知原理的振动数据修复方法，另外在航天^[17,18]、超宽带雷达^[19]等领域也都开展了CS 的应用研究工作。

本文在上述研究的基础上，将VMD算法和CS原理引入弱观测条件下雷达信号数据修复问题中，提出了一种基于VMD-CS的雷达信号数据重构方法。该方法首先利用VMD算法对采集信号进行降解去噪处理，然后根据VMD算法的处理结果构造CS框架下的观测矩阵，进而针对雷达信号的特点，建立稀疏表示字典矩阵，在此基础上利用正交追踪匹配(Orthogonal Matching Pursuit, OMP)算法重构完整的雷达信号^[20]。

2.   VMD-CS雷达信号数据重构

弱观测条件下，实际采集的雷达信号不仅在完整性上有损失，且伴有污染。因此，高精度重构原始雷达信号须在实际采集的残损信号中，一方面消除噪声干扰的影响，另一方面修复丢失的数据。因此本文提出了一种基于VMD-CS算法的雷达信号数据重构方法，如图1所示。假定采集信号是由残缺的实际信号叠加上噪声等干扰信号构成。首先使用VMD算法对采集信号进行降解去噪处理，进而使用CS模型处理降噪数据，高精度重构原始雷达信号。

图  1  VMD-CS雷达信号数据重构系统框图

下载: 全尺寸图片 幻灯片

3.   VMD降解去噪

VMD算法针对EMD算法存在的模态混叠、虚假分量和端点效应等问题作出改进。该方法首先假设每一个固有模态(Intrinsic Mode Function, IMF)分量有着不同的中心频率和有限带宽，然后采用交替方向乘子法(Alternate Direction Method of Multipliers, ADMM)，不断更新各模态及其中心频率，使得每一个IMF分量的估计带宽之和最小，逐步将各模态解调到相应的基频带，这样便实现了信号从低频到高频的有效分离。换言之，只需在VMD算法处理之后，选取合适的IMF即可从采集信号中去除噪声等干扰信号的影响。

假定将输入信号f(t)分解为N个IMF分量，则相应的变分问题构造过程如下：

(1) 对每个IMF分量进行Hilbert -Huang变换得到其解析信号

其中，$\delta \left(t\right)$为狄利克来函数，${u}_{k}\left(t\right)$为IMF函数。

(2) 对各解析信号预估其中心频率，将各信号的频谱变换到基带上，即

其中，${w}_{k}$为IMF的中心频率

(3) 计算式(2)解调信号的欧氏距离，估计各模态带宽。得到变分约束问题

为计算该变分约束模型的最优解，需将该约束问题转换成无约束问题。因此，引入二次罚函数项$\alpha$和Lagrange乘数$\lambda$算子进行约束。

(4) ADMM算法迭代，得到最终的IMF分量。

由于VMD算法对信号进行模态分解时，雷达信号中的随机噪声将会被分离到较高模态中，通过去除较高频率模态分量，便能抑制信号中的噪声。另一方面，本文同时使用了非相关模态去除原则，既可以抑制噪声的干扰，也可以消除实际雷达信号中混入的非相关无意干扰信号，处理结果如图2所示。可见通过VDM算法处理之后的降噪信号无论在时域上还是频域上都得到了一定改善，信号中的噪声部分被较好地分离出来了。

图  2  雷达信号VMD降噪结果

下载: 全尺寸图片 幻灯片

4.   CS重构

雷达信号采集数据可能存在连续丢失或者随机丢失，通过标签索引的方式标记所有数据，即标签1表示数据未丢失，标签0表示数据丢失，构造采集信号的索引集合。基于索引集合和单位矩阵${{{I}}_n} \in {R^{N \times N}}$可构造出残缺雷达信号的观测矩阵$\varPhi \in {R^{N \times N}}$。考虑到CS的数据恢复原理是将数据压缩过程和采样过程合并，将原始采集数据投影到只有少量测量值的稀疏空间，再采用适当的重构算法由测量值重构得到完整信号。而在稀疏空间，雷达信号数据丢失概率相对较小，通过稀疏空间的信号重构可以大概率恢复出原始雷达信号。因而，本文引入CS方法重构雷达信号，同时考虑到离散余弦变换(Discrete Cosine Transform, DCT)适用于雷达信号数据的稀疏表示，故基于DCT变换构造稀疏字典矩阵，再通过OMP算法重构雷达信号。综上所述，可将残缺雷达信号重构算法流程如图3所示，具体步骤如下：

图  3  弱观测条件下雷达信号重构方法流程图

下载: 全尺寸图片 幻灯片

(1) 观测得到采集信号$y \in {R^M}$(M<N)以及索引集合Z；

(2) 设${{{I}}_n} \in {{{R}}^{N \times N}}$为单位矩阵，将I_n中的第j行($j\in Z$)删除，其它行保持顺序和大小不变，得到观测矩阵$\varPhi \in {R^{N \times N}}$；

(3) 利用DCT变换构造稀疏字典矩阵$C \in {R^{N \times N}}$；

(4) 使用OMP算法重构稀疏向量$\varTheta$得到其估计值$\hat \varTheta$；

(5) 利用逆DCT变换得到修复信号$\hat x = {C^{ - 1}}\hat \varTheta$。

4.1   信号重构原理

设观测得到的采集信号为${{{y}}\in {{R}}}^{{M}}$(M<N)，其中M表示实际采集得到的数据点数，N表示无丢失情况下的原始数据点数，即：采集设备在实际采集过程中丢弃了原始完整信号不是特别重要的N–M点数据。令这丢失的N–M点数据的数值为零(因为该N–M点数据对应时刻未能成功采集到信号)，因此可由实际采集的M点数据和这补充的N–M个零数据，重新构建一组N点数据，记为x₀，其在数据丢失位置对应的数据值为零。将丢失位置的索引记为集合Z，则将x₀中对应集合Z索引位置的N–M个零值去掉，保持其它元素值顺序和大小不变，可得到观测信号${{{y}}\in {{R}}}^{M}$(M<N)。

设${{{I}}_n} \in {{{R}}^{N \times N}}$为单位矩阵，将${{I}}_{{n}}$中的第j行($j\in {{Z}}$)删除，其它行保持顺序和大小不变，可得到观测矩阵$\varPhi \in {{{R}}^{N \times N}}$，进而应用CS修复算法，可以重构原始信号N点数据的估计值。

根据上述定义和假设，可得

其中，$\varPhi$为观测矩阵。

原始信号x可以通过稀疏表示字典矩阵$\varPsi \in {{{R}}^{N \times N}}$进行稀疏化，即

其中，$\varTheta$为稀疏表示系数向量。

令$A = \varPhi \varPsi$,则

此时已知实际观测值(采集信号)y和测量矩阵A，需要重构稀疏系数向量$\varTheta$，如果$\varTheta$满足稀疏便可成为CS框架下的信号重构过程。

有学者证明，只需测量矩阵A满足约束等距性质(Restricted Isometry Property, RIP)或者观测矩阵$\varPhi$和稀疏表示字典矩阵$\varPsi$不相关即可实现信号的重构^[21]。

RIP的定义如下：对任意k-稀疏的信号$\varTheta$和常数${\delta }_{k}\in \left({0,1}\right)$，满足

因此，只需要保证测量矩阵A满足RIP条件，或者观测矩阵$\varPhi$和稀疏表示字典矩阵$\varPsi$不相关，则可重构出稀疏向量$\varTheta$的估计值$\hat \varTheta$，进而可以计算完整采集信号x的估计值$\widehat{{x}}$，即

4.2   雷达信号修复

由上所述，雷达信号的修复过程主要与观测矩阵${{\varPhi}}$、稀疏字典矩阵${{\varPsi}}$以及稀疏表示系数向量${{\varTheta }}$有关。其中观测矩阵${{\varPhi }} \in {{{R}}^{N \times N}}$的构造方法前文已经详叙。为实现信号的重构，测量矩阵A必须满足RIP条件或者观测矩阵${{\varPhi}}$和稀疏表示字典矩阵${{\varPsi}}$不相关。即：只需要选择合适的字典矩阵${{\varPsi}}$即可实现信号修复。一般来说能够使原始采集信号足够稀疏的字典矩阵${{\varPsi}}$并不是惟一的，能够将采集信号稀疏化的字典矩阵均可用于信号重构中。由于雷达信号的数据经过DCT后对应的系数比较稀疏，因此本文采用DCT矩阵作为稀疏字典矩阵，设其为${{C}}\in {{{R}}}^{N{\times}N}$，矩阵C的第i行、j列的元素计算方法如式(10)所示。

稀疏表示向量${{\varTheta}}$为

基于压缩感知，重构稀疏向量${{\varTheta}}$的估计值$\hat \varTheta$，进而得到原始完整雷达信号x的估计值$\widehat{{x}}$。

4.3   OMP算法

在H表示的Hilbert空间中，有一组向量{I₁ I₂···I_m}构成字典矩阵，其每一个向量可以称为原子，其尺寸与被表示的信号y均相同，且是进行过归一化处理的两两正交单位向量，即$\left\| {{{{I}}_i}} \right\| = 1,$$< {{{I}}_i},{{{I}}_j} > = 0,(i \ne j)$。如果从该字典矩阵(完备的原子库)中，选择一个与信号y最为匹配的原子，满足：$\left|\left\langle {y,{I}_{{r}_{0}}} \right\rangle \right|={\rm{sup}}_{i\in ({1,2},··· ,k)}\left|\left\langle {y,{I}_{{r}_{i}}} \right\rangle \right|$, r₀表示一个字典矩阵的列索引。这样便构建了信号y的稀疏逼近，并进而求出逼近之后的残差，$y = \left\langle {y,{I_{{r_0}}}} \right\rangle {I_{{r_0}}} +$${R_1}f$。如此反复的选择信号残差的最佳匹配原子，这样信号y便可以分解成这些原子以及最终可以忽略的残差的线性组合。因此，OMP算法的核心步骤归纳如下：

(1)输入观测矩阵${\varPhi}$，采样信号${{y}}$及稀疏度K；

(2)初始化残差${R_t}f = y$，索引集${{\varLambda }}_{0}={\varnothing },t=0$；

(3)找到残差和观测矩阵各列中内积最大值时的脚标${\lambda }_{t}$，即：${\lambda _t} = \arg \max \left| {\left\langle {{R_t}f,{\varphi _i}} \right\rangle } \right|$；

(4)更新索引集${{\varLambda }}_{\rm{t}}={\varLambda }_{t-1}\cup {\lambda }_{t}$，并记录寻找到的观测矩阵重建原子集${\varPhi _t} = \left[ {{\varPhi _{t - 1}},{\varPhi _{{\lambda _t}}}} \right]$；

(5)由最小二乘法解得：$\hat \varTheta = \arg {\rm{min}}{\left\| {y - {\varPhi _t}{\varTheta _t}} \right\|_2}$

(6)更新残差${R_t}f = y - {\varPhi _t}{\hat \varTheta _t}$, t=t+1；

(7)判断是否满足$t > k$，若满足，则迭代终止，反之则重复步骤(3)—步骤(7)。

5.   仿真实验

为验证本文方法的有效性，本文针对常见雷达信号(线性调频信号LFM、非线性调频信号NLFM、Baker码信号PSK)进行对比实验，设置信号中心频率为52 MHz采集时长为2 μs，采样频率为500 MHz，采样点数为1000。对于线性调频信号设置其调频带宽为4 MHz且为正向调频，对于非线性调频信号采用正向偶二次方方式进行调频，调频带宽为4 MHz，对于Baker码信号选用11位编码方式。分别进行连续丢失和随机丢失两种情况下的对比实验，分别从时域、频域以及时频域角度对比了传统方法与本文方法的信号重构结果。为了进一步精确分析本文方法的重构效果，特别选用平均绝对误差进行衡量，如式(12)所示

5.1   连续丢失数据信号重构实验

针对原始完整信号连续丢失的情况，分别从信号时域、频域以及时频域角度对比了原始完整信号、连续丢失数据信号、本文方法重构结果以及传统方法重构结果，如图4所示(由于图片较多，连续丢失情况只显示了LFM信号修复结果，随机丢失选择NLFM信号显示)对于多脉冲信号，其脉内信号与单脉冲信号相似，因此对于多脉冲信号情况，可采用逐个脉冲依次重构的方法处理，然后再将各脉冲的重构结果拼接，从而得到多脉冲信号的重构结果。因此，多脉冲信号的实验结果及其产生原因与单脉冲情况相似，只是幅值上有所增大。

图  4  连续丢失数据条件下修复结果

下载: 全尺寸图片 幻灯片

分析图4可知，在连续丢失10%数据的情况下，传统的插值方法对于这种情况重构效果非常差，丢失区域的时域重构结果仍全为0，频域以及时频域上与原始完整信号区别更为明显，时频域上已经丢失了原始信号的调频特性。而本文的方法无论是对于LFM信号、NLFM信号还是PSK信号，都可以较好地重构出了原始信号，重构后信号在时域很接近，频域以及瞬时频率特性与原始完整信号基本相同，特别需要指出的是，在信号在时频域重构效果图中，虽然由于缺失点数达到了100，线条上在缺失区域对应位置有波动，调频特性有所损失，但是重构信号在总体趋势依然保留了原始完整信号的调频特性。

分析原因可能是：传统的插值算法重构丢失点信号时，需要利用丢失点的邻域信息进行插值计算，这就要求其邻域信息是存在且可靠的，然而在连续丢失情况下，丢失区域中间点的邻域点也发生了信息丢失情况，这显然不满足插值重构的前提条件。相反，本文方法原理是将信号变换到稀疏域进行重构修复，虽然信号在时域连续丢失了相对多的信息，但是当将其转换到稀疏域之后，有可能只丢失了整个信号相对较少的信息，或者是丢失了原始信号的部分不重要信息，因此本文方法的重构结果相较于传统方法明显有所改善。

观察图4(c)可以发现：当使用传统插值方法进行残缺数据修复重构时，重构结果中出现了较为突出的一段低频成分，这明显是不合理的。分析产生这段低频成分的原因可能是：当原始信号某处发生信息丢失时，直接反映在时域上是该处的采样点数值为零，在连续发生数据丢失的条件下，便形成了一段时间的零采样点情况。当采用传统的插值算法进行重构修复时，由于丢失区域中间点的邻域点也发生了信息丢失情况，导致这些丢失区域中间点的插值条件基本相同，进而使得这段区域的插值结果基本相同，且数值变化平缓，这种平缓变化的数据反映在频谱中便产生了一个较低的频率成分。由于LFM信号与PSK信号的频率变化相较于NLFM信号简单，因此低频混入情况也更为明显。以至于当使用传统插值法重构LFM信号与PSK信号时，往往会出现混入的低频成分幅值上超过原信号的主要频率成分的情况，这种情况随着丢失率的进一步上升越来越明显。

同样由于以上原因，在处理LFM, NLFM和PSK信号时，本文方法相较于传统插值算法明显更好的保留了原始信号的调频特性。本文方法在处理LFM与NLFM信号时，虽然在丢失区域存在一定的波动，但其波动的幅值较小且整个时间段上依旧保留了原信号的调频特性，而传统插值算法在该区域完全丢失了原信号的调频特性；处理PSK信号的时频结果显示本文方法完全修复了丢失区域频率损失，只是带来了一定的频率波动，如图4(f)所示。而传统插值算法却完全丢失了该区域，如图4(e)所示。此外需要指出，插值算法重构结果中，非丢失区域与原信号调频特性是相同，然而由于丢失区域瞬时频率变化幅值较大，导致非丢失区域的调频特性被掩盖了。

分析图4(a)可以发现：数据丢失率不高于30%时，本文方法重构出的采集信号与原始真实信号相比平均相对误差比较小，在丢失率超过40%后误差增加明显，有可能无法有效重构出原始采集信号。即本文方法在连续数据丢失率不超过40%时，均可以较为完好的重构原始信号。

5.2   随机丢失数据信号重构实验

针对原始完整信号随机丢失的情况，分别从信号时域、频域以及时频域角度对比了原始完整信号、随机丢失数据信号、本文方法重构结果以及传统方法重构结果，结果如图5所示(由于图片较多，随机丢失情况只绘制了NLFM信号修复结果，连续丢失选择LFM信号显示)。对于多脉冲信号，其处理方法与连续丢失情况时相同。

图  5  连续丢失数据条件下修复结果

下载: 全尺寸图片 幻灯片

分析图5可知，在随机丢失50%数据的情况下，传统插值方法的重构效果相对较差，时域波动明显，频域峰值点处虽然与原始信号相同，但是出现很多杂乱的邻域点，时频域与原始完整信号偏差很大，完全丢失了原始信号的线性调频特性。而本文的方法无论是对于LFM信号、NLFM信号还是PSK信号，都可以较好地重构出了原始完整信号，重构后信号在时域很接近原始信号，频域特性与原始完整信号基本相同，时频域上更是可以明显的体现出原始信号的调频特性。

相较于连续丢失情况，随机丢失情况下，传统方法的重构结果有所改善，但是当丢失数据进一步提高时，依然会出现局部小区域连续丢失的情况，因此其依然逃避不了其自身修复数据的局限性。相反，由于本文方法是将信号变换到稀疏域进行相应处理，数据点之间存在一定的信息关联性，在稀疏域中，随机丢失情况下丢失的重要信息概率相较于连续丢失情况降低了，因此，在随机丢失情况下，可以容忍的数据丢失率相较于连续情况更高。

观察图5(c)可以发现：相较于连续丢失情况，当使用传统插值方法对随机丢失信号进行残缺数据修复重构时，重构结果中并未有出现了连续丢失情况中的那段低频成分。分析原因可能是：由于随机丢失率较低的时候，当前丢失点附近很小概率会出其他丢失点，因此插值方法重构该点时所依据的信息还是准确的，而且由于每个丢失点所依据的邻域数据点都不同，所以不会出现连续丢失情况中的大段相同数据点，即不会出现异常的低频成分。然而随着随机丢失率的进一步提升，依然会出现局部小区域连续丢失的情况，所以这种异常的低频成分会随着随机丢失率的提升而出现且越来越明显。

此外，由于随机丢失数据点位置的采样值突变成零，这导致频谱中出现一个高频分量，这样一个数据点的丢失就会出现1～2次的高频分量情况，另一方面由于各丢失点之间的差异，因此导致丢失信号频谱中出现较多的毛刺频率分量。进一步观察发现经过本文方法重构修复之后这种毛刺频率分量也明显较少，使得重构频谱与原信号频谱更为接近。

同样由于以上原因，在处理LFM，NLFM和PSK信号时，本文方法相较于传统插值算法明显更好的保留了原始信号的调频特性。本文方法重构结果相较于原信号只是多了一些频率波动，但波动幅值不大，如图5(f)所示。而传统插值算法在整个时间段内完全丢失了原信号的调频特性，如图5(e)所示。

分析图5(a)可以发现：随机丢失情况下的信号重构效果要明显优于连续丢失情况下的信号重构效果，数据随机丢失率不高于60%时，本文方法的重构信号与原始真实信号相比平均相对误差接近于0，在丢失率超过65%后误差明显增大，有可能无法有效重构出原始采集信号。即本文方法在随机数据丢失率不超过60%时，可以较为完好的重构原始信号。

6.   结束语

本文将VMD算法与CS原理引入弱观测条件下雷达信号数据修复问题中，提出了一种基于VMD-CS的雷达信号数据重构方法。该方法首先利用VMD算法对采集信号进行降解降噪处理，然后根据VMD算法处理结果构造CS框架下的观测矩阵，接着针对雷达信号的特点，建立稀疏表示字典矩阵，在此基础上利用OMP算法重构采集信号。对比实验结果显示，无论是在连续丢失还是随机丢失情况下，本方法相较于现有技术可以更好地重构原始信号，其在时域、频域以及瞬时频率等方面与原始信号相当接近。