CN111650672B - 采用时间堆栈实现空气污染物大气扩散快速预测的方法 - Google Patents

Abstract

本发明提供一种采用时间堆栈实现空气污染物大气扩散快速预测的方法，解决现有随机游走方法模拟大气污染物扩散，存在占用内存过大、计算时间过长、计算效率低的问题。该方法包括以下步骤：1)设初始化条件：模拟粒子数、截断时间、计算范围、记录网格分辨率和污染源的空间分布；2)污染源抽样；3)采用时间堆栈方法跟踪模拟标记粒子的运动轨迹，具体为在每一个时间步中，计算所有标记粒子经过一个时间步Δt后的位置，并将每一个粒子信息存储在堆栈中；计算完某一时间步内的所有粒子运动后，根据粒子驻留时间给出该时刻每个空间网格上污染物浓度值；4)根据不同时刻所有网格点的空间浓度分布，得到污染物整个扩散过程中浓度分布的时空图像。

Description

采用时间堆栈实现空气污染物大气扩散快速预测的方法

技术领域

本发明涉及一种空气污染物扩散模拟方法，具体涉及一种采用时间堆栈实现空气污染物大气扩散快速预测的方法。

背景技术

大气是人类赖以生存的最基本环境要素之一，它的状态和变化时时处处影响到人类的活动与生存。随着工业化和城市化的快速发展，大气污染已成为全球性的社会公害。另外，近年来，有毒或放射性气体的泄漏扩散事故也呈逐渐增多之势，给人民的生活造成严重影响，给国家财产造成重大损失。研究空气污染物在大气中的扩散，对污染物扩散造成的环境影响进行评价，对事故的实时环境后果评估和应急决策有积极意义。

空气污染物的散布主要是在大气边界层的湍流场中进行的，其扩散行为主要取决于大气的平均风速和湍流运动，前者使其在大气中产生空间位移，后者则使它们与周围空气相混合。最初，高斯扩散模式是污染物大气扩散模拟的主要理论，该模式物理概念清晰，比较适合于稳态流场中小尺度范围的计算，特别是其有很高的计算效率和空间分辨率，因而至今仍然是最受欢迎的模式之一，如由美国国家环保局联合美国气象学会组建法规模式改善委员会(AERMIC)开发的AERMOD模型(Cimorelli A J等发表在《Journal of AppliedMeteorology》期刊2005年第44卷第5期第682页，“AERMOD:A dispersion model forindustrial source applications.Part I:General model formulation and boundarylayer characterization”)和英国剑桥大学环境研究中心(CERC)开发的ADMS模型(CERC于2002年发表的用户手册，“ADMS-EIA User Manual”)均属于高斯型扩散模型，它们被广泛应用于研究局部区域大气环境污染问题。考虑到大气湍流具有高度的随机性，空气污染物在大气中的扩散行为可以使用随机游走方法来模拟(Ley A J发表在《AtmosphericEnvironment》期刊1982年第16卷第12期第2799页，“A random walk simulation of two-dimensional turbulent diffusion in the neutral surface layer”)。随机游走方法不但能够模拟均匀湍流场中的污染物大气扩散问题，还可以很好的模拟非均匀、非定常和有较强切变的复杂流场条件下的污染物扩散(Rotach M W等发表在《Quarterly Journal ofthe Royal Meteorological Society》期刊1996年第122卷第530期第367页，“A two-dimensional stochastic Lagrangian dispersion model for daytime conditions”)。近年来，随着大气科学和边界层理论研究的不断深入以及计算机科学的迅速发展，随机游走方法在理论研究(LIU Feng等发表在2012年Conference on environmental pollutionand public health会议报告，“Some theoretical aspects of atmospheric dispersionrandom walk model”)和应用领域(WANG Peng等发表在《Environmental Monitoring andAssessment》期刊2011年第172卷第507页，“Random-walk model simulation of airpollutant dispersion in atmospheric boundary layer in China”)都相当活跃，并取得令人鼓舞的进展。

但采用随机游走方法模拟空气污染物大气扩散时，需要跟踪模拟大量的标识粒子才能达到满足误差要求的计算结果，因此相较其它计算方法，该方法的计算效率较低，特别是在大尺度污染物扩散模拟中，由于粒子跟踪历时较长、高时空分辨率的气象数据占用内存空间大，导致计算效率非常低、计算机内存需求巨大。因而，如何提高随机游走模拟方法的计算效率，并减少计算对内存的需求是目前研究和实践应用中的一个难点。

