CN103954931A - 一种远场和近场混合信号源的定位方法 - Google Patents

Abstract

本发明公开了一种远场和近场混合信号源的定位方法，利用斜投影算子，有效地从混合信号中分离出远场信号，得到仅包含近场信号信息的传感器协方差矩阵，进而计算出近场信号的一维方位信息，另外，对于斜投影算子的计算，提出了一种更有效的估计方法，并引入交替迭代过程来消除近场信源定位中出现的饱和现象，本发明降低了计算复杂度，方法简单有效。

Description

一种远场和近场混合信号源的定位方法

技术领域

本发明属于阵列信号处理技术领域，具体涉及一种远场和近场混合信号源的定位方法(LOFNS)。

背景技术

信源定位在雷达、声纳、无线通信等领域具有重要的应用，在一些实际应用，如基于麦克风阵列的话音定位问题中，输入信号为远场和近场混合信号，需要分别估计出远场和近场信号的一维波达方向。尽管已经有很多用于解决远场和近场混合信号源定位问题的方法被提出，如基于高阶统计量的估计方法、基于矩阵差分的估计方法、基于二阶统计量(SOS)的估计方法等。但都存在一定的缺陷和不足，基于高阶统计量的估计方法需要较高的计算复杂度；基于矩阵差分的估计方法通过消除混合信号中的远场信号和噪声，进而估计出近场信号波达方向，此种方法受采样次数的影响较大，仅在大采样数据情况下较为有效；基于SOS的估计方法虽然计算复杂度较低，但是该方法依据主观准则来区分远场和近场信号类型，有一定主观性。

发明内容

本发明的目的在于提供一种计算有效的远场和近场混合信号源的定位方法。

为达到上述目的，本发明采用了以下技术方案。

下文中，对于任意变量a，表示该变量a的估计值。

利用消除噪声后的混合信号估计远场信号的波达方向角，利用斜投影算子从消除噪声后的混合信号中分离出远场信号，得到仅包含近场信号信息的协方差矩阵，根据仅包含近场信号信息的协方差矩阵估计出近场信号的波达方向角及距离，所述混合信号为入射到对称均匀线阵上的K个非相干窄带信号{s_k(n)}，所述对称均匀线阵包含2M+1个全向传感器阵元，M的取值范围为M≥K，阵元间距为d，设前K₁个入射信号为远场信号，远场信号的方位信息为后K₂个入射信号为近场信号，近场信号的方位信息为K＝K₁+K₂，θ_k表示第k个入射信号的波达方向角，所述波达方向角为第k个入射信号相对于y轴的逆时针夹角，r_k是第k个入射信号相对于坐标原点的距离。

所述远场信号的波达方向角的估计方法包括以下步骤：

1)消除阵列协方差矩阵R中的噪声部分，得到无噪协方差矩阵将分为两部分：

2)构造代价函数

$g_f\left(\mathrm{\θ}\right)=a_f^H\left(\mathrm{\θ}\right){\mathrm{\Π}}_{\hat{Q}}{a}_f\left(\mathrm{\θ}\right)---\left(1\right)$

式(1)中(·)^H表示共轭转置；

${\mathrm{\Π}}_{\hat{Q}}=\hat{Q}{\left({\hat{Q}}^H\hat{Q}\right)}^{-1}{\hat{Q}}^H=\hat{Q}(I_{2M+1-K}-{\hat{P}}^H{(\hat{P}{\hat{P}}^H+I_K)}^{-1}\hat{P}){\hat{Q}}^H$

$\hat{Q}={[{\hat{P}}^T,{-I}_{2M+1-K}]}^T$

$\hat{P}={\left({\hat{G}}^H\hat{G}\right)}^{-1}{\hat{G}}^H\hat{H}$

(·)^T表示转置，I_m表示m×m的单位矩阵；

3)根据式(1)构造多项式其中 $p\left(z\right)\stackrel{\mathrm{\Δ}}{=}{[z^M,...,z,1,z^{-1},...,z^{-M}]}^T,z\stackrel{\mathrm{\Δ}}{=}{e}^{j2\mathrm{\π}d\sin \mathrm{\θ}/\mathrm{\λ}},$ λ表示入射信号波长，j表示单位虚数，j²＝-1，通过求多项式的K₁个零相位点来估计远场信号的波达方向角。

所述步骤1)具体包括以下步骤：

a、根据对称均匀线阵接收的数据求得阵列协方差矩阵R的估计值：

$R=\frac{1}{N}\underset{n=1}{\overset{N}{\mathrm{\Σ}}}x\left(n\right)x^H\left(n\right)---\left(2\right)$

其中，N表示采样数，x(n)表示阵列接收数据；

b、将阵列协方差矩阵R划分为如下形式

c、估计噪声方差σ²

${\mathrm{\σ}}^2=\frac{\mathrm{tr}\left\{R_{22}\mathrm{\Π}\right\}}{\mathrm{tr}\left\{\mathrm{\Π}\right\}}---\left(4\right)$

其中tr{·}表示求矩阵的迹，表示广义逆；

d、无噪协方差矩阵

所述近场信号的波达方向角及距离的估计方法包括以下步骤：

2.1)估计斜投影算子

2.2)消除无噪协方差矩阵中的远场信号，得到仅包含近场信号信息的协方差矩阵R_n， $R_n=(I_{2M+1}-E_{A_f|A_n})\stackrel{\‾}{R}{(I_{2M+1}-E_{A_f|A_n})}^H;$

2.3)利用仅包含近场信号信息的协方差矩阵R_n的反对角元素信息构造一个新的协方差矩阵并将新构造的协方差矩阵分块：

2.4)构造代价函数

