李玉茹，楊勤科，王春梅，吳 江

(西北大學(xué) 城市與環(huán)境學(xué)院，陜西 西安710127)

數(shù)字地形分析研究中，地表粗糙度指真實地表面與理想地表面(如大地水準(zhǔn)面)在垂直方向上的偏離程度[1]，偏差越大，表面越粗糙，反之越平滑。地表粗糙度主要反映地表信息的高頻率(短波長)成分，可以用來表示地球和星球表面的復(fù)雜程度[2-3]。地表粗糙度與諸多地表過程(如坡面徑流、土壤侵蝕、地表災(zāi)害等)有關(guān)[1,4-6]，主要應(yīng)用于地質(zhì)災(zāi)害評價[7]、基本地形類型劃分[8]、水土流失評價[9]及人口分布研究[10]等領(lǐng)域，因而量化地球表面的粗糙程度，是數(shù)字地形分析、地表過程分析和侵蝕地貌學(xué)的重要研究內(nèi)容。

國內(nèi)外學(xué)者提出了一組量化地表粗糙度的算法，Hobson[11]曾將粗糙度算法概括為3種，分別為比表面積、高程統(tǒng)計分布不規(guī)則性和垂直地面法向量的聚集程度。諸多學(xué)者對不同粗糙度算法進(jìn)行了對比研究，Shepard等[12]利用收集的60個數(shù)據(jù)集對比分析了7種粗糙度算法，并對尺度依賴問題進(jìn)行了探討；Grohmann等[6]利用不同空間尺度和不同分辨率的數(shù)據(jù)對6種粗糙度算法進(jìn)行了評價；Berti等[13]以滑坡識別為目的，利用數(shù)學(xué)合成表面和自然表面，對10種粗糙度算法進(jìn)行了對比分析。Lane[4]和Smith[1]先后對地表粗糙度研究進(jìn)行了比較全面的總結(jié)，認(rèn)為有必要對不同研究者從不同角度提出的多種算法進(jìn)行比較分析。目前，對各種粗糙度算法進(jìn)行對比分析，并將計算的粗糙度結(jié)果用于侵蝕地形，特別是地形類型區(qū)分的研究鮮見報道。

本研究以SRTM3高程數(shù)據(jù)為基礎(chǔ)，以區(qū)分地形類型為目標(biāo)，選取文獻(xiàn)報道的11種算法提取研究區(qū)的地表粗糙度，對比分析不同算法的優(yōu)缺點，并探討不同算法區(qū)分研究區(qū)地形類型的能力，旨在為地形類型分類和分區(qū)研究提供參考。

1 研究方法

1.1 研究區(qū)概況與數(shù)據(jù)來源

研究區(qū)覆蓋了陜西、四川、甘肅、湖北、山西和寧夏6個省(區(qū))的部分區(qū)域(29°54′06″N-36°08′48″N，104°54′37″E-110°58′36″E)(圖1)，地理上包括黃土高原南部、秦嶺、巴山和四川盆地東北部，總面積約38萬km2。區(qū)內(nèi)地形類型有平原、丘陵、山地(低山到中高山)，是對比分析各種粗糙度算法，研究各種算法對不同地形類型識別能力的理想研究區(qū)。

