基于DEM的流域地貌特征分形量化研究

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

(1.西安理工大學 西北旱區生態水利工程國家重點實驗室, 西安 710048;2.陜西省水利廳, 西安 710004; 3.陜西省水利電力勘測設計研究院, 西安 710001)

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

1 原理與方法

1.1 分形布朗運動理論

分形布朗運動理論是由Mandelbrot提出的描述自然界中隨機分形的一種數學模型[11]。分形布朗運動的特點是具有統計自相似性即自仿射性(≈表示兩者的概率分布相同)：

(1)

在應用FBM模型計算分形數據對象的分形維數時，關鍵是估計模型的非規則維參數即H參數。H參數估計是否準確關系到FBM模型對應用對象的適用性。采用絕對矩估計法，由公式(1)得到：

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

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

(2)

兩邊取對數得到：

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

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

(3)

實際上等于：y=H·x+C。再由分形維數FD和Hurst指數H的關系得：

FD=DE-H

(4)

式中：DE為分形體的歐氏拓撲維數。

1.2 地貌分形量化模型

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

對于DEM來說，s=(x，y)∈E2，代表DEM上各像元點位置。t=Z(s)，代表DEM上各像元點的高程屬性值。若令A(t,s)代表DEM的統計表面積，則由分形布朗運動理論可將其推廣到二維的情形，推導出一個與表面積有關的性質為：

(5)

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

(6)

(7)

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

(8)

式中：FD為基于表面積的FBM流域地貌因子，DE為流域空間曲面的歐氏維數，對于表面積來說DE=2；r為DEM像元邊長，A(t,rs)為像元尺度r下所測得的流域空間曲面表面積。

1.3 地貌特征分形量化的計算過程

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

(9)

式中：Z(i,j)表示某DEM在尺度r下的高程值(i=1，2，…，M；j=1，2，…，N)。選定一條流域為研究對象，構建GIS實體模型，利用公式(9)就能計算DEM地形曲面的表面積A(t,s)。

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

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

2 流域概況及計算過程

2.1 流域概況

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

2.2 數據源與像元尺寸設置

在數字高程模型構建過程中，釆樣是關鍵的一環。確定最優的數據采樣密度和地表重建方法，都需要對地形表面形態特征有深刻地認識。本研究地貌形態GIS模型的創建選擇由紙質地形圖創建的方式，經過地形圖掃描、幾何糾正、影像二值處理與細化、等高線矢量化、數據接邊、構建不規則三角網(TIN)、DEM的生成等步驟。

為了確保擁有足夠的分析數據和DEM研究精度，這里將該流域地貌形態GIS實體模型的像元尺度分別設置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種規格，共生成14個小流域DEM模型。

2.3 計算過程

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

逐個計算陳家畈流域在不同像元尺度下DEM地形曲面的表面積A(t,rs)，再對變化的尺度r和對應的DEM表面積A(t,s)取對數，計算結果見表1。

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

以lnh為橫坐標軸，以lnA(t,s)為縱坐標軸將計算結果點繪在直角坐標系上，確定相關系數最高的區間即像元尺度18 m×18 m—36 m×36 m為無標度區間(R2≥0.99)，通過在該無標度區間內觀察利用最小二乘法進行一階線性擬合，見公式(10)：

(10)

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

3 結果與分析

3.1 地貌形態特征的標度不變性

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

圖2 擬和直線示意圖

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

因此，無標度區間的確定也就為不同像元尺度的流域地貌特征之間的比較與評價提供了可能性。另外由于18 m×18 m—36 m×36 m范圍剛好處于常規高精尺度地形數據的尺度范圍，這樣也為高精尺度分布式地貌量化因子的開發和推廣應用提供了新思路。

3.2 地貌形態特征的綜合性表達

通過以上分形計算過程可以看出，流域地貌特征分形量化因子能夠對流域三維地貌特征進行綜合表達。lnA(t,s)實際上反映了整個流域表面的起伏特征，是對整個流域表面復雜度的總體表征。對于具體流域而言，其像元尺寸一定時所對應的lnA(t,s)值越大，流域地貌的起伏變化越強烈，流域表面越復雜，反之亦然。因此lnA(t,s)與lnr的增量比值，即流域地貌特征分形量化因子反映的是整個流域的總體性綜合性特征。

準確有效地評價流域地貌特征，有利于了解和辨別地形地貌演變的客觀規律和整體特征。一般來說，流域地貌形態分形量化因子D的取值范圍為2～3。當D的取值為2時，流域表面趨于二維平面，流域地貌形態最簡單。當D的取值為3時，流域表面趨于三維立方體，流域地貌形態也相對簡單。只有當D=2.5時，流域地貌形態處于布朗隨機運動狀態，此時地貌特征處于最雜亂和無序狀態，即地貌形態破碎和復雜程度的臨界值，越遠離D=2.5流域地貌特征的破碎程度越低。兩種變化方式表征意義不同，當D從2.5增長到3時，地貌特征的演變處于正向堆積隆升階段，趨于形成規整的塬臺地貌。當D從2.5減小到2時，地貌特征的演變處于負向侵蝕拗陷階段，趨于形成平坦的平原地貌。

3.3 地貌形態特征與傳統地貌因子的關系

為了更加有效地對流域地貌形態特征FBM量化因子進行研究。在ArcGIS平臺下，分別計算陳家畈流域在18 m×18 m—36 m×36 m像元尺度內的溝壑密度、地形起伏度、平均坡度、流域中值高程等地貌形態單因子指標，并與FBM流域地貌因子進行對比分析。可以看出當像元尺度處于18 m×18 m—36 m×36 m范圍時，溝壑密度、地形起伏度、平均坡度、流域中值高程等地貌形態單因子指標分別呈現不同的變化形態。

隨著像元尺度的增大，溝壑密度呈現先減小后增大的趨勢，地形起伏度不斷增大，平均坡度不斷減小，流域中值高程呈現先平穩后增大最后減小的趨勢。然而在此過程中FBM量化因子的值始終穩定在2.012 7處，其幾何形態始終處于其他各因子曲線的中間區域，表現出與其他各因子曲線幾何形態的“均值”形態相似的特征。亦即FBM量化因子蘊涵了各個地貌形態單因子的“綜合疊加”特點，因此是對流域地貌形態特征的綜合性表達。

4 結 語

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