基于无人机的复杂地貌上空风场实测研究

2021-05-18吴红华胡昊辉李正农

实验流体力学 2021年2期

吴红华, 张亮, 冯豪, 胡昊辉, 李正农

湖南大学建筑安全与节能教育部重点实验室, 长沙 410082

0 引言

自然界中，近地风在其行程中会遇到不同的地形地貌，随着地形和地表粗糙度的变化，风场变化也不尽一致[1]，目前对于近地风特性的研究主要是风场沿高度方向的梯度分布和风剖面研究，而对地形地貌影响下风场的平面分布特征研究较少，并且多集中于风电领域。陈爱等[2]简化三维分析模型，用数值模拟方法研究了复杂地形影响下不同来流风向时的风场分布，确定风力发电机的最佳安装位置。Ren等[3]在某6 km×6 km的试验场内设立10个固定测点，并结合CFD模拟和测点相关性，建立了风场预测模型，并将预测值与实测值进行对比，结果表明：弱相关测点预测效果较差。Palma等[4]采用线性与非线性方法研究马德拉岛东部海岸某区域风场特性，并与实测数据对比，得到了目标区域的风场参数，用以指导复杂地形环境下风资源评估。雷若冰等[5]对多个测风塔风速数据进行相关性分析，有效地描述了整个风电场群的风速分布。现有研究多基于数值模拟和多个测风塔之间数据相关性的算法研究，当测风塔相距较远时，对局部区域的风场情况难以全面掌握，有一定局限性。

传统的风场实测主要依靠固定式测风塔[6-7]和现场搭设可移动式测风桅杆[8-10]开展，但这类设备体积大、可移动性差，往往只能测量风场沿高度方向的梯度分布，测量平面风场难度较大。随着无人机技术的发展，利用小型无人机测风引起广泛关注。1992年，澳大利亚气象局Holland等首次提出以气象无人机携带皮托静压管等仪器测风的方法[11]。马舒庆等[12]基于微型气象无人机盘旋飞行，利用水平空速归零法和解析法求得风速风向。周伟静等[13]利用固定翼无人机携带皮托-静压管等测风仪器进行风场测量，并对实测数据进行误差修正，提高了测量精度。固定翼无人机操控难度大，难以实现定点测量，且需要特定的起降跑道，对场地要求较高。多旋翼无人机具有垂直起降、机动灵活、定点悬停、精确定位等诸多优势，因此用多旋翼无人机搭载风速仪进行风场测量值得深入探讨。Bruschi等[14]在四旋翼飞行器机身上方安装二维风速传感器并进行风洞试验，结果表明：风速数据存在一定误差，而风向角的测量几乎不受影响，但该试验未对机身姿态改变时的影响进行研究。李正农等[15]对六旋翼无人机搭载风速仪测风的准确性进行了风洞试验研究，结果表明：机身倾角会对风速测量产生一定影响，对风向角测量无影响，但试验时无人机机身固定，不能模拟真实的飞行状态。现阶段对于多旋翼无人机测风的研究多基于数值模拟和风洞试验，这两种研究方式并不能模拟无人机在复杂大气环境中真实的飞行状态，因此，需要对无人机搭载测风仪实测的可行性进行深入研究。

本文通过多旋翼无人机搭载风速仪进行实地风速和风向测量，并对实测数据进行分析及误差修正，初步验证在复杂环境中利用无人机搭载风速仪进行风速风向实测的可行性；其次，基于不同地形地貌对于近地风场的影响不同，通过无人机多点测风并结合测风塔数据，获取各测点风场与测风塔数据的相互关系，从而在已知测风塔数据时，推知测风塔周围不同地貌影响下的某一高度平面的风场特性。不同来流方向地面粗糙度不同，对风场造成的影响也不尽相同，因此本文选取实测期间某一特定来流方向的数据进行分析，旨在对多旋翼无人机进行风场测量这一新方法进行初步探讨和研究。

1 实测概况

1.1 场地概况

