本發(fā)明屬于信號處理技術(shù)領(lǐng)域，尤其涉及對雷達(dá)信號、聲學(xué)信號及電磁信號的波達(dá)方向估計(jì)，具體是一種基于內(nèi)插虛擬陣列協(xié)方差矩陣toeplitz化重建的互質(zhì)陣列波達(dá)方向估計(jì)方法，可用于無源定位和目標(biāo)探測。

背景技術(shù)：

波達(dá)方向(direction-of-arrival,doa)估計(jì)是陣列信號處理領(lǐng)域的一個重要分支，它是指利用陣列天線接收空域信號，并通過現(xiàn)代信號處理技術(shù)和各類優(yōu)化方法實(shí)現(xiàn)對接收信號統(tǒng)計(jì)量的有效處理，從而實(shí)現(xiàn)信號的doa估計(jì)，在雷達(dá)、聲吶、語音、無線通信等領(lǐng)域有著重要的應(yīng)用價(jià)值。

doa估計(jì)方法的自由度是指其能夠估計(jì)的入射信號源的個數(shù)。現(xiàn)有的doa估計(jì)方法通常采用均勻線性陣列進(jìn)行信號的接收與建模，但是基于均勻線性陣列方法的自由度受限于實(shí)際天線陣元個數(shù)。具體而言，對于一個包含l個天線陣元的均勻線性陣列，其自由度為l-1。因此，當(dāng)某個空域范圍內(nèi)入射信號源的個數(shù)大于陣列中天線陣元的個數(shù)時(shí)，現(xiàn)有采用均勻線性陣列的方法將無法進(jìn)行有效的doa估計(jì)。

互質(zhì)陣列能夠在天線陣元個數(shù)一定的前提下增加doa估計(jì)的自由度，因而受到了學(xué)術(shù)界的廣泛關(guān)注。作為互質(zhì)采樣技術(shù)在空間域上的一個典型表現(xiàn)形式，互質(zhì)陣列提供了一個系統(tǒng)化的稀疏陣列架構(gòu)方案，并能夠突破傳統(tǒng)均勻線性陣列自由度受限的瓶頸，實(shí)現(xiàn)doa估計(jì)方法自由度性能的提升?，F(xiàn)有的基于互質(zhì)陣列的doa估計(jì)方法主要通過利用質(zhì)數(shù)的性質(zhì)將互質(zhì)陣列推導(dǎo)到虛擬域，并形成等價(jià)虛擬均勻線性陣列接收信號以實(shí)現(xiàn)doa估計(jì)。由于虛擬陣列中包含的虛擬陣元數(shù)大于實(shí)際的天線陣元數(shù)，自由度因此得到了有效的提升。但是由于從互質(zhì)陣列推導(dǎo)而來的虛擬陣列屬于非均勻陣列，因此很多現(xiàn)有基于均勻線性陣列的信號處理方法無法直接應(yīng)用于虛擬陣列等價(jià)接收信號的doa估計(jì)。當(dāng)前采用互質(zhì)陣列的doa估計(jì)方法常用的一個解決方案是，僅利用虛擬陣列中連續(xù)的陣元部分形成一個虛擬均勻線陣以進(jìn)行doa估計(jì)，但是這造成了部分原始信息的丟失和相關(guān)估計(jì)性能的降低。

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

本發(fā)明的目的在于針對上述現(xiàn)有技術(shù)存在的不足，提出一種基于內(nèi)插虛擬陣列協(xié)方差矩陣toeplitz化重建的互質(zhì)陣列波達(dá)方向估計(jì)方法，充分利用了非均勻虛擬陣列所提供的全部信息，從而提高了doa估計(jì)的自由度與分辨率。

本發(fā)明的目的是通過以下技術(shù)方案來實(shí)現(xiàn)的：一種基于內(nèi)插虛擬陣列協(xié)方差矩陣toeplitz化重建的互質(zhì)陣列波達(dá)方向估計(jì)方法，包含以下步驟：

(1)接收端使用m+n-1個天線，并按照互質(zhì)陣列結(jié)構(gòu)進(jìn)行架構(gòu)；其中m與n為互質(zhì)整數(shù)；

