基于集合卡爾曼濾波法的二維土壤水流狀態(tài)變量和參數(shù)聯(lián)合估計(jì)

劉 琨，黃冠華

（1.中國農(nóng)業(yè)大學(xué) 水利與土木工程學(xué)院，北京 100083；2.中國-以色列國際農(nóng)業(yè)研究培訓(xùn)中心，北京 100083）

1 研究背景

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

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

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

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

2 理論與方法

2.1 數(shù)據(jù)同化系統(tǒng)數(shù)據(jù)同化系統(tǒng)由模型算子、觀測算子和同化算法組成。首先模型算子根據(jù)上一時(shí)刻的狀態(tài)向量向前預(yù)測得到狀態(tài)向量的預(yù)測值；觀測算子將狀態(tài)向量轉(zhuǎn)化為觀測值所對(duì)應(yīng)的模型預(yù)測值；最后同化系統(tǒng)利用觀測信息，對(duì)預(yù)測值進(jìn)行同化運(yùn)算，得到狀態(tài)向量的分析值。分析值作為下一時(shí)刻的模型算子輸入變量，重復(fù)以上過程直至模型模擬結(jié)束。

2.1.1 模型算子 本研究中的模型算子為二維土壤水流運(yùn)動(dòng)模型，即二維Richards方程:

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

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

當(dāng)有觀測數(shù)據(jù)時(shí)，利用觀測數(shù)據(jù)計(jì)算得到當(dāng)前時(shí)刻狀態(tài)向量的分析值:

式中:H為觀測算子；Pef為預(yù)測值觀測誤差協(xié)方差；Re為觀測誤差協(xié)方差；di，t+1為t+1時(shí)刻各集合成員的觀測向量。詳細(xì)的EnKF計(jì)算流程請(qǐng)參見文獻(xiàn)［3，16］。

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

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

表1 土壤類型及其土壤水力特性參數(shù)

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

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

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

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

圖1 不同觀測點(diǎn)位置布置方案（實(shí)心點(diǎn)為觀測點(diǎn)，單位:cm）

3 結(jié)果分析與討論

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

參數(shù)KS標(biāo)準(zhǔn)差隨著同化次數(shù)增加逐漸減小（見圖2）。垂向布置方案的KS標(biāo)準(zhǔn)差變化趨勢類似，經(jīng)過3次同化標(biāo)準(zhǔn)差值已經(jīng)很??；水平向布置條件下，不同方案的標(biāo)準(zhǔn)差變化差異較大，隨著觀測點(diǎn)埋深的增加，參數(shù)估計(jì)效果變差，標(biāo)準(zhǔn)差降低的速度變小。這是由于當(dāng)入滲水分沒有運(yùn)動(dòng)到觀測點(diǎn)埋深時(shí)，觀測點(diǎn)的壓力水頭不變，同化系統(tǒng)尚未利用觀測信息更新狀態(tài)變量和參數(shù)，因此參數(shù)估計(jì)結(jié)果較差。圖3為不同觀測點(diǎn)布置方案進(jìn)氣值參數(shù)的演化，進(jìn)氣值參數(shù)α的演化與KS的演化趨勢一致，KS同化效果較好的方案，對(duì)α的估計(jì)效果也較好，這是由于同化計(jì)算中參數(shù)KS和α是同時(shí)更新的。與參數(shù)KS不同，8種方案參數(shù)α的標(biāo)準(zhǔn)差均下降較快，這是由于參數(shù)α的初始方差設(shè)置較小。

圖2 不同觀測點(diǎn)布置方案飽和導(dǎo)水率演化

圖3 不同觀測點(diǎn)布置方案進(jìn)氣值參數(shù)演化

表2 不同觀測點(diǎn)位置KS估計(jì)MRE

表3 不同觀測點(diǎn)位置α估計(jì)MRE

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

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

3.2 觀測點(diǎn)數(shù)量影響在觀測點(diǎn)位置影響結(jié)果的基礎(chǔ)上選擇每種土壤最優(yōu)的三種布置方案組合得到新的布置方案，觀測點(diǎn)數(shù)量分別為4、8和12，觀測點(diǎn)布置見圖5。

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

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

圖5 不同觀測點(diǎn)數(shù)量布置方案（實(shí)心點(diǎn)為觀測點(diǎn) 單位:cm）

圖6 不同觀測點(diǎn)數(shù)量條件下飽和導(dǎo)水率和進(jìn)氣值參數(shù)演化

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

與傳統(tǒng)的參數(shù)估計(jì)方法相比，EnKF方法在處理不確定性方面有著明顯的優(yōu)勢。例如在Levenberg-Marquardt算法中，觀測數(shù)據(jù)被認(rèn)定為是衡量預(yù)測精度的唯一標(biāo)準(zhǔn)，模型的不確定性完全歸因于參數(shù)，并沒有考慮觀測數(shù)據(jù)誤差。而在模型構(gòu)建過程中對(duì)物理過程的認(rèn)識(shí)不足和簡化以及觀測手段的局限，模型結(jié)構(gòu)和觀測數(shù)據(jù)都具有不確定性，因此僅將不確定歸因于參數(shù)是不恰當(dāng)?shù)?。EnKF方法綜合觀測數(shù)據(jù)和模型預(yù)測的不確定性，可以得到參數(shù)的優(yōu)化值。本研究中我們已通過數(shù)值實(shí)驗(yàn)從理論上證明了狀態(tài)變量與參數(shù)聯(lián)合估計(jì)效果，在未來的研究中還需結(jié)合田間實(shí)驗(yàn)進(jìn)一步檢驗(yàn)該方法的精度。

4 結(jié)論

本文采用EnKF方法開展二維土壤水流運(yùn)動(dòng)模型狀態(tài)變量和參數(shù)聯(lián)合估計(jì)研究，采用數(shù)值實(shí)驗(yàn)方法研究了粉壤土、壤土和砂壤土在線源入滲條件下EnKF方法對(duì)飽和導(dǎo)水率和進(jìn)氣值參數(shù)估計(jì)以及壓力水頭的同化效果，分析了觀測點(diǎn)布置方式和觀測點(diǎn)數(shù)量對(duì)同化效果的影響。主要結(jié)論如下:

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

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

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