0 引言

甲烷(CH₄)作为第二重要的温室气体，由于其相对二氧化碳(CO₂)具有更高的增温潜势以及更短的大气寿命(约10 a)，对其排放进行控制被认为是短期内减缓气候变暖的有效手段^[1-3]。然而，全球大气CH₄背景值浓度从1750年的722 × 10^-9增长至2023年的1 920 × 10^-9[4]，升高了近2倍，且过去数年更是呈现加速增长的态势，2014—2017年，每年增长7×10^-9 ~ 12.7 × 10^-9[5]，2021年增幅达到了有记录以来的最高值(17.56 × 10^-9)^[4]。为了缓解这一趋势并科学地对CH₄排放进行管控，首先应对各类来源的CH₄排放进行有效的监测与估算。CH₄来源分为自然源与人为源，现今与人类活动密切相关的如化石燃料的开采和应用、农业面积的增加、反刍动物数量的增大、垃圾填埋量的增长等，是大气CH₄浓度增加的主要原因，人为排放已占全球排放总量的60%^[6]。对于这些人为排放源的排放估算，“自下而上”的排放清单与“自上而下”的量化估算仍存在较大差异，对于量化估算的方法，目前针对从局域到区域乃至更大范围尺度，已经有多种独立的监测估算技术，包括卫星^[7-8]、机载观测^[9-10]、高塔及站点等形式的地基固定监测^[11-13]、地面移动平台^[14]等。然而各类平台与技术，均存在一定的局限性与不确定性，在地面和100 m高度之间仍存在大量采样空隙^[15]。

近年来，无人机已成为温室气体观测的新兴平台，尤其对于设施尺度的CH₄点源，它具有更高的灵活性和机动性，更快的响应时间和更高的空间分辨率，适合在低空进行活动。无人机测量CH₄的最常见方法包括：(1) 使用无人机收集环境空气样本，落地后使用地面仪器对样本进行分析^[16-21]；(2)通过连接到无人机的管子采样并在地面进行实时分析^[22-23]；(3)用搭载在无人机上的仪器进行实时分析^[24-25]。无人机可以在整个垂直和水平范围内横切羽流，与仅能捕获一小部分羽流的地面移动或固定观测平台相比，减少了对水平和垂直扩散假设的依赖^[21]。但基于无人机平台观测的不确定度依然较大^[18-19]，低空环境中极为重要的大气湍流所造成的不确定度尚无定量分析，同时，无人机不同的飞行安排及策略造成估算结果的差异也尚不清楚。

为此，引入大涡模拟(Large Eddy Simulation，LES)来模拟大气湍流环境中CH₄点源羽流的扩散。大涡模拟采用高时空分辨率数值模式，可以解析再现大气中主要的湍流运动，近年来已有研究成功地将大涡模拟下的模拟羽流应用于各种场景。Viatte等^[26]结合4个移动地基光谱仪的观测数据和实际案例中的LES，引入了优化算法，提升了对地区最大的奶牛场CH₄排放的估算效果；Gaudet等^[27]通过5层嵌套的中尺度天气预报模式(Weather Research and Forecasting，WRF) 降尺度至大涡尺度，评估了WRF大气化学模式(WRF-Chem)对CO₂传输的预测能力；Caulton等^[28]将传统的高斯方法、LES建模和控制释放结合起来，评估了天然气站点下风向CH₄排放估算的不确定度，并发现了不确定度的最大来源是由湍流造成的大气变化性(大气中湍流变量的平均时间不足造成的随机误差)；Ražnjević等^[29]将在非理想条件下进行的现场观测与LES结合起来，验证了利用LES模拟羽流模型直接进行排放速率反演的有效性。

本研究使用LES对1种特定大气条件下已知排放速率的CH₄羽流进行理想模拟，量化了强湍流混合环境下无人机对排放点源使用质量平衡方法(Mass Balance，MB)与逆高斯模型方法(Inverse Gaussian，IG)估算的准确性和不确定度，并探究了不同飞行安排及采样策略对估算结果的影响及改进。

1 实验方法 1.1 大涡模拟

