復雜地形環境下基于GAMP-STAP的低空風切變風速估計方法

李 海 謝瑞杰 謝伶莉 孟凡旺

①(中國民航大學天津市智能信號與圖像處理重點實驗室 天津 300300)

②(中國航空工業集團公司雷華電子技術研究所 無錫 214063)

1 引言

低空風切變是影響民航運輸安全和效率的危險天氣之一，其具有持續時間短、作用距離小、瞬時強度大、變化多端、危害性強等典型特征[1]。當飛機在起降過程遭遇風切變時，由于此時飛行高度低，飛行員缺乏足夠的時間和空間調整飛機姿態，從而導致飛行事故的發生，因而低空風切變檢測和預警成為民航安全領域的重要課題。而風切變風速估計是風切變檢測和預警的基礎，因此，準確估計風速至關重要。近年來，我國大力推進機場建設，中西部高原地區需建設的機場越來越多[2]。建設在高原地區的機場，由于海拔較高，機場接近大氣層的中層，亂流多，風場變化大，低空風切變明顯，同時復雜的地形環境增大風切變的探測難度，嚴重威脅航空安全。機載氣象雷達作為飛行氣象環境實時探測的重要手段，探測低空風切變時，處于下視工作模式，風切變目標淹沒在強地雜波中，因此首先需要雜波抑制然后估計風速。而高原機場復雜的地形環境，增大了地雜波的非均勻性，也會使風切變的類型變為更加復雜的非對稱型風切變，進一步增大雜波抑制風速估計的難度，從而影響風切變檢測的可靠性。

對于相控陣體制的機載氣象雷達，空時自適應處理(Space-Time Adaptive Processing, STAP)技術[3]可以在空時聯合域中獲得更好的地雜波抑制效果。STAP技術最優權矢量依賴于雜波協方差矩陣(Clutter Covariance Catrix, CCM)的準確估計[4]。但是傳統STAP技術要求獨立同分布(Independent Identically Distributed, IID)樣本達到系統自由度(Degree Of Freedom, DOF) 2倍或以上[5]才能準確的估計雜波協方差矩陣，以構造STAP濾波器。而在復雜地形環境下，由于地雜波的非均勻性和IID樣本數嚴重匱乏，無法滿足該條件。因此，探討復雜地形環境情況下的低空風切變風速估計具有十分重要的意義。

目前，對于非均勻環境下的STAP目標檢測技術主要有非均勻STAP[6]、降秩STAP[7]、降維STAP[8]、知識輔助STAP(Knowledge-Assisted STAP, KA-STAP)[9,10]、稀疏STAP[11]。單樣本直接數據域(Direct Data Domain, DDD)非均勻STAP方法有效地消除了雜波非均勻的影響，但在空時域滑窗得到IID樣本會存在孔徑損失問題且受空域和時域自由度的影響較大，影響雜波抑制的性能。對角加載樣本協方差矩陣求逆(L o a d e d Sampling Covariance Matrix Inversion, LSMI)的降秩STAP和多通道聯合自適應處理(Multiple Doppler Channels Joint Adaptive Processing, MCAP)的降維STAP方法都減少了雜波協方差矩陣估計的IID樣本數，但是，所需樣本的數量仍然較多。知識輔助STAP (Knowledge-Assisted STAP,KA-STAP)方法利用地面雜波先驗信息，減少估計雜波協方差矩陣的樣本數需求，能夠很好地估計出風速。然而復雜環境的無法預料可能導致IID訓練樣本嚴重不足，并且無法獲取外界地理環境知識的輔助，無法滿足KA-STAP的條件。稀疏STAP方法利用雜波空時功率譜的稀疏特性可在極少樣本情況下高精度地恢復出雜波協方差矩陣，因而利用稀疏恢復STAP可提高復雜非均勻地形環境下雜波抑制與目標檢測性能。然而傳統稀疏恢復類算法[12]，如加權二范數最小化(FOCal Underdetermined System Solver, FOCUSS)算法[13]、正交匹配追蹤(Orthogonal Matching Pursuit, OMP)算法[14]等方法面臨計算復雜度高且正則化參數設置難的問題；稀疏貝葉斯學習方法[15](Sparse Bayesian Learning,SBL)，無須設置正則化參數，性能穩健，但是SBL每次迭代都伴隨著矩陣逆運算，因此隨著問題維度的增加，將導致很高的計算量，從而無法進行實時處理。

針對上述稀疏恢復算法參數設置困難和運算量大的不足，本文提出了一種基于廣義近似消息傳遞(Generalized Approximate Message Passing,GAMP)STAP的方法。GAMP是近似消息傳遞算法[16]的推廣，稀疏恢復問題中，GAMP能夠有效近似計算邊緣后驗分布，將其運用于稀疏恢復算法，可以減少算法復雜度[17]。因此本文利用GAMP算法開發低復雜度稀疏貝葉斯學習算法，該方法采用高斯混合稀疏先驗模型，將GAMP嵌入到EM(Expectation-Maximization)框架中，通過GAMP迭代得到稀疏系數矢量后驗分布的近似[18]。從而規避了傳統SBL算法中矩陣求逆運算，顯著地降低了運算時間。此外，引入了飛機飛行參數和雷達系統參數輔助GAMP-STAP中空時導向字典的構造，以避免網格失配和字典間相關性太強進而導致稀疏恢復算法性能下降的問題。該方法實現了復雜地形環境下雜波協方差矩陣的準確估計，有效地提高了STAP抑制雜波估計風速的性能。

