反演同化技术在激光雷达中的应用研究

杜嵩山，赵文灿，王丽琼，左瑞亭，沈 沉

(1.江苏省环境监测中心，南京 210000；2.南京雨后地软环境技术有限公司，南京 210013；3.浙江省德清县气象局，浙江 德清 313200)

前 言

长江三角洲地区作为中国经济发展中心和城市化水平较高地区，灰霾(气溶胶细粒子浓度较高导致能见度恶化的天气现象)日数呈现逐年上升的趋势，已然成为我国主要的灰霾污染区[1-2]。

当前对颗粒物浓度预报，主要依据空气质量模式的数值预报，由于排放源、气象场以及模式本身的不确定性等因素，空气质量模式的预报结果与实况数据有较大差距。通过资料同化，耦合空气质量数值模式的模拟结果和大气观测资料，可以有效地改善数值模式的初始场，提高预报准确率，通常使用地面观测颗粒物浓度数据或卫星遥感结果进行同化[3-4]。刘汉卫等[5]利用WRF-Chem模式，对美国西部地区一次PM2.5污染过程的PM2.5总量及其主要成分进行同化实验，实验结果与实况接近，站点分布和演变趋势吻合较好；蒋维东等[6]利用卫星遥感数据，同化模式初始场AOD(Aerosol Optical Depth，气溶胶光学厚度)，更新后气溶胶分布与卫星观测值更加接近，表明同化AOD观测资料对初始场起到了较好地调整作用；靳璐滨等[7]构建本地化背景场误差协方差矩阵，利用监测站PM2.5和PM10逐小时观测数据，对南京青奥会期间进行滚动同化和预报实验，发现同化后，预报准确率明显提升，改进效果可以持续到20h以后。

作为一种新型监测仪器，激光雷达能展示出大气垂直结构的污染物分布和传播情况，多用于气溶胶浓度和颗粒物浓度监测[8]，在模式方面应用较少。Wang et al[9]首次基于单波长(355nm)激光雷达和POLAIR 3D模式开发了适用于该模式的同化方案，通过对模拟的雷达信号修正距离(PR2)直接进行同化以改进PM10和PM2.5的预报，同化后能有效地改进PM2.5预报结果；郑海涛[10]利用不同高度层气溶胶消光系数与150m处气溶胶消光系数比值关系作为PM2.5浓度垂直分布关系，以此订正原始场PM2.5垂直分布，研究过程带有一定主观性。

本文基于前人研究，构建一套操作简单、移植性强的激光雷达AOD同化系统：利用江苏省站等四个激光雷达观测站数据，针对2016年12月6日至10日一次污染过程进行逐小时模拟，通过反演同化方案建立AOD与PM2.5浓度关系，实现AOD三维变分同化。实验对比同化前后模式分析场、预报场的改变情况，定量评估模式预报效果，验证激光雷达AOD同化可行性，为后人研究提供参考和依据。

1 模式和资料

1.1 模式介绍

本文采用区域大气动力-化学耦合模式WRF-Chem V3.6进行同化对比实验。该模式该模式由美国国家大气和海洋局(NOAA)研发的新一代气象-大气化学耦合空气质量模式，是完全可压缩的非静力中尺度模式，采用Arakawa-C网格，利于在高分辨模拟中提高准确性，铅直方向采用地形追随的静力气压坐标，重点考虑1～10km左右的大气运动。与离线空气质量模式不同，该模式同步计算气象过程和大气污染过程，它的化学和气象过程使用相同的水平和垂直坐标系，相同的物理参数化方案，不存在时间上的插值，并且能够考虑化学对气象过程的反馈作用。事实上，在实际大气中化学和气象过程是同时发生的，并且能够互相影响，如气溶胶能影响地气系统辐射平衡，气溶胶作为云凝结核，能影响降水，而气温、云和降水对化学过程也有非常强烈的影响。因此，WRF-Chem实时计算污染物浓度对大气辐射的影响，考虑了污染物与气象场的相互作用机制，能够模拟再现一种更加真实的大气环境。研究表明，WRF-Chem模式可以合理模拟和预报出中国区域天气状况、污染物浓度、分布、传输和化学转化等[11-12]。

