CN104407335A - 一种3轴交叉阵列的doa估计方法 - Google Patents

Abstract

本发明提出了一种基于SLS-NC-ESPRIT算法的3轴交叉阵列的DOA估计方法，适用于严格二阶非圆(NC)信号的DOA的估计。不同于传统的DOA估计算法，本发明提出的方法首先利用了信号的非圆特性，扩大了阵列的虚拟空间，提高了估计精度并使得可检测信源数增加；其次利用了子阵列配置的重叠结构，考虑了信号子空间误差的重叠性，利用结构最小二乘(SLS)方法有效地解决了3个交叉阵列轴方向上的旋转不变方程；最后以增加约束条件的形式保证了3个轴方向具有近似相同的特征矢量，从而解决了当来波信号在某一轴方向有相同投影时引入的秩亏问题，保证了算法的有效性，提供了更精确的DOA估计值。

Description

一种3轴交叉阵列的DOA估计方法

技术领域

本发明涉及阵列信号处理技术领域，尤其涉及一种3轴交叉阵列的DOA估计方法。 

背景技术

来波信号DOA估计在雷达，声纳，移动通信和生物医学成像等多个领域获得了广泛的应用。众多用于解决高分辨率的DOA估计问题的技术已被提出，其中比较经典的技术有：极值搜索，多项式根，矩阵旋转不变等。作为一种矩阵旋转不变方法，ESPRIT算法能够有效地利用信号子空间的旋转不变特性计算出信号的DOA，并且可以给出解析解。通常，ESPRIT算法利用最小二乘(Least Squares,LS)或总体最小二乘(Total Least Squares,TLS)求解DOA估计值，由于LS和TLS方法没有考虑信号子空间重叠部分具有相同的误差这一特点，使得得到的DOA估计值误差较大。为了更好的利用模型具有的特殊结构信息，超定线性方程的结构化总体最小范数(Structured Total Least Norm,STLN)算法被提出，假设信号模型为Ax≈b，其中A和b中可能存在误差，STLN保留了A的仿射结构，例如，托普利茨、汉克尔等结构类型。由于ESPRIT算法的旋转不变方程没有STLN所需要的仿射结构，因此，STLN不适用于ESPRIT(Estimating Signal Parameters via Rotational Invariance Techniques)算法的求解。与STLN方法类似，SLS算法也忽略了残余矩阵的展开式的二次项，通过仿真实验可以发现，当涉及到重叠子阵列的配置时，SLS比LS和TLS获得的DOA估计值更加精确。ESPRIT算法最初是为确定一维的方位问题提出的，估计范围为[0,π]。为了解决二维的DOA估计问题，平面天线结构被利用，例如圆形、长方形和交叉型等阵列结构，这些阵列结构可以联合估计出信号的方位角和俯仰角。由于交叉阵列的特殊结构，相比于其他多维均匀阵列结构，它可以提供一个较大的孔径和更高的角度分辨率。 

为了进一步提高DOA估计器的性能，除了考虑信号子空间重叠部分具有相同的误差这一特点，还利用了信号在时域的特性，即NC特性。非圆信号如BPSK，OQPSK，PAM和ASK调制信号，已在许多现代通信系统中都有广泛的应用。另外，当I/Q失调时，圆信号也会变为非圆信号，所以非圆信号在实际应用中是非常普遍的，如何利用信号的非圆特性，提高 DOA估计器的性能具有较高的研究价值。通过利用所接收信号的非圆特性，一些改进的基于子空间的DOA估计已经被提出，如NC-MUSIC，NC-root-MUSIC，NC-ESPRIT和NC-unitary ESPRIT，NC-ESPRIT等算法，这类算法可以使阵列孔径加倍，提高DOA估计的精度。 

发明内容

为了克服上述问题，利用信号子空间重叠部分具有相同的误差扰动和到达信号的非圆特性，本发明提出了一种基于SLS-NC-ESPRIT算法的3轴交叉阵列DOA估计方法。由于旋转不变等式(Shift-Invariant Equations，SIEs)中重叠部分的信号子空间具有相同的误差扰动，SLS算法能够提供更精确的DOA估计。 

本发明通过如下技术方案实现： 

一种基于SLS-NC-ESPRIT算法的3轴交叉阵列的DOA估计方法，所述3轴交叉阵列由三个沿着x轴，y轴，z轴的线性子阵列以原点为几何中心组成，所有的天线阵元都是独立和等距的，3轴的阵元之间距离为δ，阵元数为M＝M_x+M_y+M_z，有d个远场窄带信号到达阵列，其中M_x,M_y,M_z分别表示在x，y，z轴方向的子阵列的包含的阵元个数；其特征在于：所述包括以下步骤： 

1)x(t)表示接收样本数据，利用信号非圆特性，得到扩张的测量数据向量 $x^{\left(\mathrm{nc}\right)}\left(t\right)=\left[\begin{array}{c}\mathrm{As}\left(t\right)\\ {\mathrm{\Π}}_M{A}^{*}{s}^{*}\left(t\right)\end{array}\right]+\left[\begin{array}{c}n\left(t\right)\\ {\mathrm{\Π}}_M{n}^{*}\left(t\right)\end{array}\right]=A^{\left(\mathrm{nc}\right)}s\left(t\right)+n^{\left(\mathrm{nc}\right)}\left(t\right),$ 其中Π_M为M×M的交换矩阵，该交换矩阵的反对角元素为1，其它元素为0，A为信号导向矢量矩阵，s(t)为信号矢量，n(t)为噪声； 

2)计算扩张的采样协方差矩阵其中，X^(nc)表示一个2M×N的接收数据矩阵，由N个采样快拍数据x^(nc)(t_n)，1≤n≤N组成； 

3)设 $\begin{array}{c}{J}_{\mathrm{\ξ}1}=[I_{M_{\mathrm{\ξ}}-1},0_{(M_{\mathrm{\ξ}}-1)\×1}]\\ J_{\mathrm{\ξ}2}=[0_{(M_{\mathrm{\ξ}}-1)\×1},I_{M_{\mathrm{\ξ}}-1}]\end{array},$ ξ∈{x,y,z}，此时3轴交叉阵列具有最大重叠度的选择矩阵如下： 

$K_{x1}^{\left(\mathrm{nc}\right)}=\left[\begin{array}{cccccc}{J}_{x1}& 0_{(M_x-1)\×M_y}& 0_{(M_x-1)\×M_z}& 0_{(M_x-1)\×M_z}& 0_{(M_x-1)\×M_y}& 0_{(M_x-1)\×M_x}\\ 0_{(M_x-1)\×M_x}& 0_{(M_x-1)\×M_y}& 0_{(M_x-1)\×M_z}& 0_{(M_x-1)\×M_z}& 0_{(M_x-1)\×M_y}& {\mathrm{\Π}}_{M_x-1}{J}_{x2}{\mathrm{\Π}}_{M_x}\end{array}\right]$

$K_{x2}^{\left(\mathrm{nc}\right)}=\left[\begin{array}{cccccc}{J}_{x2}& 0_{(M_x-1)\×M_y}& 0_{(M_x-1)\×M_z}& 0_{(M_x-1)\×M_z}& 0_{(M_x-1)\×M_y}& 0_{(M_x-1)\×M_x}\\ 0_{(M_x-1)\×M_x}& 0_{(M_x-1)\×M_y}& 0_{(M_x-1)\×M_z}& 0_{(M_x-1)\×M_z}& 0_{(M_x-1)\×M_y}& {\mathrm{\Π}}_{M_x-1}{J}_{x1}{\mathrm{\Π}}_{M_x}\end{array}\right]$

$K_{y1}^{\left(\mathrm{nc}\right)}=\left[\begin{array}{cccccc}{0}_{(M_y-1)\×M_x}& J_{y1}& 0_{(M_y-1)\×M_z}& 0_{(M_y-1)\×M_z}& 0_{(M_y-1)\×M_y}& 0_{(M_y-1)\×M_x}\\ 0_{(M_y-1)\×M_x}& 0_{(M_y-1)\×m_y}& 0_{(M_y-1)\×M_z}& 0_{(M_x-1)\×M_z}& {\mathrm{\Π}}_{M_y-1}{J}_{y2}{\mathrm{\Π}}_{M_y}& 0_{(M_y-1)\×M_z}\end{array}\right]$

$K_{y2}^{\left(\mathrm{nc}\right)}=\left[\begin{array}{cccccc}{0}_{(M_y-1)\×M_x}& J_{y2}& 0_{(M_y-1)\×M_z}& 0_{(M_y-1)\×M_z}& 0_{(M_y-1)\×M_y}& 0_{(M_y-1)\×M_x}\\ 0_{(M_y-1)\×M_x}& 0_{(M_y-1)\×m_y}& 0_{(M_y-1)\×M_z}& 0_{(M_x-1)\×M_z}& {\mathrm{\Π}}_{M_y-1}{J}_{y1}{\mathrm{\Π}}_{M_y}& 0_{(M_y-1)\×M_z}\end{array}\right]$

