六旋翼无人机测风系统实测海岛地区风剖面

2023-06-03黄斌李昊董金爽王相军

湖南大学学报（自然科学版） 2023年5期

黄斌，李昊，董金爽，王相军

（1.海南大学土木建筑工程学院，海南海口 570228；2.扬州大学电气与能源动力工程学院，江苏扬州 225127）

风受地面粗糙度影响的近地面大气层称为“大气边界层”，在大气边界层内风速随高度和地面粗糙度而变化［1］.强风天气下，各类建筑物和基础设施的风致破坏［2］、树木和农作物的倒伏［3-4］，已给人类生命财产和社会经济造成严重损失.为减少风灾的带来的损失，开展近地面风场实测研究成为首要任务.此外，大气边界层风场测量对风环境的预测和评估（城市的街区风环境和污染物扩散等）［5］、风能资源利用（风电场选址）［6］也具有重要的理论和工程意义.

传统测风方法主要利用固定式测风塔［7］、测风桅杆［8］、测风雷达［9］等进行风场测量，但该方法安装难度大、移动困难且价格昂贵.随着无人机技术的发展，利用固定翼和多旋翼无人机搭载测风设备进行测风引起了学者的关注.Holland 等［10］在1992 年提出利用无人机搭载气象传感器对难以到达地区进行气象要素观测的方法.Reuder 等［11］在小型固定翼无人机上配备传感器，组成气象观测站，实现了对大气边界层的温度、湿度、压强、风速和风向等要素的测量，但该方法难以实现定点观测.王彦杰等［12］和周树道等［13］针对固定翼无人机转动探测、迎角和航迹倾斜角对测风精确性的影响，提出了改善测量精度的测风模型，并通过数值仿真验证了模型的有效性.Borup 等［14］提出一种基于固定翼无人机的指数稳定非线性风速观测器，能够估计风速并计算相对速度，并通过数值模拟验证了其测风的准确性.虽然固定翼无人机续航时间长、速度快、测量高度大，但是起降时对场地要求高，操控难度大，灵活性较差，不能定点悬停，无法测得特定地点的垂直风剖面［15］，使得其在测风领域的应用受到一定限制.

相比而言，多旋翼无人机具有操作简单、机动灵活、垂直起降、精确定位、定点悬停、成本较低等诸多优势，在边界层测风领域的应用值得深入研究.Moyano［16］采用3 种算法将四旋翼无人机的倾角与局部风速和风向联系起来，推导出风剖面，并将结果与安装在桅杆上的风速仪实测值进行比较.Bruschi等［17］利用四旋翼无人机搭载二维风速仪在风洞中进行风速和风向测量，发现旋翼转动只对低于 10 m/s的风速有显著影响，对风向基本无影响.Prudden等［18］在小型四旋翼无人机的旋翼前方约4.5 倍旋翼直径处安装风速探头，测量了特定高度的风速和湍流强度，验证了其测风的可行性；但是测量范围仅限于机头前方90°锥形区域内的来流.李正农等［19］通过风洞试验研究了搭载超声风速仪的六旋翼无人机旋翼转动对测风精确性的影响，指出旋翼转动会影响无人机中心上方风场，风速误差随高度先增大后减小，对风向角测量没有影响.吴红华等［20］采用测风塔和一台搭载超声风速仪的六旋翼无人机，对复杂地貌下平面风场进行实测，通过对比无人机与测风塔的测量结果，分析实测场地地貌对上空风场的影响.李正农等［21］利用测风塔和一台搭载超声风速仪的六旋翼无人机进行边界层风场实测，将无人机实测风场参数与测风塔实测值对比，通过风速归一化拟合得到来流风剖面.总体而言，现有研究基于风洞试验和测风塔实测，初步验证了多旋翼无人机搭载测风仪进行风场测量的可行性.然而，应用多旋翼无人机代替测风塔，对复杂地貌的风剖面和湍流强度进行实测的研究仍很少.

