CN100541169C - 扩散物质的扩散状况预测方法及扩散状况预测系统 - Google Patents

Abstract

预测排往大气中的扩散物质的扩散状况。将扩散物质换算为粒子，预先通过计算求出排出源产生的粒子的移动位置，使移动位置和排出后经过时间相对应并记录下来。此外预先设定表示沿排出后经过时间的粒子强度的强度数据。接着根据排出后经过时间，参照上述强度数据，求出粒子发生时刻的强度，使每一粒子的移动位置、排出后经过时间及强度相对应，并记录下来。要想计算规定时间内的规定区域的浓度，可通过累计计算求出该规定区域内存在的粒子强度。

Description

扩散物质的扩散状况预测方法及扩散状况预测系统

技术领域

本发明涉及扩散物质的扩散状况预测方法。本发明是通过预测从扩散源(例如使用放射性物质的设施及烟囱)排往大气中的物质(例如放射性物质及气体)怎样在大气中扩散来预测各个地点时刻变化的物质浓度的。

背景技术

现已开发出放射性物质因事故从处理放射性物质的设施泄漏到外部时，预测放射性物质的扩散范围及各地点的放射性物质浓度，预测有可能遭受放射性物质带来的危险的地区的扩散状况预测方法。

该扩散状况预测方法不仅可适用于预测放射性物质的扩散状况，还可在工厂的烟囱中排出的气体(烟)在大气中扩散的情况下，适用于计算各个地点的气体浓度，以及在环境评估的分析中，适用于分析扩散物质的扩散状况。

要想通过计算预测排往大气中的物质的扩散状况，就得进行以下两种计算。

(1)气体状况预测计算，

(2)扩散状况预测计算。

上述(1)的气体状况预测计算是指根据气象GPV(Grid PointValue)数据及AMEDAS等气象观测数据，通过计算用于分析大气现象的偏微分方程式，求出从事实与现象发生(例如放射性物质泄露到外部)的时刻到规定的目标时刻之间一定时间段内的每一时刻的多数评价地点(方格点位置)的风向、风速，也就是指求出表示一定时段内每一时刻的风速场数据的气体状况的计算。

此外，上述(2)的扩散状况预测计算是指通过将外泄的扩散物质的浓度及性状及上述风速场数据代入用于计算物质(粒子)的扩散状态的扩散方程式，求出每个时间点的方格点位置上的扩散物质浓度的计算。

<气体状况预测计算介绍>

首先介绍气体状况预测计算的概况。气象观测数据、例如气象GPV数据由气象业务支援中心每12小时发送一次。该气象GPV数据表示在地球表面沿南北方向沿伸的同时，东西方向相互间隔距离为规定距离(2Km)的多条经度假想线和在地球表面沿东西方向沿伸的同时，南北方向相互间隔距离为规定距离(2Km)的复数条纬度假想线交叉的地点(将此称之为大方格点位置)的气象数据(包括风速矢量(风向、风速)、气压、温度、湿度)。而且气象GPV数据作为各大方格点位置的气象数据，统一发送发送时刻、以及从发送时刻起未来3小时、未来6小时、未来9小时等每次间隔3小时共计51小时的数据。

由于上述气象GPV数据的大方格点位置的气象数据从空间角度而言，大方格点位置相互间相距2Km之广，而且从时间角度而言彼此相隔3小时之久，因而仅靠该大方格点位置的气象数据所示的气体状况(风向、风速)数据，即风速场数据，无法计算扩散物质的扩散浓度。

因此有必要根据空间性与时间性均粗略的气象观测数据，通过计算用于分析大气现象的偏微分方程式，求出从空间方面及时间方面均更为细密的气体状况(风向、风速等)。

因此在待计算的计算区(在地球表面之内预先设定的特定区域)设定的大方格点位置之间设定小方格点位置。小方格点位置配置于在地球表面沿南北方向沿伸的同时，东西方向的相互间隔距离为规定距离(50m)的复数条经度假想线和在地球表面沿东西方向沿伸的同时，南北方向的相互间距为规定距离(50m)的复数条纬度假想线交叉的地点。

并且通过差分分析计算用于分析大气现象的偏微分方程式，求出从计算开始每隔一定时间段(例如每隔20秒种)的小方格点位置及大方格点位置的气象数据。作为分析大气现象的偏微分方程式，可使用科罗拉多州立大学和Mission Reseavch公司开发的，用RAMS(Regional Atmosphevie Modeling System)码表示的分析风速场的基本方程式。

该用RAMS码表示的分析风速场的基本方程式是由运动方程式、热能方程式、水分的扩散方程式及连续式构成，用下面的式(1)～式(6)表示。

(数学式1)

运动方程式：

$\frac{\&PartialD;u}{\&PartialD;t}=-u\frac{\&PartialD;u}{\&PartialD;x}-v\frac{\&PartialD;u}{\&PartialD;y}-w\frac{\&PartialD;u}{\&PartialD;z}-\mathrm{\&theta;}\frac{\&PartialD;{\mathrm{\&pi;}}^{\&prime;}}{\&PartialD;x}+\mathrm{fv}+\frac{\&PartialD;}{\&PartialD;x}\left(K_m\frac{\&PartialD;u}{\&PartialD;x}\right)+\frac{\&PartialD;}{\&PartialD;y}\left(K_m\frac{\&PartialD;u}{\&PartialD;y}\right)+\frac{\&PartialD;}{\&PartialD;z}\left(K_m\frac{\&PartialD;u}{\&PartialD;z}\right)---\left(1\right)$

$\frac{\&PartialD;v}{\&PartialD;t}=-u\frac{\&PartialD;v}{\&PartialD;x}-v\frac{\&PartialD;u}{\&PartialD;y}-w\frac{\&PartialD;v}{\&PartialD;z}-\mathrm{\&theta;}\frac{\&PartialD;{\mathrm{\&pi;}}^{\&prime;}}{\&PartialD;y}+\mathrm{fu}+\frac{\&PartialD;}{\&PartialD;x}\left(K_m\frac{\&PartialD;v}{\&PartialD;x}\right)+\frac{\&PartialD;}{\&PartialD;y}\left(K_m\frac{\&PartialD;v}{\&PartialD;y}\right)+\frac{\&PartialD;}{\&PartialD;z}\left(K_m\frac{\&PartialD;v}{\&PartialD;z}\right)---\left(2\right)$

$\frac{\&PartialD;w}{\&PartialD;t}=-u\frac{\&PartialD;w}{\&PartialD;x}-v\frac{\&PartialD;w}{\&PartialD;y}-w\frac{\&PartialD;w}{\&PartialD;z}-\mathrm{\&theta;}\frac{\&PartialD;{\mathrm{\&pi;}}^{\&prime;}}{\&PartialD;z}-\frac{g{\mathrm{\&theta;}}_v^{\&prime;}}{{\mathrm{\&theta;}}_0}+\frac{\&PartialD;}{\&PartialD;x}\left(K_m\frac{\&PartialD;w}{\&PartialD;x}\right)+\frac{\&PartialD;}{\&PartialD;y}\left(K_m\frac{\&PartialD;w}{\&PartialD;y}\right)+\frac{\&PartialD;}{\&PartialD;z}\left(K_m\frac{\&PartialD;w}{\&PartialD;z}\right)---\left(3\right)$

热能方程式：

$\frac{\&PartialD;{\mathrm{\&theta;}}_{\mathrm{il}}}{\&PartialD;t}=-u\frac{\&PartialD;{\mathrm{\&theta;}}_{\mathrm{il}}}{\&PartialD;x}-v\frac{\&PartialD;{\mathrm{\&theta;}}_{\mathrm{il}}}{\&PartialD;y}-w\frac{\&PartialD;{\mathrm{\&theta;}}_{\mathrm{il}}}{\&PartialD;z}+\frac{\&PartialD;}{\&PartialD;x}\left(K_h\frac{\&PartialD;{\mathrm{\&theta;}}_{\mathrm{il}}}{\&PartialD;x}\right)+\frac{\&PartialD;}{\&PartialD;y}\left(K_h\frac{\&PartialD;{\mathrm{\&theta;}}_{\mathrm{il}}}{\&PartialD;y}\right)+\frac{\&PartialD;}{\&PartialD;z}\left(K_h\frac{\&PartialD;{\mathrm{\&theta;}}_{\mathrm{il}}}{\&PartialD;z}\right)+{\left(\frac{\&PartialD;{\mathrm{\&theta;}}_{\mathrm{il}}}{\&PartialD;t}\right)}_{\mathrm{rad}}---\left(4\right)$

水分扩散方程式：

$\frac{\&PartialD;r_n}{\&PartialD;t}=-u\frac{\&PartialD;r_n}{\&PartialD;x}-v\frac{\&PartialD;r_n}{\&PartialD;y}-w\frac{\&PartialD;r_n}{\&PartialD;z}+\frac{\&PartialD;}{\&PartialD;x}\left(K_h\frac{\&PartialD;r_n}{\&PartialD;x}\right)+\frac{\&PartialD;}{\&PartialD;y}\left(K_h\frac{\&PartialD;r_n}{\&PartialD;y}\right)+\frac{\&PartialD;}{\&PartialD;z}\left(K_h\frac{\&PartialD;r_n}{\&PartialD;z}\right)---\left(5\right)$

连续式：

