黃河下游沖積平原典型區(qū)包氣帶土壤水分特征曲線的模型優(yōu)選

湛 江，李志萍，2，于小朋

（1.華北水利水電大學(xué) 地球科學(xué)與工程學(xué)院，河南 鄭州 450045；2.水資源高效利用與保障工程河南省協(xié)同創(chuàng)新中心，河南 鄭州 450003）

土壤水分特征曲線是土壤水分和溶質(zhì)運(yùn)移必不可少的參數(shù)，對于研究水文循環(huán)、指導(dǎo)農(nóng)業(yè)生產(chǎn)具有重要意義。 該曲線描述了基質(zhì)吸力和含水率間的關(guān)系，其測定通常費(fèi)時費(fèi)力且成本較高，因此這種關(guān)系需要一種適用性廣泛的解析式來表達(dá)。 由于基質(zhì)吸力與含水率之間的高度非線性，因此建立一個適用性廣泛的表達(dá)式并不容易［1］。

從20 世紀(jì)60年代至今，不斷有水分特征曲線的經(jīng)驗(yàn)?zāi)Ｐ蛦柺溃辽儆?0 種經(jīng)驗(yàn)?zāi)Ｐ偷玫搅顺浞謶?yīng)用和發(fā)展［2］，應(yīng)用較為廣泛的有 Brooks Corey（BC）模型［3］、Gardner 模型［4］、Campbell 模型（C）［5］、van Genuchten（ VG） 模 型［6］、 Fredlund ＆ Xing（ FX） 模型［7］和 Kosugi 模型［8］等，其中以 VG 模型最為著名。近年來，許多研究表明VG 模型的表現(xiàn)優(yōu)于其他模型［9-12］。 但也有一些研究指出，由于土壤質(zhì)地和土壤物理化學(xué)性質(zhì)的不同，因此VG 模型并非最佳選擇［13］。 比如，當(dāng)土壤中砂粒含量較高（超過 50%）時，VG 模型的表現(xiàn)不佳［14-15］。 劉暢等［16］通過對比研究認(rèn)為斥水砂土的水分特征曲線更適合Kosugi 模型而并非 VG 模型和 BC 模型。 王志超等［17］構(gòu)建了殘膜土的水分特征曲線模型，以RPF-SWCC 模型（改進(jìn)的Fredlund 模型）為最優(yōu)模型，該模型優(yōu)于VG模型、BC 模型、Kosugi 模型。 許多模型的擬合效果都很好，但考慮到方程形式的復(fù)雜性和推廣性，研究者會采用方程形式簡單的模型以便于應(yīng)用。 董義陽等［18］研究表明古爾班通古特沙漠丘間地土壤水分特征曲線的VG 模型的擬合效果略優(yōu)于Gardner 模型，卻因Gardner 模型形式簡單、參數(shù)少而認(rèn)為Gardner模型更為適用。

以上研究多集中于擾動土，對于原狀土的相關(guān)研究較少，僅有研究表明FX 模型是原狀黃土土壤水分特征曲線的最優(yōu)模型［19］，但對于黃河下游地區(qū)的包氣帶深層土壤的相關(guān)研究鮮見報(bào)道。 本研究基于粒子群優(yōu)化算法對6 種經(jīng)驗(yàn)?zāi)Ｐ瓦M(jìn)行擬合精度比較，優(yōu)選出適用于研究區(qū)的最佳土壤水分特征曲線的數(shù)學(xué)模型。

1 研究區(qū)概況與研究方法

1.1 研究區(qū)概況

自古以來，黃河下游水患頻繁，沿岸生態(tài)環(huán)境屢遭破壞。 歷史上黃河的頻繁決口改道，伴隨著洪水沖刷地表和泥沙快速沉積，使得黃河下游地區(qū)的包氣帶土壤不斷得以更新，垂向分層明顯。

河南省蘭考縣位于黃河下游地區(qū)，北部緊鄰黃河。近代以來的“花園口決堤”事件，使得蘭考縣的地貌發(fā)生巨大變化，成為歷史上著名的黃泛區(qū)。 蘭考縣地層沉積具有代表性，能夠全面反映黃河下游泥沙的沉積過程；而且蘭考縣曾飽受風(fēng)沙災(zāi)害，土地貧瘠，是著名的鹽漬區(qū)。 但經(jīng)過近30 a 的土地改良等，該地區(qū)的農(nóng)業(yè)生產(chǎn)條件得到改善，因此蘭考縣地表土壤具有一定特殊性，能夠反映人類活動對土地的影響。

