本發(fā)明涉及陣列信號(hào)處理領(lǐng)域，尤其是一種波達(dá)方向估計(jì)方法。

背景技術(shù)：

波達(dá)方向估計(jì)在雷達(dá)、聲納和無(wú)線通信中，有很廣泛的應(yīng)用，然而傳統(tǒng)的波達(dá)方向估計(jì)方法只能求解目標(biāo)數(shù)少于陣元數(shù)的情況，因此如何用少量的陣元檢測(cè)更多的目標(biāo)是一個(gè)值得研究的問(wèn)題。實(shí)際中，一般常用的N元均勻線列陣，最多只能估計(jì)出N-1個(gè)目標(biāo)的來(lái)波方向。近幾年，提出了一種新型幾何結(jié)構(gòu)的線列陣——互質(zhì)陣，使得估計(jì)的目標(biāo)數(shù)目遠(yuǎn)超過(guò)陣元數(shù)目N，事實(shí)上，一個(gè)有N+M陣元的互質(zhì)陣所能估計(jì)的目標(biāo)數(shù)目可以達(dá)到O(MN)，即與MN同階大小的目標(biāo)數(shù)目，由于互質(zhì)陣陣元位置的巧妙分布，在經(jīng)過(guò)數(shù)學(xué)運(yùn)算處理后，可以形成孔徑更大的虛擬陣。具體來(lái)說(shuō)，上述提到的N+2M-1個(gè)物理陣元的互質(zhì)陣可以虛擬出孔徑為2MN+1的均勻線列陣。

對(duì)于互質(zhì)陣，已經(jīng)提出了一種基于互質(zhì)陣的空域平滑MUSIC算法，對(duì)于使用N+2M-1個(gè)物理陣元的互質(zhì)陣，使用該方法最多可以估計(jì)出MN個(gè)目標(biāo)，盡管這個(gè)算法估計(jì)出的目標(biāo)數(shù)目已經(jīng)遠(yuǎn)大于同等數(shù)目物理陣元的均勻線列陣，但是由于使用了空域平滑方法，因此導(dǎo)致虛擬陣元減少到MN+1，檢測(cè)性能降低。最近幾年，由D.Donoho、E.Candes及華裔科學(xué)家T.Tao等人提出了一種新的信息獲取指導(dǎo)理論，即壓縮感知，該理論一經(jīng)提出，就在信息論、信號(hào)處理、模式識(shí)別、無(wú)線通信等領(lǐng)域受到高度關(guān)注，該理論認(rèn)為，對(duì)于一個(gè)稀疏信號(hào)或者是可以稀疏表示的信號(hào)，可以對(duì)該信號(hào)進(jìn)行少量觀測(cè)后，然后采用特定的重構(gòu)算法重構(gòu)原始信號(hào)，對(duì)于信號(hào)x的有限觀測(cè)y，滿足y＝Ax，式中觀測(cè)矩陣M＜N，若滿足一定條件，便可從有限觀測(cè)y中重構(gòu)原始的稀疏信號(hào)x。

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

為了克服現(xiàn)有技術(shù)的不足，本發(fā)明基于壓縮感知中的稀疏重構(gòu)思想，對(duì)于互質(zhì)陣構(gòu)造基于OMP算法的波達(dá)方向估計(jì)方法，克服了空域平滑MUSIC算法導(dǎo)致虛擬陣元減少的問(wèn)題，使得檢測(cè)目標(biāo)數(shù)目超過(guò)MN，并且仿真證明該方法在高信噪比時(shí)有很好的估計(jì)精度。

本發(fā)明解決其技術(shù)問(wèn)題所采用的技術(shù)方案具體包括如下步驟：

第一步：互質(zhì)陣結(jié)構(gòu)設(shè)計(jì)

