一種基于激光雷達探測的風場多特征反演方法

周 潔,高 航,張旭暉,李健兵

(1.國防科技大學電子科學學院,湖南 長沙 410073；2.海軍91236部隊,遼寧 葫蘆島 125100)

1 引 言

風是日常生活中最常見的一種自然現象,對人類生產生活產生巨大的影響,風中蘊含巨大的破壞能量,一旦感知不當,就將造成嚴重的人員傷亡和經濟損失。2021年3月23日,“長賜號”貨輪在蘇伊士運河河口南端6海里處疑遭瞬間強風吹襲,造成船身偏離航道,觸底擱淺。擱淺事件造成了巨大損失,全球12 %的國際貿易通道被該貨輪“切斷”,經濟損失估計超過數十億美元。2015年6月1日2130左右,重慶東方輪船公司所屬“東方之星”游輪上行至長江水域湖北荊州市監利縣大馬洲水道44號過河標水域處,遭遇嚴重垂直切變和水平切變災害性天氣,導致游輪翻沉,442人遇難[1]。2018年8月28日,首都航空JD5759航班在澳門機場著陸時出現“海豚跳”,在雙發不同程度受損、前起落架斷裂情況下復飛備降深圳,事故原因查明是突遇嚴重風切變[2]。

在眾多災害性風場氣象當中,低空風切變及湍流是國際上公認的嚴重危害飛行安全的風場現象。低空風切變通常是指近地面600 m高度以下風矢量(風向、風速)在水平和(或)垂直距離上的變化[3]；湍流是指風速的時間不規則性和空間的不均勻性由各種尺度的漩渦連續分布疊加形成[4]。目前,用于風切變和湍流探測技術主要包括測風儀、氣象雷達、風廓線雷達以及激光雷達,其中激光雷達被認為晴空條件下較優的風場探測系統。

采用激光雷達對風切變和湍流探測的掃描策略有多種,主要包括平面位置顯示器(Plan Position Indicator,PPI)、距離高度顯示器(Range Height Indicator,RHI)[5]、多普勒波束擺動(Doppler Beam Swinging,DBS)[5]、凝視[6]、滑動路徑掃描(Glide Path Scan,GPScan)[7]等方式。常見的掃描方式,只能夠實現對于單一高度角或者方位角區域風切變和湍流的探測,且存在空間覆蓋范圍和數據更新率之間的矛盾,為克服該問題,國內外開展了相關研究。2011年,陳柏緯等人基于下滑道掃描方式,利用F因子實現了對風切變的精準預測,并采用縱向結構函數實現了低空湍流預警算法,其掃描周期比傳統的PPI掃描短,從而提高了預警效率[8],但無法同時實現對于機場三維風場的探測。2015年,Sathe在常規5波束DBS掃描方式基礎上,通過最小化風場矢量測量隨機誤差矩陣構建了六波束掃描方式,并在此基礎上實現了對風場廓線和湍流參數的快速探測反演[9],但該方法采用固定高度角掃描,無法感知全域三維空間,只能提供風場廓線信息；2008年,Banakh和Smalikho基于RHI掃描方式,采用徑向速度譜反演湍流耗散率[10],2019年,采用2高度角PPI掃描方式,實現了對于非均勻湍流特征參數的探測和反演[11],該方法實現了對于各向異性湍流的探測感知,但其空間垂直采樣率不足。2017年,Norman等人位于葡萄牙開展湍流觀測實驗,利用一部采用垂直凝視方式和二部RHI掃描方式的激光雷達,利用修正的多普勒譜寬方法實現了對湍流和風切變以及低空急流的演變過程的探測[12],該實驗基于多部激光雷達協同探測,原始數據準確性提升,但經濟成本和多雷達協同方式難度增加。

