高斯扩散模型-高斯烟羽大气污染扩散模型

高斯，大名鼎鼎

约翰·卡尔·弗里德里希·高斯（德语：Johann Carl Friedrich Gauß;  ，英语：Gauss，拉丁语：Carolus Fridericus Gauss，1777年4月30日－1855年2月23日），生于布伦瑞克，卒于哥廷根。德国著名数学家、物理学家、天文学家、几何学家，大地测量学家。享有“数学王子”的美誉。

这个人太厉害了,今天用到的核心便是高斯正态分布

正态分布概念是由德国的数学家和天文学家Moivre于1733年首次提出的，但由于德国数学家Gauss率先将其应用于天文学研究，故正态分布又叫高斯分布，高斯这项工作对后世的影响极大，他使正态分布同时有了“高斯分布”的名称，后世之所以多将最小二乘法的发明权归之于他，也是出于这一工作。

连续点源扩散

    连续点源一般指排放大量污染物的烟囱、放散管、通风口等。排放口安置在地面的称为地面点源，处于高空位置的称为高架点源。

1. 大空间点源扩散

高斯扩散公式的建立有如下假设：
①风的平均流场稳定，风速均匀，风向平直；
②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. 高架点源扩散

    在点源的实际扩散中，污染物可能受到地面障碍物的阻挡，因此应当考虑地面对扩散的影响。处理的方法是，或者假定污染物在扩散过程中的质量不变，到达地面时不发生沉降或化学反应而全部反射；或者污染物在没有反射而被全部吸收，实际情况应在这两者之间。

（1） 高架点源扩散模式

点源在地面上的投影点o 作为坐标原点，有效源位于z 轴上某点， z ＝ H。高架有效源的高度由两部分组成，即H＝h＋Δ h，其中h为排放口的有效高度， Δh 是热烟流的浮升力和烟气以一定速度竖直离开排放口的冲力使烟流抬升的一个附加高度

    经过推理得出计算公式：

C—空间点（x，y，z）的污染物的浓度，mg/m3；
σy、σz—分别为水平、垂直方向的标准差，即y、x方向的扩散参数

u —为平均风速
x —为风向轴上空间点到源的距离
y —为风向轴垂直方向上空间点到源的距离
z —为空间点的高度
H—为烟囱高度

近期用高架点源扩散模式做了个效果，大概实现还有很大的优化空间

（2） 地面全部反射时的地面浓度

实际中，高架点源扩散问题中最关心的是地面浓度的分布状况，尤其是地面最大浓度值和它离源头的距离。

y 方向的浓度以x 轴为对称轴按正态分布；沿x 轴线上，在污染物排放源附近地面浓度接近于零，然后顺风向不断增大，在离源一定距离时的某处，地面轴线上的浓度达到最大值，以后又逐渐减小。
高架点源的地面浓度公式中令y＝0则可得到沿x轴线上的浓度分布：

3. 地面点源扩散

对于地面点源，则有效源高度H＝0。当污染物到达地面后被全部反射时，可令H＝0，即得出地面连续点源的高斯扩散公式:

连续线源扩散

    当污染物沿一水平方向连续排放时，可将其视为一线源，如汽车行驶在平坦开阔的公路上。线源在横风向排放的污染物浓度相等，这样，可将点源扩散的高斯模式对变量y积分，即可获得线源的高斯扩散模式。但由于线源排放路径相对固定，具有方向性，若取平均风向为x轴，则线源与平均风向未必同向。所以线源的情况较复杂，应当考虑线源与风向夹角以及线源的长度等问题。

连续面源扩散

    当众多的污染源在一地区内排放时，如城市中家庭炉灶的排放，可将它们作为面源来处理。因为这些污染源排放量很小但数量很大，若依点源来处理，将是非常繁杂的计算工作。

大气的稳定度

大气稳定度是指大气中的某一气团在垂直方向上的稳定程度。

一团空气受到某种外力作用而产生上升或者下降运动，当运动到某—位置时消除外力，此后气团的运动可能出现三种情况：
①气团仍然继续加速向前运动，这时的大气称为不稳定大气；
②气团不加速也不减速而作匀速运动，或趋向停留在外力去除时所处的位置，这时的大气称为中性大气；
③气团逐渐减速并有返回原先高度的趋势，这时的大气称为稳定大气。

分类

大气稳定度与天气现象、时空尺度及地理条件密切相关，其级别的准确划分非常困难。
目前国内外对大气稳定度的分类方法已多达10余种，应用较广泛的有帕斯奎尔（Pasquill）法和特纳尔(Turner)法。
帕斯奎尔法用地面风速(距离地面高度10m)、白天的太阳辐射状况（分为强、中、弱、阴天等）或夜间云量的大小将稳定度分为A～F六个级别

