基于高分辨率图像监督中分辨率图像的灌溉面积统计方法与流程

本发明涉及水资源管理领域，尤其涉及一种基于高分辨率图像监督中分辨率图像的灌溉面积统计方法。

背景技术：

基于遥感的灌溉面积统计是现代化农业管理不可或缺的管理手段，同时也是资源遥感的重要组成部分。我国土地辽阔，地形复杂，农户规模小，灌溉面积、农作物产量等信息依然是通过传统的全面报表和抽样调查两种方式获得的。这两种方式调查工作量大，费时费力且受人为因素影响较大，越来越难适应宏观决策和管理的需求。

随着遥感技术的发展，利用遥感手段进行大面积灌溉面积统计，进而进行蒸散发计算和作物水分利用效率评价的做法，具有一定的可行性和实用价值。不过，当需要对耕地区域内的作物进行识别时，目前大都需要以年为频率收集地表参考数据，即需要每年的地表实地采样数据才能完成对作物的识别。即便，近年来出现了使用一年数据监督多年分类结果的研究方法，但，当前年的训练数据仍然只有通过实地进行采集，仍然耗时耗力。

技术实现要素：

本发明的目的在于提供一种基于高分辨率图像监督中分辨率图像的灌溉面积统计方法，从而解决现有技术中存在的前述问题。

为了实现上述目的，本发明所述基于高分辨率图像监督中分辨率图像的灌溉面积统计方法，所述方法包括：

s1，确定研究区域，并获取所述研究区域的基于中分辨率成像光谱仪的中分辨率植被覆盖指数卫星影像包，记为modis-ndvi卫星影像包；

所述modis-ndvi卫星影像包中包括所述研究区域的一个训练年modis-ndvi卫星影像包和若干个实测年的modis-ndvi卫星影像包；

s2，获取训练年中若干日期的研究区域的地图卫星影像包；

s3，设定选取点的总数为q，其中，用地类型a的数量为α个，用地类型b的数量为β，用地类型c的数量为δ个，q＝α+β+δ；

提取训练年地图卫星影像包中各张地图卫星影像的中心点坐标，然后将任意一张地图卫星影像γ的中心点坐标及q选取点对应到训练年的modis-ndvi卫星影像包中与地图卫星影像γ的时间相同的modis-ndvi卫星影像γ中；

s4，首先，设定k个选取点作为率定点，q-k个选取点作为验证点；

然后，计算每张modis-ndvi卫星影像中的k个率定点的ndvi均值；

最后，形成训练年的ndvi均值曲线；

s5，采用五参数非对称logistic曲线对ndvi均值曲线进行最小二乘法函数拟合并得到五个参数的值；因，用地类型包括三种，故，每种用地类型分别对应一条logistic曲线，每种用地类型均具有五个参数；根据任意一种用地类型所对应的logistic曲线分别计算得到该土地类型的向量，用地类型a的向量为aa，用地类型b的向量为ab，用地类型c的向量为ac；

s6，采用五参数非对称logistic曲线计算任意一个验证点j的向量，记为avj，j＝1，2，......，q-k；

根据min(||avj-aa||2，||avj-ab||2，||avj-ac||2)获取验证点j的计算用地类型；

当最小值为||avj-aa||2；验证点j的用地类型为用地类型a；

当最小值为||avj-ab||2；验证点j的用地类型为用地类型b；

当最小值为||avj-ac||2；验证点j的用地类型为用地类型c；

s7，获取s4中验证点j的实测用地类型，判断验证点j的实测用地类型与验证点j的计算用地类型是否相同，如果是，则成功；如果否，则不成功；计算q-k个验证点的成功率，如果成功率大于等于90％，则进入s8；如果否，则成功率小于90％，则返回s5重新计算验证点的计算用地类型；

s8，对modis-ndvi卫星影像包中训练年的所有modis-ndvi卫星影像和若干实测年的所有modis-ndvi卫星影像进行s5至s6的处理，然后计算用地类型为耕地的灌溉面积。

