利用重要性采樣的時差－頻差聯合估計算法

趙勇勝，趙擁軍，趙闖

利用重要性采樣的時差－頻差聯合估計算法

趙勇勝，趙擁軍＊，趙闖

解放軍信息工程大學 導航與空天目標工程學院，鄭州 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