混合标气标校温室气体数据的计算方法与实现

2024-03-13蓝维嘉韦汉勇詹棠邱国强

广东气象 2024年1期

蓝维嘉，韦汉勇，詹棠，邱国强

（珠海市国家气候观象台，广东珠海 519000）

精确的温室气体观测数据能客观反映相关区域碳排放状况，是各类碳研究模式得以有效运行、各级政府及职能部门应对气候变化、实施节能减排提供决策的依据［1－2］。目前国内使用较多的Picarro系列温室气体在线分析仪（简称“分析仪”）是基于波长扫描光腔衰荡光谱技术（Wavelength Scan Cavity Ring Down Spectroscopy，WS-CRDS）的高精度观测系统［3－4］。但是分析仪在连续观测期间会发生微小的漂移，因此对于痕量级的CO2（单位量级parts per million by volume，ppmv）、CH4（单位量级parts per billion by volume，ppbv）等观测要素，原始数据的质量控制和标校订正必不可少。张东启等［5］通过保证观测系统稳定性和对标气进行逐级严格标定来实现CO2、CH4观测数据的质量控制；刘立新等［6］建立了科学规范的本底大气CO2采样观测过程的质量保证与质量控制方法；邓雪娇等［7］对广东省温室气体观测数据进行了对比分析，对温室气体观测质控方法进行研究，掌握了观测标准传递、标校方法流程，建立了温室气体可视化应用平台；臧昆鹏等［8］组装新型CO2和CH4混合标气标校系统，研究建立了简便高效的标校流程和方法。

基层监测站往往不具备国家级标校实验的条件，无法实时对设备进行专业的校准，混合标气标校法便是各站对数据进行后期标校订正的一种常用方法。对此，各站基本依托外部第三方开发的应用系统，但由于其源代码未公开，各个流程不可见、不可控，这不仅使应用系统维护不便，还制约着台站人员对数据的质量控制和应用。因此，开发出技术可控、运维方便的温室气体数据混合标气标校软件或工具，就显得十分必要。本研究应用JAVA语言编译了混合标气标校温室气体数据的计算机方法［9］，对各个计算流程进行详细介绍（图1）。

图1 温室气体数据混合标气、标校法流程设计示意图

1 混合标气、标校温室气体数据的流程

1.1 原始数据的质量控制及平均值计算

1）原始数据文件及关键字段。

样气（空气）经过二级过滤、二级除水、流量控制等处理，在分析仪内进行WS-CRDS分析后，生成名为CFADSaaaa-yyyyMMdd-HHmmssZ-Data-Log＿User＿Sync（aaaa为设备序列号）的原始数据文件（每5 s一行数据），关键字段如下：

（1）DATE、TIME：日期、时间（世界时）；

（2）ALARM＿STATUS：报警状态，0表示分析仪状态正常，否则为异常提示；

（3）MPVPosition：多口阀代码。以珠海为例，1为样气、4为高质量浓度标气WH、5为中质量浓度目标气T、6为低质量浓度标气WL；WH（WL）的标称值一般应高于（低于）站点样气质量浓度的最高（最低）水平。

（4）CO2＿dry＿sync、CH4＿dry＿sync：经水汽校正后，样气的CO2、CH4质量浓度实测值［10］。

2）气路切换与相关计算方法。

按照生成时间，遍历所有原始数据文件，逐行解析。首先需要判断本行的多口阀代码与上一行是否一致，条件判别语句为：

if（（mpv＿1－mpv＿0）！＝0）｛／／mpv＿1、mpv＿0：分别为本行、上一行数据的多口阀代码

多口阀代码不一致的，表明标气与样气的气路已通过多口阀切换。气路切换后2 min内的数据应剔除，因为气路中的气体是标气和样气的混合态（对已运行时间较长或者响应速率变慢的分析仪，以3 min为宜）；此后，气路中的气体质量浓度基本达到平衡（图2）。气路切换与否，是计算方法变换的重要依据。

图2 气路切换前后的CO2质量浓度变化曲线

本研究归纳了5种气路切换情形，其对应的计算步骤和计算机代码如下：

（1）情形1。气路不变，测标气。

①判断多口阀位置：mpv＿0＝＝mpv＿1、mpv＿0＞1、mpv＿1＞1。

②确保距上次气路切换时间已经超过2～3 min（120～180 s）：

（sec.ss（t＿str＿1，t＿mpv＿changed））＞＝180；其中，sec.ss（）是判断两个时间相差秒数的方法，t＿str＿1为本行数据时间，t＿mpv＿changed为最近一次气路切换时间。

③判断本行数据时间与上一行数据时间处在同一个标气观测时段（间隔小于10 s）。

