本發(fā)明屬于陣列信號處理領域。

背景技術：

在雷達、通信、聲吶、氣象等領域，陣列信號處理有廣泛而重要的應用。波達方向角(directionofarrival,doa)估計是陣列信號處理中的一種側(cè)向技術，而實際情況中，目標多處于運動狀態(tài)，doa跟蹤是對運動目標doa進行實時估計，高精度的doa跟蹤在個人定位業(yè)務、戰(zhàn)場移動通信和語音信號處理等領域中都有廣泛應用。然而，目標的運動通常導致空間譜擴展，而且運動目標的快拍信號難以長時間積累，這導致傳統(tǒng)靜止目標的doa估計方法不再適用。因此，高精度的doa跟蹤成為了研究重點。

目前，最具有代表性的doa跟蹤方法為：子空間類方法和貝葉斯類方法。

(一)子空間類方法。子空間類方法是分段估計出信號波達方向作為瞬時值，通過跟蹤或者更新子空間實現(xiàn)doa跟蹤。在yangb.anextensionofthepastdalgorithmtobothrankandsubspacetracking[j].signalprocessingletters,ieee,1995,2(9):179-182(一種擴展的實現(xiàn)秩和子空間跟蹤的pastd算法).中提出的緊縮投影近似子空間跟蹤(projectionapproximationsubspacetrackingdeflation,pastd)算法由于計算復雜度較低而被引入到doa跟蹤問題當中。在sanchez-araujoj,marcoss.anefficientpastd-algorithmimplementationformultipledirectionofarrivaltracking[j].ieeetransactionsonsignalprocessing,1999,47(8):2321-2324(一種高效的用于多目標到達角度跟蹤的pastd算法).中用pastd算法實時估計信號子空間和噪聲子空間，并結(jié)合卡爾曼濾波實現(xiàn)doa跟蹤。在低信噪比的情況下，這類方法的跟蹤精度較差。

(二)貝葉斯類方法?；谪惾~斯跟蹤方法是以貝葉斯原理為基礎，利用未知參數(shù)的先驗信息和樣本信息得出后驗信息，然后根據(jù)后驗信息進行參數(shù)估計，可以在相干信源和低信噪比的情況下跟蹤。典型的方法有序列蒙特卡洛方法，即粒子濾波(particlefilter,pf)方法。在ortonm,fitzgeraldw.abayesianapproachtotrackingmultipletargetsusingsensorarraysandparticlefilters[j].signalprocessing,ieeetransactionsonsignalprocessing,2002,50(2):216-223(使用傳感陣列和粒子濾波進行多目標跟蹤的貝葉斯方法).中，將粒子濾波算法用于一維陣列的doa跟蹤。但是粒子濾波算法需要用到狀態(tài)轉(zhuǎn)移概率和觀測噪聲概率函數(shù)，在這些條件未知的情況下，可能出現(xiàn)較大誤差。

技術實現(xiàn)要素：

本發(fā)明是為了解決低信噪比情況下，基于空間類doa跟蹤方法精度差，以及粒子濾波doa跟蹤方法需要的已知條件過多的問題，本發(fā)明提供了一種基于貝葉斯方法的目標方位跟蹤方法。

一種基于貝葉斯方法的目標方位跟蹤方法，它包括如下步驟：

步驟一：初始化當前t時刻的信號到達角θk(t)、信號功率pk(t)、噪聲功率σ²(t)和角度變化量△θk(t)，使得θk(t)＝θk(t-1)、pk(t)＝pk(t-1)、σ²(t)＝σ²(t-1)、△θk(t)＝0；

其中，θk(t-1)、pk(t-1)和σ²(t-1)均已知，k表示第k個信號源，且k＝1,2,3,...,k，k為空間信號源個數(shù)，θk(t-1)為t-1時刻的信號到達角，pk(t-1)為t-1時刻的信號功率，σ²(t-1)為t-1時刻的噪聲功率，角度變化量△θk(t)為當前時刻相對于上一時刻的角度變化量；

其中，t＝1,2,3,...,t′-1,t′，t′為總觀測時間；

步驟二：在θk(t-1)處，進行一階泰勒展開，從而獲得陣列流形矩陣a；

