相干信源DOA 估計(jì)的最大似然算法應(yīng)用及改進(jìn)

陳佳

(中國(guó)船舶集團(tuán)有限公司第七二三研究所，江蘇 揚(yáng)州 225000)

0 引言

雷達(dá)信號(hào)的波達(dá)方向（Direction of Arrival，DOA）估計(jì)在陣列信號(hào)處理領(lǐng)域中一直是一個(gè)非常重要的研究?jī)?nèi)容，主要通過(guò)提取出空間中按照一定樣式排布的天線陣列接收回波信號(hào)的特征參數(shù)來(lái)估計(jì)目標(biāo)的方位信息，因此國(guó)內(nèi)外學(xué)者提出了各種各樣的DOA 估計(jì)新算法[1-3]。其中最大似然估計(jì)法[4]很早就被提出，但由于其運(yùn)算量較大而沒(méi)有得到廣泛應(yīng)用。本文討論了一種有效的算法，通過(guò)交替投影算法[5-6]把最大似然估計(jì)中多維非線性問(wèn)題轉(zhuǎn)化為一維問(wèn)題來(lái)解決，并在搜索時(shí)采用斐波那契數(shù)列法提高搜索效率，使其運(yùn)算復(fù)雜性降低同時(shí)有效減少了計(jì)算量，滿足工程應(yīng)用的實(shí)時(shí)性要求。

1 陣列信號(hào)處理的數(shù)學(xué)模型[7-9]

假定天線陣列由m 個(gè)陣元組成，后面接有M 個(gè)接收通道，d(d ＜M)個(gè)窄帶遠(yuǎn)場(chǎng)、獨(dú)立的目標(biāo)信號(hào)，其方向角分別為θ1,θ2，…，θd，各個(gè)陣元噪聲為：n1(t)，n2(t)，…，nm(t)，假定噪聲之間相互獨(dú)立、功率相同，且是空間平穩(wěn)的高斯分布白噪聲，第i 個(gè)陣元的輸出為：

式中，sk(t)為第k 個(gè)目標(biāo)回波的信號(hào)；τi(θk)為第k 個(gè)回波信號(hào)到達(dá)第i 個(gè)陣元時(shí)相對(duì)于參考陣元的時(shí)延。有：

上述各個(gè)矢量表示分別為：

得出陣列輸出矢量形式為：

寫成矩陣形式為：

式中，x(t)為陣列的M×1 維快拍數(shù)據(jù)矢量，s(t) 為目標(biāo)信號(hào)的d×1 維數(shù)據(jù)矢量，A(θ)為M×d維陣列流型（導(dǎo)向矢量）矩陣，為信號(hào)向量張成的子空間（信號(hào)子空間），n(t)為陣列M×1 維噪聲數(shù)據(jù)矢量。此式為DOA 估計(jì)的基本數(shù)學(xué)模型表達(dá)式。DOA 估計(jì)任務(wù)就是利用輸出信號(hào)x (t) 估計(jì)出信源個(gè)數(shù)及所在方位θi，i=1,2,…,d。

2 最大似然估計(jì)算法[10-15]

令N 次采樣數(shù)據(jù)為{x1,x2,…,xN}，信號(hào)源數(shù)據(jù)為{s1,s2,…,sN}，噪聲數(shù)據(jù)為{n1,n2,…,nN}。如果噪聲信號(hào){ni}是均值為零，方差為σ2的各態(tài)歷經(jīng)復(fù)高斯過(guò)程，得出觀測(cè)矢量L 次快拍聯(lián)合概率密度函數(shù)為：

式中det[]表示矩陣的行列式，兩邊同時(shí)取負(fù)對(duì)數(shù)為：

求解上式最大似然估計(jì)問(wèn)題，固定θ 和si，（忽略常數(shù)項(xiàng)）得出σ2的最大似然估計(jì)：

將上式代入（11）可得到θ 和si的最大似然估計(jì)：

根據(jù)對(duì)數(shù)函數(shù)的單調(diào)性原則，可將上式等效為：

固定θ,求上式極小得到si的估計(jì)值為：

定義A 的投影陣為：

