CN104007421B - 基于全变差和压缩感知的罗兰c被动雷达toa估计方法 - Google Patents

Abstract

本发明属于雷达目标无源定位技术领域，特别涉及基于全变差和压缩感知的罗兰C被动雷达TOA估计方法。该基于全变差和压缩感知的罗兰C被动雷达TOA估计方法包括以下步骤：对被动雷达接收的罗兰C直达波信号、罗兰C天波信号、以及罗兰C目标回波信号进行时域自适应滤波，得出时域自适应滤波后信号；对时域自适应滤波后信号进行全变差滤波，得出全变差滤波后信号；对全变差滤波后信号进行离散傅里叶变换，得出全变差滤波后频域信号；根据全变差滤波后频域信号，采用压缩感知方法重构出罗兰C目标回波时域稀疏信号；根据所述罗兰C目标回波时域稀疏信号，得出罗兰C目标回波的到达时间。

Description

基于全变差和压缩感知的罗兰C被动雷达TOA估计方法

技术领域

本发明属于雷达目标无源定位技术领域，特别涉及基于全变差和压缩感知的罗兰C被动雷达TOA估计方法。

背景技术

被动雷达由于自身不发射电磁波而对目标进行定位和跟踪，其自身具有不可探测性，比常规雷达具有更强的隐蔽性和生存能力，一直备受关注。目前，被动雷达利用广播电台或电视台、手机通讯基站，无线网络信号等非合作辐射源作为目标探测的照射源，通过接收来自辐射源的直达波和经由目标散射的回波，提取目标回波的波达方向(DOA,DirectionofArrival)、到达时间(TOA,TimeofArrival)、多普勒频移等信息，实现对目标的定位和跟踪。众所周知，作为被动雷达的广播电视信号沿视线传播，其有效目标定位作用范围从根本上受到视距限制。当目标位于视距以外时，基于广播电视信号的被动雷达将不能对此类目标进行探测定位，然而对视距以外的海面目标进行探测定位具有重要的军事和民用价值，因此在必要时，需要利用沿地波传播的具有超视距覆盖范围的罗兰C导航信号作为被动雷达的外部照射源信号，对目标回波TOA的估计方法进行研究。

一般被动雷达通过测量目标的波达方向并根据双基地距离信息对目标进行定位，然而该方案对于采用罗兰C被动雷达来说是不现实的，因为罗兰C信号的波长为3000米，如果要在接收端使用阵列对目标回波信号进行测向，则阵元间距需要1500米，这在实际中是难以实现的，而测量回波信号的到达时间则不受此限制，故以到达时间作为测量量对目标位置进行估计更具有现实意义。罗兰C被动雷达中微弱的目标回波信号淹没在强的直达波和噪声背景之中，因此很难对目标回波的到达时间进行有效的估计。为了估计目标回波信号的到达时间，现有方法包括匹配滤波法、广义互相关法和MUSIC超分辨时延估计方法。雷达系统中广泛使用的是匹配滤波法，但是该方法仅对大的时宽带宽积信号的到达时间有较好的估计效果，比如线性调频信号，步进频信号以及相位编码信号，而罗兰C信号的时宽带宽积较小，故采用匹配滤波法不能获得较好的到达时间估计效果。广义互相关法是对匹配滤波法的改进，该方法通过频谱加窗的方式降低匹配滤波后波形的副瓣，但是罗兰C信号脉冲压缩后主瓣很宽，故通过广义互相关法估计到达时间的效果仍然很差。MUSIC超分辨时延估计方法是借鉴空域波达方向估计的方式在使用频域样本来估计到达时间，由于MUSIC算法是子空间算法，事先要对信源个数进行估计以完成信号子空间和噪声子空间的正确划分，但是实践表明采用：MDL(minimumdescriptionlength)或BIC(Bayesianinformationcriterion)准则对罗兰C信源估计的个数是不准确的，采用MUSIC算法将造成大量的虚假目标到达时间。

发明内容

本发明的目的在于提出基于全变差和压缩感知的罗兰C被动雷达TOA估计方法。

为实现上述技术目的，本发明采用如下技术方案予以实现。

基于全变差和压缩感知的罗兰C被动雷达TOA估计方法包括以下步骤：

S1：利用被动雷达接收罗兰C直达波信号、罗兰C天波信号、以及罗兰C目标回波信号，然后对被动雷达接收的罗兰C直达波信号、罗兰C天波信号、以及罗兰C目标回波信号进行时域自适应滤波，得出时域自适应滤波后信号x_sur(n)，n为离散时间变量；

S2：对时域自适应滤波后信号x_sur(n)进行全变差滤波，得出全变差滤波后信号

S3：对全变差滤波后信号进行离散傅里叶变换，得出全变差滤波后频域信号k为离散频点变量；根据全变差滤波后频域信号采用压缩感知方法重构出罗兰C目标回波时域稀疏信号x；

S4：根据所述罗兰C目标回波时域稀疏信号x，得出罗兰C目标回波的到达时间。

本发明的特点和进一步改进在于：

在步骤S1中，利用被动雷达的主天线接收罗兰C直达波信号、罗兰C天波信号、以及罗兰C目标回波信号，利用被动雷达的辅助天线接收罗兰C直达波信号、以及罗兰C天波信号；根据被动雷达的辅助天线接收罗兰C直达波信号、以及罗兰C天波信号，通过时域自适应滤波方法主天线接收的罗兰C直达波信号以及罗兰C天波信号，得出时域自适应滤波后信号x_sur(n)。

在步骤S1中，被动雷达的主天线以脉冲组信号的接收方式接收罗兰C直达波信号、罗兰C天波信号、以及罗兰C目标回波信号；被动雷达的主天线对接收的多个脉冲组信号进行相干积累。

在步骤S1中，时域自适应滤波后信号x_sur(n)表示为：

$x_{\mathrm{sur}}\left(n\right)=s\left(n\right)+e_{\mathrm{sur}}\left(n\right)=\underset{i=1}{\overset{L}{\mathrm{\Σ}}}{\mathrm{\α}}_i{x}_g(n-{\mathrm{\τ}}_i)+e_{\mathrm{sur}}\left(n\right)$