$\frac{\&PartialD;{\mathrm{\&pi;}}^{\&prime;}}{\&PartialD;t}=\frac{R{\mathrm{\&pi;}}_0}{C_v{\mathrm{\&rho;}}_0{\mathrm{\&theta;}}_0}(\frac{\&PartialD;{\mathrm{\&rho;}}_0{\mathrm{\&theta;}}_{0}u}{\&PartialD;x}+\frac{\&PartialD;{\mathrm{\&rho;}}_0{\mathrm{\&theta;}}_{0}v}{\&PartialD;y}+\frac{\&PartialD;{\mathrm{\&rho;}}_0{\mathrm{\&theta;}}_{0}w}{\&PartialD;z})---\left(6\right)$

其中，u、v、w：风速，f：哥氏参数，Km：运动量的湍流粘性系数，Kn：热量与水分的湍流扩散系数，θu：水分(冰→水)的温度等级，rn：雨、露等水分的混合比，ρ：密度，rad：辐射能(radiafion)，g：重力加速度，π′：Eoner fancfion(变动部分)，θr：临时温度等级，P：压力，下标字0为参照值，R：气体常数，Cr：定容比热。

如上所示，计算用RAMS(Regional Atmospheric ModelingSystem)码表示的分析风速场的基本方程式，即可取得从计算开始起每一规定时刻(例如每隔20秒)的，表示各大方格点位置的气象数据和各个小方格点位置的气象数据的风向矢量数据(风速场数据)。

<扩散状况预测计算的简要介绍>

下面介绍扩散状况预测计算。要想进行扩散状况预测计算，就得把用RAMS码求出的每隔20秒的各个大方格点位置及各个小方格点位置的风速场数据不断代入科罗拉多州立大学和Mission Research公司开发的HYPACT(Hybrid Particle Conlenfrafion Transport Model)码进行扩散状况预测计算。作为扩散状况的预测计算实例，采用了拉格朗日(Lagrangian)粒子扩散模型。

在该拉格朗日粒子扩散模型之中，采用下面所示的式(7)～(9)计算粒子的扩散速度(u′，v′，w′)，使各粒子移动。

(数学式2)

拉格朗日粒子模型采用式(12)～(14)计算粒子的扩散速度。

$u^{\&prime;}\left(t\right)=R_u{u}^{\&prime;}(t-\mathrm{\&Delta;t})+(1-R_u^2)r_u$

$v^{\&prime;}\left(t\right)=R_v{v}^{\&prime;}(t-\mathrm{\&Delta;t})+(1-R_v^2)r_v---\left(7\right)$

$w^{\&prime;}\left(t\right)=R_w{w}^{\&prime;}(t-\mathrm{\&Delta;t})+(1-R_w^2)r_w$

式中，Ru、Rv、Rw：拉格朗目湍流自相关函数，

u′(t)、v′(t)、w′(t)：粒子的湍流扩散速度成分，t：时间。

$R_u\left(\mathrm{\&Delta;t}\right)=\frac{\stackrel{\&OverBar;}{u^{\&prime;}\left(t\right)\&CenterDot;u^{\&prime;}(t-\mathrm{\&Delta;t})}}{{\mathrm{\&sigma;}}_u^2}=\exp (-\frac{\mathrm{\&Delta;t}}{T_{\mathrm{Lu}}})$

$R_v\left(\mathrm{\&Delta;t}\right)=\frac{\stackrel{\&OverBar;}{v^{\&prime;}\left(t\right)\&CenterDot;v^{\&prime;}(t-\mathrm{\&Delta;t})}}{{\mathrm{\&sigma;}}_v^2}=\exp (-\frac{\mathrm{\&Delta;t}}{T_{\mathrm{Lv}}})---\left(8\right)$

$R_w\left(\mathrm{\&Delta;t}\right)=\frac{\stackrel{\&OverBar;}{w^{\&prime;}\left(t\right)\&CenterDot;w^{\&prime;}(t-\mathrm{\&Delta;t})}}{{\mathrm{\&sigma;}}_w^2}=\exp (-\frac{\mathrm{\&Delta;t}}{T_{\mathrm{Lw}}})$

式中σu、σv、σw：湍流速度标准偏差，

TLu、TLv、TLw：拉格朗日时间刻度。

r_u＝σ_uη_u    r_v＝σ_vη_v    r_w＝σ_wη_w+w_d  (9)

式中：ηu、ηv、ηw：正规随机抽样数(平均值为零)

wd：重力沉降速度。

这里介绍将通过RAMS码求出的每隔20秒钟的各个大方格点位置及各个小方格点位置的风速场数据陆续地代入HYPACT(HybridParticle Concentration Transport Model)码，进行扩散状况的预测计算的具体例。

为了进行该计算，将从排出源排往大气中的物质换算为多数粒子，设定为每一计算周期Δt(此处的Δt＝20秒钟)从排出源的位置产生N个(此处为20个)粒子P。

也就是说，在计算开始的时刻，使之产生20个粒子P，从计算开始时刻起20秒后再产生20个粒子，从计算开始时刻起40秒后再产生20个粒子，以此类推，在每个计算周期Δt(20秒)内使之产生20个粒子。并且每隔计算周期Δt(20秒)通过计算求出各粒子P的位置(空间座标)。

并用以下方式表示计算开始时刻(时间0秒)产生的20个粒子P：

P00⁰¹、P00⁰²、P00⁰³、P00⁰⁴、P00⁰⁵、P00⁰⁶、P00⁰⁷、P00⁰⁸、P00⁰⁹、P00¹⁰、P00¹¹、P00¹²、P00¹³、P00¹⁴、P00¹⁵、P00¹⁶、P00¹⁷、P00¹⁸、P00¹⁹、P00²⁰。

从计算开始时刻起20秒后产生的20个粒子P表示为：

P20⁰¹、P20⁰²、P20⁰³、P20⁰⁴、P20⁰⁵、P20⁰⁶、P20⁰⁷、P20⁰⁸、P20⁰⁹、P20¹⁰、P20¹¹、P20¹²、P20¹³、P20¹⁴、P20¹⁵、P20¹⁶、P20¹⁷、P20¹⁸、P20¹⁹、P20²⁰。

从计算开始时刻起40秒后发生的20个粒子P表示为：

P40⁰¹、P40⁰²、P40⁰³、P40⁰⁴、P40⁰⁵、P40⁰⁶、P40⁰⁷、P40⁰⁸、P40⁰⁹、P40¹⁰、P40¹¹、P40¹²、P40¹³、P40¹⁴、P40¹⁵、P40¹⁶、P40¹⁷、P40¹⁸、P40¹⁹、P40²⁰。

也就是说在标号“P”后下标所示的数字为从计算开始时刻起的时间，而在标号“P”后上标所示的数字则是用来区别该时刻内产生的20个粒子。在其它时刻产生的粒子也采用同样的方法标示。

首先，在计算开始时刻从排出源产生20个粒子：P00⁰¹、P00⁰²、P00⁰³、P00⁰⁴、P00⁰⁵、P00⁰⁶、P00⁰⁷、P00⁰⁸、P00⁰⁹、P00¹⁰、P00¹¹、P00¹²、P00¹³、P00¹⁴、P00¹⁵、P00¹⁶、P00¹⁷、P00¹⁸、P00¹⁹、P00²⁰。

从计算开始时刻起20秒后从图19所示的排出源S新产生了20个粒子：P20⁰¹、P20⁰²、P20⁰³、P20⁰⁴、P20⁰⁵、P20⁰⁶、P20⁰⁷、P20⁰⁸、P20⁰⁹、P20¹⁰、P20¹¹、P20¹²、P20¹³、P20¹⁴、P20¹⁵、P20¹⁶、P20¹⁷、P20¹⁸、P20¹⁹、P20²⁰。

此时计算开始时刻产生的粒子P00⁰¹、P00⁰²、P00⁰³、P00⁰⁴、P00⁰⁵、P00⁰⁶、P00⁰⁷、P00⁰⁸、P00⁰⁹、P00¹⁰、P00¹¹、P00¹²、P00¹³、P00¹⁴、P00¹⁵、P00¹⁶、P00¹⁷、P00¹⁸、P00¹⁹、P00²⁰到达远离排出源的位置，同时进行扩散。

各粒子P的位置是采用通过RAMS码求出的每隔20秒的风速场数据，计算拉格朗日粒子扩散模型中的各粒子P的扩散速度(u′、v′、w′)，通过使各粒子的移动求出的。

从计算开始时刻40秒后，由图20所示的排出源新产生了20个粒子：P40⁰¹、P40⁰²、P40⁰³、P40⁰⁴、P40⁰⁵、P40⁰⁶、P40⁰⁷、P40⁰⁸、P40⁰⁹、P40¹⁰、P40¹¹、P40¹²、P40¹³、P40¹⁴、P40¹⁵、P40¹⁶、P40¹⁷、P40¹⁸、P40¹⁹、P40²⁰。

此时计算开始时刻产生的粒子：P00⁰¹、P00⁰²、P00⁰³、P00⁰⁴、P00⁰⁵、P00⁰⁶、P00⁰⁷、P00⁰⁸、P00⁰⁹、P00¹⁰、P00¹¹、P00¹²、P00¹³、P00¹⁴、P00¹⁵、P00¹⁶、P00¹⁷、P00¹⁸、P00¹⁹、P00²⁰到达距排出源更远的位置，同时进一步扩散。

