基于MODIS数据的河套地区引黄灌溉面积动态变化分析

王理想, 蔡明玉, 白雪莲, 季树新, 贾 飚, 冯学武, 常学礼

(1.鲁东大学 资源与环境工程学院,山东 烟台 264025; 2.烟台工程职业技术学院,山东 烟台 264000; 3.内蒙古自治区水利科学研究院,内蒙古 呼和浩特 010051)

黄河是我国北方重要水资源供给来源之一,为保障黄河干流保持最低限度行水量要求,水利部对沿黄各省(自治区)引黄水量做了严格的限定[1]。在最近的流域用水安排中,内蒙古河套灌区用水量要控制在4.0×109m3以内[2]。如果这一刚性要求正式实施,则对河套灌区种植规模(引黄灌溉面积)将起到直接影响,因为近些年来河套灌区引黄水量在4.3×109～4.8×109m3[3],加上一干渠黄河干流直接引水和分凌水量,初步估计超额用水至少5.0×108m3。因此,在引黄水量最新控制标准下,如何维持河套灌区粮食生产安全和区域生态系统安全,摸清近几十年来河套灌区有效灌溉面积成为调整灌域内种植结构、推广节水灌溉规模以及实施退耕还草(林)等国计民生重大措施的主要依据。

从相关水资源合理利用研究现状来看,针对水资源短缺、水质污染等导致的生态安全研究成为其重要组成部分[4-5],特别是在干旱、半干旱地区已成为当前生态水文学研究热点[6-7]。从研究内容特点来看,主要表现在水资源现状与变化趋势[8]、区域生态需水[9]、生态系统水资源调控模式[10]以及水资源承载力[11]等四个方面。从研究区域特色来看主要集中在: (1) 受冰山融雪影响较大,季节性变化明显的青藏高原地区; (2) 处在干旱区的西部内陆河流域; (3) 介于干旱和半干旱地区的特大型人工自流灌区[12-14]。其中,内蒙古河套灌区既是我国北方重要的粮食生产基地之一,也是亚洲最大一首制自流引水灌溉区[15],该区域生态安全问题主要就表现在水资源供给关系和灌溉与盐渍化防治等方面[16-18]。

从灌溉面积确认角度来看,整个河套灌区为灌溉农业。巴彦淖尔市近十年国民经济统计年鉴数据显示农业种植面积在610 970～730 533 hm2,中国科学院资源环境科学数据中心土地利用/覆盖数据表明同期乔木林(需要灌溉维持)地面积在20 479～27 100 hm2,二者合计面积为631 449～753 137 hm2; 而同期河套灌区管理局公布全域灌溉面积为574 000 hm2[19]。二者之间存在较大差异,说明河套灌区农业种植存在面积较大的井水、井水+引黄灌溉等多种方式[20]。因此,在多种灌溉形式存在的条件下如何识别历年引黄灌溉面积,成为水量控制条件下种植结构或规模调整的重要依据。

从引黄灌溉面积确认的技术路线来看,河套灌区灌溉制度一年分春、夏和秋季灌溉。其中夏灌是农作物生长期,在遥感影像上区分引黄灌溉和井水、井水+引黄灌溉具有难度,而春秋灌溉是在农作物播种前和收获后进行,灌溉后农田在一定时间(春灌5天左右,秋灌7天左右)被水体覆盖,在遥感识别上易于甄别。同时,在方法上基于遥感手段提取水体技术成熟[21-23],而且MODIS数据时间分辨率高的特点也支持对春、秋灌溉进程监测。

据此,本文拟以MODIS日数据为依托,以水体识别方法为手段,以统计学方法为补充,以近二十年来河套地区灌溉面积估算与变化趋势分析为目的,研究河套灌溉面积时空变异特征,为灌区水资源合理利用和农业种植结构、规模调整提供技术和信息支持。

1 研究区概况

内蒙古河套灌区位于黄河中上游内蒙古段北岸的冲积平原(图1),其行政区划属于内蒙古自治区巴彦淖尔市,主要包括磴口县、杭锦后旗、临河区、五原县、乌拉特前旗等五个旗县区。地理坐标106°18′43″～109°41′42″E,40°12′50″～41°18′25″N。该区属于温带大陆性气候,降水量150～400 mm,年均温5.6～7.4 ℃,10 ℃以上积温3 000～3 280 ℃,无霜期130～150 d,是典型的干旱、半干旱区交错区。河套灌区春灌时段在每年的4月至6月,秋灌则主要集中在每年的9月中旬至11月。因河套地区盐碱化严重,秋灌溉具有淋盐、压碱为来年春季播种提供适宜播种条件的作用。

