利用重要性采样的时差－频差联合估计算法

赵勇胜，赵拥军，赵闯

利用重要性采样的时差－频差联合估计算法

赵勇胜，赵拥军＊，赵闯

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

针对无源定位中参考信号真实值未知的时差（TDOA）－频差（FDOA）联合估计问题，构建了一种新的时差－频差最大似然（ML）估计模型，并采用重要性采样（IS）方法求解似然函数极大值，得到时差－频差联合估计。算法通过生成时差－频差样本，并统计样本加权均值得到估计值，克服了传统互模糊函数（CAF）算法只能得到时域和频域采样间隔整数倍估计值的问题，且不存在期望最大化（EM）等迭代算法的初值依赖和收敛问题。推导了时差－频差联合估计的克拉美罗下界（CRLB），并通过仿真实验表明，算法的计算复杂度适中，估计精度优于CAF算法和EM算法，在不同信噪比条件下估计误差接近CRLB。

时差；频差；联合估计；最大似然；重要性采样

无源定位是近年来备受关注的问题，在雷达［1］、声呐［2］、无线通信［3］以及传感器网络［4］等领域应用广泛。而时差（Time Difference of Arrival，TDOA）和频差（Frequency Difference of Arrival，FDOA）作为定位所需的基本参数［5］，其估计精度将直接决定对目标的定位精度。因此，研究高精度的时差－频差估计算法具有重要意义。

互 模 糊 函 数 （Cross Ambiguity Function，CAF）是处理时差－频差联合估计的经典算法［6］，本质是时差－频差的二维相关。在高信噪比和高采样率条件下，互模糊函数方法可以得到时差－频差的精确估计，但其需要在参数空间上进行网格搜索求解，效率较低，且只能得到时域和频域采样间隔整数倍的时差－频差估计值。为此，一些针对互模糊函数计算的快速算法被提出，如基于快速傅里叶变换、分数阶傅里叶变换、两步估计等的改进算法［7－10］。这些算法一定程度上减小了互模糊函数的计算量。此外，基于高阶累积量［11］、小波变换［12］以及自适应算法也被提出，在一些特定情况得到了优于传统CAF算法的估计效果。但本质上，这些改进算法仍然是时差－频差的二维相关，其估计精度仍受到时域和频域采样间隔的限制。为此，文献［13］提出了基于期望最大化（Expection Maximum，EM）的时差－频差估计算法。EM算法不受采样间隔的限制，但由于其求解过程中需多次对矩阵求逆，其计算量过大，限制了信号的实时处理，且作为一种迭代算法，EM算法存在初值依赖和局部收敛问题。

近年来，信号处理领域的一些新算法为解决时差－频差联合估计问题提供了新思路。重要性采样（Importance Sampling，IS）算法作为一种重要的Monte Carlo方法，是以概率统计理论为指导的，用随机数抽样来解决参数估计问题的一类数值计算方法，其基本思想是通过一个相对简单的分布函数的随机加权平均来近似计算目标分布函数的数学期望［14］。文献［15］采用IS算法估计均匀线阵中信号到达角度，得到了优于MUSIC算法和最小范数算法的估计精度。文献［16］考虑均匀线阵中信号到达角度和多普勒频率估计问题，利用IS方法实现了角度－多普勒频率的联合估计。文献［17－18］将IS算法应用到参考信号已知条件下的主动时延估计问题中。这类算法的思想是建立待估参数的概率模型或随机过程，然后利用IS算法对概率模型或随机过程抽样，通过对样本的统计实现参数的估计，对于解决非线性的参数估计问题具有很好的实践意义。

本文针对无源定位中参考信号真实值未知的时差－频差联合估计问题，构建了一种新的最大似然估计模型，并采用IS算法求解似然函数极大值，得到时差－频差估计。IS算法通过生成时差－频差样本，并对样本进行统计得到估计值，可以突破采样间隔整数倍的限制，具有较高的估计精度，且计算复杂度适中。

1 信号模型和似然函数

考虑如图1所示的无源双基地雷达系统。参考天线用于接收来自外辐射源的直达信号，监视天线用于接收目标回波信号［1］。

假设参考天线接收到的信号真实值为s（t），则两路接收机接收到的信号可建模为［10］

式中：w1（t）、w2（t）为零均值的高斯白噪声；a 和φ为冗余参数，分别为增益系数和相移；τ和ν为待估参数，分别为两路信号的时差和频差。本文的主要工作是要通过对x1（t）、x2（t）的观测来估计两路信号的时差和频差。

