本發(fā)明屬于水文模型技術(shù)領(lǐng)域，具體涉及一種基于柵格DEM的山丘區(qū)土壤厚度預(yù)測的方法。

背景技術(shù)：

由于土壤厚度在水文模型研究與應(yīng)用中的重要地位，低成本高精度地獲得土壤厚度制圖顯得非常重要，為了研究和應(yīng)用的需要，各種各樣的土壤厚度預(yù)測模型得到發(fā)展。大體而言，模型可以分為兩大類，即隨機(jī)模型和有物理基礎(chǔ)的模型。

新興發(fā)展起來的地理信息系統(tǒng)地形分析技術(shù)和廣泛可獲得的數(shù)字高程模型對隨機(jī)模型起了輔助促進(jìn)的作用。由于這些模型具有結(jié)構(gòu)簡單，參數(shù)較少的特點(diǎn)，多被用來建立地形因素和土壤屬性(如土壤厚度)之間的定量聯(lián)系。然而，所有的這些隨機(jī)模型的應(yīng)用都需遵從一個(gè)基本的假設(shè)，就是采樣地的土壤厚度和地形因素之間的經(jīng)驗(yàn)關(guān)系可以被拓展去推測其他地區(qū)的土壤厚度屬性，即采樣地與目標(biāo)預(yù)測區(qū)具有相似性。為了使得模型滿足這一假定，就需要大量且廣泛的田間采樣數(shù)據(jù)用于參數(shù)的校驗(yàn)，這無形中增加了模型的應(yīng)用成本。正是因?yàn)榻y(tǒng)計(jì)模型的經(jīng)驗(yàn)性，它使得這類模型僅適用于建模區(qū)，這已嚴(yán)重的限制了模型預(yù)測其他無資料或者少資料流域土壤厚度的能力。

