产水量计算

发表于 2020-10-20 更新于 2023-10-04 分类于 Ecosystem service ， Annual Water Yield 阅读次数： Waline：本文字数： 5.7k 阅读时长 ≈ 21 分钟

本文由SCY原创，转载注明出处。

产水量模型

InVEST模型的产水量：水库水电生产模块，不仅评估了景观中各次一级流域对产水的相对贡献量，而且有助于研究土地利用格局变化如何影响年地表产水量和水电生产。

模拟景观格局变化和水文过程的关系是个科学难题。用来模拟这些关系及其相关过程的复杂模型（如WEAP水资源评估和规划系统模型）要求较高的数据和资源资料，并且要求大量专业知识。为满足更多学科背景需求，使用易获取的数据，InVEST模型绘制和模拟用于景观水电生产的年平均产水量，而不是直接评估土地利用/覆被变化导致的水电减产，因为这一过程在逐日和逐月的时间尺度上主要由流入水量变化决定。因此，InVEST模型计算了景观中各子流域对产水的相对贡献量及其对应的水电产量。水库经济使用年限内的水电生产净现值也可以通过年收益还原法计算。

工作原理

模型运行基于栅格地图。模型估算了研究区各次一级流域对水电生产贡献水量及其经济价值。模块包括三个组件，按顺序运行。

模型估算了每栅格单元降水量减去实际蒸散发后的水量即水源供给量。模型不做地表水、地下水、基流的区分，而是假设每个栅格单元的产水通过上述途径汇集到信息点。然后，模型计算出次一级流域产水量的总量和平均值。栅格计算有助于确定决定流域产汇流空间异质性的关键因素，如土壤类型，降水量，植被类型等。但是，这组模型的基础理论基于次一级流域到流域尺度，对次一级流域过程的模型解释是可信的，因此产水量的总量和/或平均值结果也应当维持在次一级流域水平上。模型仅为校准和模型检验提供输出结果的栅格数据。这些栅格数据图件不能用于水文过程的解释说明，或作为任何类型的决策信息源。
模型计算了用于水电生产的水源供给量，即将水源供给总量减去除水电生产以外的其他用水量。
模型计算了到达水库水流的发电量及其水库有效使用年期内的经济价值。

NOTE: Water Yield 模块运行的前提条件是假设栅格单元的产水量都是通过地表径流或者地下径流的方式汇集到流域出口，在这个前提条件下计算每个栅格单元的产水量, 即降水量减去植被蒸腾与地表蒸散。

产水量评估模型

产水量评估模块基于Budyko水热耦合平衡假设（1974）和年平均降水量数据。首先，确定研究区每个栅格单元$x$的年产水量$Y(x)$，公式如下：

$$
Y(x)=\left(1-\frac{A E T(x)}{P(x)}\right) \cdot P(x)
$$

$AET(x)$表示栅格单元$x$的年实际蒸散量、$P(x)$表示栅格单元$x$的年降水量。

水量平衡公式中，土地利用/覆被类型的植被蒸散发$\frac{A E T(x)}{P(x)}$计算，采用Fu和Zhang等提出的Budyko水热耦合平衡假设公式：

$$
\frac{A E T(x)}{P(x)}=1+\frac{P E T(x)}{P(x)}-\left[1+\left(\frac{P E T(x)}{P(x)}\right)^{\omega}\right]^{1 / \omega}
$$

$PET(x)$表示潜在蒸散量、$ω(x)$表示自然气候-土壤性质的非物理参数。
潜在蒸散量$PET(x)$定义为：

$$
{PET}(x)=K_{c}\left(\ell_{x}\right) \cdot E T_{0}(x)
$$

式中，$ET_0(x)$表示栅格单元$x$的参考作物蒸散，$K_c(\ell_x)$表示栅格单元$x$中特定土地利用/覆被类型的植物（植被）蒸散系数。$ET_0(x)$通过参考作物蒸散量反映当地气候条件，例如苜蓿的蒸散量反映其草地生境气候。$K_c(\ell_x)$很大程度上取决于栅格单元$x$中土地利用/覆被的植被性质。在土地利用/覆被图中，$K_c$用于将$ET_0(x)$修正为栅格单元中特定作物或植被类型蒸散量。

