面临全球变化、气候变暖的威胁，全球冰川基本上呈现为大范围的消退状态［1］。山地冰川位于低纬度、高海拔的地区，与南北极和格陵兰岛冰川和冰盖相比，它们对于气候的变化更加敏感;因此对于山地冰川的监测是预测其发展变化的关键所在。监测冰川表面的流动速度也是全球气候变化研究的重点问题之一。冰川物质平衡变化是对大气条件的直接、快速响应。由于冰川物质的交换主要依赖于冰川流的速度，所以冰川表面速度的估计在研究冰川物质平衡起到了重要的作用，冰川表面速度的精确测量是研究和理解冰川动力学的必要参数之一［2］。

长期以来，冰川实地观测测量是冰川研究的主要方法，但是该方法比较浪费时间，成本高，是一种比较有难度的记录冰川演变信息的方法;而且仅有极少数的较为容易抵达的冰川才采用这种方法进行测量监测。虽然对于冰川速度的实地观测测量可以达到很高的精度，但是，在实际应用中，采用这种方法不容易实现对每个冰川的常态化测量。而且，大多数的山地冰川分布在高山地区，地形条件复杂，交通条件不便，很多地区难以到达。另外，冰川区气候条件变化强烈，气温常年在0 ℃以下，这些因素都给实地观测测量带来了较大的困难。然而，遥感技术为自动监测整个冰川系统提供了一个简单有效的方式，该技术克服了实地观测测量方式存在的缺点和不足。通过搭载在卫星平台上的传感器观测山地冰川对研究山地冰川是很有帮助的，因此光学遥感和微波遥感是一种冰川表面速度获取和研究的有效手段。

光学遥感方式在监测山地冰川变化中存在一些限制。光学数据的质量严重依赖于太阳光照情况。但是，大多数被冰川覆盖的地区都有一个相似的特点，即恶劣的天气条件，比如经常有雨雪等降水和大风天气，这些天气状况都导致了光照条件的变化，而且给野外工作带来了困难。另外，光学遥感数据的质量受到云层覆盖的严重影响，搜集在冰川区域成像的无云覆盖数据是研究和模拟冰川流动模式的必要前提，然而在实际应用中，这并非是一项容易的任务。而且，覆盖冰川区的数据经常缺失，造成冰川观测的不连贯，也使得在高海拔和偏远地区提取冰川速度场困难重重。与光学数据相比，合成孔径雷达(synthetic aperture radar，SAR)不受限于天气和光照条件，可以在任何天气条件下定期获取研究区的雷达图像数据。

雷达传感器可以同时获取强度和相位信息，所以合成孔径雷达数据拥有更有效的可用性。特别地，基于相位信息的差分干涉合成孔径雷达技术可以获取厘米级的形变量，该方法应用在冰川学中得到了一些引人注目的结果［3，4］。利用两幅在不同时间和微小差异的轨道条件下对同一地区成像的复数雷达图像进行复共轭相乘，提取两次得到的信号的相位信息的不同。相位信息不仅对地表地形敏感，而且也对在两次成像时间间隔内发生在雷达视线向的位移同样敏感。利用两幅干涉相位图差分处理，去除干涉条纹图中与地形相干的相位信息，从而得到由于地表形变导致的相位变化信息并反演得到地表形变图。然而，只有少数研究成功地利用差分干涉技术获取了山地冰川形变量。强烈的地形影响和山地冰川较小的面积是利用该技术的障碍。但是，更为重要的影响因素是两次连续获取图像之间的时间间隔大小。如果超过了1 d 或者3 d，冰川表面的形变梯度超过了干涉技术的梯度阈值，使得在时间间隔内冰川表面的相干性受到了严重的破坏，无法获取与形变相干的相位信息［3，5］。只有欧洲遥感卫星ERS-1/2 冰川模式和Tandem 任务可以用于提取冰川表面的速度场分布。随着新的雷达卫星的发射，如Radarsat-2、TerraSAR、Cosme-skype、雷达数据拥有了更高的空间分辨率和时间分辨率，较多的采用了X 波段和C 波段，使得对冰川表面积雪的穿透性降低，图像特征的稳定性降低，而且商业SAR 卫星的数据成本较高。

