主动微波光谱影像对新疆天山山脉区域地表形变的监测研究

王志伟，岳广阳，吴晓东，张 文，王普昶，宋雪莲，吴佳海

1. 中国科学院西北生态环境资源研究院，冰冻圈科学国家重点实验室，青藏高原冰冻圈观测研究站，甘肃 兰州 730020 2. 贵州省农业科学院草业研究所，贵州 贵阳 550006 3. Department of Geological Sciences, University of Texas at San Antonio, San Antonio 78249, USA

引 言

相比可见光、 红外和近红外光谱波段的遥感技术，具有全天时、 全天候优势的合成孔径雷达干涉测量(interferometry synthetic aperture radar, InSAR)方法可以方便、 快捷的实现大范围、 无接触、 面状的，精度为mm级的地表形变监测。该技术方法已成熟地应用于如地下水开采、 矿区沉降、 石油挖掘等人为因素引起[1]的，以及冰川退化、 冻土消融等全球气候变暖等自然因素导致[2]的众多地表变化研究中。特别是随着SBAS-InSAR(small baseline subset InSAR)技术的发展，地表形变的反演精度和准确度得到了进一步提升，促使InSAR技术的应用得到更加广泛的推广[3]。

多年冻土是冰冻圈的重要组成部分，指温度能够维持在零摄氏度以下状态两年及两年以上的近地表土壤或岩石层。近百年来全球气温持续升高，多年冻土也开始逐步退化，众多学者[4]指出这种变暖现象在21世纪可能仍将持续。活动层是多年冻土区夏季融化、 冬季冻结的土壤或岩石层，位于多年冻土之上。多年冻土退化会导致活动层厚度增加，从而改变土壤中的水、 热环境。多年冻土的冻胀融沉过程不仅毁坏铁路、 公路，也会破坏建筑工程结构，甚至造成极大的生态、 地质灾害。环境变暖多年冻土热融效应增加后，会促使大团聚体破碎成小团聚体，释放大量有机碳、 硝态氮等物质，进而改变土壤生物和微生物环境，引起高寒生态环境的恶化。同时因气候变暖土壤水、 热状况发生改变，同样会使植被条件发生变化，进而影响到地表的一系列特征，如反照率、 降水的渗透速度、 土壤中的蒸腾和蒸散、 以及土壤侵蚀等，从而打乱水文和气候系统的循环速率。上述近地表过程的变化，会在大气圈的下垫面底部对环境产生反作用，甚至影响多年冻土周边地区乃至全球的气候变化，进而影响到生态系统的发生、 发展进程。

新疆地处我国西北干旱区生态脆弱带，是我国五大畜牧业生产基地之一，生态资源丰富。冰冻圈要素中的冰川、 冻土是区域内水循环过程中不可或缺的重要因子。新疆地区的多年冻土分布可称为“两环”模式，多年冻土分布于阿尔泰山、 天山和昆仑山构成的“两环”脊线上，“两环”内部的准格尔盆地和塔里木盆地为季节冻土区。其中，天山深居内陆，分布于其上的多年冻土具有明显的垂直地带性[5]，为探讨新疆多年冻土区地表形变提供了天然的研究场地。

目前，利用InSAR技术针对天山地区的研究多集中于对冰川的讨论[2]，较少应用到对多年冻土的分析中，本工作试图利用SBAS-InSAR方法和39景ENVISAT影像(覆盖时间从2003年6月17日到2010年6月15日)对天山多年冻土区地表变形规律做出部分探讨，以期为后续的冰冻圈相关科学研究提供数据基础和方法借鉴。

1 研究区概括

研究区位于天山山脉中段，位于北纬43.4°—44.7°，东经83.9°—85.6°之间(图1所示)，区域内北部主要为平原地区，南部基本为山区。天山山脉深居大陆内陆，属于典型的大陆性山地气候。根据早期研究可知[5]，天山山脉具有典型的冻土分布特征，区域内高寒草甸和沼泽草甸分布广泛。

