CN111175727B - 一种基于条件波数谱密度的宽带信号方位估计的方法 - Google Patents

Abstract

本发明提出的一种基于条件波数谱密度的宽带信号方位估计的方法，包括：步骤1)用空间线阵接收远距离目标的宽带信号，将所述宽带信号在时域和空域做快速傅利叶变换，生成时空二维功率谱密度函数，再将时空二维功率谱密度与阵列接收响应的频谱之比作为条件波数谱密度，对扩展的条件波数谱密度谱图净化，得到条件波数谱的净化谱；步骤2)在条件波数谱的净化谱上选取估计的扫描方位角，计算该估计的扫描方位角对应的条带内峰值能量位置集，采用最小二乘拟合对估计的扫描方位角进行修正，得到修正方位角，将修正方位角与对应的峰值能量作为输出，从而得到方位谱图，方位谱图上的峰值位置对应于所述宽带信号的方位。本发明的方位的估计精度高。

Description

一种基于条件波数谱密度的宽带信号方位估计的方法

技术领域

本发明属于声纳数字信号处理领域，特别涉及一种基于条件波数谱密度的宽带信号方位估计的方法。

背景技术

被动声呐是根据辐射噪声来检测目标并判别目标类型的，宽带连续谱作为舰船的声纹特征之一，从频带内携带能量的角度描述了舰船辐射噪声的特性，对舰船辐射噪声的检测具有不可或缺的作用。基于传感器阵列的方位估计(Direction of Arrival,DOA)是其中重要的一环，面临着阵列孔径有限、空间分辨能力差、对扰动误差的鲁棒性差等问题。

目前关于窄带信号空间测向的算法已趋于成熟，而对于携带更多有效信息、具有更高应用价值的宽带信号，传统的窄带信号算法已不完全适用。宽带DOA估计算法包括常规波束形成(CBF)算法、宽带MVDR波束形成算法和子空间类算法等，这些方法均是水声阵列信号处理中最常用的方法。CBF算法是基于波束扫描的测向方法，具有结构简单计算量小的特点，然而受到“瑞利准则”的限制，空间方位分辨能力差；参考文献[1](Somasundaram SD.Wideband robust capon beamforming for passive sonar.IEEE Journal of OceanicEngineering,2013)提出的宽带最小方差无畸变(MVDR)波束形成方法，将处理频带划分为若干子带，对每个子带进行MVDR窄带波束形成，然后对所有子带波束输出功率求和，得到目标的方位估计，具有一定的高分辨能力，但该算法稳健性差，易因扰动误差(阵元误差、通道误差等误差)和数据协方差矩阵估计误差导致阵增益下降，方位估计性能下降；宽带信号的子空间方位估计算法分为非相干信号子空间算法(Incoherent Signal Subspace,ISS)和相干信号子空间算法(Coherent Signal Subspace，CSS)。参考文献[2](Wang H,KavehM.Coherent signal-subspace processing for the detection and estimation ofangles of arrival of multiple wide-band sources.IEEE Transactions onAcoustics,Speech and Signal Processing,1985)，这类算法通过信号空间和噪声空间的分解实现目标的超分辨估计，但对信噪比和快拍数要求较高，且需要对信源数目和目标方位进行预估。

频率-波数功率谱分析(f-k分析)是地震波数据分析的一种处理方法。参考文献[3](J.Capon.High-Resolution Frequency-Wavenumber SpectrumAnalysis.Proceedings of the IEEE.1969)在分析地震波信号过程中，通过把空间采样与时域平稳随机过程类比，采用时空二维联合功率谱密度对阵列传感器输出的空间平稳场进行分析和描述，提供了空间叠加随机场的幅度均方根值。参考文献[4](Kvaerna T andDoornbos D J.An integrated approach to slowness analysis with arrays andthree-component array data.In Semiannual Technical Summary,1985)提出了地震波的f-k分析方法，利用频率域内的聚束，在慢度空间计算信号的功率谱，得到不同慢度和不同方向波的能量分布，其最大谱值对应信号的慢度矢量，根据慢度矢量可以计算信号的方位角。目前，f-k分析方法已经广泛应用于地震波慢度测定、表面波的损伤检测和表面波色散问题研究。而在声呐领域，参考文献[5](Beerens P,Ijsselmuide S V,Volwerk C,etal.Flow noise analysis of towed sonar arrays.Udt Europe.1999)利用f-k分析了拖曳阵设计中流噪声的成因，取得了良好的效果；而利用信号在f-k空间的能量分布特征进行DOA估计的研究较少。

