基于旋轉干涉儀的多目標參數估計

申皓明, 廖桂生, 楊志偉

(西安電子科技大學雷達信號處理國家重點實驗室, 陜西 西安 710071)

0 引 言

隨著需求的增長以及科技的發展,電子偵察衛星已經成為獲取各種情報的重要手段,在現代戰爭中有著無可替代的地位[1]。目標定位是電子偵察衛星的一項重要應用,其獲得的目標位置是偵察信息中最為重要的信息之一。而如何獲得高精度的多目標定位結果是提高定位系統性能的關鍵問題。

對于傳統的處理方式:①單星測角定位,其優點是系統易于實現,但由于天線孔徑和載荷的限制,測向精度不高[2-4],在文獻[5] 提出了一種基于粒子群算法的無源定位技術,有相對較好的結果,但是計算復雜度較高,文獻[6]提出了基于脈沖信號到達時間的單衛星無源定位方法,利用網格搜索的方法對輻射源的位置進行估計,系統的定位精度較低,且需要長時間觀測,對高速運動目標定位精度較差,無法對短時目標進行定位。②相位干涉儀測向系統具有測向精度高、系統實現復雜度低、計算復雜度低、可實時處理等優勢,但該方法隨著測角精度的不斷提高會造成角度模糊范圍的不斷擴大[7-12],即干涉儀系統的基線越長,測向結果越精確,同時產生的角度模糊越嚴重,無法得到真實目標的角度信息。而基于旋轉干涉儀的測向解模糊算法能有效解決測角精度與角度模糊范圍的矛盾問題[13-18],利用基線旋轉不同的角度并測量相位差進行角度估計,但這類方法只對單目標情況,在多目標情況下,算法失效。③傳統的陣列信號處理方式可以解決多個目標角度估計的問題,但為避免角度模糊的影響,陣元間距不超過半波長。星載陣列為得到較高的角度分辨率對孔徑要求比較大,系統的實現要由很多的陣元組成。由于載荷限制,星載陣列陣元數目較少,陣列的孔徑太小,定位系統的性能受到很大的限制。④文獻[19]提出了一種基于旋轉干涉儀的多目標角度估計方法,可以對工作時長不同的同頻目標源信號進行檢測,通過各目標源在時間上的差異完成多目標角度的估計,但對于同頻同時工作的信號源無法進行分辨。

針對以上問題,本文提出了一種將多重信號分類(multiple signal classification, MUSIC)算法和霍夫變換[20-27]思想相結合的多目標參數估計算法。該方法可以在較大的陣元間距下,通過對角度模糊信息的提取,實現在較大的孔徑下利用少量陣元對多個目標進行角度估計。首先利用不同時刻的采樣數據通過MUSIC算法對信號進行一維譜估計,再將估計值按時間排列,不同目標的角度模糊隨時間按照一定的反三角函數曲線規律變化,接著通過霍夫變換的思想在二維角度參數空間中對特定曲線進行參數提取,實現多目標的二維角度估計。

1 系統模型建立

1.1 旋轉干涉儀模型

圖1為旋轉干涉儀模型,圖1中干涉儀以第一個陣元為參考陣元建立空間三維直角坐標系,兩個陣元間距為d,其中d?λ/2。設陣列天線在xOy平面內繞第一個陣元沿逆時針方向以角速度ω做勻速圓周運動。假設信號為遠場窄帶信號,以x軸與y軸為參考, 輻射源相對于坐標系的俯仰角為β,其中,俯仰角的范圍為[0,π/2) ,表示和z軸正方向的夾角;方位角為α,其中,方位角的取值范圍為[0,2π),表示和x軸正方向的夾角(逆時針方向旋轉)。在基線旋轉過程中,兩天線間的相位差隨基線與輻射源夾角的變化而變化。

圖1 旋轉干涉儀模型Fig.1 Rotating interferometer model

由圖1可知,相對與參考陣元,陣元接收到信號為

x(t)=s(t)ejφ(t)+σn

(1)

式中,s(t)表示信號源的包絡;σn為陣元1上的高斯白噪聲,服從正態分布N(0,σ2);φ(t)表示在t時刻陣元2與陣元1之間的相位差。兩個天線單元間接收到信號的相位差為

(2)

式中,λ為輻射源波長;β,α分別為信源的入射俯仰角和方位角;d為干涉儀陣元間的間距,其中間距d?λ/2。由式(9)可以看出,當俯仰角β固定不變時,相位差受天線基線旋轉的影響呈余弦規律變化。

