DOA估计算法

Posted 2023-04-16

tags:

篇首语：本文由小常识网(cha138.com)小编为大家整理，主要介绍了DOA估计算法相关的知识，希望对你有一定的参考价值。

参考技术A 学号：20000300055

姓名：王铎澎

嵌牛导读：文章对DOA算法进行了简单的介绍。

嵌牛正文：https://blog.csdn.net/zhangziju/article/details/100730081?ops_request_misc=%257B%2522request%255Fid%2522%253A%2522160689878119725222413438%2522%252C%2522scm%2522%253A%252220140713.130102334..%2522%257D&request_id=160689878119725222413438&biz_id=0&utm_medium=distribute.pc_search_result.none-task-blog-2~all~baidu_landing_v2~default-1-100730081.pc_first_rank_v2_rank_v28&utm_term=Musicsuanfa&spm=1018.2118.3001.4449

DOA估计算法

DOA（Direction Of Arrival）波达方向定位技术主要有ARMA谱分析、最大似然法、熵谱分析法和特征分解法，特征分解法主要有MUSIC算法、ESPRIT算法WSF算法等。

MUSIC (Multiple Signal Classification)算法，即多信号分类算法，由Schmidt等人于1979年提出。MUSIC算法是一种基于子空间分解的算法，它利用信号子空间和噪声子空间的正交性，构建空间谱函数，通过谱峰搜索，估计信号的参数。对于声源定位来说，需要估计信号的DOA。MUSIC算法对DOA的估计有很高的分辨率，且对麦克风阵列的形状没有特殊要求，因此应用十分广泛。

运用矩阵的定义，可得到更为简洁的表达式：

X = A S + N X=AS+NX=AS+N

式中

X = [ x 1 ( t ) , x 2 ( t ) , . . . x M ( t ) ] T X=[x_1(t),x_2(t),...x_M(t)]^TX=[x1(t),x2(t),...xM(t)]T，

S = [ S 1 ( t ) , S 2 ( t ) , . . . S D ( t ) ] T S=[S_1(t),S_2(t),...S_D(t)]^TS=[S1(t),S2(t),...SD(t)]T，

A = [ a ( θ 1 ) , a ( θ 2 ) , . . . a ( θ D ) ] T A=[a(\theta_1),a(\theta_2),...a(\theta_D)]^TA=[a(θ1),a(θ2),...a(θD)]T，

N = [ n 1 ( t ) , n 2 ( t ) , . . . n M ( t ) ] T N=[n_1(t),n_2(t),...n_M(t)]^TN=[n1(t),n2(t),...nM(t)]T。

X XX为阵元的输出，A AA为方向响应向量，S SS是入射信号，N NN表示阵列噪声。

其中 φ k = 2 π d λ s i n θ k \varphi_k=\frac2\pi d\lambdasin\theta_kφk=λ2πdsinθk有

A = [ a ( θ 1 ) , a ( θ 2 ) , . . . a ( θ D ) ] T = [ 1 1 ⋯ 1 e − j φ 1 e − j φ 2 ⋯ e − j φ D ⋮ ⋮ ⋱ ⋮ e − j ( M − 1 ) φ 1 e − j ( M − 1 ) φ 2 ⋯ e − j ( M − 1 ) φ D ] A=[a(\theta_1),a(\theta_2),...a(\theta_D)]^T=\left[

1e−jφ1⋮e−j(M−1)φ11e−jφ2⋮e−j(M−1)φ2⋯⋯⋱⋯1e−jφD⋮e−j(M−1)φD11⋯1e−jφ1e−jφ2⋯e−jφD⋮⋮⋱⋮e−j(M−1)φ1e−j(M−1)φ2⋯e−j(M−1)φD

\right]A=[a(θ1),a(θ2),...a(θD)]T=⎣⎢⎢⎢⎡1e−jφ1⋮e−j(M−1)φ11e−jφ2⋮e−j(M−1)φ2⋯⋯⋱⋯1e−jφD⋮e−j(M−1)φD⎦⎥⎥⎥⎤

对x m ( t ) x_m(t)xm(t)进行N点采样，要处理的问题就变成了通过输出信号x m ( t ) x_m(t)xm(t)的采样 x m ( i ) = 1 , 2 , . . . , M \ x_m (i)=1,2,...,M\xm(i)=1,2,...,M估计信号源的波达方向角θ 1 , θ 2 . . . θ D \theta_1,\theta_2...\theta_Dθ1,θ2...θD,由此可以很自然的将阵列信号看作是噪声干扰的若干空间谐波的叠加，从而将波达方向估计问题与谱估计联系起来。