现有的宽带高分辨方位谱估计算法计算速度较慢，可行性不高且实际应用中需要预先对信源数目和目标方位进行预估。

实际应用中，需要一种能够针对宽带目标信号快速、高分辨、对扰动误差鲁棒性高的方位谱估计实时处理方法。

发明内容

本发明目的在于，为克服现有的高分辨方位谱估计算法计算速度较慢，可行性不高，且实际应用中需要预先对信源数目和目标方位进行预估的缺点，本发明提出一种基于条件波数谱密度的宽带信号方位估计(CWSD-based)方法。包括：

步骤1)用空间线阵接收远距离目标的宽带信号，将所述宽带信号在时域和空域做快速傅利叶变换，生成时空二维功率谱密度函数，再将时空二维功率谱密度与阵列接收响应的频谱之比作为条件波数谱密度，对扩展的条件波数谱密度谱图净化，得到条件波数谱的净化谱；

步骤2)在条件波数谱的净化谱上选取估计的扫描方位角，计算该估计的扫描方位角对应的条带内峰值能量位置集，采用最小二乘拟合对估计的扫描方位角进行修正，得到修正方位角，将修正方位角与对应的峰值能量作为输出，从而得到方位谱图，方位谱图上的峰值位置对应于所述宽带信号的方位。

作为上述方法的一种改进，所述步骤1)具体包括：

步骤1-1)用空间线阵接收宽带信号，得到M个阵元的时域信号ψ[d_m,t]，在时域将数据分为N个快拍，记为快拍1,2,…,n,…,N，每个快拍ψ_n[d_m,t]长度为T，其中，m为阵元序号，取值为自然数，t为采样时间点，d_m＝m*d为第m个阵元的位置，d为阵元间距；

步骤1-2)对第n个快拍数据在时间域上做快速傅利叶变换，得到目标辐射信号时域频谱ψ[d_m,f]；

其中，T为采样点数，取值为自然数，f表示频率，ψ_n[d_m,t]为第m个阵元接收的时域信号；

步骤1-3)确定目标辐射信号ψ[d_m,f]的频带范围[f_min,f_max]，f_min为频率最小值，f_max为频带最大值，对频带范围内的每一个频点的空域数据，在空域做快速傅利叶变换，快速傅利叶变换点数为Q，得到时空二维频谱Φ[k,f]：

当Q＞M时，Q和M取值为自然数，采用空间域上补零，Φ[k,f]在空间域作快速傅利叶变换，生成时空二维功率谱密度函数S[k,f]：

其中，E[·]表示数学期望，k为波数，取值为自然数；

步骤1-4)参考阵元功率谱密度S(f)为：

其中，上标*表示复共轭；

将零频分量移至谱中心，得到条件波数谱密度函数S[k|f]：

S[k|f]＝S[k,f]/S(f) (7)；

步骤1-5)沿波数k方向对条件波数谱密度函数S[k|f]进行周期性延拓，得到扩展的条件波数谱密度函数；

步骤1-6)选取扩展的条件波数谱密度函数中局部振幅极大值点为宽带信号峰值能量，将其余点置零，计算所述条件波数谱的净化谱；在条件波数谱的净化谱上，宽带信号峰值能量位置呈与入射方位角相关的线性分布。

作为上述方法的一种改进，所述波数k为：

其中f表示频率，d是阵元间距，c为声波传播速度，θ是扫描方位角，取值范围0～180度。

