北京325米氣象塔上CO2梯度觀測(cè)數(shù)據(jù)質(zhì)量控制與評(píng)價(jià)

劉郁玨　胡非　程雪玲　馮永芳中國(guó)科學(xué)院大氣物理研究所，北京0009北京市房山區(qū)氣象局，北京0488

北京325米氣象塔上CO2梯度觀測(cè)數(shù)據(jù)質(zhì)量控制與評(píng)價(jià)

劉郁玨1, 2胡非1程雪玲1馮永芳2
1中國(guó)科學(xué)院大氣物理研究所，北京100029
2北京市房山區(qū)氣象局，北京102488

劉郁玨，胡非，程雪玲，等. 2016. 北京325米氣象塔上CO2梯度觀測(cè)數(shù)據(jù)質(zhì)量控制與評(píng)價(jià) [J]. 大氣科學(xué), 40 (2): 390?400.Liu Yujue, Hu Fei, Cheng Xueling, et al. 2016. Data processing and quality assessment of the eddy covariance system of the 325-meter meteorology tower in Beijing [J]. Chinese Journal of Atmospheric Sciences (in Chinese), 40 (2): 390?400, doi:10.3878/j.issn.1006-9895.1504.14325.

近些年渦度相關(guān)系統(tǒng)在城市通量研究中得到了廣泛應(yīng)用，因城市下墊面的特殊性和復(fù)雜性，以及系統(tǒng)觀測(cè)原理和儀器精度存在著局限性，使得結(jié)果存在10%以上偏差，故必須對(duì)原始觀測(cè)數(shù)據(jù)進(jìn)行質(zhì)量控制和評(píng)價(jià)。本文對(duì)北京325 m氣象塔七層高度上的CO2（二氧化碳）、能量等通量進(jìn)行了長(zhǎng)期觀測(cè)，研究了渦度相關(guān)技術(shù)在城市環(huán)境通量觀測(cè)中的適用性，設(shè)計(jì)出了一套適合于復(fù)雜城市下墊面上的渦度相關(guān)系統(tǒng)資料前處理和質(zhì)量控制方案，并對(duì)比了不同修正方案效果。分析結(jié)果表明，城市下墊面環(huán)境中計(jì)算CO2通量的最優(yōu)周期為30 min；二次坐標(biāo)旋轉(zhuǎn)法優(yōu)于平面擬合法；頻率響應(yīng)修正后的感熱通量提高5.21%，潛熱通量和CO2通量均提高9.42%。空氣密度脈動(dòng)修正Liu法（Liu，2005）優(yōu)于WPL法（Webb et al.，1980）；湍流譜在慣性副區(qū)滿足－2/3次方定律，協(xié)方差譜滿足－4/3次方定律。經(jīng)過(guò)這套前處理和質(zhì)量控制方案，原始數(shù)據(jù)中有79%能夠用于基礎(chǔ)研究，該質(zhì)量控制與評(píng)價(jià)體系可為復(fù)雜城市下墊面通量研究提供參考。

Founded byStrategic Priority Research Program—Climate Change: Carbon Budget and Related Issues of the Chinese Academy of Sciences (Grant XDA05040301)

1　引言

國(guó)際環(huán)境發(fā)展研究所相關(guān)研究表明，城市化和氣候變化是全球面臨的最重要的兩個(gè)問(wèn)題（Hoornweg et al.，2011），認(rèn)識(shí)城市碳通量的排放情況對(duì)評(píng)估未來(lái)溫室氣體濃度變化和幫助減輕氣候效應(yīng)具有重要意義。渦度相關(guān)法（EC法）作為可以直接測(cè)量湍流通量的方法，已廣泛應(yīng)用于生態(tài)系統(tǒng)CO2（二氧化碳）、能量通量觀測(cè)中，但在城市通量觀測(cè)中的應(yīng)用還為數(shù)不多。不同于森林、草地、農(nóng)田等均一生態(tài)系統(tǒng)下墊面，城市下墊面具有其特殊的復(fù)雜性，導(dǎo)致通量源、匯分布不均，使得直接測(cè)量和定量描述變得非常困難。近些年，眾多學(xué)者在城市密集區(qū)展開(kāi)了通量觀測(cè)實(shí)驗(yàn)（Grimmond et al.，2002；Vogt et al.，2005；Coutts et al.，2007；Matese et al.，2009；Crawford et al.，2011）。這些實(shí)驗(yàn)主要集中在發(fā)達(dá)國(guó)家，而在發(fā)展中國(guó)家開(kāi)展的城市通量觀測(cè)實(shí)驗(yàn)非常稀少（Burri et al., 2009; Song and Wang, 2012; Liu et al., 2012; 劉郁玨等，2014）。

城市湍流通量研究起步雖然較晚，但發(fā)展卻非常迅速。2012年，為滿足中國(guó)科學(xué)院戰(zhàn)略先導(dǎo)科技專項(xiàng)“應(yīng)對(duì)氣候變化的碳收支認(rèn)證及相關(guān)問(wèn)題”課題，于北京325 m氣象塔上新架設(shè)了七層渦度通量觀測(cè)系統(tǒng)，并開(kāi)展了觀測(cè)時(shí)長(zhǎng)達(dá)5年的觀測(cè)實(shí)驗(yàn)，該實(shí)驗(yàn)是目前國(guó)際上少有的位于發(fā)展中國(guó)家超大城市群內(nèi)的全天候、多梯度觀測(cè)實(shí)驗(yàn)。北京325 m氣象塔觀測(cè)設(shè)備架設(shè)在七層不同高度上，其觀測(cè)值不僅受到城市觀測(cè)的環(huán)境及氣象條件的影響，也受觀測(cè)點(diǎn)周圍的“足跡”即源區(qū)分布情況的影響（劉郁玨，2014），所以必須對(duì)渦度相關(guān)系統(tǒng)的觀測(cè)結(jié)果進(jìn)行質(zhì)量控制。這不僅要考慮儀器、傳感器的測(cè)量誤差，還需要考慮渦度相關(guān)方法建立的理論假設(shè)的滿足程度（Grimmond et al.，2002），故不同的處理方法會(huì)帶來(lái)分析結(jié)果10%或更大的差異（王介民，2007；王詠薇等，2013）。Velasco et al.（2005）建議城市CO2觀測(cè)系統(tǒng)傳感器位于慣性副層，以避免其結(jié)構(gòu)對(duì)流場(chǎng)的改變帶來(lái)的附加通量貢獻(xiàn)。Kotthaus and Grimmond（2012）利用微尺度技術(shù)將微尺度（10～100 m）排放源從局地尺度（100～10000 m）中辨別出來(lái)，但目前并沒(méi)有一套完備通用的資料處理控制方案。劉樹(shù)華等（2005a，2005b，2005c）利用EBEX-2000（International Energy Balance Experiment，2000）的湍流觀測(cè)資料，計(jì)算了不同穩(wěn)定度條件下的湍流宏觀量特征、湍流能量和熱量耗散率以及湍流結(jié)構(gòu)參數(shù)特征，并與其他湍流實(shí)驗(yàn)得到的結(jié)果進(jìn)行了比較。宋濤等（2014）設(shè)計(jì)了不同城市環(huán)境的熱通量數(shù)據(jù)處理方案，對(duì)城市適用性進(jìn)行了初步探討。以上研究為本文進(jìn)行七層渦度相關(guān)系統(tǒng)CO2通量觀測(cè)資料的前處理和質(zhì)量控制方案設(shè)計(jì)奠定了基礎(chǔ)。

