气溶胶激光雷达数据融合方法及系统

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

背景技术：

资料同化根据不确定性将外场观测与数值模式进行结合，得到更加准确的状态变量以及相应的不确定性(lahozetal.，2010)，是基于变分和卡尔曼滤波等数据融合方法的统称。中国的气溶胶激光雷达的数据融合研究刚刚起步。郑海涛(2018)开发了基于模式第二层pm2.5质量浓度与消光系数比例的观测算子，首次构建了以oi(optimalinterpolation)为融合算法，基于naqpms(thenestedairqualitypredictmodelsystem)的地基激光雷达数据融合系统，在2014年11月10日-11月20日，同化了12台地基激光雷达联合观测。结果显示，加入激光雷达观测后的同化结果改进效果显著(rmse改进约为24％)，对模式影响的持续时间也更长。项衍(2018)利用gsi(gridpointstatisticalinterpolation)的3d-var算法，通过气溶胶质量浓度和消光系数的拟合关系，基于wrf-chem(weatherresearchandforecastingmodelwiththechemistry)融合了2017年11月29日-12月5日的12台地基激光雷达和160分钟的共轭走航激光雷达观测。消光系数拟合到质量浓度的相关性平均值达到了0.86，表明了构建的观测算子的可行性。同化激光雷达使得相关性平均提高42.9％，rmse平均减少50.9％。另外，该研究还厘清了这次重污染过程的前期、污染时和污染后期的pm2.5三维质量浓度场的演变过程。chengetal.(2019a)将辐射传输模式crtm(thecommunityradiativetransfermodel)作为观测算子，融合wrf-chem和2018年3月13日12:00-18:00的位于北京的四台地基激光雷达观测，得到改善后的气溶胶初始场。结果表明，同化激光雷达消光系数垂直廓线能显著改善pm2.5预报低估情况，rmse下降25μg/m³，相关系数增加0.45。chengetal.(2019b)利用4-dletkf(afour-dimensionallocalensembletransformkalmanfilter)算法融合大气化学传输模式nicam-sprintars(non-hydrostaticicosahedralatmosphericmodel-spectralradiationtransportmodelforaerosolspecies)和2016年11月的caliop的消光系数垂直廓线以及modis的aod产品。结果显示，同化caliop和modis都可以改善模式对气溶胶的模拟。对于订正不正确的气溶胶垂直分布，同化caliop的效果要优于modis，而同化modis对模式aot的模拟改善更优。由于caliop相对稀疏的观测时空范围，同化caliop对模式aot模拟影响不大。liangetal.(2020)发展了一套针对气溶胶消光系数观测算子和气溶胶质量浓度算子的3d-var数据同化系统。该研究将同化系统应用于融合wrf-chem和2018年11月13日的5台地基激光雷达来得到更优的pm2.5预报初值场。pm2.5预报初值场误差降低了10.5μg/m³(17.6％)，预报时效能持续24小时以上。maetal.(2020)利用wrf-chem/dart(dataassimilationresearchtestbed)大气化学资料同化系统，融合地面pm2.5观测、aod和地基激光雷达消光系数研究了2018年6月和11月的两个持续7天的污染个例，研究发现同化aod和气溶胶消光系数垂直廓线对融合结果会出现负的增量。这主要是由于wrf-chem对边界层内的湿度模拟偏低，导致出现模拟消光系数偏低，而模拟的地面pm2.5偏高的两种偏差趋势不一致的情况。该研究通过偏差订正，减轻了上述偏差趋势不一致的问题，使同时同化地面pm2.5、aod和消光系数垂直廓线的结果最优。

