基于相干多普勒激光雷達的斜程湍流參數反演方法研究

陳曉敏 , 張洪瑋 , 孫康聞 , 吳松華 ,2*

( 1 中國海洋大學信息科學與工程學部海洋技術學院, 山東 青島 266100;2 青島海洋科學與技術試點國家實驗室, 區域海洋動力學與數值模擬功能實驗室, 山東 青島 266237 )

0 引 言

大氣湍流作為大氣中一種廣泛存在的運動形式, 具有尺度覆蓋范圍大、變化速度快的特點。大氣湍流對航空安全的影響不可忽視[1], 起飛或降落軌跡上的湍流 (斜程湍流) 會導致飛機顛簸, 嚴重的會造成飛機失控, 引發重大安全事故。飛機起飛或降落軌跡上大氣湍流的準確探測, 不僅能夠降低由大氣湍流導致的飛機事故率, 而且對于提高機場運行效率和研究不同空間位置的大氣湍流特征具有積極意義。

高時空分辨率和高精度的風場數據是低空斜程湍流反演的基礎。傳統的超聲風速計安放高度通常為10 m, 有限的探測高度限制了其在斜程湍流探測中的應用[2];天氣雷達是航空氣象保障中重要的氣象觀測設備, 在濕度較大的天氣下可以發揮作用, 但天氣雷達空間分辨率差、數據刷新率低且在晴空無云的天氣情況下探測能力不足[3]; 風廓線雷達只能夠獲得其所在位置上方的風廓線, 無法獲取較大范圍內的風場情況, 因此無法獲得飛機起降軌跡上的風場數據[4]。近年來, 相干多普勒測風激光雷達憑借時空分辨率高、輕便、掃描模式多樣等特點成為晴空條件下航空氣象保障的可靠設備。相干多普勒測風激光雷達能夠根據不同觀測需求設置不同的掃描模式, 進而獲取低空風場。

針對不同掃描模式采用對應的速度結構函數法反演大氣湍流參數是相干激光雷達觀測大氣湍流的難點和關鍵。近年來, 相關研究已陸續開展。2002年, Frehlich和Cornman[5]利用相干多普勒激光雷達模擬數據的徑向速度估計獲取了模擬湍流速度場的空間統計特性, 并利用徑向速度計算得到的結構函數計算出湍能耗散率ε和積分尺度。2005年, Smalikho等[6]、Frehlich等[7]利用脈沖相干激光雷達距離-高度-顯示 (RHI) 掃描模式的觀測結果從多普勒譜寬與高度結構函數兩個方向反演大氣湍流參數, 并對實驗結果采用數值模擬, 證實了這兩種方法的可靠性。2008年, Frehlich[8]采用固定俯仰角、改變方位角的平面-位置-顯示 (PPI) 掃描模式, 使用縱向結構函數方法和橫向結構函數方法反演邊界層內的湍流參數, 并將反演結果同超聲風速計的探測結果進行了對比, 數據的一致性較好。2011年, Chan和Lee[9]對PPI掃描徑向風場數據進行子扇區劃分, 使用速度結構函數計算得到每個子扇區內的ε1/3, 進而獲取了ε1/3的空間分布情況。2017 年, Smalikho 和Banakh[10]使用速度-方位-顯示 (VAD) 掃描模式下的方位角結構函數反演大氣湍流參數, 并將該方法的適用范圍由對流邊界層擴展到穩定邊界層。2017至2019年, 中國海洋大學翟曉春等[11,12]采用VAD、PPI和RHI等多種激光雷達觀測模式研究了大氣邊界層內的湍流垂直結構特征, 分析了不同粗糙度下墊面影響下的風機尾流與大氣湍流的相互作用特點, 并且探究了大氣湍流對飛機尾渦演化過程的影響。

而針對飛機起降軌跡上的斜程湍流, 香港天文臺開展了一系列飛機下滑道區域內的湍流觀測研究。2007年, Kwong和Chan[13]證明了相干多普勒測風激光雷達能夠在近實時的情況下生成沿飛機滑行路徑的ε剖面。2010 年, Chan[14]利用方位角結構函數和速度結構函數進行了下滑道上的低空湍流反演與預警, 獲取了ε, 兩種計算方法的結果一致, 并與飛行數據以及風切變與湍流預警系統 (WTWS) 的結果一致。2014年,Hon和Chan[15]將速度結構函數法應用到針對飛機起降階段沿滑行路徑區域掃描的下滑道掃描模式, 獲得了沿滑行路徑的ε1/3剖面, 與飛機的飛行員報告匹配較好。2010年, Chan[14]利用方位角結構函數和速度結構函數進行了下滑道上的低空湍流反演與預警, 獲取了ε, 兩種計算方法的結果一致, 并與飛行數據和風切變與湍流預警系統 (WTWS) 結果一致。然而對于斜程空間范圍內的湍流參數反演研究依舊處于起步階段, 同時缺乏對于斜程空間湍流的連續觀測。

