CN106895905B - 一种舰船辐射噪声检测方法 - Google Patents

Abstract

本发明提供了一种舰船辐射噪声检测方法，涉及信息信号处理领域，提出一种基于时域混沌成分差异的舰船辐射噪声检测方法，采用了Duffing振子检测系统及统计复杂度参数表征，得到了缺少先验信息的水下低信噪比舰船目标检测的技术效果，并实现了嵌入式表征，海上目标监测及嵌入式无人海上预警等具有重要参考价值。

Description

一种舰船辐射噪声检测方法

技术领域

本发明涉及信息信号处理领域，尤其是一种噪声监测方法。

背景技术

进入二十一世纪，世界各国均加强了对海洋领土的监管和控制，领海中的资源将在不远的未来对本国的发展起到至关重要的作用。要及时有效地检测和识别别国舰船是至关重要的一环。同时，准确地检测到远距离下的舰船也是非常重要的研究方向，因而，海洋环境下远距离舰船的检测和识别具有巨大的研究价值和意义。

对较远距离舰船的检测主要是通过检测舰船辐射噪声来实现，舰船辐射噪声具有特殊的线谱和连续谱组成的频谱，通过提取线谱和连续谱的特征可实现舰船辐射噪声的检测。基本思路之一为检测舰船辐射噪声中的线谱分量频率峰值来实现目标的检测。其中常用的自相关检测方法、快速傅里叶变换方法、自适应线谱增强方法需要对线谱频率峰值进行搜索，得到线谱先验信息，无此先验信息时，结果会受到很大的影响，这种搜索过程很大程度上依赖于经验。同时，由于舰船声隐身技术的高速发展，机械设备主动隔振、主动约束阻尼、声学智能结构等技术的陆续应用，舰船辐射噪声的线谱得到了很好的控制，线谱的幅值变得很小、能量大幅度降低、甚至数量可控。

连续谱是舰船辐射噪声频谱的另一重要特征和组成部分，时域具有单独的谱峰。基于连续谱的检测一般基于时域谱峰的能量，而基于能量的检测方法，在远距离下结果会受到较大影响。连续谱分量声压级也与海洋环境噪声声压级更加接近。同时，在声呐监测过程中，远距离下的目标舰船可得先验知识非常有限，线谱频率和数目以及相应的变化等无法确定。海洋环境噪声的复杂性也影响着舰船辐射噪声连续谱分量。使得利用常规方法很难实现远距离下舰船目标的有效检测。

发明内容

为了克服现有技术的不足，针对常规方法依赖先验信息、在远距离低信噪比下效果变差、海洋环境噪声假设不准确等问题，本发明提出一种基于时域混沌成分差异的舰船辐射噪声检测方法。

本发明解决其技术问题所采用的技术方案的具体步骤如下：

第一步：利用声呐采集海洋中的声信号，记为g(t)，即为输入信号；

第二步：构造信号检测输入系统

利用Duffing振子检测系统，系统公式为：

其中，式中x、y为系统输出，分别x、y的导数，k为阻尼系数，-αx+βx³为非线性恢复力，α、β为非线性恢复力系数，rcos(ωt)为内策动力，ω为内策动力频率，r为内策动力幅值，g(t)为输入信号，参数设置为k＝0.5,α＝1,β＝1,ω＝1，利用Melnikov方法取r＝0.8，使系统处于混沌临界状态；

Duffing振子系统初始状态为

第三步：将第一步中的信号g(t)输入到第二步中的式(1)中，利用四阶自适应步长龙格库塔方法对公式(1)进行求解，初值定为(1,1)，得出系统方程解(x,y)，其中x,y均为向量，以方程解中的x为横轴，y为纵轴，得到输出数据的相空间图形；

输入信号同时含有单频信号、混沌特性信号和高斯噪声信号，即

其中s₁(t)＝Acos(ωt)，A为输入信号幅值，s₂(t)为一种或多种混沌特性信号的和，s₃(t)为高斯白噪声信号或色噪声信号，其中ξ,ζ,为0或1，且ξ,ζ,不同时为零；

第四步，搜索公式(1)的解(x,y)，分别得到向量x的最大值X_max和最小值X_min以及向量y的最大值Y_max和最小值Y_min，以X_max到X_min的距离为宽度，Y_max到Y_min的距离为高度组成一个长方形区域，将此长方形区域均分为2行n列的网格；

