地球物理学报  2020, Vol. 63 Issue (7): 2516-2533   PDF    
引用本文
任营营, 王虎, 王解先, 等. 2020. 基于K-Means++的省内子块体划分及中国大陆水平相对运动速度场模型的建立与分析. 地球物理学报, 63(7): 2516-2533, doi: 10.6038/cjg2020N0401.  
Ren Y Y, Wang H, Wang J X, et al. 2020. The sub-block demarcation with K-Means++ in each province's interior and establishment analysis of the relative horizontal velocity field model in Mainland China. Chinese J. Geophys. (in Chinese), 63(7): 2516-2533, doi: 10.6038/cjg2020N0401.  
基于K-Means++的省内子块体划分及中国大陆水平相对运动速度场模型的建立与分析
任营营1, 王虎2,3, 王解先1, 连丽珍4, 朱卫东5, 王永哲6, 张守建3     
1. 同济大学测绘与地理信息学院, 上海 200092;
2. 中国测绘科学研究院, 北京 100830;
3. 武汉大学测绘学院, 武汉 4300079;
4. 中国科学院上海天文台, 上海 200030;
5. 上海海洋大学, 上海 201306;
6. 中国地震局地球物理研究所, 北京 100081
收稿日期 2019-11-11, 2020-04-08 收修定稿
基金项目: 国家自然科学基金面上（41874042），中国博士后科学基金（2016M590715）资助
第一作者简介: 任营营, 男, 硕士研究生, 研究方向为大地测量学、GNSS算法与数据处理.E-mail:ryy@tongji.edu.cn
通讯作者: 王虎, 男, 博士, 副研究员, 主要从事GNSS数据处理及地学应用的研究.E-mail:wanghu@casm.ac.cn
摘要：省市级区域CORS（Continuously Operating Reference Stations）系统是当代城市数字化、信息化和智慧化的重要组成部分，便于获取各类物体的时空信息及其相关动态变化.为了进一步实现区域框架基准的现代化与自主化，全面提升现代测绘基准综合服务水平和应急保障能力，同时为了提高中国大陆区域水平速度场的精度，并精细地刻画其自身的局部运动特征，本文利用陆态网上千站2011—2017年的连续观测数据，采用GAMIT/GLOBK软件，获得高精度的定位和速度成果，进而提出和构建了基于欧拉矢量模型的中国大陆省级块体相对运动模型和部分省内子块体相对运动模型，并与欧亚板块、大陆整体和二级板块相对运动模型进行对比分析；结果表明，欧亚板块相对运动模型仅能描述大陆的部分运动趋势，中国大陆整体板块相对运动模型能够较好地展现大陆整体运动趋势，二级板块和省级块体相对运动模型均能够较为精细地反映区域的局部运动特征，且两者水平相对速度的内符合精度均小于2 mm·a-1，其外符合精度均小于3 mm·a-1，其中前者物理意义更为明显，后者使用更为简便，但在青藏、川滇等地壳运动复杂的地区两者精度仍有欠缺.因此，本文提出利用K-Means++算法对地壳运动复杂区域的水平速度场进行聚类分析，以快速准确地对这些区域进行子块体划分；结果表明，划分成果与现有部分二级块体成果相符合.为了兼顾省内复杂地质构造与地形地貌的影响，同时提高省级块体划分的物理意义，对地壳运动复杂的省份再细分块体，进而对各子块体构建欧拉矢量模型；结果发现，该模型平均误差和中误差均小于2 mm·a-1，提高新疆、西藏、川滇等地区的速度场模型精度至2 mm·a-1左右，在确保精确度的同时，满足使用简便性的要求.
关键词: 陆态网      速度场      欧拉矢量      块体划分      相对运动     
The sub-block demarcation with K-Means++ in each province's interior and establishment analysis of the relative horizontal velocity field model in Mainland China
REN YingYing1, WANG Hu2,3, WANG JieXian1, LIAN LiZhen4, ZHU WeiDong5, WANG YongZhe6, ZHANG ShouJian3     
1. College of Surveying and Geo-informatics, Tongji University, Shanghai 200092, China;
2. Chinese Academy of Surveying and Mapping, Beijing 100830, China;
3. School of Geodesy and Geomatics, Wuhan University, Wuhan 430079, China;
4. Shanghai Astronomical Observatory, Chinese Academy of Science, Shanghai 200030, China;
5. Shanghai Ocean University, Shanghai 201306, China;
6. Institute of Geophysics, China Earthquake Administration, Beijing 100081, China
Abstract: Provincial CORS system is an important part of modern city digitization, informatization and intellectualization, which is easy to obtain the spatiotemporal information of various objects and their related dynamic changes. In order to further realize the modernization and autonomy of the regional framework benchmark, improve the comprehensive service level and emergency support capacity of the modern surveying benchmark, upgrade the accuracy of the regional velocity field and accurately depict its own local motion characteristics in the mainland of China, the velocity field is derived from continuous observation data of the Crustal Movement Observation Network of China (CMONOC) that was surveyed in 2011—2017. Both high precise coordinate and velocity of these CMONOC stations have been calculated with GAMIT/GLOBK software. In addition, the Euler vector of the whole Chinese continent and two-stage plates in the China mainland are reanalyzed based on the above velocity field. On the basis of the Euler vectors of the above-mentioned, i.e. the whole Chinese continent, two-stage plates, provincial plates, sub-blocks in some provinces' interior, as well as that of Eurasia plate derived from NNR-NUVEL1A plate motion model, the contribution of the average movement of each block to the horizontal velocity field could be available and then the residual velocity fields in Mainland China are further analyzed and compared. The results show that the NNR-NUVEL1A model only deducts the partial movement trend of the velocity field in Mainland China. In contrast, the overall movement model of Mainland China plate can better reflect the overall movement trend. The movement model of two-stage plates and provincial plates in Chinese mainland, both of whose inner and outer average precision is less than 2 mm·a-1 and 3 mm·a-1, respectively, can both depict the local movement characteristics of Mainland China more precisely, and the only difference between them is that the former is easier to understand in the sense of physics while the latter is simpler in using. However, both of them is not perfect for researching the relative motion of some regions such as Xinjiang, Tibetan Plateau, Sichuan, Yunnan where exist complex crustal movement. Hence, through cluster analysis with K-Means++ on horizontal velocity field in these complex areas, these provincial blocks are further decomposed into some sub-blocks quickly and accurately. Our results are in accordance with those of present-day blocks of second order. Consideration to the impact of complex geological structure, topography and geomorphology in these provinces and the physical significance in respect to each corresponding provincial plate, we also give the Euler vectors of these sub-blocks in some provinces' interior. The result can be concluded that not only the average error and mean square error are less than 2 mm·a-1 in Mainland China but also the precision of horizontal velocity fields can reach 2 mm·a-1 in some complex regions, that is Xinjiang, Tibet, Sichuan, Yunnan. Obviously, the scheme of the sub-block demarcation in some provinces' interior by the method of K-Means++ is effective and feasible in terms of high precision and convenience.
Keywords: CMONOC    Velocity field    Euler vector    Block demarcation    Relative motion    
0 引言