优选地，步骤s1，具体按照下述步骤实现：

s11，确定研究区域，并获取研究区域的基于中分辨率成像光谱仪的中分辨率卫星图像包；所述中分辨率卫星图像包包括一个训练年的中分辨率卫星图像包和若干实测年中分辨率卫星图像包；所述中分辨率卫星图像包中的各个图像格式为.hdf；

s12，对所述中分辨率卫星图像包的各张中分辨率卫星图像进行utm重投影，得到各张中分辨率卫星图像的平面图像；

s13，对每张平面图像进行修剪，提取得到在所述研究区域经纬度范围内的修剪平面图像；

s14，预先设定两个不同的波段band1波段和band2波段，将修剪的平面图像的格式.hdf批量转换成两个波段的.gif格式图像；

s15，在同一张修剪平面图像的两个波段的.gif格式图像的基础上，采用公式(1)计算同一张修剪平面图像的ndvi值：

ndvi值＝(ρnir-ρred)/(ρnir-ρred)(1)

其中，ρnir为近红外波段反射率，单位％，即ρnir为band2波段的反射率；ρred为红色波段反射率，单位％，即ρred为band1波段的反射率；

将所有计算得到ndvi值的修剪平面图像打包成modis-ndvi卫星影像包。

更优选地，在步骤s13中的修剪过程中设置采样精度250m。

优选地，步骤s2中，训练年中相邻日期的两幅地图卫星影像的间隔时间为10-20d。

优选地，步骤s5中，根据任意一种用地类型所对应的logistic曲线分别计算得到该土地类型的向量，用地类型a的向量为aa，用地类型b的向量为ab，用地类型c的向量为ac；具体为：

获取每种用地类型对应的logistic曲线及其一阶导数的xinf1，yinf1，xmax-xinf1，ymax，xinf2-xmax，yinf2；

其中，xinf1表示曲线一阶导数最大点对应的时间；yinf1表示曲线一阶导数最大值点对应的曲线函数值；xmax-xinf1表示曲线最大值点和曲线一阶导数最大点之间的时间差，其中，xmax是曲线最大值对应的时间；ymax表示曲线最大值；xinf2-xmax表示曲线一阶导数最小值和曲线最大值之间的时间差，其中，xinf2表示曲线一阶导数的最小值；yinf2表示曲线一阶导数最小值点对应的曲线函数值；

然后通过每种用地类型对应的logistic曲线及其一阶导数的xinf1，yinf1，xmax-xinf1，ymax，xinf2-xmax，yinf2，计算得到对应用地类型的年内物候学特征向量，将用地类型a的向量记为aa，用地类型b的向量记为ab，用地类型c的向量记为ac。

优选地，所述用地类型包括灌溉用地、自然湿地及草地、荒地城镇。

本发明的有益效果是：

本发明完全通过高分辨率图像监督中分辨率图像的遥感手段，摆脱地面实地数据采集，实现监督分类并可以方便快捷的进行以年为频率的灌溉面积统计。

与现有技术相比，本发明利用modis卫星图像与googleearth进行了灌溉面积统计。本发明所述方法最大的特点是避免了使用实地参考数据，即不需要科研人员亲赴现场对不同种类作物进行比较研究，而只需要内业操作即可，且所有数据都可以免费得到，成本很低。故，可以称其为一种快速简洁的方法。

本发明所述方法的应用意义较大，方法不仅简单，同时，又由于是监督类算法，对不同地区的适应性能力强，可以作为统计灌溉面积的一个参考工具。根据实际算例分析，本发明所述方法结果较普查结果误差较小，具有较高的参考价值，在一定程度上代替普查结果或者从结果中分析偷垦等问题，为政策的制定提供一个简便快速的参考意见。

附图说明

图1是logistic曲线拟合ndvi曲线中的物候学特征示意图；

图2是新疆博斯腾流域2012年率定点三种用地年内ndvi均值曲线示意图；

