基于鳥群算法的土壤水分特征曲線模型參數(shù)估算*

(山西省水利水電科學(xué)研究院，山西 太原 030002)

隨著科學(xué)技術(shù)不斷發(fā)展，計(jì)算機(jī)數(shù)值計(jì)算方法已經(jīng)成為解決土壤中高維非線性問(wèn)題的重要工具。水力參數(shù)差異是導(dǎo)致土壤水分時(shí)空變異的重要因素，而土壤水分特征曲線是關(guān)鍵的土壤水力參數(shù)之一。土壤水分特征曲線反映了土壤水分能量和數(shù)量之間的關(guān)系，是進(jìn)行土壤水分運(yùn)動(dòng)模擬的基礎(chǔ)條件。因此，如何精準(zhǔn)獲取土壤水分特征曲線模型參數(shù)對(duì)于準(zhǔn)確實(shí)現(xiàn)水分時(shí)空動(dòng)態(tài)預(yù)測(cè)具有十分重要的意義。

Van Genuchten模型具有較高的預(yù)測(cè)精度，并且模型參數(shù)的物理意義也最為明確，因此，該模型被廣泛用于土壤水分特征曲線的定量描述。目前，該模型參數(shù)的確定主要包括入滲試驗(yàn)法和數(shù)學(xué)計(jì)算法[1-2]。相對(duì)于入滲試驗(yàn)法，數(shù)學(xué)計(jì)算法更加簡(jiǎn)單方便。但采用前人的數(shù)學(xué)計(jì)算方法進(jìn)行參數(shù)估算時(shí)，也存在一定不足。李春友[3]采用單純形調(diào)優(yōu)法進(jìn)行了水力參數(shù)估算，結(jié)果表明該方法會(huì)陷入局部最優(yōu)，因此，它不是一種合適的算法。馬英杰[4]采用非線性阻尼最小二乘法進(jìn)行了水力參數(shù)的估算, 但模型預(yù)測(cè)結(jié)果受初始值影響較大，模型精度并不穩(wěn)定。因此，尋找一種更加合適的算法對(duì)于提高Van Genuchten模型預(yù)測(cè)精度是非常必要的。鳥群算法是一種基于模仿鳥類種群活動(dòng)過(guò)程(覓食、警戒和飛行)而提出的新型全局智能優(yōu)化算法。該算法具有調(diào)節(jié)參數(shù)少、收斂精度高和魯棒性能好等優(yōu)點(diǎn)，在電網(wǎng)優(yōu)化調(diào)度、水庫(kù)優(yōu)化調(diào)度和年徑流量預(yù)測(cè)等方面取得了較好的應(yīng)用效果[4-6]，但在農(nóng)業(yè)土壤環(huán)境方面未見(jiàn)相關(guān)報(bào)道。該算法對(duì)Van Genuchten模型的參數(shù)優(yōu)化能力及效果也尚不清楚，有待進(jìn)一步研究。為此，本文將采用該算法對(duì)Van Genuchten模型參數(shù)進(jìn)行求解計(jì)算和精度評(píng)價(jià)，以期為Van Genuchten模型參數(shù)的估算提供一種新途徑。

1 Van Genuchten模型參數(shù)求解方法構(gòu)建

1.1 Van Genuchten模型簡(jiǎn)介

本文采用Van Genuchten模型對(duì)土壤水勢(shì)與含水率的關(guān)系進(jìn)行描述，具體數(shù)學(xué)表達(dá)式如下[7]：

(1)

式中θ——土壤含水率，cm3/cm3；

h——土壤水吸力，cm；

θs——土壤飽和含水率，cm3/cm3；

θr——土壤殘余含水率，cm3/cm3；

α、m、n——土壤水分特征曲線形狀參數(shù)，m=1-1/n，n>1。

1.2 構(gòu)造目標(biāo)函數(shù)

由式(1)可知，θs，θr，α和n為4個(gè)獨(dú)立參數(shù)。為此，本文構(gòu)造的目標(biāo)函數(shù)見(jiàn)式(2):

(2)

式中θs——實(shí)測(cè)含水率，cm3/cm3；

hi——實(shí)測(cè)含水率θi對(duì)應(yīng)的土壤水吸力，cm；

