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

技术领域

本发明涉及气溶胶激光雷达技术领域，特别涉及一种气溶胶激光雷达数据质控的方法及系统。

背景技术

气溶胶激光雷达(Light Detection And Ranging，Lidar)作为传统无线电或微波雷达(Radio Detection And Ranging，Radar)向光学频段的延伸，能够提供高时空分辨率垂直分布特性，是探测气溶胶垂直结构的有力工具。

如图1所示，气溶胶激光雷达通过激光器向大气发射激光，利用望远镜系统接收大气分子和大气气溶胶的后向散射信号，借助激光雷达方程进行回波信号的反演，来主动遥感大气气溶胶垂直结构。气溶胶激光雷达基于弹性散射Mie理论，雷达方程为：

其中，r表示当前回波信号所对应的高度，P(r)为r高度处的回波信号功率，P

根据激光雷达的部署位置，可以分为空基激光雷达、航测激光雷达和地基激光雷达。空基激光雷达主要有CALIPSO(Cloud-Aerosol Lidar and Infrared PathfinderSatellite Observation)和CATS(Cloud-Aerosol Transport System)。CALIPSO是由美国NASA(the National Aeronautics and Space Administration)的ESSP(Earth SystemScience Pathfinder)和法国CNES(Centre National d’Etudes Spatiales)共同合作开展，旨在填补观测上对气溶胶和云的全球分布和理化性质的空白(Winker，2003)。2018年9月，CALIPSO从“A-train”轨道下降约16.5km，和CloudSat卫星共同组成“C-train”。CALIPSO为16日循环的太阳同步轨道卫星。CALIOP(TheCloud-Aerosol Lidar with OrthogonalPolarization)是搭载在CALIPSO上的双波长(532nm和1064nm)弹性散射激光雷达，以天底视角直接测量总(气溶胶与空气分子)衰减后向散射的垂直廓线(Grigas et al.，2015)。CALIOP在地表的移动速度约为7km/s，脉冲重复频率约为20Hz，因此，在低层CALIOP的水平分辨率约为333m。此外，在低层垂直分辨率为30m。CATS激光雷达是大在在ISS(theInternational Space Station)上的双波长(532nm和1064nm)激光雷达(McGill et al.，2015)。CATS的高度约为415km，倾斜角为51°，以大约4天的重复周期观测热带和中纬度地区(Yorks et al.，2016)。目前，有很多地基激光雷达网已经建设完成并投入使用，例如AD-NET(Asian Dust Network)(Sugimoto et al.，2014)、EARLINET(the EuropeanAerosolResearch Lidar Network)(Sicard et al.，2015)、MPLNET(the Mirco Pulse LidarNetwork)(Welton，2002)和CORALNet(the Canadian Operational Research AerosolLidar 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 etal.，2020a，2020b)等大气污染研究中。

