基于ICESat/GLAS高度计数据的SRTM数据精度评估——以青藏高原地区为例

万 杰，廖静娟，许 涛，沈国状

(1.中国科学院遥感与数字地球研究所数字地球科学重点实验室，北京 100094;2.中国科学院大学，北京 100049)

0 引言

数字高程模型(DEM)在地球科学及其他相关学科(如重力场建模、水文学、地形图制图学和航空图像正射纠正等)研究中有广泛的应用[1]。航天飞机雷达测图计划(shuttle Radar topography mission，SRTM)是第一个公开发布全球高精度数字高程模型的遥感平台。国内外研究采用不同的参照数据(如 GPS[2－3])，参考 DEM[4]和高度计数据[5－6]评估SRTM的精度。由于SRTM的版本、高度计数据的版本、实验样区和参照数据的变化，上述研究及其结论各有差异。Carabajal等[7]通过对美国西部、南美洲(亚马孙流域)、非洲、亚洲和澳大利亚部分区域的SRTM的分析，发现不同区域SRTM的高程误差呈现出一定规律的变化。Rodriguez等[3]通过比较SRTM和车载KGPS导航仪数据，得到SRTM的水平和垂直精度分别为11.9 m和9 m。此外，植被覆盖度和地形起伏的增大也会明显降低SRTM的数据质量[6]。

美国国家航空航天管理局(NASA)的冰、云和陆地高程卫星(ICESat)携载的激光高度计(geoscience laser altimeter system，GLAS)数据在垂直方向上的测量精度达到了cm级，远远超过了SRTM数据的精度[2]，因此可以利用ICESat/GLAS激光高度计数据(以下简称ICESat高度计数据)评估SRTM高程数据的精度。此外，相对于通常处于视野开阔和平坦区域的GPS数据和控制点不能遍及整个区域的局限性，ICESat脚印点(footprint point)可分布在整个实验样区。因此，本文采用精度更高的ICESat高度计数据分析了SRTM的系统误差和高程精度，并探讨了SRTM高程数据的误差特征与地形因子(坡度和坡向)间的关系。

1 研究区概况与数据源

1.1 研究区概况

本文的研究区域位于中国的青藏高原。青藏高原东西方向长约2 945 km，南北方向长约1 532 km，覆盖了西藏自治区全境和云南、四川、甘肃、青海和新疆的部分区域，总面积约为2 570 000 km2。除两极区域之外，青藏高原是地球上最寒冷的区域，也被称为“世界第三极”。青藏高原是一个地形条件极为复杂的区域，有着众多的山脉及占地面积大约49 873 km2的36 800 条冰川和 1 500 多个湖泊[8]。青藏高原的高程变化在60～8 792 m之间，平均海拔在4 000 m以上(图1)。

图1 研究区DEMFig.1 DEM of study area

1.2 数据源

本文实验使用的数据包括SRTM高程数据、ICESat高度计数据和 MODIS地表类型产品(MCD12Q1)。

1)SRTM高程数据。SRTM于2000年2月发射，获取了覆盖地球N 60°～S 57°之间地表的数字高程。SRTM免费提供30 m分辨率(美国境内)和90 m分辨率(其他地区)的高程数据产品。SRTM DEM数据的高程基准是EGM96的大地水准面，平面基准是WGS84[6]。SRTM数据设计的水平和垂直精度分别高于20 m和16 m[9]。本文使用的SRTM高程数据是由国际农业研究咨询组(consultative group of international agricultural research，CGIAR)提供的版本4.1的数据。

2)ICESat高度计数据。ICESat卫星是NASA在2003年1月发射升空的，卫星运行轨道高度为600 km，轨道倾角为94°。ICESat卫星携载的GLAS在地面上的脚印约为65 m，脚印间距约为172 m[10]。ICESat的高程数据是基于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地表覆盖类型产品提供了全球的陆地地表类型数据，分辨率为500 m。本文采用的是由国际地圈－生物圈计划(international geosphere－biosphere programme，IGBP)定义的地表覆盖类型。地表类型数据的有效值为0～16，其中0和15分别代表水域和冰雪区域，254和255分别代表未分类和填充值。

2 研究方法

为了评估青藏高原地区SRTM高程数据的精度，本文首先提取和筛选了青藏高原地区的ICESat高度计数据;然后采用双线性插值算法提取了与ICESat高度计数据对应位置上的SRTM高程值，通过计算评估了青藏高原地区SRTM数据的精度;最后分析了SRTM数据与地形因子间的关系。

2.1 高度计数据处理

