Contents

16.2. 使用 SAGA 地理信息系统处理多光谱卫星图像

16.2.1. 方法

在本节中，我们将介绍如何在QGIS中使用 SAGA 地理信息系统根据高分辨率Landsat 8图像构建土地利用地图。 本部分的研究区域位于塞尔-庞松湖周围（图4.5）。 这个湖位于法国阿尔卑斯山的杜兰斯河上（北纬44°31'47';北纬6°23'35'E）。

_images/image116_x89.jpeg

图4.5.法国阿尔卑斯山塞尔-庞松湖的位置

我们将使用决策树来构建面积层的土地利用。使用这种技术，用户可以精确地定义具体的分类规则。 该决策树将基于原始陆地卫星频段的特定阈值以及土壤调整植被指数（SAVI）或自动水提取指数（AWEI）等衍生指数。 我们将遵循的步骤总结如下：

  • data acquisition;

  • radiometric corrections of the satellite images;

  • cropping of the satellite images according to the extent of the study area;

  • building of color composites and exploration of the study area;

  • calculation of the indices SAVI and AWEI;

  • definition of the decision tree;

  • implementation of the decision tree;

  • filtering of the subsequent classification;

  • simplification of the land use by reclassification of the pixels.

16.2.2. 已使用数据的获取和呈现

在本练习中，我们使用2016年8月23日拍摄的Landsat 8多光谱图像。 Landsat 8是始于20世纪70年代初的美国Landsat计划的最新卫星。 该传感器自2013年2月起进入轨道，平均每16天拍摄同一陆地部分的图像。 它是一个多光谱传感器，这意味着它能够记录11个不同光谱带中表面反射的能量（表4.1）。 空间分辨率根据15 m到100 m的频谱带而异。陆地卫星图像可从许多网站免费获取。 最方便的门户之一是由 美国地质调查局 管理的： EarthExplorer（ https://earthexplorer.usgs.gov/ ）。 一旦注册，就可以从美国宇航局甚至欧洲航天局下载所有Landsat档案和其他产品。

Spectral band

Wavelength (µm)

Domain

Resolution (m)

1

0.433–0.453

Deep blue

30

2

0.450–0.515

Blue

30

3

0.525–0.600

Green

30

4

0.630–0.680

Red

30

5

0.845–0.885

Near IR

30

6

1.560–1.660

Middle IR 1

30

7

2.100–2.300

Middle IR 2

30

8

0.500–0.680

Panchromatic

15

9

1.360–1.390

IR

30

10

10.6–11.2

Thermal IR 1

100

11

11.5–12.5

Thermal IR 2

100

表4.1. Landsat 8频段特征

除了Landsat频段外，我们还提供了一个具有研究区域范围的形状文件。它仅包含一个矩形形状。 我们将使用它来裁剪Landsat频段，以便仅在我们的研究区域工作。

16.2.3. 卫星图像的修正

在处理卫星图像之前，需要进行一些预处理。第一步是根据投影系统对图像进行地理参考， 并使用与控制点相结合的数字模型对其进行正射纠正。幸运的是，大多数时候， 这一步已经在我们从不同门户下载的图像上完成了。 在我们的案例中，我们使用的投影系统是UTM区31 N系统（EPSG代码32631）。 我们将在本次练习中保留该系统。

第二步是将传感器记录的数字值转换为大气顶部发射的反射率。 在此步骤之后，图像的每个像素的值对应于“大气顶部”（OTA）反射率[OSE 16]。 这种方法可以消除从大气到传感器的过程中可能添加到反射率中的“噪音”。 这种修正通常被称为“辐射修正”。

第三步，称为“大气修正”，包括消除大气的“噪音”，以获得表面反射率。 这种噪音可能是由于大气湿度或气溶胶造成的。有很多方法可以执行这种修正。 然而，在我们的情况下，此更正不是强制性的，因为我们只处理一个日期和一个小区域[SON 01]。 事实上，通过应用此修正，存在进一步增加图像中噪音的风险。 因此，在我们的例子中，我们不会应用大气修正。

现在我们将了解如何执行我们唯一需要的修正，即热辐射学修正，以获得“大气顶部”反射率（TOA）。 我们可以手动或自动执行此更正。在这里，我们将使用QGIS中的补充模块来执行。

