《ERDAS分类后处理与ArcGIS数据交换.doc》由会员分享，可在线阅读，更多相关《ERDAS分类后处理与ArcGIS数据交换.doc（11页珍藏版）》请在taowenge.com淘文阁网|工程机械CAD图纸|机械工程制图|CAD装配图下载|SolidWorks_CaTia_CAD_UG_PROE_设计图分享下载上搜索。
1、【精品文档】如有侵权，请联系网站删除，仅供学习与交流ERDAS分类后处理与ArcGIS数据交换.精品文档.基于专家知识的决策树分类可以将多源数据用于影像分类当中，这就是专家知识的决策树分类器，本专题以ENVI中Decision Tree为例来叙述这一分类器。概述基于知识的决策树分类是基于遥感影像数据及其他空间数据，通过专家经验总结、简单的数学统计和归纳方法等，获得分类规则并进行遥感分类。分类规则易于理解，分类过程也符合人的认知过程，最大的特点是利用的多源数据。如图1所示，影像+DEM就能区分缓坡和陡坡的植被信息，如果添加其他数据，如区域图、道路图土地利用图等，就能进一步划分出那些是自然生长的植
2、被，那些是公园植被。图1 专家知识决策树分类器说明图专家知识决策树分类的步骤大体上可分为四步：知识（规则）定义、规则输入、决策树运行和分类后处理。1.知识（规则）定义规则的定义是讲知识用数学语言表达的过程，可以通过一些算法获取，也可以通过经验总结获得。2.规则输入将分类规则录入分类器中，不同的平台有着不同规则录入界面。3.决策树运行运行分类器或者是算法程序。4.分类后处理这步骤与监督/非监督分类的分类后处理类似。知识（规则）定义分类规则获取的途径比较灵活，如从经验中获得，坡度小于20度，就认为是缓坡，等等。也可以从样本中利用算法来获取，这里要讲述的就是C4.5算法。利用C4.5算法获取规则可分
3、为以下几个步骤：（1）多元文件的的构建：遥感数据经过几何校正、辐射校正处理后，进行波段运算，得到一些植被指数，连同影像一起输入空间数据库；其他空间数据经过矢量化、格式转换、地理配准，组成一个或多个多波段文件。（2）提取样本，构建样本库：在遥感图像处理软件或者GIS软件支持下，选取合适的图层，采用计算机自动选点、人工解译影像选点等方法采集样本。（3）分类规则挖掘与评价：在样本库的基础上采用适当的数据挖掘方法挖掘分类规则，后基于评价样本集对分类规则进行评价，并对分类规则做出适当的调整和筛选。这里就是C4.5算法。4.5算法的基本思路,如下：从树的根节点处的所有训练样本D0开始，离散化连续条件属性。
4、计算增益比率，取GainRatio（C0）的最大值作为划分点V0，将样本分为两个部分D11和D12。对属性C0的每一个值产生一个分支，分支属性值的相应样本子集被移到新生成的子节点上，如果得到的样本都属于同一个类，那么直接得到叶子结点。相应地将此方法应用于每个子节点上，直到节点的所有样本都分区到某个类中。到达决策树的叶节点的每条路径表示一条分类规则，利用叶列表及指向父结点的指针就可以生成规则表。图2 规则挖掘基本思路算法描述如下：算法：从空间数据集（多波段文件）中挖掘分类规则输入：训练样本输出：分类规则表方法：一、读取数据集名字二、读取所有的训练样本A、读取属性信息C、原始类E、样本值A，并将样
5、本划分为训练样本（2/3）和评价样本（1/3）。B、属性信息C可以是连续（DISCRETE）或离散（CONTINUOUS）的，分别将属性注上这两种标记；若属性是DISCERTE，读取其可能取得值，并都存储在一个列表中；每一个属性都有一个标记，一个给定的属性编号及初始化的取值列表均存储于一个属性的数据结构中，并将数据结构存储在一个哈希表中。C、原始类E当作一个附加属性信息储存在属性列表中。D、以增量方式读取每一个样本A，将所有的样本储存在一个表中，每一行代表一个样本。三、利用数据集构建树A、离散化连续条件属性C DISCRETE，获得的分割点集T（t1，t2）作为条件属性C的新的取值。B、分别计
6、算所有条件属性的增益比率GainRatio（C），取增益比率值最大的条件属性作为树的划分节点，其值或范围作为划分值V（v1，v2）来生成树的分枝。C、判断该层与每一个等价子集的原始类类别是否一致。若一致，生成叶子结点。否则，继续计算增益比率GainRatio（C）和选择条件属性C，得到树的节点和划分值V，直至所有的样本已分类完毕。四、测试生成树将测试样本C带入树中，当某一测试样本的分类预测错误时，记录分类错误的计数，并将测试样本添加到训练样本中，转向步骤三，重新构建树。否则，输出分类树五、抽取分类规则到达树的叶节点的每条路径表示一条分类规则从树中抽取分类规则，打印规则和分类的详细信息C4.5网
7、上有源代码下载，vc和c+版本都能获得。Decision Tree的使用一、规则获取选取Landsat TM5影像和这个地区对应的DEM数据，影像和DEM经过了精确配准。规则如下描述：Class1（朝北缓坡植被）:NDVI0.3, slope20, aspect270Class2（非朝北缓坡植被）：NDVI0.3, slope20, 90=aspect0.3, slope=20,Class4（水体）：NDVI=0.3, 0b420Class5（裸地）：NDVI=20Class6（无数据区，背景）: NDVIclassification-Decision Tree-Build New Decis
8、ion Tree,如图3所示，默认显示了一个节点。图3 Decision Tree界面首先我们按照NDVI的大小划分第一个节点，单击Node1，跳出图4对话框，Name为NDVI0.3，在Expression中填写：ndvi gt 0.3。图4 添加规则表达式点击OK后，会提示你给ndvi指定一个数据源，如图5所示，点击第一列中的变量，在对话框中选择相应的数据源，这样就完成第一层节点规则输入。图5 指定数据源Expression中的表达式是有变量和运算符（包括数学函数）组成，支持的运算符如表1所示表达式部分可用函数基本运算符 +、-、*、/ 三角函数 正弦Sin(x)、余弦cos(x)、正切t
9、an(x)反正弦Asin(x)、反余弦acos(x)、反正切atan(x) 双曲线正弦Sinh(x)、双曲线余弦cosh(x)、双曲线正切tanh(x) 关系/逻辑 小于LT、小于等于LE、等于EQ、不等于NE、大于等于GE、大于GTand、or、not、XOR最大值（）、最小值 ()其他符号 指数（）、自然指数exp自然对数对数alog(x) 以10为底的对数alog10(x)整形取整round(x)、ceil(x)平方根（sqrt）、绝对值（adb）表1 运算符ENVI决策树分类器中的变量是指一个波段的数据或作用于数据的一个特定函数。变量名必须包含在大括号中，即变量名；或者命名为bx，x代
10、表数据，比如哪一个波段。如果变量被赋值为多波段文件，变量名必须包含一个写在方括号中的下标，表示波段数，比如pc2表示主成分分析的第一主成分。支持特定变量名如表2，也可以通过IDL自行编写函数。变量作用slope 计算坡度 aspect 计算坡向（南坡北坡或者东西方向）ndvi 计算归一化植被指数 Tascap n穗帽变换，n表示获取的是哪一分量。 pc n主成分分析，n表示获取的是哪一分量。lpc n局部主成分分析，n表示获取的是哪一分量。 mnf n 最小噪声变换，n表示获取的是哪一分量。Lmnfn局部最小噪声变换，n表示获取的是哪一分量。Stdev n波段n的标准差 lStdev n波段n
11、的局部标准差Mean n波段n的平均值 lMean n波段n的局部平均值Min n、max n波段n的最大、最小值 lMin n、lmax n波段n的局部最大、最小值表2变量表达式第一层节点根据NDVI的值划分为植被和非植被，如果不需要进一步分类的话，这个影像就会被分成两类：class0和class1。对NDVI大于0.3，也就是class1，根据坡度划分成缓坡植被和陡坡植被。在class1图标上右键，选择Add Children。单击节点标识符，打开节点属性窗口，Name为SlopeExecute，执行决策树，跳出图7所示对话框，选择输出结果的投影参数、重采样方法、空间裁剪范围（如需要）、输
12、出路径，点击OK之后，得到如图8所示结果。在决策树运行过程中，会以不同颜色标示运行的过程。图7 输出结果图8 决策树运行结果回到决策树窗口，在工作空白处点击右键，选择Zoom In，可以看到每一个节点或者类别有相应的统计结果（以像素和百分比表示）。如果修改了某一节点或者类别的属性，可以左键单击节点或者末端类别图标，选择Execute，重新运行你修改部分的决策树。图9 运行决策树后的效果分类后处理和其他计算机分类类似的过程。ERDAS分类后处理与ArcGIS数据交换分类后，得到一张专题图。由于ERDAS是基于像素的分类，在分类结果中会出现很多细小的杂点以及小的图斑，如何去掉这些杂点和小图斑（称为
13、ERDAS后处理），以及如何将专题图转换成ArcGIS下的shp格式，是本文需要解决的工作。一、ERDAS后处理1、首先给出原始图像和经过监督分类得到的分类图像，如下图：原始影像监督分类结果影像2、开始聚类使用Interpreter模块 - GIS Analysis - Clump，打开聚类界面，输入影像为分类后影像，给一个输出影像，聚类邻域可以选择8或者4，点击OK，即可进行聚类，结果如下：聚类结果3、去除杂点和小图斑使用Interpreter模块 - GIS Analysis - Eliminate，也可以使用Sieve，后者容易出现一些属性值为0的像素，不太好处理，所以这里使用Elimi
14、nate，输入影像为上一步聚类后的影像，给定一个输出影像，Minimum处设置杂点和小图斑的像素大小，这里可以根据实际情况设定，这里选择100，点击OK，结果如下。通过比对，可以发现分类结果中很多杂点和小图斑被有效地去除了。去除小图斑去除杂点和小图斑 4、形态学处理也可以使用形态学的腐蚀、膨胀、开运算、闭运算来消除杂点，毛刺，空洞等。使用Interpreter模块 - Utilities - Morphological，选择输入输出影像，结构元素kernel definition可以选择3*3，5*5或7*7，也可以自定义，function中Erode是腐蚀，Dilate是膨胀，Open是开运
15、算，Close是闭运算，选择Open，然后点击OK即可，只是得到的结果不是专题图，需要重新分类得到专题图。不过这个时候的分类，将会非常简单，不再赘述。形态学开运算 二、ArcGIS数据交换将分类结果转换成ArcGIS的矢量格式，即可在ArcGIS上进行后续的处理。使用Vector模块 - Raster to Vector，将分类结果转换成ArcInfo的Arc Coverage格式。该格式可以在ArcMap中直接打开，如下图：ArcMap中打开Arc Coverage格式也可以使用ArcToolBox中的Conversion Tools - To Shapefile - Feature Cla
16、ss to Shapefile(multiple)，将Arc Coverage格式转换成shp格式。ArcToolBoxArcMap中打开的shp格式，图中黑点是注记层剩下的事情，就是通过ArcGIS提供的强大的空间分析功能，提取各种地物的各种信息了。ARCGIS转移矩阵土地利用转移矩阵生成的几种方法查阅相关的资料，也没有得到土地利用类型转换矩阵确切的定义，我理解为不同时间段内同一区域内土地利用类型的相互转换关系，一般用二维表来表达，从二维表中可以快速查看各个地类间相互转化的具体情况。比如某一类别的土地有百分之多少（或者面积）分别转化成了其他的土地类型，现在某类型的土地分别是由过去的哪些类别转
17、化而来的等等。还可以生成变化统计栅格图（掩膜图像），它描述了前后两幅土地分类图之间的地类发生转变的位置和类别。土地利用类型转换矩阵可以从两幅栅格图中计算得到，也可以从两个矢量文件中计算获得。下面介绍在ENVI下从两幅分类结果的栅格图中计算土地利用类型转换矩阵。1、准备数据两个时相的土地利用分类结果，它是单波段、专题类型的伪彩色图像（ENVI Classification）。2、计算转换矩阵打开两个土地利用分类结果。（1） 在主菜单中，选择Basic Tools Change Detection Change Detection Statistics。（2）分别在Initial State对话框
18、和final state对话框中选择前一时相和后一时相的土地利用结果。（3） 在Define Equivalent Classes对话框中（图1），如果两个土地利用分类名称一致，系统自动将Initial State Class和Final State Class对应，否则手动选择，单击Add Pair。（4） 选择对应的地物类型之后，单击OK按钮，出现图2对话框。选择生成图表表示单位（Report Type）：像素（Pixels）、百分比（Percent）和面积（Area）。选择Output Classification Mask Images?为YES，输出掩膜图像，选择输入路径及文件名。（
19、5） 单击OK，执行土地利用类型转换矩阵计算过程。图1 Define Equivalent Classes对话框图2 选择数据参数3、查看结果（1）如图3为得到的土地利用类型转换矩阵结果。横字段表示前一时间段（Initial State）的土地利用类别，纵字段为后一时间段（Final State）的土地利用类别。横字段和纵字段交叉处表示变化值，如有2520900平方米林地用地变化为草地。图3 土地利用类型转换矩阵（2）还可以为每一个地类生成一个变换掩膜图像，图4所示为其中一个地类的掩膜图像。掩膜图像的灰度值表示变化类型，如这里的2草地表示林地变化为草地的像元。图4 变化掩膜图像根据你的数据类型
20、选用不同的数据生成方法ERDAS中vector数据转移矩阵方法数据是Vector格式(shp数据格式)1 Erdas Imagine-Interpreter-Gis Analysis-Matrix，输入两个时相的Vector数据即可此时注意：输出栅格大小不应设的太小 要不一运行就会提示你的空间不足做这一步之前，请做好前期的地理编码。2 ArcView3.3加载 spatial analysis模块把两时相的Vector图转成grid格式（当然中间有一些单位的设置根据你做的图的分辨率来设置即可）analysis-map calculate 直接计算即可。3 把两期解译完的Vector文件在arc
21、toolbox-overlay-union中叠加，注意：两个文件不能用同一个字段名，比如一个用93Type，另一个时相则用 00Type叠加后的文件在Arcmap中打开，选中文件，然后点右键Property空间查询，输入条件语句，比如：93Type=1And 00Type=2；查询结果即为第一种类型转化为第二种类型的图形，可以另建一图层比如：12，把查询结果复制到12图层上。统计出面积，依进行，就可以得到土地利用类型转移矩阵 。一、数据准备（图1）准备两幅不同时相的土地利用现状图（shp格式），每幅图的属性表都要有一个表示土地利用类型的字段，并且要使用不同的名称加以区分，如Type1995，T
22、ype2000。土地利用类型名称必须统一，并且完整，如都使用“城镇用地”、“有林地”等。二、数据融合（图2）DISSOLVE在ArcMap里分别打开两个时相的图层，打开ArcToolbox，选择Data Management Tools | Generalization | Dissolve工具。Input Feature选择要融合的图层，Output Feature Class选择输出结果存储的位置及名称，Dissolve Field(s)选择土地利用类型字段（如Type1995），然后勾选Creat multipart features选项，点击OK完成。重复此过程，对另一时相数据进行融合
23、。此步骤使相同利用类型的记录融合为一个记录，以提高后面步骤的计算速度。三、叠置分析（图3）OVERLAY在ArcMap中打开两个时相融合后的数据，在ArcToolbox中选择Analysis Tools | Overlay | Intersect工具，Input Features选择两个时相的图层，Output Feature Class选择叠加结果存储的位置及名称，其余选项可以忽略，单击【OK】完成。四、计算面积并导出属性表（图4-6）在ArcMap中打开叠加后的图层数据，在该图层上右键打开属性表，选择Option | Add field 新建一个字段，命名为NewArea。在Editer工
24、具条中选择Editer | Start Editing，然后在属性表中NewArea字段上单击右键选择Calculate Geometry ，在打开的Calculate Geometry对话框中，Property选择Area，Units选择要使用的面积单位，单击【OK】完成图斑面积计算。依次选择Editer | Save Edits / End Editing保存和退出编辑状态。在属性表中选择Option | Export 将属性表保存为dbf文件。五、制作转移矩阵（图7-10）（以Excel2007为例）在Excel中打开上一步保存的dbf，另存为Excel格式并打开。在Excel中选中所有
25、数据（不要点左上角，只选择有效数据），点击【插入】选项卡，选择【数据透视表】|【数据透视表】，点击【确定】。在打开的数据透视表中按图示将字段拖入相应区域。Excel自动计算矩阵，将该表稍事整饰就得到美观的土地利用转移矩阵。矩阵中r(I, j)就表示i类型向j类型转移的土地面积，空值表示i类型向j类型没有转移。Fulu土地利用转移矩阵生成的几种方法根据你的数据类型选用不同的数据生成方法若你的数据是Raster格式：则有如下方法1 Erdas Imagine-Interpreter-GIS Analysis-Matrix，输入两个时相的Raster数据即可做这一步之前记得：先对两时相的数据进行重编
26、码（Interpreter-GIS Analysis-Recode）一般运行如果出现错误 肯定是重编码没做好，请继续查证。2 先在 Erdas中利用 Modeler 计算 如下公式NC(I,J)=NC(I)*10+NC(J),(JI)其中：NC(I,J)表示i，j 两年份的土地利用变化图；NC（i）表示i年份遥感分类影像；NC（j）表示j年份的遥感分类影像。在此计算的基础上，将以上变化影像图转化为BIL格式，再利用ARC/INFO GRID模块将影像转为GRID格式，然后利用GRID模块中的属性表（vat）查看命令对影像灰度值进行统计，最后得出土地利用转化举证。（注:此方法本人尚未实现过，不知
27、可行否）。若数据是Vector格式1 Erdas Imagine-Interpreter-Gis Analysis-Matrix，输入两个时相的Vector数据即可，此时注意：输出栅格大小不应设的太小要不一运行就会提示你的空间不足。做这一步之前，请做好前期的地理编码。2 ArcView3.3加载 spatial analysis模块把两时相的Vector图转成grid格式（当然中间有一些单位的设置根据你做的图的分辨率来设置即可）analysis-map calculate 直接计算即可。3 把两期解译完的Vector文件在arctoolboxoverlayunion中叠加，注意：两个文件不能用
28、同一个字段名，比如一个用93Type，另一个时相则用 00Type叠加后的文件在Arcmap中打开，选中文件，然后点右键Property空间查询，输入条件语句，比如：93Type=1And 00Type=2；查询结果即为第一种类型转化为第二种类型的图形，可以另建一图层比如：12，把查询结果复制到12图层上。统计出面积，依进行，就可以得到土地利用类型转移矩阵 。最后输出土地利用变化图 ，如下图所示：集美大学GIS（）FLAASH大气校正Flaash大气校正（IRSP6-08.3.24）大气校正的目的是消除大气和光照等因素对地物反射的影响，获得地物反射率和辐射率、地表温度等真实物理模型参数，用来消
29、除大气中水蒸气、氧气、二氧化碳、甲烷和臭氧对地物反射的影响，消除大气分子和气溶胶散射的影响。FLAASH 可以处理任何高光谱数据、卫星数据和航空数据(860nm/1135nm)，这些数据是由HyMAP、AVIRIS、CASI、 HYDICE、HYPERION(EO-1)AISA、HARP、DAIS、Probe-1、TRWIS-3、SINDRI、MIVIS、 OrbView-4、NEMO等传感器获得的。FLAASH还可以校正垂直成像数据和侧视成像数据。Flaash大气校正使用了 MODTRAN 4+ 辐射传输模型的代码，基于像素级的校正，校正由于漫反射引起的连带效应，包含卷云和不透明云层的分类图
30、，可调整由于人为抑止而导致的波谱平滑。FLAASH可对Landsat, SPOT, AVHRR, ASTER, MODIS, MERIS, AATSR, IRS等多光谱、高光谱数据、航空影像及自定义格式的高光谱影像进行快速大气校正分析。能有效消除大气和光照等因素对地物反射的影响，获得地物较为准确的反射率和辐射率、地表温度等真实物理模型参数。校正过程点击ENVIBasic Tools Preprocessing Calibration Utilities FLAASH Spectral FLAASH.或者点击ENVI-spectral- FLAASH1、 输入数据必须是辐射定标后的数据（表观反射
31、率），对辐射校正数据转成BIL或BIP格式（Basic Tools -Convert Data）；2、 对输入数据进行头文件编辑，主要是对波长wavelenth（即每一波段的波长中心值）和波长宽度fwhm（每一波段的波长范围）的编辑。不是高光谱数据可以不对fwhm进行编辑。（ENVIfileEdit Envi Header）3、 输入数据后，弹出如下对话框共有两种选择，如果输入影像不同波段有不同的转换因子，那选择第一种，反之第二种。我用的是irs影像所有波段都为同一因子，所以选用第二种，因子的值根据输入数据的单位与ENVI标准单位的转换尺度：Radiance Scale Factors是一个单
32、位转换因子，如果你的radiance（光谱灵敏度）是标准单位w/m2 *um *rad ，而flaash要求输入的是uw/cm2*sr*nm，则该因子为10。1m=103mm=106m=109nm=1012pm(皮米)1w=103mw=106w 1兆瓦=106瓦Rad平面角弧度 sr 立体角球面度4、 设置输出参数，包括：Output Reflectance File.、Output Directory for FLAASH Files、和Output Directory for FLAASH Files5、 输入成像和传感器的参数Scene center lacation 影像的中心点的经纬
33、度，可以将影像打开，查看中心点的经纬度（通过在一下窗口输入中心点的行列号即可）sensor altitude 传感器高度（轨道高度），选择正确的传感器后就可以显示了。Ground Elevation 平均海拔（所选区域的）单位是km6、atmospheric model 地球大气模型 和气溶胶模型6种标准大气模型根据以下表选择所校正区域的大气模型数据经纬度与获取时间决定选用的大气模型水气反演设置 （Water Retrieval），采用两种方式对水气进行去除a. 利用水气去除模型恢复影像中每个像元的水气量使用水气反演模型，数据必须具有 15nm以上波谱分辨率，且至少覆盖以下波谱范围之一：105
34、01210nm(优先考虑)，770870nm，8701020nm。 对于大多传感器，水气反演默认显示的是 NO，因为大多数传感器没有适当的波段来补偿水气的影响。b.单一的水气因数用于整体影像，默认是1，多光谱数据使用水气反演模型，可以在多光谱设置中手动设置水气波段气溶胶模型 （Aerosol Retrieval ）用气溶胶模型要求数据波段覆盖 660nm和 2100nm波谱。a. 提供四种标准 MODTRAN 气溶胶模型 Rural（乡村）、Urban（城市）、Maritime（海洋）、Tropospheric b. 两种气溶胶反演方法 2-Band（K-T）方法(类似模糊减少法)，如果没有找
35、到适应的黑值（一般是阴影区或者水体），系统将采用能见度值来计算；所以即使选择了该选项也要给。天气情况与能见度的关系 7.光谱打磨(高光谱) Spectral Polishing 光谱打磨(高光谱数据) 使波谱曲线更加近似于真实地物的波谱曲线 对波谱曲线进行微调8.多光谱数据参数设置当基本设置里设置了水气反演模型和气溶胶模型时，相应的在改多光谱设置框中设置参数 水气去除模型参数 气溶胶模型参数设置（用气溶胶模型要求数据波段覆盖 660nm和 2100nm波谱.） 设置值见下表所示：9高光谱数据参数设置自动选择通道定义(推荐) 设置通道定义10高级设置光谱定义文件：内置 AVIRIS、HYMAP、HYDICE、HYPERION、CASI、AISA 气溶胶高度 CO2 混合比率：390ppm 使用领域纠正 使用以前的 MODTRAN 模型计算结果 设置 MODTRAN 模型的光谱分辨率（推荐值 5 cm-1) 设置 MODTRAN 多散射模型 天顶角方位角（针对非星下点传感器）注意： 经过FLAASH大气校正后，得到的是反射率图像，且反射率值被放大10000倍。运行中如出现问题可以参考FLAASH大气校正常见错误及解决方法：