其中，s(n)表示不含噪声的罗兰C目标回波信号，n为离散时间变量，n为整数且n取0至N-1，N表示时间长度，e_sur(n)表示被动雷达的主天线接收到的噪声信号分量；x_g(n-τ_i)表示第i个罗兰C目标回波信号，α_i表示第i个罗兰C目标回波信号的幅度，τ_i表示待估计的第i个罗兰C目标回波信号的到达时间，i取1至L，L表示罗兰C目标回波信号的个数；

在步骤S2中，对时域自适应滤波后信号x_sur(n)进行全变差滤波的过程为：求解以下优化问题：

$\hat{s}\left(n\right)=\underset{s\left(n\right)}{\arg \min }\underset{n=0}{\overset{N-1}{\mathrm{\Σ}}}{|x_{\mathrm{sur}}\left(n\right)-s\left(n\right)|}^2+\mathrm{\μ}\underset{n=1}{\overset{N-1}{\mathrm{\Σ}}}|s\left(n\right)-s(n-1)|$

其中，表示待求解的全变差滤波后信号，|·|表示取模，μ为设定的大于0的常数；

通过求解以上优化问题得出全变差滤波后信号

在步骤S2中，全变差滤波后信号表示为：

$\hat{s}\left(n\right)=s\left(n\right)+{\hat{e}}_{\mathrm{sur}}\left(n\right)=\underset{i=1}{\overset{L}{\mathrm{\Σ}}}{\mathrm{\α}}_i{x}_g(n-{\mathrm{\τ}}_i)+{\hat{e}}_{\mathrm{sur}}\left(n\right)$

其中，s(n)表示不含噪声的罗兰C目标回波信号，n为整数且n取0至N-1，N表示时间长度，x_g(n-τ_i)表示第i个罗兰C目标回波信号，α_i表示第i个罗兰C目标回波信号的幅度，τ_i表示第i个罗兰C目标回波信号的到达时间，i取1至L，L表示罗兰C目标回波信号的个数；表示经全变差滤波后被动雷达的主天线接收到的噪声信号分量；

在步骤S3中，对全变差滤波后信号进行离散傅里叶变换，得出全变差滤波后频域信号k为离散频点变量；全变差滤波后频域信号表示为：

$\hat{S}\left(k\right)=X_g\left(k\right)\underset{i=1}{\overset{L}{\mathrm{\Σ}}}{\mathrm{\α}}_i{e}^{-j\frac{2\mathrm{\π}}{M}k{\mathrm{\τ}}_i}+{\hat{E}}_{\mathrm{sur}}\left(k\right)$

其中，为经离散傅里叶变换后得出的频域信号，k为整数且k取1至M，M为离散傅里叶变换的点数；X_g(k)为x_g(n)经离散傅里叶变换后得出的频域信号，x_g(n)为：

x_g(n)＝ax₀(n)*G(n)

其中，

$G\left(n\right)=\underset{m=0}{\overset{7}{\mathrm{\Σ}}}[\mathrm{\δ}(n-1000m)\mathrm{pca}\left(m\right)+\mathrm{\δ}(n-I-1000m)\mathrm{pcb}\left(m\right)]$

其中，a为常数，x₀(n)表示罗兰C发射台的发射信号，A为设定的常数，*表示卷积；m取0至7，I表示以微秒为单位的罗兰C脉冲组信号的脉冲组重复周期，δ(·)表示罗兰C脉冲组信号的单位脉冲函数，pca(m)表示罗兰C脉冲组信号的奇数脉冲组重复周期内子脉冲的二相码，pcb(m)表示罗兰C脉冲组信号的偶数脉冲组重复周期内子脉冲的二相码；

x₀(n)经离散傅里叶变换后得出的频域信号表示为X₀(k)，将全变差滤波后频域信号除以X₀(k)，得出：

$F\left(k\right)=\hat{S}\left(k\right)/X_0\left(k\right)=\underset{i=1}{\overset{L}{\mathrm{\Σ}}}{\mathrm{\β}}_i{e}^{-j\frac{2\mathrm{\π}}{M}k{\mathrm{\τ}}_i}+{\hat{E}}_{\mathrm{sur}}\left(k\right)/X_0\left(k\right)$

其中，β_i＝aα_i；k为整数且k取1至M，M为离散傅里叶变换的点数；

然后，设置矩形窗函数W(k)，当90kHz≤k≤110kHz时，W(k)＝1；当k<90kHz或k>110kHz时，W(k)＝0；将F(k)加上矩形窗函数W(k)，得出：

$Y\left(k\right)=F\left(k\right)W\left(k\right)=\underset{i=1}{\overset{L}{\mathrm{\&Sigma;}}}{\mathrm{\&beta;}}_i{e}^{-j\frac{2\mathrm{\&pi;}}{M}k{\mathrm{\&tau;}}_i}+P\left(k\right)$

其中， $P\left(k\right)=W\left(k\right){\hat{E}}_{\mathrm{sur}}\left(k\right)/X_0\left(k\right);$

将Y(1)至Y(M)组成矢量y，y＝[Y(1),...,Y(k),...,Y(M)]^T，其中，T表示矩阵或向量的转置；则矢量y表示为：

y＝Φx+p

其中，p＝[P(1),...,P(k),...,P(M)]^T，x罗兰C目标回波时域稀疏信号，x为N×1维的稀疏矢量，Φ为M×N维的测量矩阵，Φ的矩阵展开形式为：

$\mathrm{\&Phi;}=\left[\begin{array}{cccccc}{e}^{-j\frac{2\mathrm{\&pi;}}{M}{\mathrm{\&tau;}}_i}& e^{-j\frac{2\mathrm{\&pi;}}{M}{\mathrm{\&tau;}}_2}& ...& e^{-j\frac{2\mathrm{\&pi;}}{M}{\mathrm{\&tau;}}_L}& ...& e^{-j\frac{2\mathrm{\&pi;}}{M}{\mathrm{\&tau;}}_N}\\ e^{-j\frac{2\mathrm{\&pi;}}{M}2{\mathrm{\&tau;}}_1}& e^{-j\frac{2\mathrm{\&pi;}}{M}2{\mathrm{\&tau;}}_2}& ...& e^{-j\frac{2\mathrm{\&pi;}}{M}2{\mathrm{\&tau;}}_L}& ...& e^{-j\frac{2\mathrm{\&pi;}}{M}2{\mathrm{\&tau;}}_N}\\ .& .& ...& .& ...& .\\ .& .& & .& & .\\ .& .& & .& & .\\ e^{-j\frac{2\mathrm{\&pi;}}{M}M{\mathrm{\&tau;}}_1}& e^{-j\frac{2\mathrm{\&pi;}}{M}M{\mathrm{\&tau;}}_2}& ...& e^{-j\frac{2\mathrm{\&pi;}}{M}M{\mathrm{\&tau;}}_L}& ...& e^{-j\frac{2\mathrm{\&pi;}}{M}M{\mathrm{\&tau;}}_N}\end{array}\right]$

