一种基于统计量分析的船舶辐射噪声线谱精确提取方法与流程

1.本发明涉及船舶辐射噪声提取技术领域，具体为一种基于统计量分析的船舶辐射噪声线谱精确提取方法。


背景技术：

2.船舶辐射噪声包括机械噪声、螺旋桨噪声、水动力噪声，其经过短时傅里叶变换得到的lofar图上存在丰富的线谱成分，是一种较为稳定的特征信息。这些线谱成分的精确提取对水下目标探测与识别具有重要的工程意义。
3.常用的舰船辐射噪声线谱提取方法，其要经过信号的预处理、短时傅里叶变换得到lofar图，再通过人工手动或者设计阈值提取线谱。但是由于海洋环境过于复杂，有许多其他噪声的干扰，如其他舰船目标的噪声、海洋环境噪声、海洋动物叫声等，都会导致lofar图中的线谱时明时暗，出现断点，甚至出现交叉重叠。这些因素导致通过人工手动或者设计阈值提取出线谱的方法，泛化性弱，可能存在伪线谱，不利于后续的水声目标探测或识别任务的进行，会严重影响水声目标探测或识别的精度。
4.lofar图不同频率位置的能量统计量，如能量最大值、能量最小值、能量平均值、能量中位数、能量极差、能量方差等，能帮助研究人员筛选出有价值的信息，去除混淆在真实线谱中的伪线谱成分。在舰船辐射噪声lofar图中线谱表现为在整个时间段中该频率位置的都具有较高能量值，经过时间上累积后总能量较大，而伪线谱是由于在该频率位置存在某个或某几个时间点有能量特别大的“异值”，导致时间上累积后总能量较大。线谱和伪线谱经过时间累积后能量都较高，通过人工设计阈值或者人工手动提取出线谱的方法无法区分出线谱与伪线谱。但是由于伪线谱是“异值”导致的，故线谱在整个时间段内能量分布比伪线谱在整个时间段内能量分布更为稳定。


技术实现要素：