1.2 多通道旋轉干涉儀模型

圖2為多通道旋轉干涉儀模型,圖2中干涉儀變為多個接收通道,陣元間距為d。設陣列天線在xOy平面內繞第一個陣元沿逆時針方向以角速度ω做勻速圓周運動。

圖2 多通道旋轉干涉儀模型Fig.2 Multichannel rotating interferometer model

由圖2可知,相對與參考陣元,剩余天線單元接收到信號為

xi(t)=s(t)ejφi(t)+σn

(3)

式中,s(t)表示信號源的包絡;σn為陣元的高斯白噪聲,服從正態分布N(0,σ2);φi(t)為t時刻陣元i與第一個陣元間的相位差。兩天線接收到信號的相位差為

(4)

假設旋轉干涉儀有N個陣元,空間中有L個信號,其中信號的方位角,俯仰角分別為α1,α2,…,αL,β1,β2,…,βL,可以得到干涉儀陣元間相位差矩陣為

Φ(t)=[Φα1,β1(t),Φα2,β2(t),…,ΦαL,βL(t)]

(5)

(6)

將式(6)化簡為

X(t)=Φ(t)S(t)+N(t)

(7)

式中,X(t)為陣列采樣得到的接收數據;Φ(t)為N×L維陣元間相位差變化矩陣;S(t)表示空間信號;N(t)表示噪聲矩陣。

1.3 目標角度模糊分布模型

根據第1.2節所述模型,多通道旋轉干涉儀天線陣元數目為N,信源個數為L,信號俯仰角為β,方位角為α,陣元間距為d?λ/2,在陣列信號處理中會出現角度模糊,假設d=m(λ/2),則兩個陣元間相位差為

mπsinβcos(ωt-α)

(8)

由式(4)可知，距旋轉線陣的陣列流型陣隨時間呈余弦變化,因為衛星天線旋轉的角速度ω遠大于采樣周期Ts,所以利用MUSIC算法對目標角度估計時,在快拍積累時間內認為陣流型是不變的。

在ti時刻對目標的俯仰角進行一維的譜峰搜索,搜索導向矢量只考慮俯仰角的變化,即

(9)

對exp(-j2πdsinβcos(ωti-α)/λ),其中sinβcos(ωti-α)的取值范圍為(-1,1),所以θ∈(-90°,90°)。譜峰搜索變為求解式(10)。

(10)

將d=m(λ/2)代入得

(11)

exp(-jmπsinβcos(ωti-α))

(12)

在對譜峰搜索的過程中,θ∈(-90°,90°),則

(13)

式中,(πsinθ±2kπ)∈[-mπ,mπ],即k∈[(-m-sinθ)/2,(m-sinθ)/2]。代入式(7)得

exp(-j(πsinθ±2kπ))=exp(-jmπsinβcos(ωti-α))

(14)

化簡得

(15)

最終得到的一維譜估計值為

(16)

式中,θ(ti,k)表示在ti時刻譜峰搜索中搜索到的第k個峰值;ω為天線旋轉角速度;2k/m為曲線基線位置,式(16)說明在一次MUSIC角度估計中,由于陣元間距d=m(λ/2)?2/λ,譜峰搜索過程中會出現多個峰值,其中，由式(13)可知,峰值個數M=2d/λ=m,圖3為當d=8(λ/2),信號源數目為2時,利用MUSIC算法進行俯仰維的一維譜峰搜索。

圖3 MUSIC方法角度估計結果Fig.3 Angle estimation results by MUSIC method

整個天線旋轉過程,譜峰隨時間的分布為

(17)

由式(17)可知,峰值位置隨時間變化關系為反三角函數與三角函數構成的復合函數。由于θ(t,0)是非線性函數,簡便運算復雜度,只分析k=0的情況。在θ(t,0)中,當k=0時,θ(t,0)∈[-β,+β],所以β為θ(t,0)曲線的幅度,α為曲線的初相。如圖4所示,信號源數目為2,陣元間距分別為d=λ/2和d=8(λ/2),將天線旋轉一個周期內多次估計的譜峰值提取出,按時間排列在一起。

圖4 天線旋轉一周的MUSIC估計Fig.4 MUSIC estimation of antenna rotation with one cycle

