青藏高原东南缘是我国现今地壳构造形变与地震活动非常强烈的地区之一^[1-2]，区域内发育多条深大活动断裂^[3-5]，曾发生过多次中强地震，其中危害最严重的是2008年汶川M_S8.0强震^[6-7]。开展青藏高原东南缘现今地壳构造形变，特别是强震活动影响下的地壳形变与应变场变化研究，对于深入认知青藏高原现今地壳构造活动特征与地震危险性评估具有重要意义。

大范围尺度下的高精度GNSS观测资料为研究地壳构造活动与动力学机制提供了有力的基础参考信息。相关学者基于青藏高原东南缘GNSS监测资料，通过计算库仑应力变化获取汶川地震影响范围^[8-9]，并分别利用块体形变模型^[10]、位错模型^[11]、块体模型结合最小二乘配置模型^[12]及块体应变模型^[13]等探讨汶川强震对周边块体运动特征与龙门山断裂的影响，为认知汶川强震同震形变及震后短期(1~2 a内)效应对青藏块体东南缘地壳活动影响提供了重要参考。

本文利用青藏高原东南缘在汶川地震前8 a(1999~2007年)与震后6 a(2011~2017年)的高精度GNSS监测资料，通过构建区域地壳运动与应变多尺度球面小波模型，获得研究区现今地壳形变与多尺度应变率场特征，分析汶川强震前后不同空间尺度下研究区应变率场差异性特征，并进一步深入探讨在何种尺度范围下模型计算结果是合理的及在何种尺度下模型计算结果最能反映出区域构造形变的真实情况。本文成果可为研究青藏高原东南缘构造活动及其演化规律提供有价值的参考信息。

1 青藏高原东南缘现今地壳运动及变化特征

本文GNSS监测数据来自中国地壳运动观测网络与中国大陆构造环境监测网络。数据处理流程概述如下：采用GAMIT软件单日松弛解处理基线^[12]，采用IGS提供的sp3精密卫星星历削弱轨道误差，采用模型改正固体潮和海潮对GNSS台站的影响，引入中国大陆及周边IGS站进行联合解算获得基线单日松弛解h文件，并将h文件代入GLOBK软件进行统一平差计算^[13]，获得ITRF2008框架下速度场^[14]，并进一步计算获得研究区域相对于欧亚板块的速度场(图 1)。需要说明的是，为顾及7.0级以上强震对地壳速度场的影响，采用相应同震模型对2001年可可西里M_W8.1^[14]、2011年日本宫城M_W9.0^[15]与2013年芦山M_W7.0地震^[16]进行处理。

图 1(a)和1(b)均显示，青藏高原东南缘地壳运动整体呈现出顺时针旋转的特征，并在龙门山、鲜水河、安宁河、则木河与小江断裂带处形成明显的速度差异梯度带，且该速度梯度带西侧GNSS地壳运动速度场量值显著大于东侧。研究区2011~2017年GNSS速度场量值整体大于1999~2007年，1999~2007年研究区EW向平均速率约为6.7 mm/a，NS向平均速率约为－5.8 mm/a；2011~2017年研究区EW向平均速率约为8.0 mm/a，NS向平均速率约为－6.3 mm/a。图 1(a)和1(b)最显著的差异是汶川强震后龙门山断裂带附近西北侧地壳运动GNSS速度场量值显著增大，1999~2007年速率平均值约为8.7 mm/a，2011~2017年速率平均值约为14.3 mm/a；龙门山断裂带东南侧地壳运动速度量值则变化不大。2011~2017年羌塘、巴颜喀拉与川滇地块呈现出加速向南东运移并推动华南块体的趋势。上述结果表明，汶川强震后研究区内构造块体地壳运动差异性特征明显增强。

2 地壳运动与应变多尺度球面小波模型

球面小波技术具有对空间和频率局部化的功能，将其用于分析地壳运动时，不同的分解尺度会反映出不同空间范围内的地壳运动和形变特征信息，其中大空间范围反映的是区域的整体特征，小空间范围反映的是局部细节特征^[17]。本文选用高斯函数差(difference of Gaussians, DOG)构建的球面DOG小波作为球面小波模型基函数^[18]，其表达式为：

