CN102353947A - 一种基于csa-mwf的无源雷达目标回波信号子空间的估计方法 - Google Patents

Abstract

本发明提供一种基于CSA-MWF的无源雷达目标回波信号子空间的估计方法，包括：步骤一：在无源雷达接收系统中提取观测数据矢量，并将其赋值给CSA-MWF的初始观测数据，并初始化期望信号；步骤二：推导目标回波子空间的估计方法表达式；步骤三：计算CSA-MWF中本级的前向滤波器；步骤四：计算CSA-MWF中本级的期望信号；步骤五：计算CSA-MWF中更新后的观测数据；步骤六：进行门限判决；步骤七：计算得出目标回波信号子空间。本发明将CSA-MWF这种有效的降维方法应用于无源雷达中，可以避免估计观测数据的协方差矩阵，避免对其进行特征值分解，可以有效降低计算量，非常适合信号多变的复杂环境。

Description

一种基于CSA-MWF的无源雷达目标回波信号子空间的估计方法

技术领域

本发明属于无源雷达领域，具体涉及一种基于CSA-MWF的无源雷达目标回波信号子空间的估计方法。

背景技术

无源雷达是指雷达本身不发射电磁波信号而只用目标辐射电磁波信号(外辐射源)进行目标探测和跟踪的雷达，它具有良好的“四抗性能”，并具有造价低、隐蔽性强、机动性高等优点。目标辐射的电磁信号可能是目标自身发射的信号，或者是第三方电磁波信号经目标反射后的电磁信号。因此，根据目标辐射信号源的类型，无源雷达可分为两类：一是利用目标自身辐射源的无源雷达，包括待观测的目标自身携带的辐射源，如雷达、通信、应答机、有源干扰和导航等电子设备；二是利用第三方发射信号经目标反射的信号的无源雷达，这类发射源包括地面广播电台、电视台、通信台站、直播电视卫星和卫星导航定位系统等。从雷达体制上讲，它的接收机和发出外辐射源的发射机是异地配置，因此它也属于双(多)基雷达。

无源雷达具有很多传统有源雷达没有的优点：(1)无源雷达接收机由于没有大功率的器件，且不受发射机功率泄漏的影响，从而具有较高的灵敏度；(2)多个接收机工作，可以对干扰源进行无源定位；(3)由于收、发分置，接收机是静默的，也可以机动，所以系统具有抗电子侦查、抗干扰和抗反辐射导弹摧毁的能力；(4)在接收机前置时，可利用空中、空间的照射源，可以探测到发射机的视线以下和远区超低空的目标，抗超低空突防能力强；(5)具有反隐身的效果，由于隐身目标只在鼻锥正负三十度范围内有极小的雷达截面积(RCS)，而侧向及顶部散射和绕射并没有减小，还很强，使得无源雷达对隐身目标的探测具有很好的效果；(6)由于自身不发射电磁波，无源雷达系统具有隐蔽性和突然性，系统生存能力强。

子空间类算法是一类重要的方法，该算法将观测空间分解为信号子空间和噪声子空间。传统的子空间分解方法是估计阵列协方差矩阵并对其进行特征分解(EVD)，大特征值对应的特征向量张成信号子空间，小特征值对应的矢量张成噪声子空间，需要的总运量为O(M²N)+O(M³)，其中M、N分别为阵列个数和采样快拍数。若阵元数M较多，则计算量很大，不利于实时处理。无源雷达的核心技术是无源相干定位技术，其基本的思想是以外辐射源发射的直达波信号作为参考，检测分析目标反射辐射源发射的信号能量，估计出目标反射信号的到达方向、到达时间和多普勒频移等参数，从而实现对目标的定位和跟踪。无源雷达系统在进行信号处理的过程中需要进行长时间的相干累积，计算量很大，尤其对其进行空时二维联合处理时，计算量更大。因此要考虑在空时二维联合处理时，对其进行降秩运算以减少计算量。为了降低计算复杂度，众多学者提出了一系列降秩算法，如主分量法、交叉谱法和Lanczos迭代算法等。但是这些算法都需要对观测数据的协方差进行分解，使得计算量仍旧很大。Goldstein等人提出的多级维纳滤波算法(GRS-MWF)是一种新的降秩处理方法，但是其前向分解滤波器相互并不正交，数值稳健性不好，用它进行无源雷达目标子空间估计时效果不好，无法达到无源雷达性能要求。

发明内容

针对现有技术中存在的不足，本发明提供一种基于CSA-MWF的无源雷达目标回波信号子空间的估计方法，在无源雷达接收系统中，在避免计算量过大的同时对接收信号进行有效地快速地目标回波子空间分解，从而检测到目标回波的方法。

本发明提出的一种基于CSA-MWF的无源雷达目标回波信号子空间的估计方法，具体包括以下几个步骤：

步骤一：在无源雷达接收系统中提取观测数据矢量，并将其赋值给相关相减结构的多级维纳滤波器(CSA-MWF)的初始观测数据，初始化期望信号d₀。

无源雷达的接收阵元是阵元数为M的等距线阵(LUA)，则该接收阵元在k时刻接收的M维观测数据矢量x(k)为：

x(k)＝[a(θ₁)，a(θ₂)，…，a(θ_P)]s(k)+n(k)    (1)

＝A(θ)s(k)+n(k)

其中，s(k)表示目标回波信号复振幅矢量，为P×1阶矩阵，n(k)表示空时白噪声复矢量，n(k)为M×1阶矩阵，A(θ)表示目标回波的方向矩阵，为M×P阶矩阵，k为采样时刻k＝0，1，…，N-，N是快拍数，P是目标回波信号的个数，θ₁，…，θ_P分别为目标回波信号1至P的入射角，a(θ₁)，a(θ₂)，…，a(θ_P)分别为目标回波信号1至P的导向矢量。

