1. RULSE模型介绍
RUSLE(Revised Universal Soil Loss Equation),是修正的通用土壤流失方程。修正的通用土壤流失方程(RUSLE)是一种公认的、已得到广泛应用的经验土壤侵蚀估算模型。基于5个因子即:(1)降雨侵蚀力因子(R)、(2)土壤侵蚀力因子(K)、(3)坡长和陡度因子(LS)、(4)植被覆盖管理因子(C)、(5)水土保持实践因子(P)。
式中:1)降雨侵蚀力因子(R)以动能的形式评价降雨影响的效果,并预测与该降水事件直接相关的径流量和径流量;
2)土壤可蚀性因子(K)表示在一个标准条件下,土壤抵抗雨滴影响所产生的侵蚀的能力,以及降雨影响所产生的径流量和径流量;
3)斜率长度(L)和陡度(S)因素反映了区域地形对土壤侵蚀率的影响;
4)作物经营管理因子(C)直接受植被类型、植被生长阶段和植被覆盖度的影响;
5)土壤保持措施因子(P)显示了实施措施将减少径流量和径流量的效果,从而降低了土壤侵蚀的数量和速率。
接下来将利用GEE计算R和C因子。
2. GEE计算R和C因子并导出计算LS因子所需的高程数据
此处参考推文:基于GEE的RULSE
//导入研究区
var roi = ee.FeatureCollection('projects/your_roi')
Map.centerObject(roi,12)
// 导入Landsat影像、降雨及高程数据集
//Landsat image
var Landsat8 = ee.ImageCollection("LANDSAT/LC08/C02/T1_L2");
//CHIRPS precipitation dataset -- daily (mm/day)
var CHIRPS_daily = ee.ImageCollection("UCSB-CHG/CHIRPS/DAILY")
.select('precipitation') //select precipitation
.map(function(image){return image.clip(roi.geometry())});
//ALOS DEM
var ALOS = ee.ImageCollection("JAXA/ALOS/AW3D30/V3_2").select('DSM')
.map(function(image){return image.clip(roi.geometry())});
//Map.addLayer(ALOS.median().randomVisualizer())
//Landsat 8去云及波段定标
var cloudMaskL8 = function (image) {
var opticalBands = image.select('SR_B.').multiply(0.0000275).add(-0.2)
// Bits 3 and 5 are cloud shadow and cloud, respectively.
var cloudShadowBitMask = (1 << 4);
var cloudsBitMask = (1 << 3);
var dilatedcloud = (1 << 1);
// Get the pixel QA band.
var qa = image.select('QA_PIXEL');
// Both flags should be set to zero, indicating clear conditions.
var mask = qa.bitwiseAnd(cloudShadowBitMask).eq(0)
.and(qa.bitwiseAnd(cloudsBitMask).eq(0))
.and(qa.bitwiseAnd(dilatedcloud).eq(0));
return image.addBands(opticalBands,null,true)
.updateMask(mask).clip(roi.geometry());
}
//Calculate NDVI
function findndvi_l8(image){
var vi = image.normalizedDifference(['SR_B5','SR_B4']).rename('ndvi')
return vi.updateMask(vi.gt(-1).and(vi.lt(1)))
}
//选择时间段
var img2020 = Landsat8.filterDate('2019-12-31','2021-12-31')
.filterBounds(roi)
.map(cloudMaskL8).median();
Map.addLayer(img2020.select(['SR_B5','SR_B4','SR_B3']),{min:0,max:0.3})
var ndvi_img2020 = findndvi_l8(img2020).select(['ndvi']);
//Map.addLayer(ndvi_img2020.randomVisualizer())
//因子C的计算
// C ------- Vegetation Cover Factor, a=2, b=1,
// The higher the C, the worst the growth
var cf_2020 = ((ndvi_img2020.multiply(-2)).divide(ndvi_img2020.multiply(-1).add(1))).exp();
//R因子的计算
// Calculate R factor based on Flampouris'
//phd thesis "Study of the effect of rainfall factor R on Rusle's law (in Greek)"
//(2008) (No. GRI-2008-1712),Aristotle University of Thessaloniki.
//He used gauge stations and created some stable values
//Retrieve CHIRPS percipitation data
var jan_chirps2019 = CHIRPS_daily.filterDate('2020-01-01','2020-02-01');
var feb_chirps2019 = CHIRPS_daily.filterDate('2020-02-01','2020-03-01');
var mar_chirps2019 = CHIRPS_daily.filterDate('2020-03-01','2020-04-01');
var apr_chirps2019 = CHIRPS_daily.filterDate('2020-04-01','2020-05-01');
var may_chirps2019 = CHIRPS_daily.filterDate('2020-05-01','2020-06-01');
var jun_chirps2019 = CHIRPS_daily.filterDate('2020-06-01','2020-07-01');
var jul_chirps2019 = CHIRPS_daily.filterDate('2020-07-01','2020-08-01');
var aug_chirps2019 = CHIRPS_daily.filterDate('2020-08-01','2020-09-01');
var sep_chirps2019 = CHIRPS_daily.filterDate('2020-09-01','2020-10-01');
var oct_chirps2019 = CHIRPS_daily.filterDate('2020-10-01','2020-11-01');
var nov_chirps2019 = CHIRPS_daily.filterDate('2020-11-01','2020-12-01');
var dec_chirps2019 = CHIRPS_daily.filterDate('2020-12-01','2021-01-01');
// Calculate annual precipitation (mm) for 2020 year
var annual_precip = (jan_chirps2019.sum()).add(feb_chirps2019.sum())
.add(mar_chirps2019.sum()).add(apr_chirps2019.sum())
.add(may_chirps2019.sum()).add(jun_chirps2019.sum())
.add(jul_chirps2019.sum()).add(aug_chirps2019.sum())
.add(sep_chirps2019.sum()).add(oct_chirps2019.sum())
.add(nov_chirps2019.sum()).add(dec_chirps2019.sum());
var flap = annual_precip.multiply(0.8);
Map.addLayer(flap.randomVisualizer())
var ALOS_roi = ee.ImageCollection("JAXA/ALOS/AW3D30/V3_2").select('DSM').median().clip(roi)
//Map.addLayer(ALOS_roi.randomVisualizer())
//cf_2020 for C factor; flap for R factor; ALOS_roi for elevation data
// export A value for each year (e.g. 2020);
Export.image.toDrive({
image: ALOS_roi,
description: '****',
folder:'****',
scale: 30,
crs:'EPSG:4326',
maxPixels: 1e13,
region:roi.geometry()
});
3. ArcGIS中计算LS因子
此处参考推文:GEE系列(2)——GEE处理 与arcgis 结合在土壤侵蚀上的应用
首先,通过3D Analysis中的Slope工具,利用上述下载的高程数据计算坡度roi_slope;其次,通过栅格计算器计算,公式如下:
4. K因子与P因子
K因子参考:「GIS数据」分享2023年全球土壤可蚀性数据(tif格式) | 麻辣GIS
P因子参考:忘记哪里看到的,也可以根据P = 0.2 +0.03×坡度百分比进行计算
5. 总结
最后,参考相应文献,把需要归一化的参数归一化,再将五个参数栅格数据利用ArcGIS里的栅格计算器进行计算得到结果,再根据文献中的分类等级得到研究区土壤侵蚀等级。由于数据保密,此处不放图。