基于六旋翼無人機搭載風速儀的邊界層風剖面實測研究

李正農，馮 豪，蒲 鷗，沈義俊

(1.湖南大學建筑安全與節能教育部重點實驗室，湖南，長沙410082；2.海南大學土木建筑工程學院，海南，海口570228)

大氣邊界層是大氣層靠近地面約1 km～2 km的薄層，在該區域氣流受地面粗糙度的影響嚴重[1]。目前，在邊界層海拔100 m 以下的低空區域，通常使用測風塔進行定點的風速風向測量[2?3]。測風塔可以提供長時間的穩定測量，但其安裝難度大、成本高，且設塔位置受限，對于多位置測量(如建筑物周圍和風力發電機現場)有一定難度[4]。另外，在邊界層海拔1 km～2 km 的高度內可以利用風廓線雷達進行遠程風測量[5]，但其價格昂貴，且空間分辨率有限[6]。

隨著無人機技術的發展，由于無人機所具備的空間靈活度高、操作簡單、成本低等關鍵優勢已經越來越明顯，研究者們逐漸將無人機投入到測風研究中，或能替代測風塔及風廓線雷達等方式成為新的精確廣泛的測風方法。在1992年，澳大利亞氣象局的Holland 等[7]首次提出了利用搭載皮托管及相關儀器的小型無人機測量風場的方法；2008年，Reuder 等[8]利用小型固定翼無人氣象觀測站(SUMO)對邊界層氣象要素及風速風向進行了實測，然而，使用固定翼無人機進行風場測量時無人機處于水平飛行狀況，無法測量某一地點的垂直風剖面[9]。相較之下，多旋翼無人機憑借可垂直起降和定點懸停的優勢，成為更理想的測量垂直風剖面與獲取固定位置風場參數的平臺。2015年，Marino等[10]提出了一種新的傳感方法，可以確定迎風面的流速大小和方向；2017年，Palomaki 等[11]使用四旋翼無人機的姿態數據間接估算了10 m 高度處的風速和風向并與固定測風塔進行對比，測量誤差分別為風速誤差小于0.5 m/s，風向誤差小于30°；2018年，Prudden 等[12]在四旋翼無人機上的多個方向安裝了多孔壓力探針(MHPP)，通過在一定高度范圍內的實測證實其風速測量的可行性，但是其探測范圍僅局限于90°的迎風面錐形區。

應用六旋翼無人機搭載風速儀進行邊界層風場的風剖面及脈動風特性的測量與研究開展較少，本文利用搭載風速儀的六旋翼無人機進行測量，并以北京延慶某測風塔為參考，通過對比無人機測量結果與測風塔測量結果，來探究無人機測量風場的可行性；通過對數據的分析及擬合，實現無人機對風剖面及脈動風特性的擬合測算；并將實測及分析結果與各國規范[13?17]進行對比，其結論可為無人機測量風場的研究提供參考。

1 試驗概況

1.1 實測場地概況

本試驗在北京市延慶區中科院電工所進行，如圖1所示(實測場地為星標標記處)，實測場地西北側與東南側為高山，西南側較為開闊。圖2給出了實測場地周圍的地貌，可見西面遠處低矮建筑較多，而場地周圍近處基本為空曠農田，按照《建筑結構荷載規范》(GB 50009?2012)[13]有關地面粗糙度的分類標準，該場地為典型的B類場地(地面粗糙度值為0.15)。另外，該實測場地在冬季盛行西風及西北風，風場較穩定，4級、5級風較為常見，試驗條件比較理想。

圖1 北京延慶實測地點地形地貌圖Fig.1 Topographic map of Yanqing sitein Beijing

圖2 北京延慶實測地點衛星圖Fig.2 Satellite map of Yanqing site in Beijing

在圖2中的實測場地星標標記處架立了一個如圖3所示的40 m 高的測風塔，分別在10 m、20 m、30 m、40 m 高度處各安裝有一個風速儀。風速儀正北向安裝，此時風向角定義為0°，正東向來風風向角為90°，按順時針方向依此類推。測風塔的實測結果可以校準與優化無人機所測的風速風向值，以提高實驗結果的精確性。

圖3 北京延慶實測測風塔Fig.3 Wind tower in Yanqing,Beijing

