陈 健, 田丽芳

(1. 山东服装职业学院, 山东 泰安 271000； 2. 黄淮学院 信息工程学院, 河南 驻马店 463000)



基于重要性采样的最大似然时延估计算法

陈 健1, 田丽芳2

(1. 山东服装职业学院, 山东 泰安 271000； 2. 黄淮学院 信息工程学院, 河南 驻马店 463000)

针对多径环境中的时延估计，提出了一种基于重要性采样概念的算法。该算法利用蒙特卡罗算法(MC)对未知参数的分布函数抽样，获得简化似然函数的全局最优解，进而通过计算样本均值直接得到参数估计结果。该方法避免了耗时较长的多维网格搜索和对初值较为敏感的迭代算法，且能够无条件收敛至全局最优值。仿真结果表明，在相同样本条件下，该算法相比于EM、MUSIC算法，不仅消除了对初值的依赖性，也获得了更接近克拉美罗界(CRLB)的仿真结果。将该算法与其他多种算法进行计算复杂度分析后发现，IS-based算法较其他算法更为简单，计算量更低，具有较为重要的工程应用价值。

时延估计； 最大似然； 多径环境； 重要性采样

0 引 言

近年来，对多径时延估计的经典处理方法一般由Bell等提出的匹配滤波技术[1]解决，该方法通过求接收信号对于发射信号的相关峰值来得到时延估计，能在较低信噪比下获得较高估计精度，但其时延分辨能力受限于发射信号自相关函数的主瓣宽度，对发射和接收机要求过高。MUSIC算法[2-4]将多径时延估计问题在频域等效成正弦频率估计问题，获得较高分辨率的估计算法，但其性能仅限于在频谱平坦或近似平坦的情况下较优。最大似然(ML)[5]准则用于多径时延的估计，并且可以实现时延的高分辨率，对信号形式也没有附加要求。其缺点是时延估计的处理过程中涉及到复杂的多维优化问题，因此在多径时延问题中直接使用ML估计方法并不可行。出于降低优化维数的目的，Feder等提出了基于ML准则的EM算法[6]，将多维优化的复杂问题分解成一系列一维优化问题的迭代，但是它对初值非常敏感，且没有系统的初始化方法。随着技术发展，新的信号处理算法应用到时延估计中。蒙特卡罗算法(MC)[8-9]作为一类非常重要的数值计算方法应用于时延估计中。

本文从参数最大似然估计出发，构建信号频域噪声概率密度函数，推导时延估计的似然函数并构造了一个新的代价函数，从而保证该算法总可以收敛到全局点。同时针对新的代价函数采用基于MC的重要性采样方法来求解该函数的最优解，并通过仿真试验将计算结果与已有的EM算法、MUSIC算法作对比。结果表明，该算法不仅能更好地估计时间延迟，同时其算法复杂度较其他算法也更低，证明了该算法的优越性。

1 信号模型

多径传播的一般性信号模型为

(1)

式中：x(t)为t时刻发射信号；参数αi和τi分别表示第i条传播路径的增益系数和时间延迟；y(t)为接收信号；ω(t)为t时刻的观测噪声；P为总的多径个数。

对接收信号离散采样，并对采样信号进行离散傅里叶变换(DFT)，得到接收后的信号模型：

(2)

将式(2)转化为向量的形式：

(3)

(4)

式中：

进而，构建参数τ和α的似然函数，表示如下：

(5)

式中：p(Y;τ,α)为Y的概率密度函数；σ2N为噪声功率的方差。因此，为求参数的最优估计值，可以通过求式(5)的最大值，进而代入式(2)，求得时延估计。

然而，p(Y;τ,α)是关于参数τ和α联合的非线性函数，且是α的二次函数，较难直接求解。同时，α不是此问题所关注的参数，因此α关于参数τ的估计值的解析式可以通过求偏导的方法得到，即假设在τ已知的情况下，α的估计值为：

(6)

将式(5)式代入(6)得：

(7)

(8)

2 算法简介

2.1 简化似然函数及全局寻优