LES可以直接解析含能涡的演变特征，有效捕捉主要的湍流动能，被广泛应用于大气边界层与大气湍流的理论研究中^[30-32]。本研究中，基于WRF模式V3.9.1.1版本中的LES模块，添加了恒定排放速率的被动示踪物代表CH₄排放。该示踪物仅受大气的影响被动传输和扩散，进而形成羽流，对大气无反馈作用。

1.1.1 实验设置

模拟区域为5 km × 3 km × 1.5 km(x，y，z方向)，水平分辨率为10 m，垂直方向网格大小由近地的7.9 m平滑过渡到9.0 m，共设181层，被动示踪物点源被安置在(1 500 m，1 500 m，50 m)位置，排放速率设为恒定的1 000 kg/h。在该分辨率下，除了贴地数个垂直层外，次网格尺度的影响基本可以忽略。模式水平方向采用周期性边界条件，但对于被动示踪物，其作为标量，采用WRF开放式侧边界条件，避免羽流从一侧流出后又从另一侧再次进入模拟区域。模式的次网格湍流方案采用Moeng^[33]的1.5阶湍动能闭合方案。参照Zhou等^[32]以及Shin等^[34]的设置，对模拟施加5 m/s地转风(x轴方向)作为初始条件，地表感热通量$ \overline{w^{\prime} \theta^{\prime}}$为0.2(K·m)/s，以模拟日间对流边界层中的特征，初始位温在925 ~1 075 m存在强逆温，有效限制了边界层的增长。为了使模式中湍流运动发展并达到相对稳定状态，共进行了4 h的模拟，前2 h作为启动时间，后2 h每10 s输出1次模拟结果用于分析，为方便表示，对应的模拟时间设置为00：00—04：00。

1.1.2 模拟CH₄羽流

瞬时态羽流及垂直扰动速度空间分布见图 1(a)(b)(c)。由图 1(a)可见，羽流在50 m排放高度上保持了数十米紧凑状态后，在x轴方向的1 600和1 800 m位置经历2次向上混合，紧接着又急速向下混合，对应图 1(b)水平分布中x轴方向1 700和1 900 m处2次羽流的破裂。结合图 1(c)中100 m高度处垂直扰动速度的水平剖面，可以看到该时刻的羽流受到了(2 000 m，1 500 m)位置处的结构化的下沉气流束及左侧辐条状上升气流的影响；在羽流位置之外，还可以看到范围更广的下沉气流。该结构与热力主导下的日间对流边界层结构特征相一致^[35]，代表了在晴朗日间进行观测的一种典型环境。该条件下充斥了不同尺度的湍涡，湍流混合强盛，使得羽流产生了宏观的蜿蜒行为与微观的相对扩散^[36]。这无疑对无人机采样的空间位置提出了很高的要求，Andersen等^[18-19]的工作已经证实，多数无效的飞行来源于采样轨迹截面的“幕布”未能有效捕捉到羽流，关于采样轨迹截面的位置对估算的影响将在2.2节中讨论。

1.2 基于无人机的主动式AirCore原位观测与模拟采样过程

本研究选择基于无人机的主动式AirCore观测作为研究对象，其详细介绍见文献[17]，这里只作简要概述。主动式AirCore系统是1种空气采样工具，通过1个微型泵将空气吸入到1根长的不锈钢盘管，沿着无人机飞行的轨迹收集样气。泵为KNF020L微型泵，它在1个45 μm针孔的孔口下游使压强低于环境压强的一半，以促成临界流的形成。因此，AirCore的采样流量仅取决于针孔上游的压强(环境压强)，以及针孔和AirCore采集样品的温度。AirCore系统封装在1个碳纤维盒里，固定在无人机螺旋桨的下方，进气口位于碳纤维盒的一侧，因此，采样到AirCore中的空气实际上来自螺旋桨上方，距离螺旋桨不到0.5 m^[37]，由于无人机大部分时间都以1~2 m/s的速度稳定前进，可以用采集到的空气代表飞行轨迹上的样本。获得空气样本后，通过G2401 m光腔衰荡光谱仪(CRDS，美国Picarro公司)对样本进行分析，获取CO₂、CH₄和CO的摩尔分数。在实际采样及分析过程中，分析仪腔内空气样品的混合会造成样气浓度峰值减弱，但可通过对测量信号进行去卷积以获得未受影响的浓度峰值^[17]。在本研究中，不考虑气体混合和扩散造成的浓度稀释。

