王理想, 蔡明玉, 白雪蓮, 季樹新, 賈 飚, 馮學武, 常學禮

(1.魯東大學 資源與環(huán)境工程學院,山東 煙臺 264025; 2.煙臺工程職業(yè)技術學院,山東 煙臺 264000; 3.內(nèi)蒙古自治區(qū)水利科學研究院,內(nèi)蒙古 呼和浩特 010051)

黃河是我國北方重要水資源供給來源之一,為保障黃河干流保持最低限度行水量要求,水利部對沿黃各省(自治區(qū))引黃水量做了嚴格的限定[1]。在最近的流域用水安排中,內(nèi)蒙古河套灌區(qū)用水量要控制在4.0×109m3以內(nèi)[2]。如果這一剛性要求正式實施,則對河套灌區(qū)種植規(guī)模(引黃灌溉面積)將起到直接影響,因為近些年來河套灌區(qū)引黃水量在4.3×109～4.8×109m3[3],加上一干渠黃河干流直接引水和分凌水量,初步估計超額用水至少5.0×108m3。因此,在引黃水量最新控制標準下,如何維持河套灌區(qū)糧食生產(chǎn)安全和區(qū)域生態(tài)系統(tǒng)安全,摸清近幾十年來河套灌區(qū)有效灌溉面積成為調(diào)整灌域內(nèi)種植結(jié)構(gòu)、推廣節(jié)水灌溉規(guī)模以及實施退耕還草(林)等國計民生重大措施的主要依據(jù)。

從相關水資源合理利用研究現(xiàn)狀來看,針對水資源短缺、水質(zhì)污染等導致的生態(tài)安全研究成為其重要組成部分[4-5],特別是在干旱、半干旱地區(qū)已成為當前生態(tài)水文學研究熱點[6-7]。從研究內(nèi)容特點來看,主要表現(xiàn)在水資源現(xiàn)狀與變化趨勢[8]、區(qū)域生態(tài)需水[9]、生態(tài)系統(tǒng)水資源調(diào)控模式[10]以及水資源承載力[11]等四個方面。從研究區(qū)域特色來看主要集中在: (1) 受冰山融雪影響較大,季節(jié)性變化明顯的青藏高原地區(qū); (2) 處在干旱區(qū)的西部內(nèi)陸河流域; (3) 介于干旱和半干旱地區(qū)的特大型人工自流灌區(qū)[12-14]。其中,內(nèi)蒙古河套灌區(qū)既是我國北方重要的糧食生產(chǎn)基地之一,也是亞洲最大一首制自流引水灌溉區(qū)[15],該區(qū)域生態(tài)安全問題主要就表現(xiàn)在水資源供給關系和灌溉與鹽漬化防治等方面[16-18]。

從灌溉面積確認角度來看,整個河套灌區(qū)為灌溉農(nóng)業(yè)。巴彥淖爾市近十年國民經(jīng)濟統(tǒng)計年鑒數(shù)據(jù)顯示農(nóng)業(yè)種植面積在610 970～730 533 hm2,中國科學院資源環(huán)境科學數(shù)據(jù)中心土地利用/覆蓋數(shù)據(jù)表明同期喬木林(需要灌溉維持)地面積在20 479～27 100 hm2,二者合計面積為631 449～753 137 hm2; 而同期河套灌區(qū)管理局公布全域灌溉面積為574 000 hm2[19]。二者之間存在較大差異,說明河套灌區(qū)農(nóng)業(yè)種植存在面積較大的井水、井水+引黃灌溉等多種方式[20]。因此,在多種灌溉形式存在的條件下如何識別歷年引黃灌溉面積,成為水量控制條件下種植結(jié)構(gòu)或規(guī)模調(diào)整的重要依據(jù)。

從引黃灌溉面積確認的技術路線來看,河套灌區(qū)灌溉制度一年分春、夏和秋季灌溉。其中夏灌是農(nóng)作物生長期,在遙感影像上區(qū)分引黃灌溉和井水、井水+引黃灌溉具有難度,而春秋灌溉是在農(nóng)作物播種前和收獲后進行,灌溉后農(nóng)田在一定時間(春灌5天左右,秋灌7天左右)被水體覆蓋,在遙感識別上易于甄別。同時,在方法上基于遙感手段提取水體技術成熟[21-23],而且MODIS數(shù)據(jù)時間分辨率高的特點也支持對春、秋灌溉進程監(jiān)測。