本文利用安放于某機場的相干多普勒測風激光雷達, 獲取下滑道掃描模式下的高時空分辨率風場數據,將獲取的徑向風速數據進行5°方位角和60 m距離的扇區劃分, 進一步基于Kolmogorov局地均勻各向同性湍流理論, 并使用橫向結構函數法估算斜程湍流的ε1/3分布情況。具體分析了某機場南段跑道2018年12月15日凌晨 01:38―01:50 以及 02:22―02:35 兩個時間段的觀測數據, 獲得的ε1/3斜程空間分布情況。最后, 在同區域、同時段內, 使用相同數據計算了與ε1/3相同量綱的風切變強度, 并將二者進行了對比驗證。

1 理論方法

大氣湍流在產生、發展和消散的過程中會形成各種尺度的湍渦, 某些小尺度湍渦具有各向同性的特性,可以進行統計模擬, 因此在實際研究中, 小尺度湍渦的研究是目前湍流研究領域的關注點。在Kolmogorov局地均勻各向同性湍流理論的假設中, 當雷諾數足夠高時, 存在湍流強度僅由ε確定的慣性子區。國際民用航空組織 (ICAO) 規定, 使用ε1/3來衡量大氣湍流強度: 0.3～0.5 m2/3/s為中度湍流, 大于0.5 m2/3/s為嚴重湍流。

速度結構函數法可利用相干多普勒激光雷達的徑向風速數據估算滿足Kolmogorov局地均勻各向同性理論假設下的ε值, 是當前相干多普勒激光雷達獲取大氣湍流ε的常用方法。該方法將由實測數據計算的速度結構函數與模型結構函數進行擬合, 獲得ε值。在擬合過程中, 需要根據Kolmogorov局地均勻各向同性理論對擬合出的湍流特征尺度進行限制, 以確保其處于慣性子區內。

1.1 Kolmogorov局地均勻各向同性理論

根據Kolmogorov定律, 當雷諾數足夠大時, 湍流獲得充分發展。在湍渦的發展過程中, 大尺度湍渦的能量傳遞給次級的湍渦, 顯然是各向異性的。但在串級傳輸的過程中, 存在一部分在統計特性上是各向同性的小尺度湍渦, 即在局地均勻各向同性區域內存在一個僅由ε確定的慣性子區, 其尺度l滿足L≥l≥η, 其中L為含能尺度, 一般來說研究湍流時不考慮比含能尺度L更大的尺度;η為湍流內尺度, 最小尺度約為1 mm。在湍流的慣性子區內, 速度結構函數頻譜只與ε有關[16,17]。

1.2 速度結構函數法

速度結構函數法需要對實測數據計算的速度結構函數同模型結構函數進行擬合來獲取湍流參數, 在實際湍流研究中, 主要關注的是慣性子區內湍渦的貢獻。速度結構函數根據數據的擴展方向分為橫向結構函數和縱向結構函數兩種。在本研究中激光雷達下滑道掃描模式的方位角分辨率高, 因此使用橫向結構函數計算方法可以更加精確地反演大氣湍流參數。橫向結構函數Dv(s) 表示為

式中v′(r，θ，φ) 為實測風場減掉平均風場獲得的風速脈動,r為觀測點距激光雷達的距離,θ為俯仰角,φ為方位角,s為兩項風速脈動之間的距離間隔。根據式 (1) 可獲得基于實測數據計算的結構函數。

近地面層區域的模型結構函數可以簡化表示為

式中σ2為徑向風速方差, Λ(x) 為通用函數,L為含能尺度。

當L＜＜R(R表示激光光束方向距激光雷達的距離), Von Karman 模型的二維空間湍流模型為

式中q=(r2+s2)1/2,K1/3(x) 為修正的1/3階貝塞爾函數。

當s＜＜L且滿足各向均勻同性假設時, 滿足Kolmogorov定理, 結構函數可以表示為

由于激光雷達光束始終存在一定的脈沖寬度, 實際測量中探測體積平均效應無法避免, 這將影響結構函數計算。Frehlich 等[7]基于上述理論推導出探測體積平均效應下的結構函數模型。

