Google earth engine 基于面向对象遥感影像分割 SNIC分割算法

摘要被导师嫌弃N+1次，又摘要阴影了，懒得写摘要，这篇主要是在GEE上基于面向对象以Sentinel-2数据做的无监督分类算法示例，嗯呢，就是这样，您接着往下看有没有你需要的。

文章目录

• 一、Google earth engine简介
• 二、面向对象遥感影像分析方法
• 三、SNIC分割算法
• 四、总结

---

Google earth engine

    GEE是一个专门处理卫星影像和其他地球观测数据的云端运算平台（https://earthengine.google.com/），由谷歌、卡内基梅隆大学、美国地质调查局（USGS）共同开发，能够支持PB级的数据运算与地理信息数据可视化。此平台存储了大量的遥感数据，据统计，近40年来的公共数据集超过了200多个，包括MODIS、Landsat和Sentinel等系列影像产品和气象、地理、土地利用以及人口等数据集。该平台为每个用户提供免费的250G容量谷歌资源空间（Google Earth Engine Assets）和15G容量谷歌硬盘空间（Google Drive）。用户可以根据研究需要，升级个人使用空间。基于GEE平台，用户可以将本地数据上传到Google Earth Engine Assets或者Google Drive中。之后，借助GEE在线编辑器（https://code.earthengine.google.com/）进行线上数据的处理和空间分析，如影像的导入、导出、镶嵌、裁剪、信息提取与变化分析等。GEE在线编辑器支持JavaScript语言开发，通过JavaScript应用程序编程接口（API）直接实现线上数据的处理与空间分析，不需要进行环境的搭建。整个工作台可以分为脚本文件。

       与ArcG IS、ENVI等遥感影像处理软件相比，GEE平台最大的优势是不受空间和时间的约束，支持大数据的并行运算。长时间序列的遥感影像获取和分析不仅要求计算机具有高容量的存储能力，还需具备良好的数据处理性能，以实现影像筛选、辐射校正、几何校正、消除云影响、特征计算、分类处理以及优化等大批量工作。一般本地计算机难以满足这样高容量存储和运算的需要。GEE平台则能支持PB级的数据处理与分析。GEE平台另一个明显的优势是，用户能够根据研究的需要，自由编码创造新的算法或者重组算法处理遥感数据，不受固有算法的约束。 GEE平台实现了数据获取、处理、分析与应用于一体，支持并行运算，极大提高遥感图像处理效率。目前，此平台已广泛应用于旱情、森林、耕地以及水资源提取与变化监测等方面。例如，Liu等[54]基于GEE平台，提取了1990-2010年全球尺度的不透水面；陈军等[55]提取了2000和2010年全球的耕地。而本文研究基于GEE平台，使用JavaScript语言，在GEE在线编辑器中实现2009-2018年云南省的耕地覆盖信息的提取与变化监测。

如果大家想要系统的学习GEE，推荐知乎“无形的风”以及B站吴秋生老师的讲解视频，很详细哦。

无形的风：GEE开发 - 知乎 (zhihu.com)

吴秋生老师：吴秋生老师的个人空间_哔哩哔哩_bilibili

---

提示：该图来源于知乎“无形的风”博主，强烈安利！！！

基于面向对象的遥感影像分析

基于像素的遥感影像扥西方法是以像元为分类单元，通过分析地物在遥感影像上呈现的光谱、时间和空间特性，识别像元的过程。

而面向对象的遥感影像分析方法以一组像元为处理基元，综合考虑了与空间地理实体对应的光谱波段统计特征、拓扑与相邻关系等一系列因素（Diaz-Varela et al., 2014; Zhao et al.,2017），可以有效消除基于像素分类中的“椒盐现象” (Reyes，Solla和Lorenzo，2017)，有利于在高分辨率遥感中对地物边缘信息的提取（Wenyi Sun et al.,2019）。

