陆地生态系统具有吸收二氧化碳和缓解气候变化等重要功能(Janssens et al.，2003)。为了评价陆地生态系统对气候变化的贡献，开展陆地生态系统与大气之间碳、水和能量交换监测研究至关重要，涡度相关技术是目前应用较广泛的监测方法(于贵瑞等，2006a；Running，2008；Beer et al.，2010；Zhao et al.，2011)。现有的通量观测研究网络已涵盖全球各个地区，如美国通量网、欧洲通量网、亚洲通量网等。由于受系统故障和降雨等因素影响，通量塔观测数据存在缺失现象，缺失比率高达20%～60%，夜间数据丢失比例较高(Moffat et al.，2007)，上述问题使得直接获取连续完整的日、月、季、年尺度净生态系统交换量(NEE)、总初级生产力(GPP)和生态系统呼吸速率(Re)碳通量各分量困难。为了获取完整和可靠的通量数据，需要采取合理的插补方法对缺失数据进行插补。

在通量观测数据插补和分量分解上已做了大量研究(Falge et al.，2001；Moffat et al.，2007；Desai et al.，2008)。常用的插补和分量分解方法有平均昼夜变化、查表法和非线性回归等(Falge et al.，2001；于贵瑞等，2006b)。非线性回归方法是根据一定时间段内(窗口)有效数据，建立环境因子与有效数据之间的模型，然后据此插补缺失的数据。Falge等(2001)指出建立一套标准的数据插补方法可以增强区域和全球通量网络不同站点之间的可比性。缺失数据插补方法之间的差异及其对碳通量各分量估算结果的影响已得到深入研究(Moffat et al.，2007；Desai et al.，2008；刘敏等，2010；黄昆等，2013)。通过模拟数据缺失情景，比较15种插补方法，结果表明，插补方法的精度与缺失时间长度和数据累计时间尺度(半小时或日尺度)等因素密切相关，数据缺失时间越长，插补难度越大及精度越低(Moffat et al.，2007)。但是，参数估计窗口对插补结果影响研究相对较少。在确定合适的参数估计窗口时，需要综合考虑通量塔下垫面季节变化特征和数据缺失时间长度。当通量数据季节变化明显时，通常采用15天、1月、双月或季度作为参数估计窗口拟合模型参数(Hollinger et al.，2004；Richardson et al.，2006；Noormets et al.，2007；Desai et al.，2008)。采用较小参数估计窗口可以得到较好的插补效果(Reichstein et al.，2005；Lasslop et al.，2010)。但在数据缺失时间长度较大时，过小的参数估计窗口将导致模型参数无法合理估计以及插补结果远远偏离实际值等问题。因此，如何确定合适的参数估计窗口对插补结果的可靠性至关重要。

本研究按Reichstein等(2005)和Lasslop等(2010)对2011年安吉毛竹(Phyllostachys edulis)林生态系统碳通量数据进行插补和分量分解： 基于夜间数据方法(NB)和基于白天数据方法(DB)，分析不同参数估计窗口对插补结果的影响，进而比较上述2种方法在应用于毛竹林生态系统通量观测数据插补和分量分解上的差异，为获取详实和连续完整的毛竹林NEE，Re和GPP数据及其时空变化研究提供方法。

安吉县位于浙江省西北部(119°14′—119°52′ E，30°23′—30°53′ N)，全县面积1 886.45 km²。属亚热带海洋性季风气候，光照充足、气候温和、雨量充沛、四季分明，年均降水量1 400 mm，年均气温15.6 ℃。森林覆盖率达到69.4%，拥有山林13.2万 hm²，其中竹林面积6.33万 hm²。毛竹林面积为4.99万 hm²，占林地总面积的37.8%，为著名的“中国竹乡”。土壤主要是发育于酸性岩浆岩和沉积岩的红壤土类(谢莉等，2012)。毛竹林通量观测塔位于安吉县东南部的山川与天荒坪2个乡镇交界处。以观测塔为中心1 000 m范围内主要森林类型为毛竹林，其平均树高为11 m，林分平均密度为3 235株 ·hm^-2。此外，还有少量针阔混交林、农田、城镇和道路等地类。通量观测塔概况参见孙成等(2013)。

