大家好，又见面了，我是你们的朋友全栈君。

标题阵列信号DOA估计系列(三).MVDR/Capon波束形成

MVDR算法得基本思路是在频域/空间形成一个窄带滤波器，从此出发，可见MVDR不但对噪声有抑制作用，来对观察频率/角度之外的信号有抑制作用，所以MVDR的分辨率远高于常规的FFT/DBF算法。
算法的推导和理解都比较简单，所以下面稍微推导一下。需要代码的请直接拉到最后。

1. 波束形成

波束形成，可以简单的理解为加权，或者滤波也可以。基本解释可以从这个图来看在这里插入图片描述对各个阵元接收到的信号进行加权，权系数为 $\mathbf{w} = [w_{_0},w_{_1},\cdots,w_{_M-1}]^T$ ，使得输出 $y (n)$ 达到我们想要的目的。
比如说：只让 $30\degree$ 方向的信号进来，其他方向尽量抑制掉；或者，只是抑制 $30\degree$ 方向的信号。
当然，对于权系数的设计，是一个比较重要的话题，此处不作展开。有兴趣的朋友，可以自行查阅DBF相关的资料。

2. MVDR波束形成基本推导

在前面已介绍了阵列接收信号模型和导向矢量，在此假设均匀线阵（ULA）含有 $M$ 个阵元，当有 $N$ 个信号 $s_{_1}(n),s_{_2}(n),\cdots,s_{_N}(n)$ 分别从 $\theta_1,\theta_2,\cdots,\theta_N$ 入射到阵列时，按照叠加的思维，接收信号可以表述为 $\mathbf{x}(n) _{_{M\times 1}}= \mathbf{A_{_{M\times N}}s_{_{N\times 1}}}=[\mathbf{a(\theta_1)},\mathbf{a(\theta_2)},\cdots,\mathbf{a(\theta_N)}]_{_{M\times N}}\times [s_{_1}(n), s_{_2}(n),\cdots, s_{_N}(n)]^T_{_{1\times N}}\tag1$
现在假设我们需要设计一种权值 $\mathbf{w}$ ，使得 $\theta_1$ 方向的信号完全通过，其余的信号连同噪声被最大可能抑制。于是输出 $y (n)$ 可以写作 $\mathbf{w}^H\mathbf{x}(n)\tag2$ 同时据此可以计算出输出信号的功率为 $\triangleq E[|y(n)|^2]=E[\mathbf{w}^H\mathbf{x}(n)\mathbf{x}^H(n)\mathbf{w}]=\mathbf{w}^H\mathbf{R}\mathbf{w}\tag3$

将信号 $s_{_1}(n)$ 单独分离出来，其余部分写作 $\mathbf{z(n)}$ ，于是，上式变为 $\mathbf{w}^H\mathbf{x}(n) =\mathbf{w}^H\mathbf{a(\theta_1)}s_{_1}(n) + \mathbf{w}^H\mathbf{z(n)}\tag4$ 按照我们的预想，权向量要 $\mathbf{w}$ ，使得 $s_{_1}$ 方向的信号完全通过，同时，想要 $\mathbf{z(n)}$ 被最大可能抑制，就应该满足：

条件1： $\mathbf{w}^H\mathbf{a(\theta_1)}=1\tag5$
条件2：输出功率 $=\mathbf{x}(n)\mathbf{x}(n)^H=\mathbf{w}^H\mathbf{R}\mathbf{w}$ 最小。

条件1很好理解。条件2的可以这样来理解，在条件1已经满足的前提下，如果其他信号和噪声都被抑制了，那么接收的信号功率自然就是最低的了，因此条件2是 $\mathbf{z(n)}$ 被最大可能抑制的直接结果。

把 $\theta_1$ 推广到任意角度 $\theta$ ，那么上面所有的东西，可以写成这么一个以权向量 $\mathbf{w}$ 为变量的优化问题 $\min \mathbf{w}^H\mathbf{R}\mathbf{w}，st.\mathbf{w}^H\mathbf{a(\theta)}=1\tag6$

