鹽堿土壤Philip入滲模型參數(shù)的非線性預報模型

沈 婧，樊貴盛

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

0 引 言

據(jù)糧農(nóng)組織不完全統(tǒng)計，全世界鹽堿地面積為9.544 億hm2，其中我國約為0.991 億hm2，是世界第三大鹽堿地國家。我國鹽堿地廣泛分布在華北、東北、西北及濱海地區(qū)，其中鹽堿化耕地和鹽堿荒地總面積達0.333 億hm2，約占耕地總面積的27.5%，作為一種重要的可開發(fā)利用的土地資源，采取合理措施改良鹽堿土以提高土地利用率和生產(chǎn)率對緩解我國緊張的人地矛盾意義重大。為此，國內(nèi)眾多學者在鹽堿地改良方面進行了諸多探索[1-4]，包括灌水洗鹽、平整深耕、選種育苗、間作方式、施用改良劑、田間灌排等，取得了顯著成果。

鹽堿地的形成受到氣候條件、地理條件、土壤條件和地下水等多種因素的綜合影響，但其形成的實質(zhì)是各種可溶性鹽分的隨水運移，地表集聚。目前，國內(nèi)外鹽堿地的改良多采用灌排工程技術(shù)，其原理是利用灌溉水入滲實現(xiàn)對土壤鹽分的有效淋洗。鹽堿地高含鹽量尤其是高鈉離子含量使土體結(jié)構(gòu)性差，水力傳導度低，水分入滲緩慢，一方面造成了灌溉水利用率低，洗鹽效果不佳，同時也加劇了干旱半干旱地區(qū)的水資源短缺。因此，合理確定鹽堿地水分入滲參數(shù)對提高灌水利用率，實現(xiàn)節(jié)水型鹽堿地改良意義重大。呂殿青等[5]在室內(nèi)進行了一維積水入滲試驗，分析了土壤鹽分特征量與水分入滲量之間的關(guān)系，最終得出兩者的函數(shù)關(guān)系。張莉等[6]進行了夾砂層土壤入滲過程水鹽運動規(guī)律的室內(nèi)模擬試驗，觀測了在灌水條件下夾砂層和均質(zhì)土壤的水分入滲過程并對比分析了夾砂層土壤的鹽分運動規(guī)律，該研究為濱海鹽堿地改良和開發(fā)利用提供了理論依據(jù)。PATEL等[7]通過對印度古吉拉特邦北部的鹽堿土壤進行的水分入滲試驗發(fā)現(xiàn)，入滲率受鈉離子吸附、土壤黏粒含量的影響，建議要依據(jù)鹽堿地的類型選擇合理的種植作物。任長江等[8]通過建立水分累計入滲量、入滲率與初始含鹽量的多元回歸方程，分析了土壤初始含鹽量對入滲特性的影響，結(jié)果發(fā)現(xiàn)：相同時間初始含鹽量越高，累計入滲量和入滲率越小。縱觀現(xiàn)有研究成果，在鹽堿地水鹽運移規(guī)律及土壤水分入滲機理方面研究較多，成果顯著，但試驗法獲取鹽堿地土壤水分入滲參數(shù)工作量大，耗時費力，因此，本文以Philip[9]半理論半經(jīng)驗入滲模型參數(shù)為研究對象，采用土壤傳輸函數(shù)的方法，對晉北鹽堿地區(qū)的土壤水分入滲參數(shù)進行非線性預測，為間接獲取鹽堿地土壤水分入滲參數(shù)和鹽堿地改良提供了一種新思路。

1 材料與方法

1.1 試驗區(qū)概況及土壤條件

土壤水分原位入滲試驗在山西省北部朔州市應(yīng)縣鹽堿地進行，地勢東南高，西北低，屬北溫帶大陸性季風氣候，四季分明，年際溫差大，據(jù)應(yīng)縣氣象局多年觀測資料，日溫差可達17.1 ℃，年溫差高達74.1 ℃，年均氣溫7 ℃左右，降雨分布極不均勻，主要集中在6-9月，年均降水量為380 mm，年均蒸發(fā)量為1 918.6 mm(20 cm蒸發(fā)皿)，年均蒸發(fā)量遠大于年均降雨量。