本研究中，缺失数据时间窗口指数据连续缺失的时间长度；参数估计窗口指用于拟合模型参数所需数据的时间长度；移动窗口指参数估计窗口移动的时间长度。以2011年安吉县毛竹林生态系统通量塔半小时尺度观测数据为研究数据，数据集共17 520行。依据净太阳辐射能( R_n)将数据分为白天和夜间数据，当R_n>0时为白天数据，R_n≤0时为夜间数据，其中白天数据7 822个，占44.65%；夜间数据9 698个，占55.35%。通过异常值剔除(Li et al.，2005)及夜间摩擦风速校正等处理后，夜间有效数据2 202个，仅占夜间总数据量的22.7%。白天有效数据4 953个，占白天总数据量的63.3%。白天数据的有效比例明显高于夜间数据。白天缺失数据的时间窗口平均长度为3.9 h，长度大于24 h的有15个，主要集中在6—9月份，最大缺失数据的时间窗口长度为4.6天。夜间缺失数据时间窗口平均长度为4.2 h，长度大于24 h的有28个，最大缺失数据的时间窗口长度为6天。当忽略两个缺失窗口间隔内小于3 h的有效数据，最长缺口达24天。可见，缺失数据量和缺失时间长度较大，数据插补面临挑战性。

采用夜间半小时尺度观测数据对Arrhenius呼吸方程(公式1)进行参数化，然后根据参数化后的呼吸方程插补夜间生态系统呼吸速率(Ren)，并将白天空气温度代入公式1中估算白天生态系统呼吸速率(Red)(Reichstein et al.，2005；Lasslop et al.，2010)。

式中：R_ref为当参考温度T_ref等于15 ℃时的呼吸量(Lasslop et al.，2010)；E₀为活化能值；T_air为气温，℃；T₀为常数-46.02 ℃(Lloyd et al.，1994)。

采用白天生态系统净交换量(NEEd)和光合有效辐射量对Michaelis-Menten光响应曲线模型(公式2)参数化，然后根据参数化后的模型对NEEd缺失数据进行插补。

式中： α为光响应曲线斜率，即光能利用效率；β为光饱和时最大光合速率；PAR为半小时间隔内光合有效辐射平均值。

将Red作为光响应曲线模型的截距，采用白天数据对模型(式3)进行参数化，然后根据参数化后的模型对NEEd和夜间呼吸速率缺失数据进行插补

当缺失数据窗口长度超过参数估计窗口或窗口内有效数据过少时，会出现模型参数拟合值不合理等问题(Reichstein et al.，2005；Lasslop et al.，2010)，此时采用如下处理方法进行参数拟合： 1)固定生物学含义不明确的参数，避免参数之间高度相关，提高模型的稳定性。参考文献Reichstein等(2005)，采用90天参数估计窗口和15天移动窗口估计E₀值，选择满足1倍标准差范围内所有E₀值并计算均值，将参数E₀设置为常数(181.45K)；2)当参数估计窗口内样本个数少于6时，模型所有参数采用离该窗口最近的上一个窗口参数值代替；3)当参数估计窗口样本个数大于6个时，采用非线性曲线拟合方法对模型进行参数化，当参数值等于给定的最大或最小值即临界值(表 1)时，则采用上一个窗口参数值代替，其他参数重新拟合，如果仍有参数值等于临界值，则所有参数采用上一个窗口参数值代替。参数估计窗口内的有效数据用于参数拟合，然后根据参数化后的模型对移动窗口内的缺失数据进行插补。考虑夜间数据缺失量较多，NB方法的呼吸模型插补时，移动窗口设定为15天；光合响应曲线模型插补时，移动窗口设定为2天。DB方法模型插补时，移动窗口设定为2天。为了分析不同参数估计窗口大小对插补结果的影响，本研究比较了以下15个参数估计窗口： 2，4，15，30，60，90，120，150，180，210，240，270，300，330和360天。模型参数初始值及临界值见表 1。采用参数化后的模型对缺失数据进行插补，计算预测值与实测值的均方根误差RMSE(公式4)和相对系统偏置Biasr(公式5)，分析碳通量各分量时间序列变化趋势，评价参数估计窗口对插补结果的影响。

式中： N为样本总数；${{{\hat y}_i}}$，y_i和${\bar y}$分别为第i个样本的预测值、实测值和实测值均值。

表 1 模型参数初始值和临界值 Tab. 1 Initial values and their range of model parameters

表 1 Initial values and their range of model parameters Tab. 1 模型参数初始值和临界值 参数 Parameters 初始值 Initial value 区间 Range E0 100 [50，450] Rref 0.25 [0，3] α 0.003 [0，0.5] β 0.5 [0，5]

采用实测数据对模型进行参数化，根据参数化后的模型获取预测值，然后与实测值相比计算模型插补精度。当参数估计窗口为15天时，对Ren的插补结果误差最小，RMSE为0.215 mg CO₂ ·m^-2 s^-1，Biasr为-0.4%(图 1)。随着窗口增大，Ren插补结果的RMSE有所增加但不明显，Biasr呈先降后升的趋势(图 1)。

	图 1 NB方法在不同参数估计窗口条件下对Ren缺失数据插补精度变化 Fig.1 Gap-filling accuracy of Ren with different window sizes for parameter estimation based on NB method

