謝寶良,胡軍國,李燁斐,陳 芳,毛國平

(浙江農林大學信息工程學院,林業(yè)感知技術與智能裝備國家林業(yè)和草原局重點實驗室,浙江省林業(yè)智能監(jiān)測與信息技術研究重點實驗室,杭州 311300)

土壤有機碳庫是陸地碳庫的重要組成部分,總儲量達到1 394 PgC,因此其庫容的微小變化都會對大氣CO2濃度及全球氣候產生巨大影響[1]。然而,土壤呼吸的動力學仍在不斷探究中,導致全球碳通量分析受到很大限制[2-3]。因此針對土壤呼吸的精確測定也成為研究生態(tài)系統(tǒng)碳循環(huán)和地球溫暖化的關鍵問題之一[4]。

土壤CO2濃度擴散包括兩個主要過程:土壤中二氧化碳的產生及其在土壤中的傳輸,并通過土壤表面碳從土壤擴散到大氣的這個過程[5],主要方法涉及以下兩大類,一類是微氣象學法[6],另一類是氣室測定法[7]。微氣象法主要利用渦度相關法,但是易受到較低風速脈動的影響[8],不能準確測定土壤呼吸。氣室法是在土壤表面安裝用金屬或樹脂制作的氣室,根據(jù)氣室內從土壤表面向大氣擴散CO2的速率進而算出土壤呼吸速率的方法[9]。這種測定方法的優(yōu)點是能觀測到小范圍的土壤呼吸特性及其細微的變化。但在受到空間不均一性的影響下,進行大尺度擴展同樣有一定的困難。目前針對氣室法精確測定土壤呼吸的儀器有美國LI-COR公司的LI-8100開路式土壤碳通量測量系統(tǒng)[10],但是價格昂貴,并且隨著測量室對氣體長時間的監(jiān)測獲取會影響土壤CO2的擴散,影響土壤呼吸的監(jiān)測值。Bekku Y等[11]發(fā)現(xiàn)氣室內的環(huán)境與外界環(huán)境會產生一定差異,對測定土壤呼吸產生一定的噪聲。而一般的開放型土壤呼吸監(jiān)測儀在計算土壤呼吸時易受到氣壓、溫度、傳感器質量、濃度梯度、擴散系數(shù)、算法等綜合因素的干擾,導致監(jiān)測獲取的土壤CO2時間序列噪聲過大,而一般的去噪方法無法很好的對非線性土壤CO2時間序列進行去噪分析,導致最終計算得到的土壤呼吸值與LI-8100開路式土壤碳通量測量系統(tǒng)測得的值相差太大,非常不精確,無法用來研究土壤呼吸。因此,若要精確測得土壤呼吸,必須要對收集得到的土壤CO2濃度進行去噪分析。王國英等[12]利用自主研制的設備監(jiān)測獲取了土壤CO2并對其進行小波變換去噪分析,從而獲得較高質量的碳通量數(shù)據(jù),最后再結合合適的算法實現(xiàn)了對土壤呼吸的動態(tài)監(jiān)測。

本文研制了一種開放型的土壤呼吸監(jiān)測儀,如圖一所示,主要由三個傳感器分層組合而成,并結合擴散原理設計了一種新的算法,用來獲取高質量的土壤CO2數(shù)據(jù)。首先利用Fick第二定律計算分析出濃度梯度在垂直方向擴散的滯后時間,并以此對時間序列進行時移,然后運用小波包變換對土壤CO2濃度進行分析去噪,得到的結果利用距離度量[13]來評價CO2濃度時間序列與實際時間序列的相似性。接下來利用Fick第一定律求出碳通量值與LI-8100開路式土壤碳通量測量系統(tǒng)測得的碳通量值進行對比,綜合得出了以下幾個結論:①利用數(shù)據(jù)計算分析發(fā)現(xiàn)了CO2在垂直方向上的擴散規(guī)律與Fick擴散定律吻合;②Fick定律與小波包變換結合處理得到了氣室內任意時刻任意位置的高質量CO2數(shù)據(jù),并利用處理后的數(shù)據(jù)計算得到了與 LI-8100 相近的碳通量值;③通過本文算法及研制的設備可以在未來替代LI-8100,從而節(jié)約大量研究成本。

1 方法介紹

1.1 數(shù)據(jù)來源