互質(zhì)陣由兩個(gè)間距分別為Md和Nd的均勻線列陣嵌套組成，其中M和N為兩個(gè)互為質(zhì)數(shù)的常數(shù)，d表示接收信號(hào)的半波長(zhǎng)，一個(gè)均勻線列陣有N個(gè)陣元，陣元間距為Md，另外一個(gè)均勻線列陣有2M個(gè)陣元，陣元間距為Nd，兩個(gè)均勻線列陣的第一個(gè)陣元重合，共有N+2M-1個(gè)陣元，由于M和N都是質(zhì)數(shù)，所以稱這種線列陣為互質(zhì)陣，互質(zhì)陣是一種特殊結(jié)構(gòu)的非均勻線列陣；

第二步：估計(jì)互質(zhì)陣接收信號(hào)的協(xié)方差矩陣

互質(zhì)陣接收信號(hào)的協(xié)方差矩陣R_xx通過(guò)公式(1)估計(jì)：

式(1)中，表示互質(zhì)陣接收信號(hào)向量，T表示陣列快拍數(shù)；

第三步：構(gòu)造虛擬均勻線列陣的接收向量y

構(gòu)造大小為(N+2M-1)×(N+2M-1)標(biāo)識(shí)矩陣B，B中第i行第j列的元素b_ij＝i-j，標(biāo)識(shí)矩陣B和協(xié)方差矩陣R_xx分別按列拉直，得到標(biāo)識(shí)向量從標(biāo)識(shí)向量b中按照從小到大的順序依次尋找元素值為-MN到MN的位置，并按順序記錄下2MN+1個(gè)位置信息，然后從向量z中依次將對(duì)應(yīng)位置的元素提取出來(lái)，構(gòu)造出虛擬均勻線列陣的接收向量

第四步：待檢測(cè)的角度區(qū)域剖分網(wǎng)格，構(gòu)造陣列流形矩陣A

將待檢測(cè)的角度區(qū)域離散化剖分網(wǎng)格，形成網(wǎng)格向量θ＝[θ₁,θ₂,…,θ_D]^T，其中θ_k表示離散的網(wǎng)格角度，k＝1,…,D，D表示網(wǎng)格數(shù)目，且使得剖分的網(wǎng)格數(shù)目D大于信號(hào)數(shù)目，由公式(2)和公式(3)構(gòu)造陣列流形矩陣A：

A(θ₁,θ₂,…,θ_D)＝[a(θ₁),a(θ₂),…a(θ_k),...,a(θ_D)] (2)

式(2)中，a(θ_k)為對(duì)應(yīng)陣列的陣列流形向量，k＝1，2，...，D，式(3)中，j表示虛數(shù)單位，d表示接收信號(hào)的半波長(zhǎng)，λ表示接收信號(hào)的波長(zhǎng)，θ_k表示離散的網(wǎng)格角度；

第五步：利用OMP算法估計(jì)稀疏信號(hào)

陣列流形矩陣A(θ₁,θ₂,…,θ_D)和虛擬均勻線列陣的接收向量y滿足如下方程：

式(4)中p為對(duì)應(yīng)于剖分網(wǎng)格θ＝[θ₁,θ₂,…,θ_D]^T的稀疏向量，表示噪聲向量，稀疏向量p采用壓縮感知中的OMP算法求解可得稀疏向量估計(jì)值

第六步：由稀疏向量得到信號(hào)波達(dá)方向

若稀疏向量估計(jì)值中的第i項(xiàng)非0，則表示對(duì)應(yīng)的θ_i方向具有信號(hào)，否則表示沒(méi)有信號(hào)。

本發(fā)明相比傳統(tǒng)的波達(dá)方向估計(jì)采用均勻線列陣，由于采用了互質(zhì)陣這種新型的嵌套陣，使得波達(dá)方向估計(jì)性能顯著的提高了估計(jì)的目標(biāo)數(shù)量，在陣列多目標(biāo)估計(jì)時(shí)具有很好的估計(jì)性能；在同樣使用互質(zhì)陣時(shí)，該方法相比于空間平滑的MUSIC算法而言，由于使用了壓縮感知的OMP算法，使得估計(jì)的目標(biāo)數(shù)目超過(guò)NM個(gè)目標(biāo)，在同樣陣元數(shù)量下，進(jìn)一步提高了互質(zhì)陣的估計(jì)目標(biāo)數(shù)目；由于使用了OMP算法，可以簡(jiǎn)單高效的對(duì)稀疏信號(hào)進(jìn)行重構(gòu)，最終快速高效的給出波達(dá)方向估計(jì)結(jié)果；在高信噪比時(shí)，本發(fā)明方法的估計(jì)精度優(yōu)于空間平滑的MUSIC算法，總體而言，本發(fā)明的性能不低于空間平滑的MUSIC算法。

