郭 華，樊貴盛

(太原理工大學水利科學與工程學院，太原 030024)

土壤水分入滲過程[1]是指降水或灌溉水通過地表進入土壤的過程，是田間水循環(huán)的重要組成部分，也是灌水技術參數優(yōu)化、解決地下水補給等問題的重要依據。為了反映土壤水分入滲過程，國、內外學者提出了許多入滲模型[2-4]，如Horton模型、Smith模型、蔣定生公式等，其中Horton模型和Kostiakov模型是早期提出的經驗型模型，但沒有明確的物理意義；Green-Ampt模型有明確的物理基礎，但對參數精度要求低；而Philip入滲模型不僅有明確的物理基礎，對參數精度要求高，而且能較好的反映土壤水分入滲過程。邵明安等人[5]曾用Philip模型預測了土壤水分運動參數，武世量[6]研究了Philip公式標定因子與土壤理化參數的關系，任尚崗等人[7]基于Philip公式分析了單環(huán)入滲的三維動態(tài)過程，原林虎[8]采用Philip模型描述了水分入滲過程，但目前而言，對Philip入滲模型參數的預報尚鮮見。

土壤水分入滲參數[9]是表征土壤入滲能力的指標，直接影響著大田土壤的灌水質量和灌溉效果，然而由于土壤條件空間變異性和其他限制因素的存在，很難通過實測獲得大量的入滲參數?；诖?，學者們通過土壤傳遞函數法[10](Pedo-Transfer Functions，PTFs法)，即建立試驗較易獲得的土壤理化參數(如土壤含水率、干密度、質地等)與入滲參數之間的函數關系，來間接獲得入滲參數，以解決獲取入滲參數困難這一問題。其中馮錦萍[11]用線性回歸分析法對入滲參數進行了預報，其方法簡單，但預測精度低；原林虎[8]采用人工神經網絡法對入滲參數進行了預測，其精度較高，但易出現局部極小等不良現象，不能保證全局最優(yōu)。因此，本研究基于土壤傳遞函數法，以黃土高原區(qū)大田耕作土壤為試驗對象，將試驗較易獲得的土壤體積含水率、干密度、粉粒含量、黏粒含量和有機質含量等土壤理化參數作為輸入變量，建立土壤水分入滲參數為輸出變量的土壤傳遞函數，以實現對入滲參數的多元非線性預報，從而提高預測精度，達到全局最優(yōu)。

1 材料與方法

1.1 試驗條件

試驗對黃土高原區(qū)大田耕作土壤進行，主要在大同、呂梁、應縣等地進行。大同位于山西最北端，平均海拔700～1 400 m，干寒多風，溫差較大，由于多期地貌構造的變動，形成了山地、丘陵、平川等地形；呂梁地處山西中部西側，平均海拔1 000～2 000 m，四季分明，降雨量分布不均，由于呂梁山脈的存在，形成了由南到北，由高山逐漸變?yōu)槠酱ǖ牡貏菪螒B(tài)；應縣位于山西北部的桑干河中游，海拔在1 000～2 300 m之間，全年氣候偏低，其黃土丘陵區(qū)主要集中于東北、西北地帶，其余為平川地區(qū)。試驗區(qū)均位于山西境內，是被黃土廣泛覆蓋的山地形高原，屬大陸性半干旱季風氣候，在自然作用和長期農業(yè)活動的共同影響下，土壤條件在空間上發(fā)生很大變化，形成了復雜多樣的土壤類型。試驗區(qū)有壤土、沙壤土、黏壤土、粉沙質黏壤土、粉沙壤土、黏沙壤土、輕沙壤土等土壤類型，其中耕作層土壤的沙粒含量為9.7%～88%，粉粒含量為10.88%～76%，黏粒含量為0.84%～23.56%；犁底層土壤沙粒、粉粒、黏粒含量分別變化在9.66%～93%、5.41%～85%、0.91%～33.15%。耕作層土壤體積含水率為2.08%～70.67%，有機質含量為0.416～37.54 g/kg，其中地表以下0～10 cm土層的干密度的范圍是0.854～1.426 g/cm3，地表以下10～20 cm土壤干密度變化在0.891～1.523 g/cm3；犁底層土壤體積含水率為3.03%～75.87%，其土壤干密度為0.774～1.762 g/cm3。

1.2 試驗設備與方法

試驗在非凍融期進行。試驗前將內、外環(huán)直徑分別為26、64.4 cm，內、外環(huán)高度均為25 cm的雙套環(huán)入滲儀埋于試驗地塊，試驗時對內、外環(huán)均進行灌水，并用自制的水位控制裝置將內、外環(huán)水位控制在2～3 cm，以消除積水水頭對勢梯度的影響；根據Philip入滲公式，通過定量描述累積入滲水量與灌水時間的關系，便可間接獲得土壤水分入滲參數。大量試驗資料表明，土壤水分達到穩(wěn)滲狀態(tài)的時間在60 min，考慮到試驗資料的可靠性和一致性，將水分穩(wěn)定入滲時間定為90 min。

