深圳地區(qū)多普勒測(cè)風(fēng)激光雷達(dá)的湍流觀測(cè)

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

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

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

0 引言

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

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

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

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

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

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

Wind3D 6000 三維掃描型相干多普勒測(cè)風(fēng)激光雷達(dá)由中國(guó)海洋大學(xué)與青島鐳測(cè)創(chuàng)芯科技有限公司聯(lián)合研制,它基于光學(xué)脈沖相干多普勒頻移檢測(cè)的原理,具有DBS(Doppler beam swinging)和VAD(Velocity azimuth display)等多種掃描模式,已應(yīng)用于低空急流、海陸風(fēng)等相關(guān)研究[12?14]。表1 為Wind3D 6000 型多普勒測(cè)風(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)速測(cè)量,將獲得的五條視線的一組連續(xù)測(cè)量定義為一個(gè)掃描序列,即激光雷達(dá)測(cè)量的時(shí)間分辨率。設(shè)備旋轉(zhuǎn)掃描一個(gè)序列所用的時(shí)間小于10 s,每一個(gè)掃描序列完成后更新三維風(fēng)矢量反演結(jié)果。

表1 Wind3D 6000 多普勒測(cè)風(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 尖峰檢測(cè)

使用多普勒測(cè)風(fēng)激光雷達(dá)獲取的距離激光雷達(dá)出光口60 m 高度處的風(fēng)速數(shù)據(jù)研究湍流的日變化特征,該高度與超聲風(fēng)速儀安裝位置高度大體相當(dāng)。由于測(cè)風(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)動(dòng);VE、VW、VS、VN和VZ分別表示東西南北和垂直方向上激光雷達(dá)測(cè)量的徑向速度,V>0 表示遠(yuǎn)離激光雷達(dá)的正向運(yùn)動(dòng),V<0 表示朝向激光雷達(dá)的負(fù)向運(yùn)動(dòng)。

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

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

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

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

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

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

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

坐標(biāo)旋轉(zhuǎn)法通過(guò)兩次坐標(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ù)的分析。原始信號(hào)經(jīng)過(guò)小波分解得到低頻近似部分和高頻細(xì)節(jié)部分,近似部分反映信號(hào)的整體趨勢(shì),細(xì)節(jié)部分反映信號(hào)在整體趨勢(shì)項(xiàng)上的細(xì)節(jié)波動(dòng)。

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

圖3 10 月15 日激光雷達(dá)60 m 高度處滑動(dòng)平均與小波分解獲得的平均風(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)速脈動(dòng)強(qiáng)度的物理量。垂直方向上湍流強(qiáng)度Iw的物理定義為10 min 時(shí)距內(nèi)的脈動(dòng)風(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)度幾乎沒(méi)有變化,表現(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 湍流動(dòng)能結(jié)果分析

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

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

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

圖5 10 月15 日激光雷達(dá)60 m 高度處垂直方向湍流強(qiáng)度(a)和垂直方向湍流動(dòng)能的小波系數(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 湍流動(dòng)能耗散率結(jié)果分析

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

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

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

圖6 10 月15 日激光雷達(dá)60 m 高度處湍流動(dòng)能與湍流動(dòng)能耗散率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],從湍流能譜密度圖中可以看出湍流動(dòng)能在頻域上的分布特征。

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

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