CN111665489B - 一种基于目标特性的线谱提取方法 - Google Patents

Abstract

本发明提出一种基于目标特性的线谱提取方法，所述方法将宽带目标信号转换生成时空二维功率谱密度函数，通过对频域功率谱密度函数的处理，提取频域功率谱线信息，得到目标信号强线谱；然后在去强线谱影响的时空功率谱密度上选取估计的扫描方位角，计算条带内能量分布沿波数向的均值，得到该扫描方位目标信号的功率谱估计图；在所述功率谱估计图上对各个方位的线谱信息进行提取，并将得到的扫描方位线谱和目标信号强线谱在方位向复合，得到特定目标的线谱特性。本发明的方法无需事先预知目标功率谱信息；适用范围广，可囊括强线谱成分和弱目标线谱成分的提取；速度快，实时性好，算法稳健性高；仿真和海试数据证明本发明有效。

Description

一种基于目标特性的线谱提取方法

技术领域

本发明属于声纳数字信号处理领域，特别涉及一种基于目标特性的线谱提取方法。

背景技术

利用舰船辐射噪声对目标进行远程探测是水声领域重要的研究内容。水声工作者的分析表明，舰船噪声中含有大量的声纹信息，以唯一且稳定存在、传播距离远和可检测性高为特点的低频线谱声纹特征为舰船目标远程探测、目标识别提供了一条可行的思路。

舰船螺旋桨转动，船尾不均匀伴流场中空泡崩溃的周期性与准周期性成分以及机械震动等产生调制包络，尾流沿周向变化产生按叶片频率强幅度调制效应，而桨轴安装的细微偏差可能产生轴频调制，这些调制使得舰船具有明显的节奏信息。由于舰船噪声本身是一种非平稳过程，大量文献分析中，通常将其近似为一种循环局部平稳过程，文献1[舰船噪声节奏的研究(I)——数学模型及功率谱密度.声学学报，1983；8(2):65-76]中陶笃纯采用具有随即幅度、相同形状、相同重复周期的高斯脉冲模型对舰船噪声进行建模，分析了叶频和轴频的成因，并从线谱的频域、时域信息中给出了线谱的典型形状，文献2[舰船噪声识别(I)—总体框架、线谱分析和提取.声学学报.1998；23(5):394:400]中吴国清等人利用双重谱的幅度和相位信息重建了周期调制信号的波形，并通过设定线谱识别规则将双重谱去趋势、线谱净化等方式提取线谱特征，谢俊等从线谱是叶频和轴频的谐波分量对宽带随机噪声的调制出发，提出了多子带融合算法，利用各个频带的融合信息的高阶统计量进行自适应线谱增强，实现线谱的提取；文献3[数字式声呐设计原理.合肥:安徽教育出版社,2002]中李启虎等人通过分析大量舰船的目标特性，提出了不变特征量用于目标识别。目前，小波理论、高阶统计量、神经网络等提取目标线谱特征提取。以上算法在线谱提取前均已假定目标信号已完全提取，利用信号频谱进行线谱特征提取，但对目标信号的提取主要依赖于波束扫描加权的方式，而实际应用中对信号信息知之甚少。

频率-波数功率谱分析(f-k分析)是地震波中求解波传播速度和方位角的一种处理方法。文献4[High-Resolution Frequency-Wavenumber SpectrumAnalysis.Proceedingsof the IEEE.1969]中J.Capon.将阵列采样看作是空间平稳随机场采样，利用时空二维相关函数的傅氏变换得到时空功率谱密度，用以描述阵列接收的宽平稳地震波信号；f-k分析方法已广泛应用于地震波传播速度测定，文献5[王芳,王宝善.频率-波数域方法的发展及其在台阵数据分析中的应用.中国地震.2017；2:191-202]、表面波的损伤检测文献6[Tian Z,Yu L.Lamb wave frequency–wavenumber analysis anddecomposition.Journal of Intelligent Material Systems and Structures,2014；25(9):1-17]和表面波色散问题文献[7Poggi V.Time–Frequency–Wavenumber Analysis ofSurface Waves Using the Continuous Wavelet Transform.Pure AppliedGeophysics.2013；170(3):319–335]研究，在声呐领域，文献8[Flow noise analysis oftowed sonar arrays.Udt Europe.1999]中Beerens P等人利用f-k分析了拖曳阵流噪声的成因，取得了良好的效果；而利用时空二维功率谱密度在频率-波数参数空间的能量分布特征进行线谱提取的研究较少。

发明内容

本发明目的在于，为克服现有的线谱提取算法假定目标信号已知，依赖于波束扫描加权方式，而实际应用中对信号信息知之甚少的缺点。

