周 宇，黃 輝，張勁松**，孟 平，孫守家

（1.中國林業(yè)科學(xué)研究院林業(yè)研究所，北京 100091；2.國家林業(yè)和草原局林木培育重點(diǎn)實(shí)驗(yàn)室，北京 100091；3.南京林業(yè)大學(xué)南方現(xiàn)代林業(yè)協(xié)同創(chuàng)新中心，南京 210037）

渦度相關(guān)法是目前測定大氣與植被群落間凈CO2交換（NEE）最直接、理論與技術(shù)發(fā)展最為迅速的一種微氣象學(xué)方法，得到微氣象學(xué)和生態(tài)學(xué)界的廣泛接受和認(rèn)可，廣泛應(yīng)用于全球不同陸地生態(tài)系統(tǒng)的物質(zhì)與能量交換研究中，其觀測數(shù)據(jù)也經(jīng)常應(yīng)用于各種模型的檢驗(yàn)和驗(yàn)證[1]。由于渦度相關(guān)觀測儀器需要定期進(jìn)行校準(zhǔn)和維護(hù)，觀測過程中經(jīng)常遇到難以避免的系統(tǒng)故障（儀器故障或供電系統(tǒng)故障等）和外界干擾（暴雨、連續(xù)陰雨、霧霾等惡劣天氣條件等），渦度相關(guān)觀測系統(tǒng)的原始數(shù)據(jù)往往會(huì)出現(xiàn)缺失；而在對(duì)通量數(shù)據(jù)進(jìn)行數(shù)據(jù)質(zhì)量控制后，如進(jìn)行異常值剔除、低湍流條件（如夜間大氣層結(jié)穩(wěn)定等）下數(shù)據(jù)的剔除等，觀測數(shù)據(jù)又會(huì)出現(xiàn)新的缺失。最終，年通量觀測數(shù)據(jù)往往存在20%～65%的缺失，其中還可能出現(xiàn)較長時(shí)段的連續(xù)缺失（長達(dá)半個(gè)月，甚至一個(gè)月）[2?4]。為了獲取完整和可靠的通量數(shù)據(jù)，需要采取合理的插補(bǔ)方法對(duì)缺失數(shù)據(jù)進(jìn)行插補(bǔ)[5]。

目前，常用的數(shù)據(jù)插補(bǔ)方法主要有平均晝夜變化法（Mean Diurnal Variation，MDV）、查表法（Look-Up Table，LUT）、非線性回歸法（Non-Linear Regression，NLR）、邊際分布采樣法（Marginal Distribution Sampling，MDS）和人工神經(jīng)網(wǎng)絡(luò)法（Artificial Neural Network，ANN）等[3,6?15]。選擇不同的插補(bǔ)方法會(huì)產(chǎn)生不同的結(jié)果[1]。MDV 方法不依賴通量和環(huán)境變量的函數(shù)關(guān)系，插補(bǔ)夜間數(shù)據(jù)適合7d 的窗口，而日間適合14d 的窗口，在極端晴天或陰天條件下容易產(chǎn)生估算偏差[6]。LUT 方法通過氣象變量建立索引表估計(jì)NEE，日間最適合使用氣溫，夜間適合使用土壤溫度，而NLR 方法則利用生態(tài)系統(tǒng)呼吸方程和光響應(yīng)曲線估算NEE。LUT 和NLR 方法的各個(gè)環(huán)境變量的分組越詳細(xì)，插補(bǔ)效果越好，但均受到建立索引表或回歸方程時(shí)所選取的時(shí)間階段、站點(diǎn)光響應(yīng)曲線的離散性、云量和干旱等因素的影響[6]。MDS 方法綜合了MDV 與LUT方法，自動(dòng)延長插補(bǔ)窗口，不需假設(shè)環(huán)境響應(yīng)方程和預(yù)設(shè)初始值，被歐洲通量網(wǎng)和FLUXNET（全球長期通量觀測網(wǎng)絡(luò)）用作標(biāo)準(zhǔn)化處理[10]。ANN 方法在處理不同數(shù)量級(jí)和非連續(xù)變量方面具有巨大的優(yōu)勢，可通過更多的變量對(duì)NEE 進(jìn)行插補(bǔ)，但在某些情況下輸出結(jié)果可能不穩(wěn)定[7]。是否進(jìn)行u*校正對(duì)計(jì)算結(jié)果影響很大，u*校正通常使年總NEE更偏正[6]。保留NEE 與氣象變量間基本生態(tài)學(xué)響應(yīng)關(guān)系的方法能夠取得較為認(rèn)同的插補(bǔ)結(jié)果。綜上，以往針對(duì)通量數(shù)據(jù)缺失插補(bǔ)方法的研究，大多關(guān)注不同插值方法在較短時(shí)間缺失情景下的插補(bǔ)性能，對(duì)不同方法在長時(shí)段且連續(xù)缺失情景下的插補(bǔ)精度和穩(wěn)定性及其能夠取得較好插補(bǔ)效果的缺失范圍關(guān)注較少，且一直未形成通用的缺失數(shù)據(jù)插補(bǔ)方案。

本研究以華北低丘山地栓皮櫟人工混交林生態(tài)系統(tǒng)為例，以基于渦度相關(guān)法觀測得到的NEE數(shù)據(jù)為研究對(duì)象，在隨機(jī)生成連續(xù)1、3、7、15和31d 數(shù)據(jù)缺失的情景下，選取平均晝夜變化法（MDV）、查表法（LUT）、非線性回歸法（NLR）、邊際分布采樣法（MDS）和人工神經(jīng)網(wǎng)絡(luò)法（ANN）分別對(duì)50 個(gè)缺失數(shù)據(jù)集進(jìn)行NEE 數(shù)據(jù)插補(bǔ)，探究不同插值方法在不同缺失情景下的插補(bǔ)精度和穩(wěn)定性，并評(píng)估不同插值方法對(duì)缺失片段長度的敏感性，探討不同插補(bǔ)方法所適用的連續(xù)數(shù)據(jù)缺失時(shí)長，以期為山地森林生態(tài)系統(tǒng)渦度相關(guān)法通量觀測數(shù)據(jù)插補(bǔ)方法的選擇提供參考，為準(zhǔn)確估算區(qū)域碳收支、預(yù)測氣候變化對(duì)碳儲(chǔ)存及碳匯的影響、深入量化凈生態(tài)系統(tǒng)碳交換等提供理論基礎(chǔ)。

