萬 杰，廖靜娟，許 濤，沈國狀

(1.中國科學院遙感與數(shù)字地球研究所數(shù)字地球科學重點實驗室，北京 100094;2.中國科學院大學，北京 100049)

0 引言

數(shù)字高程模型(DEM)在地球科學及其他相關學科(如重力場建模、水文學、地形圖制圖學和航空圖像正射糾正等)研究中有廣泛的應用[1]。航天飛機雷達測圖計劃(shuttle Radar topography mission，SRTM)是第一個公開發(fā)布全球高精度數(shù)字高程模型的遙感平臺。國內外研究采用不同的參照數(shù)據(jù)(如 GPS[2－3])，參考 DEM[4]和高度計數(shù)據(jù)[5－6]評估SRTM的精度。由于SRTM的版本、高度計數(shù)據(jù)的版本、實驗樣區(qū)和參照數(shù)據(jù)的變化，上述研究及其結論各有差異。Carabajal等[7]通過對美國西部、南美洲(亞馬孫流域)、非洲、亞洲和澳大利亞部分區(qū)域的SRTM的分析，發(fā)現(xiàn)不同區(qū)域SRTM的高程誤差呈現(xiàn)出一定規(guī)律的變化。Rodriguez等[3]通過比較SRTM和車載KGPS導航儀數(shù)據(jù)，得到SRTM的水平和垂直精度分別為11.9 m和9 m。此外，植被覆蓋度和地形起伏的增大也會明顯降低SRTM的數(shù)據(jù)質量[6]。

美國國家航空航天管理局(NASA)的冰、云和陸地高程衛(wèi)星(ICESat)攜載的激光高度計(geoscience laser altimeter system，GLAS)數(shù)據(jù)在垂直方向上的測量精度達到了cm級，遠遠超過了SRTM數(shù)據(jù)的精度[2]，因此可以利用ICESat/GLAS激光高度計數(shù)據(jù)(以下簡稱ICESat高度計數(shù)據(jù))評估SRTM高程數(shù)據(jù)的精度。此外，相對于通常處于視野開闊和平坦區(qū)域的GPS數(shù)據(jù)和控制點不能遍及整個區(qū)域的局限性，ICESat腳印點(footprint point)可分布在整個實驗樣區(qū)。因此，本文采用精度更高的ICESat高度計數(shù)據(jù)分析了SRTM的系統(tǒng)誤差和高程精度，并探討了SRTM高程數(shù)據(jù)的誤差特征與地形因子(坡度和坡向)間的關系。

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

1.1 研究區(qū)概況

本文的研究區(qū)域位于中國的青藏高原。青藏高原東西方向長約2 945 km，南北方向長約1 532 km，覆蓋了西藏自治區(qū)全境和云南、四川、甘肅、青海和新疆的部分區(qū)域，總面積約為2 570 000 km2。除兩極區(qū)域之外，青藏高原是地球上最寒冷的區(qū)域，也被稱為“世界第三極”。青藏高原是一個地形條件極為復雜的區(qū)域，有著眾多的山脈及占地面積大約49 873 km2的36 800 條冰川和 1 500 多個湖泊[8]。青藏高原的高程變化在60～8 792 m之間，平均海拔在4 000 m以上(圖1)。

圖1 研究區(qū)DEMFig.1 DEM of study area

1.2 數(shù)據(jù)源

本文實驗使用的數(shù)據(jù)包括SRTM高程數(shù)據(jù)、ICESat高度計數(shù)據(jù)和 MODIS地表類型產品(MCD12Q1)。

1)SRTM高程數(shù)據(jù)。SRTM于2000年2月發(fā)射，獲取了覆蓋地球N 60°～S 57°之間地表的數(shù)字高程。SRTM免費提供30 m分辨率(美國境內)和90 m分辨率(其他地區(qū))的高程數(shù)據(jù)產品。SRTM DEM數(shù)據(jù)的高程基準是EGM96的大地水準面，平面基準是WGS84[6]。SRTM數(shù)據(jù)設計的水平和垂直精度分別高于20 m和16 m[9]。本文使用的SRTM高程數(shù)據(jù)是由國際農業(yè)研究咨詢組(consultative group of international agricultural research，CGIAR)提供的版本4.1的數(shù)據(jù)。

