ERDAS分类后处理与ArcGIS数据交换.docx

ERDAS分类后处理与ArcGIS数据交换.docx

《ERDAS分类后处理与ArcGIS数据交换.docx》由会员分享，可在线阅读，更多相关《ERDAS分类后处理与ArcGIS数据交换.docx（41页珍藏版）》请在冰豆网上搜索。

ERDAS分类后处理与ArcGIS数据交换

基于专家知识的决策树分类

可以将多源数据用于影像分类当中，这就是专家知识的决策树分类器，本专题以ENVI中

DecisionTree为例来叙述这一分类器。

概述

基于知识的决策树分类是基于遥感影像数据及其他空间数据，通过专家经验总结、简单的数

学统计和归纳方法等，获得分类规则并进行遥感分类。

分类规则易于理解，分类过程也符合人的认知过程，最大的特点是利用的多源数据。

如图1所示，影像+DEM就能区分缓坡和陡坡的植被信息，如果添加其他数据，如区域图、道路图土地利用图等，就能进一步划分出那些是自然生长的植被，那些是公园植被。

图1专家知识决策树分类器说明图

专家知识决策树分类的步骤大体上可分为四步：

知识（规则）定义、规则输入、决策树运行和分类后处理。

1•知识（规则）定义

规则的定义是讲知识用数学语言表达的过程，可以通过一些算法获取，也可以通过经验总结

获得。

2•规则输入

将分类规则录入分类器中，不同的平台有着不同规则录入界面。

3•决策树运行

运行分类器或者是算法程序。

4•分类后处理这步骤与监督/非监督分类的分类后处理类似。

知识（规则）定义

分类规则获取的途径比较灵活，如从经验中获得，坡度小于20度，就认为是缓坡，等等。

也可以从样本中利用算法来获取，这里要讲述的就是C4.5算法。

利用C4.5算法获取规则可分为以下几个步骤：

（1）多元文件的的构建：

遥感数据经过几何校正、辐射校正处理后，进行波段运算，得到一些植被指数，连同影像一起输入空间数据库；其他空间数据经过矢量化、格式转换、地理配准，组成一个或多个多波段文件。

（2）提取样本，构建样本库：

在遥感图像处理软件或者GIS软件支持下，选取合适的图层，

采用计算机自动选点、人工解译影像选点等方法采集样本。

（3）分类规则挖掘与评价：

在样本库的基础上采用适当的数据挖掘方法挖掘分类规则，后

基于评价样本集对分类规则进行评价，并对分类规则做出适当的调整和筛选。

这里就是C4.5算法。

4.5算法的基本思路，如下：

从树的根节点处的所有训练样本DO开始，离散化连续条件属性。

计算增益比率，取GainRatio

（CO）的最大值作为划分点V0,将样本分为两个部分D11和D12。

对属性CO的每一个值产生一个分支，分支属性值的相应样本子集被移到新生成的子节点上，如果得到的样本都属

于同一个类，那么直接得到叶子结点。

相应地将此方法应用于每个子节点上，直到节点的所有样本都分区到某个类中。

到达决策树的叶节点的每条路径表示一条分类规则，利用叶列表

及指向父结点的指针就可以生成规则表。

瘵始数据”

第一次划分

第二次划分

图2规则挖掘基本思路

算法描述如下：

算法：

从空间数据集（多波段文件）中挖掘分类规则

输入：

训练样本输出：

分类规则表

方法：

一、读取数据集名字

二、读取所有的训练样本

A、读取属性信息C、原始类E、样本值A,并将样本划分为训练样本（2/3）和评价样本（1/3）。

B、属性信息C可以是连续（DISCRETE）或离散（CONTINUOUS）的，分别将属性注上这两种标记；若属性是DISCERTE，读取其可能取得值，并都存储在一个列表中；每一个属

性都有一个标记，一个给定的属性编号及初始化的取值列表均存储于一个属性的数据结构中，并将数据结构存储在一个哈希表中。

C、原始类E当作一个附加属性信息储存在属性列表中。

D、以增量方式读取每一个样本A，将所有的样本储存在一个表中，每一行代表一个样本。

三、利用数据集构建树