1表示灌溉用地的ndvi均值曲线，2表示自然湿地及草地的ndvi均值曲线，3表示荒地城镇的ndvi均值曲线。

具体实施方式

为了使本发明的目的、技术方案及优点更加清楚明白，以下结合附图，对本发明进行进一步详细说明。应当理解，此处所描述的具体实施方式仅仅用以解释本发明，并不用于限定本发明。

modis表示中分辨率成像光谱仪；ndvi表示植被覆盖指数；utm的英文全称universaltransversemercatorprojection，中文简称通用横轴墨卡托投影；在utm投影的技术中，地球不同经度的地区有不同的utmzone；

更详细的解释说明为：

(一)步骤s1，免费获取打包卫星图像数据后分离并合成ndvi，然后再打包成modis-ndvi卫星影像包，具体为：

s11，确定研究区域，并获取研究区域的基于中分辨率成像光谱仪的中分辨率卫星图像包；所述中分辨率卫星图像包包括一个训练年的中分辨率卫星图像包和若干实测年中分辨率卫星图像包；所述中分辨率卫星图像包中的各个图像格式为.hdf；

从美国nasa官网查找对应研究区域的modis中分辨率卫星图像，下载得到研究区域的mod09gq系列图像数据，mod09gq系列图像数据的时间分辨率为日、格式为.hdf，下载时间跨度包括一个训练年与若干个实测年。

s12，利用modistool软件，对遥感数据进行utm重投影，utm重投影避免获取的modis中分辨率卫星图像数据因地球表面弯曲而产生失真的问题，通过utm重投影将modis中分辨率卫星图像精准地投影到一个平面，得到平面图像数据；查询研究区域的国际通用标准utmzone，并设置投影基准为wgs84。

s13，对平面图像数据进行修剪，提取在研究区域经纬度范围内的修剪平面图像，并保持修剪平面图像250m的采样精度，得到修剪平面图像数据包；

s14，预先设定.hdf格式的修剪平面图像数据包中需要的两个不同波段，band1及band2。利用批处理程序对修剪平面图像数据包中的数据进行格式转换处理，将.hdf格式的文件转化为.gif格式，该批处理程序集成于modistool中，反复调用modistool进行格式转化，将任意一张修剪平面图像转换成band1波段和band2波段的.gif00图像各一幅。利用matlab编写代码对图像的两个波段进行波段运算，计算ndvi值，最终得到研究区域的分辨率为250m×250m的modis-ndvi卫星遥感图像包，即modis-ndvi卫星影像包。ndvi的计算方法如下所示：

ndvi值＝(ρnir-ρred)/(ρnir-ρred)(1)

其中，ρnir为近红外波段反射率，单位％，即ρnir为band2波段的反射率；ρred为红色波段反射率，单位％，即ρred为band1波段的反射率。

utm投影过程中的基准，是为gps全球定位系统使用而建立的坐标系统。通过遍布世界的卫星观测站观测到的坐标建立。

(二)获取研究区域的训练年中若干日期的地图卫星影像，打包成地图卫星影像包。所述地图卫星影像包中只包括训练年中的若干日期的地图卫星影像，不需包括实测年；

根据当地的气候条件，训练年的选取会有不同，应选择可用图像多、关键日期附近云层覆盖少的年份进行训练，每幅图像的间隔时间大约在10-20d并尽量保持均匀，以增加后续插值过程中的准确性。

(三)任意选取一张训练年的地图卫星影像，从中选q个选取点，包括容易辨别的灌溉用地α个，自然湿地与草地β个，以及城市及荒地δ个。地图卫星影像包中的地图卫星影像仅仅是时间不同，分辨率和区域均相同，故，在任意一张地图卫星影像上设定q个选取点，则在其他地图卫星影像上均能找到一一匹配的q个选取点。

然后计算地图卫星影像包中每张地图卫星影像的中心点坐标，并将任意一张地图卫星影像γ的中心点坐标及其上的q选取点对应到训练年的modis-ndvi卫星影像包中与地图卫星影像γ的时间相同的modis-ndvi卫星影像γ中。

