兩種模型對土壤水分特征曲線擬合的比較分析

朱蔚利，肖自幸，牛健植，邵文偉，張由松，李 想，武曉麗，趙玉麗

（北京林業(yè)大學(xué)水土保持與荒漠化防治教育部重點實驗室，北京林業(yè)大學(xué)水土保持學(xué)院，北京 100083）

土壤水分特征曲線是描述土壤水的基質(zhì)勢或土壤水吸力隨土壤含水率變化的關(guān)系曲線，可以表示出土壤水的能量和數(shù)量之間的關(guān)系[1-2]，是研究土壤水分的保持和運動所用到的反映土壤水分基本特征的曲線，在土壤研究中有重要的實用價值：一是可利用他進(jìn)行土壤水吸力和含水率之間的換算；二是可以間接地反映出土壤孔隙大小的分布；三是可用來分析不同質(zhì)地土壤的持水性和土壤水分的有效性；四是應(yīng)用數(shù)學(xué)物理方法對土壤中水運動進(jìn)行定量分析的必不可少的重要參數(shù)[2]。近年來，土壤水分特征曲線在土壤水管理、土壤-植物-大氣連續(xù)體系研究中也有應(yīng)用，且具有重要意義[3-5]。

目前土壤水分特征曲線主要通過直接測定法和間接推導(dǎo)法獲得，直接測定法如張力計法、離心機法、壓力膜儀法等。但是，土壤水分特征曲線的影響因素復(fù)雜，僅憑這些實測數(shù)據(jù)還難以從理論上得到水吸力與土壤含水率的確切關(guān)系式。于是，在大量試驗的基礎(chǔ)上，人們提出一些經(jīng)驗公式來擬合實測數(shù)據(jù)[6]，其中常見的有：Broods-Corey模型[7]，Gardner模型[8]，Campbell模型[9]，van-Genuchten模型[10]和Gardner-Russo模型[11]等。Van-Genuchten模型由于適用土壤質(zhì)地范圍廣[12]，其線型與實測數(shù)據(jù)的擬合程度好而被廣泛應(yīng)用；土壤水分單一參數(shù)模型則是Gregson等依據(jù)模型參數(shù)的內(nèi)在關(guān)系建立的[13]。Williams利用該模型只有一個參數(shù)的特點，提出了通過測定土壤水吸力與土壤含水量關(guān)系曲線上一個點的數(shù)據(jù)來計算整個土壤水分特征曲線[14]。本文用Van-Genuchten模型和單一參數(shù)模型擬合鷲峰國家森林公園4個土樣的實測土壤水分特征數(shù)據(jù)，并進(jìn)行比較，為模型的選用提供參考。

1 材料與方法

1.1 試驗材料

試驗所用土樣取自北京市海淀區(qū)北安河鄉(xiāng)鷲峰國家森林公園，在山腳與山腰各選一個樣地，每個樣地的取樣土層為0～20 cm、20～40 cm，每層取兩個土樣，作為重復(fù)。試驗所用離心機為落地式，土壤用高速冷凍離心機H-1400 pF，與該離心機配套的環(huán)刀的內(nèi)直徑和高度均為5 cm。

1.2 采樣方法

在野外用與離心機配套的環(huán)刀取原狀土樣，將土樣帶回實驗室在蒸餾水中浸泡48 h，使其達(dá)到飽和狀態(tài)，取出稱重。然后用離心機對土樣進(jìn)行脫濕處理，從小到大共設(shè)置了10個不同轉(zhuǎn)速，分別為500、1 000、1 500、2 000、2 500、3 000、4 000、5 000、6 000、7 000 r/min，每次旋轉(zhuǎn)完成后都進(jìn)行稱重。根據(jù)所得數(shù)據(jù)可計算出不同轉(zhuǎn)速下所對應(yīng)的土壤體積含水率，取兩個重復(fù)土樣的平均值作為最終土壤體積含水率，從而得到實測的土壤水分特征曲線。

1.3 研究方法

1.3.1 Van-Genuchten模型 Van-Genuchten（1980）模型既連續(xù)又有連續(xù)斜率，得到的曲線光滑，對絕大多數(shù)土壤在相當(dāng)寬的水勢或含水量范圍內(nèi)具有普遍適用性，并可得到相對導(dǎo)水率的解析解，因而在土壤物理領(lǐng)域得到了最為廣泛的應(yīng)用，其公式表示為：

