郯庐断裂带是一条位于我国东部地区的最大规模的断裂带，其主要表现为右旋走滑^[1]。其内部构造复杂，历史上多次发生5.0级以上地震，其中1668年发生的郯城地震达到8.5级。为探究2011年日本宫城地震是否对郯庐断裂带南段的运动性质产生影响，对该地区的地壳形变与断层活动特征展开研究。李彦川等^[2]采用块体负位错模型分析郯庐断裂带中南段的地壳形变特征；李腊月等^[1]对不同时期的速度场利用负位错模型反演郯庐断裂带中南段的闭锁程度和滑动亏损速率，认为日本宫城地震缓解了郯庐断裂带中南段的应变积累。以往的研究中，采用的站点基本来源于中国大陆构造环境监测网络项目(CMONOC)，然而郯庐断裂中南段站点密度有限，无法揭示区域断层活动的细节特征。随着安徽连续运行参考站(Anhui continuously operating reference station，AHCORS)的建设与发展，增加了郯庐断裂带南段地区的站点密度，可为揭示区域活动断层的变形特征和地震危险性提供更详细的数据约束。

利用全球导航卫星系统(GNSS)的观测资料建立适当的数学模型，解算某一地区的地壳应变特征是研究地震危险性的重要手段。多尺度球面小波函数可以根据GNSS测站的密集程度和分布特征选择不同的尺度因子，揭示不同空间尺度下的应变特征^[3]。苏小宁等^[4]将多尺度球面小波解算GPS应变场的结果与位错模型结果对比，验证多尺度球面小波解算应变的有效性和稳健性；徐克科等^[5]提出多尺度球面小波的尺度和位置的选取、模型正则化因子以及参数估计的方法；李承涛等^[6]采用多尺度球面小波解算不同空间尺度下巴颜喀拉块体东北缘的应变率场，并分析其与2017年九寨沟M_S7.0地震的联系。本文基于2013-01~2018-06 AHCORS和CMONOC的同期共170个测站的观测数据，利用GAMIT/GLOBK解算其速度场，通过多尺度球面小波函数估计郯庐断裂带南段在不同空间尺度下的应变参数，揭示郯庐断裂带南段及邻近地区的地壳变形特征。研究结果可为该区域未来地震危险性评估提供参考。

1 GNSS数据处理

本文使用的数据主要是2013-01~2018-06期间37个AHCORS参考站的观测文件，利用GAMIT/GLOBK(Ver.10.7)软件对每隔7 d的观测数据进行处理。解算时引入研究区域以外的9个历史同期国际全球卫星导航定位服务(IGS)站(BJFS、DAEJ、IRKT、KIT3、LHAZ、POL2、TCMS、URUM、WUHN)，从而与AHCORS参考站构成一个区域网。

数据处理主要分为3个步骤：首先利用GAMIT软件解算AHCORS参考站以及周边IGS站的观测数据，获取单日松弛解，GAMIT单日解的解算策略如表 1所示；然后利用GLOBK联合SOPAC(Scripps Orbit and Permanent Array Center)机构发布的全球h文件进行网平差处理；最后通过坐标转换获取AHCORS参考站在ITRF2008框架下的坐标和速度。

采用双差相位解算方法，利用IGS的精密星历和地球旋转参数对其进行强约束，同时对测站坐标和卫星轨道进行解算。

采用GLOBK解算时间序列和速度场的处理策略如下：1)引入SOPAC的7个全球h文件(igs1~7)，将其与GAMIT基线解算得到的AHCORS参考站的h文件按照1.0的权重因子合并；2)利用GLOBK进行网平差，以IGS站在ITRF2008框架下的坐标和速度为基准，求得AHCORS参考站的速度。

