高原冰川表面温度反演研究——以祁连山老虎沟12号冰川为例

余寒，阳昌霞，刘汉湖，曾敏

(成都理工大学 国土资源部地学空间信息技术重点实验室，成都 610059)

冰川温度是冰川最基本的物理特征[1]，而冰川表面温度反映了冰川表面与外界能量交换，从而影响着冰川运动和变化[2-5]。19世纪，国外学者开始开展对冰川温度的研究，20世纪中期中国开始对冰川温度研究[1]。目前，对冰川温度的研究多数是实测[6-9]，耗时耗力且获取的数据有限，遥感数据具有更新周期快，易获取等优点，被广泛应用于城市和海洋表面温度的研究中[10-11]，如金旭峰等利用AMSR2遥感数据对地表温度进行了反演[12]，郑贵洲等利用MODIS遥感数据对海洋表面温度进行了反演[13]，但对于冰川表面温度的遥感反演研究甚少。

目前地表温度反演算法主要有单通道算法、劈窗算法和多通道算法三大类，运用比较成熟的是单通道算法，具有代表性的有辐射传输方程、Jiménez-Muoz(JM)模型算法[14]和覃志豪单窗算法[15]。其中JM算法只需求取地表比辐射率和大气含水量便可得到地表温度，算法简单，参数少，被广泛应用于温度反演，而Landsat系列卫星参数易获取，是良好的温度反演源数据。本研究利用Landsat 8遥感数据反演了祁连山脉老虎沟12号冰川的表面温度，并进行了冰川温度分布及规律的相关性分析。

1 研究区概况

老虎沟12号冰川(39°25′N～39°30′N，96°30′E～96°34′E，图1) 位于青藏高原祁连山西段北坡，面积约为20.42 km2，长约9.7 km，海拔跨度约1 200 m，平均海拔4 943 m，是祁连山最大的山谷多温型冰川，降水集中在5—9月[16-18]。老虎沟12号冰川分东西两支，汇合于4 550 m的消融区。过去50年间，老虎沟12号冰川萎缩显著，冰川长度减少了403 m( 总长度的4%) ，面积减少了1.54 km2( 总面积的7%)[19]，平衡线高度也在不断增加，冰川表面流速也较1959 年减慢了11%[20]。此外，老虎沟12 号冰川东、西支的平均厚度分别为190 m 和150 m，冰川槽谷形态具有空间差异，东、西支冰川槽谷形态近似于对称的V型,4 550 m汇合区中部冰川厚度增加到162 m，槽谷底部变宽，边坡变缓，发育有不对称槽谷[21-22]。

图 1 研究区位置示意及RGB影像(B2-B4)Figure 1 Location of study area, RGB images (B2-B4)

表 1 遥感影像信息 Table 1 Image information

表 2 大气参数 Table 2 Atmospheric parameters

2 数据与方法

2.1 数据来源

