DOA估计

DOA估计（精选7篇）

DOA估计 篇1

智能天线系统是利用共享同一信道的用户之间的空间差异来实现通信容量倍增的, 因此波达方向 (DOA) 估计技术作为获取用户空间位置的重要方法得到了广泛的研究, 取得了一系列高效的处理方法。其中, 经典的有基于子空间技术的多信号分类 (MUSIC) 算法、DOA估计的最优解方法及最大似然 (ML) 算法等。然而这些方法主要是基于一维DOA估计的, 而通信中的用户位于三维空间, 显然二维DOA比一维DOA (方位角) 能更有效地开发信道冗余, 含有更多的信息量, 因此更需要对有效的二维DOA估计算法, 即时空DOA矩阵方法进行研究[1]。

1一种有效的二维DOA估计算法——时空DOA矩阵方法

1.1 阵列信号模型

图1是二维虚拟阵列技术中的阵列结构。该阵列是由位于x轴上阵元间距为dx的M-1元均匀线阵和一个位于y轴上的导引阵元xM组成, xM与x轴的间距为[1]dy。

假设有D个己知波长为λ的窄带信号以DOA{ (α1, β1) , (α2, β2) , …, (αD, βD) }入射, 则各阵元上的观测信号为:

$\begin{array}{l}x{}_i(t)={\displaystyle \sum _{k=1}^{d}s}{}_k(t)\exp \left[-\text{j}\frac{2\text{π}}{\lambda }(i-1)\mathit{d}{}_x\cos \mathit{\alpha }{}_k\right]+n{}_i(t),\\ i=1,2,\cdots ,{\rm M}-1(1)\\ x{}_{{\rm M}}(t)={\displaystyle \sum _{k=1}^{d}s}{}_k(t)\exp \left[\text{j}\frac{2\text{π}}{\lambda }\mathit{d}{}_y\cos \mathit{\beta }{}_k\right]+n{}_{{\rm M}}(t)(2)\end{array}$

式中:sk (t) 为参考阵元x1上接收的信号;ni (t) 代表均值为零, 方差为σ2的加性高斯白噪声[2]。假设sk (t) 之间互不相关, 各阵元间的噪声ni (t) 与sk (t) 之间互不相关, 将式 (1) , 式 (2) 写为矩阵形式:

X (t) =AS (t) +N (t) (3)

A=[α1, α2, …, αD]T, αk=[α1k, α2k, …, αik, …, αMk]T (4)

$\alpha {}_{ik}=\exp \left[-\text{j}\frac{2\text{π}}{\lambda }(i-1)\mathit{d}{}_x\cos \mathit{\alpha }{}_k\right],i=1,2,\cdots ,{\rm M}(5)$

$\alpha {}_{{\rm M}k}=\exp \left[\text{j}\frac{2\text{π}}{\lambda }(i-1)\mathit{d}{}_y\cos \mathit{\beta }{}_k\right](6)$

注意这时αk的结构是Vandermonde矢量[3]。如果得到αk, 就可以求得对应的二维角 (αk, βk) 。

1.2 时空DOA矩阵方法处理[4]

为估计A, 定义xM (t) 与各阵元的相关函数为:

$R{}_{x{}_{i}x{}_{{\rm M}}}=E\{x{}_i(t+\tau /2)x{}_{{\rm M}}^{*}(t-\tau /2)\},i=1,2,\cdots ,{\rm M}(7)$

其中:RX (τ) , RY (τ) 矩阵形式为:

$\begin{array}{l}\mathit{R}{}_X(\tau )=\mathit{A}R{}_S(\tau )\\ \mathit{R}{}_Y(\tau )=\mathit{A}\mathit{\Phi }R{}_S(\tau )\end{array}(8)$

由式 (8) 得:

$\mathit{X}=\mathit{A}\mathit{S},\mathit{Y}=\mathit{A}\mathit{\Phi }\mathit{S}(9)$

其中:

$\begin{array}{l}\mathit{S}=[R{}_S({\rm T}{}_S),R{}_S(2{\rm T}{}_S),\cdots ,R{}_S({\rm N}{\rm T}{}_S)](10)\\ \mathit{R}{}_{{\rm T}{}_S}=\mathit{Y}\cdot [\mathit{X}]{}^{-}(11)\end{array}$

其中:[·]-表示伪逆, 即:

$[\mathit{X}]{}^{-}=\mathit{X}{}^{\text{Η}}(\mathit{X}\mathit{X}{}^{\text{Η}}){}^{-1}(12)$

由$\mathit{R}{}_{{\rm T}{}_S}\overline{\mathit{A}}=\overline{\mathit{A}}\mathit{\Phi }$, 就得到了信号的二维DOA。

1.3 仿真结果与分析

采用如图1所示阵列结构M=4, 阵元间距dx, dy均为0.5波长, 实际快拍数为100, 伪快拍数为50;每一次实验均分别在信噪比为 20 dB下, 进行 50次Monte-Carlo实验。

(1) 有3个来波方向分别为 (50, 50) , (60, 80) , (90, 120) 的窄带信号入射, 仿真结果如图2所示。

这个实验确定该估计方法有效准确。从图2中可以看出该方法准确有效。

(2) 有3个来波方向分别为 (51, 71) , (52, 72) , (53, 73) 的窄带信号入射, 估计结果如图3所示。

这个实验确定该估计方法的分辨率。从图3中可以看出该方法大致可以分辨间隔为1°的两个角。

(3) 有3个来波方向分别为 (40, 50) , (40, 80) , (70, 80) 的窄带信号入射, 估计结果如图4所示。

这个实验确定该估计方法抗角度兼并能力。从图4中可以看出该方法可以处理角度兼并问题[5]。

通过仿真结果可以看出, 时空DOA矩阵是一种非常高效的估计算法, 具有很强的分辨能力, 且能解决二维估计中的角度兼并问题。但是仿真中由于阵元数的限制, 所以该方法也只能分辨比阵元数少一个的信号数[6]。

2 算法的改进

从上述分析可以看出, 虚拟阵列变换的方法可以弥补时空DOA矩阵方法的缺点, 经过分析可将两种方法有效地结合在一起。步骤如下:

(1) 根据时空DOA矩阵方法得到伪快拍数据矩阵X=AS, Y=AΦS;

(2) 根据虚拟阵列变换方法对期望的区间进行阵列变换, 求得变换矩阵$\mathit{B}=\overline{\mathit{A}}\mathit{A}{}^{-1}$;

(3) 可以得到虚拟阵列的伪快拍数据矩阵$\overline{\mathit{X}}=\mathit{B}\mathit{A}\mathit{S}‚\overline{\mathit{Y}}=\mathit{B}\mathit{A}\mathit{\Phi }\mathit{S};$

(4) 得到时空虚拟DOA矩阵, 求得需要的DOA。

根据这一步骤进行了算法仿真[7]:

采用如图1所示阵列结构, M=4, 阵元间距为一个波长, 对于不同的信号个数和变换区域进行了2次仿真实验, 实际快拍数为100, 伪快拍数为50;每一次实验均在信噪比为20 dB下, 进行了10次Monte-carlo仿真[8], 其结果如下:

(1) 5个窄带源分别以入射角度为 (40, 50) , (40, 60) , (43, 70) , (46, 80) , (50, 80) 入射, 变换区域为40°～50°, 虚拟变换后的$\overline{{\rm M}}=7$, 阵元间距为半个波长。DOA估计结果如图5所示。

(2) 6个窄带源分别以入射角度为 (30, 50) , (36, 60) , (42, 70) , (48, 80) , (56, 90) , (62, 100) 入射, 变换区域为30°～70°, 虚拟变换后的$\overline{{\rm M}}=9$, 阵元间距为0.735个波长。DOA估计结果如图6所示。

3 结 论

通过对时空DOA矩阵法数学分析、模型建立及算法的仿真处理, 由仿真结果可以看出[9]:

(1) 时空DOA矩阵法是一种非常高效的估计算法, 具有很强的分辨能力, 且能解决二维估计中的角度兼并问题。

(2) 时空DOA矩阵法的仿真中由于阵元数的限制, 只能分辨比阵元数少一个的信号数。

(3) 对时空DOA矩阵法进行有效的改进, 提出了时空虚拟DOA矩阵方法。该方法可以弥补时空DOA矩阵方法的缺点, 可以估计超过阵元数的信号, 具有较高的分辨性能和较好的稳健性。

参考文献

[1]LIBERTI Joseph C, RAPPAPORT Theodore S.无线通信中的智能天线:IS-95和第3代CDMA应用[M].马凉, 译.北京:机械工业出版社, 2002.

[2]SHANBAZPANAHI Shahram, VABEE Shahrokh, BAS-TANI Mohammad Hasan.Distributed source localizationusing ESPRIT algorithm[J].IEEE Trans.on SignalProessing, 2001, 49 (10) :2164-2178.

