基于 landsat 图像用NDWFI来确认地表水范围

引用: Cai Yaotong, Qian Shi, Xiaoping Liu. Spatiotemporal mapping of surface water using Landsat images and spectral mixture analysis on Google Earth Engine. (2024) Journal of Remote Sensing
NDWFI
landsat图像

整体概览

1. 应用缩放因子：applyScaleFactors 函数用于调整 Landsat 影像的光学波段和热波段的比例。

2. 获取 Landsat 影像数据：从 GEE 的 Landsat 8 影像集合中筛选特定时间和地点的影像，并应用缩放因子。

3. 去除山体阴影：使用 SRTM（Shuttle Radar Topography Mission）数据和太阳方位角、高度角来生成山体阴影遮罩。

4. 光谱混合分析（SMA）：定义了五种地表类型（绿色植被、非光合植被、土壤、暗物质、亮物质）的光谱端元，并使用这些端元对 Landsat 影像进行光谱解混。

5. 计算 NDWFI：使用光谱解混的结果来计算 NDWFI，这是一种用于区分水体和非水体的指数。代码中还提到了使用地形校正和不透水数据来提高映射精度，以及使用 GAUD（Global Artificial Impervious Area Dataset）作为不透水源。

6. 可视化：设置了用于 Landsat 影像和 NDWFI 的可视化参数，并将它们添加到地图上进行显示。

第一部分：缩放因子

// Function to apply scaling factors.
function applyScaleFactors(image) {
  var opticalBands = image.select('SR_B.').multiply(0.0000275).add(-0.2);
  var thermalBands = image.select('ST_B.*').multiply(0.00341802).add(149.0);
  return image.addBands(opticalBands, null, true)
              .addBands(thermalBands, null, true);
}

1. 用途：将 Landsat 影像的原始数字量化值（DN值）转换为反射率和亮温值
2. 热波段与光学波段：
   • 热波段（ST_B.*）：指的是热红外波段，用于分析地表温度等。
   • 光学波段（SR_B.）：包括可见光和近红外波段，用于分析地表特征如植被、水体等。
3. ST_B.* 和 SR_B. 的含义：

   • ST_B.*：选择所有热红外波段，* 是通配符，表示匹配所有以 ST_B 开头的波段。

   • SR_B.：选择所有光学波段，. 在这里足以表示选择以 SR_B 开头的所有波段。

4. . 和 * 的区别：
   • * 是通配符，用于匹配任何字符。
   • . 在这个上下文中用作选择特定前缀的所有项，功能类似于通配符。

第二部分：获取landsat影像数据

// Get Landsat imagery 
var dataset = ee.ImageCollection('LANDSAT/LC08/C02/T1_L2')
    .filterDate('2021-05-01', '2021-06-01')
    .filterBounds(ee.Geometry.Point(-122.0808, 37.3947))
    .sort('CLOUD_COVER')
    .map(applyScaleFactors)
    .first();

这一串代码非常常用：先选择影像集合，再筛选日期【filterdate】，以及选中ROI【filterbounds】

1. .sort('CLOUD_COVER')：对筛选出的影像按照云覆盖度进行排序。云覆盖度较低的影像排在前，优先选择云量较少的影像，以便获得更清晰、更准确的地表信息。
   • sort('CLOUD_COVER', false) // 修改这里，将云覆盖度从高到低排序
2. .map(applyScaleFactors)：这一步是对每个影像应用 applyScaleFactors 函数。这个函数的作用是对影像中的光学波段和热波段数据进行缩放处理，将原始的数字量化值转换为反射率和温度等更有用的物理量。
   • map 函数: 是对影像集合中的每个影像执行指定的操作或函数。
3. .first()：从经过排序和处理的影像集合中选择第一个影像。这里选择的是云覆盖度最低的那个影像。

第三部分：去除山体阴影

// remove hillshadow using SRTM, azimuth and zenith
var dem = ee.Image("USGS/SRTMGL1_003").select('elevation')

//hillshadow函数
function hillshadow (img){
var shadowMap = ee.Terrain.hillShadow({
  image: dem.clip(img.geometry()),
  azimuth: img.get('SUN_AZIMUTH'),
  zenith: img.get('SUN_ELEVATION'),
  neighborhoodSize: 200,
  hysteresis: true
})
return shadowMap
}

