Research on airport boundary layer height based on Doppler LiDAR

国务院和中央军事委员会发布的相关文件指出，各种低空空域的垂直范围为距地面高度1000 m以下。对飞行活动而言，低空复杂多变的飞行环境是影响飞行品质和飞行安全的重要因素之一，尤其是活动范围多为低空的通用航空，由于飞行高度低、机载设备简陋、且多为轻小型飞机，其抵御复杂天气条件的能力明显弱于运输航空飞机，这也是国内外通用航空事故率远高于运输航空的重要原因之一。另外，航空业对天气有着更为严苛的标准要求，以飞行训练为例，大众认为的晴好天气却不一定适合学员执行单飞科目，因加热不均引起的低空湍流会严重干扰飞行状态的保持，甚至会威胁到训练安全。通常是在单飞开飞前安排有经验的机组起飞测试湍流情况，但由于低空气象条件变化快，该方法只能作为临时措施，这也反映出目前气象部门的精细化服务能力无法充分满足飞行训练的保障需求。结合气象学中对边界层高度的定义可知，低空空域与边界层相重合。一方面，边界层顶是湍流与非湍流运动的分界面，其高度的变化反映了边界层内湍流混合和对流发展等物理过程^[1]，这一系列物理过程的变化，会引起飞机颠簸和风切变; 另一方面，边界层内污染物、水汽的扩散和输送很大程度上依赖边界层的结构和厚度，而污染物、水汽的聚集会影响地面和空中能见度，进而影响飞行安全。因此，细致刻画机场区域边界层高度对提升航空气象保障能力有重要意义。

边界层高度存在明显的日变化特征和空间变化特征^[2]，也是判断湍流混合、大气污染物扩散和分析大气环境容量的主要特征量^[3]。边界层高度无法直接通过观测获得，需要借助气象要素(温度、风速和湿度等)的垂直剖面进行判定^[4]。所以在实际应用中，边界层高度的反演方法主要从湍流、热力、动力、风速和物质分布角度进行反演^[5]。从湍流运动角度反演边界层高度，通常将湍流能量接近消失的高度视为边界层顶^[6]。由于湍流能量的垂直分布观测困难，所以在实际应用中多从热力和物质垂直分布来反演边界层高度^[7]；从热力角度反演边界层高度，通常将温度梯度明显不连续定义为边界层高度^[6]；动力角度通常将风速逼近的转风高度定义为边界层高度^[8]。从物质角度计算边界层高度，通常是根据示踪物(水汽、气溶胶等)在空间中分布，确定边界层的高度^[7]。获取温度、湿度、风速等要素廓线资料的手段有很多，例如无线电探空^{[1, 9]}、飞机观测^[10]等，但这些设备具有一定的时空局限，无法提供连续的观测。随着遥感技术近年来不断的发展，观测手段逐渐增多，如激光雷达(light detection and ranging, LiDAR)^[11]、微波辐射计和风廓线雷达等，利用这些探测数据能够很好地确定边界层高度^[12-13]，其中激光雷达由于设备轻巧、时间和空间分辨率高、探测空间范围大、刷新率快的优点，已成为了国际公认的检测边界层高度的新技术手段^[14]。

利用激光雷达对气溶胶的后向散射信号反演边界层高度，国内外学者提出了多种识别方法。WANG等人^[15]基于波长为532 nm的米散射激光雷达数据，采用梯度法进行了计算，结果显示, 北京夏季边界层日变化特征明显。LI等人^[16]基于激光雷达数据，采用小波协方差变换的方法，计算了边界层高度，结果发现该方法只适用于过渡层比较明显的情况。MENUT等人^[17]基于激光雷达数据，针对巴黎一次重污染过程，利用标准方差法确定了巴黎地区城市的边界层高度。除此以外，基于激光雷达探测技术还发展了目测法^[18]、阈值法^[19]、理想廓线法^[20]以及基于这些方法的衍生算法。总体而言，大多数方法都能准确地反演边界层高度，但针对不同情况下的廓线结构，都有自己的局限性和适用范围。本文中研究利用多普勒激光雷达于机场开展边界层高度的反演，通过无线电探空仪和飞机数据验证了该雷达反演边界层高度的可行性，有助于加深对机场边界层演变特征认识，提升多普勒激光雷达的利用率，可为航空高影响天气的预警预报提供有力支撑。