2)ICESat高度計數(shù)據(jù)。ICESat衛(wèi)星是NASA在2003年1月發(fā)射升空的，衛(wèi)星運行軌道高度為600 km，軌道傾角為94°。ICESat衛(wèi)星攜載的GLAS在地面上的腳印約為65 m，腳印間距約為172 m[10]。ICESat的高程數(shù)據(jù)是基于TOPEX/POSEIDON橢球的大地高程，相對于WGS84橢球約有70 cm的差距[11]。ICESat的工作周期為2003年2月—2009年10月。本文采用的是版本33的Level－2級別的全球陸表高度計(global land surface altimetry，GLA14)在 2004年2—3月(L2B)的產品。

3)MODIS地表類型產品(MCD12Q1)。MODIS地表覆蓋類型產品提供了全球的陸地地表類型數(shù)據(jù)，分辨率為500 m。本文采用的是由國際地圈－生物圈計劃(international geosphere－biosphere programme，IGBP)定義的地表覆蓋類型。地表類型數(shù)據(jù)的有效值為0～16，其中0和15分別代表水域和冰雪區(qū)域，254和255分別代表未分類和填充值。

2 研究方法

為了評估青藏高原地區(qū)SRTM高程數(shù)據(jù)的精度，本文首先提取和篩選了青藏高原地區(qū)的ICESat高度計數(shù)據(jù);然后采用雙線性插值算法提取了與ICESat高度計數(shù)據(jù)對應位置上的SRTM高程值，通過計算評估了青藏高原地區(qū)SRTM數(shù)據(jù)的精度;最后分析了SRTM數(shù)據(jù)與地形因子間的關系。

2.1 高度計數(shù)據(jù)處理

使用NASA提供的GLAS高度計高程提取工具(NASA GLAS altimetryelevation extractortool，NGAT)從ICESat高度計數(shù)據(jù)中提取必要的參量。首先，剔除明顯異常的數(shù)據(jù)點，這部分數(shù)據(jù)點是受云影響產生的異常值;其次，鑒于陸地表面的高自然反射率和低功率發(fā)射激光脈沖會產生過度飽和的問題，故對數(shù)據(jù)進行飽和度改正;再次，將反射率值大于1的數(shù)據(jù)剔除，以避免不可靠的高程值;最后，因為由云和水蒸氣等大氣成分引起的散射會使陸地表面的信號減小、返回波形變寬，這部分數(shù)據(jù)為高增益數(shù)據(jù)，故需要將增益值大于30的數(shù)據(jù)剔除[12]。

2.2 SRTM高程值提取

2.2.1 基準轉換

ICESat高度計數(shù)據(jù)和與ICESat腳印對應位置的SRTM數(shù)據(jù)基于不同參考橢球和高程基準，因此，需先對ICESat高度計數(shù)據(jù)進行橢球基準和高程基準轉換。這2個橢球在經線方向的差距為0，在緯線方向的差距也非常小，故只需要考慮其在高程方向的差距?；鶞兽D換的公式為

式中:Hellip為大地高程;Hortho為正射高程;Geoid為大地水準面差距;△h為2個橢球的高程之差。其中Geoid和△h的值可以直接從原始數(shù)據(jù)中提取。

2.2.2 高程值提取

根據(jù)MODIS數(shù)據(jù)(MCD12Q1)，剔除位于水域和冰雪區(qū)域的數(shù)據(jù)點(約9 214個)，按ICESat衛(wèi)星的每一個腳印點的位置P，通過雙線性插值算法獲取與該腳印點相對應位置上的SRTM高程值。雙線性插值的算法如下:

已知函數(shù) f在 Q11=(x1，y1)，Q12=(x1，y2)，Q21=(x2，y1)以及 Q22=(x2，y2)這 4 個點的值(圖2)，求未知函數(shù)f在點P=(x，y)的值。

圖2 雙線性插值算法Fig.2 Bilinear interpolation algorithm

圖2 中黑色點Q11，Q12，Q21和Q22為已知的4個數(shù)據(jù)點。通過X方向的線性插值，得到

通過Y方向的插值，得到

這樣就可以求得點P的函數(shù)值，也就是對應于ICESat腳印點位置的SRTM上的高程值。

