一种空时二维滤波的快速宽带波束形成方法与流程

本发明属于水下声传感器阵列信号处理领域，具体涉及一种空时二维滤波的快速宽带波束形成方法。

背景技术：

在水下被动定位中，目标声源往往为宽带声信号。快速、准确的宽带波束形成技术对提升水下武器装备的性能具有重要意义。

常规的宽带数字波束形成建立在量化时延的基础上，因而输出波束的扫描范围、角度步进值等参量均受到采样率的限制。为了保证输出波束的质量，往往需要较高的采样率。当接收阵为标准半波长均匀线阵时，描述了宽带时域波束形成的过程。其中，n表示线阵中一共存在n个阵元，xi(n)表示第i个通道的接收信号时间序列。τs为时延补偿，该补偿必须为采样周期的整数倍。即τs＝kts,k＝0,±1,±2,…，其中ts为采样周期，与采样率fs的关系为

由于k只能取离散值，所以τs只能为离散值。若希望使用更小数值的τs进行精确度更高的扫描，则只能使用更小的ts值，即需要提高系统的采样率。根据采样定理，为了能够不失真地恢复信号，采样率应当为信号最高频率的2倍，即萘奎斯特(niquest)采样定理。在数字多波束中，为了保证增益不明显下降，通常要求采样率为信号最高频率的3-5倍。而这种要求只能保证增益下降不多，若希望波束指向性图失真很小，则一般要求采样率为信号最高频率的10倍以上，由此导致高采样率。高采样率对系统整体性能提出了更高的要求。在数据采集阶段，高采样率要求使用高转换速率的ad芯片，造成了硬件实现上的困难，提高了系统成本。在数据存储和传输阶段，高采样率导致了大的数据吞吐量，对数据存储设备、数据传输过程都提出了更高的要求。在数据处理阶段，更大的数据量要求处理器具有更强的计算能力。因此，需要寻求一种方法，既能够保证输出波束的质量，又能够保持较低的采样率。

采样定理表明，只要按照萘奎斯特采样定律对模拟信号进行采样，就可以无失真地将其恢复。所以，可以按照较低的采样率进行采样，在处理阶段提高采样率，这样便可以在使用低样率的同时保证输出波束的质量。通常的方法是，对每个通道采集的信号进行升采样操作，然后在高采样率下进行时延波束形成。对于存在n个通道的声呐接收阵，该方法需要进行n次插值操作。考虑到插值是一种计算量较大的数值计算方法，所以这种方法虽然降低了对硬件设备的要求，但具有较大的计算量。利用时延波束形成和升采样均为线性过程这一特点，还可以先对接收信号进行补零，然后进行时延波束形成，最后再通过插值的方法获得高采样率下的波束。对于存在n个通道的接收阵，这种方法只需要进行1次插值操作，所以降低了计算量。但由于补零过程的存在，使得在波束形成阶段的计算量成倍地增加。当希望进行m倍升采样时，计算量则提升m倍。这限制了宽带时延波束形成的快速实现。

技术实现要素：

针对升采样导致计算量成倍地增加的问题，本发明提出了一种基于空时二维滤波的快速宽带时延波束形成方法。一方面，该方法能够得到任意原采样率整数倍数的采样率下的输出波束，从而能够保证输出波束的精度。另一方面，该方法在保证输出波束质量的前提下几乎没有增大计算量。实际上，对于传统方法而言，升采样m倍时计算量也提升了m倍。而本文所述方法为传统方法的即本发明所述方法在不增加计算量的前提下完成了升采样m倍的任务，计算量与升采样倍数无关且保持恒定。

下面对本发明方法的原理进行阐述。

一种空时二维滤波的快速宽带波束形成方法，包括以下步骤：

(1)以符合niquest定理的采样率进行采样，每个通道采集的数据按时间顺序排列成列向量，组成一个t×n维的数据矩阵，列的编号依次为j＝0,1,…,n-1，行的编号依次为i＝0,1,…,t-1；

