本發(fā)明涉及一種基于互質陣列的寬帶信號DOA估計方法，屬于陣列信號處理領域。

背景技術：

陣列信號處理技術被廣泛應用于國防和人們生活的方方面面。寬帶陣列信號的波達方向DOA估計一直是陣列信號處理領域中的一個引起廣泛興趣的問題，在雷達、通信、智能天線和無源定位等方面具有較多的應用。現(xiàn)有的寬帶DOA估計方法主要分為兩大類：一類主要是最大似然(maximum likelihood，ML)參數(shù)估計方法，該方法不僅需要已知或者預先估計出信源的個數(shù)，而且其代價函數(shù)往往呈現(xiàn)多峰特性，全局收斂性差，其求解的過程需要涉及多維搜索，運算量較大；另一類算法是基于信號子空間的方法，該類算法因為具有相對較低的運算復雜度，一直得到陣列信號處理領域的學者和工程技術人員的廣泛關注，目前已發(fā)展出眾多的基于信號子空間的寬帶陣列信號DOA估計算法。另外，基于信號子空間的這一類方法，在低信噪比情況下，由于子空間估計得不準確，導致測向誤差較大、精度較低；此外，在相同陣元數(shù)目的情況下，陣列的可識別信源數(shù)目也較少。

在實際中，并不是在整個空間內任意位置處都存在信號源，這就意味著信號在整個空域中的能量分布具有稀疏性，即空間譜是稀疏的，正如單個正(余)弦信號在整個頻域范圍內具有稀疏性一樣。正是由于信源空間分布的稀疏性才使得利用稀疏信號表示理論實現(xiàn)空間譜估計提供了可能。由于基于信號稀疏表示的DOA估計方法在理論上具有不需要信號的統(tǒng)計特性、無需初始化、解相干等諸多優(yōu)點，因而利用信源的空域稀疏性，可將信源的DOA估計轉化為一個稀疏信號重構的問題。稀疏信號重構可以為信號處理領域中的許多問題提供一種新的解決方法，目前越來越引起信號處理學界的廣泛興趣，特別是近年來新興的壓縮感知的發(fā)展進一步推動了稀疏信號重構的應用。目前，將稀疏重構應用于窄帶DOA估計已經獲得了較多的研究。比如申請?zhí)枮椤癈N201410172438.8”的專利申請公開了“一種基于樣本協(xié)方差矩陣稀疏性的波達方向估計方法”，涉及陣列信號處理領域。主要解決現(xiàn)有技術中需要已知噪聲功率，而且計算量很大的問題。申請?zhí)枮椤癈N201410280744.3”公開了一種基于空間平滑協(xié)方差矩陣稀疏表示的波達方向角估計方法，主要解決現(xiàn)有技術運算量大，處理相干信號源性能差，造成無源定位估計誤差大的問題。但是這些DOA估計算法都只能適用于窄帶信號；而且都要求陣列滿足空間Nyquist采樣定理，即相鄰陣元間距不大于入射信號半波長，以保證DOA估計的唯一性，不過這樣的陣列又會存在陣列孔徑小、測向精度低和分辨率差等缺點，已不能滿足日益復雜的現(xiàn)代戰(zhàn)爭環(huán)境和競爭日益激烈 的通信市場。

由于稀疏陣列的高自由度，而且優(yōu)于均勻陣列的特點，稀疏陣列在近幾年的發(fā)展步伐很快，尤其是互質陣列?；ベ|陣列是一種陣元間距大于半波長非均勻陣列系統(tǒng)。它擁有更大的陣列孔徑，甚至包含多尺度陣元間距?；ベ|陣列的這些特點為其提高測向性能奠定了物理基礎，使得進一步改善測向性能或節(jié)約成本成為可能。和常規(guī)陣列相比，互質陣具有更多的優(yōu)勢，主要體現(xiàn)為如下三個方面：1)陣元數(shù)相同時擁有更大的陣列孔徑，甚至更多的自由度，改善了測向算法在測向精度、測向分辨率和最大可處理信源數(shù)等方面的性能；2)陣列孔徑相同時所需陣元數(shù)更少，更少的陣元數(shù)不僅意味著更少的天線單元，而且還意味著更小規(guī)模的接收系統(tǒng)和信號處理系統(tǒng)等，極大地降低設備成本；3)減小陣元間互耦，削弱其對測向性能的影響?？傊@些優(yōu)勢使互質陣列成為重點研究內容。如申請?zhí)枮椤癈N201510126041.X”的專利申請公開了“一種基于稀疏重構的互質陣列波達方向角估計方法”，該方法通過在空域波達方向角觀測空間進行網格劃分，構造超完備基Φ(θ)，并定義一個空域稀疏向量s，然后將DOA估計轉化為一種稀疏約束方程的求解，主要解決現(xiàn)有技術運算量大、需要先驗知識多、識別信源數(shù)目少以及DOA估計誤差大的問題，實現(xiàn)了在低先驗知識的情況下提高了無源測向的精度和運算速度，在相同陣元數(shù)目時，增加了陣列可識別的信源數(shù)目，提高了低信噪比下對信號方向角估計精度，可用于目標偵察和無源定位。但是該方法中進行一次差分操作后，僅僅利用兩個空穴之間這段連續(xù)的虛擬線性陣列進行DOA估計，導致在相同陣元數(shù)目時，陣列可識別的信源數(shù)目以及DOA估計的精確仍有待提高；另外，該方法也只是針對窄帶信號，并不能適用于寬帶信號，已經不能夠支持現(xiàn)在越來越復雜的電子通信環(huán)境。隨著通信技術的發(fā)展，跳頻信號、擴頻信號、線性調頻信號等寬帶信號在通信系統(tǒng)中的應用越來越普遍，另外，在自然界中還有些多信號本質上就屬于寬帶信號，如聲音信號、地震波等。因此，研究寬帶信號測向就顯得越來越重要。另外無論是在民用領域還是軍用領域，普通陣列(ULA)也已經不能滿足復雜的電子系統(tǒng)環(huán)境，且其高成本在實施的過程中也帶來了一定的困難。因此，仍需人們繼續(xù)進行研究改進。

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

本發(fā)明的目的在于，提供一種基于互質陣列的寬帶信號DOA估計方法，它可以有效解決現(xiàn)有技術中存在的問題，尤其是在低信噪比情況下，進行寬帶DOA估計時，測向誤差較大、精度較低的問題。