本研究選取集黃泛區(qū)、鹽漬區(qū)和農(nóng)業(yè)區(qū)于一身的蘭考 縣 閆 樓 鄉(xiāng) （ 114° 57′ 23″ E—115° 00′ 02″ E，34°54′27″N—34°52′31″N）作為黃河下游沖積平原典型區(qū)。 研究區(qū)位于蘭考縣中部，原為黃河故道［20］；南北長4 km、東西長4 km，總面積約16 km2；土地利用類型為農(nóng)田和村莊。

1.2 采樣與試驗(yàn)

1.2.1 采樣方案

試驗(yàn)采用人工掘進(jìn)法，依每個鉆孔按不同巖性分層取樣。 考慮地表土壤受土地利用（農(nóng)業(yè)、村莊和工程建設(shè)等）影響較大，本研究按4 層采集樣品，從地表向下依次為土層、粉土層、粉質(zhì)黏土層、粉砂層。 由于研究區(qū)北部農(nóng)田實(shí)為洼地，致使粉質(zhì)黏土層出露地表，因此第2 層粉土層樣品不足64 個，其余巖性均取64個樣品，共計(jì)取得233 個樣品，每個樣品取若干份原狀環(huán)刀樣品以供試驗(yàn)。

本研究采用網(wǎng)格法確定采樣點(diǎn)鉆孔，網(wǎng)格單元大小為500 m×500 m。 為避開建筑物和街道，采樣點(diǎn)位置稍偏離網(wǎng)格中心（見圖1）。 圖1 中顯示了4 個鉆孔的剖面柱狀圖，表層與第2 層均以地表以下30 cm 為界。 可以看出，東北部和南部高地的包氣帶土層共4層；位于研究區(qū)西北部洼地的LK0104 號鉆孔，由于粉土層的風(fēng)蝕剝離，因此包氣帶土層只有3 層。

圖1 研究區(qū)采樣點(diǎn)分布和代表鉆孔剖面柱狀圖

1.2.2 試驗(yàn)方法

采用QT-2012 型激光粒度儀測得土壤顆粒組成，以確定土壤質(zhì)地，質(zhì)地分類標(biāo)準(zhǔn)參考USDA（美國農(nóng)業(yè)部）土壤分類表［21］。

土壤水分特征曲線的測定方法很多，但考慮壓力膜儀吸力范圍廣，對于黏土等細(xì)質(zhì)地土壤可以測得完整曲線，因此本研究選擇壓力膜儀法測定土壤水分特征曲線。

1.3 研究方法

1.3.1 經(jīng)驗(yàn)?zāi)Ｐ瓦x擇

按照模型的普適性以及發(fā)展歷程，選取6 種代表性經(jīng)驗(yàn)?zāi)Ｐ瓦M(jìn)行比較和優(yōu)選，主要考慮這些模型的適用性、泛化性、方程類型和形式（參數(shù)數(shù)量），其表達(dá)式見表1。

表1 水分特征曲線經(jīng)驗(yàn)?zāi)Ｐ?/p>

BC 模型和C 模型提出的年代較早，其方程形式簡單，便于應(yīng)用，但對于非均質(zhì)土壤預(yù)測精度不高。VG 模型的優(yōu)點(diǎn)是適用于絕大部分土壤，其預(yù)測精度高，并且在全吸力段光滑連續(xù)，但是該方程缺點(diǎn)為：首先，與所有經(jīng)驗(yàn)?zāi)Ｐ鸵粯樱锢硪饬x不夠明確［22］；其次，對于非均質(zhì)明顯的擾動土壤（如雙峰、多峰結(jié)構(gòu)），其擬合效果一般。 Fredlund ＆ Xing 基于土壤孔隙分布模型建立了水分特征曲線經(jīng)驗(yàn)公式，因而其具有較為完備的理論基礎(chǔ)，參數(shù)對土壤理化性質(zhì)敏感性很弱［19］，預(yù)測精度高，不過該模型對于砂性土、擾動土壤以及對于高吸力段曲線預(yù)測結(jié)果有待檢驗(yàn)［23］。 Omuto提出的Biexponential 模型有5 個參數(shù)，并將土壤分為結(jié)構(gòu)性孔隙（structural pore-space）和質(zhì)地孔隙（textural pore-space），對應(yīng)不同物理意義的參數(shù)，該模型在全吸力段連續(xù)可導(dǎo)，精度優(yōu)良。 Matlan 等［24］改進(jìn)了適用性廣泛的Gardner 模型，將其發(fā)展為4 個參數(shù)的Modified Gardner（MG）模型，其優(yōu)點(diǎn)是適用土壤類型廣泛且參數(shù)的物理意義較為明確。