土壤呼吸排放的CO2會由高濃度向低濃度擴散,考慮到氣體擴散的規(guī)律,采用描述物質擴散的Fick定律作為建模的理論基礎。該模型設計的測量儀器為一個開放式測量儀,如圖1所示。將測量儀器放置在土壤表面,土壤呼吸產生的氣體首先聚集在儀器底部,當儀器底部的濃度比上部的濃度高時,氣體發(fā)生擴散,最后擴散出儀器,擴散過程如圖1中紅色箭頭所示。儀器內部傳感器每隔一段時間采集一次數(shù)據(jù),然后將采集的數(shù)據(jù)發(fā)送到服務器端,最后利用Fick定律建模并計算出土壤碳通量。

圖1 土壤呼吸監(jiān)測儀

1.2 Fick定理

Fick定律是描述物質擴散現(xiàn)象的宏觀規(guī)律,這是生物學家Fick于1855年發(fā)現(xiàn)的。Fick定律包括第一定律和第二定律。Fick定律在各個行業(yè)發(fā)揮著重要作用[14-16]。第一定律描述的是在單位時間內通過垂直于擴散方向的單位截面積的擴散物質流量(擴散通量),用J表示與該截面處的濃度梯度成正比,即

(1)

Fick第二定律描述的是在非穩(wěn)態(tài)擴散過程中,在距離x處濃度隨時間的變化率等于該處的擴散通量隨距離變化率的負值,即

(2)

式中:C為擴散物質的濃度,t為擴散時間,x擴散的距離,D為擴散系數(shù),如式(3)所示。

(3)

式中:T為熱力學溫度,P為大氣壓強,μA、μB為氣體的分子量,本文研究CO2氣體在空氣中擴散,因此μA=44,μB=29。VA、VB為氣體A、B在正常沸點時液態(tài)克摩爾容積,VA=34,VB=29.9。從式(3)可以知道擴散系數(shù)D只與測量環(huán)境的溫度和大氣壓強有關,跟擴散氣體的濃度無關。測量儀器可以實時測量當前環(huán)境的溫度和大氣壓強,實時更新擴散系數(shù)D。

本文研究的氣體擴散問題屬于非穩(wěn)態(tài)擴散,因此利用Fick第二定律對土壤呼吸排放的二氧化碳進行計算。得到如下方程:

(4)

式中:C為CO2濃度;t為CO2在監(jiān)測儀內的擴散時間;x為CO2在監(jiān)測儀內的擴散距離;D為CO2在監(jiān)測儀內的擴散系數(shù);C0為擴散到監(jiān)測儀頂部的CO2;Cs為土壤擴散到監(jiān)測儀內的CO2。

解得:

(5)

式中:erf(·)為誤差函數(shù),

(6)

1.3 小波包變換

由于本裝置傳感器精度較低,在實驗過程中易受外界環(huán)境影響,測得的CO2濃度值存在著較大的噪聲。因此通過小波包變換對非線性時間序列進行分析去噪,來獲取更高質量的CO2濃度時間序列,提高土壤呼吸監(jiān)測儀的監(jiān)測精度。

(7)

式中:S為信號;A為低頻;D為高頻。

圖2 小波包分解示意圖

(8)

式中:h,g為濾波器系數(shù);d為小波包分解系數(shù);p,t為分解層數(shù);j,k為小波包節(jié)點號。

1.4 評價指標

利用第2節(jié)中的2017年10月份采集的浙江省杭州市臨安區(qū)武肅街666號的6個傳感器采樣間隔為10 s的二氧化碳數(shù)據(jù),進行小波包變換處理。為了說明小波包變換在土壤二氧化碳去噪分析上的優(yōu)越性,利用下列公式來說明方法的優(yōu)越性:

①斜率距離

首先對時間序列進行分段線性化,主要由重要點來確定分段。本文將極值點作為重要點對時間序列進行分段線性化,步驟如下:

若給定時間序列X=(x1,x2,…,xn),參數(shù)R>1,若xm∈X(1

上述步驟給出了xm是重要極大值點滿足的條件:xm為子序列X[i:j]中的最大值點,且xm與子序列兩個端點的比值大于參數(shù)R。通過參數(shù)R可以控制時間序列分段的精細度,R值越大選中重要點越少,分段越粗,反之分段就越細。類似的,重要極小值點定義如下:

若給定時間序列X=(x1,x2,…,xn),參數(shù)R>1,若xm∈X(1

對于獲得的極值點,我們將其進行按順序排列,得到X=(x1,s1,s2,…,sk,xn),對極值點間進行線性擬合,得到如下一元線性分段函數(shù)模型:

(9)

式中:bi為直線的斜率估計值,i=1,2,…,k,代表分段趨勢大小;ai為直線的截距估計值,i=1,2,…,k;ei(t)為零均值白噪聲,i=1,2,…,k,表示一段時間內時間序列與它的分段線性表示之間的誤差。

通過以上步驟,可以獲得如下公式:

S=([b1,s1-x1],[b2,s2-s1],…,[bk,sk-sk-1],
[bn,xn-sk])

(10)

這里我們將公式稍作修改方便后面的計算,令:

(11)

則可得:

S=([b1,l1],[b2,l2],…,[bk,lk],[bn,ln])

(12)

接下來,利用上述公式計算斜率距離來分析時間序列A和時間序列B的相似性。在本文中,斜率距離越小,代表兩時間序列越接近。公式如下:

(13)

②信噪比SNR(Signal-to-Noise Ratio),公式定義為:

(14)

③原始信號與降噪信號之間的均方根誤差RMSE(Root Mean Squared Error),定義如下:

(15)

2 Fick-WPT方法介紹

在利用Fick定律計算由開放型氣室法測得的土壤CO2值時,存在著以下幾個問題:①由開放型氣室測得的土壤CO2時間序列含噪聲過大,會影響后續(xù)的Fick定律計算得到的碳通量值;②氣室中各個高度的CO2值不相同,并且隨著土壤CO2擴散,各層傳感器在何時接收到CO2也不相同,導致無法辨識時間序列的去噪時間段。針對上述兩個問題,本文利用小波包變換在處理時間序列噪聲的優(yōu)勢上,結合Fick第二定律能計算出CO2隨著時間和空間的擴散,提高了開放型氣室法測得的CO2濃度值的質量。

步驟1首先將圖一中的a層和c層傳感器測得的CO2濃度時間序列μ(a,t)和ψ(c,t)代入式(5)中,得到:

(16)

令式中x=b,則可得到b層傳感器的理論CO2濃度時間序列,即

(17)

步驟2通過Fick第二定律計算出CO2濃度在豎直高度上的擴散變化,通過連接轉折處,分析得到CO2濃度在豎直高度上隨時間變化的關系圖,擬合得到CO2濃度在豎直高度上隨時間變化的關系式,即

x=kt+b

(18)

式中:x為高度,t為時間。

利用式(18)作為門限標準選取需要去噪的時間段,即

x(t)=x(T)-x(t′)

(19)

式中:T∈[0,N],t′∈[0,(x-b)/k]。

步驟3將新獲得的時間序列作為去噪對象,對每一高度處的CO2濃度時間序列C(x,t)進行小波包變換。步驟如下:①時間序列的小波包分解。選擇一個小波并確定一個小波分解的層次N,然后對時間序列C(x,t)進行N層小波包分解。②計算最佳樹(即確定最佳小波包基)。對于一個給定的熵標準來計算最佳樹。③小波包分解系數(shù)的閾值量化。對每一個小波包分解系數(shù),選擇一個適當?shù)拈撝挡ο禂?shù)進行閾值量化。④小波包重構。根據(jù)第N層的小波包分解系數(shù)和經過量化處理系數(shù),進行小波包重構。

步驟4通過重構后的時間序列計算氣室內的濃度梯度,并將其代入Fick第一定律里,計算得到土壤碳通量。

3 實驗結果與分析

3.1 Fick第二定律對CO2濃度時間序列的分析

數(shù)據(jù)選取2017年12月份同一研究區(qū)域的CO2濃度時間序列,將a層傳感器和c層傳感器測得的CO2濃度值代到Fick第二定律中,求得隨著時間變化,CO2濃度從a傳感器擴散到c傳感器的變化規(guī)律,如圖3所示。

圖3 隨時間變化的CO2濃度在豎直高度上的擴散變化圖

圖4 二氧化碳濃度隨時間在豎直高度上的變化關系圖

從圖3可得,由Fick第二定律計算出的CO2在垂直高度擴散時,需要一定的時間。而濃度梯度則隨著時間的增加不斷的變小。這里,我們將圖三曲線各個轉折處的點坐標重新描繪到一個平面坐標系中,如圖4所示,可以得到CO2濃度隨時間在豎直高度上的關系式為x=0.001t+0.31,說明CO2在開放型監(jiān)測儀內是較為均勻擴散的,擴散10 cm所需要的時間大約為90 s左右。在實驗中,由于土壤呼吸監(jiān)測儀放置在土壤上方后不是立即測量,存在著一定的系統(tǒng)誤差。并且測量儀內本身含有CO2,導致3個傳感器在同時測量時已經發(fā)生擴散,與理論值發(fā)生一定的偏差,如圖五所示。因此,本文將由Fick第二定律計算分析得到的CO2濃度隨時間在豎直高度上的變化關系作為門限進行時移來確定去噪時間段。

