無跡加權集合卡爾曼濾波及其在土壤濕度同化中的應用

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

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

0 引 言

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

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

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

1 方 法

1.1 Richards方程

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

(1)

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

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

1.2 無跡加權集合卡爾曼濾波UWEnKF

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

(2)

本文在EnKF的基礎上，結合無跡變換[26,27]理論，提出了無跡加權集合卡爾曼濾波UWEnKF，該方法增加了影響性較大的集合成員的權重，且集合成員關于均值對稱。

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

(3)

UWEnKF的預測和更新過程如下：

預測:

每個集合成員的狀態(tài)預測過程為：

(4)

更新:

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

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 評價準則

本研究采用均方根誤差RMSE、絕對誤差AE和平均絕對誤差MAE三個標準對UWEnKF的有效性進行評估，具體公式如下：

(20)

AE=|Xsim,i-Xobs,i|

(21)

(22)

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

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

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

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

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

2 實驗結果與分析

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

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

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

圖2 2007年9月1日至30日(第273 d)梅林實驗區(qū)域觀測站點A的土壤濕度模擬/同化結果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

本文計算了一維Richards方程、EnKF方法和UWEnKF方法在不同土壤層的均方根誤差RMSE、絕對誤差AE和平均絕對誤差MAE來評價方法對土壤濕度的模擬/同化效果，圖3顯示了模擬/同化結果與觀測值之間的均方根誤差RMSE，圖4和圖5分別顯示了土壤濕度模擬/同化值與觀測值之間的絕對誤差AE變化過程和平均絕對誤差MAE。

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

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

圖3 2007年9月1日(第244 d)至30日梅林實驗區(qū)域觀測站點A不同深度的土壤濕度同化/模擬結果與觀測值之間的均方根誤差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日梅林實驗區(qū)域觀測站點A的土壤濕度同化/模擬結果的絕對誤差變化過程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的結論，即EnKF和UWEnKF同化的MAE在各層要明顯小于模擬結果的MAE，且UWEnKF方法的MAE值在各土壤層均明顯小于EnKF方法的MAE。

圖5 2007年9月1日至30日梅林實驗區(qū)域觀測站點A不同深度的土壤濕度同化/模擬結果與觀測值之間的平均絕對誤差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 結 論

本文基于無跡轉移方程和集合卡爾曼濾波，提出了無跡加權集合卡爾曼濾波UWEnKF。該濾波方法中集合成員的影響性越大，其相應的權重越大，且集合成員的分布關于均值對稱。結合一維Richards方程在梅林流域開展了同化表層土壤濕度的實時同化實驗。具體的結論如下：

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

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

(3)通過土壤濕度模擬的同化實驗可知，相較于EnKF方法，無跡加權集合卡爾曼濾波UWEnKF能夠更顯著地改善土壤濕度的模擬結果，提高土壤濕度的模擬精度。

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