发明内容

为了解决现有随机游走方法模拟大气污染物扩散时，模拟所使用的气象数据占用内存过大、计算时间过长、计算效率非常低的技术问题，本发明提供了一种采用时间堆栈实现空气污染物大气扩散快速预测的方法。

为实现上述目的，本发明提供的技术方案是：

一种采用时间堆栈实现空气污染物大气扩散快速预测的方法，其特殊之处在于，包括以下步骤：

1)计算条件初始化处理

设置初始化条件，所述初始化条件包括计算条件参数和污染源条件参数；

所述计算条件参数包括模拟粒子数N_S、截断时间T_CUT、计算范围和记录网格分辨率；

所述污染源条件参数包括污染源的空间分布；

2)污染源的抽样

根据污染源的空间分布抽样源粒子在初始时刻t₀的位置X_j(t₀)，并对抽样得到的源粒子进行标记为j，j＝1，2……N_S；

其中，X_j(t₀)的位置坐标为(x_j0，y_j0，z_j0)，定义经度方向为x方向，纬度方向为y方向，竖直方向为z方向；

3)采用时间堆栈方法跟踪模拟标记粒子的运动轨迹；

3.1)读入t₀时刻的气象参数，所述气象参数包括边界层厚度H_b、下垫面粗糙度r₀和平均风速V；根据初始时刻t₀标记粒子的位置坐标，插值计算标记粒子所处位置的气象参数，并用大气科学和行星边界层理论计算边界层湍流扩散参数σ与拉格朗日时间尺度T_L；

3.2)计算标记粒子t₁时刻位置和污染物浓度值

3.2.1)根据大气输运和湍流扩散，分别计算初始时刻t₀的所有标记粒子经过时间步Δt后的t₁时刻位置X_j(t₁)，并将t₁时刻每一个标记粒子信息存储在堆栈中；

3.2.2)计算t₁时刻每个空间网格上的污染物浓度值，根据每个空间网格上的污染物浓度值，确定并输出t₁时刻污染物的浓度空间分布；

污染物浓度值计算公式如下：

式中，Q为释放污染物总量，i代表记录网格编号，ΔV代表网格体积，T_i,j为该时间步内第j个粒子在第i号记录网格中的停留时间；

3.3)读入t₁时刻的气象参数，所述气象参数包括边界层厚度H_b、下垫面粗糙度r₀和平均风速V；根据t₁时刻标记粒子的位置坐标，插值计算标记粒子所处位置的气象参数，并用大气科学和行星边界层理论计算边界层湍流扩散参数σ与拉格朗日时间尺度T_L；

3.4)以标记粒子t₁时刻位置X_j(t₁)为下一时间步Δt的初始位置，利用步骤3.2)的计算方法，计算t₁时刻的所有标记粒子经过时间步Δt后的t₂时刻位置X_j(t₂)和t₂时刻每个空间网格上的污染物浓度值，并将t₂时刻每一个标记粒子信息存储在堆栈中，以及输出t₂时刻污染物的浓度空间分布；

3.5)重复利用步骤3.3)和3.4)，直至跟踪时间达到截断时间T_CUT或者计算出的标记粒子t_k时刻位置X_j(t_k)逸出计算边界；

4)得到污染物扩散过程

根据不同时刻所有网格点的空间浓度分布，得到污染物整个扩散过程中浓度分布的时空图像。

进一步地，步骤3.1)和步骤3.3)中，用大气科学和行星边界层理论计算边界层湍流扩散参数σ与拉格朗日时间尺度T_L具体为：

设大气稳定度为中性，则σ与T_L在顺风方向的分量为σ₁和T_L1、水平侧风方向的分量为σ₂和T_L2、垂直方向的分量σ₃和T_L3，计算公式如下：

σ₁＝2.3×u^*

σ₂＝2.0×u^*

其中，z为标记粒子离地高度；u^*为摩擦速度，

式中，κ为卡曼常数，u₀为地面风速，z₀为对应测风高度。

进一步地，步骤3.2.1)中，t₁时刻位置X_j(t₁)计算公式如下：

