基于Logistic大气消光廓线识别混合层厚度方法与流程

本发明属于大气和环境技术领域，尤其涉及一种基于Logistic大气消光廓线识别混合层厚度方法。

背景技术：

大气边界层是地气之间物质和能量交换的桥梁，对其结构特征及变化成因的深入认知是天气气候模式精度提升的基础。作为大气边界层的一种重要表现形式，混合层定义为湍流特征不连续界面以下的大气，表征污染物在垂直方向被热力湍流稀释的范围，即低层空气热力对流与动力湍流输送所能达到的高度，相关研究一直是大气科学和环境科学共同关注的前沿问题(胡非，2003；杨辉，2006；丁辉，2015)。Nozaki等早在1973年就提出了一种利用地面常规气象资料估算混合层厚度的方法(程水源，1997)；随着计算机水平的飞速发展，数值模式已经成为模拟边界层复杂物理过程的重要工具(Lilly，1968；Betts，1973；蒋维楣，2010；伍大洲，2006)。此外，由于Mie散射激光雷达具备对大气光学特性垂直分布进行高时空分辨率探测的优点，近年来已成为反演大气边界层结构及颗粒物垂直分布的另一手段(何涛，2013)。就目前常用的反演算法而言，梯度法把激光雷达距离订正信号(RCS)的梯度最小值所在的高度作为边界层高度(Hoff R M，1996)；拐点法则是将混合层顶高定义为滤波处理后的距离订正信号对高度的二阶导数最小值所在高度(Menut L，1999)；SBH99算法通过将实测后向散射系数的垂直分布拟合成理想曲线，据此确定混合层厚度(Steyn，1999)；Brooks(2003)等基于复合阶跃函数(Haar函数)提出了一种用于探测边界层信号阶跃变化的方法，即为小波协方差变换法。虽然Mie散射激光雷达已被广泛应用于边界层结构的诊断，但仍困扰于上述算法稳定性和普适性差以及人为任意性强等不足，例如：梯度法和拐点法均容易受到信号局部结构的影响，对信噪比要求较高，反演结果受噪声干扰影响较大(王琳，2012；李红，2015)；小波协方差变换法、SBH99算法对自身输入参数较为敏感(贺千山，2005；王琳，2012)。

综上所述，现有大气光学特性垂直分布进行高时空分辨率探测存在稳定性和普适性差以及人为任意性强，对信噪比要求较高，反演结果受噪声干扰影响较大，对自身输入参数较为敏感。目前利用激光雷达反演混合层厚度的方法均是建立在边界层消光系数(或后向散射系数)的垂直分布特征的基础之上，由于边界层内部活动复杂多变以及激光雷达受噪声等因素的干扰，导致探测结果出现一定的波动，在此情况下，对于混合层厚度的确定带有很强的主观色彩，稳定性和普适性得不到保证。

技术实现要素：

本发明的目的在于提供一种基于Logistic大气消光廓线识别混合层厚度方法，旨在解决目前利用激光雷达反演混合层厚度的方法均是建立在边界层消光系数(或后向散射系数)的垂直分布特征的基础之上，由于边界层内部活动复杂多变以及激光雷达受噪声等因素的干扰，导致探测结果出现一定的波动，在此情况下，对于混合层厚度的确定带有很强的主观色彩，稳定性和普适性得不到保证的问题。

本发明是这样实现的，一种基于Logistic大气消光廓线识别混合层厚度方法，所述基于Logistic大气消光廓线识别混合层厚度方法包括以下步骤：

步骤一，基于Mie散射激光雷达反演的大气消光系数廓线，确定混合层顶附近自下而上大气消光系数历经的急剧下降区、总体缓变区以及二者之间的过渡区；将消光系数和垂直高度分别作为横纵坐标，构建直角坐标系；

步骤二，用半程Logistic曲线拟合研究区内消光系数随高度的变化特征，构建半程Logistic曲线；

步骤三，根据最小二乘原理求解半程Logistic曲线参数，求解半程Logistic曲线曲率最大值对应高度。

进一步，构建直角坐标系，记急剧下降区的起始点为(σ₀，z₀)。

