CN101098179A - 一种宽带频域数字波束形成方法 - Google Patents

Abstract

本发明涉及一种宽带声信号的数字波束形成方法，该方法包括正交采样、傅立叶变换、复乘移相和波束加等步骤，其中复乘移相步骤中的移相因子由两部分组成：第一部分用于完成正交波束形成中所需的相位旋转，第二部分用于完成复包络信号的波前对齐。本发明的波束形成方法着重于频域精确时延补偿的实现，利用分数点DFT循环移位方法可以对宽带信号进行精确延时，这种延时精度与信号的采样率无关，从而可以最大限度降低采样频率，减少冗余数据，大大提高波束形成实现的效率和减少系统复杂度。

Description

一种宽带频域数字波束形成方法

技术领域

本发明涉及数字信号处理领域的波束形成技术，特别涉及宽带声信号的数字波束形成方法。

背景技术

波束形成是广泛应用于雷达、声纳、超声成像以及无线通讯中的一项技术；波束形成器的主要功能是将传感器阵列接收信号加以组合得到来自某一具体方向的信号分量，或通过辐射阵列将发射信号约束在某一空间角度范围内定向发射；这里考虑接收端波束形成过程。

基本的波束形成方法通过对各阵元接收信号进行一定的延迟处理，以补偿某一方向入射波因阵元空间位置不同带来的时延差，并将延时后的数据进行简单加和(或增加幅度加权以改善旁瓣性能)即可得到指向该方向的波束输出(参见图3a)。在数字系统中，这种“延时-加和”的波束形成方法通过序列的移位迭加实现，而延时的精度与数据采样率成正比。为了得到良好的波束形成性能(较低的波束副瓣)需要对阵元序列进行精确的延时补偿。一般情况下，为了得到良好的波束形成性能，系统的采样率需要约10倍的奈奎斯特频率(即采样定理所要求的信号采样频率)。随着信号频率的提高和阵列规模的扩大，波束形成系统的复杂度(包括数据存储、传输及运算速率等)将急剧上升，这是制约数字波束形成系统实现的主要障碍。

现有的波束形成方法包括时域波束形成和频域波束形成两大类，时域法除基本的“时延-迭加”方法外，还包括“部分-加和”、内插方法以及相移边带波束形成方法等；频域方法基于频域相移等效于时域延迟这一事实，包括基于离散傅立叶变换的应用、相移波束形成以及补零2维FFT近似方法等。如何根据实际应用场合，选择合适的时延补偿方法是波束形成系统设计的关键。

对于主动声纳，声成像和水声通讯等带通系统，最常用的方法为正交采样波束形成方法，包括四个基本的处理步骤：

1)正交采样过程：将带通波形变换到复基带，得到复包络信号(可以在基带或中频)；

2)相位旋转(Phase Rotation)：对一定波束方向，补偿频谱搬移(正交采样)所产生的相位偏差；

3)时延复包络信号，补偿复信号波前延迟；

4)迭加各通道信号，得到波束输出(复包络)。

其中，第2步和第3步的先后次序可以颠倒。这一方法可以将高频带通信号变换到基带进行，由于复基带信号的奈奎斯特频率等于信号的带宽，对于信号中心频率与带宽之比较大的系统，可以显著降低波束形成系统的采样率和运算复杂度。

基于第2、3步的具体实现，现有技术有：

a)相移波束形成

这是现有的雷达，通讯系统中普遍采用的技术，对于窄带系统，信号复包络在阵列尺度上不会发生明显变化，这时可以去掉上述第3步过程简化波束形成计算。对于宽带信号，这种方法将使得时域信号发生畸变和动目标方位估计的偏差，如对短脉冲，其信号畸变的过程可参见图3b。因此，对于现代宽带成像系统和高速宽带通讯系统，这种方法已不能满足需求。

b)直接延时方法

直接对复包络序列进行移位延时补偿(美国专利US4290127)可以改善波束波形，然而复包络信号波前补偿精度仍由采样速率决定，为了满足宽带波束形成性能的要求，采样率仍需要数倍于基带奈奎斯特频率。为了降低采样率，减少冗余数据，现有技术采用时域内插方法改善延时精度，这样对内插后的数据序列进行波束形成仍然需要在较高的频率上进行。

