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

趙擁軍 趙勇勝 趙 闖

?

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

趙擁軍*趙勇勝 趙 闖

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

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

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

1 引言

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

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

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

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

2 信號模型和似然函數(shù)

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

圖1 無源雙基地雷達系統(tǒng)配置

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

3 時差-頻差聯(lián)合估計

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

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

3.1 MCMC方法

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

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

(3)重復(直至馬氏鏈達到平穩(wěn)狀態(tài))：

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

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

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

那么最終構(gòu)建馬爾科夫鏈的平穩(wěn)分布函數(shù)為

圖2 對平穩(wěn)分布函數(shù)形狀的影響

綜上，利用MCMC方法進行時差-頻差聯(lián)合估計的具體實現(xiàn)過程總結(jié)如下：

4 CRLB分析

式中，

那么，F(xiàn)isher信息矩陣為

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

5 仿真實驗

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

圖3給出了信號信噪比為10 dB時利用MCMC方法抽取的時差頻差樣本?？梢钥闯觯闃舆^程開始后，樣本很快收斂至平穩(wěn)分布，然后圍繞著時差頻差真實值上下波動。統(tǒng)計不同樣本數(shù)量下MCMC算法估計的均方根誤差(Root Mean Square Error, RMSE)，結(jié)果如圖4所示?？梢钥闯觯S著樣本數(shù)量的增加，時差和頻差的估計精度均不斷提高，但提高的速度變慢，在樣本數(shù)量增加至2000后，基本不再提高。且樣本數(shù)量的增加意味著計算復雜度的增加，因此，作為折中，在后續(xù)仿真中樣本數(shù)量設(shè)置為。

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

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

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

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

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

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

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

6 結(jié)論

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

差-頻差最大似然估計模型，并采用MCMC方法求解最大似然模型，得到時差-頻差估計。MCMC方法通過生成時差-頻差的樣本，進而通過統(tǒng)計樣本均值得到時差-頻差估計。算法的計算復雜度與利用FFT的CAF快速計算方法基本相同，但是克服了CAF算法只能得到時域和頻域采樣間隔整數(shù)倍的時差-頻差估計問題，可以得到采樣間隔非整數(shù)倍的時差-頻差估計，因此估計精度高于CAF算法。而與EM算法相比，本文算法不存在迭代的初值依賴和收斂問題，且計算復雜度遠低于EM算法。推導了時差-頻差聯(lián)合估計的CRLB，并通過仿真實驗表明，算法的估計精度優(yōu)于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的線性調(diào)頻信號最大似然參數(shù)估計[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)，國家高技術(shù)研究發(fā)展計劃(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年生，教授，主要研究方向為雷達信號處理.