随着毫米级现代化大地测量技术的不断发展，高精度速度场成为大地测量工作中不可或缺的基础性要素之一，有助于解释和反映中国大陆整体或局部区域的内部构造的变化特征.由于测站坐标受地壳运动和形变等各种因素的影响会随时间而发生变化，因此必须通过高精度速度场模型推算测站在用户所需历元时刻的高精度点位坐标才具有工程实用价值，那么如何建立简便的高精度速度场模型就显得尤为重要.

通常，采用ITRF框架下的速度场对于描述大区域的运动特征具有一定的便利和优势，但对于较小的区域，必须要扣除相应的背景速度场才能更真实地反映其自身内部的运动特征.欧拉矢量模型是将某一地表区域当作刚性块体，可以很好地描述该块体的整体运动特征.为了构建中国大陆高精度速度场模型，近年来，诸多学者利用不同时期的观测数据建立了不同的基于欧拉矢量的速度场模型.将欧亚板块当作一个块体，通过构建欧亚板块运动模型来描述中国大陆的运动特征，这可以扣除大陆的部分运动趋势(Argus and Gordon, 1991; Zhao et al., 2015)；将中国大陆整体作为一个块体，可以分离整体运动趋势，但这无法描述局部区域的运动特征(Qu et al., 2018; Yu et al., 2018; 魏子卿等, 2011)；将大陆分为若干不同的块体，虽然可以得到较高精度的产品，但这需要软件支持(姚宜斌等, 2002; 刘经南等, 2002; 谢方等, 2015)；通常划分块体越小，其精度越高(游新兆等, 2012; 杨元喜等, 2011; 苗岳旺等, 2019)，由于中国大陆内部没有明确的具体块体划分，因此利用局部移动窗口模型计算的速度场精度最高，但其较为复杂且实用性不强(吴富梅等, 2012).

由于可用的数据较少、数据连续性不佳、数据现势性不强以及缺乏明确的应用背景等各种因素，使得现有研究成果的实用性价值有一定的局限性.此外，仅仅以省级边界作为块体划分界线(简称省级块体)，不具有确切的物理意义，而且基于此建立的速度场模型产品在地壳运动较为复杂地区的精度欠佳.因此，本文利用陆态网的上千站2011—2017年观测数据解算得到的高精度速度场，建立了基于欧拉矢量的中国大陆欧亚板块相对运动模型、整体板块相对运动模型和二级板块相对运动模型，提出了省级块体相对运动模型，且在省级块体基础上，利用K-Means++方法划分部分省内部的子块体，同时建立省内子块体运动模型，并对比分析了这5种模型.其中，结合中国大陆省内地壳运动特征再细分子块体的省级块体运动模型具有较强的科学研究意义和工程实用价值，进一步提升了省级区域基准框架的自主维护和应用水平.

1 中国大陆高精度水平速度场的获取

中国大陆构造环境监测网络(Crustal Movement Observation Network of China，CMONOC，简称陆态网)主要由均匀分布在全国的260余连续运行基准站和2000余不定期观测区域站构成.其中，基准站观测大多起始于2010—2011年，其数据总体时间跨度约7—8年；区域站观测大多始于1999年，在2008年以前，进行不定期复测，在2009年以后，陆态网二期流动站进一步加密，新增1078个流动站，复测频率定为2年1次.

为了将国内GNSS观测站同IGS站进行联合解算，本文选取国内及周边的17个IGS站进行基准约束.图 1中三角形表示基准站，圆点表示区域站，五角星表示IGS框架站.

图 1 陆态网及周边IGS站点分布图 Fig. 1 Distribution of CMONOC and its surrounding IGS stations