设定第一个接收阵元为基准阵元，则任一目标i的导向矢量a(θ_i)具有如下的结构：

$a\left({\mathrm{\θ}}_i\right)=\frac{1}{\sqrt{M}}{[1,e^{j\frac{2\mathrm{\πd}}{\mathrm{\λ}}\sin \left({\mathrm{\θ}}_i\right)},\·\·\·,e^{j(M-1)\frac{2\mathrm{\πd}}{\mathrm{\λ}}\sin \left({\mathrm{\θ}}_i\right)}]}^T---\left(2\right)$

其中，目标i的入射角

d表示阵元间距，λ表示载波波长，T表示矩阵转置，M为阵元数，

表示每个阵元相对基准阵元的相移。

设定雷达接收阵元数大于目标回波信号个数，即M＞P。加性噪声为独立同分布的满足(0，σ²)的空时高斯白噪声矢量，即：

E[n(k)ⁿH(l)]＝σ²I_M    (3)

E[n(k)n^T(l)]＝0        (4)

其中，(·)^H表示共轭转置，E[·]表示求取数学期望。σ²表示空时高斯白噪声的方差，n(k)表示空时高斯白噪声，n^H(l)表示空时高斯白噪声的共轭转置，n^T(l)表示空时高斯白噪声的转置，I_M表示M维的单位矩阵。

步骤二：推导目标回波子空间的估计方法的表达式。

设定目标回波信号和加性噪声是不相关的，则观测数据的协方差矩阵R_x为：

R_x＝E[x(k)x^H(k)]＝A(θ)R_sA^H(θ)+σ²I_M    (5)

其中，R_s为目标回波信号的协方差矩阵；A(θ)表示目标回波信号的方向矩阵，A^H(θ)表示A(θ)的共轭转置，x(k)表示观测数据矢量，x^H(k)表示观测数据矢量的共轭转置，σ²表示空时高斯白噪声的方差，E[·]表示求取数学期望，I_M表示M维的单位矩阵。

对观测数据的协方差矩阵做特征值分解：

$R_x=\underset{i=1}{\overset{M}{\mathrm{\Σ}}}{\mathrm{\λ}}_i{v}_i{v}_i^H=V_s{\mathrm{\Λ}}_s{V}_s^H+{\mathrm{\σ}}^2{V}_n{V}_n^H---\left(6\right)$

其中，特征值λ₁＞λ₂＞…＞λ_P＞λ_P+1＝…＝λ_M＝σ²，其中较大的特征值λ₁，…，λ_p对应目标回波信号，较小的特征值λ_P+1，…，λ_M对应噪声信号，即V_s＝[v₁，v₂，…，v_P]，V_n＝[v_P+1，v_P+2，…，v_M]。M为阵元个数，同时也是特征值的个数。V_s的列数等于目标回波信号协方差矩阵R_s的秩P，从而张成A(θ)的P维子空间。v_i表示观测数据协方差矩阵的特征向量，

表示观测数据协方差矩阵的特征向量的共轭转置，V_s表示目标回波信号子空间的特征向量组成的矩阵，

表示V_s的共轭转置，Λ_s表示目标回波信号子空间的特征值组成的对角矩阵，V_n表示噪声子空间的特征向量组成的矩阵，

表示V_n的共轭转置。

由公式(5)和公式(6)得到：

V_s＝A(θ)Q    (7)

其中，Q＝R_sA^H(θ)V_s(Λ_s-σ²I_M)^-1，Q为P维满秩矩阵。

CSA-MWF相当于在最小均方意义下得到Wiener-Hopf方程R_xw＝r_xd的渐进最优解(其中w表示滤波器权值，r_xd表示观测数据x与期望信号d的互相关函数)，CSA-MWF的各级前向滤波器h_i，i＝{1，2，…，P}是相互正交的(i表示第i级前向滤波器，P表示前向滤波器的级数)，所以级数为P的CSA-MWF相当于Wiener-Hopf方程在Krylov子空间 ${\mathrm{\κ}}^{\left(P\right)}(R_{x_0},r_{x_0{d}_0})=\mathrm{span}\{r_{x_0{d}_0},R_{x_0}{r}_{x_0{d}_0},\·\·\·,R_{x_0}^{(P-1)}{r}_{x_0{d}_0}\}$ 的解，其中x₀表示输入的观测向量，d₀表示期望向量，

表示x₀与d₀的互相关函数，表示x₀的自相关函数，

表示的(P-1)次方，κ^(P)表示P级的Krylov子空间，span{·}表示将括号中的向量张成空间。目标回波信号子空间按照如下进行估计：

$\mathrm{span}\{h_1,h_2,\·\·\·,h_P\}=\mathrm{span}\{r_{x_0{d}_0},R_{x_0}{r}_{x_0{d}_0},\·\·\·,R_{x_0}^{(P-1)}{r}_{x_0{d}_0}\}---\left(8\right)$

因此，存在一个P阶满秩矩阵K，使得公式(9)成立：

$[h_1,h_2,\·\·\·,h_P]=[r_{x_0{d}_0},R_{x_0}{r}_{x_0{d}_0},\·\·\·,R_{x_0}^{(P-1)}{r}_{x_0{d}_0}]K---\left(9\right)$

令T_s＝[h₁，h₂，…，h_P]，T_n＝[h_P+1，h_P+2，…，h_M]。由于

和其中I_P和I_M-P分别表示维数为P和(M-P)的单位阵，由公式(6)可得下式成立，即

