基于L型阵列的2维DOA估计

摘要：在实际的雷达目标DOA估计过程中，局限于1维DOA估计显然是不能满足现代雷达系统探测目标的要求，例如主流的4D雷达成像，不仅要估计目标的距离、速度、方位(水平)角外，还要估计目标的俯仰角，然而我们要同时获得目标的水平-俯仰向信息，在硬件层面，就需要拥有两个维度的天线。因此，本期内容通过L型阵列，并采用Music算法进行谱峰扫描实现目标的水平-俯仰向信息的估计，希望对各位读者带来帮助。

1 . L形阵列

        如图1所示为L形阵列的阵列结构，有M+N-1个阵元。此L型阵列由x轴上阵元数为N的均匀线阵和y轴上阵元数为M的均匀线阵构成，阵列间距为d。假设空间有K个信源照射到此阵上。其二维波达方向为（θk,ϕk），k = 1,2,3,…,K，其中,θk和ϕk分别表示第k个信源的俯仰角和方位角。

图 1 L形阵列的阵列结构

假设入射到此阵列上的信源数为K，则x轴上N个阵元对应的导向矩阵为 

y轴上M个阵元对应的导向矩阵为 

2. MUSIC算法 

        关于MUSIC算法是一种非常经典的超分辨率算法，有关于它的原理在此不在赘述。我们直接给出MUSIC算法表达式：

式中，a(ϕ,θ)表示二维搜索的导向矢量，Un表示噪声子空间。2D MUSIC算法就是在水平向和方位向两个维度进行谱峰扫描，从而确定(θk,ϕk)的值。 

3. 仿真代码 

       仿真环境：Matlab2021b

       阵列单元总数：Nx = 12,Ny = 8

       波源方向： (15,60)、(40,17)、(34,27)

       信噪比(SNR): 40 dB 

       仿真代码：

%% Author  ： Poulen
%% DAta    ： 2023.12.18

%%  简介：该代码主要分享如何实现2维DOA估计问题。截至目前为止，基于1维DOA估计问题的理论已经相对成熟，因此在这下面，基于MUSIC算法完成L型阵列目标的2维DOA估计，
%%        同时这也是做2维DOA估计的基础，希望可以为读者带来启示作用。

clear
close all;
clc

Nx  =   12;
Ny  =    8;
c   =  physconst('Lightspeed');
lambda = c/77e9;    %% 波长
deta   =  lambda/2;    %% 天线间隔
dx     = 0:deta:(Nx-1)*deta;
dy     = 0:deta:(Ny-1)*deta;
snapshot = 512;     %% 信号采样长度

source = [  15,60;
            40,17;
            34,27]; %% 波达方向(方位向，俯仰向)

source_num = size(source,1);
Ax = zeros(Nx,source_num);
Ay = zeros(Ny,source_num);
for i = 1:source_num
    Ax(:,i) = exp(1j*2*pi/lambda*dx(:)*cosd(source(i,1))*sind(source(i,2)));
    Ay(:,i) = exp(1j*2*pi/lambda*dy(:)*sind(source(i,1))*sind(source(i,2)));
end

S = randn(source_num,snapshot);  %产生随机信号
A = [Ax.' Ay.'].';               %合并x、y轴方向上的导向矢量
Z = A*S;                         %L型阵列接收的数据
% X = Ax*S;                      %%x轴方向阵列接收的数据
% Y = Ay*S;                      %%y轴方向阵列接收的数据
% Z = [X;Y];


%%  加入噪声
SNR = 40; %单位dB
Z = Z +(randn(size(Z)).*std(Z))/db2mag(SNR);

R = (1/snapshot)*Z*conj(Z.');
[V,D] = eig(R);             %在matalb2021版本以上，特征值已经从小到大排列
Un = V(:,1:end-source_num);          %提取小特征对应的特征向量作为噪声子空间 M*(M-K)

Azimuth = 0:0.1:60;
Pitch   = 0:0.1:90;
P2D = zeros(length(Azimuth),length(Pitch));
for n = 1:length(Azimuth)
    for m = 1:length(Pitch)
        ax = exp(1j*2*pi/lambda*dx(:)*cosd(Azimuth(n))*sind(Pitch(m)));
        ay = exp(1j*2*pi/lambda*dy(:)*sind(Azimuth(n))*sind(Pitch(m)));
        a  = [ax.' ay.'].';
        P2D(n,m) = 1/(a'*(Un*Un')*a);
    end
end
Pazimuth = sum(abs(P2D),2); 
Ppitch   = sum(abs(P2D),1);

%%  方位-俯仰向2D图
figure;
mesh(Pitch,Azimuth,db(P2D));
xlabel('Pitch(°)');
ylabel('Azimuth(°)');
title('方位-俯仰向图');

%%  绘制方位向剖面
figure;
plot(Azimuth,db(Pazimuth),'LineWidth',1.5);
title('目标方位角剖面');
xlabel('Azimuth(°)');
ylabel('幅度(dB)');

%%  绘制俯仰向剖面
figure;
plot(Pitch,db(Ppitch),'LineWidth',1.5);
title('目标俯仰角剖面');
xlabel('Azimuth(°)');
ylabel('幅度(dB)');

仿真结果：

          根据仿真结果，可以看出，不管是联合水平-方位向结果，还是水平向和俯仰向剖面，都得到了很好的估计效果。通过该例子，主要是为了提醒大家，要学会举一反三，不管是什么类型的阵列，主要找准他们的方向矩阵，并通过谱峰扫描就可以估计得到目标的水平-俯仰信息。另外，希望本期的内容可以帮助到迷茫的大家，谢谢！