相對而言，有物理基礎(chǔ)的模型(如地貌演化模型)，則側(cè)重于描述地質(zhì)年代尺度上的土壤生成、輸移的過程，即土壤厚度演化的過程。例如，通過假設(shè)土壤輸移的速率與坡度成正比，即假定土壤運(yùn)動(dòng)為簡單的蠕動(dòng)模型，Dietrich et al.[1995]建立了基于柵格數(shù)字高程模型的預(yù)測流域尺度土壤厚度的數(shù)值模型。Roering et al.[2008]則進(jìn)一步發(fā)展了Dietrich et al.[1995]的工作，他在已有的土壤厚度演化模型中引入三個(gè)線性、非線性土壤輸移公式。但是，傳統(tǒng)上這些土壤厚度演化模型通常被認(rèn)為是地貌演化動(dòng)力學(xué)模型的一部分，它的求解往往依賴于復(fù)雜且精確的求解技術(shù)，例如采用數(shù)值方法模擬地質(zhì)時(shí)間尺度上氣候和構(gòu)造作用力影響下的地貌演化[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.

顯然，在水文模擬中應(yīng)用地貌演化數(shù)學(xué)模型非常專業(yè)、代價(jià)非常之大。因此，為了便于水文領(lǐng)域的應(yīng)用，我們急需一些簡單而通用的土壤厚度表達(dá)式，這些表達(dá)式一般為地形的函數(shù)?；谶@一認(rèn)識，許多研究往往采用解析或者降低數(shù)值求解難度的簡化方法。例如，土壤生成和侵蝕的平衡態(tài)假設(shè)，Bertoldi et al.[2006]給出一個(gè)簡單的土壤厚度預(yù)測的解析公式，其采用土壤蠕動(dòng)模型來描述土壤的輸移。Pelletier和Rasmussen[2009]則延續(xù)了Bertoldi的平衡穩(wěn)定方法，其分別用三個(gè)非線性土壤輸移方程來描述土壤顆粒運(yùn)動(dòng)。然而平衡穩(wěn)定狀態(tài)只在滿足特定條件下的地區(qū)有效，并非廣泛適用的。

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.

土壤厚度是山坡水文過程的一個(gè)關(guān)鍵性的控制因素。早期研究表明，不同的土壤厚度空間分布模式在很大程度上影響降雨徑流的比率，即徑流系數(shù)。這是由于土壤厚度和土壤中的孔隙決定了山坡土壤的蓄水能力。因此，在很多有物理基礎(chǔ)的水文模型中，土壤厚度都是一個(gè)很重要的變量。在概念性的水文模型中，土壤厚度也可用來導(dǎo)出或指示其中的重要參數(shù)，如土壤蓄水能力等。因此，水文研究和實(shí)踐日益需要越來越高質(zhì)量的土壤厚度制圖。然而，到目前為止，水文學(xué)家在進(jìn)行水文模擬和相關(guān)規(guī)劃研究時(shí)很大程度上仍需依賴以往的土壤調(diào)查數(shù)據(jù)庫，如U.S.Department of Agriculture’s national soils databases和the Soil Information System of China，來獲取土壤厚度圖。但是，標(biāo)準(zhǔn)的土壤調(diào)查不能提供高分辨率的土壤厚度圖，這顯然限制了分布式水文模型的有效運(yùn)用。

技術(shù)實(shí)現(xiàn)要素：

土壤厚度是山坡水文過程的一個(gè)關(guān)鍵性的控制因素。早期研究表明，不同的土壤厚度空間分布模式在很大程度上影響降雨徑流的比率，即徑流系數(shù)。這是由于土壤厚度和土壤中的孔隙決定了山坡土壤的蓄水能力，因此，在很多有物理基礎(chǔ)的水文模型中，土壤厚度都是一個(gè)很重要的變量。在概念性的水文模型中，土壤厚度也可用來導(dǎo)出或指示其中的重要參數(shù)，如土壤蓄水能力等。故水文研究和實(shí)踐日益需要越來越高質(zhì)量的土壤厚度制圖。然而，到目前為止，水文學(xué)家在進(jìn)行水文模擬和相關(guān)規(guī)劃研究時(shí)很大程度上仍需依賴以往的土壤調(diào)查數(shù)據(jù)庫，如U.S.Department of Agriculture’s national soils databases和the Soil Information System of China，來獲取土壤厚度圖。但是，標(biāo)準(zhǔn)的土壤調(diào)查不能提供高分辨率的土壤厚度圖，這顯然限制了分布式水文模型的有效運(yùn)用。

因此，本發(fā)明提供了一種基于柵格DEM的山丘區(qū)土壤厚度預(yù)測的方法。該方法采用一個(gè)簡單的基于非穩(wěn)態(tài)假設(shè)的土壤厚度預(yù)測模型，即一個(gè)依據(jù)土壤演化動(dòng)力學(xué)方程推導(dǎo)出的解析表達(dá)式；該模型運(yùn)行在柵格數(shù)字高程模型之上，采用簡化的解析方程預(yù)測土壤厚度時(shí)空演化及分布，最后通過數(shù)據(jù)可視化的技術(shù)手段生成可供廣泛應(yīng)用的土壤厚度制圖。

為實(shí)現(xiàn)上述目的，本發(fā)明采用了以下技術(shù)方案：

一種基于柵格DEM的山丘區(qū)土壤厚度預(yù)測的方法，包括：

S1、根據(jù)土壤生成和輸移的質(zhì)量守恒，描述土壤演化的動(dòng)力學(xué)方程寫為：

h是土壤厚度，e是高程，q是土壤輸移通量，η是巖石密度和土壤密度之比，即η＝ρ_r/ρ_s，方程(1)右邊的第一項(xiàng)是由于風(fēng)化作用，基巖上土壤生成的速率，

其中，P₀是裸露基巖上土壤生成率(即土壤厚度為零時(shí)的土壤生成率)，h₀是一個(gè)特征侵蝕深度，為經(jīng)驗(yàn)參數(shù)；

S2、使用兩個(gè)線性沉積物輸移法則，即一個(gè)土壤蠕動(dòng)輸移模型q_d和一個(gè)徑流輸移模型q_t：

q＝q_d+q_t (3)

土壤蠕動(dòng)q_d用一個(gè)線性蠕變函數(shù)表示：

其中k_d是擴(kuò)散系數(shù)，z是土壤表面高程；

降雨徑流驅(qū)動(dòng)下的土壤輸移一般方程寫為：

k_t是沉積物輸移系數(shù)，A是集水面積，m和n是一個(gè)給定地形的指數(shù)常數(shù)；

S3、將方程(2)，(4)和(5)代入方程(1)，得到用來描述土壤厚度演化的土壤輸移動(dòng)力學(xué)方程：

其中：c＝ηP₀，且：

其中f是土壤蠕動(dòng)和降雨徑流土運(yùn)的下坡土壤輸移通量散度，利用基于數(shù)字高程模型的數(shù)字水系，根據(jù)方程(6)和(7)計(jì)算f；根據(jù)方程(4)，方程(7)中的可表示為其中的值等于地形曲率，C；

通過三乘三的柵格網(wǎng)絡(luò)計(jì)算土壤輸移通量，柵格i的C_i計(jì)算如下：

其中z_i是研究柵格i的高程，z₁-z₈是柵格i周圍八個(gè)柵格的高程，Δx系柵格的邊長大??；方程(7)中的可用一個(gè)坡度和上游集水面積的函數(shù)來表達(dá)：

和q_tok分別是柵格i的土壤入流和出流通量，是比例系數(shù)；

方程(6)和(7)通過數(shù)值方法求解，采用顯式差分法來求解微分方程(6)和(7)；

S4、提出兩個(gè)假設(shè)：首先，假定研究的對象僅限于濕潤和半濕潤的地區(qū)，這里壤中流較為豐富，是主要的徑流成分，且下墊面基巖機(jī)械性質(zhì)足夠強(qiáng)；其次，假設(shè)在相對較短的地質(zhì)年代尺度內(nèi)，沒有構(gòu)造隆升或者基巖下沉的情況，即地表地形在土壤厚度演化過程中相對穩(wěn)定或穩(wěn)態(tài)變化。因此，可以推定地形特征(如曲率)的演化受巖石風(fēng)化和土壤輸移作用的驅(qū)動(dòng)，且認(rèn)為其速率比土壤厚度演化的速率慢得多。這樣，方程(6)可被看作一個(gè)一階非線性、非齊次、常微分方程，方程(6)的解析解可以通過如下導(dǎo)出：

首先，方程改寫為：

其中：

之后，改寫成一個(gè)線性常微分方程：

Z′＝aZ+b (13)

其中

方程(13)的兩邊都乘以e^-at并重組為：

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

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

對方程(14)進(jìn)行積分，得到線性常微分方程的一般解：

其中R₁是積分常量；

方程(12)代入方程(15)中，土壤厚度的一般表達(dá)就可導(dǎo)出為：

在方程(16)中，系數(shù)R₁可以通過設(shè)定初始土壤厚度來決定；所以，如果初始土壤厚度假設(shè)為：h(0)＝hi (17)

那么積分常量R₁導(dǎo)出為：

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

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

h(0)＝0 (20)

那么常數(shù)R₁為：

從而，方程(16)可改寫成一個(gè)更簡化的表達(dá)，如下：

其中h₀,c是土壤生成參數(shù)，f是整體土壤侵蝕率，可以根據(jù)方程(7)決定；

由于土壤厚度一般在流域的低洼地增長，根據(jù)方程(2)土壤生成率為零，那么方程(6)重寫成：

對方程(23)積分，積分時(shí)間從t_s到t，土壤厚度從h_s到h，得到一個(gè)可以預(yù)測土壤厚度的方程：

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

其中t_s是土壤生成率為零的時(shí)間，h_s是時(shí)刻t_s的土壤厚度。

優(yōu)選的，土壤生成方程即方程(2)中的參數(shù)P₀根據(jù)放射性同位素測量。

優(yōu)選的，方程(4)中的擴(kuò)散系數(shù)k_d范圍值是10^-5到10^-2m²/yr(平方米每年)。

優(yōu)選的，基于Fagherazzi的定義，m和n的值分別是1.4和1.2。

優(yōu)選的，方程(9)中的是D∞method確定的比例系數(shù)。

進(jìn)一步的，對于出流的系數(shù)和可定義為：

其中α₁,α₂分別是與基方向以及斜方向的夾角。

本發(fā)明的有益效果在于：

1)本發(fā)明開發(fā)了一個(gè)簡單的基于非穩(wěn)態(tài)假設(shè)的土壤厚度預(yù)測模型，該模型是依據(jù)土壤演化動(dòng)力學(xué)方程推導(dǎo)出的解析表達(dá)式，像Dietrich et al.[1995]和Pelletier及Rasmussen[2009]的模型一樣，該模型也是運(yùn)行在柵格數(shù)字高程模型之上，但是本發(fā)明采用簡化的解析方程預(yù)測土壤厚度時(shí)空演化及分布，最后通過數(shù)據(jù)可視化的技術(shù)手段生成可供廣泛應(yīng)用的土壤厚度制圖，具有模型簡單、廣泛適用、易于實(shí)踐、提供高質(zhì)量土壤厚度制圖的優(yōu)勢。

