推进器被动声学探测是水下目标识别的重要问题^[1-2]，其辐射噪声包括推进噪声、水动力噪声和机械噪声. 推进噪声蕴含着与桨叶旋转有关的重要信息^[3-6]，包括纯音成分和宽带成分2部分. 其中，纯音成分（即低频线谱噪声）主要源自厚度噪声、叶片-涡流相互作用噪声、载荷噪声及高速脉冲噪声；宽带成分主要源自空化噪声、湍流吸入噪声、边界层湍流噪声、涡流脱落噪声，其信号能量受到旋转桨叶的调制. 这2类信号成分可被频谱分析、信号解调技术所提取识别，故成为水声目标识别领域的关注重点.

目前主流的解调方法包括传统DEMON（detection of envelope modulation on noise）、窄带包络解调和循环平稳分析. 传统DEMON分析通过从辐射噪声中定位推进噪声所在频带来提取有效信息，较依赖于操作人员的主观判断^[7]，为此，学者们提出利用模糊逻辑代替人工选取^[8]并通过信号统计来抑制高斯噪声^[9] ，该方法已取得广泛应用^[10]，但频带优选准确度仍存在较大的发展空间. 窄带包络解调同样作为推进器调制线谱提取的主流技术之一，其核心在于选取包含推进噪声的最佳解调频带以提高解调性能^[11-12]. Kurtogram^[13]及FK^[14]利用信号冲击特性，借助峭度指标自动识别最佳频带，然而该方法易受非高斯噪声干扰而失效；Infogram^[15]利用香农熵作为指标，同时在时、频域中捕捉重复瞬态特征，达到稳健的解调性能，但作为全局性的统计量，熵值在目标频率冲击微弱时难以表征其特性；Autogram^[16]通过计算自相关平方包络的峭度来选取最佳解调频带，该方法无需目标频率信息，但易受到与目标频率无关的周期性干扰噪音的影响；IESFOgram^[17-18]在已知目标频率的前提下能有效定位最佳解调频带，但过于依赖先验信息，且其单频带输出方式容易在循环平稳干扰（如电磁激励噪声）下失效. 此外，窄带包络解调检测范围不能超过解调频带的带宽^[19]，否则可能导致轴叶频调制信息丢失，在低信噪比及含复杂噪声（冲击噪声、循环平稳噪声）的情况下识别效果不佳. 相较而言，循环平稳分析具有高抗噪性的优点. 由于叶片式水下推进器辐射噪声由具有一阶循环平稳性的纯音噪声和具有二阶循环平稳性的宽带噪声组成，而宽带噪声主导了推进器辐射噪声的循环平稳特性，故学者们^[6,20]利用增强包络谱（enhanced envelope spectrum, EES）来提取螺旋桨特征频率. EES的解调效果受限于积分频带上下限的选择，为了提升循环平稳分析在水下目标探测领域的适应性，学者们引入加权函数，对轴叶频调制线谱进行增强提取. 自适应加权包络谱（adaptive weighted envelope spectrum, AWES）^[21]可在无轴频先验信息条件下实现推进器轴叶频的增强解调提取，但该方法在循环平稳噪声干扰下易于失效. 在此基础上，周期稀疏包络谱（periodic sparsity envelope spectrum, PSES）^[3]可在已知轴频先验信息的条件下进一步增强对轴叶频线谱的解调提取性能，而该方法计算所需的经验参数较多，并且要求精准的轴频先验信息.

针对以上研究现状，本研究分析揭示水下推进器辐射噪声的循环平稳特性，设计频域相关峭度指标，构建有轴频先验信息下的谱相干加权函数，进而提出遍历加权包络谱（traversal weighted envelope spectrum，TWES）. 在此基础上，本研究构造自适应统计阈值曲线，提出特征峰噪比（characteristic peak-to-noise peak ratio，CNR）来量化评价解调方法，支撑推进器调制线谱的提取和识别.

水下推进器辐射噪声$ x(t) $可视为目标信号组分$ m(t) $和噪声干扰组分$ n(t) $的加和. 目标信号组分$ m(t) $是由推进器叶轮旋转引发的循环平稳信号，包括纯音成分（即一阶循环平稳成分$ {m}_{\text{CS1}}(t) $）和宽带成分（即二阶循环平稳成分$ {m}_{\text{CS2}}(t) $）. 噪声干扰组分$ n(t) $包含平稳噪声$ {n}_{\text{s}}(t) $（主要来自于测量环境和测试系统）、一阶循环平稳噪声$ {n}_{\text{CS1}}(t) $、二阶循环平稳噪声$ {n}_{\text{CS2}}(t) $和冲击噪声$ {n}_{\text{P}}(t) $，这些干扰组分主要由电磁激励噪声、机械部件激励、流场压力脉动等因素导致. 据此，可以建立推进器辐射噪声的循环平稳信号模型：

(1)$ \begin{split}x(t)=&{m}_{\text{CS1}}(t)+{m}_{\text{CS2}}(t)+{n}_{\text{CS1}}(t)+{n}_{\text{CS2}}(t)+\\ & {n}_{\text{P}}(t)+{n}_{\text{s}}(t).\end{split} $

将一阶循环平稳成分$ {m}_{\text{CS1}}(t) $和二阶循环平稳成分$ {m}_{\text{CS2}}(t) $转化成指数形式：

(2)$ \left.\begin{split}{m}_{\text{CS1}}(t)=&\displaystyle\sum\limits_{i=-I}^{I}{a}_{i}{\text{e}}^{\text{j}2{\text{π}} i\omega t},\\{m}_{\text{CS2}}(t)=&\displaystyle\sum\limits_{k=-K}^{K}{b}_{k}{\text{e}}^{\text{j}2{\text{π}} k\omega t}\cdot v\left(t\right). \end{split}\right\} $

式中：$ \omega $代表推进器轴频，$ {a}_{i} $代表$ {m}_{\text{CS1}}(t) $中第i个正弦信号成分的幅值，I代表$ {m}_{\text{CS1}}(t) $所包含的正弦信号成分最高阶次，$ {b}_{k} $代表$ {m}_{\text{CS2}}(t) $中第k个调制信号成分的幅值，K代表$ {m}_{\text{CS2}}(t) $中所包含的调制信号成分最高阶次，$ v\left(t\right) $代表随机载波.

