遥感影像分类方法实验报告

./实验报告实验报告目录TOC\o"1-3"\h\u224241实验目的4139662实验数据475933实验容475094实验步骤5318804.1对人口矢量数据〔shapefile进行投影转换532574.1.1Census.shp文件投影坐标的检查 5147144.1.2将投影坐标转换为WGS_1984_UTM_Zone_16N 6100314.2对遥感影像进行几何精校正〔以经过投影变换的人口矢量数据为基准677634.2.1Census.shp在ENVI软件的加载 6254854.2.2对遥感影像进行几何精校正〔以矢量数据为基准 745154.2.3用矢量图层对遥感影像进行裁剪 10217624.3将Pan波段和多光谱波段进行融合,并对融合效果进行定性和定量评价1133784.3.1两种融合方法的原理11235964.3.2进行Gram-SchmidtSpectralSharpening融合11197924.3.4融合效果进行定性评价14208894.3.5融合效果进行定量评价<软件提供的计算方法>15101924.3.6融合效果进行定量评价<Matlab编程计算>16321034.3.7遥感影像融合定量分析代码20121924.4生成住房密度栅格影像2333074.4.1两表的连接 23128924.4.2计算房屋密度 24323854.4.3直接栅格化 25182614.4.4IDW插值 25319944.4.5对房屋密度图进行重分类 2691804.5将住房密度栅格影像作为额外的通道与ETM+多光谱波段进行叠加26120864.6监督分类〔融合方法为HSV,波段为5,4,327181294.6.1打开GoogleEarth影像作为监督分类的参照 27261064.6.2建立兴趣区 29240054.6.3训练样区的选择 3084444.6.4训练样区的评价 3162064.6.5执行监督分类 33280864.6.6分类后处理 354574.6.7评价结果分析 37300894.6.8分类结果面积统计 385244.6.9分类结果 41323044.7分类结果评价与分析41230254.7.1未加入房屋密度图层的分类结果评价与分析 41303574.7.2加入IDW插值房屋密度图层的分类结果评价与分析 42285764.7.3加入直接栅格化房屋密度图层的分类结果评价与分析 43152394.7.4加入重分类后IDW插值房屋密度图层的分类结果评价与分析 4417714.7.5从总精度与Kappa系数对分类结果进行评价 45232654.7.6分类结果总体评价 46128744.7.7与其他训练样区的分类精度和Kappa系数的计算 4843824.8决策树分类49267894.8.1决策树分类原理 49236724.8.2数据预处理 49200884.8.3指数的计算 5166514.8.4执行决策树 54299794.8.5不同参数设置的对比 5710775实验体会60105385.1实验中存在的问题60267335.2软件平台使用63290375.3实验总结631实验目的①掌握ArcGIS10和ENVI4.7对遥感图像处理的基本操作与原理②熟悉几何精校正的方法,掌握ENVI软件对遥感影像进行几何精校正③掌握全色波段与多光谱波段的融合方法和原理,学会对融合效果进行定性定量分析④熟悉掌握ArcGIS的栅格化方法和IDW插值方法⑤熟悉监督分类的方法和基本原理,掌握ENVI软件中进行监督分类⑥了解监督分类后评价过程,对分类结果进行精度评价和分析⑦掌握Erdas的空间建模方法以及原理⑧了解RuleGen算法,掌握决策树分类方法2实验数据①带属性数据的shapefile：Census.shp②带有陆地面积字段的矢量图层：③GoogleEarth_原始拼接：GE1005211134.jpg④研究区域的多光谱波段数据：Stack_b1-6162-7.img⑤研究区域的全色波段数据：b8.img⑥监督分类参照影像：GoogleEarth3实验容①对人口矢量数据〔shapefile进行投影转换：WGS_1984_UTM_Zone_16N②对遥感影像进行几何精校正〔以经过投影变换的人口矢量数据为基准：〔1对多光谱波段〔30m空间分辨率进行几何精校正〔小于0.25个像元；〔2对Pan波段〔15m空间分辨率进行几何精校正〔小于0.25个像元；③将Pan波段和多光谱波段进行融合〔自选至少一种融合算法,并对融合效果进行定性和定量评价；④生成住房密度栅格影像：〔1直接栅格化；〔2IDW插值；⑤将住房密度栅格影像作为额外的通道〔或波段与ETM+多光谱波段进行叠加；⑥进行监督分类和分类后处理〔Post-Classification,ExpertRules⑦利用ERDAS软件的空间建模〔SpatialModeler进行水体信息〔MNDWI指数和植被信息〔NDVI指数的提取；⑧利用"自动阈值决策树分类算法"进MarionCounty的土地利用/覆盖分类信息提取〔使用的数据：原始各波段+MNDWI+NDVI+ISODATA等,或其他有益的波段组合①探讨"自动阈值决策树分类算法"中的各个参数意义及如何设置更合理②对分类结果进行评价与分析⑨对分类结果进行精度评价和分析；4实验步骤4.1对人口矢量数据〔shapefile进行投影转换4.1.1Census.shp文件投影坐标的检查根据实验要求,人口矢量数据〔shapefile进行投影坐标应为：WGS_1984_UTM_Zone_16N在ArcGIS软件的图层右击Properties,在LayerProperties的Source下查看投影信息,如图1。得到Census.shp的投影坐标为：GCS_North_American_1983,与实验要求不符合,需进行投影转换。图14.1.2将投影坐标转换为WGS_1984_UTM_Zone_16N在Arcgis软件的工具箱中的DefineProjection工具,设置输入数据为：Census.shp,坐标系统为GCS_WGS_1984,在工具箱中的工具,设置输入数据：Census.shp,导入遥感影像的投影坐标系,即GCS_WGS_1984〔如图2。图24.2对遥感影像进行几何精校正〔以经过投影变换的人口矢量数据为基准4.2.1Census.shp在ENVI软件的加载在ENVI软件中,File->OpenVectorFile,选择Census.shp,设置好参数,生成evf文件〔如图3。图34.2.2对遥感影像进行几何精校正〔以矢量数据为基准在ENVI分别打开遥感影像img和矢量文件vef,选择Map->Registration->SelectGCPs:imagetomap,以矢量数据为基准,设置好投影,如图4。图4在ENVI中,在zoom窗口下采集控制点,这次实验采集的控制点数为13个,控制点的主要定位在道路与道路之间的交叉点,如图5,其控制点的RMSError为0.246390,如图6,如图7为20个控制点的采集情况。图5图6图7选择校正参数输出结果,在GroundControlPointsSelection窗口选择Option->WarpFile,如图8；数学模型为Polynomial,设定参数为2,从采样方法为最临近法,如图9。图8图9如图10、11为几何校正前后,矢量图层与遥感影像吻合度的对比,可以明显看出,经过几何校正后的遥感影像与矢量图层吻合程度有明显的改善,有部分水体边界不吻合,这主要是由于水体会随时间而改变；而街区与道路吻合程度良好。图10图11如图12为对多光谱波段〔30m空间分辨率进行几何精校正后的遥感影像；如图13为对Pan波段〔15m空间分辨率进行几何精校正后的遥感影像。图12图134.2.3用矢量图层对遥感影像进行裁剪在ENVI中打开Census.shp文件,将该图层转换成ROI,如图14。图14通过ROI进行裁剪遥感影像,选择BasicTools->SubsetDataviaROIs；选择转换好的ROI进行裁剪,如图15。图15同样,对pan波段的遥感影像进行裁剪,得到遥感影像如图16。图164.3将Pan波段和多光谱波段进行融合,并对融合效果进行定性和定量评价4.3.1两种融合方法的原理Gram-Schmidt可以对具有高分辨率的高光谱数据进行锐化。①从低分辨率的波谱波段中模拟出一个全色波段。②对该全色波段和波谱波段进行Gram-Schmidt变换,其中模拟的全色波段被作为第一个波段。③用Gram-Schmidt变换后的第一个波段替换高空间分辨率的全色波段。④应用Gram-Schmidt反变换构成pan锐化后的波谱波段。Gram-SchmidtSpectralSharpening方法进行图像增强能够比较好的保留原多光谱图像的光谱信息,使遥感影像的融合保留多光谱影像的增强效果。用PC可以对具有高空间分辨率的光谱图像进行锐化。①先对多光谱数据进行主成分变换。②用高分辨率波段替换第一主成分波段,在此之前,高分辨率波段已被缩放匹配到第一主成分波段,从而避免波谱信息失真。③进行主成分逆变换。函数自动地用最近邻、双线性或三次卷积技术将多光谱数据重采样到高分辨率像元尺寸。4.3.2进行Gram-SchmidtSpectralSharpening融合在ENVI软件中,选择Transform->ImageSharpening->Gram-SchmidtSpectralSharpening,在SelectLowSpatialResolutionMultiBandInputFile中选择Stack_b1-6162-7_CJ.img多光谱波段,SpatialSubset为FullScene,SpectralSunset为8个波段,如图17。图17在SelectHighSpatialResolutionPanInputFile窗口中选择b8_CJ.img全色波段,如图18。图18SelectMethodforLowResoutionPan选择AverageofLowResolutionMultispectralFile：利用多光谱波段的平均值来模拟低分辨率的全色波段。Resampling选择NearestNeighbor,OutputResult选择保存路径,如图19。图194.3.3进行PCSpectralSharpening融合

