0 引　言

风是使单体树木和森林系统受灾的主要非人为因素^[1-4]。巨大的强风能量对树木造成危害，不仅影响生态系统的稳定性^[4]，还会给林业生产带来严重的经济损失^[5-8]。保障树木风致安全的首要问题是掌握强风下树木-风相互作用特性。特别是，树木的风致响应。

针对树木-风相互作用特性这一科学问题，通常采用的研究手段主要包括理论分析^[9-10]、数值模拟^{[3, 11-12]}和现场实测。陶嗣巍利用计算机辅助设计软件Pro/E（Pro/Engineer）中的簇表与UDF（User Defined Function）功能，基于参数化建模思路和分形理论，提出了一种快速建立树木结构几何模型的方法并分析了白毛杨的风振动特性^[13]。宋晓鹤基于非线性假设，将东北鱼鳞云杉模拟成弹性变截面圆锥杆，建立了东北云杉树风倒非线性力学模型，并对其进行了静力学分析^[14]。王琳将东北云杉树干近似为横截面随着高度成幂函数变化的弹性杆，一端固定，自由端连有质量团，建立了力学模型，分析了云杉的受迫振动响应^[15]。崔云静通过Pro/E三维绘图软件建立了新疆杨的实体模型，并运用有限元软件Ansys workbench对模型分别进行静力和动力分析，得到新疆杨在模拟风载荷作用下的应力、应变和位移变化情况^[16]。相比较其他两种研究手段，现场实测具有真实反映树木及风场的实际状况、数据可靠的优点，是开展相关研究较好的方式。在强风天气下获得树木风致振动实测结果具有一定难度，到目前为止，大多数的研究主要是基于低风速条件下的现场实测^[17-23]，且针叶树的风致振动行为是目前现场实测研究关注的主要对象^{[17-20, 23-25]}。例如Dupont等人基于现场实测结果发现了边缘流中湍流风场频率峰值与海岸松（Pinus pinaster Ait）基本振动模式相关性较大^[23]。

落叶阔叶树相对针叶树而言结构要复杂的多，关于其在强风条件下的风致运动实测研究目前较少。但事实上，落叶阔叶树与针叶树一样都具有很高的风敏感性，其覆盖面积在各树种中也占有较高的比重。在欧洲的许多森林中，针叶树正逐渐被落叶阔叶树所取代^[5]；且大多数欧洲城市的行道树均为落叶阔叶树种^[26]。因此，开展强风天气下落叶阔叶树的振动响应分析十分必要。

综上，开展强风天气下落叶阔叶林的风场特性及其风致运动特性研究具有重要意义。本文以黑龙江帽儿山森林生态系统国家野外科学观测研究站的落叶阔叶林为研究对象，实测研究了三次强风时段的林间风场特性，并研究了某一白桦树的风致动力行为，采用谱分析方法分析了脉动风特性、树木脉动风荷载及风致响应。本文包括引言在内共分为4个部分：第1节描述了样本树、现场实测的具体位置及监测系统，还简要介绍了功率谱分析的方法；第2节介绍并讨论了湍流风场特性、白桦树的风致响应和风树相互作用的结果；第3节总结了我们的研究。

1 试验概况 1.1 现场实测地点和样本树

观测地点选为黑龙江帽儿山森林生态系统国家野外科学观测研究站老爷岭主站点（东经127°40′，北纬45°24′）附近。观测点位于东南-西北的山谷中，周围的平均坡度约为9°^[27]。其周围的主要树种包括榆树、水曲柳、白桦、山杨和胡桃楸等^[27]。观测对象选为该站点的一棵白桦树（Betula platyphylla Suk），如图1所示。其中，高50 m的通量塔位于样本树西南方向10 m处，塔与样本树之间无其他树木遮挡。

在2016年到2018年期间，老爷岭主站点树木冠层的平均叶面积指数（mean leaf area index，LAI）和森林平均林分高度分别约为6.95±0.65 和18～20 m^[27-28]，平均断面积为24.16 m²/hm^2[29]。所选的白桦树观测对象高23 m，胸径（即，地面以上1.3 m处树干的直径）0.28 m，树冠中心高度16.8 m。