${R_x}^{\left(i\right)}=V_s{\mathrm{\Λ}}_s^{\left(i\right)}{V}_s^H+{\mathrm{\σ}}^{2i}{V}_n{V}_n^H---\left(10\right)$

其中，R_x ⁽ⁱ⁾表示R_x的i次方。i＝1，2，…，P-1。由CSA-MWF性质可知，

落在目标回波信号子空间中，从而有把公式(10)带入公式(9)中，同时

和V_s＝A(θ)Q，即可得到公式(11)，其中K表示一个P阶满秩矩阵，

表示目标子空间的特征值组成的对角矩阵Λ_s的(P-1)次方，Γ表示表示一个P阶矩阵具体如公式(12)所示，H表示一个P阶矩阵具体如公式(13)所示，Q＝R_sA^H(θ)V_s(Λ_s-σ²I)^-1。

$T_s=[V_s{V}_s^H{r}_{x_0{d}_0},V_s{\mathrm{\Λ}}_s{V}_s^H{r}_{x_0{d}_0},\·\·\·,V_s{\mathrm{\Λ}}_s^{(P-1)}{V}_s^H{r}_{x_0{d}_0}]K$

$=V_s[V_s^H{r}_{x_0{d}_0},{\mathrm{\Λ}}_s{V}_s^H{r}_{x_0{d}_0},\·\·\·,{\mathrm{\Λ}}_s^{(P-1)}{V}_s^H{r}_{x_0{d}_0}]K$

$=A\left(\mathrm{\θ}\right)Q[V_s^H{r}_{x_0{d}_0},{\mathrm{\Λ}}_s{V}_s^H{r}_{x_0{d}_0},\·\·\·,{\mathrm{\Λ}}_s^{(P-1)}{V}_s^H{r}_{x_0{d}_0}]K---\left(11\right)$

$=A\left(\mathrm{\θ}\right)\mathrm{Q\ΓK}$

$=A\left(\mathrm{\θ}\right)H$

其中：

$\mathrm{\Γ}=[V_s^H{r}_{x_0{d}_0},{\mathrm{\Λ}}_s{V}_s^H{r}_{x_0{d}_0},\·\·\·,{\mathrm{\Λ}}_s^{(P-1)}{V}_s^H{r}_{x_0{d}_0}]---\left(12\right)$

H＝QΓK    (13)

Γ是满秩矩阵，又因为Q和K均是非奇异矩阵，因而H也是非奇异矩阵，所以由式(11)可得P阶的目标回波信号子空间的表达式：

$\mathrm{span}\{h_1,h_2,\·\·\·,h_P\}=\mathrm{col}\left\{A\left(\mathrm{\θ}\right)\right\}\stackrel{\mathrm{\Δ}}{=}{\mathrm{\Φ}}_s^{\left(P\right)}---\left(14\right)$

其中，col{}表示列空间，即为其列矢量所有线性组合集合成的空间，

表示P阶的目标回波信号子空间。

由CSA-MWF的M个相互正交的匹配滤波器构成的预滤波矩阵为T_M＝[h₁，h₂，…h_P，h_P+1，…，h_M]，由于其所有列矢量h_j，j＝1，2，…，M均相互正交，则h_k⊥col{A(θ)}，k＝P+1，P+2，…，M，所以h_k位于由A(θ)的各列向量张成的列空间col{A(θ)}的正交补子空间，即(M-P)维的噪声子空间的表达式：

$\mathrm{span}\{h_{P+1},h_{P+2},\·\·\·,h_M\}=\mathrm{null}\left\{A\left(\mathrm{\θ}\right)\right\}\stackrel{\mathrm{\Δ}}{=}{\mathrm{\Φ}}_n^{(M-P)}---\left(15\right)$

其中null{·}表示括号中空间的正交补子空间；表示(M-P)维的噪声子空间。

步骤三：按照下式计算CSA-MWF中前向分解的第j级前向滤波器h_j，i＝{1，2，…，P，…，M}：

$h_j=\frac{E\left[d_{j-1}^{*}\left(k\right)x_{j-1}\left(k\right)\right]}{{\left|\right|E\left[d_{j-1}^{*}\left(k\right)x_{j-1}\left(k\right)\right]\left|\right|}_2}---\left(16\right)$

x_j-1(k)表示第(j-1)级前向滤波器的观测数据；表示第(j-1)级前向滤波器的期望信号d_j-1(k)的共轭信号。

步骤四：按照下式计算CSA-MWF中前向分解的第j级期望信号d_j(k)：

$d_j\left(k\right)=h_j^H{x}_{j-1}\left(k\right)---\left(17\right)$

步骤五：按照下式计算CSA-MWF中前向分解的各级前向滤波器更新后的观测数据：

x_j(k)＝x_j-1(k)-h_jd_j(k)    (18)

其中x_j(k)表示第j级前向滤波器的观测数据；d_j(k)表示第j级前向滤波器的期望信号。

步骤六：进行门限判决，若|x_j(k)|²≤2MNσ²即j＝P，则进行步骤七，其中M为阵元数，N为快拍数；否则令j＝j+1，返回步骤三计算下一级前向滤波器。

步骤七：将计算得出的CSA-MWF中的各级前向滤波器h₁，h₂，…，h_P带入公式(14)中，计算得出目标回波信号子空间 ${\mathrm{\Φ}}_s^{\left(P\right)}=\mathrm{span}\{h_1,h_2,...,h_P\}.$

本发明在完成步骤七之后还优选的包含由计算噪声子空间的步骤具体为：

步骤八：(1)令j＝P+1，返回执行步骤三，计算CSA-MWF的第j级前向滤波器，计算得出前向滤波器h_P+1；

(2)计算CSA-MWF中前向分解的第j级期望信号

(3)计算CSA-MWF中前向分解的第j级前向滤波器更新后的观测数据：

