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

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

1. 波束形成

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

2. MVDR波束形成基本推导

在前面已介绍了阵列接收信号模型和导向矢量，在此假设均匀线阵（ULA）含有 MMM 个阵元，当有 NNN 个信号 s1(n),s2(n),⋯,sN(n)s_{_1}(n),s_{_2}(n),\cdots,s_{_N}(n)s1(n),s2(n),⋯,sN(n) 分别从 θ1,θ2,⋯,θN\theta_1,\theta_2,\cdots,\theta_Nθ1,θ2,⋯,θN 入射到阵列时，按照叠加的思维，接收信号可以表述为x(n)M×1=AM×NsN×1=[a(θ1),a(θ2),⋯,a(θN)]M×N×[s1(n),s2(n),⋯,sN(n)]1×NT(1)\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 x(n)M×1=AM×NsN×1=[a(θ1),a(θ2),⋯,a(θN)]M×N×[s1(n),s2(n),⋯,sN(n)]1×NT(1)
现在假设我们需要设计一种权值 w\mathbf{w}w，使得 θ1\theta_1θ1 方向的信号完全通过，其余的信号连同噪声被最大可能抑制。于是输出 y(n)y(n)y(n) 可以写作 y(n)=wHx(n)(2)y(n)= \mathbf{w}^H\mathbf{x}(n)\tag2y(n)=wHx(n)(2)同时据此可以计算出输出信号的功率为P≜E[∣y(n)∣2]=E[wHx(n)xH(n)w]=wHRw(3)P \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}\tag3P≜E[∣y(n)∣2]=E[wHx(n)xH(n)w]=wHRw(3)

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

条件1：wHa(θ1)=1(5)\mathbf{w}^H\mathbf{a(\theta_1)}=1\tag5wHa(θ1)=1(5)
条件2：输出功率 P=x(n)x(n)H=wHRw{P} =\mathbf{x}(n)\mathbf{x}(n)^H=\mathbf{w}^H\mathbf{R}\mathbf{w}P=x(n)x(n)H=wHRw最小。

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

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

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

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

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

３.MVDR波束形成计算步骤

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

4. 结论和思考

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

5.仿真结果

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

5.1 用MVDR波束形成方法进行DOA估计

从仿真结果可以看出，信号非相干时，此方法具有较高的分辨率；但当信号相干是，虽然还是在　10°,20°10\degree,20\degree10°,20° 方向依然有两个峰值，但是其对应的纵坐标较小，且其余地方也有峰值，这就给后续的检测算法带来了难度。

作为对比，也将空间FFT的结果放在这里，可以看到，MVDR波束形成方法的分辨率要高很多

6.仿真代码

代码在此处，请点击下载

PS：如有错误，请大家不吝指正。谢谢！

阵列信号DOA估计系列(三).MVDR/Capon波束形成(附代码)相关推荐

阵列信号DOA估计系列(二).导向矢量与空间FFT(附代码)
阵列信号DOA估计系列(二).导向矢量在DOA估计里面,经常会看到导向矢量这个名词,也有的地方叫方向矢量,方向矩阵,基本上都是array steering vector 的翻译. 本文首先对均匀线阵 ...
阵列信号DOA估计系列(一).概述
阵列信号DOA估计系列之概述 1.从相位差说起 2.空间相位差的来源 3.从"空间相位差"到"DOA估计" 3.1 时域 3.2 空域 3.3 DOA估计 ...
【压缩感知】基于matlab压缩感知理论的窄带信号DOA估计【含Matlab源码 2616期】
⛄一.压缩感知理论阵列信号波达方向(Direction ofArrival,DOA)估计是阵列信号处理领域中主要研究内容之一,广泛应用于军事及民用领域.基于压缩感知理论的稀疏重构算法的阵列信号DOA ...
基于matlab的相干信号的doa 估计,基于空间平滑MUSIC算法的相干信号DOA估计(1)
基于空间平滑MUSIC算法的相干信号DOA估计(1) 基于空间平滑MUSIC算法的相干信号DOA估计(1) 空间平滑MUSIC算法(1) 在上一篇博客中有提到,当多个入射信号相干时,传统MUSIC算法 ...
常见传统算法实现DOA估计总结CBF、Capon、MUSIC、ESPRIT、OMP
常见传统算法DOA估计总结 CBF算法传统时域傅里叶谱估计方法在空域中简单拓展形式,空间分辨能力会受到"瑞利限"的限制 Capon算法通过对与信号协方差矩阵以及阵列方向矢量相关 ...
宽带信号doa matlab,宽带信号DOA估计处理方法研究
学术研究 DOI:10． 3969/j． issn． 1001-3824． 2012． 06． 008 宽带信号 DOA 估计处理方法研究收稿日期:2012-06-14 闫杰1,周围1,2,杜晓 ...
基于空间平滑MUSIC算法的相干信号DOA估计(1)
空间平滑MUSIC算法(1) 1. 前言在上一篇博客中有提到,当多个入射信号相干时,传统MUSIC算法的效果就会不理想.具体原因是多个入射信号相干时,有部分能量就会散发到噪声子空间,使得MUSIC算 ...
基于空间平滑MUSIC算法的相干信号DOA估计(2)
空间平滑MUSIC算法(2) 继续上一篇博客,继续讲后向空间平滑和前/后向空间平滑MUSIC算法. 基于空间平滑MUSIC算法的相干信号DOA估计(1) 2.3 后向空间平滑算法后向空间平滑更准确的 ...
DOA估计基于互质阵列的DOA估计
前言传统阵列的配置方式是均匀线阵,该阵列要求相邻阵元的间距为半波长,易产生耦合效应,影响 DOA 估计精度.而稀疏阵列利用协方差矩阵构建差分共阵方式在虚拟域上生成虚拟阵列,并利用虚拟阵列实现波达方向 ...

阵列信号DOA估计系列(三).MVDR/Capon波束形成(附代码)