2　研究區(qū)域簡(jiǎn)介及觀測(cè)實(shí)驗(yàn)概況

2.1區(qū)域簡(jiǎn)介

中國(guó)科學(xué)院大氣物理研究所氣象塔建成于1979年8月，氣象塔位于（39°58′N，116°22′E），海拔高度49 m，塔高325 m，塔上設(shè)有15層觀測(cè)平臺(tái)，并裝有風(fēng)速儀、風(fēng)向儀和溫度儀。為減小塔體阻擋氣流對(duì)測(cè)風(fēng)儀造成的影響，分別在西北和東南兩個(gè)盛行風(fēng)向各裝一臺(tái)測(cè)風(fēng)儀。其主要功能為首都空氣污染以及大氣邊界層、大氣湍流擴(kuò)散等提供研究服務(wù)。目前氣象塔半徑2 km內(nèi)地勢(shì)平坦，地表覆蓋以建筑為主，塔的北面和南面為密集居民區(qū)，覆蓋率為0.65；其次為植被和道路。距離氣象塔西部約0.3 km處為小型公園，由樹(shù)木、草坪和河流組成。塔的東面為京藏高速，塔的北面為北辰路，每天有大量車輛行駛，覆蓋率分別為0.21、0.14 （Miao et al., 2012），圖略。

2.2實(shí)驗(yàn)概況

為滿足中國(guó)科學(xué)院戰(zhàn)略先導(dǎo)科技專項(xiàng)“應(yīng)對(duì)氣候變化的碳收支認(rèn)證及相關(guān)問(wèn)題”課題（簡(jiǎn)稱碳專項(xiàng)）中“碳衛(wèi)星驗(yàn)證系統(tǒng)與綜合觀測(cè)”的第一子課題“高塔驗(yàn)證系統(tǒng)和綜合觀測(cè)”（簡(jiǎn)稱高塔驗(yàn)證）的需要，2012年于氣象塔上新架設(shè)了七層渦動(dòng)通量觀測(cè)系統(tǒng)（分別位于8、16、47、80、140、220和280 m處），如圖1所示。系統(tǒng)主要由開(kāi)路快速響應(yīng)紅外水汽—二氧化碳分析儀（Infrared Gas Analyzer，IRGA；Model Li-7500，Li-Cor Inc，USA）、三維超聲風(fēng)速儀（Wind Master，Gill，USA）和數(shù)據(jù)采集器（Model CR5000，Campbell，USA）組成。儀器設(shè)定的基本參數(shù)如表1所示，表1中u、v、w為三維方向的風(fēng)速，Ts為氣溫，p為氣壓，ρv為水汽密度，ρc為二氧化碳密度。2011年8月底，碳專項(xiàng)高塔驗(yàn)證課題確定了觀測(cè)方案，完成了觀測(cè)儀器申請(qǐng)及招標(biāo)工作。2012年7月1日，七層數(shù)據(jù)已經(jīng)正式完備。本文中所使用的數(shù)據(jù)均來(lái)自該觀測(cè)系統(tǒng)2012年12月1日至2014年2月14日間各高度層的觀測(cè)數(shù)據(jù)。

圖1　七層渦度通量觀測(cè)系統(tǒng)示意圖Fig. 1　The EC (Eddy Covariance) system of the 325-m meteorology tower in Beijing

表1　儀器參數(shù)表Table 1　Parameter table of instruments

3　數(shù)據(jù)質(zhì)量控制與評(píng)價(jià)方案設(shè)計(jì)

3.1方案設(shè)計(jì)

對(duì)觀測(cè)點(diǎn)的通量數(shù)據(jù)進(jìn)行質(zhì)量控制與評(píng)價(jià)是獲得可信的地—?dú)饨粨QCO2通量的必要步驟，也是后續(xù)分析過(guò)程的前提。然而對(duì)近地層通量數(shù)據(jù)進(jìn)行評(píng)價(jià)存在困難，目前還沒(méi)有通用的處理方案，且不同數(shù)據(jù)處理軟件采用的修正方法也略有不同，需研究者根據(jù)自身需求、實(shí)驗(yàn)數(shù)據(jù)的特點(diǎn)進(jìn)行反復(fù)的校正和評(píng)價(jià)分析，從而得到最適合的校正方案。本方案的處理步驟如圖2所示。

3.2輔助軟件選擇

目前國(guó)際上已有一些較為實(shí)用的渦度相關(guān)儀器數(shù)據(jù)的質(zhì)量控制軟件（徐自為，2008），如愛(ài)丁堡大學(xué)的EdiRe（Eddy Reprocessing）軟件、德國(guó)拜羅伊特TK2/TK3（Turbulence Knight 2/3）、荷蘭瓦赫寧根大學(xué)的ECPACK軟件、Campbell的自帶處理軟件等。其中EdiRe功能相對(duì)更加豐富，適合科研人員靈活組合運(yùn)用，本文推薦使用這款免費(fèi)處理軟件來(lái)處理地表微氣象通量數(shù)據(jù)，其相關(guān)教程及幫助文件參見(jiàn)（http://www.geos.ed.ac.uk/abs/ researchh/micromet/ EdiRe [2014-11-10]）。

4　數(shù)據(jù)質(zhì)量控制

4.1確定平均周期