$ \begin{gathered} \psi_{o, a}\left(o^{\prime}\right)=\chi_a(\gamma)^{1 / 2} \exp \left[-\frac{\tan ^2(\gamma / 2)}{a^2}\right]- \\ \frac{1}{\alpha} \chi_{a a}(\gamma)^{1 / 2} \exp \left[-\frac{\tan ^2(\gamma / 2)}{(\alpha a)^2}\right] \end{gathered} $

(1)

式中，a=1/2^q, q为尺度因子；γ为观测点o′与小波中心极o的夹角；α为常数(α>1)；$\chi_a(\gamma)=\frac{4 a^2}{\left[\left(a^2-1\right) \cos \gamma+a^2+1\right]^2} $。由式(1)可知，尺度因子q越大，球面小波的基函数涵盖的空间范围越小。由球面小波基函数的线性组合构成的球面小波框架可以更加灵活地表示信号，将此球面小波框架用于估计GNSS速度场，则速度场v(λ，φ)可表示为：

$ \begin{gathered} v(\lambda, \varphi)=\sum\limits_{k=1}^M\left[a_k g_k(\lambda, \varphi) \hat{r}+\right. \\ \left.b_k g_k(\lambda, \varphi) \hat{\lambda}+c_k g_k(\lambda, \varphi) \hat{\varphi}\right] \end{gathered} $

(2)

式中，λ为经度；φ为纬度；M为框架函数个数；g_k(λ，φ)为第k个球面小波基函数；$ \hat{r} $为垂向；$ \hat{\lambda} $为东方向；$ \hat{\varphi} $为北方向；a_k、b_k、c_k分别为3个方向的待定系数。由于g_k(λ，φ)为非正交基，容易导致法方程系数矩阵出现秩亏，使得待定系数的解不唯一，故采用正则化方法进行处理^[19]。

地壳运动监测中，GNSS站点通常分布不均，因此需根据站点疏密程度选择相应的尺度因子^[20-21]。基于多尺度球面小波模型构建原理，根据研究区域范围的球面积得到区域半径L，最小尺度因子下球面小波基函数涵盖的空间范围应小于2L，不同尺度因子q下球面小波的基函数涵盖的空间范围如表 1所示。根据本文研究区域GNSS站点分布范围，可计算得出L≈1 400 km，因此取最小尺度因子q=3；根据满足建立球面小波模型的基本条件，即最大尺度因子下球面小波模型的每个小波点涵盖的空间范围内应至少包括3个GNSS测站点，可确定最大尺度因子q=8。图 2为根据本文研究区域内GNSS站点的空间分布情况确定的不同尺度因子q的情况。

相对于欧亚板块下的区域GNSS速度场包含区域整体旋转运动，这一整体旋转运动可通过欧拉旋转矢量进行模拟^[22]。本文在进行地壳应变解算前，先将研究区域GNSS水平速度场进行欧拉旋转拟合，得到欧拉旋转角速率为0.276°/Ma，欧拉极为(14.603°S, 89.170°W)，再将两期GNSS水平速度场减去该欧拉旋转矢量，得到研究区域无旋转GNSS水平速度场，并将实测GNSS速度场与模型模拟值进行对比分析(图 3)。从图 3(a)和3(d)可看出，无旋转GNSS水平速度场可较好地反映出区域地壳相对运动特征，羌塘、川滇块体保持顺时针旋转运动趋势，华南块体整体表现出北北西向运动趋势，巴颜喀拉块体整体运动趋势一致性较差；安宁河-则木河-小江断裂带两侧地壳运动GNSS速度场量值在2011~2017年明显增强；龙门山断裂两侧地壳运动GNSS速度场在2011~2017年不仅在量值上有显著增大趋势，且在运动方向上也表现出明显的差异性。图 3(b)、3(c)、3(e)、3(f)分别展示了东向与北向速度分量模拟值与实测值残差统计直方图。可以看出，绝大部分残差趋于0且近似服从正态分布，1999~2007年东向与北向模拟值与观测值之差标准差分别仅为0.63 mm/a和0.64 mm/a，2011~2017年分别仅为0.53 mm/a和0.50 mm/a, 均小于各时间段GNSS速度场观测误差。