為解決上述技術問題，本發(fā)明采用如下的技術方案：基于互質陣列的寬帶信號DOA估 計方法，包括以下步驟：

S1，利用天線設計互質陣列結構，每個天線稱為一個陣元；

S2，對互質陣列中的天線接收到的寬帶信號進行抽樣及離散傅里葉變換，獲得頻域信號輸出模型X[l]；

S3，計算頻域信號輸出模型X[l]的自相關矩陣并將其矢量化，得新的信號模型Z[l]；

S4，對所述新的信號模型Z[l]進行去冗余、去相關處理，并獲得寬帶信號的空間平滑協(xié)方差矩陣R′；

S5，劃分空域網絡，構造字典，并利用寬帶信號的多個頻率點處的字典對所述的空間平滑協(xié)方差矩陣R′進行稀疏表示，形成寬帶信號多個字典的多測量矢量稀疏表示模型(即MMV模型)；

S6，通過對多字典稀疏表示系數(shù)的聯(lián)合稀疏約束以求解稀疏反問題的形式實現(xiàn)寬帶信號的波達方向估計。

前述的基于互質陣列的寬帶信號DOA估計方法中，步驟S1中所述的互質陣列為傳統(tǒng)互質陣列。

具體的，所述的傳統(tǒng)互質陣列采用下述方法進行設計：

(1)用N個天線形成均勻線性陣列1，簡稱線陣1；用M-1個天線形成均勻線性陣列2，簡稱線陣2；其中，線陣1的陣元間距為Md，線陣2的陣元間距為Nd，N＞M≥2且M與N互質，0＜d≤λ/2，λ為入射到互質陣列的寬帶信號的波長；

(2)組合線陣1和線陣2為互質陣列：設線陣1的第一個陣元為互質陣列的陣元0，將線陣2的第一個陣元放置于與陣元0相距為Nd的位置，依次進行處理，使得線陣2的所有陣元依次插于線陣1中；從互質陣列的陣元0開始，從頭至尾依次命名各個陣元為陣元0，陣元1，……，陣元M+N-1。

優(yōu)選的，步驟S1還包括：對形成互質陣列的各個天線的位置進行二次差分操作，形成大孔徑均勻線性虛擬陣列，從而可以實現(xiàn)在相同陣元數(shù)目時，獲得更大的陣列自由度，增加陣列可識別的信源數(shù)目，實現(xiàn)更精確的DOA估計。具體的說，現(xiàn)有技術中通過對上述方法中由線性陣列1和線性陣列2構成的互質陣列進行一次差分操作后形成的虛擬線性陣列自由度較低，因而在DOA估計時，相應的譜分辨率也較低。而本發(fā)明通過對經過一次差分后形成的虛擬線性陣列進行第二次差分操作可以獲得更大的自由度，從而可以實現(xiàn)在相同陣元數(shù)目和相同的快拍數(shù)下，增加陣列可識別的信源數(shù)目，降低DOA估計的誤差，同時提高DOA估計時的譜分辨率，而且成本較低。其進一步形成的虛擬均勻線性陣列天線位置為 (-MN-M-N+1)d到(MN+M+N-1)d。

步驟S2中，假設K個寬帶信號入射到由M+N-1個天線構成的互質陣列中，其來波方向為θ＝[θ₁，θ₂，...，θ_K]^T，用一個抽樣頻率f_s對線陣1和線陣2接收到的信號進行抽樣，得到兩個線陣接收信號的離散形式(因為每個線陣接收的信號在時間上是連續(xù)的，當采用一個抽樣頻率f_s對線陣1和線陣2接收到的信號進行抽樣時，在時間上就是離散的，故為離散形式)；對所述離散形式的信號進行L點DFT變換(由于寬帶測向不能在時域進行，經過DFT變換后，天線接收到的信號由時域變換到了頻域；其中，L的值與寬帶信號的帶寬、抽樣頻率的設定相關)，獲得頻域的信號輸出模型X[l]。

優(yōu)選的，步驟S3具體包括：計算頻域信號輸出模型X[l]的自相關矩陣R_xx[l]并利用矢量化vec(·)函數(shù)對該自相關矩陣進行按列堆棧，得到一個新的向量Z[l]，再將Z[l]的協(xié)方差矩陣進行矢量化，獲得二次差分互質陣列的信號模型Z₁[l]。

優(yōu)選的，步驟S3包括：計算頻域信號輸出模型X[l]的自相關矩陣并將其矢量化，得新的信號模型Z[l]，再將Z[l]的協(xié)方差矩陣進行矢量化，獲得二次差分互質陣列的信號模型Z₁[l]；步驟S4包括：利用空間平滑和子陣選取技術對所述二次差分互質陣列的信號模型Z₁[l]進行去冗余、去相關處理，并獲得寬帶信號的空間平滑協(xié)方差矩陣R′；具體包括以下步驟：

S41，從新的信號模型Z₁[l]的流型矩陣B_l^*e B_l中建立一個(2MN+2M+2N-1)×K的新的流型矩陣B_1l，也即從流型矩陣B_l^*e B_l中取出不同行對應的2MN+2M+2N-1個連續(xù)的差分值，并對取出的連續(xù)的差分值進行排序；(這相當于從新信號模型Z₁[l]中移除了相應的重復行并進行排序得到一個新的向量Z₂[l]，新形成的流型矩陣B_1l等于一個帶有2MN+2M+2N-1根天線的虛擬均勻線性陣列的流型矩陣)其中，M表示形成互質陣列的子陣1的天線數(shù)目，N表示形成互質陣列的子陣2的天線數(shù)目，K表示入射到互質陣列的信號數(shù)；

S42，將新的流型矩陣B_1l中的2MN+2M+2N-1根天線劃分成NM+M+N個子陣列，劃分后的每個子陣列中有NM+M+N根天線；其中，第i個子陣列天線的位置位于{(-i+1+n)d,n＝0,1,...,NM+M+N-1}；

S43，設第i個子陣列對應新向量Z₂[l]的第(NM+M+N+1-i)行到第(2MN+2M+2N+1-i)行，從而獲得流型矩陣B_1l中每個子陣列對應的新信號模型Z_2i[l]，計算新信號模型Z_2i[l]的協(xié)方差矩陣，進而獲得寬帶信號的空間平滑協(xié)方差矩陣R′。