帕斯奎尔法虽然可以利用常规气象资料确定大气稳定度等级，简单易行，应用方便，但这种方法没有确切地描述太阳的辐射强度，云量的确定也不准确，较为粗略，为此特纳尔作了改进与补充。
特纳尔方法首先根据某地、某时及太阳倾角的太阳高度θ h 和云量( 全天空为10 分制) ，确定太阳辐射等级，再由太阳的辐射等级和距地面高度10m的平均风速确定大气稳定度的级别。
我国采用特纳尔方法

模型计算:

  //获取稳定度对应Qy，x为距离，Q大气稳定度
  function GetQy(x, Q) {
    var a = 0.0;
    var r = 0.0;
    if (Q == "A") {
      if (x <= 1000) {
        a = 0.901074;
        r = 0.425809;
      }
      else {
        a = 0.850934;
        r = 0.602052;
      }
    }
    else if (Q == "B") {
      if (x <= 1000) {
        a = 0.914370;
        r = 0.281846;
      }
      else {
        a = 0.865014;
        r = 0.396353;
      }
    }
    else if (Q == "B-C") {
      if (x <= 1000) {
        a = 0.919325;
        r = 0.229500;
      }
      else {
        a = 0.875086;
        r = 0.314238;
      }
    }

    else if (Q == "C") {
      if (x <= 1000) {
        a = 0.924279;
        r = 0.177154;
      }
      else {
        a = 0.885157;
        r = 0.232123;
      }
    }
    else if (Q == "C-D") {
      if (x <= 1000) {
        a = 0.926849;
        r = 0.143940;
      }
      else {
        a = 0.886940;
        r = 0.189396;
      }
    }

    else if (Q == "D") {
      if (x <= 1000) {
        a = 0.929418;
        r = 0.110726;
      }
      else {
        a = 0.888723;
        r = 0.146669;
      }
    }

    else if (Q == "D-E") {
      if (x <= 1000) {
        a = 0.925118;
        r = 0.0985631;
      }
      else {
        a = 0.892794;
        r = 0.124308;
      }
    }

    else if (Q == "E") {
      if (x <= 1000) {
        a = 0.920818;
        r = 0.0864001;
      }
      else {
        a = 0.896864;
        r = 0.101947;
      }
    }

    else {
      if (x <= 1000) {
        a = 0.929418;
        r = 0.0553634;
      }
      else {
        a = 0.888723;
        r = 0.0733348;
      }
    }

    var Qz = Math.pow(x, a) * r;
    return Qz;
  }

  //获取稳定度对应Qz，x为距离，Q大气稳定度
  function GetQz(x, Q) {
    var a = 0.0;
    var r = 0.0;
    if (Q == "A") {
      if (x <= 300) {
        a = 1.12154;
        r = 0.0799904;
      }
      else if (x > 300 || x <= 500) {
        a = 1.51360;
        r = 0.00854771;
      }
      else {
        a = 2.10881;
        r = 0.000211545;
      }

    }

    else if (Q == "B") {
      if (x <= 500) {
        a = 0.964435;
        r = 0.127190;
      }
      else {
        a = 0.109356;
        r = 0.057025;
      }
    }
    else if (Q == "B-C") {
      if (x <= 500) {
        a = 0.941015;
        r = 0.114682;
      }
      else {
        a = 1.00770;
        r = 0.0757182;
      }
    }
    else if (Q == "C") {
      a = 0.917595;
      r = 0.106803;
    }

    else if (Q == "C-D") {
      if (x <= 2000) {
        a = 0.838628;
        r = 0.126152;
      }
      else if (x > 2000 || x <= 10000) {
        a = 0.756410;
        r = 0.235667;
      }
      else {
        a = 0.815575;
        r = 0.136659;
      }
    }
    else if (Q == "D") {
      if (x <= 1000) {
        a = 0.826212;
        r = 0.104634;
      }
      else if (x > 1000 || x <= 10000) {
        a = 0.632023;
        r = 0.400167;
      }
      else {
        a = 0.55536;
        r = 0.810763;
      }
    }
    else if (Q == "D-E") {
      if (x <= 2000) {
        a = 0.776864;
        r = 0.111771;
      }
      else if (x > 2000 || x <= 10000) {
        a = 0.572347;
        r = 0.5289922;
      }
      else {
        a = 0.499149;
        r = 1.03810;
      }
    }
    else if (Q == "E") {
      if (x <= 1000) {
        a = 0.788370;
        r = 0.0927529;
      }
      else if (x > 1000 || x <= 10000) {
        a = 0.565188;
        r = 0.433384;
      }
      else {
        a = 0.414743;
        r = 1.73241;
      }
    }
    else {
      if (x <= 1000) {
        a = 0.784400;
        r = 0.0620765;
      }
      else if (x > 1000 || x <= 10000) {
        a = 0.525969;
        r = 0.370015;
      }
      else {
        a = 0.322659;
        r = 2.40691;
      }
    }


    var Qy = Math.pow(x, a) * r;
    return Qy;

  }