θ(hi,X)——含水率計(jì)算值，cm3/cm3；

X——待求參數(shù)向量(θs，θr，α,m)；

N——樣本數(shù)。

1.3 鳥群(GSA)算法基本思想

a.覓食行為。通常，鳥類在覓食過(guò)程中都需要隨時(shí)保持警惕性，即覓食行為和警戒行為需要進(jìn)行切換，并且該過(guò)程為一個(gè)隨機(jī)發(fā)生過(guò)程。這可以用數(shù)學(xué)語(yǔ)言描述為：在(0,1)范圍內(nèi)，當(dāng)均勻分布的隨機(jī)小數(shù)小于P[P∈(0,1)]時(shí)，鳥群進(jìn)行覓食，否則，保持警惕。在覓食過(guò)程中，每只鳥都能對(duì)食物及群體以往最好經(jīng)驗(yàn)進(jìn)行及時(shí)記錄和更新，并且這種群體信息是共享的。這些可以采用數(shù)學(xué)模型進(jìn)行描述，見(jiàn)式(3):

(3)

式中j——1,2,…,D;

C——非負(fù)常數(shù)，認(rèn)知行為的加速因子；

S——非負(fù)常數(shù)，是群體行為的加速因子；

Pi，j——第i只鳥最佳位置；

gj——群體當(dāng)前最佳位置。

b.警戒行為。在警戒行為發(fā)生時(shí)，個(gè)體均會(huì)向群體靠攏。由于個(gè)體存在警覺(jué)性差異，它們向群體中心移動(dòng)的難易程度也會(huì)各不相同，并且在這種移動(dòng)過(guò)程中還會(huì)伴隨著相互競(jìng)爭(zhēng)。該警戒行為可采用下式進(jìn)行描述：

(4)

(5)

(6)

式中k(k≠i)——(1,N)之間的正整數(shù)；

a1和a2——(0,2)之間的實(shí)數(shù)；

PFiti——第i只鳥最佳適應(yīng)度；

sumFit——鳥群適應(yīng)度之和；

ε——常數(shù)；

meanj——群體第j維平均位置。

c.飛行行為。在鳥群尋找食物過(guò)程中必然會(huì)發(fā)生飛行行為。由于個(gè)體差異，個(gè)體間的生存能力會(huì)存在差異。部分群體能夠及時(shí)準(zhǔn)確發(fā)現(xiàn)食物，成為生產(chǎn)者；而另外群體只能跟隨生產(chǎn)者獲取食物，成為食物索取者。生產(chǎn)者和隨機(jī)者的出現(xiàn)為隨機(jī)過(guò)程。飛行行為可采用下式進(jìn)行描述：

(7)

(8)

式中randn(0,1)——服從方差為0，均值為1的Gaussian分布的隨機(jī)數(shù)；

k——(1,N)之間的正整數(shù)；

FL[FL∈(0,2)]——索取者會(huì)跟隨生產(chǎn)者去尋找食物。

另外假設(shè)每只鳥會(huì)在FQ時(shí)間間隔內(nèi)飛往另一個(gè)地方，F(xiàn)Q為正整數(shù)。

1.4 Van Genuchten模型參數(shù)尋優(yōu)流程

Van Genuchten模型參數(shù)尋優(yōu)主要包括如下步驟：

步驟1：初始化鳥群規(guī)模N、迭代次數(shù)M、飛行頻率FQ、覓食概率P和常量參數(shù)。

步驟2：根據(jù)適應(yīng)度函數(shù)計(jì)算個(gè)體適應(yīng)度值，并由此確定個(gè)體和群體的最佳位置。

步驟3：迭代開(kāi)始，判斷t%FQ是否有余數(shù)。

a.有余數(shù)，判斷鳥群發(fā)生覓食行為還是警戒行為。若為覓食行為，采用式(3)更新位置；若為警戒行為，采用式(4)更新位置。

b.若沒(méi)有余數(shù)，判斷個(gè)體屬于生產(chǎn)者還是索取者。若為生產(chǎn)者，采用式(7)更新位置；若為索取者，采用式(8)更新位置。

