0 引　言

水声应用中需要检测水声脉冲信号。由于各项参数具有未知性与随机性，准确探测脉冲信号仍是难题^[1]。频域能量聚集性强的信号，如单频（CW）和线性调频（LFM）脉冲信号等，其检测方法已无法满足对频域能量聚集性弱、形式多变脉冲信号的检测要求，如直扩信号（DSSS）等，因此需要深入研究^[2]。

常用的信号检测方法，如能量检测法，可以不依赖信号的先验信息，但低信噪比下性能下降。熵检测法对于检测非合作脉冲信号优势突出，具有普适性和稳健性，且运算量小、计算速度快、对水声环境条件及信号形式不敏感等^[3]。熵值类型一般分为时域熵和频谱熵，在熵检测法计算过程中量化方式又分为均匀量化和非均匀量化^[4]。噪声与信号的熵值不同，因此可以通过频谱熵实现对脉冲信号的检测，而对于频域能量聚集性弱的信号，循环谱检测法对此类信号具有良好的检测效果^[5]。

本文针对7种常用声呐脉冲信号，从理论上对熵检测法进行分析，提出频谱熵与循环谱联合检测方法，该方法对信号形式不敏感，具有普适性与稳定性，对频域能量聚集性弱的信号形式检测性能提升；并针对检测曲线滤波的需求，分析高斯滤波及双边滤波原理，提出改进双边高斯滤波方法；根据奈曼皮尔逊准则，在给定虚警率条件下进行蒙特卡洛仿仿真，验证了检测与滤波方法的合理性与可行性。

1 频谱熵与循环谱联合检测原理 1.1 非均匀量化频谱熵法原理

当信号中只有白噪声时，频谱谱线随机分布，偏于无序，熵值最大。当存在特定信号时，谱线分布偏向有序，熵值较小，所以能够根据熵值判断信号是否存在。信息熵定义为，一离散型随机变量$X$，取值空间为$R$，概率密度函数为：

$ p\left( x \right) = \Pr \left( {X = x} \right)。$

(1)

则随机变量$X$的熵值为：

$ {H_a}\left( X \right) = - \sum\limits_{x \in R} {p\left( x \right){{\log }_a}p\left( x \right)}。$

(2)

式中：底数$a = 2$时，熵的单位为比特（bit），本文熵值均以该底数进行计算。规定当$p\left( x \right) = 0$时，$p\left( x \right) \times {\log _a}p\left( x \right) = 0$。

已知接收信号二元假设模型为：