x_j(k)＝x_j-1(k)-h_jd_j(k)

其中x_j(k)表示第j级前向滤波器的观测数据；d_j(k)表示第j级前向滤波器的期望信号。

(4)判断j＝M是否成立，若成立，将计算得到的各级前向滤波器h_P+1，h_P+2，…，h_M代入公式(15)中，计算得到噪声子空间

否则令j＝j+1，返回步骤八(1)，计算下一级前向滤波器。

本发明的优点在于：

(1)本发明提供一种基于CSA-MWF的无源雷达目标回波信号子空间的估计方法，将CSA-MWF这种有效的降维方法应用于无源雷达中，可以避免估计观测数据的协方差矩阵，避免对其进行特征值分解，可以有效降低计算量，非常适合信号多变的复杂环境；

(2)本发明提供一种基于CSA-MWF的无源雷达目标回波信号子空间的估计方法，CSA-MWF的维数可以小于目标子空间的实际维数时，收敛速度快，收敛需要的快拍低；

(3)本发明提供一种基于CSA-MWF的无源雷达目标回波信号子空间的估计方法，可以有效准确地完成无源雷达目标子空间的快速估计，保持优良的抗干扰性能；

(4)本发明提供一种基于CSA-MWF的无源雷达目标回波信号子空间的估计方法，其中各级前向滤波器均正交，具有很好的数值稳健性。

附图说明

图1：本发明提供一种基于CSA-MWF的无源雷达目标回波信号子空间的估计方法的流程图；

图2：本发明提供一种基于CSA-MWF的无源雷达目标回波信号子空间的估计方法中基于外辐射源无源雷达系统结构示意图。

具体实施方式

下面将结合附图对本发明作进一步的详细说明。

本发明提出的一种基于CSA-MWF的无源雷达目标回波信号子空间的估计方法，如图1所示，具体包括以下几个步骤：

步骤一：在无源雷达接收系统中提取观测数据矢量，并将其赋值给相关相减结构的多级维纳滤波器(CSA-MWF)的初始观测数据，初始化期望信号d₀。

如图2所示，基于外辐射源的无源雷达系统属于多基地无源雷达探测系统，由M个接收阵元的等间距d＝λ/2的均匀线阵组成无源雷达接收系统。本发明采用相关相减结构的多级维纳滤波器CSA-MWF的前向分解特性，相关相减结构可有效的降低多级维纳滤波器前向分解的计算量，并且CSA-MWF为酉多级维纳滤波器，具有良好的降秩性能。

无源雷达的接收阵元是阵元数为M的等距线阵(LUA)，则该接收阵元在k时刻接收的M维观测数据矢量x(k)为：

x(k)＝[a(θ₁)，a(θ₁)，…，a(θ_P)]s(k)+n(k)    (1)

＝A(θ)s(k)+n(k)

其中，s(k)表示目标回波信号复振幅矢量，为P×1阶矩阵，n(k)表示空时白噪声复矢量，n(k)为M×1阶矩阵，A(θ)表示目标回波的方向矩阵，为M×P阶矩阵，k为采样时刻k＝0，1，…，N-，N是快拍数，P是目标回波信号的个数，θ₁，…，θ_P分别为目标回波信号1至P的入射角，a(θ₁)，a(θ₂)，…，a(θ_P)分别为目标回波信号1至P的导向矢量。

设定第一个接收阵元为基准阵元，则任一目标i的导向矢量a(θ_i)具有如下的结构：

$a\left({\mathrm{\θ}}_i\right)=\frac{1}{\sqrt{M}}{[1,e^{j\frac{2\mathrm{\πd}}{\mathrm{\λ}}\sin \left({\mathrm{\θ}}_i\right)},\·\·\·,e^{j(M-1)\frac{2\mathrm{\πd}}{\mathrm{\λ}}\sin \left({\mathrm{\θ}}_i\right)}]}^T---\left(2\right)$

其中，目标i的入射角

d表示阵元间距，λ表示载波波长，T表示矩阵转置，M为阵元数，

表示每个阵元相对基准阵元的相移。

设定雷达接收阵元数大于目标回波信号个数，即M＞P。加性噪声为独立同分布的满足(0，σ²)的空时高斯白噪声矢量，即：

E[n(k)n^H(l)]＝σ²I_M    (3)

E[n(k)n^T(l)]＝0        (4)

其中，(·)^H表示共轭转置，E[·]表示求取数学期望。σ²表示空时高斯白噪声的方差，n(k)表示空时高斯白噪声，n^H(l)表示空时高斯白噪声的共轭转置，n^T(l)表示空时高斯白噪声的转置，I_M表示M维的单位矩阵。

步骤二：推导目标回波子空间的估计方法的表达式。

设定目标回波信号和加性噪声是不相关的，则观测数据的协方差矩阵R_x为：

R_x＝E[x(k)x^H(k)]＝A(θ)R_sA^H(θ)+σ²I_M    (5)

其中，R_s为目标回波信号的协方差矩阵；A(θ)表示目标回波信号的方向矩阵，A^H(θ)表示A(θ)的共轭转置，x(k)表示观测数据矢量，x^H(k)表示观测数据矢量的共轭转置，σ²表示空时高斯白噪声的方差，E[·]表示求取数学期望，I_M表示M维的单位矩阵。

对观测数据的协方差矩阵做特征值分解：

$R_x=\underset{i=1}{\overset{M}{\mathrm{\Σ}}}{\mathrm{\λ}}_i{v}_i{v}_i^H=V_s{\mathrm{\Λ}}_s{V}_s^H+{\mathrm{\σ}}^2{V}_n{V}_n^H---\left(6\right)$