(四)用五参数非对称logistic曲线对ndvi均值曲线进行最小二乘法函数拟合并得到五个参数的值，灌溉用地，自然湿地与草地，城市及荒地分别对应一条logistic曲线，故，最后要生成3个参数组，每个参数组中包括五个参数。

本领域技术人员公知：logistic曲线的特点是开始增长缓慢，而在以后的某一范围内迅速增长，达到某限度后，增长又缓慢下来。logistic曲线略呈拉长的“s”型。而五参数logistic曲线的特点是不仅有曲线的上升段，也有曲线的下降段，所以适合于模拟作物在一个年周期中的生长成熟凋零过程。logistic曲线方程的另一个优点是其参数具有实际的物候学意义，这使得logistic曲线的应用与作物年变化规律十分契合。logistic曲线方程的表达形式如方程组(2)：

式中，x代表天数，即day(或称doy，dateofyear，指日期在一年中对应为第几天)，y代表函数的拟合值，共有a、b、c、d、e五个参数。附图1展示了一个logistic曲线。

对logistic曲线进行一阶求导，可以得到导函数，方程组(3)。

植被生长阶段中，有一些生长阶段的转折点。比如：生长最快的点，可以认为是ndvi增长最快的点，也就是一阶导数的最大值点；而生长成熟的点，可以认为是ndvi的最大值点。那么植物的生长时间可以定义为从生长最快的点到成熟点之间的时间，还可以得到从成熟到凋零最快的时期。这些值也是识别不同作物类型的关键数据。

另外，如前所述，logistic曲线的一个重要优点是其参数本身具有一定的生物学意义，以本申请中的非对称五参数logistic曲线为例，令logistic曲线的一阶导数为0，可以得到ndvi最大值为a+b，而ndvi取得最大值的时间即植被成熟的时间为c。另外，ndvi的最小值为a，即可以认为裸土的ndvi值为a；所以导出logistic曲线拟合的参数之后凭肉眼就可以轻松识别参数的合理性；参见附图1。

(五)通过得到的参数，计算训练年中用地类型对应的logistic曲线的最大值与最小值及对应时间点；计算logistic曲线一阶导数的最大值与最小值，及最大值与最小值分别对应的时间点。

获取每种用地类型对应的logistic曲线及其一阶导数的xinf1，yinf1，xmax-xinf1，ymax，xinf2-xmax，yinf2；

其中，xinf1表示曲线一阶导数最大点对应的时间；yinf1表示曲线一阶导数最大值点对应的曲线函数值；xmax-xinf1表示曲线最大值点和曲线一阶导数最大点之间的时间差，其中，xmax是曲线最大值对应的时间；ymax表示曲线最大值；xinf2-xmax表示曲线一阶导数最小值和曲线最大值之间的时间差，其中，xinf2表示曲线一阶导数的最小值；yinf2表示曲线一阶导数最小值点对应的曲线函数值；之所以使用(xmax-xinf1)与(xinf2-xmax)是为了排除年纪间水热变化带来的影响。

然后通过每种用地类型对应的logistic曲线及其一阶导数的xinf1，yinf1，xmax-xinf1，ymax，xinf2-xmax，yinf2，计算得到对应用地类型的年内物候学特征向量，将灌溉用地的向量记为aa，自然湿地及草地的向量记为ab，荒地城镇的向量记为ac。

