☰
InVEST产水模型基岩深度栅格数据获取与预处理全流程指南
2026/9/30 3:41:35 网站建设 项目流程

做生态系统服务评估的人,大概率都绕不开InVEST。尤其是产水模型(Water Yield),几乎所有水源涵养、水生态功能评价项目里都会碰到它。可很多人在数据准备阶段就被一个看似不起眼的参数卡住:root restricting layer depth,也就是根系限制层深度。这个参数本质上反映的是土壤剖面里阻碍根系向下发展的那一层离地表有多远,绝大多数情况下可以直接用基岩深度栅格数据来表征。我项目里常用的是一套覆盖中国全境、250米分辨率的基岩深度栅格,数据源自ISRIC World Soil Information的SoilGrids产品,经过裁剪、重投影、单位换算之后,能非常顺滑地接入InVEST产水模型。这篇文章就把从数据源头、预处理到模型配置的完整链条拆开讲清楚,顺便把ArcMap里栅格数据导出Excel这种高频操作一并交代。

这套数据的价值不只是"填一个参数"那么简单。基岩深度会直接影响模型对土壤储水、实际蒸散和产水量的模拟,差之毫厘,结果就可能从"水源涵养功能高"变成"中等偏低",进而影响生态补偿、国土空间规划里的决策判断。所以别小看这张栅格图,跑模型之前把它处理明白,比纠结Z参数怎么调更值得投入时间。

1. 为什么产水模型离不开基岩深度栅格数据

1.1 InVEST产水模型的输入清单里,最容易被低估的一环

打开InVEST里的Annual Water Yield模型,摆在你面前的是这么一串输入项:年均降水量栅格、参考蒸散量栅格、土地利用/覆被栅格及配套属性表、根系限制层深度栅格、植物可利用水分比例栅格、流域边界矢量,外加一个Z常数。

大多数教程会花大量篇幅讲土地利用怎么分类、Kc系数怎么查文献、Z值怎么率定,很少有人认真对待根系限制层深度这一项。不少人随手填一个常数,或者干脆拿土壤厚度数据顶上去,运气好能跑通,但结果是否贴近真实情况就不好说了。实际上,这一项直接进入模型的水量平衡方程,决定了土壤剖面中可用于植被蒸散的水分总量。基岩深度数据能在多大程度上影响产水模拟,下面这一段展开讲。

1.2 基岩深度在产水量计算中的真实作用

InVEST年产水模型的核心公式并不复杂:

Y(x) = (1 - AET(x) / P(x)) × P(x)

其中AET是实际蒸散量,P是降水量。模型用Budyko曲线估算AET/P的比值,而这个比值的形状受一个无量纲参数ω控制:

ω(x) = Z × AWC(x) / P(x) + 1.25

注意这个AWC,植物可利用水含量。在模型里,它由植物可利用水比例(PAWC,通常0到1之间)和根系限制层深度共同决定。根系限制层越深,土壤剖面能存储的可用水越多,旱季植物就有更充足的水分维持蒸散,AET/P比值相应偏高,最终算出来的产水量就会偏低。反过来,基岩很浅的地区,土壤储水空间有限,降水一来就形成径流,产水系数自然更高。

这里有个反直觉的细节值得单独说:基岩深度并不是"越深越好",也不是"越深产水越多"。它对产水量是间接作用——通过调节土壤水库的容量来影响蒸散过程。你甚至可以把它想象成一个蓄水池:基岩深度决定了池子的深度,PAWC决定了池子里水的“有效比例”,降水量是进水管,蒸散是出水管。池子越深,旱季能供给植物的水越多,出水管持续开着,最终从溢流口排出的水就越少。理解了这张比喻,你就明白为什么产水模型特别看重基岩深度栅格数据的质量了。

1.3 250米分辨率:全国尺度评估的“甜点”

既然基岩深度这么重要,分辨率怎么选?250米这个尺度,不是拍脑袋定的。