第五步：统计每个网格中的数据量k_i，利用

计算每个网格出现数据的概率p_i，得到2*n个格子的分布列，其中C为系统相空间中的总数据量；

第六步：利用Bandt-Pompe算法及统计复杂度计算方法得到接收信号的统计复杂度，根据统计复杂度值大小进行目标检测，统计复杂度的计算方法如下：

Shannon熵表示概率分布为P＝{pi,i＝1,…,N}的物理过程的不确定程度，表述为

Shannon熵的最大值为概率分布为均匀分布P_e＝{1/N,…,1/N}时的取值，由此，标准Shannon熵为

其中S_max＝S[P_e]＝lnN，P_e＝{1/N,…,1/N}表示均匀分布，0≤H_S[P]≤1；

失衡度K[P]度量系统任一状态T时的概率分布P与均匀分布P_e之间的距离D_S表述为：

K[P]＝K₀·D_S[P,P_e] (7)

其中K₀是归一化常数，则0≤K≤1，D_S选用Jensen-Shannon散度J_S进行刻画，即对概率空间中任意两个分布P₁和P₂，表述为：

J_S[P₁，P₂]＝{S[(P₁+P₂)/2]-S[P₁]/2-S[P₂]/2} (8)

那么，失衡度表述为

K_J[P]＝K₀·J_S[P,P_e] (9)

其中，归一化常数K₀为J_S[P,P_e]取最大值时的倒数，完全有序状态和均匀分布之间的距离为J_S[P,P_e]取值的最大值；

则由如式(6)所示的标准Shannon熵和式(9)所示的失衡度，可得统计复杂度为：

C_JS[P]＝K_J[P]·H_S[P] (9)

定义阈值为0.05，若统计复杂度C_JS[P]大于0.05则判定有目标，小于0.05则判定无目标。

本发明由于采用了Duffing振子检测系统及统计复杂度参数表征，得到了缺少先验信息的水下低信噪比舰船目标检测的技术效果，并实现了嵌入式表征，海上目标监测及嵌入式无人海上预警等具有重要参考价值。

附图说明

图1是本发明同宿轨道状态相空间图。

图2是本发明分叉状态相空间图。

图3是本发明混沌状态相空间图。

图4是本发明大尺度周期状态相空间图。

图5是本发明输入Lorenz信号系统输出相空间3。

图6是本发明输入Lorenz信号系统输出相空间4。

图7是本发明加入高斯白噪声特殊相空间。

图8是本发明是参数r＝0.8时相空间。

图9是本发明中1型舰船辐射噪声在输入信号进行连续三点平移后的输出相空间形态。

图10是本发明2型舰船辐射噪声在输入信号进行连续三点平移后的输出相空间形态。

图11是本发明3型舰船辐射噪声在输入信号进行连续三点平移后的输出相空间形态。

图12是本发明4型舰船辐射噪声在输入信号进行连续三点平移后的输出相空间形态。

图13是本发明海洋环境噪声时域波形。

图14是本发明海洋环境噪声输入系统相空间1。

图15是本发明海洋环境噪声输入系统相空间2。

图16是本发明检测效果图。

图中x1和x2分别代表两列数据。

具体实施方式

下面结合附图和实施例对本发明进一步说明。

第一步：利用声呐采集海洋中的声信号，记为g(t)，即为输入信号；

第二步：构造信号检测输入系统

利用Duffing振子检测系统，系统公式为：

其中，式中x、y为系统输出，分别x、y的导数，k为阻尼系数，-αx+βx³为非线性恢复力，α、β为非线性恢复力系数，rcos(ωt)为内策动力，ω为内策动力频率，r为内策动力幅值，g(t)为输入信号，参数设置为k＝0.5,α＝1,β＝1,ω＝1，利用Melnikov方法取r＝0.8，使系统处于混沌临界状态；

Duffing振子系统初始状态为

式(2)所示系统在固定k时，随着r从0逐渐增加，系统相空间将出现不同形态的变化，如图1～图4所示，其中以混沌状态到大尺度周期状态相空间变化最为明显，将待测信号以加项形式输入系统，由于系统对单频信号敏感和对随机信号免疫，则根据相空间变化可以完成单频信号的检测；

第三步：将第一步中的信号g(t)输入到第二步中的式(1)中，利用四阶自适应步长龙格库塔方法对公式(1)进行求解，初值定为(1,1)，得出系统方程解(x,y)，其中x,y均为向量，以方程解中的x为横轴，y为纵轴，得到输出数据的相空间图形；