使用NASA提供的GLAS高度计高程提取工具(NASA GLAS altimetryelevation extractortool，NGAT)从ICESat高度计数据中提取必要的参量。首先，剔除明显异常的数据点，这部分数据点是受云影响产生的异常值;其次，鉴于陆地表面的高自然反射率和低功率发射激光脉冲会产生过度饱和的问题，故对数据进行饱和度改正;再次，将反射率值大于1的数据剔除，以避免不可靠的高程值;最后，因为由云和水蒸气等大气成分引起的散射会使陆地表面的信号减小、返回波形变宽，这部分数据为高增益数据，故需要将增益值大于30的数据剔除[12]。

2.2 SRTM高程值提取

2.2.1 基准转换

ICESat高度计数据和与ICESat脚印对应位置的SRTM数据基于不同参考椭球和高程基准，因此，需先对ICESat高度计数据进行椭球基准和高程基准转换。这2个椭球在经线方向的差距为0，在纬线方向的差距也非常小，故只需要考虑其在高程方向的差距。基准转换的公式为

式中:Hellip为大地高程;Hortho为正射高程;Geoid为大地水准面差距;△h为2个椭球的高程之差。其中Geoid和△h的值可以直接从原始数据中提取。

2.2.2 高程值提取

根据MODIS数据(MCD12Q1)，剔除位于水域和冰雪区域的数据点(约9 214个)，按ICESat卫星的每一个脚印点的位置P，通过双线性插值算法获取与该脚印点相对应位置上的SRTM高程值。双线性插值的算法如下:

已知函数 f在 Q11=(x1，y1)，Q12=(x1，y2)，Q21=(x2，y1)以及 Q22=(x2，y2)这 4 个点的值(图2)，求未知函数f在点P=(x，y)的值。

图2 双线性插值算法Fig.2 Bilinear interpolation algorithm

图2 中黑色点Q11，Q12，Q21和Q22为已知的4个数据点。通过X方向的线性插值，得到

通过Y方向的插值，得到

这样就可以求得点P的函数值，也就是对应于ICESat脚印点位置的SRTM上的高程值。

尽管SRTM和ICESat数据已经转化到了同一椭球基准和高程基准，本文仍然对两者的基准一致性进行了检查:沿经线和纬线方向分别平移ICESat高度计数据，然后计算两者的相关系数。经计算比较，相关系数没有明显提高。因此，本文认为ICESat和SRTM数据的基准之间没有偏差，并据此分析了ICESat高度计数据和SRTM数据的相关性。

2.3 高差计算及SRTM精度评估

首先计算ICESat和SRTM数据之间的高差hdiff，即

式中:hSRTM为SRTM上的高程值;hICESat为ICESat高度计数据的高程值。

然后对数据剔除高差超过100 m的数据点，再按照3倍标准差的准则剔除剩余的异常值(这些异常点主要是受到云雾影响的数据点)。类似的方法在文献[5－7]中也用到过。前者使用100 m作为阈值，后者使用50 m作为阈值。

目前常用的DEM高程精度评价标准是DEM高程中误差(即RMSE)模型。本文中的RMSE计算公式为

式中n为数据点总数。

2.4 地形因子与SRTM数据关系分析

坡度和坡向是地形的2个重要特征。本文使用ArcGIS软件，用SRTM的DEM数据生成坡度和坡向数据。坡度表示单元格与其相邻单元格之间的最大变化率;坡向则表示该单元格与其相邻单元格的最大变化方向。坡向沿顺时针方向从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°)。

把坡度数据分为10类，依此类别分析了高程精度、高差(hICESat－hSRTM)和坡度的关系，并分析了SRTM误差值偏大或偏小的数据在坡向数据的8个类别中所占比例。

3 结果和讨论

3.1 ICESat和SRTM 高程数据比较

研究区共有GLAS原始数据点319 869个。本文首先剔除了位于湖泊和冰雪区域的数据点，然后采用3倍标准差准则剔除掉异常数据，剩余GLAS数据点276 394个;并分别找出与之对应位置的SRTM像元，统计结果如表1所示。

表1 ICESat和SRTM高程数据基本统计量Tab.1 Basic statistics of ICESat and SRTM elevation data

经计算，2个数据集的相关系数为0.999 8。图3示出青藏高原地区的ICESat和SRTM数据的散点图及线性回归图。2种数据间存在很强的相关性，斜率和判定系数R2分别为1.000 2和0.999。

图3 ICESat与SRTM高程数据关系Fig.3 Relationship between ICESat and SRTM elevation data

计算了2种数据之间的高差。图4是青藏高原地区高差数据的空间分布图。