激光雷达数据融合有主要的三个不确定性来源：大气化学传输模式、激光雷达观测和融合算法。大气化学传输模式的不确定性来源主要有气象场、气溶胶初边值、排放源、模式参数化方案和观测算子等(maetal.，2019)。作为连接模式变量和观测变量的“桥梁”，模式观测算子是影响气溶胶激光雷达数据融合的重要因素。目前的研究中多采用单一的方法作为观测算子，例如假设球形粒子下的米散射(yumimotoetal.，2007)、比值法(郑海涛，2018)、基于观测的improve(theinteragencymonitoringofprotectedvisualenvironments)方法(liangetal.，2020)、拟合法(项衍，2018)和简化的辐射传输模式法(chengetal.，2019a)等。但由于气溶胶在不同的粒径谱、不同的组分、不同的复折射指数、不同环境湿度下的消光表现不同，有必要考虑集合的方式应对不同的气溶胶，例如针对球形粒子采用mie理论，针对椭球体粒子采用tmm(t-matrixmethod)或igom(improvedgeometricopticsmethods)方法，针对不规则粒子应用dda(discretedipoleapproximation)方法(gasteigerandwiegner，2018)。融合算法主要分为变分方法和集合卡尔曼滤波方法。集合卡尔曼滤波能处理模式模拟观测需要用到复杂的非线性观测算子，还能体现模式背景场的“流依赖”特性，更适合用于激光雷达数据融合研究(houtekamerandzhang，2016)。

技术实现要素：

本发明要解决的技术问题是提供一种气溶胶激光雷达数据融合方法及系统，以获得气溶胶激光雷达数据更加准确的状态变量以及相应的不确定性。

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

一种气溶胶激光雷达数据融合方法，包括以下步骤：

s1，n组模式初值分别作为嵌套网格空气质量预报模式子系统的输入值，而后得到n组预报场结果；

s2，n组预报场结果与雷达观测数据作为平行数据同化子系统的输入值，而后得到n组分析场结果；

s3，n组分析场结果作为下一个时刻的n组模式初值进行下一个时刻的预报与数据融合。

优选的，所述嵌套网格空气质量预报模式子系统包括源扰动分析处理模块、平流输送分析处理模块、湍流扩散分析处理模块、液相和非均相分析处理模块、重力沉降和干沉降分析处理模块、气相化学分析处理模块和湿清除分析处理模块中的一个或者任意多个；n组模式初值分别经过所述嵌套网格空气质量预报模式子系统中的一个或者任意多个模块后，转换成n组预报场结果。

优选的，所述平行数据同化子系统包括集合卡尔曼滤波模块、水平垂直局地花模块、垂直深度开发模块、remove观测算子模块、mpi二级并行模块和在线耦合模块中的一个或者任意多个；n组预报场结果各自与雷达数据经过所述平行数据同化子系统中的一个或者任意多个模块后，转换成n组分析场结果。

优选的，雷达观测数据的获得包括以下步骤：

s100，扫描气溶胶分布后获得地基激光雷达初始数据；

s200，对地基激光雷达初始数据极性数据质控管理后获得地基激光雷达质控数据；

s300，对地基激光雷达质控数据分别进行归一化处理和虚拟观测；

s400，对归一化后的数据以及虚拟观测后的数据进行不确定性分析，从而获得步骤s2中的雷达观测数据。

优选的，所述平流输送分析处理模块采用高精度正定质量守恒差分格式方案计算。

优选的，所述湍流扩散分析处理模块采用梯度输送理论。

优选的，所述重力沉降和干沉降分析处理模块中干沉降过程基于阻力模型处理，气体部分采用wesely方案，而气溶胶部分采用slinn方案。

优选的，所述气相化学分析处理模块中气相化学机制为cbm-z。

本发明还公开了一种气溶胶激光雷达数据融合系统，包括嵌套网格空气质量预报模式子系统和平行数据同化子系统，用于以上任一项所述的气溶胶激光雷达数据融合方法。

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