输入信号为混沌信号时，Duffing振子输出相空间为如图5和图6所示形态；

输入信号同时含有单频信号、混沌特性信号和高斯噪声信号，即

其中s₁(t)＝Acos(ωt)，A为输入信号幅值，s₂(t)为一种或多种混沌特性信号的和，s₃(t)为高斯白噪声信号或色噪声信号，其中ξ,ζ,为0或1，且ξ,ζ,不同时为零；输入如式(3)所示信号时，输出相空间依然如图3、图4、图5和图6所示形态，当输入信号为舰船辐射噪声信号时，输出相空间出现如图7所示特殊相空间形态，且这种特殊形态相空间会随着输入信号的平移连续出现；

第四步：搜索公式(1)的解(x,y)，分别得到向量x的最大值X_max和最小值X_min以及向量y的最大值Y_max和最小值Y_min，以X_max到X_min的距离为宽度，Y_max到Y_min的距离为高度组成一个长方形区域，将此长方形区域均分为2行n列的网格，本实施例中取n＝5；

第五步：统计每个网格中的数据量k_i，利用

计算每个网格出现数据的概率p_i，得到2*n个格子的分布列，其中C为系统相空间中的总数据量；

第六步：利用Bandt-Pompe算法及统计复杂度计算方法得到接收信号的统计复杂度，根据统计复杂度值大小进行目标检测，统计复杂度的计算方法如下：

统计复杂度是在Shannon熵的基础上进行定义的，Shannon熵表示概率分布为P＝{p_i,i＝1,…,N}的物理过程的不确定程度，表述为

Shannon熵的最大值为概率分布为均匀分布P_e＝{1/N,…,1/N}时的取值，由此，标准Shannon熵为

其中S_max＝S[P_e]＝lnN，P_e＝{1/N,…,1/N}表示均匀分布，0≤H_S[P]≤1；

失衡度K[P]度量系统任一状态T时的概率分布P与均匀分布P_e之间的距离D_S表述为：

K[P]＝K₀·D_S[P,P_e] (7)

其中K₀是归一化常数，则0≤K≤1，D_S选用Jensen-Shannon散度J_S进行刻画，即对概率空间中任意两个分布P₁和P₂，表述为：

J_s[P₁，P₂]＝{S[(P₁+P₂)/2]-S[P₁]/2-S[P₂]/2} (8)

那么，失衡度表述为

K_J[P]＝K₀·J_S[P,P_e] (9)

其中，归一化常数K₀为J_S[P,P_e]取最大值时的倒数，完全有序状态和均匀分布之间的距离为J_S[P,P_e]取值的最大值；

则由如式(6)所示的标准Shannon熵和式(9)所示的失衡度，可得如下述统计复杂度：

C_JS[P]＝K_J[P]·H_S[P] (9)

定义阈值为0.05，若统计复杂度C_JS[P]大于0.05则判定有目标，小于0.05则判定无目标。

利用仿真及海上实测数据，进行最低检测信噪比的换算和对比，并分析实时性：

根据以下信噪比的计算公式，可得利用Duffing振子检测舰船辐射噪声的最低检测信噪比：

其中P_signal为信号的平均功率，P_noise为噪声的平均功率，最低检测信噪比为-9.1329，由于本发明并非以统计理论为基础得出的，故而无法得出在不同虚惊概率下的检测概率曲线。为了得到检测概率，根据本发明所述方法做1000次蒙特卡洛实验，结果并未出现检测失败的情况，而若输入信号信噪比低于最低检测信噪比，检测概率则是0％，此方法失效。值得说明的是，声纳收集到的舰船辐射噪声本身就含有海洋环境噪声，故而实际的最低检测信噪比计算得到的还要低。

本发明主要涉及的计算包括归一化、四阶龙格-库塔方法解Duffing振子方程、Bandt-Pompe算法和统计复杂度的计算。其中，归一化和统计复杂度涉及简单的乘法和除法，Bandt-Pompe算法主要为分段和除法，而龙格-库塔方法比较复杂，占据了主要的运算时间，将实时性以运算时间衡量，得到表1：

表1 Duffing振子检测方法实时性判断

