0 引言

据中国地震台网测定，2021年5月21日21时48分云南省大理白族自治州漾濞县(25.67°N，99.87°E)发生M_S6.4地震，震源深度8km。通过对GNSS观测资料数据进行处理，可获取地震前后应变率场，观测、研究地壳运动应变场及其随时间的变化过程，对探索震源的孕育过程和地震的监测预报具有重要意义(徐克科等，2017)，目前国内外诸多学者都积极致力于利用GPS观测资料探索地壳运动与形变的机制，研究地震孕育发生的形变过程，捕捉与认识孕震信息(杨国华等，2008；魏文薪等，2012)。杨少敏等(2003)利用GPS区域观测资料计算分析了2011年昆仑山8.1级地震前中国大陆水平构造应变场空间分布特征；魏文薪等(2018)利用GPS同震位移及反演分析得到鲁甸M_S6.5地震沿NE向表现出拉张应变释放的同震特征；朱治国等(2016)通过GPS观测资料分析认为精河M_S6.6地震前应变率场反映出区域构造动力环境。

“十一五”期间，由中国地震局牵头的国家重大科学工程项目“中国大陆构造环境监测网络工程”(简称陆态网络)启动，项目在网络工程一期的基础上布设了高密度、大范围的GNSS观测，为GNSS观测资料应用于地壳活动性分析、地球动力学研究及地震预测预报奠定了基础。本研究所采用的GNSS数据包括网络工程位于云南及邻区的55个连续观测站和近300个GNSS区域站的多期观测数据。采用GAMIT/GLOBK软件对观测数据进行解算，在获取各个测点时间序列和速度场的基础上，采用克里金插值方法分时段估计该区域2009—2015年和2015—2020年两个时间区域应变率场，了解漾濞地震前区域应变场背景信息；同时对2014年以来云南省内GNSS连续数据进行处理，以连续基准站时间序列为约束，根据李延兴等(2001、2004、2009)提出的块体运动计算模型，获取漾濞M_S6.4地震近场区域的块体时间序列。通过漾濞地震前后多种应变参数的时空变化特征，分别从区域应变场背景和发震区域局部应力动态变化两个角度分析地震前后地壳运动的动态演变过程，以期捕捉与认识孕震信息，提取可作为震情跟踪的预测判据。云南及邻区构造背景及周边GNSS观测站具体分布见图 1。

1 区域应变场背景分析 1.1 数据解算

本文采用MIT(Massachusetts Institute of Technology)和SIO(Scripps Institution of Oceanography)研制的GAMIT/GLOBK软件对数据进行处理。选取BJFS、LHAZ、PIMO、USUD、WUHN等十几个中国及周边的IGS跟踪站与研究区域内GNSS站点联合解算。每24h作为一个分析时段，利用GAMIT获取包含测站位置、大气参数、地球定向参数以及相应的方差-协方差矩阵的单日松弛解；利用GLOBK卡尔曼滤波程序将GAMIT处理得到的单日松弛解和SOPAC(Scripps Orbital and Permanent Array Center)给出的全球IGS站的单日松弛解(h-文件)合并，选取全球范围内均匀分布、时间序列稳定的IGS连续观测站作为框架约束，获取待解算测站在该框架下的点位结果。综合近几年来的单日解，可得到ITRF2014参考框架下点位时间序列，通过对不同时间段时间序列线性拟合获得站点在ITRF2014参考框架下的速度场。

图 2给出了漾濞地震前云南及邻区相对于稳定欧亚板块的GNSS水平速度场，红色速度矢量为2009—2015年速度场结果，蓝色速度矢量为2015—2020年速度场结果，并绘制了两期速度结果75% 置信水平的误差椭圆，后期的应变率场分析均以此为约束。

1.2 区域应变率场计算方法

在获得各期扣除了大震同震影响的速度场的基础上，将速度场转化为区间位移场，再采用克里金插值方法以经纬度0.2°×0.2°为单元对其进行格网化。区域应变场计算采用GLOBK软件的NEU站心坐标系。站心坐标系是一种地方空间直角坐标系，其坐标原点在一个选定的测站O点上，其北向坐标轴(N坐标)为过O点的子午线切线，指北为正；其东向坐标轴(E坐标)为过O点的椭球的平行圈切线，指东为正。在实际计算过程中，通过坐标N、E分量计算网格中心点经纬度信息，并通过大地反算计算任意两个网格之间的大地线长度以及大地线长度的变化量，并计算网格中心点之间的大地方位角。每个格网与其他各个相邻格网点通过线应变与方位角之间的关系联立方程组，再通过最小二乘法求解得到应变状态分量，进一步可以计算其他应变参数，包括最大剪应变、面应变等(Hong et al，2018)。

