基于DEM 與Logistic 函數(shù)的灌區(qū)渠系工作制度模擬
——以周橋灌區(qū)為例

劉 永，繳錫云,2,3*，程明瀚，李 江

（1.河海大學(xué) 農(nóng)業(yè)工程學(xué)院，南京 210098；2.水文水資源與水利工程國(guó)家重點(diǎn)實(shí)驗(yàn)室， 南京210098；3.河海大學(xué) 水安全與水科學(xué)協(xié)同創(chuàng)新中心，南京 210098）

0 引 言

【研究意義】農(nóng)業(yè)生產(chǎn)是我國(guó)的用水主體，灌溉是農(nóng)業(yè)生產(chǎn)的前提。根據(jù)《第一次全國(guó)水利普查公報(bào)》，全國(guó)大中型灌區(qū)的灌溉面積占總灌溉面積的50.20%[1]，因此在水資源緊缺的大背景下，大中型灌區(qū)的水資源高效利用具有重要的意義?！狙芯窟M(jìn)展】目前，對(duì)于大中型灌區(qū)來(lái)說(shuō)，渠道級(jí)別多、數(shù)量大，渠道輸水效率、干支渠水位穩(wěn)定性、末級(jí)渠道水位高低等對(duì)于灌區(qū)的高效、穩(wěn)定運(yùn)行具有重要的影響[2-4]，因此渠系調(diào)度是大中型灌區(qū)高效運(yùn)行與節(jié)水的重要環(huán)節(jié)。渠系工作制度優(yōu)化的目標(biāo)主要分為2 大類(lèi)[5]：一是以灌區(qū)某次配水增產(chǎn)效益最大為目的優(yōu)化渠道配水。前人基于人工魚(yú)群算法[6]、“Generalized Wolfe”算法[7]、作物水分生產(chǎn)函數(shù)[8]等模型對(duì)灌區(qū)渠系配水進(jìn)行模擬優(yōu)化，實(shí)現(xiàn)對(duì)灌區(qū)灌溉效益的優(yōu)化。二是以輸水損失量最小為目標(biāo)對(duì)渠道進(jìn)行優(yōu)化配水；宋靜茹等[9]、馬孝義等[10]基于遺傳算法，以輪灌期間所有配水時(shí)段的渠道輸水損失最小為目標(biāo)建立了渠道優(yōu)化配水模型，解決了渠道引水流量不穩(wěn)定情況下的渠系優(yōu)化配水決策問(wèn)題。Suryavanshi 等[11]建立0-1 線性規(guī)劃配水模型，對(duì)渠道配水及輪灌方式進(jìn)行優(yōu)化，減少了灌區(qū)渠道輸水損失及工程成本。宋松柏等[12]、呂宏興等[13]提出渠道配水均一化方法，應(yīng)用遺傳算法對(duì)配水渠道流量?jī)?yōu)化求解，實(shí)現(xiàn)同級(jí)渠道的同時(shí)關(guān)閉，解決了配水渠道閘門(mén)多次調(diào)節(jié)的問(wèn)題。

【切入點(diǎn)】對(duì)于平原灌區(qū)，自流灌溉面積也是衡量灌區(qū)渠系工作制度優(yōu)劣的一個(gè)重要指標(biāo)，不同渠系工作制度會(huì)產(chǎn)生不同的渠道水位，進(jìn)而對(duì)灌區(qū)自流灌溉面積產(chǎn)生巨大影響。然而，上述渠系工作制度優(yōu)化研究主要是以灌區(qū)效益最大化、輸水損失小為目標(biāo)函數(shù)，沒(méi)有將自流灌溉面積作為優(yōu)化目標(biāo)考慮進(jìn)去?！緮M解決的關(guān)鍵問(wèn)題】周橋灌區(qū)為洪澤湖東岸的大型平原自流灌區(qū)，灌溉配水制約因素較多。灌區(qū)自流灌溉面積受水位影響較大，往往由于10～30 cm 的水頭差，導(dǎo)致整個(gè)灌區(qū)中大面積農(nóng)田無(wú)法自流灌溉。因此，本文以周橋灌區(qū)為例，以ASTER GDEM 數(shù)據(jù)[14-15]及實(shí)測(cè)高程為基礎(chǔ)，結(jié)合Logistic 函數(shù)建立渠首閘后水位與灌區(qū)自流灌溉面積的關(guān)系模型，模擬不同渠系工作制度下周橋灌區(qū)的自流灌溉情況，確定最優(yōu)渠系工作制度，以提高灌區(qū)自流灌溉面積，為周橋灌區(qū)灌溉決策提供參考。

