资讯详情

【GEE笔记】最大类间方差法（otsu、大津法）算法实现——计算阈值、图像二值化分割

主要内容

1.概述最大类间方差法的原理

2、GEE频率分布统计，直方图绘制

3.具体实现算法GEE JavaScript版本为例

4.以遥感图像提取水为例，提取目标像元

算法原理

概念

最大类间方差法(也称otsu、日本学者大津法OTSU一种高效算法是在1979年提出的图像二值化。算法假设图像根据频率分布直方图包含(前景像元和背景像元)，计算可类可以分开，使其类内方差最小；由于两两平方距恒定，即其类间方差最大。

计算简单快捷，不受图像亮度和对比度的影响。

对图像噪声敏感；如果图像中的前景像元和背景像元的面积差异很大，直方图没有明显的双峰，或者两峰的大小差异很大，分割效果不好。在实际应用中，通常与其他方法相结合。

原理

借助穷举法搜索，大津法可以最大限度地降低类内方差的阈值，计算两类方差的加权和权值。前人证明最小化类内方差和最大化类间方差是一样的。

图像的总平均灰度为：M=P1(t) * M0 P2(t) * M1

前景与背景像元类之间的方差：S=P1(t) * (M1 - M) * (M1 - M) P2(t) * (M2 - M) * (M2 - M)

t前景像元占图像的比例是前景和背景的划分阈值P1(t)，平均灰度为M1.背景像元占图像的比例P2(t)，平均灰度为M1。

算法实现

数据准备

1.原始图像：定义示例矢量区域geometry(山东省潍坊市峡山水库周边，便于提取水体)，时间范围，云量低于20%，筛选出合格的Landsat8影像集dataset，获得中值合成l8_image

var geometry =      ee.Geometry.Polygon(         [[[119.3140376290338, 36.559328749628065],           [119.3140376290338, 36.263933411986294],           [119.62234146204162, 36.263933411986294],           [119.62234146204162, 36.559328749628065]]], null, false);  var dataset = ee.ImageCollection('LANDSAT/LC08/C02/T1_L2')     .filterDate('2021-10-01', '2021-12-01')     .filterBounds(geometry)     .filter(ee.Filter.lt('CLOUD_COVER', 20));   function applyScaleFactors(image) {     var opticalBands = image.select('SR_B.').multiply(0.0000275).add(-0.2);     var thermalBands = image.select('ST_B.*').multiply(0.00341802).add(149.0);     return image.addBands(opticalBands, null, true)         .addBands(thermalBands, null, true); }  //获取图像 中值合成 var l8_image = dataset     .map(applyScaleFactors)     .median().clip(geometry);  //可视化 var visualization = {     bands: ['SR_B4', 'SR_B3', 'SR_B2'],     min: 0.0,     max: 0.3, }; Map.centerObject(geometry); Map.addLayer(l8_image, visualization, 'True Color (432)');

原始图像真色显示：

2、提取NDWI：后续用于计算阈值，识别水和非水

NDWI(Normalized Difference Water Index,归一化水指数)用遥感图像的特定波段进行归一化差处理，以突出图像中的水信息。表达式为:

NDWI =(p(Green)-p(NIR))/(p(Green) p(NIR))

它是基于绿色波段外波段的归一化比值指数，一般用于提取图像中的水信息，效果较好。

var ndwi = l8_image.normalizedDifference(['SR_B3', 'SR_B5']).float().rename('l8_NDWI'); print("ndwi", ndwi) var visParams1 = {     min: 0, max: 1, palette: ['FFFFFF', 'CE7E45', 'DF923D', 'F1B555', 'FCD163', '99B718',         '74A901', '66A000', '529400', '3E8601', '207401', '056201',         '004C00', '023B01', '012E01', '011D01', '011301'] }; Map.addLayer(ndwi, visParams1, "l8_NDWI");

NDWI灰度显示：

数据查看

NDWI计算并显示，可见有明显的双峰，且两类无重叠，使用大津法计算分割阈值非常方便，其中参数： maxBuckets minBucketWidth设置更为关键，其他具体参数参考官方文件

GEE右上角“Download”可以将csv数据下载到本地进行分析

//频数分布直方图 var chart = ui.Chart.image.histogram({     image: ndwi,     region: geometry,     scale: 30,     maxBuckets: 1000,///最大组数     minBucketWidth: 0.01, ///最小组距离     // maxRaw,      maxPixels: 1e13 }) print(chart)

计算阈值

首先，计算频数分布数据，并注意上述数据，然后用穷举法计算类内方差，得到类间方差结果表，按方差排序，得到方差最大对应值为最佳阈值，

阈值=0.0650840364171252

详见代码注释

//计算OTSU阈值  var yuzhi = otsu(ndwi) print("阈值", yuzhi)  //OTSU function otsu1(histogram) {     // 各组频数     var counts = ee.Array(ee.Dictionary(histogram).get('histogram'))     // 各组的值     var means = ee.Array(ee.Dictionary(histogram).get('bucketMeans'))     // 组数     var size = means.length().get([0])     // 总像元数量     var total = counts.reduce(ee.Reducer.sum(), [0]).get([0])     // 所有组的值之和     var sum = means.multiply(counts).reduce(ee.Reducer.sum(), [0]).get([0])     // 整个图像的平均值     var mean = sum.divide(total)      // 索引的长度与组数相同     var indices = ee.List.sequence(1, size)     // 穷举计算类内方差     var bss = indices.map(function (i) {         // 当 i = 1, aCounts = [counts[0]], 当 i = 2, aCounts = [counts[0], counts[1]]          //从i分为两类A、B 计算A方差         var aCounts = counts.slice(0, 0, i)         var aCount = aCounts.reduce(ee.Reducer.sum(), [0]).get([0])         var aMeans = means.slice(0, 0, i)         // 类别A均值         var aMean = aMeans.multiply(aCounts)             .reduce(ee.Reducer.sum(), [0]).get([0])             .divide(aCount)          var bCount = total.subtract(aCount)         // 类别B均值        var bMean = sum.subtract(aCount.multiply(aMean)).divide(bCount)
        //类间方差公式
        return aCount.multiply(aMean.subtract(mean).pow(2)).add(
            bCount.multiply(bMean.subtract(mean).pow(2)))
    })

    print('类间方差', ui.Chart.array.values(ee.Array(bss), 0, means))

    // 排序选出最适阈值
    return means
        .sort(bss)
        .get([-1])
}
function otsu(image) {
    var histogram = image.reduceRegion({
        reducer: ee.Reducer.histogram(1000, 0.01),// 自行修改合适的最大组数，最小组距
        geometry: geometry,
        scale: 30,
        bestEffort: true,
        // tileScale:16
    });
    print("频数分布", histogram)
    return otsu1(histogram.get(histogram.keys().get(0)));
}