其中，L<N；然后，求解以下l₀范数最优化问题：

min||x||₀

s.t.y＝Φx+p

其中，||x||₀表示x的l₀范数；通过求解以上l₀范数最优化问题重构出罗兰C目标回波时域稀疏信号x。

在步骤S4中，将所述l₀范数最优化问题替换为以下l₁范数最优化问题：

min||x||₁s.t.y＝Φx+p

其中，||·||₁表示求l₁范数。

在步骤S4中，在所述罗兰C目标回波时域稀疏信号x中，找出各个非零元素在x中的位置；根据每个非零元素在x中的位置、以及被动雷达的采样频率，得出对应的罗兰C目标回波的到达时间。

在步骤S4中，根据所述罗兰C目标回波时域稀疏信号x，对罗兰C目标回波进行恒虚警检测，得出经恒虚警检测后稀疏矢量，然后根据所述恒虚警检测后稀疏矢量，找出各个非零元素在恒虚警检测后稀疏矢量中的位置；根据每个非零元素在恒虚警检测后稀疏矢量中的位置、以及被动雷达的采样频率，得出对应的罗兰C目标回波的到达时间。

本发明的有益效果为：

(1)本发明利用沿地波传播的具有超视距覆盖范围的罗兰C导航信号作为被动雷达的外部照射源信号，有效目标定位作用范围相当可观,接收站自身具有不可探测性，比常规雷达具有更强的隐蔽性和生存能力，具有重要的军事和民用价值。

(2)本发明利用全变差滤波方法去除罗兰C频带内的色噪声,从而使得虚假目标回波到达时间估计对应的模值较小的尖峰全部消失，提升了算法的性能，既而通过恒虚警检测后，可以得到目标回波到达时间的精确估计。

(3)本发明利用压缩感知方法在频域对接收信号进行稀疏重构，通过测量目标回波到达时间的方法进行目标定位，解决了通过测量目标的波达方向对目标进行定位的传统方法在遇到罗兰C这种波长长达3000m的信号时无法实现的难题。

附图说明

图1为本发明的应用场景示意图；

图2为本发明的基于全变差和压缩感知的罗兰C被动雷达TOA估计方法的流程图；

图3a为本发明仿真实验1得出的重构结果示意图；

图3b为图3a的局部放大示意图；

图4a为本发明仿真实验2得出的重构结果示意图；

图4b为图4a的局部放大示意图；

图5为本发明仿真实验3得出的正确估计罗兰C目标回波的到达时间的成功概率与信噪比的关系曲线图；

图6为本发明仿真实验4得出的正确估计两个目标回波到达时间的成功概率与两个目标的回波到达时间差的变化关系示意图。

具体实施方式

下面结合附图对本发明作进一步说明：

参照图1，为本发明的应用场景示意图。在该应用场景中，包括罗兰C发射台、待检测目标、电离层、以及被动雷达。被动雷达接收的信号包括罗兰C直达波信号、罗兰C天波信号(经电离层发射的信号)、以及罗兰C目标回波信号。罗兰C发射台1向外发射信号，利用被动雷达的主天线3接收罗兰C直达波信号、罗兰C天波信号、以及罗兰C目标回波信号，利用被动雷达的辅助天线4接收罗兰C直达波信号、以及罗兰C天波信号，在图1中，2表示目标。

参照图2，为本发明的基于全变差和压缩感知的罗兰C被动雷达TOA估计方法的流程图。该基于全变差和压缩感知的罗兰C被动雷达TOA估计方法包括以下步骤：

S1：利用被动雷达接收罗兰C直达波信号、罗兰C天波信号、以及罗兰C目标回波信号，然后对被动雷达接收的罗兰C直达波信号、罗兰C天波信号、以及罗兰C目标回波信号进行时域自适应滤波，得出时域自适应滤波后信号x_sur(n)，n为离散时间变量。具体说明如下：

在步骤S1中，针对罗兰C发射台，天线接收到的单个目标散射的脉冲组重复周期(grouprepetitioninterval,GRI)罗兰C脉冲组信号x_g(n-τ)可以表示为：

x_g(n-τ)＝ax₀(n)*G(n-τ)

其中，

x₀(n)＝An²exp[-2n/65]sin(0.2πn)

$G\left(n\right)=\underset{m=0}{\overset{7}{\mathrm{\&Sigma;}}}[\mathrm{\&delta;}(n-1000m)\mathrm{pca}\left(m\right)+\mathrm{\&delta;}(n-I-1000m)\mathrm{pcb}\left(m\right)]$

其中，a为常数(与目标散射截面积相关的常数)，x₀(n)表示罗兰C发射台的发射信号，n为离散时间变量，其单位为微秒，*表示卷积；A为设定的常数(与信号峰值有关的常数)，m取0至7，I表示以微秒为单位的罗兰C脉冲组信号的脉冲组重复周期，δ(·)表示罗兰C脉冲组信号的单位脉冲函数，pca(m)表示罗兰C脉冲组信号的奇数脉冲组重复周期内子脉冲的二相码(取1或-1)，pcb(m)表示罗兰C脉冲组信号的偶数脉冲组重复周期内子脉冲的二相码(取1或-1)。

在步骤S1中，利用被动雷达的主天线接收罗兰C直达波信号、罗兰C天波信号、以及罗兰C目标回波信号，利用被动雷达的辅助天线接收罗兰C直达波信号、以及罗兰C天波信号。为了提升接收信号的信噪比，首先对被动雷达的主天线接收到的多个脉冲组信号进行相干积累，通过长时间的相干积累可以获得较大的相干积累增益(比如在基于南海罗兰C台链的被动雷达中，信号带宽20kHz，积累13.56s对应的相干积累增益为64dB)。