盡管SRTM和ICESat數(shù)據(jù)已經轉化到了同一橢球基準和高程基準，本文仍然對兩者的基準一致性進行了檢查:沿經線和緯線方向分別平移ICESat高度計數(shù)據(jù)，然后計算兩者的相關系數(shù)。經計算比較，相關系數(shù)沒有明顯提高。因此，本文認為ICESat和SRTM數(shù)據(jù)的基準之間沒有偏差，并據(jù)此分析了ICESat高度計數(shù)據(jù)和SRTM數(shù)據(jù)的相關性。

2.3 高差計算及SRTM精度評估

首先計算ICESat和SRTM數(shù)據(jù)之間的高差hdiff，即

式中:hSRTM為SRTM上的高程值;hICESat為ICESat高度計數(shù)據(jù)的高程值。

然后對數(shù)據(jù)剔除高差超過100 m的數(shù)據(jù)點，再按照3倍標準差的準則剔除剩余的異常值(這些異常點主要是受到云霧影響的數(shù)據(jù)點)。類似的方法在文獻[5－7]中也用到過。前者使用100 m作為閾值，后者使用50 m作為閾值。

目前常用的DEM高程精度評價標準是DEM高程中誤差(即RMSE)模型。本文中的RMSE計算公式為

式中n為數(shù)據(jù)點總數(shù)。

2.4 地形因子與SRTM數(shù)據(jù)關系分析

坡度和坡向是地形的2個重要特征。本文使用ArcGIS軟件，用SRTM的DEM數(shù)據(jù)生成坡度和坡向數(shù)據(jù)。坡度表示單元格與其相鄰單元格之間的最大變化率;坡向則表示該單元格與其相鄰單元格的最大變化方向。坡向沿順時針方向從0°增加到360°，N方向定義為0°。將坡向分為8個類別:N[0°，22.5°)和[337.5°，360°]，NE[22.5°，67.5°)，E[67.5 °，112.5°)，SE[112.5°，157.5°)，S[157.5°，202.5°)，SW[202.5 °，247.5°)，W[247.5°，292.5°)，NW[292.5°，337.5°)。

把坡度數(shù)據(jù)分為10類，依此類別分析了高程精度、高差(hICESat－h(huán)SRTM)和坡度的關系，并分析了SRTM誤差值偏大或偏小的數(shù)據(jù)在坡向數(shù)據(jù)的8個類別中所占比例。

3 結果和討論

3.1 ICESat和SRTM 高程數(shù)據(jù)比較

研究區(qū)共有GLAS原始數(shù)據(jù)點319 869個。本文首先剔除了位于湖泊和冰雪區(qū)域的數(shù)據(jù)點，然后采用3倍標準差準則剔除掉異常數(shù)據(jù)，剩余GLAS數(shù)據(jù)點276 394個;并分別找出與之對應位置的SRTM像元，統(tǒng)計結果如表1所示。

表1 ICESat和SRTM高程數(shù)據(jù)基本統(tǒng)計量Tab.1 Basic statistics of ICESat and SRTM elevation data

經計算，2個數(shù)據(jù)集的相關系數(shù)為0.999 8。圖3示出青藏高原地區(qū)的ICESat和SRTM數(shù)據(jù)的散點圖及線性回歸圖。2種數(shù)據(jù)間存在很強的相關性，斜率和判定系數(shù)R2分別為1.000 2和0.999。

圖3 ICESat與SRTM高程數(shù)據(jù)關系Fig.3 Relationship between ICESat and SRTM elevation data

計算了2種數(shù)據(jù)之間的高差。圖4是青藏高原地區(qū)高差數(shù)據(jù)的空間分布圖。

圖4 高差的空間分布圖Fig.4 Spatial distribution of elevation differences

從圖4可以看出在青藏高原地區(qū)2種數(shù)據(jù)間的高差分布情況。圖5是青藏高原地區(qū)高差直方圖，高差數(shù)據(jù)呈左偏、尖峰分布。

圖5 高差直方圖統(tǒng)計分析Fig.5 Histogram analysis of elevation differences

經計算，整個區(qū)域2種數(shù)據(jù)的高差最大值和最小值分別為46.91 m和－51.54 m。SRTM高程數(shù)據(jù)的系統(tǒng)誤差為2.36±16.48 m，研究區(qū)SRTM 數(shù)據(jù)的高程精度約為16.65 m。

