青藏高原地區3種全球DEM精度對不同地形因子的響應

高志遠, 謝元禮,3, 王寧練,3, 蔣廣鑫, 周 鵬

(1.西北大學 城市與環境學院, 陜西 西安 710127; 2.陜西省地表系統與環境承載力重點實驗室,陜西 西安 710127; 3.西北大學 地表系統與災害研究院, 陜西 西安 710127)

在地理學研究中，數字高程模型(digital elevation model, DEM)已成為一個重要且可靠的分析數據源，在土壤侵蝕、水土流失、流域分析、地貌分析以及地形因子提取方面有著十分廣泛的應用。針對全球尺度DEM產品的精度評估工作已有不少進展。Rodriguez等人[1]利用絕對地理誤差、絕對高程誤差和相對高程誤差等精度指標對全球6大洲的SRTM DEM數據進行精度評估；Shortridge等人[2]針對全美格局的景觀結構對SRTM DEM精度影響進行了大范圍的評估分析。但是對于特定地區DEM精度評價可能面臨以下問題：第一，全球大尺度DEM精度評價結果不一定適用于特定區域的數字高程模型；第二，DEM精度評價主要依賴于地面測量數據如GPS測量作為驗證標準[3-4]，而地面測量數據對于特殊地理環境地區通常較難獲取。高精度星載雷達(ICESAT/GLAS)數據的水平精度為±20 cm，高程精度為±18 cm，已有不少學者以ICESAT/GLAS高程數據為參考進行了相關的精度分析。En?le等人[5]利用GLAH14數據進行了樹冠頂層、地面高程和植被高度的精度評價。Yue[6]利用ICESAT/GLAS高程數據對ASTER GDEM數據進行改正并與SRTM DEM數據的進行無縫融合。Liu等人[7]利用經過ICESAT/GLAS數據糾正后的ASTER數據計算蘭伯特冰川的變化率，結果表明其速率符合當期冰川的物質平衡。但是全球DEM在青藏高原地區的精度情況及其地形影響因素研究的比較少。為此，本文中選取3種的全球免費中分辨率DEM數據SRTM DEM(shuttle radar topography mission digital elevation model)，ASTER GDEM(the advanced spaceborne thermal emission and reflection radiometer global digital elevation model)以及HydroSHEDS，以星載激光雷達ICESAT/GLAS數據中的GLAH14數據為高程參考標準，針對青藏高原地區具體情況，對SRTM DEM，ASTER GDEM以及HydroSHEDS DEM這3種DEM進行整體的精度評價，并分別從坡度、坡向以及地形粗糙度等地形因子進行精度變化趨勢及規律的探究。

1 研究區概況與方法

1.1 研究區概況

青藏高原是世界上面積最大的高原，同時也是世界上平均海拔最高的高原，有“世界屋脊”和“第三極”之稱[8]，青藏高原總面積為2.60×106km2，平均海拔4 500 m，全球海拔超過8 000 m的14座山峰全部分布在該地區[9]。青藏高原是中國3個凍土主要聚集區域之一，分布有大量季節性凍土和部分多年凍土[10]，因為其凍土面積廣泛，地表覆蓋特征受植被季節性變化影響小。萬杰等人[11]提出時間間隔為4 a的青藏高原ICESAT/GLAS數據和SRTM DEM數據的高程線性擬合的決定系數高達0.999 8，因此在青藏高原地區利用多年ICESAT/GLAS數據進行DEM精度評價十分合適。本文研究區為青藏高原非冰雪覆蓋區域。

1.2 研究數據