本研究使用的數(shù)據(jù)是從CGIAR-CSI網(wǎng)站(http://srtm.csi.cgiar.org/)下載的SRTM3高程數(shù)據(jù)。對下載的數(shù)據(jù)進(jìn)行拼接并將投影轉(zhuǎn)換為6度分帶的高斯投影，中央經(jīng)線為108°E，分辨率設(shè)置為90 m。

圖1 研究區(qū)位置示意圖Fig.1 Location of study area

1.2 地表粗糙度算法

本研究從相關(guān)文獻(xiàn)中選取了以下11種粗糙度計算方法進(jìn)行分析。參考Berti等[13]的研究，考慮到其計算結(jié)果均為粗糙度的量化表達(dá)，所以下文統(tǒng)一將這些計算方法稱為算法。

(1)坡度算法[2]，簡稱“Slope”，計算公式為：

式中：S為坡度，p為中心像元在x方向上的變化率(高程梯度)，q為中心像元在y方向上的變化率(高程梯度)。

(2)局地高差算法[14]，簡稱為“Range”，計算公式如下：

R=Zmax-Zmin。

式中：R為局地高差，Zmax為移動窗口內(nèi)像元最大高程值，Zmin為移動窗口內(nèi)像元最小高程值。

(3)高程均方差算法[13]，簡稱為“RMSH”，計算公式如下：

(4)高程均方根偏差算法[13]，簡稱為“RMSD”，計算公式如下：

式中：RMS_D為高程均方根偏差，Zc為移動窗口中心像元高程值，Zb為移動窗口邊緣像元高程值。

(5)均方根坡度算法[13]，簡稱為“RMSS”，計算公式如下：

式中：RMS_S為均方根坡度,Δxb為中心像元至邊緣像元的距離。

(6)絕對坡度算法[13]，簡稱為“AS”，計算公式為：

式中：A_S為絕對坡度，W為柵格大小。

(7)坡度標(biāo)準(zhǔn)差算法[13]，簡稱為“SDS”，計算公式如下：

(8)面積比算法[6]，簡稱為“AR”，計算公式為：

式中：A為面積比。

(9)二維變異性算法[15]，簡稱為“γ”，計算公式為：

式中：V為二維變異性，(xi,yi)為第i個像元的空間坐標(biāo)，h為滯后距。

(10)矢量離差算法[11]，簡稱為“VD”，計算公式如下：

式中：D為矢量離差，Ai為移動窗口內(nèi)第i個像元的坡向。

(11)余弦特征值算法[16]，簡稱為“DCE”，計算公式如下：

C=ln(E1/E2)。

式中：C為余弦特征值，E1和E2為3×3方向余弦矩陣的歸一化特征值。

除DCE外的其他算法都是數(shù)值越高代表地表越粗糙，對于DCE算法，利用計算出來的最大值減去計算結(jié)果，使其與其他算法保持一致。

1.3 數(shù)據(jù)分析方法

1.3.1 提取結(jié)果比較方法 在Python腳本語言的支持下，用3×3的矩形窗口使用鄰域分析方法提取地表粗糙度。對于計算的各種粗糙度結(jié)果，利用以下方法進(jìn)行分析。

(1)圖像表面格局對比。研究所采用的11種算法計算的粗糙度結(jié)果數(shù)量級不同，因此利用極差標(biāo)準(zhǔn)化方法將提取結(jié)果統(tǒng)一拉伸到0～100。拉伸后的結(jié)果圖像采用相同色帶顯示使其具有可比性。觀察不同算法在不同地形類型區(qū)的圖像表面，直觀對比算法間的不同。

(2)相關(guān)性分析。采用ArcGIS中的Correlation函數(shù)實現(xiàn)，若2種算法提取結(jié)果的空間相關(guān)性越接近1，表明2種算法的相關(guān)性越強，可以相互替代。

(3)統(tǒng)計分布分析。提取每種粗糙度計算結(jié)果(極差標(biāo)準(zhǔn)化后)的統(tǒng)計特征值(平均值、標(biāo)準(zhǔn)差)和頻率分布曲線，從統(tǒng)計分布角度對比算法的異同。

1.3.2 粗糙度計算結(jié)果低通濾波 粗糙度的計算結(jié)果主要反映地表的局地信息(高頻成分)，但也包含了一部分相對低頻的全局信息?？紤]到地形類型是一個比較宏觀的概念，因此為了區(qū)分地形類型，需要對粗糙度計算結(jié)果進(jìn)行低通濾波，以剔除高頻信息，保留相對宏觀的地形信息；并對初步篩選出的粗糙度提取結(jié)果進(jìn)行均值濾波處理。

1.3.3 地形類型區(qū)分能力分析 對濾波后的粗糙度結(jié)果從以下2個方面進(jìn)行分析：

(1)典型斷面分析。在研究區(qū)布設(shè)一條自北向南縱跨黃土丘陵、渭河平原、秦嶺山地的斷面，提取不同粗糙度濾波結(jié)果沿斷面線的形態(tài)分布曲線，討論不同算法區(qū)分地形類型的能力。

