陳佳

(中國船舶集團有限公司第七二三研究所，江蘇 揚州 225000)

0 引言

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

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

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

式中，sk(t)為第k 個目標(biāo)回波的信號；τi(θk)為第k 個回波信號到達第i 個陣元時相對于參考陣元的時延。有：

上述各個矢量表示分別為：

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

寫成矩陣形式為：

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

2 最大似然估計算法[10-15]

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

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

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

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

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

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

定義A 的投影陣為：

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

可得最大似然準則以跡的形式表示為：

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

AP 算法具體流程如下：

（1）估計第一個信源的情形，為：

（2）估計第二個信源，假定第一個信源的方位為θ1

（3）進行迭代得：

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

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

具體流程如圖1 所示。

圖1 兩維交替投影過程

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

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

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

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

具體步驟如下：

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

計算試探點p1和q1。令k=1。

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

3 仿真結(jié)果

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

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

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

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

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

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

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

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

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

本文利用最大似然估計算法測DOA 的前提是信號源個數(shù)已知。

4 結(jié)論

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