本發明屬于陣列天線信號處理技術領域，特別是涉及一種基于粒子濾波的相干分布源波達方向跟蹤方法。

背景技術：

波達方向估計是陣列天線信號處理領域的重要研究內容。針對波達估計方法，相關研究人員已經提出了capon波束形成、最大似然估計、music譜以及esprit等方法。但是這些方法主要針對點目標源，即信號由一個遠場質點目標發出。這種點目標源模型在雷達、聲吶以及無線通信等領域并不適用，這些場景中的目標源通常在角度空間上存在一定的擴展，因此利用空間分布源模型作為這類場景中目標的波達方向估計模型更為適合。根據分布源的散射特性，可以將其分為相干分布源和非相干分布源兩類。

針對相干分布源的波達方向估計方法，研究人員提出了廣義capon波束形成方法、dpse方法、廣義esprit方法等。但是由于相干分布源參數包括中心波達方向及角分布參數等至少兩個參數，所以上述方向均需要進行多維搜索處理，這使得其計算量大大增加。

上述針對點目標源及分布源的波達估計方法只適用于目標位置保持不變的情況，對于目標運動的非平穩信號的波達方向估計則需要對于每一個快拍的陣列天線接收信號進行一次波達方向估計，并對目標源的波達方向進行跟蹤。目前這類方法主要針對點目標源，包括投影近似子空間跟蹤算法(past)、正交投影近似子空間跟蹤算法(opast)等。

對于分布源的波達方向跟蹤方法僅有fapi-tls-esprit方法，該方法利用fapi算法實現噪聲子空間在每個快拍的更新，再利用tls-esprit算法實現每個快拍的分布源的波達方向估計。該方法僅能實現中心波達方向的跟蹤，不能進行角分布參數的跟蹤。

粒子濾波算法是一種貝葉斯遞推估計方法，能夠實現較強噪聲條件下的不確定參數準確估計。但到目前為止尚未發現利用粒子濾波算法實現相干分布源波達方向跟蹤方面的報道。

技術實現要素：

為了解決上述問題，本發明的目的在于提供一種基于粒子濾波的相干分布源波達方向跟蹤方法。

為了達到上述目的，本發明提供的基于粒子濾波的相干分布源波達方向跟蹤方法包括按順序進行的下列步驟：

2)建立相干分布源的廣義陣列流型及陣列天線接收信號模型的s1階段；

2)基于上述相干分布源的廣義陣列流型及陣列天線接收信號模型，建立基于相干分布源參數譜的觀測似然函數的s2階段；

3)構造相干分布源參數狀態方程的s3階段；

4)基于步驟2)獲得的觀測似然函數和步驟3)獲得的相干分布源參數的狀態方程，利用粒子濾波算法實現相干分布源波達方向跟蹤的s4階段。

在步驟1)中，所述的建立相干分布源的廣義陣列流型及陣列天線接收信號模型的方法是：在等距線陣情況下，首先將分布源的陣列天線接收信號模型描述為角信號密度函數在分布空間積分的形式；然后對于相干分布源模型，則角信號密度函數可以表示為隨機信號幅度與分布源的空間分布函數乘積的形式；假設相干分布源的空間分布函數為高斯分布，則可以建立起相干分布源的廣義陣列流型及陣列天線接收信號模型。

在步驟2)中，所述的基于上述相干分布源的廣義陣列流型及陣列天線接收信號模型，建立基于相干分布源參數譜的觀測似然函數的方法是：在步驟1)建立的相干分布源的廣義陣列流型及陣列天線接收信號模型的基礎上，利用dspe算法估計出接收天線信號的協方差矩陣，并對接收天線信號的協方差矩陣進行特征值分解，選取小特征值對應的特征向量構建噪聲子空間，利用陣列天線導向矢量與噪聲子空間估計相干分布源參數譜，并以該相干分布源參數譜作為觀測似然函數。