④将本行标气值写入标气原始观测样本数组（每个数组包含一个标气观测时段的后2 min的数据，每5 s一个观测数据，共计24个数据）：

co2＿mpv［k＿5min－1］＝Double.valueOf（zz.substring（416，428））*100；

ch4＿mpv［k＿5min－1］＝Double.valueOf（zz.substring（468，480））*1000；

其中，co2＿mpv［］、ch4＿mpv［］为CO2、CH4的标气原始观测样本数组。

（2）情形2。气路不变，测样气，且要求本行与上一行数据的时间在同一个5 min时段——“5 min时段”，指的是正5 min时刻（含）之前的300 s。例如，00：00：01—00：05：00的60个样气有效数据的平均值，表示为“00：05：00的平均值”；若因气路切换、数据中断等原因造成5 min时段的样气有效数据不满60个，如00：00：01—00：04：35或00：02：30—00：05：00的样气有效数据的平均值，均表示为“00：05：00的平均值”。

①判断多口阀位置：mpv＿0＝＝1、mpv＿1＝＝1；

②气路切换时间已经超过3 min（180 s）；

③判断本行与上一行数据在同一个5 min时段：sec.ss（w5＿1，w5＿0）＝＝0；其中：w5＿1、w5＿0分别为本行、上一行数据所对应的正5 min时刻。

④5 s样气实测值累加：

co2＿5min＿all＝co2＿5min＿all＋Double.valueOf（zz.substring（416，428））*100；

ch4＿5min＿all＝ch4＿5min＿all＋Double.valueOf（zz.substring（468，480））*1000；

k＿5min＋＋。

（3）情形3。气路不变，测样气，但是本行与上一行数据的时间已经不在同一个5 min时段。

①判断多口阀位置：mpv＿0＝＝1、mpv＿1＝＝1；

②判断本行与上一行数据不在同一个5 min时段：sec.ss（w5＿1，w5＿0）！＝0；

③计算样气5 min平均值：

data＿co2＿ave＝Double.valueOf（String.format（″%.4f″，co2＿5min＿all／k＿5min））；

data＿ch4＿ave＝Double.valueOf（String.format（″%.4f″，ch4＿5min＿all／k＿5min））；

④将上一行数据对应的正5 min时刻作为样气5 min平均值的探测时间。

（4）情形4。气路从样气切换至标气。

①判断多口阀位置：mpv＿0＝＝1、mpv＿1＞1；

②计算样气5 min平均值；

③将上一行数据对应的正5 min时刻作为样气5 min平均值的探测时间。

（5）情形5。气路从标气切换至样气。

①判断多口阀位置：mpv＿0＞1、mpv＿1＝＝1；

②某标气测定时段内（通常为5 min时长，每5 s 1个、共计60个实测值），标气测定值序列为｛c1、c2、……c59、c60｝。取后2 min的观测序列，作为1.1节情形1中提出的标气原始观测样本数组co2＿mpv［］、ch4＿mpv［］，依次计算1 min时长的数据序列｛c49、c50…c59、c60｝、｛c48、c49…c58、c59｝、……｛c38、c39…c48、c49｝、｛c37、c38…c47、c48｝的标准差σ；

③对CO2（CH4），选取σ最小、且小于0.1 ppmv（或2 ppbv）的1 min数据序列，其平均值为该时段的标气实测平均值。关于标气采集样本的选取原则，不在于数量多，而在于质量精，可保证最终所选的标气样本集合的平衡态最佳［11］；

④若有效测定序列为｛ck－11、ck－10……ck－1、ck｝，则选取ck的测定时间作为标气实测平均值的对应时间。

3）计算程序与生成文件。

利用Data＿5min＿ave.java来实现样气5 min平均值和标气实测平均值的计算，并形成“STA＿5min＿ave.txt”文件（表1）。；STA为站名，下同。

表1 STA＿5min＿ave.txt1）文本格式示例

1.2 标校参数的获取（基于混合标气标校法）

1）标校参数计算。

设温室气体浓度的实测值为X、标称值为Y，对符合线性条件的设备［12］，X、Y符合最小二乘法原则，即X、Y满足线性回归方程Y＝aX＋b（a、b为标校参数）。以珠海为例，分析仪观测一轮样气、标气的时长分别为2 h和5 min，顺序依次为：WH、样气、T、样气、WL、样气、WH……。WH、WL用于计算标校参数；T作为目标气，用于检验标校参数。

举例假设分析仪的观测时间和要素如下：07：55—08：00观测WH、08：00—10：00观测样气、10：00—10：05观测T、10：05—12：05观测样气、12：05—12：10观测WL……。联立方程组：

