水流强度指数（SPI）地形湿度指数（TWI）计算

指数计算公式：

        本文介绍如何在ArcGIS中利用DEM数据进行水流强度指数（stream power index，SPI）和地形湿度指数（topographic wetness index，TWI）的计算。

DEM数据可在【数据分享】系列推文中免费获取。

TWI是区域地形对径流流向和蓄积影响的物理指标。公式如下：

TWI=ln[SCA / Tan(slope)]

其中，SCA代表单位面积的汇流量，slope代表坡度（以度为单位）。

---

SPI值越高，表明径流集中可能导致土壤侵蚀。SPI公式形式较多，常见公式如下：

SPI=SCA * Tan(slope)

SPI=ln[SCA * Tan(slope)]

SPI=ln[SCA * Tan(slope)*100]

其中，SCA和slope意义与上述相同。本文基于公式2进行计算。

由公式可知，SPI和TWI计算均基于SCA和slope，因此本文先介绍如何获得这两个数据。再依据公式进行SPI和TWI的计算。

---

一、SCA和slope计算

1、对DEM数据进行填洼

ArcToolbox→Spatial Analyst Tools (Spatial Analyst工具)→Hydrology (水文分析)→Fill (填洼)

2、填洼后的数据进行坡度计算，获得slope

ArcToolbox→Spatial Analyst Tools (Spatial Analyst工具)→Surface (表面分析)→Slope(坡度)

3、填洼后的数据计算水流方向

ArcToolbox→Spatial Analyst Tools (Spatial Analyst工具)→Hydrology (水文分析)→Flow Direction(流向)

4、水流方向数据计算汇流累积量

ArcToolbox→Spatial Analyst Tools (Spatial Analyst工具)→Hydrology (水文分析)→Flow Accumulation(流量)

5、汇流累积量数据计算单位面积的汇流量，获得SCA

ArcToolbox→Spatial Analyst Tools (Spatial Analyst工具)→Map Alqebra (地图代数)→Raster Calculator (栅格计算器)

公式为：

Con("Flow_Acc.tif"==0,1,"Flow_Acc.tif")*900/Con("Flow_Dir.tif"==1,30,Con("Flow_Dir.tif"==4,30,Con("Flow_Dir.tif"==16,30,Con("Flow_Dir.tif"==64,30,Con("Flow_Dir.tif"==2,30*SquareRoot(2),Con("Flow_Dir.tif"==8,30*SquareRoot(2),Con("Flow_Dir.tif"==32,30*SquareRoot(2),Con("Flow_Dir.tif"==128,30*SquareRoot(2))))))))) 

其中，900为栅格大小，即栅格分辨率的平方，需要根据自己的数据像元大小进行修改（本文为30m*30m,因此为900），且当"Flow Accumulation"为0时没意义，至少应为1。对于流向为1、4、16、64 的栅格，除以栅格尺寸，这里是30；对于流向为2、8、32、128的栅格，除以栅格尺寸的√2倍，这里是√2倍的30。

 二、SPI计算

根据公式计算SPI   

SPI=ln[SCA * Tan(slope)]

ArcToolbox→Spatial Analyst Tools (Spatial Analyst工具)→Map Alqebra (地图代数)→Raster Calculator (栅格计算器)

公式为：

Ln("SCA.tif"*Tan(Con("Fill_slope.tif<=0,0.00001,Con("Fill_slope.tif">0,"Fill_slope.tif"*3.1415926/180))))

其中，求Tan（slope）时需要先将角度化为弧度；且当角度为0时将其自定义为0.00001，否则无意义。

 三、TWI计算

根据公式计算TWI   

TWI=ln[SCA / Tan(slope)]

ArcToolbox→Spatial Analyst Tools (Spatial Analyst工具)→Map Alqebra (地图代数)→Raster Calculator (栅格计算器)

公式为：

Ln("SCA.tif"/Tan(Con("Fill_slope.tif<=0,0.00001,Con("Fill_slope.tif">0,"Fill_slope.tif"*3.1415926/180))))

同理，求Tan（slope）时需要先将角度化为弧度；且当角度为0时将其自定义为0.00001，否则无意义。

---

参考文献

【1】Gessler, P.E., Peterson, G.A., 1993. Soil Attribute Prediction Using Terrain Analysis. Soil Sci. Soc. Am. J. 57. https://doi.org/10.2136/sssaj1993.572npb

【2】Mhiret, D.A., Dagnew, D.C., Assefa, T.T., Tilahun, S.A., Zaitchik, B.F., Steenhuis, T.S., 2019. Erosion hotspot identification in the sub-humid Ethiopian highlands. Ecohydrol. Hydrobiol. 19, 146–154. https://doi.org/10.1016/j.ecohyd.2018.08.004

【3】Mohseni, N., Salar, Y.S., 2021. Terrain indices control the quality of soil total carbon stock within water erosion-prone environments. Ecohydrol. Hydrobiol. 21, 46–54.https://doi.org/10.1016/j.ecohyd.2020.08.006

【4】刘欣. 基于地形指标的黄土高原坡面侵蚀沟发生与空间分异研究[D].西北大学,2021.DOI:10.27405/d.cnki.gxbdu.2021.001838.

---

微信公众号【GIS攻略】同步更新，其中高程数据可在公众号【数据分享】中获得，供大家练习。