对x1（t）、x2（t）以间隔Ts采样，得到信号的离散形式为

由于估计过程中参考信号的真实值s（n）是未知的，因此难以直接根据式（2）得到时差和频差的似然函数。为此，从频域出发，推导时差和频差的似然函数。首先对x1（t）、x2（t）做离散傅里叶变换（Discrete Fourier Transform，DFT），得到信号的频域形式为

式中：k＝0，1，…，N－1。

式（4）可以表示为

式中：

由于DFT变换具有周期性，因此m的范围仍为0，1，…，N－1。

将式（3）代入式（5），并整理化简，可以消去未知量S（k），得到仅含有观测量和待估参量的表达式为

式中：

那么，令Xν＝［Xν（0） Xν（1） … Xν（N－1）］T，Xτ＝［Xτ（0） Xτ（1） … Xτ（N－1）］T，W＝［W（0） W（1） … W（N－1）］T，式（7）可以表示为

由于w1（n）、w2（n）为高斯白噪声，经过 DFT等一系列线性变换后，向量W仍满足高斯性。那么根据W 的概率统计特性，给定参数（b，τ，ν）条件下X1、X2的概率密度函数为

式中：σ2为噪声W（k）的方差。

对式（10）中的概率密度函数取对数并去掉与（b，τ，ν）无关的常数项，得到对数似然函数为

对数似然函数L（b，τ，ν）关于τ、ν高度非线性，但却是关于冗余参数b的二次函数。因此可以通过对L（b，τ，ν）关于b求偏导，并令偏导为零，来消除冗余参数b。

对式（12）求解，得到b关于（τ，ν）的表达式为

将式（13）代入式（11），得到仅含有（τ，ν）的对数似然函数为

式中：θ＝［τ ν］T为时差和频差构成的待估向量。

2 时差－频差联合估计算法

Lc（θ）是关于θ的非线性函数，无法利用解析方法得到时差－频差的估计值。常用方法是通过多维网格搜索或利用迭代方法求解。但网格搜索方法效率低，且其只能得到时域和频域采样间隔整数倍的时差和频差。而迭代方法往往需要一个较为准确的初始解，否则难以收敛至全局最优解。为此，本文采用重要性采样的方法估计时差和频差。

2．1 重要性采样方法

根据Pincus的理论［19］，使得代价函数Lc（θ）取得全局最大值的变量θ可以通过式（15）得到

式中：θi为向量θ 中 的 第i 个 变 量；＾θi为 其 估计值。

定义exp［ρLc（θ）］的归一化函数为

式（17）中的积分难以直接通过解析方法计算。但是如果能够得到足够多服从L′c，ρ0（θ）分布的θ样本，则可以通过数值计算方法估计式（17）中的积分。假设已经得到了R个θ的样本，那么式（17）中积分可通过统计样本均值近似得到

IS方法是处理复杂多维积分的有效方法。对于式（17）中的积分，有

式中：g′（θ）为一个容易采样的概率密度函数，称为归一化重要性函数。通过生成服从g′（θ）分布的样本，式（19）中的积分可以利用如下加权平均计算得到：

显然，式（20）的估计效果受样本数量R和选取的归一化重要性函数g′（θ）的影响。增加样本数量R可以提高估计精度，但同时意味着更大的计算量。而实际上，R同时也受到g′（θ）的影响。选取的g′（θ）与L′c，ρ0（θ）越相似，则式（20）的估计方差越小，但同时g′（θ）还应考虑易于采样。因此，重要性函数的选取是利用IS方法估计时差－频差的关键。

2．2 重要性函数选取

首先，根据Lc（θ）初步构建重要性函数为

式中：ρ1为常数，用于调整重要性函数的形状。ρ1的不同取值对应重要性函数形状如图2所示。从图2可以看出，ρ1越大，重要性函数形状越尖锐，反之则越平缓。选择较大的ρ1可以消除重要性函数峰值附近的旁瓣，从而减少采样过程中对参数估计无益的样本数，使得样本集中于峰值附近，提高估计精度。但ρ1并非越大越好，过大的ρ1会使得构造的重要性函数丢失一些有用信息，同时过于尖锐的函数形状也会造成采样困难。ρ1的不同选取参见文献［19］。

归一化重要性函数构建为

但在实际应用中，式（22）中的归一化重要性函数只能采用离散形式

式中：L和K分别为时差和频差取值区间划分的总点数。

2．3 利用重要性采样估计时差－频差