其中，WH、WL为标气实测平均值；WH1、WL1为标称值；解得标校参数a、b，用于标校08：00—12：05的样气5 min平均值。

2）标校参数检验。

分析仪在连续观测期间会发生微小的漂移，因此，不同时段内的标校参数不同。1.2节例子中，令目标气计算误差k＝｜（aT＋b）－T1｜（T为目标气实测平均值，T1为目标气标称值）。如果对于CO2（CH4），k＜0.1 ppmv（或2 ppbv），则a、b值是有效的［10］；否则，应采用上一个样气观测时段的a、b值进行替换。

3）计算程序与生成文件。

利用Data＿calibration.java计算标校参数，并生成名为“STA＿a＿b.xlxs”的标校参数文件。

（1）标气实测平均值：从“STA＿5 min＿ave.txt”中读取的标气实测平均值；

（2）标气标称值：对应时间点的标气标称值；

（3）差值：标气实测平均值与标称值的差；

（4）目标气计算误差：将目标气实测平均值T、标称值T1代入公式｜（a×T＋b）－T1｜的结果；

（5）a＿CO2、b＿CO2、a＿CH4、b＿CH4：通过公式（1）和（2）计算所得、且经过检验或替换后的标校参数。

标校参数与时间点的对应关系如下：

（1）t1、t3为相邻的WH、WL标气观测时间点，根据t1、t3的标气实测平均值、标称值计算出时间t1至t3的标校参数a＿CO2、b＿CO2；之后再用t1、t3之间的目标气T观测时间点（即t2）的实测平均值、标称值来验证a＿CO2、b＿CO2。

（2）时间t3至t4没有目标气T观测，则选取t4之后、距离t4最近的一次目标气T观测时间点（即t5）的实测平均值、标称值来验证时间t3至t4的标校参数。

（3）t1对应行的a＿CO2、b＿CO2是时间t1至t3的标校参数，t3对应行的a＿CO2、b＿CO2是时间t3至t4的标校参数，依此类推。

4）算法的可靠性评估。

利用混合标气标校温室气体数据的计算工具，计算珠海市某温室气体监测站2014年全年的CO2、CH4目标气计算误差k。该监测站在2013年12月底正式运行，2014年业务运行期间满足以下条件：（1）新投入使用的分析仪性能稳定、漂移较小；（2）WH、T、WL等3瓶标气均由中国气象科学研究院温室气体及相关微量成分实验室标定，定值结果可溯源至世界气象组织CO2、CH4中心标校实验室维持的国际一级标气，期间标气质量浓度较稳定且不更换。2014年，该站的绝大部分CO2目标气计算误差k＜0.1 ppmv，所有CH4目标气计算误差k＜2 ppbv（表2），这说明编译出的混合标气标校温室气体数据的计算方法，其计算精度满足WMO的要求。

表2 珠海市某温室气体监测站2014年CO2、CH4目标气计算误差k值分布

1.3 样气5 min平均值的标校

利用Data＿output.java计算标校值，将“STA＿5min＿ave.txt”文件的样气5 min平均值，代入相应时段的标校公式Y＝aX＋b，可得出标校值，生成“STA＿output.txt”文件。

2 算法集成与运行

运行算法集成模块QC＿manual.java（图3），一键执行Data＿5 min＿ave.java、Data＿calibration.java、Data＿output.java，完成多年历史温室气体原始数据的质量控制和标校订正。该计算机工具使温室气体数据标校更方便快捷，根据标校法计算工具的测试情况，完成一个站点7年原始5 s数据（超过4 000万行）的质量控制、标校订正，并生成样气5 min平均值、标校参数、样气标校值等数据文件，时长约为8 min。

图3 算法集成模块QC＿manual.java（手动质控）界面

该工具还提供更多排查系统故障的途径，以标校参数文件“STA＿a＿b.xlxs”为例。当标气实测平均值与标气标称值差异过大时，提示标气质量浓度发生改变、标气更换记录未输入电脑等可能（图4），可为维护人员更换标气、完善标气记录提供参考依据。

图4 标气实测平均值与标称值差异过大（红色方框处）界面（提示可能未及时录入标气更换记录）

当标气实测平均值为0时，提示可能出现以下情况：（1）标气质量浓度稳定，但分析仪出现震荡性漂移，从而导致对标气的观测达不到规定标准；如果分析仪对所有标气的测量都无法达到平衡态，则应对分析仪进行人工校准甚至检查激光发生器、光腔等关键组件是否故障，以提高分析仪运行效率和观测数据精度。（2）分析仪不存在漂移，但是标气质量浓度出现震荡性变化；若分析仪只是对某种标气的测量无法达到平衡态，则需要对相应的标气进行重新标定或更换标气。