1.   引言

辐射源的位置可以通过测量多个传感器接收信号的到达时差(Time Difference Of Arrival, TDOA)来确定^[1]。经典的基于相关的TDOA估计器可以使用几种广义互相关技术中的一种来构造^[2]。在许多传播环境中，每个传感器的接收信号中除了包含直达辐射源信号外还包括一个或多个辐射源信号的多径副本信号。由于受到多径干扰的影响，传感器间的互相关呈现多个峰值。互相关峰值之间的延迟分辨能力、相对幅度和相位偏移决定了多径副本信号对直达辐射源信号TDOA估计精度的影响。多径环境下的TDOA估计技术有很多种，在辐射源信号已知的情况下匹配滤波提供了估计多径信号TDOA的近似最优方法^[3]。但是这种方法中多径TDOA的分辨能力受到辐射源信号相关时间的限制，当多径TDOA之间的时间间隔小于辐射源信号的相关时间时就必须应用高分辨率的多径TDOA估计算法。

高分辨率的多径TDOA估计算法主要可以分为以下几类：第1类是在辐射源信号已知的情况下将多径信号TDOA估计的多维优化问题转化为多个1维优化问题，从而降低优化问题求解的复杂度，这类算法主要有WRELAX (Weighted RELAXation)算法^[4,5]和期望最大化(Expectation Maximization, EM)算法^[6,7]。第2类是在辐射源信号未知的情况下，将时域的多径TDOA估计问题转化为对频域多个正弦信号的频率估计问题，然后通过子空间技术实现对TDOA的估计^[8,9]。第3类是利用时延值与信道冲激响应的稀疏支撑集之间的关系，将多径信号的TDOA估计问题转化为稀疏重构问题，通过合理地构造测量矩阵并结合不同的稀疏重构算法来估计TDOA^[10,11]。还有一些其他应用于高分辨率多径时间延迟估计的方法，如重要性采样方法^[12]和马尔可夫链蒙特卡罗方法^[13]等。

上述多径信号TDOA估计算法都是假设接收传感器处于高斯分布噪声环境中，但在实际的无线通信和水声通信环境中，自然现象或者人为因素的影响使得传感器的输入噪声具有明显的尖峰脉冲特性。这类尖峰脉冲特性噪声的概率密度分布曲线与高斯分布的概率密度曲线类似，但比高斯分布概率密度曲线具有更厚的拖尾。对称α-稳定分布(Symmetric α-Stable distribution, SαS) 是描述这类具有尖峰脉冲性的非高斯噪声的一种非常广泛的数学模型^[14]。SαS分布的概率密度可用3个参数来描述：特征指数α、分散系数γ和位置参数a。其中特征指数α($0 < \alpha \le 2$)表示分布函数的拖尾厚度，α取值越小，表示SαS分布的拖尾越厚，噪声脉冲性越强；α取值越大，表示SαS分布的拖尾越薄，噪声脉冲性越弱。当α=2时，SαS分布退化为高斯分布。当$0 < \alpha < 2$时，SαS分布不具有α阶及其以上阶矩。由于SαS分布的这个特性，使得经典的基于2阶统计量提出的各种TDOA估计算法的性能发生退化，甚至不能得到正确的估计结果。于是Shao等人^[15]提出了分数低阶统计量的概念并将传统的2阶和高阶统计量信号处理方法拓展到分数低阶统计量信号处理方法。