1 材料與方法

1.1 研究區(qū)概況

周橋灌區(qū)位于江蘇省淮安市洪澤區(qū)中北部，東經(jīng)118°47′40″—119°07′35″，北緯33°11′01″—33°23′25″，西靠洪澤湖，南以草澤河為界、與淮安洪金灌區(qū)毗鄰，東連白馬湖，北連蘇北灌溉總渠?？绺吡紳尽ⅫS集、朱壩、岔河4 鎮(zhèn)（街道）全部及東雙溝鎮(zhèn)部分區(qū)域，耕地面積2.433 萬(wàn)hm2，約占全區(qū)耕地面積的62%，設(shè)計(jì)灌溉面積2.133 萬(wàn)hm2。周橋灌區(qū)主要通過(guò)渠首周橋洞引用洪澤湖水及黃集洞、薛橋洞從蘇北灌溉總渠補(bǔ)水進(jìn)行灌溉。灌區(qū)現(xiàn)有總干渠、硯臨干渠、潯南干渠、潯北干渠等4 條干渠，總長(zhǎng)48.2 km；流量1 m3/s以上的渠道、灌排結(jié)合渠道96 條，總長(zhǎng)302 km。周橋灌區(qū)渠系分布見(jiàn)圖1。

圖1 周橋灌區(qū)渠系分布圖 Fig.1 Distribution of channel system in Zhouqiao irrigation district

1.2 研究方法

1.2.1 數(shù)字高程模型（DEM）構(gòu)建

茲通過(guò)實(shí)測(cè)高程數(shù)據(jù)對(duì)ASTER GDEM 數(shù)據(jù)進(jìn)行擬合，構(gòu)建周橋灌區(qū)DEM。ASTER GDEM 數(shù)據(jù)為NASA 于2009 年發(fā)布的覆蓋全球98%以上陸地DEM，是目前使用廣泛的數(shù)字高程模型，具有覆蓋度廣、精度高的特點(diǎn)[16-17]，DEM 是模擬地表高程變化特征的主要方式，灌區(qū)田面高程分析建立在該模型的基礎(chǔ)上。結(jié)合周橋灌區(qū)實(shí)測(cè)高程點(diǎn)數(shù)據(jù)，圖形法建立周橋灌區(qū)DEM 可分為2 步：第一步，對(duì)灌區(qū)內(nèi)多個(gè)離散點(diǎn)的高程數(shù)據(jù)進(jìn)行空間差值從而構(gòu)成不規(guī)則網(wǎng)絡(luò)（TIN）模型；第二步，在不規(guī)則網(wǎng)絡(luò)模型的基礎(chǔ)上生成規(guī)則網(wǎng)格（Grid）模型。

1.2.2 自灌率模型構(gòu)建

本文將水源水位一定條件下的灌區(qū)自流灌溉面積占灌溉總面積的百分比定義為自灌率[18]。干渠渠首水位扣除沿程水頭損失，以及沿程各種渠系建筑物的局部水頭損失，推算出農(nóng)渠的控制點(diǎn)水位（閘后水位），根據(jù)控制點(diǎn)水位結(jié)合DEM 進(jìn)行空間分析，提取控制單元內(nèi)各點(diǎn)高程得出相鄰等高線之間柵格數(shù)目，計(jì)算出各渠道控制范圍內(nèi)對(duì)應(yīng)的耕地面積，即自流灌溉面積。本文以DEM 為基礎(chǔ)結(jié)合Logistic 函數(shù)，得到周橋灌區(qū)不同地塊的自灌率模型。

1.2.3 渠系工作制度模擬