本文采用GAMIT/GLOBK软件分别处理陆态网基准站和区域站数据，首先按区域将基准站和区域站进行子网分区，并将17个IGS站加入每个分区数据进行单天基线解算(基线解算配置见表 1)，得到单日松弛解，然后将全球IGS站的各分区的单日松弛解进行合并，合并时仅保留稳定的IGS核心站，再将合并好的全球IGS单日解与陆态网各分区的单日松弛解进行合并，得到包含全球IGS站和陆态网测站的坐标、卫星轨道和极移等参数的单日整体松弛解，最后基于所有的单日整体松弛解文件提供的法方程进行时序平差和综合处理，结合相似变换法解算陆态网测站的坐标和速度等信息.在求解测站坐标和速度过程中，通过加入eq重命名文件，已将受地震运动和仪器更换等影响而发生阶跃变化的测站进行了分段处理(测站分段名称如XIAM_GPS、XIAM_2PS、XIAM_3PS…)，其中仪器更换等因素仅影响某一时刻的测站位移发生变化，但测站分段前后的速度保持一致，若测站周围发生地震，测站会受到地震同震位移和震后形变带来的影响，分段前后测站的速度可能发生变化.因此在进行时序分析求取测站速度时，剔除粗差后加入了地震同震位移改正，考虑了震后形变对坐标序列的影响，经分段处理，综合求解最终得到陆态网测站高精度水平速度解.

表 1 基线解算方案设置 Table 1 Baseline solution settings

为了保证数据解算结果的质量，解算过程中剔除观测时间跨度小于3年和坐标时间序列残差较大的测站.图 2所示为ITRF2014框架下陆态网基准站和区域站的速度场，IGS框架站和陆态网基准站水平速度精度优于1 mm·a-1，陆态网区域站水平速度精度优于2 mm·a-1.

图 2 ITRF2014框架下陆态网测站水平速度场 Fig. 2 Velocity field of CMONOC stations in ITRF2014
2 欧拉矢量模型的建立与评估 2.1 欧拉矢量的基本平差原理

板块的运动通常采用欧拉矢量来描述.如果把地球近似为一个规则球体，地球质心保持不变，地球表面的块体当作一个刚体，那么，根据欧拉定理可将地球表面块体的运动视为过地心的定轴转动(任雅奇, 2012)，均满足公式(1)：

(1)

式中，V为块体内某站点的速度，Ω为该站点所处板块的欧拉矢量，R为该站点的位置矢量.可改写为便于计算的公式(任雅奇, 2012)：

(2)

式中，Ve、Vn为某测站的水平方向速度；R为地球平均半径，约为6371393 m；λ、φ为该测站的大地经纬度，单位rad；Ωx、Ωy、Ωz为欧拉矢量的三个分量.

根据间接平差原理，设欧拉矢量的参数为向量X，站点的水平速度观测向量L，其观测误差向量为Δ，可组成观测方程

(3)

设向量X的估值为，则有误差方程

(4)

上述误差方程的最小二乘解为

(5)

为欧拉矢量三参数的协因数阵，则参数的中误差估值为

(6)

式中，为的主对角元素，为单位权中误差估值，有

(7)

式中，r为多余观测量，V为观测值改正数，与Δ大小相等，符号相反，可由计算得到的欧拉参数估值代入误差方程得到.Δ由东、北方向的速度误差向量Δve和Δvn组成，则东、北方向的速度中误差估值为

(8)

式中，n为观测站个数.

2.2 欧拉极的数值转换模型

笛卡尔坐标系中的Ωx、Ωy、Ωz在球坐标中的分量表示为λ、φ、ω，其中λ、φ为板块旋转的经纬度，ω为旋转角速度，其相互对应关系为

(9)

(10)

式中，λ取值范围为[-90, 90]、φ取值范围为[-180, 180].

3 中国大陆水平相对运动速度场模型的建立与评估 3.1 基于欧拉矢量的欧亚板块相对运动模型

地球表面是由若干不断运动和变化的板块组成的，地质学家依据近百万年的地球物理观测资料推导出各个板块运动的平均速度.目前国际推荐使用的是NNR-NUVEL1A板块运动模型，该模型将全球划分为16个板块，中国大陆位于欧亚板块之中.表 2给出NNR-NUVEL1A模型欧亚板块的欧拉矢量(Argus and Gordon, 1991).

表 2 NNR-NUVEL1A模型欧亚板块的欧拉矢量(单位：rad/Ma) Table 2 Euler vector of Eurasian plate based on NNR-NUVEL1A model

将陆态网各区域站的水平速度值和经纬度以及表 2中的欧拉矢量参数值代入公式(2)，可得到NNR-NUVEL1A模型下各区域站点的运动速度，代入公式(4)进而得到相对于欧亚板块的中国大陆区域站剩余速度场，如图 3所示.

图 3 相对欧亚板块的中国大陆水平运动速度场 Fig. 3 Horizontal velocity field of Chinese mainland relative to Eurasian plate

根据计算得到的区域站的剩余速度场(即相对于欧亚框架背景速度场的运动场)，代入公式(8)计算得到NNR-NUVEL1A模型的中误差.表 3给出东、北方向相对速度场的绝对值的平均值(Mve和Mvn)以及剩余速度的中误差(Sve和Svn).

表 3 欧亚板块相对运动模型精度统计结果(单位：mm·a-1) Table 3 Accuracy of relative motion model of Eurasian plate

分析图 3可得：相对于欧亚板块运动，中国大陆内部具有西部地区整体向东北方向、东部地区整体向东方向的运动趋势，且剩余速度场较大，显然西部地区的速度场精度低于东部地区.结合表 3分析可得：NNR-NUVEL1A模型的水平方向的中误差接近于8 mm·a-1，其整体误差较大.显然，NNR-NUVEL1A模型仅能消除中国大陆速度场的部分运动趋势，无法表现中国大陆整体的运动趋势.