[3]金梁, 殷勤业.时空DOA矩阵方法[J].电子学报, 2000, 28 (6) :8-12.

[4]GREEN Peter J, TAYLOR Desmond P.Dynamic signalenumeration algorithm for smart antennas[J].IEEETrans.on Signal Processing, 2002, 50 (6) :1307-1314.

[5]吴伟陵.移动通信中的关键技术[M].北京:北京邮电大学出版社, 2001.

[6]GODARA L C.无线通信天线手册[M].左群声, 金林, 胡明春, 等译.北京:国防工业出版社, 2004.

[7]杨大成.现代移动通信中的先进技术[M].北京:机械工业出版社, 2005.

[8]谢显中.TD-SCDMA第三代移动通信系统技术与实现[M].北京:电子工业出版社, 2004.

[9]杨维, 陈俊仕, 李世明, 等.移动通信中的阵列天线技术[M].北京:清华大学出版社, 2005.

[10]王磊, 束坤.一种子空间识别DOA估计算法的研究[J].现代电子技术, 2010, 33 (11) :29-31.

[11]GHAVAI Mohammad.Wideband smart antenna theoryusing rectangular structures[J].IEEE Trans.on SignalProcessing, 2002, 50 (9) :2143-2151.

利用盲源分离算法实现DOA估计 篇2

波达方向(DOA)是信号源的一个重要参数,利用传感器阵列接收数据并实现波达方向估计,其在雷达、通信、声纳、地震探测等阵列系统中具有极其重要的地位[1,2]。盲源分离所要解决的问题是[3],在源信号和传输信道参数等先验知识未知的情况下,仅根据输入源信号的统计特性估计信道参数、恢复源信号[3,4,5,6,7,8]。由于通过盲源分离算法解出的信道参数中,蕴含着波达方向信息,因此,可以考虑利用盲源分离算法实现DOA估计[9,10]。

Pourrostam等人将解盲源分离问题的传统二阶盲辨识(SOBI)算法[4]引入DOA处理[9](记作SOBI-DOA),实现了DOA估计,但是SOBI算法本身所需要的预白化操作,降低了算法的应用范围并导致了误差[7]。聂卫科等人利用不需要预白化操作的三迭代(TIA)算法[5]估计DOA[10](记作TIA-DOA),估计性能有了提高,但是该算法要求阵元数严格大于信源数,且三迭代(TIA)算法本身存在左右混迭矩阵估计不一致的情况,从而会带来DOA的估计误差[6]。为了克服以上缺陷,本文将一种快速复数域盲源分离算法[8]引入阵列信号处理,实现DOA估计(记作FBSS-DOA)。算法充分利用所构造相关矩阵组的对角结构,引入盲源分离领域常用的联合对角化代价函数[3,4,8,9],并利用一种快速的复数域乘性迭代算法求解代价函数,实现矩阵组的联合对角化,求得混迭矩阵。进一步通过深入挖掘混叠矩阵中蕴含的波达方向信息,求解DOA参数。与SOBI-DOA[9] 算法相比,FBSS-DOA算法具有不需要预白化操作,避免了预白化的缺点;而与严格限制阵元数大于信源数的TIA-DOA[10]算法相比,FBSS-DOA算法只限制阵元数不小于信源数。可见,本文提出的FBSS-DOA算法具有更广的适用性及更精确的DOA估计性能。

1 数据模型

考虑有M个全向阵元的等距线阵,阵元间隔为d(d≤λ/2,λ表示信号波长),接收N个(M≥N)入射方向为θ=(θ1,θ2,…,θN)的远场窄带信号,第u个阵元接收信号可表示为:

$\begin{array}{l}x{}_u(t)={\displaystyle \sum _{v=1}^{{\rm N}}s}{}_v(t)\text{e}{}^{\text{j}\frac{2\text{π}}{\lambda }(u-1)d\text{s}\text{i}\text{n}\theta {}_v}+n{}_u(t)\\ u=1,2,\cdots ,{\rm M}(1)\end{array}$

将式(1)表示为矩阵形式:

$\mathit{x}(t)=\mathit{A}\mathit{s}(t)+\mathit{n}(t)(2)$

式中:阵列流型矩阵A=[a(θ1),a(θ2),…,a(θN)]是列满秩的,矩阵A的列向量具有形式a(θv)=[1,ej2π(d/λ)sin θv,…,ej2πM(d/λ)sin θv]T,v=1,2…,N,很显然的,矩阵A中蕴含着波达方向θ=(θ1,θ2,…,θN)信息。式(2)中的s(t)=[s1(t),s2(t),…,sN(t)]T,表示均值为零,两两互不相关的源信号,即:

$\begin{array}{l}E[\mathit{s}(t)\mathit{s}{}^{\text{Η}}(t+\tau )]=\text{d}\text{i}\text{a}\text{g}[\rho {}_1(\tau ),\rho {}_2(\tau ),\cdots ,\rho {}_{{\rm N}}(\tau )]\\ =\mathit{R}{}_{\mathit{s}}(\tau )(3)\end{array}$

式中diag[·]表示由向量组成其对角线元素的对角阵, ρv(τ)=E[sv(t)sHv(t+τ)],v=1,2,…,N,τ表示时延。n(t)=[n1(t),n2(t),…,nM(t)]T表示均值为零方差为σ${}_n^2$的阵元噪声,即:

$R{}_n=E[\mathit{n}{}_p(t)\mathit{n}{}_q^{\text{Η}}(t)]=\sigma {}_n^2\delta {}_{p,q}‚p,q=1,2,\cdots ,{\rm M}(4)$

为了分析的直观,忽略噪声的影响,取接收信号x(t)在K个不同的非零时延kΔτ(Δτ为时延步长,k=1,2,…,K)下的二阶相关矩阵Rkx:

$\begin{array}{l}\mathit{R}{}_{\mathit{x}}^k=E[\mathit{x}(t)\mathit{x}{}^{\text{Η}}(t+k\text{Δ}\tau )]=\mathit{A}\mathit{R}{}_{\mathit{s}}^k\mathit{A}{}^{\text{Η}}=\mathit{A}\mathit{\Lambda }{}^k\mathit{A}{}^{\text{Η}},\\ k=1,2,\cdots ‚{\rm K}(5)\end{array}$

式中:Rks的含义可由式(3)类推;Λk是对角矩阵Rks的简化表示,通过取k=1,2,…,K表征不同的非零延时,将Rkx简记为Ck,则可得到K个具有对角化结构的矩阵组,简记为:C={Ck=AΛkAH,k=1,2,…,K},称为目标矩阵组。在实际应用中,DOA估计所要解决的问题是:根据线阵传感器输出的T个样本{x(t)}${}_{t=1}^{{\rm T}}$,估计信源波达方向角θ=(θ1,θ2,…,θN)。由于噪声的影响以及样本数目的有限性,导致矩阵组C中的相关矩阵,只具有近似的对角化结构。

观察式(2)并考虑相应的若干假设可知,式(2)表征的数据模型与盲源分离问题的数据模型完全一致,并且,在解盲源分离问题的诸多算法中,很多针对复数域问题的盲源分离算法[3,4,5,7,8]都在此数据模型下,产生一组具有近似对角化结构的矩阵组{Ck=AΛkAH,k=1,2,…,K},并通过复数域的联合对角化思想求解阵列流型矩阵A(在盲源分离领域称为混迭矩阵)的估计。因而,面对DOA估计问题,在得到具有近似对角化结构的目标矩阵组C后,可以利用解盲源分离思路,实现对混迭矩阵A的估计,然后再挖掘求解A中蕴含的波达方向信息,进而估计DOA,这正是本文的主要研究思路。

考虑到当M>N时,A为非方阵,而盲源分离中的联合对角化算法一般假设混迭矩阵A为方阵,为了便于计算,并降低运算量,根据文献[3]所述方法对目标矩阵组C进行降维处理:令$\mathit{\psi }={\displaystyle \sum _{k=1}^{{\rm K}}\mathit{C}}{}^k(\mathit{C}{}^k){}^{\text{Η}}$,对ψ进行奇异值分解ψ=UDVH,取降维矩阵Γ为ψ的N个主特征矢量组成的矩阵,即Γ=[u1,u2,…,uN]∈M×N。令$\overline{\mathit{A}}=\mathit{\Gamma }{}^{\text{Η}}\mathit{A}$,则可得到降维后的矩阵组$\overline{\mathit{C}}=\{\overline{\mathit{C}}{}^k=\overline{\mathit{A}}\mathit{\Lambda }{}^k\overline{\mathit{A}}{}^{\text{Η}},k=1,2‚\cdots ,{\rm K}\}$。显然地,降维处理并没有破坏目标矩阵组的对角化结构,并且,据此方法得到的降维矩阵Γ还具有抑制噪声的作用[3]。

