基于多維傅里葉變換的陣列空時信號處理與應用

石 榮,胡 柱,謝佳霖

(1.電子信息控制重點實驗室,四川 成都 610036;2.電磁空間安全全國重點實驗室,四川 成都 610036)

0 引 言

傳統的相控陣天線工作過程中首先進行空域處理,即通過各個單元天線所接收信號的加權求和形成指向目標信號方向的主波束,以此合成一路信號,然后再對該路信號進行后續的時域處理[1-2]。近年來隨著技術的發展,對陣列天線接收的信號進行空時聯合處理在工程上的應用越來越廣泛[3-4],在此情況下,其數據通常以矩陣形式來表達,以一維線陣產生的二維空時采樣數據為例,矩陣的列數與陣列中接收天線單元的個數一致,矩陣的行數與觀察時段內每個接收通道的時域采樣點個數一致。這意味著每一列數據對應了一個天線單元所接收到的信號,每一行數據對應了同一個采樣時刻對所有單元天線接收信號的采樣結果,所以該二維矩陣的2個維度分別對應了時間與空間。截止目前,所發展的大量空時聯合處理算法都以二維采樣矩陣為基礎,通過各種運算來完成時空域目標信號的檢測、來波方向估計與濾波等功能[5-8]。在工程實現中,對于大中型陣列來講,由于空時聯合處理所涉及的采樣信號通常為二維甚至三維矩陣數據,矩陣維數較大,巨大的計算量使其難以在實時性工程任務中應用。

針對上述問題,本文將圖像處理中的多維傅里葉變換方法引入到空時信號處理問題中,構建了陣列空時信號處理的工程實現新框架。圖像處理中的灰度圖像數據是二維矩陣形式,多光譜圖像數據是由不同光譜接收傳感器形成的高維矩陣形式,而多維傅里葉變換是圖像數據進行頻譜分析、頻率分量增強與濾波的常用處理手段。借鑒圖像處理領域的分析方法,將陣列天線接收的信號采樣數據看成是一幅圖像,一維線陣輸出類似于二維圖像的數據,二維面陣輸出類似于三維圖像的數據,通過多維傅里葉變換對上述空時采樣數據進行處理,同樣能夠提取出信號的時頻特征與來波方向等信息。所以本文在對信號時域與空域對偶關系概要總結之后,將傅里葉變換從時域推廣至空域,指出了空域與角度域之間存在的傅里葉變換關系,然后分別對線陣與面陣的空時聯合處理的二維與三維傅里葉變換的表達及其物理意義進行了詳細的闡述,并通過示例性仿真驗證了前述理論分析的正確性與有效性。

1 將傅里葉變換從時域推廣至空域

通過傅里葉變換將時域信號變換至頻域而得到信號的頻譜,對此大家早已熟悉。但是對空域采樣信號實施傅里葉變換,并清晰地解釋其物理意義,這還需要從時域與空域之間的對偶關系談起。

1.1 時域采樣定理與空域采樣定理

對一個連續時間信號x(t)進行時域離散采樣,當采樣頻率大于信號最高頻率的2倍時,采樣之后的離散信號完整地保留了原始信號中的所有信息。從理論上講,通過采樣信號能夠完整地恢復原始信號x(t),這就是大家十分熟悉的時域采樣定理,又稱為奈奎斯特采樣定理。

實際上在空域中也有相應的采樣定理,即對空間中傳播的信號x(t,d)在某一時刻t0的空間樣本x(t0,d)進行采樣,只有當空間上的采樣間隔小于信號波長的一半,即λ/2時,采樣之后的空間信號才能完整地保留原始信號x(t0,d)中與空間傳播方向相關的所有信息。換句話講,空域采樣定理要求在1個信號波長的空間距離上至少要有2個采樣點。如果將單位長度上的波數1/λ看成是一種空間頻率,那么空域采樣定理同樣要求空間采樣頻率大于信號空間頻率的2倍。

由上可見,信號的時域采樣定理與空域采樣定理具有對偶性,這也為將時域信號處理方法推廣應用于空域提供了條件。

1.2 時域與頻域之間的傅里葉變換關系

對一個信號進行時域處理一般要求信號帶寬與伴隨的噪聲帶寬基本一致,這樣才能獲得較好的效果,所以通常采用與信號帶寬基本相同的濾波器對采樣信號進行濾波,濾除帶外噪聲與其他干擾之后就能滿足這一條件。于是將信號由時域變換至頻域,通過頻譜分析得到信號的帶寬等參數。具體來講,在給定的觀察時段內,時域采樣信號x(n),n=0,1,…,N-1,與其頻域頻譜X(k1)(k1=0,1,…,N-1)之間滿足如下傅里葉變換關系:

(1)

在工程上,式(1)通常由快速傅里葉變換來實現,這也是信號分析處理工程應用中最常見的手段。如果將傅里葉變換應用于空域,其變換結果所展現的物理意義,下面以一維線陣為例進行分析。

1.3 空域與角度域之間的傅里葉變換關系

空間中一維均勻線陣包含N個單元天線,記單個信號入射到陣列上時所輸出的信號為ys(n,m):

ys(n,m)=yt(nTs-2πmdscosθs/λ,0)

(2)

式中:第1個變量n表示時域采樣序號,n=0,1,…,N-1,N為時域中總的采樣點數;第2個變量m表示空域采樣序號,m=0,1,…,M-1,M為線陣中單元天線的個數;yt(t,0)=xt(t),表示第0號單元天線接收到的連續時間域的信號;Ts表示時間采樣間隔;ds表示空間采樣間隔,即相鄰2個單元天線之間的距離;λ為信號波長;θs表示信號來波方向與線陣基線之間的夾角,且θs∈[0,π]。

在實際工程應用中的大多數情況下,信號帶寬遠小于信號載頻,滿足窄帶接收條件,于是式(2)可簡化為:

ys(n,m)≈xt(nTs)·exp(-j2πmdscosθs/λ)

(3)

如果對某一時域采樣時刻ns0的M個單元天線采集的信號ys(ns0,m)進行空域傅里葉變換,其結果記為Y(ns0,k2),k2=0,1,…,M-1,于是可得:

(4)

參照時域傅里葉逆變換表達方式,得到與式(4)相關的空域傅里葉逆變換:

(5)

空域傅里葉變換與逆變換中變量k2的物理意義主要體現在信號來波方向的角度信息上,所以空域與角度域是1對傅里葉變換關系式。

時域信號通過傅里葉變換變換至頻域,得到信號的頻譜;將空域信號通過傅里葉變換變換至角度域,可以得到信號的角度譜。信號頻譜反映了信號所包含的不同頻率分量的分布與構成情況;而信號的角度譜則反映了信號入射到一維線陣上的來波方向角度的特征信息。將式(3)代入式(4)可得:

Y(ns0,k2)=xt(ns0Ts)·

(6)

式(6)實際上描述了變量k2與信號來波方向角θs之間的對應關系,由于k2取值[0,M-1]范圍內的整數,而來波方向角θs的取值范圍為[0,π]。所以,在對上述空域傅里葉變換后所得到的角度譜的物理意義解釋,需要特別注意坐標單位的平移對應關系。

(1)當θs∈[π/2,π]時,有:

(7)

由式(6)可知:在k2對應的角度θs方向上相干匯聚了整個一維均勻線陣M個單元天線接收的來波信號,并且會在角度θs上出現一個峰值。當k2從0變化到M/2時,由式(7)可解得:

θs=arccos(-k2λ/(Mds))

(8)

由式(8)可知,當ds=λ/2時,k2從0變化到M/2對應了θs從π/2變化到π。

(2)當θs∈[0,π/2]時,有:

(9)

同樣由式(6)可知,在k2所對應的角度θs方向上相干匯聚了整個一維均勻線陣M個單元天線接收的來波信號,并且會在角度θs上出現1個峰值。于是當k2從M/2變化到M-1時,由式(9)可得:

θs=arccos((M-k2)λ/(Mds))

(10)

由式(10)可知,k2從M/2變化到M-1對應了θs從0變化到接近π/2。

由于離散傅里葉變換具有周期性,如果以k2=0時對應的來波方向角θs=π/2為中心,按照周期性,k2從M/2變化到M-1等效對應了k2從-M/2變化到接近于0,于是來波方向角θs從0變化到π就與k2從-M/2變化到M/2完全對應起來了。

在以上空域與角度域的變換關系分析過程中采用的是單個信號,如果有多個信號從不同角度方向同時入射到一維線陣上,針對每一個信號而言,式(2)均成立,在窄帶信號應用條件下,陣列中所有單元天線輸出的信號y(n,m)可表示為:

(11)

式中:Ls表示同時入射到陣列中的信號個數;θs,i表示第i個信號的來波方向線與陣列基線之間的夾角。

將式(11)代入式(4),可得某一時域采樣時刻ns0的M個單元天線采集的信號y(ns0,m)的傅里葉變換結果Y(ns0,k2)。時域信號的傅里葉變換滿足線性疊加原理,空域信號的傅里葉變換同樣滿足線性疊加原理,所以式(4)與式(5)所表示的空域與角度域之間的傅里葉變換與逆變換表達式對同時入射到陣列上的多個信號仍然成立。