1.2 試驗裝置

測風塔上安裝的風速儀為GILL 三維超聲風速儀，采樣頻率為1 Hz，并在風洞中進行了標定，誤差小于±0.5%。如圖4(a)所示，無人機型號為大疆M600PRO，具有穩定性強、懸停精度高的特點，其飛行垂直誤差±0.5 m，水平誤差±0.5 m，單次航行約30 min。無人機上搭載的風速儀為賽能超聲風速儀，通過改裝可以由無人機直接供電。該風速儀的采樣頻率為1 Hz，與測風塔風速儀相同，并且與氣象部門所使用風速儀的采樣頻率一致。將小型電臺與風速儀連接，實現無線遠程實時傳輸(見圖4(b))。

圖4 實測所用系統Fig.4 System used in the measurement

由于旋翼轉動引起的擾流會對安置在中央架的風速儀測風的準確性產生影響，需要對風速儀的安裝高度進行探究，經過李正農等[18]的研究，將風速儀安裝在機身中央架上方0.53倍旋翼直徑高度處，可以減小槳葉轉動產生的氣流對風速儀測量數據的影響。另外，以上六旋翼的型號與參數選擇及風速儀安裝的具體參數皆只適用此實驗，其普適性暫未深入探究。

2 風特性參數分析方法

2.1 平均風速和風向角

風速風向時程由風速儀直接測得，對風速時程按式(1)分解，得到x、y方向的分量ux(t)、uy(t)(如圖5所示)：

圖5 風速、風向示意圖Fig.5 Schematic diagram of wind speed and direction

式中：ux(t)為N方向分解風速時程；uy(t)為E方向分解風速時程；u(t)為風速時程；θ (t)為風向角時程。

則在測量時距內平均風速U與平均風向角φ可根據式(2)～式(5)計算：

2.2 平均風剖面

在大氣邊界層中，風速隨高度發生變化，該規律可形成平均風速剖面。目前多使用指數率與對數率模型[19?20]來擬合風剖面，本文將采用中國規范使用的指數率模型式(7)對風剖面進行擬合：

式中：z為離地高度；Uz為z高度處的風速；α 為地面粗糙度的參數。

2.3 湍流度

湍流度是對氣流速度脈動程度的一種度量指標，定義為平均時距內脈動風速均方根與平均風速的比值，即：

式中：i指代順風向u及橫風向v；Ii為i風向湍流度；σi為i風向脈動風速均方根；U為平均時距內的平均風速。

2.4 陣風因子

陣風因子是考慮到瞬時風較平均風大而乘的系數，一般定義為陣風持續時間內平均風速的最大值與平均時距內的平均風速的比值，如式(9)所示：

式中：Gu為順風向陣風因子；Gv為橫風向陣風因子；tg為陣風持續時間，本文取tg=3 s；u(tg)為tg內的順風向平均風速；v(tg)為tg內的橫風向平均風速。

2.5 湍流積分尺度

湍流積分尺度是對脈動風旋渦尺寸的度量，反映了湍流空間中兩點脈動風速的相關性。本文根據Taylor 假設自相關函數積分法來計算，即：

式中：Li為i風向湍流積分尺度；Ri(τ)為i風向脈動風速自相關函數。

2.6 脈動風速功率譜密度函數

脈動風速功率譜描述了脈動風速能量隨頻率變化的分布情況，反映了脈動風中不同頻率成分對湍流脈動總動能的貢獻。基于Kolmogorov 理論的統一表達式為：

目前的研究中常使用的順風向脈動風速功率譜經驗譜有Davenport 譜、Kaimal 譜和Karman 譜等[21]。本文也將以這些模型作為標準進行對比，并根據擬合結果給出了Karman 譜表達式：

3 實測結果分析

3.1 平均風特性

3.1.1平均風速和風向角