2 機載氣象雷達回波信號模型

機載雷達前視陣由N元均勻線陣組成，其探測低空風切變時如圖1，機載平臺高度為H，速度為V，θ0,φl分別表示風切變目標的方位角和俯仰角，θ和φ分別表示地面散射體的方位角和俯仰角。相干脈沖數為K，脈沖重復頻率為fr，相干處理時間有L個距離單元，第l(1,2,...,L)個距離單元空時快拍數據yl表示為

圖1 機載氣象雷達視探測低空風切變幾何模型

其中，nl為 高斯白噪聲，sl為 第l個距離單元的風切變信號，cl為非均勻地雜波信號，本文采用Ward模型融合DEM和NLCD數據仿真高保真地雜波信號[19]。

sl表示為

其中，Γ為風切變信號復幅度，A(ψ0，fl)為風切變信號空時導向矢量，ψ0為風切變信號空間錐角，cosψ0=cosθ0cosφl，fl(fl ∈[?1,1])為風切變信號歸一化多普勒頻率，As(θ0，?l)、At(fl)分別為風切變信號的空域導向矢量和時域導向矢量[20]

3 基于GAMP-STAP的低空風切變風速估計方法

基于GAMP-STAP的低空風切變風速估計方法結合飛機飛行參數和雷達系統參數，利用雜波脊的先驗信息構造空時導向字典；然后在空時導向字典上通過GAMP算法迭代獲得雜波幅度計算協方差矩陣；進而設計STAP濾波器實現雜波抑制并估計風速，該方法原理框圖如圖2所示。

圖2 GAMP-STAP的低空風切變風速估計原理框圖

3.1 空時導向字典構造

機載雷達回波信號由不同方向和不同多普勒頻率的回波信號疊加而成。如果將空域和多普勒頻率分別均分成Ns=ρsN和Nt=ρtK個網格點(ρs，ρt分別表示離散化系數)，整個空時平面被分成J=Ns×Nt個 點，則距離單元l的回波接收數據稀疏模型[22]可以表示為

目前，稀疏恢復STAP方法中字典主要是根據以上均勻劃分處理得到，雖然簡單易處理，但是存在網格失配問題。因此本文利用雜波脊的先驗信息構造空時導向字典[23]可以減小上述影響。

機載前視陣雷達下雜波的歸一化多普勒頻率與歸一化空間頻率關系為

3.2 基于GAMP估計雜波協方差矩陣

基于GAMP算法估計雜波協方差矩陣，首先建立GAMP模型，給出稀疏系數矢量的稀疏先驗分布，通過GAMP迭代得到稀疏系數矢量的近似真實后驗分布，然后利用EM算法更新超參數，最后獲得空時導向字典每個網格點的雜波幅度以計算雜波協方差矩陣。

(1) GAMP模型。GAMP算法中采用通用的D項高斯混合(GM)[23]對待估計的稀疏向量x=[x1,x2,...,xg,...,xG]T建模，其先驗分布為式(8)所示的復高斯分布：

根據文獻[24]，GAMP算法中有兩個重要的近似：

一是使用輸出端的邊緣后驗分布來近似真正的邊緣后驗分布p(zw|y;q)( 其中w=1,2,...,W)

其中，pY|Z(yw|zw;q)是 輸出端的似然函數，zw?aTwx表示無噪聲輸出量，aTw表 示字典A的 第w行 ，?pw和μzw是GAMP算法的迭代中間變量。

另一個是在輸入端用下式來近似真實的后驗分布pX|Y(xg|y;q)( 其中g=1,2,...,G)

此時，根據貝葉斯準則，結合稀疏向量先驗式(8)和后驗分布式(10)，得到稀疏向量的邊緣分布，根據高斯分布乘積法則可得到pX|Y(xg|y;q)服從高斯分布，具體式(11)所示。

GAMP算法估計待重建信號的過程可用圖3的因子圖來表示，其中輸出變量xg(紅色圓圈表示)代表空時平面上第g(g=1,2,...,G)個網格點復幅度，即所求的稀疏向量；輸出變量yw(藍色方框表示)代表接收數據向量分量，即雷達回波數據IID樣本；首先是紅色的小箭頭，表示信息傳遞方向，代表的是xg傳 遞給yw節 點的信息，即xg結合自身先驗概率更新的后驗概率；藍色的小箭頭表示yw傳 遞給xg節點的信息，即對xg的后驗概率的更新推斷；物理意義即在已知雷達接收數據、稀疏字典以及輸入后驗分布、輸出端的先驗分布基礎上，通過輸入端和輸出端的不停迭代直到算法收斂[24]，得到基于MMSE準則下xg的估計。對于GAMP算法，最重要的是尺度函數的求解。

