侯 雯，寇飛周，周穎明(陜西華地勘察設(shè)計(jì)咨詢有限公司，西安 710018)

地下滴灌是把灌水毛管及其灌水器埋入土壤中，將水或水肥的混合液緩慢滲入到作物根區(qū)土壤，再借助毛細(xì)管作用或重力作用將水分?jǐn)U散到根系層供作物吸收利用的一種灌水方法[1,2]。

王曉愚[3]等以灌水器工作壓力、土壤初始含水率和土壤容重為影響因素，在試驗(yàn)的基礎(chǔ)上，建立了地下滴灌灌水器流量計(jì)算經(jīng)驗(yàn)公式，但公式經(jīng)驗(yàn)參數(shù)是針對(duì)某一種土壤試驗(yàn)獲得的。其實(shí)對(duì)不同質(zhì)地的土壤，地下滴灌灌水器的流量大小也不盡相同[4]。Lazarovitch et al.[5]基于正壓地下滴灌灌水器穩(wěn)定流量計(jì)算公式[6]和毛管的沿程水頭損失計(jì)算公式相結(jié)合，同時(shí)考慮了土壤的空間變異性，提出了地下滴灌毛管水力計(jì)算方法，這一方法首次考慮了土壤空間變異性對(duì)地下滴灌毛管水力要素的影響，但其中反映土壤空間變異性的有關(guān)參數(shù)是隨機(jī)生成的，因此這一方法缺乏試驗(yàn)驗(yàn)證。Gil et al.[7]在試驗(yàn)的基礎(chǔ)上，綜合考慮了滴頭制造偏差等因素，提出了地下滴灌毛管水力計(jì)算方法，這一方法具有一定的試驗(yàn)基礎(chǔ)，但試驗(yàn)裝置與工程實(shí)際相差太大。

滴頭埋于土中，地下滴灌毛管水力要素的空間分布必然受到土壤物理要素的影響，因此土壤的空間變異性會(huì)影響地下滴灌毛管水力要素的空間分布規(guī)律。本文以此為研究目的，結(jié)合試驗(yàn)建立空間變異性影響下的地下滴灌毛管水力要素計(jì)算模型。為揭示土壤空間變異性和地下滴灌毛管水力特性的內(nèi)在聯(lián)系，地下滴灌系統(tǒng)的規(guī)劃、設(shè)計(jì)、管理和產(chǎn)品開發(fā)提供理論依據(jù)。

1 地下滴灌條件下土壤物理參數(shù)的空間變異性

1.1 試驗(yàn)方案

試驗(yàn)區(qū)位于北京市平谷區(qū)東高村鎮(zhèn)的東高村，試驗(yàn)布置如圖1所示。

圖1 試驗(yàn)布置示意圖Fig.1 Experimental arrangement diagram

試驗(yàn)區(qū)地下水埋深為83 m，采用機(jī)井作為灌溉水源。試驗(yàn)地塊為200 m×35 m大小。支管和毛管的埋深均為0.4 m。支管選用φ60PE管，毛管(滴灌管)選用φ16PE管，支管兩側(cè)各布置8條毛管，毛管長度為100 m，毛管上滴頭間距為0.3 m，毛管之間的間距為5 m。壓力表用來測(cè)量毛管進(jìn)口處壓力值，水表用來測(cè)量毛管進(jìn)口處水量，根據(jù)灌水時(shí)間，可計(jì)算得出毛管入口處流量值。土壤初始含水率和土壤黏粒含量測(cè)定沿著毛管鋪設(shè)方向從右往左進(jìn)行。沿著毛管方向每隔5 m設(shè)置一個(gè)取樣點(diǎn)，在取樣點(diǎn)土壤表層下深度為60 cm處取土樣，采用烘干法[8]測(cè)定土壤初始含水率；沿著毛管方向每隔5 m設(shè)置一個(gè)取樣點(diǎn)，在每個(gè)取樣點(diǎn)土壤表層下深度30、60、90 cm取土樣，采用激光粒度儀[9]測(cè)定土壤顆粒級(jí)配。對(duì)該試驗(yàn)地塊的土壤初始含水率和土壤顆粒級(jí)配的取樣值進(jìn)行統(tǒng)計(jì)特征分析和空間分布特征分析，探求整片試驗(yàn)地塊上土壤初始含水率和土壤黏粒含量的空間分布格局，為研究土壤物理參數(shù)(土壤初始含水率和土壤黏粒含量)空間變異性對(duì)地下滴灌毛管水力特性的影響奠定基礎(chǔ)。

