王 正，孫兆軍，王 旭

(1.寧夏大學(xué)土木與水利工程學(xué)院，銀川 750021；2.寧夏大學(xué)新華學(xué)院，銀川 750021)

0 引 言

近年來，隨著人們對環(huán)境質(zhì)量重要性的認(rèn)識的提高，研究水分和可溶性鹽在非飽和土壤中的運(yùn)移機(jī)理已成為眾多領(lǐng)域的熱點(diǎn)[1]。而可溶性鹽則是通過水分在土壤中運(yùn)移來進(jìn)入地下水中的，因此，土壤中的水分運(yùn)動(dòng)對其有著至關(guān)重要的影響。為了更準(zhǔn)確地從定量的方面來揭示水分和鹽分的運(yùn)移規(guī)律，就需要獲取更為精確的輸入?yún)?shù)，而土壤水分特征曲線(Water retention characteristics)是其中最為重要的參數(shù)之一[2]。

對于刻畫土壤水分運(yùn)動(dòng)與含水率的土壤水分特征曲線問題，國內(nèi)外學(xué)者先后提出了許多數(shù)學(xué)模型和計(jì)算方法，其中吻合度最好的就是Van Genuchten模型[3](簡稱VG模型)，而VG方程是一個(gè)高度非線性的數(shù)學(xué)模型，含有眾多的未知參數(shù)，其表達(dá)式如下：

(1)

另外,求解方程(1)中的參數(shù)問題，本質(zhì)就是利用實(shí)驗(yàn)實(shí)測數(shù)據(jù)來擬合這些參數(shù)的問題。一般情況下，實(shí)測數(shù)據(jù)的樣本容量往往較小，故采用如下以回歸估計(jì)值與實(shí)際觀測值的偏差平方和最小為原則構(gòu)造魯棒性較好的優(yōu)化正則函數(shù)[4]：

(2)

式中：θi為實(shí)測土壤含水率；N為實(shí)驗(yàn)次數(shù)。

從本質(zhì)上講，式(2)隸屬于非線性問題的范疇，很多方法都可解決該問題。一般地，有直接法和間接法。通過實(shí)驗(yàn)數(shù)據(jù)擬定土壤水分特征曲線，也就是直接法，是一種最快捷的方式。而直接法中又有許多實(shí)驗(yàn)室方法和田間方法，例如測定土壤基水勢就有負(fù)壓力法、砂型漏斗法、壓力儀法等。雖然這些實(shí)驗(yàn)方法可以使用，但會(huì)大量消耗人力物力，實(shí)驗(yàn)時(shí)間漫長，并且不能保證測定范圍和實(shí)驗(yàn)精度[5]。間接法是指利用實(shí)驗(yàn)測得的有限的樣本，在已知模型上結(jié)合智能算法和原先假定的參數(shù)、源匯項(xiàng)和邊界條件等約束求得較優(yōu)參數(shù)。最近，將待求解問題與智能算法相結(jié)合的方式被廣泛應(yīng)用，而且計(jì)算準(zhǔn)確性較高。

1 參數(shù)反演算法現(xiàn)狀

隨著計(jì)算機(jī)和數(shù)學(xué)理論的不斷發(fā)展，對于參數(shù)優(yōu)化的問題，很多學(xué)者提出了不同的優(yōu)化方法，例如仿生算法、層次貝葉斯方法、集合卡爾曼濾波算法等。本部分將著重介紹21世紀(jì)以來用于土壤水分特征曲線參數(shù)反演的不同算法，并研究它們在實(shí)驗(yàn)精度和算法性能方面的異同。

1.1 仿生算法

1.1.1 螢火蟲算法

螢火蟲算法(Firefly Algorithm，F(xiàn)A)是根據(jù)螢火蟲獨(dú)特的生理構(gòu)造和生物習(xí)性總結(jié)出來的一種算法，其主要思想就是螢火蟲之間會(huì)利用發(fā)光來相互吸引，它們聚集的地方就應(yīng)是最優(yōu)解。但是由于原始的算法會(huì)在最優(yōu)值附近出現(xiàn)振蕩的可能，隨后付強(qiáng)[6]等人提出了一種改進(jìn)的螢火蟲算法。建立了一種約束判定條件：每兩個(gè)螢火蟲的距離lij小于某一特定值時(shí)，減小其步長λ。

