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

阵列信号DOA估计系列(三).MVDR/Capon波束形成(附代码)本文主要介绍Capn波束形成算法,又名最小方差无失真相应(MinimumVarianceDistortionlessResponse,MVDR),并将其方法应用于DOA估计。

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

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

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

1. 波束形成

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

2. MVDR波束形成基本推导

在前面已介绍了阵列接收信号模型和导向矢量,在此假设均匀线阵(ULA)含有 M M M 个阵元,当有 N N N 个信号 s 1 ( n ) , s 2 ( n ) , ⋯   , s N ( 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 = A M × N s N × 1 = [ a ( θ 1 ) , a ( θ 2 ) , ⋯   , a ( θ N ) ] M × N × [ s 1 ( n ) , s 2 ( n ) , ⋯   , s N ( n ) ] 1 × N T (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 ) = w H x ( n ) (2) y(n)= \mathbf{w}^H\mathbf{x}(n)\tag2 y(n)=wHx(n)(2)同时据此可以计算出输出信号的功率为 P ≜ E [ ∣ y ( n ) ∣ 2 ] = E [ w H x ( n ) x H ( n ) w ] = w H R w (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}\tag3 PE[y(n)2]=E[wHx(n)xH(n)w]=wHRw(3)

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

条件1: w H a ( θ 1 ) = 1 (5) \mathbf{w}^H\mathbf{a(\theta_1)}=1\tag5 wHa(θ1)=1(5)
条件2:输出功率 P = x ( n ) x ( n ) H = w H R w {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 ⁡ w H R w , s t . w H a ( θ ) = 1 (6) \min \mathbf{w}^H\mathbf{R}\mathbf{w},st.\mathbf{w}^H\mathbf{a(\theta)}=1\tag6 minwHRwst.wHa(θ)=1(6)

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

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

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

3.MVDR波束形成计算步骤

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

4. 结论和思考

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

5.仿真结果

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

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

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

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

6.仿真代码

代码在此处,请点击下载

PS:如有错误,请大家不吝指正。谢谢!

版权声明:本文内容由互联网用户自发贡献,该文观点仅代表作者本人。本站仅提供信息存储空间服务,不拥有所有权,不承担相关法律责任。如发现本站有涉嫌侵权/违法违规的内容, 请发送邮件至 举报,一经查实,本站将立刻删除。

发布者:全栈程序员-用户IM,转载请注明出处:https://javaforall.cn/148861.html原文链接:https://javaforall.cn

【正版授权,激活自己账号】: Jetbrains全家桶Ide使用,1年售后保障,每天仅需1毛

【官方授权 正版激活】: 官方授权 正版激活 支持Jetbrains家族下所有IDE 使用个人JB账号...

(0)
blank

相关推荐

发表回复

您的电子邮箱地址不会被公开。

关注全栈程序员社区公众号