一种高时空分辨率大气污染扩散迁移模拟方法

技术领域

本发明涉及大气污染监测技术领域，具体涉及中小尺度下大气污染扩散迁移的数值模拟，更具体地涉及一种高时空分辨率大气污染扩散迁移模拟方法。

背景技术

由于大气污染物及其前体物的大量排放，使众多城市时常遭受严重的大气污染问题。尤其是近年来大气中的PM

对于大气污染物扩散规律的数值研究模型通常可分为三类：高斯扩散模型，欧拉扩散模型和拉格朗日粒子扩散模型。其中高斯扩散模型发展的最早，也是相对应用较为成熟的一种经验模型(如《高斯模型在中小城市多点源大气扩散模拟中的应用研究》，哈尔滨师范大学硕士学位论文，2020年；《基于变化风场下燃气泄漏扩散GIS系统研究》，重庆大学硕士学位论文，2018年)。高斯模型是一种经典便捷的解析模型，但该模型为定常模型，对于气象及地形(下垫面)要素的考虑较为简单，因此很难适用于复杂非定常条件下的数值模拟。欧拉扩散模型将计算区域进行三维空间离散，通过联立求解不同离散网格的守恒方程来获得污染物在不同网格内的组分特性，较流行的计算模式有WRF-CMAQ中尺度模式(《基于WRF-CMAQ模式的沈阳市大气污染特征研究》，沈阳航空航天大学硕士学位论文，2020年)及CFD模式(《基于多相质点网格法的大气污染仿真模拟》，科学技术与工程，2020年第20卷，第14期)等。欧拉扩散模型的优势在于可对某些复杂的物理过程进行参数化处理，并可进行精细化建模，但欧拉模型普遍较为耗费计算资源，且其数值计算收敛性上存在一定的限制，因此该模型无法用于快速评估。拉格朗日粒子模式的核心是追踪粒子物理行为的时空变化，其中比较有代表性的有HYSPLIT模式(《基于HYSPLIT模型的大气CO2传输研究》，中国地质大学(北京)硕士学位论文，2014年)以及FLEXPART模式(《北京地区FLEXPART模式适用性初步研究》，环境科学学报，2010年8月，第30卷第8期)等。

其中，FLEXPART模型作为一种发展较成熟、移植方便的拉格朗日粒子扩散模型，其复杂程度介于简单的轨迹计算和大气化学模式之间。与欧拉扩散模型相比，FLEXPART模型具有时空分辨率高、计算时间短以及后向轨迹追踪等优点。此外，FLEXPART模型还考虑了湍流和中小尺度对流过程，因此可较为准确地模拟示踪物的扩散、沉降以及衰减过程。

然而，目前FLEXPART模型对于街区尺度下(计算区域约10×10km量级，空间分辨率100m量级)的大气污染物扩散过程的模拟仍存在不足，其中一个主要原因在于用于驱动FLEXPART进行污染物扩散模拟的中尺度气象模式WRF的空间分辨率通常只能达到3-5km量级，无法驱动FLEXPART进行街区尺度的污染物扩散模拟。因此，研究如何利用FLEXPART模型对街区尺度下大气污染物扩散的时空分布特征进行高时空分辨率的捕捉，由此对大气污染扩散、沉降等过程进行深入分析，具有其重要意义。

发明内容

本发明的目的是解决现有大气污染物扩散规律的数值研究采用FLEXPART模型时，因用于驱动FLEXPART进行污染物扩散模拟的中尺度气象模式WRF的空间分辨率低，无法驱动FLEXPART进行街区尺度的污染物扩散模拟的技术问题，而提供一种高时空分辨率大气污染扩散迁移模拟方法。

本发明的构思为：在中尺度三维气象模拟的基础上，首先采用降尺度气象诊断模型进行动力学降尺度诊断模拟，将中尺度气象模型得到的3-5km空间分辨率的逐小时气象场数据在空间上降尺度至百米空间分辨率。随后利用耦合接口程序，将街区尺度三维气象场数据转换为FLEXPART模型可读取的数据格式，以该街区尺度气象场数据驱动FLEXPART模型进行下一步大气污染物扩散过程的数值运算，从而实现街区尺度下污染物扩散的高时空分辨率(空间分辨率百米，时间分辨率分钟量级)模拟。

为完成上述发明目的、实现上述方面构思，本发明采用的技术方案为：

一种高时空分辨率大气污染扩散迁移模拟方法，其特殊之处在于，包括以下步骤：

