基于Rao检测器的舰船轴频电场滑动门限检测方法

喻鹏，程锦房，张伽伟，姜润翔

(1.海军工程大学 兵器工程学院，湖北 武汉 430033；2.海军工程大学 电气工程学院，湖北 武汉 430033)

0 引言

舰船在航行过程中，会在海水中产生轴频电场，其基频与螺旋桨的转速基本一致(大约处于1～7 Hz之间)[1]，具有较明显的线谱特征，可作为对目标探测和识别的重要物理场信息[2-4]。

由于轴频电场在海水中衰减较快，且目前轴频电场消减技术也在逐步发展[5-6]，如何在更低信噪比下，实现更可靠的检测，成为轴频电场应用研究的重点。文献[7]首先采用一种基于小波包分解的微弱轴频电场检测算法。文献[8]基于海洋环境电场和轴频电场具有不同自回归(AR)模型参数的特点进行轴频电场检测。文献[9-10]采用经验模态分解(EMD)、集合经验模态分解(EEMD)和改进功率谱熵的轴频电场检测方法。上述检测方法计算量较大。文献[11]采用广义似然比检验(GLRT)法对船舶的静电场进行探测和跟踪，但仅利用水平电偶极子进行了仿真计算，后期未进一步研究。文献[12]基于舰船通过时轴频电场特征频段内能量增加特性，提出一种滑动功率谱检测方法，较好地解决了检测概率与虚警问题，计算量小，目前应用较多。

上述研究多重点关注信号滤波处理算法，然后基于能量变化实现目标检测，较少考虑信号模型和噪声模型本身的特性，即没有充分利用目标的信息。本文基于统计信号处理理论，提出一种基于Rao检测器[13]的轴频电场滑动门限检测方法(简称滑动门限Rao检测方法)。首先建立轴频电场信号的检测模型，然后利用混合高斯模型(GMM)对环境电场噪声进行建模，并采用期望最大化(EM)迭代算法实现GMM参数的实时估计，在上述信号模型和噪声模型估计基础上计算Rao检测统计量，并采用前一时刻的Rao检测统计量均值作为滑动门限实现轴频电场检测。

为验证滑动门限Rao检测方法的有效性，本文首先利用仿真计算对比Rao检测方法与基于能量的检测方法的性能差异，然后利用不同船舶的轴频电场数据对比了滑动门限Rao检测方法与文献[12]提出的滑动功率谱检测方法的性能。

1 基本原理

在高斯背景噪声中，当没有待检测信号的任何先验信息时，能量检测器是一种最优检测器[14]，但是在非高斯背景噪声中，能量检测器的检测性能会随着噪声非高斯特性的增强而逐渐降低[15]。在实际环境现场，噪声往往表现出非高斯特性，如近岸海洋环境电场容易受到自然及人为因素的影响，往往表现出非高斯特性。在这种情况下，采用GLRT法将实现较好的效果，由于GLRT法需要估计出H0(目标不存在)和H1(目标存在)情况下未知参数的最大似然估计(MLE)，而Rao检测方法是GLRT法的改进算法，它只需要在H0情况下进行参数估计，大大降低了算法的复杂度和运算量，在通讯领域被广泛采用[16-17]。

为实现基于Rao检测器的轴频电场检测，需要建立信号模型和噪声模型。

1.1 信号模型

轴频电场信号属于零均值、未知频率(已知频率范围)和幅值的周期信号，并且还存在一定数量的谐波分量，基于傅里叶级数可以用如(1)式线性模型来表示该轴频电场信号：

(1)

式中：s[n]为轴频电场信号(不包含噪声)；M为谐波分量的个数，根据谐波能量分布特性，一般取M=3；ak和bk分别为相应谐波的幅值；f0为轴频电场信号的基频，f0∈[fmin,fmax],fmin、fmax分别为轴频电场基频的最小值和最大值；n为离散时刻值；N为信号总时长。

将信号模型转换为矩阵形式：

(2)

(2)式可简写为

S=H(f0)θ，

(3)

式中：H(f0)和θ分别代表(2)式右侧两项，θ=[a1,a2,a3,b1,b2,b3]T.

1.2 噪声模型

海洋环境电场长期统计特性近似呈高斯分布规律[18]，但是对于轴频电场实时性检测应用背景下(0.5～30 Hz频率范围)，由于检测数据时长短，且受到外界因素的影响，环境噪声的高斯特性很难满足。这里根据文献[19]的理论，对噪声的非高斯特性进行判别。对于零均值概率密度函数(PDF)的信号w0[n],n=1,2,…,N，其非高斯程度通常由相对高斯PDF的峰态γ表示，如(4)式所示：

(4)

式中：E(·)表示求平均。

