一套GEE数据插补平滑完整方案2：使用S-G和Whittaker滤波
KeyPoints
-
遥感数据缺失和平滑是一个问题。 -
数据缺失存在两种情况，一种是局部数据缺失，一种是整张影像缺失，影像后续分析。 -
数据平滑有助于去除噪声，使得数据更为清晰、准确。此外，遥感数据中常常包含着地形、云层等导致数据的波动较大，数据平滑可以改善图像质量。 -
对于数据处理，采取“先插补，再平滑”的策略。上期推文解决插补问题，这次解决平滑问题。 -
本文利用Savitzky-Golay滤波和Whittaker滤波解决上述两种情况。
S-G滤波
原理
S-G滤波（Savitzky-Golay滤波）是一种常用于数据平滑和拟合的数字滤波器。它可以在不降低信号分辨率的情况下去除信号中的噪声和粗糙部分，适用于各种数据类型，如光谱数据、生物信号、地球物理数据等。
S-G滤波的原理是基于局部多项式拟合的思想，其核心是对信号进行局部加权多项式拟合。在滤波过程中，通过对一段信号进行多项式拟合，将该段信号近似为一个多项式函数。在拟合过程中，多项式的阶数和拟合窗口大小是两个关键参数，它们决定了拟合的复杂度和拟合的精度。一般来说，阶数和窗口大小应该根据数据的特征进行调整，以达到最佳的滤波效果。
S-G滤波可以通过矩阵运算的方式进行实现，这种方法计算速度较快，并且可以使用线性代数库来实现。S-G滤波常用于对光谱数据进行平滑处理，也可以用于对其他类型的数据进行平滑和去噪处理。
公式推导
以拟合一个三次多项式为例，窗口是5，根据三次多项式公式：
那么这五个窗口的值可以代入进去，写成矩阵形式：
不妨把上式写为：
那么我们只需要求解系数A，对应a0-a3，就能算出所有平滑的Y了，根据线性代数最小二乘的推导，A的解为：
只要在GEE函数中实现上述矩阵运算就行了
代码
第一步，加载数据
var geometry = /* color: #d63000 */ee.Geometry.Point([105.77243064871553, 13.54490072286009])
var filterLAI = function(img){
var newimg = img.copyProperties(img, img.propertyNames());
return newimg
}
var start_date = ee.Date('2019-09-27');
var end_date = ee.Date('2022-11-18');
var latter_lai = ee.ImageCollection('MODIS/006/MCD15A3H').select('Lai').filter(ee.Filter.date('2019-09-27', '2022-11-18')).map(filterLAI)
// getting rid of masked pixels
var latter_lai = latter_lai.map(function(img){return img.unmask(latter_lai.mean())});
/****************************************************************************
// Add predictors for SG fitting, using date difference
// We prepare for order 3 fitting,
// but can be adapted to lower order fitting later on
****************************************************************************/
//to remove cloud
function mosaicByDate(imcol){
// convert the ImageCollection into List
var imlist = imcol.toList(imcol.size());
// print(imlist)
// Obtain the distinct image dates from the ImageCollection
var unique_dates = imlist.map(function(im){
return ee.Image(im).date().format("YYYY-MM-dd");
}).distinct();
// print(unique_dates);
// mosaic the images acquired on the same date
var mosaic_imlist = unique_dates.map(function(d){
d = ee.Date(d);
var im = imcol.filterDate(d, d.advance(1, "day")).mosaic();
//print(im)
// return the mosaiced same-date images and set the time properties
return im.set(
"system:time_start", d.millis(),
"system:id", d.format("YYYY-MM-dd")
);
});
return ee.ImageCollection(mosaic_imlist);
}
var modis_res = latter_lai
.map(function(img){
var dstamp = ee.Date(img.get('system:time_start'));
var ddiff = dstamp.difference(ee.Date(start_date), 'day');
img = img.select(['Lai']).set('date', dstamp)
.clip(table).set('footprint',table);
return img.addBands(ee.Image(1).toFloat().rename('constant'))
.addBands(ee.Image(ddiff).toFloat().rename('t'))
.addBands(ee.Image(ddiff).pow(ee.Image(2)).toFloat().rename('t2'))
.addBands(ee.Image(ddiff).pow(ee.Image(3)).toFloat().rename('t3'));
});
print("modis_res",modis_res);
Map.addLayer(modis_res.max().select('Lai').clip(table),
{'min':0,'max':30,'palette':['black','green']}, 'modis_res');
第二步，建立多项式：
这里设置window_size为9，多项式次数order为3，结合平滑效果请自行修改。
/****************************************************************************
// Step 2: Set up Savitzky-Golay smoothing
****************************************************************************/
// set the window size
var window_size = 9;
var half_window = (window_size - 1)/2;
// Define the axes of variation in the collection array.
var imageAxis = 0;
var bandAxis = 1;
// Set polynomial order
var order = 3;
var coeffFlattener = [['constant', 'x', 'x2', 'x3']];
var indepSelectors = ['constant', 't', 't2', 't3'];
// Convert the collection to an array.
var array = modis_res.toArray();
// For the remainder, use modis_res as a list of images
modis_res = modis_res.toList(modis_res.size());
var runLength = ee.List.sequence(0, modis_res.size().subtract(window_size));
// Run the SG solver over the series, and return the smoothed image version
var sg_series = runLength.map(function(i) {
var ref = ee.Image(modis_res.get(ee.Number(i).add(half_window)));
return getLocalFit(i).multiply(ref.select(indepSelectors))
.reduce(ee.Reducer.sum()).copyProperties(ref);
});
print("sg_series",sg_series);
第三步，求解，利用matrixSolve求解矩阵
// Solve
function getLocalFit(i){
// Get a slice corresponding to the window_size of the SG smoother
var subarray = array.arraySlice(imageAxis, ee.Number(i).int(), ee.Number(i).add(window_size).int());
var predictors = subarray.arraySlice(bandAxis, 1, 1 + order + 1);
var response = subarray.arraySlice(bandAxis, 0, 1); // NDVI
var coeff = predictors.matrixSolve(response);
coeff = coeff.arrayProject([0]).arrayFlatten(coeffFlattener);
return coeff;
}
第四步，根据求解多项式的系数，生成结果：
// Part 3: Generate some output
// 3A. Get an example original image and its SG-ed version
var NDVI = (ee.Image(modis_res.get(50)).select('Lai'));
var NDVIsg = (ee.Image(sg_series.get(50))); // half-window difference in index
Map.addLayer(NDVI,{min:0,max:30},"NDVI");
Map.addLayer(NDVIsg,{min:0,max:30},"NDVIsg");
最后结果如图所示：
完整的代码可以后台回复获取，见后文
Whittaker滤波
原理
Whittaker滤波是一种基于时间序列的平滑方法，通过对时间序列进行加权平均来平滑信号。其原理基于以下两个假设：
-
时间序列中的数据点存在一个随机误差，该误差是服从正态分布的，并且误差是独立同分布的。 -
时间序列中数据点之间并没有显著的周期性或趋势性变化。
Whittaker滤波的核心思想是对每个数据点进行加权平均，权重系数是根据时间和空间距离计算的。具体来说，Whittaker滤波算法将每个数据点表示为一个二维网格（即空间-时间网格），并且使用高斯核函数来计算每个数据点在该网格内的权重。然后，通过加权平均来计算每个数据点的平滑值。
Whittaker滤波可以应用于多种不同类型的时间序列数据，包括光谱数据、图像数据和地理空间数据等。它被广泛应用于信号处理、数据降噪和数据插值等领域。
根据Nishanta等人发表在Remote Sensing上的文章：
Whittaker能用于GEE，并产生比SG更好的效果
代码
Nishanta等人在论文附录公布了代码，然而他的代码只是对于土地覆盖，因此对其进行了一些修改，使之适用于NDVI、LAI等变量：
由于篇幅有限，这里不再展示完整代码，仅展示结果：
S-G与Whittaker对比
可以看到Whittaker滤波总体上比三次多项式S-G好。
在细节信息上，S-G能较大程度保留一些细节信号（但这些细节信号更有可能是异常信号，如云量，噪声等等），反而产生错误的结果。
在遥感数据，植被信号等，Whittaker滤波可能比S-G滤波更适合。
Whittaker通过对原始数据进行重构来消除噪声和不规则性。它的效果通常取决于所选择的平滑因子和曲线拟合程度。
与之相比，S-G滤波是使用多项式拟合来平滑信号，并且可以控制拟合的阶数和窗口大小，该方法较为粗糙。
最近越来越多研究也关注改进后的Whittaker滤波，而很少用S-G处理这些遥感信号。
如Kong等人发表在ISPRS的研究，提出了一种加权后的Whittaker方法，并在空间尺度上设置动态参数lambda（称为wWHd方法）。显著改善了精度。
效果如图：
如有相关需求，可以引用系列文章：
Kong, D., Zhang, Y., Gu, X., & Wang, D. (2019). A robust method for reconstructing global MODIS EVI time series on the Google Earth Engine. ISPRS Journal of Photogrammetry and Remote Sensing, 155, 13-24.
Kong, D., McVicar, T. R., Xiao, M., Zhang, Y., Peña‐Arancibia, J. L., Filippa, G., ... & Gu, X. (2022). phenofit: An R package for extracting vegetation phenology from time series remote sensing. Methods in Ecology and Evolution, 13(7), 1508-1527.
代码
上述S-G代码和Whittaker代码可以后台回复【GEE滤波】，wWHd相关请参考孔老师GitHub，wWHd：https://github.com/gee-hydro/gee_whittaker
S-G：https://code.earthengine.google.com/c73b17534179363adb578558d46f78a7
Whittaker：https://code.earthengine.google.com/67448f422e8a55132e87cff1653a67ee
本文由 mdnice 多平台发布
994
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