本发明通过嵌套网格空气质量预报模式子系统和平行数据同化子系统获得了更加准确的气溶胶激光雷达数据；嵌套网格空气质量预报模式子系统采用包括源扰动、平流输送、湍流扩散、液相和非均相、重力沉降和干沉降、气相化学和湿清除在内的方法来订正气溶胶激光雷达数据；平行数据同化子系统采用包括集合卡尔曼滤波、水平垂直局地化、垂直维深度开发、improve观测算子、mpi二级并行以及在线耦合嵌套网格空气质量预报模式子系统输出数据在内的方法来进一步提高气溶胶激光雷达数据的准确度，比如：气溶胶激光雷达数据融合能显著改善消光系数的模拟，以及在特定条件下的模式高估和模式低估都可以通过激光雷达数据融合得到订正。

附图说明

图1为本发明naqpms-pdaf数据融合系统运行流程图；

图2为本发明naqpms-pdaf数据融合系统运行流程的详解图；

图3为本发明naqpms-pdaf数据融合系统的应用图；

图4为本发明naqpms和pdaf的二级并行示意图；

图5为本发明300m处的先验rmse和先验总传播时间序列图；

图6为本发明470m处的先验rmse和先验总传播时间序列图；

图7为本发明680m处的先验rmse和先验总传播时间序列图；

图8为本发明930m处的先验rmse和先验总传播时间序列图；

图9为本发明激光雷达数据融合前后由消光系数换算的aod与aeronet的beijingpku站aod验证结果对比图；

图10为本发明激光雷达数据融合前后由消光系数换算的aod与aeronet的xianghe站aod验证结果对比图；

图11为本发明2019年4月12日激光雷达数据融合前后消光系数与caliop验证结果对比图；

图12为本发明2019年4月13日激光雷达数据融合前后消光系数与caliop验证结果对比图。

具体实施方式

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

实施例一

如图1所示，一种气溶胶激光雷达数据融合方法，包括以下步骤：sl，n组模式初值分别作为嵌套网格空气质量预报模式子系统的输入值，而后得到n组预报场结果，从而能得到当前时刻下的集合预报结果；s2，n组预报场结果与雷达观测数据作为平行数据同化子系统的输入值，而后得到n组分析场结果，从而得到了当前时刻下的集合分析场；s3，n组分析场结果作为下一个时刻的n组模式初值进行下一个时刻的预报与数据融合。

本实施例中个，第1组模式初值作为第1嵌套网格空气质量预报模式子系统的输入值，而后得到第1组预报场结果；第1组预报场结果和雷达观测数据作为平行数据同化子系统的输入值，而后得到第1组分析场结果；第2组模式初值至第n组模式初值同理，可以分别得到第2组至第n组分析场结果。

如图2所示，嵌套网格空气质量预报模式子系统包括源扰动分析处理模块、平流输送分析处理模块、湍流扩散分析处理模块、液相和非均相分析处理模块、重力沉降和干沉降分析处理模块、气相化学分析处理模块和湿清除分析处理模块中的一个或者任意多个；n组模式初值分别经过嵌套网格空气质量预报模式子系统中的一个或者任意多个模块后，转换成n组预报场结果。

本实施例中，嵌套网格空气质量预报模式子系统是一种双向嵌套的三维欧拉大气化学传输模式，简称为naqpms(nestedairquailitypredictionmodelsystem)。naqpms包含污染物排放、平流输送、湍流扩散、干湿沉降、气相、液相及非均相反应等物理与化学过程。其采用地形追随坐标作为垂直坐标，水平结构为多重网格嵌套，能同时模拟区域和城市尺度沙尘、pm2.5、pm10、so2、nox、co、o3、nh3等多种污染物。

化学成分欧拉输送方程的一般形式为：

在高度地形追随球坐标系下的三维化学成分欧拉输送方程为：

其中ci为第i种污染物浓度，t为时间，θ、分别为纬度和经度，r为地球半径，kθ、和kσ分别为经向、纬向、垂直湍流扩散系数，u、v为水平风速、p为化学转化项，s为源排放源速率，rd为干沉降项，wash为湿清除项，地形对污染物输送有一定的影响，σ为地形坐标，满足：

其中为模式顶高度，为地形高度，z为格点的海拔高度。在地形追随坐标系下，低层水平面是一个地形面，对于讨论污染物的分布是非常方便的。w为等效垂直运动速度，由空气质量守恒方程得出：