設(shè)定λ的起始值為0.3，當(dāng)lij大于0.3時(shí)，可用下式計(jì)算其位置：

(3)

其中，符號意義同參考文獻(xiàn)[6]中表示。

其主要步驟如下：

(1)設(shè)定算法參數(shù)的初始值；

(2)隨機(jī)生成種群的初始位置；

(3)決定螢火蟲移動(dòng)方向；

(4)用式(3)修正螢火蟲位置，找尋最佳位置；

(5)在更改位置后，重新計(jì)算相對亮度；

(6)結(jié)果與判定條件進(jìn)行對比，決定是否繼續(xù)進(jìn)行算法；

(7)算法結(jié)束，輸出結(jié)果。

實(shí)驗(yàn)結(jié)果表明，用這種算法得到的計(jì)算結(jié)果與實(shí)測值有很好的擬合度，尤其是在較低的土壤水吸力條件下，誤差相當(dāng)小。但是螢火蟲算法的理論支撐還不夠，而且也沒有一個(gè)完備的數(shù)學(xué)模型，所以仍需要完善現(xiàn)有的理論，與其他算法糅合也能夠提高算法的性能。

1.1.2 粒子群算法

粒子群算法是一種基于進(jìn)化的計(jì)算方法(Particle Swarm Optimization，PSO)，該算法是建立在已發(fā)現(xiàn)的飛鳥集群活動(dòng)規(guī)律基礎(chǔ)之上，然后利用人工智能得到的簡化模型。該算法的基本思想是由個(gè)體與種群共享其獲得的“食物”信息，使整個(gè)群體都向求解空間中的“食物”所處位置進(jìn)行運(yùn)動(dòng)，從而產(chǎn)生從無序到有序的演化過程以獲得最優(yōu)解。

劉利斌[4]等、馬亮[7]完善了最初的算法，雖二人提出的算法結(jié)構(gòu)有些不同，但其大致步驟是相類似的，其主要步驟總結(jié)如下：

(1)參數(shù)初始化。計(jì)算各粒子對應(yīng)的目標(biāo)向量、確定種群規(guī)模、隨機(jī)產(chǎn)生θs、θr、α、m等參數(shù)的初始值。

(2)從初始值開始進(jìn)行迭代，得到位置更優(yōu)的初始值；

(3)計(jì)算粒子適應(yīng)度；

(4)持續(xù)更新適應(yīng)度和位置；

(5)對比適應(yīng)度值，得到全局最優(yōu)位置；

(6)比對結(jié)果和判定條件，決定是否繼續(xù)進(jìn)行。

計(jì)算結(jié)果表明，改進(jìn)后的粒子群算法有了更快的搜索速度，并且其精度也有很大的提高，因此，該算法具有較強(qiáng)的全局優(yōu)化能力和計(jì)算結(jié)果準(zhǔn)確的特點(diǎn)；另外，迭代次數(shù)少、魯棒性強(qiáng)、效率高也是該算法最主要的優(yōu)點(diǎn)。

1.2 層次貝葉斯方法

層次貝葉斯方法(Hierarchical Bayesian Method，HBM)于1996年由Berliner[8]提出，將復(fù)雜問題利用條件概率理論，逐級分解為不同層次的簡單后驗(yàn)概率求解問題[9]。而且各層次之間并不是孤立的，它們是由條件概率聯(lián)系起來的。層次貝葉斯方法在數(shù)據(jù)分析方面具有很大的潛力，現(xiàn)已經(jīng)成功應(yīng)用于其他領(lǐng)域[10-13]。