通過上述方法進行去冗余、去相關處理，從而可以增加從新的向量Z[l]中獲得的半定矩陣的秩，實現(xiàn)更精確地對寬帶信號DOA進行估計。

上述去冗余、去相關處理方法中，與原來的寬帶信號頻域輸出模型X[l]相比，新的信號模型Z₁[l]也類似于一個天線陣列的接收信號模型，其流型矩陣通過B_l^*e B_l給出，等源信號向量表示為q,它們類似于相關源，其中流型矩陣B_l^*e B_l中不同的行對應著一個擁有更多陣元數(shù)目的天線陣列接收信號模型的流型矩陣所對應的行(所述的擁有更多陣元數(shù)目的天線陣列的天線位置在所有d的整數(shù)倍位置上，范圍從(-NM-M-N+1)d到(NM+M+N-1)d)。

前述方法中，步驟S5具體包括以下步驟：

首先，利用寬帶信號的頻率w_l處的字典對所述的空間平滑協(xié)方差矩陣R′的第q列進行稀疏表示，得單個快拍矢量(SMV矢量)模型的稀疏表示為：

式中，Θ＝{θ₁,θ₂,...,θ_Q}集合表示整個DOA空間域的搜索網格，共有Q(Q>>K，K表示入射到互質陣列的信號數(shù))個柵格點；B_1l(Θ)表示整個空間內單個頻率點w_l處的流型矩陣，稱為過完備基或字典；e_q是一個(NM+M+N)×1維的列向量，它表示噪聲，它的第q行元素是1，其余位置元素均是0；b_lq是Q×1維列稀疏向量，表示字典B_1l(Θ)的稀疏表示系數(shù)，為單個列稀疏向量；為噪聲功率；(若稀疏表示系數(shù)b_lq的第q行為非零元素，則意味著對應θ_q＝θ_k(k＝1,2,...,K)處存在一個信號源，若θ_q不在柵格點上，則θ_q≈θ_k)

其次，將NM+M+N個快拍矢量(SMV矢量)模型合成一個矩陣，即得整個空間平滑協(xié)方差矩陣R′的稀疏表示，形成寬帶信號多個字典的多測量矢量稀疏表示模型(即MMV模型)：

這里只是寫單個頻率點的稀疏表示，因為l的值是變化的，實際上是利用寬帶信號的多個頻率點處的字典對所述的空間平滑協(xié)方差矩陣R′進行稀疏表示；也即所述的SMV稀疏表示模型只是對寬帶信號空間平滑協(xié)方差矩陣的第q列進行稀疏表示，而將NM+M+N個SMV稀疏表示模型合成一個矩陣，形成MMV稀疏表示模型，即是對整個寬帶信號空間平滑協(xié)方差矩陣的稀疏表示；

其中，G_l＝[b_l1,b_l2,...,b_l(NM+M+N)]，由于G_l中的每列稀疏向量都具有共同的稀疏結構， 因此G_l具有聯(lián)合稀疏性，I₁表示一個(NM+M+N)×(NM+M+N)的單位矩陣。

前述方法中，步驟S6具體包括以下步驟：

S61，將寬帶信號DOA參數(shù)估計問題表述為單個頻率點w_l處的多個單字典稀疏約束優(yōu)化模型(求解此稀疏問題相當于求解L個單字典稀疏優(yōu)化模型進行聯(lián)合稀疏處理(由下面的公式可知))：

式中，是G_l的估計；μ(G_l)表示矩陣G_l的非零行對應的索引集合，|μ(G_l)|＝K表示G_l中非零行的個數(shù)(即稀疏度)，K表示入射到互質陣列的信號數(shù)；

S62，根據所述的稀疏約束優(yōu)化模型，將其轉化為凸優(yōu)化問題進行求解，從而獲得寬帶信號的DOA估計。

本發(fā)明從壓縮感知稀疏恢復理論出發(fā)，將稀疏約束優(yōu)化模型轉化為最終的凸優(yōu)化問題，通過CVX進行求解，從而可以獲得更好的DOA估計結果以及更高的分辨率。

優(yōu)選的，步驟S61中，當多個字典存在時，利用所述的稀疏約束優(yōu)化模型構造一個多字典聯(lián)合優(yōu)化問題：

式中，為無混疊的約束條件。

優(yōu)選的，步驟S62具體包括以下步驟：(所述的多字典聯(lián)合優(yōu)化問題相當于組合優(yōu)化問題，其為NP難問題)將所述的多字典聯(lián)合優(yōu)化問題轉化為1范數(shù)凸優(yōu)化BP問題進行求解(將NP難問題進行轉化)：

首先，(約束條件相當于同一稀疏結構約束，為了對其進行聯(lián)合稀疏約束)設待估計變量的行之間具有稀疏性，而列之間不具有稀疏性，同時假設b^ο中的每一個元素為矩陣對應該行向量的2范數(shù)，即 而||b^ο||₁表示只對在空域θ上進行稀疏約束；因此，對應的BP優(yōu)化模型為：

其中，||b^ο||₁表示對所有字典的表示系數(shù)的聯(lián)合行稀疏約束(這是根據壓縮感知信號恢復理論而得，由前面知，待估計變量同時b^ο中的每一個元素為矩陣對應該行向量的2范數(shù)，即這里b^ο中包含所有字典的稀疏表示系數(shù))；

其次，將所述的BP優(yōu)化模型用對應的正則化模型表示為：

式中，δ為正則化參數(shù)，表示Frobenius范數(shù)；如果能求得或(繼續(xù)轉化為下面標準的SOCP形式，根據CVX求解)，則通過其非零行對應的位置就可估計參數(shù)θ；

最后，所述的正則化模型為凸優(yōu)化問題，對其進行求解，獲得寬帶信號的DOA估計。

前述的基于互質陣列的寬帶信號DOA估計方法中，將所述的正則化模型轉化為標準的SOCP形式(便于求解)：

式中，1_Q是一個所有元素為1的Q×1維列向量；z和η為引入的輔助變量，γ是一個1×Q的向量，其中第q個元素是γ_q，目標函數(shù)為線性函數(shù)；

采用數(shù)值軟件包sedumi或者CVX對所述的標準的SOCP形式進行求解，進而獲得寬帶信號的DOA估計。