首先是和数据源匹配。ISRIC SoilGrids原生产品就是250米分辨率,直接用这套栅格,省去大量重采样的麻烦和由此引入的误差。重采样不是不能做,只是每做一次就会丢失一部分空间异质性,尤其在基岩起伏明显的山区,硬把250米升采样到90米,并不会带来额外真实信息,只是增加了文件体积和运算时间。

其次,250米在全国尺度上恰好处于“够细”和“算得动”的交界处。一个中国全境的250米栅格,像元数大约在千万级,对InVEST这种单机运行的模型来说是能扛住的。如果用30米分辨率,文件动辄几个GB,裁剪、重投影、读取每一步都磨人,而且大多数流域尺度的生态服务评估也用不到这么细。反过来,如果用1公里分辨率,很多峡谷、断层带附近的浅基岩区域会被平滑掉,产水空间格局会失真。

所以我的建议很直接:如果你做的是省级、国家级或者大流域范围的产水评估,250米基岩深度栅格就是性价比最高的选择,没有之一。

2. 数据源头与处理链条:从ISRIC到可用的中国栅格

2.1 ISRIC World Soil Information能提供什么级别的数据

先说ISRIC。全称International Soil Reference and Information Centre,国际土壤参考和信息中心,跟联合国粮农组织、各国土壤调查机构都有长期合作。普通人接触它,多半是因为SoilGrids——一个全球土壤属性网格数据集。

SoilGrids基于全球大量土壤剖面实测数据和协变量(气候、植被、地形、母岩等),用机器学习模型预测出多种土壤属性的空间分布,分辨率250米。我们要的基岩深度,在SoilGrids里的正式名称通常是Depth to Bedrock,不同版本字段代码可能有差异,有的叫bdticm,有的直接写成depth_to_bedrock,单位基本都是厘米,值域一般被截断到0到200厘米。也就是说,凡是基岩深度超过2米的地方,产品返回的数值大多是200。这个截断对产水模型的影响其实有限,因为绝大多数农作物和草本的根系深度远不到2米,乔木类虽然根系深,但真正参与水分吸收的活跃根层也多在上层1米到1.5米以内。用200厘米作为上限值,在生态服务评估里完全可以接受。

现在想获取SoilGrids数据,路径很多。新手可以直接去官网的交互式地图逐块下载,懂一点自动化的话,用REST API或WCS服务按经纬度范围请求更高效。下载的时候注意选对图层和坐标系,SoilGrids默认是WGS84地理坐标系,后续到中国区域使用一般要重投影。

2.2 全球栅格裁剪到中国区域的完整操作

拿到全球或亚洲范围的深度到基岩栅格后,常规处理流程分四步:裁剪、投影、换算单位、统一NoData。这一步在ArcMap或QGIS里都能完成,我这里用命令行工具的写法给你一套可复现的流程,方便批量处理。

# 1. 先按中国范围的经纬度范围裁剪,大致范围可设为 # 东经73至135度、北纬18至54度,留出一点余量 gdal_translate -projwin 73 54 135 18 -of GTiff \ depth_to_bedrock_raw.tif china_dtb_raw.tif # 2. 重投影到适合中国区域的等面积投影,比如Albers Conical Equal Area # 目标EPSG可选用自定义或当地常用投影,关键是保证后续所有栅格一致 gdalwarp -t_srs "+proj=aea +lat_1=25 +lat_2=47 +lat_0=0 +lon_0=105 +datum=WGS84" \ -r bilinear china_dtb_raw.tif china_dtb_aea.tif # 3. 单位换算:SoilGrids基岩深度单位是cm,InVEST产水模型要求mm # cm转mm,乘10即可。这一步简单但极其容易漏掉。 gdal_calc.py -A china_dtb_aea.tif \ --outfile=china_dtb_250m_mm.tif \ --calc="A*10" --NoData=-9999

这里有几个细节必须交代。