$ω(x)$是一个经验参数，通常用$\frac{A W C \times N}{P}$线性函数表示，式中$N$表示每年的降水事件数，$AWC$表示植物可利用水含量。虽然基于全球数据的$ ω(x) $公式亟待进一步研究， InVEST模型采用Donohue等人提出的公式表达，定义为：

$$
$$

$AWC(x)$表示土壤有效含水量（mm），由土壤质地和土壤有效深度决定，用来确定土壤为植物生长储存和提供的总水量。由植物利用水分含量$(PAWC)$，以及土壤的最大根系埋藏深度和植物根系深度的最小值决定：

$$
A W C(x)=\operatorname{Min}( Re st.layer.depth, root.depth ) \cdot PAWC
$$

土壤的最大根系埋藏深度是指由于环境的物理和化学特征不同，植物根系在土壤中能够延伸的最大深度（也叫土壤深度）。植物根系深度通常指特定植物类型的根系生物量为95%的土层深度。$PAWC$表示植物利用水分含量，即田间持水量和萎蔫点之间的差值。

$Z$为经验常数，又称季节常数，能够代表区域降水分布及其他水文地质特征。$Z$与$N$正相关，$N$是每年降水发生次数。1.25为$ω(x)$基数，即裸地（根系深度为0）的植被年需水量和年降水量比值。$ω(x)$上限为5。

其他土地利用/覆被类型（开放水域，城市，湿地）的实际蒸散发通过参考作物蒸散$ET_0(x)$直接计算，由降水量决定其最大值：

$$
A E T(x)=\operatorname{Min}\left(K_{c}\left(\ell_{x}\right) \cdot E T_{0}(x), P(x)\right)
$$

$ET_0(x)$表示参考作物蒸散，而$K_c(\ell_x)$表示特定土地利用/覆被类型蒸腾作用的影响因子。

数据需求

模型使用的数据需求列表如下，关于数据来源和预处理的详细信息见附录。所有数据输入前，应先定义栅格数据投影，栅格单位为米（m）。

土壤的最大根系埋藏深度（必需）：

每个栅格对应一个土壤的最大根系埋藏深度平均值的GIS栅格数据集。土壤的最大根系埋藏深度是指由于环境的物理和化学特征不同，植物根系在土壤中能够延伸的最大深度，单位毫米（mm）。根系限制层深度可从一些土壤图中获得。如果无法获得根系限制层深度或按土壤类型划分的可扎根深度，可使用土壤深度作为替代。如果有几个土壤层是详细的，那么限制根系层的深度就是非限制性土壤层的深度之和。

*命名：*用户自定义，但若为ESRI GRID格式，文件名不能有空格并且少于13个字，若为TIF或IMG格式，命名可能更长。

格式： GIS标准栅格文件（如：ESRI GRID，TIF或IMG），每个栅格对应一个土壤的最大根系埋藏深度平均值，单位毫米（mm）。

数据获得：

http://globalchange.bnu.edu.cn/research/cdtb.jsp （单位为m，需转化为mm，已下载处理，上传百度☁️）
原文 https://doi.org/10.1038/s41597-019-0345-6
SoilGrids250m 2017-03 - Depth to bedrock (R horizon) （单位为cm，需转化为mm）https://data.isric.org/geonetwork/srv/eng/catalog.search#/metadata/bfb01655-db81-4571-b6eb-3caae86c037a

年降水量（必需）：

每个栅格对应一个非空值的年平均降水量的GIS栅格数据集，单位毫米（mm）。

命名：用户自定义，但若为ESRI GRID格式，文件名不能有空格并且少于13个字，若为TIF或IMG格式，命名可能更长。

格式：GIS标准栅格文件（如：ESRI GRID，TIF或IMG），每个栅格对应一个年平均降水量。

数据获得：

利用日值（v3）数据进行差值。
国家地球系统科学数据中心：中国1km分辨率年降水量数据（2001-2020年）已下载直接使用，单位为0.1mm http://www.geodata.cn/data/datadetails.html?dataguid=113786088533256&docid=3549
WorldClim (v 2.1) 数据范围：1970-2000，月值，最大精度：30s ，单位:mm https://www.worldclim.org/data/worldclim21.html
GEE （“UCSB-CHG/CHIRPS/DAILY”)数据集，完整教程见编程分享-利用GEE得到研究区年、月、日降水量数据

