基于聚类分析的川滇地区块体划分与应变研究

张小青　许才军　王帅

摘要 本文研究川滇地区块体划分与区域地壳应变.首先，对川滇地区GPS水平速度场采用K-medoids聚类法直接聚类;然后，基于改进的K-medoids方法，利用GPS站点地理位置与水平速度场对初步聚类结果联合约束进行第2次聚类分析，并结合F统计检验法确定各块体边界;最后，运用整体旋转线性应变模型计算各块体的应变参数，分析水平运动、应变场空间分布与区域构造变形的关系.研究结果表明：川滇地区的GPS聚类分析结果与地质学提出的板块划分结果相吻合，为块体划分提供另一种方式;整体旋转线性应变模型计算的应变结果显示川滇菱形块体表现为最大剪应变和面应变率梯度高值区，表明该区域积累了一定的弹性应变能，具有发生地震的危险性.

关键词 GPS速度场;川滇地区;聚类分析;应变

中图分类号 P315

文献标志码 A

0 引言

川滇地区位于青藏高原东缘，为印度板块向青藏高原东北向挤压、碰撞作用的前沿地带，区内构造活动显著和强震频发，是中国大陆构造变形最复杂和最强烈的区域之一.近年来，随着GPS技术的进一步发展以及GPS观测网络的密集布设，利用GPS观测数据研究块体运动和地壳变形，已经成为地壳形变研究的热点[1-5].利用GPS 观测到的地壳形变速度场，结合地质、地球物理学等观测资料，通过构建地表变形的场源数学物理模型，可进一步求解地壳应力-应变参数，进而深入掌握区域地壳形变的内在机制特征[1-5].党亚民等[1]根据两期GNSS实测速度场，运用平移-旋转-应变模型，研究了川滇地区地震前后块体运动与应变特征.温扬茂等[2]联合GPS资料、重力资料反演了川滇地区的现今地壳形变，结果表明川滇地区地壳运动模式更趋向于连续变形.丁开华等[3]根据川滇地区1998—2004 年间的GPS观测结果，基于地质学给出的块体划分模型，反演了川滇地区各活动块体的现今地壳运动和变形参数，认为川滇地区的地壳变形模式以连续和渐变的剪切为主.依据刘峡等[4]、党亚民等[1]和廖思佩等[6]对川滇地区活动块体划分的研究结果，川滇地区活动块体可分为巴颜喀拉块体（Ⅰ）、川滇菱形块体（Ⅱ）、滇东块体（Ⅲ）、藏东块体（Ⅳ）和滇西南块体（Ⅴ）.其中巴颜喀拉块体（Ⅰ）以龙日坝断裂为界又分为阿坝（Ⅰ1）和龙门山（Ⅰ2）两个次级块体;川滇菱形块体（Ⅱ）以丽江—小金河断裂为界可进一步划分为川西北次级块体（Ⅱ1）和滇中（Ⅱ2）两个次级块体，各板块分布如图1所示.

合理有效地划分构造活动块体是构造运动及地震活动性研究的重要基础.许多研究者致力于板块，尤其是板内块体及构造亚板块的划分[3-10].石耀霖等[7]对基于GPS 观测资料计算的欧拉极进行聚类分析，将中国大陆划分为 9 个活动块体，该方法在计算欧拉极时存在研究范围大、计算精度低的不足，且研究区域内GPS测站较少，分布不均匀，板块划分难于精细.Savage等[8]对加利福尼亚东南部莫哈韦地区的GPS水平速度进行聚类分析，将莫哈韦地区分成4个聚类地块，分类结果与Meade[9]的划分结果类似.

本文首先利用K-medoids聚类算法对川滇地区的GPS水平速度场进行聚类分析，获得初始聚类中心;然后，联合GPS站点地理位置与速度场数据对初始聚类结果进行约束，重新判断各点与初始聚类中心的距离，进行第2次聚类，并运用F统计检验方法进一步筛选GPS测站，得到最终的聚类分析结果;最后，基于块体划分结果，利用整体旋转线性应变模型计算川滇地区地壳运动与变形参数.

1 聚类算法与数据

1.1 K-medoids聚类算法

