使用Google Earth Engine计算1994-2000年Landsat 5 TM数据的Forel-Ule指数（FUI）

一、什么是Forel-Ule指数（FUI）？

Forel-Ule指数（FUI）是一种用于描述水体颜色的分类系统，最初由Forel和Ule于19世纪末提出。该指数将水体颜色分为21个等级，从蓝绿色（1级）到棕色（21级），可以有效反映水体的光学特性和水质状况。

在遥感应用中，FUI通常基于多光谱卫星影像的反射率数据计算得出，能够帮助研究人员大规模、长期监测水体颜色变化，进而分析水质变化趋势、富营养化状况等环境问题。

二、研究数据与区域

1. 数据源

本研究使用的是Landsat 5 TM传感器的Level 2级表面反射率数据（LANDSAT/LT05/C02/T1_L2），时间范围为1994-2000年。Landsat 5卫星携带的TM传感器具有7个波段，其中我们主要使用了以下4个光学波段进行FUI计算：

SR_B1：蓝光波段（0.45-0.52μm）
SR_B2：绿光波段（0.52-0.60μm）
SR_B3：红光波段（0.63-0.69μm）
SR_B4：近红外波段（0.76-0.90μm）

2. 研究区域

研究区域可以通过修改代码中的geometry变量进行设置，默认使用的是"projects/ee-hllutlu2024/assets/DT2"路径下的矢量数据。实际应用时，只需替换为自己的研究区域矢量数据即可。

三、完整代码解析

1. 研究区域加载

var geometry = ee.FeatureCollection("projects/ee-hllutlu2024/assets/DT2"); // 替换为实际研究区域
Map.centerObject(geometry, 9);

这段代码的作用是加载研究区域的矢量数据，并将地图中心定位到该区域，缩放级别为9级。

2. Landsat 5 TM数据预处理函数

function preprocessL5(image) {
  var qaPixel = image.select('QA_PIXEL');
  var qaRadsat = image.select('QA_RADSAT');
  
  // 云阴影掩码
  var cloudShadowMask = qaPixel.bitwiseAnd(parseInt('11111', 2)).eq(0);
  // 饱和度掩码
  var saturationMask = qaRadsat.eq(0);
  
  // 光学波段反射率计算（将DN值转换为实际反射率）
  var opticalBands = image.select(['SR_B1', 'SR_B2', 'SR_B3', 'SR_B4'])
                          .multiply(0.0000275)
                          .add(-0.2);
  
  // 返回处理后的影像
  return image.addBands(opticalBands, null, true)
              .updateMask(cloudShadowMask)
              .updateMask(saturationMask)
              .clip(geometry);
}

预处理函数主要完成以下工作：

利用QA_PIXEL波段创建云阴影掩码，去除云阴影影响
利用QA_RADSAT波段创建饱和度掩码，去除饱和像元
将光学波段的DN值转换为实际反射率（应用缩放因子0.0000275和偏移量-0.2）
对影像进行裁剪，只保留研究区域内的部分

3. FUI计算核心函数

function calculateFUI(meanImage) {
  // 计算归一化水体指数(NDWI)并提取水体
  var ndwi = meanImage.normalizedDifference(['SR_B2', 'SR_B4']);
  var waterMask = ndwi.gte(0.25).selfMask();
  var waterReflectance = meanImage.updateMask(waterMask);

  // 计算X、Y、Z三刺激值
  var X = waterReflectance.expression(
    '11.053 * b1 + 6.950 * b2 + 51.135 * b3 + 34.457 * b4',
    {b1: waterReflectance.select('SR_B1'), b2: waterReflectance.select('SR_B2'),
     b3: waterReflectance.select('SR_B3'), b4: waterReflectance.select('SR_B4')}
  ).rename('X');

  var Y = waterReflectance.expression(
    '1.32 * b1 + 21.053 * b2 + 66.023 * b3 + 18.034 * b4',
    {b1: waterReflectance.select('SR_B1'), b2: waterReflectance.select('SR_B2'),
     b3: waterReflectance.select('SR_B3'), b4: waterReflectance.select('SR_B4')}
  ).rename('Y');

  var Z = waterReflectance.expression(
    '58.038 * b1 + 34.931 * b2 + 2.606 * b3 + 0.016 * b4',
    {b1: waterReflectance.select('SR_B1'), b2: waterReflectance.select('SR_B2'),
     b3: waterReflectance.select('SR_B3'), b4: waterReflectance.select('SR_B4')}
  ).rename('Z');

  // 计算色度坐标
  var chrx = X.divide(X.add(Y).add(Z)).rename('chrx');
  var chry = Y.divide(X.add(Y).add(Z)).rename('chry');
  
  // 计算色调角度并转换为度数
  var hueRadians = chrx.subtract(0.3333).atan2(chry.subtract(0.3333));
  var hueDegrees = hueRadians.multiply(180).divide(Math.PI);
  var normalizedHue = hueDegrees.divide(100);

  // 色调校正
  var correction = normalizedHue.pow(5).multiply(25.81)
                   .add(normalizedHue.pow(4).multiply(-177.4))
                   .add(normalizedHue.pow(3).multiply(476.69))
                   .add(normalizedHue.pow(2).multiply(-653.3))
                   .add(normalizedHue.multiply(463.3))
                   .add(-94.41);
  var correctedHue = hueDegrees.add(correction);

  // FUI分配逻辑（覆盖所有范围并转换为整数）
  var fui = correctedHue
    .where(correctedHue.lte(19.0), 1)
    .where(correctedHue.gt(19.0).and(correctedHue.lte(22.741)), 1)
    .where(correctedHue.gt(22.741).and(correctedHue.lte(26.337)), 2)
    .where(correctedHue.gt(26.337).and(correctedHue.lte(30.439)), 3)
    .where(correctedHue.gt(30.439).and(correctedHue.lte(34.906)), 4)
    .where(correctedHue.gt(34.906).and(correctedHue.lte(39.769)), 5)
    .where(correctedHue.gt(39.769).and(correctedHue.lte(45.129)), 6)
    .where(correctedHue.gt(45.129).and(correctedHue.lte(50.665)), 7)
    .where(correctedHue.gt(50.665).and(correctedHue.lte(56.435)), 8)
    .where(correctedHue.gt(56.435).and(correctedHue.lte(62.186)), 9)
    .where(correctedHue.gt(62.186).and(correctedHue.lte(67.957)), 10)
    .where(correctedHue.gt(67.957).and(correctedHue.lte(74.572)), 11)
    .where(correctedHue.gt(74.572).and(correctedHue.lte(83.346)), 12)
    .where(correctedHue.gt(83.346).and(correctedHue.lte(94.037)), 13)
    .where(correctedHue.gt(94.037).and(correctedHue.lte(109.054)), 14)
    .where(correctedHue.gt(109.054).and(correctedHue.lte(132.999)), 15)
    .where(correctedHue.gt(132.999).and(correctedHue.lte(163.084)), 16)
    .where(correctedHue.gt(163.084).and(correctedHue.lte(190.779)), 17)
    .where(correctedHue.gt(190.779).and(correctedHue.lte(209.994)), 18)
    .where(correctedHue.gt(209.994).and(correctedHue.lte(220.977)), 19)
    .where(correctedHue.gt(220.977).and(correctedHue.lte(227.168)), 20)
    .where(correctedHue.gt(227.168), 21)
    .round()  // 强制转换为整数
    .rename('FUI');

  return fui;
}