作为上述方法的一种改进，所述步骤1-3)的采用空间域上补零，具体包括采用直接在原数据后补零，或在原数据中插值补零。

作为上述方法的一种改进，所述步骤1-3)中，单信号的时空二维功率谱密度为：

作为上述方法的一种改进，所述步骤1-3)的傅利叶变换点数Q为阵元数M的整数倍。

作为上述方法的一种改进，步骤1-4)中单信号的条件波数谱密度函数为：

作为上述方法的一种改进，所述步骤1-5)的沿波数k方向对条件波数谱密度函数S[k|f]进行周期性延拓为沿波束k方向向左向右同时复制多条条件波数谱密度函数S[k|f]。

作为上述方法的一种改进，所述步骤2)具体包括：

步骤2-1)在方位角范围(0～180度)内，设定初始扫描方位角θ＝θ₀，并计算扫描斜率μ：

步骤2-2)以过(f,k)＝(0,0)、斜率为μ的直线l₀为中心，宽度2Q/M+1的条带L作为参数域方位，计算局部峰值能量位置集P＝{(k_i,f_i)}；k_i为第i个位置的波数，f_i为第i个位置的频率；

步骤2-3)对条带L内条件波数谱密度函数S[k|f]的局部峰值能量位置集P＝{(k_i,f_i)}做最小二乘直线拟合，得到修正斜率μ′，拟合过程为：

根据修正斜率μ′，得到修正波达方位cosθ_s：

θ_s为修正方位角；

步骤2-4)改变扫描方位角θ作为估计的扫描方位角，θ为按照指定的顺序选取的方位角；

θ＝θ₀+lΔθ (12)

θ₀为初始扫描方位角，l为第l次扫描，取值为自然数，Δθ为方位角扫描间隔，取值为0～180度。

判断所述估计的扫描方位角θ是否在有效声锥范围内，如果“是”，返回步骤2-1)继续执行；

如果“否”，则对条件波数谱密度函数在有效声锥范围内的扫描斜率完成遍历；

步骤2-5)以每一个修正方位角为横坐标，修正方位角对应幅度值为纵坐标绘制方位谱曲线，曲线峰值处为能量输出的最大值，峰值处对应修正方位角为精确的宽带信号方位角，从而实现对目标方位的精确定位。

作为上述方法的一种改进，所述有效声锥范围为条件波数谱密度S[k|f]上斜率μ的一个锥形区域内，即斜率μ满足：

本发明的优势在于：

1、本发明的基于条件波数谱密度的宽带信号方位估计的方法利用条件波数谱密度函数将宽带目标信号转换到f-k空间，宽带目标信号能量在该空间呈与入射角相关的线性分布，通过直线检测原理，实现方位谱估计，对宽带信号方位的估计精度和分辨能力明显优于常规波束形成和宽带最小方差无畸变算法，具有较高的主副瓣比，对微弱目标的检测能力更优异；

2、本发明的基于条件波数谱密度的宽带信号方位估计的方法对阵元位置误差具有很好的宽容性；

3、本发明的基于条件波数谱密度的宽带信号方位估计的方法方位估计均方误差小，在宽带处理频带最高频为f时，理论分辨率可达2c/Mfd；

4、本发明的基于条件波数谱密度的宽带信号方位估计的方法采用两重快速傅立叶变换，对信号频带范围内的条件波数谱密度进行直线扫描处理，速度快，实时性好，算法稳健性高；

5、海试数据证明本发明的方法有效。

附图说明

图1为本发明的基于条件波数谱密度的宽带信号方位估计的方法所用线阵模型原理示意图；

图2为现有技术中由舰船拖曳的常规线列阵示意图；

图3为本发明的基于条件波数谱密度的宽带信号方位估计的方法流程图；

图4为本发明的基于条件波数谱密度的宽带信号方位估计的方法采用的条件波数谱密度图及其沿波数k方向的扩展图；

图5为本发明的基于条件波数谱密度的宽带信号方位估计的方法采用的条件波数谱密度的净化谱图；

