一種新的地下滴灌毛管水力計算模型
——基于土壤初始含水率和土壤黏粒含量空間變異性研究

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

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

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

滴頭埋于土中，地下滴灌毛管水力要素的空間分布必然受到土壤物理要素的影響，因此土壤的空間變異性會影響地下滴灌毛管水力要素的空間分布規律。本文以此為研究目的，結合試驗建立空間變異性影響下的地下滴灌毛管水力要素計算模型。為揭示土壤空間變異性和地下滴灌毛管水力特性的內在聯系，地下滴灌系統的規劃、設計、管理和產品開發提供理論依據。

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

1.1 試驗方案

試驗區位于北京市平谷區東高村鎮的東高村，試驗布置如圖1所示。

圖1 試驗布置示意圖Fig.1 Experimental arrangement diagram

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

1.2 土壤物理參數的統計特征分析

采用SPSS19.0軟件包K-S檢驗[10]或z的標準正態Q-Q圖[11]來分析土壤初始含水率和土壤黏粒含量在土壤中分布的正態性。根據SPSS軟件對樣本均值及標準差的分析來確定各因素的變異性質。

1.2.1土壤初始含水率的統計特征分析

按經典統計方法對土壤初始含水率的統計特征值進行分析：剔除異常點后的樣本數為308個，土壤初始含水率的變化范圍為12.43%～23.74%，其平均值為17.34%。土壤初始含水率的變異系數為11.31%，屬于中等變異[12]的范圍。土壤初始含水率用K-S正態分布檢驗概率(PK-S)[13]進行檢驗，取顯著水平α=0.05，土壤初始含水率PK-S=0.200>0.05，符合標準正態分布，滿足地質統計學方法分析的前提。

1.2.2土壤黏粒含量的統計特征分析

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

根據K-S檢驗，土壤黏粒含量不符合標準正態分布。對其進行對數轉換，經對數轉換后，消除了可能存在的比例效應。經對數轉換后的數據進行K-S值檢驗，土壤黏粒含量對數轉換后的PK-S=0，仍不滿足標準正態分布。采用標準正態Q-Q圖來檢驗數據的正態性，若樣本數據近似于正態分布，在標準Q-Q圖上這些點近似成一條直線，由圖2所示可知，土壤黏粒含量取對數后滿足近似正態分布。

圖2 土壤黏粒含量對數轉換后標準Q-Q圖Fig.2 Standard Q-Q plot of soil clay content after natural logarithmic transformation

1.3 土壤物理參數的空間分布特征分析

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

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

經檢驗，土壤初始含水率服從標準正態分布，可直接進行半方差函數模型的擬合。其變異函數理論模型和擬合參數如表1所示。

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

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

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

經檢驗，土壤黏粒含量經過對數轉換后滿足近似正態分布，因此可進行半方差函數模型的擬合，其變異函數理論模型擬合參數見表2。

土壤黏粒含量經過對數轉換后的半方差函數理論模型的塊金系數為0.640，具有中等的空間相關性。

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

圖3、圖4分別為土壤初始含水率和對數轉換后的土壤黏粒含量的空間分布插值圖。

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

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

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

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

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

qi=kHxiθaiMbi

(1)

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

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

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

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

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

平臺功能：1）檢索，提供標題、關鍵詞、長詞等檢索，通過語義分析給予最符合的檢索結果，并在詳情頁推薦相似文獻；II）訂閱服務，讀者可根據作者、主題、關鍵字或出版標準來選擇出版物提醒服務。

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

(2)

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

(3)

dz=IL

(4)

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

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

(5)

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

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

ΔH=αhf

(6)

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

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

(7)

(8)

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

由式(1)可得：

(9)

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

(10)

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

2.3 地下滴灌毛管水力要素計算模型參數求解

毛管水力要素計算數學模型中，除滴頭流量qi外，a、b、x、k、α也為未知量。為求解模型中的相關未知參數，可將其轉化為目標函數的優化問題。以Matlab遺傳算法工具箱[17]為主要工具，解決有約束的最優化問題。

根據試驗方案，毛管的進口流量值已知，因此可建立目標函數為：毛管進口流量的均方差最小。即：

(12)

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

目標函數的求解與參數a、b、x、k、α的取值有關。假設已知a、b、x、k、α的取值范圍，則目標函數是尋找一組最優的(a，b，x，k，α)的組合，使得毛管進口流量的均方差達到最小。

根據相關文獻[18]，可用a、b、x、k和α這5個參量的上下限來約束目標函數的值，令目標函數達到最小。a、b、x、k和α這5個參量的上下限取值如表3所示。

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

對14組毛管試驗數據通過遺傳算法進行迭代，可求出a、b、x、k、α這5個參量的最優解如表4所示。

表4 優化參量最優解Tab.4 optimal solution of optimizing variables

為驗證模型的可靠性和正確性，將優化目標函數的方法求得的參量a、b、k、x、α代入模型計算得出毛管入口流量Q0計算，通過比較試驗測量值Q0和模型計算值Q0計算的大小，利用另外2組毛管的試驗數據進行校核和驗證。對比結果如表5所示。

表5 地下滴灌毛管水力要素試驗結果與水力計算模擬結果Tab.5 Experiment and simulation results of hydraulic calculation of the lateral under SDI

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

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

(13)

3 結 語

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

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

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

□

[1] 寇 丹,蘇德榮,吳 迪，等.地下調虧滴灌對紫花苜蓿耗水、產量和品質的影響[J].農業工程學報,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] 王曉愚,白 丹,李占斌,等.土壤物理特性對地下滴灌灌水器流量影響分析[J].干旱區資源與環境,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] 陳立新. 土壤實驗實習教程[M].哈爾濱：東北林業大學出版社,2005.

[9] 程 鵬,高 抒,李徐生.激光粒度儀測試結果及其與沉降法、篩析法的比較[J].沉積學報,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].農業機械學報,2005,36(10):121-124.

[12] 洪 楠.Statistic for Windows 統計與圖表分析教程[M].北京:清華大學出版社,2002:152-159.

[13] 王政權.地質統計學及在生態學中的應用[M]. 北京：科學出版社,1999.

[14] 龔振平.土壤學與農作學[M]. 北京：中國水利水電出版社,2008.

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

[16] GB/T 50485-2009，微灌工程技術規范[S].

[17] 吳祈宗.運籌學與最優MATLAB編程[M]. 北京：機械工業出版社,2009:195-205.

[18] 李 剛,王曉愚,白丹.地下滴灌中毛管水力計算的數學模型與試驗[J].排灌機械工程學報,2011,29(1):87-92.