李 海, 程偉杰, 謝瑞杰

(中國(guó)民航大學(xué)天津市智能信號(hào)與圖像處理重點(diǎn)實(shí)驗(yàn)室, 天津 300300)

0 引 言

低空風(fēng)切變是一種高度在600 m以下、氣流在一個(gè)較小的距離內(nèi)突然改變其方向或速度的大氣現(xiàn)象,它具有持續(xù)時(shí)間短、作用區(qū)域小、瞬間強(qiáng)度大、危害性強(qiáng)等特點(diǎn)。當(dāng)飛機(jī)在起降階段突遇低空風(fēng)切變時(shí),飛行員往往缺乏充足的時(shí)間和空間對(duì)飛機(jī)姿態(tài)進(jìn)行控制,從而導(dǎo)致飛行事故的發(fā)生,因此低空風(fēng)切變檢測(cè)和預(yù)警成為當(dāng)前民航領(lǐng)域的一項(xiàng)重要課題,而低空風(fēng)切變風(fēng)速估計(jì)決定了低空風(fēng)切變檢測(cè)的準(zhǔn)確程度,因而準(zhǔn)確估計(jì)低空風(fēng)切變風(fēng)速顯得至關(guān)重要[1]。機(jī)載氣象雷達(dá)能夠?qū)铰飞系奈ｋU(xiǎn)氣象進(jìn)行實(shí)時(shí)探測(cè)并發(fā)出預(yù)警,已經(jīng)成為保障民航飛機(jī)飛行安全必備的航空電子設(shè)備。但其屬于下視工作模式,所面臨的雜波分布范圍廣、強(qiáng)度大,同時(shí)由于飛機(jī)運(yùn)動(dòng)致使雜波譜嚴(yán)重展寬,從而導(dǎo)致風(fēng)切變這類分布式氣象目標(biāo)常被淹沒(méi)于強(qiáng)地雜波之中,故地雜波抑制效果將直接影響風(fēng)速估計(jì)結(jié)果[2]。

在機(jī)載相控陣?yán)走_(dá)系統(tǒng)中,空時(shí)自適應(yīng)處理(space-time adaptive processing, STAP)技術(shù)能夠在空時(shí)聯(lián)合域有效抑制雜波并達(dá)到目標(biāo)檢測(cè)的目的。STAP技術(shù)的最優(yōu)權(quán)矢量依賴于雜波協(xié)方差矩陣(clutter covariance matrix, CCM)的準(zhǔn)確估計(jì)[3]。但是傳統(tǒng)STAP技術(shù)要求在獨(dú)立同分布(independent identically distributed, IID)環(huán)境下,IID訓(xùn)練樣本達(dá)到系統(tǒng)自由度(degree of freedom, DOF)2倍或以上時(shí),方可保證輸出信噪比的損失不超過(guò)3 dB[4]。然而在實(shí)際情況中,傳統(tǒng)STAP對(duì)于IID訓(xùn)練樣本往往難以達(dá)到要求,從而無(wú)法得到準(zhǔn)確的雜波協(xié)方差矩陣,此時(shí)STAP在實(shí)際環(huán)境中的性能急劇下降。尤其是在機(jī)載前視陣下,地雜波空時(shí)二維譜隨距離變化呈現(xiàn)差異性,這種雜波特性大大減少了可用于估計(jì)雜波協(xié)方差矩陣的IID訓(xùn)練樣本數(shù)量,加大了實(shí)現(xiàn)STAP的難度。因此,探討在少量鄰近距離單元的情況下進(jìn)行低空風(fēng)切變風(fēng)速估計(jì)具有意義。