为实现上述目的，本发明提出一种基于目标特性的线谱提取方法，利用信号特征在时空二维功率谱上的分布信息，对目标线谱特征进行提取，并在方位谱图上呈现，仿真和海试数据验证了算法的有效性。

本发明提出的一种基于目标特性的线谱提取方法，所述方法一种基于目标特性的线谱提取方法，所述方法将宽带目标信号转换生成时空二维功率谱密度函数，通过对频域功率谱密度函数的处理，提取频域功率谱线信息，得到目标信号强线谱；然后在去强线谱影响的时空功率谱密度上选取估计的扫描方位角，计算条带内能量分布沿波数向的均值，得到该扫描方位目标信号的功率谱估计图；在所述功率谱估计图上对各个方位的线谱信息进行提取，并将得到的扫描方位线谱和目标信号强线谱在方位向复合，得到特定目标的线谱特性。作为所述方法的一种改进，具体包括：

步骤1)将空间线阵接收的远距离目标信号，在时域和空域做快速傅利叶变换，生成时空二维功率谱密度函数，再对时空二维功率谱密度沿空间频率方向积分得频域功率谱密度；

步骤2)通过对频域功率谱密度进行平滑处理得到整个接收平面的趋势项，对频域功率谱密度进行去趋势后，计算各频域功率谱线的高度，再通过设置的第一阈值提取频域功率谱线局部最大值，对伪峰进行剔除，得到目标信号强线谱；

步骤3)对所述目标信号强线谱的频率置零，得到去强线谱影响的时空功率谱密度，在去强线谱影响的时空功率谱密度上选取估计的扫描方位角，计算该估计的扫描方位角对应的条带，计算条带内能量分布沿波数向的均值，得到目标信号的功率谱估计；对所述目标信号的功率谱估计通过去趋势和设置的第二阈值进行方位线谱提取，对伪峰进行剔除，得到所述扫描方位线谱；在一定的角度容差范围内，对所述扫描方位线谱和目标信号强线谱在方位向复合，得到特定目标的线谱特性。

作为所述方法的一种改进，所述步骤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个快拍数据在时间域上做快速傅利叶变换，得到目标信号时域频谱ψ_n[d_m,f]：

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

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

采用空间域上补零，Φ[k,f]在空间域作快速傅利叶变换，生成时空二维功率谱密度函数S[k,f]：

当目标信号为单信号时，单信号的时域频谱

Ψ(f)为单信号的频谱，则目标信号为单信号时的时空二维功率谱密度为：

其中，E[·]表示数学期望，k为波数，取值为自然数，上标H表示复共轭；波数k为：

其中f表示频率，θ是扫描方位角，取值范围0～180度；

步骤1-4)在时空二维功率谱密度函数S[k,f]上，将零频分量移至谱中心，频率f下的空间互谱密度函数H[χ,f]为：

由S(f)＝H[0,f]得到频域功率谱密度S(f)为：

目标信号为单信号时的频域功率谱密度为：

在时空二维功率谱密度S[k,f]上，宽带信号峰值能量位置呈与入射方位角相关的线性分布。

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

步骤2-1)对频域功率谱密度S(f)取P个数据点，进行平滑处理，得到整个接收平面的趋势项，则趋势项F_trend(f)为：

F_trend(f)＝S(f)*H(f) (10)

符号*为卷积运算，H(f)为低通滤波器，包括平滑低通滤波器：

步骤2-2)在频域功率谱S(f)上去趋势项，得到去趋势的频域功率谱密度L(f)：

L(f)＝S(f)-F_trend(f) (12)

步骤2-3)计算去趋势的频域功率谱密度L(f)中各谱线的高度，利用设置的第一阈值对去趋势的频域功率谱线局部最大值进行提取，得到目标信号强线谱；

步骤2-4)在目标信号强线谱中，找出频率非零的最大峰值，用此作为标准，在一定容差范围内计算其他峰值与其的最大公约数，由结果判断是否具有谐波关系，若无则剔除，经过筛选得到S个强线谱频率f_i，i＝1,2,...S；得到对应的目标信号强线谱频率集为Fre＝{f_i}；

步骤2-5)对于目标频域功率谱强线谱频率集Fre中每一个强线谱频率f_i，从S[k,f_i]曲线中提取超过设定固定阈值的局部最大值，得到该最大值对应的坐标(k,f_i)，求解强线谱频率f_i对应的目标方位θ_r，满足：

sinθ_r＝kc/f_iQd (13)。

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

步骤2-3-1)计算去趋势的频域功率谱密度L(f)中各谱线的高度为：去趋势谱线中的局部最大值谱峰高度；

