1 引言(Introduction)

近年来, 随着经济的快速增长和城市规模的不断扩大, 我国大气环境污染问题日趋加剧, 大中型城市重霾天气频发, 而细颗粒物PM_2.5是主要污染物之一(吴兑, 2012；吴兑等, 2014；Ao et al., 2016).目前, 大多数大气污染过程的实验观测和数值模拟研究主要集中在与污染过程有关的天气条件变化(蒋伊蓉等, 2015；Ye et al., 2016；Li et al., 2018；Wang et al., 2018；吴战平等, 2019)、边界层高度变化(王耀庭等, 2012；杜川利等, 2014；Luan et al., 2018；Miao et al., 2018)和气象要素变化(Zhang et al., 2015；吴序鹏等, 2018；Zhong et al., 2018；尹承美等, 2019)等方面, 少有针对重霾天气过程大气细颗粒物的湍流垂直输送和湍流通量的研究.而细颗粒物湍流垂直通量的准确估算是影响空气质量预报准确与否的重要参量之一, 有利于空气质量和气候数值模式中排放参数化方案的发展与改进.

由于目前缺乏较为成熟可靠且使用广泛的颗粒物快速响应测量仪器, 而利用传统的涡动相关法(Eddy Covariance, EC)直接计算颗粒物湍流通量存在一定的难度, 因此, 基于条件采样的松弛涡旋累积法(The Relaxed Eddy Accumulation, REA)逐渐成为一种实用、有效的获取气态物质或颗粒物湍流通量的方法.REA法利用采样时间区间内气态物质或颗粒物的平均浓度计算湍流通量, 其显著优点在于不需要气态物质或颗粒物的快速响应测量仪器, 计算湍流通量时不需要进行复杂的仪器学修正, 后期的数据处理过程较为简单(Ueyama et al., 2013), 可有效解决缺少气态物质或颗粒物快速涨落测量手段的情形或测量精度低、误差大的问题.因此, 在大气气溶胶和大气成分分析领域, REA法比EC法适用范围更为广泛, 目前经常用于CH₄(葛红星等, 2016；Ueyama et al., 2012；2013；2018)、NO₂(Desjardins et al., 2010)、氨(Nelson et al., 2017；2019)、挥发性有机物(VOC)(Arnts et al., 2013)等气态物质和气溶胶(Held et al., 2008)湍流通量的获取.

由于目前缺少EC法所需的颗粒物快速响应测量仪器, 因而难以获取细颗粒物的湍流通量.因此, 本文采用REA法进行模拟条件采样, 利用山东省平原县大气环境实验站观测资料, 计算2018年12月27日—2019年1月8日持续重霾天气过程的细颗粒物湍流通量, 并与细颗粒物湍流通量测量装置给出的结果进行对比, 以验证REA法计算细颗粒物湍流通量的适用性.

2 数据获取与处理(Data acquisition and processing)

平原县大气环境实验站位于山东省德州市平原县气象局(37.154144°N, 116.476774°E, 见图 1), 地处中国华北平原, 德州市南部, 观测场海拔高度为23.8 m, 周围是相对平坦的农田, 无高大建筑.本文选取大气环境加强实验期间, 即2018年12月27日—2019年1月8日13 d重霾天气过程的大气湍流和颗粒物实验观测数据.与本文分析有关的主要观测项目包括:涡动相关系统(IRGASON, Campbell Scientific, Inc., USA)获取风速、温度、湿度和CO₂快速涨落信息, 超声风温仪朝向正北方向, 采样频率为10 Hz；颗粒物连续测量仪(E-Sampler, Met One, Inc., USA)获取细颗粒物质量浓度的高频信息, 采样频率为1 Hz.细颗粒物湍流通量测量系统(一种细颗粒物湍流通量测量方法, 专利受理号:201811376291.9)是一种利用大气消光特征获取颗粒物湍流通量的测量装置(Ren et al., 2020).仪器安装高度均为2.8 m.