該縣屬海河流域和桑干河水系，洪積傾斜平原地貌特征明顯，地下水資源豐富，地下水位最高為1.5 m，礦化度落在0.5～1.5 g/L之間。試驗區(qū)土壤堿化程度相對較輕而鹽漬化程度相對較重，pH值介于7.0～8.3之間，電導率范圍為0.123 7～1.916 6 S/m2(0～20 cm)和0.125 0～3.158 1 S/m2(20～40 cm)，試驗區(qū)土壤理化參數(shù)情況見表1。

1.2 試驗設(shè)備與方法

鹽堿地土壤原位入滲試驗采用雙套環(huán)入滲儀法，雙套環(huán)為2個同心圓環(huán)(內(nèi)、外環(huán)直徑分別為26 、60 cm，兩環(huán)高度均為25 cm)，內(nèi)環(huán)控制試驗的入滲面積，外環(huán)起保證內(nèi)環(huán)水體垂直下滲、減少側(cè)向滲流的作用。試驗開始前，將雙套環(huán)打入地表以下20 cm的深度。用自制的水位控制器對內(nèi)外環(huán)水位進行控制，且水位基本保持在地表以上2 cm處。試驗過程中，準確記錄內(nèi)環(huán)加水的時間和注水量，以便日后擬合入滲參數(shù)。試驗研究表明，土壤水分入滲在60 min左右基本達到穩(wěn)定入滲，本試驗選定入滲時間為90 min，以確保大田土壤水分入滲都達到穩(wěn)定入滲狀態(tài)。

在進行野外試驗的同時，還需在室內(nèi)對土壤理化指標進行測定，包括土壤全鹽量、土壤容重、體積含水率、土壤有機質(zhì)、土壤質(zhì)地。土壤全鹽量采用電導法測定；土壤容重通過蠟封法進行測定；體積含水率測定用傳統(tǒng)烘干稱重法獲得；試驗點土壤有機質(zhì)量的測定是利用化學方法通過重鉻酸鉀滴定法得到；土壤質(zhì)地通過比重計法得到篩分曲線，然后分析土壤的顆粒級配，進而確定土壤質(zhì)地；氣溫等其他信息通過試驗站氣象設(shè)施觀察得到。

表1 試驗區(qū)土壤理化參數(shù)范圍Tab.1 Soil Physical and chemical parameters in the test area

1.3 Philip入滲模型

1957年，在系統(tǒng)研究Richard方程的基礎(chǔ)上，Philip提出了能較好地描述土壤水分一維垂直入滲情況的入滲模型，該模型是半經(jīng)驗半理論性的模型，具有明確的物理意義，入滲公式中只含有兩個入滲參數(shù)，形式簡單，是被人們認可的、表述土壤水分入滲過程的模型之一，在世界范圍內(nèi)得到較廣泛的應(yīng)用[10]。Philip入滲公式如下：

I=St0.5+At,

(1)

式中：I為累積入滲量：自入滲開始至t時刻的入滲總量，cm；S為吸滲率：入滲開始后第一個單位時段末的累積入滲量與穩(wěn)定入滲率之差，在數(shù)值上等于入滲開始至一分鐘末的累積入滲量減去穩(wěn)滲率，cm/min0.5；t為入滲歷時，min；A為穩(wěn)滲率，在單位勢梯度下飽和土壤的入滲速度或非飽和土壤入滲達到相對穩(wěn)定階段的入滲速度[11]，cm/min。

1.4 建模樣本

本文建立了200組由Philip入滲模型參數(shù)和鹽堿土壤基本理化參數(shù)構(gòu)成的樣本，其中197組作為建模訓練樣本，3組用于檢驗所建預報模型的精度。

表2 部分建模樣本數(shù)據(jù)Tab.2Partial sample data

2 輸入因子和輸出因子的確定

2.1 鹽堿地土壤理化特性對入滲參數(shù)的影響分析