在步驟3)中，所述的構造相干分布源參數狀態方程的方法是：以相干分布源的波達方向、波達方向的速度、角分布參數、角分布速度構成相干分布源的狀態向量，然后利用勻速運動模型和上述相干分布源的狀態向量建立相干分布源參數的狀態方程。

在步驟4)中，所述的基于步驟2)獲得的觀測似然函數和步驟3)獲得的相干分布源參數的狀態方程，利用粒子濾波算法實現相干分布源波達方向跟蹤的方法是：使用粒子濾波算法，利用波達方向變化的時間相關性，結合每個快拍的陣列天線接收信號，實現波達方向跟蹤。

本發明提供的基于粒子濾波的相干分布源波達方向跟蹤方法首先建立相干分布源的陣列天線接收信號模型，然后在此模型基礎上利用dspe算法建立相干分布源參數譜的觀測似然函數，并構造相干分布源參數狀態方程，進而使用粒子濾波算法實現相干分布源的波達方向跟蹤。本發明方法具有不需要搜索處理，同時估計中心波達方向及角分布參數，估計精度高、估計性能好等優點。

附圖說明

圖1為本發明提供的基于粒子濾波的相干分布源波達方向跟蹤方法流程圖。

圖2為snr＝3db時本發明方法獲得的相干分布源中心波達方向估計值。

圖3為snr＝3db時本發明方法獲得的相干分布源中心波達方向估計值的rmse。

圖4snr＝3db時本發明方法獲得的相干分布源角分布參數估計值。

圖5為相干分布源中心波達方向估計值隨信噪比變化的平均rmse。

具體實施方式

下面結合附圖和具體實施例對本發明提供的基于粒子濾波的相干分布源波達方向跟蹤方法進行詳細說明。

圖1為本發明提供的基于粒子濾波的相干分布源波達方向跟蹤方法流程圖。其中的全部操作都是在計算機系統中完成的，操作的主體均為計算機系統。

如圖1所示，本發明提供的基于粒子濾波的相干分布源波達方向跟蹤方法包括按順序進行的下列步驟：

3)建立相干分布源的廣義陣列流型及陣列天線接收信號模型的s1階段：

本階段是在等距線陣情況下，利用分布源陣列天線接收信號模型，結合高斯角分布函數，建立相干分布源的廣義陣列流型及陣列天線接收信號模型，然后進入下一步s2階段。

在此階段中，首先將如式(1)所示的分布源的陣列天線接收信號模型描述為角信號密度函數在分布空間積分的形式，即：

其中，y(t)為分布源的陣列天線接收信號，且y(t)＝[y1(t),y2(t),…,yl(t)]^t；a(θ)為陣列天線導向矢量；si(θ-θi,t)為第i個分布源的角信號密度函數，nt為觀測噪聲。

對于相干分布源模型，則角信號密度函數可以表示為隨機信號幅度與分布源的空間分布函數乘積的形式，即：

si(θ-θi,t)＝si(t)gi(θ-θi)(2)

其中，si(t)為隨機信號幅度；gi(θ-θi)為分布源的空間分布函數。

將公式(2)帶入公式(1)，并假設相干分布源的空間分布函數為高斯分布，即：

其中，δi為未知角分布參數，則有：

其中，

將公式(3)帶入公式(5)則可以建立相干分布源的廣義陣列流型及陣列天線接收信號模型，即：

yt＝btst+nt(6)

其中，bt＝[b1(θ1),…,bq(θq)]，且

2)基于上述相干分布源的廣義陣列流型及陣列天線接收信號模型，建立基于相干分布源參數譜的觀測似然函數的s2階段：

本階段是利用dspe(distributedsourceparameterestimation，分布源參數估計)算法，基于上述相干分布源的廣義陣列流型及陣列天線接收信號模型，先估計陣列天線接收信號的協方差矩陣，并對接收信號的協方差矩陣進行特征值分解，選取小特征值對應的特征向量構建噪聲子空間，利用陣列天線導向矢量與噪聲子空間估計相干分布源參數譜，并以此作為觀測似然函數，然后進入下一步s3階段。