除此之外，差分雷达干涉技术的成功应用同样受限于相位噪声，其质量一般可以反映在相干性信息中。如果没有在整幅图的配准中考虑到冰川区局部的形变因素影响，得到的冰川区的相干性就会受到破坏。在干涉处理中，整个冰川的相干性随着两次数据获取时间的间隔增加而减少。在间隔期内，由于冰雪的融化，降雪和大风影响，冰川表面的冰雪的分布发生了变化，导致了冰川表面的失相干。而且冰川相对于雷达视线向的方向在山地环境中也至关重要。理想境况下，冰川位于背向雷达卫星的斜坡，而且是平行于雷达视线向滑动。另外，干涉测量技术得到的相位结果是缠绕的，相位解缠是必要的一步。然而，强度跟踪算法可以直接获取方位向和距离向确定的形变分量，可以参考相对稳定且表面没有冰川覆盖的不动区域来得到冰川移动的绝对量，因此长时间间隔的雷达图像依然可以使用该方法。尽管利用ERS 和Envisat 卫星有些成功监测冰川的例子，但是自从1990年开始运行以来，这些卫星获取了大量的雷达数据。不幸的是，在一些区域，如被浓密植被或雪覆盖，C 波段干涉往往不能成功应用。另外一个比较合适的选择是使用搭载在高级陆地观测卫星上的相阵式L 波段雷达传感器获取的L 波段的SAR 数据。与C 波段相比，L 波段微波具有更强的穿透性，可以到达地表以下更深的部分，所以在L 波段图像中各种特征能在相对较长的时间内保持稳定，不受表面临时降雪的影响。基于L波段数据的强度跟踪算法相比C 波段得到的结果精度更高［6］。

由于快速流动和其他气象条件造成雷达图像对之间失相干，所以相比干涉技术，基于雷达图像的偏移跟踪技术是另外一个颇受欢迎的用于估计冰川表面速度场的方法。两种不同的技术用于计算影像的偏移，分别是强度跟踪和相干跟踪，这两种方法在文献［7］进行了详细的介绍。强度跟踪算法最大的优势在于它可以在不管雷达图像是否相干都可以应用。在很多区域，对于长时间间隔获取的图像对成功分析都是利用且仅限于强度跟踪算法，因为相干性在冰川表面不能继续保持超过几天的时间间隔。

冰川表面的特征，例如裂缝、冰碛物等，它们在冰川表面保持相对稳定，并随着冰川一起流动。冰川表面特征的变化往往可以通过配准两幅图像得到，并对比共同特征，自动追踪移动的特征从而得到冰川流速。基于强度的位移估计方法是通过获取得到图像特征的偏移量，拟合去掉整体形变，从而得到位移估计。基于强度的方法提供了另外一个可以完全利用在干涉处理中因为失相干而无法使用的数据。

本文利用1999年1月份和3月份的ALOS/PALSAR 数据和基于图像强度信息的归一化互相关算法，获取慕士塔格峰地区间隔44 d 的冰川速度分布图。

1 研究区与数据

慕士塔格峰位于新疆维吾尔族自治区塔里木盆地西部边缘，地理位置为北纬38°15'，东经75°8'，地处帕米尔高原的东南缘和青藏高原的西北缘，是两大高原的连接处，平均海拔高度4 000 m 以上，其主峰高度为7 546 m［8］。该区常年受西风带影响，气候寒冷，终年以固体降水为主，因此非常有利于冰川发育。

慕士塔格峰是一座呈浑圆形的断块山，围绕其形成了很多山地冰川，它们以慕士塔格峰冰原为中心，呈放射状向外分布，见图1。该地区有现代冰川128 处，冰川总面积为377.21 km2，面积超过10 km2的有8 条，其中面积约为85 km2 的科克萨依冰川是该地区最大的冰川，长度约有15 km，宽度约为1.5 km。