$K_{z1}^{\left(\mathrm{nc}\right)}=\left[\begin{array}{cccccc}{0}_{(M_z-1)\×M_x}& 0_{(M_z-1)\×M_y}& J_{z1}& 0_{(M_z-1)\×M_z}& 0_{(M_z-1)\×M_y}& 0_{(M_z-1)\×M_x}\\ 0_{(M_z-1)\×M_x}& 0_{(M_z-1)\×M_y}& 0_{(M_z-1)\×M_z}& {\mathrm{\Π}}_{M_z-1}{J}_{z2}{\mathrm{\Π}}_{M_z}& 0_{(M_z-1)\×M_y}& 0_{(M_z-1)\×M_x}\end{array}\right]$

$K_{z2}^{\left(\mathrm{nc}\right)}=\left[\begin{array}{cccccc}{0}_{(M_z-1)\×M_x}& 0_{(M_z-1)\×M_y}& J_{z2}& 0_{(M_z-1)\×M_z}& 0_{(M_z-1)\×M_y}& 0_{(M_z-1)\×M_x}\\ 0_{(M_z-1)\×M_x}& 0_{(M_z-1)\×M_y}& 0_{(M_z-1)\×M_z}& {\mathrm{\Π}}_{M_z-1}{J}_{z1}{\mathrm{\Π}}_{M_z}& 0_{(M_z-1)\×M_y}& 0_{(M_z-1)\×M_x}\end{array}\right];$

4)对扩张的采样协方差矩阵做SVD分解得到信号子空间估计 非圆信号的旋转不变等式可表示为： 

$K_{\mathrm{\ξ}1}^{\left(\mathrm{nc}\right)}{\hat{U}}_s^{\left(\mathrm{nc}\right)}{\mathrm{\γ}}_{\mathrm{\ξ}}=K_{\mathrm{\ξ}2}^{\left(\mathrm{nc}\right)}{\hat{U}}_s^{\left(\mathrm{nc}\right)},\mathrm{\ξ}\∈\{x,y,z\},$

对上式进行最小二乘法求解，得到ξ∈{x,y,z}作为γ_ξ的初始估计，信号空间估计的初始值为

5)改进的信号子空间估计可以表示为：得到残差矩阵 $R_{\mathrm{\ξ}}({\tilde{U}}_s^{\left(\mathrm{nc}\right)},{\mathrm{\γ}}_{\mathrm{\ξ}})=K_{\mathrm{\ξ}1}^{\left(\mathrm{nc}\right)}{\tilde{U}}_s^{\left(\mathrm{nc}\right)}{\mathrm{\γ}}_{\mathrm{\ξ}}-K_{\mathrm{\ξ}2}^{\left(\mathrm{nc}\right)}{\tilde{U}}_s^{\left(\mathrm{nc}\right)},\mathrm{\ξ}\∈\{x,y,z\};$

6)利用3轴分叉阵列的阵列矩阵γ_ξ,ξ∈{x,y,z}具有相同的特征向量的特点建立矩阵F₁,F₂和F₃： 

F₁＝γ_xγ_y-γ_yγ_x＝0_d×d

F₂＝γ_yγ_z-γ_zγ_y＝0_d×d

F₃＝γ_zγ_x-γ_xγ_z＝0_d×d； 

7)对于第k次迭代，残余矩阵,F₁,F₂和F₃可以表示为： 

$R({\tilde{U}}_{\mathrm{sk}}^{\left(\mathrm{nc}\right)},{\mathrm{\γ}}_{\mathrm{\ξk}})=K_{\mathrm{\ξ}1}^{\left(\mathrm{nc}\right)}{\tilde{U}}_{\mathrm{sk}}^{\left(\mathrm{nc}\right)}{\mathrm{\γ}}_{\mathrm{\ξ}}-K_{\mathrm{\ξ}2}^{\left(\mathrm{nc}\right)}{\tilde{U}}_{\mathrm{sk}}^{\left(\mathrm{nc}\right)}$

F_1k＝γ_xkγ_yk-γ_ykγ_xk

F_2k＝γ_ykγ_zk-γ_zkγ_yk

F_3k＝γ_zkγ_xk-γ_xkγ_zk； 

8)第k+1次的迭代中残余矩阵，F₁,F₂,F₃可写为： 

$\begin{array}{c}R({\tilde{U}}_{\mathrm{sk}}^{\left(\mathrm{nc}\right)},{\mathrm{\γ}}_{\mathrm{\ξk}})=R({\tilde{U}}_{\mathrm{sk}}^{\left(\mathrm{nc}\right)}+\mathrm{\Δ}{\hat{U}}_{\mathrm{sk}}^{\left(\mathrm{nc}\right)},{\mathrm{\γ}}_{\mathrm{\ξk}}+\mathrm{\Δ}{\mathrm{\γ}}_{\mathrm{\ξk}})\\ \≈R({\tilde{U}}_{\mathrm{sk}}^{\left(\mathrm{nc}\right)},{\mathrm{\γ}}_{\mathrm{\ξk}})+K_{\mathrm{\ξ}1}^{\left(\mathrm{nc}\right)}{\tilde{U}}_{\mathrm{sk}}^{\left(\mathrm{nc}\right)}\mathrm{\Δ}{\mathrm{\γ}}_{\mathrm{\ξk}}+K_{\mathrm{\ξ}1}^{\left(\mathrm{nc}\right)}\mathrm{\Δ}{\tilde{U}}_{\mathrm{sk}}^{\left(\mathrm{nc}\right)}{\mathrm{\γ}}_{\mathrm{\ξk}}-K_{\mathrm{\ξ}2}^{\left(\mathrm{nc}\right)}\mathrm{\Δ}{\tilde{U}}_{\mathrm{sk}}^{\left(\mathrm{nc}\right)}\end{array}$

F_1(k+1)＝(γ_xk+Δγ_xk)(γ_yk+Δγ_yk)-(γ_yk+Δγ_yk)(γ_xk+Δγ_xk) 

≈F_1k+γ_xkΔγ_yk+Δγ_xkγ_yk-γ_ykΔγ_xk-Δγ_ykγ_xk

F_2(k+1)＝(γ_yk+Δγ_yk)(γ_zk+Δγ_zk)-(γ_zk+Δγ_zk)(γ_yk+Δγ_yk) 

≈F_2k+γ_ykΔγ_zk+Δγ_ykγ_zk-γ_zkΔγ_yk-Δγ_zkγ_yk

F_3(k+1)＝(γ_zk+Δγ_zk)(γ_xk+Δγ_xk)-(γ_xk+Δγ_xk)(γ_zk+Δγ_zk) 

≈F_3k+γ_zkΔγ_xk+Δγ_zkγ_xk-γ_xkΔγ_zk-Δγ_xkγ_zk； 

9)建立代价函数： 

$\{{\widehat{\mathrm{\γ}}}_x,{\widehat{\mathrm{\γ}}}_y,{\widehat{\mathrm{\γ}}}_z\}=\underset{\mathrm{\Δ}{\mathrm{\γ}}_{\mathrm{\ξk}},\mathrm{\Δ}{\tilde{U}}_{\mathrm{sk}}^{\left(\mathrm{nc}\right)}}{\arg \min }\left|\right|{\left[\begin{array}{c}\mathrm{vec}\left\{R({\tilde{U}}_{\mathrm{sk}}^{\left(\mathrm{nc}\right)}+\mathrm{\Δ}{U}_{\mathrm{sk}}^{\left(\mathrm{nc}\right)},{\mathrm{\γ}}_{\mathrm{xk}})\right\}\\ \mathrm{vec}\left\{R({\tilde{U}}_{\mathrm{sk}}^{\left(\mathrm{nc}\right)}+\mathrm{\Δ}{U}_{\mathrm{sk}}^{\left(\mathrm{nc}\right)},{\mathrm{\γ}}_{\mathrm{yk}})\right\}\\ \mathrm{vec}\left\{R({\tilde{U}}_{\mathrm{sk}}^{\left(\mathrm{nc}\right)}+\mathrm{\Δ}{U}_{\mathrm{sk}}^{\left(\mathrm{nc}\right)},{\mathrm{\γ}}_{\mathrm{zk}})\right\}\\ \mathrm{vec}\left\{F_{1k}\right\}\\ \mathrm{vec}\left\{F_{2k}\right\}\\ \mathrm{vec}\left\{F_{3k}\right\}\\ \widetilde{\mathrm{\κ}}\·\mathrm{vec}\left\{\mathrm{\Δ}{\tilde{U}}_{\mathrm{sk}}^{\left(\mathrm{nc}\right)}\right\}\end{array}\right]+H_k\·\left[\begin{array}{c}\mathrm{vec}\left\{\mathrm{\Δ}{\mathrm{\γ}}_{\mathrm{xk}}\right\}\\ \mathrm{vec}\left\{\mathrm{\Δ}{\mathrm{\γ}}_{\mathrm{yk}}\right\}\\ \mathrm{vec}\left\{\mathrm{\Δ}{\mathrm{\γ}}_{\mathrm{zk}}\right\}\\ \mathrm{vec}\left\{\mathrm{\Δ}{\hat{U}}_{\mathrm{sk}}^{\left(\mathrm{nc}\right)}\right\}\end{array}\right]\left|\right|}_F,$