與現(xiàn)有技術相比，本發(fā)明具有以下優(yōu)點：

1、本發(fā)明通過對每個頻率點處的寬帶信號空間平滑協(xié)方差矩陣建立字典進行稀疏表示，建立多個字典的多測量矢量(Multiple-Measurement Vectors，MMV)稀疏表示模型，然后利用多個字典的聯(lián)合稀疏約束實現(xiàn)寬帶信號的波達方向估計，不僅可以使用稀疏陣列直接處理更多數(shù)目的寬帶信號，提高其波達方向估計的分辨率，而且可以獲得較好的檢測性能(進行寬帶DOA估計時，測向的均方根誤差較低)；

2、本發(fā)明中通過對由M+N-1根天線形成的互質陣列進行二次差分操作，形成的虛擬均勻線性陣列結構中，不同的連續(xù)差分值數(shù)目為2MN+2M+2N-1個，其陣列自由度也為2MN+2M+2N-1，因此可以在相同陣元數(shù)目下，更多的增加陣列可識別的信源數(shù)目；另 外，使用子陣選取和空間平滑技術后，在形成的2MN+2M+2N-1根天線的虛擬均勻線性陣列中可利用的天線數(shù)目是MN+M+N，其自由度也為MN+M+N(說明采用本發(fā)明的方法進行DOA估計時，可利用的陣列自由度為MN+M+N，明顯多于互質陣列經過一次差分后獲得的自由度)；從而可以識別更多的天線陣源數(shù)目，同時獲得更好的檢測性能(尤其在低信噪比情況下，進行寬帶DOA估計時，測向的均方根誤差更低；另外，可以準確的估計所有的來波方向，既無偽峰，又能表現(xiàn)出更高的譜分辨特性；具有更高的角度分辨率，可以提高低信噪比下對寬帶信號方向角估計精度，降低測向誤差)；

3、在陣元數(shù)目相同的情況下，使用本發(fā)明的方法可以獲得更大的自由度，檢測相同信號數(shù)目時，可更多的增加陣列可識別的信源數(shù)目，而且具有更好的檢測性能(隨著信噪比的變化表現(xiàn)出了更低的均方根誤差特性，特別是在低SNR條件下，這種優(yōu)勢更加明顯)；另外，在相同的陣元數(shù)目和相同的快拍數(shù)下，使用本發(fā)明的均方根誤差明顯低于其他方法，同時，由于本發(fā)明的方法可以獲得更大的陣列自由度，因而在DOA估計時，具有更好的譜分辨率，且系統(tǒng)成本也更低(因為沒有進行互質陣列陣元數(shù)目的擴展，而僅僅使用互質陣列的數(shù)目是M+N-1)，此外，本發(fā)明需要的先驗知識也較少。

4、發(fā)明人經大量試驗論證，相對于廣義互質陣列，本發(fā)明采用傳統(tǒng)互質陣列進行二次差分處理獲得的效果更好，可以處理更多數(shù)目的寬帶信號，提高其波達方向估計的分辨率，而且可以獲得更好的檢測性能。

為了驗證上述效果，發(fā)明人進行了以下試驗論證：

首先，假設構造的互質陣列結構，M＝3,N＝5，給定的信源數(shù)目K＝10，其信號的來波方向分別為θ＝[-60°,-45°,-30°,-20°,-10°,5°,20°,30°,45°,60°]，其信噪比(SNR)用公式表示為和分別表示信號功率和噪聲功率，這里設定SNR＝5dB，快拍數(shù)為512，實驗選取d＝0.05m滿足空域采樣無模糊條件，θ＝[-90°,90°]以步長1°為間隔，本次實驗方法選取寬帶空間平滑MUSIC算法與本發(fā)明中的方法進行比較，兩種方法的測向實驗結果如圖7所示，由圖7可以看出，在給定的來波方向中，寬帶空間平滑MUSIC算法基本可以分辯大部分的來波方向，但是也存在一些方向上分辨率不好，不能真實估計出的來波方向，致使真實DOA參數(shù)的位置無法判別；而采用本發(fā)明中的方法(即互質陣列二次差分操作和寬帶信號空間平滑協(xié)方差矩陣稀疏表示方法)則可以準確的估計所有的來波方向，既無偽峰，又能表現(xiàn)出更高的譜分辨特性。

其次，假設有兩個遠場寬帶信號從波達方向-30°和-30°+Δθ入射到由8陣元組成的互質陣列中，快拍數(shù)為256，信噪比為5dB，其中Δθ從3°到10°，以1°為間隔進行變化，每個 角度間隔點做100次蒙特卡羅實驗。圖8給出了寬帶信號空間平滑MUSIC方法和本發(fā)明的方法(即互質陣列二次差分操作和寬帶信號空間平滑協(xié)方差矩陣稀疏表示方法)隨角度間隔變化的均方誤差曲線。從圖8可以看出，在給定的仿真條件下，相對于寬帶信號空間平滑MUSIC算法，本發(fā)明的方法具有更高的角度分辨率。

最后，比較使用三種不同的陣列結構實現(xiàn)寬帶信號DOA估計的均方根誤差。實驗中使用1024次快拍進行協(xié)方差矩陣的估計，信號源的數(shù)目設定為10，使用的陣列結構分別為：(1)均勻線性陣列結構，使用12根天線組成一個均勻線性陣列結構，其陣列自由度是11；(2)擴展互質陣列結構，使用12根天線組成一個擴展的互質陣列，因為總數(shù)目是2M+N-1，這里設定M＝4,N＝5，其陣列自由度是2MN+1＝41；(3)本發(fā)明中的大孔徑均勻線性虛擬陣列(即對形成傳統(tǒng)互質陣列結構的各個天線的位置進行二次差分操作后形成的陣列)，所述的互質陣列結構為使用12根天線形成的傳統(tǒng)互質陣列，因為天線總數(shù)目是M+N-1，這里設定M＝6,N＝7，大孔徑均勻線性虛擬陣列的自由度為2MN+2M+2N-1＝109，經過空間平滑使用的自由度為MN+M+N＝55。比較三種方法的DOA估計均方根誤差隨信噪比變化的情況，其實驗結果如圖9所示，從圖9中可以看出，由于在陣元數(shù)目相同的情況下，使用本發(fā)明的方法可以獲得更大的自由度，檢測相同信號數(shù)目時，可更多的增加陣列可識別的信源數(shù)目，而且具有更好的檢測性能(隨著信噪比的變化表現(xiàn)出了更低的均方根誤差特性，特別是在低SNR條件下，這種優(yōu)勢更加明顯)。另外，比較三種方法的DOA估計均方根誤差隨快拍數(shù)變化情況，其實驗結果如圖10所示，從圖10中可以看出，在相同的陣元數(shù)目和相同的快拍數(shù)下，使用本發(fā)明的均方根誤差明顯低于其他方法，同時，由于本發(fā)明的方法可以獲得更大的陣列自由度，因而在DOA估計時，具有更好的譜分辨率，且系統(tǒng)成本也更低(因為沒有進行互質陣列陣元數(shù)目的擴展，而僅僅使用互質陣列的數(shù)目是M+N-1)。

