劉佳鑫,云龍,邵士勇,程雪玲,宋小全,5?

(1 中國海洋大學(xué)信息科學(xué)與工程學(xué)部海洋技術(shù)學(xué)院,山東 青島 266100;

2 深圳市環(huán)境監(jiān)測中心站,廣東 深圳 518049;3 中國科學(xué)院合肥物質(zhì)科學(xué)研究院安徽光學(xué)精密機械研究所,中國科學(xué)院大氣光學(xué)重點實驗室,安徽 合肥 230031;4 中國科學(xué)院大氣物理研究所大氣邊界層物理和大氣化學(xué)國家重點實驗室,北京 100029;5 青島海洋科學(xué)與技術(shù)試點國家實驗室區(qū)域海洋動力學(xué)與數(shù)值模擬功能實驗室,山東 青島 266237)

0 引言

受到熱力和動力因素作用,當雷諾系數(shù)超過臨界值時,大氣的分層不再穩(wěn)定,表現(xiàn)出如漩渦式紊亂、不規(guī)則的運動形式,這被稱為湍流[1]。湍流影響著包括空氣質(zhì)量、風(fēng)能、航空和數(shù)值天氣預(yù)報等在內(nèi)的許多領(lǐng)域。由湍流引起的下墊面與大氣之間動量、熱量和物質(zhì)交換的垂直輸送直接影響著人們生活和生產(chǎn)活動,在大氣運動和天氣演變中起著至關(guān)重要的作用[2]。因此研究湍流相關(guān)特性及其時空分布特征具有重要的意義。

安裝于氣象梯度塔上的風(fēng)杯式風(fēng)速計是觀測風(fēng)速湍流的傳統(tǒng)儀器,但其受到梯度塔塔高的限制,并存在風(fēng)速傳感器啟動速度等問題。超聲風(fēng)速儀比風(fēng)杯式風(fēng)速計更精確地測量湍流,能夠以較高的采樣頻率記錄三維風(fēng)速信息,但只能做到定點觀測,并受周圍環(huán)境影響,導(dǎo)致風(fēng)速和湍流的測量不準確。針對這些問題,聲雷達和激光雷達等遙感設(shè)備被應(yīng)用于風(fēng)湍流的研究。近年來,大氣湍流參數(shù)及其垂直廓線的測量方法及設(shè)備不斷發(fā)展完善,有微波雷達、聲雷達、激光雷達測量和探空式測量等[3?6]。微波雷達及聲雷達測量受到空間分辨率的限制,且可移動性較差;探空式測量易受風(fēng)力影響,定點性和實時性較差。激光雷達具有高時空分辨率的特點[7],使其在反演大氣湍流參數(shù)等方面具有獨特優(yōu)勢,它能夠獲得遠高于傳統(tǒng)氣象梯度塔高度處的廓線信息,且對所測量的大氣環(huán)境無影響。湍流強度廓線激光雷達在大氣湍流的研究中有較多應(yīng)用[8]。此外,作為目前風(fēng)場測量最先進的設(shè)備,相干多普勒激光雷達的測風(fēng)技術(shù)迅速成熟[9],進而較多地被應(yīng)用于湍流的監(jiān)測及研究[10],通過檢測回波信號的多普勒頻移可實現(xiàn)包括大氣邊界層在內(nèi)的對流層三維風(fēng)場的精細化探測。

本文使用多普勒測風(fēng)激光雷達數(shù)據(jù)并結(jié)合地面氣象站的常規(guī)氣象資料,對深圳楊梅坑大氣湍流運動的基本日變化特征展開探索性研究,其結(jié)果為激光雷達探測和半島地區(qū)湍流運動的研究提供重要參考。

1 實驗簡介

2019 年10 月14–21 日,利用相干多普勒測風(fēng)激光雷達在深圳大鵬半島上的楊梅坑監(jiān)測站樓頂平臺(22?33.3′N,114?35.2′E,海拔約60 m)進行大氣風(fēng)場和湍流觀測實驗。如圖1 所示,半島位于珠江口東側(cè)大鵬灣與大亞灣之間,監(jiān)測站位于半島丘陵地帶的山腳處,面朝大亞灣[11]。監(jiān)測站測風(fēng)塔(海拔約120 m)上安裝的測風(fēng)儀器為中國科學(xué)院大氣物理研究所研制的UAT-3 型超聲風(fēng)速儀,與多普勒測風(fēng)激光雷達的水平距離約100 m,垂直方向高度差約60 m。

圖1 楊梅坑地理位置(a)、實驗場地(b)與儀器(c)Fig.1 Geographical location(a),experimental site(b)and equipment(c)

