基于集合卡爾曼濾波法的二維土壤水流狀態變量和參數聯合估計

劉 琨，黃冠華

（1.中國農業大學 水利與土木工程學院，北京 100083；2.中國-以色列國際農業研究培訓中心，北京 100083）

1 研究背景

土壤水流模型的預測精度很大程度上依賴于土壤水力特性參數的準確估計。常用的模型參數反求方法有Levenberg-Marquardt算法［1］，Shuffled complex evolution算法［2］等，它們均以觀測值與模擬值差異最小為目標，得到一組最優的參數為估計值。然而這些反求優化算法將預測結果的不確定性歸因于模型參數，可能導致參數的次優解，進而降低了模型的預測精度。

集合卡爾曼濾波（EnKF）方法［3］可以顯式地考慮模型輸入、輸出以及模型結構等因素的不確定性，與傳統的卡爾曼濾波相比，EnKF方法引入了集合預報技術，通過蒙特卡洛方法計算預報誤差的協方差，因此可以處理高維和非線性的問題。此外EnKF方法對觀測數據集進行時間序貫處理，降低了數據的存儲量和計算量。因此，EnKF方法被廣泛應用于水文模型參數估計的研究中［4-6］。

非飽和帶水力參數與狀態變量之間的高度非線性關系增加了其參數估計的難度，近年來基于En-KF方法已開展了一些非飽和帶水力模型的參數估計研究。Brandhorst等［7］研究了在模型參數高度不確定且不可識別的條件下非飽和區模型的集合卡爾曼濾波系統的數據同化的潛力。Li和Ren［8］采用En-KF方法研究了不同土壤條件下狀態變量同化以及非均勻土壤條件下的多參數估計問題，并分析了樣本數量和觀測誤差等因素對同化效果的影響。Song等［9］評估了三種迭代集合卡爾曼濾波算法的同化效果，分析了不同觀測數據類型，觀測誤差方差，集合大小和阻尼因子等因素對同化效果的影響。

目前關于集合卡爾曼濾波法的非飽和帶參數估計的研究大多考慮一維問題［10-12］，然而在實際問題中由于復雜的邊界條件（例如溝灌、滴灌等）以及土壤剖面非均質性，土壤水流運動多為二維問題。此外與一維問題相比，二維土壤水流模型模擬區域離散后通常產生較多的計算節點，而在實際中考慮成本問題，布置觀測點的數量有限，因此如何合理的布置觀測網以實現觀測點信息的高效利用具有重要意義。總之，EnKF方法在二維土壤水流運動條件下的土壤水力特性參數估計還需進一步研究。本研究基于EnKF方法，開展二維土壤水流模型狀態和參數聯合估計研究，探討在線源入滲條件下EnKF方法對粉壤土、壤土和砂壤土的飽和導水率和進氣值參數估計精度以及對土壤壓力水頭的同化效果，分析了觀測點布置方式和觀測點數量對同化效果的影響，以期指導田間觀測網點的合理布置。

2 理論與方法

2.1 數據同化系統數據同化系統由模型算子、觀測算子和同化算法組成。首先模型算子根據上一時刻的狀態向量向前預測得到狀態向量的預測值；觀測算子將狀態向量轉化為觀測值所對應的模型預測值；最后同化系統利用觀測信息，對預測值進行同化運算，得到狀態向量的分析值。分析值作為下一時刻的模型算子輸入變量，重復以上過程直至模型模擬結束。

2.1.1 模型算子 本研究中的模型算子為二維土壤水流運動模型，即二維Richards方程:

式中:θ為體積含水率，cm3/cm3；h為壓力水頭，cm；K為非飽和導水率，cm/min；KijA為無量綱各向異性張量KA的分量；KizA為垂向方向無量綱張量KA的分量；xi和xj是空間坐標，cm；S為源匯項。非飽和導水率函數采用van Genuchten-Mualem模型表示［13-14］。本研究中使用HYDRUS-2D［15］的水流運動模塊模擬二維土壤水流運動。由于本研究觀測數據和模型預測值均為壓力水頭，因此觀測算子為對應觀測點位置處為1，其余元素為0的矩陣。

