付曉雷，余鐘波，丁永建，蔣曉蕾，楊傳國(guó)，鞠 琴

(1.福州大學(xué)土木工程學(xué)院，福州 350116;2.中國(guó)科學(xué)院西北生態(tài)環(huán)境資源研究院冰凍圈科學(xué)國(guó)家重點(diǎn)實(shí)驗(yàn)室，蘭州 730000；3.水文水資源與水利工程科學(xué)國(guó)家重點(diǎn)實(shí)驗(yàn)室，南京 210098;4.河海大學(xué)水文水資源學(xué)院，南京 210098)

0 引 言

土壤濕度是氣象、水文、農(nóng)業(yè)等領(lǐng)域重要的物理狀態(tài)變量，它對(duì)地表能量[1]以及降雨下滲[2,3]的再分配和水資源管理、農(nóng)業(yè)灌溉等[4,5]方面具有重要影響。盡管土壤濕度在不同學(xué)科領(lǐng)域的重要性越來(lái)越受到關(guān)注，但是由于缺乏長(zhǎng)期的全球高密度的土壤濕度觀測(cè)數(shù)據(jù)，使其土壤濕度的研究受到了影響[6]。隨著衛(wèi)星遙感技術(shù)的發(fā)展，極大地提高了全球土壤濕度數(shù)據(jù)庫(kù)的空間分辨率。然而，遙感反演產(chǎn)品僅限于距地表幾厘米[7,8]，且易受土壤類型等影響，精度不高。此外，陸面水文模型可以利用氣象、土壤、植被等信息模擬土壤水分的動(dòng)態(tài)過(guò)程，且能夠得到高時(shí)空分辨率的區(qū)域土壤濕度分布數(shù)據(jù)[9]。但是，由于模型物理過(guò)程概化的不精確性(包括大氣模型和陸面模式及參數(shù)化方案)、土壤參數(shù)和初始狀態(tài)的不準(zhǔn)確性以及氣象信息、陸面模式本身的不確定性，使得土壤濕度的模擬精度受到限制。

為了滿足對(duì)高精度高分辨率土壤濕度數(shù)據(jù)的需求，綜合觀測(cè)、反演、模擬數(shù)據(jù)的優(yōu)點(diǎn)，數(shù)據(jù)同化方法逐漸在土壤濕度研究領(lǐng)域得到了應(yīng)用發(fā)展[10-15]。尤其是Kalman[16]提出針對(duì)線性系統(tǒng)的卡爾曼濾波方法之后，數(shù)據(jù)同化方法得到了快速有效的發(fā)展以及在土壤濕度研究中的應(yīng)用。Jazwinski[16]提出了擴(kuò)展卡爾曼濾波EKF，有效解決了非線性系統(tǒng)問(wèn)題。Yeh等[10]基于擴(kuò)展卡爾曼濾波評(píng)估氣候和水文變化對(duì)土壤濕度的影響；Lü等[5]利用擴(kuò)展卡爾曼濾波同化表層土壤含水量來(lái)預(yù)測(cè)根域土壤含水量。但是，擴(kuò)展卡爾曼濾波中需要建立的動(dòng)力模型切線性算子，在強(qiáng)非線性的動(dòng)力系統(tǒng)中可能并不存在。為了更好地解決非線性問(wèn)題，Evensen[18]提出了集合卡爾曼濾波(EnKF)，并取得了巨大成功。Brandhorst等[19]通過(guò)處理土壤水力參數(shù)的不確定性，結(jié)合EnKF來(lái)預(yù)測(cè)土壤濕度；Liu等[20]利用EnKF和支持向量機(jī)分析土壤濕度模擬不確定性的因素。

然而，在利用集合卡爾曼濾波對(duì)土壤濕度進(jìn)行同化時(shí)，也存在一些不足：由于集合成員取樣的隨機(jī)性，使得每次取樣的隨機(jī)誤差不同，導(dǎo)致每次的同化結(jié)果存在一定差別；另外，每個(gè)集合成員賦予相同的權(quán)重，降低了影響性較大的集合成員的作用。本文提出了無(wú)跡加權(quán)集合卡爾曼濾波UWEnKF，增加了較大影響性集合成員的權(quán)重，并結(jié)合一維Richards方程進(jìn)行了土壤濕度同化實(shí)驗(yàn)，驗(yàn)證其有效性和適用性。

