陳 兵, 彭 芳, 李 鵬, 袁水龍

(1.西安理工大學(xué) 西北旱區(qū)生態(tài)水利工程國(guó)家重點(diǎn)實(shí)驗(yàn)室, 西安 710048;2.陜西省水利廳, 西安 710004; 3.陜西省水利電力勘測(cè)設(shè)計(jì)研究院, 西安 710001)

小流域是黃土高原生態(tài)環(huán)境恢復(fù)重建與治理的基本單元，其侵蝕產(chǎn)沙規(guī)律和水土流失防治一直是土壤侵蝕與水土保持學(xué)界研究的重點(diǎn)[1-3]。地貌形態(tài)是黃土高原流域降雨侵蝕產(chǎn)沙過(guò)程的重要下墊面影響因素之一，它在降雨動(dòng)力、人為活動(dòng)和地質(zhì)構(gòu)造運(yùn)動(dòng)等因素的共同作用下，通過(guò)水沙運(yùn)移、能量消耗等方式實(shí)現(xiàn)其形態(tài)的不斷演化，最終形成地形破碎、溝壑縱橫的侵蝕地貌景觀[4-5]。在已有的流域水土流失評(píng)價(jià)和預(yù)報(bào)模型中，諸多學(xué)者將流域高差比、溝壑密度、平均坡度、平均長(zhǎng)度及溝谷切割深度等因子中的一個(gè)或多個(gè)作為流域水土流失預(yù)報(bào)模型中地形因子的量化參數(shù)[6-9]。但是由于地貌系統(tǒng)是一個(gè)非線性動(dòng)態(tài)系統(tǒng)(Nonlinear Dynamical System，NDS)[10]，傳統(tǒng)的地貌因子難以有效地表達(dá)這個(gè)非線性系統(tǒng)的復(fù)雜性，如地貌演化過(guò)程中的自組織性、突現(xiàn)性、自相似性、多尺度性及時(shí)空耦合性等特征，因此在流域水土流失評(píng)價(jià)和預(yù)報(bào)研究領(lǐng)域存在明顯的局限性。分形理論的形成與發(fā)展為流域地形因子的整體性描述和地貌過(guò)程的非線性研究開(kāi)辟了新的思路。但就目前而言，流域地貌的分形大多集中于水系分形研究，而對(duì)整個(gè)流域地貌形態(tài)特征的綜合描述、地貌因子分形建模及其尺度范圍(即無(wú)標(biāo)度區(qū))等方面問(wèn)題有待進(jìn)一步研究。本文以地理信息系統(tǒng)為平臺(tái)，建立基于DEM表面積的流域分形布朗運(yùn)動(dòng)FBM地貌因子計(jì)算模型，并進(jìn)行相關(guān)的GIS算法開(kāi)發(fā)，以期能為黃土高原小流域降雨侵蝕產(chǎn)沙預(yù)報(bào)模型中的地貌因子量化研究提供新思路與新角度。

1 原理與方法

1.1 分形布朗運(yùn)動(dòng)理論

分形布朗運(yùn)動(dòng)理論是由Mandelbrot提出的描述自然界中隨機(jī)分形的一種數(shù)學(xué)模型[11]。分形布朗運(yùn)動(dòng)的特點(diǎn)是具有統(tǒng)計(jì)自相似性即自仿射性(≈表示兩者的概率分布相同)：

(1)

在應(yīng)用FBM模型計(jì)算分形數(shù)據(jù)對(duì)象的分形維數(shù)時(shí)，關(guān)鍵是估計(jì)模型的非規(guī)則維參數(shù)即H參數(shù)。H參數(shù)估計(jì)是否準(zhǔn)確關(guān)系到FBM模型對(duì)應(yīng)用對(duì)象的適用性。采用絕對(duì)矩估計(jì)法，由公式(1)得到：

E[|BH(t+hs)-BH(t)|]=

‖h‖H·E[|BH(t+s)-BH(t)|]

(2)

兩邊取對(duì)數(shù)得到：

lnE[|BH(t+hs)-BH(t)|]=

H·ln‖h‖+lnE[|BH(t+s)-BH(t)|]

(3)

實(shí)際上等于：y=H·x+C。再由分形維數(shù)FD和Hurst指數(shù)H的關(guān)系得：