进一步，构建半程Logistic曲线为：

其中，z表示研究区高度，σ₀、z₀为研究区下端的消光系数以及对应的高度，r是半程Logistic曲线待定参数，与消光系数垂直递减速率有关。

进一步，混合层厚度的表达式为z_m：

进一步，所述基于Logistic大气消光廓线识别混合层厚度方法包括以下步骤：

(1)基于Mie散射激光雷达反演的大气消光系数廓线，确定混合层顶附近自下而上大气消光系数历经的急剧下降区、总体缓变区以及二者之间的过渡区，考虑到激光雷达探测盲区的影响，记急剧下降区的起始点为(σ₀，z₀)；

(2)用半程Logistic曲线拟合研究区内消光系数随高度的变化特征，见下式：

其中，z表示研究区高度，σ₀、z₀为研究区下端的消光系数以及对应的高度，r是半程Logistic曲线待定参数，为Logistic大气消光廓线；

(3)根据最小二乘原理求解半程Logistic曲线参数r，进一步求解半程Logistic曲线曲率最大值对应高度，由此得到混合层厚度z_m的表达式为：

本发明提供的基于Logistic大气消光廓线识别混合层厚度方法，根据实例分析表明，相对于基于Mie散射激光雷达计算混合层厚度的现有方法而言，半程logistic曲线拟合新算法相应的计算结果能合理地反映边界层结构的日变化特点及其污染效应，且分析精度更高。

本发明基于成都市2013年6月至2014年2月Mie散射激光雷达探测数据的分析表明，混合层顶附近大气消光系数自下而上均存在急剧下降区、总体缓变区以及二者之间的过渡区。利用半程logistic曲线对上述变化特征进行拟合，通过计算该曲线曲率最大值对应的高度，据此提出Mie散射激光雷达识别混合层厚度的新算法。针对新算法的进一步分析表明：(1)半程logistic曲线拟合算法着眼于混合层附近大气消光系数垂直演变的普遍规律，受人为因素和信噪比影响小，计算的混合层即为湍流特征不连续界面以下湍流充分发展的气层，具有物理意义明确、普适性好、客观性强以及计算稳定等优点。(2)针对不同空气质量等级样本的分析表明，半程logistic曲线拟合算法与位温探空曲线计算的混合层厚度总体一致，同时也取得了比梯度法、拐点法和SBH99算法更好的效果。(3)将半程logistic曲线拟合算法应用于成都地区一次重霾污染过程的诊断，准确地揭示了区域在该时段主要的污染气象特征，计算的混合层厚度与地面颗粒物质量浓度保持了高度的相关。

附图说明

图1是本发明实施例提供的基于Logistic大气消光廓线识别混合层厚度方法流程图。

图2是本发明实施例提供的混合层厚度求解的示意图。

图3是本发明实施例提供的6种不同空气质量状况下新方法和位温法计算的混合层厚度对比示意图；

图中：(a)2014年1月31日07时<严重>；(b)2013年12月2日07时<重度>；(c)2013年11月30日07时<中度>；(d)2013年11月26日07时<轻度>；(e)2013年8月8日07时<良>；(f)2013年7月10日07时<优>。

图4是本发明实施例提供的成都市2014年1月21日-2月5日AQI及地面PM2.5质量浓度逐日变化示意图。

图5是本发明实施例提供的成都市1月27日-1月31日大气消光系数廓线、混合层高度和PM2.5质量浓度逐时变化示意图。

图6是本发明实施例提供的污染时段内(1月23日-2月4日)平均混合层高度与地面PM2.5质量浓度的日变化示意图。

图7是本发明实施例提供的污染时段内平均混合层高度与PM2.5质量浓度散点示意图。

具体实施方式

为了使本发明的目的、技术方案及优点更加清楚明白，以下结合实施例，对本发明进行进一步详细说明。应当理解，此处所描述的具体实施例仅仅用以解释本发明，并不用于限定本发明。

下面结合附图对本发明的应用原理作详细的描述。