2 線陣空時聯合處理的二維傅里葉變換表達及物理意義

在給定的觀察時段內,一維線陣采樣得到的空時數據的二維矩陣,其元素記為y(n,m),其中自變量n與m分別代表時域與空域采樣序號,矩陣大小為N×M。對y(n,m)進行二維傅里葉變換的結果記為Y2(k1,k2):

Y2(k1,k2)=

(12)

如前所述,Y2(k1,k2)展示了時域與空域相對應的頻域與角度域的信號分量分布情況,其中第1個自變量k1代表頻域分量序號,反映了信號頻率分量的分布,第2個自變量k2代表角度域分量序號,反映了不同信號的一維來波方向,如圖1所示。

圖1 線陣數據的二維傅里葉變換中各物理量的轉換關系

在計算實施方面,將式(12)改寫成如下形式:

(13)

由式(13)可知,空時采樣數據矩陣的二維傅里葉變換能夠通過(N+M)次一維傅里葉變換來實現,前N次一維傅里葉變換是針對y(n,m)的M列數據進行,分別將各個單元天線接收的時域數據變換至頻域,從而得到1個相同維度的中間形式的二維矩陣,其元素記為Yy(k1,m),然后再對Yy(k1,m)的N行數據進行N次一維傅里葉變換,從而最終得到Y2(k1,k2)。而上述N+M次一維傅里葉變換均能夠通過快速傅里葉變換來實現,所以在工程應用中線陣信號的二維傅里葉變換的計算量并不大。特別是對于大中型陣列天線來講,本方法相對于傳統的基于信號互相關矩陣分解的空間譜之類的方法,計算量更小,并具有可實時計算的優勢。

在物理意義上,變換結果的模值平方‖Y2(k1,k2)‖2完整地反映了不同頻率信號從不同角度到達一維線陣的能量分布情況。通過‖Y2(k1,k2)‖2不僅可實現信號的頻域檢測、頻率與帶寬參數的估計,而且還能估計各個信號的來波方向角度參數。

式(12)和式(13)的運算對信號并沒有附加特定的限制條件,即便是頻率相同的多個信號同時入射到一維線陣中時,只要這些信號的來波方向不同,在線陣陣元個數足夠多的情況下,也能夠將各個同頻信號的來波方向信息逐一提取出來。這一點是線陣空時聯合處理中二維傅里葉變換相對于傳統的基于陣列信號互相關矩陣分解的MUSIC類空間譜估計算法的優勢所在。在MUSIC等算法中,為了避免同頻多信號所引入的互相關矩陣的秩虧損等問題,需要采用子陣平滑等操作,不僅計算繁雜,而且性能還會下降;而本方法直接通過快速傅里葉變換即可分辨同時到達的多個同頻信號。

3 面陣空時聯合處理的三維傅里葉變換表達及物理意義

在工程應用中,陣列形式有線陣和面陣,面陣相對于線陣增加了一個空間維度,其空時采樣數據由二維矩陣變成了三維矩陣,包含1個時間維度和2個空間維度,可記為y(n,m,l),其中第1個自變量n代表時域采樣序號,第2個自變量m代表橫向空域采樣序號,第3個自變量l代表縱向空域采樣序號,整個矩陣大小為N×M×L。對y(n,m,l)進行三維傅里葉變換的結果記為Y3(k1,k2,k3):

Y3(k1,k2,k3)=

(14)

式中:k1代表頻域分量序號,反映了信號頻率分量的分布;k2與k3分別代表了二維角度域分量序號,反映了不同信號的二維來波方向,如圖2所示。

圖2 面陣數據的三維傅里葉變換中各物理量的轉換關系

作為前一節相關理論分析結果的推廣,對面陣空時采樣數據進行三維傅里葉變換同樣可以通過(NL+LM+MN)次一維快速傅里葉變換來實現,在工程上可實時實現。2個空間維度變換之后也對應了2個角度維度的數據,從物理意義上講,分別表示了各個信號來波方向與面陣相互垂直的2條基線之間的夾角。二維平面上的信號來波方向用1個角度參數即可描述,而三維空間中的信號來波方向需要用2個角度參數來描述,例如工程上常常用來波方向的方位角與俯仰角來描述,在此所使用的2個角度參數分別是來波方向與面陣2條正交測量基線之間的夾角,由立體幾何可知該夾角與方位角、俯仰角之間可以相互換算。于是三維傅里葉變換結果的模值平方‖Y3(k1,k2,k3)‖2同樣完整地反映了不同頻率的信號從三維空間不同角度到達二維面陣的能量分布情況。通過‖Y3(k1,k2,k3)‖2同樣能夠實現信號的頻域檢測、頻率與帶寬參數的估計、對各個信號的來波方向角度參數的估計。而且對于同頻多信號同時達到二維面陣的情況在陣元個數足夠多的情況下也可進行分辨。