1.2 现场观测系统

观测系统包括：风速仪、加速度传感器以及数据采集系统。风速仪为英国Gill公司生产的Gill Wind Master三维超声波风速仪，风速测量范围为0～45 m/s，分辨率为0.01 m/s，风速小于12 m/s时精度小于1.5%；风向测量范围为0°～359°，分辨率为0.1°，风速小于12 m/s时精度为2°。本次试验采样频率设置为10 Hz，长期离线采集。Gill Wind Master三维超声波风速仪如图2所示。在通量塔塔身的5、10、20、35、50 m处各安装一个Gill Wind Master 三维超声波风速仪，通量塔如图3所示。

采用扬州晶明科技有限公司生产的三向加速度传感器（型号：AIT2500）及JM5981A数据采集系统测量并记录树木的摇摆，如图4和图5所示。加速度传感器量程为±2.5g，最大冲击极限±50g，横向灵敏度小于5%，外形尺寸45 mm×45 mm×50 mm，重量200 g。本试验在执行中，加速度传感器布设在白桦树树干上，距地面高度为7.97 m。试验时，采集频率设置为50 Hz，长期离线采集。试验现场整体测试仪器布置如图6所示。

1.3 三次强风事件

选取帽儿山实测地点2019年11月2日、11月5日和11月8日这三个强风天气中所测的风特性和树木响应观测数据进行后续的研究与分析。在不到十天的时间里实测地发生了三次强风，每个强风天气的一般气象条件见表1。三次强风中11月5日的强风最强，其次为11月8日的。塔顶风速仪（z/h = 2.17，z为测点高度，h为树高）记录的11月5日强风最大平均风速为11.42 m/s，最大阵风为19.82 m/s。

1.4 风速及风向的处理计算方法

本试验所用的Gill超声波风速仪可以获得笛卡尔坐标系下三个方向的风速矢量数据，U_x、U_y、U_z。其中，U_x正方向定义为与方位角系统中正北方向一致的方向，U_y正方向定义为从正北方向逆时针旋转90°的方向，U_z正方向为平行于风速仪安装轴垂直向上。

依照中国规范，在计算分析中设置基本时距为10 min。则在基本时距内的三维平均风速为：

$ {\overline {{U}}} =\sqrt {{{{\overline {{U}}_x }}^2} + {{{{\overline {{U}}_y }} }^2}{\text{ + }}{{ {{\overline {{U}}_z }}}^2}} $

(1)

式中 ${{\overline {{U}}_x }}$ 、 ${{\overline {{U}}_x }}$ 和 ${{\overline {{U}}_x }}$ 分别为10 min基本时距内的三维风速平均值，通过下式计算得到：

$ {{\overline {{U}}_x }} =\frac{1}{N}\sum\limits_{i=1}^N {{U_x}(i)} $

(2)

$ {{\overline {{U}}_y }} =\frac{1}{N}\sum\limits_{i=1}^N {{U_y}(i)} $

(3)

$ {{\overline {{U}}_z }} =\frac{1}{N}\sum\limits_{i=1}^N {{U_z}(i)} $

(4)

式中N为基本计算时距内的风速样本数。

风向角定义为风向与正北方的夹角，沿顺时针增加。则瞬时风向角及10 min基本计算时距内的风向角，即瞬时风向及平均风向可用下式表示：

$ \beta '{\text=}\arctan \frac{{{U_x}}}{{{U_y}}} $

(5)

$ \beta {\text=}\arctan \frac{{\overline {{U}}_y }}{{\overline {{U}}_x }} $

(6)

基于前述10 min基本计算时距内的风向角，可计算得到顺风向风速、横风向风速及竖向风速分别为：

$ U={U_x}\cos \beta + {U_y}\sin \beta $

(7)

$ V={U_y}\cos \beta + {U_x}\sin \beta $

(8)

$ W={U_z} $

(9)

则顺风向脉动风速、横风向脉动风速及竖向脉动风速可分别表示为：

$ u={U_x}\cos \beta + {U_y}\sin \beta - {{\overline {U}}} $

(10)

$ v={U_y}\cos \beta + {U_x}\sin \beta $

