基于 GEE 计算核植被指数 kNDVI

基于 GEE 计算核植被指数 kNDVI

引用

CSDN

1.

https://blog.csdn.net/ZK180531/article/details/145551198

kNDVI（核植被指数）是近年来提出的一种新型植被指数，它通过引入机器学习的核方法理论，解决了传统NDVI（归一化差异植被指数）和EVI（增强型植被指数）在高植被覆盖度区域的饱和问题。本文将详细介绍kNDVI的计算原理，并提供完整的GEE（Google Earth Engine）代码实现，帮助读者更好地理解和应用这一新的植被指数。

1 核植被指数（kNDVI）

目前，NDVI依然是比较流行的植被指数。但由于NDVI自身是“非线性”的，当植被覆盖持续增加时，比值NIR/RED持续增加，这种非线性的转换过程也使得NDVI对RED波段过度敏感，因此，在植被覆盖较高的或者生物含量较高地区，NDVI不可避免的伴有不同程度的饱和现象。此外，NDVI对红光(叶绿素)高度敏感，当叶绿索浓度超过一定浓度时，红光对叶绿素不再具有敏感性，很快出现饱和现象。

为了调整NDVI以适应大气和土壤噪声，特别是在植被茂密地区，减缓NDVI饱和性问题，在NDVI中加入了蓝光波段，提出了增强型植被指数EVI，但是EVI中饱和现象依旧存在。

有学者为解决NDVI（归一化差异植被指数）和EVI（增强型植被指数）的饱和度问题上存在的缺陷，通过采用机器学习的核方法理论，将NDVI线性化，转变为了kNDVI（核植被指数）。

kNDVI指数计算公式如下：

n为NIR；r为Red。核函数k表示这两个波段之间的相似性。在任何情况下均通过径向基函数（RBF）再生核，sigma控制近红外和红外波段之间的距离概念。

对于径向基函数（RBF），自相似性k（n，n）为1，可先简化为下式。

最后计算公式可以简化为：

kNDVI的数值含义与NDVI类似，但提供了更细致的反映植被状况的能力，尤其是在高植被覆盖度的地区。kNDVI的数值范围通常在-1到1之间，但实际应用中，它主要取值在0到1的区间。下面是kNDVI数值的一些含义：

（1）接近0：kNDVI值接近0通常表示植被覆盖度较低，可能是由于植被稀疏、植被处于非生长季节、或者是非植被区域（如裸土、岩石、水体等）。

（2）接近1：kNDVI值接近1表示植被覆盖度很高，通常与健康、茂盛的植被相关。这可能意味着该地区有丰富的植被生物量和活跃的光合作用。

（3）中间值：kNDVI值在0到1之间的中间值表示不同程度的植被覆盖度，具体解释可能需要结合具体的生态系统和环境条件。

（4）负值：虽然理论上kNDVI可以取负值，但在实际应用中很少见，因为kNDVI的设计是为了减少负值的出现，并且提高在低植被覆盖度区域的敏感性。

kNDVI的一个关键优势是它通过核方法减少了NDVI在高植被覆盖度区域的饱和现象，使得它能够更准确地反映植被的生物物理和生理参数，如叶面积指数（LAI）、总初级生产力（GPP）和太阳诱导叶绿素荧光（SIF）等。这使得kNDVI在植被监测、物候和绿化研究、参数升级等领域具有重要的应用潜力。在实际应用中，kNDVI的数值需要结合具体的遥感影像、时间序列分析、地面验证数据以及其他环境参数来综合解读，以获得更准确的植被状况评估。

2 完整代码

// Define region
var empty = ee.Image().toByte();
var outline = empty.paint({
    featureCollection: geometry,
    color: 0,
    width: 3
});
// Draw a red border
Map.addLayer(outline, {palette: "ff0000"}, "geometry");
// Function to extract QA bits
var getQABits = function(image, start, end, newName) {
    var pattern = 0;
    for (var i = start; i <= end; i++) {
        pattern += Math.pow(2, i);
    }
    return image.select([0], [newName])
                .bitwiseAnd(pattern)
                .rightShift(start);
};
// Function to mask out cloudy pixels
var maskClouds = function(image) {
    var QA = image.select('state_1km');
    var internalCloud = getQABits(QA, 0, 1, 'cloud_state').eq(1).or(getQABits(QA, 0, 1, 'cloud_state').eq(2));
    image = image.addBands(internalCloud);
    return image.updateMask(internalCloud.eq(0));
};
// Function to mask out non-land pixels
var masknoLand = function(image) {
    var QA = image.select('state_1km');
    var internalLand = getQABits(QA, 3, 5, 'Land/water_flag').gt(1).or(getQABits(QA, 3, 5, 'Land/water_flag').eq(0));
    image = image.addBands(internalLand);
    return image.updateMask(internalLand.eq(0));
};
// Function to mask out pixels with no observations
var maskEmptyPixels = function(image) {
    var withObs = image.select('num_observations_1km').gt(0);
    return image.updateMask(withObs);
};
var collection = ee.ImageCollection('MODIS/006/MOD09GA')
                  .filterDate('2018-01-01', '2019-01-11');
// Clip the geometry using a map function
var clipImage = function(image) {
    return image.clip(geometry);
};
collection = collection.map(clipImage)
                      .map(maskEmptyPixels)
                      .map(masknoLand);
var collectionCloudMasked = collection.map(maskClouds);
// Select relevant bands and scale reflectance values
collectionCloudMasked = collectionCloudMasked.select([
    'sur_refl_b01', 'sur_refl_b02', 'sur_refl_b03', 'sur_refl_b04'
]).map(function(image) {
    return image.multiply(ee.Image.constant(1e-4))
                .set({'system:time_start': image.get('system:time_start')});
});
var vizParams = {
    bands: ["sur_refl_b01", "sur_refl_b04", "sur_refl_b03"],
    min: 0,
    max: 0.15
};
Map.centerObject(geometry, 9);
// Function to compute kNDVI
var addKNDVI = function(image) {
    var red = image.select('sur_refl_b01');
    var nir = image.select('sur_refl_b02');
    var D2 = nir.subtract(red).pow(2).rename('d2');
    var sigma = 0.15;
    var kndvi = D2.divide(ee.Number(sigma)).divide(ee.Number(sigma)).divide(4.0).tanh().rename('kndvi');
    return image.addBands(kndvi);
};
// Function to compute NDVI
var addNDVI = function(image) {
    var red = image.select('sur_refl_b01');
    var nir = image.select('sur_refl_b02');
    var ndvi = nir.subtract(red).divide(nir.add(red)).rename('ndvi');
    return image.addBands(ndvi);
};
var collectionNDVI = collectionCloudMasked.map(addNDVI);
var collectionKNDVI = collectionCloudMasked.map(addKNDVI).max();
print(collectionKNDVI);
var visParam = {
    palette: 'FFFFFF, CE7E45, DF923D, F1B555, FCD163, 99B718, 74A901, 66A000, 529400,' +
             '3E8601, 207401, 056201, 004C00, 023B01, 012E01, 011D01, 011301'
};
Map.addLayer(collectionNDVI.select("ndvi"), visParam, 'NDVI', false);
Map.addLayer(collectionKNDVI.select("kndvi"), visParam, 'kNDVI', false);

3 运行结果

研究区归一化植被指数NDVI

研究区核植被指数kNDVI