本次研究使用的數據主要包括DEM數據和ICESAT/GALS數據。DEM數據包括SRTM DEM，ASTER GDEM和HydroSHEDS。SRTM DEM數據選用的是SRTM V4版本，SRTM V4版本數據是針對SRTM V2版本的缺失數據進行插值補充；ASTER GDEM為對地觀測衛星Terra的近紅外波段獲取的立體像對生成的DEM；HydroSHEDS DEM是在3弧秒的SRTM DEM數據基礎上，采用新改進和新開發的算法所生產的數據。3類DEM數據對比如表1所示。ICESAT/GLAS數據采用的是GLAH14數據，GLAH14數據是由GLAH05數據和GLAH06數據再生產的陸地表面二級測高數據。GLA14數據為散點記錄方式，相鄰點的距離為172 m，與目前常見的數字高程模型相比GLAH14數據有更高的精度，因此可以當成DEM的驗證標準來進行精度驗證。ICESAT/GLAS工作周期為2003—2009年，共18個工作時間段，數據來源于國家冰雪數據中心(National Snow and Ice Data Center， http:∥nsidc.org/)，試驗中在GLAH14的18期數據序列中選取300 000余個數據點進行分析。

1.3 數據預處理

由于3種DEM數據的參考橢球體與GLAH14數據不同，要對GLAH14數據進行預處理工作。預處理工作包括高程基準處理、參考橢球體處理和粗差剔除3部分。

1.3.1 高程基準處理 在GLAH14的數據集中，直接利用HDFView工具提取出大地水準面差距的數值t，再利用公式(1)進行計算，把橢球體高程轉化為正高(地面點沿過此點的重力線到大地水準面的距離)。

HT/P=h-t

(1)

式中：HT/P——TOPHEX/Poseidon橢球體的正高；t——水準面差距。

1.3.2 參考橢球體處理 3種DEM的參考橢球體都是WGS84橢球體，而GLAH14的參考橢球體是TOPHEX/Poseidon橢球體(T/P橢球體)。橢球體之間的差異主要來源于水平差異和垂直差異兩方面，兩者在水平方向上的差異只有厘米極[12]，因此這里主要針對高程差異進行轉換。表2為TOPHEX/Poseidon橢球體和WGS84橢球體參數對比，經計算可知，兩種橢球體的高程差異約在70 ～71 cm之間，一般取70 cm[13]。最后再利用公式(2)進行計算，求出WGS84的正高高程。

表2 T/P橢球體和WGS84橢球體參數對比

HWGS84=HT/P-0.70=h-t-0.70

(2)

式中：HWGS84——WGS84橢球體的正高；HT/P——TOPHEX/Poseidon橢球體的正高。

1.3.3 粗差剔除 先剔除NODATA的數據點(數據異常值點)，然后計算GLAH14數據點上對應的3種DEM的差值，高程差超過標準差3倍的數據點不參與統計。

1.4 精度指標

本文利用坡度、坡向以及地形粗糙度等地形因子，結合不同的精度指標，綜合對3種DEM進行精度評價以及地形因子的響應分析。本文利用高程的誤差d，平均誤差Mean及中誤差RSME作為評價精度的指標。其公式如式(3)—(5)所示。

d=hDEM-hICESAT/GLAS

(3)

(4)

(5)

其中高程誤差d表示了數據集中每個DEM數據點高程與GLAH14數據之間的差值；平均誤差 表示數據集中DEM數據點高程與GLAH14數據的高程誤差的平均值；中誤差 代表了數據集中DEM高程數據與相應GLAH14數據之間的離散程度，是評價精度的直接指標。綜合研究3種DEM在不同的精度指標與各個地形因子之間的關系及變化規律。

針對DEM分形維數的計算方法，本文采用表面積—尺度法，表面積—尺度法也稱作投影覆蓋法，此方法由Clarke等人[14]于1986年首次提出，而后謝和平、周銀軍等人[15-16]又在此方法的基礎之上做出了改進。首先，需要對原始的DEM柵格進行重采樣，以分出不同的空間分辨率尺度Ri，并且計算每一種空間分辨率尺度下DEM的表面積Si。而后，計算各個尺度下DEM分辨率尺度Ri的對數值以及其對應表面積Si的對數值，并將其進行線性擬合，再在雙對數坐標軸上找尋一段擬合較好的區域作為無標度區。雖然分形維數具有尺度不變性，但當尺度大到一定值時，粗糙表面根本不表現出分形性質[17]。因此只有當觀測尺度介于一定范圍內時，粗糙表面才能出現分形性質，因此需要尋找一個合適的尺度范圍來計算分形維數D。

