近几年来，东北冬季重雾霾污染过程频发.东北地区大气污染正朝着区域性跨界污染发展, 随着冬季供暖期的到来，由于东北三省地区较为封闭的地形和相似的气候气象条件, 供暖期持续时间较长，燃煤量很大，冬季地表气温低，逆温天气频繁.烟霾污染又使垂向热力对流交换进一步减弱，逆温天气往往伴随着静风，此时城市犹如一个“箱体”，污染物在箱体内不断累积，形成了污染气象的恶性循环，致使大气污染物持续累积；而逆温层下的水平向风风速很小，几乎为零，污染物的衰减以及传输能力很差，加重了区域内污染的程度^[1-2].

哈尔滨市作为典型的东北城市，具有东北寒冷地区城市大气污染物排放的主要特征，即采暖期哈尔滨市供暖燃煤使用量增多，机动车保有量不断增加、采暖期冰雪路面使机动车行驶时间增加怠速时间长导致尾气污染加剧，农村秋收后秸秆回收利用率低、经济效益差，秸秆焚烧问题无法彻底根治^[3-5].国内学者针对区域大气污染物源强的计算往往通过利用排放因子法构建区域大气排放清单方式来开展相关研究，主要集中于京津翼、长三角和珠三角等地区.潘月云等^[6]根据收集的广东省人为源活动水平数据采用排放因子法建立该区域2010年人为源大气污染物排放清单.宁亚东等^[7]在考虑移动源大气污染物排放标准变化上利用经验公式法估算了道路排放源和非道路排放源的大气污染物排放量.虽然可以通过国内学者建立的区域排放清单以及经验公式法来获得主要人为源的污染物总量，但存在着受排放源类别不全，排放清单的不确定性，人力物力财力因素制约，以及企业上报数据不准确，存在偷排多排的情况等问题.

本研究将以哈尔滨为研究区域，以该市2016年9月—2016年12月起始日期为时间节点，采集哈尔滨市11个大气国控监测点位大气污染物小时浓度值，得到采暖期大气污染物时间分布特征.利用大气箱式模型，计算出满足箱式模型假设特征的时段大气污染物排放源强，跳过企业大气污染物排放数据监控与上报步骤，符合排污许可证制度理念要求，为源强估算提供新的思路与较为简易的新方法，为东北城市大气排放控制提供科学依据.

1 数据来源与研究方法 1.1 数据来源与研究范围

大气污染物浓度数据由黑龙江省环境监测中心站提供，气象数据与地理信息数据获取来自中国空气质量在线监测分析平台和谷歌地图信息共享.

本研究区域空间范围为哈尔滨市区，包括动力和平路、南岗学府路、呼兰师专、太平宏伟公园、岭北、平房东轻厂、松北商大、省农科院、道外承德广场、道里建国路、阿城会宁、香坊红旗大街11座大气监测点所在分区.时间范围为2016年9月—2016年12月.

1.2 研究方法 1.2.1 哈尔滨市基于ArcGis和mapinfo的泰森多边形的构建

泰森多边形是荷兰气象学家A.H.Thiessen提出的，通过离散雨量点的雨量计算区域平均雨量，即将相邻的离散雨量点连成三角形，再做出三角形各边的垂线，每个离散点周围被垂线围成一个多边形，用这个泰森多边形中唯一雨量观测的降雨数据来代表这个多边形范围内的降水量.在泰森多边形内，每个泰森多边形都对应着一个雨量观测点，这个雨量观测点的降水数据就代表对应的泰森多边形内的降水数据.除在雨量计算方面的应用，方琪璐等^[8]建立了应用泰森多边形法计算行人密度的算法，宋洁华等^[9]以Gis为平台利用泰森多边形分析临高县居民点的空间分布特征.将此方法归纳到本研究中，以哈尔滨市区域的大气国控监测点作为泰森多边形的离散点，将相邻离散大气监测点连成三角形，做垂线围成多边形.在哈尔滨市泰森多边形中，每个泰森多边形对应着一个大气监测点，此大气监测点的大气污染物监测浓度就代表此泰森多边形内的大气污染物浓度^[10].由此可以将哈尔滨市分成多个泰森多边形集合，分别计算得到每个泰森多边形中大气污染物的排放源强，由此计算得出研究区间哈尔滨市采暖期的排放总源强.哈尔滨市市区地势较为平坦、低洼，大气监测点分布较为分散，每个站点同一时段大气浓度相差较大，大气污染物时空分布不均，因此，采用泰森多边形法对研究区域进行处理适用于哈尔滨市大气源强的估算.