naqpms中的平流输送过程采用高精度正定质量守恒差分格式方案计算。该方案是基于通量法的一个修正和简化版，可保证化学物种的浓度为正和质量守恒，减少计算伪扩散等问题，该方案在对流层化学模式中得到广泛应用。湍流扩散采用梯度输送理论，即将大气扩散问题处理成在一定边界条件下求解扩散方程的问题。在naqpms中，由于水平扩散远小于水平平流的贡献，水平湍流扩散系数取常数。垂直扩散采用byunanddennis方案计算，该方案在cmaq和raqm等对流层化学模式中得到了广泛的应用。naqpms中的干沉降过程基于阻力模型处理，气体部分采用wesely方案，而气溶胶部分采用slinn等方案。干沉降通量f＝vd·c，其中vd为干沉降速率，c为气溶胶质量浓度。在阻力模型中，气体干沉降速率vd的计算如下：

其中，ra表示空气动力学，表征由于湍流扩散导致污染物向地面的迁移；rb代表次边界层阻力，表征由于分子扩散作用穿过薄层空气直接到达地面；rc为表面阻力，主要表征建筑物、植被等对污染物的捕获、吸附。对于气溶胶颗粒物，还会考虑重力沉降作用。另外，通常假定粒子可粘附到表面，即表面阻力rc＝0，其计算公式如下：

其中vs为重力沉降速率，主要与粒子的大小和密度相关。naqpms中湿清除模块主要包含云下(1600m以下)清除过程，该方案已成功用于东亚nox的湿清除模拟。云下雨水的冲刷过程定义为：

wash＝wa·ci

其中，wash为冲刷量，ci为污染物浓度，wa为冲刷系数，与雨强、雨滴谱分布及各种污染物的浓度、气溶胶粒子谱和气溶胶化学成分有关。naqpms中，云对流、液相化学、云内云下湿清除过程基于radm2机制处理，该机制已广泛应用于大气化学的数值模拟。利用wrf模式提供的云水、雨水资料和网格内的气象要素，考虑污染物的云雨吸收、溶解、电离、化学反应等液相过程，降水条件下的云内清除和云下冲刷，另外诊断积云对流是否存在，对污染物进行垂直再分布。naqpms的气相化学机制为cbm-z。cbm-z是基于cbm-iv发展起来的按结构分类的一个新的归纳化学机理。其化学反应方程由81个增加至176个，反应物种由32个增加至67个，其中涉及光化学反应的物种有20个。与cbm-iv机制相比，cbm-z在活性长寿命物种及中间产物的化学反应、无机物的化学反应、活性烷烃、石蜡以及芳香烃的化学反应、若干自由基以及异戊二烯的化学反应等方面的考虑更为全面。同时，该机制设定了包括背景条件、城市、远郊和生物区以及海洋等四个反应场景，可适用于全球、区域和城市尺度的研究。根据反应场景的不同，参与反应的物种和反应方程也不相同。在背景大气下，化学反应有74个；在城市大气中，化学反应为118个；在有生物源的城市大气中，化学反应增加为134个；在远离人为源的海洋上空，化学反应为109个；在近海，化学反应增加为153个；在受到生物源影响的近海地区，化学反应达到169个。naqpms中考虑的气溶胶过程主要包括无机气溶胶热力学过程、二次有机气溶胶生成、气溶胶表面非均相化学反应等。naqpms中采用气溶胶热力学模式isorropia来计算无机气溶胶的组分等于前体物的热力平衡与分配。气溶胶表面非均相化学反应模块，主要侧重对硫氧化物和氮氧化物的影响。在一定温度和湿度下，水汽在气溶胶粒子表面凝结增长形成薄水层，气态污染物向气溶胶表面扩散可被其表面溶液吸收、溶解。