其中，为在第k步骤中估计的信号子空间的误差矩阵，γ_ξk＝γ_ξk-1+Δγ_ξk-1,ξ∈{x,y,z}，vec{·}是向量化函数，

$H_k=\left[\begin{array}{cccc}{I}_d\⊗\left(K_{x1}^{\left(\mathrm{nc}\right)}{\tilde{U}}_{\mathrm{sk}}^{\left(\mathrm{nc}\right)}\right)& 0_{2(M_y-1)d\×d^2}& 0_{2(M_z-1)d\×d^2}& {\mathrm{\γ}}_{\mathrm{xk}}^T\⊗K_{x1}^{\left(\mathrm{nc}\right)}-I_d\⊗K_{x2}^{\left(\mathrm{nc}\right)}\\ 0_{2(M_x-1)d\×d^2}& I_d\⊗\left(K_{y1}^{\left(\mathrm{nc}\right)}{\tilde{U}}_{\mathrm{sk}}^{\left(\mathrm{nc}\right)}\right)& 0_{2(M_z-1)d\×d^2}& {\mathrm{\γ}}_{\mathrm{yk}}^T\⊗K_{y1}^{\left(\mathrm{nc}\right)}-I_d\⊗K_{y2}^{\left(\mathrm{nc}\right)}\\ 0_{2(M_x-1)d\×d^2}& 0_{2(M_y-1)d\×d^2}& I_d\⊗\left(K_{z1}^{\left(\mathrm{nc}\right)}{\tilde{U}}_{\mathrm{sk}}^{\left(\mathrm{nc}\right)}\right)& {\mathrm{\γ}}_{\mathrm{zk}}^T\⊗K_{z1}^{\left(\mathrm{nc}\right)}-I_d\⊗K_{z2}^{\left(\mathrm{nc}\right)}\\ {\mathrm{\γ}}_{\mathrm{yk}}^T\⊗I_d-I_d\⊗{\mathrm{\γ}}_{\mathrm{yk}}& I_d\⊗{\mathrm{\γ}}_{\mathrm{yk}}-{\mathrm{\γ}}_{\mathrm{yk}}^T\⊗I_d& 0_{d^2\×d^2}& 0_{d^2\×2\mathrm{Md}}\\ 0_{d^2\×d^2}& {\mathrm{\γ}}_{\mathrm{zk}}^T\⊗I_d-I_d\⊗{\mathrm{\γ}}_{\mathrm{zk}}& I_d\⊗{\mathrm{\γ}}_{\mathrm{yk}}-{\mathrm{\γ}}_{\mathrm{yk}}^T\⊗I_d& 0_{d^2\×2\mathrm{Md}}\\ I_d\⊗{\mathrm{\γ}}_{\mathrm{yk}}-{\mathrm{\γ}}_{\mathrm{yk}}^T\⊗I_d& 0_{d^2\×d^2}& {\mathrm{\γ}}_{\mathrm{xk}}^T\⊗I_d-I_d\⊗{\mathrm{\γ}}_{\mathrm{xk}}& 0_{d^2\×2\mathrm{Md}}\\ 0_{2\mathrm{Md}\×d^2}& 0_{2\mathrm{Md}\×d^2}& 0_{2\mathrm{Md}\×d^2}& \widetilde{\mathrm{\κ}}{I}_{2\mathrm{Md}}\end{array}\right];$

10)判断代价函数是否满足ξ∈{x,y,z}，如果满足则认为收敛到最优解，停止迭代，其中ε为允许的最大误差范围；如果不满足，则令k＝k+1，返回步骤9； 

11)将得到的最优解做EVD分解，得到μ_xp,μ_yp,μ_zp,p＝1,...,d的最优估计值p＝1,...,d，其中， 

${\mathrm{\μ}}_x({\mathrm{\θ}}_i,{\mathrm{\φ}}_i)=\frac{2\mathrm{\π}}{{\mathrm{\λ}}_0}\mathrm{\δ}\cos \left({\mathrm{\φ}}_i\right)\sin \left({\mathrm{\θ}}_i\right)$

${\mathrm{\μ}}_y({\mathrm{\θ}}_i,{\mathrm{\φ}}_i)=\frac{2\mathrm{\π}}{{\mathrm{\λ}}_0}\mathrm{\δ}\sin \left({\mathrm{\φ}}_i\right)\sin \left({\mathrm{\θ}}_i\right)$

${\mathrm{\μ}}_z({\mathrm{\θ}}_i,{\mathrm{\φ}}_i)=\frac{2\mathrm{\π}}{{\mathrm{\λ}}_0}\mathrm{\δ}\cos \left({\mathrm{\φ}}_i\right),$

λ₀为波长，θ_i,φ_i分别表示第i个信号的俯仰角和方位角； 

根据d个θ_i,φ_i得到d个信源的DOA的估计值 

本发明的有益效果是：对于3轴交叉阵列，当两个或多个到达信号在x轴，y轴，z轴具有相同的投影时，会出现秩亏问题，此时LS，TLS和SLS方法将失效。为了解决秩亏问题，本发明以增加约束条件的形式保证了3个轴方向具有近似相同的特征矢量，从而保证了算法的有效性，可以提供更精确的DOA估计值。因为前者有效地利用了信号的非圆特性，本发明提出的基于SLS-NC-ESPRIT算法的3轴交叉阵列DOA估计方法能够提供比SLS-ESPRIT算法更好的性能。 

附图说明

图1是3轴交叉阵列的结构示意图； 

图2是本发明提供的基于SLS-NC-E ESPRIT算法的3轴交叉阵列DOA估计方法流程图； 

图3是本发明的方法与传统算法的空间频率估计值的RMSE与信噪比的关系曲线图(θ₁＝θ₂,φ₁≠φ₂,N＝300)； 

图4是本发明的方法与传统算法的空间频率估计值的RMSE与信噪比的关系曲线图(θ₁＝θ₂,φ₁≠φ₂,N＝20)； 

图5是本发明的方法与传统算法的空间频率估计值的RMSE与快拍数的关系曲线图(θ₁＝θ₂,φ₁≠φ₂)； 

图6是本发明的方法与传统算法的空间频率估计值的RMSE与信噪比的关系曲线图(θ₁≠θ₂,φ₁≠φ₂,N＝300)； 

图7是本发明的方法与传统算法的空间频率估计值的RMSE与信噪比的关系曲线图(θ₁≠θ₂,φ₁≠φ₂,N＝20)； 

图8是本发明的方法与传统算法的空间频率估计值的RMSE与快拍数的关系曲线图(θ₁≠θ₂,φ₁≠φ₂)； 

图9是本发明的方法与传统算法的空间频率估计值的RMSE与方位间 隔φ的关系曲线图； 

图10是本发明的方法与传统算法的空间频率估计值的RMSE与仰角间隔θ的关系曲线图。 

具体实施方式

下面结合附图说明及具体实施方式对本发明进一步说明。 

3轴交叉阵列结构如附图1所示，由三个沿着x轴，y轴，z轴的线性子阵列(以原点为几何中心)组成。假设所有的天线阵元都是独立和等距的，3轴的阵元之间距离δ是相同的。假设3轴交叉阵列的阵元数为M，有d个远场窄带信号到达阵列，这里M＝M_x+M_y+M_z，其中M_x,M_y,M_z分别表示在x，y，z轴方向的子阵列的包含的阵元个数。将阵列接收信号建模为： 

x(t)＝As(t)+n(t)    (1) 

其中s(t)＝[s₁(t),...,s_d(t)]^T为信号矢量，n(t)＝[n₁(t),...,n_M(t)]^T包含加性的传感器噪声；为信号导向矢量矩阵，包含d个信号导向矢量 i＝1,...,d，其中θ_i,φ_i分别表示第i个信号的俯仰角和方位角；信源数d一般假定是预先已知的或者可以利用信息理论准则对其进行估计。3轴交叉阵列的导向矢量如下： 

其中， 

矢量p_n,n＝1,...,M表示天线位置，μ_i可以表示为 

其中： 

$\begin{array}{c}{\mathrm{\μ}}_x({\mathrm{\θ}}_i,{\mathrm{\φ}}_i)=\frac{2\mathrm{\π}}{{\mathrm{\λ}}_0}\mathrm{\δ}\cos \left({\mathrm{\φ}}_i\right)\sin \left({\mathrm{\θ}}_i\right)\\ {\mathrm{\μ}}_y({\mathrm{\θ}}_i,{\mathrm{\φ}}_i)=\frac{2\mathrm{\π}}{{\mathrm{\λ}}_0}\mathrm{\δ}\sin \left({\mathrm{\φ}}_i\right)\sin \left({\mathrm{\θ}}_i\right)\\ {\mathrm{\μ}}_z({\mathrm{\θ}}_i,{\mathrm{\φ}}_i)=\frac{2\mathrm{\π}}{{\mathrm{\λ}}_0}\mathrm{\δ}\cos \left({\mathrm{\φ}}_i\right)\end{array}---\left(5\right)$