其中，特征值λ₁＞λ₂＞…＞λ_P＞λ_P+1＝…＝λ_M＝σ²，其中较大的特征值λ₁，…，λ_p对应目标回波信号，较小的特征值λ_P+1，…，λ_M对应噪声信号，即V_s＝[v₁，v₂，…，v_P]，V_n＝[v_P+1，v_P+2，…，v_M]。M为阵元个数，同时也是特征值的个数。V_s的列数等于目标回波信号协方差矩阵R_s的秩P，从而张成A(θ)的P维子空间。v_i表示观测数据协方差矩阵的特征向量，

表示观测数据协方差矩阵的特征向量的共轭转置，V_s表示目标回波信号子空间的特征向量组成的矩阵，表示V_s的共轭转置，Λ_s表示目标回波信号子空间的特征值组成的对角矩阵，V_n表示噪声子空间的特征向量组成的矩阵，

表示V_n的共轭转置。

由公式(5)和公式(6)得到：

V_s＝A(θ)Q    (7)

其中，Q＝R_sA^H(θ)V_s(Λ_s-σ²I_M)^-1，Q为P维满秩矩阵。

CSA-MWF相当于在最小均方意义下得到Wiener-Hopf方程R_xw＝r_xd的渐进最优解(其中w表示滤波器权值，r_xd表示观测数据x与期望信号d的互相关函数)，CSA-MWF的各级前向滤波器h_i，i＝{1，2，…，P}是相互正交的(i表示第i级前向滤波器，P表示前向滤波器的级数)，所以级数为P的CSA-MWF相当于Wiener-Hopf方程在Krylov子空间 ${\mathrm{\κ}}^{\left(P\right)}(R_{x_0},r_{x_0{d}_0})=\mathrm{span}\{r_{x_0{d}_0},R_{x_0}{r}_{x_0{d}_0},\·\·\·,R_{x_0}^{(P-1)}{r}_{x_0{d}_0}\}$ 的解，其中x₀表示输入的观测向量，d₀表示期望向量，表示x₀与d₀的互相关函数，

表示x₀的自相关函数，

表示

的(P-1)次方，κ^(P)表示P级的Krylov子空间，span{·}表示将括号中的向量张成空间。目标回波信号子空间按照如下进行估计：

$\mathrm{span}\{h_1,h_2,\·\·\·,h_P\}=\mathrm{span}\{r_{x_0{d}_0},R_{x_0}{r}_{x_0{d}_0},\·\·\·,R_{x_0}^{(P-1)}{r}_{x_0{d}_0}\}---\left(8\right)$

因此，存在一个P阶满秩矩阵K，使得公式(9)成立：

$[h_1,h_2,\·\·\·,h_P]=[r_{x_0{d}_0},R_{x_0}{r}_{x_0{d}_0},\·\·\·,R_{x_0}^{(P-1)}{r}_{x_0{d}_0}]K---\left(9\right)$

令T_s＝[h₁，h₂，…，h_P]，T_n＝[h_P+1，h_P+2，…，h_M]。由于

和

其中I_P和I_M-P分别表示维数为P和(M-P)的单位阵，由公式(6)可得下式成立，即

${R_x}^{\left(i\right)}=V_s{\mathrm{\Λ}}_s^{\left(i\right)}{V}_s^H+{\mathrm{\σ}}^{2i}{V}_n{V}_n^H---\left(10\right)$

其中，R_x ⁽ⁱ⁾表示R_x的i次方。i＝1，2，…，P-1。由CSA-MWF性质可知，

落在目标回波信号子空间中，从而有

把公式(10)带入公式(9)中，同时

和V_s＝A(θ)Q，即可得到公式(11)，其中K表示一个P阶满秩矩阵，表示目标子空间的特征值组成的对角矩阵Λ_s的(P-1)次方，Γ表示表示一个P阶矩阵具体如公式(12)所示，H表示一个P阶矩阵具体如公式(13)所示，Q＝R_sA^H(θ)V_s(Λ_s-σ²I)^-1。

$T_s=[V_s{V}_s^H{r}_{x_0{d}_0},V_s{\mathrm{\Λ}}_s{V}_s^H{r}_{x_0{d}_0},\·\·\·,V_s{\mathrm{\Λ}}_s^{(P-1)}{V}_s^H{r}_{x_0{d}_0}]K$

$=V_s[V_s^H{r}_{x_0{d}_0},{\mathrm{\Λ}}_s{V}_s^H{r}_{x_0{d}_0},\·\·\·,{\mathrm{\Λ}}_s^{(P-1)}{V}_s^H{r}_{x_0{d}_0}]K$

$=A\left(\mathrm{\θ}\right)Q[V_s^H{r}_{x_0{d}_0},{\mathrm{\Λ}}_s{V}_s^H{r}_{x_0{d}_0},\·\·\·,{\mathrm{\Λ}}_s^{(P-1)}{V}_s^H{r}_{x_0{d}_0}]K---\left(11\right)$

$=A\left(\mathrm{\θ}\right)\mathrm{Q\ΓK}$

$=A\left(\mathrm{\θ}\right)H$

其中：

$\mathrm{\Γ}=[V_s^H{r}_{x_0{d}_0},{\mathrm{\Λ}}_s{V}_s^H{r}_{x_0{d}_0},\·\·\·,{\mathrm{\Λ}}_s^{(P-1)}{V}_s^H{r}_{x_0{d}_0}]---\left(12\right)$

H＝QΓK    (13)

Γ是满秩矩阵，又因为Q和K均是非奇异矩阵，因而H也是非奇异矩阵，所以由式(11)可得P阶的目标回波信号子空间的表达式：