图6为采用本发明的基于条件波数谱密度的宽带信号方位估计的方法和采用其他算法得到的方位谱输出对比；

图7-1为采用本发明的基于条件波数谱密度的宽带信号方位估计的方法和采用其他算法算法时间复杂度随阵元数变化情况；

图7-2为采用本发明的基于条件波数谱密度的宽带信号方位估计的方法和采用其他方法时间复杂度随频带划分个数变化情况；

图8为采用宽带最小方差无畸变算法得到的方位历程图；

图9为采用本发明的基于条件波数谱密度的宽带信号方位估计的方法得到的方位历程图。

具体实施方式

下面结合附图和具体实施例对本发明的基于条件波数谱密度的宽带信号方位估计的方法进行详细的说明。

本发明提出的基于条件波数谱密度的宽带信号方位估计的方法直接利用宽带目标信号的特性，根据直线检测原理，将宽带目标信号在条件波数谱密度上呈现的直线转化为方位-能量谱参数域，然后利用最小二乘拟合对估计的方位进行修正，得到修正方位的能量输出，从而快速精确得到远距离目标的方位。

为了实现上述目的，本发明所述的方法所采用的一种线列阵声纳装置，可以是一条拖曳阵，也可以是舷侧阵，所述线列阵声纳由多个水听器组成。设实际阵元数目M，阵元间距d；目标入射方向θ，第m个阵元接收的信号表示为ψ[d_m,f]；快拍长度为T。

本发明包括如下步骤：

步骤1)用空间线阵接收远距离目标的宽带信号，将所述宽带信号在时域和空域做快速傅利叶变换，生成时空二维功率谱密度函数，再将时空二维功率谱密度与阵列接收响应的频谱之比作为条件波数谱密度，根据生成的时空二维功率谱密度函数得到条件波数谱密度，对扩展的条件波数谱密度谱图净化，在条件波数谱的净化谱上，宽带信号峰值能量呈与入射方位角相关的条带状线性分布；

步骤2)在条件波数谱的净化谱上选取估计的扫描方位角，获得该估计的扫描方位角对应的条带内峰值能量位置集，采用最小二乘拟合对估计的扫描方位角进行修正，得到修正方位角，将修正方位角对应的峰值能量作为输出，从而得到算法的方位谱图，方位谱图上峰值位置对应于所述宽带信号的方位。

所述步骤1)具体包括：

步骤1-1)用线阵接收空间信号，得到M个阵元的时域信号ψ[d_m,t]，在时域将数据分为N个快拍，记为快拍1,2,…,n,…,N，N取自然数，每个快拍ψ_n[d_m,t]长度为T，m为阵元序号，t为时间，d_m＝m*d为第m阵元的位置，d为阵元间距；阵列模型如图1所示。

步骤1-2)对第n个快拍数据在时间域上做快速傅利叶变换，得到目标辐射信号时域频谱ψ[d_m,f]：

T为采样点数，取值为自然数，f表示频率，ψ[d_m,t]为第m阵元接收的时域信号；

步骤1-3)确定目标辐射信号的频带范围[f_min,f_max]，f_min为频率最小值，f_max为频带最大值，对频带范围内的每一个频点的空域数据，在空域做变换点数为Q的快速傅利叶变换，得到时空二维频谱Φ[k,f]：

所述的空域快速傅利叶变换点数Q为阵元数M的整数倍，当Q＞M时，采用空间域补零法，再对其在空间域作FFT变换，生成时空二维功率谱密度函数S[k,f]：

其中，E[·]表示数学期望；k为波数，取值为自然数。

所述波数k为：

其中f_i表示位置i处的频率，d是阵元间距，c为声波传播速度，θ是扫描方位角，取值范围0～180度。

将所述的目标辐射信号ψ[d_m,f]做空域FFT变换，能取得更好的方向分辨率和技术效果。

所述的空间域补零法，具体采用直接在原数据后直接补零，或在原数据中插值补零；只要补零后的数据长度满足要求即可。

