一种基于栅格DEM的山丘区土壤厚度预测的方法与流程

本发明属于水文模型技术领域，具体涉及一种基于栅格DEM的山丘区土壤厚度预测的方法。

背景技术：

由于土壤厚度在水文模型研究与应用中的重要地位，低成本高精度地获得土壤厚度制图显得非常重要，为了研究和应用的需要，各种各样的土壤厚度预测模型得到发展。大体而言，模型可以分为两大类，即随机模型和有物理基础的模型。

新兴发展起来的地理信息系统地形分析技术和广泛可获得的数字高程模型对随机模型起了辅助促进的作用。由于这些模型具有结构简单，参数较少的特点，多被用来建立地形因素和土壤属性(如土壤厚度)之间的定量联系。然而，所有的这些随机模型的应用都需遵从一个基本的假设，就是采样地的土壤厚度和地形因素之间的经验关系可以被拓展去推测其他地区的土壤厚度属性，即采样地与目标预测区具有相似性。为了使得模型满足这一假定，就需要大量且广泛的田间采样数据用于参数的校验，这无形中增加了模型的应用成本。正是因为统计模型的经验性，它使得这类模型仅适用于建模区，这已严重的限制了模型预测其他无资料或者少资料流域土壤厚度的能力。

相对而言，有物理基础的模型(如地貌演化模型)，则侧重于描述地质年代尺度上的土壤生成、输移的过程，即土壤厚度演化的过程。例如，通过假设土壤输移的速率与坡度成正比，即假定土壤运动为简单的蠕动模型，Dietrich et al.[1995]建立了基于栅格数字高程模型的预测流域尺度土壤厚度的数值模型。Roering et al.[2008]则进一步发展了Dietrich et al.[1995]的工作，他在已有的土壤厚度演化模型中引入三个线性、非线性土壤输移公式。但是，传统上这些土壤厚度演化模型通常被认为是地貌演化动力学模型的一部分，它的求解往往依赖于复杂且精确的求解技术，例如采用数值方法模拟地质时间尺度上气候和构造作用力影响下的地貌演化[Saco et al.，2006；Pelletier et al.，2011]。

Dietrich,W.E.,R.Reiss,M.-L.Hsu,and D.R.Montgomery(1995),A process-based model for colluvial soil depth and shallow landsliding using digital elevation data,Hydrol.Process.,9,383–400,doi:10.1002/hyp.3360090311.

Roering,J.J.(2008),How well can hillslope evolution models“explain”topography？,Geol.Soc.Am.Bull.,120,1248-1262,doi:10.1130/B26283.1.

Saco,P.M.,G.R.Willgoose,and G.R.Hancock(2006),Spatial organization of soil depths using a landform evolution model,J.Geophys.Res.,111,F02016,doi:10.1029/2005JF000351.

显然，在水文模拟中应用地貌演化数学模型非常专业、代价非常之大。因此，为了便于水文领域的应用，我们急需一些简单而通用的土壤厚度表达式，这些表达式一般为地形的函数。基于这一认识，许多研究往往采用解析或者降低数值求解难度的简化方法。例如，土壤生成和侵蚀的平衡态假设，Bertoldi et al.[2006]给出一个简单的土壤厚度预测的解析公式，其采用土壤蠕动模型来描述土壤的输移。Pelletier和Rasmussen[2009]则延续了Bertoldi的平衡稳定方法，其分别用三个非线性土壤输移方程来描述土壤颗粒运动。然而平衡稳定状态只在满足特定条件下的地区有效，并非广泛适用的。

Bertoldi,G.,R.Rigon,and T.M.Over(2006),Impact of watershed geomorphic characteristics on the energy and water budgets,J.Hydrometeorol.,7,389–403.

Pelletier,J.D.,and C.Rasmussen(2009),Geomorphically based predictive mapping of soil thickness in upland watersheds,Water Resour.Res.,45,W09417,doi:10.1029/2008WR007319.

