气溶胶激光雷达数据质控的方法及系统

lidar network)(mckendry et al.，2011)等。气溶胶激光雷达被广泛应用于表征气溶胶垂直结构(huang et al.，2012；liu et al.，2013；sun et al.，2013；tian et al.，2018)、评估污染管控效果(yang et al.，2010)、评估大气化学传输模式(li et al.，2012，2013；wang et al.，2014)、确定边界层高度(liu et al.，2015；yang et al.，2017b)、反演气溶胶组分(burton et al.，2012；sugimoto et al.，2014)和计算气溶胶液态含水量(tan et al.，2020a，2020b)等大气污染研究中。
7.随着对气溶胶激光雷达数据定量化程序的要求越来越高，气溶胶激光雷达数据质量控制的研究也逐渐开展起来。项衍(2018)基于激光雷达硬件层面，对光机结构、电子学模块中的信号采集方式和采集累加次数，定标高度和雷达比的选取、能见度仪标定、瑞利散射标定、太阳光度计标定等方面进行了讨论。即使保证硬件质控，还需要进一步进行软件层面的质量控制。气溶胶激光雷达的数据产品是消光系数垂直廓线的时间序列，相比于地面常规观测，具有更加复杂的信息量和维度，因此对地面常规观测的质控方案无法直接应用于激光雷达上(durre et al.，2010；wu et al.，2018)。同时，常规的异常值检测方法例如孤立森立(liu et al.，2008)和k阶邻近法(ramaswamy et al.，2000)无法直接处理气溶胶激光雷达数据产品，因为检测出的异常值是代表污染烟羽的消光系数大值。earlinet的scc(single calculus chain)是一套针对激光雷达数据质控的自动化工具，包括硬件质控模块elpp(earlinet lidar pre
‑
processor)和软件质控模块elda(earlinet lidar dataanalyzer)(d’amico et al.，2015)。elda将elpp质控后的激光雷达信号作为输入来反演气溶胶光学性质，并在反演中和反演后进行软件层面的质量控制(mattis et al.，2016)。elda包含了许多经典的算法，例如蒙特卡洛不确定性分析、误差传播法，同时也包含了一些专门开发的算法，例如垂直平滑算法和时间平均算法。尽管对elda的适用性和准确性已有了一些研究(d’amico et al.，2015)，但elda仅仅是针对欧洲的激光雷达数据软件质控方案，对于具有大气复合污染(he et al.，2002)的中国的适用性还有待研究。另外，elda目前没有针对弹性散射激光雷达的平滑方案。这可能是因为elda中仅有的平滑方案是为拉曼雷达设计，拉曼雷达反演时需要对原始信号求导。除了适用于欧洲的elda方案，综合的软件层面的激光雷达数据质控研究相对较少。因此，亟需对气溶胶激光雷达软件层面的质控开展相关研究。


技术实现要素：