c)基于离散傅立叶变换(DFT)的频域时延法

通过对频域序列乘上一定的相移因子，等效于信号的时域延时，基本的处理过程参见文献(Mucci，R.A.，AComparison of Efficient Beamforming Algorithms，IEEETrans.ASSP，Vol.32，No.3，pp548-558，June 1984.)。这种方法的时延操作精度与采样率无关，可以极大减少冗余数据，而且算法结构简单(仅包括FFT和复乘运算)，易于实现大规模并行处理。但是，这种方法基于DFT循环移位定理，在分数点移位条件下不能得到正确的时域输出。因此，这种方法通常用于只需要对波束频域输出进行处理，或对时域波形输出的连续性和精度要求不高的场合。

d)2维傅立叶变换频域波束形成方法

这种方法首先对每一通道应用FFT，然后对每一频率的各通道分量应用FFT，相当于2维FFT，具体过程可参见文献(Weber，M.E.and Heisler，R.，A.Frequency-domain beamforming algorithm for wideband coherent signal processing，J.Acoust.Soc.Amer.，Vol.76，No.4，pp 1132-1144，1984.)。对单频信号，空间场的FFT可得到远场波束特征；在宽带信号条件，不同频率的一定波数分量对应着空间不同方位，因此这种方法需要通过在波数域内插以得到近似波束输出，通常要求空域FFT内插倍数为8-16倍，同时这种频域2维FFT方法隐含了传感器阵元必须在空间上均匀分布。显然，这种方法结构复杂，只在阵元数目很大的条件下具有运算速度的优势。

近年来，现代声成像、高速水声通讯以及被动声定位等技术领域为了得到更高的性能(高成像分辨率，高数据传输率等)，普遍采用宽带信号体制。在布阵形式上，更倾向于采用稀疏阵、曲面阵等非规则阵型，例如稀疏阵形的引入可以使阵列规模极大下降，降低系统成本。在对波束输出的后处理上，更加强调波束输出波形的质量，如水声通讯和语音处理中的说话人声音的跟踪和识别等。对这些新的应用需求，现有的波束形成方法，包括一些近似技术，在系统复杂度或性能上已难以满足实际需求。

发明内容

本发明的目的是为了克服现有的数字波束形成方法中冗余数据多，精确延时代价高的缺点，利用分数点DFT循环移位方法对宽带信号进行精确延时，从而提供一种宽带频域数字波束形成方法，该方法在保证波束形成精确时延补偿的前提下，可以简化信号处理系统在数据传输，存储量，运算量和算法等方面的复杂度。

为了实现上述目的，本发明提供的宽带频域数字波束形成方法包括如下步骤：

1)对各阵元的接收信号进行前置放大；

2)在信号处理系统中，对前置放大后的信号进行正交采样得到复包络信号，该正交采样过程包括前置滤波和次级放大，高速模数转换，正交解调，低通滤波以及数据重采样；

3)对多通道的连续数据流进行分块处理，将数据流分解为相互重叠，大小为N_E×N_F的数据块，其中符号N_E表示阵元数，N_F表示进行傅立叶变换的数据长度；

4)对每个数据块沿行方向进行FFT变换，得到各通道数据的频域形式X(n，k)，其中n为阵元标号，k为频率标号；

5)计算各阵元在不同方向上各频率的相移因子；将不同方向的波束编号，成像区域的波束方向θ满足sinθ＝0，

，...，

...；其中m为波束编号(m≤Z)，Z为正整数，则对第m号波束，其相移因子P(n，k，m)＝(η(n，k))^m；

$\mathrm{\η}(n,k)=\exp \{-j\frac{2\mathrm{\πd}\left(n\right)}{\mathrm{cZ}}\·(f_c-\frac{k^{\′}}{F_s{N}_F})\};$