图1 研究区概括(a): 研究区地理位置；(b): 高程图Fig.1 Location of study areas(a): Distribution of study area; (b): Digital elevation model

2003年至2010年期间，区域内年均温和年降水分别为(-2.7±8.27) ℃和(340±85) mL，如图2所示(依据《中国区域高时空分辨率地面气象要素驱动数据集》[6-7]计算获取)，其基本分布规律为北部平原地区温度高、 降水少，南部山地温度低、 降水多，以上雨、 热分布特点为多年冻土的发育提供了十分有利的条件。以3 000 m海拔为界，研究区的山地和平原地区2003年至2010年间的年均温和年降水分别为(-10.87±2.98) ℃、 (420±38) mL和(0.54±7.34) ℃，(310±77) mL，具有明显的差异。

图2 2003年—2010年间多年平均温度(℃)(a)和降水(mm·yr-1)(b)分布图Fig.2 The spatial distribution of annual average temperature (℃) (a) and annual average precipitation (mm·yr-1) (b) from 2003 to 2010

2 地表形变反演研究算法

2.1 研究数据

如表1所示，本研究中涉及的ASAR (advanced synthetic aperture radar)数据主要有39景。经初始处理[8]后，因基线和多普勒质心差的原因，其中6景影像没有参与配对。该ENVISAT ASAR数据集属于单视复数影像(single look complex, SLC)，为C波段、 IS2模式遥感数据，入射角23°。研究数据获取方式为降轨，数据来源于欧空局(European Space Agency, ESA)。方位向分辨率和距离向分辨率分别为22.6和4.05 m。本文中所用DEM资料为STRM V4版本数据，该数据在赤道的分辨率约为90 m[9]。《中国区域高时空分辨率地面气象要素驱动数据集》[6-7]下载于寒区旱区科学数据中心，覆盖时间范围为2003年1月1日到2010年12月31日，包括月、 日、 年三种时间尺度数据，其空间分辨率为0.1°。此外，文中所有地图投影方式均为WGS84坐标系。

表1 本研究中ENVISAT ASAR影像集的详细信息(连接图生成时主影像为20090106)Table 1 ENVISAT ASAR datasets products applied in this study (The image of 20090106 was used as the super master-imagine generating a connected graph)

2.2 算法原理介绍

星载InSAR技术在地表形变监测领域应用已经十分成熟，其工作原理主要是集合成孔径雷达技术与干涉测量技术于一体，计算距离信息的经典干涉相位计算公式如式(1)[3]

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

(1)

式(1)中,φflat指平地相位，φtopo指地形相位，φdef指地表形变相位；φatmo指大气延迟相位；φnoise指噪声引起的相位。首次将该技术应用于厘米级地表形变的研究方法称为雷达差分干涉技术(D-InSAR)[10-11]。在多年冻土区的研究指出早期InSAR方法的主要限制因素是时空失相关和冻土地表的非线性形变[9, 12-13]，而SBAS-InSAR方法则因奇异值分解法(singular value decomposition, SVD)可以很好地克服上述问题[14]，SBAS-InSAR方法的具体原理如式(2)所示[15]。首先，生成差分干涉图。干涉图的数量如式(2)所示

(N+1)/2≤M≤N(N+1)/2

(2)

式(2)中，N+1景的ASAR研究区影像获取于具体的研究时期，M对干涉图则需满足特定的时空基线阈值设置规则，假设每期ASAR影像都至少有1景影像与之干涉，则如式(2)，N+1景的ASAR影像生成的M对干涉图数量应在(N+1)/2至N(N+1)/2之间。时空基线的阈值设置规则一般通过多普勒质心差来控制[13]。之后，通过Ta和Tb时刻获取的ASAR影像来计算像元上的解缠相位，假设该像元在方位向和距离向的数值为(x,r)，则生成第i景的干涉图公式为

d(Ta，x，r)]+Δφi，topo(x，r)+Δφi,orb+Δφi，res(x，r)

(3)

