李 海，湛 蕾，吳仁彪

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

大氣的運動形式千變?nèi)f化，對航空運輸?shù)陌踩?、正點及效益有著不容忽視的影響。大氣湍流是飛機在飛行過程中時常會遇到的一種危險氣象目標，它是由大氣中急速變換的氣流引起的，具有隨機性、間歇性和不規(guī)則性。湍流區(qū)域內(nèi)氣流運動速度向量有很大差異，并且變化都相當急劇。因此飛機一旦飛進湍流區(qū)域，輕則會產(chǎn)生顛簸，重則可能會短暫失控，發(fā)生難以想象的后果。根據(jù)民航組織（ICAO）的飛行事故和事故癥候報告統(tǒng)計，每9次天氣因素造成的飛行事故中，就有1次是由大氣湍流造成的[1]。所以，準確檢測出大氣湍流并有效避開湍流航線成為近半個世紀以來航空界關注的重要研究課題之一[2]。

機載氣象雷達具有湍流檢測、雨區(qū)探測和風切變預警等功能，可為飛行員實時提供危害天氣的方位及強度，極大地提高了飛機飛行的安全性[3]。對于機載氣象雷達而言，大氣湍流的運動特征體現(xiàn)在風場粒子的風速分布特性上，對湍流的檢測關鍵在于如何從回波中獲取風場的譜矩信息，這些信息包含功率、平均速度以及譜寬，其中譜寬是實現(xiàn)湍流目標檢測最為重要的參數(shù)。大量實際檢測數(shù)據(jù)表明，湍流回波的多普勒頻譜具有較大的譜寬，譜寬越大，湍流越強[1]。根據(jù)國際民航界統(tǒng)一標準，規(guī)定速度偏差即譜寬大于5 m/s的氣象回波目標為湍流目標[4]。

目前，大氣湍流的譜寬估計普遍采用脈沖對處理法（PPP）。傳統(tǒng)的脈沖對處理法通過處理信號回波采樣序列的互相關函數(shù)來對湍流回波的平均速度和譜寬進行估計，該方法運算復雜度低，能夠滿足機載氣象雷達的實時性探測，但其受噪聲的影響較大，僅適用于回波信噪比高的情況[5]?？焖俑道锶~變換（FFT）估計法也是湍流回波譜寬估計的常用方法之一。文獻[6-7]分析了基于FFT的湍流譜寬估計方法及其改進方法，將時域的回波信號轉(zhuǎn)換到頻域?qū)ψV寬進行估計，但是FFT法受回波數(shù)據(jù)長度的限制，在脈沖數(shù)較少或信噪比較低的情況下，譜寬估計性能變差。對于含水量低的大氣湍流，雷達反射率弱，導致回波信噪低，此時傳統(tǒng)的PPP法和FFT法都無法有效估計湍流回波譜寬。針對低信噪比下傳統(tǒng)方法估計性能下降的問題，可借用參數(shù)化模型進行檢測。文獻[8]提出了一種單天線體制下基于RELAX的參數(shù)化譜寬估計方法，該方法利用非線性最小二乘準則將回波信號的自相關序列和自相關序列的理想模型進行擬合，然后通過RELAX算法獲得譜寬的準確估計。然而該算法在脈沖數(shù)少的情況下估計性能會變差。文獻[9]提出了一種基于回波幅度的譜寬估計方法，通過構建AR模型的自適應檢測器對譜寬進行估計，該算法在低信噪比下譜寬估計精度高；然而其需要極其大量的數(shù)據(jù)樣本，并不適用于機載氣象雷達的實際應用。文獻[10-11]提出一種將多普勒頻譜分析與雷達極化信息處理技術相結合的方法，對湍流目標進行分類檢測，但該方法沒有對湍流回波進行詳細的參數(shù)估計研究。

