基于WRF模型的计算机程序化计算大气环境容量的方法与流程

本发明涉及环境影响评价
技术领域：
，尤其是涉及一种基于wrf模型的计算机程序化计算大气环境容量的方法。
背景技术：
：环境容量(environmentcapacity)又称环境负载容量、地球环境承载容量或负荷量，是在人类生存和自然生态系统不致受害的前提下，某一环境所能容纳的污染物的最大负荷量。大气环境容量主要与大气混合层厚度、风速、沉降速率、降水、污染物标准浓度等因素相关，以往的大气环境容量研究和实际应用中大气环境容量计算主要采用a值法，依据国家颁布的相关标准和及再次基础上进行改良的公式进行计算，但由于采用人工估算，公式中相关参数的取值多采用较大地理范围内的年平均值和经验数值，因此计算结果存在时空精度差、经验依赖性较大的问题。本发明以wrf模型为例，该模型可利用在互联网上公开的基础数据经计算输出精确到小时的各种气象参数，且计算地理范围可精确到数平方公里，通过编制计算机程序实现自动化计算，最终可实现任意地理范围内时间精确到小时的大气环境容量计算；因此，本发明提供一种基于wrf模型的计算机程序化计算大气环境容量的方法。技术实现要素：针对上述情况，为克服现有技术之缺陷，本发明提供一种基于wrf模型的计算机程序化计算大气环境容量的方法。有效的解决了在计算大气环境容量时采用人工估算，公式中相关参数采用较大地理范围内年平均值和经验数值所带来的计算结果时空精度差、经验依赖性大的问题。基于wrf模型的计算机程序化计算大气环境容量的方法，其特征在于，计算步骤如下：步骤一、计算某一选定区域w污染物的大气环境容量，选定所要计算区域中心坐标为北纬a，东经b，边长为x*y的范围，对于所要计算区域进行网格化处理，将其格网化为x*y个网格，选取计算范围以单个网格为例进行计算，计算时间长度为k年l月m日0时-23时；根据公式δ＝[0.006918-0.399912cosθ0+0.070257sinθ0-0.006758cos2θ0+0.000907sin2θ0-0.002697cos3θ0+0.001480sin3θ0]×180/π式中：θ0＝360dn/365,deg；δ—太阳倾角，deg；dn—一年中日期序数，0，1，2……，364得出所要计算区域k年l月m日的太阳倾角δ步骤二、将步骤一中求得的太阳倾角代入公式h0＝arcsin[sinφsinδ+cosφcosδcos(15t+λ-300)]式中：h0—太阳高度角，deg；φ—当地纬度，deg；t—北京时间，h；λ—当地纬度，deg；计算网格点中心坐标为：北纬a，东经b，得出0时-23时各个时间段内对应的太阳高度角h0，根据太阳高度角可知：k年l月m日昼间为c-d时，夜间为g-h时，e-f时；步骤三、查询地面气象sam文件得出网格点中心坐标为：北纬a，东经b，0时-23时各个时间段内对应的总云量/低云量；步骤四、由太阳高度角h0、总云量/低云量、以及昼夜情况查询太阳辐射等级表得出所要计算网格点中心坐标为：北纬a，东经b，0时-23时各个时间段内对应的太阳辐射等级；步骤五、将以上数据参数由wrf模型运算完成输出所要计算网格点中心坐标为：北纬a，东经b，0时-23时各个时间段内的风速u、干沉降速率ud、降水速率r；步骤六、由太阳辐射等级与风速的对应关系，查询大气稳定度等级表得出网格点中心坐标为：北纬a，东经b，0时-23时各个时间段内的大气稳定度等级；步骤七、结合地区序号表确定网格点中心坐标为：北纬a，东经b，所处的地区序号，并且结合大气稳定度等级查询我国不同地区和在不同大气稳定度等级下所对应的as/bs数值，求得网格点中心坐标为：北纬a，东经b，0时-23时各个时间段内对应的as/bs数值；步骤八、综合网格点中心坐标为：北纬a，东经b，0时-23时各个时间段内对应的风速以及0时-23时各个时间段内对应的as/bs数值，根据公式：在大气稳定度为a、b、c和d级时：在大气稳定度为e和f级时：式中：h—混合层厚度，m；u—10m高度上平均风速，m/s；大于6m/s时取为6m/s；as，bs—混合层系数；f—地转参数；ω—地转角速度，由《制定地方大气污染物排放标准的技术方法》gb/t13201-91里规定取为7.29×10-5rad/s；—地理纬度；求得网格点中心坐标为：北纬a，东经b，地转参数f以及0时-23时各个时间段内对应的大气混合层厚度；步骤九、大气环境容量q的表达式表示为：其中：h:大气混合层厚度(m)u：风速(m/s)ud:干沉降速度(m/s)r：降水速率(mm/a)ωr：清洗比，由《制定地方大气污染物排放标准的技术方法》gb/t13201-91里规定取值为1.9×10-5cs：污染物标准浓度(mg/m3)s：区域面积(km2)q：大气污染物容量(104t/a)区域内存在不同大气功能区时各功能区ai计算公式为：ai＝aci式中ci为大气环境质量标准规定的第i功能区类别对应的年日平均浓度限值(mg/m3)qi＝ai×csi×si/s1/2式中：csi：大气环境质量标准规定的第i功能区类别对应的年平均浓度限值(mg/m3)si：第i功能区面积ci、csi由查询气象文件及相关国家标准获得；计算得出网格点中心坐标为：北纬a，东经b，0时-23时各个时间段内对应的qi值，由于公式中qi值量纲为104t/a，乘以系数1.14最终得出小时值qhi对qhi进行求和得出网格点中心坐标为：北纬a，东经b，0时-23时内w污染物的大气环境容量。(1)本发明以wrf模型为基础，由于wrf模型可以利用互联网上公开的基础数据经计算输出精确到小时的各种气象参数，且计算地理范围可精确到数平方公里，通过编制计算机程序从wrf计算结果文件中提取计算大气环境容量所需的参数实现自动化计算，最终可实现任意地理范围内时间精确到小时的大气环境容量计算；(2)本发明可以快速准确的计算所关注区域的大气环境容量，为企业生产和大气环境相关政策的制定提供科学依据。附图说明图1为本发明所选计算区域示意图。具体实施方式有关本发明的前述及其他技术内容、特点与功效，在以下配合参考附图至图对实施例的详细说明中，将可清楚的呈现。以下实施例中所提到的结构内容，均是以说明书附图为参考。下面将参照附图描述本发明的各示例性的实施例。实施例1，基于wrf模型的计算机程序化计算大气环境容量的方法，其特征在于，计算步骤如下：步骤一、计算某一选定区域so2污染物的大气环境容量，选定所要计算区域中心坐标为北纬33.3°，东经113.5°，边长为50km×50km的范围，对于所要计算区域进行网格化处理，将其格网化为50×50个网格，选取计算范围以单个网格为例进行计算，计算时间长度为2012年1月1日0时-23时；计算区域示意图为选取计算网格根据公式δ＝[0.006918-0.399912cosθ0+0.070257sinθ0-0.006758cos2θ0+0.000907sin2θ0-0.002697cos3θ0+0.001480sin3θ0]×180/π式中：θ0＝360dn/365,deg；δ—太阳倾角，deg；dn—一年中日期序数，0，1，2……，364dn＝0,θ0＝0，计算得出δ＝-23.06步骤二、将步骤一中求得的太阳倾角代入公式h0＝arcsin[sinφsinδ+cosφcosδcos(15t+λ-300)]式中：h0—太阳高度角，deg；φ—当地纬度，deg；t—北京时间，h；λ—当地纬度，deg；计算网格点中心坐标为：北纬33.3°，东经113.5°，得出0时-23时各个时间段内对应的太阳高度角h0t012345678910111213141516h0-78-77-67-55-42-30-18-651524303333302314t17181920212223h04-7-19-31-44-56-68根据太阳高度角可知：2012年1月1日昼间为8-17时，夜间为0-7时，18-23时；步骤三、查询地面气象sam文件得出网格点中心坐标为：北纬33.3°，东经113.5°，0时-23时各个时间段内对应的总云量/低云量t012345678910111213141516总云55555555555555555低云55555555555555555t17181920212223总云5555555低云5555555注：云量(全天空十分制)观测规则见中央气象局编定的《地面气象观测规范》第3.3节。步骤四、由太阳高度角h0、总云量/低云量、以及昼夜情况查询太阳辐射等级表得出所要计算网格点中心坐标为：北纬33.3°，东经113.5°，0时-23时各个时间段内对应的太阳辐射等级；太阳辐射等级网格点中心坐标为：北纬33.3°，东经113.5°，0时-23时各个时间段内对应的太阳辐射等级：t012345678910111213141516b100000000000000000t17181920212223b10000000步骤五、将以上数据参数由wrf模型运算完成输出所要计算网格点中心坐标为：北纬33.3°，东经113.5°，0时-23时各个时间段内的风速u、干沉降速率ud、降水速率rt012345678910111213141516u2.43.03.53.60.54.54.95.04.64.85.35.20.54.94.33.52.5t17181920212223u1.81.41.00.70.40.30.6t012345678910111213141516r00000000000000000t17181920212223r0000000步骤六、由太阳辐射等级与风速的对应关系，查询大气稳定度等级表得出网格点中心坐标为：北纬33.3°，东经113.5°，0时-23时各个时间段内的大气稳定度等级大气稳定度等级表注：1)地面风速(m/s)系指离地面10m高度处10分钟平均风速网格点中心坐标为：北纬33.3°，东经113.5°，0时-23时各个时间段内的大气稳定度等级：t012345678910111213141516b2dddddddddddddddddt17181920212223b2ddddddd步骤七、结合地区序号表确定网格点中心坐标为：北纬33.3°，东经113.5°，所处的地区序号，并且结合大气稳定度等级查询我国不同地区和在不同大气稳定度等级下所对应的as/bs数值，求得网格点中心坐标为：北纬33.3°，东经113.5°，0时-23时各个时间段内对应的as/bs数值；由表1确定地区序号为“3”表1结合地区序号和大气稳定度等级b2，由表e1查出as,bs值：表e1我国各地区as和bs值网格点中心坐标为：北纬33.3°，东经113.5°，0时-23时各个时间段内对应的as/bs数值：t012345678910111213141516as0.0190.0190.0190.0190.0190.0190.0190.0190.0190.0190.0190.0190.0190.0190.0190.0190.019t17181920212223as0.0190.0190.0190.0190.0190.0190.019步骤八、综合网格点中心坐标为：北纬33.3°，东经113.5°，0时-23时各个时间段内对应的风速以及0时-23时各个时间段内对应的as/bs数值，根据公式：在大气稳定度为a、b、c和d级时：在大气稳定度为e和f级时：式中：h—混合层厚度，m；u—10m高度上平均风速，m/s；大于6m/s时取为6m/s；as，bs—混合层系数，见表e1；f—地转参数；ω—地转角速度，由《制定地方大气污染物排放标准的技术方法》gb/t13201-91里规定取值为7.29×10-5rad/s；—地理纬度；得出网格点中心坐标为：北纬33.3°，东经113.5°的地转参数为f＝7.96×10-5；求得网格点中心坐标为：北纬33.3°，东经113.5°，0时-23时各个时间段内对应的大气混合层厚度：t012345678910111213141516h573716835859119107411701193109811461265124111911701026835597t17181920212223h4303342391679572143综合以上结果，a值计算公式中所需参数已全部求得，计算所需的各参数分别为：步骤九、大气环境容量q的表达式表示为：其中：h:大气混合层厚度(m)u：风速(m/s)ud:干沉降速度(m/s)r：降水速率(mm/a)ωr：清洗比，由《制定地方大气污染物排放标准的技术方法》gb/t13201-91里规定取值为1.9×10-5；cs：污染物标准浓度(mg/m3)s：区域面积(km2)q：大气污染物容量(104t/a)由此计算出2012年1月1日0-23时逐时a值：t012345678910111213141516a3.845.998.168.640.1713.4915.9916.6514.1015.3518.7118.010.1715.9912.318.164.16t17181920212223a2.161.310.670.330.110.060.24区域内存在不同大气功能区时各功能区ai计算公式为：ai＝aci式中ci为大气环境质量标准规定的第i功能区类别对应的年日平均浓度限值(mg/m3)qi＝ai×csi×si/s1/2式中：csi：大气环境质量标准规定的第i功能区类别对应的年平均浓度限值(mg/m3)si：第i功能区面积ci、csi由查询气象文件及相关国家标准获得；根据环境空气质量标准(gb3095-2012)中相关规定：环境空气质量功能区的分类和标准分级环境空气质量功能区分类：一类区为自然保护区、风景名胜区和其它需要特殊保护的地区；二类区为居住区、商业交通居民混合区、文化区、工业区和农村地区；环境空气质量标准分级：一类区执行一级浓度限值二类区执行二级浓度限值各项污染物的浓度限值：设定计算网格点执行二级标准csi＝0.06mg/m3ci＝0.15mg/m3由公式qi＝ai×csi×si/s1/2ai＝aci由此计算出2012年1月1日0-23时逐时ai值：t012345678910111213141516ai0.580.901.221.300.032.022.402.502.112.302.812.700.032.401.851.220.62t17181920212223ai0.320.200.100.050.020.0090.04若出现u＝0公式中分母为0的情况，该时间点计算结果无效舍弃计算得出网格点中心坐标为：北纬33.3°，东经113.5°，0时-23时各个时间段内对应的qi值，由于公式中qi值量纲为104t/a(104吨/年)，乘以系数1.14最终得出小时值qhi对qhi进行求和得出网格点中心坐标为：北纬33.3°，东经113.5°，0时-23时内w污染物的大气环境容量；1.14为万吨/年与吨/小时的换算系数，由10000(吨)/8760(小时)计算得到；求解qi值，由于公式中qi值量纲为104t/a，乘以系数1.14最终得出小时值qhi(t/h)对于qhi进行求和，得出2012年1月1日一天时间内选定计算网格内so2的环境容量为0.038t。本发明以wrf模型为基础，由于wrf模型可以利用互联网上公开的基础数据经计算输出精确到小时的各种气象参数，且计算地理范围可精确到数平方公里，通过编制计算机程序从wrf计算结果文件中提取计算大气环境容量所需的参数实现自动化计算，最终可实现任意地理范围内时间精确到小时的大气环境容量计算；本发明可以快速准确的计算所关注区域的大气环境容量，为企业生产和大气环境相关政策的制定提供科学依据。以上所述的实施例并非对本发明的范围进行限定，在不脱离本发明设计构思的前提下，本领域所述技术人员对本发明的技术方案作出的各种变形和改进，均应纳入本发明的权利要求书确定的保护范围内。当前第1页12