Google Earth Engine谷歌地球引擎图像自动地理配准与空间位置校正

  本文主要对 GEE 中栅格图像的 地理配准 （空间坐标位置校正）操作加以介绍。

  本文是谷歌地球引擎（Google Earth Engine，GEE）系列教学文章的第十四篇，更多 GEE 文章请参考专栏：GEE学习与应用_疯狂学习GIS的博客-CSDN博客。

  在之前的两篇文章基于ENVI实现遥感影像栅格图层手动地理配准 与ENVI自动产生地面控制点GCP实现栅格遥感影像的自动地理配准 中，我们介绍了基于 ENVI 软件的遥感影像 手动、自动 地理配准方法；本文则介绍在 GEE 中， 自动 实现遥感影像间地理配准的方法。

  在 GEE 中实现遥感影像地理配准，总体分为两个步骤。首先，需要设置一景“ 位移影像 ”，其包含三个波段：首先前两个波段为 dx 与 dy ——二者分别为图像中每一个像元在 x 和 y 方向的偏移量；第三个波段为上述偏移量的置信度。其次，将设置好的“位移影像”应用于需要进行地理配准的影像中即可。

  首先，依据第九篇GEE教学博客（Google Earth Engine谷歌地球引擎GEE中JavaScript脚本语言代码基础规则与函数语句_疯狂学习GIS的博客-CSDN博客）中介绍的代码导入遥感影像方法，我们将两景 SkySat 多光谱校正公开图像集 Public Ortho Imagery 产品导入 GEE 中，并在地图中显示。 SkySat 是由美国商业公司 Planet 拥有并运行的亚米级地球遥感观测卫星星座，其多光谱波段空间分辨率可达 2 m 左右，全色波段空间分辨率更可高达 0.8 m 左右。

  随后，依据第十二篇GEE教学博客（Google Earth Engine谷歌地球引擎GEE图层投影信息与参考坐标系及其空间转换参数获取_疯狂学习GIS的博客-CSDN博客）中介绍的 Map.centerObject() 函数，将GEE地图缩放至第一景遥感影像的中心经、纬度处。

var image_1=ee.Image("SKYSAT/GEN-A/PUBLIC/ORTHO/MULTISPECTRAL/s01_20150305T081019Z");
var image_2=ee.Image("SKYSAT/GEN-A/PUBLIC/ORTHO/MULTISPECTRAL/s01_20150502T082736Z");
Map.addLayer(image_1);
Map.addLayer(image_2);
Map.centerObject(image_1);

  执行代码，得到结果如下。

此时如果对两景图像分别设置不同的可视化参数，仔细观察可以看到两景图像之间的地物位置有一定细微的差异。这里需要注意，是同一地物在两景遥感影像中的位置具有一定差异，而不是两景图像的形状与成像范围不一致（例如两景图像一个形似向北开的坦克，另一个形似向南开的坦克）；这里二者形状与成像范围不一致是由于在不同时刻成像时，轨道方向与参数不一致导致的，是正常情况。

接下来，我们需要对进行地理配准过程中遥感影像的重采样方法进行规定。在第十三篇GEE教学博客（Google Earth Engine谷歌地球引擎GEE地理坐标系与投影坐标系变换与重投影_疯狂学习GIS的博客-CSDN博客_谷歌坐标系）中我们提到，进行投影变换，也包括本文介绍的地理配准等遥感影像的空间操作时，GEE默认的重采样方法为最邻近插值方法；如果我们不想用这一默认方法，可以通过 .resample() 函数或 .reduceResolution() 函数加以实现。这里我们就用 .resample() 函数来实现。

var image_1_resample=image_1.resample("bicubic");
var image_2_resample=image_2.resample("bicubic");

其中， "bicubic" 表示通过双三次插值方法进行重采样；此外，这一参数还可以换为 "bilinear" ，表示双线性插值方法。

接下来，依据第六篇GEE教学博客（Google Earth Engine谷歌地球引擎GEE栅格代数与NDVI波段计算手动求取_疯狂学习GIS的博客-CSDN博客_gee波段计算）中介绍的 .select() 函数，选出两景图像的红光波段进行后续操作。这里需要注意，对于两景需要进行地理配准的遥感影像，二者只需要保持相互之间的波段数量一致，且两景图像的每一组波段（如图像1的红光波段与图像2的红光波段）在视觉上尽可能比较近似即可，并不一定分别只能含有一个波段；这里选出一个红光波段，仅仅是为了后续进行确立两景遥感影像之间的偏移量时方便。

var image_1_red=image_1_resample.select("R");
var image_2_red=image_2_resample.select("R");

接下来，便可开始两景遥感影像之间的地理配准；我们以将第二景遥感影像配准到第一景遥感影像为例进行介绍。

首先进行本文开头提及的第一个步骤，即设置一景“位移影像”。