目前,在小樣本的情況下實(shí)現(xiàn)分布式目標(biāo)參數(shù)估計(jì)的手段主要有降秩STAP[5]、降維STAP[6]、稀疏STAP[7]等。文獻(xiàn)[5]提出了對(duì)角加載樣本協(xié)方差矩陣求逆(loaded sampling covariance matrix inversion, LSMI)方法,該方法直接利用雷達(dá)回波數(shù)據(jù)估計(jì)濾波器權(quán)向量,能將IID訓(xùn)練樣本數(shù)目減少至雜波空間維數(shù)的兩倍多,但是仍需利用較多的IID訓(xùn)練樣本估計(jì)雜波協(xié)方差矩陣,實(shí)際應(yīng)用中難以得到滿足;文獻(xiàn)[6]提出了基于多通道聯(lián)合自適應(yīng)處理(multiple Doppler channels joint adaptive processing, M-CAP)的低空風(fēng)切變風(fēng)速估計(jì)方法,該方法采用了降維結(jié)構(gòu),將全局系統(tǒng)自由度降至局域系統(tǒng)自由度,但I(xiàn)ID訓(xùn)練樣本數(shù)量仍然較大。稀疏恢復(fù)理論可利用極少觀測(cè)樣本高精度恢復(fù)信號(hào)[8-9],將該特點(diǎn)與STAP技術(shù)相結(jié)合可提高在小樣本條件下雜波抑制與目標(biāo)檢測(cè)性能。然而傳統(tǒng)的稀疏恢復(fù)方法,如正交匹配追蹤(orthogonal matching pursuit,OMP)算法[10]、加權(quán)二范數(shù)最小化(focal underdetermined system solver,FOCUSS)算法[11]等存在參數(shù)設(shè)置的問(wèn)題,其中最重要的參數(shù)就是正則化參數(shù)和稀疏度。如果稀疏恢復(fù)算法對(duì)參數(shù)比較敏感,則參數(shù)設(shè)置會(huì)影響到其算法的性能。正則化參數(shù)是L1范數(shù)最優(yōu)化的問(wèn)題,在實(shí)際中難以設(shè)計(jì),而稀疏度盡管在正側(cè)陣下與雜波秩有關(guān),但是風(fēng)切變檢測(cè)是處于機(jī)載前視陣情況下,此種情況下稀疏度理論研究暫時(shí)沒(méi)有文獻(xiàn)報(bào)道。因此,研究參數(shù)設(shè)置簡(jiǎn)單的稀疏STAP算法,并將其應(yīng)用至低空風(fēng)切變風(fēng)速估計(jì)問(wèn)題中具有重要意義。

針對(duì)上述小樣本、稀疏恢復(fù)算法參數(shù)設(shè)置難的問(wèn)題,提出了一種基于同倫稀疏STAP的低空風(fēng)切變風(fēng)速估計(jì)方法。同倫是一種參數(shù)設(shè)置簡(jiǎn)單的非參數(shù)化稀疏恢復(fù)方法,該方法僅需要少量鄰近距離單元就能夠獲得較為精確的雜波協(xié)方差矩陣,然后構(gòu)建STAP處理器實(shí)現(xiàn)對(duì)雜波的抑制和歸一化多普勒頻率估計(jì),最終得到風(fēng)場(chǎng)速度的準(zhǔn)確估計(jì)。該方法能夠很好地提高小樣本條件下雜波協(xié)方差矩陣的估計(jì)性能,從而有效抑制雜波并對(duì)風(fēng)速進(jìn)行準(zhǔn)確估計(jì)。

1 回波信號(hào)模型描述

本文中,xl表示第l個(gè)待檢測(cè)距離單元的NK×1維空時(shí)二維快拍數(shù)據(jù),其表達(dá)式為

xl=sl+cl+nl

(1)

式中:sl為第l個(gè)待檢測(cè)距離單元內(nèi)低空風(fēng)切變風(fēng)場(chǎng)產(chǎn)生的雷達(dá)回波信號(hào);cl為第l個(gè)待檢測(cè)距離單元內(nèi)的地雜波,在此假設(shè)地雜波沒(méi)有起伏也不存在模糊現(xiàn)象[12];nl為加性高斯白噪聲。

1.1 低空風(fēng)切變信號(hào)模型

第l個(gè)待檢測(cè)距離單元內(nèi)的低空風(fēng)切變風(fēng)場(chǎng)回波信號(hào)數(shù)據(jù)sl為

(2)

(3)

(4)

1.2 地雜波信號(hào)模型

(5)

(6)

式中:α為載機(jī)飛行方向與天線陣面之間的夾角,大小為90°。

則第l個(gè)待檢測(cè)距離單元的地雜波回波信號(hào)數(shù)據(jù)cl為

(7)

(8)

2 基于同倫稀疏STAP的風(fēng)速估計(jì)方法