(2)假設(shè)有k個來自θ1,θ2,…,θk方向的遠(yuǎn)場窄帶非相干信號源，則(m+n-1)×1維互質(zhì)陣列接收信號x(t)可建模為：

其中，sk(t)為信號波形，n(t)為與各信號源相互獨(dú)立的噪聲分量，a(θk)為θk方向的導(dǎo)引矢量，表示為：

其中，pid,i＝1,2,…,m+n-1表示互質(zhì)陣列中第i個物理天線陣元的實(shí)際位置，且p1＝0；d為入射窄帶信號波長λ的一半，即d＝λ/2，[·]^t表示轉(zhuǎn)置操作。共采集t個采樣快拍，得到采樣協(xié)方差矩陣

這里，(·)^h表示共軛轉(zhuǎn)置；

(3)計(jì)算互質(zhì)陣列接收信號所對應(yīng)的等價(jià)虛擬信號：矢量化互質(zhì)陣列接收信號的采樣協(xié)方差矩陣獲得虛擬陣列等價(jià)接收信號v：

其中，為(m+n-1)²×k維虛擬陣列導(dǎo)引矩陣，包含k個入射信號源的功率，為噪聲功率，iv＝vec(im+n-1)。這里，vec(·)表示矢量化操作，即把矩陣中的各列依次堆疊以形成一個新的矢量，(·)^*表示共軛操作，表示克羅內(nèi)克積，im+n-1表示(m+n-1)×(m+n-1)維單位矩陣。矢量v對應(yīng)的虛擬陣列中各虛擬陣元的位置為

去除集合中各位置上重復(fù)的虛擬陣元，得到一個非均勻的虛擬陣列其對應(yīng)的等價(jià)虛擬信號vc可通過選取矢量v中相對應(yīng)的元素獲得；

(4)構(gòu)造內(nèi)插虛擬陣列及其接收信號并建模：首先對于非均勻的虛擬陣列在保留其原有虛擬陣元位置不變的前提下，向其中非連續(xù)的位置插入若干虛擬陣元，從而將非均勻虛擬陣列轉(zhuǎn)化為間距為d、陣列孔徑與互質(zhì)陣列相同、且虛擬陣元數(shù)目增加的均勻虛擬陣列該內(nèi)插均勻虛擬陣列共包含個虛擬陣元，其中|·|表示集合的勢，其對應(yīng)的等價(jià)虛擬信號vi可通過往矢量vc中插入0獲得，插入0的位置與中插入的虛擬陣元的位置相對應(yīng)；

(5)構(gòu)造內(nèi)插虛擬陣列多采樣快拍信號及其采樣協(xié)方差矩陣：將切割為li個長度為li的連續(xù)子陣列，其中

相應(yīng)地，內(nèi)插虛擬陣列的多采樣快拍信號可通過截取矢量vi中對應(yīng)的元素獲得，即：vi,l,l＝1,2,…,li由vi中第li+1-l到第2li-l個元素組成。接著，vi的采樣協(xié)方差矩陣rv可以由如下方式得到：

其中，<vi>i表示位置為id的虛擬陣元所對應(yīng)的等價(jià)接收信號；

(6)構(gòu)造投影矩陣并定義投影運(yùn)算：投影矩陣p的維度與rv相同，如果矩陣rv中某個元素為0，則投影矩陣p中相同位置的元素值也為0；否則投影矩陣p中的元素值為1。定義為投影運(yùn)算，其中括號內(nèi)變量為與p維度相同的矩陣，投影運(yùn)算通過變量矩陣中每一個元素與投影矩陣p中相應(yīng)位置上的元素一一相乘實(shí)現(xiàn)，得到一個與矩陣p維度相同的矩陣；

(7)設(shè)計(jì)基于內(nèi)插虛擬陣列協(xié)方差矩陣toeplitz化重建的優(yōu)化問題并求解：根據(jù)接收信號理論協(xié)方差矩陣的toeplitz性，利用(5)中得到的內(nèi)插虛擬陣列協(xié)方差矩陣rv作為參考值，尋找一個與其差異最小的低秩toeplitz矩陣作為接收信號的協(xié)方差矩陣，可構(gòu)建如下以矢量z為變量的優(yōu)化問題：

