基于小波變換模極大值的LiDAR風切變預警算法

張千千，史緯恒，伍 波*，萬家碩，成家豪，龔 靖，趙青虎

(1.成都信息工程大學 光電工程學院，成都 610225；2.中國氣象局 大氣探測重點開放實驗室，成都 610225；3.中國人民解放軍 32368部隊，北京 100042)

引 言

低空風切變指的是在距地面600m范圍內，風矢量(風速、風向)在水平或垂直距離上發(fā)生明顯改變的現(xiàn)象[1]，具有時間短、類型多變、破壞力強等特點，且常伴隨著有極端天氣產生，為飛行器的起飛和降落帶來了極大的困擾[2]。根據(jù)國際民航組織(International Civil Organization, ICAO)的規(guī)定：風切變強度分為輕度、中度、重度、嚴重4個等級，可用風切變因子來表示。風切變因子是指空間兩點距離(水平距離或垂直距離)為30m時風矢量差的大小[3]。由于大部分飛行事故都發(fā)生在飛機起飛和進近階段，因此，機場的風場探測技術尤為重要。相干激光測風雷達作為一種新型的探測裝置，體積小、重量輕、抗干擾能力強、分辨率高，能夠彌補其它探測手段的不足，尤其是在晴空條件下最有效的風場探測手段，對保障飛機起降安全、預警低空風切變有重要指導作用，因此，利用激光測風雷達，準確地對低空風場預測是保證飛行的重要工具。

20世紀70年代，著名氣象學家FUJITA在調查了3次飛機飛行事故后，認為始作俑者為強烈的下沉氣流，從此開啟風切變的研究[4]。1983年，WOODFIELD和WOODS提出的S因子算法[5]，在某些特殊情況下僅通過最大最小值之間的差值會產生較大誤差。2011年，香港天文臺CHEN等人提出香港機場沿用至今的單斜坡低空風切變自動預警算法[6]和一種新的下滑道掃描策略[7]，但是該算法由于受到斜坡長度的限制僅對小尺度風切變有更好的效果。2012年，CHEN提出F因子算法[8]，但直接計算由激光雷達得到的逆風梯度數(shù)據(jù)可能會導致快速波動且無法計算垂直分量的風切變。2014年，他們又提出利用渦流耗散率[9]預警低空風切變，但因為渦流耗散率的閾值在國際中沒有明確的規(guī)定，該方法還需要在實際應用中繼續(xù)探索。2016年，JIANG等人提出單雙斜坡結合的算法[10]來預警風切變，不足的是該算法只能預警大尺度風切變。2019年，MA等人提出基于主成分分析(principal component analysis,PCA)和相位差校正法的風切變預警算法[11]，但是該算法存在一定的問題且結果不夠準確，多普勒頻移和風切變閾值還需要進一步通過經驗確定。

鑒于現(xiàn)有算法存在不足，本文中提出了一種基于小波變換模極大值的激光雷達風切變預警算法。該算法通過選取一定的小波函數(shù)確定重組逆風廓線上的模極大值，再根據(jù)模極大值的閾值確定是否有風切變發(fā)生，最終確定風切變發(fā)生的位置和時刻。

1 算法原理

通過檢測小波變換模極大值來檢測信號奇異點的方法最早是由MALLAT等人提出[12]，后人對此加以引申并應用到了電力系統(tǒng)故障檢測[13]、心電圖異常信號檢測[14-15]、橋梁裂縫檢測[16]。在風切變的檢測中，風切變可以看成是徑向速度信號的突變信號。信號發(fā)生的突變時刻被認為是信號的奇異點，而小波變換模極大值通常正好對應著信號的突變點，因此，基于小波變換模極大值檢測方法可用于低空風切變預警，通過檢測速度徑向數(shù)據(jù)的小波變換模極大值確定信號的突變位置，即檢測飛機下滑道上風切變發(fā)生的時刻和位置。

MALLAT系統(tǒng)地論述了如何利用小波變換的局部化特性檢測信號奇異點位置[12]，并對小波變換的定義采用了具有濾波意義的卷積形式[17]。

Wψ(a,x)=f(x)*ψa(x)=

(1)

式中，τ為時間。

(2)

(3)