式(3)中，φ是干涉相位信息；i是干涉图的序号，i∈(1, 2, …,M)；λ是传输信号的中心波长(本研究中ENVISAT ASAR数据的中心波长是5.67 cm)；d(Tb，x，r)和d(Ta，x，r)是雷达视线方向(line-of-sight，LOS)的积累量，Ta和Tb时刻主要相比T0时刻而言(T0时刻是参考影像的获取时间)；Δφi,topo是来自于DEM误差的相位信息；Δφi,orb是来自于卫星轨道信息误差的相位信息；Δφi,res是来自于残余相位数据的相位信息。

假设vj是p到q时刻的平均相位速度，则相位和时间满足如式(4)

vj×(Tp-Tq)=Δ(φp-φq)

(4)

式(4)中，p和q是获取ASAR影像的具体时刻，φ是相位信息。考虑到不同的时间分段，T0到Tj时刻的第j幅干涉图的相位信息可以通过如式(5)计算得到

(5)

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

Av=Δφ

(6)

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

2.3 数据处理流程

研究中的数据处理流程如图3所示。其中生成连接图时，选取空间基线距小于500 m和时间基线距小于550 d的ASAR影像，最终生成M幅(本研究中M为126)差分干涉图，图4显示结果为生成连接图时的时空基线图。

图3 地表形变反演流程Fig.3 Processing of surface deformation calculation

图4 时空基线图(a): 空间基线距；(b): 时间基线距；(c): 3D空间基线距黄色、 绿色和红色点分别代表超级主影像、 从影像及未参与配对影像Fig.4 Plots of space and time baselines(a): Plot of position baseline; (b): Plot of time baseline distance; (c): Delaunay 3D plot of position baseline distance Blue lines represent interferometric pairs；Yellow, green and red diamonds denote super-master, slave and no-matched images, respectively

在生成连接图后，配合DEM数据，经过干涉图去平、 自适应滤波、 相干系数图生成、 相位解缠后，再对M对干涉图进行编辑，保留处理较优的结果，例如图5(a，b)所示分别为20090421和20081202干涉的结果)，对处理较差的结果(例如图6(a，b)所示分别为20050308和20060221干涉的结果)进行剔除处理，此时移除的干涉图数量即为L(本研究中剔除的L干涉图为52对)。

图5 较优干涉结果(a): 较优结果；(b): 较差结果Fig.5 High quality results of interferograms(a): High quality result; (b): Low quality result

图6 较差干涉结果(a): 20050308; (b): 20060621Fig.6 Low quality results of interferograms(a): 20050308; (b): 20060621

之后利用地面控制点，进行轨道精炼和重去平处理，然后估算形变速率和残余地形，经过相关系数阈值控制和SVD方法计算，再结合空间低通滤波和时间高通滤波方法，完成最终的研究区时间序列地表形变反演，如图7所示。

3 反演结果与分析

研究区反演结果共包括33期形变影像，除第一景20040113(参考影像)的所有数值为0外，其余32期形变结果都是相对参考影像的数值[14]。本文以其中4期形变量影像为例，包括20041228(与20040113相比，季节接近)、 20050726(与20040113相比，季节相反)、 20090106(与20040113相比，季节接近)和20100615(与20040113相比，季节相反)，如图7(a,b,c,d)所示。根据图7可知，研究区内反演质量较佳区域为平原地区，山区的地表变形反演稀少而零散。不过如有需要，可通过各种空间插值方法来获取完整的面状数据，从而用以评估感兴趣区域。为保证更准确地探索地表形变与温度、 降水的关系，不做插值处理。此外，仅通过从影像表现来看，20041228和20050726在山区西部表现出地表抬升状况，而在20090106和20100615两期影像中则转变为沉降现象，四期影像在山区东部的区域基本都呈现出抬升的态势。而在平原地区，地表的沉降和抬升在4期影像中变化不明显，大部分表现为抬升，在接近城市的区域多表现为沉降。

图7 研究区4期示例形变结果(a): 2004/12/28; (b): 2005/07/26; (c): 2009/01/06; (d): 2010/06/15Fig.7 Deformation results of 4 example imagines in the study area(a): 2004/12/28; (b): 2005/07/26; (c): 2009/01/06; (d): 2010/06/15

