2. 中国科学院新疆生态与地理研究所, 新疆 乌鲁木齐 830011
2. Xinjiang Institute of Ecology and Geography, Chinese Academy of Sciences, Urumqi 830011, China
1 引 言
由于容易受到积雪及山体阴影的影响，在利用遥感技术提取冰川信息的过程中，遥感数据的时相选择、观测角度等会对冰川提取的准确度和精度造成较大影响。2002年全球陆地冰空间监测计划(global land ice measurements from space,GLIMS)发布了由遥感技术提取的全球陆地冰川编目数据，然而用于提取新疆山区冰川的遥感数据在时间跨度上长达30多年，冰川边界的精度和时效性存在较大问题[1]。目前，遥感冰川信息提取方法主要有分类法和阈值分割法[2]。分类方法适用于较小局部区域精细分类，效果较好，但算法复杂，难以推广；分割方法主要采用基于雪盖指数或比值图像[3]的单一阈值提取，然而针对大区域的冰川制图，简单的阈值分割无法适用于高山冰川区域复杂多样的局部环境。“全域-局部”迭代的阈值分割方法可实现大范围，复杂环境下地物提取阈值的快速选择[4,5]。
新疆地区冰川分布在地形陡峭的高山地区，受太阳入射角度的影响，部分冰川往往被山体阴影遮挡。文献[6]在提取喜马拉雅山地区冰川时采用具有较高太阳高度角的影像来减少山体阴影的影响。另外，由于积雪与冰川很难区分，而高山区域终年降雪，很难获取无积雪覆盖的遥感影像[7]。因此，采用单一时相的遥感影像很难获取精确的冰川边界。采用一定时段内的多时相遥感影像对冰川边界进行观测，综合利用多角度信息以减少地形对冰川的影响，可获取某个时段内冰川的最佳边界信息。为此，本文拟采用一定时段内的Landsat多时相遥感影像提取天山托木尔峰西部的冰川，通过不同太阳入射角度下的冰雪信息叠加分析，消除山体阴影和积雪对冰川信息的干扰，以获得较精确的冰川边界信息。 2 高山冰川提取的特点
新疆天山地区的冰川属于温带大陆型冰川，发育于海拔4000m以上的地形陡峭山区，其北坡由于受崎岖地形的遮挡，冰川发育良好，不少冰川位于巨大山体形成的阴影区域内[8]。在不同的太阳高度条件下，山体形成的阴影分布是不同的，可以通过多角度遥感的方法减少遮挡冰川的阴影面积。
通过对天山冰川变化最为剧烈地区的冰川观测记录发现，1972-2005年间哈密地区的庙尔沟冰川平均退缩量为5m/a，庙尔沟冰帽年退缩量为2.3m/a，乌鲁木齐东部的博格达冰川长度平均退缩量为4.11m/a，四工河4号冰川末端2006-2009年最大退缩量为13.3m/a[9]。如果采用空间分辨率为30m的Landsat TM数据(30m)监测冰川的变化，在冰川变化最剧烈的区域，可观测到冰川变化的时间尺度约为2.26a。因此，可认为Landsat提取的冰川边界在2年之内是相对稳定的。
为此，本文考虑采用某个时段内的多期遥感数据联合确定冰川的边界，其主要思想是合理利用不同时相遥感影像成像时的不同太阳入射角度来减小山体阴影对冰川遮挡的影响，而多期遥感数据联合观测可以获取冰川上的最小积雪覆盖面积。该方法基于以下两个假设条件：
(1) 利用多时相遥感数据的多角度信息可以最大程度减小甚至消除阴影对冰川提取的影响。
(2) 采用一定时段内(2年内)的多时相Landsat影像提取冰川边界时，可认为这个时间段内的冰川边界不变。
这样就可以用多时相遥感提取的冰川边界进行叠置分析，从而逼近冰川的真实边界。本文的多角度遥感冰川提取的技术路线如图 1所示，包括山体阴影提取、冰雪信息提取、多角度阴影剔除和多时相冰川边界叠置分析。
 |
