利用OICC-SBAS技术提取岗普冰川表面流速时空变化特征

欧海沨,字城岱,滕兴发,关舒丹

(西南交通大学地球科学与环境工程学院,四川 成都 611756)

冰川受人类活动影响较小,受气候变化影响的同时,也对其作出快速的响应[1-2]。在全球气候变暖的影响下,我国大部分冰川呈现出退缩减弱趋势[3],导致由冰川变化引发的区域淹没、冰湖溃决、滑坡和碎屑流等自然灾害的规模增大、数量增加。提取并分析冰川表面运动信息,不仅能揭示全球气候变暖背景下冰川的动态变化信息,同时也能为冰川变化诱导的灾害预警提供重要依据[4-6]。由于冰川多分布在海拔较高、地势严峻、条件恶劣的地区,传统的监测技术受到诸多限制。而空间遥感技术体现出更多的优势,包括覆盖范围大、数据源广、监测周期长、无须现场勘测等。目前,在冰川表面流速的监测中,普遍采用的遥感手段包括3种:合成孔径雷达(SAR)影像的差分干涉[7-8]、SAR像素偏移追踪技术(POT)[9-10]及光学影像互相关技术(OICC)[11-12]。

其中,光学影像具有更加丰富的纹理信息,更有利于监测地表的几何变化,同时,接近垂直的入射角也可避免陡坡之间的相互遮挡。虽然光学影像质量易受云层的限制,但随着高性能光学遥感卫星数量的不断增加,光学影像的存档数量也在不断积累,为及时、准确获取冰川动态变化信息提供了可能。本文以加拉白垒峰北侧的岗普冰川为研究对象,位于雅鲁藏布江大拐弯地带,该地区冰川多属于海洋型冰川,运动速度较快;借鉴短基线集(SBAS)思想,利用OICC-SBAS技术提取冰川表面运动状态,并分析其时空变化特征。

1 研究区与数据

1.1 研究区概况

加拉白垒峰位于念唐古拉山脉,属于高原温带半湿润季风气候,具有雨量丰沛、干湿季节分明、昼夜温差大等特点,受降雨和温度的影响,该区域内冰崩和雪崩等灾害十分频繁。研究对象岗普冰川为加拉白垒峰周边规模较大的冰川之一,地势南高北低,东侧发育着一条大型的支流(如图1所示),区域内最高海拔为4169 m,最低海拔为3438 m。

图1 研究区地理位置及冰川地貌

1.2 研究数据

Sentinel-2包含A、B两颗卫星,已收集了约7年的存档数据,且具有重访周期短、中等分辨率和免费开放等特点。Sentinel-2传感器共提供13个光谱波段,空间分辨率位于10～60 m之间,本文采用空间分辨率最高的近红外波段。从欧洲航天局(European Space Agency,ESA)网站获取了2015年12月6日至2022年3月9日的12景Sentinel-2遥感影像,详细参数信息见表1。

表1 Sentinel-2影像参数

2 研究方法与数据处理

为获取岗普冰川的时序表面运动速度及时空变化特征,本文采用OICC-SBAS时间序列反演算法,如图2所示。该算法共包含5个模块:影像配对策略、数据预处理、COSI-Corr批处理、误差后处理及时序反演。

图2 数据处理流程

2.1 影像配对策略与预处理

不同影像之间存在着太阳光照条件、轨道位置、采集时间等因素的差异,这些因素与图像相关性测量值之间存在着一定的关系,被称为基线变量[13],包括太阳方位角差值(SAD)、太阳高度角差值(SED)、时间基线(TB)和空间基线(SB)。通过设置基线变量的阈值,可优化影像对的选择,具体见表2。在后续处理中手动去除质量较差的影像对,最终获取14组影像对。由于获取的Sentinel-2数据为L1C级数据,即经过了粗略几何校正和正射校正[13-14],可直接用于影像互相关处理。预处理只需提取覆盖研究区域的范围,即裁剪出包含相同地理覆盖的子图像,并按时间顺序排列,形成图像堆栈。

表2 基线阈值与批处理相关参数

2.2 偏移量估计与误差后处理

以COSI-Corr软件为数据处理平台,经多次试验,具体参数设置为:初始搜索窗口和最终搜索窗口均设置为64×64,步长设置为1,其他参数(稳健性迭代次数和掩膜阈值)均采用系统默认值,确保在降低测量不确定度的同时具有较高的计算效率。