(2)确定采样倍数m以及相邻通道时延补偿点数t0；

(3)计算坐标矩阵主值区间p⁽⁰⁾(t0)，规定坐标矩阵的列编号从左至右依次为j＝0,1,…,n-1，行编号从上至下依次为i＝0,1,…,m-1，则p⁽⁰⁾(t0)的第i行第j列元素为

(4)去除p⁽⁰⁾(t0)中非整数点，只保留整数点；p⁽⁰⁾(t0)的每列有且只有1个整数点元素，设定循环因子k＝0；

(5)依次取出p^(k)(t0)中第0～m-1行的整数点，根据其数值从数据矩阵中取出对应的数据点进行累加，分别得到波束的第km,km+1,…,km+m-1个点；

(6)对p^(k)(t0)中所有整数点数值加1，得到下一个周期的坐标矩阵p^(k+1)(t0)；更新循环因子：k←k+1；

(7)判断km≥t是否成立，若成立，则说明已经遍历数据矩阵中所有数值点，执行步骤(8)，否则执行步骤(5)；

(8)计算滤波器通带范围；

(9)根据步骤(8)中计算的参数设计滤波器对进行滤波，得到输出波束y(i)。

所述计算坐标矩阵主值区间p⁽⁰⁾(t0)，规定坐标矩阵的列编号从左至右依次为j＝0,1,…,n-1，行编号从上至下依次为i＝0,1,…,m-1，则p⁽⁰⁾(t0)的第i行第j列元素为包括：

由步骤(2)确定的采样倍数m以及相邻通道时延补偿点数t0得输出波束表达式为：

其中，表示第j个通道中序号为的采样点参与了累加，得到了输出波束的第i个数据点，

显然，必为整数，由niquest采样定理可知，只要采样率大于信号最高频率的2倍，则通过插值的方式恢复信号而无信息丢失；阵列信号中，目标信号最高频率对接收数据进行m倍降采样处理，利用补零和时延求和得到波束然后再通过低通滤波无失真地恢复输出波束，具体过程表示为：

其中，为降采样m倍的信号，为第j个通道参与波束形成的数据序号；

与pij(t0)关系为：

其中，非整数的pij(t0)表示序号为的数据点不参与降采样后的累加过程，即在降采样后，这部分数据点不参与运算；

定义p(t0)为：

pij(t0)＝{p(t0)}ij,(i＝0,1,…,t-1,j＝0,1,…,n-1)

p(t0)为t×n的矩阵，称为坐标矩阵，将坐标矩阵的行从上至下依次称为第0行，第1行，…，第t-1行，同时将坐标矩阵的列从左至右依次称为第0列，第1列，…，第n-1列，于是，p(t0)的第i行，第j列元素为pij(t0)。

所述去除p⁽⁰⁾(t0)中非整数点，只保留整数点；p⁽⁰⁾(t0)的每列有且只有1个整数点元素，设定循环因子k＝0，包括：

由于坐标矩阵表示降采样后数据点的序号，且只有当坐标矩阵为整数的元素表示的数据才能参与波束形成运算，所以只需要找出坐标矩阵中的整数元素，然后找到该元素表示的数据点，对其进行累加即完成了降采样后的波束形成；

下面考虑坐标矩阵的第0行，其元素依次为：

p00(t0)必为整数，且值为0；设下一个为整数的元素为p0j(t0),则

其中，表示整数组成的集合，k表示某个整数；若s0为t0和m的最小公倍数，则均为整数，此时，p0j(t0)为整数；若j'＝kj0,k＝0,1,2,…，则p0j'(t0)亦为整数；所以，坐标矩阵的第0行中存在整数的列编号呈周期性变化，周期为

所述依次取出p^(k)(t0)中第0～m-1行的整数点，根据其数值从数据矩阵中取出对应的数据点进行累加，分别得到波束的第km,km+1,…,km+m-1个点，包括：

