基于四階累量的嵌套陣列波達方向角估計方法
【專利摘要】本發明公開了一種基于四階累量的嵌套陣列波達方向角估計方法，主要解決現有技術中陣元利用率低，信號識別數量少的問題，其實現步驟是：1)分別構建一層嵌套陣列和二層嵌套陣列；2)獲取一層嵌套陣列的輸出信號和二層嵌套陣列的輸出信號；3)根據一層嵌套陣列輸出信號與二層嵌套陣列輸出信號計算所有四階累量；4)根據所有四階累量構造四階累量矩陣；5)計算四階累量矩陣的噪聲子空間；6)根據四階累量的噪聲子空間和陣列流型矩陣計算空間譜；7)根據空間譜繪制幅度譜圖，得到波達方向角。本發明在陣元數量有限的情況下大大提高了陣列可識別的信源數目，適用于目標偵察和無源定位。
【專利說明】
基于四階累量的嵌套陣列波達方向角估計方法
技術領域
[0001] 本發明屬于信號處理技術領域，特別設及一種電磁信號的陣列信號波達方向角估 計方法，可用于對飛機、艦船運動目標的偵察與無源定位。
【背景技術】
[0002] 信號的波達方向角DOA估計是陣列信號處理領域的一個重要分支，它是指利用天 線陣列對空間聲學信號、電磁信號進行感應接收，再運用現代信號處理方法快速準確的估 計出信號源的方向，在雷達、聲納、無線通信等領域具有重要應用價值。隨著科技的不斷進 步，對陣列在進行信號波達方向估計時達到的自由度也有越來越高的要求。
[0003] 針對該問題的研究中，出現較早、應用較為廣泛的是多重信號分類MUSIC子空間的 模型，對于一個L陣元的典型線性均勻陣列，傳統的MUSIC類計算方法可檢測的信源數目是 L-I個。之后的大部分算法都是利用該模型生成的，例如信號參數估計旋轉不變技術 ESPRIT。運些算法采用典型的線性均勻陣列，造成估計的信號數目低于陣元數目，目標個數 很多時甚至無法識別，導致目標捕獲失敗。
[0004] 為了在少的陣元條件下得到盡量大的角度自由度，檢測更多的信源，一些新的陣 列結構被提出，比較有代表性的是嵌套陣列W及互質陣列。Piya Pal等人在其發表的論文 "Nested Arrays=A Novel Approach to Array Processing With Enhanced Degrees of Freedom"(《IEEE transactions on signal processing》，VOL 58,NO.8,August 2010)中 公開了 一種基于嵌套陣列的DOA估計方法，該方法能夠使用M+N個陣元，生成2MN+2N-1個虛 擬陣元，可檢測MN+N-1個信號。該方法具有估計多于陣元數目的信號數的能力，但是，該方 法仍然存在的不足之處是，在陣元數量一定的情況下，能估計的信號數量有時仍不能滿足 實際應用。
[0005] 為解決上述問題，近年來，出現了使用四階累量增加陣列虛擬陣元的方法，可使信 號估計數量得到較大提升。但已有的基于四階累量的波達角度估計方法多數是針對均勻陣 列，而對于嵌套陣列的算法，由于沒有合理設置陣元位置而不能使陣元達到最大利用率。在 實際應用中，給定一定數量的陣元，如果不能合理利用運些陣元獲得足夠多的虛擬陣元，就 不能估計足夠多的信號，造成偵察和定位資源的浪費。
【發明內容】
[0006] 本發明的目的在于針對上述現有技術存在的不足，提出一種基于四階累量的嵌套 陣列波達方向角估計方法，W在陣元數量一定的情況下，提高能夠進行估計的信號數量，避 免因不能合理利用陣元造成的資源浪費。
[0007] 為實現上述目的，本發明技術方案包括如下：
[000引（1)用M+N個天線接收機形成一層嵌套陣列：
[0009] (Ia)將每個天線接收機稱為一個陣元，用M個天線接收機形成第一均勻線性陣列 曰，其陣元間距為d;用N個天線接收機形成第二均勻線性陣列b，其陣元間距為(M+l)d;定義 第一均勻線性陣列a的第一個陣元為起始陣元，其中，1>1，1^>1，0<(1《^/2，^為入射到陣列 的窄帶信號波長；
[0010] (化)將第一均勻線性陣列a與第二均勻線性陣列b組合為一層嵌套陣列:將第二均 勻線性陣列b的第一個陣元放置于與起始陣元相距為Md的位置;將第二均勻線性陣列b的所 有陣元依次插于第一均勻線性陣列a中，形成一層嵌套陣列；
[0011] (Ic)計算一層嵌套陣列的虛擬陣元數A = 2(M+1)N-1;
[0012] (2)用P+Q個天線接收機形成二層嵌套陣列：
[0013] (2a)用P個天線接收機形成第S均勻線性陣列C，其陣元間距為AcU用Q個天線接收 機形成第四均勻線性陣列f，其陣元間距為(P+l)Ad，其中，Q>1，P>1;
[0014] (2b)將第S均勻線性陣列C與第四均勻線性陣列f組合為二層嵌套陣列，即將第S 均勻線性陣列C的第一個陣元放置于與起始陣元相距為(A-l)d的位置，將第四均勻線性陣 列f的第一個陣元放置于與起始陣元相距為[(P+l)A-l]d的位置，并將第四均勻線性陣列f 的所有陣元依次插于第=均勻線性陣列C中，組合為二層嵌套陣列；
[0015] (2c)計算二層嵌套陣列的虛擬陣元數B = 2(P+1)Q-1;
[0016] (3)由一層嵌套陣列和二層嵌套陣列的天線接收機對空間目標信號進行采樣，分 別得到一層嵌套陣列輸出信號Y(t)和二層嵌套陣列輸出信號Z(t)，其中，Y(t) = [yo (t)，…，yi(t)，???，yM+^t)]，yi(t)表示一層嵌套陣列第i個陣元的輸出信號，i的取值范圍是 0《i《M+N，Z(t) = [zo(t)，???，zj(t)，…，zp+Q(t)]，zj(t)表示二層嵌套陣列的第j個陣元的 輸出信號J的取值范圍是0刮《P+Q;
[0017] (4)計算一層嵌套陣列輸出信號Y(t)與二層嵌套陣列輸出信號Z(t)的所有四階累 量c(kl ,k2,k3,k4) = cum(yki(t) ,yk2(t)*,zk3(t) ,zk4(t)*)，
[001引其中，0《kl ,k2《M+N, l《k3 ,k4《P+Q，cum表示求解四階累量操作，（?）*表示向 量的共輛運算；
[0019] (5)定義中間變量D=(AB+2A-l)/2,從所有四階累量C化1，42,43,44)中依次找出 kl，k2，k3，k4滿足kl+k3A-k2-k4A = -D，…，-1，0，1，…，D條件的四階累量，并依次定義為四 階累量元素g(-D)，…，旨(-1)，旨(0)，旨(1)，-'，旨化），將運些四階累量元素重新排列形成四階 累量矩陣G:
[0020]
;
[0021] (6)計算四階累量矩陣G的噪聲子空間Un;
[0022] (7)根據陣列流型矩陣a(0)和四階累量矩陣G的噪聲子空間Un,計算空間譜P(0);
[0023] (8) W波達方向角范圍0的值為X軸坐標，W空間譜P(0)的幅度值為y軸坐標，繪制 幅度譜圖，從該幅度譜圖中按照從高到低的順序尋找幅值較大的前K個譜峰，運些譜峰的峰 值點所對應的X軸坐標即為目標的波達方向角度值，其中，K表示入射到一層嵌套陣列與二 層嵌套陣列的空間目標信號個數，K>1。
[0024] 本發明與現有技術相比具有W下優點：
[0025] 1)本發明采用了雙層嵌套陣列模型進行波達方向角度估計，克服了現有技術中采 用典型的線性均勻陣列，造成估計的信號數目低于陣元數目的缺點，提高了在陣元數目相 同的條件下的陣列可識別信源數目。
[0026] 2)本發明將四階累量應用到嵌套陣列的DOA估計中，通過使用四階累量，設置陣元 位置，可使用M+N+P+Q+1個陣元獲得(2M+2N+1) (2PQ+2Q+1)個連續虛擬陣元，大大提高了陣 列利用率，同時通過獲得更多的虛擬陣元，進一步增加了陣列可識別的信源數目。
[0027] 3)本發明采用了雙層嵌套陣列模型進行波達方向角度估計，相比于使用其他新型 陣列模型對陣元的數目要求更低，提高了陣元數目使用的靈活性。
【附圖說明】
[0028] 圖1是本發明的實現流程圖；
[0029] 圖2是本發明中一層嵌套陣列與二層嵌套陣列的結構示意圖。
【具體實施方式】
[0030] W下參照附圖，對本發明的技術方案和效果作進一步的詳細說明。
[0031 ]參附圖1，本發明的具體步驟如下：
[0032] 步驟1:用M+N個天線接收機形成一層嵌套陣列，得到其一層嵌套陣列虛擬陣元數 A。
[0033] (Ia)將每個天線接收機稱為一個陣元，用M個天線接收機形成第一均勻線性陣列 曰，其陣元間距為d;用N個天線接收機形成第二均勻線性陣列b，其陣元間距為(M+l)d;定義 第一均勻線性陣列a的第一個陣元為起始陣元，其中，1>1，1^>1，0<(1《^/2，^為入射到陣列 的窄帶信號波長；
[0034] (化)將第一均勻線性陣列a與第二均勻線性陣列b組合為一層嵌套陣列:將第二均 勻線性陣列b的第一個陣元放置在與起始陣元相距為Md的位置;將第二均勻線性陣列b的所 有陣元依次插于第一均勻線性陣列a中，形成一層嵌套陣列；
[0035] (Ic)計算一層嵌套陣列的虛擬陣元數A = 2(M+1)N-1。
[0036] 步驟2:用P+Q個天線接收機形成二層嵌套陣列，得到二層嵌套陣列的虛擬陣元數 B。
[0037] (2a)用P個天線接收機形成第S均勻線性陣列C，其陣元間距為AcU用Q個天線接收 機形成第四均勻線性陣列f，其陣元間距為(P+l)Ad，其中，Q>1，P>1;
[0038] (2b)將第S均勻線性陣列C與第四均勻線性陣列f組合為二層嵌套陣列，即將第S 均勻線性陣列C的第一個陣元放置于與起始陣元相距為(A-l)d的位置，將第四均勻線性陣 列f的第一個陣元放置于與起始陣元相距為[(P+l)A-l]d的位置，并將第四均勻線性陣列f 的所有陣元依次插于第=均勻線性陣列C中，組合為二層嵌套陣列；
[0039] (2c)計算二層嵌套陣列的虛擬陣元數B = 2(P+1)Q-1。
[0040] 步驟3:獲得一層嵌套陣列輸出信號和二層嵌套陣列輸出信號。
[0041] 由一層嵌套陣列和二層嵌套陣列的天線接收機對空間目標信號進行采樣，分別得 到一層嵌套陣列輸出信號Y(t)和二層嵌套陣列輸出信號Z(t)，其中，Y(t) = [yo(t)，…，yi (t)，???，yM+N(t)]，yi(t)表示一層嵌套陣列第i個陣元的輸出信號，i的取值范圍是0《i《M+ N，Z(t) = [Z0(t)，…，Zj(t)，…，ZP+Q(t)]，Zj(t)表示二層嵌套陣列的第j個陣元的輸出信號， j的取值范圍是〇《j《P+Q。
[0042] 步驟4:計算輸出信號的所有四階累量。
[0043] 計算一層嵌套陣列輸出信號Y(t)與二層嵌套陣列輸出信號Z(t)的所有四階累量
[0044] c(kl ,k2,k3,k4) = cum(yki(t) ,yk2(t)*,zk3(t) ,zk4(t)*)，
[0045] 其中，0《kl ,k2《M+N, l《k3 ,k4《P+Q,cum表示求解四階累量操作，（?）*表示向 量的共輛運算。
[0046] 步驟5:構造四階累量矩陣G。
[0047] 定義中間變量D= (AB+2A-1 )/2，從所有四階累量C化1，k2，k3，k4)中依次找出kl， k2，k3，k4滿足kl+k3A-k2-k4A =-D，…，-I，0，I，…，D條件的四階累量，并依次定義為四階累 量元素g(-D)，…，g(-l)，g(0)，g(l)，…，g(D)，將運些四階累量元素重新排列，形成四階累 量矩陣G:
[004引
[0049] 步驟6:計算四階累量矩陣的噪聲子空間Un。
[0050] (6a)對四階累量矩陣G進行如下特征分解：
[0051] G = U-A-IjH,
[0052] 其中，A為四階累量矩陣（.）H的特征值矩陣，U為四階累量矩陣G的特征值所對應 的特征向量矩陣，（?）H表示矩陣的共輛轉置遠算；
[0053] (6b)將特征值矩陣A中的特征值按從大到小排序，取后D-K個較小特征值對應的 特征向量矩陣作為噪聲子空間Un。
[0054] 步驟7:計算空間譜P(目）。
[0055] (7a)將觀測空域[-90°，90° ]按等間隔L劃分成E個角度，定義波達方向角范圍為0 =[目1，目2，'。，目6,…目E]，其中，目e表示第e個目標波達方向角，e = l，2,…，E，E〉〉max{M，N，P, QhL的取值根據期望達到的角度估計精度進行設定，網格劃分間隔越小則最終得到的角度 估計值精度越高；
[0056] (7b)計算第e個目標波達方向角06對應的導向矢量a(0e):
[0化7]
[005引其中，（?）嗦示矩陣轉置運算，j為虛數單位；
[0059] (7c)計算所有的目標波達方向角對應的導向矢量，得到陣列流型矩陣a(0):
[0060] a(目）= [a(目 1),a(目e),a(目E)];
[0061] (7d)根據陣列流型矩陣a(0)和四階累量矩陣G的噪聲子空間Un，計算空間譜P(0):
[0062]
[0063] 其中，（?）H表示矩陣的共輛轉置遠算。
[0064] 步驟8:繪制幅度譜圖。
[0065] (8a似波達方向角范圍0的值為X軸坐標，W空間譜P(0)的幅度值為y軸坐標，繪制 幅度譜圖；
[0066] (8b)從繪制的幅度譜圖中，按照從高到低的順序將所有譜峰幅度值進行排序，將 前K個譜峰的峰值點所對應的X軸坐標做為目標的波達方向角度值，其中，K表示入射到一層 嵌套陣列與二層嵌套陣列的空間目標信號個數，且假設空間目標信號在傳播過程中加入了 均值為零的復高斯白噪聲，K^l。
[0067] 下面結合仿真實例對本發明的效果做進一步的描述。
[0068] 仿真實例1.構建雙層嵌套陣列結構。
[0069] 1.1)利用第一均勻線性陣列a與第二均勻線性陣列b形成一層嵌套陣列，設N = 2，M =2,第一均勻線性陣列a含有3個陣元，陣元間距為d，第二均勻線性陣列b含有3個陣元，陣 元間距為3d，第一均勻線性陣列a的第一個陣元放置在Id位置；
[0070] 1.2)計算第一均勻線性陣列a與第二均勻線性陣列b形成的一層嵌套陣列虛擬陣 元數 A = 2MN+2N-1 = 11;
[0071] 1.3)利用第S均勻線性陣列C與第四均勻線性陣列f形成二層嵌套陣列，設P = 2，Q =2,第S均勻線性陣列C含有2個陣元，陣元間距為lld，第四均勻線性陣列f含有2個陣元， 陣元間距為33d，第S均勻線性陣列C的第一個陣元在Ild的位置，計算虛擬陣元數B = 2PQ+ 2Q-1 = 11。
[0072] 其中第一均勻線性陣列a、第二均勻線性陣列b和第S均勻線性陣列c、第四均勻線 性陣列f形成了雙層嵌套陣列，其陣列結構圖如圖2所示。
[0073] 由圖2可見，嵌套陣列結構在陣元數量一定的情況下，可W得到更多的陣元位置信 息，從而增加陣列可識別信源數目，同時，嵌套陣列結構相比于使用其他新型陣列模型對陣 元的數目要求更低，提高了陣元數目使用的靈活性。
[0074] 仿真實例2,根據四階累量計算雙層嵌套陣列所有虛擬陣元。
[0075] 2.1)設第一均勻線性陣列a與第二均勻線性陣列b形成的一層嵌套陣列陣元位置 為[l，2,3,6]d，；
[0076] 2.2)設第S均勻線性陣列C與第四均勻線性陣列f形成的二層嵌套陣列陣元位置 為[ll，22,33,66]d;
[0077] 2.3)計算得到一層嵌套陣列每個陣元分別與二層嵌套陣列中每個陣元形成的虛 擬陣列為[12,13,15，17]d、[23,24,26，28]d、[34,35,36，39]d、[67,68,70，72]d;
[0078] 2.4)根據四階累量計算的方法和原則，用上述運些虛擬陣列生成最終的所有虛擬 陣元，結果如表1。
[0079] 表1雙層嵌套陣列所有虛擬陣元
[00801
[0081 ] 從表1可見，該雙層嵌套陣列用M+N+P+Q = 8個原始陣元生成了一個從-60d到60d共 含有AB = 121個虛擬連續陣元的線性陣列。
[0082]綜上，本發明解決了現有技術陣元利用率低，識別信源數目少，無源定位估計誤差 大的問題，降低了對陣元數目的要求，保證了陣元數目使用的靈活性，提高了一定陣元數情 況下陣列可識別的信源數目W及低信噪比下對信號方向角的估計性能。
【主權項】
1. 一種基于四階累量的嵌套陣列的波達方向角估計方法，其特征在于，包括以下步驟： (1) 用M+N個天線接收機形成一層嵌套陣列： (la) 將每個天線接收機稱為一個陣元，用M個天線接收機形成第一均勻線性陣列a，其 陣元間距為d;用N個天線接收機形成第二均勻線性陣列b，其陣元間距為(M+l)d;定義第一 均勻線性陣列a的第一個陣元為起始陣元，其中，M多1，N多1，0〈(1<λ/2，λ為入射到陣列的窄 帶信號波長； (lb) 將第一均勻線性陣列a與第二均勻線性陣列b組合為一層嵌套陣列:將第二均勻線 性陣列b的第一個陣元放置于與起始陣元相距為Md的位置;將第二均勻線性陣列b的所有陣 元依次插于第一均勻線性陣列a中，形成一層嵌套陣列； (lc) 計算一層嵌套陣列的虛擬陣元數A = 2 (M+1)N-1; (2) 用P+Q個天線接收機形成二層嵌套陣列： (2a)用P個天線接收機形成第三均勻線性陣列c，其陣元間距為AcU用Q個天線接收機形 成第四均勻線性陣列f，其陣元間距為(P+l)Ad，其中，Q多1，P彡1; (2b)將第三均勻線性陣列c與第四均勻線性陣列f組合為二層嵌套陣列，即將第三均勻 線性陣列c的第一個陣元放置于與起始陣元相距為(A-l)d的位置，將第四均勻線性陣列f的 第一個陣元放置于與起始陣元相距為[(P+l)A-l]d的位置，并將第四均勻線性陣列f的所有 陣元依次插于第三均勻線性陣列c中，組合為二層嵌套陣列； (2c)計算二層嵌套陣列的虛擬陣元數B = 2(P+1)Q-1; (3) 由一層嵌套陣列和二層嵌套陣列的天線接收機對空間目標信號進行采樣，分別得 到一層嵌套陣列輸出信號Y(t)和二層嵌套陣列輸出信號Z(t)，其中，Y(t) = [y〇(t)，H_，yi (t)，···，yM+N(t)]，yi(t)表示一層嵌套陣列第i個陣元的輸出信號，i的取值范圍是(Xi彡M+ N，Z(t) = [ZO(t)，···，Zj(t)，…，ZP+Q(t)]，Zj(t)表示二層嵌套陣列的第j個陣元的輸出信號， 的取值范圍是〇彡j彡P+Q; (4) 計算一層嵌套陣列輸出信號Y(t)與二層嵌套陣列輸出信號Z(t)的所有四階累量c (kl，k2，k3，k4) = cum(yki(t)，yk2(t)*，zk3(t)，zk4(t)*)，其中，0<kl，k2<M+N，Kk3，k4<P +Q，cum表示求解四階累量操作，（· Γ表示向量的共輒運算； (5) 定義中間變量D = (AB+2A-1)/2，從所有四階累量c (kl，k2，k3，k4)中依次找出kl， k2，k3，k4滿足kl+k3A-k2-k4A =-D，…，-I，O，I，…，D條件的四階累量，并依次定義為四階累 量元素 g(_D)，…，g(-l)，g(0)，g(l)，~，g(D)，將這些四階累量元素重新排列形成四階累量 矩(6) 計算四階累量矩陣G的噪聲子空間Un; (7) 根據陣列流型矩陣α(θ)和四階累量矩陣G的噪聲子空間Un，計算空間譜Ρ(θ); (8) 以波達方向角范圍Θ的值為X軸坐標，以空間譜Ρ(θ)的幅度值為y軸坐標，繪制幅度 譜圖，從該幅度譜圖中按照從高到低的順序尋找幅值較大的前K個譜峰，這些譜峰的峰值點 所對應的X軸坐標即為目標的波達方向角度值，其中，K表示入射到一層嵌套陣列與二層嵌 套陣列的空間目標信號個數，K多1。2. 根據權利要求1所述的基于四階累量的嵌套陣列波達方向角估計方法，其中步驟(6) 中計算四階累量矩陣G的噪聲子空間Un，按如下步驟進行： (6a)對四階累量矩陣G進行如下特征分解： 其中，Λ為四階累量矩陣（·）H的特征值矩陣，U為四階累量矩陣G的特征值所對應的特 征向量矩陣，（·）H表示矩陣的共輒轉置遠算； (6b)將特征值矩陣Λ中的特征值按從大到小排序，取后D-K個較小特征值對應的特征 向量矩陣作為噪聲子空間Un。3. 根據權利要求1所述的基于四階累量的嵌套陣列波達方向角估計方法，其中步驟(7) 中的陣列流型矩陣α (Θ )，按如下步驟構造： (7a)將觀測空域[-90°，90°]等間隔劃分成E個角度，定義為波達方向角范圍Q = M1, Θ2,…，9e，…ΘΕ]，其中，0e表示第e個目標波達方向角， e = l，2,…，E，E?max{M，N，P，Q}; (7b)計算第e個目標波達方向角對應的導向矢量Ci(Qe3):其中，（·）T表示矩陣轉置運算，j為虛數單位； (7c)計算所有的目標波達方向角對應的導向矢量，得到陣列流型矩陣α(θ): α(θ) = [α(θι) ,---,Q(Oe)4. 根據權利要求1所述的基于四階累量的嵌套陣列波達方向角估計方法，其中步驟(7) 中計算空間譜Ρ(θ)，按如下形式進行：其中，α(θ)表示陣列流型矩陣，（·）Η表示矩陣的共輒轉置運算。
【文檔編號】G01S3/14GK106019215SQ201610560482
【公開日】2016年10月12日
【申請日】2016年7月15日
【發明人】蔡晶晶, 宗汝, 趙晗希, 武斌, 李鵬, 蘇瑤
【申請人】西安電子科技大學