趙擁軍 趙勇勝 趙 闖

?

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

趙擁軍*趙勇勝 趙 闖

(解放軍信息工程大學(xué)導(dǎo)航與空天目標(biāo)工程學(xué)院 鄭州 450001)

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

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

1 引言

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

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

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

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

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

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

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

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

3 時(shí)差-頻差聯(lián)合估計(jì)

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

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

3.1 MCMC方法

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

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

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

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

3.2 基于MCMC的時(shí)差-頻差估計(jì)

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

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

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

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

4 CRLB分析

式中，

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

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

5 仿真實(shí)驗(yàn)

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

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

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

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

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

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

圖5 對(duì)算法性能的影響???圖6 對(duì)算法性能的影響???圖7 對(duì)算法性能的影響

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

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

6 結(jié)論

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

差-頻差最大似然估計(jì)模型，并采用MCMC方法求解最大似然模型，得到時(shí)差-頻差估計(jì)。MCMC方法通過(guò)生成時(shí)差-頻差的樣本，進(jìn)而通過(guò)統(tǒng)計(jì)樣本均值得到時(shí)差-頻差估計(jì)。算法的計(jì)算復(fù)雜度與利用FFT的CAF快速計(jì)算方法基本相同，但是克服了CAF算法只能得到時(shí)域和頻域采樣間隔整數(shù)倍的時(shí)差-頻差估計(jì)問(wèn)題，可以得到采樣間隔非整數(shù)倍的時(shí)差-頻差估計(jì)，因此估計(jì)精度高于CAF算法。而與EM算法相比，本文算法不存在迭代的初值依賴(lài)和收斂問(wèn)題，且計(jì)算復(fù)雜度遠(yuǎn)低于EM算法。推導(dǎo)了時(shí)差-頻差聯(lián)合估計(jì)的CRLB，并通過(guò)仿真實(shí)驗(yàn)表明，算法的估計(jì)精度優(yōu)于CAF算法和EM算法。

表1不同算法的計(jì)算復(fù)雜度對(duì)比

算法計(jì)算復(fù)雜度BPSK信號(hào)計(jì)算復(fù)雜度比 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] 曲付勇, 孟祥偉. 基于約束總體最小二乘方法的到達(dá)時(shí)差到達(dá)頻差無(wú)源定位算法[J]. 電子與信息學(xué)報(bào), 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] 林彥, 王秀壇, 彭應(yīng)寧, 等. 基于MCMC的線性調(diào)頻信號(hào)最大似然參數(shù)估計(jì)[J]. 清華大學(xué)學(xué)報(bào)(自然科學(xué)版), 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] 李晶, 趙擁軍, 李冬海. 基于馬爾科夫鏈蒙特卡羅的時(shí)延估計(jì)算法[J]. 物理學(xué)報(bào), 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

國(guó)家自然科學(xué)基金(61401469, 41301481, 61501513)，國(guó)家高技術(shù)研究發(fā)展計(jì)劃(2012AA7031015)

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

趙擁軍： 男，1964年生，教授，博士生導(dǎo)師，主要研究方向?yàn)槔走_(dá)信號(hào)處理、陣列信號(hào)處理.

趙勇勝： 男，1990年生，碩士生，研究方向?yàn)闊o(wú)源定位.

趙 闖： 男，1978年生，教授，主要研究方向?yàn)槔走_(dá)信號(hào)處理.