蔣伊琳,唐三強,陸滿君,張莉婷

(1． 哈爾濱工程大學 信息與通信工程學院, 黑龍江 哈爾濱 150001;2. 哈爾濱工程大學 先進船舶通信與信息技術(shù)工業(yè)和信息化部重點實驗室, 黑龍江 哈爾濱 150001;3. 上海無線電設備研究所, 上海 200000)

受硬件限制,單一子帶系統(tǒng)觀測的頻帶范圍受到限制,回波所含目標信息有限,由這些信息反演所得到的圖像分辨率也相應受到制約。為了提高距離分辨率需要提高信號帶寬,而受制于Nyquist采樣定理這將導致發(fā)射、接收系統(tǒng)更加復雜,給系統(tǒng)硬件實現(xiàn)增加了很大的難度。隨著雷達技術(shù)和信息處理技術(shù)的發(fā)展,為了克服單一子帶的局限,多子帶融合技術(shù)的綜合使用逐漸成為一種趨勢[1]。

多子帶融合的本質(zhì)是外推預測高低子帶之間的空余頻段[2]。傳統(tǒng)的多頻段融合方法可以分為兩類:

第一類是非參數(shù)方法,不需要目標的先驗信息。文獻[3]提出了放大間隙數(shù)據(jù)幅度和相位估計(gapped-data amplitude and phase estimation, GAPES)方法。該方法采用最小二乘法迭代估計未知頻譜。仿真和實測數(shù)據(jù)證實了該方法的有效性,但不適當?shù)某跏蓟赡軙乖摲椒ㄏ萑刖植孔顑?yōu)。文獻[4]將不同雷達的相位偏差建模為線性相位和恒定相位。該方法采用全相位快速傅里葉變換(all-phase fast Fourier transform, apFFT)對脈沖壓縮后的圖像進行相關(guān)處理,估計恒定相位和線性相位,但是,apFFT脈沖壓縮方式失去了信號初相信息,只能顯示小角度范圍內(nèi)的相位信息,在面對大角度相位差距時,會出現(xiàn)跨距離窗現(xiàn)象,造成估計誤差。文獻[5]提出了一種以稀疏表示為基礎的融合算法,通過相干補償與高分辨率成像相結(jié)合的方法來補償非相干相位,獲得高分辨率的逆合成孔徑雷達融合圖像。

第二類是參數(shù)化方法,建立參數(shù)化模型,求解相關(guān)參數(shù)。與非參數(shù)方法相比,參數(shù)化方法利用了豐富的先驗信息,具有更優(yōu)越的性能。文獻[6]提出了一種基于全極點模型的融合方法。該方法分別為較低和較高子帶建立前向預測矩陣,然后采用奇異值分解和Akaike信息準則(Akaike information criterion, AIC)估計極點數(shù),并采用最小二乘法求解模型參數(shù)。補償相位偏移后,得到積分頻率信號。然而,這種方法在低信噪比下難以確定模型階數(shù),模型階數(shù)錯誤對線性模型估計結(jié)果會造成極大偏差。文獻[7]研究了復合制導體制下多傳感器異步信息融合的時間同步和空間配準問題,并提出一種自適應無跡卡爾曼濾波算法,該算法采用預測殘差構(gòu)造狀態(tài)模型誤差統(tǒng)計量,通過自適應因子調(diào)整狀態(tài)模型信息對狀態(tài)參數(shù)估值的貢獻,有效控制狀態(tài)模型噪聲異常對狀態(tài)參數(shù)估值的影響。文獻[8]提出了一種基于廣義似然比(generalized likelihood ratio, GLR)檢測的多階段方法,用于使用稀疏子帶測量數(shù)據(jù)估計目標距離像。文獻[9]采用root-MUSIC算法和最小二乘法估計非相干相位(incoherent phase, ICP),然后利用相干處理對不同雷達之間的ICP進行補償,成功合成了四種光子雷達的雷達回波。文獻[10]提出了一種復雜噪聲環(huán)境下的多雷達子帶融合新算法,該算法基于噪聲概率密度函數(shù),引入懲罰函數(shù)來抑制不同類型的噪聲。文獻[11]提出了一種基于改進RELAX算法的多頻帶融合方法,該算法將最大差分準則應用于回波Hankel矩陣的奇異值,以提高計算散射中心數(shù)的精度。文獻[12]提出了一種基于合成實矩陣奇異值的模型階數(shù)估計方法,合成數(shù)據(jù)矩陣包括復雜觀測數(shù)據(jù)及其共軛數(shù)據(jù),充分利用了現(xiàn)有數(shù)據(jù)信息,該方法比傳統(tǒng)的基于復雜觀測數(shù)據(jù)矩陣的模型階數(shù)估計方法有更好的效果,可以在較低的信噪比下估計模型階數(shù)。文獻[13]提出了一種基于多輸出變量高斯過程模型的多波段圖像融合方法。