我们使用的模块称为“半自动分类插件”。该模块主要用于对卫星图像进行分类，但也可以进行一些预处理。 要安装此模块，请转至“插件-管理并安装插件”。通过在搜索栏中输入插件名称来搜索插件， 然后单击“安装”。

_images/image117_xb6.png

安装模块“半自动分类插件：SCP”

安装此模块后，名为“SCP”的新菜单即可使用。我们现在可以对地球资源卫星图像进行辐射修正。 要执行此操作，请转至菜单“SCP '预处理' Landsat”。在新窗口中，转到“包含Landsat频段的目录”行， 然后搜索包含Landsat频段的目录。此目录还应包含与Landsat图像关联的元数据文件（.MTL）。 此元数据文件是执行辐射修正所必需的。 设置此路径后，模块将自行检测（由于元数据文件）与图像相关且需要纠正的数据 （传感器、采集数据、光谱带等）。

_images/image118_xdt.png

SCP模块中的辐射修正设置。

有关该图的彩色版本，请参阅 www.iste.co.uk/baghdadi/qgis1.zip

虽然我们在这里不会使用它，但在这个阶段，可以将热红外带转换为摄氏度， 通过NOS 1方法执行大气修正并执行全景细化。我们可以保留默认参数。 要运行该过程，请单击窗口右下方的黄色和绿色按钮。我们必须定义要存储更正的频段的目录。 我们可以称之为“L8_2016-08- 23_DOE”。此更正可能需要几分钟。

执行纠正后，新图像就会出现在QGIS的主窗口中。我们可以使用“识别特征”工具检查一些像素的值。 在纠正后的图像上，像素值的范围从0到1。这些值对应于反射率值。 给定光谱带中反射率值为1的像素表明该像素反射了该光谱带中所有进入的能量。 该像素将是白色的。相反，反射率值为0的像素不反射任何东西。该像素将是黑色的。

16.2.4. 根据研究区域裁剪图像

我们现在将根据提供的形状文件“study_area_32631.shp”裁剪Landsat图像， 以便仅关注我们的研究区域。首先，我们在QGIS中加载该层。 该层已经在我们正在使用的投影系统中（UTM 31 N-EPSG 32631）。 我们组织我们的层，将研究区域的层带到顶部，然后排列Landsat频段，第一个频段在顶部， 最后一个频段在底部（图4.8）。

_images/image119_xk9.jpeg

图4.8. QGIS中不同层的组织

为了裁剪我们的图像，我们使用了名为“带多边形的裁剪栅格”的 SAGA 模块。 该模块位于“处理工具箱”中的菜单“佐贺”和“矢量<->栅格”中。 我们也可以通过在“处理工具箱”的搜索栏中输入它的名称来找到它。 此时将出现一个名为“Clip Rster with Polygons”的新窗口。 由于我们想要裁剪具有相同范围的多个栅格，因此我们选择“Run as Batch Process…”选项。 在新窗口中，我们设置了所有要裁剪的栅格和要用来裁剪的图形文件。 让我们从点击图标“…”开始列“输入”，然后点击“从打开的层中选择”。 在“多选”窗口中，我们选择要剪裁的图像。我们只想剪辑前七个乐队，我们不会与剩下的四个乐队合作。 在“多边形”一栏中，我们选择要用来剪裁陆地卫星波段的shapefile。 因为我们对所有的波段都使用相同的多边形，所以我们在每行设置了相同的文件。 在“已裁剪”一栏中，我们设置要保存已裁剪的频带的路径。 我们将它们存储在名为“L8_2016-08-23_TOA”的目录中， 我们将这些频段称为“BX_L8_2016-08-23_TOA_SP”，用当前频段的编号替换“X”。 “Sp”用于“Serre-Ponçon”。每次，选择“不自动填充”选项。我们得到了一个如图4.9所示的窗口。 我们可以直接将生成的图像加载到QGIS中。点击“Run”(运行)启动该过程。 一旦治疗完成，就会出现一个窗口通知我们。

_images/image120_x1d.jpeg

图4.9.使用“用多边形剪裁栅格”模块剪裁Landsat频段

处理完成后，修剪层就会出现在QGIS中。为了更好地组织我们的练习，我们关闭了所有层并导入剪辑的图像。 我们将第一个乐队放在顶部对乐队进行排序。我们现在拥有了前七个陆地卫星频段， 它们以大气顶部的反射率表示，并根据研究区域进行剪辑（图4.10）。

