利用蒙特卡罗的最大似然时延估计算法

巴斌,郑娜娥,朱世磊,胡捍英

(解放军信息工程大学导航与空天目标工程学院,450001,郑州)



利用蒙特卡罗的最大似然时延估计算法

巴斌,郑娜娥,朱世磊,胡捍英

(解放军信息工程大学导航与空天目标工程学院,450001,郑州)

针对最大似然时延估计算法的峰值搜索计算复杂度较高且容易陷入局部收敛,造成估计误差较大的问题。提出了一种利用蒙特卡罗的最大似然时延估计(MCML)算法。首先利用信道频域响应估计矢量建立似然函数;然后把时延估计问题转化为求解随机变量的期望问题,将采用指数化似然函数构造的标准化概率密度函数趋近于冲激函数,使得随机变量的方差趋近于零;最后采用蒙特卡罗方法对随机变量进行抽样,从而利用抽样的均值估计出时延。较之传统方法,蒙特卡罗方法避免了网格搜索,降低了计算复杂度,保证了全局收敛性和估计精度。仿真结果表明:在信噪比0～25 dB的条件下,MCML算法均能始终逼近克拉美罗界;当信噪比为25 dB时,MCML算法的时延估计分布范围缩小为马尔科夫链蒙特卡罗算法的34%。

最大似然;蒙特卡罗;时延估计;克拉美罗界

在定位系统中,测量端通过接收辐射源信号精确估计到达时间(time of arrival, TOA),从而实现高精度目标定位[1]。因此,TOA估计一直都是国内外十分关注的研究热点,其中超分辨TOA估计技术是研究的重点。该项技术被广泛应用于声呐、雷达[2]和通信中的目标定位[3]等各个领域。

超分辨TOA估计算法主要包括两类:子空间类估计算法和最大似然估计算法。子空间类估计算法包括最小范数谱估计算法[4]、多重信号分类(multiple signal classification, MUSIC)算法[5-6]、传播算子算法(propagator method, PM)[7]、求根MUSIC算法[8]、旋转不变技术估计信号参数(estimation of signal parameters via rotational invariance techniques, ESPRIT)算法[9]等。这些算法在较为理想的应用环境中具有较好的估计性能,但是在恶劣环境下,特别是低信噪比条件下,算法的性能会有所降低。最大似然估计算法被公认为最佳估计算法,不仅在高信噪比时性能逼近克拉美罗界(Cramer-Rao bound, CRB),并且在低信噪比时也具有很好的性能。然而,最大似然估计算法需要进行网格搜索,计算量巨大,工程上难以实现。

为了解决最大似然时延估计算法计算量巨大的问题,必须采用新的信号处理方法。蒙特卡罗(Monte Carlo, MC)方法是一类非常重要的数值计算方法,但与一般的数值计算方法有本质区别,它依据概率论中的大数定律,采用统计抽样理论近似求解数学问题。文献[10]利用MC方法对热辐射输运问题进行隐式求解。文献[11]利用可逆跳转马尔科夫链蒙特卡罗(reversible jump Markov Chain Monte Carlo, RJMCMC)算法可以同时解决检测和估计问题,该算法需要在所有的参数子空间上跳转,因此效率低、计算量巨大。文献[12]和文献[13]利用蒙特卡罗重要性采样的方法实现多径条件下的最大似然时延估计,进而将蒙特卡罗思想引入到了到达时间差(time difference of arrival, TDOA)估计中。文献[14]将马尔科夫链蒙特卡罗(Markov chain Monte Carlo, MCMC)应用于无源雷达系统中的最大似然TDOA估计,对于单快拍条件下的TOA估计具有重要的借鉴意义。然而,文献[10-14]中的蒙特卡罗方法虽然解决了相应的问题,但是约束条件较多。

本文利用信道频域响应估计建立信号模型,受约束的条件更少,并且充分利用了单径条件下似然函数可以简化计算的特性。相比于文献[14]的MCMC方法,本文蒙特卡罗直接抽样方法所需的抽样点数更少,时延估计精度更高。

