前言
在使用遥感模型对陆地生态系统碳循环的模拟研究中，研究者一般以植被净初级生产力（NPP）或净生态系统生产力（NEP）作为基本指示因子。NEP指生态系统在单位时间和单位面积上碳的收支状况，它不仅是反映植物活动的重要变量，而且是鉴定生态系统碳源、碳汇以及调节生态过程的主要因子。本文使用遥感模型，参考CASA模型的计算方法，基于GEE平台实现区域尺度的生态系统碳收支分布的估算。
遥感模型理论基础
遥感模型，也称为光能利用率模型（LUE模型），主要基于Monteith提出的光能利用率理论，计算植被层吸收的光合有效辐射（PAR）以及植被将其转化为植物有机碳的效率，从而估算植被的净初级生产力（NPP）。该过程的基本机制如下：
首先，绿色植物通过光合作用吸收光能，将二氧化碳和水转化为有机物。假设太阳总辐射为A，但植物主要利用波长在0.4-0.7μm范围内的辐射，这部分辐射称为植被可利用的入射光合有效辐射（PAR）。PAR通常通过将太阳总辐射乘以一定比例系数k（一般取0.5）来计算。接下来，植物吸收一定比例的入射光合有效辐射，这部分被称为植被吸收的光合有效辐射（APAR），而吸收的比例称为植被冠层对光合有效辐射的吸收比率（FPAR）。最终，植被吸收的光合有效辐射（APAR）中的一部分会转化为植被的有机物质，这一转化过程的效率称为光能利用率（LUE），所转化出的有机物质即为净初级生产力（NPP）。计算公式如下：
可以类比为我们吃饭的过程：整桌菜代表太阳辐射，符合口味的菜肴相当于入射光合有效辐射（PAR），我们吃下的饭菜则相当于植被吸收的光合有效辐射（APAR），而最终转化为我们身体营养物质的部分，类似于光能利用率（LUE）转化为净初级生产力的过程。
流程图如下:
由上述可知，光能利用率模型具有一定的生理基础、可以引入遥感数据、可以进行大尺度 NPP模拟等，但光能利用率模型只能模拟 NPP 和 GPP(总初级生产力)，无法模拟衡量碳源/汇的重要指标 NEP。
由上篇推文可知NEP可以由净初级生产力减去异养生物呼吸(Rh)消耗有机物质得到。
因为这里的异养生物呼吸主要是指微生物对植被凋落物分解，所以很多研究通过计算土壤呼吸量表征Rh。
那么也就是说如果我要计算NEP，先通过光能利用率模型计算NPP，再用土壤数据计算Rh，相差得到NEP。
估算NEP
反演思路
若反演月累计NEP，先得获取这个月的太阳总辐射，在计算FPAR，其取决于植被覆盖状况和植被类型可由ndvi计算获得，光能转化率主要受温度、降水等的影响，引入月均温度等因子带入到已研究的公式计算光能转化率。最后各提取的参数相乘得到NPP，最后由土壤数据计算Rh，相差得到NEP。具体可参考《近30年中国陆地生态系统碳收支时空变化模拟研究》。公式如下：
数据准备
- 太阳总辐射数据，太阳辐射到达地球大气层的部分被称为总辐射，也叫短波辐射。这里我们可以使用。太阳总辐射数据、降雨、空气温度数据都可在NASA/FLDAS/NOAH01/C/GL/M/V001数据集中找到
- 土地分类数据，不同的植被类型具有不同的光能利用率，使用modis的土地分类数据
- FPAR计算，FPAR取决于植被覆盖状况和植被类型，使用modis的ndvi数据计算。
- 光能利用率，计算公式太长了，具体可参考《近30年中国陆地生态系统碳收支时空变化模拟研究》。图方便的话直接使用已有的光能利用率研究成果。
- 土壤呼吸Rh计算，由月平均气温和土壤含碳量数据计算得出。
利用GEE平台实现计算
该代码我已经写到计算npp，最后导入soc数据就能计算NEP了，晚上有局，我得吃席去了，不写了，拜拜。
var startDate=ee.Date.fromYMD(2023,6,1)
var endDate=ee.Date.fromYMD(2023,7,2)
// 导入数据
var dataset=ee.ImageCollection('NASA/FLDAS/NOAH01/C/GL/M/V001')
.filter(ee.Filter.date(startDate, endDate)).first().clip(geometry)
var addIndices = function(image) {
// 计算月累计雨量
var rain = image.expression(
'rain * (60 * 60 * 24 * 30)', {
'rain': image.select('Rainf_f_tavg')
}).rename('rain');
// 将温度从开尔文转换为摄氏度
var temp = image.expression(
'temp - 273.15', {
'temp': image.select('Tair_f_tavg')
}).rename('temp');
// 转换为月累计辐射能量 (J/m²)，转换为 MJ/m²，调整为光合有效辐射
var par = image.expression(
'rad * (60 * 60 * 24 * 30) / 1000000 * 0.5', {
'rad': image.select('Swnet_tavg')
}).rename('par');
// 添加计算结果为新波段
return image
.addBands(rain)
.addBands(temp)
.addBands(par);
};
// 调用函数并选择结果波段，同时重投影
var rain_temp_par = addIndices(dataset)
.select(['rain', 'temp', 'par'])
.reproject({crs: 'EPSG:4326', scale: 1000});
var ndvi= ee.ImageCollection('MODIS/061/MOD13A2')
.filter(ee.Filter.date(startDate, endDate))
.select('NDVI').first().clip(geometry);
var constant = ee.Image(1).clip(geometry)
.reproject({crs: 'EPSG:4326',scale: 1000});
var ndvi_pro = ndvi.reproject({crs: 'EPSG:4326', scale: 1000}).multiply(0.0001)
var SR=ndvi_pro.add(constant).divide(constant.subtract(ndvi_pro)).rename("SR")
// 计算 SR 的最小值和最大值
var SR_min = SR.reduceRegion({
reducer: ee.Reducer.min(),
geometry: geometry,
scale: 1000,
maxPixels: 1e13
}).get('SR');
var SR_max = SR.reduceRegion({
reducer: ee.Reducer.max(),
geometry: geometry,
scale: 1000,
maxPixels: 1e13
}).get('SR');
// print(SR_min,SR_max)
// // 计算 FPAR
var FPAR = SR.expression(
'min((SR - SR_min) / (SR_max - SR_min), 0.95)', {
'SR': SR,
'SR_min': ee.Number(SR_min),
'SR_max': ee.Number(SR_max)
}
);//这里的 min是限制最大值不超过0.95
//----------------------
//计算光能利用率
var Evapotranspiration = ee.ImageCollection('MODIS/061/MOD16A2')
.filter(ee.Filter.date(startDate, endDate))
.first().clip(geometry)
.select(['ET', 'PET']);
// 计算水分胁迫影响系数 WE
var WE = Evapotranspiration.expression(
'0.5 * ((0.5 + E) / Ep)', {
'E': Evapotranspiration.select('ET'),
'Ep': Evapotranspiration.select('PET')
}
).rename('WE')
var epsilon_max = ee.Image(0.95).clip(geometry).reproject({crs: 'EPSG:4326', scale: 1000});
var T_opt = rain_temp_par.select('temp'); // 示例：以当前温度为最适温度
Map.addLayer(WE)
// 低温限制因子 Te1
var Te1 = ee.Image(0.8)
.add(ee.Image(0.027).multiply(T_opt))
.subtract(ee.Image(0.00005).multiply(T_opt.pow(2)))
.rename('Te1');
// 高温限制因子 Te2
var Te2 = ee.Image(1.184)
.divide(
ee.Image(1).add((ee.Image(0.2).multiply(T_opt.subtract(10).subtract(rain_temp_par.select('temp')))).exp())
).multiply(
ee.Image(1).divide(
ee.Image(1).add((ee.Image(0.3).multiply(T_opt.subtract(10).add(rain_temp_par.select('temp')))).exp())
)
).rename('Te2');
// 计算最终光能利用率 epsilon
var epsilon = Te1.multiply(Te2).multiply(WE).multiply(epsilon_max).rename('epsilon');
print(epsilon)
// 可视化光能利用率
// Map.addLayer(epsilon, {min: 0, max: 0.5, palette: ['purple', 'blue', 'green', 'yellow', 'red']}, '光能利用率');
var NPP = rain_temp_par.select("par").multiply(FPAR).multiply(epsilon)
// // 显示结果
Map.addLayer(NPP, {min: 0, max: 10}, 'FPAR');
//----------------------
//计算土壤呼吸
var SOC =ee.Image("OpenLandMap/SOL/SOL_ORGANIC-CARBON_USDA-6A1C_M/v02")
.select('b10').clip(geometry)