(11)

$ w={U_z} - {{\overline {{U}}_x }} $

(12)

1.5 谱分析方法

脉动风荷载作用下的树木风致响应分析方法包括时程法和功率谱密度法。时程法通常需要知道树木气动力参数，如：树木阻力系数C_D和树冠垂直于气流的迎风面积A等，这些参数的获得较为困难。功率谱密度法可直接进行树木结构实测响应分析，是常采用的分析方法^{[17-18, 20, 24, 30-33]}。功率谱密度法的一般过程如图7所示^{[24, 31]}。

对于一阶线性系统，图7所示的传递关系可通过方程式^[24]表达为：

$ {P_l}\left( f \right)={\left( {\rho {C_D}AU{T_a}\left( f \right)} \right)^2}{P_u}\left( f \right) $

(13)

式中： ${P_l}\left( f \right)$ 为风荷载的功率谱密度； $\rho$ 为空气密度，kg/m³； $U$ 为顺风向平均风速，m/s； ${T_a}\left( f \right)$ 为气动导纳函数； ${P_u}\left( f \right)$ 为风速的功率谱密度。

$ {P_t}\left( f \right)={\left| {{T_m}\left( {{\rm{i}}f} \right)} \right|^2}{P_l}\left( f \right) $

(14)

式中：P_t（f）为树响应（位移或者加速度）的功率谱密度； ${\left| {{T_m}\left( {{\rm{i}}f} \right)} \right|^2}$ 为机械传递函数，其中i为虚数单位。

在实际中由于C_D和A并不容易获得^[34]，且风荷载的作用也并不仅仅是一维的，所以风荷载的功率谱并不能由式（13）得到，在研究树木在风载下的响应的问题通常采用雷诺应力 $ \tau $ 作为风荷载的度量。雷诺应力是湍流脉动速度引起的动量传递，由法向和切向压力构成^{[18, 24]}。平均雷诺应力 $\bar{\tau }$ 可表达为：

$ \bar \tau = - \rho \overline {uw} $

(15)

式中u、w为顺风向和竖向的脉动风速，m/s；同理，下文中的v为横风向脉动风速，m/s。

由图7可知，顺流向脉动风速功率谱和树响应功率谱之间的关系可通过机械传递函数 ${\left| {{T_m}\left( {{\rm{i}}f} \right)} \right|^2}$ 表示， ${\left| {{T_m}\left( {{\rm{i}}f} \right)} \right|^2}$ 提供了树木对不同频率湍流的响应的信息，同时 ${\left| {{T_m}\left( {{\rm{i}}f} \right)} \right|^2}$ 也是对湍流风场到树木的能量传递效率的一种测量^{[17-18, 20]}。通过式（15）的平均雷诺应力表达式以及功率谱估计可得到雷诺应力功率谱 $ {P_{\bar \tau }}\left( f \right) $ 。将 $ {P_{\bar \tau }}\left( f \right) $ 代替 ${P_l}\left( f \right)$ 代入式（14），则机械传递函数 ${\left| {{T_m}\left( {{\rm{i}}f} \right)} \right|^2}$ 可通过下式得到：

$ {\left| {{T_m}\left( {{\rm{i}}f} \right)} \right|^2}=\frac{{{P_t}\left( f \right)}}{{{P_{\bar \tau }}\left( f \right)}} $

(16)

以上所有功率谱估计均应用MATLAB软件（8.3 R2014a，Mathworks，纳提克马，马萨诸塞州，美国）中的Pwelch函数实现。通过对近5小时数据历程的约30个基于10 min计算时距的功率谱密度进行平均，得到了z = 0.35h高度处的树干加速度、z = 2.17h、1.52h、0.87h、0.43h和0.22h高度处的风速以及雷诺应力的时间平均功率谱密度结果。其中，树响应的功率谱分别计算了水平纵向（东北-西南方向，a_l）和水平横向（东南-西北方向，a₂）加速度的功率谱密度。其中，用于各谱时间平均计算的具体10 min计算时距个数见表1。采用双对数坐标绘出各个功率谱密度图，其中纵坐标功率谱密度乘以频率再除以方差进行归一化，即 ${S_i}=\dfrac{{n{S_{ii}}(n)}}{{\sigma _i^2}}$ (i = u, v, w, a)，横坐标采用莫宁坐标，即 $f=\dfrac{{nz}}{U}$ 。