1 方 法

1.1 Richards方程

一維Richards方程描述了非飽和區(qū)域土壤水分的運(yùn)動(dòng)過(guò)程[5,21]，該方程可以表達(dá)為：

(1)

式中：θ是土壤水分含量，m3/m3；t是時(shí)間，s；z是土壤深度，mm；ψ是土水勢(shì)，mm；K是土壤水力傳導(dǎo)度，mm/s；e是源匯項(xiàng)，即根系吸水后的蒸散發(fā)損失率，mm/s。

源匯項(xiàng)的求解參見(jiàn)Lü等[5]，對(duì)方程離散求解過(guò)程中，上邊界條件取決于下滲率、雨強(qiáng)(當(dāng)雨強(qiáng)小于下滲率，按雨強(qiáng)下滲，反之按下滲率下滲)；下邊界條件采用最底層的土壤水力傳導(dǎo)度。方程的時(shí)間和空間離散分別采用Crank-Nicolson方法和向后的歐拉方法。

1.2 無(wú)跡加權(quán)集合卡爾曼濾波UWEnKF

集合卡爾曼濾波EnKF[8,11,18,21-25]于1994年由Evensen[18]提出，之后便被廣泛用于提高土壤濕度的模擬精度。該方法包含預(yù)測(cè)和分析兩個(gè)階段，在分析階段，通過(guò)引入噪聲擾動(dòng)，生成狀態(tài)變量的初始集合，利用模型算子得到每個(gè)集合成員的預(yù)測(cè)值，以及狀態(tài)變量的誤差協(xié)方差矩陣；在分析階段，利用觀測(cè)值對(duì)狀態(tài)變量的誤差協(xié)方差矩陣和預(yù)測(cè)值進(jìn)行更新，得到狀態(tài)變量的分析值。許多學(xué)者對(duì)其理論過(guò)程進(jìn)行了詳細(xì)的闡述[8,11,18,21-25]。因此，本文不在對(duì)EnKF進(jìn)行描述。

(2)

本文在EnKF的基礎(chǔ)上，結(jié)合無(wú)跡變換[26,27]理論，提出了無(wú)跡加權(quán)集合卡爾曼濾波UWEnKF，該方法增加了影響性較大的集合成員的權(quán)重，且集合成員關(guān)于均值對(duì)稱。

基于公式(2)和隨機(jī)誤差粒子，可以得到UWEnKF中對(duì)稱性的 個(gè)Sigma集合成員：

(3)

UWEnKF的預(yù)測(cè)和更新過(guò)程如下：

預(yù)測(cè):

每個(gè)集合成員的狀態(tài)預(yù)測(cè)過(guò)程為：

(4)

更新:

k時(shí)刻狀態(tài)變量X的誤差協(xié)方差矩陣的計(jì)算如下：

PXk|k-1=A(I·×WT)AT

(5)

(6)

(7)

(8)

(9)

(10)

(11)

PYk=S(I·×WT)ST

(12)

(13)

PXYk=A(I·×WT)ST

(14)

Kk=PXYk(PYk+R)-1

(15)

(16)

(17)

(18)

PXk=E(I·×WT)ET

(19)

1.3 評(píng)價(jià)準(zhǔn)則

本研究采用均方根誤差RMSE、絕對(duì)誤差A(yù)E和平均絕對(duì)誤差MAE三個(gè)標(biāo)準(zhǔn)對(duì)UWEnKF的有效性進(jìn)行評(píng)估，具體公式如下：

(20)

AE=|Xsim,i-Xobs,i|

(21)

(22)

式中：Xsim,i和Xobs,i分別是第i時(shí)刻各層土壤濕度的模擬值/同化值和觀測(cè)值；m為總模擬時(shí)刻數(shù)。

1.4 研究區(qū)和實(shí)驗(yàn)數(shù)據(jù)