(2)頻率曲線分析。單一地形類型條件下，坡度的統(tǒng)計分布呈單峰(正態(tài))分布，而復(fù)合地形區(qū)(包含幾種地形類型)的統(tǒng)計分布，是幾個單一類型之和。據(jù)此，對比濾波后的粗糙度結(jié)果的頻率曲線，根據(jù)頻率曲線的波峰波谷分布分析不同算法區(qū)分地形的能力。

1.4 地表粗糙度計算與分析的技術(shù)路線

技術(shù)路線如圖2所示。

圖2 地表粗糙度計算與分析的技術(shù)路線圖Fig.2 Procedure for surface roughness calculation and analysis

首先下載SRTM3高程數(shù)據(jù)，對高程數(shù)據(jù)進(jìn)行拼接、投影轉(zhuǎn)換等預(yù)處理，通過編程實現(xiàn)對研究區(qū)上述11種粗糙度的計算，然后從空間格局分析、統(tǒng)計分布分析、地形類型區(qū)分能力分析3個角度，對各粗糙度算法進(jìn)行對比，并給出各算法區(qū)分地形類型的適宜性評價。

2 結(jié)果與分析

2.1 不同地表粗糙度算法提取結(jié)果的空間格局分析

2.1.1 圖像表面格局對比 11種地表粗糙度算法提取結(jié)果經(jīng)極差標(biāo)準(zhǔn)化處理后(值域設(shè)為0～100)的結(jié)果如圖3所示。

圖3 不同地表粗糙度算法提取結(jié)果的比較Fig.3 Comparison of results extracted by different surface roughness algorithms

圖3表明：(1)Slope、RMSH、Range、AS、RMSS、RMSD 6種算法提取的結(jié)果高度相似，均表現(xiàn)為秦巴山區(qū)粗糙度值高，黃土丘陵和川中丘陵粗糙度值居中，渭河平原和漢中平原粗糙度值低，從圖3可直觀地區(qū)分山地、丘陵、平原這3種地形；(2)VD、AR和γ 3種算法粗糙度數(shù)值整體較低，高值集中在秦巴山地高程急劇變化的地方，難以區(qū)分丘陵和平原，其中γ算法的結(jié)果對丘陵和平原的區(qū)分能力最差；(3)SDS和DCE 2種算法均難以區(qū)分丘陵和山地，但DCE算法在平原地區(qū)表現(xiàn)出比其他算法更豐富的紋理特征。

2.1.2 空間相關(guān)性 根據(jù)11種算法提取結(jié)果之間的相關(guān)性(表1)，可以將其分為2組，第1組包括Slope、RMSH、Range、AS、RMSS、RMSD、AR、γ 8種算法，這8種算法提取結(jié)果之間的相關(guān)性均高于0.88，相關(guān)性較強，其中Slope、RMSH、Range、AS、RMSS、RMSD 6種算法提取結(jié)果之間的相關(guān)性均高于0.97，相關(guān)性很高；第2組包括SDS、DCE、VD 3種算法，這3種算法與其余8種算法提取結(jié)果之間的相關(guān)性均低于0.73。

表1 不同地表粗糙度算法提取結(jié)果的相關(guān)性Table 1 Correlation of results extracted by different surface roughness algorithms

2.2 不同地表粗糙度算法提取結(jié)果的統(tǒng)計分布分析

統(tǒng)計分布是對地理數(shù)據(jù)表面值整體特征進(jìn)行分析的最常用方法，本研究從統(tǒng)計特征值和頻率分布曲線兩方面對各算法進(jìn)行分析。

2.2.1 統(tǒng)計特征值分析 極差標(biāo)準(zhǔn)化(值域設(shè)為0～100)后，不同算法提取粗糙度的統(tǒng)計特征值如圖4所示。

圖4 不同地表粗糙度算法提取結(jié)果的統(tǒng)計特征值Fig.4 Arithmetic mean and STD of results extracted by different surface roughness algorithms

