青藏高原多年冻土区不同高寒草地类型地表形变特征

王志伟,岳广阳,吴晓东

1 贵州省农业科学院草业研究所,贵阳 550006 2 中国科学院西北生态环境资源研究院 冰冻圈科学国家重点实验室 青藏高原冰冻圈观测研究站, 兰州 730020 3 美国德克萨斯州大学圣安东尼奥分校地质学系,美国圣安东尼奥 78249

青藏高原作为“世界第三极”,不仅在亚洲季风系统中发挥着重要作用,对全球气候变化也异常敏感[1- 2]。它不仅是我国气候系统的重要组成因子,还影响到全球尺度的气候变化[3- 4]。被称为“气候变化指示器”的多年冻土[5]广泛分布于青藏高原[6],面积约为1.06×106km2[7]。多年冻土指温度能够维持在零摄氏度以下状态两年及两年以上的近地表土壤或岩石层,活动层是多年冻土区位于多年冻土之上,夏季融化、冬季冻结的地表部分。伴随全球气候变暖[8],多年冻土在逐年退化[9],活动层厚度相应增加[10],影响到土壤特别是近地表的水、热循环过程和生态环境[11]。多年冻土的热融效应增长会促使大团聚体破碎成小团聚体,释放大量有机碳、硝态氮等物质,进而造成地表植被发生改变,影响到地表的一系列特征,如反照率、降水的渗透速度、土壤中的蒸腾和蒸散、以及土壤侵蚀等,从而打乱水文和气候系统的循环速率,造成高寒地区的生态环境恶化。

青藏高原多年冻土区蕴含世界上海拔最高、类型最为独特的高寒草地生态系统[12],其分布面积分别占青藏高原多年冻土区和全国草地总面积的80%[13]和38%[14],具有防风固沙、涵养水源、固氮储碳、调节碳循环及气候变化、维护生物多样性等诸多生态服务功能[15- 16],是国家生态安全的重要屏障[14]。已有研究表明,高纬度和高海拔地区的生态系统对气候变暖的响应更加敏感[17]。而青藏高原的高寒草地作为高原气候和冻土环境双重属性的生态系统,其对全球气候变化的响应更迅速、更超前[18],其生态系统的稳定性也更脆弱[19]。综上所述,对青藏高原多年冻土区高寒草地生态系统地表环境变化的特征分析极为重要,可以为揭示高寒草地生态系统在全球变化中的生态价值和贡献提供重要的理论基础和科学依据[12]。

随着星载合成孔径雷达干涉测量技术(Interferometric Synthetic Aperture Radar, InSAR)的发展[20],可以实现对多年冻土区近地表冻融循环变化特征的大范围、高精度和高分辨率监测[21],为寒区生态环境的安全、平稳和可持续发展提供了十分有效的监测工具和技术手段[22]。近些年利用InSAR技术分析多年冻土区的相关研究日益成熟[23- 24],Liu等[25]在阿拉斯加多年冻土区利用ERS 1/2 SAR数据成功探测到厘米级的地表沉降量,之后[26]进一步反演该区域的活动层厚度。Chen等[21]使用ALOS PALSAR和Envisat ASAR数据,获得青藏高原北麓河的地表形变量,结果显示部分区域每年的变化量甚至为±2 cm。Zhao等利用Envisat ASAR数据[24]对青藏高原当雄至羊八井的冻土区地表形变进行研究,发现自然地表的形变量在3.6—5.0 cm,而铁路和公路的形变量则在2.8—3.7 cm。上述研究多只是单纯针对地表形变(Ground surface deformation, GSD)和活动层厚度进行探讨[26],而发掘不同高寒草地生态环境条件下地表形变特征变化的研究还较为薄弱[27]。鉴于此,本研究以青藏高原多年冻土区五道梁区域的不同类型高寒草地为研究对象,通过分析其地表形变量和遥感生态指数(Remote Sensing Based Ecological Index, RSEI)的关系(本研究中的生态指数包括绿度、湿度和热度三项[28]),试图揭示：(1)不同高寒草地类型条件下的地表形变量变化特征如何？(2)地表形变量变化与绿度、湿度和热度指标存在何种关系？旨在对不同高寒草地类型地表形变特征进行挖掘,以探讨青藏高原多年冻土区不同高寒草地生态系统对气候变化的响应机制。

1 研究区概况

