基于馬爾科夫鍵蒙特卡洛抽樣的最大似然時差-頻差聯合估計算法

趙擁軍 趙勇勝 趙 闖

?

基于馬爾科夫鍵蒙特卡洛抽樣的最大似然時差-頻差聯合估計算法

趙擁軍*趙勇勝 趙 闖

(解放軍信息工程大學導航與空天目標工程學院 鄭州 450001)

該文針對無源定位中參考信號真實值未知的時差-頻差聯合估計問題，構建了一種新的時差-頻差最大似然估計模型，并采用馬爾科夫鏈蒙特卡洛(MCMC)方法求解似然函數的全局極大值，得到時差-頻差聯合估計。算法通過生成時差-頻差樣本，并統計樣本均值得到估計值，克服了傳統互模糊函數(CAF)算法只能得到時域和頻域采樣間隔整數倍估計值的問題，且不存在期望最大化(EM)等迭代算法的初值依賴和收斂問題。推導了時差-頻差聯合估計的克拉美羅界，并通過仿真實驗表明，算法在不同信噪比條件下的估計精度優于CAF算法和EM算法，且計算復雜度較低。

無源定位；時差；頻差；聯合估計；最大似然；馬爾科夫鏈蒙特卡洛方法

1 引言

無源定位是近年來備受關注的問題，在雷達[1]、聲吶[2]、無線通信[3]、傳感器網絡[4]等領域應用廣泛。而時差(Time Difference Of Arrival, TDOA)和頻差(Frequency Difference Of Arrival, FDOA)作為無源定位所需的基本參數[5]，其估計精度將直接決定對目標的定位精度。因此，研究高精度的時差-頻差估計算法具有重要意義。

互模糊函數(Cross Ambiguity Function, CAF)方法是處理時差-頻差聯合估計的經典算法[6]，本質是時差-頻差的2維相關。在高信噪比和高采樣率條件下，互模糊函數方法可以得到時差-頻差的精確估計，但其需要在參數空間上進行網格搜索求解，效率較低，且只能得到時域和頻域采樣間隔整數倍的時差-頻差估計值。為此，一些針對互模糊函數計算的快速算法被提出，如基于快速傅里葉變換，分數階傅里葉變換，兩步估計等的改進算法。這些算法一定程度上減小互模糊函數的計算量。此外，基于高階累積量[11]，小波變換[12]，以及自適應算法也被提出，在一些特定情況得到了優于傳統CAF算法的估計效果。但本質上，這些改進算法仍然是時差-頻差的2維相關，其估計精度仍受到時域和頻域采樣間隔的限制。為此，文獻[13]提出了基于期望最大化(Expection Maximum, EM)的時差-頻差估計算法。EM算法不受采樣間隔的限制，但由于其求解過程中需多次對矩陣求逆，計算量過大，限制了信號的實時處理，且作為一種迭代算法，EM算法存在初值依賴和局部收斂問題。

馬爾科夫鏈蒙特卡洛方法(Markov Chain Monte Carlo, MCMC)方法作為一種重要的Monte Carlo方法，是以概率統計理論為指導的，用隨機數抽樣來解決參數估計問題的一類數值計算方法，因其較高的靈活性，以及在復雜高維、高度非線性等問題中表現出的優異性能[14]，近年來在參數估計領域中得到了廣泛的應用。文獻[15]利用MCMC方法估計降水-徑流模型中的未知參數；文獻[16]采用MCMC方法實現合成孔徑雷達中運動目標的線性調頻(chirp)回波信號的參數估計。文獻[17]針對陣列信號測角問題，通過引入可逆跳轉馬爾科夫鏈蒙特卡羅(Reversible Jump Markov Chain Monte Carlo, RJMCMC)方法實現了真正意義上的信號源數和到達角度的聯合估計。文獻[18]將MCMC方法應用到時延估計問題中，得到了優于傳統算法的估計性能。這類算法的思想是建立待估參數的概率模型或隨機過程，然后利用MCMC方法對概率模型或隨機過程抽樣，通過對樣本的統計實現參數的估計，具有估計精度高、計算復雜度低的優點。

本文針對無源定位中參考信號真實值未知的時差-頻差聯合估計問題，構建了一種新的時差-頻差最大似然估計模型，并采用MCMC方法求解似然函數的全局極大值，得到時差-頻差估計。MCMC方法通過生成時差-頻差樣本，并對樣本進行統計得到估計值，可以突破傳統算法只能得到采樣間隔整數倍的限制，具有較高的估計精度，且計算復雜度較低。

2 信號模型和似然函數

考慮如圖1所示的無源雙基地雷達系統。參考天線用于接收來自外輻射源的直達信號，監視天線用于接收目標回波信號[1]。