图1 研究区位置图与解译精度野外检查点

2 数据与方法

2.1 数据源

研究选取MODIS和资源三号卫星数据分别进行灌溉面积及研究区范围识别和确定。MODIS数据采用美国USGS遥感数据库提供的MOD09GQ 250 m地表反射率产品(轨道号h26v04)。经过对河套灌区管理局灌溉调度记录查询,研究区春灌、秋灌起始时间为每年4月1日至5月31日和9月15日至11月15日,在2000-2018年春灌秋灌期间河套灌区无云有效数据共774景。

2.2 研究方法

2.2.1 地面反射光谱测定与阈值确定 地面反射光谱测定是在2018年9月29日至10月16日之间进行,选择晴朗无云、日照充足且稳定的时段(10:00～14:00)进行典型地物光谱反射测定[24]。仪器采用波长范围为325～1 075 nm的手持便携式光谱仪(HandHeld-2),测定前对仪器充分预热和校光处理,测定时探头视场角设定为25°,探头垂直向下并与实测地物保持0.67 m的高度差,每隔15 min重新校光。野外光谱测定工作时,对同一地物连续测10组数据,取10次所测光谱平均值为最终光谱。测定下垫面分别为灌溉积水、湿润裸地、盐碱地表、植被和各种留茬地表(表1)。

表1 各类典型地物NDVI阈值

2.2.2 MODIS数据计算流程 首先下载有效MODIS数据(250 m×250 m),剔除研究区有条带噪音或云覆盖的影像,经MRT数据转换并导出1、2波段; 然后在ArcGIS平台内计算NDVI指数,根据地面遥感野外实测的典型地物NDVI阈值提取水体。根据资源三号卫星高分影像和河套灌区管理局提供的2017年内蒙古河套灌区现状图,确定研究区范围及黄河边界并进行掩膜处理。最后,利用逐年非灌溉期MODIS 13Q1数据叠加提取湖泊、水渠等水域并在对应年份剔除。

2.2.3 精度检查 1960年,Cohen提出Kappa系数来评价生成图件与地物的一致性,随后在此基础上发展了标准Kappa系数、随机Kappa系数、位置Kappa系数等用来评价遥感图像解译的正确程度[25]。本文采用随机Kappa系数(K)分析解译后影像的精度,公式为

(1)

其中:p1为相同位置经遥感解译整合后的生成图件与标记的水浇地类型一致的点数;p2与p1相反,为存在差异的点数;p3为所有样点点数总和。在精度检查中,检查点为2017年10月22日至10月25日和2018年 10月20日至10月24日野外GPS标定点,随机选择在水浇地边缘,样点数为76个,水浇地分类精度检查结果见表 2,K值为0.914,大于0.75(表2)。

表2 灌溉农田解译精度检查

2.2.4 灌溉面积估算 由于影像获取条件所限(有云数据较多)导致研究时段内有效影像缺失,遥感监测获取面积必定偏小的事实,本文将通过统计学方法对缺失时段内面积进行补充修正。具体方法为在灌溉时段求算所有相邻两天有效数据灌溉面积增量,然后从增量值中提取50个样本值求平均,分别计算出春、秋灌溉每天平均增量值(Ac j),最后采用下列公式计算年引黄灌溉面积

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

(2)

其中:Si j为春秋引黄灌溉面积;Tj为春秋灌溉天数,Nj为有效遥感影像天数;Ai j为有效遥感影像识别所获取的春灌秋浇控制面积;i=1,2,…,19,分别代表2000年-2018年;j=1,2,分别代表春灌和秋灌。

此外,为了更好验证统计学方法补充河套灌区引黄灌溉面积的准确性以及对变化特征进行分析,分别选取研究时段最小(2003年)、最接近平均(2010年)和最大(2016年)年份进行分析。

3 结果与分析

3.1 河套灌区灌溉面积变化与时间尺度特征分析

3.1.1 春季灌溉面积 2000-2018年河套灌区春灌面积呈增加趋势,斜率为9 852.9(Sig<0.001,图2a1)。变化特点以2008年为分界点,分为平稳波动和迅速升高两个阶段。在平稳波动阶段,灌溉面积变化幅度小,极差为16 824 hm2,在迅速升高阶段,极差为166 665 hm2,是前一阶段的10倍。为了弥补影像缺失导致对春季灌溉面积估算偏小的事实,采用公式(1)计算补充的结果表明(图2a2),补充后春灌面积变化特征与遥感获取面积变化趋势一致,但变化特点有所区别。在平稳波动阶段,灌溉面积极差为52 364.3 hm2; 而在迅速升高阶段春灌面积极差为152 041.9 hm2,是前一阶段的2.9倍。比较遥感监测(图2a1)和遥感监测+统计补充(图2a2)结果可以看出,年平均春灌面积相差149 111 hm2。

