李 海，唐 芳，李雙雙

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

0 引 言

在航空氣象領(lǐng)域，通常將飛機(jī)高度在600 m 以下，風(fēng)向風(fēng)速急劇變化的這一大氣現(xiàn)象稱為低空風(fēng)切變，其具有變化快、強(qiáng)度大、危害性高等典型特征［1］。如果在起飛或降落階段飛機(jī)遭遇風(fēng)切變，而此時飛行員缺乏充足的時間來調(diào)整飛機(jī)狀態(tài)，可能會導(dǎo)致嚴(yán)重飛行事故發(fā)生［2］。因此，有關(guān)風(fēng)切變檢測的問題一直以來都是民航領(lǐng)域關(guān)注熱點(diǎn)。作為對于低空風(fēng)切變檢測流程的關(guān)鍵，風(fēng)場速度估計的結(jié)果將直接影響整個風(fēng)切變檢測結(jié)果的準(zhǔn)確程度，因此，準(zhǔn)確估計風(fēng)速至關(guān)重要［3-5］。

機(jī)載氣象雷達(dá)是飛行中不可或缺的設(shè)備，可以對飛機(jī)航路上氣象狀況進(jìn)行實(shí)時監(jiān)測，及時獲取氣象數(shù)據(jù)發(fā)現(xiàn)危險天氣，以幫助飛行員提前規(guī)避風(fēng)切變等危險大氣現(xiàn)象，從而保障飛行過程的安全［6］。隨著航空領(lǐng)域事業(yè)的發(fā)展和進(jìn)步，航空安全的實(shí)際使用需求對機(jī)載氣象雷達(dá)的探測范圍、掃描精度等性能提出了更高的要求。目前，與相對成熟的平面相控陣?yán)走_(dá)相比，安裝在飛機(jī)表面與飛行器貼合的共形陣機(jī)載氣象雷達(dá)具有許多優(yōu)點(diǎn)，如共形陣天線孔徑的增大能提高目標(biāo)探測時方位向的分辨率，天線掃描范圍的擴(kuò)大有助于實(shí)現(xiàn)全方位掃描等［7］。共形陣機(jī)載雷達(dá)所具備的優(yōu)點(diǎn)可以更好地適應(yīng)航空氣象領(lǐng)域發(fā)展的需求。因此，對共形陣機(jī)載氣象雷達(dá)體制下目標(biāo)檢測的課題意義重大。

機(jī)載氣象雷達(dá)在檢測風(fēng)切變過程中風(fēng)速估計是關(guān)鍵步驟，但由于地雜波分布范圍較廣，強(qiáng)度大，導(dǎo)致雷達(dá)接收到的風(fēng)切變信號完全被地雜波信號所覆蓋，進(jìn)而導(dǎo)致風(fēng)速估計結(jié)果的準(zhǔn)確性受到影響。目前，針對低空風(fēng)切變風(fēng)速估計方法主要基于色加載知識輔助的空時自適應(yīng)處理（Space-Time Adaptive Processing, STAP）風(fēng)速估計方法［8］、基于直接數(shù)據(jù)域-擴(kuò)展因子化法（Direct Data Domain-Extended Factored Approach, DDD-EFA）的風(fēng)速估計方法［9］、基于同倫稀疏STAP 的風(fēng)速估計方法［10］和基于改進(jìn)輔助通道的風(fēng)速估計方法［11］等。以上方法在傳統(tǒng)線陣或面陣體制下都能準(zhǔn)確估計風(fēng)速，但對于陣元在空間上呈現(xiàn)兩維或三維分布的半球共形陣，還未有文獻(xiàn)對其風(fēng)切變檢測技術(shù)和風(fēng)速估計方法進(jìn)行研究報道。

基于此，本文提出了一種基于D3AR 風(fēng)速估計方法。該方法將空時自回歸算法（Space-Time Autoregression，STAR）思想引入到直接數(shù)據(jù)域（Direct Data Domain, DDD）算法中，首先將待檢測單元的數(shù)據(jù)分別從空域、時域以及空時域進(jìn)行信號對消處理，然后將處理后的數(shù)據(jù)矩陣描述為空時自回歸模型并估計模型參數(shù)，再通過構(gòu)造與雜波子空間正交的空間來實(shí)現(xiàn)對雜波的抑制，最后通過提取待檢測單元的最大多普勒頻率來估計風(fēng)場速度。根據(jù)仿真結(jié)果顯示，該方法有效地實(shí)現(xiàn)了地雜波抑制，并且能夠精確估計風(fēng)速。