研究区位于青藏高原腹地、青海省的西部,大片连续多年冻土区(图1),地处北纬34.4°—35.9°,东经92.4°—93.9°之间,平均海拔(4600±190) m,属高原山地气候,夏季降水多,冬季降水少[29]。2005—2010年期间,年均植被指数为0.085±0.037(-),属于植被稀疏区,主要以高寒草地植被类型为主,包括高寒草甸和高寒草原两种。

图1 青藏高原多年冻土分布及研究区位置Fig.1 Permafrost distribution of Qinghai-Tibet Plateau and location of study area

2 数据来源与研究方法

2.1 数据来源及预处理

本研究中涉及的ASAR (Advanced Synthetic Aperture Radar,ASAR)数据有28景。因基线和多普勒质心差的原因,其中4景影像没有参与配对,剩余24景参与运算的影像获取年月日(YYYY-MM-DD)分别为：2005-4-7、2005-9-29、2006-4-27、2006-10-19、2007-5-17、2007-10-4、2008-5-1、2008-10-23、2009-1-1、2009-2-5、2009-3-12、2009-4-16、2009-5-21、2009-6-25、2009-7-30、2009-9-3、2009-10-8、2009-12-17、2010-1-21、2010-2-25、2010-4-1、2010-5-6、2010-6-10和2010-7-15。该ENVISAT ASAR数据集属于单视复数影像(Single Look Complex, SLC),为C波段、IS2模式遥感数据,入射角23°。研究数据获取方式为降轨,数据来源于欧空局(European Space Agency, ESA)。方位向分辨率和距离向分辨率分别为22.6 m和4.05 m。

DEM资料为STRM V4版本数据,该数据在赤道的分辨率约为90 m[24]。

本文使用的MODIS产品有09A1、11A2和13A1三类,每类产品选取Terra(MOD产品)卫星结果。09A1产品为500 米8 天合成数据集,选取其1、2、3、4、6和7六个波段的数值(表1)计算生态指数。地表温度(Land Surface Temperature,LST)11A2和增强型植被指数(Enhanced Vegetation Index,EVI)13A1产品分别为1 km 8天和500 m 16天合成数据集,运算过程中对LST数据执行降尺度处理至500 m。植被类型数据结果从全国植被类型图中提取[30],重采样至500 m。

此外,文中所有地图投影方式都为WGS84坐标系。

表1 MODIS 09A1与TM各波段波长

2.2 地表形变反演方法

地表形变反演,采用短基线集合成孔径雷达干涉测量技术(Small Baseline Subset Interferometric Synthetic Aperture Radar, SBAS-InSAR)。该技术使用的微波数据,不受或少受云、雪和稀疏植被的影响,不需要光照条件,可全天候、全天时地获取。SBAS-InSAR技术能够获取长时间序列的地表形变特征,成熟应用于青藏高原多年冻土区[21, 24]。其原理是通过构建具有较短时-空基线的影像数据对,计算获得干涉图,然后应用奇异值分解(Singular Value Decomposition,SVD)方法生成形变时间序列和平均形变速率结果。SBAS-InSAR干涉时序分析方法计算获取影像形变结果基本计算模型如下式所示：

φint=φdef+φflat+φtopo+φatmo+φnoise

(1)

式中,φint指干涉相位,φdef指地表形变相位,φflat指平地相位,φtopo指地形相位,φatmo指大气延迟相位；φnoise指噪声引起的相位。差分干涉测量主要关注的是地表形变相位,该结果并非借助高程数据(高程信息会作为辅助数据参与运算)反演的高程差,而是相对于第一景影像的形变差(因此第一景影像的形变值为0)。

具体来讲,在研究区内参与运算的影像共有N+1景,获取时间按照顺序排列,依次为t0、t1、…、tN。以其中任意一幅影像为例,该影像会与其他影像进行匹配,生成至少一幅差分干涉图。假设第k幅差分干涉图由ti和tj时刻获取的影像生成,则其计算过程如下式所示：

(2)

获取具有物理意义的形变时间序列结果,其计算过程如下：

假设vk是i到j时刻的平均相位速度,则其相位和时间关系满足如下公式：

vk×(tj-ti)=Δ(φj-φi)

(3)

式中,φ指相位信息。

考虑到不同的时间分段,T0到Tk时刻的第k幅干涉图的相位信息可以通过如下公式计算得到：

(4)

式中,Δφk是影响不同时间间隔的相位速度积分信息,改写成矩阵表达式为

Av=Δφ

(5)