步骤2-3-2)设置的第一阈值包括：阈值A和阈值B：

a)设所有谱线高度的平均值的2倍为阈值A；

b)设包络谱能量最高点的谱强度为S₀，基线强度S_a，则阈值B为0.3(S₀-S_a)；

步骤2-3-3)若谱线的高度大于等于阈值A或阈值B之一，则判定为目标信号频域功率谱强线谱的局部最大值，基于dB谱进行频率局部最大值提取，得到目标信号强线谱。

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

步骤3-1)根据目标信号强线谱集Fre中每一个频率f₀，将时空二维功率谱上任意波数k对应的S[k,f₀]置零，得到去强线谱时空功率谱密度S′[k,f₀]；

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

步骤3-3)以过(f,k)＝(0,0)、斜率为μ的直线l₀为中心，宽度2Q/M+1作条带L；对条带L内的去强线谱时空功率谱密度S′[k,f]沿波数k方向求均值，得θ＝θ₀方向的信号功率谱估计E(f|l₀,θ₀)；

步骤3-4)对信号功率谱估计E(f|l₀,θ₀)取R个数据点进行平滑处理，得到整个接收平面的趋势项，则趋势项F^E _trend(f)：

F^E _trend(f)＝E(f|l₀,θ₀)*H'(f) (15)

符号*为卷积运算，H′(f)为低通滤波器，包括平滑低通滤波器：

步骤3-5)在θ₀方向的信号功率谱估计E(f|l₀,θ₀)上去趋势项，得到去趋势的信号功率谱密度L′(f)：

L′(f)＝E(f|l₀,θ₀)-F^E _trend(f) (17)

步骤3-6)计算去趋势的信号功率谱密度L′(f)中的各谱线的高度，根据设置的第二阈值进行提取方位线谱；

步骤3-7)在方位线谱的超过第二阈值的谱峰中，找出频率非零的最大峰，用此峰值作为标准，在一定容差范围内求解其他峰与其的最大公约数，由结果判断是否具有谐波关系，若无则剔除，经过筛选得到对应的X个线谱频率f_j，j＝1,2,...X；获得方向的去强线谱的信号功率谱所对应扫描方位线谱频率集为F_sig＝{f_j}；

步骤3-8))对于扫描方位线谱频率集Fre中每一个强线谱频率f_j′，求解f_j′对应的方位θ_r′；满足：

sinθ′_r＝kc/f′_jQd

θ_r′为估计的方位角θ₀所对应的准确方位；改变扫描方位角θ作为估计的扫描方位角，θ为按照指定的顺序选取的方位角；

θ＝θ₀+lΔθ (18)

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

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

如果“否”，则对时空二维功率谱在有效声锥范围内的扫描斜率完成遍历；

将每一个方位角θ_r′对应的方位线谱复合得到扫描方位线谱，执行步骤3-10)；

步骤3-10)在给定的角度容差范围内，对步骤2)得到的目标信号强线谱和扫描方位线谱在方位向复合，得到特定目标的线谱特性。

作为所述方法的一种改进，步骤3-6)包括：

步骤3-6-1)计算去趋势的L′(f)中各谱线的高度为：去趋势谱线中的局部最大值谱峰高度；

步骤3-6-2)设置第二阈值包括：阈值C和阈值D：

a)设所有谱线高度的平均值的2倍为阈值C；

b)设包络谱能量最高点的谱强度为S₀，基线强度S_a，则阈值D为0.3(S₀-S_a)；

步骤3-6-3)若谱线的高度大于等于阈值C或阈值D之一，则判定为方位谱线局部最大值，基于dB谱进行方位谱线局部最大值提取，得到方位线谱。

作为所述方法的一种改进，步骤3-6)还包括：设置第二阈值为一固定值，取值为6～10dB或为目标信号的所有频率对应的L′(f)方差的2～3倍。

作为所述方法的一种改进，步骤3-9)的有效声锥范围为时空二维功率谱密度S[k,f]上斜率μ的一个锥形区域内，斜率μ满足：

本发明的优势在于：

1、本发明的基于目标特性的线谱提取方法利用时空二维功率谱密度函数将宽带目标信号转换到f-k空间，宽带目标信号能量在该空间呈与入射角相关的线性分布，利用信号特征在时空二维功率谱上的分布特性，对给定方位的目标功率谱信息进行提取，然后对目标线谱特征进行提取，并在方位谱图上呈现，得到目标线谱特性的方位谱估计的同时，对各个方位的线谱信息进行提取。

2、本发明的基于目标特性的线谱提取方法无需事先预知目标功率谱信息；

3、本发明的基于目标特性的线谱提取方法适用范围广，可囊括强线谱成分和弱目标线谱成分的提取；