p(t0)的第0行第j列出现了整数坐标，则考察p(t0)的第(j+c)列第i行的整数坐标的分布情况，由周期性可知，只需考察第j列到第2j列之间的范围即可，故j＜j+c＜2j，pi(j+c)(t0)为：

由于若希望则令i＝km-ct0，又由于0≤i≤m-1，得到：

若p(t0)的第0行第j列出现了一个整数坐标点，则p(t0)的第(j+c)列第i＝km-ct0行也将出现一个整数坐标点；

若p(t0)的第(j+c)列第i行出现了一个整数点坐标，下面考察第(j+c)列中整数点坐标出现的情况，设以第(j+c)列第i行为起点，向上或者向下移动m行出现下一个整数点，则该整数点的坐标为第(i+m)行，第(j+c)列，于是此处数值为：

由于为整数，故只需要求为整数；所以：

m＝km,k＝0,±1,±2,…

m的取值说明，在第(j+c)列中整数点出现的位置是周期的，周期为m，则：

定义坐标矩阵p(t0)的第km～km+m-1行为坐标矩阵的子矩阵，记为p^(k)(t0)，称k＝0时的子矩阵p⁽⁰⁾(t0)为坐标矩阵的主值区间。

所述对p^(k)(t0)中所有整数点数值加1，得到下一个周期的坐标矩阵p^(k+1)(t0)；更新循环因子：k←k+1，包括：

主值区间内有且只有n个整数点，可求得p(t0)中相邻整数点的数值为：

即同一列中相邻整数点之间的数值总是相差1；

在坐标矩阵中，每行出现整数点是周期的，周期为每列出现整数点的位置是周期的，周期为m，主值区间内每列有且必有1个整数点；由于坐标矩阵是周期的，并且相邻整数点数值相差1，故只需要计算主值区间内的数值，即能得到坐标矩阵全部区域的值。

所述计算滤波器通带范围，包括：

滤波器的下截止频率为wfl，上截止频率为wfh，则wfl和wfh应当满足以下约束：对于实信号和复信号wfl应当满足：对于实信号的上截止频率应当满足：对于复解析信号wfh应当满足：

其中，w0l和w0h是采样率为fs时信号通带的上截止频率和下截止频率。

本发明的有益效果在于：

本发明提出了坐标矩阵计算方法，并利用坐标矩阵直接对有效数据进行计算，略去了大部分的冗余计算量，从而达到快速计算的目的。

附图说明

图1是本发明的流程图；

具体实施方式

下面结合附图对本发明做进一步描述。

针对升采样导致计算量成倍地增加的问题，本发明提出了一种基于空时二维滤波的快速宽带时延波束形成方法。一方面，该方法能够得到任意原采样率整数倍数的采样率下的输出波束，从而能够保证输出波束的精度。另一方面，该方法在保证输出波束质量的前提下几乎没有增大计算量。实际上，对于传统方法而言，升采样m倍时计算量也提升了m倍。而本文所述方法为传统方法的即本发明所述方法在不增加计算量的前提下完成了升采样m倍的任务，计算量与升采样倍数无关且保持恒定。

下面对本发明方法的原理进行阐述。

假设接收系统中共存在n个传感器(又称n个通道)，依次编号为0,1,…,n-1。第j个传感器的输出为i＝0,1,…,t-1，j＝0,1,…,n-1。其中，i表示接收数据的序号，每个通道中共采集了t个数据点。时延波束形成通过在相邻通道之间插入时延补偿的方式完成波束的扫描。在数字系统中，插入的时延只能为1个采样周期的整数倍。设系统采样率为mfs，其中m为整数。设相邻通道之间需要插入个采样点，则输出波束表达式为

其中，表示第j个通道中序号为的采样点参与了累加，得到了输出波束的第i个数据点。的表达式为