X_j(t₁)＝X_j(t₀)+V(X_j(t₀),t₀)Δt+V′(X_j(t₀),t₀)Δt

其中，X为标记粒子的坐标，三维分量分别为x,y,z；

V为平均风速，三维分量分别为u,v,w；

V′为湍流脉动速度，三维分量分别为u′,v′,w′；

湍流脉动速度V′遵从马尔科夫假定，即

式中，ξ为符合正态分布的随机数；R_u，R_v，R_w为u，v，w三个速度分量风向的拉格朗日自相关函数，其计算公式如下：

其中，σ_u、σ_v、σ_w为x、y、z三个方向的湍流扩散参数，T_Lu、T_Lv、T_Lw为x、y、z三个方向的拉格朗日时间尺度，σ_u、σ_v、σ_w和T_Lu，T_Lv，T_Lw可通过坐标变换从σ₁，σ₂，σ₃和T_L1，T_L2，T_L3计算得到。

进一步地，步骤3.2.1)和步骤3.4)中，标记粒子信息包括位置坐标和速度信息。

与现有技术相比，本发明的优点是：

1、本发明方法通过采用时间堆栈的方式，突破了现有模拟算法计算内存空间需求过大、计算时间过长的局限性；解决了采用随机游走方法模拟大气扩散时，计算过程中涉及大计算内存需求及海量I/O数据传输而导致计算效率低下的问题。

2、本发明方法能够用于不同尺度范围的空气污染物大气扩散情形，可应用于有毒有害气体区域扩散预测领域，为事故发生后快速制定应急方案提供有效技术支持，也可用于全球范围的污染物输运模拟，为全球环境保护、污染治理提供理论依据。

附图说明

图1为本发明实施例采用时间堆栈实现空气污染物大气扩散快速预测的方法流程图；

图2为本发明实例中10月27日0时的空气污染物分布情况示意图；

图3为距离释放源约1000km位置上处空气中的污染物浓度随时间变化示意图，其中，图a为采用传统模拟方法计算，图b为采用本实施例时间堆栈方法计算；

图4为距离释放源约1000km位置上处空气中的污染物浓度随时间变化采用传统模拟方法和本实施例方法的比较示意图。

具体实施方式

以下结合附图和具体实施例对本发明的内容作进一步详细描述。

本发明采用随机游走方法建立空气污染物扩散计算模型，通过释放并跟踪大量标记粒子的轨迹来模拟空气污染物的大气扩散过程，粒子在流场中按平均风输送，同时又用一系列随机位移来模拟湍流扩散，最后由这些粒子在空间和时间上的总体分布统计给出污染物的时空分布。

一种采用时间堆栈实现空气污染物大气扩散快速预测的方法，包括以下步骤：

1)计算条件初始化处理

设置初始化条件，初始化条件包括模拟粒子数N_S、截断时间T_CUT、计算范围、记录网格分辨率、污染源的空间分布；

2)污染源的抽样

根据污染源的空间分布、泄漏出射速度分布抽样源粒子在初始时刻t₀的位置X_j(t₀)和运动状态，并对抽样得到的源粒子进行标记为j，j＝1，2……N_S；

其中，X_j(t₀)的位置坐标为(x_j0，y_j0，z_j0)，定义经度方向为x方向，纬度方向为y方向，竖直方向为z方向；

3)采用时间堆栈方法跟踪模拟标记粒子的运动轨迹；

在每一个时间步中，计算所有标记粒子经过一个时间步Δt后的位置，并将每一个粒子的信息(时间步结束时粒子位置和速度)存储在堆栈中；计算完某一时间步内的所有粒子运动后，根据粒子驻留时间给出该时刻每个空间网格上的污染物浓度值；

3.1)湍流参数计算

读入t₀时刻的气象参数，气象参数包括边界层厚度H_b、下垫面粗糙度r₀和平均风速V；根据初始时刻t₀标记粒子的位置坐标，插值计算标记粒子所处位置的气象参数，并用大气科学和行星边界层理论计算边界层湍流扩散参数σ与拉格朗日时间尺度T_L；具体的计算方法与大气的稳定度密切相关，当大气稳定度为中性时，则σ与T_L在顺风方向的分量为σ₁和T_L1、水平侧风方向的分量为σ₂和T_L2、垂直方向的分量σ₃和T_L3，计算公式如下：

