基于SiB2模型的低丘紅壤區(qū)花生根域土壤水分動態(tài)模擬與數(shù)據(jù)同化

張　乾，景元書

(1.氣象災害預報預警與評估協(xié)同創(chuàng)新中心/江蘇省農(nóng)業(yè)氣象重點實驗室，南京 210044；2.南京信息工程大學，南京 210044)

0　引　言

土壤水分是研究作物生長與地氣相互作用之間聯(lián)系的重要指標，由于其在空間和時間上變化非常大，因此，建立準確的陸面過程模型評估陸面條件尤其是土壤水分變化在陸地大氣相互作用中的表現(xiàn)，研究大尺度和流域尺度的水分和能量循環(huán)，獲取完備的地球表層系統(tǒng)的時空信息[1]，對提升區(qū)域尺度天氣預報精確性有著積極的意義。根據(jù)基本氣象資料對若干段時間內(nèi)作物根域土壤水分的變化進行準確的預估，有利于更深層次地研究和理解地球表層生物物理過程。同時，陸面模型在進行模擬時，模擬誤差會隨著時間不斷積累，因此可以建立數(shù)據(jù)同化方案，在觀測時間段內(nèi)對模型的模擬進行校正以提高模擬精度。

Sellers 等人于1996年建立的SiB2(Simple Biosphere Model version 2，簡單生物圈模型)，是目前研究單層植被中水、能量以及CO2通量變化較成熟的陸面模型之一。SiB2在全球尺度下使用衛(wèi)星圖像，獲得了多種類型植被的物候?qū)W資料作為模型運行的重要參數(shù)，應用能量、動量和質(zhì)量守恒定律, 能夠方便有效地模擬出土壤、大氣和生物圈諸多參數(shù)以及陸面(植被和裸地)和大氣之間的相互作用，該模式已被嵌入其他區(qū)域和全球模式中[2]。

南方低丘紅壤區(qū)是我國水土流失較為嚴重的地區(qū)之一，由于紅壤性質(zhì)上的酸、瘦、黏等弱點，紅壤分布區(qū)域降水時空分布的不均勻，高溫伏旱頻繁，以及不合理開發(fā)利用造成的水土流失、土壤退化、土壤污染等，導致紅壤地區(qū)的生態(tài)環(huán)境惡化，紅壤資源潛在的生產(chǎn)能力得不到應有的發(fā)揮，使整個地區(qū)農(nóng)業(yè)及經(jīng)濟持續(xù)發(fā)展受到嚴重影響[1]。花生是我國南方紅壤地區(qū)重要的經(jīng)濟作物之一，利用模型以及相應的數(shù)據(jù)同化方案準確地模擬花生根域土壤水分變化，有利于指導當?shù)剞r(nóng)業(yè)生產(chǎn)，趨利避害。

1　研究地區(qū)與研究方法

1.1　研究區(qū)概況

研究區(qū)位于江西省鷹潭市余江縣下轄劉墾場三分場的紅壤生態(tài)試驗站( 116°55′E，28°15′N) 。該地區(qū)地形以低丘崗地為主，南北有少量丘陵，平均海拔37.56 m，坡度3°～8°，土壤多為紅砂巖或紅黏土發(fā)育的普通紅壤[4]，是典型的我國南方低丘紅壤區(qū)之一。屬于亞熱帶濕潤季風氣候，四季分明、氣溫偏高。年平均氣溫17.6 ℃，年平均降水量1 788.8 mm，主要集中在春末夏初時期，年平均日照時數(shù)1 739.4 h。研究區(qū)域內(nèi)，經(jīng)濟作物以水稻、花生、柑橘、板栗為主，植被分布主要以樟科、茶科等常綠闊葉林以及部分竹類、藤本、蕨類等構成。

1.2　數(shù)據(jù)來源及處理

土壤濕度和地表溫度數(shù)據(jù)來源于鷹潭市余江縣劉墾場三分場的FDR水分傳感器觀測的2015年花生地的水分數(shù)據(jù)。進行統(tǒng)計分析前，對土壤水分數(shù)據(jù)進行單樣本k-s檢測，其中6-7月的數(shù)據(jù)采用平方根反正切函數(shù)進行轉(zhuǎn)化，最終使所有數(shù)據(jù)均符合正態(tài)分布，可以進行統(tǒng)計分析[5]。根據(jù)對花生地實地測量，得到適用于試驗區(qū)土壤特性和植被條件的 SiB2所需要的全部參數(shù)。