3.1.2 秋季灌溉面积 从遥感监测秋灌面积变化特点来看(图2b1),在研究时期内增加趋势不明显(斜率为1 690.3,Sig=0.488),最小为91 350 hm2(2002年),最大为268 448 hm2(2004年),近20年间极差为177 098 hm2。同样采用公式(1)对遥感数据缺失天数进行补充计算(图2b2)。结果表明,在研究时期内变化呈下降趋势,斜率为-419.06(Sig=0.842)。从秋灌面积变化特征来看,平均值为367 623 hm2,最小为292 463 hm2(2015年),最大为471 097.3 hm2(2012年),变幅为178 635 hm2。从遥感获取面积(图2b1)和统计补充面积(图2b2)的比较来看,年平均相差185 431 hm2。此外,在叠加分析中发现每年春、秋灌溉都有重复灌溉,重复面积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 河套灌区灌溉面积变化特征 (a 春灌面积, b 秋灌面积; 1 监测面积, 2 统计补充面积)

3.2 河套灌区灌溉面积空间分异

河套灌区从西向东由沈乌灌域、解放闸灌域、永济灌域、义长灌域和乌拉特灌域组成(图1)。为了探究引黄灌溉面积空间分异特征,分别选取灌溉面积最小、最接近平均和最大的2001年、2010年以及2013年三个典型年份进行分析(图3,表3)。

图3 典型年份河套灌区引黄灌溉区域分布