FD=DE-H

(4)

式中：DE為分形體的歐氏拓?fù)渚S數(shù)。

1.2 地貌分形量化模型

數(shù)字高程模型(Digital Elevation Model，DEM)是描述地面三維信息的常用模型，其數(shù)據(jù)結(jié)構(gòu)是一種基于柵格模型的空間數(shù)據(jù)結(jié)構(gòu)[12]。它將整個(gè)空間分割成大小相等緊密相鄰的正方形網(wǎng)格，稱為像元(cell)，高程信息以屬性值的形式儲(chǔ)存在相應(yīng)的像元中。依托DEM進(jìn)行數(shù)字地形分析，能夠有效挖掘更深層的地形特征和地學(xué)現(xiàn)象，獲取地形屬性信息。隨著地理信息系統(tǒng)的應(yīng)用和普及，數(shù)字地形分析技術(shù)的優(yōu)勢(shì)和應(yīng)用前景逐步為人們所認(rèn)識(shí)，已成為進(jìn)行地形模擬和地學(xué)分析的核心技術(shù)之一，在測(cè)繪、交通、軍事、規(guī)劃等相關(guān)領(lǐng)域被廣泛應(yīng)用。

對(duì)于DEM來(lái)說(shuō)，s=(x，y)∈E2，代表DEM上各像元點(diǎn)位置。t=Z(s)，代表DEM上各像元點(diǎn)的高程屬性值。若令A(yù)(t,s)代表DEM的統(tǒng)計(jì)表面積，則由分形布朗運(yùn)動(dòng)理論可將其推廣到二維的情形，推導(dǎo)出一個(gè)與表面積有關(guān)的性質(zhì)為：

(5)

lnA(t,rs)=H·lnr+lnF(t,s)

(6)

(7)

綜上提出基于表面積的地貌特征分形量化模型為：

(8)

式中：FD為基于表面積的FBM流域地貌因子，DE為流域空間曲面的歐氏維數(shù)，對(duì)于表面積來(lái)說(shuō)DE=2；r為DEM像元邊長(zhǎng)，A(t,rs)為像元尺度r下所測(cè)得的流域空間曲面表面積。

1.3 地貌特征分形量化的計(jì)算過(guò)程

規(guī)則網(wǎng)格DEM將地形曲面劃分成一系列的規(guī)則網(wǎng)格單元，每個(gè)網(wǎng)格單元對(duì)應(yīng)一個(gè)高程值。因此可以將其理解成由大量相同的小立方體單元“堆積”而成，小立方體邊長(zhǎng)等于像元的邊長(zhǎng)r，小立方體的高度為高程的屬性值Z。計(jì)算每個(gè)小立方體的頂部水平面積及4個(gè)側(cè)面的垂直面積之和就得到了整個(gè)DEM的表面積A(t,s)。為了避免重復(fù)計(jì)算，對(duì)每個(gè)像元只計(jì)算前面和左面的垂直面。r為像元尺度，AH為相鄰像元高程之差，其中AH=r×r；AV1，AV2為相鄰像元高程值之差與尺度r的乘積。設(shè)DEM的尺寸為M×N(像元個(gè)數(shù))，DEM的表面積A(t,s)的計(jì)算如下[13]：

(9)

式中：Z(i,j)表示某DEM在尺度r下的高程值(i=1，2，…，M；j=1，2，…，N)。選定一條流域?yàn)檠芯繉?duì)象，構(gòu)建GIS實(shí)體模型，利用公式(9)就能計(jì)算DEM地形曲面的表面積A(t,s)。

式中基于表面積的FBM地貌因子的計(jì)算，實(shí)質(zhì)上是將DEM像元尺寸逐步趨近于零計(jì)算極限的過(guò)程。然而在實(shí)際計(jì)算過(guò)程中所選定的DEM精度總是有限的，因此必須對(duì)所研究流域DEM的像元尺寸進(jìn)行適當(dāng)?shù)倪x擇(圖1)。在尺度r下，求取表面積A(t，s)；變換像元尺度r，求出若干個(gè)相應(yīng)的A(t，s)；將A(t，s)與r同取對(duì)數(shù)，并對(duì)兩數(shù)列進(jìn)行曲線擬合；利用統(tǒng)計(jì)軟件求出lnA(t，s)相對(duì)于lnr的一階導(dǎo)數(shù)，得到FBM地貌因子的Hurst指數(shù)H；最后由公式(8)計(jì)算得到基于表面積的FBM流域地貌因子FD。

