垃圾焚烧厂经济赔偿问题——基于高斯模型.doc

垃圾焚烧厂的经济补偿问题基于高斯模型摘要：目前我国三分之二以上的城市面临“垃圾围城”问题，垃圾焚烧这种具有高效益的方法正逐步成为中国垃圾处理的主要手段之一。然而，近年来由于垃圾焚烧电厂产生环境污染，遭到了居民的极大阻挠，因此研究成为民众所信服的全方位垃圾焚烧厂环境监控体系并制定合理的经济补偿方案对于推行“垃圾焚烧”具有关键性意义。本文便针对这个问题建立科学模型并提出合理化的建议。问题一：运用高斯烟羽模型，得到各因素对污染物扩散分布的影响及各类污染物浓度的初步函数。以垃圾场为圆心，5千米为半径画出废气污染范围，并将此圆按风向分为8个扇形区域，再以1千米为间距，在每个风向上取5个监测点，共40个。继而运用层次分析法计算颗粒物、的权重，分别为0.1765，0.4706，0.3529。则每个监测点的总污染浓度为：。本文是按照区域来对居民进行赔偿，因此需要根据监测点来求出每个区域的平均污染浓度，定义为每个区域的4个顶点的总污染浓度的平均值。假设政府共有赔偿金W元，则每个区域可按该区域的平均污染浓度得到相应的赔偿金。分发到每个区域的赔偿金按人口平均分配。问题二：由于焚烧炉的除尘装置出现故障或损坏，各污染物的指标将严重超标，因此，各个监测区的污染物浓度也会大大增加。利用SPSS画出附件2中的各污染物三月里焚烧炉排放浓度散点图。观察后发现，只有颗粒物的浓度在某些天内超标，而浓度和浓度均符合在国家标准。因此本文由颗粒物的故障条件来类比故障时的浓度和浓度。定义颗粒物浓度为25以上为故障浓度，因此对附件2中所提供的超过25的数据求平均，得到30，超过国家标准50%。则故障时的各监测点的浓度为：。最后，按照问题一中给出的补偿金额计算方式，得出新的补偿方案。关键词：高斯烟羽模型 层次分析法 垃圾焚烧污染物 一、问题的背景与提出“垃圾围城”是世界性难题，在今天的中国显得尤为突出。2012年全国城市生活垃圾清运量达到1.71亿吨，比2010年增长了1300万吨。数据显示，目前全国三分之二以上的城市面临“垃圾围城”问题，垃圾堆放累计侵占土地75万亩。因此，垃圾焚烧正逐步成为中国垃圾处理的主要手段之一。城市垃圾经过分类处理，剔除可回收垃圾和有害垃圾后将剩余垃圾在焚烧炉中焚烧处理，既可避免垃圾填埋侵占大量的土地，又可利用垃圾焚烧产生的能量进行发电等获得可观的经济效益。然而，由于政府监管不力、投资者目光短浅等多方面的原因，致使前些年各地建设的垃圾焚烧电厂在运营中出现了环境污染问题，给垃圾焚烧技术在我国的推广造成了很大阻力，许多城市的新建垃圾焚烧厂选址都出现因居民反对而难以落地的局面。事实上目前建厂选址尤为困难，国内各大城市目前均倾向于采用新型大型焚烧炉的焚烧厂取代分散的小型焚烧炉的举措。然而大型焚烧厂又存在需要考虑垃圾运输成本与道路建设成本等问题，因此对于不同城市来说，究竟该把大型焚烧厂的建设规模控制在什么水平，这是一个值得研究的课题。在垃圾焚烧厂运行监管方面，目前主要是在垃圾焚烧厂内进行测量监控，缺少从周边环境视角出发的外围动态监控，因而难以形成为民众所信服的全方位垃圾焚烧厂环境监控体系。因此，从构建环境动态监控体系、并根据潜在污染风险对周围居民进行合理经济补偿的需求出发，综合考虑垃圾焚烧厂对周围带来环境污染以及其他危害的多种因素（例如，焚烧炉的污染物排放量、居住点离开垃圾焚烧厂的距离、风力和风向及降雨等气象条件、地形地貌以及建筑物的遮挡程度等等），在进行科学定量分析的基础上，确立一套可行的垃圾焚烧厂环境影响动态监控评估方法，并针对潜在环境风险制定出合理的经济补偿方案是迫切需要的。二、基本假设1. 该厂是该地区唯一一个高架污染点源。 2. 每个监控点的监控能力是相同的，即每个监控点对同一种污染物的监控效果相同。 3. 监测的成本应尽量低，即监控点要尽可能少。 4. 假设垃圾焚烧厂周围居民风险承担经济补偿只与总污染程度有关。三、主要变量符号说明为了便于描述问题，我们用一些符号来代替问题中涉及的一些基本变量，如表1所示。表1：主要变量符号说明一览表政府用于赔偿总金额每个区域的赔偿金额每个区域的平均污染程度，第个监测点的颗粒物浓度第个监测点的浓度第个监测点的浓度每个监测点的总污染程度第个监测点故障时的浓度【注】其余没有列出的符号，我们将在文章第一次出现时给出具体说明四、问题的分析针对问题一，环境检测体系包括污染浓度函数的建立与监测点的设置。为了建立对垃圾焚烧厂周边的环境动态监控体系,首先需要建立污染物排放扩散模型,实现对垃圾焚烧厂外围污染物浓度的监控。通过分析，污染物扩散应考虑源强、风力风速及高度等因素的影响。通过附件4将风向分为8个方向，并能求出每个风向的年平均风速。本文将总污染程度定义颗粒物、的浓度总和，因此需要确定颗粒物、对人体的危害程度以用层次分析法确定权重。最后利用区域周围监测点的总污染物浓度确定该区域的平均污染浓度，并以此为依据对该区域居民进行合理的赔偿。针对问题二：由于焚烧炉的除尘装置出现故障或损坏，各污染物的指标将严重超标，因此，各个监测区的污染物浓度也会大大增加。利用SPSS画出附件2中的各污染物三月里焚烧炉排放浓度散点图。观察后发现，只有颗粒物的浓度在某些天内超标，而浓度和浓度均符合在国家标准。因此本文由颗粒物的故障条件来类比故障时的浓度和浓度。定义颗粒物浓度为25以上为故障浓度，因此对附件2中所提供的超过25的数据求平均，得到30，超过国家标准50%。则故障时的各监测点的浓度为：。最后，按照问题一中给出的补偿金额计算方式，得出新的补偿方案。五、问题一的求解5.1高斯模型5.1.1坐标系高斯模型的坐标系如图1所示，原点为排放点（若为高架源，原点为排放点在地面的投影），x轴正向为风速方向，y轴在水平面上垂直于x轴，正向在x轴的左侧，z轴垂直于水平面xoy，向上为正向。在此坐标系下烟流中心线或烟流中心线在xoy面的投影与x轴重合。图1：高斯模式的坐标系5.1.2模型假设(1)污染物的浓度在y、z轴上的分布是高斯分布（正态分布）的；(2)污染源的源强是连续且均匀的，初始时刻云团内部的浓度、温度呈均匀分布；(3)扩散过程中不考虑云团内部温度的变化，忽略热传递、热对流及热辐射；(4)泄漏气体是理想气体，遵守理想气体状态方程；(5)在水平方向，大气扩散系数呈各向同性；(6)取x轴为平均风速方向，整个扩散过程中风速的大小、方向保持不变，不随地点、时间变化而变化；(7)地面对泄漏气体起全反射作用，不发生吸收或吸附作用；(8)整个过程中，泄漏气体不发生沉降、分解，不发生任何化学反应等。5.1.3模型公式推导由正态分布假设可以导出下风向任意一点X（x,y,z）处泄漏气体浓度的函数为： （1） 由概率统计理论可以写出方差的表达式为： （2）由假设可以写出源强的积分公式： （3）式中：、为泄漏气体在y、z方向分布的标准差，单位为 m；X（x,y,z）为任一点处泄漏气体的浓度，单位为 kg/m；u为平均风速，单位为 m/s；Q为源强（即泄漏速度），单位为 kg/s；将（1）式代入（2）式，积分可得： （4）将（1）式和（4）式代入（3）式，积分可得： （5）再将（4）式和（5）式代入（1）式，可得： （6）上式为无界空间连续点源扩散的高斯模型公式，然而在实际中，由于地面的存在，烟羽的扩散是有界的。根据假设可以把地面看作一个镜面，对泄漏气体起全反射作用，并采用像源法处理，原理如图2所示。可以把任一点p处的浓度看作两部分的贡献之和：一部分是不存在地面时所造成的泄漏物浓度；一部分是由于地面反射作用增加的泄漏物浓度。该处的泄漏物浓度即相当于不存在地面时由位于（0,0,H）的实源和位于（0,0,-H）的像源在P点处所造成的泄漏物浓度之和。图2：像源法原理示意图其中，实源的贡献为： （7）其中，像源的贡献为： （8）则该处的实际浓度为： （9）由以上条件公式可得到高架连续点源扩散的高斯烟羽模型公式为： （10）其中：X（x,y,z）为下风向x米、横向y米、地面上方z米处的扩散的气体浓度，单位为kg/m；Q为源强（即源释放速率），单位为kg/s；u为平均风速，单位为m/s；为水平扩散参数，单位为m；为垂直扩散参数，单位为 m；t为泄漏后是时间，单位为s；H为泄漏源有效高度，单位为m；y为横向距离，单位为m；z为垂直方向距离，单位为m。如：式中，令z=0，即可得到地面气体浓度计算公式： (11)令y=z=0，即可得到地面轴线气体浓度计算公式： （12）其中，X（x,y,0）为下风向x米、横向y米处的地面扩散气体浓度，单位为kg/m；若令y=0，则可以得到下风向中心线上的浓度分布。5.1.4泄漏源有效高度以上式中的泄漏源有效高度是指泄漏气体形成的气云基本上变成水平状的时候气云中心的离地高度。实际上，泄漏源有效高度就等于泄漏源几何高度加泄漏烟云抬升高度。影响烟云抬升高度的因素有很多，主要包括：泄漏气体的初始速度和方向、初始温度、泄漏口直径、环境风速及风速岁高度的变化率、环境温度及大气稳定度。有效源高：。其中，为泄漏源几何架高，为烟云抬升高度。实验表明，泄漏源抬升高度可以用下面公式近似计算：其中，是烟囱出气口速度，单位为m/s；是出口直径，单位为 m；V为环境风速，单位为m/s。经查阅资料，出气口速度可以确定为15m/s，出口直径为4 m。上式是20世纪70年代末、80年代初，Wilson根据管道破裂泄漏实验所得的实验公式。实验时气体的喷射方向与风速垂直并且气体喷射路径上无障碍物。5.1.5扩散系数的选取扩散系数、的大小与大气湍流结构、离地高度、地面粗糙度、泄漏持续时间、抽样时间间隔、风速以及离开泄漏源的距离等因素有关。大气的湍流结构和风速在大气稳定度中考虑。大气稳定度由10米高度以上的风速、白天的太阳辐射或夜间的云量等参数决定。按照Pasquill的分类方法，随着气象条件稳定性的增加，大气稳定度可以分为A、B、C、D、E、F六类。从 A 到 F 的大气稳定性 依次增强。A 类表示气象条件极不稳定，B 类表示气象条件中度稳定，C 类表示气象条件略不稳定。大气稳定度具体分类方法见下表2、表3。表2：Pasquill大气稳定度的确定地面风速（m/s）白天日照夜间条件强中等弱阴天且云层薄，或低空云量为4/8天空云量为3/8AA-BBCCA-BBB-CC-DDBCCDDEDDDFEDD表3：日照强度的确定天空云层的情况和ABCDEF 通过对深圳市气象局对 10m 高度以上的风速、白天的太阳辐射、或夜间的云量等参数进行估计，确定该区域日照为强，并确定大气稳定度，即可以测量出该地区污染物扩散程度 ， 。详见下表。表4：垃圾场地区大气稳定度风向风速（m/s）大气稳定度东1.5A南3.1B西3B北1.9A东北1.8A东南1.6A西南2.9B西北2.6B5.2层析分析模型层次分析法(AHP)基于数学和心理学，是一种组织和分析复杂决策的结构化技术,。它是由Thomas L. Saaty在1970年提出的，自那时以来,一直得到广泛的研究和细化。层次分析法的基本原理是排序的原理，即最终将各方法（或措施）排出优劣次序，作为决策的依据。具体可描述为：层次分析法首先将决策的问题看作受多种因素影响的大系统，这些相互关联、相互制约的因素可以按照它们之间的隶属关系排成从高到低的若干层次，叫做构造递阶层次结构。然后请专家对各因素两两比较重要性，再利用数学方法，对各因素层层排序，最后对排序结果进行分析，辅助进行决策。它的主要特点是定性与定量分析相结合，将人的主观判断用数量形式表达出来并进行科学处理，因此，更能适合复杂的社会科学领域的情况，较准确地反映社会科学领域的问题。同时，这一方法虽然有深刻的理论基础，但表现形式非常简单，容易被人理解、接受，因此，这一方法得到了较为广泛的应用。5.2.1建立递阶层次结构颗粒物浓度各监测点浓度浓度浓度监测点监测点目标层:准则层:方案层:5.2.2建立两两比较矩阵通过一个和另一个元素的比较来评估它的各种因素。在得到比较矩阵之前，先给出表5。表 5: 重要性标度含义表重要性标度含义1表示两个元素相比，具有同等重要性3表示两个元素相比，前者比后者稍重要5表示两个元素相比，前者比后者明显重要7表示两个元素相比，前者比后者强烈重要9表示两个元素相比，前者比后者极端重要2,4,6,8表示上述判断的中间值根据表3，通过比较一个元素和另一个元素的重要性来评估各种因素。每个因素的重要性由来自某高校大学生物化学系教授和另一高校生科院教授共同判断，得到判断矩阵M。如果用表来呈现判断矩阵，结果如下表6所示。表 6: 重要性的比较颗粒物颗粒物1115.2.3计算每个指标的权重通过判断矩阵M，用matlab计算最大特征值，最大特征值的向量，得到每个因素的权重如下表7所示。表7:每个因素的权重颗粒物权重0.17650.47060.35295.2.4一致性检验假设一致性指标为C.I.，有。随机一致性指标R.I.可以在下表8中找到。判断矩阵的阶数是3。因此本文选择。表 8: 平均随机一致性指标R.I.表矩阵阶数12345678000.520.891.121.261.361.41因此，由公式，可得，即结果可被接受。5.3模型求解按照上述模型利用高斯公式计算任意坐标点的浓度，本文在以垃圾场为中心，边长为10千米的正方形中选取500个点，利用matlab编程计算每个风向的每种污染物的浓度，并画出图像。下图为风向取北时，随距离垃圾场远近的不同，浓度变化情况的分布图。图3：风向为北时浓度岁距离变化图由上图可知，在距离垃圾场5千米以外的区域基本不受废气污染的影响。因此为了全面准确地监控垃圾场附近居民住所的废气污染程度，本文以垃圾场为圆心，5千米为半径画出废气污染范围，并将此圆按8个风向分为8个扇形区域，再以1千米为间距，在每个风向上取5个监测点，即40个监测点（交点），40个区域，见图4，图5。 图4：扇形区域图 图5：各区域划分图利用matlab计算出上述40个监测点的各污染物浓度，结果如表9，表10，表11所示。表9：颗粒物污染浓度监测点（西北）12345各污染物浓度和(mg/m3)0.0248 0.0146 0.0063 0.0035 0.0022 监测点（北）678910各污染物浓度和(mg/m3)0.0287 0.0122 0.0050 0.0027 0.0017 监测点（东北）1112131415各污染物浓度和(mg/m3)0.0288 0.0136 0.0056 0.0031 0.0019 监测点（东）1617181920各污染物浓度和(mg/m3)0.0289 0.0142 0.0059 0.0032 0.0021 监测点（东南）2122232425各污染物浓度和(mg/m3)0.0250 0.0144 0.0062 0.0034 0.0022 监测点（南）2627282930各污染物浓度和(mg/m3)0.0211 0.0154 0.0068 0.0038 0.0024 监测点（西南）3132333435各污染物浓度和(mg/m3)0.0212 0.0156 0.0069 0.0038 0.0024 监测点（西）3637383940各污染物浓度和(mg/m3)0.0211 0.0163 0.0073 0.0041 0.0026 表10：污染浓度监测点（西北）12345各污染物浓度和(mg/m3)0.0993 0.0585 0.0253 0.0140 0.0089 监测点（北）678910各污染物浓度和(mg/m3)0.1148 0.0490 0.0201 0.0109 0.0069 监测点（东北）1112131415各污染物浓度和(mg/m3)0.1152 0.0542 0.0225 0.0123 0.0078 监测点（东）1617181920各污染物浓度和(mg/m3)0.1154 0.0569 0.0237 0.0130 0.0082 监测点（东南）2122232425各污染物浓度和(mg/m3)0.1001 0.0577 0.0247 0.0136 0.0086 监测点（南）2627282930各污染物浓度和(mg/m3)0.0845 0.0617 0.0272 0.0151 0.0096 监测点（西南）3132333435各污染物浓度和(mg/m3)0.0846 0.0625 0.0277 0.0153 0.0098 监测点（西）3637383940各污染物浓度和(mg/m3)0.0844 0.0654 0.0291 0.0162 0.0104 表11：污染浓度监测点（西北）12345各污染物浓度和(mg/m3)0.3104 0.1827 0.0791 0.0437 0.0278 监测点（北）678910各污染物浓度和(mg/m3)0.3588 0.1530 0.0628 0.0341 0.0216 监测点（东北）1112131415各污染物浓度和(mg/m3)0.3599 0.1695 0.0704 0.0384 0.0243 监测点（东）1617181920各污染物浓度和(mg/m3)0.3607 0.1777 0.0742 0.0405 0.0257 监测点（东南）2122232425各污染物浓度和(mg/m3)0.3127 0.1803 0.0772 0.0425 0.0270 监测点（南）2627282930各污染物浓度和(mg/m3)0.2642 0.1929 0.0851 0.0472 0.0301 监测点（西南）3132333435各污染物浓度和(mg/m3)0.2644 0.1954 0.0864 0.0479 0.0306 监测点（西）3637383940各污染物浓度和(mg/m3)0.2638 0.2043 0.0911 0.0506 0.0324 用层次分析法计算每个监测点的总污染浓度为：。由于本文是按照区域来对居民进行赔偿，因此需要根据监测点来求出每个区域的平均污染浓度。在此，作出以下规定，每个区域的4个顶点的总污染浓度取平均即为该区域的平均污染浓度。下表即为40个区域的污染程度。表12：各区域平均污染浓度区域（西北）12345各污染物浓度和(mg/m3)0.1606730.0945470.0409550.0225990.014403区域（北）678910各污染物浓度和(mg/m3)0.1857060.079210.0324880.0176520.011179区域（东北）1112131415各污染物浓度和(mg/m3)0.1863070.0877340.0364410.0198690.012601区域（东）1617181920各污染物浓度和(mg/m3)0.18670.0919720.03840.0209670.013305区域（东南）2122232425各污染物浓度和(mg/m3)0.1618560.0933450.0399720.0219810.013989区域（南）2627282930各污染物浓度和(mg/m3)0.1367360.0998510.0440720.0244350.015603区域（西南）3132333435各污染物浓度和(mg/m3)0.1368640.1011660.0447350.0248150.01585区域（西）3637383940各污染物浓度和(mg/m3)0.1365580.1057580.0471440.026210.016756经观察google地图后发现，某些区域为无人居住区，因此，对这些区域不考虑予以赔偿（图中用阴影表示的区域为无人居住区）。假设政府共有赔偿金W元，则每个区域可按该区域的平均污染浓度得到相应的赔偿金。分发到每个区域的赔偿金按人口平均分配。下表为每个区域所应得的赔偿。表13：各区域赔偿金额区域（西北）12345赔偿金额（元）0.0355W0.0154W0.0085W0.0054W区域（北）678910赔偿金额（元）0.0698W0.0298W0.0066W区域（东北）1112131415赔偿金额（元）0.07W0.033W0.0137W0.0075W区域（东）1617181920赔偿金额（元）0.0144W0.0079W区域（东南）2122232425赔偿金额（元）0.0083W0.0053W区域（南）2627282930赔偿金额（元）0.0092W0.0059W区域（西南）3132333435赔偿金额（元）0.038W0.0168W0.0093W0.006W区域（西）3637383940赔偿金额（元）六、问题二的求解由于焚烧炉的除尘装置出现故障或损坏，各污染物的指标将严重超标，因此，各个监测区的污染物浓度也会大大增加，甚至会因气体重力、浓度等各种因素致 使其他地区被污染。利用SPSS画出附件二中的各污染物三月里焚烧炉排放浓度波动图，如图6，图7，图8。 图6：颗粒物浓度随时间变化散点图 图7：浓度随时间变化散点图图8：浓度随时间变化散点图由上述三幅图可以看出，只有颗粒物的浓度在某些天内超标，而浓度和浓度均符合在国家标准。因此本文由颗粒物的故障条件来类比故障时的浓度和浓度。定义颗粒物浓度为25以上为故障浓度，因此对附件2中所提供的超过25的数据求平均，得到30，超过国家标准50%。则故障时的各监测点的浓度为：最后，按照问题一中给出的补偿金额计算方式，得出新的补偿方案。七、模型的检验和评价检验： 我们采用相对误差值来对模型进行检验。首先检验烟口平均环境风速，研究的变化对浓度的影响，可以考虑如果的相对变化，则发生的相对变化是的多少倍，即定义对的灵敏度为：如果增加1%，则变化的百分比是1%的倍。如果很小，则对不灵敏；反之则对灵敏，的微小变化会带来较大的变化。其他参数检验同理可得。优点：1. 输入的参数易于从常规气象观测中获取，简便易行；计算费用低，并且可以到达一定精度。2. 可以把比较多的变量因素简化处理从而简单方便来设计监测体系。3. 全面性。不仅考虑纵切和横切的情况，还考虑了单面和双面打印文件的碎纸片拼接复原问题，从而使模型的应用范围更广泛。4. 本文所建模型与实际紧密联系，可以应用于现实生活中，通用性强，具有现实意义。缺点：1. 忽略了地形降雨等很多因素影响，这种比较与实际有点误差差距。八、参考文献1 程勇，于林，姚安林，采用高斯模型分析输气管道泄露后气体的扩散J， 内蒙古石油化工，2010,（16）：49-512 杨洪,林大伟,张宏礼,MATLA 中三个适于数据拟合命令的比较分析J,大庆 师范学院学报,2010,(6)：62-653 方成贤，董兴玲，龚光辉，姜宗顺，高虹，江晶，垃圾焚烧厂的环境补偿机 制探讨J，环境工程, ,2009,（27）（增刊）：603-605.4 王文，于雷，裴文文，杨方，基于高斯线源模式的主要尾气扩散模型综述J， 交通环保，2004,(05)期：8-105 徐田伟，矿产资源开发生态补偿机制初探J，环境保护与循环经济， 2009,(07)：58-616 孔波，对于环境监测数据的具体分析，科技资讯J，2011,(24):34-367 徐晓敏, 层次分析法的运用J,统计与决策, 2008,(01)：11-13附录附录一：高斯模型syms Sy Sz H%Q源强（即泄漏速度），u为平均风速，H烟囱高度,Vs排烟速度Q=15*2*pi*4*0.08;u=1.9;Vs=15;H=80+2.4*Vs*4/u;%A扩散系数x,y=meshgrid(-5000:20:5000);z=zeros(size(x);Sy=0.22*x./sqrt(1+0.0001*x);Sz=0.2*x;Bei=Q./(2*pi*u*Sy.*Sz).*exp(-y.2./(2*Sy.2).*(exp(-(z-H).2./(2*Sz.2)+exp(-(z+H).2./(2*Sz.2);Bei(:,1:251)=0;mesh(x,y,Bei)view(0 -1 0)syms Sy Sz H%Q源强（即泄漏速度），u为平均风速，H烟囱高度,Vs排烟速度Q=15*2*pi*4*0.02;u=3.1;Vs=15;H=80+2.4*Vs*4/u;%B扩散系数x,y=meshgrid(-5000:20:5000);Sy=0.16*x./sqrt(1+0.0001*x);Sz=0.12*x;z=zeros(size(x);Bei=Q./(2*pi*u*Sy.*Sz).*exp(-y.2./(2*Sy.2).*(exp(-(z-H).2./(2*Sz.2)+exp(-(z+H).2./(2*Sz.2);Bei(:,1:250)=0;mesh(x,y,Bei)附录二：层次分析法A=1 3/8 1/2 8/3 1 4/3 2 3/4 1;v,d=eig(A)v=v(:,1)/sum(v(:,1)CI=(max(max(d)-3)/2CR=CI/0.52附录三：单方向污染物浓度计算程序function P=pollution(d,c,x,y)%Q源强（即泄漏速度），u为平均风速，H烟囱高度，t时间，c浓度Q=15*2*pi*4*c;Vs=15;z=zeros(size(x);if d=1%北 u=1.9;H=80+2.4*Vs*4/u;Sy=0.22*x./sqrt(1+0.0001*x);Sz=0.2*x;P=Q./(2*pi*u*Sy.*Sz).*exp(-y.2./(2*Sy.2).*(exp(-(z-H).2./(2*Sz.2)+exp(-(z+H).2./(2*Sz.2); elseif d=2%东u=1.5;H=80+2.4*Vs*4/u;Sy=0.22*x./sqrt(1+0.0001*x);Sz=0.2*x;P=Q./(2*pi*u*Sy.*Sz).*exp(-y.2./(2*Sy.2).*(exp(-(z-H).2./(2*Sz.2)+exp(-(z+H).2./(2*Sz.2); elseif d=3%东北u=1.8;H=80+2.4*Vs*4/u;Sy=0.22*x./sqrt(1+0.0001*x);Sz=0.2*x;P=Q./(2*pi*u*Sy.*Sz).*exp(-y.2./(2*Sy.2).*(exp(-(z-H).2./(2*Sz.2)+exp(-(z+H).2./(2*Sz.2); elseif d=4%东南u=1.6;H=80+2.4*Vs*4/u;Sy=0.22*x./sqrt(1+0.0001*x);Sz=0.2*x;P=Q./(2*pi*u*Sy.*Sz).*exp(-y.2./(2*Sy.2).*(exp(-(z-H).2./(2*Sz.2)+exp(-(z+H).2./(2*Sz.2); else