由于ITRF2008框架以地球质心为基准，该框架下的点位运动速度包含了板块的运动速度，难以准确表达AHCORS参考站间的相对运动和形变，因此需要将ITRF2008框架下的运动速度转换到欧亚框架下。本文采用Wang等^[7]提出的相对于ITRF2008框架下的欧拉矢量(－0.087 mas/a，－0.514 mas/a，0.741 mas/a)，将AHCORS参考站在ITRF2008框架下的速度减去欧拉矢量所得到的速度，获得其在欧亚框架下的速度。AHCORS参考站在ITRF2008框架下的速度场如图 1(a)所示，在欧亚框架下的速度场如图 1(b)所示。从图 1(b)可以看出，安徽区域在欧亚框架下的速度主要往东偏南方向运动，与袁鹏等^[8]的结论一致。表 2(单位mm/a)为欧亚框架下AHCORS参考站在N、E、U方向的速度估值和中误差统计。

2 联合AHCORS参考站的多尺度球面小波应变场解算 2.1 多尺度球面小波基本原理

球面小波具有将某一给定空间进行局部化的特点，可以用来建立速度场和应变场的模型。其主要原理为，在半径为1的单位球上定义一个小波母函数，然后将母函数通过平移和伸缩得到小波基函数^[4]。本文选择高斯函数的差(difference of gaussians，DOG)形成的球面DOG小波作为基函数，其表达式为：

$ \begin{array}{l} {\psi _{x, a}}X' = {\lambda _a}{\left( \gamma \right)^{\frac{1}{2}}}{\rm{exp}}\left[ { - \frac{{{{\tan }^2}\left( {\frac{\gamma }{2}} \right)}}{{{a^2}}}} \right] - \\ \;\;\;\;\;\;\;\frac{1}{a}{\lambda _{{\alpha ^a}}}{\left( \gamma \right)^{1/2}}{\rm{exp}}\left[ { - \frac{{{{\tan }^2}\left( {\frac{\gamma }{2}} \right)}}{{\alpha {a^2}}}} \right] \end{array} $

(1)

式中，尺度a=1/2^q，q为尺度因子，q的值越大，尺度越小；γ为观测点与球面小波中心间的夹角，一般在0°~180°之间；α是一个大于1的常数，本文根据实验结果将α取值为1.25。

2.2 多尺度球面小波估计速度

在使用多尺度球面小波解算郯庐断裂带南段应变场前，需要根据GNSS测站分布的疏密情况来选择分解尺度。多尺度球面小波方法具有可将给定空间局部化的特点，可有效消除测站分布不均匀带来的影响。选择分解尺度的原则为：在测站分布密集的区域选择较大的分解尺度，在测站分布相对稀疏的区域选择较小的分解尺度。

本文根据AHCORS参考站和郯庐断裂带南段附近的CMONOC站点的疏密程度来选择分解尺度，如图 2所示，图中黑色实心小圆点为测站。从图 2可以看出，小波分解的最大分解尺度为6~9级，测站间距为21.8 km的9级分解尺度区域共有8个，其余大部分区域是测站间距为43.6 km的8级分解尺度。

在球面上二维GPS速度场的表达式为：

$ \nu \left( {\theta , \varphi } \right) = {\nu _\theta }\left( {\theta , \varphi } \right)\hat \theta + {\nu _\varphi }\left( {\theta , \varphi } \right)\hat \varphi $

(2)

式中，$ {\hat \theta }$、$ {\hat \varphi }$分别表示南北向、东西向测站速度的方向，θ和φ分别表示测站的纬度和经度。将平方可积空间的标量函数代入式(2)得：

$ \nu \left( {\theta , \varphi } \right) = \mathop \sum \limits_{k = 1}^n \left[ {{a_k}{g_k}\left( {\theta , \varphi } \right)\hat \theta + {b_k}{g_k}\left( {\theta , \varphi } \right)\hat \varphi } \right] $

(3)

式中，ν(θ, φ)为测站速度，g_k(θ, φ)为DOG小波基函数，a_k和b_k为待估参数。由于g_k(θ, φ)是经过平移和伸缩得到的，属于非正交基，存在着法方程系数矩阵秩亏的问题，因此解算出的待估参数不唯一。本文可采用广义交叉验证法(generalized cross-validation，GCV)确定正则化参数，结果如图 3所示。