var ROI = ee.FeatureCollection("users/rice20220411/AH");
Map.addLayer(ROI,{},'ROI');
for(var i = 2015;i<=2021;i++){
var CHIRPS_Daily = ee.ImageCollection("UCSB-CHG/CHIRPS/DAILY")
.filterDate(i+'-01-01', i+'-12-31')
 .select('precipitation')
var CHIRPS_Year_mean = CHIRPS_Daily.mean().clip(ROI)
var precipitationVis = {
  min: 1.0,
  max: 17.0,
  palette: ['001137', '0aab1e', 'e7eb05', 'ff4a2d', 'e90000'],
};
print(CHIRPS_Year_mean)
Map.addLayer(CHIRPS_Year_mean, precipitationVis, i+'_CHIRPS_Year_mean');
// Map.addLayer(CHIRPS_Daily.first().clip(ROI), precipitationVis, 'CHIRPS_Year_mean_first');
Export.image.toDrive({
      image: CHIRPS_Year_mean,
      description: i+'year_mean',
      region: ROI,
      maxPixels: 1e13,
      folder: 'CHIRPS'
    }) 
}

植物可利用水量（必需）：

每个栅格对应一个植物可利用水的GIS栅格数据集。植物可利用水（PAWC）是指土壤土层中为植物生长提供的水量所占比例。PAWC是[0,1]的小数。

命名：用户自定义，但若为ESRI GRID格式，文件名不能有空格并且少于13个字，若为TIF或IMG格式，命名可能更长。

格式：GIS标准栅格文件（如：ESRI GRID，TIF或IMG），每个栅格对应一个植物可利用水含量百分比。

数据获得：（下载好处理好后传百度☁️）