附圖說明

圖1是本發(fā)明的方法流程圖；

圖2是傳統(tǒng)的互質陣列結構圖；

圖3是傳統(tǒng)的互質陣列差分集合陣列圖；

圖4是擴展的互質陣列結構圖；

圖5是擴展的互質陣列差分集合陣列圖；

圖6是本發(fā)明使用傳統(tǒng)的互質陣列結構二次差分集合陣列圖；

圖7是本發(fā)明方法的寬帶信號DOA估計譜；

圖8是本發(fā)明方法的均方根誤差隨角度間隔變化圖；

圖9是本發(fā)明方法與其他兩種陣列結構的均方根誤差隨信噪比變化圖；

圖10是本發(fā)明方法與其他兩種陣列結構的均方根誤差隨快拍數(shù)變化圖。

下面結合附圖和具體實施方式對本發(fā)明作進一步的說明。

具體實施方式

本發(fā)明的實施例：基于互質陣列的寬帶信號DOA估計方法，如圖1所示，包括以下步驟：

S1，利用天線設計互質陣列結構，每個天線稱為一個陣元；可選的，所述的互質陣列為廣義互質陣列或傳統(tǒng)互質陣列；具體的，若為傳統(tǒng)互質陣列時，則可通過以下方法設計而成：

(1)用N個天線形成均勻線性陣列1，簡稱線陣1；用M-1個天線形成均勻線性陣列2，簡稱線陣2；其中，線陣1的陣元間距為Md，線陣2的陣元間距為Nd，N＞M≥2且M與N互質，0＜d≤λ/2，λ為入射到互質陣列的寬帶信號的波長；

(2)組合線陣1和線陣2為互質陣列：設線陣1的第一個陣元為互質陣列的陣元0，將線陣2的第一個陣元放置于與陣元0相距為Nd的位置，依次進行處理，使得線陣2的所有陣元依次插于線陣1中；從互質陣列的陣元0開始，從頭至尾依次命名各個陣元為陣元0，陣元1，……，陣元M+N-1；

對形成互質陣列的各個天線的位置進行二次差分操作，形成大孔徑均勻線性虛擬陣列；

S2，對互質陣列中的天線接收到的寬帶信號進行抽樣及離散傅里葉變換，獲得頻域信號輸出模型X[l]；

S3，計算頻域信號輸出模型X[l]的自相關矩陣R_xx[l]并利用矢量化vec(·)函數(shù)對該自相關矩陣進行按列堆棧，得到一個新的向量Z[l]；再將Z[l]的協(xié)方差矩陣進行矢量化，獲得二次差分互質陣列的信號模型Z₁[l]；

S4，利用空間平滑和子陣選取技術對所述二次差分互質陣列的信號模型Z₁[l]進行去冗余、去相關處理，并獲得寬帶信號的空間平滑協(xié)方差矩陣R′，具體包括以下步驟：

S41，從新的信號模型Z₁[l]的流型矩陣B_l^*e B_l中建立一個(2MN+2M+2N-1)×K的新的流型矩陣B_1l，也即從流型矩陣B_l^*e B_l中取出不同行對應的2MN+2M+2N-1個連續(xù)的差分值，并對取出的連續(xù)的差分值進行排序；(這相當于從新信號模型Z₁[l]中移除了相應的重復行并進行排序得到一個新的向量Z₂[l]，新形成的流型矩陣B_1l等于一個帶有2MN+2M+2N-1根天線的虛擬均勻線性陣列的流型矩陣)其中，M表示形成互質陣列的 子陣1的天線數(shù)目，N表示形成互質陣列的子陣2的天線數(shù)目，K表示入射到互質陣列的信號數(shù)；

S42，將新的流型矩陣A_l中的2MN+2M+2N-1根天線劃分成NM+M+N個子陣列，劃分后的每個子陣列中有NM+M+N根天線；其中，第i個子陣列天線的位置位于{(-i+1+n)d,n＝0,1,...,NM+M+N-1}；

S43，設第i個子陣列對應新向量Z₂[l]的第(NM+M+N+1-i)行到第(2MN+2M+2N+1-i)行，從而獲得流型矩陣B_1l中每個子陣列對應的新信號模型Z_2i[l]，計算新信號模型Z_2i[l]的協(xié)方差矩陣，進而獲得寬帶信號的空間平滑協(xié)方差矩陣R′；

S5劃分空域網絡，構造字典，并利用寬帶信號的多個頻率點處的字典對所述的空間平滑協(xié)方差矩陣R′進行稀疏表示，形成寬帶信號多個字典的多測量矢量稀疏表示模型(即MMV模型)；具體包括以下步驟：

首先，利用寬帶信號的頻率w_l處的字典對所述的空間平滑協(xié)方差矩陣R′的第q列進行稀疏表示，得單個快拍矢量(SMV矢量)模型的稀疏表示為：

式中，Θ＝{θ₁,θ₂,...,θ_Q}集合表示整個DOA空間域的搜索網格，共有Q(Q>>K，K表示入射到互質陣列的信號數(shù))個柵格點；B_1l(Θ)表示整個空間內單個頻率點w_l處的流型矩陣，稱為過完備基或字典；e_q是一個(NM+M+N)×1維的列向量，它表示噪聲，它的第q行元素是1，其余位置元素均是0；b_lq是Q×1維列稀疏向量，表示字典B_1l(Θ)的稀疏表示系數(shù)，為單個列稀疏向量；為噪聲功率；(若稀疏表示系數(shù)b_lq的第q行為非零元素，則意味著對應θ_q＝θ_k(k＝1,2,...,K)處存在一個信號源，若θ_q不在柵格點上，則θ_q≈θ_k)