式中,A是系数矩阵；速度v是速度矢量,可以通过系数矩阵A计算获取。因为在SBAS-InSAR处理过程中会有多个主影像,这就有可能导致系数矩阵A秩亏。通常利用SVD方法来处理,首先会生成一个逆矩阵,然后获得速度矢量的最小范数解,最终通过各个时间段内的速度积分获得各个时间段的形变量。

2.3 生态指标计算方法

已有研究[28]指出,遥感生态指数包括绿度、湿度、热度和干度4个指标,分别可以用植被指数NDVI、湿度分量Wet、地表温度LST 和土壤指数NBSI来代表。

考虑到青藏高原多年冻土区植被稀疏的特点[31],本研究中使用MODIS 13A1的EVI产品替代NDVI作为绿度指标。

湿度指标利用TM近似波段的MODIS 09A1产品计算获得,其计算公式如下：

Wet(MODIS)=0.0315b3+0.2021b4+0.3102b1+0.1594b2-0.6806b6- 0.6109b7

(6)

式中,b3、b4、b1、b2、b6和b7分别为MODIS 09A1的第3、4、1、2、6和7波段的反射率,分别对应TM数据的第1、2、3、4、5和7波段数据,如表1所示。

热度指标在此直接利用MODIS 11A2 LST产品。

因研究对象为青藏高原多年冻土区的高寒草地,在此有植被区域不适合用裸土指数。同时,研究区内较少有建筑存在,故在本研究中也不考虑土壤指数。因此本研究中未对干度指标进行计算和分析。

最后,为了后续的统计分析,除了干度指标不做考虑外,其他3个指数都实行归一化(Normalized Differential)操作,具体计算公式如下：

NDRSEI=(RSEIx- RSEImin)/(RSEImax-RSEImin)

(7)

式中,NDRSEI、RSEIx、RSEImin和RSEImax分别对应绿度、湿度和热度3个指标的归一化值、当前值、最小值和最大值,计算后的结果包括归一化绿度(Normalized differential green index,NDGI)、归一化湿度(Normalized differential wet index,NDWI)和归一化热度(Normalized differential land surface temperature index,NDLI)。

3 结果与分析

3.1 研究区地表形变时间序列反演结果

地表形变的长时间序列数据结果主要覆盖2005—2010年时间段,如图2所示。利用SBAS-InSAR方法反演研究区24个时间段(共采集28景影像,4景没有参与配对的影像无结果)的形变结果。其中,于2005年4月7日获取的第一景运算影像被设置为参考影像(即t0时刻影像),该影像所有数值皆为0,故未在图2中展示。此外,图2中红色区域代表地表存在抬升现象；绿色则相反,代表地面的沉降现象。

3.2 研究区地表形变率结果

研究区在2005年至2010年间的地表形变率如图3所示,所有位置的形变率位于15 mm/a之内。图中无数据区域是形变率为0 mm/a和数据质量较差的像元,其余区域地表形变率绝大部分都在±8 mm/a内。如图4所示,其中形变率在±4 mm/a以内的像元占所有形变像元的89.24%,变形率为正的区域占57.70%,地表形变为负的区域占42.30%。

3.3 不同植被类型条件下地表形变结果

研究区24个时间点在不同植被类型条件下的地表形变结果如图5所示,高寒草地整体表现地表下沉的现象,同全球变暖冻土退化进而导致冰融化为水造成土层体积减少的规律一致。

3.4 不同植被类型条件下生态指数结果

利用公式7中的方法,计算归一化生态指数(Remote sensing based ecological indexes,NDRSEI)获取结果如图6所示。图中归一化热度(Normalized differential land surface temperature index,NDLI)体现高寒草甸类的NDLI值大于高寒草原类,但是利用LST计算获取的NDLI结果中并没有当即体现出2008年极端气温低年,而是在2009年3月后开始出现一次极端低地温现象。归一化绿度(Normalized differential green index,NDGI)显示高寒草甸类的NDGI值要大于高寒草原类,特别是在2008年下半年时,存在NDGI显著变小的现象。归一化湿度(Normalized differential wet index,NDWI)不仅较好的体现出湿度季节性变化的特点,而且曲线在2005年4月至2010年7月期间呈现一种湿度稳定增长的状态。

图3 研究区2005年到2010年地表形变率图 Fig.3 Deformation rate of ground surface in the study area from 2005 to 2010 红色为抬升区域,绿色为沉降区域

图4 研究区地表形变率频率分布直方图 Fig.4 Frequency distribution histogram of ground deformation rate in study areas