步骤一：采用中尺度气象模式，在拟研究时段以设定的时间间隔对拟研究区域进行中尺度气象场模拟，获得在拟研究时段包含拟研究区域在内的中尺度范围内三维气象场数据；

所述三维气象场数据包括气温、气压、湿度、风速、风向和降水量；

步骤二：根据步骤一获得的中尺度范围内三维气象场数据，利用降尺度气象诊断模式在拟研究区域进行动力学降尺度诊断模拟，得到精细化三维气象场数据；

步骤三：保持时间间隔和空间分辨率不变，采用耦合接口程序，将步骤二得到的精细化三维气象场数据的结果与FLEXPART模型运行所需的气象场变量进行对照映射；

步骤四：在FLEXPART模型的参数脚本中设定污染的源项参数；

所述源项参数包括输出文件路径、三维气象数据路径、计算域起止时间、时间参数、网格参数、排放物种参数、污染物释放位置、污染物释放时间、污染物释放强度；

步骤五：结合步骤三对照映射结果和步骤四设定的源项参数，驱动FLEXPART模型进行迁移扩散模拟，得到大气污染物不同时刻的浓度分布、沉降量信息。

进一步地，步骤三中，所述采用耦合接口程序，将步骤二得到的精细化三维气象场数据的结果与FLEXPART模型运行所需的气象场变量进行对照映射，具体为：

3.1、使用CALMET2NETCDF程序，将步骤二得到的精细化三维气象场数据的结果文件转换为三个NETCDF格式文件：CALGRID2D，CALMET2D，CALMET3D；其中CALGRID2D为网格信息文件，CALMET2D为地面气象场数据，CALMET3D为三维气象场数据；

3.2、采用耦合接口程序，将CALGRID2D，CALMET2D，CALMET3D三个文件中的气象场变量与FLEXPART模式运行所需的气象场变量进行对照映射。

进一步地，步骤3.2中，所述CALGRID2D，CALMET2D，CALMET3D三个文件中的气象场变量与FLEXPART模式运行所需的气象场变量进行对照映射的变量映射关系具体为：

变量纬度：FLEXPART模式中的XLAT对应CALGRID2D中的lat；

变量经度：FLEXPART模式中的XLONG对应CALGRID2D中的lon；

变量基础位势：FLEXPART模式中的PHB＝0对应CALMET3D中0；

变量扰动位势：FLEXPART模式中的PH对应CALMET3D中global attributes:VGLVLS以及CALGRID2D中ELEV；

变量基础气压：FLEXPART模式中的PB＝0对应CALMET3D中0；

变量风速分量U：FLEXPART模式中的U对应CALMET3D中U；

变量风速分量V：FLEXPART模式中的V对应CALMET3D中V；

变量风速分量W：FLEXPART模式中的W对应CALMET3D中W；

变量气温T：FLEXPART模式中的T对应CALMET3D中T；变量10m风速分量U10：FLEXPART模式中的U10对应CALMET3D中U插值，CALMET设置时设置2m和10m高度层；

变量10m风速分量V10：FLEXPART模式中的V10对应CALMET3D中V插值，CALMET设置时设置2m和10m高度层；

变量2m气温：FLEXPART模式中的T2对应CALMET3D中T插值，CALMET设置时设置2m和10m高度层；

变量2m比湿：FLEXPART模式中的Q2对应CALMET2D中IRH：相对湿度(％)计算获取；CALMET设置时设置2m和10m高度层；

变量地表气压：FLEXPART模式中的PSFC对应CALMET2D中RHO：空气密度p＝RHO×Rd×Tv计算获取；

变量大尺度降水：FLEXPART模式中的RAINNC对应CALMET2D中RMM；

变量陆地使用类型：FLEXPART模式中的LU_INDEX对应CALGRID2D中ILANDU；

变量摩擦速度：FLEXPART模式中的UST对应CALMET2D中USTAR；

变量边界层高度：FLEXPART模式中的PBLH对应CALMET2D中ZI；

变量短波辐射：FLEXPART模式中的SWDOWN对应CALMET2D中QSW。

进一步地，步骤五中，所述迁移扩散模拟时，采用湍流参数化方法来估算对流边界层中的湍流混合作用及边界层高度；采用基于三维湍动能的对流参数化方法来计算中小尺度对流过程。

进一步地，步骤一中，所述中尺度模式的气象场模拟采用三层嵌套方式，且最内层需覆盖拟研究区域，其中最内层空间分辨率为3km。