流域面積約70 hm2的梅林小流域(31°20′ N,119°51′ E)屬于封閉型小流域，地處江蘇省宜興市東南，位于太湖以西約9 km處，如圖1[29,30]。流域?qū)贊駶?rùn)氣候區(qū)域，多年平均降雨量約1 150 mm，年平均氣溫15.5 ℃，土壤類型主要為紅壤和黃棕壤，植被主要為林地、茶樹、草地[30]。本文收集了該流域2007年8月1日(第244 d)至9月30日(第273 d)的實(shí)時(shí)(步長(zhǎng)為1 h)氣象數(shù)據(jù)，期間最大日降雨量65 mm，日平均氣溫23.6 ℃，平均相對(duì)濕度86.5%，平均風(fēng)速8.1 m/s。收集了流域內(nèi)楊家山水文實(shí)驗(yàn)站A站點(diǎn)2007年8月10日至22日、以及8月26日至9月30日的深度分別為5、30、50、100 cm的土壤濕度實(shí)時(shí)觀測(cè)數(shù)據(jù)等，觀測(cè)步長(zhǎng)為1 h，以及土壤孔隙率、殘差土壤含水量等參數(shù)值。根據(jù)Lü等[5]關(guān)于梅林流域不同深度的土壤屬性數(shù)值，本文采用的各層土壤參數(shù)見(jiàn)表1。

圖1 梅林流域及土壤濕度觀測(cè)站點(diǎn)[28,29]Fig.1 Meilin study area and soil moisture observation sites

表1 梅林流域各層土壤參數(shù)的設(shè)定[5]Tab.1 The values of soil parameters at different depth in Meilin area

2 實(shí)驗(yàn)結(jié)果與分析

在實(shí)驗(yàn)設(shè)置中，首先利用2007年8月1日至8月31日的實(shí)時(shí)氣象數(shù)據(jù)對(duì)一維Richards方程進(jìn)行校正，然后，利用9月1日至30日的數(shù)據(jù)進(jìn)行土壤濕度數(shù)據(jù)同化實(shí)驗(yàn)。在同化過(guò)程中，僅利用EnKF和UWEnKF同化土壤表層(5 cm)的站點(diǎn)觀測(cè)土壤濕度，進(jìn)而更新不同深度的土壤濕度模擬值。在同化實(shí)驗(yàn)中，隨機(jī)產(chǎn)生的誤差粒子數(shù)N設(shè)為80；方程模擬了4層不同土壤層的土壤濕度，所以，狀態(tài)變量的維數(shù)n為4，則Sigma集合成員數(shù)為2nN+1=641。根據(jù)氣象數(shù)據(jù)和土壤濕度數(shù)據(jù)的觀測(cè)步長(zhǎng)，同化時(shí)間步長(zhǎng)和模型模擬步長(zhǎng)均設(shè)為1 h。圖2顯示了不同深度的土壤濕度同化及模擬值，其中模擬值是基于表1的土壤參數(shù)和氣象數(shù)據(jù)等，單獨(dú)求解Richards方程得到。

由圖2可知，不同深度(5、30、50、100 cm)的土壤濕度變化趨勢(shì)基本一致，都在降雨量較大的第246 d(9月3日)和第262 d(9月19日)有明顯增高，其中，30 cm深的土壤濕度變化最為劇烈，100 cm深的土壤濕度變化最為緩慢。校正后的一維Richards方程在5、30、50 cm土壤層的模擬結(jié)果與實(shí)測(cè)變化趨勢(shì)一致，在強(qiáng)降雨發(fā)生時(shí)(第246 d和第262 d)，一維Richards方程可以反映土壤濕度的突增變化，但對(duì)增長(zhǎng)幅度的模擬精度不夠。盡管如此，經(jīng)過(guò)一段時(shí)間(第248～261 d，第263～273 d)的穩(wěn)定期(無(wú)強(qiáng)降雨)后，土壤濕度的模擬結(jié)果與觀測(cè)趨于一致。在底層(100 cm)，一維Richards方程的模擬結(jié)果無(wú)法反映土壤濕度的變化過(guò)程，尤其是在降雨量較大的第262 d(9月19日)，原因可能是下邊界采用最底層的土壤水力傳導(dǎo)度，而非實(shí)際下邊界的土壤水力傳導(dǎo)度，導(dǎo)致土壤濕度模擬結(jié)果與實(shí)際情況存在偏差。