本文采用两台搭载超声风速仪的六旋翼无人机测风系统对海岛地区近海岸、农业温室基地以及密集建筑群等典型地貌的风场进行现场实测.基于风洞试验得到的修正系数，对无人机实测数据进行修正；通过归一化计算，拟合出多类典型地貌的来流风剖面和湍流强度剖面；分析不同测点上空风场受地表粗糙度的影响，并将实测结果与相关风荷载规范和标准进行对比；在验证无人机测风可靠性的基础上，获得热带海岛地区多类典型地貌的风场参数.

1 实测概况

1.1 实测场地

如图1 所示，实测场地位于海口市海甸岛，实测区域包括一处近海岸、一处农业温室基地和两处密集建筑群，分别代表几类较为典型的地貌类别，以研究地面粗糙度对风场的影响.

图1 实测场地Fig.1 Test site

实测海岸为海甸岛西海岸，西面临海，东面为开阔平坦空地，选取风向较稳定（变化幅度在7°内）的偏西风对海面风场进行实测.依据《建筑结构荷载规范》（GB 50009—2012）［22］对地面粗糙度类别的分类标准，实测场地为典型A类场地.

实测农业温室基地位于海南大学海甸校区内，如图2 所示，基地尺寸为175 m×125 m，分布有密集温室群，各类温室高度为5～7 m，呈带状分布，相邻温室之间相距2～3 m.依据《建筑结构荷载规范》（GB 50009—2012）［22］对地面粗糙度类别的分类标准，大致属于B 类场地.基地东侧为田径场；北部为海洋学院五层办公楼和高为5～7 m的树林；西侧为池塘，地貌空旷平坦；南部零星分布有3 m高的工作间、2 m 高的果树和0.5 m 高的蔬菜地.实测期间风向较稳定，为偏东风（变化幅度在10°内）.

图2 农业温室基地及周围环境Fig.2 Agricultural greenhouse base and surrounding environment

考虑实测期间主导风向为偏东风（风向较稳定，变化幅度在7°内），选取以下两处密集建筑群进行实测：一是密集温室群以北约200 m 处，高层建筑农科楼以东的区域；二是密集温室群东南方向约460 m处，高层建筑研发楼以东的区域.如图3 所示，农科楼和研发楼以东2 km 范围内的建筑物分布较为密集，且基本为多层和高层建筑.依据《建筑结构荷载规范》（GB 50009—2012）［22］对地面粗糙度类别的分类标准，大致属于C类或以上场地.

图3 农科楼和研发楼以东地貌Fig.3 The landscape on the east side of agricul science building and R&D building

1.2 实测系统

如图4 所示，实测系统由两台型号为 DJI M600PRO 的六旋翼无人机、两台超声风速仪、无线电台和地面电脑接收端组成，能够进行风数据的实时采集、传输和存储.无人机飞行时的外形尺寸为1 668 mm×1 518 mm×727 mm，载重为6 kg，续航时间约为30 min，能够抵抗7 级风力，悬停精度高（垂直±0.5 m，水平±1.5 m），稳定性强，这些均为无人机搭载风速仪进行测风提供了保障.

图4 实测系统Fig.4 The measurement system

为尽可能减小无人机旋翼转动对风速仪测量造成的干扰，根据李正农等［19］的试验结果，在无人机中心上方20.2 cm（0.53 倍旋翼直径）处安装风速仪，且风速仪与旋翼中心水平距离为53 cm.规定正北方向来风时，风向角为0°，正东方向来风时，风向角为90°，按顺时针方向类推.无人机A 和B 分别搭载SA210、FT205 二维超声风速仪，采样频率分别为1 Hz、10 Hz，风速仪的相关参数见表1.

表1 风速仪相关参数Tab.1 Related parameters of anemometer

为确保风速仪测风的准确性，在实测前，通过风洞试验对两台超声风速仪进行二次标定（出厂已标定）.如图5 所示，将TFI Series 100 眼镜蛇探头和超声风速仪分别固定在同一高度，选取4 m/s、6 m/s、8 m/s和10 m/s共4个风洞控制风速进行测量，由表2可知，两台超声风速仪与眼镜蛇探头所测风速值吻合较好，误差（与眼镜蛇探头所得结果相比）均小于0.5%，认为两台超声风速仪的测量数据准确.