K中心点算法（K-medoids）的基本思想是：对含有n个数据的数据集，欲分k个簇，首先任意选择k个不同的数据作为初始簇中心，再根据每个数据对象到簇中心的距离，将它们重新分配到距离最近的簇;然后通过反复地用非中心点代替初始簇中心，而选取的标准就是当该样本点成为新的中心點后能提高类簇的聚类质量，使得类簇更紧凑.K-medoids相比k-means算法采用均值计算法获得中心点，能改善k-means的“噪声”敏感特性[11].

1.2 川滇地区GPS水平运动速度场

本文利用Liang等[12]、Gan 等[13]公布的川滇地区1999—2013年间GPS观测数据.数据处理主要有以下4步：

1）利用GAMIT软件处理GPS数据，获得站坐标和卫星轨道的单日松弛解;

2）将获得的单日松弛解与SOPAC提供的全球解用GLOBK软件进行融合;

3）用QOCA软件根据卡尔曼滤波方法计算全球参考框架ITRF2000站坐标和速度;

4）将基于全球参考框架ITRF2000的测站速度转换到稳定欧亚参考框架的速度.

川滇地区的水平运动速度场由北向南、西向东逐渐减弱，呈顺时针涡旋式运动[1-3]（图1）.

2 聚类结果

2.1 K-medoids直接聚类

首先，利用K-medoids聚类分析方法对川滇地区的水平运动速度场进行初步聚类分析.在聚类分析过程中，设定K取不同的值，以顾及不同聚类数对聚类分析结果的影响.图2a—2d分别是聚类数K取3、4、5和6时的聚类分析结果，其中同一集群的GPS站点用相同的符号表示.

图2a为聚类数K取3时的分析结果，可以看出川滇菱形块体与滇东块体能够很好地被划分出来;但巴颜喀拉块体中的龙门山块体与滇东块体被划分到了一起，其中滇东块体与滇西南块体也被分在一起.值得说明是，龙门山块体中的一个点被错误地划分到川滇菱形块体中.

图2b为聚类数分成4类的分析结果，可以看出与3类聚类分析结果相似，川滇菱形块体与滇东块体能够被很好地识别并划分出来.除此之外，滇东块体与滇西南块体也被成功地划分出来，巴颜喀拉块体中的阿坝块体与藏东块体被分在一起，并且龙门山断裂带未被检测到，玉树-甘孜断裂带也未被检测到.其中滇东块体的3个点被错误地划分到滇西南块体中（图2b粉色三角形表示的点），同时巴颜喀拉块体和滇西块体各有1个点和2个点被错误地分到滇东块体中（蓝色方块点所示）.

分5類的结果如图2c所示，与分4类聚类分析结果相比，巴颜喀拉块体中的阿坝块体与滇西块体已经被成功分开，并且能够检测到玉树-甘孜断裂带，但龙门山块体与滇东块体被分在了一起.虽然能够识别出滇西块体与川滇菱形块体，但滇东块体中有一些点被错误地划分到滇西南块体中.

相比前几类分析结果，分6类的聚类分析能够很好地识别出川滇地区的各个块体，如图2d，川滇菱形块体可分为川西北和滇中两个次级块体.

2.2 K-medoids聚类结果分析

因K-medoids算法聚类中心选择的随意性，即使两次聚类分析取相同的聚类数，聚类分析的结果也存在较大的差异，如图3a—3b所示.所以，需要一定的指标来评定聚类分析的结果的可靠性.目前聚类分析效果的评定指标主要有轮廓系数、戴维森堡丁指数和偏差估计等[14].本文采用轮廓系数对聚类分析结果进行评定.将待分类数据分成k个簇，对簇中的各点分别计算轮廓系数S（i）：

其中，a（i）为样本i到同簇其他样本的平均距离，b（i）为样本i到其他簇的所有样本的平均距离.根据式（1），簇内距离越小，簇间距离越大，也即Sk值越大，聚类效果越好.图3中，聚类数为5的图3a和3b中，图3a聚类轮廓系数小于3b的轮廓系数值，也就是说图3b聚类效果更好;聚类数是6的图3c和3d中，图3c轮廓系数小于图3d，显然图3d聚类效果更好.

2.3 改进的K-medoids聚类分析