8.本发明要解决的技术问题是提供一种气溶胶激光雷达数据质控的方法，以从综合的软件层面对激光雷达数据进行优化。
9.为了解决上述技术问题，本发明的技术方案为：
10.一种气溶胶激光雷达数据质控的方法，包括：
11.s1，利用气溶胶激光雷达系统获取探测大气气溶胶的激光雷达数据，当数值为0～1时，进入步骤s2；当数值＞1时，将数值设为1，然后进入步骤s3；当数值＜0时，将＜0的数值以及标记数据设为0，然后进入步骤s4；
12.s2，根据雷达比和雷达信号对激光雷达数据进行蒙特卡洛不确定性分析，当存在高不确定性值时，则对高不确定性值剔除，并进入步骤s4；当存在低不确定性值时，则直接进入步骤s3；
13.s3，采用图像去噪算法对激光雷达数据进行去噪处理，当存在椒盐噪声时，剔除椒盐噪声并进入步骤s4；否则直接进入步骤s4；
14.s4，采用双边滤波图像处理算法对激光雷达数据进行平滑，进而获得质控后的消光系数。
15.进一步地，步骤s1中所述利用气溶胶激光雷达系统获取探测大气气溶胶的激光雷达数据包括：
16.s11，计算每个高度的信号均值和标准差；
17.s12，基于均值和标准差，在每个高度层建立高斯分布；
18.s13，在每个高度层的高斯分布中随机抽取n个样本，从而构建n组激光雷达信号垂直廓线；
19.s14，将这n组雷达信号连同原始雷达信号一同进入气溶胶激光雷达系统得到n+1组消光系数垂直廓线；
20.s15，通过其中的n组消光系数廓线得到每个高度层的消光系数均值和标准值，从而构建消光系数相对不确定性垂直廓线。
21.进一步地，步骤s2中：大气边界层高度设为1km。
22.大气边界层以下的相对不确定性阈值设为20％，以上的不确定性阈值设有30％。
23.进一步地，步骤s2中：雷达比被假定为50sr。
24.进一步地，雷达比的相对不确定性设为10％。
25.进一步地，步骤s3中：每个格点周围非0格点小于等于3则被去除。
26.本发明还公开了一种控制气溶胶激光雷达数据质量的系统，使用以上任一项所述的方法。
27.与现有技术相比，本发明的有益效果是：
28.本发明包含蒙特卡洛不确定性分析、图像去噪算法以及双边滤波，其中蒙特卡洛不确定性分析基于对激光雷达信号的随机采样从而区分高不确定性和低不确定性的消光系数；图像去噪算法用于去除激光雷达数据中的脉冲噪声；双边滤波对激光雷达数据进行平滑，并且能够保持数据边缘结构；本发明主要去除了位于气溶胶上部的异常值，提高了气溶胶层的平滑度，维持了原有的气溶胶垂直结构。
附图说明
29.图1为气溶胶激光雷达硬件的结构图；
30.图2为本发明一种气溶胶激光雷达数据质控的方法的流程图；
31.图3为本发明2017年1月、4月、7月和10月的消光系数垂直平均(vaec)表征的散点图，其中图3a表示transfer质控前与mca质控后vaec的关系；图3b表示transfer质控前与dn质控后vaec的关系；图3c表示transfer质控前与transfer质控后vaec的关系；图3d表示mca质控后与dn质控后vaec的关系；图3e表示mca质控后与transfer质控后vaec的关系；图3f表示dn质控后与transfer质控后vaec的关系；
32.图4为本发明2017年1月、4月、7月和10月以及4个月耦合风矢量后transfer质控前后的垂直平均消光系数比较的示意图；
33.图5为本发明2017年1月、4月、7月和10月以及4个月耦合能见度结果后transfer质
控前后的垂直平均消光系数比较的示意图；
34.图6为本发明2017年1月、4月、7月和10月以及4个月耦合rh后transfer质控前后的垂直平均消光系数比较的示意图；
35.图7为本发明2017年1月、4月、7月和10月以及4个月耦合pm
2.5
质量浓度后transfer质控前后的垂直平均消光系数比较的示意图；
36.图8为本发明2017年1月、4月、7月和10月以及4个月耦合o3质量浓度后transfer质控前后的垂直平均消光系数比较的示意图；
37.图9为本发明2017年7月1日至3日nies激光雷达系统反演的消光系数时间
‑
高度截面图原始数据(a)，经过mca质控后(b)、经过dn质控后(c)、经过bf质控后(d)和质控增量(e)的效果图；
38.图10为本发明2017年7月1日至3日的消光系数垂直平均(vaec)表征的散点图。
具体实施方式
39.下面结合附图对本发明的具体实施方式作进一步说明。在此需要说明的是，对于这些实施方式的说明用于帮助理解本发明，但并不构成对本发明的限定。此外，下面所描述的本发明各个实施方式中所涉及的技术特征只要彼此之间未构成冲突就可以相互组合。
40.本发明公开了一种气溶胶激光雷达数据质控的方法，是一套专门针对气溶胶激光雷达数据质量控制的集成方案，也可以称为transfer(a data quality assurance and control scheme for aerosol lidar，transfer)，如图2所示，主要包含蒙特卡洛不确定性分析和双边滤波，具体包括以下步骤：
41.s1，利用气溶胶激光雷达系统获取探测大气气溶胶的激光雷达数据，当数值为0～1时，进入步骤s2；当数值＞1时，将数值设为1，然后进入步骤s3；当数值＜0时，将＜0的数值以及标记数据设为0，然后进入步骤s4；
42.s2，根据雷达比和雷达信号对激光雷达数据进行蒙特卡洛不确定性分析，当存在高不确定性值时，则对高不确定性值剔除，并进入步骤s4；当存在低不确定性值时，则直接进入步骤s3；
43.s3，采用图像去噪算法对激光雷达数据进行去噪处理，当存在椒盐噪声时，剔除椒盐噪声并进入步骤s4；否则直接进入步骤s4；
44.s4，采用双边滤波图像处理算法对激光雷达数据进行平滑，进而获得质控后的消光系数。
45.蒙特卡洛方法通过统计采样能够提供数理问题的近似解(fishman，1996)，是处理非线性反演问题和同时考虑多种不确定性因子的常用方法(唐晓，2010)。transfer中采用的蒙特卡洛分析法(monte carlo uncertaintyanalysis，mca)基于对激光雷达信号的随机采样从而区分高不确定性和低不确定性的消光系数。雷达信号的垂直分辨率为6m，通过气溶胶激光雷达系统中的双线性内插(yang et al.，2017a)，输出的消光系数垂直廓线的分辨率为30m。为了表征不确定性从输入到输出的传播，利用对应激光雷达产品垂直分辨率的5格滑动窗口来计算每个高度的信号均值和标准差。基于均值和标准差，在每个高度层建立高斯分布。这个建立高斯分布的方式同样用在elda中(mattis et al.，2016)。在每个高度层的高斯分布中随机抽取50个样本，从而构建了50组激光雷达信号垂直廓线。这50组雷达
信号连同原始雷达信号一同进入气溶胶激光雷达系统得到51组消光系数垂直廓线。通过其中的50组消光系数廓线可以得到每个高度层的消光系数均值和标准差，从而构建了消光系数相对不确定性垂直廓线。由原始信号反演而得到的消光系数垂直廓线通过以上构建的消光系数相对不确定性垂直廓线而进行数据质量筛选。考虑到大气污染主要发生在大气边界层内，于是在边界层上下分别设置不同的消光系数相对不确定性阈值。2006
‑
2017年华北平原的大气边界层高度约为1km(su et al.，2018)，同时为了减少基于激光雷达计算的边界层高度而引入新的误差，本发明中的大气边界层高度设为1km。大气边界层以下的相对不确定性阈值设为20％，这与elda中的设置相同(mattis et al.，2016)。考虑到利用fernald反演法的定标高度以下的积分误差相对较小以及本发明相对较低的边界层高度设置，边界层以上的阈值设为30％。在气溶胶激光雷达系统中，雷达比被假定为50sr，这是城市气溶胶和沙尘气溶胶的典型值(liu et al.，2002；shimizu，2004；shimizu et al.，2016)，同时也是在532nm不易受到相对湿度影响的值(takamura and sasano，1987)。但是雷达比取决于气溶胶空间的非均一性(kovalev，1995)且会随高度而变化。雷达比的不确定性被广泛研究，包括不同的气溶胶组分(ackerman and toon，1981；m
ü
ller et al.，2007)、不同地理位置(tesche et al.，2007)和不同高度(et al.，2012)。由于transfer的关注重点在于数据质控方案对激光雷达数据产品的改善程度，因此参考最近相关研究(xie et al.，2010)，雷达比的相对不确定性设为10％。
46.经过蒙特卡洛不确定性分析后，激光雷达数据中仍然会有一些高不确定性值存在。消光系数垂直廓线的时间序列可以认为是一种二维数据，一维的地面常规观测质控方法无法直接应用。因此，transfer采借用一种简单但是实用的图像去噪算法(denoising method，dn)来去除激光雷达数据中的脉冲噪声，即基于每个格点周围非0格点小于等于3则被去除的原则。
47.在气溶胶激光雷达系统中，每5min内的3000个回波信号被平均后得到每15min的雷达信号垂直廓线。为了延长使用寿命和减少维护成本，其间的10min内雷达的激光器停止工作(shimizu，2004)。尽管气溶胶激光雷达产品的时间分辨率为15min，但是气溶胶在激光器关闭的10min内的变化特征无法被捕捉。这会导致反演的消光系数垂直廓线原则上无法完全表征气溶胶的变化特征。因此，这种情况下需要考虑保持现有激光雷达观测的大气污染烟羽边缘，而对烟羽内部进行平滑处理。高斯低通滤波通过计算周围点的加权值来平滑数据，其中的权重值随距离衰减。高斯低通滤波仅仅考虑了距离上的相关性，而忽略了数值本身之间的相似性，会使得激光雷达数据中的烟羽边缘结构被破坏。因此，transfer采用能保持数据边缘结构的同时进行数据平滑的双边滤波图像处理算法(bilateral filtering，bf)来对激光雷达数据进行平滑(tomasi and manduchi，1998)。这种算法同时考虑了距离邻近度和值域相似性。bf对消光系数f(x)滤波计算公式如下：
[0048][0049][0050]
其中，c(ξ，x)表征邻近点ξ与中心点x的距离邻近度，s(f(ξ)，f(x))表示两点之间的消光系数值相似程度。距离邻近函数c(ξ，x)和值域相似性函数s(f(ξ)，f(x))计算公式如下：
[0051][0052][0053]
其中，d(ξ，x)表示两点之间的距离，δ(f(ξ)，f(x))表示两邻近点之间的消光系数值的差异。σ
d
和σ
r
分别表示距离和消光系数的标准差，transfer中分别设为0.2个格距和0.2km
‑1。另外，离散化积分直径设为5个格距。
[0054]
下面以nies激光雷达为例，简单介绍transfer对激光雷达数据处理的实际效果与改善程度。nies激光雷达的观测时间段为2017年1月、4月、7月和10月。基准数据采用aeronet产品。通过相较于angstrom插值具有更高准确性和可靠性的二次多项式插值(eck et al.，1999)，将aeronet产品中波长为440nm、500nm和670nm的aod转化成532nm的aod。将原始的和通过transfer质控后的消光系数在垂直上积分得到基于激光雷达的aod。通过transfer质控前后的消光系数转化为aod后与aeronet的aod产品的统计检验结果如表1所示。如表1所示，经过transfer质控，rmse(the root mean square error)由0.87降至0.51，相对下降了41％，同时mfe(the mean fractional error)相对下降了17％，而相关系数(pearson correlation coefficient，r)改善有限。由此可知，transfer通过mca、dn和bf集成方法改善了消光系数的绝对偏差而对原始数据的线性趋势影响不大。
[0055]
表1：2017年1月、4月、7月和10月transfer质控前后的消光系数转化的aod与aeronet的aod的rmse、mef和相关系数
[0056][0057]
注：bt表示transfer质控前，at表示transfer质控后，ai表示transfer质控前后的绝对增量.ri表示transfer质控前后的相对增量
[0058]
在证明了transfer质控方案对整个观测期间的消光系数改善有正贡献后，接下来对transfer方案内部的子方案之间进行两两比较，以定量评估不同子方案对消光数据的影响程度。图3为transfer各步骤质控前后的垂直平均消光系数(vertical average extinction coefficients，vaecs)散点图。图3中每行的横坐标分别为原始vaec、经过mca的vaec和经过dn的vaec，每列的纵坐标分别为经过mca的vaec、经过dn的vaec和经过整个transfer方案的vaec。从各个子方案单独对数据产生的影响角度来看(图3a、图3d和图3f)，对原始数据影响最大的子方案是mca，具有最大的rmse(0.08km
‑1)和最低的相关系数(0.91)；其次是dn，rmse为0.04km
‑1，相关系数为0.96；bf的影响最小。沿着图3a、图3b和图3c可以看出transfer中子方案的逐步实施相较于原始数据对数据的变化。如图所示，散点图
中的散点主要集中在坐标(0.05km
‑1，0.045km
‑1)附近，与斜率为1的直线接近，消光系数原始数据整体呈现出轻微高估的趋势。图中依然会存在若干低估的散点，这将会在后文中进行讨论。但是在坐标(0.55km
‑1，0.35km
‑1)和坐标(0.5km
‑1，0.1km
‑1)附近分别存在两团散点簇。第二处的散点簇在质控前具有相对较大的消光系数值，而在质控后具有相对较低的消光系数值，这主要是原始数据中这些数据点具有较大的不确定性导致。在本发明中将这个散点簇范围称为“虚高区”，同理，称第一处散点簇为“实高区”。
[0059]
将图3c中transfer质控前后的vaec与气象因子进行耦合来进一步探讨transfer在不同气象条件下的表现(图4a、图5a、图6a)。如图4a所示，大多数风矢量集中在坐标(0.05km
‑
1，0.045km
‑
1)附近，以强偏北风为主。风速相对较小和风向为偏南方向的风矢量主要在“实高区”，而“虚高区”的风矢量风速较低且无统一的风向。如图5a所示，高能见度主要发生在有强偏北风时。能见度小于5km的情况同时出现在“实高区”和“虚高区”。另外，“实高区”还包含一些能见度接近15km的个例。相对湿度(relative humidity，rh)低于40％主要出现在坐标(0.05km
‑1，0.0045km
‑1)附近。“实高区”既有rh大于80％的个例，也存在rh约等于60％的个例，而“虚高区”的rh都大于80％。通过耦合气象因子可知，vaecs的高估，特别是出现在“实高区”，主要发生在弱南风、低能见度、高rh等的条件下，而这些条件恰是华北平原典型不利气象条件的特征(wang et al.，2014)。
[0060]
气体和颗粒物的吸收和散射共同贡献了消光系数，其中气溶胶对总消光系数贡献显著(singh and dey，2012)。最近研究发现ho
x
和no
x
自由基的存在使得pm
2.5
和o3具有很强的负相关(li et al.，2019a，2019b)。因此，地面pm
2.5
和o3分别对消光产生直接和间接的影响。用与气象因子相同的方法，将这两种污染因子配置到vaec散点图上(图7a和图8a)。在“实高区”的pm
2.5
质量浓度大于250μg/m3，坐标(0.05km
‑1，0.045km
‑1)附近和“虚高区”的pm
2.5
浓度低于75μg/m3。o3与pm
2.5
具有相似的分布结构，这可能是由于当地面o3浓度较高时，对消光贡献较高的气溶胶污染烟羽主要位于高空，而当地面pm
2.5
浓度较高时，气溶胶污染团主要集中在地面。
[0061]
尽管结合气象因子和污染因子对观测期间transfer的总体表现进行了分析，但是不同季节的主导因子也各不相同。由表1可知，rmse从绝对增量来看，在夏季、冬季、春季和秋季呈现下降趋势，而从相对变化来看呈现完全相反的顺序。以上结果取决于aeronet的aod的偏离程度的，也许会根据选取不同校准数据而改变。而为了研究transfer在不同季节对消光系数数据本身的影响，rmse和相关系数被用来表征质控前后数据的改变程度，而不是与aeronet的aod的符合程度。如图4b所示，在冬季，风矢量存在显著差异：强偏北风出现在坐标(0.05km
‑1，0.045km
‑1)附近，弱风出现在实高区。冬季的能见度和rh分布规律相似(图5b和图6b)。在坐标(0.05km
‑1，0.045km
‑1)附近出现能见度大于10km，甚至达到30km，rh小于20％的情况，而在“实高区”出现能见度低于5km，rh高于80％的情况。在“实高区”pm
2.5
质量浓度达到250μg/m3以上，而在坐标(0.05km
‑1，0.045km
‑1)附近pm
2.5
质量浓度低于75μg/m3。冬季，o3质量浓度在各个散点簇区域都很低(图8b)。四季中，冬季vaecs的变化最小，具有最大的相关系数(0.93)和最小的rmse(0.07μg/m3)，而夏季vaecs的变化最大，具有最小的相关系数(0.8)和最大的rmse(0.15μg/m3)。在坐标(0.05km
‑1，0.045km
‑1)出现大于15km的能见度，而在“实高区”和“虚高区”都出现了小于5km的能见度。夏季rh的分布结构与能见度相似。但是，部分高能见度和低rh的散点出现在“实高区”，这主要是由于气溶胶不同的垂直分
布。夏季的o3散点分布与冬季pm
2.5
散点分布相似，“实高区”的臭氧质量浓度大于200μg/m3，而在坐标(0.05km
‑1，0.045km
‑1)附近的o3浓度低于50μg/m3。与春季和夏季的比较发现，图7a和图8a中“实高区”的大值分别为pm
2.5
在冬季的高浓度和o3在夏季的高浓度的贡献。
[0062]
如果vaecs能提供大气柱的丰富信息，那消光系数的水平平均能够给出transfer对气溶胶垂直结构的订正作用。如图9所示，经过transfer质控后的消光系数平均廓线为负，表明气溶胶激光雷达系统消光系数的系统性高估得到订正，这与前文消光系数垂直平均的结论一致。夏季transfer质控前后的增量绝对值在0.5～6km高度范围内为四个季节最大，表明transfer对夏季的消光系数改善程度最大。在2km以下，冬季的增量绝对值大于春季，而在2km以上则相反。除了春季，夏季、秋季和冬季的消光系数时间平均廓线在大气边界层内存在极大值，分别为0.09km
‑1、0.09km
‑1和0.05km
‑1，出现的高度分别为0.95km、0.6km和0.4km。差异出现的原因是大气边界层高度的变化受热力对流、机械湍流和气象条件的共同作用(guinot et al.，2006)。如图9a所示，夏季的消光系数平均垂直廓线在0.5～2km具有极大值，这是由于夏季发生的气溶胶区域输送(he et al.，2009)以及高层pm
2.5
具有较高吸湿性(sun et al.，2013)。
[0063]
尽管验证了transfer质控方案在长时间激光雷达观测的数据质量控制表现，但仍需在污染个例下探究transfer对激光雷达数据质量的改善程度。2017年7月1日至7月3日的污染个例被用来检验transfer的质控效果。结果表明，经过transfer质控后，该污染个例期间消光系数换算的aod与aeronet的aod的rmse由1.4降至0.84，mfe也从0.58降至0.43。
[0064]
图9为此次污染个例期间transfer质控方案中各个子方案的质控效果，图10与图3类似，为污染个例期间的vaec在各个质控子方案中的两两比较结果。如图所示，激光雷达成功捕获了这次污染过程中的多次消光系数大值。图9a显示在2～6km范围内，有很多脉冲噪声的存在。尽管这些脉冲噪声可以被人为去除，但是去除的情况会随着不同操作者而改变。同时有些脉冲噪声无法直接被分辨出来。图9b为原始数据经过mca质控后的结果，可见大量高不确定性值被去除，这些值有些能直接被分辨出来而有些不能，其中2～4km范围内的效果最显著。激光雷达探测的盲区也因为较高的相对不确定性而被一定程度去除。如图10a所示，经过mca质控之后的消光系数数据显示，大量vaec集中在坐标(0.1km
‑1，0.09km
‑1)附近，随后是“实高区”，仅有若干散点出现在“虚高区”。本发明中采用的dn算法是mca的一种补充，能在mca质控之后继续去除剩余的脉冲噪声。在4km以上，dn算法的作用更为显著。dn的补充作用也能在图10d中表现出来，经过mca质控之后的vaec和经过dn质控后的vaec的rmse和相关系数分别为0.05km
‑1和0.96，分别小于和大于原始数据与经过mca质控后数据的rmse和相关系数。bf是transfer质控方案的最后一步，对经过mca质控和dn质控后的消光系数数据进行“保边去噪”，即不破坏消光系数大值区边缘而对大值区内部进行平滑。由图10f可知，经过bf质控后，消光系数数据趋势并未发生改变，相关系数约为1.0。如图9e可知，transfer主要去除了位于气溶胶层上部的异常值，提高了气溶胶层的平滑度，维持了原有气溶胶垂直结构。原始数据的高估主要是由于一些参数设置中存在的不确定性(例如剔除云的阈值)。
[0065]
本发明的保护范围不限于本发明还公开了一种气溶胶激光雷达数据质控的方法，当一种控制气溶胶激光雷达数据质量的系统利用了上述的方法或者基本原理，也属于本发明的保护范围内。
[0066]
以上结合附图对本发明的实施方式作了详细说明，但本发明不限于所描述的实施方式。对于本领域的技术人员而言，在不脱离本发明原理和精神的情况下，对这些实施方式进行多种变化、修改、替换和变型，仍落入本发明的保护范围内。