此外，从计算开始时刻起20秒后产生的20个粒子：P20⁰¹、P20⁰²、P20⁰³、P20⁰⁴、P20⁰⁵、P20⁰⁶、P20⁰⁷、P20⁰⁸、P20⁰⁹、P20¹⁰、P20¹¹、P20¹²、P20¹³、P20¹⁴、P20¹⁵、P20¹⁶、P20¹⁷、P20¹⁸、P20¹⁹、P20²⁰到达远离排出源的位置的同时进行扩散。

各粒子P的位置是采用通过RAMS码求出的每隔20秒的风速场数据，计算拉格朗日粒子扩散模型中的各粒子P的扩散速度(u′v′w′)，通过使各粒子移动求出的。

从计算开始时刻起60秒后，由图21所示的排出源新产生20个粒子：P60⁰¹、P60⁰²、P60⁰³、P60⁰⁴、P60⁰⁵、P60⁰⁶、P60⁰⁷、P60⁰⁸、P60⁰⁹、P60¹⁰、P60¹¹、P60¹²、P60¹³、P60¹⁴、P60¹⁵、P60¹⁶、P60¹⁷、P60¹⁸、P60¹⁹、P60²⁰。

此时，计算开始时刻产生的粒子P00⁰¹、P00⁰²、P00⁰³、P00⁰⁴、P00⁰⁵、P00⁰⁶、P00⁰⁷、P00⁰⁸、P00⁰⁹、P00¹⁰、P00¹¹、P00¹²、P00¹³、P00¹⁴、P00¹⁵、P00¹⁶、P00¹⁷、P00¹⁸、P00¹⁹、P00²⁰到达距排出源更远的位置，同时进一步扩散。

此外，从计算开始时刻起20秒后产生的20个粒子P20⁰¹、P20⁰²、P20⁰³、P20⁰⁴、P20⁰⁵、P20⁰⁶、P20⁰⁷、P20⁰⁸、P20⁰⁹、P20¹⁰、P20¹¹、P20¹²、P20¹³、P20¹⁴、P20¹⁵、P20¹⁶、P20¹⁷、P20¹⁸、P20¹⁹、P20²⁰到达距排出源更远的位置同时进一步扩散。

此外，从计算开始时刻起40秒后产生的20个粒子P40⁰¹、P40⁰²、P40⁰³、P40⁰⁴、P40⁰⁵、P40⁰⁶、P40⁰⁷、P40⁰⁸、P40⁰⁹、P40¹⁰、P40¹¹、P40¹²、P40¹³、P40¹⁴、P40¹⁵、P40¹⁶、P40¹⁷、P40¹⁸、P40¹⁹、P40²⁰到达远离排出源的位置，同时进行扩散。

各粒子P的位置是采用通过RAMS码求出的每隔20秒的风速场数据，计算拉格朗日粒子扩散模型中的各粒子P的扩散速度(u′、v′、w′)，通过使各粒子移动求出的。

正如上述，每隔一个计算周期Δt(20秒)使之接连产生20个粒子的同时，接连求出每个计算周期Δt(20秒)内的粒子位置，即空间座标(xi(t)、yi(t)zi(t))。

而在从计算开始到经过规定时间时，正如图22所示，在距排出源规定距离的单位空间(预测地区的单位体积)内存在粒子P的情况下，即可根据该粒子的数量计算出该单位空间内的物质浓度。

也就是说，设在该排出源S每秒排出Q(m³)物质，由于粒子P在20秒钟内产生20个(如进行换算则为每秒钟一个)，则各个粒子P每个粒子具有Q/1(m³)的排出源强度。因此通过将该单位空间内存在的粒子P的数量乘以排出源强度Q/1(m³)即可求出单位空间内的物质浓度。

如果将上述具体例一般性表示，则如下所示。将排出源排出的气体等物质换算为多数粒子，而且从排出源每秒排出N个粒子。这种情况下，计算上的粒子排出量为N/秒。实际由排出源排出的物质的排出量为Q(m³/秒)时，各粒子成为具有Q/N(m³)的排出源强度。

通过针对每个粒子进行非稳定性数值计算运动方程式，也就是说将通过RAMS码求出的风速场数据代入粒子运动方程式式HYPACT码，采用拉格朗日粒子扩散模型计算各粒子P的扩数速度(u′、v′、w′)，通过使各粒子的移动，即可非稳定性地决定各粒子的座标。也就是可决定每一计算周期Δt各粒子的空间座标。而通过拉格朗日粒子扩散模型求出后，记录在记录装置中的各粒子的数据仅为各粒子的空间座标(xi(t)、yi(t)、zi(t))。

粒子(物质)的运动方程式HYPACT码是用来表现粒子的平流、扩散、重力沉降现象的。在此，粒子的平流现象依赖于大气的时间平均速度，扩散现象依赖于大气的湍流速度，重力沉降依赖于粒子的质量、重力加速度、空气的粘度系数等因素(参照图23)。

当空气中的单位体积内的粒子个数为n个时，则该空间内的气体浓度(物质浓度)为n×Q/N(气体m³/空气m³)。也就是该单位空间内存在的粒子数n乘以各粒子具有的排出源强度Q/N。

发明内容

该环境浓度(单位体积内的物质浓度)依赖于排出物质的排出量的经时变化。所以在排出量与时间同为变化的条件的情况下，扩散计算必须针对每个排出条件实施。因此，当设想的排出条件很多的情况下，就得针对每例排出量进行扩散计算，其结果是需要极其大量的计算时间。

也就是说，正如图24所示，例如当从排出源S(例如烟囱)中排出气体时，处于下风地点F的气体浓度的经时变化随排出源S排出的物质的经时变化而变化。

也就是说，如图25(a)所示，当物质的排出量与时间共同变化时，地点F的物质浓度如图25(b)所示，随时间一同变化，当如图26(a)所示，物质的排出量固定的情况下，地点F的物质浓度如图26(b)所示，上升到一定值之后即维护一定的浓度，当如图27(a)所示，物质为瞬间排出的情况下，地点F的物质浓度如图27(b)所示，短暂上升之后变为零。

象这样，物质的排出量随时间一同变化的情况下，有必要使粒子的产生个数与物质的排出量结合，使之随时间一同变化。并且象这样求出与时间的经过一同改变其产生个数的粒子的移动位置，从该粒子的移动位置计算出物质的浓度。因此，必须针对排出量的变化不同的各种物质进行扩散计算，从而需要相当庞大的计算结果。

例如，当处理放射性物质的设施发生了放射物质泄露事故的情况下，可排出种类极多的的物质(例如100种左右的物质)。而且每种物质的排出量又因时间而各自不同。因此需要针对每种物质，使粒子的发生个数与物质的排出量结合，随时间一同变化，这样求出发生个数改变的粒子的移动位置，再从该粒子的移动位置计算出物质的浓度。因此在这种情况下，必须进行与100种物质相对应的100种扩散计算。

本发明的目的在于鉴于以上现有技术，提供一种扩散物质的扩散状况预测方法及扩散物质的扩散状况预测系统。该方法及系统即使在排出多种物质，同时每种物质的排出量随时间一道变化的情况下，仍然可在短时间内预测计算出物质的扩散状况。

解决上述课题的本发明的扩散物质的扩散状况预测方法的特征在于：为了预测从排出源排往大气中的物质在大气中的扩散状况，将上述物质换算为多数粒子，设定为每一计算周期从排出源位置产生预先设定个数的粒子；同时通过将包括排出源位置在内的区域内的多数地点的表示沿时间过程变化的风向、风速的风速场数据代八计算粒子扩散状态的扩散方程式，求出各粒子的扩散速度，根据该扩散速度，求出表示各粒子在每一计算周期内所在的空间位置的空间座标的同时，计测上述粒子最初产生的时刻起的经过时间——排出后经过时间，使各粒子在各计算周期内的空间座标和各粒子的排出后经过时间相对应，记录到记录装置之中；按照排出的物质排出后经过时间的时间过程所伴随的排出量的变化比例，预先设定与沿排出后经过时间的时间过程的粒子相对应的排出源强度数据；读出记录在上述数据记录装置中的，各粒子在每一计算周期内的空间座标和各粒子的排出后经过时间，同时参照读出的排出后经过时间求出各粒子产生的时刻，根据上述排出源强度数据求出该时刻的各粒子的排出源强度，使排出源强度与各粒子在每一计算周期内的空间座标和各粒子的排出后经过时间相对应，重新记录到上述数据记录装置之中；规定的计算周期内的规定区域的上述物质浓度，通过累计计算该规定的计算周期内，存在于该规定区域内的所有粒子的排出源强度求出。

此外，本发明的扩散物质的扩散状况预测方法的特征还在于：为了预测从多个排出源排往大气中的物质在大气中的扩散状况，将上述物质换算为多数粒子，设定为从各个排出源位置，在每一计算周期内各自产生预先设定个数的粒子；同时通过将包括排出源位置在内的区域内的多数地点的表示沿时间过程变化的风向、风速的风速场数据代入计算粒子扩散状态的扩散方程式，求出各粒子的扩散速度，根据该扩散速度，求出每一计算周期内表示各粒子所在的空间位置的空间座标，同时计测上述粒子最初产生的时刻起的经过时间——排出后经过时间，使各计算周期内的各粒子的空间座标和各粒子的排出后经过时间以及识别排出源的排出源识别信息相对应，并记录到数据装置之中。