本文提出一種基于空時協(xié)方差矩陣擬合的湍流譜寬估計方法，從大氣湍流的分布式氣象特性出發(fā)，利用湍流信號服從高斯分布的特點，給出湍流目標協(xié)方差矩陣參數(shù)化模型。然后分析回波信號協(xié)方差矩陣與湍流目標協(xié)方差矩陣之間的關系，從而構造相應的代價函數(shù)。最后通過對代價函數(shù)的求解實現(xiàn)湍流譜寬的精確估計。

1 湍流模型

1.1 回波模型

設載機以勻速V直線飛行，機載雷達天線為N元前視線陣，陣元間距為d（d=λ/2，λ為雷達脈沖波長）。圖1為湍流幾何觀測圖，圖中X-Y平面為湍流散射體所在平面。假設1個距離單元的寬度為ΔR，其中包含Q 個湍流散射粒子，對于第 q（q=1，2，…，Q）個散射粒子，其相對于飛機的徑向距離為Rq，方位角、俯仰角以及空間錐角分別為 θq、φq和 ψq，且有cos θqcos φq=cos ψq。

圖1 大氣湍流觀測圖Fig.1 Geometry of turbulence observation

一個相干處理間隔（CPI，coherent processing interval）內(nèi)雷達發(fā)射的脈沖數(shù)為K，則對第l（l=1，2，…，L，L為距離單元個數(shù)）個待檢測距離單元，雷達接收的空時快拍數(shù)據(jù)可表示為

其中：xl（n，k）表示雷達天線第n（n=1，2，…，N）個陣元發(fā)射的第 k（k=1，2，…，K）個脈沖對該距離單元所接收的采樣數(shù)據(jù)；sl表示該距離單元的湍流目標信號；nl表示加性高斯白噪聲。由于機載氣象雷達在檢測湍流時巡航高度較高，故在此不考慮地雜波對雷達回波的影響。

大量實際檢測數(shù)據(jù)表明，大氣湍流是一種高斯型的隨機過程，即風場中湍流散射粒子的速度大小服從正態(tài)分布，在此基礎上對湍流風場建模[12]，得到大氣湍流的速度場，用sl（n，k）表示第n個陣元的第k個脈沖接收的第l個待檢測距離單元內(nèi)湍流目標的數(shù)據(jù)，則有

其中：Aq為第q個湍流散射粒子的回波幅度；Rq表示該散射粒子相對于飛機的距離；vq為第q個湍流散射粒子與飛機的相對速度；ωs（θq，φq）和ωt（vq）分別為第q個湍流散射粒子的空間角頻率和時間角頻率，且有

1.2 湍流目標協(xié)方差矩陣模型

本文算法需要求解湍流目標的協(xié)方差矩陣模型，在此利用湍流信號的分布特性對其建立空時二維協(xié)方差矩陣參數(shù)化模型。由于大氣湍流屬于分布式氣象目標，其功率譜服從高斯譜，據(jù)此推導出湍流目標協(xié)方差矩陣Rs的參數(shù)化模型[12]為

其中：P為回波信號功率；σf為信號多普勒譜寬；f0為信號多普勒中心頻率，且有

其中：Bt為K×K維矩陣，且有Bs為N×N維全1矩陣；ψ為雷達對湍流場進行探測時天線的空間錐角。

2 基于協(xié)方差矩陣擬合的湍流譜寬估計方法

第l個待檢測距離單元回波信號的協(xié)方差矩陣可由該距離單元的回波數(shù)據(jù)樣本xl估計[13]，根據(jù)式（1）可知，一個CPI內(nèi)該距離單元的采樣數(shù)據(jù)為M（M=N×K）個，求解序列樣本 xl的自相關函數(shù) r（m1-m2）=E{xl（m1）xl（m2）}，m1，m2=1，…，M，則第l個待檢測距離單元回波信號的協(xié)方差矩陣可表示為

其中：M為數(shù)據(jù)樣本個數(shù)。