SoilGrids250m 2017-03 - “Derived available soil water capacity (volumetric fraction) until wilting point” (https://data.isric.org/geonetwork/srv/eng/catalog.search#/metadata/e33e75c0-d9ab-46b5-a915-cb344345099c) SoilGrids 2.0版目前不提供AWC。SoilGrids 2017提供7个土壤深度区间的AWC层。所有7个深度区间都需要下载，然后合并成一个单一的图层，以便在模型中使用。

栅格值以整数百分比的形式给出（如25，表示AWC值为25%）。
（标准）深度区间的平均值，如0-5厘米或0-30厘米，可以通过使用数字积分，如梯形规则，对深度区间内的预测值进行加权平均来得出：

$$
\left(\frac{1}{(b-a)}\right)\left(\frac{1}{2}\right) \sum_{k=1}^{N-1}\left(x_{k+1}-x_{k}\right)\left(f\left(x_{k}\right)+f\left(x_{k+1}\right)\right)
$$

其中，$N$是深度数，$x_k$是第$k$个深度，$f(x_k)$是目标变量（即土壤属性）在深度$x_k$的值。

操作步骤：

从ISRIC网站上下载所有可用的深度区间。深度区间为0cm-200cm。注意，每个栅格的大小为1.5GB。
使用 GIS 缓冲区工具，在你要建模的流域/感兴趣的区域周围创建一个缓冲区。由于SoilGrids数据的分辨率为250米，因此缓冲区的宽度为250或500米。这样做是为了确保土壤数据完全覆盖你正在建模的流域，边界周围没有漏洞。
使用缓冲流域，将所有原始的ISRIC AWC 栅格数据剪辑到你感兴趣的区域。在ArcGIS中，这可以通过空间分析工具掩模提取来完成。在这个例子中，我们将把剪下的图层称为AWC_sl1_clip.tif、AWC_sl2_clip.tif … AWC_sl7_clip.tif。
使用GIS栅格计算器工具来计算合并后的AWC层。将其代入上面的方程，我们可以得到：
$$(1/(200-0)) * (1/2) * ( ((5-0) * (AWC_sl1_clip.tif + AWC_sl2_clip.tif)) + ((15-5) * (AWC_sl2_clip.tif + AWC_sl3_clip.tif)) + ((30-15) * (AWC_sl3_clip.tif + AWC_sl4_clip.tif)) + ((60-30) * (AWC_sl4_clip.tif + AWC_sl5_clip.tif)) + ((100-60) * (AWC_sl5_clip.tif + AWC_sl6_clip.tif)) + ((200-100) * ( AWC_sl6_clip.tif + AWC_sl7_clip.tif)) )$$
将此公式输入光栅计算器，根据需要调整文件名。
得到的栅格应该包含0-100范围内的数值，代表整数百分比。该模型要求AWC以分数的形式给出，因此将步骤4中计算的栅格除以100。
重新投影AWC分数层，使其具有与其他模型输入相同的投影坐标系统。这个栅格现在可以用作模型的可用水含量输入。

年平均潜在蒸散发（必需）：

每个栅格对应一个年平均潜在蒸散发的GIS栅格数据集。潜在蒸散发是指水分充足的情况下，通过土壤蒸发和植物（如苜蓿或其他草类等健康植被）蒸散作用可能散逸的水量，单位毫米（mm）。

命名：用户自定义，但若为ESRI GRID格式，文件名不能有空格并且少于13个字，若为TIF或IMG格式，命名可能更长。

格式：GIS标准栅格文件（如：ESRI GRID，TIF或IMG），每个栅格对应一个年平均潜在蒸散发。

数据获得：

FA0工具计算，详见使用EToCalculator 计算潜在蒸散发详细教程 https://blog.sciencenet.cn/blog-3459054-1308734.html
Global Aridity Index and Potential Evapotranspiration (ET0) Climate Database v2 (Penman Monteith Evapotranspiration equation，30弧秒，1970-2000，暂时未使用，有待验证。) https://figshare.com/articles/dataset/Global_Aridity_Index_and_Potential_Evapotranspiration_ET0_Climate_Database_v2/7504448/3
大部分研究采用，Modified-Hargreaves 法对潜在蒸散量进行计算。
式中：$ET_0$代表栅格单元上的潜在蒸散量（mm）；$RA$为太阳大气顶层辐射（$MJ·m^{-2}·d^{-1}$）；$T_a$为不同小流域的日最高温均值和日最低温均值的平均值（°C）； $T_d$为小流域上的日最高温均值和日最低温均值的差值（°C）；$P$为小流域的降雨量（mm）。

$$
ET_{0}=0.0013 \times 0.408 \times R A \times\left(T_{\mathrm{a}}+17\right) \times\left(T_{\mathrm{d}}-0.0123 P\right)^{0.76}
$$

GEE获取 (Penman-Monteith法)（代码传GEE和GitHub）已经下载全国区域的2000-2020年数据压缩上传百度云

原文链接 https://doi.org/10.1016/j.rse.2018.12.031

//完美运行，输入坐标即可，获取每年的平均值。
//加入矢量边界之后直接运行var 
var ROI = table.geometry();


//Map.addLayer(ROI,{},'ROI');
var imgcol_PMLV2_v016_8d = ee.ImageCollection('projects/pml_evapotranspiration/PML/OUTPUT/PML_V2_8day_v016');
var imgcol_PMLV2_v016_8d = ee.ImageCollection(imgcol_PMLV2_v016_8d.toList(2000));
print(imgcol_PMLV2_v016_8d);

var pkg_export = require('users/kongdd/pkgs:pkg_export.js');
var options = {
    type: "drive",
    range: [110,34,114,40],//ROI, // [73, 25, 105, 40], 
    cellsize: 1/240,
    // crsTransform : [463.312716528, 0, -20015109.354, 0, -463.312716527, 10007554.677], // prj.crsTransform;
    // scale        : 463.3127165275, // prj.scale
    crs: 'EPSG:4326', // 'SR-ORG:6974', // EPSG:4326
    folder: 'ET0'
};


var imgcol_years = aggregateToYearly(imgcol_PMLV2_v016_8d, 2010, 2019, 0.01);
function aggregateToYearly(imgcol, year_begin, year_end, scale_factor) {
    var bands = ['ET_pot']; //,'qc'
    var years = ee.List.sequence(year_begin, year_end);
    var imgcol_years = years.map(function(year) {
        var date_begin = ee.Date.fromYMD(year,1,1);
        var date_end   = ee.Date.fromYMD(year,12,31);
        var ydays = date_begin.advance(1, 'year').difference(date_begin, 'day');
        var imgcol_year = imgcol.filterDate(date_begin, date_end);
        var scale = ydays.multiply(scale_factor);
        return imgcol_year.select(bands)
            //.multiply(scale_factor).
            .mean()
            .multiply(scale)
            .toFloat()
            .set('system:time_start', date_begin.millis())
            .set('system:id', date_begin.format());
    });
    imgcol_years = ee.ImageCollection(imgcol_years);
    // pkg_export.ExportImg(img_year, task, range, cellsize, type, folder_yearly, crs, crsTransform);
    return imgcol_years;
}
var imgcol_years = imgcol_years.filterBounds(ROI);
print(imgcol_years);
//var batch = require('users/fitoprincipe/geetools:batch');
//batch.Download.ImageCollection.toDrive(imgcol_years,"ETO", {
//scale: 1/240})


pkg_export.ExportImgCol(imgcol_years.limit(10), 'ET0', options);

GEE批量导出方法：

    var runButtons = document.querySelector('#task-pane').shadowRoot.querySelectorAll(".run-button")
     runButtons.forEach(function(e) {e.click()})
}
runTaskList()
setTimeout( 
    function(){
       var taskDialog = document.querySelectorAll("ee-image-config-dialog") //table的话-image-改成-table-
        taskDialog.forEach(function(e) {e.shadowRoot.querySelector("ee-dialog").shadowRoot.querySelector("paper-dialog").querySelector(".ok-button").click()})
 },5 * 1000 );