$f_n={\stackrel{\‾}{a}}_n^H\left(\mathrm{\θ}\right){\mathrm{\Π}}_{\widehat{\stackrel{\‾}{Q}}}{\stackrel{\‾}{a}}_n\left(\mathrm{\θ}\right)---\left(12\right)$

式(12)中， ${\mathrm{\Π}}_{\widehat{\stackrel{\‾}{Q}}}=\widehat{\stackrel{\‾}{Q}}{\left({\widehat{\stackrel{\‾}{Q}}}^H\widehat{\stackrel{\‾}{Q}}\right)}^{-1}{\widehat{\stackrel{\‾}{Q}}}^H=\widehat{\stackrel{\‾}{Q}}(I_{M+1-K_2}-{\widehat{\stackrel{\‾}{P}}}^H{(\widehat{\stackrel{\‾}{P}}{\widehat{\stackrel{\‾}{P}}}^H+I_{K_2})}^{-1}\widehat{\stackrel{\‾}{P}}){\widehat{\stackrel{\‾}{Q}}}^H,\widehat{\stackrel{\‾}{P}}={\left({\widehat{\stackrel{\‾}{G}}}^H\widehat{\stackrel{\‾}{G}}\right)}^{-1}{\widehat{\stackrel{\‾}{G}}}^H\widehat{\stackrel{\‾}{H}},$ $\widehat{\stackrel{\‾}{Q}}={[{\widehat{\stackrel{\‾}{P}}}^T,-I_{M+1-K_2}]}^T;$

2.5)根据式(12)构造多项式 $p_N\left(z\right)\stackrel{\mathrm{\Δ}}{=}{z}^M{p}^H\left(z\right){\mathrm{\Π}}_{\widehat{\stackrel{\‾}{Q}}}p\left(z\right),$ 其中 $p\left(z\right)\stackrel{\mathrm{\Δ}}{=}{[1,z^{-1},...,z^{-M}]}^T,$ $z\stackrel{\mathrm{\Δ}}{=}{e}^{j4\mathrm{\π}d\sin \mathrm{\θ}/\mathrm{\λ}},$ 通过求多项式 $p_N\left(z\right)\stackrel{\mathrm{\Δ}}{=}{z}^M{p}^H\left(z\right){\mathrm{\Π}}_{\widehat{\stackrel{\‾}{Q}}}p\left(z\right)$ 的K2个零相位点来估计近场信号的波达方向角；

2.6)构造代价函数

$\begin{array}{c}{\stackrel{\‾}{f}}_N\left(r\right)=a_n^H(r,\widehat{\mathrm{\θ}}){\mathrm{\Π}}_{\hat{Q}}{a}_n(r,\widehat{\mathrm{\θ}})\\ ={\mathrm{\μ}}_n^H(\widehat{\mathrm{\θ}},r)D^H\left(\widehat{\mathrm{\θ}}\right){\mathrm{\Π}}_{\hat{Q}}D\left(\widehat{\mathrm{\θ}}\right){\mathrm{\μ}}_n(\widehat{\mathrm{\θ}},r)\end{array}---\left(13\right)$

其中 ${\mathrm{\μ}}_n(\widehat{\mathrm{\θ}},r)={[e^{M^2\mathrm{\φ}},e^{{(M-1)}^2\mathrm{\φ}},...,1]}^T,\mathrm{\φ}\stackrel{\mathrm{\Δ}}{=}{\mathrm{\πd}}^2{\cos }^2\widehat{\mathrm{\θ}}/\left(\mathrm{\λr}\right),D\left(\widehat{\mathrm{\θ}}\right)\stackrel{\mathrm{\Δ}}{=}{[\mathrm{\Φ}\left(\widehat{\mathrm{\θ}}\right),{\mathrm{\Φ}}^{*}\left(\widehat{\mathrm{\θ}}\right)J]}^T,$ 是块对角矩阵， $\mathrm{\Φ}\left(\widehat{\mathrm{\θ}}\right)=\mathrm{diag}\{e^{-j\widehat{\mathrm{\ω}}M},e^{-j\widehat{\mathrm{\ω}}(M-1)},...,1\},$ (·)^*表示取共轭， $\widehat{\mathrm{\ω}}\stackrel{\mathrm{\Δ}}{=}-2\mathrm{\π}d\sin \widehat{\mathrm{\θ}}/\mathrm{\λ},$ J是一个(M+1)×M的选择矩阵，定义为J＝[e_M,e_M-1,…,e₁]，e_m表示单位矩阵I_M+1的第m列，m＝1,2,…,M；

2.7)根据式(13)构造多项式 $p_r\left(z\right)\stackrel{\mathrm{\Δ}}{=}{z}^{M^2}{p}^H\left(z\right)D^H\left(\widehat{\mathrm{\θ}}\right){\mathrm{\Π}}_{\hat{Q}}D\left(\widehat{\mathrm{\θ}}\right)p\left(z\right),$ 其中 $p\left(z\right)\stackrel{\mathrm{\Δ}}{=}{[1,z,...,z^{M^2}]}^T,z\stackrel{\mathrm{\Δ}}{=}{e}^{\mathrm{\π}{d}^2{\cos }^2\widehat{\mathrm{\θ}}/\left(\mathrm{\λr}\right)},$ 通过求多项式 $p_r\left(z\right)\stackrel{\mathrm{\Δ}}{=}{z}^{M^2}{p}^H\left(z\right)D^H\left(\widehat{\mathrm{\θ}}\right){\mathrm{\Π}}_{\hat{Q}}D\left(\widehat{\mathrm{\θ}}\right)p\left(z\right)$ 的K₂个零相位点来估计近场信号的距离；

2.8)更新斜投影算子返回步骤2.2)，继续顺序执行，如此往复循环，直至满足精度要求时结束。