將（15）代入（16）中可得到確定性最大似然的準(zhǔn)則，即：

可得最大似然準(zhǔn)則以跡的形式表示為：

綜上所述，就可以得到關(guān)于θ 的最大似然估計(jì)器，上式一般用修正的牛頓多維優(yōu)化方法求解。由于在計(jì)算最大似然估計(jì)中需要多維非線性搜索，并且在每次的迭代搜索都要計(jì)算A 的投影矩陣PA，能夠看出這種算法的計(jì)算量相當(dāng)大，因此提出了交替投影迭代（Alternating Projection iterative algorithm，AP）算法應(yīng)用于最大似然估計(jì)算法研究中用以減少計(jì)算量。本文以兩個(gè)信源的情況描述AP 算法，很容易就能夠擴(kuò)展到多個(gè)信源的情況。

AP 算法具體流程如下：

（1）估計(jì)第一個(gè)信源的情形，為：

（2）估計(jì)第二個(gè)信源，假定第一個(gè)信源的方位為θ1

（3）進(jìn)行迭代得：

（4）繼續(xù)迭代得：

（5）重復(fù)（3）和（4）直到收斂。

具體流程如圖1 所示。

圖1 兩維交替投影過(guò)程

采用交替投影法將多維搜索轉(zhuǎn)化到一維搜索上，但是隨著測(cè)角精度提升，搜索次數(shù)也隨之增大，再此基礎(chǔ)上采用一維斐波那契數(shù)列法搜索極值點(diǎn)。

斐波那契數(shù)列法，其基本思想是通過(guò)試探點(diǎn)函數(shù)值比較，使包含極大點(diǎn)的搜索區(qū)間不斷縮小。該方法僅需要計(jì)算函數(shù)值，適用范圍廣，使用方便。

斐波那契數(shù)列法是建立在區(qū)間消去法原理基礎(chǔ)上的試探方法，即在搜索區(qū)間[θmin,θmax]內(nèi)適當(dāng)插入兩點(diǎn)θ1，θ2，并計(jì)算其估計(jì)值。

θ1，θ2將整個(gè)搜索區(qū)間分為三段，應(yīng)用函數(shù)的單峰性質(zhì)，通過(guò)函數(shù)值大小的比較，刪去其中一段，使搜索區(qū)間得以縮小。然后再在保留下來(lái)的區(qū)間上作同樣的處理，如此迭代下去，搜索區(qū)間無(wú)限縮小，從而得到極大點(diǎn)的數(shù)值近似解。

具體步驟如下：

（1）選定初始區(qū)間[θmin,θmax]，及ε＞0，利用式（b1-a1）求出計(jì)算函數(shù)值的次數(shù)n。并設(shè)eps=1×106。接著由下式：

計(jì)算試探點(diǎn)p1和q1。令k=1。

（2）如果f(pk)＜f(qk)，轉(zhuǎn)步驟（3）；否則轉(zhuǎn)步驟（4）。

3 仿真結(jié)果

本文以8 陣元為例，進(jìn)行最大似然算法測(cè)向仿真。均勻線陣，陣元數(shù)為M=8，陣元之間間隔為λ/2，信號(hào)頻率為f=24.3 GHz (λ=12.3 mm)，快拍點(diǎn)數(shù)為1 個(gè)，信源數(shù)為d=2，信噪比為5 dB，設(shè)噪聲服從高斯分布的白噪聲。信源入射角為-5°和5°，搜索角范圍為（-18°，18°），以0.01°為一個(gè)步進(jìn)。采用AP 方法估計(jì)來(lái)波信號(hào)方向角的仿真結(jié)果如圖2 所示。

圖2 列出4 次迭代效果圖，由圖2 可知第一次估計(jì)目標(biāo)1 的初始角度為-0.93°，然后根據(jù)-0.93°估計(jì)出目標(biāo)2 的初始角度為1.75°，第一次迭代根據(jù)目標(biāo)2 的1.75°估計(jì)出目標(biāo)1 角度為-4.56°，再根據(jù)-4.56°估計(jì)出目標(biāo)2 為5.08°，直到估計(jì) 目標(biāo)1 為-4.98°，目標(biāo)2 為5.01°滿足收斂條件結(jié)束迭代，最終測(cè)出角度為-4.98°和5.01°。