由于对GPS速度场直接聚类时未考虑站点的地理空间分布连贯性，分类结果中仍有一些点被错分.为获得更好的聚类结果，本文对K-medoids聚类算法进行改进.首先利用K-medoids进行初步聚类，获得初始聚类中心;然后联合各点地理位置与速度场数据对初始聚类中心进行约束，重新判断各点与原始聚类中心的距离，进行第2次聚类.利用改进后的K-medoids算法获得的聚类结果如图4所示，可以看出聚类效果得到了明显的改善，不再有错分的点.对比图1与图4b，可看出改进的聚类算法获得的结果与基于地质学提出的块体划分结果基本吻合.

2.4 GPS 速度解的筛选

值得注意的是，虽然改进的K-medoids算法能够明显改善聚类的结果，但仍有部分测站的速度值与周边测站明显不符，出现这种速度奇异点的原因很复杂，可能是由观测期间测站的不稳定性等原因造成的.在板块构造划分和地壳形变研究中，常使用F统计检验法来判定点的相符性[15]，从而研究块体间的独立性.采用同时考虑块体刚性运动和内部应变的整体旋转模型来计算各块体的应变参数[16]，进而根据块体模型和应变参数预测、估计各测站速度，拟合速度值见图1中红色箭头，黑色箭头为GPS观测速度.针对块体边界点采用F统计检验，若该测站与其他测站相互独立，则将此站点从该块体去除，划到地理相邻的块体中.F统计量：

式中n为所选用的测站数，χ2n为所有测站速度共同运动的欧拉极拟合后的速度残差，

χ2n-1

为去除某一测站拟合后的速度残差.

3 块体运动特征分析

3.1 块体应变分析

通过上述聚类算法与测站筛选法确定各块体边界之后，采用同时考虑块体刚性运动和内部应变的整体旋转模型来计算各块体的欧拉矢量（表1）与应变参数（图5），藏东、巴颜喀拉、滇西南、滇中和滇东块体的旋转速率与肖卓辉等[17]的对应区域大致相同.

由于2008年汶川地震后川滇地区震后形变机制复杂，本文使用的1999—2013年GPS速度场未考虑2008以后川滇地区的形变，因此本文只对2008年汶川地震之前的应变状态进行分析.图5a为川滇地区最大剪应变结果，结果表明红河断裂带、安宁河断裂带、小江断裂带和丽江—小金河断裂带处的最大剪应变比较显著，尤其在鲜水河断裂带北段以及玉树—甘孜断裂带最为显著，与江在森[18]的研究结果相吻合，表明川滇地区受到较明显的外力.

从图5b可以看出，川滇地区内部以面膨胀为主，而沿川滇块体边界以面压缩为主，并且在红河断裂带、安宁河断裂带、小江断裂带和丽江—小金河断裂带处的面应变比较显著，与江在森[18]以及魏文薪等[19]的研究结果相吻合.

图5c是面应变梯度结果，可以看出川滇地区西部的面应变梯度比东部大，同样表现在红河断裂带、安宁河断裂带、小江断裂带和丽江—小金河断裂带，面应变梯度最大可达4.8×10-10 km.龙门山断裂带面应变梯度为中值区，应变特征不明显.滇西南地区的面应变梯度高值区转换为梯度中值区，这是由于震后地壳的松弛和变形调整所致.面应变梯度高低值交汇处主要沿区内主干走滑断裂带分布，1900年以来发生的中强以上地震多出现在面应变梯度高低值交汇处.

3.2 块体间相对运动

在获取块体欧拉运动矢量后，本文以GPS点的外包络线作为块体边界，计算了边界点的相对运动速度，计算结果如表2所示.从表2可看出，巴颜喀拉块体与藏东块体、滇东块体以挤压为主，巴颜喀拉块体与川西北块体以扩张为主，川西北块体与藏东块体、滇中块体以挤压为主，滇中块体与滇西南块体以挤压为主，滇中与滇东块体以扩张为主，这与党亚民等[1]的研究结果基本吻合.

4 结论

本文基于1999—2004年以及1999—2008年间川滇地区的GPS观测资料，在K-medoids初步聚类的基础上，联合各点地理位置与速度场数据对初始

聚类中心进行约束，进行第2次聚类，能有效改善被错分的站点，最终将川滇地区划分为6个次级块体，包括巴颜喀拉块体、川滇菱形块体、滇东块体、藏东块体和滇西南块体.其中川滇菱形块体以丽江—小金河断裂为界可进一步划分为川西北和滇中两个次级块体.块体划分结果与地质上块体划分模型相吻合，表明聚类分析用于块体划分的合理性，可与地质上块体划分方法形成互补.值得一提的是，利用GPS观测资料进行聚类分析，块体划分精度依赖于GPS站点的密集度.由于各站点速度代表邻近区域的平均速度，因此，GPS站点越密集，其所代表的邻近区域越小，块体划分精度越高.

