GEE：应用遥感影像时空插值技术的实践（插值填补去云空洞）

作者：_养乐多_

本文介绍了几个 Google Earth Engine (GEE) 平台中常用的处理遥感数据中的缺失值代码片段，这些代码可以用于在时间序列中对遥感图像进行线性插值，提供更加连续和完整的时间序列。

第一段代码是一个线性插值函数，它能够在一个时间序列图像集合中对缺失数据进行插值。该函数使用了 GEE 平台中的图像集合操作函数和图像合成函数，并采用了线性插值方法对缺失值进行处理，最后返回一个插值后的图像集合。(该代码可以应用与去云后通过时间序列数据进行插值填补影像空洞）

第二段代码是一个掩模替换函数，它的作用是将原始影像中的某些值（如缺失值）用一个新值替换。它使用了 GEE 平台中的图像操作函数，并使用掩膜对影像进行处理，最后返回替换后的影像。

第三段代码是一个添加时间属性的函数，它会在影像中添加一个时间波段，用于后续时间序列分析。该函数同样使用了 GEE 平台中的图像操作函数和元数据属性获取函数。

最后一段代码是一个提取质量控制 (QA) 位的函数，它能够从一个多波段图像中提取指定位置的 QA 位数据。它使用了位运算函数和移位函数，并使用了 GEE 平台中的图像操作函数。

一、线性插值函数

二、掩模替换函数

三、添加时间属性的函数

四、提取质量控制 (QA) 位函数

一、线性插值函数

//线性插值填补无云像素
function linearInterp(imgcol, frame, nodata, bandname){frame  = frame  || 32;nodata = nodata || 0;var time   = 'system:time_start';imgcol = imgcol.map(addTimeBand);var maxDiff = ee.Filter.maxDifference(frame * (1000*60*60*24), time, null, time);var cond    = {leftField:time, rightField:time};// 待插影像之前的所有影像// Images before, 按时间升序排序 (最近的是排序后最后一个)var f2 = ee.Filter.and(maxDiff, ee.Filter.greaterThanOrEquals(cond));var c2 = ee.Join.saveAll({matchesKey:'before', ordering:time, ascending:true}).apply(c1, imgcol, f2);// 待插影像之后的所有影像// Images after 按时间降序排列 (最近的是排序后最后一个)var f1 = ee.Filter.and(maxDiff, ee.Filter.lessThanOrEquals(cond));var c1 = ee.Join.saveAll({matchesKey:'after', ordering:time, ascending:false}).apply(imgcol, imgcol, f1);var interpolated = ee.ImageCollection(c2.map(function(img) {img = ee.Image(img);var before = ee.ImageCollection.fromImages(ee.List(img.get('before'))).mosaic();var after  = ee.ImageCollection.fromImages(ee.List(img.get('after'))).mosaic();img = img.set('before', null).set('after', null);//替换空值，确认线性Interp有结果before = replace_mask(before, after, nodata);after  = replace_mask(after , before, nodata);//计算图像的时间之间的比率var x1 = before.select('time').double();var x2 = after.select('time').double();var now = ee.Image.constant(img.date().millis()).double();//待插影像时间减去前面的合成时间，除以后面的合成影像时间和前面的影像时间的差var ratio = now.subtract(x1).divide(x2.subtract(x1));  //x1 = x2时，为零//计算插值图像before = before.select(0); //移除时间波段after  = after.select(0);img    = img.select(0); var interp = after.subtract(before).multiply(ratio).add(before).toFloat();var qc = img.mask().not().rename('qc');interp = replace_mask(img, interp, nodata).rename(bandname);//return interp.addBands(qc).copyProperties(img, img.propertyNames());return interp.copyProperties(img, img.propertyNames());}));return interpolated;
}

这段代码实现了对图像集合进行线性插值，生成时间连续的图像序列，从而对时间序列进行进一步分析和处理。下面是代码的详细介绍：

linearInterp(imgcol, frame, nodata, bandname)：该函数接受四个参数，分别为：

imgcol：图像集合，即待插值的图像序列。
frame：时间窗口，即选择多少天内的图像插值，默认值为32。
nodata：缺失值，即需要插值的像素点的默认值，默认值为0。
bandname：插值后的波段名称，默认值为原始波段名称。

addTimeBand(img)：该函数用于添加时间属性，即为图像添加一个名为“time”的新波段，该波段的值为图像的起始时间。
ee.Filter.maxDifference(frame * (1000*60*60*24), time, null, time)：该语句用于生成一个时间过滤器，用于对图像集合进行分组，以便进行插值操作。frame * (1000*60*60*24) 表示时间差异的阈值，单位为毫秒。
ee.Join.saveAll({matchesKey:'after', ordering:time, ascending:false}).apply(imgcol, imgcol, f1)：该语句用于将图像集合按照时间降序排序，并为每个图像找到它之后的最近的图像。在该语句中，matchesKey 表示保存关联图像的键值，ordering 表示排序关键字，ascending 表示是否升序排序。
ee.Join.saveAll({matchesKey:'before', ordering:time, ascending:true}).apply(c1, imgcol, f2)：该语句用于将图像集合按照时间升序排序，并为每个图像找到它之前的最近的图像。在该语句中，matchesKey、ordering、ascending 的含义与上述相同。
interpolated：该语句是一个图像集合，其中的每个图像都是通过线性插值生成的新图像。具体来说，对于两个相邻的图像 before 和 after，首先用 replace_mask 函数替换缺失值，并计算出它们之间的时间比率，然后根据时间比率，通过线性插值生成一张新图像 interp。最后，将新图像的波段名称重命名为 bandname，并复制原始图像的属性。

综上所述，该代码实现了对图像序列进行线性插值，并返回一系列时间连续的图像序列，为后续的时间序列分析提供了重要的基础。

二、掩模替换函数

//掩模替换
function replace_mask(img, newimg, nodata) {nodata   = nodata || 0;var mask = img.mask();img = img.unmask(nodata);img = img.where(mask.not(), newimg);img = img.updateMask(img.neq(nodata));return img;
}

这段代码实现的是掩模替换的功能。该函数的输入参数包括一个原始图像img、一个用于替换无效数据的新图像newimg和一个可选的无效数据值nodata。输出为经过掩模替换后的图像。

函数中的第一行nodata = nodata || 0表示如果未传入无效数据值nodata，则默认值为0。
接下来，该函数通过img.mask()获取了原始图像的掩模（mask），然后通过img.unmask(nodata)将原始图像的无效数据（即值等于nodata的像素）转换为NaN，以便后续处理。
然后，函数调用img.where(mask.not(), newimg)，该函数会将掩模中值为false的像素替换为新图像newimg中对应像素的值。这个操作就是所谓的掩模替换。
最后，函数调用img.updateMask(img.neq(nodata))，该函数会将原始图像中值不等于nodata的像素重新赋上掩模，即原始图像中的有效数据区域。

最终函数返回经过掩模替换和更新掩模后的图像。

三、添加时间属性的函数

//添加时间属性
function addTimeBand(img) {/** 确保掩模一致 */var mask = img.mask();var time = img.metadata('system:time_start').rename("time").mask(mask);return img.addBands(time);
}

这段代码实现的是为图像添加时间属性的功能。该函数的输入参数为一个图像img，输出为添加时间属性后的图像。

函数中的第一步操作是获取图像的掩模，以确保添加的时间属性的掩模和原始图像一致。具体而言，函数通过img.mask()获取了原始图像的掩模（mask）。
接下来，函数通过img.metadata('system:time_start')获取原始图像的起始时间，即该图像的时间属性。获取的时间属性被重命名为time，然后通过.mask(mask)将其掩模与原始图像的掩模一致。这样就保证了时间属性的掩模和原始图像的掩模一致。
最后，函数通过img.addBands(time)将时间属性作为一维带（band）添加到原始图像中，从而得到了一个带有时间属性的新图像。

这个函数通常在处理遥感影像时使用，由于遥感影像一般都包含时间属性，因此为遥感影像添加时间属性可以更好地进行时间序列分析，例如通过时间属性对影像进行分类、监测变化等。

四、提取质量控制 (QA) 位函数

//返回提取的QA位的单波段图像
var getQABits = function(image, start, end) {var pattern = 0;for (var i = start; i <= end; i++) {pattern += 1 << i;}return image.bitwiseAnd(pattern).rightShift(start);
};

这段代码实现的是从多波段图像中提取QA（Quality Assurance）位的单波段图像的功能。该函数的输入参数包括一个多波段图像image，一个表示QA位的起始位置start，和一个表示QA位的终止位置end。输出为提取出来的QA位的单波段图像。

函数中的核心操作是通过位运算提取QA位。

首先，函数定义了一个变量pattern，并将其初始化为0。然后，函数通过循环遍历从start到end的所有位数，将每个位的二进制表示中的1向左移动相应的位数，最后将它们相加，得到一个表示QA位的二进制模式。具体而言，这个过程就是将一个二进制的"1"左移若干位（位数由i决定），然后与变量pattern相加。
接下来，函数通过image.bitwiseAnd(pattern)将二进制模式应用到输入的多波段图像image上，得到一个新的多波段图像，其中仅保留了QA位上的像素值。
最后，函数通过.rightShift(start)将提取出来的QA位向右移动start位，得到一个单波段图像，即为提取出来的QA位的单波段图像。

该函数通常在遥感影像处理中使用，例如对Landsat影像的QA信息进行处理，以确定哪些像素是有用的（如云、阴影、水等）并排除掉无用的像素。