本发明属于声呐信号处理领域，特别涉及水下声信号检测及参数估计分析处理，提出了一种基于恒虚警检测的水下声信号时频端点参数估计方法，为水下目标检测、水声试验数据分析、目标识别分类等提供相关参数。
背景技术
随着水声领域的迅速发展，水声信号检测及分析问题成为该领域专家日益关切的问题之一。海洋中存在种类繁多的信号，例如水声通信信号、声呐探测信号、导航信号、舰船辐射噪声、海洋生物叫声、海洋地质活动产生的声波等。如果能够将水声信号的时频端点等特征提取出来，则能对后续的水声试验数据分析、水下信号分类、水下目标识别提供有效信息，从而深化对水声信号的认识，提高信号处理精度。
目前，有关水下信号时间和频率端点参数估计的公开方法极难找到，可供参考的是语音信号处理中的时间端点参数估计。在语音信号处理中，需要通过对时间端点参数的估计来截取说话区，进一步对数据进行处理，例如自然语义分析等。常用的语音信号时间端点参数估计方法大致分为两类：
第一类是基于阈值的方法，基本思路是首先提出语音信号每一帧的特征值，将该特征值与预先设定的特征值进行比较，从而判断该帧是否包含信号。典型的特征值包括：时域能量、过零率、线性预测参数等。其中，低信噪比情况下，能量、过零率、线性预测参数不够可靠，且过零率对加性噪声很敏感。该方法的优点在于简单快速，但准确率受噪声影响较大，且只能估计时间端点参数。
第二类是基于模式匹配的方法。典型算法包括：k最近邻分类法、支持向量机法、神经网络法等。该方法随着数据逐渐完善，准确率不断提升，但需要训练数据，这就需要已知待检测信号的先验信息，且训练时间长、算法复杂。
目前，常用的信号端点参数估计方法中存在的问题如下：
(1)常用的基于阈值的语音信号时间端点参数估计方法，估计结果受信噪比影响较大。
(2)常用的基于模式匹配的语音信号时间端点参数估计方法，需要已知待检测信号的先验信息，不适用于未知条件下的端点参数估计系统。
(3)目前常用的端点参数估计方法主要是用来估计时间起止点。
对信号的时间端点和频率端点同时进行有效估计是实际工程中迫切需要的一种技术，本发明通过对信号的时间端点估计，可以对水声信号进行时域截取，提高信噪比，为下一步信号分析提供更好的素材；实际工程中，因为不同类型的水声信号处在不同的频带范围内，因此通过对信号的频率端点进行估计，可以为信号粗分类提供依据，亦可通过对信号进行频域截取为下一步详细处理(例如使用机器学习方法进行信号分类等)提供有用的素材。
技术实现要素：
要解决的技术问题
为了解决水下信号时间和频率端点参数估计问题，本发明提出了一种基于恒虚警检测的水下信号时频端点参数估计方法，用于对水声信号进行时频端点参数估计，提取信号的时间起止点、频率起止点等相关信息特征。不需要先验信息，适用于未知条件下的端点参数估计系统，受信噪比影响较小且检测虚警率恒定，能够对水声信号的时间端点和频率端点同时进行估计。
技术方案
一种基于恒虚警检测的水下信号时频端点参数估计方法，其特征在于步骤如下：
步骤1：读取数据，对数据进行分帧处理：
x'(n)＝x(t1*fs+n)(0≤n＜t*fs)
其中，x(n)为接收数据，t1为时域截取的开始时刻，fs为数据采样频率，t为截取长度；n为分帧后的点数；
步骤2：对分帧后的数据进行分频带处理：频带划分采用非等间隔的划分方法，低频段频带划分相对高频段频带划分较窄；对于在高频段和低频段中间部分的信号，则采用非等间隔的频带重叠划分方法进行处理；
根据划分后的频带范围，得到每个频带对应的滤波器下限频率fl和上限频率fh；使用傅里叶变换算法将x'(n)变换到频域，得到频谱x'(f)；根据滤波器下限频率fl和上限频率fh，设计频域滤波器，对信号进行频域滤波，得到进行分频段处理后的信号x*(n)；
步骤3：对经过分频带处理的数据进行时频分析：
对分频带后的数据进行时频分析，采用短时傅里叶变换的方法，得到信号的声谱图：
式中，x*(n)为经过分频段后的时域数据，n的值等于ts×fs，ts为移动的时间段；w(n)为窗函数，其数学表示为：
sx(m,r)中的m代表对信号进行截取的时刻所对应的点数，等于a×τ×fs，其中a为大于等于0的整数且满足a×τ＜t，τ为每次进行帧移的长度，r代表短时傅里叶变换后的频率；
步骤4：对时频分析结果进行恒虚警检测：
提取sx(m,r)中的第i行、第j列元素作为检测单元的数据tij，按照一维恒虚警检测窗结构提取检测单元两侧的参考单元rij；
对rij中的结果进行由大到小排序得到
式中，m为参考单元个数，针对检测单元tij的检测阈值t(i,j)的表达式如下：
t(i,j)＝α×ttemp(i,j)
其中，α的值和虚警率pfa之间关系如下：
pfa＝1-iα(n,m)
式中，n为检测单元个数，m为参考单元个数，iα(n,m)为不完全归一化beta函数，其定义如下：
式中，γ(n)为伽马函数；
对所有检测单元tij平方求和：
γij＝(tij)2
进行如下比较
如果γij大于t(i,j)，说明该检测单元中有信号，xcfar(i,j)＝1；否则没有，xcfar(i,j)＝0；
经过恒虚警检测后，得到尺寸和sx(m,r)矩阵一致的矩阵xcfar，其每个元素只等于1或是0，其中1代表该单元有信号，0代表该单元没有信号；
步骤5：计算恒虚警检测频率边际谱并判决：
计算恒虚警检测结果xcfar的频率边际谱xpf(f)，即对xcfar进行时间积分，其表达式为：
对凹凸不平的频率边际谱进行滤波平滑，将其与阈值tf＝a×k×pfa进行比较，其中，pfa为虚警概率，k为时间积分长度，a为调整系数；
如果在一段频率范围内只有一个频点上的能量超过了阈值，则认为该频点上存在线谱信号，该频点认为是线谱信号的中心频率f0；如果在一段频率范围内的频点都超过了阈值，则认为该频率范围内存在块状信号，该频率范围的上下限认为是块状信号的频率上下限f′h和f′l；
步骤6：对恒虚警检测结果进行频域滤波，计算时间边际谱并判决：
根据步骤5中得到的上下限频率f′h和f′l，设计带通滤波器h(f)：
将恒虚警检测结果xcfar(t,f)通过该滤波器进行处理，得到x'cfar(t,f)，其表达式为：
x'cfar(t,f)＝xcfar(t,f)*h(f)
计算x'cfar(t,f)的时间边际谱，即对x'cfar(t,f)进行频率积分，其表达式为：
对xpt(t)进行滤波平滑，将其与阈值tt＝b×l×pfa进行比较，其中，pfa为虚警概率，l为频率积分长度，b为调整系数；
将xpt(t)中第一个超过阈值和最后一个超过阈值的点为信号的起止点tstart和tend。
步骤3中的τ＝0.5×ts。
步骤3中的ts，对低频段的信号选择1-2s；对高频段的信号选择0.01-0.05s。
步骤4中的h为6。
有益效果
本发明提出的一种基于恒虚警检测的水下信号时频端点参数估计方法，与现有技术相比具有的有益效果是：
(1)本发明提出了对水下信号同时进行时间端点和频率端点参数估计的方法：在对信号进行分帧、分频段处理的基础上，首先对信号进行时频变换得到声谱图，在此基础上对信号进行恒虚警检测，对检测结果分别计算频率边际谱和时间边际谱并判决，从而得到信号的频率和时间端点估计结果。
(2)从技术角度来看，本发明能够同时估计信号的时间和频率端点，较传统方法有更高检测概率，能够提取声谱图中难以发觉的信息，从而获得对应的目标信息；检测过程不需要先验信息，也不需要大规模数据训练，适用于未知条件下的端点参数估计系统。
(3)从经济效益角度看，一方面，借助计算机、声学传感器等硬件设备，本发明可以构建一个水声信号现场分析仪，为海洋科考、海军作战提供有效信息；另一方面，本发明与后续的信号分析、目标分类识别、机器学习等相关水下信号处理方法相结合，能够为水声信号分析提供一套完整的软件分析工具。
附图说明
图1本发明流程图
图2短时傅里叶变换说明图
图3一维恒虚警检测窗示意图
图4本算法中恒虚警检测流程图
图5常见水声信号经过本发明中恒虚警检测的结果示意图：a)线谱信号、b)lfm信号、c)hfm信号、d)复杂信号；
图6常见水声信号的频域边际谱示意图：a)线谱信号、b)lfm信号、c)hfm信号、d)复杂信号；
图7常见水声信号的时域边际谱示意图：a)线谱信号、b)lfm信号、c)hfm信号、d)复杂信号。
具体实施方式
现结合实施例、附图对本发明作进一步描述：
本发明核心技术如下：
1.非等间隔的频率带宽划分
对于种类繁多的水声信号、信号成分复杂的水声数据，通常低频段信号的频率随时间变化较慢，例如线谱信号，而高频段信号的频率随时间变化较快，例如调制信号、瞬态信号、伪随机信号等。根据这个现象，本发明中，需要依据水声信号、水声实验数据的频率范围信息，对信号进行频带划分，采用非等间隔的频率带宽划分分析信号，以满足不同频段信号对时频分辨率的要求。在非等间隔的频率带宽划分时，低频部分的频带划分较窄，采用频域分辨率高、时域分辨率低的时频分析方式对低频信号或信号的低频成分进行表示；高频部分的频带划分较宽，采用频域分辨率低、时域分辨率高的时频分析方式表示信号或信号成分；对于在高频和低频中间部分的信号，则采用非等间隔的频带重叠划分方法进行分析，以便分析重叠部分可能包含的低频、高频信号及信号成分，满足不同信号对时频分辨率的要求。
频带划分准则依赖于水声信号、水声实验数据的先验信息。目前，水声信号主要包括水声通信信号、声呐探测信号、导航信号、舰船辐射噪声、海洋生物叫声、海洋地质活动产生的声波等，其所处的频段、信号参数为非等间隔频带划分提供了依据。对于完全未知的水声信号，则可采用粗、细相结合的频带划分方式进行分析研究，获取未知信号的信息。
2.恒虚警检测的边际谱
通常，在非等间隔带宽时频分析的恒虚警检测结果中，如果在某时频区域存在信号，那么就会有较大的能量集中在该信号所在的时段和频段。根据这个原理，本发明中，对非等间隔带宽时频分析声谱图的恒虚警检测结果，引入恒虚警检测的边际谱，通过边际谱检测估计时频端点参数。
在频率端点参数估计时，将恒虚警检测结果对时间积分，得到恒虚警检测的频域边际谱。对凹凸不平的频域边际谱滤波平滑，从频域边际谱判决信号的频率范围，得到信号的频率起止点。
在时间端点参数估计时，则首先根据频率端点参数估计结果，设置滤波器提取该频率起止点范围内的恒虚警检测结果。然后，按照频率端点参数估计的方法估计时间端点参数，即将滤波提取的恒虚警检测结果对频率积分，得到该频率起止点范围内信号的时间边际谱。滤波平滑时间边际谱，从边际谱判决信号的时间范围，得到信号的时间起止点。
3.由边际谱求端点参数的阈值
在存在信号的时频段，由恒虚警检测得到的边际谱幅度较大；在不存在信号的时频段，边际谱幅度较小；无噪声时，在不存在信号的时频段，边际谱幅度等于零。边际谱幅度的均值与恒虚警检测的虚警概率、求边际谱的频率(时间)积分长度成正比。根据这个原理，本发明中，门限值由恒虚警检测的虚警概率和求边际谱的频率(时间)轴积分的长度之积确定。在实际工程应用中，可适当的对门限值进行调整。
其具体技术方案如下：
步骤一：读取数据，对数据进行分帧处理，每帧取时长为10秒的信号。
步骤二：对分帧后的数据进行分频带处理，其中，频带划分采用非等间隔的划分方法，低频段频带划分相对较窄，高频段频带划分相对较宽。对于在高频段和低频段中间部分的信号，则采用非等间隔的频带重叠划分方法进行处理。频带划分准则依赖于水声信号、水声实验数据的先验信息。对于完全未知的水声信号，则可采用粗、细相结合的频带划分方式进行分析研究，获取未知信号的信息。
步骤三：对分频带后的数据进行时频分析，采用短时傅里叶变换的方法，得到信号的声谱图。其中，对于每个划分后的频带采用不同的时频分辨率进行时频分析，低频部分采用较高的频率分辨率，高频部分采用较高的时间分辨率。
步骤四：对时频分析后的数据进行恒虚警检测，得到对声谱图的恒虚警检测结果。需要说明的是，声谱图的尺寸和恒虚警检测结果的尺寸一致。
步骤五：计算恒虚警检测的频率边际谱并检测判决。恒虚警检测结果的横轴和纵轴分别代表时间和频率，通过对恒虚警检测结果进行时间积分，得到基于恒虚警检测结果的频率边际谱xpf(f)，滤波平滑频率边际谱曲线。再根据检测虚警率和进行时间积分的长度确定合适的阈值，与频率边际谱进行比较，当边际谱的点大于阈值时，说明该点对应频率点处存在信号，由此得到信号的频率上下限。
步骤六：对恒虚警检测结果进行频域滤波，再计算时间边际谱并判决。根据上一步得到的频率范围，设计频域滤波器，提取频带范围内的恒虚警检测结果。将该结果进行频率积分，得到时间边际谱xpt(t)，滤波平滑时间边际谱曲线，根据频率积分长度和虚警率确定合适的门限，与时间边际谱xpt(t)进行比较，超过门限的点认为该时刻有信号，由此得到信号的时间起止点。
实施例：
1)读取数据，对数据进行分帧处理。
假设x(n)为接收到的数据，其点数为fs*t，其中fs为数据采样频率，t为接收数据时长。假如对从t1时刻开始的数据进行时域截取，截取长度为t，例如t为10秒。得到的时域截取信号为
x'(n)＝x(t1*fs+n)(0≤n＜t*fs)
2)对分帧后的数据进行分频带处理
本发明中，需要依据水声信号、水声实验数据的频率范围信息，对信号进行频带划分，采用非等间隔的频率带宽划分分析信号。在非等间隔的频率带宽划分时，低频段部分的频带划分较窄(≤5khz)，高频段部分的频带划分较宽(≥15khz)，对于在高频和低频中间部分的信号，则采用非等间隔的频带重叠划分方法进行分析。取值介于低频段频带范围和高频段频带范围之间。
频带划分准则依赖于水声信号、水声实验数据的先验信息。目前，水声信号主要包括水声通信信号、声呐探测信号、导航信号、舰船辐射噪声、海洋生物叫声、海洋地质活动产生的声波等，其所处的频段、信号参数为非等间隔频带划分提供了依据。对于完全未知的水声信号，则可采用粗、细相结合的频带划分方式进行分析研究，获取未知信号的信息。
举例说明：假设采集到了一段水声数据，数据的采样频率为100khz，根据奈奎斯特采样定理，数据中所包含的信号有效频率成分不超过50khz。那么，频带范围可以如下划分：第一频带为0-5khz；第二频带为4-12khz；第三频带为10-27khz；第四频带为25-50khz。实际工程中可根据情况进行适当调整。
假设x'(n)为分帧后的数据。首先根据划分后的频带范围，得到每个频带对应的滤波器下限频率fl和上限频率fh。使用傅里叶变换算法将x'(n)变换到频域，得到频谱x'(f)；根据滤波器下限频率fl和上限频率fh，设计频域滤波器，对信号进行频域滤波，得到进行分频段处理后的信号x*(n)。
3)对经过分频带处理的数据进行时频分析
短时傅里叶变换是一种工程中常用的分析非平稳信号的手段，其具体步骤见图2，说明如下：
利用窗函数截取从t0时刻到t0+ts时刻信号，假定信号在窗内是平稳的，采用傅里叶变换来分析窗内信号，以确定在那个时段中信号所包含的频率信息，将得到的结果保存在一个矩阵sx(m,r)的列中(矩阵的每一列物理意义是ts长的信号所包含的频率信息，每一行的物理意义为该段信号在某个特定的频点上的时间分布)，即图2中矩阵的阴影部分，然后沿着信号时间轴将窗函数移动τs截取下一段ts时间的信号(这段信号和上一次截取的信号之间可以重叠，本发明中重叠率为50％，即τ＝0.5×ts)，得到下一段信号频率随时间的变换关系，保存在矩阵sx(m,r)中的下一列，以此类推，直到这段信号被完全处理。
其定义为：
式中，x*(n)为经过分频段后的时域数据，n的值等于ts×fs，w(n)为窗函数，其数学表示为：
sx(m,r)中的m代表对信号进行截取的时刻所对应的点数，等于a×τ×fs，其中a为大于等于0的整数且满足a×τ＜t，τ为每次进行帧移的长度，本发明中τ＝0.5×ts，r代表短时傅里叶变换后的频率。
需要注意的是，本发明中，对于不同频带的信号，在时频变换中的ts是不同的，ts直接影响了对该频带信号的时频分析分辨率。一般，对低频段信号的分析选择较大的ts(1-2s)，从而具有较高的频率分辨率；对高频段信号的分析选择较小的ts(0.01-0.05s)，从而具有较高的时间分辨率。
4)对时频分析结果进行恒虚警检测。
图3为一维恒虚警检测所用检测窗示意图，检测窗包括三部分：检测单元、保护单元、参考单元。其中检测单元xij位于检测窗中央，相邻的左右两边为保护单元(即图3中斜线标出部分)，左右两边保护单元的外侧分布着参考单元rij(即图3中十字阴影部分)。
恒虚警检测的基本思路是：按照一维恒虚警检测窗的结构对数据进行处理，得到相关参考单元、保护单元、检测单元。当某一单元中含有信号和噪声时，该单元所包含的能量大于其他仅含有噪声的单元中的能量。因此，使用参考单元中的能量对背景噪声能量进行估计，估计值和检测单元能量进行比较，从而确定检测单元中是否含有信号。本发明中恒虚警检测的具体流程见图4，详细说明如下：
设置恒虚警检测所用的参数，包括：虚警率、检测单元个数、参考单元个数、保护单元个数。本例中，虚警率为10-4，参考单元个数为60，检测单元个数为1，保护单元个数为4。
根据一维恒虚警检测窗提取上述短时傅里叶变换结果sx(m,r)中参考单元rij、检测单元的数据tij，具体方法是：提取sx(m,r)中的第i行，第j列元素作为检测单元的数据tij，按照一维恒虚警检测窗结构提取检测单元两侧的参考单元rij(下标i和j分别代表检测单元在矩阵sx(m,r)中的行数和列数)和保护单元。对这些数据进行如下处理：
对rij中的结果进行由大到小排序得到
式中，m为参考单元个数。针对检测单元tij的检测阈值t(i,j)的表达式如下：
t(i,j)＝α×ttemp(i,j)
其中，α的值和虚警率pfa之间关系如下：
pfa＝1-iα(n,m)
式中，n为检测单元个数，m为参考单元个数，iα(n,m)为不完全归一化beta函数，其定义如下：
式中，γ(n)为伽马函数。
对所有检测单元tij平方求和，本发明中因为检测单元只有1个，有：
γij＝(tij)2
进行如下比较
如果γij大于t(i,j)，说明该检测单元中有信号，xcfar(i,j)＝1；否则没有，xcfar(i,j)＝0。
经过恒虚警检测后，得到尺寸和sx(m,r)矩阵一致的矩阵xcfar，其每个元素只等于1或是0，其中1代表该单元有信号，0代表该单元没有信号。
相关常见的水声信号经过恒虚警检测后形状见图5，图中黑色点代表该点有信号。
5)计算恒虚警检测频率边际谱并判决。
计算恒虚警检测结果xcfar的频率边际谱xpf(f)，即对xcfar进行时间积分，其表达式为：
对凹凸不平的频率边际谱进行滤波平滑，选择合适的阈值tf，进行如下判断：
式中，门限值tf由虚警概率pfa和时间积分长度(即图2)矩阵sx(m,r)的列数)k决定，具体表达式为
tf＝a×k×pfa
根据实际工程中的情况，通过调整a适当对tf进行放大，保证检测有效性。
如果在一段频率范围内只有一个频点上的能量超过了阈值，则认为该频点上存在线谱信号，该频点认为是线谱信号的中心频率f0。如果在一段频率范围内的频点都超过了阈值，则认为该频率范围内存在块状信号，该频率范围的上下限认为是块状信号的频率上下限f′h和f′l。
相关水声信号的频率边际谱示意图见图6。
6)对恒虚警检测结果进行频域滤波，计算时间边际谱并判决。
根据第五步中得到的上下限频率f′h和f′l，设计带通滤波器，本发明中滤波器的传递函数h(f)表达式如下：
将恒虚警检测结果xcfar(t,f)通过该滤波器进行处理，得到x'cfar(t,f)，其表达式为：
x'cfar(t,f)＝xcfar(t,f)*h(f)
计算x'cfar(t,f)的时间边际谱，即对x'cfar(t,f)进行频率积分，其表达式为：
对xpt(t)进行滤波平滑。类似频率端点参数估计，选择合适的阈值tt对xpt(t)进行如下检测判决：
式中，门限值tt由恒虚警检测的虚警概率pfa和频率积分长度(即图2)矩阵sx(m,r)的行数)l决定，具体表达式为
tt＝b×l×pfa
根据实际工程中的情况，通过调整b适当对tt进行放大，保证参数估计的有效性。
认为xpt(t)中第一个超过阈值和最后一个超过阈值的点为信号的起止点tstart和tend，图7为常见水声信号时间边际谱示意图。