图4 高差的空间分布图Fig.4 Spatial distribution of elevation differences

从图4可以看出在青藏高原地区2种数据间的高差分布情况。图5是青藏高原地区高差直方图，高差数据呈左偏、尖峰分布。

图5 高差直方图统计分析Fig.5 Histogram analysis of elevation differences

经计算，整个区域2种数据的高差最大值和最小值分别为46.91 m和－51.54 m。SRTM高程数据的系统误差为2.36±16.48 m，研究区SRTM 数据的高程精度约为16.65 m。

3.2 地形因子对高差的影响

表2示出整个实验区SRTM数据的地形因子和高差数据的基本统计量。

表2 地形因子和高差的基本统计量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

根据研究区域的坡度特征，把坡度分为10类， 并计算了各类别坡度对应的RMSE值(图7)。

图7 坡度和RMSE的散点图Fig.7 Scatter diagram of slope and RMSE

从图7可以看出，当坡度低于区间[25°，30°]，RMSE随坡度增加而显著增加;当坡度超过该区间，RMSE随坡度增加而降低，但RMSE依然明显超过坡度低于20°的地区;另外，当坡度超过了25°，RMSE变化不大。在地形平坦区域(坡度小于1°)，高差的平均值和标准差为－2.49±2.20 m，RMSE=3.33 m，远远优于SRTM标称的16 m的精度。然而，当坡度位于区间[25°，30°]，RMSE 会增加到约29 m，大约是地形平坦区域的9倍。

Jarvis等[13]发现，当坡向的朝向为 N，NE 和 E时，SRTM的高程值比当地的数字化地形图的高程值要高;当坡向的朝向为S，SW和W时，则相反;并认为这和雷达成像时的入射角有关。Zhao等[14]通过比较在华北平原地区的实测高程值和SRTM，发现SRTM的测量值在N方向的数值明显高于实测值，而在其他方向上则不明显。Shortridge等[4]通过比较美国地质调查局的国家高程数据(NED)和SRTM数据发现，与NED相比，SRTM在NW方向偏高，在SE方向偏低;并认为这与SRTM在升轨和降轨时的航向密切相关。为验证SRTM数据在青藏高原是否与坡向相关，将坡向分类，并统计SRTM和ICESat数据之间差异较大的数据在各个坡向类别中的个数。图8示出276 394个数据点对应位置的坡向。

图8 数据点在8个坡向上的分布Fig.8 Distribution of data in eight direction categories

从图8可以看出，除了NW和N方向的数据比其他方向的数据少外，另外6个坡向内的数据个数则基本一致。将 SRTM和 ICESat数据差异超过±20 m的数据点定义为大误差数据，并统计大误差数据在坡向类别中所占的百分比(图9)。

图9 SRTM大误差数据点在各坡向类别中所占的百分比Fig.9 Percentage of SRTM big error data in each direction categories

从图9可以看出，SRTM数据在坡向为N，NW和NE方向的测量值偏高，最大偏高在N向，达0.31左右;而在坡向为S，SE和SW方向的测量值偏低，S向的偏低最显著，平均可达0.2。

4 结论

SRTM提供了一种免费的高精度的全球DEM数据，现已被广泛应用在众多科学研究中。本文用ICESat/GLAS激光高度计数据检验了中国青藏高原地区SRTM的高程精度，并讨论了SRTM数据与地形因子之间的关系。得出结论如下:

1)中国青藏高原地区的SRTM和ICESat数据高度相关。尽管在本文实验中这2种数据的获取时间有4 a的间隔，但两者的相关系数依然达到了0.999 8。

2)SRTM和ICESat这2种数据之间的系统差为2.36±16.48 m，其中SRTM的测量值偏高。整个实验区的SRTM的精度达到RMSE=16.65 m，比数据标称的16 m略大。

3)地形因子对SRTM的数据质量有一定影响。当坡度低于25°时，坡度和SRTM的高程精度呈现明显的负相关;在地形平坦地区，SRTM的垂直精度约为3.33 m;然而在地形陡峭区域，SRTM的精度只有 29.00 m。

4)SRTM数据在青藏高原地区坡向为N，NW和NE方向的测量值偏高，而在坡向为S，SE和SW方向的测量值偏低。

尽管本文中ICESat高度计数据覆盖了各种地貌类型，且精度较高，但与实测数据尚有一定误差。在今后的研究中，可用更高精度的高度计数据(如ICESat－2等)评估并改善现有SRTM的数据质量。

[1]Hirt C，Filmer M S，Featherstone 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，Fotopoulos 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卫星观测数据确定极地冰盖变化[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.