2)土壤厚度在上游水文過程中扮演著重要的控制作用。然而，在流域水文研究中，人們對其土壤厚度的空間分布仍然知之甚少，本發(fā)明基于構(gòu)造穩(wěn)定的假設(shè)，即在較近的地質(zhì)年代，流域未發(fā)生構(gòu)造隆升或下沉，導(dǎo)出一個(gè)基于地貌演化動(dòng)力方程的土壤厚度預(yù)測解析解能幫助水文學(xué)家理解土壤厚度隨地質(zhì)年代變化的規(guī)律，也為相應(yīng)的水文模擬研究提供了一個(gè)土壤厚度空間分布預(yù)測的實(shí)用工具。

附圖說明

圖1為三乘三表格的3D結(jié)構(gòu)圖。

圖2為計(jì)算曲率的柵格排序圖。

圖3為土壤通量f導(dǎo)出的計(jì)算網(wǎng)格。其中，q代表土壤輸移通量，灰色箭頭代表柵格i的三角刻面方向α₁,α₂分別是與基方向以及斜方向的夾角。

圖1、2、3即為本發(fā)明三乘三表格中土壤輸移通量演算圖。

具體實(shí)施方式

下面將結(jié)合本發(fā)明實(shí)施例中的附圖，對本發(fā)明中的技術(shù)方案進(jìn)行清楚、完整地描述。以下實(shí)施例僅用于更加清楚地說明本發(fā)明的技術(shù)方案，而不能以此來限制本發(fā)明的保護(hù)范圍。

