羅小清 楊屹杰 唐輝亮

摘要 研究了土壤水分特征曲線對探討土壤水分運動機理及溶質(zhì)運移動力學機制具有重要科學意義。由于測定土壤水分特征曲線的長歷時性和時間限制，在從相關參考文獻中提取土壤水分特征曲線測定數(shù)據(jù)的基礎上，建立土壤水分特征曲線的資料庫，采用van Genuchten模型對土壤水分特征曲線進行描述，分別應用vanfit軟件和MATLAB軟件下的lsqcurvefit函數(shù)2種方法對土壤進行van Genuchten模型下的4個參數(shù)擬合分析。結(jié)果表明，2個程序擬合結(jié)果大體相近且優(yōu)良，均方差較小，擬合曲線大致相同，且土壤質(zhì)地不同，其土壤水分特征曲線差別較大。

關鍵詞 土壤水分特征曲線；van Genuchten模型；參數(shù)擬合；vanfit軟件；lsqcurvefit函數(shù)

中圖分類號：S152.7 文獻標識碼：B 文章編號：2095–3305（2023）07–0254-03

土壤水分特征曲線是描述土壤水基質(zhì)勢或土壤吸力隨土壤含水量變化關系的曲線，用于表征土壤水分在非飽和狀態(tài)下能態(tài)與數(shù)量之間的關系，反映土壤持水和釋水特性，可作為最重要和最基本的工具，用于研究土壤水分運動，調(diào)節(jié)利用土壤水，改良土壤等，該曲線的基本運動趨勢為：隨著含水量的減少，土壤基質(zhì)勢降低或水吸力升高。受到土壤在吸水和釋水過程中土壤空氣的作用和固、液而接觸角不同的影響，土壤水分特征曲線并不是單值函數(shù)曲線。土壤物理性質(zhì)、土溫、土壤中有機質(zhì)含量會對土壤水分特征曲線造成一定的影響，土壤種植作物也會影響土壤物理性質(zhì)及有機質(zhì)含量，不同的作物，其耕作制度和施肥方法也存在顯著的差異，且還會影響土壤的性狀，因此相同層次的不同作物表現(xiàn)出的土壤水分特征曲線也有明顯的不同[1]。

1 資料收集及資料庫的建立

鑒于測定土壤水分特征曲線實驗耗時較長，擬從有測定土壤水分特征曲線的參考文獻中收集相關資料，提取相關數(shù)據(jù)，建立土壤水分特征曲線資料庫。通過資料的收集和數(shù)據(jù)的提取整理，得到相關資料庫（表1）。

2 土壤水分特征曲線的模擬

2.1 土壤水分特征曲線模型選定

土壤持水曲線是研究土壤水動力學性質(zhì)必不可少的重要參數(shù)。土壤水分特征曲線數(shù)學模型種類較多，較常應用的是Brooks-Corey模型、van?Genuchten模型、Gardner-Simmons模型、Broadbridge-White模型。大量研究得出，在眾多構(gòu)建出的數(shù)學模型中，van Genuchten模型以其與實測數(shù)據(jù)擬合程度最好。

Van Genuchten（1980）模型具有連續(xù)特點且有連續(xù)斜率，得到的曲線光滑，對大多數(shù)水勢或含水量范圍較寬的土壤具有普遍適用性，并可得到相對導水率的解析解，因此被廣泛應用于土壤物理領域。該公式為：

式中，θ表示體積含水量（cm3/cm3）；θr指殘留含水量（cm3/cm3）；θs為飽和含水量（cm3/cm3）；h為吸力（cm H2O）；

土壤呈飽和狀態(tài)時，可認為土壤吸力h為0；a、n、m為擬合參數(shù)。參數(shù)a是尺度函數(shù)，與平均孔隙直徑呈反比，a=1/hb，hb為進氣吸力，向飽和的土壤中施加吸力，吸力小時，土壤中無水排出，土壤含水率處于飽和值，若吸力增加至某一臨界值以上時，土壤孔隙中排出水分，該臨界負壓值即為進氣值，土壤水由飽和轉(zhuǎn)變?yōu)榉秋柡蜁r的負壓值。不同質(zhì)地土壤進氣值存在差異，重質(zhì)粘性土壤通常進氣值較大，輕質(zhì)土或良好結(jié)構(gòu)的土壤則進氣值小，或無明顯進氣值，不易直接從土壤水分特征曲線圖上獲取a[2]。參數(shù)n的大小表征土壤水分特征曲線的坡度，當n比較大時曲線較緩，當n比較小時則曲線較陡。參數(shù)m可由參數(shù)n表達，m=1-1/n （n>1）。