圖2還表明，經(jīng)過(guò)同化后(EnKF和UWEnKF)的土壤濕度比模擬值更接近觀測(cè)值，而UWEnKF方法的同化效果比EnKF方法更好。在上述4種不同深度的土壤層中，UWEnKF同化值比EnKF同化值更趨近實(shí)測(cè)土壤濕度。當(dāng)有強(qiáng)降雨輸入時(shí)(第246 d和第262 d)，雖然Richards方程能夠較好地反映土壤濕度的突然升高變化，但是，土壤濕度模擬值與觀測(cè)值之間相差較大，導(dǎo)致同化值與觀測(cè)值之間的誤差也較大。在第264～271 d的第3層(50 cm)，同化值遠(yuǎn)離觀測(cè)值，表明在此期間，同化方法在該層沒(méi)有改善土壤濕度模擬結(jié)果。主要原因可能在于經(jīng)過(guò)9月19日的強(qiáng)降雨后，土壤濕度增量較大，模擬結(jié)果雖然反映出土壤濕度變化過(guò)程，但是增量幅度較小(第3層)，而EnKF和UWEnKF使得同化結(jié)果大于模擬結(jié)果，以至于在264～273 dEnKF和UWEnKF同化后的結(jié)果在此期間大于模擬結(jié)果，但是，UWEnKF同化后的結(jié)果與模擬結(jié)果越來(lái)越接近。在底層(100 cm)，由于Richards方程的模擬結(jié)果較穩(wěn)定，導(dǎo)致同化結(jié)果也趨于穩(wěn)定，但是，整體上UWEnKF同化值比EnKF同化值更接近觀測(cè)值。綜上所述：在土壤濕度同化中，UWEnKF方法比EnKF方法更有效。

圖2 2007年9月1日至30日(第273 d)梅林實(shí)驗(yàn)區(qū)域觀測(cè)站點(diǎn)A的土壤濕度模擬/同化結(jié)果Fig.2 The soil moisture simulations and assimilations from Sep.1 (day 244)to Sep.30 (day 273),2007 at Point A in Meilin watershed

本文計(jì)算了一維Richards方程、EnKF方法和UWEnKF方法在不同土壤層的均方根誤差RMSE、絕對(duì)誤差A(yù)E和平均絕對(duì)誤差MAE來(lái)評(píng)價(jià)方法對(duì)土壤濕度的模擬/同化效果，圖3顯示了模擬/同化結(jié)果與觀測(cè)值之間的均方根誤差RMSE，圖4和圖5分別顯示了土壤濕度模擬/同化值與觀測(cè)值之間的絕對(duì)誤差A(yù)E變化過(guò)程和平均絕對(duì)誤差MAE。

由圖3可知，EnKF和UWEnKF同化的結(jié)果與觀測(cè)值之間的RMSE在前3層要明顯小于模擬結(jié)果與觀測(cè)值之間的RMSE，說(shuō)明EnKF和UWEnKF可以改善土壤濕度的模擬結(jié)果，提高土壤濕度的模擬精度。然而，在底層，EnKF同化后的結(jié)果與觀測(cè)值之間的RMSE與模擬結(jié)果與觀測(cè)值之間的RMSE相差不大，表明EnKF僅同化表層土壤含水量時(shí)，在底層改善土壤濕度模擬結(jié)果的效果不明顯；UWEnKF在該層仍可顯著提高土壤濕度的模擬精度。另外，UWEnKF方法的RMSE值在各土壤層均明顯小于EnKF方法，表明UWEnKF方法比EnKF方法更好地提高了土壤濕度的模擬精度。這是因?yàn)镋nKF方法中集合成員等權(quán)重，而在UWEnKF方法中，影響性大的集合成員被賦予了更高的權(quán)重，這提高了濾波方法的有效性。