然后，通过时域自适应滤波方法主天线接收的罗兰C直达波信号以及罗兰C天波信号，得出时域自适应滤波后信号x_sur(n)。其具体过程为：由于被动雷达的主天线的接收信号中，直达波信号和天波信号的能量较强，目标回波信号能量较弱，因此首先根据被动雷达的辅助天线接收的罗兰C直达波信号以及罗兰C天波信号，通过时域自适应滤波的方法消除主天线接收的罗兰C直达波信号以及罗兰C天波信号。

经过相干积累和自适应滤波后，得出时域自适应滤波后信号x_sur(n)，x_sur(n)表示为：

$x_{\mathrm{sur}}\left(n\right)=s\left(n\right)+e_{\mathrm{sur}}\left(n\right)=\underset{i=1}{\overset{L}{\mathrm{\&Sigma;}}}{\mathrm{\&alpha;}}_i{x}_g(n-{\mathrm{\&tau;}}_i)+e_{\mathrm{sur}}\left(n\right)$

其中，s(n)表示不含噪声的罗兰C目标回波信号，n为离散时间变量其单位为微秒，n为整数且n取0至N-1，N表示时间长度(即罗兰C目标回波信号的数据长度)，e_sur(n)表示被动雷达的主天线接收到的噪声信号分量；x_g(n-τ_i)表示第i个罗兰C目标回波信号，α_i表示第i个罗兰C目标回波信号的幅度，τ_i表示第i个罗兰C目标回波信号的到达时间(TOA)，i取1至L，L表示罗兰C目标回波信号的个数。

S2：对时域自适应滤波后信号x_sur(n)进行全变差滤波，得出全变差滤波后信号具体说明如下：

实践中发现时域自适应滤波后信号中携带一些虚假TOA信息，这些虚假TOA信息是由于被动雷达的主天线接收到的噪声信号分量e_sur(n)在罗兰C频带内的色噪声的分量引起的扰动所致，然而此色噪声不容易用传统的滤波方法去除。此问题可在用时域自适应滤波后信号x_sur(n)做TOA估计之前先对x_sur(n)做全变差滤波解决，具体而言就是先求解以下优化问题得出全变差滤波后信号(为不含噪声的罗兰C目标回波信号s(n)的的波形估计)：

$\hat{s}\left(n\right)=\underset{s\left(n\right)}{\arg \min }\underset{n=0}{\overset{N-1}{\mathrm{\&Sigma;}}}{|x_{\mathrm{sur}}\left(n\right)-s\left(n\right)|}^2+\mathrm{\&mu;}\underset{n=1}{\overset{N-1}{\mathrm{\&Sigma;}}}|s\left(n\right)-s(n-1)|$

其中，表示待求解的全变差滤波后信号，|·|表示取模，μ为设定的大于0的常数；n为整数且n取0至N-1。

通过求解以上优化问题得出全变差滤波后信号实现对时域自适应滤波后信号x_sur(n)的全变差滤波。全变差滤波能够保证在消除噪声扰动的同时保证处理后的波形不产生相位移动和大幅度失真。

以上优化问题可以表示为如下矢量形式：

$\hat{s}=\underset{s}{\arg \min }{\left|\right|x_{\mathrm{sur}}-s\left|\right|}_2+\mathrm{\&mu;}{\left|\right|\mathrm{Ds}\left|\right|}_1$

其中D为(N-1)×N阶的矩阵。为全变差滤波后信号的矢量形式，x_sur为时域自适应滤波后信号x_sur(n)的矢量形式，s为不含噪声的罗兰C目标回波信号s(n)的矢量形式，||·||₂表示求l₂范数，||·||₁表示求l₁范数。该优化问题可利用增广拉格朗日算法求解。

全变差滤波之所以能滤除上述虚假TOA信息的原因有三：第一，由于噪声随机变换比信号快很多，在目标函数中令||Ds||₁最小可以抑制快速变换的噪声；第二，全变差滤波有保持信号转折点的特性，故滤波后波形失真较小；第三，全变差滤波没有FIR(finiteimpulseresponse)滤波或IIR(infiniteimpulseresponse)滤波产生的波形位置移动。

求解得出的全变差滤波后信号表示为：

$\hat{s}\left(n\right)=s\left(n\right)+{\hat{e}}_{\mathrm{sur}}\left(n\right)=\underset{i=1}{\overset{L}{\mathrm{\&Sigma;}}}{\mathrm{\&alpha;}}_i{x}_g(n-{\mathrm{\&tau;}}_i)+{\hat{e}}_{\mathrm{sur}}\left(n\right)$

其中，s(n)表示不含噪声的罗兰C目标回波信号，n为整数且n取0至N-1，N表示时间长度，x_g(n-τ_i)表示第i个罗兰C目标回波信号，α_i表示第i个罗兰C目标回波信号的幅度，τ_i表示第i个罗兰C目标回波信号的到达时间，i取1至L，L表示罗兰C目标回波信号的个数；表示经全变差滤波后被动雷达的主天线接收到的噪声信号分量。

S3：对全变差滤波后信号进行离散傅里叶变换，得出全变差滤波后频域信号k为离散频点变量；根据全变差滤波后频域信号采用压缩感知方法重构出罗兰C目标回波时域稀疏信号x。具体说明如下：

根据压缩感知理论，对于任一信号x，需要自身是稀疏的或者在某个基下进行稀疏表示(例如在稀疏基Ψ下可以稀疏表示为x＝Ψb，b是稀疏信号)时，才能使用压缩感知根据信号x的观测矢量y＝Φx恢复出稀疏信号x或x在基下的稀疏表示。因此首先需要找到表征目标回波TOA信息的稀疏信号x，然后才能使用压缩感知根据该信号的观测矢量y对x进行稀疏重构，以达到超分辨TOA估计的目的。

考虑到x_g(n-τ_i)＝x_g(n)*δ(n-τ_i),i＝1,2,...,L，由于时域的卷积相当于频域的乘积，故在频域中分析第i个罗兰C目标回波信号的到达时间τ_i的估计问题更方便。具体地说，首先对全变差滤波后信号进行离散傅里叶变换，得出全变差滤波后频域信号k为离散频点变量；全变差滤波后频域信号表示为：