模式网格分辨率设置为5km，模式投影中心点取为(118.5°E，32.5°N)，经向细分为115个网格，纬向细分为90个网格，包括江苏大部分地区和安徽东部地区，垂直方向分为30层，模式顶层气压为50hPa。表1为模式设置的参数化方案：气溶胶和化学过程参数化方案选取GOCART气溶胶方案和RADM2-KPP气相化学方案，RADM2机制是第二版本的区域酸性沉降机制，该机制综合考虑了化学过程、预报精度和计算机资源，被广泛应用与各种大气化学模式当中，RADM2中无机物包括14个稳定物种，4个活泼的中间物种和3个足量稳定物种(氧气、氮气和水汽)，有机物包括26类稳定物种和16类过氧自由基，使用一种准稳态近似(QSSA)的方法求解22个诊断物种、38个预报物种的浓度；选用TUV光化学过程方案，该方案考虑了21个物种的光分解反应，对于不同波长和物种的光分解系数，考虑反应物分子特性(吸收系数和量子产率)；陆面过程选用Noah方案，考虑城市五岛效应影响(热岛效应，干岛效应，湿岛效应，雨岛效应和浑浊岛效应)，同时采用UCM城市冠层方案，该方案对二维结构的街区进行参数化，来计算街区的热力特性，考虑建筑物的阴影和建筑物的反射作用，区分路面、墙面和屋顶的热力作用，对城市热力和动力过程模拟效果较好；考虑气溶胶光化学反应，模式采用RRTMG方案长短波辐射反馈，触发气溶胶消光系数输出[13-14]；同时启动Lin方案云微物理反馈机制，考虑水汽、云水、雨、云冰、雪和霰的预报，模拟云中短波和光解过程，其他过程如表1所示。

表1 WRF-Chem模式参数的本地化设置

1.2 数据资料

观测数据：江苏省站、苏州站、徐州站和连云港站4站532nm波长激光雷达垂直消光系数，观测时段从2016年1月～2017年1月，时间分辨率5min，垂直高度分辨率为7.5m，观测范围从地面到15km。

污染物监测数据：颗粒物浓度(PM10、PM2.5)和AQI数据来自国控站逐小时监测数据，数据经过质量控制，可信度较高。

模式资料：气象初始场及侧边界资料来自GFS 0.25°× 0.25°再分析资料；模式人为排放源使用过去40年全球对流层再分析排放数据(RETRO)和全球大气研究排放数据(EDGAR)，两套数据都提供几种温室气体(CO2，CH4，N2O等)和一些初级粒子的年排放数据，分辨率分别为0.5°× 0.5°和1°× 1°；模拟区域地形和地表信息分别来自模式自带GTOPO30和MODIS数据，分辨率为1km × 1km[23-24]。

1.3 预报评估方法

分别采用相关系数(R)、平均绝对误差(MD)和风险评分(Ts)进行预报结果评估[25]，计算公式如下：

(1)

(2)

(3)

式中，N表示样本数，Oi为观测值，Mi为模型反演值，i=1,2,3,...N；n11表示预报正确次数，n21表示空报次数，n12表示漏报次数，n11+n21为总预报次数，n11+n12为观测总次数。

2 同化方案及实验设计

2.1 同化方案设计

本系统基于空气质量模式WRF-Chem，利用江苏省激光雷达观测数据、监测站污染物浓度数据和常规气象数据，基于AOD反演算法进行三维变分同化。其技术实现涉及两个基本方面：(1)观测信息的输入(观测增量)；(2)观测信息的权重与传播平滑。针对这两部分，本文设计了相应的同化系统，如图1所示，系统分为数据预处理模块和同化模块，前者通过读取激光雷达消光系数、污染监测数据和常规气象观测数据分别进行质量控制，并以年-月-日命名存入对应文件夹，对消光系数进行垂直积分获取AOD，通过反演计算，利用相关性、绝对误差等评估指标，建立权重系数，综合得到相应的PM2.5浓度，实现同化所需观测信息的格式化载入；后者即同化模块，导入气象再分析、地理信息、源清单等数据，进行模式预处理，通过模式冷启动，生成未同化初始场，启动三维变分同化技术，载入观测数据订正模式初始场，随后更新侧边界，启动模式进行积分运算。两个模块通过反演PM2.5浓度，实现实况数据载入和模式同化。

图1 激光雷达同化系统流程图

2.1.1 数据预处理模块设计

激光雷达消光系数、污染监测数据和常规气象观测数据的质量控制包括缺失值和异常值剔除、数据整合和规范化排列、保存等，这里不做赘述。

