1 在线模拟示例

在线示例https://simu.cool/vue/

2 高斯模型

2.1 高斯烟团模型

   突发性泄漏事故中，经常发生污染源在短时间内突然释放大量的有害气体，此时对地面污染浓度的求解适合采用高斯烟团模型。烟团模型假定污染气云的体积沿水平和垂直方向增长，模拟污染气云在时间和空间上的变化。

2.2 高斯烟羽模型

  高斯烟羽模式是计算释入大气中的气载污染物下风向浓度的应用最广的方法。此模式假定烟羽中污染物浓度分布在水平方向和垂直方向都遵循高斯分布。对于在恒定气象条件(指风向、风速、大气稳定度不随时间而变)高架点源的连续排放，在考虑了烟羽在地面的全反射后，下风向任一点的污染物浓度C (r，y，z)可由高斯烟羽公式进行模拟。

高斯扩散公式的建立有如下假设：
①风的平均流场稳定，风速均匀，风向平直；
②y、z轴方向符合正态分布；
③污染物在输送扩散中质量守恒；
④污染源的源强均匀、连续

    如图所示，有效源位于坐标原点o处，平均风向与x轴平行，并与x轴正向同向。假设点源在没有任何障碍物的自由空间扩散，不考虑下垫面的存在。大气中的扩散是具有y与z两个坐标方向的二维正态分布，当两坐标方向的随机变量独立时，分布密度为每个坐标方向的一维正态分布密度函数的乘积。
    经过推理得出计算公式：

C—空间点（x，y，z）的污染物的浓度，mg/m3；
σy、σz—分别为水平、垂直方向的标准差，即y、x方向的扩散参数
u —为平均风速
x —为风向轴上空间点到源的距离
y —为风向轴垂直方向上空间点到源的距离
z —为空间点的高度

σy、σz与大气稳定度和水平距离x有关，并随x的增大而增加。通过理论或经验的方法可得σ＝f（x），可求出最大浓度点离源的距离x，具体可查阅我国GB384091《制定地方大气污染物排放标准的技术方法》。

2.2.1 方程

式中c为污染物浓度(单位：mg/m3)
Q为源强(单位：mg/s)
u为泄漏高度的平均风速(单位：m/s)
y、z分别用浓度标准偏差表示的y轴及z轴上的扩散参数
H为泄漏有效高度(单位：m) 

2.2.2 适用条件

(1)污染物浓度在y、z轴上的分布符合高斯分布（正态分布）；

(2)在全部空间中风速是均匀的、稳定的；

(3)源强是连续均匀的；

(4)在扩散过程中污染物质量是守恒的（不考虑转化）。

3 高斯羽烟模型参数

3.1 坐标系

高斯模式的坐标系如下图所示，其原点为排放点（无界点源或地面源）或高架源排放点在地面的投影点，x轴正向为平均风向，y轴在水平面上垂直于x轴，正向在x轴的左侧，z轴垂直于水平面xoy，向上为正向，即为右手坐标系。

3.2 烟气抬升

有效源高H,等于烟囱的几何高度与烟气抬升高度△H之和，即

对于一确定的烟囱，其几何高度H是确定值，只要计算出抬升高度△H值，就可以求出有效源高H。我国“制订地方大气污染物排放标准的技术方法”(GB/T13201-91)中规定选用Briggs 公式作为适用计算模型，请参考相关文献。

3.3 大气稳定度

大气稳定度指大气中某一高度上的气团在垂直方向上的相对稳定程度。如果给一团空气一个初始作用力，使其作向上的垂直运动，垂直运动的气块在外力消失后，又逐渐回到原来的位置，这种状况的大气是稳定的;当外力消失后，气块仍继续上升，甚至加速前进，这种状况的大气是不稳定的;当外力消失后，气块停留在其已到达的位置，既不上升也不下降，这种状况的大气处于中性状态。

我国《环境影响评价技术导则》中推荐了用常规地面观测资料划分大气稳定度的方法。大气稳定度的分类方法采用经过修正的帕斯奎尔(Pasquill)稳定度分级法(Ps)，将大气扩散稳定度分为强不稳定、不稳定、弱不稳定、中性、弱稳定和稳定六级，分别用A、B、C、D、E、F来表示。确定等级时首先根据云量与太阳高度角按下表查出太阳辐射等级数，再由太阳辐射等级数与地面风速按下表查找稳定度等级。

 

3.4 有风条件下的扩散系数

    1. 有风时扩散参数、:
    扩散参数是大气稳定度和距离的函数，其表达式可写为以下幂指数的

 

     如无实测值时,上述各指数、系数的值可按下述方法选取:

     (1) 平原地区农村及城市远郊区的扩散参数如下:A、B、C级稳定度直接由表2-9和表2-10查算，D、E、F级稳定度则需向不稳定方向提半级后由表2-9和表2-10查算。
     (2) 工业区或城区中的点源，其扩散参数选取方法如下:
A、B级不提级，C级提到B级，D、E、F级向不稳定方向提一级，再按表2-9和表2-10查算。