圖5 二氧化碳的理論值與實測值之間的對比

3.2 Fick-WPT方法與其他方法之間的對比

為了更好的說明本文算法的優(yōu)越性,實驗選取了同一天內不同時間段以及不同環(huán)境下的CO2濃度時間序列進行測試。首先為了驗證儀器的可行性,將LI-8100儀器橫放置于封閉室內地面,同時將本實驗裝置土壤呼吸監(jiān)測儀也橫放置于封閉室內地面進行監(jiān)測。LI-8100儀器顯示CO2濃度均值477.4×10-6,碳通量值為-0.03 μmol/(m2/s),土壤呼吸監(jiān)測儀的三個傳感器的CO2濃度均值分別為505.77×10-6,501.26×10-6,500.41×10-6,利用Fick定理計算出來的碳通量值為0.01 μmol/(m2/s)?？梢园l(fā)現(xiàn),雖然通過本實驗裝置計算得到的碳通量值相近,但可能具有一定的偶然性,而且傳感器測得的數(shù)據(jù)質量并不是很理想,因此為了克服這一問題,選擇合適的算法,可以使本裝置應用于土壤呼吸的監(jiān)測,具有非常大的意義。本文選取了2017年12月份的一組實驗的數(shù)據(jù),利用MATLAB軟件及其中的wavemenu進行處理。根據(jù)第三節(jié)所述算法,這里我們的小波包變換最終選取的最佳小波包基為sym5,并且對時間序列進行6層分解,閾值選取的是threshold,數(shù)值為2。首先,我們對b層傳感器測得的CO2濃度時間序列進行小波包分析,得到如圖6所示,并且得到了如圖7所示的剩余誤差。接下來我們以處理之后的b層數(shù)據(jù)作為標準與其他方法計算得到的數(shù)據(jù)進行對比。然后,我們將a層傳感器和c層傳感器進行分析計算,分別利用Fick定理、Fick定理算出之后由小波包分析直接處理和Fick定理與小波包分析結合的方法算出的b層傳感器的CO2濃度時間序列與標準進行對比。從圖8和表1可以看到,無論從圖像、RMSE、信噪比和時間序列相似性上看,本文算法都比其他方法優(yōu)秀。從圖8看,利用本文算法得到的時間序列非常平滑,與標準時間序列及實際監(jiān)測得到的時間序列都非常接近。而利用其他方法得到的時間序列,雖然在700 s左右之后與實際時間序列和標準時間序列接近,但是在0～700 s左右與實際和標準時間序列相差甚遠。說明,若沒有利用Fick第二定律分析得到的時差,會導致前期的時間序列與實際不符。從表1可以看出,本文算法得到時序信噪比為8.340 4,接近標準時序信噪比14.351 1,而其他方法得到的時序則出現(xiàn)了負值,明顯噪聲過大,數(shù)據(jù)質量不利于后續(xù)的土壤呼吸研究;RMSE上,本文算法僅為6.279 9與標準時序3.143 5接近,而其他算法計算得到的RMSE超過了20,與標準相差過大;時序相似性上,本文算法與先小波后Fick計算得到的時間序列與標準相比分別為0.217 4和0.221 6,而Fick計算得到的為0.440 9,與標準相似度過低。綜上,利用本文算法,可以更好的結合Fick第二定律,精確算出氣室內各個位置各個時間點的CO2濃度值。

表1 驗證結果

圖6 二氧化碳濃度時間序列對比

圖7 處理之后得到的剩余誤差

圖8 算法比較

3.3 Fick-WPT方法在不同時間段計算得到的碳通量值與傳統(tǒng)LI-8100對比

為了探究算法得到的數(shù)據(jù)是否具有更加實際的意義,本文分別在2018年8月9號的上午9點到10點半,中午12點到1點半,下午6點到7點半三個時間段對研究區(qū)域進行了土壤呼吸監(jiān)測的對比實驗。本次實驗為了提高數(shù)據(jù)的質量,分別將3個傳感器放置于開放型監(jiān)測儀內的0.1 m,0.2 m,0.3 m處進行監(jiān)測土壤CO2濃度時間序列,每一個傳感器看作一個氣室,之間的距離為1 cm。首先通過本文算法,計算出0.2 cm處的理論值,并和實際值及用Fick定理算出來的理論值進行對比,如圖9所示。

