第十二章 综合案例：咸海 40 年水体变迁

本文整理自我的原创教程《通过遥感历年水体变化监测看咸海40年变迁！》，内容基本保持原貌。 这是前十一章所有方法的综合实战：数据获取（第三、四、五章）→ 预处理（第八、九章）→ 指数提取（第七章）→ 精度意识（第十章），最终完成一项完整的遥感变化监测。

前言：咸海危机

咸海，曾经世界第四大湖泊。自上世纪60年代起，由于农业灌溉需求的急剧增加， 阿姆河和锡尔河流域开始大规模调水，导致咸海面积不断缩小，水位不断下降， 原本一整片的大湖也开始走向一分为二。

渐渐地，干涸的湖底变成了盐碱度极高的盐漠，不再适宜耕作； 草原上吹拂的狂风也变成了有害的盐尘暴，肆意威胁着生命。

90年代后，中亚国家意识到咸海危机，开始着手计划治理咸海。 但由于咸海地区每年高于十倍降水的蒸发量和超强度的工农业用水， 远超水循环负荷，让咸海消亡的趋势愈发不可收拾。

在这篇案例中，将用 ArcGIS Pro 和 Landsat 系列影像， 带你看看咸海近40年来是如何从大湖一步步走向消亡的！

1 确定研究区基本信息

1.1 投影坐标系

为了统计咸海水体面积，选择了 Albers 投影（等面积圆锥投影）作为地图的投影坐标系， 并根据咸海的地理位置，设置中央经线为 59.6°E，两条标准纬线分别为 44°N 和 46°N。 统计面积必须用等面积投影，这一点不能省。

1.2 研究区行列号

也就是确定研究区所需影像的位置。对于Landsat系列卫星影像， 官方采用 WRS（Worldwide Reference System，全球参考系统） 来标记影像所覆盖的 地理位置，通过 Path（路径）和 Row（行）两个数值，我们就能知道要下载哪些位置的影像。

首先在 ArcGIS Pro 中叠加研究区范围、WRS（daytime）与卫星底图图层， 确定研究区需要哪些位置的影像。要拼接出咸海研究区的卫星图像， 需要以下行列号的影像：

为了减轻工作量，这里计划只下载 Row 为 27、28、29 行的影像。

1.3 历时范围与月份

计划分析咸海近40年每隔5年的变化，即从1984年到2024年，每5年一景影像， 所以选取 1984、1989、1994、1999、2004、2009、2014、2019、2024 作为观测年份。

由于咸海地处中亚干旱地区（哈萨克斯坦与乌兹别克斯坦交界处）， 云量特征属于季节性类型，夏秋季节云量较少，所以选取7、8、9月份的少云影像 作为观测数据。如果这三个月的影像都有严重的云层干扰， 就尽量选取时间靠近这几个月份的少云影像，以保持影像色调大致相同。

2 逐一下载影像

根据年份建立任务文件夹，将每个年份下指定月份的影像下载到对应年份文件夹中， 每个年份文件夹都包含研究区所有行列号的影像。

输入行列号（Row/Path）确定影像位置，设置年份和月份 （为了检索出结果，可以不用限制云量），下载 Landsat Collection 2 Level-2 数据集中的影像——因为 Level-2 中包含经过辐射定标和大气校正处理的地表反射率产品， 到后面可以直接用于地表分析，不需要额外的大气校正步骤。 具体下载方法见第五章（教程原文链接在其参考资料中）。

下载结果中有几点值得注意：

1. 1984年有部分影像缺失，用1985年的影像代替；
2. 个别年份的影像在7、8、9月份时云量较多，选取了其它较接近的月份代替；
3. 2009年的影像依然用的是Landsat 5影像——因为在此期间Landsat 7因传感器故障， 影像普遍有明显的缺失条纹；
4. 2004年的影像来自Landsat 7（没检索到Landsat 5的结果），后续需要填充缺失条纹。

问题反思：由于不太会利用 Earth Explorer 的 API，只能采取笨方法一幅幅地下载。 这意味着要弄到9个年份、每个年份9个行列号的影像，即便过程顺利， 也要重复检索八十一次。实际还不止——下载请求太频繁，IP直接被限制了， 只能换台电脑下，光下载影像就花费了一晚上和一早上。 这也是为什么值得学习第六章的 GEE 方案。

3 拼接影像

解压数据：将每个年份文件夹中的数据都解压到当前位置。