在圖4(a)中,陣元間距為半波長,測向過程中不會產生角度模糊,此時曲線的k=0,即角度曲線的基線為零,圖4(a)中兩條曲線分別對應相應的兩個信號源,其中每條曲線的幅度和初相均為相應信號源的俯仰角β,方位角α。在圖4(b),陣元間距d=8(λ/2),故在除了基線位置為零的曲線外,仍出現多條反三角函數曲線,圖4(b)中這些基線不為零的曲線均為測向過程中角度模糊對應的峰值曲線,因為在單次譜峰搜索過程中無法區分所得譜峰是真實值還是模糊值,因此真實值不能單獨提取出,只能將模糊值與真實值全部列出。綜上只需要對所有基線位置為零的曲線進行參數估計,得到相應曲線的幅度和相位信息,即可完成所有信號源的俯仰角β、方位角α的估計,實現多目標角度解模糊。

2 本文算法原理

提取的峰值曲線表達式為

(18)

在參數估計中,化簡參數積累過程,只需考慮2k/m=0的情況。則式(18)變為

θ(t)=arcsin(sinβcos(ωt-α))

(19)

式中,因為天線旋轉角速度ω已知,在曲線參數估計中只需對俯仰角β、方位角α進行估計。

下面對峰值曲線的Hough變換原理作以說明。式(19)中β,α為常數,分別代表曲線的幅度和初相。對曲線上任意一點(ti,θ(ti))來說,在θ(t)-t構成的二維平面中有無數條曲線可以通過該點。若將θ(ti),ti視為常數,β,α視為變量,其中β∈(0～90°),α∈(0～360°),則式(19)變為

(20)

通過式(19)的變換,將峰值曲線上的點(ti,θ(ti))變換為α-β的二維參數空間中的一條曲線,即完成對(ti,θ(ti))點的Hough變換。Hough變換原理如圖5所示,圖5中θ(t)-t的二維坐標系中峰值曲線為θ(t)=arcsin(sinβ1cos(ωt-α1)),在峰值曲線中任取3個點:(t1,θ(t1))、(t2,θ(t2))和(t3,θ(t3)),將這3個點分別代入式(20)中,可得到參數空間中對應的曲線,3條曲線有共同的交點(α1,β1)。因此將峰值曲線上所有點都投影到參數空間,每個點在參數空間都有對應的曲線,這些曲線在參數空間有共同的交點(α1,β1),對參數空間中曲線上每個點的出現次數進行統計,統計次數最多的點對應的坐標位置便是峰值曲線的參數。

圖5 Hough變換原理Fig.5 Hough transform principle

假設空間中L個信號源的峰值曲線向參數空間投影:①凡是k/m=0的峰值曲線在向式(17)構造的參數空間投影后都可以在相應參數點(αl,βl)處積累出一個極大值,其中,l=1,2,…,L,由此實現多目標的檢測;②各極大值點的坐標(αl,βl)分別為各峰值曲線的初相和幅度,即對應各信號源的方位角和俯仰角,對于一個信號源而言,得到一個極大值點的坐標即可同時得到方位角和俯仰角的估計,避免了信號源間的參數配對問題;③對于k/m≠0的曲線向式(17)構造的參數空間中投影沒有固定的交點,所以角度模糊產生的峰值曲線在參數空間中無法得到有效積累,由此避免了角度模糊影響,實現目標的解模糊。圖6給出了多目標情況下的Hough變換示意圖。

圖6 多曲線Hough變換示意圖Fig.6 Schematic diagram of multiple curve Hough transform

本文算法具體檢測方案流程如圖7所示。

圖7 數據處理流程圖Fig.7 Data processing flow chart

算法步驟如下:

步驟1設參數空間中的角度精度為P,參數空間累加矩陣CA[β][α],維數為(90/P)×(360/P),初始時刻所有元素置零。

步驟2令方位角α從0到360°以角度精度P變化,利用峰值點θ(t)和對應時刻t代入式(17)計算不同α取值下俯仰角β的值:

步驟3若β在(0～90°)范圍以內,則在CA[α/P][β/P]位置處累加1。若β不在(0～90°)范圍以內,則舍棄該次計算結果。直到干涉儀旋轉過程中所有峰值點都被計算過。

步驟4對CA[β][α]矩陣進行遍歷,搜索得到最大的L個峰值,每個峰值的坐標位置代表該目標的方位角,俯仰角信息。第l個峰值的所處位置的值分別該目標的方位角αl,俯仰角βl。完成多信號源的二維角度解模糊,實現目標定位。