1.5 試驗方案

利用以上幾種精度評價指標，先計算總體精度進行統計，然后結合坡度、坡向以及地表粗糙度等地形因子進行3種DEM精度的分析。在坡度因子分析中，對坡度進行離散化分級，級差根據地形特點選取1°～5°，在每個坡度級內統計3種DEM的平均誤差和中誤差，分析坡度對精度的影響。在坡向因子分析中，將坡向按照45°范圍分為8個坡向區域，在每個坡向區域內分別統計3種DEM的平均誤差和中誤差，分析坡向對精度的影響。分形維數的計算是在青藏高原研究區選取600余個ICESAT/GLAS數據點較多分布的區域，每個區域大小為60 km×60 km(圖1)，因為在60 km的空間尺度下，研究區地形有較為明顯的分維特征。為了保證精度，地形因子采用SRTM DEM進行計算，將每個區域的DEM柵格進行重采樣，采樣分辨率范圍為100～600 m，以10 m作為間隔，最后利用求得的分形維數D與每個區域的中誤差擬合，分析地形破碎度對精度的影響。

圖1 青藏高原分形維數試驗分區

2 結果及分析

經過統計計算，青藏高原地區SRTM DEM平均誤差為0.49 m，中誤差為14.29 m。ASTER GDEM平均誤差為1.71 m，中誤差為18.81 m。Hydro SHEDS DEM平均誤差為4.17 m，中誤差為31.08 m。其中，SRTM DEM數據精度與萬杰等人[11]在青藏高原地區研究得到的16.65 m高程中誤差相當，ASTER GDEM數據精度與杜小平等人[13]在中國西南地區高海拔山區得到的20.42 m高程中誤差相接近。很明顯地，在青藏高原地區，SRTM DEM精度最高，HydroSHEDS DEM精度最低。

2.1 坡度對DEM精度的影響

將DEM與GLAH14高程的高程誤差d(hDEM-hICESAT/GLAS)按照點所在的坡度級進行統計，坡度級差設定為：<1°，1°～5°，5°～10°，10°～15°，15°～20°，20°～25°，25°～30°，30°～35°，35°～40°，40°～45°，45°～50°以及>50°分別求出在各個坡度級中3種DEM的平均誤差Mean和中誤差RSME，計算結果詳見表3。

表3 各坡度區間3類DEM精度指標

3種DEM平均誤差隨坡度的變化及中誤差隨坡度的變化如圖2所示。SRTM DEM的平均誤差在坡度為40°～45°時較大，HydroSHEDS DEM平均誤差在坡度為0°～10°及35°～45°出現兩個峰值，而ASTER GDEM的平均誤差隨坡度增大而增大。從中誤差分析來看，SRTM DEM數據高程精度在各個坡度上均優于ASTER GDEM；HydroSHEDS DEM中誤差隨坡度變化規律呈明顯的對數關系，與SRTM DEM和ASTER GDEM的規律有很大的不同。在地形較為平坦地區(坡度小于5°)，HydroSHEDS DEM精度高于ASTER GDEM精度，但在坡度較大地區(坡度大于5°)，則ASTER GDEM精度高于HydroSHEDS DEM精度。

圖2 3類DEM中誤差和平均誤差與坡度的關系

2.2 坡向對DEM精度影響