本文中使用WINDCUBE V2 3-D扫描式多普勒激光雷达，发射的激光脉冲波长为1540 nm，扫描径向线上的距离分辨率为50 m，探测盲区为200 m以下，理论探测距离200 m~6000 m，包括距离高度指示器(range height indication, RHI)、平面位置指示器(plane position indication, PPI)。该雷达采用的是外差法的多普勒频移激光探测和测距技术，探测得到的数据参数包括时间、方位角、仰角、距离、径向风速、载噪比(carrier-to-noise ratio, CNR)、谱宽、可信度等, 系统原理如图 1所示。图中，CW(continuous wave)为连续波，ADC(analog-to-digital converter)为模数转换器。

Figure 1.  Schematic of a scanning Doppler LiDAR system

边界层内的气溶胶粒子的分布浓度明显高于边界层之上的“自由大气”部分的气溶胶浓度^[21]，气溶胶浓度越大，信号越强；反之，信号越弱。基于此原理，本文中使用了RHI扫描模式下的90°朝天顶方向扫描接收到的载噪比数据，载噪比的理论公式可表示为^[22-23]：

式中，η为系统效率，E为激光器发射能量，β为气溶胶后向散射系数，λ为波长, T为传输R处的大气透射率, R为探测距离，D为望远镜孔径, ℏ为普朗克常量, B为探测电路带宽。

作为对比，本文作者还使用了温江站的L波段雷达探空数据。探空数据包括温度、湿度、风速、位温、压强等参数。探测时次(北京时间)为每天3次，分别为02、08和20时刻，采样间隔小于1.5 s。温江站与广汉机场都位于成都平原，直线距离约50 km，边界层高度在一定范围内差距较小，可以用温江站的探空反映广汉机场地区的边界层特征。为了验证及补充探空资料时次少的问题，本文作者还使用了飞机传感器获得的1次/s的高度、温度数据。由于训练飞行高度所限，只能获得1200 m以下的温度廓线信息。

为了了解当天的天气情况，本文作者还使用了广汉机场气象台提供的逐小时的观探测资料，包括自动观测的气温、相对湿度、露点温度和风速，以及人工观测的云量、云底高信息。

从热力作用角度，边界层高度可以定义为位温或温度梯度变化不连续的高度^{[6, 24]}。本文中依据上述定义，利用位温廓线数据，将位温明显不连续的高度定义为边界层高度。

利用激光雷达CNR数据反演边界层高度之前，首先对雷达数据进行了质量控制^[7]，以剔除不可信数据，然后运用梯度法、小波协方差法和标准偏差法计算了边界层高，3种反演方法具体见下。

(a) 梯度法。激光雷达的距离平方校正信号R_CNR(Z²)的梯度D(Z)可表示为：

式中, Z表示高度。为了消除气溶胶浓度微小变化造成的梯度值大幅变化的问题，用Savitzjy-Golay方法对每条CNR数据进行了平滑处理后再求取D(Z)，最后将D(Z)的负值最小值确定为边界层高度^[25]。

(b) 小波协方差法。小波协方差变换是用于检测信号突变的方法，当函数值越大时，表明信号函数与小波函数相似性越高^[26]。小波协方差法多使用哈尔小波作为小波基，哈尔小波基由一个复合函数构成：

式中，a为空间范围或者计算步长，b为哈尔小波基的函数中心，利用哈尔小波函数定义的协方差如下：

式中，f(Z)是后向散射信号，Z_t和Z_b分别是信号高度的上下限。经整理可得：

小波协方差变换函数W(a, b)反映的是在高度b±a/2的范围内f(Z)与哈尔函数的相似程度。所以，当W(a, b)越大，表明f(Z)与哈尔函数越相似，即其阶跃变化越明显，由此可将W(a, b)取得最大值的高度作为边界层高度。