圖1 DEM像元尺度變化示意圖

2 流域概況及計(jì)算過(guò)程

2.1 流域概況

陳家畈流域位于黃土高原中部黃土丘陵溝壑區(qū)第一副區(qū)，是大理河流域的一級(jí)子流域。該流域干流全長(zhǎng)約為28 km，流域面積282 km2。流域地勢(shì)西高東低，平均海拔1 100～1 700 m。流域地形復(fù)雜破碎，植被稀疏，水土流失嚴(yán)重。經(jīng)過(guò)長(zhǎng)期侵蝕，形成梁峁突起，溝壑縱橫，鑲嵌有墹地的黃土梁峁丘陵溝壑地貌。本研究擬以大理河陳家畈流域?yàn)槔?，?duì)流域地貌特征量化因子的計(jì)算方法進(jìn)行研究，并對(duì)FBM地貌因子的特征及應(yīng)用進(jìn)行探討。

2.2 數(shù)據(jù)源與像元尺寸設(shè)置

在數(shù)字高程模型構(gòu)建過(guò)程中，釆樣是關(guān)鍵的一環(huán)。確定最優(yōu)的數(shù)據(jù)采樣密度和地表重建方法，都需要對(duì)地形表面形態(tài)特征有深刻地認(rèn)識(shí)。本研究地貌形態(tài)GIS模型的創(chuàng)建選擇由紙質(zhì)地形圖創(chuàng)建的方式，經(jīng)過(guò)地形圖掃描、幾何糾正、影像二值處理與細(xì)化、等高線矢量化、數(shù)據(jù)接邊、構(gòu)建不規(guī)則三角網(wǎng)(TIN)、DEM的生成等步驟。

為了確保擁有足夠的分析數(shù)據(jù)和DEM研究精度，這里將該流域地貌形態(tài)GIS實(shí)體模型的像元尺度分別設(shè)置18 m×18 m，20 m×20 m，22 m×22 m，24 m×24 m，26 m×26 m，28 m×28 m，30 m×30 m，32 m×32 m，34 m×34 m，36 m×36 m，38 m×38 m，40 m×40 m，45 m×45 m，50 m×50 m共14種規(guī)格，共生成14個(gè)小流域DEM模型。

2.3 計(jì)算過(guò)程

首先選取陳家畈流域18 m×18 m像元尺度DEM為研究對(duì)象，利用公式(9)求取表面積A(t,s)；然后依次變換像元尺度從20 m×20 m經(jīng)過(guò)14次尺度變化到50 m×50 m，求出相應(yīng)的A(t,s)；將所得表面積A(t,s)與像元尺度r同取對(duì)數(shù)，并對(duì)兩數(shù)列進(jìn)行線性擬合；確定相關(guān)系數(shù)最高的區(qū)間為無(wú)標(biāo)度區(qū)間(R2≥0.99)；利用統(tǒng)計(jì)軟件SPSS求出lnA(t,s)相對(duì)于lnr的一階導(dǎo)數(shù)，得到地貌因子的Hurst指數(shù)H；最后由公式(8)計(jì)算得到陳家畈流域基于表面積的流域地貌特征分形量化因子 。

逐個(gè)計(jì)算陳家畈流域在不同像元尺度下DEM地形曲面的表面積A(t,rs)，再對(duì)變化的尺度r和對(duì)應(yīng)的DEM表面積A(t,s)取對(duì)數(shù)，計(jì)算結(jié)果見(jiàn)表1。

表1 陳家畈流域地貌特征分形量化的計(jì)算成果

以lnh為橫坐標(biāo)軸，以lnA(t,s)為縱坐標(biāo)軸將計(jì)算結(jié)果點(diǎn)繪在直角坐標(biāo)系上，確定相關(guān)系數(shù)最高的區(qū)間即像元尺度18 m×18 m—36 m×36 m為無(wú)標(biāo)度區(qū)間(R2≥0.99)，通過(guò)在該無(wú)標(biāo)度區(qū)間內(nèi)觀察利用最小二乘法進(jìn)行一階線性擬合，見(jiàn)公式(10)：