表2 风速标定结果Tab.2 Results of wind speed calibration

图5 标定风速仪的风洞试验Fig.5 Wind tunnel test of calibrated anemometer

1.3 实测方案

采用双无人机测风的基本原理为：在相同风向下，无人机A 悬停于某一高度作为参考点，无人机B移动测量其他各个高度处的风数据，通过对同一时距内两台无人机的实测结果进行归一化处理，近似得到相同时刻的风剖面［21］.为探究双无人机实测系统对特定区域风场测量的可行性，选取海岛地区几类具有典型地面粗糙度的地貌进行测风研究.实测高度范围均为近地面10～120 m，以10 m 为间距，总共12 个高度，各高度采集10 min 数据.实测时，无人机A 悬停于100 m 处作为参考点，无人机B 移动测量.为保证实测结果可靠性，选取风向稳定阶段实测.在海岸实测时，选取从海面吹向海岸的偏西风对海面风场进行实测.农业温室基地和两处密集建筑群在实测期间主导风向均为偏东风.如图2 所示，为细致考虑地面粗糙度对近地面风场的影响，在温室基地均匀布置了9 个测点，分别位于基地迎风上方、中部和尾部.在实测农科楼和研发楼东侧密集建筑群的来流风剖面时，应使无人机从楼前空地起飞以避免其与高层建筑相撞.

2 数据前处理

2.1 实测风速修正

多旋翼无人机测风时，机身的抖动与倾斜会影响风速仪的测量精度.李正农等［19］和吴红华等［20］通过推算不同来流风速下无人机的机身倾角开展风洞试验研究，结果表明：机身倾斜使得无人机测风系统测得的风速偏大，需要对无人机实测风速进行修正；而机身倾斜对风向角基本无影响.总结风洞试验结果，得到各风速区间内的修正系数β，见表3，实测风速的修正公式如下：

表3 机身倾角引起的风速修正系数Tab.3 The wind speed correction coefficient caused by inclination of fuselage

式中：为修正的风速时程；uz(t)为无人机系统测得的原始风速时程；β为修正系数.

为说明无人机实测数据的修正过程，以无人机B 在70 m 高度所测得的10 min 风速时程为例，利用表3 中修正系数对实测风速数据进行机身倾斜修正.由图6 可知，修正后的风速时程与原始风速时程的趋势大体相似，但是峰值和部分瞬时值明显减小.表4 给出了修正前后的相关参数，结果表明：机身倾斜导致平均风速、X方向分量均值和Y方向分量均值的相对误差为2.046%和2.463%；机身倾斜导致平均风向的绝对误差仅为0.048°，相对误差仅为0.047%，即机身倾斜对风向的影响很小.

表4 修正前后相关风参数Tab.4 Related wind parameters before and after correction

图6 机身倾斜修正前后风速时程Fig.6 Wind speed time history before and after correction of wind speed caused by fuselage tilt

2.2 实测湍流强度修正

在实测过程中，一方面在风速突变时，无人机为维持飞行稳定性会实时进行必要的姿态调整，使个别瞬时风场数据偏差较大，增大数据时程的波动性；另一方面信号的无线传输有可能产生个别“野点”，导致数据时程的波动性增大.此外，吴红华等［20］对比无人机测量结果与测风塔实测结果，发现无人机测得的湍流强度明显大于测风塔实测值.因此，需要修正无人机测得的湍流强度.五点滑动平均法是一种算法简便且计算量小的平滑预测方法［23］，因此，采用五点滑动平均法修正实测湍流强度.

对已经过机身倾斜修正和分解的X、Y方向风速分量时程进行滑动平均处理，结果如图7 和表4 所示.表4中Uz为z高度的平均风速；φz为z高度的平均风向角；为z高度x方向风速分量平均值；为z高度y方向风速分量平均值；σz，x为z高度x方向风速分量时程标准差；σz，y为z高度y方向风速分量时程标准差；Iz，x为z高度x方向湍流强度；Iz，y为z高度y方向湍流强度.从图7 和表4 可知，经过滑动平均法处理后，风速时程的波动性显著减小；滑动平均法对平均风速、平均风向以及风速分量均值均没有影响，但是能够明显削弱湍流强度，有效消除风速时程样本中部分“野点”的影响.

