DEM地表坡向變率的向量幾何計算法

胡光輝，熊禮陽，湯國安

1. 南京師范大學(xué)地理科學(xué)學(xué)院，江蘇 南京 210023； 2. 南京師范大學(xué)虛擬地理環(huán)境教育部重點實驗室，江蘇 南京 210023； 3. 江蘇省地理信息資源開發(fā)與利用協(xié)同創(chuàng)新中心，江蘇 南京 210023

坡向(aspect)是重要的地形因子之一，決定著地表接收的太陽能量強度與地表徑流方向，對于山地生態(tài)具有重要的作用[1-3]。坡向的變化也代表了地表接收能量的改變和坡面形態(tài)的轉(zhuǎn)換。坡向變率(slope of aspect,SOA)的概念即由此提出，指坡面指向的變化程度[4]。坡向變率在地學(xué)中具有極其重要的科學(xué)與實踐意義。從地貌學(xué)角度來看，坡面指向的轉(zhuǎn)折映射著不同坡面形態(tài)的轉(zhuǎn)換，并據(jù)此產(chǎn)生眾多地貌特征線，如山脊線，山谷線等[5-6]。從生態(tài)學(xué)角度來看，坡面指向的變化表征著地面接收太陽輻射能量的轉(zhuǎn)換與地表植被的變化，即陰陽坡產(chǎn)生的植被生態(tài)變化[7]。此外，風(fēng)積過程中迎風(fēng)坡與背風(fēng)坡的差異也體現(xiàn)出坡面指向的變化[8]。可見，科學(xué)準確地計算坡向變率是地貌學(xué)、生態(tài)學(xué)等眾多學(xué)科相關(guān)研究的重要內(nèi)容。

20世紀以來，數(shù)字高程模型(DEM)與數(shù)字地形分析(DTA)的提出與應(yīng)用，為傳統(tǒng)的地學(xué)分析方法帶來革命性的變化[9-13]。現(xiàn)階段，基于DEM數(shù)據(jù)，一階地形因子能夠較為方便地計算[14]。但是，坡向變率為二階地形因子，其計算方法更為復(fù)雜[15]。起初，坡向變率的計算是參考坡度的計算方法(對高程數(shù)值矩陣求偏導(dǎo))，也就是對坡向數(shù)值矩陣求取偏導(dǎo)數(shù)，即計算坡向的坡度。據(jù)此，隨著對坡向變化理解的深入，前人依次提出直接法、正反DEM法、六分法實現(xiàn)了標量條件的坡向變率求解，并取得了廣泛的應(yīng)用[16-19]?？梢姡叭藢ζ孪蜃兟实挠嬎阕髁顺醪降难芯?，取得了一定的研究成果。但是，坡向變率的計算是一個基于坡向數(shù)據(jù)基礎(chǔ)二階地形因子的提取[20-21]。該坡向數(shù)值矩陣與原始的高程數(shù)值矩陣存在本質(zhì)的不同，即高程數(shù)值在很大程度上被看成是一個沒有方向特征的常規(guī)數(shù)學(xué)標量，標量值只有大小關(guān)系，而無須考慮方向?qū)傩?。此時對其進行一階求偏導(dǎo)可直接進行高程求差，即一階地形因子。而相對二階地形因子坡向變率而言，坡向數(shù)值矩陣的每個像元都代表著其獨有的坡面指向，這種指向代表著其特有的方向性[22]。因此，坡向數(shù)值矩陣不是一個標量矩陣，簡單數(shù)學(xué)標量的代數(shù)運算方式用在坡向數(shù)值矩陣上不可避免地違背了數(shù)學(xué)運算機制。