由圖4可見：(1) Slope、RMSH、Range、RMSD、RMSS、AS 6種算法的標(biāo)準(zhǔn)差與平均值相近，與2.1節(jié)中這6種算法提取結(jié)果的空間格局相似、相關(guān)性極高的結(jié)果一致，且這6種算法的標(biāo)準(zhǔn)差高于其他算法，說明其粗糙度結(jié)果的離散程度高于其他算法，表達(dá)了更豐富的表面信息；(2) AR、VD、γ 3種算法的平均值和標(biāo)準(zhǔn)差都很低，與2.1節(jié)中這3種算法結(jié)果數(shù)值整體偏低、層次不分明的圖像特征一致；(3) DCE算法的平均值最高，但標(biāo)準(zhǔn)差居中，這與其提取結(jié)果在山地和丘陵區(qū)粗糙度值均較高的特點一致;(4) SDS算法平均值與標(biāo)準(zhǔn)差均居中，與其提取結(jié)果數(shù)值整體不高，且難以區(qū)分山地和丘陵的特征一致。

2.2.2 頻率分布曲線 各粗糙度結(jié)果極差標(biāo)準(zhǔn)化(值域設(shè)為0～100)后的頻率分布曲線如圖5所示。由圖5可見：(1) Slope、RMSH、Range、RMSD、RMSS、AS 6種算法的頻率曲線幾乎重疊，分別在粗糙度值為0.2和5.5附近有2個波峰，SDS算法的頻率曲線雖然也有2個波峰，但第2個波峰位置在粗糙度值為3.2附近；(2) DCE算法的頻率曲線比較特殊，只有1個明顯的波峰，且波峰對應(yīng)的粗糙度集中在高值(60～80)；(3) VD、AR、γ 3種算法沒有明顯波峰，且粗糙度集中在低值。不同算法頻率分布曲線的變化與表面值的分布特征和統(tǒng)計特征值變化一致。

圖5 不同地表粗糙度算法提取結(jié)果的頻率分布曲線Fig.5 Distribution of results extracted by different surface roughness algorithms

2.3 不同地表粗糙度算法的地形類型區(qū)分能力分析

由不同算法空間格局和統(tǒng)計分布特征分析結(jié)果可以看出，Slope、RMSH、Range、RMSD、RMSS、AS 6種算法有較強的相似性，且對地形特征的表達(dá)都較好；DCE算法具有一定特殊性，與其他算法相比，其在平原區(qū)有更好的紋理特征。綜合以上分析，且考慮到Slope和Range算法成熟且應(yīng)用廣泛，故選用Slope、Range和DCE 3種算法，分析其區(qū)分地形類型的能力。

2.3.1 濾波窗口選取 不同矩形窗口均值濾波損失的高頻信息統(tǒng)計結(jié)果如圖6所示。

圖6 不同矩形窗口均值濾波損失的高頻信息統(tǒng)計特征值Fig.6 Maximum,minimum and local variances of high frequency information lost by mean filter of different rectangular windows

由圖6可見，當(dāng)矩形窗口大小達(dá)到13×13時，高頻信息的最小值、最大值和局地方差基本趨于穩(wěn)定，說明濾掉的高頻信息不再增加，根據(jù)這一結(jié)果，本研究最終選擇13×13的矩形窗口對粗糙度結(jié)果進(jìn)行均值濾波。

2.3.2 典型斷面分析 研究布設(shè)的斷面線如圖7所示，Slope、Range、DCE 3種算法濾波后的結(jié)果(值域拉伸到0～100)沿斷面線的形態(tài)分布如圖8所示，這3種算法濾波后的結(jié)果在不同地形區(qū)的分布情況見圖9。

圖7 斷面線位置示意Fig.7 Location of the profile

圖8 不同算法地表粗糙度濾波后結(jié)果的斷面Fig.8 Profiles of different algorithms based roughness with low pass filtering

箱線圖從下至上依次為頻率10%,25%,50%,75%,90%的個數(shù)，十字絲代表平均值The boxplots from bottom to top refer to the frequency of 10%,25%,50%,75%,90% and the crosshairs refer to mean values