在此階段中，陣列天線接收信號的協方差矩陣可以表示為：

對上述協方差矩陣rt進行特征值分解，選取l-q個較小特征值對應的特征向量構成噪聲子空間，l為陣元數，q為分布源個數，即：

un＝[u1,…,ul-q](9)

則dspe相干分布源參數譜可以表示為：

以該此相干分布源參數譜作為觀測似然函數。

3)構造相干分布源參數狀態方程的s3階段：

假設相干分布源波達方向估計服從馬爾可夫過程，以相干分布源的波達方向θk、波達方向的速度角分布參數δk、角分布速度構成相干分布源的狀態向量，即：

利用勻速運動模型和上述相干分布源的狀態向量建立相干分布源參數的狀態方程，即：

xk+1＝axk+vk(12)

其中，a表示狀態轉移矩陣，vk表示隨機擾動噪聲，其滿足零均值高斯分布。

4)基于步驟2)獲得的觀測似然函數和步驟3)獲得的相干分布源參數的狀態方程，利用粒子濾波算法實現相干分布源波達方向跟蹤的s4階段：

假設在k-1快拍時，相干分布源參數的狀態向量可以由一組粒子表示，其權值為

則k快拍時的相干分布源參數的狀態向量可以預測為：

其權值更新為：

其中，為觀測似然函數，本發明以公式(10)所獲得的dspe相干分布源參數譜作為觀測似然函數，即：

并對更新后的權值進行歸一化處理，即：

對更新后的粒子進行重采樣處理，即可獲得相關分布源參數的k快拍時估計值，即：

本發明提供的基于粒子濾波的相干分布源波達方向跟蹤方法的效果可以通過以下仿真結果進一步說明。

仿真數據描述：陣列天線為32個陣元組成的均勻線陣，陣元間距為半波長。仿真相干分布源目標波達方向由-1.5°開始運動，并在100個快拍后運動至-1.0°，相干分布源角分布參數為δ＝2.0°。陣列接收信號的信噪比snr＝3db。

圖2為snr＝3db時本發明方法獲得的相干分布源中心波達方向估計值。其中‘—’為中心波達方向真值，‘---’為本發明方法獲得的中心波達方向估計值。從圖中可以看出本發明方法獲得的中心波達方向估計值更加接近真值，而fapi-tls-esprit方法的估計值則有較大波動，特別是在收斂前(前20個快拍)波動更為明顯。這是由于傳統的分布源中心波達方向估計方法(如fapi-tls-esprit方法)僅能利用當前時刻的陣列天線接收信號進行估計，而忽略了運動的分布源中心波達方向的時間相關性。本發明方法利用了分布源中心波達方向的時間相關性，提高了估計精度。

圖3為snr＝3db時本發明方法獲得的相干分布源中心波達方向估計值的rmse。從圖中可以看出，本發明方法在大約10個快拍后便收斂到rmse小于1°，而fapi-tls-esprit方法的收斂速度明顯較本發明方法更慢。

圖4為snr＝3db時本發明方法獲得的相干分布源角分布參數估計值。其中‘—’為角分布參數真值，‘---’為本發明方法獲得的角分布參數估計值。fapi-tls-esprit方法僅能對分布源對中心波達方向進行跟蹤，分布源的角分布參數則需要在獲得中心波達方向估計值后，利用搜索方法獲得。而本發明方法在對分布源對中心波達方向進行跟蹤的同時，可以對分布源的角分布參數進行估計，更加適合實時性要求較高的分布源跟蹤場景。

圖5為相干分布源中心波達方向估計值隨信噪比變化的平均rmse(200次蒙特卡洛實驗)。由圖中可以看出，本發明方法在信噪比較低的情況下(-10db、-5db)rmse仍然小于10°，fapi-tls-esprit方法的rmse則明顯大于本發明方法。