綜上所述,現有各種掃描方式只能實現特定的風場特征獲取,為同時實現大范圍三維風場、湍流強度、風切變的反演,亟需提出一種兼顧空間覆蓋范圍、時空分辨率及魯棒性的激光雷達體積掃描方法。本文提出一種晴空條件下的多用途的測風激光雷達體積掃描策略及風場特征獲取方法,采用交替多個特定高度角的圓錐掃描方式,基于局部線性風場假設條件,在郎需興等人二維風場反演方法基礎上[13],通過徑向風速的處理實現三維風場的反演(后稱為3D-VPP方法)；通過計算特定高度角下飛機下滑道區域F因子強度,實現對風切變強度的評估和預警[7,14-15]；利用35.3°高度角圓錐掃描,采用部分傅里葉分解算法,快速計算湍流動能強度廓線分布情況；應用均勻湍流假設,通過計算徑向速度結構函數,建立湍流耗散率與徑向速度縱向結構函數關系式,從而計算該區域湍流耗散率分布情況[11]。通過仿真比對實驗,驗證了本方法對于機場周邊風切變和湍流的探測反演具有較好的性能。

2 方法介紹

2.1 體積掃描方式

在晴空條件下,激光雷達通過大氣氣溶膠粒子后向散射信號的多普勒頻移獲得掃描空間中風場的徑向速度。激光雷達體積掃描示意如圖1(a)所示,在多個高度角作PPI(Plan Position Indicator)掃描,該掃描方式可以獲得更大的掃描波束覆蓋范圍,可為精細的風場反演算法提供足夠的數據支撐。

圖1 激光雷達掃描策略和VPP算法分析單元的示意圖Fig.1 Schematic diagram of lidar scanning strategy and schematic diagram of VPP algorithm

激光雷達應部署在飛機降落點附近,確保激光雷達與飛機下滑道之間無障礙物遮擋,激光雷達有效探測距離Rmax能完全覆蓋飛機下滑道所有區域。進行體積掃描時交替采用多個不同高度角PPI掃描的方式,其中包含φ=1°,3°,35.3°,90°這4個特定的高度角,1°高度角主要實現對于近地面風場的探測,3°高度角主要實現對于飛機下滑道區域風切變的探測,35.3°高度角主要實現對于湍流的探測和湍流動能強度的快速計算,90°高度角作為垂直風速參考。其余高度角可按照φ=1°、3°、8°、15°、25°、35.3°、45°、90°進行配置,若特定區域需進行重點掃描,可進行加密配置,但需滿足2個相鄰高度角之差Δφ<10°,最大高度角根據所測空域高度需求及Rmax而確定(45°滿足600 m以下空域探測需求)。激光雷達旋轉角速度、掃描扇區大小可根據關注的探測空間區域而調整,在每個高度角采取完整PPI或者扇區PPI掃描方式獲得各高度角的徑向速度信息。

基于該掃描方式,在保證探測精度的同時,兼顧數據更新率,可以同時實現風場三維速度反演、下滑道的風切變識別、湍流動能強度和耗散率廓線的有效探測。

2.2 風場反演算法

在2.1節提出的體積掃描方式基礎上,采用3D-VPP方法實現對于探測空域的三維風場的反演。建立以激光雷達為原點的笛卡爾坐標系(x,y,z)和球坐標系(φ,θ,r),如圖1(a)所示,x,y,z分別為垂直機場跑道方向、平行跑道方向和垂直方向上距激光雷達的距離,φ為高度角,θ為方位角,r為到激光雷達的徑向距離。在激光雷達掃描空間中,將高度角跨度插值為φi∈[φl-Iδφ,φl+Iδφ],方位角跨度為θj∈[θm-Jδθ,θm+Jδθ],徑向距離跨度為rk∈[rn-Kδr,rn+Kδr]的連續區域作為分析體積單元Almn,其中O(φl,θm,rn)是該分析單元的中心,δφ,δθ,δr分別為激光雷達的高度角、方位角和徑向距離分辨率,lmn表示掃描空間中的第lmn個探測單元,示意圖如圖1(b)所示。2I+1,2J+1,2K+1分別為分析體積單元內高度角、方位角和徑向分辨單元的個數,φi,θj,rk分別為分析單元內任意一點的高度角、方位角、徑向距離。

