Google Earth Engine:对NDVI进行惠特克平滑算法进行长时序分析

2024-08-31 10:36

本文主要是介绍Google Earth Engine:对NDVI进行惠特克平滑算法进行长时序分析,希望对大家解决编程问题提供一定的参考价值,需要的开发者们随着小编来一起学习吧!

目录

简介

函数

ee.Array.identity(size)

Arguments:

Returns: Array

transpose(axis1, axis2)

Arguments:

Returns: Array

matrixMultiply(image2)

Arguments:

Returns: Image

matrixSolve(image2)

Arguments:

Returns: Image

arrayFlatten(coordinateLabels, separator)

Arguments:

Returns: Image

arrayReduce(reducer, axes, fieldAxis)

Arguments:

Returns: Image

代码

结果


简介

惠特克(GEE)平滑算法是一种用于时间序列预测的统计方法,特别适用于非线性、非平稳和非高斯的数据。该算法基于广义估计方程,通过最小化残差的平方和来拟合数据并找到最佳的平滑曲线。

GEE平滑算法的主要思想是在时间序列数据中引入一个平滑函数来描述数据的趋势和周期性变化。该平滑函数由一系列基函数的线性组合组成,其中每个基函数具有不同的频率和振幅。通过调整基函数的权重,可以得到最佳的平滑曲线,以最大程度地拟合数据。

在实际应用中,GEE平滑算法通常与其他统计方法结合使用,例如自回归移动平均模型(ARIMA)或指数平滑法。通过将GEE平滑算法与其他方法相结合,可以进一步提高时间序列的预测准确度和稳定性。

总的来说,GEE平滑算法是一种针对非线性、非平稳和非高斯数据的时间序列预测方法,通过引入一个平滑函数来描述数据的趋势和周期性变化,以最大程度地拟合数据。它在实际应用中通常与其他统计方法结合使用,以进一步提高预测的准确度和稳定性。

函数

ee.Array.identity(size)

Creates a 2D identity matrix of the given size.

创建一个给定大小的二维标识矩阵。

Arguments:

size (Integer):

The length of each axis.

Returns: Array

transpose(axis1axis2)

Transposes two dimensions of an array.

平移数组的两个维度。

Arguments:

this:array (Array):

Array to transpose.

axis1 (Integer, default: 0):

First axis to swap.

axis2 (Integer, default: 1):

Second axis to swap.

Returns: Array

matrixMultiply(image2)

Returns the matrix multiplication A * B for each matched pair of bands in image1 and image2. If either image1 or image2 has only 1 band, then it is used against all the bands in the other image. If the images have the same number of bands, but not the same names, they're used pairwise in the natural order. The output bands are named for the longer of the two inputs, or if they're equal in length, in image1's order. The type of the output pixels is the union of the input types.

返回图像 1 和图像 2 中每对匹配波段的矩阵乘法 A * B。如果图像 1 或图像 2 中只有一个波段,则该波段将与另一幅图像中的所有波段相对应。如果图像中的条带数量相同,但名称不同,则按自然顺序成对使用。输出波段以两个输入波段中较长的一个命名,如果两个输入波段长度相等,则按图像 1 的顺序命名。输出像素的类型是输入类型的组合。

Arguments:

this:image1 (Image):

The image from which the left operand bands are taken.

image2 (Image):

The image from which the right operand bands are taken.

Returns: Image

matrixSolve(image2)

Solves for x in the matrix equation A * x = B, finding a least-squares solution if A is overdetermined for each matched pair of bands in image1 and image2. If either image1 or image2 has only 1 band, then it is used against all the bands in the other image. If the images have the same number of bands, but not the same names, they're used pairwise in the natural order. The output bands are named for the longer of the two inputs, or if they're equal in length, in image1's order. The type of the output pixels is the union of the input types.

求解矩阵方程 A * x = B 中的 x,如果 A 对图像 1 和图像 2 中每对匹配的波段都是过确定的,则找到最小二乘法解。如果图像 1 或图像 2 中只有一个波段,则使用该波段与另一幅图像中的所有波段进行比对。如果图像中的波段数相同,但名称不相同,则按自然顺序成对使用。输出波段以两个输入波段中较长的一个命名,如果两个输入波段长度相等,则按图像 1 的顺序命名。输出像素的类型是输入类型的组合。

Arguments:

this:image1 (Image):