(1)全鹽量對入滲參數(shù)的影響：鹽堿地土壤的可溶性鹽濃度若高于作物的耐鹽度，則會影響作物的正常生長，輕則作物出苗率低，重則作物生長受到抑制甚至不生長。大量試驗資料表明：在土壤結(jié)構(gòu)、質(zhì)地、含水率、有機質(zhì)等條件相近的情況下，鹽堿地土壤的入滲能力明顯小于非鹽堿地土壤的入滲能力，根本原因在于：與非鹽堿地的水力傳導度相比，鹽堿地土壤的水力傳導度變小[13]，全鹽量越高，土體團聚性越差，土粒越分散，結(jié)構(gòu)性越差[14]，且在入滲過程中，土粒隨著水流逐漸下移沉積的可能性越大[15]，可能導致大孔隙孔徑減小、小孔隙孔徑更小甚至被堵塞，水流下滲受到的阻力增大，吸滲率S越小；在單位勢梯度下，單位土壤面積通過的流量減小，水力傳導度減小，入滲量小，穩(wěn)滲率A越小。

(2)有機質(zhì)對入滲參數(shù)的影響：有機質(zhì)具有很好的膠結(jié)作用，比砂粒的黏結(jié)力大，促進砂土團粒狀結(jié)構(gòu)的形成，比黏粒的黏結(jié)力小，降低黏粒間的黏結(jié)力，使土壤的大、小孔隙分配合理，從而改善了土壤的通透性，增強孔隙的穩(wěn)定性，改良土壤結(jié)構(gòu)。在其他條件一定的條件下，有機質(zhì)量低的土壤，較大、較小孔隙居多，水分進入土壤初期，較大孔隙先充滿水，并在重力作用下向下層土壤入滲，S值較大；反之，有機質(zhì)量高的土壤，中等孔隙發(fā)育較好，在水分入滲初期，S值較小。穩(wěn)滲率A的大小主要由孔隙的大小和多少決定，有機質(zhì)量高的土壤，中等孔隙居多且土壤團粒結(jié)構(gòu)穩(wěn)定性較好，A值較大；有機質(zhì)量低的土壤，大孔隙和小孔隙居多且土壤團粒結(jié)構(gòu)穩(wěn)定性差，水在入滲過程中，土粒容易隨著水流移動，使小孔隙孔徑減小，甚者堵塞小孔隙，迫使入滲水流改變方向，下滲路徑曲折蜿蜒，A值較小。

(3)土壤結(jié)構(gòu)對入滲參數(shù)的影響：土壤結(jié)構(gòu)是指土壤的松緊程度、孔隙程度和板結(jié)程度，用土壤干容重(土壤相對體積質(zhì)量)作為反映土壤結(jié)構(gòu)的指標[16]。若土壤干密度大，土壤越密實，土壤孔隙就小，孔隙連通性差，水分入滲路徑越彎曲，水流受到的阻力越大，單位勢梯度下，水分通量小，土壤的水力傳導度就小，則S和A的值越小；反之，越大。

(4)體積含水率對入滲參數(shù)的影響：土壤含水率決定著土壤固、液相的比例，常用體積含水率表示。在土壤質(zhì)地、結(jié)構(gòu)相同或相近的情況下，在入滲開始后的很短時間內(nèi)，表層土壤很快達到飽和，若土壤含水量大，勢梯度就小，單位過水斷面的入滲通量就小，S值就小；在土壤結(jié)構(gòu)一定的條件下，含水率越大，孔隙相對減少，穩(wěn)滲率A就越小。

(5)土壤質(zhì)地對入滲參數(shù)的影響：土壤質(zhì)地的分類和劃分標準有很多種，本文采用國際制對鹽堿土壤質(zhì)地進行劃分，并以黏粒含量、粉粒含量和砂粒含量作為數(shù)學表征值。在土壤結(jié)構(gòu)、含水率相同的條件下，質(zhì)地越重，黏粒含量越多，固體相比表面積越大，顆粒表面的吸附能力越大，粒間孔隙孔徑越小，在單位勢梯度下，入滲水流路徑嚴重曲折，導致水分通量較小，水力傳導度低，入滲參數(shù)S和A都減小。

2.2 預報模型輸出和輸入因子的確定

預報模型的輸出因子為Philip入滲模型的入滲參數(shù),即吸滲率S和穩(wěn)滲率A。

預報模型的輸入因子根據(jù)2.1分析的影響因素初步確定，然后由輸入變量的顯著性檢驗決定取舍。在輸入因子的初步確定中，考慮了以下三點：