实测场地位于北京市延庆区中国科学院某试验基地内，如图1所示。

图1 实测地点Fig.1 Test site

红色星标位置为测风塔。为便于描述，将场地划分为4个象限平面。如图所示，第一象限地貌较为复杂，靠近测风塔处有2组东西走向的槽式聚光镜(高约7 m)，较远处为低矮松树(2 m左右)和灌木丛；第二象限中树木较多，靠近测风塔处为低矮松树(2 m左右)，距离测风塔较远处有厂房和大片高大的杨树林；第三象限相对平坦，多为低矮松树(2 m左右)和灌木丛，其对上空风场的影响可能较小；第四象限地貌最为复杂，有2组南北走向的槽式聚光镜(高约7 m)和定日镜群组(单个约12 m高)，还有一些试验厂房。复杂的场地状况可能会造成上空风场的波动更大。该场地冬季西风盛行，来风方向场地较为平坦空旷，4、5级风较为常见，是进行实测研究的理想场所。

1.2 实测系统

实测系统由40 m高的测风塔和大疆M600 PRO六旋翼无人机组成，如图2所示。

图2 现场实测系统Fig.2 The field measurement system

在测风塔10、20、30和40 m高度处各安装了一个WindSonic二维超声风速仪，风速仪采样频率为1 Hz。无人机机身上部搭载SA210二维超声风速仪，采样频率也为1 Hz。风速仪测得数据通过无线电台实时传输到电脑。两超声波风速仪的具体参数见表1。测风塔风速仪和无人机风速仪均为正北向安装，规定正北方向来风时，风向角θ为0°；正东方向来风时，风向角θ为90°，其余风向角按照顺时针方向类推。

表1 超声波风速仪参数Table 1 Parameters of ultrasonic anemometer

为保证风速仪测量结果的准确性，实测前在风洞中利用澳大利亚TFI Series 100眼镜蛇三维脉动风速探头对两超声风速仪进行二次标定(出厂时已标定)。标定时，将眼镜蛇风速探头分别与两超声波风速仪在同一高度同一位置处进行测量，标定结果见表2。由表2可知，两超声风速仪与眼镜蛇风速探头所测风速基本一致，误差(相对于眼镜蛇风速探头数据)均小于0.5%，可以认为两超声风速仪测风数据准确，可以用其所测数据直接进行对比分析。

表2 风速标定结果Table 2 Wind speed calibration results

1.3 实测方案

本次试验分2阶段进行。第一阶段探究复杂环境中无人机搭载风速仪实测的可行性与准确性:为避免无人机与测风塔相互影响，保证无人机飞行安全，将无人机放置在与测风塔水平距离为5 m处垂直起飞，与风场实测相对应，在20 m高度处进行风速风向测量，每次测量时间10 min，将测得的无人机数据与测风塔数据进行对比分析。第二阶段进行平面风场的测定：通过无人机多点测风，并结合测风塔数据作为参考，获取各测点位置风场与测风塔数据的相互关系，推知测风塔周边不同地貌影响下的平面风场特性。

由于地貌中定日镜约为12 m高，综合考虑地貌对风场可能的影响以及无人机飞行安全，此次风场实测高度选为20 m。图3中红点位置为测风塔，将其作为参考点，设置为原点(0，0)；其他测点位置通过坐标x、y确定，x指示东西方向，东方为正，y指示南北方向，北方为正，x、y的范围均为[-100，100]。测点x、y方向间隔均为20 m，总共设置11×11=121个测点。无人机在每个测点悬停10 min。

图3 实测测点坐标系Fig.3 Coordinate system of measured points

2 无人机数据前处理

2.1 数据前处理基本理论

矢量风速与风向是密切相关的，风向不同时直接比较风速意义不大。后续风场测量过程中，受地貌影响，不同测点处无人机所测风向与测风塔风向不可能完全一致，因此需要将风速分解到x、y方向上，以便于比较。测风塔和无人机上安装的风速仪所实测到的数据为风速时程u(t)和风向时程φ(t)两列，规定正北方向(y轴正向)为0°，正东方向(x轴正向)为90°，如图4。通过风向φ(t)可以将风速时程u(t)分解为分量ux(t)和uy(t)。实测风速时程沿x、y方向的分量ux(t)和uy(t)为：