在實測中，由于無人機的自身抖動、傾角等機身姿態的影響，其搭載風速儀所測數據與測風塔風速儀所測數據存在一定誤差。因此，將無人機定位在測風塔正南方向水平距離5 m 的定點處，在該定點豎向的4個垂直高度(10 m、20 m、30 m、40 m)上進行了風場測量，每次測量時長為10 min，并進行多次實測與對比分析。并結合相關論文的研究，李正農等[18]通過推算傾角進行了風洞試驗，試驗表明無人機所測風速隨機身的傾斜角度增大而增大，且當控制風速為8 m/s時，無人機測得風速相比水平時測得風速增大約3%；當控制風速為10 m/s時，無人機測得風速相比水平時測得風速增大約6%。根據該結論，由于實測中風速不恒定，風場狀態與風洞中有區別，且無人機處于懸停狀態，機身抖動同樣會造成一定影響，因此結合李正農等的研究與實測對比分析對無人機搭載風速儀所測量得到的風速時程進行了如式(13)的修正，同時具體修正系數λ，見表1。

表1 無人機機身修正系數表Table 1 MUA fuselage correction coefficient

式中：u(t)為機身姿態修正后的風速時程；uw(t)為無人機實測風速時程。

表2給出了一組40 m 高度處的風速數據，記錄了無人機機身姿態修正前后的相關風速數據與測風塔數據的對比。與測風塔原始數據進行對比可以發現，無人機機身姿態修正后，平均風速誤差由原始數據誤差的3.04%減小到0.76%，平均風速得到了較好的修正。另外，對比x、y方向的分解風速，所得結果誤差也有所減小，并控制在3%以內；兩者風向角基本吻合，誤差較小。

表2 機身姿態修正后無人機與測風塔風速風向對比Table 2 Comparison of wind speed and direction of drone and wind tower after airframe attitude correction

表3中記錄了10 m～40 m 各高度處經過機身姿態修正后的無人機所測的風速及風向角均值以及測風塔所測的對應值。由表3可見，經過修正后，基于無人機所測的風速與測風塔對應值的誤差均控制在±3%以內，風向角兩者誤差在±0.5%內。由此認為，機身姿態修正能對基于無人機所測的風速起到較好的校準作用，使無人機的測量值更加精確，貼近測風塔值。

表3 不同高度處無人機與測風塔風速風向對比Table 3 Comparison of wind speed and direction between UAV and wind tower at different heights

3.1.2平均風剖面及粗糙度

在實測中由于僅有一臺無人機，單臺無人機無法同時在用一時距內完成10 m～100 m 各個高度風速的測量，而在不同時刻由于風速及湍流度的不同，某一固定高度的風速風向會有所區別，如果利用不同時刻無人機所測得的風速來擬合風剖面將會產生較大誤差，因此需將所有高度處的風速歸一化。每次測量為10 m～100 m 中的10個高度(以10 m 為增量)，這里將每次在10 m 高度測量的時刻稱為t1時刻，并將其余高度及時刻的風速統一以t1時刻為基準進行歸一化處理，從而近似得到在t1時刻10 m～100 m 各個高度處的風速。風速歸一化需要將測風塔40 m 高度作為參考點(40 m 高度的風速風向相對穩定準確)，則在t1時刻，測風塔40 m 高度處參考點風速Uc,40?1即為參考歸一化風速，將此后任意高度測量點的測量時刻稱為t2時刻，則t2時刻z高度處基于無人機所測的風速表示為Uw,z?2，t2時刻測風塔40 m 高度處參考點的風速表示為Uc,40?2，由此，無人機在z高度處相對t1時刻的歸一化風速Uw,z為：

利用歸一化后的10 m～100 m 風速來進行風剖面擬合(如圖6所示為其中一組的擬合結果)，由圖6可見實測值的擬合較好，且其擬合表達式在圖中給出，粗糙度α 為0.147。同時表4中記錄了不同時間各組實測值的擬合相關系數R2值，亦可見各組的擬合精度較高，并且隨著風速的增大，R2值增大，說明風速越大，各高度所測的風速越穩定準確。將實測擬合曲線與中國規范B類地貌，美國規范C類地貌，歐洲規范Ⅱ類地貌及澳大利亞規范B類地貌[13?16]的風剖面曲線進行對比，其中中國規范與美國規范采用指數率模型，歐洲規范及澳大利亞規范采用對數率模型，由圖6可知實測擬合曲線與中國規范B類地貌風剖面最為接近。

圖6 風剖面擬合(B類地貌)Fig.6 Wind profile fitting (Type B geomorphology)