随着对气溶胶激光雷达数据定量化程序的要求越来越高，气溶胶激光雷达数据质量控制的研究也逐渐开展起来。项衍(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 LidarDataAnalyzer)(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方案，综合的软件层面的激光雷达数据质控研究相对较少。因此，亟需对气溶胶激光雷达软件层面的质控开展相关研究。

发明内容

本发明要解决的技术问题是提供一种气溶胶激光雷达数据质控的方法，以从综合的软件层面对激光雷达数据进行优化。

为了解决上述技术问题，本发明的技术方案为：

一种气溶胶激光雷达数据质控的方法，包括：

S1，利用气溶胶激光雷达系统获取探测大气气溶胶的激光雷达数据，当数值为0～1时，进入步骤S2；当数值＞1时，将数值设为1，然后进入步骤S3；当数值＜0时，将＜0的数值以及标记数据设为0，然后进入步骤S4；

S2，根据雷达比和雷达信号对激光雷达数据进行蒙特卡洛不确定性分析，当存在高不确定性值时，则对高不确定性值剔除，并进入步骤S4；当存在低不确定性值时，则直接进入步骤S3；

S3，采用图像去噪算法对激光雷达数据进行去噪处理，当存在椒盐噪声时，剔除椒盐噪声并进入步骤S4；否则直接进入步骤S4；

S4，采用双边滤波图像处理算法对激光雷达数据进行平滑，进而获得质控后的消光系数。

进一步地，步骤S1中所述利用气溶胶激光雷达系统获取探测大气气溶胶的激光雷达数据包括：

S11，计算每个高度的信号均值和标准差；

S12，基于均值和标准差，在每个高度层建立高斯分布；

S13，在每个高度层的高斯分布中随机抽取N个样本，从而构建N组激光雷达信号垂直廓线；

S14，将这N组雷达信号连同原始雷达信号一同进入气溶胶激光雷达系统得到N+1组消光系数垂直廓线；

S15，通过其中的N组消光系数廓线得到每个高度层的消光系数均值和标准值，从而构建消光系数相对不确定性垂直廓线。

进一步地，步骤S2中：大气边界层高度设为1km。

大气边界层以下的相对不确定性阈值设为20％，以上的不确定性阈值设有30％。

进一步地，步骤S2中：雷达比被假定为50sr。

进一步地，雷达比的相对不确定性设为10％。

进一步地，步骤S3中：每个格点周围非0格点小于等于3则被去除。

本发明还公开了一种控制气溶胶激光雷达数据质量的系统，使用以上任一项所述的方法。

与现有技术相比，本发明的有益效果是：

本发明包含蒙特卡洛不确定性分析、图像去噪算法以及双边滤波，其中蒙特卡洛不确定性分析基于对激光雷达信号的随机采样从而区分高不确定性和低不确定性的消光系数；图像去噪算法用于去除激光雷达数据中的脉冲噪声；双边滤波对激光雷达数据进行平滑，并且能够保持数据边缘结构；本发明主要去除了位于气溶胶上部的异常值，提高了气溶胶层的平滑度，维持了原有的气溶胶垂直结构。

附图说明

图1为气溶胶激光雷达硬件的结构图；

图2为本发明一种气溶胶激光雷达数据质控的方法的流程图；

图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的关系；

图4为本发明2017年1月、4月、7月和10月以及4个月耦合风矢量后TRANSFER质控前后的垂直平均消光系数比较的示意图；

图5为本发明2017年1月、4月、7月和10月以及4个月耦合能见度结果后TRANSFER质控前后的垂直平均消光系数比较的示意图；

图6为本发明2017年1月、4月、7月和10月以及4个月耦合RH后TRANSFER质控前后的垂直平均消光系数比较的示意图；

图7为本发明2017年1月、4月、7月和10月以及4个月耦合PM

图8为本发明2017年1月、4月、7月和10月以及4个月耦合O

图9为本发明2017年7月1日至3日NIES激光雷达系统反演的消光系数时间-高度截面图原始数据(a)，经过MCA质控后(b)、经过DN质控后(c)、经过BF质控后(d)和质控增量(e)的效果图；

图10为本发明2017年7月1日至3日的消光系数垂直平均(VAEC)表征的散点图。

具体实施方式

下面结合附图对本发明的具体实施方式作进一步说明。在此需要说明的是，对于这些实施方式的说明用于帮助理解本发明，但并不构成对本发明的限定。此外，下面所描述的本发明各个实施方式中所涉及的技术特征只要彼此之间未构成冲突就可以相互组合。

本发明公开了一种气溶胶激光雷达数据质控的方法，是一套专门针对气溶胶激光雷达数据质量控制的集成方案，也可以称为TRANSFER(a data quality assurance andcontrol scheme for aerosol lidar，TRANSFER)，如图2所示，主要包含蒙特卡洛不确定性分析和双边滤波，具体包括以下步骤：

S1，利用气溶胶激光雷达系统获取探测大气气溶胶的激光雷达数据，当数值为0～1时，进入步骤S2；当数值＞1时，将数值设为1，然后进入步骤S3；当数值＜0时，将＜0的数值以及标记数据设为0，然后进入步骤S4；

S2，根据雷达比和雷达信号对激光雷达数据进行蒙特卡洛不确定性分析，当存在高不确定性值时，则对高不确定性值剔除，并进入步骤S4；当存在低不确定性值时，则直接进入步骤S3；

S3，采用图像去噪算法对激光雷达数据进行去噪处理，当存在椒盐噪声时，剔除椒盐噪声并进入步骤S4；否则直接进入步骤S4；

S4，采用双边滤波图像处理算法对激光雷达数据进行平滑，进而获得质控后的消光系数。

蒙特卡洛方法通过统计采样能够提供数理问题的近似解(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)和不同高度(

经过蒙特卡洛不确定性分析后，激光雷达数据中仍然会有一些高不确定性值存在。消光系数垂直廓线的时间序列可以认为是一种二维数据，一维的地面常规观测质控方法无法直接应用。因此，TRANSFER采借用一种简单但是实用的图像去噪算法(DeNoisingmethod，DN)来去除激光雷达数据中的脉冲噪声，即基于每个格点周围非0格点小于等于3则被去除的原则。

在气溶胶激光雷达系统中，每5min内的3000个回波信号被平均后得到每15min的雷达信号垂直廓线。为了延长使用寿命和减少维护成本，其间的10min内雷达的激光器停止工作(Shimizu，2004)。尽管气溶胶激光雷达产品的时间分辨率为15min，但是气溶胶在激光器关闭的10min内的变化特征无法被捕捉。这会导致反演的消光系数垂直廓线原则上无法完全表征气溶胶的变化特征。因此，这种情况下需要考虑保持现有激光雷达观测的大气污染烟羽边缘，而对烟羽内部进行平滑处理。高斯低通滤波通过计算周围点的加权值来平滑数据，其中的权重值随距离衰减。高斯低通滤波仅仅考虑了距离上的相关性，而忽略了数值本身之间的相似性，会使得激光雷达数据中的烟羽边缘结构被破坏。因此，TRANSFER采用能保持数据边缘结构的同时进行数据平滑的双边滤波图像处理算法(Bilateral Filtering，BF)来对激光雷达数据进行平滑(Tomasi and Manduchi，1998)。这种算法同时考虑了距离邻近度和值域相似性。BF对消光系数f(x)滤波计算公式如下：

其中，c(ξ，x)表征邻近点ξ与中心点x的距离邻近度，s(f(ξ)，f(x))表示两点之间的消光系数值相似程度。距离邻近函数c(ξ，x)和值域相似性函数s(f(ξ)，f(x))计算公式如下：

其中，d(ξ，x)表示两点之间的距离，δ(f(ξ)，f(x))表示两邻近点之间的消光系数值的差异。σ

下面以NIES激光雷达为例，简单介绍TRANSFER对激光雷达数据处理的实际效果与改善程度。NIES激光雷达的观测时间段为2017年1月、4月、7月和10月。基准数据采用AERONET产品。通过相较于Angstrom插值具有更高准确性和可靠性的二次多项式插值(Ecket 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集成方法改善了消光系数的绝对偏差而对原始数据的线性趋势影响不大。

表1：2017年1月、4月、7月和10月TRANSFER质控前后的消光系数转化的AOD与AERONET的AOD的RMSE、MEF和相关系数

注：BT表示TRANSFER质控前，AT表示TRANSFER质控后，AI表示TRANSFER质控前后的绝对增量.RI表示TRANSFER质控前后的相对增量

在证明了TRANSFER质控方案对整个观测期间的消光系数改善有正贡献后，接下来对TRANSFER方案内部的子方案之间进行两两比较，以定量评估不同子方案对消光数据的影响程度。图3为TRANSFER各步骤质控前后的垂直平均消光系数(Vertical AverageExtinction Coefficients，VAECs)散点图。图3中每行的横坐标分别为原始VAEC、经过MCA的VAEC和经过DN的VAEC，每列的纵坐标分别为经过MCA的VAEC、经过DN的VAEC和经过整个TRANSFER方案的VAEC。从各个子方案单独对数据产生的影响角度来看(图3a、图3d和图3f)，对原始数据影响最大的子方案是MCA，具有最大的RMSE(0.08km

将图3c中TRANSFER质控前后的VAEC与气象因子进行耦合来进一步探讨TRANSFER在不同气象条件下的表现(图4a、图5a、图6a)。如图4a所示，大多数风矢量集中在坐标(0.05km-1，0.045km-1)附近，以强偏北风为主。风速相对较小和风向为偏南方向的风矢量主要在“实高区”，而“虚高区”的风矢量风速较低且无统一的风向。如图5a所示，高能见度主要发生在有强偏北风时。能见度小于5km的情况同时出现在“实高区”和“虚高区”。另外，“实高区”还包含一些能见度接近15km的个例。相对湿度(Relative Humidity，RH)低于40％主要出现在坐标(0.05km

气体和颗粒物的吸收和散射共同贡献了消光系数，其中气溶胶对总消光系数贡献显著(Singh and Dey，2012)。最近研究发现HO

尽管结合气象因子和污染因子对观测期间TRANSFER的总体表现进行了分析，但是不同季节的主导因子也各不相同。由表1可知，RMSE从绝对增量来看，在夏季、冬季、春季和秋季呈现下降趋势，而从相对变化来看呈现完全相反的顺序。以上结果取决于AERONET的AOD的偏离程度的，也许会根据选取不同校准数据而改变。而为了研究TRANSFER在不同季节对消光系数数据本身的影响，RMSE和相关系数被用来表征质控前后数据的改变程度，而不是与AERONET的AOD的符合程度。如图4b所示，在冬季，风矢量存在显著差异：强偏北风出现在坐标(0.05km

如果VAECs能提供大气柱的丰富信息，那消光系数的水平平均能够给出TRANSFER对气溶胶垂直结构的订正作用。如图9所示，经过TRANSFER质控后的消光系数平均廓线为负，表明气溶胶激光雷达系统消光系数的系统性高估得到订正，这与前文消光系数垂直平均的结论一致。夏季TRANSFER质控前后的增量绝对值在0.5～6km高度范围内为四个季节最大，表明TRANSFER对夏季的消光系数改善程度最大。在2km以下，冬季的增量绝对值大于春季，而在2km以上则相反。除了春季，夏季、秋季和冬季的消光系数时间平均廓线在大气边界层内存在极大值，分别为0.09km

尽管验证了TRANSFER质控方案在长时间激光雷达观测的数据质量控制表现，但仍需在污染个例下探究TRANSFER对激光雷达数据质量的改善程度。2017年7月1日至7月3日的污染个例被用来检验TRANSFER的质控效果。结果表明，经过TRANSFER质控后，该污染个例期间消光系数换算的AOD与AERONET的AOD的RMSE由1.4降至0.84，MFE也从0.58降至0.43。

图9为此次污染个例期间TRANSFER质控方案中各个子方案的质控效果，图10与图3类似，为污染个例期间的VAEC在各个质控子方案中的两两比较结果。如图所示，激光雷达成功捕获了这次污染过程中的多次消光系数大值。图9a显示在2～6km范围内，有很多脉冲噪声的存在。尽管这些脉冲噪声可以被人为去除，但是去除的情况会随着不同操作者而改变。同时有些脉冲噪声无法直接被分辨出来。图9b为原始数据经过MCA质控后的结果，可见大量高不确定性值被去除，这些值有些能直接被分辨出来而有些不能，其中2～4km范围内的效果最显著。激光雷达探测的盲区也因为较高的相对不确定性而被一定程度去除。如图10a所示，经过MCA质控之后的消光系数数据显示，大量VAEC集中在坐标(0.1km

本发明的保护范围不限于本发明还公开了一种气溶胶激光雷达数据质控的方法，当一种控制气溶胶激光雷达数据质量的系统利用了上述的方法或者基本原理，也属于本发明的保护范围内。

以上结合附图对本发明的实施方式作了详细说明，但本发明不限于所描述的实施方式。对于本领域的技术人员而言，在不脱离本发明原理和精神的情况下，对这些实施方式进行多种变化、修改、替换和变型，仍落入本发明的保护范围内。