这里λ₀为波长，假设噪声为高斯白噪声过程，其均值为零，协方差矩阵为 与信号s(t)是不相关的。 

假设信号是严格非圆信号，也就是说，在星座图中，信号的所有可能状态可以表示在一条线上，这样的信号在实际应用中广泛存在，例如，BPSK，OQPSK和PAM等信号。将信号矢量分解为: 

s(t)＝ψs₀(t)    (6) 

其中φ_i，i＝1,…,d为源信号的初始相位，s₀(t)是实值的信号矢量，信号s_i(t)的非圆系数定义为： 

${\mathrm{\γ}}_{\mathrm{sl}}=\frac{E\left\{s_l^2\left(t\right)\right\}}{E\left\{{\left|s_l\left(t\right)\right|}^2\right\}},l=1,...,d---\left(7\right)$

根据前述假设，接收信号矢量x(t)的协方差矩阵可以表示为： 

$R_x=E\left\{x\left(t\right)x{\left(t\right)}^H\right\}=A{R}_s{A}^H+{\mathrm{\σ}}_n^2{I}_M---\left(8\right)$

其中，是每个阵列天线的接收的噪声功率，R_s＝E{s(t)s(t)^H}为信号协方差矩阵，对R_x进行SVD分解可得： 

R_x＝UΣU^H    (9) 

其中Σ＝diag{λ₁,...,λ_M}，λ₁,…,λ_M为接收信号协方差矩阵R_x的特征值，并且满足U＝[u₁,...,u_M],i＝1,...,M为对应的特征向量。有 其中Σ_s＝diag{λ₁,...,λ_d}为主特征值，U_s＝[u₁,...,u_d]为信号子空间，U_n＝[u_d+1,...,u_M]为噪声子空间，根据式(7)，可得： 

$A{R}_s{A}^H+{\mathrm{\σ}}_n^2{I}_M=U_s{\mathrm{\Σ}}_s{U}_s^H+{\mathrm{\σ}}_n^2{U}_n{U}_n^H---\left(10\right)$

将利用到(10)中，得到： 

$U_s=A{R}_s{A}^H{U}_s{[{\mathrm{\Σ}}_s-{\mathrm{\σ}}_n^2{I}_M]}^{-1}=\mathrm{AT}---\left(11\right)$

其中，是非奇异矩阵，U_s为信号子空间。式(11)是ESPRIT算法的能够给出高分辨率DOA估计值的关键所在。 

假设ESPRIT算法使用具有最大重叠的子阵列，此时选择矩阵为： 

$\begin{array}{c}{J}_{\mathrm{\ξ}1}=[I_{M_{\mathrm{\ξ}}-1},0_{(M_{\mathrm{\ξ}}-1)\×1}]\\ J_{\mathrm{\ξ}2}=[0_{(M_{\mathrm{\ξ}}-1)\×1},I_{M_{\mathrm{\ξ}}-1}]\end{array},\mathrm{\ξ}\∈\{x,y,z\}---\left(12\right)$

此时3轴交叉阵列的选择矩阵如下： 

$\begin{array}{c}{K}_{\mathrm{xl}}=[J_{\mathrm{xl}},0_{(M_x-1)\×M_y},0_{(M_x-1)\×M_z}],l=\mathrm{1,2}\\ K_{\mathrm{yl}}=[0_{(M_y-1)\×M_x},J_{\mathrm{yl}},0_{(M_y-1)\×M_z}],l=\mathrm{1,2}\\ K_{\mathrm{zl}}=[0_{(M_z-1)\×M_x},0_{(M_z-1)\×M_y},J_{\mathrm{zl}}],l=\mathrm{1,2}\end{array}---\left(13\right)$

由于所有子阵列具有旋转不变的特性，可得到 

K_ξ1AΦ_ξ＝K_ξ2A,ξ∈{x,y,z}    (14) 

其中根据式(11)可得: 

K_ξ1U_sγ_ξ＝K_ξ2U_s    (15) 

其中 

γ_ξ＝T^-1Φ_ξT    (16) 

为了利用信号的非圆信息，应用广泛的线性处理技术，定义扩张的测量数据向量x^(nc)(t)为： 

其中Π_M为M×M的交换矩阵，其反对角元素为1，其它元素为0。将式(1)和(6)代入(17)，得： 

$\begin{array}{c}{x}^{\left(\mathrm{nc}\right)}\left(t\right)=\left[\begin{array}{c}\mathrm{As}\left(t\right)\\ {\mathrm{\Π}}_M{A}^{*}{s}^{*}\left(t\right)\end{array}\right]+\left[\begin{array}{c}n\left(t\right)\\ {\mathrm{\Π}}_M{n}^{*}\left(t\right)\end{array}\right]\\ =A^{\left(\mathrm{nc}\right)}s\left(t\right)+n^{\left(\mathrm{nc}\right)}\left(t\right)\end{array}---\left(18\right)$

其中A^(nc)可以视为扩展的阵列导向矩阵，与原导向矢量相比，它使阵列虚拟孔径变为原来的两倍，n^(nc)(t)代表扩展的噪声矢量。应当指出的是，阵列虚拟孔径加倍能够提高DOA估计的分辨率，并使得最大可检测信源数加倍。 

设X^(nc)表示一个2M×N的接收数据矩阵，由N个采样快拍数据x^(nc)(t_n)，1≤n≤N组成。类似于ESPRIT算法，这里选择矩阵的结构与扩张采样矩阵X^(nc)的定义有关系，并且这些矩阵式(19)在本发明中是首次提出。 

$\begin{array}{c}{K}_{x1}^{\left(\mathrm{nc}\right)}=\left[\begin{array}{cccccc}{J}_{x1}& 0_{(M_x-1)\×M_y}& 0_{(M_x-1)\×M_z}& 0_{(M_x-1)\×M_z}& 0_{(M_x-1)\×M_y}& 0_{(M_x-1)\×M_x}\\ 0_{(M_x-1)\×M_x}& 0_{(M_x-1)\×M_y}& 0_{(M_x-1)\×M_z}& 0_{(M_x-1)\×M_z}& 0_{(M_x-1)\×M_y}& {\mathrm{\Π}}_{M_x-1}{J}_{x2}{\mathrm{\Π}}_{M_x}\end{array}\right]\\ K_{x2}^{\left(\mathrm{nc}\right)}=\left[\begin{array}{cccccc}{J}_{x2}& 0_{(M_x-1)\×M_y}& 0_{(M_x-1)\×M_z}& 0_{(M_x-1)\×M_z}& 0_{(M_x-1)\×M_y}& 0_{(M_x-1)\×M_x}\\ 0_{(M_x-1)\×M_x}& 0_{(M_x-1)\×M_y}& 0_{(M_x-1)\×M_z}& 0_{(M_x-1)\×M_z}& 0_{(M_x-1)\×M_y}& {\mathrm{\Π}}_{M_x-1}{J}_{x1}{\mathrm{\Π}}_{M_x}\end{array}\right]\\ K_{y1}^{\left(\mathrm{nc}\right)}=\left[\begin{array}{cccccc}{0}_{(M_y-1)\×M_x}& J_{y1}& 0_{(M_y-1)\×M_z}& 0_{(M_y-1)\×M_z}& 0_{(M_y-1)\×M_y}& 0_{(M_y-1)\×M_x}\\ 0_{(M_y-1)\×M_x}& 0_{(M_y-1)\×M_y}& 0_{(M_y-1)\×M_z}& 0_{(M_x-1)\×M_z}& {\mathrm{\Π}}_{M_y-1}{J}_{y1}{\mathrm{\Π}}_{M_y}& 0_{(M_y-1)\×M_z}\end{array}\right]\\ K_{y2}^{\left(\mathrm{nc}\right)}=\left[\begin{array}{cccccc}{0}_{(M_y-1)\×M_x}& J_{y2}& 0_{(M_y-1)\×M_z}& 0_{(M_y-1)\×M_z}& 0_{(M_y-1)\×M_y}& 0_{(M_y-1)\×M_x}\\ 0_{(M_y-1)\×M_z}& 0_{(M_y-1)\×M_y}& J_{z1}& 0_{(M_z-1)\×M_z}& {\mathrm{\Π}}_{M_y-1}{J}_{y1}{\mathrm{\Π}}_{M_y}& 0_{(M_y-1)\×M_z}\end{array}\right]\\ K_{z1}^{\left(\mathrm{nc}\right)}=\left[\begin{array}{cccccc}{0}_{(M_z-1)\×M_x}& 0_{(M_z-1)\×M_y}& J_{z1}& 0_{(M_z-1)\×M_z}& 0_{(M_z-1)\×M_y}& 0_{(M_z-1)\×M_x}\\ 0_{(M_z-1)\×M_x}& 0_{(M_z-1)\×M_y}& 0_{(M_z-1)\×M_z}& {\mathrm{\Π}}_{M_z-1}{J}_{z2}{\mathrm{\Π}}_{M_z}{}_{}& 0_{(M_z-1)\×M_y}& 0_{(M_z-1)\×M_x}\end{array}\right]\\ K_{z2}^{\left(\mathrm{nc}\right)}=\left[\begin{array}{cccccc}{0}_{(M_z-1)\×M_x}& 0_{(M_z-1)\×M_y}& J_2& 0_{(M_z-1)\×M_z}& 0_{(M_z-1)\×M_y}& 0_{(M_z-1)\×M_x}\\ 0_{(M_z-1)\×M_x}& 0_{(M_z-1)\×M_y}& 0_{(M_z-1)\×M_x}& {\mathrm{\Π}}_{M_z-1}{J}_{z1}{\mathrm{\Π}}_{M_z}& 0_{(M_z-1)\×M_y}& 0_{(M_z-1)\×M_x}\end{array}\right]\end{array}---\left(19\right)$

