版权声明:本文为博主原创文章,未经博主允许不得转载。 https://blog.csdn.net/qq_23947237/article/details/82502848
二维 DOA 估计算法
0. 前言
阵列信号处理基础本质上属于参数估计问题,和信道估计知识基本上别无二致。末学在这里整理了阵列信号处理的基础知识,包括公式推导,以及代码。一方面为了节省同行人士寻找资料和整理吸收的时间,开方便之门。另一方面为了和大家多多交流这方面的知识,寻找研究灵感。
如果有任何问题或者有相关的 MATLAB 代码,本着交流的态度请分享到我的邮箱:[email protected]。
愿以此功德,庄严佛净土。上报四重恩,下济三途苦。
若有见闻者,悉发菩提心。尽此一报身,同生极乐国。
南无大智文殊菩萨。
在实际通信中,用户位于三维空间,因此用户的 DOA 需要由二维角度来刻画。如果只是估计出用户的一维 DOA,那只是可以确定用户存在于一个平面,还不能精确的定向。因此利用二维 DOA 比只利用一维 DOA 方位角信息能更有效的开发信道的冗余信息,更好的描述信号特征。
二维 DOA 估计算法可以分为两大类:
利用天线阵列来完成空间二维DOA估计;
利用空时二维处理的方法来完成空间二维DOA估计。
1. 空间阵列二维 DOA 估计
利用空间天线阵列来完成空间二维 DOA 估计的算法有两大类:
而常用的空间天线阵列模型有:平行线阵和 L 型天线阵列。 基于 ESPRIT 类算法,下面提出一种解决参数配对问题的二维 DOA 估计算法,其不仅解决了空间方位角和俯仰角参数配对问题,而且相对于其他二维来波方向估计算法,具有较少的计算复杂度。该二维 DOA 估计算法是基于平行天线的阵列模型。
快速二维 DOA 估计
系统模型如图所示,天线阵列为双平行天线阵列结构。
x
轴阵元间距为
dx
,两平行天线阵列间距为
dy
。第一个均匀线阵由
(N+1)
个阵元组成,第二个均匀线阵由
N
个阵元组成。假设位于远场的具有相同波长的
P
个统计独立的窄带信号从二维方向
(rk,φk)
入射到阵列。
其中
(rk,φk)
分别为方位角和俯仰角。方位角表示在
xoy
面上的投影与
x
轴夹角,俯仰角为投影与入射信号之间的夹角。
为处理方便,下面用
(αk,βk)
表示
(rk,φk)
,其中对应关系有
cosαkcosβk=cosφk⋅cosrk=cosφk⋅sinrk
从双平行线阵列结构构造 3 个子阵,其中第一个子阵为第一个均匀线阵的前
N
个阵元组成,第二个子阵为第一个均匀线阵的后
N
个阵元组成,第 3 个子阵为第二个均匀线阵
N
个阵元组成。分别用
X,Y,Z
表示子阵 1,2,3 接收的来波信号矩阵。子阵 1,2,3 的输出噪声为
n1,n2,n3
,假设它们为零均值且方差为
σ2
的高斯白噪声,且与信号矩阵不相关。需要注意的是,
n1
和
n2
相关,它们和
n3
不相关。
n1n2n3=[n11(t),n12(t),⋯,n1N(t)]T=[n21(t),n22(t),⋯,n2N(t)]T=[n12(t),n13(t),⋯,n1,N+1(t)]T=[n31(t),n32(t),⋯,n3N(t)]T
假设来波信号
si(t), (i=1,2,⋯,P)
与
x
轴夹角为
αi
,与
y
轴夹角为
βi
,则子阵所接收的来波信号矩阵为
XYZ=AS+n1=AΦ1S+n2=AΦ2S+n3
其中,
X=[x1(t),x2(t),⋯,xN(t)]T
,
Y,Z
类似。注意两个阵元之间的波程差如下图,同一个信号到达相邻阵元存在一个波程差
τ=dsinθc
,波程差导致了相位差
ψ=e−jwτ
。
ψ=e−jwτ=e−jwdsinθc=e−j(2πf)dsinθλf=e−j2πdsinθλ
于是根据下图在一条线上的阵列几何关系,有
cτ=dycosβ
成立。
再考虑以原点为参考点,另外一个阵元位于
(dx,dy)
,则有关系
cτ=dxcosα+dycosβ
成立。
其中,
a(αi)=[1,e−j2πdsinαi/λ,⋯,e−j(M−1)2πdsinαi/λ]T
。因此有
SAΦ1Φ2=[s1(t),s2(t),⋯,sP(t)]T=[a(α1),a(α2),⋯,a(αP)]=diag(exp(j2πdxcosα1/λ),⋯,exp(j2πdxcosαP/λ))T=diag(exp(j2πdycosβ1/λ),⋯,exp(j2πdycosβP/λ))T
利用上述数据模型,可以得到下列数据矩阵
RXXRXYRXZRYYRYZRZZ=E(XXH)=ARSAH+σ2I=E(XYH)=ARSΦH1AH+σ2Q=ARSΦH1AH+σ2[0IN−100]=E(XZH)=ARSΦH2AH=E(YYH)=AΦ1RSΦH1AH+σ2I=E(YZH)=AΦ1RSΦH2AH=E(ZZH)=AΦ2RSΦH2AH+σ2I
计算数据矩阵
RXX
的特征分解,由特征值大小排列估计噪声功率
σ2
。接下来我们令
CXXCXYCXZCYYCYZCZZ=RXX−σ2I=ARSAH=RXY−σ2Q=ARSΦH1AH=RXZ=ARSΦH2AH=RYY−σ2I=AΦ1RSΦH1AH=RYZ=AΦ1RSΦH2AH=RZZ−σ2I=AΦ2RSΦH2AH
利用上面的形式构造矩阵
C
,得到
C=⎡⎣⎢CXXCHXYCHXZCXYCYYCHYZCXZCYZCZZ⎤⎦⎥=⎡⎣⎢⎢ARSAHAΦ1RSAHAΦ2RSAHARSΦH1AHAΦ1RSΦH1AHAΦ2RSΦH1AHARSΦH2AHAΦ1RSΦH2AHAΦ2RSΦH2AH⎤⎦⎥⎥=⎡⎣⎢AAΦ1AΦ2⎤⎦⎥RS[AHΦH1AHΦH2AH]=BRSBH
如果
P
个信号源来波方向
α,β
不同时接近,则有:
Rank(B)=P
。则可以对矩阵
C
进行特征值分解得到子空间
US
,显然信号子空间
US
和
B
张成的子空间一致,因此存在唯一的可逆矩阵
T
,满足下面的关系:
US=⎡⎣⎢U1U2U3⎤⎦⎥=⎡⎣⎢AAΦ1AΦ2⎤⎦⎥T=BT
根据以上关系可知,
U2=U1T−1Φ1T=U1Ψ1
,
U3=U1T−1Φ2T=U1Ψ2
。我们定义
Ψ1Ψ2=T−1Φ1T=U†1U2=T−1Φ2T=U†1U3
其中,
U†1=(UH1U1)−1UH1
。显然对
Ψ1,Ψ2
进行特征分解可得到
Φ1,Φ2
的特征值,就可得到入射角度
(αi,βi)
。由于
Ψ1,Ψ2
单独处理,所以存在角度配对问题。我们做以下处理
Ψ1Ψ2=T−11Φ1T1=T−12Φ2T2
这里
T1,T2
是由矩阵
T
的列组成的。通过比较特征向量矩阵
T1,T2
的列向量可以决定二维 DOA 估计
(αi,βi)
的参数配对问题。在实际计算中 2 个特征值分解是独立进行的,可能特征向量的排列顺序是不同的,但同一信号的特征向量是
强相关的。所以可构造排序矩阵
G
,对特征值的顺序进行调整以实现参数配对。计算
G=TH2T1
根据
G
中每一列元素中绝对值最大的一个元素的矩阵坐标来调整
Φ
中对角元素的顺序,实现特征向量的对应关系。进而可以求取信号的二维到达角。
matlab 仿真
clear variables;
close all; clc;
degrad=pi/180;
radeg=180/pi;
N=5;
P=3;
snr=10;
dd=0.5;
k=1024;
d=0:dd:(N-1)*dd;
a = [30 50 70];
b = [20 40 60];
figure;
plot(a,b,'r*');
hold on;
twpi=2*pi;
A1 = [exp(-1i*twpi*d.'*(sin(a*degrad)))];
ph1= [exp(1i*twpi*dd.'*(cos(a*degrad)))];
ph2= [exp(1i*twpi*dd.'*(cos(b*degrad)))];
A2 = A1*diag(ph1);
A3 = A1*diag(ph2);
N0 = sqrt(1/2)*(randn(N+1,k)+1i*randn(N+1,k));
N1 = N0(1:N,:);
N2 = N0(2:N+1,:);
N3 = sqrt(1/2)*(randn(N,k)+1i*randn(N,k));
S = randn(P,k);
X = A1*S+N1;
Y = A2*S+N2;
Z = A3*S+N3;
Rxx=X*X'/k;
Rxy=X*Y'/k;
Rxz=X*Z'/k;
Ryy=Y*Y'/k;
Ryz=Y*Z'/k;
Rzz=Z*Z'/k;
% 估计噪声功率
[V,D]=eig(Rxx); % 特征值 对角阵D和特征向量V
D=diag(D)'; % 把特征值对角阵变换成行矩阵
[D,~]=sort(D);
sigma2 = mean(D(1:N-P));
Cxx=Rxx-sigma2*eye(size(Rxx));
Cxy=Rxy-sigma2*diag([ones(1,N-1)],-1);
Cxz=Rxz;
Cyy=Ryy-sigma2*eye(size(Ryy));
Cyz=Ryz;
Czz=Rzz-sigma2*eye(size(Rzz));
C = [Cxx, Cxy, Cxz;
Cxy',Cyy, Cyz;
Cxz',Cyz',Czz];
% 估计噪声功率
[Ug,Sg]=eig(C); % 特征值分解
Sg=diag(Sg)'; % 行向量
[Sg,I]=sort(Sg);
Sg=fliplr(Sg);
Ug=fliplr(Ug(:,I));
U=Ug(:,1:P);
U1 = U(1:N,:);
U2 = U(N+1:2*N,:);
U3 = U(2*N+1:3*N,:);
Psi1 = (U1'*U1)\U1'*U2;
Psi2 = (U1'*U1)\U1'*U3;
[T1,S1]=eig(Psi1);
[T2,S2]=eig(Psi2);
G=T2'*T1;
G = abs(G);
[max_g,index]=max(G,[],2);
[Row,Col]=size(G);
for i=1:Row
for j=1:Col
if j == index(i)
G(i,j)=1;
else
G(i,j)=0;
end
end
end
T1=T1*inv(G);
S1=inv(T1)*Psi1*T1;
S1=diag(S1)'; % 把特征值对角阵变换成行矩阵
S2=diag(S2)';
an = acos(abs(angle(S1))/(twpi*dd))*radeg;
bn = acos(abs(angle(S2))/(twpi*dd))*radeg;
plot(an,bn,'.')
axis([10 90 10 90]);
set(gca, 'XTick',[10:10:90]);
grid on;
仿真估计结果如图所示:
2. 空时二维 DOA 估计
一般的二维 DOA 特征结构方法都会遇到二维谱峰搜索、非线性优化、分维处理及参数配对等难题,而且二维 DOA 估计必须借助平面阵或立体阵,阵元数较多。因此计算量和存贮量变大,降低了算法的实用性。另外,还需要通过特殊的阵列结构才能实现二维 DOA 估计,增加了阵列孔径损失,阵元利用率降低,阵元间互耦加剧,设备的成本、体积和复杂性增加,阵元的校正也十分困难。
因此,在利用特征结构法进行空域处理的同时引入适当的时域处理,更加充分地利用信号中的信息。这样可以降低对阵列结构的约束,提高算法的抗噪能力。将同时利用信号空域信息,以及时域信息的方法称为空时二维 DOA 估计法。该算法不需要双平行线阵,也不存在角度兼并等问题,对噪声不敏感,并可推广到任意形状阵列且阵元大部分未校准的二维 DOA 估计问题。
另外,在无线通信信道中,同信道干扰和码间干扰是限制天线容量和数据率的主要因素,利用时域一维处理或者利用空域一维处理都存在很多局限性,因此选用空时二维处理的方法则能取到较好的效果。典型的空时接收器是接收调制解调器,它同时对所有天线工作,并在时间和空间的两个域处理接收的信号。下面首先对空时二维处理方法数据进行系统建模。
上图所示均匀线阵模式,在
M
个传感器的输出端加上
L
个延迟单元,每个延迟单元的时延为
τ
。其中
l
为延迟单元的序号,当
l=0
时表示第
m
个阵元的信号。则第
m
个阵元连同其后的第
l
个延迟单元上的信号可以表示为
xml(t)=∑i=1Psi(t+(m−1)dsinθic−lτ)+nml(t)
考虑所有的阵元,写成矩阵的形式为
Xl=AlS+nl
其中
Xl
为第
l
个延迟单元上的信号矩阵,
nl
为零均值、方差为
σ2
统计独立的高斯白噪声,且与信号矩阵不相关。其中
Al=A0Φl
。
当不存在延迟时,有
X0=A0S+n0
这里
A0=[a(θ1),⋯,a(θp)]
,其中
a(θp)=[1,e−jw(dsinθp/c),⋯,e−jw((M−1)dsinθp/c)]T
。
可知,通过时域的延时,可以得到虚拟阵元形式,使得算法能处理的入射信源个数大于等于阵元个数。由于实际延迟单元的延迟时间可知,所以可以结合时间域上所得到虚拟天线阵列和空间域上的天线阵列进行空时二维 DOA 估计。
二维虚拟 ESPRIT 算法
数据模型如下图所示的天线接收阵列,它由三个子线性阵列组成。每个子阵列有
N+1
个阵元,子阵列相互平行,且完全相同,子阵列里每个阵元的间距均为
dx
,子阵列之间的距离为
dy
。
定义
Y0
为
1~N
阵元输出序列,
Y1
和
Y2
为
2~(N+1)
阵元输出序列。信号为窄带远场平面波时,假设来波信号
si(t),(i=1,2,⋯,p)
与
x
轴夹角为
αi
,与
y
轴夹角为
βi
。子阵
Y0,Y1,Y2
的输出噪声为
n0,n1,n2
,假设它们为零均值的方差为
σ2
统计独立的高斯白噪声,且与信号矩阵不相关。则有
Y0Y1Y2Y0Y1Y2SA=AS+n0=AΦ1S+n1=AΦ2S+n2=[y01(t),y01(t),⋯,y0N(t)]T=[y11(t),y11(t),⋯,y1N(t)]T=[y21(t),y21(t),⋯,y2N(t)]T=[s1(t),s2(t),⋯,sp(t)]T=[a(α1),a(α2),⋯,a(αp)]
由上面定义可以看到,
Y1,Y2
并不是简单的平行于
Y0
,而是滞后
Y0
于一个子阵列的间隔单元
dx
。则有
Φ1Φ2uivi=diag(ej(u1+v1),ej(u2+v2),⋯,ej(up+vp))=diag(ej(u1−v1),ej(u2−v2),⋯,ej(up−vp))=2πdxcosαiλ=2πdycosβiλ
利用上述数据模型,令
Y⟹Y=Y1+Y22=A(Φ1+Φ2)S+(n1+n2)2=AΦS+(n1+n2)2=AΦS+N2
复指数
ejx=cosx+jsinx
。利用和差化积公式
sinx+sinycosx+cosy12(ejx+ejy)=2sinx+y2cosx−y2=2cosx+y2cosx−y2=cosx−y2[cosx+y2+jsinx+y2]=cosx−y2exp(jx+y2)
可得,
Φ=diag(cosv1eju1,cosv2eju2,⋯,cosvpejup)
。
由
Φ
的形式可知,利用以前的 ESPRIT 算法,可以求出
Φ
的模值和相位,就能求出相应的方位角和俯仰角,故能达到计算量小、估计精度高、对子阵列要求不高的优点。为了叙述方便,令
X=Y0=AS+n0⟹X=AS+N1
于是可得
RXXRXY⟹CXXCXY⟹C=CXX−γCXY=E(X(n)XH(n))=ARSAH+σ2I=E(X(n)YH(n))=ARSΦHAH=RXX−σ2I=ARSAH=RXY=ARSΦHAH=ARSAH−ARS(γΦH)AH=ARS(I−γΦH)AH
若矩阵
A,B
,满足
Au=ξBu
,则称
ξ
是
{A,B}
的广义特征值,
u
是广义特征向量。
(ξ,u)
是广义特征对。
当且仅当
ξ
是
{A,B}
的广义特征值时,
A−ξB
奇异。广义特征值是广义特征多项式
∣∣A−ξB∣∣=0
的根。
假设信号角度不相关,所以
A,RS
的秩都为
p
,故奇异值为
∣∣I−γΦH∣∣=0
的根。 注意在这里求多个
γ
值组成的对角阵的对角元素,且
ΦH=diag(cosv1e−ju1,⋯,cosvpe−jup)
。于是得到
γ=diag(1cosv1eju1,⋯,1cosvpejup)
到目前为止,建模为
XY=AS+N1=AΦS+N2
通过前面的分析,可以知道我们只要获得
Φ
或者
γ
,就可以求出角度。一种方法是利用求
{CXX,CXY}
的广义特征值得到
γ
。另一种方法是按照 TLS-ESPRIT 的思路求
Φ
,可以进行如下操作:
Z=[XY]=[AS+N1AΦS+N2]=[AAΦ]S+[N1N2]=A¯S+N
然后求
RZ=E(ZZH)=1K∑Kt=1Z(t)Z(t)H
,进行分解得到信号子空间
US
。根据
US=[UXUY]=[AAΦ]T
记
UST=[UX|UY]
,然后有
UHSTUST=EΣEH=[ESEN][ΛSΛN][EHSEHN]
最后得到
EENΨTLS=[E11E21E12E22]=[E12E22]=−E12E−122
通过
ΨTLS
最终得到
Φ
。利用
Φ=diag(cosv1eju1,cosv2eju2,⋯,cosvpejup)
求出
u
和
v
。
u = angle(phi);
cos_v = abs(phi);
v = acos(cos_v);
a = acos(lamda*u/(2*pi*dx));
b = acos(lamda*v/(2*pi*dy));
虚拟波达方向矩阵法
改进的二维虚拟 ESPRIT 算法:虚拟波达方向矩阵法中,其数据模型与二维虚拟 ESPRIT 算法一致。根据二维虚拟 ESPRIT 算法,求出
Φ
即可。虚拟波达方向矩阵法与虚拟 ESPRIT 算法的不同就在于求
Φ
的过程不同。将波达矩阵法的求解思想用在求解
Φ
上,从而得到虚拟波达方向矩阵法。其相对于二维虚拟ESPRIT算法,减少了计算复杂度,更适合于低信噪比,多信源的通信环境。
系统模型为
XY=AS+N1=Y0=AS+n0=AΦS+N2=AΦS+(n1+n2)2
协方差以及估计噪声功率等相关操作
RXXRYX⟹CXX=E(X(n)XH(n))=ARSAH+σ2I=E(Y(n)XH(n))=AΦRSAH=RXX−σ2I=ARSAH=VΣVH
式中
Σ=[μ1,⋯,μN]
是特征值,而
{V1,⋯,VN}
是相应的特征向量,由前面 MUSIC 算法的信号子空间和噪声子空间的分析可知:
-
CXX
的特征值
μ1⩾μ2⩾⋯⩾μp>μp+1=⋯=λN=0
。
- 特征向量组
V=[Vs,Vn]
。
span{Vp+1,Vp+2,⋯,VN}⊥span{A}
,因此
AHVn=0
。
于是有
AHVnVHnAAHVVHA=0=AH(VsVHs)A+AH(VnVHn)A=AHVsVHsA=AHA
有了上述性质,推出波达方向矩阵法,记波达方向矩阵
R
:
R=RYX∗C†XX
其中
C†XX=VΣ−1VH
。
有以下推导
CXXRSAHRYXRA=ARSAH=A†CXX=(AHA)−1AHCXX=AΦRSAH=AΦ(AHA)−1AHCXX=RYX∗C†XXA=AΦ(AHA)−1AHCXX∗C†XXA=AΦ
可得结论:
RA=AΦ
。对波达方向矩阵
R
进行特征分解得到
Φ
,根据空间方位夹角与
Φ
之间对应关系得出来波方向角
(α,β)
。
虚拟波达方向矩阵法与二维虚拟 ESPRIT 算法的不同点在于,虚拟波达方向矩阵法只需要一次
(N×N)
的互协方差运算,一次
(N×N)
互协方差的求逆,以及一次
(N×N)
的特征值分解。而二维虚拟 ESPRIT 算法需要一次
(2N×2N)
阶,一次
(2p×2p)
阶和一次
(p×p)
阶的特征值分解,所以当信源数较大时使用虚拟波达方向矩阵法可以大大节省算法的运算量。