虽然当参数估计窗口大小为15天时误差最小，但是Ren预测值变化趋势出现异常波动(图 2)，随着窗口从15天增大到360天(图 2)，Ren预测值局部波动幅度降低，呈单峰变化，冬季Ren速率较低和夏季高的常态趋势。 随着窗口增大，夏季Ren预测值不断增大而冬季降低(图 2)。另一方面，窗口越小，局部Ren值波动幅度越大，表明减小窗口更有利于反映短时间区间内呼吸速率的差异性。在确保Biasr尽量小(图 1)和体现短时间区间内Ren波动性的情况下，用于插补Ren缺失数据和预测Red的最优参数估计窗口为90天(图 2)，此时Ren缺失数据插补结果的Biasr为-8.41%，从Biasr为负值可推出其年尺度累计值存在系统低估。

	图 2 NB方法在不同参数估计窗口(15，60,90,360天)条件下Ren插补结果 Fig.2 Ren with different window sizes (15,60,90,360 d) for parameter estimation based on NB method

根据Arrhenius呼吸模型预测Red，利用公式2对白天缺失数据进行插补。移动窗口为2天，不同参数估计窗口的插补效果见图 3。当参数估计窗口为4天时，插补误差最小，RMSE为0.253 mg CO₂ · m^-2 s^-1及Biasr为-1.89%。当移动窗口为2天和参数估计窗口为4天时插补结果最优。随着参数估计窗口增大，RMSE增加(图 3)。NEEd预测值及其5日滑动平均值的变化趋势(图 4)表明随着窗口从4天增大到360天(图 4)，季节性波动越不明显。在体现季节波段性的前提下，参数估计窗口为4天时的NEEd插补结果相对较好(图 4)。

	图 3 NB方法在不同参数估计窗口条件下对NEEd缺失数据插补的精度变化 Fig.3 Gap-filling accuracy of NEEd with different window sizes for parameter estimation based on NB method

	图 4 NB方法在不同参数估计窗口(4，15，90，360天)条件下NEEd插补结果 Fig.4 NEEd with different window sizes (4,15,90,360 d) for parameter estimation based on NB method

当移动窗口为2天及参数估计窗口为4天时，模型误差最小，NEEd插补结果的RMSE和Biasr分别为0.256 mg CO₂ ·m^-2 s^-1和-3.15%(图 5)。随着窗口增大，NEEd插补结果的RMSE增加并趋向稳定，但Biasr呈起伏变化，总系统偏差存在高估或低估(图 5)。DB方法在不同参数估计窗口条件下Ren插补结果变化趋势(图 6)表明参数估计窗口过小时，Ren预测值出现明显异常(图 6)，相应的Red也出现明显异常。DB方法在不同参数估计窗口条件下NEEd插补结果变化趋势(图 7)分析表明当参数估计窗口大于60天时，DB方法的NEEd预测值变化趋势较合理(图 7)，窗口越大预测值变异程度越小。综合考虑模型预测精度和更好地保留变异信息，参数估计窗口为60天时插补结果最优，此时的NEEd插补结果RMSE和Biasr分别为0.277 mg CO₂ ·m^-2 s^-1和-3.06%(图 5)，从Biasr为负值可推出其年尺度净交换量存在系统低估。

	图 5 DB方法在不同参数估计窗口条件下对NEEd缺失数据插补的精度 Fig.5 Gap-filling accuracy of NEEd with different window sizes for parameter estimation based on DB method

	图 6 DB方法在不同参数估计窗口(4，15，60，360天)条件下Ren插补结果 Fig.6 Ren with different window sizes (4,15,60,360 d) for parameter estimation based on DB method

	图 7 DB方法在不同参数估计窗口(4,15,60,360天)条件下NEEd插补结果 Fig.7 NEEd with different window sizes (4,15,60,360 d) for parameter estimation based on DB method

根据NB和DB方法最优参数估计窗口插补和分量分解得到的日尺度的NEE，GPP和Re散点图分析表明，2种方法的NEE和GPP结果非常接近，R²分别为0.73和0.83，无明显的系统偏差。NB方法的GPP稍高于DB方法的GPP，主要原因是由于NB方法估算的Re明显高于DB方法(图 8)。

	图 8 NB和DB方法插补的日尺度碳通量分量散点图 Fig.8 Scatter plots between different carbon flux components by using NB and DB methods

NB和DB方法插补和分量分解得到的年尺度碳通量各分量分析表明，NB方法的GPP和Re高于DB方法的结果，分别高出13.8%和26.8%，而NEE低于DB方法32.2%(表 2)。NB和DB方法2种方法插补得到的NEEd相差最小，NB方法约高出DB方法4.6%(56 g C · m^-2 a^-1)，而对于Ren，NB的结果高出DB方法23.7%(190.5 g C ·m^-2a^-1)。