由于3轴交叉阵列的每个子阵列都是旋转不变的，因此同样A^(nc)具有旋转不变特性，故满足下式： 

$K_{\mathrm{\ξ}1}^{\left(\mathrm{nc}\right)}{A}^{\left(\mathrm{nc}\right)}{\mathrm{\Φ}}_{\mathrm{\ξ}}=K_{\mathrm{\ξ}2}^{\left(\mathrm{nc}\right)}{A}^{\left(\mathrm{nc}\right)},\mathrm{\ξ}\∈\{x,y,z\}---\left(20\right)$

根据式(17)，通过对接收数据矩阵X^(nc)进行SVD分解，提取d个最大的奇异值对应的特征矢量，组成信号子空间的估计此时A^(nc)和具有近似相同的列空间，即其中T^(nc)∈C^d×d是非奇异矩阵。与ESPRIT相类似，非圆信号的SIEs可表示为： 

$K_{\mathrm{\ξ}1}^{\left(\mathrm{nc}\right)}{\hat{U}}_s^{\left(\mathrm{nc}\right)}{\mathrm{\γ}}_{\mathrm{\ξ}}=K_{\mathrm{\ξ}2}^{\left(\mathrm{nc}\right)}{\hat{U}}_s^{\left(\mathrm{nc}\right)},\mathrm{\ξ}\∈\{x,y,z\}---\left(21\right)$

求解式(21)最常见的方法是基于LS的方法，利用LS方法计算可得γ_ξ的估计值为ξ∈{x,y,z}，然而，LS算法假设是已知的且无误差，误差只存在于中，这显然与实际应用中的情况是不相符的。与LS的算法不同，TLS方法假设和都存在误差，通过最小化存在于和的干扰和从而获得比LS更精确的DOA估计值。但是，只有在两个子阵列的阵元没有重叠的情况下，TLS可以提供最优的DOA估计值。反之，当沿x，y和z轴方向的每个线性阵列的两个子阵列有重叠的阵元时，SLS方法由于考虑了子阵列间的在结构上存在的具体关系，因此相对于TLS方法，可以提高估计值的精确度。所以本发明讨论了使用SLS算法来求解非圆信号的SIEs，仿真结果证明，在子阵列之间重叠度最大的情况下，SLS算法大大提高了俯仰角和方位角的估计准确度。 

因为只张成了包含噪声的信号子空间，假设在信号子空间中存在一个小的扰动改进的信号子空间估计可以表示为 定义残差矩阵为： 

$R_{\mathrm{\ξ}}({\tilde{U}}_s^{\left(\mathrm{nc}\right)},{\mathrm{\γ}}_{\mathrm{\ξ}})=K_{\mathrm{\ξ}1}^{\left(\mathrm{nc}\right)}{\tilde{U}}_s^{\left(\mathrm{nc}\right)}{\mathrm{\γ}}_{\mathrm{\ξ}}-K_{\mathrm{\ξ}2}^{\left(\mathrm{nc}\right)}{\tilde{U}}_s^{\left(\mathrm{nc}\right)},\mathrm{\ξ}\∈\{x,y,z\}---\left(22\right)$

SLS算法在使得的Frobenius范数最小化的基础上，同时保持的Frobenius范数尽可能小，综上所述，建立如下代价函数，通过令下面的表达式最小化，进一步提高DOA估计的精度： 

$\underset{\mathrm{\Δ}{\hat{U}}_s^{\left(\mathrm{nc}\right)},\mathrm{\Δ}{\hat{U}}_s^{\left(\mathrm{nc}\right)}}{\min }{\left|\right|\left[\begin{array}{c}{R}_{\mathrm{\ξ}}({\tilde{U}}_s^{\left(\mathrm{nc}\right)},{\mathrm{\γ}}_{\mathrm{\ξ}})\\ \mathrm{\κ}\·\mathrm{\Δ}{\hat{U}}_s^{\left(\mathrm{nc}\right)}\end{array}\right]\left|\right|}_F,\mathrm{\ξ}\∈\{x,y,z\}---\left(23\right)$

其中是调整和残差矩阵的之间的关系的调节因子。如果α的值大于1，在迭代过程中的的每个元素扰动范围相对于中的元素要大，但是已有文献验证SLS算法对于α的选择并不敏感。定义和中的第(i,j)个元素分别为r_ij和e_ij，式(23)可写为： 

$\underset{\mathrm{\Δ}{\hat{U}}_s^{\left(\mathrm{nc}\right)},{\mathrm{\γ}}_{\mathrm{\ξ}}}{\min }\underset{j=1}{\overset{d}{\mathrm{\Σ}}}(\frac{\mathrm{\α}}{2(M-3)d}\underset{i=1}{\overset{2(M-3)}{\mathrm{\Σ}}}{\left|r_{\mathrm{ij}}\right|}^2+\frac{1}{2\mathrm{Md}}\underset{i=1}{\overset{2M}{\mathrm{\Σ}}}{\left|e_{\mathrm{ij}}\right|}^2)---\left(24\right)$

因为3轴交叉阵列的每个轴上的子阵列都是线性的，当两个或多个到达的信号在轴ξ，ξ∈{x,y,z}上有相同的投影时，矩阵l＝1,2可能出现秩亏或者近似秩亏的情况，因此，LS，TLS和SLS算法在特定的ξ轴上将无法提供正确的空间频率估计值。但是，因为矩阵γ_ξ,ξ∈{x,y,z}具有相同的特征向量组成矩阵T，因而只存在一个有效的可行解。根据矩阵分析的知识可知，矩阵A和B有相同的特征向量集的充分必要条件是AB＝BA。利用这个特性，本发明的方法能够通过运用矩阵γ_ξ,ξ∈{x,y,z}间的关系进一步提高算法的性能，定义矩阵F₁,F₂和F₃，建立如下等式： 

F₁＝γ_xγ_y-γ_yγ_x＝0_d×d

F₂＝γ_yγ_z-γ_zγ_y＝0_d×d    (25) 

F₃＝γ_zγ_x-γ_xγ_z＝0_d×d

式(25)保证了矩阵γ_ξ,ξ∈{x,y,z}有相同的特征向量集，可以解决当两个或多个到达的信号在轴ξ，ξ∈{x,y,z}上有相同的投影时导致的秩亏问题，因而估计的DOA的精度得到了提升。 

通过对式(23)进行迭代求解，同时最小化F₁,F₂,F₃的Frobenius范数就可以获得本发明所采用的改进的SLS-NC-ESPRIT算法来求解基于3轴交叉阵列的俯仰角和方位角估计的问题。 

在第k次的迭代中，残差矩阵F₁,F₂和F₃可以表示为： 

$\begin{array}{c}R({\tilde{U}}_{\mathrm{sk}}^{\left(\mathrm{nc}\right)},{\mathrm{\γ}}_{\mathrm{\ξk}})=K_{\mathrm{\ξ}1}^{\left(\mathrm{nc}\right)}{\tilde{U}}_{\mathrm{sk}}^{\left(\mathrm{nc}\right)}{\mathrm{\γ}}_{\mathrm{\ξ}}-K_{\mathrm{\ξ}2}^{\left(\mathrm{nc}\right)}{\tilde{U}}_{\mathrm{sk}}^{\left(\mathrm{nc}\right)}\\ F_{1k}={\mathrm{\γ}}_{\mathrm{xk}}{\mathrm{\γ}}_{\mathrm{yk}}-{\mathrm{\γ}}_{\mathrm{yk}}{\mathrm{\γ}}_{\mathrm{xk}}\\ F_{2k}={\mathrm{\γ}}_{\mathrm{yk}}{\mathrm{\γ}}_{\mathrm{zk}}-{\mathrm{\γ}}_{\mathrm{zk}}{\mathrm{\γ}}_{\mathrm{yk}}\\ F_{3k}={\mathrm{\γ}}_{\mathrm{zk}}{\mathrm{\γ}}_{\mathrm{xk}}-{\mathrm{\γ}}_{\mathrm{xk}}{\mathrm{\γ}}_{\mathrm{zk}}\end{array}---\left(26\right)$

其中ξ∈{x,y,z}.因此，在第k+1次迭代中残差矩阵F₁,F₂和F₃为： 

