利用ENVI、eCognition Developer（易康）和Excel和ArcMap软件对南京江北新区统计分析土地利用变化的具体信息

1.研究区概述

南京的江北新区是2015年提出的在南京长江北设置的新区南京城市扩张的新发展极，经过这些年的发展，土地利用类型变化较大。

2.数据源介绍与数据获取过程

Landsat 8卫星，卫星高度705m，携带陆地成像仪（OLI）和热红外传感器（TIRS），OLI陆地成像仪包括9个波段（Band 1 Coastal（海岸波段）Band 2 Blue（蓝波段）Band 3 Green（绿波段）Band 4 Red（红波段）Band 5 NIR（近红外波段）Band 6 SWIR 1（短波红外1）Band 7 SWIR 2（短波红外2）Band 8 Pan（全色波段）Band 9 Cirrus（卷云波段）），空间分辨率为30米，其中包括一个15米的全色波段。

步骤：

数据获取首先登录“地理空间数据云”网站，点击高级检索，在检索框中填入检索信息，这里选择Landsat8卫星的影像数据，根据研究区的位置，填入空间位置信息；选择时间这里选择南京江北新区的设立时间2015年成立，想下载2015—2022年的影像。但是根据实际的影像数据我选择了2014年和2021年的影像数据；云量小于10%。点击检索，检索框下面出现了符合条件的所有数据。

图 1数据下载

将其加载到ENVI中可以更细致的查看下载的影像，确定影像合格，开始对遥感影像进行处理。

图2 ENVI加载图像

3.影像预处理

3.1 2014影像辐射处理

进行辐射定标，使用工具栏“辐射定标”工具，先选择文件。选择影像的可见光-红外组。

图3 辐射定标文件选择

单击确定，进入下一步，在“辐射定标”对话框设置如下参数，因为后续要进行FLASS大气校正，点击面板上的“应用FLASS设置”按钮可以快速自动填充可用于大气校正的参数。设置输出路径，保存定标数据。

图4辐射定标参数设置

辐射定标后影像结果：

图5 辐射定标结果

定标完成后进行大气校正。大气校正模块中需要输入成像区域的平均高程，先从 DEM 中计算成像区域的平均高程。

在菜单栏，选择文件->打开世界数据->高程（GMTED2010），打开ENVI自带高程数据。用作业影像的范围，用其统计研究区的平局高程，得到研究区平均高程值38.76米。

图6得到平均高程数据

接下来正式进行FLAASH大气校正，在工具栏选择“FLAASH”工具。打开对话框，单击 Input Radiance Image 按钮，在 FLAASH Input File 一栏选择辐射定标后的多光谱影像文件，确定。在弹出的 Radiance Scale Factors 对话框，选择 Use single scale factor for all bands，并将 Single scale factor 设为 1。

图7 Radiance Scale Factors设置

回到 FLAASH Atmospheric Correction Model Input Parameters 对话框，选择传感器类型为Landsat-8 OLI，填写卫星高度705km，平均高程0.038km。成像季节根据ENVI帮助里的介绍选择MLS，气溶胶模型根据影像特点选择Urben。

图8  ENVI帮助文档

设置 FLAASH 参数如下：

图9  FLAASH 参数设置界面

点击底部的 Multispectral Setting 按钮，在Multispectral Settings 对话框中，选择 Kaufman-Tanre Aerosol Retrieval→ Defaults→ Over-Land-Retrieval Standard（660：2100nm）

图10 Multispectral Settings设置

点击确定，回到 FLAASH Atmospheric Correction Model Input Parameters 对话框，选择 Advanced Setting，设置采用文件分块处理方式，设置分块文件大小为 100M，点击确定回到 FLAASH 对话框，点击 Apply 按钮进行大气校正。得到大气校正后图像。

图11大气校正影像结果

3.2 影像融合与裁切

先对2014影像的全色波段进行融合。辐射定标波段选全色波段，步骤与多光谱一样。得到全色波段辐射定标结果。

图12 全色波段辐射定标结果

全色影像定标完成后开始进行图像融合。这里使用NNDiffuse变换融合

使用工具栏“NNDiffuse图像融合”工具。在对话框中，填入低分辨率的多光谱图像和高分辨全色影像。

图13 融合设置界面

点击“确定”得到融合影像。

图14   NNDiffuse变换融合方法融合结果

融合完毕后对图像进行裁切，使用矢量数据进行裁切，将江北新区的矢量数据导入ENVI中

图15  导入矢量数据

使用工具栏里的“subset data from ROIS”进行裁剪。选择被裁影像和裁剪要素。

图16裁剪影像

裁剪结果。

图17裁剪结果

3.3   2021年的影像数据预处理

2021年的影像数据预处理的步骤与2014影像的处理步骤相同。得到2021年的研究区裁剪影像。

图18 2021年影像裁剪结果

4. 变化区域获取与变化信息分析

4.1 2014年影像分类

使用eCognition 软件中的面向对象分类；添加影像，对波段进行重新命名，方便观看。

图19 添加影像

点击界面上的“No Data”按钮，打开对话框后，勾选“Use single value for all layers(union)”去除的值为“0”将遥感影像的黑边去掉，后续分割就不会将黑边分割了。

图20 去除影像黑边

影像创建完成后，打开View Setting，设置显示波段为彩红外显示。再打开处理界面在“process tree”界面的空白处，右键，选择 Append New,在出现添加处理的界面，在Algorithm 里选择“multisolution segmentation”， 在 Level Name 里输入分割层的名字，这里分割尺度准备设成50，就命名为 L50，设置的参数界面。

图21 设置分割处理界面