需要指出的是,TLS-ESPRIT算法是一种传统而经典的DOA估计方法[2],而SOBI-DOA算法[9],TIA-DOA算法[10]以及本文的FBSS-DOA算法,都采用了利用盲源分离算法实现DOA估计的思路。因此,比较最终的估计效果是体现各算法性能的最直接方式。具体到SOBI-DOA算法,TIA-DOA算法以及本文的FBSS-DOA算法而言,从上文所述的数据分析角度看,SOBI-DOA算法要求阵列流型矩阵A是酉矩阵,否则便需要根据矩阵组C求得白化矩阵,实现所求混迭矩阵A的预白化,由于采样点数目的有限性导致矩阵组C中的相关矩阵的结构本身就有误差存在,因而导致据此求得的白化矩阵具有误差,而白化过程中引入的误差无法在后续的算法中得以纠正[7],因此降低了算法的应用范围及估计性能;与上述的通过直接求解相关矩阵的方法获得可联合对角化目标矩阵组C的方式不同,TIA-DOA算法需要划分子阵并利用两个子阵的旋转不变性,以得到该所需的目标矩阵组C,因此必需要求阵元数目严格大于信源数目(M>N),另外,在求解混迭矩阵A的估计的过程中,由于算法本身的限制,需要将式(5)中的混迭矩阵A分别看作左混迭矩阵与右混迭矩阵,并需要将左右混迭矩阵看作两个不同的变量分别求解,这就可能导致最终的左右混迭矩阵估计结果不一致的情况,产生估计误差。本文算法克服了SOBI-DOA算法及TIA-DOA算法的上述缺陷,因而保证了算法极广的适用性及良好的估计性能。

2 FBSS-DOA算法描述

为了求得$\overline{\mathit{A}}$,先求其逆矩阵$\mathit{V}=\overline{\mathit{A}}{}^{-1}$。求解V通过对矩阵组$\overline{\mathit{C}}$的联合对角化实现:希望求得的V能使得矩阵组$\overline{\mathit{C}}$中的各矩阵最大程度的对角化,即使得所有的$\mathit{V}\overline{\mathit{C}}{}^k\mathit{V}{}^{\mathit{{\rm H}}},k=1,2,\cdots ,{\rm K}$尽量为对角化矩阵。为了实现这一目标,首先建立表征联合对角化近似程度的代价函数[3,4,8]:

$L(V)=\underset{V\in \mathit{C}{}^{{\rm N}\times {\rm N}}}{{\displaystyle \min }}{\displaystyle \sum _{k=1}^{{\rm K}}\text{o}}\text{f}\text{f}(\mathit{V}\overline{\mathit{C}}{}^k\mathit{V}{}^{\text{Η}})(6)$

式中:off(V$\overline{\mathit{C}}$kVH) 表示矩阵V$\overline{\mathit{C}}$kVH非对角线元素的F-范数。

引入我们新近提出的一种采用了乘性迭代机制的快速复数域盲源分离算法[8]求解代价函数(6),并实现DOA估计,算法简称为FBSS-DOA,描述如下:

(1) 初始化V(0)=I;

(2) 进入迭代流程,在第n步迭代时,令V(n)具有形式V(n)=(I+W(n))V(n-1),其中I为单位阵,设定更新矩阵W(n)对角线元素为零并被限定为:

$\mathit{W}{}_{(n)}\leftarrow \frac{\alpha }{\parallel \mathit{W}{}_{(n)}\parallel {}_{\infty }}\mathit{W}{}_{(n)},0<\alpha <1$

以保证I+W(n)具有对角占优的特性,根据对角占优定理[11],I+W(n)的可逆性得以保证,进而保证了V(n)的可逆性,防止算法收敛到平凡解V=0或其他非可逆解;第n步迭代的目标是,求解更新矩阵W(n),将$\overline{\mathit{C}}$${}_{(n)}^k$,k=1,2,…,K转化成$\overline{\mathit{C}}{}_{(n+1)}^k,k=1,2‚\cdots ,{\rm K}$,并使得K个$\overline{\mathit{C}}$${}_{(n+1)}^k$都尽量对角化;更新矩阵W(n)的求解,可以通过分步求解其所有关于对角线对称的成对元素实现,求得W(n)后可求得V(n)。

(3) 算法收敛后得到V,在此基础上,可进一步得到阵列流型矩阵A的估计:$\tilde{\mathit{A}}=(\mathit{\Gamma }{}^{\text{Η}}){}^{+}\mathit{V}{}^{-1}‚(\cdot ){}^{+}$表示伪逆。

(4) 设au,v为$\tilde{\mathit{A}}$的第(u,v)元素,av为$\tilde{\mathit{A}}$的第v列向量,取$\tilde{\mathit{A}}\leftarrow \left[\frac{\mathit{a}{}_1}{a{}_{1,1}},\frac{\mathit{a}{}_2}{a{}_{1,2}},\cdots ,\frac{\mathit{a}{}_{{\rm N}}}{a{}_{1,{\rm N}}}\right]$,则$\tilde{\mathit{A}}$便具有了与式(2)中的A完全一致的标准形式。文献[8]提出只利用$\tilde{\mathit{A}}$的最后一行元素求解各个信源的波达角度,经过实验验证表明,采用文献[9]提出的求解第j个信源波达方向的方法更具准确性和稳健性:

$\begin{array}{l}\widehat{\mathit{\theta }}{}_v=\frac{180}{({\rm M}-1)\text{π}}{\displaystyle \sum _{u=1}^{{\rm M}-1}\left\{\arcsin \left[\frac{\lambda }{2\text{π}d}\text{a}\text{n}\text{g}\text{l}\text{e}(a{}_{u,v}^{*}a{}_{u+1,v})\right]\right\}}\\ j=1,2,\cdots ,{\rm N}(7)\end{array}$

式中:angle(·)表示复角主值。基于此,根据矩阵A的各列,可得波达方向估计结果$\widehat{\mathit{\theta }}=(\widehat{\theta }{}_1,\widehat{\theta }{}_2,\cdots ,\widehat{\theta }{}_{{\rm N}})$。

3 仿 真

本节通过三个实验来分析FBSS-DOA算法的性能。首先给出三个性能指标。显然,矩阵组$\overline{\mathit{C}}=\{\overline{\mathit{C}}{}^k=\overline{\mathit{A}}\mathit{\Lambda }{}^k\overline{\mathit{A}}{}^{\text{Η}},k=1,2,\cdots ,{\rm K}\}$的联合对角化程度能直接表征所使用的盲源分离算法的有效性[7]。用对角化误差(Diagonalization Error,DE)来描述目标矩阵组的联合对角化程度:

$\text{D}\text{E}={\displaystyle \sum _{k=1}^{{\rm K}}\text{o}}\text{f}\text{f}(\mathit{V}\overline{\mathit{C}}{}^k\mathit{V}{}^{\text{Η}})={\displaystyle \sum _{k=1}^{{\rm K}}{\displaystyle \sum _{u\ne v}|}}\mathit{V}\overline{\mathit{C}}{}^k\mathit{V}{}^{\text{Η}}|{}_{uv}^2(8)$

可见,DE的定义方式与代价函数(6)一致。

同样地,阵列流型矩阵A的估计精度直接影响波达方向的估计精度,用全局拒噪水平(Global Rejection Level,GRL)来描述A的估计误差,定义为[4,5,7,8]:

$\begin{array}{l}\text{G}\text{R}\text{L}={\displaystyle \sum _{u=1}^{{\rm N}}\left({\displaystyle \sum _{v=1}^{{\rm N}}\frac{\left|g{}_{uv}\right|}{\max {}_k\left|g{}_{uk}\right|}}-1\right)}+\\ {\displaystyle \sum _{v=1}^{{\rm N}}\left({\displaystyle \sum _{i=1}^{{\rm N}}\frac{\left|g{}_{uv}\right|}{\max {}_k\left|g{}_{kv}\right|}}-1\right)}(9)\end{array}$

式中:guv为矩阵G=VΓHA的第u行第v列元素。

第三个性能指标为角度估计的均方根误差(RMSE):设独立实验次数为Q次,定义$\widehat{\mathit{\theta }}$v(q)为第q次试验对第v个信源波达方向的估计,θv为第v个信源的真实波达方向,定义角度估计的均方根误差为[10]:

$\text{R}\text{Μ}\text{S}\text{E}=\sqrt{\frac{1}{Q}{\displaystyle \sum _{q=1}^Q\left[\frac{1}{{\rm N}}{\displaystyle \sum _{v=1}^{{\rm N}}(}\widehat{\mathit{\theta }}{}_v(q)-\mathit{\theta }{}_v){}^2\right]}}(10)$