其次，將NM+M+N個快拍矢量(SMV矢量)模型合成一個矩陣，即得整個空間平滑協(xié)方差矩陣R′的稀疏表示，形成寬帶信號多個字典的多測量矢量稀疏表示模型(即MMV模型)：

(這里只是寫單個頻率點的稀疏表示，因為l的值是變化的，實際上是利用寬帶信號的多個頻率點處的字典對所述的空間平滑協(xié)方差矩陣R′進行稀疏表示)

其中，G_l＝[b_l1,b_l2,...,b_l(NM+M+N)]，由于G_l中的每列稀疏向量都具有共同的稀疏結構，因此G_l具有聯(lián)合稀疏性，I₁表示一個(NM+M+N)×(NM+M+N)的單位矩陣；

S6，通過對多字典稀疏表示系數(shù)的聯(lián)合稀疏約束以求解稀疏反問題的形式實現(xiàn)寬帶信號的波達方向估計；具體包括以下步驟：

S61，將寬帶信號DOA參數(shù)估計問題表述為單個頻率點w_l處的多個單字典稀疏約束優(yōu)化模型(求解此稀疏問題相當于求解L個單字典稀疏優(yōu)化模型進行聯(lián)合稀疏處理)：

式中，是G_l的估計；μ(G_l)表示矩陣G_l的非零行對應的索引集合，|μ(G_l)|＝K表示G_l中非零行的個數(shù)(即稀疏度)，K表示入射到互質陣列的信號數(shù)；

當多個字典存在時，利用所述的稀疏約束優(yōu)化模型構造一個多字典聯(lián)合優(yōu)化問題：

式中，為無混疊的約束條件；

S62，根據所述的稀疏約束優(yōu)化模型，將其轉化為凸優(yōu)化問題進行求解，從而獲得寬帶信號的DOA估計；具體包括以下步驟：(所述的多字典聯(lián)合優(yōu)化問題相當于組合優(yōu)化問題，其為NP難問題)將所述的多字典聯(lián)合優(yōu)化問題轉化為1范數(shù)凸優(yōu)化BP問題進行求解：

首先，(約束條件相當于同一稀疏結構約束，為了對其進行聯(lián)合稀疏約束)設待估計變量的行之間具有稀疏性，而列之間不具有稀疏性，同時假設bο中的每一個元素為矩陣對應該行向量的2范數(shù)，即 而||bο||₁表示只對在空域θ上進行稀疏約束；因此，對應的BP優(yōu)化模型為：

其中，||b^ο||₁表示對所有字典的表示系數(shù)的聯(lián)合行稀疏約束(這是根據壓縮感知信號恢復理論而得，由前面知，待估計變量同時b^ο中的每一個元素為矩陣對應該行向量的2范數(shù)，即這里b^ο中包含所有字典的稀疏表示系數(shù))；

其次，將所述的BP優(yōu)化模型用對應的正則化模型表示為：

式中，δ為正則化參數(shù)，表示Frobenius范數(shù)；如果能求得或(繼續(xù)轉化為下面標準的SOCP形式，根據CVX求解)，則通過其非零行對應的位置就可估計參數(shù)θ；

最后，所述的正則化模型為凸優(yōu)化問題，對其進行求解，獲得寬帶信號的DOA估計；

可選的，可將所述的正則化模型轉化為標準的SOCP形式：

式中，1_Q是一個所有元素為1的Q×1維列向量；z和η為引入的輔助變量，γ是一個1×Q的向量，其中第q個元素是γ_q，目標函數(shù)為線性函數(shù)；

采用數(shù)值軟件包sedumi或者CVX對所述的標準的SOCP形式進行求解，進而獲得寬帶信號的DOA估計。

實驗例：基于互質陣列的寬帶信號DOA估計方法，包括以下步驟：

S1，利用天線設計互質陣列結構，每個天線稱為一個陣元；具體的，所述的互質陣列可通過以下方法設計而成(如圖2所示)：

(1)用N個天線形成均勻線性陣列1，簡稱線陣1；用M-1個天線形成均勻線性陣列2，簡稱線陣2；其中，線陣1的陣元間距為Md，線陣2的陣元間距為Nd，N＞M≥2且M與N互質，0＜d≤λ/2，λ為入射到互質陣列的寬帶信號的波長；

(2)組合線陣1和線陣2為互質陣列：設線陣1的第一個陣元為互質陣列的陣元0，將線陣2的第一個陣元放置于與陣元0相距為Nd的位置，依次進行處理，使得線陣2的所有陣元依次插于線陣1中；從互質陣列的陣元0開始，從頭至尾依次命名各個陣元為陣元0，陣元1，……，陣元M+N-1；傳統(tǒng)的互質天線陣列位置用集合表示為：S＝{Mnd,0≤n≤N-1}U{Nmd,1≤m≤M-1}；

假設M＝3,N＝4，對該互質陣列進行一次差分操作，如圖3所示，差分集合陣列擴展形成虛擬陣列后存在空穴，即位置7處，因而阻礙了其應用到許多實際的應用中，包括使用空間平滑技術的一些DOA估計情況；

如圖4所示為一種擴展互質陣列，所述的擴展互質陣列與傳統(tǒng)的互質陣列有相同的結構，擴展互質陣列天線位置集合表示為：S_E＝{Mnd,0≤n≤N-1}U{Nmd,0≤m≤2M-1}， 只是相比傳統(tǒng)的互質陣列，其第二個子陣列的天線數(shù)目加倍，即多出M根天線，所以擴展互質陣列總的天線數(shù)目是2M+N-1。圖5為由擴展互質陣列進行差分集合陣列擴展形成的虛擬陣列，同樣假設M＝3,N＝4。由圖5可以看出，擴展互質陣列進行差分集合陣列擴展后可以形成一個更大的連續(xù)的虛擬線性陣列，連續(xù)的天線位置至少位于C＝{kd,-MN≤k≤MN}，這種情況下，使用2M+N-1個天線就可以獲得O(MN)自由度，和傳統(tǒng)互質陣列的差分集合陣列擴展相比，盡管擴展的互質陣列形成更大的虛擬陣列是連續(xù)的，且增加了陣列自由度、增大了陣列孔徑，但是它需要更多的天線。