图4 风速、风向示意图Fig.4 Sketch map of wind speeds and directions

(1)

为确定最大平均风速而规定的时间间隔称为平均时距。我国规范将平均时距取为10 min，平均时距内2正交方向的平均风速为：

(2)

式中，n为采样频率与时距的乘积。相应平均时距内水平平均风速U和平均风向θ分别为：

(3)

(4)

定义Ii(i=x，y)为x、y方向上的湍流度分量。Ii反映了风的脉动程度，为某一高度处在给定持续时间内风速分量时程标准差与相应高度处水平平均风速U的比值:

Ii=σi/U, (i=x、y)

(5)

式中，σx、σy分别为ux(t)、uy(t) 的标准差。

2.2 无人机实测风速风向修正

风洞试验结果表明，无人机机身倾角会对无人机风速测量产生影响，因此需要先根据风洞试验结果对无人机倾角造成的误差进行修正。将无人机实测风速时程数据u(t)按照每30 s时距求取平均值，根据平均值的大小选取修正系数，进行不同倾角下风速修正，具体修正公式和修正系数为：

u1(t)=u(t)/α

(6)

式中，u1(t)为修正后风速时程，u(t)为原始风速时程，α为修正系数。李正农等[15]根据无人机相关参数推算出不同风速下机身倾角，并进行风洞试验，试验结果表明：机身倾斜时无人机所测风速偏大，控制风速为8 m/s时，无人机测得风速比水平时约增大2%；控制风速为10 m/s时，无人机测得风速比水平时约增大5%，其余风速范围按照插值法确定，具体修正系数如表3所示。

表3 风速修正系数Table 3 Correction coefficient of wind speed

利用上述修正方法对实测无人机风速时程进行修正，然后按式(1)～(5)分解计算，将计算结果与测风塔数据进行对比，如表4所示。可以看出，无人机原始实测风速较测风塔偏大，通过修正，可以很大程度上消除机身倾角对风速测量的影响。修正后的无人机平均风速数据与测风塔数据比较接近，但是由于风向角存在误差，分解到x、y两个方向的风速分量仍然存在较大误差。

表4 无人机风速修正前后数据对比Table 4 Data comparison before and after MUA wind speed correction

风洞试验结果还表明，无人机风向的测量不会受到机身姿态变化的影响，但通过实测数据可以看出无人机和测风塔存在大约3.5°的风向测量误差。这是由于测风塔和无人机风速仪均以正北方向为0°方向角，在实际操作过程中，无人机风速仪0°角使用GPS和电子指南针辅助对北，精度高，几乎不会产生误差；而测风塔风速仪安装在20 m高度处，高空作业难度大，在手动进行0°角对北的过程中不可避免地可能产生一定的偏差，造成此处3.5°的安装误差。风向误差的存在导致分解到x、y两个方向的风速分量误差较大，因此以测风塔数据为基准，修正无人机风向数据，可以消除风向误差对风速的影响，从而更好地探究无人机实测的准确性。将无人机所测风向时程每一瞬时点风向修正3.5°以消除安装误差，然后利用修正后的风向时程进行计算，结果如表5所示。

从表5可以看出，经过风向修正后，x方向风速分量增大，y方向风速分量减小，2个风速分量和平均风速都更接近测风塔风速数据，风速误差显著减小。在后续风场测量中，为保证无人机数据的准确性，将无人机起飞点统一设定于距离测风塔5 m处，利用测风塔数据修正无人机风速仪误差后，再进行不同测点处的风场测量。

表5 无人机风向修正前后数据对比Table 5 Data comparison before and after MUA wind direction correction