在以下3个实验中,设定有6个相距半波长的均匀线阵接收3个波达方向为θ=(-5°,0°,5°)T的远场窄带信号,取16个非零时延相关矩阵,乘性迭代算法的初值取单位阵V(0)=I,并令正常数α=0.9。

实验1:本实验验证FBSS-DOA算法的有效性及快速收敛性。设定信噪比SNR=5 dB,快拍数T=200,运行20次独立实验。对角化误差随迭代次数变化的性能曲线如图1所示。全局拒噪水平随迭代次数变化的性能曲线如图2所示。并可求得角度估计的均方根误差RMSE=-15.82 dB。图1和图2所示,以及角度估计的均方根误差说明,FBSS-DOA算法收敛快速,且有好的估计性能。

实验2:本实验分别考察信噪比,快拍数对FBSS-DOA迭代算法收敛性能的影响。表1所示为T=200,不同信噪比情况下,算法收敛所需的迭代次数及GRL收敛值。表2所示为SNR=5 dB,不同快拍数情况下,算法收敛所需的迭代次数及GRL收敛值。表1及表2所示均为100次Monte Carlo独立实验的平均结果(由于是平均值,因而迭代次数产生了小数)。表1表2所示结果表明了迭代算法的良好性能。

实验3:本实验分别考察信噪比,快拍数对FBSS-DOA算法DOA估计性能的影响,并与TLS-ESPRIT[2]算法,SOBI-DOA算法[4]以及TIA-DOA算法[9]的估计性能作出比较。图3为T=200时,四种算法的角度估计的均方根误差随信噪比变化曲线,图4为SNR=5 dB时,角度估计的均方根误差随快拍数变化曲线,所有结果均为100次Monte Carlo独立实验的平均值。图3图4所示结果说明,FBSS-DOA算法的估计性能优于其余三种算法。

4 结 语

本文首先构造一组具有对角化结构的空时相关目标矩阵组,然后利用现有的科学方法实现降维,降低了运算量并在一定程度上抑制了噪声。随后,引入解盲源分离的联合对角化思路,建立代价函数并采用了一种快速的复数域乘性迭代算法以求解代价函数,得到混迭矩阵的估计,进而实现DOA估计。仿真实验表明,所提出的FBSS-DOA算法与同类算法相比,具有更精确的估计性能。

参考文献

[1]DMOCHOWSKI J,BENESTY J,AFFES S.An informa-tion-theoretic view of array processing[J].IEEE Transac-tions on Audio,Speech,and Language Processing,2009,17(2):392-401.

[2]TRIPATHY P,SRIVASTAVA S C,SINGH S N.A mo-dified TLS-ESPRIT-based method for low-frequency modeidentification in power systems utilizing synchrophasormeasurements[J].IEEE Transactions on Power Systems,2011,26(2):719-727.

[3]FENG D Z,ZHENG W X,CICHOCKI A.Matrix-groupalgorithm via improved whitening process for extracting sta-tistically independent sources from array signals[J].IEEETransactions on Signal Processing,2007,55(3):962-977.

[4]BELOUCHRANI A,ABED-MERAIM K,CARDOSO J F,et al.A blind source separation technique using second-or-der statistics[J].IEEE Transactions on Signal Processing,1997,45(2):434-444.

[5]张华,冯大政,聂卫科,等.非正交联合对角化盲源分离算法[J].西安电子科技大学学报,2008,35(1):27-31.

[6]LATHAUWER L D,CASTAING J.Blind identification ofunderdetermined mixtures by simultaneous matrix diagona-lization[J].IEEE Transactions on Signal Processing,2008,56(3):1096-1105.

[7]SOULOUMIAC A.Nonorthogonal joint diagonalization bycombining givens and hyperbolic rotations[J].IEEETransactions on Signal Processing,2009,57(6):2222-2231.

[8]XU X F,FENG D Z,ZHENG W X.a fast algorithm fornon-unitary joint diagonalization and its application to blindsource separation[J].IEEE Transactions on signal proce-ssing,2011,59(7):3457-3463.

[9]POURROSTAM J,ZEKAVAT S A,POURKHAATOUNM.Super-resolution direction-of-arrival estimation via blindsignal separation methods[C]//IEEE Radar Conference.Boston,MA:IEEE,2007:614-617.

[10]聂卫科,冯大政,张华.空时矩阵组三迭代DOA估计新方法[J].系统工程与电子技术,2008,30(3):404-408.

DOA估计 篇3

ESPRIT方法[1]现在已经成为信号DOA和极化参数估计的重要方法,与经典的MUSIC算法相比,ESPRIT算法不需要进行谱峰搜索,所以大大降低了计算和存储的代价。ESPRIT算法的关键是对接收数据的协方差矩阵进行特征分解,进而估计出信号的参数。然而在实际运用中,信号的参数往往是时变的,因此必须对协方差矩阵进行实时更新,但若对每一个样本都进行特征分解,所需运算量太大,难以实现。

Liu等人提出的自适应ESPRIT算法[2]是建立在秩显URV分解[3]的基础上的,使得运算量大大减少,且可以有效地并行实现,王雪松等人提出的极化域-空域联合谱估计[4]较传统的MUSIC方法运算量大大减少,但是仍需进行谱峰搜索。

本文在信号序贯处理的基础上,提出了信号相关矩阵的递推估计,并且根据矩阵扰动理论,利用矩阵特征分解的一阶修正法,实现了特征值和特征向量的递推估计,从而使得ESPRIT算法能够自适应地估计信号的参数,并且由于不需要每次都进行特征分解,其运算量也大大减少了。

2 ESPRIT 算法