圖1 無源雙基地雷達系統配置

對式(10)中的概率密度函數取對數并去掉常數項，得到對數似然函數為

3 時差-頻差聯合估計

文獻[19]提出的處理非線性問題全局最優解的基本理論為尋找一種無需搜索且不存在初值依賴的全局最優算法奠定了基礎。根據文獻[19]的理論，使得對數似然函數取得全局最大值的變量可以通過式(15)得到

式(17)中的積分難以直接通過解析方法計算。但是如果能夠得到足夠多服從分布的樣本，則可以通過數值計算方法估計式(17)中的積分。假設已經得到了個的樣本，那么式(17)中積分可通過統計樣本均值近似得到。

3.1 MCMC方法

MCMC方法的第1個“MC”, Markov Chain，表示利用Markov Chain從目標分布中抽取樣本，第2個“MC”, Monte Carlo，則表示在抽取的樣本下利用Monte Carlo方法對積分進行計算。它的基本思想是通過建立一個平穩分布為的馬爾可夫鏈來得到服從分布的樣本，然后通過對樣本的統計來估計參數值。對于本文最大似然估計模型，平穩分布即為的后驗分布。

Metropolis-Hasting(M-H)抽樣是構造馬爾科夫鏈的常用方法。MCMC方法的精髓在于構造合適的馬爾科夫鏈，因此算法的主要目的是對馬爾科夫鏈，在給定一個所處狀態下，產生下一步的狀態。M-H算法構造馬爾科夫鏈的主要步驟總結如下：

(3)重復(直至馬氏鏈達到平穩狀態)：

為此，本文采用自適應的隨機游走采樣方法，自適應地控制增量方差的大小，使其隨著抽樣次數的增加取值不斷減小，即游走的范圍不斷縮小。在抽取第個樣本時，

3.2 基于MCMC的時差-頻差估計

同時為避免式中指數運算的數值過大，令

那么最終構建馬爾科夫鏈的平穩分布函數為

圖2 對平穩分布函數形狀的影響

綜上，利用MCMC方法進行時差-頻差聯合估計的具體實現過程總結如下：

4 CRLB分析

式中，

那么，Fisher信息矩陣為

CRLB是無偏算法估計誤差的理論下限，其等于Fisher信息矩陣的逆。那么時差和頻差的估計誤差方差滿足式(32)中的不等式

5 仿真實驗

選取一段BPSK信號作為輻射源信號，進行仿真實驗。BPSK信號參數設置為：碼元速率，信號載頻。采樣頻率，信號快拍數，兩路信號之間的時差，頻差。信號的信噪比初步設置為10 dB, MCMC方法的參數初步設置為,,,。

圖3給出了信號信噪比為10 dB時利用MCMC方法抽取的時差頻差樣本。可以看出，抽樣過程開始后，樣本很快收斂至平穩分布，然后圍繞著時差頻差真實值上下波動。統計不同樣本數量下MCMC算法估計的均方根誤差(Root Mean Square Error, RMSE)，結果如圖4所示。可以看出，隨著樣本數量的增加，時差和頻差的估計精度均不斷提高，但提高的速度變慢，在樣本數量增加至2000后，基本不再提高。且樣本數量的增加意味著計算復雜度的增加，因此，作為折中，在后續仿真中樣本數量設置為。

圖3 MCMC方法抽取樣本圖????????圖4 不同樣本數量下算法的RMSE

為評估算法估計性能，在不同信噪比條件下利用算法進行蒙特卡羅仿真。算法的估計誤差為1000次仿真的RMSE。為了突出本文基于MCMC的ML(MCMC-ML)算法性能，將算法的RMSE與基于FFT的CAF(FFT-CAF)算法[10]、EM算法[13]和CRLB對比。仿真結果如圖8所示。

從圖8(a)可以看出，隨著信噪比的增加，幾種算法的時差估計精度均有提高。但FFT-CAF算法在信噪比增加至5 dB后，估計精度基本保持不變，不再隨信噪比的增加而提高。原因在于FFT-CAF只能得到時域和頻域采樣間隔整數倍的時差-頻差估計，估計精度受到限制。EM算法和MCMC-ML算法均可得到頻域和時域采樣間隔非整數倍的時差-頻差估計，因此在-5 dB至20 dB信噪比范圍內，估計精度可隨著信噪比的增加而不斷提高，但EM算法的估計精度對初值依賴嚴重。初值較差時，EM算法的估計精度低于CAF算法。而在給定較為準確的初值時，EM算法的估計精度高于CAF算法，較高信噪比條件下估計精度與MCMC-ML算法相近，但在信噪比較低時的估計精度低于本文MCMC-ML算法。圖8(b)表明，MCMC-ML算法在頻差估計中性能同樣優于FFT-CAF算法和EM算法，但與時差估計相比不同的是，幾種算法對頻差估計的精度均相對較高，這主要由信號的互模糊特性決定。