式中：θ為體積含水量（cm3/cm3）；θr為殘留含水量（cm3/cm3）；θs為飽和含水量（cm3/cm3）；h為吸力(cm H2O)，當(dāng)土壤處于飽和狀態(tài)時，可認(rèn)為土壤吸力h等于零；a、n、m為擬合參數(shù)。其中參數(shù)a為尺度函數(shù)，與平均孔隙直徑成反比，a=1/hb，hb為進(jìn)氣吸力，在飽和土壤中施加吸力，當(dāng)吸力較小時，土壤中尚無水排出，土壤含水率維持飽和值，但當(dāng)吸力增加超過某一臨界值時，土壤孔隙中的水分開始排出，該臨界負(fù)壓值就稱為進(jìn)氣值，即土壤水由飽和轉(zhuǎn)變?yōu)榉秋柡蜁r的負(fù)壓值。不同質(zhì)地土壤的進(jìn)氣值不同，一般來說，重質(zhì)粘性土壤進(jìn)氣值較大；輕質(zhì)土或結(jié)構(gòu)良好的土壤進(jìn)氣值較小，或者沒有明顯的進(jìn)氣值，因此很難直接從土壤水分特征曲線圖上獲得a[6]。參數(shù)n的大小則決定著土壤水分特征曲線的坡度，當(dāng)n較大時，曲線較緩，當(dāng)n較小時，曲線較陡。參數(shù)m可以用參數(shù)n表達(dá)，m=1-1/n（n>1）。

1.3.2 單一參數(shù)模型 Gardner模型是土壤水分特征曲線模型中參數(shù)比較少、形式比較簡單的模型，應(yīng)用也比較廣泛，其雙對數(shù)形式為：

式中：h為負(fù)壓水頭（kPa）；θ為體積含水量（cm3/cm3）；A、B為大于0的經(jīng)驗參數(shù)。且A和B之間具有線性關(guān)系：

式中p、q為常數(shù)。

將式（3）代入式（2）得：

式（4）中只有一個變量B，因此稱作土壤水分特征曲線的單一參數(shù)模型[15-16]。不同土壤的水分特征曲線B值不同，則不同的B值就能表示不同的土壤水分特征曲線，所以應(yīng)用該參數(shù)可以分析土壤水分特征曲線的空間變異性。

1.4 數(shù)據(jù)處理

主要用Excel2010處理基礎(chǔ)數(shù)據(jù)，用專業(yè)函數(shù)繪圖軟件Origin 8.1對Van-Genuchten模型和土壤水分特征曲線單一參數(shù)模型進(jìn)行擬合，并求相關(guān)參數(shù)值。

2 結(jié)果與分析

2.1 實測的土壤水分特征曲線

根據(jù)試驗結(jié)果，做出了4條土壤水分特征曲線，如圖1所示。為了便于說明，分別用A、B、C、D來表示：A表示山腳樣地的0～20 cm土樣，B表示山腳樣地的20～40 cm土樣，C表示山腰樣地的0～20 cm土樣，D表示山腰樣地的20～40 cm土樣。由于圖1（a）中1 000 cm H2O吸力一些變化情況表現(xiàn)不明顯，故特將該部分做成圖1（b），從圖1（b）中可看到B曲線呈“S”型。由圖1可知，山腳土樣的土壤飽和含水率大于山腰的土壤飽和含水率，且0～20 cm土層土壤飽和含水率大于20～40 cm土層的土壤飽和含水率。其原因可能是山腳土壤發(fā)育較山腰厚，發(fā)育更為成熟，腐殖質(zhì)等有機質(zhì)含量較高，并且根系、枯枝落葉等對0～20 cm土層土壤結(jié)構(gòu)的改良作用大于對20～40 cm土層土壤結(jié)構(gòu)的改良作用，使其持水性較好。

2.2 Van-Genuchten擬合模型

圖1 土壤水分特征曲線

將實測土壤含水率與對應(yīng)的土壤水吸力輸入origin8.1軟件，設(shè)置擬合參數(shù)的初始值為：θs=0.5，θr=0.14，a=0.08，n=3。擬合結(jié)果見表1，可見擬合曲線的相關(guān)系數(shù)均大于0.99，說明擬合的精度很高。

表1 土壤水分特征曲線Van-Genuchten模型擬合參數(shù)

將擬合參數(shù)帶入Van-Genuchten模型，計算出10個土壤水吸力值所對應(yīng)的土壤含水率，作為擬合含水率，與實測含水率進(jìn)行比較，如圖2（a）、（b）、（c）、（d）所示。由圖2可知，絕大部分實測點落在擬合曲線上，近于重合。為了進(jìn)一步比較擬合值與實測值的差別，計算了它們的相關(guān)系數(shù)和均方根誤差，見表2。均方根誤差公式[17]見公式（5）。

圖2 Van-Genuchten模型模擬土壤水分特征曲線與實測土壤水分特征曲線對比

表2 Van-Genuchten模型含水率模擬值與實測值的比較

式中：RMSE為均方根誤差；N為實測含水率的個數(shù)；θmi為實測土壤含水率；θei為由Van-Genuchten模型計算得到的土壤含水率。