城市碳通量的計(jì)算方法與自然生態(tài)系統(tǒng)相似，CO2通量（Fc）均通過(guò)計(jì)算垂直風(fēng)速脈動(dòng)w'（單位：m s?1）和CO2濃度脈動(dòng)c'（單位：μmol m?3）乘積在平均周期內(nèi)的平均值得到，如式（1）所示：

基于渦度相關(guān)技術(shù)觀測(cè)地氣間的CO2、H2O通量，需要考慮通量變化特征和微氣象相關(guān)原則來(lái)確定適宜的數(shù)據(jù)平均周期，其原則是盡可能包含各種頻率的湍流成分（脈動(dòng)量），數(shù)據(jù)越長(zhǎng)所包含的采樣周期內(nèi)的長(zhǎng)周期脈動(dòng)量越多，但數(shù)據(jù)平均長(zhǎng)度太長(zhǎng)則會(huì)削弱通量明顯的小時(shí)變化及日變化等時(shí)間尺度上的特征，因此平均周期的選取需考慮城市下墊面特征、人為活動(dòng)、植被光合作用等對(duì)CO2通量影響的時(shí)空特征。

在計(jì)算通量值之前，首先需分析低頻成分對(duì)北京城市觀測(cè)結(jié)果的影響，這里用到Ogive函數(shù)（Foken and Wichura, 1996），它是一個(gè)協(xié)譜從高頻至低頻的累積積分函數(shù)，標(biāo)量協(xié)譜定義為：FA、FB，分別為同時(shí)觀測(cè)的標(biāo)量A和B的離散傅里葉變化，其中FA= FAr+ FAi，F(xiàn)B= FBr+ FBi，下標(biāo)r、i代表離散傅氏變化實(shí)部、虛部。協(xié)譜表示為Co=FArFBr+ FAiFBi。Ogive函數(shù)表示如下：

式中w為垂直風(fēng)速，c為某一標(biāo)量，f_low為可分辨的最低頻率（2T?1，T為所選時(shí)間序列長(zhǎng)度），f_high為尼奎斯特（Nyquist）頻率，大小為采樣頻率一半，Cowc(t) 為垂直風(fēng)速w和標(biāo)量c密度的協(xié)譜。Ogive函數(shù)的曲線由高頻向低頻方向呈現(xiàn)漸進(jìn)的形式，在某一低頻頻率時(shí)則逐漸收斂為一定值，該頻率即為計(jì)算通量時(shí)所需要的平均時(shí)間。

圖2　渦度相關(guān)通量數(shù)據(jù)質(zhì)量處理與評(píng)價(jià)步驟圖Fig. 2　Flow chart of the data processing and quality assessment of the EC system

圖3a為2013年7月1日06:00～18:00（協(xié)調(diào)世界時(shí)，下同），47 m高度處每2 h平均周期后的CO2通量的Ogive函數(shù)曲線圖，圖3b為2013年7 月1日10:00～12:00，七層高度處2 h平均周期的Ogive函數(shù)曲線圖?？v向?qū)嵕€從左至右，分別對(duì)應(yīng)于15 min、30 min、60 min、120 min。從圖中不難發(fā)現(xiàn)，無(wú)論是同一層不同時(shí)段還是不同高度同一時(shí)段，CO2通量每120 min連續(xù)采樣的Ogive區(qū)間在15 min處，部分時(shí)間段通量的函數(shù)已經(jīng)收斂，余部分仍有上升趨勢(shì)；而在30 min的平均周期則基本己經(jīng)全部收斂，變化幅度很小。15 min到120 min的平均周期都是可以接受的，基本能獲取全部的低頻信息，沒(méi)有明顯的低頻信號(hào)缺失，這也證明了理論條件下計(jì)算的經(jīng)驗(yàn)平均時(shí)間30 min對(duì)于城市下墊面也是適用的。

圖3　基于2 h平均周期的CO2通量Ogive函數(shù)圖：（a）47 m高度，取自2013年7月1日06:00～18:00；(b) 七層高度，取自2013年7月1日10:00～12:00；從左至右縱向?qū)嵕€分別代表15 min、30 min、60 min、120 min的周期Fig. 3　Ogive function of CO2flux based on a 2-h average cycle: (a) 0600–1800 UTC 1 July 2013, at 47 m; (b) 1000–1200 UTC I July 2013, at 7 levels. Theblack lines from left to right represent periods 15 min, 30 min, 60 min and 120 min, respectively

4.2資料前處理

4.2.1野點(diǎn)剔除

原始的10 Hz采樣數(shù)據(jù)為觀測(cè)得到的直接數(shù)據(jù)，其中渦動(dòng)相關(guān)觀測(cè)系統(tǒng)除了觀測(cè)到有效信息外，也包含了噪聲信號(hào)，因此必須對(duì)原始數(shù)據(jù)進(jìn)行整理，對(duì)噪聲信號(hào)進(jìn)行剔除。使用的方法如下：（1）利用方差分析原理采用4倍標(biāo)準(zhǔn)差為閉值（徐自為，2008），對(duì)湍流數(shù)據(jù)進(jìn)行檢驗(yàn)并剔除噪聲；（2）利用觀測(cè)站點(diǎn)平均場(chǎng)資料，對(duì)降水前后一小時(shí)及降水期間的數(shù)據(jù)進(jìn)行剔除，剔除由雪、塵粒等環(huán)境因子對(duì)傳感器聲光程的干擾、瞬間斷電等因素造成的尖點(diǎn)假峰值，并利用超聲風(fēng)速儀Gill和二氧化碳水汽分析儀Licor7500的狀態(tài)異常標(biāo)志diag≠0，剔除電子電路、電纜電源不穩(wěn)定原因造成的連續(xù)假峰值；（3）由于冬季寒冷容易造成儀器結(jié)冰，儀器結(jié)冰后會(huì)直接影響觀測(cè)效果。Li-7500儀器觀測(cè)光路被若被冰晶所阻擋，儀器探測(cè)值失真，使得光路分析數(shù)據(jù)值A(chǔ)GC（Automatic Gain Control）明顯偏大。故可用AGC值來(lái)檢驗(yàn)儀器是否存在結(jié)冰，當(dāng)AGC值大于90時(shí)可認(rèn)為儀器結(jié)冰嚴(yán)重。