大气湍流观测数据采用EddyPro软件(Advanced 6.2.1, LI-COR Biosciences, Inc., USA)进行数据质量控制和预处理.具体过程包括:野点剔除(Vickers et al., 1997)、二次坐标旋转(Wilczak et al., 2001)、仪器学订正等, 同时剔除了明显不合理的数值, 数据平均时间区间为30 min.后续研究主要使用EddyPro软件Advanced模式中高级设置后输出的经过以上修正的10 Hz高频Level 6数据, 并对垂直风速数据分段每10 s取一个平均得到1 Hz的数据, 用于模拟条件采样计算细颗粒物湍流通量.

细颗粒物质量浓度脉动数据的预处理仅包括野点剔除, 平均时间为30 min, 通过计算其平均值和标准差对野点进行判别, 即:若脉动数据中某一数值相对其平均值的偏差超过4倍的标准差, 即: C-C ≥4σ_C, 则认为该数据点为一个“野点”, 将“野点”用其前后相邻两数据的线性内插代替.

3 理论和计算方法(Theory and calculation method)

松弛涡旋累积法是一种基于条件采样的气态物质或颗粒物湍流通量获取方法(Businger et al., 1990), 是在涡旋累积法(Eddy Accumulation, EA)基础上, 放宽采样条件的简化方法(Desjardins, 1977).REA法不需要EA法中进行条件采样所必需的复杂阀门装置, 仅需用恒定速率抽气即可.一般来说, 湍流通量的计算方法都存在一定的前提假设, EC法的相关假设也适用于EA法和REA法(Bowling et al., 1998).首先, 大气中的标量气体应具有无化学反应的或惰性的属性, 如甲烷、二氧化碳等(Businger et al., 1990)；其次, 假设某一实验站垂直输送的全部通量代表实验站所在区域的湍流通量的总体平均值(Businger, 1986)；最后, 大气运动比较平稳, 流场水平均匀(忽略平流的影响), 准平稳(定常)湍流, 近地层满足常通量层假设, 影响湍流通量的各种尺度涡旋都能被测量到(徐自为等, 2008).REA法基于垂直风速进行条件采样, 当垂直风速向上(w>0)时, 将气态物质收集至向上的采样罐；当垂直风速向下(w < 0)时, 将气态物质收集至向下的采样罐, 从而获取一定时间区间内垂直风速向上和向下两个采样罐内待测气态物质或颗粒物的平均浓度和, 则待测气态物质或颗粒物的湍流通量可以表示为(Businger et al., 1990):

(1)

式中, b为REA法的经验系数, σ_w为垂直风速标准差.但条件采样获取的气态物质或颗粒物的平均浓度差通常很小, 甚至可能小于气态物质或颗粒物测量仪的精度, 使得计算得到的湍流通量误差较大.为此, 有研究发现在条件采样过程中引入一定的采样阈值w_d可有效增大获取的平均浓度差, 提高计算精度(Bowling et al., 1998；Myles et al., 2007).即:当|w|≤w_d时, 认为颗粒物浓度C=0.通常引入归一化采样阈值H_REA, 以使不同平均时间的结果具有可比性, 即:

(2)

在没有进行实际条件采样的情况下, 也可以通过EC系统观测获取的垂直风速w和气态物质或颗粒物质量浓度的高频观测数据进行模拟条件采样来获取式(1)中的浓度差(Andreas et al., 1998；Ammann et al., 2002；Sakabe et al., 2012；Tsai et al., 2012), 同时选取一定的采样阈值w_d.在每个采样平均时间区间内, 分别对垂直风速w>w_d和w < -w_d对应的浓度C的高频数据取平均, 得到该时间区间内的平均浓度和, 进一步得到平均浓度差.