Definiens公司出品的面向对象软件eCognition v.9.0中的多尺度分割算法是当前最受欢迎的遥感影像分割算法,（之后会写在博文里）但是该分割算法局限于小区域遥感影像分割，不适用于大区域的遥感影像分割，而云计算平台GEE的出现为大区域遥感影像分割的实现提供了可能，下面介绍景点简单非迭代聚类（SNIC）超像素分割算法。

简单非迭代聚类（SNIC）超像素分割

SNIC算法从规则网格上创建的种子开始，然后根据像素与簇质心的距离（考虑CIELAB颜色和空间坐标的归一化五维空间）将像素分组为超像素簇。随后，基于所有图像像素与现有超像素簇的连接性，将所有图像像素添加到优先级队列。每次添加到聚类质心的距离最小的像素时，为每个增长的聚类生成更新的质心值和优先级队列。

SNIC分割中需要设置的参数是种子距离、分割紧密度、连通性和邻域大小。尽管不可取，但许多将分割应用于LULC标测的研究依赖于试错法和目视检查来找到最佳分割参数。基于这种方法，SNIC参数在各种实验后选择，基于视觉检查，以在实验测试部位内产生最均匀和最有意义的对象，并符合“良好分割是显示很少过分割和没有欠分割的分割”的概念。

现在，我们将在几个编号的部分中构建一个脚本，让您有机会了解它是如何构建的，并在进行过程中观察中间结果和对比结果。我们将首先定义一个函数，用于获取图像并将其分解为一组无监督类。当调用该函数时，该函数将图像划分为指定数量的类，而不直接使用图像的任何空间特征。

将下面的代码复制到新的脚本中

// 1.1 Unsupervised k-Means classification

// This function does unsupervised clustering classification 
// input = any image. All bands will be used for clustering.
// numberOfUnsupervisedClusters = tunable parameter for how 
//        many clusters to create.
var afn_Kmeans = function(input, numberOfUnsupervisedClusters,
    defaultStudyArea, nativeScaleOfImage) {

    // Make a new sample set on the input. Here the sample set is 
    // randomly selected spatially. 
    var training = input.sample({
        region: defaultStudyArea,
        scale: nativeScaleOfImage,
        numPixels: 1000
    });

    var cluster = ee.Clusterer.wekaKMeans(
            numberOfUnsupervisedClusters)
        .train(training);

    // Now apply that clusterer to the raw image that was also passed in. 
    var toexport = input.cluster(cluster);

    // The first item is the unsupervised classification. Name the band.
    var clusterUnsup = toexport.select(0).rename(
        'unsupervisedClass');
    return (clusterUnsup);
};

我们还需要一个函数来将频带值归一化为从0到1的公共刻度。当我们创建对象时，这将是最有用的。此外，我们还需要一个函数将平均值添加到波段名称中。将以下函数粘贴到代码中。注意，代码编号故意从1.2跳到1.4;稍后我们将添加第1.3节。

// 1.2 Simple normalization by maxes function.
var afn_normalize_by_maxes = function(img, bandMaxes) {
    return img.divide(bandMaxes);
};

// 1.4 Simple add mean to Band Name function
var afn_addMeanToBandName = (function(i) {
    return i + '_mean';
});

我们将创建一个小节，定义您可以调整的变量。一个重要的可调参数是集群器使用的集群数量。将以下代码添加到函数定义下面。

// 2. Parameters to function calls

 
// 2.1. Unsupervised KMeans Classification Parameters
var numberOfUnsupervisedClusters = 4;

该脚本将允许您放大到指定的区域以更好地查看，并且已经存在于代码存储库检查点中。添加下面的代码。

// 2.2. Visualization and Saving parameters
// For different images, you might want to change the min and max 
// values to stretch. Useful for images 2 and 3, the normalized images.
var centerObjectYN = true;

