【双评价学习笔记】生态—防风固沙功能重要性（1）：挑战一下模型法

之前做生态评价都是基于NPP法，今天想挑战下模型法。

==嗯，就挑战下吧，也不确定对不对。

防风固沙是生态系统（如森林、草地等）通过其结构与过程减少由于风蚀所导致的土壤侵蚀的作用，是生态系统提供的重要调节服务之一。防风固沙功能主要与风速、降雨、温度、土壤、地形和植被等因素密切相关。以防风固沙量（潜在风蚀量与实际风蚀量的差值）作为生态系统防风固沙功能的评估指标。

（1）评估模型

参考《指南6月版》以及《生态保护红线划定指南》中的模型评估法。采用修正风蚀方程来计算防风固沙量，公式如下：

SR、Sl、Sp的量纲单位为t/km2·a。Ql、Qlmax分别为潜在风沙转移量和最大转移量，Qp、Qpmax分别为实际风沙转移量和最大转移量，量纲单位为kg/m；z 为最大风蚀出现距离（m）；WF、EF、SCF、RS、Vc分别为气候因子、土壤可蚀因子、土壤结皮因子、地表糙度因子和植被覆盖因子。

公式很复杂，说到底就是要先求出WF、EF、SCF、RS、Vc这五个因子。此外还有一个z，为最大风蚀出现距离（m），参考发表在生态学报的论文《防风固沙型重点生态功能区防风固沙服务的评估与受益区识别》，本次计算取50m。

徐洁，肖玉，谢高地，王洋洋，江源，陈文辉．防风固沙型重点生态功能区防风固沙服务的评估与受益区识别．生态学报，2019，39(16) : 5857-5873．DOI:10．5846/stxb201901040042

有一些具体的指南中没给出的公式，参考了下面的文章：

张路,肖燚,郑华, 等.2010年中国生态系统服务空间数据集[J/OL].中国科学数据, 2017. (2018-01-02). DOI: 10.11922/csdata.180.2017.0145.

（2）气候因子WF

气候因子WF表征了在考虑土壤湿度、雪盖等因素下风力对土壤颗粒的搬运能力，单位为kg/m，12 个月WF总和得到多年年均WF，其计算公式如下:

Wf为各月多年平均风力因子，p为空气密度，g为重力加速度，SW为各月多年平均土壤湿度因子，无量纲；SD 为雪盖因子，无量纲。

多年平均风力因子Wf

Wf为各月多年平均风力因子。参考相关论文，特别是生态学报的那篇论文，给出了风力因子（m/s）的计算公式：

u1为2m处起沙风速，本次计算取5m/s（参考多篇论文后的取值均为5m/s）; u2为气象站2m处月均风速值(m/s) ; Nd为各月2m处风速大于5m/s的天数。

由于气象台站测得的风速数据是10m高度处的风速值，所以风蚀力模拟时所用2m处的风速值是由10m高度处的风速值转换而得，方法采用Eiliot的转换方法：

此处的u2为z2（2m）高度处的风速，此处的u1为z1（10m）高度处的风速。