從數(shù)學(xué)上看，標量運算法則針對的是沒有方向?qū)傩詳?shù)據(jù)的簡單數(shù)值計算，而向量運算法則恰好是針對帶有方向?qū)傩詳?shù)據(jù)的數(shù)學(xué)向量計算。向量既有大小，也有方向，其運算法則遵循向量運算法則，運算結(jié)果是一個向量(圖1(a))，而不是一個標量。因此，坡向變率的計算應(yīng)該遵循向量運算法則而不是標量運算法則。針對該思路，本文提出基于數(shù)學(xué)向量幾何的坡向變率計算方法。該方法首先以DEM數(shù)據(jù)為基礎(chǔ)提取坡向矩陣數(shù)據(jù)。其次，對坡向矩陣數(shù)據(jù)進行數(shù)學(xué)向量表達。最后，依據(jù)該向量化的坡向矩陣來計算坡向變率。并將結(jié)果與傳統(tǒng)方法進行對比，以期更為科學(xué)準確地認識坡向以及計算坡向變率。

1 研究方法

1.1 坡向變率計算基礎(chǔ)

坡向變率，即坡向的變化率。變化率在DEM上應(yīng)用最廣泛的是坡度的概念，即求一階偏導(dǎo)。坡向變率的計算即采用坡向矩陣代替高程矩陣，求取地表某一個點在鄰域內(nèi)的坡向變化率。以高程矩陣為基礎(chǔ)的DEM數(shù)據(jù)，高程變化率提取算法眾多。其中，由于三階反距離平方權(quán)差分算法相對合理準確，在研究中已被廣泛采用[23]。將此方法借鑒到坡向變率計算過程中，可以得到坡向變率的計算公式

(1)

式中，fx與fy分別是地表某一個點的坡向在鄰域內(nèi)東西和南北方向上的變化率。在三階反距離平方權(quán)差分算法中，鄰域柵格單元內(nèi)格網(wǎng)數(shù)值的代數(shù)做差運算是很重要的部分。由于傳統(tǒng)方法進行差分運算使用的是未經(jīng)向量化轉(zhuǎn)換的坡向數(shù)值數(shù)據(jù)，運算法則使用的是標量之間的代數(shù)運算法則，忽略了方向?qū)\算的影響，其實質(zhì)是一維運算(圖1(b))。在本文研究中，對坡向隱式表達的坡向數(shù)值矩陣數(shù)據(jù)進行向量化表達。坡向矩陣中每個坡向值都有其方向。此時，坡向做差運算不是一維運算，而是符合坡向定義的二維向量運算(圖1(a))。

1.2 坡向矩陣向量化表達

坡向變率向量計算方法的核心問題是對隱式的坡向方向特征進行向量表達以及向量運算。本文對坡向數(shù)據(jù)進行向量化表達的方法具體為以下步驟：

(1) 坡向在坡向極坐標系下的向量表達。坡向向量可以認為處在一個以正北方向為起始方向，順時針為旋轉(zhuǎn)方向的特殊極坐標系下。本文稱之為坡向極坐標系，如圖2(a)所示。坡向值應(yīng)當首先在坡向極坐標系下表達為向量。其表示為

A=(r,θ)

(2)

式中，r是長度，這里是柵格單元的大??；θ是旋轉(zhuǎn)角度，這里是坡向值；A是坡向在坡向極坐標系下的表示。

圖2 兩種不同的極坐標系Fig.2 Two different polar coordinate systems

(2) 坡向在普通極坐標系下的向量表達。如圖2(b)所示，普通極坐標系是較為常用的一種坐標系，并且可以與平面直角坐標系進行快速轉(zhuǎn)換。坡向極坐標轉(zhuǎn)為普通極坐標的方法如下

(3)

式中，θ是坡向值。轉(zhuǎn)換完成之后，每一個地面點的坡向向量在普通極坐標系下表示為

(4)

式中，r是柵格單元大??；θ是坡向值；AP是坡向在普通極坐標系下的表示。

(3) 坡向在平面直角坐標系下的向量表達。平面直角坐標系更適合于向量之間的運算，因此需要將極坐標系下的坡向向量轉(zhuǎn)換到平面直角坐標系下進行表達。轉(zhuǎn)換公式如下

(5)

式中，r是柵格單元大?。沪潦茿P的旋轉(zhuǎn)角度。坡向向量最終在平面直角坐標系下表示為AR=(x,y)。