继续参考图2，平行数据同化子系统包括集合卡尔曼滤波模块、水平垂直局地花模块、垂直深度开发模块、remove观测算子模块、mpi二级并行模块和在线耦合模块中的一个或者任意多个；n组预报场结果各自与雷达数据经过平行数据同化子系统中的一个或者任意多个模块后，转换成n组分析场结果。

为了实现模式数据与地基观测及空基观测的融合，利用pdaf(paralleldataassimilationframework)下的lestkf(localerrorsubspaceensembletransformkalmanfilter)方法来进行数据融合，搭建了naqpms-pdaf在线同化系统。pdaf是一款致力于简化部署集合同化的免费开源软件，为实施并行集合滤波和平滑提供了统一的框架，例如letkf、netf和estkf等。pdaf基于fortran编写，同时可通过mpi和openmp来并行运算。

estkf作为一种enkf，利用一组ne维的状态变量的集合i＝1，...，ne来表征模式系统的状态变量x^k和误差协方差矩阵pk。状态变量的维度nx包含模式变量和参数。其由集合平均来表征：

误差协方差矩阵为：

其中，为集合扰动矩阵。

集合同化的算法分为预报步和分析步，分别用下标f和a表示。预报步将模式系统的状态变量和误差协方差矩阵从t＝tk-1传播到下一个观测时间t＝tk。在模式运行期间，数值模式(mk)对每一个集合成员进行状态变量和误差协方差的传播：

在t＝tk，大小为ny的观测向量通过与真值系统x^f建立联系。其中，假设观测误差ε^k服从协方差矩阵r^k的正态分布，h是将模式空间投影至观测空间的观测算子。

分析步通过融合预报场和观测信息来提供新状态的估计。由于分析步是在同一个时间下进行，因此下面的介绍略去时间下标。estkf利用卡尔曼滤波更新方程，从预报场、观测场和这两者的不确定性信息，来计算分析场状态变量和误差协方差矩阵：

其中，k为卡尔曼增益：

k＝p^fh^t(hp^fh^t+r)^-1

对于集合：

y^f＝hx^f

其中，为预报场的集合矩阵，对应于模式场投影到观测空间的观测集合矩阵。

由于状态变量维度较大，以大小为ne的整体近似的样本协方差矩阵只是真实协方差矩阵的低阶近似，其秩最多为ne-1。estkf利用这一性质，将分析场写入集合张成的ne-1维的子空间中，这个空间被称为误差子空间。于是预报误差协方差矩阵p^f被写为：

其中，为：

l＝x^fω

矩阵定义为：

将集合矩阵x^f投影至误差子空间。这之后，p^a可以写为：

p^a＝lal^t

其中，转化矩阵为

a^-1＝ρ(ne-1)i+(y^fω)^tr^-1yω

其中，ρ∈[0，1]为遗忘因子。最后，更新步将预报集合转换为分析集合：

其中，为预报集合平均矩阵，的每一列为：

w为：

其中，c是通过奇异值分解获得的a的对称平方根。最后，分析集合x^a作为下一次集合预报的初值，滚动地进行下去直至同化窗口结束。

如图3所示，雷达观测数据的获得包括以下步骤：s100，扫描气溶胶分布后获得地基激光雷达初始数据；s200，对地基激光雷达初始数据极性数据质控管理后获得地基激光雷达质控数据；s300，对地基激光雷达质控数据分别进行归一化处理和虚拟观测；s400，对归一化后的数据以及虚拟观测后的数据进行不确定性分析，从而获得步骤s2中的雷达观测数据。

本实施例中，雷达观测数据在输入pdaf之前可以先进行处理，处理方法包括对地基激光雷达初始数据进行质控管理、归一化、虚拟观测以及不确定性分析。本发明气溶胶激光雷达数据融合系统在经过pdaf和naqpms在线耦合后，输出集合分析场，这些更加准确的气溶胶激光雷达数据可以应用到以下技术领域：数据融合系统不确定性分析、分析长时间pm2.5垂直结构变化特征和评估naqpms对pm2.5模拟/预报改善。本发明可以根据融合应用的效果来评估对地基激光雷达初始数据的处理是否妥当，是否需要调整。同时，本发明也可以根据融合应用的效果来对naqpms的方案进行改进。