一般地，層次貝葉斯方法有數(shù)據(jù)S、過程Q和參數(shù)C等三層，數(shù)據(jù)S包括實(shí)地實(shí)測數(shù)據(jù)、反演同化數(shù)據(jù)和模型計(jì)算數(shù)據(jù)等；過程Q表示土壤濕度等參數(shù)的時(shí)空分布。其主要原理為：定義數(shù)據(jù)、過程和參數(shù)三層為隨機(jī)變量，并構(gòu)建條件概率模型，在已有數(shù)據(jù)條件下，將數(shù)據(jù)同化問題轉(zhuǎn)化為推理過程和參數(shù)的后驗(yàn)概率分布：

(4)

式中：p(×)為某事件的概率；p(a|b)表示b條件下a發(fā)生概率。

根據(jù)式(4)，層次貝葉斯方法數(shù)據(jù)模型為第一層，過程模型為第二層，參數(shù)模型為第三層，因此，定義了這三層模型之后，就能推導(dǎo)后驗(yàn)概率。另外，這三個(gè)模型還能夠進(jìn)行拆分，這樣做的目的是用其他的觀測數(shù)據(jù)或過程就能完全確定這些模型，這就實(shí)現(xiàn)了用簡單條件概率確定復(fù)雜的模型，層次貝葉斯方法也可以用于確定復(fù)雜模型的參數(shù)[14]。

(1)數(shù)據(jù)模型。數(shù)據(jù)模型表示X和θD在真實(shí)過程條件下與數(shù)據(jù)Y的概率關(guān)系，其模型為：

(5)

式中，數(shù)據(jù)Y=(Y1,Y2,…,Yn)表示不同數(shù)據(jù)源數(shù)據(jù)，X=(X1,X2,…,Xn)為每個(gè)數(shù)據(jù)源對應(yīng)尺度的過程，θD為參數(shù)，p(Yi|Xi,θD)表示為第i個(gè)數(shù)據(jù)源的數(shù)據(jù)模型。當(dāng)由不同分辨率的多源數(shù)據(jù)構(gòu)成輸入數(shù)據(jù)時(shí)，層次貝葉斯方法可分別建立相對應(yīng)的數(shù)據(jù)模型，并依賴于真實(shí)過程和參數(shù)，這就是說層次貝葉斯方法在解決多源數(shù)據(jù)問題方面具有一定優(yōu)勢。

(2)過程模型。層次貝葉斯方法中最為關(guān)鍵一步就是構(gòu)建過程模型，而過程模型需要對物理過程有大量的經(jīng)驗(yàn)信息和深刻的理解，最后建立過程X的時(shí)空分布模型p(X|θp)。通常來說，可將過程模型拆分為若干層次的子模型，其形式為：

(6)

式中：X=(X1,X2,…,Xn)為多尺度不同源數(shù)據(jù)所對應(yīng)的過程；θp為參數(shù)。

(3)參數(shù)模型。假定參數(shù)之間均相互獨(dú)立，則模型可定義為：

p(θD,θp)=p(θD)p(θp)

(7)

將式(5)～(7)代入式(4)中，數(shù)據(jù)同化要推理的概率為：

(8)

用層次貝葉斯方法來解決數(shù)據(jù)同化問題是一種非常有效的方法[15]，層次貝葉斯方法有以下幾個(gè)明顯區(qū)別于其他數(shù)據(jù)同化算法的特點(diǎn)：

(1)以貝葉斯理論為基礎(chǔ)，將估計(jì)狀態(tài)的后驗(yàn)概率作為目標(biāo)，同時(shí)完全考慮給定的觀測數(shù)據(jù)和模型數(shù)據(jù)；

(2)同化的數(shù)據(jù)中難免存在不確定性，而且對結(jié)果有重要影響，該方法能克服到這些缺點(diǎn)，也能脫離線性和高斯假設(shè)的束縛；

(3)分層考慮各個(gè)模型，并整合為一體，這樣處理的方式便于計(jì)算，而且在理論和實(shí)際應(yīng)用中更加合理方便。

1.3 集合卡爾曼濾波算法