將坡向劃分為北方向(337.5°～22.5°)、東北方向(22.5°～67.5°)、東方向(67.5°～112.5°)、東南方向(112.5°～157.5°)、南方向(157.5°～202.5°)、西南方向(202.5°～247.5°)西方向(247.5°～292.5°)以及西北方向(292.5°～337.5°)8個坡向級。圖3為3種DEM高程誤差隨坡向分異圖，圖中上下邊緣代表了除數據異常值(大于或小于上下四分位數1.5倍四分位數差的數據)外的最大和最小值。3種DEM的高程誤差分布均具有不同程°的坡向分異性特征，其中SRTM DEM的分異幅°最小，不同坡向的高程誤差數據中位數接近，且均大于0，高程誤差隨坡向的分異特征不明顯，南方向和東南方向SRTM DEM的高程誤差相對較大。相較于SRTM DEM，ASTER GDEM和HydroSHEDS DEM數據高程誤差隨坡向的分異性強，規律顯著，且高程誤差有明顯的正負分布性。

其中ASTER GDEM高程誤差正值最大值分布在東北方向和西南方向，而高程誤差負值在東南方向最大；HydroSHEDS DEM的高程誤差則呈現明顯的“東北—西南對稱”分布，坡向為337.5°～157.5°時，高程誤差總體傾向為正值，且誤差在東北方向和北方向值最大，當坡向為157.5°～337.5°時，高程誤差總體傾向為負值，誤差在西南方向和南方向值最大。可見不同DEM測量值在不同坡向上的正負分布有一定差異性。

SRTM DEM ASTER GDEM HydroSHEDS DEM

圖3 青藏高原3種DEM高程誤差隨坡向分異

為了進一步研究這種差異規律特征，本文分別對測量值大于正2倍中誤差和小于負2倍中誤差數據點(下文中稱為正負2倍中誤差數據點，即di>2RSME和di<-2RSME)按照其在某一坡向數據點所占百分比進行計算統計，嘗試探究DEM測量正偏離值、負偏離值隨坡向的分異特征，結果如圖4所示。

3種DEM測量偏離值的分布呈現不同的特征，且程度不同。SRTM DEM正負2倍中誤差數據點主要呈“南—北”分布，正2倍中誤差數據點主要集中在南方向、西南方向和東南方向，其中南方向占比達3.23%，負2倍中誤差數據點則主要出現在西北方向上，占比為3.07%，ASTER GDEM正負2倍中誤差數據點呈現“西北—東南”分布，正2倍中誤差數據點主要出現在西方向，占比為1.46%，負2倍中誤差數據點主要集中在南方向，占比為0.89%。HydroSHEDS DEM正負2倍中誤差數據點呈現“東北—西南”分布，正2倍差數據點峰值分布在東方向上，占比為21.52%，負2倍中誤差數據點峰值位于西南方向，占比為8.69%。

研究表明，SRTM DEM測量正偏離值主要分布在南方向，而負偏離值主要分布在北方向。ASTER GDEM測量正偏離值主要分布在西北方向，負偏離值分布在東南方向，HydroSHEDS DEM測量正偏離值分布在東北方向，負偏離值分布在西南方向。3種DEM測量偏離值的離散程度不盡相同，HydroSHEDS DEM的偏離程度最大，ASTER GDEM的偏離程度最小。3種DEM出現高程測量偏離值隨坡向的分異特征的原因可能與衛星傳感器在上升軌道和下降軌道的航向[2]以及SRTM傳感器雷達與地表的入射角度有關。

圖4 3類DEM各坡向正負2倍中誤差數據點占比

Zhang[18]等人在研究地表覆蓋對SRTM DEM精度的影響時將坡度進行分級固定，以排除坡度因子的干擾，本文利用這一思想進行坡度、坡向的控制變量分析。不同DEM的中誤差隨坡向變化規律不同，且隨著坡度的增加，中誤差隨坡向變化幅度增大，呈現更加明顯的坡向分異特征。