如图4所示，naqpms模式采用“一级并行”方案，即设置一种通信域world来实现mpi分布式计算。如图1所示，为了在线耦合naqpms和pdaf，将通信域world拆分成三个子通信域：model、filter和couple，分别起模式集合模拟、集合卡尔曼滤波和数据交换的作用。此外，为了融合垂直观测，扩展了pdaf中原有的二维数据框架，增加了垂直维度，并增加了观测资料的垂直局地化半径参数来方便后续研究。

基于京津冀7台微脉冲激光雷达于2019年4月10日07:00-4月14日10：00的观测，对naqpms-pdaf数据融合系统进行测试。7台激光雷达分别布设于北京、保定(2台)、石家庄、沧州、衡水和邢台。在此期间有共有四条过境的caliop观测轨道。本研究时空范围内包含了2个有效的aeronet站点：beijing_pku站和xianghe站，还包含了72个pm2.5国控监测站。本研究的融合算法采用lestkf，局地化半径设为40格点(约为200km)，局地化权重函数为5-th多项式。通过扰动初始排放源来产生20组集合，so2和voc以及其余物种的排放源不确定性分别设为20％、50％和30％(maetal.，2019；zhangetal.，2006)。气溶胶激光雷达消光系数的观测误差设为0.02km^-1。

集合卡尔曼滤波类算法中，集合被用来表征模式背景场的演变特征，常用先验rmse和先验总传播来检验产生的集合传播是否充分。先验rmse的公式为先验总传播为(houtekameretal.，2009)。其中，表示第i个观测的观测值，表示观测空间先验集合平均，表示观测误差方法，表示观测空间先验观测误差方差。如图5-8所示，300m、470m、680m和930m高度处的先验rmse和先验总传播在量级和变化趋势上基本平衡，表明集合模拟传播较为充分。

根据京津冀7台激光雷达数据融合的分析场和预报场结果得知，经过数据融合后，7台激光雷达消光系数的rmse由0.25km^-1降至0.11km^-1，bias由0.13km^-1降至0.05km^-1，mae由0.17km^-1降至0.09km^-1，而相关系数由于0.28上升至0.81，bias也更向0附近集中。同时，消光系数1小时预报场的rmse由0.25km^-1降至0.18km^-1，bias由0.13km^-1降至0.09km^-1，mae(meanabsoluteerror)由0.17km^-1降至0.14km^-1，相关系数由0.28变为0.41。由此可见，分析场结果优于1小时预报场结果。经过内部验证表明，气溶胶激光雷达数据融合能显著改善消光系数的模拟。图9和图10分别为激光雷达数据融合前后由消光系数换算的aod与aeronet的beijing_pku站与xianghe站的aod验证结果。如图9所示，经过激光雷达数据融合后，尽管4月12日18：00出现偏离观测的情况，但总体上，aod时间序列与观测更为接近，rmse减少0.14，相关系数增加0.1。图10显示激光雷达数据融合后，基于消光系数换算的aod与aeronet的aod变化不大(rmse减少0.02)，这是由于naqpms模式相较于beijing_pku站点，在这一站点的aod模拟与观测更为接近。图11和图12分别是4月12日和4月13日激光雷达数据融合前后与caliop结果的验证图。如图所示，4月12日1000-2000m和4月13日1000m左右的模式高估，以及4月13日2000m以上出现的模式低估都通过激光雷达数据融合得到了订正。

实施例二

一种气溶胶激光雷达数据融合系统，包括嵌套网格空气质量预报模式子系统和平行数据同化子系统，用于实施实施例一中任一项所述的气溶胶激光雷达数据融合方法。

应当理解的是，任何系统或者设备，只要是利用本发明的基本原理，用来实施本发明上述任一项所述的气溶胶激光雷达数据融合方法，均可以认为是在本发明的保护范围之内。

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