Evensen[16]基于卡爾曼濾波提出了集合卡爾曼濾波(Ensemble Kalman Filter，En-KF)算法，用集合卡爾曼濾波研究土壤水分時(shí)，主要關(guān)注的是狀態(tài)變量的估計(jì)和預(yù)測[17-21]。

En-KF是基于蒙特卡洛的方法形成每一個(gè)元素的數(shù)據(jù)向前積分來預(yù)報(bào)誤差的統(tǒng)計(jì)量集合。一般用均值來近似最佳估計(jì)，用集合的取值范圍來表征參數(shù)估計(jì)的不確定性。

用EnKF來估計(jì)參數(shù)的主要步驟如下：

(1)形成初始變量和初始猜想集合。

(2)預(yù)報(bào)。集合中的每個(gè)元素向前積分，得到土壤持水率以矩陣形式存在的預(yù)報(bào)結(jié)果，并存儲在搜索空間中。

(3)產(chǎn)生觀測值的集合。

(4)分析。利用觀測算子分析其結(jié)果，將分析的結(jié)果存儲到搜索空間中。

(5)判定結(jié)果是否滿足要求。如果不滿足，返回第2步，否則停止并輸出。

Evensen[22,23]等給出更詳細(xì)的程序。

EnKF方法優(yōu)化的結(jié)果不受初始條件的影響，且具有很強(qiáng)的適應(yīng)性和更快的收斂性，還能夠同時(shí)考慮土壤水分函數(shù)模型以及數(shù)據(jù)觀測不確定性，大大降低了獲得偽優(yōu)參數(shù)值的風(fēng)險(xiǎn)。

2 存在的問題

進(jìn)入21世紀(jì)以來，我國在土壤水分滲透和水鹽運(yùn)移方面取得了許多成績，但是在國際科技不斷發(fā)展的大環(huán)境下，我國在土壤水分?jǐn)?shù)值模擬方面仍不能跟上世界的步伐。急需一大批科技工作者投入到這項(xiàng)事業(yè)當(dāng)中來，同時(shí)，也需要從更高的角度冷靜地思考仍未解決的問題。存在的主要問題歸納如下：

(1)嚴(yán)重依賴仿真技術(shù)，缺乏地質(zhì)考察經(jīng)驗(yàn)和理論。由于實(shí)驗(yàn)不僅可以發(fā)展理論還可以驗(yàn)證理論，因此，本文認(rèn)為實(shí)驗(yàn)是提升理論最有效的方式之一。建議國家有關(guān)部門要加大投入，要讓更多的人才參與到這項(xiàng)事業(yè)中的，同時(shí)讓他們戒除浮躁，安心工作。大量的經(jīng)驗(yàn)表明，只有通過大量的試驗(yàn)，獲取更多的數(shù)據(jù)樣本，才有利于發(fā)掘新的認(rèn)識。

(2)輕視具體地質(zhì)條件，隨意主觀地刪改數(shù)據(jù)，為了達(dá)到較好的擬合效果，隨意忽視一些現(xiàn)象或邊界條件。這種潮流，在現(xiàn)在的學(xué)術(shù)界尤為盛行。沒有深入考察當(dāng)?shù)氐牡刭|(zhì)情況下，就生搬硬套，必然只會(huì)產(chǎn)生片面的、帶主觀臆想的甚至是錯(cuò)誤的結(jié)論來。無數(shù)實(shí)踐證明，擬合效果不佳往往都是模型有缺失，而這又是不認(rèn)真考察地質(zhì)條件的結(jié)果。所以，結(jié)合當(dāng)?shù)卣鎸?shí)情況建立符合實(shí)際的模型，對客觀反映地區(qū)具體條件大有裨益。

(3)盲目追求參數(shù)擬合軟件?，F(xiàn)在有很多商用軟件都可以進(jìn)行參數(shù)反演的研究，但是不能過分依賴這些軟件，不是說有了軟件就萬無一失了，而是要在軟件的基礎(chǔ)上加入符合當(dāng)?shù)靥厣囊蛩?，只有這樣擬合出來的結(jié)果才能令人信服。