(c) 标准偏差法。标准偏差(standard deriation, STD)反映了激光雷达后向散射信号在某高度的离散程度，该值越大, 离散性越强。由于在自由大气与边界层的交界处总是存在强烈的夹卷，信号在边界层顶处会存在剧烈的信号变化，因此可将标准偏差最大值的高度视为边界层高度^[27]。标准偏差的公式可表示为：

式中, N表示数据点。本文作者将D_STD(Z)取得最大值的高度定义为边界层高度。

查阅广汉机场逐时天气纪要后，选取了天气晴好且激光雷达探测效能高的2020-07-19进行对比分析。当天的日出时间为06:12，日落时间为20:06。图 2展示了逐时的云量H_b(见图 2a)、温度T、露点温度T_d、相对湿度h_r(见图 2b)以及风速v(见图 2c)。图中黑线表示有L波段探空, 虚线表示有飞机探测。从广汉机场观探测资料可以看到，07-19夜间，广汉机场云量接近8分(八分制云量，即满天云)，但为云底较高的中云。从日出前的05时(本文中所用时间均为北京时间)起至日落时分的20时，天空为无云或少量中云，云底进一步抬高，期间1000 m以下有少量低云。日落后，又变为满天云，云底与日出前一样(见图 2a)。结合温度来看(见图 2b)，05时气温达到全天最低值，为19.9 ℃，相对湿度接近100%。日出后，地面温度上升，空气含水能力逐渐增大，导致相对湿度下降。至16时，气温达到全天最高的34.2 ℃，而后气温开始缓慢下降。结合云量来看，天空无云的06时~08时气温上升更快。随着地面升温，热力对流增强，地面风速也略有上升，08时~14时的地面风速为2 m/s，而其它时刻的风速接近不大于1 m/s的静风(见图 2c)。

Figure 2.  The height of cloud bottom H_b, temperature T, dew point temperature T_d, relative humidity h_r and wind velocity v of airport weather

从02时温江站探空曲线可以看到，因夜间持续向外放出长波辐射而降温，100 m左右形成了逆温，此外，在1500 m附近还存在另一个明显的逆温层(见图 3a)。从位温能够更加清晰地看到(见图 3b)，300 m以下大气为稳定层结，300 m~1500 m大气更接近于中性，1500 m处大气稳定度又再次增强。由大气稳定度变化可以判断出稳定边界层高度约为300 m左右，300 m~1500 m之间为白天对流边界层遗留形成的残余层，残余层顶高约为1500 m。受稳定边界层的影响，300 m上下相对湿度有明显差异，边界层内相对湿度随高度增加而减小，而在边界层之上，相对湿度随高度增加而增大。日出后的08时，由于太阳短波辐射增强，地面迅速增温，大气由稳定边界层向不稳定边界层转化。近地面逆温层被破坏，逆温层被抬升至100 m~300 m之间(见图 3a)。对比08时前后的飞机探测数据来看，尽管广汉机场温度整体高于温江，但二者曲线形状吻合，这表面利用温江探空数据能够反映广汉机场边界层特征。位温廓线上可以看到(见图 3b)，地面附近更稳定，随着高度增大，大气稳定度减小，约在1000 m以上，大气层结接近中性，同样相对湿度在1000 m~1200 m上下发生了明显变化(见图 3c)。由此可以确定08时的热力边界层高度在1000 m~1200 m之间。从其它时次的飞机探测曲线可以看到，随着地面气温的进一步升高，近地面的逆温结构完全消失，对流边界层快速发展，但由于训练飞行的高度限制，无法获得对流边界层顶的具体高度。20时的温度和位温廓线与02时非常相似，但此时近地面形成的稳定层结高度要低于02时，高度不足200 m。对比20时飞机的温度廓线可以看到，广汉机场逆温层高度要略高于温江站，这与广汉机场地面气温更高有关。另外，由于太阳刚落山不久，白天对流边界层残余部分的影响依然显著，相较于08时，20时残余层高度更高，达到了2000 m左右。