據(jù)此,本文擬以MODIS日數(shù)據(jù)為依托,以水體識別方法為手段,以統(tǒng)計學方法為補充,以近二十年來河套地區(qū)灌溉面積估算與變化趨勢分析為目的,研究河套灌溉面積時空變異特征,為灌區(qū)水資源合理利用和農(nóng)業(yè)種植結(jié)構(gòu)、規(guī)模調(diào)整提供技術和信息支持。

1 研究區(qū)概況

內(nèi)蒙古河套灌區(qū)位于黃河中上游內(nèi)蒙古段北岸的沖積平原(圖1),其行政區(qū)劃屬于內(nèi)蒙古自治區(qū)巴彥淖爾市,主要包括磴口縣、杭錦后旗、臨河區(qū)、五原縣、烏拉特前旗等五個旗縣區(qū)。地理坐標106°18′43″～109°41′42″E,40°12′50″～41°18′25″N。該區(qū)屬于溫帶大陸性氣候,降水量150～400 mm,年均溫5.6～7.4 ℃,10 ℃以上積溫3 000～3 280 ℃,無霜期130～150 d,是典型的干旱、半干旱區(qū)交錯區(qū)。河套灌區(qū)春灌時段在每年的4月至6月,秋灌則主要集中在每年的9月中旬至11月。因河套地區(qū)鹽堿化嚴重,秋灌溉具有淋鹽、壓堿為來年春季播種提供適宜播種條件的作用。

圖1 研究區(qū)位置圖與解譯精度野外檢查點

2 數(shù)據(jù)與方法

2.1 數(shù)據(jù)源

研究選取MODIS和資源三號衛(wèi)星數(shù)據(jù)分別進行灌溉面積及研究區(qū)范圍識別和確定。MODIS數(shù)據(jù)采用美國USGS遙感數(shù)據(jù)庫提供的MOD09GQ 250 m地表反射率產(chǎn)品(軌道號h26v04)。經(jīng)過對河套灌區(qū)管理局灌溉調(diào)度記錄查詢,研究區(qū)春灌、秋灌起始時間為每年4月1日至5月31日和9月15日至11月15日,在2000-2018年春灌秋灌期間河套灌區(qū)無云有效數(shù)據(jù)共774景。

2.2 研究方法

2.2.1 地面反射光譜測定與閾值確定 地面反射光譜測定是在2018年9月29日至10月16日之間進行,選擇晴朗無云、日照充足且穩(wěn)定的時段(10:00～14:00)進行典型地物光譜反射測定[24]。儀器采用波長范圍為325～1 075 nm的手持便攜式光譜儀(HandHeld-2),測定前對儀器充分預熱和校光處理,測定時探頭視場角設定為25°,探頭垂直向下并與實測地物保持0.67 m的高度差,每隔15 min重新校光。野外光譜測定工作時,對同一地物連續(xù)測10組數(shù)據(jù),取10次所測光譜平均值為最終光譜。測定下墊面分別為灌溉積水、濕潤裸地、鹽堿地表、植被和各種留茬地表(表1)。

表1 各類典型地物NDVI閾值

2.2.2 MODIS數(shù)據(jù)計算流程 首先下載有效MODIS數(shù)據(jù)(250 m×250 m),剔除研究區(qū)有條帶噪音或云覆蓋的影像,經(jīng)MRT數(shù)據(jù)轉(zhuǎn)換并導出1、2波段; 然后在ArcGIS平臺內(nèi)計算NDVI指數(shù),根據(jù)地面遙感野外實測的典型地物NDVI閾值提取水體。根據(jù)資源三號衛(wèi)星高分影像和河套灌區(qū)管理局提供的2017年內(nèi)蒙古河套灌區(qū)現(xiàn)狀圖,確定研究區(qū)范圍及黃河邊界并進行掩膜處理。最后,利用逐年非灌溉期MODIS 13Q1數(shù)據(jù)疊加提取湖泊、水渠等水域并在對應年份剔除。