_images/image121_x9f.jpeg

图4.10. Landsat频段根据研究区域（44°37'N，5°57'E）进行剪辑

16.2.5. 色彩合成和该地区的探索

为了制作颜色合成，我们使用GDAL 工具构建虚拟栅格，将Landsat频段合并在一起。 为此，请转至“栅格'其他'构建虚拟栅格（目录）.”。在新窗口中，我们勾选“使用可见的栅格层进行输入”框。 在“输出文件”行设置存储结果文件的路径，我们将其称为“L8_2016-08-23_TOA_stack”。 我们勾选“September”框（图4.11），以便能够单独管理虚拟栅格的每个频段。

_images/image122_x9d.jpeg

图4.11.使用GDAL 构建虚拟栅格。

我们可以在窗口底部查看（并编辑）GDAL 语法

该过程完成后，虚拟栅格直接出现在QGIS中。我们可以使用文本编辑器打开这个虚拟栅格。 它只是一个指向每个单独乐队的ML文件。使用这个虚拟栅格，我们构建颜色合成物。 为了突出该地区的植被，我们构建了一个假彩色合成体。我们将绿色带涂成蓝色，红色带涂成绿色， 近红外带（植被反射最多的一个）涂成红色。 右键单击虚拟栅格图层，转至“Properties”和“Style”。对于“渲染样式”，选择“多频段颜色”。 对于红色频段，我们选择频段5（近红外），对于绿色频段，我们选择频段4（红色），对于蓝色频段， 我们选择频段3（绿色）。我们通过选择“拉伸到最小最大值”选项来提高对比度， 然后单击“加载最小/最大值”面板上的“加载”，以获取每个频段的极端值（图4.12）。

_images/image123_x01.png

图4.12.从虚拟栅格构建颜色合成

我们获得以下颜色合成（图4.13）。在这个颜色合成物上，植被根据其状态或多或少地呈现出强烈的红色。 事实上，“健康状况良好”的植被会反射很多红外线。像素越红，包含的植被就越多。

_images/image124_xik.png

图4.13.伪彩色和近红外红色的复合色

有关该图形的彩色版本，请访问见 www.iste.co.uk/baghdadi/qgis1.zip

我们的研究区域是一个山区，山脉西北侧翼的黑暗区与阴影相对应。 不同频段的反射率受到这种现象的影响。与平坦地区相比，此类地区的土地使用分类更难执行。 然而，一旦我们处理将多个乐队组合在一起的指数，我们就会看到这个伪影往往会消失。

16.2.6. 计算指数以提取植被和水面

在这部分中，我们将计算两个指标。我们首先计算SAVI，这对于提取不同类别的植被很有用。 然后我们将计算AWEI以提取水面。

使用SAVI提取植被

最著名、最广泛使用的指数之一是归一化植被指数。 然而，它并不是唯一的指数，而且确实存在许多其他指数，例如SAVI。 该指数考虑了通过树冠不同程度可见的土壤的反射率。SAVI可以使用以下公式计算：

(1+L)(NIR-R)
SAVI = [4.1]
NIR+R+L

其中近红外带，R是红色带，L是考虑土壤亮度的因子。根据文献，L的值为0.5适合大多数情况。

SAVI的计算可以直接在 SAGA 中使用模式“植被指数（基于斜坡）”进行。 我们在“Processing Toolbox”的 SAGA 菜单的子菜单“图像分析”中找到了这个模块。 当我们执行此模块时，会出现一个名为“植被指数（基于斜坡）”的新窗口。 我们首先需要设置与近红外和红色相对应的频段。在我们的例子中（Landsat 8）， 这些频段分别是频段5和4（图4.14）。然后我们选择要计算的索引以及要存储它们的目录。 在这里，我们仅计算SAVI并将其存储在目录“L8_201608- 23_TOM”中， 将其命名为“savi_L8_2016-08-23”。

_images/image125_xci.png

图4.14.使用 SAGA 模块“植被指数（基于斜坡）”计算SAVI