分割图像

借助前文得到的的阈值进行影像二值化分割，得到水体示意图

var result = ndwi.gt(yuzhi)
Map.addLayer(result.randomVisualizer(), "", "water");

结果图：

代码附录

链接：https://code.earthengine.google.com/4456dc4c799a8673d0f1aec1431250f4

var geometry = 
    /* color: #d63000 */
    /* shown: false */
    /* displayProperties: [
      {
        "type": "rectangle"
      }
    ] */
    ee.Geometry.Polygon(
        [[[119.3140376290338, 36.559328749628065],
          [119.3140376290338, 36.263933411986294],
          [119.62234146204162, 36.263933411986294],
          [119.62234146204162, 36.559328749628065]]], null, false);

var dataset = ee.ImageCollection('LANDSAT/LC08/C02/T1_L2')
    .filterDate('2021-10-01', '2021-12-01')
    .filterBounds(geometry)
    .filter(ee.Filter.lt('CLOUD_COVER', 20));


function applyScaleFactors(image) {
    var opticalBands = image.select('SR_B.').multiply(0.0000275).add(-0.2);
    var thermalBands = image.select('ST_B.*').multiply(0.00341802).add(149.0);
    return image.addBands(opticalBands, null, true)
        .addBands(thermalBands, null, true);
}

