土地利用分类详细教程——以高分一号影像为例（上）

1 前言

相信跟我一样刚刚接触遥感的计算机人士来讲，突然转行开始弄遥感，刚开始肯定会一头雾水。更别说什么土地利用分类的制作了。在这里，我将会用最为直接的图形流程操作来给像我一样初识遥感的童鞋们讲解一下制作土地利用的整个过程。

首先，我先用一幅流程图来概括一下制作土地利用的整个过程。如下图1所示：



图1土地利用分类流程图

这个流程操作适用于一般影像制作土地利用过程。然而对于高分影像数据而言，制作土地利用流程需要进行一些修改。在文章的最后，我会将高分影像数据制作土地利用的流程图附上。现在，让我们开始高分影像数据土地利用制作的过程。

注：在教程中，我有些地方使用了envi5.3。如若想获取envi5.3的请在评论中留下邮箱，届时我将发送链接分享。该软件分享仅作为科研使用，如若进行商用，请使用正版，否则后果自负。

2 裁剪

先对影像进行裁剪，这里我对影像进行裁剪的目的有两个：

1、减少接下来工作量

2、这样可以选择较高分辨率的DEM（DEM可能不能覆盖影像，这里我通过shp图来进行裁剪影像，便可以确定使得下载的DEM更小些）

裁剪步骤如下图所示：



图 2 裁剪工具



图3 选择裁剪影像



图4 结果输出

这样就可以将影像裁剪出来了，展示一下我裁剪的影像图：



图5 裁剪影像图

同理，将研究区所有的影像进行裁剪（直接对影像拼接会出现问题，具体说明看正射校正部分）。

3 正射校正

这里我们之所以将正射校正放在拼接前边，是因为我们正射校正的结果需要作为大气校正的初始文件。这里需要说明其实我们刚开始说的流程图也可以先进行正射校正，也就是说正射校正和大气校正之间是没有依赖关系的。高分影像数据的制作可以算做一个特例。

这里还需要进行说明一点：正射校正之前可能无法进行图像镶嵌（拼接）的，无法进行操作会弹出如下图2所示的说明（其实我们影像并没有问题，下图3-5显示两幅影像都是一样的）。我们之所以先进行拼接是因为可能研究区有多幅影像组成，那么如果正射校正的话需要一幅一幅影像的进行校正，然后再进行拼接。这样将会花费大量的时间。因此，我们往往会想到先进行影像的拼接，然后对拼接后的一幅影像进行校正。如果你跟我一样未果的话，那么可以先进行校正，然后再进行拼接操作。



图6 拼接失败提示



图7 影像数据信息



图8 影像坐标信息



图9 影像光谱波段信息

这个问题目前我还没有想到为什么。这里我们就先不纠结这个问题为什么会出现了，那就先用麻烦点的操作进行，先进行正射校正。

这里由于我没有控制点，没有参考的栅格影像，于是采用了RPC文件进行正射校正。如果你有现成的栅格影像或者控制点信息，你可以采用RPC Orthorectification Using Reference Image 工具（envi5.3版本，这里因为我使用的envi5.1，因此放弃了这种方式，不过这种方式要比我所做的方式正射效果优秀，建议使用这个种正射校正方式）。具体操作步骤参考博文：自动采集控制点的RPC正射校正工具 。我采用的方式是使用RPC Orthorectification Workflow工具，具体过程如下图所示：



图10 选择RPC Orthorectification Workflow工具



图11 文件选择

这里DEM File默认是2010年的900m分辨率的DEM，这里根据实际情况，如果有更高分辨率的DEM的话，选择高分辨率的DEM效果更好。这里我已经对影像进行了裁剪，因此影像范围肯定小于DEM范围了。这里我选择的是30m分辨率的DEM。



图12 DEM跟裁剪影像图

中间白色边框内是两幅未拼接的裁剪影像图。接下来进行选择裁剪影像和DEM进行文件选择。