其中，是以矢量z為第一列的厄米特對稱toeplitz矩陣；∈為閾值常數(shù)，用于約束協(xié)方差矩陣的重建誤差；保證了重建的協(xié)方差矩陣滿足半正定的條件；‖·‖f表示frobenius范數(shù)；rank(·)表示矩陣的秩。將上述非凸優(yōu)化問題轉(zhuǎn)化為凸優(yōu)化問題，并求得最優(yōu)化值相應(yīng)地，重建的toeplitz矩陣為內(nèi)插虛擬陣列協(xié)方差矩陣；

(8)根據(jù)重建的內(nèi)插虛擬陣列協(xié)方差矩陣進(jìn)行波達(dá)方向估計(jì)。

進(jìn)一步地，步驟(1)所述的互質(zhì)陣列結(jié)構(gòu)可具體描述為：首先選取一對互質(zhì)整數(shù)m、n；然后，構(gòu)造一對稀疏均勻線性子陣列，其中第一個子陣列包含m個間距為nd的天線陣元，其位置為0,nd,…,(m-1)nd，第二個子陣列包含n個間距為md的天線陣元，其位置為0,md,…,(n-1)md；接著，將兩個子陣列按照首個陣元重疊的方式進(jìn)行子陣列組合，獲得實(shí)際包含m+n-1個天線陣元的非均勻互質(zhì)陣列架構(gòu)。

進(jìn)一步地，步驟(7)中所構(gòu)建的非凸優(yōu)化問題可通過凸松弛技術(shù)，將優(yōu)化問題目標(biāo)函數(shù)中的矩陣的秩最小化操作替換為矩陣的跡最小化操作，得到以下以矢量z為變量的凸優(yōu)化問題：

其中，tr(·)表示矩陣的跡。

進(jìn)一步地，步驟(7)中所構(gòu)建的非凸優(yōu)化問題可轉(zhuǎn)化為如下以矢量z為變量的凸優(yōu)化問題：

其中μ為正則化參數(shù)，用于在最小化過程中權(quán)衡矩陣重建誤差和矩陣的跡。

進(jìn)一步地，步驟(8)中的波達(dá)方向估計(jì)，可采用以下方法：多重信號分類方法、旋轉(zhuǎn)不變子空間方法、求根多重信號分類方法、協(xié)方差矩陣稀疏重建方法等。

進(jìn)一步地，步驟8中，通過多重信號分類方法進(jìn)行波達(dá)方向估計(jì)，具體為：畫出虛擬域空間譜pmusic(θ)：

其中d(θ)是li×1維內(nèi)插虛擬陣列導(dǎo)引矢量，對應(yīng)于位置為由0到(li-1)d的一段虛擬均勻陣列；en是li×(li-k)維矩陣，表示內(nèi)插虛擬陣列協(xié)方差矩陣的噪聲子空間；θ是假定的信號波達(dá)方向；通過譜峰搜索尋找空間功率譜pmusic(θ)上的峰值，并將這些峰值所對應(yīng)的響應(yīng)值從大到小排列，取前k個峰值所對應(yīng)的角度方向，即為波達(dá)方向估計(jì)結(jié)果。

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

(1)本發(fā)明在互質(zhì)陣列等價(jià)虛擬域上引入了陣列內(nèi)插的思想，充分利用了虛擬陣列提供的全部信息。通過在非均勻虛擬陣列中內(nèi)插虛擬陣元的方式構(gòu)建出均勻線性陣列，保留了由原始非均勻虛擬陣列接收到的全部信息，提高了doa估計(jì)的自由度及分辨率；

(2)本發(fā)明基于內(nèi)插虛擬陣列協(xié)方差矩陣toeplitz化重建的思想設(shè)計(jì)優(yōu)化問題。由于均勻線性陣列的理論協(xié)方差矩陣滿足toeplitz結(jié)構(gòu)，因此利用它的toeplitz特性進(jìn)行協(xié)方差矩陣的重建可以使得重建結(jié)果與真實(shí)值差異更小，從而提高doa估計(jì)方法的性能。

附圖說明

圖1是本發(fā)明的方法總體流程框圖。

圖2是本發(fā)明中組成互質(zhì)陣列的一對稀疏均勻子陣列結(jié)構(gòu)示意圖。

