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

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

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

0 引言

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

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

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

1 實(shí)驗(yàn)簡介

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

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

Wind3D 6000 三維掃描型相干多普勒測風(fēng)激光雷達(dá)由中國海洋大學(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)激光雷達(dá)與UAT-3 超聲風(fēng)速儀的技術(shù)與性能指標(biāo)。實(shí)驗(yàn)中激光雷達(dá)工作在改進(jìn)型DBS 掃描模式下,使用具有高指向精度的光學(xué)掃描轉(zhuǎn)鏡使激光雷達(dá)光束分別沿東、北、南、西方向以71.38?仰角傾斜出射以及沿垂直方向出射。激光雷達(dá)以30 m 的空間分辨率沿每條視線進(jìn)行徑向風(fēng)速測量,將獲得的五條視線的一組連續(xù)測量定義為一個(gè)掃描序列,即激光雷達(dá)測量的時(shí)間分辨率。設(shè)備旋轉(zhuǎn)掃描一個(gè)序列所用的時(shí)間小于10 s,每一個(gè)掃描序列完成后更新三維風(fēng)矢量反演結(jié)果。

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

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

2.1 尖峰檢測

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

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

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

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

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

圖2 10 月15 日激光雷達(dá)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ù)進(jìn)行湍流統(tǒng)計(jì)分析方法有滑動平均、坐標(biāo)旋轉(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)取平均消除了由于湍流的影響造成的瞬時(shí)風(fēng)速相對于平均值的偏離。

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

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

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

圖3 10 月15 日激光雷達(dá)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 湍流強(qiáng)度結(jié)果分析

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

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

圖4 10 月15–20 日實(shí)驗(yàn)期間60 m 高度處的湍流強(qiáng)度(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)隨時(shí)間的變化是衡量湍流發(fā)展或衰退的指標(biāo)。平均湍流動能ETK對湍流的發(fā)展更具有代表性[18],其通過對湍流動能瞬時(shí)值求平均計(jì)算得到,即

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

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

圖5 10 月15 日激光雷達(dá)60 m 高度處垂直方向湍流強(qiáng)度(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)化為分子熱運(yùn)動動能的速率。激光雷達(dá)計(jì)算湍流動能耗散率的常用方法有回波多普勒譜寬法、結(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的兩點(diǎn)間的速度差。垂直氣流w的二階結(jié)構(gòu)函數(shù)Dww與湍流動能耗散率ε 之間的關(guān)系為

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

圖6 10 月15 日激光雷達(dá)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ū)三個(gè)部分[18],從湍流能譜密度圖中可以看出湍流動能在頻域上的分布特征。

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

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

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