按照从各排出的物质的排出源排出后经过时间的时间过程相伴随的排出量的变化比例，预先分别设定与沿排出后经过时间的时间过程的粒子相对应的排出源强度数据；读出记录在上述记录装置中的，各粒子在每一计算周期内的空间座标和各粒子排出后经过时间及各粒子的排出源识别信息，同时参照读出的排出后经过时间，求出各粒子的产生时刻，参照读出的排出源识别信息，根据与产生该粒子的排出源相对应的上述排出源强度数据，求出粒子产生时刻的各粒子的排出源强度，使各粒子在每一计算周期的空间座标和各粒子的排出后经过时间及排出源强度相对应，由上述数据记录装置重新记录；规定的计算周期内的规定区域的上述物质的浓度，通过累计计算该规定的计算周期内的存在于该规定区域内的所有粒子的排出源强度求出。

此外，本发明的扩散物质的扩散状况预测方法的特征在于：既可以通过实测求出从上述排出源实际排出的物质的浓度后设定上述排出源强度数据，也可根据在上述排出源周围的观测点实测的物质浓度的经时变化设定上述排出源强度数据。

此外，本发明的扩散物质的扩散状态预测系统，其特征在于：该系统包括当扩散物质排往大气中时，实测扩散物质的浓度，并发送表示扩散物质的排出量的数据的企业；发送气象观测数据的气象数据发送设施；向上述企业及上述企业周围的居民发布避难通告的政府监察部门；进行扩散物质的扩散状况预测计算处理，计算规定区域内的物质浓度的安全分析中心；上述安全分析中心通过信息传输手段接收上述企业发送的表示扩散物质的排出量的数据，以及通过信息传输手段接收上述气象数据发送设施发送的气象观测数据；上述政府监察部门通过信息传输手段接收上述安全分析中心发送的物质浓度，并根据收到的物质的浓度报告发布避难通告。

附图说明

图1是表示本发明的第1实施方式的计算流程的流程图。

图2是表示本发明的第1实施方式的粒子扩散状态的说明图。

图3是表示本发明的第1实施方式的粒子扩散状态的说明图。

图4是表示本发明的第1实施方式的粒子扩散状态的说明图。

图5是表示物质的排出量的经时变化示例的特性图。

图6是表示与物质的排出量的经时变化相对应的排出源强度示例的特性图。

图7是表示在规定的方格区域内粒子分布的说明图。

图8是表示排出源与方格区域的说明图。

图9是当排出量固定时，表示排出量与浓度关系的特性图。

图10当排出量经时变化时，表示排出量与浓度关系的特性图。

图11是表示两个排出源与方格区域的说明图。

图12是当排出量固定时，表示排出量与浓度关系的特性图。

图13是当排出量经时变化时，表示排出量与浓度关系的特性图。

图14是表示第3实施方式的的说明图。

图15是表示第4实施方式的说明图。

图16是表示第4实施方式中的计算流程的流程图。

图17是表示第5实施方式中的计算流程的流程图。

图18是表示第6实施方式涉及的系统的系统构成图。

图19是表示现有技术中的粒子的扩散状态的说明图。

图20是表示现有技术中的粒子的扩散状态的说明图。

图21是表示现有技术中的粒子的扩散状态的说明图。

图22是表示规定的方格区域内的粒子分布的说明图。

图23是表示粒子扩散模型的功能的说明图。

图24是表示排出源和方格区域的说明图。

图25是当排出量经时变化时，表示排出量与浓度关系的特性图。

图26是当排出量固定时，表示排出量与浓度关系的特性图。

图27是当排出量为瞬间时，表示排出量与浓度关系的特性图。

具体实施方式

下面根据附图，详细介绍本发明的具体实施方式。

<第1实施方式(排出源为一个)>

下面参照图1～图10，介绍本发明的第1实施方式涉及的扩散物质的扩散状况预测方法。

在第1实施方式的第1步骤(参照计算流程图图1)中，实施下述处理。即，无论从排出源S(参照表示粒子扩散状况的图2～图4)实际排出的物质的排出量是固定值，还是沿时间过程排出量有变化，首先将物质的排出量Q(m³/秒)设为固定值(＝1.0)，采用现用的拉格朗日扩散模型，计数粒子的轨迹。进而除各粒子具有的信息——空间座标(xi(t)、yi(t)、zi(t))之外，将粒子最初产生的时刻起的经过时间——排出后经过时间Ti(t)，在每个计算周期内记录到数据记录装置1之中。

对该第1步骤的处理的具体说明。如下述。在该计算中，每一计算周期Δ(t)(此处Δt＝20秒)，产生20个粒子，同时在每一计算周期Δt(20秒)内，计算粒子P的位置(空间座标)。

首先，在计算开始时刻，从排出源S产生20个粒子P00⁰¹、P00⁰²、P00⁰³、P00⁰⁴、P00⁰⁵、P00⁰⁶、P00⁰⁷、P00⁰⁸、P00⁰⁹、P00¹⁰、P00¹¹、P00¹²、P00¹³、P00¹⁴、P00¹⁵、P00¹⁶、P00¹⁷、P00¹⁸、P00¹⁹、P00²⁰。

从计算开始时刻起20秒钟之后，从图2所示的排出源S新产生20个粒子P20⁰¹、P20⁰²、P20⁰³、P20⁰⁴、P20⁰⁵、P20⁰⁶、P20⁰⁷、P20⁰⁸、P20⁰⁹、P20¹⁰、P20¹¹、P20¹²、P20¹³、P20¹⁴、P20¹⁵、P20¹⁶、P20¹⁷、P20¹⁸、P20¹⁹、P20²⁰。

此时，计算开始时刻产生的粒子P00⁰¹、P00⁰²、P00⁰³、P00⁰⁴、P00⁰⁵、P00⁰⁶、P00⁰⁷、P00⁰⁸、P00⁰⁹、P00¹⁰、P00¹¹、P00¹²、P00¹³、P00¹⁴、P00¹⁵、P00¹⁶、P00¹⁷、P00¹⁸、P00¹⁹、P00²⁰已到达远离排出源S的位置，同时在扩散。

各粒子P00⁰¹～P00²⁰的位置，是采用通过RAMS码求出的每20秒的风速场数据，计算拉格朗日粒子扩散模型中的各粒子P00⁰¹～P00²⁰的扩散速度(u′、v′、w′)，通过使各粒子移动求出的。

还有，计算开始时刻产生的粒子P00⁰¹～P00²⁰，从计算开始时刻(粒子最初产生的时刻)起，经过了20秒钟。因此分别使排出后经过时间Ti(t)＝20对应于各粒子P00⁰¹～P00²⁰的各空间座标(xi(t＝20)、yi(t＝20)、zi(t＝20))之后，记录到记录装置1之中(参照图1、图2)。

从计算开始时刻起40秒种之后，从图3所示的排出源S新产生20个粒子P40⁰¹、P40⁰²、P40⁰³、P40⁰⁴、P40⁰⁵、P40⁰⁶、P40⁰⁷、P40⁰⁸、P40⁰⁹、P40¹⁰、P40¹¹、P40¹²、P40¹³、P40¹⁴、P40¹⁵、P40¹⁶、P40¹⁷、P40¹⁸、P40¹⁹、P40²⁰。

此时计算开始时刻产生的粒子P00⁰¹、P00⁰²、P00⁰³、P00⁰⁴、P00⁰⁵、P00⁰⁶、P00⁰⁷、P00⁰⁸、P00⁰⁹、P00¹⁰、P00¹¹、P00¹²、P00¹³、P00¹⁴、P00¹⁵、P00¹⁶、P00¹⁷、P00¹⁸、P00¹⁹、P00²⁰到达距排出源更远的位置，同时进一步扩散。

此外，从计算开始时刻起20秒钟后产生的20个粒子P20⁰¹、P20⁰²、P20⁰³、P20⁰⁴、P20⁰⁵、P20⁰⁶、P20⁰⁷、P20⁰⁸、P20⁰⁹、P20¹⁰、P20¹¹、P20¹²、P20¹³、P20¹⁴、P20¹⁵、P20¹⁶、P20¹⁷、P20¹⁸、P20¹⁹、P20²⁰到达远离排出源的位置，同时在扩散。

各粒子P00⁰¹～P00²⁰、P20⁰¹～P20²⁰的位置，是采用通过RAMS码求出的每隔20秒钟的风速场数据，计算拉格朗日粒子扩散模型中的各粒子P00⁰¹～P00²⁰、P20⁰¹～P20²⁰的扩散速度(u′、v′、w′)，通过使各粒子移动求出的。

还有，计算开始时刻产生的粒子P00⁰¹～P00²⁰，从计算始时刻(粒子最初产生的时刻)起，已过了40秒钟。因此使排出后经过时间Ti(t)＝40秒分别对应于各粒子P00⁰¹～P00²⁰的各空间座标(xi(t＝40)、yi(t＝40)、zi(t＝40))之后记录到数据记录装置1之中(参照图1、图3)。

此外，从计算开始时刻起20秒后产生的粒子P20⁰¹～P20²⁰从计算开始时刻(粒子最初产生的时刻)起已经过了20秒钟。因此使排出后经过时间下Ti(t)＝20秒分别对应于各粒子P20⁰¹～P20²⁰的空间座标(xi(t＝40)、yi(t＝40)、zi(t＝40))之后，记录到记录装置1之中(参照图1、图3)。