对阵列输出X做相关处理，得到其协方差矩阵

R x = E [ X X H ] R_x=E[XX^H]Rx=E[XXH]

其中H HH表示矩阵的共轭转置。

根据已假设信号与噪声互不相关、噪声为零均值白噪声，因此可得到：

R x = E [ ( A S + N ) ( A S + N ) H ] = A E [ S S H ] A H + E [ N N H ] = A R S A H + R N R_x=E[(AS+N)(AS+N)^H] =AE[SS^H]A^H+E[NN^H]=AR_SA^H+R_NRx=E[(AS+N)(AS+N)H]=AE[SSH]AH+E[NNH]=ARSAH+RN

其中R s = E [ S S H ] R_s=E[SS^H]Rs=E[SSH]称为信号相关矩阵

R N = σ 2 I R_N=\sigma^2IRN=σ2I是噪声相关阵

σ 2 \sigma^2σ2是噪声功率

I II是M × M M\times MM×M阶的单位矩阵

在实际应用中通常无法直接得到R x R_xRx，能使用的只有样本的协方差矩阵：

R x ^ = 1 N ∑ i = 1 N X ( i ) X H ( i ) \hatR_x=\frac1N \sum_i=1^NX(i)X^H (i)Rx^=N1∑i=1NX(i)XH(i)，R x ^ \hatR_xRx^是R x R_xRx的最大似然估计。

当采样数N → ∞ N\to\inftyN→∞，他们是一致的，但实际情况将由于样本数有限而造成误差。根据矩阵特征分解的理论，可对阵列协方差矩阵进行特征分解，首先考虑理想情况，即无噪声的情况：R x = A R s A H R_x=AR_sA^HRx=ARsAH，对均匀线阵，矩阵A由

A = [ a ( θ 1 ) , a ( θ 2 ) , . . . a ( θ D ) ] T = [ 1 1 ⋯ 1 e − j φ 1 e − j φ 2 ⋯ e − j φ D ⋮ ⋮ ⋱ ⋮ e − j ( M − 1 ) φ 1 e − j ( M − 1 ) φ 2 ⋯ e − j ( M − 1 ) φ D ] A=[a(\theta_1),a(\theta_2),...a(\theta_D)]^T=\left[

1e−jφ1⋮e−j(M−1)φ11e−jφ2⋮e−j(M−1)φ2⋯⋯⋱⋯1e−jφD⋮e−j(M−1)φD11⋯1e−jφ1e−jφ2⋯e−jφD⋮⋮⋱⋮e−j(M−1)φ1e−j(M−1)φ2⋯e−j(M−1)φD

\right]A=[a(θ1),a(θ2),...a(θD)]T=⎣⎢⎢⎢⎡1e−jφ1⋮e−j(M−1)φ11e−jφ2⋮e−j(M−1)φ2⋯⋯⋱⋯1e−jφD⋮e−j(M−1)φD⎦⎥⎥⎥⎤

所定义的范德蒙德矩阵，只要满足θ i ≠ θ j , i ≠ j \theta_i\neq \theta_j,i\neq jθi=θj,i=j，则他的各列相互独立。

若R s R_sRs为非奇异矩阵R a n k ( R s ) = D Rank(R_s)=DRank(Rs)=D，各信号源两两不相干，且M > D M>DM>D，则r a n d ( A R s A H ) = D rand(AR_sA^H)=Drand(ARsAH)=D，

由于R x = E [ X X H ] R_x=E[XX^H]Rx=E[XXH]，有：

R s H = R x R_s^H=R_xRsH=Rx

即R s R_sRs为Hermite矩阵，它的特性是都是实数，又由于R s R_sRs为正定的，因此A R s A … … H AR_sA……HARsA……H为半正定的，它有D个正特征值和M − D M-DM−D个零特征值。

再考虑有噪声存在的情况

R x = A R s A H + σ 2 I R_x=AR_sA^H+\sigma^2IRx=ARsAH+σ2I

由于σ 2 > 0 \sigma^2>0σ2>0，R x R_xRx为满秩阵，所以R x R_xRx有M个正实特征值λ 1 , λ 2 . . . λ M \lambda_1,\lambda_2...\lambda_Mλ1,λ2...λM