所述步骤2.1)中斜投影算子的估计方法包括以下步骤：

a、构造一个(2M+1)×(2M+1)矩阵：

$R_o\stackrel{\mathrm{\Δ}}{=}\stackrel{\‾}{R}{\mathrm{\Π}}_{A_f}^{\⊥}---\left(14\right)$

其中A_f是远场信号方向矩阵，由估计出的远场信号的波达方向角计算得出；

b、对R_oΔ进行QR分解得斜投影算子由下述方式获得

其中是矩阵的后2M+1-K₂列，Δ是(2M+1)×(2M+1)的置换矩阵，不改变矩阵R_o的列相关性，表示广义逆。

所述步骤2.3)具体包括以下步骤：

仅包含近场信号信息的协方差矩阵R_n的反对角元素r_n(i)定义为

$r_n\left(i\right)\stackrel{\mathrm{\Δ}}{=}{\left(R_n\right)}_{i,2M+2-i}=\underset{k=K_1+1}{\overset{K}{\mathrm{\Σ}}}{r}_{\mathrm{sk}}^2{e}^{-j2(M+1-i){\mathrm{\ω}}_k}---\left(16\right)$

i＝1,2,…,2M+1，是第k个近场信号的功率；

更新斜投影算子的方法为：

$E_{A_f|A_n}=A_f{\left(A_f^H{\mathrm{\Π}}_{A_n}^{\⊥}{A}_f\right)}^{-1}{A}_f^H{\mathrm{\Π}}_{A_n}^{\⊥}---\left(25\right)$

A_f是远场信号方向矩阵，由估计出的远场信号的波达方向角计算得出。

本发明的有益效果是：

本发明对远场和近场信号方位信息分别进行估计，利用斜投影算子，有效地从混合信号中分离出远场信号，得到仅包含近场信号信息的传感器协方差矩阵，进而计算出近场信号的波达方向角及距离，相比于已有的远场和近场混合信号源定位算法，本发明计算复杂度低，方法简单有效，若以matlab每秒浮点运算次数来描述算法复杂度，本发明的复杂度为每秒运算10(2M+1)²+30(2M+1)³次浮点数。

另外，本发明提出了一种可以更有效地计算斜投影算子的方法，并引入交替迭代方法来消除近场信源定位中出现的饱和现象。

附图说明

图1为计算复杂度随采样数及阵元个数变化曲线：(a)为随采样数，阵元个数设为7，(b)为随阵元数，采样数设为100；虚线：LOFNS(无迭代)；实线：EAMA；实线带有“×”：MSLA；虚线加“o”：TSA；点划线：EMPL。

图2为波达方向角的估计均方根误差随信噪比变化曲线：(a)为远场信号，(b)为近场信号；点线加“Δ”：LOFNS；虚线：LOFNS(无迭代)；实线：EAMA；实线带有“×”：MSLA；虚线加“o”：TSA；点划线：EMPL；点线：CRB。

图3为波达方向角的估计均方根误差随采样数变化曲线：(a)为远场信号，(b)为近场信号；点线加“Δ”：LOFNS；虚线：LOFNS(无迭代)；实线：EAMA；实线带有“×”：MSLA；虚线加“o”：TSA；点划线：EMPL；点线：CRB。

图4为近场信号距离的估计均方根误差随信噪比及采样数的变化曲线：(a)为随信噪比，(b)为随采样数；点线加“Δ”：LOFNS；虚线：LOFNS(无迭代)；实线：EAMA；实线带有“×”：MSLA；虚线加“o”：TSA；点划线：EMPL；点线：CRB。

图5为阵列结构图。

具体实施方式

下面结合附图和实施例对本发明做详细描述。

一种远场和近场混合信号源的定位方法，具体实现步骤概括如下：

1)计算阵列协方差矩阵R的估计值

2)消除阵列协方差矩阵中的噪声部分，得到无噪协方差矩阵

3)由无噪协方差矩阵构造代价函数，通过求解优化问题估计出远场信号的波达方向角 ${\left\{{\mathrm{\θ}}_k\right\}}_{k=1}^{K_1};$

4)计算斜投影算子的初始值；

5)利用斜投影算子消除无噪协方差矩阵中的远场信号部分，得到仅包含近场信息的新的协方差矩阵

6)由新的协方差矩阵构造代价函数，通过求解优化问题估计出近场信号的波达方向角 ${\left\{{\mathrm{\θ}}_k\right\}}_{k=K_1+1}^K$ 及距离 ${\left\{r_k\right\}}_{k=K_1+1}^K;$

7)将上述求得的远场信号和近场信号的波达方向角以及近场信号的距离作为初始值，重新估计斜投影算子，返回步骤5)，反复迭代，直到满足条件。

下面进行具体描述。

K个非相干窄带信号{s_k(n)}入射到对称均匀线阵ULA上，该对称均匀线阵包含2M+1个全向传感器阵元，M的取值范围为M≥K，阵元间距为d。设前K₁个入射信号为远场信号，其方位信息为后K₂个入射信号为近场信号，其方位信息为并且K＝K₁+K₂。θ_k表示第k个入射信号(远场信号或近场信号)相对于y轴的逆时针夹角(波达方向角)，r_k是第k个入射信号(仅近场信号)相对于坐标原点的距离，参见图5。

令对称均匀线阵ULA的中心为参考阵元，阵列输出信号为

x(n)＝A_fs_f(n)+A_ns_n(n)+w(n)

                   (1)

＝As(n)+w(n)