Wind3D 6000 三維掃描型相干多普勒測風(fēng)激光雷達由中國海洋大學(xué)與青島鐳測創(chuàng)芯科技有限公司聯(lián)合研制,它基于光學(xué)脈沖相干多普勒頻移檢測的原理,具有DBS(Doppler beam swinging)和VAD(Velocity azimuth display)等多種掃描模式,已應(yīng)用于低空急流、海陸風(fēng)等相關(guān)研究[12?14]。表1 為Wind3D 6000 型多普勒測風(fēng)激光雷達與UAT-3 超聲風(fēng)速儀的技術(shù)與性能指標。實驗中激光雷達工作在改進型DBS 掃描模式下,使用具有高指向精度的光學(xué)掃描轉(zhuǎn)鏡使激光雷達光束分別沿東、北、南、西方向以71.38?仰角傾斜出射以及沿垂直方向出射。激光雷達以30 m 的空間分辨率沿每條視線進行徑向風(fēng)速測量,將獲得的五條視線的一組連續(xù)測量定義為一個掃描序列,即激光雷達測量的時間分辨率。設(shè)備旋轉(zhuǎn)掃描一個序列所用的時間小于10 s,每一個掃描序列完成后更新三維風(fēng)矢量反演結(jié)果。

表1 Wind3D 6000 多普勒測風(fēng)激光雷達和UAT-3 超聲風(fēng)速儀技術(shù)指標Table 1 Technical specifications of Wind3D 6000 Doppler wind lidar and UAT-3 ultrasonic anemometer

2 風(fēng)速湍流的獲取與分析

2.1 尖峰檢測

使用多普勒測風(fēng)激光雷達獲取的距離激光雷達出光口60 m 高度處的風(fēng)速數(shù)據(jù)研究湍流的日變化特征,該高度與超聲風(fēng)速儀安裝位置高度大體相當。由于測風(fēng)激光雷達獲取的是徑向風(fēng)速信息,所以首先需使用徑向風(fēng)速反演三維風(fēng)矢量,其計算公式為

式中:遵循標準氣象慣例,u表示自然坐標系下風(fēng)矢量在東西方向分量,u>0 表示西風(fēng);v表示風(fēng)矢量在南北方向分量,v>0 表示南風(fēng);w是垂直氣流,w>0 表示向上運動;VE、VW、VS、VN和VZ分別表示東西南北和垂直方向上激光雷達測量的徑向速度,V>0 表示遠離激光雷達的正向運動,V<0 表示朝向激光雷達的負向運動。

因為湍流參數(shù)的值依賴于時間序列中風(fēng)速測量的不確定性,實驗過程中的干擾信號產(chǎn)生的異常尖峰會導(dǎo)致信號局部出現(xiàn)異常高的湍流強度。為了減輕這種影響,對三維風(fēng)矢量進行尖峰檢驗,以消除異常的高風(fēng)速值[15]。i時刻的風(fēng)速預(yù)測值ufcst,i的表達式為

式中:cov(τ)表示滯后時間τ 的信號自協(xié)方差;uL,i?1是i?1 時刻風(fēng)速的觀測值;表示平均風(fēng)速。上述參數(shù)的計算使用每個測量時刻前10 min 的數(shù)據(jù)點。當獲得新的風(fēng)速值時,在假設(shè)ufcst,i?uL,i符合高斯分布的前提下,使用如下準則檢測可能存在的異常值,即

式中:Cspike是尖峰檢測的閾值,第一次檢測取3,每次迭代后增加0.1;σ 是標準偏差。尖峰檢測的窗口在整個數(shù)據(jù)集上移動,檢測到的尖峰將使用相鄰的非尖峰值的線性插值代替,并重復(fù)尖峰檢測的過程,直到不再發(fā)現(xiàn)尖峰。尖峰檢測的結(jié)果如圖2 所示,檢測出的異常值使用三角符號標記,圓形符號表示經(jīng)三次保形樣條插值處理獲得的校正值,從中可以看到尖峰檢測已較好地剔除了原始數(shù)據(jù)中的異常值。實驗中使用的數(shù)據(jù)有效率達到99.1%。

圖2 10 月15 日激光雷達60 m 高度處尖峰檢測異常值剔除Fig.2 Elimination of outlier values during peak detection at 60 m height on October 15

2.2 風(fēng)速湍流獲取