QGIS中出现与SAVI对应的新层。我们可以通过转到该层的属性和“样式”面板来更改该层的样式。 作为“渲染类型”，我们选择“单带伪色”，对于“颜色”，我们选择“绿色”查找表（图4.15）。

_images/image126_xvh.jpeg

图4.15.设置SAVI层的样式

然后单击“分类”和“应用”并获取SAVI层的色调（图4.16）。

_images/image127_x9r.jpeg

图4.16.根据2016年8月23日的Landsat 8图像计算出Serre-Ponçon地区的SAVI。 有关该图的彩色版本，请参阅 www.iste.co.uk/baghdadi/qgis1.zip 。

在这张图像上（图4.16），深绿色的像素是接近1（最大值）的像素，白色的像素是SAVI值非常低的像素， 接近0甚至负值。SAVI值高的像素是植被茂密的像素。SAVI值较低的像素是没有植被的像素。 最后，不同的绿色深浅，即从0到1的不同值，代表不同的植被密度。

使用AWEI提取水面

我们可以使用许多不同的指数来从多光谱卫星图像中提取水面。 AWEI将在本练习中使用，因为它强大且使用方便 [FEY 14] 。 它是基于蓝色、绿色、近红外和中红外频段的经验指数。 使用AWEI，可以非常方便地将水像素从陆地像素中数字化。 事实上，AWEI旨在对值大于0的水像素和值小于0的陆地像素进行编码。该索引存在两个版本。 第一个可以用于没有阴影的图像，第二个可以用于有阴影的图像。以下公式适用于第二版本[FEY 14]：

AWEIℎ = B + 2.5 ×G − 1.5 × (NIR +MIR1) − 0.25 ×MIR2 [4.2]

AWEIsh是用于具有阴影区域的图像的AWEI版本。 B对应于蓝色频段（Landsat 8的B2），G对应于绿色频段（B3），近红外频段（B5）， MIR 1和MIR 2对应于中红外频段（B6和B7）。

在QGIS或 SAGA GIS中无法自动计算该指数。 我们需要使用 Raster → Raster calculator 中的QGIS的栅格计算器进行计算。 应考虑正确的频段作为输入，手动编写公式。 我们将生成的层存储在工作目录中，并将其命名为“awei_L8_2016-08-23”（图4.17）。

_images/image128_xba.jpeg

图4.17.使用栅格计算器功能计算AWEI

生成的层直接出现在QGIS中。我们可以按照与SAVI层相同的步骤更改此层的样式。 这次我们没有选择“绿色”查找表，而是选择“蓝色”。我们得到图4.18。

_images/image129_x2m.png

图4.18.根据2016年8月23日的Landsat 8图像计算出塞尔-庞松地区的AWEI。有 关该图的彩色版本，请参阅 www.iste.co.uk/baghdadi/qgis1.zip

在此图像上，深蓝色像素是AWEI值为正的像素。因此，这些像素显示水面。 白色像素是AWEI值为负的像素，因此是陆地像素。 然而，由于我们研究区域的地形，阴影和雪带往往会出现正值。在这里，我们面临着该指数的局限性之一。 当我们构建用于对土地利用进行分类的决策树时，我们将在以下部分中处理这个限制。

土地利用分类

决策树的定义

我们现在将使用原始陆地卫星频段、衍生指数和决策树对该地区的土地利用进行分类。 如前所述，该技术允许用户准确定义每个频段和每个土地利用类别的反射率值。 这种方法可能很耗时，并且可能需要用户的专业知识来实施。 但是，最终我们可以准确地定义不同的土地使用类别。

在这里，我们想根据8个类别对塞尔-庞松地区的土地利用进行分类： 干雪、融雪、水面、裸地/人造土壤、稀疏植被、草地/文化、林地和茂密的森林。

让我们从定义雪开始。探测到这一现象的第一张有用的图像是来自AWEI的那张。 在这一层上，我们可以看到雪地的像素大于0.12。但是，水像素也呈现大于0.12的AWEI值。 因此，我们需要找到一个标准来区分自由水和雪。 通过观察其他层，我们可以看到SAVI可能有助于区分这两个类别。 事实上，大多数时候，雪地的SAVI大于-0.25，而自由水域像素的SAVI低于此值。 最后，我们得到了我们的两个标准来表示这两类之间的差异。 如果某个特定像素的AWEI大于0.12，并且其SAVI大于-0.25，则该像素将被雪覆盖。 如果另一个像素的AWEI也大于0.12，但如果其SAVI小于-0.25，则该像素对应于开阔水域。 请注意，我们选择AWEI值0.12而不是0来检测水。此阈值允许用户避开阴影区域中的某些像素。