由于选取的g′ρ′1（θ）与L′c，ρ0（θ）存在偏差，因此生成的θ样本的统计特性也是有偏差的，但是这种偏差可以通过 加权因 子 L′c，ρ0（θ）／g′ρ′1（θ）消除。按照g′ρ′1（θ）生成θ的样本，则IS方法得到的时差和频差的估计为

由于时差和频差存在取值范围［θmin，θmax］，因此对于式（26）中的加权平均，采用角度平均［16］方法可以减小计算量，并改善算法在低信噪比和信号快拍数较少条件下的估计效果［20］。首先，将时差和频差归一化至［0，1］内

那么珔θ的角度加权均值为

式中：∠表示取弧度角；珋θ（t）i表示按照g′ρ′1（珔θ）生成的第t个珋θi样本；F（珔θ（t））为加权因子，其表达式为

对于式（28）中的角度均值计算，只需要确定一个复数的相位角，而L′c，ρ0（珔θ）和g′ρ′1（珔θ）中的归一于相位角并无影响，且难以计算，因此将其省略。同时，为了使计算F（珔θ（t））过程中指数运算不至于过大而溢出，将其归一化为

式中：I（珔θ（t））＝I（珋τ（t），珋ν（t））。

综上，利用IS方法估计时差－频差的步骤总结如下：

步骤1 根据接收信号x1（n）、x2（n），按照式（25）构建离散形式的归一化重要性函数

步骤2 按照g′ρ′1（珔θ）生成R个珔θ样本。

步骤3 利用式（30）计算样本的权重，并按照式（28）计算珔θ的角度加权均值。

3 克拉美罗下限（CRLB）分析

假设两路接收机接收到的噪声方差相同，记为var［W1（k）］＝var［W2（k）］＝，那么参考接收机接收到的信号x1（n）的信噪比为

式中：S为参考信号真实值s（n）的DFT变换构成的向量。由式（32）可得

根据噪声的统计特性，W（k）的均值为零，即

对于时差－频差联合估计问题，Fisher信息矩阵为2×2的方阵。根据Fisher信息矩阵的定义，对式（10）中的概率密度函数p（Xν，Xτ；b，θ）取对数，并关于θ求偏导：

式中：

那么，Fisher信息矩阵为

CRLB是无偏算法估计误差的理论下限，其等于Fisher信息矩阵的逆。那么时差和频差的估计误差方差满足：

式中：（J－1）ii为矩阵J－1主对角线上的第i个元素。

4 仿真实验

选取一段二进制相移键控（Binary Phase Shift Keying，BPSK）信号作为辐射源信号，进行仿真实验。BPSK信号参数设置为：码元速率RB＝1MHz，信号载频fc＝50MHz。采样频率Fs＝1／Ts＝128MHz，信号快拍数 N＝2 048，两路信号之间的时差τ＝150．4Ts，频差ν＝100．6Fs／N。IS方法的参数设置为ρ0＝20，ρ1＝15，归一化重要性函数在时差和频差取值区间划分的总点数为L＝K＝2 048。

为了分析样本数量对算法估计精度的影响，利用不同数量的样本进行蒙特卡罗仿真。信号信噪比设置为10dB，仿真结果如图3所示。图3给出了不同样本数量下算法估计的均方误差（Mean Square Error，MSE）。可以看出，随着样本数量的增加，时差和频差的估计精度均不断提高，但提高的速度变慢，在样本数量增加至4 000后，基本不再提高。且样本数量的增加意味着计算复杂度的增加，因此，作为折中，在后续仿真中样本数量设置为R＝4 000。

为评估算法的估计性能，在不同信噪比条件下的利用算法进行蒙特卡罗仿真。为了突出算法的性能，将算法的MSE与基于FFT的CAF算法［10］、EM 算 法［13］和 CRLB 对 比。 仿 真 结 果 如图4所示。从图4（a）可以看出，随着信噪比的增加，几种算法的时差估计精度均有提高。但CAF算法在信噪比增加至5dB后，估计精度基本保持不变，不再随信噪比的增加而提高。原因在于CAF只能得到时域和频域采样间隔整数倍的时差－频差估计，估计精度受到限制。EM算法和IS算法均可得到频域和时域采样间隔非整数倍的时差－频差估计，因此在－5～20dB信噪比范围内，估计精度可随着信噪比的增加而不断提高，但EM算法的估计精度对初值依赖严重。初值较差时，EM算法的估计精度低于CAF算法。而在给定较为准确的初值时，EM算法的估计精度高于CAF算法，较高信噪比条件下估计精度与IS算法相近，但在信噪比较低时估计精度低于本文IS算法。图4（b）表明，IS算法在频差估计中性能同样优于CAF算法和EM算法，但与时差估计相比不同的是，几种算法对频差估计的精度均相对较高，这主要由信号的互模糊特性决定。