根据似然函数的物理含义可知，式(8)以τ为变量生成的伪谱最大值对应的τ即为时延的ML估计[10-12]。但由于简化似然函数Lc(τ)是关于自变量τ的非线性函数，其理论解是非常难求的。在关于此问题的研究中，最常用的即是迭代算法，用以替代更为复杂耗时的网格搜索算法。但是,迭代算法的性能取决于是否有一个好的初始值，而且极有可能收敛于局部最优点，无法找到全局最优值。在这里我们使用Pincus的全局最大值定理[13]，它不需要初始化并且保证总能收敛到全局最大值点。其基本原理：如果Lc(τ)存在全局最大值且最大值为正，则最大值点可近似为：

(9)

如果定义ρLc(τ)的归一化函数为

(10)

(11)

2.2 重要性采样技术

IS[14-16]作为一种高效的针对稀有事件的降方差算法，是计算多重积分问题的有效手段。其基本思想是通过一个相对简单的分布函数生成一定数量的随机数，并将其加权平均来近似计算目标分布函数的数学期望。其中，相对简单的分布函数被称为重要性函数(Importance Function，IF)[17]，权重值近似正比于这两种分布的似然比。通过修改重要性函数，并引入重要性权值可以大幅度地减少仿真样本数，从而在较短的运行时间内得到给定精确度的仿真结果。

(12)

称为权重函数。

对g′(τ)采样，得到τ的R个独立分布的样本点，故式(11)的积分可以用下面的样本平均来逼近，即

(13)

权重函数w(x)用来调整采样点输出，以确保对τ的重要性采样估计为无偏估计。

由此可以看出，重要性采样的关键在于选取合适的偏置函数，以增加对于待估值参数影响更大的变量的采样。好的偏置函数可以极大地缩短仿真时间；相反，偏置函数一旦选取不好，就会大大延长仿真时间，甚至比传统的蒙特卡罗方法还要长。

重要性函数的选取对估值算法的性能起着举足轻重的作用。因此，我们面临的关键问题则是选取合适的重要性函数，在减轻计算量的同时，也要能够保证良好的估值性能。

(14)

经统计验证表明,ΦH(τ)Φa(τ)的非对角线上的值要远小于其对角线上的值，因此，可以将ΦH(τ)Φa(τ)近似为对角矩阵，即

(15)

式中：Ip为P×P的单位矩阵。此时，定义重要性函数为：

(16)

经一系列代数运算化简后，得：

(17)

式中，

并且ρ1是一个不同于ρ0的常数。最终，重要性采样函数gρ1(τ)写为如下形式：

(18)

式中，

(19)

一旦确定了重要性采样函数，多径时延估计就可以由下式给出：

(20)

式中：τk是τ的第k个向量，且τk(i)是τk的第i个元素。

该算法整体流程如下：

form=1,R

forn=1,P

令ξ中不再含有τ1,τ2,…,τn；

end for

end for

由之产生的数组即为τ矢量的变现，进而可以由之得到时间延迟的估计值。

3 结果分析与讨论

3.1 采样参数

由推导过程可以看出，ρ0、ρ1和样本数R的对估值的准确度和计算量有着重要影响。其中，参数R的大小决于参数τ的具体取值，即R的取值大小完全取决于重要性采样函数和简化似然函数之间的相似度。同时由式(10)可以看出，当参数ρ0趋向于无穷大时，简化似然函数呈现P维的Dirac-delta函数形状，使全局最大值和其他极值点有很大的区别。实际上，如果当函数的全局最大点和其他极值点已经分开，则ρ0的取值不必太大 ；否则只有ρ0较大时，才能满足计算要求。参数ρ1能够保证重要性采样函数与最大函数形状的类似程度，并且描述的伪概率密度函数具有“重尾”特性，因而应取小值。

图1、2给出了参数ρ0和ρ1的选取对算法的估计性能的影响。从图1可以看出，估值的标准差(RMSE)随着ρ0的增大，先缓慢减小，后迅速下降，并随着ρ0的继续增大几乎保持不变。这表明，当全局最大点和其他极值点已经分开，此时的ρ0的取值已经能够满足计算要求，继续增大ρ0对结果改善不大。与之相反，从图2可以看出，随着ρ1的增大，估值标准差先缓慢变化，后迅速增加，并在之后随着ρ1的继续增大几乎保持不变。这说明，当ρ1的值过大时，会导致算法性能急剧变差，可能会导致某些延时峰消失。在本文的仿真中，令ρ0和ρ1分别为16和4。