图7 滑动平均处理前后的风速风向时程Fig.7 Wind speed and direction time history before and after moving average correction of UAV data

采用五点滑动平均法处理后的风向时程可由公式（2）反演得到，如图7（c）所示.

式中：为滑动平均处理后的反演风向时程；、分别为滑动平均处理后X、Y方向风速分量时程.

2.3 风剖面和湍流强度剖面计算理论

2.3.1 风剖面计算理论

不同时刻某一高度处的风速和风向可能不同，因此采用不同时刻多个高度的风速测量值拟合风剖面会产生较大误差.然而，在同一时刻，两台无人机只能分别测量两个高度的风速.为确保风剖面测量的准确性，选取实测期间风向稳定阶段，并利用第1.3 节中双无人机测风的基本原理，将无人机测得的10～120 m 范围内各个高度风速值进行归一化处理，近似得到同一时刻下所有高度的风速值.具体为：设t1时刻无人机A 在100 m 高度处测得风速UA，100-1为参考归一化风速，之后任意高度的实测时刻为t2，则在t2时刻高度z处无人机B 测得的风速为UB，z-2，无人机A 在100 m 参考高度处测得的风速为UA，100-2，从而在高度z处无人机B 测得相对t1时刻的归一化风速UB，z为：

2.3.2 湍流强度剖面计算理论

由于无法同时测量10～120 m范围内所有高度的湍流强度，因此，采用与风剖面类似的归一化处理方法将无人机A、B 在各个高度测得的湍流强度归一化到同一时刻.设t1时刻无人机A 在100 m 高度处测得湍流强度IA，100-1为参考归一化湍流强度，之后任意高度的实测时刻为t2，则在t2时刻高度z处无人机B 测得的湍流强度为IB，z-2，无人机A 在100 m 参考高度处测得的湍流强度为IA，100-2，从而在高度z处无人机B测得相对t1时刻的归一化湍流强度IB，z为：

针对顺风向湍流强度的描述，《建筑结构荷载规范》（GB 50009—2012）［22］给出的计算公式为：

式中：I10为10 m 高度处的名义湍流强度；α为地面粗糙度指数.

3 实测结果分析

3.1 风剖面

指数律是一种描述风剖面的简便方法，《建筑结构荷载规范》（GB 50009—2012）［22］、AIJ-2004［24］、ASCE/SEI 7-16［25］等规范均采用指数律，其通过地面粗糙度指数来评价场地的地表粗糙程度.在气象学中，对数律模型常用来描述大气边界层风剖面，欧洲规范［26］也采用对数律模型，其通过地面粗糙度长度来反映场地的粗糙程度.将无人机测风系统获得的实测数据按第2 节中的相关方法进行修正和归一化处理后，分别采用指数律模型和对数律模型进行海岛典型地貌的风剖面拟合.

以农业温室基地测点1为例，表5列出了无人机实测得到的12 个高度的平均风速、风向、湍流强度以及所有高度在同一时刻的归一化值.结果表明，实测期间平均风向基本在86°～96°，风向比较稳定，接近东风.如图8 所示，将各个高度归一化风速值代入指数律和对数律公式中，分别拟合出测点1 的地面粗糙度指数α=0.170 6，摩擦速度u*=0.414 8 m/s，地面粗糙度长度z0=0.106 9 m.根据同样的方法依次求得海岛地区典型地貌各测点的地面粗糙度指数、摩擦速度和地面粗糙度长度，如表6 所示，所有测点的拟合相关系数R²均大于0.96，拟合效果很好.

表5 测点1实测结果Tab.5 The measured results at point 1

表6 风剖面拟合结果Tab.6 Fitting results of wind profile

图8 测点1实测结果与拟合结果Fig.8 Measured results and fitting results at point 1

3.1.1 指数律拟合