Mapinfo的绘图功能、地图显示、属性导入及图像编辑功能较为直观、方便，可以直接导入Excel表的地理坐标数据，但无法显式表达图形对象间联系的数据结构，在图形中每个对象都是独立存在独立表达，空间分析功能相对单一，数据结构简单，无拓扑关系. Arcgis具有强大的空间分析能力，拓扑分析功能十分强大，除了支持传统的点线面简单要素模型外，也引入了面向对象的空间数据模型(GeoDatabase)，可以利用这个模型来定义和操作不同用户或应用的具体模型.同时，Mapinfo中的格式转换工具箱可将Mapinfo所用的TAB文件转换为Arcgis所使用的SHP格式文件，因此，实现Mapinfo和Arcgis的联用从而提高地理信息的处理效率是可行的.

采用Mapinfo以及ArcGis构建、优化、裁剪得到哈尔滨市的泰森多边形^[11-12].哈尔滨市泰森多边形集合如图 1所示.

通过ArcGis中计算几何功能可以快速得到各个多边形的面积，计算各个多边形的面积与区域总面积的比即可得到大气排放源强的权重系数，如表 1所示.通过箱式模型计算得到的各个多边形区域的排放源强乘以权重系数进行加权平均，就可以计算得到该区域的大气平均排放源强.

1.2.2 模型所选定目标时段的确定

综合哈尔滨市2016年9月—12月气象条件可知，冬季风速小，降水量很低，逆温天气出现频繁，气象条件不利于大气污染物的平流运输，导致哈尔滨市大气污染具有时间性变化特征.不同月份环境大气中的污染物富集与扩散能力不同，呈现出污染扩散从难到易排序依次为11月份、12月份、10月份、9月份，一天污染状况从重到轻依次为傍晚、上午、下午、凌晨早上的特点.2016年哈尔滨市9月—12月AQI变化曲线如图 2所示.

哈尔滨市冬季采暖期逆温气象出现较为频繁，原因多为辐射逆温，在逆温过程中，贴地逆温持续时间长且逆温层厚度较厚.以空气质量最严重的2016年11月为例，从相对湿度角度来看，空气中水汽条件是形成大雾天气的关键因素，能见度随相对湿度的增加而下降，风力风速较小，地面湿度大，地面气压场弱，遂形成大雾^[13]；其次，大气处于稳定状态，垂直和水平方向均无流动，相对湿度较小，低层空气比较干燥致使空气中的污染物无法扩散，导致雾霾天气产生.从风速和温度角度来看，天气现象表现为烟雾转大雾，但从2016年11月4日开始，哈尔滨市此时弱的偏南风，几乎全部表现为静风，为污染最为严重的时期，在雾霾天气过程期间，大气能见度随风速的减弱而持续降低.哈尔滨市最高温度从3日起突然回暖，这个气温的回升递直接诱导了逆温天气的产生，为霾的形成创造了条件，这种形势不利于形成大雾而利于形成霾，天气状况随温度的逐日上升由大雾转为雾霾，导致11月4日哈尔滨市出现严重污染天气.到6日，随冷空气的入侵，哈尔滨上空的逆温层被破坏，风速也逐时加大，大气层结由稳定状态转变为不稳定状态，污染颗粒物流动加速，污染物在对流层中传输能力增强，雾霾天气逐渐结束，空气质量指数也逐渐降低^[14]，如图 3所示.以岭北研究区为例，在时间区间2016/11/4 15:00—2016/11/4 20:00中SO₂、NO_x和CO质量浓度为近平滑曲线上升，由第二节中逆温天气分析可知，此时气象条件基本满足大气箱式模型的假设，因此，选定此时段和大气污染物浓度数据作为模型的输入.如图 4所示.

1.2.3 基于箱式模型法的排放强度估算方法

单箱模型是计算一个区域或者城市大气质量最简单的模型，箱子的平面尺寸就是研究区域的平面，其高度是由地面计算的混合层高度h.其基本模型为

$ \frac{{{\rm{d}}C}}{{{\rm{d}}t}}lbh = ubh\left( {{C_0} - C} \right) + lbQ - kClbh. $