随后，在大涡模拟资料中模拟这一原位采样过程。首先，须确定采样位置。利用排放源位置100 m高度内的2 h平均风向确定飞行轨迹线的角度，轨迹线与平均风向垂直且交点位于其中心，并在下风向距离点源50，100，150，200 m处进行采样。2 h平均羽流的垂直累积分布和飞行轨迹垂直截面见图 2(a)(b)。

由图 2(a)可见，充足的时间平均几乎消除了湍流混合在水平方向上的随机性，距点源较近的2条采样线较完整地覆盖了平均羽流的扩散范围。之后，确定飞行轨迹及其垂直截面。根据Andersen等^[18-19]及Vinkovic等^[38]的实际采样经验，设定无人机飞行速度恒定为2 m/s，参考实际应用中的续航能力，共飞行10 min 50 s，构建出图 2(b)所示的飞行轨迹，其涵盖的宽度为200 m，高度为100 m(分布在5~105 m之间)的垂直截面。飞行轨迹每次横切的高度间距为20 m，共有6条横切线。浓度的模拟采样频率为1 Hz，即在飞行路径上每隔2 m显示1次数值。以100 m距离为例，图 2(b)展示了采样垂直截面内的平均态羽流的分布，发现该羽流浓度高值在截面中心偏右，但非对称分布，浓度范围中心大致位于截面中心。经统计，以1 min 1次的统计频率，在2 h共121个截面中，仅有5个时刻羽流中心不在该截面中，这表明该采样策略的设置能有效捕获羽流。将上述飞行路径结合下风向100 m采样距离定义为标准飞行策略，后文中以此为基础进行分析，关于不同下风向距离的讨论见2.2节。

关于单次模拟飞行中的采样过程，根据每10 s输出1次的高精度LES三维场，做出以下处理，由1个时刻的瞬时场通过插值确定该时刻及其后9 s(9个时刻)共10个采样数据点位置的值，1次标准飞行过程(10 min 50 s)，由起始时间开始的66个逐10 s的瞬时场信息合成而来，定义为合成路径场，其后通过合成路径场中的浓度及气象要素信息来对排放速率进行估算。需要说明的是，上述处理做出了羽流在10 s内不会发生大幅度蜿蜒行为的假设，作为对比，当模拟采样频率为10 s 1次(0.1 Hz)时，不需要对羽流行为做出假设，但采样点会变得稀疏，经计算二者估算结果类似，0.1 Hz采样时不确定度更大，因此在之后的实验中仍采用1 Hz模拟采样频率。

1.3 量化方法 1.3.1 排放速率

IG方法使用主动式AirCore系统采样得来的三维空间上的CH₄浓度与高斯羽流模型，对排放速率进行反演，该模型计算公式见式(1)。

$ \begin{aligned} & \qquad C(x, y, z)=\frac{Q}{2 {\rm{ \mathsf{ π} }} \sigma_y \sigma_z v \cdot \cos (\theta)} \exp \left[-\frac{1}{2} \cdot\left(\frac{y}{\sigma_y}\right)^2\right] \cdot \\ & \left\{\exp \left[-\frac{1}{2} \cdot\left(\frac{h-z}{\sigma_z}\right)^2\right]+\exp \left[-\frac{1}{2} \cdot\left(\frac{h+z}{\sigma_z}\right)^2\right]\right\} \cdot\\ & \frac{V}{M_{\mathrm{CH}_4}} \end{aligned} $

(1)

式中：C(x，y，z)——给定位置(x，y，z)处的干摩尔分数，mol/mol；(x，y，z)——羽流流向、水平穿过羽流和垂直穿过羽流的投影位置坐标；Q——排放速率，kg/s；σ_y和σ_z——污染物在水平和垂直方向上的扩散参数，m；v——风速，m/s；v·cos(θ)——垂直飞行轨迹剖面的有效速度；h——有效排放高度，m，该高度随监测到的羽流中心位置变化而变化；V——干摩尔体积，m³/mol；M_CH4——CH₄的摩尔质量，为0.016 kg/mol。