现在，有了这些函数、参数和标志之后，让我们定义一个新图像，并设置有助于分析图像的特定图像值。我们将把它放在代码的一个新部分中，其中包含来自多个传感器的图像的“if”语句。我们像这样设置代码，因为我们将在接下来的部分中使用来自不同传感器的几张图像，因此它们是预加载的，所以你所要做的就是改变参数“whohimage”。在这张特别的Sentinel-2图像中，专注于区分美国华盛顿普吉特海湾的森林和非森林地区。脚本将自动缩放到感兴趣的区域

//
// 3. Statements
//

// 3.1  Selecting Image to Classify 
var whichImage = 1; // will be used to select among images
if (whichImage == 1) {
    // Image 1. 
    // Puget Sound, WA: Forest Harvest
    // (April 21, 2016)
    // Harvested Parcels
    // Clear Parcel Boundaries
    // Sentinel 2, 10m
    var whichCollection = 'COPERNICUS/S2';
    var ImageToUseID = '20160421T191704_20160421T212107_T10TDT';
    var originalImage = ee.Image(whichCollection + '/' +
    ImageToUseID);
    print(ImageToUseID, originalImage);
    var nativeScaleOfImage = 10;
    var threeBandsToDraw = ['B4', 'B3', 'B2'];
    var bandsToUse = ['B4', 'B3', 'B2'];
    var bandMaxes = [1e4, 1e4, 1e4];
    var drawMin = 0;
    var drawMax = 0.3;
    var defaultStudyArea = ee.Geometry.Polygon(
        [
            [
                [-123.13105468749993, 47.612974066532004],
                [-123.13105468749993, 47.56214700543596],
                [-123.00179367065422, 47.56214700543596],
                [-123.00179367065422, 47.612974066532004]
            ]
        ]);
    var zoomArea = ee.Geometry.Polygon(
        [
            [
                [-123.13105468749993, 47.612974066532004],
                [-123.13105468749993, 47.56214700543596],
                [-123.00179367065422, 47.56214700543596],
                [-123.00179367065422, 47.612974066532004]
            ]
        ], null, false);
}
Map.addLayer(originalImage.select(threeBandsToDraw), {
    min: 0,
    max: 2000
}, '3.1 ' + ImageToUseID, true, 1);

 现在，让我们将图像剪辑到我们感兴趣的研究区域，然后提取用于分类过程的波段。

// 4. Image Preprocessing 

var clippedImageSelectedBands = originalImage.clip(defaultStudyArea)
    .select(bandsToUse);
var ImageToUse = afn_normalize_by_maxes(clippedImageSelectedBands,
    bandMaxes);

Map.addLayer(ImageToUse.select(threeBandsToDraw), {
        min: 0.028,
        max: 0.12
    },
    '4.3 Pre-normalized image', true, 0);

 现在，让我们查看使用k-means分类器生成的逐像素无监督分类。请注意，正如前面所做的，我们跳过了代码编号的一部分(从第4节移动到第6节)，我们将在稍后进一步开发脚本时填充该部分。

// 6. Execute Classifications


// 6.1 Per Pixel Unsupervised Classification for Comparison
var PerPixelUnsupervised = afn_Kmeans(ImageToUse,
    numberOfUnsupervisedClusters, defaultStudyArea,
    nativeScaleOfImage);
Map.addLayer(PerPixelUnsupervised.select('unsupervisedClass')
    .randomVisualizer(), {}, '6.1 Per-Pixel Unsupervised', true, 0
);
print('6.1b Per-Pixel Unsupervised Results:', PerPixelUnsupervised);

// 7. Zoom if requested

if (centerObjectYN === true) {
    Map.centerObject(zoomArea, 14);
}

然后插入此代码，以便在需要时可以缩放。 

 运行脚本。它将以真彩色视图绘制研究区域，在那里您可以检查景观的复杂性，就像2016年拍摄图像时您所看到的那样。

 S2A原始影像

 SNIC分割结果