本發(fā)明采用的是：對形成互質陣列的各個天線的位置進行二次差分操作，形成大孔徑均勻線性虛擬陣列；其中，所述的互質陣列和傳統(tǒng)的互質陣列結構相同，如圖2所示，由M+N-1個天線組成傳統(tǒng)的互質陣列。圖6為一個傳統(tǒng)互質陣列結構進行二次差分集合陣列擴展后形成的虛擬均勻線性陣列，其中，M＝3,N＝4，圖中只給出了二次差分集合陣列擴展后的正數(shù)部分。由此可以看出：傳統(tǒng)的互質陣列結構經過二次差分后形成了一個更大的連續(xù)虛擬均勻線性陣列；而且和擴展互質陣列技術相比，省去了M根天線，節(jié)約了成本。實際上，本發(fā)明使用互質陣列二次差分技術擴展后形成連續(xù)的虛擬均勻線性陣列的天線位置范圍是(-MN-M-N+1)d到(MN+M+N-1)d，相比傳統(tǒng)的互質陣列結構和擴展的互質陣列結構經過一次差分集合陣列擴展后，明顯的獲得了一個自由度的提高。

S2，假設有K個寬帶信號源s_k(t),k＝1,2,...,K位于遠場入射到互質天線陣列，其來波方向為θ＝[θ₁，θ₂，...，θ_K]^T；線陣1和線陣2接收信號表達式如下：

其中，τ_1,n(θ_k)，τ_2,m(θ_k)分別表示角度為θ_k的信號到達線陣1和線陣2的時間延遲， 分別表示線陣1和線陣2上天線的觀測噪聲；

用一個頻率f_s對其進行抽樣，線陣1和線陣2接收信號的離散形式表示為：

x₁[i]＝[x_1,0[i],x_1,1[i],...,x_1,N-1[i]]^T

x₂[i]＝[x_2,1[i],x_2,2[i],...,x_1,M-1[i]]^T

其中，x₁[i]表示線陣1中各個天線接收信號的離散形式，x₂[i]表示線陣2中接收信號的離散形式。

對每個天線接收的信號采用L點DFT變換，相鄰的DFTs沒有重疊，線陣1和線陣2的第l個 頻帶樣本可以寫成一個向量形式，即

x₁[l]＝[x_1,0[l],x_1,1[l],...,x_1,N-1[l]]^T

x₂[l]＝[x_2,1[l],x_2,2[l],...,x_1,M-1[l]]^T

其中，x₁[l]表示線陣1上第l個頻帶的接收信號表示形式，x₂[l]表示線陣2上第l個頻帶的接收信號表示形式。

式中，x₁[i]表示線陣1中各個天線接收信號的離散形式，L表示DFT點數(shù)，也是總的頻帶數(shù)目，l表示第l個頻帶；

經過DFT變換后，輸出信號頻域模型可以表示為：

其中，A₁(l,θ)＝[a₁(l,θ₁),a₁(l,θ₂),...,a₁(l,θ_K)]和A₂(l,θ)＝[a₂(l,θ₁),a₂(l,θ₂),...,a₂(l,θ_K)]是線陣1和線陣2對應的第l個頻帶在頻率f_l處的指向矩陣，S[l]＝[S₁[l]，S₂[l]，L,S_K[l]]^T表示所有信號在第l個頻帶的表示形式，和表示線陣1和線陣2對應的噪聲向量。將上述線陣1和線陣2頻域輸出信號模型進行合并，有如下定義

S3，計算頻域信號輸出模型X[l]的自相關矩陣R_xx[l]：

式中，其中，表示第l頻帶第k個信號的功率，I表示一個單位矩陣，為噪聲功率，每個頻帶上噪聲均為

然后利用矢量化函數(shù)vec(·)矢量化上述的自相關矩陣R_xx[l]，得到下面的向量：

式中，B_l＝A^*(l,θ)e A(l,θ)，符號e表示Khatri-Rao(KR)積，與陣列觀測向量X[l]相比，可以說Z[l]也類似一個天線陣列的接收信號，其流型矩陣通過B_l給出，等源信號向量表示為p，噪聲變成確定性的向量

為了獲得二次差分陣列的接收數(shù)據模型，進行如下操作：

定義

由于p和均為確定性向量，則有

讓pp^H＝B＝R_P+B,其中，B中的元素遠小于R_P中的元素，R_P是一個對角矩陣，對角元素表示信號功率的自相關，B中的元素表示信號功率的互相關，而且diag(B)＝diag(R_P)成立。所以,

利用克羅內克積的性質：

矢量化R_z[l]，得

通過簡化得

Z₁[l]＝B_l^*e B_lq+e

其中，q＝diag(R_P)，表示噪聲，這里已經不是白噪聲。注意到Z₁[l]可以認為是表示另一個陣列信號模型，其流型矩陣為B_l^*e B_l，源信號矢量是q，類似于相關源，噪聲為e。如果原始陣列天線位置可以通過集合{x_i,1≤i≤M+N-1},那么經過二次差分后形成的虛擬天線位置表示為集合形式{x_i+x_j-x_k-x_l,1≤i,j,k,l≤M+N-1}。正如上面所述，利用傳統(tǒng)的互質陣列天線位置經過二次差分集合陣列擴展后形成連續(xù)的差分值一共包括2NM+2M+2N-1個，其范圍是從-NM-M-N+1到NM+M+N-1。所以對互質陣列進行二次差分技術處理進行DOA估計明顯地給出了2NM+2M+2N-1個自由度，而原始的陣列天線結構僅僅提供M+N-1個自由度。

為了簡化計算，我們移除噪聲e中的兩部分，即信號功率的互相關部分 信號功率與噪聲功率的互相關部分，即則形成的虛擬陣列模型為

其中，e₁＝vec(11^H)；

S4，利用空間平滑和子陣選取技術對所述的二次差分互質陣列的信號模型Z₁[l]進行去冗余、去相關處理，并獲得寬帶信號的空間平滑協(xié)方差矩陣R′，具體包括以下步驟：