| 图 1 多角度冰川提取流程 Fig. 1 Multi-angle glacier delineation flow |
本文基于高程数据和遥感成像的太阳入射角度，利用朗伯特余弦定律[10]建立山体阴影提取的物理模型，通过计算数字高程模型(digital elevation model,DEM)的中的坡面法线，结合太阳角度得到朗伯特余弦定律的入射角，最终计算出山体阴影图像。本文在地形建模工具中输入DEM数据和对应时相的太阳角度参数，生成地形晕渲图像，通过分析得到山体阴影。
图 2以天山西段托木尔峰东北部的几条冰川[11]为例，通过DEM(图 2(b))和Landsat TM影像(图 2(a)波段组合为7/4/2)对应的太阳角度求得该区域的地形晕渲图，提取出山体阴影(图 2(d))，所用Landsat TM影像成像时太阳方位角147.22°，高度角47.56°，成像日期为2011-09-13。图 2(c)是单纯利用光谱信息提取阴影区域叠加在Landsat TM影像上显示。由图 2可以看出，由高程数据结合卫星成像时的太阳提取的山体阴影效果较光谱信息提取的阴影好。与单纯利用遥感影像的光谱信息提取阴影相比，利用太阳角度结合DEM提取山体阴影能够减少因“同物异谱”或“异物同谱”造成的错分现象。
 |
| 图 2 山体阴影提取 Fig. 2 Extraction of mountain shadow |
本文采用“全域-局部”分步迭代的阈值提取方法[12]提取遥感图像上的冰雪信息。这种方法先用单一阈值对整景遥感图像进行全域阈值分割，获取初步的目标与背景信息，然后针对每个目标单元，采用分步迭代的方法对缓冲区内的目标进行精细分割。以图 3为例，“全域-局部”迭代法提取冰雪的技术流程如下：
 |