土地利用/覆被（必需）：

每个栅格对应一个土地利用类型的GIS栅格数据集。土地利用类型代码定义为整数。

*命名：*用户自定义，但若为ESRI GRID格式，文件名不能有空格并且少于13个字，若为TIF或IMG格式，命名可能更长。

*格式：*GIS标准栅格文件（如：ESRI GRID，TIF或IMG），每个栅格对应一个土地利用类型代码（如：1表示森林，3表示草地，等）。地类代码必须与生物物理系数表中的地类代码一致。

数据获取得：

Land cover classification gridded maps from 1992 to present derived from satellite observations https://cds.climate.copernicus.eu/cdsapp#!/dataset/satellite-land-cover?tab=overview
ESA全球10米https://esa-worldcover.org/en
GLOBELAND30 包含2000，2010，2020三期数据。http://www.globallandcover.com/defaults.html?src=/Scripts/map/defaults/download.html&head=download&type=data
note:汾河流域在图幅N49_35,山西省全域需要图幅N49_30,N49_35,N49_40,N50_35,N50_40
Sentinel-2 Land Use/ Land Cover Downloader 2017-2021年数据，10米分辨率。 https://www.arcgis.com/apps/instant/media/index.html?appid=fc92d38533d440078f17678ebc20e8e2
1990-2021全国30米数据（每年更新） https://zenodo.org/record/5816591#.YzQUF-xBwbk

流域（必需）：

用多边形表示流域的图形文件（shapefile）。即与研究区水电生产研究相关的所有小流域图层。

命名： 用户定义，但文件名不能有空格

*格式： *图形文件（.shp）

*属性表横列：*每一行表示一个小流域

*属性表纵列：*必须包含定义为整数的“ws_id”字段，每个小流域赋予唯一数值。

次一级小流域（必需）：

流域图层中选定区域中用多图形文件（shapefile）。生成次一级小流域的工具和方法，详见 “DEM数据处理”章节。

*格式：*图形文件（.shp）

*属性表横列：*每一行表示一个次一级小流域

*属性表纵列：*必须包含定义为整数的“subws_id”字段，每个次一级小流域赋予唯一数值。

生物物理系数表（必需）：

土地利用/覆被（LULC）类型表，包括用于该工具使用的生物物理系数数据。注意事项：这些数据主要针对每种土地利用类型属性而非栅格图栅格单元属性。

*命名：*文件名由字母、数字和下划线组成，不能有空格。

*格式：*ArcGIS模型的*.dbf或*.mdb格式，独立运行模块要求一个*.csv文件。