如图1所示，本发明实施例提供的基于Logistic大气消光廓线识别混合层厚度方法包括以下步骤：

S101：基于Mie散射激光雷达反演的大气消光系数廓线，确定混合层顶附近自下而上大气消光系数历经的急剧下降区、总体缓变区以及二者之间的过渡区；将消光系数和垂直高度分别作为横纵坐标，据此构建直角坐标系；

S102：用半程Logistic曲线拟合研究区内消光系数随高度的变化特征，据此构建半程Logistic曲线；

S103：根据最小二乘原理求解半程Logistic曲线参数，进一步求解半程Logistic曲线曲率最大值对应高度。

下面结合具体实施例对本发明的应用原理作进一步的描述。

本发明利用成都市2013年6月至2014年2月Mie散射激光雷达的探测数据，通过对边界层大气光学特征垂直演变共性的挖掘，结合物理意义的分析和多重实例的验证，据此提出了基于半程Logistic曲线拟合的混合层厚度识别新算法。

1资料与方法

1.1资料来源

本发明采用的Mie散射雷达数据与颗粒物浓度资料由四川省环境监测总站提供。激光雷达观测仪器为北京怡孚和融科技有限公司生产的EV-激光雷达激光雷达(偏振Mie散射微脉冲激光雷达)，探测点设置于西南交通大学九里堤校区土木馆楼顶，距地高度35m，于2013年6月正式投入使用，可实现全天自动数据采集，数据采集频率为1次/3min，主要观测范围可达30km，探测盲区为60m，测距分辨率为15m；颗粒物质量浓度使用美国热电子环境设备公司生产的大气污染物环境监测分析仪(1400a环境颗粒物检测仪)进行观测，监测点位于成都市人民南路四段，数据采集频率为1次/5min，精度为±1.5μgn。用于对比的探空数据则来源于成都市温江气象站，数据为一日两次(探测时间分别为北京时间07:00、19:00)的探空观测资料，其垂直空间上的分辨率为5m，主要包含气温、气压、湿度和露点温度等气象要素。

激光雷达与颗粒物质量浓度监测点均位于主城区，两者距离10km左右，因而两地之间颗粒物浓度差异很小。气象探空站与激光雷达布设位置距离不超过20km，在当今城市化进程迅猛的背景条件下，二者对边界层的探测结果虽受到下垫面差异的影响，但在一定的误差范围内仍具有可比性。

1.2方法介绍

1.2.1原理

大气边界层(ABL)也称为行星边界层(PBL)，是大气与下垫面直接发生相互作用的层次，从地面向上厚度可达几百至上千米(胡非，2003)。当大气层结处于不稳定时，从贴地层向上的大气湍流活动强烈，处于这种运动状态的边界层称对流边界层(CBL，也称混合层)。在自由大气与边界层间有一个过渡区域(卷夹层)，从边界层中受热上升的气块可以穿透卷夹层进入自由大气，湍流重力波亦可使自由大气中具有较高位温的气块进入边界层。不稳定边界层和中性边界层对应的ABL高度位于逆温层的层底，而稳定边界层的ABL高度则位于逆温层的层顶。受边界层热力和动力因素的共同影响，湍流扩散导致颗粒物浓度由下往上呈现出明显的递减趋势(Wu，2002)。由于受逆温“顶盖”的制约，上述颗粒物浓度显著递减的趋势会在混合层顶附近出现转折，由此产生颗粒物浓度变化特征在垂直方向上的不连续现象。边界层以上颗粒物浓度含量骤减，气体分子的消光逐渐占据主导地位，大气消光因而表现为总体缓变的特征。颗粒物浓度急剧下降区和总体缓变区本质上体现的是两种不同的湍流状态，这种状态的转折对应于消光系数曲率的最大值，因而可将消光系数曲率的最大值对应的高度视为混合层顶。

综上，尽管边界层气象条件呈现出非常复杂的表现形式，混合层顶附近自下而上大气消光系数均存在急剧下降区、总体缓变区以及二者之间的过渡区，整体表现出与半程Logistic曲线一致的变化特征。因此，利用半程Logistic曲线拟合消光系数的上述变化特征，进一步求取消光系数曲率最大点对应的高度，据此提出基于Mie散射激光雷达识别混合层厚度的新算法。