圖3是本發(fā)明中互質(zhì)陣列的結(jié)構(gòu)示意圖。

圖4是本發(fā)明中內(nèi)插虛擬陣列的結(jié)構(gòu)示意圖。

圖5是本發(fā)明中內(nèi)插虛擬陣列分割方法的示意圖。

圖6是用于體現(xiàn)本發(fā)明所提方法自由度性能的空間功率譜示意圖。

圖7是用于體現(xiàn)本發(fā)明所提方法分辨率性能的歸一化空間譜示意圖。

具體實(shí)施方式

以下參照附圖，對本發(fā)明的技術(shù)方案和效果作進(jìn)一步的詳細(xì)說明。

對于doa估計(jì)方法在實(shí)際系統(tǒng)中的應(yīng)用，互質(zhì)陣列由于其可以通過等價(jià)虛擬陣列信號的計(jì)算和統(tǒng)計(jì)信號處理，突破物理陣元數(shù)量對自由度的限制而備受青睞。但是受限于虛擬陣列的非均勻性，目前很多方法都會選擇利用其中連續(xù)的部分進(jìn)行doa估計(jì)，從而造成了信息損失。為了充分利用非均勻虛擬陣列的所有信息，本發(fā)明提供了一種基于內(nèi)插虛擬陣列協(xié)方差矩陣toeplitz化重建的互質(zhì)陣列波達(dá)方向估計(jì)方法，參照圖1，本發(fā)明的實(shí)現(xiàn)步驟如下：

步驟一：在接收端使用m+n-1個天線陣元架構(gòu)互質(zhì)陣列；首先，選取一組互質(zhì)整數(shù)m、n；然后，參照圖2，構(gòu)造一對稀疏均勻線性子陣列，其中第一個子陣列包含m個間距為nd的天線陣元，其位置為0,nd,…,(m-1)nd；第二個子陣列包含n個間距為md的天線陣元，其位置為0,md,…,(n-1)md；單位間距d取為入射窄帶信號波長λ的一半，即d＝λ/2；接著，將兩個子陣列的首個天線陣元視為參考陣元，參照圖3，將兩個子陣的參考陣元重疊以實(shí)現(xiàn)子陣列組合，獲得實(shí)際包含m+n-1個天線陣元的非均勻互質(zhì)陣列架構(gòu)。

步驟二：采用互質(zhì)陣列接收信號并建模。假設(shè)有k個來自θ1,θ2,…,θk方向的遠(yuǎn)場窄帶非相干信號源，采用步驟一架構(gòu)的非均勻互質(zhì)陣列接收入射信號，得到(m+n-1)×1維互質(zhì)陣列接收信號x(t)，可建模為：

其中，sk(t)為信號波形，n(t)為與各信號源相互獨(dú)立的噪聲分量，a(θk)為θk方向的互質(zhì)陣列導(dǎo)引矢量，表示為

其中，pid,i＝1,2,…,m+n-1表示互質(zhì)陣列中第i個物理天線陣元的實(shí)際位置，且p1＝0；[·]^t表示轉(zhuǎn)置操作。共采集t個采樣快拍，得到采樣協(xié)方差矩陣

其中，(·)^h表示共軛轉(zhuǎn)置。

步驟三：計(jì)算互質(zhì)陣列接收信號所對應(yīng)的等價(jià)虛擬信號。矢量化互質(zhì)陣列接收信號的采樣協(xié)方差矩陣獲得虛擬陣列等價(jià)接收信號v：

其中，為(m+n-1)²×k維虛擬陣列導(dǎo)引矩陣，包含k個入射信號源的功率，為噪聲功率，iv＝vec(im+n-1)。這里，vec(·)表示矢量化操作，即把矩陣中的各列依次堆疊以形成一個新的矢量，(·)^*表示共軛操作，表示克羅內(nèi)克積，im+n-1表示(m+n-1)×(m+n-1)維單位矩陣。矢量v對應(yīng)的虛擬陣列中各虛擬陣元的位置為其中

去除集合中各位置上重復(fù)的虛擬陣元，得到一個非均勻的虛擬陣列其對應(yīng)的等價(jià)虛擬信號vc可通過選取矢量v中相對應(yīng)的元素獲得。