4、本发明的基于目标特性的线谱提取方法采用两重快速傅立叶变换，对信号频带范围内的时空二维功率谱密度进行直线扫描处理，速度快，实时性好，算法稳健性高；

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

附图说明

图1为本发明的基于目标特性的线谱提取方法所用线阵模型原理示意图；

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

图3为本发明的基于目标特性的线谱提取方法流程图；

图4为本发明的基于目标特性的线谱提取方法采用的频域功率谱示意图；

图5(a)为本发明的基于目标特性的线谱提取方法得到的去趋势频域功率谱示意图；

图5(b)为本发明的基于目标特性的线谱提取方法得到的去趋势频域功率谱阈值选取方式；

图6为采用本发明的基于目标特性的线谱提取方法得到线谱提取结果；

图7为采用本发明的基于目标特性的线谱提取方法得到时空二维功率谱密度上线谱检测结果；

图8(a)为海上试验20～100Hz方位历程图；

图8(b)为采用本发明的基于目标特性的线谱提取的线谱提取海上试验线谱方位历程图。

具体实施方式

下面结合附图和具体实施例对本发明的基于目标特性的线谱提取方法进行详细的说明。

本发明提出的基于目标特性的线谱提取方法利用时空二维功率谱密度函数将宽带目标信号转换到f-k空间，宽带目标信号能量在该空间呈与入射角相关的线性分布，利用信号特征在时空二维功率谱上的分布特性，对给定方位的目标功率谱信息进行提取，然后对目标线谱特征进行提取，并在方位谱图上呈现，得到目标线谱特性的方位谱估计的同时，从而快速对各个方位的线谱信息进行提取。

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

本发明包括如下步骤：

步骤1)用空间线阵接收远距离目标的宽带信号，将所述宽带信号在时域和空域做快速傅利叶变换，生成时空二维功率谱密度函数，再对时空二维功率谱密度沿空间频率方向积分得频域功率谱密度；

步骤2)通过对频域功率谱密度平滑得到整个接收平面的趋势项，对频率功率谱密度进行去趋势后得到去趋势的频域功率谱密度，提取去趋势后的频域功率谱密度的局部最大值，通过给定的阈值的方式进行线谱提取，通过谐波知识对伪峰进行剔除，得到信号强线谱频率成分；

步骤3)对所述目标信号强线谱的频率置零，得到去强线谱影响的时空功率谱密度，在去强线谱影响的时空功率谱密度上选取估计的扫描方位角，计算该估计的扫描方位角对应的条带，计算条带内能量分布沿波数向的均值，得到目标信号的功率谱估计；对所述目标信号的功率谱估计通过去趋势和设置的第二阈值进行方位线谱提取，对伪峰进行剔除，得到所述扫描方位线谱；在一定的角度容差范围内，对所述扫描方位线谱和目标信号强线谱在方位向复合，得到特定目标的线谱特性。

所述步骤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为波数，取值为自然数,上标H表示复共轭；

并将S[k,f]的零频分量移至谱中心。

所述波数k为：

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

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

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

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

作为一个例子，单信号的阵列响应

Ψ(f)为单信号的频谱，则时空二维功率谱密度为：

步骤1-4)频域功率谱密度S(f)可由空间互谱密度函数导出，具体如下：

其中，E[·]表示数学期望，频率f下的空间互谱密度函数

由Winner-Khintchine定理

频域功率谱密度S(f)为：

作为一个例子，单信号频谱为Ψ(f)，利用定义式(8)求解单信号频域功率谱密度为：

所述步骤2)具体包括：

步骤2-1)对频域功率谱密度S[f]平滑得到整个接收平面的趋势项，平滑长度取P个数据点，这一过程类似于设计合适的低通滤波器H(f)，则趋势项

F_trend(f)＝S(f)*H(f) (10)

符号*为卷积运算，H(f)为低通滤波器，包括平滑低通滤波器：

所述的低通滤波器包括但不限于平滑滤波器，滤波器阶数越高，值越均匀，滤波效果越好。

步骤2-2)在频域功率谱S(f)上去趋势项，得到去趋势的频域功率谱密度

L(f)＝S(f)-F_trend(f) (12)

步骤2-3)提取去趋势的频域功率谱密度L(f)中的局部最大值，利用设置的第一阈值采取卡门限的方法进行目标信号强线谱提取；

所述的利用设定的第一阈值的方法进行线谱提取，第一阈值的提取不采用绝对阈值，而采用相对大小，这是考虑到相对阈值有利于检测弱调制的情况。具体方法为：首先求出L(f)中各谱线的高度(即去趋势谱线中的局部最大值谱峰高度)，然后定义第一阈值为两个相对阈值：