var my_displacement=image_2_red.displacement({
  referenceImage:image_1_red,
  maxOffset:50.0,
  patchWidth:100.0
});

其中，需要用到 .displacement() 函数进行“位移影像”的设立； referenceImage 表示基准遥感影像； maxOffset 表示两景遥感影像之间的最大偏移量，这一参数的具体含义稍后会有讲解； patchWidth 表示用以确定两景遥感影像之间偏移量的集合大小，对于 patchWidth 这一参数我的个人理解就是，其相当于就是用以确定两景遥感影像之间偏移量的点对个数——即GEE自动在两景遥感影像范围内选定 patchWidth 个点对，随后计算偏移量；和ENVI、ArcMap等软件中地理配准操作时选定的校准点类似，只不过GEE中是系统自动选定，ENVI、ArcMap等软件中需要我们手动选定。得到的变量 my_displacement 即为“位移影像”。

用两景遥感影像的红光波段确定了两景遥感影像之间的偏移量后，即可开始对两景遥感影像的全部波段进行配准了，也就是本文开头提及的第二个步骤。

var my_registered=image_2_resample.displace(my_displacement);

其中，用 .displace() 函数实现两景遥感影像之间的配准操作，其参数 my_displacement 便是我们刚刚进行“位移影像”设立时建立的图像变量（也就是“位移影像”自身）。

完成了地理配准，我们可以将参考图像 image_1_resample 、配准前图像 image_2_resample 与配准后图像 my_registered 分别加以可视化显示，从而对比地理配准操作的结果。

var vis_parameter={
  min:1765,
  max:3725,
  bands:["R","G","B"]
};
Map.addLayer(image_1_resample,vis_parameter,"Reference");
Map.addLayer(image_2_resample,vis_parameter,"Before Registration");
Map.addLayer(my_registered,vis_parameter,"After Registration");

执行代码，得到结果如下。

为了更好看清楚地理配准的结果，我们可以将图像缩放，找到一个地物比较清晰的位置，首先对比参考图像 image_1_resample 与配准前图像 image_2_resample 之间的地物空间位置差异，随后对比参考图像 image_1_resample 与配准后图像 my_registered 之间的地物空间位置差异；通过上述对比，会非常清晰地看到地理配准后的图像与参考图像之间的地物空间位置非常一致，即地理配准的效果还是非常明显的。

我们还可以看一下“位移影像” my_displacement 长什么样子。

Map.addLayer(my_displacement);

可以通过Inspector查看 my_displacement 图层的属性。可以看到其三个波段正如本文开头所述，前两个波段为 dx 与 dy ——二者分别为图像中每一个像元在 x 和 y 方向的偏移量；第三个波段为上述偏移量的置信度。

再将第一个波段 dx 显示为灰度图像，可以看到，其呈现出一个渐变、阶梯状的色彩分布情况；这是由于其自身存储的是 x 方向的偏移量，因此其相邻两个像元之间的数值肯定是比较接近的，从而呈现出这种色彩分布情况。

此外，还可以看到前两个波段 dx 与 dy 的像元数值取值范围均为 [-50,50] ，这就是 .displacement() 函数中 maxOffset 参数的效果。

以上代码是将本文开头提及的GEE中实现遥感影像地理配准的两个步骤分开来实现的过程，我们还可以将这两个步骤合并为一个，提高地理配准的效率。

var my_registered_2=image_2_resample.register({
  referenceImage:image_1_resample,
  maxOffset:50.0,
  patchWidth:100.0
});
Map.addLayer(image_1_resample,vis_parameter,"Reference");
Map.addLayer(image_2_resample,vis_parameter,"Before Registration");
Map.addLayer(my_registered,vis_parameter,"After Registration");
Map.addLayer(my_registered_2,vis_parameter,"After Registration 2");

其中，和第一种地理配准方法相比，第二种方法中我们仅需要将 .displacement() 函数换为 .register() 函数，其它参数不变，即可实现一步到位的地理配准，省去了执行 .displace() 函数这一步；不过这种方法相信大家也可以看到，就是不能再将“位移影像”给显示出来了。如果需要显示或保存“位移影像”，还是只能用第一种地理配准方法。

执行代码，得到结果如下。

可以看到，经过两种地理配准方法操作后所得到的配准后图像 my_registered 与 my_registered_2 ，较之配准前图像 image_2_resample ，均更接近参考图像 image_1_resample ；而两种地理配准方法操作后所得到的配准后图像 my_registered 与 my_registered_2 之间也并不是完全一致的。个人认为 my_registered 与 my_registered_2 二者不完全一致是因为：GEE在选取 patchWidth 个点对时也是随机选取的，因此其得到的偏移量肯定也是不完全一致的，从而导致所得地理配准结果图层也是不完全一致的。