由圖8和圖9可以看出，Slope和Range濾波后粗糙度沿斷面線的形態(tài)分布幾乎重疊，3種地形對應(yīng)的縱坐標(biāo)分別集中在3個層次上；這2種算法在3種地形區(qū)至少80%的粗糙度值完全不相交，可以明顯區(qū)分平原、丘陵、山地。DCE濾波后粗糙度的形態(tài)分布與Slope、Range不同，在平原區(qū)表現(xiàn)出更復(fù)雜的地表粗糙度特征，縮小了平原與丘陵的差距，丘陵與山地粗糙度值重疊度高，混淆程度大；從箱線圖中也可看出，DCE方法在3種地形區(qū)粗糙度值相交程度大，說明DCE方法區(qū)分地形的能力較弱。

2.3.3 濾波后粗糙度的頻率分布曲線 不同地形類型的粗糙度頻率曲線呈現(xiàn)不同的特征，其峰值和值域范圍均有較大差異。當(dāng)研究區(qū)存在多種地形類型時，粗糙度頻率曲線呈多峰分布。通過分析粗糙度頻率曲線的波峰-波谷序列，可辨識粗糙度算法對各種地形類型的區(qū)分能力。Slope、Range和DCE 3種算法濾波后的頻率曲線如圖10所示。由圖10可見，Slope和Range算法濾波后粗糙度結(jié)果的頻率曲線相似，都有4個波峰-波谷序列，說明用這2種算法提取的粗糙度可以辨識平原、臺地、丘陵和山地4種地形類型；DCE只有2個明顯波峰，說明其只能辨識部分地形類型，這與2.3.2節(jié)的結(jié)果一致。

圖10 不同地表粗糙度算法濾波后結(jié)果的頻率分布曲線Fig.10 Results after filtering by different algorithms

3 討 論

根據(jù)地表粗糙度算法的地理學(xué)原理，將其分為3類，包括基于高程變異性的算法、基于坡度變異性的算法和基于地表法向量的算法。

基于高程變異性的算法包括 Slope、RMSH、Range、AS、RMSS、RMSD、AR和γ。本研究中，Slope、RMSH、Range、AS、RMSS、RMSD 6種算法之間的相關(guān)性很高，在分辨率和計算窗口固定的前提下，其基本地學(xué)含義都是高差。Range、RMSH、RMSD 3種算法都是直接對窗口內(nèi)像元高程的變異性進(jìn)行分析得到的，其中Range是計算窗口內(nèi)的最大高差，RMSH是計算窗口內(nèi)的高程均方差，RMSD與RMSH相似，但只計算中心像元與邊緣像元高差的均方根。RMSS和AS實質(zhì)上是Slope的變形，都是根據(jù)高差與水平距離之比來計算的。因此，Slope、RMSH、Range、AS、RMSS、RMSD 6種算法在提取地表粗糙度時表現(xiàn)出相似的特征。AR、γ與坡度的相關(guān)性也很高，但提取效果卻不如上述6種算法。AR基于比表面積(斜面面積與投影面積之比)，具體計算時AR是坡度的余弦函數(shù)，會壓縮低坡度值，拉伸高坡度值，導(dǎo)致突出山地，模糊丘陵和平原。γ描述了移動窗口中心像元與鄰域像元的相關(guān)性，地表越粗糙，地表高程值的空間差異性越大，對應(yīng)的空間依賴范圍越短，變異性越高。固定窗口下，AR和γ算法實質(zhì)也是計算局部地表高程的高差，因此與坡度的相關(guān)性也較高。但γ計算了高差的平方，會放大地表高程變化大的地域如山地，從而導(dǎo)致提取結(jié)果只突出起伏較大的山區(qū)和丘陵。

基于坡度變異性的算法只有SDS，其與坡度的變化率高度相關(guān)。SDS實質(zhì)為窗口像元坡度的標(biāo)準(zhǔn)差，表達(dá)了地表坡度的變異程度，這與Evans[17]、李新艷等[18]的認(rèn)識一致。這一原因?qū)е耂DS算法不易區(qū)分坡度變率相似的丘陵和山地，所以在地形類型分區(qū)時，不是重點考慮的算法。