如圖5所示，當坡向固定時，3種DEM中誤差隨著坡度的增加，呈現增加趨勢，其中SRTM DEM和ASTER GDEM在坡度大于30°時中誤差迅速增加，HydroSHEDS DEM在坡度大于8°時，精度迅速變差。而當坡度處于不同的特定區間時，3種DEM精度的坡向分異規律也不同，其中SRTM DEM在各個坡度區間內中誤差隨坡向的分異規律均不明顯，ASTER GDEM在坡度小于16°時，規律不顯著，而當坡度大于16°時，坡向分異明顯，且隨著坡度的增加，規律愈發突出，其中最大中誤差集中在西南方向，而在西北方向和北方向上精度最好，受到坡度的影響較小。HydroSHEDS DEM在坡度小于8°時，各個坡向上的精度大致相當，且中誤差均在30 m以內，當坡度大于8°時，HydroSHEDS DEM中誤差在不同坡向上呈現明顯的“雙峰”分布，即在東北方向和西南方向中誤差往往較大，而在東南方向，精度相對較好。由以上分析可以發現，HydroSHEDS DEM受到地形坡度的波動精度變異性更明顯，并且相對SRTM DEM和ASTER GDEM其在更小的的坡度區間內就出現了明顯的坡向分異性。

SRTM DEM ASTER GDEM HydroSHEDS DEM

圖5 3類DEM各坡向中誤差分市

2.3 DEM精度與地形粗糙度的關系

按照分形維數的計算步驟進行逐個樣方的計算，每個樣方中樣點的數量從168到756，656個樣方的分形維數D計算統計結果如表4所示。

表4 分形維數D統計結果

由圖6可以看出，樣方分形維數的分布主要集中在2.00～2.06之間，當分形維數大于2.06時，樣方數迅速減少，而在2.00～2.01區間內的樣方數最多，為262個，達到總數的40%。

通過對3種DEM中誤差與分維數擬合，發現它們均呈現二次多項式關系，如圖7所示。SRTM DEM精度受到地形破碎度的影響較小，ASTER GDEM次之，HydroSHEDS DEM精度隨地形破碎度影響最大。當分形維數較大時(D>2.04)時，SRTM DEM和ASTER GDEM中誤差與分維數的相關性變弱，而HydroSHEDS DEM中誤差與分維數一貫的相關性，可見HydroSHEDS DEM精度對地形破碎度響應較為敏感。

圖6 分形維數數頻數

SRTM DEM ASTER GDEM HydroSHEDS DEM

圖7 3種DEM分形維數和中誤差的關系

3 結 論

(1) 在青藏高原地區，SRTM DEM的精度最高，中誤差為14.29 m，ASTER GDEM的精度次之，中誤差為18.81 m，HydroSHEDS DEM精度最差，中誤差為31.08 m。

(2) SRTM DEM精度隨坡度變化趨勢與ASTER GDEM趨勢相似，但是在各坡度區間內略優于ASTER GDEM。HydroSHEDS DEM中誤差隨坡度變化呈明顯的對數關系，在地形較為平坦地區(坡度小于5°)，HydroSHEDS DEM精度高于ASTER GDEM精度；但在坡度較大地區(坡度大于5°)，ASTER GDEM精度高于HydroSHEDS DEM精度。

(3) SRTM DEM誤差正偏離主要位于南坡向，而負偏離主要分布在北坡向；ASTER GDEM誤差正偏離主要分布在西北坡向，負偏離分布在東南坡向；HydroSHEDS DEM誤差正偏離主要分布在東北坡向，負偏離主要分布在西南坡向。HydroSHEDS DEM測量偏離值的離散程度最大，ASTER GDEM離散程度最小。

(4) SRTM DEM中誤差在不同坡度、各個坡向上分異幅度較小，具有一定的同質性；ASTER GDEM在坡度大于16°時分異性較強，中誤差峰值出現在西南坡向；HydroSHEDS DEM在坡度大于8°時分異性較強，最大中誤差分布在西南和東北坡向。

(5) 分形維數D和3種DEM中誤差擬合均呈較為明顯的二次多項式關系，分形維數越小，分形維數對DEM精度的影響越大。其中分形維數對HydroSHEDS DEM精度影響最大，對SRTM DEM精度影響最小。