設分析單元中心點速度為(ulmn,vlmn,wlmn),ulmn表示徑向速度,vlmn表示水平切向速度,wlmn表示ulmn,vlmn所在平面法向的垂直速度。對于較小尺度的分析單元,可以認為單元內所有速度相等,則分析單元內任一點P(φi,θj,rk)的徑向速度Vobs(φi,θj,rk)與中心點風速在其徑向上的風速投影F(φi,θj,rk)可以表示為式(1)。

使分析單元內觀測到的徑向風場Vobs(φi,θj,rk)與各徑向的投影風場F(φi,θj,rk)之間差的平方和最小建立目標函數式。目標函數分別關于vlmn,wlmn作變分處理,并令?H/?vlmn=0和?H/?wlmn=0,可以得到使H取極小值時的vlmn,wlmn速度的大小,ulmn,vlmn,wlmn的表達式可以表示為式(2),將反演速度轉為笛卡爾坐標系中速度(u,v,w)。

F(φi,θj,rk)=ulmn[cosφlcosφicos(θm-θj)+sinφlsinφi]-vlmncosφisin(θm-θj)+

wlmn[cosφlsinφi-sinφlcosφicos(θm-θj)]

(1)

(2)

2.3 下滑道風切變算法

假設飛機下滑道是具有一定尺度的管道區域,以3°飛機下滑道為軸選取附近30 m范圍內的風場數據,示意圖如圖2所示。激光雷達到跑道的距離為dT,著陸點與激光雷達到跑道垂足的距離為dL,激光光束與跑道的夾角為Θ,下滑道上徑向數據點到y軸的距離為dH,且dH=(y-dL)tan3°。

圖2 飛機下滑道區域的示意圖Fig.2 Schematic diagram of the aircraft glide path

2.4 湍流特征參數反演算法

基于均勻湍流假設條件,在湍流慣性副區當中,可以采用統計學方法進行相關湍流特征參數的計算和反演,其中湍流動能強度和湍流耗散率是最主要關注的量。本方法對機場上空完成一次掃描時間僅需幾分鐘,且PPI掃描水平投影半徑滿足遠大于湍流場外尺度條件,因此認為湍流特征沒有發生明顯變化,符合靜態、均勻條件。

2.4.1 湍流動能強度算法

(3)

其中,〈·〉θ表示按方位角取平均值,使2/tan2φ=1,φ=35.3°,可快速地獲取湍流動能強度計算公式(Turbulence Kinetic Energy,TKE):

(4)

2.4.2 湍流耗散率算法

(5)

其中,tδr=R1-R2要滿足tδr≤Lv;δr為徑向距離間隔;t=0,1,2,…,Lv為湍流外尺度;R1,R2為到激光雷達的徑向距離;k=0,1,2,…,K;N為有效速度波動總次數;θm=m·δθ為方位角,δθ為方位角精度,m=0,1,2,…,M;E(R1,R2)為測量誤差的無偏估計,本文采用脈動速度差的協方差計算測量誤差。

基于Kolmogorov的局地均勻各向同性湍流理論,當徑向距離間隔r?Lv時,其結構函數可以近似表示為[4]:

Dv(r)=Cvε2/3r2/3

(6)

其中,Cv≈2為Kolmogorov常數;ε為湍流耗散率,將計算得到的速度結構函數公式(5)與公式(6)聯立,可得到相應的湍流耗散率強度。

3 比對結果和分析

3.1 仿真設置

仿真激光雷達在每個仰角上作完整的PPI掃描,設置掃描仰角為φ=1°、3°、8°、15°、25°、35.3°、45°、90°,在每個仰角上,雷達掃描的方位角范圍為θ0+m·δθ,其中θ0=0°,δθ=5°,m=0,…,72,探測單元到激光雷達的徑向距離為r0+k·δr,其中r0=0 m,δr=10 m,k=0,…,100,掃描示意圖如圖1所示。為了充分驗證本方法對于在復雜風場探測性能,仿真風場包含均勻項、湍流項及風切變項3部分,具體仿真風場為:

(7)