附圖說(shuō)明

圖1是本發(fā)明進(jìn)行波達(dá)方向估計(jì)的方法流程圖。

圖2是本發(fā)明互質(zhì)陣的幾何結(jié)構(gòu)。

圖3是本發(fā)明與空間平滑MUSIC算法的波達(dá)方向估計(jì)結(jié)果。

圖4是本發(fā)明對(duì)16個(gè)來(lái)波信號(hào)的波達(dá)方向估計(jì)結(jié)果。

圖5是本發(fā)明與空間平滑MUSIC算法的單目標(biāo)的估計(jì)精度對(duì)比結(jié)果。

圖5中MUSIC表示空間平滑MUSIC算法，OMP表示本發(fā)明提出的基于OMP算法的波達(dá)方向估計(jì)方法。

具體實(shí)施方式

下面結(jié)合附圖和實(shí)施例對(duì)本發(fā)明進(jìn)一步說(shuō)明。

第一步：互質(zhì)陣結(jié)構(gòu)設(shè)計(jì)

選取兩個(gè)互為質(zhì)數(shù)的常數(shù)M和N，用常數(shù)d表示接收信號(hào)的半波長(zhǎng)，互質(zhì)陣由兩個(gè)間距分別為Md和Nd的均勻線列陣嵌套組成，其中一個(gè)均勻線陣有N個(gè)陣元，陣元間距為Md，另外一個(gè)線陣有2M個(gè)陣元，陣元間距為Nd，兩個(gè)線列陣第一個(gè)陣元重合，總共有N+2M-1個(gè)陣元，由于M和N都是質(zhì)數(shù)，所以稱這種線列陣為互質(zhì)陣，互質(zhì)陣的結(jié)構(gòu)如圖1所示，互質(zhì)陣是一種特殊結(jié)構(gòu)的非均勻線列陣，因此可以在線列陣的模型上進(jìn)一步推導(dǎo)；

首先考慮包含N個(gè)陣元的線陣，陣元位置為l_i，假設(shè)K個(gè)目標(biāo)方位為并且設(shè)K個(gè)目標(biāo)信號(hào)為窄帶信號(hào)，那么陣列接收信號(hào)表示為：

式中1≤t≤T，s_k(t)是第k個(gè)來(lái)波信號(hào)，n(t)是獨(dú)立同分布的白噪聲，是陣列流形向量，且

的第i個(gè)元素代表第k個(gè)信號(hào)在第i個(gè)陣元的時(shí)間延遲所帶來(lái)的相位變化。

第二步：估計(jì)互質(zhì)陣接收信號(hào)的協(xié)方差矩陣

互質(zhì)陣接收信號(hào)的協(xié)方差矩陣R_xx通過(guò)公式(1)估計(jì)：

式中，表示互質(zhì)陣接收信號(hào)向量，T表示陣列快拍數(shù)；

設(shè)信號(hào)s_k(t)服從方差為的獨(dú)立高斯分布，考慮每個(gè)陣元接收數(shù)據(jù)的二階統(tǒng)計(jì)量，求解陣列接收信號(hào)x(t)的協(xié)方差矩陣R_xx

式中σ²是噪聲功率，將矩陣R_xx按列拉直，根據(jù)矩陣Kronecker積(也稱直積)的運(yùn)算性質(zhì)，由式(7)得到向量z：

其中，表示第k個(gè)發(fā)射信號(hào)的方差，表達(dá)式如下：