空间域上补零的目的是保证波数数目Q远大于阵元数M，当波数数目Q足够大时，θ_s能够平均分布在不同的k值上，这样根据空间域FFT运算得到的波束输出也就越精确。

单信号的时空二维功率谱密度为：

步骤1-4)参考阵元功率谱密度S(f)为：

其中，E[·]表示数学期望，上标*表示复共轭。

利用定义求解条件波数谱密度函数S[k|f]，并将零频分量移至谱中心：

S[k|f]＝S[k,f]/S(f) (7)

单信号的条件波数谱密度函数为：

步骤1-5)沿波数k方向对条件波数谱密度函数进行周期延拓，即沿波束k方向向左向右同时复制多条条件波数谱密度函数，得到扩展的条件波数谱密度谱；

步骤1-6)仅选取扩展的条件波数谱密度函数中局部极大值点，将其余点能量置零，计算条件波数谱的净化谱；

所述的仅选取扩展的条件波数谱密度函数中局部极大值点，将其余点能量置零，净化了条件波数谱，降低了S[k|f]非主轴外的能量对方位谱输出的影响，减小了主瓣宽度，同时降低了非目标方位的能量输出，提高了主副瓣比。

所述步骤2)具体包括：

步骤2-1)采用直线扫描思想对条件波数谱密度进行方位谱扫描，在方位角范围(0～180度)内，设定初始扫描方位角θ＝θ₀，并计算扫描斜率μ：

步骤2-2)以过(f,k)＝(0,0)、斜率为μ的直线l₀为中心，宽度2Q/M+1的条带L作为参数域方位，计算局部峰值能量位置集P＝{(k_i,f_i)}；k_i为第i个位置的波数，f_i为第i个位置的频率；

步骤2-3)对条带L内条件波数谱密度函数S[k|f]的局部峰值能量位置集P＝{(k_i,f_i)}做最小二乘直线拟合，将L内净化谱能量累积P作为修正方位角的能量输出，得到修正斜率μ′，拟合过程为：

根据修正斜率μ′，得到修正波达方位cosθ_s：

θ_s为修正方位角；

步骤2-4)：改变扫描方位角θ作为估计的扫描方位角，θ为按照指定的顺序选取的方位角：

θ＝θ₀+lΔθ (12)

θ₀为初始扫描方位角，l为第l次扫描，取值为自然数，Δθ为方位角扫描间隔，取值为0～180度，数值越小，扫描越精细。

判断所述估计的扫描方位角θ是否在有效声锥范围内，如果“是”，返回步骤2-1)继续执行；

所述有效声锥范围，为条件波数谱密度S[k|f]上斜率

所夹的一个锥形区域内，即斜率μ满足

如果“否”，则对条件波数谱密度函数在有效声锥范围内的扫描斜率完成遍历，执行步骤2-5)：

步骤2-5)以每一个修正方位角为横坐标，修正方位角对应幅度值为纵坐标绘制方位谱曲线，峰值处为能量输出的最大值，峰值处对应修正方位角为精确的宽带信号方位角，从而实现对目标方位的精确定位。

实施例

下面结合某次海试数据和附图对本发明的具体实施方式做进一步的详细描述。

试验参数：拖曳阵水听器数目M＝24，水听器间距d＝0.1875m；信号采样率fs＝20000Hz。目标频带范围：1k～4k Hz，方位125度，声速c＝1516m/s，快拍长度N＝20000。

采用基于条件波数谱密度的宽带信号方位估计的方法，具体步骤如下：

步骤S1)：如图3所示，其中的301、302、303，用线阵接收空间信号，得到M个阵元的时域信号x_m[t]，取第l帧快拍的信号在时间域上做20000点FFT运算，得到快拍信号在不同频率的响应。如下式，N表示快拍长度，行表示时间采样，列表示阵元：

步骤S2)：如图3所示，其中对应的304，确定目标辐射信号的频带范围[f_min,f_max]，对频带范围内的每一个频点的空域数据，在空域做Q点快速傅利叶变换如公式(2)：

