2. 长安大学公路学院, 陕西 西安 710064
2. College of Highway Engineering, Chang'an University, Xi'an 710064, China
2008年5月12日四川省汶川县(31.0°N,103.4°E)MS8.0地震,是近年来破坏程度和强烈程度最为严重的一次[1]。地壳运动的实质是地球内部构造应力作用下引起的地壳一些相对运动(包括垂直运动和水平运动等)。从预变应变建立到同震移位,进一步到后震弹性恢复,地震周围地表变速变形,甚至达几十厘米[2]。本文利用小基线技术对汶川地区震前进行变形监测研究,得到一年之内汶川地区地面的形变速率,并将研究结果与前人结果进行对比,探讨地震前地表形变呈现的规律。
汶川地震震区海拔高,地势险峻,利用水准测量和布设GPS形变观测点的方法较困难,很难获取研究区域内长时间的形变情况,即使获取到部分图像,也不是连续的[3]。曾有人通过分析获取到的断层和GPS观测资料,得到了与经验和理解相悖的结果,因为在地震爆发前期,无明显的变化结果,且越靠近震中,变化越不明显[4]。显然,这是与现实不符的。
D-InSAR是近几年发展起来的一种用于地形和地面变形测绘的强大技术,有着全天候和昼夜成像的能力以及宽空间覆盖高分辨率和高测量精度的特点[5]。2004年, 李震等采用该技术对青藏公路沿线进行形变监测,所得结果与实地监测控制点形变数据非常吻合[6];但是,该技术因为监测精度受时间和空间相关性以及大气延迟的影响,只能得到单个变形成果,不能得到研究区域表面变形的时间演化情况[7]。
基于此,部分国际学者提出了基于D-InSAR的时序InSAR技术(multi-temporal InSAR),通过对大量数据在高相位点目标的分析,获得了研究区地表面变形速率和变形时间序列[8]。短基线集(SBAS-InSAR)技术通过多干涉对组合对多个解缠相位图进行最小二乘求解,可以有效消除或削弱解缠粗差和大气延迟误差的影响,高精度地获取地面形变量[9]。此外,该方法无需大量SAR影像,是一种非常有效的InSAR处理方法,在大城市地面沉降、火山缓慢变化、断层蠕动等方面表现出了极大的潜力和优势[10]。本文尝试利用该技术提取汶川震前震中位置的形变时间序列,以探究发震前期活动的状态。
1 SBAS-InSAR技术 1.1 SBAS-InSAR技术原理SBAS-InSAR是由Berardino和Lanari等研究人员提出的一种不同于PS-InSAR的时间序列分析方法,可以较好地克服空间去相干[11]。SBAS-InSAR采用了奇异值分解(SVD)法求解形变速率,使较大空间基线分开成孤立SAR数据集连接起来,大大提高了观测采样率[12]。
图 1为SBAS-InSAR数据处理流程。
SBAS-InSAR方法是基于时间序列SAR数据生成的小基线差分干涉图集,它可以利用平时通用的InSAR处理软件获得的差分干涉图来得到地表相干目标的形变速率和时间序列[13]。小基线技术的关键特征在于充分利用获取的SAR数据,因此提高了形变量测的时间采样频率和研究区域的空间覆盖率。SBAS-InSAR应用于影像覆盖范围内所有表现出足够高相干性的像元,因此算法是很稳健的[14]。
1.2 Stacking原理Stacking方法数学模型为
式中,ph_rate为年沉降相位速率;wi为加权因子,wi=Δt-1,Δt为干涉对时间基线;phi为干涉对解缠相位值。
Stacking方法实质是:设定大气影响为时间上的随机信号,而形变则为线性变化,当
本次试验选用由美国宇航局NASA测得的SRTM数字高程模型(如图 2所示)作为地形数据。SRTM系统覆盖面积之广、采集数据量之大、精度之高在测绘史上是前所未有的。该模型的地面分辨率为3 m,垂直和水平精度分别为16和20 m,选用WGS-84坐标系为平面基准,高程基准是EGM96。
本文选择23.6 cm的L波段数据,较常规5.6 cm的C波段数据信号更好,获得的数据信噪比更高。由于卫星上直接获得的raw数据不能够在Gamma终端下运行,必须将通过互联网获得的卫星精密星历,转换成软件可以处理的单视复影像图(SLC Image)。
为达到最佳配准精度,本次试验选取时间居中的一景影像为主影像,其余7幅影像与其进行配准,查看配准结果发现,其中两幅数据由于大气和噪声影响,配准效果较差,达不到研究要求,将这两幅剔除。故本次试验能够完全配准的只有6幅影像。表 1是所用影像的基本参数。
编号 | 摄像日期 | Track | Sensor | Incident angle | Mode | 垂直基线距离/m |
1 | 2007-06-20 | 475 | PALSAR | 38.7373 | FBD | -165.809 6 |
2 | 2007-08-05 | 475 | PALSAR | 38.7370 | FBD | 217.873 7 |
3 | 2007-09-20 | 475 | PALSAR | 38.7272 | FBD | 0 |
4 | 2007-12-21 | 475 | PALSAR | 38.7396 | FBS | 447.638 7 |
5 | 2008-02-05 | 475 | PALSAR | 38.7323 | FBS | 1 157.938 1 |
6 | 2008-05-07 | 475 | PALSAR | 38.7330 | FBD | 1 702.306 6 |
选中20070920为主影像,6幅影像两两干涉生成12幅干涉图,但由于其中3幅不满足研究要求,故而只有9幅。接着对干涉图多视处理来增强图像相干性和检测较大尺度的形变信息。多视因子设为12:28。采用两次自适应滤波法,去除由于失相关引起的微小形变,从而达到提高信噪比的目的[14]。两次滤波窗口分别设为256、128,滤波因子设为0.5。数据处理过程中,许多失相干因素的影响会间断地获取干涉相位,使得相位周期性和趋势性受到破坏,增加相位解缠的难度,更甚的是可能造成相位解缠出错。而小基线多干涉对组合可剔除相位解缠误差。截至目前,相位解缠方法多达30多种。本次经过几种方法对比,最后所用的是效果较好的最小费用流法,解缠时相干系数阈值设为0.7。
所选地区由部分山区和扬子地块组成,而山区部分由于植被覆盖率较高,大气(水汽等)影响最后结果的真实性,故而为了使最后生成的时间序列和年平均形变速率更加精确,对影像进行了一次去大气操作和两次去趋势操作。最后,用stacking技术获取年平均形变速率图。
2.3 InSAR监测震前变形剖面分析龙门山断裂带由龙门山后山断裂带(汶川-茂县-平武-青川)、主中央断裂带(映秀-北川-关庄)、龙门山主山前边界大断裂(都江堰-汉旺-安县)3条主断裂带组成[15]。为了进一步对该地区的震前变形活动进行空间域的分析,根据断裂带走向和分布,提取了3条断裂中变形较为明显的汶川-茂汶断裂带以及主中央断裂带,作为南北剖面线(如图 3 H1H2所示)。
H1、H2为汶川-茂汶段裂带的剖面线,L1、L2为主中央断裂带剖面线,图 4(a)为北侧H1H2的剖面线形变图,最大沉降值达到14 cm,最大隆升值达19 cm,相比于南侧的L1L2(如图 4(b)所示),变化较为凸出,而L1L2整体变化较为平缓,活动区域基本在10 cm左右。图 4(c)为两条剖面的叠加图,可以看出,汶茂断裂带比映秀-北川断裂带更为活跃,这也解释了汶茂断裂带为主要断裂带,发生地震的严重程度要远大于都江堰一带。
V1-V6为垂直于H1H2的6条剖面线,其形变如图 5所示,从中可以看到V1剖面线地形沉降幅度最大,大约在24 cm内左右,并且在V1、V4、V6处出现了斑状沉降区,而沿着汶茂断裂带出现了“凸”状隆升区,最大隆升值甚至达到22 cm。在绵虒镇附近沿着汶川一带有大范围的隆起,最大幅度可达18 cm,这也可能预示着此区域附近为地震初始破裂位置,与单新建等人在2015年的调查结果一致。汶川地震震前形变场的分布格局显示,沿着断层线呈现出“凸”状形变变化特征, 即以发震断层为中心的隆起形变变化,远离发震断层两盘均出现下沉。
2.4 时间序列分析针对同一剖面V3做时间序列,如图 5所示,得到6幅不同时间的形变图。对每个时间段的V3剖面做形变图,得到图 6。
从图 5可以看出,远离断裂带处的地区有斑状下沉区,而离断层越近,上升特征越为明显。在图 6(c)中,即2007年9月20日,达到了沉降最大值,大约24 cm。至V3西北方向剖面起点15 km处开始,出现了上升的趋势。在大约剖面线和断裂带的交点,出现了急剧的上升。2008年2月,上升值达到了最大,大约在20 cm左右。随后在断裂带的东南侧,开始出现下降的趋势。从图 6可以看出,是在距离西北端20 km处开始出现下降的特征。
从图 6(b)和(c)中可以看出,2007年8月和9月断裂带的东南处也有大面积的隆起。而在2007年12月以后一直到地震前夕,变形则缓和了很多,推测与断裂带起伏运动有关系。从整体上看,同比雅安地震,汶川震中在地震爆发前一年变形较大,这也与后面汶川大地震的大爆发力以及大破坏力遥相呼应。
3 结论本文总结运用了InSAR的小基线技术,通过对汶川震中地区时间间隔为320 d的数据分析,提取了汶川一个条带在震前一年的形变场,得到的震前垂直形变集中在为10~15 cm之内,个别严重变形区域可达20 cm。汶川-茂县断裂带为主要隆升区域,绵虒镇的西北端和汶川的东南段为主要沉降区域。较前人的研究结果变形略大,可能与选取数据包含了2/3的山区部分有关,山区植被多、水汽大,对变形结果有较大的干预。此外,对于其他地区震前活动(如北川、安县、绵阳区域),在InSAR形变图上则不太明显,推测与地震孕震规律有关。还需采用可以覆盖龙门断裂带区较多的数据,进行长时间序列监测。而对于地面垂直变形原因与地震活动的具体关系,还需结合构造地质学与地震活动资料等进一步加以佐证。
[1] | 江在森, 方颖, 武艳强, 等. 汶川8.0级地震前区域地壳运动与变形动态过程[J]. 地球物理学报, 2009, 52(2): 505–518. |
[2] | SHAN X J, ZHANG G H. A Characteristic Analysis of the Dynamic Evolution of Preseismic-coseismic-postseismic Interferometric Deformation Fields Associated with the MS7.9 Earthquake of Mani, Tibet in 1997[J]. Acta Geologica Sinica, 2007, 81(4): 587–592. DOI:10.1111/acgs.2007.81.issue-4 |
[3] | 单新建, 宋小刚, 韩宇飞, 等. 汶川MS8.0地震前In-SAR垂直形变场变化特征研究[J]. 地球物理学报, 2009, 52(11): 2739–2745. |
[4] | 李媛, 牛安福, 刘希康, 等. 汶川MS8.0地震震前近震源区地壳形变机制探究[J]. 地震学报, 2015, 37(6): 959–972. DOI:10.11939/jass.2015.06.007 |
[5] | 许才军, 林敦灵, 温扬茂. 利用InSAR数据的汶川地震形变场提取及分析[J]. 武汉大学学报(信息科学版), 2010, 35(10): 1138–1142. |
[6] | 李珊珊, 李志伟, 胡俊, 等. SBAS-InSAR技术监测青藏高原季节性冻土形变[J]. 地球物理学报, 2013, 56(5): 1476–1485. DOI:10.6038/cjg20130506 |
[7] | 梁涛. 利用短基线集InSAR技术监测矿区地表形变[J]. 测绘通报, 2014(S1): 82–84. |
[8] | 季灵运, 王庆良, 崔笃信, 等. 利用SBAS-DInSAR技术提取腾冲火山区形变时间序列[J]. 大地测量与地球动力学, 2011, 31(4): 149–153. |
[9] | 喻永平, 林鸿, 王会强. 利用时序InSAR技术监测广花盆地地面沉降[J]. 测绘通报, 2015(S1): 157–159. |
[10] | OSMANOGLU B, SUNAR F, WDOWINSKI S, et al. Time Series Analysis of InSAR Data:Methods and Trends[J]. ISPRS Journal of Photogrammetry and Remote Sensing, 2016(115): 90–102. |
[11] | 严剑锋, 杨俊凯, 赵伟颖, 等. 基于SBAS的采动区高速路沉降监测与分析[J]. 测绘通报, 2016(1): 148–149. |
[12] | 尹宏杰, 朱建军, 李志伟, 等. 基于SBAS的矿区形变监测研究[J]. 测绘学报, 2011, 40(1): 52–57. |
[13] | 程滔, 葛春青, 陶舒, 等. 小基线集合成孔径雷达干涉测量算法及其应用[J]. 测绘科学, 2015, 40(11): 96–99. |
[14] | 付政庆, 刘国林, 陶秋香, 等. 雷达干涉测量图像的梯度自适应光滑子函数滤波法[J]. 测绘学报, 2014, 43(3): 263–266. |
[15] | 李勇, 黄润秋, 周荣军, 等. 龙门山地震带的地质背景与汶川地震的地表破裂[J]. 工程地质学报, 2009, 17(1): 4–17. |