显然，必为整数。由萘奎斯特采样定理可知，只要采样率大于信号最高频率的2倍，则可以通过插值的方式恢复信号而无信息丢失。若阵列信号中，目标信号最高频率则可以对接收数据进行m倍降采样处理，利用补零和时延求和得到波束然后再通过低通滤波无失真地恢复输出波束，具体过程表示为

此处，为

其中，为降采样m倍的信号，为第j个通道参与波束形成的数据序号。与pij(t0)关系为

显然，pij(t0)可能不为整数。非整数的pij(t0)表示序号为的数据点不参与降采样后的累加过程，即在降采样后，这部分数据点可以不参与运算。

定义p(t0)为

pij(t0)＝{p(t0)}ij,(i＝0,1,…,t-1,j＝0,1,…,n-1)

则p(t0)为t×n的矩阵，称为坐标矩阵。为描述方便，将坐标矩阵的行从上至下依次称为第0行，第1行，…，第t-1行，同时将坐标矩阵的列从左至右依次称为第0列，第1列，…第n-1列。于是，p(t0)的第i行，第j列元素为pij(t0)。

由于坐标矩阵表示降采样后数据点的序号，且只有当坐标矩阵为整数的元素表示的数据才能参与波束形成运算，所以只需要找出坐标矩阵中的整数元素，然后找到该元素表示的数据点，对其进行累加即可完成降采样后的波束形成。由于坐标矩阵中只有部分元素为整数，故只有部分数据点参与了运算，从而节约了计算量，提高了计算速度。

下面考虑坐标矩阵的第0行。其元素依次为

显然p00(t0)必为整数，值为0。设下一个为整数的元素为p0j(t0),则

其中，表示整数组成的集合，k表示某个整数。若s0为t0和m的最小公倍数，则均为整数。此时，p0j(t0)为整数。若j'＝kj0,k＝0,1,2,…，则p0j'(t0)亦为整数。所以，坐标矩阵的第0行中存在整数的列编号呈周期性变化，周期为

下面利用反证法证明在第0列和第列之间不存在另一个整数坐标。设p0j'(t0)为整数，且0＜j'＜j0，则

于是有

s＝j't0＝k'm

所以，s为t0和m的公倍数，于是应当有

这与假设矛盾。故p(t0)矩阵的第0行的整数坐标周期性地出现，周期为

设p(t0)的第0行第j列出现了整数坐标，下面考察p(t0)的第(j+c)列第i行的整数坐标的分布情况。由周期性可知，只需考察第j列到第2j列之间的范围即可，故j＜j+c＜2j。pi(j+c)(t0)为

由于若希望则当令：

于是有

i＝km-ct0

由于0≤i≤m-1，所以得到

这表示，若p(t0)的第0行第j列出现了一个整数坐标点，则p(t0)的第(j+c)列第i＝km-ct0行也将出现一个整数坐标点。

下面证明，若p0j(t0)为整数，则对于j＜j+c＜2j，0≤i≤m-1，一定存在pi(j+c)(t0)为整数。

证明：

由于则只需证设ct0＝km+r，其中0≤r≤m-1。取i＝m-r，则有

ct0+i＝(k+1)m

证毕。

若p(t0)的第(j+c)列第i行出现了一个整数点坐标，下面考察第(j+c)列中整数点坐标出现的情况。设以第(j+c)列第i行为起点，向上或者向下移动m行出现下一个整数点，则该整数点的坐标为第(i+m)行，第(j+c)列。于是此处数值为

由于为整数，故只需要求为整数。所以，可以得到

m＝km,k＝0,±1,±2,…

m的取值说明，在第(j+c)列中整数点出现的位置是周期的，周期为m。结合上一结论可知，在0～m-1行中，第(j+c)列有且只有1个整数点。为了叙述方便，作出如下定义：