窄带包络解调的核心思想在于选取一个或若干个最佳解调频带，提高信噪比进而实现调制线谱的增强提取. 这一思想可扩展到循环平稳分析中：对谱相关或谱相干在一个最佳频带或者若干个最佳频带内进行积分操作，从而增强其解调效果. 本研究提出更普适的增强解调策略：将谱相关或谱相干乘以一个权重函数，该权重函数针对目标信号组分的载波频带赋予较高值，对其他区域赋予较低值. 因此，如何确定信号增强解调的最佳加权函数成为本研究所关注的核心问题. 为了解决该问题，通过理论推导揭示推进器辐射噪声的循环平稳特性.

首先，时域信号$ x\text{(}t\text{)} $的瞬时自相关函数定义如下：

(3)$ R_x(t, \tau)={E}\left[x\left(t+\frac{\tau}{2}\right) x^*\left(t-\frac{\tau}{2}\right)\right] . $

式中：$ {E }[\cdot ] $代表集合平均算子，$ t $和$ \tau $分别代表时间和时延，$ \ast $代表复共轭. 在此基础上，可以进一步获取循环自相关函数：

(4)$ {R}_{x}\text{(}\alpha ,\tau \text{)=}\underset{T\rightarrow \mathrm{\infty }}{\lim }\;\frac{1}{T}\displaystyle{\int}_{-T/2}^{T/2}{R}_{x}\text{(}t,\tau \text{)}{\text{e}}^{-\text{j}2{\text{π}} \alpha t}\text{d}t. $

式中：$ \alpha $代表循环频率，即调制频率. 最后，可获取谱相关函数：

(5)$ {{\mathrm{SC}}}_{x}\text{(}\alpha ,f\text{)=}\displaystyle{\int}_{-\mathrm{\infty }}^{+\mathrm{\infty }}{R}_{x}\text{(}\alpha ,\tau \text{)}{\text{e}}^{-\text{j}2{\text{π}} f\tau }\text{d}\tau . $

式中：$ f $代表谱频率，即载波频率.

本研究利用推进器辐射噪声的循环平稳特性来实现水下推进器被动声学探测. 如1.1节所述，目标信号组分包括一阶循环平稳成分$ {m}_{\text{CS1}}(t) $和二阶循环平稳成分$ {m}_{\text{CS2}}(t) $，因此接下来对这2部分的谱相关理论值进行理论推导.

考察含有平稳噪声的一阶循环平稳成分$ {x}_{\text{CS1}}(t) $和二阶循环平稳成分$ {x}_{\text{CS2}}(t) $，即

(6)$ \left.\begin{split} {x}_{\text{CS1}}(t)=&\displaystyle\sum\limits_{i=-I}^{I}{a}_{i}{\text{e}}^{\text{j}2{\text{π}} i\omega t}+{n}_{\text{s}}(t),\\ {x}_{\text{CS2}}(t)=&\displaystyle\sum\limits_{k=-K}^{K}{b}_{k}{\text{e}}^{\text{j}2{\text{π}} k\omega t}\cdot v\left(t\right)+{n}_{\text{s}}(t).\end{split}\right\} $

将式（6）代入到式（3）~（5），可得到$ {x}_{\text{CS1}}(t) $和$ {x}_{\text{CS2}}(t) $的谱相关理论值为

(7)$ \left.\begin{split} {\rm{S C}}_{x_{\mathrm{CS} 1}}(\alpha, f)= & \displaystyle\sum_{i_1, i_2} a_{i_1} a_{i_2}^* \delta\left[\alpha-\left(i_1-i_2\right) \omega\right]\times \\&\delta\left[f - \frac{\left(i_1 + i_2\right) \omega}{2}\right] + {A} + {B} + P_{n_{\mathrm{s}}}(f) \delta(\alpha),\\{\rm{S C}}_{x_{\mathrm{CS} 2}}(\alpha, f) =&\displaystyle\sum_{k_1, k_2} b_{k_1} b_{k_2}^* \delta\left[\alpha-\left(k_1-k_2\right) \omega\right]\times \\& P_v\left[f-\frac{\left(k_1+k_2\right) \omega}{2}\right] +P_{n_{\mathrm{s}}}(f) \delta(\alpha).\end{split}\right\} $

式中：A、B分别代表$ \displaystyle\sum\limits_{{i}_{1}}{u}_{{{n}_{\text{s}}}}{a}_{{{i}_{1}}}\delta (\alpha -{i}_{1}\omega )\delta (f-{i}_{1}\omega /2) $、$ \displaystyle\sum\limits_{{i}_{2}}{u}_{{{n}_{\text{s}}}}a_{{i}_{2}}^{\ast }\delta (\alpha +{i}_{2}\omega )\delta (f-{i}_{2}\omega /2) $，$ {u}_{{{n}_{\text{s}}}} $和$ {P}_{{{n}_{\text{s}}}}(f) $分别代表平稳噪声$ {n}_{\text{s}}(t) $的期望值和功率谱密度；$p_v(f) $表示信号v(t)的功率谱密度，v(t)表示零均值的平稳随机信号；$ \delta (\alpha ) $代表克罗内克符号. $ \delta (\alpha ) $满足

(8)$ \delta (\alpha )=\begin{cases} 1, & \alpha =0;\\0, & \alpha \neq 0.\end{cases} $

式（7）表明谱相关函数只在$ \alpha =0 $和$ \alpha =l\omega $调制频率处取非零值，且能描述二阶循环平稳成分$ {m}_{\text{CS2}}(t) $的载波分布. 在工程现场环境的低信噪比条件下，谱相关容易忽略掉一些关键的微弱调制特征，为此进一步提出谱相干函数：

(9)$ {\gamma }_{x}{(}\alpha ,f{)=}\dfrac{{\rm{S C}}_{x}(\alpha ,f)}{\sqrt{{\rm{S C}}_{x}(0,f){\rm{S C}}_{x}(0,f-\alpha )}}. $

谱相干函数能对不同区域的调制特征进行幅值归一化，进而增强微弱的调制特征. 根据式（7）可推导得到一阶和二阶循环平稳成分的谱相干函数在$ \alpha =0 $和$ \alpha =l\omega $调制频率处的理论值：