圖3 GAMP算法估計待重建信號因子圖

3.3 風速估計

對雷達回波信號進行濾波處理后，當前待檢測距離單元的地雜波被抑制，且該距離單元風場目標的多普勒頻率可由下式搜索得到

4 方法流程

本文提出的基于GAMP-STAP的風切變風速估計方法流程如圖4所示。

圖4 基于GAMP-STAP的風切變風速估計方法流程圖

步驟總結如下：

步驟1 初始化參數值；

步驟2 結合飛機飛行參數和雷達系統參數構造空時導向字典；

步驟3 給出稀疏系數向量的稀疏先驗模型，并利用GAMP算法進行迭代，求解稀疏系數向量；

步驟4 將GAMP算法得到的稀疏系數向量和中間量按照式(16)、式(17)EM更新規則，進行迭代更新，并判斷是否達到迭代停止條件，若無，則重復步驟3，若有，則輸出最終的稀疏系數向量估計值；

步驟5 計算雜波協方差矩陣；構建STAP濾波器，進行雜波抑制以及低空風切變風速估計。

5 實驗結果分析

5.1 仿真條件描述

本文中，設定非對稱低空風切變場的入射傾角為20°，入口速度為40 m/s ，風場范圍為飛機前方8.5～16.5 km，雷達系統仿真參數如表1所示。

表1 雷達系統仿真參數

5.2 仿真結果分析

圖5展示了地雜波的距離多普勒圖，圖5(a)和圖5(b)分別為均勻環境和復雜地形環境下地雜波距離多普勒圖，從圖中可看出，復雜地形環境下雜波的非均勻性和非平穩性現象十分明顯，表現了高原地區回波數據特征。

圖5 地雜波距離多普勒圖

圖6為前視陣機載氣象雷達回波信號空時2維譜，雜波譜呈半圓形分布，體現了雜波譜前視陣結構的空時耦合特性，低空風切變信號在空時域呈現為有一定的展寬的一條“窄帶”，并且從圖中可以看出地雜波信號的回波功率遠大于低空風切變信號的回波功率，風切變回波多普勒信息淹沒在地雜波的多普勒信息中，嚴重影響風切變速度估計與檢測。

圖6 機載氣象雷達回波空時2維譜

圖7為本文基于GAMP算法對第60號距離單元雜波功率譜的恢復結果，空間尺度ρs=10,ρt=8，當采用傳統的空時字典時，算法性能不穩定并不能保證每次算法都能收斂；當網格劃分致密導致字典的列相關性太強，算法性能損失時，雜波譜恢復結果如圖7(a)所示，圖7(b)為在相同的空間尺度下引入雜波脊先驗信息構造的字典情況下恢復的結果，能夠保證算法每次都能收斂，得到較好的信號恢復結果。同時說明本文算法利用8個訓練樣本能夠較好地恢復雜波的功率譜。

圖7 GAMP算法恢復的雜波功率譜

圖8為GAMP-STAP方法與傳統STAP, KASTAP, SBL-STAP的風速估計結果對比圖，可以看出本文方法在復雜地形環境下能夠準確估計風速，而傳統STAP, KA-STAP方法性能很差。相較于SBL-STAP方法，本文方法準確度更高一些，并且運算量和運算速度得到很大提升。

圖8 不同方法風速估計結果對比

表2為本文方法與傳統SBL-STAP方法完成單個距離門風速估計的運行時間對比，其中W=NK為接收信號維度，其中N為陣元數，K為脈沖數；G為空時2維平面的離散點數。計算復雜度為1次迭代復雜度，運行時間為估計1個距離單元風速的時間。取同樣的空間分辨尺度(ρs=10,ρt=8)，可以看出由于傳統SBL-STAP方法涉及矩陣逆運算，計算復雜度高，運行時間慢，而本文方法在貝葉斯框架下使用GAMP算法來迭代估計雜波幅度，計算過程僅涉及乘法運算，大大減少了運行時間。

表2 算法運行時間對比

6 結束語

本文針對復雜地形環境造成的非均勻雜波情況和IID樣本數不足問題，提出一種基于GAMP-STAP的風切變風速估計方法。在稀疏背景下，該方法首先利用雜波脊的先驗信息構造空時導向字典，接著用GAMP算法迭代計算稀疏系數向量，并用EM算法更新學習超參數集合，最后得到雜波協方差矩陣的估計結果，構造STAP濾波器抑制雜波實現風切變風速的準確估計。實驗仿真結果表明該方法能使用少量IID樣本實現風速較準確估計，并且算法計算復雜度低，較于傳統SBL-STAP大大減小了運行時間，對于復雜地形環境下的低空風切變風速估計具有更好的實用價值。