本文為了解決全極點模型在低信噪比下模型階數(shù)不準、性能較差的問題,將一種噪聲抑制方法引入模型,通過在建立模型過程中對信號Hankel矩陣主奇異值各分量加權(quán)和歸置等方法,達到抑制噪聲的目的,并采用整體前向預測矩陣對全極點模型參數(shù)進行整體估計,從而減少計算量。

1 噪聲抑制方法

針對文獻[6]在低信噪比情況下的不足進行改進,有效提高低信噪比條件下模型參數(shù)的估計精度。

將雷達回波信號與發(fā)射信號共軛相乘可以得到基帶信號,其頻率分量與目標上散射中心之間的相對距離成比例。對基帶信號進行采樣和逆傅里葉變換,以提供目標的距離分辨輪廓。靠近雷達的散射中心具有隨頻率衰減的散射振幅,而遠離雷達的散射中心具有隨頻率增加的散射振幅。以如下基帶信號s(n)[14]為例:

(1)

其中:Ak表示散射點散射強度;fk表示子帶起始頻率;fn=fk+n·Δf,n∈[0,Nk),Nk為頻率步進數(shù);am表示第m個散射中心類型;θm表示第m個散射點的角度。在帶寬不是很寬時可以將式(1)表示為全極點模型[15]。以信號s1為例抑制噪聲,首先構(gòu)建如式(2)所示Hankel矩陣:

(2)

其中,L1表示相關(guān)窗長度,Hankel矩陣與線性時不變系統(tǒng)的瞬態(tài)響應相關(guān)聯(lián)。子空間分解方法利用Hankel矩陣的特征結(jié)構(gòu)來估計線性時不變信號模型的參數(shù)[16]。為得到精確的模型參數(shù),通常采用L1=N1/3,當L1較大時雖然能得到更好的分辨率,但是對噪聲的適應能力較弱。

對H1矩陣進行奇異值分解[17],得到式(3):

(3)

S1=diag(σ1,σ2,…,σr)

(4)

其中,σi(i=1,2,…,r)為矩陣H1的所有奇異值,且σ1≥σ2≥…≥σr≥0,即S1為奇異值依次減小的對角矩陣。根據(jù)奇異值分解原理,分解后得到奇異值依次減小的對角奇異值矩陣S1。其中σi(i=1,2,…,r)代表不同的信號分量。通常接近于零的較小的σi代表噪聲分量,較大的σi代表目標散射點分量。將噪聲分量對應的σi根據(jù)需要取一個較小值或者歸置為零,即可實現(xiàn)對噪聲的抑制。如此可以得到新的對角奇異值S′1。

再重構(gòu)恢復出H′1矩陣:

(5)

利用式(5),對矩陣H′1求取矩陣H′1的反對角線元素的平均,即可重構(gòu)出抑制噪聲后的雷達信號序列:

(6)

其中,mn為求解的第n點時H′1矩陣中的反對角線上元素個數(shù),即符合i+j-1=n的H′1(i,j)元素個數(shù),∑H′1(i,j)表示每次對矩陣H′1的反對角線元素求和。

奇異值分解法抑制噪聲的步驟總結(jié)如下:

1)將雷達回波信號序列s(n)重排成Hankel矩陣;

2)對Hankel矩陣進行奇異值分解,即H=USVΗ,得到奇異值對角矩陣S;

3)根據(jù)奇異值跟蹤目標和噪聲等不同信號分量,然后將噪聲分量對應的奇異值歸置,即為抑制噪聲,得到新的奇異值對角矩陣S′;