1 數(shù)據(jù)與方法

1.1 基準(zhǔn)數(shù)據(jù)集

通量觀測地點(diǎn)位于黃河小浪底森林生態(tài)系統(tǒng)國家野外定位觀測研究站南山觀測區(qū)（35°01'45''N，112°28'08''E，海拔410m）。該站位于河南省濟(jì)源市轄區(qū)，地處太行山南麓與黃河中游的交錯(cuò)帶。觀測對(duì)象為栓皮櫟（Quercus vari abilisBl.）、側(cè)柏[Platycladus orie ntalis(L.) Franco]和刺槐（Robinia pseudoacacia）人工混交林生態(tài)系統(tǒng)，造林時(shí)間分別為1972年、1974年和1976年，平均株高分別為10.5、8.2 和9.3m，主要樹種為栓皮櫟（約占80%）。研究區(qū)屬暖溫帶亞濕潤季風(fēng)氣候，年平均氣溫為12.4℃。歷年平均降水量為641.7mm，受季風(fēng)氣候影響，降水季節(jié)分配不均勻，5?9月平均降水量為438.0mm，占全年的68.3%[16]。生長季風(fēng)向以東北偏東、西南方向?yàn)橹鳌?/p>

渦度相關(guān)通量觀測系統(tǒng)由CSAT3 型三維超聲風(fēng)速儀（USA）、LI-7500A 型開路式CO2/H2O 氣體分析儀（USA）組成，安裝高度距地面36m。數(shù)據(jù)采集器為SmartFlux 型（USA），原始數(shù)據(jù)采樣頻率為10Hz[17]。

小氣候梯度觀測系統(tǒng)包括7 層HMP155 型溫濕度傳感器（安裝高度分別為5、8、11、14、18、26 和32m），6 層WindSonic 型二維超聲風(fēng)速計(jì)（安裝高度分別為5、8、11、14、18 和26m），CNR1 型四分量輻射計(jì)（安裝高度為17m），3 層土壤溫濕度傳感器（安裝深度分別為 5、10 和20cm），土壤熱通量計(jì)2 個(gè)（分別位于塔的東西兩側(cè)）[18]。數(shù)據(jù)采集器為CR1000 型，數(shù)據(jù)存儲(chǔ)間隔為10min。

選取2017年3月1日?11月30日10Hz 通量觀測數(shù)據(jù)，經(jīng)渦度協(xié)方差數(shù)據(jù)分析軟件EddyPro（USA）處理后，剔除數(shù)據(jù)質(zhì)量標(biāo)記為不可用的數(shù)據(jù)[19?20]，使用R Package ‘REddyProc’（https://github.com/ bgctw/REddyProc）計(jì)算摩擦風(fēng)速閾值[10,21?22]，剔除低湍流條件（摩擦風(fēng)速小于0.22m·s?1）下的觀測數(shù)據(jù)[23]，將最終獲得的 0.5h 尺度凈生態(tài)系統(tǒng)碳交換量（NEE）數(shù)據(jù)作為基準(zhǔn)數(shù)據(jù)集。NEE 基準(zhǔn)數(shù)據(jù)集共有13200 條數(shù)據(jù)，數(shù)據(jù)缺失率為39%，其中白天（Rg＞20W·m?2）數(shù)據(jù)缺失率為32%，夜間（Rg≤20W·m?2）數(shù)據(jù)缺失率為48%。將小氣候觀測系統(tǒng)10min 氣象數(shù)據(jù)聚合為30min，數(shù)據(jù)缺失率為9%，小于2h 的缺失數(shù)據(jù)用線性內(nèi)插法插補(bǔ)，大于2h 的缺失數(shù)據(jù)采用平均晝夜變化法（MDV）進(jìn)行插補(bǔ)，選擇18m 位置處氣溫（Ta）、飽和水汽壓差（VPD）、17m 處總輻射（Rg）作為氣象基準(zhǔn)數(shù)據(jù)集，用于數(shù)據(jù)插補(bǔ)。

1.2 數(shù)據(jù)缺失情景

為評(píng)估不同插值方法對(duì)缺失片段長度的敏感性，隨機(jī)生成缺失片段重復(fù)次數(shù)不同，但缺失總數(shù)大致相同，以連續(xù)1、3、7、15 和31d 數(shù)據(jù)為缺失片段的5 類數(shù)據(jù)缺失情景（表1），各缺失片段分布均服從隨機(jī)分布，每類缺失約占基準(zhǔn)數(shù)據(jù)集的10%。為提高數(shù)據(jù)利用率，保證計(jì)算結(jié)果的準(zhǔn)確性和穩(wěn)定性，各類情景隨機(jī)重復(fù)10 次，并與基準(zhǔn)數(shù)據(jù)集相疊加，生成50 個(gè)缺失數(shù)據(jù)集。

表1 五類數(shù)據(jù)缺失情景Table 1 Description of five gap scenarios

1.3 數(shù)據(jù)插補(bǔ)方法

1.3.1 平均晝夜變化法

平均晝夜變化法（MDV）使用鄰近一段時(shí)間內(nèi)相同時(shí)段的觀測平均值來代替缺失值，分為獨(dú)立窗口法和滑動(dòng)窗口法[6]。獨(dú)立窗口法使用特定窗口內(nèi)任一時(shí)間點(diǎn)在該時(shí)刻所有有效觀測數(shù)據(jù)的平均值來代替缺失值（一般要求窗口內(nèi)至少有4 個(gè)有效觀測數(shù)據(jù)），滑動(dòng)窗口法則使用缺失數(shù)據(jù)周圍指定大小窗口內(nèi)的所有有效觀測數(shù)據(jù)的平均值來代替缺失值。窗口大小通常為4～15d。采用獨(dú)立窗口法，即