σ₁＝2.3×u^*

σ₂＝2.0×u^*

其中，z为标记粒子离地高度；u^*为摩擦速度，

式中，κ为卡曼常数，u₀为地面风速，z₀为对应测风高度，z₀通常为10m。

3.2)计算标记粒子t₁时刻位置和污染物浓度值

3.2.1)根据大气输运和湍流扩散，分别计算初始时刻t₀的所有标记粒子经过时间步Δt后的t₁时刻位置X_j(t₁)，并将t₁时刻每一个标记粒子信息(时间步结束时粒子位置和速度)存储在堆栈中，同时记录粒子在记录网格内的驻留时间；

标记粒子在流场运动时，一方面要受到平均风场的作用做空间平移，同时又受到湍流作用而产生一系列随机位移，运动轨迹可表示为

X_j(t₁)＝X_j(t₀)+V(X_j(t₀),t₀)Δt+V′(X_j(t₀),t₀)Δt

其中，X为标记粒子的坐标，三维分量分别为x,y,z；

V为平均风速，三维分量(x，y，z三个方向的平均风速)分别为u,v,w；

V′为湍流脉动速度，三维分量分别为u′,v′,w′；

不同时间步上的湍流脉动速度V′遵从马尔科夫(Markov)假定(t+Δt时刻的运动状态仅与t时刻的运动状态有关)，即

式中，右边第二项代表速度涨落量中的随机部分，ξ为符合正态分布(平均值为0，标准差为1)的随机数；R_u，R_v，R_w为u，v，w三个速度分量风向的拉格朗日自相关函数，其计算公式如下：

其中，σ_u、σ_v、σ_w为x、y、z三个方向的湍流扩散参数，T_Lu、T_Lv、T_Lw为x、y、z三个方向的拉格朗日时间尺度，σ_u、σ_v、σ_w和T_Lu，T_Lv，T_Lw可通过坐标变换从σ₁，σ₂，σ₃和T_L1，T_L2，T_L3计算得到。

3.2.2)计算完t₀到t₁时间步内的所有粒子运动后，根据粒子驻留时间给出t₁时刻每个空间网格上的污染物浓度值，根据每个空间网格上的污染物浓度值，确定并输出t₁时刻污染物的浓度空间分布；

污染物浓度值计算公式如下：

式中，Q为释放污染物总量，N为模拟粒子数，i代表记录网格编号，ΔV代表网格体积，T_i,j为该时间步内第j个粒子在第i号记录网格中的停留时间；

3.3)读入t₁时刻的气象参数，气象参数包括边界层厚度H_b、下垫面粗糙度r₀和平均风速V；根据t₁时刻标记粒子的位置坐标，插值计算标记粒子所处位置的气象参数，并用大气科学和行星边界层理论计算边界层湍流扩散参数σ与拉格朗日时间尺度T_L；

3.4)以标记粒子t₁时刻位置X_j(t₁)为下一时间步Δt的初始位置，利用步骤3.2)的计算方法，计算t₁时刻的所有标记粒子经过时间步Δt后的t₂时刻位置X_j(t₂)和t₂时刻每个空间网格上的污染物浓度值，并将t₂时刻每一个标记粒子信息存储在堆栈中，以及输出t₂时刻污染物的浓度空间分布；

3.5)重复利用步骤3.3)和3.4)，在下一个时间步计算时，从上一时间步的堆栈中取出该标记粒子信息，根据该粒子信息以及当前时间步的气象数据，进行当前时间步的粒子输运模拟计算，并将本时间步结束时的粒子信息存到另一个堆栈中，以此往复，直至跟踪时间达到截断时间T_CUT或者计算出的标记粒子t_k时刻位置X_j(t_k)逸出计算边界；

4)得到污染物扩散过程

根据不同时刻所有网格点的空间浓度分布，得到污染物整个扩散过程中浓度分布的时空图像。

本发明方法对传统模拟方式计算流程进行了改进，从而实现了空气污染物大气扩散的随机游走方法的快速模拟，从而提高空气污染物扩散的模拟效率。该方法主要针对空气污染物随机游走模拟的情形，也可针对涉及高精度时空网格剖分的拉格朗日粒子模式计算模拟。

实施例