圖4　兩種坐標(biāo)旋轉(zhuǎn)修正后（a）摩擦速度、（b）CO2通量、（c）感熱通量和（d）潛熱通量比較圖Fig. 4　Comparison of (a) friction velocity, (b) CO2flux, (c) sensible heat flux, and (d) latent heat flux after DR (Double coordinate Rotation) and PF (Plane Fitting) coordinate rotation corrections

4.2.2坐標(biāo)旋轉(zhuǎn)

在湍流觀測(cè)過(guò)程中，實(shí)際觀測(cè)條件絕不是理想的。地形傾斜或起伏、儀器不是絕對(duì)水平固定等問(wèn)題會(huì)直接對(duì)垂直風(fēng)速觀測(cè)產(chǎn)生影響。尤其在城市這種極為復(fù)雜的下墊面上，絕對(duì)水平很難保證，目前常用的方法有二次坐標(biāo)旋轉(zhuǎn)法DR（Double coordinate Rotation）和平面擬合法PF（Plane Fitting）（Finnigan，2004），具體的操作可以使用EdiRe程序進(jìn)行編輯計(jì)算。為了比較這兩種方法在城市下墊面的適用性，取七層高度上，每個(gè)月的1日與15日2天的半小時(shí)平均值（摩擦風(fēng)速、CO2通量、感熱通量與潛熱通量）進(jìn)行兩種方法的修正，修正結(jié)果如圖4所示。

通過(guò)二次旋轉(zhuǎn)修正的數(shù)據(jù)相對(duì)于平面擬合修正后略小，兩種修正方法對(duì)感熱通量的修正最為接近，差異不足2%；而潛熱通量差異最大，平面擬合后比二次旋轉(zhuǎn)后減小4%；摩擦速度和CO2通量的修正差異介于二者之間。

圖5給出了不同風(fēng)向上，兩種方法修正后的湍流統(tǒng)計(jì)量相關(guān)性隨風(fēng)向變化規(guī)律。其相關(guān)系數(shù)隨風(fēng)向變化趨勢(shì)較為一致，相關(guān)系數(shù)達(dá)到0.97。在0～120°和210°～360°的風(fēng)向范圍內(nèi)相關(guān)系數(shù)最高，主要由于該風(fēng)向上，風(fēng)向與地形變換較為一致；另一方面該風(fēng)向上的數(shù)據(jù)較多，能夠保證平面擬合順利進(jìn)行。在120°～180°風(fēng)向范圍里兩種方法的修正結(jié)果差異較大，相關(guān)性降低到0.96，對(duì)能量通量修正差異可達(dá)1.4 W m?2。究其原因認(rèn)為差異主要來(lái)自于平面擬合修正對(duì)平均風(fēng)向選擇的理論基礎(chǔ)。氣象塔北面由于受到高層樓房的影響，很少有西南氣流，故該風(fēng)向內(nèi)的數(shù)據(jù)較少，導(dǎo)致平面擬合方法出現(xiàn)偏差。在城市下墊面對(duì)通量值進(jìn)行修正，若數(shù)據(jù)較少應(yīng)避免進(jìn)行平面擬合修正。如需精確分析數(shù)據(jù)，需將兩種方法結(jié)合使用。

4.2.3頻率響應(yīng)修正

即便排除人為和環(huán)境造成的損失，超聲風(fēng)速儀和CO2/H2O分析儀本身之間有一個(gè)間距d，使儀器只能測(cè)到渦旋尺度大于2d的湍渦，在低風(fēng)速下高頻出現(xiàn)嚴(yán)重衰減。進(jìn)行頻率損失修正常常可使不同通量分別增加5%到30%左右，夜間修正會(huì)更加明顯。修正公式參照姜明等（2011）。

表2給出不同大氣穩(wěn)定度下通量修正損失率，其中Uw為風(fēng)速u分量和w分量的協(xié)方差，Wt為風(fēng)速w分量和溫度t的協(xié)方差、Wc為風(fēng)速w分量和二氧化碳濃度的協(xié)方差。整體數(shù)據(jù)感熱通量的損失率大約為5.21%，潛熱通量和CO2通量的損失率為9.42%。該通量損失百分比與王介民等（2007）關(guān)于黑河地區(qū)的研究結(jié)果類似。白天的通量損失要明顯小于夜間，即大氣不穩(wěn)定條件下的損失率要小于大氣穩(wěn)定條件下。這是由于白天湍流發(fā)展旺盛，響應(yīng)尺度較大，對(duì)于渦動(dòng)相關(guān)觀測(cè)系統(tǒng)則更容易獲取湍流信號(hào)；夜間則多屬于大氣穩(wěn)定層結(jié)，穩(wěn)定條件下湍流發(fā)展緩慢且物質(zhì)垂直擴(kuò)散緩慢，不易捕捉湍流高頻信號(hào)。

表2　不同大氣穩(wěn)定度條件下通量損失率和修正系數(shù)表Table 2　Comparison of flux loss and correction coefficient under different atmospheric stabilities

選取的數(shù)據(jù)來(lái)自47 m處，每月1日和15日的全天半小時(shí)平均通量來(lái)分析通量損失與平均風(fēng)速和大氣穩(wěn)定度。在圖6中，可以看到在不穩(wěn)定層結(jié)下，即z/L＜0時(shí)，通量損失率隨風(fēng)速的增大而增大，損失率與風(fēng)速呈線性關(guān)系；穩(wěn)定層結(jié)下，即z/L≥0時(shí)，通量損失率與z/L呈現(xiàn)出二階關(guān)系。具體的擬合公式如下：

其中，F(xiàn)Loss為通量損失率，u為水平風(fēng)速，z為觀測(cè)高度，L為莫寧—奧布霍夫長(zhǎng)度。擬合函數(shù)曲線的相關(guān)系數(shù)t檢驗(yàn)?zāi)軌驖M足置信區(qū)間90%，證明該曲線能夠較好的表達(dá)環(huán)境變量與頻率損失的關(guān)系，可為用來(lái)簡(jiǎn)易估算北京城市下墊面上的通量頻率損失，為其他研究提供科學(xué)參考。

4.2.3空氣密度脈動(dòng)修正

空氣密度脈動(dòng)修正是為了修正由于空氣膨脹、壓縮等物理過(guò)程而導(dǎo)致測(cè)量的通量值增大的假象，最常用的有Liu法（Liu，2005）和WPL法（Webb et al.，1980）。為了比對(duì)這兩種方法對(duì)城市下墊面上的適用性，將七層高度上2012年12月至2013年12月的CO2通量值分別用兩種方法進(jìn)行修正，圖7為7月和12月的修正結(jié)果。