S41，新流型矩陣B′對應一個更多陣元數(shù)目的均勻線性陣列形成的流型矩陣，其陣列天線位置在所有間距d的整數(shù)倍位置，范圍從(-NM-M-N+1)d到(NM+M+N-1)d，等源信號向量q類似于相關源。因此，使用子陣選取和空間平滑技術應用于新信號模型Z₁[l]，同時增加從這個信號模型獲得的半定矩陣的秩。為了采用空間平滑，首先從新流型矩陣B_l^*e B_l建立一個(2NM+2M+2N-1)×K的矩陣B_1l，也就是精確地取出對應的2NM+2M+2N-1個連續(xù)的差分值位置對應的行，將取出的不同差分值進行排序；這相當于從新觀測向量Z₁[l]中移除了重復的行并進行排序得到一個新的向量表示如下：

式中，e₂是一個(2NM+2M+2N-1)×1列向量，表示噪聲。

S42，將新的流型矩陣B_1l中的2NM+2M+2N-1根天線劃分成NM+M+N個子陣列，劃分后的每個子陣列中有NM+M+N根天線；其中，第i個子陣列天線的位置位于{(-i+1+n)d,n＝0,1,...,NM+M+N-1}；

S43，設第i個子陣列的輸出信號流型對應新向量Z₂[l]的第(NM+M+N+1-i)到第(2NM+2M+2N+1-i)行，把它表示為：

Z_2i[l]＝B_1liq+e_2i

式中，B_1li表示第i個子陣列對應的流型矩陣，它是一個由新形成的流型矩陣B_1l的第(NM+M+N+1-i)到第(2NM+2M+2N+1-i)行組成的一個(NM+M+N)×K維矩陣，e_2i表示第i個子陣列對應的噪聲。將上式Z_2i[l]表達為另外一種形式

Z_2i[l]＝B_1lΦ^i-1q+e_2i

式中，

而且其中

計算Z_2i[l]求其協(xié)方差矩陣，具體如下所述：

對于所有i取R_2i[l]的平均值，得到：

又矩陣R_s可以表示為其中，

另外可以把上式簡寫成

式中，I₁是一個(NM+M+N)×(NM+M+N)的單位矩陣，把矩陣R′稱為寬帶信號的空間平滑協(xié)方差矩陣，它與一個由NM+M+N個天線組成的均勻線性陣列的輸出信號模型的協(xié)方差矩陣具有相同的表達形式，由此可知，可以使得僅僅使用M+N-1個天線組成的陣列結構執(zhí)行NM+M+N-1個寬帶源的DOA估計。

S5，劃分空域網絡，構造字典，并利用寬帶信號的多個頻率點處的字典對所述的空間平滑協(xié)方差矩陣R′進行稀疏表示，形成寬帶信號多個字典的多測量矢量稀疏表示模型(即MMV模型)；具體包括以下步驟：

假設集合Θ＝{θ₁,θ₂,...,θ_Q}表示整個DOA空間域的搜索網格，共有Q(Q＞＞K)個柵格點，于是整個空間內頻率w_l處的流型矩陣可以表示為B_1l(Θ)。因為Q＞＞K，所以把稱為過完備基或者字典。那么寬帶信號的輸出模型經過空間平滑后求得的協(xié)方差矩陣R′的第q列是：

式中，是字典的稀疏表示系數(shù)。的第q行元素是1，其余位置元素均是0。如果稀疏表示系數(shù)b_lq的第q行為非零元素，則意味著對應θ_q＝θ_k(k＝1,2,...,K)處存在一個信號源，若θ_q不在柵格點上，則θ_q≈θ_k，b_lq中非零元素的個數(shù)為K。因b_lq是單個列稀疏向量，所以稱R′_q為單個快拍矢量(SMV)模型的稀疏表示。將MN+M+N個SMV矢量合成一個矩陣，得到下述稀疏表示問題：

其中，G_l＝[b_l1,b_l2,...,b_l(NM+M+N)]，易知G_l中的每列稀疏向量都具有共同的稀疏結構，因此G_l具有聯(lián)合稀疏性，上式的合成矩陣R′稱為多測量(MMV)模型的稀疏表示；

S6，通過對多字典稀疏表示系數(shù)的聯(lián)合稀疏約束以求解稀疏反問題的形式實現(xiàn)寬帶信號的波達方向估計；具體包括以下步驟：

假設μ(G_l)表示矩陣G_l的非零行對應的索引集合，|μ(G_l)|＝K表示G_l中非零行的個數(shù)(即為稀疏度)，則對應的寬帶信號DOA參數(shù)估計問題可以表述為如下約束最優(yōu)化問題：

式中，是G_l的估計。上式的最優(yōu)化問題是單個頻率點w_l處的單字典稀疏優(yōu)化模型，求解此稀疏問題相當于求解L個單字典稀疏優(yōu)化模型進行聯(lián)合稀疏處理。

當多個字典存在時，將利用上述的稀疏約束優(yōu)化模型構造一個多字典聯(lián)合優(yōu)化問題：

多字典聯(lián)合優(yōu)化問題相當于組合優(yōu)化問題，其為NP難問題，可以將其轉化為1范數(shù)凸優(yōu)化BP問題求解。約束條件相當于同一稀疏結構約束，為了對其進行聯(lián)合稀疏約束，設待估計變量明顯地，的行之間具有稀疏性，而列之間不具有稀疏性。設b^ο中的每一個元素為矩陣對應該行向量的2范數(shù)，即 而||b^ο||₁表示只對在空域θ上進行稀疏約束。因此，對應的BP優(yōu)化模型為

其中，||b^ο||₁表示對所有字典的表示系數(shù)的聯(lián)合行稀疏約束。

其次，將所述的BP優(yōu)化模型用對應的正則化模型表示為：

式中，δ為正則化參數(shù)，表示Frobenius范數(shù)；如果能求得或(繼續(xù)轉化為下 面標準的SOCP形式，根據CVX求解)，則通過其非零行對應的位置就可估計參數(shù)θ。上面的正則化模型是凸優(yōu)化問題，可以采用二階錐規(guī)劃(SOCP)求解，為了將其轉化為標準的SOCP形式，引入輔助變量z和η使得目標函數(shù)是線性函數(shù)，SOCP標準形式是：

式中，1_Q是一個所有元素為1的Q×1維列向量；z和η為引入的輔助變量，γ是一個1×Q的向量，其中第q個元素是γ_q，目標函數(shù)為線性函數(shù)；

所述的標準的SOCP問題可用數(shù)值軟件包sedumi或CVX求解，進而可以獲得DOA估計。