关于σ_y和σ_z值的计算，首先根据稳定度参数L(Obukhov长度)和粗糙度参数z₀确定环境的稳定度等级A—F^[39]，在本实验中，L的均值为29.2 m，z₀为0.1 m，确认其稳定度等级为C。其后采用Klug^[40]提出的经验公式$ \sigma_i(x)=R_i x^{r_i}$(i=y，z)确定系数R_i和r_i的值，进而计算不同距离上的σ_y和σ_z。设定高斯模型中先验排放速率(Q_ref)，进而根据公式(1)计算出二维网格平面上的高斯羽流摩尔分数分布(C_IG)，结合模拟采样所得的摩尔分数(C_Observation)，利用公式(2)计算对应的排放速率^[28]。

$ Q=\frac{\sum C_{\text {Observation }}}{\sum C_{\mathrm{IG}}} \times Q_{\text {ref }} $

(2)

MB方法沿用Andersen等^[18]的描述，使用克里金插值法将飞行轨迹上的浓度及气象数据插值到垂直的二维网格平面，之后利用公式(3)求得排放速率。MB方法的计算公式见式(3)。

$ Q=\frac{\Delta X \cdot M_{\mathrm{CH}_4}}{R} \sum\limits_i^{k_i} \sum\limits_j^{k_j} \frac{C_{i, j} \cdot P_{i, j} \cdot v_{i, j}}{T_{i, j}} $

(3)

式中：

ΔX——每个网格的面积，为2 m × 2 m；R——通用气体常数，为8.314 5(kg·m²)/(s²· K·mol)；k_i，k_j——克里金平面中水平网格和垂直网格的数量；C_i，j——格点处CH₄的摩尔分数，mol/mol；P_i，j——每个格点的压力，Pa；v_i，j——格点处风速，m/s；T_i，j——格点处温度，K。

CH₄摩尔分数分布见图 3(a)(b)(c)，其中IG方法的浓度数值为当Q_ref设置为100 g/s时C_IG的分布。

1.3.2 不确定度

如1.1节所述，可供分析的2 h LES资料每10 s保存1次，因此共有721个时刻的三维场资料得以保存。每次模拟飞行采样均能获取1个合成路径场，并通过IG与MB方法分别获取1个排放速率量化值。从02：00时刻起，每隔1 min进行1次模拟飞行实验，采用1.2节中所描述的标准飞行策略，由于每次飞行耗时10 min 50 s，至03：49时刻，共进行110次独立模拟飞行实验，构成集合平均的样本集，以该序列标准差表示为不确定度，为方便比较，用均值对其做标准化，用百分数表示。

2 结果与讨论 2.1 排放速率及不确定度的量化

110次独立模拟飞行实验通过IG和MB 2种方法得到的估算结果见图 4。

由图 4可见，IG与MB方法均能有效量化排放速率，在标准飞行条件下，所估算的排放速率值均值分别达到了模拟排放真值的95.3%和86.1%，2种方法的不确定度相似，分别为56.6%和56.9%。

IG和MB 2种方法估算结果减去均值后的概率密度分布见图 5(a)(b)。由图 5可见，二者均基本符合正态分布，并通过了K-S检验，这也验证了2种方法的估算结果集合具有代表性，意味着如果在该环境中进行任意1次飞行观测，根据正态分布的概率密度分布，估算值有大约68%的概率会落在(953 ± 539)kg/h范围内(以IG方法为例)，这与图 4观察到的连续多次实验结果的变化趋势一致。IG与MB估算结果具有高度相关性，相关系数达0.95，但可以发现MB量化结果在大多数情况下低于IG。该情形的1种可能解释是MB方法对飞行范围更为敏感，如果羽流只有部分处于飞行截面中，由于MB方法的机制会造成结果的低估，但这无法解释采样时期羽流蜿蜒行为未超出飞行范围时估算值依然低于IG的情况，其涉及的具体机制有待进一步探究。同时，估算值的不确定度较高、极大值与极小值间差异较大，体现出了在IG与MB 2种方法的计算过程中所涉及物理量较强的变化性，该变化性可能与稳定度及湍流混合强度有关。针对较强的变化性，对物理量进行长时间平均可以有效消除高频信号。