从计算开始时刻起60秒钟后，由图4所示的排出源S新产生20个粒子P60⁰¹、P60⁰²、P60⁰³、P60⁰⁴、P60⁰⁵、P60⁰⁶、P60⁰⁷、P60⁰⁸、P60⁰⁹、P60¹⁰、P60¹¹、P60¹²、P60¹³、P60¹⁴、P60¹⁵、P60¹⁶、P60¹⁷、P60¹⁸、P60¹⁹、P60²⁰。

此时，计算开始时刻产生的粒子P00⁰¹、P00⁰²、P00⁰³、P00⁰⁴、P00⁰⁵、P00⁰⁶、P00⁰⁷、P00⁰⁸、P00⁰⁹、P00¹⁰、P00¹¹、P00¹²、P00¹³、P00¹⁴、P00¹⁵、P00¹⁶、P00¹⁷、P00¹⁸、P00¹⁹、P00²⁰到达距排出源更远的位置，同时进一步扩散。

此外，从计算开始时刻起20秒后产生的20个粒子P20⁰¹、P20⁰²、P20⁰³、P20⁰⁴、P20⁰⁵、P20⁰⁶、P20⁰⁷、P20⁰⁸、P20⁰⁹、P20¹⁰、P20¹¹、P20¹²、P20¹³、P20¹⁴、P20¹⁵、P20¹⁶、P20¹⁷、P20¹⁸、P20¹⁹、P20²⁰到达距排出源更远的位置同时进一步扩散。

此外，从计算开始时刻起40秒后产生的20个粒子P40⁰¹、P40⁰²、P40⁰³、P40⁰⁴、P40⁰⁵、P40⁰⁶、P40⁰⁷、P40⁰⁸、P40⁰⁹、P40¹⁰、P40¹¹、P40¹²、P40¹³、P40¹⁴、P40¹⁵、P40¹⁶、P40¹⁷、P40¹⁸、P40¹⁹、P40²⁰到达远离排出源的位置，同时在扩散。

各粒子P00⁰¹～P00²⁰、P20⁰¹～P20²⁰、P40⁰¹～P40²⁰的位置，是采用通过RAMS码求出的每隔20秒钟的风速场数据，计算拉格朗日粒子扩散模型中的各粒子P00⁰¹～P00²⁰、P20⁰¹～P20²⁰、P40⁰¹～P40²⁰的扩散速度(u′、v′、w′’)，通过使各粒子移动求出的。

还有，计算开始时刻产生的粒子P00⁰¹～P00²⁰、从计算开始时刻(粒子最初产生的时刻)起，已过了60秒钟，因此使排出后经过时间Ti(t)＝60秒分别与各粒子P00⁰¹～P00²⁰的空间座标(xi(t＝60)、yi(t＝60)、zi(t＝60))相对应之后，记录到数据记录装置1之中(参照图1、图4)。

此外，从计算开始时刻起20秒后产生的粒子P20⁰¹～P20²⁰，从计算开始时刻(使粒子最初产生的时刻)起，已过了40秒钟。因此，使排出后经过时间Ti(t)＝40秒分别与各粒子P20⁰¹～P20²⁰的空间座标(xi(t＝60)、yi(t＝60)、zi(t＝60))相对应之后，记录到数据记录装置1之中(参照图1、图4)。

此外，从计算开始时刻起40秒后产生的粒子P40⁰¹～P40²⁰，从计算开始时刻(使粒子最初产生的时刻)起，已经过20秒钟，因此，使排出后经过时间Ti(t)＝20秒分别与各粒子P40⁰¹～P40²⁰的各空间座标(xi(t＝60)、yi(t＝60)、zi(t＝60))相对应之后，记录到记录装置1之中(参照图1、图4)。

如上所述，使之在每个计算周期Δt(20秒)接连产生20个粒子的同时，接连求出在每个计算周期Δt(20秒)内的粒子位置，即空间座标(xi(t)、yi(t)、zi(t))。此外，预先计测出各计算周期Δt内的排出后经过时间Ti(t)，使各计算周期内的各粒子的空间座标和各粒子的排出后经过时间相对应，陆续记录到数据记录装置1之中。

下面具体介绍第2步骤(参照图1)的处理。在上述的第1步骤中将物质的排出量Q(m³/秒)设为固定值(＝0.1)，进行了计数。然而从实际的排出源S排出的物质的排出量正如图5所示，大多数随着排出后经过时间Ti(t)的经过而改变。因此在此类排出量随时间变化的情况下，设定与该图5所示的物质的排出量相对应的，如图6所示的，表示与沿排出后经过时间Ti(t)的时间过程的粒子相对应的排出强度的数据。

图6中的排出源强度数据，例如排出后经过时间Ti(t)为0秒、20秒、60秒时，排出源强度分别为0.3、0.9、0.6。

接着，在读出记录在数据记录装置1中的，每一计算周期的各粒子的空间座标和各粒子的排出后经过时间Ti(t)的同时，参照每一粒子排出后经过时间Ti(t)，求出该粒子的产生时刻，根据图6所示的排出源强度数据，求出该时刻的各粒子的排出源强度。进而在每一计算周期，使各粒子的空间座标和各粒子的排出后经过时间Ti(t)及排出源强度相对应，重新记录到数据记录装置1之中。

具体而言，作为排出后经过时间Ti(t)为20秒(第1个计算周期)时的数据，各粒子P00⁰¹～P00²⁰的各空间座标(xi(t＝20)、yi(t＝20)、zi(t＝20))和排出后经过时间Ti(t)＝20秒相对应，由数据记录装置1记录(参照图2)。

因此，读出该各粒子P00⁰¹～P00²⁰的各空间座标(xi(t＝20)、yi(t＝20)、zi(t＝20))和排出后经过时间Ti(t)＝20，通过从现在时刻t＝20秒减去排出后经过时间Ti＝20秒，求出各粒子P00⁰¹～P00²⁰产生时刻的排出后经过时间Ti(t)＝0秒。并且根据图6所示的排出强度数据，求出粒子P00⁰¹～P00²⁰产生时的排出后经过时间Ti(t)＝0秒时的排出源强度0.3。

并使各粒子P00⁰¹～P00²⁰的各空间座标(xi(t＝20)、yi(t＝20)、zi(t＝20))和排出后经过时间Ti(t)＝20秒、及各粒子P00⁰¹～P00²⁰的排出源强度0.3相对应之后，重新记录到数据记录装置1之中。

此外，作为排出后经过时间Ti(t)为40秒(第2个计算周期)时的数据，

使各粒子P00⁰¹～P00²⁰的各空间座标(xi(t＝40)、yi(t＝40)、zi(t＝40))和排出后经过时间Ti(t)＝40秒相对应，以及

各粒子P20⁰¹～P20²⁰的各空间座标(xi(t＝40)、yi(t＝40)、zi(t＝40))和排出后经过时间Ti(t)＝20秒相对应后，由数据记录装置1记录(参照图3)。

因此，读出各粒子P00⁰¹～P00²⁰的各空间座标(xi(t＝40)、yi(t＝40)、zi(t＝40))和排出后经过时间Ti(t)＝40秒，通过从现在的时刻t＝40秒中减去排出后经过时间Ti(t)＝40秒，求出各粒子P00⁰¹～P00²⁰产生时刻的排出后经过时间Ti(t)＝0秒。并根据图6所示的排出强度数据，求出粒子P00⁰¹～P00²⁰产生时的排出后经过时间Ti(t)＝0秒时的排出源强度0.3。

与此相同，读出各粒子P20⁰¹～P20²⁰的各空间座标(xi(t＝40)、yi(t＝40)、zi(t＝40))和排出后经过时间Ti(t)＝20秒，通过从现在的时刻t＝40秒中减去排出后经过时间Ti(t)＝20秒，求出各粒子P20⁰¹～P20²⁰产生时刻的排出后经过时间Ti(t)＝20秒。并根据图6所示的排出强度数据，求出粒子P20⁰¹～P20²⁰产生时的排出后经过时间Ti(t)＝20秒时的排出源强度0.5。

并且使各粒子P00⁰¹～P00²⁰的各空间座标(xi(t＝40)、yi(t＝40)、zi(t＝40))和排出后经过时间Ti(t)＝40秒及各粒子P00⁰¹～P00²⁰的排出源强度0.3相对应，重新记录到数据记录装置1之中。

此外，使各粒子P20⁰¹～P20²⁰的各空间座标(xi(t＝40)、yi(t＝40)、zi(t＝40))和排出后经过时间Ti(t)＝20秒及各粒子P20⁰¹～P20²⁰的排出源强度0.5相对应，重新记录到数据记录装置1之中。

此外，作为排出后经过时间Ti(t)为60秒(第3次计算周期)时的数据，

使各粒子P00⁰¹～P00²⁰的各空间座标(xi(t＝60)、yi(t＝60)、zi(t＝60))和排出后经过时间Ti(t)＝60秒，

各粒子P20⁰¹～P20²⁰的各空间座标(xi(t＝60)、yi(t＝60)、zi(t＝60))和排出后经过时间Ti(t)＝40秒，以及

各粒子P40⁰¹～P40²⁰的各空间座标(xi(t＝60)、yi(t＝60)、zi(t＝60))和排出后经过时间Ti(t)＝20秒相对应后，由数据记录装置1记录(参照图4)。

因此，读出各粒子P00⁰¹～P00²⁰的各空间座标(xi(t＝60)、yi(t＝60)、zi(t＝60))和排出后经过时间Ti(t)＝60秒，通过从现在的时刻t＝60秒中减去排出后经过时间Ti(t)＝60秒，求出各粒子P00⁰¹～P00²⁰产生时刻的排出后经过时间Ti(t)＝0秒。并根据图6所示的排出强度数据，求出粒子P00⁰¹～P00²⁰产生时的排出后经过时间Ti(t)＝0秒时的排出源强度0.3。