$\mathrm{span}\{h_1,h_2,\·\·\·,h_P\}=\mathrm{col}\left\{A\left(\mathrm{\θ}\right)\right\}\stackrel{\mathrm{\Δ}}{=}{\mathrm{\Φ}}_s^{\left(P\right)}---\left(14\right)$

其中，col{}表示列空间，即为其列矢量所有线性组合集合成的空间，表示P阶的目标回波信号子空间。

由CSA-MWF的M个相互正交的匹配滤波器构成的预滤波矩阵为T_M＝[h₁，h₂，…h_P，h_P+1，…，h_M]，由于其所有列矢量h_j，j＝1，2，…，M均相互正交，则h_k⊥col{A(θ)}，k＝P+1，P+2，…，M，所以h_k位于由A(θ)的各列向量张成的列空间col{A(θ)}的正交补子空间，即(M-P)维的噪声子空间的表达式：

$\mathrm{span}\{h_{P+1},h_{P+2},\·\·\·,h_M\}=\mathrm{null}\left\{A\left(\mathrm{\θ}\right)\right\}\stackrel{\mathrm{\Δ}}{=}{\mathrm{\Φ}}_n^{(M-P)}---\left(15\right)$

其中null{·}表示括号中空间的正交补子空间；

表示(M-P)维的噪声子空间。

步骤三：按照下式计算CSA-MWF中前向分解的第j级前向滤波器h_j，i＝{1，2，…，P，…，M}：

$h_j=\frac{E\left[d_{j-1}^{*}\left(k\right)x_{j-1}\left(k\right)\right]}{{\left|\right|E\left[d_{j-1}^{*}\left(k\right)x_{j-1}\left(k\right)\right]\left|\right|}_2}---\left(16\right)$

x_j-1(k)表示第(j-1)级前向滤波器的观测数据；

表示第(j-1)级前向滤波器的期望信号d_j-1(k)的共轭信号。

步骤四：按照下式计算CSA-MWF中前向分解的第j级期望信号d_j(k)：

$d_j\left(k\right)=h_j^H{x}_{j-1}\left(k\right)---\left(17\right)$

步骤五：按照下式计算CSA-MWF中前向分解的各级前向滤波器更新后的观测数据：

x_j(k)＝x_j-1(k)-h_jd_j(k)    (18)

其中x_j(k)表示第j级前向滤波器的观测数据；d_j(k)表示第j级前向滤波器的期望信号。

步骤六：进行门限判决，若|x_j(k)|²≤2MNσ²即j＝P，则进行步骤七，其中M为阵元数，N为快拍数；否则令j＝j+1，返回步骤三计算下一级前向滤波器。

步骤七：将计算得出的CSA-MWF中的各级前向滤波器h₁，h₂，…，h_P带入公式(14)中，计算得出目标回波信号子空间 ${\mathrm{\Φ}}_s^{\left(P\right)}=\mathrm{span}\{h_1,h_2,...,h_P\}.$

本发明在完成步骤七之后还优选的包含由计算噪声子空间的步骤具体为：

步骤八：(1)令j＝P+1，返回执行步骤三，计算CSA-MWF的第j级前向滤波器，计算得出前向滤波器h_P+1；

(2)计算CSA-MWF中前向分解的第j级期望信号

(3)计算CSA-MWF中前向分解的第j级前向滤波器更新后的观测数据：

x_j(k)＝x_j-1(k)-h_jd_j(k)

其中x_j(k)表示第j级前向滤波器的观测数据；d_j(k)表示第j级前向滤波器的期望信号。

(4)判断j＝M是否成立，若成立，将计算得到的各级前向滤波器h_P+1，h_P+2，…，h_M代入公式(15)中，计算得到噪声子空间

否则令j＝j+1，返回步骤八(1)，计算下一级前向滤波器。

本发明提出的一种基于CSA-MWF的无源雷达目标回波信号子空间的估计方法，在给定某一期望信号的训练信号的条件下，无源雷达目标回波子空间和噪声子空间可以分别估计如公式(14)和公式(15)所示，因此只需要得出CSA-MWF的前向分解前向滤波器就可以估计出无源雷达的目标回波信号子空间。采用该方法应用于无源雷达空时联合处理过程中，可以有效达到降秩减少计算量的作用。

Claims (2)

1.一种基于CSA-MWF的无源雷达目标回波信号子空间的估计方法，其特征在于：具体包括以下几个步骤：

步骤一：在无源雷达接收系统中提取观测数据矢量，并将其赋值给相关相减结构的多级维纳滤波器CSA-MWF的初始观测数据，初始化期望信号d₀：

无源雷达的接收阵元是阵元数为M的等距线阵，则该接收阵元在k时刻接收的M维观测数据矢量x(k)为：

x(k)＝[a(θ₁)，a(θ₂)，…，a(θ_P)]s(k)+n(k)    (1)

＝A(θ)s(k)+n(k)

其中，s(k)表示目标回波信号复振幅矢量，为P×1阶矩阵，n(k)表示空时白噪声复矢量，n(k)为M×1阶矩阵，A(θ)表示目标回波的方向矩阵，为M×P阶矩阵，k为采样时刻k＝0，1，…，N-，N为快拍数，P为目标回波信号的个数，θ₁，…，θ_P分别为目标回波信号1至P的入射角，a(θ₁)，a(θ₂)，…，a(θ_P)分别为目标回波信号1至P的导向矢量；

设定第一个接收阵元为基准阵元，则任一目标i的导向矢量a(θ_i)具有如下的结构：

$a\left({\mathrm{\θ}}_i\right)=\frac{1}{\sqrt{M}}{[1,e^{j\frac{2\mathrm{\πd}}{\mathrm{\λ}}\sin \left({\mathrm{\θ}}_i\right)},\·\·\·,e^{j(M-1)\frac{2\mathrm{\πd}}{\mathrm{\λ}}\sin \left({\mathrm{\θ}}_i\right)}]}^T---\left(2\right)$