代表除了第i個(gè)位置為1其余元素都為0的列向量，根據(jù)(8)式將向量z看作是信源向量q經(jīng)過(guò)矩陣Φ觀測(cè)后接收到的信號(hào)，矩陣中存在一個(gè)陣列流形矩陣，并且這個(gè)矩陣是一個(gè)具有更多虛擬陣元的線列陣的陣列流形矩陣。

第三步：構(gòu)造虛擬均勻線列陣的接收向量y

構(gòu)造大小為(N+2M-1)×(N+2M-1)標(biāo)識(shí)矩陣B，B中第i行第j列的元素b_ij＝i-j，標(biāo)識(shí)矩陣B和協(xié)方差矩陣R_xx分別按列拉直，得到標(biāo)識(shí)向量從標(biāo)識(shí)向量b中從小到大的順序依次尋找元素值為-MN到MN的位置，按順序記錄下共2MN+1個(gè)位置信息，假設(shè)b＝[-3,-1,1,5,0]，從b中以小到大的順序依次尋找元素值為-1到1的位置，記錄下的對(duì)應(yīng)位置信息就是[2,5,3]，然后從向量z中依次將對(duì)應(yīng)位置的元素提取出來(lái)，構(gòu)造出虛擬均勻線列陣的接收向量

第四步：待檢測(cè)的角度區(qū)域剖分網(wǎng)格，構(gòu)造陣列流形矩陣A

將待檢測(cè)的角度區(qū)域離散化剖分網(wǎng)格，形成網(wǎng)格向量θ＝[θ₁,θ₂,…,θ_D]^T，其中θ_k(k＝1,…,D)表示離散的網(wǎng)格角度，D表示網(wǎng)格數(shù)目，且使得剖分的網(wǎng)格數(shù)目D大于信號(hào)數(shù)目，取D為信號(hào)數(shù)目的10倍以上，由下面公式構(gòu)造陣列流形矩陣A：

A(θ₁,θ₂,…,θ_D)＝[a(θ₁),a(θ₂),…,a(θ_D)] (2)

式(2)中，a(θ_k)為對(duì)應(yīng)陣列的陣列流形向量，k＝1，2，...，D，式(3)中，j表示虛數(shù)單位，d表示接收信號(hào)的半波長(zhǎng)，λ表示接收信號(hào)的波長(zhǎng)，θ_k表示離散的網(wǎng)格角度；

第五步：利用OMP算法估計(jì)稀疏信號(hào)

陣列流形矩陣A(θ₁,θ₂,…,θ_D)和虛擬均勻線列陣的接收向量y滿足如下方程：

式(4)中p為對(duì)應(yīng)于剖分網(wǎng)格θ＝[θ₁,θ₂,…,θ_D]^T的稀疏向量，表示噪聲向量，稀疏向量p采用壓縮感知中的OMP算法求解可得稀疏向量估計(jì)值

第六步：由稀疏向量得到信號(hào)波達(dá)方向

若稀疏向量估計(jì)值中的第i項(xiàng)非0，則表示θ_i方向具有信號(hào)，否則表示沒(méi)有信號(hào)。

在本發(fā)明中，更具體的假設(shè)這個(gè)線列陣為如圖1所示的互質(zhì)陣，由兩個(gè)間距分別為Md和Nd的均勻線列陣嵌套組成，一共有N+2M-1個(gè)物理陣元，其中一個(gè)均勻線陣有N個(gè)陣元，另外一個(gè)線陣有2M個(gè)陣元，兩個(gè)線陣的第一個(gè)陣元重合。假設(shè)信號(hào)的協(xié)方差矩陣R_xx已知，矩陣大小為(N+2M-1)×(N+2M-1)，根據(jù)公式(4)，由直積的性質(zhì)可以知道，矩陣Φ共有(N+2M-1)²個(gè)行向量，已經(jīng)證明了，如果N,M都是質(zhì)數(shù)，這些行向量的子集對(duì)應(yīng)于一個(gè)孔徑更大的均勻線列陣的陣列流形矩陣A，其陣元位置從-MNd到MNd，陣元間距為d，這時(shí)虛擬均勻線列陣的孔徑可以達(dá)到(2MN+1)d，下面由矩陣構(gòu)造均勻線列陣的陣列流形矩陣