与此相同，读出各粒子P20⁰¹～P20²⁰的各空间座标(xi(t＝60)、yi(t＝60)、zi(t＝60))和排出后经过时间Ti(t)＝60秒，通过从现在的时刻t＝60秒中减去排出后经过时间Ti(t)＝60秒，求出各粒子P20⁰¹～P20²⁰产生时刻的排出后经过时间Ti(t)＝20秒。并根据图6所示的排出强度数据，求出粒子P20⁰¹～P20²⁰产生时的排出后经过时间Ti(t)＝20秒时的排出源强度0.5。

与此相同，读出各粒子P40⁰¹～P40²⁰的各空间座标(xi(t＝60)、yi(t＝60)、zi(t＝60))和排出后经过时间Ti(t)＝20秒，通过从现在的时刻t＝60秒中减去排出后经过时间Ti(t)＝20秒，求出各粒子P40⁰¹～P40²⁰产生时刻的排出后经过时间Ti(t)＝40秒。并根据图6所示的排出强度数据，求出粒子P40⁰¹～P40²⁰产生时的排出后经过时间Ti(t)＝40秒时的排出源强度0.9。

并使各粒子P00⁰¹～P00²⁰的各空间座标(xi(t＝60)、yi(t＝60)、zi(t＝60))和排出后经过时间Ti(t)＝60秒，及各粒子P00⁰¹～P00²⁰的排出源强度0.3相对应之后，重新记录到数据记录装置1之中。

此外，使各粒子P20⁰¹～P20²⁰的各空间座标(xi(t＝60)、yi(t＝60)、zi(t＝60))和排出后经过时间Ti(t)＝40秒，及各粒子P20⁰¹～P20²⁰的排出源强度0.5相对应之后，重新记录到数据记录装置1之中。

此外，使各粒子P40⁰¹～P40²⁰的各空间座标(xi(t＝60)、yi(t＝60)、zi(t＝60))和排出后经过时间Ti(t)＝20秒，及各粒子P40⁰¹～P40²⁰的排出源强度0.9相对应之后，重新记录到数据记录装置1之中。

在以下的计算周期中仍进行同样的处理计算，使各粒子的各空间座标和排出后经过时间Ti(t)及各粒子的排出源强度相对应，重新记录到数据记录装置之中。

下面具体介绍第3步骤(对照图1)的处理。例如在排出后经过时间Ti(t)＝120秒时，要想计算在离图7所示的排出源S有规定距离的规定的方格区域(构成单位体积的单位空间)I、J、K的物质浓度，从数据记录装置i中读出排出后经过时间Ti(t)＝120秒的，存在于该方格区域内的粒子。读出后，若如图7所示存在粒子的情况下，通过累计计算这些粒子各自具有的排出源强度，即可计算出该单位空间内的物质浓度。

也就是说，在图7所示的方格区域内存在

强度为0.3的4个粒子，P00⁰¹、P00⁰⁵、P00¹⁰、P00²⁰；

强度为0.5的3个粒子P20⁰¹、P20⁰⁷、P20¹⁷；

强度为0.9的两个粒子P40⁰⁸、P40¹⁰；以及

强度的0.6的1个粒子P60¹⁷。

因此，通过如下所示累计计算这些粒子的排出源强度，即可计算出该单位空间内的物质浓度为5.1。

(0.3×4)+(0.5×3)+(0.9×2)+(0.6×1)＝5.1

若用普通方式(数学方式)说明上述第1实施方式，则如下所示。在第1实施方式之中，如图8所示，当从排出源S排出物质(气体等)时可根据时间变化，预测处于排出源下风处的方格区域I、J、K中的物质浓度(气体浓度)。而且正如图9(a)所示，当物质的排出量Q固定时，自然能如图9(b)所示预测方格区域I、J、K的浓度经时变化，即使是如图10(a)所示，物质的排出量为随时间变化的排出量Q(t)，也能够如图10(b)所示，预测计算出方格区域的浓度经时变化。

在第1实施之中，不论从排出源S实际排出的物质的排出量是固定值，还是排出量沿时间过程变化，首先把物质的排出量Q(m³/秒)设定为固定值(＝1.0)，采用现有的拉格朗日粒子扩散模型，将物质换算为粒子，使从排出源S每秒种产生N个粒子，将各粒子的轨迹计数后求出表示粒子位置的空间座标(xi(t)、yi(t)、zi(t))。进而除各粒子具有的信息——空间座标(xi(t)、yi(t)、zi(t))之外，将粒子最初产生的时刻起的经过时间——排出后经过时间Ti(t)在每一计算周期内记录到数据记录装置之中。这样一来即可利用现有的拉格朗日粒子扩散模型，计算与测量与经时变化的所有排出量q(t)相对应的浓度分布的经时变化。

为此，须设定与经时变化的物质的排出量q(t)成比例的，与沿排出后经过时间Ti(t)的时间经过的粒子相对应的表示排出源强度的数据。并在某个时刻(t)从数据存储装置中读出表示各粒子位置的空间座标(xi(t)、yi(t)、zi(t))和排出后经过时间Ti(t)。

当排出量固定(Q＝1.0)的情况下，各粒子的排出源强度Q/N(m³)＝1/N，但在排出量q(t)为经时变化的情况下，各粒子的排出源强度则成为q(t-Ti)/N(m³)。

除了各粒子具有的信息——空间座标(xi(t)、yi(t)、zi(t))之外，再次将排出后经过时间Ti(t)以及各粒子的排出源强度qi(t-Ti)/N(m³)重新记录到数据记录装置中。

在利用了拉格朗日粒子模型、排出后经过时间、以及对应于经时变化的排出量q(t)的排出源强度qi(t-Ti)/N(m³)的本实施方式之中，由于空间内的单位体积(＝1m³空气)的各粒子具有的排出源强度不同，因而累计计算各粒子的排出源强度qi(t-Ti)/N(m³)的∑qi(t-Ti)/N(m³)即成为该单位体积内存在的气体量。因此，该空间内的气体浓度即为∑qi(t-Ti)/N(m³)/N(m³)/N(气体m³/空气m³)。

<第2实施方式(排出源为复数)>

当存在多个(J个)排出源，从各个排出源以经时变化不同的排出量(qi(t))排出物质时，在第1实施方式采用的拉格朗日粒子扩散模型的粒子信息(位置、排出后经过时间)的基础之上，通过使各粒子具有排出源识别信息(Si)，即可发挥与第1实施方式相同的功能。

例如图11所示，当有两个排出源S1、S2的情况下，使从排出源S1排出的粒子具有排出源识别信息S1，使从排出源S2排出的粒子具有排出源识别信息S2。并使用拉格朗日粒子扩散模型，将物质换算为粒子，使之从两个排出源S1、S2每秒钟分别排出N个粒子，计数各粒子的轨迹，求出表示粒子位置的空间座标(xi(t)、yi(t)、zi(t))。进而除各粒子具有的信息——空间座标(xi(t)、yi(t)、zi(t))之外，使之具有从粒子最初产生的时刻起的经过时间——排出后经过时间Ti(t)和排出源识别信息S1或S2。

如图12(a)所示，从排出源S1、S2排出的排出量Q，分别设定为固定值(＝1)，采用现有的拉格朗日粒子扩散模型，求出图12(b)所示的方格区域I、J、K内的物质浓度经时变化后记录下来。

接着，若设从排出源S1排出的物质的排出量如图13(a)中的虚线所示，为经时变化的q1(t)，则对于从该排出源S1排出的粒子，用与第1实施方式相同的方法，将粒子的排出源强度设定为q1(t)。并参照该粒子产生时刻的排出源强度q1(t)，针对每一个粒子设定排出源强度。其结果，在所求的方格区域I、J、K内，通过累计计算具有排出源识别信息S 1的粒子的排出源强度，即可求出从排出源S1排出的物质浓度(所求的方格区域I、J、K之内的物质浓度)的经时变化。

同样，若设从排出源S2排出的物质的排出量为q2(t)，则对于从该排出源S2排出的粒子，用与第1实施方式相同的方法，将粒子的排出源强度设定为q2(t)。并针对每一个粒子，参照该粒子产生时刻的排出源强度q2(t)，设定排出源强度。其结果，在所求的方格区域I、J、K之内，通过累计计算具有排出源识别信息S2的粒子的排出源强度，即可求出从排出源S2排出的物质浓度(所求的方格区域I、J、K之内的物质浓度)的经时变化。

如上所述，即可通过将从排出源S1排出的物质浓度的经时变化和从排出S2排出的物质浓度的经时变化相加，求出所求的方格区域I、J、K之内的物质浓度。

<第3实施方式>

第3实施方式是一种发生了有毒气体泄漏事故之后，根据该有毒气体的排出量的实测结果，在短时间内计算当前及将来的浓度分布的方法。

环境评估之类的扩散计算，因不需要马上看到计算结果，所以在确定沿时间经过变化的气体的排出量q(t)之后，花费数日到数月的时间实施气体的扩散状况的预测计算。但在发生有毒气体泄露事故之类的紧急情况下，由于有必要采取紧急避难周围居民的措施，因而必须在尽可能短的时间内得出扩散计算结果。