其中,X,Y,Z∈L×M×N分別為探測單元在以激光雷達為原點的笛卡爾坐標系(如圖1所示)中的x,y,z坐標;Uturb是采用偽小波方法[18],給定湍流耗散率為5×10-4m2/s3仿真所得到的湍流速度,如圖3(a)所示。Ushear是根據風切變的流體力學特征,結合小尺度氣象學通過函數擬合方法[19],給定垂直初始風速為10 m/s強度下所獲得的三維低空風切變速度,如圖3(b)所示。

圖3 湍流和風切變仿真示意圖Fig.3 Simulation of turbulence and wind shear

3.2 仿真結果分析

根據仿真的背景風場,采用本方法分別對反演的三維風場、風切變F因子、湍流動能強度以及湍流耗散率進行計算,并與真值進行了比較。

3.2.1 三維風場反演算法驗證

圖4給出了采用本方法反演的三維風場與仿真風場的對比情況,使用相關系數和均方根誤差用于定量表征風速分量的反演誤差。從圖4中可以看出,采用3D-VPP方法反演的水平風場與仿真風場吻合度較高,相關系數可達0.9以上。因該方法直接使用徑向速度作三角變換計算其他風場分量,因此經反演后合成的徑向速度與初始徑向速度相同,但使用該方法反演垂直風速時誤差較大,修改背景風場中風切變類型為線性風場,其他條件保持不變時,此時各高度角垂直風速之間具有較強相關性,再次使用本方法反演垂直風速時性能較好,均方根誤差可控制在0.7 m/s,性能相對于非線性風場具有較大提升。為驗證本方法的反演性能,采用VAD方法基于線性風場假設對各高度層二維風速進行反演,并與仿真風速進行比對,圖4(c)、(d)所示,VAD方法與3D-VPP方法相比,反演精度較低,相關系數均小于0.7,對于非線性風場反演較差。

圖4 仿真獲取的風場水平分量和3D-VPP、VAD方法反演的結果的比較Fig.4 Comparison of wind field obtained by 3D-VPP and VAD method

2020年9月到11月,利用部署于長沙黃花機場的激光雷達開展了風場探測實驗,激光雷達的系統參數及掃描策略見表1,按本掃描方式完成一次體積掃描需要約3 min。為了對本方法反演的水平風速進行驗證,通過與在同一時間段、同一場地進行探測的探空氣球觀測數據進行比較,追蹤氣球軌跡獲得不同高度上的水平風速以獲得風速的廓線信息。

表1 外場實驗中激光雷達的系統參數及測量參數Tab.1 System parameters and measurement parameters of Lidar

通過與探空氣球測量得到的水平風速進行比較,驗證反演算法反演水平速度的性能。圖5中給出了探空氣球和本方法得到的x,y方向上的水平風速的比較。橫坐標為探空氣球獲得的水平風速的大小,縱坐標表示本方法反演得到的風速的大小,其中反演得到的u與探空氣球數據的比較結果見圖5(a),v的比較結果見圖5(b)。可以發現對于水平速度的反演,采用3D-VPP方法可得到較理想的結果(與探空氣球獲取風速的一致性高),反演水平速度與探空氣球實測速度相關系數均可達0.97以上(u方向為0.9866,v方向為0.9748),為進一步驗證算法反演水平風速的性能,還需積累更多的探空氣球數據。

圖5 本算法反演得到的水平速度與探空氣球數據的比較結果Fig.5 Comparison between the retrieved horizontal wind speed(u,v)and the pilot balloon observed data

3.2.2 風切變反演算法驗證

激光雷達-降落點連線與飛機跑道夾角Θ小于30°,取跑道所在方位角上的3°高度角的徑向速度作為飛機的迎頭風。由于激光雷達只能直接獲取徑向風速數據,只有通過VVP、VPP以及變分法反演才能得到垂直風速,因此在過去關于下滑道風切變F因子的計算中,大多是采用忽略垂直風速,只考慮迎風梯度近似為F因子整體數值。如表2所示,若垂直風速與飛機進場速度(通常為75 m/s)比值大于0.1時,若此時仍忽略垂直項將導致誤差增大,計算結果極易偏離真值,一旦超出國際慣例F因子±0.105風切變閾值,將會出現“虛警”或“漏警”情況,此時不能只考慮迎風梯度計算相應結果。