步驟4：若第i只鳥當(dāng)前位置優(yōu)于先前自身最佳位置，則將當(dāng)前位置更新為自身最佳位置，同時(shí)更新鳥群當(dāng)前最佳位置。

步驟5: 判斷是否達(dá)到最大迭代次數(shù)，若迭代結(jié)束轉(zhuǎn)至步驟6; 否則迭代數(shù)加1，轉(zhuǎn)回步驟3。

步驟6: 輸出鳥群最佳適應(yīng)度值及鳥群最佳位置所包含的θs，θr，α和n為4個(gè)獨(dú)立參數(shù)。

1.5 模型評(píng)價(jià)

本文主要從模型預(yù)測(cè)值與實(shí)測(cè)值的線性一致性、平均相對(duì)誤差(MAPE)、均方根誤差(RMSE)三個(gè)方面對(duì)模型預(yù)測(cè)效果進(jìn)行評(píng)價(jià)。MAPE和RMSE的計(jì)算式分別為：

(9)

(10)

式中WR——試驗(yàn)觀測(cè)值,mg/kg;

WL——模型預(yù)測(cè)值,mg/kg;

N——樣本數(shù)。

2 數(shù)據(jù)來(lái)源

本文采用室內(nèi)土壤水分特征曲線試驗(yàn)資料對(duì)Van Genuchten模型參數(shù)進(jìn)行估算。土壤取樣點(diǎn)位于山西省水利水電科學(xué)研究院節(jié)水高效示范基地。該基地位于太原市小店區(qū)，屬暖溫帶大陸性氣候區(qū)，年平均降水量為495.5mm，年平均氣溫為9.5℃。土壤經(jīng)風(fēng)干、碾壓、混合、過(guò)篩(篩孔徑為2mm)后備用。土壤質(zhì)地為黏壤土。首先通過(guò)壓力膜儀法獲得土壤水分特征曲線數(shù)據(jù)資料，然后采用單純形調(diào)優(yōu)法、非線性最小二乘法和鳥群算法分別對(duì)Van Genuchten模型參數(shù)進(jìn)行優(yōu)化求解。

3 參數(shù)估算分析

3.1 Van Genuchten模型參數(shù)估算結(jié)果分析

