六氟化硫(SF₆)具有良好的化学稳定性和绝缘性，自1980年代以来被广泛应用于电力行业，主要作为电气绝缘介质和灭弧剂，因而排放量大大增加.研究发现，SF₆具有长达3200年的大气寿命和极强的增温潜势(百年尺度的GWP值是CO₂的23900倍)(IPCC，2007)，因此，被列为《京都议定书》限排的6种/类温室气体之一.近年来，国内外学者也对我国SF₆等卤代温室气体源排放、变化趋势、辐射强迫及全球增温潜能进行了大量的研究(石广玉，1992；张仁健等，2000； 王明星等，2000；2002；IPCC，2007；吴金秀等，2009；张华等，2013).

采用传统的“自下而上”清单方法估算SF₆排放量时，在所需排放因子、活动水平数据的准确性和时效性等方面尚有较大局限.因此，国际上一些研究机构采用“自上而下”的方法，利用在线观测浓度数据和气象资料，结合各种模型估算区域及全球的排放源分布及源强，其中，较具代表性的有：利用1978—1996年大气CCl₄浓度观测资料研究了全球CCl₄变化趋势和排放状况；Mulqiey等(1998)采用澳大利亚国立大学开发的化学传输模式初步反演估算了CFCl₃通量；Ryall等(1998；2001)采用拉格朗日粒子扩散模式探讨了Mace Head地区的CFC-11传输，并结合观测资料初步估计了欧洲的排放量和源分布；Manning等(2003)利用Mace Head站6年观测数据及NAME模式，反演推算了欧洲地区CFC-11、CFC-12及CH₄、N₂O排放源的分布；Stohl等(2009；2010)通过改进反演方法，利用FLEXPART模式结合多站观测资料，反演了区域和全球HFC-134a、HFC-152a和HCFC-22排放量，以及东亚地区氢氟氯烃(HCFCs)和氢氟烃(HFCs)的排放量；Vollmer(2009)、Kim(2010)、Li等(2011)采用改进的反演方法结合FLEXPART模式，尝试反演了中国及东亚地区部分卤代温室气体排放量.但Vollmer等(2009)的研究方法还不是很完善，模拟网格比较粗，有些物种(如CFC-11、HCFC-22)明显偏高，而Kim等(2010)和Li等(2011)的研究则是基于韩国Jeju站点的观测资料来反演，对于中国区域的代表性差一些.目前，国内利用FLEXPART模式反演卤代温室气体排放源及评估验证的工作才刚刚起步(李岩等，2010；An et al., 2012).

因此，本文拟采用拉格朗日粒子扩散模式FLEXPART，模拟得到排放源对测站浓度的敏感性系数.同时，结合地面测站观测资料，建立观测浓度值与源排放的线性关系，并利用最小二乘法或其他方法优化该线性关系求解源排放，通过调整合理的参数，使得反向模式在利用调整后的源排放进行模拟时，模拟结果与观测数据匹配最好，这样所得到的调整优化后的源即为反演源.最后，利用2009年北京上甸子区域大气本底站SF₆浓度观测资料，结合FLEXPART模式尝试对中国区域SF₆排放量进行初步反演.

FLEXPART模式是由挪威大气研究所(Norwegian Institute for Air Research，NILU)开发的一种拉格朗日粒子扩散模式(Stohl et al., 2005).通过计算点、线、面或体积源释放的大量粒子的轨迹，来描述示踪物在大气中长距离、中尺度的传输、扩散、干湿沉降和辐射衰减等过程.该模式可通过时间的前向运算来模拟示踪物由源区向周围的扩散，也可通过后向运算来确定对固定站点有影响的潜在源区的分布，尤其当研究区域内观测站数量少于排放源数量时，后向运算更具有优势(Seibert et al., 2004).

FLEXPART模式最初开发应用于研究长距离、中尺度的点源空气污染扩散问题(Stohl，1996).Stohl等(1998)根据3个大型示踪实验CAPTEX(Cross-Appalachian Tracer Experiment)、ANATEX(Across North America Tracer Experiment)、ETEX(European Tracer Experiment)的数据对FLEXPART模式早期版本(Version 2.0)的模拟效果进行了评估，结果显示，FLEXPART 是当前最好的扩散模式之一.后来，该模式经过发展，成为了模拟和分析大气传输的一个广泛性工具，应用也从空气污染扩展到了大气传输领域(Stohl et al., 1999; 2003; Cristonfanelli et al., 2003; Avey et al., 2007).近年来，该模式逐渐被用于估算全球及区域卤代温室气体排放量及排放源的研究中.