第一,裁剪时最好留出缓冲区,别把国境线刚好切在像元边界上,否则后续跟流域边界叠加时边缘容易出细缝。第二,重采样方法我建议用双线性,基岩深度是连续变量,最近邻法会让地形台阶感特别强,尤其在山谷区域会有明显的锯齿。第三,投影选择上,如果项目后续还要算面积、算体积,建议用Albers或UTM这类等面积/等距离投影,别用WGS84经纬度直接交到InVEST里跑,虽然有时候能跑通,但模型在计算像元面积的时候会基于投影后的几何结果,投影不合适会引入系统偏差。

2.3 拿到栅格后的三步质检

数据别急着用。我处理这套中国基岩深度数据时,每次拿到新的版本都会做三件事,十分钟不到,能避免后面跑模型时很多莫名其妙的报错。

第一,统计极值。用ArcMap的栅格统计工具或Python的numpy读一遍,正常基岩深度范围应该在0到200厘米,换算成毫米就是0到2000。如果出现负值,基本可以断定是浮点运算时的NoData边界处理出了问题;如果出现几千的巨值,就要怀疑是不是数据拼接时没处理好。

第二,检查NoData空洞。SoilGrids在积雪覆盖、水体、难以取样的区域经常会有缺失值。这些空洞如果在流域内部,InVEST会直接报错或者把该像元当0处理,导致产水量出现黑洞状异常值。我的做法是先做一个空洞掩膜,再用周边像元填补,或者干脆把空洞区域标记出来,最后跟流域边界比对看是否会造成实质影响。

第三,地形逻辑校验。基岩深度和地形有强相关性:高山陡坡处基岩通常浅,冲积平原和谷地通常深。随便抽几个点,叠加DEM看一眼。如果发现高山顶上写着200厘米,河边反而是0厘米,就要小心投影偏移或者数据源波段选错了。这一步能帮你拦截掉至少一半的低级错误。

3. 基岩深度数据接入InVEST产水模型的实操配置

3.1 单位换算与NoData归一

上一节提到了厘米到毫米的换算,这里单独拎出来说,是因为我见过太多次因为单位问题导致结果量级离谱的例子。InVEST产水模型里,降水量的单位是毫米/年,根系限制层深度的单位也是毫米。你把SoilGrids的基岩深度直接填进去,数值会差10倍。表面看产水量还能跑出来,但模型内部在算植物可利用水总量时,AWC直接被低估或高估10倍,ω指数跟着偏移,最终产水量的空间格局可能完全变形。

NoData同样要在进模型前统一。InVEST对NoData有自己的约定,一般要求所有栅格的NoData一致,且在属性设置里明确指定。我习惯把NoData设成-9999,在gdal_calc里用--NoData=-9999写好,这样不管是后续在ArcMap里检查,还是拿到Python里做统计分析,都不容易误伤有效像元。

还有一点,InVEST要求输入的所有栅格:范围一致、投影一致、分辨率一致。这三条看起来是常识,实际项目里却总是翻车。降水数据可能是1公里分辨率,蒸散数据可能是500米,基岩深度是250米,土地利用是30米。我的经验是先把所有栅格统一到同一个基准栅格上,用栅格重采样工具配准。基准建议就选250米基岩深度栅格,省一次重采样。

3.2 InVEST面板参数填写实操

在InVEST图形界面里,模型列表选择Water Yield (Annual),然后会出现完整的参数面板。我按顺序讲一遍我的填法,都是实测过的。

第一个是降水栅格。年均降水量,单位毫米/年。这里要特别注意,InVEST要的是多年平均,不是某一年的降水,别把单年数据填进去,否则结果只能解释那一年。

第二个是参考蒸散量栅格,单位同样是毫米/年。你可以用FAO Penman-Monteith公式结合气象站数据插值得到,也可以直接用全球或区域的蒸散产品。