为减少约束条件,并简化计算,本文推导了TOA估计的似然函数,构建了标准化概率密度函数,利用MC方法对标准化概率密度函数进行直接抽样,并通过计算样本均值得到TOA估计值,最后,通过仿真对比MCML算法与MUSIC算法及MCMC算法的估计性能并分析了计算复杂度。

1 信号模型

为了分析接收信号的频谱,实现时延估计,接收端必须准确已知发送端发送信号的时间。在空旷的室外环境下,收发两端可以通过GPS实现时间同步。在城市峡谷或者室内环境,由于GPS信号微弱,因此收发两端无法通过GPS实现时间同步。时延估计的目的主要是用来完成无线电测距,而无线电测距可以通过测量飞行时间的方式来实现。例如,发送端在发送的信号中标记本地时间戳t1,接收端收到信号并在反馈信号中标记本地接收时间戳t2和发送反馈信号时间戳t3,发送端收到反馈信号后记录时间t4,那么[(t4-t1)-(t3-t2)]/2即为信号由发送端到接收端的空中飞行时间,飞行时间乘以光速即为两端的距离。t1和t4、t2和t3分别共时间基准,飞行时间的计算利用时间差避免了收发两端的时间同步。

在定位系统中,无线信道脉冲响应可以建模为

(1)

式中:a=|a|ejφ为无线信道的复衰落系数,|a|为幅度,φ为相位并在区间(0,2π)上服从均匀分布[6],记为φ～U(0,2π);τ为传播时延,即TOA;δ表示冲激函数。对式(1)进行傅里叶变换,信道响应的频域表示形式为

(2)

式中:f为频率域。

信号模型利用频域信道响应进行TOA估计。信道频域响应的离散采样在不同的系统中可以采用不同的方法获得,如在正交频分复用(orthogonal frequency division multiplexing, OFDM)系统中采用多载波解调技术,在直接序列扩频系统中采用接收信号解卷积方法等。

对频域信道响应H(f)进行K个等频率间隔采样,考虑到测量过程中的加性高斯白噪声,频域信道响应的离散采样可以表示为

x(k)=H(k)+w(k)=ae-j2π(fc+kΔf)τ+w(k)

(3)

式中:k=0,1,…,K-1;fc为载波频率;Δf为频域采样间隔;w(k)表示均值为0、方差为σ2的加性复高斯白噪声,记为w(k)～N(0,σ2)。信号模型的矢量形式可以表示为

(4)

式中:x=[x(0),x(1),…,x(K-1)]T和H=[H(0),H(1),…,H(K-1)]T分别为频域信道估计矢量和频域信道响应矢量。为简化计算,并在下文的求导过程中消除无关参数,定义α=ae-j2πfcτ、V(τ)=[1,e-j2πΔfτ,…,e-j2πΔf(K-1)τ]T为导向矢量,w=[w(0),w(1),…,w(K-1)]T为加性复高斯白噪声矢量。

由式(4)和噪声的相关假设,似然函数可以表示为

(5)

(6)

为了简化时延估计过程,首先消除无关参数α,得到近关于时延的似然函数。由于p(x|τ,α)是τ和α的联合分布函数,并且α为二次形式,可以通过求偏导的方法得到α关于τ的估计值解析式

(7)

将式(7)代入式(5),取对数,并去掉常数部分,可得τ的似然函数

(8)

由于

(9)

式(8)可以化简为

(10)

式(10)是关于时延τ的表达式,从而简化了似然函数。通过以上分析,信号的处理流程如图1所示。

图1 信号处理流程图

2 蒙特卡罗时延估计算法

为了计算似然函数L(τ)的全局最大值,首先对L(τ)取指数得到指数似然函数,使得指数似然函数更加尖锐;其次利用指数似然函数构造标准化概率密度函数,标准化概率密度函数趋近于冲激函数,并且其峰值所对应的时延(等价于标准化概率密度函数的期望)即为时延估计;然后进行随机抽样并利用统计平均代替积分得到时延估计,由于标准化概率密度函数趋近于冲激函数,因此抽样方差趋近于零,从而保证全局收敛和高精度的时延估计。

