郭豪　劉文祥　王超然　黃智剛

摘要：以山東榮成地區(qū)棕壤、風(fēng)沙土為對(duì)象，測(cè)定分析土壤剖面不同層次的容重、質(zhì)地、有機(jī)質(zhì)含量、水勢(shì)和含水量等狀況;運(yùn)用Matlab軟件和神經(jīng)網(wǎng)絡(luò)法對(duì)土壤水分特征曲線van-Genuchten（簡(jiǎn)稱VG）模型進(jìn)行擬合，將2種擬合結(jié)果與實(shí)測(cè)數(shù)據(jù)進(jìn)行對(duì)比分析。結(jié)果表明，研究土壤不同土層含水量均隨著水吸力的增加呈“快速下降-緩慢下降-基本平穩(wěn)”的變化趨勢(shì);由于受到土壤顆粒組成、容重等物理因素和有機(jī)質(zhì)含量差異的影響，棕壤和風(fēng)沙土不同土層的飽和含水量和殘留含水量差異較大，其中棕壤60～70 cm和風(fēng)沙土10～20 cm的土壤殘留含水量最大;Matlab軟件擬合的殘差平方和比神經(jīng)網(wǎng)絡(luò)預(yù)測(cè)方法低幾個(gè)數(shù)量級(jí)，說(shuō)明采用Matlab對(duì)水分特征曲線VG模型擬合精度高，擬合效果較好。

關(guān)鍵詞：土壤水分特征曲線;VG模型;Matlab;神經(jīng)網(wǎng)絡(luò)預(yù)測(cè)

中圖分類號(hào)： S152.7? 文獻(xiàn)標(biāo)志碼： A? 文章編號(hào)：1002-1302（2019）03-0232-04

土壤水分特征曲線表示土壤水能量與土壤含水量之間的關(guān)系，是研究水分保持和運(yùn)動(dòng)的基本特征曲線[1]。土壤水分特征曲線測(cè)定方法主要有張力計(jì)法、離心機(jī)法和壓力膜法等直接測(cè)定法[2]。但由于影響土壤水分特征曲線的因素復(fù)雜，通過(guò)這些實(shí)測(cè)數(shù)據(jù)還難以從理論上獲得土壤水吸力與土壤含水量的關(guān)系式。因此，科研人員基于大量試驗(yàn)研究，提出了一些經(jīng)驗(yàn)公式來(lái)擬合實(shí)測(cè)數(shù)據(jù)[3]，其中，常見(jiàn)的有 Broods-Corey 模型[4]、Gardner模型[5]、Campbell模型[6]、van-Genuchten 模型[7]和Gardner-Russo模型[8]等。van-Genuchten （簡(jiǎn)稱VG）模型由于適用土壤質(zhì)地范圍廣、模擬線型與實(shí)測(cè)曲線吻合程度高而被廣泛應(yīng)用[9]。

朱蔚利等通過(guò)VG模型和土壤水分特征曲線單一參數(shù)模型擬合了土壤水分特征曲線，認(rèn)為VG模型擬合的精度比單一參數(shù)模型高一個(gè)數(shù)量級(jí)[10]。劉洪波等比較了RETC軟件和Matlab軟件對(duì)土壤水分特征曲線VG模型擬合效果，表明Matlab軟件擬合效果更好[11]。相關(guān)研究表明，采用Matlab對(duì)VG求參既快捷又準(zhǔn)確，是切實(shí)可行的方法[12]。

前人已經(jīng)做了大量關(guān)于土壤含水量與土壤水吸力擬合的研究，特別是不同經(jīng)驗(yàn)?zāi)Ｐ椭g的比較。針對(duì)同一經(jīng)驗(yàn)?zāi)Ｐ筒捎貌煌瑪M合方法的研究報(bào)道較少。本研究以棕壤和風(fēng)沙土為研究對(duì)象，采用離心機(jī)法測(cè)定土壤剖面不同深度土壤的水勢(shì)和對(duì)應(yīng)含水量，運(yùn)用Matlab軟件對(duì)土壤水分特征曲線VG模型參數(shù)進(jìn)行擬合計(jì)算，同時(shí)基于土壤物理性質(zhì)參數(shù)采用神經(jīng)網(wǎng)絡(luò)預(yù)測(cè)方法對(duì)VG模型求參，然后將2種擬合值與實(shí)測(cè)數(shù)據(jù)進(jìn)行對(duì)比，分析2種方法擬合效果差異，確定擬合研究區(qū)域土壤水分特征曲線VG模型的合適方法，以期為該區(qū)域的放射性核素在土壤中遷移數(shù)值模擬研究提供基礎(chǔ)數(shù)據(jù)。