其中： $k^{\′}=\left\{\begin{array}{cc}k& k=1,...,N_F/2;\\ k-N_F& k=(N_F/2)+1,...,N_F;\end{array}\right.$

式中，n为阵元标号，k为频率标号，d(n)为第n阵元相对于阵中心的距离，c为水中声速，f_c为阵元接收信号的中心频率，F_s为正交采样后的复包络信号的采样率，N_F为傅立叶变换的数据长度；

6)利用步骤5)得到的相移因子与步骤4)得到的频域数据块进行复数乘法运算，得到各波束方向上的完成补偿的频域数据，对第m号波束，补偿后的频域数据

X_m′(n，k)＝X(n，k)×P(n，k，m)；

7)对步骤6)得到的结果进行加和得到该波束方向的频域输出Y(k)，对第m号波束，其频域输出

$Y_m\left(k\right)=\underset{n=1}{\overset{N_E}{\mathrm{\Σ}}}{X}_m^{\′}(n,k),$  k＝1,2,...,N_F；m＝1,2,...,N_B；

其中N_B为波束数目；

8)判断是否所有的波束形成完毕，若判断为否，则回到步骤5)，进行下一波束的相移因子和波束形成计算，若判断为是，则进入步骤9)；

9)对所有频域波束输出数据进行逆快速傅立叶变换(IFFT)，拼接相邻数据块，得到连续时域波束输出。

上述技术方案中，所述步骤2)中前置滤波的滤波通带选择1.5-4倍的阵元接收信号带宽。

上述技术方案中，所述步骤2)中的模数转换由同步时钟统一控制。

上述技术方案中，所述步骤2)中，每一通道的正交解调均由2个乘法器完成，其参考周期序列为sin(-2πf_cn/F_B)和cos(2πf_cn/F_B)，其中f_c为带通信号中心频率，n为整数，F_B为模数转换采样率。

上述技术方案中，所述步骤2)中，数据重采样后的复包络信号采样率为1.2～2倍奈奎斯特频率(复包络信号的奈奎斯特频率即为该复包络信号的带宽)。

上述技术方案中，所述步骤5)中，在存储器中存储中间参数η(n，k)，通过累乘运算实时计算获得相应波束的相位因子P(n，k，m)。

本发明的波束形成方法着重于频域精确时延补偿的实现，利用分数点DFT循环移位方法可以对宽带信号进行精确延时，这种延时精度与信号的采样率无关，从而可以最大限度降低采样频率，减少冗余数据，大大提高波束形成实现的效率和减少系统复杂度。

附图说明

图1是本发明一个实施例的水声成像系统示意图；

图2是本发明方法的结构框图；由正交采样、傅立叶变换、复乘移相和波束加和等几个部分组成；

图3a为基本的时延波束形成方法示意图，接收信号波前差异被延迟线补偿；

图3b是相位补偿方式实现波束形成的示意图；

图4是分数点DFT移位方法数值仿真结果；图中延时量为0.8个采样周期；

图5是本发明一个优选实施例中系统前级解调部分框图；由前级调理、高速AD、正交采样以及数据传输部分组成；

图6是本发明一个优选实施例的波束形成流程图；首先将上传数据分块，然后变换至频域进行移相复乘，波束加和，最后变换到时域得到波束输出；

具体实施方式

下面结合附图和具体实施方式对本发明作进一步详细描述。

实施例：

本发明的波束形成方法可用水声成像系统来实现，该系统可以采用现有技术中主动式声扫描工作方式。图1示例性地给出了该系统的结构示意图，如图1所示，该水下探测系统包括主动声源1和接收阵2，主动声源发射探测脉冲(在垂直方向扫描)，接收阵为一等间隔均匀线阵，测量的回波信号通过信号线传输至信号处理系统进行波束形成和时域处理完成水下目标的成像。为了有效区分水下目标，需要波束形成系统具有良好的性能(波束旁瓣应小于-20dB)。图2是本发明方法的结构框图；由正交采样、傅立叶变换、复乘移相和波束加和等几个部分组成。

图6示出了本发明波束形成方法的一个优选实施例流程，具体步骤如下所述：

在步骤101中，首先对声阵阵元输出信号进行前置放大后输入信号处理系统。