其中，Α是阵列响应矩阵， $A\stackrel{\mathrm{\Δ}}{=}[A_f,A_n],A_f\stackrel{\mathrm{\Δ}}{=}[a_f\left({\mathrm{\θ}}_1\right),a_f\left({\mathrm{\θ}}_2\right),...,a_f\left({\mathrm{\θ}}_{K_1}\right)],$ $A_n\stackrel{\mathrm{\Δ}}{=}[a_n(r_{K_1+1},{\mathrm{\θ}}_{K_1+1}),a_n(r_{K_1+2},{\mathrm{\θ}}_{K_1+2}),...,a_n(r_K,{\mathrm{\θ}}_K)],$ 远场和近场的方向向量分别定义为 $a_f\left({\mathrm{\θ}}_k\right)\stackrel{\mathrm{\Δ}}{=}{[e^{-\mathrm{jM}{\mathrm{\ω}}_k},...,e^{-j{\mathrm{\ω}}_k},1,e^{j{\mathrm{\ω}}_k},...,e^{\mathrm{jM}{\mathrm{\ω}}_k}]}^T,a_n(r_k,{\mathrm{\θ}}_k)\stackrel{\mathrm{\Δ}}{=}{[e^{j{\mathrm{\τ}}_{-\mathrm{Mk}}},...,e^{j{\mathrm{\τ}}_{-k}},1,e^{j{\mathrm{\τ}}_k},...,e^{j{\mathrm{\τ}}_{\mathrm{Mk}}}]}^T,$ (·)^T表示转置。远场信号的相位延迟ω_k定义为λ是入射信号波长，j表示单位虚数，j²＝-1。当近场距离满足r_f∈(0.62(D³/λ)^1/2,2D²/λ)，其中D是阵列孔径，根据菲涅尔近似，近场信号的相位延迟τ_mk定义为其中m＝-M,…,-1,0,1,…,M， ${\mathrm{\φ}}_k\stackrel{\mathrm{\Δ}}{=}{\mathrm{\πd}}^2{\cos }^2{\mathrm{\θ}}_k/\left(\mathrm{\λ}{r}_k\right).$

接收数据的阵列协方差矩阵为

$\begin{array}{c}R\stackrel{\mathrm{\Δ}}{=}E\left\{x\left(n\right)x^H\left(n\right)\right\}\\ =A_f{R}_{\mathrm{sf}}{A}_f^H+A_n{R}_{\mathrm{sn}}{A}_n^H+{\mathrm{\σ}}^2{I}_{2M+1}\\ =\stackrel{\‾}{R}+{\mathrm{\σ}}^2{I}_{2M+1}\end{array}---\left(2\right)$

其中是无噪协方差矩阵，定义为R_sf和R_sn分别是远场和近场信号协方差矩阵，定义为 $R_{\mathrm{sf}}\stackrel{\mathrm{\Δ}}{=}E\left\{s_f\left(n\right)s_f^H\left(n\right)\right\},R_{\mathrm{sn}}\stackrel{\mathrm{\Δ}}{=}E\left\{s_n\left(n\right)s_n^H\left(n\right)\right\}.$ 另外E{w(n)w^H(n)}＝σ²I_2M+1，blkdiag(·)表示构建分块对角矩阵，(·)^H表示矩阵共轭转置，I_m是一个m×m单位阵。

下面描述远场信号波达方向角的估计方法，包括以下步骤：

1)消除阵列协方差矩阵R中的噪声部分，得到无噪协方差矩阵

2)无噪协方差矩阵可分为两部分，

3)构造代价函数

$g_f\left(\mathrm{\θ}\right)=a_f^H\left(\mathrm{\θ}\right){\mathrm{\Π}}_{\hat{Q}}{a}_f\left(\mathrm{\θ}\right)---\left(3\right)$

4)根据上式构造多项式其中 $p\left(z\right)\stackrel{\mathrm{\Δ}}{=}{[z^M,...,z,1,z^{-1},...,z^{-M}]}^T,z\stackrel{\mathrm{\Δ}}{=}{e}^{j2\mathrm{\π}d\sin \mathrm{\θ}/\mathrm{\λ}},$ 通过求该多项式的K₁个零相位点来估计远场信号的波达方向角θ_k；

所述步骤1)消除阵列协方差矩阵R中的噪声部分，得到无噪协方差矩阵的过程包括：

a、根据阵列(对称均匀线阵)接收数据求得阵列协方差矩阵R的估计值：

$R=\frac{1}{N}\underset{n=1}{\overset{N}{\mathrm{\Σ}}}x\left(n\right)x^H\left(n\right)---\left(4\right)$

其中，N表示采样数，x(n)表示阵列接收数据；

b、阵列协方差矩阵R可以划分为如下形式

c、估计噪声方差σ²

${\mathrm{\σ}}^2=\frac{\mathrm{tr}\left\{R_{22}\mathrm{\Π}\right\}}{\mathrm{tr}\left\{\mathrm{\Π}\right\}}---\left(6\right)$

其中tr{·}表示求矩阵的迹，表示广义逆；

d、无噪协方差矩阵

所述步骤3)中

${\mathrm{\Π}}_{\hat{Q}}=\hat{Q}{\left({\hat{Q}}^H\hat{Q}\right)}^{-1}{\hat{Q}}^H=\hat{Q}(I_{2M+1-K}-{\hat{P}}^H){(\hat{P}{\hat{P}}^H+I_K)}^{-1}\hat{P}){\hat{Q}}^H---\left(7\right)$

$\hat{Q}={[{\hat{P}}^T,{-I}_{2M+1-K}]}^T---\left(8\right)$

$\hat{P}={\left({\hat{G}}^H\hat{G}\right)}^{-1}{\hat{G}}^H\hat{H}---\left(9\right)$

是一个线性变化矩阵，阵列响应矩阵A可以被分为两个子阵，如下

