波束形成:最小方差无畸变响应波束形成器(MVDR)

MVDR波束形成器

1. 问题描述

  假设远场空间有一个感兴趣的信号 d ( t ) d\left(t\right) d(t),称为期望信号,其波达方向为 θ d \theta_d θd,和 J J J个不感兴趣的信号 i j , ( j = 1 , 2 , ⋯   , J ) i_j,\left( j=1,2,\cdots,J\right) ij,(j=1,2,,J),称为干扰信号,其波达方向分别为 θ i j \theta_{ij} θij,每个阵元上的加性白噪声为 n m ( t ) n_m\left(t\right) nm(t),则第 m m m个阵元接收到的信号为
x m ( t ) = a m ( θ d ) d ( t ) + ∑ j = 1 J a m ( θ i j ) i j ( t ) + n m ( t ) x_m\left(t\right)=a_m\left(\theta_d\right)d\left(t\right)+\sum_{j=1}^{J}a_m\left(\theta_{ij}\right)i_j\left(t\right)+n_m\left(t\right) xm(t)=am(θd)d(t)+j=1Jam(θij)ij(t)+nm(t)
对于整个阵列,写成矩阵形式为
x ( t ) = a ( θ d ) d ( t ) + ∑ j = 1 J a ( θ i j ) i j ( t ) + n ( t ) \mathbf{x}\left(t\right)=\mathbf{a}\left(\theta_d\right)d\left(t\right)+\sum_{j=1}^{J}\mathbf{a}\left(\theta_{ij}\right)i_j\left(t\right)+\textbf{\textit{n}}\left(t\right) x(t)=a(θd)d(t)+j=1Ja(θij)ij(t)+n(t)
式中, a ( θ k ) = [ a 1 ( θ k ) , a 2 ( θ k ) , ⋯   , a M ( θ k ) ] \mathbf{a}\left(\theta_k\right)=\left[a_1\left(\theta_k\right),a_2\left(\theta_k\right),\cdots,a_M\left(\theta_k\right)\right] a(θk)=[a1(θk),a2(θk),,aM(θk)],表示来自波达方向 θ k ( k = d , i 1 , i 2 , ⋯   ) \theta_k\left(k=d,i_1,i_2,\cdots\right) θk(k=d,i1,i2,)的发射信源的方向向量。

   N N N个快拍的波束形成器输出 y ( t ) = w H ( θ ) x ( t ) y\left(t\right)=\mathbf{w}^H\left(\theta\right)\mathbf{x}\left(t\right) y(t)=wH(θ)x(t)的平均功率为
