時啟軍
(中國南水北調集團中線有限公司河北分公司,石家莊 050000)
作為華北農田糧食的重要生產基地的河北省,灌區斗渠以下(田間)的土質農毛渠基本不具備水力學測流條件,大部分未安裝計量設施,計量設施配備率僅為2%,且完好率僅為17.6%,“按畝收費”“按時收費”“交平攤水費”的現象普遍存在,灌溉面積依靠人工統計,費時費力,且不具備客觀性,亟待利用先進技術,對灌溉面積進行統計提取。
遙感技術是獲取地表信息的一種有效手段,可從定性到定量、從宏觀到微觀地反映出地面物體(如土壤)的特性及其變化(如含水率),具有覆蓋范圍廣、時效性強、數據客觀、經濟效益高等優點[1]。因此,可以借助遙感測量,利用不同土壤濕度對光譜反射率的不同,將遙感數據進行解譯,提取灌溉面積,作為用水計量和末級渠道水費收取的參考依據。
根據光學理論研究,不同地表含水率在紅光波段及近紅外波段,表現出強吸收或強發射的特征,基于灌溉前后土壤光譜反射率有所不同,遙感數據能夠在一定層面上反映出土壤含水率的變化情況[2]。可以利用灌溉前后土壤濕度的變化,判定地塊灌溉與否,進而進行實際灌溉面積提取和灌溉進度監測[3]。
垂直干旱指數PDI (perpendicular drought index)就是利用影像光譜對不同土壤濕度的響應、基于近紅光反射率和紅光反射率構建特征空間,來表征土壤含水量的變化[4]。在Nir-Red 散點圖上,遙感影像各像素點的分布近似于一個三角形(如圖1):L 為土壤基線的法線;PDI 為直線L 的垂線,用以描述在Nir-Red 光譜特征空間上,水分含量的分布規律,離直線L 越近代表水分含量越高,越遠代表土壤水分含量越低;任取一點E,從點E 到直線L 的距離EF,即為垂直干旱指數PDI。即垂直干旱指數計算公式為:

圖1 Nir-Red 特征空間和PDI 示意圖
式中PDI 為垂直干旱指數;RRED為紅光波段反射率;RNIR為近紅外光波段的反射率;M 為土壤基線BC 的斜率。
當土壤濕度狀態未發生變化時,M 值也不會發生變化。在灌溉時,土壤會隨著灌溉,含水量信息發生變化。
計算時,為了減小或消除植被光譜對土壤含水量的影響,引入植被覆蓋度概念[5],改進PDI,從特征空間混合元素中去除植被反射率等信息的影響,減從而形成改進的垂直干旱指數MPDI(modified perpendicular drought index)。改進的垂直干旱指數計算公式:
式中 RV,RED為植被紅光波段反射率;RV,NIR為植被近紅光波段反射率;M 為土壤線斜率;fv 為植被覆蓋度。
(1)數據來源。針對灌區實際灌溉面積監測信息的需求,本文使用的遙感影像來自環境減災衛星HJ-1A/1BCCD,可由中國資源衛星應用網站下載,進而進行遙感指數反演,計算改進的垂直干旱指數。
(2)計算過程。對下載的遙感數據,經過去除熱噪聲、輻射定標[6]、相干斑濾波、地形校正預處理后得到水體和高植被區掩膜,大氣校正后,利用不同波段的反射率,選取目標區域,進行波段合成,結合實際實際灌溉面積的監測情況,構建Nir-Red 特征空間,利用線性回歸求出土壤線斜率M,根據M 和歸一化植被指數,分別計算區域灌溉前后的改進垂直干旱指數,在Erdasmodel 中構建MPDI 閾值提取區域實際灌溉面積監測模型:
式中I 為研究區域某個像元的垂直干旱指數差值;MPDIt1為灌溉前土壤垂直干旱指數;MPDIt2為灌溉后土壤垂直干旱指數。
當灌溉結束后,發生灌溉的區域土壤含水量增加,MPDI 減小;未發生灌溉的區域土壤含水量減小,MPDI 增大,可利用灌溉前MPDI 值減去灌溉后MPDI 值判斷土壤含水量變化情況,構建基于干旱指數差異閾值的灌溉面積監測模型。I 值越大,代表區域發生灌溉的可能性越大。如上所述,土壤濕度發生變化時,其光譜特征信息也會改變,因此理論上只要某一像元前后MPDI 存在差值,即I>0,該像元就有發生灌溉可能性。但土壤含水量的變化可能不是因為灌溉所造成的,也可能是由于風速等其他因素干擾,差值I 可能存在一定無效空間,因此需要設定一個閾值[7],只有I 變化程度大于該閾值時,才認為該像元的土壤水分變化是由灌溉造成的。
河北省某灌區始建于1959 年,設計灌溉面積5300 hm2,有效灌溉面積3000 hm2,經過多年發展,形成了1 總、4 干,16 支、92 斗渠的灌溉工程系統。灌區實行渠系專業三級組織管理:灌區管理處→渠道管理所→村農民用水戶協會。
根據現狀調查,在斗口以下的田間,該灌區未配套計量設施,目前按畝或按時收費,存在灌溉面積統計不準確、農戶認可度較低等缺點,一定程度上增加了水費征收的難度,不能滿足精確配水、合理灌溉的需要,管理部門調配用水和用水戶之間經常有矛盾發生,亟待采用客觀公正的方法,提取、識別灌溉面積。該灌區2020 年4 月2~11 日、5 月1~9 日,分別進行冬小麥拔節期、抽穗揚花期灌水。
采用灌區矢量數據進行裁剪,利用IDL(Interactive Data Language)語言建立MPDI 計算模型,依據公式分別計算灌溉前(2020 年5 月1 日)、后(2020 年5 月9 日)2 期影像的MPDI 差值,利用無人機航飛,結合GPS 點位信息,進行控制點布設,最終獲取灌溉區域的點位等信息,確定差值閾值0.17,對冬小麥灌溉面積進行提取,形成矢量圖及灌溉前后信息分布如圖2~圖5。