进一步地，步骤二中，所述精细化三维气象场数据的街区尺度为10km×10km，空间水平网格分辨率为100m，垂直方向分为十层以上。

进一步地，步骤一中，所述设定的时间间隔为1小时。

进一步地，步骤五中，所述FLEXPART模型进行迁移扩散模拟得到的结果的空间分辨率为100m，垂直最高分辨率20m，时间分辨率达到分钟量级。

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

1、本发明提供的高时空分辨率大气污染扩散迁移模拟方法，基于降尺度耦合技术，该方法在中尺度三维气象模拟的气象数据的基础上，采用采用降尺度气象诊断模型进行动力学降尺度诊断模拟，获得街区尺度下的三维气象场分布参数，实现了中尺度气象模型的逐小时3-5km空间分辨率的气象场数据在空间上降尺度至百米空间分辨率。

2、本发明提供的高时空分辨率大气污染扩散迁移模拟方法，利用自编程耦合接口，耦合降尺度三维气象模型和FLEXPART模型，实现街区尺度下大气污染物扩散的高时空分辨率的数值模拟，具有时空分辨率高的特点，在针对街区尺度的物理问题模拟上，时间分辨率可以达到分钟量级，空间水平分率可达到百米量级，空间垂直分辨率可达20m。

3、本发明在污染物的扩散迁移数值模拟过程中，考虑模型湍流和中小尺度对流，采用湍流参数化方法来估算对流边界层中的湍流混合作用及边界层高度，既兼顾了垂直速度中湍流的倾斜度，也兼顾了空气密度的垂直梯度，使得大气污染扩散迁移模拟湍流的因素更为完备，可以更好地计算对流边界层中的混合状态；采用基于三维湍动能的对流参数化方法来计算中小尺度对流过程，可以更为精确捕捉污染物的大气扩散过程。

4、本发明方法耗费计算资源小，计算时间少。在已获取降尺度三维气象场数据的条件下，单机运行5分钟内即可完成全部扩散迁移计算过程，因此可以用于事故后果的快速评价。

附图说明

图1为本发明高时空分辨率大气污染扩散迁移模拟方法的流程图；

图2为本发明实施例中降尺度模型计算得到的分辨率为100m×100m的风场分布情况；其中(a)为污染物释放时刻计算区域地表风矢量图；(b)为污染物释放后1h计算区域地表风矢量图；(c)为污染物释放后2h计算区域地表风矢量图；(d)为污染物释放后3h计算区域地表风矢量图；(e)为污染物释放后4h计算区域地表风矢量图；(f)为污染物释放后5h计算区域地表风矢量图；

图3为本发明实施例中使用耦合接口转换后的FLEXPART读入气象文件；

图4为本发明实施例中FLEXPART模型输出的街区尺度下不同时刻污染物的浓度分布云图；其中(a)为污染物释放时刻计算区域地表浓度分布云图；(b)为污染物释放后1h计算区域地表浓度分布云图；(c)为污染物释放后2h计算区域地表浓度分布云图；(d)为污染物释放后3h计算区域地表浓度分布云图；(e)为污染物释放后4h计算区域地表浓度分布云图；(f)为污染物释放后5h计算区域地表浓度分布云图；

图5为本发明实施例中FLEXPART模型输出的街区尺度下不同时刻污染物的沉降量分布云图；其中(a)为污染物释放时刻计算区域地表沉降分布云图；(b)为污染物释放后1h计算区域地表沉降分布云图；(c)为污染物释放后2h计算区域地表沉降分布云图；(d)为污染物释放后3h计算区域地表沉降分布云图；(e)为污染物释放后4h计算区域地表沉降分布云图；(f)为污染物释放后5h计算区域地表沉降分布云图。

具体实施方式

为使本发明的目的、优点和特征更加清楚，以下结合附图和具体实施例对本发明提出的一种高时空分辨率大气污染扩散迁移模拟方法作进一步详细说明。本领域技术人员应当理解的是，这些实施方式仅仅用来解释本发明的技术原理，目的并不是用来限制本发明的保护范围。

如图1所示，本实施例提供的高时空分辨率大气污染扩散迁移模拟方法，具体包括五个步骤：

第一步：采用中尺度气象模式，在拟研究时段对拟研究区域的进行中尺度气象场模拟，获得在拟研究时段包含拟研究区域在内的中尺度范围内三维气象场数据；

本实施例中，中尺度模式的气象场模拟采用三层嵌套方式，最内层空间分辨率为3km，空间分辨率分别是27km，9km和3km，最内层需覆盖拟研究区域。