2.1 似然函数的全局最大值

为了使样本均值逼近全局最大值,根据文献[13],对似然函数L(τ)取指数能够使分布函数更加尖锐,从而使估计更加精确。定义指数似然函数

(11)

式中:ρ为常数。

(12)

式中:J表示时延搜索区间。

定义标准化概率密度函数

(13)

(14)

由于fρ0(τ)趋近于冲激函数,因此对fρ0(τ)的抽样方差趋近于0,从而不仅保证了全局收敛,还保证了时延估计的高精度。

(15)

2.2 随机抽样方法

(16)

式中:n=1,2,…,N。

对式(16)进行随机抽样,首先利用均匀分布U(0,1)生成矢量u=[u1,u2,…,uM],计算τi=F-1(ui)。其中,F-1(·)是概率累积分布函数的逆函数,F(x)定义如下

(17)

根据式(16),概率累积分布函数可以表示为

(18)

式中:n=1,2,…,N。

逆函数F-1(x)的闭式表达式不易获得,因此第m次采样的抽样公式为

(19)

2.3 算法步骤

由以上推导与分析,MCML算法的流程可以归纳如下:

(2)利用均匀分布U(0,1)生成矢量u=[u1,u2,…,uM];

(3)根据式(19)得到M次抽样τ(m),其中m=1,2,…,M;

（1）加快企业现金流出。根据权责发生制，企业赊销所获得的收入全部记入当期收入，企业通过赊销获得了相应的经济收益，企业当期纳税义务随之产生，企业当期并没有获得现金的流入，却要帮客户垫付各项税费，这就导致了现金流出企业，应收账款的产生增加了企业的费用，企业管理这些应收账款所需要支付的成本都加快了企业现金的流出。

3 克拉美罗界

克拉美罗界给出了无偏估计子的均方误差下界,下面给出模型相应的克拉美罗界。为了便于计算,同时不失一般性,对公式(5)取对数,从而将指数形式转化为求和形式

lnplnp(x|τ,α)=-Kln(2π)-Kln(σ2/2)-

(20)

定义αRe和αIm分别表示α的实部和虚部,即αRe=Reα,αIm=Imα。克拉美罗界可以利用Fisher信息矩阵(Fisher information matrix, FIM)求得。FIM可以定义为E(γγT),γ=∂lnp/∂[σ2,αRe,αIm,τ]。因此,在求克拉美罗界时,需要首先计算出lnp对σ2、αRe、αIm和τ的一阶偏导数,可得

(21)

(22)

(23)

(24)

式中

(25)

根据文献[16],利用式(21)～式(24),可得E(γγT)中的元素

(26)

(27)

(28)

(29)

(30)

(31)

(32)

根据FIM矩阵及文献[16],τ的克拉美罗界Bcrb(τ)满足下式

(33)

4 仿真实验及性能分析

4.1 仿真实验

为了评估MCML算法的性能,将在OFDM系统下,本节对MCML算法与文献[6]中的MUSIC时延估计算法、文献[14]中的MCMC时延估计算法以进行对比分析。最后对以上算法的计算复杂度进行比较分析。仿真采用高斯信道,信道复衰落系数的幅度|a|设置为1,传播时延τ设置为200 ns。OFDM系统仿真参数设置如表1所示。

表1 OFDM系统参数设置

首先定义信噪比(signal to noise ratio, SNR)

(34)

定义均方误差(mean square error, MSE):

(35)

仿真1 似然函数与指数似然函数性能对比。考虑到ρ0取值较大时计算溢出,分别取ρ0=1和ρ0=5时的指数似然函数与似然函数的归一化结果进行对比,2种函数的性能对比结果如图2～图4所示。

从图2～图4可以看出,指数似然函数全局最大值明显突出,曲线非常尖锐。对比ρ0=1和ρ0=5时的归一化指数似然函数可以发现,随着ρ0的增大,归一化指数似然函数变得更加尖锐,逼近冲激函数。归一化指数似然函数越尖锐,表示时延随机变量的方差越小,从而利用随机抽样平均进行时延估计的误差越小。

