联合GPS和InSAR研究海原-六盘山断裂现今的地壳变形特征
蒋锋云, 季灵运, 朱良玉, 刘传金
摘要：
海原-六盘山构造区为青藏高原东北部构造变形最为显著的区域之一, 历史强震活动频繁, 是研究青藏高原NE向扩展的重要窗口和地震孕育过程的理想场所。文中处理了跨海原-六盘山断裂2014-2020年期间2个轨道的时序Sentinel-1A/B SAR数据, 获得了该区域InSAR视线向现今的地壳形变场。融合公开发表的近十多年时间尺度的水平GPS地壳运动速度场, 获得了研究区高密度地壳水平形变场。对比GPS、水准和InSAR观测结果, 以及GPS-InSAR融合的高密度水平形变场, 分析讨论了该区域的地壳形变、应变场特征及其与构造之间的对应关系。主要结论如下: 1)GPS和InSAR观测表明, 1920年海原8.5级大地震的震后黏弹性松弛效应在海原断裂南侧至今仍较为明显; 2)GPS-InSAR高分辨率水平形变场表明, 狭义海原断裂左旋滑动速率的递减主要发生在中东段, 而中西段递减并不显著, 可能与海原断裂向六盘山断裂之间由左旋走滑向逆冲推覆构造转换有关; 3)六盘山断裂中-南段的地壳垂直形变和水平形变场特征均显示, 该段断裂可能处于强震孕育的中晚期, 根据反演得到的断层运动参数和地震地质资料, 估算六盘山断裂中-南段发生强震的最大矩震级约达7.5级; 4)研究区应变积累较快的区域主要集中在海原断裂附近和海原断裂-香山-天景山断裂之间的左旋剪切区, 香山-天景山断裂东南段的应变率场和周围相比明显偏小, 存在应变不匹配现象, 可能与强震孕育有关。
0 引言
海原-六盘山断裂位于中国青藏高原东北部, 是该区域主要的地貌边界与构造边界之一。该断裂带自西向东由近EW走向的冷龙岭断裂、金强河断裂、毛毛山-老虎山断裂、狭义海原断裂和近SN走向的六盘山断裂组成(张培震等, 2003; 郑文俊等, 2019)。新生代以来, 该断裂的构造活动极其强烈, 地壳变形极为显著。从历史地震活动来看, 海原断裂及周边地区历史强震非常活跃, 自1900年以来, 曾发生1920年海原8.5级和1927年古浪8.0级2次8级以上地震及多次7级以上地震(Gaudemer et al., 1995)(图1)。从历史构造演化来看, 该断裂的形成和演化主要受新生代以来印度板块与欧亚板块的会聚以及之后青藏高原的快速隆升和向外扩展所控制。该区域现今主要表现为以左旋走滑为主的海原构造系统向以逆冲推覆为主的六盘山构造系统转换的变形特征(王伟涛等, 2014)。
HYF 海原断裂(S1、S2、S3、S4); LPSF 六盘山断裂; MMSF 毛毛山断裂; LHSF 老虎山断裂; GLF 古浪断裂; XS-TJSF 香山-天景山断裂; GG-BJF 固关-宝鸡断裂
随着空间对地观测技术的不断发展, 以GPS和InSAR为代表的空间对地观测结果在地壳形变研究中得到了广泛的应用, 在构造变形、地震孕育过程等方面取得了重要成果。近年来, 将InSAR和GPS 2种观测资料结合起来, 发挥各自的优势, 获取高精度三维地壳形变场的方法日益成熟和完善(Gudmundsson et al., 2002; Samsonov et al., 2007; Hergert et al., 2010; Guglielmino et al., 2011; Hu et al., 2012; Walters et al., 2014; Fuhrmann et al., 2015; Luo et al., 2016; Song et al., 2019; Wang et al., 2019; Shen et al., 2020; Zhang et al., 2022)。其中, Song等(2019)采用2003-2010年期间覆盖青藏高原东北部广义海原断裂带的6个Envisat ASAR降轨干涉合成孔径雷达(InSAR)数据, 融合GPS地壳水平运动速度场, 获得了该区域高精度、高分辨率的地壳水平运动速度场, 并讨论和分析了水平形变场、应变场特征与构造之间的关系。但未涉及与海原构造紧密相连的六盘山地区, 也没有涉及与之相关的垂向变形方面的内容。而水准观测数据(Hao et al., 2014)表明, 六盘山地区现今的地壳垂直形变具有和汶川地震之前龙门山以西地壳的垂直形变(王庆良等, 2010; 蒋锋云等, 2011)类似的隆升特征(图2a)。
图2 海源-六盘山地区的GPS速度场
a 六盘山地区现今的水准观测垂直形变速率(据Hao等(2014)的数据绘制); b 平行于海原断裂中西段的GPS速度场。图中的颜色表示速率大小, 箭头方向表示速度的方向
这是否表明六盘山断裂的强震危险性较高? 鉴于水准和GPS垂向观测的空间分辨率较低, InSAR观测到的高密度形变场是否也存在这种特征, 值得进一步研究。另外, 我们从Wang等(2020)给出的中国大陆GPS形变场中可以看出, 在1920年海原大地震发震断层的南侧存在一个明显的地壳快速运动区域, 将GPS速度分解到平行于海原断裂带的方向之后(图2b), 这一特征更为明显, 初步推测该特征可能是1920年海原大地震震后地壳介质黏弹性松弛效应引起的变形。但Song等(2019)给出的跨海原断裂带的InSAR形变场并未显示出这一特征, 其结合GPS观测融合后的高密度地壳水平形变场也未显示这一形变特征。在地壳形变监测方面, 与Song等(2019)使用的Envisat卫星相比, Sentinel-1卫星具有更大的优势, 后者重访周期短、获取的数据量更大、观测值增多, 使得监测结果的点位密度和精度都有所提高。此外, Sentinel-1卫星采用先进的控轨技术, 使干涉图的基线更小, 更有利于开展形变监测。
为此, 本研究将处理覆盖海原-六盘山断裂的2014-2020年期间2个轨道的Sentinel-1A/B SAR数据, 获得该构造区的InSAR视线向变形场, 并结合该区域的GPS和水准观测结果讨论上述2个问题。在此基础上, 联合公开发表的青藏高原及其周缘通过GPS观测得到的近十多年尺度的地壳水平运动速度场(Wang et al., 2020), 采用Shen等(2020)给出的GPS和InSAR数据融合方法, 获取海原-六盘山地区现今高精度的地壳形变场及应变率场, 并进一步探讨海原-六盘山断裂现今的滑动速率、应变特征的空间分布与构造之间的对应关系。
1 InSAR数据及其处理过程
我们调查并收集了覆盖研究区域的2个轨道(D62和D157)2014-2020年的Sentinel-1A/B SAR数据, 主要覆盖了狭义海原断裂和六盘山断裂。对于金强河、毛毛山和老虎山断裂, 已有的研究较多(Jolivet et al., 2012; 乔鑫等, 2019; 吴东霖等, 2020; Li et al., 2021), 本文将不再涉及。对每个轨道的Frame进行上下拼接, 共获得208景影像, 数据信息见表1。
采用GAMMA商业软件平台处理InSAR数据。为获取可靠的形变场结果, 我们主要解决了以下关键问题:
(1)采用增强频谱分集方法(Yagüe-Martínez et al., 2016)并引入美国NASA发布的30m空间分辨率的SRTM DEM作为外部地形信息, 以保证配准精度(优于千分之一像元)。
(2)在选取干涉图时需要同时限定空间基线和时间基线, 本文限定垂直基线<100m, 时间间隔≥300d(图3)。
(3)为了提高干涉图在植被覆盖区域的相干性, 对干涉图进行2次多视处理及Goldstein滤波(Goldstein et al., 1998)。第1次采用的多视比为8︰2, 目的是有效利用30m分辨率SRTMDEM数据中的地形信息以消除地形相位; 第2次采用的多视比在此基础上又扩大了15倍(即120︰30), 最终得到的像元大小约为450m, 其目的是提高干涉图的相干性。
(4)采用基于Delaunay三角网的最小费用网络流(Minimum Cost Flow)算法进行相位解缠, 解缠的起始点选取在断裂附近。
(5)对于大气相位延迟引起的误差, 首先使用英国纽卡斯尔大学李振洪教授团队开发的GACOS(Generic Atmospheric Correction Online Service for InSAR)在线气象模型估计大气延迟并从干涉图中去除。对于残留的大气湍流层延迟, 根据其空间相关、时间不相关的时空域特性, 采用时间域高通滤波和空间域低通滤波进行估计并削弱(Yu et al., 2018)。
(6)对于轨道误差及残余的长波长大气误差, 采用“移去恢复法”削弱长波长非构造信号。移去恢复法的主要思想是使用外部的形变数据通过插值投影到InSAR的LOS方向, 从而得到初始形变模型。在估计InSAR轨道误差之前, 移除初始形变模型相位, 通过二次多项式拟合估计残余相位的轨道误差, 最后将初始形变模型加回到去除轨道误差的解缠图中。通过这种方式可以有效地分离轨道误差和构造形变。这里的外部数据主要指GPS和水准数据观测到的形变结果。
最后利用干涉图堆叠(Stacking)技术处理获取区域平均地壳形变速率场图像。为验证InSAR观测的可靠性, 我们将其与研究区观测时间尺度较长的GPS稳定地壳运动变形数据进行了对比分析。
2 利用GPS和InSAR获得高分辨率地壳形变场
获取精确的高密度地壳水平形变场对认识区域地壳形变细节具有重要意义。然而, InSAR观测只能获得一维的LOS向形变。因此, 仅仅使用升轨和降轨InSAR观测数据不可能获得真实可靠的三维形变场。虽然GPS方法可以以毫米级的精度测量离散位置的三维位移, 但观测站点分布稀疏且观测周期较长, 无法捕捉到三维形变场的细节。很明显, 这2种方法是互补的, 因此我们将两者结合以进行高空间分辨率的三维形变监测。
这里, 我们基于Shen等(2020)开发的一种算法, 采用D062升轨和D157降轨的InSAR数据, 以及公开发表的约十多年观测获得的GPS数据(Wang et al., 2020), 获取研究区高精度三维形变场。首先对InSAR速率图进行多视处理, 将数据压缩到1km的距离, 以提高计算效率。此外, 考虑到GPS网络的密度和数据加权配置, 基于之前提出的算法(Shen et al., 2015)对基于点的离散GPS速度采用高斯距离加权函数和Voronoi单元空间加权函数进行插值, 生成与多视InSAR数据相同网格间隔的连续三维矢量图。我们将预先确定的公共加权阈值参数(W)设置为5, 将不确定性平滑常数(σ0)设置为120km, 由此实现对解决方案平滑程度的总体控制。我们根据海原-六盘山断裂周围的GPS网络站点空间密度, 通过试错确定GPS水平速度数据输入的不确定阈值为1mm/a。
由于不同轨道相对测量选择的参考点不同, 轨道之间的InSAR视线向速度通常存在偏差, 采用GPS速度及其内插值和InSAR视线向速率, 通过最小二乘回归求解各网格单元的三维速度。
3 结果
3.1 InSAR视线向震间速率特征及其可靠性分析
我们的结果自西向东、从南到北覆盖了狭义海原断裂带和六盘山断裂带的整个区域。根据已有的地质构造、大地测量知识可知, 狭义海原断裂带作为左旋走滑构造带, 现今主要表现为近EW向水平变形, 垂向变形并不显著。而六盘山断裂带作为逆冲推覆构造带, 其现今地壳主要以垂直断裂方向的水平向会聚和垂向变形为主, 走滑运动并不显著。鉴于狭义海原断裂带和六盘山断裂带2种不同的构造变形模式, 我们分别对海原断裂和六盘山断裂的震间InSAR视线向速率及其可靠性进行简单分析。
3.1.1 狭义海原断裂带
图4a为处理得到的狭义海原断裂带的InSAR观测雷达视线向LOS方向速率。结合卫星几何图可以看出, 该区域在InSAR观测时间段表现为显著的左旋走滑构造运动, 和已有的认识一致。其他断裂两侧的差异运动相比狭义海原断裂要小得多。值得关注的是, 在狭义海原断裂带中西段西南侧的景泰、白银、靖远、海原之间有一个明显的速率相对较大的区域, 处于该区域的几个GPS站点的运动速率相对周围其他站点也偏大。
由于狭义海原断裂带主要以水平左旋走滑为主, 且垂向变形(Hao et al., 2014)相比水平向变形要小得多, 因此忽略其垂向运动分量, 通过雷达入射角(θ)、雷达视线方位角(α)和断层走向(β), 可将雷达视线方向形变速率
转换为与断层(狭义海原断裂中西段)平行方向的速率(Dfault)(Lyons et al., 2003), 计算时考虑雷达入射角随经度的变化。
从结果来看(图4b), 相比InSAR视线向速率, 研究区的速率整体变大, 同时跨狭义海原断裂带两侧的速率差异也显著增加。狭义海原断裂带中西段西南侧的速率相对较大的区域也变得更为显著。与将GPS速度场沿着平行于断层方向分解后的内插结果(图2b)进行对比可以看出, 虽然由于基准不一致, 无法在量值上和InSAR结果进行直接比较, 但其反映的相对形变特征, 特别是狭义海原断裂带西南侧速率较大的区域和InSAR结果具有一定的相似性。2种形变结果的不谋而合, 进一步证实了这一形变特征的可靠性。
在震后和震间期断层可以发生蠕变行为, 一般可以通过断层近场位移梯度监测。为此, 我们沿着海原断裂走向每隔5km选取一条垂直于断层走向的剖面(图4a中沿断层的圆圈), 选取剖面内的InSAR视线向速率进行投影。通过结果发现, 除了海原断裂西段与老虎山断裂相邻的部位存在明显的蠕滑之外(图4a中红色圆圈), 狭义海原断裂现今整体上并不存在明显的蠕滑。同时我们还发现, 跨海原断裂带的震间变形总体上可以描述为走滑断层弹性半无限空间的二维反正切螺旋位错模型(Savage et al., 1973)。另外, 通过反复核对以及同GPS形变场进行对比, 认为图4a中局部的深蓝色变形为盆地内局部沉降变形所致, 对研究较大尺度构造变形的影响不大。
为了进一步比较狭义海原断裂带平行于断层方向的InSAR变形和GPS结果之间的差异, 我们设计了横跨狭义海原断裂带中段的剖面AA'(图4a中的黑色方框), 通过平行于断层方向的GPS速率和InSAR速率剖面(AA')来看(图6), GPS结果和InSAR结果总体上具有相当的一致性, 更加证实了InSAR观测结果的可靠性。
图b中灰色小点为平行于断层方向的InSAR观测速率, 带误差棒的红色圆点为将沿剖面走向每隔5km方框内的数据去除大于2倍均方根误差的点后的数据平均值, 带误差棒的蓝色圆点为剖面内平行于断层方向的GPS速率
3.1.2 六盘山断裂带
六盘山断裂走向近SN, 与InSAR卫星的飞行方向几乎平行, 因此很难监测到平行于断层方向(SN向)的变形。六盘山以逆冲推覆运动为主, 兼具一定的左旋走滑分量, 构造变形主要以垂直断层方向的水平会聚和断裂以西的地壳垂直隆升变形为主, 现今水准资料研究表明跨六盘山断裂垂直形变速率梯度相当显著(Hao et al., 2014), 而GPS观测则显示平行断层方向的走滑速率较小。由于卫星飞行方向和视线方向的限制, InSAR本身对垂直方向变形信息更为敏感。因此, 六盘山地区的InSAR视线方向变形应当主要由垂直方向的变形引起。
从InSAR视线方向变形速率来看(图7), 六盘山断裂以西表现为1~2mm/a的轻微上升, 而西秦岭和陇县-宝鸡断裂带的交会部位以及鄂尔多斯和渭河盆地的交界处则为2~3mm/a的下沉。鄂尔多斯地块内部为显著的上升, 相对渭河盆地达到6~7mm/a, 这与邓起东等(1999)提出的鄂尔多斯相对周围盆地平均呈6~8mm/a的速度上升一致。与3~10a时间尺度(站点观测时间的差异较大)GPS流动站垂向观测速率进行比较, 尽管存在点位密度分布的巨大差异和局部形变的不匹配性, 但我们仍然可以看出, 总体而言, 六盘山断裂以西的青藏高原一侧呈1~2mm/a的速度上升与InSAR结果具有一致性。和40a时间尺度的水准观测结果进行比较可知, 水准观测得到六盘山断裂以西呈显著的局部地壳隆升, 而这在我们解算的InSAR观测结果中并没有出现。由于地壳垂向运动本身具有波动性, 时间尺度不完全一致可能会造成不同结果之间的差异。特别是在地震孕育晚期, 由于震源区地壳介质的弹性形变已趋于极限, 形变会显著减弱。
3.2 InSAR和GPS融合后的高密度三维形变场
对比InSAR和GPS融合后的高密度水平形变场EW向速率和SN向速率(图8a, b)以及相应单独使用GPS水平速度场内插后的形变场(图8c, d)可以看出, InSAR-GPS融合水平速度场的长波信号总体和GPS结果一致, 这能够为和构造变形相关的距断层中长距离的水平变形提供可靠约束。由于InSAR近极轨道观测对SN向变形通常不敏感, 因此, 融合后速度场的SN向分量主要吸收了GPS数据的绝大部分信息。由于EW向分量融合了高分辨率的InSAR数据, 可对在研究区内, 特别是断层附近的一些变形细节表现得更为明显。为了进一步与前述平行于海原断裂方向的GPS形变场进行比较, 将融合后的速度场沿着狭义海原断裂的走向进行投影, 与平行于断层方向的GPS速度相比可以看出, 融合后的速度场继承了GPS速度场和InSAR变形的显著特征, 同时形变场受单个GPS站点的影响也显著减小(图9a)。在三维形变场融合时, 仅以GPS水平形变场作为约束, 而不加入GPS垂向约束, 则在六盘山构造区得到的垂向变形大体和InSAR视线向结果一致, 如果同时加入Liang等(2013)得到的GPS垂向速率结果作为约束, 则得到的垂直形变场在六盘山以西呈现出明显的局部隆升特征(图9b)。
4 讨论
4.1 沿着狭义海原断裂带现今的断层滑动速率变化
断层滑动速率的研究对于认识区域构造变形、分析断裂强震危险性等均具有重要意义。针对海原断裂带现今的滑动速率, 国内外大量学者借助GNSS观测获得了相当多的结果, 但由于GPS观测密度相对稀疏及所用模型不同, 得到的结果之间存在一些明显的差异。例如, 戴洪宝等(2016)通过剖面分析法认为海原断裂西、中和东段的滑动速率分别为4.0mm/a、4.0mm/a和2mm/a; 崔笃信等(2009)和张晓亮等(2011)分别采用Smith 3D体力模型和三维非线性有限元方法反演得到海原断裂西、中和东段的滑动速率分别为5.0mm/a、4.0mm/a和1.0mm/a; 李强等(2014)和胡亚轩等(2009)分别借助块体模型和位错模型得到海原断裂西、中和东段的滑动速率分别为4.0mm/a、5.6mm/a和1.0mm/a。也有学者利用时序InSAR资料计算得到海原断裂带现今的滑动速率(乔鑫等, 2019), 但由于所采用的影像覆盖范围有限, 并未得到海原断裂整体、精细的滑动速率结果。为此, 我们以InSAR和GPS融合后的高密度水平形变速度作为约束, 借助二维弹性半空间螺旋位错模型, 沿着海原断裂带的走向进行不同分段连续、精细的滑动速率反演。
由于海原断裂现今以走滑为主, 本文将海原断裂震间期的变形假设为一般倾角的弹性半空间走滑断层螺旋位错模型, 借助InSAR和GPS融合后的高密度水平形变场, 用该模型来研究海原断裂沿断层走向的滑动速率。
根据已有的研究结果, 考虑断层深、浅部特征, 将海原断裂的倾角范围限定在70°~90°之间。采用上述公式, 通过垂直海原断裂走向剖面内的InSAR-GPS融合水平形变场数据计算现今沿海原断裂断层走向的滑动速率。如图1所示, 根据断裂实际展布, 将狭义海原断裂分为4段(S1-S4), 沿着每个分段断层的走向提取剖面, 走向每变化一次则提取一个剖面, 剖面垂直于断层走向, 长200km(距断层两侧各100km), 每个剖面平均宽约5km。
在计算中发现, 由于在海原断裂西南侧相对断层近场(距断裂约50km)处存在可能因1920年海原地震导致的震后变形, 使得模型对海原断裂西南侧的拟合效果不佳。为此, 我们采用2种方案来考虑这一问题: 第1种方案采用所提取剖面内距断层(南侧)50km以内的数据进行拟合, 以获得断层滑动速率; 第2种方案采用所提取剖面内距断层(南侧)50km以外的数据进行拟合, 以获得断层滑动速率。从2种方案对观测数据的拟合效果来看, 尽管有个别剖面的拟合效果不佳, 但拟合残差总体较小, 基本能反映沿海原断裂走向的断层滑动速率变化的细节特征。
从反演得到的滑动速率来看(图10), 整体上采用第2种方案获得的断层滑动速率比采用第1种方案获得的断层滑动速率更小。但二者每个分段的均值差异≤1mm, 与不同学者仅利用GPS形变场研究得到的每个分段断层滑动速率之间的差异相差不大。从利用2种方案获得的海原断裂沿断层走向的滑动速率来看, 总体上自西向东存在滑动速率递减的现象, 且这种现象不仅存在于不同分段断层之间, 也存在于同一分段断层之内, 但在同一分段断层之内速率自西向东递减的特征不如不同分段断层之间的明显。从第2种方案的结果来看, 海原断裂不同分段的平均滑动速率自西向东依次为4.9mm/a(S1)、4.8mm/a(S2)、3.4mm/a(S3)和1.4mm/a(S4段), 滑动速率的递减主要发生在中东段, 而西段递减并不显著。对比仅利用GPS形变资料得到的海原断裂带的滑动速率结果可知, 尽管存在断裂分段特征和研究模型的差异, 但总体上无论是分段滑动速率大小, 还是速率自西向东的递减特征, 均与本文的研究结果具有相似性, 区别仅在于由于本文采用高密度InSAR-GPS融合水平形变场, 得到的滑动速率在空间尺度上更为精细。
4.2 六盘山断裂带现今的断层运动与强震危险性分析
杜方等(2018)的研究表明, 六盘山断裂带的中-南段至少于距今约1 400a内未发生过M≥6.5的强震, GPS形变显示已有显著的应变积累, 出现小地震稀疏或空缺的部位及低b值区, 表明六盘山中-南段已经闭锁并具备高应力积累, 是未来潜在的强震发震区。另外, 跨六盘山断裂中-南段的精密水准观测也显示该区的地壳垂直形变具有和汶川地震前跨龙门山断裂类似的特征(Hao et al., 2014)。
GPS和水准观测资料表明, 六盘山断裂现今以逆冲运动为主, 兼具少量的左旋走滑分量。因此, 该断裂现今的变形特征以垂直断裂方向的水平挤压变形和垂向变形为主。本文将六盘山断裂的震间水平变形假设为一般倾角的弹性半空间倾滑断层螺旋位错模型, 使用跨六盘山断裂中-南段的InSAR-GPS融合后的高密度水平形变场垂直断裂方向的剖面数据(图7中的剖面BB')研究六盘山断裂中-南段现今的倾滑速率和闭锁深度。考虑到逆冲断裂垂直断层方向的水平变形主要集中在上盘, 即六盘山以西的青藏高原一侧, 我们在断裂以西取距断裂100km范围内的数据; 而在以东的鄂尔多斯一侧, 考虑到鄂尔多斯地块西南部地壳本身存在差异变形, 仅取距断裂30km以内的数据。
我们利用InSAR-GPS融合后的高密度水平形变场, 通过螺旋位错模型反演得到六盘山中-南段断层现今的逆冲倾滑速率为2.1mm/a, 搜索得到断层倾角为41°, 断层的闭锁深度为20km, 已到达上地壳孕震层的下界(小地震分布层深约20km), 与已有的关于六盘山断裂现今运动和闭锁特征的研究结果基本一致(李强等, 2014)。根据断层面的滑动速率S、断层倾角和断层长度(计算断层面积A),可由公式
计算潜在发震断层的地震矩累积率。其中, μ为剪切模量, 一般取30GPa。进一步通过地震矩和矩震级的关系MW=
logM0-6.03估计潜在发震断层未来潜在地震的震级。根据上述反演得到的六盘山断裂现今的滑动速率、倾角和闭锁深度, 以及杜方等(2018)给出的6.5级以上地震的离逝时间(约1 400a)和断层长度范围(70~80km), 估算认为六盘山断裂中-南段未来潜在地震的最大矩震级MW约为7.5, 与杜方等(2018)仅根据断层闭锁长度经验公式得到的约 MW7.3 相当。
4.3 1920年海原8.5级地震现今的震后变形
相关研究(Rice et al., 1983; Lorenzo-Martín et al., 2006; Jiang et al., 2020; 刘雷等, 2021)表明, 一次大地震由于介质黏弹性效应引起的震后变形持续的空间范围和时间长度由地球介质的流变结构控制。通过现今的GPS和InSAR形变场可以看出, 虽然1920年海原大地震距今已有100a, 但在海原断裂南侧仍然可以观测到类似震后地壳介质黏弹性松弛引起的变形特征。从跨海原断裂平行于断层方向的地壳运动速率中可以看出, 这种震后变形主要存在于断裂南侧, 具有明显的非对称性。从半无限空间弹性螺旋位错模型的拟合效果来看, 南侧形变速率的拟合效果较差, 特别是靠近断层附近, 而北侧的拟合效果则较好。其中, 南侧的变形特征和Savage(2000)与Devries等(2013)指出的中下地壳黏弹性介质震后松弛引起的变形非常类似。至于产生断裂两侧变形非对称性的原因可能包括以下2个方面: 1)断裂两侧中-下地壳的地球介质流变参数不同, 即海原断裂南、北2侧的中-下地壳介质波速结构存在显著差异(莘海亮等, 2020)。2)考虑到海原断裂南侧青藏高原一侧相对于北侧阿拉善一侧是主动盘, 故与2001年昆仑山口西8.1级地震类似(李海兵等, 2004), 1920年海原大地震的破裂方式应该也以单侧块体运动为主。断层面的南、北2盘相对震前的同震静态位错是不对称的(Jiang et al., 2020), 从而使得震后变形也具有非对称性。
如果认为海原断裂南侧现今的地壳形变叠加了1920年海原大地震由于地壳介质黏弹性效应引起的震后变形, 采用现今的观测结果研究海原断裂长期滑动速率就可能存在问题。但从使用相对断层的远场数据(与断层的距离>50km)和相对断层的近场数据分别拟合的结果来看, 得到的断层滑动速率相差不大, 在1mm/a范围内。而已有的利用GPS水平速度场对海原断裂现今的滑动速率开展的研究中, 大部分结果虽然平均拟合残差较小, 但对前述海原断裂南侧运动相对较大的GPS站点的拟合效果并不好(崔笃信等, 2009; 石富强等, 2018), 残差相对其他站点更大, 这也从侧面说明仅采用震间位错模型并不能完全解释海原断裂远、近场的变形特征, 尽管这对断层滑动速率的研究影响有限。
由地壳介质黏弹性效应引起的震后变形不仅与同震断层面静态滑动、破裂方式有关, 也与断层两侧地壳介质的性质有关, 很难用简单的二维模型(Savage, 2000)获得较好的模拟效果。另外, 由于观测技术发展的客观原因, 至今仍无法获得1920年海原大地震的同震断层面精细破裂模型, 因此上述有关1920年海原大地震引起的震后变形仍然处于推测阶段, 需要更多的观测资料和模拟研究加以证实。
4.4 海原-六盘山构造区现今的应变率场特征
与速度场相比, 应变率场能够直接给出地壳变形量的大小及空间差异信息, 能更好地反映断层应变率的积累特征, 从而揭示变形与构造之间的对应关系。由于GPS站点空间分布密度较低, 计算得到的应变率场空间的分辨率较低, 无法更多地揭示地壳变形的细节。因此, 我们依据上面得到InSAR-GPS高密度地壳水平运动场, 采用Shen等(2015)提出的优化最小二乘配置法计算了研究区的高分辨率应变率场, 借助不同应变率场参数的空间分布研究海原-六盘山构造区现今构造变形和断层应变积累的细节特征。需要说明的是, 由于我们的观测时间与地震复发周期相比很短, 除当前区域构造荷载外, 断层当前的地震周期也会影响断层周围的应变率分布。1920年发生在海原断裂上的8.5级大地震的震后变形虽距今已有100a, 但其震后影响仍然可能对现今的应变率场产生影响。
图 12 分别为利用InSAR-GPS融合形变场得到的海原-六盘山地区的面应变率、最大剪应变率、主应变率场和旋转应变率场的结果, 与单独利用GPS速度计算的结果相比, 应变率场的空间分布细节更为清楚。研究区的最大剪应变率(图12a)高值区主要集中在海原断裂及其附近, 表明海原断裂是该区域现今构造变形最为活跃的走滑断裂, 是调节青藏高原NE向扩张的重要边界断裂, 这与大多数地质研究的认识一致(Zhang et al., 1990)。海原断裂与香山-天景山断裂之间为青藏高原与阿拉善地块之间的过渡地带, 其内最大剪应变率也较为明显, 形成了一个显著的左旋剪切带。而六盘山断裂附近的最大剪应变率相比海原断裂附近要小得多, 表明六盘山断裂现今的走滑剪切变形相对较弱。整体上, 研究区的面应变率(图12b)在青藏高原一侧大部分表现为压性特征, 而在鄂尔多斯一侧表现为轻微张性特征, 这和基于块体的应变研究具有一致性(蒋锋云等, 2013), 而面应变率变化显著的部位则主要沿着海原-六盘山构造带展布。需要强调的是, 在六盘山断裂的中段和南段以西分别存在2个显著的面应率压性高值区, 与前文提到的六盘山断裂以西的地壳挤压缩短隆升存在一定的对应关系。另外, 海原断裂带中东段、六盘山断裂带以西2个压缩区之间及断层附近呈现局部张性特征。从主应变率特征来看, 与最大剪应变率类似, 主应变率的高值区主要分布于海原断裂附近与香山-天景山断裂之间, 且主压应变率方向与海原断裂、香山-天景山断裂中西段走向的夹角近45°, 表现为典型的左旋走滑特征。而对于香山-天景山断裂东南段、烟筒山断裂、云雾山断裂等近SN走向的断裂, 主压应变率的方向近垂直于断裂走向, 主要表现为逆冲特征。相比海原断裂, 六盘山地区的主应变速率整体较弱, 仅在青藏高原一侧靠近六盘山断裂附近以西的位置相对较强, 且主压应变率方向沿着断裂走向具有明显的差异, 六盘山断裂中段以垂直断裂方向为主, 而在六盘山断裂南段主压应变率的方向较为凌乱, 但整体上并非垂直于断裂, 和仅利用GPS获得的大尺度主压应变率方向近垂直于六盘山断裂的情况存在差异(江在森等, 2003; 武艳强等, 2020)。另外, 香山-天景山断裂东南段附近的面应变率、主应变率和最大剪应变率相对周围明显更小, 这是否与该断层的闭锁导致周围介质硬化有关(刘雷等, 2020), 需要结合其他资料开展进一步的深入研究。邹镇宇(2015)经过公式推导显示, 微元的旋转与剪应变不同, 不会随着所取方向改变而发生变化, 因此地壳旋转率在反映走滑断层变形特征方面具有独特优势。从旋转率场(图12c)结果来看, 研究区分为多个旋转变形单元, 这些旋转变形单元以走滑断层为界。其中, 陇西地块以顺时针旋转为主, 鄂尔多斯地块西南部包括六盘山构造带以轻微的逆时针旋转为主, 鄂尔多斯西缘右旋剪切带以顺时针旋转为主, 海原断裂与香山-天景山断裂之间的左旋剪切带则以逆时针旋转为主。
5 结论
本文在分析InSAR、GPS观测得到的海原-六盘山地区现今的地壳形变场特征的基础上, 进一步采用二者融合后的速度场研究了海原-六盘山地区现今的地壳形变细节, 主要获得以下认识:
(1)1920年海原8.5级地震的发震断裂南侧存在明显的地壳快速运动区域。该区域可能是1920年海原8.5级地震引起的地壳介质长期黏弹性松弛效应所致。
(2)狭义海原断裂带现今的左旋走滑速率自西向东存在减少特征, 但滑动速率的递减主要发生在中东段, 而中西段递减并不显著。从西段和中西段的近5mm/a到中东段的3.5mm/a, 再到最东段的1.4mm/a, 变化幅度较大。这可能与海原断裂的左旋构造系统与六盘山断裂的逆冲推覆构造系统之间的构造转换有关。
(3)从水准和GPS反映的跨六盘山断裂中-南段的垂直形变速率剖面来看, 六盘山断裂中-南段具有和汶川地震之前龙门山地区地壳垂直隆升类似的变形特征。从InSAR和GPS融合后的高密度地壳水平形变跨六盘山断裂中-南段垂直于断裂的会聚速率剖面拟合结果来看, 六盘山断裂目前处于完全闭锁状态, 闭锁深度约达20km, 远场会聚速率约为2.1mm/a, 估计最大矩震级约达7.5级。
(4)从高密度形变场、应变率场和旋转变形场结果来看, 海原断裂表现为以左旋走滑为主的运动特征, 跨断裂变形差异显著, 最大剪切应变率尤为显著, 是青藏高原NE向扩张的重要边界断裂。海原断裂与香山-天景山断裂之间为青藏高原与阿拉善地块之间的过渡地带, 其内应变率也较为显著, 左旋剪切变形强烈, 需要注意的是香山-天景山断裂东南段出现和周围应变场不协调的特征。六盘山断裂则主要表现为以逆冲推覆构造为主的运动特征, 面应变率在六盘山断裂以西整体上表现为明显的挤压变形特征。从旋转变形来看, 研究区表现出多个旋转变形单元, 其中逆时针旋转对应左旋走滑变形(海原断裂与香山-天景山断裂间的左旋剪切带), 而顺时针旋转对应右旋走滑变形(鄂尔多斯西缘右旋剪切带和陇西地块)。
（内容有删减，可文后↙“阅读原文”）
图10
引用格式
蒋锋云, 季灵运, 朱良玉, 刘传金. 联合GPS和InSAR研究海原-六盘山断裂现今的地壳变形特征[J]. 地震地质, 2023, 45(2): 377-400 DOI:10.3969/j.issn.0253-4967.2023.02.005
JIANG Feng-yun, JI Ling-yun, ZHU Liang-yu, LIU Chuan-jin. THE PRESENT CRUSTAL DEFORMATION CHARACTERISTICS OF THE HAIYUAN-LIUPANSHAN FAULT ZONE FROM INSAR AND GPS OBSERVATIONS[J]. Seismology and Geology, 2023, 45(2): 377-400 DOI:10.3969/j.issn.0253-4967.2023.02.005
作者简介:蒋锋云，男，1978年生，高级工程师，主要从事地壳形变机理与地震长期预测研究，E-mail:jfy267862@163.com。
基金资助:
国家自然科学基金(42104061);国家自然科学基金(41874117);陕西省自然科学基础研究计划项目(2023-JC-QN-0292)
转自：“测绘学术资讯”微信公众号
如有侵权，请联系本站删除！