The image from which the left operand bands are taken.

image2 (Image):

The image from which the right operand bands are taken.

Returns: Image

arrayFlatten(coordinateLabels, separator)

Converts a single-band image of equal-shape multidimensional pixels to an image of scalar pixels, with one band for each element of the array.

将等形多维像素的单波段图像转换为标量像素图像,阵列中的每个元素对应一个波段。

Arguments:

this:image (Image):

Image of multidimensional pixels to flatten.

coordinateLabels (List):

Name of each position along each axis. For example, 2x2 arrays with axes meaning 'day' and 'color' could have labels like [['monday', 'tuesday'], ['red', 'green']], resulting in band names'monday_red', 'monday_green', 'tuesday_red', and 'tuesday_green'.

separator (String, default: "_"):

Separator between array labels in each band name.

Returns: Image

arrayReduce(reducer, axes, fieldAxis)

Reduces elements of each array pixel.

减少每个阵列像素的元素。

Arguments:

this:input (Image):

Input image.

reducer (Reducer):

The reducer to apply.

axes (List):

The list of array axes to reduce in each pixel. The output will have a length of 1 in all these axes.

fieldAxis (Integer, default: null):

The axis for the reducer's input and output fields. Only required if the reducer has multiple inputs or outputs.

Returns: Image

代码