1.3.2 粒子群優(yōu)化算法

以往水分特征曲線的擬合多借助Matlab、RETC、Vanfit、HypropFit 等軟件，采用非線性最小二乘回歸或者非線性擬合［26］。 試驗(yàn)中低吸力段和高吸力段的數(shù)據(jù)點(diǎn)較多，為了降低極值對擬合參數(shù)的影響［27］，本研究采用粒子群優(yōu)化算法（Particle Swarm Optimization，PSO）對模型參數(shù)尋優(yōu)。 算法是一種自啟發(fā)式的全局優(yōu)化算法，它通過模擬鳥群個體的覓食行為，同時傳遞信息給其他同伴，來協(xié)助鳥群整體覓食。 該算法的優(yōu)勢在于簡單容易實(shí)現(xiàn)，并且沒有許多參數(shù)的調(diào)節(jié)。 有關(guān)該算法的介紹和應(yīng)用較多，其原理和程序設(shè)計(jì)詳見文獻(xiàn)［28］。

1.3.3 精度評價(jià)標(biāo)準(zhǔn)

均方根誤差（RMSE）為模型預(yù)測值與實(shí)測值方差均值的平方根，均方根誤差越小，表明預(yù)測值越接近實(shí)測值。

式中：N為樣本點(diǎn)數(shù)；ym為實(shí)測值；yp為預(yù)測值。

決定系數(shù)R2反映了預(yù)測值變異在總變異中的比率，可以從整體上反映擬合優(yōu)度，R2越接近于1，表明整體上預(yù)測值越接近于實(shí)測值。

AIC信息準(zhǔn)則是由日本統(tǒng)計(jì)學(xué)家赤池弘次創(chuàng)立的，用于權(quán)衡所估計(jì)模型的精度與復(fù)雜性的指標(biāo)［29］。它不僅考慮模型本身的精度，而且考慮模型參數(shù)的數(shù)量，參數(shù)越多意味著模型越精確，但增加了復(fù)雜度。 式（3）所確定的AIC指數(shù)表明，AIC越小模型越優(yōu)良。AIC的簡化表達(dá)式［30］如下：

式中：np為模型參數(shù)數(shù)量。

2 結(jié)果與討論

2.1 土壤的質(zhì)地分類

由圖2 看出，研究區(qū)土壤有黏土、粉質(zhì)黏土、粉質(zhì)黏壤土、粉質(zhì)壤土、粉土、砂壤土、壤砂土和砂土8 個種類。 這說明研究區(qū)土壤具有較強(qiáng)變異性，土壤質(zhì)地種類繁多，土壤結(jié)構(gòu)變化復(fù)雜。

圖2 研究區(qū)土壤質(zhì)地分類三角圖

2.2 擬合結(jié)果分析

根據(jù)粒子群優(yōu)化算法，對實(shí)測水分特征曲線采用BC 模型、C 模型、VG 模型、FX 模型、Biexponential（雙指數(shù)）模型以及MG 模型分別擬合。 所有模型中，殘余含水率θr和飽和含水率θs的約束條件為

式中：θ2kPa為壓力在 2 kPa 時的含水率；θ1400kPa為壓力在1 400 kPa 時的含水率。

以最小RMSE為目標(biāo)函數(shù)，對擬合參數(shù)進(jìn)行非線性規(guī)劃。 在參數(shù)設(shè)置中，學(xué)習(xí)因子c1＝c2＝2，慣性權(quán)重w＝0.5。 粒子群數(shù)目和迭代次數(shù)影響擬合精度和運(yùn)算速度，經(jīng)反復(fù)調(diào)試，其最優(yōu)設(shè)置見表2。

表2 粒子群優(yōu)化算法的參數(shù)設(shè)置

經(jīng)過粒子群優(yōu)化算法尋優(yōu)，得到總計(jì)233 個樣品的精度評價(jià)情況。 受文章篇幅所限，為突出擬合結(jié)果的差異性，本文選擇8 個樣品的擬合結(jié)果（對應(yīng)8 種質(zhì)地），對6 種模型的擬合結(jié)果進(jìn)行簡要對比分析（見圖3）。