式中，h（1，…，48）為一天中每0.5h 的索引，i[1，…，interger（d/n）+1]為平均窗口的索引，n 為窗口大小，d 為一年的天數(shù)，k 為一中間變量，上劃線表示排除缺失數(shù)據(jù)后對(duì)該上滑線下的子集進(jìn)行算術(shù)平均，下劃線表示消除子集中缺失值后的索引。窗口大小選擇：（1）固定窗口大小為15d，在連續(xù)缺失長達(dá)31d 情景下擴(kuò)大為30d；（2）可變窗口，先以7d 為窗口插補(bǔ)數(shù)據(jù)，若還有數(shù)據(jù)缺失，則逐步擴(kuò)大窗口為14d、21d、28d…，直至完成全部缺失數(shù)據(jù)插補(bǔ)。分別使用 Excel 和 R Package ‘REddyProc’sEddyProc_sFillMDC 命令基于窗口（1）和（2）進(jìn)行缺失NEE 數(shù)據(jù)插補(bǔ)，分別記為MDV 和MDC。

1.3.2 查表法

查表法（LUT）基于6 個(gè)雙月或4 個(gè)季節(jié)時(shí)段，建立特定站點(diǎn)各種氣象條件下的NEE 索引表，根據(jù)缺失數(shù)據(jù)時(shí)段的氣象條件在NEE 索引表中查找相似環(huán)境下的NEE 來代替缺失數(shù)據(jù)[6]。通常選取光強(qiáng)和溫度作為索引因子，光強(qiáng)以100μmol·m?2·s?1的間隔從0 漸增至2200μmol·m?2·s?1，溫度以2℃的間隔從可能最低溫度到可能最高溫度，缺失NEE 數(shù)據(jù)用線性內(nèi)插法生成。使用 R Package ‘REddyProc’sEddyProc_sFillLUT 命令建立季節(jié)索引表進(jìn)行數(shù)據(jù)插補(bǔ)，記為LUT。

1.3.3 非線性回歸法

非線性回歸法（NLR）基于參數(shù)化非線性方程，通過建立一定時(shí)間內(nèi)有效NEE 觀測值與相關(guān)環(huán)境因子（如溫度和光強(qiáng)）的經(jīng)驗(yàn)方程，根據(jù)缺失時(shí)段的環(huán)境因子估算缺失的NEE[3,6,14]。通常，將日間與夜間的NEE 數(shù)據(jù)分開處理，使用呼吸方程來計(jì)算夜間NEE（等于夜間生態(tài)系統(tǒng)呼吸），使用光響應(yīng)方程結(jié)合日間生態(tài)系統(tǒng)呼吸來計(jì)算日間NEE；進(jìn)行回歸分析的時(shí)段無明確限定。使用R Package ‘Nonlinear Least Squares’ nls 命令擬合Lloyd & Taylor 方程和Michaelis-Menten 方程進(jìn)行數(shù)據(jù)插補(bǔ)，記為NLR。

呼吸方程，即Lloyd & Taylor 方程為[6,24]

式中，F(xiàn)RE,night為夜間生態(tài)系統(tǒng)呼吸（等于夜間的NEE，μmol·m?2·s?1，以CO2物質(zhì)的量計(jì)）；E0是常量，常設(shè)為309K；Tref是參考溫度（K），一般為298.16K；FRE,Tref為參考溫度下的生態(tài)系統(tǒng)呼吸（μmol·m?2·s?1，以CO2物質(zhì)的量計(jì)）；T0是生態(tài)系統(tǒng)呼吸為0 時(shí)的溫度（K）； T0為空氣或土壤溫度（K）。參數(shù)T0和FRE,Tref通過觀測數(shù)據(jù)回歸擬合得到。

光響應(yīng)方程，即Michaelis-Menten 方程為[6,25]

式中，QPPFD為光量子通量密度（μmol·m?2·s?1）；α '為生態(tài)系統(tǒng)量子效率（μmolCO2·μmol?1quanta）；FGPP,opt是最佳光照條件下的總初級(jí)生產(chǎn)力（μmol·m?2·s?1，以CO2物質(zhì)的量計(jì)）；FRE,day是日間的生態(tài)系統(tǒng)呼吸（μmol·m?2·s?1，以CO2物質(zhì)的量計(jì)）。參數(shù) α '和FGPP,opt通過觀測數(shù)據(jù)回歸擬合得到。

1.3.4 邊際分布采樣法邊際分布采樣法（MDS）是平均晝夜變化法和查表法的綜合使用，在氣溫（Ta）、總輻射（Rg）和飽和水汽壓差（VPD）觀測數(shù)據(jù)均可用時(shí)，在一定的時(shí)間窗口（缺失數(shù)據(jù)前后14～28d）內(nèi)，分別以2.5℃、50W·m?2和0.5kPa 的變異范圍約束Ta、Rg和VPD 進(jìn)行數(shù)據(jù)插補(bǔ)；在僅Rg 數(shù)據(jù)可用時(shí)，將插補(bǔ)的時(shí)間窗口縮小至前后14d；在氣象觀測要素全部缺失時(shí)，則使用平均晝夜變化法對(duì)缺失數(shù)據(jù)進(jìn)行插補(bǔ)；若仍有缺失則擴(kuò)大時(shí)間窗口重復(fù)上述步驟，直至完成對(duì)全部缺失數(shù)據(jù)的插補(bǔ)[10]。使用R Package‘REddyProc’ sEddyProc_sMDSGapFill 命令對(duì)缺失數(shù)據(jù)進(jìn)行插補(bǔ)，記為MDS。

1.3.5 人工神經(jīng)網(wǎng)絡(luò)法