运用拉格朗日乘子法，构造代价函数 $J(\mathbf{w})$ 为 $J(\mathbf{w}) = \mathbf{w}^H\mathbf{R}\mathbf{w} – \lambda(1-\mathbf{w}^H\mathbf{a(\theta)})\tag7$ 以权向量 $\mathbf{w}$ 为变量求代价函数 $J(\mathbf{w})$ 的梯度并令其等于 $0$ ，即$ $\nabla J(\mathbf{w}) = 2\mathbf{R}\mathbf{w}-2\lambda\mathbf{a(\theta)}\triangleq0\tag8$ 可以解得 $\mathbf{w} = \lambda\mathbf{R}^{-1} \mathbf{a(\theta)}\tag9$ 两端取共轭转置(H)，并右乘 $\mathbf{a(\theta)}$ ，同时注意到约束条件 $\mathbf{w}^H\mathbf{a(\theta)}=1$ ，可以得到 $\lambda =\frac{1}{\mathbf{a(\theta)}^H(\mathbf{R}^{-1})^H\mathbf{a(\theta)}}\tag{10}$ 将 $(10)$ 带入 $(9)$ ，即可解得的权向量为 $\mathbf{w} = \frac{\mathbf{R}^{-1} \mathbf{a(\theta)}}{\mathbf{a(\theta)}^H(\mathbf{R}^{-1})^H\mathbf{a(\theta)}}\tag{11}$

到此，权向量就计算出来了。回望一下，这个权向量的目的：仅使得 $(11)$ 式中给定的 $\mathbf{a(\theta)}$ 对应方向的信号通过，其余方向的信号和噪声最大程度的抑制。这时候，输出的平均功率由 $(3)$ 式给定，稍微计算一下可得 $P_{\theta}=\frac{1}{\mathbf{a(\theta)}^H\mathbf{R}^{-1}\mathbf{a(\theta)}}\tag{12}$

$(12)$ 式可以这样理解：
任意给一个角度，就有其对应的导向矢量 $\mathbf{a(\theta)}$ ，如果这个方向有信号，那么 $(12)$ 式计算出来的就是此方向信号平均功率，其值应该较大；如果没有信号只有噪声，那么计算出来的就是噪声平均功率，其值应该较小。
因此，将 $\theta$ 在自己想观察的角度扫描一圈，每个地方都计算一遍功率，选出其中较大的值，其对应的 $\theta$ 就是DOA估计结果。

３.MVDR波束形成计算步骤

Step1：由接收到的快拍信号 $\mathbf{x}(n)$ 估计自相关矩阵 $\mathbf{R}$ ;
Step2：计算自相关矩阵 $\mathbf{R}$ 的逆矩阵 $\mathbf{R}^{-1}$ ;
Step3：按照阵列几何形状，构造对应的导向矢量 $\mathbf{a(\theta)}$ ;
Step4：使 $\theta$ 按照一定的步进，在自己想观察的角度扫描，逐次计算 $P_{\theta}$ ;
Step5：对 $P_{\theta}$ 进行谱峰搜索，找出峰值点对应的 $\theta$ ;

4. 结论和思考

MVDR波束形成方法只能处理非相干信号。
在对 $(8)$ 式的求解中，对自相关矩阵 $\mathbf{R}$ 进行了求逆运算。这就要求 $\mathbf{R}$ 满秩，即信号之间是不相干的。若存在相干信号，那么上面的推导，到 $(8)$ 式就无法继续进行了。
那么，如果信号相干又该怎么办？
MVDR波束形成具有通用性，而不限于线阵。
从通篇推导可以看出，没有应用到 $\mathbf{a(\theta)}$ 的具体结构。对于其他形式的阵列，修改 $\mathbf{a(\theta)}$ 的形式即可；
用MVDR波束形成方法进行DOA估计，不用知道信源个数。MUSIC、ESPRIT算法等均需要进行信源个数估计；
用MVDR波束形成方法进行DOA估计，分辨率比空间FFT高得多，这一点从下面的仿真中可以看出来。

5.仿真结果

假设一个均匀线阵有16阵元， $\lambda/2$ 布阵；取1024个快拍估计自相关矩阵 $\mathbf{R}$ ，两个信号分别从 $10\degree,20\degree$ 方向入射大阵列，信噪比均为 $10 d B$ 。取信号相干和不相干的情况，用本文所述的MVDR波束形成方法和空间FFT及进行DOA估计，结果如下。