当Q＞M时，采用空间域上补零，再对其在空间域作FFT变换，生成时空二维功率谱密度函数S[k,f]如公式(3)：

Φ[k,f]为阵列的时空二维频谱，E[·]表示数学期望。

需要说明的是，目标频带范围[f_min,f_max]是1k～4k Hz，采样率20000Hz，FFT长度20000点，那么目标信号对应的离散频点是：(1000～4000)/20000*20000＝(1000～4000)点，因此整个过程只需要对这一段频率范围进行处理即可；空域FFT点数Q＝16·M。

空间域上补零的目的是保证波数数目Q远大于阵元数M，使得空间域FFT运算结果和阵列波束输出的对应关系更加均匀，波束输出也就越精确。空间域补零可以有多种选择。既可以直接在原数据后直接补零，也可以在原数据中插值补零，只要补零后的数据长度满足要求即可。本实例直接在原数据后延长补零。

步骤S3)：如图3所示，其中对应的305，求解条件波数谱密度函数S[k|f]，并将零频分量移至谱中心如公式(7)：

S[k|f]＝S[k,f]/S(f)

具体做法是，首先沿波数k方向对时空二维功率谱积分，求得参考阵元功率谱密度S(f)；后利用公式(3)和公式(6)计算条件波数谱密度S[k|f]；

采用时空二维FFT进行时空二维功率谱密度函数求解时，如公式(5)所示，单信号的时空二维功率谱密度S[k,f]：

将零频分量移至谱中心是一种常用的处理方法，在许多文献中都有说明，对本领域技术人员来说，理解和实现该处理是可以胜任的。

步骤S4)：如图3所示，其中对应的306，沿波数k方向对条件波数谱密度函数进行周期延拓，得扩展的条件波数谱密度；周期延拓方式为沿波束k方向向左向右同时适当复制条件波数谱密度，如图4所示；

如公式(8)所示，单信号条件波数谱密度函数为：

步骤S5)：仅选取扩展的条件波数谱密度函数中局部极大值点，将其余点能量置零，计算条件波数谱的净化谱；净化结果如图5所示。

步骤S6)：选择初始扫描方位角，根据公式(9)计算扫描斜率μ，以过(f,k)＝(0,0)、斜率为μ的直线l₀为中心，宽度2Q/M+1的条带L作为参数域方位，计算局部峰值能量位置集P＝{(k_i,f_i)}；

步骤S7)：对条带L内S(k|f)局部峰值能量位置集P＝{(k_i,f_i)}做最小二乘直线拟合得到修正斜率和修正波达方位，并将此方位作为修正的扫描方位角，实现方位校正，将L内净化谱能量累积作为修正方位的能量输出，拟合过程如公式(10)：

步骤S8)：改变扫描方位角，并根据公式(9)重新计算扫描斜率，重复步骤5、6、7，进行K次波束扫描，直到对有效声锥范围完成遍历；

在本发明的实施例中，步骤S8)扫描方位角是按一定的间隔角度顺序选取的，即将有效声锥范围均分为K段，按从小到大或从大到小的顺序，依次选取一个角度作为扫描角。

如图6所示，对比了采用本发明和其他算法，海试数据的波束输出。蓝色线表示采用本发明得到的波束输出结果；红色线表示采用宽带最小方差无畸变算法得到的波束输出结果，黄色线代表采用常规波束形成算法得到的波束输出结果。从图中可以清楚看出，采用本发明得到的目标波峰的半功率波束宽度远小于其他2中宽带方位谱估计算法，也就说本发明能够有效提高目标方位的分辨率。

为了说明本发明的实时处理优势，专门统计了matlab程序处理海试数据需要的时间长度。信号采样率fs＝20000Hz，数据快拍长度为20000点，直接对24阵元信号做方位谱估计。

采用常规波束形成算法和宽带最小方差无畸变算法和本发明提出的CWSD-based算法时间复杂度随阵元数变化如图7-1所示和频段划分个数变化如图7-2所示的情况。可以看出，采用本发明在阵元数少于25元时，算法耗时低于0.15s,远低于宽带最小方差无畸变算法的1.55s,完全可以实时处理海试数据。

