基于集合卡爾曼濾波和生物圈模型的陸面數(shù)據(jù)同化

基于集合卡爾曼濾波和生物圈模型的陸面數(shù)據(jù)同化

1對(duì)土壤水分的數(shù)據(jù)同化64數(shù)據(jù)的同源性特征（cdda，foys，簡(jiǎn)稱(chēng)數(shù)據(jù)同源性）是指在考慮數(shù)據(jù)的時(shí)空分布、觀測(cè)場(chǎng)地和背景場(chǎng)誤差的基礎(chǔ)上，在數(shù)值模擬的動(dòng)態(tài)執(zhí)行過(guò)程中整合新的觀測(cè)數(shù)據(jù)的方法。陸面數(shù)據(jù)同化系統(tǒng)(LDAS,LandDataAssimilationSystem)是近年來(lái)將四維同化方法應(yīng)用到地球表層科學(xué)和水文學(xué)中而迅速發(fā)展起來(lái)的,它作為獨(dú)立的領(lǐng)域出現(xiàn)是在1995年之后。當(dāng)前,陸面數(shù)據(jù)同化的研究主要對(duì)給定的陸面模型和水文模型,采用不同的數(shù)據(jù)同化算法同化地表觀測(cè)資料、衛(wèi)星和雷達(dá)數(shù)據(jù),優(yōu)化地表和根區(qū)土壤水分、溫度、地表能量通量等的估算。李新等\,Pathmathevan等利用模擬退火算法對(duì)陸面過(guò)程模型(SiB2)進(jìn)行了土壤水分和地表溫度的同化實(shí)驗(yàn);Houser等、Schuurmans等把數(shù)據(jù)同化算法與水文模型結(jié)合,進(jìn)行了土壤水分和潛熱通量的同化試驗(yàn);Galanlowicz等\,Hoeben等、Reichle等、Walker等也對(duì)土壤水分的同化方法進(jìn)行了研究和比較;Crow等、Kumar、Lakshmi進(jìn)行了地表溫度的同化方法的研究。集合卡爾曼濾波是Evensen等在1994年根據(jù)Epstein的隨機(jī)動(dòng)態(tài)預(yù)報(bào)理論提出的順序數(shù)據(jù)同化算法。它將模型狀態(tài)預(yù)報(bào)看成近似隨機(jī)動(dòng)態(tài)預(yù)報(bào),用一個(gè)狀態(tài)總體(設(shè)數(shù)目為N)去代表隨機(jī)動(dòng)態(tài)預(yù)報(bào)中的概率密度函數(shù),通過(guò)向前積分,狀態(tài)總體很容易計(jì)算不同時(shí)間的概率密度函數(shù)所對(duì)應(yīng)的統(tǒng)計(jì)特性(如均值與方差)。集合卡爾曼濾波的最大特點(diǎn)是它克服了卡爾曼濾波要求線(xiàn)性化的模型算子和觀測(cè)算子的缺點(diǎn)。經(jīng)過(guò)近十年的研究和發(fā)展,集合卡爾曼濾波算法逐漸發(fā)展成熟,也逐步受到數(shù)據(jù)同化研究人員的重視,不同作者[22,23,24,25,26,27,28]對(duì)集合卡爾曼濾波算法的理論和應(yīng)用進(jìn)行了廣泛探討。土壤水分是地氣相互作用、水文循環(huán)等研究的關(guān)鍵變量,它影響地表能量通量、徑流、輻射平衡、物質(zhì)遷移等。土壤水分的準(zhǔn)確估計(jì)對(duì)于研究和理解地球表層生物物理過(guò)程起著重要作用。土壤水分在空間和時(shí)間上變化非常大,建立稠密的地球表層變量的觀測(cè)網(wǎng)幾乎不現(xiàn)實(shí)。此外,陸面過(guò)程或水文模型雖然可以模擬土壤水分的連續(xù)變化,但隨著時(shí)間的變化,誤差不斷累積,導(dǎo)致模擬結(jié)果較差。通過(guò)建立陸面數(shù)據(jù)同化系統(tǒng)可以把相關(guān)的多源信息(地表觀測(cè)、衛(wèi)星、雷達(dá)等)融入陸面過(guò)程模型或水文模型,在有觀測(cè)的時(shí)間段校正模型對(duì)土壤水分的預(yù)報(bào)值,降低誤差的累積,從而可以提高對(duì)土壤水分的估算精度。2方法陸面數(shù)據(jù)同化系統(tǒng)主要由數(shù)據(jù)同化算法、陸面過(guò)程模型、數(shù)據(jù)(驅(qū)動(dòng)數(shù)據(jù)、參數(shù)、觀測(cè)數(shù)據(jù)、輸出數(shù)據(jù))等構(gòu)成(圖1)。2.1層土壤的水分控制方程陸面過(guò)程是指發(fā)生在地表和大氣之間水分、熱量和動(dòng)量交換的作用過(guò)程,包括地面上的熱力過(guò)程、水文過(guò)程和生物過(guò)程,地氣間的能量和物質(zhì)交換以及地面以下土壤中的熱傳導(dǎo)和水熱輸送過(guò)程等。本文采用了發(fā)展比較成熟的SiB2模型模擬土壤水分的變化。在SiB2中定義了表層、根區(qū)和深層3層土壤。公式(1)~(3)為每層土壤的水分控制方程。?θ1?t=1D1[Q1-Q12-Egρw],(1)?θ2?t=1D2[Q12-Q23-Ectρw],(2)?θ3?t=1D3[Q23-Q3],(3)?θ1?t=1D1[Q1?Q12?Egρw],(1)?θ2?t=1D2[Q12?Q23?Ectρw],(2)?θ3?t=1D3[Q23?Q3],(3)式中θ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為上邊界進(jìn)入土壤表層的水流(m/s);Q3為重力排水(m/s);Eg,Ect分別是土壤表層的蒸發(fā)率(m/s)和通過(guò)氣孔的植被蒸騰率(m/s)。2.2數(shù)據(jù)同化算法2.2.1卡爾曼濾波的計(jì)算步驟1960年,Kalman等針對(duì)隨機(jī)過(guò)程狀態(tài)估計(jì)提出Kalman濾波的思想,它包括預(yù)報(bào)和分析兩個(gè)步驟。在預(yù)報(bào)階段,根據(jù)前一時(shí)刻的模式狀態(tài)生成當(dāng)前時(shí)刻模式狀態(tài)的預(yù)報(bào)值。在分析階段,引入觀測(cè)數(shù)據(jù),利用最小方差估計(jì)方法對(duì)模式狀態(tài)進(jìn)行重新分析。隨著模式狀態(tài)預(yù)報(bào)的繼續(xù)進(jìn)行和新的觀測(cè)數(shù)據(jù)的陸續(xù)輸入,這個(gè)過(guò)程不斷向前推進(jìn)。在卡爾曼濾波算法中,認(rèn)為模型誤差方差矩陣(Q)和觀測(cè)誤差方差矩陣(R)是已知的?？柭鼮V波的計(jì)算步驟如下:(1)設(shè)定初始時(shí)刻t0的狀態(tài)變量Xa0a0和它的誤差方差矩陣Pa0a0。(2)根據(jù)第k時(shí)刻(k=0,1,2,…,)的分析場(chǎng)的狀態(tài)變量Xakak和背景場(chǎng)誤差方差矩陣Pakak計(jì)算第k+1時(shí)刻預(yù)報(bào)場(chǎng)的狀態(tài)變量Xfk+1fk+1和背景場(chǎng)誤差方差矩陣Pfk+1fk+1。Xfk+1=ΜXak,(4)Ρfk+1=ΜΡakΜΤ+Q,(5)Xfk+1=MXak,(4)Pfk+1=MPakMT+Q,(5)式中M為線(xiàn)性的模型算子,MT為模型算子的轉(zhuǎn)置,Q為模型誤差方差矩陣。(3)計(jì)算k+1時(shí)刻的卡爾曼增益矩陣Kk+1。Κk+1=Ρfk+1ΗΤ(ΗΡfk+1ΗΤ+R)-1?(6)式中H為線(xiàn)性的觀測(cè)算子,HT為觀測(cè)算子的轉(zhuǎn)置,R為觀測(cè)誤差方差矩陣。(4)計(jì)算k+1時(shí)刻分析場(chǎng)的狀態(tài)變量Xak+1和它的誤差方差矩陣Pak+1。Xak+1=Xfk+1+Κk+1(Ζk+1-ΗXfk+1),(7)Ρak+1=(Ι-Κk+1Η)Ρfk+1,(8)式中Zk+1為k+1時(shí)刻的觀測(cè)值,I為單位矩陣。(5)進(jìn)入下一時(shí)刻,再返回步驟(2)。如此循環(huán)往復(fù)地執(zhí)行預(yù)報(bào)和分析步驟。從卡爾曼濾波的計(jì)算公式我們可以看出,在計(jì)算過(guò)程中該算法需要線(xiàn)性的模型算子和觀測(cè)算子。對(duì)于復(fù)雜的非線(xiàn)性系統(tǒng)(特別是地球系統(tǒng)科學(xué)),給出線(xiàn)性的模型算子是非常困難的。模型算子和觀測(cè)算子線(xiàn)性化問(wèn)題成為制約卡爾曼濾波廣泛應(yīng)用的瓶頸。2.2.2模型狀態(tài)變量的計(jì)算集合卡爾曼濾波是Evensen在1994年根據(jù)Epstein的隨機(jī)動(dòng)態(tài)預(yù)報(bào)理論提出的順序數(shù)據(jù)同化算法。它利用了蒙特卡羅方法的思想,用符合高斯分布的一組隨機(jī)變量(設(shè)數(shù)目為N)去代表隨機(jī)動(dòng)態(tài)預(yù)報(bào)中的概率密度函數(shù),通過(guò)向前積分,計(jì)算下一時(shí)刻狀態(tài)總體的概率密度函數(shù),并得到該時(shí)刻的統(tǒng)計(jì)特性(如均值與協(xié)方差)。在本研究中,使用Houtekamer提出的算法。集合卡爾曼濾波的計(jì)算步驟如下:(1)初始化背景場(chǎng)。給定N個(gè)符合高斯分布的隨機(jī)變量Xi(i=1,…,N),在本文中,狀態(tài)變量為地表、根區(qū)、深層土壤體積含水量,Xi=[θ1,θ2,θ3]Τi。(2)計(jì)算每個(gè)隨機(jī)變量在第k+1時(shí)刻的預(yù)報(bào)值Xfi,k+1。Xfi,k+1=Μk(Xai,k)+wi,kwi,k~Ν(0,Qk),(9)式中Xfi,k+1為隨機(jī)變量i在k+1時(shí)刻的預(yù)報(bào)值,Xai,k為隨機(jī)變量i在k時(shí)刻的分析值,Mk(·)為非線(xiàn)性的模型算子,在本文中為SiB2模型;Qk為模型誤差方差矩陣,wi,k為期望為0,方差為Qk的高斯白噪聲。(3)計(jì)算k+1時(shí)刻的卡爾曼增益矩陣Kk+1。Kk+1=Pfk+1HT(HPfk+1HT+Rk)-1,(10)Ρfk+1=1Ν-1Ν∑i=1(Xfi,k+1-ˉXfk+1)?(Xfi,k+1-ˉXfk+1)Τ(11)Ρfk+1ΗΤ=1Ν-1Ν∑i=1[Xfi,k+1-ˉXfk+1][Η(Xfi,k+1)-Η(ˉXfk+1)]Τ,(12)ΗΡfk+1ΗΤ=1Ν-1Ν∑i=1[Η(Xfi,k+1)-Η(ˉXfk+1]?[Η(Xfi,k+1)-Η(ˉXfk+1)]Τ,(13)ˉXfk+1=1ΝΝ∑i=1Xfi,k+1,(14)式中Xfi,k+1為第i個(gè)隨機(jī)變量在k+1時(shí)刻的預(yù)報(bào)值,ˉXfk+1為k+1時(shí)刻狀態(tài)變量預(yù)報(bào)值的平均值,Rk為觀測(cè)誤差方差矩陣,Hk(·)為觀測(cè)算子,本文中的狀態(tài)變量和觀測(cè)均為3層土壤體積含水量,所以Hk(·)為3×3的單位方陣。(4)計(jì)算k+1時(shí)刻分析場(chǎng)的狀態(tài)變量平均值Xak+1和背景場(chǎng)誤差方差矩陣Pak+1。Xai,k+1=Xfi,k+1+Kk+1[yk+1-H(xfi,k+1)+vi,k+1]vi,k~N(0,Rk),(15)ˉXak+1=1ΝΝ∑i=1Xai,k+1,(16)Ρak+1=1Ν-1Ν∑i=1(Xai,k+1-ˉXak+1)(Xai,k+1-ˉXak+1)Τ,(17)式中Xai,k+1為第i個(gè)隨機(jī)變量在k+1時(shí)刻的分析值,ˉXak+1為k+1時(shí)刻狀態(tài)變量分析值的平均值,vi,k為期望為0,方差為Rk的高斯白噪聲。(5)進(jìn)入下一時(shí)刻,返回步驟(2)。從公式(10)~(15)中可以看出,集合卡爾曼濾波在計(jì)算卡爾曼增益矩陣和k+1時(shí)刻分析場(chǎng)的狀態(tài)變量時(shí),利用蒙特卡羅方法有效地避免了卡爾曼濾波需要線(xiàn)性化的模型算子和觀測(cè)算子的難題。2.3驅(qū)動(dòng)數(shù)據(jù)和觀測(cè)數(shù)據(jù)“全球能量與水循環(huán)亞洲季風(fēng)之青藏高原試驗(yàn)”(GAME-Tibet)已經(jīng)積累了大量觀測(cè)資料并且已被用于青藏高原的水熱平衡模擬研究,本文所需的驅(qū)動(dòng)數(shù)據(jù)和觀測(cè)數(shù)據(jù)均取自GAME-Tibet試驗(yàn)區(qū)MS3608觀測(cè)站1998年7月6日至8月9日的觀測(cè)數(shù)據(jù),時(shí)間分辨率為1h。該站點(diǎn)位于31°13.6′N(xiāo),91°47.0′E,海拔高度4610m,地表覆蓋類(lèi)型為矮草。2.3.1輻射和降水條件SiB2模型所需的驅(qū)動(dòng)數(shù)據(jù)包括短波輻射、氣溫、風(fēng)速、水汽壓、長(zhǎng)波輻射和降水。其中短波輻射、氣溫、風(fēng)速資料直接來(lái)自于自動(dòng)氣象站(AWS)觀測(cè),長(zhǎng)波輻射使用SiB2中的標(biāo)準(zhǔn)方法-Brunt公式計(jì)算,水汽壓由濕度換算,降水資料取自自動(dòng)記錄的雨量筒。2.3.2errs等和李新等的定義青藏高原MS3608站點(diǎn)的土壤、植被等參數(shù)主要參照Sellers等和李新等的定義。不同土壤的飽和水力傳導(dǎo)系數(shù)變化范圍很大(1.2×10-4~6×10-6m/s),對(duì)于模擬結(jié)果影響很大,本文采用6×10-5m/s。2.3.3深度4.10,5.1土壤水分?jǐn)?shù)據(jù)由時(shí)域反射儀(TDR)測(cè)量,深度分別為4,20,60,100,160和196cm。在同化實(shí)驗(yàn)中只利用了4,20,100cm的土壤水分觀測(cè)數(shù)據(jù)。2.4otmoanstraft安氏原理對(duì)于模擬、同化、觀測(cè)結(jié)果,我們采用了兩種誤差分析方法:(1)均方根誤差(RMSE,RootMeanSquareError);(2)平均誤差(AVE,AverageError)。RΜSE=√1ΝΝ∑t=1(Οbst-Xt)2,(18)AVE=1ΝΝ∑t=1(Οbst-Xt),(19)式中N為整個(gè)觀測(cè)期,Obst為t時(shí)刻的觀測(cè)值,Xt為t時(shí)刻的SiB2模擬值或同化結(jié)果。3結(jié)果分析3.1亞風(fēng)場(chǎng)模擬結(jié)果有6個(gè)月的情況,自然集合數(shù)取48,利用1998年7月6日、7日逐時(shí)4,20,100cm的土壤水分?jǐn)?shù)據(jù)作為初始狀態(tài)變量的集合,同化間隔為6h(利用00,06,12和18時(shí)的4,20,100cm土壤體積含水量觀測(cè)數(shù)據(jù)同化SiB2模型輸出的表層、根區(qū)、深層土壤水分),模型誤差方差矩陣(Q)和觀測(cè)場(chǎng)誤差方差矩陣(R)分別為Qi=[2.5×10-5,0,00,1.44×10-5,00,0,2.5×10-5],Rk=[1×10-4,0,00,4×10-4,00,0,1×10-4].圖2a~c是1998年7月6日至8月10日亞洲季風(fēng)實(shí)驗(yàn)青藏高原實(shí)驗(yàn)區(qū)MS3608站點(diǎn)表層、根區(qū)、深層土壤水分同化結(jié)果。圖2a是表層土壤水分的同化結(jié)果。由圖可看到,降水對(duì)于表層土壤水分的影響很大。當(dāng)有降水時(shí),SiB2模擬結(jié)果與4cm觀測(cè)值相比有些偏大;無(wú)降水時(shí),SiB2模擬結(jié)果接近觀測(cè)值。對(duì)于同化結(jié)果,有降水時(shí),同化降低了SiB2對(duì)土壤水分過(guò)高的估計(jì),使結(jié)果接近觀測(cè)值;無(wú)降水時(shí),同化結(jié)果與觀測(cè)相比略偏低,在部分時(shí)段同化效果比SiB2直接模擬的效果還要差。在整體上,同化還是取得了較好的效果,但擾動(dòng)較大。圖2b是根區(qū)土壤水分的同化結(jié)果。SiB2對(duì)根區(qū)土壤水分的模擬結(jié)果整體偏大,降水對(duì)根區(qū)土壤水分的影響仍然比較大。有大的降水時(shí),模擬的根區(qū)土壤水分急劇增大。在根區(qū),同化對(duì)土壤水分的估計(jì)有比較明顯的效果。圖2c是深層土壤水分的同化結(jié)果。降水對(duì)深層土壤水分影響很小,土壤水分的變化幅度不大。經(jīng)歷一次較強(qiáng)的降水后,深層土壤水分增大的幅度較小。SiB2模擬結(jié)果略低于100cm的觀測(cè)值,同化效果顯著,同化結(jié)果幾乎和觀測(cè)值重合。3.2同化土壤水分模型圖3為表層、根區(qū)、深層土壤水分同化后誤差方差隨時(shí)間的變化。表層、根區(qū)、深層土壤水分方差都比較穩(wěn)定,在有觀測(cè)時(shí),三層土壤水分的誤差明顯降低,所以同化觀測(cè)可以消除誤差的累積效應(yīng),抑制誤差的增長(zhǎng)。表層、根區(qū)、深層誤差方差分別為<5×10-5,1.5×10-4,2.5×10-5。由表1可以看到,在SiB2模型運(yùn)行過(guò)程中同化土壤水分的觀測(cè)值降低了土壤水分模擬值的誤差。表層、根區(qū)、深層土壤水分模擬值的均方根誤差約分別降低為原來(lái)的1/2及1/9,平均誤差也分別降低為原來(lái)的1/4及1/17。在同化試驗(yàn)中,考慮到土壤水分觀測(cè)數(shù)據(jù)的代表性問(wèn)題,觀測(cè)誤差給定的較大,所以同化前后相比,表層、根區(qū)的誤差減少幅度不大。4應(yīng)用結(jié)果分析我們發(fā)展了基于集合卡爾曼濾波和SiB2的試驗(yàn)型單點(diǎn)陸面數(shù)據(jù)同化方案,利用GAME-Tibet土壤水分觀測(cè)數(shù)據(jù)進(jìn)行了同化試驗(yàn),得到以下結(jié)論:(1)基于集合卡爾曼濾波的數(shù)據(jù)同化方法可以降低模擬的誤差,同化后和同化前相比,表層、根區(qū)、深層土壤水分模擬值的均方根誤差分別降低為原來(lái)的1