第三个是土地利用/覆被栅格和它的属性表。属性表里必须有两个关键字段:植被蒸散系数Kc和根系深度root_depth。Kc一般在0到1.5之间,跟作物系数类似;root_depth是不同地类的根系深度,单位也是毫米。我通常通过连接字段把LULC属性表接好,确保每个地类唯一编码都有对应的Kc和root_depth。

第四个就是我们的主角:根系限制层深度栅格,直接选处理好的china_dtb_250m_mm.tif。

第五个是植物可利用水比例栅格,PAWC,取值范围0到1。这个一般由土壤质地数据换算得到,常用公式或者直接查表。

第六个是流域边界矢量,注意别跟栅格坐标系打架,统一投影后再导入。

最后是Z常数。Z值是季节降雨特征参数,一般在1到10之间,湿润地区偏低、干旱地区偏高。很多项目直接填默认值,但如果你研究区季节差异大,最好用实测径流来率定。我这套基岩深度数据不会直接改变Z的最优取值,但基岩深度影响AWC,而AWC进入ω公式后跟Z耦合,所以原则上Z也要跟着数据源一起重新校准,不能照搬别人写的区域参数表。

参数填完后,指定输出目录就可以运行了。模型运行时间一般在几分钟到几十分钟,看流域大小和像元数量,250米分辨率全国尺度通常十几分钟能跑完。

3.3 运行结果合理性校验

模型跑完别急着出图写报告,先做三个合理性检查。

一是看产水量的数值范围。湿润区年产水量在几百毫米合理,严重干旱区接近0也合理,但如果湿润区算出来上千米,或者干旱区算出来几百毫米,就要回头查输入栅格单位了。

二是看空间格局。把产水量图层跟DEM、土地利用图层做视觉对比。山区降水多、植被少的区域产水应该偏高,平原农田区蒸散大产水应该偏低。如果出现完全反着的格局,多半是土地利用的Kc配错了。

三是按流域统计总量,跟实测径流资料对一对。市县级水文站往往有多年平均径流深数据,用流域边界做zonal统计,算出的产水总量跟实测径流深对比,误差在20%以内就算很理想。如果差得远,优先调Kc和Z,而不是动基岩深度。基岩深度数据作为底层输入,除非确信数据本身有错,否则不建议为了拟合径流去修改它,那属于掩耳盗铃。

4. 常见问题与排查技巧实录

4.1 栅格范围、分辨率不一致导致的报错

InVEST跑产水模型最常见的问题,就是所有栅格放进面板后,运行日志抛出一堆“extent mismatch”或者“grid aligned”之类的报错。很多新手第一时间以为是模型坏了,实际上就是图层之间范围、像元对齐差了一丁点。

解决办法不是去面板里找开关,而是在进入模型之前统一好所有输入。我在ArcMap里的标准操作是:选中基岩深度栅格作为基准,用Environment设置里的Snap Raster锁定它,然后分别对其他栅格做重采样。这样所有输出栅格的像元位置完全对齐,范围也一致,InVEST那边就安静了。另外注意投影,如果你用了WGS84地理坐标系,InVEST处理起来会偏慢且可能出现奇怪的单位问题,尽量在预处理阶段转成米制投影。

4.2 基岩深度为0和异常值怎么处理

SoilGrids数据里,基岩深度为0意味着基岩直接出露,现实中常见于裸岩、戈壁、高海拔流石滩。但问题在于,0值在InVEST产水模型里会直接导致土壤水库容量为0,AWC为0,进而让该像元的实际蒸散被压到极低,产水量被推到接近降水量的极端数值。这在地貌上说得通,但如果你的研究区里零值区域占比很大,就要想想这个极端产水信号是否掩盖了周边区域的格局。

我的处理策略是,先统计零值像元的数量和空间分布。如果只是在陡峭山脊线零星出现,直接保留即可;如果是大面积戈壁,建议把基岩深度下限设为10到20厘米,避免数值奇点导致模型数值不稳定。还有一类异常是栅格边缘的环状低值,通常是裁剪范围没留缓冲区、重采样插值时在边界处拉出过渡带,把边缘100到200米内的像元检查一下,必要时用掩膜把国境线外的值清成NoData,再让InVEST自己忽略。