本文利用日本高级陆地观测卫星携带的相控阵型L 波段合成孔径雷达(advanced land observing satellite，phase array L-band synthetic aperture radar，ALOS/PALSAR)获取的雷达图像对该区域进行研究，考虑到该地区地形复杂，高差较大，容易使得雷达图像产生较大的几何畸变，如顶底倒置和阴影等，从而影响雷达图像质量，使得基于雷达图像的结果被噪声淹没，无法得到与冰川区运动相关的有效信号。考虑到以上原因，在数据选取阶段重点考虑了雷达影像对的空间基线大小，最终选择了水平和垂直基线均小于300 m 的一对PALSAR 数据作为研究的基础数据，图像对的参数特征如表1 所示。

图1 墓士塔格峰地区Landsat 743 波段组合图像(叠加了等高线，小图为该区域位置示意图)
Fig.1 Image of Landsat 743 band combination of Muztag Ata adding the contours(the small map on the upper left shows the location of the region in China)

表1 匹配影像对的参数特征
Table 1 Characteristics of the parameters of the matching images

2 方法与结果

采用了基于雷达图像强度的归一化互相关算法来提取与冰川移动相关的像素偏移信息，它是一种基于相似性的量度，通过获取子模块与搜索窗的最大相似度所对应的欧几里的距离［9］，反演得到研究区冰川流动和变化的一种有效手段。基于归一化互相关算法的图像跟踪方法可以得到冰川区雷达图像在方位向和距离向的偏移。在预配准完成后，便可以在雷达图像对之间使用该方法得到两个方向的偏移量，通过矢量合成最后得到每一采样点处总的位移大小。

对归一化的互相关算法进行一个简单的介绍，如公式(1)所示，进行归一化的目的主要是消除不同图像亮度和对比度对配准结果的影响［10］，使该算法具有更好的适用性和更强的鲁棒性。

式(1)中γ 是相关性系数，f 是在辅图像中的搜索窗口，t 是在主图像中与之相对应的模板窗口。搜索窗口f 的大小通常大于模板窗口t 的大小。μ 和ν分别是两个图像子块在x 和y 方向的偏移量。模块窗口和搜索窗口的大小选择对最后偏移结果的有着重要的影响，采用的模块窗口大小应该能够保证较大的信噪比，而且可以反映出速度梯度较小变化［11］。搜索窗口大小应该足够大从而可以包括偏移较大的模块窗口，但同时也要考虑计算效率，所以搜索窗口也不能太大，否则会造成数据处理效率低下［9］。

通过计算主图像中t 窗口和辅图像中的搜索窗口f 的归一化相关系数，从而在辅图像中找到与t 模板窗口相对应的图像窗口。通过式(1)计算得到归一化的互相关系数γ，并在搜索窗口中每次偏移一个像素循环计算得到相关系数γ 组成的相关矩阵，其峰值γ(u，v)的位置(u，v)对应着两个窗口图像的偏移量。峰值的大小表示两个图像窗口的相关性大小，其取值分为［0，1］，γ 为0 时表示两幅窗口图像完全失相关，γ 为1 时表示两幅窗口图像完全相关。图像之间即使没有真实相对应的地物，该方法也会得到一定的相关峰值。在实际应用当中，不同时间的成像条件会发生一定的变化，相关性较强的图像块之间得到的相关系数峰值明显，而相关性较差的图像窗口之间的相关系数矩阵峰值不够突出，见图2，有时由于噪声的影响也可能导致图像窗口之间的失配，因此，必须对γ 设置一个阈值，认为小于该阈值的γ 对应的偏移量(u，v)为无效值。在本文中，阈值设为0.3。

图2 归一化互相关系数矩阵
Fig.2 Normalized cross-correlation coefficient matrix