1架无人机通过落地后重复飞行采样可获取更长的观测时间，多次复飞视为1次整体飞行安排，可对观测和估算结果进行时间平均以接近稳态，实际实验中可通过落地后更换电池再进行复飞实现。本研究对比了飞行1次到重复飞行5次的量化结果。对于多次重复飞行，飞行起飞间隔为12 min，n(2~5)次重复飞行视为1次整体的模拟飞行实验，在2 h可用的模拟时间内，相应共进行了98，86，74和62次实验，以该样本的估算值进行分析。结果表明，连续进行5次飞行可使不确定度从56.6%(IG)和56.9%(MB)降低至＜30%，分别为28.3%(IG)和26.4%(MB)。对于多次飞行的数据，有2种处理方案，第1种直接采用图 4中单次模拟飞行的连续估算数据，1组多次重复飞行的估算结果为对应的多个飞行次数估算值的平均值；第2种方案对多次采样得到的浓度、风速等必要物理量计算平均值，之后再进行量化估算。经对比，2种方法结果无显著差异，因此选择第1种更易实施的方案进行实验。重复飞行的量化结果见表 1，重复飞行估算结果的不确定度随飞行次数的变化见图 6。由表 1和图 6可见，经过多次飞行，2种方法的均值几乎不变，但不确定度显著减小，2次飞行相对1次有较大的降幅，分别为14.4%(IG)和16.6%(MB)，其后下降速度趋缓，整体上，平均每多1次飞行，不确定度下降7.1%(IG)和7.6%(MB)。

2.2 飞行策略对估算的影响

本节旨在探究空间上不同的飞行策略对排放估算的影响，为此首先考虑了基础的单架无人机观测，从而同时在空间上获取更多的采样信息，之后评估了2架无人机的同步观测。

由于羽流随下风向逐渐扩散，不同的下风向距离会使得获取的浓度分布出现差异，进而影响估算结果。在此开展关于飞行轨迹距排放源距离的敏感性实验，量化结果见图 7(a)(b)。由图 7可见，4种采样距离均出现总体分布相对模拟排放真值偏低的情形，当距离达到150 m时，真值高于估算数据分布的上四分位数，距离200 m时整体分布进一步降低。而对于50和100 m距离，二者分布类似，但100 m处2种量化方法的均值及中位数均都更加接近真值，以此为指标，认为100 m相对最优。这解释了标准飞行策略采用100 m距离的原因。适当的距离与风速有直接关联，本研究中10 ~100 m距离间的风速在3~4 m/s。Andersen等^[18]在实际观测中发现，当风速＜2 m/s时会造成估算结果极大的不确定度，其原因与风速较小使得羽流更为弥散，难以有效观测且风向变化性强有关，关于其定量分析有待于进一步探究。在确定了合适的距离(100 m)之后，针对无人机有限的采样范围，设计了向右、向左及向上偏移的垂直采样截面，对应图 2(b)中截面在水平方向上沿图 2(a)中原方向线向两侧各移动100 m，在水平方向上向上移动50 m，仅覆盖一部分平均羽流。由图 7(b)可见，在截面中心偏移平均风向半个截面长度或高度后，IG方法估算均值能达到标准状态的71.8%~ 83.8%，而MB方法估算均值仅能达到标准状态的51.7%~64.3%。

相比2.1节中单架无人机连续多次飞行拓展了时间范围，多架次无人机同步飞行会缓解无人机非瞬时采样造成的对蜿蜒羽流的捕获缺失，理论上会有利于获取更好的估算结果^{[18, 38]}，本节对此假想进行定量分析。为此设计了3种基于2架无人机同步采样的飞行策略，包括与标准飞行方法相同的横切距离(100_20策略)、与标准飞行近似采样时间内加密的横切距离(110_10策略)，以及为了便于对比，1种理想化的2架无人机近距离同步采样(110_10 a策略)，其中下划线前首个数字表示截面高度范围，第2个数字表示横向飞行轨迹间距，需要说明的是110_10 a策略由于2架无人机距离较近(10 m)，在实际操作中会产生较大的相互干扰，可以通过其中1架稍滞后飞行消除干扰。标准采样及3种双无人机飞行策略的飞行路径见图 8(a)—(d)，不同飞行策略的量化结果见表 2。