基于同倫稀疏STAP的風(fēng)速估計(jì)方法首先將少量鄰近距離單元進(jìn)行距離依賴性矯正得到IID訓(xùn)練樣本,然后利用同倫算法進(jìn)行雜波協(xié)方差矩陣估計(jì),接著通過(guò)構(gòu)造STAP處理器求解最優(yōu)權(quán)矢量,最后實(shí)現(xiàn)雜波抑制和低空風(fēng)切變風(fēng)速估計(jì)。其中距離依賴性矯正、基于同倫算法的雜波協(xié)方差矩陣估計(jì)和低空風(fēng)切變風(fēng)速估計(jì)是同倫稀疏STAP方法的關(guān)鍵步驟,下面分別進(jìn)行論述。

2.1 距離依賴性矯正

由于機(jī)載前視陣下雜波分布具有距離依賴性,從而導(dǎo)致鄰近距離單元的雜波空時(shí)二維分布特性與待檢測(cè)距離單元存在差異,因此需要先進(jìn)行距離依賴性矯正。本文采用多普勒頻移(Doppler warping,DW)法進(jìn)行距離依賴性矯正,該方法簡(jiǎn)單并且易于工程實(shí)現(xiàn)[14]。該方法依據(jù)雜波空時(shí)二維譜的空間角頻率和多普勒頻率之間具有耦合特性計(jì)算得到各鄰近距離單元雜波的多普勒頻率,然后將鄰近距離單元雜波進(jìn)行一維多普勒頻率平移,使得鄰近距離單元與待檢測(cè)距離單元雜波主瓣重合,從而減少IID訓(xùn)練樣本間的雜波空時(shí)二維分布的差異性,即降低雜波距離依賴性。

(9)

式中：

(10)

則經(jīng)過(guò)DW矯正后的第p個(gè)鄰近距離單元為

(11)

2.2 基于同倫算法的雜波協(xié)方差矩陣估計(jì)

雜波協(xié)方差矩陣估計(jì)需要超完備基矩陣和空時(shí)分布向量,利用同倫算法恢復(fù)在超完備基矩陣下的空時(shí)分布向量,就能夠獲得高性能的雜波協(xié)方差矩陣估計(jì)。

(1) 超完備基矩陣

(12)

(2) 空時(shí)分布向量

(13)

由于式(13)是一個(gè)欠定方程,即存在無(wú)數(shù)可能的解[15-16]。因此式(13)可以被轉(zhuǎn)化為凸優(yōu)化中的L1范數(shù)最優(yōu)化問(wèn)題[17],也稱為L(zhǎng)ASSO(least absolute shrinkage and selection operator)問(wèn)題來(lái)進(jìn)行求解:

(14)

式中:Γ(αp)是目標(biāo)函數(shù);β是正則化參數(shù)。

式(14)中的最優(yōu)化問(wèn)題可以轉(zhuǎn)化為求目標(biāo)函數(shù)Γ(αp)次微分等于零向量的解,即?Γ(αp)/?αp=0[18-19]。目標(biāo)函數(shù)Γ(αp)的次微分表示為

(15)

根據(jù)KKT(Karush-Kuhn-Tucker)條件,?Γ(αp)/?αp=0可以等價(jià)轉(zhuǎn)化[21-22]為

(16)

式中:β≥0;Λ、u分別表示αp的支撐集、αp在其支撐集Λ上的符號(hào)序列,支撐集是αp中非零元素所對(duì)應(yīng)單元網(wǎng)格索引的集合;AΛ表示由支撐集Λ對(duì)應(yīng)單元網(wǎng)格的空時(shí)導(dǎo)向矢量所組成的集合;Λc表示支撐集Λ的補(bǔ)集,補(bǔ)集是αp中零元素所對(duì)應(yīng)單元網(wǎng)格索引的集合。

同倫算法通過(guò)下式進(jìn)行迭代更新,表達(dá)式為

(17)

式中:δ為更新方向;τ為步長(zhǎng)。

(18)

(19)

(20)

(21)

(22)

(23)

則其對(duì)應(yīng)的雜波協(xié)方差矩陣可表示為

(24)

2.3 低空風(fēng)切變風(fēng)速估計(jì)

根據(jù)式(23)得出的雜波協(xié)方差矩陣來(lái)設(shè)計(jì)STAP處理器從而消除雜波并估計(jì)風(fēng)速。由線性約束最小方差準(zhǔn)則求解最優(yōu)權(quán)矢量,該問(wèn)題可表示為

(25)

利用最優(yōu)權(quán)矢量wl對(duì)雷達(dá)回波信號(hào)進(jìn)行處理后,待檢測(cè)距離單元xl中的地雜波被抑制,且該距離單元低空風(fēng)切變信號(hào)的歸一化多普勒頻率估計(jì)值可由下式搜索得到:

(26)

則第l個(gè)待檢測(cè)距離單元的風(fēng)場(chǎng)目標(biāo)速度估計(jì)結(jié)果為

(27)

3 方法流程

基于同倫稀疏STAP的低空風(fēng)切變風(fēng)速估計(jì)方法流程如圖2所示。

所提方法可以有效地抑制地雜波并進(jìn)行低空風(fēng)切變風(fēng)速估計(jì),其關(guān)鍵步驟如下:

步驟 1利用DW補(bǔ)償法矯正雜波距離依賴性得到IID訓(xùn)練樣本;

步驟 2通過(guò)同倫算法計(jì)算待檢測(cè)距離單元的空時(shí)分布向量;

步驟 3對(duì)待檢測(cè)距離單元的雜波協(xié)方差矩陣進(jìn)行估計(jì)并求解STAP處理器的最優(yōu)權(quán)矢量,接著對(duì)地雜波進(jìn)行抑制;

步驟 4估計(jì)待檢測(cè)距離單元內(nèi)低空風(fēng)切變的歸一化多普勒頻率,從而得出待檢測(cè)距離單元的風(fēng)場(chǎng)中心風(fēng)速估計(jì)值。

4 仿真結(jié)果及分析

4.1 仿真條件設(shè)置

本文中假設(shè)低空風(fēng)切變場(chǎng)中心位于載機(jī)平臺(tái)60°方向距離12.5 km處,天線主瓣對(duì)準(zhǔn)低空風(fēng)切變場(chǎng),低空風(fēng)切變場(chǎng)寬度為8 km,高度為600 m,風(fēng)場(chǎng)水平風(fēng)速范圍為-40～40 m/s,垂直風(fēng)速范圍為-20～20 m/s。主要仿真參數(shù)如表1所示。

表1 系統(tǒng)仿真參數(shù)Table 1 System simulation parameters

4.2 仿真結(jié)果分析

圖3為機(jī)載氣象雷達(dá)前視陣下雷達(dá)回波信號(hào)的空時(shí)二維譜,由于機(jī)載前視陣下地雜波具有空時(shí)耦合特性,因此地雜波在空時(shí)二維譜中呈現(xiàn)半圓形特征,低空風(fēng)切變信號(hào)功率譜在歸一化多普勒頻率域和空間錐角余弦域具有一定的展寬。通過(guò)分析可得,地雜波的功率比低空風(fēng)切變信號(hào)的功率要強(qiáng)得多,使得地雜波完全掩蓋低空風(fēng)切變信號(hào)的歸一化多普勒頻率信息,這是造成低空風(fēng)切變估計(jì)性能不準(zhǔn)確的重要因素。

圖4為同倫稀疏STAP算法恢復(fù)地雜波空時(shí)二維譜,其中圖4(a)表示雷達(dá)回波信號(hào)中的地雜波空時(shí)二維譜,圖4(b)為直接使用1個(gè)IID訓(xùn)練樣本得出的空時(shí)二維譜,圖4(c)～圖4(f)為在IID訓(xùn)練樣本數(shù)量為1、5、10和15時(shí),利用同倫稀疏STAP方法處理得到的空時(shí)二維譜。從中可以看出所提方法可以僅使用少量IID訓(xùn)練樣本數(shù)據(jù)就能夠高分辨地恢復(fù)出地雜波空時(shí)二維譜;隨著IID訓(xùn)練樣本數(shù)量的增加,所提方法性能并沒(méi)有較大改善,說(shuō)明所提方法對(duì)IID訓(xùn)練樣本的依賴性比較低。

圖5為所提方法、LSMI、M-CAP和OMP方法雜波協(xié)方差矩陣的雜波特征譜圖。從中可以看出,所提方法和OMP方法同屬于稀疏恢復(fù)類SATP方法,其雜波協(xié)方差矩陣大特征值個(gè)數(shù)明顯要比LSMI和M-CAP兩種非稀疏恢復(fù)類STAP方法少,即雜波自由度低,更有利于STAP算法性能的提高;同時(shí)所提方法的大特征值個(gè)數(shù)比OMP方法少,說(shuō)明所提方法得到的雜波協(xié)方差矩陣更加準(zhǔn)確。圖6為所提方法、LSMI、M-CAP和OMP方法在IID訓(xùn)練樣本數(shù)量為10時(shí)的改善因子曲線圖,以第60號(hào)距離單元為例,從中可以看出所提方法相比于其他3種方法,改善因子曲線在主雜波區(qū)凹口更深和更窄,可以更為有效抑制地雜波。