圖7中可以看出，經(jīng)WPL法和Liu法修正后，CO2通量分別減小27%和20%，WPL法的修正量更大，通量減小得更多。通量修正較大主要集中在白天，夜間修正量不大。其中峰值和谷值的部分修正幅度更大，尤其在對(duì)應(yīng)的交通早、晚高峰時(shí)段和白天正午前后波谷時(shí)段。圖7中，第一高峰峰值修正前為27.12 mg m?2s?1，經(jīng)過(guò)WPL法和Liu法修正后明顯減小，WPL法的修正幅度要大于Liu的修正幅度，最大幅度為4.29 mg m?2s?1和5.17 mg m?2s?1。夜間經(jīng)過(guò)修正后，碳通量的輸出也有所減小，WPL法的修正量仍然要比Liu方法大，白天和夜間修正前平均CO2通量值分別比修正后的大1.2倍及1.05倍。因WPL方法不考慮水汽的作用，而北京城市環(huán)境下的水汽含量不可忽略，故WPL修正方法偏大，不完全適用。Liu法從理論上對(duì)于水汽的估計(jì)是比較準(zhǔn)確的，故Liu的方法更適合用于北京城市下墊面通量修正。

5　湍流數(shù)據(jù)分析與評(píng)價(jià)

5.1功率譜和協(xié)方差譜檢驗(yàn)

中等大小的渦度既沒(méi)有粘滯效應(yīng)，也沒(méi)有湍流動(dòng)能的生成。這些渦由慣性從較大尺度的渦中獲得能量，又以同樣的方式傳遞給較小的渦，能量向譜下游傳遞的串級(jí)率通常是－2/3（Stull，1988；胡非，1995）。協(xié)譜作為頻率的函數(shù)能夠反映兩個(gè)變量與其各自的譜強(qiáng)的位相關(guān)系，因此協(xié)譜可以用來(lái)檢測(cè)渦度相關(guān)觀測(cè)中與頻率有關(guān)的一些可以導(dǎo)致系統(tǒng)觀測(cè)誤差的變化（Velasco et al., 2005）。確定不同觀測(cè)量的觀測(cè)譜在慣性副區(qū)的斜率對(duì)于判斷渦度相關(guān)儀器的響應(yīng)能力是非常重要的。圖8為47 m處（zm=z－d=6.88 m；z＝47 m，d＝40.12 m）2013年7月1日12:00～12:30經(jīng)過(guò)前處理的數(shù)據(jù)。其中，圖8a中的縱坐標(biāo)為水汽濃度與二氧化碳濃度的功率譜與二氧化碳濃度的比值，圖8b中的縱坐標(biāo)為垂直風(fēng)速w與二氧化碳濃度的協(xié)方差譜與二氧化碳濃度的比值；兩個(gè)圖的橫坐標(biāo)同時(shí)采用歸一化頻率f(zm／u)。

圖5　兩種坐標(biāo)旋轉(zhuǎn)修正后湍流統(tǒng)計(jì)量相關(guān)系數(shù)隨風(fēng)向變化Fig. 5　Comparison of the coefficient wind turbulence statistics in different wind directionsafter two coordinate rotation correction

圖6?。╝）不穩(wěn)定層結(jié)下通量損失隨風(fēng)速變化；（b）穩(wěn)定層結(jié)下通量損失隨大氣穩(wěn)定度變化Fig. 6　(a) Flux loss rate vs wind speed under unstable stratification; (b) flux loss rate vs atmospheric stabilities under stable stratification

從圖8中左圖可以看出，在高頻區(qū)[f(zm／u)＞2]處，CO2濃度和H2O濃度功率譜、T的協(xié)譜均有相同的模態(tài)。在相同的頻率范圍內(nèi)，功率譜變化模態(tài)符合慣性副區(qū)－2/3斜率理論值。圖8中右圖可以看出，在渦度測(cè)量系統(tǒng)中，無(wú)論是超聲風(fēng)速儀還是CO2分析儀，T和CO2濃度都與垂直風(fēng)速w的協(xié)譜表現(xiàn)出一致性，這樣的相似性說(shuō)明渦度相關(guān)測(cè)量設(shè)備沒(méi)有系統(tǒng)性的相移和失真。

5.2湍流穩(wěn)定性檢驗(yàn)

湍流穩(wěn)態(tài)檢驗(yàn)（Stationarity Test）是測(cè)試大氣湍流特征變化的重要方法，它是指湍流統(tǒng)計(jì)量不隨時(shí)間而發(fā)生變化，在觀測(cè)過(guò)程中的非穩(wěn)態(tài)情況會(huì)影響湍流觀測(cè)質(zhì)量。通過(guò)對(duì)平均時(shí)段內(nèi)數(shù)據(jù)進(jìn)行分割，比較整體時(shí)間段內(nèi)的通量與分割段內(nèi)通量單獨(dú)計(jì)算后的平均值的差異的方法（Foken and Wichura，1996；Al-Jiboori et al., 2005）。

舉例說(shuō)明：對(duì)于一個(gè)垂直風(fēng)速為w和某變量x的協(xié)方差時(shí)間序列（30 min），則數(shù)據(jù)長(zhǎng)度為N＝18000（10 Hz）；將這18000個(gè)數(shù)據(jù)分割成4到8個(gè)間斷（通常取5 min為間隔，6個(gè)間斷）；在每個(gè)小時(shí)間短時(shí)間段內(nèi)，協(xié)方差表示為

再對(duì)6個(gè)平均時(shí)段內(nèi)的協(xié)方差進(jìn)行算術(shù)平均計(jì)算，計(jì)算公式為

同樣，全部時(shí)間里的協(xié)方差的平均值為

穩(wěn)態(tài)檢驗(yàn)值（IST）計(jì)算式為

其中，測(cè)試的變量通常為垂直速度w和CO2（或H2O）密度。當(dāng)計(jì)算出來(lái)的穩(wěn)態(tài)測(cè)試比值小于30%認(rèn)為湍流通量處于穩(wěn)態(tài)條件。Foken et al.（2004） 根據(jù)IST進(jìn)行了分類，將湍流穩(wěn)態(tài)測(cè)試劃分為9級(jí)標(biāo)準(zhǔn)，并規(guī)定第1級(jí)數(shù)據(jù)質(zhì)量最好，第9級(jí)數(shù)據(jù)質(zhì)量最差。在這里規(guī)定，1～3級(jí)數(shù)據(jù)為較好數(shù)據(jù)；4～6級(jí)為中等質(zhì)量數(shù)據(jù)；7～9級(jí)數(shù)據(jù)為質(zhì)量較差數(shù)據(jù)，詳見(jiàn)表3。