在步骤102中，对输入回波信号进行正交采样，现有技术中，正交采样可以通过模拟方式或数字方式实现，考虑到模拟方式在通道间有2°～3°的相位偏差，影响成像质量，这里采用数字方式。图5给出数字正交采样的实现过程，由前置滤波和次级放大、高速模数转换，正交解调，低通滤波以及数据重采样四个部分组成。接收信号首先进行前置滤波，滤波通带选择1.5-4倍的信号带宽，这样各通道滤波器的相位一致性较好。然后对各通道信号进行同步模数转换(ADC)，采样率为F_B赫兹，模数转换由同步时钟统一控制。每一通道的正交解调部分由2个乘法器组成，相乘的参考周期序列为sin(-2πf_cn/F_B)和cos(2πf_cn/F_B)，其中f_c为阵元接收到的带通信号的中心频率，n为整数。显然，参考周期序列可通过循环查表方式获得，这样频谱搬移的精度与模数转换同步时钟源精度相同；通过数字低通滤波器(LPF)后以N_d取1的方式对同相分量I和正交分量Q进行降采样抽取，降采样的倍数N_d的选择应当使复包络信号采样率接近其奈奎斯特采样率(一般可取1.2～2倍奈奎斯特频率)。这样正交采样后复信号采样率为F_s＝F_B/N_d。这一部分算法结构规整，可多路并行处理，尤其适合在大规模FPGA阵列中通过定点乘法实现。

在步骤103中，由于频域处理以数据块方式进行，需要对多通道的连续数据流进行分块处理。将数据流分解为相互重叠，大小为N_E×N_F的数据块，其中符号N_E表示阵元数，N_F表示需要进行傅立叶变换的数据长度。最大时延点数在最大波束偏角θ_MAX时得到

N_T＝[sin(θ_MAX)·L·F_s/c]

其中L表示水听器线阵长度，c为水中声速。

为了减小因循环移位过程中数据段边缘引起的冲激，需要对数据段两端进行加窗，窗的形式可以取半个哈明窗，窗长度N_win可取为8-16。

因此，数据块与块之间相互重叠部分的长度应大于等于N_T+N_win。

在步骤104中，直接对每个数据块沿行方向(本实施例中数据块的行表示不同的通道，列表示采样时间序列)进行FFT变换，得到各通道数据的频域形式，这里以X(n，k)表示，其中n为阵元标号，k为频率标号。

在步骤105中，需要计算对某一波束方向对应的相移因子，可以采用两种方法得到相移因子，一种方法是存储所有波束的相移因子，这可以使波束形成过程的运算量最小，但缺点是存储量大，对于阵元及需要生成波束数目众多的系统，直接存储相移因子的方法难以实现；第二种方法需要实时计算这些相移因子，如何快速地计算该相移因子是本发明的一个关键。首先分析移相因子的构成：

该相移因子由两部分组成：第一部分用于完成正交波束形成中所需的相位旋转，其公式为：p₁(n，θ)＝exp[- j2πf_cd(n)sinθ/c]

其中，θ为波束方向，d(n)为第n阵元相对于阵中心的距离。

第二部分用于完成复包络信号的波前对齐(时间延迟)，根据分数点DFT移位定理(该定理为本发明对整数点DFT移位定理的推广，在下文中有相应叙述)，可以得到这部分的相移因子为：

