当前位置：首页 > article >正文

DOA估计算法——ESPRIT算法

article 2025/1/14 0:53:51

1 简介

ESPRIT（Estimation of Signal Parameters via Rotational Invariance Techniques）最早是由Roy等人于1986年提出，是一种广泛应用于高分辨率方向到达（DOA）估计和频率估计的子空间方法。其核心思想基于信号子空间的旋转不变性，能够通过有限的快拍数据实现对信号参数的高精度估计。与经典的MUSIC算法相比，ESPRIT不依赖于复杂的空间谱搜索过程，而是通过结构化的阵列设计（如双阵列或具有重复结构的阵列）直接解算信号的参数，具有较低的计算复杂度。由于在参数估计等方面的优越性，ESPRIT算法近年来得到了广泛应用，并出现了许多变化算法，如LS-ESPRIT算法、TLS-ESPRIT算法、SVD-ESPRIT算法、多重不变ESPRIT算法；波束空间ESPRIT算法和酉ESPRIT算法等。

ESPRIT算法的主要优势在于其对阵列结构要求相对较低，同时能够在较小样本数和较高噪声条件下保持优良的性能。这使得其在雷达、声纳、通信、天文观测等领域得到了广泛应用。此外，由于其求解过程中采用了特征值分解和最小二乘方法，ESPRIT具有良好的稳健性和解析能力，为多信号处理提供了可靠的工具。随着对信号处理需求的不断提高，ESPRIT算法的研究和扩展也在持续进行，例如推广到二维或三维DOA估计、结合稀疏优化技术以提升分辨能力，以及适应动态场景的实时处理等。这些进展进一步拓宽了ESPRIT算法的应用范围，增强了其实际价值。

2 ESPRIT算法原理

如图1所示的具有N个阵元的均匀线阵。假设通过该线阵在一定观测时长接收到的阵列数据为X，则其无偏估计的协方差矩阵可以表示为：

其中L表示快拍数。对协方差矩阵进行特征分解将其分解为噪声子空间Un可信号子空间Us ：

其中分别表示信号子空间和噪声子空间对应的特征值对角矩阵。根据阵列信号模型与矩阵分解理论，具有两个重要的引理：

信号子空间与导向矢量张成的子空间属于同一个子空间，即span{Us}=span{A }.
信号子空间与噪声子空间正交，即.

值得一提的是，MUSIC算法正是基于信号子空间与噪声子空间正交的特性推导的，关于MUSIC算法的讲解我将另外一篇文章中讲解。现在我们知道Us和A属于同一个子空间，基于该引理，Roy首先假设存在一个唯一的满秩矩阵T，使得下式成立：

紧接着他又将一个含有N个阵元的线性阵列拆解为两个子阵（如图1所示），其中子阵1有前N-1个阵元组成，子阵2由后N-1个阵元组成。于是分别得到子阵1和子阵2的信号子空间Ex和Ey，它们满足下式关系：

Ey=ExΦ (4)

其中是信号相位变化矩阵，与待估计参数相关。式(4)正是利用了旋转不变的特性，因此结合(3)、(4)可以得到：

(5)

根据式(5)进一步可可得

(6)

其中。至此可知，Ex和Ey张成相似的子空间，且矩阵Φ的对角线元素为Ψ的特征值。

总结：关于ESPRIT算法，它的实际过程归结于求解式(6)中的Ψ,显然这个问题可以等价于求解方程AX=B的问题，求解该问题可以通过LS、TLS、SLS等方法求解，然后根据ESPRIT思想，旋转因子Φ与Ψ相似，并根据矩阵相似的理论，故它们有相同的特征值，因此通过特征值可以估计DOA。

3 算法仿真

根据第2部分中原理，我们知道求解ESPRIT算法归结为欠定方程求解问题，因此在本demo案例中使用TLS（总体最小二乘）进行求解，当然读者如何感兴趣，可以使用LS、SLS等其它一些方法进行求解，值得注意的是，每种方法各有优缺点：以下是TLS-ESPRIT算法求解步骤：

Step1 :根据阵列接收的矩阵X求协方差矩阵R.