图2 归一化后的似然函数

图3 ρ0=1时归一化后的指数似然函数

图4 ρ0=5时归一化后的指数似然函数

仿真2 MCML算法与文献[6]和文献[14]算法性能对比。参数设置为:样本数M=2 000,做1 000次独立统计实验。考虑到ρ0取值较大时计算溢出,取ρ0=5。这是由于当ρ0=5时,归一化指数似然函数已经足够尖锐。

由图5可以看出,随着信噪比的增大,MUSIC算法和MCML算法估计的均方误差迅速减小,而MCMC算法在高信噪比时的均方误差变化趋于平缓。

图5 时延估计性能比较

对于MUSIC算法,由于单快拍下协方差矩阵的秩为1,无法分解出信号子空间和噪声子空间,因此需要采用频域平滑使协方差矩阵满秩。频域平滑导致有效带宽减小,从而导致时延估计方差始终与CRB有一定距离。

对于MCMC算法,由于马尔科夫链具有一定的拒绝概率,因此,在高信噪比下为了提高时延的估计精度,需要增大马尔科夫链的样本点。MCMC算法的样本点条件与MCML算法相同。

MCML算法由于利用直接抽样的方法,始终能够保证抽样均匀分布于整个归一化后的指数似然函数时延区间。并且ρ0=5时,归一化后的指数似然函数逼近冲激函数,因此抽样平均的方差更小,时延估计精度更高。因此当信噪比在0dB到25dB的条件下,MCML算法均能始终逼近克拉美罗界。

仿真3MCML算法与文献[14]算法时延估计分布。仿真参数设置与仿真2一致,对比MCML算法与MCMC算法时延估计的散布图。

图6给出了MCML算法与MCMC算法的时延估计散布图。从图中可以看出,随着信噪比的提高,MCML算法与MCMC算法的时延估计值分布趋于紧密,估计误差逐渐缩小;在信噪比较小时,MCML算法时延估计的分布范围与MCMC算法相当;在信噪比较大时,MCML算法的时延分布范围明显小于MCMC算法。当信噪比为25dB时,MCML算法的时延估计分布范围缩小为MCMC算法的34%。因此MCML算法的时延估计精度优于MCMC算法。

图6 MCML算法与MCMC算法的时延估计散布图

4.2 复杂度分析

通过仿真实验与性能分析可知,虽然MCML算法的计算复杂度略高于文献[14]算法,但其估计精度却大幅度提高,并且信噪比在0 dB到25 dB的条件下,均能始终逼近克拉美罗界。当信噪比为25 dB时,MCML算法的时延估计分布范围缩小为MCMC算法的34%。

5 总 结

在定位系统中,最大似然时延估计的均方误差具有逼近克拉美罗界的性能,但是计算复杂度较大。针对该问题,给出了一种基于蒙特卡罗的最大似然时延估计算法。并给出了标准化概率密度函数、随机抽样方法、模型的克拉美罗界以及计算复杂度分析。算法利用直接抽样蒙特卡罗方法将统计平均代替积分,实现了最大似然时延估计,从而使估计结果均方误差始终逼近克拉美罗界。

[1] 王云龙, 吴瑛. 联合时延与多普勒频率的直接定位改进算法 [J]. 西安交通大学学报, 2015, 49(4): 1-7. WANG Yunlong, WU Ying. An improved direct localization algorithm with combined time delay and Doppler [J]. Journal of Xi’an Jiaotong University, 2015, 49(4): 1-7.

[2] REN H P, LI W C, LIU D. Hopf bifurcation analysis of Chen circuit with direct time delay feedback [J]. Chinese Physics: B, 2010, 19(3): 030511.

[3] LUO B W, DONG J J, YU Y, et al. Novel loop-like solitons for generalized Vakhnenko equation [J]. Chinese Physics: B, 2013, 22(3): 023201.