周橋灌區(qū)的渠系調(diào)度方式不同于常規(guī)灌區(qū)以某一級(jí)渠道作為輪灌組，而是通過(guò)節(jié)制閘劃分輪灌組，本文根據(jù)灌區(qū)內(nèi)設(shè)立的22 個(gè)節(jié)制閘將灌區(qū)劃分為16個(gè)灌溉單元（圖1），每個(gè)灌溉單元內(nèi)實(shí)行續(xù)灌。根據(jù)灌區(qū)灌溉單元，選擇分條輪灌、分段輪灌、分片輪灌以及現(xiàn)行輪灌4 種渠系工作制度，通過(guò)自灌率模型模擬不同渠系工作制度下灌區(qū)自流灌溉情況，尋求周橋灌區(qū)最優(yōu)渠系工作制度。

2 結(jié)果與分析

2.1 DEM 率定與驗(yàn)證

根據(jù)統(tǒng)計(jì)學(xué)理論，采用相關(guān)系數(shù)R、平均偏差MBE 等2個(gè)指標(biāo)對(duì)周橋灌區(qū)ASTER GDEM 數(shù)據(jù)進(jìn)行評(píng)估，通過(guò)最小二乘法[19]用實(shí)測(cè)數(shù)據(jù)率定ASTER GDEM 數(shù)據(jù)，計(jì)算式為：

式中：H0為ASTER GDEM 數(shù)據(jù)（m）；H1為實(shí)測(cè)數(shù)據(jù)（m）；a、b 為回歸系數(shù)。

據(jù)統(tǒng)計(jì)，周橋灌區(qū)實(shí)測(cè)高程在5.4～12.8 m 之間，ASTER GDEM 數(shù)據(jù)在4.0～12.7 m 之間，二者擬合關(guān)系見(jiàn)圖2。經(jīng)計(jì)算，ASTER GDEM 數(shù)據(jù)與實(shí)測(cè)數(shù)據(jù)的相關(guān)系數(shù)R 達(dá)到0.96，MBE 為0.27 m，2 組數(shù)據(jù)偏差較小，相關(guān)性高。本次驗(yàn)證的樣本數(shù)據(jù)共349 個(gè)，其中有327 個(gè)數(shù)據(jù)的ASTER GDEM 高程值小于實(shí)測(cè)值，表明周橋灌區(qū)ASTER GDEM 數(shù)據(jù)整體偏小，與郭笑怡等[20]研究結(jié)果相似，但周橋灌區(qū) ASTER GDEM 數(shù)據(jù)與實(shí)測(cè)值線性相關(guān)系數(shù)R 達(dá)到0.96，說(shuō)明二者之間存在緊密的相關(guān)關(guān)系。通過(guò)最小二乘法構(gòu)建ASTER GDEM 與實(shí)測(cè)值的線性關(guān)系，計(jì)算式為：

根據(jù)此關(guān)系式可得到率定后的周橋灌區(qū)數(shù)字高程模型，其與周橋灌區(qū)實(shí)際高程基本吻合，能夠反映出不同地塊的高程情況。因此，DEM 可用于模擬周橋灌區(qū)自流灌溉面積與渠首閘后水位的變化關(guān)系。

圖2 ASTER GDEM 與實(shí)測(cè)數(shù)據(jù)線性回歸 Fig.2 Linear regression of Aster GDEM and measured data

2.2 各灌溉單元自灌率模型

以地塊#1 為例計(jì)算該灌溉單元的自灌率，自灌率與渠首閘后水位關(guān)系見(jiàn)圖3。自灌率隨閘后水位的上升呈先慢后快又變慢的趨勢(shì)，符合Logistic 曲線分布，Logistic 曲線的數(shù)學(xué)模型[21]如下：

式中：y（H）為自灌率（%）；H 為閘后水位（m）；a、r、k 為根據(jù)實(shí)際地形信息確定的參數(shù)（k>0、r>0、a∈R）。

圖3 地塊#1 自灌率與閘后水位關(guān)系曲線 Fig. 3 Relation curve between self-irrigation rate and water level behind sluice in plot 1

表1 自灌率模型相關(guān)參數(shù) Table 1 Relevant parameters of self-irrigation rate model