2.1.2 同化算法 采用擴展狀態變量法將模型狀態變量和土壤水力特性參數包含在一個矩陣中形成狀態向量，使用EnKF方法同時更新狀態變量和待估計參數。EnKF方法中通過向模型狀態變量和參數加入高斯白噪聲生成狀態向量集合，利用集合計算出模型預報誤差協方差矩陣，再基于蒙特卡洛抽樣方法來估計高維非線性動力學模型中的狀態變量和參數。在具體的實施過程中，EnKF方法包括預測和分析兩個步驟，首先根據前一時刻的模型狀態向量生成當前時刻狀態向量的預報值:式中:Xi,tf+1為t+1時刻第i個集合成員模型向量的預測值；Xi,at為t時刻第i個集合成員模型向量的分析值；M為模型算子，即二維土壤水流運動模型；ei為模型誤差向量，其滿足均值為0、方差為Q的高斯分布。

當有觀測數據時，利用觀測數據計算得到當前時刻狀態向量的分析值:

式中:H為觀測算子；Pef為預測值觀測誤差協方差；Re為觀測誤差協方差；di，t+1為t+1時刻各集合成員的觀測向量。詳細的EnKF計算流程請參見文獻［3，16］。

基于EnKF土壤水力參數估計流程如下:（1）模型初始設置，包括模型參數，初始條件和邊界條件設置；（2）根據集合數量通過隨機數函數抽樣得到待估計參數集合；（3）每一時刻，HYDRUS模型向前計算，得到模型向量的預測值；（4）當該時刻有觀測數據時，對觀測值進行擾動得到觀測向量；（5）利用公式（3）計算得到狀態向量的分析值，即更新壓力水頭和水力參數，作為下一時刻模型預測的初始條件；重復步驟（3）—（5）直至模擬時間結束。

2.2 數值試驗采用數值實驗的方法探究EnKF方法的同化效果。設置人工假想算例，實驗方案如下:地表設置長度為20cm的線形入滲源，模擬區域為50 cm×50 cm的均質土壤剖面?？紤]粉壤土、壤土和砂壤土三種土壤類型，三種土壤的土壤水力特性參數根據Carsel和Parrish取值［17］（見表1）。粉壤土和壤土的模擬時間設為400 min；由于砂壤土的導水率較大，入滲時間較短，模擬時間設為200 min。模型的初始條件為均勻分布的土壤壓力水頭剖面，水頭值為-100cm；上邊界條件為定水頭邊界，壓力水頭為2 cm；下邊界條件為自由排水邊界。

表1 土壤類型及其土壤水力特性參數

待估計參數選取土壤水力特性參數中變異性最大的飽和導水率和進氣值參數［17］。為了探究EnKF方法的參數估計效果，參數的初值設置均大于參數的真值，假定待估計參數服從正態分布，初始統計特征為lnKS～N（1.5a，0.01）、α～N（1.2a，0.0001）。其中，a為參數的真值。由于參數KS的方差較大，為保證KS均為正值，因此對KS進行了取對數處理。觀測數據通過以參數真值模擬水分入滲過程，每隔20 min記錄觀測點的壓力水頭，添加一定的擾動作為具有觀測誤差的觀測值，假設擾動服從正態分布N（0，0.001）。研究中樣本數目設置為150以保證同化效果。

2.3 觀測點布置方案為了研究觀測點位置對同化效果的影響，設置了不同觀測點布置方案，觀測點布置見圖1。其中方案1—4中觀測點垂向布置，分別距離模擬區域中心軸0、5、10和15 cm；方案5—8中觀測點水平向布置，分別距離地表10、20、30和40 cm。