4 仿真與應用示例

(1)一維均勻線陣采樣信號的二維傅里葉變換

由128個單元天線構成的一維均勻線陣中,相鄰單元天線的間距為0.12 m,有3個相等功率的1.25 GHz的同頻信號分別從不同的方向入射,第1個為脈寬0.5 μs的單頻雷達脈沖信號,第2、3個為連續波瞄頻干擾信號,上述3個信號的來波方向線與線陣基線之間的夾角分別為57°,106°,113°。以3 GHz采樣頻率在信號持續存在期間進行連續采樣0.2 μs,從而得到一個600×128的二維采樣信號矩陣。

按照前述的線陣空時聯合處理的二維傅里葉變換方法,首先進行128次列方向上的一維時域傅里葉變換,然后對變換后的數據再進行600次行方向上的一維空域傅里葉變換,從而將時域與空域的采樣信號矩陣變換成了頻域與角度域的二維矩陣,各元素求模之后如圖3所示。

圖3 二維傅里葉變換后展示的頻域與角度域結果

圖3中的3個峰值分別對應了3個信號,在頻率維度3個信號的頻率相同,均為1.25 GHz,在角度維度3個信號的來波方向分別為56.8°,106.3°,113.0°。由此可見,從一維線陣采樣的空時信號的二維傅里葉變換中不僅可以獲得各個信號的頻率,而且還能對同頻多信號的不同來波方向進行估計。

(2)二維均勻面陣采樣信號的三維傅里葉變換

將前述一維線陣擴展成128×128的二維均勻面陣,相鄰單元天線的間距仍然為0.12 m,同樣是前述3個信號入射,入射角度與面陣的2條正交基線的夾角分別是:(57°,120°),(106°,45°),(113°,98°)。同樣以3 GHz采樣頻率在信號持續存在期間進行連續采樣0.2 μs,從而得到一個600×128×128的三維采樣信號矩陣。

按照前述的面陣空時聯合處理的三維傅里葉變換方法,首先進行128×128次列方向上的一維時域傅里葉變換,然后對變換后的數據再進行600次二維空域傅里葉變換,從而將時域與空域的采樣信號矩陣變換成了頻域與角度域的三維矩陣。每一次二維空域傅里葉變換還可以繼續分解成256個一維傅里葉變換,所以總的計算量并不大。由于三維矩陣難以在平面上展示,所以用切分方式進行局部觀察,以頻域頻率采樣序號等于250的頻點對三維矩陣進行切片,該切片對應的信號頻率為1.25 GHz,從而得到一個128×128的二維矩陣,該矩陣的2個維度分別是來波方向線與2條正交基線之間的夾角,該二維矩陣的各元素求模之后如圖4所示。

圖4中的3個峰值分別對應了3個信號,在頻率維度3個信號的頻率相同,均為1.25 GHz,在角度維度3個信號的二維來波方向分別為(56.8°,120.0°),(106.3°,45.3°),(113.0°,98.1°)。由此可見,從二維面陣采樣的空時信號的三維傅里葉變換中不僅可以估計出各個信號的頻率,而且還能對同頻多信號的不同來波方向進行準確估計。

5 結束語

在大中型陣列天線所接收信號的空時信號處理過程中,由于陣元數量眾多,導致空時采樣數據矩陣的維數很大,高維矩陣運算所需巨大計算量給實時處理類的工程應用帶來了極大的困難。本文借鑒圖像處理中的多維傅里葉變換方法,將時空域采樣數據進行多維傅里葉變換,而這些變換均能夠分解并利用更低維數的快速傅里葉變換來實現,從而為工程上的實時處理創造了條件。而且空域與角度域的傅里葉變換關系也使得該方法將大中型陣列的空時數據轉換成了頻域與角度域的數據,實現了信號的頻域參數提取與來波方向估計。最后通過仿真驗證了理論分析的正確性,同時也展現了線陣與面陣各自空時聯合處理的二維與三維傅里葉變換的本質物理意義,從而為大中型陣列的空時信號實時處理工程應用提供了新的技術手段。