圖7　修正后CO2通量日變化曲線：（a）12月；（b）7月Fig. 7　Revised CO2flux diurnal variation curves: (a) December, (b) July

圖8 （a）CO2、H2O濃度功率譜圖；（b）垂直風(fēng)速w和CO2濃度、溫度T歸一化協(xié)譜圖。橫坐標(biāo)同時(shí)采用歸一化頻率f(zm／u)Fig. 8　(a) Power spectrum of CO2, H2O concentrations; (b) normalized spectra association of vertical wind speed w, CO2concentration, and temperature T. x-axis: f(zm／u) denotes normalized frequency

表3　穩(wěn)態(tài)檢驗(yàn)值（IST）分類表Table 3　Classification of the value of stationarity test (IST)

利用上述方法對(duì)北京2012年至2013年的CO2濃度數(shù)據(jù)進(jìn)行穩(wěn)態(tài)檢驗(yàn)，并按表3對(duì)數(shù)據(jù)進(jìn)行質(zhì)量等級(jí)標(biāo)志。圖9給出l～9級(jí)各級(jí)數(shù)據(jù)量的分布情況，其中1～3級(jí)數(shù)據(jù)占總量的86%，4～6級(jí)數(shù)據(jù)占數(shù)據(jù)總量的16%，7～9級(jí)數(shù)據(jù)占數(shù)據(jù)總量的3%?？梢钥闯鰯?shù)據(jù)大部分質(zhì)量較好，能滿足穩(wěn)態(tài)檢驗(yàn)的要求。

圖9　穩(wěn)態(tài)檢驗(yàn)統(tǒng)計(jì)圖Fig. 9　The steady-state test

5.3湍流發(fā)展性檢驗(yàn)

湍流的發(fā)展性檢驗(yàn)（總體湍流特征指數(shù)ITC，Integrated Turbulence Characteristics）又稱為湍流積分統(tǒng)計(jì)特性檢驗(yàn)，其物理意義是探討衡量湍流的發(fā)展情況。該檢驗(yàn)以近地層湍流方差相似性理論為基礎(chǔ)，通過(guò)計(jì)算湍流方差相似性來(lái)判斷大氣湍流是否能夠符合近地層Monin-Obukhov相似理論，從而判斷湍流發(fā)展情況。但由于缺少CO2方差相似性的函數(shù)，因?yàn)樵诜讲钕嗨菩詸z驗(yàn)中主要還是針對(duì)風(fēng)速分量u和w，其普適函數(shù)表達(dá)為

式中，x可以是三維風(fēng)速分量u、v、w或溫度T或其他標(biāo)量；則為選取的x量的標(biāo)準(zhǔn)差（可以是、、和或其他標(biāo)量的標(biāo)準(zhǔn)差）為標(biāo)量的相似性參數(shù)，一般選取摩擦速度u*；z為觀測(cè)高度；L為莫寧—奧布霍夫長(zhǎng)度，xf為擬合參數(shù)。利用渦度相關(guān)觀測(cè)值計(jì)算的與相似理論提供的模型值進(jìn)行對(duì)比，就可以計(jì)算湍流數(shù)據(jù)特性積分測(cè)試結(jié)果，稱為ITC值，計(jì)算公式為

同樣參考Foken提出的劃分標(biāo)準(zhǔn)（Foken et al., 2004），1表示最高質(zhì)量數(shù)據(jù)，9表示需要剔除的數(shù)據(jù)，如表4所示。

圖10　湍流積分檢驗(yàn)統(tǒng)計(jì)圖Fig. 10　The turbulence integral test

進(jìn)行湍流發(fā)展特性檢驗(yàn)后，各質(zhì)量等級(jí)的統(tǒng)計(jì)圖如圖10，部分?jǐn)?shù)據(jù)集中在前3等級(jí)，第1等級(jí)的數(shù)據(jù)量占整體數(shù)據(jù)的78%，質(zhì)量差的數(shù)據(jù)較少，不超過(guò)1%。這表明大部分?jǐn)?shù)量能夠符合湍流的方差相似性規(guī)律，質(zhì)量較為理想。值得注意的是9級(jí)的高級(jí)標(biāo)志位代表的質(zhì)量數(shù)據(jù)極度不好，需要剔除的數(shù)據(jù)也存在著不可忽略的一小部分，大多數(shù)異常數(shù)據(jù)都集中在這級(jí)別中，異常數(shù)據(jù)能夠利用分級(jí)標(biāo)志位來(lái)解釋，但也并不代表這一級(jí)別的所有數(shù)據(jù)都是異常數(shù)據(jù)。

表4　湍流積分檢驗(yàn)（ITC）分類表Table 4　Classification of the value of integrated turbulence characteristics (ITC)

5.4綜合等級(jí)評(píng)定

總體質(zhì)量評(píng)價(jià)（QE）是結(jié)合湍流平穩(wěn)性和發(fā)展性檢驗(yàn)對(duì)湍流通量資料做出的質(zhì)量分級(jí)標(biāo)志，F(xiàn)oken et al.（1996）建議的質(zhì)量評(píng)價(jià)標(biāo)準(zhǔn)，如表5。本文研究數(shù)據(jù)的質(zhì)量評(píng)價(jià)結(jié)果顯示，可用于長(zhǎng)期觀測(cè)資料處理研究的數(shù)據(jù)占總數(shù)據(jù)的79%左右（表6）。

表5　綜合湍流數(shù)據(jù)質(zhì)量等級(jí)Table 5 Comprehensive quality grade of turbulence data

表6　質(zhì)量評(píng)價(jià)結(jié)果Table 6　Quality evaluation results

6　結(jié)論

本文主要北京325 m氣象塔2012年12月至2013年12月7層湍流數(shù)據(jù)進(jìn)行前期處理，使用處理后得數(shù)據(jù)進(jìn)行通量計(jì)算，并對(duì)計(jì)算結(jié)果進(jìn)行修正。在此過(guò)程中比對(duì)了不同的處理方案和修正方案。所得結(jié)論如下：