在大田土壤水分入滲試驗進行的同時，還需對試驗地塊的土壤理化參數進行測定。土壤理化參數包括土壤體積含水率、干密度、質地以及有機質，其中土壤體積含水率采用烘干、稱重的方法測定；土壤干密度用環(huán)刀法獲得；通過篩分曲線法分析各個土層的顆粒級配，從而獲得各層土壤質地；土壤有機質主要存在于耕作層(0～20 cm)土壤中，通過重鉻酸鉀滴定法獲得。

1.3 Philip入滲模型

在對Richard方程系統(tǒng)研究的基礎上，Philip[12]提出了適合土壤一維垂直入滲的入滲模型，該模型形式簡單，計算方便，參數容易確定，具有明確的物理意義，且能較好反映土壤水分入滲過程，在多年的土壤入滲研究中應用廣泛，其具體表達式為：

I=St0.5+At

(1)

式中：I為單位面積上的累積入滲量，cm；t為入滲歷時，min；S為吸濕率，在數值上等于第一個單位時段末的累積入滲量，決定著入滲初期入滲率的大小，cm/min0.5；A為穩(wěn)定入滲率,接近土壤的水力傳導度， cm/min。

2 土壤入滲參數的預報模型

2.1 輸出因子與輸入因子

隨著農業(yè)用水的浪費、灌溉水利用率低、土壤環(huán)境退化等問題的出現，如何合理灌溉成為需要急切解決的問題。為此，學者們提出了土壤水分入滲模型，其中土壤水分入滲參數作為地面灌溉水流運動模擬模型的輸入參數，其取值的合適與否，直接決定著灌水質量和灌溉效果，但通過實測獲取大量的入滲參數相當困難，若建立以試驗較易獲得的土壤理化參數為輸入因子，以入滲參數為輸出因子的土壤傳遞函數，便可方便獲得入滲參數S和A。因此，本文著力研究Philip入滲模型參數S和A的預測問題，即模型的輸出因子為入滲參數S和A。

S在數值上等于第一個單位時段末的累積入滲量，只與表層土壤的理化性質有關；而A主要取決于土壤結構狀態(tài)，即土壤孔隙的數量和連通性，受到整個濕潤土層理化參數的影響。相關研究[8]表明，影響入滲參數的主要因素為土壤體積含水率、干密度、質地和有機質含量，土壤體積含水率越高的土壤，其水吸力越小，土水勢越大，即水勢差越小，土壤達到飽和的速度越快，S值越小；土壤結構以干密度來表征，能夠反映土壤的孔隙率和密實度，土壤越密實，干密度越大，土壤孔隙越小，連通性越差，入滲通量也越小，入滲參數S和A隨之變??；土壤質地反映不同大小礦物顆粒的組成狀況，其中以黏粒含量作為劃分土壤類型的標準，土壤質地越重，黏粒含量越高，土壤孔隙度越小，孔隙越彎曲，連通性越差，對水分入滲的阻力越大，導致水力傳導度越小，并且固體顆粒表面積相應增大、吸附能力增強，S和A隨之變?。煌寥烙袡C質是由動、植物殘體分解合成的產物，可促進土壤團聚體的形成，增強土壤的穩(wěn)定性，有機質含量少的土壤，其結構不穩(wěn)定，孔隙不均勻，大孔隙相對較多，使得入滲初期的入滲通量較大，隨著水分的入滲，裂隙兩邊的土壤發(fā)生崩塌、擠壓，堵塞土壤中的大孔隙，入滲路徑嚴重彎曲，導致相對穩(wěn)定入滲率減小，因而，有機質含量少的土壤，吸濕率S值較大，相對穩(wěn)定入滲率A較小。

由上述分析可知，S主要受耕作層土壤物理性質的影響，A主要受耕作層和犁底層土壤理化性質的共同影響。土壤體積含水率以θ0(0～20 cm)、θ1(20～40 cm)作為指標；土壤結構以干密度γ0(0～10 cm)、γ1(10～20 cm)、γ2(20～40 cm)作為指標；考慮到土壤沙粒含量、粉粒含量和黏粒含量之間的線性關系，將耕作層(0～20 cm)土壤的粉粒含量w1(0.05～0.002 mm)、黏粒含量w2(<0.002 mm)以及犁底層(20～40 cm)土壤的粉粒含量w3、黏粒含量w4作為土壤質地指標；土壤有機質以耕作層土壤的有機質含量G作為指標。因此，S的回歸模型以土壤理化參數θ0、γ0、γ1、w1、w2、G作為輸入因子，A的回歸模型以土壤理化參數θ0、θ1、γ0、γ1、γ2、w1、w2、w3、w4、G作為輸入因子。試驗區(qū)部分土壤條件如表1所示。