下面对AOD反演算法进行详细说明。大气气溶胶光学深度表示垂直气柱内各个高度气溶胶消光系数的积分，由于大气中不同高度气溶胶浓度受大气相对湿度、辐射强度等大气参数影响较大，导致大气垂直方向上气溶胶浓度呈现不均、复杂的廓线分布，因此直接利用AOD估算地面颗粒物浓度必然导致准确率下降。本文采用标高法获取地面气溶胶消光系数进行颗粒物浓度反演，具体公式如下[26]：

(4)

由上述公式变换，可得地面消光系数σ0：

(5)

本文利用地面消光系数σ0建立PM2.5浓度的单要素拟合模型：y=-52.08x+111.62 ，y为地面PM2.5浓度，x为地面消光系数σ0。

为考虑其他气象要素的影响[27-28]，本文同时基于多种气象要素，包括温度、相对湿度、风速等，同步建立地面消光系数与PM2.5浓度的多元拟合模型：

y=-35.24x1+131.07x2+2.93x3-108.36x4-172.88

(6)

其中y为地面PM2.5浓度，x1至x4分别为地面消光系数、近地面温度和相对湿度、近地面风速。

对比两模型的拟合结果，如图2所示，单要素模型得到的PM2.5与实况相关性为0.27，回归方差占总方差比重仅为29%，表明单要素线性回归模型的数据损耗较为明显。多要素多元拟合模型反演的PM2.5浓度值相关性为0.71，且拟合后数据围绕斜率为1的直线分布，主要落区分布在原始数据的1/2至2倍之间，平均绝对误差大幅下降，表明多元拟合反演模型能够更好地反映激光雷达数据和PM2.5的关系。因此，本文采用多元拟合模型方案作为PM2.5的预测模型，建立观测算子。

图2 单要素线性模型和多要素线性模型反演PM2.5的分布

2.1.2 同化模块设计

资料同化就是在功能上实现按照某种算法通过计算机程序(无需人为干预)，将不规则分布的观测资料结合天气动力学和观测资料特征而得到某一预定规则网格点上最可能的值，即表示尽可能准确的大气状态的分析场，来提供数值天气预报模式的初值条件[29]。

通过数据预处理模块，建立PM2.5的多元拟合预测模型，实现实时观测数据反演、订正输出，作为同化模块的观测算子。基于贝叶斯定理，在先验背景场误差和观测误差是不相关的假定下，由分析场似然概率最大的解转化为等价于目标泛函数最小的解：

-y)TR-1(H(x)-y)

(7)

式中：x是同化分析场，xb是背景场；B是背景误差协方差；H是观测算子；y是观测向量；R是观测误差协方差。本项目中的控制变量包括流函数、非平衡速度势、非平衡温度、非平衡地表气压、假相对湿度等资料。

通过三维变分同化，将观测信息载入模式初始场，改善模式初始条件，随后更新侧边界条件，防止模式运行过程溢出，最后启动模式进行积分运算。

2.2 实验方案设计

2016年12月7日～10日，南京地区出现一次以细颗粒物为主要污染物的中度至重度污染过程：7日开始颗粒物浓度逐渐增大，出现轻度污染；8日到9日凌晨PM2.5浓度达到最大，且浓度值维持稳定；9日下午前后PM2.5浓度逐渐下降，至10日污染过程结束。

以该污染过程为模拟对象，本文设计了2组模拟实验方案，具体如下。

方案1为控制实验：即不采用同化方案，从2016年12月6日20时向后积分WRF-Chem模式做84小时预报，到2016年12月10日08时结束，其中前12小时作为spin-up，对气象和化学过程进行充分适应，取后72小时用于研究分析。

方案2为同化实验：利用本文前述反演算法和高度订正，得到同化时刻PM2.5浓度，再将PM2.5浓度作为同化观测场进行同化。以12月6日20时起WRF-Chem的12小时预报场作为背景场，同化7日08时的PM2.5资料，将同化得到的分析场作为初始场，利用restart进行热启动，更新模式的边界条件，向后积分72小时，其结果用于模拟分析。

3 结果分析

3.1 同化后分析场分析

图3是反演同化后PM2.5浓度的分析场与同化前原始场的对比。对比发现，苏北地区原始场模拟偏高，经过同化订正后徐州地区PM2.5浓度降幅达到50μg/m3，连云港地区也有30μg/m3降幅，同化有效半径达到50km左右。苏州地区模拟的原始场与实况接近，同化订正效果不明显。南部地区由于缺少观测，同化对原始场的调整幅度相对较小。

(图a：同化后，图b：未同化，单位：μg/m3)

(沿32.1°N做纬向剖面，单位：μg/m3)