添加数据：打开 ArcGIS Pro，进入地图，然后【添加数据】并以类型排序， 将年份文件夹中所有的栅格产品（MTL.txt）添加进去。 添加后都默认显示的是地表反射率产品，以 SurfaceReflection 名称前缀或 SR 名称后缀标识。

镶嵌数据：利用【镶嵌至新栅格】工具将所有行列号的影像拼接成一幅完整的影像。 【像素类型】选择与原影像一样的16位无符号，【波段数】保持与原影像一致， 【镶嵌运算符】选择平均数（Mean），即对重叠区域的像素值按平均值处理。

去除黑色背景：栅格镶嵌或裁剪后都会产生黑色背景， 如果觉得影响观察可打开影像的符号系统，在掩膜中显示背景值并设为透明。 同时在这里也可以调整 RGB 对应的波段分别为 Band3、Band2、Band1 并标准拉伸， 使影像色彩更接近自然。

4 填充影像缺失条纹

由于Landsat 7的机载扫描行校正器（Scan Lines Corrector，SLC）在2003年5月31日 发生了故障，所以下载的Landsat 7影像普遍有缺失条纹的情况。 而这些缺失的部分将会导致后续统计的水体面积值偏小。

针对Landsat 7影像，缺失的部分可采取邻域平均值插值的方式进行填充。 其它Landsat系列的影像一般没有这种情况，可跳过此步。

将背景值设为 Nodata：利用【栅格计算器】工具， 输入 SetNull("目标栅格" == 0, "目标栅格") 函数公式， 将值为0的背景（包括缺失条纹）设为Nodata。

填充 Nodata：在栅格函数中利用【统计分析】， 以邻域平均值的方式仅填充Nodata像素。如果用【栅格计算器】， 则通过形如以下的表达式实现：

Con(IsNull("目标栅格"), FocalStatistics("目标栅格", NbrRectangle(3,3,"CELL"),"MEAN"), "目标栅格")

导出栅格：填充后的影像图层还并不是能进行地理处理的栅格数据， 需要将填充的图层导出为栅格，像素类型与原影像保持一致， 这样后续才能利用工具裁剪栅格。注意，不要勾选使用渲染器， 否则导出的栅格就是一张RGB图像，将丢失原始影像的波段。

5 裁剪影像

利用【裁剪栅格】或【按掩膜提取】工具，裁剪出研究区的影像。 或者更方便一点，在导出栅格时就选择用裁剪几何导出，省去裁剪步骤。 其中【裁剪类型】选择【外部】，即导出的栅格不保留研究区以外的栅格区域。

6 提取水体

在近红外波段，水体几乎只吸收而不反射，与其它地物具有明显区别 （原理见第一、四章）。要从遥感影像中识别并提取水体， 可采用 MNDWI（改进的归一化水体指数）法，指数值越大， 说明越可能是水体，这也是最便捷常用的方法之一。

与 NDWI 相比，MNDWI 考虑了近红外波段中土壤、建筑信息的干扰， 可以更有效地区分出水体：

MNDWI = (Green - SWIR) / (Green + SWIR)

计算 MNDWI 的参数只有两个：绿色波段值和 SWIR（短波红外）波段值， 其中 Landsat 影像的 SWIR 波长范围通常在 1.57～1.65 微米之间 （对应 Landsat 5/7/8/9 的 Band 5 / Band 6，注意查各传感器的波段表，见第四章）。

6.1 计算 MNDWI

在 ArcGIS Pro 中计算MNDWI非常方便，因为影像处理栏中提供了许多常见指数的工具， 无需输入公式：选择 MNDWI 指数，然后选择影像波段中对应的绿色波段和短波红外波段， 就会得到 MNDWI 结果（拉伸类型为标准差）：

其中云量会干扰 MNDWI 结果，后期需要去除。

6.2 导出 MNDWI 栅格

生成的 MNDWI 结果是32位浮点型，数值范围非常大，远超实际 MNDWI 值范围， 使得水体与非水体在直方图上显得几乎没什么区别，不方便之后的手动阈值划分； 而且32位的栅格值范围会超出 ArcGIS 的分类运算限制，导致程序崩溃。

所以还要再将 MNDWI 结果导出栅格，相当于重新计算栅格值范围， 让栅格值范围变成实际的 MNDWI 值范围。导出的像素类型保持与 MNDWI 结果一致， 都为浮点型。

6.3 阈值划分

接下来就是将一定阈值范围内的 MNDWI 值归于一个值，作为水体的标识； 而处于阈值范围外的 MNDWI 值则都归于另一个值，作为非水体的标识。