表1 試驗區(qū)部分土壤條件Tab.1 Some soil condition of test area

2.2 模型的建立與分析

回歸分析是建立實測數據間的相互依賴關系，分析數據的內在規(guī)律，以用于預報、控制等問題，其中非線性回歸分析難度相對較大，為此，MATLAB語言得到廣泛應用。MATLAB軟件[13]是高級的數值分析、處理與計算軟件，具有優(yōu)秀的計算能力和卓越的數據化可視功能，在矩陣運算方面有其他語言難以比擬的優(yōu)越性，可節(jié)省大量的編程時間。土壤入滲參數與各理化參數之間的關系本質上是非線性的，由于入滲參數受到多個理化參數的共同影響，用多個理化參數的最優(yōu)組合共同來預測入滲參數相比只用一個理化參數更符合實際。因此，本研究以MATLAB軟件為計算工具，建立土壤理化參數與入滲參數之間的多元非線性傳遞函數，對試驗區(qū)100組樣本進行回歸，以更加精準的預測入滲參數，使預測模型更符合實際 。

(1)單因素函數形式的確定?；诖罅康膶崪y資料，在其他因素相同時，選取多組單因素試驗樣本，采用MATLAB中的Cftool工具，繪制散點圖和回歸曲線圖，根據曲線走勢，確定單因素最優(yōu)回歸方程，如表2所示。

由單因素最優(yōu)回歸方程的形式可知：在其他因子相同或接近相同時，入滲參數與體積含水率既存在對數關系，又存在線性關系，干密度與入滲參數呈線性關系，耕作層和犁底層土壤的粉粒、黏粒含量與入滲參數均呈對數關系，耕作層土壤的有機質含量與入滲參數為對數關系。

表2 入滲參數與各土壤理化參數間的單一定量關系Tab.2 Single quantitative relationship between each infiltration parameter and soil physicochemical parameters

(2)多元非線性函數的確定?；趯Ω饔绊懸蛩嘏c入滲參數內在關系的研究，確定單因素回歸方程，以各單因素回歸方程中的因子作為回歸模型的輸入因子，建立入滲參數的多元非線性回歸模型，并對各因子進行T檢驗，剔除每次檢驗中T值最小對應的因子，直至所有因子∣T∣≥T0.05/2為止(T0.05/2=2.014 4)，檢驗結果如表3所示。

表3 S回歸方程和A回歸方程自變量的確定Tab.3 Independent variables' determination of S and A regression equations

注：*表示每次進行T檢驗時的最小值，將其對應的因子剔除，不作為回歸方程的自變量。

對回歸方程輸入因子的檢驗結果可知，lnw2在S的回歸方程中不顯著，lnθ0、θ1、lnw3在A的回歸方程中不顯著。最終確定入滲參數回歸方程的自變量，如下所示：S回歸方程的自變量為lnθ0、θ0、γ0、γ1、lnw1、lnG；A回歸方程的自變量為θ0、lnθ1、γ0、γ1、γ2、lnw1、lnw2、lnw4、lnG。

基于土壤傳遞函數----非線性回歸分析法，以上述自變量為輸入參數，建立以S和A為輸出參數的多元非線性回歸方程，如下所示：

S=1.519 7-0.092 9 lnθ0+0.001 1θ0-0.635γ0+

0.040 8γ1+0.050 9 lnω1+0.140 7 lnG

(2)

A=0.321 5+0.000 2θ0-0.029 2 lnθ1-0.103 5γ0-

0.024 4γ1+0.036 4γ2-0.013 9 lnω1+0.015 4 lnω2-

0.012 lnω4+0.018 2 lnG

(3)

(3)灰色關聯(lián)分析。為反映自變量與預測變量的關聯(lián)程度，以Philip模型參數S和A為參考序列，各自變量為比較序列，依據灰色關聯(lián)理論[14,15]，對自變量進行排序，排序結果見表4。

表4 自變量的排序結果Tab.4 Sorted Results of independent variables

由表4可知，不同回歸模型中，自變量與預測值之間的關聯(lián)度不同，S回歸模型中自變量的關聯(lián)度順序為θ0>lnG>γ1>lnw1>γ0> lnθ0，A回歸模型中自變量的關聯(lián)度順序為lnw2>γ1> lnG>θ0>lnw1>γ0>γ2> lnθ1>lnw4。

(4)模型檢驗。在自變量檢驗的基礎上，對回歸方程進行顯著性、誤差以及擬合度分析，以保證方程整體的可靠性，分析結果如表5所示。