土壤厚度是山坡水文过程的一个关键性的控制因素。早期研究表明，不同的土壤厚度空间分布模式在很大程度上影响降雨径流的比率，即径流系数。这是由于土壤厚度和土壤中的孔隙决定了山坡土壤的蓄水能力。因此，在很多有物理基础的水文模型中，土壤厚度都是一个很重要的变量。在概念性的水文模型中，土壤厚度也可用来导出或指示其中的重要参数，如土壤蓄水能力等。因此，水文研究和实践日益需要越来越高质量的土壤厚度制图。然而，到目前为止，水文学家在进行水文模拟和相关规划研究时很大程度上仍需依赖以往的土壤调查数据库，如U.S.Department of Agriculture’s national soils databases和the Soil Information System of China，来获取土壤厚度图。但是，标准的土壤调查不能提供高分辨率的土壤厚度图，这显然限制了分布式水文模型的有效运用。

技术实现要素：

土壤厚度是山坡水文过程的一个关键性的控制因素。早期研究表明，不同的土壤厚度空间分布模式在很大程度上影响降雨径流的比率，即径流系数。这是由于土壤厚度和土壤中的孔隙决定了山坡土壤的蓄水能力，因此，在很多有物理基础的水文模型中，土壤厚度都是一个很重要的变量。在概念性的水文模型中，土壤厚度也可用来导出或指示其中的重要参数，如土壤蓄水能力等。故水文研究和实践日益需要越来越高质量的土壤厚度制图。然而，到目前为止，水文学家在进行水文模拟和相关规划研究时很大程度上仍需依赖以往的土壤调查数据库，如U.S.Department of Agriculture’s national soils databases和the Soil Information System of China，来获取土壤厚度图。但是，标准的土壤调查不能提供高分辨率的土壤厚度图，这显然限制了分布式水文模型的有效运用。

因此，本发明提供了一种基于栅格DEM的山丘区土壤厚度预测的方法。该方法采用一个简单的基于非稳态假设的土壤厚度预测模型，即一个依据土壤演化动力学方程推导出的解析表达式；该模型运行在栅格数字高程模型之上，采用简化的解析方程预测土壤厚度时空演化及分布，最后通过数据可视化的技术手段生成可供广泛应用的土壤厚度制图。

为实现上述目的，本发明采用了以下技术方案：

一种基于栅格DEM的山丘区土壤厚度预测的方法，包括：

S1、根据土壤生成和输移的质量守恒，描述土壤演化的动力学方程写为：

h是土壤厚度，e是高程，q是土壤输移通量，η是岩石密度和土壤密度之比，即η＝ρ_r/ρ_s，方程(1)右边的第一项是由于风化作用，基岩上土壤生成的速率，

其中，P₀是裸露基岩上土壤生成率(即土壤厚度为零时的土壤生成率)，h₀是一个特征侵蚀深度，为经验参数；

S2、使用两个线性沉积物输移法则，即一个土壤蠕动输移模型q_d和一个径流输移模型q_t：

q＝q_d+q_t (3)

土壤蠕动q_d用一个线性蠕变函数表示：

其中k_d是扩散系数，z是土壤表面高程；

降雨径流驱动下的土壤输移一般方程写为：

k_t是沉积物输移系数，A是集水面积，m和n是一个给定地形的指数常数；

S3、将方程(2)，(4)和(5)代入方程(1)，得到用来描述土壤厚度演化的土壤输移动力学方程：

其中：c＝ηP₀，且：

其中f是土壤蠕动和降雨径流土运的下坡土壤输移通量散度，利用基于数字高程模型的数字水系，根据方程(6)和(7)计算f；根据方程(4)，方程(7)中的可表示为其中的值等于地形曲率，C；

通过三乘三的栅格网络计算土壤输移通量，栅格i的C_i计算如下：

其中z_i是研究栅格i的高程，z₁-z₈是栅格i周围八个栅格的高程，Δx系栅格的边长大小；方程(7)中的可用一个坡度和上游集水面积的函数来表达：

和q_tok分别是栅格i的土壤入流和出流通量，是比例系数；