Figure 3.  Wenjiang sounding data and aircraft data

图 4展示了从激光雷达测得的随时间变化的气溶胶浓度垂直分布。红色线为小波协方差法计算的边界层高度，蓝色点和红色点分别为梯度法和标准偏差法计算的边界层高度。02时前后，由于稳定边界层的存在，大气湍流被抑制，大部分气溶胶粒子被限制300 m以下，该高度内回波信号强烈。此外，在500 m~1000 m之间还有一层强度更弱的回波信号区，这是残余层保持的对流边界层混合到该层的气溶胶粒子回波。以上回波特征与探空资料所反映的稳定边界层和残余层特征一致。而从20时前后的回波可以看到，500 m以下有一个强回波区，1000 m左右形成了第2个强回波区，与02时不同，下层的回波强度要弱于上层。这是由于太阳刚落山，地面稳定边界层刚开始形成，地面气溶胶刚开始积累，浓度较低。而在稳定边界层之上，由于白天湍流活动充分发展，气溶胶粒子被输送至更高高度，并充分混合，使得日落后残余层依然比较厚，气溶胶粒子浓度较近地面更高。通过对比还能看出，在气温快速上升的06时~15时，由于逆温层被破坏，热力湍流快速发展并向上扩展，使得气溶胶粒子的层顶逐渐抬高。由于气溶胶粒子被扩散至更大空间，浓度下降使得回波强度较其他时刻弱。综合以上分析可知，随着边界层高度的变化，边界层内的气溶胶粒子扩散高度、浓度都有明显的日变化。

Figure 4.  Fill plot of LiDAR CNR data

从反演结果来看，梯度法、小波协方差法和标准偏差法都能捕获到边界层的变化特征，但对边界层的不同发展阶段，3种方法识别的边界层高度有一定的差异。从反演结果看，三者都能较好地捕获中午前后对流边界层高度的变化特征且一致性较好，三者的差异主要表现在对日出前和日落后稳定边界层的识别上。日出前，由于稳定边界层高度过低，小波协方差法受波长限制，无法识别出高度较低的稳定边界层顶，而将2000 m左右的残余层顶识别为了边界层顶。03时~04时，随着稳定边界层的抬升，小波协方差法也能够识别出稳定边界层，使得小波协方差法识别的边界层高度在残余层顶和稳定边界层顶不断跳变。比较来看，梯度法和标准方差法则能够识别出稳定边界层逐渐抬升的过程，但标准偏差法的连续性要低于梯度法。日出后至日落前，小波协方差法和梯度法识别的边界层顶高度的变化特征比较一致，区别在于小波协方差法识别的高度要低于梯度法。15时之前，标准偏差法识别的边界层高度与另外两种方法接近。但在15时后，随着边界层顶高度的降低，近地面附近气溶胶浓度随时间变化明显，标准偏差法会间断性地将这种变化识别为边界层高度，使得边界层高度发生跳变。日落后，随着近地面稳定边界层的形成，小波协方差法已不能识别出其高度，而梯度法依然能够反映出稳定边界层的高度。值得注意的是，此时空气中气溶胶的上下分布不同于日出前，梯度法也能在个别时间识别出残余层顶高度。

结合探空数据识别的边界层高度来看，02时梯度法确定的稳定边界层高度略高于探空确定的高度。08时3种方法观测的边界层顶高度约为1400 m~1600 m之间，高于温江站的1000 m~1200 m。20时前后，梯度法和标准方差法识别的稳定边界层高度约400 m左右，高于探空数据得到的边界层高度。梯度法和小波协方差法识别的残余层顶高度约为1500 m左右，要低于探空数据确定的高度。结合飞机探测数据来看，由于广汉机场多个时间点低空的气温要高于温江站，故激光雷达获得的边界层高度与探空得到的边界高度存在差异可能是实际情况的反映。

总体比较来看，标准方差法识别边界层高度的时间分辨率低，特别是午后识别出的边界层顶与实际有明显偏差。小波协方差法对日出前和日落后高度较低的稳定边界层顶无法识别。梯度法无论是准确性、连续性还是稳定性上都表现出了明显优势。

