遥感图像常用植被参数计算系统的设计与实现
2018-12-20黄义忠谭荣建马义超温煜未
钟 凯,黄义忠*,谭荣建,马义超,温煜未
遥感图像常用植被参数计算系统的设计与实现
钟 凯1,黄义忠1*,谭荣建1,马义超1,温煜未2
(1. 昆明理工大学 国土资源工程学院,云南 昆明 650093;2. 忻州市静乐县气象局,山西 忻州 035100)
植被遥感利用了植被对于各个波段的波谱曲线特征,广泛应用于农林作物的估算中,包括在一定范围农林作物种类的探测,种植生长区域面积的估算,生产和生态评价等。为了方便农业、林业及生态等非遥感专业领域研究人员快速获取植被参数,利用Visual Studio开发环境和ArcEngine开发工具建立了植被指数计算系统,通过归一化植被指数模型(NDVI)、比值植被指数模型(RVI),并结合矢量化的行政区图限定研究的范围,对该范围内的植被进行提取和对覆盖度进行分级。简化了遥感图像处理的繁琐操作过程,较为专业地进行计算、提取、分析,给出了具体的成果,能够在实际问题中为使用者提供更为便捷的植被参数计算方式。
遥感图像;植被参数;NDVI;RVI
0 引言
植被指数的引入基于绿色植物的叶片对于不同波段的电磁波的吸收与反射的特性,在波谱的红光波段范围与近红外波段范围内与其他地物的特性区别最大,由此得到植被指数的概念[1]。西方的发达国家于上世纪初就关注到环境监测方面的问题,发展较早,科学引用植被指数来治理环境,A Stepčenko[2]等人经过长期研究发现叶绿素的绿色质量、植物相对密度与植物的健康状态与归一化植被指数密不可分,说明了NDVI是植被遥感估算的重要指标。我国的李斌斌[3]等人利用分形布朗运动理论,结合归一化植被指数的空间分布,提出并计算了植被覆盖分形布朗运动分形维数。因此,可以看出,常用植被参数在林业、农业、土壤等生态方面具有重大研究意义[4-6],而本文的遥感图像常用植被参数计算系统正是基于快速获取常用植被参数的需求所开发设计,并经过实验证明可以快速获取所需的植被参数以及进行植被提取,能够给提取植被参数过程提供极大的便捷。
1 技术路线
系统功能的实现遵循如下图1所示的技术路线:
(1)对于需要计算植被参数的遥感原始数据进行加载,为了进行遥感数据的植被遥感估算,首先实现对典型植被指数的计算功能,典型植被指数可以反映出区域内植被信息,设置不同的植被指数的临界值提取出所需求的植被区域,最后对植被区域结果渲染显示以提供更为直观的信息。
(2)对加载的遥感影像利用数学模型进行处理,计算出植被覆盖度的栅格影像,对植被覆盖度成果图进行科学分类获取植被覆盖度的分级地图。
(3)按照植被遥感需求,对需求范围提取。
2 数学模型
2.1 常用植被指数
植被指数的计算原理大致可以概述为:植物叶片对电磁波谱的近红外波段与红光波段的吸收与反射具有一定的特性,利用遥感影像这两个波段的遥感栅格数据进行处理得到不同的植被指数[7]。利用典型的植被指数可以获取相应的植被信息,典型植被指数计算模型如下[8]:
比值植被指数(RVI),是近红外NIR波段与红光R波段段的比值,计算公式为[9]
RVI的取值范围大于0,在这个指数中有一个分界值2,在2以上代表了植被,且值越大植被的覆盖程度越高;在0-1之间一般是由于水体,冰雪的缘故,或者是遥感影像的自身原因比如云的影响;在1-2之间表示的是没有植被的区域,比如裸岩,建筑等等。
归一化植被指数(NDVI),计算公式为[10]
NDVI的取值范围在-1到1之间,分界值为0,在0以上表示的是我们需要的信息,0以下的是没有植物的区域;与比值植被指数相似,植被覆盖程度与指数呈正比。
差值植被指数(DVI),计算公式为
不同的植被指数有不同的局限性。比值植被指数对于大气的影像较为敏感,因此需要对得到的成果进行大气校;归一化植被指数在较高的植被覆盖区域内对于植被的识别度较弱,对于高植被覆盖区域的显示不如其他影像的效果;差值植被指数则受到土壤背景的影响较大。
通过科学实验表明,归一化植被指数(NDVI)尽管在高植被覆盖区域内敏感性不强,但这个指数对于土壤背景的变化相比其他指数具有优势,能够对于植被类型、长势、覆盖形状进行较好的反映。NDVI反映的植被信息较广,拥有良好的适应性,因此,利用NDVI来得到植被信息是本系统的重要标准。此外,为了弥补NDVI的部分局限性,系统也可以计算其他的植被指数,并利用指数相关的特性来获取植被信息。
2.2 像元二分模型
像元二分模型是一种简单实用的遥感估算模型,它假设一个像元的地表由有植被覆盖部分地表与无植被覆盖部分地表组成,而遥感传感器观测到的光谱信息也由这2个组分因子线性加权合成,各因子的权重是各自的面积在像元中所占的比率,如其中植被覆盖度可以看作是植被的权重[11-13]。李苗苗[14-15]等人利用NDVI并结合像元二分模型得到估算植被覆盖度的公式如下:
式中代表NDVI很低的区域,例如泥土、建筑等区域;代表NDVI很高有很多植物的区域。其具体计算公式如下:
因此,求得植被覆盖度的关键便是获取这两个值,如果近似的取覆盖度的最大值与覆盖度最小值为栅格影像的最高与最低,即和分别为100%和0。那么公式(4)可以化简为:
实际操作中不可避免有噪声存在,因此与的取值一般要选取一个置信度的范围内的值,置信度是依据不同的栅格影像来得到的。
在获取实际的最大与最小值可以借助ENVI来操作,获取最大与最小值按NDVI累积的百分比需要利用Compute Static工具来统计在栅格数据的集中方式。依据不同的遥感影像得到在累积概率较高或者较低的栅格影像点的范围,最后取得平均值便可以得到与。
2.3 重分类
重分类在地理数据分析中的使用规则为:基于地理数据如栅格数据中每个像素的原有像素值进行重新分配新值,对这些存在的像素值进行重新分配,并可以对原有的像素值用新的值来替代,新值与旧值的替换有一定的规则。本文经多次实验发现采用如图2所示的重分类规则,渲染显示的效果最佳,系统采用的重分类的原理如下。
3 功能模块
3.1 影像载入
系统处理的Landsat8遥感影像与Modis影像性质上为栅格数据,可以对于栅格数据加载,处理与渲染的方法操作应用于遥感影像,此外还需要添加openfiledialog类,openfiledialog拥有多个返回值,当点击OK后,可以根据这个返回值设立IF判断语句实现对文件的加载功能,此后为了打开栅格数据集必须有worksapce工作空间来打开,而工作空间的创建还需要的是workspacefactory工作空间工厂,而这两种类对于栅格数据还有专门的rasterworkspace与rasterworkspacefactory,还需强制转换。获取了栅格工作空间之后可以打开拥有路径的栅格数据。
以TIFF格式的Landsat8河南省郑州市遥感数据为例,利用openfiledialoge的返回值设置的if语句将打开文件的完整路径赋予定义字符串变量,设置整型数据标记“//”的位置,最后获取文件路径与文件路径,利用栅格工作空间打开数据集并且加载到图层之中。得到影像如图3所示。
该系统影像载入模块能够打开TIFF和HDF两种数据格式的文件,分别对应数据源Landsat8与Modis影像,但其载入原理类似且均在上面说明,故只以其中一种TIFF格式为例介绍具体代码:
OpenFileDialog openTiff = new OpenFileDialog();
openTiff.Filter = "Map Documents(*.TIF)|*. TIF";
openTiff.Title = "打开TIFF地图";
string Tiffpath;
if (openTiff.ShowDialog() == DialogResult.OK)
{
Tiffpath = openTiff.FileName;
int pos = Tiffpath.LastIndexOf("\");
string Filepath = Tiffpath.Substring (0, pos);
string Tiffname = Tiffpath.Substring (pos + 1);
IWorkspaceFactory workspacefac = new RasterWorkspaceFactory();
IRasterWorkspace rasterworkspc;
IRasterDataset rasterdatst = new RasterDataset();
IRasterLayer rasterlay = new RasterLayerClass();
rasterworkspc = workspacefac.OpenFromFile(Filepath, 0) as IRasterWorkspace;
rasterdatst = rasterworkspc.OpenRasterDataset(Tiffname);
rasterlay.CreateFromDataset(rasterdatst);
axMapControl1.AddLayer(rasterlay);
}
图3 影像载入
3.2 常用植被指数计算
植被指数计算功能涉及到IMathop类的方法调用,在自己定义的命名空间中定义了相关的函数进行计算,主要方法在于输入栅格数据集后对IMathop的方法引用,最后返回一个Geodatabase类的数据存储在相应的路径之中并加载至图层。以归一化植被指数NDVI计算为例,得到如图4、5所示的结果。
植被指数计算中,最关键的一步便是栅格数据相减,是在自定义的命名空间下(Function)编写的栅格数据计算的静态函数,首先if语句判断输入的参数是否为栅格数据进行计算的有效类型,接着实例化mathop对象并调用该对象的.minus()方法,最后返回一个类型为栅格数据集的对象。关键代码如下:
图4 计算过程
图5 NDVI成果图
public static IGeoDataset MathOpPlusRaster (IGeoDataset geoDataset_1, IGeoDataset geoDataset_2)
{
if ((geoDataset_1 is ESRI.ArcGIS.Geodatabase. IRaster | geoDataset_1 is ESRI.ArcGIS.Geodatabase. IRasterDataset | geoDataset_1 is ESRI.ArcGIS.DataSourcesRaster.IRasterBand | geoDataset_1 is ESRI. ArcGIS.GeoAnalyst.IRasterDescriptor) & (geoDataset_2 is ESRI.ArcGIS.Geodatabase.IRaster | geoDataset_2 is ESRI.ArcGIS.Geodatabase.IRasterDataset | geoDataset_2 is ESRI.ArcGIS.DataSourcesRaster. IRasterBand | geoDataset_2 is ESRI.ArcGIS.GeoAnalyst. IRasterDescriptor))
{
IMathOp mathOp = new ESRI.ArcGIS. SpatialAnalyst.RasterMathOpsClass();
IGeoDataset geodataset_output = mathOp. Plus(geoDataset_1, geoDataset_2);
return geodataset_output;
}
else
{
return null;
}
}
3.3 植被提取
植被的提取功能利用归一化植被指数与比值植被指数来确定。根据用户所需要的范围来区分植被与其他区域,点击信息提取栏的植被提取可以实现以上功能,首先是对用户需求的设置临界值的窗体设计,新建窗体如下图6所示。这里我们以NDVI值大于0.3为例进行提取,输入0.3点击确定,可以得到如图7所示的植被提取图结果。
设置临界值页面确定按钮的click事件中的植被提取代码,在计算栅格的过程中设置了遍历栅格像素的IRasterCursor类,这个类类似于指针,指向栅格图像上的每一个像素点,计算涉及的类还有IRasterEdit类,这个类的对象用以编辑栅格图像中的栅格属性值,以及IPixelBlock3数据块类,通过设置for循环语句判断栅格属性值并进行判断与临界值的关系设置新值Nan来实现提取依据用于的需求得到的栅格图像,最后调用myfunction的存储功能存储栅格数据集。关键步骤代码如下:
图6 设置临界值
图7 植被提取结果
IWorkspaceFactory OpenWorkspaceFac = new RasterWorkspaceFactory();
IRasterWorkspace OpenRasterWorkspc = OpenWorkspaceFac.OpenFromFile(Filepath, 0) as IrasterWorkspace;
IRasterDataset OpenRasterDataset = OpenRasterWorkspc.OpenRasterDataset(HDFname);
IRasterBandCollection OpenRasterBandCol = OpenRasterDataset as IRasterBandCollection;
IRasterBand OpenRasterBand = OpenRasterBandCol.Item(0);
OpenRasterBand.ComputeStatsAndHist();
IRasterStatistics pRasterStatistic = OpenRasterBand.Statistics;
Single dMaxValue = Convert.ToSingle(pRasterStatistic.Maximum);
Single dMinValue = Convert.ToSingle(pRasterStatistic.Minimum);
Single MaxMinusMinValue = Convert.ToSingle (dMaxValue - dMinValue);
IRaster OpenRater1 = OpenRasterDataset.CreateDefaultRaster();
IRaster2 OpenRaster = (IRaster2)OpenRater1;
IRasterDataset OutRasterDataset = OpenRasterDataset;
IRaster OutRaster1 = OutRasterDataset.CreateDefaultRaster();
IRaster2 OutRaster = (IRaster2)OutRaster1;
try
{
IRasterCursor rasterCursor = OpenRaster.CreateCursorEx(null);
IRasterEdit rasterEdit = OutRaster as IRasterEdit;
IPixelBlock3 pixelblock3 = null;
int blockwidth = 0;
int blockheight = 0;
System.Array pixels;
IPnt tlc = null;
Single pixval;
do
{
pixelblock3 = rasterCursor.PixelBlock as IPixelBlock3;
blockwidth = pixelblock3.Width;
blockheight = pixelblock3.Height;
pixels = (System.Array)pixelblock3. get_PixelData(0);
for (long i = 0; i < blockwidth; i++)
{
for (long j = 0; j < blockheight; j++)
{
if (Convert.ToSingle(pixels.GetValue(i, j)) == Single.MinValue || Convert.ToSingle(pixels.GetValue(i, j)) == Single.MaxValue)
continue;
else
{
pixval = (Convert.ToSingle(pixels.GetValue(i, j)) - dMinValue) / MaxMinusMinValue;
pixels.SetValue(pixval, i, j);
}
}
}
pixelblock3.set_PixelData(0, pixels);
tlc = rasterCursor.TopLeft;
rasterEdit.Write(tlc, (IPixelBlock) pixelblock3);
}
while (rasterCursor.Next() == true);
rasterEdit.Refresh();
System.Runtime.InteropServices.Marshal.ReleaseComObject(rasterEdit);
3.4 重分类渲染
重分类渲染的主要作用是将提取图渲染为直观的表现图,但提取的功能在于表示数据的属性,能够将大于临界值的所有栅格数据提取出来并且保留它们的像元。以提取二分为例,通过调用GP工具来实现,除了输入输出栅格图层的参数外,设置重分类标准与覆盖度的分类标准不一致,只需要设置两个分类0与1即可将植被与裸土提取出来。得到的提取二分结果如图8所示。
图8 重分类渲染结果
代码设计为先调用GP工具并实例化重分类工具对象,设置好输入、输出、分类标准与字段名等参数,最后执行Excute。在其他窗体中的事件中只需要将这些参数设置好,直接调用命名空间中的这个函数即可。调用GP工具的具体代码如下:
Geoprocessor tGp = new Geoprocessor();
Reclassify RasterRecla = new Reclassify();
RasterRecla.in_raster = inRaster;
RasterRecla.out_raster = outRaster;
RasterRecla.reclass_field = "Value";
RasterRecla.remap = "-1 0.3 0;0.3 1 1;";
IGeoProcessorResult tGeoResult = (IGeoProcessorResult)tGp.Execute(RasterRecla, null);
return outRaster;
3.5 裁剪
系统的成果裁剪所采用的方法是调用GP工具,GP工具的使用的参数为输入的栅格数据,输入的掩膜数据以及输出的路径,首先新建一个GP工具,创建新的临时gbg文件,在运算结束后会自动删除,设置临时的工作环境。最后输入要裁剪的数据与矢量数据。得到的裁剪后的结果如图9所示。
调用GP工具实例化Extractbymask工具的关键代码如下:
Geoprocessor extractGeoprocessor = new Geoprocessor();
string gdb = Myfunction.CreateFileGDB(inRaster + "random");
extractGeoprocessor.SetEnvironmentValue("workspace", gdb);
extractGeoprocessor.OverwriteOutput = true;
if (File.Exists(extractByMask))
File.Delete(extractByMask);
ExtractByMask newExtractByMask = new ExtractByMask()
{
in_raster = inRaster
out_raster = extractByMask
in_mask_data = inMask
}
图9 郑州市及周边植被提取结果
此时,我们可以将郑州行政区划图的图层添加至图9的植被提取图,得到如下图10所示的郑州市植被比较分析图。据此,可以做出相关的植被资源评价:郑州市的植被覆盖最少,需要重视植被树林的栽种;新郑、中牟、新密的植被覆盖较高,需要重视保护森林并预防森林火灾;登封、巩义、荥阳覆盖面积居中,需要坚持植被森林的浇灌和维护。
4 结语
本文开发设计出的遥感图像常用植被参数计算系统,经过实验证明,能够载入Landsat8与Modis遥感影像产品,通过常用植被指数计算功能快速计算出RVI、NDVI和DVI的数值,并通过植被提取和重分类渲染的功能将植被区域筛选出来,结合郑州行政区划图得到了郑州市植被比较分析图,得出了郑州市的植被覆盖最少,新郑、中牟、新密的植被覆盖较高,登封、巩义、荥阳覆盖面积居中的植被覆盖状况,能够进行相关的植被资源评价并提供一些决策依据。
[1] 梅安新. 遥感导论[M].北京: 高等教育出版社, 2001.
[2] Stepčenko A, Čižovs, Jurijs. Markov Chain Modelling for Short-Term NDVI Time Series Forecasting[J]. Information Technology & Management Science, 2016, 19(1).
[3] 李斌斌, 李占斌, 宇涛, 李鹏. 基于归一化植被指数的流域植被覆盖分形维数研究[J]. 农业工程学报, 2014, 30(15): 239-247.
[4] 徐斌, 杨秀春, 金云翔著. 草原植被遥感监测[M]. 北京: 科学出版社, 2016.
[5] 谭昌伟. 农业遥感技术[M]. 北京: 中国农业出版社, 2017.
[6] 林辉. 林业遥感[M]. 北京: 中国林业出版社, 2011.
[7] 孙家抦. 遥感原理与应用[M]. 武汉: 武汉大学出版社, 2013.
[8] 张学艺. 基于EOS/MODIS的几种植被指数[J]. 安徽农业科学, 2009, 37(26): 12842-12845.
[9] Jordan C F. Derivation of Leaf‐Area Index from Quality of Light on the Forest Floor[J]. Ecology, 1969, 50(4): 663-666.
[10] Deering D W. Rangeland reflectance characteristics measured by aircraft and spacecraft sensors[J]. Ph.d.thesis Texas A & M Univ.college Station, 1978, 338.
[11] Catherine Leprieur, Michel M. Verstraete, Bernard Pinty. Evaluation of the performance of various vegetation indices to retrieve vegetation cover from AVHRR data[J]. Remote Sensing Reviews, 1995, 10(4): 265-284.
[12] Zribi M, LE HE’GARAT-MASCLES, TaconetO.D, BoussemaMR. et al. Derivation of wild vegetation cover density in semi-arid regions: ERS2 SAR evaluation[J] .International Journal of Remote Sensing , 2003 ,(24):1335-1352.
[13] 陈晋, 陈云浩, 何春阳, 等. 基于土地覆盖分类的植被覆盖率估算亚像元模型与应用[J]. 遥感学报, 2001, 5(6): 416-422.
[14] 李苗苗. 植被覆盖度的遥感估算方法研究[D]. 中国科学院研究生院(遥感应用研究所), 2003.
[15] 李苗苗, 吴炳方, 颜长珍, 等. 密云水库上游植被覆盖度的遥感估算[J]. 资源科学, 2004, 26(4): 153-159.
Design and Achieve of Common Vegetation Parameter Calculation System for Remote Sensing Images
ZHONG Kai1, HUANG Yi-zhong1*, TAN Rong-jian1, MA Yi-chao1, WEN Yu-wei2
(1. Faculty of Land Resource Engineering, Kunming University of Science and Technology, Kunming 650093, China; 2. Jingle County Meteorological Bureau of Xinzhou City, Xinzhou 035100, China)
Vegetation remote sensing utilizes the spectral characteristics of vegetation for each band, and is widely used in the estimation of agricultural and forestry crops, including the detection of a range of agricultural and forestry crops, the estimation of planting and growing area, production and ecological evaluation. In order to facilitate the rapid acquisition of vegetation parameters by researchers in non-remote sensing fields such as agriculture, forestry and ecology, a vegetation index calculation system was established using Visual Studio development environment and ArcEngine development tools, through the normalized difference vegetation index (NDVI) model or the ratio of vegetation index (RVI) model, combined with the vector of the administrative region defining the scope of the study, extract the vegetation within the range and grade the coverage. It simplifies the cumbersome operation process of remote sensing image processing, performs calculation, extraction and analysis more professionally, and gives specific results. It can provide users with more convenient calculation methods of vegetation parameters in practical problems.
Remote sensing image; Vegetation parameters; NDVI; RVI
TP751
A
10.3969/j.issn.1003-6970.2018.11.003
国家自然基金地区基金,滇中城市群国土空间格局多尺度演化模拟及优化配置(KKGD201721098)
钟凯(1995-),男,昆明理工大学硕士研究生,研究方向为国土资源遥感;谭荣建(1965-),男,昆明理工大学副教授、硕士生导师,研究方向为测绘与土地资源管理;马义超(1993-),男,昆明理工大学硕士研究生,研究方向为地理信息工程;温煜未(1994-),男,华北水利水电大学毕业,山西省忻州市静乐县气象局从事工作,研究方向为地理信息科学。
黄义忠(1972-),男,昆明理工大学副教授、硕士生导师,研究方向为土地资源管理。
钟凯,黄义忠,谭荣建,等. 遥感图像常用植被参数计算系统的设计与实现[J]. 软件,2018,39(11):11-17