图1 参数ρ0对算法性能的影响曲线

图2 参数ρ1对算法性能的影响曲线

3.2 仿真结果与分析

本文采用高斯分布的白噪声叠加8PSK信号对该算法的估值性能进行仿真分析。为了更好地对新算法性能做出评价，本文针对上述问题，将该算法与EM、MUSIC算法的计算结果作比较，同时也将上述结果与文献[7]的CRLB理论界限作对比，用作评价不同算法性能的指标。其中参数设置如下：样本数N=100，采样间隔Ts=1×10-4s,多径路数P=3，假定各路径增益相同且3个路径时延分别为9Ts，12Ts，20Ts。

图3所示为本文所提算法的重要性函数(IF)随τ/Ts的变化。从图中可以看出，重要性函数在不同位置出现3个延时峰，其对应位置分别为τ/Ts=9,12,20，这与前面的时延假设完全吻合，表明本文所取重要性函数是可信、可靠的。

图3 SNR=10 dB时重要性函数分布

在图4中，固定SNR为0 dB ,样本数R从30～150，采用1 000次Monte Carlo 仿真来计算估值标准差。由图4可以看出，随着样本数的增加，CRLB理论界限不断下降，但是下降速度逐渐变慢并存在极限，这也表明所有算法性能也存在极限值。

如图4(a)所示，对于对初值有强烈依赖性的EM算法，可以看出该算法对于较好的初值，能够得到较低，甚至接近于CRLB理论界限的标准差。然而，当初值不够好时，则无论样本数如何增大，都无法接近CRLB界，这说明此时该算法收敛于局部最优解。而本文提出的新算法，则总能够随着样本数的增加迅速接近CRLB界，表明了该算法的优越性。

(a)

(b)

同时，从图4(b)可以看出，尽管MUSIC算法消除了对初值的依赖性，其估值性能依然逊色于本文提出的IS-based算法。

(a)

(b)

从图5可以看出，随着SNR的增大，不同算法的RMSE随之减小，同时CRLB界在指数空间随SNR线性减小。

如图5(a)所示，EM算法对初值依赖性较大。对于较好的初值，不仅能够快速收敛，同时能够保证收敛于CRLB界；对于较差的初值，则收敛较慢且不能收敛到CRLB界。同时，本文提出的IS-based算法与EM算法对比可以看出，该算法能够无条件收敛到全局最优值，消除了对初值的依赖性，并且该算法收敛速度和精度与有良好初值的EM算法类似。同理，由图5(b)可知，IS-based算法性能要优越于MUSIC算法。

3.3 算法复杂度分析

衡量一个算法的优劣，除了以其求解精准度作为判别标准外，还必须将该算法的计算复杂程度考虑在内。表1对比了本文提出的算法与网格搜索、EM、MUSIC算法的复杂度。其中，Niter代表对应算法的迭代次数；M为算法引入参数；K为相关常数。可以看出，本文所提算法更为简单，计算量更小。

表1 算法复杂度对比

4 结 语

本文针对多径环境下的时延估计问题，提出了基于重要性采样的最大似然估计算法。本文对比了EM算法、MUSIC算法与本文所提IS-based算法，结果表明，该算法不仅消除了初值依赖性，也能够快速收敛至全局最优解，并且使用样本数较其他算法更少。本文将该算法与其他多种算法进行计算复杂度分析后发现，IS-based算法较其他算法更为简单，计算量更低。

综上可知，本文提出的IS-based算法计算精准度高，且算法简单，计算量较低，具有良好的实际应用价值。

[1] Bell B M, Ewart T E. Separating multipaths by global optimization of multidimensional mathched filter [J]. IEEE Trans Acoust, Speech, and Signal processing, 1986, 34(5) :1029-1037.