[4] 李晶, 裴亮, 曹茂永, 等. 一种用于多径环境的超分辨率TOA定位算法 [J]. 电波科学学报, 2006, 21(5): 771-776. LI Jing, PEI Liang, CAO Maoyong, et al. Super-resolution TOA algorithm in multi-path environments [J]. Chinese Journal of Radio Science, 2006, 21(5): 771-776.

[5] LI X, MA X, YAN S, et al. Super-resolution time delay estimation for narrowband signal [J]. IET Radar Sonar & Navigation, 2012, 6(8): 781-787.

[6] LI X, PAHLAVAN K. Super-resolution TOA estimation with diversity for indoor geolocation [J]. IEEE Transactions on Wireless Communication, 2004, 3(1): 224-234.

[7] JIANG Hong, CAO Fucheng, DING Rui. Propagator method-based TOA estimation for UWB indoor environment in the presence of correlated fading amplitudes [C]∥Proceedings of IEEE International Conference on Circuits and Systems for Communications. Piscataway, NJ, USA: IEEE, 2008: 535-538.

[8] 王方秋, 张小飞, 汪飞. IR-UWB系统中基于root-MUSIC算法的TOA和DOA联合估计 [J]. 通信学报, 2014, 35(2): 137-145. WANG Fangqiu, ZHANG Xiaofei, WANG Fei. Root-MUSIC-based joint TOA and DOA estimation in IR-UWB [J]. Journal on Communications, 2014, 35(2): 137-145.

[9] OH D, KIM S, YOON S H, et al. Two-dimensional ESPRIT-like shift-invariant TOA estimation algorithm using multi-band chirp signals robust to carrier frequency offset [J]. IEEE Transactions on Wireless Communication, 2013, 12(7): 3130-3139.

[10]李树, 李刚, 田东风, 等. 热辐射输运问题的隐式蒙特卡罗方法求解 [J]. 物理学报, 2013, 62(24): 249501. LI Shu, LI Gang, TIAN Dongfeng, et al. An implicit Monte Carlo method for thermal radiation transport [J]. Acta Physica Sinica, 2013, 62(24): 249501.

[11]NG W, REILLY J P, KIRUBARAJAN T, et al. Wideband array signal processing using MCMC methods [J]. IEEE Transactions on Signal Processing, 2005, 53(2): 411-426.

[12]MASMOUDI A, BELLILI F, AFFES S, et al. A non-data-aided maximum likelihood time delay estimator using importance sampling [J]. IEEE Transactions on Signal Processing, 2011, 59(10): 4505-4515.

[13]MASMOUDI A, BELLILI F, AFFES S, et al. A maximum likelihood time delay estimator in a multipath environment using importance sampling [J]. IEEE Transactions on Signal Processing, 2013, 61(1): 182-193.

[14]李晶, 赵拥军, 李东海. 基于马尔科夫链蒙特卡罗的时延估计算法 [J]. 物理学报, 2014, 63(13): 130701. LI Jing, ZHAO Yongjun, LI Donghai. Time delay estimation using Markov chain Monte Carlo method [J]. Acta Physica Sinica, 2014, 63(13): 130701.

[15]PINCUS M. A closed form solution for certain programming problems [J]. Operations Research, 1968, 16(3): 690-694.

[16]STOICA P, ARYE N. MUSIC, maximum likelihood, and Cramer-Rao bound [J]. IEEE Transactions on Acoustics, Speech, Signal Processing, 1989, 37(5): 720-741.

[本刊相关文献链接]

王云龙,吴瑛.联合时延与多普勒频率的直接定位改进算法.2015,49(4):123-129.[doi:10.7652/xjtuxb201504020]

熊涛,江桦,崔鹏辉,等.应用基扩展模型的混合信号单通道盲分离算法.2015,49(6):60-66.[doi:10.7652/xjtuxb201506 010]

吴一全,孟天亮,吴诗婳.人工蜂群优化的非下采样Shearlet域引导滤波图像增强.2015,49(6):39-45.[doi:10.7652/xjtuxb201506007]