基于A是满秩矩阵这一基本假设，矩阵A的前K行线性独立，因此，存在一个K×(2M+1-K)的线性变化矩阵P，使得

A₂＝P^HA₁，i.e.Q^HA＝O_(2M+1-K)×K  (11)

其中O_p×q是一个p×q的零矩阵，由(11)式可以很容易得到

Q^Ha_n(r,θ)＝0_(2M+1-K)×1  (12)

Q^Ha_f(θ)＝0_(2M+1-K)×1  (13)

下面描述近场信号波达方向角及距离估计。

1)估计斜投影算子

2)消除无噪协方差矩阵中的远场信号，得到仅包含近场信号信息的协方差矩阵R_n， $R_n=(I_{2M+1}-E_{A_f|A_n})\stackrel{\‾}{R}{(I_{2M+1}-E_{A_f|A_n})}^H;$

3)利用仅包含近场信号信息的协方差矩阵R_n的反对角元素信息，构造新的协方差矩阵

4)新的协方差矩阵可分为两部分：

5)构造代价函数

$f_n={\stackrel{\‾}{a}}_n^H\left(\mathrm{\θ}\right){\mathrm{\Π}}_{\widehat{\stackrel{\‾}{Q}}}{\stackrel{\‾}{a}}_n\left(\mathrm{\θ}\right)---\left(14\right)$

6)根据上式构造多项式 $p_N\left(z\right)\stackrel{\mathrm{\Δ}}{=}{z}^M{p}^H\left(z\right){\mathrm{\Π}}_{\widehat{\stackrel{\‾}{Q}}}p\left(z\right),$ 其中 $p\left(z\right)\stackrel{\mathrm{\Δ}}{=}{[1,z^{-1},...,z^{-M}]}^T,$ 通过求该多项式的K₂个零相位点来估计近场信号的波达方向角θ_k；

7)构造代价函数

$\begin{array}{c}{\stackrel{\‾}{f}}_N\left(r\right)=a_n^H(r,\widehat{\mathrm{\θ}}){\mathrm{\Π}}_{\hat{Q}}{a}_n(r,\widehat{\mathrm{\θ}})\\ ={\mathrm{\μ}}_n^H(\widehat{\mathrm{\θ}},r)D^H\left(\widehat{\mathrm{\θ}}\right){\mathrm{\Π}}_{\hat{Q}}D\left(\widehat{\mathrm{\θ}}\right){\mathrm{\μ}}_n(\widehat{\mathrm{\θ}},r)\end{array}---\left(15\right)$

其中为步骤6)估计的近场信号波达方向角， $\mathrm{\φ}\stackrel{\mathrm{\Δ}}{=}{\mathrm{\πd}}^2{\cos }^2\widehat{\mathrm{\θ}}/\left(\mathrm{\λr}\right),$ 定义为 $D\left(\widehat{\mathrm{\θ}}\right)\stackrel{\mathrm{\Δ}}{=}{[\mathrm{\Φ}\left(\widehat{\mathrm{\θ}}\right),{\mathrm{\Φ}}^{*}\left(\widehat{\mathrm{\θ}}\right)J]}^T,$ 是块对角矩阵， $\mathrm{\Φ}\left(\widehat{\mathrm{\θ}}\right)=\mathrm{diag}\{e^{-j\widehat{\mathrm{\ω}}M},e^{-j\widehat{\mathrm{\ω}}(M-1)},...,1\},\widehat{\mathrm{\ω}}\stackrel{\mathrm{\Δ}}{=}-2\mathrm{\π}d\sin \widehat{\mathrm{\θ}}/\mathrm{\λ},$ J是一个(M+1)×M的选择矩阵，定义为J＝[e_M,e_M-1,…,e₁]，e_m表示单位矩阵I_M+1的第m列，m＝1,2,…,M；

8)根据上式构造多项式其中 $p\left(z\right)\stackrel{\mathrm{\Δ}}{=}{[1,z,...,z^{M^2}]}^T,z\stackrel{\mathrm{\Δ}}{=}{e}^{\mathrm{\π}{d}^2{\cos }^2\widehat{\mathrm{\θ}}/\left(\mathrm{\λr}\right)},$ 通过求该多项式的K₂个零相位点来估计近场信号的距离；

9)更新斜投影算子返回步骤2)，继续顺序执行，如此往复循环，直至满足精度要求时结束。

所述步骤1)中斜投影算子的估计方法包括以下步骤：

a、构造一个(2M+1)×(2M+1)矩阵如下

$R_o\stackrel{\mathrm{\Δ}}{=}\stackrel{\‾}{R}{\mathrm{\Π}}_{A_f}^{\⊥}=A_n{R}_{\mathrm{sn}}{A}_n^H{\mathrm{\Π}}_{A_f}^{\⊥}---\left(16\right)$

其中A_f是远场信号方向矩阵，由以上估计出的远场信号波达方向角计算得出；

b、对R_oΔ进行QR分解得斜投影算子可以由下述方式获得

其中是矩阵的后2M+1-K₂列，Δ是(2M+1)×(2M+1)的置换矩阵，不改变矩阵R_o的列相关性。

所述步骤3)具体叙述如下：

仅包含近场信号信息的协方差矩阵R_n的反对角元素r_n(i)定义为

$r_n\left(i\right)\stackrel{\mathrm{\Δ}}{=}{\left(R_n\right)}_{i,2M+2-i}=\underset{k=K_1+1}{\overset{K}{\mathrm{\Σ}}}{r}_{\mathrm{sk}}^2{e}^{-j2(M+1-i){\mathrm{\ω}}_k}---\left(18\right)$