//加载研究区
var geometry = /* color: #d63000 *//* displayProperties: [{"type": "rectangle"}] */ee.Geometry.Polygon([[[113.44773227683973, 38.6708907304602],[113.44773227683973, 38.64783484482313],[113.47588474266004, 38.64783484482313],[113.47588474266004, 38.6708907304602]]], null, false);// 将 qa 位图像转换为标志的辅助函数
function extractBits(image, start, end, newName) {// 计算我们需要提取的比特。var pattern = 0;for (var i = start; i <= end; i++) {pattern += Math.pow(2, i);}// 返回提取的质量保证位的单波段图像,并为该波段命名。return image.select([0], [newName]).bitwiseAnd(pattern).rightShift(start);
}// 在输入矩阵上获取指定阶次的差分矩阵的函数。将矩阵和阶次作为参数
function getDifferenceMatrix(inputMatrix, order){var rowCount = ee.Number(inputMatrix.length().get([0]));var left = inputMatrix.slice(0,0,rowCount.subtract(1));var right = inputMatrix.slice(0,1,rowCount);if (order > 1 ){return getDifferenceMatrix(left.subtract(right), order-1)}return left.subtract(right);
};// 将数组图像解包为图像和波段
// 以数组图像、图像 ID 列表和乐队名称列表为参数
function unpack(arrayImage, imageIds, bands){function iter(item, icoll){function innerIter(innerItem, innerList){return ee.List(innerList).add(ee.String(item).cat("_").cat(ee.String(innerItem)))}var temp = bands.iterate(innerIter, ee.List([]));return ee.ImageCollection(icoll).merge(ee.ImageCollection(ee.Image(arrayImage).select(temp,bands).set("id",item)))}var imgcoll  = ee.ImageCollection(imageIds.iterate(iter, ee.ImageCollection([])));return imgcoll}// 用于计算回归结果的反对数比率并转换回百分比单位的函数
function inverseLogRatio(image) {var bands = image.bandNames();var t = image.get("system:time_start");var ilrImage = ee.Image(100).divide(ee.Image(1).add(image.exp())).rename(bands);return ilrImage.set("system:time_start",t);
}function whittakerSmoothing(imageCollection, isCompositional, lambda){// 快速配置以设置默认值if (isCompositional === undefined || isCompositional !==true) isCompositional = false;if (lambda === undefined ) lambda = 10;// 程序启动  var ic = imageCollection.map(function(image){var t = image.get("system:time_start");return image.toFloat().set("system:time_start",t);});var dimension = ic.size();var identity_mat = ee.Array.identity(dimension);var difference_mat = getDifferenceMatrix(identity_mat,3);var difference_mat_transpose = difference_mat.transpose();var lamda_difference_mat = difference_mat_transpose.multiply(lambda);var res_mat = lamda_difference_mat.matrixMultiply(difference_mat);var hat_matrix = res_mat.add(identity_mat);// 备份原始数据var original = ic;// 获取原始图像属性var properties = ee.List(ic.iterate(function(image, list){return ee.List(list).add(image.toDictionary());},[]));var time = ee.List(ic.iterate(function(image, list){return ee.List(list).add(image.get("system:time_start"));},[]));// 如果数据是合成的// 计算图像在 0 到 100 之间的对比率。首先// 夹在 delta 和 100-delta 之间,其中 delta 是一个很小的正值。if (isCompositional){ic = ic.map(function(image){var t = image.get("system:time_start");var delta = 0.001;var bands = image.bandNames();image = image.clamp(delta,100-delta);image = (ee.Image.constant(100).subtract(image)).divide(image).log().rename(bands);return image.set("system:time_start",t);});}var arrayImage = original.toArray();var coeffimage = ee.Image(hat_matrix);var smoothImage = coeffimage.matrixSolve(arrayImage);var idlist = ee.List(ic.iterate(function(image, list){return ee.List(list).add(image.id());},[]));var bandlist = ee.Image(ic.first()).bandNames();var flatImage = smoothImage.arrayFlatten([idlist,bandlist]);var smoothCollection = ee.ImageCollection(unpack(flatImage, idlist, bandlist));if (isCompositional){smoothCollection = smoothCollection.map(inverseLogRatio);}// 通过添加后缀fitted获得新的乐队名称var newBandNames = bandlist.map(function(band){return ee.String(band).cat("_fitted")});// 重新命名平滑图像中的波段smoothCollection = smoothCollection.map(function(image){return ee.Image(image).rename(newBandNames)});// 一个非常笨的方法,可以flatten谷歌地球引擎生成的 ID,这样就可以将两张图片合并为图表了var dumbimg = arrayImage.arrayFlatten([idlist,bandlist]);var dumbcoll = ee.ImageCollection(unpack(dumbimg,idlist, bandlist));var outCollection = dumbcoll.combine(smoothCollection);var outCollectionProp = outCollection.iterate(function(image,list){var t = image.get("system:time_start")return ee.List(list).add(image.set(properties.get(ee.List(list).size())));},[]);var outCollectionProp = outCollection.iterate(function(image,list){return ee.List(list).add(image.set("system:time_start",time.get(ee.List(list).size())));},[]);var residue_sq = smoothImage.subtract(arrayImage).pow(ee.Image(2)).divide(dimension);var rmse_array = residue_sq.arrayReduce(ee.Reducer.sum(),[0]).pow(ee.Image(1/2));var rmseImage = rmse_array.arrayFlatten([["rmse"],bandlist]);return [ee.ImageCollection.fromImages(outCollectionProp), rmseImage];
}var ndvi =ee.ImageCollection("NOAA/VIIRS/001/VNP13A1").select('NDVI').filterDate("2019-01-01","2019-12-31");
// 去除屏蔽像素
ndvi = ndvi.map(function(img){return img.unmask(ndvi.mean())});var ndvi =  whittakerSmoothing(ndvi)[0];// 添加图表
print(ui.Chart.image.series(ndvi.select(['NDVI', 'NDVI_fitted']), geometry, ee.Reducer.mean(), 500).setSeriesNames(['NDVI', 'NDVI_fitted']).setOptions({title: 'smoothed',lineWidth: 1,pointSize: 3,
}));

结果

这篇关于Google Earth Engine:对NDVI进行惠特克平滑算法进行长时序分析的文章就介绍到这儿,希望我们推荐的文章对编程师们有所帮助!



http://www.chinasem.cn/article/1123663

相关文章

不懂推荐算法也能设计推荐系统

本文以商业化应用推荐为例,告诉我们不懂推荐算法的产品,也能从产品侧出发, 设计出一款不错的推荐系统。 相信很多新手产品,看到算法二字,多是懵圈的。 什么排序算法、最短路径等都是相对传统的算法(注:传统是指科班出身的产品都会接触过)。但对于推荐算法,多数产品对着网上搜到的资源,都会无从下手。特别当某些推荐算法 和 “AI”扯上关系后,更是加大了理解的难度。 但,不了解推荐算法,就无法做推荐系

康拓展开(hash算法中会用到)

康拓展开是一个全排列到一个自然数的双射(也就是某个全排列与某个自然数一一对应) 公式: X=a[n]*(n-1)!+a[n-1]*(n-2)!+...+a[i]*(i-1)!+...+a[1]*0! 其中,a[i]为整数,并且0<=a[i]<i,1<=i<=n。(a[i]在不同应用中的含义不同); 典型应用: 计算当前排列在所有由小到大全排列中的顺序,也就是说求当前排列是第

