技术领域
本发明属于声纳数字信号处理领域,特别涉及一种实时频域超分辨方位估计方 法及装置。
背景技术
声纳
数字信号处理的一个基本问题就是方位估计,即确定感兴趣的空间信号的 方向,又称目标测向。高分辨方位估计,它是指采用某种
算法,对阵列接受到的目 标
辐射或反射的
时空二维信号进行处理,从而得到目标方位的高分辨估计。超分辨 方位估计技术能超越传统高分辨技术性能,在声纳、雷达及通信等领域有着广泛的 应用。
目前典型的用于声纳的高分辨方位估计算法有:Capon法、窄带二维频域波束 形成法、MUSIC(多重信号分类)法、自适应波束形成法以及ESPRIT(旋转不变量信 号参数估计方法)法。
其中,窄带二维频域波束形成法(可参考文献:“孙长瑜等,二维频域波束形成 方法,应用声学,1995”),是指在窄带目标信号输入的条件下,在时域和空间域对 输入阵列数据进行FFT(快速
傅立叶变换),从而在频域得到波束形成。所述窄带二 维频域波束形成法相对于Capon法等其它几种典型方法,具有计算高效、快速,便 于DSP工程实现,能实现实时处理等优点。所述窄带二维频域波束形成法的缺点是: 只能处理窄带信号,与测向时经常遇到的宽带信号不匹配,实用性差。
ESPRIT法(可参考文献:“Roy R等,ESPRIT-A Subspace Rotation approach to estimation of parameters of cissoids in noise,IEEE Trans.ASSP,1986”),是指一种利用 阵列
传感器分组所得的信号子空间中隐含的旋转不变量实现测量来波方向的方法。 子空间的不变量是通过建立自相关和互相关矩阵来实现的。将原阵列分解为两个子 阵,通过最小二乘(Least Square,LS)或总体最小二乘(Total Least Square,TLS)拟 合来估计子阵列的
相移矩阵,从而得到目标信号方位估计。所述ESPRIT法的
缺陷是: 在计算中涉及元素间互相关值的估计时,由于环境的不确定性,噪声方差值的估计 经常不准确,所以从自相关和互相关矩阵中扣除估得的噪声方差有时会出现致命错 误,从而造成系统的
稳定性差,实用性能较差。
除此以外,近年来,由于现代
信号处理理论的迅速发展,高阶累积量、时-频分 析、小波分析、循环平稳信号分析与处理等理论和方法在
水下目标方位估计中也得 到了广泛的研究,但这些方法要么计算量巨大,系统难以负荷,要么稳健性不高, 难以在实际的应用中得到采用。
总体而言,由于声纳工作环境的恶劣性(噪声模型易失配,海洋传播信道的复 杂,阵形易失配,
信噪比低等),以上算法在现阶段的实际应用中,效果不是很理想, 实际应用中迫切需要一种既能快速实现又具有高
精度的方位估计方法及装置。
发明内容
本发明的目的是克服
现有技术的不足,解决现有技术中计算速度慢、稳健性差、 难以处理宽带目标信号及分辨
力不高的问题,提供一种用于线列阵的快速稳健的高 精度的方位估计方法及装置。本发明可以扩展到平面阵等其它接收阵。
为了实现上述目的,本发明提供的实时频域超分辨方位估计方法包括以下步骤:
1)对线列阵进行划分,得到M1个子阵,所述M1至少为2;
2)对每个子阵的时空二维信号分别进行波束形成,得到每个子阵在扫描方位上 的输出波束;
3)在扫描方位上对各子阵的输出波束进行代数组合得到多子阵合成处理后的波 束。
上述技术方案中,还包括步骤4)在波束域内通过
能量检测器扫描得到目标方位。
上述技术方案中,所述步骤1)中,所述子阵可以是对单个线列阵在逻辑上进行 分割得到的,也可以是物理上区隔的多个单阵。
上述技术方案中,所述步骤1)包括如下子步骤:
11)在探测
角度和探测频带内,对任意角度任意频点产生
频率-
波数网格;
12)产生子阵间延时表,所述子阵间延时表记录子阵间每个频点
位置在每个扫 描角下的延时。
上述技术方案中,所述步骤2)中,所述输出波束可以是时间-方位数据矩阵, 也可以是频率-方位数据矩阵。
上述技术方案中,所述步骤2)包括如下子步骤:
21)用各子阵接收时空二维信号,对各子阵进行波束形成处理,所述波束形成 处理包括步骤22)至25);
22)对子阵的各阵元上的接收数据在时间域上做
快速傅立叶变换,获得各阵元 的频域数据;
23)对子阵获得的频域数据,在空间域上补零,然后在目标频段内的每一个频 率分量上,做空间域快速傅立叶变换,得到一组与各频点和波束相对应的数据;
24)在处理带宽内的每一个频点上,依据所述频率-波数网格和子阵间延时表对 步骤23)得到的数据进行校正,得到频率-方位数据矩阵;
25)对步骤24)获得的频率-方位数据矩阵在频域做反快速傅立叶变换,得到时 间-方位数据矩阵,并存储子阵处理最终数据。
上述技术方案中,所述步骤3)包括如下子步骤:
31)对每一个子阵输出的波束,在探测角度范围内的任意角度上,得到各子阵 输出波束的绝对值和两两子阵间输出波束差的绝对值;
32)对上述波束绝对值进行代数组合,得到新的波束;
上述技术方案中,所述步骤4)包括如下子步骤:
41)求步骤32)得到波束的能量,得到方位-幅度谱,在角度域内幅度最大的位 置即认为是目标信号位置,输出估计的目标方位;
上述技术方案中,所述步骤1)中,对线列阵在逻辑上进行分割,子阵的划分方 式有很多种,有三种常见的子阵划分方式,一种为子阵不重叠的划分方式;一种为 子阵重叠的划分方式;一种为子阵交叉的划分方式。各子阵包含阵元数可以相等, 也可以不等。
上述技术方案中,所述步骤23)中,在做空间域快速傅立叶变换前,首先需要 在空间域对原数据进行补零。
上述技术方案中,所述在空间域对原数据进行补零既可以直接在原数据后直接 补零,也可以在原数据中插值补零。
上述技术方案中,所述步骤31)中,对各子阵输出波束进行代数组合的公式如 下:
B=(B_sumn-B_diffn)1/n
其中n称为控制因子,取值范围为正实数集;Q值的选择区间为[1,M1-1];Bi 是第i个子阵输出波束;wi为正实数。
为了实现本发明的另一目的,本发明提供的实时频域超分辨方位估计装置,包 括:
子阵划分单元,用于对线列阵进行划分,得到M1个子阵,所述M1至少为2;
波束形成单元,用于对每个子阵的时空二维信号分别进行波束形成,得到每个 子阵在扫描方位上的输出波束;
合成处理单元,用于在扫描方位上对各子阵的输出波束进行合成处理得到多子 阵合成波束;以及
目标方位判定单元,用于根据合成处理单元得到的多子阵合成波束获取目标方 位。
上述技术方案中,所述子阵划分单元包括:
频率-波数网格生成单元,用于在探测角度和探测频带内,对任意角度任意频点 产生频率-波数网格;以及
子阵间延时表生成单元,用于产生子阵间延时表,所述子阵间延时表记录子阵 间每个频点位置在每个扫描角下的延时。
上述技术方案中,所述波束形成单元包括:
信号输入单元,用于输入各子阵所接收到的时空二维信号;
时间域快速傅立叶变换单元,用于对子阵的各阵元上的接收数据在时间域上做 快速傅立叶变换,获得各阵元的频域数据;
空间域快速傅立叶变换单元,用于对子阵获得的频域数据,在空间域上补零, 然后在目标频段内的每一个频率分量上,做空间域快速傅立叶变换,得到一组与各 频点和波束相对应的数据;
数据校正单元,用于在处理带宽内的每一个频点上,依据所述频率-波数网格和 子阵间延时表,对所述空间域快速傅立叶变换单元得到的数据进行校正,得到频率- 方位数据矩阵;以及
反快速傅立叶变换单元,用于对所述数据校正单元输出的频率-方位数据矩阵在 频域做反快速傅立叶变换,得到时间-方位数据矩阵,并存储子阵处理最终数据。
上述技术方案中,所述目标方位判定单元包括:
能谱计算单元,用于计算所述合成处理单元输出的波束的能量,得到方位-幅度 谱;以及
判定单元,用于根据所述方位-幅度谱,找出在角度域内幅度最大的位置并判定 该位置为目标方位。
与现有技术相比,本发明具有如下技术效果:
(1)采用频域波束形成,充分利用快速傅立叶变换的计算高效性,计算速度快, 便于DSP工程实现,能实现实时处理;
(2)采用快速傅立叶变换,在频域上运算,将宽带信号分解为多个频点信号, 对每一个频点进行处理,适应于宽带噪声目标测向;
(3)采用子阵处理,在子阵内保证了信号的相关性,同时对各子阵输出波束进 行代数组合,可以得到高
分辨率的波束;
(4)采用本发明后,与现有技术中的频域波束形成算法相比较,可以得到主瓣 更窄,旁瓣强度更弱的波束输出,从而得到更高分辨力的目标方位估计。
附图说明
以下,结合附图来详细说明本发明的
实施例,其中:
图1是现有技术中的常规线列阵示意图;
图2是本发明中线列阵子阵分割方案示意图;
图3是本发明中频率-波数网格的示意图;
图4是本发明实时高分辨方位估计方法的总体
流程图;
图5是本发明实时高分辨方位估计方法的具体流程图;
图6是本发明实施例的声纳
数据处理和显控系统;
图7是在一个仿真目标下本发明实施例与现有技术的分辨效果对比图;
图8是在两个仿真目标下本发明实施例与现有技术的分辨效果对比图;
图9是采用频域波束形成算法对实际数据处理后的效果图;图中白色表示目 标的轨迹;
图10是本发明实施例用于对实际数据处理后的效果图;图中白色表示目标的 轨迹。
具体实施方式
本发明的基本构思为:利用快速傅立叶变换计算的高效性,结合频域波束形成 算法、子阵处理技术和子阵波束代数组合技术,即可以大大降低计算的复杂性,同 时也能获得高精度的宽带目标方位估计。本发明的总体流程如图4所示。
下面结合附图和具体实施方式对本发明做进一步的详细描述:
实施例1
如图5所示,本实施例的用于线列阵的实时高分辨方位估计方法及装置的具体 步骤如下:
步骤501:用单条线列阵接收时空二维信号,从而获得该线列阵的各阵元时域数 据。
步骤502:对该线列阵进行逻辑划分,得到M1个子阵;
本步骤中,子阵的划分方式主要有三种,如图2所示,201为子阵不重叠的划分 方式,这种划分方式计算量小,适合于不是很长的线列阵;当各个子阵中心之间的 距离过大时,不重叠的子阵划分方式会导致输出波束中出现栅瓣;202为子阵重叠的 划分方式,这种方式由于阵元的重叠,可以避免栅瓣的出现,但是增加计算量;203 为子阵交叉的划分方式,这
种子阵划分的方法,保留了整条线阵的分辨率,但是如 果子阵数过多,将不可避免出现栅瓣现象。本实施例中采用201子阵不重叠的划分 方式。
步骤503:对各子阵中的每个阵元,分别取一段时间内的接收数据,在时间域上 做快速傅立叶变换(FFT),得到各阵元的频域数据,从而组成一个频域-空间域数据 矩阵。
本步骤中,需要注意的是:因为取的数据长度越长,即数据矩阵维数越大,则 运算速度会相应减慢。为了保证运算速度能够满足实时处理的要求,数据长度不宜 过大。一般FFT运算点数为2的幂次方,通常取1024点、2048点等。
步骤504:在空间域上对上一步的数据补零,然后对处理带宽内的每一个频点, 在空间域上做快速傅立叶变换,得到一组与各频点和波束相对应的数据,即频率-波数 数据矩阵;
由于线阵总阵元数M一般不大,为了在空间域上做FFT后,频率-波数网格(如 图3所示)中有更“稠密”的分布,使得样本点能接近任意的导向角,在空间域上对 数据补零至Np点,Np远远大于M,一般情况下,Np至少为M的8倍。
空间域补零可以有多种选择。既可以直接在原数据后直接补零,也可以在原数 据中插值补零,只要补零后的数据长度满足要求即可。
步骤505:为下一步的校正,预先产生频率-波数网格。
具体地说,本步骤中,预先产生的频率-波数网格由以下各式定义:
t(k,m)=round(NP*fk*τ(θm)+0.5),fk=k*fs/N, τ(θm)=d*sin(θm)/c,
这里k是频点,t(k,m)表示频率-波数网格中频点k和波数m位置的数据,fs,fk 分别为
采样频率和第k个频点频率,d是阵元间距,θm是扫描角,NP是空间域补零 后点数,N是FFT点数,round(x)表示取最接近x的整数,x(n,:)表示取矩阵x的第 n行,c是
声波在水中的传播速度。
步骤506:为了聚焦各子阵输出波束,预先产生子阵延时表。本实施例做法是选 定一个基准子阵,计算其它子阵与该基准子阵间的延时表。不失一般性,把第一个 子阵作为所述基准子阵,所述子阵延时表记录其它子阵与第一个子阵间,每个频点 位置在每个扫描角下的延时;
具体的说,本步骤
中子阵间延时表的生成由下式决定:
上式中,fs为
采样频率,D是子阵阵间距(即子阵中心之间的间距),i是处理 带宽内的频点位置,j是扫描角。
步骤507:在每一个频点上依据预先产生的频率-波数网格对由步骤503所得的 频域-波数矩阵进行校正。
采用FFT实现频域波束形成所对应的角度是波数域上的值,不是真正的角度, 为了得到对应真实的角度的值,一般需要对其进行插值,插值一般意味着巨大的计 算量,本发明采用预先生成用来校正的频率-波数网格的方法进行校正,该方法可参 考文献“Brian Maranda,Efficient digital beamforming in the frequency domain,1989, J.Acoustical Society of America”。
具体地说,本步骤中,频率-波数网格校正由下式定义:
P2(k,:)=P(k,t(k,:)+NP/2)
这里k是频点,P是空间域补零后FFT得到的数据,P2是经频率-波数网格校正 后数据,NP是空间域补零后点数,x(n,:)表示取矩阵x的第n行。
步骤508:依据子阵延时表,对各子阵输出波束进行聚焦,输出聚焦后的频率- 波数数据矩阵P3;频率-波数数据矩阵就是前文中所述的频率-方位数据矩阵。
具体的说,本步骤中,频率-波数网格校正后数据乘以一个延时量由下式定义:
P2(k,:)=P(k,t(k,:)+NP/2),
P2(k,:)=P2(k,:)*exp(j*pp*delay(k,:))
这里pp是子阵序号,k是频点,P是空间域补零后FFT得到的数据,P2是经频 率-波数网格校正后数据,P3是经频率-波数网格校正后数据乘以一个延时量后的数 据,NP是空间域补零后点数,x(n,:)表示取矩阵x的第n行。
步骤509:对步骤508)获得的频率-波数数据矩阵在频域做反快速傅立叶变换, 获得时间-方位数据矩阵;
步骤510:各子阵输出波束代数组合;所述各子阵输出波束既可以是步骤508输 出聚焦后的频率-波数数据矩阵,也可以是步骤509输出的时间-方位数据矩阵。
本步骤是核心步骤,对各子阵输出波束进行代数组合的公式如下:
B=(B_sumn-B_diffn)1/n
其中n称为控制因子,取值范围为正实数集;Q值的选择区间为[1,M1-1];Bi是第 i个子阵输出波束;wi为正实数,一般可取为Q/M2。
在实际应用中,控制因子n一般取为0.5,1,2,这样便于实现。当线列阵的总 阵元数32时,n取0.5。
步骤511:将步骤510获得的波束共轭转置与其自身相乘,结果即为在搜索方位 上的输出功率;
本步骤中时间-波数数据矩阵的一列表示一个搜索方位,一行表示一个时间点的 信号。
搜索角度的步长越小,输出的波束越多,探测的精度也越高。
步骤512:输出目标方位估计;
本步骤中对各方位上的能量进行比较,取能量最大的方位位置为估计的目标方 位。
对长度很大的数据,可以采取分段处理的方法,对每个分段数据进行上述用于 线列阵的快速稳健的高精度测向,再积累输出各个分段的检测结果。
图7、图8、图9和图10,可以看出,本发明与频域常规波束形成算法相比,有 更窄的主瓣宽度,更高的主副瓣高度比和分辨率更高的性能。子阵数为4,每个子阵 阵元数为8,子阵间距为1.2m,控制因子等于0.5,Q值为3,公式510中的ωi均等 于1/2。
需要说明的是,本发明的子阵即可以是对单个线列阵在逻辑上进行分割得到的, 也可以是物理上区隔的多个单阵。采用物理上区隔的多个单阵时,在划分子阵后的 处理步骤与上述步骤503至512完全一致。
下面对本实施例的
硬件环境的一个示例进行描述。
本实施例可以直接使用现有技术中的常规线列阵,如图1所示,用于接收空间 信号的线列阵声纳102的线阵是由32个市场所售常规的接收信号中心频率为 5000KHz的无
指向性水听器101组成,两两水听器间的间隔是半个
波长,从而线阵 声振段的总长度是4.65米。线列阵102安装在潜艇或无人潜航器上。
图6是本发明实施例的声纳数据处理和显控系统。拖曳线列阵102和前置
电路 601从每个水听器接收到
模拟信号。前置电路601包括前放、滤波和其它常规的电路。 每个通道的模拟信号输入到A/D转换器602得到数字信号。从A/D转换器602出来 的是多通道的数字数据流,每个通道的数据流对应于一个水听器接收到的模拟信号。 将这些数据流输入到微型处理器603。经微型处理器603处理后的输出信息可以存储 在数据
存储器件605,比如磁盘存储设备中,或直接输出到显示设备606上显示。
微型处理器603首先将接收到的数据流存储到动态存取区604,在输入满足处理 要求数量的数据流后就开始处理。图5是算法具体流程示意图,它包括时间域FFT 503,空间域补零后做FFT504,频率-波数网格505,产生子阵间延时表506,频率- 波数校正507,子阵聚焦508,频域IFFT 509,各子阵输出波束代数组合510,能量 检测511及目标方位估计512。这些程序储存在动态存取区604中。
因为本发明中的A/D转换器602输出的是多通道的数据流,因此可以采用多片 的微型处理器来并行处理。能够实现图6功能的其它一些硬件设备,比如专用硬件、 基于应用的集成电路(ASIC)、DSP、ARM等都可以用来代替微型处理器603。
最后所应说明的是,以上实施例仅用以说明本发明的技术方案而非限制。尽管 参照实施例对本发明进行了详细说明,本领域的普通技术人员应当理解,对本发明 的技术方案进行
修改或者等同替换,都不脱离本发明技术方案的精神和范围,其均 应涵盖在本发明的
权利要求范围当中。