在聚類分析的基础上，采用F统计检验法能有效去除部分与周边测站的速度值明显不符的测站，估算得到的各块体欧拉矢量与其他研究者得到的结果相吻合.给出的应变特征表明：红河断裂带、安宁河断裂带、小江断裂带和丽江—小金河断裂带处的最大剪应变比较显著，尤其在鲜水河断裂带处最为显著，川滇地区西部是最大剪应变作用最广泛、最强烈的区域.川滇地区以面膨胀为主，整个川滇菱形块体为应变能和面应变率梯度高值区，为不同活动方式的构造断裂共同作用的结果.

参考文献

References

[1] 党亚民，杨强，梁诗明，等.川滇区域活动块体运动与应变特征地震影响分析[J].测绘学报，2018，47（5）：559-566

DANG Yamin，YANG Qiang，LIANG Shiming，et al.Block movement and strain characteristics effected by earthquake in Sichuan-Yunnan region[J].Acta Geodaetica et Cartographica Sinica，2018，47（5）：559-566

[2] 温扬茂，许才军.联合GPS与重力资料反演分析川滇地区现今地壳形变[J].武汉大学学报（信息科学版），2009，34（5）：568-572

WEN Yangmao，XU Caijun.Contemporary crustal deformation in Sichuan-Yunnan region from GPS and gravity data[J].Geomatics and Information Science of Wuhan University，2009，34（5）：568-572

[3] 丁开华，许才军，邹蓉，等.利用GPS分析川滇地区活动地块运动与应变模型[J].武汉大学学报（信息科学版），2013，38（7）：822-827

DING Kaihua，XU Caijun，ZOU Rong，et al.Crustal movement and strain model of active blocks analyzed by GPS in Sichuan-Yunnan region[J].Geomatics and Information Science of Wuhan University，2013，38（7）：822-827

[4] 刘峡，马瑾，杜雪松，等.川滇主要断裂带近期运动变化及与地震活动关联性[J].中国科学（地球科学），2016，59（6）：1267-1282

LIU Xia，MA Jin，DU Xuesong，et al.Recent movement changes of main fault zones in the Sichuan-Yunnan region and their relevance to seismic activity[J].Science China Earth Sciences，2016，59（6）：1267-1282

[5] 金双根，朱文耀.太平洋板块运动和形变及其边缘现今相对运动[J].大地测量与地球动力学，2002，22（2）：57-60

JIN Shuanggen，ZHU Wenyao.The motion and deformation of pacific plate and contemporary relative motion at its boundaries[J].Crustal Deformation and Earthquake，2002，22（2）：57-60

[6] 廖思佩，侯强，杜永超.基于GPS形变资料的川滇地区应力场数值模拟研究[J].大地测量与地球动力学，2016，36（7）：645-649

LIAO Sipei，HOU Qiang，DU Yongchao.Numerical simulation of the stress field in Sichuan-Yunnan region with GPS deformation data[J].Journal of Geodesy and Geodynamics，2016，36（7）：645-649

[7] 石耀霖，朱守彪.利用GPS观测资料划分现今地壳活动块体的方法[J].大地测量与地球动力学，2004，24（2）：1-5

SHI Yaolin，ZHU Shoubiao.Method for division of present active crustal blocks by gps survey data[J].Crustal Deformation and Earthquake，2004，24（2）：1-5

[8] Savage J C，Simpson R W.Clustering of GPS velocities in the mojave block，southeastern california[J].Journal of Geophysical Research（Solid Earth），2013，118（4）：1747-1759

[9] Meade B J.Block models of crustal motion in southern California constrained by GPS measurements[J].Journal of Geophysical Research，2005，110（B3）：B03403

[10] Aranganayagi S，Thangavel K.Clustering categorical data using silhouette coefficient as a relocating measure[C]∥International Conference on Conference on Computational Intelligence & Multimedia Applications.IEEE，2007

[11] 朱纯，吴建华，潘毅.关于优化K-medoids聚类算法搜索策略研究[J].计算机仿真，2016，33（10）：244-248，277