為研究觀測點數量對同化效果的影響，選擇每種土壤最優的三種觀測點布置方案進行組合確定最終的布置方案。采用平均相對誤差（MRE）評價參數估計效果，MRE值越小，表示參數越快收斂于參數真值。采用絕對誤差（AE）評價土壤壓力水頭分布剖面的同化效果。

式中:Ya為同化后的參數樣本均值；Y為參數真值；N為同化次數。

圖1 不同觀測點位置布置方案（實心點為觀測點，單位:cm）

3 結果分析與討論

3.1 觀測點位置影響圖2為不同觀測點布置方案飽和導水率的演化。由圖2可知粉壤土條件下，垂向布置觀測點的方案比水平向布置的方案對參數KS估計效果更好，最優的3種布置方案為1、2和3（見表2）。壤土條件下，除了方案8不能較好地估計參數KS，其他方案參數KS都可以收斂于真值，最優的3種布置方案為1、5和6。砂壤土條件下，8種布置方案均可以較好地估計KS，最優的3種布置方案為2、3和4。這可能是由于粉壤土的導水率較小，水平布置觀測點時當水分運動還未到達觀測點深度，同化系統未能利用觀測信息更新參數，參數估計偏差較大，尤其隨著觀測點埋深的增加，估計效果較差。而對于壤土和砂壤土由于其導水率較大，入滲水分可以很快到達觀測點深度，因此水平向布置方案也可以得到較好的參數估計結果。

參數KS標準差隨著同化次數增加逐漸減?。ㄒ妶D2）。垂向布置方案的KS標準差變化趨勢類似，經過3次同化標準差值已經很小；水平向布置條件下，不同方案的標準差變化差異較大，隨著觀測點埋深的增加，參數估計效果變差，標準差降低的速度變小。這是由于當入滲水分沒有運動到觀測點埋深時，觀測點的壓力水頭不變，同化系統尚未利用觀測信息更新狀態變量和參數，因此參數估計結果較差。圖3為不同觀測點布置方案進氣值參數的演化，進氣值參數α的演化與KS的演化趨勢一致，KS同化效果較好的方案，對α的估計效果也較好，這是由于同化計算中參數KS和α是同時更新的。與參數KS不同，8種方案參數α的標準差均下降較快，這是由于參數α的初始方差設置較小。

圖2 不同觀測點布置方案飽和導水率演化

圖3 不同觀測點布置方案進氣值參數演化

表2 不同觀測點位置KS估計MRE

表3 不同觀測點位置α估計MRE

由于3種土壤條件下壓力水頭誤差結果類似，我們僅給出了壤土條件下不同觀測點布置方案土壤壓力水頭誤差分布（見圖4）。垂向布置條件下，方案1—4隨著觀測點到中軸距離的增加，誤差分布也隨之偏移，觀測點兩側部分的AE較大。水平布置條件下，土壤剖面0～30 cm區域的AE較大，其中方案7和8的同化效果較差。王文等［10］和Montzka等［18］也得到類似的研究結果，其研究表明使用土壤表層觀測信息，同化效果僅限于一定深度的土壤。隨著觀測點埋深增大，參數估計偏差增大，導致土壤壓力水頭的同化結果較差，并且這種參數估計偏差引起的誤差不會在同化過程中被消除［19］。以上結果表明有觀測點的區域誤差較小，說明同化系統可以有效地利用觀測信息改善土壤壓力水頭的預測結果。

圖4 壤土不同觀測點布置方案土壤壓力水頭誤差分布（單位:cm）

3.2 觀測點數量影響在觀測點位置影響結果的基礎上選擇每種土壤最優的三種布置方案組合得到新的布置方案，觀測點數量分別為4、8和12，觀測點布置見圖5。