4 参考代码

 /**
     * 获取横向扩散参数回归系数
     *
     * @param windSpeed 风速
     * @param Q         大气稳定度
     * @return
     */
    public double getGamma01Light(double windSpeed, final String Q) {
        double gamma = 0;
        if (Q.equalsIgnoreCase("A")) {
            if (windSpeed < 0.5) {
                gamma = 0.93;
            } else if (windSpeed < 1.5) {
                gamma = 0.76;
            }
        } else if (Q.equalsIgnoreCase("B")) {
            if (windSpeed < 0.5) {
                gamma = 0.76;
            } else if (windSpeed < 1.5) {
                gamma = 0.56;
            }
        } else if (Q.equals("C")) {
            if (windSpeed < 0.5) {
                gamma = 0.55;
            } else if (windSpeed < 1.5) {
                gamma = 0.35;
            }
        } else if (Q.equals("D")) {
            if (windSpeed < 0.5) {
                gamma = 0.47;
            } else if (windSpeed < 1.5) {
                gamma = 0.27;
            }
        } else if (Q.equals("E")) {
            if (windSpeed < 0.5) {
                gamma = 0.44;
            } else if (windSpeed < 1.5) {
                gamma = 0.24;
            }
        } else if (Q.equals("F")) {
            if (windSpeed < 0.5) {
                gamma = 0.44;
            } else if (windSpeed < 1.5) {
                gamma = 0.24;
            }
        }

        return gamma;
    }

    /**
     * 获取铅直扩散参数回归系数
     *
     * @param windSpeed 风速
     * @param Q         大气稳定度
     * @return
     */
    public double getGamma02Light(double windSpeed, final String Q) {
        double gamma = 0;
        if (Q.equalsIgnoreCase("A")) {
            if (windSpeed < 0.5) {
                gamma = 0.15;
            } else if (windSpeed < 1.5) {
                gamma = 1.57;
            }
        } else if (Q.equalsIgnoreCase("B")) {
            if (windSpeed < 0.5) {
                gamma = 0.47;
            } else if (windSpeed < 1.5) {
                gamma = 0.47;
            }
        } else if (Q.equals("C")) {
            if (windSpeed < 0.5) {
                gamma = 0.21;
            } else if (windSpeed < 1.5) {
                gamma = 0.21;
            }
        } else if (Q.equals("D")) {
            if (windSpeed < 0.5) {
                gamma = 0.12;
            } else if (windSpeed < 1.5) {
                gamma = 0.12;
            }
        } else if (Q.equals("E")) {
            if (windSpeed < 0.5) {
                gamma = 0.07;
            } else if (windSpeed < 1.5) {
                gamma = 0.07;
            }
        } else if (Q.equals("F")) {
            if (windSpeed < 0.5) {
                gamma = 0.05;
            } else if (windSpeed < 1.5) {
                gamma = 0.05;
            }
        }

        return gamma;
    }

package com.planet.engine.air.param;

import lombok.Data;

/**
 * 运算参数
 */
@Data
public class AirProcParam {
    /**
     * 经纬度坐标系：污染点原点x
     */
    private double lon = 108.14D;

    /**
     * 经纬度坐标系：污染点原点y
     */
    private double lat = 36.27D;

    /**
     * 起算点即原点高程
     */
    private Integer startEle = 0;

    /**
     * 投影横坐标
     */
    private double x = 0;

    /**
     * 投影纵坐标
     */
    private double y = 0;

    /**
     * 当前计算点z坐标
     * 不做地形改正时，为计算点离地面高
     * 做地形改正时，为计算点的高程改正值
     */
    private double z = 20;

    /**
     * 模拟时间
     * 单位秒
     */
    private int timeSec = 600;

    /**
     * 扇形开度，算法要求不超过45度
     */
    private double sectorAngle = 90;

    /**
     * 扇形弧段加密点数
     */
    private int pointsNumSector = 6;

    /**
     * 扇形弧段加密点数
     */
    private int pointsNumCircle = 24;

    /**
     * 高斯烟团或者高斯烟羽
     * 分别用0和1表示
     */
    private Integer type = 0;

    /**
     * 网格点间距，单位：米
     */
    private float sideLen = 200f;

    /**
     * 密度阈值，单位mg/m3
     * 超过此值的需要做叠加分析
     */
    private Integer densityThr = 3;
}

package com.planet.controller;

import com.planet.common.annotations.LimitRequest;
import com.planet.engine.air.param.GaussAirParam;
import com.planet.model.SecurityUser;
import com.planet.model.enums.SimuType;
import com.planet.model.enums.StatusType;
import com.planet.service.GaussAirService;
import com.planet.service.SimuHistService;
import io.swagger.annotations.Api;
import io.swagger.annotations.ApiOperation;
import org.springframework.beans.factory.annotation.Autowired;
import org.springframework.web.bind.annotation.*;


/**
 * Created by 之乎者也
 * 
 */
@RestController
@RequestMapping("/api/air")
@Api(value = "AirSimuApi", description = "空气污染扩散模拟模型")
public class AirSimuApi extends MyBaseApi {

    @Autowired
    private GaussAirService gaussAirService;

    @Autowired
    private SimuHistService simuHistService;

    @LimitRequest
    @ResponseBody
    @RequestMapping(value = "/gauss", method = RequestMethod.POST)
    @ApiOperation(value = "高斯烟羽扩散模拟接口")
    public Result gaussSimu(@RequestBody GaussAirParam params) {
        // 验证参数有效性，例如风速
        String result = this.gaussAirService.calc(params);

        // 异步保存历史纪录
        SecurityUser securityUser = this.getLoginUser();
        String paramsStr = this.getReqParams(params);
        if(params.getProcParams().getType().equals(1)){
            this.simuHistService.save(securityUser, SimuType.GAUSS_PLUME, paramsStr, StatusType.SUCCESS.code(), result);
        }else{
            this.simuHistService.save(securityUser, SimuType.GAUSS_BALL, paramsStr, StatusType.SUCCESS.code(), result);
        }


        return this.getResult(ResultCode.SUCCESS, result);
    }
}