为了解决SαS分布脉冲噪声环境下的多径信号时延估计问题，文献[16]利用分数低阶统计量理论对高分辨率多径时延估计EM算法进行改进，提出一种基于分数低阶相关的P-EM算法，仿真结果表明P-EM算法对尖峰脉冲噪声具有更好的适应性。刘文红等人^[17]用分数低阶统计量理论中的最小分散系数准则代替最小二乘准则，并利用波形松弛思想将高维优化问题分解为多个1维优化问题，提出了一种韧性的高分辨率多径时延估计方法P-WRELAX。文献[18]将两路接收信号的分数低阶互协方差谱序列当作时间序列，利用子空间旋转不变技术实现了脉冲噪声环境下多径时延的高分辨率估计。基于分数低阶统计量理论的信号处理算法要求选取的分数低阶统计量的阶数p小于SαS分布的特征指数α，而在实际系统中SαS分布的特征指数是不易获得的，从而限制了基于分数低阶统计量理论的信号处理算法的应用。相关熵作为随机变量间局部相似性的度量不仅具有良好的抗噪性能，而且不依赖于脉冲噪声的先验知识，近年来受到了国内外学者的极大关注^[19,20]。文献[21]首先证明了在特定条件下最大相关熵准则和分数低阶统计量理论中的最小分散系数准则是等价的，然后提出了一种基于最大相关熵准则的自适应时延估计算法。文献[22]将最大相关熵准则分别与Farrow非整数时延滤波器和Lagrange非整数时延滤波器相结合解决了脉冲噪声环境下的非整数时延估计问题。文献[23]通过对两路接收信号的互相关熵谱密度进行白化加权处理提高了相关熵时延估计算法抗尖峰脉冲噪声的能力。文献[24]提出一种基于相关熵希尔伯特差值的时延估计算法，克服了窄带射频信号时延估计的精度受相对带宽限制和脉冲噪声影响的问题。文献[25]基于最大相关熵准则提出了一种针对时变线性系统状态估计的扩展递归最大相关熵自适应算法，并利用所提出的算法解决了ADS-B报文的缓变时延估计问题。文献[26]研究了脉冲噪声环境下振幅衰减信号的时延估计问题，基于最大相关熵准则和自适应时延估计思想提出了两步法时延估计算法；为了避免两步估计的误差累积和估计偏差，进一步提出了一种一步法自适应时延估计算法。受相关熵是随机变量相似性度量思想的启发，本文对期望最大化高分辨率多径TDOA估计算法进行改进，提出了基于相关熵的期望最大化算法(Correntropy based Expectation Maximization, CEM)。仿真实验表明，该算法在不依赖于脉冲噪声先验信息的情况下具有比P-EM算法和P-WRELAX算法更好的脉冲噪声适应性和估计精度。

2.   相关熵基本理论

两个随机变量X和Y之间的相似性可以用相关熵来度量，其定义为^[19]

其中，${\rm{E}}\left[\cdot\right]$表示期望值，${k_\sigma }( \cdot ) = \dfrac{1}{{\sqrt {2\pi \sigma } }}\exp \left[ - \dfrac{{{{( \cdot )}^2}}}{{2{\sigma ^2}}}\right]$表示高斯核函数，σ>0为核长参数。在工程应用中，两个随机变量X和Y的联合概率密度通常情况下是未知的，可以通过N组观测数据$\{ ({x_i},{y_i})\} _{i = 1}^N$近似估计相关熵，其近似估计表达式为

相关熵具有对称性和有界性，即${V}_{\sigma }\left(X,Y\right)={V}_{\sigma }\left(Y,X\right)$和${V}_{\sigma }\left(X,Y\right)\le {V}_{\sigma }\left(0\right)$。由相关熵可以诱导出一个距离测度CIM(Correntropy Induced Metric)。与常规的距离测度不同，CIM距离测度表现出由核长σ控制的混合范数特征。即当核长σ较大时，CIM距离测度的L₂范数距离区域也较大；当核长σ较小时，其L₂范数距离区域缩小而L₁范数距离区域和L₀范数距离区域扩大。基于CIM距离测度，可以定义最大相关熵准则为

3.   CEM多径TDOA估计算法

3.1   TDOA估计模型

在多径信号传播模型中，一般假设源信号$s(k)$是已知的或是可以通过某种方式获得的，传感器的接收信号可以表示为

式中，${\left\{{a}_{i}\right\}}_{i=1}^{D}$和$\{ {\tau _i}\} _{i = 1}^D$分别为各路径信号的衰减系数和时间延迟，D表示多径数目。$w(k)$为加性背景噪声，假设$w(k)$服从SαS分布并且与源信号互不相关。