如图4所示，沿32.1°N做纬向剖面，对比同化前后分析场在垂直方向的浓度变化。同化后，垂直方向调整明显，在119.5°E东西两侧分别形成差值中心，西部为负值区，表明同化后分析场减小，东部为正值区，表明同化后分析场增大；西部调整幅度较大，垂直影响较高，达到800hPa(2 000m左右)，东部调整较小，垂直方向达到750hPa(1 500m左右)；高度升高，浓度调整的影响也随之减弱。

3.2 预报结果对比分析

图5为南京市PM10、PM2.5和AQI 72小时预报结果和观测值对比时间序列图(徐州和苏州地区时间序列图略)，模式模拟出污染物浓度变化趋势，相关性达到0.55以上，但未经过同化的颗粒物浓度普遍偏高，最高达到2倍以上，主要因为排放源清单分辨率较低和未进行动力降尺度造成此现象；利用近地面反演数据订正初始场，同化后结果得到较大改善。

实验从7日08时开始restart，气象和化学模式需要磨合，所以前几小时模拟偏差较大；随着预报时长的增加，初始场自身的误差导致预报不确定性逐渐降低，同化效果逐渐显现，且前36小时预报与实况走势差异较小；随后排放源、气象背景场及空气质量模式本身的不确定性逐渐成为导致预报偏差的主要原因，同时近地层同化效果逐渐减弱，初始场垂直订正的同化对预报时长的影响放大模式本身不确定，导致后期模拟偏差增大；这种偏差突出表现在PM10浓度预报，PM2.5同化效果较好，但峰值偏差依然较大。

此次过程以PM2.5污染为主，AQI整体走势与PM2.5浓度变化一致。反演同化后，预报与观测相关性提高了0.1左右，准确率提高了10%左右，PM2.5的平均绝对误差下降了20，PM10下降了50左右；南京与徐州地区颗粒物浓度变化幅度较大，整个过程出现多波峰波谷，模拟结果较差，同化后预报改善明显；苏州地区前期污染严重，后期污染消退，整个过程污染物浓度变化单一，波动幅度较小，所以模拟结果优于其他两个地区，同化后平均绝对误差明显减小，但相关性和准确率调整较小(表2所示)。

图5 南京市区2016年12月7日～10日PM2.5(a)、PM10(b)和AQI(c)72小时预报结果与观测浓度对比时间序列图

表2 实验结果评估对比

续表2

为客观检验实验对污染事件的预报准确度，根据空气质量浓度等级(表3所示)，将颗粒物浓度和AQI按照阈值划分为优良、轻度、中度和重度以及严重污染5个等级，采用Ts(Threat score)方法，对不同等级污染预报能力进行分类检验。

表3 空气质量分级浓度

如表4所示，使用污染物浓度等级来定义污染事件，计算过程预报检验数据。统计结果显示，优良等级天气预报评分最高；PM2.5和PM10的中度、重度污染等级评分高于轻度污染预报，表明模式预报对中度和重度污染预报能力较强；另一方面，同化后各污染物的不同等级评分都有显著提升：徐州地区轻度、中度和重度污染的预报改善效果提升明显，幅度达到0.2左右，苏州轻度污染改善明显，达到0.1左右，南京地区也有小幅提升；不同污染物的同化改善效果也有不同，PM2.5浓度改善最为明显，其次是AQI，最后是PM10，表明反演同化直接对PM2.5浓度初始场产生影响，间接影响PM10浓度，此次过程主要为PM2.5污染，AQI浓度收PM2.5影响较大，所以预报结果优于PM10。整体而言，同化后模式对于污染事件的预报准确度有明显提高。

表4 不同等级污染物浓度的分类检验

续表4

4 结 论

4.1 采用高度订正和多元拟合算法，建立激光雷达AOD与近地面PM2.5关系，引入观测算子，通过同化订正模式初始场，影响模式预报结果。

4.2 同化对PM2.5预报改善优于PM10，反演同化分析场，对PM2.5浓度进行订正，间接影响PM10，所以PM2.5预报改善效果明显优于PM10；同化后AQI预报也有显著改善，相关性提高了0.05～0.1，预报准确率提高了10%～20%。

4.3 不同地区，由于污染过程的特征不同，同化改进效果不同，总体而言，模拟结果越偏离监测值、污染过程越复杂，同化改善效果越明显。

4.4 本文通过实验对比，发现同化后分析场调整明显，对颗粒物和AQI预报改善效果显著，表明激光雷达数据在同化方面表现良好，在模式同化领域有广阔的应用空间。