表1為各地塊的自灌率模型相關(guān)參數(shù)。由表1可知，各地塊自灌率與閘后水位的決定系數(shù)R2均在0.74 以上，擬合度較高。因此，利用Logistic 曲線建立自灌率模型可模擬灌區(qū)不同渠系工作制度下自流灌溉情況。

2.3 渠系工作制度模擬及比較

洪澤湖水位變化情況：水稻灌溉季節(jié)最高水位13.59 m，最低水位12.01 m，平均水位為12.80 m；本文以洪澤湖的平均水位12.80 m 為例，灌水周期為6～8 d，起始時(shí)間約為7 月21 日—8 月10 日，該時(shí)間段周橋灌區(qū)用水限額為0.583 億m3，流量上限為30 m3/s。設(shè)計(jì)分條輪灌、分片輪灌以及分段輪灌等3種渠系工作制度，代入模型進(jìn)行自灌率模擬并與現(xiàn)行輪灌進(jìn)行對(duì)比。分條輪灌以支渠或分干渠控制區(qū)域?yàn)閯澐忠罁?jù)，即以每條分干渠/支渠為一個(gè)輪灌組將灌區(qū)劃分為6 個(gè)輪灌組，見(jiàn)圖4（a）；分片輪灌按分干渠或支渠將灌區(qū)劃分為5 個(gè)集中片區(qū)，每個(gè)片區(qū)作為1 個(gè)輪灌組，見(jiàn)圖4（b）；分段輪灌是按分干渠/支渠將灌區(qū)分為3 段輪灌，每段1 個(gè)輪灌組，見(jiàn)圖4（c）；現(xiàn)行輪灌是周橋灌區(qū)現(xiàn)行輪灌方式，灌區(qū)劃分為5 個(gè)片區(qū)，每個(gè)片區(qū)作為1 個(gè)輪灌組，見(jiàn)圖4（d）。

不同渠系工作制度的自灌率模擬結(jié)果見(jiàn)表2。基于自流灌溉面積并結(jié)合輸水損失、渠系水利用系數(shù)等輔助指標(biāo)進(jìn)行對(duì)比，確定最優(yōu)渠系工作制度。其中渠系輸水損失、渠系水利用系數(shù)可表示為：

式中：Vs為該渠道毛配水流量（m3/s）；t1、t0分別為灌水結(jié)束時(shí)間和灌水開(kāi)始時(shí)間（h）；ηs為支渠及以上渠系水利用系數(shù)（無(wú)量綱）；q支凈為支渠的凈流量（m3/s）；q干毛為支渠的凈流量（m3/s）。

由表2 可知，4 種渠系工作制度的自流灌溉面積占比均達(dá)到96%以上。對(duì)比4 種渠系工作制度，現(xiàn)行輪灌的自流灌溉面積占比為96.59%，分片輪灌的自流灌溉面積和自灌率最大，為20.67 萬(wàn)hm2和97.11%，分段輪灌和分條輪灌自流灌溉面積基本相同，基于自灌率最大的目標(biāo)，4 種渠系工作制度以分片輪灌為最優(yōu)。對(duì)比渠道輸水損失，相對(duì)現(xiàn)行輪灌，分片輪灌和分段輪灌的渠道輸水損失分別減少了 45.37%、13.56%，而分條輪灌的輸水損失增加了3.34%。對(duì)比渠系水利用系數(shù)，分條輪灌和現(xiàn)行輪灌渠系水利用系數(shù)差別不大，均在0.65 左右，分段輪灌、分片輪灌的渠系水利用系數(shù)分別為0.69、0.78，較現(xiàn)行輪灌分別提高了4.55%、18.18%，均優(yōu)于現(xiàn)行輪灌，其中分片輪灌最優(yōu)，分片輪灌輸水損失小且提高了灌區(qū)的渠系輸水效率。

圖4 4 種渠系工作制度概化圖 Fig.4 Generalized diagram of channel work system

表2 渠系工作制度模擬結(jié)果 Table 2 Simulation results of channel work system

3 討 論