其中，目标i的入射角

d表示阵元间距，λ表示载波波长，T表示矩阵转置，M为阵元数，

分别表示每个阵元相对基准阵元的相移；

设定雷达接收阵元数大于目标回波信号个数，M＞P；加性噪声为独立同分布的满足(0，σ²)的空时高斯白噪声矢量：

E[n(k)n^H(l)]＝σ²I_M    (3)

E[n(k)n^T(l)]＝0        (4)

其中，(·)^H表示共轭转置，E[·]表示求取数学期望，σ²表示空时高斯白噪声的方差，n(k)表示空时高斯白噪声，n^H(l)表示空时高斯白噪声的共轭转置，n^T(l)表示空时高斯白噪声的转置，I_M表示M维的单位矩阵；

步骤二：推导目标回波子空间的估计方法的表达式：

目标回波信号和加性噪声不相关，观测数据的协方差矩阵R_x为：

R_x＝E[x(k)x^H(k)]＝A(θ)R_sA^H(θ)+σ²I_M    (5)

其中，R_s为目标回波信号的协方差矩阵；A(θ)表示目标回波信号的方向矩阵，A^H(θ)表示A(θ)的共轭转置，x(k)表示观测数据矢量，x^H(k)表示观测数据矢量的共轭转置，σ²表示空时高斯白噪声的方差，E[·]表示求取数学期望，I_M表示M维的单位矩阵；

对观测数据的协方差矩阵做特征值分解：

$R_x=\underset{i=1}{\overset{M}{\mathrm{\Σ}}}{\mathrm{\λ}}_i{v}_i{v}_i^H=V_s{\mathrm{\Λ}}_s{V}_s^H+{\mathrm{\σ}}^2{V}_n{V}_n^H---\left(6\right)$

其中，特征值λ₁＞λ₂＞…＞λ_P＞λ_P+1＝…＝λ_M＝σ²，V_s＝[v₁，v₂，…，v_P]，V_n＝[v_P+1，v_P+2，…，v_M]；V_s的列数等于目标回波信号协方差矩阵R_s的秩P，张成A(θ)的P维子空间；v_i表示观测数据协方差矩阵的特征向量，

表示观测数据协方差矩阵的特征向量的共轭转置，V_s表示目标回波信号子空间的特征向量组成的矩阵，

表示V_s的共轭转置，Λ_s表示目标回波信号子空间的特征值组成的对角矩阵，V_n表示噪声子空间的特征向量组成的矩阵，表示V_n的共轭转置；

由公式(5)和公式(6)得到：

V_s＝A(θ)Q    (7)

其中，Q＝R_sA^H(θ)V_s(Λ_s-σ²I_M)^-1，Q为P维满秩矩阵；

CSA-MWF的各级前向滤波器h_i，i＝{1，2，…，P}是相互正交的，i表示第i级前向滤波器，P表示前向滤波器的级数，所以级数为P的相关相减结构的多级维纳滤波器CSA-MWF为Wiener-Hopf方程在Krylov子空间 ${\mathrm{\κ}}^{\left(P\right)}(R_{x_0},r_{x_0{d}_0})=\mathrm{span}\{r_{x_0{d}_0},R_{x_0}{r}_{x_0{d}_0},\·\·\·,R_{x_0}^{(P-1)}{r}_{x_0{d}_0}\}$ 的解，其中x₀表示输入的观测向量，d₀表示期望向量，

表示x₀与d₀的互相关函数，

表示x₀的自相关函数，表示

的(P-1)次方，κ^(P)表示P级的Krylov子空间，span{·}表示将括号中的向量张成空间，目标回波信号子空间按照公式(8)进行估计：

$\mathrm{span}\{h_1,h_2,\·\·\·,h_P\}=\mathrm{span}\{r_{x_0{d}_0},R_{x_0}{r}_{x_0{d}_0},\·\·\·,R_{x_0}^{(P-1)}{r}_{x_0{d}_0}\}---\left(8\right)$

存在一个P阶满秩矩阵K，使得公式(9)成立：

$[h_1,h_2,\·\·\·,h_P]=[r_{x_0{d}_0},R_{x_0}{r}_{x_0{d}_0},\·\·\·,R_{x_0}^{(P-1)}{r}_{x_0{d}_0}]K---\left(9\right)$

令T_s＝[h₁，h₂，…，h_P]，T_n＝[h_P+1，h_P+2，…，h_M]；由于和

其中I_P和I_M-P分别表示维数为P和(M-P)的单位阵，由公式(6)得到：

${R_x}^{\left(i\right)}=V_s{\mathrm{\Λ}}_s^{\left(i\right)}{V}_s^H+{\mathrm{\σ}}^{2i}{V}_n{V}_n^H---\left(10\right)$

其中，R_x ⁽ⁱ⁾表示R_x的i次方，i＝1，2，…，P-；

落在目标回波信号子空间中，

将公式(10)带入公式(9)中，同时

和V_s＝A(θ)Q，得到公式(11)：

$T_s=[V_s{V}_s^H{r}_{x_0{d}_0},V_s{\mathrm{\Λ}}_s{V}_s^H{r}_{x_0{d}_0},\·\·\·,V_s{\mathrm{\Λ}}_s^{(P-1)}{V}_s^H{r}_{x_0{d}_0}]K$

$=V_s[V_s^H{r}_{x_0{d}_0},{\mathrm{\Λ}}_s{V}_s^H{r}_{x_0{d}_0},\·\·\·,{\mathrm{\Λ}}_s^{(P-1)}{V}_s^H{r}_{x_0{d}_0}]K$