1 半球共形陣?yán)走_(dá)信號模型

1.1 地雜波信號

基于Ward 模型［12］仿真半球共形陣體制下的雜波回波信號。首先依據(jù)載機(jī)雷達(dá)參數(shù)確認(rèn)雷達(dá)所照射區(qū)域內(nèi)單個距離單元的寬度，然后將單個距離單元按照方位角劃分為多個雜波塊并計算每個雜波塊的回波信號，通過將單個距離單元內(nèi)的所有雜波塊相加，得到雜波回波信號。

圖1為半球共形陣機(jī)載氣象雷達(dá)幾何模型，假設(shè)半球共形陣共有M個半徑分別為Rm的圓?。旤c(diǎn)不放陣元），每一層圓弧上均勻分布著N個陣元，且相鄰兩個圓弧之間的間距為d，θ為方位角，φ為俯仰角，ψ為空間錐角，V為載機(jī)速度，由幾何關(guān)系可知波束矢量K(θ,φ)=[cosθcosφsinθcosφsinφ]T。

圖1 半球共形陣機(jī)載氣象雷達(dá)幾何模型

如圖1（b）所示，令右側(cè)第一段圓弧為半球形陣列的第一層，每一層的第一個陣元位于X軸正方向，則由半球共形陣機(jī)載雷達(dá)陣列幾何模型可以推導(dǎo)出第m(m= 1,2,…,M)行第n(n= 1,2,…,N)列的陣元位置坐標(biāo)Pm,n為

式中，Rm為第m層圓弧的半徑，，β為共形陣機(jī)載雷達(dá)圓弧所對應(yīng)的圓心角，α為圖1（c）YOZ平面所示的夾角。

根據(jù)陣元坐標(biāo)以及波束矢量推導(dǎo)得到半球共形陣?yán)走_(dá)地雜波回波信號的空間角頻率wsc：

式中，φl為第l(l= 1,2,…,L)個距離單元雜波回波的俯仰角，L為回波總距離門數(shù)。

半球共形陣?yán)走_(dá)地雜波回波信號的時間角頻率wtc為

式中，fr為脈沖重復(fù)頻率，λ為波長。

假設(shè)半球形陣列天線的第m行第n列個陣元在第k個脈沖下接收的第l個距離單元的雜波回波數(shù)據(jù)為Cl(m,n,k)，則

式中，Rl為第l個距離單元的斜距，θ的取值范圍為[ 0,π ]。

1.2 低空風(fēng)切變信號

由于低空風(fēng)切變是體分布式目標(biāo)，所包含的氣象粒子具有微粒性、疊加性和隨機(jī)性等特點(diǎn)，因此采用撒點(diǎn)法［13］對風(fēng)切變場回波數(shù)據(jù)進(jìn)行仿真，通過將單個距離單元內(nèi)的所有散射點(diǎn)相加，得到該距離單元的風(fēng)場回波信號。則半球形陣列天線的第m行第n列陣元在第k個脈沖下接收的第l個距離單元的風(fēng)切變信號回波數(shù)據(jù)為Sl(m,n,k)，且有

式中，Q表示一個距離單元內(nèi)散射點(diǎn)總數(shù)，Rq為第q個散射點(diǎn)的斜距，Aq為由雷達(dá)方程推導(dǎo)得到的第q個散射點(diǎn)回波的幅度［14］：

式中，Z為風(fēng)場反射率因子，Pt為機(jī)載雷達(dá)的發(fā)射功率，G為天線的增益。

式（5）中，wss(θq,φq)表示第q個散射點(diǎn)的空間角頻率，wts(θq,φq)為時間角頻率，且有

式中，θq為第q個散射點(diǎn)的方位角，φq為俯仰角，vq為徑向速度。

1.3 回波信號

當(dāng)半球共形陣機(jī)載氣象雷達(dá)用于探測風(fēng)切變時，其回波信號由地雜波C、風(fēng)切變S和噪聲組成，則半球形陣列天線接收的第l個距離單元的總回波數(shù)據(jù)Xl為

假設(shè)在雷達(dá)工作時，在一個CPI 內(nèi)發(fā)射K個脈沖，用Xl表示第l個距離單元的半球共形陣機(jī)載氣象雷達(dá)所接收的回波數(shù)據(jù)矩陣：

2 基于D3AR 低空風(fēng)切變風(fēng)速估計方法