$\begin{array}{c}R({\tilde{U}}_{\mathrm{sk}}^{\left(\mathrm{nc}\right)},{\mathrm{\γ}}_{\mathrm{\ξk}})=R({\tilde{U}}_{\mathrm{sk}}^{\left(\mathrm{nc}\right)}+\mathrm{\Δ}{\hat{U}}_{\mathrm{sk}}^{\left(\mathrm{nc}\right)},{\mathrm{\γ}}_{\mathrm{\ξk}}+\mathrm{\Δ}{\mathrm{\γ}}_{\mathrm{\ξk}})\\ \≈R({\tilde{U}}_{\mathrm{sk}}^{\left(\mathrm{nc}\right)},{\mathrm{\γ}}_{\mathrm{\ξk}})+K_{\mathrm{\ξ}1}^{\left(\mathrm{nc}\right)}{\tilde{U}}_{\mathrm{sk}}^{\left(\mathrm{nc}\right)}\mathrm{\Δ}{\mathrm{\γ}}_{\mathrm{\ξk}}+K_{\mathrm{\ξ}1}^{\left(\mathrm{nc}\right)}\mathrm{\Δ}{\tilde{U}}_{\mathrm{sk}}^{\left(\mathrm{nc}\right)}{\mathrm{\γ}}_{\mathrm{\ξk}}-K_{\mathrm{\ξ}2}^{\left(\mathrm{nc}\right)}\mathrm{\Δ}{\tilde{U}}_{\mathrm{sk}}^{\left(\mathrm{nc}\right)}\\ F_{1(k+1)}=({\mathrm{\γ}}_{\mathrm{xk}}+\mathrm{\Δ}{\mathrm{\γ}}_{\mathrm{xk}})({\mathrm{\γ}}_{\mathrm{yk}}+\mathrm{\Δ}{\mathrm{\γ}}_{\mathrm{yk}})-({\mathrm{\γ}}_{\mathrm{yk}}+\mathrm{\Δ}{\mathrm{\γ}}_{\mathrm{yk}})({\mathrm{\γ}}_{\mathrm{xk}}+\mathrm{\Δ}{\mathrm{\γ}}_{\mathrm{xk}})\\ \≈F_{1k}+{\mathrm{\γ}}_{\mathrm{xk}}\mathrm{\Δ}{\mathrm{\γ}}_{\mathrm{yk}}+\mathrm{\Δ}{\mathrm{\γ}}_{\mathrm{xk}}{\mathrm{\γ}}_{\mathrm{yk}}-{\mathrm{\γ}}_{\mathrm{yk}}\mathrm{\Δ}{\mathrm{\γ}}_{\mathrm{xk}}-\mathrm{\Δ}{\mathrm{\γ}}_{\mathrm{yk}}{\mathrm{\γ}}_{\mathrm{xk}}\\ F_{2(k+1)}=({\mathrm{\γ}}_{\mathrm{yk}}+\mathrm{\Δ}{\mathrm{\γ}}_{\mathrm{yk}})({\mathrm{\γ}}_{\mathrm{zk}}+\mathrm{\Δ}{\mathrm{\γ}}_{\mathrm{zk}})-({\mathrm{\γ}}_{\mathrm{zk}}+\mathrm{\Δ}{\mathrm{\γ}}_{\mathrm{zk}})({\mathrm{\γ}}_{\mathrm{yk}}+\mathrm{\Δ}{\mathrm{\γ}}_{\mathrm{yk}})\\ \≈F_{2k}+{\mathrm{\γ}}_{\mathrm{yk}}\mathrm{\Δ}{\mathrm{\γ}}_{\mathrm{zk}}+\mathrm{\Δ}{\mathrm{\γ}}_{\mathrm{yk}}{\mathrm{\γ}}_{\mathrm{zk}}-{\mathrm{\γ}}_{\mathrm{zk}}\mathrm{\Δ}{\mathrm{\γ}}_{\mathrm{yk}}-\mathrm{\Δ}{\mathrm{\γ}}_{\mathrm{zk}}{\mathrm{\γ}}_{\mathrm{yk}}\\ F_{3(k+1)}=({\mathrm{\γ}}_{\mathrm{zk}}+\mathrm{\Δ}{\mathrm{\γ}}_{\mathrm{zk}})({\mathrm{\γ}}_{\mathrm{xk}}+\mathrm{\Δ}{\mathrm{\γ}}_{\mathrm{xk}})-({\mathrm{\γ}}_{\mathrm{xk}}+\mathrm{\Δ}{\mathrm{\γ}}_{\mathrm{xk}})({\mathrm{\γ}}_{\mathrm{zk}}+\mathrm{\Δ}{\mathrm{\γ}}_{\mathrm{zk}})\\ \≈F_{3k}+{\mathrm{\γ}}_{\mathrm{zk}}\mathrm{\Δ}{\mathrm{\γ}}_{\mathrm{xk}}+\mathrm{\Δ}{\mathrm{\γ}}_{\mathrm{zk}}{\mathrm{\γ}}_{\mathrm{xk}}-{\mathrm{\γ}}_{\mathrm{xk}}\mathrm{\Δ}{\mathrm{\γ}}_{\mathrm{zk}}-\mathrm{\Δ}{\mathrm{\γ}}_{\mathrm{xk}}{\mathrm{\γ}}_{\mathrm{zk}}\end{array}---\left(27\right)$

注意(27)中最后的近似式是忽略二次项得到的。根据矩阵分析可知，对于已知矩阵有： 

$\mathrm{vec}\left\{C_1{C}_2{C}_3\right\}=(C_3^T\⊗C_1)\mathrm{vec}\left\{C_2\right\}---\left(28\right)$

将矢量化运算运用到式(27)中，得： 

$\begin{array}{c}\mathrm{vec}\{R({\tilde{U}}_{s(k+1)}^{\left(\mathrm{nc}\right)},{\mathrm{\γ}}_{\mathrm{\ξ}(k+1)}\}\≈\mathrm{vec}\left\{R({\tilde{U}}_{\mathrm{sk}}^{\left(\mathrm{nc}\right)},{\mathrm{\γ}}_{\mathrm{\ξk}})\right\}\\ +[I_d\⊗\left(K_{\mathrm{\ξ}1}^{\left(\mathrm{nc}\right){\tilde{U}}_{\mathrm{sk}}^{\left(\mathrm{nc}\right)}}\right)]\mathrm{vec}\left\{\mathrm{\Δ}{\mathrm{\γ}}_{\mathrm{\ξk}}\right\}]\\ +[{\mathrm{\γ}}_{\mathrm{\ξk}}^T\⊗K_{\mathrm{\ξ}1}^{\left(\mathrm{nc}\right)}-I_d\⊗K_{\mathrm{\ξ}1}^{\left(\mathrm{nc}\right)}]\mathrm{vec}\left\{\mathrm{\Δ}{\hat{U}}_{\mathrm{sk}}^{\left(\mathrm{nc}\right)}\right\}\\ \mathrm{vec}\left\{F_{1(k+1)}\right\}\≈\mathrm{vec}\left\{F_{1k}\right\}\\ +[{\mathrm{\γ}}_{\mathrm{yk}}^T\⊗I_d-I_d\⊗{\mathrm{\γ}}_{\mathrm{yk}}]\mathrm{vec}\left\{\mathrm{\Δ}{\mathrm{\γ}}_{\mathrm{xk}}\right\}\\ +[I_d\⊗{\mathrm{\γ}}_{\mathrm{xk}}-{\mathrm{\γ}}_{\mathrm{xk}}\⊗I_d]\mathrm{vec}\left\{\mathrm{\Δ}{\mathrm{\γ}}_{\mathrm{yk}}\right\}\\ \mathrm{vec}\left\{F_{2(k+1)}\right\}\≈\mathrm{vec}\left\{F_{2k}\right\}\\ +[{\mathrm{\γ}}_{\mathrm{zk}}^T\⊗I_d-I_d\⊗{\mathrm{\γ}}_{\mathrm{zk}}]\mathrm{vec}\left\{\mathrm{\Δ}{\mathrm{\γ}}_{\mathrm{yk}}\right\}\\ +[I_d\⊗{\mathrm{\γ}}_{\mathrm{yk}}-{\mathrm{\γ}}_{\mathrm{yk}}\⊗I_d]\mathrm{vec}\left\{\mathrm{\Δ}{\mathrm{\γ}}_{\mathrm{zk}}\right\}\\ \mathrm{vec}\left\{F_{3(k+1)}\right\}\≈\mathrm{vec}\left\{F_{3k}\right\}\\ +[{\mathrm{\γ}}_{\mathrm{xk}}^T\⊗I_d-I_d\⊗{\mathrm{\γ}}_{\mathrm{xk}}\}\mathrm{vec}\left\{\mathrm{\Δ}{\mathrm{\γ}}_{\mathrm{zk}}\right\}\\ +[I_d\⊗{\mathrm{\γ}}_{\mathrm{zk}}-{\mathrm{\γ}}_{\mathrm{zk}}\⊗I_d]\mathrm{vec}\left\{\mathrm{\Δ}{\mathrm{\γ}}_{\mathrm{xk}}\right\}\end{array}---\left(29\right)$

