均匀圆阵UCA到虚拟线阵变换的DOA估计

SignalProcessingNWU

已于 2025-01-05 22:35:53 修改

阅读量1.3k

点赞数 22

文章标签：信号处理

于 2025-01-05 22:32:27 首次发布

本文链接：https://blog.csdn.net/qq_51604300/article/details/144931595

版权

1.引言

参考文章：Eigenstructure techniques for 2-D angle estimation with uniform circular arrays .https://ieeexplore.ieee.org/document/317861

梁国龙,田蕴琦,付进,等.基于均匀圆阵中心对称性的相干源方位估计[J].哈尔滨工程大学学报,2019,40(12):1958-1964.

https://kns.cnki.net/kcms2/article/abstract?v=2G5Y6l-6Cy_-_pwwqunUS3qQyxehVS9kDpysxfwHCyr0ouhpSr2a7GWs4sk1Kc5rW0Yn0FGTG71_-mjUT8u0pFwBv1L-iWitoYHRSFC_MzeOwD0FfmQRFw08QgC8LeNHvWITbsN6JEhNAmWN8WVUePC0WpHoPXlSiPTAQN0jp0-j0xn0NKX0Vtq-UsSyS4nh&uniplatform=NZKPT&language=CHS

均匀圆阵的阵列流形不具有Vandermonde结构性质，无法使用ULA阵列的一些DOA算法。本文学习虚拟阵变换法，通过模式空间变换将均匀圆阵转换为均匀线阵，其阵列流形转换成具有范德蒙结构。

2.UCA阵列结构

文章把入射信源与阵元看作在同一平面上，即俯仰角为0°，方位角 $\Theta$ 。 UCA阵半径为R，阵元数M，窄带信号数N，快拍数为L。

第m个阵元与信源的夹角： $\Theta -\Theta _{m}=\Theta -\frac{2\pi (m-1)}{M}$

阵列导向矢量为： $a=[e^{-j\omega cos(\Theta -\Theta _{1})r/c},...,e^{-j\omega cos(\Theta -\Theta _{m})r/c}]^{T}$

t时刻第m个阵元的接收信号表达式为： $x_{m}=A(\Theta )x(t)+n(t)$

其中x为窄带信号，n为噪声。

3.模式空间变换

图源自参考论文

虚拟均匀线阵如图，d为阵元间距，虚拟线阵的输出信号Y(t)由原UCA接收信号经转换矩阵T相乘得到。

$Y(t)=TX(t)=A_{v}\cdot s(t)+T\cdot n(t)$

其中 $A_{v}$ 是经转换后的虚拟线阵导向矢量。

模式空间变换矩阵 $T=\frac{1}{M}J^{-1}F$

F是离散傅里叶变换子矩阵，J是由k阶第1类贝塞尔 $J_{k}(\cdot )$ 构成的对角阵。

$J=diag\left \{ j^{-k}J_{-k} (-\beta ),...,j^{k}J_{k} (-\beta )\right \}$

$\beta =2\pi r/\lambda$

$F=[\omega _{-k},\omega _{-k+1},...,\omega _{k-1},\omega _{k}]^{H}$

$\omega _{k}=[1,e^{-j2\pi k/M},..,e^{-j2\pi k(M-1)/M}]$

那么带代入上式，虚拟线阵阵列流形为：

$a_{nv}=[e^{-jk\Theta _{n}},e^{-j(k-1)\Theta _{n}},...,e^{j(k-1)\Theta _{n}},e^{jk\Theta _{n}}]^{T}$

$A_{v}=[a_{(\Theta _{1})},a_{(\Theta _{2})},...,a_{(\Theta _{N})}]$

转换后的虚拟线阵阵列流形具有范德蒙结构。

此时可以算出虚拟线阵输出信号Y(t)。

由文章推导

UCA可激励最大模式数满足 $2K< M$ 。

4.MATLAB仿真实验

实验加入了俯仰角估计，使用MUSIC算法进行虚拟线阵DOA估计。

clc;
clear all;
M=24;
c = 3*10^8;%光速
f = 3*10^9;%接受信号频率
lamda = c/f;%波长
theta=[22]*pi/180;
phi=[60]*pi/180;
P=length(theta);
kp=150;
% R=lamda*sin(pi-pi/2/M)/sin(2*pi/M)/2;
R = lamda/(4*sin(pi/M));%不模糊最大阵列半径
%R=lamda*2/1;
beita=2*pi*R/lamda;
% N=floor(beita);
N=floor(M/2-1);
PN=2*N+1;
dtheta = 0.01; % 方位角搜索步长
dphi = 0.01; % 俯仰角搜索步长
%%%%%%%%%%%%%%%%%%%%生成预处理矩阵
% for in=1:PN
%     for ip=1:P
%        Aba(in,ip)=exp(i*(in-N-1)*theta(ip))*sin(phi(ip));
%     end
% end
fain=2*pi*[0:M-1]/M;
for ia=1:P
    for ifai=1:M
        a_r(ifai,ia)=exp(i*beita*cos(theta(ia)-fain(ifai))*sin(phi(ia)));
        %a_r(ifai,ia)=exp(i*beita*cos(theta(ia)-fain(ifai)));
    end
end
m=-N:N;
D=diag(i.^m.*besselj(m,beita));
for fi=1:PN
      for mi=1:M
          F(mi,fi)=exp(-i*2*pi*(fi-N-1)*(mi-1)/M);   
      end
end
T=pinv(D)*F'/M;%预处理矩阵

%shengcheng xinhao
fw=[0.01,0.05]*pi;
for is=1:1:P
    for t=1:1:kp
         if randn>0%%%%%%%%%%%%%产生不相干信号
             S(is,t)=cos(fw(is)*t)*sin(phi(is));
         else S(is,t)=-cos(fw(is)*t)*sin(phi(is));
         end
            %%%%%%%%%%%%产生相干信号
%              S(is,t)=cos(fw(is)*t);       
    end
end
zn=1:PN;
P_noise=0.1;
noise=sqrt(P_noise/2)*(randn(M,kp)+i*randn(M,kp));
Av=T*a_r;
%y=T*a_r*S+T*noise;  %预处理信号
y1=a_r*S+noise;
y=T*y1;
% MUSIC 算法估计方位角和俯仰角
Rxx = y * y' / kp; 

% 特征值分解
[V,D] = eig(Rxx);
[lambda_sorted, index] = sort(diag(D), 'descend');
V_sorted = V(:, index);
U_S = V_sorted(:, 1:P);
U_N = V_sorted(:, P+1:end);

% MUSIC算法搜索角度
for k = 1:ceil(2*pi/dtheta)
    theta1 = (k-1)*dtheta;
    for m = 1:ceil(pi/(2*dphi))
        phi1 = (m-1)*dphi;
        for n = 1:M
            phi_n = (2*pi*R/lamda)*cos(theta1 - fain(n))*sin(phi1);
            a1(n) = exp(-1i*phi_n);
        end
        a=T*a1';
        P_MUSIC(k,m) = 1/(a'*U_N*U_N'*a);
    end
end

% 找到谱峰位置
[max_value, index] = max(P_MUSIC(:));
[index_theta, index_phi] = ind2sub(size(P_MUSIC), index);
theta_est = (index_theta-1)*dtheta*180/pi;
phi_est = (index_phi-1)*dphi*180/pi;

disp(['估计的方位角：', num2str(theta_est)]);
disp(['估计的俯仰角：', num2str(phi_est)]);
disp(['真实的方位角：', num2str(theta*180/pi)]);
disp(['真实的俯仰角：', num2str(phi*180/pi)]);

估计结果