根据多尺度球面小波法可以模拟郯庐断裂带南段区域的速度场。图 4为多尺度球面小波估计速度场的矢量图，图 4(a)为实际速度场与估计速度场的对比，图 4(b)为残差矢量图。图 5为多尺度球面小波估计速度场的标量图，图 5(a)为实际的速度场，图 5(b)为利用球面小波理论估计的速度场，图 5(c)为二者之间的残差。从图 4和图 5可以看出，利用球面小波法估计得到的速度场与实际速度场的误差在2 mm/a以内，整体运动情况也大致相同。因此，多尺度球面小波法可以有效估计速度场。

2.3 应变场解算

在利用GCV法确定模型的正则化参数后，可利用Savage等^[9]、孟国杰等^[10]提出的在球面坐标系下应变与位移的微分公式，来求得目标点的应变张量。由于垂直方向上产生的形变对于水平方向的应变影响较小，且将不处于同一球面的测站投影到同一球面来进行应变解算的过程较为复杂，本文不予考虑。应变张量的计算表达式为：

$ \begin{array}{l} \;\;\;\;\;\;\;\;\;\;\;\;\;\;D = \left[ {\begin{array}{*{20}{c}} {{e_{\theta \theta }}}&{{e_{\theta \varphi }}}\\ {{e_{\theta \varphi }}}&{{e_{\varphi \varphi }}} \end{array}} \right] = \\ \left[ {\begin{array}{*{20}{c}} {\frac{1}{r}\cdot\frac{{\partial {u_\theta }}}{{\partial \theta }}}&{\frac{1}{{2r}}\left( {\frac{{\partial {u_\varphi }}}{{\partial \theta }} + \frac{1}{{\sin \theta }}\cdot\frac{{\partial {u_\theta }}}{{\partial \theta }} - {u_\varphi }\cot \theta } \right)}\\ {{D_{23}}}&{\frac{1}{r}\left( {\frac{1}{{\sin \theta }}\cdot\frac{{\partial {u_\varphi }}}{{\partial \varphi }} + {u_\theta }\cot \theta } \right)} \end{array}} \right] \end{array} $

(4)

式中，r表示地球半径，$ \frac{{\partial {u_\theta }}}{{\partial \theta }}、\frac{{\partial {u_\theta }}}{{\partial \varphi }}、\frac{{\partial {u_\varphi }}}{{\partial \theta }}、\frac{{\partial {u_\varphi }}}{{\partial \varphi }}$表示速度的梯度张量，θ、φ分别表示余纬和经度。由上文所求得的应变张量可得到主应变率、面膨胀率以及剪应变率。

3 郯庐断裂带南段地壳形变与断层活动分析 3.1 应变场分析

利用多尺度球面小波法计算郯庐断裂带南段的面膨胀率结果如图 6所示，(a)~(d)4幅图分别表示空间分解尺度6~9下的面膨胀率, 图中面膨胀率为正值表示膨胀，负值表示压缩。当q_max=6时，要求最小站点间的距离不超过174.4 km。该图可以明显反映出整个郯庐断裂带南段大区域的构造特点，面膨胀率最大值为9.995×10^－9/a。当q_max=9时，要求最小站点间的距离不超过21.8 km，此时面膨胀率最大值为1.976×10^－8/a。总体上看，郯庐断裂带南段以西除阜阳、徐州、枣庄等地，大部分地区呈面膨胀状态，以东地区大多表现为面压缩状态。比较不同空间尺度下的面膨胀率可以发现，某些局部区域存在差异。因此，在小的分解尺度下可以反映出大尺度的应变积累特征，在大的分解尺度下能够反映出小区域范围内的应变积累特征。

比较最大分解尺度q_max=8和q_max=9时的应变特征结果可以发现，两者间的差异较小，故本文将针对最大分解尺度q_max=8时的应变特征进行分析。