(10)$ \left.\begin{split} & {\gamma }_{{{x}_{\text{CS1}}}}\text{(}\alpha ,f\text{)}=\\ &\begin{cases} 1, & \alpha =0;\\ & \\\dfrac{\displaystyle\sum\limits_{{i}_{1},{i}_{2}}{a}_{{{i}_{1}}}a_{{i}_{2}}^{\ast }\delta \left[f-\frac{\left({i}_{1}+{i}_{2}\right)\omega }{2}\right]+{C}+{D}}{{{Q}}_{{{x}_{\text{CS1}}}}\text{(}f\text{)}\times {{Q}}_{{{x}_{\rm{CS1}}}}\text{(}f-\alpha \text{)}}, & \alpha =l\omega . \\ & \end{cases}\\ & {\gamma }_{{{x}_{\text{CS2}}}}\text{(}\alpha ,f\text{)}=\\ & \begin{cases} 1, & \alpha =0;\\ & \\\dfrac{\displaystyle\sum\limits_{{k}_{1},{k}_{2}}{b}_{{{k}_{1}}}b_{{k}_{2}}^{\ast }{P}_{v}\left[f-\frac{{k}_{1}+{k}_{2}}{2}\omega \right]}{{{Q}}_{{{x}_{\text{CS2}}}}\text{(}f\text{)}\times {{Q}}_{{{x}_{\text{CS2}}}}\text{(}f-\alpha \text{)}}, & \alpha =l\omega. \\ & \end{cases}\end{split}\right\} $

式中：i₁ − i₂ = 1；i₁ = − I + l ，$ \cdots $, I； i₂ = − I, $\cdots,$ I − l；$ {Q}_{x_{\mathrm{CS} 1}}$($\lambda $) = $ \displaystyle\sum_i\left|a_i\right|^2 \delta[\lambda\;-\;i \omega]\;+\;$$ u_{n_{\mathrm{s}}} $a₀δ($\lambda $)+$u_{n_{\mathrm{s}}} a_0^* $δ($\lambda $)+$P_{n_{\mathrm{s}}}(\lambda) $；C=$ \;\;\;u_{n_{\mathrm{s}}} a_l $δ(f−lω/2)；D=$\;\;\;u_{n_{\mathrm{s}}} a_{-l}^*$δ(f+lω/2)；k₁−k₂=l；k₁=−K+l, $ \cdots, $K；k₂=−K, $ \cdots, K-l$；$ \mathrm{Q}_{x_{\mathrm{CS} 2}}(\mu)=$$\displaystyle\sum_k\left|b_k\right|^2 P_v[\mu-k \omega]+P_{n_{\mathrm{s}}}(\mu)$.

根据式（10），一阶和二阶循环平稳成分的谱相干函数理论值如图1所示，基于此，可以对推进噪声的循环平稳特性进行总结：

1）$ {\gamma }_{{{x}_{\text{CS1}}}}\text{(}\alpha ,f\text{)} $在$ \alpha =0 $处恒等于1. 在其他区域，$ {\gamma }_{{{x}_{\text{CS1}}}}\text{(}\alpha ,f\text{)} $是关于循环频率$ \alpha $和谱频率$ f $的离散函数. $ {\gamma }_{{{x}_{\text{CS2}}}}\text{(}\alpha ,f\text{)} $在$ \alpha =0 $和$ \alpha =l\omega $调制频率处取非零值，并且在$ \alpha =l\omega $处非零值分布在载波频率范围$ [f_{{\mathrm{l}}}-\left(K-l/2\right)\omega ,\;f_{{\mathrm{u}}}+\left(K-l/2\right)\omega ] $内. 与$ \;\gamma_{x_{\mathrm{CS} 1}}(\alpha, f)$不同，$ {\gamma }_{{{x}_{\text{CS2}}}}\text{(}\alpha ,f\text{)} $是循环频率$ \alpha $的离散函数和谱频率$ f $的连续函数，并且$ {\gamma }_{{{x}_{\text{CS2}}}}\text{(}\alpha ,f\text{)} $在$ \alpha =0 $处恒等于1，在其他区域取值分布在$ [0,1.0] $.

2）通过对比式（10）的上下2个公式可以发现，推进噪声的谱相干函数在$ \alpha =l\omega $调制频率处的取值由$ {\gamma }_{{{x}_{\text{CS2}}}}\text{(}\alpha ,f\text{)} $所主导，这是因为$ {\gamma }_{{{x}_{\text{CS1}}}}\text{(}\alpha ,f\text{)} $只在某些特定的离散点处有非零值，而$ {\gamma }_{{{x}_{\text{CS2}}}}\text{(}\alpha ,f\text{)} $则会在目标信号组分的载波范围内有一系列连续的非零值.

3）相较于谱相关，谱相干显著地增强了$ \alpha =l\omega $处的微弱调制特征，这使得谱相干成为构建加权包络谱的优先选择.

循环平稳分析中一般对谱相干在某一特定频带进行积分得到平方包络谱$ S_{x}^{\text{SES}}\text{(}\alpha \text{)} $或增强包络谱$ S_{x}^{\text{EES}}\text{(}\alpha \text{)} $：

(11)$ S_{x}^{\text{SES}}\text{(}\alpha \text{)=}\left| \displaystyle{\int}_{{f}_{1}}^{{f}_{2}}{\gamma }_{x}\text{(}\alpha ,f\text{)}\right| \text{d}f, $

(12)$ S_{x}^{\text{EES}}\text{(}\alpha \text{)=}\displaystyle{\int}_{{f}_{1}}^{{f}_{2}}\left| {\gamma }_{x}\text{(}\alpha ,f\text{)}\right| \text{d}f. $

式中：$ {f}_{1} $和$ {f}_{2} $分别代表积分频带的下限和上限. 不难发现，积分频带的选择对于增强包络谱的性能有决定性作用，且这一过程具有较强的主观性. 为了降低人为选择积分频带的影响，本研究提出将谱相干的幅值乘以一个权重函数$ w(f) $，再进行积分操作得到加权包络谱，聚焦性地实现轴叶频调制线谱的增强提取：