4)用新的奇異值對角矩陣S′重構(gòu)H′矩陣,H′=US′VΗ;

5)將重構(gòu)的H′矩陣再進行重排,得到抑制噪聲后的信號序列s′(n)。

在實際多子帶融合過程中,只需對Hankel矩陣進行加權(quán)和對噪聲分量進行歸置即可,不用再恢復出抑制噪聲后的信號序列,這里是為了將噪聲抑制后信號與原信號進行對比。

2 多子帶融合方法

超寬帶處理要求每個子帶具有一致的頻譜信號,即每個子帶的全極點模型必須一致。當子帶測量值由獨立工作的寬帶雷達收集時,會出現(xiàn)互相干問題。本文應用相干函數(shù)的方法補償任意數(shù)量的子帶之間缺乏一致性的問題。

全極點模型是一種用指數(shù)函數(shù)近似冪指函數(shù)的方法,其原理是通過線性模型參數(shù)估計逼近非線性信號,算法整體流程如圖1所示。

圖1 算法整體流程Fig.1 Overall algorithm flow chart

可以將式(1)表示的基帶信號表示為如下所示的全極點模型:

(7)

其中,P表示極點數(shù)目,ak表示極點幅值系數(shù),pk表示模型的極點。通過將單獨的全極點模型擬合到每個子帶,并調(diào)整模型直到它們一致,可以使子帶相互一致。

設計多部同視角觀測的獨立工作的雷達為例,其起始頻率為fi(i=1,2,…,n),頻率步進數(shù)分別為Ni(i=1,2,…,n),各子帶表示為si。根據(jù)式(2)構(gòu)建各子帶的Hankel矩陣:

(8)

為了估計各個子帶的全極點模型參數(shù),對各個子帶的Hankel矩陣按式(3)進行奇異值分解,得到:

(9)

通過分解Hi,可以通過以下四步過程來估計每個子帶的全極點模型參數(shù):

步驟1：根據(jù)奇異值矩陣Si估計各個子帶的模型階數(shù)Pi。

步驟2：Pi將Vi劃分為相互正交的信號子空間和噪聲子空間。利用root-MUSIC算法估計每個子帶的信號極點。

步驟3：根據(jù)最小線性二乘擬合來確定全極點模型的幅值系數(shù)。

步驟4：根據(jù)極點值與極點幅值系數(shù)調(diào)整各個子帶互相干。

在步驟1中,Si中的奇異值用于估計兩個子帶的適當模型階數(shù)。相對較大的奇異值對應于強信號分量,而較小的奇異值通常對應于噪聲。對于低噪聲水平,在大的和小的奇異值之間有一個急劇的轉(zhuǎn)變,如圖2所示。這個轉(zhuǎn)變點可以用來估計模型階數(shù)。但是在較高的噪聲水平下,從大奇異值到小奇異值的過渡是非常平滑的,因此在低信噪比環(huán)境下精確估計模型階次會相當困難。

圖2 主對角奇異值曲線Fig.2 Main diagonal singular value curve

估計模型的方法有很多,其中比較經(jīng)典的方法有采用AIC和最終預測誤差(final prediction error, FPE )準則等,本文采用一種歸一化比值法進行模型估計。設σi(i=1,2,…,m)為矩陣S的奇異值,且有σ1≥σ2≥…≥σm=0,定義歸一化比值為

(10)

預先設定一個閾值(如0.990),當p是ρ(k)的值大于或等于該閾值的最小整數(shù)的時候,即可認為在矩陣的奇異值中,前p個為主要奇異值,即對應于強的信號成分的奇異值,而后面的奇異值即為次奇異值,對應于噪聲奇異值,從而可以將p定為該信號模型的階數(shù)。

(11)

Ai的第一行為dij,表示各個子帶的第一行的j個元素,構(gòu)成如下多項式:

(12)

Ai(z)的根即為各子帶的極點值。一般來講,信號模型的變化會導致極點在復平面上偏離單位圓,但是在每個子帶上,這個變化會非常小,所以可以認定主導信號成分的極點更靠近單位圓。于是可以選取最靠近單位圓的Pi個點作為各個子帶的極點,得到pi=[p1,p2,…,pPi]。