從Φ中選擇相應(yīng)的行向量，并對(duì)這些行向量進(jìn)行排序得到A，使得陣列流形矩陣A中的第(n,k)個(gè)元素為k＝1,2,…,K,n＝-MN,…,0,…,MN，最終形成大小為(2NM+1)×K的矩陣A，此時(shí)A便對(duì)于一個(gè)孔徑為(2MN+1)d的均勻線列陣的陣列流形矩陣，對(duì)于式(8)左端的向量z，按照步驟3中所述構(gòu)造向量y的方法選擇相應(yīng)的項(xiàng)并且排序得到向量y，y便是對(duì)應(yīng)于此虛擬均勻線列陣的接收向量，那么A、y滿足

式中向量是向量1_n去掉相應(yīng)項(xiàng)排序得到的向量，并且第NM+1項(xiàng)為1，其余均為0。

下面對(duì)壓縮感知的稀疏重構(gòu)思想進(jìn)行波達(dá)方向估計(jì)進(jìn)行說(shuō)明：

首先對(duì)待檢測(cè)的角度區(qū)域離散化：θ＝[θ₁,θ₂,…,θ_D]^T，且D＞＞K，即剖分的網(wǎng)格數(shù)目大于信號(hào)數(shù)目，構(gòu)造陣列流形矩陣A(θ₁,θ₂,…,θ_D)：

A(θ₁,θ₂,…,θ_D)＝[a(θ₁),a(θ₂),…,a(θ_D)] (2)

式中是一個(gè)陣元數(shù)為2NM+1的均勻線列陣的陣列流形向量。因此可以將公式(4)改寫為

其中，y是根據(jù)互質(zhì)陣的接收信號(hào)x(t)計(jì)算得到的虛擬均勻線列陣的接收向量，如果第k個(gè)信號(hào)的來(lái)波方向?yàn)棣?sub>i，那么向量p的第i個(gè)位置假設(shè)所有目標(biāo)的位置恰好位于剖分的網(wǎng)格點(diǎn)上，那么只在K個(gè)來(lái)波方向上是非零的，此時(shí)，p為K稀疏向量。對(duì)于公式(9)，如果給定y和A，通過(guò)壓縮感知中的稀疏重構(gòu)算法來(lái)求解方程。

根據(jù)壓縮感知理論，需要求得滿足公式(10)的最稀疏解因此稀疏重構(gòu)問(wèn)題表示為：

式中||·||₀表示向量非零項(xiàng)的數(shù)目，為了表示方便A＝A(θ₁,θ₂,…,θ_D)，假設(shè)噪聲方差σ²是未知的，為了求解式(10)，本發(fā)明采用OMP算法求解，OMP算法的迭代步驟如下：

輸入：陣列流形矩陣A，虛擬均勻線列陣的接收向量y，目標(biāo)信號(hào)數(shù)量K；

輸出：p的K稀疏逼近誤差向量E；

初始化：余量E⁰＝y(tǒng)，重構(gòu)信號(hào)p⁰＝0，索引集迭代次數(shù)n＝0；

步驟1：計(jì)算余量和構(gòu)造陣列流形矩陣A(θ₁,θ₂,…,θ_D)的每一列之間的內(nèi)積Gⁿ＝A^TE^n-1，E^n-1是第n次迭代的余量，n＝1時(shí)，E^n-1＝E⁰；

步驟2：找出Gⁿ中絕對(duì)值最大的元素對(duì)應(yīng)的位置，即

步驟3：更新索引集Γⁿ＝Γ^n-1∪{k}及原子集合

步驟4：利用最小二乘法求得近似解

步驟5：更新余量