當激光雷達的橫向分辨率Rsin Δφ＜＜Lv時, 則模型結構函數為

式中F(x，μ) 為高斯分布脈沖和方波窗函數下的濾波函數, 其計算公式為

式中Dmeasure(kΔs)為實測數據計算的結構函數,kΔs為結構函數距離分辨率的倍數。根據式 (10) 進行最小二乘法擬合, 即可求得湍流參數徑向風速方差σ和含能尺度L, 最終通過式 (6) 求得ε。

1.3 斜程湍流反演算法

圖1為斜程湍流參數反演算法流程圖。使用速度結構函數法反演局地均勻各向同性假設下的斜程湍流參數, 首先需要對激光雷達獲取的風速數據進行質量控制, 通過設置信噪比閾值的方法排除奇異數據。下滑道掃描模式的俯仰角與方位角同時變換, 因此獲取的是空間傾斜平面上的風速數據。假設一定高度范圍內的風場是均勻的, 對實測風場數據進行高度平均和時間平均獲取平均風場, 進而獲得風速脈動, 本研究采用15 min時間平均和2 m高度平均。圖2展示了2018年12月15日01:41:17的平均風場分布和風速脈動分布, 其中x軸為距激光雷達北-南方向的距離,y軸為距激光雷達東-西方向的距離,z軸為距激光雷達的垂直高度。

圖1 算法流程圖Fig.1 Flow chart of algorithm

圖2 2018年12月15日01:41:17的平均風場分布圖 (a) 和風速脈動分布圖 (b)Fig.2 Mean radial velocity (a) and radial velocity fluctuation (b) of wind field at 01:41:17 on 15 December 2018

為獲取傾斜下滑道區域內的不同空間位置的湍流參數并提高運算效率, 采用空間區域劃分方法將掃描區域劃分為多個子扇區 (每個子扇區尺度為 5°方位角、2個距離庫), 計算每個子扇區內的橫向結構函數并使用自協方差法進行系統隨機誤差校正。后將計算的結構函數與模型結構函數進行最小二乘法擬合, 獲得大氣湍流參數ε。圖3為距離激光雷達210 m處的計算結構函數與修正后的模型結構函數的擬合效果圖, 此時對應的橫向分辨率為 0.2 m, 遠小于擬合出的湍流積分尺度74.28 m, 符合橫向結構函數法的適用前提, 反演出的湍流參數數值可信。

圖3 2018年12月15日01:44:28距激光雷達210 m處的結構函數估計Fig.3 Structure function estimation of turbulence at 210 m from the lidar at 01:44:28 on 15 December 2018

2 結果與分析

2.1 實驗設置

中國海洋大學于2018年11月至2019年8月在某機場開展湍流觀測實驗。機場所在區域于秋冬季節盛行西北風, 且風速較大, 跑道端周圍多建筑與高大樹木, 下墊面情況復雜, 易產生地形誘導風切變。實驗期間, 將青島鐳測創芯科技有限公司研制的Wind3D 6000-AP相干測風激光雷達 (詳細參數見表1) 安放于機場跑道一端, 安放位置示意圖如圖4所示。

圖4 機場下滑道湍流觀測模式示意圖Fig.4 Schematic diagram for glide path scanning mode of lidar

表1 Wind3D 6000-AP相干測風激光雷達技術指標Table 1 Technical specifications of Wind3D 6000-AP coherent wind lidar

實驗中激光雷達的觀測模式為下滑道掃描模式。本次實驗中激光雷達安放位置距離跑道邊緣線約為70 m, 以保證激光光束與下滑道夾角不大于30°, 減小側風分量的影響, 使激光雷達徑向風速可以近似代替飛機真實經歷的風速(真實風速在跑道方向上的風速分量)。下滑道掃描模式下的激光雷達掃描參數為: 俯仰角范圍3°～4°, 方位角范圍148°～174°, 方位角的掃描速度約為0.2 °/s～0.3 °/s, 距離分辨率為30 m, 風速測量精度約0.1 m/s, 風向測量精度約5°。

2.2 數據分析和討論

2.2.1 穩定大氣條件下的觀測數據

安放在機場跑道南端的激光雷達于2018年12月15日01:38―01:50獲取的徑向風速分布如圖5所示, 其中x軸為距激光雷達北-南方向的距離,y軸為距激光雷達東-西方向的距離,z軸為距激光雷達的垂直高度。圖中顏色表示風速大小, 暖色表示風向遠離激光雷達, 冷色表示風向朝向激光雷達。通過圖5的序列圖可以發現該時段內近地面的風速普遍高于1.5 m/s, 且不存在明顯的大氣擾動現象, 揭示該段時間內大氣狀態穩定。