2 结果与分析

选取3次强风过程中z = 2.17h高度处风速传感器所测得的平均风速大于5 m/s的时间段进行后续的所有分析。每个强风天气内所选取的时间段约为5 h。

2.1 平均风速与风向

如图8所示，为在三次强风天气里测量并处理得到的顺风向平均风速廓线，其中高度（z）采用z = 0.87h（近冠顶）处的高度作为参考高度（z_r）进行归一化，风速（u）以对应的z = 0.87h处的风速作为参考风速（u_r）进行归一化。图中虚线代表样本树冠顶高度。可以发现在三次强风天气内有着几乎相同的顺风向平均风速廓线。在冠层顶部以上平均风速沿测点高度近似呈现对数分布。

绘制了3个强风天气的平均风速和风向的时程曲线，如图9所示，图中阴影区域为用于本文分析的强风时段。

从图9可以看出，三个强风展示了相似的风特性。从z = 1.25h到0.87h的高度，即到达冠层时，平均风速急剧下降（图9（a）），这可能是由于冠层树叶和树枝等从气流中提取动量导致的。为了说明该结论，我们选取2021年5月7日研究区域附近某一空旷区域（无树木影响区域）的风场数据与本文数据进行对比分析。2021年5月7日研究区域附近某一空旷区域的20 m、35 m和50 m的平均风速时程曲线如图10所示。

通过对比图10和图9（a）可以发现，没有树木影响的区域下，随着测点高度的降低平均风速稍有减小，幅值很小；从z = 1.25h到0.87h，风速并没有急剧减小的现象。而在有树木的影响下，从z = 1.25h到0.87h，即到达冠层时，平均风速急剧下降（图9（a）），可知这是由于冠层树叶和树枝等从气流中提取动量导致的。

在冠层下面的两个高度（即z = 0.43h和0.22h）的风速几乎相同。通过表1和图9可以发现对于同一测点高度，在风速较大的强风期段，平均风向变动范围较小，基本保持稳定；而在z = 2.17h～1.52h高度处平均风速近似小于2.0 m/s（z = 0.87h高度处平均风速近似小于1.0 m/s，z = 0.43h～0.22h高度处平均风速近似小于0.5 m/s）的低风速时段内，风向变化范围很大。还可以发现平均风速方向随高度方向发生偏转，冠层高度范围以内，偏转梯度大，冠层高度以外，偏转梯度较小；最大风偏转角约为67°，从z = 2.17h处的西北风转变为z = 0.22h处的西南风，这是由于树木的扰动作用对风特性产生了影响。三次强风时段的平均风向廓线如图11所示。

2.2 脉动风速功率谱的分析

脉动风速功率谱反映脉动风能量随频率变化的分布情况，是非常重要的脉动风特性之一。当脉动风频率接近或等于高耸结构的固有频率时，湍流脉动分量的大小将会显著影响结构的响应。脉动风速功率谱包括水平纵向（顺风向）脉动风速功率谱、水平横向（横风向）脉动风速功率谱及竖向脉动风速功率谱。在这项研究中不同测点高度处顺风向、横风向和竖向的归一化时间平均脉动风速功率谱如图12所示。

三次强风事件中的三个风速分量（u、v和w）在z = 1.52h和2.17h高度处的时间平均功率谱展现出与常见的大气表层谱相似的形状，即功率谱含能区的斜率近似以 $ {f^{ + 1}} $ 的函数上升，在惯性子区近似以 $ {f^{ - 2/3}} $ 的函数下降（图12）。三次强风事件中，相比其他测点的u-谱，z = 0.22h测点的u-谱在耗散区的斜率最陡。这是因为冠层以下植被背后的尾迹结构发展，导致湍流耗散加速^[23]。