(13)$ S_{x}^{\text{WES}}(\alpha )=\displaystyle{\int}_{0}^{{F}_{{\mathrm{s}}}/2}w(f)\left| {\gamma }_{x}\left(\alpha ,f\right)\right| \text{d}f. $

式中：$ S_{x}^{\text{WES}}(\alpha ) $代表加权包络谱，$ {F}_{{\mathrm{s}}} $代表信号的采样频率.

为了获得合适的加权函数，须对谱相干在每一谱频率处进行切片操作，并且度量每一个谱频率切片的调制强度进而获得该谱频率对应的加权系数，其核心思想在于：对于调制强度大的谱频率赋予较高的加权系数，对于调制强度小的谱频率赋予较低的加权系数，将所有谱频率对应的加权系数组合起来就得到加权函数$ w(f) $. 本研究提出频域相关峭度这一指标度量谱频率切片的调制强度：

(14)$\left.\begin{split} &{\text{FCK}}_{1}({\alpha }_{\rm{T}},f)=\frac{\displaystyle{\sum}_{n=1}^{N}{\left[\gamma ({\alpha }_{n},f)\cdot \gamma ({\alpha }_{n}-{\alpha }_{\rm{T}},f)\right]}^{2}}{{\left[\displaystyle{\sum}_{n=1}^{N}{\gamma }^{2}({\alpha }_{n},f)\right]}^{2}},\\&{\text{FCK}}_{M}({\alpha }_{\rm{T}},f)=\frac{\displaystyle{\sum}_{n=1}^{N}{\left[\displaystyle{\prod}_{m=0}^{M}\gamma ({\alpha }_{n}-m{\alpha }_{\rm{T}},f)\right]}^{2}}{{\left[\displaystyle{\sum}_{n=1}^{N}{\gamma }^{2}({\alpha }_{n},f)\right]}^{M+1}}.\end{split}\right\}$

式中：$ {\text{FCK}}_{M}({\alpha }_{\rm{T}},f) $代表相关阶次为$ M $时的频域相关峭度；$ {\alpha }_{\rm{T}} $代表调制基频（本研究中特指推进器轴频）；$ \gamma ({\alpha }_{n}-{\alpha }_{\rm{T}},f) $代表$ \gamma ({\alpha }_{n},f) $沿着循环频率轴向右平移$ {\alpha }_{\rm{T}} $后的版本；$ m $代表相关阶次变量，其取值范围是$ 0 \sim M $.

为了基于频域相关峭度构造谱相干加权函数，须先设置相关阶次$ M $的遍历范围$ \left\{{M}_{\min },\cdots ,{M}_{\max }\right\} $和目标频率的遍历范围$ \left[\alpha _{\rm{T}}^{\rm{L}},\alpha _{\rm{T}}^{\rm{U}}\right] $. 在完成参数设置后，提出调制基频$ {\alpha }_{\rm{T}} $自适应确定准则：对于目标频率备选范围$ \left[\alpha _{\rm{T}}^{\rm{L}},\alpha _{\rm{T}}^{\rm{U}}\right] $中的每一循环频率，均计算累加后的频域相关峭度值$ {\sum}_{f}{\text{FCK}}_{M}(\alpha ,f) $.在此基础上，选取最大的$ {\sum}_{f}{\text{FCK}}_{M}(\alpha ,f) $对应的循环频率作为调制基频$ {\alpha }_{\rm{T}} $：

(15)$\left. \begin{split} {\text{FCK}}_{M}(\alpha )=&{\sum}_{f}{\text{FCK}}_{M}(\alpha ,f), \\ {\alpha }_{\rm{T}}=& \arg \max \;{\text{FCK}}_{M}(\alpha ). \end{split}\right\} $

在此基础上，可以得到基于频域相关峭度的谱相干加权函数：

(16)$ w(f)={\text{FCK}}_{M}({\alpha }_{\rm{T}},f). $

遍历加权包络谱的计算流程主要包括3个步骤：计算参数设置、加权函数构造和遍历加权包络谱计算.

1）计算参数设置. 计算参数的设置主要包括相关阶次$ M $遍历范围$ \left\{{M}_{\min },\cdots ,{M}_{\max }\right\} $和目标频率遍历范围$ \left[\alpha _{\rm{T}}^{\rm{L}},\alpha _{\rm{T}}^{\rm{U}}\right] $.

当谱相干中存在较多的干扰线谱成分时，$ M $应取较大的值，如3~5；当谱相干只有少数的调制线谱时，$ M $应取较小的值，如1. 本研究中将相关阶次$ M $遍历范围设为$ \left\{1,\cdots ,5\right\} $，在此基础上，本研究通过遍历操作自适应地确定最优的相关阶次.

对于水下目标识别而言，可根据拟识别的目标类型预估目标频率遍历范围$ \left[\alpha _{\rm{T}}^{\rm{L}},\alpha _{\rm{T}}^{\rm{U}}\right] $. 假设拟提取的目标频率为$ {\alpha }_{\rm{T}} $，那么目标频率遍历范围$ \left[\alpha _{\rm{T}}^{\rm{L}},\alpha _{\rm{T}}^{\rm{U}}\right] $一般取为$ \left[0.8{\alpha }_{\rm{T}},1.2{\alpha }_{\rm{T}}\right] $.

2）加权函数构造. 加权函数计算的核心内容是频域相关峭度的遍历计算，其包含3层的嵌套循环结构. 第1层循环结构是相关阶次$ M $的遍历寻优，其目的是从$ \left\{{M}_{\min },\cdots ,{M}_{\max }\right\} $中自适应地确定最优的相关阶次$ M $；第2层循环结构是谱频率$ f $的遍历计算，其遍历范围是$ \left[0,{F}_{{\mathrm{s}}}/2\right] $，其目的是确定每一个循环频率对应的加权系数；第3层循环结构是循环频率$ \alpha $的遍历计算，其目的是从遍历范围$ \left[\alpha _{\rm{T}}^{\rm{L}},\alpha _{\rm{T}}^{\rm{U}}\right] $中自适应地确定调制基频$ {\alpha }_{\rm{T}} $. 上述嵌套循环结构具体可见表1.