在 ArcGIS Pro 中，针对两类的阈值划分，可通过符号系统中的【分类】 或地理处理工具中的【重分类】手动划分；也可利用【二进制阈值】函数自动划分。

手动阈值划分：利用【重分类】工具，点击分类， 综合水体与非水体样本的 MNDWI 值和直方图中波峰两端点值，确定水体的阈值范围。 例如将水体的阈值范围确定为 0～0.338，并归为新建的 2（类别）， 这样之后栅格转面后，属性表中水体面的标识就是 2。

自动阈值划分：利用【二进制阈值】栅格函数，将 MNDWI 结果划分为水体与非水体， MNDWI 值大的区域将归为 1（水体）；值小的区域将归为 0（非水体）。 虽然【二进制阈值】划分简单方便，但没有参数控制， 生成的结果有时候会与实际偏差非常大、不够精确，就只能手动重分类了。

6.4 栅格转面与导出

利用【栅格转面】工具将阈值划分的结果转为面要素；打开属性表， 选择 gridcode 为 1（水体）的部分，按所选部分导出要素—— 这样就得到咸海的水体要素了。

6.5 排除云干扰

由于影像中的云会干扰 MNDWI 值，这也就导致栅格转面后， 裸地上空的云可能会被归于水体面要素中，水面上空的云可能会阻碍水体提取； 如果不去除的话，将会导致计算的水域面积值出现偏差。

所以对于由裸地上空云转换而来的水面要素，要将其删除； 而对于因水面上空云阻碍而缺失的水面要素则补齐： 对照着彩色影像，在编辑中选择由裸地上空云转换而来的水面要素，然后删除。

7 计算水体面积

融合：利用【融合】工具将水面要素中的所有部件融合成一个整体，方便计算总面积。

投影：目前要素的坐标系只有地理坐标系，单位是经纬度，不是米，还无法计算面积。 需要利用【投影】工具赋予要素第1.1节设置好的 Albers 投影坐标系。 这样就得到正确的 Shape_Area 了。例如 2004 年的咸海水体面积 为 25061356431.77 m² ≈ 25061.35 km²。

计算几何：如果要素属性表中没有自动计算好的面积字段， 则要在投影后打开要素属性表，手动计算几何面积。

8 模型批处理

基于以上流程方法，利用【模型构建器】创建批处理工具， 完成对其余年份影像的拼接、填充、裁剪与面积计算。 其中MNDWI 提取水体的步骤没办法用模型处理，估计要用上 Python—— 这也是从软件操作走向编程处理的一个自然动机。

9 统计汇总

最后导出布局，汇总面积至 Excel：

（以上所有数据仅供演示，不作实际参考依据。）

从数据可以读出两个阶段：1984-2004 年是匀速萎缩，每5年减少一到两成； 2004-2014 年是崩塌式退缩，2009 年单期暴跌 54.7%—— 南湖（咸海南部）在这一时期基本干涸。2019 年之后面积趋于稳定（+0.7%）， 但这不是治理的胜利，而是剩下的北湖已接近水量平衡的底部。

结语与补充

1. 有时候因为影像质量问题或地物信息干扰，MNDWI 并不能明显区分水体与非水体， 可以试试将自然彩色影像导出为灰度图，再重分类，最后转矢量。 （照这么说，其实有造假嫌疑——拿谷歌历史影像和 PS 就能复刻整个过程……）
2. 虽然有批处理，但依然有许多重复操作，而且下载数据过程重复繁琐、耗时过长， 用软件批处理还是有点不合适，或许遥感云计算（例如GEE，见第六章） 或编程处理才是遥感历年变化监测的不二之选。
3. 由于我并非地信遥感专业，只是感兴趣想试试能不能只用 GIS 就能完成遥感地学分析，所以水平有限，文中难免有不合理之处， 还请各位读者火眼金睛，多多建议！

参考资料

1. 王志信, 林友明, 黄鹏, 等. 基于时间序列的中亚地区云量特征分类及云量变化趋势[J]. 遥感技术与应用, 2014, 29(05): 839-845.
2. Landsat ETM+ SLC-off 遥感影像条带修复 - 开源地理空间基金会中文分会
3. Moradi, M., Sahebi, M.R., & Shokri, M. (2017). Modified Optimization Water Index (MOWI) for LANDSAT-8 OLI/TIRS. ISPRS Archives, 185-190.
4. ArcGIS 填补栅格空缺值 Nodata - CSDN博客
5. Landsat 数据集合集（Landsat 5/7/8/9）- CSDN博客