1 材料與方法

1.1 試驗(yàn)材料

本研究選擇山東省榮成市石島灣核電廠周邊為研究區(qū)域，研究土壤主要是棕壤與風(fēng)沙土。土樣于2014年7月取自榮成市大疃李家村（地理位置：122°30′55″E、37°13′45″N）和后海崖村（地理位置：122°31′26″E、36°55′10″N），土壤類型分別對(duì)應(yīng)棕壤和風(fēng)沙土。每個(gè)樣地的取樣深度依次為0～10、10～20、…、60～70 cm，每種土壤類型取3個(gè)剖面。由于風(fēng)沙土土層較淺，60 cm以下是母質(zhì)層，取樣最大深度為50～60 cm，采集樣品的環(huán)刀是與高速冷凍離心機(jī)H-1400 pF配套的環(huán)刀，樣品采集后用離心機(jī)測(cè)量土壤含水量和土壤水吸力。

1.2 測(cè)定方法

用環(huán)刀采集原狀土樣，在實(shí)驗(yàn)室于蒸餾水中浸泡48 h，使其達(dá)到飽和狀態(tài)，取出稱質(zhì)量。用離心機(jī)對(duì)土樣進(jìn)行脫濕處理，共設(shè)置了8個(gè)不同轉(zhuǎn)速，分別為500、1 200、2 000、2 200、3 800、4 900、6 900、8 500 r/min，每次旋轉(zhuǎn)完成后都進(jìn)行稱質(zhì)量，根據(jù)稱質(zhì)量數(shù)據(jù)計(jì)算出不同轉(zhuǎn)速下對(duì)應(yīng)的土壤含水量。每個(gè)轉(zhuǎn)速對(duì)應(yīng)一個(gè)土壤水吸力，然后繪制土壤含水量與土壤水吸力關(guān)系曲線，即實(shí)測(cè)土壤水分特征曲線。土壤容重、有機(jī)質(zhì)含量等理化性質(zhì)測(cè)定根據(jù)土壤農(nóng)化分析方法[13]，土壤顆粒組成采用激光粒度儀（Microtrac S3500）測(cè)量。

1.3 分析方法

1.3.1 經(jīng)驗(yàn)?zāi)Ｐ?van-Genuchten模型：

式中：θ（h）為壓力水頭為h（kPa，或表示為基質(zhì)勢(shì)cm，1 cm=0.1 kPa）時(shí)的土壤體積含水量，cm3/cm3;α為尺度參數(shù)，kPa;n為曲線形態(tài)參數(shù)。

1.3.2 土壤水吸力估算 在離心機(jī)法測(cè)定土壤水分特征曲線過(guò)程中，土壤吸力可根據(jù)角速度和與之對(duì)應(yīng)的離心半徑確定。壓力水頭以及與其對(duì)應(yīng)的轉(zhuǎn)速和平衡時(shí)間見(jiàn)表1。

式中：h為土壤吸力（壓力水頭），cm;R1為離心機(jī)軸心到土樣中心的徑向距離，cm;R2為離心機(jī)軸心到土樣底部（即自由水面處）的徑向距離，cm;ρw為水密度，g/cm3;ω為角速度，r/s;g為重力加速度，cm/s2。

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

采用函數(shù)繪圖軟件Sigmaplot 12.0繪制土壤含水量和水勢(shì)關(guān)系曲線，采用Matlab 2012a和Hydrus 1D中神經(jīng)網(wǎng)絡(luò)預(yù)測(cè)對(duì)van-Genuchten模型進(jìn)行擬合，并求解相關(guān)參數(shù)值和殘差平方和。

2 結(jié)果與分析

2.1 2種土壤的水分特征曲線比較