初始相关结果中仍然存在大量的噪声和伪影,包括失相关噪声误差、轨道误差、条带误差、卫星姿态角误差,需要进一步后处理。针对异常值和失相关噪声,将信噪比值低于0.96的像素和超过TB阈值估计的最大预期位移的像素进行掩膜。对于轨道误差,采用一次多项式仿射变换模型解决;对于条带误差,通过均值相减法去除;对于卫星姿态角误差,采用改进的均值相减法[15]去除;最后再对相关图进行非局部均值滤波,获得进行了误差处理的东西向形变图和南北向形变图。

2.3 时序偏移量计算

基于互相关技术的位移时间序列反演借鉴了InSAR时序反演中的小基线时序分析技术(SBAS)[15]。假设有N+1幅光学影像,获取时间序列为[t0,t1,…,tn],M为互相关影像对。以影像对在每个连续时间间隔上的平均速率v作为未知参数,构建观测方程为

B·v=δd

(1)

式中,δd为影像对的位移序列;B为由相邻影像的时间间隔组成系数矩阵。当系数矩阵B满秩时,可采用最小二乘求解代求参数v,即

v=(BTB)-1BTδd

(2)

但当系数矩阵B秩亏时,需要采用SVD分解法求解参数v。为了更好地考虑去相关随时间的变化,根据主影像和从影像之间的时间间隔,加入权重矩阵P,P=diag(w1,w2,…,wM),其中

(3)

Rl为原始测量值和反演位移之间差值,即Rl=δd-Bv。高残差代表高偏差或低质量的相关值,这些值会直接影响反演的精度。为了避免该问题,利用残差Rl对反演进行第2次迭代,即

(4)

式中,R0为一个常数,近似于“成功”相关性的噪声标准差。因此,权重Wil限制了非常高的残差异常值对反演结果的影响,但不会人为地高估残差接近于零的相关值。获取反演结果后,可以通过掩膜和滤波对结果进一步改进。首先,根据文献[16]中的方法,过滤方向标准差较大的像素值,并去除运动方向与实际斜坡方向大于90°的值;然后,掩膜后对导出的位移时间序列进行方向滤波[17]。

3 结果与分析

3.1 冰川运动结果可靠性分析

为验证SAD与SED对其相关性结果的影响程度,在研究区域内选取3个稳定区域(如图1(b)中的Roi1、Roi2、Roi3),对3个区域东西向和南北向的形变均值进行精度评定,最终获得的相互关系图如图3所示。通过统计发现,当SAD与SED较大时,其均值较高;当SAD与SED较小时,其均值较低。这充分说明在适当的时间基线下,选择较小的SAD与SED,能够有效提高结果精度。

图3 太阳高度角差值和太阳方位角差值与稳定区域形变均值的关系

由于缺乏冰川的实测数据,在精度评定时,通常采用稳定区域的均值和标准差等先验信息进行辅助。理论上稳定区域的速率为零,因此提取3个稳定区域的残余运动,并绘制相应的残余运动频率分布图(如图4所示)。其中,Roi1、Roi2靠近冰川,地形起伏较大,且在一定程度上受积雪影响,残余运动偏大;Roi3离冰川较远,残余运动明显减小。误差统计结果见表3。结果显示,误差值相对于冰川表面最大流速范围为0.04%～0.7%,稳定区域平均运动速度远小于冰川实际运动速度,保证了冰川时序流速结果的可靠性与准确性。

表3 稳定区域运动速度统计结果(m/a)

图4 稳定区域残余运动频率分布

3.2 冰川表面流速空间分布特征

本文通过2015—2022年获取的12幅Sentinel-2光学影像,基于OICC-SBAS反演算法提取了岗普冰川的年平均表面流速(如图5所示)。由于两侧山体挤压,冰川中间的表面流速大于两侧,符合冰川运动规律。同时,利用东西向和南北向偏移量获取的偏移方向,与冰川实际运动方向相符(如图6(a)所示)。冰川主流以南北向形变为主,支流以东西向形变为主,这主要是由坡度、坡向等地形因子联合导致的。为了进一步分析冰川表面流速的空间分布特征,在岗普冰川上提取了3条剖面线(剖面线位置如图5(b)中的AA′、BB′、CC′),并联合高程、坡度信息进行分析,结果如图6所示。通过提取冰川剖线与地形因子的关系,得到如下结论。