FLEXPART模式的核心内容是研究大气污染物的源-受体关系(Source-receptor relationship)，简称SRR系数：污染排放为“源”，观测站点为受体.通过研究污染物水平输送、扩散、对流、干、湿沉降、辐射衰减和一阶化学反应等过程，可以得到随时间序列变化的格点污染浓度(正向模拟)或格点驻留时间(也称敏感性系数或印痕函数(后向模拟))(蔡旭晖，2008).

通过源-受体关系转换，计算格点驻留时间τ_i的公式可以表示为：

式中，T为时间分辨率(s)；N为T时间范围内采样的数量(个)；J为释放的粒子总数(个)；f_ijn是一个函数，它决定了对于指定格点有“贡献”的粒子的多少.

本文中选取的模拟时间为2009年2月10—12月31日，时间分辨率为3 h，格点分辨率为0.125°×0.125°，垂直方向分为100、300、3000、20000 m 4层.为了更好地与近地面观测资料匹配，选取模式第一层输出结果.模式输出结果为驻留时间(Residence times，也称为敏感性系数)，单位为s，表示单位质量污染气体在水平格点上的停留时间.

上甸子站(东经117°07′，北纬40°39′，海拔293.9 m)位于北京城区东北方向，坐落在密云县高岭镇上甸子村，距北京城区约120 km，是全国最早的3个区域大气本底站之一(图 1)，代表了我国华北经济发达区域和京津冀经济圈大气本底状况.该站地处暖温带半湿润气候区，夏季炎热多雨，冬季寒冷干燥，春、秋季短促，具有典型的华北区域气候特征.

图 1(Fig. 1)

图 1 北京上甸子区域大气本底站 Fig. 1 The location of Shangdianzi regional atmospheric background station