此外，为防止构建的多尺度球面小波模型在确定球面小波函数空间分布时出现过拟合现象，基于研究区域GNSS监测站点分布分别计算1999~2007年与2011~2017年多尺度球面小波模型东向与北向最优正则化参数(图 4)。

进一步将无旋转GNSS水平速度场作为解算球面小波模型应变参数的输入数据，采用球坐标系计算水平向应变张量D：

$ \begin{array}{l} \mathit{\boldsymbol{D}} = \left[ {\begin{array}{*{20}{l}} {{\varepsilon _{\theta \theta }}}&{{\varepsilon _{\theta \lambda }}}\\ {{\varepsilon _{\lambda \theta }}}&{{\varepsilon _{\lambda \lambda }}} \end{array}} \right] = \\ \left[ {\begin{array}{*{20}{c}} {\frac{1}{r}\frac{{\partial {v_\theta }}}{{\partial \theta }}}&{\frac{1}{{2r}}\left( { - {v_\lambda }\cot \theta + \frac{1}{{\sin \theta }}\frac{{\partial {v_\theta }}}{{\partial \lambda }} + \frac{{\partial {v_\lambda }}}{{\partial \theta }}} \right)}\\ {{D_{12}}}&{\frac{1}{r}\left( {{v_\theta }\cot \theta + \frac{1}{{\sin \theta }}\frac{{\partial {v_\lambda }}}{{\partial \lambda }}} \right)} \end{array}} \right] \end{array} $

(3)

式中，v_θ=－v_N，v_N为速度场北向分量；v_λ=v_E，v_E为速度场东向分量；r为地球半径；θ为余纬；λ为经度；D₁₂=D₂₁表示矩阵元素对称。由于在球坐标系中式(3)以南向为正，而GNSS水平速度场常以北向为正，因此通过坐标旋转分别得到北向、北东向、东北向和东向应变分量ε_NN、ε_NE、ε_EN、ε_EE，其中，ε_NE=ε_EN，进一步可计算出最小主应变率、最大主应变率、最小主应变率方位角、面膨胀率和最大剪应变率^[2]。

3 青藏高原东南缘现今地壳多尺度应变率场变化特征

为揭示青藏高原东南缘在汶川强震前后不同尺度下的应变场时空演化特征，利用球面小波模型分别计算获得研究区在1999~2007年和2011~2017年多尺度应变特征参数。

3.1 主应变率和面膨胀率变化特征

图 5为研究区1999~2007年与2011~2017年在尺度因子q=3~5、6、7、8时的主应变率与面膨胀率分布。可以看出，不同尺度因子下研究区应变特征表现出较显著的差异，不同空间与尺度范围内应变积累与变化特征也反映出青藏高原东南缘复杂的地壳构造运动与形变特性。

图 5(a)和5(b)为尺度因子q=3~5时的计算结果，主要反映的是区域整体应变场特征。可以看出，研究区主压应变方向整体以NW-SE向为主，在川滇块体南部主压应变率呈近EW向；羌塘、巴颜喀拉与川滇块体内部主应变率量值显著大于华南块体内部。对比1999~2007年与2011~2017年的结果看出，羌塘与川滇块体内部呈现出面膨胀特征，而面压缩特征主要位于巴颜喀拉块体内，不同的是，2011~2017年巴颜喀拉块体与华南块体交界的龙门山断裂带处也表现出一定的面压缩特征。

图 5(c)和5(d)为尺度因子q=6时的计算结果。可以看出，q=6时应变率结果较q=3~5时体现出更多的细节差异，最显著的特征为1999~2007年主压应变率主要位于鲜水河(与龙门山断裂交界处)-安宁河-则木河-小江断裂带处，最大量值约－21×10^－9/a，方向以NW-SE向为主，且该区域也表现出显著的面压缩状态；2011~2017年主压应变率虽仍然沿鲜水河-安宁河-则木河-小江断裂带处展布，但主压应变最显著的区域则位于龙门山断裂带及其附近区域，最大量值约－29×10^－9/a，同时龙门山断裂带也是面压缩最为显著的区域，量值高达约30×10^－9/a。