由圖1、圖2可知，不同深度土壤水分特征曲線形態(tài)相似，隨著土壤水吸力的增大，土壤含水量先急劇降低，隨后緩慢減小，最終達(dá)到平穩(wěn)狀態(tài)。當(dāng)土壤水吸力為零時(shí)，土壤含水量達(dá)到飽和狀態(tài);當(dāng)土壤水吸力達(dá)到1 200 cm H2O后，棕壤不同深度土壤含水量出現(xiàn)差異;隨著土壤水吸力的增大，7個(gè)層次土壤含水量減小速率出現(xiàn)明顯差異，當(dāng)水吸力增大到一定程度時(shí)土壤含水量達(dá)到穩(wěn)定狀態(tài)，各層土壤含水量大小順序表現(xiàn)為40～50 cm>50～60 cm>20～30 cm>30～40 cm>10～20 cm>60～70 cm>0～10 cm，最終土壤殘留含水量介于0.09～0.15 cm3/cm3之間。

當(dāng)土壤水吸力在0～2 000 cm H2O時(shí)，風(fēng)沙土不同深度土壤含水量急劇下降（圖2），水吸力大于2 000 cm H2O后，土壤含水量隨水吸力的變化變緩。當(dāng)水吸力增大到一定程度時(shí)，土壤含水量達(dá)到穩(wěn)定狀態(tài)，各層土壤含水量大小順序表現(xiàn)為20～30 cm>10～20 cm>40～50 cm>30～40 cm>0～10 cm>50～60 cm。除20～30 cm和50～60 cm深度土壤含水量差異較大外，其他深度土壤含水量變化基本重合。

2.2 影響土壤水分特征曲線的因素

2.2.1 土壤理化性質(zhì) 由表2可知，除棕壤30～40 cm深度土壤外，2種土壤的容重均隨土壤深度的增加而增大，與表層土壤（0～10 cm）相比，2種土壤深層（50～60 cm）容重分別增加11.72%、11.56%。而土壤有機(jī)質(zhì)含量隨深度的增大呈減小趨勢(shì)，但是2種類型土壤有機(jī)質(zhì)含量均在20～30 cm深度時(shí)減少幅度最大（與10～20 cm土壤比較，均減少了約 5 g/kg）。2種類型土壤相比，相同深度土壤的容重和有機(jī)質(zhì)含量差異均較小。

2.2.2 土壤顆粒組成 從表3可以看出，棕壤和風(fēng)沙土顆粒組成以粉粒為主，其中棕壤不同深度的粉粒含量均高于80%，風(fēng)沙土粉粒含量也高達(dá)75%～84%。2種類型土壤中的黏粒含量均小于10%，表明研究區(qū)域黏粒含量偏低，土壤偏沙。比較2種土壤顆粒組成發(fā)現(xiàn)，棕壤沙粒含量總體比風(fēng)沙土砂粒低，而粉粒含量高于風(fēng)沙土，兩種土壤黏粒含量基本相同。根據(jù)美國(guó)土壤質(zhì)地分類標(biāo)準(zhǔn)[24]，2種土壤類型的質(zhì)地均是壤質(zhì)沙土。

2.3 VG模型擬合求參比較

基于實(shí)測(cè)的土壤含水量和水吸力數(shù)據(jù)，采用Matlab軟件中非線性擬合函數(shù)lsqcurvefit對(duì)VG模型參數(shù)求解。Matlab擬合程序?yàn)椋?/p>

在Matlab擬合VG模型求參的同時(shí)，采用神經(jīng)網(wǎng)絡(luò)預(yù)測(cè)方法對(duì)VG進(jìn)行擬合求參。神經(jīng)網(wǎng)絡(luò)預(yù)測(cè)是Hydrus 1D中Rosetta Lite窗口的方法，擬合時(shí)選用輸入方式是土壤顆粒組成和容重，然后根據(jù)實(shí)測(cè)土壤顆粒組成和容重?cái)?shù)據(jù)求解VG模型中的參數(shù)。通過(guò)Matlab和神經(jīng)網(wǎng)絡(luò)預(yù)測(cè)擬合值與實(shí)測(cè)值的殘差平方和比較，確定選擇合適的擬合方法。