3.2   基于相关熵的期望最大化算法

期望最大化算法是一种迭代优化策略，由于它的每一次迭代都分两步，其中一个为期望步(Expectation step, E-step)，另一个为极大步(Maximization step, M-step)，所以算法简称为EM算法。EM算法最初作为一种数据添加算法主要用于解决在数据量非常大的情况下存在数据缺失或者数据不可用的问题。随着理论的发展，EM算法已经不仅仅用于处理数据缺失问题，而成为一种用于解决复杂问题的简化策略。利用EM算法的基本思想，通过引入恰当的潜在数据可以将复杂的多维优化问题分解为若干个独立的1维优化问题。

EM算法用于解决多径信号TDOA估计问题时，首先引入潜在数据 ${\left\{{x}_{i}\left(k\right)\right\}}_{i=1}^{D}$

其中，$\{ {\hat a_i}\} _{i = 1}^D$和$\{ {\hat \tau _i}\} _{i = 1}^D$是待估计的参数且潜在数据中的噪声${w}_{i}\left(k\right)$是由$w(k)$分解而来的，即${w}_{i}\left(k\right)\!=\!{\beta }_{i}w\left(k\right)$。参数${\beta }_{i}$满足$\displaystyle\sum\nolimits_{i = 1}^D {{\beta _i}} = 1$，它的取值可以控制EM算法的收敛速度，并阻止算法收敛到局部极值点。通过EM算法的迭代优化最终要使潜在数据与传感器接收信号$r(k)$之间满足

因此EM算法收敛后的$\{ {\hat a_i}\} _{i = 1}^D$和$\{ {\hat \tau _i}\} _{i = 1}^D$即为各路径信号的衰减系数和时间延迟的估计值。

在给定参数${\hat a_i}$, ${\hat \tau _i}$和${\beta }_{i}$的起始值后，EM算法第m次迭代的过程为

E-step：对$i = 1,2, ··· ,D$计算

M-step：对$i = 1,2, ··· ,D$计算

式(8)是一个求解两个参数的最优化问题。它的求解过程可以首先根据期望步给定的$\hat a_i^{(m)}$求得$\hat \tau _i^{(m + 1)}$，然后再由$\hat \tau _i^{(m + 1)}$求解新的$\hat a_i^{(m + 1)}$。根据式(8)可得$\hat \tau _i^{(m + 1)}$的估计值为

然后由$\hat \tau _i^{(m + 1)}$可以得到$\hat a_i^{(m + 1)}$的估计值为

由于SαS分布脉冲噪声不存在有限的2阶矩，因此文献[16]利用最小分散系数准则代替了式(8)中的最小均方误差准则提出了一种适用于脉冲噪声环境下的鲁棒的EM算法，由于在该算法中迭代公式(9)和式(10)中的2阶相关被p-阶分数低阶相关所代替，因此称该算法为P-EM算法，其中参数1<p<α。

本文用最大相关熵准则代替EM算法中的最小均方误差准则，即在第m次迭代的M-step中

对$i = 1,2, ··· ,D$计算

可得$\hat \tau _i^{(m + 1)}$的估计值为

$\hat a_i^{(m + 1)}$的估计值为

观察式(7)可以发现$x_i^{(m)}(k)$是第i路多径信号的估计值，而式(12)相当于求${x}_{i}^{\left(m\right)}\left(k\right)$和第i路多径信号的相关熵，因此称该算法为基于相关熵的EM算法，记为CEM。

4.   算法仿真