$p_2(n,k,\mathrm{\θ})=\left\{\begin{array}{cc}\exp \{j2\mathrm{\πkd}\left(n\right)\sin \mathrm{\θ}/c{F}_s{N}_F& k=1,...,N_F/2;\\ \exp \{j2\mathrm{\π}(k-N_F)d\left(n\right)\sin \mathrm{\θ}/c{F}_s{N}_F& k=(N_F/2)+1,...,N_F;\end{array}\right.$

总的相移因子可以写为：

P(n，k，θ)＝p₁(n，θ)p₂(n，k，θ)＝exp{sinθ·ξ(n，k)}

           ＝(exp{ξ(n，k)})^sinθ

令成像区域的波束方向这样选择，使得sinθ＝0，

，...，

...；这里m为波束编号，生成波束数目为2Z+1。其中Z为决定波束在空间方向上的密集程度的整数，Z越大则波束分布越密集，反之则越稀疏；通常，Z的取值范围在0.5×N_E到2×N_E之间。

总的相移因子可以以乘幂形式给出：

P(n，k，m)＝((exp(ξ(n，k))^1/Z)^m

          ＝(η(n，k))^m

$\mathrm{\η}(n,k)=\exp \{-j\frac{2\mathrm{\πd}\left(n\right)}{\mathrm{cZ}}\·(f_c-\frac{k^{\′}}{F_s{N}_F})\}$

其中： $k^{\′}=\left\{\begin{array}{cc}k& k=1,...,N_F/2\\ k-N_F& k=(N_F/2)+1,...,N_F\end{array}\right.$

这样需要存储的系数为η(n，k)，通过累乘可以得到需要波束的相位因子；注意到有P(n，k，m)＝P^*(n，k，-m)

其中符号(^*)表示复数共扼，这样可以节省一半相移因子的计算。

综上所述，相移因子的计算可以按波束编号通过累乘运算得到：首先计算m＝0时的0号波束，所需的相移因子P(n，k，0)≡1；完成波束形成后计算m＝1时的1号波束，所需的相移因子P(n，k，1)＝η(n，k)；然后计算m＝-1时的-1号波束，相移因子P(n，k，-1)＝P^*(n，k，1)；对m＝2时的2号波束，所需的相移因子P(n，k，2)＝P(n，k，1)_η(n，k)，符号_表示矩阵对应元素相乘；重复以上过程，直至完成所有波束的相移因子计算。这一过程可以在现有的数字信号处理系统中方便的实现。

在本步骤中，利用到了分数点的DFT循环移位。基于整数点循环移位的DFT定理套用到分数点移位时会导致错误的结果，如果对频域相移因子进行适当的修正，就可以将整数点循环移位的DFT定理推广到分数点移位的场合，分数点DFT循环移位的(即m可为分数)移相因子为：

$S(k,m)=\left\{\begin{array}{cc}\exp (-j2\mathrm{\πkm}),& k=\mathrm{0,1},...,[N/2];\\ \exp [-j2\mathrm{\π}(N-k)m],& k=[N/2]+1,...,N;\end{array}\right.$

这里符号[·]为取整运算，m为移位的点数，k为频率标号，N为数据段点数。需要注意的是，频域移位方法实际上是对数据段的周期延拓序列进行移位的，由于数据段的两端点取值不同，可能会在移位后的序列中出现端点冲激现象，这是与整数点移位定理不同的地方。为了消除端点冲激，本发明采用对数据段两端加窗平滑的方法。

在步骤106中，直接利用105步骤得到的移相复系数与104步骤得到的频域数据块进行复乘运算，这样对任一波束方向需要进行N_E×N_F次复乘运算。

X_m′(n，k)＝X(n，k)×P_m(n，k)；

这里P_m(n，k)与前面的P(n，k，m)和附图2中的P_n，m(k)含义相同，只是写法上的区别。

在步骤107中，对步骤106得到的结果进行加和得到该波束方向的频域输出，如对第m号波束：

$Y_m\left(k\right)=\underset{n=1}{\overset{N_E}{\mathrm{\Σ}}}{X}_m^{\′}(n,k),$   k＝1，2，...，N_F；m＝1，2，...，N_B；

在步骤108中，判断是否所有的波束形成完毕，若判断为否，则回到步骤105，进行下一波束的相移因子和波束形成计算。若判断为是，则进入步骤109。

在步骤109中，对频域波束输出数据(大小为N_B×N_F，N_B为波束数目)进行逆傅立叶变换(IFFT)变换，拼接相邻数据块，得到连续时域波束输出。至此，一段数据的实时波束形成过程结束。

图4给出了根据分数点移位方法得到的延时0.8个采样周期后的时域波形，可以看到延时后波形与原始波形符合很好。由于可以得到精确的分数点移位序列，所以时延补偿精度与信号采样率没有关系，这是本发明可以采用较低采样率实现精确延时补偿的根本原因。

本发明的波束形成方法着重于频域精确时延补偿的实现，利用分数点DFT循环移位方法可以对宽带信号进行精确延时，这种延时精度与信号的采样率无关，从而可以最大限度降低采样频率，减少冗余数据，大大提高波束形成实现的效率和减少系统复杂度。对于相移操作需要的复相位因子，可以利用前述的对称波束分布特点通过乘幂运算得到，运算量约增加0.5倍。

Claims (6)

1、一种宽带频域数字波束形成方法，其特征在于，包括如下步骤：

1)对接收阵各阵元接收到的信号进行前置放大；

2)在信号处理系统中，对前置放大后的信号进行正交采样得到复包络信号，该正交采样过程包括前置滤波和次级放大，高速模数转换，正交解调，低通滤波以及数据重采样；

3)对多通道的连续数据流进行分块处理，将数据流分解为相互重叠，大小为N_E×N_F的数据块，其中符号N_E表示阵元数，N_F表示进行傅立叶变换的数据长度；

4)对每个数据块沿行方向进行快速傅里叶变换，得到各通道数据的频域形式X(n，k)，其中n为阵元标号，k为频率标号；

5)计算各阵元在不同波束方向上的相移因子；将不同方向波束编号，成像区域的波束方向θ满足 $\sin \mathrm{\θ}=0,\±\frac{1}{Z},\±\frac{2}{Z},...,\±\frac{m}{Z},...;$ 其中m为波束编号，m≤Z，Z为正整数，则对第m号波束，其相移因子P(n，k，m)＝(η(n，k))^m；

$\mathrm{\η}(n,k)=\exp \{-j\frac{2\mathrm{\πd}\left(n\right)}{\mathrm{cZ}}\·(f_c-\frac{k^{\′}}{F_s{N}_F})\};$

其中： $k^{\′}=\left\{\begin{array}{cc}k& k=1,...,N_F/2;\\ k-N_F& k=(N_F/2)+1,...,N_F;\end{array}\right.$

式中，n为阵元标号，k为频率标号，d(n)为第n阵元相对于阵中心的距离，c为介质中波速，f_c为阵元接收信号的中心频率，F_s为正交采样后的对复包络信号的采样率，N_F为傅立叶变换的数据长度。

6)利用步骤5)得到的相移因子与步骤4)得到的频域数据块进行复数乘法运算，得到各波束方向上的完成补偿的频域数据，对第m号波束，补偿后的频域数据