方程(6)和(7)通过数值方法求解，采用显式差分法来求解微分方程(6)和(7)；

S4、提出两个假设：首先，假定研究的对象仅限于湿润和半湿润的地区，这里壤中流较为丰富，是主要的径流成分，且下垫面基岩机械性质足够强；其次，假设在相对较短的地质年代尺度内，没有构造隆升或者基岩下沉的情况，即地表地形在土壤厚度演化过程中相对稳定或稳态变化。因此，可以推定地形特征(如曲率)的演化受岩石风化和土壤输移作用的驱动，且认为其速率比土壤厚度演化的速率慢得多。这样，方程(6)可被看作一个一阶非线性、非齐次、常微分方程，方程(6)的解析解可以通过如下导出：

首先，方程改写为：

其中：

之后，改写成一个线性常微分方程：

Z′＝aZ+b (13)

其中

方程(13)的两边都乘以e^-at并重组为：

Z′e^-at-aZe^-at＝(Ze^-at)′＝be^-at

得到：(Ze^-at)′＝be^-at (14)

对方程(14)进行积分，得到线性常微分方程的一般解：

其中R₁是积分常量；

方程(12)代入方程(15)中，土壤厚度的一般表达就可导出为：

在方程(16)中，系数R₁可以通过设定初始土壤厚度来决定；所以，如果初始土壤厚度假设为：h(0)＝hi (17)

那么积分常量R₁导出为：

将方程(18)代入(16)，得到：

方程(19)是土壤厚度演化的普通形式，土壤厚度演化模型对初始土壤厚度的敏感度不高，所以，初始土壤厚度假定为：

h(0)＝0 (20)

那么常数R₁为：

从而，方程(16)可改写成一个更简化的表达，如下：

其中h₀,c是土壤生成参数，f是整体土壤侵蚀率，可以根据方程(7)决定；

由于土壤厚度一般在流域的低洼地增长，根据方程(2)土壤生成率为零，那么方程(6)重写成：

对方程(23)积分，积分时间从t_s到t，土壤厚度从h_s到h，得到一个可以预测土壤厚度的方程：

h(t)＝h_s+f(t-t_s) (24)

其中t_s是土壤生成率为零的时间，h_s是时刻t_s的土壤厚度。

优选的，土壤生成方程即方程(2)中的参数P₀根据放射性同位素测量。

优选的，方程(4)中的扩散系数k_d范围值是10^-5到10^-2m²/yr(平方米每年)。

优选的，基于Fagherazzi的定义，m和n的值分别是1.4和1.2。

优选的，方程(9)中的是D∞method确定的比例系数。

进一步的，对于出流的系数和可定义为：

其中α₁,α₂分别是与基方向以及斜方向的夹角。

本发明的有益效果在于：

1)本发明开发了一个简单的基于非稳态假设的土壤厚度预测模型，该模型是依据土壤演化动力学方程推导出的解析表达式，像Dietrich et al.[1995]和Pelletier及Rasmussen[2009]的模型一样，该模型也是运行在栅格数字高程模型之上，但是本发明采用简化的解析方程预测土壤厚度时空演化及分布，最后通过数据可视化的技术手段生成可供广泛应用的土壤厚度制图，具有模型简单、广泛适用、易于实践、提供高质量土壤厚度制图的优势。

2)土壤厚度在上游水文过程中扮演着重要的控制作用。然而，在流域水文研究中，人们对其土壤厚度的空间分布仍然知之甚少，本发明基于构造稳定的假设，即在较近的地质年代，流域未发生构造隆升或下沉，导出一个基于地貌演化动力方程的土壤厚度预测解析解能帮助水文学家理解土壤厚度随地质年代变化的规律，也为相应的水文模拟研究提供了一个土壤厚度空间分布预测的实用工具。

附图说明

图1为三乘三表格的3D结构图。

图2为计算曲率的栅格排序图。