中尺度模式中的物理参数化方案选择如表1所示。

表1WRF模式物理过程设置

中尺度范围内逐小时三维气象场数据包括气温、气压、湿度、风速、风向、降水量。

第二步：利用降尺度气象诊断模式，将中尺度气象模式输出的气象数据结果在空间上进行动力学降尺度诊断模拟。

降尺度气象诊断模式中设置中尺度气象模式输出的气象场为初猜气象场，计算域(街区尺度)需覆盖大气污染物的排放源，计算域范围取10km×10km，水平网格分辨率取100m，垂直方向分为十层以上，最终输出逐小时的空间水平分辨率为100m左右的精细化风场，如图2所示，污染物释放0～5小时计算区域地表风矢量图，(a)～(f)分别为污染物释放时刻以及释放后每小时的表风矢量图。

降尺度气象诊断模式的核心参数设置如表2所示。

表2降尺度气象诊断模式核心参数设置列表

/>

第三步：使用CALMET2NETCDF程序，将降尺度气象诊断模式输出结果文件转换为三个NETCDF格式文件：CALGRID2D，CALMET2D，CALMET3D。其中CALGRID2D为网格信息文件，CALMET2D为地面气象场数据，CALMET3D为三维气象场数据。

进一步使用自编写的模型耦合接口程序模块，将CALGRID2D，CALMET2D，CALMET3D三个文件中的气象场变量与FLEXPART模式运行所需的气象场变量进行对照映射，即将CALMET模式输出的降尺度气象场信息转换成FLEXPART模式所需的输入文件(如图3所示)，以驱动FLEXPART进行污染扩散迁移模拟。

该耦合过程需要进行CALMET模式输出气象场与FLEXPART模式所需气象场间的变量映射，两种模型间的变量映射关系如表3所示：

表3FLEXPART模式运行所需变量与CALMET结果中变量映射关系

/>

第四步：在flexwrf.input文件中进行计算参数的设置。具体设置内容包括输出文件路径、CALMET气象数据路径、计算域起止时间、时间相关参数设定、网格相关参数设定、排放物种相关参数设定、污染物释放有关信息(释放位置、释放时间、释放强度等)的设定。

第五步：使用FLEXPART模式读入转换格式后的降尺度气象场，基于定义好的源项参数信息进行污染物的扩散迁移数值模拟。

在污染物的扩散迁移数值模拟过程中，湍流方面采用湍流参数化方法来估算对流边界层中的湍流混合作用及边界层高度，该方法既考虑了垂直速度中湍流的倾斜度，也考虑了空气密度的垂直梯度，因此该方法考虑湍流的因素更为完备，可以更好计算对流边界层中的混合状态；对流方面采用基于三维湍动能的对流参数化方法来计算中小尺度对流过程，因此可以更为精确捕捉污染物的大气扩散过程。

该步骤需设置FLEXPART模式参数文件flexwrf.input中的LDIRECT参数为1，以控制FLEXPART模拟示踪物由源区向周围扩散的过程。

FLEXPART模拟输出的变量有：浓度(CONC，污染物浓度，单位：ng·m

图4为污染物释放0～5小时计算区域地表浓度分布云图，(a)～(f)分别为污染物释放时刻以及释放后每小时计算区域地表浓度分布云图。

图5为污染物释放0～5小时计算区域地表沉降分布云图，(a)～(f)分别为污染物释放时刻以及释放后每小时计算区域地表沉降分布云图。

通过结果的后处理作图，可获取源排放污染物在大气中的扩散、沉降等规律。该FLEXPART模式结果的计算域范围及空间分辨率与CALMET模式相同，计算域范围取10km×10km，水平方向空间分辨率为100m，而时间分辨率则达到分钟量级。

本发明可以克服现有大气污染扩散模型中时空分辨率低的问题，从而在对中尺度气象模式的模拟结果进行动力学降尺度的基础上，实现FLEXPART模型在工业园区、城市街区等环境下对大气污染物的扩散、传输、沉降等过程的精确模拟和捕捉。

最后应说明的是：以上实施例仅用以说明本发明的技术方案，而非对其限制；尽管参照前述实施例对本发明进行了详细的说明，本领域的普通技术人员应当理解：其依然可以对前述实施例所记载的技术方案进行修改，或者对其中部分或者全部技术特征进行等同替换；而这些修改或者替换，并不使相应技术方案的本质脱离本发明技术方案的范围。