var shadow_mask = hillshadow(dataset);

SRTM 高程数据

美国地质调查局（USGS）的 SRTM（Shuttle Radar Topography Mission）全球1弧秒数字高程模型（DEM）。这个模型提供了地表的高程信息。
L1_003：

• L1：这通常表示数据的处理级别。可能表示 Level 1，即一定级别的处理已经被应用于原始数据。
• 003：表示数据集的版本或更新号。

hillshadow 函数：利用地形数据和太阳位置生成一个山体阴影图像

首先在SRTM中选择了高程数据，命名为dem
传入img作为参数，设置一个变量为shadowmap
ee.Terrain.hillShadow 是 Google Earth Engine (GEE) 中的一个方法。这个方法用于根据给定的参数（如数字高程模型、太阳方位和高度角）生成一个表示山体阴影的影像。这里的 hillShadow 是 GEE 的一个内置函数，用于特定的地形分析任务。

• dem.clip(img.geometry()) 中，img.geometry() 被用于获取 img 影像的几何形状。然后，这个几何形状被用作 clip 方法的参数，以裁剪 dem 影像。
• azimuth: 太阳的方位角，从输入图像 img 中获取。
• zenith: 太阳的天顶角，也是从输入图像 img 中获取。
• neighborhoodSize: 用于计算阴影的邻域大小。精度
• hysteresis: 一个布尔值，用于指定是否应用滞后效应以改善阴影的连续性。稳定性

hillshadow函数返回shadowmap

ee.Terrain.hillShadow 函数的参数是一个 JavaScript 对象，它在语法和用途上类似于字典。在 JavaScript 中，对象是一种复合数据类型，它存储的是键值对（key-value pairs）。这些键值对可以包含各种数据类型，如字符串、数字、数组、甚至是其他对象或函数。

shadow_mask

传入dataset作为hillshadow函数的参数，最后返回一个遮罩

第四部分：光谱混合分析

光谱解混：用于从混合的像素光谱中估计地表覆盖类型的比例。这种方法基于假设一个像素的光谱可以表示为几种基本地表类型（称为端元）光谱的线性组合。

/ /********  parameters for spectral mixture analysis (SMA)**************
// 'gv', 'npv', 'soil', 'dark', 'bright' are the endmembers. 
// ******************************************************************/
var endmembers = ee.Dictionary({
  'gv': [0.019,	0.036,	0.0222,	0.303,	0.144, 0.054],
  'npv': [0.028, 0.068, 0.0044, 0.403, 0.224, 0.095],
  'soil': [0.034,	0.069,	0.074,	0.335,	0.274, 0.133],
  'dark': [0.064, 0.097, 0.087, 0.088, 0.046, 0.031],
  'bright': [0.066, 0.087, 0.088, 0.128, 0.123, 0.104],
});

var gv = endmembers.get('gv');
var npv = endmembers.get('npv');
var soil = endmembers.get('soil');
var dark = endmembers.get('dark');
var bright = endmembers.get('bright');
// /********  parameters for spectral mixture analysis (SMA)**************
// 'gv', 'npv', 'soil', 'dark', 'bright' are the endmembers. 
// ******************************************************************/
var endmembers = ee.Dictionary({
  'gv': [0.019,	0.036,	0.0222,	0.303,	0.144, 0.054],
  'npv': [0.028, 0.068, 0.0044, 0.403, 0.224, 0.095],
  'soil': [0.034,	0.069,	0.074,	0.335,	0.274, 0.133],
  'dark': [0.064, 0.097, 0.087, 0.088, 0.046, 0.031],
  'bright': [0.066, 0.087, 0.088, 0.128, 0.123, 0.104],
});

var gv = endmembers.get('gv');
var npv = endmembers.get('npv');
var soil = endmembers.get('soil');
var dark = endmembers.get('dark');
var bright = endmembers.get('bright');

// spectral mixture analysis (SMA)光谱解混分析
var unmixi = ee.Image(dataset)//从dataset的影像集合中创建unmixi
            .select(['SR_B1', 'SR_B2', 'SR_B3', 'SR_B4', 'SR_B5', 'SR_B7'])
            .unmix([gv, dark, npv, soil, bright],true,true)
            .rename(['GV','Dark', 'NPV','Soil','Bright']);