(4) 向量在平面直角坐標系下的坐標表示。分別取與X軸、Y軸方向相同的兩個單位向量i、j作為基底向量。任作一個向量a，由平面向量基本定理可知，有且只有一對實數(shù)x、y，可以使得等式a=xi+yj成立，于是將(x,y)叫作向量a的直角坐標表示，記作a=(x,y)，如圖3所示。

圖3 平面直角坐標系下向量表示特征Fig.3 Vector representation method of plane Cartesian coordinate system

1.3 坡向變率向量計算

在平面直角坐標系下，平面向量間的運算法則如下：

(1) 已知a=(x1,y1),b=(x2,y2)，則a+b=(x1+x1,y2+y2)。

(2) 已知a=(x1,y1),b=(x2,y2)，則a-b=(x1-x2,y2-y2)。

(3) 已知a=(x,y)和實數(shù)λ,則λa=(λx,λy)。

通過以上向量之間的4個運算法則，研究中可完整、準確地將三階反距離平方權(quán)差分算法采用向量運算法則進行求解，并應(yīng)用于坡向變率的計算。圖4為基于向量幾何法計算坡向變率的完整技術(shù)流程。

2 試驗樣區(qū)與數(shù)據(jù)

本文采用數(shù)學(xué)高斯曲面和陜北黃土高原3種典型地貌塬、墚、峁所在地區(qū)為研究樣區(qū)，用于測試基于向量幾何法計算坡向變率。其中，高斯曲面為模擬數(shù)學(xué)地形曲面，高斯合成曲面參數(shù)方程的定義如式(6)所示[22]

(6)

式中，A、B、C為地勢起伏參數(shù)；m、n為范圍控制參數(shù)。使用基于數(shù)學(xué)公式建立起來的DEM數(shù)據(jù)具有無誤差的優(yōu)點，并且可以根據(jù)公式直接求取x和y方向上的偏導(dǎo)數(shù)fx和fy。結(jié)合式(7)有利于求取坡向，為下一步坡向變率的求取奠定基礎(chǔ)

(7)

取A=3，B=10，C=1/3，m=500，n=500，標準差為1.326 0，均值為0.624 7，格網(wǎng)分辨率為5×5構(gòu)建地形曲面如圖5(a)所示。針對高斯數(shù)學(xué)曲面分別求取偏導(dǎo)fx和fy，再根據(jù)式(7)求得坡向結(jié)果作為原始坡向數(shù)值矩陣。將此坡向數(shù)值矩陣按照2.2節(jié)向量化后的坡向矢量場圖如圖5(b)所示，圖5(c)、(d)、(e)為局部放大圖。

圖5 模擬高斯曲面及曲面坡向Fig.5 Simulated Gaussian surface and its aspect

真實地形所選擇的3個樣區(qū)均位于陜北黃土高原區(qū)，其位置如圖6(a)所示，從南至北依次為宜君樣區(qū)(圖6(d))、吳起樣區(qū)(圖6(c))和綏德樣區(qū)(圖6(b))。宜君樣區(qū)位于陜西省宜君縣城東北部，溝谷溯源侵蝕強烈，重力侵蝕活躍[24]。但是其本身屬于殘塬地貌，因此又具有部分黃土塬地貌特征。吳起樣區(qū)位于陜西省吳起縣北部，區(qū)內(nèi)以黃土墚狀丘陵為主，墚坡上面蝕、細溝和切溝侵蝕處于加速階段，墚地間的沖溝，河溝下切加深[25]。綏德樣區(qū)位于綏德縣無定河中游，區(qū)內(nèi)丘陵起伏，溝壑縱橫，土壤侵蝕極為劇烈[26]。表1為3個樣區(qū)的基本信息。本試驗分別采用3個樣區(qū)5 m分辨率的DEM數(shù)據(jù)作為原始數(shù)據(jù)，設(shè)計相關(guān)算法展開地表坡向變率計算試驗。

圖6 試驗樣區(qū)位置及其地貌暈眩Fig.6 Locations of study areas and the hill shade in sample areas

表1 樣區(qū)概況

3 試驗結(jié)果與分析

3.1 高斯曲面坡向變率計算結(jié)果與分析