利用上述归一化的互相关算法只能得到像素级的图像配准精度，然而图像窗口匹配块之间的偏移不一定都刚好是整数值。为了获取亚像素级的精度的位置偏移量，通常采用两种不同的方法提高配准精度:一是对相关图像窗口进行插值处理，另外一种是拟合二维相关矩阵并计算其峰值点对应的坐标值。本文主要是通过快速傅里叶变换插值方法对图像进行插值处理，从而得到图像窗口配准所需的亚像素级精度。

3 结果与分析

首先利用ROI_PAC 软件包对ALOS/PALSAR原始信号数据进行成像处理，得到研究区的单视复图像(SLC)，其方位向分辨率约为3.5 m，距离向分辨率约为4.7 m。接着基于强度的互相关算法被用来精确估计主辅图像之间的亚像素级偏移量，从而准确地配准主辅图像。根据经验，选择了128×128大小的主图像模块进行粗配准，32×32 大小的主图像模块进行精配准，相应的辅图像搜索窗口大小分别设定为256×256 和64×64，这样既满足了获取最大位移的需要，同时又降低了计算时间，提高了数据处理效率，而且粗配准得到的偏移量作为精配准操作时的初始值，有利于提高配准结果的准确度。

由于精确配准值对应的γ 值在0.2～0.8 之间，均值为0.45，另外一个有效的配准质量评价量信噪比SNR 也同时进行了计算，其为相关矩阵γ 峰值与非峰值的均值的比值，它也可以作为判断配准质量好坏的一个重要参数，在本文中，经验性的将其阈值设置为5，信噪比在3 到90 之间变动，均值约为20。因此，利用相关系数大于0.3，信噪比大约5 所对应的偏移量进行最小二乘拟合，去除由于卫星轨道、地形和卫星姿态变动导致的部分位移量，得到与冰川运动相关的位移信号。

通过上述方法，分别提取到了研究区雷达图像中距离向和方位向的偏移量，见图3。由图可知，慕士塔格峰地区相对较大的冰川在结果图中均有与其冰川速度相关的信号反应。图像中部颜色较深、变化比较杂乱的区域为噪声或者地形影响造成的失相关区，无法提取到有效的偏移量值，因此无法获取慕士塔格峰峰顶附近的位移量。在慕士塔格峰的右侧，主要分布着一条大型冰川，即该区域最大的冰川—科克萨依冰川，图3 中比较详细的反应出了该区域沿距离向和方位向的变化情况:在距离向，冰川整体上朝远离主峰的方向移动;在方位向，由于冰川运动相对复杂，但是图像中依然可以有效地反应出其运动特点，同时说明该方法的有效性。在主峰左侧，分布了较多相对科克萨依冰川较小的冰川，它们方位向和距离向的运动特征在图像中也有较好的表现。值得注意的是，在距离向位移图中，位于主峰上下部的冰川流动方向接近于方位向，因此与它们相关的信号较弱，但是在方位位移图中可以看出它们的信号反应较强。

最后利用得到的距离向和方位向的速度位移图，通过矢量合成的方法，获取到了该区域冰川速度大小分布图，见图4。非冰川区由较浅的蓝色覆盖，其表示的数值大多在0.5 m 以下，通过统计这些非运动区域的偏移量残差，可以进行本文所用方法的精度分析，结果表明，结果精度约为0.5 m。在平坦区域误差较小，地形变化较大的区域误差相对较大，因为地形对偏移量提取结果和最后整体形变的去除都有比较大的影响。冰川整体趋势呈现为，底部区域冰川位移量较小，随着海拔的升高，冰川位移矢量会进一步增加。在坡度比较小的冰川区，其相应速度也比较小，在坡度较大的地区其冰川流速也较大，速度大小与地形表现出了比较强的相关性。通过图4 可以得到，在数据获取间隔时间内，该地区冰川速度大小基本在7 m 以下，部分区域由于较大的地形坡度，速度大小到达了8～10 m，甚至超过了10 m，如图4 中A 方框所示区域。提取到了位于主峰右侧的科克萨依冰川表面的运动特征，除了靠近主峰的部分区域失相关外，其余区域的运动特征细致的体现在了图4 中的B 区域。该区域的最大值出现在科克萨依上部冰川的一个分支上，大小约8 m。

