ESA SNAP水色反演算法在GEE中的实现

最新推荐文章于 2024-11-18 11:39:54 发布
原创 最新推荐文章于 2024-11-18 11:39:54 发布 · 5.3k 阅读 · 48 ·
CC 4.0 BY-SA版权
版权声明：本文为博主原创文章，遵循 CC 4.0 BY-SA 版权协议，转载请附上原文出处链接和本声明。
文章标签：

#google earth #gis

该文章已生成可运行项目，

1. ESA SNAP 与FUI水色指数

SNAP (Sentinel Application Platform) 哨兵数据应用平台，是欧空局针对哥白尼计划中哨兵数据任务开发的免费开源的软件。详细了解点击这里。

SNAP 软件中集成了多种算法的实现，我要介绍的是 FUI 水色指数，它被集成在 Optical / Thematic Water Processing / FU Classification 模块中。与水温，盐度和透明度一起，水体颜色属于最古老的时间序列水质数据之一 ，水中的三种主要光学物质—叶绿素、非藻类悬浮物和 CDOM 以及水分子本身的吸收和散射作用共同决定了水体呈现出的颜色，基于 Forel-Ule 比色表将水体颜色从深蓝 色到红棕色划分为 21 个等级。研究表明，由于涵盖非常广泛的自然水体光学特征，FUI 水色指数在区域和全球长时间尺度上监测反演水体水质方面具有强大潜力和优势。

image-20211121165440539

基于卫星图像的 FUI 水色指数和色度角度的系统研究始于 2012 年欧洲Citclops 项目。在该项目中，Wernand 和 Novoa 等逐步完善了 Forel-Ule 比色液的颜色光谱分析，构建了完整的 FUI 指数 CIE 色度坐标点。以下为基于卫星影像的 FUI 水色指数计算流程。

[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-6hrDKkGQ-1637637812526)(https://raw.githubusercontent.com/cyandee/Image/main/imgFUI.png)]

针对不同遥感影像，根据不同波段系数计算得到颜色三刺激值 XYZ，将其归一化为二维点坐标 （chrx，chry），计算其与等能白光点（x=y=1/3) 沿着 X 正轴方向逆时针旋转所得的夹角α，称其为色度角，色度角按照 FUI 指数查找表重分类为 FUI 水色指数的 21 个等级，完成 FUI 水色指数的计算（值得注意的是，不同传感器由于自身波段离散的特点，R、G、B 三波段中心波长与 CIE-XYZ 系统定义的 R(645)、G(555)、B(469)并不是一一对应，因此，多光谱图像计算所得色度角与高光谱波谱的积分结果相比必然会有一定的系统误差，需要进行色度角校正）。

上述算法在 ESA SNAP 中已经内置，直接输入影像运算可直接得到计算结果。而我们在对某湖库的水色状况进分析时，一般需要计算其长时间序列的水色指数，使用 SNAP 虽然可以便捷的计算单景影像的 FU 结果，但是长时序分析计算需要下载处理多景影像，本地计算耗时费力。因此，考虑将 ESA SNAP 中内置的算法移植到 GEE 平台中。 SNAP的源代码在 GitHub 开源，可以找到 FU 计算的核心代码：

package org.esa.s3tbx.fu;

import com.bc.ceres.core.Assert;

class FuAlgo {

    private static final double HUE_MIN = 45.0;
    private static final double HUE_MAX = 234.0;
    public static final byte MAX_FU_VALUE = 21;
    public static final byte MIN_FU_VALUE = 0;
	
    //等能白光点
    private static final double CONST_WHITE_POINT = 0.333333;
    //FUI指数查找表
    private static final double[] ANGLE_OF_TRANSITIONS = new double[]{
            232.0, 227.168, 220.977, 209.994, 190.779, 163.084, 132.999,
            109.054, 94.037, 83.346, 74.572, 67.957, 62.186, 56.435,
            50.665, 45.129, 39.769, 34.906, 30.439, 26.337, 22.741, 19.0, 19.0
    };
    private DominantLambdaLookup lambdaLookup;
    
    //三刺激值
    private double[] x3Factors;
    private double[] y3Factors;
    private double[] z3Factors;
    //色度角校正值
    private double[] polyCoeffs;

    FuAlgo(Instrument instrument, boolean includeDominantLambda) {
        x3Factors = instrument.getXFactors();
        y3Factors = instrument.getYFactors();
        z3Factors = instrument.getZFactors();
        polyCoeffs = instrument.getPolynomCoefficients();
        if (includeDominantLambda) {
            lambdaLookup = new DominantLambdaLookup();
        }
    }