通過圖3 看出，8 種質(zhì)地土壤的VG 模型和MG 模型曲線與實(shí)測曲線幾乎重合，表明這兩種模型具有較高精度，但是模型的表現(xiàn)在不同質(zhì)地和不同吸力段有一定差別。 因此，可以得到結(jié)論：①無論何種質(zhì)地，BC模型和C 模型在低吸力段明顯高估含水率，在中吸力段的擬合效果好于高吸力段。 ②BC 模型和C 模型更適合粉土和壤砂土。 FX 模型的表現(xiàn)不穩(wěn)定，對于絕大部分土壤質(zhì)地， FX 模型表現(xiàn)尚可； 對于粉質(zhì)黏土和砂土，該模型明顯偏離實(shí)測曲線，原因可能是修正系數(shù)的選擇問題導(dǎo)致FX 模型表現(xiàn)不佳。 ③Biexponential模型在8 種質(zhì)地土壤的表現(xiàn)穩(wěn)定，在中吸力段的表現(xiàn)優(yōu)于在低吸力段和高吸力段的表現(xiàn)，尤其對于黏土和粉黏壤土，該模型明顯低估了低吸力段和高吸力段含水率。

2.3 精度比較

圖3 未能定量反映233 個樣品的擬合精度。 因此，根據(jù)表3 對233 個樣品的整體擬合情況進(jìn)行分析，優(yōu)選最佳水分特征曲線模型。 233 個樣品的土壤水分特征曲線擬合精度表顯示，233 個樣品的決定系數(shù)為0.855 9～0.999 8，這表明不同模型的表現(xiàn)有較大差距。通過R2可以看出VG 模型和MG 模型顯示出明顯的優(yōu)越性，F(xiàn)X 模型和 Biexponential 模型次之，C 模型再次之，BC 模型最差。 從RMSE均值來看，VG 模型和 MG模型相似，其余模型顯示了與R2相反的大小關(guān)系，表明RMSE反映出的模型優(yōu)劣與R2是一致的。

模型參數(shù)數(shù)目越少，表明其越具有推廣能力。AIC均值（赤池信息準(zhǔn)則）兼顧了擬合精度與參數(shù)數(shù)目，是綜合評判模型性能的指標(biāo)，已有研究者利用AIC均值評價(jià)了不同水分特征曲線經(jīng)驗(yàn)?zāi)Ｐ偷男阅埽?0-31］，因此本研究選擇AIC均值作為衡量模型優(yōu)劣的最終指標(biāo)。根據(jù)AIC均值越小模型越優(yōu)的原則（見表3），6 種模型的優(yōu)劣為： VG 模型＞MG 模型 ＞FX 模型 ＞Biexponential 模型＞C 模型＞BC 模型。 進(jìn)一步根據(jù)AIC均值大小，可以將6 種模型的性能分為3 類，第一類為VG 模型和MG 模型，具有很高精度；第二類FX 模型和Biexponential 模型具有較高精度；第三類C 模型和BC 模型具有一般精度。 因此，研究區(qū)包氣帶土壤水分特征曲線的最優(yōu)模型為VG 模型，但是VG 模型的AIC均值僅比MG 模型高了2.91%，MG 模型與VG 模型的精度相差無幾，表明MG 模型同樣具有一定的應(yīng)用價(jià)值。

表3 土壤水分特征曲線模型的擬合精度評價(jià)

3 結(jié) 論

為了篩選適用于黃河下游地區(qū)包氣帶土壤的最佳水分特征曲線模型，本研究以典型區(qū)233 個包氣帶土壤樣品的實(shí)測水分特征曲線數(shù)據(jù)，基于粒子群優(yōu)化算法，對比6 種經(jīng)驗(yàn)?zāi)Ｐ偷臄M合結(jié)果，得出以下結(jié)論：

（1）BC 模型和C 模型明顯高估了低吸力段的含水率，相對而言，BC 模型和C 模型更適合粉土和壤砂土；對于粉質(zhì)黏土和砂土，F(xiàn)X 模型表現(xiàn)不佳；Biexponential 模型在8 種質(zhì)地土壤的表現(xiàn)穩(wěn)定，對于黏土和粉黏壤土，Biexponential 模型明顯低估了低吸力段和高吸力段含水率。

（2）以赤池信息準(zhǔn)則AIC均值為最優(yōu)模型的評選指標(biāo)，結(jié)論表明6 種模型的優(yōu)劣程度為：VG 模型＞MG模型＞FX 模型＞Biexponential 模型＞C 模型＞BC 模型。因此，適宜研究區(qū)包氣帶土壤樣品的最佳水分特征曲線為VG 模型。

（3）本研究的經(jīng)驗(yàn)?zāi)Ｐ瓦x擇存在一定主觀性，擬合結(jié)果的差異性表明經(jīng)驗(yàn)?zāi)Ｐ团c不同土壤質(zhì)地和吸力段之間的關(guān)系仍需進(jìn)一步研究。