DOA估计算法从原理到仿真——CBF、Capon、MUSIC算法

从原理出发理解并实现经典的DOA估计算法

m0_55845774

1143人浏览 · 2024-11-17 14:18:33

m0_55845774 · 2024-11-17 14:18:33 发布

本人第一篇CSDN博客~欢迎关注！

DOA是指Direction Of Arrival，是利用电磁波信号来获取目标或信源相对天线阵列的角度信息的方式，也称测向、空间谱估计。主要应用于雷达、通信、电子对抗和侦察等领域。

一、阵列信号模型

如上图所示，假设有 $M$ 个阵元以 $d$ 等间距线性排列组成接收阵列（这里 $d$ 的选取我们在后文进一步讨论），设目标距离第一个阵元距离为 $R$ ，方向为 $\theta$ ，其中第 $m$ 个阵元对该目标的回波信号为 :

远场条件下，信号传输模型近似与平面波，而近场则可视作球面波。

下面的DOA估计都将按照远场条件来近似。我们先要重构对于阵列接收和发射信号的表达方式，采用矩阵的形式对各个阵元收发的信号进行建模：

首先在连续时域上，假设存在M个阵元，有P个目标回波（分别来自P个不同方向），那么各个阵元接收到的信号可以表示为:

其中，矢量X代表阵元接收回波的信号，A代表空域导向矢量，S代表发射信号矢量，N代表噪声信号，经过采样转换到离散域下有：

即：

另外值得一提的是，这里的角度变量 $\theta$ 也是一个连续值，在仿真分析中我们也要将其离散化（这里在之后的仿真中我们具体考虑如何实现）。

针对这样的接收信号，我们再构建下面这样的阵列接收系统模型：

其中， $W_n$ 称作阵元的权矢量，这个权矢量可以理解为一个空域的“滤波器”，它决定了各个阵元所接收到的信号经过怎样的增益值后再进行加权得到整个阵列的输出，经过权矢量的“滤波”效果，我们希望阵列接收到的信号在目标来向角度下有更大的增益，而在别的角度则产生衰减效果，这种“滤波”效果也实际上附带了目标的来向角度信息，从而实现DOA估计的效果，其阵元输出表达式：

二、算法原理

（1）常规波束形成算法CBF

CBF（Conventional BeamForming）是最简单的一种DOA估计算法。其做法是将权矢量矩阵W定义为空域导向矢量矩阵，即考虑最简单的情况。从计算上来看，它表征了不同方向角下信号的功率分布情况。类比于功率谱密度函数和自相关函数的关系（维纳-辛钦定理），这样的算法有点像对信号作傅里叶变换，教材和课件上把这个称作空域的傅里叶变换。其表达式：

这种算法优点是实现起来非常简单，但是它受限于瑞利极限，其角分辨率并不高，并且从计算上我们可以看出，这种类似傅里叶变换的计算方式得到的会是一个类似sinc函数的表达式，这会在谱线上引入一系列的栅瓣。它相当于是一种简单的，被动的DOA估计算法。

（2）Capon算法

在解释Capon算法的原理前，我们先了解一些数学知识：

首先回想我们求解多元函数极值问题时的步骤：通常来说，我们求解一个函数的极值时，需要先考虑它是否存在约束条件，如果不存在约束条件，使用求导的方法就可以解决，但是如果存在约束条件，常规多元函数求极值问题可以抽象为：

保证约束条件的同时，找到目标函数的极值。

在针对波束成形/DOA估计时，我们希望解决的问题是：让阵列能够正常接收到目标来向上的信号，并且让非目标来向角度下的增益变得更小。根据这个，Capon算法把波束成形当做了一个抽象化的求极值问题：

这种方法也叫做最小方差无失真响应(Minimum Variance Distortion-less Response,MVDR)法。不同于CBF，它是一种具备自适应调节能力的超分辨率的DOA方法，不受瑞利极限限制，具有更好的角度分辨率，并且无需知道信源个数作为先验知识，这种思路其实和维纳滤波器的原理一模一样。