一種基于柵格DEM的山丘區(qū)土壤厚度預(yù)測的方法

S1、根據(jù)土壤生成和輸移的質(zhì)量守恒，描述土壤演化的動(dòng)力學(xué)方程寫為：

h是土壤厚度，t為時(shí)間，e是高程，q是土壤輸移通量，η是巖石密度和土壤密度之比，即η＝ρ_r/ρ_s，方程(1)右邊的第一項(xiàng)是由于風(fēng)化作用，基巖上土壤生成的速率，Heimsath et al.[2000]研究表明土壤生成率隨土壤厚度的增加大致呈指數(shù)下降：

P₀是裸露基巖上土壤生成率(即土壤厚度為零時(shí)的土壤生成率)，h₀是一個(gè)特征侵蝕深度，為經(jīng)驗(yàn)參數(shù)；

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、根據(jù)McKean et al.[1993]，Small et al.[1999]，Dietrich et al[1995]等人的報(bào)告，在許多濕潤，半濕潤流域，坡度是土壤輸移主要的主導(dǎo)因素，即土壤蠕動(dòng)主導(dǎo)土壤輸移的模型，此外，這些線性輸移法則會(huì)得出質(zhì)量守恒方程的一個(gè)簡單解析解。所以，在本發(fā)明中，使用兩個(gè)線性沉積物輸移法則，即一個(gè)土壤蠕動(dòng)輸移模型q_d和一個(gè)徑流輸移模型q_t：

q＝q_d+q_t (3)

土壤蠕動(dòng)q_d用一個(gè)線性蠕變函數(shù)表示：

其中k_d是擴(kuò)散系數(shù)，z是土壤表面高程；

根據(jù)Willgoose et al.[1991b]and Saco et al.[2006]，降雨徑流驅(qū)動(dòng)下的土壤輸移一般方程寫為：

k_t是沉積物輸移系數(shù)，A是集水面積，m和n是一個(gè)給定地形的指數(shù)常數(shù)?；贔agherazzi的定義，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)，得到用來描述土壤厚度演化的土壤輸移動(dòng)力學(xué)方程：

其中：c＝ηP₀，且：