算法的計算復雜度也是衡量算法優劣的重要指標。為此，這里比較本文MCMC-ML算法，基于網格搜索的ML(Grid search-ML)算法，FFT-CAF算法及EM算法的計算復雜度。由于實際運算過程中運算量主要由復數乘法運算次數決定，因此將算法復數乘法的次數作為衡量算法計算復雜度的指標。為便于統計，這里將共軛運算和指數運算均作為一次復數乘法運算參與統計。結果如表1所示。其中，分別為Grid search-ML算法和FFT-CAF算法在時差和頻差取值區間劃分點數。為MCMC算法的樣本數。為EM算法的迭代次數。

圖5 對算法性能的影響???圖6 對算法性能的影響???圖7 對算法性能的影響

圖8 不同信噪比條件下算法的估計誤差

從表1可以看出，4種算法中，Grid search-ML算法的計算計算復雜度最高，難以滿足實時處理的要求。而與之相比，MCMC-ML算法的計算復雜度很低，與FFT-CAF算法的計算復雜度相當。EM算法計算復雜度較高，原因在于EM算法在期望最大化的迭代過程中需多次對的矩陣求逆，造成算法計算復雜度的急劇增加。從計算復雜度的表達式可以看出，本文MCMC-ML算法的計算復雜度主要由信號快拍數和樣本數量決定，計算復雜度隨著生成樣本數的增加而增加。對于仿真BPSK信號情況，MCMC-ML算法的計算復雜度低于FFT-CAF快速計算方法，可以滿足信號實時處理的要求。

6 結論

針對無源雙基地定位中參考信號真實值未知的時差-頻差聯合估計問題，本文構建了一種新的時

差-頻差最大似然估計模型，并采用MCMC方法求解最大似然模型，得到時差-頻差估計。MCMC方法通過生成時差-頻差的樣本，進而通過統計樣本均值得到時差-頻差估計。算法的計算復雜度與利用FFT的CAF快速計算方法基本相同，但是克服了CAF算法只能得到時域和頻域采樣間隔整數倍的時差-頻差估計問題，可以得到采樣間隔非整數倍的時差-頻差估計，因此估計精度高于CAF算法。而與EM算法相比，本文算法不存在迭代的初值依賴和收斂問題，且計算復雜度遠低于EM算法。推導了時差-頻差聯合估計的CRLB，并通過仿真實驗表明，算法的估計精度優于CAF算法和EM算法。

表1不同算法的計算復雜度對比

算法計算復雜度BPSK信號計算復雜度比 MCMC-ML 1.000 Grid search-ML 2095.900 FFT-CAF 1.102 EM 808.070

[1] HIGGINS T, WEBSTER T, and MOKOLE E L. Passive multistatic radar experiment using WiMAX signals of opportunity. Part 1: Signal processing[J].,&, 2016, 10(2): 238-247. doi: 10.1049/iet-rsn. 2015.0020.

[2] LI Ruiyang and HO K. Efficient closed-form estimators for multistatic sonar localization[J]., 2015, 51(1): 600-614. doi: 10.1109/TAES.2014.140482.

[3] ZEMMARI R, BROETJE M, BATTISTELLO G,. GSM passive coherent location system: Performance prediction and measurement evaluation[J].,&, 2014, 8(2): 94-105. doi: 10.1049/iet-rsn.2013.0206.

[4] DECARLI N, GUIDI F, and DARDARI D. A novel joint RFID and radar sensor network for passive localization: Design and performance bounds[J]., 2014, 8(1): 80-95. doi: 10.1109 /JSTSP.2013.2287174.

[5] 曲付勇, 孟祥偉. 基于約束總體最小二乘方法的到達時差到達頻差無源定位算法[J]. 電子與信息學報, 2014, 36(5): 1075-1081. doi: 10.3724/SP.J.1146.2013.01019.

QU Fuyong and MENG Xiangwei. Source localization using TDOA and FDOA measurements based on constrained total least squares algorithm[J].&, 2014, 36(5): 1075-1081. doi: 10.3724 /SP.J.1146.2013.01019.

[6] STEIN S. Algorithms for ambiguity function processing[J].,,, 1981, 29(3): 588-599. doi: 10.1109/TASSP. 1981.1163621.

[7] TOLIMIERI R and WINOGRAD S. Computing the ambiguity surface[J].,,, 1985, 33(5): 1239-1245. doi: 10.1109/ TASSP.1985.1164688.