图 5(e)和5(f)为尺度因子q=7时的计算结果。可以看出，研究区域应变场特征较q=6时反映出更丰富的局部细节特征。2个时间段内主应变率和面膨胀率整体分布特征较为类似，地壳面膨胀和压缩区域整体表现出成对轮流出现的特点，但各部分形变特征并不完全一致。该尺度下的应变特征主要反映研究区域内局部连续形变特性，进一步展示出青藏块体东南缘构造形变的广泛分布和复杂性。1999~2007年和2011~2017年主应变率与面膨胀率最显著的差异特征仍主要体现在龙门山断裂带附近，1999~2007年龙门山断裂带北西侧主要表现为轻微面膨胀，东南侧则表现出微弱面压缩状态；而2011~2017年龙门山断裂带两侧均表现出显著的NW-SE向主压应变与面压缩状态。

图 5(g)和5(h)为尺度因子q=8时的计算结果。可以看出，其所反映的区域主应变和面膨胀分布与尺度因子q=3~5、6、7时差异较大，有些地方甚至出现相反的结果。如当q=8时，2011~2017年龙门山断裂带处并未呈现出较显著的面压缩(压应变)特征，反而呈现出一定的面膨胀特征。而根据区域构造活动特性，2008年汶川地震后，在龙门山断裂带附近应呈现出较显著的面压缩(压应变)特征。据此，本文舍弃尺度因子q=8，保留合理的最大尺度因子q=7。

从上述分析可看出，当尺度因子q=6时，即可较好地揭示出汶川强震后龙门山断裂带附近显著的主压应变与面压缩高值现象。从图 5(c)也可发现，研究区域的主应变中面膨胀和面压缩区的高值区主要沿川滇块体的东西边界交替分布，这一交替区域也是地壳运动速度场显著差异带，易造成区域应变积累。而汶川地震后巨大的应力释放打破了区域原有的应力平衡，在应力再平衡的过程中，研究区面膨胀和压缩区高值区转换为沿龙门山-安宁河-则木河-小江断裂带对称分布(图 5(d))，表明在尺度因子q=6时，研究区地壳形变受边界断裂带控制。同时，图 5(e)、5(f)中主应变率与面膨胀率正负值轮流出现的特征也精细地反映出研究区复杂的地壳运动形变特性与地质构造背景、地貌特征；面膨胀与面压缩相间区域反映出研究区地壳岩性不均匀的特征，可能与区域内广泛发育的波浪形褶皱有关，反映出了局部区域应力平衡的力学机制。

进一步，为充分发挥出多尺度球面小波模型在表达区域复杂地壳形变特征中的优势，更加清晰地综合揭示出区域地壳形变大尺度(整体)与局部形变特征，分别计算尺度因子q=3~5、3~6、3~7下球面小波模型模拟GNSS速度场与实测GNSS速度场残差的RMSE量值(表 2)。表 2(单位mm/a)显示，当尺度因子q=3~7时，模型模拟与实测速度场RMSE量值最小，可较好地综合揭示出区域地壳大尺度(整体)与局部形变特征，即最佳尺度因子为q=3~7。

图 6为q=3~7时研究区主应变率与面膨胀率分布结果。可以看出，2个时间段中研究区均呈现出以龙门山、安宁河、则木河、小江断裂带为界，主应变率总体表现出西强东弱的趋势；甘孜-玉树断裂南段、鲜水河、安宁河、则木河、小金河-丽江、红河、金沙江和小江断裂带处是主应变率与面膨胀率高值区；甘孜-玉树断裂南段、鲜水河、安宁河、则木河、金沙江断裂带处均呈现出NW-SE向挤压应变与NE-SW向拉张应变的特征；小金河-丽江、红河、小江断裂带处呈现出近EW向挤压应变与近SN向拉张应变的特征。