表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)灌溉区域分布居中,永济灌域(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 典型年份灌溉进程

为准确理解研究时段内引黄灌溉面积年内累积进程,以及统计补充计算灌溉面积与遥感监测面积之间关系,同样选取上述三个典型年份,分别对其遥感监测面积和统计补充计算面积累积进程趋势以及二者之间关系进行一元回归分析,选择灌溉进程最佳表达模型,解释河套灌区年灌溉面积估测的可信性和变化特点。

首先从最小引黄灌溉年分析结果来看(图4a),随灌溉时间增加遥感监测引黄灌溉面积和统计补充计算灌溉面积变化趋势表现出较高的相似性(图4a1,a2)。一元回归模型最佳表达为指数函数。其中遥感监测累积灌溉面积进程模型为y=3 536.3 e0.043 3x(R=0.960>R0.001,33=0.532),统计补充计算灌溉面积累积进程模型为y=21 546 e0.045 2x(R=0.957>R0.001,33=0.532)。相关系数显著性都达到了极显著水平。从变化特点来看,引黄灌溉累积面积在90天以前缓慢增加,90天以后累积面积迅速上升; 统计补充计算灌溉面积也具有相同的变化趋势。

图4 年灌溉面积进程趋势(a 最小灌溉年份(2001), b 最接近平均灌溉年份(2010),c 最大灌溉年份(2013); 1 监测面积, 2 修正面积, 3 监测与统计补充面积关系)

其次从最接近平均灌溉面积年分析结果来看(图4b),随灌溉时间增加遥感监测引黄灌溉面积和统计补充计算灌溉面积变化趋势同样表现出较高的相似性(图4b1,b2)。一元回归模型最佳表达为指数函数。其中遥感监测累积引黄灌溉面积进程模型为y=26 171 e0.027 8x(R=0.956>R0.001,34=0.525),统计补充计算灌溉面积累积进程模型为y=25 597 e0.045 9x(R=0.976>R0.001,34=0.525)。相关系数显著性都达到了极显著水平。从变化特点来看,遥感监测灌溉累积面积和统计补充计算灌溉面积具有相同的变化趋势,均呈现增加趋势。

最后从最大引黄灌溉年分析结果来看(图4c),随灌溉时间增加遥感监测引黄灌溉面积和统计补充计算灌溉面积变化趋势表现出较高的相似性(图4c1,c2)。一元回归模型最佳表达为幂函数。其中遥感监测累积灌溉面积进程模型为y=24 361x0.801 2(R=0.986>R0.001,41=0.484),统计补充计算灌溉面积累积进程模型为y=16 477x1.424 6(R=0.985>R0.001,41=0.484)。相关系数显著性都达到了极显著水平。从变化特点来看,遥感监测引黄灌溉累积面积和统计补充计算灌溉面积均表现为阶梯式增加趋势,在59天以前呈快速增加趋势,59至108天时段内增速较慢,108天以后再次呈现快速增长趋势。

此外对三个典型年份遥感监测引黄灌溉面积和统计补充计算灌溉面积二者之间关系分析表明(图4a3,b3,c3),它们之间都存在极显著相关,相关系数显著性检查分别为R=0.975>R0.001,33=0.532、R=0.974>R0.001,34=0.525和R=0.970>R0.001,41=0.484。综上所述,河套灌区年引黄灌溉面积进程无论采用统计补充计算灌溉面积还是遥感监测面积其变化趋势规律具有相似性(图4)。三个典型年份都可以用显著性极高的一元一次回归模型表达,其中最小年和最接近平均年为指数函数,最大年为幂函数。三个典型年份引黄灌溉面积和统计补充计算灌溉面积回归分析表明,统计补充计算灌溉面积变化趋势具有较高可信性,即运用公式(2)可以较好地补充影像缺失导致面积偏小的问题。

4 讨论

河套灌域农业种植主要依赖引黄灌溉,其中当地土壤盐碱化水平偏高的自然现象,使农作物收获的秋季灌溉不仅具有增加农田土壤水分的作用,而且对于土壤盐碱含量通过土壤下渗排除也具有决定性作用[16,19]。与其他自然条件相似地区相比较,河套灌域秋季灌溉规模也是最大,这既与黄河水文特征有关,也与消减河套地区土壤盐碱化对农业种植影响密切相关[17-18]。从本文研究结果来看,多年秋灌面积平均为367 623 hm2,是春灌面积(219 579 hm2)的1.67倍。从变化趋势来看这一差距逐渐变小,从2010年开始春灌面积呈明显增加趋势(图2a1,a2)。从导致春秋灌溉变化趋势的原因来看,一是特定土壤盐碱含量较高需要灌溉排盐排碱[18],二是市场需求导致的近十年来向日葵种植面积增加[15]。因为在向日葵种植品种中,生长期短的品种(油葵和西葫芦,春灌导致的播种期后延不影响产量)种植面积逐年增加。

从河套灌区年引黄灌溉面积进程来看,无论采用统计补充计算灌溉面积还是遥感监测面积,其变化趋势规律具有相似性(图4)。在三个典型年份分析中,都可以用显著性极高的一元一次回归模型表达,其中最小年和最接近平均年为指数函数,最大年为幂函数。遥感监测引黄灌溉面积和统计补充计算灌溉面积关联分析表明(图4a3,b3,c3),二者年灌溉进程具有极高相似性(R>Rs,n-1),即运用公式(2)可以较好地补充影像缺失导致面积偏小的问题。

在年灌溉面积叠加分析中发现,每年春、秋灌溉都有重复灌溉现象发生,重复面积为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倍。重复灌溉面积逐年增加的原因与生长期短的油葵种和西葫芦植面积增加有关,因为春季灌溉导致的播种期滞后,对它们的产量无影响。

5 结论

(1) 在2000-2018年间河套灌区春季引黄灌面积波动大于秋季,二者变异系数分别为24.4%和13.1%。引黄灌溉面积研究期间内变化在103 652～385 071 hm2,其中多年春灌面积平均为70 468 hm2,秋灌平均为182 192 hm2,秋灌大于春灌111 742 hm2。从统计补充计算面积变化特点来看,在研究时期内变化在486 040～709 655 hm2,其中多年春灌面积平均为219 579 hm2,秋灌面积为367 622 hm2,秋灌大于春灌148 044 hm2。

(2) 灌溉面积变化特点在时间尺度上具有显著差异,春季引黄灌面积以2008年为分界点分平稳波动和迅速升高个阶段。2000-2008年间小幅波动,2008年以后呈快速增加趋势,前后增幅相差105 315 hm2,总体变异系数为10.9%。秋季引黄灌面积在2000-2018年间面积变化呈减少趋势,变异系数为13.1%。

(3) 春、秋灌溉面积空间分异显著,义长和乌拉特灌域春、秋灌溉面积分布都比较广; 沈乌灌域春、秋灌溉面积分布较小,面积均不超过15 000 hm2; 永济灌域以秋灌为主,春灌分布最少,2010年春、秋灌溉面积相差约34 000 hm2。

(4) 典型年份累积监测灌溉和统计补充计算灌溉面积进程具有相似变化趋势,二者关系可用一元一次线性回归模型估测,相关系数检验都在0.001水平上呈极显著相关。