| 图 3 “全域－局部”冰雪提取以及其NDSI统计直方图 Fig. 3 Snow/ice extraction based on “global-local” segmentation and their NDSI histogram |
(1) 归一化雪覆盖指数NDSI的提取[13]。冰雪在短波红外波段有强吸收特性，在绿波段有强的反射特性，能很好地区分冰雪与背景信息[14]。以Landsat TM影像为例，归一化雪覆盖指数NDSI的计算公式如下
式中，ρB2、ρB5分别代表Landsat TM影像第2、第5波段的反射率。图 3第2行为4个冰川样例图中的归一化雪盖指数。
(2) 全域阈值分割。采用初始阈值T0对NDSI图像进行单阈值分割，获取初步的冰雪单元[16]。已有研究表明，当NDSI>0.4时，能够有效提取绝大部分的冰川和积雪[15]，这里也采用T0= 0.4进行阈值分割。由图 3可知，不同形态与状态的冰川，采用单一阈值提取冰川均取得了一定的效果，明显的冰川边界均能精确地提取，部分区域的背景信息也有较大的NDSI值，被划分为冰川单元。
(3) 局域阈值分割：在全域阈值分割的每个冰雪单元，计算与其面积相等的缓冲区，作为局域分割单元(图 3第2行)，并对局域单元的NDSI进行直方图统计，可以看到冰雪与背景的双峰分布图，通过双峰分布直方图阈值选择方法[12]确定新的分割阈值(图 3第3行)。图 3中，红色边界为全域分割冰雪边界；绿色边界为缓冲区边界，NDSI直方图中红色竖线为局部分割阈值；黄色边界为局部分割的冰雪边界；图 3(d)中紫色边界为水体边界。4个冰川的最终的局部分割阈值分别为：0.376、0.324、0.553、0.475，这就使得不同类型的冰雪单元依据自身不同的特点实现精细的提取。
因此，针对不同形态和特征的冰川，采用不同的阈值方法提取可得到较好的效果。然而，图 3(a)、图 3(b)、图 3(d) 的冰雪均有阴影的遮挡，且无法获取阴影区域的冰雪边界，部分地区的山体阴影也被视为冰川[17] 3.3 阴影区域冰雪边界提取
为了识别阴影区域内冰雪边界，这里采用一定时期内不同时相的遥感数据进行联合分析。不同时相数据对同一地区进行观测时具有不同的太阳角度，所形成的山体阴影也不同。通过对多期阴影信息的空间分析获取阴影分布最小的冰雪范围的方法如下：
(1) 提取出多期影像受阴影影响的总和，确定多角度遥感影像总的阴影区域Smax
设有N个时相的遥感影像用于多角度剔除山体阴影，S(i)表示时相i影像的山体阴影区域，i∈[1,N]。
(2) 为了使总阴影区域提取的冰川不受积雪干扰，确定总阴影区内的冰川边界。先求出时相i对应的总的阴影区域的冰雪为GS(i)
Gt(i)为3.2节方法提取的时相i冰雪。
然后，剔除总阴影区各个时相冰川受积雪的影响，得到总阴影区的冰川GS
(3) 确定i时相剔除阴影后冰雪G(i)
G(i)表示时相i影像确定阴影区域冰川后的冰雪边界，即G(i)剔除了山体阴影的影响。判断多角度方法能完全剔除山体阴影的条件是
，若无法完全剔除阴影，标记未剔除的阴影。图 4选取了4幅不同时相的Landsat数据进行阴影剔除试验，成像日期分别为2011-02-17、2011-10-22、2011-09-29、2011-06-25，太阳高度角分别为31.20°、34.56°、42.32°、63.53°，太阳方位角分别为150.42°、157.50°、152.28°、126.81°。图 4(a)、(b)、(c)、(d)的底图分别为4个时相Landsat数据的波段7、4、2合成影像，图中无填充和有填充的多边形区域分别为利用其余时相依次剔除图 4(a)积雪和山体阴影后的结果。由图 4可见，试验区在春、秋、冬季的影像受山体阴影面积较大，而夏季影像的山体阴影面积最小。阴影区的面积大小跟太阳高度角成反比。通过对4个阴影区的空间分析，可以获取山体阴影的最小范围(图 4 (d))。虽然图 4 (d)中仍有未剔除的4段阴影，但阴影对冰雪提取的影响降低到最低程度。在有更理想的夏季时相数据的情况下，可以完全消除阴影的影响。
 |