_images/image130_xp1.png

图4.19.用于对塞尔地区土地利用进行分类的决策树。有关该图形的彩色版本，请访问 www.iste.co.uk/baghdadi/qgis1.zip

通过对所有类使用类似的方法，我们获得以下决策树（图4.19）。 例如，AWEI值低于0.12、 第一带（深蓝色）反射率低于0.17且SAVI值在0.3和0.7之间的像素是稀疏植被的像素。

决策树的实施

定义此决策树后，我们将使用 SAGA GIS栅格计算器实现它。 我们可以通过在“处理工具箱”的搜索栏中输入该模块的名称来找到该模块。 在此模块中，我们首先需要设置将用作输入的层。在“主输入层”行，我们选择要使用的第一层。 例如，我们选择“awei_L8_2016-08-23”，这是我们之前计算的AWEI。 在我们将使用的公式中，该层将被命名为“a”。 然后，我们通过点击行末尾的图标“…”来访问子菜单“附加层”。 在这里，我们选择将在决策树中使用的其他层， 即SAVI（“savi_L8_2016-08-23_TOA_SP”）和第一个Landsat频段（“B1_L8_201608-23”）。 根据选择的顺序，这些层将分别命名为“b”和“c”。

主要困难是将我们的决策树翻译成 SAGA GIS语法。

我们使用逻辑 ifelse 指令：

ifelse(test, value if True, value if False)

例如，如果我们处理AWEI的决策树的第一层，我们会写这样的内容：

ifelse(Is the AWEI greater than 0.12, if yes I test if the SAVI is less
than -0.25, if no I test if the band 1 is greater than 0.17)

通过使用几个重叠的 ifelse 指令，我们获得以下决策树。我们将其写入计算器的“公式”部分。

ifelse(a>0.12,ifelse(b<(-0.25),1,2),ifelse(c>0.17,3,ifelse(b<0.3,4,ifelse(b<0.7,5,
ifelse(b<1,6, ifelse(b<1.17,7,8))))))

请注意，当我们处理负值时，我们需要将它们写在括号之间。 因此，我们测试的是“b<（-0.25）”，而不是“b<-0.25”。当所有设置完成后，我们得到图4.20。

在“Calculated”行中，我们将生成的分类命名为“landuse_ 8_classes_2016-08-23.tif”， 并将其存储在我们的工作目录中。过程结束后，分类直接出现在QGIS中。 我们可以通过其属性和“常规”面板更改其名称。 我们现在拥有Serre-Ponçon地区的土地利用分为8类（表4.2）。

_images/image131_xhc.png

图4.20.使用 SAGA GIS栅格计算器实现决策树

Pixel value

Land use

Hexa code

Color

1

Water

#001dff

2

Dry snow

#257cff

3

Melted snow

#13d0ff

4

Bare soil/artificialized soil

#8c8c8c

5

Sparse vegetation

#ffbe3b

6

Grassland/culture

#00db36

7

Woodland

#02a200

8

Dense forest

#055c0f

表4.2. Serre-Ponçon地区的土地利用分为8类。有关该表格的彩色版本， 请参阅seewww.iste.co.uk/baghdadi/qgis1.zip

我们可以通过其“属性”和“波段渲染”面板更改该层的风格。 可以通过将“渲染类型”设置为“单带伪色”来手动定义样式。 在我们的示例中，我们将通过加载文件“L8_landuse_8_classes.qml”来使用预定义的样式。 加载此文件后，我们单击“Bandrendering”面板底部的“Style”图标，并选择“Load Style…”菜单。 我们获得了图4.21所示的分类。

_images/image132_x24.png

图4.21. Serre-Ponçon地区的土地利用分为8类，摘自2016年8月23日Landsat 8图像。 有关该图的彩色版本，请参阅 www.iste.co.uk/baghdadi/qgis1.zip

使用多数过滤器简化土地利用