（3）MUSIC算法

同样需要先补充一些数学知识（这边参考SYSU阵列信号处理课程课件以及科学出版社高等工程数学第一章的部分内容）：

由于在前文中，我们已经推导出阵列接收信号 $y(t)$ 的自相关矩阵 $R$ 实际上是一个Hermite矩阵，根据Schur定理，则存在酉矩阵 $U$ 可将其进行酉分解，即

$\Sigma = diag(\lambda _1,....,\lambda _n)=U^H R U$

变换得到

$R = U \Sigma U^H = \Sigma _{i=1} ^{n}\lambda _ie_ie^H_i$

其中 $U=[e_1,e_2,...e_M]$ 为特征矢量构成的酉矩阵，且特征值满足从大到小排序。参考上述对线性空间和子空间的定义，我们知道 $U$ 实际上可以考虑作一个由特征矢量张成的线性空间，而对于自相关矩阵 $R$ ，定义上有：

可以发现这个自相关矩阵其实是由信号 $s(n)$ 与噪声 $n(n)$ 两种成分组成，并且对于噪声（白噪声）这种随机过程，其与任意其他函数是不相关的（相互独立），因此在直觉上看 $R$ 是可以通过某种方法分解成信号与噪声两部分的。那么该如何分解呢？

这时候我们再考虑特征值分解（EVD）的几何意义，实际上EVD在几何上实现的就是将一个空间中的所有 “经过线性变换后不改变'形状'，只改变长度” 的所有向量找出来，并且特征值对应的就是其“改变长度”的比例，或者叫做“拉伸/缩放比例”，而在酉矩阵的EVD中，重要特性是它保持向量的范数和内积不变，因此可以看作是纯旋转或反射矩阵，换句话说，对于一个酉矩阵，其特征值总是复数模长为1的复数，因此它只会改变特征向量的方向，而不会改变其大小。

在上图对 $R$ 的分解中可以看出， $R_s$ 是信号子空间的矩阵，其EVD后对角元素是较大的特征值，反映信号的能量。剩下的部分是噪声子空间的特征值对角矩阵，其EVD对角元素是较小的特征值，通常接近噪声功率 $\sigma^2$ 。

由于信号子空间 $\mathbf{E}_s$ 对应的是由信号源贡献的协方差矩阵 $\mathbf{A} \mathbf{R}_s \mathbf{A}^H$ ，其能量通常较大，因此信号空间的特征值 $\mathbf{\Lambda}_s$ 较大。这是因为信号子空间承载了信号的能量，其特征值远大于噪声的功率。而噪声子空间 $\mathbf{E}_n$ 对应的是噪声协方差矩阵 $\sigma^2 \mathbf{I}$ ，其特征值接近噪声功率 $\sigma^2$ 。由于噪声是均匀分布的，其特征值较小，通常接近 $\sigma^2$ 并且彼此相等。

至此我们就完成了对接收信号的信号空间与噪声空间的分解，即：

$R = U_s \Sigma_s U_s^H +U_N \Sigma_N U_N^H = AR_sA^H + \sigma ^2I$

对上式同时右乘 $U_N$

$RU_N =U_s \Sigma_s U_s^HU_N +U_N \Sigma_N U_N^H U_N = AR_sA^H U_N + \sigma ^2IU_N$

由于信号子空间与噪声子空间相互正交，其内积为0，则上式

$=0+U_N \Sigma_N I_N = AR_sA^H U_N + \sigma ^2U_N$

$\Rightarrow \sigma^2 U_N= AR_sA^H U_N + \sigma ^2U_N$

$\Rightarrow AR_sA^H U_N =0$

左乘 $U_N^H$ 得到：

$U_N^HAR_sA^H U_N =0$

令 $X = A^HU_N$ :

$X^H R_s X = 0$

由于 $X$ 非奇异时，当且仅当 $X=0$ 时上式才成立，故得到：

$A^H U_N = 0$