由图12（a）可知，顺风向脉动风速功率谱（S_u）的归一化峰值频率（f_pu，即最大功率对应的频率）随风暴强度的增加而减小；对于较强的风暴（11月5日和11月8日）来说，顺风向脉动风速功率谱（S_u）的最大功率随测点高度的减小而增加，结合图12（a）的顺风向湍流强度廓线图可以发现这是由于树木的扰动作用增加了冠层内风场的湍流强度造成的。李正农等人通过风洞试验也提出树木的背风面存在风影响区，此区域的流场特性会发生显著改变，体现为风速降低和湍流强度增大^[35-36]。由图12（b）可知，横风向脉动风速功率谱（S_v）的归一化峰值频率（f_pv）随风暴强度的变化规律与S_u的相同；S_v的最大功率，对于较强的风暴（11月5日和11月8日），除了测点z = 0.22h，随测点高度的减小而增加；对于较弱的风暴（11月2日），各测点的最大功率较为接近，无明显规律。各个风暴事件的竖向脉动风速功率谱（S_w）（图12（c））的最大功率在各测点之间较为接近，无明显规律。

总的来说，观察三个强风事件的u-谱和v-谱发现，相邻两个测点之间的最大功率的差值随风暴强度的增加而增加。同时，对于同一强风事件的同一速度分量在不同高度处的功率谱归一化峰值频率几乎相同。11月2日、5日和8日强风的u、v和w谱的归一化峰值频率分别为f = nz/u = 0.33、0.35±0.01和0.60±0.09、f = nz/u = 0.28±0.02、0.28±0.02和0.67±0.10和f = nz/u = 0.30±0.02、0.29和0.57±0.01。基于脉动风速谱及归一化频率公式整理可得各测点高度的脉动风峰值频率，发现三次强风期间各测点的顺风向和横风向脉动风速（u和v）的峰值频率约为0.03 Hz，竖向脉动风速（w）的峰值频率在0.05～0.07 Hz之间，不同测点高度之间有细微差别。

对同一强风事件同一测点高度下三个风速分量（u、v和w）的功率谱密度图进行了绘制，如图13所示。

对于11月2日的强风事件，5个测点高度的功率谱具有相似的变化规律，由u → v → w，峰值频率越来越大，最大功率越来越小。对于11月5日和8日的强风时段，所有测点高度的峰值频率变化规律具有相似变化规律，即f_pu≈ f_pv<f_pw。11月5日的最大功率，除z = 0.87h测点外，由u → v → w，越来越小。11月8日的最大功率，对于测点z = 0.22h和0.43h，由u → v → w，越来越小；对于测点z = 0.87h、1.52h和2.17h，v-谱最大，u-谱次之。三个风速分量的不同的功率谱密度表明了林分树冠上部、顶部及冠层内湍流的各向异性特征。总体看来，随风暴强度的增加，湍流各向异性增加；对于同一风暴，湍流各向异性随测点高度的增加而减小。

2.3 树响应及风荷载的分析

用三次强风期间树干加速度分量（东北-西南和东南-西北方向）来表征风荷载作用下样本树的运动，如图14所示为三次强风期间样本树加速度响应的时程曲线。通过图14可以发现，在三次强风期间，白桦树均呈现出不规则的复杂的摇摆运动。比较三个强风天气里样本树的运动大小可以发现，随着风强度的增加（“11月2日”<“11月8日”<“11月5日”），白桦树在两个方向的运动变大。且在强风的影响下，白桦在东北-西南和东南-西北两个方向的运动大小相近，可知样本树在强风期间的运动轨迹近似为圆形。

z = 2.17h、1.52h、0.87h、0.43h和0.22h的5个测点处的雷诺应力功率谱、基于z = 0.87h测点处的雷诺应力功率谱和树木加速度响应功率谱得到的机械传递函数以及z = 0.35h高度处的树木水平纵向和横向加速度响应功率谱分别如图15～图17所示。

通过图15可以发现，总体而言，在三次强风期间各测点低频区雷诺应力功率谱密度增加，到达功率谱密度的最大值后，又逐渐减小，但随着风暴强度的增加，减小的斜率越来越平缓。但在11月5日的强风暴期间，在z = 0.87h测点处（近冠层顶部），当功率谱密度值达到最大后，在一定频率范围内，呈现出一段功率谱密度几乎在相同水平的频率区间。11月2日、5日和8日的τ-谱的归一化峰值频率分别为f = nz/u = 0.58 ± 0.08、0.84 ± 0.30和0.55 ± 0.09。