$\hat{S}\left(k\right)=X_g\left(k\right)\underset{i=1}{\overset{L}{\mathrm{\&Sigma;}}}{\mathrm{\&alpha;}}_i{e}^{-j\frac{2\mathrm{\&pi;}}{M}k{\mathrm{\&tau;}}_i}+{\hat{E}}_{\mathrm{sur}}\left(k\right)$

其中，为经离散傅里叶变换后得出的频域信号，k为整数且k取1至M，M为离散傅里叶变换的点数；X_g(k)为x_g(n)经离散傅里叶变换后得出的频域信号，x_g(n)为：

x_g(n)＝ax₀(n)*G(n)

其中，

x₀(n)＝An²exp[-2n/65]sin(0.2πn)

$G\left(n\right)=\underset{m=0}{\overset{7}{\mathrm{\&Sigma;}}}[\mathrm{\&delta;}(n-1000m)\mathrm{pca}\left(m\right)+\mathrm{\&delta;}(n-I-1000m)\mathrm{pcb}\left(m\right)]$

其中，a为常数(与目标散射截面积相关的常数)，x₀(n)表示罗兰C发射台的发射信号，A为设定的常数(与信号峰值有关的常数)，*表示卷积；m取0至7，I表示以微秒为单位的罗兰C脉冲组信号的脉冲组重复周期，δ(·)表示罗兰C脉冲组信号的单位脉冲函数，pca(m)表示罗兰C脉冲组信号的奇数脉冲组重复周期内子脉冲的二相码(取1或-1)，pcb(m)表示罗兰C脉冲组信号的偶数脉冲组重复周期内子脉冲的二相码(取1或-1)。

x₀(n)经离散傅里叶变换后得出的频域信号表示为X₀(k)，将全变差滤波后频域信号除以X₀(k)，得出：

$F\left(k\right)=\hat{S}\left(k\right)/X_0\left(k\right)=\underset{i=1}{\overset{L}{\mathrm{\&Sigma;}}}{\mathrm{\&beta;}}_i{e}^{-j\frac{2\mathrm{\&pi;}}{M}k{\mathrm{\&tau;}}_i}+{\hat{E}}_{\mathrm{sur}}\left(k\right)/X_0\left(k\right)$

其中，β_i＝aα_i；k为整数且k取1至M，M为离散傅里叶变换的点数；表示对应的频域噪声。

由于罗兰C信号的能量集中在90kHz～110kHz的频带内，为了抑制此频带外的频域噪声，设置矩形窗函数W(k)，当90kHz≤k≤110kHz时，W(k)＝1；当k<90kHz或k>110kHz时，W(k)＝0；将F(k)加上矩形窗函数W(k)，得出：

$Y\left(k\right)=F\left(k\right)W\left(k\right)=\underset{i=1}{\overset{L}{\mathrm{\&Sigma;}}}{\mathrm{\&beta;}}_i{e}^{-j\frac{2\mathrm{\&pi;}}{M}k{\mathrm{\&tau;}}_i}+P\left(k\right)$

其中，Y(k)表示F(k)加上矩形窗函数W(k)后形成的频谱。

将Y(1)至Y(M)组成矢量y，y＝[Y(1),...,Y(k),...,Y(M)]^T，其中，T表示矩阵或向量的转置；则矢量y表示为：

y＝Φx+p

其中，p＝[P(1),...,P(k),...,P(M)]^T，x罗兰C目标回波时域稀疏信号，x为N×1维的稀疏矢量，Φ为M×N维的测量矩阵，Φ的矩阵展开形式为：

$\mathrm{\&Phi;}=\left[\begin{array}{cccccc}{e}^{-j\frac{2\mathrm{\&pi;}}{M}{\mathrm{\&tau;}}_i}& e^{-j\frac{2\mathrm{\&pi;}}{M}{\mathrm{\&tau;}}_2}& ...& e^{-j\frac{2\mathrm{\&pi;}}{M}{\mathrm{\&tau;}}_L}& ...& e^{-j\frac{2\mathrm{\&pi;}}{M}{\mathrm{\&tau;}}_N}\\ e^{-j\frac{2\mathrm{\&pi;}}{M}2{\mathrm{\&tau;}}_1}& e^{-j\frac{2\mathrm{\&pi;}}{M}2{\mathrm{\&tau;}}_2}& ...& e^{-j\frac{2\mathrm{\&pi;}}{M}2{\mathrm{\&tau;}}_L}& ...& e^{-j\frac{2\mathrm{\&pi;}}{M}2{\mathrm{\&tau;}}_N}\\ .& .& ...& .& ...& .\\ .& .& & .& & .\\ .& .& & .& & .\\ e^{-j\frac{2\mathrm{\&pi;}}{M}M{\mathrm{\&tau;}}_1}& e^{-j\frac{2\mathrm{\&pi;}}{M}M{\mathrm{\&tau;}}_2}& ...& e^{-j\frac{2\mathrm{\&pi;}}{M}M{\mathrm{\&tau;}}_L}& ...& e^{-j\frac{2\mathrm{\&pi;}}{M}M{\mathrm{\&tau;}}_N}\end{array}\right]$

其中，L<N；Φ的N列对应着一个脉冲重复周期内的N点等间隔采样点造成的相移，x为L稀疏的N×1维的矢量。x的一般项为β_j,j＝0,1,...,N-1，代表第j个罗兰C目标回波信号的复振幅。由于中观测矩阵Φ是M×N阶的离散傅里叶矩阵，已经有学者证明该矩阵高概率满足有限等距性质(RestrictedIsometryProperty,RIP)，只要y的长度M≥cLlog₂(N/L)即可以根据y以高概率不失真的重构出x，其中c一般取常数1即可。由于y不是实际测量而是通过离散傅里叶变换计算得到的，因此条件M≥cLlog₂(N/L)很容易满足。根据矢量y，x可以通过求解以下l₀范数最优化问题得到：

min||x||₀

s.t.y＝Φx+p

其中，||x||₀表示x的l₀范数；通过求解以上l₀范数最优化问题重构出罗兰C目标回波时域稀疏信号x。

虽然可以通过(OrthogonalMatchingPursuit,OMP)等贪婪算法求解以上l₀范数最优化问题，但是由于以上l₀范数最优化问题是一个非凸优化问题且在本问题中x的维度较高，使用贪婪算法效率非常低。有学者已经证明(参见Robustuncertaintyprinciples:exactsignalreconstructionfromhighlyincompletefrequencyinformation[J].InformationTheory,IEEETransactionson,2006,52(2):489-509)可以将以上l₀范数最优化问题转换为以下l₁范数最优化问题：