| 图 4 山体阴影剔除与冰川提取示意图 Fig. 4 The diagram to eliminate mountain shadow and to extract glacier |
上节利用多角度信息确定了山体阴影区域的冰川边界，从而得到每个时相剔除山体阴影影响的冰雪信息。由于没有气象数据的辅助，无法确定单一时相遥感影像中哪些冰川边界不受积雪影响。本文求出多时相冰雪最小边界为理想的冰川边界，在3.3节得到时相i剔除阴影影响后的冰雪为G(i)，则多时相冰川边界Gmulti的计算公式如下
G(i)表示时相i影像确定阴影区域冰川边界后的冰雪边界。图 4(a)、(b)、(c)、(d)中相近时间段内各时相的冰雪分布不同，通过3.3节中的方法分别剔除四个时相的山体阴影，得到对应时相完整的冰雪边界G(1)、G(2)、G(3)、G(4)，利用公式(5)求得冰川边界Gmulti(见图 4(d))，最大限度地排除了积雪的影响。 4 试验与讨论
为验证本文提出方法的有效性，本文使用天山托木尔峰西部地区Landsat TM数据和SRTM 高程数据进行多角度冰川提取试验，试验区的范围为79.33°E至80.16°E，42.12°N至42.55°N。遥感数据采用了2009年至2010年4个时相Landsat TM数据(见表 1)，轨道号为147/031。
| 时相n | 日期 | 太阳 高度角/(°) | 太阳 方位角/(°) | 山体阴影下 冰雪面积/km2 | 积雪面积 /km2 | 山体阴影下的冰雪 比例A′n/(%) | 积雪比例 SC′n/(%) |
| 1 | 2010-10-22 | 34.56 | 157.49 | 229.9 | 945.4 | 16.5 | 67.7 |
| 2 | 2010-10-03 | 41.02 | 153.97 | 174.6 | 480 | 18.8 | 51.5 |
| 3 | 2010-07-15 | 61.77 | 128.54 | 25.8 | 235 | 3.8 | 34.2 |
| 4 | 2009-09-14 | 47.23 | 148.16 | 2.9 | 1919.8 | 0.1 | 81.0 |
利用第3节方法提取试验区冰川的步骤为：①利用30m分辨率的DEM与表 2中各时相Landsat TM数据观测的太阳角度提取该区域的各个时相的山体阴影，分析得4个时相山体阴影的交集为空集；②通过NDSI “全域-局部”阈值迭代的方法提取单一时相对应的冰雪边界；③多时相的冰雪边界与阴影边界叠加运算，计算出每个时相山体阴影区域内的最小的冰雪边界，确定山体阴影区域的冰川边界，得到多角度剔除山体阴影后的冰雪结果(图 5(a)、(b)、(c)、(d))；④用4个时相的剔除山体阴影后的冰雪边界求最小值得到冰川边界，得到的冰川边界叠加在该区域2009年的RapidEye高分辨率(分辨率为5m)影像上(图 5(e))。图 5(a)、(b)、(c)、(d)的底图分别为4个时相遥感数据的波段7、4、2合成影像，多边形区域为对应时相剔除阴影后的冰雪边界。图 5(e)的底图为2009年高分辨率影像，多边形区域为利用图 5(a)、(b)、(c)、(d)的冰雪边界提取的冰川边界。
 |
| 图 5 多角度提取冰川边界 Fig. 5 Delineation of glacier boundary using multi-angle glaciers |
山体阴影区域的冰雪的比例表示对应时相的冰雪受山体阴影影响的程度[18,19]，计算方法为某一时相山体阴影区域冰雪面积除以该时相的总冰雪面积。
时相n的冰雪面积Gn减去冰川面积Gmulti，得到该时相的积雪面积，再除以Gn，得到该时相的积雪比例。本文方法提取2009至2010年该区域的冰川面积为Gmulti=568.7km2。由表 2分析可得试验区2009年9月14日山体阴影影响最小，2010年7月15日的积雪范围最小(见表 1)。
结合纹理、形态、邻近度与DEM数据的高度和坡度信息，采用面向对象方法可以快速提取高分辨率影像中的冰川[20]。但是对比高分辨率方法和本文方法(见图 6)。
 |