A、离散化连续条件属性CDISCRETE，获得的分割点集T（t1,t2……）作为条件属性C的新的取值。

B、分别计算所有条件属性的增益比率GainRatio（C）,取增益比率值最大的条件属性作为

树的划分节点，其值或范围作为划分值V（v1,v2……）来生成树的分枝。

C、判断该层与每一个等价子集的原始类类别是否一致。

若一致，生成叶子结点。

否则，继

续计算增益比率GainRatio（C）和选择条件属性C,得到树的节点和划分值V，直至所有的

样本已分类完毕。

四、测试生成树

将测试样本C'带入树中，当某一测试样本的分类预测错误时，记录分类错误的计数，并将测试样本添加到训练样本中，转向步骤三，重新构建树。

否则，输出分类树

五、抽取分类规则

到达树的叶节点的每条路径表示一条分类规则从树中抽取分类规则，打印规则和分类的详细

信息

C4.5网上有源代码下载，vc和C++版本都能获得。

DecisionTree的使用

一、规则获取

选取LandsatTM5影像和这个地区对应的DEM数据，影像和DEM经过了精确配准。

规则如下描述：

Classi（朝北缓坡植被）：

NDVI>0.3,slope<20,aspect<90andaspect>270

Class2（非朝北缓坡植被）：

NDVI>0.3,slope<20,90<=aspect<=270

Class3（陡坡植被）：

NDVI>0.3,slope>=20,

Class4（水体）：

NDVI<=0.3,0

Class5（裸地）：

NDVI<=0.3,b4>=20

Class6（无数据区，背景）：

NDVI<=0.3,b4=0

也可以按照二叉树描述方式：

第一层，将影像分为两类，NDVI大于0.3，NDVI小于或等于

0.3；第二层，NDVI高的，分为坡度大于或等于20度和坡度小于20度。

以此往下划分。

二、输入决策树规则

打开主菜单->classification-'DecisionTree->BuildNewDecisionTree,如图3所示，默认显示了一个节点。

图3DecisionTree界面

首先我们按照NDVI的大小划分第一个节点，单击Nodel，跳出图4对话框，Name为

NDVI>0.3，在Expression中填写：

{ndvi}gt0.3。

图4添加规则表达式

点击0K后，会提示你给｛ndvi｝指定一个数据源，如图5所示，点击第一列中的变量，在对话框中选择相应的数据源，这样就完成第一层节点规则输入。

聖Variable/FilePairings口

图5指定数据源

Expression中的表达式是有变量和运算符（包括数学函数）组成，支持的运算符如表1所示

表达式

部分可用函数

基本运算符

+、-、*、/

三角函数

正弦Sin（x）、余弦cos（x）、正切tan（x）

反正弦Asin（x）、反余弦acos（x）、反正切atan（x）双曲线正弦Sinh（x）、双曲线余弦cosh（x）、双曲线正

切tanh（x）

小于LT、小于等于LE、等于EQ、不等于NE、大于

关系/逻辑

等于GE、大于GT

and、or、not、XOR

最大值（>）、最小值（<）

指数（人）、自然指数exp

自然对数对数alog（x）

其他符号

以10为底的对数alog10（x）整形取整round（x）、ceil（x）

平方根（sqrt）、绝对值（adb）

表1运算符

ENVI决策树分类器中的变量是指一个波段的数据或作用于数据的一个特定函数。

变量名必

须包含在大括号中，

即｛变量名｝；或者命名为bx,x代表数据，比如哪一个波段。

如果变量

被赋值为多波段文件，

变量名必须包含一个写在方括号中的下标，表示波段数，比如｛pc[2]｝

表示主成分分析的第

•主成分。

支持特定变量名如表2，也可以通过IDL自行编写函数。

变量

作用

slope

计算坡度

aspect

计算坡向（南坡北坡或者东西方向）

ndvi

计算归一化植被指数

Tascap[n]

穗帽变换，n表示获取的是哪一分量。

pc[n]

主成分分析，n表示获取的是哪一分量。

lpc[n]

局部主成分分析，n表示获取的是哪一分量。

mnf[n]

最小噪声变换，n表示获取的是哪一分量。

Lmnf[n]

局部最小噪声变换，n表示获取的是哪一分量。