本研究在1.3.2节中构建的频域相关峭度指标倾向于关注少数调制强度较高的谱频率切片，存在遗漏轴叶频调制线谱的风险. 因此，进一步提出基于频域相关峭度的改进指标SFCK（soft FCK）以保留更多有用的谱相干切片，增强抗噪性能：

(17)$ {\text{SFCK}}_{M}(\alpha ,f)=\sqrt{{\text{FCK}}_{M}(\alpha ,f)}. $

式中：$ {\text{SFCK}}_{M}(\alpha ,f) $代表特定的循环频率$ \alpha $和谱频率$ f $下的SFCK. 在以上3层嵌套循环结构中，重复计算对应的SFCK.

在完成循环频率$ \alpha $和谱频率$ f $的SFCK遍历计算后，须确定$ \left\{{M}_{\min },\cdots ,{M}_{\max }\right\} $中所有备选相关阶次$ M $对应的调制基频. 采用1.3.2节所提出的方法：对于目标频率备选范围$ \left[\alpha _{\rm{T}}^{\rm{L}},\alpha _{\rm{T}}^{\rm{U}}\right] $中的每一循环频率，均计算累加后的SFCK值$ {\sum}_{f}{\text{SFCK}}_{M}(\alpha ,f) $，在此基础上，选取最大的$ {\sum}_{f}{\text{SFCK}}_{M}(\alpha ,f) $对应的循环频率作为调制基频$ {\alpha }_{\rm{T}} $：

(18)$ \left.\begin{split} {\text{SFCK}}_{M}(\alpha )=&{\sum}_{f}{\text{SFCK}}_{M}(\alpha ,f), \\{\alpha }_{\text{TWES},M}=&\arg \max \;{\text{SFCK}}_{M}(\alpha ). \end{split}\right\} $

在确定了所有相关阶次$ M $对应的调制基频$ {\alpha }_{\text{TWES},M} $后，对应的加权函数为

(19)$ {w}_{\text{TWES},M}(f)={{\mathrm{SFCK}}}_{M}({\alpha }_{\text{TWES},M},f). $

在此基础上，可以获取所有相关阶次$ M $对应的遍历加权包络谱：

(20)$ S_{x}^{\text{TWES},M}\left(\alpha \right)=\displaystyle{\int}_{0}^{{F}_{{\mathrm{s}}}/2}{w}_{\text{TWES},M}(f)\left| {\gamma }_{x}\left(\alpha ,f\right)\right| \text{d}f. $

进一步地，对以上遍历加权包络谱进行归一化处理：

(21)$ S_{x-N}^{\text{TWES},M}\left(\alpha \right)=\dfrac{S_{x}^{\text{TWES},M}\left(\alpha \right)-\min \left[\mathrm{S}_{x}^{{{\mathrm{TWES}}},M}\left(\alpha \right)\right]}{\max \left[S_{x}^{\text{TWES},M}\left(\alpha \right)\right]-\min \left[S_{x}^{\text{TWES},M}\left(\alpha \right)\right]}. $

3）遍历加权包络谱计算. 在获取所有相关阶次$ M $对应的归一化遍历加权包络谱$ S_{x-N}^{\text{TWES},M}\left(\alpha \right) $后，用谐波强度来度量其解调效果优劣，谐波强度指标的定义为

(22)$ H(M)=\left(\prod_{l=1}^L p_{M, l \omega}\right)^{{1}/{L}}. $

式中：$ {p}_{M,l\omega } $为$ S_{x-N}^{\text{TWES},M}\left(\alpha \right) $在特征频率$ \alpha =l\omega $处的取值，$ L $为谐波结构中的最高阶次. 谐波强度$ H\left(M\right) $取得最大值时对应的$ M $为最优相关阶次$ {M}_{\text{TWES-opt}} $，其对应的加权函数为$ {w}_{\text{TWES},{{M}_{\text{TWES-opt}}}}\left(f\right) $.

将$ {w}_{\text{TWES},{{M}_{\text{TWES-opt}}}}\left(f\right) $转化为联合加权形式，利用信息阈值（$ T=\overline{w}(f)+\sigma (w(f)) $，$ \overline{w}(f) $和$ \sigma (w(f) $，分别为加权函数$ w(f) $的平均值和标准差）滤除噪声干扰组分所对应的谱频率加权系数，进一步提升解调抗噪性：

(23)$ \begin{split} &{w}_{\text{TWES}}(f)=\\&\begin{cases} 0, & {w}_{\text{TWES},{{M}_{\text{TWES-opt}}}}\left(f\right) \lt T;\\{w}_{\text{TWES},{{M}_{\text{TWES-opt}}}}\left(f\right), & {w}_{\text{TWES},{{M}_{\text{TWES-opt}}}}\left(f\right)\geqslant T.\end{cases}\end{split} $

最终计算得到最优的遍历加权包络谱：

(24)$ S_{x}^{\text{TWES}}\left(\alpha \right)=\displaystyle{\int}_{0}^{{F}_{\rm{s}}/2}{w}_{\text{TWES}}(f)\left| {\gamma }_{x}\left(\alpha ,f\right)\right| \text{d}f. $

遍历加权包络谱的计算流程如算法1所示.

算法1

开始

% 计算参数设置 %

设置相关阶次$ M $遍历范围$ \left\{{M}_{\min },\cdots ,{M}_{\max }\right\} $和目标频率遍历范围$ \left[\alpha _{\rm{T}}^{\rm{L}},\alpha _{\rm{T}}^{\rm{U}}\right] $

% 遍历加权包络谱计算 %

For $ M={M}_{\min }\colon {M}_{\max } $

For $ f=0\colon {F}_{\rm{s}}/2 $

For $ \alpha ={\alpha }_{n1}\colon {\alpha }_{n2} $

计算$ {\text{SFCK}}_{M}\left(\alpha ,f\right) $

End

End

计算不同$ \alpha $对应的$ {\sum}_{f}{\text{SFCK}}_{M}(\alpha ,f) $

找到$ {\sum}_{f}{\text{SFCK}}_{M}(\alpha ,f) $最大值, 确定$ {\alpha }_{\text{TWES},M} $

End

计算所有$ M $对应的$ {w}_{\text{TWES},M}\left(f\right) $

构建所有$ M $对应的$ S_{x}^{\text{TWES},M}\left(\alpha \right) $

获取所有$ M $对应的归一化$ S_{x-N}^{\text{TWES},M}\left(\alpha \right) $

计算所有$ S_{x-N}^{\text{TWES},M}\left(\alpha \right) $的谐波强度

找到谐波强度最大值对应的相关阶次$ {M}_{\text{TWES-opt}} $

将w_TWES、M_TWES-opt(f)转化为联合加权形式w_TWES(f)

得到最优的遍历加权包络谱$ S_{x}^{\text{TWES}}\left(\alpha \right) $

结束

基于信号解调的推进器被动声学探测可形式化为一个二元假设检验，原假设$ {H}_{0} $表示未探测水下推进器（包络谱中只含背景噪声），备择假设$ {H}_{1} $表示探测到水下推进器（包络谱中检测到轴叶频调制线谱）：

(25)$ \left.\begin{aligned}{H}_{0}\colon {{\boldsymbol{d}}}(n)&={\boldsymbol{n}}(t);\\{H}_{1}\colon {{\boldsymbol{d}}}(n)&={{\boldsymbol{m}}}_{\text{CS}}(t)+{\boldsymbol{n}}(t).\end{aligned} \right\}$