图 6(a)显示，1999~2007年研究区内最大挤压应变率高值区位于甘孜-玉树断裂南段、鲜水河、安宁河、则木河和小江断裂带处，量值约为－70×10^－9/a。小金河-丽江与红河断裂带处挤压应变率平均量值约为－25×10^－9/a。面压缩特征也主要沿上述几条主干断裂带展布。巴颜喀拉和华南块体内主应变率与面膨胀率数值相对均较小。龙门山断裂带整体主应变率量值较小，在龙门山断裂带最南端与鲜水河及安宁河断裂带交界处呈现出一定NW-SE向挤压应变与面压缩状态。图 6(b)显示，2011~2017年研究区内整体主应变率与面膨胀率量值较1999~2007年均有所增加，特别是沿着甘孜-玉树、鲜水河、安宁河、则木河、龙门山、小江、小金河-丽江与红河断裂带处尤为突出。2011~2017年研究区域内最大挤压应变率高值区位于龙门山断裂带及其附近区域，断裂带区NW-SE向挤压应变率量值高达约－91×10^－9/a，该区域也是整个研究区域内面压缩最大量值区，面压缩量值高达－72×10^－9/a。甘孜-玉树断裂南段、鲜水河、安宁河、则木河和小江断裂带处呈现出更显著的挤压应变特征，挤压应变量值平均约为－45×10^－9/a。巴颜喀拉块体东部与华南块体西部也呈现出较强地震前显著的挤压应变特征。

3.2 最大剪应变率变化特征及与地震活动性关系

图 7为研究区1999~2007年和2011~2017年在尺度因子q=3~5、6、7时的最大剪应变率分布特征。可以看出，研究区最大剪应变率在不同尺度因子下虽呈现出细节差异，但最大剪应变率高值区几乎都沿主干断裂展布。

图 7(a)和7(b)为研究区在尺度因子q=3~5时的计算结果，反映的是大尺度范围内区域最大剪应变率整体分布特征。可以看出，最大剪应变率高值区在整个研究区内未呈现出较多的细节性特征，但在1999~2007年区域内最大剪应变率最大值位于鲜水河断裂带北东段与甘孜-玉树断裂带交界处，量值约26×10^－9/a；而2011~2017年最大值则位于鲜水河与安宁河以及龙门山断裂带南段的交界处，量值约29×10^－9/a。

图 7(c)和7(d)为研究区在尺度因子q=6时的计算结果。可以看出，最大剪应变率较q=3~5时的结果有所不同，最显著的特征为相较于1999~2007年，2011~2017年时段龙门山断裂带呈现出明显的最大剪应变率高值区，量值约17×10^－9/a；在巴颜喀拉块体东部以及华南块体西部靠近龙门山断裂带处也呈现出较明显的最大剪应变率高值区，量值约15×10^－9/a。上述特征表明，汶川强震后龙门山断裂带及其附近区域构造剪切运动显著增强。

图 7(e)和7(f)为研究区在尺度因子q=7时的计算结果。可以看出，1999~2007年与2011~2017年研究区内主干断裂带剪切活动在汶川强震后有显著的增强趋势。

图 8为青藏高原东南缘1999~2007年与2011~2017年在尺度因子q=3~7时的最大剪应变率分布特征。可以看出，研究区最大剪应变率高值区主要分布于甘孜-玉树断裂南段、鲜水河、安宁河、则木河、小江、小金河-丽江与红河断裂带及其附近区域。表明主干断裂及其附近区域是青藏块体东南缘地壳构造活动最为剧烈的区域。同时，对比2个时间段研究区发生的M_S>4.0地震空间分布特征也可看出，上述主干断裂及其附近区域也是地震活动的高频带。

如图 8(a)所示，1999~2007年研究区最大剪应变率最大值位于甘孜-玉树断裂带与鲜水河断裂带交界处，量值约79×10^－9/a；除研究区域内几条主干断裂处呈现出最大剪应变率高值区，在巴颜喀拉地块中部也呈现出部分高量值区，但在龙门山断裂带及其附近区域最大剪应变量值很小。如图 8(b)所示，2011~2017年研究区主干断裂及其附近区域最大剪应变率量值大小和范围整体均较1999~2007年有所增加。虽然该时间段内最大剪应变最大值仍位于甘孜-玉树断裂带与鲜水河断裂带交界处及鲜水河断裂带南段，量值高达72×10^－9/a，但相较于1999~2007年该时间段最显著的差异特征是龙门山断裂带及其附近呈现出最大剪应变率高值区，量值达约50×10^－9/a。同时，羌塘、巴颜喀拉、川滇甚至华南地块西部均呈现出较1999~2007年更显著的地壳活动性。甘孜-玉树断裂南段、鲜水河、安宁河、则木河、小江、小金河-丽江与红河断裂带活动性进一步增强，此外，小金河-丽江与红河断裂带南部区域地壳活动性也进一步增强。从地震空间活动性来看，汶川强震对青藏高原东南缘区域内主干断裂带及其附近区域的地震活动性影响最为显著。