4.3 ArcMap里把栅格数据导出Excel的两种实用办法

做产水模型前,我经常需要把基岩深度栅格抽出来跟气象站点、土壤采样点的实测数据做回归分析,这就绕不开一个高频操作:把栅格数据转化导出为Excel。

第一种办法适合点位明确的情况。先在ArcMap里准备一个点要素层,比如你的采样点或水文站坐标点,然后打开ArcToolbox里的Spatial Analyst Tools,找到Extract Multi Values to Points。输入栅格可以是好几个,它会自动在每个点上提取所有栅格的像元值,写到点属性表里。提取完打开点图层的属性表,你会发现多了几列字段,值已经挂好。然后,关键一步:找到Conversion Tools下的Table to Excel工具,输出成xls或xlsx,Excel里直接就能用。这个流程是我最推荐的,因为一次能把降水、基岩深度、PAWC全部提取到一张表里,后续做相关性分析特别顺手。

第二种办法适合没有点文件的纯栅格抽样。比如你就想看几个像元对应的基岩深度值,可以先用Create Fishnet生成一个采样网格点,然后同样用Extract Values to Points提取。或者更简单一点,如果你的栅格是整型且像元数量不算太爆炸,直接右键图层打开属性表,选择按像元值统计,再把表导出。不过连续型浮点栅格一般不显示属性表,这时候必须先转成点或整型才能导出。实话说,80%的项目需求用第一种办法就够,不用纠结第二种。

还有个细节提示:Table to Excel工具在老版本ArcMap里不是默认加载,需要确保安装了Data Interoperability或相应扩展模块。如果提示找不到工具,你可以把点属性表导出成dbf,用Excel直接打开dbf文件也能看到数据,只是格式上不如xlsx方便。

4.4 避坑清单速查表

下面这张表是我这几套流程跑下来整理的避坑清单,每一条都对应一次真实翻车现场:

坑位点症状处理办法
单位忘换算产水量整体偏低/偏高近10倍基岩深度cm转mm,全体乘10后重导
投影混用栅格范围错位、结果变形全流程锁定同一米制投影
分辨率不一致InVEST报grid aligned错误以250米基岩深度为基准栅格,统一重采样
NoData空洞流域内部出现产水黑洞先填补空洞,再统一NoData为-9999
基岩深度大面积0值产水极端值掩盖区域格局对裸岩区设下限或单独掩膜
Kc和root_depth表漏字段模型找不到属性列直接报错跑模型前先Check Model Inputs
把单年降水当多年平均结果年份特异性强确认降水栅格为多年均值
边界裁剪无缓冲国境线边缘出现锯齿细缝预留半个像元以上的缓冲区

表格里的问题我几乎都踩过,尤其是单位和NoData这两个,属于最不起眼但破坏力最大的两类。

最后想多说一句

在InVEST产水模型这条技术线路上摸爬滚打这几年,我最大的体会是:真正决定模型上限的往往不是算法,而是数据预处理得干不干净。基岩深度栅格数据看起来只是流程里一个普普通通的输入项,但它在土壤储水、蒸散控制、产水格局这条链条上扮演了关键角色。把这套250米、源自ISRIC World Soil Information的基岩深度数据用好,你的产水模型结果就稳定了一大半。

我实际使用中还有一个习惯:把预处理后的基岩深度栅格跟其他关键图层一起打包保存,连同处理日志记清楚投影、重采样方法、NoData设置。这样过半年再看这个项目,或者换同事接手,都能快速回溯,不用重新踩一遍坑。下次如果你也在InVEST里被某个栅格参数磨得没脾气,不妨回头检查一下数据底子,多半会有惊喜。

需要专业的网站建设服务?

联系我们获取免费的网站建设咨询和方案报价,让我们帮助您实现业务目标

立即咨询