$ \left\{ {\begin{aligned} &{{H_0}:x\left( n \right) = w\left( n \right)}，\\ &{{H_1}:x\left( n \right) = s\left( n \right) + w\left( n \right)}。\end{aligned}n = 0,1, \cdots N - 1}，\right. $

(3)

式中：$w\left( n \right)$表示均值为0、方差为$\sigma _0^2$服从高斯分布的噪声；$x\left( n \right)$为接收到方差为$\sigma _x^2$的信号，$s\left( n \right)$表示${\mu _s}$为均值$\sigma _s^2$为方差的信号；n = 0, 1, ··· N-1，N为采样点数。对式(3)作离散傅里叶变换，并非均匀量化频谱幅值区间，量化后各区间频谱幅值服从等概率分布，即每个子区间中的幅值个数出现概率相同。接收信号为纯噪声时，$X\left( k \right)$服从瑞利分布，其累计分布函数为：

$ {F_X}\left( x \right) = 1 - \exp \left( { - \frac{{{y^2}}}{{2{\sigma ^2}}}} \right)。$

(4)

通过$ {F_X}\left( {{x_i}} \right) = \left( {i - 1} \right)/L $得到各子区间分界点：

$ {x_i} = \sqrt { - 2\ln \left( {1 - \frac{{i - 1}}{L}} \right)\sigma } ,i = 1,2, \cdots ,L。$

(5)

假设${H_0}$成立时，频谱$X\left( k \right)$幅值的熵值最大；假设${H_1}$成立时，由于频谱幅值偏离等概率分布，熵值变小。

熵值${H_L}\left( X \right)$作为检验统计量，对比门限$\gamma $，当${H_L}\left( X \right) \leqslant \gamma $时，${H_1}$成立，则信号存在；当${H_L}\left( X \right) > \gamma $时，${H_0}$成立，则信号不存在。

1.2 循环谱法原理

相关函数随时间成周期或多个周期变化的信号统称为循环平稳信号^[6]。与平稳信号相似，循环平稳信号同样具有相关函数和谱密度函数，其相关函数的傅里叶序列可表示为：

$ {R_x}\left( {t + \frac{\tau }{2},t - \frac{\tau }{2}} \right) = \sum\limits_\alpha {R_x^a\left( \tau \right){e^{j2{\text{π}} \alpha t}}}。$

(6)

式中：$R_x^\alpha \left( \tau \right)$为循环自相关函数，表征时延为$\tau $时${R_x}$在频率$\alpha $处的强度，称为周期频率集合。循环频率$\alpha = 0$时，$R_x^0\left( \tau \right)$为普通自相关函数，体现平稳特性；$\alpha \ne 0$时，体现循环平稳特性。循环平稳信号$R_x^\alpha \left( \tau \right)$为：

$ R_x^\alpha \left( \tau \right) = \mathop {\lim }\limits_{T \to \infty } \frac{1}{T}\int_{ - \frac{T}{2}}^{\frac{T}{2}} {x\left( {t + \frac{\tau }{2}} \right){x^*}} \left( {t - \frac{\tau }{2}} \right){e^{ - j2{\text{π}} \alpha t}}\mathrm{d}t。$

(7)

令$ u\left( t \right) = x\left( t \right){e^{ - j{\text{π}} \alpha t}} $、$ v\left( t \right) = x\left( t \right){e^{j{\text{π}} \alpha t}} $，可将$R_x^\alpha \left( \tau \right)$看作$ u\left( t \right) $和$ v\left( t \right) $的互相关函数：

$ {R_{uv}} = R_x^\alpha \left( \tau \right) = \mathop {\lim }\limits_{T \to \infty } \frac{1}{T}\int_{ - \frac{T}{2}}^{\frac{T}{2}} {u\left( {t + \frac{\tau }{2}} \right){v^*}} \left( {t - \frac{\tau }{2}} \right)\mathrm{d}t。$

(8)

$u\left( t \right)$和$v\left( t \right)$的互谱密度函数为${S_{UV}}\left( f \right)$，再由互谱原理可得：

$ {S_x^\alpha = \mathop {\lim }\limits_{T \to \infty } \mathop {\lim }\limits_{\Delta t \to \infty } {S_{U{V_T}}}{\left( f \right)_{\Delta t}} = \mathop {\lim }\limits_{T \to \infty } \mathop {\lim }\limits_{\Delta t \to \infty } \dfrac{1}{{\Delta t}}\displaystyle\int_{ - \frac{{\Delta t}}{2}}^{\frac{{\Delta t}}{2}} {S_{{X_T}}^\alpha \left( {t,f} \right)\mathrm{d}t},} $

(9)

$ S_{{X_T}}^\alpha \left( {t,f} \right) \triangleq \frac{1}{T}{X_T}\left( {t,f + \frac{\alpha }{2}} \right)X_T^*\left( {t,f - \frac{\alpha }{2}} \right)，$

(10)

$ {X_T}\left( {t,f} \right) \triangleq \int_{t - \frac{T}{2}}^{t + \frac{T}{2}} {x\left( u \right){e^{ - j2{\text{π}} fu}}\mathrm{d}u}。$

(11)

式中：$ S_{{X_T}}^\alpha \left( {t,f} \right) $为时变周期的周期图，${X_T}\left( {t,f} \right)$为短时傅里叶变换频谱。$ S_x^\alpha \left( f \right) $相当于在频率$\left( {f + \alpha /2} \right)$和$\left( {f - \alpha /2} \right)$处求得信号$x\left( t \right)$的频谱分量之间的谱相关程度。

相比常规功率谱，循环谱在周期频率处能够表现信号特征。在非零周期频率处，平稳信号和噪声的谱相关值恒为0，但有用信息不为0。这使得即使有干扰和噪声存在，依旧能有效判断信号存在。

1.3 联合检测法原理

接收到采样率为${f_s}$的时域信号$x\left( t \right)$，需确定短时滑动窗长${\tau _\Delta }$，频率分辨率${\Delta _f}$，根据${\Delta _f}$计算快速傅里叶变换长度。对滑动窗内时域信号${x_i}\left( t \right)$补零，补零后信号为${x_{e{x_i}}}\left( t \right)$，而后进行快速傅里叶变换，得到当前信号段频谱${X_i}\left( f \right)$。

确定非均匀量化阶数$L$，计算当前信号段频谱标准差${\sigma _i}$，根据式(5)、$L$和${\sigma _i}$计算得到非均匀量化边界点向量${{\boldsymbol{L}}_i}$，进行归一化，得到归一化边界点向量$ {{\boldsymbol{\hat L}}_i} $。通过频谱幅值最大值${X_{i\max }}$、最小值${X_{i\min }}$以及归一化边界点向量$ {{\boldsymbol{\hat L}}_i} $计算频谱实际边界点向量${{\boldsymbol{L}}_{{X_i}}}$：

$ {{\boldsymbol{L}}_{{X_i}}} = {X_{i\min }} + \left( {{X_{i\max }} - {X_{i\min }}} \right){L_{ij}}。$

(12)

计算当前滑动窗信号频谱熵值${H_i}\left( X \right)$：

$ {H_i}\left( X \right) = - \sum\limits_{j = 1}^L {{p_j}{{\log }_2}{p_j}} ,j = 1,2, \cdots ,L。$

(13)

重复上述步骤，得到整段信号熵值向量$ \boldsymbol{{H}} $。对整段信号频谱熵计算结果进行高度门限${T_h}$和宽度门限${T_w}$判定，判定方式如下：

${y_{{h_i}}}\left[ H \right] = \left\{ {\begin{aligned} &{1,{H_{bg{f_i}}} \leqslant {T_h}}，\\ &{0,{H_{bg{f_i}}} > {T_h}} 。\end{aligned}} \right. $

(14)

对于小于高度门限的位置，在判定结果向量对应位置写入1，对所有连续置1段长度${L_{{y_h} = 1}}$进行宽度门限判定：

$ {y}_{{w}_{i}}\left[H\right]= \left\{\begin{aligned} &1,{T}_{w\mathrm{min}}\leqslant {L}_{{y}_{h}=1}\leqslant {T}_{w\mathrm{max}}，\\ &0.5,\beta {T}_{w\mathrm{min}}\leqslant su{m}_{{y}_{h}=1} < {T}_{w\mathrm{min}}，\\ &0,{L}_{{y}_{h}=1} > {T}_{w\mathrm{max}}或su{m}_{{y}_{h}=1} < \beta {T}_{w\mathrm{min}}。\end{aligned} \right.$

(15)

式中：${T_{w\min }}$为宽度门限下限；${T_{w\max }}$为宽度门限上限；$\beta $为系数且$0 < \beta < 1$；$su{m_{{y_h} = 1}}$为置1的总个数。若${y_{{w_i}}}\left[ H \right] = 1$，为当前脉宽符合宽度门限，存在脉冲信号，结束检测。若${y_{{w_i}}}\left[ H \right] = 0$，为当前脉宽小于宽度门限下限，不存在脉冲信号，结束检测。

若${y_{{w_i}}}\left[ H \right] = 0.5$，则表示熵值无法连续稳定小于高度门限，因此出现熵值在高度门限上下浮动的情况，从而导致在判定结果间断出现1的情况，累计检测结果置1的总数，经过二次宽度门限判定，若符合二次宽度门限，则认为疑似存在脉冲信号，需对当前段接收信号进行二次检测。

确定循环频率范围$\left[ { - \alpha ,\alpha } \right]$（一般为$\left[ { - {f_s}/2,{f_s}/2} \right]$）以及循环频率步长${\alpha _0}$，则需计算$N = 2\alpha /{\alpha _0} + 1$个循环频率处的频谱。时延当前滑动窗内时域信号${x_i}\left( t \right)$，将时间轴调整为从0时刻开始，得到该段信号从0时刻开始的时间长度${t_i}$，计算当前信号段的循环谱。

当前循环频率为${\alpha _k}$，可得到当前时域信号频移${\alpha _k}$的信号${u_i}\left( t \right)$及其复共轭信号${v_i}\left( t \right)$：

$ \left\{ {\begin{aligned} &{{u_i}\left( t \right) = {x_i}\left( t \right){e^{ - j{\text{π}} {\alpha _k}{t_0}}}}，\\ &{{v_i}\left( t \right) = {x_i}\left( t \right){e^{j{\text{π}} {\alpha _k}{t_0}}}}。\end{aligned}} \right. $

(16)

计算当前互相关函数$R_{u{v_i}}^{{\alpha _k}}\left( \tau \right)$：

$ R_{u{v_i}}^{{\alpha _k}}\left( \tau \right) = \mathop {\lim }\limits_{T \to \infty } \frac{1}{T}\int_{ - \frac{T}{2}}^{\frac{T}{2}} {{u_i}\left( {t + \frac{\tau }{2}} \right)v_i^*} \left( {t - \frac{\tau }{2}} \right)\mathrm{d}t 。$

(17)

傅里叶变换得到互谱密度函数$S_{_{U{V_i}}}^{{\alpha _k}}\left( f \right)$：

$ S_{_{U{V_i}}}^{{\alpha _k}}\left( f \right) = \int_{ - \infty }^\infty {R_{u{v_i}}^{{\alpha _k}}\left( \tau \right) \cdot {e^{ - j2{\text{π}} {\alpha _k}\tau }}\mathrm{d}\tau }。$

(18)

重复式(16)～式(18)，计算完整循环频率轴上的循环谱互谱密度函数，得到当前信号段的三维循环谱密度函数$ {S_{{C_i}}}\left( {f,\alpha } \right) $，$x$和$y$坐标为频率和循环频率，沿频率轴对循环谱密度函数求和，并取正频率轴范围，即可得到当前信号段二维循环频谱$ {S_{{C_i}}}\left( f \right) $：

$ {S_{{C_i}}}\left( f \right) = \sum\limits_{i = 1}^N {{S_{{C_i}}}\left( {f,\alpha } \right)} ,i = 1,2, \cdots ,N。$

(19)

计算频谱特征平滑度$S{S_{{C_i}}}$和峰度${K_{{C_i}}}$：

$ {S{S_{{C_i}}} = \displaystyle\sum\limits_{j = 1}^{N - 1} {{{\left( {{S_{{C_{i,j + 1}}}}\left( f \right) - {S_{{C_{i,j}}}}\left( f \right)} \right)}^2}} ,j = 1,2, \cdots ,M - 1，} $

(20)

$ {K_{{C_i}}} = \frac{{E\left[ {{{\left( {{S_{{C_i}}}\left( f \right) - {\mu _{{C_i}}}} \right)}^4}} \right]}}{{{{\left\{ {E\left[ {{{\left( {{S_{{C_i}}}\left( f \right) - {\mu _{{C_i}}}} \right)}^2}} \right]} \right\}}^2}}}。$

(21)

式中：$i$为当前滑动窗；$j$为当前循环谱频谱幅值；$M$为幅值总数。若$S{S_{{C_i}}} \geqslant {K_{{C_i}}}$，认为当前信号段为噪声，则计算当前信号段常规频谱熵值；若$S{S_{{C_i}}} < {K_{{C_i}}}$，认为当前信号段存在脉冲信号，则对当前循环谱频谱进行频谱优化。

循环谱频谱极大值中的最大值${M_{\max }}$，即当前循环谱频谱最大峰值。最大峰值邻域内左右最邻近极小值${m_l}$和${m_r}$，即循环谱频谱完整最大峰。

计算当前循环谱频谱均值${\mu _{{S_{Ci}}}}$和标准差${\sigma _{{S_{Ci}}}}$，确定系数$\gamma $取值。对完整最大峰外大于$\gamma \cdot {M_{\max }}$的频谱幅值${S_{{C_i}}}\left( {{f_j}} \right)$进行更新：

$ {S_{{C_i}}}\left( {{f_j}} \right) = {\mu _{{S_{Ci}}}} + U \cdot {\sigma _{{S_{Ci}}}}。$

(22)

式中：$U$为均匀分布在$\left[ {0,1} \right]$区间的随机变量${S_{{C_i}}}\left( {{f_j}} \right)$更新为当前循环谱均值与标准差缩放随机数的和。

使用归一化边界点向量$ {{\boldsymbol{\hat L}}_i} $计算当前循环谱频谱幅值的边界点向量${{\boldsymbol{L}}_{{C_i}}}$，从而得到当前滑动窗信号循环谱熵值${H_C}_i\left( X \right)$：

$ {H_C}_i\left( X \right) = - \sum\limits_{j = 1}^L {{p_C}_j{{\log }_2}{p_C}_j} ,j = 1,2, \cdots ,L。$

(23)

重复式(16)～式(23)，计算得到整段信号二次检测循环谱熵值向量${{\boldsymbol{H}}_C}$。进行高度门限${T_C}_h$和宽度门限${T_C}_w$判定，判定方式如下：

$ {y_{{h_i}}}\left[ {{H_C}} \right] = \left\{ {\begin{aligned} &{1,{H_C}_{bg{f_i}} \leqslant {T_C}_h}，\\ &{0,{H_C}_{bg{f_i}} > {T_C}_h}。\end{aligned}} \right. $

(24)

对于小于高度门限的位置，判定结果置1，对所有连续置1段长度${L_{{y_h} = 1}}$进行宽度门限判定：

$ {y}_{{w}_{i}}\left[{H}_{C}\right]=\left\{\begin{aligned} &1,{T}_{C}{}_{w\mathrm{min}}\leqslant {L}_{{y}_{h}={1}_{k}}\leqslant {T}_{C}{}_{w\mathrm{max}}，\\ &0,其他。\end{aligned} \right.$

(25)

式中：${T_C}_{w\min }$为宽度门限下限，${T_C}_{w\max }$为宽度门限上限。若${y_{{w_i}}}\left[ {{H_C}} \right] = 1$，为当前脉宽符合宽度门限，熵检测结果中存在脉冲信号，结束检测。若${y_{{w_i}}}\left[ {{H_C}} \right] = 0$，为当前脉冲信号宽度小于宽度门限下限，熵检测结果中不存在所需脉冲信号，结束检测。

随着信噪比逐渐降低，频谱熵检测法对频域能量聚集性弱的信号检测性能下降，频谱熵与循环谱联合检测方法能够提高其检测性能。

2 改进双边高斯滤波原理

短时频谱熵检测曲线表示熵值随时间的变化特性，反映了脉冲信号的起止时刻及熵值大小，理想的滤波效果为检测曲线平滑且保留信号边缘突变特征。信噪比降低会加剧熵值数值抖动，难以设置门限。

根据中心极限定理，采样点数较大时，噪声近似为高斯噪声，使用高斯滤波器进行低通滤波可以得到较好的平滑效果。高斯滤波常用于消除高斯噪声，为线性滤波。整段待滤波数据每一点的值，由其自身及其邻域内其他点加权平均后得到。多使用滤波窗长较长的高斯滤波器，以达到较好滤波效果。

双边滤波是一种非线性滤波器，具有局部性、非迭代以及设计简单等特点。对局部待滤波数据加权平均，待滤波点一维邻域内的加权系数由2个因子相乘组成，一是由待滤波数据的空间距离决定的空间邻近度因子${\sigma _s}$（或称空间权重）；二是待滤波数据数值之差决定的值域相似度因子${\sigma _r}$（或称灰度权重）。灰度权重无法由一个固定值满足不同信噪比的滤波需求。滤波窗长较小时，空间权重接近1，未充分发挥其作用；滤波窗长较大时，则计算量增加。双边滤波能有效保留边缘特征，滤波过程依赖于局部信息，适用于处理动态变化的信号^[7]。

将高斯滤波与双边滤波进行结合，得到改进双边高斯滤波方法。具体滤波步骤如下：

对接收信号进行短时频谱熵检测，得到检测结果频谱熵${\boldsymbol{H}}$，该一维向量中的每一个元素代表当前滑动窗内信号的频谱熵值。

设定滤波函数中的各项参数初始值，包括基础滤波窗长${L_0}$、高斯滤波器基础标准差${\sigma _{g0}}$、双边滤波器中的基础空间标准差${\sigma _{s0}}$、基础值域标准差${\sigma _{r0}}$。

根据短时熵检测窗长${L_t}$确定局部时间尺度因子$ {\tau _\Delta } = {\alpha _1}{L_t} $，$ {\alpha _1} $为系数。

确定改进双边高斯滤波窗长$L$范围：

$ \left\{ {\begin{aligned} &{{L_{\min }} = {L_0}\left( {1 + {\alpha _2} \cdot {a^{ - {\tau _\Delta }}}} \right)} ，\\ &{{L_{\max }} = {\alpha _3}{L_{\min }} + c}。\end{aligned}} \right. $

(26)

式中：${\alpha _2}$为系数，$a$为底数，${L_{\min }}$为滤波窗长最小值，滤波窗长需为奇数，则取${L_{\min }}$最邻近奇数；${\alpha _3}$为系数，$c$为常数，取${L_{\max }}$最邻近奇数。

设置初始滤波窗口大小$L = {L_{\min }} + 2$，使用$L$滑动计算整段信号熵值${{H}}$的局部样本方差$ {\mathbf{\sigma }}_\Delta ^2 $，${\mu _{{{\boldsymbol{H}}_i}}}$为第$i$个窗口内熵值${{H}}$的均值，${H_{ij}}$表示第$i$个窗口内的第$j$个熵值，$\sigma _{\Delta i}^2$表示第$i$个窗口的方差，共$i$个方差$\sigma _{\Delta i}^2$组成整段熵值的局部样本方差向量$ {\mathbf{\sigma }}_\Delta ^2 $，并计算$ {\mathbf{\sigma }}_\Delta ^2 $的中值平均值$ {\mu _{{\sigma ^2}}} $，对$ \sigma _{\Delta 1}^2,\sigma _{\Delta 2}^2, \cdots ,\sigma _{\Delta L}^2 $进行升序排序后取中值平均数。$ {\mathbf{\sigma }}_\Delta ^2 $和$ {\mu _{{\sigma ^2}}} $用来评价整段熵值的平滑度。

计算整段熵值${{H}}$的一阶局部变化率${{D}}$：

$ {{D}} = {\left[ {\left| {{H_2} - {H_1}} \right|, \cdots ,\left| {{H_i} - {H_{i - 1}}} \right|, \cdots ,\left| {{H_L} - {H_{L - 1}}} \right|} \right]^\mathrm{T}}。$

(27)

基础参数向量为${{\mathbf{\sigma }}_0} = {\left[ {{\sigma _{g0}},{\sigma _{s0}},{\sigma _{r0}}} \right]^\mathrm{T}}$，根据${{\mathbf{\sigma }}_0}$及$ {\mu _{{\sigma ^2}}} $计算各滤波器参数，包括高斯滤波器标准差的最大值${\sigma _{g\max }}$和最小值${\sigma _{g\min }}$、双边滤波器空间标准差的最大值${\sigma _{s\max }}$和最小值${\sigma _{s\min }}$、双边滤波器值域标准差的最大值${\sigma _{r\max }}$和最小值${\sigma _{r\min }}$。得到最小值参数向量${{\mathbf{\sigma }}_{\min }} = {\left[ {{\sigma _{g\min }},\;{\sigma _{s\min }},\;{\sigma _{r\min }}} \right]^\mathrm{T}}$，最大值参数向量${{\mathbf{\sigma }}_{\max }} = [ {{\sigma _{g\max }}, {\sigma _{s\max }},{\sigma _{r\max }}} ]^\mathrm{T}$。共分为以下3种情况：①$ 0 \leqslant {\mu _{{\sigma ^2}}} \leqslant {m_1} $、②$ {m_1} < {\mu _{{\sigma ^2}}} \leqslant {m_2} $、③$ {\mu _{{\sigma ^2}}} > {m_2} $。

$i$为第$i$种情况，滤波器参数最值可由式(28)、式(29)计算：

$ \left[ {\begin{array}{*{20}{c}} {{\sigma _{g\min }}} \\ {{\sigma _{s\min }}} \\ {{\sigma _{r\min }}} \end{array}} \right] = \left[ {\begin{array}{*{20}{c}} {{\beta _{i1}}}&0&0 \\ 0&{{\beta _{i2}}}&0 \\ 0&0&{{\beta _{i3}}} \end{array}} \right] \cdot \left[ {\begin{array}{*{20}{c}} {{\sigma _{g0}}} \\ {{\sigma _{s0}}} \\ {{\sigma _{r0}}} \end{array}} \right] + \left[ {\begin{array}{*{20}{c}} {{\gamma _{i1}}} \\ {{\gamma _{i2}}} \\ {{\gamma _{i3}}} \end{array}} \right]，$

(28)

$ {\left[ {\begin{array}{*{20}{c}} {{\sigma _{g\max }}} \\ {{\sigma _{s\max }}} \\ {{\sigma _{r\max }}} \end{array}} \right] = \left[ {\begin{array}{*{20}{c}} {{\beta _{i4}}}&0&0 \\ 0&{{\beta _{i5}}}&0 \\ 0&0&{{\beta _{i6}}} \end{array}} \right] \cdot \left[ {\begin{array}{*{20}{c}} {{\sigma _{g\min }}} \\ {{\sigma _{s\min }}} \\ {{\sigma _{r\min }}} \end{array}} \right] + \left[ {\begin{array}{*{20}{c}} {{\gamma _{i4}}} \\ {{\gamma _{i5}}} \\ {{\gamma _{i6}}} \end{array}} \right]。} $

(29)

式中：${\boldsymbol{\beta }}$为系数矩阵；${\boldsymbol{\gamma }}$为常数向量。

为防止窗长变化过大，需减缓窗长增长速率，对熵值${{H}}$的局部变化率$D$取对数${D_{\log }} = \lg \left( {1 + D} \right)$。

根据${D_{\log }}$、${L_{\min }}$、${\sigma _{g\min }}$、${\sigma _{g\max }}$、${\sigma _{s\min }}$、${\sigma _{s\max }}$、${\sigma _{r\min }}$、${\sigma _{r\max }}$，计算每个熵值点$ {H_i} $的滤波窗长向量${{\boldsymbol{L}}_w}$、高斯标准差向量${{\mathbf{\sigma }}_g}$、双边空间标准差向量${{\mathbf{\sigma }}_s}$和值域标准差向量${{\mathbf{\sigma }}_r}$，计算方法如下：

$ \left\{ {\begin{aligned} &{{{\boldsymbol{L}}_w} = {L_{\min }} + \left( {{L_{\max }} - {L_{\min }}} \right) \cdot \left( {1 - \displaystyle\frac{{{D_{\log }}}}{{{D_{{{\log }_{\max }}}}}}} \right)}，\\ &{{{\mathbf{\sigma }}_g} = {\sigma _{g\min }} + \left( {{\sigma _{g\max }} - {\sigma _{g\min }}} \right) \cdot \left( {1 - \displaystyle\frac{{{D_{\log }}}}{{{D_{{{\log }_{\max }}}}}}} \right)} ，\\ &{{{\mathbf{\sigma }}_s} = {\sigma _{s\min }} + \left( {{\sigma _{s\max }} - {\sigma _{s\min }}} \right) \cdot \left( {1 - \displaystyle\frac{{{D_{\log }}}}{{{D_{{{\log }_{\max }}}}}}} \right)}，\\ &{{{\mathbf{\sigma }}_r} = {\sigma _{r\min }} + \left( {{\sigma _{r\max }} - {\sigma _{r\min }}} \right) \cdot \displaystyle\frac{{{D_{\log }}}}{{{D_{{{\log }_{\max }}}}}}} 。\end{aligned}} \right. $

(30)

取当前滤波窗的高斯标准差${\sigma _{gi}}$，构建高斯滤波器$ {\boldsymbol{g}}{{\boldsymbol{f}}_i} $：

$ {\boldsymbol{g}}{{\boldsymbol{f}}_i} = \Bigg[ {{e^{ - \frac{{L_{wi1}^2}}{{2\sigma _{gi}^2}}}},{e^{ - \frac{{L_{wi2}^2}}{{2\sigma _{gi}^2}}}}, \cdots ,{e^{ - \frac{{L_{wij}^2}}{{2\sigma _{gi}^2}}}}, \cdots ,{e^{ - \frac{{L_{wiN}^2}}{{2\sigma _{gi}^2}}}}} \Bigg]。$

(31)

式中：$i$为当前滤波窗，窗长${L_{wi}} = N$，${L_{wij}}$为第$i$个窗内的第$j$个值。计算当前高斯滤波器的总值$ {S_{{\boldsymbol{g}}{{\boldsymbol{f}}_i}}} $，对当前高斯滤波器每一个点计算权重，得到当前滤波窗的高斯权重向量${{\boldsymbol{w}}_{g{f_i}}}$

取与当前高斯滤波器长度${L_{wi}}$相等的第$i$点熵值${H_i}$邻域内的$N$个点组成当前滤波窗口熵值向量${{\boldsymbol{H}}_i}$，对该熵值向量进行高斯滤波得到高斯滤波后的熵值向量$ {{\boldsymbol{H}}_{g{f_i}}} $。

重复上述步骤，对所有滤波窗都进行高斯滤波，得到整段熵值向量高斯滤波后的结果向量${{\boldsymbol{H}}_{gf}}$。

在高斯滤波的基础上进行双边滤波，双边滤波时需对熵值向量边缘突变特征进行保留，因此需确定边缘检测门限。计算局部变化率向量${\boldsymbol{D}}$的均值${\mu _{\mathbf{D}}}$及标准差${\sigma _{\mathbf{D}}}$。根据${\mu _{\mathbf{D}}}$及${\sigma _{\mathbf{D}}}$得到边缘检测门限${T_{edge}} = {\mu _{\mathbf{D}}} + {\sigma _{\mathbf{D}}}$。若$ D_i > T_{edge} $，则该点需要边缘保护。

当前滤波窗对应的空间标准差${\sigma _{si}}$和值域标准差${\sigma _{ri}}$。若$ D_i > T_{edge} $，${\sigma _{ri}}$更新为${\sigma _{r\min }}$。构建当前滤波窗空间滤波器$ {\boldsymbol{bf}}{_{{s_i}}} $：

$ {\boldsymbol{b}}{{\boldsymbol{f}}_{{s_i}}} = \Bigg[ {{e^{ - \frac{{L_{wi1}^2}}{{2\sigma _{si}^2}}}},{e^{ - \frac{{L_{wi2}^2}}{{2\sigma _{si}^2}}}}, \cdots ,{e^{ - \frac{{L_{wij}^2}}{{2\sigma _{si}^2}}}}, \cdots ,{e^{ - \frac{{L_{wiN}^2}}{{2\sigma _{si}^2}}}}} \Bigg]。$

(32)

对未滤波熵值向量${\boldsymbol{H}}$进行滤波窗长最大值${L_{\max }}$半窗长$\left\lfloor {{L_{\max }}/2} \right\rfloor $扩展（$ \lfloor ·\rfloor $为向下取整），即首尾对称填充$\left\lfloor {{L_{\max }}/4} \right\rfloor $个元素，以保留可能存在首尾的突变特征和边界信息，减少边缘效应。得到扩展后的扩展熵值向量${{\boldsymbol{H}}_{ex}}$。

找出当前滤波窗对应高斯滤波后的熵值${{\boldsymbol{H}}_{g{f_i}}}$，当前滤波点位$i$对应扩展熵值向量${{\boldsymbol{H}}_{ex}}$的同点位$i$，从$i$点取当前滤波窗长${L_{wi}}$的扩展熵值向量${{\boldsymbol{H}}_{e{x_i}}}$，结合当前值域标准差${\sigma _{ri}}$构建当前值域滤波器$ {\boldsymbol{bf}}{_{{r_i}}} $：

$ {{{\boldsymbol{b}}{{\boldsymbol{f}}_{{r_i}}} = \Bigg[ {{e^{ - \frac{{{{\left( {{H_{ex1}} - {H_{g{f_i}}}} \right)}^2}}}{{2\sigma _{ri}^2}}}}, \cdots ,{e^{ - \frac{{{{\left( {{H_{exj}} - {H_{g{f_i}}}} \right)}^2}}}{{2\sigma _{ri}^2}}}}, \cdots ,{e^{ - \frac{{{{\left( {{H_{exN}} - {H_{g{f_i}}}} \right)}^2}}}{{2\sigma _{ri}^2}}}}} \Bigg] }。}$

(33)

根据当前空间滤波器$ {\boldsymbol{bf}}{_{{s_i}}} $和值域滤波器$ {\boldsymbol{b}}{{\boldsymbol{f}}_{{r_i}}} $计算双边滤波器$ {\boldsymbol{b}}{{\boldsymbol{f}}_{s{r_i}}} $：

$ \boldsymbol{b}\boldsymbol{f}_{sr_i}=\boldsymbol{b}\boldsymbol{f}_{s_i}\cdot\boldsymbol{bf}_{_{r_i}}^{\mathrm{T}}。$

(34)

计算当前双边滤波器$ {\boldsymbol{b}}{{\boldsymbol{f}}_{s{r_i}}} $的总值$ {S_{{\boldsymbol{b}}{{\boldsymbol{f}}_i}}} $，得到权重向量${{\boldsymbol{w}}_{s{r_i}}}$。

取与当前双边滤波器长度${L_{wi}}$相等的扩展熵值矩阵${{\boldsymbol{H}}_{ex}}$的第$i$点熵值${H_{e{x_i}}}$邻域$N$个点组成当前待滤波熵值向量${{\boldsymbol{H}}_{e{x_i}}}$，对该熵值向量进行双边滤波得到滤波后的熵值向量$ {{\boldsymbol{H}}_{b{f_i}}} $。

重复高斯滤波后的步骤，对所有高斯滤波后的熵值进行双边滤波，得到熵值经改进双边高斯滤波后的结果${{\boldsymbol{H}}_{bgf}}$。

3 仿真分析

仿真信号为7种不同形式的常用声呐信号，分别为单频矩形（CW）脉冲、线性调频（LFM）脉冲、双曲调频（HFM）脉冲、二进制相移键控调制（BPSK）脉冲、跳频（FHSS）脉冲、直接序列扩频（DSSS）脉冲、组合波形脉冲，其中组合波形脉冲为CW叠加LFM信号。脉冲信号脉宽在30～350 ms随机取值。

如图1所示为7种信号形式的时频图，仅展示信号类型及时频特性。以脉宽为300 ms示例，前5种信号频域能量聚集性强，后2种信号频域能量聚集性弱。后文信号形式出现顺序与图1一致。

对以上7种声呐信号使用不同信号检测方法进行检测，分别为非均匀量化频谱熵检测法、均匀量化频谱熵检测法、时域熵检测法、能量检测法、短时傅里叶检测法。对不同方法的检测统计量归一化后取幅度，使用固定参数高斯滤波对归一化幅度进行平滑处理，得到不同方法检测统计量归一化幅度对比结果。

如图2所示，可以看出非均匀量化频谱熵法检测性能最好，对于频域能量聚集性弱的信号，即直扩信号和组合波形信号，信号与噪声叠加段与纯噪声段检测量的差值小于频域能量聚集性强的信号差值，因此检测效果下降。在低信噪比条件下，时域信号被噪声淹没，时域熵检测性能较差，后续不再使用时域熵检测法进行对比验证。

针对于该种形式信号，使用上述提出的频谱熵与循环谱联合检测方法在−20～0 dB信噪比下进行验证，恒虚警率${P_f}$为0.001进行1000次蒙特卡洛实验得到不同形式脉冲信号的检测率曲线。

如图3所示，通过对7种形式脉冲信号的仿真验证结果对比，在相同恒虚警概率的条件下，循环谱频谱与频谱熵联合检测法的检测性能均优于其他检测方法，5种检测方法分别为：①循环谱与频谱熵联合法、②非均匀量化频谱熵法、③均匀量化频谱熵法、④能量检测法、⑤短时傅里叶变换法，当每种检测方法检测概率达到1时信噪比的值如表1所示。

可知，对比5种检测方法的检测能力，同信噪比下，循环谱与频谱熵联合检测法检测率最高。证明该信号检测方法的有效性与优势性。

在脉冲信号检测流程中，短时频谱熵检测步骤计算结果如图4所示，仿真条件为信噪比−6 dB，每种脉冲信号脉宽为300 ms，滑动窗长为4 ms，非均匀量化，量化阶数为16。

可以看出，短时频谱熵检测结果在一高值和低值交替浮动，由于纯噪声段只有噪声存在，混乱度高，因此熵值较大，在3.7 bit附近浮动，低于纯噪声熵值的7处为脉冲信号段。由于检测结果抖动较大，难以确定高度门限和宽度门限，因此需要对检测结果进行滤波处理。

使用3种滤波方法对熵值曲线进行滤波。滤波器参数设置为：固定参数高斯滤波，滤波窗长为41，高斯标准差${\sigma _g} = 5$。固定参数双边滤波，滤波窗长为41，空间标准差${\sigma _s} = 5$，值域标准差${\sigma _r} = 0.2$。改进双边高斯滤波，基础滤波窗长${L_0} = 5$、高斯滤波器基础标准差${\sigma _{g0}} = 2$，双边滤波器基础空间标准差${\sigma _{s0}} = 2$，基础值域标准差${\sigma _{r0}} = 0.1$，系数$ {\alpha _1} = 1000 $，系数${\alpha _2} = 1.5$，底数$a = 5$，系数${\alpha _3} = 5/3$，常数$c = 16$。参数分别为：

$ \left\{\begin{aligned} & diag\left(\beta_{11},\beta_{12},\beta_{13}\right)=diag\left(1,1,1\right)，\\ & \left[\gamma_{11},\gamma_{12},\gamma_{13}\right]^{\mathrm{T}}=\left[0,0,0\right]^\mathrm{T}，\\ & diag\left(\beta_{14},\beta_{15},\beta_{16}\right)=diag\left(2,1.5,4\right)，\\ & \left[\gamma_{14},\gamma_{15},\gamma_{16}\right]^\mathrm{T}=\left[0,0,0\right]^\mathrm{T}。\end{aligned}\right. $

(35)

$ \left\{\begin{aligned} & diag\left(\beta_{21},\beta_{22},\beta_{23}\right)=diag\left(2,1,1.5\right)，\\ & \left[\gamma_{21},\gamma_{22},\gamma_{23}\right]^{\mathrm{T}}=\left[0,0,0\right]^{\mathrm{T}}，\\ & diag\left(\beta_{24},\beta_{25},\beta_{26}\right)=diag\left(2,1.5,4\right)，\\ & \left[\gamma_{24},\gamma_{25},\gamma_{26}\right]^{\mathrm{T}}=\left[0,0,0\right]\mathrm{^T}。\end{aligned}\right. $

(36)

$ \left\{\begin{aligned} & diag\left(\beta_{31},\beta_{32},\beta_{33}\right)=diag\left(6,1,2\right)，\\ & \left[\gamma_{31},\gamma_{32},\gamma_{33}\right]\mathrm{^T}=\left[0,2,0\right]\mathrm{^T}，\\ & diag\left(\beta_{34},\beta_{35},\beta_{36}\right)=diag\left(2,1.5,4\right)，\\ & \left[\gamma_{34},\gamma_{35},\gamma_{36}\right]^{\mathrm{T}}=\left[0,0,0\right]\mathrm{^T}。\end{aligned}\right. $

(37)

3种方法滤波处理结果如图5所示，可以明显分辨脉冲信号段与纯噪声段，固定参数高斯滤波对于熵检测结果有良好的平滑效果，熵值数值大小较好保留。但脉冲信号起止时间产生模糊，加长脉宽，宽度门限判定时会产生较大偏差。检测曲线首末位置产生边界效应，边界处滤波效果不佳，与滤波前熵值不一致或失真。

固定参数双边滤波能够较好保留突变特征，边缘处理效果良好，但脉宽产生较大偏差。且滤波后熵值数值均变小，与滤波前熵值无法保持一致性，难以设置高度检测门限。

改进双边高斯滤波结合以上2种滤波方法优点，舍弃缺点，既达到良好的平滑效果，又保留了突变特征。脉冲起止时间以及熵值数值大小保持一致，消除了边界效应。

4 结　语

本文对7种常用声纳信号进行恒虚警率的检测性能仿真。对于频域能量聚集性弱的信号形式，频谱熵与循环谱联合检测方法性能提升较大；其他信号形式的检测性能也均有提升。针对检测曲线滤波问题，使用改进双边高斯滤波方法，既对曲线产生较好的平滑效果，又保留脉冲信号突变特征，达到了理想的滤波效果。仿真结果表明，相同信噪比条件下，本文所提方法的检测率均高于其他检测方法，验证了该信号检测方法的可行性与有效性。