根據(jù)步驟3估計各子帶的全極點幅值系數(shù),可以歸結(jié)為一個線性最小二乘問題,根據(jù)式(13)即可得到各子帶極點幅值系數(shù)。

(13)

隨后可以根據(jù)式(7)表示各個子帶的值Mi(fn),根據(jù)步驟4調(diào)整各個子帶至互相干位置,根據(jù)相干函數(shù)

(14)

可以得到極點旋轉(zhuǎn)角度Δθ和復振幅稀疏B的最小值。于是,較低子帶數(shù)據(jù)可以由式(15)給出的互相干數(shù)據(jù)代替:

(15)

這樣就完成了各個子帶的互相干。在完成各個子帶互相干之后,可以進行多子帶融合步驟,各子帶可以看作根據(jù)目標散射點格式由多個正弦波疊加的混合信號,如圖3所示。

圖3 多子帶信號示意圖Fig.3 Schematic diagram of multi-subband signal

根據(jù)Cuomo等[6]提出的全極點模型,提出一種計算更加簡便的方式:將多個子帶的Hankel矩陣疊加,根據(jù)一個Hankel矩陣計算出各個子帶需要預測頻段的參數(shù),根據(jù)將要進行預測的所用子帶的Hankel矩陣構(gòu)建整體前向預測矩陣:

(16)

根據(jù)整體前向預測矩陣可以直接求解得到多子帶融合信號的極點值和極點幅值,相較于兩兩子帶構(gòu)造預測矩陣可以減少很大的運算量。按照式(9)進行奇異值分解,得到U、S、V矩陣。然后按照步驟1～3進行計算,得到整體極點值與極點幅值系數(shù)。

根據(jù)式(7)即可求出預測子帶的值,得到全頻段信號子帶表示如式(17)所示:

(17)

其中,Ni表示各子帶頻點數(shù),Nmi表示待預測子帶頻點數(shù)。這樣就可以得到一個由多子帶融合而成的全頻段子帶。

3 實例仿真

以三段子帶信號為例﹐其起始頻率分別為f1、f2、f3,頻率步進數(shù)分別為N1、N2、N3,高低子帶的跳頻間隔均為Δf,其中f2=f1+ΔB1,f3=f2+ΔB2, ΔB>Δf·N1,ΔB為高子帶與低子帶起始頻率間隔。根據(jù)式(1)以兩個散射點目標為例,按表1的各子帶數(shù)據(jù)參數(shù)進行實驗得到各子帶信號s1、s2、s3。

表1 各子帶數(shù)據(jù)參數(shù)

圖4 降噪前主對角奇異值Fig.4 Singular value of main diagonal before noise reduction

當信噪比越低時,奇異值曲線逐漸趨于平滑,對模型精度會有一定影響。對主對角奇異值矩陣按第2節(jié)中方法對SNR=-15 dB和SNR=-20 dB的信號進行噪聲抑制處理,將噪聲分量歸置為零或是極小值,對目標主分量乘以權(quán)重。得到噪聲抑制后的對角奇異值如圖5所示。

圖5 噪聲抑制后的主對角奇異值Fig.5 Singular value of main diagonal after noise suppression

噪聲抑制后的SNR=-15 dB與SNR=-20 dB信號主對角奇異值與SNR=10 dB基本重合,且噪聲分量低于SNR=-10 dB。根據(jù)式(5)重構(gòu)Hankel矩陣,然后根據(jù)式(6)重構(gòu)信號。對重構(gòu)后的信號進行傅里葉變換得到距離包絡,如圖6所示。根據(jù)圖6可知在距離方向上噪聲抑制明顯,達到10 dB以上。

(a) SNR=-15 dB噪聲抑制后信號距離包絡(a) SNR=-15 dB signal distance envelope after noise suppression

為驗證算法可靠性,選取SNR=-15 dB進行100次蒙特卡羅實驗,并計算均方根誤差(root mean square error, RMSE),結(jié)果如圖7所示。從圖7(a)中看出,以信噪比為10 dB的信號脈壓值為真值,100次實驗均方根誤差在0～10 dB范圍內(nèi)浮動,即可認為該算法在噪聲抑制方面具有可靠性。為驗證在信噪比范圍內(nèi)算法的可行性,選取信噪比范圍-20～10 dB,對每個信噪比取值都進行100次蒙特卡羅實驗,并以信噪比為10 dB的信號脈壓值為真值,計算100次實驗的RMSE平均值。由圖7(b)可知,各信噪比在噪聲抑制處理后相較于10 dB信噪比都有所減小,且與信噪比成反比,即算法在各信噪比情況下都有適用性和可靠性。