對三維風(fēng)速數(shù)據(jù)進行湍流統(tǒng)計分析方法有滑動平均、坐標旋轉(zhuǎn)和小波分解等[16]。在近地面風(fēng)速脈動能譜中,湍流對應(yīng)的微尺度渦具有10 min 的周期?；瑒悠骄ㄒ罁?jù)Reynolds 分解的原理,將風(fēng)速分解為平均部分和湍流部分的疊加形式,即

式中:w′表示垂直氣流的脈動,表示平均垂直氣流,滑動平均周期T取10 min。風(fēng)速數(shù)據(jù)在10 min 內(nèi)取平均消除了由于湍流的影響造成的瞬時風(fēng)速相對于平均值的偏離。

坐標旋轉(zhuǎn)法通過兩次坐標旋轉(zhuǎn)獲取風(fēng)速湍流[17]。第一次坐標旋轉(zhuǎn)以z軸為中心軸,將坐標系x軸旋轉(zhuǎn)至水平風(fēng)的合成風(fēng)方向;第二次旋轉(zhuǎn)是以旋轉(zhuǎn)后的y軸為中心軸,使旋轉(zhuǎn)后的坐標系x軸與三維風(fēng)矢量的合成風(fēng)方向平行。地形復(fù)雜情形時一般不進行兩次旋轉(zhuǎn),以避免由地形引起的垂直氣流被消掉。

本研究將小波分解的原理應(yīng)用于獲取風(fēng)速湍流并進行湍流參數(shù)的分析。原始信號經(jīng)過小波分解得到低頻近似部分和高頻細節(jié)部分,近似部分反映信號的整體趨勢,細節(jié)部分反映信號在整體趨勢項上的細節(jié)波動。

式中:aj+1和bj+1分別為低頻分量和高頻分量,j為分解的層數(shù);H、G分別為低頻和高頻分解函數(shù),類似濾波器。小波分解得到的脈動風(fēng)速數(shù)據(jù)在統(tǒng)計上應(yīng)滿足在10 min 的時間尺度上積分為0,因此使用小波分解法獲得的第六層脈動風(fēng)速來獲得湍流參數(shù)?；瑒悠骄ㄅc小波分解得到的平均風(fēng)速如圖3 所示,平均垂直氣流出現(xiàn)負值表示平均垂直氣流向下運動?；瑒悠骄ㄝ^為簡單,但存在邊緣效應(yīng),即在邊緣部分平均數(shù)據(jù)的變化趨勢不準確,而小波分解能夠較好地處理邊緣數(shù)據(jù)。

圖3 10 月15 日激光雷達60 m 高度處滑動平均與小波分解獲得的平均風(fēng)速Fig.3 Average wind speed obtained by moving average and wavelet decomposition at 60 m height on October 15

3 湍流參數(shù)反演結(jié)果分析

3.1 湍流強度結(jié)果分析

湍流強度(Turbulence intensity,TI)是表征風(fēng)速脈動強度的物理量。垂直方向上湍流強度Iw的物理定義為10 min 時距內(nèi)的脈動風(fēng)速標準差σw和平均風(fēng)速U的比值,其表達式為

實驗期間深圳楊梅坑地區(qū)日出時刻為06:20 左右,日落時刻為17:57 左右。日平均湍流強度與平均溫度的變化如圖4 所示,湍流強度的變化具有典型的帽式特征,形態(tài)上表現(xiàn)為“單峰單谷”。日出前湍流強度幾乎沒有變化,表現(xiàn)為谷值;日出后,湍流強度基本隨著溫度的升高而增高,午后02:00 湍流強度達到最大值;日落后,湍流強度隨著溫度的降低而逐漸降低。

圖4 10 月15–20 日實驗期間60 m 高度處的湍流強度(a)和溫度(b)日平均變化Fig.4 Daily average change of turbulence intensity(a)and temperature(b)at 60 m height from October 15 to 20

3.2 湍流動能結(jié)果分析

湍流動能(Turbulent kinetic energy,TKE)隨時間的變化是衡量湍流發(fā)展或衰退的指標。平均湍流動能ETK對湍流的發(fā)展更具有代表性[18],其通過對湍流動能瞬時值求平均計算得到,即

式中:u′、v′、w′分別為東西方向風(fēng)分量u、南北方向的風(fēng)分量v、垂直氣流w的脈動量。在理想條件下,湍流動能的計算應(yīng)通過同一時間、同一觀測點的多個測量值的系綜平均來獲得。在實際情況中無法達到系綜平均的要求,通常在滿足局部均勻假設(shè)的條件下,使用一段時間的平均值來代替[19],這里借助具有多分辨率分析的小波變換來分析垂直方向上的湍流動能。小波變換通過伸縮和平移對信號逐次進行多尺度細化,特別適用于如風(fēng)速等非平穩(wěn)信號的分析。