1.2.2流程

(1)基于Mie散射激光雷达反演的大气消光系数廓线，确定混合层顶附近自下而上大气消光系数历经的急剧下降区、总体缓变区以及二者之间的过渡区。将消光系数和垂直高度分别作为横纵坐标，据此构建直角坐标系，记急剧下降区的起始点为(σ₀，z₀)，见图2。

(2)用半程Logistic曲线拟合研究区内消光系数随高度的变化特征σ_z，据此构建半程Logistic曲线如下：

其中，z表示研究区高度，σ₀、_z0为研究区下端的消光系数以及对应的高度，r是半程Logistic曲线待定参数，与消光系数垂直递减速率有关。

(3)根据最小二乘原理求解半程Logistic曲线参数r，进一步求解半程Logistic曲线曲率最大值对应高度，由此可得混合层厚度的表达式为z_m：

图2为实际消光系数曲线以及应用半程Logistic曲线拟合得到的理想曲线。图中黑色虚线以上部分为拟合区所在，起始点为(σ₀，z₀)，根据公式(1)和(2)将该区域消光系数廓线分布特征拟合得到理想曲线，并进一步求解半程Logistic曲线曲率最大值对应高度，即红色虚线表示高度，其上为消光系数总体缓变区，其下为消光系数急剧下降区。

2实例验证

2.1新算法与位温探空曲线计算混合层厚度的对比分析

依据空气质量指数(AQI)的大小，在研究时段内选取优、良、轻度污染、中度污染、重度污染、严重污染六类不同的样本，对logistic曲线拟合新算法与位温探空曲线计算的混合层厚度(马金，2011)进行比对，据此分析新算法的适用性，结果见图3。以图3(a)为例，该图反映的是2014年1月31日(空气质量状况为严重污染，对应AQI为432)07:00时的位温曲线与消光系数廓线。从图中可以看出，从地面至239m位温虽随高度增大而有所增加，但增长率较小，表明该气层处于中性偏稳定的状态；从239m到295m位温增长率较大，气层转为稳定状态，由于混合层之上通常会覆盖一层稳定的气层(逆温层)，据此可判断混合层厚度约为239m。作为比较，该图也给出了同一时刻Mie散射激光雷达反演的大气消光系数垂直廓线，利用logistic曲线拟合新算法得出对应该时刻的混合层厚度为224m；进一步分析表明，大气消光系数随高度变化趋势在此高度附近发生明显转折，其上下分别对应于为整体缓变区和急剧下降区。针对图3中其余5类样本的分析也同时表明，两种方法得到的混合层厚度基本一致。

为进一步验证新算法的计算精度，按不同污染等级各自等量增选2个样本，分别利用梯度法、拐点法和SBH99算法计算了对应的混合层高度，见表1。由表1可见，以位温法计算的混合层高度作为判定的依据，新算法取得了最为精准的结果。

表1利用Mie散射激光雷达相关算法与位温法计算混合层高度的结果对比(单位：m)

注：时间分别为每日07:00时，相对误差均以位温法结果为标准。

综上分析，以位温探测结果作为混合层高度的判定标准，梯度法和拐点法易受单点信息的干扰，往往不能很好地反映边界层的真实演变过程；SBH99算法由于对自身输入参数较敏感，计算不稳定；而新算法基于对混合层附近大气消光系数垂直变化共性的挖掘，计算结果具有物理意义明确、客观性强以及精度高等优点。

2.2新算法分析一次灰霾过程混合层厚度的演变