2.2.3 精度檢查 1960年,Cohen提出Kappa系數(shù)來評價生成圖件與地物的一致性,隨后在此基礎上發(fā)展了標準Kappa系數(shù)、隨機Kappa系數(shù)、位置Kappa系數(shù)等用來評價遙感圖像解譯的正確程度[25]。本文采用隨機Kappa系數(shù)(K)分析解譯后影像的精度,公式為

(1)

其中:p1為相同位置經(jīng)遙感解譯整合后的生成圖件與標記的水澆地類型一致的點數(shù);p2與p1相反,為存在差異的點數(shù);p3為所有樣點點數(shù)總和。在精度檢查中,檢查點為2017年10月22日至10月25日和2018年 10月20日至10月24日野外GPS標定點,隨機選擇在水澆地邊緣,樣點數(shù)為76個,水澆地分類精度檢查結(jié)果見表 2,K值為0.914,大于0.75(表2)。

表2 灌溉農(nóng)田解譯精度檢查

2.2.4 灌溉面積估算 由于影像獲取條件所限(有云數(shù)據(jù)較多)導致研究時段內(nèi)有效影像缺失,遙感監(jiān)測獲取面積必定偏小的事實,本文將通過統(tǒng)計學方法對缺失時段內(nèi)面積進行補充修正。具體方法為在灌溉時段求算所有相鄰兩天有效數(shù)據(jù)灌溉面積增量,然后從增量值中提取50個樣本值求平均,分別計算出春、秋灌溉每天平均增量值(Ac j),最后采用下列公式計算年引黃灌溉面積

Si j=(Tj-Nj)×Ac j+Ai j,

(2)

其中:Si j為春秋引黃灌溉面積;Tj為春秋灌溉天數(shù),Nj為有效遙感影像天數(shù);Ai j為有效遙感影像識別所獲取的春灌秋澆控制面積;i=1,2,…,19,分別代表2000年-2018年;j=1,2,分別代表春灌和秋灌。

此外,為了更好驗證統(tǒng)計學方法補充河套灌區(qū)引黃灌溉面積的準確性以及對變化特征進行分析,分別選取研究時段最小(2003年)、最接近平均(2010年)和最大(2016年)年份進行分析。

3 結(jié)果與分析

3.1 河套灌區(qū)灌溉面積變化與時間尺度特征分析

3.1.1 春季灌溉面積 2000-2018年河套灌區(qū)春灌面積呈增加趨勢,斜率為9 852.9(Sig<0.001,圖2a1)。變化特點以2008年為分界點,分為平穩(wěn)波動和迅速升高兩個階段。在平穩(wěn)波動階段,灌溉面積變化幅度小,極差為16 824 hm2,在迅速升高階段,極差為166 665 hm2,是前一階段的10倍。為了彌補影像缺失導致對春季灌溉面積估算偏小的事實,采用公式(1)計算補充的結(jié)果表明(圖2a2),補充后春灌面積變化特征與遙感獲取面積變化趨勢一致,但變化特點有所區(qū)別。在平穩(wěn)波動階段,灌溉面積極差為52 364.3 hm2; 而在迅速升高階段春灌面積極差為152 041.9 hm2,是前一階段的2.9倍。比較遙感監(jiān)測(圖2a1)和遙感監(jiān)測+統(tǒng)計補充(圖2a2)結(jié)果可以看出,年平均春灌面積相差149 111 hm2。