取源信号$s\left(k\right)$为BPSK信号，载波频率为${f_{\rm{c}}} = 0.25{f_{\rm{s}}}$, ${f_{\rm{s}}}$为采样频率且${f_{\rm{s}}} = {10^5}$ Hz, ${T_{\rm{s}}} = {1 / {{f_{\rm{s}}}}}$为采样间隔，码元持续时间为$20{T_{\rm{s}}}$，信号带宽$B = {{0.1} / {{T_{\rm{s}}}}}$。分别考虑多径数目$D = 2$和$D = 3$两种情况，在$D = 2$的情况下，两路信号的TDOA和衰减系数分别是${\tau _1} = 28{T_{\rm{s}}}$, ${a_1} = 0.8$和${\tau _2} = 32{T_{\rm{s}}}$, ${a_2} = 0.6$。$D = 3$是在$D = 2$的情况下增加了${\tau _3} = 36{T_{\rm{s}}}$, ${a}_{3}=0.7$的第3路多径信号。各路多径信号之间的TDOA之差为$\Delta \tau = \left| {{\tau _3} - {\tau _2}} \right| = \left| {{\tau _2} - {\tau _1}} \right| = 4{T_{\rm{s}}} < {1 / B} = 10{T_{\rm{s}}}$，属于基本相关法中的多径不可分辨情况。

仿真中3路信号的TDOA和衰减系数初值分别设为0及0.5, ${\beta }_{1}$, ${\beta _2}$和${\beta _3}$分别设为0.25, 0.40和0.35，实验数据长度为3000点。由于SαS分布不具有有限的2阶矩，定义广义信噪比${\rm{GSNR = 10lg}}({{\sigma _{\rm{s}}^2} / {{\gamma _{\rm{n}}}}})$，其中${\gamma _{\rm{n}}}$为SαS分布的分散系数，${\sigma }_{{\rm{s}}}^{2}$为信号的方差。采用均方根误差(Root Mean Square Error, RMSE)和平均准确率${{{P}}_{\rm{A}}}$来衡量EM算法、P-EM算法、P-WRELAX算法和CEM算法的估计性能。它们的定义分别为

其中，${\hat \tau _i}(l)$是${\tau }_{i}$的第l次仿真实验的估计值，$l'(i)$是在L次总仿真实验次数中第i径TDOA得到与设定的时延值完全一致估计结果的次数。实验中总仿真实验次数L=200, P-EM算法和P-WRELAX算法中的参数p=α–0.1，CEM算法中的核长σ=10。

实验1　给定SαS分布的特征指数α=1.4，图1给出了EM算法、P-EM算法、P-WRELAX算法和CEM算法分别在2径和3径情况下的均方根误差和平均准确率随着广义信噪比GSNR变化的曲线。从图1可以明显看出，4种算法在2径情况下的估计性能都要优于3径情况下的估计性能，并且无论在2径还是3径情况下采用了脉冲噪声抑制技术的P-EM算法、P-WRELAX算法和CEM算法的性能都要优于没有采用脉冲噪声抑制技术的EM算法。特别是在2径情况下P-EM算法、P-WRELAX算法和CEM算法在GSNR>–2 dB时平均准确率为100%，均方根误差为零，而EM算法要到GSNR>–8 dB时才会得到平均准确率为100%，均方根误差为零。与基于分数低阶矩的P-EM算法和P-WRELAX算法相比，基于相关熵的CEM算法在低信噪比下的性能更有优势，例如在2径情况下当GSNR=–8 dB时P-EM算法的准确率为58.5%, P-WRELAX算法的准确率为64.5%，而CEM算法的准确已经达到90%以上。

图 1  不同广义信噪比下各算法的估计性能

下载: 全尺寸图片 幻灯片

实验2　给定广义信噪比GSNR=–3 dB，图2给出了EM算法、P-EM算法、P-WRELAX算法和CEM算法的均方根误差和平均准确率随着SαS分布噪声的特征指数α变化的曲线。从图中可以看出在2径和3径的情况下P-EM 算法、P-WRELAX算法和CEM算法的平均准确率和均方根误差在各种脉冲噪声环境下的性能都要优于EM算法，说明P-EM算法、P-WRELAX算法和CEM算法均具有对脉冲噪声的抑制作用。但是基于相关熵的CEM算法对脉冲噪声的适应性要优于基于分数低阶矩的P-EM算法和P-WRELAX算法，特别是在脉冲噪声的特征指数α<1的强脉冲噪声环境下，例如在3径情况下当α=0.6时，P-EM算法和P-WRELAX算法的平均准确率只有32.2%和60.9%，而CEM算法的平均准确率为99.1%，这是因为这两种算法中利用的p阶分数低阶矩虽然在α<1的情况下可以通过求时间平均的方法来计算，但理论上是没有定义的；而CEM算法中利用的相关熵在特征指数$0 < \alpha \le 2$的情况下都是有定义的。

