引言

激光在大气中的传输特性受到各种大气参数的影响，其中光路上的风速和风向对热晕效应的分析具有重要的意义。高能激光通过大气传输时，一部分能量被空气中的分子和气溶胶吸收，加热空气形成负透镜，使光束发散、弯曲、变形，这种效应称作光的热晕效应。热晕效应计算中，一般假定风速为一定值，但实际大气中的风场比较复杂，其对热晕的影响将随外界条件的改变而发生很大变化，若要分析实时实地的激光传输效果，则相应风速的分析是必不可少的。不仅如此，边界层中的水汽、热量、动量和污染物等的输送也受到风速、风向的支配。因此，研究风速廓线特点具有十分重要的价值。

众多研究人员对近地层风速廓线做了多方面的研究。刘学军等^[1]分析了天津城市建筑群对近地层风速廓线的影响，发现不同性质下垫面上近地层的平均风速廓线存在显著差异。赵鸣等^[2]利用天津250 m塔层资料得到风切变的若干特征。申华羽等^[3]应用江苏和吉林的测风塔资料，探讨了一种推算近地层内任意高度上风能参数的新方法。业已证明，在中性大气层结条件下，近地面层风速随高度基本上呈现对数变化规律^[4-6]，而地面粗糙度是对数风速廓线中一个重要的参量，它反映了下垫面的形态结构对风速的影响。因此，粗糙度的确定对获得风速廓线特点至关重要。在小尺度上，粗糙度的一个度量是各个障碍物即粗糙元的空间密度，Lettau^[7]提出了以根据粗糙元的平均垂直范围、平均轮廓及地域大小来计算粗糙度长度的关系式，适用于粗糙元间距均匀、不太靠近、高度和形状相似的情况。Kondo等^[8]提出了类似的表达式，并考虑了各个粗糙元的变化，适用于城市建筑。Chamberlain^[9]把应用于海面粗糙度长度计算的查诺克关系式扩展到用于吹沙和吹雪条件下粗糙度的计算。近年来许多研究人员对不同下垫面的粗糙度进行了研究分析。张强等^[10]讨论和分析了控制粗糙度的物理因子，给出确定城市粗糙度长度的简化公式。韩素芹等^[11]利用天津255 m气象塔近10年观测资料，定量计算空气动力学粗糙度随下垫面的改变情况，研究了城市化对边界层风场结构的影响。何清等^[12]利用塔克拉玛干沙漠北部肖塘地区梯度自动站观测数据，计算出该区域的粗糙度，发现其随下垫面性质变化明显，理论结果对生态环境研究及治理具有重大意义。

本文根据35 m铁塔测量系统测得的2007—2009年4月的风速及其他一些大气参数，采用空气动力粗糙度长度z₀(定义为近地层风速向下递减到零时的高度)，研究了位于合肥郊区某非均一地形的地面粗糙度及在春季中性大气层结条件下的近地层风速廓线，这可为分析特定大气条件下的激光大气传输效果提供一些参考。

1 风速廓线理论

风速测量设备及周边地形环境如图 1所示。北面是四层楼建筑物，东面靠近大片树林，南面为水库，西面是草坪、少数树木和比较低矮的建筑物。因此不同方向的风受阻挡的程度不同，这为粗糙度的确定带来一定困难。35 m铁塔在五层上(离地面高度从低到高分别为5 m、10 m、15 m、23 m、31 m)装有风速风向、温度、湿度等多种传感器，可以测量5个高度的风速风向及温度等，第二层每10分钟存储一组数据，其余四层每20秒存储一组数据。

在近地面层(边界层下部10%左右，高度在50~100 m^[5])风速随高度是否遵从对数规律分布主要由当时的大气层结决定。大气层结可分为不稳定、稳定和中性三种，在静力中性条件下，风速随高度呈现对数变化规律。因此首先需要判断当时的大气稳定性，再利用有关关系式进行粗糙度的分析，进一步推出风速廓线。

本文利用梯度Richardson数(Ri)法^[12-14]来判断大气稳定度。有关文献给出了Richardson数大气稳定度分类标准^[12-13]，若Ri满足-0.0807＜Ri＜0.0918，根据其标准认为大气处于静力中性状态。

在描述近地面层结构时最常用的就是Monin-Obukhov相似性理论^[5]，由此理论可得近地面层的无量纲风速：

$ \frac{{kz}}{{{u_{\rm{*}}}}}\frac{{\partial \bar u}}{{\partial z}} = {\varphi _m}\left( {\frac{z}{L}} \right) $

(1)