点击“Execute”执行，得到一个名为 L50的分割结果.

图22 分割结果

可以看到分割还是比较细碎，进一步区域合并得到更大尺度的分割结果。再添加处理的界面，在L50的基础上继续合并，所以“Domain”这里应该选择对象层，同时这里选择“Create above”尺度设为150，给分割层起名为L150，尺度权值比较大时，形状权值设的大一些，设为0.5 。

图23 设置分割处理界面

分割结果。

图24 分割结果

分割结果还较好，接下来进行特征计算，在 Feature View 里，点击打开界面可以看到有非常多的特征，我们主要用的对象特征，波段数比较多的时候，每个波段的对象灰度均值特征也需要创建，依次将八个波段的特征都创建出来。

图25 每个波段的对象灰度均值特征创建结果

点击用户自定义特征，在出现的界面里输入特征的四则运算，计算 NDVI。

图26 计算 NDVI

点击确定后，建立了一个新特征NDVI，双击这个新建的特征，可以查看特征效果。

图27 计算NDVI结果

建立NDWI的新特征。

图28 计算 NDWI

点击确定后，建立了一个新特征NDWI，双击这个新建的特征，可以查看特征效果。

图29 计算NDWI结果

打开样本管理工具，在窗口上右击，创建类别，创建4个类别：植被（zhibei）、水体（shuiti）、裸地（luodi）、建成区（jianchengqu）。

图30 创建分类结果

点击菜单栏“classification”-> “samples”->“select sample”，打开样本选择模式对每类地物进行样本选择，选择结果如下。

图31 样本选择

为所有别分配一个最近邻（NN）分类器，分类完毕后，点击任意一个类别，可以看到最近邻分类器的存在。双击选择分类使用的特征，这里我们选择各波段特征、NDVI、NDWI、长宽比、矩形度、椭圆度、紧致性等特征。

图32 分配一个最近邻（NN）分类器，选择分类使用的特征

然后进行分类过程，在“process tree”右击添加一个“classification”处理，分类层选择 L150，在 Active Classes 里选择所有类别。

图33 在 Active Classes 里选择所有类别

分类结果如下图。

图34 分类结果

完成分类后利用样本进行精度评价，一般用于精度评价的样本应该比较多，所以重新打开样本选择模式，为各类添加增多的样本。

图35 添加分类样本

然后从样本中创建一个 TTA。

图36 创建TTA结果

点击 tool –>Accuracy Assessment，打开精度评价工具，设置参数如下所示：

图37 精度评价界面设置

精度报告如下。

图38 精度评价报告

导出分类结果为矢量数据。

4.2 2021年影像分类

加载影像数据的步骤与2014年影像相同，分类直接使用150为尺度进行分割。

图39 2021年影像设置分割处理界面

分割结果。

图40 2021年影像分割结果

进行特征计算，依次将八个波段的特征都创建出来，创建 ndvi 特征， 双击这个新建的特征，可以查看特征效果。

图41 2021年影像NDVI计算结果

建立NDWI的新特征，双击这个新建的特征，可以查看特征效果。

图42 2021年影像NDWI计算结果

创建类别，为每个类别选择样本。

图43 2021年影像样本添加结果

与2014影像的步骤一样，为所有别分配一个最近邻（NN）分类器，再选择分类使用的特征，选择的特征与2014影像相同。然后进行分类过程，在“process tree”右击添加一个“classification”处理，分类层选择 L150，在 Active Classes 里选择所有类别。进行分类，得到分类结果。

图44 分类结果

完成分类后利用样本进行精度评价，一般用于精度评价的样本应该比较多，所以重新打开样本选择模式，为各类添加增多的样本。

图45 添加样本数量

然后从样本中创建一个 TTA，创建TAT结果：

图46 创建TAT

点击 tool –>Accuracy Assessment，打开精度评价工具，设置参数，得到精度报告。

图47 2021年影像分类精度报告

4.3 变化信息提取分析与成图

导出分类结果的矢量数据。与2014影像的分类结果矢量数据一起导入ARCMAP（易康导出的矢量数据没有分类字段，可能是我操作的不当，这里导出的栅格数据在arcmap中转换为矢量）根据分类字段对分类矢量数据进行“融合”处理。

图48 影像分类矢量结果融合结果

根据融合的结果在属性表中添加字段，进行计算几何的操作，计算出各类地物的面积。将计算后的属性表转换为excel表格，在excel表格中对两年的同一地物类型进行图表对比。

图49导出excel表格结果

图50 南京江北新区用地类型变化图

图51 南京江北新区各类用地类型面积变化图

在ARCMAP中对两年矢量数据的融合结果再进行“相交”操作，将相交后的属性表导出为excel表，在excel表中建立数据透视表得到土地利用类型转移矩阵（单位km2）。

图52 南京江北新区土地利用变化矩阵

回到ARCMAP中依据相加结果，制作重要变化地类的矢量结果图。

图53 南京江北新区2014-2021年地物变化图

5.总结

完成分类后对分类结果进行分析，可以看到分类结果建成区和水体的变化都还比较合理，但是裸地和植被的变化我觉得有些太剧烈，通过地物变化图和变化矩阵，我发现裸地变为植被的比例很高。通过对原始影像的观察，我觉得原因是，虽然南京是常绿植被带，但是其中2014年的影像选用的时间为3月，这个时间段是农田的农作物轮种的时间段，有的农田的农作物被收获了或者新种的农作物还没长出，这造成出现了大量裸地，而2021年的影像上农田里是有农作物的。此次的操作结果不能完全意义上反映江北新区的严格意义上的地类变换。可以选择两年的相同月份遥感影像再进行操作。