圖2 灌區矢量化圖形

圖3 5 月1 日(灌溉前)灌溉信息分布

圖5 灌溉面積提取分布
2.3.1 土壤含水量地面監測點驗證
兩期灌溉結束(4 月11 日、5 月9 日),在灌區選擇20 個土壤含水量實測點0~20 cm 深度平均土壤含水量的變化情況得出該時期各點是否灌溉。將提取的實際灌溉范圍與土壤含水量實測值點上灌溉情況相比,可在點尺度上進行驗證,如表1。兩輪灌溉驗證結果:4 月11 日灌溉范圍與點上實測結果的精度達0.9;5 月9 日灌溉范圍與點上實測結果的精度達0.85。

表1 提取結果點驗證
2.3.2 土壤灌溉面積驗證
4 月2 日,土壤干燥,冬小麥缺水,開始提閘灌溉;至4 月11 日,灌區完成冬小麥拔節期灌溉,對灌溉面積進行提取;5 月1 日陸續提閘放水,至5月9 日,完成抽穗揚花期灌水,再對灌溉面積進行提取。同時對灌區灌溉系統上報灌溉面積進行查詢,得到結果如表2。

表2 遙感提取灌溉面積與統計記錄灌溉面積比較
2.3.3 結果分析
實地調查發現,灌溉面積的統計都是由人工上報所得,其結果存在一些人為誤差,使統計灌溉面積會略小于真實灌溉面積,所以監測灌溉面積大于記錄灌溉面積是符合實際情況的。同時,在5月3 日,灌區部分區域有小雨,對提取的灌溉面積造成一定影響。但總體來看,遙感提取灌溉面積與統計記錄灌溉面積基本接近,可以作為統計灌溉面積的一項手段,為灌溉面積的計算節省大量的人力、物力、財力,提升灌區信息化管理水平,節省逐級上報的中間環節,提高管理運行效率。同時,可將遙感監測情況以一定方式,向用水單元發送,實現信息雙向流動和透明化。
(1)針對傳統灌溉面積上報、計算方式存在監測站點少、耗時費力、更新周期長等不足,充分發揮“3S”技術可快速采集與處理海量空間信息數據的優勢,以多源多時相遙感數據為基礎,以GIS 技術為支撐,進行灌溉地塊識別,基于灌溉前后改進的垂直干旱指數變化規律,結合地面少量調查點,實現實際灌溉時間和灌溉面積的監測,提出一種基于低空間分辨率遙感數據的低成本灌區灌溉面積提取的精確計量技術,為實現末級渠道分水計量提供技術支撐。
(2)結合現場低空攝像機拍攝實際灌溉信息,確定灌溉前后提取閾值,對研究區域灌溉面積進行提取,最終得到實際灌溉面積,用于校核上報的灌溉面積,計量末級渠道分水情況,提高農田中灌區面積提取方法的精度與效率。在灌區進行應用,其精度在0.85 左右,可在大中型灌區推廣應用。