    // just for testing
    FuAlgo() {
    }


    void setPolyCoeffs(double[] polyCoeffs) {
        this.polyCoeffs = polyCoeffs;
    }

    void setZ3Factors(double[] z3Factors) {
        this.z3Factors = z3Factors;
    }

    void setY3Factors(double[] y3Factors) {
        this.y3Factors = y3Factors;
    }

    void setX3Factors(double[] x3Factors) {
        this.x3Factors = x3Factors;
    }

    public FuResult compute(double[] spectrum) {
        final double x3 = getTristimulusValue(spectrum, x3Factors);
        final double y3 = getTristimulusValue(spectrum, y3Factors);
        final double z3 = getTristimulusValue(spectrum, z3Factors);

        final double denominator = x3 + y3 + z3;

        final double chrX = x3 / denominator;
        final double chrY = y3 / denominator;

        double hue = getHue(chrX, chrY);
        final double polyCorr;
        if (hue < HUE_MIN || hue > HUE_MAX) {
            hue = Double.NaN;
            polyCorr = Double.NaN;
        } else {
            final double hue100 = hue / 100;
            polyCorr = getPolyCorr(hue100, polyCoeffs);
        }

        FuResultImpl result = new FuResultImpl();
        result.x3 = x3;
        result.y3 = y3;
        result.z3 = z3;
        result.chrX = chrX;
        result.chrY = chrY;
        result.hue = hue;
        result.polyCorr = polyCorr;
        result.hueAngle = hue + polyCorr;
        result.fuValue = getFuValue(result.hueAngle);
        if (lambdaLookup != null) {
            result.domLambda = lambdaLookup.getDominantLambda(result.hueAngle);
        }
        return result;
    }

    static byte getFuValue(final double hueAngle) {
        if (Double.isNaN(hueAngle)) {
            return MIN_FU_VALUE;
        }
        for (byte i = 0; i < ANGLE_OF_TRANSITIONS.length; i++) {
            if (hueAngle > ANGLE_OF_TRANSITIONS[i]) {
                return i;
            }
        }
        return MAX_FU_VALUE;
    }

    double getTristimulusValue(double[] spectrum, double[] factors) {
        if (spectrum.length != factors.length) {
            throw new IllegalArgumentException("The spectrum must have equal length as factors.");
        }
        double summation = 0;
        for (int i = 0; i < spectrum.length; i++) {
            summation = (spectrum[i] * factors[i]) + summation;
        }
        return summation;
    }

    double getHue(double chrX, double chrY) {
        double atan2 = Math.atan2((chrY - CONST_WHITE_POINT), (chrX - CONST_WHITE_POINT));
        final double hue = (180 * atan2) / Math.PI;
        return hue < 0 ? hue + 360 : hue;
    }

    double getPolyCorr(double hue100, double[] constPolyHue) {
        Assert.argument(constPolyHue.length == 6, "constPolyHue.length == 6");
        double value = 0.0;
        for (int i = 0; i < constPolyHue.length; i++) {
            if ((i + 1) == constPolyHue.length) {
                value += constPolyHue[i];
                break;
            }
            value += constPolyHue[i] * Math.pow(hue100, 5 - i);
        }
        return value;
    }
}

2. FUI 反演算法在 GEE 中的实现

参考论文中针对不同传感器的参数，可在 GEE 中实现基于 Landsat 系列，Sentinel-2 的 FUI 水色反演。下图为不同传感器在计算三刺激值时所对应的波段系数。

image-20211122225849895

下图为不同传感器的色度角校正系数。

image-20211122230154749

根据论文中的波段系数，结合 SNAP 开源代码，可以方便的将 FUI 计算流程移植到 GEE 平台中，借助 GEE 平台海量的数据和强大的运算能力，轻松计算某地长时间序列水色结果。