在绘制真彩色图像时，Earth Engine还执行k-means分类，并将其添加到图层集中。打开6.1每像素无监督层的可见性，它显示了使用随机选择的颜色的每像素四类分类结果。结果应该如图F3.3.2所示。

 

 

 您在基于像素的分类中注意到的噪声现在将使用基于对象的图像分析的两步方法进行改进。首先，您将使用SNIC算法分割图像，然后使用k-means无监督分类器对其进行分类。返回到脚本，我们将在其中添加一个新函数。请注意，代码的部分是编号的，找到代码部分1.2并将下面的函数添加到它下面。

// 1.3 Seed Creation and SNIC segmentation Function
var afn_SNIC = function(imageOriginal, SuperPixelSize, Compactness,
    Connectivity, NeighborhoodSize, SeedShape) {
    var theSeeds = ee.Algorithms.Image.Segmentation.seedGrid(
        SuperPixelSize, SeedShape);
    var snic2 = ee.Algorithms.Image.Segmentation.SNIC({
        image: imageOriginal,
        size: SuperPixelSize,
        compactness: Compactness,
        connectivity: Connectivity,
        neighborhoodSize: NeighborhoodSize,
        seeds: theSeeds
    });
    var theStack = snic2.addBands(theSeeds);
    return (theStack);
};

 如您所见，该函数组装了运行SNIC所需的参数(Achanta和Süsstrunk 2017, Crowley等人，2019)，该函数在图像中描绘对象。对SNIC的调用需要几个参数，我们将探讨这些参数。在代码节2.2下面添加以下代码。

// 2.3 Object-growing parameters to change
// Adjustable Superpixel Seed and SNIC segmentation Parameters:
// The superpixel seed location spacing, in pixels.
var SNIC_SuperPixelSize = 16;
// Larger values cause clusters to be more compact (square/hexagonal). 
// Setting this to 0 disables spatial distance weighting.
var SNIC_Compactness = 0;
// Connectivity. Either 4 or 8. 
var SNIC_Connectivity = 4;
// Either 'square' or 'hex'.
var SNIC_SeedShape = 'square';

// 2.4 Parameters that can stay unchanged
// Tile neighborhood size (to avoid tile boundary artifacts). Defaults to 2 * size.
var SNIC_NeighborhoodSize = 2 * SNIC_SuperPixelSize;

 现在添加一个对SNIC函数的调用。您将注意到，它接受代码第2节中指定的参数，并将它们发送给SNIC算法。将下面的代码放入脚本中，作为代码的第5节，位于第4节和第6节之间。

// 5. SNIC Clustering


// This function returns a multi-banded image that has had SNIC
// applied to it. It automatically determine the new names 
// of the bands that will be returned from the segmentation.
print('5.1 Execute SNIC');
var SNIC_MultiBandedResults = afn_SNIC(
    ImageToUse,
    SNIC_SuperPixelSize,
    SNIC_Compactness,
    SNIC_Connectivity,
    SNIC_NeighborhoodSize,
    SNIC_SeedShape
);

var SNIC_MultiBandedResults = SNIC_MultiBandedResults
    .reproject('EPSG:3857', null, nativeScaleOfImage);
print('5.2 SNIC Multi-Banded Results', SNIC_MultiBandedResults);

Map.addLayer(SNIC_MultiBandedResults.select('clusters')
    .randomVisualizer(), {}, '5.3 SNIC Segment Clusters', true, 1);

var theSeeds = SNIC_MultiBandedResults.select('seeds');
Map.addLayer(theSeeds, {
    palette: 'red'
}, '5.4 Seed points of clusters', true, 1);