（1）利用累計(jì)次數(shù)Ogive函數(shù)對(duì)數(shù)據(jù)平均時(shí)間進(jìn)行了分析，結(jié)果表明采用30 min作為10 Hz高頻數(shù)據(jù)的平均周期，可以捕獲絕大部分的低頻信號(hào)；

（2）分別使用二次旋轉(zhuǎn)和平面擬合的方法，對(duì)觀測(cè)點(diǎn)非常復(fù)雜的城市下墊面的地形狀況進(jìn)行了校正。結(jié)果表明，二次旋轉(zhuǎn)結(jié)果較好，而平面擬合方法因?qū)δ扯螘r(shí)間內(nèi)的風(fēng)向進(jìn)行平均計(jì)算導(dǎo)致對(duì)地形反應(yīng)較差。對(duì)湍流數(shù)據(jù)進(jìn)行頻率響應(yīng)修正，修正后整體數(shù)據(jù)的感熱通量提高5.21%，潛熱通量和CO2通量均提高9.42%；應(yīng)用WPL法及Liu法對(duì)通量進(jìn)行修正，Liu法更適合北京城市地區(qū)；

（3）從數(shù)據(jù)質(zhì)量的整體來(lái)看，呈現(xiàn)高層優(yōu)于低層、白天好于夜間、夏季好于冬季的變化特點(diǎn)；通過(guò)數(shù)據(jù)評(píng)價(jià)方案的選取，選擇了穩(wěn)態(tài)性檢驗(yàn)和發(fā)展性檢驗(yàn)測(cè)試，分析結(jié)果表明總數(shù)據(jù)中有79%的數(shù)據(jù)為良好，前處理對(duì)數(shù)據(jù)質(zhì)量的提高有幫助。

（References）

Al-Jiboori M H, Fei H. 2005. Surface roughness around a 325-m meteorological tower and its effect on urban turbulence [J]. Adv. Atmos. Sci., 22 (4): 595–605, doi:10.1007/BF02918491.

Burri S, Fre C, Parlow E, et al. 2009. CO2fluxes and concentrations over an urban surface in CAIRO/EGYPT [C]// The Seventh International Conference on Urban Climate. Yokohama.

Coutts A M, Beringer J, Tapper N J. 2007. Characteristics influencing the variability of urban CO2fluxes in Melbourne, Australia [J]. Atmos. Environ., 41 (1): 51–62, doi:10.1016/j.atmosenv.2006.08.030.

Crawford B, Grimmond C S B, Christen A. 2011. Five years of carbon dioxide fluxes measurements in a highly vegetated suburban area [J]. Atmos. Environ., 45 (4): 896–905, doi:10.1016/j.atmosenv.2010.11.017.

Grimmond C S B, King T S, Cropley F D, et al. 2002. Local-scale fluxes of carbon dioxide in urban environments: Methodological challenges and results from Chicago [J]. Environ. Pollut., 116 (S1): S243–S254, doi:10.1016/S0269-7491(01)00256-1.Finnigan J J. 2004. A re-evaluation of long term-flux measurement techniques part II : Coordinate systems [J]. Boundary-Layer Meteorology, 113 (1): 1–41, doi:10.1023/B:BOUN.0000037348.64252.45.

Foken T, Wichura B.1996. Tools for quality assessment of surface-based flux measurements [J]. Agric. Forest Meteor., 78 (1–2): 83–105.

Foken T, Leclerc M Y. 2004. Methods and limitations in validation of footprint models [J]. Agricultural and Forest Meteorology, 127 (3): 223–234.

Hoornweg D, Sugar L, Gómez C L T. 2011. Cities and greenhouse gas emissions: Moving forward [J]. Environment and Urbanization, 23 (1): 207–227, doi:10.1177/0956247810392270.

胡非. 1995. 湍流、間歇性與大氣邊界層 [M]. 北京: 科學(xué)出版社, 288 pp. Hu Fei. 1995.The Intermittent Turbulence and Atmospheric Boundary Layer (in Chinese) [M]. Beijing: Science Press, 288pp.

姜明, 景元書(shū), 郭建俠, 等. 2011. 三維超聲風(fēng)速儀觀測(cè)中風(fēng)向角計(jì)算方法 [J]. 氣象科技, 39 (5): 615–619.Jiang Ming, Jing Yuanshu, Guo Jianxia, et al. 2011. Method for calculating wind direction in measurement using 3D sonic anemometer [J]. Meteor. Sci. Technol. (in Chinese), 39 (5): 615–619, doi:10.3969/j.issn.1671-6345.2011.05.014.

Kotthaus A, Grimmond C S B. 2012. Identification of micro-scale anthropogenic CO2, heat and moisture sources-processing eddy covariance fluxes for a dense urban environment [J]. Atmos. Environ., 57: 301–316, doi:10.1016/j.atmosenv.2012.04.024.

Liu H P. 2005. An alternative approach for CO2flux correction caused by heat and water vapor transfer [J]. Boundary-Layer Meteorology, 115: 151–168.

Liu H Z, Feng J W, J?rvi J, et al. 2012. Four-year (2006–2009) eddy

covariance measurements of CO2flux over an urban area in Beijing [J].

劉樹(shù)華, 劉和平, 李潔, 等. 2005b. 在EBEX-2000實(shí)驗(yàn)資料中湍流耗散率、長(zhǎng)度尺度和結(jié)構(gòu)參數(shù)特征 [J]. 大氣科學(xué), 29 (3): 475–481.Liu Shuhua, Liu Heping, Li Jie, et al. 2005b. Characteristics of turbulence dissipation rates, characteristic length scales, and structure parameters in EBEX-2000 [J]. Chinese Journal of Atmospheric Sciences (in Chinese), 29 (3): 475–481, doi:10.3878/j.issn.1006-9895.2005.03.15.

劉樹(shù)華, 李潔, 劉和平, 等. 2005c. 在EBEX-2000實(shí)驗(yàn)資料中湍流譜和局地各向同性特征 [J]. 大氣科學(xué), 29 (2): 213–224.Liu Shuhua, Liu Heping, Li Jie, et al. 2005c. Characteristics of turbulence spectra and local isotropy in EBEX-2000 [J]. Chinese Journal of Atmospheric Sciences (in Chinese), 29 (2): 213–224, doi:10.3878/j.issn.1006-9895. 2005.02.06.