a)所有谱线高度的平直至的2倍为阈值A；

b)设包络谱能量最高点的谱强度为S₀，基线强度S_a，阈值B为0.3S₀-S_a)，此阈值的物理意义对应于调制深度。

若谱线高度大于等于阈值A或阈值B之一，则判定为目标信号强线谱。这里所有的处理都是基于已化为dB的谱上进行的。

步骤2-4)利用调制线谱通常具有的谐波关系，在目标信号强线谱中，找出频率非零的最大峰，用此峰值作为标准，在一定容差范围内求解其他峰与其的最大公约数，由结果判断是否具有谐波关系，若无则剔除，经过筛选得到若干个强线谱频率f_i，i为自然数，得到对应的目标信号强线谱频率集为Fre＝{f_i}；

步骤2-5)对于目标频域功率谱强线谱频率集Fre中每一个强线谱频率f_i，从S[k,f_i]曲线中提取超过设定固定阈值的局部最大值，得到该最大值对应的坐标(k,f_i)，求解强线谱频率f_i对应的目标方位θ_r，满足：

sinθ_r＝kc/f_iQd (13)。

其中，Q为空域傅里叶变换点数，d为阵元间距，c为声音在介质中传播速度。

所述步骤3)具体包括：

步骤3-1)根据步骤2)中求得的强线谱频率集Fre中每一个频率f₀，将时空二维功率谱上任意波数k对应的S[k,f₀]置零，得到去强线谱时空功率谱密度S′[k,f]；

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

步骤3-3)以过(f,k)＝(0,0)、斜率为μ的直线l₀为中心，宽度2Q/M+1作条带L，条带L是信号能量分布的一种形式，对条带L内的去强线谱时空功率谱密度S′[k,f]沿波数k方向求均值，得θ＝θ₀方向的信号功率谱，记作E(f|l₀,θ₀)

步骤3-4)对信号功率谱估计E(f|l₀,θ₀)取R个数据点进行平滑处理，得到整个接收平面的趋势项，则趋势项F^E _trend(f)：

F^E _trend(f)＝E(f|l₀,θ₀)*H'(f) (15)

符号*为卷积运算，H′(f)为低通滤波器，包括平滑低通滤波器：

所述的低通滤波器包括但不限于平滑滤波器，滤波器阶数越高，值越均匀，滤波效果越好。

步骤3-5)在θ₀方向的信号功率谱E(f|l₀,θ₀)上去趋势项，得到去趋势的信号功率谱密度

L′(f)＝E(f|l₀,θ₀)-F^E _trend(f) (17)

步骤3-6)提取去趋势的信号功率谱密度L′(f)中的局部最大值，利用设置的第二阈值采用卡门限的方式进行线谱提取；

所述的利用设置的第二阈值进行线谱提取，第二阈值取值包括绝对值阈值6～10dB或所有频率成分方差的2～3倍，此外还有方法为：首先求出L′(f)中各谱线的高度(即去趋势谱线中的局部最大值谱峰高度)，然后定义第二阈值为两个相对阈值：

a)所有谱线高度的平直至的2倍为阈值C；

b)设包络谱能量最高点的谱强度为S₀，基线强度S_a，阈值D为0.3(S₀-S_a)，此阈值的物理意义对应于调制深度。

若谱线高度大于等于阈值C或阈值D之一，则判定为方位谱线局部最大值，基于d B谱进行方位谱线局部最大值提取，得到方位线谱。

步骤3-7)利用方位线谱通常具有的谐波关系，在提取的方位谱线局部最大值中，找出频率非零的最大峰，用此峰值作为标准，在一定容差范围内求解其他峰与其的最大公约数，由结果判断是否具有谐波关系，若无则剔除，获得扫描方位线谱频率集为F_sig＝{f_j}；j为自然数；

步骤3-8))对于扫描方位线谱频率集Fre中每一个强线谱频率f_j′，求解f_j′对应的方位θ_r′；满足：sinθ′_r＝kc/f′_jQd；

θ_r′为估计的方位角θ₀所对应的准确方位；改变扫描方位角θ作为估计的扫描方位角，θ为按照指定的顺序选取的方位角；

θ＝θ₀+lΔθ (18)

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

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

如果“否”，则对时空二维功率谱在有效声锥范围内的扫描斜率完成遍历；所述有效声锥范围为时空二维功率谱密度S[k，f]上斜率μ的一个锥形区域内，即斜率μ满足：

将每一个方位角θ_r′对应的的方位线谱复合得到扫描方位线谱，执行步骤3-10)；