(4)模型雖多，但是沒有證明模型的可信度。早期一些學(xué)者提出來的模型和現(xiàn)如今的實(shí)際情況已相去甚遠(yuǎn)，原因就在于當(dāng)時(shí)的地質(zhì)條件已發(fā)生了變化，而且在主觀認(rèn)識上也有了不同，因此，以早期的模型來擬合現(xiàn)在的參數(shù)，恐怕也只能是刻舟求劍、朝花夕拾了。所以，要敢于挑戰(zhàn)權(quán)威，在已有的模型基礎(chǔ)上加以改進(jìn)優(yōu)化，這樣才是科學(xué)發(fā)展的規(guī)律。

3 發(fā)展方向與展望

近年來，已將眾多智能算法和技術(shù)應(yīng)用到土壤水分曲線參數(shù)反演的研究當(dāng)中來，但是隨著更多更好的、更新的方法的提出，在這方面的探索仍然任重道遠(yuǎn)。下面提出幾點(diǎn)發(fā)展的方向與展望：

(1)微波遙感技術(shù)可實(shí)時(shí)測量觀測，而且對地表環(huán)境的要求不高，能避免其他方式在土壤水分監(jiān)測應(yīng)用中的缺點(diǎn)和短板，這就為研究土壤水分提供了一個(gè)全新的方法和途徑。但是，由于目前的研究尚未成熟，所以有許多問題仍未完全解決，例如，一旦有植被附著在土壤的表面，微波會(huì)發(fā)生時(shí)空變異，會(huì)直接影響水分模型的參數(shù)；再如，雷達(dá)監(jiān)測雖然具有一定的穿透能力，但是也只是針對淺層土壤而言，對于深層土壤水分的探測仍然無能為力，無法探測到作物的根系，這樣會(huì)得到片面的結(jié)果。

(2)影響土壤含水率的因素很復(fù)雜，包括溫度、濕度、有機(jī)質(zhì)含量、礦物類型、滯后、陽離子交量和陽離子類型等，還需要進(jìn)一步考慮這些因素對水分曲線的影響機(jī)理。

(3)土壤結(jié)構(gòu)影響土壤儲水性能。土壤入滲能力、持水能力及土壤水分特征參數(shù)直接受土壤孔隙數(shù)量與質(zhì)量的影響，而團(tuán)聚體含量又影響土壤孔隙數(shù)量與質(zhì)量，所以要進(jìn)一步來研究土壤結(jié)構(gòu)性與土壤水分特征曲線的關(guān)系，也就要建立基于土壤傳遞函數(shù)的可以反映土壤入滲、蓄水能力的計(jì)算模型，為后續(xù)的研究提供簡便的方法。

4 結(jié) 語

預(yù)測土壤水分特征曲線的間接方法在今后的研究中有很廣泛的前景。而且，尤其在較大范圍的模型研究中，間接方法的優(yōu)勢就突顯出來。就目前來說，土壤水分特征曲線參數(shù)反演主要依靠的是智能算法，例如仿生算法、數(shù)據(jù)同化算法等，這些算法在不同程度上都具有很好的反演精度，但是在模型的使用和算法的結(jié)構(gòu)上仍然需要進(jìn)一步提高。另外，其他新技術(shù)的摸索和使用，例如遙感技術(shù)，在很大程度上也能提供更廣闊的平臺。

[1] 李 峰，繳錫云，李盼盼，等.田間土壤水分特征曲線參數(shù)反演[J].河海大學(xué)學(xué)報(bào)(自然科學(xué)版)，2009，(4)：373-376.

[2] 劉建立，徐紹輝，劉 慧.估計(jì)土壤水分特征曲線的間接方法研究進(jìn)展[J].水利學(xué)報(bào)，2004，(2)：68-73.

[3] Van Genuchten M.Th. A closed form equation for predicting the hydraulic conductivity of unsaturated soils [J].Soil Sci.Soc.Am.J.1980，44(5)：892-898.