1.3　試驗設計

在試驗區(qū)代表性地段，選擇上坡、下坡建造的20 m×5 m(100 m2)的試驗小區(qū)，種植的花生品種為粵油551[6]，小區(qū)由水泥板圍成，水泥板高出地面30 cm，以阻擋小區(qū)外徑流流入[4]。FDR水分傳感器每30 min測定一次10～100 cm深度的土壤溫度以及水分含量。分別選取5月10-19日、6月8-17日、8月17-26日3個均為期10 d的時間段作為花生幼苗期、開花期、飽果成熟期生長階段的代表時間段，將3個時間段的逐小時深度分別為10、15、20 cm的土壤濕度作為觀測值，分別在SiB2中代入測定好的3個時間段的花生形態(tài)學參數(shù)、模型主程序運行參數(shù)等，選擇5月10日、6月8日、8月17日3個日期0時刻儀器測量得到的植被溫度、地表溫度、深土溫度、冠層空氣溫度、參照高度溫度、土壤表層水分含量、植被根域水分含量、深層土壤水分含量數(shù)據(jù)分別作為花生幼苗期、開花期、飽果成熟期SiB2運行的初始值，最后代入相應的自動氣象站氣象資料。其中，修正SiB2長波輻射的計算，得到土壤水分模擬值后，引入集合卡爾曼濾波進行數(shù)據(jù)同化以調(diào)高模擬精度。

1.4　誤差統(tǒng)計

對于模型模擬和同化結(jié)果，采用均方根誤差(RMSE, Root Mean Square Error)、平均誤差(AVE, Average Error)和相對誤差(RE, Relative Error)三種誤差分析方法[7]。

(1)

(2)

(3)

式中：N為整個觀測期的時長；Obst為t時刻的觀測值；Xt為t時刻的模型模擬和同化結(jié)果[7]。

2　SiB2和數(shù)據(jù)同化方案

2.1　SiB2簡介

陸面過程作為氣候-陸地系統(tǒng)間重要的反饋機制之一，是天氣、氣候、水文和環(huán)境預報的關鍵因素。針對不同的陸地下墊面，具有的地表反照率、比輻射率、地面粗糙度等理化性質(zhì)以及相應的植被類型各不相同，局地氣候受其影響的敏感性也不同的陸面過程難點，1996年，Seller等人在SiB1的基礎上建立了SiB2，融入了更科學的植被冠層光合作用傳導模型，分別描述CO2和水汽在葉片內(nèi)部和外部的傳輸過程，克服了SiB1單純使用經(jīng)驗公式的缺陷。Sellers還對水文子模型進行了改進，以確切地描述基本徑流，并對土壤剖面內(nèi)的層間交換進行了更可靠的計算，模式定義了9種植被類型和7種土壤類型，在模擬地表能量收支、徑流傳輸?shù)热〉昧撕玫男Ч?/p>

SiB2 利用氣象驅(qū)動數(shù)據(jù)，可以模擬一定時間一定區(qū)域內(nèi)的凈輻射、感熱、潛熱、地表熱通量、地表有效輻射溫度、土壤濕度變化。在模擬地面以下土壤的熱傳導和水熱輸送過程中，SiB2定義了表層、根域和深層3 層土壤。公式(4)～(6)為每層土壤的水分控制方程[7]。

(4)

(5)

(6)

式中：θi(i=1,2,3)分別為表層、根區(qū)、深層土壤中液態(tài)水的體積含水量，m3/m3；ρw為液態(tài)水的密度，kg/m3；Di(i=1, 2, 3)為每層土壤厚度，m；Qi,i+1(i=1,2)為第i和i+1 層之間的水流，m/s；Q1為上邊界進入土壤表層的水流，m/s；Q3為重力排水，m/s；Eg為土壤表層的蒸發(fā)速率，m/s；Ect為通過氣孔的植被蒸騰速率(m/s)[7]。

2.2　參數(shù)設置

根據(jù)試驗地條件，將SiB2中的植被類型設定為9(農(nóng)田/C3草地)，土壤類型設定為3(壤土)[9]，在實地測量的基礎上確定了在花生幼苗期、開花期和飽果成熟期內(nèi)低丘紅壤區(qū)農(nóng)田小流域適于SiB2運行的所有參數(shù)，見表1～表3。