[2] Schmidt R O. Multiple emitter location and signal Parameter estimation [J]. IEEE Trans. Antenna and Propagation, 1986, 34(3):276-280.

[3] 陈 栋，查代奉，杨保海. 基于FLOC-MUSIC算法的3D-DOA估计[J].计算机应用研究，2015(10):32.

[4] 钟林波，吉建华，文小军,等.基于角度预估计的声源定位MUSIC算法实验研究[J].计算机工程与应用，2015,51(2)：205-208.

[5] 陈华伟,赵俊渭,郭业才. 一种频域自适应最大似然时延估计算法[J]. 系统工程与电子技术,2003,11:1355-1357,1361.

[6] Feder M, Weinstein E. Parameter estimation of superimposed signals using the EM algorithm[J]. IEEE Trans Acoust, Speech Signal Process, 1988,36:477-489.

[7] 黄 龙，郭树人，徐 昂，等.基于脉冲串信号的无源定位系统频差估计误差CRBL分析[J].国防科技大学学报，2013，35(2)：109-114.

[8] 王惠刚,刘 强. 多普勒方位联合估计的蒙特卡罗算法[J]. 电子学报,2009,33(2):427-431.

[9] 李 晶, 赵拥军,李冬梅. 基于马尔科夫链蒙特卡罗的时延估计算法[J]. 物理学报,2014,63(13):1-7.

[10] 胡 德，郭刚正.最小二乘法、矩法和最大似然法的应用比较[J]. 统计与决策，2015(9)：20-24.

[11] 朱娟娟，王 伟.OFDM中基于CP的最大似然估计同步技术[J]. 信息技术，2015(1)：125-128.

[12] 周 成，黄高明，单鸿昌，等.基于最大似然估计的TDOA/FDOA无源定位偏差补偿算法[J]. 航空学报，2015,36(3)：979-986.

[13] PINCUS M M. A closed form solution for certain programming problems [J]. Oper Research,1968(18):1225-1228.

[14] Masmoudi A, Bellili F, Affes S. A non-data-aided maximum likelihood time delay estimator using imortance sampling[J]. IEEE Trans Signal Process, 2011,59(10):4505-4515.

[15] 张玲霞，刘志仓，王 辉，等.非线性系统故障诊断的粒子滤波方法[J].电子学报，2015,43(3)：615-619.

[16] 刘 璞，王春平，徐 艳.基于Mean Shift重要性采样的粒子滤波跟踪算法[J].计算机工程与设计，2014,35(3)：909-913.

[17] 邱有恒，邓 力，李百文，等.集中重要性函数在粒子输运蒙特卡罗模拟中的应用[J].原子能科学技术，2013,47：673-677.

Research on an Importance-Sampling Based Maximum Likelihood Time Delay Estimator

CHENJian1,TIANLi-fang2

(1. Shandong Vocational Institute of Clothing Technology, Taian 271000, China;2. School of Information Engineering, Huanghuai University, Zhumadian 463000, China)

To estimate the time delay in multipath environments, a new algorithm was proposed based on importance sampling (IS). This algorithm applied Monte Carlo method into the sampling of unknown parameters. The global maximum of the compressed likelihood function could be found, thus, the parameters were approximated by the average of samples. Therefore, the estimator proposed avoided the complex multidimensional grid search and iterative methods which depended seriously on the initial guess. The algorithm guaranteed convergence to the global maximum. Results of simulation showed that, this algorithm held better performance than EM and MUSIC algorithms, under the conditions of the same samples. In addition, the algorithm also showed superiority in terms of estimation complexity. In addition, IS-based algorithm proposed demonstrated superiority of other ones in terms of algorithm complexity.

time-delay estimation; maximum likelihood (ML); multipath environments; importance sampling (IS)

2015-03-12

河南省科技厅发展计划(142102110088)；河南省科技攻关项目(No.122102210430；No.142102110088)

陈 健(1975-)，男，山东泰安人，硕士，副教授，主要从事计算机人工智能算法及其相关算法的理论和应用研究。

Tel.: 0538-6959656; E-mail: 274912359@qq.com

TN 911.7

A

1006-7167(2015)12-0016-05