图3 慕士塔格峰地区冰川位移分布图
Fig.3 Map of glacier displacement distribution:slant-range displacement

慕士塔格峰地区冰川整体上表现为主峰左侧速度变化较大，主峰右侧科克萨依冰川速度变化相对较小。整体误差在0.5 m 以下，但是其大小分布与地形变化相关，因此由于高差起伏对最终结果造成了较强的噪声。但是由于SAR 的成像几何条件以及缺乏该地区的高精度的DEM，使得对雷达图像地形校正难以实现，因此获取的大小变化是在雷达斜视坐标系统内，最终的平面二维形变度量还需要借助外部DEM 和雷达成像参数完成，在本文中对此并没有涉及，进一步地研究工作还在进行中。

图4 慕士塔格峰地区冰川流速大小分布图(m)
Fig.4 Map of glacier surface flow velocity field(m)

4 总结

利用2009年1月和3月获取的ALOS/PALSAR雷达图像和基于强度信息的归一化的互相关算法成功地提取到了慕士塔格峰地区冰川流动所对应偏移量，并得到了令人满意的效果，结果精度达到了1/10个像素。对比主峰东西两侧的冰川变化，发现，在观测期内，左侧冰川流速明显大于右侧冰川。本文所述方法可以应用在雷达干涉处理无法得到有效结果的情况下，从而为更好地更有效地利用类似数据提供了一种补充手段，提高了存档数据的利用率。总之，基于强度信息的归一化互相关算法为冰川动力学研究提供了一种有效的工具，特别是在冰川流速过快，SAR 数据获取时间间隔较长，干涉SAR 方法失效的情况下，该方法可能是唯一一种能够应用在长时间基线SAR 数据对中，获取山地冰川流速信息的手段。

参考文献

1 Magnusson J，Jonas T，Lopez-Moreno I，et al.Snow cover response to climate change in a high alpine and half-glacierized basin in switzerland.Hydrology Research，2010;41 (3—4):230—240.doi:10.2166/nh.2010.115

2 Pedersen V K，Egholm D L，Nielsen S B.Alpine glacial topography and the rate of rock column uplift:a global perspective.Geomorphology，2010;122(1—2)，129—139

3 Goldstein R M，Engelhardt H，Kamb B，et al.Satellite radar interferometry for monitoring ice-sheet motion，application to an antarctic ice stream.Science，1993;262(5139):1525—1530

4 Joghin I，Winebrenner D，Fahnestock M，et al.Measurement of icesheet topography using satellite-radar interferometry.Journal of Glaciology，1996;42(140):10—22

5 Gourmelen N，Kim S W，Shepherd A，et al.Ice velocity determined using conventional and multiple-aperture inSAR.Earth and Planetary Science Letters，2011;307(1—2):156—160

6 Nakamura K，Wakabayashi K，Doi K，et al.Ice flow estimation of shirase glacier by using JERS-1/SAR image correlation.Geoscience and Remote Sensing Symposium，Barcelona，2007:4213—4216

7 Strozzi T，Luckman A，Murray T，et al.Glacier motion estimation using SAR offset-tracking procedures.Geoscience and Remote Sensing，2002;40:2384—2391

8 蔡迪花，马金辉，年雁云，等.慕士塔格峰冰川变化遥感研究.兰州大学学报，2006;42(1):13—16 Cai D H，Ma J H，Nian Y Y，et al.The study of glacier change using remote sensing in Mt Muztag Ata.Journal of Lanzhou University，2006;42(1):13—16

9 Gilo M，Kääb A.Sub-pixel precision image matching for measuring surface displacements on mass movements using normalized cross-correlation.Remote Sensing of Environment，2011;15(1):130—142

10 Lewis J P.Fast normalized cross-correlation.Vision Interface，1995;10(1):120—123

11 Kanade T，Okutomi M.A stereo matching algorithm with an adaptive window，theory and experiment.IEEE Transactions on Pattern Analysis and Machine Intelligence，1994;16(9):920—932