全球核安全咨询组织与世界气象组织合作，于1994年在欧洲进行了两次大规模的粒子示踪实验。本实例以对第一次欧洲粒子示踪实验模拟为例，空气污染物释放总量为340kg；释放起始时间为1994年10月23日16时整(国际标准时)，释放时间持续11小时50分；污染物出口流速为7.98g/s；释放点为于东经2°，北纬48°；计算时间步长为1hr；计算截断时间为120hr；计算范围为西经10°～东经40°，北纬20°～70°；记录网格为0.1°×0.1°；模拟粒子数为10⁶个；计算中气象数据的空间分辨率为1.125°×1.125°；气象数据的时间分辨率为6hr。

本实施例采用时间堆栈实现空气污染物大气扩散快速预测的方法流程图见图1所示，首先是建立模型并对计算参数初始化和记录网格剖分，然后是主程序对污染物粒子进行输运计算，包括每一时间步初始时刻粒子状态的获取、时间步内粒子运动轨迹的计算、时间步结束时粒子信息的堆栈存储，以及当每个时间步内的所有粒子计算完后输出计算结果，当计算完当前时间步内的所有粒子运动后就开始计算下一时间步的粒子输运过程，依次往复，直至达到截断时间或污染物粒子全部逸出计算区域为止。具体包括以下步骤：

1)计算条件初始化处理：对模拟粒子数，截断时间，计算范围，记录网格剖分等计算条件参数进行初始化设置；根据污染源的空间分布、泄漏出射速度信息抽样标记粒子初始时刻位置和运动状态。

2)根据计算区域设置：将当前时刻指定区域的气象数据从硬盘读入到缓存中。

3)本时刻粒子初始运动状态：根据当前所处时刻的不同，粒子的初始信息获取方式分以下两种情况：

i.当t＝t₀时，根据污染源的空间分布信息抽样标记粒子初始时刻位置；

ii.当t≠t₀时，从上一时间步的粒子堆栈中取出一条粒子信息，以获取当前时刻粒子初始位置。

4)湍流参数计算：根据抽样粒子的位置坐标，插值计算其所在位置处的气象参数，并用大气科学和行星边界层理论计算湍流参数。

5)粒子运动计算：计算抽样粒子在大气平均风场和湍流的作用下，经过一个时间步Δt后的位置。

6)将当前时间步结束时的粒子位置和运动信息存入堆栈中。

7)重复步骤(3)～(6)直至模拟抽样粒子数达到预定模拟粒子数N_S。

8)统计输出当前时刻污染物空间分布计算结果。

9)下一时刻粒子模拟：计算粒子计数器重置N＝1，时间增长t＝t+Δt。

10)重复步骤(2)～(9)，直至粒子时间t达到截断时间T_CUT或粒子全部逸出边界为止。

图2为计算得到的10月27日0时污染物分布结果。图3给出一个距源释放点约1000km监测点的处地表空气中污染物浓度随时间的变化情况示意图，作为比较，图4中还给出了传统模拟方法和本申请时间堆栈方法的比较结果。从图3和图4中可以看出，采用本发明模拟所得的地表污染物浓度变化曲线与传统模拟方法符合很好。而且根据模拟可知，测点污染物峰值浓度的相对差为0.26％。综上可知，本实施例模拟所得到的结果与传统模拟方法的结果符合很好。此外，由于每个时间步的粒子计算过程中只需读入当前时刻的气象数据，极大节省了计算内存和I/O数据交换耗时，与传统模拟相比，本发明能节省约70％的计算时间、减少约80％的计算占用内存。

以上仅是对本发明的优选实施方式进行了描述，并不将本发明的技术方案限制于此，本领域技术人员在本发明主要技术构思的基础上所作的任何公知变形都属于本发明所要保护的技术范畴。

Claims (2)

1.采用时间堆栈实现空气污染物大气扩散快速预测的方法，其特征在于，包括以下步骤：

1)计算条件初始化处理

设置初始化条件，所述初始化条件包括计算条件参数和污染源条件参数；

所述计算条件参数包括模拟粒子数N_S、截断时间T_CUT、计算范围和记录网格分辨率；

所述污染源条件参数包括污染源的空间分布；

2)污染源的抽样

根据污染源的空间分布抽样源粒子在初始时刻t₀的位置X_j(t₀)，并对抽样得到的源粒子进行标记为j，j＝1，2……N_S；