其中：λm和em為協(xié)方差矩陣的特征值和與特征值對應的特征向量；Λs=diag（λ1，λ2，…，λr）為協(xié)方差矩陣的大特征值，Es=（e1，e2，…，er）是與大特征值對應的特征向量張成的子空間（r為大特征值個數(shù)），即信號子空間；Λn=diag（λr+1，λr+2，…，λM）為小特征值，En=（er+1，er+2，…，eM）是由小特征值對應的特征向量張成的噪聲子空間。其中為湍流目標協(xié)方差矩陣，為噪聲協(xié)方差矩陣。

根據(jù)信號子空間與噪聲子空間的正交性，有

則可推出

式（10）兩邊同時右乘湍流目標協(xié)方差矩陣的參數(shù)化模型Rs（f0，σf），得

因此可利用式（11）定義如下代價函數(shù)，求解出湍流回波的多普勒譜寬f和多普勒中心頻率為

3 算法步驟

在相控陣雷達體制下，構造符合湍流特性的湍流目標空時協(xié)方差矩陣模型，與回波協(xié)方差矩陣進行擬合，實現(xiàn)湍流譜寬的準確估計。圖2為實現(xiàn)所提方法的流程圖，具體步驟描述如下。

步驟1設置載機速度V、飛機飛行高度H以及天線掃描的相關參數(shù)，初始化所有待檢測單元的雷達回波數(shù)據(jù)；

步驟2由陣列接收的待檢測距離單元回波數(shù)據(jù)估計出回波的協(xié)方差矩陣；

步驟3確定待檢測距離單元空間錐角，構造該距離單元湍流目標的協(xié)方差矩陣模型Rs（f0，σf）；

步驟4構造代價函數(shù)，在目標參數(shù)范圍內(nèi)進行二維搜索找到使代價函數(shù)最大值的參數(shù)組合由所求多普勒譜寬得到該距離單元湍流目標的譜寬估計值；

步驟5重復上述步驟，檢測下一個距離單元。

圖2 本文方法流程圖Fig.2 Flowchart of proposed method

4 仿真結果及分析

仿真條件描述：飛機速度V為200 m/s，飛行高度H為8 000 m；湍流場分布于飛機前方8.5～21.5 km處；信噪比為5 dB，天線掃描的主瓣波束等雷達參數(shù)設置如表1所示。

表1 雷達參數(shù)Tab.1 Radar simulation parameters

圖3給出了第78號距離單元湍流回波的參數(shù)估計圖。圖中Z軸為式（11）代價函數(shù)值，峰值點所對應的歸一化多普勒譜寬和多普勒中心頻率即為該距離單元回波譜矩的估計值，估計出的歸一化多普勒譜寬為根據(jù)速度譜寬與歸一化多普勒譜寬之間的關系推算出速度譜寬估計值m/s，與設定速度譜寬接近，可見在低信噪比SNR=5 dB的條件下，參數(shù)估計值較準確。

圖3 多普勒譜寬估計（二維搜索）Fig.3 Doppler spectral width estimation（2D search）

由于多普勒譜寬與多普勒平均頻率非耦合[12]，基于這一事實，對多普勒中心頻率和多普勒譜寬的二維搜索轉(zhuǎn)化為兩個一維搜索。在求解多普勒譜寬時，可先假設多普勒譜寬為某個大于0的常數(shù)σf=C（C為大于0的常數(shù)），對代價函數(shù)進行一維搜索估計出多普勒中心頻率得到后，再估計多普勒譜寬來降低運算量。圖4為第78號距離單元湍流回波的一維搜索估計圖，歸一化多普勒估計結果與二維搜索一致，說明快速算法在保持估計精度的基礎上降低了運算量。

圖4 多普勒譜寬估計（一維搜索）Fig.4 Doppler spectral width estimation（1D search）