2.3　長波輻射修正

運行SiB2所需要的驅(qū)動數(shù)據(jù)包括向下短波輻射、向下長波輻射、水汽壓、氣溫、風速、降水6個基本量，其中鷹潭市余江縣紅壤生態(tài)試驗站自動氣象站提供了短波輻射、相對濕度、氣溫、風速、降水量的逐小時觀測資料。

水汽壓使用相對濕度換算得到，T為氣溫。

(7)

向下長波輻射是來自整層大氣的輻射，即大氣逆輻射。SiB2中可以根據(jù)已有驅(qū)動數(shù)據(jù)換算向下長波輻射，但我們選擇利用江西省鷹潭市氣象局自動氣象站的逐小時云量數(shù)據(jù)，修正brunt公式從而計算得到地面有效長波輻射。

表1　花生不同生長期形態(tài)學參數(shù)

表2　花生和土壤光學特性參數(shù)

表3　試驗區(qū)作物生理參數(shù)及土壤特性參數(shù)

(8)

則向下長波輻射為：

(9)

式中：ε為地面比輻射率，5、6、8月的值分別為0.97、0.97、0.975[10]；σ為斯蒂芬玻爾茲曼常數(shù)；T、d為地面溫度；C為春夏季江西地區(qū)云對地面有效輻射的減弱系數(shù)(0.59)[13]；N為逐小時云量。

2.4　數(shù)據(jù)同化方案

數(shù)據(jù)同化的基本思想是假設模型的模擬結(jié)果和觀測都有一定的不確定性，通過一定的數(shù)據(jù)同化方法融合這兩種來源的數(shù)據(jù)可以得到一個更優(yōu)的模型狀態(tài)估計結(jié)果[16]。傳統(tǒng)的卡爾曼濾波算法是一種利用觀測數(shù)據(jù)和預測值(兩者都含有噪聲和誤差)來逼近隱藏變量真實值的算法，模型預測值被視為先驗值，觀測值和模型預測值一起用于獲得隱藏變量的更新估計值，即為濾波過程。由于在利用氣象模型進行數(shù)值模擬的過程中，模擬誤差會不斷積累，最終導致時間序列越長，結(jié)果越不可靠，解決誤差積累的核心問題是如何在模型運行過程中引入更準確的初值作為校正。所以我們引入集合卡爾曼濾波算法，在整個模擬過程中設置多個同化點，每次同化引入觀測值進行濾波，得到的更逼近真實值的土壤水分數(shù)據(jù)作為下一階段的模型初值。

集合卡爾曼濾波是順序數(shù)據(jù)同化算法的典型代表之一，用狀態(tài)樣本的集合來代表模型狀態(tài)后驗概率密度函數(shù)，通過這些樣本的向前積分計算狀態(tài)的均值和方差，克服了卡爾曼濾波線性化的缺點[17]，從而可以應用到陸面過程等非線性系統(tǒng)模型中，其計算步驟如下：

(10)

ωi,k～N(0,Qk)

(3)計算k+1時刻的卡爾曼增益矩陣Kk+1：

(11)

(12)

(13)

(14)

(15)

vi,k～N(0,Rk)

(16)

(17)

(18)

(5)進入下一時刻返回步驟(2)。

集合卡爾曼濾波主要包括預報和分析兩個過程。首先進行模式狀態(tài)的預報，接著引入觀測資料，然后根據(jù)觀測資料對模式狀態(tài)進行重新分析，公式(10)為變量的集合預報，公式(13)、(14)為預測誤差的分析，公式(16)為對變量的分析，公式(18)為對分析變量的誤差分析。設定集合數(shù)為48，同化間隔為6 h，SiB2模型誤差方差矩陣(Q)和觀測誤差方差矩陣(R)分別為[7,18]：

3　結(jié)果與分析

3.1　模擬與同化結(jié)果

各生育期模擬與同化結(jié)果見圖1-圖3。

圖1　幼苗期根域土壤水分模擬與同化結(jié)果

圖2　開花期根域土壤水分模擬與同化結(jié)果

圖3　飽果成熟期根域土壤水分模擬與同化結(jié)果

幼苗期降雨量多且集中，對花生根域土壤濕度影響較大，短期的大量突發(fā)性降水使SiB2對根域土壤濕度估計過高，在最初階段就積累了大量誤差，導致模擬結(jié)果整體偏高。同化結(jié)果雖然在有降水時也產(chǎn)生了過高估計，但隨著時間推移可以不斷消減誤差累積，在無降水或降水偏少時，同化結(jié)果接近觀測值。