从经典边界层理论来看，边界层高度的升高主要是源于太阳辐射对地面加热，近地层产生的热通量增加使得湍流混合加强，导致了边界层高度的上升^[20]。从垂直运动可以看到(见图 5a)，从02时起，出现了上升和下沉运动交替出现，这种运动由地面延伸至边界层顶，表现出重力波的特征。随着时间的推移，升降运动的强度增加，垂直运动的范围也变得更为深厚，边界层高度也随之抬升。日出前后，垂直运动的强度最强，可达1 m/s，厚度也达到了该时段的最高值。谱宽数据能够反映小尺度的大气湍流情况，由图 5b可以看到, 日出前边界层内垂直运动活跃期，湍流活动非常旺盛，这也反映了动力性湍流对日出前边界层抬升有重要贡献。07时以后，随着稳定边界层逐渐被破坏，重力波维持的条件不复存在，对流边界层快速发展，原来深厚的升降气流，变为在垂直方向上交替出现，表明由于热力湍流加强，原来深厚的垂直运动已破碎为更小尺度的湍流。11时以后，边界层内热量和动量混合均匀，边界层内以上升运动为主。值得注意的是，在边界层顶上下湍流活动有明显差异，这是由于对流边界层在发展的过程中会伴随大量的湍流运动，这些湍流运动与自由大气的混合形成了夹卷，而在边界层顶以下，湍流活动明显减弱。经比较可以发现，湍流强度的分界和梯度法识别的边界层顶高度相吻合，进一步说明了梯度法识别边界层高度的优势。

Figure 5.  Analysis of turbulent characteristics in boundary layer

从12°仰角PPI(见图 6a)和RHI(见图 6b)的径向风来看，日出前在1000 m左右，存在7 m/s以上的低空西南风大风，配合低空逆温层，为重力波的形成提供了有力条件。正是由于这种深厚的垂直运动存在，使得动力性湍流不断向上发展，边界层高度也不断抬升。图中蓝色表示靠近激光雷达，红色表示远离激光雷达。

Figure 6.  Radial wind chart

本文作者通过利用温江站L波段的探空数据、飞机数据和多普勒激光雷达数据对广汉机场上空的边界层高度进行了反演和对比分析，并通过激光雷达对机场上的边界层进行了持续性的观测。

(a) 对比探空仪和飞机的低空数据来看，广汉机场的温度整体要略高于温江，但二者曲线形状相吻合，说明了温江站的探空数据对广汉机场边界层高度的反演具有一定的参考价值。

(b) 从探空数据确定的边界层高度与激光雷达数据持续观测到的机场上空边界层高度的对比分析来看，两者反演的结果基本一致，这验证了该多普勒激光雷达反演边界层高度的可行性。另外，机场上空当天边界层高度呈现出明显的日变化特征，其中对流边界层和残余层高度基本维持在2000 m左右，在日出前和日落后稳定边界层高度在100 m~200 m。

(c) 从不同方法反演的边界层高度来看，3种方法都能有效估算对流边界层高度，且变化趋势基本一致，但针对稳定边界层的反演，小波协方差法由于受到步长限制，无法识别稳定边界层；标准偏差法虽然能够在一定程度上识别出稳定边界层，但连续性较差。相比之下，梯度法无论是在对流边界层的识别上，还是在稳定边界层的识别上都优于小波协方差法和标准偏差法。

(d) 从多普勒激光雷达的低空风数据来看，日出前物质边界层的抬升主要是由于在重力波的影响下地面垂直气流运动的发展。另外，从多普勒激光雷达的谱宽数据能够看出，湍流强弱的分界与梯度法反演的边界层顶一致性较好。

本次观测中只选用了放单考试当天机场的数据，对边界层高度进行了探测研究，后续会针对不同的天气背景对边界层进行研究，探究不同条件下边界层的变化特征, 以便能够为训练中的飞行学员提供一些预警信息，从而更好地保障飞行安全。