算法流程圖如圖8所示。

算法復雜度分析,在算法計算過程中,P為角度搜索精度,L為信號源個數,m為陣列天線的陣元間距對半波長的倍數,設天線旋轉過程中角度估計頻率為fE,天線旋轉時間為t,在系統結構確定時,L,m,fE,P這4個參數為常數,則算法計算復雜度為O(LmfE/P·t),即復雜度是O(t)，為線性階復雜度。算法復雜度較低,天線旋轉過程中,在每完成一次一維角度估計后即可將結果向參數空間投影做參數積累。

圖8 算法流程圖Fig.8 Algorithm flow chart

3 仿真實驗與分析

本小節在零均值高斯白噪聲背景下進行,通過不同目標二維角度估計仿真實驗來驗證本文方法對多目標角度解模糊的有效性。

3.1 曲線參數檢測算法有效性仿真分析

在本次仿真中,主要驗證算法在曲線檢測過程中只能對基線為零的曲線進行有效檢測,對于基線非零的曲線在參數空間中無法進行有效積累。本仿真分為兩部分,在第一部分中根據式(19)產生一條基線為零的峰值曲線,并對其進行曲線參數檢測,零基線處峰值曲線參數如表1所示;第二部分根據式(18)產生3條基線不為零的峰值曲線,這3條曲線與第一部分曲線的幅度和初相一致,并進行曲線參數檢測,非零基線處峰值曲線參數如表2所示。仿真圖如圖9和圖10所示。仿真第一部分中由圖9和圖10 可知:對于基線為零的曲線,算法在參數空間對其可進行有效積累,圖10(a)為參數空間積累結果,在曲線對應參數處有明顯的峰值出現,將圖10(a)向幅度維和相位維投影可得到圖10(b)和圖10(c),可知,曲線的幅度和初相檢測結果與表1中仿真參數一致。

表1 零基線曲線仿真參數

表2 非零基線曲線仿真參數

仿真第二部分中由圖11和圖12可知:對基線非零處的峰值曲線,算法在參數空間對其無法進行有效積累,在圖12參數空間中沒有峰值出現。

圖9 零基線處峰值曲線仿真Fig.9 Simulation of peak curve at zero baseline

圖10 零基線處峰值曲線參數空間處理結果Fig.10 Processing result of parameter space in peak curve at zero baseline

圖11 非零基線處峰值曲線仿真Fig.11 Simulation of peak curve at non zero baseline

圖12 非零基線處峰值曲線參數空間處理結果Fig.12 Processing result of parameter space in peakcurve at non zero baseline

根據第2節分析可知,基線為零的峰值曲線對應的是信號源的真實角度信息,基線不為零的峰值曲線對應的是信號源的模糊角度信息,因此算法能出檢測出真實角度信息,同時消除了角度模糊的影響。

3.2 陣元間距半波長仿真分析

本仿真主要驗證算法對多信號源參數檢測的有效性,系統參數及目標參數如表3和表4所示。

表3 系統仿真參數

表4 目標參數

圖13、圖14給出了算法在陣元間距為半波長時對多目標的角度估計。一維譜估計使用MUSIC算法實現,快拍數為200。如圖13所示,3條峰值曲線的基線2k/m=0。參數空間積累結果如圖14(a)所示,參數空間積累中出現3個信號的峰值,各峰值對應坐標為(20°,150°)、(30°,40°)、(55°,210°),其中橫坐標為峰值曲線的幅度對應信號源的俯仰角,縱坐標為曲線的初相對應信號源的方位角,由此根據坐標點位置即可得到信號源二維角度信息,避免了參數匹配問題。圖14(b)和圖14(c)是圖14(a)在方位維和俯仰維的投影。

圖13 d=λ/2時天線轉動一周的譜估計峰值Fig.13 Spectral estimation peak of antenna rotation at array element d=λ/2

3.3 陣元間距多倍半波長仿真分析

本小節分析算法在測向結果有角度模糊時在參數空間積累的有效性,仿真利用Matlab對整個信號處理流程進行仿真,系統參數與目標參數如表5和表6所示。

圖14 陣元間距半波長時參數空間處理結果Fig.14 Processing result of parameter space in half wavelength of array element