分别对应于M个特征向量v 1 , v 2 . . . v M v_1,v_2...v_Mv1,v2...vM。又由于R x R_xRx为Hermite矩阵，所以各特征向量是正交的，即：v i H v j = 0 , i ≠ j v_i^Hv_j=0,i\neq jviHvj=0,i=j与信号有关的特征值只有D个,分别等于矩阵A R s A H AR_sA^HARsAH的各特征值与σ 2 \sigma^2σ2之和，其余M − D M-DM−D个特征值为σ 2 \sigma^2σ2，即σ 2 \sigma^2σ2为R RR的最小特征值，它是M − D M-DM−D维的，对应的特征向量v i , i = 1 , 2 , . . . , M v_i,i=1,2,...,Mvi,i=1,2,...,M中，也有D个是与信号有关的，另外M − D M-DM−D个是与噪声有关的，可利用特征分解的性质求出信号源的波达方向θ k \theta_kθk。

MUSIC算法的原理及实现

通过对协方差矩阵的特征值分解，可得到如下结论：

将矩阵R x R_xRx的特征值进行从小到大的排序，即λ 1 ≥ λ 2 ≥ . . . ≥ λ M > 0 \lambda_1 \geq \lambda_2\geq...\geq\lambda_M>0λ1≥λ2≥...≥λM>0，其中D个较大的特征值对应于信号，M − D M-DM−D个较小的特征值对应于噪声。

矩阵R x R_xRx的属于这些特征值的特征向量也分别对应于各个信号和噪声，因此可把R x R_xRx的特征值（特征向量）划分为信号特征（特征向量）与噪声特征（特征向量）。

设λ i \lambda_iλi为R x R_xRx的第i ii个特征值，v i v_ivi是与λ i \lambda_iλi个相对应的特征向量，有：

R x v i = λ i v i R_xv_i=\lambda_iv_iRxvi=λivi

再设λ i = σ 2 \lambda_i=\sigma^2λi=σ2是R x R_xRx的最小特征值R x v i = σ 2 v i i = D + 1 , D + 2... M R_xv_i=\sigma^2v_i i=D+1,D+2...MRxvi=σ2vii=D+1,D+2...M，

将R x = A R s A H + σ 2 I R_x=AR_sA^H+\sigma^2IRx=ARsAH+σ2I代入可得σ 2 v i = ( A R s A H + σ 2 I ) v i \sigma^2v_i=(AR_sA^H+\sigma^2I)v_iσ2vi=(ARsAH+σ2I)vi，

将其右边展开与左边比较得：

A R s A H v i = 0 AR_sA^Hv_i=0ARsAHvi=0

因A H A A^HAAHA是D ∗ D D*DD∗D维的满秩矩阵，( A H A ) − 1 (A^HA)^-1(AHA)−1存在；

而R s − 1 R_s^-1Rs−1同样存在，则上式两边同乘以R s − 1 ( A H A ) − 1 A H R_s^-1(A^HA)^-1A^HRs−1(AHA)−1AH，

有：

R s − 1 ( A H A ) − 1 A H A R s A H v i = 0 R_s^-1(A^HA)^-1A^HAR_sA^Hv_i=0Rs−1(AHA)−1AHARsAHvi=0

于是有

A H v i = 0 ， i = D + 1 , D + 2 , . . . , M A^Hv_i=0，i=D+1,D+2,...,MAHvi=0，i=D+1,D+2,...,M

上式表明：噪声特征值所对应的特征向量（称为噪声特征向量）v i v_ivi，与矩阵A AA的列向量正交，而A AA的各列是与信号源的方向相对应的，这就是利用噪声特征向量求解信号源方向的出发点。

用各噪声特征向量为例，构造一个噪声矩阵E n E_nEn:

E n = [ v D + 1 , v D + 2 , . . . v M ] E_n=[v_D+1,v_D+2,...v_M]En=[vD+1,vD+2,...vM]

定义空间谱P m u ( θ ) P_mu(\theta)Pmu(θ):

P m u ( θ ) = 1 a H ( θ ) E n E n H a ( θ ) = 1 ∥ E n H a ( θ ) ∥ 2 P_mu(\theta)=\frac1a^H(\theta)E_nE_n^Ha(\theta)=\frac1\Vert E_n^Ha(\theta)\Vert^2Pmu(θ)=aH(θ)1EnEnHa(θ)=∥EnHa(θ)∥21

该式中分母是信号向量和噪声矩阵的内积，当a ( θ ) a(\theta)a(θ)和E n E_nEn的各列正交时，该分母为零，但由于噪声的存在，它实际上为一最小值，因此P m u ( θ ) P_mu(\theta)Pmu(θ)有一尖峰值，由该式，使θ \thetaθ变化，通过寻找波峰来估计到达角。