(10)

將計(jì)算結(jié)果代入公式(8)中“基于表面積的地貌特征分形量化模型”，就可以計(jì)算得到陳家畈流域基于DEM表面積的流域地貌特征分形量化因子為FD=DE-H=2.0127。

3 結(jié)果與分析

3.1 地貌形態(tài)特征的標(biāo)度不變性

標(biāo)度不變性是分形量化模型應(yīng)用的最大優(yōu)勢(shì)。從圖2可以看出，當(dāng)像元尺度大致處于18 m×18 m—36 m×36 m區(qū)間時(shí)，lnA(r)-lnr擬合的一階線性方程的決定系數(shù)R2均達(dá)到0.99以上。說(shuō)明在此區(qū)間內(nèi)ln 與lnr表現(xiàn)出良好的相關(guān)性，并且呈現(xiàn)出明顯的線性分布特征。此時(shí)lnA(r)與lnr的增量比值不隨lnA(r)的變化而變化，呈現(xiàn)出標(biāo)度不變性特征。lnA(r)與lnr的比值具有統(tǒng)計(jì)意義上的自相似性，地貌形態(tài)表現(xiàn)出了分形性質(zhì)，此區(qū)間應(yīng)處在無(wú)標(biāo)度區(qū)之內(nèi)。當(dāng)像元尺度處在36 m×36 m以上時(shí)，線形的總體走向趨于下彎，lnA(r)與lnr的比值不能保持穩(wěn)定。lnA(r)-lnr的線性分布特性消失，此時(shí)地貌形態(tài)的分形性質(zhì)也隨之消失，此區(qū)間應(yīng)處于無(wú)標(biāo)度區(qū)之外。

圖2 擬和直線示意圖

無(wú)標(biāo)度區(qū)間的確定是分形維數(shù)應(yīng)用的重要前提條件。在進(jìn)行分形維數(shù)比較時(shí)應(yīng)特別注意，無(wú)標(biāo)度區(qū)間的確定方法必須是統(tǒng)一的，另外只能對(duì)處于同一個(gè)無(wú)標(biāo)度區(qū)間內(nèi)的分維進(jìn)行比較。由圖2可知，在18 m×18 m—36 m×36 m范圍內(nèi)(R2≥0.99)，陳家畈流域基于DEM表面積的FBM流域地貌因子為2.012 7，即當(dāng)像元尺度處于這個(gè)范圍時(shí)，流域地貌特征可以用一個(gè)分形維數(shù)表征。

因此，無(wú)標(biāo)度區(qū)間的確定也就為不同像元尺度的流域地貌特征之間的比較與評(píng)價(jià)提供了可能性。另外由于18 m×18 m—36 m×36 m范圍剛好處于常規(guī)高精尺度地形數(shù)據(jù)的尺度范圍，這樣也為高精尺度分布式地貌量化因子的開(kāi)發(fā)和推廣應(yīng)用提供了新思路。

3.2 地貌形態(tài)特征的綜合性表達(dá)

通過(guò)以上分形計(jì)算過(guò)程可以看出，流域地貌特征分形量化因子能夠?qū)α饔蛉S地貌特征進(jìn)行綜合表達(dá)。lnA(t,s)實(shí)際上反映了整個(gè)流域表面的起伏特征，是對(duì)整個(gè)流域表面復(fù)雜度的總體表征。對(duì)于具體流域而言，其像元尺寸一定時(shí)所對(duì)應(yīng)的lnA(t,s)值越大，流域地貌的起伏變化越強(qiáng)烈，流域表面越復(fù)雜，反之亦然。因此lnA(t,s)與lnr的增量比值，即流域地貌特征分形量化因子反映的是整個(gè)流域的總體性綜合性特征。