[8] AUSLANDER L and TOLIMIERI R. Computing decimated finite cross-ambiguity functions[J].,,, 1988, 36(3): 359-364. doi: 10.1109/29.1532.

[9] OZDEMIR A K and ARIKAN O. Fast computation of the ambiguity function and the Wigner distribution on arbitrary line segments[J]., 2001, 49(2): 381-393. doi: 10.1109/78.902121.

[10] TAO R, ZHANG W Q, and CHEN E Q. Two-stage method for joint time delay and Doppler shift estimation[J].,&, 2008, 2(1): 71-77. doi: 10.1049 /iet-rsn:20060014.

[11] SHIN D C and NIKIAS C L. Complex ambiguity functions using nonstationary higher order cumulant estimates[J]., 1995, 43(11): 2649-2664. doi: 10.1109/78.482115.

[12] NIU X, CHING P C, and CHAN Y T. Wavelet based approach for joint time delay and Doppler stretch measurements[J].and, 1999, 35(3): 1111-1119. doi: 10.1109/7. 784079.

[13] BELANGER S P. Multipath TDOA and FDOA estimation using the EM algorithm[C]. IEEE International Conference on Acoustics, Speech, and Signal Processing, Minneapolis, USA, 1993: 168-171. doi: 10.1109/ICASSP.1993.319621.

[14] GILAVERT C, MOUSSAOUI S, and IDIER J. Efficient Gaussian sampling for solving large-scale inverse problems using MCMC[J]., 2015, 63(1): 70-80. doi: 10.1109/TSP.2014.2367457.

[15] BATES B C and CAMPBEL E P. A Markov chain Monte Carlo scheme for parameter estimation and inference in conceptual rainfall-runoff modeling[J]., 2001, 37(4): 937-947. doi: 10.1029/2000WR900363.

[16] 林彥, 王秀壇, 彭應寧, 等. 基于MCMC的線性調頻信號最大似然參數估計[J]. 清華大學學報(自然科學版), 2004, 44(4): 511-514. doi: 10.3321/j.issn:1000-0054.2004.04.020.

LIN Yan, WANG Xiutan, PENG Yingning,. Maximum likelihood parameter estimation of chirp signals based on MCMC[J].(), 2004, 44(4): 511-514. doi: 10.3321/j.issn:1000- 0054.2004.04.020.

[17] NG W, REILLY J P, KIRUBARAJAN T,. Wideband array signal processing using MCMC methods[J]., 2005, 53(2): 411-426. doi: 10.1109/TSP.2004.838934.

[18] 李晶, 趙擁軍, 李冬海. 基于馬爾科夫鏈蒙特卡羅的時延估計算法[J]. 物理學報, 2014, 63(13): 67-73. doi: 10.7498/aps.63. 130701.

LI Jing, ZHAO Yongjun, and LI Donghai. Time delay estimation using Markov chain Monte Carlo method[J]., 2014, 63(13): 67-73. doi: 10.7498/aps.63. 130701.

[19] PINCUS M. A closed form solution of certain programming problems[J]., 1968, 16(3): 690-694. doi: 10.1287/opre.16.3.690.

Maximum Likelihood TDOA-FDOA Estimator Using Markov Chain Monte Carlo Sampling

ZHAO Yongjun ZHAO Yongsheng ZHAO Chuang

(,,450001,)

This paper investigates 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

ignal is unknown. A novel Maximum Likelihood (ML) estimator of TDOA and FDOA is constructed, and Markov Chain Monte Carlo (MCMC) method is applied to finding the global maximum of likelihood function by generating the realizations of TDOA and FDOA. Unlike the Cross Ambiguity Function (CAF) algorithm or the Expectation Maximization (EM) algorithm, the proposed algorithm can also 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 for different SNR conditions with higher accuracy and lower computational complexity.

Passive location; Time Difference Of Arrival (TDOA); Frequency Difference Of Arrival (FDOA); Joint estimation; Maximum Likelihood (ML); Markov Chain Monte Carlo (MCMC) method

TN971

A

1009-5896(2016)11-2745-08

10.11999/JEIT160050

2016-01-13；改回日期：2016-06-08；

2016-09-01

趙擁軍 zhaoyongjuntg@126.com

國家自然科學基金(61401469, 41301481, 61501513)，國家高技術研究發展計劃(2012AA7031015)

The National Natural Science Foundation of China (61401469, 41301481, 61501513), The National High Technology Research and Development Program of China (2012AA7031015)

趙擁軍： 男，1964年生，教授，博士生導師，主要研究方向為雷達信號處理、陣列信號處理.

趙勇勝： 男，1990年生，碩士生，研究方向為無源定位.

趙 闖： 男，1978年生，教授，主要研究方向為雷達信號處理.