| 图 6 结果对比 Fig. 6 Comparison results |
本文方法对于冰川提取显示出很好的优势，具体分析如下：
(1) 本文方法采用的Landsat数据时相较为齐全，而高分辨率影像的时相单一。尤其是在积雪影响较为严重的时期，单一时相的高分辨率影像提取的冰川可以通过形状、光谱、地形等特征将冰川、积雪分离，虽然其描述冰川的轮廓极为精细，但是一些与冰川特征相似的积雪却无法分离。而连续的多时相的Landsat影像，可以综合不同时相的冰川边界，最大限度地区分冰川和积雪的边界。图 6(a)显示了多角度提取的冰川边界与受积雪影响较小的高分影像的冰川边界线非常一致，而且在描述冰川内部的裸地时，多角度方法要比面向对象的高分辨率方法准确。高分辨率影像在描述有很多细小分支的冰川时极为精细。但是，对于较长时间(30年或者更长时间)的冰川监测而言，本文方法所用数据的时间的长度和密度都显示出很大优势。
(2) Landsat影像与高分辨率影像的成像条件有一定的差异，尤其是在几何形变较为严重高山地区表现尤为明显，因此针对地势起伏剧烈的高山地区冰川的提取，选择近似正射投影(几何形变较小)的本文冰川提取方法较好(图 6(b))。
(3) 图 6(c)通过面向对象方法提取高分辨率影像的阴影区域的冰川，非冰川信息与冰川信息在光谱和地形特征等方面极为相似，不容易区分。图 6(d)显示本文方法在相同山体阴影区域结合不同太阳角度下的冰川分布信息，能够很好地提取出阴影区域的冰川边界，并且能够保证提取的阴影区域冰川边界受积雪影响最小。 5 结 论
针对冰川提取容易受到积雪和山体阴影的影响，本文在“全域－局部”阈值分割提取冰雪的基础上，用多角度多时相遥感数据对短时期内的冰川进行提取，并对该方法提取冰川的精度进行了验证。结论如下：
(1) 本文用DEM和遥感成像的太阳角度信息建立山体阴影的物理模型。这种方法提取的山体阴影准确度高，与单纯利用遥感影像的光谱信息提取阴影相比，减少了因“异物同谱”造成的错分现象。
(2) 冰川退缩是天山地区冰川变化的总体趋势，冰川提取中把冰雪变化最剧烈的情况作为理想的冰川边界[6]。从本文多角度山体阴影剔除方法可以看出，山体阴影的分布与太阳的高角度有关，在选择多时相数据进行冰川提取时，可根据太阳高度信息来选择合适的遥感数据进行冰川提取。
(3) 多期遥感影像可以剔除积雪的影响。由结果分析可得，多角度遥感提取有效地融合了积雪范围最小和山体阴影影响最小两个最佳遥感观测条件，提取的冰川最接近真实情况。
| [1] | LI X, CHENG G D, JIN H J, et al. Cryospheric Change in China[J]. Global and Planetary Change, 2008, 62(3): 210-218. |
| [2] | PAN B T, ZHANG G L, WANG J, et al. Glacier Changes from 1966—2009 in the Gongga Mountains, on the South-Eastern Margin of the Qinghai-Tibetan Plateau and Their Climatic Forcing[J]. The Cryosphere, 2012, 6: 1087-1101. |
| [3] | LI Zhen, SHI Jiancheng. Snow Mapping Algorithm Development and Validation Using Hyperspectral Data[J]. Acta Geodaetica et Cartographica Sinica, 2001, 30(1): 67-72.(李震, 施建成. 高光谱遥感积雪制图算法及验证[J]. 测绘学报, 2001, 30(1): 67-72.) |
| [4] | ZHU Changming, LUO Jiancheng, SHEN Zhanfeng, et al. River Liner Water Adaptive Auto-extraction on Remote Sensing Image Aided by DEM[J]. Acta Geodaetica et Cartographica Sinica, 2013, 42(2): 277-283. (朱长明, 骆剑承, 沈占锋, 等. DEM辅助下的河道细小线性水体自适应迭代提取[J]. 测绘学报, 2013, 42(2): 277-283.) |
| [5] | HU Xiaodong, LUO Jiancheng, XIA Liegang, et al. Adaptive Water Body Information Extraction Using RS TUPU Computing Model [J]. Acta Geodaetica et Cartographica Sinica, 2011, 40(5): 544-550. (胡晓东, 骆剑承, 夏列钢. 图谱迭代反馈的自适应水体信息提取方法[J]. 测绘学报, 2011, 40(5): 544-550.) |
| [6] | FREY H, PAUL F, STROZZI T. Compilation of a Glacier Inventory for the Western Himalayas from Satellite Data: Methods, Challenges, and Results[J]. Remote Sensing of Environment, 2012, 124(9): 832-843. |
| [7] | PAUL F, KAAB A, MAISCH M; et al. The New Remote Sensing Derived Swiss Glacier Inventory: I: Methods[J]. Annals of Glaciology, 2002, 34(1): 355-361. |
| [8] | BOLCH T. Climate Change and Glacier Retreat in Northern Tienshan (Kazakhstan/Kyrgyzstan) Using Remote Sensing Data[J]. Global and Planetary Change, 2007, 56(1): 1-12. |
| [9] | WANG P Y, Li Z Q, Li H L, et al. Comparison of Glaciological and Geodetic Mass Balance at Urumqi Glacier No.1, Tian shan, Central Asia[J]. Global and Planetary Change, 2014, 114: 14-22. |
| [10] | KATZIL Y, DOYTSHER Y. A Logarithmic and Sub-pixel Approach to Shaded Relief Representation[J]. Computers & Geosciences, 2003, 29(9): 1137-1142. |
| [11] | RAUP B, RACOVITEANU A, KHALSA S J S, et al. The GLIMS Geospatial Glacier Database: A New Tool for Studying Glacier Change[J]. Global and Planetary Change, 2007, 56(1): 101-110. |
| [12] | LI J L, SHENG Y W. An Automated Scheme for Glacial Lake Dynamic Mapping with Landsat Imagery and Digital Elevation Models: A Case Study in the Himalayas[J]. International Journal of Remote Sensing, 2012, 33(16): 5194-5213. |
| [13] | XIAO J Y, JI N, LI X. A New Index for Steel Framed Roof Information Extraction Based on Remote Sensing[J]. Advanced Materials Research, 2013, 726: 4682-4685. |
| [14] | YIN D, CAO X, CHEN X, et al. Comparison of Automatic Thresholding Methods for Snow-cover Mapping Using Landsat TM Imagery[J]. International Journal of Remote Sensing, 2013, 34(19): 6529-6538. |
| [15] | BISHOP M P, BONK R, KAMP J U, et al. Terrain Analysis and Data Modeling for Alpine Glacier Mapping[J]. Polar Geography, 2001, 25(3): 182-201. |
| [16] | YANG Shuwen, XUE Chongsheng, LIU Tao, et al. A Method of Small Water Information Automatic Extraction from TM Remote Sensing Images[J]. Acta Geodaetica et Cartographica Sinica, 2010, 39(6): 611-616. (杨树文, 薛重生, 刘涛, 等. 一种利用TM影像自动提取细小水体的方法 [J]. 测绘学报, 2010, 39(6): 611-616.) |
| [17] | CHEN W B, FUKUI H, DOKO T, et al. Improvement of Glacial Lakes Detection under Shadow Environment Using ASTER Data in Himalayas, Nepal[J]. Chinese Geographical Science, 2013, 23(2): 217-226. |
| [18] | LIANG Dong, KONG Jie, HU Gensheng, et al. The Removal of Thick Cloud and Cloud Shadow of Remote Sensing Image Based on Support Vector Machine [J]. Acta Geodaetica et Cartographica Sinica, 2012, 41(2): 225-231, 238.(梁栋,孔颉, 胡根生, 等. 基于支持向量机的遥感影像厚云及云阴影去除[J]. 测绘学报, 2012, 41(2): 225-231, 238.) |
| [19] | CHENG Wei, WANG Liming, TIAN Qingjiu. A Method of Atmospheric Correction Based on Shadow-pixel for Optical Satellite Data[J]. Acta Geodaetica et Cartographica Sinica, 2008, 37(4): 469-475. (程伟, 王黎明, 田庆久. 一种基于阴影像元的光学遥感大气校正方法[J]. 测绘学报, 2008, 37(4): 469-475.) |
| [20] | ZHONG Zhenwei, YE Qinghua. Integrated Method of Glacial Information Extraction around the Mount Naimonacnyi[J]. Journal of Glaciology and Geocryology, 2009, 13(4): 717-723. (仲振维, 叶庆华. 纳木那尼峰地区冰川信息的综合提取方法[J]. 冰川冻土, 2009, 13(4): 717-723.) |