參數數值載頻/GHz3天線孔徑/m5陣元間距/m0.5(10倍半波長)陣元個數10采樣頻率/MHz40角度估計頻率/Hz300天線旋轉周期/s1天線旋轉時長/s1

表6 目標參數

圖15、圖16給出了本文方法在陣元間距為多倍半波長情況下對多目標二維角度的估計驗證。一維角度估計使用MUSIC算法實現,其中MUSIC算法快拍數為200。

圖15 d=10(λ/2)時天線轉動一周的譜估計峰值Fig.15 Spectral estimation peak of the antenna rotation at array element d=10(λ/2)

圖16 陣元間距多倍半波長時參數空間處理結果Fig.16 Processing results of parameter space in multiple wavelength half space of array elements

由圖15可見,天線相對目標角度成多條三角復合函數曲線變化,其中,陣元間距為10倍半波長,因此對于同一個目標會產生多個角度模糊。由式(18)可知,對于同一個信號源而言,峰值曲線有兩種,分別為真實峰值曲線與角度模糊峰值曲線,這兩種曲線的初相和幅度是相同的,區別在于基線分布的位置不同,真實峰值曲線分布在基線為零的位置,即2k/m=0,而角度模糊峰值曲線的基線分布在非零處,即2k/m≠0。

在曲線參數檢測過程中:①對真實峰值曲線進行了有效檢測,參數空間積累結果由圖16(a)所示,參數空間積累中出現3個信號峰值,3個峰值所處位置的坐標分別為:(20°,150°)、(30°,40°)、(55°,210°),可以準確分辨各目標的角度信息,無需對信號源進行角度匹配。圖16(b)和圖16(c)是圖16(a)分別在方位維和俯仰維的投影。②對于基線位置非零的角度模糊峰值曲線,在投影到角度參數空間時不能產生有效積累,在參數空間中以小毛刺的形式存在,形成不了峰值,由此避免了角度模糊的影響。綜上本文所提處理方案可以有效解決測向過程中出現的角度模糊問題。

3.4 性能分析

在本小節中,利用本文所提方法與互質基線解模糊方法、傳統旋轉干涉儀方法和陣元間距為半波長的均勻線陣MUSIC測向方法進行對比,為了便于分析,通過單目標蒙特卡羅仿真實驗進行說明,多通道旋轉干涉儀系統參數、互質基線干涉儀系統參數、傳統旋轉干涉儀系統參數、陣元間距為半波長的均勻陣列系統參數和目標參數、分別如表7～表11所示。

表7 多通道旋轉干涉儀系統仿真參數

表8 互質基線系統仿真參數

表9 旋轉干涉儀系統仿真參數

表10 半波長均勻線陣系統仿真參數

表11 目標參數

輸入信噪比(signal-to-noise ratio, SNR)變化范圍為-5～20 dB,步長為1 dB,對每一個信噪比進行1 000次蒙特卡羅實驗。

圖17為4種算法的測向均方根誤差(root mean square error,RMSE)的比較,由仿真結果可得:①基于多通道旋轉干涉儀的方法利用子空間類測向算法進行處理性能優于互質基線解模糊方法和傳統旋轉干涉儀方法,且對信噪比的要求較低;②本文所提方法為稀疏布陣,陣元個數少,因此性能與陣元間距為半波長的均勻線陣相比較差。③本文所提方法為均勻稀疏布陣,對一定頻率范圍內的信號均可處理,對于低頻信號避免了由于陣元間距過小造成的互耦現象,對于高頻信號帶來的測向模糊,本文所提處理流程可有效解決。本文所提方案可適應電子偵察中對寬頻程信號偵測的要求。

圖17 4種方法在不同SNR下的性能對比Fig.17 Performance comparison of four methods under different SNR

4 結束語

本文采用的將傳統測角方式與圖像檢測相結合的目標檢測算法解決了在星載平臺下大口徑天線測向模糊的問題,同時可以應用到衛星編隊中。優點如下:①通過陣列信號的處理方式,實現了多個目標的分辨,準確度高,對SNR的要求大大降低;②利用旋轉基線的方式角度模糊呈三角曲線規律變化,避免不同目標間的角度匹配問題;③利用曲線參數積累的方式估計目標角度信息,運算復雜度低,可以進行并行處理,具有一定的實時性,便于工程實現。仿真實驗證明本文方法的有效性。針對稀疏天線角度解模糊問題,將在接下來的科研工作中進一步展開研究。