根據(jù)Fourier變換的微分性質，ζ(x)，η(x)均滿足容許性條件，均可作為小波母函數(shù)。對其分別做小波變換，有：

(4)

(5)

式中，Wζ(a,x)是f(x)通過濾波器θa(x)濾波后的1階導數(shù)。由于θa(x)是一個平滑濾波函數(shù)，所以f(x)經過θa(x)濾波后，f(x)的噪聲得到了抑制；而1階導數(shù)，即微分運算，反映了f(x)的變化率，當存在突變點時，它的變化率就很大，達到模極大值，所以Wζ(a,x)取極值點的地方就是f(x)的邊沿位置。Wη(a,x)是通過濾波器θa(x)濾波后的2階導數(shù)，2階導數(shù)的過零點對應1階導數(shù)的極值點，所以也常用Wη(a,x)的過零點來檢測信號的突變點，但是對于受到強噪聲污染的信號，Wη(a,x)的過零點很多，由此很難真正的判斷信號的邊沿。此外，Wη(a,x)的過零點只能給出拐點的位置信息，不能給出變換的快慢，所以本文中用Wζ(a,x)的模極大值來檢測徑向速度信號的風切變。若對x0的任意鄰域內的任意點x有|W(a0,x)|≤|W(a0,x0)|，則稱點(a0,x)為小波變換的模極大值。

圖1是本文中提出的新算法的流程圖。

Fig.1 Algorithm flow chart

2 算法驗證

首先討論小波變換階數(shù)及插值對檢測結果的影響，由于橫坐標表示下滑道上距離飛機降落點的距離，每個點間隔為100m，即0.1km，故插值的步長選為0.05。階數(shù)影響的結果如圖2所示，插值影響的結果如圖3所示。

圖2是在不同分解階數(shù)下使用相同小波得到的結果。其中圖2a～圖2f中的階數(shù)依次為1～6。由圖3a可知，最大的兩個速度突變點發(fā)生在x=38和x=55處，而圖2中只有圖2e即分解階數(shù)為5階時是準確的。再通過其它數(shù)據(jù)驗算，表明此算法檢測激光雷達測量反演得到的重組逆風廓線最適合的分解階數(shù)為5階，故后面的驗證將全部采用5階小波變換分解。

圖3是插值對檢測結果的影響。其中圖3a～圖3c分別為原始重組逆風廓線、首尾鏡像延長各10個點的重組逆風廓線、先插值再首尾延長各10個點的重組逆風廓線，圖3d～圖3f是分別對圖3a～圖3c做小波變換、求模極大值并歸一化的結果。根據(jù)圖3的結果，顯然，直接處理數(shù)據(jù)和首尾鏡像延長后再處理數(shù)據(jù)沒有區(qū)別，而先插值后延長的處理方式使結果有10個數(shù)據(jù)點左右的偏差，考慮到插值采用的是樣條插值，會額外數(shù)據(jù)個數(shù)，即增加了徑向速度數(shù)據(jù)中本不存在的速度值，所以導致了結果的偏移。再驗證了其它數(shù)據(jù)之后可以確定，在對激光雷達測量得到的重組逆風廓線處理時不需要插值。而對于首尾鏡像延長的步驟，是有必要的，這一步是為了防止首尾突變的風速信息缺失。經過驗證，首尾鏡像延長不會對小波變換模極大值的結果造成影響。

Fig.2 The normalized results of modulus maxima of different orders

Fig.3 Influence of interpolation on detection results

2.1 數(shù)值驗證

為了驗證上述算法的可行性和準確性，首先構造了幾組比較有代表性的徑向速度數(shù)據(jù)。數(shù)據(jù)點數(shù)均為120個，數(shù)據(jù)如圖4所示。其中圖4a為脈沖數(shù)據(jù)，圖4b為階躍數(shù)據(jù)，圖4c為斜坡數(shù)據(jù)，圖4d為明顯未發(fā)生風切變數(shù)據(jù)。4組數(shù)據(jù)的橫坐標表示下滑道上距離飛機降落點的距離，每個點間隔100m；縱坐標表示風速，單位為m/s。突變點均設置在距離下滑道上距離飛機降落點7.6km處，即橫坐標x=76。利用以下4組數(shù)據(jù)分別采用雙正交(biorthogonal， Bior)小波、多貝西(Daubechies，Db)小波對算法進行數(shù)值驗證。