min||x||₁s.t.y＝Φx+p

其中，||·||₁表示求l₁范数，以上l₁范数最优化问题可以等价为以下简化形式的优化问题：

min||y-Φx||₂+λ||x||₁

其中，λ为大于0的常数。以上l₁范数最优化问题和简化形式的优化问题都是凸优化问题，本发明实施例中采用内点法、梯度投影方法和迭代阈值法去求解这两个凸优化问题。

S4：根据所述罗兰C目标回波时域稀疏信号x，得出罗兰C目标回波的到达时间。具体说明如下：

在所述罗兰C目标回波时域稀疏信号x，找出各个非零元素在x中的位置；根据每个非零元素在x中的位置、以及被动雷达的采样频率，得出对应的罗兰C目标回波的到达时间。具体地，如果罗兰C目标回波时域稀疏信号x中第i'行为非零元素，i'取1至N；则对应的第i'个罗兰C目标回波信号的到达时间τ_i'为：一个脉冲组重复周期内第i'个采样点的位置，根据采样频率f_s可以得出第i'个罗兰C目标回波信号的到达时间τ_i'：τ_i'＝i/f_s。

采用上述方法得出的罗兰C目标回波信号的到达时间会出现虚假TOA信号，因此，作为本发明实施例的一种改进，在在步骤S4中，根据所述罗兰C目标回波时域稀疏信号x，对罗兰C目标回波进行恒虚警检测(单元平均恒虚警检测)，得出经恒虚警检测后稀疏矢量，然后根据所述恒虚警检测后稀疏矢量，找出各个非零元素在恒虚警检测后稀疏矢量中的位置；根据每个非零元素在恒虚警检测后稀疏矢量中的位置、以及被动雷达的采样频率，得出对应的罗兰C目标回波的到达时间。

本发明的效果可以通过以下仿真实验进一步说明：

仿真条件：

为了验证本发明的正确性、有效性和稳定性，需要进行性能分析。在仿真实验中，直接对时域自适应滤波后信号进行处理。设定罗兰C脉冲组(GRI)重复周期为67800μs，经过相干积累后，单个脉冲重复周期的数据长度取1000点，采样率为1MHz。本仿真实验中，只对2个罗兰C目标回波信号超分辨TOA估计性能分析：

仿真内容：

仿真实验1，以直达波的到达时间作为回波到达时间的零基准点。设两个目标回波到达时间分别为50μs和80μs，信噪比均为-36dB，噪声为高斯白噪声。经过13.56s(I值)脉冲组间和脉冲组内的相干积累后先不对此信号做全变差滤波处理，然后分别使用匹配滤波方法和本发明的方法估计罗兰C目标回波的到达时间。在使用本发明的方法对罗兰C目标回波的到达时间进行估计时，首先按照本发明的步骤对时域自适应滤波后信号进行1024点(M＝1024)的离散傅里叶变换，在重构罗兰C目标回波时域稀疏信号x的过程中，N＝1024。然后，对罗兰C目标回波的到达时间进行估计(共两种方法，一种方法为直接对罗兰C目标回波的到达时间进行估计，另一种方法为在对罗兰C目标回波进行恒虚警检测之后，对罗兰C目标回波的到达时间进行估计)，当需要对罗兰C目标回波进行恒虚警检测之时，将虚警率设为0.0281。参照图3a，为本发明仿真实验1得出的重构结果示意图。在图3a中，横轴表示相对于直达波的到达时间的时延，单位为微秒，纵轴表示罗兰C目标回波信号的幅度。图3a中，CA-CFAR检测阈值表示恒虚警检测阈值，压缩感知表示本发明的方法，匹配滤波表示现有的匹配滤波方法。参照图3b，为图3a的局部放大示意图。在图3a和图3b中，两个最高的尖峰分别表示2个罗兰C目标回波信号。

仿真实验2，其仿真条件与仿真实验1相同，采用全变差滤波后，对罗兰C目标回波的到达时间进行估计。在采用全变差滤波后，仿真实验2的内容与仿真实验1相同。参照图4a，为本发明仿真实验2得出的重构结果示意图。在图4a中，横轴表示相对于直达波的到达时间的时延，单位为微秒，纵轴表示罗兰C目标回波信号的幅度。图4a中，CA-CFAR检测阈值表示恒虚警检测阈值，压缩感知表示本发明的方法，匹配滤波表示现有的匹配滤波方法。参照图4b，为图4a的局部放大示意图。在图4a和图4b中，两个最高的尖峰分别表示2个罗兰C目标回波信号。

仿真实验3，仍然设定两个目标回波的到达时间为50μs和80μs，分析随信噪比变化对罗兰C目标回波的到达时间的估计性能的影响。设两个目标回波信噪比相同，令此信噪比从-46dB到-24dB依次变化，每种情况非别做150次独立蒙特卡罗实验，重构后能正确估计罗兰C目标回波的到达时间的成功概率与信噪比的关系如图5所示。参照图5，为本发明仿真实验3得出的正确估计罗兰C目标回波的到达时间的成功概率与信噪比的关系曲线图。在图5中，横轴表示罗兰C目标回波信噪比，单位为dB，纵轴表示正确估计罗兰C目标回波的到达时间的成功概率，单位为％。

仿真实验4，设定目标1的到达时间为50μs，分析两个目标回波的到达时间相对间隔的变化对目标回波TOA估计性能的影响。设定两个目标回波信噪比均为-36dB，令目标2回波的到达时间相对于目标1回波的到达时间从1μs到25μs变化，每种情况非别做150次独立蒙特卡罗实验，重构后能正确估计两个目标回波到达时间的成功概率与目标2相对目标1的回波到达时间差的变化关系如图6所示。参照图6，为本发明仿真实验4得出的正确估计两个目标回波到达时间的成功概率与两个目标的回波到达时间差的变化关系示意图。在图6中，横轴表示目标2相对目标1的回波到达时间差，单位为微秒，纵轴表示确估计两个目标回波到达时间的成功概率，单位为％。

仿真结果分析：

