基于ICESat/GLAS高度計(jì)數(shù)據(jù)的SRTM數(shù)據(jù)精度評(píng)估——以青藏高原地區(qū)為例

萬(wàn) 杰，廖靜娟，許 濤，沈國(guó)狀

(1.中國(guó)科學(xué)院遙感與數(shù)字地球研究所數(shù)字地球科學(xué)重點(diǎn)實(shí)驗(yàn)室，北京 100094;2.中國(guó)科學(xué)院大學(xué)，北京 100049)

0 引言

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

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

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

1.1 研究區(qū)概況

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

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

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

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

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

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

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

2 研究方法

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

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

使用NASA提供的GLAS高度計(jì)高程提取工具(NASA GLAS altimetryelevation extractortool，NGAT)從ICESat高度計(jì)數(shù)據(jù)中提取必要的參量。首先，剔除明顯異常的數(shù)據(jù)點(diǎn)，這部分?jǐn)?shù)據(jù)點(diǎn)是受云影響產(chǎn)生的異常值;其次，鑒于陸地表面的高自然反射率和低功率發(fā)射激光脈沖會(huì)產(chǎn)生過(guò)度飽和的問(wèn)題，故對(duì)數(shù)據(jù)進(jìn)行飽和度改正;再次，將反射率值大于1的數(shù)據(jù)剔除，以避免不可靠的高程值;最后，因?yàn)橛稍坪退魵獾却髿獬煞忠鸬纳⑸鋾?huì)使陸地表面的信號(hào)減小、返回波形變寬，這部分?jǐn)?shù)據(jù)為高增益數(shù)據(jù)，故需要將增益值大于30的數(shù)據(jù)剔除[12]。

2.2 SRTM高程值提取

2.2.1 基準(zhǔn)轉(zhuǎn)換

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

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

2.2.2 高程值提取

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

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

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

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

通過(guò)Y方向的插值，得到

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

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

2.3 高差計(jì)算及SRTM精度評(píng)估

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

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

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

目前常用的DEM高程精度評(píng)價(jià)標(biāo)準(zhǔn)是DEM高程中誤差(即RMSE)模型。本文中的RMSE計(jì)算公式為

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

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

坡度和坡向是地形的2個(gè)重要特征。本文使用ArcGIS軟件，用SRTM的DEM數(shù)據(jù)生成坡度和坡向數(shù)據(jù)。坡度表示單元格與其相鄰單元格之間的最大變化率;坡向則表示該單元格與其相鄰單元格的最大變化方向。坡向沿順時(shí)針?lè)较驈?°增加到360°，N方向定義為0°。將坡向分為8個(gè)類別: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)和坡度的關(guān)系，并分析了SRTM誤差值偏大或偏小的數(shù)據(jù)在坡向數(shù)據(jù)的8個(gè)類別中所占比例。

3 結(jié)果和討論

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

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

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

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

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

計(jì)算了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)計(jì)分析Fig.5 Histogram analysis of elevation differences

經(jīng)計(jì)算，整個(gè)區(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 地形因子對(duì)高差的影響

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

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

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

圖6 坡度與高差的散點(diǎn)圖Fig.6 Scatter diagram of slope and elevation differences

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

圖7 坡度和RMSE的散點(diǎn)圖Fig.7 Scatter diagram of slope and RMSE

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

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

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

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

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

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

4 結(jié)論

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

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

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

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

4)SRTM數(shù)據(jù)在青藏高原地區(qū)坡向?yàn)镹，NW和NE方向的測(cè)量值偏高，而在坡向?yàn)镾，SE和SW方向的測(cè)量值偏低。

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

[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)星觀測(cè)數(shù)據(jù)確定極地冰蓋變化[D].武漢:武漢大學(xué)，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.