經過大量小波對4組數(shù)據(jù)的驗證，表明脈沖型數(shù)據(jù)使用Bior系中雙數(shù)小波檢結果較準確；階躍型和斜坡型數(shù)據(jù)用Db系中Db5檢測結果較準確。在驗證過程中數(shù)據(jù)首尾延長各10個數(shù)據(jù)點防止信息缺失，由于本文中所說的模極大值實際意義是模的極大值和模的極小值，為了避免出現(xiàn)負值，生成的模極大值結果又經過歸一化處理，最終如圖5所示。

圖6是經過風切變判斷后保留的模極大值結果，即為風切變發(fā)生的位置。由于圖4中1～3組數(shù)據(jù)設置的突變點位置均為下滑道上距離飛機降落點7.6km處，即第76個數(shù)據(jù)點發(fā)生突變，數(shù)據(jù)首端經鏡像延長10個數(shù)據(jù)點后，應在第86個數(shù)據(jù)點發(fā)生突變。圖5展示的都是準確找到了第86個數(shù)據(jù)點的小波變換模極大值的結果，其中圖5a、圖5d中選取的小波為Bior6.8小波，圖5b、圖5c中選取的小波為Db5小波。經過風切變判斷后，圖6a～圖6c中也僅僅只保留了第86個數(shù)據(jù)點的模極大值，而圖6d中的模極大值全部置為0，預設值的結果完全一致，預警率達到100%。至此，基于小波變換模極大值的激光雷達風切變預警算法在數(shù)值驗證上通過。

2.2 外場測試

2.2.1 鄂西北鄖西測試 2017-08-08，鄂西北鄖西縣馬安鄉(xiāng)發(fā)生了局地大暴雨，使用中國氣象局的CINRAD/SA型S波段多普勒雷達對暴雨進行了監(jiān)測，該雷達一個平面位置顯示器(plane position indicator， PPI)掃描周期為6min,可以測量920個距離庫230km范圍內的天氣情況。結果顯示，20:24:00時最大風速為20.1m/s，20:30:06時最大風速為27.5m/s，20:36:13時最大風速降為19m/s，此處粗略判斷有風切變發(fā)生。

Fig.4 Four groups of reconstructed upwind profiles of the structure

圖7是該時間段內(20:24:00～20:36:13)使用天氣雷達測量得到的基本速度。由于直接截取的雷達界面圖片較模糊，因此將數(shù)據(jù)導入MATLAB并重新生成了如圖7所示的基本速度。圖7的看圖順序是從左下角到右上角，圖中各顏色色塊表示不同的風速和方向，其中白色是速度模糊的部分；左下角淺綠色到深綠色再到中間淺灰色，表明風速是逐漸減小，從-15m/s降到0m/s；右上角由深灰到橙黃色交替，表明風速逐漸增大，最大約為25m/s雷達圖的色塊反映了有極端強對流導致的風切變發(fā)生。

Fig.5 Wavelet transform modulus maximum result

與天氣雷達同時測量的激光雷達經數(shù)據(jù)反演和質量控制后得到如圖8所示的重組逆風廓線，激光雷達仰角為0.5°。由于圖8所示的重組逆風廓線前半段類似脈沖函數(shù)(脈沖較寬)，后半段類似斜坡函數(shù)，因此根據(jù)數(shù)值仿真的結果，數(shù)據(jù)點1～23選取Bior6.8小波，數(shù)據(jù)點23～120選取Db5小波，分解層數(shù)均為5階，首尾各鏡像延長10個數(shù)據(jù)點，對數(shù)據(jù)進行處理，得到如圖9所示歸一化的模極大值。圖9中前半段為Bior6.8小波處理的模極大值，后半段為Db5小波處理的模極大值，中間為了區(qū)分，用紅色虛線加以隔開。為了圖片的美觀，作者將圖9中Db5小波變換模極大值的結果整體平移到與Bior6.8小波變換模極大值的結果在同一水平高度，但是在后續(xù)閾值的處理上兩部分有不同的閾值。