仿真实验1，从图3a和图3b中可以看出，使用匹配滤波方法(在图3a和图3b中，标记为匹配滤波)时，无法分辨2个目标回波的到达时间，进行全变差滤波(在图3a和图3b中，标记为压缩感知)后，能够有效分辨2个目标回波的到达时间，此时，两个最大尖峰的位置分别为50μs和80μs，与两个目标回波的到达时间一致。最大尖峰两侧周围存在一些副瓣，这些副瓣可以通过适当的选取合适的虚警率通过恒虚警检测剔除。另外，重构结果中存在一些模值较小的尖峰，这些尖峰采用恒虚警检测无法剔除，会造成虚假目标到达时间A估计。

仿真实验2，从图4a和图4b中可以看出，使用匹配滤波方法(在图4a和图4b中，标记为匹配滤波)时，无法分辨2个目标回波的到达时间，进行全变差滤波(在图4a和图4b中，标记为压缩感知)后，能够有效分辨2个目标回波的到达时间，进行全变差滤波后，图3a和图3b中出现的虚假目标TOA估计对应的模值较小的尖峰全部消失，全变差滤波处理提升了目标回波的到达时间的估计的正确率，此外通过对全变差滤波后的结果进行恒虚警检测，可以得到目标回波的到达时间的精确估计。

仿真实验3，从图5可以看出，正确估计罗兰C目标回波的到达时间的成功概率随着罗兰C目标回波信噪比的增加而增加。当罗兰C目标回波的信噪比大于-36dB时，采用本发明的方法正确估计罗兰C目标回波的到达时间的成功概率为100％。

仿真实验4，从图6可以看出，正确估计两个目标回波到达时间的成功概率随着两个目标的回波到达时间差增加而增加，当两个目标的回波到达时间差大于22μs时，正确估计两个目标回波到达时间的成功概率大于98％。

综合以上分析过程可以得出如下结论：针对以罗兰C信号作为照射源的被动雷达中目标回波时宽带宽积较小，匹配滤波结果分辨率较差的问题，本发明先利用自适应滤波消除接收信号中的直达波和天波，然后对其做全变差滤波处理，最后根据频域样本采用压缩感知重构携带目标TOA信息的时域稀疏信号，并结合恒虚警检测对目标回波的个数进行检测并对其进行超分辨TOA估计。仿真结果表明了本发明的有效性。

显然，本领域的技术人员可以对本发明进行各种改动和变型而不脱离本发明的精神和范围。这样，倘若本发明的这些修改和变型属于本发明权利要求及其等同技术的范围之内，则本发明也意图包含这些改动和变型在内。

Claims (6)

1.基于全变差和压缩感知的罗兰C被动雷达TOA估计方法，其特征在于，包括以下步骤：

S1：利用被动雷达接收罗兰C直达波信号、罗兰C天波信号以及罗兰C目标回波信号，然后对被动雷达接收的罗兰C直达波信号、罗兰C天波信号以及罗兰C目标回波信号进行时域自适应滤波，得出时域自适应滤波后信号x_sur(n)，n为离散时间变量；

S2：对时域自适应滤波后信号x_sur(n)进行全变差滤波，得出全变差滤波后信号

S3：对全变差滤波后信号进行离散傅里叶变换，得出全变差滤波后频域信号k为离散频点变量；根据全变差滤波后频域信号采用压缩感知方法重构出罗兰C目标回波时域稀疏信号x；

S4：根据所述罗兰C目标回波时域稀疏信号x，得出罗兰C目标回波的到达时间。

2.如权利要求1所述的基于全变差和压缩感知的罗兰C被动雷达TOA估计方法，其特征在于，在步骤S1中，时域自适应滤波后信号x_sur(n)表示为：

$x_{sur}\left(n\right)=s\left(n\right)+e_{sur}\left(n\right)=\underset{i=1}{\overset{L}{\&Sigma;}}{\mathrm{\&alpha;}}_i{x}_g(n-{\mathrm{\&tau;}}_i)+e_{sur}\left(n\right)$

其中，s(n)表示不含噪声的罗兰C目标回波信号，n为离散时间变量，n为整数且n取0至N-1，N表示时间长度，e_sur(n)表示被动雷达的主天线接收到的噪声信号分量；x_g(n-τ_i)表示第i个罗兰C目标回波信号，α_i表示第i个罗兰C目标回波信号的幅度，τ_i表示待估计的第i个罗兰C目标回波信号的到达时间，i取1至L，L表示罗兰C目标回波信号的个数；

在步骤S2中，对时域自适应滤波后信号x_sur(n)进行全变差滤波的过程为：求解以下优化问题：

$\hat{s}\left(n\right)=\underset{s\left(n\right)}{\mathrm{argmin}}\underset{n=0}{\overset{N-1}{\&Sigma;}}|x_{sur}\left(n\right)-s\left(n\right){|}^2+\mathrm{\&mu;}\underset{n=0}{\overset{N-1}{\&Sigma;}}|s\left(n\right)-s(n-1)|$

其中，表示待求解的全变差滤波后信号，|·|表示取模，μ为设定的大于0的常数；

通过求解以上优化问题得出全变差滤波后信号

3.如权利要求1所述的基于全变差和压缩感知的罗兰C被动雷达TOA估计方法，其特征在于，在步骤S2中，全变差滤波后信号表示为：

$\hat{s}\left(n\right)=s\left(n\right)+{\hat{e}}_{sur}\left(n\right)=\underset{i=1}{\overset{L}{\&Sigma;}}{\mathrm{\&alpha;}}_i{x}_g(n-{\mathrm{\&tau;}}_i)+{\hat{e}}_{sur}\left(n\right)$

其中，s(n)表示不含噪声的罗兰C目标回波信号，n为整数且n取0至N-1，N表示时间长度，x_g(n-τ_i)表示第i个罗兰C目标回波信号，α_i表示第i个罗兰C目标回波信号的幅度，τ_i表示第i个罗兰C目标回波信号的到达时间，i取1至L，L表示罗兰C目标回波信号的个数；表示经全变差滤波后被动雷达的主天线接收到的噪声信号分量；

在步骤S3中，对全变差滤波后信号进行离散傅里叶变换，得出全变差滤波后频域信号k为离散频点变量；全变差滤波后频域信号表示为：

$\hat{S}\left(k\right)=X_g\left(k\right)\underset{i=1}{\overset{L}{\&Sigma;}}{\mathrm{\&alpha;}}_i{e}^{-j\frac{2\mathrm{\&pi;}}{M}{\mathrm{k\&tau;}}_i}+{\hat{E}}_{sur}\left(k\right)$