数据处理前后的风速和风向时程如图5所示，从图中可以看出，修正后无人机风速和风向时程数据与测风塔风速风向时程数据吻合良好。需要注意的是，y方向风速分量时程原本为负值，为便于绘图，将y分量同时乘以-1处理，图6同理。

图5 测风塔、修正后无人机风速风向时程图Fig.5 Wind speed and direction time history of wind tower and modified UAV

2.3 无人机湍流度修正

修正后的无人机风向与测风塔风向一致，无人机平均风速、风速分量与测风塔数据相比误差较小，但是2个分量方向的无人机湍流度仍然较大，原因主要有2点：一是当风速发生变化时，无人机为保持飞行的稳定性会立刻进行横滚角和俯仰角的姿态调整，造成部分瞬时点误差偏大，增大了数据的波动程度；二是无人机在实测过程中由于信号传输等原因产生少量野值点，也使测得的数据波动变大，因此需要进一步处理计算湍流度误差。

(i=3,…,n-2)

(7)

式(7)无法计算4个端点(开始2个点和最后2个点)，这种情况称为“端部效应”。4个端点值可以采用原始数据补齐，即：

(8)

将上述经过无人机风速和风向修正后的数据进行分解，然后对分解后的x、y方向风速分量时程进行5点滑动平均处理，结果如表6所示。从表6中可以看出，滑动平均对平均风速、平均风向、风速分量没有影响，但是可以有效地减小湍流度。其原因是采用5点滑动平均后，部分野值点的影响被消除，数据更为平滑，时程数据的波动程度明显减小。

表6 无人机数据滑动平均修正前后结果对比Table 6 Comparison of results before and after MUA data moving average correction

对风速分量滑动平均之后，风速分量时程的每个瞬时点会发生变化，因此与之对应的风向时程瞬时数据也发生了改变，利用式(9)可以推算变化之后的风向时程数据。

或

(9)

图6 测风塔、滑动平均后无人机风速风向时程图Fig.6 Wind speed and wind direction time history of wind tower and UAV after moving average

3 实测风场数据分析

由于地貌的影响，不同测点位置的风向数据与测风塔存在不同，风场数据无法直接进行比较，需要将实测数据进行分解处理。通过第2节所述无人机数据前处理方法对各测点无人机实测风速风向数据进行修正分解，然后取对应相同时间段的测风塔分量数据，利用式(10)求取各测点位置无人机和测风塔分量数据的比值C，比值C表示地貌影响下20 m高度水平面不同测点位置风场参数与测风塔风场参数的相对关系。通过比值C和测风塔风场数据可以推知不同测点在地貌影响下的风场状况；此外，比值C能在一定程度上反映地貌对测点风场的影响程度，C越大表示测点位置地貌对风场的影响程度越大。比值C计算公式为：

(10)

3.1 风场风速结果分析

实测后得到全部测点的风速比CS见表7和8，为了更加直观地得到整个风场平面风速实测值的变化趋势，通过Matlab对所有实测值风速比绘制三维图和等值线图，并将等值线图与场地图进行对应，详见图7和8。

从表7、8和图7、8可以看出，平面内风速比有一定的波动，这种波动情况能够反映不同地貌影响下风场风速变化情况。通过x、y方向风速分量比三维图可以看出，第二、四象限风速比波动较大，第一、三象限波动较小。需要注意的是实测风向接近270°，分解后y方向风速分量较小，导致部分比值大，而x方向更接近实测风的来流方向，更具有代表性。

表7 x方向风速分量比CSxTable 7 The ratio of wind speed component in x direction CSx

图7 风速比CS三维图Fig.7 Three-dimensional graph of wind speed ratio CS

第一象限中地势平坦，风场受地貌影响较小，因此风速比有一定变化但是变化不剧烈。

第二象限中，距测风塔远处的风场出现了较大的风速比波动，此象限x、y方向风速分量比最大值均位于(-60，80)测点，分别为1.045、3.540，原因可能是受到附近成片杨树林的干扰，风速发生较大变化。

第三象限的地貌较为平坦，当西侧来风时，整个场地的树木对于20 m高度处的风场影响较小，故此象限内风速比的整体变化较小，仅有个别测点数值较大。