圖6給出了不同觀測點數量條件下參數KS和α的演化。隨著觀測點數量增加，參數更快趨近于其真值。觀測點數量為12時參數的MRE值最小值，粉壤土、壤土和砂壤土KS的MRE分別為0.097、0.065和0.021，α的MRE分別為0.075、0.073和0.059。這是由于同化系統利用更多的觀測信息更新參數，提高了參數的估計精度。圖7給出了土壤壓力水頭誤差分布剖面。土壤壓力水頭誤差隨著觀測點數量的增加而減小，這是由于同化系統利用更多的觀測信息實時校正土壤壓力水頭剖面；此外增加觀測點數量可以得到更好的參數估計值，進而提高了土壤剖面壓力水頭的預測精度。這與蘭天等［20］結果一致，其研究表明增加觀測點數量可以有效地改善同化效果，參數收斂速度也更快。

本研究在數值試驗中僅考慮了3種土壤類型。這是由于質地更粗的土壤如砂土，土壤導水率較大導致入滲時間較短，當整個土壤剖面飽和時樣本之間無差異，導致在同化過程中無法進行矩陣求逆運算，導致計算過程終止；而對于細質土壤如黏土，在長時間的積水入滲條件下，如果參數的初始均值和方差以及觀測誤差方差等變量設置不合理會導致模擬計算失敗，這是由于同化過程中部分集合的土壤參數偏差較大，導致模型求解不收斂。對于細質土壤宜采用較小的參數方差以及較大的觀測值方差以保證同化過程中模型的順利運行。而本研究重點關注觀測點布置方式和數量對同化效果的影響，因此為了保證參數統計變量和觀測誤差方差一致以及模型的正常運行，選取了粉壤土、壤土和砂壤土3種土壤類型進行試驗。在之后的研究中我們將考慮更多的同化影響因素，如集合數量和同化時間間隔等，并通過設置合理的參數統計變量和觀測誤差方差等方法以涵蓋更多土壤類型開展進一步研究。

圖5 不同觀測點數量布置方案（實心點為觀測點 單位:cm）

圖6 不同觀測點數量條件下飽和導水率和進氣值參數演化

圖7 不同觀測點數目條件下土壤壓力水頭誤差分布（單位:cm）

與傳統的參數估計方法相比，EnKF方法在處理不確定性方面有著明顯的優勢。例如在Levenberg-Marquardt算法中，觀測數據被認定為是衡量預測精度的唯一標準，模型的不確定性完全歸因于參數，并沒有考慮觀測數據誤差。而在模型構建過程中對物理過程的認識不足和簡化以及觀測手段的局限，模型結構和觀測數據都具有不確定性，因此僅將不確定歸因于參數是不恰當的。EnKF方法綜合觀測數據和模型預測的不確定性，可以得到參數的優化值。本研究中我們已通過數值實驗從理論上證明了狀態變量與參數聯合估計效果，在未來的研究中還需結合田間實驗進一步檢驗該方法的精度。

4 結論

本文采用EnKF方法開展二維土壤水流運動模型狀態變量和參數聯合估計研究，采用數值實驗方法研究了粉壤土、壤土和砂壤土在線源入滲條件下EnKF方法對飽和導水率和進氣值參數估計以及壓力水頭的同化效果，分析了觀測點布置方式和觀測點數量對同化效果的影響。主要結論如下:

（1）觀測點的布置方式應捕捉二維土壤水流運動過程，才能得到較好的參數估計值。粉壤土條件下，宜采用垂向的觀測點布置方式；壤土和砂壤土條件下，在0～30 cm深土壤中水平向布置觀測點可以得到較好的參數估計值。

（2）入滲條件下觀測點布置應盡量靠近地表，使得同化系統可以盡早利用觀測信息更新狀態向量，參數可以更快地收斂于真值，但壓力水頭的同化效果僅限于一定深度的土壤。

（3）增加觀測點數量可以使參數更快地收斂于其真值，進而提高土壤剖面壓力水頭的預測精度。