NNR-NUVEL1A模型是利用地球物理观测数据推导出来的，反映的是地球在几百万年内的平均板块运动，由其所推算的欧亚板块的整体运动不能代表其近期运动情况，因此必须建立适合中国大陆板块整体欧拉运动模型，以便能够更好地展现其自身内部的运动趋势.

3.2 基于欧拉矢量的大陆块体相对运动模型

将中国大陆整体当作一个刚性块体，使用陆态网2000余个区域站的实测速度，根据欧拉矢量模型的平差原理，求取中国大陆的整体欧拉矢量.为避免粗差对解算的影响，需进行迭代计算，每次解算后以三倍中误差准则判断数据集中是否含有粗差并剔除，直至不再含有粗差结束迭代，最终得到中国大陆整体板块运动的欧拉矢量.利用解算得到的欧拉矢量参数，将其代入公式(4)即可得到相对中国大陆整体板块的区域站速度场，如图 4所示.

图 4 相对中国大陆整体板块自身的水平运动速度场 Fig. 4 Horizontal velocity field of Chinese mainland relative to itself

为了验证本文计算的中国大陆整体板块相对运动的欧拉矢量模型(记为CHINA-2017)的精度和可靠性，特与国内学者魏子卿等(2011)利用“地壳网络工程”1999—2008年的观测数据解算的速度场模型(记为CHINA-2008)进行比较.表 4、5给出两种模型的欧拉参数和速度残差的统计结果.

表 4 两种大陆整体欧拉矢量参数对比(单位：rad/Ma) Table 4 Comparison of Euler vector parameters of two models in Chinese mainland
表 5 两种大陆整体欧拉矢量模型的精度统计(单位：mm·a-1) Table 5 Accuracy of two Euler vector models in the Chinese mainland

分析表 4、5可得：两种模型的欧拉参数值基本一致，欧拉三参数的最大差值小于0.001 rad/Ma；水平方向平均误差和中误差大致相同，但CHINA-2017模型可以达到5 mm·a-1的精度，略优于CHINA-2008模型的，出现这种差异的主要原因是两者采用的数据不一致，前者是依据1999—2017年的观测数据，而后者主要依据1999—2008年的观测数据，且随着近年来陆态网的不断建设与发展，陆态网测站的分布更均匀、数量更多、数据质量更佳，因此CHINA-2017较CHINA-2008更具现势性，且精度更加可靠，适用时间范围更广.

分析图 4可得：中国大陆相对于自身整体运动具有明显的区域运动特征，其中运动最为剧烈的区域主要集中在青藏活动地块边界，新疆区域整体向西北方向运动，西藏南部由南向北运动，青藏边界区域和川滇地区由西向东呈90°顺时针旋转，华南华北平原地区的活动效应相对较弱，华南地区主要由西向东南方向呈180°逆时针微小旋转，华北地区主要由东向西南方向呈180°方向微小旋转，东北地区整体由东南向西北方向逆时针旋转.综合分析可知，中国大陆内部的运动呈现三大态势：东部地区以环渤海区域为中心呈现逆时针旋涡状的运动趋势；西南部地区以青藏地区为中心呈顺时针旋转趋势；西北部地区以准噶尔和天山的分界线为边界，分别向西南、西北方向运动.

相较于NNR-NUVEL1A模型，中国大陆块体相对运动模型能更好地描述大陆内部的运动特征.中国大陆东西部的相对运动速度差别较大，西部地区的运动速度明显大于东部地区；其中，华南华北地区的相对运动速度最小.将中国大陆整体作为刚性块体求取的欧拉矢量，仅可以较好地描述华南华北等地区的水平运动特征，但在东北地区、西部大部分地区仍有较大的误差，因此中国大陆整体欧拉矢量模型仍无法精细描述各区域的运动特征.

3.3 基于欧拉矢量的二级板块相对运动模型

中国大陆整体板块欧拉矢量法和局部欧拉矢量法的主要区别在于参与计算的站点所占的区域大小.前者是依据中国大陆整体区域进行计算，而后者是依据现有大陆20个二级板块(杨元喜等, 2011)进行计算，每一个块体都有一组对应的欧拉矢量.中国大陆二级板块划分成果见图 5.

图 5 中国大陆二级板块划分示意图 Fig. 5 Secondary block demarcation in Chinese mainland

将中国大陆任一二级板块当作一个刚性块体，首先采用射线法判断陆态网各区域站点分别位于哪一个块体之中，并分类提取每一块体的区域站信息，然后使用对应块体的区域站点，根据欧拉矢量模型的平差原理迭代计算，最终得到中国大陆的各二级块体的欧拉矢量参数.利用解算得到的各块体欧拉矢量估值，可得到相对于二级块体的中国大陆区域站速度场，如图 6所示.表 6给出了中国大陆各二级板块欧拉矢量参数和水平方向速度场的精度.

图 6 相对二级板块的中国大陆水平运动速度场 Fig. 6 Horizontal velocity field of Chinese mainland relative to the secondary plate
表 6 中国大陆二级板块的欧拉参数和水平速度场精度 Table 6 Euler parameters and the accuracy of the horizontal velocity field of the secondary plate in Mainland China

为了评估中国大陆二级块体相对运动模型(记为CHINA-PLATE)的区域站速度场整体转换精度，表 7统计了中国大陆二级板块相对运动模型的整体区域站速度场精度.