1.2 土壤物理參數(shù)的統(tǒng)計(jì)特征分析

采用SPSS19.0軟件包K-S檢驗(yàn)[10]或z的標(biāo)準(zhǔn)正態(tài)Q-Q圖[11]來分析土壤初始含水率和土壤黏粒含量在土壤中分布的正態(tài)性。根據(jù)SPSS軟件對(duì)樣本均值及標(biāo)準(zhǔn)差的分析來確定各因素的變異性質(zhì)。

1.2.1土壤初始含水率的統(tǒng)計(jì)特征分析

按經(jīng)典統(tǒng)計(jì)方法對(duì)土壤初始含水率的統(tǒng)計(jì)特征值進(jìn)行分析：剔除異常點(diǎn)后的樣本數(shù)為308個(gè)，土壤初始含水率的變化范圍為12.43%～23.74%，其平均值為17.34%。土壤初始含水率的變異系數(shù)為11.31%，屬于中等變異[12]的范圍。土壤初始含水率用K-S正態(tài)分布檢驗(yàn)概率(PK-S)[13]進(jìn)行檢驗(yàn)，取顯著水平α=0.05，土壤初始含水率PK-S=0.200>0.05，符合標(biāo)準(zhǔn)正態(tài)分布，滿足地質(zhì)統(tǒng)計(jì)學(xué)方法分析的前提。

1.2.2土壤黏粒含量的統(tǒng)計(jì)特征分析

選用0.01 mm以下粒徑的顆粒作為研究對(duì)象[14]。對(duì)不同取樣點(diǎn)處各深度土樣的黏粒含量測(cè)定結(jié)果取平均值進(jìn)行分析。土壤黏粒含量的樣本數(shù)為327個(gè)，變化范圍為19.33%～60.71%，平均值為39.25%。土壤黏粒含量的變異系數(shù)為21.27%，屬于中等變異的范圍。

根據(jù)K-S檢驗(yàn)，土壤黏粒含量不符合標(biāo)準(zhǔn)正態(tài)分布。對(duì)其進(jìn)行對(duì)數(shù)轉(zhuǎn)換，經(jīng)對(duì)數(shù)轉(zhuǎn)換后，消除了可能存在的比例效應(yīng)。經(jīng)對(duì)數(shù)轉(zhuǎn)換后的數(shù)據(jù)進(jìn)行K-S值檢驗(yàn)，土壤黏粒含量對(duì)數(shù)轉(zhuǎn)換后的PK-S=0，仍不滿足標(biāo)準(zhǔn)正態(tài)分布。采用標(biāo)準(zhǔn)正態(tài)Q-Q圖來檢驗(yàn)數(shù)據(jù)的正態(tài)性，若樣本數(shù)據(jù)近似于正態(tài)分布，在標(biāo)準(zhǔn)Q-Q圖上這些點(diǎn)近似成一條直線，由圖2所示可知，土壤黏粒含量取對(duì)數(shù)后滿足近似正態(tài)分布。

圖2 土壤黏粒含量對(duì)數(shù)轉(zhuǎn)換后標(biāo)準(zhǔn)Q-Q圖Fig.2 Standard Q-Q plot of soil clay content after natural logarithmic transformation

1.3 土壤物理參數(shù)的空間分布特征分析

在GS+ 9.0平臺(tái)上，土壤初始含水率、土壤黏粒含量分別用線性、球狀、高斯和指數(shù)等理論模型進(jìn)行擬合，根據(jù)擬合誤差大小選出最適模型，得出土壤物理參數(shù)的半方差函數(shù)。用普通克里格法(Ordinary Kriging)內(nèi)插獲得空間變異分布圖。通過Kriging插值圖進(jìn)行土壤物理參數(shù)空間變異特征的分析。

1.3.1土壤初始含水率的空間分布特征分析

經(jīng)檢驗(yàn)，土壤初始含水率服從標(biāo)準(zhǔn)正態(tài)分布，可直接進(jìn)行半方差函數(shù)模型的擬合。其變異函數(shù)理論模型和擬合參數(shù)如表1所示。

表1 土壤初始含水率變異函數(shù)理論模型擬合參數(shù)Tab.1 Semivariogram theoretical models and contents and fitting parameters of the soil initial moisture content

土壤初始含水率的塊金系數(shù)為0.510，說明其具有中等的空間相關(guān)性[15]。

1.3.2土壤黏粒含量的空間分布特征分析