FUI计算函数是整个代码的核心，主要包含以下步骤：

水体提取：利用NDWI（归一化水体指数）提取水体区域，NDWI≥0.25的区域被视为水体
三刺激值计算：基于四个光学波段的反射率计算X、Y、Z三刺激值
色度坐标计算：将三刺激值转换为色度坐标(chrx, chry)
色调角度计算：根据色度坐标计算色调角度并转换为度数
色调校正：对色调角度进行校正，以更准确地匹配Forel-Ule颜色标准
FUI等级分配：根据校正后的色调角度，将水体划分为21个FUI等级

4. 批量处理1994-2000年数据

var startYear = 1994;
var endYear = 2000;

for (var year = startYear; year <= endYear; year++) {
  // 丰水期处理
  var wetSeason = ee.ImageCollection("LANDSAT/LT05/C02/T1_L2")
                      .filterDate(ee.Date.fromYMD(year, 4, 1), ee.Date.fromYMD(year, 9, 30))
                      .filterBounds(geometry)
                      .map(preprocessL5);
  var meanWet = wetSeason.mean();
  var fuiWet = calculateFUI(meanWet);
  Export.image.toDrive({
    image: fuiWet,
    description: year + "_FUI_wet",
    folder: "FUI_Export",
    scale: 30,
    region: geometry,
    maxPixels: 1e13,
    fileFormat: "GeoTIFF"
  });

  // 枯水期处理
  var drySeason = ee.ImageCollection("LANDSAT/LT05/C02/T1_L2")
                      .filterDate(ee.Date.fromYMD(year, 10, 1), ee.Date.fromYMD(year + 1, 3, 31))
                      .filterBounds(geometry)
                      .map(preprocessL5);
  var meanDry = drySeason.mean();
  var fuiDry = calculateFUI(meanDry);
  Export.image.toDrive({
    image: fuiDry,
    description: year + "_FUI_dry",
    folder: "FUI_Export",
    scale: 30,
    region: geometry,
    maxPixels: 1e13,
    fileFormat: "GeoTIFF"
  });
}

这段代码实现了1994-2000年数据的批量处理，将每年分为丰水期和枯水期分别处理：

丰水期：每年4月1日至9月30日
枯水期：每年10月1日至次年3月31日

对于每个时期，代码会：

加载对应时期的Landsat 5 TM影像集
应用预处理函数对影像进行处理
计算该时期的平均反射率影像
基于平均反射率影像计算FUI
将计算得到的FUI结果导出到Google Drive

导出参数说明：

scale: 30，Landsat TM数据的空间分辨率为30米
fileFormat: "GeoTIFF"，导出为GeoTIFF格式，便于后续分析
maxPixels: 1e13，设置最大像元数，避免因区域过大而导出失败

四、使用方法

登录Google Earth Engine平台
创建新的脚本，将上述代码复制粘贴到编辑器中
将geometry变量替换为自己的研究区域
根据需要调整时间范围（修改startYear和endYear）
点击运行按钮，代码会自动创建导出任务
在"Tasks"面板中点击每个任务的"Run"按钮，开始导出数据

五、注意事项

研究区域不宜过大，否则可能导致计算时间过长或导出失败
确保Google账号有足够的存储空间，导出的GeoTIFF文件会占用一定空间
Landsat 5数据在不同年份的质量可能有差异，部分年份可能存在数据缺失
水体提取的NDWI阈值（0.25）可能需要根据研究区实际情况进行调整
FUI计算中的参数（如三刺激值计算公式、色调校正公式等）是基于特定研究的，对于不同区域可能需要进行验证和调整

通过上述代码，我们可以高效地计算长时间序列的Forel-Ule指数，为水体环境监测、水质评估等研究提供有力的数据支持。在这里插入图片描述