图13 裁剪文件选择

点击“Next”。接下进行RPC设置。



图 14 RPC设置

这里因为我没有控制点，因此对于GPC的设置我没有进行处理。直接进行接下来选项的设置。



图 15 Advanced选项设置

Advanced设置，其中辐射校正算法中分别是：最近邻法、线性内插、三次卷积法，这里我们可以选择三次卷积方法。



图 16 Statistics选项设置

Statistics设置我们这里选择All。



图17 辐射校正输出选项

（这里明天进行一下对比验证）

4 辐射定标和大气校正

在开始之前，我这里先说明一下辐射定标和大气校正。同样，这里引用上边所引用的文章几何校正；正射校正；几何配准；影像配准，空间配准；辐射定标；大气校正；辐射校正的概念解释一下各定义。

4.1 辐射定标说明

正如很多资料所说，不同学者解释辐射定标的定义不同：

①定标是将传感器所得的测量值变换为绝对亮度为与地表反射率、表面温度等物理量有关的相对值的处理过程（赵英时等《遥感应用分析原理与方法》）

②遥感器定标就是建立遥感传感器的数字量化输出值DN与其所对应的视场中的辐射亮度值之间的定量关系（陈述彭，也是度娘的解释）

③辐射定标是将传感器记录的电压或数字值转换成绝对辐射亮度的过程（梁顺林《定量遥感》）

④辐射定标就是将记录的原始DN值转换为大气外层表面反射率，目的是消除传感器本身产生的误差

辐射定标有很多方法：实验室定标、星上定标、场地定标。具体参见ESRI中国社区ENVI版块，dsbin传感器定标http://bbs.esrichina-bj.cn/ESRI/viewthread.php?tid=56191

另特别说一下：将DN值转化为辐射亮度和大气表观反射率的公式：

1.

将初始DN值转化为辐射亮度，其中Lb是辐射亮度值，单位是： W/cm2.μm.sr（瓦特/平方厘米.微米.球面度），Gain和Bias是增益和偏差，单位和辐射亮度值相同，可以看出，辐射亮度和DN值是线性关系。

2.

将辐射亮度值转换为大气表观反射率，式中Lλ为辐射亮度值，d为天文单位的日地距离，ESUNλ为太阳表观辐射率均值，θs是以度为单位的太阳高度角。不过总的来说，这部分的工作基本上不需要用户自己做，相关的系数都包含在数据的头文件或者元数据中了。

4.2 大气校正说明

大气校正的目的是消除大气和光照等因素对地物反射的影响，广义上讲获得地物反射率、辐射率或者地表温度等真实物理模型参数；狭义上是获取地物真实反射率数据。用来消除大气中水蒸气、氧气、二氧化碳、甲烷和臭氧等物质对地物反射的影响，消除大气分子和气溶胶散射的影响。大多数情况下，大气校正同时也是反演地物真实反射率的过程。

什么情况下需要做大气校正，我们购买或者其他途径获取的影像是否做过大气校正。通俗来讲，如果我们需要定量反演或者获取地球信息、精确识别地物等，需要使用影像上真实反映对太阳光的辐射情况，那么需要做大气校正。一般购买的影像，说明文档会注明经过辐射校正，其实这个辐射校正是指的粗的辐射校正，只是做了系统大气校正，就跟系统几何校正的意义是一样的。

大气校正跟辐射校正是相互联系的，在envi中做大气校正前提就是需要辐射定标后的结果，实际过程就是将影像的辐射亮度值，最终转换为地表反射率。

4.3 辐射校正说明

是指对由于外界因素，数据获取和传输系统产生的系统的、随机的辐射失真或畸变进行的校正，消除或改正因辐射误差而引起影像畸变的过程。（包括辐射定标和大气校正）三者关系可以归结为下图：



图 18 辐射校正、辐射定标、大气校正关系图