Stdev[n]

波段n的标准差

lStdev[n]

波段n的局部标准差

Mean[n]

波段n的平均值

lMean[n]

波段n的局部平均值

Min[n]、max[n]

波段n的最大、最小值

lMin[n]、lmax[n]

波段n的局部最大、最小值

表2变量表达式

第一层节点根据NDVI的值划分为植被和非植被，如果不需要进一步分类的话，这个影像就

会被分成两类：

classO和classl。

对NDVI大于0.3，也就是classl，根据坡度划分成缓坡植被和陡坡植被。

在classl图标上右键，选择AddChildren。

单击节点标识符，打开节点属性窗口，Name为

Slope<20，在Expression中填写：

{Slope}It20。

同样的方法，将所有规则输入，末节点图标右键EditProperties，可以设置分类结果的名称

和颜色，最后结果如图6所示。

图6规则输入结果图

三、执行决策树

选择Options->Execute，执行决策树，跳出图7所示对话框，选择输出结果的投影参数、重采样方法、空间裁剪范围（如需要）、输出路径，点击OK之后，得到如图8所示结果。

在决策树运行过程中，会以不同颜色标示运行的过程。

聖Dgwo貝TreeExecutionPau…

0K

图7输出结果

图8决策树运行结果

回到决策树窗口，在工作空白处点击右键，选择ZoomIn，可以看到每一个节点或者类别

有相应的统计结果（以像素和百分比表示）。

如果修改了某一节点或者类别的属性，可以左

键单击节点或者末端类别图标，选择Execute，重新运行你修改部分的决策树。

图9运行决策树后的效果

分类后处理和其他计算机分类类似的过程。

ERDAS分类后处理与ArcGIS数据交换

分类后，得到一张专题图。

由于ERDAS是基于像素的分类，在分类结果中会出现很多细小

的杂点以及小的图斑，如何去掉这些杂点和小图斑（称为ERDAS后处理），以及如何将专

题图转换成ArcGIS下的shp格式，是本文需要解决的工作。

一、ERDAS后处理

1首先给出原始图像和经过监督分类得到的分类图像，如下图：

原始影像

监督分类结果影像

2、开始聚类

使用Interpreter模块->GISAnalysis->Clump，打开聚类界面，输入影像为分类后影像，给一个输出影像，聚类邻域可以选择8或者4，点击0K,即可进行聚类，结果如下：

聚类结果

3、去除杂点和小图斑

使用Interpreter模块->GISAnalysis->Eliminate，也可以使用Sieve，后者容易出现一些属性值为0的像素，不太好处理，所以这里使用Eliminate，输入影像为上一步聚类后的影像，给定一个输出影像，Minimum处设置杂点和小图斑的像素大小，这里可以根据实际情况设

定，这里选择100，点击0K，结果如下。

通过比对，可以发现分类结果中很多杂点和小图斑被有效地去除了。

去除小图斑

Tievtcf2=惠鳶樂乱:

（sLaycr_I）.L*