步驟四：構(gòu)造內(nèi)插虛擬陣列及其接收信號建模。參照圖4，對于非均勻的虛擬陣列在保留其原有虛擬陣元位置不變的前提下，向其中存在孔洞的位置插入若干虛擬陣元(如圖4中的空心圓所示)，從而將非均勻虛擬陣列轉(zhuǎn)化為間距為d、陣列孔徑與互質(zhì)陣列相同、且虛擬陣元數(shù)目增加的均勻虛擬陣列內(nèi)插虛擬陣列共包含個虛擬陣元，其中|·|表示集合的勢。內(nèi)插虛擬陣列對應(yīng)的等價(jià)虛擬信號vi可通過往矢量vc中孔洞的相應(yīng)位置填入0獲得。

步驟五：構(gòu)造內(nèi)插虛擬陣列多采樣快拍信號及其采樣協(xié)方差矩陣。參照圖5，將切割為li個長度為li的連續(xù)子陣列，其中

由于中的虛擬陣元以零位對稱分布，始終為奇數(shù)，故li為實(shí)數(shù)。相應(yīng)地，內(nèi)插虛擬陣列的多采樣快拍信號可通過截取矢量vi中對應(yīng)的元素獲得，即：其中vi,l,l＝1,2,…,li由vi中第li+1-l到第2li-l個元素組成。接著，vi的采樣協(xié)方差矩陣rv可以由如下方式得到：

其中，<vi>i表示位置為id的虛擬陣元所對應(yīng)的等價(jià)接收信號。

步驟六：構(gòu)造投影矩陣并定義投影運(yùn)算。由于步驟五所得的協(xié)方差矩陣rv中包含有在步驟四中插入的0，因此其相應(yīng)位置的對角線上元素全部為0。根據(jù)這樣的結(jié)構(gòu)定義一個與rv維度相同的投影矩陣p，如果rv中某一位置上的元素是0，則投影矩陣p中相同位置的元素值也為0；否則投影矩陣p中的元素值為1。定義為投影運(yùn)算，其中括號內(nèi)變量為與p維度相同的矩陣，投影運(yùn)算通過變量矩陣中每一個元素與投影矩陣p中相應(yīng)位置上的元素一一相乘實(shí)現(xiàn)，得到一個與矩陣p維度相同的矩陣。

步驟七：設(shè)計(jì)基于內(nèi)插虛擬陣列協(xié)方差矩陣toeplitz化重建的優(yōu)化問題并求解。根據(jù)接收信號理論協(xié)方差矩陣的toeplitz性，利用步驟五得到的內(nèi)插虛擬陣列協(xié)方差矩陣rv作為參考值，尋找一個與其差異最小的低秩toeplitz矩陣作為接收信號的協(xié)方差矩陣，可構(gòu)建如下以矢量z為變量的優(yōu)化問題：

其中，表示以矢量z為第一列的厄米特對稱toeplitz矩陣；為閾值常數(shù)，用于約束協(xié)方差矩陣的重建誤差；保證了重建的協(xié)方差矩陣滿足半正定的條件；‖·‖f表示frobenius范數(shù)，rank(·)表示矩陣的秩。求解上述非凸優(yōu)化問題可得到最優(yōu)化值由于上述優(yōu)化問題包含求解矩陣的秩最小化這一非凸項(xiàng)，這將導(dǎo)致求解困難；為了獲得優(yōu)化解，我們考慮引入凸松弛技術(shù)，將上述優(yōu)化問題目標(biāo)函數(shù)中矩陣的秩最小化操作替換為矩陣的跡最小化操作，得到以下以矢量z為變量的凸優(yōu)化問題：

其中tr(·)表示矩陣的跡。上述凸優(yōu)化問題可等價(jià)寫為以下以矢量z為變量的優(yōu)化問題：

其中μ為正則化參數(shù)，用于在最小化過程中權(quán)衡矩陣重建誤差和矩陣的跡。求解上述凸優(yōu)化問題可得到最優(yōu)化值相應(yīng)地，重建的toeplitz矩陣為內(nèi)插虛擬陣列協(xié)方差矩陣。

