DOA估计算法——Capon算法

1.波速形成基本思想

在理解Capon算法之前，我们有必要先了解波束形成的基本思想以及原理到底是什么。这有助于我们更好的理解Capon算法的思想。

如图1展示了均匀阵列波束导向的示意图。图中wm表示加权值，波速形成(DBF)的基本思想就是将各阵元输出进行加权求和，在一定时间内将天线阵列“导向”到一个方向上，对期望信号得到最大输出功率的导向位置，同时这个位置也表征了目标或波达方向。

此外，DBF的基本原理简单来说，就是利用阵元直接相干叠加而获得输出，其缺点在于只有垂直于阵列平面方向的入射波在阵列输出端才能同相叠加，从而形成方向图中主瓣的极大值。反过来说，如果阵列可以围绕它的中心轴旋转，那么当阵列输出最大时，空间波必然由垂直于阵列平面的方向入射而来。但有些天线阵列是很庞大的，且是不能转动的。因此，设法设计一种相控天线法(或称常规波束形成法)，这是最早出现的阵列信号处理方法。这种方法中，阵列输出选取一个合适的加权向量以补偿各个阵元的传播延时，从而使在某一期望方向上阵列输出可以同相叠加，进而使阵列在该方向上产生一个主瓣波束，而对其他方向上产生较小的响应，用这种方法对整个空间进行波束扫描就可确定空中待测信号的方位。因此，计算权值wm是波束形成这类方法的关键，目前DBF权重主要基于以下准则进行计算：

最大信噪比准则(MSNR)：使期望信号分量功率与噪声分量功率之比最大，但是必须知道噪声的统计量和期望信号的波达方向。
最大信干噪比准则(MSINR)：使期望信号功率与干扰功率及噪声分量功率之和的比最大。
最小均方误差准则(MMSE)：在非雷达应用中，阵列协方差矩阵中通常都含有期望信号，基于此种情况提出的准则。使得阵列输出与某期望响应的均方误差最小，不需要知道期望信号的波达方向。
最大似然比准则(MLH)：在对有用信号完全先验未知的情况下，参考信号无法设置，因此，在干扰噪声背景下，首先要取得对有用信号的最大似然估计。
线性约束最小方差准则(LCMV)：对有用信号形式和来向完全已知，在某种约束条件下使阵列输出的方差最小。

2.Capon算法

Capon算法属于一种在线性约束最小方差准则(LCMV)下的波束形成算法，所谓波束形成即在某些准则约束下，求解阵列输出的最优权。实质上波束形成属于一种空域滤波器。波束形成的“导向”作用是通过调整加权系数完成的,阵列的输出是对各阵元的接收信号量x(n)在各阵元上的加权和，令权向量为w,则输出可写作为：

可见对不同的权向量，上式对来自不同方向的信号有不同的响应从而形成不同方向的空间波束。假设空间远场有一个感兴趣的信号d(t) (其波达方向为θd)和J个干扰信号,1,..J((其波达方向为θi),令每个阵元上的加性白噪声为nk(t),它们具有相同的方差。在这些假定的条件下，第k个阵元上的接收信号可以表示为

等式右边的三项分别表示信号，干扰和噪声。假设有M个阵元如果用矩阵的形式表示式(2),则有

其中，分别为M个阵元上接收的数据；

表示波达方向来自的方向向量。假设接收端得到了N个快拍数据，则根据式(1)可得到波形形成器的输出的平均功率为：

当式(4)可以表示为：

如上式子所示，输出功率关于波达角的函数通常被称为空间谱。为了保证来自方向的信号正常接收，同时完全抑制掉其它J个干扰，很容易根据式子(5)得到权向量的约束条件为：

综上所述，我们可以得出Capon算法所要求解的优化问题可以表述为：

可以采用拉格朗日乘子法求解。求解过程如下：

令：

L分别对w和λ求偏导并使各自的偏导数为零，由此可得：

将式(9)中的第二个式子左乘并将第一个式子代入可得

式(10)右乘，可得

于是得

式(12)右乘可得

于是得

将式(14)代入式(12)并取复共轭转置可得：

将式(15)代入式(5)即可得到Capon的空间谱为:

根据式(16)进行谱峰搜索，峰值所在的索引表征了目标波达方向(DOA)。

3.Capon算法仿真

仿真环境：Matlab2021b；

波达方向分别为：10°、-15°、21°；

信源快拍数、阵元个数：1024、32;

分别对信噪比为:-10dB、20dB、30dB

仿真代码如下：

%%  Author：Poulen
%%  Data：2023.5.29
%%  Capon算法仿真
clear 
close all;
clc;

%%  产生信号
M=32;                       %阵元单元
c=3e8;                      %光速
f0=77e9;                    %初始频率
lambda=c/f0;                %波长
slope=30e12;                %调频斜率
time=60e-6;                 %60us
d=0:lambda/2:(M-1)*lambda/2;%阵列天线
thita=[-15,21,10];          %波达方向
K = length(thita);
N=1024;                     %信号长度
t=linspace(0,time,N);
A=zeros(M,K);               %导向向量空间 M*K
S=zeros(K,N);               %信号空间

f = 100+f0; 
for i = 1:K
    A(:,i) = exp(-1j*2*pi/lambda*d(:)*sind(thita(i)));
    S(i,:) = exp(1j*2*pi*f*t(:));
    f = 1000+f;
end
S = A*S;                    %产生阵列接收数据

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

%%  计算信号协方差矩阵
R = (1/N).*S*conj(S).';
R_inv = inv(R);

%%  Capon 算法
scale = -60:0.1:59;         %扫面范围
P_Capon = zeros(length(scale),1);
idx = 1;
for i=scale
    a = exp(-1j*2*pi/lambda*d(:)*sind(i));
    P_Capon(idx) = 1/(conj(a.')*R_inv*a); 
    idx = idx + 1;
end
figure;
plot(scale,db(P_Capon),'LineWidth',1.6,'Color',[0 0 0]);
xlabel('AngleRange(°)');
ylabel('Amplitude(dB)');
title('CaponAlgorithm');
legend(' SNR = 60dB');

仿真效果如下图所示：