如图8所示，是直接对24阵元信号做宽带最小方差无畸变算法的方位历程图，如图9所示，是采用本发明得到的海试数据方位历程图，可以看出，本发明得到的目标波束非常精细，更加有利于检测远距离微弱目标。

本发明提出基于条件波数谱密度(the Conditional Wavenumber SpectralDensity based，CWSD-based)的宽带信号方位估计的方法，所述方法利用条件波数谱密度函数将阵列信号转换到f-k空间，宽带目标信号能量在该空间呈与入射角相关的线性分布，通过直线检测原理，实现具有超低旁瓣和高分辨性能的DOA估计；能够快速实时得到高分辨的目标方位功率输出。

本发明所述的基于条件波数谱密度的宽带信号方位估计(CWSD-based)方法，所述方法对宽带信号方位的估计精度和分辨能力明显优于常规波束形成方法和宽带最小方差无畸变算法，对阵元位置误差具有很好的宽容性。所述方法无需预估角度和信源数目等先验信息，当宽带处理频带最高频为f时，理论分辨率可达2c/Mfd，均方估计误差小，稳健性高，并具有时间复杂度低、运算效率高的优点，便于实时处理。

最后所应说明的是，以上实施例仅用以说明本发明的技术方案而非限制。尽管参照实施例对本发明进行了详细说明，本领域的普通技术人员应当理解，对本发明的技术方案进行修改或者等同替换，都不脱离本发明技术方案的精神和范围，其均应涵盖在本发明的权利要求范围当中。

Claims (9)

1.一种基于条件波数谱密度的宽带信号方位估计的方法，包括：

步骤1)用空间线阵接收远距离目标的宽带信号，将所述宽带信号在时域和空域做快速傅利叶变换，生成时空二维功率谱密度函数，再将时空二维功率谱密度函数与参考阵元功率谱密度之比作为条件波数谱密度，对扩展的条件波数谱密度谱图净化，得到条件波数谱的净化谱；

步骤2)在条件波数谱的净化谱上选取估计的扫描方位角，计算该估计的扫描方位角对应的条带内峰值能量位置集，采用最小二乘拟合对估计的扫描方位角进行修正，得到修正方位角，将修正方位角与对应的峰值能量作为输出，从而得到方位谱图，方位谱图上的峰值位置对应于所述宽带信号的方位；

所述步骤2)具体包括：

步骤2-1)在方位角范围0～180度内，设定初始扫描方位角θ＝θ₀，并计算扫描斜率μ：

其中c为声波传播速度；

步骤2-2)以过(f,k)＝(0,0)、斜率为μ的直线l₀为中心，宽度2Q/M+1的条带L作为参数域方位，计算局部峰值能量位置集P＝{(k_i,f_i)}；k_i为第i个位置的波数，f_i为第i个位置的频率；

步骤2-3)对条带L内条件波数谱密度函数S[k|f]的局部峰值能量位置集P＝{(k_i,f_i)}做最小二乘直线拟合，得到修正斜率μ′，拟合过程为：

根据修正斜率μ′，得到修正波达方位cosθ_s：

θ_s为修正方位角，c为声波传播速度；d为阵元间距，Q为空域快速傅利叶变换点数；

步骤2-4)改变扫描方位角θ作为估计的扫描方位角，θ为按照指定的顺序选取的方位角；

θ＝θ₀+lΔθ (12)

θ₀为初始扫描方位角，l为第l次扫描，取值为自然数，Δθ为方位角扫描间隔，取值为0～180度；

判断所述估计的扫描方位角θ是否在有效声锥范围内，如果“是”，返回步骤2-1)继续执行；如果“否”，则对条件波数谱密度函数在有效声锥范围内的扫描斜率完成遍历，

执行步骤2-5)；

步骤2-5)以每一个修正方位角为横坐标，修正方位角对应幅度值为纵坐标绘制方位谱曲线，曲线峰值处为能量输出的最大值，峰值处对应修正方位角为精确的宽带信号方位角，从而实现对目标方位的精确定位。