基于地表法向量的算法包括DCE和VD。這2種算法用垂直于地面的法向量在空間上的聚集程度刻畫地表粗糙度，法向量越聚集，地表越平滑；越分散，地表越粗糙。這2種算法在一定程度上表達(dá)了地表坡度和坡向的局地變異，因而不適合對宏觀地形類型進(jìn)行劃分。

本研究以地形類型劃分為目的，對比分析不同地表粗糙度算法的特點，所以在認(rèn)識上與已有的研究結(jié)論不完全相同。在地質(zhì)災(zāi)害應(yīng)用研究中，Mckean等[7]基于高分辨率地形數(shù)據(jù)，通過計算地表粗糙度來識別滑坡，認(rèn)為DCE算法在局地范圍內(nèi)用于滑坡范圍的識別，是較為適宜的。本研究則認(rèn)為，DCE算法在區(qū)域尺度上進(jìn)行地形類型的劃分并不適宜。在區(qū)域尺度的研究中，Grohmann等[6]基于分辨率5 m DEM的分析認(rèn)為，VD算法很難表達(dá)區(qū)域性地形起伏，本研究也認(rèn)為VD算法未能有效劃分區(qū)域尺度地形區(qū)。地貌類型制圖研究中，常用起伏度(Range)來劃分基本地形形態(tài)(平原、臺地、丘陵和山地)[8,19-20]。本研究通過分析不同算法區(qū)分地形類型的能力，進(jìn)一步認(rèn)為Range是區(qū)域尺度地形類型劃分的較優(yōu)算法。

本研究討論的11種算法中，除坡度算法以外，其余算法都與計算時移動窗口的尺寸有關(guān)。窗口尺寸是一個不可回避的問題，有研究對此進(jìn)行過探討，包括尋找最佳窗口[9,21]、避免窗口尺寸影響[22]等。但該問題還在探討中，窗口尺寸對粗糙度的影響有待進(jìn)一步研究。

另外，近年來更復(fù)雜的基于頻譜分析的傅里葉算法[23]和小波算法[24]也被用于粗糙度的提取。Berti等[13]認(rèn)為，簡單算法更能用于檢測滑坡信息，因此本研究中未考慮這些復(fù)雜計算方法，這些算法也有必要進(jìn)一步分析探討。

4 結(jié) 論

本研究以辨識和區(qū)分地形類型為目標(biāo)，探討了坡度(Slope)、局地高差(Range)、高程均方差(RMSH)、高程均方根偏差(RMSD)、均方根坡度(RMSS)、絕對坡度(AS)、坡度標(biāo)準(zhǔn)差(SDS)、面積比(AR)、矢量離差(VD)、余弦特征值(DCE)、二維變異性(γ)這11種地表粗糙度算法的特點及其區(qū)分地形類型的能力，可得出以下結(jié)論：

(1)Slope、Range、RMSH、RMSD、RMSS、AS 6種算法的空間格局特征和統(tǒng)計分布特征都極其相似，其中具有代表性的Slope和Range算法經(jīng)過低通濾波消除局部高頻信息后，可有效區(qū)分研究區(qū)地形類型?？梢哉J(rèn)為以Slope和Range為代表的這6種算法是表達(dá)地表粗糙程度和區(qū)分地形類型的較佳算法。

(2)AR和γ算法提取結(jié)果難以區(qū)分平原和丘陵，雖與坡度相關(guān)性較高，但因為經(jīng)過非線性變化，反而不能有效反映地表粗糙度的空間變異。SDS算法表達(dá)了坡度空間變異，提取結(jié)果難以區(qū)分坡度變率相近的丘陵和山地。DCE算法提取結(jié)果雖在平原地區(qū)有更好的紋理特征，但同樣難以區(qū)分丘陵和山地。VD算法提取結(jié)果在空間格局特征和統(tǒng)計分布特征兩方面的表現(xiàn)最差?？梢哉J(rèn)為AR、γ、SDS、DCE、VD 5種算法都難以表達(dá)研究區(qū)地表粗糙程度的宏觀變異，不能有效區(qū)分研究區(qū)的地形類型。