MUSIC算法实现的步骤

1.根据N个接收信号矢量得到下面协方差矩阵的估计值：

R x = 1 N ∑ i = 1 N X ( i ) X H ( i ) R_x=\frac1N \sum_i=1^NX(i)X^H(i)Rx=N1∑i=1NX(i)XH(i)

对上面得到的协方差矩阵进行特征分解

R x = A R s A H + σ 2 I R_x=AR_sA^H+\sigma^2IRx=ARsAH+σ2I

2.按特征值的大小排序将与信号个数D DD相等的特征值和对应的特征向量看做信号部分空间，将剩下的M − D M-DM−D个特征值和特征向量看做噪声部分空间，得到噪声矩阵E n E_nEn:

A H v i = 0 ， i = D + 1 , D + 2 , . . . , M A^Hv_i=0，i=D+1,D+2,...,MAHvi=0，i=D+1,D+2,...,M

E n = [ v D + 1 , v D + 2 , . . . v M ] E_n=[v_D+1,v_D+2,...v_M]En=[vD+1,vD+2,...vM]

3.使θ \thetaθ变化，按照式

P m u ( θ ) = 1 a H ( θ ) E n E n H a ( θ ) P_mu(\theta)=\frac1a^H(\theta)E_nE_n^Ha(\theta)Pmu(θ)=aH(θ)EnEnHa(θ)1

来计算谱函数，通过寻求峰值来得到波达方向的估计值。

clear; close all;

%%%%%%%% MUSIC for Uniform Linear Array%%%%%%%%

derad = pi/180; %角度->弧度

N = 8; % 阵元个数

M = 3; % 信源数目

theta = [-30 0 60]; % 待估计角度

snr = 10; % 信噪比

K = 512; % 快拍数

dd = 0.5; % 阵元间距

d=0:dd:(N-1)*dd;

