Google Earth Engine谷歌地球引擎遥感影像自动化地理配准、空间校正

  本文主要对GEE中栅格图像的地理配准（空间坐标位置校正）操作加以介绍。本文是谷歌地球引擎（Google Earth Engine，GEE）系列教学文章的第十四篇。

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

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

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

  随后，依据Google Earth Engine谷歌地球引擎获取投影信息、坐标系及其转换参数中介绍的Map.centerObject()函数，将GEE地图缩放至第一景遥感影像的中心经、纬度处。

1var image_1=ee.Image("SKYSAT/GEN-A/PUBLIC/ORTHO/MULTISPECTRAL/s01_20150305T081019Z");

2var image_2=ee.Image("SKYSAT/GEN-A/PUBLIC/ORTHO/MULTISPECTRAL/s01_20150502T082736Z");

3Map.addLayer(image_1);

4Map.addLayer(image_2);

5Map.centerObject(image_1);

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

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

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

1var image_1_resample=image_1.resample("bicubic");

2var image_2_resample=image_2.resample("bicubic");

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

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

1var image_1_red=image_1_resample.select("R");

2var image_2_red=image_2_resample.select("R");

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

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

1var my_displacement=image_2_red.displacement({

2  referenceImage:image_1_red,

3  maxOffset:50.0,

4  patchWidth:100.0

5});

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

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

1var my_registered=image_2_resample.displace(my_displacement);

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

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

1var vis_parameter={

2  min:1765,

3  max:3725,

4  bands:["R","G","B"]

5};

6Map.addLayer(image_1_resample,vis_parameter,"Reference");

7Map.addLayer(image_2_resample,vis_parameter,"Before Registration");

8Map.addLayer(my_registered,vis_parameter,"After Registration");

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

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

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

1Map.addLayer(my_displacement);

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

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

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

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

1var my_registered_2=image_2_resample.register({

2  referenceImage:image_1_resample,

3  maxOffset:50.0,

4  patchWidth:100.0

5});

6Map.addLayer(image_1_resample,vis_parameter,"Reference");

7Map.addLayer(image_2_resample,vis_parameter,"Before Registration");

8Map.addLayer(my_registered,vis_parameter,"After Registration");

9Map.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个点对时也是随机选取的，因此其得到的偏移量肯定也是不完全一致的，从而导致所得地理配准结果图层也是不完全一致的。