高斯噪声对应γ=0，非高斯噪声对应γ偏离0. 选取实测环境电场数据进行评估，选择第3节中水面浮动平台测量背景下目标不存在时刻对应的环境电场数据，每次选择50 s数据长度，共11组数据，以基本覆盖全天的环境电场特性。这11组环境电场噪声的峰态γ分别为0.272 6，0.280 2，0.207 5，-0.046 4，0.064 2，-0.249 2，0.217 0，-0.314 0，0.600 9，-0.543 6，-0.357 7. 可以看出峰态γ偏离0的程度较大，说明在轴频电场范围内，该水面浮动平台下海洋环境电场噪声呈现出一定的非高斯特性，可以采用2阶0均值GMM进行拟合。

(5)

混合高斯分布的4阶矩[20]为

(6)

(5)式代入(4)式，可以得到混合高斯噪声的峰态值γ为

(7)

1.3 检测方法实现

在1.1节信号模型和1.2节噪声模型分析的基础上，完整的数据模型为

x[n]=s[n]+w[n],

(8)

式中：x[n]表示轴频电场信号(包含噪声)。

定义如下检测统计量[19]：

(9)

(10)

(11)

由此可求得y[n]：

y[n]=

(12)

(10)式代入(11)式，得

(13)

由(13)式可以进一步看出i(A)能够用[g(w)]2的数学期望值代替，则在无目标条件下可以用高斯化滤波结果的数学平均来近似i(A)，从而大幅度降低计算复杂度[21]。

另外，当基频f0≫1/N时(可通过调节数据长度实现)，H(f0)的列向量近似正交，则HT(f0)H(f0)≈(N/2)I，I为单位矩阵，则(9)式可化简为

(14)

由于

(15)

则(14)式可化简为

TRao[y]=

(16)

1.4 具体流程

理论上，Rao检测方法的步骤为：首先在H0情况下估计出GMM的PDF，然后求解检测统计量TRao，最后依据阈值对检测结果进行判断。值得注意的是，在实际水下目标检测应用领域，为降低虚警概率，往往并不会依据虚警率指标来设定一个固定门限，而是采用滑动门限的方式进行目标检测，从而在降低虚警概率的基础上提升检测概率。另外，利用GMM对环境噪声建模时，由于环境噪声特性也在不断变化，需要不断更新GMM的参数。

对于GMM的参数估计，采用EM迭代算法计算[22]，该方法计算量小，对于本文的2阶0均值GMM参数估计只需数十次迭代(计算量很小)，即可实现较好的估计效果，其基本步骤如下：定义噪声数据中第p个样本值的指示向量为zp,q，当该样本来自第q个高斯源时zp,q=1，否则zp,q=0，其中高斯源总数Mg=2. 根据(17)式计算zp,q的期望值〈zp,q〉，即样本值来自第q个高斯源的后验概率密度：

(17)

(18)

(19)

值得注意的是，在估计GMM参数时，所选择的样本应该为本时刻之前Δt时刻某固定时长的环境电场数据，这样就可以实现自适应GMM参数估计，即(10)式具有自适应限幅器的功能。具体实施流程如图1所示。图1中Δn为时延，m为选取的数据长度。

图1 滑动门限Rao检测方法流程图

实测信号首先进行带通滤波(0.5～30 Hz)并估计GMM参数，基于本时刻之前Δt时刻的GMM参数并利用(16)式计算检测统计量TRao；同时计算本时刻n之前n-Δn-m+1时刻至n-Δn时刻的平均TRao值，之后乘以倍数κ作为本时刻的门限值，即阈值Γ[n]为

(20)

值得注意的是，由于舰船通过时TRao会大幅增加，从而导致在舰船通过Δn时长后，阈值Γ会大于实际环境对应的阈值，从而影响目标通过后这段时刻的检测效果。因此在具体计算阈值时，当确认检测到目标后，会将确认检测到目标时刻附近对应的TRao值作为阈值Γ在本阶段的上界，从而降低这段时间的漏检概率，该部分阈值设定并未在图1中进行描述，但是在后文实船检测时可以看到其效果。

2 仿真分析

为评估Rao检测方法相比于能量检测的优势，进行如下仿真计算。设定信号s的基频f0为3 Hz，并带有一定的谐波分量。设定每段数据长度N为500个(模拟采样率100 Hz情况下5 s时长的检测数据长度)，仿真次数共1 000次。注意这里仿真计算仅评估Rao检测方法的效果，而与第1节的滑动门限无关，滑动门限属于实际应用而非仿真计算。

由于仿真分析时确切知道信号和噪声存在时的PDF，可利用Neyman-Pearson(NP)准则求出最佳检测效果，即NP上界作为参考，其检测统计量为

(21)

式中：x为待检测信号。

另外，为对比分析，这里也求出基于能量检测的检测效果，其检测统计量为

(22)