性能分析之MySQL索引实战案例

文章目录 一、前言二、准备三、MySQL索引优化四、MySQL 索引知识回顾五、总结 一、前言 在上一讲性能工具之 JProfiler 简单登录案例分析实战中已经发现SQL没有建立索引问题,本文将一起从代码层去分析为什么没有建立索引? 开源ERP项目地址:https://gitee.com/jishenghua/JSH_ERP 二、准备 打开IDEA找到登录请求资源路径位置

csu 1446 Problem J Modified LCS (扩展欧几里得算法的简单应用)

这是一道扩展欧几里得算法的简单应用题,这题是在湖南多校训练赛中队友ac的一道题,在比赛之后请教了队友,然后自己把它a掉 这也是自己独自做扩展欧几里得算法的题目 题意:把题意转变下就变成了:求d1*x - d2*y = f2 - f1的解,很明显用exgcd来解 下面介绍一下exgcd的一些知识点:求ax + by = c的解 一、首先求ax + by = gcd(a,b)的解 这个

综合安防管理平台LntonAIServer视频监控汇聚抖动检测算法优势

LntonAIServer视频质量诊断功能中的抖动检测是一个专门针对视频稳定性进行分析的功能。抖动通常是指视频帧之间的不必要运动,这种运动可能是由于摄像机的移动、传输中的错误或编解码问题导致的。抖动检测对于确保视频内容的平滑性和观看体验至关重要。 优势 1. 提高图像质量 - 清晰度提升:减少抖动,提高图像的清晰度和细节表现力,使得监控画面更加真实可信。 - 细节增强:在低光条件下,抖

【数据结构】——原来排序算法搞懂这些就行,轻松拿捏

前言:快速排序的实现最重要的是找基准值,下面让我们来了解如何实现找基准值 基准值的注释:在快排的过程中,每一次我们要取一个元素作为枢纽值,以这个数字来将序列划分为两部分。 在此我们采用三数取中法,也就是取左端、中间、右端三个数,然后进行排序,将中间数作为枢纽值。 快速排序实现主框架: //快速排序 void QuickSort(int* arr, int left, int rig

【Prometheus】PromQL向量匹配实现不同标签的向量数据进行运算

✨✨ 欢迎大家来到景天科技苑✨✨ 🎈🎈 养成好习惯,先赞后看哦~🎈🎈 🏆 作者简介:景天科技苑 🏆《头衔》:大厂架构师,华为云开发者社区专家博主,阿里云开发者社区专家博主,CSDN全栈领域优质创作者,掘金优秀博主,51CTO博客专家等。 🏆《博客》:Python全栈,前后端开发,小程序开发,人工智能,js逆向,App逆向,网络系统安全,数据分析,Django,fastapi

poj 3974 and hdu 3068 最长回文串的O(n)解法(Manacher算法)

求一段字符串中的最长回文串。 因为数据量比较大,用原来的O(n^2)会爆。 小白上的O(n^2)解法代码:TLE啦~ #include<stdio.h>#include<string.h>const int Maxn = 1000000;char s[Maxn];int main(){char e[] = {"END"};while(scanf("%s", s) != EO

业务中14个需要进行A/B测试的时刻[信息图]

在本指南中,我们将全面了解有关 A/B测试 的所有内容。 我们将介绍不同类型的A/B测试,如何有效地规划和启动测试,如何评估测试是否成功,您应该关注哪些指标,多年来我们发现的常见错误等等。 什么是A/B测试? A/B测试(有时称为“分割测试”)是一种实验类型,其中您创建两种或多种内容变体——如登录页面、电子邮件或广告——并将它们显示给不同的受众群体,以查看哪一种效果最好。 本质上,A/B测

秋招最新大模型算法面试,熬夜都要肝完它

💥大家在面试大模型LLM这个板块的时候,不知道面试完会不会复盘、总结,做笔记的习惯,这份大模型算法岗面试八股笔记也帮助不少人拿到过offer ✨对于面试大模型算法工程师会有一定的帮助,都附有完整答案,熬夜也要看完,祝大家一臂之力 这份《大模型算法工程师面试题》已经上传CSDN,还有完整版的大模型 AI 学习资料,朋友们如果需要可以微信扫描下方CSDN官方认证二维码免费领取【保证100%免费