利用指數率模型擬合計算地面粗糙度 α，表4記錄了在不同時間、不同平均風速下擬合得到的地面粗糙度及置信區間，分析所測數據，數據基本穩定在0.137～0.156，范圍為0.019，均值為0.148。從地貌來看，西南西北方向雖然有低矮房屋但是大部分不超過9 m，因此值約為0.15，由此判斷該地貌為B類地貌是合適的。圖7則給出了粗糙度與風速的關系，由圖中趨勢可知，隨著風速的增大，所擬合的地面粗糙度的離散性有所減小，產生的原因是隨著風速的增大，各高度所測的風速更穩定且準確，則計算得到的地面粗糙度結果偏差更小，離散性減小。

圖7 地面粗糙度αFig.7 Ground roughnessα

表4 地面粗糙度擬合Table 4 Ground roughness fitting

3.2 脈動風特性

脈動風特性參數包括湍流度、陣風因子、湍流積分尺度及脈動風速功率譜等，其中基于無人機與測風塔所測的湍流度均由式(8)計算得出，兩者陣風因子根據式(9)計算得到，湍流積分尺度則根據式(10)得出。

表5記錄了10 m～40 m 的4個高度處基于無人機與測風塔所測順風向湍流度均值、順風向陣風因子均值和順風向湍流積分尺度均值的統計結果及兩者對比結果。可知，10 m～40 m 各高度處基于無人機所測的順風向湍流度平均值表現出與測風塔的對應數值近似一致的特性，由表5中測量組中的最大誤差組可以發現各高度上兩者最大誤差控制在±5%以內，由此可見，基于無人機所測的湍流度能夠較為準確真實地反映測風塔所測湍流度。同樣，基于無人機所測順風向陣風因子均值與測風塔所測值也較為接近，且變異系數較小(均小于0.1)，數據較為穩定；基于無人機所測順風向湍流積分尺度均值與測風塔所測均值對比誤差控制在±4%內，最大誤差值在±5%以內。可見無人機能夠較為真實準確地反映測風塔的實測結果，證實了其測量脈動風特性的可行性。

表5 無人機與測風塔風場參數對比Table 5 Comparison of UAV and wind tower wind field parameters

表6給出了無人機在10 m～100 m 范圍內的10個高度處所測得的順風向和橫風向各項湍流參數的均值，同時給出測風塔在10 m～40 m 范圍內4個高度處所測得的對應數據的均值進行比較。由表6可知，隨著高度的增加，無人機所測順風向湍流度Iw,u及橫風向湍流度Iw,v均表現為減小趨勢，另外Iw,u/Iw,v在1∶0.82～1∶0.88，與我國橋梁抗風設計指南規定的1∶0.88 接近，較為準確地反映了良態風的湍流特性。同樣，順風向陣風因子Gw,u及橫風向陣風因子Gw,v也隨著高度的增大而減小；而順風向與橫風向湍流積分尺度則是隨高度增加而增大，且10 m～100 m 范圍內Lw,u/Lw,v在1∶0.68～1∶0.75，這與文獻[22]的比值接近，說明該測點處10 m 范圍內與10 m～100 m 范圍內的湍流積分尺度的比值差異不大。

表6 不同高度脈動風特性參數Table 6 Pulsating wind characteristic parametersat different heights

圖8(a)給出了基于無人機所測順風向湍流度與風速的變化關系(由于篇幅原因只給出順風向參數，橫風向的結論與順風向相同)，可知各高度處順風向湍流度隨著風速的增大呈減小的趨勢，這是由于大氣穩定性隨風速增大而增加，從而減小了湍流度[23]；圖8(b)給出基于無人機所測順風向陣風因子與風速的變化關系，可以看到各高度的陣風因子均隨著風速的增大而減小；圖8(c)則給出了基于無人機所測順風向湍流積分尺度與風速的變化關系，可見隨著風速的增大，各高度的湍流積分尺度都呈現增長趨勢。

圖8 順風向各項湍流參數隨風速的變化Fig.8 Changesof turbulence parameterswith wind speed

3.2.1湍流度與陣風因子

湍流度與陣風因子之間存在一定相關性，一般認為湍流度越大，陣風因子也越大。多年來，人們對湍流強度與陣風因子的關系進行了研究，Ishizaki[24]和Choi[25]分別提出的公式可以組合成一個方程式(15)來表示湍流強度和陣風因子之間的關系：