式中：$ {{\boldsymbol{d}}}(n)\in {\bf{R}}^{N} $为待分析的包络谱，$ {\boldsymbol{n}}(t) $代表背景噪声，$ {{\boldsymbol{m}}}_{\text{CS}}(t) $代表推进噪声. 为了通过以上假设检验实现推进器被动声学探测，须设计对应的假设检验指标.

提出利用自适应阈值统计曲线来实现假设检验指标的构建：对于包络谱$ S_{x}\left(\alpha \right) $，其统计阈值$ {\lambda }_{1-{{p}_{\rm{fa}}}} $定义为有$ 1-{p}_{\rm{fa}} $比例的取值在其以下的界限值. 包络谱的整个循环频率范围可以划分为一系列窗口$ {W}_{k} $（其重叠比例为$ R $）. 对于每个窗口，计算其统计阈值$ \lambda _{1-{p}_{\rm{fa}}}^{k} $. 在此基础上，将所有窗口阈值连接起来并用三次样条曲线进行插值，得到自适应阈值曲线$ {\lambda }_{1-{{p}_{\rm{fa}}}}(\alpha ) $，其点数与$ {{\boldsymbol{S}}}_{x}\left(\alpha \right) $完全一致. 自适应阈值曲线的构造过程如图2所示. 其中，$ {{A}}_{\text{norm}} $为归一化幅值，$ \alpha $为循环频率，n₀为点数.

基于此，可以利用特征峰噪比（characteristic peak-to-noise peak ratio，CNR）整体性地评价调制线谱相对于背景噪声和干扰线谱的突出程度：

(26)$ \left.\begin{split}\text{CNR=}&{\text{CPR}}/{\text{NPR}},\\ \text{CPR}=&{\left(\prod\limits_{k=1}^{K}\dfrac{{S}_{x}(k\omega )}{{\lambda }_{1-{{p}_{\rm{fa}}}}(k\omega )}\right)}^{{1}/{K}}, \\ \text{NPR}=&{\left(\prod\limits_{l=1}^{{L}}\dfrac{{S}_{x}({\alpha }_{{{n}_{l}}})}{{\lambda }_{1-{{p}_{\rm{fa}}}}({\alpha }_{{{n}_{l}}})}\right)}^{{1}/{L}}.\end{split}\right\} $

式中：$ \text{CPR} $代表特征峰值比（characteristic peak ratio），其可以评价调制线谱的凸显程度；$ \text{NPR} $代表噪声峰值比（noise peak ratio）；$ {\alpha }_{{{n}_{l}}} $代表包络谱$ {S}_{x}\left(\alpha \right) $中干扰线谱对应的循环频率. 本研究关注轴叶频调制线谱，因此参与CNR指标构建的调制线谱为轴频及其谐波（所关注的最高谐波为叶频，即K为叶片数）. 在此基础上，可利用峰噪比实现推进器被动声学探测的二元假设检验：

(27)$\left. \begin{split} &接受{H}_{0}并拒绝{H}_{1}\colon \text{CNR}\leq 1;\\&接受{H}_{1}并拒绝{H}_{0}\colon \text{CNR} \gt 1. \end{split}\right\}$

利用仿真信号验证遍历加权包络谱的有效性和优越性，所用到的仿真信号表达式如下：

(28)$ \left.\begin{split} x(t)=& m_{\mathrm{CS} 1}(t)+m_{\mathrm{CS} 1}(t) \cdot v(t)+n_{\mathrm{CS} 1}(t)+ \\& n_{\mathrm{CS} 1}(t) \cdot v_N(t)+n_{\mathrm{p}}(t)+n_{\mathrm{s}}(t), \\ m_{\mathrm{CS} 1}(t)=&\sum_i A_i \cos \;\left(2 {\text{π}} i f_{\mathrm{SF}} t+\phi_i\right)+ \\& \sum_j A_j \cos\; \left(2 {\text{π}} j f_{\mathrm{BPF}} t+\phi_j\right) , \\ n_{\mathrm{CS} 1}(t)=&\sum_k B_k \cos \;\left(2 {\text{π}} k \alpha_N t+\varphi_k\right) .\end{split} \right\}$

式中：$ {m}_{\text{CS1}}(t) $和$ {m}_{\text{CS1}}(t)\cdot v(t) $分别代表一阶和二阶循环平稳成分，其中$ v(t) $为二阶循环平稳成分的载波信号（首先生成指定功率的高斯白噪声，然后进行带通滤波得到指定频率范围的载波信号）；$ {n}_{\text{CS1}}(t) $和$ {n}_{\text{CS1}}(t)\cdot {v}_{N}(t) $分别代表一阶和二阶循环平稳噪声，其中$ {v}_{N}(t) $为二阶循环平稳噪声的载波信号（生成方式同上）；$ {n}_{\text{p}}(t) $代表冲击噪声（首先生成若干个指定功率的高斯白噪声，信号时长控制为0.1 s，然后将其随机添加至仿真信号中）；$ {n}_{\text{s}}(t) $代表平稳噪声（通过生成指定功率的高斯白噪声得到）.