i可取i＝1,2,…,2M+1，是第k个近场信号的功率，定义为(·)^*表示取共轭。由(18)式可以看到反对角元素r_n(i)只包含近场信号波达方向角度信息。类似于远场Toepliz协方差矩阵，近场信号协方差矩阵如下

将(18)式带入(19)式，经过简单的代数运算，重写协方差矩阵如下式

${\stackrel{\‾}{R}}_n={\stackrel{\‾}{A}}_n\left(\mathrm{\θ}\right)R_{\mathrm{sn}}{\stackrel{\‾}{A}}_n^H\left(\mathrm{\θ}\right)---\left(20\right)$

其中 ${\stackrel{\‾}{A}}_n\left(\mathrm{\θ}\right)\stackrel{\mathrm{\Δ}}{=}[{\stackrel{\‾}{a}}_n\left({\mathrm{\θ}}_1\right),{\stackrel{\‾}{a}}_n\left({\mathrm{\θ}}_2\right),...,{\stackrel{\‾}{a}}_n\left({\mathrm{\θ}}_{K_1}\right)],{\stackrel{\‾}{a}}_n\left({\mathrm{\θ}}_k\right)\stackrel{\mathrm{\Δ}}{=}{\left[1,e^{j2{\mathrm{\ω}}_k},e^{j4{\mathrm{\ω}}_k},...,e^{j2M{\mathrm{\ω}}_k}\right]}^T.$

所述步骤5)中

${\mathrm{\Π}}_{\hat{Q}}=\hat{Q}{\left({\hat{Q}}^H\hat{Q}\right)}^{-1}{\hat{Q}}^H=\hat{Q}(I_{2M+1-K}-{\hat{P}}^H{(\hat{P}{\hat{P}}^H+I_K)}^{-1}\hat{P}){\hat{Q}}^H---\left(21\right)$

$\hat{Q}={[{\hat{P}}^T,{-I}_{2M+1-K}]}^T---\left(22\right)$

$\widehat{\stackrel{\‾}{P}}={\left({\widehat{\stackrel{\‾}{G}}}^H\widehat{\stackrel{\‾}{G}}\right)}^{-1}{\widehat{\stackrel{\‾}{G}}}^H\widehat{\stackrel{\‾}{H}}---\left(23\right)$

是一个线性变化矩阵，可以被分为两个子阵，如下

其中是一个满秩矩阵，的行向量可以由的行向量线性表示，因此，存在一个K₂×(M+1-K₂)的线性算子使得

${\stackrel{\‾}{A}}_{n2}\left(\mathrm{\θ}\right)={\stackrel{\‾}{P}}^H{\stackrel{\‾}{A}}_{n1}\left(\mathrm{\θ}\right),i.e.{\stackrel{\‾}{Q}}^H{\stackrel{\‾}{A}}_n\left(\mathrm{\θ}\right)=O_{(M+1-K_2)\×K_2}---\left(25\right)$

其中O_p×q是一个p×q的零矩阵。

在步骤8)完成后更新斜投影算子

$E_{A_f|A_n}=A_f{\left(A_f^H{\mathrm{\Π}}_{A_n}^{\⊥}{A}_f\right)}^{-1}{A}_f^H{\mathrm{\Π}}_{A_n}^{\⊥}---\left(26\right)$

更新斜投影算子的目的在于得到更为精确的的估计值。步骤1)中的估计值有一定误差，因为该估计过程需要利用信号协方差矩阵R_s的块对角结构和仅包含近场信号信息的协方差矩阵R_n的反对角元素。然而，在采样数有限的情况下，R_s并非块对角结构，R_n的反对角元素也并非仅包含近场信号信息，可以由下式看到

${\hat{R}}_o=A_n{\hat{R}}_{\mathrm{sn}}{A}_n^H{\mathrm{\Π}}_{A_f}^{\⊥}+A_f{\hat{R}}_{\mathrm{sfn}}{A}_n^H{\mathrm{\Π}}_{A_f}^{\⊥}---\left(27\right)$

其中由于的存在，的列空间与A_n的列空间不完全相等。因此步骤1)中斜投影算子的估计并不准确，并且不受叠加噪声的影响，这意味着尽管信噪比趋于无穷大，的影响仍然存在，估计误差并不随信噪比SNR的增加而减小，称为“饱和现象”。而利用第一次计算得到的远场信号和近场信号波达方向角以及近场信号距离，重新计算斜投影算子，可以进一步优化结果，并且通过迭代，消除“饱和现象“。

下面通过以下不同情形对上述方法的效果进行说明：

空间有两个波达方向未知的入射信号，其方向信息分别为(∞,-5°)，(1.7λ,25°)，对称均匀线阵含有7个阵元，阵元间隔为d＝λ/4。仿真中加入基于矩阵差分的算法如EMPL、TSA，基于SOS的算法如EAMA、MSLA与LOFNS对比，同时给出了CRB界。每一个仿真结果都是经由1000次独立重复实验得到的。

由图1说明LOFNS计算复杂度最低。

由图2及图4(a)说明对于远场信号波达方向角的估计，LOFNS接近EAMA且优于EAMA，特别是在高信噪比情况下。并且通过迭代，消除了饱和现象。

由图3及图4(b)说明LOFNS在远场信号和近场信号的估计性能都优于其他方法。

Claims (7)

1.一种远场和近场混合信号源的定位方法，其特征在于，包括以下步骤：利用消除噪声后的混合信号估计远场信号的波达方向角，利用斜投影算子从消除噪声后的混合信号中分离出远场信号，得到仅包含近场信号信息的协方差矩阵，根据仅包含近场信号信息的协方差矩阵估计出近场信号的波达方向角及距离，所述混合信号为入射到对称均匀线阵上的K个非相干窄带信号{s_k(n)}，所述对称均匀线阵包含2M+1个全向传感器阵元，M的取值范围为M≥K，阵元间距为d，设前K₁个入射信号为远场信号，远场信号的方位信息为后K₂个入射信号为近场信号，近场信号的方位信息为K＝K₁+K₂，θ_k表示第k个入射信号的波达方向角，所述波达方向角为第k个入射信号相对于y轴的逆时针夹角，r_k是第k个入射信号相对于坐标原点的距离。