图 2  不同特征指数α下各算法的估计性能

下载: 全尺寸图片 幻灯片

实验3　给定广义信噪比GSNR=–8 dB，图3给出了CEM算法在SαS稳定分布噪声的特征指数α=0.5, 1.0, 1.5和2.0时的估计性能随着核长参数σ的变化曲线。从图中可以看出，在某一给定特征指数的SαS稳定分布噪声环境下，为了获得最佳的估计性能存在着最佳核长选取问题，且最佳核长值与SαS稳定分布噪声的特征指数密切相关。当给定的核长较小时(如σ=1)，随着噪声脉冲性的增加，CEM算法的估计性能越来越好；当给定的核长较大时(如σ=14)，随着噪声脉冲性的增加，CEM算法的估计性能却越来越差。这是因为当核长较小时，相关熵诱导的CIM距离测度的L₁和L₀范数距离空间增大，L₁和L₀范数对脉冲噪声具有很好的抑制作用，所以随着噪声脉冲性的增加，CEM算法的估计性能也越来越好。当核长较大时，相关熵诱导的CIM距离测度的L₁和L₀范数距离空间被压缩，与信号2阶相关等价的L₂范数距离空间增大，也就是说此时相关熵对脉冲噪声的抑制能力减弱却更多地聚集了信号中的2阶相关信息，所以随着噪声脉冲性的增加，CEM算法的估计性能下降。CEM算法中核长参数的选择虽然与脉冲噪声的特性有关，但它与P-EM算法中分数低阶参数p(p<α)必须依赖于脉冲噪声的先验信息是不同的。

图 3  不同核长σ和不同特征指数α下CEM算法的估计性能

下载: 全尺寸图片 幻灯片

实验4　给定SαS分布的特征指数α=0.5，图4给出了CEM算法在广义信噪比GSNR=–12 dB, –8 dB, –4 dB和0 dB时的估计性能随着核长参数σ的变化曲线。从图中可以看出，在4种广义信噪比的环境下，CEM算法能够获得最佳估计性能的核长参数的取值范围都在1.4～1.8范围内，结合实验3可以说明相关熵函数中核长参数的最优选取与SαS分布的脉冲特性密切相关，而与表征脉冲噪声强弱的广义信噪比无关。在相同脉冲特性不同脉冲强度的SαS分布噪声环境下，当核长参数的取值在大数范围内时，为了获得相同的估计性能，对于强度较大的脉冲噪声应该选取较小的核长。例如图4中在σ>4的范围内，为了获得99.5%的平均准确率，当GSNR=0 dB时，应选取的核长为σ=10，而当脉冲噪声的强度增加到GSNR=–4 dB时，核长参数降为σ=6。

图 4  不同核长σ和不同广义信噪比GSNR下CEM算法的估计性能

下载: 全尺寸图片 幻灯片

5.   结束语

为了解决基于EM算法的多径信号TDOA估计算法在脉冲噪声环境下性能退化的问题，本文基于最大相关熵准则提出了基于相关熵的期望最大化高分辨率TDOA估计算法(CEM)。仿真实验结果表明CEM算法和基于分数低阶矩的P-EM算法和P-WRELAX算法一样对脉冲噪声具有很好的抑制作用，并且在强脉冲噪声和低信噪比的环境下CEM算法的估计性能要优于基于分数低阶矩的算法，而且CEM算法中参数的选取不依赖于脉冲噪声的先验信息。