X_m′(n，k)＝X(n，k)×P(n，k，m)；

7)对步骤6)得到的结果进行加和得到该波束方向的频域输出Y(k)，对第m号波束，其频域输出

$Y_m\left(k\right)=\underset{n=1}{\overset{N_E}{\mathrm{\Σ}}}{X}_m^{\′}(n,k),k=\mathrm{1,2},...,N_F;m=\mathrm{1,2},...,N_B;$

其中N_B为波束数目；

8)判断是否所有的波束形成完毕，若判断为否，则回到步骤5)，进行下一波束的相移因子和波束形成计算，若判断为是，则进入步骤9)；

9)对所有频域波束输出数据进行逆快速傅立叶变换，拼接相邻数据块，得到连续时域波束输出。

2、按权利要求1所述的宽带频域数字波束形成方法，其特征在于，所述步骤2)中前置滤波的滤波通带选择1.5-4倍的阵元接收信号带宽。

3、按权利要求1所述的宽带频域数字波束形成方法，其特征在于，所述步骤2)中的模数转换由同步时钟统一控制。

4、按权利要求1所述的宽带频域数字波束形成方法，其特征在于，所述步骤2)中，每一通道的正交解调均由2个乘法器完成，其参考周期序列为sin(-2πf_cn/F_B)和cos(2πf_cn/F_B)，其中f_c为阵元接收信号的中心频率，n为整数，F_B为模数转换采样率。

5、按权利要求1所述的宽带频域数字波束形成方法，其特征在于，所述步骤2)中，数据重采样后的复包络信号采样率为1.2～2倍奈奎斯特频率。

6、按权利要求1所述的宽带频域数字波束形成方法，其特征在于，所述步骤5)中，在存储器中存储中间参数η(n，k)，通过累乘运算实时计算获得相应波束的相位因子P(n，k，m)。

Images