表 7 二级板块相对运动模型的精度(单位：mm·a-1) Table 7 Accuracy of the relative motion model of the secondary plate

由图 6分析可得：采用中国大陆二级块体相对运动模型可以更好地描述各块体的运动趋势，除青藏活动地块边界处与华北-东北活动地块边缘处个别站存在较大残差外，其余各地区的残差均较小.块体内部活动最为剧烈的是川滇和青藏地区，川滇地区整体由西向北逆时针旋转，滇西南地区呈逆时针旋涡状运动趋势，青藏地区由南向西、向北、向东呈发散状.

结合表 6和表 7分析可得：中国大陆二级板块相对运动模型的区域站速度场残差整体中误差小于2 mm·a-1，各板块内部区域站残差中误差大部分不超过3 mm·a-1，其中速度残差超过2 mm·a-1的块体主要集中在西部地区，分别为：拉萨、羌塘、巴颜喀拉、川滇和滇西南块体，其中羌塘块体的速度场残差最大，主要是由于所在区域测站数量太少，其余四个块体速度场残差较大是由于自身内部构造运动较剧烈所造成的，与分析图 6得出的结论相吻合.

相较于中国大陆整体板块相对运动模型，二级板块相对运动模型可以更好地反映大陆各块体区域的水平运动，且整体精度更加可靠.但是一般用户较难确定点位所属的二级板块，尤其在板块边界地区，因此必须有相应的软件支持，工程实用性不高.

3.4 基于欧拉矢量的省级块体相对运动模型

在考虑中国大陆二级块体的局部运动特征的同时，为提高模型的工程实用性，把中国大陆省(直辖市、自治区)划分为省级块体，除港澳台地区，各省级块体计算一组对应的欧拉矢量用以描述其局部运动特征.省级块体欧拉矢量模型的求解方法与二级块体模型求解方法一致.首先提取每一省级块体的区域测站信息，然后利用对应省级块体的区域测站信息，根据欧拉矢量模型的平差原理进行迭代计算，最终得到中国大陆的各省级块体的欧拉矢量参数.利用解算得到的欧拉矢量参数，可得相对于省级块体的中国大陆各区域站点水平速度场，如图 7所示.表 8给出中国大陆各省级块体欧拉矢量参数和水平速度精度.

图 7 相对省级块体的中国大陆水平运动速度场 Fig. 7 Horizontal velocity field of Chinese mainland relative to the provincial block
表 8 中国省级块体欧拉参数和精度 Table 8 Euler parameters and its accuracy of Chinese provincial blocks

为了评估中国大陆省级块体相对运动模型(记为CHINA-PRO)的区域站速度场整体转换精度，表 9统计了中国大陆省级块体运动模型的区域站点水平速度场的精度.

表 9 中国省级块体相对运动模型的精度(单位：mm·a-1) Table 9 Accuracy of relative motion model in the Chinese provincial block

由图 7分析可得：采用中国大陆省级块体相对运动模型可以明显减弱局部地区的运动趋势，以东经105°为分界线，大陆西部地区的构造运动更为剧烈，东部地区除个别测站剩余速度较大，其他测站的剩余速度都较小，西部地区以新疆、西藏、青海、四川和云南5个省份的区域站剩余速度较大，其中新疆地区由东北向东南逆时针旋转，西藏地区由南向西、向北、向东呈发散状，青海甘肃边界地区整体由东向西运动，四川北部地区由西向南顺时针旋转，四川南部由东向北逆时针旋转，云南西部由北向西南顺时针旋转，云南、贵州和广西边缘地区由西南向东北方向运动.

由表 8、9分析可得：中国大陆省级块体相对运动模型的剩余速度场整体中误差在2 mm·a-1左右，各块体内部区域站剩余速度中误差大部分不超过2 mm·a-1，其中剩余速度超过2 mm·a-1的省份主要集中在四川、云南、西藏、新疆四个省份，四个省份块体相对于自身块体速度较大是由于内部构造运动较剧烈且变形不一致性造成，与图 7得出的结论一致.

相较于中国大陆二级板块相对运动模型，虽然省级块体相对运动模型的整体精度略次，但各省份块体的区域精度更优，且用户可以很方便地确定测站的所属省份，因此在保证精度的同时，其应用更方便.

3.5 基于欧拉矢量的省内子块体相对运动模型

虽然省级块体欧拉矢量运动模型使用方便，且满足大部分省市块体的高精度速度场转换要求，但是由于川滇、青藏地区地壳内部活动较为复杂(如图 8所示)，导致该区域的欧拉矢量模型精度较低.为了提高部分地壳运动比较活跃的省份的速度场模型精度，有必要对省级块体的二次划分作进一步研究，以提高川滇、青藏等地区的速度场模型精度.

图 8 新疆、青藏和川滇地区省级速度场 Fig. 8 Velocity field of provincial blocks in Xinjiang, Tibetan Plateau, Sichuan and Yunnan
3.5.1 基于K-Means++聚类分析的省内子块体划分

对于区域较小的GNSS测站网，HAC法或K-Means聚类法均可有效地实现基于GPS速度场的块体划分(Özdemir et al., 2019).其中K-Means算法属于非监督分类模型，使用广泛且简单有效，但是由于该算法的初始点选取具有随机性，进而可能导致错误的分类结果，因此本文采用改进后的K-Means++(Arthur and Vassilvitskii, 2007)算法.

对于包含N个测站点的省级区域速度场Vi(i=1, 2, …, N)，使用K-Means++进行块体划分的算法可以描述如下：