2.根据权利要求1所述的基于条件波数谱密度的宽带信号方位估计的方法，其特征在于，所述步骤1)具体包括：

步骤1-1)用空间线阵接收宽带信号，得到M个阵元的时域信号ψ[d_m,t]，在时域将数据分为N个快拍，记为快拍1,2,…,n,…,N，每个快拍ψ_n[d_m,t]长度为T，其中，m为阵元序号，取值为自然数，t为采样时间点，d_m＝m*d为第m个阵元的位置，d为阵元间距；

步骤1-2)对第n个快拍数据在时间域上做快速傅利叶变换，得到目标辐射信号时域频谱ψ[d_m,f]；

其中，T为采样点数，取值为自然数，f表示频率，ψ_n[d_m,t]为第m个阵元接收的时域信号；

步骤1-3)确定目标辐射信号ψ[d_m,f]的频带范围[f_min,f_max]，f_min为频率最小值，f_max为频带最大值，对频带范围内的每一个频点的空域数据，在空域做快速傅利叶变换，快速傅利叶变换点数为Q，得到时空二维频谱Φ[k,f]：

当Q＞M时，Q和M取值为自然数，采用空间域上补零，Φ[k,f]在空间域作快速傅利叶变换，生成时空二维功率谱密度函数S[k,f]：

其中，E[]表示数学期望，k为波数，取值为自然数；

步骤1-4)参考阵元功率谱密度S(f)为：

其中，上标*表示复共轭；

将零频分量移至谱中心，并通过时空二维功率谱密度函数S[k,f]与参考阵元功率谱密度S(f)之比得到条件波数谱密度函数S[k|f]：

S[k|f]＝S[k,f]/S(f) (7)；

步骤1-5)沿波数k方向对条件波数谱密度函数S[k|f]进行周期性延拓，得到扩展的条件波数谱密度函数；

步骤1-6)选取扩展的条件波数谱密度函数中局部振幅极大值点为宽带信号峰值能量，将其余点置零，计算所述条件波数谱的净化谱；在条件波数谱的净化谱上，宽带信号峰值能量位置呈与入射方位角相关的线性分布。

3.根据权利要求2所述的基于条件波数谱密度的宽带信号方位估计的方法，其特征在于，所述波数k为：

其中f表示频率，d是阵元间距，c为声波传播速度，θ是扫描方位角，取值范围0～180度。

4.根据权利要求2所述的基于条件波数谱密度的宽带信号方位估计的方法，其特征在于，所述步骤1-3)的采用空间域上补零，具体包括采用直接在原数据后补零，或在原数据中插值补零。

5.根据权利要求2所述的基于条件波数谱密度的宽带信号方位估计的方法，其特征在于，所述步骤1-3)中，单信号的时空二维功率谱密度为：

其中，c为声波传播速度，θ为扫描方位角，θ取值范围为0～180度。

6.根据权利要求2所述的基于条件波数谱密度的宽带信号方位估计的方法，其特征在于，所述步骤1-3)的傅利叶变换点数Q为阵元数M的整数倍。

7.根据权利要求2所述的基于条件波数谱密度的宽带信号方位估计的方法，其特征在于，步骤1-4)中单信号的条件波数谱密度函数为：

其中c为声波传播速度，θ为扫描方位角，θ取值范围为0～180度。

8.根据权利要求2所述的基于条件波数谱密度的宽带信号方位估计的方法，其特征在于，所述步骤1-5)的沿波数k方向对条件波数谱密度函数S[k|f]进行周期性延拓为沿波束k方向向左向右同时复制多条条件波数谱密度函数S[k|f]。

9.根据权利要求1所述的基于条件波数谱密度的宽带信号方位估计的方法，其特征在于，所述有效声锥范围为条件波数谱密度S[k|f]上斜率μ的一个锥形区域内，即斜率μ满足：

其中c为声波传播速度。

Images