從表4可以看出，2種方法求解VG模型參數(shù)結(jié)果差異較大。土壤殘留含水量（θr）和飽和含水量（θs）擬合結(jié)果顯示，預(yù)測(cè)值隨土壤深度變化規(guī)律與實(shí)測(cè)值相同，不同深度的θr和θs不同。Matlab軟件對(duì)2種土壤θr和θs預(yù)測(cè)值均大于神經(jīng)網(wǎng)絡(luò)預(yù)測(cè)值。比較2種擬合方法的預(yù)測(cè)值與實(shí)測(cè)值之間的殘差平方和可知，在相同深度土層，采用神經(jīng)網(wǎng)絡(luò)預(yù)測(cè)的殘差平方和比Matlab擬合的大2～3個(gè)數(shù)量級(jí)，表明Matlab對(duì)VG模型擬合效果更好。

3 討論與結(jié)論

研究發(fā)現(xiàn)，2種土壤不同土層土壤含水量均隨著水吸力的增加呈“快速下降-緩慢下降-基本平穩(wěn)”的變化趨勢(shì)。但在同一水吸力下，各土層所保持土壤水分含量有所不同。由于受到土壤顆粒組成、容重等物理因素和有機(jī)質(zhì)含量差異的影響，棕壤和風(fēng)沙土不同土層的飽和含水量和殘留含水量差異較大，其中棕壤60～70 cm和風(fēng)沙土10～20 cm的土壤殘留含水量最大。

土壤水動(dòng)力學(xué)參數(shù)特征受土壤結(jié)構(gòu)、土壤質(zhì)地及土壤有機(jī)質(zhì)含量等因素的影響，土壤容重的增大導(dǎo)致土壤含水量減少，持水能力下降[14];土壤黏粒含量增多，細(xì)小孔隙數(shù)量亦增多，土壤毛管持水作用增強(qiáng)，土壤顆粒的比表面積增大，對(duì)水分吸附力增強(qiáng)，從而提高了土壤持水能力[15];土壤有機(jī)質(zhì)對(duì)自然及風(fēng)干狀態(tài)下的土壤含水量、土壤水分特征曲線均有顯著影響[16]。本研究表明，2種研究土壤的容重均隨土壤深度增加而增大，且棕壤在土層深度達(dá)到30～40 cm以后、風(fēng)沙土在土層達(dá)到20～30 cm以后土壤容重大小基本趨于穩(wěn)定，土壤有機(jī)質(zhì)含量隨深度增加呈減小的趨勢(shì)，而土壤的粒徑組成則隨著土層深度變化呈現(xiàn)不規(guī)律變化。影響土壤質(zhì)地的決定因素是成土母質(zhì)類型[1]，產(chǎn)生這些變化的原因可能是受到農(nóng)田管理措施特別是耕作的影響。表層土壤經(jīng)翻耕后，容重變小，孔隙增多，持水能力增強(qiáng)，這與傳統(tǒng)耕作和免耕土壤水分特征曲線特性相同[17]。此外，作物類型特別是根系的機(jī)械作用，也會(huì)對(duì)土壤持水特性產(chǎn)生不同程度的影響[18]。本研究采樣區(qū)域位于耕作區(qū)，主要種植作物是冬小麥和夏玉米，農(nóng)田的耕作和作物輪作導(dǎo)致不同土層的土壤殘留含水量和飽和含水量發(fā)生改變。

與2種土壤含水量的實(shí)測(cè)值相比，神經(jīng)網(wǎng)絡(luò)方法預(yù)測(cè)的殘差平方和較大，這是因?yàn)楸狙芯克痔卣髑€是基于土壤質(zhì)地和容重估算的，與考慮土壤飽和含水量等土壤其他理化因素的人工神經(jīng)網(wǎng)絡(luò)ANN預(yù)測(cè)效果[19]相比，本研究中神經(jīng)網(wǎng)絡(luò)預(yù)測(cè)法在對(duì)VG模型擬合求參時(shí)精度相對(duì)于Matlab較低。因此，在比較Matlab和神經(jīng)網(wǎng)絡(luò)法對(duì)土壤水分特征曲線VG模型擬合求參時(shí)，當(dāng)土壤含水量和水吸力易于測(cè)量時(shí)，采用Matlab軟件擬合VG模型精度高且簡(jiǎn)單高效，反之，當(dāng)考慮土壤飽和含水量等其他理化因素時(shí)，則可以采用神經(jīng)網(wǎng)絡(luò)預(yù)測(cè)方法獲取VG模型參數(shù)，但在采用神經(jīng)網(wǎng)絡(luò)預(yù)測(cè)時(shí)需要考慮影響土壤結(jié)構(gòu)的諸多因素。