3.1.2 秋季灌溉面積 從遙感監(jiān)測秋灌面積變化特點來看(圖2b1),在研究時期內(nèi)增加趨勢不明顯(斜率為1 690.3,Sig=0.488),最小為91 350 hm2(2002年),最大為268 448 hm2(2004年),近20年間極差為177 098 hm2。同樣采用公式(1)對遙感數(shù)據(jù)缺失天數(shù)進行補充計算(圖2b2)。結(jié)果表明,在研究時期內(nèi)變化呈下降趨勢,斜率為-419.06(Sig=0.842)。從秋灌面積變化特征來看,平均值為367 623 hm2,最小為292 463 hm2(2015年),最大為471 097.3 hm2(2012年),變幅為178 635 hm2。從遙感獲取面積(圖2b1)和統(tǒng)計補充面積(圖2b2)的比較來看,年平均相差185 431 hm2。此外,在疊加分析中發(fā)現(xiàn)每年春、秋灌溉都有重復灌溉,重復面積2 406～44 075 hm2,平均為18 495 hm2。從趨勢來看,近二十年重復灌溉面積呈擴大趨勢,2007年以前重復面積較小,均不超過9 000 hm2,從2008年開始重復灌溉明顯增加,平均在10 000 hm2以上。從典型年份來看最小年份(2004年)重復灌溉面積2 406 hm2,最大年份(2012年)達到40 119 hm2,是最小年份的16.7倍。

圖2 河套灌區(qū)灌溉面積變化特征 (a 春灌面積, b 秋灌面積; 1 監(jiān)測面積, 2 統(tǒng)計補充面積)

3.2 河套灌區(qū)灌溉面積空間分異

河套灌區(qū)從西向東由沈烏灌域、解放閘灌域、永濟灌域、義長灌域和烏拉特灌域組成(圖1)。為了探究引黃灌溉面積空間分異特征,分別選取灌溉面積最小、最接近平均和最大的2001年、2010年以及2013年三個典型年份進行分析(圖3,表3)。

圖3 典型年份河套灌區(qū)引黃灌溉區(qū)域分布

表3 典型年份分灌域引黃灌溉面積

從春季引黃灌面積特點來看(圖3,表3),在最小引黃灌溉年,依次為烏拉特灌域(6 469 hm2)、義長灌域(6 356 hm2)、沈烏灌域(5 338 hm2)、永濟灌域(1 181 hm2)、解放閘灌域(894 hm2),其中引黃灌溉面積最大的烏拉特灌域較最小的解放閘灌域多5 575 hm2。在最接近平均年,義長灌域和烏拉特灌域的引黃灌溉面積最大,分別為33 988 hm2和26 656 hm2,沈烏灌域(5 750 hm2)與解放閘灌域(3 294 hm2)灌溉區(qū)域分布居中,永濟灌域(2 275 hm2)最小(僅為義長灌域7%左右)。在引黃灌溉面積最大年,各個灌域間的空間差異進一步擴大,灌溉面積最大的義長灌域(74 088 hm2)已接近永濟灌域(3 838 hm2)的19倍,灌溉面積第二的烏拉特灌域較第三的沈烏灌域灌溉面積差異也擴大至37 031 hm2。

從秋季引黃灌溉面積排序特點來看(圖3,表3),在最小年份,灌溉面積依次為義長灌域(48 388 hm2)、解放閘灌域(31 594 hm2)、烏拉特灌域(14 369 hm2)、永濟灌域(11 713 hm2)、沈烏灌域(7 019 hm2),分別占13%、43%、10%、28%、6%,灌域間面積存在較大差別。在平均年份,灌溉面積最大的義長灌域與最小的沈烏灌域相差達到90 719 hm2(義長灌域為115 175 hm2,沈烏灌域為6 325 hm2),排序居中的烏拉特灌域、永濟灌域和解放閘灌域的引黃灌溉面積為29 000～42 000 hm2。

3.3 典型年份灌溉進程

為準確理解研究時段內(nèi)引黃灌溉面積年內(nèi)累積進程,以及統(tǒng)計補充計算灌溉面積與遙感監(jiān)測面積之間關系,同樣選取上述三個典型年份,分別對其遙感監(jiān)測面積和統(tǒng)計補充計算面積累積進程趨勢以及二者之間關系進行一元回歸分析,選擇灌溉進程最佳表達模型,解釋河套灌區(qū)年灌溉面積估測的可信性和變化特點。