考虑一个由M个阵元组成的极化敏感阵列,如图1所示。每个阵元都由相互正交的偶极子对组成,阵元间距为δ。假设有K个连续波信号入射到阵列(K

信号k的极化矢量表达式为:

undefined

定义相位延迟因子:undefined。

则阵列接收信号的矢量表达式为:

undefined

其中A=[a1,a2,…aK]。

undefined

各信道噪声相互独立,且与信号不相关。协方差阵为:

undefined

其中R0=ARsA,Rs为信号自相关阵,σ2为噪声功率,对R进行特征分解,假设其特征值为:

undefined

对应的特征向量为e1,e2,…e2M。

由此可得信号子空间为Es=[e1e2 … eK],因此存在非奇异矩阵T,满足Es=AT。

让Aq1,Aq2分别表示矩阵A的前2(M-1)和后2(M-1)行组成的子矩阵,则:Aq2=Aq1φ,其中φ=diag(q1,q2,…qK)。

Eq1,Eq2分别表示信号子空间的前2(M-1)和后2(M-1)行组成的子空间,则:

undefined

则φ中的元素就是矩阵Ψ=T-1φT的特征值,Ψ满足Eq2=Eq1Ψ。

假设undefined,让Ar1,Ar2分别表示矩阵A的偶数行和奇数行组成的子矩阵,则Ar2=Ar1φr,φr=diag(r1,r2,…rK)。

Er1,Er2分别表示信号子空间的偶数行和奇数行组成的子矩阵,则:

undefined

则φr的元素就是Ψr=T-1φrT的特征值,Ψr满足Er2=Er1Ψr。

在多个信号的情况下,我们还需要判断Ψ,Ψr中特征值的对应关系,假设Ψq=Ψr-1Ψ,则Ψq的特征值是Ψ和Ψr特征值的比值。这样Ψ和Ψr的对应特征值应满足使下式最小:

undefined是Ψq的特征值

在估计出(qk,rk)后,就可以估计信号的参数了。

undefined

3 阵列协方差的更新

根据阵列采样数据{x(k),k=1,2,…,N}可得阵列协方差矩阵的极大似然估计为:

undefined

对于非平稳情形,利用当前观测数据x(k)对前一时刻阵列协方差矩阵R(k-1)更新得到当前时刻阵列协方差矩阵R(k)。更新过程可以表示如下:

undefined

其中:ε 为遗忘因子(0<ε<1),ε提供了指数窗有效长度的一个粗略测度。上式重写为:

undefined

其中:Rk为当前观测数据带来的新息矩阵,定义为:

undefined

可以看出,对于足够小的ε,当前的阵列协方差矩阵R(k)可以解释为R(k-1)经过微小扰动εRk得到,由于扰动项由ε线性表出,所以称为一阶扰动。

4 矩阵特征分解的一阶修正

假定已经获得R(k-1)的特征值和特征向量对{λi(k-1),ui(k-1),i=1,2,…,2M},根据k-1时刻的特征分解可以得到:

R(k-1)ui(k-1)=λi(k-1)ui(k-1)

ui(k-1)Huj(k-1)=δij

其中:δij为第拉克函数。假定R(k)的特征值和特征向量对为{λi(k),ui(k),i=1,2,…,2M},则根据矩阵扰动理论[5]{λi(k),ui(k),i=1,2,…,2M}可以用ε的幂级数展开,并且当ε→0 时,他们将分别收敛到{λi(k-1),ui(k-1),i=1,2,…2M}。存在幂级数:

undefined

根据k时刻的特征分解可知:

undefined

则由上面三个式子可以得到undefined的表达式。下面不加推导的给出其表达式。具体推导可参见文献[6]。

undefined

其中yi=ui(k-1)Hx(k),为k时刻的观测向量x(k)在k-1时刻特征向量ui(k-1)上的投影坐标。

undefined

其中,

undefined

以上修正式是在没有重根的情况下得到的,对于有重根的情况必须进一步修正[7]。

undefined

其中,undefined,噪声子空间的第一个特征矢量集中了噪声的全部能量,所以只需考虑第一个特征向量为:undefined,并且当j=K+2,…,2M时,bij=0。

5 仿真分析

(1) 阵列由3个阵元组成,间距为δ=λ/2,排列方式为沿y轴。

假设有一个信号入射到阵列,其参数为:

undefined

持续时间为500 s,采样周期为1 s,SNR为10 dB,信道噪声为0均值高斯噪声,与信号不相关。仿真结果如图2所示。

从图1可以看出在经过30次左右的适应之后,便能很好地跟踪信号参数的变化,且信号的DOA和极化参数都能够较好地估计出来。

(2) 假设由两个信号入射到阵列,信号的参数为:

undefined

(SNR1=10 dB,SNR2=20 dB)。其余条件和图1中一样,仿真结果如图3,图4所示。

从图3,图4中可以看出文中方法的有效性,信号2参数的仿真结果较好,但与(1)的仿真结果相比,误差较大。这主要是因为阵列的阵元数太少,如果加大阵元同样能够取得较好的结果。对于(2)中极化参数的仿真误差较大,这主要是由于信号DOA 估计误差增大的原因,同样也有阵元数较少的原因,这主要是为了与(1)的结果相比较。

6 结 语

ESPRIT算法可以有效地估计出信号的DOA 和极化参数,但是面临着矩阵特征分解的问题。如果阵元较多,计算量将非常大,特别是需要对信号进行实时估计的情况。文中针对这一问题,利用矩阵特征分解一阶修正法,使得矩阵特征值和特征向量能够进行递推,从而实现了自适应估计,且计算量大大减少,有利于实时估计。经过总结认为还有以下方面可以完善:运用二阶修正应该能够减小误差;加大阵元数应该能够取得更好的结果。

参考文献

[1]Roy R,Kailath T.ESPRIT-Estimation of Signal Parametersvia Rotational Invariance Techniques[J].IEEE Trans.onAcoust.,Speech,Signal Processing,1989,37:984-995.

[2]Liu KJ R,Leary D P O,Stewart GW,et al.URV Esprit forTracking Time-Varying Signals[J].IEEE Trans.on SignalProcessing,1994,42:3 441-3 448.

[3]Stewart G W.An Updating Algorithm for Subspace Track-ing[J].IEEE Trans.on Signal Processing,1992,40(6):1 535-1 541.

[4]徐振海,王雪松,肖顺平.极化域-空域联合谱估计[J].国防科技大学学报,2004,26(3):63-67.

[5]张贤达.信号处理中的线性代数[M].北京:科学出版社,1997.

[6]徐振海,王雪松,肖顺平.矩阵特征分解二阶修正算法[J].信号处理,2004,20(6):600-604.

DOA估计 篇4

目前,精确定位技术是武器系统的关键技术之一,而目标方位角和俯仰角的确定首当其冲[1]。在信号处理过程中,MUSIC算法由于具有较高的分辨力、估计精度和稳定性而备受青睐[2]。另外,作为信号处理领域内的一个重要分支,阵列信号处理在近30年来得到迅速发展,广泛应用在雷达、通信、声纳、地震勘探、射电天文以及生物医学工程等众多军事及国民经济领域。利用阵列信号处理技术对波达方向的估计(也称空间谱估计)也一直是人们研究的热点。从20世纪70年代末开始,在空间谱估计方面涌现出了大量的研究成果和文献,其中以美国的Schmidt R.O等人提出的多重信号分类(MUSIC)算法最为突出。MUSIC算法的特点就是通过对阵列接收数据的数学分解,将接收的数据划分为相互正交的信号子空间和噪声子空间,利用其正交特性构造出“针状”空间谱峰,从而提高算法的分辨力。

1 信号模型及MUSIC算法

天线阵列有多种,均匀线性阵列结构简单,分析和工程实现都比较容易。均匀线阵的模型[3](见图1),空间信号的全部信息包含在阵列信号矢量x(t)中,MUSIC算法的基本思想是先对信号协方差矩阵进行特征分解,将信号子空间和噪声子空间分离,然后进行信号子空间匹配完成空间谱估计(角度估计)。具体算法为:

假设空间有N个信号(s1(t),s2(t),…,sN(t))入射到m元均匀线阵,各信号源入射角分别为:θ1,θ2,…,θN,各阵元噪声为: n1(t),n2(t),…,nm(t),假定噪声是相互独立、具有相同的功率、空间平稳的高斯白噪声,则第i个阵元输出为:

$x{}_i(t)={\displaystyle \sum _{k=1}^{{\rm N}}s}{}_k(t+(i-1)\tau {}_k)+n{}_i(t)$。

式中:i=1,2,…,m,$\tau {}_k=\frac{d\text{s}\text{i}\text{n}\theta {}_k}{c}$为各信号在相邻阵元之间的相对延迟;d为相邻阵元间距;c为光速。

在窄带假设下有:

$x{}_i(t)={\displaystyle \sum _{k=1}^{{\rm N}}s}{}_k(t)e{}^{j\omega {}_0(i-1)\tau {}_k}+n{}_i(t)$。

信号向量张成的子空间为:

A(θ)=[a(θ1),a(θ2),…,a(θN)]。

其中:信号向量a(θk),也称信号流型(k=1,2,…,N):

a(θk)=[1,e-jω0τk,…,e-j(m-1)ω0τk]T。

各阵元输出(以矩阵形式表示)为:

x(t)=A(θ)s(t)+n(t)。

因此,阵列协方差矩阵R为:

R=E[x(t)x(t)H]=A(θ)RsA(θ)H+σ2I。

式中:Rs是信号协方差矩阵,并且Rs=E[s(t)s(t)H],σ${}_{}^2$为高斯白噪声方差值,I为单位矩阵 ,{}H为复共扼转置运算。

对阵列协方差矩阵R进行特征分解就可以得到N个较大的特征值和m-N个等于σ${}_{}^2$的特征值。按照特征结构的性质,若m>N,将特征值按降序排列,前N个最大特征值对应的特征向量构成的矩阵其张成的空间称为信号子空间,后m-N个最小特征值对应的特征向量构成的矩阵其张成的空间称为噪声子空间。假设各信号非相关,则信号子空间和A(θ)的列向量a(θk)所张成的子空间L(A)相同。信号特征向量张成的空间和噪声向量张成的空间正交,亦信号空间的任意向量和噪声向量正交,有:

aH(θk) vj=0。

式中:j=N+1,N+2,…,m,利用这一特性,可以对R矩阵进行特征空间分解,得出信号空间及噪声空间,根据信号参数进行谱峰搜索,极大值点所对应的角度就是信号入射方向。

2 误差分析

影响目标方位角估计的误差因素大体上有三个部分:建模引起的误差、信噪比引起的误差、阵元数引起的误差等等。

2.1 建模误差

实际系统中的模型误差总是存在的,测向定位技术对模型误差的敏感性使得系统不容易实现。模型误差主要有通道失配误差、阵元间互耦误差和阵元位置误差,任一误差因素都会对系统精度产生很大的影响。建模引起的误差是不可避免的,在实际中要多加考虑,只能对误差做全面的补偿。另外,实际过程阵列的校正也不一定精确,各通道的幅相也不一定一致。

2.2 信噪比引起的误差

信噪比降低的影响有两个方面:第一,信噪比降低导致信号源数目估计与真实信号源数目不一致时,也就是信号子空间、噪声子空间估计不准使得信号子空间和噪声子空间的划分错误而造成虚警或漏报,也就会在估计信号方向时有偏差;第二,实际计算时,只能得到协方差估计矩阵,而噪声电平的提高会破坏由协方差估计矩阵分解得到的信号子空间和噪声子空间的正交性,同样会使角度估计的分辨率下降,估计偏移量及方差增加。

在信噪比分别-20 dB、-15 dB、-10 dB、-5 dB、0 dB、5 dB、10 dB、20 dB时进行仿真,仿真计算采用8元均匀线阵,空间存在4个互不相关的窄带源,其方位角分别为:16°、0°、-17°、39°。为了仿真计算方便,将窄带信号近似等同于正弦信号,其表达式为:x(t)=cos ω0t,根据雷达接收的实际情况,X波段[(8~12.5) GHz]是机载雷达的常用工作范围,工作最短波长λmin =0.024 m,取阵元间距$d=\frac{\lambda {}_{\min }}{2}=0.012\text{m}$,取信号回波的频率分别为8 GHz、9 GHz、10 GHz、12 GHz,空间采样数为NX=1 280。对于噪声信号,可取均值为零,方差为σ2的高斯白噪声。

MUSIC算法估计角度(见图2),峰值对应的横坐标就是目标的方位角,从图中可以看出在SNR=-10 dB的时候,角度的分辨率还是很高的。

信噪比误差的仿真结果(见图3),由图可以看出,MUSIC算法随着信噪比的增加分辨率也相应逐渐提高。在仿真过程中,在信噪比为-20 dB的时候,两信号源的谱峰出现了叠加的现象,几乎不能分辨。由此可见,信噪比在MUSIC算法中是一个很重要的因素。

2.3 阵元数引起的误差

阵元数作为模型的基本物理参数,在波达方向估计中也有一定的影响。取阵元数分别为8、9、10、11、12、13,信噪比为-10 dB,窄带源的角度、频率、阵元间距以及空间采样数均同2.2中所列,进行仿真,结果见图4,由图4可以看出,随着阵元的增加,算法的估计精度也逐步提高;但阵元的增加加大了计算量,效率降低。因此,在保证精度的前提下应选用合适天线阵阵元数。

除此之外,还有其它因素影响波达方向的估计精度,比如阵元间距[4](通常取工作波长的一半)、两信号源之间的角度差、阵列孔径等等。因此,在误差允许的情况下,应权衡各个因素以达到最佳效果。

3 结论

通过对MUSIC算法以及在波达方向估计中误差的分析与研究,从理论和仿真两方面证明将此法用于目标方位角估计的实用价值。尽管MUSIC算法在搜索过程中使用了所有的噪声特征向量导致了较大的计算量,但是作为空间谱估计理论体系中的标志性算法,MUSIC算法仍具有很高的分辨力、估计精度及稳定性,应用领域和前景十分广阔。民用上如海上救援,军用上如空中预警,智能制导武器系统,尤其是用于精确打击的机载雷达探测定位系统。另外,MUSIC算法的许多改良版[5],如基于相干的MUSIC算法、基于波速空间的MUSIC算法、求根MUSIC算法等等,也大大丰富了MUSIC算法的内容。

参考文献

[1]毕兰金,刘勇志.精确制导武器在现代战争中的应用与发展趋势.战术导弹技术,2004;23(6):54—57

[2]王永良,陈辉,彭应宁,等.空间谱估计理论与算法.北京:清华大学出版社,2004

[3]张光义,赵玉洁.相控阵雷达技术.北京:电子工业出版社,2006

[4]弋稳.雷达接收机技术.北京:电子工业出版社,2005

DOA估计 篇5

以MUSIC算法为代表的子空间类高分辨DOA估计算法是利用信号子空间与噪声子空间的正交性来实现DOA估计。此类算法通过对阵列接收数据协方差矩阵的特征分解来实现子空间的分离, 所以在实现参数估计高分辨力的同时也带来了巨大的运算量 (m×m矩阵的特征分解的运算量为O (m3) ) 。为了使此类算法减小运算量, 具备实时性, 已有的研究主要集中在以下2个方面:文献[1]给出的PM算法避开了矩阵的特征分解和奇异值分解, 有效地减小了算法的计算量 (运算量为O (m2n) ) ;文献[2]提出的子空间快速分解法采用Lanczos方法求取阵列接收数据协方差矩阵的正交基, 然后利用求得的正交基将协方差矩阵变换为一个三角对角阵, 再对该对角阵进行特征分解, 从而减小计算量 (运算量为O (m2n) , n<m, n为信源个数) 。本文通过详细分析PM算法机理, 利用阵列分割去噪原理对PM算法中提取噪声子空间基的数据矩阵进行了修改, 在不增加算法复杂度的前提下提高了PM算法的分辨力, 得到了较好的估计效果, 并通过计算机仿真试验验证了改进后算法的有效性。

1 PM算法

PM算法的前提假设:

① 窄带远场信号;② m>n, m为阵列的阵元数, n为信源数;③ 点信源假设, 且相互独立;④ 各阵元的噪声为加性白噪声, 且与信源不相关;⑤ 阵列流形矩阵A的前n行线性独立。

在上述假设条件下, PM算法可表述如下[1]:

阵列的流形矩阵是以导向矢量为列, 信源个数为列数, 阵元个数为行数构成的矩阵:

对于均匀线阵:

式 (1) ～式 (3) 中, d为阵元间隔, λ为信号波长, θi为第i个信源对应的角度, (·) T表示矩阵 (·) 的转置。

将流形矩阵作如下划分:

式中, A1为n×n的矩阵, A2为 (m-n) ×n的矩阵。因为A的前n行为线性独立且其秩为n, 所以A2为A1的线性变换:

式中, PH为 (m-n) ×n的矩阵, (·) H表示矩阵 (·) 的共轭转置; (·) -1表示矩阵 (·) 的求逆操作。构造m× (m-n) 的矩阵Q如下:

I为秩为m-n的单位阵, 则有:

A的列空间与信号子空间相同, 所以Q成了噪声子空间。获得矩阵Q的关键是求得矩阵P。假设阵列各阵元接收信号时没有噪声, 即阵列接收信号如下:

式 (9) 中, S为n×1的信号矢量;R1为m×n的矩阵;R2为m× (m-n) 的矩阵:

阵列实际的接收数据是含有噪声的, 则P的估计值$\widehat{\mathit{{\rm P}}}$为[1]:

从阵列接收的数据协方差矩阵中估计出矩阵$\widehat{\mathit{{\rm P}}}‚$结合式 (6) 可得矩阵Q, 从而可得噪声子空间, 再利用式 (11) 进行峰值搜索, 完成DOA估计:

2 消噪预处理的PM算法

由上述分析可知, 经典PM算法是在阵列接收数据协方差矩阵中提取噪声子空间的基。由于数据中噪声的存在, 从协方差矩阵中提取噪声子空间的算法为理论分析算法的近似 (式 (11) 为式 (10) 的最小二乘近似) , 从而影响了算法对噪声子空间逼近的精度。针对以上算法的不足, 本文提出通过将阵列分割为2个子阵来消除接收数据中的噪声, 在经过消噪预处理的数据基础上运用PM算法。这样就使从实际数据中提取噪声子空间的算法与理论分析的算法一致 (均利用式 (10) 进行求解) , 避免了由于噪声而产生的算法模糊, 从而提高算法的估计精度。消噪PM算法的表述如下:

消噪PM算法的前提假设:

① 窄带远场信号;② m>n, m为阵列的阵元数, n为信源数;③ 点信源假设, 且相互独立;④ 阵元的噪声为加性噪声与信源线性独立, 两子阵的噪声相互独立;⑤ 对阵列分割所得的两子阵的流行矩阵前n行线性独立。

在上述假设前提下, 将阵列分割为含有m-n个阵元的子阵1与含有n个阵元的子阵2, 二者的流形矩阵分别为B1=[b1 (θ1) , b1 (θ2) , …, b1 (θn) ], B2=[b2 (θ1) , b2 (θ2) , …, b2 (θn) ]。且

两子阵接收的数据为:

作2个子阵接收数据的相关:

由于任何矩阵C左乘满列秩矩阵或者右乘满行秩矩阵后, 矩阵C的秩保持不变。

式 (18) 可以变换为:

式中, B11为n×n的矩阵, 秩为n;B12为 (m-2n) ×n的矩阵;Rs为n×n的信号相关矩阵, 其为满秩;BH2为满行秩矩阵, 秩为n。由式 (19) 可得, 矩阵T1=B11BsBH2的秩与B14相同为n, T1为n×n的矩阵;T2为 (m-2n) ×n的矩阵。所以RY的前n行线性独立, 且其列空间与B1的列空间相同, 即为信号子空间, 则:$\mathit{{\rm P}}{}^{´\text{Η}}\mathit{{\rm T}}{}_1=\mathit{{\rm T}}{}_2‚\mathit{{\rm P}}{}^{´\text{Η}}=\mathit{{\rm T}}{}_2\mathit{{\rm T}}{}_1^{-1}‚{\mathit{Q}}^{\prime }=[\frac{{\mathit{{\rm P}}}^{\prime }}{-{{\rm I}}^{\prime }}]。$

I′为m-2n秩为的单位矩阵, 则有:

综合上述讨论可得消噪的PM算法的峰值搜索公式:

3 仿真验证

为了充分的说明本文对PM算法的改进效果, 对于均匀线阵分别运用PM算法与消噪PM算法进行DOA估计, 仿真过程与结果如下:

① 仿真前提:16阵元的均匀线阵;信源个数n=3;阵元间隔d=0.5λ;信噪比snr=5 dB;快拍数N=1 500;信源俯仰角θ=0, 1/4, 4/19, 单位为π;② 仿真结果:如图1所示。

由文献[1]可知PM算法随着快拍数的增大, 估计的偏差越来越小, 但计算量也越来越大。由图1可清楚地看出在大快拍数、低信噪比条件下, 消噪PM算法分辨力明显好于经典的PM算法。同时表明, 阵列接收数据协方差矩阵的噪声对PM算法的影响要大于由于阵列分割而产生的阵列孔径损失对PM算法的影响。

4结束语

本文基于文献[1]提出的快速子空间类DOA估计算法—PM算法, 利用阵列分割去噪原理, 对原算法进行修改, 在消去噪声的数据矩阵中提取噪声子空间的基, 实现对噪声子空间的估计。与原算法相比, 改进后的算法由于阵列的分割虽然带来了阵列孔径损失, 但在大快拍数条件下可以达到较好的去噪效果, 从而更好地逼近噪声子空间。相同条件下的DOA估计仿真试验表明改进的算法在大快拍数时较之原算法具有更好的分辨力, 同时算法的复杂度并未增加, 与原算法相比不影响算法的实时性。

参考文献

[1]MARCOS S, BENIDIR M.Source bearing estimation and sensor positioning with the propagator method[C].Signal Processing Alg, Arch, and Implementations, SPIE.1990:312-323.

[2]XUG, KAILATHT.Fast subspace decomposition[J].IEEE Trans.On Signal Processing, 1994, 42 (3) :535-551.

DOA估计 篇6

由于MUSIC算法在特定的条件下具有很高的分辨力、估计精度及稳定性。它是建立在以下假设基础上的:阵列形式为线性均匀阵;信号源数小于阵元的数目;处理器的噪声为加性高斯分布;空间信号为零均值平稳随机过程;信号源通常为窄带远场信号。

1、多重信号分类算法研究

1.1 经典MUSIC

假设有P个信号入射到阵列中, 则N元阵列的接收向量可表示为P个入射波与噪声的线性组合, 即

式中n (t) 为噪声向量, A为满列秩, s (t) =[s1 (t) , s 2 (t) …sp (t) ]T为入射信号向量,

α (θi) 是第i个信号波达方向的阵列导引向量。求得s (t) 协方差矩阵为:

Rxx是N*N的矩阵, 所以能分解为N个特征值和特征向量iλ, 则有:|Rxx-λiI|=0即|ARssAH- (λi-σn2) I|=0, ARssA的特征值为ηi, ηi=λi-σn2。ARssA为半正定, 秩为P。Rxx的特征值中有N-P个值等于σn2。如果将这N个特征值从大到小排列, 即λp=λp+1=…=λN-1=σn2。当最小特征值的重数K确定时, 信号的估计个数就可由下式确定.设特征值λi对应特征向量为νi, 与N-P个最小特征值对应的特征向满足:

可得:

设Un=[UP, UP+1, …UN-1]为噪声子空间, 则从式 (4) 可看出A正交于Un, 表明对应于波达方向的导引向量位于信号子空间内。因此, 对于多个入射信号的DOA可通过确定M U S I C空间谱的峰值作出估计, 其峰值为:

1.2 Eular MUSIC

在现代通信系统以及卫星系统中, BPSK和MASK信号是常用的信号, 这两种调制信号都具有信号为实信号的特点。根据信号的这种特点, 结合欧拉公式 (Euler's equation) , 将接收数据转换为实值的正弦和余弦数据, 然后将互不相等的新数据串联起来, 从而扩展了数据的维数, 相当于在原阵列阵元个数的基础上增添了同样个数的虚拟阵元。

对于B P S K和M A S K等实信号有s k (t) =sk* (t) , 令:

构造接收数据Yr (t) 为

重新构造的Yr (t) 是由Yc (t) 和Ys (t) 拼接而成的, 其维数为2N*M, 相当加倍了可利用的阵元个数。因为Ar相当于A去除虚数因子后的实部与虚部直接拼接而成, 且Ar的第N+1行元素均为0, 所以容易证明Ar对于2N-1>p其列满秩, 即提出算法最多可以处理2 (N-1) 个信号。此时的数据已经转化为实值矩阵, 所以下面的计算均是在实数的基础上进行的。

对Yr (t) 的协方差矩阵进行特征值分解, 得到表示信号对应的特征值与特征向量Us, 噪声对应的特征值与特征向量Us, 从而对Eular_MUSIC空间谱的峰值作出估计:

其中:

2、matlab仿真结果

取等间距线阵, 对MUSIC算法进行仿真, 仿真条件为:阵元间距1/2波长, 快拍数1024, 阵元数目6, SNR=30。图1中3个信号源的初始方位角为30°、50°、70°, 然后减小信源来角间隔, 信号源的方位角变为30°、50°、58°, 降低信噪比到SNR=10, 减小阵元数目m=6。由以上仿真结果可以得出, 在阵元数目减少、低信噪或角度相隔比较近的情况下, 经典MUSIC算法的准确度下降。图2是在3个信源来角方向为30°、50°、58°, SNR=10, m=4和8个信号源, 来波方向分别从30°到135°, 角度间隔为15°, 阵元数目为6, SNR=20d B的仿真条件采用Eular_MUSIC算法进行仿真。

3、结语

MUSIC算法是一种信号参数估计算法, 给出的信息包括入射信号的数目、各个方向的波达方向、强度以及入射信号和噪声间互相关等。采用经典的MUSIC算法对于高信噪比, 角度相隔较远的信号的DOA估计较为理想, 但是对于小信噪比和角度相隔比较近的信号时, MUSIC算法并不能准确地估计信号的DOA。而通过与欧拉公式的结合, 扩展了数据的维数, 降低了运算的复杂性, 提高了经典MUSIC算法的精度高和分辨力, 并且可以处理比阵元个数多的信号的能力。

摘要：在阵列信号处理领域, 波达方向 (DOA) 估计一直是研究的重点之一。在波达方向 (DOA) 估计中, 利用多重分类算法 (MUSIC) 对来波方向进行估计是最常用的方法。本文概述了经典MUSIC算法, 针对现代通信中常用的BPSK和MASK信号都是实信号的特点, 结合Eular公式对MUSIC算法进行了改进, 使用matlab进行了仿真及对比。

关键词：波达方向估计,多重分类算法,matlab仿真

参考文献

[1]YEH C C.Bearing estimations with mutual coupling present[J].IEEE Trans on Antennas and Propagation, 1989, 37 (10) ;1332-1335.

[2]STOICA P, NEHORAI A.MUSIC, maximum likelihood, and Cramer-Raobound[J].IEEE Trans on ASSP, May 1989, 37 (5) :720-741.

[3]何子述, 黄振兴.修正MUSIC算法对相关信号源的DOA估计性能.通信学报, 2000 (10) :15-17.

DOA估计 篇7

关键词：MIMO雷达,稀疏信号恢复,迭代适应法 (IAA) ,到达角 (DOA) 估计,降维

最近几年来, MIMO雷达受到国内外研究者的广泛关注[1]。相比传统相控阵雷达, MIMO雷达发射相互正交的信号, 在接收端进行匹配滤波, 形成独立的信号通道, 因此拥有更高的系统自由度, 更高的空间分辨力和参数识别能力[1—3]。

MIMO雷达按阵元配置方式, 可以分为单基地、双基地、统计MIMO雷达三种。方位角估计是MIMO雷达参数估计中的一个关键问题。文献[4]研究了基于Capon算法的MIMO雷达目标角度估计算法。文献[5]利用接收阵列与发射阵列的旋转不变性, 运用ESPRIT算法对到达角 (DOA) 与离开角 (DOD) 进行估计;文献[6, 7]应用Parafac模型实现了双基地MIMO雷达DOA与DOD的联合估计;文献[8]提出一种降维Capon算法, 减少了运算量。文献[9]研究了波形未知情形下, 单基地MIMO雷达的波达方向估计。上述算法均是将传统算法扩展到MIMO雷达方位角估计问题中, 在低快拍情形下, 这些算法均不能很好地分辩。近几年提出的基于稀疏信号表示的DOA估计算法, 在低快拍情形下也具有较高分辩力[10,11], 但这些只是针对传统相控阵列雷达, 针对MIMO雷达体制, 利用稀疏信号表示的DOA估计算法并未深入研究。提出一种基于稀疏信号表示的单基地MIMO雷达DOA估计方法, 该算法利用目标在空域的稀疏性, 建立了MIMO雷达回波信号的稀疏模型, 同时结合接收导向向量与发射导向向量Kronecker积的特性, 对阵列流形过完备字典进行降维预处理, 然后构造加权代价函数, 以最小化该代价函数的方式, 通过迭代获得波达方向估计。仿真实验验证了所提算法的有效性, 并分析了信噪比、快拍数与阵元个数对算法性能的影响。

1 信号模型

考虑一单基地MIMO雷达, 收发同置。设发射阵元个数为M, 阵元间距为dt, 接收阵元个数为N, 阵元间距为dr。各个发射阵元发射波长为λ的窄带正交信号, 发射阵与接收阵均为均匀线阵, 且阵元间的间距均大于半波长。设远场区域存在P个不相关目标, 发射信号经过目标反射, 通过匹配滤波器后的输出信号为

式 (1) 中, θi为第i个目标的到达角, ar (θi) 和αt (θi) 分别为接收导向向量、发射导向向量。s (t) 表示目标向量

式 (2) 中, sp (t) 为第p个目标的多普勒频率, βp为幅度, n (t) 表示高斯白噪声向量, 均值为0, 协方差矩阵Rn=σ2IMN, σ2为噪声方差, IMN表示MN维单位阵, 表示Kronecker积, (·) T和 (·) H表示矩阵转置和共轭转置。

1.1 接收信号稀疏模型

远场区域内存在P个目标, 在空域内可以认为其具备稀疏性。因此, 原问题可以转化为稀疏信号恢复问题。将空域从-90°到90°等分成K个区域, PK, 若K足够大, 则离散观察点十分接近真实的角度, 接收信号y (t) 可近似表示为

式 (3) 中, B∈RMN×K是离散阵列流形, x (t) ∈RK×1是与B中列向量对应系数构成的系数向量。

式 (4) 中, br (θ1) 与bt (θ1) 分别为接收方向向量与发射方向向量。角度是等分观察区域后, 离散观察点的角度。若方位k处存在目标, 则向量x (t) 中相对应元素不为0;反之, 则向量中相对应元素为0。因此, 从y (t) 中求出一个稀疏度为P的解, 便可求得目标到达角。

1.2 降维预处理

发射方向向量与接收方向向量的Kronecker乘积形成了B中的每一列Bi∈RMN×1, 但是Bi中出现了许多重复元素, 因此可以通过矩阵变换进行降维, 消除重复元素。定义降维矩阵T∈RMN× (M+N-1) 、降维导向向量b∈R (M+N-1) ×1和W∈R (M+N-1) × (M+N-1) 。

式中, 表示整除, mod () 表示取余。将式 (3) 左乘W-12TH, 得到降维后的接收信号向量

n~ (t) 的协方差矩阵Rn~ (t) 为

由式 (12) 可以看出, 降维预处理并没有引入其他的误差, 各阵元噪声的统计特性也没有改变。考虑式 (9) , 阵列流形C左乘了一对角矩阵, 这相当于对C的每一行进行加权, 但稀疏性并未改变, 因此, 仍然可以采用稀疏恢复算法来求得方向角。

2 算法描述

通过前面的降维预处理, 接收信号向量已由原来的MN×1维变成 (M+N-1) ×1维, 扫描阵列流形也由原来的MN×K变成了 (M+N-1) ×K维。为了从中恢复出稀疏信号, 现提出一种非参数的迭代算法。考虑一加权最小方差代价函数

式 (13) 中, y~ (l) 为降维后, 第l次快拍接收信号向量, xk (l) 为系数向量x (l) 的第k个元素。定义矩

阵C~

式 (14) 中, 为矩阵的第k列向量。假设目标之间相互独立, 定义对角阵P∈RK×K,

求代价函数f[xk (l) ]的最小值, 得

根据矩阵求逆公式Q-1≈Rs-1, 将其代入式 (18) 得

算法步骤如下:

(1) 按式 (5) 、式 (7) 、式 (8) 分别构造矩阵T、C和W。

(2) 按式 (10) 、式 (14) 求得。

(3) 计算pk初值

(5) 对每个k值计算:

重复步骤 (4) 、 (5) , 直到收敛。

算法迭代15次之后, xk^的值几乎不再变化, 将该值作为xk的估计值。

3 仿真结果

本节通过蒙特卡罗实验验证文中所提算法的有效性。定义DOA估计均方根误差 (RMSE) 。

式 (20) 中, θp为第p个目标的真实到达角, θp, n为第p个目标第n次实验的到达角估计值。

考虑3个等功率非相关目标s1、s2、s3其方位角分别是θ1=-10o、θ2=0o、θ3=20o。空间角度扫描范围-90o≤θ≤90o, 空间采样点数K=181, 迭代次数I=15。

3.1 低信噪比、低快拍下算法性能

信噪比SNR=5 d B, 发射阵元数M=6, 接收阵元数N=6, 快拍数L=2。分别采用降维Capon[8]、IAA-R[12]和本文所提算法进行10次蒙特卡罗仿真实验, 所得结果分别由图1、图2、图3所示。由图1可知, 降维Capon算法在低快拍情形下已不能分辩。由图2、图3可知, 尽管降低了扫描阵列流形的维数, 但是仍能很好地分辩出目标到达角。

信噪比-10 d B≤SNR≤30 d B, 发射阵元数M=5, 接收阵元数N=5, 快拍数L=2。分别使用IAA-R[12]和本文算法进行100次蒙特卡罗仿真实验, 所得RMSE如图4所示。图4中, 本文算法与IAA-R算法RMSE变化曲线很相近, 表明了降维后的接收数据和扫描阵列流形对算法性能几乎没有影响。

3.2 信噪比与快拍数对算法性能的影响

信噪比-10 d B≤SNR≤30 d B, 发射阵元数M=6, 接收阵元数N=6, 快拍数L=2、8、16。进行100次蒙特卡罗仿真实验。由图5可知, 随着快拍数与信噪比增大, 估计值越接近真实值。

3.3 接收阵元数、发射阵元数对算法性能影响

信噪比0 d B≤SNR≤30 d B, 发射、接收阵元数{M, N}分别取5、6、7, 快拍数L=2, 进行100次蒙特卡罗实验。由图6可知, 随着快拍数与发射阵元和接收阵元数增大, 估计值越接近真实值。

4 结论

利用目标在空域的稀疏性, 将目标方向角估计转化成稀疏信号恢复问题进行求解。根据MIMO雷达接收导向向量与发射导向向量kronecker乘积的特性, 首先对接收数据进行降维预处理, 构造新的阵列流形, 然后再利用稀疏恢复算法完成DOA估计。本文算法在低快拍数下能很好地分辨。同传统稀疏信号恢复算法相比, 该算法运算量要低于一般稀疏信号恢复算法。仿真实验验证了新算法的有效性, 具有一定的工程参考价值。

参考文献

[1] Fishler E, Haimovich A, Blum B, et al.MIMO radar:an idea whose time has come.Proc of IEEE Radar Conf, 2004:71—78

[2] Fishler E, Haimovich A, Blum R, et al.Performance of MIMO radar systems:advantages of angular diversity.Proc of 38th Asilomar Conf Signal Syst, Comput, Pacific Grove, CA, 2004;1:305—309

[3] Haimovich A M, Blum R S, Lenard J, et al.MIMO radar with widely separated antennas.IEEE Signal Processing Magazine, 2008;25 (1) :116—129

[4] Xia Wei, He Zishu.Multiple-target localization and estimation of MIMO radars using capon and APES tech-niques.Proc of IEEE Radar Conference, Rome, 2008:1—6

[5] Duofang C, Baixiao C, Guodong Q.Angle estimation using ESPRIT in MIMO radar.Electron Lett, 2008;44 (12) :770—771

[6] Nion D, Sidiropoulos N D.A PARAFAC-based techn-ique for detection and localization of multiple targets in a MIMO radar system.2009IEEE Int Conf on Acoustics Speech and Signal Processing, 2009:2077—2080

[7] Zhang X, Xu Z, Xu L, et al.Trilinear decomposition-based transmit angle and receive angle estimation for multiple-input multiple-output radar.IET Radar Sonar Navig, 2011:5 (6) :626—631

[8] Zhang X, Huang Y, Chen C, et al.Reduced-complexity Capon for direction of arrival estimation in a monostatic multiple-input multipleoutput radar.IET Radar, Sonar and Navigation, 2012;6 (8) :796—801

[9]陈晨, 张小飞.波形相关矩阵未知情况下单基地MIMO雷达中一种改进MUSIC的DOA估计算法.电子与信息学报, 2012;34(12) :2966—2971

[10] Donoho D L, Elad M, Temlyakvo V.Stable recovery of sparse overcomplete representations in the presence of noise.IEEE Trans Inf Theory, 2006;52 (1) :6—8

[11] Chunmei Z, Zhongke Y, Xiangdong C, et al.Signal overcomplete representation and sparse decomposition based on redundant dictionaries.Chinese Sci Bull, 2005;50 (23) :2672—2677