定义坐标矩阵p(t0)的第km～km+m-1行为坐标矩阵的子矩阵，记为p^(k)(t0)。称k＝0时的子矩阵p⁽⁰⁾(t0)为坐标矩阵的主值区间。

于是，主值区间内有且只有n个整数点。可以求得p(t0)中相邻整数点的数值为

即，同一列中相邻整数点之间的数值总是相差1。

综上所述，在坐标矩阵中，每行出现整数点是周期的，周期为每列出现整数点的位置是周期的，周期为m。主值区间内每列有且必有1个整数点。

利用此结论，可以计算出利用坐标矩阵进行波束形成的计算量。由于坐标矩阵是周期的，并且相邻整数点数值相差1，故只需要计算主值区间内的数值，即可得到坐标矩阵全部区域的值。

对于经典方法需要在采样率为mfs情况下对所有数据点进行求和，所以主值区间内共需要计算m×n次加法。若按照本文方法，只需要在采样率为fs时对坐标矩阵中的整数点进行求和，则共需要计算n次加法。于是，利用坐标矩阵，可以将时延波束形成的计算量降低至原方法的并且可以看到，利用坐标矩阵进行波束形成的方法，其计算量与m无关。这个特点允许将采样率提升至任意倍数而不增加计算量。

所以，只要确定了m、t0与n的值，便可以计算得到完整的坐标矩阵中的整数点。在实现过程中，考虑到数据存储的成本，可以利用坐标矩阵的周期性与相邻整数点数值相差1的特点，只需要存储主值区间内的数值，便可推知坐标矩阵任意位置的数值。据此，便可以设计基于空时二维联合滤波的快速时延波束形成算法。

下面考察从得到y(i)的方法。本文所述方法虽然没有对数据补零的操作，但本文方法实际上利用坐标矩阵进行波束形成相当于同时完成了对数据补零并对其累加的操作。所以，也需要设计一个滤波器对波束进行恢复。下面将通过计算得出滤波器应当满足的设计指标。

对于离散信号，当采样率为时fs，设接收信号为x(n)，于是在z域有

频谱为

其中，j为复数单位，满足关系为j²＝-1。

利用坐标矩阵进行时延波束形成实际上相当于在x(n)的相邻采样点中插入(m-1)个零，得到序列可以表示为

于是，序列的z变换为

所以，的频谱为

结合上述公式，可以得到补零前后信号频谱的关系为

观察到式中存在关系：

由于x(w)为周期函数，最小正周期为2π，则

也是周期的，最小正周期为上式表明了补零前后信号的频谱变化关系以及补零后信号频谱的周期性质。利用这两个等式可以实现对补零后频谱的清晰刻画，同时能够为滤波器设计提供指导。

若x(n)存在于频带w0l～w0h之间，则补零后频带变为

同时，由于具有周期性，最小正周期为则中也将出现周期性延拓，频谱范围表达式为

对于实信号，由于频谱存在对称性，所以x(n)在频率-w0h～(-w0l)之间也存在能量分布，且与w0l～w0h之间关于w＝0对称。此时，中也将出现另一种周期性延拓，频谱范围表达式为

的频谱分布范围为上式的叠加。不妨称k＝0时表示的信号为基频信号，|k|＞0时表示的信号称为谐波。可以看到，补零操作使信号出现了谐波，且谐波的频率取决于参数m。在宽带时延波束形成中，基频信号为希望得到的信号，而谐波信号则为干扰，应当避免。

本发明中，采用滤波器将谐波滤除。下面设计滤波器的指标参数。

由于目标信号存在于k＝0表示的区间内，而谐波信号存在于|k|＞0表示的区间内。所以，可以设计一个低通或者带通类型的滤波器，将|k|＞0的部分都滤除。设滤波器的通带范围是wfl～wfh,不失一般性，可以限制wfl≥0。对于wfl，为了使信号不失真，则应当要求

下面考察wfh的选择。下表列出了信号类型与通带的关系。