Ok，这里感谢一下“倒斗我不会”给予的解释。接下来我们继续我们的教程——辐射校正。

4.4 辐射定标

接下来我们开始进行辐射校正。首先先进行辐射定标。在envi中使用Radiometric Correction—>>Radiometric Calibration（辐射校正）。具体如下图所示：



图 19 辐射定标

当我选择所选影像进行下一步时会出现如下图所示警示：



图20 警示

？？？这是什么意思呢？经过查阅相关文献，我从中得以将此问题解决。实际上辐射定标的目的就是将DN值转换为辐射亮度值（前边已经提及到了）。因此需要将传感器的Gain值和Offset值输入到envi中进行转换成辐射亮度值。envi中通过编辑envi 头文件来进行编辑Gain和Offset值。Envi 5.3中（具体获取方式参照前言）通过Raster Management中的Edit ENVI Header工具进行对两个值修改，具体操作如下：



图21 Set Raster Metadata

通过点击Add可以添加Gain和Offset工具:



图22 Add Gain and Offset工具

添加完后，我们可以发现在界面中出现了对Gain和Offset值修改的界面。如下图所示：



图23 Gain and Offset界面

接下来我们就开始对两个值进行添加，这里其实还有一种修改方式（手动），具体的操作方式，不在叙述，具体操作请参照博客ENVI下高分一号PMS相机多光谱数据大气校正，在这里我采用自动选择文件的方式进行快速添加。首先，我在这里列一下2015年国产陆地观测卫星外场绝对辐射了定标系数表供手动方式使用：



图24 2015各国产卫星辐射定标系数

在接下来就是对两个值的修改，这里我们以修改Gain值为例，影像的波谱函数等参数获取，可以通过中国资源卫星应用中心获取。



图25 添加Gain文件

然后我们可以发现各波段的Gain值已经输入到了界面中，如下：



图26 添加Gain值

同理，再修改Offset值。通过View Metadata可以来确认一下各波段是否已经赋值给了Gain和Offset。



图27 检查各波段参数值

通过检查光谱信息我们发现，影像需要大气校正还缺条件，大气校正需要以下几个条件：

1、数据是经过定标后的辐射亮度（辐射率）数据（uW/cm2*nm*sr）

2、数据有中心波长（wavelenth）值，如果是高光谱还需要波段宽度（FWHM）。

3、数据类型支持四种：浮点型、长整型、整型、无符号整型

4、波谱范围：400~2500nm

可以发现我们的影像中并没有中心波长，因此我们还需要进行输入中心波长信息。接下来我们进行中心波长的添加，对于中心波长的添加这里我采用手动方式进行添加。具体的方式我可以通过使用envi打开影像的光谱曲线来进行确定。具体步骤如下：

使用Spectral Library进行对光谱曲线文件的添加。



图28 Spectral Library Viewer



图29 添加光谱曲线函数文件

通过点击各波段最高点可以获取中心波长（Data Value=1处，需要注意这里的单位是nm）。



同理获取其它波段的波长，这里进行一下高分一号各波段中心波长的整理：

中心波长（um）

band1： 0.485

band2： 0.56

band3： 0.696

band4： 0.797

这样我们就完成了大气校正的条件，接下来可以进行辐射定标的其它操作了。需要注意的是我需要对辐射定标的文件的格式进行限制一下，因为大气校正所识别的文件格式为BIL或者BIP格式，因此需要进行选择。



图30 文件格式

（随后修改：其实这里只是对影像文件进行了头文件的参数编辑，因此文件格式这里也可以不进行选择）

接下来时使用Radiometric Correction > Radiometric Calibration工具，进行辐射定标，具体步骤如图所19，选择文件后发现上述的错误没了，弹出辐射校正对话框，如下图所示：



图31 辐射定标输出

这样，辐射定标操作就完成了！

4.5 大气校正–FLAASH工具