根据五参数非对称logistic曲线，获取每个验证点的xinf1，yinf1，xmax-xinf1，ymax，xinf2-xmax，yinf2，然后计算得到任意一个验证点j的向量，记为avj(其中v代表verification，j代表是第j个验证点，j＝0，1，...，q-k。

根据min(||avj-aa||2，||avj-ab||2，||avj-ac||2)获取验证点j的计算用地类型；

当最小值为||avj-aa||2；验证点j的用地类型为用地类型a；

当最小值为||avj-ab||2；验证点j的用地类型为用地类型b；

当最小值为||avj-ac||2；验证点j的用地类型为用地类型c。

实施例1

本发明以新疆博斯腾流域为实施例，按照本申请所述方法进行流域灌溉面积统计，最终得到的识别结果。

博斯腾湖流域深居欧亚大陆腹地，远离海洋，呈现出明显的干旱大陆性气候特征。总的特点是干旱少雨，蒸发量大，多晴天，日照时间长，光热资源丰富。盆地平原区与海拔高程1700m以上的山区具有显著性差异，平原区四季分明，但冬长夏短，冬季长达110～115天，夏季56～69天。山区终年无夏，春秋相连，只有冷半年与暖半年之分，雨雪量远大于平原区，高寒风大。湖区气候受平原区大气候控制，但也有其独特性，由于大水体的储温效应，昼夜温差和最高最低温差较陆地小，局部风向风速也与陆地有所不同。博斯腾湖流域主要的气象站有巴音布鲁克、和静、和硕和焉耆气象站等。开都河是流入焉耆盆地的最大常年性河流，发源于天山中部艾尔宾、伊连哈尔嘎、那拉提、科克铁克等山脉，河源高度4292～4812m，河长560m，最终汇入博斯腾湖。

步骤一：使用2012年的卫星遥感数据进行监督训练

步骤二：根据以上分类模型，对2000至2012年中的灌溉面积进行统计，提取其中属于灌溉地的像元，计算数量，再乘以每个像元的面积得到总面积。

步骤三：完成地图区域识别后，在gis环境下完成后续统计工作。用给定区划的四个小灌区的掩膜进行裁剪，并最终合成为完整的灌区统计图，计算结果。

由于条件限制，只得到了2012年的当地普查统计值，这也是选取2012年为训练年的一个原因。对比可知，误差约为6.57％，经分析，误差的来源主要包含以下三个因素：

(1)模型的系统误差，由于与工作量妥协，率定点数量不足，可能导致将少量自然植被误识别为灌溉地；也可能既将部分灌溉地错误识别为自然植被，但同时，也有更多的自然植被被识别为灌溉地。

(2)当地普查统计过程中可能对灌溉面积有遗漏，比如城市中的绿地或其他植被。

(3)当地存在一定程度的偷垦现象，使得实际灌溉面积较普查值相对较高。

而通过分析从2000年到2012偶数年的计算结果，我们可以知道，在这12年间，博斯腾流域的灌溉面积呈稳步提高状态，其中，从2006到2008年有一个较大的增长。(五个偶数年的结果统计分别为303.82，344.35，365.39，385.39，484.12，505.84，564.03万亩)。

本发明通过高分辨率卫星遥感图像监督中分辨率卫星遥感图像的遥感监测手段，实现了监督分类并可以方便快捷的进行以年为频率的灌溉面积统计。

通过采用本发明公开的上述技术方案，得到了如下有益的效果：

本发明完全通过高分辨率图像监督中分辨率图像的遥感手段，摆脱地面实地数据采集，实现监督分类并可以方便快捷的进行以年为频率的灌溉面积统计。

与现有技术相比，本发明利用modis卫星图像与googleearth进行了灌溉面积统计。本发明所述方法最大的特点是避免了使用实地参考数据，即不需要科研人员亲赴现场对不同种类作物进行比较研究，而只需要内业操作即可，且所有数据都可以免费得到，成本很低。故，可以称其为一种快速简洁的方法。

本发明所述方法的应用意义较大，方法不仅简单，同时，又由于是监督类算法，对不同地区的适应性能力强，可以作为统计灌溉面积的一个参考工具。根据实际算例分析，本发明所述方法结果较普查结果误差较小，具有较高的参考价值，在一定程度上代替普查结果或者从结果中分析偷垦等问题，为政策的制定提供一个简便快速的参考意见。

以上所述仅是本发明的优选实施方式，应当指出，对于本技术领域的普通技术人员来说，在不脱离本发明原理的前提下，还可以做出若干改进和润饰，这些改进和润饰也应视本发明的保护范围。