由圖4可知，在上述4個(gè)土壤層，同化值(EnKF和UWEnKF)的AE整體上而言均小于模擬值，且模擬結(jié)果與觀測(cè)值之間的AE越大，基本上同化結(jié)果與模擬結(jié)果之間的AE也越大，即模型誤差越大，同化結(jié)果的誤差越大。但在第3層(50 cm)，EnKF和UWEnKF同化值的AE在264～273 d要大于模擬值的AE；在第四層(100 cm)，EnKF同化值的AE在第268～273 d高于模擬值的AE(見(jiàn)前文分析)。另外，UWEnKF同化結(jié)果的AE整體上小于EnKF同化結(jié)果的AE。其中在第2～4層，除了降雨量較大的前后幾天，相比EnKF，UWEnKF同化結(jié)果的AE明顯小于0.02 m3/m3，甚至小于0.01 m3/m3。因此，UWEnKF方法比EnKF方法更好地改善了土壤濕度的模擬結(jié)果。

圖3 2007年9月1日(第244 d)至30日梅林實(shí)驗(yàn)區(qū)域觀測(cè)站點(diǎn)A不同深度的土壤濕度同化/模擬結(jié)果與觀測(cè)值之間的均方根誤差RMSEFig.3 The RMSE between the soil moisture assimilations/simulations and observations Sep.1 to Sep.30,2007 at Point A in Meilin watershed

圖4 2007年9月1日至30日梅林實(shí)驗(yàn)區(qū)域觀測(cè)站點(diǎn)A的土壤濕度同化/模擬結(jié)果的絕對(duì)誤差變化過(guò)程Fig.4 The AE changes of soil moisture assimilations/simulations from Sep.1 to Sep.30,2007 at Point A in Meilin watershed

由圖5可知，同化和模擬的MAE類似于圖3中RMSE的結(jié)論，即EnKF和UWEnKF同化的MAE在各層要明顯小于模擬結(jié)果的MAE，且UWEnKF方法的MAE值在各土壤層均明顯小于EnKF方法的MAE。

圖5 2007年9月1日至30日梅林實(shí)驗(yàn)區(qū)域觀測(cè)站點(diǎn)A不同深度的土壤濕度同化/模擬結(jié)果與觀測(cè)值之間的平均絕對(duì)誤差MAEFig.5 The MAE between the soil moisture assimilations/simulations and observations Sep.1 to Sep.30,2007 at Point A in Meilin watershed

3 結(jié) 論

本文基于無(wú)跡轉(zhuǎn)移方程和集合卡爾曼濾波，提出了無(wú)跡加權(quán)集合卡爾曼濾波UWEnKF。該濾波方法中集合成員的影響性越大，其相應(yīng)的權(quán)重越大，且集合成員的分布關(guān)于均值對(duì)稱。結(jié)合一維Richards方程在梅林流域開展了同化表層土壤濕度的實(shí)時(shí)同化實(shí)驗(yàn)。具體的結(jié)論如下：

(1)一維Richards方程可以較好地模擬出土壤濕度的動(dòng)態(tài)變化過(guò)程，但在底層，該方程的模擬效果不佳，無(wú)法較好地反應(yīng)土壤濕度的變化過(guò)程。

(2)在土壤濕度模擬過(guò)程中，采用數(shù)據(jù)同化方法(UWEnKF和EnKF)對(duì)表層(5 cm)土壤濕度進(jìn)行同化，可以提高不同土壤層的土壤濕度模擬精度。同時(shí)，同化效果與模型聯(lián)系緊密，模型誤差越大，同化結(jié)果的誤差越大。

(3)通過(guò)土壤濕度模擬的同化實(shí)驗(yàn)可知，相較于EnKF方法，無(wú)跡加權(quán)集合卡爾曼濾波UWEnKF能夠更顯著地改善土壤濕度的模擬結(jié)果，提高土壤濕度的模擬精度。

總之，無(wú)跡加權(quán)集合卡爾曼濾波UWEnKF比EnKF能夠更顯著提高土壤水分模擬精度，是一種有效的、實(shí)用的數(shù)據(jù)同化方法。