[4] 劉利斌，歐陽艾嘉，樂光學(xué)，等.基于混合粒子群的土壤水分特征曲線參數(shù)優(yōu)化[J].計(jì)算機(jī)工程與應(yīng)用，2011,47(35)：218-221.

[5] 唐亞莉，董文明，趙晶晶，等.優(yōu)化算法確定土壤水分特征曲線的分析[J].新疆大學(xué)學(xué)報(bào)(自然科學(xué)版)，2006，23(2)：240-243.

[6] 付 強(qiáng)，蔣睿奇，王子龍，等.基于改進(jìn)螢火蟲算法的土壤水分特征曲線參數(shù)優(yōu)化[J].農(nóng)業(yè)工程學(xué)報(bào)，2015,31(11)：117-122.

[7] 馬 亮.基于改進(jìn)粒子群的土壤水分特征曲線參數(shù)優(yōu)化研究[J].節(jié)水灌溉，2014，(11):17-20.

[8] Berliner L M，Hanson K，Silver R. Hierarchical bayesian time series models[M]. Maximum Entropy and Bayesian Method.Dordrecht：Kluwer Academic Public，1996:15-23.

[9] Cressie N，Wikle C K. Statistics for Spatio-temporal Data[M]. New Jersey：Wiley，John & Sons，2010.

[10] Gelman A，Carlin J B，Stern H S，et al.Bayesian Data Analysis[M].New York：Chapman & Hall/CRC，2003.

[11] Gelfand A E, Diggle P, Guttorp P, et al. Handbook of Spatial Statistics[M].New York：Chapman and Hall /CRC，2010.

[12] Gelfand A E, Sahu S K. Combining monitoring data and computer model output in assessing environmental exposure[M].The Oxford Handbook of Applied Bayesian Analysis. Oxford：Oxford University Press，2009.

[13] Sahu S K，Yip S，Holland D M. Improved space-time forecasting of next day ozone concentrations in the eastern US[M]. Atmospheric Environment，2009，43( 3)：494-501.

[14] Liu Y，Gupa H V. Uncertainty in hydrologic modeling toward an integrated data assimilation framework[J].Water Resources Research，2007，43：10.1029 /2006WR005756.

[15] Wikle C K，Berliner L M.A Bayesian tutorial for data assimilation[J]. Physica D，2006，230( 1 /2) ： 1-16.

[16] Evensen G. Sequential data assimilation with a nonlinear quasi-geostrophic model using monte carlo methods to forecast error statistics [J]. Journal of Geophyscial Research，1994，99：10 143-10 162.

[17] Reichle R，Walker J，Koster R，et al.Extended versus ensemble Kalman filtering for land data assimilation [J]. Journal of Hydrometeorology，2002，(3)：728-740.

[18] Crow W T. Correcting land surface model predictions for the impact of temporally sparse rainfall rate measurements using an ensemble Kalman filter and surface brightness temperature observation [J].Journal of Hydrometeorology, 2003，(4)： 960-973.

[19] Zhang S W，Li H R，Zhang W D，et al. Estimating the soil moisture profile by assimilating near-surface observations with the Ensemble Kalman Filter (EnKF) [J]. Advances in Atmospheric Sciences，2006，22(6)： 936-945.

[20] 黃春林，李 新.土壤水分同化系統(tǒng)的敏感性試驗(yàn)研究[J].水科學(xué)進(jìn)展， 2006，17(4)： 457-465.

[21] 茍浩鋒，劉彥華，張述文，等.評估集合卡曼濾波反演土壤濕度廓線的性能[J].地球科學(xué)進(jìn)展，2010，25(4)：400-407.

[22] Evensen G.The ensemble Kalman filter：Theoretical formulation and practical implementation [J].Ocean Dynamic，2003，53：343-367.

[23] Evensen G. Sampling strategies and square root analysis schemes for the EnKF[J].Ocean Dynamic，2004，54：539-560.