經(jīng)檢驗(yàn)，土壤黏粒含量經(jīng)過對(duì)數(shù)轉(zhuǎn)換后滿足近似正態(tài)分布，因此可進(jìn)行半方差函數(shù)模型的擬合，其變異函數(shù)理論模型擬合參數(shù)見表2。

土壤黏粒含量經(jīng)過對(duì)數(shù)轉(zhuǎn)換后的半方差函數(shù)理論模型的塊金系數(shù)為0.640，具有中等的空間相關(guān)性。

表2 土壤黏粒含量變異函數(shù)理論模型擬合參數(shù)Tab.2 Semivariogram theoretical models and contents and fitting parameters of the soil clay content

圖3、圖4分別為土壤初始含水率和對(duì)數(shù)轉(zhuǎn)換后的土壤黏粒含量的空間分布插值圖。

2 地下滴灌毛管水力要素研究

2.1 地下滴灌毛管上單滴頭公式的建立

考慮土壤初始含水率和土壤黏粒含量的空間變異性，建立單滴頭公式如下：

圖3 土壤初始含水率(%)的Kriging插值圖Fig.3 Kriging interpolation map of the soil initial moisture content(%)

圖4 對(duì)數(shù)轉(zhuǎn)換后的土壤黏粒含量(%)的Kriging插值圖Fig.4 Kriging interpolation map of the soil clay content after natural logarithmic transformation(%)

qi=kHxiθaiMbi

(1)

式中：qi為地下滴灌毛管上第i個(gè)滴頭流量，L/h；Hi為地下滴灌毛管上第i個(gè)滴頭壓力，m；θi為地下滴灌毛管上第i個(gè)滴頭處土壤初始含水率(重量)，%；Mi為地下滴灌毛管上第i個(gè)滴頭處土壤黏粒含量，%；k為經(jīng)驗(yàn)系數(shù)；a、b和x為經(jīng)驗(yàn)指數(shù)。

2.2 地下滴灌毛管水力要素計(jì)算模型的建立

地下滴灌毛管沿程壓力水頭分布如圖5所示。毛管上有n個(gè)滴頭，0點(diǎn)為毛管入口，在毛管入口處安裝水表和壓力表，可測(cè)得毛管入口流量Q0和入口端壓力H0。

圖5 毛管上沿程壓力水頭分布示意圖Fig.5 Pressure head distribution of along multiple-outlet pipeline

將毛管從入口0處開始按滴頭位置分為n個(gè)管段，由能量方程可知對(duì)于第i段管段有：

平臺(tái)功能：1）檢索，提供標(biāo)題、關(guān)鍵詞、長詞等檢索，通過語義分析給予最符合的檢索結(jié)果，并在詳情頁推薦相似文獻(xiàn)；II）訂閱服務(wù)，讀者可根據(jù)作者、主題、關(guān)鍵字或出版標(biāo)準(zhǔn)來選擇出版物提醒服務(wù)。

Hi-1+dZ-ΔHi-Hi=0

(2)

由圖5和式(2)可列出毛管上各管段能量方程組成的非線性方程組：

(3)

dz=IL

(4)

式中：Hi為毛管上第i個(gè)滴頭處的壓力水頭，m；dZ為相鄰滴頭的高差，m；ΔHi為第i段管段的水頭損失，m；I為地面坡度，I<0為上坡，I=0為平坡，I>0為下坡；L為相鄰滴頭的間距，m。

管道水頭損失由沿程水頭損失和局部水頭損失組成。沿程水頭損失可按下式進(jìn)行確定：

(5)

式中：hf為管道沿程水頭損失，m；f為摩阻系數(shù)，毛管為聚乙烯管，管徑16 mm，查規(guī)范[16]，f=0.505；Q為管道流量，L/h；D為管道內(nèi)徑，mm；L為管段長度(本處即指滴頭間距)，m；m為流量系數(shù)，毛管為聚乙烯管，管徑16 mm，查規(guī)范[16]，m=1.75；B為管徑指數(shù)，毛管為聚乙烯管，管徑16 mm，查規(guī)范[16]，B=4.75。

管道的局部水頭損失可按沿程水頭損失的一定比例估算。因此，管段水頭損失的計(jì)算可按下式確定：

ΔH=αhf

(6)

式中：α為毛管局部水頭損失加大系數(shù)。

將式(5)代入式(6)，即可得各管道水頭損失計(jì)算公式：

(7)

(8)

式中：Qi為毛管第i段管段流量，L/h；qi為毛管上第i個(gè)滴頭的流量，L/h。

由式(1)可得：

(9)