此外，图4（b）中，在低信噪比时，算法的估计误差逼近CRLB，而随着信噪比的增加，算法的估计误差虽然降低，但偏离CRLB的程度反而有增大的趋势。原因在于本文IS算法通过生成样本并统计样本加权均值的方式来近似参数的最大似然解。在低信噪比条件下，算法估计误差主要受噪声影响，这种近似引起的误差相比之下很小，可以忽略，而在高信噪比条件下，噪声对算法估计误差的影响较小，这种近似引起的误差相比之下变大，导致算法的估计误差偏离CRLB。

算法的计算复杂度也是衡量算法优劣的重要指标。为此，这里比较基于网格搜索的 ML算法、基于FFT的CAF算法、EM算法及本文算法的计算复杂度。由于实际运算过程中运算量主要由乘法运算决定，因此将算法实数乘法的次数作为衡量算法计算复杂度的指标。为便于统计，这里将共轭运算和指数运算均作为乘法运算参与统计。结果如表1所示。表中：Niter为EM算法的迭代次数。

表1 不同算法的计算复杂度对比Table 1 Computational complexity comparison among different algorithms

从表1可以看出，4种算法中，基于FFT的CAF快速计算算法计算复杂度最低。文献［13］中的EM算法计算复杂度最高，原因在于文献［13］在期望最大化的迭代过程中需多次对N×N的矩阵求逆，造成算法计算复杂度的急剧增加。而基于网格搜索的ML算法的计算效率同样较低，难以满足实时处理的要求。从计算复杂度的表达式可以看出，本文IS算法的计算复杂度主要由构建式（25）中归一化重要性函数的计算量决定，同时会随着生成样本数的增加而有所增加，对于仿真BPSK信号情况，IS算法的计算复杂度略高于基于FFT的CAF快速计算算法，可以满足信号实时处理的要求。

5 结 论

针对无源双基地定位中参考信号真实值未知的时差－频差联合估计问题，提出了一种基于重要性采样的最大似然估计算法。该算法具有如下优势：

1）算法的计算复杂度与利用FFT的CAF快速计算算法基本相同，但是克服了CAF算法只能得到时域和频域采样间隔整数的时差－频差估计问题，可以得到采样间隔非整数倍的时差－频差估计。

2）与EM算法相比，本文算法不存在迭代的初值依赖和收敛问题，且计算复杂度远低于EM算法。

3）推导了时差－频差联合估计的CRLB，并通过仿真实验表明，算法的估计精度优于CAF算法和EM算法，在不同信噪比条件下估计误差接近CRLB。

［1］ LIU J，LI H，HIMED B．On the performance of the cross－correlation detector for passive radar applications［J］．Signal Processing，2015，113：32－37．

［2］ CORALUPPI S．Multistatic sonar localization［J］．IEEE Journal of Oceanic Engineering，2006，31（4）：964－974．

［3］ CAFFERY J J，STUBER G L．Overview of radio location in CDMA cellular systems［J］．IEEE Communications Magazine，1998，16（4）：38－45．

［4］ GEZICI S，ZHI T，GIANNAKIS G B，et al．Localization via ultra－wideband radios：A look at positioning aspects for future sensor networks［J］．IEEE Signal Processing Magazine，2005，22（1）：70－84．

［5］ 刘洋，杨乐，郭福成，等．基于定位误差修正的运动目标TDOA／FDOA无源定位方法［J］．航空学报，2015，36（5）：1617－1626．LIU Y，YANG L，GUO F C，et al．Moving targets TDOA／FDOA passive localization algorithm based on localization error refinement［J］．Acta Aeronautica et Astronautica Sinica，2015，36（5）：1617－1626（in Chinese）．

［6］ STEIN S．Algorithms for ambiguity function processing［J］．IEEE Transactions on Acoustics，Speech，and Signal Processing，1981，29（3）：588－599．

［7］ TOLIMIERI R，WINOGRAD S．Computing the ambiguity surface［J］．IEEE Transactions on Acoustics，Speech，and Signal Processing，1985，33（5）：1239－1245．

［8］ AUSLANDER L，TOLIMIERI R．Computing decimated finite cross－ambiguity functions［J］．IEEE Transactions on Acoustics，Speech，and Signal Processing，1988，36（3）：359－364．

［9］ OZDEMIR A K，ARIKAN O．Fast computation of the ambiguity function and the Wigner distribution on arbitrary line segments［J］．IEEE Transactions on Signal Processing，2001，49（2）：381－393．