2.如权利要求1所述一种远场和近场混合信号源的定位方法，其特征在于，所述远场信号的波达方向角的估计方法包括以下步骤：

1)消除阵列协方差矩阵R中的噪声部分，得到无噪协方差矩阵将分为两部分：

2)构造代价函数

$g_f\left(\mathrm{\θ}\right)=a_f^H\left(\mathrm{\θ}\right){\mathrm{\Π}}_{\hat{Q}}{a}_f\left(\mathrm{\θ}\right)---\left(1\right)$

式(1)中(·)^H表示共轭转置；

${\mathrm{\Π}}_{\hat{Q}}=\hat{Q}{\left({\hat{Q}}^H\hat{Q}\right)}^{-1}{\hat{Q}}^H=\hat{Q}(I_{2M+1-K}-{\hat{P}}^H{(\hat{P}{\hat{P}}^H+I_K)}^{-1}\hat{P}){\hat{Q}}^H$

$\hat{Q}={[{\hat{P}}^T,{-I}_{2M+1-K}]}^T$

$\hat{P}={\left({\hat{G}}^H\hat{G}\right)}^{-1}{\hat{G}}^H\hat{H}$

(·)^T表示转置，I_m表示m×m的单位矩阵；

3)根据式(1)构造多项式其中 $p\left(z\right)\stackrel{\mathrm{\Δ}}{=}{[z^M,...,z,1,z^{-1},...,z^{-M}]}^T,z\stackrel{\mathrm{\Δ}}{=}{e}^{j2\mathrm{\π}d\sin \mathrm{\θ}/\mathrm{\λ}},$ λ表示入射信号波长，j表示单位虚数，j²＝-1，通过求多项式的K₁个零相位点来估计远场信号的波达方向角。

3.如权利要求2所述一种远场和近场混合信号源的定位方法，其特征在于，所述步骤1)具体包括以下步骤：

a、根据对称均匀线阵接收的数据求得阵列协方差矩阵R的估计值：

$R=\frac{1}{N}\underset{n=1}{\overset{N}{\mathrm{\Σ}}}x\left(n\right)x^H\left(n\right)---\left(2\right)$

其中，N表示采样数，x(n)表示阵列接收数据；

b、将阵列协方差矩阵R划分为如下形式

c、估计噪声方差σ²

${\mathrm{\σ}}^2=\frac{\mathrm{tr}\left\{R_{22}\mathrm{\Π}\right\}}{\mathrm{tr}\left\{\mathrm{\Π}\right\}}---\left(4\right)$

其中tr{·}表示求矩阵的迹，表示广义逆；

d、无噪协方差矩阵

4.如权利要求2所述一种远场和近场混合信号源的定位方法，其特征在于，所述近场信号的波达方向角及距离的估计方法包括以下步骤：

2.1)估计斜投影算子

2.2)消除无噪协方差矩阵中的远场信号，得到仅包含近场信号信息的协方差矩阵R_n， $R_n=(I_{2M+1}-E_{A_f|A_n})\stackrel{\‾}{R}{(I_{2M+1}-E_{A_f|A_n})}^H;$

2.3)利用仅包含近场信号信息的协方差矩阵R_n的反对角元素信息构造一个新的协方差矩阵并将新构造的协方差矩阵分块：

2.4)构造代价函数

$f_n={\stackrel{\‾}{a}}_n^H\left(\mathrm{\θ}\right){\mathrm{\Π}}_{\widehat{\stackrel{\‾}{Q}}}{\stackrel{\‾}{a}}_n\left(\mathrm{\θ}\right)---\left(12\right)$

式(12)中， ${\mathrm{\Π}}_{\widehat{\stackrel{\‾}{Q}}}=\widehat{\stackrel{\‾}{Q}}{\left({\widehat{\stackrel{\‾}{Q}}}^H\widehat{\stackrel{\‾}{Q}}\right)}^{-1}{\widehat{\stackrel{\‾}{Q}}}^H=\widehat{\stackrel{\‾}{Q}}(I_{M+1-K_2}-{\widehat{\stackrel{\‾}{P}}}^H{(\widehat{\stackrel{\‾}{P}}{\widehat{\stackrel{\‾}{P}}}^H+I_{K_2})}^{-1}\widehat{\stackrel{\‾}{P}}){\widehat{\stackrel{\‾}{Q}}}^H,\widehat{\stackrel{\‾}{P}}={\left({\widehat{\stackrel{\‾}{G}}}^H\widehat{\stackrel{\‾}{G}}\right)}^{-1}{\widehat{\stackrel{\‾}{G}}}^H\widehat{\stackrel{\‾}{H}},$ $\widehat{\stackrel{\‾}{Q}}={[{\widehat{\stackrel{\‾}{P}}}^T,-I_{M+1-K_2}]}^T;$

2.5)根据式(12)构造多项式 $p_N\left(z\right)\stackrel{\mathrm{\Δ}}{=}{z}^M{p}^H\left(z\right){\mathrm{\Π}}_{\widehat{\stackrel{\‾}{Q}}}p\left(z\right),$ 其中 $p\left(z\right)\stackrel{\mathrm{\Δ}}{=}{[1,z^{-1},...,z^{-M}]}^T,$ $z\stackrel{\mathrm{\Δ}}{=}{e}^{j4\mathrm{\π}d\sin \mathrm{\θ}/\mathrm{\λ}},$ 通过求多项式 $p_N\left(z\right)\stackrel{\mathrm{\Δ}}{=}{z}^M{p}^H\left(z\right){\mathrm{\Π}}_{\widehat{\stackrel{\‾}{Q}}}p\left(z\right)$ 的K₂个零相位点来估计近场信号的波达方向角；

