阵列信号处理（一）常规波束形成算法

1、问题描述

\qquad 阵列信号处理在雷达、声呐、无线通信、医学成像、地质勘探、射电天文学等多种领域具有广泛的应用。波束形成是阵列信号处理的一个重要的任务。常规波束形成算法（Conventional Beamforming, CBF）又称为时延求和（Delay And Sum, DAS）波束形成，其对各阵元数据通过简单的延迟求和以达到空间滤波的效果，具有最好的稳健性。

2、窄带时域快拍模型

\qquad 考虑单平面波输入的情况，假设在参考点接收的窄带信号 $$s(t)=\bar{s}(t)e^{j\omega t}$$其中，$\bar{s}(t)$ 为接收信号的复包络，$\omega$ 为载波角频率。在第 $m$ 个阵元接收的窄带信号表示为$$s_m(t)={\bar{s}}_m(t)e^{j\omega t}$$显然，第 $m$ 个阵元接收信号 $s_m(t)$ 是参考点信号 $s(t)$ 经过时间 ${\tau }_m$ 延迟得到的。因此，$$s_m(t)=s(t-{\tau }_m)=\bar{s}(t-{\tau }_m)e^{j\omega (t-{\tau }_m)}$$在窄带假设条件下，有$$\bar{s}(t-{\tau }_m)\approx \bar{s}(t)$$则 $$s_m(t)=\bar{s}(t)e^{j\omega (t-{\tau }_m)}=s(t)e^{-j\omega {\tau }_m}$$ \qquad 对于标准均匀线阵，时延 ${\tau }_m$ 可表示为$${\tau }_m=\frac{mdsin(\theta )}{c},m=0,1,\cdots ,M$$其中，$M$ 为阵元个数，$d$ 为阵元间距，$\theta$ 为入射角度，$c$ 为信号传播速度。此时，第 $m$ 个阵元接收信号 $s_m(t)$ 可表示为$$s_m(t)=s(t)e^{-j\omega \frac{mdsin(\theta )}{c}}=s(t)e^{-j2\pi f\frac{mdsin(\theta )}{c}}$$其中，$f$ 为载波频率（$\omega =2\pi f$）。
\qquad 将所有 $M$ 个阵元在 $t$ 时刻接收的信号排列成一个列矢量，$$\mathbf{x}(t)=s(t)[1,e^{-j2\pi f\frac{dsin(\theta )}{c}},\cdots ,e^{-j2\pi f\frac{(M-1)dsin(\theta )}{c}}{]}^T+\mathbf{n}(t)=s(t)\mathbf{a}(\theta )+\mathbf{n}(t)$$其中，$a(\theta )$ 称为阵列流形向量，$\mathbf{n}(t)=[n_1(t),\cdots ,n_M(t){]}^T$ 为基阵接收的 $M\times 1$ 维的噪声向量。
\qquad 对上式进行采样获得时域快拍模型，即$$\mathbf{x}(n)=s(n)\mathbf{a}(\theta )+\mathbf{n}(n),n=1,\cdots ,N$$式中，$N$ 为数据快拍长度。

3、算法描述

\qquad 波束形成一般指的是对基阵各阵元采集数据进行线性时不变滤波再求和，得到波束输出。在窄带情况下，可以采用一个复数加权代替线性时不变滤波。对窄带快拍数据进行加权求和，得到波束形成器得输出快拍为$$y(n)={\mathbf{w}}^H\mathbf{x}(n)$$其中，$\mathbf{w}=[w_1,w_2,\cdots ,w_M{]}^T\in {\mathbb{C}}^M$。常规波束形成算法的加权向量为$$\mathbf{w}=\mathbf{a}(\theta )/M$$ \qquad 下面再介绍两个性能指标，分别是波束响应和方位谱。波束响应指的是波束形成器对某方位单位功率平面波信号的响应，其只与加权向量有关，而与该方向是否存在信号无关。其定义如下：$$p(\theta )={\mathbf{w}}^{H}a(\theta ),\theta \in \Theta$$其中，$\Theta$ 为所有可能的信号到达方位的集合。对于标准均匀线阵，$\Theta =[-90^{\circ },90^{\circ }]$。
\qquad 根据波束响应能量相对于方位的函数，即 $20\text{lg}|p(\theta )|$，便可得到波束形成器的指向性图，也称为波束图。波束图反映的是基阵对不同方位到达信号响应情况，可以用来评估其他方向干扰与噪声对感兴趣方向信号产生的影响大小。
\qquad 方位谱也称为波束扫描方位谱，指的是波束观察方向扫描得到的波束输出功率相对于方位的函数。其定义为$$P(\theta )={\mathbf{w}}^H(\theta ){\mathbf{R}}_x\mathbf{w}(\theta ),\theta \in \Theta$$其中，${\mathbf{R}}_x$ 为数据协方差矩阵。在实际中，协方差矩阵是未知的，可以由一段数据快拍样本的空间相关矩阵来估计，即 $${\mathbf{R}}_x=\frac{1}{N}\sum _{n=1}^N[\mathbf{x}(n){\mathbf{x}}^H(n)]$$其中，$N$ 为数据快拍长度。
\qquad 根据方位谱可以绘制方位谱图，一般对方位谱取对数，即 $10\text{lg}|P(\theta )|$。方位谱图可以作为信号到达方位的估计值，体现的是某一信号源对其他方位输出功率的影响。

4、仿真结果

\qquad 假设一个均匀线阵有 256 个阵元，阵元间距为 0.5 m。取 6000 个快拍长度，源信号信号从 $10^{\circ }$ 方向入射基阵，信号传播速度为1500 m/s，频率为 500 HZ。利用常规波束形成算法可以得到下列波束图和方位谱图。

5、matlab 代码

clc
clear all
close all
simpleFrequency = 6e3;
t = 1/simpleFrequency:1/simpleFrequency:1;
snapNums = 6e3;
expectAngle = [10];
arrayNums = 256;
frequency = 500;
speed = 1500;
arraySpace = 0.5;
steeringVector = exp(-1i*2*pi*frequency*arraySpace*(0:arrayNums-1)'*sind(expectAngle)/speed);
Signal = exp(1j*2*pi*frequency*t);
snapData = steeringVector*Signal;
noise = 1/sqrt(2)*(randn(arrayNums, snapNums) + 1j*randn(arrayNums, snapNums));
snapData = snapData + noise;
R = snapData*snapData'./snapNums;
theta = [-90:1:90];
steerMatrix = exp(-1j*2*pi*arraySpace*frequency/speed*(0:arrayNums-1)'*sind(theta));

% 绘制CBF波束响应

beamResponse = steeringVector'*steerMatrix./arrayNums;
figure(1)
plot(theta,20*log10(abs(beamResponse)),'k')
title('波束图')
xlabel('角度/degree')
ylabel('波束/dB')
xlim([-90 90])
xticks([-90 -60 -30 0 30 60 90])
grid on
% 绘制CBF空间谱
Pcbf = zeros(1,length(theta));
for ii = 1:length(theta)
    Pcbf(ii) = steerMatrix(:,ii)'*R*steerMatrix(:,ii)/(arrayNums*arrayNums);
end
figure(2)
plot(theta,20*log10(abs(Pcbf)),'k')
title('方位谱图')
xlabel('角度/degree')
ylabel('方位谱/dB')
xlim([-90 90])
xticks([-90 -60 -30 0 30 60 90])

grid on

参考资料

[1] 鄢社锋,马远良.传感器阵列波束优化设计及应用[M].科学出版社,2009.