ZHU Chun，WU Jianhua，PAN Yi.Research about optimizing search strategy of K-medoids clustering algorithm[J].Computer Simulation，2016，33（10）：244-248，277

[12] Liang S M，Gan W J，Shen C Z，et al.Three-dimensional velocity field of present-day crustal motion of the Tibetan Plateau derived from GPS measurements[J].Journal of Geophysical Research（Solid Earth），2013，118（10）：5722-5732

[13] Gan W J，Zhang P Z，Shen Z K，et al.Present-day crustal motion within the Tibetan Plateau inferred from GPS measurements[J].Journal of Geophysical Research，2007，112（B8）：B08416

[14] 贾虎，邓力珲.基于流线聚类人工智能方法的水驱油藏流场识别[J].石油勘探与开发，2018，45（2）：312-319

JIA Hu，DENG Lihui.Oil reservoir water flooding flowing area identification based on the method of streamline clustering artificial intelligence[J].Petroleum Exploration and Development，2018，45（2）：312-319

[15] 郭東美，许厚泽.基于统计检验法研究Amurian板块相对独立性[J].地球物理学报，2013，56（4）：1246-1252

GUO Dongmei，HSU Houtse.Research on the existence of Amurian plate based on statistical tests[J].Chinese Journal of Geophysics，2013，56（4）：1246-1252

[16] 李延兴，张静华，何建坤，等.由空间大地测量得到的太平洋板块现今构造运动与板内形变应变场[J].地球物理学报，2007，50（2）：437-447

LI Yanxing，ZHANG Jinghua，HE Jiankun，et al.Current-day tectonic motion and intraplate deformation-strain field obtained from space geodesy in the Pacific Plate[J].Chinese Journal of Geophysics，2007，50（2）：437-447

[17] 肖卓辉，许才军，江国焰，等.汶川地震前十年间龙门山区域顾及断层闭锁的地壳应变场[J].地球物理学报，2017，60（3）：953-961

XIAO Zhuohui，XU Caijun，JIANG Guoyan，et al.Crustal strain in the Longmenshan region considering fault locking during ten years before the 2008 Wenchuan earthquake[J].Chinese Journal of Geophysics，2017，60（3）：953-961

[18] 江在森.利用动态大地测量资料研究中国大陆构造形变与强震关系[D].武汉：武汉大学，2012

JIANG Zaisen.Research on the relationship between the tectonic deformation and the strong earthquake in the Chinese mainland using dynamic geodetic data[D].Wuhan：Wuhan University，2012

[19] 魏文薪，江在森，刘晓霞，等.川滇地区应变率场分布及其变化特征研究[J].地震，2015，35（4）：11-20

WEI Wenxin，JIANG Zaisen，LIU Xiaoxia，et al.Distribution and variation characteristics of strain rate field in the Sichuan：Yunnan region[J].Earthquake，2015，35（4）：11-20

Block partitioning and strain study in Sichuan-Yunnan

region based on cluster analysis

ZHANG Xiaoqing1 XU Caijun2 WANG Shuai2

1 Department of Hydraulic Engineering，Fujian College of Water Conservancy and Electric Power，Yongan 366000

2 School of Geodesy and Geomatics，Wuhan University，Wuhan 430079

AbstractThis paper focuses on block partitioning and strain calculation in the Sichuan-Yunnan region.First，the K-medoids clustering method is used to cluster the GPS horizontal velocity field in Sichuan-Yunnan.Next，the initial classification results are further clustered based on constraints from the geographical location and velocity field of GPS stations by using the improved K-medoids method.Subsequently，the F statistical test method is used to determine the boundary of each block.Finally，the strain parameters of each block are calculated using an integral rotation linear strain model and the relationship between horizontal motion，spatial distribution of strain field，and regional tectonic deformation is analyzed.The results of the velocity field clustering analysis of GPS stations in the Sichuan-Yunnan region are consistent with the results of the alternative block division methodology based on geological theory of plate division.The strain results，which are calculated using the integral rotation linear strain model，show that the Sichuan-Yunnan rhombic block has the maximum shear strain and the highestsurface strain rate gradient.In addition，this region has accumulated some elastic strain energy，an indicator of earthquake risk.

Key wordsGPS velocity field;the Sichuan-Yunnan region;cluster analysis;strain