式中z为距离地面的高度；u为z高度上的平均水平风速；k称为Karman常数，Tennekes和Lumley指出该常数取值范围为0.33~0.4^[6]，通常取值0.4；u_*是摩擦速度，表征下垫面摩擦应力的大小；L为Monin-Obukhov长度。利用Businger-Dyer关系式^[4-5]对式(1) 进行高度积分可以得到风速廓线^{[4, 15]}：

$ \bar u = \frac{{{u_*}}}{k}\left[{\ln \left( {\frac{z}{{{z_0}}}} \right) + {\psi _m}\left( {\frac{z}{L}} \right)} \right] $

(2)

在不同层结条件下，稳定函数${{\psi _m}\left( {\frac{z}{L}} \right)}$有具体的表达关系。在静力中性气流(z/L=0) 极限情况下，Monin-Obukhov长度L趋于无限长，${{\psi _m}\left( {\frac{z}{L}} \right)}$等于零，上式便简化为常见的对数风速廓线：

$ \bar u = \frac{{{u_*}}}{k}\ln \frac{z}{{{z_0}}} $

(3)

式中z₀即空气动力学粗糙度长度。一般在边界层中很少出现完全中性，但是在强风和阴天时边界层近似中性。

利用两个高度的风速差根据式(3) 可以推出：

$ \ln {z_0} = \frac{{{u_2}\ln {z_1}-{u_1}\ln {z_2}}}{{{u_2}-{u_1}}} $

(4)

利用上式由两层高度和风速可以求出粗糙度z₀。

剔除恶劣天气下测量结果，共得到2007—2009年4月77天的数据，根据Ri法选出满足中性大气层结条件的风速。本文同时选取第一层和第三层风速、温度等及第一层和第五层风速、温度等分别计算对应的Ri数，选出两个Ri数都满足-0.0807＜Ri＜0.0918的风速，按照小时进行平均，同时将第二层风速也按照小时进行平均。为排除外界干扰导致的某层风速测量结果偏离其他塔层过大的情况，将第i层和第i+1层(i=1, 2, 3, 4) 的风速和高度分别代入式(4)，求得4个粗糙度值，进行两两比较，若任一差值超过4 m，则这一组是无效数据。这样便得到了逐时平均符合中性大气层结条件的五层有效风速值。

由式(3) 可以得到每层风速：

$ {\bar u_i} = \frac{{{u_*}}}{k}\ln \frac{{{z_i}}}{{{z_0}}}, \;\;\;\;i = 1, 2, \cdots, 5$

(5)

式中u_i即第i层的风速，z_i为第i层离地面的高度。将第一层风速作为参照风速，可以由上式推出：

$ {\overline u _i} = {k_i}{\overline u _1}, \;\;\;\;{k_i} = \frac{{\ln \left( {{z_i}/{z_0}} \right)}}{{\ln \left( {{z_1}/{z_0}} \right)}} $

(6)

式中k_i为第i层风速u_i和第一层风速u₁之间关系的比例因子，可以看出若粗糙度z₀和第一层高度z₁确定，则比例因子仅是第i层高度的函数，这样点绘第i层风速与第一层风速的散点图应当是线性相关的^[16-17]。

由于随着层数的增加，第i层与第一层风速的散点图愈加分散，可以点绘第i层与第i-1层风速之间的散点图^[16]。若将第i-1层风速作为参照风速，第i层的风速可以表示为：

$ {\overline u _i} = {c_i}{\overline u _{i-1}}, \;\;\;\;{c_i} = \frac{{\ln \left( {{z_i}/{z_0}} \right)}}{{\ln \left( {{z_{i-1}}/{z_0}} \right)}} $

(7)

式中c_i为第i层风速u_i和第i-1层风速u_i-1的比例因子。由式(6) 和式(7) 可以得到k_i和c_i之间的关系为：

$ {k_i} = \prod\limits_{i = 1}^i {{c_i}} $

(8)

则k₁=c₁=1.0，k₂=c₂。

2 统计结果与分析 2.1 粗糙度的拟合分析

粗糙度与下垫面的形态结构有关，和粗糙元存在着一一对应关系，主要由粗糙元的形状大小、空间密度等几何特性决定。前面提到的Lettau、Kondo和Yamazawa提出的计算粗糙度的方法主要应用在粗糙元比较均匀的情况，而35 m铁塔所处的地形比较复杂，可以利用35 m铁塔测量的风速及其他大气参数根据上述理论分析其粗糙度。