式（28）中相关参数的含义和取值详见表1，其中轴频$ {f}_{\text{SF}} $设置为21 Hz，叶片数设置为7，即叶频$ {f}_{\text{BPF}} $为147 Hz. 为了充分验证遍历加权包络谱的优越性，将其与经典及前沿的信号解调方法进行对比，具体包括3类对比方法：基于快速峭度谱的增强包络谱（fast kurtogram-based EES, FK EES）^[14]、增强包络谱（enhanced envelope spectrum, EES）^[22]和改进包络谱（improved envelope spectrum, IES）^[17]. 所有对比方法以及本研究方法的结果均被归一化至$ [0,1.0] $，其提取到的调制基频及其谐波用红圈标出. 对比方法的具体介绍如下：1）FK EES：通过快速峭度谱计算不同分解频层中每个频带信号的峭度，选择最佳解调频带（最佳积分频带），并将谱相干函数在该频带内进行积分操作，得到FK EES；2）EES：将谱相干函数在$ [0,{F}_{\rm{s}}/2] $的频带范围内（即全频带范围）进行积分操作，得到EES；3）IES：通过IESFOgram选择最佳解调频带（最佳积分频带），并将谱相干函数在该频带内进行积分操作，得到IES.

当信噪比为−10 dB时，仿真信号的分析结果如图3、4所示. 解调频带选择结果如图3(a)~(d)所示. 其中，$ {n}_{\text{d}} $为分解频层，f为谱频率，SCoh表示谱相干，$ {{A}}_{\text{w}} $为加权系数，SK为谱峭度，DF为诊断指标. 图3(a)和3(c)中的纵坐标n_d代表了滤波器组对信号频谱的划分精细程度，层数越高，对应的频带划分越细，带宽越窄. 由图3(a)可知，FK选择的最佳解调频带为8 333.3~10 000.0 Hz，这与CS2噪声的载波范围以及脉冲噪声的频率范围一致. 不难发现，FK未能捕捉到目标组分的载波频带. 图3(b)验证了谱相干函数能在特征频率$ l\omega $处捕捉到目标组分的载波分布（如红框所示），此外谱相干函数可提取到噪声干扰成分及其载波分布（如白框所示），即谱相干会将所有的二阶循环平稳信号提取出来. 如图3(c)所示，IESFOgram选择的最佳解调频带是14 375~15 000 Hz，该频带落在目标组分的第2个载波频带内，从而确保了后续IES能有效提取调制频率. 遍历加权包络谱的加权函数（即联合加权形式下的$ {w}_{\text{TWES},{{M}_{\text{TWES-opt}}}}\left(f\right) $）如图3(d)所示，该加权函数捕捉到目标组分位于3 500~6 000 Hz和14 000~16 500 Hz的2个载波频带，对于目标组分的提取能力优于其他对比方法.

信号解调分析结果如图4(a) ~(d)所示. 由图4(a)可知，FK EES提取到噪声调制频率而未能提取到目标调制频率，$ {\text{CNR}}_{\text{FK EES}}=0.61 \lt 1.00 $也证实FK失效. 如图4(b)所示，EES可以提取到目标调制频率，$ {\text{CNR}}_{\text{EES}}=1.63 \gt 1.00 $也支持这一发现，但噪声调制频率仍在EES中占据主导地位. 如图4(c)所示，IES不仅成功提取到目标调制频率，同时抑制了噪声调制频率，$ {\text{CNR}}_{\text{IES}}=2.99 $也明显高于前2种方法. 如图4(d)，TWES不仅聚焦性地提取到目标调制频率，并有效地抑制了循环平稳噪声及其他干扰成分，而$ {\text{CNR}}_{\text{TWES}} = 3.46 \gt {\text{CNR}}_{\text{IES}} \gt {\text{CNR}}_{\text{EES}} \gt {\text{CNR}}_{\text{FK EES}} $进一步证实了TWES相较于对比方法展示出更加优异的解调性能.

为了进一步验证遍历加权包络谱的抗噪性能，对−20~0 dB（以1 dB为间隔）下的仿真信号进行分析处理. 对于每一信噪比，生成1 000组仿真信号进行分析，即有 $ {N}_{{\mathrm{a}}}=1\;000 $. 对于TWES及各对比方法，计算特定信噪比下所有仿真信号的峰噪比平均值来评价其解调性能：

(29)$ \text{CNR}=\displaystyle\sum\limits_{n=1}^{1\;000}{\text{CNR}}_{n}/{N}_{{\mathrm{a}}}. $

式中：$ {\text{CNR}}_{n} $为待评价的解调方法在某特定信噪比下对第$ n $个仿真信号分析后的峰噪比. TWES及各对比方法的蒙特卡洛模拟结果如图5所示，其中，SNR为信噪比. 可以看出，FK解调性能不佳，$ {\text{CNR}}_{\text{FK EES}} $在−20~0 dB的信噪比范围内均保持着较低值；IESFOgram的解调性能相对较差，$ {\text{CNR}}_{\text{IES}} $在信噪比＜−11 dB时已小于1.5，处于较低水平；谱相干函数的解调性能优于前2种方法，$ {\text{CNR}}_{\text{EES}} $在信噪比≥−14 dB时保持大于1.5，证明其在该信噪比范围内能有效检测到目标调制频率，这主要得益于其在全频段内寻找循环平稳特征的能力；本研究所提出的TWES呈现最优的解调性能，$ {\text{CNR}}_{\text{TWES}} $在−20~0 dB的信噪比范围内均高于各对比方法，在信噪比≥−16 dB时保持大于1.5，证明该方法具有较好的抗噪性.

为了验证遍历加权包络谱的解调提取能力，对某泵喷推进器的辐射噪声数据进行处理分析. 实验在上海船研所大型水洞中进行^[23]. 实验对象为某一自研的前置导叶泵喷推进器，设计参数由某一航行体推进需求确定，并进行缩比处理得到，如表2所示. 单组水声信号的测量时长为10 s，采样频率为20 960 Hz，由于额定转速为1 260 r/min，对应地，其调制基频为21 Hz，在所有方法的分析结果中均将提取到的调制基频及其谐波用红圈标出.