这种情况下，根据各个时刻的三维空间风速分布，如图14(a)、(b)、(c)所示，用与第1实施方式相同的拉格朗日粒子扩散模型，24小时连续实施固定排放率(Q＝1)的扩散计算。

发生有毒气体泄漏事故后，用设置在排出源的烟囱出口处的气体浓度监测装置等，测定实测气体排放量q(t)。通过按该实测气体排放量q(t)设定各粒的排出源强度q(t)，用与第1实施方式相同的方法，即可修正计算与该有毒气体排放量相对应的浓度经时变化(参照图14(d)(e))。也就是说，针对每个粒子，参照该粒子产生时刻的排出源强度q(t)，设定排出源强度(图14(d)。并在所求的方格区域I、J、K之中，通过累计计算该区域内存在的粒子的排出源强度，即可求出从排出源S排出的物质浓度(所求的方格区域I、J、K之内的物质浓度)的经时变化。

而将来的有毒气体排放量或设定为与当前的实测排出量q(t)等值，或根据另行设定的预测式计算。

作为计算当前及将来的三维空间风速分析的方法，可利用气象局等定期发送的大范围方格(20km)的气象数据(GPV：Grid PointValue)，用大范围气象预测模型(RAMS、MM5等)计算出小范围方格(数公里～数十米)的气象数据(风向、风速、气温)的经时变化。

<第4实施方式>

第4实施方式是一种发生了有毒气体泄漏事故之后，根据泄露地点周边的浓度实测结果在短时间内计算出有毒气体当前排放量的方法。

在发生了有毒气体泄漏事故时，大多数情况下，很难测定从图15(a)所示的排出源——有毒气体泄露源S排出的排出量q(t)。此种情况下，可根据泄漏地点周边的观测点(xk)实测出的浓度(ck(t)观测)的经时变化，采用第1实施方式的拉格朗日粒子扩散模型，短时间内推算出排出量q(t)。

在事故发生前即采用第1实施方式的拉格朗日粒子扩散模型，以固定排出量(Q＝1)(参照图15(b))，预先计算出某个周边观测点(xk)的浓度(Ck(t)计算(参照图15(c))。

该浓度(Ck(t)计算)是根据以观测点(xk)为中心的三维空间体积(Δx×Δy×Δz)内存在的粒子(n个)，用下面的式(10)计算出来的。

Ck(t)计算＝n×∑Q(t-Ti)/N/(Δx×Δy×Δz)...    (10)

此种情况下，排出量Q(t-Ti)为固定值(＝1)。

接着一经测定观测点(xk)的浓度经时变化(Ck(t)观测)，利用式(10)即可求出式(11)

(Ck(t)观测＝n0×∑q0(t-0)/N/(Δx×Δy×Δz)

+n1×∑q1(t-ΔT)/N/(Δx×Δy×Δz)

-n2×∑q2(t-2·ΔT)/N/(Δx×Δy×Δz)

+n3×∑q3(t-3·ΔT)/N/(Δx×Δy×Δz)

+n4×∑q4(t-4·ΔT)/N/(Δx×Δy×Δz)

+......

+n1×∑q1(t-t)/N/(Δx×Δy×Δz)

+C0(t)                      ......(11)

式中的n0、n1、n2、n3、以及n1分别为排放开始后，0秒前、Δt秒前、(Δt×2)秒前、(Δt×3)秒前以及t秒前排放出的粒子，是存在于以观测点(xk)为中心的三维空间体积(Δx ×Δy ×Δz)内的粒子总数。

此外，式中的C0(t)被称之为背景浓度，是从构成计算对象的排出源以外的地方排出，与观测地点无关的存在的浓度，为固定值或随时间而变化。

式中的q0、q1、q2、q3以及q1分别为从t秒时起0秒前、ΔT秒前、(ΔT×2)秒前、(ΔT×3)秒前以及t(＝Δt×1)秒前的排放率。

由于Ck t)观测是24小时连续测定的，因而可测定从泄漏开始时刻到0秒后、ΔT秒后、(ΔT ×2)秒后、(ΔT ×3)秒后、以及t(ΔT×1)秒后(1+1)个以上。

此外，当有k个观测点(xk)的情况下，可取得k×(1+1)个观测数据。

由于式(11)中的未知数为q0、q1、q2、q3以及q1的(1+1)个，已知数Ck(t)观测有k×(1+1)个以上，因而未知数少于已知数。这种情况下采用最小二乘法将与观测浓度Ck t)观测对应的最小平方误差设为最小，即可决定未知数q0、q1、q2、q3及q1。图16即是表示第4实施方式的计算状态的流程图。

<第5实施方式>

第5实施方式是一种在发生有毒气体泄漏之类的事故之后，根据泄漏地点周边的浓度实测结果，在短时间内推算有毒气体的当前排出量、计算其浓度分布的方法。

根据采用第4实施方式，从在泄漏地点周边的观测点(xk)实测的浓度(Ck(t)观测)的经时变化，采用第1实施方式的拉格朗日粒子扩散模型，即可在短时间内推算出排出量q(t)。

正如图17所示，在泄漏事故发生之前，即采用第1实施方式的拉格朗日粒子扩散模型，以固定排放率(Q＝1)预先计算出某周边观测点(xk)的浓度(Ck(t)计算)。

该浓度(Ck(t)计算)是根据存在于以观测点(xk)为中心的三维空间体积(Δx×Δy×Δz)内的(n个)粒子，用式(12)计算出来的。

Ck(t)计算＝n×∑Q(t-Ti)/N/(Δx ×Δy×Δz)...(12)

这种情况下排放率Q(t-Ti)为固定值(＝1)。

接着一经测定观测点(xk)的浓度经时变化(Ck(t)观测)，即可利用式(12)求出式(13)。

(Ck(t)观测＝n0×∑q0(t-0)/N/(Δx×Δy×Δz)

+n1×∑q1(t-ΔT)/N/(Δx×Δy×Δz)

+n2×∑q2(t-2·ΔT)/N/(Δx×Δy×Δz)

+n3×∑q3(t-3·ΔT)/N/(Δx×Δy×Δz)

+n4×∑q4(t-4·ΔT)/N/(Δx×Δy×Δz)

+......                         (13)

+n1×∑q1(t-t)/N/(Δx ×Δy×Δz)

+C0(t)

式中的n0、n1、n2、n3以及n1分别为排放开始后0秒前、ΔT秒前、(ΔT ×2)秒前、(ΔT ×3)秒前以及t秒前排放出的粒子，是存在于以观测点为中心的三维空间(Δx×Δy×Δz)内的粒子总数。

此外，式中的C0(t)被称之为背景浓度，是从构成计算对象的排出源以外的地方排出，与观测地点无关的存在的浓度，为固定值或随时间而变化。

式中的q0、q1、q2、q3以及q1分别为从t秒时起，0秒前，ΔT秒前、(ΔT×2)秒前、(ΔT×3)秒前以及t(＝ΔT×1)秒前的排放率。

由于Ck(t)观测是24小时连续测定的，因而可测定出从泄漏开始时刻到0秒后、ΔT秒后、(ΔT×2)秒后、(ΔT×3)秒后以及t(＝ΔT×1)秒后(1+1)个以上。

此外，当有k个观测点(xk)的情况下，可取得k×(1+1)个观测数据。

由于式(13)中的未知数为q0、q1、q2、q3以及q1的(1+1)个，而已知数Ck(t)观测为k×(1+1)个以上，所以未知数小于已知数。这种情况下，采用最小二乘法将与观测浓度Ck(t)观测对应的最小平方误差设为最小，即可决定未知数的q0、q1、q2、q3及q1。

若采用该推算值q0、q1、q2、q3及q1，以第1实施方式的方法进行浓度经时变化的修正计算，即可计算出各自的每一经过时间的浓度分布。

<第6实施方式>

第6实施方式是在发生有毒气体泄漏事故之后，根据排放量实测结果，在短时间内用因特网提供当前及将来的浓度分布计算结果的系统。

由于在发生有毒气体泄漏之类事故的紧急情况下，必须紧急采取避难周边居民的措施，因而要求在尽可能短的时间内得出扩散计算结果。但由于此类泄漏事故不知会在何时发生，所以在政府监察部门——消防、警察以及自治团体和企业的各工厂必须以24小时体制管理运作该计算机系统。

由于此种管理运作需要花费很大的人力与财力以及高超的计算机操作技术，因此除常设大规模危机管理系统的机构之外，很难管理运作。

作为该问题的补救办法，提出了利用最近的因特网的信息提供系统的方案。

在该第6实施之中，正如图18所示，在与政府监察部门的消防、警察以及自治团体等职能部门10及企业的各工厂11不同的其它地方设置安全分析中心12。安全分析中心12通过因特网等信息传输手段接收从气象厅13等的气象数据发送设施发送的气象观测数据，平常时通过使用大型机算机进行计算，根据各个时刻的三维空间风速分布，用与第1实施方式相同的拉格朗日粒子扩散模型，24小时连续实施固定排放率(Q＝1)的扩散计算。

发生有毒气体泄漏事故后，若能从设置在企业11的烟囱出口处的气体浓度测定装置等处了解到实测的有毒气体排出量q(t)，则把该有毒气体排放量q(t)通过因特网等信息传输手段传输给安全分析中心12。这样即可在安全分析中心12采用与第1实施方式相同的方法，修正计算与该有毒气体排放率相对应的浓度经时变化。

安全分析中心12将计算出的浓度计算结果发送给消防、警察及自治团体等政府监察部门10。职能部门10根据其浓度，对企业11及工厂周边居民14发出避难通告。

而将来的有毒气体排放量或设定为与当前的实测排出量q(t)等值，或根据另行设定的预测试计算。在有多家企业存在的情况下，在安全分析中心执行第二实施方式的计算，来预测浓度经时变化。

作为计算当前及将来的三维空间风速分布的方法，可利用气象厅13等处定期发送的大范围方格(20km)的气象观测数据(GPV：GridPointValue)，用大范围气象预测模型(RAMS、MM5等)计算出小范围方格(数公里～数十米)的气象数据(风向、风速、气温)的经时变化。

(发明效果)

正如以上的实施方式及具体说明中所示，本发明涉及的扩散物质的扩散状态预测方法，为了预测从排出源排往大气中的物质在大气中的扩散状况，将上述物质换算为多数粒子，设定出每一计算周期内从排出源位置产生预先设定个数的粒子；同时通过将包括排出源位置的区域内的多个地点的，表示沿时间过程变化的风向、风速的风速场数据代入计算粒子扩散状态的扩散方程式，求出各个粒子的扩散速度，并在每个计算周期内根据该扩散速度求出表示各个粒子所在的空间位置的空间座标，同时计算测量出从上述粒子最初产生的时刻起的经过时间——排出后经过时间，使各个计算周期内的各个粒子的空间座标与各个粒子的排出后经过时间相对应，然后记录到数据记录装置之中；按照排出物质的排出后经过时间的伴随时间过程的排出量的变化比例，预先设定与排出后经过时间的沿时间过程的粒子相对应的排出源强度数据；读出记录在上述数据记录装置中的，每一计算周期的各个粒子的空间座标和各个粒子的排出后经过时间，同时参照读出的排出后经过时间求出各个粒子产生的时刻，根据上述排出源强度数据，求出该时刻的各个粒子的排出源强度，使排出源强度每一计算周期的各个粒子的空间座标和各粒子的排出后经过时间相对应之后由上述数据记录装置重新记录；规定计算周期内的规定区域的上述物质浓度，则设定为通过累计计算该规定计算周期内的存在于该规定区域的所有粒子的排出源强度求出。

因此，即使从排出源排出的量不同，也能在短时间内计算出特定区域内的物质浓度。

此外，本发明的扩散物质的扩散状况预测方法，为了预测从多个排出源排往大气中的物质在大气中的扩散状况，将上述物质换算为多数粒子，设定出每一计算周期从排出源位置分别产生预先设定个数的粒子；同时通过将包括排出源的位置的区域内的多数地点的，表示沿时间过程变化的风向、风速的风速场数据，代入计算粒子扩散状态的扩散方程式，求出各粒子的扩散速度，根据该扩散速度，求出每一计算周期内表示各个粒子所在的空间位置的空间座标，同时计算测量出上述粒子最初产生的时刻起的经过时间——排出后经过时间，使各个计算周期内的各个粒子的空间座标和各个粒子的排出后经过时间及识别排出源的排出源识别信息相对应之后，记录到数据记录装置之中；按照从排出源排出的物质的排出后经过时间的时间过程伴随的排出量的变化比例，分别预先设定每个排出源与沿排出后经过时间的时间过程的粒子相对应的排出源强度数据；读出记录在上述记录装置中的，各个计算周期的各个粒子的空间座标和各粒子的排出后经过时间以及各个粒子的排出源识别信息，同时参照读出的排出源识别信息，根据与该粒子产生的排出源所对应的上述排出源强度数据，求出粒子产生时刻的各个粒子的排出源强度，使每一计算周期的各个粒子的空间座标和各个粒子的排出后经过时间及排出源强度相对应，由上述数据记录装置重新记录；规定的计算周期内的规定区域的上述物质浓度设定为通过累计计算该规定的计算周期内的存在于该规定区域的所有粒子的排出源强度求出。

因此，即使从排出源排出的量不同，也能在短时间内计算出特定区域内的物质浓度。此外，即使从多个排出源排出物质，也能准确计算出物质的浓度。

此外，在本发明的扩散物质的扩散状况预测方法之中，既可以通过实测求出从上述排出源实际排出的物质的浓度后设定上述排出源强度数据，也可根据在上述排出源周围的观测点实测的物质浓度的经时变化设定上述排出源强度数据。

因此，即使预先未能求出排出源排出的物质的浓度数据，也能用实测数据进行浓度计算。

此外，本发明的扩散物质的扩散状况预测系统包括：扩散物质排往大气中时，实测扩散物质的浓度，发送表示扩散物质的排出量数据的企业；发送气象观测数据的气象数据发送设施；向上述企业及周边居民发布避难通告的政府监察部门；进行扩散物质的扩散状况预测计算处理，计算规定区域内的物质浓度的安全分析中心。设定为上述安全分析中心通过信息传输手段接收上述企业发送的表示扩散物质排出量的数据，接收上述气象数据发送设施发送的气象观测数据。上述政府监察部门通过信息传输手段接收上述安全分析中心发送的物质浓度。上述政府监察部门根据接收到的浓度报告发布避难通告。

因此，政府监察部门可根据安全分析中心计算出的信息迅速发出避难通告，紧急采取避难周边居民的措施。

Claims (5)

1.一种扩散物质的扩散状况预测方法，其特征在于：

为了预测从排出源排往大气中的物质在大气中的扩散状况，将上述物质换算为多数粒子，当每个计算周期内从排出源位置产生预定个数的粒子时进行设定；

同时通过将包括排出源位置在内的区域内的多数地点的表示随时间变化的风向、风速的风速场数据代入计算粒子扩散状况的扩散方程式，求出各个粒子的扩散速度，根据该扩散速度，求出在每一计算周期内表示各粒子所在的空间位置的空间座标，同时计测从最初产生上述粒子的时刻起的经过时间即排出后经过时间，使各计算周期内的各粒子的空间座标和各粒子的排出后经过时间相对应，并记录到数据记录装置之中；

与排出物质的排出量随排出后经过时间的变化成比例，预先设定随排出后经过时间变化的与粒子相对应的排出源强度数据；

读出记录在上述数据记录装置中的，每一计算周期的各粒子的空间座标和各粒子的排出后经过时间，同时参照读出的排出后经过时间，求出各粒子的发生时刻，根据上述排出源强度数据，求出该时刻的各粒子的排出源强度，使排出源强度与每一计算周期的各粒子的空间座标和各粒子的排出后经过时间相对应，然后由上述数据记录装置重新记录；

规定的计算周期内的规定区域的上述物质的浓度，通过累计计算该规定计算周期内的存在于该规定区域的所有粒子的排出源强度求出。

2.一种扩散物质的扩散状况预测方法，其特征在于：为了预测从多个排出源排往大气中的物质在大气中的扩散状况，将上述物质换算为多数粒子，当在每个计算周期内从各个排出源位置各自产生预定个数的粒子时进行设定；

同时通过将包括排出源位置在内的区域内的多数地点的表示随时间变化的风向、风速的风速场数据代入计算粒子扩散状态的扩散方程式，求出各粒子的扩散速度，根据该扩散速度，求出每一计算周期内表示该粒子所在的空间位置的空间座标，同时计测上述粒子最初产生的时刻起的经过时间即排出后经过时间，使各计算周期内的各粒子的空间座标和各粒子的排出后经过时间及识别排出源的排出源识别信息相对应，并记录到数据记录装置之中；

与从各个排出源排出的物质的排出量随排出后经过时间的变化成比例，预先分别设定随排出后经过时间变化的粒子的排出源强度数据；

读出记录在上述数据记录装置中的，每一计算周期的各粒子的空间座标和各粒子的排出后经过时间以及各粒子的排出源识别信息，同时参照读出的排出后经过时间，求出各粒子产生时刻，参照读出的排出源识别信息，根据与产生该粒子的排出源相对应的上述排出源强度数据，求出产生该粒子的时刻的各粒子的排出源强度，使排出源强度与每一计算周期的各粒子的空间座标和各粒子的排出后经过时间相对应，由上述数据记录装置重新记录；

规定的计算周期内的规定区域的上述物质的浓度，通过累计计算该规定的计算周期内的存在于该规定区域的所有粒子的排出源强度求出。

3.根据权利要求1或2所述的扩散物质的扩散状况预测方法，其特征在于：通过实测求出从上述排出源实际排出的物质的浓度后设定上述排出源强度数据。

4.根据权利要求1或2所述的扩散物质的扩散状况预测方法，其特征在于：上述排出源强度数据依据在上述排出源周围的观测点实测的物质浓度的经时变化进行设定。

5.一种扩散物质的扩散状况预测系统，其特征在于，该系统包括：

当扩散物质排往大气中时，实测扩散物质的浓度，并发送表示扩散物质排出量的数据的企业；

发送气象观测数据的气象数据发送设施；

向上述企业及上述企业周围的居民发布避难通告的政府监察部门；以及

安全分析中心，其根据权利要求1或2的预测方法进行计算处理以计算规定区域内的物质浓度，上述企业发送的表示扩散物质排出量的数据以及上述气象数据发送设施发送的气象观测数据通过信息传输手段传送到上述安全分析中心，上述安全分析中心的物质浓度报告通过信息传输手段传送到上述政府监察部门，

上述政府监察部门根据传送来的物质浓度报告发布避难通告。

Images