基于 GEE 利用 MODIS 数据计算核植被指数 kNDVI
基于 GEE 利用 MODIS 数据计算核植被指数 kNDVI
核植被指数（kNDVI）是为了解决传统NDVI和EVI在高植被覆盖区域的饱和问题而提出的。通过采用机器学习的核方法理论，kNDVI能够更准确地反映植被的生物物理和生理参数。本文将详细介绍kNDVI的原理、计算公式，并通过Google Earth Engine（GEE）平台基于MODIS数据进行实际计算。
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");
var getQABits = function(image, start, end, newName) {
// Compute the bits we need to extract
var pattern = 0;
for (var i = start; i <= end; i++) {
pattern += Math.pow(2, i);
}
// Return a single band image of the extracted QA bits, giving the band a new name
return image.select([0], [newName])
.bitwiseAnd(pattern)
.rightShift(start);
};
// A function to mask out cloudy pixels
var maskClouds = function(image) {
// Select the QA band
var QA = image.select('state_1km');
// Get the cloud_state bits and find cloudy areas
var internalCloud = getQABits(QA, 0, 1, 'cloud_state');
var cloudMask = internalCloud.eq(1).or(internalCloud.eq(2));
image = image.addBands(internalCloud);
return image.updateMask(cloudMask.eq(0));
};
// A function to mask out no Land pixels
var masknoLand = function(image) {
// Select the QA band
var QA = image.select('state_1km');
// Get the land/water flag bits
var internalCloud = getQABits(QA, 3, 5, 'Land/water flag');
var landMask = internalCloud.gt(1).or(internalCloud.eq(0));
image = image.addBands(internalCloud);
return image.updateMask(landMask.eq(0));
};
// A function to mask out pixels that did not have observations
var maskEmptyPixels = function(image) {
// Find pixels that had observations
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
function clipImage(image) {
return image.clip(geometry);
}
var collection = collection.map(clipImage);
// Mask out no land areas and areas that were not observed
collection = collection.map(maskEmptyPixels).map(masknoLand);
// Map the cloud masking function over the collection
var collectionCloudMasked = collection.map(maskClouds);
// This is important to have some intuition about the sigma parameter in the RBF kernel
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);
// 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).select([0], ['d2']);
// Fix sigma parameter to a reasonable value
// Recommended values range between [0.1-0.3] for normalized reflectance data
var sigma = 0.15;
var kndvi = D2.divide(ee.Number(sigma))
.divide(ee.Number(sigma))
.divide(4.0)
.tanh()
.select([0], ['kndvi']);
return image.addBands(kndvi);
};
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))
.select([0], ['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计算结果