人工神經(jīng)網(wǎng)絡(luò)（ANN）是基于計(jì)算機(jī)網(wǎng)絡(luò)還原人腦或生物神經(jīng)的網(wǎng)絡(luò)結(jié)構(gòu)和激勵(lì)行為的并行非線性計(jì)算系統(tǒng)[[7,26?28]，通過建立經(jīng)驗(yàn)非線性回歸模型進(jìn)行數(shù)據(jù)插補(bǔ)。目前，較常使用的人工神經(jīng)網(wǎng)絡(luò)是在監(jiān)督訓(xùn)練程序下能夠還原各環(huán)境因子間復(fù)雜關(guān)系的誤差反傳、信息前饋神經(jīng)網(wǎng)絡(luò)（Feed-forward backpropagation neural network，BP 網(wǎng)絡(luò)）。在氣象樣本數(shù)據(jù)信號(hào)（通常經(jīng)過標(biāo)準(zhǔn)化處理轉(zhuǎn)換為[0,1]）由輸入層進(jìn)入BP 神經(jīng)網(wǎng)絡(luò)各結(jié)點(diǎn)后，各輸入變量都被乘以分配給該結(jié)點(diǎn)的權(quán)重值，并經(jīng)轉(zhuǎn)換函數(shù)反復(fù)修正各層神經(jīng)元之間的連接權(quán)重，使誤差評(píng)價(jià)函數(shù)最優(yōu)，實(shí)現(xiàn)網(wǎng)絡(luò)實(shí)際輸出與期望輸出的誤差最小化。以氣溫（Ta）、飽和水汽壓差（VPD）、總輻射（Rg）作為輸入變量，以70%的數(shù)據(jù)作為訓(xùn)練集（training set）、15%的數(shù)據(jù)作為驗(yàn)證集（validation set）、15%的數(shù)據(jù)作為檢驗(yàn)集（testing set），使用Matlab Neural Net Fit 工具箱構(gòu)建具有10 個(gè)隱含神經(jīng)元、基于Bayesian Regularization 的BP 神經(jīng)網(wǎng)絡(luò)來進(jìn)行數(shù)據(jù)插補(bǔ)，記為ANN。

1.4 統(tǒng)計(jì)參數(shù)

使用決定系數(shù)（R2）、相對(duì)均方根誤差（RRMSE）和平均絕對(duì)誤差（MAE）評(píng)估各插補(bǔ)方法的性能。

式中， pi為各插補(bǔ)方法預(yù)測的 NEE（μmol·m?2·s?1），oi為實(shí)際觀測的NEE（μmol·m?2·s?1），為各插補(bǔ)方法預(yù)測NEE 的均值（μmol·m?2·s?1），為實(shí)際觀測NEE 的均值（μmol·m?2·s?1），均以CO2物質(zhì)的量計(jì)；N 為樣本數(shù)。

由于基準(zhǔn)數(shù)據(jù)集中晝夜缺失比率不同，在對(duì)插補(bǔ)方法進(jìn)行比較時(shí)，若將日間和夜間整合分析，白天和黑夜對(duì)統(tǒng)計(jì)指標(biāo)的貢獻(xiàn)權(quán)重的差異會(huì)影響最終的結(jié)果，帶來一定的偏差。因此，單獨(dú)計(jì)算和分析日間和夜間的統(tǒng)計(jì)參數(shù)。

使用 SPSS 25 對(duì)統(tǒng)計(jì)參數(shù)進(jìn)行方差分析（ANOVA）和多重比較。使用Origin 2021 作圖。

2 結(jié)果與分析

2.1 不同方法日間數(shù)據(jù)插補(bǔ)效果的比較

如圖1a 所示，在5 類缺失情景下，不同方法插補(bǔ)所得日間NEE 與實(shí)測NEE 間的R2均在0.5 以下；NLR（非線性回歸法）的R2最低，始終在0.2 以下；LUT（查表法）在連續(xù)1d 和3d 缺失情景下極差較?。O差在0.12 以下），在連續(xù)缺失達(dá)到7d 時(shí)極差明顯增大（極差始終大于0.18）。MDV（固定窗口平均晝夜變化法）和MDC（可變窗口平均晝夜變化法）的R2在連續(xù)1、3 和7d 缺失情景下波動(dòng)較小，結(jié)果相對(duì)穩(wěn)定，在連續(xù)15 和31d 缺失情景下變異系數(shù)顯著增大；MDS（邊際分布采樣法）和ANN（人工神經(jīng)網(wǎng)絡(luò)法）的R2在連續(xù)1、3、7 和15d 缺失情景下波動(dòng)較為穩(wěn)定（極差在0.2 左右），在連續(xù)31d 缺失情景下極差開始增大，分布更為離散，獲得結(jié)果的穩(wěn)定性趨于減弱。

圖1 五類缺失情景下不同方法插補(bǔ)所得日間NEE 數(shù)據(jù)集的統(tǒng)計(jì)參數(shù)（2017年3月1日?11月30日，Rg＞20W·m?2，0.5h 數(shù)據(jù)間隔）Fig.1 Statistical parameters of the daytime NEE datasets obtained by six gap?filling methods in five types of gap scenarios(Mar.1st-Nov.30th,2017,Rg＞20W· 0.5h data interval)

由表2可知，在連續(xù)1d 缺失情景下，ANN 方法的R2最高，NLR 方法的R2最低，與其它方法存在顯著差異（P＜0.05）；在連續(xù)3 和連續(xù)7d 缺失情景下，NLR 方法的R2仍顯著最低，ANN 方法的R2較高，但與MDV 方法差異不顯著（P＜0.05）；當(dāng)缺失達(dá)到連續(xù)15d 時(shí)，NLR 方法的R2最低，與其它方法差異顯著（P＜0.05），而其它方法間R2差異不顯著。隨著缺失片段長度的增加，除MDV 方法外，各方法的R2均呈下降趨勢，MDS 方法的R2在連續(xù)15d缺失與連續(xù)31d 缺失情景下的R2差異顯著（P＜0.05）；MDC 和NLR 方法的R2在連續(xù)7d 缺失與連續(xù)31d 缺失情景下的R2差異顯著（P＜0.05）；LUT和ANN 方法的R2在連續(xù)3d 缺失與連續(xù)31d 缺失情景下的R2差異顯著（P＜0.05）；MDV 方法的R2始終無顯著差異。

表2 五類缺失情景下不同方法插補(bǔ)所得日間NEE 數(shù)據(jù)集R2 多重比較表（Duncan 法）Table 2 The multiple comparison of R2 of the daytime NEE datasets obtained by six gap-filling methods in five types of gap scenarios (Duncan method)

如圖1b 所示，在連續(xù)缺失≤15d 時(shí)，LUT 方法插補(bǔ)所得日間NEE 與實(shí)測值NEE 間的RRMSE（相對(duì)均方根誤差）明顯低于其它方法，平均值在0.4 左右，極差相對(duì)較小，分布較為集中，但在連續(xù)31d缺失情景下，極差顯著增大，分布較為離散；MDV、MDC 和MDS 方法的RRMSE 在連續(xù)1、3 和7d 缺失情景下波動(dòng)較小，在連續(xù)15d 和連續(xù)31d 缺失情景下極差顯著增大；NLR 和ANN 方法的RRMSE 波動(dòng)范圍相對(duì)平穩(wěn)，在連續(xù)31d 缺失情景下，極差有增大趨勢，分布開始更為離散。

由表3可知，在連續(xù)1d 缺失情景下，LUT 方法的RRMSE 最低，NLR 方法的RRMSE 最高，與其它方法存在顯著差異（P＜0.05）；在連續(xù)3d 和連續(xù)7d 缺失情景下，LUT 方法的RRMSE 仍最低，NLR方法的RRMSE 較高，但與MDS 方法差異不顯著，與MDV、ANN 方法差異顯著（P＜0.05）；當(dāng)缺失達(dá)到連續(xù)15d 時(shí)，LUT 方法的RRMSE 較低，與其它方法有顯著差異（P＜0.05），而其它方法間RRMSE差異不顯著；當(dāng)缺失達(dá)到連續(xù)31d 時(shí)，各方法間RRMSE 無顯著差異。隨著缺失片段長度的增加，MDV 和MDS 方法的RRMSE 呈增大趨勢，連續(xù)1d缺失與連續(xù)31d 缺失情景下的RRMSE 差異顯著（P＜0.05）；其它方法的RRMSE 差異相對(duì)不顯著。

表3 五類缺失情景下不同方法插補(bǔ)所得日間NEE 數(shù)據(jù)集RRMSE 多重比較表（Duncan 法）Table 3 The multiple comparison of RRMSE of the daytime NEE datasets obtained by six gap-filling methods in five types of gap scenarios (Duncan method)

如圖1c 所示，在連續(xù)缺失＜15d 時(shí)，各插補(bǔ)方法所得日間NEE 與實(shí)測值NEE 間的MAE（平均絕對(duì)誤差）無明顯差異，分布較為集中；在連續(xù)缺失31d 情景下，MDV 方法的MAE 出現(xiàn)較多異常值，各方法之間的MAE 開始出現(xiàn)分化的趨勢。

2.2 不同方法夜間通量插補(bǔ)效果的比較

如圖2a 所示，與日間通量插補(bǔ)結(jié)果不同，在5類缺失情景下，夜間通量插補(bǔ)所得NEE 與實(shí)測NEE間的R2普遍較小，不同插補(bǔ)方法的R2始終在0.2 以下。LUT（查表法）的R2始終很小，分布集中且趨近于0；MDV（固定窗口平均晝夜變化法）的R2相對(duì)穩(wěn)定，隨連續(xù)缺失片段的增加無明顯變化；ANN（人工神經(jīng)網(wǎng)絡(luò)法）的R2均值相對(duì)較高，但始終存在較大波動(dòng)，穩(wěn)定性較差，隨著連續(xù)缺失片段的增加分布趨于集中；在連續(xù)缺失片段增加到7d 時(shí)，各方法的平均R2趨于相等，差異越來越??；在連續(xù)缺失片段增加到15d 時(shí)，MDC（可變窗口平均晝夜變化法）、MDS（邊際分布采樣法）和NLR（非線性相關(guān)法）的R2趨于穩(wěn)定。

由表4可知，在各缺失情景下，ANN 方法的R2較高，LUT 方法的R2較低，二者之間差異顯著（P＜0.05）；在連續(xù)缺失片段長度大于3d 的情景下，ANN 與MDV、NLR 方法的R2的差異不再顯著；在連續(xù)缺失片段達(dá)到15d 的情景下，ANN 與MDS 方法的R2出現(xiàn)顯著差異（P＜0.05）；而在連續(xù)缺失片段達(dá)到31d 時(shí)，ANN 與MDC 方法的R2出現(xiàn)顯著差異（P＜0.05）。隨著缺失片段長度的增加，MDC、MDS 和ANN 方法的R2呈下降趨勢，MDC 和MDS方法在連續(xù)7d 缺失與連續(xù)31d 缺失情景下的R2差異顯著（P＜0.05），ANN 方法在連續(xù)1d 和連續(xù)3d 缺失與連續(xù)15 和連續(xù)31d 缺失情景下的R2差異顯著（P＜0.05），MDV 與LUT 方法的R2始終無顯著差異。

如圖2b 所示，在各類缺失情景下，LUT 方法插補(bǔ)所得夜間NEE 與實(shí)測NEE 間的RRMSE（相對(duì)均方根誤差）明顯高于其它方法，平均值始終在0.9 以上，且極差相對(duì)較大，分布較為離散；NLR 和MDS方法的RRMSE 在連續(xù)1d 和連續(xù)3d 缺失情景下波動(dòng)較小，在連續(xù)缺失達(dá)到7d 時(shí)分布開始趨于離散；與MDV 方法不同，MDC 和ANN 方法的RRMSE 波動(dòng)相對(duì)平穩(wěn)，在各缺失情景下具有相似的極差。

由表5可知，在5 類缺失情景下，LUT 方法的RRMSE 最高，與其它方法存在顯著差異（P＜0.05）。除LUT 方法外，在連續(xù)缺失≤15d 的情景下，ANN方法的RRMSE 顯著低于MDC 方法（P＜0.05），而在連續(xù)缺失大于31d 的情景下，各方法的RRMSE 差異均不顯著。隨著缺失片段長度的增加，各方法的RRMSE 差異無顯著變化。

如圖2c 所示，夜間不同方法插補(bǔ)所得夜間NEE與實(shí)測NEE 間的MAE（平均絕對(duì)誤差）遠(yuǎn)低于日間；LUT 方法的MAE 始終為正，明顯高于其它方法，且分布較為離散，在連續(xù)缺失1d 情景下極差約為其它方法的2 倍；在連續(xù)缺失片段長度達(dá)到7d 時(shí)，除LUT方法外，其它方法的MAE 開始出現(xiàn)較大的波動(dòng)，分布趨于離散；在5 類缺失情景下，除LUT 方法外，其它方法的MAE 無明顯差異；與其它方法相比，MDV 方法的MAE 更趨近于0，NLR、MDS 和ANN方法的MAE 趨近于負(fù)值。

圖2 五類缺失情景下不同方法插補(bǔ)所得夜間NEE 數(shù)據(jù)集的統(tǒng)計(jì)參數(shù)（2017年3月1日?11月30日，Rg≤20W·m?2，0.5h 數(shù)據(jù)間隔）Fig.2 Statistical parameters of the nighttime NEE datasets obtained by six gap-filling methods in five types of gap scenarios(Mar.1st-Nov.30th,2017,Rg≤20W·m?2,0.5h data interval)

表4 五類缺失情景下不同方法插補(bǔ)所得夜間NEE 數(shù)據(jù)集R2 多重比較表（Duncan 法）Table 4 The multiple comparison of R2 of the nighttime NEE datasets obtained by six gap-filling methods in five types of gap scenarios(Duncan method)

表5 五類缺失情景下不同方法插補(bǔ)所得夜間NEE 數(shù)據(jù)集RRMSE 多重比較表（Duncan 法）Table 5 The multiple comparison of RRMSE of the nighttime NEE datasets obtained by six gap-filling methods in five types of gap scenarios (Duncan method)

2.3 不同方法典型晴天日變化數(shù)據(jù)插補(bǔ)效果的比較

為探究不同方法在0.5h 尺度上日變化的插值效果，以在5 類缺失情景中均涉及的生長穩(wěn)定的典型晴天2017年5月25日NEE 數(shù)據(jù)（數(shù)據(jù)缺失率為8%）為例，對(duì)比插補(bǔ)NEE 與實(shí)測NEE 的差異，分析各方法在不同缺失情景下對(duì)NEE日變化趨勢的還原效果（圖3）。

圖3 典型晴天下不同插補(bǔ)方法在五類缺失情景下NEE 的日變化（2017年5月25日，0.5h 數(shù)據(jù)間隔）Fig.3 Variation curve of the daily change of NEE filled by six gap-filling methods under five types of gap scenarios in a typical sunny day (May 25,2017,0.5h data interval )

由圖3可見，MDV（固定窗口平均晝夜變化法）在1d、3d 連續(xù)缺失情景下，對(duì)0.5h 尺度NEE 插值效果相對(duì)較好，但當(dāng)連續(xù)缺失達(dá)到7d 時(shí)，對(duì)夜間0.5h尺度NEE 開始出現(xiàn)較為明顯的低估，對(duì)日間0.5h 尺度NEE 開始出現(xiàn)較為明顯的高估，不能準(zhǔn)確預(yù)測日間CO2的吸收峰值。

MDC（可變窗口平均晝夜變化法）對(duì)0.5h 尺度NEE 插值效果相對(duì)最優(yōu)，在不同連續(xù)缺失情景下均能較好地還原NEE 的日變化。

LUT（查表法）對(duì)夜間0.5h 尺度NEE 存在明顯且穩(wěn)定的低估，在不同連續(xù)缺失情景下無明顯差異，對(duì)日間0.5h 尺度NEE 預(yù)測效果相對(duì)較好，在6：00?9：30 存在一定程度的低估，在10：00?16：00 存在一定程度的高估，使得日變化曲線在12：00?15：00時(shí)出現(xiàn)不明顯的“上凸”現(xiàn)象。

NLR（非線性回歸法）在不同缺失情景下，插補(bǔ)效果相對(duì)穩(wěn)定，但對(duì)夜間0.5h 尺度NEE 存在明顯低估，對(duì)日間0.5h 尺度NEE 存在明顯高估，“U”型曲線不明顯，對(duì)NEE日變化趨勢的還原效果最差。

與MDV 和MDC 方法類似，MDS（邊際分布采樣法）能夠在一定程度上還原夜間0.5h 尺度NEE 的波動(dòng)，但還原效果會(huì)隨著連續(xù)缺失片段長度的增加而降低；與MDV 和MDC 方法不同，MDS 方法對(duì)日間0.5h 尺度NEE 曲線的還原具有明顯的左偏現(xiàn)象，在9：00 左右會(huì)出現(xiàn)一個(gè)極小值。

與LUT 和NLR 方法類似，ANN（人工神經(jīng)網(wǎng)絡(luò)法）對(duì)夜間0.5h 尺度NEE 預(yù)測相對(duì)較差，存在明顯的低估，不能捕捉到夜間NEE 的波動(dòng)。同時(shí)，ANN 方法對(duì)日間NEE 存在高估，預(yù)測NEE 最小值出現(xiàn)的時(shí)間比實(shí)際觀測提前4h（由14：00 提前至10：00），與實(shí)際NEE 曲線相比在10：00?16：00存在明顯“上凸”現(xiàn)象。隨著缺失片段長度的增加，ANN 方法對(duì)夜間NEE 的低估和對(duì)日間NEE 的高估程度呈增加趨勢。

3 討論與結(jié)論

3.1 討論

3.1.1 日間和夜間NEE 插補(bǔ)效果的差異

在對(duì)缺失NEE 數(shù)據(jù)進(jìn)行插補(bǔ)時(shí)，日間和夜間插補(bǔ)效果存在顯著差異，日間不同方法插補(bǔ)所得NEE數(shù)據(jù)與實(shí)測數(shù)據(jù)的R2明顯高于夜間，RRMSE（相對(duì)均方根誤差）明顯低于夜間，插補(bǔ)效果明顯優(yōu)于夜間。

日間大氣層結(jié)不穩(wěn)定，熱量交換頻繁，對(duì)流作用較強(qiáng)，適合湍流交換，基本滿足利用渦度相關(guān)觀測系統(tǒng)進(jìn)行通量觀測的要求，經(jīng)過數(shù)據(jù)質(zhì)量控制后數(shù)據(jù)缺失較少。而夜間輻射冷卻導(dǎo)致大氣層結(jié)穩(wěn)定，對(duì)流較弱，抑制湍流混合，渦度相關(guān)法通量觀測結(jié)果較差，在剔除不滿足通量觀測條件的值后，往往會(huì)產(chǎn)生較多的缺失。同時(shí)，夜間經(jīng)常發(fā)生平流、泄流效應(yīng)，垂直方向上湍流運(yùn)動(dòng)傾向于向高頻移動(dòng)，以小渦運(yùn)動(dòng)占優(yōu)勢，開路式渦度相關(guān)系統(tǒng)傳感器的分離等會(huì)造成觀測儀器響應(yīng)的不足，進(jìn)而引起對(duì)夜間通量的低估，引入選擇性系統(tǒng)誤差[1]。在對(duì)NEE數(shù)據(jù)進(jìn)行插補(bǔ)時(shí)，由于夜間可用樣本數(shù)遠(yuǎn)低于日間，再加上夜間渦度相關(guān)法更傾向于低估NEE 的特性，使得夜間NEE 插補(bǔ)效果遠(yuǎn)低于日間。

3.1.2 不同插補(bǔ)方法NEE 插補(bǔ)效果的差異

不同插補(bǔ)方法的通量插補(bǔ)效果存在差異。不管是日間還是夜間，ANN（人工神經(jīng)網(wǎng)絡(luò)法）往往總能取得相對(duì)較好的插補(bǔ)結(jié)果，而NLR（非線性相關(guān)法）則相對(duì)表現(xiàn)較差；LUT（查表法）在日間的表現(xiàn)明顯優(yōu)于夜間；MDV（固定窗口平均晝夜變化法）、MDC（可變窗口平均晝夜變化法）和MDS（邊際分布采樣法）之間差異不顯著。

與Moffat 等[29?30]的研究結(jié)果相比，本站點(diǎn)各方法對(duì)NEE 數(shù)據(jù)插補(bǔ)結(jié)果的R2相對(duì)較低。這可能與站點(diǎn)因素有關(guān)，本研究所用通量觀測數(shù)據(jù)來源于山地丘陵林區(qū)，下墊面相對(duì)較為復(fù)雜，與農(nóng)田和草地相比，對(duì)通量觀測和數(shù)據(jù)插補(bǔ)的要求更高。

與Ooba 等[12,29,31]的結(jié)果類似，ANN 方法在連續(xù)缺失小于7d 時(shí)，數(shù)據(jù)插補(bǔ)結(jié)果的R2較穩(wěn)定，高于其它方法，RRMSE（相對(duì)均方根誤差）也相對(duì)較低，表現(xiàn)出較好的插值效果。ANN 方法對(duì)0.5h 尺度NEE的日變化還原較差，傾向于高估中午時(shí)分的NEE，一方面源于統(tǒng)計(jì)平均效應(yīng)，另一方面也與輸入?yún)?shù)的數(shù)量有關(guān)。除將總輻射、氣溫、飽和水汽壓作為輸入?yún)?shù)外，進(jìn)一步增加土壤溫濕度、風(fēng)速等氣象因子，并引入季節(jié)、生長期等非連續(xù)變量，可以進(jìn)一步提高ANN 方法的精度。