由表2知，Van-Genuchten模型含水率模擬值與實測值的相關(guān)率都大于0.995，均方根誤差都小于0.01，擬合效果很好，所以可用Van-Genuchten模型參數(shù)的物理含義來描述實際的土壤水分特征曲線的特征，這也正是Van-Genuchten模型的價值體現(xiàn)。由表1可知：所測原狀土樣的飽和含水率θs基本在40%～60%之間，山腳的土壤含水率比山腰的高出10%，可能是因為山腳的土層比山腰的厚，土壤發(fā)育更為成熟，孔隙結(jié)構(gòu)良好，所以持水能力更強。所有土樣的殘留含水率θr都接近10%，該模型中殘留含水率θr指土壤水分特征曲線導(dǎo)數(shù)等于0時的土壤含水率，即dθ/dh=0時的土壤含水率。在實際應(yīng)用中，一般取凋萎點處的含水率，即用凋萎系數(shù)作為殘留含水率。由此可知所測土樣的凋萎系數(shù)約為0.1；參數(shù)a一般認(rèn)為是進(jìn)氣值ha的倒數(shù)，由a計算ha得到A、B、C、D四條曲線的進(jìn)氣值，分別為6.35、39.68、37.51、9.19 cm H2O；參數(shù)n的大小表示著土壤水分特征曲線的坡度。由表1和圖1可看出，曲線B的n值最大，曲線也最緩，曲線D的n值最小，其曲線也最陡。

2.3 單一參數(shù)擬合模型

為了求得土壤特征單一參數(shù)模型的常數(shù)p、q的值，先利用Gardner模型的雙對數(shù)形式對實測數(shù)據(jù)進(jìn)行擬合，以便得到A、B值，擬合結(jié)果見表3。

表3 土壤水分特征曲線Gardner對數(shù)模型的擬合參數(shù)

由表3知，擬合結(jié)果的相關(guān)系數(shù)均大于0.9，擬合效果較好。A、B之間呈線性關(guān)系，故將表3中的4對值進(jìn)行擬合，得到p=2.227 58，q=1.618 15，所以A=2.227 58 B+1.618 15，R2=0.981 59。將p、q值代入式（4），得到所取土樣的土壤水分特征曲線單一參數(shù)模型：

在4組實測土樣中各隨機選取一對值（土壤水吸力及對應(yīng)土壤含水率）代入單一參數(shù)模型，求出B值，并以此計算各個水吸力下對應(yīng)的土壤含水率，從而繪制土壤水分特征曲線，并與實測土壤水分特征曲線作對比，如圖3（a）、（b）、（c）、（d）。

圖3 單一參數(shù)模型模擬土壤水分特征曲線與實測土壤水分特征曲線對比

由圖3可以看出，單一參數(shù)模型模擬土壤水分特征曲線與實測土壤水分特征曲線有一些差別。為了確定這個差別對單一參數(shù)模型應(yīng)用的影響，計算了模擬值與實測值的相關(guān)系數(shù)及均方根誤差，見表4。由表4可知，模擬值與實測值的相關(guān)系數(shù)均在0.96以上，均方根誤差在0.01～0.04之間，誤差較小，可用于實踐。宋孝玉等[16]測得的單一參數(shù)模型土壤含水率模擬值與實測值的相關(guān)系數(shù)為0.957 4～0.992 1，均方根誤差為0.017 5～0.052 3，其土壤水分特征曲線單一參數(shù)模型模擬結(jié)果在甘肅省慶陽市境內(nèi)的南小河溝流域不同地貌以及不同土壤質(zhì)地的應(yīng)用中效果良好。

單一參數(shù)模型中的B參數(shù)，沒有很明確的物理意義，只是不同的土樣有不同的B值，反映了土壤水分特征曲線的空間變異情況。

2.4 Van-Genuchten模型與單一參數(shù)模型模擬精度的比較

比較表2和表4可知，Van-Genuchten模型與實測值的相關(guān)性大于單一參數(shù)模型，且Van-Genuchten模型的均方根誤差比單一參數(shù)模型小了一個數(shù)量級，說明Van-Genuchten模型對實測值模擬的精度遠(yuǎn)遠(yuǎn)高于單一參數(shù)模型。這跟它們的模擬原理有關(guān)，Van-Genuchten模型是直接對實測值進(jìn)行擬合，而單一參數(shù)模型則是通過測若干土壤樣點，先得到Gardner的對數(shù)模型參數(shù)A、B，再由A、B進(jìn)行線性擬合，求出p、q，在隨機選取一對土壤水吸力以及對應(yīng)的土壤含水率求出參數(shù)B，然后才進(jìn)行擬合，由于每一步擬合都可能會有誤差，從而導(dǎo)致了誤差的累積，使得其模擬精度下降。

表4 單一參數(shù)模型土壤含水率模擬值與實測值比較

3 結(jié)論