模型分析結果顯示：S和A回歸模型的F值均大于F0.05、平均相對誤差分別為7.53%和7.64%、負相關系數分別為0.882、0.933，累積入滲量I的綜合平均相對誤差為7.81%，則入滲參數S和A的多元非線性預測模型顯著性強、精度高、擬合度高。

表5 模型分析Tab.5 Model analysis

注：R為復相關系數，R越接近1，表明方程的擬合度越高。

2.3 預報實例

利用以下2組樣本，采用上述100組樣本建立的非線性預測模型對Philip模型參數進行預報，預報樣本的土壤條件如表6所示，預報結果如表7所示。

表6 預報樣本的土壤條件Tab.6 Soil condition of forecasting samples

表7 非線性預報結果Tab.7 Nonlinear forecasting results

由非線性預報結果可知，入滲參數S和A的預測相對誤差分別變化在5%～4.63%、5.46%～7.81%，累積入滲量I的綜合影響相對誤差為3.42%～4.32%，以上所得精度均可滿足非凍融條件下Philip模型參數的預報要求。預報結果表明用土壤體積含水率、干密度、粉粒含量、黏粒含量以及有機質含量對Philip模型參數進行非線性預報具有一定的可靠性。

3 結 語

(1)不同土壤條件下，用土壤理化參數對Philip入滲模型參數進行預報是可行的，即采用土壤傳遞函數----非線性回歸分析法，以土壤理化參數作為輸入變量，對Philip入滲模型參數S、A進行預報是可行的，其綜合相對誤差可控制在的8%以下，實現了對入滲參數的高度預測，解決了獲取入滲參數困難的問題。

(2)不同入滲參數的輸入變量不同，入滲參數S回歸模型的最終輸入變量為耕作層土壤(0～20 cm)的土壤體積含水率、干密度、地表以下0～10 cm土壤的干密度、粉粒含量和有機質含量；入滲參數A回歸模型的最終輸入變量為 耕作層土壤的土壤體積含水率、干密度、地表以下0～10cm土壤的干密度、粉粒含量、黏粒含量、有機質含量以及犁底層土壤的土壤體積含水率、干密度、黏粒含量。

本文所建立的Philip入滲模型參數的非線性預報模型，體現了非凍融期土壤理化參數與入滲參數之間的非線性本質，較好地反映了土壤水分入滲過程，克服了預測精度低、不能保證全局最優(yōu)等缺陷，但Philip模型僅適用均質一維垂直入滲的情況，還需進一步研究非均質土壤的入滲過程。

[1] 邵明安，王全九，黃明斌. 土壤物理學[M]. 北京：高等教育出版社，2006.

[2] 王全九，來劍斌，李 毅. Green-Ampt模型與Philip入滲模型的對比分析[J]. 農業(yè)工程學報， 2002,18(2):13-16.

[3] Green W H, Ampt G A.Studies on soil physics, flow of air and water through soils [J].J. Agr. Sci.，1911,76(4)：1-24.

[4] Horton RE．An approach toward a physical interpretation of infiltration-capacity[J]．Soil Sci., 1957,84(4):257-264．

[5] 邵明安，王全九. 推求土壤水分運動參數的簡單入滲法[J]. 土壤學報， 2000, 37(1):1-8.

[6] 武世亮，聶衛(wèi)波，馬孝義. 土壤特性與Philip入滲公式標定因子的空間變異性關系[J]. 排灌機械工程學報， 2014,32(8):730-736.

[7] 任尚崗，張振華，楊潤亞，等. Philip公式在三維入滲及參數測算中的應用[J]. 干旱地區(qū)農業(yè)研究， 2011,29(5):192-196.

[8] 原林虎. PHILIP入滲模型參數預報模型研究與應用[D]. 太原：太原理工大學， 2013.

[9] 雷志棟，楊詩秀，謝傳森.土壤水動力學[M].北京：清華大學出版社，1988．

[10] 朱安寧，張佳寶，陳效民，等．封丘地區(qū)土壤傳遞函數的研究[J]. 土壤學報， 2003,40(1):54-59.

[11] 馮錦萍，樊貴盛. 土壤入滲參數的線性傳輸函數研究[J]. 中國農村水利水電， 2014,(9):8-11.

[12] Philip J R．The theory of infiltration about sorptivity and algebraic infiltration equations[J]. Soi Sci., 1957，84(4)：257-264．

[13] 彭 勇，牛俊峰. 用MATLAB回歸非線性模型參數[J]. 化學通報， 2007,(11):880-883.

[14] 鄧聚龍.灰色系統(tǒng)控制系統(tǒng)[M].武漢：華中理工大學出版社， 1993.

[15] 劉思峰，黨耀國，方志耕.灰色系統(tǒng)理論及其應用[M].北京：科學出版社， 2004.