(a) 100次實驗的均方根誤差(a) Root mean square error of 100 experiments

按照第2節(jié)中步驟1和步驟2進行計算得到各子帶極點值,可以觀察噪聲抑制前的極點分布,如圖8所示。在沒有抑制噪聲時,低信噪比情況下,極點值不能準確估計。

圖8 噪聲抑制前子帶1極點分布Fig.8 Subband 1 pole distribution before noise suppression

各子帶未相干前極點分布并不能重合,在單位圓上存在角度偏差,在距離包絡也有一定差距,如圖9(a)、(c)所示,值與Δθ1與Δθ2一致。完成極點值計算后,按照第2節(jié)中步驟3進行極點幅值計算,即可根據(jù)式(7)得到子帶估計值,根據(jù)相干函數(shù)即可估計極點旋轉(zhuǎn)角度和復振幅系數(shù)值。根據(jù)式(15)完成各子帶相干,得到相干后的極點分布如圖9(b)、(d)所示。對比各子帶距離包絡,未相干時各子帶距離包絡無法對齊,相干后各子帶距離包絡保持一致。

(a) 噪聲抑制后各子帶極點分布(未相干)(a) Distribution of poles of each sub-band after noise suppression (not coherent)

完成相干后即可對各子帶根據(jù)式(16)構(gòu)建整體全極點模型前向預測矩陣,再根據(jù)第2節(jié)中步驟1、 2、 3進行計算,得到整體極點值與極點幅值系數(shù)。然后根據(jù)式(7),按照式(17)子帶形式進行外推,即可得到全頻段信號。若不進行相干處理,兩子帶外推信號將在頻域上不重合,融合信號會在距離包絡上產(chǎn)生分量。未相干直接對子帶1和子帶2進行外推,結(jié)果如圖10(a)所示;相干后子帶1與子帶2的外推信號在頻域上重合,不會產(chǎn)生多余的距離分量,如圖10(b)所示。

(a) 未相干子帶1與子帶2外推(a) Extrapolation of uncorrelated subband 1 and subband 2

最終可以獲得一個合成寬帶信號,當由多個子帶融合而成時,信號帶寬變大。根據(jù)距離分辨單元公式(式中c為光速,B為信號帶寬)

(18)

信號帶寬與距離分辨單元成反比,帶寬越大,能分辨的最小距離單元越小,距離分辨率越高。

多頻段融合合成寬帶信號如圖11所示。由圖可知,合成信號與原始信號之間的誤差很小,這樣拓展帶寬后,一些不能分辨的距離也能分辨出。

圖11 多頻段融合合成寬帶信號Fig.11 Multi-band fusion synthesis of broadband signals

通過蒙特卡羅實驗,對-30～20 dB范圍內(nèi)信噪比進行100次實驗并取無噪聲信號為真值計算RMSE的平均值。與不進行噪聲抑制直接使用算法進行子帶融合相比,進行噪聲抑制后,子帶模型更加精確,如圖12所示。

圖12 合成寬帶信號蒙特卡羅實驗RMSEFig.12 Monte Carlo experiment RMSE for synthesizing broadband signal

4 結(jié)論

通過上述實驗結(jié)果可知,通過在全極點模型中對Hankel矩陣進行奇異值分解這一運算過程中,加入噪聲抑制處理,對主對角奇異值進行加權(quán)、對噪聲分量進行歸置,可以大大減弱噪聲在全極點模型的模型階數(shù)中對精度的影響,且能更加準確地獲得信號檢測結(jié)果,仿真實驗中對信噪比抑制能達到10 dB以上。

在抑制噪聲之后,對全極點模型各參數(shù)計算更加精確。利用整體全極點模型前向預測矩陣可以將已經(jīng)相干后的子帶直接求得多子帶融合信號的極點值與極點幅值系數(shù),減小運算重復量。