计算机配置不同、数据不同等会得到不同的运算时间，因而仅利用以上表格观察耗时趋势。首先，由表格可以看出总体耗时均在十秒以下，相对于海上舰船的运动速度，可以达到实时性要求；其次，输入信号点数对耗时没有很大影响，主要取决于龙格-库塔方法步长的选择，考虑到相空间的完整性，选取自适应步长方法。

由以上分析可知，本发明可实现真实海洋环境下，缺少先验信息时，对低信噪比下舰船目标的有效检测，对海上目标监测及嵌入式无人海上预警等具有重要参考价值。

利用海上实测数据，描述目标辐射噪声检测的过程。

对于式(1)所示系统，采用自适应步长四阶龙格-库塔方法进行数值计算，取k＝0.5,α＝1,β＝1,ω＝1，初值为(1,1)，经测量，系统内策动力幅值r＝0.863～0.864之间时，系统处于混沌临界状态，取r＝0.8，调整系统至混沌到大尺度周期的临界状态，如图8所示。

选用四种舰船，分别称为1型～4型舰船，分别以15节航速航行时实测到的舰船辐射噪声数据为对象进行分析，声呐收集信号时的采样频率为24000Hz，将舰船辐射噪声信号归一化，分别记为s_J1(t)，s_J2(t)，s_J3(t)，s_J4(t)，分别对四种信号取s1_J1(t)＝[s_J1(1),s_J1(2),…,s_J1(n-1),s_J1(n)]，s2_J1(t)＝[s_J1(2),s_J1(3),…,s_J1(n-2),s_J1(n-1)]，s3_J1(t)＝[s_J1(3),s_J1(4),…,s_J1(n-3),s_J1(n-2)]，即将信号进行连续三点平移，将得到的平移后信号分别输入公式(2)所示的Duffing振子检测系统，得到四种舰船辐射噪声的相空间形态：

由以上Duffing系统输出相空间可知，连续相移的舰船辐射噪声信号输入后，相空间均呈特殊形态，且这种特殊相空间形态有两种，可以由其排列的不同进行目标的识别。

四级海况下，相对平稳平台上以102.4kHz采样率采集某海域环境噪声，如图13所示为未经归一化的海洋环境噪声时域波形。

将采集到的海洋环境噪声归一化，输入式(1)所示系统，系统输出相空间如图14和图15所示。

取此信号中不同时段部分分别输入系统，仍然得到如图14和图15所示相空间形态，只是左右势阱高度略有差别；收集其他海域3～5级海况环境噪声数据，做如上处理，得到结果相同。故而，连续相移的海洋环境噪声背景下的舰船辐射噪声输入式(1)所示系统后，得到如图9～图12的相空间形态，而海洋环境背景噪声下的其它信号则不具备这种输出相空间，利用这种复杂相空间形态可实现真实海洋环境噪声背景下舰船辐射噪声的检测。

首先计算图3所示的混沌临界状态和图4所示大尺度周期状态相空间形态的统计复杂度，其次计算输出相空间对应的分布，由于相空间会出现图4所示的中空形态，故将y轴分为上下两层，即将相空间分割为2×n的小块，取n＝5，取有限个数分割的前提下，均为离散型随机变量，离散型随机变量的分布以分布列的形式表示。分割块数取的越大，对应的分布越精确，分割块数为无穷时，即得到概率密度曲线。故而可得图3和图4对应的分布列，分别如表和表3所示，其中u为块数标号，p_i为每块对应的概率值，保留小数点后三位。

表2混沌临界状态相空间分布列

u	1	2	3	4	5	6	7	8	9	10
											p<sub>i</sub>	0.104	0.137	0.060	0.101	0.077	0.078	0.118	0.072	0.144	0.109

表3大尺度周期状态相空间分布列

u	1	2	3	4	5	6	7	8	9	10
											p<sub>i</sub>	0.119	0.166	0.055	0.092	0.069	0.067	0.091	0.055	0.166	0.120

根据以上分布列可得混沌临界状态统计复杂度为0.0169，大尺度周期状态统计复杂度为0.0349。

同样可得出图9～图12所示两种特殊形态相空间的分布列：

表4特殊相空间1分布列

u	1	2	3	4	5	6	7	8	9	10
											p<sub>i</sub>	0.001	0.001	0.090	0.034	0.195	0.086	0.343	0.068	0.131	0.051

表5特殊相空间2分布列

u	1	2	3	4	5	6	7	8	9	10
											p<sub>i</sub>	0.075	0.225	0.135	0.555	0.001	0.003	0.001	0.003	0.001	0.001