图5 2005年到2010年间研究区不同高寒草地类型地表形变Fig.5 Ground deformation of different alpine grassland types from 2005 to 2010

图6 2005—2010年间不同高寒草地类型生态指数Fig.6 Ecological indexes of different alpine grassland types from 2005 to 2010

3.5 地表形变与生态指数的相关性

为分析研究区2005年4月至2010年7月3种生态指数对地表形变的贡献与作用,本文利用研究区地表形变(Ground surface deformation,GSD)和3种归一化生态指数,包括NDLI(热度)、NDGI(绿度)和NDWI(湿度)结果,计算获得图7所示的相关图。

图7 高寒草甸和高寒草原条件下地表形变与3种生态遥感参数的相关性Fig.7 Correlations between ground surface deformation and different remote sensing based ecological index in alpine meadow and steppe

4 讨论

4.1 SBAS-InSAR技术的优势及局限性

InSAR技术提供了一种大范围地表形变估算方法,特别是SBAS-InSAR技术尤为适用于青藏高原多年冻土区,同时多种SAR影像传感器的在轨运行,也为长时间序列的地表形变研究提供了必要的数据基础[26]。同时借助现有的植被类型图,甚至可以在青藏高原多年冻土区腹地开展不同植被类型条件下的地表形变特征研究。而高原冻土区腹地一般而言多位于海拔高、人类生存环境恶劣的区域。因此,本文尝试利用现有的SAR影像和植被类型分类结果,利用SBAS-InSAR技术在五道梁地区(青藏高原多年冻土区腹地)开展对不同高寒草地地表形变特征的分析,并得到一些初步结果。

虽然本研究成功的利用InSAR技术、植被类型图集和生态遥感植被对青藏高原多年冻土区进行了地表形变技术特征分析,并得到一些初步的结论,但仍存在以下几方面需要继续加强和改进。

(1)SAR数据过于零散。特别是在图6NDGI曲线中,2007年下半年的点甚至未表现出植被季节性变化的特点。在以后的研究中,需要增加每一年的数据获取量,以此来更准确的反映地表形变的客观规律。

(2)冻土区的研究数据存在时间滞后现象。如图6NDLI曲线中,LST产品对2008年极端低气温年的滞后表现。时间滞后性的研究,在多年冻土区的重要性日益凸显,在下一步的研究中,将会利用合理、有效的模型、算法来挖掘这种规律。

4.2 生态遥感指数对不同高寒草地类型地表形变的影响分析

如图2所示,呈现大规模抬升的区域主要位于研究区东北部,该区域包括3个湖泊,从西至东依次为库塞湖、海丁诺尔湖和盐湖。研究区东北部的大范围抬升可能受区域内湖泊影响,而湖泊区域外,大部分区域表现地面沉降的特点,该结果同当雄至羊八井段的地表形变特征一致[10]。此外,通过图5可知,在2008年附近,地表形变存在短期抬升的现象,可能同2008年我国的极端低温事件发生有关。同时,图5中也展示出高寒草原的地表沉降现象明显强于高寒草甸地区,符合植被条件较好的区域,冻土更加稳定,地下冰融化慢的规律。在图7中,对于高寒草甸而言,地表形变量和3种生态指数的负相关性略大于高寒草原。其中在高寒草甸条件下地表形变与NDLI的负相关性最大,数值为-0.0285,与NDWI和NDGI的负相关性数值依次为-0.0241和-0.0214。而在高寒草原条件下,地表形变量则与NDWI的负相关性最大为-0.0121,然后为NDLI和NDGI,负相关性数值分别为-0.0111和-0.0107。以上研究结果表明,针对不同草地类型,其主导因子存在差异。高寒草甸的地表形变有可能更多的受限于温度变化(NDLI热度指标,由地表温度LST产品计算获得),而高寒草原的地表形变则可能更多的由水分条件所影响(NDWI湿度指标,在遥感生态指数中为表征湿度的分量)。

5 结论

以上研究利用成熟的SBAS-InSAR技术,对青藏高原多年冻土区的五道梁研究区进行了地表形变特征分析,准确的反演出研究区高分辨率的地表形变信息。同时借助生态遥感指标对研究区高寒草甸和高寒草原两种不同的植被类型进行分析,发现高寒草原的地表沉降状况要严重于高寒草甸。以上结果反映出在青藏高原多年冻土区,植被生长越差的区域,其生态系统的脆弱性更为强烈。