如表6 所示，近海岸所处地形平坦，实测地面粗糙度指数为0.117 9，与《建筑结构荷载规范》（GB 50009—2012）［22］中典型A 类场地的α值和场地特征均吻合，验证了无人机测风系统的可靠性和准确性.农业温室基地各测点的地面粗糙度指数在0.170 6～0.189 6，均值为0.180 6，处于《建筑结构荷载规范》（GB 50009—2012）［22］中B 类和C 类场地之间.农业温室基地中测点5 的α值最大，为0.189 6，测点1 的α值最小，为0.170 6，这与测点所处位置的环境有关.结合图2 分析可知，测点5 位于温室基地中央，地形较为复杂，来流风受地面温室建筑和障碍物干扰较大；测点1 位于温室基地迎风边缘，东侧为田径场，东侧来流经过田径场所受干扰较小.在温室基地实测范围内，东部迎风边缘测点1、2、3的α平均值为0.172 3，中部测点4、5、6 的α平均值为0.185 6，尾部边缘测点7、8、9的α平均值为0.184 0，即温室基地中部的地面粗糙度最大，尾部位置次之，东部迎风边缘粗糙度最小.两处密集建筑群（农科楼和研发楼以东区域）的地面粗糙度指数达到0.27 以上，处于《建筑结构荷载规范》（GB 50009—2012）［22］中C 类和D 类场地之间.结合图3 分析可知，农科楼以东建筑群分布有教学楼、实验楼和办公楼等多、高层建筑，研发楼以东建筑群主要分布有多层学生宿舍楼和教师公寓等建筑，因此，两处密集建筑群来流风剖面的α值比其余两处地貌明显更大.

图9 所示为各类地貌下的实测风剖面和指数律拟合曲线.为便于对比，同时给出采用指数律的我国［22］、日本［24］和美国风荷载规范［25］在相应地貌下的风剖面.结果表明，近海岸实测风剖面与日本规范Ⅰ类地貌在10～60 m高度内偏差较大，与美国规范D类地貌在20～50 m高度内存在一定偏差，整体与我国规范A 类地貌最为吻合.为描述温室基地内风场特性，选取基地中心测点5 作为代表性测点分析.温室基地内实测风剖面处于我国和美国规范的B 类和C 类地貌之间，与美国规范的B 类地貌偏差最大，整体与日本规范的Ⅲ类地貌最为吻合.研发楼以东实测风剖面处于我国规范的C 类和D 类地貌之间，与日本规范的Ⅳ类地貌在20～100 m 高度内吻合较好，整体与美国规范的B 类地貌最为吻合.农科楼以东实测风剖面与我国规范的D类地貌在30～120 m高度内吻合较好，与日本规范的Ⅳ类地貌在20～120 m 高度内吻合较好，整体与美国规范的B类地貌较接近.

图9 风剖面实测结果、拟合结果与规范对比Fig.9 Comparison of measured results，fitting results and standard values of wind profiles

3.1.2 对数律拟合

如表6 所示，摩擦速度u*和地面粗糙度长度z0随场地粗糙度增大而提高，即两处密集建筑群的拟合系数最大，温室基地次之，近海岸最小.近海岸实测摩擦速度u*为0.170 4 m/s，地面粗糙度长度z0为0.006 8 m，处于欧洲规范0 类（z0=0.003 m）和Ⅰ类（z0=0.01 m）场地之间，处于澳大利亚/新西兰规范［27］1类（z0=0.002 m）和2 类（z0=0.02 m）场地之间，符合平坦开阔地貌的粗糙度.农业温室基地各测点u*在0.414 8～0.518 0 m/s 之间，均值为0.450 0 m/s；z0在0.106 9～0.183 9 m 之间，均值为0.153 3 m，处于欧洲规范Ⅱ类（z0=0.05 m）和Ⅲ类场地（z0=0.3 m）之间，处于澳大利亚/新西兰规范［27］2.5 类（z0=0.063 m）到3 类（z0=0.2 m）场地之间.和指数律拟合规律类似，测点5位于温室基地中央，摩擦速度和地面粗糙度长度最大，测点1 位于温室基地迎风边缘，摩擦速度和地面粗糙度长度最小.在温室基地实测范围内，东部迎风边缘测点1、2、3的u*平均值为0.418 8 m/s，z0平均值为0.117 9 m，中部测点4、5、6的u*平均值为0.485 5 m/s，z0平均值为0.170 3 m，尾部边缘测点7、8、9 的u*平均值为0.445 6 m/s，z0平均值为0.171 7 m，即温室基地东部迎风边缘的摩擦速度和地面粗糙度长度较小，中部和尾部较大.两处密集建筑群（农科楼和研发楼以东区域）u*均值为0.645 8 m/s，z0均值为0.867 6 m，接近欧洲规范Ⅳ类（z0=1.0 m）场地，处于澳大利亚/新西兰规范［27］3 类（z0=0.2 m）和4 类（z0=2.0 m）场地之间.