*属性表横列：*每一行表示一个土地利用类型。

*属性表纵列：*每一列包含每种土地利用类型的不同属性，属性命名如下：

Lucode（土地利用类型代码）：每种土地利用类型地类代码（如：1表示森林， 3表示草地，等），必须与上述土地利用类型栅格图保持一致。

LULC_desc：土地利用类型的描述性命名（可选填）。

LULC_veg：包括使用的实际蒸散发AET计算公式。植被覆盖地类（不包括湿地）赋值为1，其他土地利用类型（包括湿地、城市用地、水体）赋值为0。

root_depth：植被覆盖地类的最大根系深度，单位毫米（mm），取整数。植物根系深度通常指特定植物类型的根系生物量为95%的土层深度。对不适用一般 Budyko干燥指数（即应使用公式2计算实际蒸散发AET）的土地利用类型而言，不需要根系深度数据，设为N/A。

Kc ：每种土地利用类型的植物蒸散系数，通过将植物生理学特性与苜蓿相比较，将苜蓿的参考作物蒸散修正为特定土地利用类型的潜在蒸散量。因此土地利用类型的植物蒸散系数取值为[0，1.5]的小数（在某些非常潮湿的热带区域，水分充足）。

季节常数Z（必需）：

是根据季节性降水分布定义的从1 到30排序浮动值。

用水需求表（必需）：

土地利用类型表，表示不同土地利用类型的消耗性用水量。消耗性用水量是指提供植物和作物生长，被人类和畜牧消耗，或其他应从流域水量平衡中扣除的水量。

*命名：*文件名由字母、数字和下划线组成，不能有空格。

*格式：*ArcGIS模型的*.dbf或*.mdb格式，独立运行模块要求一个*.csv文件属性表横列：每一行表示一个土地利用类型，并且必须包含土地利用栅格图中所有土地覆被属性值。

*属性表纵列：*每一列包含每种土地利用类型的不同用水需求属性，属性命名如下：

lucode（土地利用类型代码）：每种土地利用类型地类代码（如：1表示森林， 3表示草地，等），必须与上述土地利用类型栅格图保持一致。

demand：每种土地利用类型的预测平均消耗性用水量。土地利用类型图中的用水量用立方米/年/栅格单元表示。注意事项：由于区域越大，相同土地覆被类型消耗的水量可能，因此用水量的栅格计算方法十分重要。

结果分析

output\per_pixel\fractp（分数）：模型估算每个栅格单元降水量的实际蒸散发占比（实际蒸发量/降水量）。这是栅格单元的实际蒸散发占降水量的平均值。

output\per_pixel\aet（mm）：模型估算栅格单元实际蒸散发。

output\per_pixel\wyield（mm）：模型估算栅格单元产水量。

output\subwatershed_results_wyield.shp 和 output\subwatershed_results_wyield.csv：包含模型估算次一级流域生物物理学参数值的shapefile文件和表格，属性包括：
- precip_mn（mm）：次一级流域栅格单元的平均降水量。
- PET_mn（mm）：次一级流域栅格单元的平均潜在蒸散发。
- AET_mn（mm）：次一级流域栅格单元的平均实际蒸散发。
- wyield_mn（mm）：次一级流域栅格单元的平均产水量。
- num_pixel：次一级流域栅格单元的数量。
- wyield_vol（$ m^3 $）：次一级流域产水量体积。
- wyield_ha（$ m^3 $）：次一级流域每公顷产水量体积。
output\watershed_results_wyield.shp 和 output\watershed_results_wyield.csv ：包含模型估算每个小流域生物物理学参数值的shapefile文件和表格：运行产水量评估模块时，输出结果的生物物理学参数如下：
- precip_mn（mm）：每个小流域栅格单元的平均降水量。
- PET_mn（mm）：每个小流域栅格单元的平均潜在蒸散发。
- AET_mn（mm）：每个小流域栅格单元的平均实际蒸散发。
- wyield_mn（mm）：每个小流域栅格单元的平均产水量。
- num_pixel：每个小流域栅格单元的数量。
- wyield_vol（$ m^3 $）：每个小流域内产水量体积。
- wyield_ha（$ m^3 $）：每个小流域每公顷产水量体积。