根据以上分布列，可得两种特殊相空间分布的统计复杂度分别为0.1907和0.2800，由以上计算可知，特殊相空间统计复杂度值与正常相空间值差别较大，将阈值定为0.05，将输入信号为g(t)每隔i个时刻，进行连续三次平移，分别输入Duffing振子系统，i＝2时得到以下检测效果，如图16所示。

图16中，虚线为阈值大小，实线为检测值即统计复杂度值，图中1～2时刻系统处于临界状态为检测准备时刻，3～5时刻为第一个检测时段，第4时刻有一单频干扰信号进入系统，使得系统进入了大尺度周期状态，此处假设单频信号已与内策动力同步的理想情况；8～10时刻为第二个检测时段，目标信号为舰船辐射噪声；13～15时刻为第三个检测时段，目标信号为与第二时段不同的舰船辐射噪声；18～20时刻为第四个检测时段，目标信号为具有一定复杂度的干扰信号。

Claims (1)

1.一种舰船辐射噪声检测方法，其特征在于包括下述步骤：

第一步：利用声呐采集海洋中的声信号，记为g(t)，即为输入信号；

第二步：构造信号检测输入系统

利用Duffing振子检测系统，系统公式为：

其中，式中x、y为系统输出，分别为 x、y的导数，k为阻尼系数，-αx+βx³为非线性恢复力，α、β为非线性恢复力系数，rcos(ωt)为内策动力，ω为内策动力频率，r为内策动力幅值，g(t)为输入信号，参数设置为k＝0.5,α＝1,β＝1,ω＝1，利用Melnikov方法取r＝0.8，使系统处于混沌临界状态；

Duffing振子系统初始状态为

第三步：将第一步中的信号g(t)输入到第二步中的式(1)中，利用四阶自适应步长龙格库塔方法对公式(1)进行求解，初值定为(1,1)，得出系统方程解(x,y)，其中x,y均为向量，以方程解中的x为横轴，y为纵轴，得到输出数据的相空间图形；

输入信号同时含有单频信号、混沌特性信号和高斯噪声信号，即

其中s₁(t)＝Acos(ωt)，A为输入信号幅值，s₂(t)为一种或多种混沌特性信号的和，s₃(t)为高斯白噪声信号或色噪声信号，其中ξ,ζ,为0或1，且ξ,ζ,不同时为零；

第四步，搜索公式(1)的解(x,y)，分别得到向量x的最大值X_max和最小值X_min以及向量y的最大值Y_max和最小值Y_min，以X_max到X_min的距离为宽度，Y_max到Y_min的距离为高度组成一个长方形区域，将此长方形区域均分为2行n列的网格；

第五步：统计每个网格中的数据量k_i，利用

计算每个网格出现数据的概率p_i，得到2*n个格子的分布列，其中C为系统相空间中的总数据量；

第六步：利用Bandt-Pompe算法及统计复杂度计算方法得到接收信号的统计复杂度，根据统计复杂度值大小进行目标检测，统计复杂度的计算方法如下：

Shannon熵表示概率分布为P＝{p_i,i＝1,…,N}的物理过程的不确定程度，表述为

Shannon熵的最大值为概率分布为均匀分布P_e＝{1/N,…,1/N}时的取值，由此，标准Shannon熵为

其中S_max＝S[P_e]＝lnN，P_e＝{1/N,…,1/N}表示均匀分布，0≤H_S[P]≤1；

失衡度K[P]度量系统任一状态T时的概率分布P与均匀分布P_e之间的距离D_S表述为：

K[P]＝K₀·D_S[P,P_e] (7)

其中K₀是归一化常数，则0≤K[P]≤1，D_S选用Jensen-Shannon散度J_S进行刻画，即对概率空间中任意两个分布P₁和P₂，表述为：

J_S[P₁，P₂]＝{S[(P₁+P₂)/2]-S[P₁]/2-S[P₂]/2} (8)

那么，失衡度表述为

K_J[P]＝K₀·J_S[P,P_e] (9)

其中，归一化常数K₀为J_S[P,P_e]取最大值时的倒数，完全有序状态和均匀分布之间的距离为J_S[P,P_e]取值的最大值；

则由如式(6)所示的标准Shannon熵和式(9)所示的失衡度，可得统计复杂度为：

定义阈值为0.05，若统计复杂度大于0.05则判定有目标，小于0.05则判定无目标。