图3为土壤通量f导出的计算网格。其中，q代表土壤输移通量，灰色箭头代表栅格i的三角刻面方向α₁,α₂分别是与基方向以及斜方向的夹角。

图1、2、3即为本发明三乘三表格中土壤输移通量演算图。

具体实施方式

下面将结合本发明实施例中的附图，对本发明中的技术方案进行清楚、完整地描述。以下实施例仅用于更加清楚地说明本发明的技术方案，而不能以此来限制本发明的保护范围。

一种基于栅格DEM的山丘区土壤厚度预测的方法

S1、根据土壤生成和输移的质量守恒，描述土壤演化的动力学方程写为：

h是土壤厚度，t为时间，e是高程，q是土壤输移通量，η是岩石密度和土壤密度之比，即η＝ρ_r/ρ_s，方程(1)右边的第一项是由于风化作用，基岩上土壤生成的速率，Heimsath et al.[2000]研究表明土壤生成率随土壤厚度的增加大致呈指数下降：

P₀是裸露基岩上土壤生成率(即土壤厚度为零时的土壤生成率)，h₀是一个特征侵蚀深度，为经验参数；

Heimsath,A.M.,J.Chappell,W.E.Dietrich,K.Nishiizumi,and R.C.Finkel(2000),Soil production on a retreating escarpment in southeastern Australia,Geology,28,787-790,doi:10.1130/0091-7613.

S2、根据McKean et al.[1993]，Small et al.[1999]，Dietrich et al[1995]等人的报告，在许多湿润，半湿润流域，坡度是土壤输移主要的主导因素，即土壤蠕动主导土壤输移的模型，此外，这些线性输移法则会得出质量守恒方程的一个简单解析解。所以，在本发明中，使用两个线性沉积物输移法则，即一个土壤蠕动输移模型q_d和一个径流输移模型q_t：

q＝q_d+q_t (3)

土壤蠕动q_d用一个线性蠕变函数表示：

其中k_d是扩散系数，z是土壤表面高程；

根据Willgoose et al.[1991b]and Saco et al.[2006]，降雨径流驱动下的土壤输移一般方程写为：

k_t是沉积物输移系数，A是集水面积，m和n是一个给定地形的指数常数。基于Fagherazzi的定义，m和n的值分别是1.4和1.2。

McKean,J.A.,W.E.Dietrich,R.C.Finkel,J.R.Southon,and M.W.Caffee(1993),Quantification of soil production and downslope creep rates from cosmogenic 10Be accumulations on a hillslope profile,Geology,21,343–346.

Small,E.E.,R.S.Anderson,G.S.Hancock,and R.C.Finkel(1999),Estimates of regolith production from 10Be and 26Al:Evidence for steady state alpine hillslopes,Geomorphology,27,131–150.

Willgoose,G.,R.L.Bras,and I.Rodriguez-Iturbe(1991b),A physical explanation of an observed link area-slope relationship,Water Resour.Res.,27,1697-1702.

Fagherazzi,S.,A.D.Howard,and P.L.Wiberg(2002),Am implicit finite difference method for drainage basin evolution,Water Resour.Res.,38(7),1116,doi:10.1029/2001WR000721.

S3、将方程(2)，(4)和(5)代入方程(1)，得到用来描述土壤厚度演化的土壤输移动力学方程：

其中：c＝ηP₀，且：

其中f是土壤蠕动和降雨径流土运的下坡土壤通量散度，为了根据方程(6)和(7)计算f，要利用基于数字高程模型(DEM)的数字水系。如图1所述，演示了怎样通过栅格网络计算土壤输移通量。根据方程(4)，方程(7)中的可表示为其中的值等于地形曲率，C；

图2中栅格i的C_i根据Heimsath et al.[1999]可计算如下：

其中z_i是研究栅格i的高程，z₁-z₈是图2中栅格i周围八个栅格的高程，Δx系栅格的边长大小；方程(7)中的可用一个坡度和上游集水面积的函数来表达：

和q_tok分别是栅格i的土壤入流和出流通量，是D∞method[Tarboton,1997]确定的比例系数，例如对于出流的系数和可定义为：