(1)

式中:C为箱内污染物浓度；l为箱的长度；b为箱的宽度；h为箱子的高度；C₀为初始条件下污染物的本底浓度；k为污染物的衰减速度常数；Q为污染源的源强；u为平均风速；t为时间坐标.

如果不考虑污染物的衰减，k=0，当污染物浓度稳定排放时，即视源强Q为定值，可以得到式(1)的解：

$ C = {C_0} + \frac{{Ql}}{{uh}}(1 - {{\rm{e}}^{ - \frac{{ut}}{l}}}). $

(2)

2016年11月4日哈尔滨市出现严重污染天气，空气质量为全供暖期最差，温度明显回暖，逆温严重，几乎全部表现为静风天气，外部污染源横向扩散速度过缓(静风条件下以分子扩散为主)，对哈尔滨市辖区空气污染贡献较小，即此时大气污染物的积累以本地排放源贡献为主，定义此日期为哈尔滨市采暖期的“特征日期”，在此日期可以把哈尔滨市视为一个箱体，箱内污染物混合均匀.因此，利用大气箱式模型反算哈尔滨市是可行的，此时模型可信度最高.

考虑到研究区域呈现静风天气的特点，区域内采暖期平均风速较小，可以先忽略大气污染物的区域间输送，认为箱内污染物浓度主要由区域本身的污染源排放引起，即C₀=0.则式(2)可简化为

$ C = \frac{{Ql}}{{uh}}(1 - {{\rm{e}}^{ - \frac{{ut}}{l}}}). $

(3)

由式(3)可得源强表达式：

$ Q = \frac{{Cuh}}{{l\left( {1 - {{\rm{e}}^{ - \frac{{ut}}{l}}}} \right)}}. $

(4)

根据环境影响评价技术导则，混合层高度H的计算公式如下：

当大气稳定度为A、B、C和D时

$ H = \frac{{{a_{\rm{s}}}{u_{10}}}}{f}. $

(5)

当大气稳定度为E和F时

$ H = {b_{\rm{s}}}\sqrt {{u_{10}}/f} , $

(6)

$ f = 2\mathit{\Omega }\;{\rm{sin}}\varphi . $

(7)

式中：a_s、b_s为混合层系数；u₁₀为10 m高度处平均风速(m/s)；f为地转参数，Ω为地转角频率，取为7.29× 10^-5rad/s; φ为地理纬度(°).

在计算顺风向l时，顺风向长度l可以设为研究区等效直径，若研究区域的面积为A，则该区域的等效直径为$l = {\rm{2}}\sqrt {\frac{A}{{\rm{ \mathsf{ π} }}}} $^[15].混合层高度计算时根据环境影响评价技术导则，重污染天气期间哈尔滨市大气稳定度普遍在E和F，气象条件处于静风天气，b_s分别为1.66和0.7，风速为0.3~2 m/s，采暖期哈尔滨市大气混合层采用采暖期经验数值范围91.5~ 216.9 m较为准确合理.

1.2.4 最优化估参法

最优化方法适合于具有单峰极值的非线性模型的参数估计，通常采用迭代法求解，主要通过以下步骤实现：1)需要设定参数的初值和允许迭代误差; 2)计算出目标函数的初值; 3)计算目标函数对参数的梯度; 4)计算参数的修正步长; 5)计算参数的修正值; 6)计算新的目标函数值并与旧目标值进行比较，若不符合要求，返回第3步重新迭代，直到迭代误差小于预定数值.

本研究中采用迭代法对模型参数进行估计^[16]，由于模型估计参数的准确度较为依赖初值，在估参前应对模型参数的大致范围有所了解.鉴于本研究为单变量的曲线拟合，可以使用Matlab中的单一变量曲线拟合工具箱.本模型中主要参数风速u和混合层高度的范围确定方法已在1.2.3中给出.

2 结果与分析 2.1 各区域源强的计算

最重污染时段即2016年11月4日各个区域源强计算结果如表 2所示.可以看出，大气污染排放源强于道外区、动力和平路、松北商大、省农科院等距离市中心较近的区域处较大，而阿城会宁及呼兰师专等距离市中心较远的区域大气污染排放源强较小，这与中重度污染时城区区域污染程度现状较为接近.

2.2 源强估算模式的验证 2.2.1 估算模式的适用性分析