式中：T為平均時距；tg為陣風持續時間；k1、k2為無量綱參數，其中Ishizaki[24]建議k1=0.5，k2=1.0，Choi[25]建議k1=0.62，k2=1.27。圖9給出了40 m高度處湍流度與陣風因子的相關性，利用式(15)對基于無人機所測的各高度處的數據進行擬合，擬合結果與Ishizaki和Choi的經驗表達式一起給出，所得k1=0.49，k2=0.98，并且在其他高度上，10 m 高度處k1=0.51，k2=0.99；20 m 高度處k1=0.53，k2=1.00；30 m 高度處k1=0.49，k2=0.95，這與Ishizak 擬合的曲線非常接近，說明陣風因子隨湍流強度呈近似線性增長的趨勢。

圖9 湍流度與陣風因子相關曲線Fig.9 Turbulence and gust factor correlation curve

3.2.2湍流度剖面及湍流積分尺度剖面

由于無法在同一時間測得各個高度的湍流度及湍流積分尺度，因此湍流度剖面及湍流積分尺度剖面的各高度均值與平均風速剖面中的計算相似，為了反映同一時間各高度處的湍流度及湍流積分尺度情況，利用式(16)和式(17)將所得值進行轉換：

式中：Iw,z為最終所得同一時間下各高度處無人機所測湍流度；Ic, 40?1為t1時刻測風塔40 m參考點湍流度；Iw,z?2為t2時刻z高度處無人機測得的湍流度；Ic,40?2為t2時刻參考點的湍流度。同理，式中Lw,z為最終歸一化的湍流積分尺度，Lc, 40?1為t1時刻測風塔參考點湍流積分尺度，Lw,z?2為t2時刻z高度處無人機測得的湍流積分尺度，Lc,40?2為t2時刻參考點的湍流積分尺度。

圖10給出了基于無人機所測的順風向湍流度均值隨高度的變化關系及湍流度剖面的擬合。由圖8可知湍流度與風速相關，隨著風速增大湍流度呈減小趨勢，因此根據圖8的分布將風速分成3個風速段來進行湍流度剖面的擬合，分別是4 m/s～6 m/s風速段、6 m/s～8 m/s風速段、8 m/s～10 m/s風速段，結果在圖10中呈現。圖10中還給出了同地貌條件(B類地貌)及高度下的中國、美國、日本及歐洲規范中的湍流度剖面，由圖可知各風速段順風向湍流度均值的擬合較好，而擬合曲線與各國規范都有差異，較中國規范偏大，而相對于美國、日本及歐洲規范數值偏小。通過式(16)換算后的湍流度變異系數較大，經分析，由于不同高度的脈動風特性較平均風特性的差異性更大，即t2時刻Iw，z-2與Ic，40-2相關性較低，因此在利用式(16)時Iw，z-2/Ic，40-2比值關系將產生一定誤差，導致得到的歸一化湍流度也會產生誤差。

圖10 順風向湍流度隨高度變化Fig.10 Downwind turbulence variation with altitude

圖11(a)給出了基于無人機所測的順風向湍流積分尺度均值與高度的變化關系。同樣將風速分成3段進行擬合，并與同地貌條件下(B類地貌)各國規范進行對比。由圖可知湍流積分尺度的數值大小相比于各國規范都要小，初步分析可能原因是測量時距較短，而根據Yu 等[26]的研究，隨著實測時距的減小，湍流積分尺度會減小，并存在如式(18)的關系：

式中：T表示測量時距。將實測數據按式(18)轉換為10 min 時距相對值，得到圖11(b)中的修正后的湍流積分尺度剖面，發現在修正后不同風速段的湍流積分尺度有一定差異，總體均值可以認為與日本規范較為接近。誤差的產生可能是由于時距對湍流積分尺度的影響仍然存在，且式(17)存在與式(16)相似的計算誤差，另外場地與風場的不同對湍流積分尺度的修正也有較大影響。

圖11 順風向湍流積分尺度隨高度變化Fig.11 Turbulence integral scale variation with height

3.2.3脈動風速功率譜密度函數