合理的渠系工作制度對(duì)自流灌區(qū)實(shí)現(xiàn)自流灌溉面積最大化十分重要。本文利用DEM 和Logistic 函數(shù)建立周橋灌區(qū)自灌率模型，將自流灌區(qū)自灌率與水源水位建立數(shù)學(xué)關(guān)系，為模擬灌區(qū)渠系工作制度提供了重要方法。

數(shù)字高程模型是實(shí)現(xiàn)灌區(qū)自流灌溉面積定量估算的有效手段。ASTER GDEM 數(shù)據(jù)在周橋灌區(qū)整體偏小，其值與實(shí)測(cè)數(shù)據(jù)相關(guān)系數(shù)達(dá)到0.96，呈極顯著相關(guān)。夏濤等[22]、惠鳳鳴等[23]研究也表明，ASTER GDEM 數(shù)據(jù)與實(shí)際高程存在一定誤差，但是ASTER GDEM 數(shù)據(jù)與實(shí)測(cè)數(shù)據(jù)相關(guān)性高，故利用ASTER 數(shù)據(jù)提取DEM 經(jīng)率定后可用于模擬不同閘后水位的灌區(qū)自流灌溉情況。相對(duì)傳統(tǒng)單點(diǎn)測(cè)量的形式，無(wú)人機(jī)測(cè)量能提供面上的數(shù)據(jù)，且無(wú)人機(jī)的高程測(cè)量精度最高可達(dá)到±5 cm 的水平，另外無(wú)人機(jī)的特性也保證了其空間分辨率較衛(wèi)星數(shù)據(jù)有巨大的提升。因此為了提高自灌率模型的精準(zhǔn)度，無(wú)人機(jī)將會(huì)是下一步灌區(qū)灌溉合理化的研究方向。

通過(guò)周橋灌區(qū)數(shù)字高程模型，結(jié)合Logistic 函數(shù)對(duì)閘后水位和自灌率的關(guān)系進(jìn)行擬合，建立周橋灌區(qū)自灌率模型，發(fā)現(xiàn)各地塊閘后水位與自灌率的決定系數(shù)R2均在0.74 以上，擬合度高，該模型可用于模擬自流灌區(qū)不同渠系工作制度的自流灌溉情況，這與宋靜茹等[18]研究結(jié)果一致。將分條輪灌、分段輪灌、分片輪灌以及現(xiàn)行輪灌等4 種渠系工作制度分別代入該模型進(jìn)行模擬，確定不同工況下周橋灌區(qū)的自灌率，并結(jié)合輸水損失和渠系水利用系數(shù)等參數(shù)進(jìn)行比選，其中分條輪灌和現(xiàn)行輪灌的3 個(gè)指標(biāo)之間無(wú)顯著差異，分段輪灌的輸水損失以及渠系水利用系數(shù)優(yōu)于現(xiàn)行輪灌，但其自流灌溉面積較現(xiàn)行輪灌降低1.46%，而分片輪灌的自流灌溉面積占比達(dá)到97.11%，輸水損失為172.74 萬(wàn)m3，渠系水利用系數(shù)為0.78，較現(xiàn)行輪灌工作制度的自灌率和渠系水利用系數(shù)提高了5.84%、18.18%，輸水損失降低了45.37%，提高了灌區(qū)自流灌溉面積且輸水損失小、渠系水利用系數(shù)高，為最優(yōu)的灌區(qū)渠系工作制度。

4 結(jié) 論

1）ASTER GDEM 數(shù)據(jù)經(jīng)實(shí)測(cè)值率定后，與周橋灌區(qū)的實(shí)際高程基本吻合，可建立周橋灌區(qū)DEM。

2）不同地塊自灌率模型擬合度較高，決定系數(shù)均達(dá)到0.74 以上，可用自灌率模型模擬周橋灌區(qū)不同渠系工作制度自流灌溉情況。

3）將4 種渠系工作制度分別帶入該模型進(jìn)行模擬，4 種渠系工作制度中分片輪灌為最優(yōu)渠系工作方式，可將分片輪灌結(jié)合自灌率模型應(yīng)用到周橋灌區(qū)灌溉系統(tǒng)中，提高灌區(qū)的灌溉水平。