基于D3AR 低空風(fēng)切變風(fēng)速估計方法首先對待檢測單元的回波數(shù)據(jù)分別從空域、時域、空時域進(jìn)行信號對消處理，然后將處理后的樣本數(shù)據(jù)通過AR 模型構(gòu)造與雜波子空間正交的空間來實(shí)現(xiàn)對雜波的抑制，最后通過提取待檢測單元的最大多普勒頻率來估計風(fēng)場速度。下面詳細(xì)討論以下3個部分：信號對消處理、估計STAR 模型系數(shù)矩陣以及低空風(fēng)切變風(fēng)速估計。

2.1 信號對消處理

DDD 方法僅利用待檢測距離單元的數(shù)據(jù)信息進(jìn)行樣本數(shù)據(jù)的獲取，當(dāng)待檢測距離單元含有低空風(fēng)切變信號目標(biāo)時，此時通過AR 模型來構(gòu)造的空時濾波器在抑制雜波的同時會造成風(fēng)切變信號損失。因此，為了防止風(fēng)切變信號被抑制導(dǎo)致無法準(zhǔn)確估計風(fēng)速，首先需要對待檢測距離單元數(shù)據(jù)分別從時域、空域以及空時域進(jìn)行信號對消處理，具體過程如下：

定義風(fēng)切變信號Sl(m,n,k) 為

根據(jù)載機(jī)速度、脈沖重復(fù)頻率等先驗(yàn)信息，假設(shè)風(fēng)切變信號方位已知，fl為歸一化多普勒頻率（fl的范圍為[]-1,1）。此時，對于低空風(fēng)切變信號，相鄰兩脈沖間的相位差為

相鄰兩陣元之間的相位差為

式中，ΔPm,n為相鄰兩陣元位置坐標(biāo)差。因此待檢測單元可通過利用兩陣元（脈沖）之間的相位差分別從空域、時域以及空時域進(jìn)行信號對消處理。

對第l個待檢測單元的和回波數(shù)據(jù)從空域利用兩陣元相消濾除風(fēng)切變信號時，可以得到如下數(shù)據(jù)矩陣：

式中，Yls為第l個待檢測距離單元從空域利用兩陣元相消濾除風(fēng)切變信號后的數(shù)據(jù)矩陣，其維度為(MN- 1) ×K，“*”表示共軛運(yùn)算。

對第l個待檢測距離單元的回波數(shù)據(jù)從時域進(jìn)行對消處理時，當(dāng)待更新的fl在[]-1,1 中取某一值時，利用兩脈沖對消后得到如下數(shù)據(jù)矩陣：

式中，Ylt為第l個待檢測距離單元從時域利用兩脈沖相消進(jìn)行信號濾除處理后的數(shù)據(jù)矩陣，其維度為MN×(K- 1) 。

同理，利用兩脈沖、兩陣元相消對第l個待檢測距離單元回波數(shù)據(jù)從空時域進(jìn)行信號對消處理后的數(shù)據(jù)矩陣為

式中，Ylst為第l個待檢測距離單元從空時域利用兩陣元兩脈沖相消進(jìn)行信號濾除處理后的數(shù)據(jù)矩陣，其維度為(MN- 1) ×(K- 1) 。

更新fl的取值并執(zhí)行式（13）～（15），當(dāng)待檢測距離單元利用信號對消后的樣本數(shù)據(jù)矩陣經(jīng)過AR 模型所構(gòu)造的空時自適應(yīng)處理器的輸出功率最大時，此時Yls、Ylt和Ylst3 個數(shù)據(jù)矩陣中的風(fēng)切變信號被有效消除。

2.2 估計STAR模型系數(shù)矩陣

假設(shè)當(dāng)待更新的fl在[-1,1 ]范圍中所取某值是該待檢測距離單元風(fēng)場目標(biāo)的歸一化多普勒頻率，此時由該單元通過信號對消處理后得到的3個樣本數(shù)據(jù)矩陣Yls、Ylt以及Ylst中就只包含雜波和噪聲，為了方便計算，選取Yls中的(K- 1) 列數(shù)據(jù)、以及Ylt中的(MN- 1) 行數(shù)據(jù)，再將3 個樣本數(shù)據(jù)矩陣排列成(MN- 1 )(K- 1) × 1 維矢量并分別記為Y1、Y2和Y3。由于任意隨機(jī)過程都可描述為AR 模型，則Y1、Y2和Y3滿足如下自回歸模型［15］：