圖5 為10 月15 日垂直方向的湍流強度與湍流動能的小波系數(shù),從中可以看出,二者峰值所對應(yīng)的時間段幾乎相一致,由湍流引起的垂直輸運對應(yīng)的時間尺度主要集中在午后至傍晚這一段時間,夜晚的湍流輸送較弱,這與一天中溫度的變化有較為一致的趨勢。

圖5 10 月15 日激光雷達60 m 高度處垂直方向湍流強度(a)和垂直方向湍流動能的小波系數(shù)(b)Fig.5 Vertical turbulence intensity(a)and wavelet coefficients of vertical turbulent kinetic energy(b)at 60 m height above lidar on October 15

3.3 湍流動能耗散率結(jié)果分析

湍流動能耗散率(Turbulent kinetic energy dissipation rate)是表征湍流動能在分子的粘性作用下轉(zhuǎn)化為分子熱運動動能的速率。激光雷達計算湍流動能耗散率的常用方法有回波多普勒譜寬法、結(jié)構(gòu)函數(shù)估算法和相似理論估算法等[20,21]。由于記錄回波信號多普勒譜寬的數(shù)據(jù)量較大,給數(shù)據(jù)的存儲帶來了較大的困難。這里采用二階結(jié)構(gòu)函數(shù)估算法估算湍流動能耗散率。垂直氣流w的二階結(jié)構(gòu)函數(shù)Dww定義為

式中:w(x+r)?w(x)表示距離為r的兩點間的速度差。垂直氣流w的二階結(jié)構(gòu)函數(shù)Dww與湍流動能耗散率ε 之間的關(guān)系為

圖6 為10 月15 日激光雷達在60 m 高度處湍流動能與湍流動能耗散率。由圖可知,湍流動能耗散率的數(shù)值位于10?2～10?4之間,這符合Banakh 等[10]的研究結(jié)果。垂直方向湍流動能與湍流動能耗散率有較為一致的變化趨勢,根據(jù)公式[17],二者的相關(guān)系數(shù)為0.77,呈顯著的相關(guān)性。

圖6 10 月15 日激光雷達60 m 高度處湍流動能與湍流動能耗散率Fig.6 Turbulent kinetic energy and energy dissipation rate at 60 m height on October 15

3.4 湍流速度功率譜密度結(jié)果分析

湍流能譜密度(Turbulent velocity power spectral density)一般分為含能渦區(qū)、慣性副區(qū)和耗散區(qū)三個部分[18],從湍流能譜密度圖中可以看出湍流動能在頻域上的分布特征。

依據(jù)Kolmogorov 提出的相似理論,在慣性副區(qū)中湍流特性僅由湍流耗散率決定,湍流的能譜密度表示為

式中:aw是Kolmogorov 常數(shù),取值0.5;ε 為湍流動能耗散率,本研究使用一天中能耗率的平均值0.0022。對10 min 數(shù)據(jù)的時間相關(guān)函數(shù)做傅氏變換得到湍流能譜密度,找到慣性副區(qū)所對應(yīng)的頻率范圍,比較慣性副區(qū)內(nèi)的擬合斜率與Kolmogorov 能譜密度的斜率,如圖7 所示,二者的關(guān)系基本符合“-5/3”次律[22]。

圖7 10 月15 日激光雷達60 m 高度處的譜密度在慣性副區(qū)的斜率Fig.7 The slope of the spectral density in the inertial subregion at 60 m height on October 15

4 結(jié)論

依據(jù)測風(fēng)激光雷達獲取的徑向風(fēng)速信息可以獲得描述湍流動力特征的湍流強度和湍流動能耗散率等相關(guān)參數(shù),綜合分析了湍流的時空變化特征。小波分解能較好地從風(fēng)速數(shù)據(jù)中識別湍流脈動。湍流強度的日平均變化結(jié)果基本符合帽式結(jié)構(gòu),在時間上的變化與溫度的變化呈現(xiàn)較高相關(guān)性。使用小波變換獲得湍流動能在時域及頻域上的分布,湍流強度與湍流動能在時間尺度上的變化相對應(yīng)。使用二階結(jié)構(gòu)函數(shù)計算得到的湍流耗散率與湍流動能顯著相關(guān)。湍流能譜密度與頻率的關(guān)系在慣性副區(qū)范圍內(nèi)較為符合Kolmogorov 的“-5/3”冪律。實驗結(jié)果驗證了Wind3D 6000 型相干多普勒測風(fēng)激光雷達獲取三維風(fēng)速觀測湍流的可行性。