郝雯洁,齐春.一种鲁棒的稀疏信号重构算法.2015,49(4):98-103.[doi:10.7652/xjtuxb201504016]

屈鉴铭,刘志镜,贺文骅.结合有向场景运动模式的粒子滤波行人跟踪方法.2014,48(12):74-79.[doi:10.7652/xjtuxb 201412012]

尚佳栋,王祖林,周丽娜,等.采用随机共振增强的混合扩频信号跳频参数估计.2014,48(10):42-48.[doi:10.7652/xjtuxb201410007]

储颖,牟轩沁,洪伟.采用形状一致性特征的盲图像质量评价方法.2014,48(8):12-17.[doi:10.7652/xjtuxb201408003]

吴仁斌,贾维敏,等.采用幅度响应约束的鲁棒自适应波束形成算法.2014,48(4):109-114.[doi:10.7652/xjtuxb201404 019]

郝红侠,刘芳,焦李成,等.采用结构自适应窗的非局部均值图像去噪算法.2013,47(12):71-76.[doi:10.7652/xjtuxb 201312013]

刁瑞朋,孟庆丰.一种衰减信号加窗频域插值算法.2013,47(7):85-90.[doi:10.7652/xjtuxb201307016]

穆为磊,高建民,王昭,等.考虑人眼视觉特性的射线检测数字图像质量评价方法.2013,47(7):91-95.[doi:10.7652/xjtuxb201307017]

郑纪彬,任爱锋,苏涛,等.多分量Chirp信号相位参数的精确估计算法.2013,47(2):69-74.[doi:10.7652/xjtuxb2013 02012]

吴俊峰,牟轩沁,张砚博.一种快速迭代软阈值稀疏角CT重建算法.2012,46(12):24-29.[doi:10.7652/xjtuxb201212 005]

杨旭,李兵兵,常义林.引入拉格朗日算子的纹理深度联合比特分配方法.2012,46(12):42-48.[doi:10.7652/xjtuxb 201212008]

朱磊,水鹏朗,章为川,等.利用区域划分的合成孔径雷达图像相干斑抑制算法.2012,46(10):83-88.[doi:10.7652/xjtuxb201210015]

谢航,廖与禾,林京.一种含噪声多模态调频调幅信号的提取方法.2012,46(7):93-97.[doi:10.7652/xjtuxb201207017]

(编辑 刘杨)

A Maximum Likelihood Time Delay Estimation Algorithm Using Monte Carlo Method

BA Bin,ZHENG Na’e,ZHU Shilei,HU Hanying

(Institute of Navigation and Aerospace Engineering, Information Engineering University, Zhengzhou 450001, China)

A maximum likelihood time delay estimation algorithm using Monte Carlo method (MCML) is proposed to solve the problems that the maximum likelihood delay algorithm has high computational complexity due to peak searching and is easy to fall into local convergence. A likelihood function is constructed by using the channel response estimation vector in frequency domain. Then the MCML translates the time delay estimation into the expectation of a random variable, and a standardization probability density function is built from the index likelihood function to approximate an impulse function, and to make the variance of the random variable approach zero. Finally, the random variable is sampled using the Monte Carlo method, and the time delay is estimated from sampling mean. Compared with the traditional methods, the MCML avoids the grid search, reduces the computational complexity, and ensures the global convergence and estimation accuracy. Simulation results show that the MCML is always close to Cramer-Rao bound, and the time delay estimation range of the MCML is 34% of the MCMC’s range when the signal to noise ratio is from 0 dB to 25 dB.

maximum likelihood; Monte Carlo; time delay estimation; Cramer-Rao bound

2015-01-16。 作者简介:巴斌(1987—),男,博士生;胡捍英(通信作者),男,教授,博士生导师。 基金项目:国家高技术研究发展计划资助项目(2012AA01A502,2012AA01A505)。

时间:2015-07-22

10.7652/xjtuxb201508005

TN911.7

A

0253-987X(2015)08-0024-07

网络出版地址:http:∥www.cnki.net/kcms/detail/61.1069.T.20150722.1638.001.html