unmix方法

1. 线性混合模型：unmix 方法基于线性混合模型，该模型假设每个像素的光谱是其组成端元光谱的加权和。这意味着一个像素的反射率可以表示为不同端元反射率的线性组合。
2. 端元和权重：端元是预先定义的，代表特定地表类型的光谱特征（如绿色植被、土壤、水体等）。权重是解混过程中计算出来的，表示每个端元在像素中的相对比例。

接受几个参数：

1. 端元数组：包含了一组端元（endmembers）。端元是指在光谱解混分析中用作参考的光谱特征。在您的代码中，这个数组是 [gv, dark, npv, soil, bright]，分别代表不同的地表类型（如绿色植被、暗物质、非光合植被、土壤、亮物质）
2. 非负约束：这是一个布尔值（true 或 false），指定是否对解混结果应用非负约束。如果设置为 true，则所有解混得到的端元比例都不会小于零。
3. 和为一约束：这也是一个布尔值，指定是否对解混结果应用和为一的约束。如果设置为 true，则每个像素的所有端元比例之和将等于一。

rename

如果在 unmix 方法之后不使用 rename 方法，解混后的影像波段将保留默认的命名或没有明确的命名。

第五部分：计算NDWFI

NDWFI的原理是利用水体的暗端元比例较高
公式

/ ************************Calculate NDWFI**************************
// 1. Please use terrain corrections and impervious data to improve mapping accuracy.
// 2. The GAUD was selected as the impervious source in our study.
// 3. Please refer to terrain corrections in Poortinga et al. (2019).
// 4. Zero is recommended as the threshold for discriminating water and non-water.
// 5. You can also use the OTSU method to determine the optimal threshold when conducting a regional study.
// ******************************************************************/
var NDWFI = unmixi.expression(
'(dark-bright-soil-gv-npv)/(dark+bright+soil+gv+npv)',
{
'dark':unmixi.select('Dark'),
'soil':unmixi.select('Soil'),
'npv':unmixi.select('NPV'),
'gv':unmixi.select('GV'),
'bright':unmixi.select('Bright'),
}).rename('NDWFI')
.updateMask(shadow_mask);

var visualization = {
  bands: ['SR_B5', 'SR_B4', 'SR_B3'],
  min: 0.0,
  max: 0.2,
};
var palette = ['e8ebef','fd6f7c','fc8129','b3c120','4da910','1e7d83','001c83']
Map.setOptions('SATELLITE')
Map.setCenter(-122.0077, 37.4454, 12);
Map.addLayer(dataset, visualization, 'Landsat imagery');
Map.addLayer(NDWFI,{min:0,max:1,palette:palette},'NDWFI')

1. expression 方法：

• 用途：expression 方法用于对影像应用一个数学表达式。允许您对影像中的像素值进行复杂的数学运算。

• 参数：

  • 一个字符串：表示要应用的数学表达式。表达式是 '(dark-bright-soil-gv-npv)/(dark+bright+soil+gv+npv)'。
  • 一个字典（键值对），用于定义表达式中变量的来源。这些变量（如 dark, soil, npv, gv, bright）分别对应于 unmixi 影像中的不同波段。

• 作用：在这个例子中，expression 方法用于计算标准化差异水体分数指数（NDWFI）。这个指数是通过特定波段的线性组合计算得出的，用于区分水体和非水体区域。

2. updateMask 方法：

• 用途：用于应用一个掩膜（mask）到影像上。掩膜决定了影像的哪些部分应该被保留，哪些部分应该被忽略或隐藏。

• 参数：接受一个影像作为参数，这个影像的像素值决定了原始影像的哪些部分被保留。像素值为 1（或真）的地方被保留，像素值为 0（或假）的地方被忽略。

• 作用：在您的代码中，updateMask(shadow_mask) 使用 shadow_mask 作为掩膜来更新 NDWFI 影像。这意味着只有在 shadow_mask 中值为 1 的区域（即没有山体阴影的区域）的 NDWFI 值会被保留，而其他区域（山体阴影区域）则会被忽略。