//获取图像 中值合成
var l8_image = dataset
    .map(applyScaleFactors)
    .median().clip(geometry);

//可视化
var visualization = {
    bands: ['SR_B4', 'SR_B3', 'SR_B2'],
    min: 0.0,
    max: 0.3,
};
Map.centerObject(geometry);
Map.addLayer(l8_image, visualization, 'True Color (432)');

var ndwi = l8_image.normalizedDifference(['SR_B3', 'SR_B5']).float().rename('l8_NDWI');
print("ndwi", ndwi)
var visParams1 = {
    min: 0, max: 1
};
Map.addLayer(ndwi, visParams1, "l8_NDWI");

//频数分布直方图
var chart = ui.Chart.image.histogram({
    image: ndwi,
    region: geometry,
    scale: 30,
    maxBuckets: 1000,//最大组数
    minBucketWidth: 0.01, //最小组距
    // maxRaw, 
    maxPixels: 1e13
})
print(chart)

//计算OTSU阈值 分割图像
var yuzhi = otsu(ndwi)
print("阈值", yuzhi)
var result = ndwi.gt(yuzhi)
Map.addLayer(result.randomVisualizer(), "", "water");
//
OTSU/

function otsu1(histogram) {
    // 各组频数
    var counts = ee.Array(ee.Dictionary(histogram).get('histogram'))
    // 各组的值
    var means = ee.Array(ee.Dictionary(histogram).get('bucketMeans'))
    // 组数
    var size = means.length().get([0])
    // 总像元数量
    var total = counts.reduce(ee.Reducer.sum(), [0]).get([0])
    // 所有组的值之和
    var sum = means.multiply(counts).reduce(ee.Reducer.sum(), [0]).get([0])
    // 整幅影像的均值
    var mean = sum.divide(total)

    // 与组数相同长度的索引
    var indices = ee.List.sequence(1, size)
    // 穷举法计算类内方差
    var bss = indices.map(function (i) {
        // 当 i = 1, aCounts = [counts[0]], 当 i = 2, aCounts = [counts[0], counts[1]] 
        //从i分割为两类A、B 计算A方差
        var aCounts = counts.slice(0, 0, i)
        var aCount = aCounts.reduce(ee.Reducer.sum(), [0]).get([0])
        var aMeans = means.slice(0, 0, i)
        // 类别A均值
        var aMean = aMeans.multiply(aCounts)
            .reduce(ee.Reducer.sum(), [0]).get([0])
            .divide(aCount)

        var bCount = total.subtract(aCount)
        // 类别B均值
        var bMean = sum.subtract(aCount.multiply(aMean)).divide(bCount)
        //类间方差公式
        return aCount.multiply(aMean.subtract(mean).pow(2)).add(
            bCount.multiply(bMean.subtract(mean).pow(2)))
    })

    print('类间方差', ui.Chart.array.values(ee.Array(bss), 0, means))

    // 排序选出最适阈值
    return means
        .sort(bss)
        .get([-1])
}
function otsu(image) {
    var histogram = image.reduceRegion({
        reducer: ee.Reducer.histogram(1000, 0.01),
            // .combine('mean', null, true)
            // .combine('variance', null, true),
        geometry: geometry,
        scale: 30,
        bestEffort: true,
        // tileScale:16
    });
    print("频数分布", histogram)
    return otsu1(histogram.get(histogram.keys().get(0)));
}

标签： d01一体化瓦振变送器

锐单商城拥有海量元器件数据手册、 IC替代型号，打造 电子元器件IC百科大全！

锐单商城 - 一站式电子元器件采购平台