采用NCEP/NCAR每日4次(0：00、6：00、12：00、18：00 UTC)的GFS(Global Forecast System)模式分析资料，垂直分26 层，水平分辨率1°×1°.时间范围为2009年1月20—12月31日，在NCEP网站通过FTP实时下载(ftp://ftpprd.ncep.noaa.gov/pub/data/nccf/com/gfs/prod).

上甸子站大气SF₆浓度高精度在线观测资料，时间范围为2009年2月10—12月31日，原始观测数据为每分钟1个观测值，为方便与模式结果进行比较，处理成每3小时一次的平均浓度资料，有效样本数为1893.

反演排放量的核心是求解线性超定方程组.超定方程组即控制方程的个数大于未知量的个数，方程的解不唯一.在数学上，一般采用最小二乘法，通过逐步迭代使可得解与真实解的平均标准偏差达到最小，得到方程最优解.结合FLEXPART模拟结果与上甸子测站SF₆在线观测资料，求解该超定方程组的最优解.通过设置合理的参数，使方程及方程的解具有实际的物理意义.方程组表达式为：

式中，y表示模拟得到的浓度值，x是未知量，需要反演求得的排放源，M是m×n矩阵.用于实际的反演问题时，公式(2)可以写成：

式中，m为观测资料的时间序列；n为排放网格数； M _mn矩阵为不同格点排放对上甸子测站SF₆浓度产生贡献的系数矩阵，由FLEXPART模式结果按时间序列排列而成； Y _m为m行矩阵，由上甸子测站SF₆观测浓度-时间序列排列而成，需要注意的是，这里的观测值并非台站数据的直接应用，而是剔除了本底浓度(也称为基线浓度)后的观测值(即非本底浓度)，因为本底浓度代表滤除了局地条件和人为活动直接影响、已混合均匀的大气组成特征，非本底浓度常用于排放源反演(Reimann et al., 2008)； X _n为n列矩阵，数值代表了研究区域内SF₆的格点排放量，是未知量，也是本文需要求的方程的解.

解方程过程中需要输入初始参考值，由此引入了先验源的信息，公式(2)可变换为：

式中，y_o表示测站浓度观测值；x_a为先验源，格点排列与 X _n完全一致，数值可由人口数据和初步估计的排放总量计算得出，方程(4)可缩写为：

式中，.

参考最小二乘法，建立目标函数(6)，通过求微分使目标函数最小化，具体见式(7).

式中，σ₀代表来自观测的误差，下标0代表观测；σ_x代表来自先验源的误差.由此求解出的，

加上先验值x_a即为SF₆的格点排放量数据.反演流程如图 2所示.

图 2(Fig. 2)

图 2 FLEXPART模式反演SF6排放量流程图 Fig. 2 The flow chart of inversion way based on FLEXPART

大气SF₆浓度在线观测系统位于上甸子站实验室内，采样塔(高度10 m)位于实验室正北约60 m 处，采样口位于采样塔8 m高处.系统由样品采集、样品分离、样品检测3部分组成，具体组成、观测原理、标定和浓度处理方法等见文献(姚波等，2010).SF₆分析精度约为1%，并可溯源至Scripps Institute of Oceanography(SIO)—2005标准尺度.为与模式输出结果时次对应，选择同一时间段内(2009年2月10—2009年12月31日)的观测资料作为研究对象.

运用REBS(Robust extraction of background signal)算法(Ruckstuhl et al., 2012)，将上甸子测站SF₆浓度-时间序列筛分出本底浓度(也称为基线浓度)和非本底浓度(图 3).图 3b中，白色柱状图代表所有观测数据的概率密度分布，灰色代表本底浓度的概率密度函数分布，图 3右上方显示了估算的SF₆本底浓度的不确定性，约有40%的观测样本属于本底观测数据.

图 3(Fig. 3)

图 3 2009年上甸子大气本底站SF6浓度变化图(黑色圆点代表观测数据，实线为基线浓度值) Fig. 3 Baseline estimate for SF6 at Shangdianzi

反演方程(2)中，M_mn矩阵是FLEXPART模式结果按时间序列排列成的敏感性系数矩阵.原则上选取的反演区域应该包括所有时次FLEXPART印痕结果覆盖到的区域，这样才能避免遗漏重要的排放源对测站SF₆浓度的影响.

图 4为2009年2月10日—12月31日期间FLEXPART模式模拟的平均印痕结果.为了求解反演方程的方便(提高观测数据样本数与反演网格数比例)和提高计算效率，将反演区域内敏感性系数(SRR)小于一定数值的网格合并，图 4b显示的是合并后的不规则网格分布.从图 4可以看出，印痕主要分布在我国的华北区域，敏感性数值以上甸子为中心向周围递减.敏感性数值的大小代表了潜在源区对上甸子测站影响的程度.

图 4(Fig. 4)

图 4 2009年FLEXPART模拟平均印痕结果(a.规则网格,b.不规则网格; 国审字(2014)第612号) Fig. 4 Total surface footprint for Shangdianzi on a regular grid of 0.125°×0.125° resolution (a) and total surface residence time on irregular grid used in inversion (b)

先验源为初步估计的SF₆格点排放量，作为解方程的初始值，先验源设定越接近真值，方程迭代步数越少，程序运行越快，得到的解也越合理，因此，如何设定合理的先验源在反演过程中至关重要.

本文中先验源的计算用到两种数据，一为研究区域内预估排放总量，二为研究区域内格点分配系数.中国地区SF₆的预估排放量为参考EDGAR使用bottom-up方法得到的估计值(http://edgar.jrc.ec.europa.eu)，约为1.23×10³ t · a^-1.由于SF₆排放量与人为活动密切相关，主要广泛应用于铝镁冶炼及在电力行业的气体绝缘体和高压转换器上，可近似把我国人口密度作为格点分配系数来分配研究区域内的预估排放总量.图 5 给出的是2009年我国的人口密度分布，数据来自于CIESIN中心(the Center for International Earth Science Information Network，http://sedac.ciesin.columbia.edu/gpw)，分辨率为0.5°×0.5°，单位为人 · km^-2.从图中可看到，我国华北、河南、四川地区及长三角、珠三角等沿海地区人口分布密集，西部及内蒙地区人口分布较少.

图 5(Fig. 5)

图 5 研究区域人口密度分布(国审字(2014)第612号) Fig. 5 Population density distribution of the study area

图 6给出的是根据人口密度和不规则网格分派的SF₆先验源不规则网格分布.首先根据人口密度作为格点分配系数，摊派到均匀分布的网格上，然后再分配到根据敏感性系数大小合并的不规则网格上.从图中可以看出，SF₆先验源较大值分布在我国华北、四川及长三角、珠三角等地区，人口密集地区排放量大是符合实际情况的.

图 6(Fig. 6)

图 6 研究区域 SF6 先验源分布图(黑点为上甸子测站, 国审字(2014)第612号) Fig. 6 Maps of the irregular grid emissions of a-priori of SF6 in research area

利用上述观测浓度和先验源等信息，得到求解线性超定方程组的基础数据.利用3.1节中的反演方法，通过逐步迭代，求方程组的解，得到中国区域SF₆排放量的初步反演结果.图 7为反演源分布图及反演源与先验源的差值分布图.从图中可以看出，反演源的较高值分布在华北、四川、长三角及珠三角地区.对比先验源，发现北京、天津、河北地区的排放量有所增强.反演结果基本可以看出高SF₆排放主要分布在人口较密集、工业化程度较高的区域.

图 7(Fig. 7)

图 7 研究区域SF6反演源分布图(a)及反演源与先验源差值分布图(b)(黑点为上甸子站, 国审字(2014)第612号) Fig. 7 Maps of the irregular grid emissions of a-posteriori of SF6(a) and changes (a-priori minus a-posteriori) of SF6(b) as estimated from FLEXPART simulations

统计反演前后中国区域SF₆排放量发现，反演后SF₆排放量从1.23×10³ t · a^-1增加到1.25×10³ t · a^-1，排放量的不确定性从1.05×10³ t · a^-1 减小到0.72×10³ t · a^-1.本文反演的中国区域SF₆排放量1.25×10³(0.53×10³~1.97×10³)t · a^-1，其他学者反演的中国区域SF₆排放量分别为0.80×10³(0.53×10³~1.10×10³)t · a^-1(Vollmer et al., 2009)、1.3×10³(0.93×10³~1.7×10³)t · a^-1(Kim et al., 2010)和1.2 ×10³(0.9×10³~1.7×10³)t · a^-1(Li et al., 2011)，本文结果与之接近.将反演源、先验源模拟得到的SF₆浓度与观测的SF₆浓度比较，发现反演源的模式模拟结果比先验源的结果有了一定的改善，相关系数从0.37 提高到0.43，归一化均方误差(NMSE)从0.056减小到0.053，均方根减小了 2.64%.

本文利用FLEXPART模式结合上甸子测站观测资料，初步反演了中国区域SF₆排放量，结果发现，观测到的2009年上甸子区域本地站SF₆的浓度数据，约有40%的观测样本属于本底浓度.基于人口密度分配排放总量得到先验源的分布，体现了人口分布的特性，在华北、四川、长三角及珠三角等人口密集地区，排放量相对较大.利用反演方法得到的2009年中国区域SF₆的排放量为(1.25±0.72)×10³ t · a^-1，与他人研究结果相当.模拟结果与观测结果比较表明：反演源改善了模式模拟结果，相关系数从0.37 提高到0.43，归一化均方误差(NMSE)从0.056减小到0.053，均方根误差减小了 2.64%.

总体来说，运用模式结合观测反演排放源的原理较为成熟，国外已有不少相关研究，但国内对于该反演方法的研究尚处于起步阶段.本文的模拟结果基本反映出了中国地区SF₆排放量和排放源分布.但本文的研究尚存在一些不足及不确定性，文中先验源的设定是将排放量按人口密度数据分配到模式格点上的，SF₆排放与人为活动密切相关，而本文只考虑了人口数量这一因素，这可能会忽略其他因素(如大型电厂或大型镁铝冶炼厂的点源排放等)对排放量的影响，从而使得先验源的估计存在一些不可知的偏差.

本文反演中用到的敏感性系数矩阵来自于单站(上甸子站)的FLEXPART模拟结果，印痕覆盖区域有限，使得部分地区排放源不能被考虑到反演方程的实际计算中，在数学上，这不利于方程的解.国外一般采用多站点、长时间序列浓度资料做反演研究.因此，需要积累不同站点及更长时间的观测资料，这将有利于得到更符合实际意义的方程解.

总之，运用模式结合观测反演SF₆排放量的理论是合理的，但应用于国内排放研究时，具体算法及一些参数设置还需要进一步调整，需要积累更多站点更长时间序列的观测资料使得反演结果更完善及接近实际情况.