P ( w ) = 1 N ∑ t = 1 N ∣ y ( t ) ∣ 2 = 1 N ∑ t = 1 N ∣ w H x ( t ) ∣ 2 = ∣ w H a ( θ d ) ∣ 2 1 N ∑ t = 1 N ∣ d ( t ) ∣ 2 + ∑ j = 1 J [ 1 N ∑ t = 1 N ∣ i j ( t ) ∣ 2 ] ∣ w H a ( θ i j ) ∣ 2 + 1 N ∥ w ∥ 2 ∑ t = 1 N ∥ n ( t ) ∥ 2 \begin{aligned} P\left(w\right) & = \dfrac{1}{N}\sum_{t=1}^{N}\left| y\left(t\right)\right| ^2=\dfrac{1}{N}\sum_{t=1}^{N}\left| \mathbf{w}^H\mathbf{x}\left(t\right)\right| ^2 \\ & = \left| \mathbf{w}^H\mathbf{a}\left(\theta_d\right)\right| ^2\dfrac{1}{N}\sum_{t=1}^{N}\left| d\left(t\right)\right| ^2 +\sum_{j=1}^{J}\left[\dfrac{1}{N}\sum_{t=1}^{N}\left| i_j\left(t\right)\right| ^2\right]\left| \mathbf{w}^H\mathbf{a}\left(\theta_{ij}\right)\right| ^2 + \dfrac{1}{N}\left\|\mathbf{w}\right\|^2 \sum_{t=1}^{N}\left\| \textbf{\textit{n}}\left(t\right)\right\| ^2\\ \end{aligned} P(w)=N1t=1Ny(t)2=N1t=1NwHx(t)2=wHa(θd)2N1t=1Nd(t)2+j=1J[N1t=1Nij(t)2]wHa(θij)2+N1w2t=1Nn(t)2
这里忽略了不同用户之间的相互作用项,即交叉项 i j ( t ) i k ∗ ( t ) i_j\left(t\right)i^*_k\left(t\right) ij(t)ik(t),当 N → ∞ N\to\infty N时,上式可写为
P ( w ) = E [ ∣ y ( t ) ∣ 2 ] = w H E [ x ( t ) x H ( t ) ] w = w H Rw = E [ ∣ d ( t ) ∣ 2 ] ∣ w H a ( θ d ) ∣ 2 + ∑ j = 1 J E [ ∣ i j ( t ) ∥ 2 ] ∣ w H a ( θ i j ) ∣ 2 + σ n 2 ∥ w ∥ 2 \begin{aligned} P\left(w\right) & =\mathrm{E}\left[\left| y\left(t\right)\right| ^2\right]=\mathbf{w}^H\mathrm{E}\left[\mathbf{x}\left(t\right)\mathbf{x}^H\left(t\right)\right]\mathbf{w}=\mathbf{w}^H\textbf{\textit{Rw}} \\ & =\mathrm{E}\left[\left|d\left(t\right) \right|^2 \right]\left|\mathbf{w}^H\mathbf{a}\left(\theta_d\right) \right| ^2+\sum_{j=1}^{J}\mathrm{E}\left[\left|i_j\left(t\right) \right\| ^2\right]\left|\mathbf{w}^H\mathbf{a}\left(\theta_{ij}\right) \right| ^2+\sigma^2_n\left\| \mathbf{w}\right\|^2 \\ \end{aligned} P(w)=E[y(t)2]=wHE[x(t)xH(t)]w=wHRw=E[d(t)2]wHa(θd)2+j=1JE[ij(t)2]wHa(θij)2+σn2w2
式中, R = E [ x ( t ) x H ( t ) ] \mathbf{R} = \mathrm{E}\left[\mathbf{x}\left(t\right)\mathbf{x}^H\left(t\right)\right] R=E[x(t)xH(t)]为阵列输出的协方差矩阵。为了保证来自方向 θ d \theta_d θd期望信号的正确接收,并完全已知其他 J J J个干扰,很容易根据上式得到关于权重的约束条件
{ w H a ( θ d ) = 1 w H a ( θ i j ) = 0 \left\lbrace\begin{aligned} & \mathbf{w}^H\mathbf{a}\left(\theta_d\right)=1 \\ & \mathbf{w}^H\mathbf{a}\left(\theta_{ij}\right)=0 \\ \end{aligned}\right. {wHa(θd)=1wHa(θij)=0
该约束条件被称为波束“置零条件”,因为它强迫接收阵列波束方向图的“零点”指向所有 J J J个干扰信号。在满足约束条件时,保证噪声最小,则优化目标函数为
min ⁡ w E [ ∣ y ( t ) ∣ 2 ] = min ⁡ w w H R w \min_w \mathrm{E}\left[\left|y\left(t\right) \right|^2\right]=\min_w \mathbf{w}^H\mathbf{Rw} wminE[y(t)2]=wminwHRw

  则波束形成问题可以表述为
min ⁡ w w H Rw s . t . { w H a ( θ d ) = 1 w H a ( θ i j ) = 0 \begin{aligned} & \min_w \mathbf{w}^H\textbf{\textit{Rw}} \\ & \mathrm{s.t.} \left\lbrace\begin{aligned} & \mathbf{w}^H\mathbf{a}\left(\theta_d\right)=1 \\ & \mathbf{w}^H\mathbf{a}\left(\theta_{ij}\right)=0 \\ \end{aligned}\right. \end{aligned} wminwHRws.t.{wHa(θd)=1wHa(θij)=0

2. 最小方差无畸变响应波束形成器

  最小均方无畸变响应(Minimum Variance Distortionless Response, MVDR)波束形成器就是求解各个阵元权重的结果。

采用Lagrange乘子法,
L ( w ) = w H R w + λ [ w H a ( θ d ) − 1 ] L\left(\mathbf{w}\right) = \mathbf{w}^H\mathbf{R}\mathbf{w}+\lambda\left[\mathbf{w}^H\mathbf{a}\left(\theta_d\right)-1\right] L(w)=wHRw+λ[wHa(θd)1]
对其求导,并令其为 0 0 0,则得到
∂ L ( w ) ∂ w = 2 R w + λ a ( θ d ) = 0 \dfrac{\partial L\left(\mathbf{w}\right)}{\partial \mathbf{w}}=2\mathbf{R}\mathbf{w}+\lambda \mathbf{a}\left(\theta_d\right)=0 wL(w)=2Rw+λa(θd)=0
可以解出
w = μ R − 1 a ( θ d ) \mathbf{w}=\mu \mathbf{R}^{-1}\mathbf{a}\left(\theta_d\right) w=μR1a(θd)
式中, μ \mu μ是比例常数。

  约束条件 w H a ( θ d ) = 1 \mathbf{w}^H\mathbf{a}\left(\theta_d\right)=1 wHa(θd)=1也可以写为 a ( θ d ) H w = 1 \mathbf{a}\left(\theta_d\right)^H\mathbf{w}=1 a(θd)Hw=1,则上式两遍同时乘以 a H ( θ d ) \mathbf{a}^H\left(\theta_d\right) aH(θd),可以得到
μ = 1 a H ( θ d ) R − 1 a ( θ d ) \mu = \dfrac{1}{\mathbf{a}^H\left(\theta_d\right)\mathbf{R}^{-1}\mathbf{a}\left(\theta_d\right)} μ=aH(θd)R1a(θd)1

  则根据MVDR准则求得的最优权重为
w M V D R = R − 1 a ( θ d ) a H ( θ d ) R − 1 a ( θ d ) \mathbf{w}_{\rm{MVDR}}=\dfrac{\mathbf{R}^{-1}\mathbf{a}\left(\theta_d\right)}{\mathbf{a}^H\left(\theta_d\right)\mathbf{R}^{-1}\mathbf{a}\left(\theta_d\right)} wMVDR=aH(θd)R1a(θd)R1a(θd)

3. Matlab代码

对12阵元均匀线阵进行仿真,频率 f = 10 G H z f=10\rm{GHz} f=10GHz,期望信号角度为 ϕ d = 1 0 ∘ \phi_d=10^\circ ϕd=10,干扰信号角度为 ϕ j = 3 0 ∘ \phi_j=30^\circ ϕj=30

MVDR波束形成

% Title:    最小方差无畸变响应波束形成器,MVDR
% Author:   Xu Zhe
% Date:     2020-07-25
% Brief:    一维均匀线阵波束形成,MVDR波束形成器

% 系统初始化 Initialization
close all
clear
clc

% 参数 Parameters
M   = 12;               % 阵元数目
c   = 3e8;              % 光速
f   = 10e9;             % 频率
l   = c/f;              % 波长
d   = l/2;              % 阵元间隔
phi = [10,30]*pi/180;   % 信号角度
Kt  = 1;                % 目标信号数目
snr = 20;               % 信噪比
inr = 10;               % 干噪比
Ac  = 3601;             % 角度采样数

% 计算阵元位置
R = d * (0:M-1)';

% 计算方向向量
K = length(phi);        % 信号数目
A = zeros(M,K);         % 初始化方向矩阵 Steering Matrix
for k = 1:K
    for m = 1:M
        A(m,k) = exp(1j * 2 * pi * R(m,1)* sin(phi(k))/l);
    end
end

% 构造信号
Ts = 1024;              % 快拍数
Kj = K - Kt;            % 干扰信号数目
s = zeros(K,Ts);        % 信号
sd = zeros(1,Ts);       % 期望信号
for k = 1:Kt
    for i = 0:Ts-1
        s(k,i+1) = sqrt(10^(snr/10)) * exp(1j * 2 * pi * f * i/(2*Ts)); 
    end
    sd = sd + s(k,:);
end
sd = sd/Kt;
if Kj ~= 0
    for k = Kt+1:K
        for i = 0:Ts-1
            s(k,i+1) = sqrt(10^(inr/10)/2) * (randn(1) + 1j*randn(1)) * exp(1j * 2 * pi * f * i/(2*Ts));
        end
    end
end

% 接收数据向量
n = (randn(M,Ts)+1j*randn(M,Ts))/sqrt(2);
x = A * s+n;

% MVDR
J = A(:,Kt+1:K) * s(Kt+1:K,:);
Rnj = 1/Ts * (J*J'+n*n');
Rpriv = pinv(Rnj);
w = Rpriv*A(:,1)/(A(:,1)'*Rpriv*A(:,1));

% 方向图
theta = linspace(-pi,pi,Ac);
y = zeros(1,Ac);
for i = 1:Ac
    a = exp(1j * 2 * pi * R(:,1)* sin(theta(i))/l);
    y(i) = w' * a;
end

ydB = 20*log10(abs(y)/max(abs(y)));

% 绘图
figure('name','波束方向图','color',[1,1,1]);
set(gcf,'position',[458,342,290,200])
plot(theta * 180/pi,ydB,'color','b','linestyle','-','linewidth',1);
grid on
xlim([-90,90]);
ylim([-60,0]);
set(gca,'xtick',-90:30:90);
set(gca,'FontName','Times New Roman','FontSize',10.5);
xlabel('Direction of Arrival of the Signal (Degree)','FontSize',10.5);
ylabel('Magnitude Response (dB)','FontSize',10.5);

% 保存
% print('-depsc','E:\07 Latex\Paper Note\figures\MVDR_1.eps');
©️2020 CSDN 皮肤主题: 大白 设计师:CSDN官方博客 返回首页