Step2 :对协方差矩阵进行特征分解，取K个最大特征值对应的特征向量构成信号子空间，并分为Ex

和Ey两部分。

Step3 : 计算下式的特征值分解：

并分解成K×K的子矩阵

Step4 :计算的特征值

Step5 : 计算到达角估计值

3.1 仿真代码示例

仿真环境：Matalb2021b，Windows11

阵元数量：M=16

目标来波方向：40°，10°，20°的非相干信号源

信号的中心频率为：f = 77GHz，信号能量幅度默认为1。

信噪比：SNR=15dB

快拍数：N = 1024

%%  Author：Poulen
%%  Data：2024/12/15
%%  TLS-ESPRIT算法仿真
%%  总体最小二乘ESPRIT算法，考虑到信号子空间存在估计误差的问题，利用总体最小二乘法来求解ESPRIT方程
clear 
close all;
clc;

%产生信号
%仿真初始值设定
M=16;         %阵元单元
c=3e8;       %光速
f0=77e9;     %初始频率
lambda=c/f0; %波长
slope=30e12; %调频斜率
time=60e-6;  %60us
d=0:lambda/2:(M-1)*lambda/2;%阵列天线
thita=[-40];%波达方向
K = length(thita);

%产生原始信号
N=1024;%信号长度
t=linspace(0,time,N);

A=zeros(M,K);%导向向量空间 M*K
S=zeros(K,N);%信号空间

% St = exp(1j*2*pi*(f0*t(:)+1/2*slope*t(:).^2));
f = 100+f0; %非相干信号的频率
for i = 1:K
    A(:,i) = exp(-1j*2*pi/lambda*d(:)*sind(thita(i)));
%     S(i,:) = St;
    S(i,:) = exp(1j*2*pi*f*t(:));
    f = 10000+f;
end
S = A*S; %产生阵列接收数据

%向数据添加白噪声
SNR = 30; %单位dB
S = S +(randn(size(S)).*std(S))/db2mag(SNR);

%TLS-ESPRIT算法
%step1：计算阵列输出的协方差矩阵
%step2: 对改进协方差矩阵进行特征分解，求解信号子空间Us
%step3: 划分信号子空间为两个子阵对应的子空间，Esx、Esy
%step4: 合并子阵1和子阵2的信号子空间，并求P=Usxy'*Usxy
%step5: 对P进行特征分解，并将2K*2K的特征矩阵E进行分块，分成K*K矩阵，即E = [E11 E12;E21 E22];
%step5: 根据F=-E12 * inv(E22)，并对F特征分解求取特征值
%step6: 根据特征值，估计DOA

%协方差估计
R = (1/N)*S*conj(S.');

%重构协方差矩阵
[V,D] = eig(R);   %对于实数来说已经从小到大排列，复数则随机排列

%重新排列特征值大小（按照实部）
EVA = real(diag(D)');  
[EVA,I] = sort(EVA);
V=V(:,I);


Us = V(:,(end-K+1):end);%信号子空间

Usx = Us(1:end-1,:);
Usy = Us(2:end,:);
% Uxy = [Usx Usy];

%对Uxy' * Uxy 进行特征分解
% Exy1 = Uxy' * Uxy;

Exy =  [Usx Usy]'*[Usx Usy];

[E,~] = eig(Exy);

E11 = E(1:K,1:K);
E12 = E(1:K,K+1:2*K);
E21 = E(K+1:2*K,1:K);
E22 = E(K+1:2*K,K+1:2*K);

F = -E12 * (inv(E22)) ;
F_eig = eig(F);

%DOA估计
DOAEs = asind(-angle(F_eig)/pi);
disp(sort(DOAEs));


figure;
theta=deg2rad(DOAEs); %转化为弧度值
polarplot(theta,8,'Marker','^','MarkerSize',6,'Color',[1 0 0],'LineWidth',1.2);
%设置极坐标属性
ax = gca;
ax.ThetaLim = [-60 60];
ax.ThetaDir = 'clockwise';
ax.ThetaZeroLocation = 'top';
ax.RLim = [0 10];
ax.RGrid = 'on' ;
Targetnum = length(DOAEs);
title(['共搜寻出 ' num2str(Targetnum) ' 个目标 ']);