由于 $A$ 为目标的导向矢量矩阵，即目标的导向矢量所张成的子空间与噪声子空间正交。此时将目标导向矢量矩阵中的目标方向替换为扫描的形式，得到：

$A^H(\theta) U_N = 0$

其中 $\theta$ 代表对角度纲量进行“扫描”的自变量

为了得到最后标量形式的DOA估计，定义一种类似于功率谱的函数

$P(\theta) = \frac{1}{A^H(\theta)U_NU_N ^H A(\theta)}$

则为最后的MUSIC方法得到的DOA估计表达式。

三、MATLAB仿真实现

仿真实现了三种算法的DOA估计，并且探究了不同快拍数、不同阵元数量与不同信噪比下的三种方法DOA估计效果。

clc;close all;clear;

%DOA algorithms : CBF_Capon_MUSIC
%Author:ChrisZhou 2024.10 @ SYSU

function [doa_cbf,doa_capon,doa_music,angleIndex] = myDOA(arrayNum,snapshotLength,SNR)
%Parameters:
c=3e8;%光速
fc=1.5e9;    %雷达载频
lambda=c/fc; %波长
d=lambda/2;  %阵元间距
theta_T=[10,20].'; %目标来向角
targetNum=length(theta_T);   %目标个数
%arrayNum=16;                 %阵元个数
arrayPos = 0:d:(arrayNum-1) * d; %阵元位置
sampleRate= 48e5;            %采样率
ts = 1/sampleRate;           %采样间隔
L = snapshotLength;          %快拍数

angleIndex = asin((-512:1:512-1)/512) * 180 /pi; %角度纲量
angleIndexLength = length(angleIndex);
timeIndex = (0:1:L-1)*ts;                        %时间纲量

%Signal Genration