2.2 van Genuchten模型中參數(shù)的獲取

運用van Genuchten模型的關鍵是如何獲得其中的4個參數(shù)。本研究采用2種不同方法獲取van Genuchten模型中的參數(shù)。

第一，利用vanfit軟件獲取van Ge-nuchten模型中的參數(shù)。vanfit軟件是中文軟件，根據(jù)壓力膜儀實驗數(shù)據(jù)，即可計算得出土壤持水曲線參數(shù)。

利用vanfit軟件擬合土壤水分特征曲線參數(shù)的方法如下：

從上述步驟得到的Excel表格中提取土壤吸力與體積含水量數(shù)據(jù)。由于vanfit擬合軟件內(nèi)置的土壤吸力為土壤負壓h，單位cm水柱；體積含水率為土壤含水量θ，單位g/g。因此有些土壤吸力以Pa、kPa等表示的數(shù)據(jù)需要單位換算：100 cm水柱=10 kPa=10 000 Pa。體積含水率也以小數(shù)形式表示。

根據(jù)實測點數(shù)目，在van格式文件中輸入土樣數(shù)據(jù)，如：

6（實測點數(shù)）

1.0 30.7 83.0 220.1 493.3 691.4

（土壤吸力）

0.313 0.242 0.197 0.178 0.160 0.146

Q（體積含水率）

然后在軟件中打開van數(shù)據(jù)文件，根據(jù)土壤性質(zhì)在軟件迭代參數(shù)框中輸入?yún)?shù)初始值后，“擬合”“圖形”選項從灰選狀態(tài)變成可選狀態(tài)，點擊“擬合”，可得到一組模型參數(shù)；點擊“圖形”，得到2個擬合曲線結(jié)果（圖1）。

左圖為一般坐標下擬合的土壤水分特征曲線，右圖為半對數(shù)（橫軸為對數(shù)坐標）下擬合的土壤水分特征曲線。若擬合結(jié)果較好，則取軟件所求得擬合結(jié)果參數(shù)；否則進行參數(shù)迭代，直至最優(yōu)。

第二，利用MATLAB軟件lsqcurvefit函數(shù)獲取van Genuchten模型中的參數(shù)。lsqcurvefit函數(shù)擬合步驟如下：

（1）建立fun1.m文件：

function ydata=fun1（p，xdata）

ydata=p（1）+[p（2）-p（1）]./{1+

[p（3）*xdata].^p（4）}.^[1-1./p（4）]；

% xdata，ydata為實驗所得數(shù)據(jù)；p

（ ）為土壤持水曲線參數(shù)。

（2）在Command Window輸入初始數(shù)據(jù)：

clear

%輸入實驗數(shù)據(jù) xdata，ydata

xdata = [實驗數(shù)據(jù)] ；

ydata = [實驗數(shù)據(jù)] ；

%設定迭代初始值Parameters0，分別代表θr，θs，α，n

Parameters0 =[0.1，0.1，0.01，1]；

（3）輸入初始數(shù)據(jù)后，代入lsqcurvefit

擬合函數(shù)擬合，繼續(xù)Command Window 輸入：

[Parameters， resnorm]=lsqcurvefit （@fun1，Parameters0，xdata，ydata）

% Parameters 為擬合的土壤持水曲線參數(shù)

% resnorm 為均方差

（4）在工作區(qū)內(nèi)，便會顯示出土壤持水曲線參數(shù)Parameters的擬合結(jié)果，使用非常方便。

（5）如果之后要顯示擬合圖形，首先要根據(jù)模擬參數(shù)計算出模擬曲線的若干點的坐標值，其命令如下：

x1 = [0：0.1：1000]；

y1=Parameters（1）+（Parameters（2）-Parameters（1））./（1+（Parameters（3）*x1）.^Parameters（4））.^（1-1./Parameters（4））;