首先從最小引黃灌溉年分析結(jié)果來看(圖4a),隨灌溉時間增加遙感監(jiān)測引黃灌溉面積和統(tǒng)計補充計算灌溉面積變化趨勢表現(xiàn)出較高的相似性(圖4a1,a2)。一元回歸模型最佳表達為指數(shù)函數(shù)。其中遙感監(jiān)測累積灌溉面積進程模型為y=3 536.3 e0.043 3x(R=0.960>R0.001,33=0.532),統(tǒng)計補充計算灌溉面積累積進程模型為y=21 546 e0.045 2x(R=0.957>R0.001,33=0.532)。相關系數(shù)顯著性都達到了極顯著水平。從變化特點來看,引黃灌溉累積面積在90天以前緩慢增加,90天以后累積面積迅速上升; 統(tǒng)計補充計算灌溉面積也具有相同的變化趨勢。

圖4 年灌溉面積進程趨勢(a 最小灌溉年份(2001), b 最接近平均灌溉年份(2010),c 最大灌溉年份(2013); 1 監(jiān)測面積, 2 修正面積, 3 監(jiān)測與統(tǒng)計補充面積關系)

其次從最接近平均灌溉面積年分析結(jié)果來看(圖4b),隨灌溉時間增加遙感監(jiān)測引黃灌溉面積和統(tǒng)計補充計算灌溉面積變化趨勢同樣表現(xiàn)出較高的相似性(圖4b1,b2)。一元回歸模型最佳表達為指數(shù)函數(shù)。其中遙感監(jiān)測累積引黃灌溉面積進程模型為y=26 171 e0.027 8x(R=0.956>R0.001,34=0.525),統(tǒng)計補充計算灌溉面積累積進程模型為y=25 597 e0.045 9x(R=0.976>R0.001,34=0.525)。相關系數(shù)顯著性都達到了極顯著水平。從變化特點來看,遙感監(jiān)測灌溉累積面積和統(tǒng)計補充計算灌溉面積具有相同的變化趨勢,均呈現(xiàn)增加趨勢。

最后從最大引黃灌溉年分析結(jié)果來看(圖4c),隨灌溉時間增加遙感監(jiān)測引黃灌溉面積和統(tǒng)計補充計算灌溉面積變化趨勢表現(xiàn)出較高的相似性(圖4c1,c2)。一元回歸模型最佳表達為冪函數(shù)。其中遙感監(jiān)測累積灌溉面積進程模型為y=24 361x0.801 2(R=0.986>R0.001,41=0.484),統(tǒng)計補充計算灌溉面積累積進程模型為y=16 477x1.424 6(R=0.985>R0.001,41=0.484)。相關系數(shù)顯著性都達到了極顯著水平。從變化特點來看,遙感監(jiān)測引黃灌溉累積面積和統(tǒng)計補充計算灌溉面積均表現(xiàn)為階梯式增加趨勢,在59天以前呈快速增加趨勢,59至108天時段內(nèi)增速較慢,108天以后再次呈現(xiàn)快速增長趨勢。

此外對三個典型年份遙感監(jiān)測引黃灌溉面積和統(tǒng)計補充計算灌溉面積二者之間關系分析表明(圖4a3,b3,c3),它們之間都存在極顯著相關,相關系數(shù)顯著性檢查分別為R=0.975>R0.001,33=0.532、R=0.974>R0.001,34=0.525和R=0.970>R0.001,41=0.484。綜上所述,河套灌區(qū)年引黃灌溉面積進程無論采用統(tǒng)計補充計算灌溉面積還是遙感監(jiān)測面積其變化趨勢規(guī)律具有相似性(圖4)。三個典型年份都可以用顯著性極高的一元一次回歸模型表達,其中最小年和最接近平均年為指數(shù)函數(shù),最大年為冪函數(shù)。三個典型年份引黃灌溉面積和統(tǒng)計補充計算灌溉面積回歸分析表明,統(tǒng)計補充計算灌溉面積變化趨勢具有較高可信性,即運用公式(2)可以較好地補充影像缺失導致面積偏小的問題。

4 討論