用離心機法對從鷲峰國家森林公園選取的4個土樣進(jìn)行脫濕處理，得到這4個土樣的土壤水分特征曲線。由實測土壤水分特征曲線可知，土壤水分特征曲線可以比較準(zhǔn)確地反映出土壤結(jié)構(gòu)的優(yōu)劣、植物根系及枯枝落葉對0～20 cm土層土壤結(jié)構(gòu)的改良效果優(yōu)于20～40 cm土層土壤結(jié)構(gòu)，從而使其土壤持水性較好。并用Van-Genuchten模型和土壤水分特征曲線單一參數(shù)模型進(jìn)行擬合，通過對擬合結(jié)果的分析比較，得出如下結(jié)論：

（1）從模型參數(shù)包含的物理意義來看，Van-G enuchten模型的參數(shù)可以反映土樣的飽和含水率、殘留含水率及進(jìn)氣吸力，具有直觀的物理意義；而單一參數(shù)模型中的參數(shù)表示空間變異情況，比較抽象。

（2）從工作量大小來看，Van-Genuchten模型需要測定每個土樣的若干實際含水率，并且實測值不能太少，否則影響精度；而單一參數(shù)模型每個土樣只需知道一個土壤水分特征曲線上的點就可以計算出參數(shù)B的值，從而繪出土壤水分特征曲線，省時省力，尤其對某一大型樣地進(jìn)行大規(guī)模土壤水分特征曲線測定時，優(yōu)勢更加明顯。

（3）從模型應(yīng)用范圍來看，Van-Genuchten模型由于參數(shù)較多，對鷲峰國家森林公園土樣的水分特征曲線擬合結(jié)果與實測結(jié)果吻合度較高，能用于高精度要求的情況；而單一參數(shù)模型與實測結(jié)果的吻合度相對較低，更適用于精度要求稍低的情況。

[1] 雷志棟，楊詩秀，謝森傳.土壤水動力學(xué)[M].北京：清華大學(xué)出版社，1988.18-24.

[2] 黃昌勇.土壤學(xué)[M].北京：中國農(nóng)業(yè)出版社，2000.108-110.

[3] 孫 迪，夏靜芳，關(guān)德新，等.長白山闊葉紅松林不同深度土壤水分特征曲線[J].應(yīng)用生態(tài)學(xué)報，2010，21（6）：1405-1409.

[4] 楊啟紅，陳麗華，王 宇.基于采樣密度不同土壤的水分特征參數(shù)的預(yù)測及誤差分析[J].安徽農(nóng)業(yè)科學(xué)，2009，（27）：13189-13191.

[5] 李小剛.影響土壤水分特征曲線的因素[J].甘肅農(nóng)業(yè)大學(xué)學(xué)報，1994，29（3）：273-278.

[6] 彭建平，邵愛軍.用MatLab確定土壤水分特征曲線參數(shù)[J].土壤，2007，39（3）：433-438.

[7] 楊靖宇，屈忠義.河套灌區(qū)區(qū)域土壤水分特征曲線模型的確定與評價[J].干旱區(qū)資源與環(huán)境，2008，22（5）：155-159.

[8] Gardner W R,Hillel D,Benyamini Y.Post irrigation movement of soil water：I.Redistribution[J].Water Resources Research,1970,6（3）:851-861.

[9] Campell G S.A simple method for determining unsaturated conductivity from moisture retention data[J].Soil Science,1974,117:311-314.

[10]Van Genuchten M Th.A closed-form equation for predicting the hydraulic conductivity of unsaturated soils[J].Soil Science Society of America Journal,1980,44:892-898.

[11]Russo D.Determining soil hydraulic properties by parameter estimation:on the selection of a model for the hydraulic properties[J].Water Resour,1988,24:453-459.

[12]來劍斌，王全九.土壤水分特征曲線模型比較分析[J].水土保持學(xué)報，2003，17(1)：137-140.

[13]Gregson K,Hector D J,McGowan M.A one-parameter model for the soil water characteristic[J].Journal of Soil Science,1987,38:483-486.

[14]Williams R D,Ahujia L R.Evaluation of similar-media scaling and a one-parameter model for estimating the soil the water characteristic[J].Journal of Soil Science,1992,43:237-248.

[15]賈宏偉，康紹忠，張富倉.土壤水力參數(shù)的單一參數(shù)模型[J].水利學(xué)報，2006，37（3）：272－277.

[16]宋孝玉，李亞娟，李懷有，等.土壤水分特征曲線單一參數(shù)模型的建立及應(yīng)用[J].農(nóng)業(yè)工程學(xué)報，2008，24（12）：12-15.

[17]Willmott C J,Ackleson S G,Davis R E,et al.Statistics for the evaluation and comparison of models[J].Journal of Geophysical Research,1985,90:8995-9005.