(1)入滲參數(shù)受土壤層次的影響。從吸滲率S的概念可知，入滲開始較短時間內(nèi)，水分還未能到達犁底層，所以影響S的因素僅限于耕作層；而影響A的因素既有耕作層的也有犁底層的。

(2)黏粒、粉粒和砂粒的數(shù)學關(guān)系：由于黏粒含量(%)+粉粒含量(%)+砂粒含量(%)=100%，且三者對入滲參數(shù)的影響高度相關(guān)，故本文只選擇黏粒含量和粉粒含量作為入滲參數(shù)的輸入因子[17]。

(3)鹽分的集聚性：鄒長明等[18]研究認為：鹽分主要積聚在0～15 cm土層，尤其是0～5 cm土層。馬獻發(fā)等[19]更是提出：0～10 cm土層的鹽分含量遠遠高于其下各土層。因而，只選擇耕作層的全鹽量作為輸入因子。

那么，初步確定吸滲率S的輸入因子為：μ1、μ2、G、ω1、ω2、γ1、γ2、θ1；初步確定穩(wěn)滲率A的輸入因子為：μ1、μ2、G、ω1、ω2、ω*1、ω*2、γ1、γ2、γ3、θ1、θ2。

3 入滲參數(shù)預報模型的建立與檢驗

3.1 建立模型

建立模型的目的是找到一個數(shù)學表達式，以測定的土壤常規(guī)理化參數(shù)為自變量，計算出不易測定的入滲參數(shù)。本文基于鹽堿地197組田間入滲試驗樣本，借助MATLAB[12]軟件，建立鹽堿地Philip入滲模型參數(shù)與其土壤理化參數(shù)間的關(guān)系模型。其具體操作步驟和結(jié)果如下：

(1)確定入滲參數(shù)與單個輸入因子的函數(shù)關(guān)系。運用控制變量法，選擇一個因子作為唯一變量，要求其他因子相同或相近，對數(shù)據(jù)進行篩選，選出5～10組滿足要求的數(shù)據(jù)，借助MATLAB軟件的Cftool功能[12]，輸入單個因子及入滲參數(shù)的實測值，由離散點的走勢，選擇或定義函數(shù)形式，參考擬合度的大小，最終確定函數(shù)的最佳形式。表3即入滲參數(shù)與各個輸入因子的函數(shù)關(guān)系。

由表3可知：全鹽量μ與入滲參數(shù)的關(guān)系既包含對數(shù)關(guān)系又包含線性關(guān)系；有機質(zhì)G與入滲參數(shù)呈指數(shù)關(guān)系；容重γ與入滲參數(shù)呈線性關(guān)系；含水率θ與入滲參數(shù)的關(guān)系既包含對數(shù)關(guān)系又包含線性關(guān)系；黏、粉粒含量ω與入滲參數(shù)呈線性關(guān)系。

表3 入滲參數(shù)與單個輸入因子的函數(shù)形式Tab.3 The function form of infiltration parametersand single input factors

(2)確定回歸方程的初步形式。將表2給出的各函數(shù)進行機械相加，并添加一個常數(shù)項φ，便得出回歸方程的初步形式：

S=α1lnμ1+α2μ2+α3eG+α4γ1+α5γ2+

α6lnθ1+β6θ1+α7ω1+α8ω2+φ

(2)

A=α1lnμ1+β1μ1+α2lnμ2+β2μ2+α3eG+

α4γ1+α5γ2+α6γ3+α7lnθ1+α8lnθ2+α9ω1+

α10ω2+α11ω*1+α12ω*2+φ

(3)

(3)回歸方程中自變量的顯著性檢驗(T檢驗)。將樣本數(shù)據(jù)和預報模型的初步函數(shù)形式輸入到已編好的程序中，運行后可得到初步模型各項輸入變量前的回歸系數(shù)及其進行T檢驗后的結(jié)果，若存在某個|T|≤T0.05/2的值(顯著水平α=0.05，T0.05/2=2.014 1)，則剔除掉最小T值所對應(yīng)的自變量，繼續(xù)對剩余自變量進行T檢驗，直至所有自變量對應(yīng)的|T|值都大于T0.05/2后，停止檢驗。表4是所有自變量進行T檢驗的過程和結(jié)果。