劉郁玨, 胡非, 程雪玲, 等. 2014. 北京城市通量足跡及源區(qū)分布特征分析 [J]. 大氣科學(xué), 38 (6): 1044–1054.Liu Yujue, Hu Fei, Cheng Xueling, et al. Distribution of the source area and footprint of Beijing [J]. Chinese Journal of Atmospheric Sciences (in Chinese), 38 (6): 1044–1054, doi:10.3878/j.issn.1006-9895.2013.13237.

Matese A, Gioli B, Vaccari F P, et al. 2009. Carbon dioxide emissions of the city center of Firenze, Italy: Measurement, evaluation, and source partitioning [J]. J. Appl. Meteor. Climatol., 48 (9): 1940–1947, doi:10.1175/2009JAMC1945.1.

Miao S G, Dou J X, Chen F, et al. 2012. Analysis of observations on the urban surface energy balance in Beijing [J]. Science China Earth Sciences, 55 (11): 1881–1890, doi:10.1007/s11430-012-4411-6.

Song T, Wang Y S. 2012. Carbon dioxide fluxes from an urban area in Beijing [J]. Atmos. Res., 106: 139–149, doi:10.1016/j.atmosres.2011. 12.001.

宋濤, 修天陽(yáng), 孫楊, 等. 2014. 城市復(fù)雜環(huán)境下渦度相關(guān)通量觀測(cè)的適用性分析 [J]. 環(huán)境科學(xué)學(xué)報(bào), 35 (5): 1089–1098.Song Tao, Xiu Tianyang, Sun Yang, et al. 2014. Applicability of eddy covariance technique to the measurement of turbulent fluxes in complex urban context [J]. Acta Scientiae Circumstantiae (in Chinese), 34 (5): 1089–1098, doi:10.13671/j.hjkxxb.2014.0173.

Stull R B. 1988. An Introduction to Boundary Layer Meteorology [M]. Netherlands: Springer, 401–405.

Velasco E, Pressley S, Allwine E, et al. 2005. Measurements of CO2fluxes from the Mexico city urban landscape [J]. Atmos. Environ., 39 (38): 7433–7446, doi:10.1016/j.atmosenv.2005.08.038.

Vogt R, Christen A, Rotach M W, et al. 2006. Temporal dynamics of CO2fluxes and profiles over a central European city [J]. Theoret. Appl. Climatol., 84 (1–3): 117–126, doi:10.1007/s00704-005-0149-9.

王介民, 王維真, 奧銀煥, 等. 2007. 復(fù)雜條件下湍流通量的觀測(cè)與分析[J]. 地球科學(xué)進(jìn)展, 22 (8): 791–797.Wang Jiemin, Wang Weizhen, Ao Yinhuan, et al. 2007. Turbulence flux measurements under complicated conditions [J]. Adv. Earth Sci. (in Chinese), 22 (8): 791–797, doi:10.3321/j.issn:1001-8166.2007.08.004.

王詠薇, 王雪倩, 胡楠, 等. 2013. 城市地表湍流觀測(cè)數(shù)據(jù)的后處理及質(zhì)量控制 [J]. 氣象科學(xué), 33 (2): 153–159.Wang Yongwei, Wang Xueqian, Hu Nan, et al. 2013. The post-processing and quality control of the urban boundary layer turbulence observation [J]. J. Meteor. Sci. (in Chinese), 33 (2): 153–159, doi:10.3969/2012jms.0116.

Webb E K, Pearman G I, Leuning R. 1980. Correction of flux measurements for density effects due to heat and water vapour transfer [J]. Quart J. Roy. Meteor. Soc., 106 (447): 85–100, doi:10.1002/qj.49710644707.

徐自為, 劉紹民, 宮麗娟, 等. 2008. 渦動(dòng)相關(guān)儀觀測(cè)數(shù)據(jù)的處理與質(zhì)量評(píng)價(jià)研究 [J]. 地球科學(xué)進(jìn)展, 23 (4): 357–370.Xu Ziwei, Liu Shaomin, Gong Lijuan, et al. 2008. A study on the data processing and quality assessment of the eddy covariance system [J]. Adv. Earth Sci. (in Chinese), 23 (4): 357–370, doi:10.11867/j.issn.1001-8166.2008.04.0357.

Data Processing and Quality Assessment of the Eddy Covariance System of the 325-m Meteorology Tower in Beijing

LIU Yujue1, 2, HU Fei1, CHENG Xueling1, and FENG Yongfang2
1 Institute of Atmospheric Physics, Chinese Academy of Sciences, Beijing 100029
2 Fangshan Meteorological Bureau, Beijing 102488

Advances in instrumentation, notably the eddy covariance (EC) technique, offer a tool for the direct measurement of representative flux data in urban areas. The complexity of the urban underlying surface and limitations of the EC system cause the obtained surface fluxes to not be true values. If no correction is applied, the error will be greater than 10%. The applicability of the EC technique, quality assessment of the EC system, and processing of its data, in urban Beijing, were analyzed using the CO2flux data and energy measured at seven levels of the 325-m meteorology tower in the city. Analysis of the data from the tower showed that the best time period to calculate the CO2flux for the urban underlying surface is 30 minutes. The secondary coordinate rotation method is superior to plane fitting. Afterfrequency response correction, the sensible heat flux increased by 5.21%, and the latent heat flux and CO2flux were increased by 9.42%. Liu’s method (Liu, 2005) was found to be better than the Webb E K, Pearman G I, and Leuning R’s (WPL) method (Webb et al., 1980) in the air density pulse correction step. The turbulent spectral checks of the data quality evaluation were satisfactory. After rigorous data screening, the final result showed that approximately 79% of the flux data were good.

Eddy covariance system, Quality control, Quality assessment, Turbulent flux

渦度相關(guān) 質(zhì)量控制 質(zhì)量評(píng)價(jià) 湍流通量

1006-9895(2016)02-0390-11

P442

A

10.3878/j.issn.1006-9895.1504.14325

2014-11-22；網(wǎng)絡(luò)預(yù)出版日期 2015-04-27

劉郁玨，女，1988年出生，博士，主要從事大氣邊界層物理研究。E-mail: lyjsa@163.com

中國(guó)科學(xué)院戰(zhàn)略性先導(dǎo)科技專項(xiàng)——應(yīng)對(duì)氣候變化的碳收支認(rèn)證及相關(guān)問(wèn)題XDA05040301