步驟三：采用貝葉斯方法對步驟一中的信號功率pk(t)、噪聲功率σ²(t)和陣列流形矩陣a進行處理，獲得后驗均值μ和后驗方差∑；

步驟四：根據(jù)后驗均值μ和后驗方差∑更新信號功率pk(t)、噪聲功率σ²(t)和角度變化量△θk(t)；

步驟五：采用em算法追蹤m次更新后的信號功率pk(t)、噪聲功率σ²(t)和角度變化量△θk(t)，判斷m次更新后的信號功率pk(t)、噪聲功率σ²(t)和角度變化量△θk(t)的數(shù)值是否收斂，結(jié)果為是，執(zhí)行步驟六，結(jié)果為否，令m＝m+1，執(zhí)行步驟一；

其中，m為重復執(zhí)行步驟一至步驟四的次數(shù)；

步驟六：收斂狀態(tài)下，第m次更新后的信號功率pk(t)為當前t時刻信號功率pk(t)的估計值，更新后的噪聲功率σ²(t)為當前t時刻噪聲功率σ²(t)的估計值，更新后的角度變化量△θk(t)為當前t時刻角度變化量θk(t-1)的估計值，

根據(jù)當前t時刻角度變化量△θk(t)的估計值，獲得當前t時刻的信號到達角θk(t)的估計值，其中，θk(t)＝θk(t-1)+△θk(t)；

步驟七：判斷t是否等于t′，結(jié)果為否，令t＝t+1，執(zhí)行步驟一；結(jié)果為是，執(zhí)行步驟八；

步驟八：完成了總觀測時間t′內(nèi)，各時刻的信號功率pk(t)、噪聲功率σ²(t)和信號到達角θk(t)的估計，從而實現(xiàn)對目標方位的實時追蹤。

所述的步驟二中，在θk(t-1)處，進行一階泰勒展開，從而獲得陣列流形矩陣a的具體過程為：

步驟二一：根據(jù)泰勒展開式，將t時刻的流型向量a(θk(t))，在θk(t-1)處進行一階泰勒展開得到：

a(θk(t))＝a(θk(t-1))+a'(θk(t-1))△θk(t)(公式一)，

其中，

a(θk(t-1))為t-1時刻的流型向量，a′(θk(t-1))表示t-1時刻的流型向量的導數(shù)，e為自然常數(shù)，j為虛數(shù)單位，d為陣元間距，λ為信號波長，m為陣元個數(shù)，[]^t表示矩陣轉(zhuǎn)置；

步驟二二：根據(jù)流型向量a(θk(t))，獲得t時刻的陣列流形矩陣a；

a＝a1+bδ(公式二)，

其中，

a1＝[a(θ1(t-1)),a(θ2(t-1)),a(θ3(t-1)),...,a(θk-1(t-1)),a(θk(t-1))]，