表4 T檢驗過程及結(jié)果Tab.4 Inspection process and results

注：“___”表示T檢驗中的最小值。

由表4知，經(jīng)過4次T檢驗后，lnμ1、μ2、eG、γ1、lnθ1和ω1對應(yīng)的|T|都大于T0.05/2，即在入滲參數(shù)S的回歸方程中顯著，那么S最終的輸入因子：μ1、μ2、G、γ1、ω1、θ1；經(jīng)過5次檢驗后，lnμ1、μ1、lnμ2、eG、γ1、γ2、lnθ1、ω2、lnθ2和ω*1對應(yīng)的|T|都大于T0.05/2，即在入滲參數(shù)A的回歸方程中顯著，那么A最終的輸入因子：μ1、μ2、G、γ1、γ2、θ1、ω2、γ3、θ2、ω*1。

(4)回歸方程的顯著性檢驗(F檢驗)。除了要對方程自變量進行顯著性檢驗，還要對整個回歸方程自身進行顯著性檢驗，即F檢驗，以確保回歸方程顯著。給定顯著水平α=0.05，由樣本情況查表可知Fα=0.05值，若計算得到的F值大于Fα=0.05，則方程顯著；反之，不顯著。F檢驗結(jié)果見表5。

表5 F檢驗結(jié)果Tab.5 F test results

由表5可知：FS和FA都大于Fα=0.05,即入滲參數(shù)S和A的回歸方程顯著，兩者的平均誤差分別為8.062%、9.043%。則入滲參數(shù)的回歸方程的最終形式為：

S=0.033 6 lnμ1-0.000 03μ2+0.049eG-

0.637 2γ1-0.038 5 lnθ1-0.009 9ω1+1.265 7

(4)

A=-0.040 6 lnμ1+0.000 03μ1-0.010 4 lnμ2+0.025eG-

0.037 9γ1-0.021 4γ2+0.009 1 lnθ1-0.002 3 lnθ2-

0.000 6ω2-0.002 4ω*1+0.400 3

(5)

3.2 預報模型檢驗

3.2.1 單參數(shù)預測結(jié)果檢驗

表6是3組檢驗樣本(未參與建模的樣本)的土壤條件，采用以上所建模型對檢驗樣本的入滲參數(shù)進行預報，預報結(jié)果見表7。

表6 檢驗樣本的土壤條件Tab.6 Soil conditions for test samples

表7 預報結(jié)果Tab.7 Forecast result

由表7可知，預測吸滲率S的相對誤差范圍是6.977%～8.358%；預測穩(wěn)滲率A的相對誤差范圍是6.25%～11.688%，計算I90與實測I90的相對誤差介于1.969%～11.331%之間，預報結(jié)果的精度在可接受范圍。

3.2.2 Philip入滲模型二參數(shù)綜合檢驗

將197組建模樣本的預測值S和A帶入到Phillip入滲公式中，得到對應(yīng)的90 min累積入滲量I90的計算值，與197組建模樣本的實測值I90相比較，兩者的擬合結(jié)果見圖1。

圖1 實測值I90與計算值I90擬合效果Fig.1 Fitting effect ofmeasured values and calculated value

由圖1所示：197組實測I90與計算I90的擬合程度較高(R2=0.985)，說明由入滲參數(shù)S、A計算得到的I90比較接近實測I90，進一步說明所建模型的預報精度較高，模型較可靠。

綜上所述，無論是預報模型的單參數(shù)預測結(jié)果檢驗，還是Philip入滲模型二參數(shù)綜合結(jié)果檢驗，都表明預報精度較高。因此，所建立的多元非線性預報模型是可行的。

4 結(jié) 語

(1)以常規(guī)土壤理化參數(shù)為輸入因子的多元非線性模型對鹽堿地Philip入滲模型參數(shù)A和S進行預測，兩參數(shù)平均誤差都小于10%，在建模誤差范圍之內(nèi)，預報模型的參數(shù)檢驗誤差也在可接受范圍之內(nèi)，故所建模型是可行的，可實現(xiàn)對鹽堿地入滲參數(shù)的準確預報。研究結(jié)果也為鹽堿地土壤灌水技術(shù)參數(shù)的合理確定和鹽堿地改良提供了理論與技術(shù)支撐。