開花期由于降水減少，同時高溫伏旱增多，SiB2對根域土壤濕度估計過低，誤差隨著時間逐漸增大。同化則取得了很好的效果，整體略低于觀測值。

飽果成熟期降水等氣象因素對模擬結(jié)果的干擾明顯減小，模擬和同化結(jié)果與觀測值都保持了較好的一致性，同化結(jié)果更接近觀測值。

3.2　誤差統(tǒng)計

對試驗區(qū)幼苗期、開花期、飽果成熟期花生地分別為10、15、20 cm深度的根域土壤水分的模擬結(jié)果顯示(見表4)，SiB2對飽果成熟期的花生根域土壤濕度模擬精度最好，整體偏低；對開花期的模擬精度次之，整體偏低，其離散程度低于飽果成熟期模擬結(jié)果；對幼苗期的模擬結(jié)果較差，其整體偏高，離散程度較高。在利用集合卡爾曼濾波算法進行數(shù)據(jù)同化后，均方根誤差分別降至原先的1/4、1/4、1/2；平均誤差分別降至原先的1/5、1/5、1/2；相對誤差分別降至原先的1/10、1/5、1/2。

表4　花生不同生長期模擬結(jié)果誤差比較

4　結(jié)　語

本文模擬了低丘紅壤區(qū)花生幼苗期、開花期、飽果成熟期3個生長階段花生根域深度土壤濕度10 d內(nèi)的變化情況，并發(fā)展了基于集合卡爾曼濾波數(shù)據(jù)同化方案，得出如下結(jié)論：

(1)針對江西低丘紅壤區(qū)花生地土壤性質(zhì)、作物生物學特性等進行的實地測量，以及對SiB2長波輻射計算的改進，使得SiB2可以初步適用于對當?shù)赜酌?、開花、飽果成熟等關鍵生長期花生作物根域土壤水分的短期動態(tài)變化模擬，在經(jīng)過數(shù)據(jù)同化后，3個生長期根域土壤濕度模擬結(jié)果均方根誤差分別降至原先的1/4、1/4、1/2；平均誤差分別降至原先的1/5、1/5、1/2；相對誤差分別降至原先的1/10、1/5、1/2。

(2)模型的模擬和同化結(jié)果對降水等氣象因素都具有較好的響應能力，短期強降水會使得模擬結(jié)果偏高，長期的無降水或伏旱會使得模擬結(jié)果逐漸變低，類似的極端天氣現(xiàn)象如果發(fā)生在早期或者持續(xù)較長時間，則會使得模擬結(jié)果積累大量誤差，在進行數(shù)據(jù)同化后則可以有效地避免這種情況。

(3)在花生不同生長期內(nèi)，根域土壤水分動態(tài)變化都表現(xiàn)了不同的特點：幼苗期花生生長發(fā)育迅速，生物學參數(shù)不穩(wěn)定影響模型計算，且由于花生枝干發(fā)育不完全，保持水土的能力不強，受氣象因素干擾較大，春夏之交降水多發(fā)，因而模擬結(jié)果精度和一致性較差，進行長時間的模擬不現(xiàn)實；開花期花生生長發(fā)育速度變慢，植株枝干基本發(fā)育完全，有一定的水土保持能力，受氣象因素影響較小，模擬精度以及一致性均好于幼苗期，但夏季多發(fā)伏旱，降水缺乏使得模擬結(jié)果雖然一致性較高，但整體結(jié)果偏低；飽果成熟期花生發(fā)育完全，植株枝干生長茂密、葉面積指數(shù)最大，具有較強的水土保持能力，不易受到氣象因素的影響，模擬精度最佳，一致性次于開花期。所以對幼苗期和開花期的花生根域土壤濕度進行模擬時，有必要進行數(shù)據(jù)同化。