var bandMeansToDraw = threeBandsToDraw.map(afn_addMeanToBandName);
print('5.5 band means to draw', bandMeansToDraw);
var clusterMeans = SNIC_MultiBandedResults.select(bandMeansToDraw);
print('5.6 Cluster Means by Band', clusterMeans);
Map.addLayer(clusterMeans, {
    min: drawMin,
    max: drawMax
}, '5.7 Image repainted by segments', true, 0);

 现在运行脚本。它将绘制几个图层，图F3.3.4所示的图层在最上面。

 

 这显示了SNIC在发送给它的图像上的工作——在这种情况下，在波段8、11和12的合成上。如果你仔细观察彩色图层，你可以看到红色的小“种子”像素。为了启动这个过程，需要创建这些种子，并按照传递给函数的参数所给出的间隔形成正方形或六角形的“超级像素”。然后，这些块的边缘被推拉，并在输入图像的边缘处停止。作为算法的一部分，一些超级像素被合并成更大的块，这就是为什么你会发现一些形状包含两个或更多的种子像素。通过改变第5层的透明度来探索结构，以判断对于给定的参数值集图像分割的表现。你也可以比较5.7层和3.1层。5.7层是3.1层的重新解释，其中5.3层中给定形状中的每个像素都被分配了形状内像素的平均值。当以一种对给定项目目标有用的方式进行参数化时，图像中同质的部分将获得相同的颜色，而高异质性的区域将获得多种颜色。

现在你可以花些时间研究控制代码参数的效果来回答下面的问题。

问题1：改变参数SNIC_SuperPixelSize对SNIC集群有什么影响?

问题2：改变参数SNIC_Compactness会有什么影响?

问题3：改变参数SNIC_Connectivity和SNIC_SeedShape会有什么影响?

本教程中使用的k-means分类器没有意识到我们通常更喜欢将相邻像素分组到同一类中——它没有物理空间的感觉。这就是为什么你会在无监督分类中看到噪音。但是，因为我们重新着色了SNIC集群中的像素以共享完全相同的频带值，k-means将每个集群中的所有像素分组为具有相同的类。在最好的情况下，这允许我们从基于像素增强我们的分类，以显示景观中干净和明确的对象。在本节中，我们将对这些对象进行分类，探索在此图像中寻找对象的优势和局限性。返回SNIC设置的第一个值，即:

// The superpixel seed location spacing, in pixels.
var SNIC_SuperPixelSize = 16;
// Larger values cause clusters to be more compact (square/hexagonal). 
// Setting this to 0 disables spatial distance weighting.
var SNIC_Compactness = 0;
// Connectivity. Either 4 or 8. 
var SNIC_Connectivity = 4;
// Either 'square' or 'hex'.
var SNIC_SeedShape = 'square';

作为代码6.2节，添加这段代码，它将调用SNIC函数并绘制结果。 

// 6.2 SNIC Unsupervised Classification for Comparison
var bandMeansNames = bandsToUse.map(afn_addMeanToBandName);
print('6.2 band mean names returned by segmentation', bandMeansNames);
var meanSegments = SNIC_MultiBandedResults.select(bandMeansNames);
var SegmentUnsupervised = afn_Kmeans(meanSegments,
    numberOfUnsupervisedClusters, defaultStudyArea,
    nativeScaleOfImage);
Map.addLayer(SegmentUnsupervised.randomVisualizer(), {},
    '6.3 SNIC Clusters Unsupervised', true, 0);
print('6.3b Per-Segment Unsupervised Results:', SegmentUnsupervised);

 当您运行脚本时，新函数将在5.7层中对图像进行分类，这是根据第5层中显示的片段对原始图像进行重新着色。比较超级像素的分类(6.3)和逐像素值的无监督分类(6.1)。您应该能够更改这两个层的透明度，以便直接比较它们。

总结

SNIC分割算法有效解决了大区域计算数量的难题，关于其中最优参数的调节要根据具体数据的使用情况来看，如果研究区域太大，可以将研究区分为几个小瓦片进行同样也可以完成分割，时隔半年才来写的总结，前面的代码如果哪里有问题可以在评论区留言。