(1) 随机选取1个速度矢量作为初始聚类中心V1；首先计算每个站点与当前已有聚类中心之间的最短二次型距离dij——两速度矢量差值平方和的开方值，接着计算每个站点被选为下一个聚类中心的概率；最后按照轮盘法选择下一个聚类中心；重复上一步直至选择K个聚类中心Vk(k=1, 2, …, K).

(11)

式中，vie和vin是第i个测站的水平方向速度矢量，vje和vjn是第j个测站的水平方向速度矢量.

(2) 以K个聚类中心为参考点，逐一计算其余站点速度矢量与参考点的距离dij并将其作为两者相似度的考量指标，将该站点的速度矢量分配到距离最小的参考点Vk所在的数据簇中.

(3) 针对每一个数据簇，重新计算其聚类中心.重复步骤(2)和步骤(3)直至聚类中心位置不再发生变化为止.

关于数据簇总数K的选取问题，当K值较小可利用直观有效的肘部法则(Elbow Method)进行选取.即选取不同的K值，计算相应聚类结果的代价函数，绘制聚类个数与代价函数的折现关系图，选取折线拐点幅度最大的K值作为最佳值.K-Means算法最小化问题本质是最小化每个数据点与其所关联的聚类中心点之间的距离之和，所以K-Means算法的代价函数表示为WK(Özdemir et al., 2019)，运用代价函数便于找到更好的簇数，且避免局部最优解的现象发生.

(12)

3.5.2 川滇、青藏地区省内子块体划分结果

对于新疆、西藏、四川和云南地区，分别以聚类中心K=1, 2, …, 10进行K-Means++聚类，绘制每个地区聚类个数与代价函数的关系图，如图 9所示.分析可得，新疆、西藏和四川地区的最佳块体分类个数建议为2，而云南地区地壳运动最为复杂，最佳分类个数建议为4.由于四川与西藏地区测站数量分布稀疏，该地区最佳块体分类个数可能会随后期测站数量的增多而发生改变，本文暂依据现有测站速度场进行块体划分.

图 9 部分省级水平速度场聚类个数与代价函数关系图 Fig. 9 The relationship between the number of clusters and the cost function using velocity fields in some provinces

基于新疆、西藏、四川和云南地区的GPS水平速度场空间聚类分析的块体划分结果如图 10和图 11所示.分析可知，速度场聚类分析识别结果与空间位置信息和地壳构造基本一致：新疆地区聚类分析识别的块体边界与二级块体中准噶尔和天山块体之间的边界基本重合，南北两块体聚类中心的速度场差约为(-2.93，13.88) mm·a-1，两组聚类间的东西向速度场差异较小，南北向速度场差异显著；西藏地区东西两块体聚类中心的速度场差约为(5.77，2.33) mm·a-1，两组聚类间的东西向速度场差异较大，南北向速度场差异较小；四川地区南北两块体聚类中心的速度场差约为(1.95，-7.36) mm·a-1，两组聚类间的东西向速度场差异较小，南北向速度场差异显著；西南地区共划分为四个块体，分别为YN1、YN2、YN3和YN4，各块体的聚类中心分别为：(25.37，-10.45)、(28.91，-14.89)、(33.03，-10.64)和(33.63，-17.54)mm·a-1，其中YN2和YN3两块体聚类中心的速度场差约为(8.25，-7.08)mm·a-1，两组聚类间的东西南北向速度场差异均较大.

图 10 部分省级水平速度场聚类结果 Fig. 10 Clustering results using horizontal velocity field of some provinces′ interior
图 11 部分省内部子块体划分结果 Fig. 11 Sub-blocks in some provinces′ interior
3.5.3 川滇、青藏地区省内子块体的欧拉矢量模型分析

根据上节划分的部分省级内部的块体结果，求取各省内部的子级块体对应的欧拉矢量以精细描述其局部运动特征，再利用解算得到的欧拉矢量参数计算对应块体的剩余速度场，其他省份按照省级块体进行求解，绘制相对部分省内部的子块体的中国大陆水平运动速度场，如图 12所示.表 10给出部分省级子块体欧拉矢量模型的参数和精度.

图 12 相对省内子块体的中国大陆水平运动速度场 Fig. 12 Horizontal velocity field of the Chinese mainland relative to some sub-blocks in provinces′ interior
表 10 部分省级子块体欧拉参数和精度 Table 10 Euler parameters and its precision of sub-blocks in some provinces′ interior

为评估中国大陆省内子块体相对运动模型(记为CHINA-PRO1)的区域站速度场整体转换精度，表 11给出了中国大陆的省内子块体相对运动模型的区域站水平速度场精度.

表 11 部分省内子块体欧拉矢量模型的精度(单位：mm·a-1) Table 11 Accuracy of the Euler vector model of sub-blocks in some provinces′ interior

由图 12分析可得：采用中国大陆省内子块体相对运动模型可以很好地描述青藏和川滇等地壳运动复杂地区的运动趋势，中国大陆剩余速度场基本小于2 mm·a-1.由表 10、11分析可得：中国大陆省内子块体运动模型的速度场残差中误差在2mm左右，其中YN3块体的速度场残差中误差略大，说明扣除省内子块体整体运动后该块体内部还是存在一定的形变.

相较于中国省级板块相对运动模型，省内子块体相对运动模型提高了青藏、川滇等地壳运动较为剧烈的地区的速度场模型精度，同时保证了模型的简便易用性.

3.6 五种速度场模型的外部精度检核