Z.jl-1J[Ll13tr£3aaAQI^aitvrHilp

◎■町田将翳器口+%丘X.瞰T肯沖

去除杂点和小图斑

4、形态学处理

也可以使用形态学的腐蚀、膨胀、开运算、闭运算来消除杂点，毛刺，空洞等。

使用Interpreter

模块->Utilities->Morphological，选择输入输出影像，结构元素kerneldefinition可以选

择3*3，5*5或7*7，也可以自定义，function中Erode是腐蚀，Dilate是膨胀，Open是开运算，Close是闭运算，选择Open，然后点击OK即可，只是得到的结果不是专题图，需要重新分类得到专题图。

不过这个时候的分类，将会非常简单，不再赘述。

形态学开运算

二、ArcGIS数据交换

将分类结果转换成ArcGIS的矢量格式，即可在ArcGIS上进行后续的处理。

使用Vector模

块->RastertoVector，将分类结果转换成Arcinfo的ArcCoverage格式。

该格式可以在

ArcMap中直接打开，如下图：

ArcMap中打开ArcCoverage格式

ArcloulbciiK

+JFromRas.ter

障0FrflmWFS+jiMetadata

+TqCAD

+ToCoverage

+％dBASE

+TaGeodatabase

+％EML

+JiToRaster

Shapefile

?

eatureClassToShapefile（multiple）

To囂I

BataIWanag-wnentTool's

*

+LinearReferencingTools

+ihiliidlirhemionTools

+Samples+ServerTooAs

-SpatialStfiiiceTools

也可以使用ArcToolBox中的ConversionTools->ToShapefile->FeatureClasstoShapefile（multiple）,将ArcCoverage格式转换成shp格式。

ArcToolBox

ArcMap中打开的shp格式，图中黑点是注记层

剩下的事情，就是通过ArcGIS提供的强大的空间分析功能，提取各种地物的各种信息了。

ARCGIS转移矩阵

土地利用转移矩阵生成的几种方法

查阅相关的资料，也没有得到土地利用类型转换矩阵确切的定义,一区域内土地利用类型的相互转换关系，一般用二维表来表达,个地类间相互转化的具体情况。

比如某一类别的土地有百分之多少

了其他的土地类型，现在某类型的土地分别是由过去的哪些类别转化而来的等等。

还可以生成变化统计栅格图（掩膜图像），它描述了前后两幅土地分类图之间的地类发生转变的位置和类别。

土地利用类型转换矩阵可以从两幅栅格图中计算得到，也可以从两个矢量文件中计算获得。

下面介绍在ENVI下从两幅分类结果的栅格图中计算土地利用类型转换矩阵。

1准备数据

两个时相的土地利用分类结果，它是单波段、专题类型的伪彩色图像（ENVIClassification）。

2、计算转换矩阵

打开两个土地利用分类结果。

（1）在主菜单中，选择BasicTools宀ChangeDetectionChangeDetectionStaCsti

（2）分别在InitialState对话框和finalstate对话框中选择前一时相和后一时相的土地利用结果。

（3）在DefineEquivalentClasses对话框中（图1）,如果两个土地利用分类名称一致，系

统自动将InitialStateClass和FinalStateClass对应，否则手动选择，单击AddPair。

（4）选择对应的地物类型之后，单击0K按钮，出现图2对话框。

选择生成图表表示单位

（ReportType）:

像素（Pixels）、百分比（Percent）和面积（Area）。

选择OutputClassificationMaskImages?

为YES，输出掩膜图像，选择输入路径及文件名。

（5）单击OK,执行土地利用类型转换矩阵计算过程。

警DefineEquivalentClasses

图1DefineEquivalentClasses对话框

图2选择数据参数

3、查看结果

（1）如图3为得到的土地利用类型转换矩阵结果。

横字段表示前一时间段（InitialState）

的土地利用类别，纵字段为后一时间段（FinalState）的土地利用类别。

横字段和纵字段交

叉处表示变化值，如有2520900平方米林地用地变化为草地。

图3土地利用类型转换矩阵

（2）还可以为每一个地类生成一个变换掩膜图像，图4所示为其中一个地类的掩膜图像。

掩膜图像的灰度值表示变化类型，如这里的2｛草地｝表示林地变化为草地的像元。

图4变化掩膜图像

根据你的数据类型选用不同的数据生成方法

ERDAS中vector数据转移矩阵方法

数据是Vector格式（shp数据格式）

1ErdasImagine----Interpreter---GisAnalysis---Matrix，输入两个时相的Vector数据即可此时注意：

输出栅格大小不应设的太小要不一运行就会提示你的空间不足

做这一步之前，请做好前期的地理编码。

2ArcView3.3加载spatialanalysis模块

把两时相的Vector图转成grid格式（当然中间有一些单位的设置根据你做的图的分辨率来设置即可）analysis-mapcalculate直接计算即可。

3把两期解译完的Vector文件在arctoolbox-overlay-union中叠加，注意：

两个文件不能用同一个字段名，比如一个用93Type，另一个时相则用OOType

叠加后的文件在Arcmap中打开，选中文件，然后点右键Property空间查询，输入

条件语句，比如：

93Type='TAn00Type='2';查询结果即为第一种类型转化为第二种类型的图形，可以另建一图层比如：

12，把查询结果复制到12图层上。

统计出面积，依进行，

就可以得到土地利用类型转移矩阵。

一、数据准备（图1）

准备两幅不同时相的土地利用现状图（shp格式），每幅图的属性表都要有一个表示土地利

用类型的字段，并且要使用不同的名称加以区分，如Type1995,Type2000。

土地利用类型

名称必须统一，并且完整，如都使用城镇用地”、有林地”等。

二、数据融合（图2）DISSOLVE

在ArcMap里分别打开两个时相的图层，打开ArcToolbox，选择DataManagementTools|Generalization|Dissolve工具。

InputFeature选择要融合的图层，OutputFeatureClass选择输出结果存储的位置及名称，DissolveField（s）选择土地利用类型字段（如Type1995）,然后勾选Creatmultipartfeatures选项，点击OK完成。

重复此过程，对另一时相数据进行融合。

此步骤使相同利用类型的记录融合为一个记录，以提高后面步骤的计算速度。

三、叠置分析（图3）OVERLAY

在ArcMap中打开两个时相融合后的数据，在ArcToolbox中选择AnalysisTools|Overlay|

Intersect工具，InputFeatures选择两个时相的图层，OutputFeatureClass选择叠加结果存储

的位置及名称，其余选项可以忽略，单击【OK】完成。

在Editer工具条中选择Editer|StartEditing,然后在属性表中NewArea字段上单击右键选择

CalculateGeometry…，•在打开的CalculateGeometry对话框中，Property选择Area,Units

选择要使用的面积单位，单击【0K】完成图斑面积计算。

依次选择Editer|SaveEdits/End

Editing保存和退出编辑状态。

在属性表中选择Option|Export…将属性表保存为dbf文件。

五、制作转移矩阵（图7-10）（以Excel2007为例）

在Excel中打开上一步保存的dbf,另存为Excel格式并打开。

在Excel中选中所有数据（不要点左上角，只选择有效数据），点击【插入】选项卡，选择【数据透视表】|【数据透视表】，点击【确定】。

在打开的数据透视表中按图示将字段拖入相应区域。

UMI

J

i*^->1rjh■—4i

«g|

!

l：

勺丨■■■■■▼・■呼

斗

祖廿佗比2料/用•耳舗

Excel自动计算矩阵，将该表稍事整饰就得到美观的土地利用转移矩阵。

矩阵中r（l,j）就表示

i类型向j类型转移的土地面积，空值表示i类型向j类型没有转移。

料I

車0

irti总踰n

KKARHU

OlOMH

ki.ni

・

««M9

-■*上

*7■"，s.—亠土】11■

亦1

■]|]-l■!