式中，t= 1,2,3，k= 1,2,…,K-P，Ai為D×(MN- 1)維的空時自回歸系數(shù)矩陣，D代表AR 模型的空域階數(shù)，P為時域階數(shù)。AR 模型的時域階數(shù)和空域階數(shù)未知，需要通過樣本數(shù)據(jù)對其進(jìn)行估計，由于已有許多文獻(xiàn)介紹了AR 模型階數(shù)估計方法［15-16］，本文不再敘述。模型階數(shù)D、P確定后，再通過訓(xùn)練樣本估計自回歸系數(shù)矩陣。

記A=[A0A1…AD-1]，A的維數(shù)為D×(MN-1)P，將處理后的訓(xùn)練樣本單元的數(shù)據(jù)重新改寫為

式中，et為(MN- 1)P×(K-P)維矩陣，則處理后的所有樣本數(shù)據(jù)矩陣為

式中，E為(MN- 1)P× 3(K-P)維矩陣，因此，根據(jù)式（17），訓(xùn)練樣本數(shù)據(jù)矩陣E滿足如下方程：

可以看出，式（20）中A的解是由矩陣E的零空間確定，但由于E滿秩不存在零空間，所以根據(jù)最小均方誤差準(zhǔn)則，對于自回歸系數(shù)矩陣A的求解可以轉(zhuǎn)換為最小二乘準(zhǔn)則來估計［16］：

為了確保不產(chǎn)生平凡解，A必須滿足AAH=I,且樣本矩陣E的D個小奇異值對應(yīng)的左奇異矢量構(gòu)成矩陣A的列。

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

完成自回歸系數(shù)矩陣A的估計后，利用A0,A1,…,AD-1構(gòu)造如下濾波器系數(shù)矩陣：

式中，B為D(K-P)×(MN- 1 )(K- 1) 維矩陣,且滿足如下方程：

矩陣B與雜波和噪聲所張成的子空間正交，根據(jù)匹配濾波理論，得到D3AR 算法的空時自適應(yīng)權(quán)矢量為［17］

式中，ws為歸一化空間角頻率，fl上文中假設(shè)的歸一化多普勒頻率。由此可以得到信號匹配結(jié)果為

通過更新fl的取值，可以得到對應(yīng)的空時自適應(yīng)處理器權(quán)矢量，進(jìn)行信號匹配后求解輸出功率。當(dāng)其最大時，此時fl的取值即為對應(yīng)距離單元的多普勒頻率估計值。這是因?yàn)楫?dāng)處理器的輸出功率最大化時，該待檢測距離單元經(jīng)過信號對消處理后的樣本數(shù)據(jù)矩陣中風(fēng)切變信號被有效濾除，再通過AR 模型所構(gòu)造的空時濾波器就可以使雜波信號得到有效抑制，同時風(fēng)切變信號得到積累；如果處理器的輸出功率不是最大，說明該待檢測距離單元在信號對消處理過程中風(fēng)切變信號沒有被有效濾除，導(dǎo)致得到的樣本數(shù)據(jù)經(jīng)過AR 模型所構(gòu)造的空時濾波器在抑制雜波的同時會造成風(fēng)切變信號損失。即可以通過式（28）來估計得到該待檢測距離單元：

最后，依照上述方法分別估計每一個距離單元的風(fēng)速值，完成整個風(fēng)場的風(fēng)速估計。

2.4 方法流程及步驟

圖2展示了基于D3AR的半球共形陣風(fēng)切變風(fēng)速估計方法的整體流程。

圖2 基于D3AR風(fēng)速估計方法基本流程

步驟1：分別從空域、時域以及空時域利用兩脈沖（陣元）進(jìn)行信號對消處理；

步驟2：將處理后的訓(xùn)練樣本單元描述為空時自回歸模型；

步驟3：根據(jù)樣本單元估計空時自回歸模型階數(shù)，并構(gòu)造濾波器模型系數(shù)矩陣；

步驟4：根據(jù)匹配濾波理論，構(gòu)造與雜波子空間正交的空間來實(shí)現(xiàn)對雜波的抑制，最后通過搜索該待檢測距離單元的輸出功率最大時所對應(yīng)的多普勒頻率，完成風(fēng)速估計；

步驟5：更新待檢測單元，重復(fù)上述步驟估計所有待檢測距離單元的風(fēng)速值，完成整個風(fēng)場速度估計。

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

3.1 仿真參數(shù)設(shè)計

假設(shè)風(fēng)場目標(biāo)處于8.5～16.5 km 處，表1 所示為仿真參數(shù)值。

表1 載機(jī)及雷達(dá)仿真參數(shù)

3.2 仿真結(jié)果分析