將式(7)～(9)和(4)代入式(3)中，即可建立在考慮土壤空間變異性的條件下毛管水力要素計(jì)算的數(shù)學(xué)模型。

(10)

毛管坡度I、毛管管徑D、滴頭個(gè)數(shù)n、滴頭間距L已知，毛管上各滴頭處土壤初始含水率和土壤黏粒含量的數(shù)據(jù)可通過Kriging插值得出，毛管進(jìn)口處壓力H0由試驗(yàn)測(cè)得。式(10)是一個(gè)由333個(gè)方程組成的非線性方程組，可求解333個(gè)滴頭流量qi，通過式(9)也可以求解毛管上的沿程壓力值。

2.3 地下滴灌毛管水力要素計(jì)算模型參數(shù)求解

毛管水力要素計(jì)算數(shù)學(xué)模型中，除滴頭流量qi外，a、b、x、k、α也為未知量。為求解模型中的相關(guān)未知參數(shù)，可將其轉(zhuǎn)化為目標(biāo)函數(shù)的優(yōu)化問題。以Matlab遺傳算法工具箱[17]為主要工具，解決有約束的最優(yōu)化問題。

根據(jù)試驗(yàn)方案，毛管的進(jìn)口流量值已知，因此可建立目標(biāo)函數(shù)為：毛管進(jìn)口流量的均方差最小。即：

(12)

式中：j為用來求解模型中參數(shù)的毛管條數(shù)；Q0計(jì)算為模型計(jì)算得出的毛管進(jìn)口流量值，即模型計(jì)算后所得毛管上所有滴頭流量之和；Q0為試驗(yàn)測(cè)得的毛管進(jìn)口流量值；n為每條毛管上的滴頭總個(gè)數(shù)；qi為滴頭流量值，由毛管水力要素計(jì)算數(shù)學(xué)模型求解得出。

目標(biāo)函數(shù)的求解與參數(shù)a、b、x、k、α的取值有關(guān)。假設(shè)已知a、b、x、k、α的取值范圍，則目標(biāo)函數(shù)是尋找一組最優(yōu)的(a，b，x，k，α)的組合，使得毛管進(jìn)口流量的均方差達(dá)到最小。

根據(jù)相關(guān)文獻(xiàn)[18]，可用a、b、x、k和α這5個(gè)參量的上下限來約束目標(biāo)函數(shù)的值，令目標(biāo)函數(shù)達(dá)到最小。a、b、x、k和α這5個(gè)參量的上下限取值如表3所示。

表3 優(yōu)化變量上下限取值Tab.3 The upper and lower bounds of optimizing variables

對(duì)14組毛管試驗(yàn)數(shù)據(jù)通過遺傳算法進(jìn)行迭代，可求出a、b、x、k、α這5個(gè)參量的最優(yōu)解如表4所示。

表4 優(yōu)化參量最優(yōu)解Tab.4 optimal solution of optimizing variables

為驗(yàn)證模型的可靠性和正確性，將優(yōu)化目標(biāo)函數(shù)的方法求得的參量a、b、k、x、α代入模型計(jì)算得出毛管入口流量Q0計(jì)算，通過比較試驗(yàn)測(cè)量值Q0和模型計(jì)算值Q0計(jì)算的大小，利用另外2組毛管的試驗(yàn)數(shù)據(jù)進(jìn)行校核和驗(yàn)證。對(duì)比結(jié)果如表5所示。

表5 地下滴灌毛管水力要素試驗(yàn)結(jié)果與水力計(jì)算模擬結(jié)果Tab.5 Experiment and simulation results of hydraulic calculation of the lateral under SDI

通過比較可知，入口流量試驗(yàn)值Q0與入口流量計(jì)算值Q0計(jì)算之間的誤差較小，說明根據(jù)試驗(yàn)建立的地下滴灌毛管水力要素計(jì)算模型式(10)是正確的。因此可確定考慮土壤物理參數(shù)(土壤初始含水率和土壤黏粒含量)空間變異性條件下的地下滴灌滴頭流量計(jì)算公式：

qi=0.520 1H0.583 7θ-0.040 4M-0.077 8

(13)

3 結(jié) 語

(1)經(jīng)統(tǒng)計(jì)特征分析，土壤初始含水率和土壤黏粒含量均屬于中等變異的范圍；由半方差函數(shù)分析可知，土壤初始含水率和土壤黏粒含量均符合球狀模型，且均具有中等相關(guān)性的空間分布特征。通過Kriging插值可以直觀地反映土壤初始含水率和土壤黏粒含量的空間分布情況。