A=exp(-1i*2*pi*d.'*sin(theta*derad)); %方向矢量

%%%%构建信号模型%%%%%

S=randn(M,K); %信源信号，入射信号

X=A*S; %构造接收信号

X1=awgn(X,snr,'measured'); %将白色高斯噪声添加到信号中

% 计算协方差矩阵

Rxx=X1*X1'/K;

% 特征值分解

[EV,D]=eig(Rxx); %特征值分解

EVA=diag(D)'; %将特征值矩阵对角线提取并转为一行

[EVA,I]=sort(EVA); %将特征值排序从小到大

EV=fliplr(EV(:,I)); % 对应特征矢量排序

% 遍历每个角度，计算空间谱

for iang = 1:361

angle(iang)=(iang-181)/2;

phim=derad*angle(iang);

a=exp(-1i*2*pi*d*sin(phim)).';

En=EV(:,M+1:N); % 取矩阵的第M+1到N列组成噪声子空间

Pmusic(iang)=1/(a'*En*En'*a);

end

Pmusic=abs(Pmusic);

Pmmax=max(Pmusic)

Pmusic=10*log10(Pmusic/Pmmax); % 归一化处理

h=plot(angle,Pmusic);

set(h,'Linewidth',2);

xlabel('入射角/(degree)');

ylabel('空间谱/(dB)');

set(gca, 'XTick',[-90:30:90]);

grid on;

实现结果

基于DOA联合TDOA时间积累的二维GDOP仿真分析

up目录

一、理论基础

二、核心程序

三、测试结果

一、理论基础

无人机(UAV)因其体积小,灵活性高,成本低等优势得到快速发展并被广泛应用于军事战争,城市管理,民用,地质,抢险救灾等各个领域,与此同时,无人机定位技术也得到了深入研究,其中无线电探测与定位技术备受众多学者关注.基于电磁波参数的无线电探测与定位技术主要包括角度估计,信号到达时间差估计(TDOA),接收信号强度估计(RSS),三种方法均可实现独立定位.信号到达时间差估计要求时钟同步,接收信号强度估计对信道环境比较敏感,而针对角度估计的联合定位算法普遍需要发射阵列和接收阵列角度配对.本文旨在提高无人机定位的估计精确度和降低联合估计过程中的运算复杂度,围绕多输入多输出(MIMO)双基地雷达下波离方向(DOD)和波达方向(DOA)联合估计技术开展研究,提出消除配对过程的联合估计算法,并结合角度和信号到达时间差实现无人机的位置定位.

首先介绍下TDOA的概念，如图所示，假设我们在空间中有一个声源(记为s(t)，其在空间的位置为S)两个麦克风(记为m1和m2，它们在空间的位置分别为M1和M2，接收到的信号为x1(t)和x2(t)

那么麦克风m1和m2收到的信号分别为：

其中τ1和τ2分别是声源到达两个麦克风的延迟时间，n1(t)和n2(t)为加性噪声。那么声源信号到达两个麦克风的TDOA为

其中c是声速。一般情况下，我们选择一个麦克风的信号作为参考信号，例如我们把M2作为参考信号，那么τ2=0。在麦克风阵列几何形状已知的情况下，声源定位问题变为对时延的估计问题。

DOA定位技术原理是利用接收机处的阵列天线和波达方向（DOA）估计技术来确定一个从接收机到信源的波达方向线，即为方向线（LOB），最后利用多个接收机估计的DOA进行三角测量，方向线的交点就是信源的估计位置。

波达方向（Direction Of Arrival）估计，又称为角谱估计（Angle spectral estimation）、波达角（Angle Of Arrival）估计。一个信源有很多可能的传播路径和到达角。如果几个发射机同时工作，每个信源在接收机处形成潜在的多径分量。因此，接收天线能估计出这些到达角就显得很重要，目的是估计出哪个发射机在工作以及发射机所处的方向，简单的说就是利用己方雷达接收来自目标发射机的来波方向进行估计；其物理原理是利用电磁波的直线传播原理。
通过测量辐射信号的波达方向（DOA全称Direction Of Arrival）或波达角（AOA）来估测辐射源的位置。理论上这种估计只需要两个接收阵元就可以确定辐射源的位置，但在实际中，由于受到角度分辨率、多径和噪声限制，所需阵元通常要多于两个。

二、核心程序

%-----------------------------------------------------------------------
%  --- 外辐射源基于DOA联合TDOA时间积累下二维平面GDOP分析 ---一发一收体制----
clc;
clear;
close all;
warning off;
addpath(genpath(pwd));

c=3e8;                         % 传播速度

 

% % 市区
xo=-0.25e3;yo=0;                 % 接收站1的位置   %市区 
xa=0.25e3;ya=0;                  % 发射站1的位置   %市区 
% 目标位置
xt=-2000:53:3000;
yt=-2000:53:2000;

da=pi/180;                     % 方位测量误差标准差
dtao=1e-7;                     % 时差测量误差标准差

R=[(da)^2 0;0 (dtao)^2];       % 测量误差协方差矩阵



for l=1:length(xt)
    for j=1:length(yt)
        
        a=(sqrt((xt(l)-xo)*(xt(l)-xo)+(yt(j)-yo)*(yt(j)-yo))+sqrt((xt(l)-xa)*(xt(l)-xa)+(yt(j)-ya)*(yt(j)-ya)))/2;
        tao=2*(a-xo)/c;
        
        afa_diff_x =-(yt(j)-yo)/(xt(l)-xo)^2/(1+(yt(j)-yo)^2/(xt(l)-xo)^2);
        afa_diff_y =1/(xt(l)-xo)/(1+(yt(j)-yo)^2/(xt(l)-xo)^2);
        
        f_diff_x =2*xt(l)/(xo+c*tao)^2;
        f_diff_y =2*yt(j)/(1/4*c^2*tao^2+c*tao*xo);
        f_diff_tao =-2*xt(l)^2/(xo+c*tao)^3*c-yt(j)^2/(1/4*c^2*tao^2+c*tao*xo)^2*(1/2*c^2*tao+c*xo);
        
        H=[afa_diff_x afa_diff_y;-f_diff_x/f_diff_tao -f_diff_y/f_diff_tao];
        
        Px=pinv(H)*R*pinv(H)';
        
        GDOP(l,j)=sqrt(Px(1,1)+Px(2,2));
        
    end
end

 
V=[0.005 0.007 0.01 0.015 0.02 0.025  0.03 0.035 0.04 0.045 0.05 0.055 0.06 0.08 0.10 0.15 0.2 0.3 0.5 1];%市区 %1e-7
figure(1);
[pic1]=contour(xt./1000,yt./1000,GDOP.'/1000,V);
clabel(pic1);xlabel('X(Km) ');ylabel('Y(Km) ');title( 'CONTOUR of valleys of GDOP in T_1R_1 (Km)' );
hold on;
plot(xo/1000,yo/1000,'gP',xa./1000,ya./1000,'rP');
axis([-1.25,1.75,-1.5,1.5]);%市区
% axis([-30,50,-40,40]);%郊区
grid on