步骤3-10)在给定的角度容差范围内，对步骤2)得到的目标信号强线谱和扫描方位线谱在方位向复合，得到特定目标的线谱特性。

实施例

下面结合某次海试数据和附图对本发明的具体实施方式做进一步的详细描述试验参数：拖曳阵水听器数目M＝32，水听器间距d＝8m；信号采样率fs＝2000Hz。目标频带范围：20～100Hz，方位150度，声速c＝1516m/s，快拍长度N＝2000。

如图3所示，基于目标特性的线谱提取方法的具体步骤如下：

步骤S301)：用线阵接收空间信号，得到M个阵元的时域信号x_m[t]；

步骤S302)：对时域信号x_m[t]取第l帧快拍的信号；

步骤S303)：在时间域上做2000点FFT运算，得到快拍信号在不同频率的响应。如下式，N表示快拍长度，行表示时间采样，列表示阵元：

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

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

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

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

步骤S305)：将S[k,f]零频分量移至谱中心，并求解频域功率谱S(f)，如公式(8)；

采用时空二维FFT进行时空二维功率谱密度函数求解时，如公式(3)所示；

如图4所示S(f)频域功率谱示例，频域功率谱S(f)计算方式如公式(8)；

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

步骤S306)：对频域功率谱密度S[f]平滑得到整个接收平面的趋势项，平滑长度取P个数据点，这一过程类似于设计合适的低通滤波器H(f)，则趋势项如公式(10)；

步骤S307)：在频域功率谱S(f)上去趋势项，得到去趋势的频域功率谱密度如公式(12)；

去趋势项的频域功率谱如图5(a)所示，图中阈值选取如图5(b)所示。

步骤S308)：提取去趋势的频域功率谱密度L(f)中的局部最大值，利用给定的阈值的方式进行线谱提取；

所得结果如图6所示。

步骤S309)：利用调制线谱通常具有的谐波关系，在超过阈值谱峰中，找出频率非零的最大峰，用此峰值作为标准，在一定容差范围内求解其他峰与其的最大公约数，由结果判断是否具有谐波关系，若无则剔除，记所有强线谱信号所对应频率集即目标信号强线谱集为Fre＝{f_i}。强线谱成分结果如图7中圆圈所示，其方位角如图中点画线所示，横轴为频率f，纵轴为波数k。

步骤S310)：对于目标信号强线谱集Fre中每一个强线谱频率f₀，从二维功率谱中提取S[k，f₀]，并在S[k，f₀]曲线上按照一定的阈值值寻找局部最大，利用波数k与频率f的关系，求解该目标信号强线谱可能对应的目标方位；目标方位θ与波数k和频率f的关系如公式(13)：

步骤S311)：根据步骤S304中求得的强线谱频率集Fre中每一个频率f₀，将时空二维功率谱上任意波数k对应的S[k，f₀]置零，得到去强线谱时空功率谱密度S′[k，f]；

步骤S312)：在方位角范围0～180度内，设定初始扫描方位角θ＝θ₀，并计算扫描斜率μ如公式(14)；

步骤S313)：以过(f，k)＝(0，0)、斜率为μ的直线l₀为中心，宽度2Q/M+1作条带L，条带L是信号能量分布的一种形式。对条带L内的去强线谱时空功率谱密度S′[k，f]沿波数k方向求均值，得θ＝θ₀方向的信号功率谱，记作E(f|l₀，θ₀)；

步骤S314)：对E(f|l₀，θ₀)平滑得到整个接收平面的趋势项，平滑长度取P个数据点，这一过程类似于设计合适的低通滤波器H(f)，则趋势项如公式(15)；

步骤S315)：在θ₀方向的信号功率谱E(f|l₀，θ₀)上去趋势项，得到去趋势的信号功率谱密度如公式(17)：

步骤S316)：提取去趋势的信号功率谱密度L′(f)中的局部最大值，利用给定的阈值的方式进行线谱提取；

步骤S317)：利用调制线谱通常具有的谐波关系，在超过阈值谱峰中，找出频率非零的最大峰，用此峰值作为标准，在一定容差范围内求解其他峰与其的最大公约数，由结果判断是否具有谐波关系，若无则剔除，得到扫描方位线谱集为F_sig＝{f¹ _i}；

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

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

判断所述估计的扫描方位角θ是否在有效声锥范围内，如果“是”，返回步骤S312)继续执行；如果“否”，则对时空二维功率谱在有效声锥范围内的扫描斜率完成遍历，执行步骤S319)；

所述有效声锥范围为时空二维功率谱密度S[k,f]上斜率μ的一个锥形区域内，即斜率μ如公式(20)满足：