如要顯示一般坐標土壤持水曲線，可直接輸入命令：

plot（xdata，ydata，ko，x1，y1，k-）

如要顯示半對數(shù)坐標土壤持水曲線，可直接輸入命令：

semilogx （xdata，ydata，ko，x1，y1，k-）

并根據(jù)需要，在圖上添加標簽及說明。

xlabel（‘土壤負壓h / cm 水柱）；

ylabel（‘土壤含水量θ / g.g-1）；

legend（‘實測值，擬合曲線）；

由于參數(shù)直接由程序計算得到，因此不需要再迭代。

2.3 土壤水分特征曲線參數(shù)擬合結(jié)果

用Vanfit軟件和MATLAB軟件按照上述方法，對土壤水分特征曲線資料庫中的土壤水吸力與土壤含水率關系數(shù)據(jù)采用van-Genuehten模型進行土壤水分特征曲線擬合，得到最終擬合效果較好的模型參數(shù)和土壤水分特征曲線圖（表2）。

3 擬合結(jié)果分析與展望

從各擬合結(jié)果來看，2個程序擬合結(jié)果大體上是相近且優(yōu)良的，4個參數(shù)值相差不大，均方差也很小，擬合的曲線大致相同。

對于砂壤，2個軟件所求出的4個參數(shù)相差較大，由于土壤吸力與土壤含水率對應關系特征不明顯，Vanfit軟件與MATLAB軟件的擬合結(jié)果均不是很好，2個軟件擬合的曲線都有較多未經(jīng)過的實測點。

而粉黏壤，2個軟件擬合的曲線有一定偏差，Vanfit軟件擬合的曲線會經(jīng)過更多的實測點，因而曲線稍向上凸；而MATLAB擬合的曲線大致走向沒問題但避開了大多數(shù)實測點。

可以發(fā)現(xiàn)，當實測土壤的吸力與含水率的對應關系明顯時，Vanfit軟件與MATLAB軟件對Van- Genuchten模型下4個參數(shù)的率定和土壤水分特征曲線的擬合結(jié)果都是很接近且準確的。

研究土壤水分特征曲線對探討土壤水分運動機理及溶質(zhì)運移動力學機制具有重要的科學意義，許多研究表明，準確的土壤性質(zhì)數(shù)據(jù)及其土壤水分特征曲線擬合的優(yōu)良程度是研究的關鍵因素，而文獻中建立的資料庫所收集的土壤樣本不多，擬合曲線僅采用了vanfit軟件與MATLAB軟件lsqcurvefit函數(shù)2種方法。因此在今后研究中需要建立起準確、可靠和全面的資料庫，以便能對土壤水分特征曲線有更深入研究，并采用更多方法進行擬合。

參考文獻

[1] 郭鳳臺，遲藝俠，程東娟，等.土壤水分特征曲線試驗研究[J].南水北調(diào)與水利科技，2006（2）：47-48.

[2] 張強，孫向陽，黃利江，等.毛烏素沙地土壤水分特征曲線和入滲性能的研究[J].林業(yè)科學研究，2004（S1）：9-14.

Simulation of Soil Moisture Characteristic Curve Based on van Genuchten Model by vanfit and MATLAB Software

Luo Xiao-qing et al（Jiangle County Meteorological Bureau of Fujian Province， Jiangle， Fujian 353300）

Abstract The study of soil water characteristic curve is of great scientific significance for the discussion of soil water movement mechanism and solute transport mechanics mechanism. Due to the long duration and time limit of measuring the soil moisture characteristic curve， established a database of soil moisture characteristic curve on the basis of extracting the measured data of soil moisture characteristic curve from relevant references， and used van Genuchten model to describe the soil moisture characteristic curve， The four parameters of the van Genuchten model were fitted and analyzed by using the two methods of vanfit software and lsqcurvefit function under MATLAB software. The results showed that the fitting results of the two programs were generally similar and excellent， and the mean square deviation was also very small. The fitting curves were basically the same. It is found that the soil moisture characteristic curves vary greatly with different soil textures.

Key words Soil water characteristic curve; van Genuchten model; Parameter fitting; vanfit software; Lsqcurvefit function