NLR 方法在夜間通量插補(bǔ)中表現(xiàn)出較高的R2和較低的RRMSE，但在日間通量插補(bǔ)中則相反，與Moffat 等[29,32]NLR 方法在日間仍表現(xiàn)相對(duì)較好的結(jié)果存在差異。這可能與環(huán)境響應(yīng)方程的選擇和擬合方程時(shí)所選擇的時(shí)段有關(guān)，不同物候期內(nèi)栓皮櫟人工林可能具有不同的光響應(yīng)曲線，使用單月或更為精確的物候分期分段擬合環(huán)境響應(yīng)方程或先將觀測數(shù)據(jù)按溫度分組后再利用環(huán)境響應(yīng)方程進(jìn)行回歸模擬可能取得更好的插補(bǔ)效果。

LUT 方法在日間通量插補(bǔ)中表現(xiàn)出較高的R2和較低的RRMSE，在夜間通量插補(bǔ)中則相反。夜間可用實(shí)測NEE 的數(shù)據(jù)量遠(yuǎn)少于日間，絕對(duì)值也較日間低，在以季節(jié)為單位建立索引表時(shí)，較多的缺失會(huì)引起LUT 方法傾向于低估夜間NEE，帶來難以避免的系統(tǒng)誤差。

MDV、MDC 和MDS 方法具有一定的相似性，3 種方法均使用缺失數(shù)據(jù)周圍一段時(shí)間內(nèi)相關(guān)數(shù)據(jù)的均值來代替缺失值。但MDV 和MDC 方法并未考慮氣象因素，MDS 方法則以氣溫、總輻射和飽和水汽壓差梯度為限制條件；MDV 方法的窗口大小是固定的，MDC 和MDS 方法的窗口則是動(dòng)態(tài)變化的。與Falge 等[6,29?30,32]的研究結(jié)果一致，在連續(xù)缺失15d以下情景時(shí)，MDV 方法技術(shù)具有較好的性能。但在統(tǒng)計(jì)平均效應(yīng)的作用下，當(dāng)連續(xù)缺失大于15d 時(shí)，MDV 方法會(huì)對(duì)日間NEE 產(chǎn)生高估。在數(shù)據(jù)連續(xù)缺失≤15d 時(shí)，MDV、MDC 和MDS 方法的R2和RRMSE 無顯著差異，都能較好地還原0.5h 尺度NEE的日變化，但MDC 方法還原性能更好。