相比于1999~2007年，2011~2017年研究区沿着鲜水河、安宁河、则木河、小江、小金河-丽江与金沙江断裂带及其附近M_S>4.0地震明显增多，特别是龙门山断裂带及其附近区域地震沿断裂走向密集展布，且在华南块体西部与龙门山、安宁河与则木河断裂带交界处地震活动也明显增多。2个时间段地震活动的显著差异性进一步揭示出汶川强震对于研究区，特别是震中附近地壳构造活动性的显著影响。

图 9为1999~2007年与2011~2017年尺度因子q=3~7下青藏高原东南缘的面膨胀率和最大剪应变率差值分布。从图 9(a)看出，在距离汶川地震震中较近的龙门山、安宁河、则木河断裂带及华南块体西侧呈现显著的面压缩状态，反映汶川强震后巴颜喀拉地块与川滇块体加速向南东运移并推动华南块体的运动趋势；在其他主干断裂带附近也呈现出面应变率的高值区。图 9(b)也显示出整个研究区内的主干断裂带附近是最大剪切应变率的高值区所在，但最显著的特征还是在汶川地震震中的龙门山断裂带与鲜水河断裂南段呈现出的明显高值区。图 9所示的这些显著的应变率差异特征进一步表明，汶川强震对于研究区，特别是龙门山断裂带附近地壳构造活动性具有显著影响。

4 结语

1) 汶川强震前后研究区地壳运动形变具有一定的继承性发展特征，地壳运动整体呈现出向南东运移同时自身作顺时针旋转运动的特征，在龙门山、安宁河、则木河与小江断裂带处形成明显的速度差异梯度带，且高应变率值也主要位于主干构造断裂及其附近区域。

2) 不同尺度因子能够揭示出研究区域在不同空间范围内的地壳应变场积累特征。本文依据满足建立球面小波模型的基本条件，首先初步定量确定出最小与最大尺度因子，再根据不同尺度因子计算结果与区域构造活动特性的符合程度确定出合理的最大尺度因子，进一步通过计算不同组合尺度因子下模型模拟与实测地壳速度场残差的RMSE大小，最终确定能够较好地综合揭示出区域地壳形变大尺度(整体)与局部形变特征的最佳组合尺度因子。结果表明，对于青藏高原东南缘区域，合理的最大尺度因子q=7；当尺度因子q=6时，可较好地揭示出区域整体构造活动特性，即清晰地揭示出汶川强震后龙门山断裂带处呈现出的主压应变率、面压缩率与最大剪应变率高值特征；当尺度因子q=3~7时(最佳组合尺度因子)，可较好地综合揭示出区域地壳大尺度(整体)形变与局部形变特征。

需要说明的是，基于GNSS观测获得的应变率结果可较好地反映出区域现今构造应力场特征，但本文聚焦于分析青藏块体东南缘在不同尺度上对汶川强震发生前后的较长期响应特征，即使细致地划分尺度因子，对研究区域实际地壳形变复杂性的表征还是有限；且考虑到区域构造与形变是一个长期的地球动力学过程，GNSS观测也仅反映了现今地壳形变活动。此外，本文主要依据地壳活动程度分析地震活动与最大剪切应变率的空间对应关系，但地震的发生和发展极为复杂，特别是对于汶川强震震后效应的细致研究，还需结合震后粘弹性松弛等机理模型进行深入分析。

致谢: 感谢野外测量工作者获取GNSS观测数据，感谢中国地震局第二监测中心提供高精度GNSS数据。