其中α₁,α₂是分别是图3中与基方向以及斜方向的夹角；

方程(6)和(7)可通过数值方法求解，采用MacCormack scheme的显式差分法来求解微分方程(6)和(7)。

Heimsath,A.M.,W.E.Dietrich,K.Nishiizumi,and R.C.Finkel(1999),Cosmogenic nuclides,topography,and the spatial variation of soil depth.Geomorphology,27,151-172.

Tarboton,D.G.(1997),A new method for the determination of flow directions and upslope areas in grid digital elevation models,Water Resour.Res.,33,309-319.D∞method引用该文献中的方法。

MacCormack R.W.(1971),Numerical solution of the interaction of a shock wave with a laminar boundary layer.Lecture Notes in Physics 8:151–163,Springer-Verlag.

S4、为了进一步简化动态方程提出两个假设：首先，假定研究的对象仅限于湿润和半湿润的地区，这里壤中流较为丰富，是主要的径流成分，且下垫面基岩机械性质足够强；其次，假设在相对较短的地质年代尺度内，没有构造隆升或者基岩下沉的情况，即地表地形在土壤厚度演化过程中相对稳定或稳态变化。因此，可以推定地形特征(如曲率)的演化受岩石风化和土壤输移作用的驱动，且认为其速率比土壤厚度演化的速率慢得多。这样，方程(6)可以看作是一个一阶非线性非齐次常微分方程。方程(6)的一般解可以通过以下的步骤导出。

首先，方程改写为：

其中：

之后，改写成一个线性常微分方程：

Z′＝aZ+b (13)

其中

方程(13)的两边都乘以e^-at并重组为：

Z′e^-at-aZe^-at＝(Ze^-at)′＝be^-at

得到：(Ze^-at)′＝be^-at (14)

对方程(14)进行积分，得到线性常微分方程的一般解：

其中R₁是积分常量；

方程(12)代入方程(15)中，土壤厚度的一般表达就可导出为：

在方程(16)中，系数R₁可以通过设定初始土壤厚度来决定；所以，如果初始土壤厚度假设为：h(0)＝hi (17)

那么积分常量R₁导出为：

将方程(18)代入(16)，得到：

方程(19)是土壤厚度演化的普通形式，根据Dietrich et al.[1995]提及，土壤厚度演化模型对初始土壤厚度的敏感度不高，所以，初始土壤厚度假定为：

h(0)＝0 (20)

那么常数R₁为：

从而，方程(16)可改写成一个更简化的表达，如下：

其中h₀,c是土壤生成参数(可以根据他们物理意义定义它的值)，f是整体土壤侵蚀率，可以根据方程(7)决定；

由于土壤厚度一般在流域的低洼地增长，根据方程(2)土壤生成率为零，那么方程(6)重写成：

对方程(23)积分，积分时间从t_s到t，土壤厚度从h_s到h，得到一个可以预测土壤厚度的方程：

h(t)＝h_s+f(t-t_s) (24)

其中t_s是土壤生成率为零的时间，h_s是时刻t_s的土壤厚度。

在此分析模型中，土壤生成方程即方程(2)中的参数P₀根据放射性同位素测量。

尽管方程(4)和(5)中的扩散系数k_d在土壤输移公式中带有经验的特性，但也可以理论地根据现场实验来确定，如果河道河床迁移速率可以评估的话。根据Chiang和Hsu[2006]已公开的研究成果，方程(4)中的参数扩散系数k_d范围值是10^-5到10^-2m²/yr(平方米每年)。

Chiang,S.H.,M.L.Hsu(2006),Parameter calibration in a process-based soil depth estimation model assuming local steady state,J.Geograph.Sci.,44,23-38[in Chinese with English Abstract].

以上所述仅是本发明的优选实施方式，应当指出，对于本技术领域的普通技术人员来说，在不脱离本发明技术原理的前提下，还可以做出若干改进和变形，这些改进和变形也应视为本发明的保护范围。