[发明专利]一种基于栅格DEM的山丘区土壤厚度预测的方法

权利要求书

1.一种基于栅格DEM的山丘区土壤厚度预测的方法，其特征在于：

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

$\frac{\∂h}{\∂t}=-\mathrm{\η}\frac{\∂e}{\∂t}-\▿\·q---\left(1\right)$

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

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

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

q＝q_d+q_t (3)

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

$q_d=-k_d\·\▿z---\left(4\right)$

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

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

$q_t=-k_t\·A^m\·{(\▿z)}^n---\left(5\right)$

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

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

$\frac{\∂h}{\∂t}={\mathrm{ce}}^{-h/h_0}+f---\left(6\right)$

其中：c＝ηP₀，且：

$f=-\▿q_d-\▿q_t---\left(7\right)$

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

通过3行×3列的栅格网络计算土壤输移通量，栅格i的C_i计算如下：

$C_i=\frac{2(z_1+z_3+z_5+z_7)+(z_2+z_4+z_6+z_8)-12z_i}{4{\mathrm{\Δx}}^2}---\left(8\right)$

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

$\begin{array}{c}-\▿q_t=\frac{1}{\▿x}(\underset{j}{\Σ}{\mathrm{\λ}}_{{\mathrm{in}}_j}{q}_{{\mathrm{tn}}_j}-\underset{k}{\Σ}{\mathrm{\λ}}_{{\mathrm{io}}_k}{q}_{{\mathrm{to}}_k})=\\ \frac{1}{\mathrm{\Δ}x}\underset{j}{\Σ}{\mathrm{\λ}}_{{\mathrm{in}}_j}{k}_t{A}_j^m{\left(\frac{z_j-z_i}{\mathrm{\Δ}x}\right)}^n-\frac{1}{\mathrm{\Δ}x}\underset{k}{\Σ}{\mathrm{\λ}}_{{\mathrm{io}}_k}{k}_t{A}_i^m{\left(\frac{z_j-z_i}{\mathrm{\Δ}x}\right)}^n\end{array}---\left(9\right)$

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

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

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

首先，方程改写为：

$\frac{1}{Z}(-c+h_0{Z}^{\′})=f---\left(11\right)$

其中：

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

Z′＝aZ+b (13)

其中

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

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

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

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

$Z=R_1{e}^{at}-\frac{b}{a}---\left(15\right)$

其中R₁是积分常量；

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

$e^{\frac{h}{h_0}}=R_1{e}^{at}-\frac{b}{a}$

$h\left(t\right)=h_{0}ln(R_1{e}^{at}-\frac{b}{a})---\left(16\right)$

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

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

$R_1=e^{h_i/h_0}+\frac{b}{a}---\left(18\right)$

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

$h\left(t\right)=h_{0}ln\left(\frac{({\mathrm{fe}}^{h_i/h_0}+c)e^{ft/h_0}-c}{f}\right)---\left(19\right)$

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

h(0)＝0 (20)

那么常数R₁为：

$R_1=1+\frac{b}{a}---\left(21\right)$

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

$h\left(t\right)=h_{0}ln\left(\frac{(f+c)e^{ft/h_0}-c}{f}\right)---\left(22\right)$

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

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

$\frac{\∂h}{\∂t}=f---\left(23\right)$

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

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

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