基于近地面层各物理量湍流输送特征的相似性, 利用REA法计算大气气态物质或颗粒物的湍流通量时, 通常以较易获取快速涨落信息的温度(T)、水汽(H₂O)和二氧化碳(CO₂)等作为代理变量, 确定经验系数b值.前人的许多实验研究发现, 以EC法获取的湍流通量F_EC为参考, REA法的经验系数b可通过式(3)计算.即, 对待测气态物质或颗粒物及某种代理变量(温度(T)、水汽(H₂O)和二氧化碳(CO₂))进行高频采样, 在每个采样平均时间区间内, 将EC法计算的代理变量的湍流通量、垂直风速标准差及模拟条件采样获取的代理变量的平均浓度差()代入公式(3), 得到REA法的经验系数b, 然后对代理变量的b值时间序列直接取中位数, 将此b值及模拟条件采样获取的待测气态物质或颗粒物的平均浓度差()代入公式(1)即可得到其湍流通量.

(3)

4 结果与讨论(Results and discussion) 4.1 代理变量法获取经验系数b值

REA法中引入一定的采样阈值可增大平均浓度差, 提高计算精度, 但过大或过小的采样阈值都会给湍流通量的计算带来很大的误差；同时, 要获取准确的待测气态物质或颗粒物湍流通量, 还需要调整经验系数b值进行补偿.因此, 本节分析了温度(T)、水汽(H₂O)和二氧化碳(CO₂)3个代理变量的b值与归一化采样阈值H_REA的关系, 从b值离散性角度确定理想归一化采样阈值所对应的代理变量的经验系数b值.假定H_REA的测试范围为0~2, 步长为0.1, 共21组数据.对每一个代理变量, 利用不同H_REA所对应的采样阈值w_d分别进行模拟条件采样, 分别得到21组b值的时间序列.

图 2分别给出了温度(T)、水汽(H₂O)和二氧化碳(CO₂)不同代理变量的b值b_T、b_H2O、b_CO2随归一化采样阈值H_REA变化的箱线图.其中, 箱体中间的实线代表H_REA对应的b值时间序列的中位数, 箱体的上、下两端分别代表第三四分位点Q₃(75%)和第一四分位点Q₁(25%), 箱体的高度为Q₃与Q₁之间的差值, 代表四分位距(Interquartile Range, IQR=Q₃-Q₁)；箱体以外两端虚线上、下延伸至的横实线为数据的上、下限, 分别定义为:上限MAX=Q₃+1.5IQR, 下限MIN=Q₁-1.5IQR.上、下限以外的数据点称为离群点, 未在图 2中体现.IQR和上、下限范围可在一定程度上体现b值的离散性.

由图 2可以看出, 不同代理变量T、H₂O和CO₂的b值的中位数均随H_REA的增大而减小, 箱线图的箱体高度及箱体上、下限范围均表现出先减小后增大的趋势, 说明b值的离散性均随H_REA的增大表现为先减小后增大.H_REA的数值范围在0.6~0.9之间对应的代理变量b值的离散性较小.当H_REA=0, 即采样阈值为0时, 不同代理变量b值的中位数分别为b_T=0.56, b_H2O=0.56, b_CO2=0.56.

对于REA法来说, 经验系数b的取值对气态物质或颗粒物湍流通量的获取十分重要.前人的理论研究发现, 当垂直风速与待测气体浓度之间的联合概率密度分布(JFD)满足理想的二维正态分布关系时, 经验系数(Baker et al., 1992；Wyngaard et al., 1992；Katul et al., 1996).同时, 前人也开展了大量的实验研究, 但目前的实验研究和应用更多的局限在农田、森林和草地下垫面.一些不同下垫面条件下的经验系数b参考值见表 1.可以看到, 经验系数b的取值范围多在0.5~0.6之间, 不同下垫面差别较小.而本文此次实验的下垫面属于农田下垫面, 取值与前人结果基本一致.未来若能对更多下垫面类型展开实验研究, 并确定经验系数b参考值将对扩大REA法的适用范围有重要意义.