準(zhǔn)確有效地評(píng)價(jià)流域地貌特征，有利于了解和辨別地形地貌演變的客觀規(guī)律和整體特征。一般來(lái)說(shuō)，流域地貌形態(tài)分形量化因子D的取值范圍為2～3。當(dāng)D的取值為2時(shí)，流域表面趨于二維平面，流域地貌形態(tài)最簡(jiǎn)單。當(dāng)D的取值為3時(shí)，流域表面趨于三維立方體，流域地貌形態(tài)也相對(duì)簡(jiǎn)單。只有當(dāng)D=2.5時(shí)，流域地貌形態(tài)處于布朗隨機(jī)運(yùn)動(dòng)狀態(tài)，此時(shí)地貌特征處于最雜亂和無(wú)序狀態(tài)，即地貌形態(tài)破碎和復(fù)雜程度的臨界值，越遠(yuǎn)離D=2.5流域地貌特征的破碎程度越低。兩種變化方式表征意義不同，當(dāng)D從2.5增長(zhǎng)到3時(shí)，地貌特征的演變處于正向堆積隆升階段，趨于形成規(guī)整的塬臺(tái)地貌。當(dāng)D從2.5減小到2時(shí)，地貌特征的演變處于負(fù)向侵蝕拗陷階段，趨于形成平坦的平原地貌。

3.3 地貌形態(tài)特征與傳統(tǒng)地貌因子的關(guān)系

為了更加有效地對(duì)流域地貌形態(tài)特征FBM量化因子進(jìn)行研究。在ArcGIS平臺(tái)下，分別計(jì)算陳家畈流域在18 m×18 m—36 m×36 m像元尺度內(nèi)的溝壑密度、地形起伏度、平均坡度、流域中值高程等地貌形態(tài)單因子指標(biāo)，并與FBM流域地貌因子進(jìn)行對(duì)比分析?？梢钥闯霎?dāng)像元尺度處于18 m×18 m—36 m×36 m范圍時(shí)，溝壑密度、地形起伏度、平均坡度、流域中值高程等地貌形態(tài)單因子指標(biāo)分別呈現(xiàn)不同的變化形態(tài)。

隨著像元尺度的增大，溝壑密度呈現(xiàn)先減小后增大的趨勢(shì)，地形起伏度不斷增大，平均坡度不斷減小，流域中值高程呈現(xiàn)先平穩(wěn)后增大最后減小的趨勢(shì)。然而在此過(guò)程中FBM量化因子的值始終穩(wěn)定在2.012 7處，其幾何形態(tài)始終處于其他各因子曲線的中間區(qū)域，表現(xiàn)出與其他各因子曲線幾何形態(tài)的“均值”形態(tài)相似的特征。亦即FBM量化因子蘊(yùn)涵了各個(gè)地貌形態(tài)單因子的“綜合疊加”特點(diǎn)，因此是對(duì)流域地貌形態(tài)特征的綜合性表達(dá)。

4 結(jié) 語(yǔ)

流域地貌特征的綜合量化指標(biāo)不僅是數(shù)字地形分析的重要參數(shù)，而且對(duì)地形數(shù)據(jù)壓縮、地形分類、測(cè)繪可視化、雪線分布、生物多樣性、地形仿真、土地利用、導(dǎo)航定位、空間分析等研究領(lǐng)域有重要的研究意義，備受國(guó)內(nèi)外研究人員的關(guān)注。本文提出的基于DEM流域地貌特征分形量化因子計(jì)算模型，利用分形維數(shù)的標(biāo)度不變性，實(shí)現(xiàn)了一定的像元尺度范圍內(nèi)流域地貌形態(tài)特征量化指標(biāo)的尺度不變，從而可以實(shí)現(xiàn)不同像元尺度下地貌形態(tài)特征的比較與評(píng)價(jià)。FBM地貌因子克服了傳統(tǒng)地貌形態(tài)單因子指標(biāo)的不足，可以對(duì)流域地貌形態(tài)特征進(jìn)行綜合性表達(dá)，為地貌因子量化研究提供了新思路與新角度。但是由于流域地貌特征分形量化計(jì)算模型的地域適應(yīng)性還需要經(jīng)過(guò)更多地貌類型的驗(yàn)證，后續(xù)研究重點(diǎn)應(yīng)該繼續(xù)放在對(duì)其他地貌類型區(qū)的檢驗(yàn)方面，并在此基礎(chǔ)上更加深入的探討分形量化模型的普適性問(wèn)題。