■

*y曲■.■帶It啊冨

p巨.

BlfL,

2f.T

_EB

■hi

*-

■■旨■■

草却l事期*揖*

KWflt力曜欽毎戒・（¥H>

<1*1*Srfftftit

IkUAl^l

OkHMIl林mm

13"

«.td

*34料訐

0如

ti0m?

4

&

0iJMtl

□■f見・

呼心.3

114««■«««_+|f■毒t早

ft”購用，

*«fflt*14M1

理启埔LM

*工*事厦i

■M

WnE»«

Jt电皙喪Rit

M4jWn*JT*t<床鼻（罕JCi理电

WW4=#减

LT±4*

Ci44L44

t.M

rf--7-<

Fulu

土地利用转移矩阵生成的几种方法

根据你的数据类型选用不同的数据生成方法

若你的数据是Raster格式：

则有如下方法

1ErdasImagine-lnterpreter-GISAnalysis-Matrix，

输入两个时相的Raster数据即可做这一步之前记得：

先对两时相的数据进行重编码

（Interpreter---GISAnalysis---Recode）

一般运行如果出现错误肯定是重编码没做好，请继续查证。

2先在Erdas中利用Modeler计算如下公式

NC（l,J）=NC（l）*10+NC（J）,（J>l）

其中：

NC（I,J）表示i，j两年份的土地利用变化图；NC（i）表示i年份遥感分类影像；NC

（j）表示j年份的遥感分类影像。

在此计算的基础上，将以上变化影像图转化为BIL格式，再利用ARC/INFOGRID模块将影

像转为GRID格式，然后利用GRID模块中的属性表（vat）查看命令对影像灰度值进行统计，