3.2 地形因子對高差的影響

表2示出整個實驗區(qū)SRTM數(shù)據(jù)的地形因子和高差數(shù)據(jù)的基本統(tǒng)計量。

表2 地形因子和高差的基本統(tǒng)計量Tab.2 Basic statistics of terrain factors and elevation differences

圖6分別為坡度與高差(圖6(a))和坡度與高差絕對值(圖6(b))的散點圖，此結果表明SRTM的精度隨著坡度的減小而提高。

圖6 坡度與高差的散點圖Fig.6 Scatter diagram of slope and elevation differences

根據(jù)研究區(qū)域的坡度特征，把坡度分為10類， 并計算了各類別坡度對應的RMSE值(圖7)。

圖7 坡度和RMSE的散點圖Fig.7 Scatter diagram of slope and RMSE

從圖7可以看出，當坡度低于區(qū)間[25°，30°]，RMSE隨坡度增加而顯著增加;當坡度超過該區(qū)間，RMSE隨坡度增加而降低，但RMSE依然明顯超過坡度低于20°的地區(qū);另外，當坡度超過了25°，RMSE變化不大。在地形平坦區(qū)域(坡度小于1°)，高差的平均值和標準差為－2.49±2.20 m，RMSE=3.33 m，遠遠優(yōu)于SRTM標稱的16 m的精度。然而，當坡度位于區(qū)間[25°，30°]，RMSE 會增加到約29 m，大約是地形平坦區(qū)域的9倍。

Jarvis等[13]發(fā)現(xiàn)，當坡向的朝向為 N，NE 和 E時，SRTM的高程值比當?shù)氐臄?shù)字化地形圖的高程值要高;當坡向的朝向為S，SW和W時，則相反;并認為這和雷達成像時的入射角有關。Zhao等[14]通過比較在華北平原地區(qū)的實測高程值和SRTM，發(fā)現(xiàn)SRTM的測量值在N方向的數(shù)值明顯高于實測值，而在其他方向上則不明顯。Shortridge等[4]通過比較美國地質調查局的國家高程數(shù)據(jù)(NED)和SRTM數(shù)據(jù)發(fā)現(xiàn)，與NED相比，SRTM在NW方向偏高，在SE方向偏低;并認為這與SRTM在升軌和降軌時的航向密切相關。為驗證SRTM數(shù)據(jù)在青藏高原是否與坡向相關，將坡向分類，并統(tǒng)計SRTM和ICESat數(shù)據(jù)之間差異較大的數(shù)據(jù)在各個坡向類別中的個數(shù)。圖8示出276 394個數(shù)據(jù)點對應位置的坡向。

圖8 數(shù)據(jù)點在8個坡向上的分布Fig.8 Distribution of data in eight direction categories

從圖8可以看出，除了NW和N方向的數(shù)據(jù)比其他方向的數(shù)據(jù)少外，另外6個坡向內的數(shù)據(jù)個數(shù)則基本一致。將 SRTM和 ICESat數(shù)據(jù)差異超過±20 m的數(shù)據(jù)點定義為大誤差數(shù)據(jù)，并統(tǒng)計大誤差數(shù)據(jù)在坡向類別中所占的百分比(圖9)。

圖9 SRTM大誤差數(shù)據(jù)點在各坡向類別中所占的百分比Fig.9 Percentage of SRTM big error data in each direction categories

從圖9可以看出，SRTM數(shù)據(jù)在坡向為N，NW和NE方向的測量值偏高，最大偏高在N向，達0.31左右;而在坡向為S，SE和SW方向的測量值偏低，S向的偏低最顯著，平均可達0.2。

4 結論

SRTM提供了一種免費的高精度的全球DEM數(shù)據(jù)，現(xiàn)已被廣泛應用在眾多科學研究中。本文用ICESat/GLAS激光高度計數(shù)據(jù)檢驗了中國青藏高原地區(qū)SRTM的高程精度，并討論了SRTM數(shù)據(jù)與地形因子之間的關系。得出結論如下:

1)中國青藏高原地區(qū)的SRTM和ICESat數(shù)據(jù)高度相關。盡管在本文實驗中這2種數(shù)據(jù)的獲取時間有4 a的間隔，但兩者的相關系數(shù)依然達到了0.999 8。

