19.3. QGIS中利用SAGA GIS提取水文网络#
为了研究河流及其流域功能的水文、生物地球化学或其他方面，通常需要提取相关的水文网络。 在这个网络上，可以计算形态指数（长度、排水密度、弯曲指数等）。 这些计算网络是以矢量或栅格格式对水文系统的简化表示。 在矢量格式中，河流表示为定向多段线，流域表示为多边形。
水文网络提取主要是通过 DEM 来完成的。 有些几乎在全球范围内提供（GTOPO 30、SRTM）， 另一些则在国家层面提供（法国的BD Alti IGN），甚至在更多地方范围内提供。 这些数字元数据库的特征在于给定的空间分辨率（即像素大小）和测高精度， 这决定了将产生的水文系统表示的质量和精度。
通过一个专注于塞尔-庞松（法国阿尔卑斯山）同一地区的示例， 我们将了解如何使用QGIS及其 SAGA 模块来生成此类数据集。 SAGA 根据河流被认为是发源的唯一区域阈值的定义，实施了一些传统方法 [OCA84] 。 SAGA 还提出了基于“流焚烧”的方法 [HUT89] 。 在本例中，我们将使用特定于 SAGA 的方法，基于Strahler的排序 [STR57] 。
19.3.1. DEM 预处理#
我们可以从前面提到的EarthExplorer网站下载部分SRTM数据集（近全球数字模型-从纬度56° S到60° N -由NASA制作） [FAR07] 。
一旦识别，就可以以"GeoTivf"格式下载对我们感兴趣的区域进行扩展的两个SRTM切片。
我们选择了最新版本的SRTM，分辨率为1弧秒，相当于赤道上约90平方米的分辨率，
在WGS 84中进行了地理参考 [SYS 08]。
在 EarthExplorer 界面中，
我们在“数据集-数字海拔-SRTM”面板中选择产品“SRTM 1 Arc-SecondGlobal”。
以下磁贴对应于我们要研究的区域：
“n44_e005_1arc_v3.tif和n44_e006_1arc_v3.tif”。
为了优化我们的计算，
我们将合并这些SRTM瓦片并将它们投影到UTM 31 N区（EPSG 32631）坐标参考中，
如专门用于土地利用分析的部分所示。
由于 SAGA 地理信息系统主要处理在公制系统中投影的数据，因此该投影是强制性的。
用户也更容易用平方米而不是平方度来思考。
我们将首先马赛克我们的SRTM磁块，以便仅操作单个文件而不是两个文件。
我们在“栅格工具”菜单中使用 SAGA 工具“马赛克栅格层”。
在弹出窗口中，我们指定要在“输入网格”框中组装的两个瓦片。
我们将其他选项保留为默认值，
并将我们想要保存马赛克的名称指定为“mos_tmp.tif”，
弹出窗口，用于马赛克我们感兴趣区域的两个SRTM磁贴如 图 19.27 。
我们的结果现在显示在QGIS主窗口中。
图 19.27 两个SRTM磁贴#
我们现在需要将这个马赛克投影到“UTM 31 N”区，以进行米制和相关投影。 为此，我们单击该层，选择“另存为.”菜单项，然后在“另存为.”输入框中指定我们要存储新层的目录。 我们将其命名为“mos_tmp2.tif”， 最后通过单击其旁边的图标指定要将该层重新投影到“CRS”下拉框列表中的“UTM Zone 31 N”中。 我们可以选择此投影，其具有EPSG代码（32631），投影到UTM区域31 N坐标参考合并的SRTM切片层的投影， 如 图 19.28 所示。
图 19.28 合并的SRTM切片层的投影#
最后，为了优化计算时间，我们将根据形状文件“study_area_32631.shp”中定义的研究区域的边界来剪辑马赛克。 首先，我们将项目的投影设置为UTM 31 N。为此， 我们单击位于QGIS主界面右下角的包含单词“EPSG”的图标。 在弹出窗口中的“CRS”选项卡中， 我们通过EPSG代码将项目的坐标参考指定为“UTM 31 N”，如 图 19.29 所示。
图 19.29 设置项目的坐标参考系UTM 31 N#
我们将“study_area_32631.shp”层加载到QGIS中，
然后使用 Clipper. 通过“Raster & Extraction”菜单中的功能，
我们根据我们的研究区域来提取数字地形图。
在弹出窗口中，我们选择“mos_tmp2.tif”马赛克数据集作为输入文件，
并命名输出文件“srtm_durance_utm31n.tif”。
我们将零值设置为“-99999”（ SAGA 中的默认值）。
我们将“剪裁模式”设置为“掩蔽层”并指定要使用的层“srtm_sp_utm31n.shp”。
最后，我们检查“将目标数据集的范围裁剪到切割线的范围”选项，根据我们的兴趣领域剪辑SRTM数据集，
如 图 19.30 所示。
图 19.30 剪辑SRTM数据集#
我们现在获得剪辑到研究区域的SRTM数据集。 我们可以删除所有中间网格，并在剩余的过程中仅保留最后一个。
19.3.2. 填充水槽#
出于水文目的使用数字模型涉及填充可能包括的水槽的第一次计算，这将使理论上的河流流量不可能。 这种情况通常发生在不是河流出口的栅格单元被海拔比其本身更高的其他单元完全包围时，从而形成一种碗。 从地形上看，实地情况可能是这种情况（例如，崩解带中的裂缝和天坑）， 但也可能是由于数字地形图的分辨率或其采集条件而造成的文物。 因此，建立象征非水道网络的定向弧线网络需要填充这些单元，以避免网络中的不连续性。 该操作称为“填充水槽”，包括将从碗边缘的值内插到相关单元。 此外，在我们稍后将使用的模块中，还考虑了一般斜坡来计算插值。
有几种方法可以实现这一点。 SAGA 实现了 [PLA01] 的方法，通过模拟 DEM 上径流通道的剩余部分来填充它们。 对于练习，我们将使用 [WAN06] 的方法，其优点是在填充的单元格中保持坡度。 此功能位于 Terrain Analysis → Hydrology → Fill sinks (Wang & Liu) 菜单中。 在弹出窗口中，我们首先指定要填充的 DEM “srtm_sp_utm31n.tif”。 我们保留“Minimumslope”的默认值。 我们指定存储已填充的栅格的路径及其名称“srtm_sp_filled_utm31n.tif”。 由于我们不会使用“FlowDirections”和“Watershed bases”栅格， 因此我们取消选中与这些栅格相关的“执行后打开文件”框，如 图 19.31 所示。
图 19.31 填充水槽#
在过程结束时，我们获得一个与第一个非常相似的新 DEM ：只有被识别为汇的像元才具有新的海拔值。
19.3.3. 海道网络提取#
除了用于提取具有最小上游流域面积的水文网络的传统方法外，
SAGA 还提供了一种基于"Strahler"对数字地形图每个像素排序的原创方法。
其逻辑与河段排序完全相同。没有上游像素的像素是1阶像素。
当一个像素位于两个1阶像素的下游时，它就会变成2阶像素。
如果这个2阶像素接收一个1阶像素，则该像素保持2阶。
只有当它接收到"forder 2"的另一个像素时，它才会传递到"order 3"，等等。
使用这种方法，我们必须在斯特拉阶上设置一个阈值，超过这个阈值，我们估计河流可以常年流动。
如果油非常不透水，那么这个阈值就会很低。
相反，当土壤具有渗透性时，这个阈值就会很高。
我们将这种方法应用于“地形分析-渠道”菜单中充满“渠道网络和流域”模块的 DEM 。
窗口如 图 19.32 所示。
图 19.32 使用Strahler排序方法提取河流网络#
我们首先指定我们希望处理的 DEM ；在我们的情况下，
我们选择填充的 DEM “srtm_sp_fill_utm31n.tif”。
默认阈值等于5。该阈值意味着任何Strahler度数等于或大于5的像素都将被视为常年河流。
我们可以在第一步保持这个阈值。
我们仅选择将“Strahler Order”栅格（我们将其命名为“srtm_sp_strahler_utm31n.tif”）和“Channels”水道网络保存为多段线的形状文件（我们将其命名为“channels 5_sp_utm31n.shp”）。
在Serre-Ponçon地区（绿色）中使用Strahler阈值5提取的水文网络秩序，出现以下结果，
如 图 19.33 所示。
图 19.33 Strahler阈值5提取的水文网络秩序#
背景图像是 DEM （SRTM.对于彩色版本图见 www.iste.co.uk/baghdadi/qgis1.Zip ）
为了充分理解 SAGA 使用的逻辑和Strahler阈值的意义，我们可以放大河流的源头，
它可以被视为常年春天，如 图 19.34 所示。
图 19.34 放大两个河流来源（红色）#
后向屏幕是逐像素的Strahler排序栅格“srtm_sp_strahler_utm31n.tif”（像素大小：30 m x 30 m）。 有关该图的彩色版本，请参阅www.iste.co.uk/baghdadi/qgis1.zip
在图4.34中，我们看到两个河流源头呈红色。 基础栅格是逐像素的Strahler排序栅格。 因此，每个像素具有与其Strahler流顺序对应的值。 最暗的像素的阶数为-3，而较亮的像素的阶数为1。 根据定义，Strahler流顺序不能为负，它应该从1开始。 事实上，在这种情况下，我们必须关注等级而不是价值。 如果我们将像素重新分类以将-3设置为1，我们就会得到"栅格的Strahler流订单和“真实”Strahler订单之间的等效性" 如 表 19.4 所示。
Raster Strahler order |
-3 |
-2 |
-1 |
0 |
1 |
“Real” Strahler order |
1 |
2 |
3 |
4 |
5 |
结果表明，河流的“真实”斯特拉阶数为5。 因此，我们与之前定义的阈值完全一致。 由于这种转变，第一批河流的水文相关斯特拉级为1。
为了了解对"Strahler"阈值的敏感性，我们可以重新提取阈值为8的水文网络。 我们将结果称为“channels8_sp_utm31n.shp”。 我们阈值为5（绿色）提取的网络和阈值为8（黄色）提取的网络之间的比较， 获得如 图 19.35 所示的网络。
图 19.35 网络之间的比较#
有关该图的彩色版本，请参阅 seewww.iste.co.uk/baghdadi/qgis1.zip
正如预期的那样，用阈值8提取的网络比用阈值5提取的网络密度低。 该阈值的定义取决于所研究的地区。 它可以通过与用户有信心的水文网络（地形图、卫星图像、法国IGN等国家数据库）进行比较来确定 [IGN14], etc.) 。
有很多方法可以从数字元中提取河流网络。可以比较它们，就像 施奈德 等人。 [SCH17] 他们提出了不仅基于海拔数据而且还基于气候和岩石数据的新方法。
19.3.4. 参考文献#
FARR T.G., ROSEN P.A., CARO E. et al., “The Shuttle RadarTopography Mission”, Reviews of Geophysics, vol. 45, 2007.
HUTCHINSON M., “A new procedure for gridding elevation andstream line data with automatic removal of spurious pits”, Journal ofHydrology, vol. 106, nos. 3–4, pp. 211–232, 1989.
IGN, Base de Données sur la Cartographie Thématique des Agencesde l’eau et du ministère chargé de l’environnement, Institut national del’information géographique et forestière, 2014.
O’CALLAHAM J.F., MARK D.M., “The extraction of drainagenetworks from digital elevation data”, Computer Vision, Graphics, andImage Processing, vol. 28, pp. 323–344, 1984.
PLANCHON O., DARBOUX F., “A fast, simple and versatilealgorithm to fill the depressions of digital elevation models”,*CATENA*, vol. 46, nos. 2–3, pp. 159–176, 2002.
SCHNEIDER A., JOST A., COULON C. et al., “Global-scale rivernetwork extraction based on high-resolution topography and constrainedby lithology, climate, slope, and observed drainage density”,*Geophysics Research Letters*, vol. 44, pp. 2773–2781, 2017.
STRAHLER A.N., “Quantitative analysis of watershedgeomorphology”, Transactions of the American Geophysical Union, vol.38, no. 6, pp. 913–920, 1957.
WANG L., LIU H., “An efficient method for identifying andfilling surface depressions in digital elevation models for hydrologicanalysis and modeling”, International Journal of GeographicalInformation Science, vol. 20, no. 2, pp. 193–213, 2006.