数据来源于中科院Landsat下载共享系统(http://ids.ceode.ac.cn/query.html)。由于研究区可利用的Landsat 8数据较少，因此选取了近年来6期过境老虎沟12号冰川云量较少的遥感影像，多光谱波段OLI分辨率为30 m，热红外波段TIRS分辨率为100 m，具体信息见表1。

2.2 Jiménez-Muoz(JM)模型算法

JM算法是通过一个热红外波段来计算地表温度单通道算法[14]，它由辐射传输方程演变而来，简化后具体形式如下式：

式(1)中:TS为地表温度，单位为K；ε为地表比辐射率；L0为热红外波段地表辐射亮度值，由ENVI软件Radiometric Calibration工具可计算得到。τ为大气透过率;L↑为上行辐射率;L↓为下行辐射率，在NASA公布的网站(https://atmcorr.gsfc.nasa.gov/)输入成像时间及中心纬度，即可查询。本研究利用的大气参数见表2，αT0、βT0为参数。

2.2.1地表比辐射率

地表比辐射率又称发射率，是指在同一温度下地表发射的辐射量与一黑体发射的辐射量的比值，其计算方法是由多光谱数据求得NDVI，再由下式求得：

(2)

ε=0.004×Pv+0.986

(3)

式(2～3)中:Pv为植被覆盖度；NDVISoil为完全是裸土或无植被覆盖区域的NDVI值；NDVIVeg为完全被植被所覆盖的像元的NDVI值，即纯植被像元的NDVI值，取经验值NDVIVeg= 0.70和NDVISoil= 0.05，即当某个像元的NDVI大于0.70时，Pv取值为1；当NDVI小于0.05，Pv取值为0。

2.2.2 冰川表面温度计算过程

由普朗克定律可得：

(4)

BT0=[L0-L↑-τ(1-ε)L↓]/(τ×ε)

(5)

(6)

式(4～6)中:T0为初始亮度温度值；BT0为初始亮度温度值下黑体辐射值；BTs为在温度等于TS时黑体辐射值；λ为波长(μm)。k1、k2为定标常数；c1、c2为普朗克辐射常数；数值在表3中给出。

由式(4～6)很难求出地表温度TS，但由高斯近似展开可得到下式：

表 3 波段参数值 Table 3 Band parameter value

(7)

(8)

(9)

将式(8～9)代入式(1)，即可得到温度反演结果TS。

2.3 相关性分析方法

本研究采用了线性回归方法研究了冰川表面温度与高程和坡向之间的关系。线性回归是利用数理统计中的回归分析来确定两种或两种以上变量间相互依赖的定量关系的一种统计分析方法，运用十分广泛[22]。采用显著性水平指标来评判高程和坡向与冰川表面温度之间的相关性。

3 结果与分析

3.1 温度反演结果与验证

利用简化后的JM算法对老虎沟12号冰川表面温度进行反演，得到研究区冰川表面温度分布(图2)。从表4得知，平均温度与2013—2015年的同时期夏季气温[22]相比差别不大，说明此方法可应用于高原地区冰川表面温度反演。表5列出了6年来老虎沟12号冰川消融区(4 450 m)，冰川多年平衡线(4 900 m)，冰川积累区(5 040 m)的反演结果和2010年、2011年夏季实测冰温[1]，虽然时期不同，但从温度上升趋势来看，反演结果与实测温度差别不大，进一步说明了该方法可用于冰川表面温度反演，但由于各种误差的存在，依然要对算法作进一步的改进。

表 4 冰川表面温度统计 Table 4 Glacier surface temperature statistics

表 5 反演结果与实测冰温比较 Table 5 Comparison of retrieval results with measured ice temperature

3.2 温度时空变化

根据图2,2013年8月11日冰川温度范围在270.13～280.98 K之间，低温主要分布在南部，北部和边界温度较高。2014年8月14日冰川温度范围在266.91～280.61 K之间，但相较于2013年，整体温度有所升高，低温依旧分布于南部，北部温度有明显升高。2015年8月17日冰川温度范围为268.48～281.93 K，整体变化相较于2014年同时期无明显变化，2016年7月2日冰川温度范围为273.16～283.05 K，最低温度较上一年有所上升。之所以从图2中看温度大部分呈现为蓝色，是因为温度极差最小，且大部分温度为273.15 K左右。2017年7月5日冰川温度范围为271.45～282.98 K，最低温度较上一年有所降低，但整体温度略微上升，2019年9月29日冰川温度范围为266.04～280.05 K，因其时间进入秋季，整体比前一年有所降低，平均温度为272.25 K。根据反演结果可以看出，研究区冰川表面温度主要集中在270～278 K之间，研究区不同区域温差较大，最大达到13.45 K，一般情况下，冰川表面温度小于273.15 K，但由于冰雪在融化过程中吸收热量并且流动速率加快[24]，导致冰川表面温度可以上升到273.15 K以上。值得注意的是南部边界冰川表面温度没有明显的升高，因为南部边界与其他冰川相邻，具体原因见讨论部分。

3.3 高程相关性分析

冰川表面温度受多种自然条件因素的影响，本研究主要从高程和坡向两个自然条件因素进行了相关性分析。

从6年的回归方程和显著性水平R2来看，冰川表面温度总体上随高程的升高而降低，且反演结果大致符合高程每升高1 000 m，温度下降6 K的规律(图3)。其中显著性水平最高的年份是2015年，达到了0.914 6，显著性水平最低的年份是2016年，数值为0.635 2，总体显著性水平较高，说明此反演方法在高原地区同样适用。2016和2017年温度整体高于2015、2014、2013年，因为前两者数据选取于7月，后三者数据选取于8月，太阳直射纬度向南移动，所以7月接受的太阳辐射相对较多，冰雪消融面积相对较大，冰川表面径流速度加快，8月接受的太阳辐射变少，因此7月的冰川表面温度相对较高，而8月相对较低。2019年9月因为季节的原因，温度整体偏低。图3中某些高程存在温度异常增高现象，原因或为混合像元或者冰雪类型不同。

3.4 坡向相关性分析

影响不同坡向温度差异的首要因素是太阳辐射[25]。从图4可以看出，由于老虎沟12号冰川88%左右的区域坡向向北，所以温度绝大多数分布在0°～90°和270°～360°之间，且受卫星过境时间和太阳辐射的影响，东北坡向(22.5°～67.5°)温度略高于西北坡向(292.5°～337.5°)温度，东南坡向(112.5°～157.5°)温度略高于西南坡向(202.5°～247.5°)温度，东部温度分布区间较西部窄。由图4和图5可知，2013年至2017年100°～250°坡向温度回归曲线有一个略微向下凹的趋势，因为在这个坡向范围之内的大部分区域位于高海拔地区，因此平均气温较低。而2019年此处曲线向上凸，表明进入秋冬季节后，由于坡向带来的太阳辐射影响减弱，东西坡向已没有明显的温度差别。

图 3 相同高程平均温度分布Figure 3 Average temperature distribution at the same elevation

图 4 不同坡向平均温度分布Figure 4 Average temperature distribution in different aspect

图 5 高程与坡向Figure 5 Elevation and aspect

4 讨论与结论

4.1 讨论

由于各种各样的原因，温度反演结果存在误差及异常，大约在-1～5 K左右，位于可接受范围内[12]。究其原因，主要表现在以下几个方面：

1)Landsat 8多光谱波段影像分辨率为30 m，热红外波段分辨率为100 m，两个波段融合计算存在一定误差，且由于不同冰雪类型的比辐射率不同，因此不同冰雪类型的混合像元在通过NDVI计算比辐射率时，会产生一定误差。

2)和其他区域相比，冰川边界和冰舌区域误差较大，这是由于冰与岩石等构成了混合像元，导致像元的辐射亮度值被高估，因此使得冰川边界部分温度反演结果偏高[26]，其中图3显示，在高程5 130 m左右有一个温度异常，原因就是因为此处岩石与冰川构成了混合像元。但是，北部边界冰川表面温度反演结果没有明显的升高，因为北部边界与其他冰川相邻，构成了连续冰川，而不是与其他地物的混合像元。同时，本研究为了方便比较，采用了相同的边界数据，但边界可能随时间会发生一定变化，在边界区域出现其他地物也会导致温度反演结果的偏差。

3)冰川温度受外界影响因素较大，会随时间变化，且周围山脉众多，太阳照射产生的阴影也会导致温度反演结果的偏差。另外，对于地表覆盖的冰雪类型未作细致划分，也可能产生一定误差。

4.2 结论

本研究利用多期Landsat 8多光谱和热红外遥感数据以及ENVI和ArcGIS软件，通过JM单通道算法对老虎沟12号冰川表面温度进行了反演，分析了不同区域的温度分布情况，并从高程、坡向对反演结果进行了相关性分析，得到以下主要结论：

1)从温度分布上看，冰川表面温度总体上随海拔升高而减小，冰川北部边界区域和冰舌等低海拔区域温度反演结果较高，低温区域主要分布在冰川南部等高海拔区域。

2)冰川表面温度反演结果范围在265～285 K之间，平均温度在274～277 K之间，冰川不同区域表面温度温差较大，最大达到了13.45 K，标准差在2 K左右，温度呈逐年上升趋势，与文献[23，27]具有一致性。

3)冰川表面温度受到坡向的影响，东北坡向温度略高于西北坡向，东南坡向温度略高于西南坡向，东部温度分布区间较西部窄。