$=A\left(\mathrm{\θ}\right)Q[V_s^H{r}_{x_0{d}_0},{\mathrm{\Λ}}_s{V}_s^H{r}_{x_0{d}_0},\·\·\·,{\mathrm{\Λ}}_s^{(P-1)}{V}_s^H{r}_{x_0{d}_0}]K---\left(11\right)$

$=A\left(\mathrm{\θ}\right)\mathrm{Q\ΓK}$

$=A\left(\mathrm{\θ}\right)H$

其中K表示一个P阶满秩矩阵，

表示目标子空间的特征值组成的对角矩阵Λ_s的(P-1)次方，Γ表示表示一个P阶矩阵：

$\mathrm{\Γ}=[V_s^H{r}_{x_0{d}_0},{\mathrm{\Λ}}_s{V}_s^H{r}_{x_0{d}_0},\·\·\·,{\mathrm{\Λ}}_s^{(P-1)}{V}_s^H{r}_{x_0{d}_0}]---\left(12\right)$

H表示一个P阶矩阵，Q＝R_sA^H(θ)V_s(Λ_s-σ²I)^-1：

H＝QΓK    (13)

Γ是满秩矩阵，由于Q和K均是非奇异矩阵，H也是非奇异矩阵；由公式(11)得到P阶的目标回波信号子空间的表达式：

$\mathrm{span}\{h_1,h_2,\·\·\·,h_P\}=\mathrm{col}\left\{A\left(\mathrm{\θ}\right)\right\}\stackrel{\mathrm{\Δ}}{=}{\mathrm{\Φ}}_s^{\left(P\right)}---\left(14\right)$

其中，col{}表示列空间，为其列矢量所有线性组合集合成的空间；

由CSA-MWF的M个相互正交的匹配滤波器构成的预滤波矩阵为T_M＝[h₁，h₂，…h_P，h_P+1，…，h_M]，由于其所有列矢量h_j，j＝1，2，…，M均相互正交，则h_k⊥col{A(θ)}，k＝P+1，P+2，…，M，所以h_k位于由A(θ)的各列向量张成的列空间col{A(θ)}的正交补子空间，为(M-P)维的噪声子空间的表达式

$\mathrm{span}\{h_{P+1},h_{P+2},\·\·\·,h_M\}=\mathrm{null}\left\{A\left(\mathrm{\θ}\right)\right\}\stackrel{\mathrm{\Δ}}{=}{\mathrm{\Φ}}_n^{(M-P)}---\left(15\right)$

其中null{·}表示括号中空间的正交补子空间；

步骤三：按照下式计算CSA-MWF中前向分解的第j级前向滤波器h_j，i＝{1，2，…，P，…，M}：

$h_j=\frac{E\left[d_{j-1}^{*}\left(k\right)x_{j-1}\left(k\right)\right]}{{\left|\right|E\left[d_{j-1}^{*}\left(k\right)x_{j-1}\left(k\right)\right]\left|\right|}_2}---\left(16\right)$

x_j-1(k)表示第(j-1)级前向滤波器的观测数据；

表示第(j-1)级前向滤波器的期望信号d_j-1(k)的共轭信号；

步骤四：按照下式计算CSA-MWF中前向分解的第j级期望信号d_j(k)：

$d_j\left(k\right)=h_j^H{x}_{j-1}\left(k\right)---\left(17\right)$

步骤五：按照下式计算CSA-MWF中前向分解的各级前向滤波器更新后的观测数据：

x_j(k)＝x_j-1(k)-h_jd_j(k)    (18)

其中x_j(k)表示第j级前向滤波器的观测数据；d_j(k)表示第j级前向滤波器的期望信号；

步骤六：进行门限判决，若|x_j(k)|²≤2MNσ²，则进行步骤七，其中M为阵元数，N为快拍数；否则令j＝j+1，返回步骤三计算下一级前向滤波器；

步骤七：将计算得出的CSA-MWF中的各级前向滤波器h₁，h₂，…，h_P带入公式 $\mathrm{span}\{h_1,h_2,\·\·\·,h_P\}=\mathrm{col}\left\{A\left(\mathrm{\θ}\right)\right\}\stackrel{\mathrm{\Δ}}{=}{\mathrm{\Φ}}_s^{\left(P\right)}$ 中，计算得出目标回波信号子空间 ${\mathrm{\Φ}}_s^{\left(P\right)}=\mathrm{span}\{h_1,h_2,\·\·\·,h_P\}.$

2.根据权利要求1所述的一种基于CSA-MWF的无源雷达目标回波信号子空间的估计方法，其特征在于：还包含步骤八，具体为：

(1)令j＝P+1，返回执行步骤三，计算CSA-MWF的第j级前向滤波器，计算得出前向滤波器h_P+1；

(2)计算CSA-MWF中前向分解的第j级期望信号

(3)计算CSA-MWF中前向分解的第j级前向滤波器更新后的观测数据：

x_j(k)＝x_j-1(k)-h_jd_j(k)

其中x_j(k)表示第j级前向滤波器的观测数据；d_j(k)表示第j级前向滤波器的期望信号；

(4)判断j＝M是否成立，若成立，将计算得到的各级前向滤波器h_P+1，h_P+2，…，h_M代入公式 $\mathrm{span}\{h_{P+1},h_{P+2},\·\·\·,h_M\}=\mathrm{null}\left\{A\left(\mathrm{\θ}\right)\right\}\stackrel{\mathrm{\Δ}}{=}{\mathrm{\Φ}}_n^{(M-P)}$ 中，计算得到噪声子空间否则令j＝j+1，返回步骤八(1)，计算下一级前向滤波器。

Images