如图16所示，在11月2日、5日和8日的强风期间a₁和a₂谱的归一化峰值振动频率分别为f = nz/u = 3.44和3.28、2.67和2.55以及3.34和3.19。同一强风期间，白桦树两个加速度分量的功率谱差异很小，这与Peltola等对黑松进行测量后得到的结果相同^[19]，也与Mayer得到的结果一致^{[24, 31]}。Mayer曾指出这种较小的差异性可能是由树木树冠结构部分的不对称性分布、平均风向以及根茎的发展导致的^[37]。基于加速度响应谱及归一化频率公式整理可得到白桦树在三次强风中a₁和a₂谱的峰值振动频率分别为0.30和0.28、0.27和0.26以及0.31和0.29 Hz，不同强风天气之间有细微差别。

在三个强风的加速度响应谱中还发现了两加速度分量基本一致的高于峰值振动频率的其他两个显著的峰值，这可能对应于白桦树的二阶和三阶谐波频率，分别约为0.76和1.72 Hz、0.69和1.52 Hz以及0.65和1.65 Hz。Peltola等人提出相比于间伐过的林分中的黑松，未间伐的林分中的黑松因具有相对大的灵活性的小直径的树木元素而使入射能量被划分为更高的共振模式^[19]。Scannell也提出这种更高阶的振动模式可能是由于枝叶等树木元素造成^[38]。

为了得到树木自由振动频率，还进行了树木拉力试验，并对相应的加速度响应数据进行傅里叶变换分析后，得到白桦树自由振动下两个加速度分量的峰值频率，分别为0.28 Hz和0.26 Hz，与强风下白桦树的峰值振动频率一致性较好。Mayer、Amtmann等对针叶树进行研究也得到了相同的结果^{[17-18, 31-32, 39]}。

如图17所示，三次强风下白桦树的机械传递函数均在白桦树加速度响应功率谱的峰值频率处达到峰值，然后在更高频率处下降。这与Gardiner、Stacey等对针叶树在风下进行测试得到的结果相同；Gardiner、Stacey等表明，通过将针叶树视为具有单个集总质量的简单阻尼谐振子并使用冠层质量、茎质量、茎弹性和阻尼比的测量值对风致树木响应进行模拟后，得到的模拟机械传递函数与通过试验得到的树响应的机械传递函数非常吻合，也就是说，对针叶树在风下进行测试得到风致树木响应的机械传递函数几乎与受迫阻尼谐振子的传递函数相同^{[17-18, 30, 32, 39]}。总的来说，森林可以被看作是一组有效地吸收树木共振频率处的湍流谱能量的谐波振荡器^[17]。此外，和加速度响应谱一样，两个加速度分量的机械传递函数之间只有微小的差异。

3 结　论

以黑龙江帽儿山森林生态系统国家野外科学观测研究站的落叶阔叶林为研究对象，实测研究了三次强风时段的林间平均风与脉动风特性，分析了林中某一白桦树的脉动风荷载及风致响应。主要得到以下几点结论：

1）树木的扰动作用导致平均风随高度发生偏转，在冠层高度范围以内偏转梯度大，冠层高度以外偏转梯度较小；最大风偏转角度近似为67°。较大的风偏转角可能对林木或者林中结构产生扭曲作用，在开展林木抗风分析、防风减灾能力评估以及林中柔性结构物抗风设计等方面应引起足够的重视。

2）树木的扰动作用增加了冠层内风场的湍流强度，这一影响随着来流风强度的增加而增加。三次强风作用下，树木在脉动风作用下呈现出复杂的随机振动，振动响应的峰值频率与树木自由振动的峰值频率一致，与来流风的脉动风特性无关。

3）阔叶树的机械传递函数的峰值频率与响应功率谱的峰值频率基本一致，树木风致响应的机械传递函数几乎与受迫阻尼谐振子的传递函数相同，这一点与针叶树具有相似的特征。