(2)考慮土壤物理參數(shù)(土壤初始含水率和土壤黏粒含量)的空間變異性，建立單滴頭公式。通過分析地下滴灌毛管沿程壓力水頭和流量分布情況，根據(jù)地下滴灌毛管管段能量方程建立考慮土壤物理參數(shù)空間變異性的地下滴灌毛管水力要素計(jì)算模型。

(3)采用MATLAB遺傳算法工具箱迭代求解地下滴灌毛管水力要素計(jì)算模型中的未知參數(shù)a、b、x、k、α的最優(yōu)解，進(jìn)而運(yùn)用MATLAB非線性方程組數(shù)值解法編程求解毛管上的滴頭流量，從而得出毛管入口流量。比較毛管入口流量的實(shí)測(cè)值和計(jì)算值，由于相對(duì)誤差較小，實(shí)測(cè)值和計(jì)算值較吻合，說明建立的地下滴灌毛管水力要素計(jì)算模型較為合理。

□

[1] 寇 丹,蘇德榮,吳 迪，等.地下調(diào)虧滴灌對(duì)紫花苜蓿耗水、產(chǎn)量和品質(zhì)的影響[J].農(nóng)業(yè)工程學(xué)報(bào),2014,30(2):116-123.

[2] Patel N, Rajput T B S. Dynamics and modeling of soil water under subsurface drip irrigated onion[J]. Agric. Water Manage., 2008,95(12):1 335-1 349.

[3] 王曉愚,白 丹,李占斌,等.土壤物理特性對(duì)地下滴灌灌水器流量影響分析[J].干旱區(qū)資源與環(huán)境,2009,23(3):126-129.

[4] Gil M, Rodríguez-Sinobas L, Sánchez R, et al. Procedures for determining maximum emitter discharge in subsurface drip irrigation[J]. J. Irrig. Drain Eng., 2011,137:287-294.

[5] Lazarovitch N, Shani U, Thompson T L, et al. Soil hydraulic properties affecting discharge uniformity of gravity-fed subsurface drip irrigation systems [J]. ASCE, Journal of Irrigation and Drainage Engineering, 2006,132(6):531-536.

[6] Warrick A W., Shani U. Soil-limiting flow from subsurface emitters. II: effect on uniformity[J].J. Irrig. Drain. Eng. ASCE,1996,122(5):296-299.

[7] Gil M, Rodriguez-Sinobas L, Juana L, et al. Emitter discharge variability of subsurface drip irrigation in uniform soils: Effect on water-application uniformity [J]. Irrigation Science,2008,26(6):451-458.

[8] 陳立新. 土壤實(shí)驗(yàn)實(shí)習(xí)教程[M].哈爾濱：東北林業(yè)大學(xué)出版社,2005.

[9] 程 鵬,高 抒,李徐生.激光粒度儀測(cè)試結(jié)果及其與沉降法、篩析法的比較[J].沉積學(xué)報(bào),19(3):449-455.

[10] Burgess T M, Webster R. Optimal interpolation and isarithmic mapping of soil properties. I. The semivariogram and punctual Kriging[J]. Soil Sci., 1980,(31):315-331.

[11] 陳翠英,江永真,袁朝春.土壤特性空間變異性研究[J].農(nóng)業(yè)機(jī)械學(xué)報(bào),2005,36(10):121-124.

[12] 洪 楠.Statistic for Windows 統(tǒng)計(jì)與圖表分析教程[M].北京:清華大學(xué)出版社,2002:152-159.

[13] 王政權(quán).地質(zhì)統(tǒng)計(jì)學(xué)及在生態(tài)學(xué)中的應(yīng)用[M]. 北京：科學(xué)出版社,1999.

[14] 龔振平.土壤學(xué)與農(nóng)作學(xué)[M]. 北京：中國水利水電出版社,2008.

[15] 古麗娜爾·托爾提,海米提·依米提,米日姑·買買提，等.伊犁河谷土壤含鹽量空間變異和格局分析[J].干旱地區(qū)農(nóng)業(yè)研究,2011,29(2):152-158.

[16] GB/T 50485-2009，微灌工程技術(shù)規(guī)范[S].

[17] 吳祈宗.運(yùn)籌學(xué)與最優(yōu)MATLAB編程[M]. 北京：機(jī)械工業(yè)出版社,2009:195-205.

[18] 李 剛,王曉愚,白丹.地下滴灌中毛管水力計(jì)算的數(shù)學(xué)模型與試驗(yàn)[J].排灌機(jī)械工程學(xué)報(bào),2011,29(1):87-92.