［10］ TAO R，ZHANG W Q，CHEN E Q．Two－stage method for joint time delay and Doppler shift estimation［J］．IET Radar，Sonar and Navigation，2008，2（1）：71－77．

［11］ SHIN D C，NIKIAS C L．Complex ambiguity functions using nonstationary higher order cumulant estimates［J］．IEEE Transactions on Signal Processing，1995，43（11）：2649－2664．

［12］ NIU X X，CHING P C，CHAN Y T．Wavelet based approach for joint time delay and Doppler stretch measurements［J］．IEEE Transactions on Aerospace ＆ Electronic Systems，1999，35（3）：1111－1119．

［13］ BELANGER S P．Multipath TDOA and FDOA estimation using the EM algorithm［C］／ICASSP IEEE Computer Society．Piscataway，NJ：IEEE Press，1993：168－171．

［14］ BEICHL I．The importance of importance sampling［J］．Computing in Science ＆Engineering，1999，1（2）：71－73．

［15］ WANG H，KAY S，SAHA S．An importance sampling maximum likelihood direction of arrival estimator［J］．IEEE Transactions on Signal Processing，2008，56（10）：5082－5092．

［16］ WANG H，KAY S．Maximum likelihood angle－Doppler estimator using importance sampling［J］．IEEE Transactions on Aerospace and Electronic Systems，2010，46（2）：610－622．

［17］ MASMOUDI A，BELLILI F，AFFES S，et al．A maximum likelihood time delay estimator using importance sampling［C］／2011IEEE Global Telecommunications Conference（GLOBECOM 2011）．Piscataway，NJ：IEEE Press，2011：1－6．

［18］ 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．

［19］ PINCUS M．A closed form solution of certain programming problems［J］．Operations Research，1968，16（3）：690－694．

［20］ KAY S M．Comments on“Frequency estimation by linear prediction”［J］．IEEE Transactions on Acoustics Speech＆Signal Processing，1979，27（2）：198－199．

TDOA－FDOA joint estimation using importance sampling method

ZHAO Yongsheng，ZHAO Yongjun＊，ZHAO Chuang
School of Navigation and Aerospace Engineering，PLA Information Engineering University，Zhengzhou 450001，China

To solve the joint estimation of time difference of arrival（TDOA）and frequency difference of arrival（FDOA）in passive location system，where the true value of the reference signal is unknown，a novel maximum likelihood（ML）estimator of TDOA and FDOA is constructed．Then importance sampling（IS）method is applied to find the maximum of likelihood function by generating the samples of TDOA and FDOA．Unlike the cross ambiguity function（CAF）algorithm or the expectation maximization（EM）algorithm，the proposed algorithm can estimate the TDOA and FDOA of non－integer multiple of the sampling interval and has no dependence on the initial estimate．The Cramer Rao lower bound（CRLB）is also derived．Simulation results show that the proposed algorithm outperforms the CAF and EM algorithm with higher accuracy and moderate computational complexity，and approaches the CRLB for different SNR conditions．

time difference of arrival；frequency difference of arrival；joint estimation；maximum likelihood；importance sampling

2015－12－31；Revised：2016－03－04；Accepted：2016－03－15；Published online：2016－03－21 13：27

URL：www．cnki．net／kcms／detail／11．1929．V．20160321．1327．008．html

s：National Natural Science Foundation of China（61401469，61501513）；National High Technology Research and Development Program of China（2012AA7031015）

V247；TN971

A

1000－6893（2017）01－319994－08

http：／hkxb．buaa．edu．cn hkxb＠buaa．edu．cn

10．7527／S1000－6893．2016．0085

2015－12－31；退修日期：2016－03－04；录用日期：2016－03－15；网络出版时间：2016－03－21 13：27

www．cnki．net／kcms／detail／11．1929．V．20160321．1327．008．html

国家自然科学基金 （61401469，61501513）；国家“863”计划 （2012AA7031015）

＊通讯作者 ．E－mail：zhaoyongjuntg＠126．com

赵勇胜，赵拥军，赵闯．利用重要性采样的时差－频差联合估计算法［J］．航空学报，2017，38（1）：319994．ZHAO Y S，ZHAO Y J，ZHAO C．TDOA and FDOA joint estimation using importance sampling method［J］．Acta Aeronautica et Astronautica Sinica，2017，38（1）：319994．

（责任编辑：苏磊）

＊Corresponding author．E－mail：zhaoyongjuntg＠126．com