辐射定标完成后，我们将开始大气校正。辐射定标的完成成功与否直接影响着大气校正的结果。这里我们使用FLAASH Atmospheric Correction工具对影像进行大气校正。如下图：



图32 FLAASH Atmospheric Correction工具

打开后弹出FLAASH模型工具的输入参数。首先我们需要编辑Input Radiance Image，该文件就是我们上边刚刚处理完成的辐射定标文件。选择文件后弹出Radiance Scale Factors，这里我们选择Use single那个选项，scale factor这里我们已经在辐射定标那里进行了转换，这里不需要再进行修改了。如下图所示：



图33 Radiance Scale Factors

接下来两个是设置大气校正的输出文件路径



图34 设置大气校正输出文件路径

然后我们再设置传感器的各个参数，Sensor Type选择GF-1，然后高分一号的参数，如传感器高度、影像分辨率都会自动填充上去，这里有个地面高程的需要我们根据研究区来进行获取。如果我们手中并没有确切的平均高程数据，在这里我们可以借助envi提供的高程影像进行获取研究区域的平均高程。具体方式如下：

打开全球高程数据（2010年）



图35 GMTED2010

然后借助Compute Statistics工具进行获取



图36 Compute Statistics工具

选择高程数据



图37 选择高程数据文件

按Stats Subset按钮，在弹出的Select Statistics Subset对话框中，点击File按钮选择研究区。



图38 选择研究区文件

弹出计算统计参数对话框，默认是基本统计



图39 Basics Stats

点击OK,将弹出基本统计的结果，这里我们看到我们想要的平均高程。



图40 Mean Elevation

这样我们就获取了我们所研究区域的平均高程是135.592195m。在大气校正填写是注意单位，需要转换为km。接下来我们需要设置Flight Data参数。

对于Flight Data影像获取时间，可以通过View Metadata来查看。不幸的是我的影像啥都没有，那怎么获取影像获取时间呢？



图41 Time

其实我们的影像数据中往往会文件进行记录数据，其中有个xml的文件，学过计算机的人都知道，xml文件主要用于数据存储传输用的，打开它，你会惊喜的发现，你所要的东西就在里面。如下图：



图42 Receive Time

这里需要注意的是，我们这里填写的是格林尼治时间，而不是北京时间。如若我们想要获取北京时间的话，需要再此基础上加8个小时。



图43 time setting

选择Atmospheric Model时，参考help给予的提示。



图44 Select Atmospheric Model

因为我们没有水汽的含量值，因此我们采用第二种方式进行选择。通过第二个表我们可以看到可以根据两个变量进行判断：一个是纬度；另一个是月份。而且我们发现该表90度北纬度没有，12月份也没有，因此通过此规律可以断定两者都是“就小原则”，这里我的研究区是35.7N~36.6N，时间是4月23，因此选择30N：March，因此模型选择MLS。

接下来需要设置气溶胶模型设置，至于做不做气溶胶，我们可以通过我们的影响是否可以做气溶胶来确定，如果KT算法设置总Defaults三个中上行下行均不能都有值，那么说明影像不能做气溶胶的反演。



图45 气溶胶KT算法

需要设置的是Advanced选项，这里我们根据自己电脑配置进行选择。



图46 Advanced设置

然后点击Apply按钮，完成大气校正。大气校正成功与否，我们可以查看一下植被的光谱曲线是否标准。

这里我粗略的画了一下常见几种地物的光谱曲线图，供大家参考。



图47 光谱曲线图

大气校正计算时间比较常，需要耐心等待，而且校正完的影响也不会自动加载到View中，需要打开数据管理进行加载。这里我们选择Load CIR假彩色加载(右击)对照一下大气校正前后的光谱曲线图，发现大气校正成功。

最后说明一点，前边好多操作其实可以省略的包括中心波长的编辑。如果你的版本是envi5.3的话直接打开xml文件就可以进行接下来的操作。由于刚开始我使用的envi5.1版本，因此手动编辑的内容要多一些。