图 3给出了不同代理变量b值的IQR随H_REA的变化关系.从图中可以看到, 不同代理变量b值的IQR均随H_REA的变化存在一个极小值, H_REA=0.8时的代理变量b值的离散性最小, 这与Tsai等(2012)和葛红星(2018)的结果一致.由图 2可见, H_REA=0.8时对应的温度(T)、水汽(H₂O)和二氧化碳(CO₂)的经验系数b值的中位数分别为b_T=0.3491、b_H2O=0.3481、b_CO2=0.3516.另外, 温度(T)的IQR最小, 说明温度对应的b值的离散性最小, 水汽(H₂O)次之, 二氧化碳(CO₂)的b值的离散性最大.所以, 温度(T)更适合作为代理变量计算细颗粒物的湍流通量.另外, 温度还具备显著优势:三维超声风温仪比气体分析仪更为成熟, 温度快速涨落信息较易获取, 观测精度较高(Bowling et al., 1998).

4.2 细颗粒物湍流通量的获取与验证

利用公式(2)可计算每个采样时间区间内的采样阈值w_d, 本文取其均值0.1 m · s^-1作为最终的采样阈值, 通过模拟条件采样获取细颗粒物平均浓度差.计算细颗粒物湍流通量时, 经验系数b的使用通常有3种方法:①“同步b值”(Myles et al., 2007), 即:由于不同时间区间内的b值不同, 利用每个采样平均时间区间内的b值计算该时间区间内的湍流通量.Oncley等(1993)使用“同步b值”计算了棉花田下垫面的二氧化碳通量；Gao(1995)使用“同步b值”计算了森林下垫面的动量通量、感热通量和潜热通量；Myles等(2007)使用“同步b值”计算了草地下垫面氨、二氧化硫等气体的湍流通量.但由于经验系数b值本身随时间变化不大, 用“同步b值”反而会引入b值的离散性带来的误差, 不建议采用.②“固定b值”(Sakabe et al., 2012；Ueyama et al., 2012), 即:对代理变量的b值时间序列取中位数计算所有时间区间内的湍流通量, 剔除b值之间的差异造成的离散性所带来的误差.Baker等(1992)使用“固定b值”计算了大豆田下垫面的感热通量；Sakabe等(2012)和Ueyama等(2012)都使用“固定b值”计算了森林下垫面的甲烷通量.③“修正b值”(Ammann et al., 2002), 即:考虑大气层结状况对经验系数b值的影响, 用稳定度参数z/L对固定b值进行修正.Ammann和Meixner(2002)研究了谷物田下垫面温度、水汽和二氧化碳的经验系数b值与稳定度参数之间的关系, 提出了“修正b值”, 并发现计算通量结果优于“同步b值”.但由于使用经验公式给出的“修正b值”本质上是对“固定b值”进行稳定度补偿, 其离散性明显小于“同步b值”, 因此, 该对比研究并不能确定通量获取结果的改善是加入稳定度修正还是离散性减小的作用, 使用“修正b值”的效果是否优于“固定b值”无法确定.可见利用经验系数b值的何种方式计算湍流通量似乎与下垫面没有太大联系.葛红星(2018)使用“固定b值”和“修正b值”计算了水稻田下垫面的甲烷通量, 发现两种方法的结果差异很小, 且加入采样阈值后不同稳定度条件下b值的差异将缩小, “固定b值”与“修正b值”获取的通量结果差异将更小.因此, 本文也沿用了其结论:“固定b值”的获取过程比较简单, 在大气层结状况变化不大的情况下是更好的选择.