3.1.3 不同插補(bǔ)方法的適用性

不同NEE 插補(bǔ)方法的插補(bǔ)效果與數(shù)據(jù)缺失的持續(xù)時(shí)間有關(guān)。隨著數(shù)據(jù)缺失時(shí)間的延長，各插補(bǔ)方法所得日間NEE 結(jié)果的穩(wěn)定性一般會(huì)越來越差。各通量插補(bǔ)方法一般有其適用范圍，即在一定限度的連續(xù)缺失時(shí)間內(nèi)具有較好的插補(bǔ)效果，超過此限度后插補(bǔ)效果會(huì)顯著降低。ANN（人工神經(jīng)網(wǎng)絡(luò)法）適用范圍較廣，即使是31d 的連續(xù)缺失，往往也能取得較好的插補(bǔ)效果，但在連續(xù)缺失少于7d 時(shí)能獲得更好的結(jié)果；LUT（查表法）、MDV（固定窗口平均晝夜變化法）、MDC（可變窗口平均晝夜變化法）和MDS（邊際分布采樣法）的適用限度均在15d 以內(nèi)；NLR（非線性相關(guān)法）表現(xiàn)相對(duì)較差，更適合7d 內(nèi)的缺失插補(bǔ)。

Moffat 等[29,31]發(fā)現(xiàn)，ANN 方法在連續(xù)缺失少于7d 時(shí)能獲得非常好的結(jié)果，在連續(xù)12d 缺失時(shí)，仍保持最優(yōu)的插補(bǔ)效果。即使缺失時(shí)間延長至31d，ANN方法仍能取得較優(yōu)的插補(bǔ)結(jié)果。與Falge 等[6,29?30,32]的研究結(jié)果一致，在缺失3～7d，甚至是12d 左右時(shí)，MDV 方法仍保持較好的性能，但在缺失時(shí)間繼續(xù)延長時(shí)，其穩(wěn)定性會(huì)大幅下降。Moffat 等[29?30]發(fā)現(xiàn)在連續(xù)缺失少于12d 時(shí)，MDS 方法具有較好的插補(bǔ)性能，與本研究MDS 方法在連續(xù)缺失達(dá)到15d 時(shí)插補(bǔ)性能出現(xiàn)下降且穩(wěn)定性開始變差略有差異。這可能與站點(diǎn)特性有關(guān)，Moffat 等[29]分析對(duì)象為歐洲典型森林，而Du 等[30]關(guān)注的則是退化草地和玉米農(nóng)田，與本研究所探討的暖溫帶落葉闊葉林存在較大差異，站點(diǎn)的差異影響了不同插補(bǔ)方法的適用性及穩(wěn)定性。