脈動風速功率譜密度函數根據式(11)計算，其中Davenport 譜不隨高度變化，而Kaimal譜和Karman 譜與高度相關。本文選取了多組10 m～40 m 無人機與測風塔在時距為10 min 下測量得到的強風時程樣本進行脈動風速功率譜的擬合，得到了4個高度處順風向脈動風速的歸一化功率譜。由于Kaimal 譜和Karman 譜與高度相關，因此將測風塔10 m～40 m 各高度實測值的順風向脈動風速功率譜依次單獨繪出并與經驗譜進行比較如圖12(a)～圖12(d)所示，由圖可見，10 m 高度處實測值的順風向脈動風速功率譜由于湍流較大，與各經驗譜的擬合都相對較差；而20 m～40 m 各高度實測值的順風向脈動風速功率譜與Karman 譜的擬合較好。

圖13則給出了在30 m 及40 m 高度處由無人機與測風塔實測值擬合的順風向脈動風速功率譜密度函數，由圖可知，無人機所擬合的結果與同高度及時距下的測風塔的擬合結果總體上較為吻合，在部分頻率段下稍有差異，高頻段無人機擬合的功率譜密度函數波動稍大但仍與測風塔擬合的函數有一致的趨勢，再將無人機與測風塔所擬合的脈動風速功率譜與經驗譜進行對比，同樣發現實測值擬合的順風向脈動風速功率譜與Karman譜有較高的一致性。圖14則給出了40 m 以上的4個高度上基于無人機所測得的順風向脈動風速功率譜密度函數，由圖可知各高度處實測值擬合的順風向脈動風速功率譜均與Karman 譜的一致性最高，這也說明基于無人機實測所擬合的各高度處的脈動風速功率譜密度函數有較高的真實性與準確性，但是無人機仍無法像測風塔在圖12(a)中一樣，在同一時間反映不同高度處的結果。

圖12 10 m～40 m 測風塔順風向脈動風速功率譜Fig.12 Measured power spectra of fluctuating wind speed along wind direction of 10 m～40 m wind tower

圖13 無人機及測風塔順風向脈動風速功率譜Fig.13 Measured power spectra of fluctuating wind speed along wind direction of UAV and wind tower

圖14 50 m～100 m 無人機順風向脈動風速功率譜Fig.14 Measured power spectra of fluctuating wind speed along wind direction of 50 m～100 m wind tower

4 結論

本文對比了北京延慶實驗基地的測風塔風速儀實測數據及六旋翼無人機搭載風速儀的實測數據，對多旋翼無人機的邊界層風場實測的可行性進行了研究，經分析得到以下結論：

(1)無人機實測平均風速與測風塔的實測平均風速存在一定差異，但通過對無人機所測風速時程進行修正，可得到與測風塔所測結果相近的數據，且兩者誤差較小(±3%以內)。

(2)利用無人機所測各高度處歸一化后的平均風速擬合的風剖面擬合精度較高，數據離散小(R2>0.87)，且與我國規范B類地貌的風剖面接近；擬合得到的地面粗糙度均值為0.148，接近規范標準值0.15，并判斷該地為B類地貌。

(3)由無人機風場實測數據計算所得的湍流度及陣風因子與由測風塔實測數據計算所得結果接近，誤差均在±5%以內，且隨著高度或風速的增大，湍流度及陣風因子均逐漸減小；基于無人機所測得的湍流度剖面受時距及換算誤差的影響與各國規范均有差異，較中國規范偏大，而相對于美國、日本及歐洲規范數值偏小；湍流度和陣風因子之間存在指數關系，陣風因子隨湍流度的增大而增大，擬合結果與Ishizak 擬合曲線接近。

(4)由無人機風場實測數據計算所得的湍流積分尺度與由測風塔實測數據計算所得結果接近，并與風速、高度均呈正相關。各風速段的湍流積分尺度隨高度的變化關系有一定差異，總體與日本規范較為接近。

(5)由無人機風場實測數據擬合的脈動風速功率譜與同高度及時距下測風塔實測數據擬合的脈動風速功率譜基本一致，在高頻部分稍有差異，且10 m～100 m 各高度所測脈動風速功率譜與經驗譜Karman 譜均較吻合。