假设某个测点的位移为u、v，其应变状态分量为ε_x、ε_y、γ_xy，其中，$\varepsilon_{x y}=\frac{1}{2} \gamma_{x y}$，d_x、d_y为两点间距离分量的变化量，ω为旋转量，α为两格网点间的坐标方位角，则与其无限接近的一点的位移分量可表示为

$ \left\{\begin{array}{l} u^{\prime}=u+\varepsilon_x d_x+\varepsilon_{x y} d_y-\omega d_y \\ v^{\prime}=v+\varepsilon_{x y} d_x+\varepsilon_y d_y+\omega d_x \end{array}\right. $

(1)

式(1)等号两边同时除以两点间距离，可转变为线应变与方位角之间的关系

$ \left\{\begin{array}{l} \frac{\Delta u}{d}=\varepsilon_x \cos \alpha+\varepsilon_{x y} \sin \alpha-\omega \sin \alpha \\ \frac{\Delta v}{d}=\varepsilon_{x y} \cos \alpha+\varepsilon_y \sin \alpha+\omega \cos \alpha \end{array}\right. $

(2)

野外观测及实验室结果均表明，由于断层的闭锁或其他震前机制的影响，往往会导致大地震前震源区域与其周围区域应变状态出现差异(马瑾等，2014)。为了解漾濞M_S6.4地震前震源区域与其周围区域应变状态，本文分别给出了漾濞地震前2009—2015年、2015—2020年2期速度场得到的应变率场结果。

1.3 应变特征分析

最大剪应变率是两个最大剪切方向的剪应变率，是应变参量中反映综合变形强度的参量(江在森等，2006)。图 3分别给出了2009—2015年和2015—2020年云南及邻区最大剪应变率场及以当前应变率场为背景的M_S5.0以上地震分布。结果显示，云南省内最大剪应变率高值区位于川滇菱形块体的整个东边界、西边界与丽江—小金河断裂带交汇区域，西南向主要表现为零星散布的腾冲、普洱等高值区域。前人研究发现，地震发生频率与最大剪应变率的大小有一定联系，地震往往发生在最大剪应变率的高值区域。本文所得最大剪应变率场也呈现出类似的特点，即研究区域内M_S5.0以上地震均发生在前期最大剪应变率高值或高低值边界区域，这与江在森等(2006)的研究结果具有一致性。

漾濞地震序列位于维西—乔后断裂的SW侧，整体呈NW—SE向展布，该断裂历史上发生过5.0~6.3级地震，地震主要发生在中段和东南段，其中最大地震为1839年的洱源两次6.3级地震。该区域构造复杂，属于多条活动断裂的交汇部位，红河断裂带北段、程海断裂西南段构成了川滇菱形块体的西南边界，与维西—乔后断裂中南段交汇，研究表明：强活动断裂带上各分支断裂的结合部分，不同方向断裂的交汇部位或断裂本身的拐点、枢纽点、交叉点和端点等闭锁构造部位是易发生中强地震的危险区。该区域曾发生1652年7月13日弥渡南7级和1925年3月16日大理7级地震。从图 3来看，川滇菱块西边界丽江—小金河断裂带交汇区域包括丽江、鹤庆、剑川等地，其一直处于最大剪应变率最高值，该区域也是云南历史地震较频繁的地区，与地震多数是由断层走滑运动而产生的结论相符。2008年汶川8.0级地震后川滇菱块西边界与维西—乔后断裂中南段交汇区出现明显的最大剪应变高梯度异常(图 3(a))，最大值约为12×10^-8/a，说明在该期间剪切活动加强，之后几年内该区陆续发生了2016年云龙M_S5.0和2017年漾濞M_S5.1地震。2015—2020年，该区域剪切应变强度有所减缓(图 3(b))。但相对于云南其他区域，该区域仍然存在应变高梯度状态。前人研究认为，地震孕育过程中断层两侧速度差异会逐渐降低(Meade et al，2005)，表明随着断层附近应变的积累，区域应变积累速率逐渐降低。当应变率低到一定程度时，标志着断层附近应变积累达到极限状态，进一步积累应变，断层将发生破裂。因此，较大地震并非发生在最大剪应变率的极大值区域，而是发生在高低值边界区域。此次漾濞地震也正好发生在前期应变高值区，发震的时间也处于区域应变积累速率逐渐降低的过程之后。

面应变率为两个方向主应变率之和，其反映的是地壳当前处于膨胀或压缩的状态，一定程度上反映逆逆断层或张性正断层的应变积累或释放速率。为研究漾濞地震前断层附近局部区域的应变积累速率的大小，并据此推测断层附近的应变积累程度，笔者给出了滇西北漾濞地震震区附近98°E~102°E、24°N~28°N区域面应变率场及主应变率场结果，图中面压缩为负值、面拉张为正值(图 4)。

研究时段内主压应变率及面应变率存在以下几个特征：汶川地震后滇西北丽江小金河断裂南段永胜、洱源、丽江附近，面应变状态长期以面拉张为主，主张应变率明显大于主压应变率，川滇菱块西边界与维西—乔后断裂中南段交汇区应变强度明显高于周围区域(图 4(a))。块体及其周缘部分变形相对较弱，主应变率方向沿交汇区边界断裂带逐渐发生旋转，相对于维西—乔后断裂，从南到北侧云龙附近逐渐变化为以主压应变率为主的应变积累状态；该区域附近出现明显的高值面应变正负交替异常，之后2年内陆续发生了2016年云龙M_S5.0、2017年漾濞M_S5.1地震，从主应变率方向与断裂走向判断，维西—乔后断裂中南段应变强度明显强于北段，中南段明显处于构造应力的转换带，反映了该地区较为复杂的构造背景。图 4(b)显示该区域仍处于拉张状态，但强度有所减缓，主应变方向明显改变，主压应变率从NW向主张应变转换成近SN向，且主张应变强度减弱，主压应变区域范围扩大。漾濞M_S6.4地震位于面应变张压转换区附近，与笔者前期研究结果吻合(王伶俐等，2020)。与最大剪应变率类似，主应变率与面应变率可以反映震前断层附近应变积累速率的大小，并通过积累速率推测断层附近的应变积累程度。区域长期的构造应力被认为是稳定不变的，而在断层失稳前，由于各种震前机制的影响，短期应力往往会偏离长期的构造应力(马瑾等，2014)。在野外观测中，这种震前偏离所表现出的就是震前主应变方向偏离长期构造主应力方向。漾濞地震序列震前主应变率方向和强度的异常可能源于地震前应变方向偏离长期构造应力方向。

2 发震区域近场应变时序分析 2.1 块体应变计算模型

根据李延兴等(2001、2004、2009)提出的块体运动计算模型REHSM，以GNSS坐标时间序列为约束，将模型转换为间接平差函数模型，按最小二乘原理，采用求自由极值的方法解出参数的最或然值，进而求得各观测值的平差值(王伶俐等，2016)。计算得到各活动块体的运动参数和应变参数为

$ \left[\begin{array}{l} v_{\mathrm{E}} \\ v_{\mathrm{N}} \end{array}\right]=\left[\begin{array}{ccc} -r \cos \lambda \sin \phi & -r \sin \lambda \sin \kappa & r \cos \phi \\ r \sin \lambda & -r \cos \lambda & 0 \end{array}\right]\left[\begin{array}{c} \omega_x \\ \omega_y \\ \omega_z \end{array}\right]+\left[\begin{array}{cc} \varepsilon_{\mathrm{E}} & \varepsilon_{\mathrm{EN}} \\ \varepsilon_{\mathrm{NE}} & \varepsilon_{\mathrm{N}} \end{array}\right]\left[\begin{array}{l} \left(\lambda-\lambda_0\right) r \cos \phi \\ \left(\phi-\phi_0\right) r \end{array}\right] $

(3)

通常，为了更清楚方便地了解块体内部应变的情况，需进一步计算应变率的最大主压应变率及其方位角等参数，计算公式为

$ \left\{\begin{array}{l} \varepsilon_1=\frac{1}{2}\left(\varepsilon_{\mathrm{E}}+\varepsilon_{\mathrm{N}}\right)-\frac{1}{2}\left[4 \varepsilon_{\mathrm{EN}}^2+\left(\varepsilon_{\mathrm{E}}-\varepsilon_{\mathrm{N}}\right)^2\right]^{\frac{1}{2}} \\ \varepsilon_2=\frac{1}{2}\left(\varepsilon_{\mathrm{E}}+\varepsilon_{\mathrm{N}}\right)+\frac{1}{2}\left[4 \varepsilon_{\mathrm{EN}}^2+\left(\varepsilon_{\mathrm{E}}-\varepsilon_{\mathrm{N}}\right)^2\right]^{\frac{1}{2}} \\ A=\arctan \left(\frac{\varepsilon_{\mathrm{EN}}}{\varepsilon_1-\varepsilon_{\mathrm{E}}}\right) \\ \mathit\gamma_{\max }=\varepsilon_2-\varepsilon_1 \\ \mathit\Delta=\varepsilon_1+\varepsilon_2 \end{array}\right. $

(4)

其中，ε₁和ε₂分别为最大主压应变率和最大主张应变率；A为最大主压应变率的方位角；γ_max为最大剪应变率；Δ为面应变率。

2.2 应变参数时间序列分析

为提高解算精度，在前期获得的包含全球IGS站和云南GNSS基准站的单日松弛解的基础上，对每4天的数据综合平差求出一个基于ITRF2014框架的点位静态解。综合近几年来各期解算结果，可得到ITRF2014参考框架下点位时间序列，图 5为部分站点的原始时间序列。

从原始时间序列来看，GNSS连续运行基准站相对全球参考框架的时间序列长期构造运动趋势显著，年速率估计的均方根误差为4~6mm，原始时间序列虽然在时间上有利于识别测站的动态变化，但在空间上并不利于观察研究区内部运动的差异；由于应变参数可以从整体上刻画研究区域的变形特征，且与参考基准无关，针对研究区域相对较小、基准站点密度不大的实际情况，更适合采用块体应变来分析区域变形特征。与GNSS时间序列相似，应变时序可以反映研究区域的变形信息，且与参考基准无关，同时应变结果还可以反映不同方向的变形特征。

为分析地震前后区域构造形变面应变的动态演变特征，本文选择震区附近9个GNSS基准站(包括云南丽江(YNLJ)、云南宾川(YNBC)、云南大姚(YNYA)、云南云龙(YNYL)、云南下关(XIAG)、云南楚雄(YNCX)、云南腾冲(YNTC)、云南施甸(YNSD)、云南景东(YNJD))构成(a)、(b)、(c)、(d)四个形变单元(图 6)。

以GNSS基准站点位时间序列为约束，参照李延兴等(2001、2004、2009)提出的块体整体旋转与均匀应变模型(REHSM)，求解发震区域的应变参数时间序列，包括最大剪应变、面应变、第一剪应变和第二剪应变4个物理量，研究漾濞地震前后多种应变参数的时空变化特征(图 7)。

第一、二剪应变率可以反映走滑断层的应变积累速率，第一剪应变率反映NE向或NW向的剪切变形，第二剪应变率反映EW向与SN向的剪切变形。从第一剪应变时间序列来看，震区周围4个区域均以NW向断层的右旋走滑应变积累为主，且大多呈现持续增强趋势。震区以东形变单元图 7(c)区域相对其他区域应变积累较弱，值得注意的是，该区域震前1年左右，即2020年9月开始出现趋势性异常变化，第一剪应变从持续性增强转为阶段性减弱，时间持续近5个月后又恢复增强趋势。震前1个月，即2021年4月又再次出现了趋势性异常变化，这种状态一直持续到地震发生，经过震后的调整又恢复到原来的线性增强趋势。

为进一步确定其趋势异常，笔者对该区域应变序列进行线性滤波分析，以4期数据为窗长进行滑动均值计算，滤去部分高频成分，再对时序曲线去趋势处理，结果如图 8所示。分析发现该区域第一剪应变在2020年9月—2021年2月超过一倍标准差幅度，虽然其他时段也出现过类似异常，但此次异常幅度与维持的时间均更为显著；由于该区域此前5级以上地震样本较少，虽不能确认其变化是否与地震孕育相关，但从第一剪应变率代表的物理意义来看，这一异常表示该区域NW向断裂右旋走滑应变积累的减弱。

从第二剪应变来看，震区西部形变单元图 7(a)持续负值且呈现持续增强趋势，反映了该区域SN向断层的右旋应变积累。震区南部形变单元图 7(d)持续正值，反映了该区域SN向断层的左旋应变积累特征，其他两个区域应变积累特征并不明显。第二剪应变在临震前后并未表现明显的趋势性异常变化。

最大剪应变是两个最大剪切方向的剪应变，是应变参量中反映综合变形强度的参量(江在森等，2006)。测区西部形变单元图 7(a)和北部形变单元图 7(b)区域最大剪应变强度最大，与前期基于区域网解算得到的漾濞M_S6.4地震前震源区域与其周围其他区域应变率结果一致。究其原因，可能是该区域维西—乔后断裂、澜沧江断裂、近SN向丽江—大理断裂系统等多条断裂带发育，区域内部GNSS站点相对位移大，使得解算所得应变值较高。震区南部形变单元图 7(d)和东部形变单元图 7(c)的最大剪应变强度相对较弱，说明块体内部较为稳定，变形速率较低，相对于其他两个区域剪切运动较弱。从整体来看，4个区域最大剪应变时间序列大多呈现持续增强趋势，与第一、第二剪应变一致，震区南部图 7(d)区域也表现出从趋势转折到转折加速、闭锁状态、发震、震后调整和恢复到原变化趋势的过程。

值得一提的是，震区附近4个区域第一剪应变和最大剪应变显现出明显同震响应，此次漾濞M_S6.4地震恰好发生在NW向红河断裂带向西北延伸段——维西—乔后断裂带(又称通甸—巍山断裂)附近。主震与前震的震源机制一致性较高，震源机制显示为走滑运动为主，这与第一剪应变反映NE向、NW向剪切应变相符，表明第一剪应变更能够提取与NW向的发震断裂关联紧密的地壳形变信息。

整体来看，4个区域面应变时间变化序列特征呈现明显分化，震区西部形变单元图 7(a)面膨胀持续增强，该区域2016年以前为负值，以面压缩为主，2016年后以面拉张为主，与前期基于区域网解算得到的面应变率演变过程一致。震区南部形变单元图 7(d)及北部形变单元图 7(b)区域以弱压缩为主，应变状态较稳定，东部形变单元图 7(c)面膨胀持续增强。震前该区域出现了短期应变趋势转折现象，这种反向加速的异常现象，或许反映了发震区挤压应变积累速率明显加大，应力-应变积累在接近临界破裂状态时的非线性调整，加速了此次漾濞地震的孕育过程。

3 结论和讨论

通过对漾濞M_S6.4地震震前区域网速度场解算获得的应变率场研究发现：川滇菱块西边界丽江—小金河断裂带交汇区域一直处于最大剪应变率最高值，面应变状态长期以面拉张为主，主张应变率明显大于主压应变率，川滇菱块西边界与维西—乔后断裂中南段交汇区应变强度明显高于周围区域。块体及其周缘部分变形相对较弱，主应变率方向沿交汇区边界断裂带逐渐发生旋转，相对于发震断层(维西—乔后断裂)，从南到北侧云龙附近逐渐变化为以主压应变率为主的应变积累状态；此次漾濞地震正好发生在前期最大剪应变高值区以及面应变高梯度带的张压转换区，发震的时间也处于区域应变积累速率逐渐降低的过程之后，这与笔者前期的研究结果相吻合(王伶俐等，2020)。

震区近场区块的应变参数时间序列反映了区域地壳形变的动态演变过程，结果显示：震区周围区域均以NW向断层的右旋走滑应变积累为主，且大多呈现持续增强趋势。最大剪应变和第一剪应变异常变化相对最为明显，漾濞地震发震断层位于维西—乔后断裂中南段与近SN向丽江—大理断裂的交汇部位，为走向135°且陡倾向WS的右旋走滑断层，这与第一剪应变反映NE和NW向剪切应变是相符的。震区东部块体震前出现了短期应变趋势转折现象，这种反向加速的异常现象，反映了发震区挤压应变积累速率明显加大，应力-应变积累在接近临界破裂状态时的非线性调整，加速了此次漾濞地震的孕育过程，与笔者前期研究结果一致(王伶俐等，2020)。

漾濞地震周围4个块体震前应变趋势并不统一，并非所有块体在震前都有明显异常。这除了与震级和地震破裂方向有关之外，可能还与块体所选站点离震中或者发震构造的距离有关，在块体应变时，尽管所用模型是应变模型，但由于测站之间的空间跨距较大，其中包含了许多断裂带，所估计的应变参数并非均匀应变结果，而是一个块体作为一个整体的应变参数，其所反映的是一个块体的整体平均应变状态，但仍不影响对整体趋势性的判断。由于GNSS观测资料积累时间尚短，加之近场GNSS基准站较少，虽然基准站点位时间序列所反映的相对运动与应变积累的动态变化在一定程度上能够客观反映地震孕育发展过程中的变化异常，但仍需继续探索产生这种异常变化现象的原因，进一步揭示地震孕育发生的动力学背景。