在對(duì)缺失NEE 數(shù)據(jù)進(jìn)行插補(bǔ)時(shí)，當(dāng)NEE 數(shù)據(jù)缺失＜15d 且氣象數(shù)據(jù)不可用或缺失嚴(yán)重時(shí)，可以使用MDV 或MDC 方法；而當(dāng)NEE 數(shù)據(jù)缺失＜15d，且氣象數(shù)據(jù)可用時(shí)，則優(yōu)先使用LUT、ANN 和MDS方法；在NEE 數(shù)據(jù)缺失≥15d 時(shí)，多次利用ANN 方法進(jìn)行插值并取均值可能是比較好的選擇。在關(guān)注NEE日變化趨勢時(shí)，可優(yōu)先使用MDC 方法。在數(shù)據(jù)缺失比較嚴(yán)重，可用數(shù)據(jù)量較少時(shí)，NLR 方法會(huì)有較大的誤差。

除站點(diǎn)因素外，不同插補(bǔ)方法選取的時(shí)間步長和窗口大小的差異等也會(huì)影響缺失通量數(shù)據(jù)插補(bǔ)效果，進(jìn)而影響各插補(bǔ)方法的適用性，本研究僅考慮了單一站點(diǎn)一年（除冬季）的通量數(shù)據(jù)，在構(gòu)建人工缺失集時(shí)忽略了實(shí)際缺失的分布，所選插補(bǔ)方法在進(jìn)行插補(bǔ)時(shí)所選用的時(shí)間步長和窗口大小也不盡相同，其結(jié)果可能并不適用于所有站點(diǎn)，但可為其它站點(diǎn)數(shù)據(jù)插補(bǔ)方法的選擇提供參考。同時(shí)，部分通量數(shù)據(jù)的缺失源于降水、露水等異常天氣的影響，通過上述方法插補(bǔ)所得的通量數(shù)據(jù)可能與實(shí)際通量有較大差異（顯著高估），尤其是不考慮氣象因素的MDV 和MDC 方法，要準(zhǔn)確估計(jì)這部分通量，還需與閉路式渦度相關(guān)觀測系統(tǒng)相結(jié)合，進(jìn)行相應(yīng)的數(shù)據(jù)校正研究。