2)SRTM和ICESat這2種數(shù)據(jù)之間的系統(tǒng)差為2.36±16.48 m，其中SRTM的測量值偏高。整個實驗區(qū)的SRTM的精度達到RMSE=16.65 m，比數(shù)據(jù)標稱的16 m略大。

3)地形因子對SRTM的數(shù)據(jù)質量有一定影響。當坡度低于25°時，坡度和SRTM的高程精度呈現(xiàn)明顯的負相關;在地形平坦地區(qū)，SRTM的垂直精度約為3.33 m;然而在地形陡峭區(qū)域，SRTM的精度只有 29.00 m。

4)SRTM數(shù)據(jù)在青藏高原地區(qū)坡向為N，NW和NE方向的測量值偏高，而在坡向為S，SE和SW方向的測量值偏低。

盡管本文中ICESat高度計數(shù)據(jù)覆蓋了各種地貌類型，且精度較高，但與實測數(shù)據(jù)尚有一定誤差。在今后的研究中，可用更高精度的高度計數(shù)據(jù)(如ICESat－2等)評估并改善現(xiàn)有SRTM的數(shù)據(jù)質量。

[1]Hirt C，F(xiàn)ilmer M S，F(xiàn)eatherstone W E.Comparison and validation of the recent freely－available ASTER－GDEM ver1，SRTM ver 4.1 and GEODATA DEM－9S ver3 digital elevation models over Australia[J].Australian Journal of Earth Sciences，2010，57(3):337－347.

[2]Braun A，F(xiàn)otopoulos G.Assessment of SRTM，ICESat and survey control monument elevations in Canada[J].Photogrammetric Engineering and Remote Sensing，2007，73(12):1333－1342.

[3]Rodríguez E，Morris C S，Belz J E.A global assessment of the SRTM performance[J].Photogrammetric Engineering and Remote Sensing，2006，72(3):249－260.

[4]Shortridge A，Messina J.Spatial structure and landscape associations of SRTM error[J].Remote Sensing of Environment，2011，115(6):1576－1587.

[5]Beaulieu A，Clavet D.Accuracy assessment of Canadian digital elevation data using ICESat[J].Photogrammetric Engineering and Remote Sensing，2009，75(1):81－86.

[6]Carabajal C C，Harding D J.ICESat validation of SRTM C－ band digital elevation models[J].Geophysical Research Letters，2005，32(22):L22S01.1－ L22S01.5.

[7]Carabajal C C，Harding D J.SRTM C－ band and ICESat laser altimetry elevation comparisons as a function of tree cover and relief[J].Photogrammetric Engineering and Remote Sensing，2006，72(3):287－298.

[8]Zhang G Q，Xie H J，Kang S C，et al.Monitoring lake level changes on the Tibetan Plateau using ICESat altimetry data(2003－2009)[J].Remote Sensing of Environment，2011，115(7):1733－1742.

[9]Farr T G，Rosen P A，Caro E，et al.The shuttle radar topography mission[J].Reviews of Geophysics，2007，45(2):RG2004.1－RG2004.33.

[10]Schutz B E，Zwally H J，Shuman C A，et al.Overview of the ICESat Mission[J].Geophysical Research Letters，2005，32(21):L21S01.1－ L21S01.4.

[11]Bhang K J，Schwartz F W，Braun A.Verification of the vertical error in C－band SRTM DEM using ICESat and Landsat－7，Otter Tail County，MN[J].IEEE Transactions on Geoscience and Remote Sensing，2007，45(1):36－44.

[12]黃海蘭.利用ICESat和GRACE衛(wèi)星觀測數(shù)據(jù)確定極地冰蓋變化[D].武漢:武漢大學，2011.Huang H L.Determination of Polar Ice Sheet Change from ICESat and Grace Satellite Observations[D].Wuhan:Wuhan University，2011.

[13]Jarvis A，Rubiano J，Nelson A，et al.Practical use of SRTM data in the tropics:Comparisons with digital elevation models generated from cartographic data[J].International Center for Tropical Agriculture，2004，198:1－32.

[14]Zhao S M，Cheng W M，Zhou C H，et al.Accuracy assessment of the ASTER GDEM and SRTM3 DEM:An example in the Loess Plateau and North China Plain of China[J].International Journal of Remote Sensing，2011，32(23):8081－8093.