其中f是土壤蠕動(dòng)和降雨徑流土運(yùn)的下坡土壤通量散度，為了根據(jù)方程(6)和(7)計(jì)算f，要利用基于數(shù)字高程模型(DEM)的數(shù)字水系。如圖1所述，演示了怎樣通過柵格網(wǎng)絡(luò)計(jì)算土壤輸移通量。根據(jù)方程(4)，方程(7)中的可表示為其中的值等于地形曲率，C；

圖2中柵格i的C_i根據(jù)Heimsath et al.[1999]可計(jì)算如下：

其中z_i是研究柵格i的高程，z₁-z₈是圖2中柵格i周圍八個(gè)柵格的高程，Δx系柵格的邊長大??；方程(7)中的可用一個(gè)坡度和上游集水面積的函數(shù)來表達(dá)：

和q_tok分別是柵格i的土壤入流和出流通量，是D∞method[Tarboton,1997]確定的比例系數(shù)，例如對于出流的系數(shù)和可定義為：

其中α₁,α₂是分別是圖3中與基方向以及斜方向的夾角；

方程(6)和(7)可通過數(shù)值方法求解，采用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引用該文獻(xiàn)中的方法。

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、為了進(jìn)一步簡化動(dòng)態(tài)方程提出兩個(gè)假設(shè)：首先，假定研究的對象僅限于濕潤和半濕潤的地區(qū)，這里壤中流較為豐富，是主要的徑流成分，且下墊面基巖機(jī)械性質(zhì)足夠強(qiáng)；其次，假設(shè)在相對較短的地質(zhì)年代尺度內(nèi)，沒有構(gòu)造隆升或者基巖下沉的情況，即地表地形在土壤厚度演化過程中相對穩(wěn)定或穩(wěn)態(tài)變化。因此，可以推定地形特征(如曲率)的演化受巖石風(fēng)化和土壤輸移作用的驅(qū)動(dòng)，且認(rèn)為其速率比土壤厚度演化的速率慢得多。這樣，方程(6)可以看作是一個(gè)一階非線性非齊次常微分方程。方程(6)的一般解可以通過以下的步驟導(dǎo)出。

首先，方程改寫為：

其中：

之后，改寫成一個(gè)線性常微分方程：

Z′＝aZ+b (13)

其中

方程(13)的兩邊都乘以e^-at并重組為：

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

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

對方程(14)進(jìn)行積分，得到線性常微分方程的一般解：

其中R₁是積分常量；

方程(12)代入方程(15)中，土壤厚度的一般表達(dá)就可導(dǎo)出為：

在方程(16)中，系數(shù)R₁可以通過設(shè)定初始土壤厚度來決定；所以，如果初始土壤厚度假設(shè)為：h(0)＝hi (17)

那么積分常量R₁導(dǎo)出為：

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

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

h(0)＝0 (20)

那么常數(shù)R₁為：

從而，方程(16)可改寫成一個(gè)更簡化的表達(dá)，如下：

其中h₀,c是土壤生成參數(shù)(可以根據(jù)他們物理意義定義它的值)，f是整體土壤侵蝕率，可以根據(jù)方程(7)決定；

由于土壤厚度一般在流域的低洼地增長，根據(jù)方程(2)土壤生成率為零，那么方程(6)重寫成：

對方程(23)積分，積分時(shí)間從t_s到t，土壤厚度從h_s到h，得到一個(gè)可以預(yù)測土壤厚度的方程：

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

其中t_s是土壤生成率為零的時(shí)間，h_s是時(shí)刻t_s的土壤厚度。

在此分析模型中，土壤生成方程即方程(2)中的參數(shù)P₀根據(jù)放射性同位素測量。

盡管方程(4)和(5)中的擴(kuò)散系數(shù)k_d在土壤輸移公式中帶有經(jīng)驗(yàn)的特性，但也可以理論地根據(jù)現(xiàn)場實(shí)驗(yàn)來確定，如果河道河床遷移速率可以評估的話。根據(jù)Chiang和Hsu[2006]已公開的研究成果，方程(4)中的參數(shù)擴(kuò)散系數(shù)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].

以上所述僅是本發(fā)明的優(yōu)選實(shí)施方式，應(yīng)當(dāng)指出，對于本技術(shù)領(lǐng)域的普通技術(shù)人員來說，在不脫離本發(fā)明技術(shù)原理的前提下，還可以做出若干改進(jìn)和變形，這些改進(jìn)和變形也應(yīng)視為本發(fā)明的保護(hù)范圍。