針對高斯曲面，將其坡向矩陣分別采用前人提出方法和本文方法求取坡向變率，即采用直接法、正反DEM法、六分法和向量法求取坡向變率，結(jié)果如圖7所示。結(jié)果表明，直接法結(jié)果出現(xiàn)了明顯的“北坡誤差”(圖7(a))，而正反DEM法出現(xiàn)了光滑曲面支離破碎以及鋸齒狀的現(xiàn)象(圖7(b))，六分法得到的結(jié)果曲面相比于正反DEM法的結(jié)果更加光滑(圖7(c))。

本文所提向量法計算得到的結(jié)果與以上3種方法得到的結(jié)果截然不同(圖7(d))。一般而言，地表地形特征是由漸變、突變和不變的地表形態(tài)組成。這種形態(tài)結(jié)構(gòu)也造成了坡向變化在地表空間上的漸變、突變和不變的統(tǒng)一，而漸變和不變的地形表面在地表空間中占據(jù)了大部分區(qū)域，突變的區(qū)域又以地形特征骨架為代表(即山脊，山谷，山頂?shù)?，造成了地形在地表空間上的高頻信息。因此，該特征相對有利于提取坡面突變的區(qū)域。此外，本研究采用各方法的結(jié)果頻率分布曲線對不同算法的特點進行進一步論述。

圖8(a)為直接法、六分法、正反DEM法和向量法計算結(jié)果的坡向變率頻率分布圖。為了更加明顯地展現(xiàn)每種方法各自的趨勢，將4種方法得到的坡向變率結(jié)果進行歸一化處理。結(jié)果可得，直接法、六分法和正反DEM法等標量法的頻率分布結(jié)果類似(幾乎重合，僅少量高值區(qū)存在差異)，與向量法頻率分布結(jié)果差異較大。

圖7 不同方法基于高斯曲面的坡向變率計算結(jié)果Fig.7 SOA results with different methods by using Gaussian surface

3.2 不同樣區(qū)坡向變率計算結(jié)果與分析

由高斯曲面結(jié)果可得，3種標量方法在坡向變率值較小處計算結(jié)果差異較小。因此，針對3個真實樣區(qū)，研究中采用廣泛應(yīng)用的正反DEM法和向量法計算坡向變率，并對所計算結(jié)果進行定量分析。圖9、10、11分別是宜君樣區(qū)、吳起樣區(qū)、綏德樣區(qū)在5 m分辨率DEM下使用這兩種方法的計算結(jié)果。

圖10 吳起樣區(qū)坡向變率不同方法的計算結(jié)果Fig.10 SOA results by different methods in Wuqi area

圖11 綏德樣區(qū)坡向變率不同方法的計算結(jié)果Fig.11 SOA results by different methods in Suide area

由3個樣區(qū)計算結(jié)果可以看出，類似于高斯曲面計算結(jié)果，本文所提向量法得到的結(jié)果明顯不同于傳統(tǒng)的標量計算方式結(jié)果。具體表現(xiàn)為，坡向變率值較高的區(qū)域呈現(xiàn)出較強的空間結(jié)構(gòu)特征，而坡向變化平緩區(qū)域的坡向變率值較小(圖9(b)、圖10(b)、圖11(b))?？梢钥闯觯诓豢紤]坡向存在方向特征時，傳統(tǒng)的基于數(shù)學(xué)標量規(guī)則的坡向變率計算方法在相當程度上夸大了地表面的坡向變化特征，使得地表空間中的不變和漸變的地形特征所表現(xiàn)出來的坡向變化難以得到合理的表達。