步驟八：根據(jù)重建的內(nèi)插虛擬陣列協(xié)方差矩陣進(jìn)行波達(dá)方向估計(jì)。通過引入經(jīng)典的方法，如多重信號分類方法、旋轉(zhuǎn)不變子空間方法、求根多重信號分類方法、協(xié)方差矩陣稀疏重建方法等，可以求得波達(dá)方向估計(jì)結(jié)果。以多重信號分類方法為例，畫出虛擬域空間譜pmusic(θ)

其中d(θ)是li×1維內(nèi)插虛擬陣列導(dǎo)引矢量，對應(yīng)于位置為由0到(li-1)d的一段虛擬均勻陣列；en是li×(li-k)維矩陣，表示內(nèi)插虛擬陣列協(xié)方差矩陣的噪聲子空間；θ是假定的信號波達(dá)方向；通過譜峰搜索尋找空間功率譜pmusic(θ)上的峰值，并將這些峰值所對應(yīng)的響應(yīng)值從大到小排列，取前k個峰值所對應(yīng)的角度方向，即為波達(dá)方向估計(jì)結(jié)果。

本發(fā)明一方面引入虛擬陣列內(nèi)插的思想，在推導(dǎo)的原始虛擬陣列基礎(chǔ)上內(nèi)插入虛擬陣元，從而將原始的非均勻虛擬陣列轉(zhuǎn)化為虛擬均勻陣列，同時(shí)保留了原始非均勻虛擬陣列上的所有信息，避免了因原始虛擬陣列的非均勻性所導(dǎo)致的統(tǒng)計(jì)信號處理模型失配及傳統(tǒng)方法截取虛擬均勻子陣所導(dǎo)致的信息損失問題；另一方面，引入了基于toeplitz特性的內(nèi)插虛擬陣列協(xié)方差矩陣重建思想，并將其應(yīng)用于虛擬域以實(shí)現(xiàn)doa估計(jì)。

下面結(jié)合仿真實(shí)例對本發(fā)明的效果做進(jìn)一步的描述。

仿真實(shí)例1：采用互質(zhì)陣列接收入射信號，其參數(shù)選取為m＝3，n＝5，即架構(gòu)的互質(zhì)陣列共包含m+n-1＝7個物理陣元。假定入射窄帶信號個數(shù)為9，且入射方向均勻分布于-50°至50°這一空間角度域范圍內(nèi)；信噪比設(shè)置為30db，采樣快拍數(shù)t＝500；正則化參數(shù)μ設(shè)置為2.5×10^-3/（（logt)²log（m+n-1））。

本發(fā)明所提出的基于內(nèi)插虛擬陣列協(xié)方差矩陣toeplitz化重建的互質(zhì)陣列波達(dá)方向估計(jì)方法空間功率譜如圖6所示，其中垂直虛線代表入射信號源的實(shí)際方向?？梢钥闯?，本發(fā)明所提方法能夠有效分辨這9個入射信號源。而對于傳統(tǒng)采用均勻線性陣列的方法，利用7個物理天線陣元最多只能分辨6個入射信號，以上結(jié)果體現(xiàn)了本發(fā)明所提方法實(shí)現(xiàn)了自由度的增加。

仿真實(shí)例2：采用互質(zhì)陣列接收入射信號，其參數(shù)同樣選取為m＝3，n＝5，即架構(gòu)的互質(zhì)陣列共包含m+n-1＝7個物理天線陣元；假定入射窄帶信號個數(shù)為2，且入射方向?yàn)?0.5°至0.5°，其余參數(shù)設(shè)置與仿真實(shí)例1保持一致。由圖7所示的歸一化空間譜可以看出，本發(fā)明所提方法可以有效地分辨出這兩個近距離的信號源波達(dá)方向，說明了本方法良好的分辨率性能。

綜上所述，本發(fā)明所提方法充分利用了非均勻虛擬陣列上的全部信息，能夠在信號源個數(shù)大于等于物理天線個數(shù)的情況下實(shí)現(xiàn)入射信號的有效估計(jì)，增加了doa估計(jì)的自由度和分辨率。此外，與傳統(tǒng)采用均勻線性陣列的方法相比，本發(fā)明所提方法在實(shí)際應(yīng)用中所需的物理天線陣元及射頻模塊也能夠相應(yīng)減少，體現(xiàn)了經(jīng)濟(jì)性和高效性。