var getFU = function(image){
  var X = image.expression('11.053 * b1 + 6.950 * b2 + 51.135 * b3 + 34.457 * b4',{
      b1:image.select('B1'),
      b2: image.select('B2'),
      b3: image.select('B3'),
      b4: image.select('B4')
    }).rename('X');
  
  var Y = image.expression('1.32 * b1 + 21.053 * b2 + 66.023 * b3 + 18.034 * b4',{
    b1:image.select('B1'),
    b2: image.select('B2'),
    b3: image.select('B3'),
    b4: image.select('B4')
  }).rename('Y');
  
  var Z = image.expression('58.038 * b1 + 34.931 * b2 + 2.606 * b3 + 0.016 * b4',{
    b1: image.select('B1'),
    b2: image.select('B2'),
    b3: image.select('B3'),
    b4: image.select('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 atan2 = chrx.subtract(0.3333).atan2((chry).subtract(0.3333)).rename('atan2');
  var hue_pre = atan2.multiply(ee.Number(180)).divide(Math.PI).rename('hue_pre');
  var hue = hue_pre < 0 ? (hue_pre + 360) : hue_pre.rename('hue');
  var a = hue.divide(100).rename('a');
  var cor = a.pow(5).multiply(-52.16).add(a.pow(4).multiply(373.81))
    .add(a.pow(3).multiply(-981.3)).add(a.pow(2).multiply(1134.19))
    .add(a.multiply(-533.61)).add(76.72).rename('cor');
  var hueAngle = hue.add(cor).rename('hueAngle');
  
  /*
    var angle_tran = ee.List([232.0, 227.168, 220.977, 209.994, 190.779, 
    163.084, 132.999,109.054, 94.037, 83.346, 74.572, 67.957, 62.186, 56.435,50.665, 
    45.129, 39.769, 34.906, 30.439, 26.337, 22.741, 19.0, 19.0]);
  */
  
  /*
    static byte getFuValue(final double hueAngle) {
        for (byte i = 0; i < ANGLE_OF_TRANSITIONS.length; i++) {
            if (hueAngle > ANGLE_OF_TRANSITIONS[i]) {
                return i;
            }
        }
        return MAX_FU_VALUE;
    }
  */
  var FU = a
    .where(hueAngle.gt(19.0).and(hueAngle.lte(22.741)),1)
    .where(hueAngle.gt(22.741).and(hueAngle.lte(26.337)),2)
    .where(hueAngle.gt(26.337).and(hueAngle.lte(30.439)),3)
    .where(hueAngle.gt(30.439).and(hueAngle.lte(34.906)),4)
    .where(hueAngle.gt(34.906).and(hueAngle.lte(39.769)),5)
    .where(hueAngle.gt(39.769).and(hueAngle.lte(45.129)),6)
    .where(hueAngle.gt(45.129).and(hueAngle.lte(50.665)),7)
    .where(hueAngle.gt(50.665).and(hueAngle.lte(56.435)),8)
    .where(hueAngle.gt(56.435).and(hueAngle.lte(62.186)),9)
    .where(hueAngle.gt(62.186).and(hueAngle.lte(67.957)),10)
    .where(hueAngle.gt(67.957).and(hueAngle.lte(74.572)),11)
    .where(hueAngle.gt(74.572).and(hueAngle.lte(83.346)),12)
    .where(hueAngle.gt(83.346).and(hueAngle.lte(94.037)),13)
    .where(hueAngle.gt(94.037).and(hueAngle.lte(109.054)),14)
    .where(hueAngle.gt(109.054).and(hueAngle.lte(132.999)),15)
    .where(hueAngle.gt(132.999).and(hueAngle.lte(163.084)),16)
    .where(hueAngle.gt(163.084).and(hueAngle.lte(190.779)),17)
    .where(hueAngle.gt(190.779).and(hueAngle.lte(209.994)),18)
    .where(hueAngle.gt(209.994).and(hueAngle.lte(220.977)),19)
    .where(hueAngle.gt(220.977).and(hueAngle.lte(227.168)),20)
    .where(hueAngle.gt(227.168),21).rename('FU');
  
  return image.addBands(FU);
}

以下为 FUI 指数查找表和不同等级 FUI 对应的 RGB数据和颜色代码：

[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-VGv0RLoJ-1637637812529)(https://raw.githubusercontent.com/cyandee/Image/main/imgFUI%E6%9F%A5%E6%89%BE%E8%A1%A8.png)]

FUIRGB十六进制颜色代码FUIRGB十六进制颜色代码
13388188#2158BC1214818296#94B660
249109197#316DC513165188118#A5BC75
350124187#327CBB14170184109#AAB86D
475128160#4B80A01517318195#ADB55F
586143150#568F9616168169101#A8A965
6109146152#6D92981717415992#AE9F5C
7105140134#698C861817916083#B3A053
8117158114#759E751917513868#AF8A44
912316684#7BA654201649610#A4600A
1012517456#7DAE3821164838#A45308
1114918269#95B645

附上之前在 GEE 里做的一个小试验：https://code.earthengine.google.com/7b0b97aa01a0f86dcf64e10e1dc4869f

另外，中科院王胜蕾的博士论文中提到了另一种色度角的计算方式和新的FUI指数查找表，并对色度角校正方法进行了详细的阐述，感兴趣的同学可自行探索。

参考论文地址

  1. https://www.mdpi.com/2072-4292/10/2/180
  2. https://www.mdpi.com/1424-8220/15/10/25663
  3. http://www.jeos.org/index.php/jeos_rp/article/view/13057
  4. https://journals.plos.org/plosone/article?id=10.1371/journal.pone.0063766
  5. https://kns-cnki-net-443.v.cnu.edu.cn/kcms/detail/detail.aspx?dbcode=CDFD&dbname=CDFDLAST2019&filename=1018995624.nh&uniplatform=NZKPT&v=8jR0NvYJM7ur0xgkduGfdTUjQK55ycqt2kmy2mmpMUjOJcydchiAcvB9pLactvP7

欢迎关注公众号：卫星环境应用
在这里插入图片描述

本文章已经生成可运行项目
GEE案例——基于landsat 8的计算特定区域（太湖和鄱阳湖）的 FUI（遥感水色指数）监测和可视化分析
此星光明博客
03-06 863
FUI（遥感水色指数，Fluorescence Utilization Index）是一种用于评估水体健康和生态状况的遥感指标。它通过分析水体的光谱特征，特别是水体中植物的光合作用和荧光信号，来提供有关水体质量的信息。
GEE案例——归一化差异水体指数丰水期、枯水期的水域面积和水深分析（青海湖为例）
此星光明博客
06-04 1054
水深反演是指利用遥感技术从航空或卫星平台上获取的数据来推断水体的深度信息。这种技术在海洋学、湖泊和河流的科学研究与管理中非常重要。以下是几种常用的水深反演方法：1. **光学遥感反演**：- 基于水体颜色和透明度的变化与水深的关系，使用光学遥感影像（如Landsat、Sentinel等）进行水深估算。- 通常适用于清澈的浅水环境，因为水体的光学特性会随着深度的增加而变化。2. **多光谱和高光谱遥感**：- 利用水体对不同波长光的吸收和散射特性，通过分析多光谱或高光谱数据来估算水深。
GEE利用哨兵2号数据进行水质监测
guaicaicexu的博客
11-18 1306
今天来简单分享下如何在GEE中利用哨兵2号数据进行水质监测总悬浮固体（TSS）总悬浮固体（TSS）是指尺寸超过2微米的水性颗粒。大部分总悬浮固体由无机材料组成；然而，藻类和细菌也可能被认为是TSS。总悬浮固体 (TSS) 过量存在对人体有害，在《清洁水法》中被列为常规污染物，它可以表示任何水样（从海水到废水）的质量。广义来讲的话，TSS可以是任何漂浮或“悬浮”在水中的东西，包括沙子、沉积物和浮游...
GEE专栏之水质反演 基于线下训练模型的在线叶绿素反演（叶绿素反演结果可视化、图例、白板底图、提取点位像元值并输出CSV）
Lbta的博客
06-25 4681
GEE专栏之水质反演 基于线下训练模型的在线预测 效果图 特色：GEE绘制图例、提取像元值、反演模型与GEE的融合。 该结果为洱海地区叶绿素水质反演结果，利用GEE ui绘制了图例，屏蔽了影像图层即增加了白色底图，提取监测点位像元值并输出CSV，代码部分注释掉了输出像元值，取消注释即可。可供参考 源码部分 var p280 = ee.Geometry.Point([100.16,25.86]) var p281 = ee.Geometry.Point([100.19,25.83]) var p283 =
GEE 反演叶绿素的app
zzme110的博客
10-12 546
GEE反演湖库叶绿素浓度
【文献阅读01】直接反演式水色遥感算法研究
qq_40784982的博客
03-06 1347
一、背景 文章题目：《直接反演式水色遥感算法研究》 作者：曾炼 2009 华中科技大学硕士论文 关键词：海洋水色遥感 大气修正算法 水色要素反演算法 神经网络算法 二、文章简要介绍 文章摘要： 介绍了主要的水色要素反演算法，用MODIS数据、MATLAB神经网络编程反演水色参数chla。 三、文章主要部分 大气校正的基本思想：I II类水体 给出了MODIS数据的光谱通道，并指出水色通道为8...
ESA 哨兵数据处理软件 SNAP
11-25
ESA（欧空局）提供的哨兵系列数据处理软件，可对哨兵系列数据、landsat系列数据、spot系列数据等光学遥感数据和SAR数据进行诸如几何校正、辐射定标、滤波、裁剪、镶嵌等一系列处理。
Google Earth Engine学习笔记（一）
qq_46610545的博客
03-27 1822
基于GEE平台，Landsat 8、Sentinel-2卫星，FUI水色指数计算
精选资源
SNAP -SLC数据双极化指导书.zip_esa snap_snap 哨兵一号_双极化_哨兵1号预处理_哨兵数据
07-14
对哨兵一号数据利用欧空局软件SNAP进行预处理的过程
精选资源
esa-snap_sentinel_windows-x64_7_0.exe
01-30
ESA SNAP（Sentinel Application Platform）是由欧洲空间局（ESA）开发的一款开源软件，主要用于处理和分析地球观测数据，特别是Sentinel系列卫星的数据。这个“esa-snap_sentinel_windows-x64_7_0.exe”文件是ESA ...
精选资源
SIF_叶绿素荧光SIF_GOME-2_叶绿素荧光反演算法_arrangedeq_GOME-2SIF_
10-01
在描述中提到的"arrangedeq"可能是指对原始数据进行预处理和排序，以适应反演算法的需求。 `SIF.py`文件很可能是实现这一反演算法的核心代码。在Python编程环境下，这个文件可能包含了数据读取、预处理、反演计算和...
esa-snap_sentinel_windows-x64_8_0.exe
06-11
esa-snap_sentinel_windows-x64_8_0.exe 欧空局官方哨兵数据处理工具、将JP2转为TIFF格式，是哨兵遥感影像下载后处理的不可少的一个插件。
GEE案例-基于landsat 8的计算特定区域（太湖和鄱阳湖）的 FUI（遥感水色指数）监测和可视化分析
最新发布
03-06
GEE案例——基于landsat 8的计算特定区域（太湖和鄱阳湖）的 FUI（遥感水色指数）监测和可视化分析 包含代码以及详细解答和全部代码
GEE随记（六）：Sentinel影像反演浊度
weixin_48012947的博客
04-22 980
Sentinel影像反演水体浊度并可视化
GEE反演近岸海域互花米草、红树林等
zzme110的博客
10-22 477
GEE滨海湿地分类识别
【项目实战】Python基于决策树多分类模型实现水色图像的水质评价
张陈亚的博客
01-08 5719
Python基于决策树多分类模型实现水色图像的水质评价
中科星图（案例）——基于Landsat8SR影像反演太湖遥感水色FUI指数
此星光明博客
02-21 2181
本文主要是利用gve云平台实现Landsat8 影像进行太湖FUI水色反演。基于Landsat 8 SR影像反演太湖遥感水色FUI指数的集体操作流程主要包括以下几个步骤：1. 数据获取与预处理：首先，需要从美国地质调查局（USGS）或其他可靠的数据源获取Landsat 8 SR影像数据。在太湖周边选择合适的日期的遥感影像，以确保太湖区域的相关信息完整。获取到影像数据后，需要进行预处理，如辐射定标、大气校正、几何校正等。
机器学习-04-基于水色图像的水质评价（包含代码和数据）
AI研习星球
11-30 3544
一般情况下，采集到的水样图片包含盛水容器，且容器颜色与水体颜色差异较大，同时水体位于图片中央，所以为了提取水色特征，就需要提取水样图片中央部分具有代表意义的图像，具体实施方式是提取水样图像中央101×101像素的图像。颜色直方图是最基本的颜色特征表示方法，它反映的是图像中颜色的组成分布即出现了哪些颜色以及各种颜色出现的概率。颜色直方图产生的特征维数一般大于颜色矩的特征维数，为了避免变量过多影响后续分类效果，选择基于颜色矩提取图像特征的方式，建立水样图像与反映该图像特征的数据信息之间的关系。
GEE学习篇：Landsat地表温度反演的方法
weixin_47787901的博客
01-15 8254
GEE平台地表温度反演方法
esa snap设置中文
03-03
要在ESA SNAP中设置中文界面，可以按照以下步骤进行操作： 1. 下载和安装ESA SNAP软件：首先，你需要从ESA SNAP官方网站（https://step.esa.int/main/download/）下载并安装最新版本的ESA SNAP软件。 2. 启动ESA ...
评论 14
成就一亿技术人!
拼手气红包6.0元
还能输入1000个字符
 
红包 添加红包
表情包 插入表情
表情包 代码片
 条评论被折叠 查看
添加红包

请填写红包祝福语或标题

个

红包个数最小为10个

元

红包金额最低5元

当前余额3.43元 前往充值 >
需支付：10.00元
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付元
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明：

1.余额是钱包充值的虚拟货币，按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载，可以购买VIP、付费专栏及课程。

余额充值