图 4为山东省平原县大气环境实验站2018年12月27日—2019年1月8日重霾天气过程的细颗粒物质量浓度、垂直风速及REA法计算的细颗粒物平均质量浓度差的时间序列.图 4b中两条白色虚线分别表示采样阈值w_d=0.1 m · s^-1和-w_d=-0.1 m · s^-1.可以看到, 当细颗粒物质量浓度较大时, 即重污染时期, 垂直风速脉动较弱, 细颗粒物平均质量浓度差越大；反之, 垂直风速脉动较强, 细颗粒物平均质量浓度差越小, 这与颗粒物的干沉降有关, 而造成干沉降的原因之一即是湍流扩散.垂直风速绝对值小于阈值w_d的时段, 忽略其对大气细颗粒物湍流通量的贡献, 取其平均浓度差和湍流通量为0.

图 5给出了山东省平原县大气环境实验站2018年12月27日—2019年1月8日重霾天气过程中的细颗粒物湍流通量时间序列.其中, 黑色实线为REA法的计算结果, 红色点线为利用大气消光系数特征的细颗粒物湍流通量测量装置的计算结果.由图 5可以看出, 整体上两种方法获取的细颗粒物湍流通量有较好的一致性.由于湿度的影响, REA法计算得到的细颗粒物湍流通量比细颗粒物湍流通量测量装置的结果略有偏大, 且在污染较为严重的1月4日和1月7日时段的偏差略有增加, 说明REA法在获取大气细颗粒物湍流通量上是可行的.

图 6给出了REA法和细颗粒物湍流通量测量装置获取的细颗粒物湍流通量的对比, 其中, 虚线斜率为1:1, 实线为两种方法计算结果的线性拟合, 拟合系数为0.93, R²=0.7987, 较分散数据组对应了细颗粒物湍流通量较大的重霾天气污染时段.

以上对比说明两种方法在计算颗粒物湍流通量上的结果是比较一致的, 但在重污染时段依旧存在微小差别.大气消光法和REA法计算垂直湍流通量各有优势, 也有不足.大气消光系数法的优点在于不需要颗粒物浓度的快速响应测量仪器, 缺点是当相对湿度较高时颗粒物质量浓度并不是影响消光系数的唯一因素, 两者之间的函数关系需要一定修正, 因此, 大气消光法往往用于相对湿度较低(RH < 80%)的情形；另外, 大气消光法对能见度测量的精度要求较高.松弛涡旋累积法的优点同样不需要颗粒物浓度的快速响应测量仪器, 缺点是计算过程中的气态物质平均浓度差为两个较大数值之差, 是颗粒物湍流通量相对误差增加的原因之一.

5 结论(Conclusions)

本文利用山东省平原县大气环境实验站加强观测期间的大气湍流和细颗粒物质量浓度脉动观测资料模拟条件采样, 采用REA法计算了2018年12月27日—2019年1月8日重霾天气过程的细颗粒物湍流通量.结果显示, 加强实验观测期间, 温度(T)、水汽(H₂O)、二氧化碳(CO₂)等不同代理变量的REA法获取的湍流通量的经验系数b分别为b_T=0.56、b_H2O=0.56、b_CO2=0.56；经验系数b值的中位数随归一化采样阈值H_REA的增大而减小；考虑四分位距IQR代表b值的离散性, 不同代理变量的IQR数值随H_REA的变化均存在一个最小值, 据此确定了b值离散性最小时对应的理想归一化采样阈值H_REA=0.8和各自的b值b_T=0.3491、b_H2O=0.3481、b_CO2=0.3516.基于近地面层不同物理量湍流输送特征的相似性, 将较易获取快速涨落信息的温度(T)、水汽(H₂O)、二氧化碳(CO₂)等作为代理变量确定经验系数b值, 选取b值离散性最小的温度(T)的经验系数b值和平均采样阈值w_d=0.1 m · s^-1, 采用REA法进行模拟条件采样, 计算了细颗粒物的湍流通量, 并与利用大气消光系数特征的细颗粒物湍流通量测量装置的结果进行对比.结果显示, 尽管在重霾天气时段, REA法计算的细颗粒物湍流通量与细颗粒物湍流通量测量装置结果的差异略有偏大, 但整体表现出了较好的一致性, 说明了利用REA法计算细颗粒物湍流通量方法的可行性.