圖3 為半球共形陣機(jī)載氣象雷達(dá)回波信號第87 號距離門空時譜仿真結(jié)果，由圖可以看出，地雜波信號的空時二維譜在空間的分布表現(xiàn)為橢圓形狀，風(fēng)切變信號在空間表現(xiàn)為一條窄帶，且地雜波信號的功率明顯大于風(fēng)切變信號，使得風(fēng)切變信號完全被淹沒，從而導(dǎo)致難以被檢測。

圖3 第87號距離單元空時譜

圖4 為D3AR 方法、最優(yōu)STAP 方法、DW-STAP方法以及STAR 方法的第87 號距離單元的改善因子對比圖。從圖可以看出，上述方法均在主雜波區(qū)域即零頻區(qū)形成凹口，造成性能損失，與其余方法相比，本文方法在主雜波區(qū)域形成的凹口更窄、更淺，雜波抑制性能明顯提高。

圖4 改善因子對比圖

圖5 為D3AR 方法與最優(yōu)STAP 方法、DWSTAP 方法、以及STAR 方法風(fēng)速估計結(jié)果對比圖。從圖中可以看出，本文所提方法在不利用參考單元樣本數(shù)據(jù)的情況下仍然能較準(zhǔn)確地估計風(fēng)速。對比圖顯示，最優(yōu)STAP 法的風(fēng)速估計結(jié)果不理想，主要是由于半球共形陣特殊陣列流型導(dǎo)致雜波非平穩(wěn)性更強(qiáng)，使得所需的獨(dú)立同分布樣本數(shù)不滿足RMB 準(zhǔn)則，最優(yōu)STAP 方法性能受損［5］；經(jīng)DW補(bǔ)償后，仍然不能得到理想效果，這是由于DW方法僅對時域進(jìn)行了多普勒頻率補(bǔ)償，沒有對空間頻率進(jìn)行補(bǔ)償，所以經(jīng)過DW 補(bǔ)償后，仍不能準(zhǔn)確估計風(fēng)速［18］；STAR 方法需要用到與待檢測單元相鄰單元的數(shù)據(jù)來構(gòu)造AR 模型建立空時二維濾波器進(jìn)行雜波抑制以及估計風(fēng)速［17］，其結(jié)果的準(zhǔn)確性依賴于IID 樣本數(shù)量。而D3AR 方法僅利用待檢測單元的數(shù)據(jù)且不需要進(jìn)行距離依賴性矯正就能準(zhǔn)確估計風(fēng)速。

圖5 風(fēng)速估計結(jié)果

本文所提方法主要包括求濾波器權(quán)矢量和空時濾波兩部分的運(yùn)算量，表2 為D3AR 方法、最優(yōu)STAP 方法、DW-STAP 方法以及STAR 方法運(yùn)算量對比，其中L1為最優(yōu)STAP 方法和DW-STAP 方法中所需的訓(xùn)練樣本數(shù)。從表中可以看出本文所提方法的運(yùn)算量比最優(yōu)STAP方法、DW-STAP法以及STAR法都要小。

表2 不同風(fēng)速估計方法運(yùn)算量對比

表3 給出了D3AR 方法、最優(yōu)STAP 方法、DWSTAP 方法以及STAR 方法的均方根誤差對比，可以看出，本文所提D3AR方法的均方根誤差最小。

表3 不同風(fēng)速估計方法均方根誤差

4 結(jié)束語

針對半球共形陣體制下進(jìn)行低空風(fēng)切變檢測時會受到強(qiáng)低雜波信號的干擾，導(dǎo)致風(fēng)切變信號難以檢測的問題，提出了一種基于D3AR 的半球共形陣風(fēng)切變風(fēng)速估計方法。該方法將空時自回歸算法思想引入到直接數(shù)據(jù)算法思想中，僅用待檢測單元濾除風(fēng)切變信號后的樣本數(shù)據(jù)來構(gòu)建空時自回歸模型，再通過構(gòu)造與雜波子空間正交的空間來實(shí)現(xiàn)對雜波的抑制，最后通過提取待檢測單元的最大多普勒頻率來估計風(fēng)場速度。仿真結(jié)果表明，該方法可以在半球共形陣體制下僅用待檢測距離單元的數(shù)據(jù)準(zhǔn)確估計風(fēng)速，且該方法僅在時域進(jìn)行滑窗處理，避免了空域滑窗對陣列結(jié)構(gòu)的嚴(yán)格要求，因此適用于任意陣列，但由于該方法中AR 模型參數(shù)的選擇會影響風(fēng)速估計結(jié)果的準(zhǔn)確性，因此需進(jìn)一步研究更為準(zhǔn)確的AR 模型參數(shù)估計算法。