从表中可以看出，当信号为实信号时，wfh应当满足的限制以及滤波器带宽为

当信号为复解析信号时，wfh应当满足的限制以及滤波器带宽为

显然，复解析信号可以允许更大的过渡带。更宽的过渡带意味着对滤波器要求更低，允许使用更低阶次的滤波器，从而尽可能降低滤波器引入的信号畸变。由于希尔伯特变换能够将实信号变换为复解析信号，则可以借用希尔伯特变换得到复解析信号，从而降低对滤波器的要求。

同时，从式(35)和(36)中也可以看到，对于同一个信号，升采样倍数不能太高。因为随着m的增大，δwr和δwi将逐渐变小。这意味着，随着m的增大，也需要逐渐提高滤波器阶数以获得更窄的过渡带才能够得到基频信号。而随着滤波器阶数的增加，滤波器引起的信号畸变也变得更加严重。所以，升采样倍数m不能设置得任意大，否则将因为滤波器过渡带过窄导致信号严重畸变，或者谐波进入滤波器通带，导致波束形成失败。

一种空时二维滤波的快速宽带波束形成方法，提出了步骤1～10，使得系统可以任意提升采样率而几乎不增加计算量。

滤波器设计指标须满足和其中，m为升采样倍数。w0l,w0h分别为x(w)的下截止角频率和上截止角频率，wfl,wfh分别为多项滤波器的下截止角频率和上截止角频率。

提出了坐标矩阵计算方法，并利用坐标矩阵直接对有效数据进行计算，略去了大部分的冗余计算量，从而达到快速计算的目的。

步骤1：以符合萘奎斯特(niquest)定理的采样率进行采样，每个通道采集的数据按时间顺序排列成列向量，组成一个t×n维的数据矩阵，列的编号依次为j＝0,1,…,n-1，行的编号依次为i＝0,1,…,t-1。

步骤2：确定升采样倍数m以及相邻通道时延补偿点数t0。

步骤3：计算坐标矩阵主值区间p⁽⁰⁾(t0)。规定坐标矩阵的列编号从左至右依次为j＝0,1,…,n-1，行编号从上至下依次为i＝0,1,…,m-1。则p⁽⁰⁾(t0)的第i行第j列元素为

步骤4：去除p⁽⁰⁾(t0)中非整数点，只保留整数点。p⁽⁰⁾(t0)的每列有且只有1个整数点元素。设定循环因子k＝0。

步骤5：依次取出p^(k)(t0)中第0～m-1行的整数点，根据其数值从数据矩阵中取出对应的数据点进行累加，分别得到波束的第km,km+1,…,km+m-1个点。由于p^(k)(t0)中共有m×n个元素，但只存在n个元素为整数，故利用p^(k)(t0)将加法次数从m×n次降至n次，计算量为传统方法的倍(对于传统方法而言，若将采样率提升m倍，则计算量亦提升m倍)。于是，本方法的在时延波束形成阶段的计算量与未升采样时的计算量相同，与升采样倍数无关。故本方法在保持计算量不变的情况下可以任意提升采样率。

步骤6：对p^(k)(t0)中所有整数点数值加1，得到下一个周期的坐标矩阵p^(k+1)(t0)。更新循环因子：k←k+1。

步骤7：判断km≥t是否成立。若成立，则说明已经遍历数据矩阵中所有数值点，执行步骤8，否则执行步骤5。

步骤8：计算滤波器通带范围。设滤波器的下截止频率为wfl，上截止频率为wfh。则wfl和wfh应当满足以下约束：对于实信号和复信号wfl应当满足：对于实信号的上截止频率应当满足对于复解析信号wfh应当满足，其中，w0l和w0h是采样率为fs时信号通带的上截止频率和下截止频率。

步骤9：根据步骤8中计算的参数设计滤波器对进行滤波，得到输出波束y(i)。

步骤10：结束。