其中，为经离散傅里叶变换后得出的频域信号，k为整数且k取1至M，M为离散傅里叶变换的点数；X_g(k)为x_g(n)经离散傅里叶变换后得出的频域信号，x_g(n)为：

x_g(n)＝ax₀(n)*G(n)

其中，

x₀(n)＝An²exp[-2n/65]sin(0.2πn)

$G\left(n\right)=\underset{m=0}{\overset{7}{\&Sigma;}}\&lsqb;\mathrm{\&delta;}(n-1000m)pca\left(m\right)+\mathrm{\&delta;}(n-I-1000m)pcb\left(m\right)\&rsqb;$

其中，a为常数，x₀(n)表示罗兰C发射台的发射信号，A为设定的常数，*表示卷积；m取0至7，I表示以微秒为单位的罗兰C脉冲组信号的脉冲组重复周期，δ(·)表示罗兰C脉冲组信号的单位脉冲函数，pca(m)表示罗兰C脉冲组信号的奇数脉冲组重复周期内子脉冲的二相码，pcb(m)表示罗兰C脉冲组信号的偶数脉冲组重复周期内子脉冲的二相码；

x₀(n)经离散傅里叶变换后得出的频域信号表示为X₀(k)，将全变差滤波后频域信号除以X₀(k)，得出：

$F\left(k\right)=\hat{S}\left(k\right)/X_0\left(k\right)=\underset{i=1}{\overset{L}{\&Sigma;}}{\mathrm{\&beta;}}_i{e}^{-j\frac{2\mathrm{\&pi;}}{M}{\mathrm{k\&tau;}}_i}+{\hat{E}}_{sur}\left(k\right)/X_0\left(k\right)$

其中，β_i＝aα_i；k为整数且k取1至M，M为离散傅里叶变换的点数；

然后，设置矩形窗函数W(k)，当90kHz≤k≤110kHz时，W(k)＝1；当k<90kHz或k>110kHz时，W(k)＝0；将F(k)加上矩形窗函数W(k)，得出：

$Y\left(k\right)=F\left(k\right)W\left(k\right)=\underset{i=1}{\overset{L}{\&Sigma;}}{\mathrm{\&beta;}}_i{e}^{-j\frac{2\mathrm{\&pi;}}{M}{\mathrm{k\&tau;}}_i}+P\left(k\right)$

其中， $P\left(k\right)=W\left(k\right){\hat{E}}_{sur}\left(k\right)/X_0\left(k\right);$

将Y(1)至Y(M)组成矢量y，y＝[Y(1),...,Y(k),...,Y(M)]^T，其中，T表示矩阵或向量的转置；则矢量y表示为：

y＝Φx+p

其中，p＝[P(1),...,P(k),...,P(M)]^T，x罗兰C目标回波时域稀疏信号，x为N×1维的稀疏矢量，Φ为M×N维的测量矩阵，Φ的矩阵展开形式为：

$\mathrm{\&Phi;}=\left[\begin{array}{cccccc}{e}^{-j\frac{2\mathrm{\&pi;}}{M}{\mathrm{\&tau;}}_1}& e^{-j\frac{2\mathrm{\&pi;}}{M}{\mathrm{\&tau;}}_2}& ...& e^{-j\frac{2\mathrm{\&pi;}}{M}{\mathrm{\&tau;}}_L}& ...& e^{-j\frac{2\mathrm{\&pi;}}{M}{\mathrm{\&tau;}}_N}\\ e^{-j\frac{2\mathrm{\&pi;}}{M}2{\mathrm{\&tau;}}_1}& e^{-j\frac{2\mathrm{\&pi;}}{M}2{\mathrm{\&tau;}}_2}& ...& e^{-j\frac{2\mathrm{\&pi;}}{M}2{\mathrm{\&tau;}}_L}& ...& e^{-j\frac{2\mathrm{\&pi;}}{M}2{\mathrm{\&tau;}}_N}\\ \&CenterDot;& \&CenterDot;& & \&CenterDot;& & \&CenterDot;\\ \&CenterDot;& \&CenterDot;& ...& \&CenterDot;& ...& \&CenterDot;\\ \&CenterDot;& \&CenterDot;& & \&CenterDot;& & \&CenterDot;\\ e^{-j\frac{2\mathrm{\&pi;}}{M}{\mathrm{M\&tau;}}_1}& e^{-j\frac{2\mathrm{\&pi;}}{M}{\mathrm{M\&tau;}}_2}& ...& e^{-j\frac{2\mathrm{\&pi;}}{M}{\mathrm{M\&tau;}}_L}& ...& e^{-j\frac{2\mathrm{\&pi;}}{M}{\mathrm{M\&tau;}}_N}\end{array}\right]$

其中，L<N；然后，求解以下l₀范数最优化问题：

min||x||₀

s.t.y＝Φx+p

其中，||x||₀表示x的l₀范数；通过求解以上l₀范数最优化问题重构出罗兰C目标回波时域稀疏信号x。

4.如权利要求3所述的基于全变差和压缩感知的罗兰C被动雷达TOA估计方法，其特征在于，在步骤S4中，将所述l₀范数最优化问题替换为以下l₁范数最优化问题：

min||x||₁s.t.y＝Φx+p

其中，||·||₁表示求l₁范数。

5.如权利要求1所述的基于全变差和压缩感知的罗兰C被动雷达TOA估计方法，其特征在于，在步骤S4中，在所述罗兰C目标回波时域稀疏信号x中，找出各个非零元素在x中的位置；根据每个非零元素在x中的位置、以及被动雷达的采样频率，得出对应的罗兰C目标回波的到达时间。

6.如权利要求1所述的基于全变差和压缩感知的罗兰C被动雷达TOA估计方法，其特征在于，在步骤S4中，根据所述罗兰C目标回波时域稀疏信号x，对罗兰C目标回波进行恒虚警检测，得出经恒虚警检测后稀疏矢量，然后根据所述恒虚警检测后稀疏矢量，找出各个非零元素在恒虚警检测后稀疏矢量中的位置；根据每个非零元素在恒虚警检测后稀疏矢量中的位置、以及被动雷达的采样频率，得出对应的罗兰C目标回波的到达时间。