As = zeros(arrayNum,targetNum);  % A(theta)
S_signal = zeros(targetNum,L);   % s(n)
%分别定义A(theta) 与 s(n)
for i=1:1:targetNum
    As(:,i) = exp(-1j*(2*pi/lambda).*arrayPos.'*(sind(theta_T(i))));
    S_signal(i,:) = exp(1j*2*pi*fc*timeIndex);
    fc=fc+1000;  %假设后一个目标回波信号的频率高1000Hz
end

S_signal = As * S_signal;      % X(n) = A(theta) * S(n)
%SNR = 20;                      % 设定信噪比
S_signal = awgn(S_signal,SNR); % 添加高斯白噪声 即x(n) = x(n) + N(n)

%定义权矢量矩阵 W
W = zeros(arrayNum,angleIndexLength);
for i=1:angleIndexLength
    W(:,i)=exp(-1j*(2*pi/lambda).*arrayPos.'*sind(angleIndex(i)));
end

%信号协方差矩阵的最大似然估计（标准化后的自相关矩阵）
R_hat = (S_signal * conj(S_signal).' ) / L;
R_hatInv = inv(R_hat);
%DOA Algorithms

%CBF方法
doa_cbf=zeros(angleIndexLength,1);

%P(w)=W^H R W
for i=1:1:angleIndexLength
    doa_cbf(i)=conj(W(:,i).') * R_hat * W(:,i);
end


%Capon方法
doa_capon = zeros(angleIndexLength,1);

%P(w) = 1 / A^H R^-1 
for i=1:1:angleIndexLength
    doa_capon(i)= 1 / (conj(W(:,i).')* R_hatInv * W(:,i));
end

%MUSIC方法
doa_music = zeros(angleIndexLength,1);
[EV,D] = eig(R_hat);    %将协方差矩阵进行特征分解
EVA = diag(D)';         %提取特征值矩阵对角线，并转为一行
[EVA,I] = sort(EVA);    %将特征值从小到大排序
EV = fliplr(EV(:,I));   %对应特征矢量的排序
Un = EV(:,targetNum+1:arrayNum);
%P(w) = 1/ A*Un*Un^H*A
for i = 1:1:length(angleIndex)
    doa_music(i)=1/(conj(W(:,i).') * Un * Un' * W(:,i));
end



% ESPRIT方法
% 首先，对信号协方差矩阵进行特征分解
[EV, D] = eig(R_hat);
EVA = diag(D)';
[EVA, I] = sort(EVA); % 特征值排序
EV = fliplr(EV(:, I)); % 特征矢量排序

% 选择信号子空间和噪声子空间
Un = EV(:, 1:targetNum); % 信号子空间
Un_plus = Un(2:end, :); % 上部分
Un_minus = Un(1:end-1, :); % 下部分

% 计算phi矩阵
Phi = pinv(Un_minus) * Un_plus;

% 计算特征值
[~, D_phi] = eig(Phi);
eigenvalues = diag(D_phi);
angles = angle(eigenvalues); % 提取相位

% 计算估计的DOA
doa_esprit = asin(lambda * angles / (2 * pi * d)) * (180/pi); % 转换为角度

% 显示结果
disp('DOA Estimates using ESPRIT:');
disp(doa_esprit);


end


%%plot

%Case 1 :SNR = 20dB时，估计 R 快拍数 = 500，阵元数为 8、16、24三种情况；
[c8_cbf,c8_capon,c8_music,angleIndex] = myDOA(8,500,20);
[c16_cbf,c16_capon,c16_music] = myDOA(16,500,20);
[c24_cbf,c24_capon,c24_music] = myDOA(24,500,20);

No=1;
figure(No);
%theta = linspace(-90, 90, 1024);

subplot(311);
plot(angleIndex,db(c8_cbf/max(c8_cbf)),'b','DisplayName', 'CBF','LineWidth',1);
hold on;
plot(angleIndex,db(c8_capon/max(c8_capon)),'r','DisplayName', 'Capon','LineWidth',1);
hold on;
plot(angleIndex,db(c8_music/max(c8_music)),'g','DisplayName', 'Music','LineWidth',1);
xlabel('Angle(°)');ylabel('dB');title('M=8 Snapshot=500 SNR = 20');
xlim([-90,90]);
legend;
grid on;

subplot(312);
plot(angleIndex,db(c16_cbf/max(c16_cbf)),'b','DisplayName', 'CBF','LineWidth',1);
hold on;
plot(angleIndex,db(c16_capon/max(c16_capon)),'r','DisplayName', 'Capon','LineWidth',1);
hold on;
plot(angleIndex,db(c16_music/max(c16_music)),'g','DisplayName', 'Music','LineWidth',1);
xlabel('Angle(°)');ylabel('dB');title('M=16 Snapshot=500 SNR = 20');
xlim([-90,90]);
legend;
grid on;

subplot(313);
plot(angleIndex,db(c24_cbf/max(c24_cbf)),'b','DisplayName', 'CBF','LineWidth',1);
hold on;
plot(angleIndex,db(c24_capon/max(c24_capon)),'r','DisplayName', 'Capon','LineWidth',1);
hold on;
plot(angleIndex,db(c24_music/max(c24_music)),'g','DisplayName', 'Music','LineWidth',1);
xlabel('Angle(°)');ylabel('dB');title('M=24 Snapshot=500 SNR = 20');
xlim([-90,90]);
legend;
grid on;

%Case 2 : 估计 R 快拍数 = 500，阵元数为 16，SNR 为 5dB、10dB 和 20dB 三种情况；
No=No+1;
figure(No);

[n5_cbf,n5_capon,n5_music] = myDOA(16,500,5);
[n10_cbf,n10_capon,n10_music] = myDOA(16,500,10);
[n20_cbf,n20_capon,n20_music] = myDOA(16,500,20);

subplot(311);
plot(angleIndex,db(n5_cbf/max(n5_cbf)),'b','DisplayName', 'CBF','LineWidth',1);
hold on;
plot(angleIndex,db(n5_capon/max(n5_capon)),'r','DisplayName', 'Capon','LineWidth',1);
hold on;
plot(angleIndex,db(n5_music/max(n5_music)),'g','DisplayName', 'Music','LineWidth',1);
xlabel('Angle(°)');ylabel('dB');title('M=16 Snapshot=500 SNR = 5');
xlim([-90,90]);
legend;
grid on;

subplot(312);
plot(angleIndex,db(n10_cbf/max(n10_cbf)),'b','DisplayName', 'CBF','LineWidth',1);
hold on;
plot(angleIndex,db(n10_capon/max(n10_capon)),'r','DisplayName', 'Capon','LineWidth',1);
hold on;
plot(angleIndex,db(n10_music/max(n10_music)),'g','DisplayName', 'Music','LineWidth',1);
xlabel('Angle(°)');ylabel('dB');title('M=16 Snapshot=500 SNR = 10');
xlim([-90,90]);
legend;
grid on;

subplot(313);
plot(angleIndex,db(n20_cbf/max(n20_cbf)),'b','DisplayName', 'CBF','LineWidth',1);
hold on;
plot(angleIndex,db(n20_capon/max(n20_capon)),'r','DisplayName', 'Capon','LineWidth',1);
hold on;
plot(angleIndex,db(n20_music/max(n20_music)),'g','DisplayName', 'Music','LineWidth',1);
xlabel('Angle(°)');ylabel('dB');title('M=24 Snapshot=500 SNR = 20');
xlim([-90,90]);
legend;
grid on;

%Case 3 :SNR = 20dB时，阵元数为 16，估计 R 快拍数为：10、20、50和 100 四种情况

No=No+1;
figure(No);

[s10_cbf,s10_capon,s10_music] = myDOA(16,10,20);
[s20_cbf,s20_capon,s20_music] = myDOA(16,20,20);
[s50_cbf,s50_capon,s50_music] = myDOA(16,50,20);
[s100_cbf,s100_capon,s100_music] = myDOA(16,100,20);

subplot(411);
plot(angleIndex,db(s10_cbf/max(s10_cbf)),'b','DisplayName', 'CBF','LineWidth',1);
hold on;
plot(angleIndex,db(s10_capon/max(s10_capon)),'r','DisplayName', 'Capon','LineWidth',1);
hold on;
plot(angleIndex,db(s10_music/max(s10_music)),'g','DisplayName', 'Music','LineWidth',1);
xlabel('Angle(°)');ylabel('dB');title('M=16 Snapshot=10 SNR = 20');
xlim([-90,90]);
legend;
grid on;

subplot(412);
plot(angleIndex,db(s20_cbf/max(s20_cbf)),'b','DisplayName', 'CBF','LineWidth',1);
hold on;
plot(angleIndex,db(s20_capon/max(s20_capon)),'r','DisplayName', 'Capon','LineWidth',1);
hold on;
plot(angleIndex,db(s20_music/max(s20_music)),'g','DisplayName', 'Music','LineWidth',1);
xlabel('Angle(°)');ylabel('dB');title('M=16 Snapshot=20 SNR = 20');
xlim([-90,90]);
legend;
grid on;

subplot(413);
plot(angleIndex,db(s50_cbf/max(s50_cbf)),'b','DisplayName', 'CBF','LineWidth',1);
hold on;
plot(angleIndex,db(s50_capon/max(s50_capon)),'r','DisplayName', 'Capon','LineWidth',1);
hold on;
plot(angleIndex,db(s50_music/max(s50_music)),'g','DisplayName', 'Music','LineWidth',1);
xlabel('Angle(°)');ylabel('dB');title('M=16 Snapshot=50 SNR = 20');
xlim([-90,90]);
legend;
grid on;

subplot(414);
plot(angleIndex,db(s100_cbf/max(s100_cbf)),'b','DisplayName', 'CBF','LineWidth',1);
hold on;
plot(angleIndex,db(s100_capon/max(s100_capon)),'r','DisplayName', 'Capon','LineWidth',1);
hold on;
plot(angleIndex,db(s100_music/max(s100_music)),'g','DisplayName', 'Music','LineWidth',1);
xlabel('Angle(°)');ylabel('dB');title('M=16 Snapshot=100 SNR = 20');
xlim([-90,90]);
legend;
grid on;

部分运行结果：