表2 風切變F因子在有、無垂直風速項的誤差Tab.2 Errors of wind shear F factor with or without vertical wind speed

為提高計算精度,本方法加入通過3D-VPP方法反演得到的垂直風速,使計算結果更加貼近真實數值。如圖6所示,該方法計算所得F因子與仿真結果高度吻合,相關系數超過0.98,均方根誤差為0.013,而無引入垂直風速計算的F因子均方根誤差大約是本方法的2倍,驗證了使用本方法計算F因子具有較高的準確性,對于風切變預警具有較強應用價值。

圖6 仿真獲取的F因子強度和3D-VPP方法反演的F因子的比較Fig.6 Comparison of the F factor obtained by simulation and 3D-VPP method

3.2.3 湍流反演算法驗證圖

圖7(a)、(b)給出了仿真風場和反演風場條件下,在35.3°高度角上采用部分傅里葉分解算法,計算所得的湍流動能強度高度廓線以及兩者相關性情況。從仿真結果來看,兩者具有較強相關性,均方根誤差控制在0.03以內,在不同高度層兩者變化趨勢基本一致,且各高度誤差均較小,因此采用本方法能夠較準確地表征湍流動能強度情況。調整背景湍流場強度,按照強中弱3個等級,采用同樣方法分析,分別得到圖7(a)、(c)、(d),從3幅圖中可以看出,無論在強中弱背景湍流場,采用本方法均能獲得較好的湍流動能強度結果,反演誤差均較小。但在強湍流條件下,如圖7(d)所示,采用本方法獲得的湍流動能強度均方根誤差要遠大于中弱湍流場環境,誤差要大一數量級,因此本方法在中弱湍流強度情況下具有更佳的性能。

圖7 仿真和反演計算所得湍流動能強度對比結果Fig.7 Comparison of turbulence energy intensity obtained by simulation and retrieval

圖8 仿真和反演計算所得湍流耗散率對比結果Fig.8 Comparison of turbulence dissipation rate obtained by simulation and retrieval method

比較本方法在不同湍流強度背景下的反演能力如表3所示,分別設置湍流耗散率為5×10-4m2/s3、1×10-3m2/s3、5×10-3m2/s3、1×10-2m2/s3背景湍流場,采用相同方法進行仿真和反演,可以看出本方法在中低湍流場條件下反演準確度都較高,相對誤差能控制在25 %以內,在強湍流場背景下,使用本方法誤差較大。

表3 仿真和反演湍流耗散率對比結果Tab.3 Comparison of turbulence dissipation rate obtained by simulation and retrieval method

4 結 論

本文提出一種晴空條件下的多用途的測風激光雷達體積掃描策略及風場特征獲取方法,交替采用多個特定高度角的PPI掃描方式,通過3D-VPP方法由徑向風速實現三維風場的反演。基于反演風場,在3°高度角下滑道區域,采用F因子強度評估風切變的存在；利用35.3°高度角圓錐掃描,采用部分傅里葉分解算法計算湍流動能強度；應用均勻湍流假設,計算徑向速度結構函數,從而計算該區域湍流耗散率分布情況。

本文通過實驗仿真,構建包含均勻風場、風切變以及湍流的背景風場,通過對比反演與仿真數據發現:本方法反演的三維風場與仿真風場吻合度較高,水平方向上反演精度更高；引入垂直風速計算風切變F因子,驗證垂直風速與飛機進場速度(通常為75 m/s)比值大于0.1時,垂直風速項將對整體產生較大影響；采用部分傅里葉分解算法計算湍流動能強度時,在中弱湍流強度情況下具有更佳的性能；采用徑向速度方位結構函數計算湍流耗散率時,在中弱湍流場時反演準確度較高,對于強湍流條件下湍流耗散率的反演方法,將在后續工作中進一步進行探索研究。

通過實驗驗證了本方法可有效兼顧空間覆蓋范圍、時空分辨率及魯棒性,既能實現對飛機下滑道區域的重點探測,也能實現對于近場湍流強度的準確感知,還可實現對大范圍三維風場的反演。