圖5 2018年12月15日01:38―01:50徑向風速分布圖Fig. 5 Distribution of radial wind velocity during 01:38―01:50 on 15 December 2018

將每次掃描反演得到的不同空間位置的大氣湍流參數按照空間位置作圖, 如圖6所示。圖中的原點位置為激光雷達的布放位置,x軸為距激光雷達北-南方向的距離,y軸為距激光雷達東-西方向的距離,z軸為距激光雷達的垂直高度。整個下滑道掃描區域被劃分為45個子扇區, 每個子扇區的方位角掃描范圍為5°、徑向距離范圍為60 m。每個子扇區反演獲得子扇區內空間的平均湍流特征參數, 將每次掃描的148.4°～174°方位角范圍、60～600 m徑向距離范圍區域中的45個子扇區內的ε1/3取均值, 代表一次掃描區域內的平均大氣湍流強度狀況,ε1/3的平均值和相應的標準差 (——-ε13±σ) 在每個圖窗中進行了標注。由圖6可知, 2018年12月15 日01:38―01:50 內, 湍流狀況較為穩定。ε1/3值主要分布在0.03～0.07 m2/3/s 之間, 每次掃描區域內的空間平均ε1/3為0.047～0.059 m2/3/s。在01:42:53, 激光雷達觀測到了發生在其南側約60～250 m、東側0～100 m范圍內一次相對較強的湍流現象,ε1/3最大值達到0.16 m2/3/s。同樣地, 在01:47:39, 激光雷達南側約60～300 m、東側0～100 m范圍處也出現了一次相對較強的湍流現象, 激光雷達估計的ε1/3最大值達到了0.12 m2/3/s, 屬于輕度湍流。

2.2.2 波動大氣狀態案例分析

2018年12月15日凌晨02:22―02:35, 跑道南端的激光雷達共進行了9次掃描, 反演獲取了約13min的徑向風速, 如圖7所示。該時間段內的ε1/3空間分布圖如圖8所示, 每次掃描空間內ε1/3的平均值和相應的標準差 (——-ε13±σ)在每個圖窗中進行了標注。圖7、圖8的坐標軸與圖5、圖6一致。結果顯示該時段的近地面徑向風速小于1 m/s, 且存在較為明顯的大氣擾動現象?？梢钥闯? 2018年12月15日02:22―02:35,ε1/3主要分布在0.04～0.08 m2/3/s之間, 每次掃描區域內的ε1/3空間平均值為0.044～0.081 m2/3/s。在02:30:51的觀測數據顯示下滑道區域內存在非常明顯的湍流強度增大現象, 位置大概在激光雷達南側150～450 m、東側0～150 m范圍內,ε1/3最大值達到0.28 m2/3/s。

圖6 2018年12月15日01:38―01:50 ε1/3 的空間位置分布圖Fig.6 Spatiotemporal distributions of ε1/3 during 01:38―01:50 on 15 December 2018

圖7 2018年12月15日02:22―02:35徑向風速分布圖Fig.7 Distribution of the radial wind velocity during 02:22―02:35 on 15 December 2018

2.2.3 湍流特征參數與風切變強度對比印證

為驗證湍流特征參數反演算法的準確性, 采用經驗證的下滑道風切變識別法[18]計算的風切變強度值進行對比印證。下滑道風切變識別法通過逆風廓線的增大或減小來判斷飛機在某一特定方向上是否遭遇風切變。該方法可以進行風切變強度值Iraw的計算, 即

式中 ΔV是一段距離內風速的總變化值;Vapp是飛機的平均進近速度, 一般取 75 m/s; ΔV/R1/3是風速在下滑道上的變化率,R是風切變發生的長度, 即飛行軌跡上風切變的長度。與ε相似, 同樣以風切變強度值Iraw的立方根進行強度等級的劃分, 當該值在0.3～0.5 m1/3/s2/3范圍時, 為中等強度風切變; 當該值大于0.5 m1/3/s2/3時,為嚴重風切變[19,20]。

對Iraw與ε進行量綱分析, 發現兩者量綱相差m/s, 同時, 飛機進近速度Vapp為固定常數, 與觀測數據無關。因此, 為更好地對ε1/3對比印證, 將式(11)中的平均進近速度項Vapp去除, 則調整后的風切變強度值I與ε的量綱相同, 即