2014年1月23日～2月4日四川盆地出现了一次持续时间长达13天的区域性灰霾天气过程，其中以成都市的污染程度最为严重。根据四川省环境监测总站监测资料的分析表明，成都地区其间严重污染天气出现了5天、重度污染6天、中度污染2天，其中1月27日到31日连续5天为严重污染，日均AQI序列演变如图4所示，首要污染物均为PM2.5。针对1月27日～1月31日天气状况的进一步分析表明，在该时段无雨以及静风的条件下，受热力调控的混合层厚度无疑是地面颗粒物变化最重要的气象成因。以此次严重污染过程为研究对象，利用半程logistic曲线拟合新算法计算其逐时混合层厚度，如图5所示。由图5可见，该时段内混合层厚度存在一定的日变化，表现为白天高和夜晚低的特点，但整体仍维持相对较低的一种状态，平均高度约为275m，其中1月31日对应的混合层厚度只有190m左右。

混合层的上述特征使得夜晚产生的污染物在白天无法得到很好的稀释扩散，而白天形成的污染物会进一步加剧夜间的污染程度。缺乏湿清除以及水平和垂直扩散能力差是此次重霾过程的主要污染气象特征，也是颗粒物得以渐进累积的气象背景。

图5同时描述了此次严重污染时段地面PM2.5质量浓度的逐时演变过程，易见混合层厚度和地面细颗粒物质量浓度存在明显的反位相关系。为深化对这一问题的认知，图6给出了研究时段内(1月23日-2月4日)逐时平均混合层高度和逐时平均细颗粒物质量浓度的日变化。由图6可知，混合层高度的日变化基本表现为单峰单谷型，从00:00时至09:00时，混合层均处于较低状态，期间有谷值出现；自09:00时开始，随着太阳辐射的增强，地面逐渐升温，下垫面与大气之间的相互作用得到加强，大气层结由夜间的稳定状态向不稳定状态过渡，混合层抬升，污染物在垂直方向上扩散能力得以加强；午后对流活动愈加旺盛，混合层高度抬升速率明显增加，并于18:00时前后达到峰值；之后，大气层结向稳定状态过渡，混合层高度也随之逐渐降低。与上述混合层的演变过程相对应，夜间至清晨，受低混合层高度的扩散制约，地面PM2.5质量浓度在此时段内均维持相对较高的状态；此后，因混合层高度快速抬升的作用，地面细颗粒物质量浓度随之迅速降低，并于19:00时左右达到最低值；在之后混合层高度波动下降的背景下，细颗粒物质量浓度逐渐上升。进一步分析表明，该污染时段内逐时平均混合层高度与逐时平均PM2.5质量浓度的相关系数R为-0.67(通过置信水平a＝0.05的显著性检验)，如图7所示。由上可知，尽管此次污染过程混合层高度相对较低，但仍呈现出以日为周期的起伏特征，并在很大程度上影响到地面颗粒物质量浓度的变化。

上述实例分析表明，相对于基于Mie散射激光雷达计算混合层厚度的现有方法而言，半程logistic曲线拟合新算法相应的计算结果能合理地反映边界层结构的日变化特点及其污染效应，且分析精度更高。

本发明基于成都市2013年6月至2014年2月Mie散射激光雷达探测数据的分析表明，混合层顶附近大气消光系数自下而上均存在急剧下降区、总体缓变区以及二者之间的过渡区。利用半程logistic曲线对上述变化特征进行拟合，通过计算该曲线曲率最大值对应的高度，据此提出Mie散射激光雷达识别混合层厚度的新算法。针对新算法的进一步分析表明：

(1)半程logistic曲线拟合算法着眼于混合层附近大气消光系数垂直演变的普遍规律，受人为因素和信噪比影响小，计算的混合层即为湍流特征不连续界面以下湍流充分发展的气层，具有物理意义明确、普适性好、客观性强以及计算稳定等优点。

(2)针对不同空气质量等级样本的分析表明，半程logistic曲线拟合算法与位温探空曲线计算的混合层厚度总体一致，同时也取得了比梯度法、拐点法和SBH99算法更好的效果。

(3)将半程logistic曲线拟合算法应用于成都地区一次重霾污染过程的诊断，准确地揭示了该时段内边界层的演变特征，计算的混合层厚度与地面颗粒物质量浓度保持了高度的相关。

以上所述仅为本发明的较佳实施例而已，并不用以限制本发明，凡在本发明的精神和原则之内所作的任何修改、等同替换和改进等，均应包含在本发明的保护范围之内。