在高斯曲面試驗中，高斯曲面是數(shù)學(xué)函數(shù)擬合DEM數(shù)據(jù)。因此，所擬合地形表面相對光滑，坡向變化結(jié)果以緩和為主，較少出現(xiàn)坡向突變的區(qū)域。而實地樣區(qū)的試驗中，樣區(qū)為地表粗糙的黃土高原丘陵溝壑區(qū)，在5 m分辨率的DEM數(shù)據(jù)基礎(chǔ)上，眾多表面細節(jié)信息能夠得到有效的表達，其中包括一定數(shù)量的坡面轉(zhuǎn)折信息。圖8(b)、圖8(c)和圖8(d)分別為宜君樣區(qū)、吳起樣區(qū)和綏德樣區(qū)采用正反DEM法和向量法計算得到的坡向變率頻率分布結(jié)果。圖8(b)為宜君樣區(qū)的頻率分布結(jié)果，可以看到坡向變率值為0區(qū)域頻率很高，這很好地反映了宜君樣區(qū)的殘塬地貌。圖8(d)為綏德樣區(qū)的頻率分布結(jié)果，坡向變率高值區(qū)域頻率明顯高于其他兩個樣區(qū)，這與該樣區(qū)的黃土峁丘陵溝壑地貌類型是相符的。圖8(c)為吳起樣區(qū)的頻率分布結(jié)果，而吳起樣區(qū)本身位于黃土墚狀丘陵溝壑區(qū)，其地形復(fù)雜度高于宜君樣區(qū)而低于綏德樣區(qū)。這樣的復(fù)雜度對比能夠較為容易地從圖8(b)、(c)、(d)的3個頻率分布結(jié)果中得到。

綜合圖8(b)、(c)、(d)、圖9—圖11可知，對于黃土地貌類型區(qū)域，山脊、山谷、山頂?shù)鹊匦翁卣鞑课辉跇訁^(qū)中表現(xiàn)為地形骨架信息，處于量少而又重要的地位。這些地形特征部位坡向變化較大，屬于地表空間地形變化的突變區(qū)域。因此，相應(yīng)的坡向變率計算結(jié)果較大。而黃土地貌類型區(qū)域地形中大部分仍然是坡向轉(zhuǎn)折較小、坡面形態(tài)改變相對平緩的區(qū)域。因此，地形特征仍屬于漸變區(qū)域，相應(yīng)的坡向變率計算結(jié)果較小。向量法得到的坡向變率結(jié)果在相當程度上更好地符合黃土地貌坡向變化的認知，坡向變率低值區(qū)域占大部分，坡向變率高值區(qū)域占小部分。而基于標量運算的正反DEM法恰恰相反，高坡向變率值柵格頻率遠高于低坡向變率值柵格頻率，在相當程度上夸大了地表坡面轉(zhuǎn)變的信息。

通過圖8(b)、(c)、(d)3個樣區(qū)的向量法結(jié)果頻率分布情況，結(jié)合向量法坡向變率的取值范圍，設(shè)置區(qū)間[0,5)、[5,10)、[10,15)、[15,20)、[20,25)、[25,30)、[30,35)、[35,40)、[40,45)、[45,50)、[50,55)對結(jié)果進行統(tǒng)計分析。P(d)為不同區(qū)間內(nèi)坡向變率計算結(jié)果柵格數(shù)占總柵格數(shù)目的百分比，E為區(qū)間內(nèi)坡向變率的均值。3個樣區(qū)的統(tǒng)計結(jié)果如表2所示。可以看出，宜君、吳起和綏德樣區(qū)在區(qū)間[0,15)內(nèi)的柵格占總柵格數(shù)的百分比分別達到80.67%、74.12%和72.56%。該結(jié)果較好地反映了向量法能夠較為真實地計算坡向變化結(jié)果，符合地表空間中地形在不變、漸變與突變的分布關(guān)系。

表2 3個樣區(qū)的向量法坡向變率值在各個區(qū)間的統(tǒng)計

3.3 不同分辨率向量法坡向變率計算結(jié)果分析