现在我们就可以计算多年平均风力因子Wf了，需要的数据是两个：月均风速和日均风速大于5m/s的天数。月均风速数据来源是“国家气象科学数据中心”（中国气象数据网）的“中国地面累年值月值数据集(1981-2010年）” ；日均风速大于5m/s的天数这个数据要计算，来源还是国家气象科学数据中心，这次是要用“中国地面气候资料日值数据集(V3.0)”。与此同时，把“中国地面气候标准值数据集台站表”下载下来。

利用经纬度，将台站表中的数据导入到gis中显示，通过字段连接将数据集中的数据与台站相连。数据集中的数据可以先在Excel中简单处理后导入gis中，直接新建字段，利用字段计算器进行运算。

利用“中国地面气候资料日值数据集(V3.0)”统计各月2m处日均风速大于5m/s的天数Nd。首先通关转换计算将10米处风速转为2米处风速，然后进行汇总统计得到年均的各月2米处风速超过5m/s的天数。

接下来计算空气密度p

式中：h 为海拔高度，km；T为绝对温度，K；其中T=t+273.15℃，t为实测摄氏温度。

论文参考：于贵瑞,徐兴良,王秋凤,等.全球变化对生态脆弱区资源环境承载力的影响研究[J].中国基础科学,2017,19(6):19-23+35.

同样直接在属性表中操作即可。

重力加速度参考下表，研究区域及周边气象站纬度都40左右，高度基本小于1km，因此取9.082m/s。

接着计算p/g，并进行插值运算得到各月因子：

详细说一下计算土壤湿度因子。

土壤湿度的数据可以通过欧空局（ESA）CCI计划下载的数据来替代。

网址：https://www.esa-soilmoisture-cci.org/

由于好像没有月度合成产品，需要自己下载后进行合成。并且数据本身剔除了雪盖部分，因此某些区域存在没数据的情况。

也可用欧洲中期天气预报中心(ECMWF)的ERA-Interim的数据来替代。

网址：Https://Apps.Ecmwf.Int/Datasets/Data/Interim-Full-Mnth/Levtype=Sfc/

也可以自己计算==

土壤湿度因子SW表达式如下：

式中，R为各月平均降雨量，单位为mm；I为灌溉量，本次取0；Rd 为各月平均降雨天数；Nrd为各月降雨平均天数；ETP为各月平均潜在蒸发量，单位为mm。ETP表达式如下：

式中，SR为各月平均太阳总辐射，单位为cal/cm2；DT为各月平均温度，单位为℃。

首先需要获取各月平均太阳总辐射的数据，在“资源学科创新平台”获取“历年各月太阳总辐射”的数据，保存为Excel格式，然后加载到ArcGIS中进行插值。

网址：http://www.data.ac.cn/

数据为全国台站数据，进行大范围插值后提取出研究区的栅格数据。

接着计算月平均温度，结合各月太阳总辐射，算出各月ETP：

最后计算各月月均降雨量和降雨天数，利用公式，算出各月土壤湿度因子SW：

各月多年平均雪盖因子SD

SD为雪盖因子，即无积雪覆盖天数/研究总天数（定义雪盖深度＜25.4mm 为无积雪覆盖）;也有文章定义为20.0.mm的。我这边引用的数据为“国家冰川冻土沙漠科学数据中心”的“中国地区MODIS雪盖产品数据集”，该数据集为2000年2月24日至2004年12月31日，分辨率为0.05度，MODIS数据，数据格式为.hdf，可以用HDFView打开，也可以直接在ArcGIS打开。像素值0代表无雪，1-100代表每个单元的积雪比例。因此在使用时，可以换算成小数，用1减去原始值后得到雪盖因子。

数据来源网址：

http://www.crensed.ac.cn/portal/metadata/6bdfe932-8c3e-4747-9600-33e578aaf12d

计算气候因子WF

根据公式计算每个月WF，求和得到总WF。

（3）土壤可蚀因子EF

在一定的物理化学环境下，土壤受到风力侵蚀的难易程度即为土壤可蚀性因子。结合《双评价指南6月版》和《生态保护红线指南》，土壤可蚀性因子EF的计算公式为：

式中，sa为土壤粗砂含量（0.2mm～2mm）（%）；si为土壤粉砂含量（%）；cl为土壤粘粒含量（%）；OM为土壤有机质含量（%）；Caco3为碳酸钙含量（%），可不予考虑。

土壤数据使用的是“世界土壤数据库（HWSD）土壤数据集(v1.2)”

文章的引用 Wieder,W.R., J. Boehnert,G.B. Bonan, andM. Langseth.2014. RegriddedHarmonized World Soil Database v1.2. Data set. Available on-line[http://daac.ornl.gov]：fromOak Ridge National Laboratory Distributed Active Archive Center, Oak Ridge,Tennessee,USA.

http://dx.doi.org/10.3334/ORNLDAAC/1247 .

有两个方案可以采取，一是直接用栅格进行连接，二是栅格转成矢量文件进行连接。个人觉得转为矢量文件更好操作。转成矢量后直接在属性表中计算出EF，再转成栅格。

回顾一下我们需要哪几个数据：

（4）土壤结皮因子SCF

土壤表面在物理作用或者生物作用下，形成的特殊土壤坚硬表层结构能够有效抵抗外力作用，减轻土壤风蚀程度，公式为：

式中，cl为土壤粘粒含量（%）；OM为土壤有机质含量（%）。

属性表中操作，然后转为栅格。

还有个地表糙度因子、植被覆盖因子，以及最后的计算==下次再说吧~

有兴趣给个赞赏吗？

END>

如需转载文章、投稿或者合作

请直接在公众号后台联系

使用关键词搜索历史文章请点击【阅读原文】