圖9 早中晚二氧化碳濃度時間序列對比

從圖9可以看到,圖9(a)(上午)和圖9(c)(傍晚)的計算結果非常優(yōu)秀,圖9(b)(中午)的計算結果相對較差,可能是因為溫度影響了最終計算結果。但總體來說最終結果還是非常出色的,說明通過本文算法,可以非常精確的計算出氣室內各個位置,各個時間段的CO2濃度值。因此,本文利用算法進一步計算出土壤的碳通量值,并與其他的評價指標來分析算法更廣泛的實際意義,如表2所示。

從表中可以看到,首先,經過本文算法處理得到的碳通量值明顯與LI-8100測得的碳通量值非常接近,并且比未處理及Fick定理計算得到的更加優(yōu)秀;其次,針對CO2濃度時間序列,經本文算法處理之后的時間序列與LI-8100監(jiān)測得到的時間序列更為相似,而未處理的時間序列及Fick定理計算得到的理論時間序列與LI-8100相比,相似度過低,從圖9中也可以發(fā)現(xiàn),這兩組時間序列起伏波動非常大,穩(wěn)定性過低;最后從平均值角度來看,未處理的時間序列,用Fick定理及本文算法處理的時間序列的濃度平均值都比LI-8100監(jiān)測得到的濃度平均值低,但相比較而言,用Fick定理及本文算法處理得到的濃度平均值略優(yōu)于未處理的濃度平均值。綜上所述,經過本文算法處理,不僅能夠精確的算出各個時間點各個氣室位置的CO2濃度值,而且算出來的碳通量值也與LI-8100測得的非常接近。因此,利用本文算法,能夠非常好的克服傳感器的質量問題以及開放型氣室所存在的隱含問題,同時所測得的土壤呼吸值又具有非常重要的實際意義,利用較廉價的實驗設備獲取了非常高的實驗結果,對土壤呼吸測量具有非常大的意義。

表2 驗證結果

其中碳通量值的單位為μmol/(m2/s),平均值單位為10-6。

4 討論

Coifman和Wickerhauser于20世紀90年代初提出了小波包WP(Wavelet Packet)分析,不僅對低通子帶進行分解,而且也對高通分量分解,從而聚焦到感興趣的任意頻段,突破了小波分析對信號頻帶進行等Q劃分的局限性。因此,利用小波包變換可以有效的對非平穩(wěn)時間序列進行去噪分析,實現(xiàn)一般去噪方法無法得到的效果,并且已經在各領域取得了一定的成就[17-21]。而小波包變換結合其他方法能夠更加合理的解決非線性時間序列數(shù)據(jù)中產生的問題,這種方法已經應用到土壤學領域,取得了比單獨一種方法更加精確有效的結果[22-23]。

Fick擴散定律由阿道夫·菲克于1855年提出,能夠在不依靠宏觀的混合作用發(fā)生的傳質現(xiàn)象時,描述分子擴散過程中傳質通量與濃度梯度之間關系的定律?；贔ick定律的預測模型,并利用其他針對性的方法對數(shù)據(jù)進行處理的方式,能夠取得比單一Fick定律更加精確有效的結論[24-25]。

本文結合上述兩種方法的優(yōu)點,首先基于Fick定律能描述非穩(wěn)態(tài)擴散過程中土壤CO2濃度時間序列的變化,計算分析出需要去噪的時間段。接下來利用小波包對時間序列的高頻信號和低頻信號逐層分解,利用軟閾值對數(shù)據(jù)進行去噪,最后重構得到質量較高的時間序列。通過分析比較時間序列以及通過時間序列計算得到的土壤碳通量值,進一步證明了該方法的有效性和可靠性。但是在小波包基的選取上和氣室間的距離,還需要進一步探究。

5 結論

本文通過結合小波包變換和Fick定理將其運用到計算土壤CO2濃度時間序列,獲得了比單獨使用小波包變換和Fick定理更加精確的時間序列。同時利用本文算法求得的時間序列,運用到Fick第一定律中,計算得到的了與LI-8100測得的相近的土壤碳量值。通過廉價的實驗設備和適當?shù)乃惴?獲得了與費用昂貴的土壤呼吸監(jiān)測儀相近的土壤呼吸值,具有非常重要的意義。但是,本文自制的實驗設備在野外測量仍會受到較大的環(huán)境因素影響,因此改進實驗設備使其具有在野外更好的適應能力在以后的研究中可以進一步考慮。