定义为在第k次迭代时相对于估计的信号子空间初始值的误差矩阵，则改进后的信号子空间估计值可表示为： 

${\tilde{U}}_{\mathrm{sk}}^{\left(\mathrm{nc}\right)}={\hat{U}}_{\mathrm{sk}}^{\left(\mathrm{nc}\right)}+\mathrm{\Δ}{\tilde{U}}_{s,k}^{\left(\mathrm{nc}\right)}---\left(30\right)$

根据以上的迭代公式，NC-SLS-ESPRIT可以用下式表达： 

$\{{\widehat{\mathrm{\γ}}}_x,{\widehat{\mathrm{\γ}}}_y,{\widehat{\mathrm{\γ}}}_z\}=\underset{\mathrm{\Δ}{\mathrm{\γ}}_{\mathrm{\ξk}},\mathrm{\Δ}{\tilde{U}}_{\mathrm{sk}}^{\left(\mathrm{nc}\right)}}{\arg \min }\left|\right|{\left[\begin{array}{c}\mathrm{vec}\left\{R({\tilde{U}}_{\mathrm{sk}}^{\left(\mathrm{nc}\right)}+\mathrm{\Δ}{U}_{\mathrm{sk}}^{\left(\mathrm{nc}\right)},{\mathrm{\γ}}_{\mathrm{xk}})\right\}\\ \mathrm{vec}\left\{R({\tilde{U}}_{\mathrm{sk}}^{\left(\mathrm{nc}\right)}+\mathrm{\Δ}{U}_{\mathrm{sk}}^{\left(\mathrm{nc}\right)},{\mathrm{\γ}}_{\mathrm{yk}})\right\}\\ \mathrm{vec}\left\{R({\tilde{U}}_{\mathrm{sk}}^{\left(\mathrm{nc}\right)}+\mathrm{\Δ}{U}_{\mathrm{sk}}^{\left(\mathrm{nc}\right)},{\mathrm{\γ}}_{\mathrm{zk}})\right\}\\ \mathrm{vec}\left\{F_{1k}\right\}\\ \mathrm{vec}\left\{F_{2k}\right\}\\ \mathrm{vec}\left\{F_{3k}\right\}\\ \widetilde{\mathrm{\κ}}\·\mathrm{vec}\left\{\mathrm{\Δ}{\tilde{U}}_{s,k}^{\left(\mathrm{nc}\right)}\right\}\end{array}\right]+H_k\·\left[\begin{array}{c}\mathrm{vec}\left\{\mathrm{\Δ}{\mathrm{\γ}}_{\mathrm{xk}}\right\}\\ \mathrm{vec}\left\{\mathrm{\Δ}{\mathrm{\γ}}_{\mathrm{yk}}\right\}\\ \mathrm{vec}\left\{\mathrm{\Δ}{\mathrm{\γ}}_{\mathrm{zk}}\right\}\\ \mathrm{vec}\left\{\mathrm{\Δ}{\hat{U}}_{\mathrm{sk}}^{\left(\mathrm{nc}\right)}\right\}\end{array}\right]\left|\right|}_F---\left(31\right)$

其中为使得矢量的优化独立于其他列矢量的调节因子，其中H_k定义如下： 

$H_k=\left[\begin{array}{cccc}{I}_d\⊗\left(K_{x1}^{\left(\mathrm{nc}\right)}{\tilde{U}}_{\mathrm{sk}}^{\left(\mathrm{nc}\right)}\right)& 0_{2(M_y-1)d\×d^2}& 0_{2(M_z-1)d\×d^2}& {\mathrm{\γ}}_{\mathrm{xk}}^T\⊗K_{x1}^{\left(\mathrm{nc}\right)}-I_d\⊗K_{x2}^{\left(\mathrm{nc}\right)}\\ 0_{2(M_x-1)d\×d^2}& I_d\⊗\left(K_{y1}^{\left(\mathrm{nc}\right)}{\tilde{U}}_{\mathrm{sk}}^{\left(\mathrm{nc}\right)}\right)& 0_{2(M_z-1)d\×d^2}& {\mathrm{\γ}}_{\mathrm{yk}}^T\⊗K_{y1}^{\left(\mathrm{nc}\right)}-I_d\⊗K_{y2}^{\left(\mathrm{nc}\right)}\\ 0_{2(M_x-1)d\×d^2}& 0_{2(M_y-1)d\×d^2}& I_d\⊗\left(K_{z1}^{\left(\mathrm{nc}\right)}{\tilde{U}}_{\mathrm{sk}}^{\left(\mathrm{nc}\right)}\right)& {\mathrm{\γ}}_{\mathrm{zk}}^T\⊗K_{z1}^{\left(\mathrm{nc}\right)}-I_d\⊗K_{z2}^{\left(\mathrm{nc}\right)}\\ {\mathrm{\γ}}_{\mathrm{yk}}^T\⊗I_d-I_d\⊗{\mathrm{\γ}}_{\mathrm{yk}}& I_d\⊗{\mathrm{\γ}}_{\mathrm{yk}}-{\mathrm{\γ}}_{\mathrm{yk}}^T\⊗I_d& 0_{d^2\×d^2}& 0_{d^2\×2\mathrm{Md}}\\ 0_{d^2\×d^2}& {\mathrm{\γ}}_{\mathrm{zk}}^T\⊗I_d-I_d\⊗{\mathrm{\γ}}_{\mathrm{zk}}& I_d\⊗{\mathrm{\γ}}_{\mathrm{yk}}-{\mathrm{\γ}}_{\mathrm{yk}}^T\⊗I_d& 0_{d^2\×2\mathrm{Md}}\\ I_d\⊗{\mathrm{\γ}}_{\mathrm{yk}}-{\mathrm{\γ}}_{\mathrm{yk}}^T\⊗I_d& 0_{d^2\×d^2}& {\mathrm{\γ}}_{\mathrm{xk}}^T\⊗I_d-I_d\⊗{\mathrm{\γ}}_{\mathrm{xk}}& 0_{d^2\×2\mathrm{Md}}\\ 0_{2\mathrm{Md}\×d^2}& 0_{2\mathrm{Md}\×d^2}& 0_{2\mathrm{Md}\×d^2}& \widetilde{\mathrm{\κ}}{I}_{2\mathrm{Md}}\end{array}\right]---\left(32\right)$

通过优化式(31)中的代价函数，获得γ_x,γ_y和γ_z的最优估计和对和进行特征分解，计算得到空间频率μ_xp,μ_yp,μ_zp,p＝1,...,d的最优估计值p＝1,...,d。注意在迭代过程中，通过对得到的矢量进行变形处理，能通过下式恢复和γ_ξk: 

$\begin{array}{c}{\tilde{U}}_{s(k+1)}^{\left(\mathrm{nc}\right)}={\tilde{U}}_{\mathrm{sk}}^{\left(\mathrm{nc}\right)}+\mathrm{\Δ}{\hat{U}}_{\mathrm{sk}}^{\left(\mathrm{nc}\right)}\\ {\mathrm{\γ}}_{\mathrm{\ξ}(k+1)}={\mathrm{\γ}}_{\mathrm{\ξk}}+\mathrm{\Δ}{\mathrm{\γ}}_{\mathrm{\ξk}},\mathrm{\ξ}\∈\{x,y,z\}\end{array}---\left(33\right)$

令信号子空间估计的初始值为SVD分解中的信号子空间估计值，即 将利用LS方法获得的作为γ_ξ的初始值。当代价函数满足以下约束条件ξ∈{x,y,z}时，可以认为本发明所采用的改进的SLS-NC-ESPRIT算法收敛到最优解，其中ε为允许的最大误差范围。 

附图2是本发明提供的基于SLS-NC-E ESPRIT算法的3轴交叉阵列DOA估计方法流程图。 

步骤1：利用信号非圆特性得到扩张的接收信号矢量x^(nc)(t)，并计算扩张的采样协方差矩阵

步骤2：利用SVD分解得到信号子空间估计假设改进的信号子空间估计可以表示为：得到残差矩阵

步骤3：利用3轴分叉阵列的阵列矩阵γ_ξ,ξ∈{x,y,z}具有相同的特征向量的特点建立矩阵F₁,F₂和F₃如式(25)所示； 

步骤4：由式(31)建立代价函数，当代价函数满足以下约束条件 ξ∈{x,y,z}时，可以认为收敛到最优解，其中ε为允许的最大误差范围； 

步骤5：由式(31)得到对其做EVD分解，根据式(5)就可以得到d个信源的DOA的估计值。 

本发明主要讨论了对来波信号的空间频率μ_xp,μ_yp,μ_zp,p＝1,...,d的估计值 的推导过程。但是由于每个轴的SVD分解都是独立进行的，所以他们计算得到的空间频率的估计值不是一一对应的。空间频率的配对问题在许多文献中都有讨论，因为不是本发明的核心内容，所以不作详细介绍。通过式(5)发现信号的到达角和俯仰角，即θ_p和p＝1,...d,是由空间频率μ_xp,μ_yp,μ_zp唯一确定的，所以本发明的方法得到的空间频率估计值可以唯一的确定来波信号的DOA。 