試驗中造成誤差的原因可能有，一是實際測量不同生長期花生的形態(tài)參數(shù)如葉面積指數(shù)、冠層覆蓋比、平均葉長葉寬等誤差較大，尤其是在花生幼苗期內(nèi)，花生植株發(fā)育迅速，葉片、枝干、根系等生長快，最初設定的形態(tài)參數(shù)隨著時間的推移會越來越不準確導致誤差偏大；二是小區(qū)周圍有較濃密高大的樹林遮蔽了部分陽光，阻礙了風的影響，引起了誤差，突發(fā)性降水、伏旱等極端天氣對模擬結(jié)果的影響顯著；三是SiB2模型存在很多局限性，集合卡爾曼濾波的同化性能受模型誤差的制約較大，同時，如何確定觀測誤差、集合的大小、同化頻率等對同化結(jié)果至關重要。江西低丘紅壤區(qū)地處亞熱帶季風區(qū)水文環(huán)境復雜、氣候多變、作物生長發(fā)育迅速，因此，在應用SiB2對低丘紅壤區(qū)農(nóng)田小流域進行長期的模擬試驗時仍然有必要對其進行改進以及發(fā)展更加精確高效的區(qū)域尺度陸面同化方案。

致謝:南京信息工程大學應用氣象學院朱博鈞博士生對研究給予了幫助，謹致謝意！

參考文獻：

[1]李新,黃春林,車濤,等.中國陸面數(shù)據(jù)同化系統(tǒng)研究的進展與前瞻[J].自然科學進展,2007,17(2):163-173.

[2]高志球,卞林根,程彥杰,等.利用生物圈模型(SiB2)模擬青藏高原那曲草原近地面層能量收支[J].應用氣象學報,2002,13(2):129-141.

[3]趙其國,黃國勤,馬艷芹.中國南方紅壤生態(tài)系統(tǒng)面臨的問題及對策[J].生態(tài)學報,2013,33(24):7 615-7 622.

[4]景元書,張斌.亞熱帶花生旱坡地的水分動態(tài)特征分析[J].中國農(nóng)學通報,2006,22(11):467-471.

[5]張川,張偉,陳洪松,等.喀斯特典型坡地旱季表層土壤水分時空變異性[J].生態(tài)學報,2015,35(19):6 326-6 334.

[6]景元書,張斌,孫秀芬.坡地對花生氣候生產(chǎn)潛力利用率的影響[J].花生學報,2003,32(1):17-20.

[7]黃春林,李新.基于集合卡爾曼濾波的土壤水分同化試驗[J].高原氣象,2006,25(4):665-671.

[8]Sellers PJ, Randall DA, Collatz GJ, et al. A revised land surface parameterization (sib2) for atmospheric gcms.part i: model formulation[J]. Journal of Climate, 1996,9(4):676-705.

[9]Sellers J P,Tucker J C,Collatz J G, et al. A revised land surface parameterization (sib2) for atmospheric gcms. part ii: the generation of global fields of terrestrial biophysical parameters from satellite data[J]. Journal of Climate, 1996,9(4):706-737.

[10]姚鎮(zhèn)海,邱新法,施國萍,等.中國比輻射率空間分布特征分析[J].干旱區(qū)地理,2015,38(6):1 234-1 240.

[11]王藝, 朱彬, 劉煜,等. 中國地區(qū)近10年地表反照率變化趨勢[J]. 氣象科技, 2011,39(2):147-155.

[12]姚彤, 張強. 我國北方不同類型下墊面地表反照率特征[J]. 物理學報, 2014,63(8):462-470.

[13]高國棟,陸渝蓉,吳必文,等.地面長波輻射量的計算及云對其影響的初步分析[J].氣象科學,1985,(2):77-84.

[14]郭元喜, 龔道溢, 汪文珊,等. 中國東部夏季云量與日氣溫統(tǒng)計關系[J]. 地理科學, 2013,33(1):104-109.

[15]Nakajima T, Higurashi A, Kawamoto K, et al. Nakajima, T. et al. A possible correlation between satellite-derived cloud and aerosol microphysical parameters. Geophysical. Research. Letters. 28, 1171-1174[J]. Geophysical Research Letters, 2001,28(7):1 171-1 174.

[16]韓旭軍, 李新. 非線性濾波方法與陸面數(shù)據(jù)同化[J]. 地球科學進展, 2008,23(8):813-820.

[17]Anderson B D O, Moore J B, Eslami M. Optimal Filtering[J]. IEEE Transactions on Systems Man & Cybernetics, 1979,12(2):235-236.

[18]李超, 惠建忠, 唐千紅,等. EnKF優(yōu)化土壤濕度方程中參數(shù)的性能研究[J]. 自然資源學報, 2015,(12):2 131-2 140.