泵喷推进器辐射噪声的解调分析结果如图6(a)~(d)所示. 由图6(a)可知，FK EES只提取到2倍轴频. 由图6(b)可知，EES有效地提取到了以轴频及其谐波为代表的目标频率. 如图6(c)所示，IES与FK EES结果比较接近，同样也只提取到2倍轴频. 如图6(d)所示，遍历加权包络谱提取到轴频及其谐波，对于泵喷推进器轴叶频调制线谱的提取能力优于其他对比方法. 峰噪比$ {\text{CNR}}_{\text{TWES}}=1.52 $>$ {\text{CNR}}_{\text{EES}} $>$ {\text{CNR}}_{\text{IES}} $>$ {\text{CNR}}_{\text{FK EES}} $也反映TWES对泵喷推进器调制线谱的检测性能最好.

为了进一步验证遍历加权包络谱的解调提取能力，对湖试环境下的某自研水下航行器辐射噪声数据进行处理分析. 实验在浙江某湖试试验场进行，水下航行器的湖试实验示意图如图7所示，主要包括水下航行器（含螺旋桨推进器）、水听器、固定装置、信号采集系统、测试平台（即测试船）. 水下航行器运行过程中围绕着某一固定中心做圆周运动，并用GPS接收器获取其位置信息. 水声信号的测量采用舟山海宇海洋科技公司的LNH-20水听器，其具体参数如表3所示. 数据采集仪型号为NI-9250，具体参数如表4所示，单组水声信号的采样频率设置为51 200 Hz.

实验过程包括以下步骤：1）测试前准备：测试船和水下航行器到达合适测量点后锚定测试船，在测试船上布置测试系统，水听器数量为2只，分别固定在水深1、2 m的位置；记录测试船GPS信息；水下航行器靠近测试船并校准GPS位置； 2）背景噪声测量：用水听器测量背景噪声，每组测试20 s，测试3组；3）水下噪声测量：启动航行器，转速控制在768 r/min，确保安全后操纵其入水，进行水下转圈航行（围绕某一中心绕圈），记录时间、转速、GPS信息，开始水下航行后进行测量，共测试3组，每组20 s；切换转速1150 r/min，进行水下噪声测量，共测试3组，每组20 s；4）测试结束：返航，回收航行器，核对下水时间点，校对数据.

为了全面验证所提方法的有效性，分别对转速为768、1 150 r/min时的航行器水下噪声进行分析，对应地，其调制基频分别为$ {f}_{\text{SF1}}\text{=12.8 Hz} $和$ {f}_{\text{SF2}}\text{=19.2 Hz} $，在所有方法的分析结果中均将提取到的调制基频及其谐波用红圈标出.

当转速为768 r/min时，解调分析结果如图8(a)~(d)所示. 由图8(a)可知，FK EES未能提取到以轴频及其谐波为代表的目标调制频率. 如图8(b)可知，EES中轴频及其谐波也几乎不可辨识. 如图8(c)所示，与FK EES一致，IES同样无法提取到目标调制频率. 如图8(d)所示，遍历加权包络谱提取到轴频及其谐波，对于推进器调制线谱的提取能力优于其他对比方法. 峰噪比结果$ {\text{CNR}}_{\text{TWES}}=2.47 $>1.00>$ {\text{CNR}}_{\text{IES}} $、$ {\text{CNR}}_{\text{EES}} $、$ {\text{CNR}}_{\text{FK EES}} $，反映只有TWES成功检测到了推进器.

对转速为768 r/min时3次测量的水声信号（测量距离依次递增）进行解调分析，对比不同方法对于轴叶频调制线谱的检测结果，其结果如表5所示. TWES均可稳定检测到轴叶频调制线谱，其对应的$ {\text{CNR}}_{\text{TWES}} $分别为2.47、1.85和1.22，而$ {\text{CNR}}_{\text{FK EES}} $、$ {\text{CNR}}_{\text{EES}} $、$ {\text{CNR}}_{\text{IES}} $始终小于1.00，证明这3种对比方法在此微弱信号场景下已失效，也进一步通过实验数据验证了本研究所提出的TWES方法对微弱信号解调检测的有效性和优越性.

转速为1 150 r/min时的水声信号的解调分析结果如图9所示. 此外还引入了电磁激励噪声（即二阶循环平稳噪声）以进一步验证TWES对于推进器调制线谱的聚焦性提取能力. 由图9(a)、(b)可知，FK EES和EES完全被100 Hz处电磁激励噪声所主导，没有提取到以轴频及其谐波为代表的目标频率. 如图9(c)所示，IES同样没有提取到目标频率. 如图9(d)所示，遍历加权包络谱提取到了以轴频及其谐波为代表的目标频率. 同样地，本研究通过峰噪比结果$ {\text{CNR}}_{\text{TWES}}=1.46 $>1.00>$ {\text{CNR}}_{\text{IES}} $、$ {\text{CNR}}_{\text{EES}} $、$ {\text{CNR}}_{\text{FK EES}} $进一步证实，TWES成功检测到推进器.

本研究提出面向水下推进器调制线谱提取的遍历加权包络谱TWES，其核心思想是对谱相干函数进行加权操作，在目标组分的载波范围内赋予较高的加权系数，而在其他区域赋予较低的加权系数. 在此基础上，本研究提出峰噪比CNR这一指标，用于实现解调方法的性能量化评估和基于信号解调的水下推进器自动检测.

仿真信号、泵喷推进器水洞实验及水下航行器湖试数据的分析结果表明，当水声信号在信噪比不低于–16 dB时，或含有一定电磁激励噪声（在先验频率范围与干扰频率不重合的前提下）的条件下，本研究所提出的遍历加权包络谱（TWES）方法能有效提取以轴频为代表的调制线谱，其解调性能优于快速峭度谱、增强包络谱及改进包络谱等现有方法.

遍历加权包络谱（TWES）方法依赖于准确的先验频率范围且假设先验频率与干扰噪声不重合，因此在实际复杂工况下的方法鲁棒性待验证. 未来研究一方面可探索自适应加权策略，将该方法拓展至多源干扰共存或更低信噪比场景，另一方面可与小样本条件下的深度学习理论相结合，发展基于调制线谱提取的推进器智能识别方法乃至桨叶数智能分类模型.