本文采用鳥群算法對(duì)Van Genuchten模型參數(shù)進(jìn)行求解，獲得的最優(yōu)模型參數(shù)見(jiàn)下表。為了進(jìn)一步評(píng)價(jià)該方法的可靠性和優(yōu)越性，采用傳統(tǒng)的單純形調(diào)優(yōu)法和非線性最小二乘法獲得的估算結(jié)果作為對(duì)比，模型參數(shù)見(jiàn)下表。由下表可以看出，經(jīng)單純形調(diào)優(yōu)法和非線性最小二乘法計(jì)算后的土壤飽和含水率(θs)數(shù)值相同，獲得的土壤殘余含水率θr數(shù)值也相同。而經(jīng)過(guò)鳥群算法處理后的θr要比前兩者結(jié)果大11.99%，θs要比前兩者結(jié)果小0.73%。由此說(shuō)明，鳥群算法對(duì)θr的估算結(jié)果影響較大。從形狀參數(shù)來(lái)看，經(jīng)單純形調(diào)優(yōu)法和鳥群算法估算后的形狀參數(shù)α數(shù)值相同，而非線性最小二乘法估算結(jié)果要比其余2種方法的估算結(jié)果高3.08%，由此說(shuō)明3種方法處理后的α數(shù)值差異不明顯。由下表還可以看出，經(jīng)3種方法處理后的形狀參數(shù)n大小表現(xiàn)為：?jiǎn)渭冃握{(diào)優(yōu)法>非線性最小二乘法>鳥群算法。相對(duì)單純形法而言，經(jīng)非線性最小二乘法和鳥群算法估算后的參數(shù)n要分別比前者小0.98%和2.00%，說(shuō)明3種算法對(duì)形狀參數(shù)n的計(jì)算結(jié)果無(wú)顯著影響。由下表還可以看出，3種方法估算方法的均方差大小表現(xiàn)為：非線性最小二乘法>單純形調(diào)優(yōu)法>鳥群算法。綜上所述，與傳統(tǒng)的單純形調(diào)優(yōu)法和非線性最小二乘法相比較，基于鳥群算法的參數(shù)估算方法主要對(duì)Van Genuchten模型中參數(shù)θr的影響較大。基于鳥群算法的Van Genuchten模型參數(shù)估算方法的結(jié)果更加準(zhǔn)確。

Van Genuchten模型參數(shù)估算結(jié)果表

3.2 Van Genuchten模型參數(shù)估算精度分析

圖1為不同估算方法下土壤含水率預(yù)測(cè)值與實(shí)測(cè)值間的線性關(guān)系。由圖1可以看出，經(jīng)3種算法處理后的估算值與實(shí)測(cè)值均呈現(xiàn)較好的線性關(guān)系。經(jīng)單純形調(diào)優(yōu)法、非線性最小二乘法和鳥群算法處理后的線性方程斜率分別為0.904、0.968和0.974，決定系數(shù)分別為0.987、0.992和0.996。由此說(shuō)明，經(jīng)3種方法處理后，實(shí)測(cè)值與預(yù)測(cè)值間線性一致性好壞表現(xiàn)為：鳥群算法>非線性最小二乘法>單純形調(diào)優(yōu)法。

圖1 不同估算方法下土壤含水率預(yù)測(cè)值與實(shí)測(cè)值間的線性關(guān)系

圖2為不同估算方法下土壤含水率預(yù)測(cè)值與實(shí)測(cè)值對(duì)比結(jié)果。由圖2可以看出，經(jīng)3種方法估算后的土壤含水率均隨土壤水吸力呈指數(shù)型衰減趨勢(shì)。由圖2可以看出，在0～2000cm和9000～12000cm范圍內(nèi)，3種算法處理后的結(jié)果較為接近，但在2000～9000cm范圍內(nèi)，3種算法處理后的結(jié)果差異較為明顯。由此說(shuō)明，不同算法主要對(duì)中間土壤水吸力段(2000～9000cm)的土壤含水率預(yù)測(cè)結(jié)果影響較大。從整體上來(lái)看，經(jīng)鳥群算法處理后的含水率預(yù)測(cè)值與實(shí)測(cè)值最為接近，單純形調(diào)優(yōu)法的處理結(jié)果次之，非線性最小二乘法的含水率預(yù)測(cè)值與實(shí)測(cè)值差異最明顯。

圖2 不同估算方法下土壤含水率預(yù)測(cè)值與實(shí)測(cè)值對(duì)比結(jié)果

為了進(jìn)一步準(zhǔn)確比較各種方法的估算效果，有必要對(duì)其進(jìn)行量化分析。圖3為不同算法估算值與實(shí)測(cè)值間的絕對(duì)百分誤差。由圖3可以看出，不同算法處理后的樣本誤差較為明顯。經(jīng)鳥群算法、非線性最小二乘法和單純形調(diào)優(yōu)法處理后的絕對(duì)百分誤差分別為1.00%～4.84%、0.74%～7.78%和0.57%～9.69%，平均相對(duì)誤差分別為3.50%、4.89%和4.31%。由此說(shuō)明，不同算法的估算精度好壞表現(xiàn)為：鳥群算法>單純形調(diào)優(yōu)法>非線性最小二乘法。鳥群算法的估算效果最為理想，單純形調(diào)優(yōu)法估算效果次之，非線性最小二乘法估算效果最差。

圖3 不同算法估算值與實(shí)測(cè)值間的絕對(duì)百分誤差

4 結(jié) 論

本文通過(guò)鳥群算法對(duì)Van Genuchten模型中的參數(shù)進(jìn)行了求解，并對(duì)參數(shù)估算精度進(jìn)行了評(píng)價(jià)。與單純形調(diào)優(yōu)法(SEM)和非線性最小二乘法(NLS)的估算結(jié)果相比較，基于鳥群算法(BSA)的參數(shù)估算方法獲得的土壤殘余含水率θr均比前兩者要大11.99%?；邙B群算法的Van Genuchten模型參數(shù)估算方法的平均相對(duì)誤差僅為3.50%，優(yōu)于傳統(tǒng)的單純形調(diào)優(yōu)法和非線性最小二乘法估算效果。