(2)本文選擇對入滲參數(shù)影響較大的土壤理化特性作為模型的輸入因子，其他如地表溫度、地下水礦化度以及土壤酸堿度等次要因素均未考慮在內(nèi)，因此，在后續(xù)的研究中應(yīng)豐富影響因素，以建立精度更高的預報模型。

□

[1] 李慶梅，侯龍魚，劉 艷，等. 黃河三角洲鹽堿地不同利用方式土壤理化性質(zhì)[J].中國生態(tài)農(nóng)業(yè)學報，2009，17(6)：1 132-1 136.

[2] 張小娟，韓 明，鄭敏娜. 不同油用向日葵品種在鹽堿地的適應(yīng)性研究[J].山西農(nóng)業(yè)科學，2012，40(4)：332- 335.

[3] QingjieWanga, CaiyunLua, Hongwen Lia,et al. The effects of no-tillage with subsoiling on soil properties and maize yield: 12-Year experiment on alkaline soils of Northeast China[J]. Soil and Tillage Research，2014，137：43-49.

[4] Chun Xia Wang, Quan Jiu Wang, Zhuang Liang. Field Experimental Study on Saline and Alkaline Soil Amelioration through Combining Chemical Method and Under-Mulch Drip Irrigation. Advanced Materials Research, 2012，599：832-837.

[5] 呂殿青，王全九，王文焰，等. 土壤鹽分分布特征評價[J].土壤學報，2002，39(5)：720-725.

[6] 張 莉，丁國棟，王翔宇，等. 夾砂層土壤入滲過程水鹽運動規(guī)律室內(nèi)模擬試驗研究[J]. 安徽農(nóng)業(yè)科學，2010，38(7)：3 605-3 606.

[7] B B PATEL, BHARAT B PATEL, R S DAVE. Studies on infiltration of saline-alkali soils of several parts of mehsana and patan districts of north gujarat[J]. Journal of Applied Technology in Environmental Sanitation.2011，(1)：87-92.

[8] 任長江，白 丹，周 文，等. 土壤初始含鹽量對水分入滲特性的影響[J].干旱區(qū)研究，2014，31(2)：222-225.

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

[10] 原林虎. PHILIP入滲模型參數(shù)預報模型研究與應(yīng)用[D]. 太原：太原理工大學，2013：6-8.

[11] 岳海晶. 樊貴盛.備耕頭水地土壤入滲參數(shù)的線性預報模型[J]. 中國農(nóng)村水利水電，2016,(2)：23.

[12] 汪嵐，黃彩虹. 基于MATLAB色差預測多元回歸模型的研究[J]. 計算機與應(yīng)用化學，2008，25(8)：1 015-1 018．

[13] 樊貴盛，李 堯，蘇冬陽，等. 大田原生鹽堿荒地入滲特性的試驗[J].農(nóng)業(yè)工程學報. 2012，28(19)：63-70.

[14] So H B,，AYLMORE L A G. How do sodic soils behave? The effects of sadicity on soil physical behavior [J]. Aust. J. Soil Res. 1993，31.

[15] 王 雪，樊貴盛. Na+含量對土壤入滲能力影響的試驗研究[J].太原理工大學學報，2009，40(4)：391-393.

[16] 樊貴盛，李雪轉(zhuǎn)，李紅星. 非飽和土壤介質(zhì)水分入滲問題的試驗研究[M]. 北京：中國水利水電出版社，2012.

[17] 郭 華. 土壤入滲模型參數(shù)的分階段非線性預報模型研究[D].太原：太原理工大學，2016：46.

[18] 鄒長明，張多姝，張曉紅，等．蚌埠地區(qū)設(shè)施土壤酸化與鹽漬化狀況測定與評價[J]．安徽農(nóng)學通報，2006，12(9):54-55．

[19] 馬獻發(fā)，白路平，張繼舟.哈爾濱市郊設(shè)施土壤積鹽規(guī)律的研究[J]．黑龍江農(nóng)業(yè)科學，2009，11(1)：50-52．