不同信噪比情况下两种检测方法的受试者工作特征曲线(ROC)对比如图2所示。由图2可以看出：在SNR=-10 dB时，Rao检测方法趋近NP上界，即趋近最佳检测效果；当信噪比逐渐降低时，Rao检测方法和能量检测方法的检测效果均下降，但是Rao检测方法的性能始终优于能量检测方法。

图2 Rao检测方法与能量检测方法检测效果对比

3 实船试验

3.1 试验情况介绍

2019年11月在大连港口附近进行试验，测量来往不同类型客货轮的轴频电场信号，验证滑动门限Rao检测方法的效果。测量电极是安装在水面浮动平台底部(在水面以下约30 cm位置)，如图3所示。

图3 电场测量装置试验现场

3.2 检测效果分析

为分析滑动门限Rao检测方法相比于滑动功率谱检测方法的优势，选择本次试验中2组信噪比较差的轴频电场信号进行检验。

滑动门限Rao检测方法参数设置为：对原始数据进行0.5～30 Hz带通滤波；GMM参数估计时延Δt=30 s；轴频电场信号基频f0取值为1～7 Hz之间的离散点，其频率间隔为0.1 Hz；Rao统计量计算时延Δn=30 s；计算阈值的数据长度m=30fs，即时长为30 s的数据；阈值计算倍数κ=3；检测数据长度N=5fs；判定检测到的持续时间为3 s.

在对比滑动门限Rao检测方法与滑动功率谱检测方法时，首先对比二者实时检测统计量的差异，暂不考虑阈值对检测效果的影响，即对二者的实时检测统计量进行归一化处理，比较两种检测方法在目标通过前、后的实时检测统计量变化。归一化基本原理为：以各检测方法的检测统计量均值为参照，分别对两种检测方法的检测量进行归一化，之后以两种检测方法中检测统计量最大值为参考，进行统一归一化处理。

图4(a)所示为“连港6号”拖船通过时的轴频电场水平分量(Ex和Ey)，其中在60 s时距离测量点最近(正横距离约100 m)。由于拖船吨位小，且仅通过牺牲阳极进行阴极保护，所以电场源强度较小，在图4中不能从时域看到明显通过特性(本文中电场时域图均进行了0.5～30 Hz带通滤波)。

图4 小型拖船轴频电场及检测效果

对Ex分量进行检测，效果如图4(b)所示。由图4(b)可以看出：在整个通过特性中，滑动功率谱检测方法的实时检测统计量始终处于波动状态，在60 s时刻虽然实时检测统计量有所上升，但并没有大幅超过目标未到达时刻对应的平均检测统计量；相比之下，本文提出的滑动门限Rao检测方法在60 s时实时检测统计量明显高于目标未达到时刻对应的平均检测统计量，具有更好的检测效果。

之后设定滑动功率谱检测方法阈值的倍数与滑动门限Rao检测方法的阈值参数一致，对一段船舶通过特性进行检测。图5所示为2艘小型拖船(45 s和70 s时刻正横距离100 m通过)、“万荣海”大型客轮(170 s时刻正横距离80 m通过)接连通过时对应的轴频电场时域图。由图5(a)可以看出，仅当大型客轮通过时，可从时域上看到明显的通过特性曲线，其他小型船舶均不能从时域上看到明显的轴频电场信号。

图5 多艘船舶连续通过时的轴频电场及检测效果

图5(b)所示为滑动门限Rao检测方法的检测效果，纵坐标为归一化检测统计量，可以看出小型拖船及大型客轮均被检测到；图5(c)所示为滑动功率谱检测方法，当其滑动门限阈值倍数κ=3时，将产生较大的虚警，主要是因为在图5(c)70～150 s和210～240 s时实时检测统计量存在较大波动，而此时并未出现目标。因此，这将滑动功率谱检测方法的门限阈值倍数κ设置为4，得到图5(c)所示的检测效果，可以看出漏检一艘拖船。

综合图4和图5的检测效果，可以看出滑动门限Rao检测方法相比于滑动功率谱检测方法的优势主要是在目标未到达时刻能够抑制噪声干扰，此时实时检测统计量处于较低水平；当目标到达时，实时检测统计量迅速增加。其主要原因是：所提方法实时估计了环境噪声的参数，即实时估计GMM的PDF，确保(10)式能够有效抑制非高斯尖峰噪声干扰；而当目标达到时，采用信号模型检测，实现最大化检测效果。

4 结论

本文针对非高斯背景噪声和低信噪比情况下船舶轴频电场实时检测问题，提出一种滑动门限Rao检测方法。首先利用仿真计算证明了Rao检测方法相比于能量检测方法在同样虚警率情况下具有更高的检测概率；然后利用水面浮动平台背景下实测轴频电场数据，对比滑动门限Rao检测方法与现有滑动门限功率谱检测方法，结果表明：针对本文实测数据，滑动门限Rao检测方法能够有效抑制噪声干扰，检测效果优于现有滑动门限功率谱检测方法。