我们使用的方法（决策树）对土地利用进行分类是一种“基于像素”的方法。 图像的每个像素根据其在不同层中的值被分类到特定的类别中。 这种方法的主要缺点之一是产生具有“盐和胡椒”效应的分类， 其中许多孤立的像素被不良分类并给结果添加了一些“噪音”。 可能很难准确地定义什么是“噪音”和什么是“现实”，尤其是对于高分辨率图像。 然而，非均匀区域中间的孤立像素可能是噪音。 减少这种噪音的一种常用方法是对结果分类应用“多数过滤器”。

多数过滤器是一个移动窗口，其大小和形状由用户定义。该算法在窗口内寻找多数值，并替换中心像素的值， 如果该值最初不同，则替换该值。有一种使用 SAGA GIS应用这种过滤器的方便方法。 我们使用名为“多数过滤器”的模块。在新窗口中，我们首先设置要过滤的层。 这里，我们选择8类土地利用分类“landuse_8_classes_2016-08-23.tif”。 然后，在“搜索模式”行中，我们设置想要使用的窗口的形状，方形或圆形。 在我们的情况下，我们选择“方形”窗口。“半径”对应于窗口的大小。 例如，我们将大小设置为5像素x 5像素。最后，在“Filtered grid”行中注明生成层的名称。 我们将其命名为“landuse_L8_2016-08-23_filter5.tif”，并将其存储在我们的工作目录中（图4.22）。

_images/image133_xb1.png

图4.22.在8类土地上应用多数过滤器 使用 SAGA GIS进行使用分类

该过程完成后，我们可以通过执行初始分类风格的“复制和粘贴”来更改过滤分类的风格。 在过滤后的分类中，噪音更少，不同土地利用斑块之间的边界更全面（图4.23）。 我们必须小心，不要过于简化初始分类。

在这个阶段，有必要计算混淆矩阵以验证我们的分类。为此，我们需要一些地面控制点。 在专门介绍清算调查的章节[OSE 16]中介绍了混淆矩阵的一个示例。

_images/image134_xhy.png

图4.23.对初始分类应用5 x 5平方多数过滤器后的土地利用分类。 有关该图的彩色版本，请参阅www.iste.co.uk/baghdadi/qgis1.zip

减少土地利用类别数量

在某些情况下，将土地利用简化为更通用的类别可能会很有用。 例如，我们将8类土地利用简化为4类土地利用。我们正在寻找以下类别：水、雪、植被和裸土。 表4.3总结了这种简化。

Initial class

Initial value

Final class

Final value

Water

1

Water

1

Dry snow

2

Snow

2

Melted snow

3

Snow

2

Bare soil/art. soil

4

Bare soil

3

Sparse vegetation

5

Vegetation

4

G rassland/culture

6

Vegetation

4

Woodland

7

Vegetation

4

Dense forest

8

Vegetation

4

表4.3.将8类土地利用重新归类为4类土地利用

我们将使用名为“重新分类值（简单）”的 SAGA GIS模块继续进行此简化。 首先，我们需要设置我们想要重新分类的栅格。

我们选择了8个类别中的过滤分类“landuse_L8_2016-0823_filter5.tif”。 在“Lookup table”行中，我们指定我们想要更改的值以及我们想要设置的新值。 为此，我们单击“Lookup table”行中的图标“…”（图4.24）。

_images/image135_xv0.png

图4.24.土地利用重新分类

我们使用表4.3中的信息填充查找表窗口（图4.25）。

_images/image136_x0m.png

图4.25.适合将8类土地利用重新归类为4类土地利用

通过单击“确定”，我们返回到上一个窗口。在“Operator”行中，我们需要设置要使用的操作员。 由于我们定义查找表的方式，我们还想重新分类每个类的最小值和最大值。 因此，我们将操作符设置为“Low value <= grid value <= high value“，这是第二个选择。 最后，我们可以取消选中“reclassing no data values”和“replace other values”框。 我们将新分类命名为“landuse_L8_2016-0823_4_classes”（图4.24），并将其存储在我们的工作目录中。 我们获得了图4.26中所示的层。

_images/image137_x7c.png

图4.26.土地利用重新分类为四个主要类别。 有关该图的彩色版本，请参阅 www.iste.co.uk/baghdadi/qgis1.zip

上一页

16.1. QGIS中的 SAGA GIS

下一页

16.3. QGIS中利用 SAGA GIS提取水文网络