調整后的風切變強度值I的立方根在進行強度等級劃分時, 也需要相應調整為I的立方根。在1.3～2.1 m2/3/s范圍時, 為中等風切變強度;大于2.1 m2/3/s時為嚴重風切變。

圖9為2.2.1節和2.2.2節所述兩個觀測個例中01:42:53與02:30:51的ε1/3空間分布情況與對應時刻的逆風廓線對比。圖9 (a)和圖9 (c)中, 原點為激光雷達的布放位置,x軸為距激光雷達的距離,y軸為逆風風速值?？梢钥闯鰞纱未髿鈹_動在逆風廓線數據和數據中都有所體現, 并且兩種方法獲得的結果在時間上較為一致。兩次擾動較弱, 均未達到輕度風切變閾值, 在逆風廓線中并無明顯體現, 但ε1/3的最大值達到輕度湍流閾值,且在ε1/3的空間分布圖中能夠被明顯觀察到。下滑道風切變探測反演方法通過衡量較短時間內風矢量的突變情況判斷風切變的發生位置與強度。相較于風切變強度值, 應用下滑道掃描模式的橫向速度結構函數法反演得到的湍流參數去除了背景平均風場, 對空間內存在的小尺度湍渦更加敏感。

圖9 2018年12月15日01:42:53逆風廓線 (a) 與ε1/3空間分布 (b) 以及02:30:51逆風廓線 (c) 與ε1/3 空間分布 (d)Fig.9 Headwind profile (a) and spatiotemporal distributions of ε1/3 (b) at 01:42:53, headwind profile (c) and spatiotemporal distributions of ε1/3 (d) at 02:30:51 on 15 December 2018

圖10 (a)和圖10 (b)分別為2018年12月15日 01:38―01:50時段內9個掃描片段與02:18―02:38時段內14個掃描片段的I1/3與單次下滑道掃描區域內的空間平均ε1/3對比。在01:38―01:50時段, 激光雷達觀測范圍內的大氣狀態穩定,I1/3與空間平均ε1/3均無明顯起伏; 在02:18―02:38時段激光雷達觀測到一次較為明顯的大氣波動現象 [02:30―02:34 時段, 如圖8 (f)―(h) 所示], 兩個強度因子隨時間的變化趨勢具有較好的一致性,I1/3與空間平均ε1/3數值均達到該時段內的最大值。此次對比基本可以印證本文下滑道掃描模式湍流參數反演結果的可信度。

圖8 2018年12月15日02:22―02:35 ε1/3 的空間位置分布圖Fig.8 Spatiotemporal distributions of ε1/3 during 02:22―02:35 on 15 December 2018

圖10 2018年12月15日 01:38―01:50 (a) 和02:18―02:38 (b) 時段的ε1/3 與風切變強度對比Fig.10 Comparison of ε1/3 and I1/3 during 01:38―01:50 (a) and 02:18―02:38 (b) on 15 December 2018

3 結 論

應用高時空分辨率相干多普勒激光雷達觀測數據, 采用橫向結構函數法估算出凌晨兩個時間段 (01:38―01:50、02:22―02:35)下滑道掃描模式的斜程湍流參數, 分析了觀測時段內大氣湍流參數的斜程空間分布情況, 捕捉到了較為穩定大氣狀態下的輕微大氣擾動。將估算的大氣湍流強度與使用相同數據、不同計算方法計算的風切變強度值進行比對, 二者的變化趨勢具有較好的一致性。然而風切變強度只代表下滑道方向上的風切變強度水平, 而空間平均后的ε1/3代表一次下滑道掃描范圍內的整體湍流強度水平, 因此兩者在時間上的變化趨勢無法完全對應。

斜程湍流參數的觀測反演研究具有廣闊的應用前景, 如在氣象觀測中獲取復雜地形條件下的湍流信息,校正天文望遠鏡在斜程觀測中因大氣湍流造成的觀測誤差等。在后期的實驗中計劃在激光雷達的掃描區域內設置多個安裝超聲風速計的測風桿, 利用超聲風速計的數據對激光雷達反演的湍流參數進行驗證對比, 修正大氣湍流反演算法, 進一步提高激光雷達的大氣湍流探測反演精度。

致謝: 感謝青島鐳測創芯科技有限公司在現場觀測實驗和數據獲取中提供的幫助， 感謝中國海洋大學激光雷達課題組的戴光耀、劉曉英在問題討論過程中提出的建議。