为了检核和评估本文计算的5种欧拉矢量模型(NNR-NUVEL1A、CHINA-2017、CHINA-PLATE、CHINA-PRO和CHINA-PRO1)的精度和可靠性，以长期连续观测的260余陆态网基准站作为外部检核点(建立速度场模型时仅使用2000余区域站，未包含260余基准站)，将陆态网基准站点的速度和经纬度以及5种模型的欧拉矢量参数代入公式(2)，进一步得到不同模型下的基准站的水平相对运动速度.以检核点的经度作为横坐标，水平相对速度场作为纵坐标，绘制基准站的水平相对速度值的对比图，如图 13所示.此外，表 12给出了5种模型的外部精度对比结果.

图 13 不同模型下基准站的水平相对速度对比图 Fig. 13 Comparison of horizontal relative speed of the base stations using different models
表 12 5种欧拉矢量模型的外符合精度比较(单位：mm·a-1) Table 12 Comparison of external accuracy of five Euler vector models

分析图 13、表 12可得：

(1) 以东经105°为分界线，5种模型中的西部地区水平相对速度场明显大于东部地区的，均说明中国大陆西部地区内部构造运动比较剧烈.

(2) 5种模型中NNR-NUVEL1A模型精度最低，其水平速度残差绝对值最大值超过20 mm·a-1，平均误差和中误差在7~9 mm·a-1；CHINA-2017模型精度次之，其水平速度残差最大值小于20 mm·a-1，平均误差和中误差在4~6 mm·a-1；CHINA-2017和CHINA-PRO模型的精度大致相当，其水平速度残差绝对值最大值小于10 mm·a-1，平均误差和中误差在2 mm·a-1左右；CHINA-PRO1模型精度最高，平均误差和中误差均小于2 mm·a-1.

(3) 相较于中国大陆二级块体的划分，省级块体的划分不具有明确的地壳构造依据，因此把省级区域当作刚性块体得到的CHINA-PRO模型精度略次于CHINA-PLATE模型的精度，其中CHINA-PRO模型东、北的中误差为2.51 mm·a-1、2.10 mm·a-1，CHINA-PLATE模型东、北方向的中误差为2.40 mm·a-1、1.72 mm·a-1，但是CHINA-PRO可以很方便地确定测站的所在省级块体，因此CHINA- PRO模型的实用性更强，同时满足高精度速度场转换的需求，CHINA-PRO1模型相较于CHINA-PRO模型物理意义更强，其在地壳运动剧烈区域精度最优.

4 结论

本文基于陆态网2011—2017年的观测数据解算得到的高精度水平速度场，使用区域站构建中国大陆水平相对运动速度场模型并利用基准站速度成果予以评估，分别建立了基于欧拉矢量的中国大陆欧亚板块相对运动模型、整体板块相对运动模型、二级板块相对运动模型、省级块体相对运动模型和省内子块体相对运动模型，分析5种模型的中国大陆相对运动特征，并统计其相对速度场内外符合精度指标，得到以下结论：

中国大陆西部地区内部构造运动较东部更为剧烈，其中以川滇、青藏地区为最.扣除相应模型的背景速度场后，中国大陆内部也表现出不尽相同的区域运动特征：以欧亚板块为准，中国大陆具有西部地区整体向东北方向、东部地区整体向东方向的运动趋势；以中国大陆整体板块为准，中国大陆东部地区以环渤海区域为中心呈现逆时针旋涡状的运动趋势、西南部地区以青藏地区为中心呈顺时针旋转趋势、西北部地区以准噶尔和天山的分界线为边界，分别向西南、西北方向运动的三大运动态势；以中国大陆局部块体为准，青藏地区由南向西、向北、向东呈发散状，川滇块体整体由西向北逆时针旋转的运动趋势.

NNR-NUVEL1A模型仅能消除中国大陆速度场的部分运动趋势，无法表现中国大陆整体的运动趋势，该模型的水平速度残余的外符合平均误差和中误差接近于10 mm·a-1，整体误差较大，精度较低.CHINA-2017模型能够更好地描述中国大陆整体运动特征，该模型的水平速度残余的外符合平均误差和中误差在6 mm·a-1左右，相对于魏子卿等建立的模型(CHINA-2008)，CHINA-2017模型参与计算的观测站分布更均匀、数量更多、数据质量更佳、观测持续时间更长，因此其精度更高且结果也更可靠，虽然使用方便，用一组欧拉矢量即可得到中国大陆区域内任意一点的速度值，但速度场转换精度无法满足工程高精度需求.CHINA-PLATE和CHINA-PRO模型分别以二级板块和省级块体为刚体求解局部块体的欧拉矢量，两者精度大致相当，其水平速度残余的外符合平均误差和中误差均小于3 mm·a-1，其中前者精度略高但使用不便，后者同时保证了高精度性和实用性要求.CHINA-PRO1模型基于K-Means++算法着重为地壳运动较为复杂的川滇、青藏地区省份划分子块体，以便更加准确地描述其局部运动特征，并提高中国大陆区域速度场模型的整体精度至2 mm·a-1左右.

总的来说，本文建立了中国大陆不同区域等级(包括欧亚块体、大陆块体、二级块体、省级块体、省内子块体)的5种欧拉矢量模型，提出和分析了基于省级的和部分省级内部子块体的速度场运动模型，即利用K-Means++算法实现了部分地壳运动较为复杂地区的子块体划分，并建立了相应的局部欧拉矢量模型，以提高省级块体速度场模型的高精度性与简便性，有利于区域基准框架的自身维护，从而进一步提升省级区域CORS站网的现代化与自主化水平.