图 2点绘了第i层与第i-1层风速之间的散点图，进行线性回归拟合，可以得到比例因子, 再由式(8) 进一步推出第i层风速u_i和第一层风速u₁之间关系的比例因子k_i。

粗糙度长度可以由关系式${k_i} = \frac{{\ln \left( {{z_i}/{z_0}} \right)}}{{\ln \left( {{z_1}/{z_0}} \right)}}$通过拟合(k_i，lnz_i)得到。图 3绘出了二者之间的散点图，进行线性回归拟合，可以看出相关性很大，相关系数R=0.998。计算出粗糙度长度为1.25 m，与粗糙度长度的标准值基本吻合^[1-2]。

2.2 粗糙度与风速关系

从空气动力学角度来讲，可以认为对数风速廓线式(3) 中的空气动力学粗糙度长度z₀是气/固界面上“无滑移”的厚度^[18]，那么外界大气扰动尤其是风速的改变会导致这个厚度随之出现起伏，所以粗糙度不仅取决于地表粗糙元的大小、形状、空间密度等几何特性，还与流经地表的流体的性质和运动状态有关。

由前面逐时平均统计的符合中性大气层结条件的五层有效风速值，利用式(3) 根据最小二乘法原理进行拟合^[19]，得到每组数据所对应的粗糙度和摩擦速度。粗糙度和第一层风速呈较好的负相关(如图 4)，随着风速的增大，空气动力学粗糙度逐渐减小。因为粗糙度与流体性质有关，所以随着外界条件的改变，不同时刻粗糙度出现起伏，变化范围为0.038~4.903 m。在利用上述方法拟合时，先求出摩擦速度，再代入一定的关系式求解粗糙度，而二者之间存在对数关系，所以即使摩擦速度有很小的变化，也会导致粗糙度出现很大的起伏，特别是当风向改变时，风速有可能会变得很小，甚至趋向于零，这时的摩擦速度也会很小，用这样的数据再去求粗糙度，因为对数关系会导致粗糙度很大，结果便出现偏差。

对3个月同一时刻的数据进行平均，得到24组风速值，同样利用式(3) 根据最小二乘法原理进行拟合，得到每组数据所对应的粗糙度和摩擦速度。图 5绘出了第一层风速与粗糙度随时间变化趋势图线(图 5a)以及二者之间的相关性(图 5b)。由图 5a可以看出，风速在夜间较小，白天逐渐增大，到中午时刻左右达到最大值，然后又逐渐减小，粗糙度和风速变化基本上呈相反变化趋势。其原因在于风速从较小的量逐渐增大时，湍流扰动加强，近地面水平动量的垂直通量增加，表征其大小的量摩擦速度（$u_*^2 = \left| {\overline {u'{{w'}_s}} } \right|$）增加，即近地面摩擦应力增大，粗糙度相应有所增加，但同时由于湍流强度加强，导致气/固界面上有更多的动量下传到“无滑移”层，引起风速为零的高度减小，即空气动力学粗糙度减小，综合效果便表现为粗糙度随风速的增大而减小。由图 5b可以得出相关系数-0.953，负相关性非常好。可见，空气动力学粗糙度长度随风速基本上呈现相反方向变化趋势。

2.3 对数风速廓线

由求出的摩擦速度和粗糙度，利用式(3) 每隔5 m计算出相对应的风速，拟合出50 m高度内的对数风速廓线，标准偏差SD最小值为0.020，最大值为0.116，大部分都集中在0.04左右；相关系数R最小值为0.9855，最大值为0.9970，相关性非常好。图 6绘出了4个时刻(02:00，10:00，15:00，23:00) 的风速廓线，可以看出，利用梯度Ri数按照中性大气层结条件进行筛选统计得到的风速廓线基本上满足对数分布规律。

3 结论

本文通过梯度Ri数筛选出合肥郊区在2007—2009年4月满足中性层结大气条件的五层有效风速，计算分析非均一地形的粗糙度，结果(1.25 m)与粗糙度长度的标准值基本吻合。由于所研究地形比较复杂，其粗糙度不仅与下垫面的形态结构有关，外界扰动也会使粗糙度出现起伏。将3个月同一时刻的风速进行平均，统计出24组数据(对应24个时刻)，分别用最小二乘法按照对数风速廓线进行拟合，分析当地不同时刻的空气动力学粗糙度长度，结果表明它不仅与地形这样固定的物理特征有关，还与流体性质有很大关系，与风速有很好的负相关性；利用统计的每组风速分别拟合50 m高度内的对数风速廓线，相关系数都很大，证明了利用梯度Ri数法筛选出的风速基本上满足对数分布规律。