步驟6：判斷迭代是否滿足停止條件，即迭代次數(shù)n是否大于K，當(dāng)?shù)螖?shù)大于K則滿足停止條件，若滿足停止條件，則E＝Eⁿ，輸出否則令n←n+1，轉(zhuǎn)步驟1，其中，是的第i項(xiàng)。

以下對(duì)本發(fā)明的方法進(jìn)一步描述：在以本發(fā)明技術(shù)方案為前提下進(jìn)行實(shí)施，給出了詳細(xì)的實(shí)施方式和具體的操作過(guò)程。

考慮包含10個(gè)物理陣元的互質(zhì)陣，即對(duì)圖2中的互質(zhì)陣取N＝5，M＝3，具體來(lái)說(shuō)，第一層陣元位置在[0,3,6,9,12]d，第二層陣元位置在[0,5,10,15,20,25]d，取d為半波長(zhǎng)值，選取若干窄帶信號(hào)，來(lái)波方向均勻分布在-60度到60度的區(qū)間內(nèi)，頻率均取為f＝1000Hz，采樣頻率fs＝8192Hz，快拍數(shù)500，信噪比為0dB。

圖3中選取了15個(gè)不同方向的窄帶信號(hào)仿真兩種算法的檢測(cè)性能。圖4中左圖是按照空間平滑矩陣MUSIC算法得到的結(jié)果。圖3中右圖展示了基于OMP算法的波達(dá)方向估計(jì)。首先從-60度到60度以0.5度為步長(zhǎng)劃分網(wǎng)格，即D＝240，然后根據(jù)式(1)到(3)，計(jì)算構(gòu)造式(4)的等式方程，最后執(zhí)行OMP算法，得到波達(dá)方向估計(jì)結(jié)果。圖4中虛線表示原始信號(hào)方向，實(shí)線表示估計(jì)結(jié)果?？梢钥闯?，兩種算法均能得到很好的估計(jì)效果，然而空間平滑MUSIC算法最多只能估計(jì)出MN＝15個(gè)信號(hào)的來(lái)波方向。對(duì)于上述提到的10元互質(zhì)陣，當(dāng)有16個(gè)來(lái)波信號(hào)時(shí)，空間平滑MUSIC算法便會(huì)出現(xiàn)問(wèn)題，MUSIC算法最多只能估計(jì)出NM＝15個(gè)信號(hào)，當(dāng)16個(gè)信號(hào)時(shí)，方法無(wú)解，若依然按照15個(gè)信號(hào)的輸入求解，估計(jì)結(jié)果與真實(shí)值相差很大，此時(shí)基于OMP算法的波達(dá)方向估計(jì)將展現(xiàn)出優(yōu)勢(shì)。

圖4給出了16個(gè)來(lái)波信號(hào)時(shí)，使用OMP算法進(jìn)行波達(dá)方向估計(jì)的結(jié)果，可以看出OMP算法的波達(dá)方向估計(jì)結(jié)果能給出準(zhǔn)確的估計(jì)結(jié)果，而這種結(jié)果是空間平滑MUSIC算法不能達(dá)到的。

在圖5給出的仿真圖中，本發(fā)明考慮單個(gè)目標(biāo)估計(jì)時(shí)，兩種算法的估計(jì)精度隨著信噪比的變化情況。不同于前面的仿真，這次在-90度到90度的范圍內(nèi)，隨機(jī)產(chǎn)生一個(gè)信號(hào)方向，橫軸表示的信噪比范圍以2dB為步長(zhǎng)，從-20dB到10dB；縱軸表示的誤差采用估計(jì)角度與真實(shí)角度差值的絕對(duì)值。每個(gè)點(diǎn)執(zhí)行5000次蒙特卡洛實(shí)驗(yàn)，其余條件與前面仿真相同。從圖5中的結(jié)果可以看出，對(duì)于單目標(biāo)的波達(dá)方向估計(jì)，在低信噪比時(shí)，兩種算法的估計(jì)性能基本一樣，但是在高信噪比時(shí)，OMP算法的估計(jì)精度優(yōu)于空間平滑的MUSIC算法。