本文采用参数估计法直接求取测站速度时未能充分考虑坐标序列的震后变形和非线性变化因素，如进行时序分析求取速度来进一步提高测站速度的精度与可靠性.本文在进行块体划分过程中是将站点的水平速度作为其基本属性，未来的研究中如能考虑测站的垂向速度和相关的地壳构造等信息特征，将会进一步提高子块体划分的准确性，基于此考虑，下一步研究中将基于K-Means++进行中国大陆块体划分，详细给出每个块体的欧拉参数及块体号，与现有的大陆块体划分成果进行对比分析，提高块体划分的科学意义和工程应用价值.

致谢  感谢“中国大陆构造环境监测网络”提供部分数据支撑.
References
Argus D F, Gordon R G. 1991. No-Net-Rotation model of current plate velocities incorporating plate motion model NUVEL-1. Geophysical Research Letters, 18(11): 2039-2042.
Arthur D, Vassilvitskii S. 2007. K-Means++: the advantages of careful seeding.//Proceedings of the Eighteenth Annual ACM-SIAM Symposium on Discrete Algorithms. New Orleans, Louisiana, USA: ACM.
Liu J N, Yao Y B, Shi C. 2002. Method for establishing the speed field model of crustal movement in China. Geomatics and Information Science of Wuhan University (in Chinese), 27(4): 331-336.
Miao Y W, Li J, Zhu W, et al. 2019. The method of local velocity field modeling in Mainland China based on Euler vector. Engineering of Surveying and Mapping (in Chinese), 28(1): 26-31.
Özdemir S, Karslıoğlu M O. 2019. Soft clustering of GPS velocities from a homogeneous permanent network in Turkey. Journal of Geodesy, 93(8): 1171-1195.
Qu W, Lu Z, Zhang Q, et al. 2018. Crustal deformation and strain fields of the Weihe Basin and surrounding area of central China based on GPS observations and kinematic models. Journal of Geodynamics, 120: 1-10.
Ren Y Q. 2012. Establishment of Chinese crustal movement velocity model based on the GPS data[Master's thesis] (in Chinese). Zhengzhou: PLA Information Engineering University.
Wei Z Q, Liu G M, Wu F M. 2011. China geodetic coordinate system 2000:velocity field in mainland China. Acta Geodaetica et Cartographica Sinica (in Chinese), 40(4): 403-410.
Wu F M, Liu G M, Wei Z Q. 2012. Velocity field model of CGCS2000 based on Euler vector of local Area. Geomatics and Information Science of Wuhan University (in Chinese), 37(4): 432-435.
Xie F, Cheng C L, Wang B, et al. 2015. Research on the modeling approach of the velocity field in mainland China. Journal of Geodesy and Geodynamics (in Chinese), 35(2): 258-260, 272.
Yang Y X, Zeng A M, Wu F M. 2011. Horizontal crustal movement in China fitted by adaptive collocation with embedded Euler vector. Science China Earth Sciences, 54(12): 1822-1829.
Yao Y B, Liu J N, Shi C, et al. 2002. Construction and application of Chinese crustal motion background field based on ITRF97 frame. Geomatics and Information Science of Wuhan University (in Chinese), 27(4): 363-366, 376.
You X Z, Yang S M, Wang Q. 2012. Progress of crustal movement observation and research in China. Earthquake (in Chinese), 32(2): 31-39.
Yu J S, Tan K, Zhang C H, et al. 2018. Present-day crustal movement of the Chinese mainland based on Global Navigation Satellite System data from 1998 to 2018. Advances in Space Research, 63(2): 840-856.
Zhao B, Huang Y, Zhang C H, et al. 2015. Crustal deformation on the Chinese mainland during 1998-2014 based on GPS data. Geodesy and Geodynamics, 6(1): 7-15.
刘经南, 姚宜斌, 施闯. 2002. 中国地壳运动整体速度场模型的建立方法研究. 武汉大学学报·信息科学版, 27(4): 331-336.
苗岳旺, 李军, 朱璇, 等. 2019. 基于欧拉矢量的中国大陆局部速度场建模方法. 测绘工程, 28(1): 26-31.
任雅奇. 2012.基于GPS数据的中国地壳运动速度场模型的建立[硕士论文].郑州: 解放军信息工程大学. http://cdmd.cnki.com.cn/Article/CDMD-90005-1013161053.htm
魏子卿, 刘光明, 吴富梅. 2011. 2000中国大地坐标系:中国大陆速度场. 测绘学报, 40(4): 403-410.
吴富梅, 刘光明, 魏子卿. 2012. 利用局域欧拉矢量法建立CGCS2000速度场模型. 武汉大学学报·信息科学版, 37(4): 432-435.
谢方, 程传录, 王斌, 等. 2015. 中国大陆速度场模型建立方法研究. 大地测量与地球动力学, 35(2): 258-260, 272.
杨元喜, 曾安敏, 吴富梅. 2011. 基于欧拉矢量的中国大陆地壳水平运动自适应拟合推估模型. 中国科学:地球科学, 41(8): 1116-1125.
姚宜斌, 刘经南, 施闯, 等. 2002. ITRF97参考框架下的中国大陆区域地壳板块运动背景场的建立及应用. 武汉大学学报·信息科学版, 27(4): 363-366, 376.
游新兆, 杨少敏, 王琪. 2012. 中国大陆地壳运动观测研究进展. 地震, 32(2): 31-39.