圖2 部分AP_DML 算法仿真結(jié)果1

其他條件與上述仿真一致，當(dāng)信源入射角為-1.5°和8°時(shí)，仿真結(jié)果如圖3 所示。

從圖3 可以得出估計(jì)第二種仿真模型的計(jì)算的測(cè)角結(jié)果目標(biāo)1 為-1.52°，目標(biāo)2 為7.91°，滿足收斂條件結(jié)束迭代。如果需要更精確的角度值，可將收斂條件提高，那么迭代次數(shù)也會(huì)增加，計(jì)算角度會(huì)更加精確。

圖3 部分AP_DML 算法仿真結(jié)果2

從上文中以第一種仿真結(jié)果為例，能輕易估算出搜索一次角度為3 600 次，本次迭代次數(shù)為6 次，那么總共需要14 次角度搜索，總搜索次數(shù)為50 400，如果迭代次數(shù)增加，搜索總次數(shù)也隨之增加，耗時(shí)量十分巨大，因此在AP 方法的基礎(chǔ)上采用斐波那契數(shù)列法極值點(diǎn)。斐波那契數(shù)列法求極值點(diǎn)的前提是整個(gè)搜索區(qū)間內(nèi)只有一個(gè)極值點(diǎn)，即一個(gè)波峰或者波谷。以估計(jì)第一個(gè)目標(biāo)的角度為例（搜索其他估計(jì)角度的極值點(diǎn)類似），根據(jù)上述仿真發(fā)現(xiàn)，當(dāng)搜索角度區(qū)間足夠大時(shí)，會(huì)出現(xiàn)多個(gè)波峰情況，因此針對(duì)這種情況本文限定了搜索區(qū)間，除此之外可以看出除初始估計(jì)第一個(gè)信源角度外，其他估計(jì)時(shí)會(huì)出現(xiàn)一個(gè)凹口，因此也不能完全滿足斐波那契數(shù)列法求極值點(diǎn)要求，但是根據(jù)算法得出凹口位置為本次迭代角度值，即初始估計(jì)第二個(gè)信源角度時(shí)凹口為初始估計(jì)第一個(gè)信源角度的值2.25°，因此在用斐波那契數(shù)列法求極值時(shí)可跳過(guò)本次迭代的角度值。具體搜索過(guò)程如圖4 所示（以第一個(gè)仿真模型初始估計(jì)第一個(gè)信源角度為例）。

圖4 中黑色點(diǎn)為整個(gè)斐波那契數(shù)列法的參考點(diǎn)，統(tǒng)計(jì)搜索出峰值點(diǎn)位置的總計(jì)算周期次數(shù)為16 次，與原先3 600 次相比大大減少周期計(jì)算次數(shù)，且最終結(jié)果與之前一致，顯著提高了整個(gè)角度估計(jì)的效率。根據(jù)仿真驗(yàn)證得出當(dāng)滿足斐波那契數(shù)列法求極值的前提條件時(shí)，該算法用于其他角度搜索的最大次數(shù)一般不超過(guò)20次，所以能夠顯著提高峰值搜索效率。

圖4 基于斐波那契數(shù)列法搜索峰值點(diǎn)效果圖

本文利用最大似然估計(jì)算法測(cè)DOA 的前提是信號(hào)源個(gè)數(shù)已知。

4 結(jié)論

本文主要討論了簡(jiǎn)化后的AP 算法用于DOA 估計(jì)。該算法適用于信號(hào)源個(gè)數(shù)已知且與實(shí)際一致，回波信號(hào)相干，低信噪比采樣點(diǎn)比較少的場(chǎng)景。給出了相關(guān)算法的公式及其詳細(xì)的推導(dǎo)過(guò)程，并在均勻線陣上應(yīng)用本文算法，仿真驗(yàn)證了算法的有效性和高效性，且證明了算法對(duì)相干信號(hào)DOA 估計(jì)的能力。