以岭北为例，根据1.2.2方法，分别选取11月10日、11月13日、11月17日、12月11日和12月17日较重污染天气时段作为大气箱式模型输入，并根据当天的气象条件分别设置风速以及混合层高度.将最重污染时段11月4日所算得的源强估算值y，与其他5个较重污染天气所算得的源强估算值y₁′, y₂′, y₃′, y₄′, y₅′对本源强估算模式进行精确度检验，其中根据1.2.3计算得到11月10日、11月13日、11月17日、12月11日和12月17日的混合层高度范围为212.3~445.2 m，如表 3所示.

通过对岭北区源强不同时段的中值误差计算，可知中值误差e_0.5＜15%，算法基本有效.为进一步验证模型的适用性，根据从黑龙江省环境监测中心站获得的2014年—2015年、2015年—2016年哈尔滨市各个监测点冬季采暖期间重污染及较重污染天气下大气污染物的时均浓度值，确定2015年—2016年采暖期最重污染且符合模型输入时段, 同样应用上述方法对本年采暖期大气平均排放源强估算模式进行精确度检验.其中，y₁′, y₂′, y₃′, y₄′, y₅′为较重污染天气按日期先后序列排序的源强估计值，计算结果如表 4所示.

综合表 3、4，岭北区2014年—2015年采暖期、2015年—2016年采暖期、2016年—2017年采暖期NO_x、SO₂的源强估算结果(10^-8t/(m²·d))分别为36.46和16.31、46.14和11.25、45.25和7.74.可以看出，2014年—2016年岭北区采暖期NO_x的排放源强逐年上升，可能与岭北区道路改造以及城市汽车保有量的增加有关^[17]；岭北区采暖期SO₂浓度逐年降低，可能与哈尔滨市2014年及之前的以“蒙煤”为主的褐煤逐年采购与燃烧量的增加，近些年出台政策监控煤质限购褐煤导致采暖期煤质转优从而二氧化硫排放量降低有关^[18].2014年—2015年及2015年—2016年源强结果中值误差e_0.5均小于30%, 说明虽然各个年份采暖期的排放源强在一定范围内波动，结合最重污染天气下的气象因素更吻合大气箱式模型的基本假设，此时在最重污染天气下估算得到的大气污染物排放源强结果可信度最高.将此最重污染天气下源强估算结果作为哈尔滨市采暖期大气污染物排放平均源强是可行的，且具备一定科学性.

2.2.2 估算模式的不确定性分析

本估算模式主要考虑数据的不确定性，包括数据的缺失和失真.以2016年—2017年采暖期岭北区最重污染天气时段11月4日15:00—20:00时均NO_x浓度和SO₂浓度为例，采用估算模式得到最优拟合结果后，分别使用坐标系中处于拟合曲线上方和下方的真实值做拟合曲线，探讨监测站点监测数据与大气实际污染源强关系模型的不确定性.各个区域NO_x和SO₂排放源强的不确定性区间如表 5所示.可以看出，绝大部分地区拟合最优条件下排放源强均处于不确定性分析区间内，占总不确定性区间个数的13.6%，是可以接受的.本模型的不确定性主要是监测数据的不确定性导致的，为提高模型对大气污染源强估算的精确度，对于各个分区监测站点应该加设大气监测点，且监测数据由时均值提高到30 min均值、15 min均值等措施，均能降低数据的不确定性对模型估参的影响.

2.3 哈尔滨市总源强计算结果

哈尔滨市区研究区域总面积为1 691.301 1 km², 通过泰森多边形法由表 2计算可得哈尔滨市区采暖期大气污染物NO_x、SO₂、CO总的估算排放源强为436.33，154.82和8 566.82 t/d.

3 结语

本文构建了哈尔滨市区大气污染物排放强度估算模式, 根据对源强算法的精确度验证中所算得中值误差可知，中值误差小于15%，说明基于大气箱式模型反算源强的算法具有很好的适用性，适合作为哈尔滨市大气污染物源强排放模拟与预测的模型理论基础，计算得到哈尔滨市区大气污染物NO_x、SO₂、CO总的估算排放源强为436.33，154.82和8 566.82 t/d.基于ArcGis的泰森多边形网格工具推广到大气污染研究区域的划分，能够避免人为主观划分区域范围的不确定性，使得大气排放源强模型更具有可靠性和说服力.