在ENVI软件中,选择Transform->ImageSharpening->PCSpectralSharpening,在SelectLowSpatialResolutionMultiBandInputFile中选择Stack_b1-6162-7_CJ.img多光谱波段,SpatialSubset为FullScene,SpectralSunset为8个波段,如图20。图20在SelectHighSpatialResolutionPanInputFile窗口中选择b8_CJ.img全色波段,如图21。图21Resampling选择NearestNeighbor,OutputResult选择保存路径,如图22。图224.3.4融合效果进行定性评价如图23为原始遥感影像与经过Gram-SchmidtSpectralShaping处理后的影像。图23从图23,我们可以得到,Gram-SchmidtSpectralShaping处理后的影像总体上来说分辨率有很大的提高,清晰度高,光谱信息比较丰富,但颜色的匹配还不是很理想,整个影像的色调基本上一致,呈现出泛红的现象,地物之间的辨别基本上是通过影像上的灰度信息,而色彩提供的信息量较少,区分度不高。如图24为原始遥感影像与经过PCSpectralSharpening处理后的影像。图24从图24,我们可以得到,PCSpectralSharpening处理后的影像清晰度高,涵盖的地物信息量大,颜色的畸变很小,和Gram-SchmidtSpectralShaping处理的影像效果类似,颜色的匹配还不是很理想,整个影像的色调基本上一致,呈现出泛红的现象,色彩不够丰富。在水域的地方颜色与原来的为黑色变成了青色,颜色变化差异较大。定性评价结论,从图23、24,我们可以了解到提高空间分辨率效果最好的是Gram-SchmidtSpectralShaping；光谱变化较小的是PCSpectralSharpening,但在水域的区域上,Gram-SchmidtSpectralShaping保持这原始影像的色彩。4.3.5融合效果进行定量评价<软件提供的计算方法>在实验中,波段的选取为4,3,2波段。在主窗口中,右击选择ZProfile,对比三个影像的情况<如图25>。图25从图25,我们可以得到,经过融合后,在Value上都有所增加,band1从90->100；说明了融合后影像的亮度信息提高了。但PC和GS两者的Value通过图上很难分辨出区别。图26XProfile图27YProfile在ENVI4.7软件中,选择BasicTools->Stayistics->ComputeStatistics,如图28。图28如表1为原始影像、PC融合后影像和GS融合后影像的部分统计参数。表1图29 图30从图29、30,两折线图上我们可以直观的了解原始影像、PC融合后影像和GS融合后影像在亮度信息上的关系。在均值上,PC融合影像和GS融合影像均有所增加,但是从幅度上来看,PC融合影像的均值增加的幅度较大,较GS融合影像更有利于目视判读。在方差上,与原始影像相比较PC融合影像的方差有所增加,而GS融合影像的方差大幅下降,不利于目视判读。因此,在亮度指标上最好的为PC融合,其次为GS融合。4.3.6融合效果进行定量评价<Matlab编程计算>此次实验从融合影像的亮度信息,清晰程度,光谱保持程度,信息丰富程度等多角度进行评价分析,相比传统的单一定量评价全面,能够减少评价的随机性,使得定量评价更加科学全面。主要通过4方面进行统计分析：①亮度信息,针对融合后影像亮度信息进行评价,主要包括均值和标准差等指标。②清晰度,评价融合影像的清晰程度,主要包括平均梯度和空间频率等指标。③光谱信息,评价融合后影像的光谱变形情况,包括扭曲程度,相关系数等指标。④信息量,评价融合后信息量的保持或增加程度,主要包括信息熵,交叉熵和互信息等指标。以下公式中,影像函数为Z<X,Y>,影像的行数和列数分别为M和N,影像的大小则为M×N。4.3.6.1亮度信息均值<z>影像均值是像素的灰度平均值,对人眼反映为平均亮度。如果均值适中,则影像效果良好,其定义为:标准差<σ>标准差反映了影像灰度相对于灰度平均值的离散情况,在某种程度上,标准差也可用来评价影像反差的大小。若标准差大,则影像灰度级分布分散,影像的反差大,可以看出更多的信息,其公式为：方差<D>方差反映了灰度相对于灰度均值的离散情况,方差越大,则灰度级分布越分散。此时,图像中所有灰度级出现的概率越趋于相等,从而包含的信息量趋于越大。方差计算公式为：4.3.6.2清晰度平均梯度<G>平均梯度可敏感地反映影像对微小细节反差表达的能力,可用来评价影像的清晰程度,还可同时反映出影像中微小细节反差和纹理变化特征,其计算公式为：空间频率<SF>空间频率反映了一幅影像空间的总体活跃程度,空间频率越大,说明融合效果越好。它包括空间行频率RF和空间列频率CF组成,其公式为：总体的空间频率值取RF和CF的均方根,即：4.3.6.3光谱信息扭曲程度<D>扭曲程度D是指融合影像F像素灰度平均值与源影像A像素灰度平均值之差,也可以说是融合影像与源影像的差值影像的灰度平均值,它的表达式为：它反映融合影像和源多光谱影像在光谱信息上的差异大小和光谱特性变化的平均程度。相关系数<ρ>融合影像与源影像的相关系数能反映两幅影像光谱特征的相似程度,其定义如下：其中,f和a分别为融合影像与源影像的均值。通过比较融合前后的影像相关系数可以看出多光谱影像的光谱信息的改变程度。4.3.6.4信息量信息熵<E>影像的信息熵值是衡量影像信息丰富程度的一个重要指标。对于一幅单独的影像,其灰度分布为,为灰度值等于i的像素数与影像总像素数之比。根据Shannon信息论的定理,一幅影像的信息熵为:融合影像的熵越大,融合影像所含的信息越丰富,融合质量越好。交叉熵<C>交叉熵可以用来测定两幅影像灰度分布的信息差异。设源影像和融合影像的灰度分布分别为：和,则交叉熵定义为:交叉熵是评价两幅影像差别的重要指标,它直接反映了两幅影像对应像素的差异。对融合影像前的源影像和融合结果影像求交叉熵,即可得到融合影像与源影像的差异。差异越小,则该融合方法从源影像提取的信息量越多。互信息<MI>互信息是信息论中的一个重要概念,它可作为两个变量之间相关性的量度,或一个变量包含另一个变量的信息量的量度。F与A,B的互信息分别表示为和：式中：,和分别是A,B和C的概率密度；和分别代表两组影像的联合概率密度。互信息的值越大,表示融合影像从源影像中获取的信息越丰富,融合效果越好。通过定量评价分析对2种融合方法进行评价。评价容包括亮度信息、清晰度、光谱信息和信息量等4个方面。试验证明这些评价参数能够很好地反映遥感影像融合效果。融合结果参数评价指标见表所示。亮度信息清晰度均值方差标准差平均梯度空间频率b873.0791488.295122.09743.46829.0288Stack90.48971488.73338.58417.645620.0627PC90.67711464.02738.26266.107315.398GS90.53151474.09938.39406.134415.4558光谱信息信息量扭曲程度相关系数信息熵交叉熵互信息b85.5133Stack6.8639PC0.18750.87336.87130.23680.0409GS0.04190.88586.90860.77940.04051>亮度信息由表第一列亮度信息可以看出2种融合方法的均值有所降低,方差和标准差有所增加；2>清晰度由表第二列可以看出,3种融合方法的平均梯度和空间频率都小于原始的多光谱影像,说明融合后的结果影像减少了细节纹理信息,使影像的清晰度降低了。3>光谱信息依据表中的第三列可以看出GS融合扭曲程度最小,而相关系数最大,因此GS融合的光谱信息保持得最好。4>信息量通过信息熵可以看出GS融合所得的结果信息量最大；从交叉熵看到,GS融合和PC融合的结果交叉熵依次减小,说明融合后的结果影像与源影像对应像素差异也依次减小,即GS融合保留原始图像信息量最多,PC融合保留最少；从互信息也可以得到同样规律。而通过定性评价分析不能得出图像携带的信息量,它仅仅是通过比较分析图像的亮度信息进行比较评价的。4.3.7遥感影像融合定量分析代码Dfusion=imread<'C:\Users\Administrator\Desktop\RS1\GS.tif'>;Dlow=imread<'C:\Users\Administrator\Desktop\RS1\Stack_b1-6162-7_CJ.tif'>;Dhigh=imread<'C:\Users\Administrator\Desktop\RS1\b8_CJ.tif'>;[rh,ch]=size<Dhigh>;Dlowh=imresize<Dlow,[rh,ch],'bicubic'>;mean=mean2<Dlow>;%均值meanf=mean2<Dfusion>;DIF=abs<meanf-mean>; %扭曲程度Std=std2<Dfusion>; %标准差Std2=std2<Dlowh>;Ds=Std^2; %方差Dl=Std2^2;D3=abs<Ds-Dl>; %求差方差p=imhist<Dfusion<:>,8>/numel<Dfusion<:>>;r=entropy<Dfusion<:>>;%信息熵c=corr2<Dfusion<:>,Dlowh<:>>; %相关系数h1=diff<Dfusion>;%求影像差分,反映图像清晰度h=mean2<h1>;g1=diff<Dlow>;g=mean2<g1>;G=h-g;%求融合后影像与原影像差分差值s=size<size<Dlow>>;%交叉熵ifs<2>==3;%判断是灰度图还是RGBf1=rgb2gray<Dlow>;f2=rgb2gray<Dfusion>;else f1=Dlow; f2=Dfusion;endG1=double<f1>;G2=double<f2>;[m1,n1]=size<G1>;[m2,n2]=size<G2>;m2=m1;n2=n1;X1=zeros<1,256>;X2=zeros<1,256>;result=0;%统计两图各灰度级像素fori=1:m1 forj=1:n1 X1<G1<i,j>+1>=X1<G1<i,j>+1>+1; X2<G2<i,j>+1>=X2<G2<i,j>+1>+1; endendfork=1:256 P1<k>=X1<k>/<m1*n1>; P2<k>=X2<k>/<m1*n1>; if<<P1<k>~=0>&<P2<k>~=0>> result=P1<k>*log2<P1<k>/P2<k>>+result; endendf=results=size<size<Dlow>>;%互信息ifs<2>==3;%判断是灰度图还是RGBa=rgb2gray<Dlow>;b=rgb2gray<Dfusion>;enda=double<Dlow>;b=double<Dfusion>;[Ma,Na]=size<a>;[Mb,Nb]=size<b>;M=min<Ma,Mb>;N=min<Na,Nb>;%初始化直方图数组hab=zeros<256,256>;ha=zeros<1,256>;hb=zeros<1,256>;%归一化ifmax<max<a>>~=min<min<a>> a=<a-min<min<a>>>/<max<max<a>>-min<min<a>>>;else a=zeros<M,N>;endifmax<max<b>>-min<min<b>> b=<b-min<min<b>>>/<max<max<b>>-min<min<b>>>;else b=zeros<M,N>;enda=double<int16<a*255>>+1;b=double<int16<b*255>>+1;%统计直方图fori=1:M forj=1:N indexx=a<i,j>; indexy=b<i,j>; hab<indexx,indexy>=hab<indexx,indexy>+1; %联合直方图 ha<indexx>=ha<indexx>+1; %a图直方图 hb<indexy>=hb<indexy>+1; %b图直方图 endend%计算联合信息熵hsum=sum<sum<hab>>;index=find<hab~=0>;p=hab/hsum;Hab=sum<sum<-p<index>.*log<p<index>>>>;%计算a图信息熵hsum=sum<sum<ha>>;index=find<ha~=0>;p=ha/hsum;Ha=sum<sum<-p<index>.*log<p<index>>>>;%计算b图信息熵hsum=sum<sum<hb>>;index=find<hb~=0>;p=hb/hsum;Hb=sum<sum<-p<index>.*log<p<index>>>>;%计算a和b的互信息mi=Ha+Hb-Hab%计算a和b的归一化互信息mi1=hab/<Ha+Hb>; img=double<img>;%平均梯度%Getthesizeofimg[r,c,b]=size<img>;dx=1;dy=1;fork=1:bband=img<:,:,k>;[dzdx,dzdy]=gradient<band,dx,dy>;s=sqrt<<dzdx.^2+dzdy.^2>./2>;g<k>=sum<sum<s>>/<<r-1>*<c-1>>;endoutval=mean2<g>;[m,n]=size<Dfusion>; %空间频率s=size<size<Dfusion>>;ifs<1>==3;%判断是灰度图还是RGBimg=rgb2gray<Dfusion>;endimg=double<img>;rf=0.00;cf=0.00;fori=1:mforj=2:nrf=rf+<img<i,j>-img<i,j-1>>*<img<i,j>-img<i,j-1>>;endendRF=rf/<m*n>;fori=2:mforj=1:ncf=cf+<img<i,j>-img<i-1,j>>*<img<i,j>-img<i-1,j>>;endendCF=cf/<m*n>;MSF=<RF+CF>.^<1/2>;4.4生成住房密度栅格影像4.4.1两表的连接在ArcGIS软件中,选择OpenAttributeTable中查看属性数据,其中,没有该图层的面积字段。因此需要进行有陆地面积的表格进行连接,如图31。图31在JoinData数据中选择两表连接的属性分别为STFID和BLKIDFP00,如图32。图324.4.2计算房屋密度创建房屋密度字段hou_den,并进行字段计算,公式为：hou_den=[HSE_UNITS]*1000000/[ALAND00]其中,[HSE_UNITS]为各街区房屋单元数,[ALAND00]为各街陆地面积,单位平方米。乘以1000000是为将房屋密度的单位转换为个/平方千米,如图33。图334.4.3直接栅格化在ArcGIS工具箱中选择工具PolygontoRaster,输入栅格数字的字段为hou_den,输出栅格大小为15m,如图34。图344.4.4IDW插值在ArcGIS工具箱中选择工具FeaturetoPoint,不勾选Inside,将选择多边形的中心点作为生成的点,如图35。图35在ArcGIS软件的空间分析模块中,打开IDW插值方法的窗口,对房屋密度的点图进行插值,选择字段为hou_den,输出栅格大小15m,如图36。图364.4.5对房屋密度图进行重分类在ArcGIS软件的空间分析模块中,打开Reclassify,如图37,将房屋密度分为0-400,400-1300,>1300三个级别。IDW插值房屋密度图和直接栅格化房屋密度图的重分类结果分别如图38、图39。图37图38图394.5将住房密度栅格影像作为额外的通道与ETM+多光谱波段进行叠加在ArcGIS10软件中,分别将两种方法生成的住房密度栅格影像数据导出,格式为dat,如图40。图40在Arcgis10软件中,分别把ras、Re_ras、idw和Re_idw导出成dat文件。在ENVI4.7软件中,分别打开影像和房屋密度栅格图,选择BasicTools->LayerStacking,如图41。图41在LayerStackingParameters窗口中,点击ImportFile将idw.dat和HSV文件导入,设置保存路径和投影坐标〔如图42。图424.6监督分类〔融合方法为HSV,波段为5,4,34.6.1打开GoogleEarth影像作为监督分类的参照使用Envi的GoogleEarthBridge将遥感影像和街区矢量图层导入GoogleEarth,以便在GoogleEarth寻找参考影像时更快定位。在ENVI4.7软件中,选择Spectral->SPEARTools->GoogleEarthBridge,如图43。图43在GoogleEarthBridge窗口中,选择添加的遥感影像E:\大三下\RS&GIS\Data\作业\数据处理\HSV,如图44,点击NEXT进入下一步,设置遥感影像显示的参数,选择5、4、3波段进行彩色显示,如图45,点击NEXT进入下一步,同样添加矢量图层,如图46,保存路径为E:\大三下\RS&GIS\Data\作业\数据处理\栅格数据\GoogleEarthBridge.kml。图44 图45图464.6.2建立兴趣区在ENVI4.7软件中,分别打开HSV.img和重分类后的直接栅格化房屋密度图,并为两个窗口建立其同步,如图47。图47打开兴趣区模版,选择Overlay->RegionofInterest,打开ROITool,如图48。图48在ROITool窗口中,我们可以对ROIName和Color进行编辑,双击ROIName的区域,可以对训练样区进行名字的修改,右击Color可以选择训练样区的颜色,如图49,windows选择Zoom窗口,即只能在Zoom中选择训练样区；在#1Zoom窗口中,我们鼠标左键选择训练样区,右击确定,第二次右击即保存该训练样区,如图50。图49图50在ROITool窗口中,点击NewRegion建立其新的ROI样区,Goto可以跳转各个训练样区,就可以对其进行修改,在实验中,我们将分成11类,进行监督分类,如：water,ransportation,forest,commercial,grass,industrial,cropland,fallow,residential-H,residential-M,residential-L,共11类,如图51。图514.6.3训练样区的选择训练样区的选择依据,如图52,〔aCommercial,〔bindustial,〔ctransportation,〔dhigh-densityresidential,<e>medium-densityresidential,<f>low-densityresidential。图524.6.4训练样区的评价4.6.4.1分离性的定量分析在ROITool窗口中,选择Options->ComputeROISeparability,计算样本的可分离性。如图53,表示各个样本类型之间的可分离性,用Jeffries-Matusita,TransformedDivergence参数表示,这两个参数的值在0~2.0之间,大于1.9说明样本之间可分离性好；小于1.8大于1.4属于合格样本；小于1.4需要重新选择样本；小于1,考虑将两类样本合成一类样本。图53图54从分离性的数据中,我们可以得到,transportation、industial和commercial这三类的分离性PairSeparation都小于1,两两分离性分别为0.66428386、0.81190866和1.27089663,可知这三类的分离效果较差,即使经过多次的训练样区的选择,效果也是一样,因此在实验中,将这三类合并为urban一类。而L_residential和H_residential的分离性只有1.03427042,H_residential和M_residential的分离性只有1.15958190,M_residential和L_residential的分离性只有1.26650863其原因主要是在住宅围干扰因素,且周围的环境类似。4.6.4.2分离度的定性分析在ROITool窗口中,选择File->ExportROIston-DVisualizer,从多维视图进行训练样区分离性的定性分析,如图55,56。图55图56从图54,在n-DControls窗口中,选择1,2,3分别为波段5,4,3,点击Start,在n-DVisualizer窗口中的像元点会自动从从多个方位来给我们展示分离情况。同样我们可以在此窗口中对其分离性进行修改,对像元进行重新归类,在n-DControls窗口中,在Class中选择将归类的颜色,在n-DVisualizer窗口中选择分离效果较低的像元点,或者混合的像元点,右击确定,如图57。图57在ROITool窗口中,选择Option->MergeRegions,将transportation、industial和commercial这三类训练样区进行合并,如图58。图584.6.5执行监督分类选择classification->Supervised->MaximumLikelihood,ENVI给我提供了四种多光谱监督分类的方法,分别为平行六面体、最小距离、马氏距离以及最大似然,在此次实验中,我们选择MaximumLikelihood〔最大似然法进行监督分类,如图59。图59在ClassificationInputFile窗口中建立Mask,选择MaskOptions->BuildMask,如图60。图60在MaskDefinition中,选择Options->ImportEVFs,添加Census_pro_.evf文件,如图61。图61在MaximumLikelihoodParameters窗口中,设置好参数和保存路径,如图62。图624.6.6分类后处理在ENVI4.7窗口中打开视图,选择Overlay->Classification,打开分类结果,并对颜色进行编辑,如图63,64。图63图64选择Classification->PostClassification->Majority/MinorityAnalysis,在ENVI中为我们提供了Majority/MinorityAnalysis、ClumpClasses、SieveClasses三类,主要是用来进行小斑点的去除,此次实验中选择Majority/MinorityAnalysis进行分类后处理,如图65,如图66为经过处理后与未处理的对比,从中我们可以看出有些小斑被周围分类给覆盖了。图65HSV融合影像未处理处理后图66在ENVI中,选择Classification->PostClassification->ConfusionMatrix->UsingGroundTruthROIs,查看混淆矩阵,并获取分类精度,生产者精度和使用者精度以及总精度和Kappa系数,如图67,图68为混淆矩阵。图67图684.6.7评价结果分析本实验先按照书本将LULC类别划为11类,分析可知交通用地、商业用地、工业用地这三类的分离效果较低,生产者精度和使用者精度都不高并且这三类地物的像元十分相近,因此,我直接将这三类合并,按照9类来执行之后的分类,此时的总精度达到93.5656%,Kappa系数为0.9222。表111类与9类精度与Kappa系数比较表11类9类LULCProd.AccUserAcc.LULCProd.AccUserAcc.water99.95%99.95%water99.95%99.95%forest97.08%99.92%forest97.08%99.92%grass96.65%92.46%grass96.65%92.61%cropland95.26%97.64%cropland95.22%97.70%fallow93.00%97.15%fallow94.53%95.50%residential-H89.97%83.25%residential-H90.62%83.43%residential-M58.00%42.43%residential-M62.52%37.72%residential-L85.18%58.25%residential-L84.96%58.27%Transportation43.28%69.49%urban86.20%97.28%commercial90.42%77.63%industrial56.82%52.47%Overallaccuracy90.81%Overallaccuracy93.57%Overallkappa0.8904Overallkappa0.9222进行分类统计,选择Classification->PostClassification->ClassStatistics,在ClassStatisticsResult窗口中,我们可以获取,我们分类的像元数、占该类像元数总数的比例以及该类地物的面积,如图66。图694.6.8分类结果面积统计1将9类的分类结果转成矢量在ENVI中,选择Classification->PostClassification->ClassificationtoVector,如图70,选择输入的图层9LULC1.img。图70在RasterToVectorParameters窗口中,选择9类分类,Output选择输出单个图层SingleLayer,如图71。图71如图72为矢量化结果图72在ENVI格式转成Shp文件,会自动生成面积,单位为m2,首先,我们应对其进行归类处理,在Table中,选择Class_Name->Summarize,在Summarize窗口中,归类的字段为Class_Name,字段Area为求和〔Sum。如图73。图73分类结果如下其中水体面积为21.4594平方千米,占总面积的2.06%,而水体的实际面积为18平方千米。所有地物中,低密度居民区面积百分比最大,达到34.23%；其次是urban,占13.31%。居民区总面积达到55.98%。其他各类地物所占比例见表表分类后各地物面积及其所占百分百地物类型面积/km2所占百分比/%water21.45942.06urban138.852713.31forest102.61269.84grass112.381210.77cropland52.77025.06fallow31.08762.98residential-H103.48569.92residential-M123.370411.83residential-L357.037634.23sum1043.0573100.00制作饼状图如下：图744.6.9分类结果图7511种分类LULC图769种分类LULC4.7分类结果评价与分析4.7.1未加入房屋密度图层的分类结果评价与分析如上图68为未加入房屋密度图层的分类结果,如下表2为未加入房屋密度图层进行分类的分离度。从表2我们可以得到,在分离度上分离度效果较差的为低、中、高密度居民区,其中高密度与中密度居民区的分离度为1.1596,中密度与低密度居民区的分离度为1.2665,低密度与高密度居民区之间的分离度为1.0343。其次的为中密度居民区与城市的分离度效果较低,只有1.0358,通过查看ROI训练样本的选择,发现主要是因为中密度的居民区周围都有道路,从而形成了干扰,故此分离效果较差。通过计算混淆矩阵,总精度为93.5656%,KappaCoefficient=0.9222。表2未加入房屋密度图层进行分类的分离度SignatureName23456789water<1>2.00002.00002.00002.00002.00002.00002.00001.9994forest<2>1.82962.00002.00001.99581.99011.77401.9999grass<3>1.99731.99991.96641.89541.72291.9922cropland<4>1.82201.99991.64751.99951.9339fallow<5>2.00001.79031.99851.8532residential-H<6>1.15961.03431.7053residential-M<7>1.26651.0358residential-L<8>1.8943urban<9>4.7.2加入IDW插值房屋密度图层的分类结果评价与分析如图77、图78和表3,从图77,我们可以直观的了解到,在加入IDW插值房屋密度图层,urban的像元数从617123points->746549points,说明urban的面积在增大,而中、高密度的居民区的面积减少,其他的训练样区的面积变化不大,分析其原因可能是因为在加入IDW插值房屋密度图,主要会影响到低、中、高密度居民区的面积,然而urban的面积急剧增加可能是因为在未加入房屋密度图时,urban与中、高密度居民区的分离效果不好而导致的。从表3中我们可以得到,高密度与中密度居民区的分离度为1.7889,中密度与低密度居民区的分离度为1.3869,低密度与高密度居民区之间的分离度为1.4155。通过与未加入房屋密度图的分离度来看,在加入IDW插值房屋密度图层的分离度有明显的提高。图77未加入房屋密度图与加入IDW插值房屋密度图层像元数比较图78加入IDW插值房屋密度图层的分类结果图表3加入IDW插值房屋密度图层进行分类的分离度SignatureName23456789water<1>2.0002.0002.0002.0002.0002.0002.0001.9994forest<2>2.0002.0002.0001.99891.99241.79051.9999grass<3>1.99901.83961.98091.92641.73531.9923cropland<4>1.82512.0001.84091.99981.9729fallow<5>2.0001.88661.99941.9427residential-H<6>1.78891.41551.8451residential-M<7>1.38691.3268residential-L<8>1.9006urban<9>4.7.3加入直接栅格化房屋密度图层的分类结果评价与分析如图79、图80和表4,从图79未加入房屋密度图与加入直接栅格化房屋密度图层像元数的比较中,得到加入直接栅格化房屋密度图层对分类面积的影像程度较大,主要为低密度居民区,总像元1586843points->1042126points,总面积急剧减少。而中、高密度居民区的面积在增加。表4的高密度与中、低密度居民区的分离达到了1.8694和1.9304,与之前相比有了明显的增加,说明了在加入直接栅格化房屋密度图层可以有效的对不同密度的居民区进行分类。图79未加入房屋密度图与加入直接栅格化房屋密度图层像元数比较图80加入直接栅格化房屋密度图层的分类结果图表4加入直接栅格化房屋密度图层进行分类的分离度SignatureName23456789water<1>2.0002.0002.0002.0002.0002.0002.0001.9999forest<2>1.83972.0002.0001.99991.99731.97271.9999grass<3>1.99942.0001.99941.93941.93871.9969cropland<4>1.84222.0001.95172.0001.9955fallow<5>2.0001.96111.99991.9864residential-H<6>1.86941.93041.8653residential-M<7>1.34561.3829residential-L<8>1.9502urban<9>4.7.4加入重分类后IDW插值房屋密度图层的分类结果评价与分析如图81、图82和表5,从分类面积来看,未加入房屋密度图与加入重分类IDW插值房屋密度图层的面积基本一样。从表5来看,其分离度与未加入房屋密度图差不多,在低、中、高密度居民区的分离效果没有加强,因此,对IDW插值房屋密度图并未能使分类效果增强。图81未加入房屋密度图与加入重分类IDW插值房屋密度图层像元数比较图82加入重分类IDW插值房屋密度图层的分类结果图表5加入重分类的IDW插值房屋密度图层进行分类的分离度SignatureName23456789water<1>2.00002.00002.000002.00002.00002.00001.9994forest<2>1.84352.000001.99951.99271.79421.9999grass<3>1.999101.98081.90181.72671.9927cropland<4>02.00001.86201.99981.9672fallow<5>0000residential-H<6>1.43371.46021.8701residential-M<7>1.31541.1672residential-L<8>1.8975urban<9>4.7.5从总精度与Kappa系数对分类结果进行评价表6分别列举了3中不同方法的分类精度比较,分别为未加入房屋密度图层,加入IDW插值房屋密度图层,加入直接栅格化房屋密度图层。从总精度和Kappa系数来看,未加入房屋密度图层的情况的分类精度最高,总精度为93.57%,Kappa系数为0.9222；加入房屋密度图进行分类的精度有了明显的提高,其中加入IDW插值房屋密度图进行分类的总精度为95.72%,Kappa系数为0.9459,加入直接栅格化房屋密度图层进行分类的总精度为94.81%,Kappa系数为0.9374。从表6中,我们可以得到加入IDW插值的房屋密度图的总精度和Kappa系数都大于加入直接栅格化房屋密度图,从每一类分类来看,两者的Water和Cropland的用户精度和生产者精度一样；forest和urban的用户者精度基本一样,其生产者精度相差较大,说明了其他类别被分为forest和urban在IDW插值和直接栅格化的密度图中的程度不一样,加入直接栅格化房屋密度被错分这两类的像元较多；grass、residential-H和residential-L的生产者精度一致,其中grass和residential-L被错分为其他类别的像元较多,而加入IDW插值房屋密度图的residential-H被错为其他类别的像元相对与加入直接栅格化房屋密度图较少；加入IDW插值房屋密度图的residential-M的生产者精度和用户者精度都高于加入直接栅格化房屋密度图。从生产者精度和用户者精度可以看出加入IDW插值房屋密度图可以得到更好的分类效果。表6不同分类方法的分类精度比较未加入房屋密度图层加入IDW插值房屋密度图加入直接栅格化房屋密度图LULCPAUAPAUAPAUAwater99.95%99.95%99.95%99.95%99.95%99.95%forest97.08%99.92%97.20%99.86%93.54%99.90%grass96.65%92.61%96.14%93.71%97.74%89.97%cropland95.22%97.70%96.17%98.74%97.82%98.26%residential-H90.62%83.43%99.24%94.50%99.35%99.02%residential-M84.96%37.72%86.27%58.04%80.61%52.69%residential-L62.52%58.27%92.26%66.19%94.47%52.98%urban86.20%97.28%91.02%98.26%87.44%98.15%Overallaccuracy93.57%95.72%94.81%Overallkappa0.92220.94590.93744.7.6分类结果总体评价本实验中,通过从总精度、Kappa系数、分类结果图三个角度进行综合评估,个人认为通过加入房屋密度作为选训练样区的依据可以很大程度上来提高分类的精度。residential_H,从图像上看在高密度住房区,IDW插值的是采取中心点位置,从而可以将周围的道路和树木的干扰,直接栅格化则是将整个区域栅格话,忽略了道路和树木的影响,如图83。图832residential_M,经过IDW插值的话,周围的Point会相互影响,从而降低了错分和被误分的概率,而通过直接栅格化则是整个区域的值都被同一化,导致了误差的增加,如图84。HSV融合影像 IDW插值 直接栅格化图84直接栅格化的误差来源直接栅格化将房屋密度都按照街区来显示,这种方案带来的弊端是,一个街区中部分区域的房屋单元数量影响整个街区的房屋密度情况。如图85,红框中的街区,实际上仅有靠街区外围部分房屋单元数量众多,中间大部分为水体,但是在直接删格化的房屋密度图中整个街区都被定为中密度住宅区,从属性值看出该区域的hou_dem1的值为681。HSV融合影像 直接栅格化 属性值图85IDW插值的误差来源虽然IDW插值中Point会相互影影响,但是在一些水域和很大面积都为一类的地物,IDW插值则会产生误判,下图中蓝色本为高密度房屋居民区,却扩展到水域的地区,导致水体部分被分为高密度居民区,如图86。图864.7.7与其他训练样区的分类精度和Kappa系数的计算在实验中,用****的训练样区来验证自己的分类总精度和Kappa系数〔实验中,只对未加入房屋密度图层的HSV进行验证。如表7为****的训练样区的分离度情况。SignatureName2345678water<1>1.99701.99761.99651.99581.99662.00001.9464forest<2>1.52911.93481.81241.60941.34601.7031grass<3>1.70011.36890.92271.09851.1832cropland<4>1.79931.49681.15481.4372residential-H<5>0.48940.86850.4947residential-M<6>0.45270.4671residential-L<7>0.4647urban<8>表7训练样区的分离度选择Classification->PostClassification->ConfusionMatrix->UsingGroundTruthImage〔如图87,选择自己的训练样区和验证训练样区,根据训练样区分类进行对比〔如图88。图87 图88如图89,为通过验证训练样区的总精度和Kappa系数,从图中,可知总精度为45.8993%,Kappa系数为0.3929,相对来说还是比较低的。表8为生产者精度和用户者精度。图89表8LULCPAUAwater66.07%99.66%forest39.88%99.93%grass64.99%46.97%cropland46.36%99.92%residential-H49.76%69.65%residential-M6.79%11.25%residential-L25.40%4.73%urban88.69%71.36%Overallaccuracy45.90%Overallkappa0.39294.8决策树分类4.8.1决策树分类原理遥感信息的提取与分类一直是遥感技术领域研究的一项重要容。遥感分类应用中,监督与非监督分类的传统分类方法以及人工神经网络分类、专家系统分类等新方法都以影像光谱特征为基础。然而由于影像本身存在"同物异谱、异物同谱"现象,这种纯粹依赖地物光谱特征的分类方法往往会出现较多的错分、漏分情况,从而导致分类精度降低。近年来蓬勃发展的决策树方法可以充分利用影像的多种特征来进行分类,因此利用影像的多特征结合决策树在影像分类领域具有很大的应用潜能。决策树是类似于流程图的树型结构,一个决策树由一个根节点,一系列部节点和分支以及若干叶节点组成,每个部节点只有一个父节点和两个或多个子节点,节点和子节点之间形成分支,决策树的每个部节点代表一个测试属性或属性的集合,每个分支对应该属性的每个可能值,叶节点则对应一个类别属性值,即分类结果,不同的叶节点可以对应相同的类别属性值。4.8.2数据预处理预处理过程主要是针对遥感影像进行大气较正,在ENVI4.7给我提供了FLAASH和QuickAtmosphericCorrection,在这两种大气校正方法QuickAtmosphericCorrection的精度要低于FLAASH。但是,在做FLAASH大气校正前需要先进行辐射定标,但是根据老师所提供的数据我们只了解到其传感器为LandsetETM+,因此在辐射定标中缺少数据,故采取精度较低的QuickAtmosphericCorrection方法进行大气校正,在执行QuickAtmosphericCorrection校正前还需要进行中心波长的录入,由于同一传感器的中心波长都是一致的,故百度到了各个波段的中心波长。1中心波长的录入,在AvailableBandsList窗口中,右击Stack,选择EditHeader〔如图90。图902在Header窗口中,点击EditAttributes,选择Wavelengths〔如图91图91通过百度得知,LandsetTM各个波段的中心波长如下图所示,在EditWavelengthsvalues窗口中编辑中心波长〔如图92。图924在ENVI中,选择Spectral->QuickAtmosphericCorrection<如图93>图935在窗口中,SensorType选择LandsetTM,点击OK〔如图94图94如图95所示,左边为大气校正后影像的ZProfile,右边为未校正影像的ZProfile,从ZProfile我们可以看出经过校正后的影像在绿色波段上的value明显高与蓝色和红色波段,在地物上能够有效的区分出植被,有利于计算NDVI。图954.8.3指数的计算1水体指数〔MNDWI在ERDASIMAINE9.2软件中,选择Modeler->ModelMaker,如图96。图96在New_Model窗口中,选择输入输出框和中间的运算函数,第一个函数框的公式为水体指数的提取MNDWI=〔band2-band5/<band2+band5+0.00001>,并输出影像,去除斑点,在View窗口中,我们先大概了解斑点的像元值,如图94,我们可得水体的像元值x小于0,因此确定水体的条件函数为EITHER1IF<$n1_stack_b1<=0>OR0OTHERWISE,输出最后提取的水体信息。图97为输入数据格式 图 98为输入运算函数图99图100在ENVI4.7软件中,选择BasicTools->LayerStacking,将在ERDASIMAINE9.2提取出来的水体信息进行叠加,如图101。图1012植被信息〔NDWI指数的提取在ENVI软件中,选择Transform->NDVI,输入校正后的影像数据〔如图102。图1023IsoData数据的提取IsoData为非监督分类中的一种方法,选择Classification->Unsupervised->IsoData,参数设置为默认值,如图103。图1034.8.4执行决策树1RuleGen插件的使用将如下四个文件分别复制到ENVI安装目录下的save_add文件夹下。2执行决策树在ENVI中,选择Classification->DecisionTree->RuleGen->Classifier,如图104。图104在RuleGen-Classifier窗口中,首先为默认缺省值进行分类,如图105。图105此时,需要重启ENVI,选择Classification->DecisionTre