為分析算法性能，將本文算法與文獻[8]所提基于RELAX的譜寬估計方法以及傳統(tǒng)脈沖對法在同等條件下進行了對比試驗。圖5為300次蒙特卡洛試驗條件下，3種算法在不同信噪比下譜寬估計均方誤差的對比。如圖5所示，在低信噪比條件下，本文提出的協(xié)方差矩陣擬合方法明顯優(yōu)于另兩種方法；信噪比SNR越低，優(yōu)勢越明顯。當SNR低于10 dB時，脈沖對法的譜寬估計性能急速下降，誤差較大；而本文方法依舊能夠保持較高精度。仿真結果證明，基于協(xié)方差矩陣擬合的湍流譜寬估計方法可在低信噪比下較準確地估計湍流譜寬。

圖5 譜寬估計結果對比Fig.5 Comparison of spectral width estimation results

5 結語

針對低信噪比下大氣湍流的檢測估計，提出一種基于空時協(xié)方差矩陣擬合的機載氣象雷達湍流回波譜寬估計方法。該方法是一種參數(shù)化的譜寬估計方法，利用信號協(xié)方差矩陣與回波協(xié)方差矩陣噪聲子空間的正交性將回波協(xié)方差矩陣與湍流目標協(xié)方差矩陣的參數(shù)化模型進行擬合，構造代價函數(shù)，通過對代價函數(shù)的求解來實現(xiàn)回波譜寬的精確估計。仿真實驗表明，相較傳統(tǒng)的譜寬估計方法，協(xié)方差矩陣擬合方法在低信噪比情況下，能夠獲得精確的風速譜寬估計。

[1]肖業(yè)倫,金長江.大氣擾動中的飛行原理[M].北京:國防工業(yè)出版社,1993:1-50.

[2]BORDEN S C,MATHEWS B D,STEVENSON J.Turbulence Radar System:US,5164731[P].1992.

[3]GEORGEWS.機載雷達導論[M].2nded.北京:電子工業(yè)出版社,2005:1-106.

[4]盧曉光,夏 冬.基于統(tǒng)計置信度的湍流檢測門限確定方法[J].中國民航大學學報,2011,29(4):27-30.

[5]BAXA E G,LEE J.The Pulse-pair Algorithm as a Robust Estimator of Turbulent Weather Spectral Parameters Using Airborne Pulse Doppler Radar[R].Final Report Clemson Univ.sc.radar Systems Lab,1991.

[6]DIAZ J M B,LEITAO J M N.Nonparametric estimation of mean Doppler and spectral width[J].IEEE Trans on Geosci Remote Sensing,2000,38(1):271-282.

[7]于瑩潔,李 勇.一種機載氣象雷達湍流信號的檢測方法及仿真[J].西北工業(yè)大學學報,2015,33(1):159-163.

[8]盧曉光,吳仁彪.基于RELAX的氣象雷達回波參數(shù)化譜矩估計方法[J].信號處理,2011,27(8):1250-1253.

[9]PINSKY M,VENTURA J F I,OTTO T,et al.Application of a simple adaptive estimator for an atmospheric Doppler radar[J].IEEE Transactions on Geoscience&Remote Sensing,2011,49(1):115-127.

[10]OSTROVSKY Y P,YANOVSKY F J,ROHLING H.Turbulence and precipitation classification based on Doppler-polarimetric radar data[J].IEEE Radar Symposium,2007,3(3):57-60.

[11]AVERYANOVA Y.Statistical algorithm for turbulence detection using polarization features of radar reflections from rain[J].IEEE Radar Symposium,2015,16(6):593-596.

[12]蔣 婷.基于空時處理的湍流檢測技術研究[D].天津:中國民航大學,2016.

[13]王永良,陳 輝,彭應寧,等.空間譜估計理論與算法[M].北京:清華大學出版社,2004:152-181.

[14]ZOUBIR A,WANG Y,CHARGE P.Efficient subspace-based estimator for localization of multiple incoherently distributed sources[J].IEEE Transactions on Signal Processing,2008,56(2):532-542.