本文方法、LSMI、M-CAP和OMP方法在不同IID訓(xùn)練樣本數(shù)量時(shí)的風(fēng)速估計(jì)結(jié)果對(duì)比如圖7所示。圖7(a)為IID樣本數(shù)為100時(shí)風(fēng)速估計(jì)結(jié)果,4種方法由于樣本數(shù)量充足均能夠有效估計(jì)風(fēng)速;圖7(b)為IID樣本數(shù)為50時(shí)風(fēng)速估計(jì)結(jié)果,LSMI方法由于樣本數(shù)減少無(wú)法得到準(zhǔn)確的風(fēng)速估計(jì)結(jié)果;圖7(c)為IID樣本數(shù)為10時(shí)風(fēng)速估計(jì)結(jié)果,從中可以看出,所提方法在IID訓(xùn)練樣本數(shù)量極少時(shí)依然可以準(zhǔn)確估計(jì)風(fēng)速且精度更高,同時(shí)在8.5～16.5 km范圍內(nèi),風(fēng)場(chǎng)風(fēng)速隨距離呈現(xiàn)反“S”型變化特征。而LSMI和M-CAP方法由于IID訓(xùn)練樣本嚴(yán)重不足導(dǎo)致無(wú)法估計(jì)風(fēng)速,OMP方法雖能大致估計(jì)風(fēng)速但存在誤差精度較低。

表2為不同樣本數(shù)的風(fēng)速估計(jì)均方根誤差對(duì)比?？梢钥闯?隨著IID樣本數(shù)的減少,LSMI、M-CAP和OMP方法的均方根誤差變化大,而本文方法的均方根誤差變化小,且均方根誤差值都小于其他3種對(duì)比方法。

表2 不同樣本數(shù)的風(fēng)速估計(jì)均方根誤差對(duì)比Table 2 Comparison of root mean square error of wind speed estimation with different number of samples

圖8為IID樣本數(shù)為100時(shí)不同雜波強(qiáng)度下4種方法的風(fēng)速估計(jì)結(jié)果對(duì)比圖。圖8(a)～圖8(c)為雜噪比分別為40 dB、50 dB和60 dB時(shí)的風(fēng)速估計(jì)結(jié)果,從中可以看出所提方法在雜噪比增加時(shí)依然能夠準(zhǔn)確估計(jì)風(fēng)速。

表3為不同雜波強(qiáng)度的風(fēng)速估計(jì)均方根誤差對(duì)比。從表中可以看出,雜噪比的增加對(duì)LSMI、M-CAP和OMP方法的均方根誤差值影響大,而本文方法的均方根誤差值受雜噪比影響較小。

表3 不同雜波強(qiáng)度的風(fēng)速估計(jì)均方根誤差對(duì)比Table 3 Comparison of root mean square error of wind speed estimation with different clutter intensities

5 結(jié) 論

本文將同倫稀疏恢復(fù)算法應(yīng)用到機(jī)載氣象雷達(dá)STAP處理中,并針對(duì)在機(jī)載氣象雷達(dá)前視陣下由于IID訓(xùn)練樣本不足而導(dǎo)致STAP性能下降的問(wèn)題,提出了一種基于同倫稀疏STAP的低空風(fēng)切變風(fēng)速估計(jì)的方法。該方法將待檢測(cè)距離單元周圍的少量鄰近距離單元經(jīng)過(guò)距離依賴性矯正后作為IID訓(xùn)練樣本,然后對(duì)其進(jìn)行同倫稀疏恢復(fù)處理得到高分辨率雜波空時(shí)二維譜,進(jìn)而計(jì)算相應(yīng)的雜波協(xié)方差矩陣,通過(guò)構(gòu)造STAP處理器對(duì)雜波進(jìn)行抑制并對(duì)低空風(fēng)切變風(fēng)速進(jìn)行估計(jì)。仿真結(jié)果表明,本文所提方法在IID訓(xùn)練樣本數(shù)量極少情況下依然可以準(zhǔn)確估計(jì)風(fēng)速,可顯著降低對(duì)IID訓(xùn)練樣本的需求。