在研究区监测范围内，2004年至2010年的地表形变速率如图8所示。研究区内平原地区的形变主要表现为抬升，仅北方接近城市的区域表现出强烈的沉降现象。研究区内山区的形变结果较为零散，大体上在西部和东部地区分别呈现出沉降和抬升现象。尽管研究区范围内有不同程度的沉降和抬升，但研究期间内的形变速率基本在±5 cm·yr-1之内，平均形变速率为(-0.07±3.38) mm·yr-1。虽然整体变形不是十分明显，但是均值体现出该区域内地表存在轻微的沉降现象，同全球气候变暖导致冻土退化的情形一致。此外，与北方西部湖泊附近沉降速率在-6和7 cm·yr-1之间相比，位于平原区北方中部的乌苏市地表沉降速率基本在-10 cm·yr-1以上。山区的整体形变基调以沉降为主，符合全球变暖冻土开始退化的规律。山区每一个有具体数值的形变点满足海拔大于3 000 m的提取条件，其具体分布位置如图9所示。利用温度和降水数据，分析这些形变点的年变化规律如图10所示。总体和不同变形速率间隔的形变样本点基本都呈现出逐步变暖的趋势，而且同样受2008年末和2009年初的极端降温事情影响。不过除地表形变速率小于-2.0 cm·yr-1的样点在2010年温度仍然有所下降外，区域形变样本点在2010年都有继续变暖的趋势。研究区的山区地表形变速率小于-2.0 cm·yr-1、 在-2.0到2.0 cm·yr-1之间和大于2.0 cm·yr-1的样本点数量依次为6 364，6 449和2 385个。

图8 平均形变速率Fig.8 Average deformation rates

图9 山区形变样点空间位置分布图Fig.9 Deformation points in mountain regions of study areas

图10 山区形变样点年均温度和降水曲线图(a): 所有样本点；(b): 小于2 cm·yr-1样本点；(c): -2～2 cm·yr-1样本点；(d): 大于2 cm·yr-1样本点Fig.10 Plots of annual precipitations and annual mean temperatures in mountain regions(a): All points；(b): Ground deformation rate <2 cm·yr-1；(c): Ground deformation rate between -2～2 cm·yr-1；(d): Ground deformation rate >2 cm·yr-1

4 结 论

近年来，因人类活动的干扰和全球气候变暖的日益严重，不同区域的多年冻土开始出现退化现象，退化严重的地区甚至会发生强烈的水土流失，对草场、 环境、 建筑及工程设置造成威胁[3]。特别是对于天山地区，不仅存在“全球气候变化指示器”的多年冻土，而且还是重要的畜牧业生产基地。无论是在环境保护，还是人类生产、 生活方面，都需要着重监测。

下列的几方面需要重点关注：

(1)地表形变反演结果从空间来分析，虽然在平原地区表现优良，但是在山区却较为稀少和零散。即使可以后期可以进行插值处理，却不一定能完全满足一些需要高精度形变数据的研究和应用。后期可以考虑利用其他软件和不同InSAR方法进行反演对比。

(2)从时间角度分析，难以反演不间断的长时间序列拟合形变结果。现有常用的InSAR数据，特别是开源数据反演结果基本上一年仅存在几期。当利用月尺度气象或其他科学数据集对比分析时(如本文中使用的温度和降水数据)，难以实现一对一匹配。如果试图利用日或者小时尺度的数据集进行分析时，则存在更大的科学挑战，因此后期地表形变的研究需要重点关注更细时间尺度的长时间序列地表形变结果反演。

(3)冻土变化的时间滞后性也需着重考虑。现有研究对时间滞后性的分析不是十分充足，但是冻土发生变化处于地表内部，自然因素对其的影响必然存在一个积累后再变化的过程。后期针对冻土变化的研究分析需要越来越多的关注时间滞后性对其影响，衡量温度或者降水等要素对其时间滞后尺度和周期的具体数值，正是一个亟需科学研究揭示的关键点。