b＝[a'(θ1(t-1)),a'(θ2(t-1)),a'(θ3(t-1)),...,a'(θk-1(t-1)),a'(θk(t-1))]，

δ＝diag([△θ1(t),△θ2(t),△θ3(t),...,△θk-1(t),△θk(t)])，

a1為t-1時刻的陣列流形矩陣，b為t-1時刻流形向量的導數(shù)構成的矩陣，δ為角度校正量，diag()表示對角矩陣函數(shù)。

所述的步驟三中，采用貝葉斯方法對步驟一中的信號功率pk(t)、噪聲功率σ²(t)和陣列流形矩陣a進行處理，獲得后驗均值μ和后驗方差∑的具體過程為：

μ＝pa^h(σ²(t)i+apa^h)^-1x(公式三)，

∑＝(p^-1+σ^-2(t)a^ha)^-1(公式四)，

其中，

p＝diag([p1(t),p2(t),p3(t),...,pk-1(t),pk(t)])，

p為信號的協(xié)方差矩陣，diag()表示對角矩陣函數(shù)，i為單位矩陣。

所述的步驟四中，根據(jù)后驗均值μ和后驗方差∑更新信號功率pk(t)、噪聲功率σ²(t)和角度變化量△θk(t)的具體過程為：

[△θ1(t),△θ2(t),△θ3(t),...,△θk-1(t),△θk(t)]^t＝u^-1v(公式七)，

其中，

(∑)(n,n)表示∑的第n行第n列元素，

(μ)n表示μ的第n行，n為正整數(shù)；

x＝as+n，

其中，上標h表示共軛轉(zhuǎn)置，表示取實部，⊙表示哈達瑪乘積，l表示快拍數(shù)，diag()表示表示對角矩陣函數(shù)，v和u均表示中間變量，x∈c^m×l為陣列接收數(shù)據(jù)，s∈c^k×l為入射信號，n∈c^m×l為觀測噪聲，c為復數(shù)集合，tr()表示矩陣的跡。

本發(fā)明考慮的是角度慢變模式下的doa跟蹤，即相鄰時刻信號到達角度值及信號功率的變化很小。本發(fā)明在信號到達角和功率緩慢變化的情形下，根據(jù)前一時刻的信號到達角、信號功率和噪聲功率的估計值，以及當前時刻的觀測值，利用em算法，對前一時刻的信息進行校正，進而實現(xiàn)對信號到達角度的追蹤。

本發(fā)明帶來的有益效果是，

在本發(fā)明中，利用泰勒展開公式，將角度慢變模式下的doa跟蹤建模為一個動態(tài)模型，并基于貝葉斯理論，將角度的跟蹤轉(zhuǎn)化為概率模型中的參數(shù)估計問題，采用em算法來求解。本發(fā)明提出的跟蹤方法的優(yōu)勢在于能夠同時實時跟蹤信號的到達角度、信號功率以及噪聲功率，且無需知道狀態(tài)轉(zhuǎn)移概率等先驗信息，而且與稀疏貝葉斯學習方法相比，計算量小。

本發(fā)明所述一種基于貝葉斯方法的目標方位跟蹤方法與子空間類doa跟蹤方法相比，本發(fā)明跟蹤方法精度提高30％以上，且在低信噪比的情況下尤為明顯。

本發(fā)明所述一種基于貝葉斯方法的目標方位跟蹤方法與粒子濾波doa跟蹤方法相比，本發(fā)明所述方法無需知道狀態(tài)轉(zhuǎn)移概率和觀測噪聲概率函數(shù)。

本發(fā)明在進行doa跟蹤的同時，本發(fā)明還可以實現(xiàn)信號功率和噪聲功率的多參數(shù)聯(lián)合跟蹤。

附圖說明

圖1為本發(fā)明所述的一種基于貝葉斯方法的目標方位跟蹤方法的流程圖；

圖2為有兩個入射信號情形下的由本發(fā)明方法得到的信號到達角(doa)追蹤圖；

圖3為利用本發(fā)明所述的一種基于貝葉斯方法的目標方位跟蹤方法，得到的信號功率追蹤圖；

圖4為利用本發(fā)明所述的一種基于貝葉斯方法的目標方位跟蹤方法，得到的噪聲功率追蹤圖；

圖5為本發(fā)明方法與pastd-kalman(緊縮投影近似子空間跟蹤-卡爾曼濾波)算法進行doa跟蹤的均方根誤差與信噪比關系曲線的對比圖。

具體實施方式

具體實施方式一：參見圖1說明本實施方式，本實施方式所述的一種基于貝葉斯方法的目標方位跟蹤方法，它包括如下步驟：

步驟一：初始化當前t時刻的信號到達角θk(t)、信號功率pk(t)、噪聲功率σ²(t)和角度變化量△θk(t)，使得θk(t)＝θk(t-1)、pk(t)＝pk(t-1)、σ²(t)＝σ²(t-1)、△θk(t)＝0；

其中，θk(t-1)、pk(t-1)和σ²(t-1)均已知，k表示第k個信號源，且k＝1,2,3,...,k，k為空間信號源個數(shù)，θk(t-1)為t-1時刻的信號到達角，pk(t-1)為t-1時刻的信號功率，σ²(t-1)為t-1時刻的噪聲功率，角度變化量△θk(t)為當前時刻相對于上一時刻的角度變化量；

其中，t＝1,2,3,...,t′-1,t′，t′為總觀測時間；

步驟二：在θk(t-1)處，進行一階泰勒展開，從而獲得陣列流形矩陣a；

步驟三：采用貝葉斯方法對步驟一中的信號功率pk(t)、噪聲功率σ²(t)和陣列流形矩陣a進行處理，獲得后驗均值μ和后驗方差∑；

步驟四：根據(jù)后驗均值μ和后驗方差∑更新信號功率pk(t)、噪聲功率σ²(t)和角度變化量△θk(t)；

步驟五：采用em算法追蹤m次更新后的信號功率pk(t)、噪聲功率σ²(t)和角度變化量△θk(t)，判斷m次更新后的信號功率pk(t)、噪聲功率σ²(t)和角度變化量△θk(t)的數(shù)值是否收斂，結(jié)果為是，執(zhí)行步驟六，結(jié)果為否，令m＝m+1，執(zhí)行步驟一；

其中，m為重復執(zhí)行步驟一至步驟四的次數(shù)；

步驟六：收斂狀態(tài)下，第m次更新后的信號功率pk(t)為當前t時刻信號功率pk(t)的估計值，更新后的噪聲功率σ²(t)為當前t時刻噪聲功率σ²(t)的估計值，更新后的角度變化量△θk(t)為當前t時刻角度變化量θk(t-1)的估計值，

根據(jù)當前t時刻角度變化量△θk(t)的估計值，獲得當前t時刻的信號到達角θk(t)的估計值，其中，θk(t)＝θk(t-1)+△θk(t)；

步驟七：判斷t是否等于t′，結(jié)果為否，令t＝t+1，執(zhí)行步驟一；結(jié)果為是，執(zhí)行步驟八；

步驟八：完成了總觀測時間t′內(nèi)，各時刻的信號功率pk(t)、噪聲功率σ²(t)和信號到達角θk(t)的估計，從而實現(xiàn)對目標方位的實時追蹤。

本實施方式，在本發(fā)明中，利用泰勒展開公式，將角度慢變模式下的doa跟蹤建模為一個動態(tài)模型，并基于貝葉斯理論，將角度的跟蹤轉(zhuǎn)化為概率模型中的參數(shù)估計問題，采用em算法來求解。本發(fā)明提出的跟蹤方法的優(yōu)勢在于能夠同時實時跟蹤信號的到達角度、信號功率以及噪聲功率，且無需知道狀態(tài)轉(zhuǎn)移概率等先驗信息，而且與稀疏貝葉斯學習方法相比，計算量小。

em算法的英文全稱為：expectationmaximization。

原理分析：首先根據(jù)角度慢變的假設，根據(jù)一階泰勒展開式，得到與上一時刻角度和角度更新量有關的陣列流形矩陣，以便有效利用上一時刻的信息。然后進行em迭代，每次迭代中根據(jù)上一輪求得的信號功率、噪聲功率及角度校正量的估計值計算信號的后驗概率密度，根據(jù)該后驗概率密度，對信號與噪聲的聯(lián)合概率密度求均值并使之最大化從而更新參數(shù)的估計值。獲得角度校正量后，即可計算出實時的到達角度。

具體實施方式二：參見圖1說明本實施方式，本實施方式與具體實施方式一所述的一種基于貝葉斯方法的目標方位跟蹤方法的區(qū)別在于，所述的步驟二中，在θk(t-1)處，進行一階泰勒展開，從而獲得陣列流形矩陣a的具體過程為：

步驟二一：根據(jù)泰勒展開式，將t時刻的流型向量a(θk(t))，在θk(t-1)處進行一階泰勒展開得到：

a(θk(t))＝a(θk(t-1))+a'(θk(t-1))△θk(t)(公式一)，

其中，

a(θk(t-1))為t-1時刻的流型向量，a′(θk(t-1))表示t-1時刻的流型向量的導數(shù)，e為自然常數(shù)，j為虛數(shù)單位，d為陣元間距，λ為信號波長，m為陣元個數(shù)，[]^t表示矩陣轉(zhuǎn)置；

步驟二二：根據(jù)流型向量a(θk(t))，獲得t時刻的陣列流形矩陣a；

a＝a1+bδ(公式二)，

其中，

a1＝[a(θ1(t-1)),a(θ2(t-1)),a(θ3(t-1)),...,a(θk-1(t-1)),a(θk(t-1))]，

b＝[a'(θ1(t-1)),a'(θ2(t-1)),a'(θ3(t-1)),...,a'(θk-1(t-1)),a'(θk(t-1))]，

δ＝diag([△θ1(t),△θ2(t),△θ3(t),...,△θk-1(t),△θk(t)])，

a1為t-1時刻的陣列流形矩陣，b為t-1時刻流形向量的導數(shù)構成的矩陣，δ為角度校正量，diag()表示對角矩陣函數(shù)。

具體實施方式三：參見圖1說明本實施方式，本實施方式與具體實施方式一所述的一種基于貝葉斯方法的目標方位跟蹤方法的區(qū)別在于，所述的步驟三中，采用貝葉斯方法對步驟一中的信號功率pk(t)、噪聲功率σ²(t)和陣列流形矩陣a進行處理，獲得后驗均值μ和后驗方差∑的具體過程為：

μ＝pa^h(σ²(t)i+apa^h)^-1x(公式三)，

∑＝(p^-1+σ^-2(t)a^ha)^-1(公式四)，

其中，

p＝diag([p1(t),p2(t),p3(t),...,pk-1(t),pk(t)])，

p為信號的協(xié)方差矩陣，diag()表示對角矩陣函數(shù)，i為單位矩陣。

具體實施方式四：參見圖1說明本實施方式，本實施方式與具體實施方式一所述的一種基于貝葉斯方法的目標方位跟蹤方法的區(qū)別在于，所述的步驟四中，根據(jù)后驗均值μ和后驗方差∑更新信號功率pk(t)、噪聲功率σ²(t)和角度變化量△θk(t)的具體過程為：

[△θ1(t),△θ2(t),△θ3(t),...,△θk-1(t),△θk(t)]^t＝u^-1v(公式七)，

其中，

(∑)(n,n)表示∑的第n行第n列元素，

(μ)n表示μ的第n行，n為正整數(shù)；

x＝as+n，

其中，上標h表示共軛轉(zhuǎn)置，表示取實部，⊙表示哈達瑪乘積，l表示快拍數(shù)，diag()表示表示對角矩陣函數(shù)，v和u均表示中間變量，x∈c^m×l為陣列接收數(shù)據(jù)，s∈c^k×l為入射信號，n∈c^m×l為觀測噪聲，c為復數(shù)集合，tr()表示矩陣的跡。

具體實施方式五：參見圖1說明本實施方式，本實施方式與具體實施方式二所述的一種基于貝葉斯方法的目標方位跟蹤方法的區(qū)別在于，所述的步驟二二中，a＝[a(θ1(t)),a(θ2(t)),a(θ3(t)),...,a(θk-1(t)),a(θk(t))]。

仿真驗證：

利用本發(fā)明所述的一種基于貝葉斯方法的目標方位跟蹤方法進行仿真試驗：

步驟一一，仿真參數(shù)設置如下：陣元個數(shù)m＝12，快拍數(shù)l＝20，陣元間距d＝λ/2，噪聲功率σ²(t)＝1，觀測時間t′＝100；

步驟一二：建立仿真信號運動狀態(tài)模型：

xk(t)＝fxk(t-1)+gvk(t)

式中，vk(t)為零均值的高斯白噪聲，且方差為0.001，其中表示θk(t)的速度；

步驟一三：設置入射信號個數(shù)k＝2，初始角度分別為40°和30°，信噪比的變化分別為8db至-3db和3db至8db，得到的doa、信號功率和噪聲功率跟蹤曲線分別如圖2、圖3和圖4所示，其中，doa為信號到達角；

步驟一四：為了更好地分析本發(fā)明方法的doa跟蹤性能，將其與pastd算法結(jié)合卡爾曼濾波的doa跟蹤方法進行比較。進行100次蒙特卡洛試驗，得到兩種算法的均方根誤差隨信噪比的關系曲線如圖5所示。

仿真結(jié)果證明：

由圖2、圖3和圖4可以看出，本發(fā)明的方法可以很好的進行信號到達角度、信號功率和噪聲功率的實時跟蹤；由圖5可以看出，本發(fā)明的方法與pastd算法結(jié)合卡爾曼濾波的doa跟蹤方法相比具有更高的跟蹤精度，在低信噪比的情況下優(yōu)勢更加明顯。