图8 风速比CS等值线图Fig.8 Contour graph of wind speed ratio CS

表8 y方向风速分量比CSyTable 8 The ratio of wind speed component in y direction CSy

第四象限地貌最为复杂，西侧来风时，由于地面槽式聚光镜和定日镜群的干扰，此象限出现了较大的风速比波动。整个风场中x、y方向风速分量比最大值均出现在此象限内，其中x方向风速分量比最大值1.098出现在两层试验用房的上方，y方向风速分量比最大值3.738出现在东西走向的槽式聚光镜上方，另外此象限内还有多个较大值出现，说明这些测点位置风速较测风塔变化较多，风场受地貌影响显著。

3.2 风场风向角分析

全部测点的风向比CD见表9。为了更加直观地得到整个风场平面风向角比的变化趋势，通过Matlab软件将所有风向角比绘制成三维图和等值线图，详见图9和10。

表9 风向比 CDTable 9 The ratio of wind direction CD

图9 风向角比值CD三维图Fig.9 Three-dimensional graph of wind direction ratio CD

图10 风向角比值CD等值线图Fig.10 Contour graph of wind direction ratio CD

测风塔风速仪和无人机风速仪均以正北向为0°风向角，试验过程中测得的来流风为西风，来流方向场地较为平坦空旷。通过分析实测数据可知，无人机与测风塔风向数据的比值非常接近1，即在相同时间段内，无人机风速仪和测风塔风速仪测得的风向角差别较小。

从风向角比值三维图可以看出，风向角比在第四象限内的波动有明显的增大，最大值1.015和最小值0.972均位于这一象限。从风向角比等值线图可以看出，风向角比波动较大的区域位于聚光镜和定日镜群上方，原因可能是槽式聚光镜和定日镜群较高，导致部分测点风场受到干扰，风向产生变化，从而增大了风向角比的波动程度；其他3个象限风向角比波动较小，仅有个别测点数据较大，说明这3个象限内地貌对20 m高度处风场风向影响有限。

3.3 风场湍流度分析

全部测点实测后计算得到的x、y方向湍流度分量比CT见表10和11。与风速和风向相同，通过Matlab软件把所有湍流度比绘制成三维图和等值线图，并将等值线图与场地图进行对应，详见图11和12。

从表10、11和图11、12可看出，x、y方向湍流度分量比的波动程度大，这说明整个风场平面湍流度变化较大。从湍流度比等值线图可以看出：

图11 计算湍流度比CT三维图Fig.11 Three-dimensional graph of calculated turbulence ratio CT

表10 x方向湍流度分量比CTxTable 10 The ratio of turbulence component in x direction CTx

表11 y方向湍流度分量比CTyTable 11 The ratio of turbulence component in y direction CTy

第一象限中部分测点风场受到下方东西走向的聚光镜的干扰，x、y方向湍流度分量比数据偏大。

第二象限远离测风塔的位置受到成片杨树的影响，湍流度有所增加，因此部分测点的湍流度分量比值较大，此象限内x方向湍流度分量比最大值达到1.365。

第三象限除个别测点湍流度分量比偏大外，整个象限湍流度比及其变化程度与其他3个象限相比较小，这是因为第三象限地貌较为平坦，地貌对上空风场影响小，风速的波动程度也更小。

结合比值表和图11、12可以看出，整个风场平面x方向湍流度分量比最大值1.431和最小值0.540均位于第四象限东西走向聚光镜与南北走向聚光镜的交汇区域，另有多个较大值出现在定日镜群上方，y方向湍流度分量比在槽式聚光镜上方也明显偏大；此外，第四象限等值线与其他3个象限相比更为密集，表明第四象限x、y方向湍流度分量比波动相对其他象限更为剧烈，原因是地面的7 m高槽式聚光镜和12 m高定日镜群对20 m高度处风场平面干扰更大，使得上空风场变化复杂，湍流度显著增大。