由图 8和表 2可见，空间同步采样会对IG方法造成较大的干扰，使其量化结果的不确定度增加，但IG方法对2架无人机的不同飞行策略并不敏感；对于MB方法则较为有效，能将不确定度范围降低至35.2%~51.9%，小于标准飞行的56.9%。关于其中原因，发现空间同步采样易使得飞行剖面出现多个浓度中心，这在很大程度上干扰了IG方法对羽流的正确拟合，也说明用IG方法模拟瞬时羽流缺少理论支持。以第2 h 7 min开始的1次模拟飞行为例，见图 9(a)(b)。

在共计5 min 20 s的采样时间内，一开始羽流中心大概位于(0 m，60 m)位置，其后多次转移均被飞行路径获取，出现多个中心，IG羽流在此例中被拟合至左上方插值结果浓度最高值处，并不能代表该采样时间内羽流的整体分布状态。但多中心并未造成MB方法的高估，这可能与每次横向飞行过程均未充分穿过整个羽流有关。

3 结论

使用基于无人机的原位观测系统对CH₄点源排放速率进行估算，使用大涡模拟来复现典型湍流混合环境中的CH₄羽流与无人机原位采样过程，经实验与分析，得出以下结论：

(1) CH₄羽流在湍流混合较强的环境中受不同尺度湍涡的影响，有较强的蜿蜒行为，为应对这一情况，应尽量选择可以涵盖羽流变化范围的飞行区域及合适的与排放源之间的距离，为此定义了标准飞行策略以供后续定量分析。

(2) 通过连续多次独立模拟飞行实验中IG和MB 2种方法的估算结果集合，得出估算的排放速率均值分别达到了模拟排放真值的95.3%和86.1%，2种方法的不确定度类似，分别为56.6%和56.9%，估算值呈现近似正态分布的特征。

(3) 为减小湍流造成的较大的不确定度，提升估算效果，1架无人机的重复飞行可获取更长的观测时间以对观测和估算结果进行时间平均，多次重复飞行通过对结果进行平均可以使不确定度显著减小，连续进行5次飞行可使不确定度从56.6%(IG)和56.9%(MB)降低至＜30%，分别为28.3%和26.4%。

(4) 相对于标准飞行策略，飞行轨迹截面中心偏移半个截面长度或高度后，估算均值为偏移前的71.8%~83.8%(IG)及51.7%~64.3%(MB)。在50~200 m 4个由近及远的下风向距离的估算均值验证了标准飞行策略下的100 m方案的准确度最高。但在实际观测中难以量化飞行轨迹截面的偏移程度。

(5) 多架次无人机同步飞行会缓解无人机非瞬时采样造成的对蜿蜒羽流的捕获缺失，在空间上可以同时获取更多采样信息。实验结果表明，空间同步采样会对于IG方法造成较大的干扰，使其量化结果的不确定度增加；对于MB方法则较为有效，能将不确定度范围降低至35.2%~51.9%，小于标准飞行的56.9%。相较于单架无人机重复飞行，该方法对不确定度的降低相对不显著。

本研究量化了1种特定湍流混合条件下基于无人机观测估算的不确定度，并通过不同的飞行安排及飞行策略对其进行了有效的限制，证明了无人机观测的可靠性与较大的可提升空间。但值得说明的是，不同程度的湍流混合强度及稳定度会在很大程度上改变羽流扩散状况及形态，进而会对估算准确度及不确定度造成影响，这对于无人机观测的飞行策略及估算模型都提出了更大的挑战。本研究仅考虑了一种典型日间对流边界层中的环境，且未对羽流参数进行量化比较，对于多种湍流混合条件下的无人机观测的适用性及相关改进仍须进一步研究。此外，本研究仅考虑CH₄的传输扩散过程，因而该结果也适用于其他被动示踪物的点源排放估算。