2.6)构造代价函数

$\begin{array}{c}{\stackrel{\‾}{f}}_N\left(r\right)=a_n^H(r,\widehat{\mathrm{\θ}}){\mathrm{\Π}}_{\hat{Q}}{a}_n(r,\widehat{\mathrm{\θ}})\\ ={\mathrm{\μ}}_n^H(\widehat{\mathrm{\θ}},r)D^H\left(\widehat{\mathrm{\θ}}\right){\mathrm{\Π}}_{\hat{Q}}D\left(\widehat{\mathrm{\θ}}\right){\mathrm{\μ}}_n(\widehat{\mathrm{\θ}},r)\end{array}---\left(13\right)$

其中 ${\mathrm{\μ}}_n(\widehat{\mathrm{\θ}},r)={[e^{M^2\mathrm{\φ}},e^{{(M-1)}^2\mathrm{\φ}},...,1]}^T,\mathrm{\φ}\stackrel{\mathrm{\Δ}}{=}{\mathrm{\πd}}^2{\cos }^2\widehat{\mathrm{\θ}}/\left(\mathrm{\λr}\right),D\left(\widehat{\mathrm{\θ}}\right)\stackrel{\mathrm{\Δ}}{=}{[\mathrm{\Φ}\left(\widehat{\mathrm{\θ}}\right),{\mathrm{\Φ}}^{*}\left(\widehat{\mathrm{\θ}}\right)J]}^T,$ 是块对角矩阵， $\mathrm{\Φ}\left(\widehat{\mathrm{\θ}}\right)=\mathrm{diag}\{e^{-j\widehat{\mathrm{\ω}}M},e^{-j\widehat{\mathrm{\ω}}(M-1)},...,1\},$ (·)^*表示取共轭， $\widehat{\mathrm{\ω}}\stackrel{\mathrm{\Δ}}{=}-2\mathrm{\π}d\sin \widehat{\mathrm{\θ}}/\mathrm{\λ},$ J是一个(M+1)×M的选择矩阵，定义为J＝[e_M,e_M-1,…,e₁]，e_m表示单位矩阵I_M+1的第m列，m＝1,2,…,M；

2.7)根据式(13)构造多项式 $p_r\left(z\right)\stackrel{\mathrm{\Δ}}{=}{z}^{M^2}{p}^H\left(z\right)D^H\left(\widehat{\mathrm{\θ}}\right){\mathrm{\Π}}_{\hat{Q}}D\left(\widehat{\mathrm{\θ}}\right)p\left(z\right),$ 其中 $p\left(z\right)\stackrel{\mathrm{\Δ}}{=}{[1,z,...,z^{M^2}]}^T,$ $z\stackrel{\mathrm{\Δ}}{=}{e}^{\mathrm{\π}{d}^2{\cos }^2\widehat{\mathrm{\θ}}/\left(\mathrm{\λr}\right)},$ 通过求多项式 $p_r\left(z\right)\stackrel{\mathrm{\Δ}}{=}{z}^{M^2}{p}^H\left(z\right)D^H\left(\widehat{\mathrm{\θ}}\right){\mathrm{\Π}}_{\hat{Q}}D\left(\widehat{\mathrm{\θ}}\right)p\left(z\right)$ 的K₂个零相位点来估计近场信号的距离；

2.8)更新斜投影算子返回步骤2.2)，继续顺序执行，如此往复循环，直至满足精度要求时结束。

5.如权利要求4所述一种远场和近场混合信号源的定位方法，其特征在于，所述步骤2.1)中斜投影算子的估计方法包括以下步骤：

a、构造一个(2M+1)×(2M+1)矩阵：

$R_o\stackrel{\mathrm{\Δ}}{=}\stackrel{\‾}{R}{\mathrm{\Π}}_{A_f}^{\⊥}---\left(14\right)$

其中A_f是远场信号方向矩阵，由估计出的远场信号的波达方向角计算得出；

b、对R_oΔ进行QR分解得斜投影算子由下述方式获得

其中是矩阵的后2M+1-K₂列，Δ是(2M+1)×(2M+1)的置换矩阵，不改变矩阵R_o的列相关性，表示广义逆。

6.如权利要求4所述一种远场和近场混合信号源的定位方法，其特征在于，所述步骤2.3)具体包括以下步骤：

仅包含近场信号信息的协方差矩阵R_n的反对角元素r_n(i)定义为

$r_n\left(i\right)\stackrel{\mathrm{\Δ}}{=}{\left(R_n\right)}_{i,2M+2-i}=\underset{k=K_1+1}{\overset{K}{\mathrm{\Σ}}}{r}_{\mathrm{sk}}^2{e}^{-j2(M+1-i){\mathrm{\ω}}_k}---\left(16\right)$

i＝1,2,…,2M+1，是第k个近场信号的功率；

7.如权利要求4所述一种远场和近场混合信号源的定位方法，其特征在于，更新斜投影算子的方法为：

$E_{A_f|A_n}=A_f{\left(A_f^H{\mathrm{\Π}}_{A_n}^{\⊥}{A}_f\right)}^{-1}{A}_f^H{\mathrm{\Π}}_{A_n}^{\⊥}---\left(25\right)$

A_f是远场信号方向矩阵，由估计出的远场信号的波达方向角计算得出。