表 2 NB和DB方法之间插补得到的年尺度碳通量各分量比较^① Tab. 2 Comparisons of carbon flux components estimated by using NB and DB methods

表 2 NB和DB方法之间插补得到的年尺度碳通量各分量比较① Tab. 2 Comparisons of carbon flux components estimated by using NB and DB methods g C ·m-2 a-1 方法Method GPP NEE Re Ren Red NEEd NB -1 899.7 -417.5 1 482.2 802.5 679.7 -1 220.0 DB -1 637.0 -552.0 1 085.0 612.0 473.0 -1 164.0 ①负值表示吸收，正值表示排放。Positive mean absorption，neqative mean emission.

缺失数据插补时，不同参数估计窗口对插补后的结果具有明显的影响，窗口越小越能体现碳通量短期变化幅度(Lasslop et al.，2010)，但在缺失数据比率较高时结果出现异常现象。采用NB方法进行Re插补时，最优移动窗口和参数估计窗口分别为15天和90天，GPP插补时最优移动窗口和参数估计窗口分别为2和4天，而DB方法的最优移动窗口和参数估计窗口分别为2和60天，本研究采用的参数估计窗口明显大于Lasslop等(2010)所采用的参数估计窗口4～15天，这主要原因是Lasslop等(2010)所用到的数据有效比率高达80%且最长缺失时间长度不超过750 h，其有效数据比率明显高于本研究采用的数据而其缺失时间长度明显小于本研究所采用的数据，可见参数估计窗口最优值与数据缺失量及其缺失时间长度密切相关。无论是NB还是DB方法，结果表明随着参数估计窗口增大，夏季呼吸速率不断增大而冬季呼吸速率降低，这主要原因是窗口越大，窗口内包含的呼吸速率实测值和温度值域变化范围越大，相比小窗口，大窗口的呼吸模型曲线陡度要大于小窗口的模型曲线，因此模型曲线斜率越大，会导致温度低处(冬季)越低，温度高处(夏季)越高的结果。白天缺失数据插补时的最优参数估计窗口明显小于夜间缺失数据插补时的最优参数估计窗口，主要是白天数据的缺失比率较低，使得在较小的参数估计窗口内有足够多的样本用于参数估计，使得模型参数拟合值更加稳定和准确，降低模型不确定性。在数据缺失比率较高的前提下，采用过小的参数估计窗口将使窗口内样本过少，进而引起呼吸方程参数不确定性大，最后造成呼吸速率预测值出现明显异常(图 2和6)。缺失数据比率对参数估计窗口大小的选择具有重要的影响，因此本研究得出的最优窗口只适用于本研究案例，同时可作为安吉县毛竹林通量塔观测数据插补时的参考值。

在碳通量各分量统计结果上，NB和DB方法存在差异。NB方法的GPP稍高于DB方法的GPP(图 8)，该结果与Lasslop等(2010)和黄昆等(2013)相一致。NB方法基于通量塔实测的夜间呼吸数据对呼吸方程进行拟合，更有利于真实反映Re对温度变化的敏感性(Lasslop et al.，2010)。 另一方面，NB方法拟合的R_ref大于DB方法的R_ref，夜间呼吸对温度的响应与白天呼吸对温度的响应方式有差异，将NB方法采用夜间呼吸拟合的参数用于预测白天呼吸，可能存在高估的问题(黄昆等，2013)。NB方法的GPP和Re高出DB方法结果的比例与黄昆等(2013)结果相似。从碳通量各分量对比分析得出，NB和DB方法之间NEEd非常接近，主要原因是白天观测数据缺失比率较低(刘敏等，2010)，因此导致2者方法之间插补结果差异的主要因素是呼吸速率分量不同，这一方面与夜间观测数据缺失比率过高造成两者方法模型参数不确定性较大密切相关。另一方面，由于昼夜生态系统呼吸温度敏感性具有差异性，采用NB方法时将高估白天生态系统呼吸(Lasslop et al.，2010)。总之，NB和DB方法各有弊端，NB方法在插补夜间呼吸速率缺失数据时采用实测夜间呼吸速率数据对呼吸速率方程参数进行拟合，Re结果更能反映实际变化，但当实测夜间呼吸速率数据缺失比率较高时，会导致拟合的参数不确定性大。而DB方法中模型需要拟合的参数过多且参数之间存在高度自相关，使得呼吸方程参数难以得到合理约束和拟合，缺乏生物学意义，得到的结果同样存在较大不确定性。在应用时，需要综合考虑数据缺失比率和下垫面的碳通量季节变化特征，确定合适的参数估计窗口，结合2种方法减少插补结果不确定(黄昆等，2013)。