其中，X_j(t₀)的位置坐标为(x_j0，y_j0，z_j0)，定义经度方向为x方向，纬度方向为y方向，竖直方向为z方向；

3)采用时间堆栈方法跟踪模拟标记粒子的运动轨迹；

3.1)读入t₀时刻的气象参数，所述气象参数包括边界层厚度H_b、下垫面粗糙度r₀和平均风速V；根据初始时刻t₀标记粒子的位置坐标，插值计算标记粒子所处位置的气象参数，并用大气科学和行星边界层理论计算边界层湍流扩散参数σ与拉格朗日时间尺度T_L；

3.2)计算标记粒子t₁时刻位置和污染物浓度值

3.2.1)根据大气输运和湍流扩散，分别计算初始时刻t₀的所有标记粒子经过时间步Δt后的t₁时刻位置X_j(t₁)，并将t₁时刻每一个标记粒子信息存储在堆栈中；

其中，t₁时刻位置X_j(t₁)计算公式如下：

X_j(t₁)＝X_j(t₀)+V(X_j(t₀),t₀)Δt+V′(X_j(t₀),t₀)Δt

其中，X为标记粒子的坐标，三维分量分别为x,y,z；

V为平均风速，三维分量分别为u,v,w；

V′为湍流脉动速度，三维分量分别为u′,v′,w′；

湍流脉动速度V′遵从马尔科夫假定，即

式中，ξ为符合正态分布的随机数；R_u，R_v，R_w为u，v，w三个速度分量风向的拉格朗日自相关函数，其计算公式如下：

其中，σ_u、σ_v、σ_w为x、y、z三个方向的湍流扩散参数，T_Lu、T_Lv、T_Lw为x、y、z三个方向的拉格朗日时间尺度，σ_u、σ_v、σ_w和T_Lu，T_Lv，T_Lw可通过坐标变换从σ₁，σ₂，σ₃和T_L1，T_L2，T_L3计算得到；

3.2.2)计算t₁时刻每个空间网格上的污染物浓度值，根据每个空间网格上的污染物浓度值，确定并输出t₁时刻污染物的浓度空间分布；

污染物浓度值计算公式如下：

式中，Q为释放污染物总量，i代表记录网格编号，ΔV代表网格体积，T_i,j为该时间步内第j个粒子在第i号记录网格中的停留时间；

3.3)读入t₁时刻的气象参数，所述气象参数包括边界层厚度H_b、下垫面粗糙度r₀和平均风速V；根据t₁时刻标记粒子的位置坐标，插值计算标记粒子所处位置的气象参数，并用大气科学和行星边界层理论计算边界层湍流扩散参数σ与拉格朗日时间尺度T_L；

3.4)以标记粒子t₁时刻位置X_j(t₁)为下一时间步Δt的初始位置，利用步骤3.2)的计算方法，计算t₁时刻的所有标记粒子经过时间步Δt后的t₂时刻位置X_j(t₂)和t₂时刻每个空间网格上的污染物浓度值，并将t₂时刻每一个标记粒子信息存储在堆栈中，以及输出t₂时刻污染物的浓度空间分布；

3.5)重复利用与步骤3.3)和3.4)同样的计算过程，直至跟踪时间达到截断时间T_CUT或者计算出的标记粒子t_k时刻位置X_j(t_k)逸出计算边界；

其中，步骤3.1)和步骤3.3)中，用大气科学和行星边界层理论计算边界层湍流扩散参数σ与拉格朗日时间尺度T_L具体为：

设大气稳定度为中性，则σ与T_L在顺风方向的分量为σ₁和T_L1、水平侧风方向的分量为σ₂和T_L2、垂直方向的分量σ₃和T_L3，计算公式如下：

σ₁＝2.3×u^*

σ₂＝2.0×u^*

其中，z为标记粒子离地高度；u^*为摩擦速度，

式中，κ为卡曼常数，u₀为地面风速，z₀为对应测风高度；

4)得到污染物扩散过程

根据不同时刻所有网格点的空间浓度分布，得到污染物整个扩散过程中浓度分布的时空图像。

2.根据权利要求1所述采用时间堆栈实现空气污染物大气扩散快速预测的方法，其特征在于：步骤3.2.1)和步骤3.4)中，标记粒子信息包括位置坐标和速度信息。

Images