图10 为各类地貌下的实测风剖面和对数律拟合曲线，为便于对比，同时给出指数律拟合曲线和采用对数律模型的欧洲风荷载规范［26］在相应地貌下的风剖面.结果表明，近海岸风剖面的指数律和对数律拟合效果均比较好，实测风剖面与欧洲规范的Ⅰ类地貌也吻合较好.在温室基地内部及研发楼和农科楼以东区域，风剖面的指数律拟合效果比对数律拟合效果要好.温室基地内部实测风剖面与欧洲规范的Ⅲ类地貌吻合较好，研发楼和农科楼以东区域实测风剖面均与欧洲规范存在一定偏差.综合指数律和对数律拟合与分析结果，近海岸风剖面采用指数律和对数律模型均可，实测风剖面与我国规范的A类地貌、欧洲规范的Ⅰ类地貌均吻合较好；温室基地内部风剖面采用指数律模型，实测风剖面与日本规范的Ⅲ类地貌吻合较好；研发楼和农科楼以东区域的风剖面采用指数律模型，实测风剖面与美国规范的B类地貌吻合较好.

图10 实测结果风剖面的对数律拟合结果Fig.10 Logarithmic law fitting results of measured wind profiles

3.2 湍流强度剖面

图11 所示为各类地貌下的实测湍流强度剖面和各国风荷载规范［22，24-26］在相应地貌下的湍流强度剖面.结果表明，实测顺风向湍流强度随高度增加而减小.近海岸实测湍流强度比现有规范值明显要小，其中，与日本规范Ⅰ类地貌偏差最大，与我国规范A类地貌偏差最小.温室基地内实测湍流强度处于我国规范C 类和D 类地貌之间，美国规范B 类和C 类地貌之间，欧洲规范Ⅱ类和Ⅲ类地貌之间，整体与日本规范Ⅲ类地貌吻合较好.研发楼和农科楼以东区域实测湍流强度均处于我国规范C 类和D 类地貌之间；与美国规范B 类地貌相比，实测湍流强度在40 m以下偏大，40 m 以上显著偏小；与欧洲规范Ⅲ类地貌相比，在80 m 以下显著偏大，80 m 以上偏小；实测湍流强度整体与日本规范Ⅳ类地貌偏差较小.

图11 湍流强度剖面拟合结果Fig.11 Fitting results of turbulence intensity profile

由以上可知，现有风荷载规范难以很好地描述各类地貌下湍流强度随高度的变化规律.因此，采用公式（5）对实测湍流强度剖面进一步拟合，拟合参数见表7.结果表明，通过湍流强度拟合的地面粗糙度指数与通过风剖面拟合的结果吻合较好，各测点的地面粗糙度指数随场地差异的变化趋势也与通过风剖面拟合的规律一致，再次验证了无人机测风系统的准确性.拟合得到近海岸的名义湍流强度Iu，10为0.115 3，与《建筑结构荷载规范》（GB 50009—2012）［22］中典型A 类场地值吻合.农业温室基地各测点的Iu，10为0.202 2～0.240 7，均值为0.228 4，接近《建筑结构荷载规范》（GB 50009—2012）［22］中C 类场地值.两处密集建筑群（农科楼和研发楼以东区域）的Iu，10均值为0.339 4，接近《建筑结构荷载规范》中D 类场地的规定值.