图5 冰川年均表面流速

图6 岗普冰川剖线表面流速与地形因子关系

(1)剖面线AA′:剖面线AA′为冰川主流中心线,以南北向形变为主。沿该剖面线的高程和坡度变化较平缓,流速波动较小,但也可以看出3个波动现象,分别位于剖面线的0.5、3.2、5.6 km处(分别为图6(a)中的T1、T2、T3)。距T1处的加速主要是由于上游支流在此汇集,从而导致冰川加速;而在T2处,虽有东侧支流的汇入,但冰川流速减弱,主要是由于该支流在汇入主冰川之前其流速已经减小至1.24 m/a,因此对主冰川流速未表现出较大的推动作用;T3位于冰舌末端,由于冰川物质积累和地形变化,也有一个加速现象。

(2)剖面线BB′:剖面线BB′位于冰川的东侧支流,以东西向流速为主。沿剖面线BB′的地形变化较明显,该支流流速随高程的减小而逐渐减小;在剖线末端,由于汇入主流,导致南北向流速有一个明显的增速现象。

(3)剖面线CC′:该剖面线横跨上游东西两条支流,地形为中间低、两头高。该支流流速在1.0、3.9 km处(分别为图6(c)中的T1、T2)随着冰碛物的积累和地形变换,其东西向流速有一段快速加速过程;后因逐渐汇入主流,开始以南北向流速为主导。冰川西侧的上游支流为表面流速高值的主要聚集区域,该区域为冰川的物质积累区域,形变量级较大且形变特征复杂,年最高流速可达82.5 m/a。

综上所述,岗普冰川主体流速受支流影响,但整体无剧烈变化,维持在40～50 m/a之间。东侧支流随着高程减小,其冰川流速也逐渐减小,位于20～40 m/a之间。岗普冰川上游为整条冰川流速高值的主要聚集区域,形变量大且形变特征较为复杂。

3.3 冰川表面流速时序变化特征

通过OICC-SBAS技术处理得到的2015—2022年的冰川时序形变结果如图7所示。可以明显看出,2017年2月—2018年3月、2020年11月—2021年12月,这两时段内的累计形变明显大于其他时期。一方面是由于这两时段的时间间隔较长,另一方面则是受气候因素的影响。

图7 时序形变处理结果

为了研究岗普冰川的变化规律,通过美国国家海洋和大气管理局(NOAA)网站获取距离研究区域最近的林芝站点降水量的逐日统计结果。在岗普冰川上选取3个区域(分别为图7中的A、B、C),分析冰川形变与气候因素的关系(如图8所示)。可以看出,研究时段内年降水量具有明显差异,2017年2月—2018年3月、2020年11月—2021年12月均具有较高的累计降水量。其中,2016年3月—2017年2月具有最高的累计降水量。该时段时间间隔较短,导致累计形变不是最大,但该时段形变仍然显著高于其他累计降水量较小的时段。总体而言,冰川的形变量与降水变化较为一致,在累计降水量高且时间间隔较长的时段内,冰川累计形变较高,而在累计降水量较低的时段内,冰川的形变量有明显减少。

图8 岗普冰川累计形变与气温降水关系

4 结 语

本文利用2015—2022年获取的12景Sentinel-2影像,采用OICC-SBAS反演算法,提取岗普冰川表面时序运动速度分布信息,并进行定性、定量分析。空间分布特征上,岗普冰川位移受地形和气候因子的共同作用,流速高值主要集中在上游的物质积累区,该区形变特征较为复杂,最高流速超过80 m/a。冰川主体运动速度主要为40～50 m/a,受地形和支流汇集影响,表现出轻微波动,但整体较为稳定。东侧支流作为岗普冰川水资源的重要补给,其运动速度随着高程减小而缓慢减小,主要维持在20～40 m/a。时序变化特征上,研究时段内冰川流速无较大变化,年际之间虽然有轻微的差距,但整体呈现出较高的相似性。即2015—2022年,岗普冰川都维持在一个稳定流动的状态。且岗普冰川的流速变化与降水变化较为一致,能较好地反映出气候变化的特征。

同时,本文利用东西向和南北向偏移量获取的最终偏移方向与冰川实际流动方向相符,偏移量误差仅占冰川最大流速的0.04%～0.7%,稳定区域的残余运动远小于冰川的实际运动速度,有力地证实了冰川时序表面流速结果具有较高的可靠性和检测精度。