Fig.7 Basic speed measured by weather radar (20:24:00～20:36:13)

Fig.8 Reorganization of the headwind profile (20:23:58～20:32:04)

Fig.9 Wavelet transform modulus maximum result

圖10是經過風切變判斷后的檢測結果。圖中，依次在1，2，3，4處有小波變換模極大值保留，由于首端有10個數(shù)據(jù)點鏡像對稱延長，且下滑道上每個數(shù)據(jù)點的距離為0.1km，因此風切變發(fā)生的實際位置為下滑道上距離飛機降落點距離1km，2km，2.6km和7.8km處，且風速差依次為-12.5m/s，18.5m/s，17.5m/s和19.5m/s，滿足兩點間距離超過30m、風速差大于7.7m/s的規(guī)定，故確認有風切變發(fā)生，風切變強度均為重度。此結果與同時測量的天氣雷達和氣象局提供的氣象資料吻合，證實由于暴雨造成強對流天氣，形成風切變，因此基于小波變換模極大值的激光雷達風切變預警算法準確。

Fig.10 Results after windshear judgment

2.2.2 攀枝花機場測試 2018年3月～11月，在青藏高原東南緣的四川省攀枝花市保安營機場削平的山頂上(經緯度為26°32′32.33″N，101°47′54.67″E)，西南技術物理研究所部署了其研制的3維激光測風雷達來檢測機場的風場情況。2018-05-06T13：30，CA4463機組報告“本場五邊有顛簸，返航成都”，部署的激光雷達正好進行下滑道掃描[18]，西南技術物理研究所提供了當時的實測下滑道側向風數(shù)據(jù)進行分析。重組逆風廓線如圖11所示。

Fig.11 The reconstructed headwind profile measured by 3-D LiDAR (13:22:29～13:31:33)

圖12是歸一化的小波變換模極大值結果。其中圖12a為Bior6.8小波變換模極大值，圖12b為Db5小波變換模極大值。圖13是經風切變判斷后保留的模極大值。其中圖13a是Bior6.8的，圖13b是Db5的。圖12中一些數(shù)值較高的模極大值點沒有保留，是因為經過風切變判斷后，這些奇異點不滿足風切變發(fā)生的速度差或者距離差，因此被置為0。對比圖13的兩幅子圖，1、2、3點的結果接近但有1～2個數(shù)據(jù)點的偏差，且速度差接近；第4個數(shù)據(jù)點，橫坐標(即發(fā)生風切變的位置)接近，但風速差差別較大；第5個數(shù)據(jù)點，兩幅子圖都是一樣的結果；第6個數(shù)據(jù)點，使用Bior6.8小波時存在，而使用Db5小波時，由于模極大值點在第91個數(shù)據(jù)點，此時徑向風速為4.9m/s，與前后均不滿足風速差大于7.7m/s的要求，因此模極大值被置為0，造成最終數(shù)據(jù)缺失。由于此數(shù)據(jù)的類型更類似于脈沖型數(shù)據(jù)，且查閱重組逆風廓線可知,Bior6.8小波的效果較Db5小波更準確。

圖13a中，有6個小波變換模極大值保留，同樣于首端有10個數(shù)據(jù)點鏡像對稱延長，因此風切變發(fā)生的實際位置為下滑道上距離飛機降落點距離3.7km，4.1km，4.8km，6.9km，8.3km和9.2km處，對應的風速差依次為-10.8m/s，11.9m/s，-18.2m/s，14m/s，-8.3m/s和9.2m/s，滿足兩點間距離超過30m、風速差大于7.7m/s的規(guī)定，故確認有風切變發(fā)生，風切變強度均為重度。對比西南技術物理研究所的激光雷達判斷的風切變和機組的報告，證明基于激光雷達的風切變預警算法是可行的且檢測結果準確。

Fig.12 Wavelet transform modulus maximum result

Fig.13 Test results after windshear judgment

3 結 論

本文中提出了一種基于小波變換模極大值的激光雷達低空風切變預警算法。該算法可以檢測各種類型的風切變，且不用考慮風切變的產生的原因和風切變的尺度,同時能夠精確檢測風切變發(fā)生的時刻和位置。通過數(shù)值驗證和外場測試證實了算法的可行性、準確性和有效性。