波束形成：最小方差无畸变响应波束形成器（MVDR）

MVDR波束形成器

1. 问题描述

  假设远场空间有一个感兴趣的信号$d\left(t\right)$，称为期望信号，其波达方向为${\theta }_d$，和$J$个不感兴趣的信号$i_j,\left(j=1,2,\cdots ,J\right)$，称为干扰信号，其波达方向分别为${\theta }_{ij}$，每个阵元上的加性白噪声为$n_m\left(t\right)$，则第$m$个阵元接收到的信号为
$$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)$$
对于整个阵列，写成矩阵形式为
$$\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)+\text{n}\left(t\right)$$
式中，$\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]$，表示来自波达方向${\theta }_k\left(k=d,i_1,i_2,\cdots \right)$的发射信源的方向向量。

  $N$个快拍的波束形成器输出$y\left(t\right)={\mathbf{w}}^H\left(\theta \right)\mathbf{x}\left(t\right)$的平均功率为
$$\begin{array}{rl}P\left(w\right)& =\frac{1}{N}\sum _{t=1}^N{\left|y\left(t\right)\right|}^2=\frac{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\frac{1}{N}\sum _{t=1}^N{\left|d\left(t\right)\right|}^2+\sum _{j=1}^J\left[\frac{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+\frac{1}{N}{\Vert \mathbf{w}\Vert }^2\sum _{t=1}^N{\Vert \text{n}\left(t\right)\Vert }^2\end{array}$$
这里忽略了不同用户之间的相互作用项，即交叉项$i_j\left(t\right)i_k^{\ast }\left(t\right)$，当$N\to \infty$时，上式可写为
$$\begin{array}{rl}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\text{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[{|i_j\left(t\right)\Vert }^2\right]{\left|{\mathbf{w}}^H\mathbf{a}\left({\theta }_{ij}\right)\right|}^2+{\sigma }_n^2{\Vert \mathbf{w}\Vert }^2\end{array}$$
式中，$\mathbf{R}=\mathrm{E}\left[\mathbf{x}\left(t\right){\mathbf{x}}^H\left(t\right)\right]$为阵列输出的协方差矩阵。为了保证来自方向${\theta }_d$期望信号的正确接收，并完全已知其他$J$个干扰，很容易根据上式得到关于权重的约束条件
$$\{\begin{array}{rl}& {\mathbf{w}}^H\mathbf{a}\left({\theta }_d\right)=1\\ & {\mathbf{w}}^H\mathbf{a}\left({\theta }_{ij}\right)=0\end{array}$$
该约束条件被称为波束“置零条件”，因为它强迫接收阵列波束方向图的“零点”指向所有$J$个干扰信号。在满足约束条件时，保证噪声最小，则优化目标函数为
$$\underset{w}{\min }\mathrm{E}\left[{\left|y\left(t\right)\right|}^2\right]=\underset{w}{\min }{\mathbf{w}}^H\mathbf{R}\mathbf{w}$$

  则波束形成问题可以表述为
$$\begin{array}{rl}& \underset{w}{\min }{\mathbf{w}}^H\text{Rw}\\ & \mathrm{s}.\mathrm{t}.\{\begin{array}{rl}& {\mathbf{w}}^H\mathbf{a}\left({\theta }_d\right)=1\\ & {\mathbf{w}}^H\mathbf{a}\left({\theta }_{ij}\right)=0\end{array}\end{array}$$

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

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

采用Lagrange乘子法，
$$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]$$
对其求导，并令其为$0$，则得到
$$\frac{\partial L\left(\mathbf{w}\right)}{\partial \mathbf{w}}=2\mathbf{R}\mathbf{w}+\lambda \mathbf{a}\left({\theta }_d\right)=0$$
可以解出
$$\mathbf{w}=\mu {\mathbf{R}}^{-1}\mathbf{a}\left({\theta }_d\right)$$
式中，$\mu$是比例常数。

  约束条件${\mathbf{w}}^H\mathbf{a}\left({\theta }_d\right)=1$也可以写为$\mathbf{a}{\left({\theta }_d\right)}^H\mathbf{w}=1$，则上式两遍同时乘以${\mathbf{a}}^H\left({\theta }_d\right)$，可以得到
$$\mu =\frac{1}{{\mathbf{a}}^H\left({\theta }_d\right){\mathbf{R}}^{-1}\mathbf{a}\left({\theta }_d\right)}$$

  则根据MVDR准则求得的最优权重为
$${\mathbf{w}}_{\mathrm{M}\mathrm{V}\mathrm{D}\mathrm{R}}=\frac{{\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)}$$

3. Matlab代码

对12阵元均匀线阵进行仿真，频率$f=10\mathrm{G}\mathrm{H}\mathrm{z}$，期望信号角度为${\varphi }_d=10^{\circ }$，干扰信号角度为${\varphi }_j=30^{\circ }$。

% 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');