河套灌域農(nóng)業(yè)種植主要依賴引黃灌溉,其中當?shù)赝寥利}堿化水平偏高的自然現(xiàn)象,使農(nóng)作物收獲的秋季灌溉不僅具有增加農(nóng)田土壤水分的作用,而且對于土壤鹽堿含量通過土壤下滲排除也具有決定性作用[16,19]。與其他自然條件相似地區(qū)相比較,河套灌域秋季灌溉規(guī)模也是最大,這既與黃河水文特征有關,也與消減河套地區(qū)土壤鹽堿化對農(nóng)業(yè)種植影響密切相關[17-18]。從本文研究結(jié)果來看,多年秋灌面積平均為367 623 hm2,是春灌面積(219 579 hm2)的1.67倍。從變化趨勢來看這一差距逐漸變小,從2010年開始春灌面積呈明顯增加趨勢(圖2a1,a2)。從導致春秋灌溉變化趨勢的原因來看,一是特定土壤鹽堿含量較高需要灌溉排鹽排堿[18],二是市場需求導致的近十年來向日葵種植面積增加[15]。因為在向日葵種植品種中,生長期短的品種(油葵和西葫蘆,春灌導致的播種期后延不影響產(chǎn)量)種植面積逐年增加。

從河套灌區(qū)年引黃灌溉面積進程來看,無論采用統(tǒng)計補充計算灌溉面積還是遙感監(jiān)測面積,其變化趨勢規(guī)律具有相似性(圖4)。在三個典型年份分析中,都可以用顯著性極高的一元一次回歸模型表達,其中最小年和最接近平均年為指數(shù)函數(shù),最大年為冪函數(shù)。遙感監(jiān)測引黃灌溉面積和統(tǒng)計補充計算灌溉面積關聯(lián)分析表明(圖4a3,b3,c3),二者年灌溉進程具有極高相似性(R>Rs,n-1),即運用公式(2)可以較好地補充影像缺失導致面積偏小的問題。

在年灌溉面積疊加分析中發(fā)現(xiàn),每年春、秋灌溉都有重復灌溉現(xiàn)象發(fā)生,重復面積為2 406～44 075 hm2,平均為18 495 hm2。從其趨勢來看,近二十年重復灌溉面積呈擴大趨勢,在2007年以前重復面積較小,均不超過9 000 hm2,從2008年開始重復灌溉明顯增加,平均在10 000 hm2以上。從典型年份來看最小年份(2004年)重復灌溉面積2 406 hm2,最大年份(2012年)達到40 119 hm2,是最小年份的16.7倍。重復灌溉面積逐年增加的原因與生長期短的油葵種和西葫蘆植面積增加有關,因為春季灌溉導致的播種期滯后,對它們的產(chǎn)量無影響。

5 結(jié)論

(1) 在2000-2018年間河套灌區(qū)春季引黃灌面積波動大于秋季,二者變異系數(shù)分別為24.4%和13.1%。引黃灌溉面積研究期間內(nèi)變化在103 652～385 071 hm2,其中多年春灌面積平均為70 468 hm2,秋灌平均為182 192 hm2,秋灌大于春灌111 742 hm2。從統(tǒng)計補充計算面積變化特點來看,在研究時期內(nèi)變化在486 040～709 655 hm2,其中多年春灌面積平均為219 579 hm2,秋灌面積為367 622 hm2,秋灌大于春灌148 044 hm2。

(2) 灌溉面積變化特點在時間尺度上具有顯著差異,春季引黃灌面積以2008年為分界點分平穩(wěn)波動和迅速升高個階段。2000-2008年間小幅波動,2008年以后呈快速增加趨勢,前后增幅相差105 315 hm2,總體變異系數(shù)為10.9%。秋季引黃灌面積在2000-2018年間面積變化呈減少趨勢,變異系數(shù)為13.1%。

(3) 春、秋灌溉面積空間分異顯著,義長和烏拉特灌域春、秋灌溉面積分布都比較廣; 沈烏灌域春、秋灌溉面積分布較小,面積均不超過15 000 hm2; 永濟灌域以秋灌為主,春灌分布最少,2010年春、秋灌溉面積相差約34 000 hm2。

(4) 典型年份累積監(jiān)測灌溉和統(tǒng)計補充計算灌溉面積進程具有相似變化趨勢,二者關系可用一元一次線性回歸模型估測,相關系數(shù)檢驗都在0.001水平上呈極顯著相關。