本文對3個樣區(qū)的5 m分辨率DEM數(shù)據(jù)，在ArcGIS軟件中，使用最鄰近點法進行了重采樣處理，分別得到10、15、20、25和30 m的DEM數(shù)據(jù)。對這些數(shù)據(jù)使用向量法進行坡向變率的求解，以期分析向量法的穩(wěn)定性。圖12為3個樣區(qū)不同DEM分辨率下向量法計算結(jié)果頻率分布。可以看出，3個樣區(qū)在5 m分辨率下地形信息最為豐富，諸多細節(jié)信息得到保留，地表的不變、漸變特征在較高的分辨率下能夠較好地表達。隨著分辨率的降低，產(chǎn)生削峰填谷現(xiàn)象，部分低等級溝谷和山脊信息得到綜合。但高等級溝谷和山脊信息仍然相對突出，因而表現(xiàn)出坡向變率值較高區(qū)域占比增大。向量法得到的計算結(jié)果在不同樣區(qū)、不同分辨率的高值區(qū)域，其頻率分布曲線形狀基本保持穩(wěn)定。圖13是基于不同DEM分辨率的3個樣區(qū)向量法計算均值(圖13(a))和中位數(shù)(圖13(b))結(jié)果?？梢钥闯?，隨著分辨率的降低，基于向量幾何的坡向變率計算方法在不同分辨率的DEM下表現(xiàn)出了較好的穩(wěn)定性。

圖12 基于不同分辨率的三樣區(qū)DEM數(shù)據(jù)的向量法計算結(jié)果頻率分布Fig.12 Probability distribution of SOA results based different resolutions using vector methods in the three sample areas

圖13 基于不同分辨率的三樣區(qū)DEM數(shù)據(jù)向量法計算結(jié)果統(tǒng)計Fig.13 Statistic result based on different resolution DEMs using vector method in the three sample areas

4 結(jié) 論

本文基于數(shù)學(xué)中的向量幾何方法，考慮到坡向因子的方向?qū)傩?，提出了基于向量方法來求解?shù)字地形分析中的坡向變率地形因子的計算方法。該方法將坡向數(shù)值矩陣進行坐標轉(zhuǎn)換、向量還原、向量運算，最終求解出更加合理的坡向變率結(jié)果。本文所提出的基于向量幾何的坡向變率方法為精準數(shù)字地形分析提供參考。向量幾何在數(shù)字地形分析中的應(yīng)用，對于涉及具有方向?qū)傩缘牡匦我蜃佑嬎愕目茖W(xué)研究具有重要的借鑒意義。從新的數(shù)學(xué)向量角度理解數(shù)字地形分析相關(guān)算法，也是借鑒數(shù)學(xué)向量幾何的方法解決數(shù)字地形分析問題的重要實踐。

本文分別采用高斯合成曲面的離散化數(shù)據(jù)和陜北宜君樣區(qū)、吳起樣區(qū)和綏德樣區(qū)的5 m分辨率DEM數(shù)據(jù)進行試驗，得出以下具體結(jié)論：①基于向量幾何的DEM地表坡向變率計算方法有效地將坡向的方向?qū)傩杂糜谟嬎?，避免了傳統(tǒng)標量方法無視方向?qū)傩栽斐善孪蜃兟视嬎憬Y(jié)果夸大的誤差；②坡向變率計算結(jié)果分布特征更加符合地表空間中坡面轉(zhuǎn)折的漸變、突變以及不變的地貌形態(tài)分布的基本認知，即不變及漸變的地形區(qū)域占據(jù)主要地位；③基于向量幾何的坡向變率計算方法具有較好的穩(wěn)定性，適用于多種分辨率DEM數(shù)據(jù)。

當今數(shù)字地形分析理論中，諸多的地形因子在概念的提出是以向量方式進行描述，帶有重要的方向?qū)傩?，這也符合真實的地貌形態(tài)描述方式。然而，在對地形因子進行計算與表達時，受制于數(shù)據(jù)模型、數(shù)據(jù)結(jié)構(gòu)，無法將其完整含義充分展現(xiàn)，使得絕大多數(shù)地形屬性計算與特征提取采用標量的形式進行求解，造成了數(shù)字地形分析中解譯算法的不確定性問題。本文立足于數(shù)學(xué)向量幾何，從一個小的側(cè)面——坡向變率的求解，將向量的思想引入到數(shù)字地形分析中?？梢钥闯?，向量幾何在數(shù)字地形分析中具有巨大的發(fā)展?jié)摿?。今后的研究將針對地形因子，在?shù)學(xué)上尋找一種更加抽象的表達方式，更加合理的運算方法，來滿足今后數(shù)字地形分析走向地貌研究，走向地形演變的過程與機理研究。