參考文獻(xiàn)：

[1]呂貽忠，李保國(guó). 土壤學(xué)[M]. 北京：中國(guó)農(nóng)業(yè)出版社，2006.

[2]van Genuchten M T，Leij F J，Lund L J. Indirect methods for estimating the hydraulic properties of unsaturated soils[J]. Proceedings of the International Workshop on Indirect Methods for Estimating the Hydraulic Properties of Unsaturated Soils，1989，157（6）：615-627.

[3]劉建立，徐紹輝，劉 慧. 估計(jì)土壤水分特征曲線的間接方法研究進(jìn)展[J]. 水利學(xué)報(bào)，2004，35（2）：68-76.

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

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

[6]Campell G S. A simple method for determining unsaturated conductivity from moisture retention data[J]. Soil Science，1974，117（6）：311-314.

[7]van Genuchten M T. A closed-form equation for predicting the hydraulic conductivity of unsaturated soils[J]. Soil Science Society of America Journal，1980，44（5）：892-898.

[8]Russo D. Determining soil hydraulic properties by parameter estimation：on the selection of a model for the hydraulic properties[J]. Water Resources Research，1988，24（3）：453-459.

[9]付 強(qiáng)，蔣睿奇，王子龍，等. 基于改進(jìn)螢火蟲(chóng)算法的土壤水分特征曲線參數(shù)優(yōu)化[J]. 農(nóng)業(yè)工程學(xué)報(bào)，2015，31（11）：117-122.

[10]朱蔚利，肖自幸，牛健植，等. 兩種模型對(duì)土壤水分特征曲線擬合的比較分析[J]. 湖南農(nóng)業(yè)科學(xué)，2011（17）：47-51.

[11]劉洪波，張江輝，虎膽·吐馬爾白，等. 土壤水分特征曲線VG模型參數(shù)求解對(duì)比分析[J]. 新疆農(nóng)業(yè)大學(xué)學(xué)報(bào)，2011，34（5）：437-441.

[12]魏義長(zhǎng)，劉作新，康玲玲，等. 土壤持水曲線van Genuchten模型求參的Matlab實(shí)現(xiàn)[J]. 土壤學(xué)報(bào)，2004，41（3）：380-386.

[13]魯如坤. 土壤農(nóng)業(yè)化學(xué)分析[M]. 北京：中國(guó)農(nóng)業(yè)出版社，1999.

[14]李 卓，吳普特，馮 浩，等. 容重對(duì)土壤水分蓄持能力影響模擬試驗(yàn)研究[J]. 土壤學(xué)報(bào)，2010，47（4）：611-620.

[15]楊欣坤，王 宇，趙蘭坡，等. 土壤水動(dòng)力學(xué)參數(shù)及其影響因素研究進(jìn)展[J]. 中國(guó)農(nóng)學(xué)通報(bào)，2014，30（3）：38-43.

[16]馮 杰，郝振純，劉方貴. 大孔隙對(duì)土壤水分特征曲線的影響[J]. 灌溉排水，2002，21（3）：4-7.

[17]王小華，賈克力，劉景輝，等. Van Genuchten模型在土壤水分特征曲線擬合分析中的應(yīng)用[J]. 干旱地區(qū)農(nóng)業(yè)研究，2009，27（2）：179-183，188.

[18]馬昌臣，王 飛，穆興民，等. 小麥根系機(jī)械作用對(duì)土壤水分特征曲線的影響[J]. 水土保持學(xué)報(bào)，2013，27（2）：105-109.

[19]胡振琪，張學(xué)禮. 基于ANN的復(fù)墾土壤水分特征曲線的預(yù)測(cè)研究[J]. 農(nóng)業(yè)工程學(xué)報(bào)，2008，24（10）：15-19.榮 榮，鄭育聲，楊林生，等. 生物炭對(duì)雞糞堆肥過(guò)程中氨氣排放的影響[J]. 江蘇農(nóng)業(yè)科學(xué)，2019，47（3）：236-240.