考虑一3轴交叉阵列包含M＝3m个阵元，沿x，y和z轴的每个子阵列均有m个阵元。根据式(34)估计RMSE 

${\mathrm{RMSE}}_p=\sqrt{\frac{1}{M_c}\underset{i=1}{\overset{M_c}{\mathrm{\Σ}}}\underset{r=x,y,z}{\mathrm{\Σ}}{({\widehat{\mathrm{\μ}}}_{p,i}^{\left(r\right)}-{\mathrm{\μ}}_p^{\left(r\right)})}^2}---\left(34\right)$

其中蒙特卡洛试验次数为M_c＝1000，p＝1,...,d为第i次试验第p个信号的估计空间频率估计值，而p＝1,...,d是相应的空间频率真实值。在所有的例子中，假设有两个相等功率的BPSK信号，非圆系数均为1，初始相位分别为5°,20°。阵元总数为M＝30,三个子阵列的阵元之间的间隔均相等为0.45λ₀，这里λ₀为来波信号波长。SLS算法的误差容限因子为ε＝10^-6。在仿真试验中为了验证本发明所采用的SLS-NC-ESPRIT算法的性能优势，将其与ESPRIT，SLS-ESPRIT，NC-ESPRIT算法进行了性能比较。为了研究附加约束条件对SLS-ESPRIT和SLS-NC-ESPRIT算法的影响，仿真结果中也提供了没有附加约束条件的SLS-ESPRIT和SLS-NC-ESPRIT的特性，在仿真图中，分别用SLS-ESPRIT(no constraint)和SLS-NC-ESPRIT(no constraint)表示。 

实验1当θ₁＝θ₂,φ₁≠φ₂时的RMSE特性 

当θ₁＝θ₂时，由式(5)可得μ_z1＝μ_z2，即两个信号在z轴上有相同的投影。根据式(21)可知当ξ＝z时会出现秩亏问题，此时LS和没有附加约束条件的SLS方法都会失效，将导致ESPRIT，NC-ESPRIT和没有附加约束条件的SLS-ESPRIT性能差，假设两个BPSK信号的方位角分别为0°和20°，俯仰角均为60°。 

由附图3可以发现，在快拍数固定为N＝300，不同方法的空间频率估计值的RMSE随SNR的变化，随着SNR的变大，没有附加约束条件的SLS-ESPRIT算法性能下降显著，ESPRIT和NC-ESPRIT算法的RMSE也较大，因为它们无法解决秩亏的问题。本发明提出的SLS-NC-ESPRIT算法在所有算法中性能最好。减小快拍数为N＝20，由附图4可知，本发 明提出的SLS-NC-ESPRIT算法的性能仍然是最优的。当固定SNR为-5dB，附图5展示了各算法输出RMSE随快拍数的变化，ESPRIT和没有附加约束条件的SLS-ESPRIT算法性能最差，快拍数增大时，性能几乎没有改善，NC-ESPRIT和没有附加约束条件的SLS-NC-ESPRIT算法虽然性能有所提升，但是仍然不能解决秩亏问题，SLS-ESPRIT和本发明提出的SLS-NC-ESPRIT算法具有较好的性能，但是由于本发明提出的SLS-NC-ESPRIT考虑了信号的非圆特性，其输出RMSE最小。 

实验2当θ₁≠θ₂,φ₁≠φ₂时的RMSE特性 

固定快拍数N＝300，两个BPSK信号的方位角范围是(-4°,4°)，仰角范围是(18°,28°)。由附图6可知，本发明提出的SLS-NC-ESPRIT算法比其他算法性能更优，特别是在信噪比小于-5dB的情况下性能优势更加明显。尽管此时不存在秩亏问题，SLS-ESPRIT和本发明提出的SLS-NC-ESPRIT算法仍然比没有附加约束条件的基于SLS算法具有更小的输出RMSE。在附图7中，减小快拍数为N＝20，此时仍可以发现本发明提出的SLS-NC-ESPRIT性能优于其他算法。附图8展示了各算法输出RMSE随快拍数的变化，利用信号的NC特性的NC-ESPRIT和SLS-NC-ESPRIT算法的性能得到有一定的提高，但是可以发现本发明提出的SLS-NC-ESPRIT算法在所有算法中性能是最优的。 

实验3RMSE随方位角和俯仰角的间隔变化的关系 

固定快拍数N＝300，信噪比SNR＝-5dB，第一个目标的DOAs为(2°,18°)，第二个目标的DOAs为(2°+Δφ,28°)，其中Δφ∈(2°,10°)。由附图9展示了输出RMSE随两信号方位角之间的间隔的变化曲线，由图可以看出,本发明提出的SLS-NC-ESPRIT算法性能优于其他算法。在附图10中，快拍数和信噪比采用与附图9相同的参数设置。固定第一个目标的DOAs为(-4°,18°)，第二个目标的DOAs为(4°,18°+Δθ),Δθ∈(2°,10°)，观察发现当俯仰角小于5°时，本发明提出的SLS-NC-ESPRIT算法的性能得到显著优于其他算法。由于ESPRIT和NC-ESPRIT忽略了子阵列之间结构上的重叠特性，因此它们提供了较差的估计精度，尽管SLS-ESPRIT算法能充分利用子阵列间的重叠结构，但由于没利用信号的NC特性，故它的特性比本发明提出的SLS-NC-ESPRIT算法差，特别是当Δθ＜5°的情况下。 

以上内容是结合具体的优选实施方式对本发明所作的进一步详细说明，不能认定本发明的具体实施只局限于这些说明。对于本发明所属技术 领域的普通技术人员来说，在不脱离本发明构思的前提下，还可以做出若干简单推演或替换，都应当视为属于本发明的保护范围。 

Claims (3)

1.一种基于SLS-NC-ESPRIT算法的3轴交叉阵列的DOA估计方法，所述3轴交叉阵列由三个沿着x轴，y轴，z轴的线性子阵列以原点为几何中心组成，所有的天线阵元都是独立和等距的，3轴的阵元之间距离为δ，阵元数为M＝M_x+M_y+M_z，有d个远场窄带信号到达阵列，其中M_x,M_y,M_z分别表示在x，y，z轴方向的子阵列的包含的阵元个数；其特征在于：所述包括以下步骤： 

1)x(t)表示接收样本数据，利用信号非圆特性，得到扩张的测量数据向量其中Π_M为M×M的交换矩阵，该交换矩阵的反对角元素为1，其它元素为0，A为信号导向矢量矩阵，s(t)为信号矢量，n(t)为噪声； 

2)计算扩张的采样协方差矩阵其中，X^(nc)表示一个2M×N的接收数据矩阵，由N个采样快拍数据x^(nc)(t_n)，1≤n≤N组成； 

3)设ξ∈{x,y,z}，此时3轴交叉阵列具有最大重叠度的选择矩阵如下： 

4)对扩张的采样协方差矩阵做SVD分解得到信号子空间估计 非圆信号的旋转不变等式可表示为： 

对上式进行最小二乘法求解，得到作为γ_ξ的初始估计，信号空间估计的初始值为

5)改进的信号子空间估计可以表示为：得到残差矩阵

6)利用3轴分叉阵列的阵列矩阵γ_ξ,ξ∈{x,y,z}具有相同的特征向量的特点建立矩阵F₁,F₂和F₃： 

F₁＝γ_xγ_y-γ_yγ_x＝0_d×d

F₂＝γ_yγ_z-γ_zγ_y＝0_d×d

F₃＝γ_zγ_x-γ_xγ_z＝0_d×d； 

7)对于第k次迭代，残余矩阵,F₁,F₂和F₃可以表示为： 

F_1k＝γ_xkγ_yk-γ_ykγ_xk

F_2k＝γ_ykγ_zk-γ_zkγ_yk

F_3k＝γ_zkγ_xk-γ_xkγ_zk； 

8)第k+1次的迭代中残余矩阵，F₁,F₂,F₃可写为： 

9)建立代价函数： 

其中，为在第k步骤中估计的信号子空间的误差矩阵，γ_ξk＝γ_ξk-1+Δγ_ξk-1,ξ∈{x,y,z}，vec{·}是向量化函数，

10)判断代价函数是否满足如果满足则认为收敛到最优解，停止迭代，其中ε为允许的最大误差范围；如果不满足，则令k＝k+1，返回步骤9； 

11)将得到的最优解做EVD分解，得到μ_xp,μ_yp,μ_zp,p＝1,...,d的最优估计值其中， 

λ₀为波长，θ_i,φ_i分别表示第i个信号的俯仰角和方位角； 

12)根据d个θ_i,φ_i得到d个信源的DOA的估计值。 

2.根据权利要求1所述的估计方法，其特征在于：利用信号的非圆特性提高了DOA估计精度并使可检测信源数增加，利用3轴分叉阵列阵列矩阵γ_ξ,ξ∈{x,y,z}具有相同的特征向量的特点有效解决秩亏问题，从而提供了更精确的DOA估计。 

3.根据权利要求1所述的估计方法，其特征在于：所述和 的迭代初始值是任给的一个值。