3.2 結(jié)論

（1）在對(duì)缺失NEE 數(shù)據(jù)進(jìn)行數(shù)據(jù)插補(bǔ)時(shí)，由于湍流穩(wěn)態(tài)因素和開路式渦度相關(guān)觀測系統(tǒng)傳感器分離等的影響，日間插補(bǔ)效果顯著優(yōu)于夜間。由于插值策略的差異，不同插補(bǔ)方法的插補(bǔ)效果存在差異。ANN（人工神經(jīng)網(wǎng)絡(luò)法）插值效果總體較好，而NLR（非線性相關(guān)法）則相對(duì)表現(xiàn)較差；LUT（查表法）在日間的表現(xiàn)明顯優(yōu)于夜間，對(duì)夜間NEE 存在低估現(xiàn)象；MDV（固定窗口平均晝夜變化法）、MDC（可變窗口平均晝夜變化法）和MDS（邊際分布采樣法）之間差異不顯著，但MDC 方法對(duì)半小時(shí)通量變化的還原性能更好。

（2）不同NEE 插補(bǔ)方法的插補(bǔ)效果與連續(xù)數(shù)據(jù)缺失的持續(xù)時(shí)間有關(guān)。隨著連續(xù)缺失時(shí)間的延長，各插補(bǔ)方法所得結(jié)果的穩(wěn)定性一般會(huì)越來越差。NLR 方法適用于氣象數(shù)據(jù)完備、NEE 數(shù)據(jù)連續(xù)缺失少于7d 的情景；MDV 或MDC 方法適用于氣象數(shù)據(jù)不可用或缺失嚴(yán)重、NEE 數(shù)據(jù)連續(xù)缺失少于15d 的情景，優(yōu)先選用MDC 方法；LUT 和MDS 方法則適用于氣象數(shù)據(jù)缺失較少、NEE 數(shù)據(jù)連續(xù)缺失少于15d的情景；ANN 方法適用性相對(duì)較廣，可用于氣象數(shù)據(jù)缺失較少、NEE 數(shù)據(jù)連續(xù)缺失長達(dá)31d 的情景。