步骤S319)在一定的角度容差范围内，对上述步骤309)得到的目标信号强线谱和步骤S317)中扫描方位线谱在方位向融合，得到特定目标(特定方位角)的线谱特性。

为了说明本发明的处理效果，专门利用了MATLAB程序处理海试数据。信号采样率fs＝2000Hz，数据快拍长度为2000点，直接利用32阵元信号对150度方位角90Hz合作单频声源进行线谱提取，处理频段20～100Hz内方位谱如图8(a)所示，线谱检测后方位谱如图8(b)所示，从图中可以清晰看出合作单频声源方位被正确检出。

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

Claims (8)

1.一种基于目标特性的线谱提取方法，所述方法将宽带目标信号转换生成时空二维功率谱密度函数，通过对频域功率谱密度函数的处理，提取频域功率谱线信息，得到目标信号强线谱；然后在去强线谱影响的时空功率谱密度上选取估计的扫描方位角，计算条带内能量分布沿波数向的均值，得到目标信号的功率谱估计图；在所述功率谱估计图上对各个方位的线谱信息进行提取，并将得到的扫描方位线谱和目标信号强线谱在方位向复合，得到特定目标的线谱特性；

所述方法具体包括：

步骤1)将空间线阵接收的远距离目标信号，在时域和空域做快速傅利叶变换，生成时空二维功率谱密度函数，再对时空二维功率谱密度沿空间频率方向积分得频域功率谱密度；

步骤2)通过对频域功率谱密度进行平滑处理得到整个接收平面的趋势项，对频域功率谱密度进行去趋势后，计算各频域功率谱线的高度，再通过设置的第一阈值提取频域功率谱线局部最大值，对伪峰进行剔除，得到目标信号强线谱；

步骤3)对所述目标信号强线谱的频率置零，得到去强线谱影响的时空功率谱密度，在去强线谱影响的时空功率谱密度上选取估计的扫描方位角，计算该估计的扫描方位角对应的条带，计算条带内能量分布沿波数向的均值，得到目标信号的功率谱估计；对所述目标信号的功率谱估计通过去趋势和设置的第二阈值进行方位线谱提取，对伪峰进行剔除，得到所述扫描方位线谱；在一定的角度容差范围内，对所述扫描方位线谱和目标信号强线谱在方位向复合，得到特定目标的线谱特性。

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个快拍数据在时间域上做快速傅利叶变换，得到目标信号时域频谱ψ_n[d_m,f]：

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

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

采用空间域上补零，Φ[k,f]在空间域作快速傅利叶变换，生成时空二维功率谱密度函数S[k,f]：

当目标信号为单信号时，单信号的时域频谱

Ψ(f)为单信号的频谱，则目标信号的时空二维功率谱密度为：

其中，E[·]表示数学期望，k为波数，取值为自然数，上标H表示复共轭；波数k为：

其中f表示频率，θ是扫描方位角，取值范围0～180度；

步骤1-4)在时空二维功率谱密度函数S[k,f]上，将零频分量移至谱中心，频率f下的空间互谱密度函数H[χ,f]为：

则，频域功率谱密度S(f)为：

在时空二维功率谱密度S[k,f]上，宽带信号峰值能量位置呈与入射方位角相关的线性分布；

当目标信号为单信号时，目标信号的频域功率谱密度为：

3.根据权利要求2所述的基于目标特性的线谱提取方法，其特征在于，所述步骤2)具体包括：

步骤2-1)对频域功率谱密度S(f)取P个数据点，进行平滑处理，得到整个接收平面的趋势项，则趋势项F_trend(f)为：

F_trend(f)＝S(f)*H(f) (10)

符号*为卷积运算，H(f)为低通滤波器，包括平滑低通滤波器：

步骤2-2)在频域功率谱S(f)上去趋势项，得到去趋势的频域功率谱密度L(f)：

L(f)＝S(f)-F_trend(f) (12)

步骤2-3)计算去趋势的频域功率谱密度L(f)中各谱线的高度，利用设置的第一阈值对去趋势的频域功率谱线局部最大值进行提取，得到目标信号强线谱；

步骤2-4)在目标信号强线谱中，找出频率非零的最大峰值，用此作为标准，在一定容差范围内计算其他峰值与其的最大公约数，由结果判断是否具有谐波关系，若无则剔除，经过筛选得到S个强线谱频率f_i，i＝1,2,...S；得到对应的目标信号强线谱频率集为Fre＝{f_i}；

步骤2-5)对于目标频域功率谱强线谱频率集Fre中每一个强线谱频率f_i，从S[k,f_i]曲线中提取超过设定固定阈值的局部最大值，得到该最大值对应的坐标(k,f_i)，求解强线谱频率f_i对应的目标方位θ_r，满足：