图 7为q_max=8时郯庐断裂带南段及周边地区的主应变率/面膨胀率、最大剪应变率的分布。从图 7(a)可以看出，郯庐断裂带泗洪至莒南段西侧的徐州附近主要表现为E-W向的拉张，商丘北部表现为NE-SW方向的拉张、NW-SE方向的压缩；商丘附近则相反，转变为NE-SW方向的压缩、NW-SE方向的拉张。宿松、岳西等地表现为S-N向压缩，黄山、婺源等地表现为近E-W向的拉张。郯庐断裂带的泗洪至莒南段东侧大多表现为S-N向的压缩。整体上看，郯庐断裂带南段西侧的主应变比东侧的大，其中最大张应变位于郯庐断裂带西侧的徐州附近，主应变的最大值为2.327×10^－8/a。商丘、岳西、婺源等应变方向变化异常地区，均处于面膨胀与面压缩的交界地带，地质构造可能较为复杂，存在着较大的应变积累，后期有发生地震的可能性，应予以关注。从图 7(b)可以看出，最大剪应变在徐州附近的表现最为异常，达到4.044×10^－8/a。其次，商丘、枣庄、阜阳、黄山等地最大剪应变率也较大。总体上看，最大剪应变率沿郯庐断裂带南段两侧呈对称分布，其西北和东南地区为高值区域，存在着应变积累，应予以关注。

3.2 速度剖面分析

为进一步分析郯庐断裂带南段及其周边地壳形变和断层活动情况，本文横跨郯庐断裂带南段设置3个速度剖面，选择的剖面位置如图 1(b)所示，剖面长度为600 km，宽度为100 km。图 8为3个速度剖面的结果，由北向南依次排列，图中圆点表示AHCORS参考站的速度，三角形表示CMONOC站点的速度。表 3为断层滑动速率的估值结果。从图 8和表 3可以看出，郯庐断裂带南段主要表现为右旋走滑，其北部存在着拉张，拉张分量向南逐渐递减，在最南部转变为压缩的状态。平行断裂的运动速率为0.10~1.21 mm/a，垂直断裂的运动速率为-0.52~0.76 mm/a。郭良迁等^[11]利用块体-应变模型反演得到的郯庐断裂带表现为右旋挤压，李彦川等^[2]利用负位错模型得到的郯庐断裂带南部的特征为右旋走滑为主兼拉张。上述结论与表 3的a、b段结果基本吻合。由于本文还结合了AHCORS的速度场，所以b段的拉张分量较弱，更能反映该区域的断裂带特征。c段剖面处于郯庐断裂带桐城-庐江段，宋传中等^[12]通过地质学分析认为，郯庐断裂带桐城-庐江段在新第三纪以来表现为逆冲挤压的特征，与本文得到的c段结论较为一致。

4 结语

本文采用GAMIT/GLOBK解算2013-01~2018-06期间AHCORS参考站在欧亚框架下的速度场，同时结合CMONOC速度场利用多尺度球面小波法计算郯庐断裂带南段及周边地区的应变率，并对横跨郯庐断裂带南段的3个速度剖面进行分析。得出以下结论：

1) 安徽地区在欧亚框架下主要朝东偏南方向运动。本文利用多尺度球面小波法模拟出的速度场与实际速度场的误差在2 mm/a以内，可有效反映该区域的地壳运动情况。

2) 总体上看，郯庐断裂带南段以西除阜阳、徐州、枣庄等地外，大部分地区呈面膨胀状态，以东地区大多表现为面压缩状态。根据最大分解尺度下的主应变率/面膨胀率和最大剪应变率的分布情况可以看出，应变方向变化异常的地区，均出现在面膨胀与面压缩的交界地带。最大剪应变率沿郯庐断裂带南段两侧呈对称分布，其西北和东南地区为高值区域，构造特征复杂，存在着较大的应变积累，有发生地震的可能，应予以关注。

3) 郯庐断裂带南段主要表现为右旋走滑，且北部存在着拉张，拉张分量向南逐渐递减，在最南部转变为压缩的状态。平行断裂的运动速率为0.10~1.21 mm/a^-1，垂直断裂的运动速率为-0.52~0.76 mm/a。