5.基于上述原因，结合大量实验分析，针对人工设计阈值或者人工手动提取出线谱存在伪线谱的问题，本发明提出一种基于统计量分析的船舶辐射噪声线谱精确提取方法，先提取船舶辐射噪声的lofar图中所有疑似线谱，然后利用6个统计量(最大值、最小值、平均值、中位数、极差、方差)对疑似线谱进行评价分析，去除伪线谱，进而精确地提取船舶辐射噪声线谱。
6.本发明的技术方案为：
7.所述一种基于统计量分析的船舶辐射噪声线谱精确提取方法，包括以下步骤：
8.步骤1：获取舰船发出的辐射噪声信号x，对舰船发出的辐射噪声信号x进行均值归一化处理：
9.对均值归一化处理后的信号进行短时傅里叶变换(short time fourier transform,stft)，得到信号lofar图对应的数据矩阵：x维度为m
×
n；x的时
间跨度为t(单位：s)，频率跨度为f(单位：hz)；
10.对lofar图的数据矩阵x进行最大最小归一化处理：
[0011][0012]
步骤2：对进行时间上的累积，得到不同频率点在t时间段内能量之和，记作向量p，向量p中第j个元素为：将不同频率点对应的能量向量p进行最大最小归一化处理：
[0013]
步骤3：取向量的极小值点p
min
及其对应的频率点位置f
min
，利用f
min
及p
min
拟合连续谱趋势曲线，得到函数φ：p
min
＝φ(f
min
)；利用得到的函数φ，代入向量对应的频率位置f，得到连续谱各频率点对应能量值，记作置f，得到连续谱各频率点对应能量值，记作
[0014]
步骤4：将向量的各个元素减去中所对应的元素，得到去除连续谱趋势后的各频率点所对应的能量值大小，记作p
d
：将去除连续谱趋势后能量向量p
d
进行最大最小归一化处理：
[0015]
步骤5：将中元素从大到小排列，形成新的向量p'
d
，取p'
d
中前1/2的元素取平均值的2倍作为提取线谱阈值s：int(n/2)为对n/2取整；
[0016]
步骤6：将中小于阈值s的元素置零：中小于阈值s的元素置零：中未置零的元素位置所对应的频率即为所提出线谱所对应的频率f
p
：f
p
＝(f
p1
,f
p2
,
…
,f
pn
)，n为所提取出线谱数量；
[0017]
步骤7：求lofar图所对应数据矩阵每一列的最大值，记作向量a：a＝(a1,a2,
…
,a
n
)，其中
[0018]
求lofar图所对应数据矩阵每一列的最小值，记作向量b：b＝(b1,b2,
…
,b
n
)，其中
[0019]
求lofar图所对应数据矩阵每一列的平均值，记作向量c：c＝(c1,c2,
…
,c
n
)，其中
[0020]
求lofar图所对应数据矩阵每一列的中位数，记作向量d：d＝(d1,d2,
…
,d
n
)，其中x'为中每一列元素从大到小排列后矩阵；
[0021]
求lofar图所对应数据矩阵每一列的极差，记作向量g：g＝a
‑
b；
[0022]
求lofar图所对应数据矩阵每一列的方差，记作向量h：h＝(h1,h2,
…
,h
n
)，其中
[0023]
步骤8：将向量a在线谱频率f
p
对应列的值取出，记作对应列的值取出，记作同理将向量b，c，d，g，h在线谱频率f
p
对应列的值取出，记作
[0024]
步骤9：将中最小的k＝int(n/5)个元素依照从小到大依次赋值为(k,k
‑
1,
…
,1)，其余元素置0；同理将中最小的k＝int(n/5)个元素依照从小到大依次赋值为(k,k
‑
1,
…
,1)，其余元素置0；将最大的k＝int(n/5)个元素依照从大到小依次赋值为(k,k
‑
1,
…
,1)；
[0025]
步骤10：将新生成向量进行加权处理并相加得到新向量，记作l，优选并设计阈值，优选阈值q＝(k*4*2 k*2)*0.6，当l(i)≥q时，则该处为伪线谱；将f
p
中伪线谱频率去除，得到实际存在线谱所对应频率f
p
'。
[0026]
有益效果
[0027]
本发明利用最大值、最小值、平均值、中位数、极差、方差作为评价指标，提出一种基于统计量分析的船舶辐射噪声线谱精确提取方法，先利用去除连续谱趋势方法提取出可能包含伪线谱的疑似线谱，再求取这些线谱在该时间段的能量统计量，并对这些统计量进行排序赋值并加权，能有效去除疑似线谱之中的伪线谱，精确提取出lofar图中的线谱成分。
[0028]
本发明的附加方面和优点将在下面的描述中部分给出，部分将从下面的描述中变得明显，或通过本发明的实践了解到。
附图说明
[0029]
本发明的上述和/或附加的方面和优点从结合下面附图对实施例的描述中将变得明显和容易理解，其中：
[0030]
图1：本发明方法的流程图。
[0031]
图2：用于本发明方法的仿真信号lofar图。
[0032]
图3：用于本发明方法的仿真信号去除伪线谱后线谱图。
[0033]
图4：用于本发明方法的实测信号lofar图。
[0034]
图5：用于本发明方法的实测信号去除伪线谱后线谱图。
具体实施方式
[0035]
下面详细描述本发明的实施例，所述实施例的示例在附图中示出，其中自始至终相同或类似的标号表示相同或类似的元件或具有相同或类似功能的元件。下面通过参考附图描述的实施例是示例性的，旨在用于解释本发明，而不能理解为对本发明的限制。
[0036]
本实施例针对辐射舰船噪声的线谱提取问题，提出了一种利用去除连续谱趋势提
取疑似线谱方法，并针对提取出的疑似线谱利用统计量中最大值、最小值、均值、中位数、极差、方差进行评价分析，以达到去除其中伪线谱的目的。本发明能更精确的提取出lofar图中疑似线谱频率位置，并针对提取出的疑似线谱的能量统计量进行评价分析，进而更为精准提取出舰船辐射噪声线谱。
[0037]
本实施例所用信号x为利用matlab仿真模拟信号，模拟时长为84s。
[0038]
(1)获取辐射噪声信号x＝x1 x2 x3 x4 x5 x6 x7 x8，
[0039]
其中信号x1＝2*sin(2*π*f1*t)，x2＝10*sin(2*π*f2*t)，x3＝5*cos(2*π*f3*t)，x4＝7*sin(2*π*f4*t)，x5＝sin(2*π*f5*t)，x6＝3*sin(2*π*f6*t)，x7＝4*cos(2*π*f7*t)，x8＝cos(2*π*f8*t)；f1＝100hz，f2＝200hz，f3＝30hz，f4＝400hz，f5＝45hz，f6＝254hz，f7＝85hz，f8＝213hz。对辐射噪声信号x进行加噪处理，是信噪比snr＝
‑
30db。
[0040]
(2)训练样本的预处理。将x依次进行均值归一化处理、短时傅里叶变换和最大最小归一化处理后得lofar图的数据矩阵处理后得lofar图的数据矩阵的维度为159
×
500，时间跨度为t＝84s，频率跨度为f＝500hz。
[0041]
(3)对进行时间上的累积，得到不同频率点在84s内能量之和为p，向量p中第j个元素为：并将p进行最大最小归一化处理，得到
[0042]
(4)取向量的极小值点p
min
及其对应的频率点位置f
min
，利用f
min
及p
min
拟合连续谱趋势曲线，得到函数φ：p
min
＝φ(f
min
)。
[0043]
(5)利用得到的函数φ，代入向量对应频率位置f，得到连续谱各频率点对应能量值，记作值，记作
[0044]
(6)将向量的各个元素减去中所对应的元素，得到去除连续谱趋势后的各频率点所对应的能量值大小，记作p
d
：并对p
d
进行最大最小归一化处理，得到
[0045]
(7)将中元素从大到小排列，形成新的向量p'
d
，取p'
d
中前1/2的元素取平均值的2倍作为提取线谱阈值s：
[0046]
(8)将中小于阈值s的元素置0：中小于阈值s的元素置0：中未置0的元素位置所对应的频率即为所提出可能包含伪线谱的疑似线谱所对应的频率f
p
：f
p
＝(30,45,85,100,200,213,254,400,461)hz，提取出线谱数量n＝9。
[0047]
(9)求lofar图所对应数据矩阵每一列的最大值a：a＝(a1,a2,
…
,a
500
)，其中a
j
＝max(x
1j
,x
2j
,
…
,x
159j
),j＝1,2,
…
,500。
[0048]
(10)求lofar图所对应数据矩阵每一列的最小值b：b＝(b1,b2,
…
,b
500
)，其中b
j
＝min(x
1j
,x
2j
,l,x
159j
),j＝1,2,
…
,n。
[0049]
(11)求lofar图所对应数据矩阵每一列的平均值c：c＝(c1,c2,
…
,c
500
)，其中
[0050]
(12)，求lofar图所对应数据矩阵每一列的中位数，记作向量d：d＝(d1,d2,
…
,d
500
)，其中d
j
＝x'
79j
，x'为中每一列元素从小到大排列后矩阵。
[0051]
(13)求lofar图所对应数据矩阵每一列的极差，记作向量g：g＝a
‑
b。
[0052]
(14)求lofar图所对应数据矩阵每一列的方差，记作向量h：h＝(h1,h2,
…
,h
500
)，其中
[0053]
(15)将向量a中提取出的线谱频率f
p
对应列的值取出，记作对应列的值取出，记作同理将向量b，c，d，g，h中提取出的线谱频率f
p
对应列的值取出，记作
[0054]
(16)将中最小的2个元素依照从小到大依次赋值(2,1)，其余元素置0；同理将中最小的2个元素依照从小到大依次赋值(2,1)，其余元素置0；将最大的2个元素依照从大到小依次赋值(2,1)，其余元素置0。
[0055]
(17)将新生成向量进行加权处理并相加得到新向量，记作l：设计阈值q＝12，当l(i)≥12时，则该处为伪线谱。
[0056]
(18)将f
p
中伪线谱频率461hz去除，得到实际存在线谱所对应频率f
p
'＝(30,45,85,100,200,213,254,400)hz。
[0057]
尽管上面已经示出和描述了本发明的实施例，可以理解的是，上述实施例是示例性的，不能理解为对本发明的限制，本领域的普通技术人员在不脱离本发明的原理和宗旨的情况下在本发明的范围内可以对上述实施例进行变化、修改、替换和变型。