sinθ_r＝kc/f_iQd (13)。

4.根据权利要求3所述的基于目标特性的线谱提取方法，其特征在于，所述步骤2-3)具体包括：

步骤2-3-1)计算去趋势的频域功率谱密度L(f)中各谱线的高度为：去趋势谱线中的局部最大值谱峰高度；

步骤2-3-2)设置的第一阈值包括：阈值A和阈值B：

a)设所有谱线高度的平均值的2倍为阈值A；

b)设包络谱能量最高点的谱强度为S₀，基线强度S_a，则阈值B为0.3(S₀-S_a)；

步骤2-3-3)若谱线的高度大于等于阈值A或阈值B之一，则判定为目标信号频域功率谱强线谱的局部最大值，基于dB谱进行频率局部最大值提取，得到目标信号强线谱。

5.根据权利要求3所述的基于目标特性的线谱提取方法，其特征在于，所述步骤3)具体包括：

步骤3-1)根据目标信号强线谱Fre中每一个频率f₀，将时空二维功率谱上任意波数k对应的S[k,f₀]置零，得到去强线谱时空功率谱密度S′[k,f₀]；

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

步骤3-3)以过(f,k)＝(0,0)、斜率为μ的直线l₀为中心，宽度2Q/M+1作条带L；对条带L内的去强线谱时空功率谱密度S′[k,f]沿波数k方向求均值，得θ＝θ₀方向的信号去强线谱功率谱估计E(f|l₀,θ₀)；

步骤3-4)对信号功率谱估计E(f|l₀,θ₀)取R个数据点进行平滑处理，得到整个接收平面的趋势项，则趋势项F^E _trend(f)：

F^E _trend(f)＝E(f|l₀,θ₀)*H'(f) (15)

符号*为卷积运算，H′(f)为低通滤波器，包括平滑低通滤波器：

步骤3-5)在θ₀方向的信号功率谱估计E(f|l₀,θ₀)上去趋势项，得到去强线谱的信号功率谱密度L′(f)：

L′(f)＝E(f|l₀,θ₀)-F^E _trend(f) (17)

步骤3-6)计算去强线谱的信号功率谱密度L′(f)中的各谱线的高度，根据设置的第二阈值进行提取方位线谱；

步骤3-7)在方位线谱的超过第二阈值的谱峰中，找出频率非零的最大峰，用此峰值作为标准，在一定容差范围内求解其他峰与其的最大公约数，由结果判断是否具有谐波关系，若无则剔除，经过筛选得到对应的X个线谱频率f_j′，j＝1,2,...X；获得方向的去强线谱的信号功率谱所对应扫描方位线谱频率集为F_sig＝{f′_j}；

步骤3-8)对于扫描方位线谱频率集Fre中每一个强线谱频率f_j′，求解f_j′对应的方位θ_r′；满足：

sinθ′_r＝kc/f′_jQd

θ′_r为估计的方位角θ₀所对应的准确方位；改变扫描方位角θ作为估计的扫描方位角，θ为按照指定的顺序选取的方位角；

θ＝θ₀+lΔθ (18)

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

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

如果“否”，则对时空二维功率谱在有效声锥范围内的扫描斜率完成遍历；

将每一个方位角θ′_r对应的方位线谱复合得到扫描方位线谱，执行步骤3-10)；

步骤3-10)在给定的角度容差范围内，对步骤2)得到的目标信号强线谱和扫描方位线谱在方位向复合，得到特定目标的线谱特性。

6.根据权利要求5所述的基于目标特性的线谱提取方法，其特征在于，步骤3-6)包括：

步骤3-6-1)计算去趋势的L′(f)中各谱线的高度为：去趋势谱线中的局部最大值谱峰高度；

步骤3-6-2)设置第二阈值包括：阈值C和阈值D：

a)设所有谱线高度的平均值的2倍为阈值C；

b)设包络谱能量最高点的谱强度为S₀，基线强度S_a，则阈值D为0.3(S₀-S_a)；

步骤3-6-3)若谱线的高度大于等于阈值C或阈值D之一，则判定为方位谱线局部最大值，基于dB谱进行方位谱线局部最大值提取，得到方位线谱。

7.根据权利要求5所述的基于目标特性的线谱提取方法，其特征在于，步骤3-6)还包括：设置第二阈值为一固定值，取值为6～10dB或为目标信号的所有频率对应的L′(f)方差的2～3倍。

8.根据权利要求5所述的基于目标特性的线谱提取方法，其特征在于，步骤3-9)的有效声锥范围为时空二维功率谱密度S[k,f]上斜率μ的一个锥形区域内，斜率μ满足：

Images