Construction of Velocity Field in Chinese Mainland Based on Local Seamless Delaunay Triangulation with Inverse Distance Weighting Method
-
摘要: 现有研究中，中国大陆速度场模型容易出现数据量小、连续性不佳、现势性不强等问题。基于近7年的中国地壳运动观测网络工程连续运行基准站GNSS(global navigation satellite system)观测数据，解算得到高精度陆态网基准站点的点位坐标和速度场，并利用提出的局部无缝Delaunay三角网反距离加权模型构建中国大陆格网速度场。相较现有的NNR-NUVEL1A(no net rotation Nubia velocity 1A)欧拉矢量模型与常用的速度场数值拟合模型，局部无缝Delaunay三角网反距离加权模型的精度最高，其水平速度场拟合精度优于1.5 mm/a; 采用陆态网约1 800个区域站进行外部检核，结果表明，东、北方向速度差值的绝对值平均值分别为1.11 mm/a、0.90 mm/a，中误差分别为1.50 mm/a、1.35 mm/a。相较于整体三角网而言，局域三角网在大陆边缘地区的拟合精度更优，水平速度场的插值精度更高，平均可以提高约0.3 mm/a，且不易产生粗差。局部无缝Delaunay三角网反距离加权模型不仅考虑了邻近点的距离和方位信息，还可以刻画出更为精细的局部特征，同时克服了边缘地区整体三角网跨度过大以及二级块体边缘处三角网不连续的缺点。Abstract:Objectives With the development of modern millimeter level geodetic technology, velocity field is one of the indispensable basic elements in Geodetic work. Under the background that current velocity models have limitations such as small data, poor continuity and timeliness update, a method based on local seamless Delaunay triangulation with inverse distance weighting is presented to establish the velocity field model in Chinese mainland.Methods Aiming at the limitations of the current velocity models, we obtain the high-precision coordinates and velocities of each CMONOC (crustal movement observational network of China) station, through processing the long-term CMONOC GNSS (global navigation satellite system) observation data of recent seven years by the GAMIT/GLOBK software. And we construct the velocity field of Chinese mainland grid with the inverse distance weighted model of locally seamless Delaunay triangulation. Compared with that calculated by NNR-NUVEL1A (no net rotation Nubia velocity 1A) Euler vector model and some other numerical fitting models, the velocity obtained by the new method in this research has the highest precision.Results The fitting precision of horizontal velocity field constructed by this new model is better than 1.5 mm/a. By using more than 1 800 regional stations for external checking, it shows that the average of absolute value of velocity differences in East and North directions are 1.11 mm/a and 0.90 mm/a, respectively, and the corresponding mean square errors are 1.50 mm/a and 1.35 mm/a. In continental margin areas, the fitting precision of horizontal velocity filed by the method of local triangulation net is much higher than that from the whole triangulation net model. Furthermore, based on the method of local triangulation net, interpolation precision of horizontal velocity filed is higher, an improvement of about 0.3 mm/a on average, and the gross errors could be largely reduced.Conclusions The new model, which is based on the method of local seamless Delaunay triangulation with inverse distance weighting, can effectively utilize both distance and orientation information of the adjacent points, depict subtler local features with high precision. And it overcomes the shortcomings of the large span of the whole triangulation network in the edge area and the discontinuity of the triangulation network at the edge of the secondary block.
-
随着现代毫米级大地测量技术的发展，速度场是大地测量工作中不可或缺的基础性要素之一。站点坐标会由于板块运动、潮汐位移等因素产生变化，因此，必须通过速度场将坐标归算到用户所需的历元时刻，才具有工程实用价值。目前，关于速度场的拟合主要有两种方法：一种是采用转换断层方位角、地震走滑矢量和磁异常得到的扩张速率等数据确定全球相对板块运动模型，如NNR-NUVEL1A(no net rotation Nubia velocity 1A)模型[1]; 另一种是使用近几十年来采集的中国大陆空间大地测量数据确定中国大陆速度场。后者能够较好地展现局部板块的整体运动趋势，精度较高，已成为目前速度场拟合的主流方案。
2011年12月中国大陆构造环境监测网络(简称陆态网)完成全部建设任务，建成了覆盖中国大陆及近海的高精度、高时空分辨率的中国地壳运动监测网络(crustal movement observation network of China，CMONOC)，它是推算中国大陆速度场的主要数据源。因此，不少学者利用网络数据开展了对中国大陆速度场的研究，常用的速度场模型建立方法分为欧拉矢量法和数值拟合法。吴富梅等[2]系统论述了基于CGCS2000(China geodetic coordinate system 2000)坐标系构建速度场模型的方法，对比了构建速度场的全局欧拉矢量法、局部欧拉矢量法和格网平均值法; 文献[3-4]研究了中国大陆区域在ITRF97(international terrestrial reference frame 1997)框架下的速度场，并采用多面函数法拟合速度场; 文献[5-7]基于地壳网络数据对比了常用的数值拟合方法的精度; 于亮等[8]基于2011—2015年的CMONOC数据并采用克里金法构建了中国大陆格网速度场模型; 苗岳旺等[9]基于1 397个区域站信息计算了中国大陆整体和各省的欧拉矢量。但是，由于使用的数据较少，或数据连续性不强，或缺乏明确的应用背景，或数据现势性不强等各种原因，现有的研究成果的实用性价值仍有一定的局限性。恰当的数值拟合方法在反映局部块体的细微变化方面较欧拉矢量法更为准确。目前，常用的数值拟合方法有克里金法、距离加权法、多元回归法、三角网插值法以及有限元法等[5-6]，其中三角网插值法结合反距离加权模型，使用灵活，且可以较好地刻画块体的局部特征，精度较高。
鉴于此，本文对陆态网基准站2011—2017年近7年的连续观测数据进行了处理与分析，采用几种常用的模型(如普通克里金法、反距离加权法和三角网内插法)构建了中国大陆速度场，在此基础上提出了基于局部无缝Delaunay三角网反距离加权模型，并对比分析了几种速度场模型的拟合精度和适用特点，证明了本文模型的高精度和实用性。
1 基于GAMIT/GLOBK软件的陆态网高精度数据解算
1.1 陆态网与周边IGS站点的选取
陆态网是以卫星导航定位系统观测为主，辅以其他多种技术手段，主要由均匀分布在全国的260个连续运行基准站和2 000个不定期观测区域站构成。其中，观测站的大部分观测设备为天宝公司的NET R8/R9型双频接收机和Dorne Margolin扼流圈天线，可以获得30 s、1 s和50 Hz采样率的观测数据。考虑到陆态网基准站建站规范、维护得当、运行稳定，且数据产出连续，其质量有一定的保障，本文采用陆态网的260个连续运行基准站2011—2017年的GNSS(global navigation satellite system)观测数据进行速度场解算。
为了将中国国内基准站同国际IGS(International GNSS Service)站进行联合解算，需选取国内及周边的IGS站进行约束。基于连续性、稳定性和高精度性等原则，挑选17个IGS站(AIRA、BJFS、DAEJ、IISC、IRKT、KIT3、KUNM、LHAZ、PIMO、POL2、SHAO、URUM、WUHN、YSSK、TIXI、ARTU和TCMS)进行数据处理。2 000个区域站、260个基准站和17个IGS站的分布见图 1。
1.2 高精度CGCS2000基准站坐标与速度的解算
1) 陆态网站点速度解算
首先，采用GAMIT(GPS analysis at massachusetts institute of technology)软件计算陆态网站点坐标和卫星轨道的单日松弛解。将陆态网260个站点与选取的17个周边IGS站进行联合解算。为了便于模糊度的分解与确定，陆态网站点和IGS站点坐标需给定适当的约束。同时，将陆态网基准站分为7个区域，每个区域单独进行解算，各区域之间以选取的17个IGS站作为公共连接点。然后，将每个分区的单日松弛解法方程文件与SOPAC(Scripps Orbit and Permanent Array Center)提供的全球IGS站的单日松弛解法方程文件进行合并，得到单日整体松弛解。最后，利用GLOBK(global Kalman filter)软件将所有的单日整体松弛解法方程文件进行合并平差，获得ITRF2014(international terrestrial reference frame 2014)框架下陆态网站点的坐标和速度信息，并根据解算的站点坐标及速度递推得到2000.0历元时刻的坐标。需要注意的是，为了给不同法方程文件选定合适的比例，需进行两次合并平差，第一次合并时各法方程文件比例因子为1，平差后提取平差结果中各文件的卡方作为比例因子，重新合并平差。
2) ITRF框架转换
顾及到ITRF的定向时间演变条件和测站高程的周期性变化等因素，将测站速度从空间直角坐标系转换为站心坐标系，且令高程方向速度为0;再将得到的速度由站心坐标系转换为空间直角坐标系; 利用表 1中的2010.0历元时刻(小数年)的转换参数与速率递推得到2000.0历元时刻ITRF2014至ITRF97的转换参数，再经过历元统一和框架统一两个步骤[10]，可以将ITRF2014框架下的坐标和速度转换到CGCS2000坐标系所参考的ITRF97框架下的坐标和速度，从而获得陆态网260个基准站在ITRF97下2000.0历元时刻的坐标和速度。
表 1 2010.0历元时刻ITRF2014至ITRF97的转换参数Table 1. Transformation Parameters from ITRF2014 to ITRF97 at Epoch 2010.0项目 Tx/mm Ty/mm Tz/mm D Rx/(″) Ry/(″) Rz/(″) 参数值 7.4 -0.5 -62.8 3.8×10-9 0 0 0.000 26 年变率 0.1 -0.5 -3.3 1.2×10-8 0 0 0.000 22 注：Tx、Ty、Tz分别为框架间3个坐标轴平移量，D为框架间的尺度变化因子，Rx、Ry、Rz分别为框架间3个坐标轴旋转角度 2 基于欧亚板块运动的欧拉矢量模型
2.1 欧拉矢量模型的转换精度
根据国际NNR-NUVEL1A模型提供的欧亚板块的欧拉矢量参数[1](wx = -0.000 81 rad/Ma，wy = -0.002 394 rad/Ma，wz = 0.003 153 rad/Ma)，将解算得到的中国大陆速度场转换至欧亚参考框架下。为直观分析NNR-NUVEL1A模型的转换效果，转换前后的速度场见图 2。
由图 2可以明显看出，中国大陆西部地区的转换速度场精度低于东部地区，说明大陆西部地区特别是青藏、川滇地区的地壳形变较大，东部地区的地壳形变较小，因此扣除NNR-NUVEL1A参数的变化量后，东部地区的速度几乎为0。计算NNR-NUVEL1A模型的E、N方向的速度误差统计值分别为11.2 mm/a、8.7 mm/a，误差较大，无法满足速度场转换要求。
2.2 欧拉矢量模型的不足之处
虽然NNR-NUVEL1A模型能够消除中国大陆速度场的部分运动趋势，但是由于NNR-NUVEL1A模型代表的是近几百万年地壳欧亚板块的平均运动，不能代表其近期运动情况，因此国内外有学者研究了欧亚板块运动的现今趋势和中国大陆地区的地壳形变特征[11-15]。另外，由于中国大陆相对于欧亚板块也有一定的运动，因此必须建立适合中国大陆板块的整体欧拉模型，以便能较好地展现其整体运动趋势。吴富梅等[2]利用网络点速度得到了大陆地区的整体欧拉矢量，虽然站点位置速度可达7.7 mm/a的精度，但仍旧无法满足高精度速度场转换的要求。从实用角度出发，本文通过数值拟合方法建立格网速度场，可以达到2 mm/a的水平速度场转换精度。
3 中国大陆速度场数值拟合方法
3.1 普通克里金法
克里金法[5]又称空间局部插值法，是以变异函数理论和结构分析为基础，在有限区域内对区域化变量进行线性无偏最优估计的一种方法，是最常用的空间插值算法。由于克里金插值公式复杂，本文仅介绍基于普通克里金法进行速度插值的主要步骤。方法如下：
1) 计算离散点数据集$ {Q}_{i}({x}_{i}, {y}_{i}, {z}_{i})(i=\mathrm{1, 2}\cdots n, {x}_{i}, {y}_{i}\mathrm{为}\mathrm{位}\mathrm{置}\mathrm{信}\mathrm{息}, {z}_{i}\mathrm{为}\mathrm{属}\mathrm{性}\mathrm{信}\mathrm{息}) $中任意两点之间的距离$ {d}_{ij} $和属性半方差值$ {r}_{ij} $，公式如下：
$$ {d}_{ij}=\sqrt[]{{\left({x}_{i}-{x}_{j}\right)}^{2}+{\left({y}_{i}-{y}_{j}\right)}^{2}} $$ (1) $$ {r}_{ij}=\frac{1}{2}{\left({z}_{i}-{z}_{j}\right)}^{2}, i, j=\mathrm{1, 2}\cdots n $$ (2) 2) 在数据对$ ({d}_{ij}, {r}_{ij}) $中寻找最优拟合函数(变异函数理论模型中有球形模型、指数模型、高斯模型等[6]，经试验，3种模型的速度场拟合结果基本一致，本文最终选用的是高斯模型)拟合距离$ d $与模型$ r $之间的函数关系，得到最优拟合函数的关系式为：
$$ r=r\left(d\right)=\left\{\begin{array}{l}0, d=0\\ {C}_{0}+C\left(1-{\mathrm{e}}^{-\frac{d}{a}}\right), d>0\end{array}\right. $$ (3) 式中，$ {C}_{0} $为函数常量; $ C $为指数系数; $ a $为衰减因子。
3) 根据最优拟合函数计算得到所有已知离散点集之间的半方差$ {r}_{ij} $。
4) 对未知点$ {z}_{0} $根据最优拟合函数计算得到未知点到已知点之间的半方差$ {r}_{i0} $。
5) 代入下列方程组利用最小二乘原理求解最优权重系数$ {\lambda }_{i} $：
$$ \left\{\begin{array}{l}{r}_{k0}-\sum\limits _{j=1}^{n}{\lambda }_{j}{r}_{kj}+u=0, k=\mathrm{1, 2}\cdots n\\ \sum\limits _{i=1}^{n}{\lambda }_{i}=1\end{array}\right. $$ (4) 式中，$ u $为拉格朗日乘子。
6) 求出权重系数$ {\lambda }_{i} $后，代入目标方程求出待定点属性值$ {\widehat{z}}_{0} $：
$$ {\widehat{z}}_{0}=\sum \limits_{i=1}^{n}{\lambda }_{i}{r}_{i0} $$ (5) 3.2 反距离加权法
反距离加权法[6]又称谢别德法，是以地理学第一定律——相似相近原理为基础，根据近邻点的属性来确定未知点属性。设空间待插值点为$ P({x}_{P}, {y}_{P}, {z}_{P}) $，P点邻域内有n个已知离散点$ {Q}_{i}({x}_{i}, {y}_{i}, {z}_{i}) $。根据离散点的值，通过反距离加权法对P点的属性$ {z}_{P} $进行插值，即：
$$ {z}_{P}=\sum\limits_{i=1}^{n}\frac{{z}_{i}}{[{d}_{i}{(x, y)]}^{k}}/\sum\limits_{i=1}^{n}\frac{1}{[{d}_{i}{(x, y)]}^{k}} $$ (6) 式中，$ {d}_{i}\left(x, y\right)=\sqrt[]{({x}_{P}-{x}_{i}{)}^{2}+({y}_{P}-{y}_{i}{)}^{2}} $，表示离散点到待定点的距离; $ k $一般取1~2，本文取$ k $=1。
反距离加权法以一定邻域内的已知点确定待定点信息，但因其仅仅以距离作为权重标准而忽视了方位信息，所以其精度有一定的限制。
3.3 Delaunay三角网反距离加权法
3.3.1 Delaunay三角网生成
在数字地形建模中，不规则三角网通过不规则分布的离散点并按照某种规则生成连续的三角网，以近似逼近地形表面。相对于格网模型，不规则三角网能在某一特定分辨率下利用更少的信息表示更复杂的表面。在所有可能的三角网中，Delaunay三角网在地形拟合中的应用较为普遍，不仅算法成熟，而且生成的三角网结构良好，数据冗余度小，存储效率高，便于信息检索、分析和提取[16-17]。
进行Delaunay三角网剖分时，必须要满足两个重要的准则：(1)Delaunay三角网具有唯一性(即任意4点不可共圆)。在Delaunay三角网中，任一三角形的外接圆内不会有其他点存在; (2)最大化最小角特性。在离散点集可能形成的三角剖分中，Delaunay三角剖分所形成的三角形的最小角最大，即Delaunay三角网是最接近规则化的三角网，具体是指两个相邻的三角形构成凸四边形的对角线相互交换后，6个内角的最小角不再增大。
根据实现过程的不同，Delaunay三角网通常分为逐点插入法、三角网生成法、分治算法和合成算法等[18-19]。其中，逐点插入法原理简单，易于编程实现，基本步骤为：(1)建立包含所有离散点集的多边形凸包; (2)在构成凸包的所有点中，初始化凸包三角网; (3)根据局部优化过程(local optimization procedure，LOP)算法优化初始凸包三角网; (4)在已生成的三角网逐点插入其余散点，每次插入一个点后，用LOP算法重新生成三角网; (5)重复步骤(4)，直至所有点插入完毕。
3.3.2 整体Delaunay三角网反距离加权法
以陆态网中260个基准站点构建中国大陆整体Delaunay三角网，如图 3所示; 通过射线法判断待定点位于哪个三角形中; 在该三角形中，利用顶点的位置和属性信息，根据反距离加权法(式(6))确定待插入点的属性信息。
![]()
图 3 中国大陆整体Delaunay三角网Figure 3. Overall Delaunay Triangulation Network in Chinese Mainland与反距离加权法相比，整体Delaunay三角网反距离加权法不仅以距离作为权重标准，而且依据三角形特性考虑了内插点各个方向的属性信息，但其明显存在的问题是，边缘区域三角网跨度过大，违反了相邻相近性原理。
3.3.3 局部Delaunay三角网反距离加权法
根据中国大陆的二级板块划分成果，以射线法判断陆态网站点分别位于哪个板块中并进行分类统计; 逐一在二级板块中利用所在区域的站点构建局部Delaunay三角网，如图 4所示; 以射线法判断待定点位于哪个三角形中; 在该三角形中，利用顶点的位置和属性信息，根据式(6)描述的反距离加权法确定待插入点的属性信息。
![]()
图 4 中国大陆二级板块Delaunay三角网Figure 4. Delaunay Triangulation Network of the Ⅱ Plates in Chinese Mainland与整体Delaunay三角网反距离加权法相比，该方法的局部刻画更精细，同时克服了大陆边缘区域三角网跨度过大的缺点，但由于站点个数有限，各块体三角网之间存在不连续的缺点，即当待定点位于板块三角网缝隙之间时，无法找到对应的三角形。
3.3.4 局部无缝Delaunay三角网反距离加权法
为了克服上述局部Delaunay三角网反距离加权法的缺点，本文提出一种结合整体和局部三角网的反距离加权法。该方法结合整体和局部三角网优点的同时，克服了三角网跨度过大、块体三角网之间不连续的缺点，即当待定点落于局部Delaunay三角网中，采用局部Delaunay三角网反距离加权法; 若待定点不在局部Delaunay三角网中，则采用整体Delaunay三角网反距离加权法，以确定位置点属性信息。采用局部无缝Delaunay三角网反距离加权法拟合速度场的流程如图 5所示。
![]()
图 5 局部无缝Delaunay三角网反距离加权法流程图Figure 5. Flowchart of Local Seamless Delaunay Triangulation Network with Inverse Distance Weighting Model4 速度场数值拟合模型的精度分析
为了便于速度场的发布和使用，将中国大陆按照1°的格网间隔分为1 000个左右的格网，分别采用普通克里金法、反距离加权法、整体Delaunay三角网反距离加权法(简称整体三角网法)和局部无缝Delaunay三角网反距离加权法(简称局部无缝三角网法)4种数值拟合方法计算出每个格网点的平均速度值。以局部无缝Delaunay三角网反距离加权法为例，计算并绘制中国大陆格网速度场见图 6。
4.1 格网速度场的内部检核
为了检验格网速度场模型的精度，利用260个基准站点进行内部检核。这里为了更直观地体现模型的整体精度，数据处理中删除了部分残差较大的点，实际使用的检核点个数少于260个。分别利用4种格网模型内插出基准站点的速度，将内插出的E、N方向的速度值与原始精确的速度值进行比较，结果如图 7所示。4种数值拟合模型得到的格网速度的内部检核精度如表 2所示，参数分别为E、N方向速度差值绝对值的平均值(MvE、MvN)及其速度差值的中误差(SvE、SvN)。
![]()
图 7 基准站点E、N方向插值速度与实测速度之差Figure 7. Differences Between the Interpolated and Measured Velocities of Reference Stations in East and North Directions表 2 4种模型的内部检核精度比较/(mm·a-1)Table 2. Comparison of Internal Checking Accuracy of Four Models/(mm·a-1)模型 MvE MvN SvE SvN 普通克里金法 0.67 0.64 1.07 1.09 反距离加权法 0.90 0.81 1.32 1.22 整体三角网法 0.59 0.62 0.89 1.00 局部无缝三角网法 0.56 0.56 0.87 0.99 4.2 格网速度场的外部检核
利用选取的中国大陆范围内观测时间较长的约1 800个区域站点进行模型的外部检核。分别用4种格网模型内插出区域站点的速度值，将内插出的区域站点E、N方向的速度值与实测得到的精确速度值进行比较，结果如图 8所示。4种数值拟合模型得到的格网速度的外部检核精度见表 3。
![]()
图 8 区域站点E、N方向插值速度与实测速度之差Figure 8. Differences Between the Interpolated and Measured Velocities of Regional Stations in East and North Directions表 3 4种模型的外部检核精度比较/(mm·a-1)Table 3. Comparison of External Checking Accuracy of Four Models/(mm·a-1)模型 MvE MvN SvE SvN 普通克里金法 1.27 0.95 1.64 1.44 反距离加权法 1.30 1.00 1.67 1.46 整体三角网法 1.15 0.92 1.52 1.38 局部无缝三角网法 1.11 0.90 1.50 1.35 经过内外部检核，对4种格网速度场数值拟合模型E、N方向的精度进行分析，得到如下结论：
1) 内部检核精度高于外部检核精度，4种格网模型的精度具有相似的对比结果，精度由低到高排序为：反距离加权法 < 克里金法 < 整体三角网法 < 局部无缝三角网法。
2) 数值拟合方法虽然物理意义不及欧拉矢量法，但生成的格网速度场使用简便，且精度整体优于欧拉矢量法，水平速度拟合精度优于2 mm/a。
3) 由于仅考虑距离作为权重标准，反距离加权法虽原理简单，但拟合精度较低; 普通克里金法作为地理统计常用的一种方法，考虑了随机数据的期望和方差等信息，因此精度有一定的提高，但其原理较为复杂，特别是半变异函数的选取会直接影响结果的精度，要想得到理想的结果，需要不断地尝试与验证。
4) 三角网法的拟合精度可达1.5 mm/a，相较于以上两种插值方法，三角网法在考虑距离和方位信息的同时，Delaunay三角剖分的特性使得采用更少的信息便可以刻画更为精细的局部，因此拟合精度较高。
5) 本文提出的局部无缝Delaunay三角网反距离加权法结合整体和局部三角网优点的同时，克服了整体三角网跨度过大、二级板块三角网不连续的缺点，其精度最高。但由于检核点大部分位于大陆区域的内部，位于跨度过大的边缘三角网处的站点较少，因此局部无缝Delaunay三角网反距离加权法的优势体现得不是特别明显，几乎与整体Delaunay三角网反距离加权法的精度处于同一水平。
6) 文献[5]中将中国大陆划分为五大块体并采用三角网进行插值，其水平速度场的外部插值精度也达到了1.5 mm/a的水平，但其检核点数量较少，且没有考虑大块体边缘三角网跨度过大和块体三角网之间不连续的问题，因此使用上有一定的局限。本文提出的局部无缝三角网法很好地解决了这一问题。
4.3 边缘地区站点整体与局域三角网模型拟合精度分析
在本文介绍的3种三角网加权方法中，整体三角网法覆盖了整个大陆区域，适应范围较广，但局部地形特征描述不足，特别是大陆边缘地区三角网跨度过大，将严重破坏相邻相近性原则，导致该区域速度场插值精度较低; 二级块体三角网可以较好地刻画局部特征，但同时存在块体边缘三角网不连续的问题; 局部无缝三角网模型结合了这两种模型的优点，在块体内部采用二级块体三角网模型且在块体边缘区域采用整体三角网模型进行拟合。
为了进一步对比研究整体与局域三角网模型在大陆边缘区域构建的速度场趋势一致性水平，本文在中国大陆边缘地区均匀选取20个测站点，分别采用两种模型内插出各区域站点的速度值，将内插出的区域站点在E、N方向上的速度值与实测速度值进行比较，结果如图 9所示。显然，局域三角网模型的插值精度明显优于整体三角网模型：整体三角网模型在E、N方向的速度差平均值分别为1.17 mm/a、0.99 mm/a，而局域三角网模型在E、N方向的速度差平均值分别为0.88 mm/a、0.72 mm/a。由此说明了局域三角网模型可以较好地克服边缘地区整体三角网跨度过大带来的影响。
![]()
图 9 大陆边缘站点E、N方向插值速度与实测速度之差Figure 9. Differences Between the Interpolated and Measured Velocities of Stations Located in Continental Margin Areas in East and North Directions另外，本文提出的局部无缝Delaunay三角网反距离加权法是基于二级块体的划分成果进行中国大陆速度场的构建，由于各块体自身内部(除地壳内部活动较为剧烈的青藏和川滇块体)的整体运动趋势保持一致，若相邻两点处于同一二级块体内部，则两处站点的速度值基本一致; 若相邻两点恰好位于二级块体边界且处于不同的块体，则有可能出现两处站点的速度值差异较大的情况。当相邻两点处于青藏或川滇地区，若两点处于同一三角网内，则两处站点的速度值基本一致; 若两点恰好处于不同三角网内，则也有可能出现两处站点的速度值差异较大的情况。
5 结语
本文以实现中国陆态网基准站CGCS2000坐标和速度场为背景，基于2011—2017年的CMONOC连续运行基准站观测数据，以GAMIT/GLOBK软件为计算平台，介绍了中国大陆速度场的处理与获取方法，得到了高精度的陆态网基准站点的定位成果和速度场数据，提出了基于局部无缝Delaunay三角网反距离加权法拟合任意点速度的方法，并对几种速度场拟合模型构建的格网速度场进行了对比。主要结论有：
1) NNR-NUVEL1A模型仅能消除中国大陆速度场的部分运动趋势，不能代表其近期运动情况，精度最差。
2) 数值拟合方法虽然缺乏物理意义，但整体精度优于欧拉矢量法，其中基于局部无缝Delaunay三角网反距离加权法构建的格网速度场使用简便，精度最高，其内、外部检核精度分别可达1 mm/a、1.5 mm/a。
3) 在4种数值拟合模型中，由于三角网法考虑了距离和方位信息的同时，可以刻画更为精细的局部特征，因此三角网法相较于反距离加权法和克里金法精度更高。局部无缝Delaunay三角网反距离加权法结合整体和局部三角网优点的同时，克服了边缘三角网跨度过大、二级板块三角网不连续的缺点，精度最高。
4) 采用不同的数值拟合方法构建中国大陆格网间距为1°的速度场，既满足了精密速度场转换的需求，也满足了工程实用的需求，对具体工程应用有一定的参考价值。局部无缝反距离加权模型主要是基于大陆二级板块进行局域地区特征的描述，但是个别二级板块内部的地壳特征仍具有多样性，像川滇地区板块内部活动较为复杂，以其为整体构建三角网可能会降低内部小块体边缘处站点的速度插值精度。因此，在本文的研究基础上可进一步细化二级块体或者以省市级为块体进行三角网建模，以获得更加理想的速度场模型。
致谢: 感谢中国大陆构造环境监测网络提供的部分数据支撑。 -
![]()
图 3 中国大陆整体Delaunay三角网
Figure 3. Overall Delaunay Triangulation Network in Chinese Mainland
![]()
图 4 中国大陆二级板块Delaunay三角网
Figure 4. Delaunay Triangulation Network of the Ⅱ Plates in Chinese Mainland
![]()
图 5 局部无缝Delaunay三角网反距离加权法流程图
Figure 5. Flowchart of Local Seamless Delaunay Triangulation Network with Inverse Distance Weighting Model
![]()
图 7 基准站点E、N方向插值速度与实测速度之差
Figure 7. Differences Between the Interpolated and Measured Velocities of Reference Stations in East and North Directions
![]()
图 8 区域站点E、N方向插值速度与实测速度之差
Figure 8. Differences Between the Interpolated and Measured Velocities of Regional Stations in East and North Directions
![]()
图 9 大陆边缘站点E、N方向插值速度与实测速度之差
Figure 9. Differences Between the Interpolated and Measured Velocities of Stations Located in Continental Margin Areas in East and North Directions
表 1 2010.0历元时刻ITRF2014至ITRF97的转换参数
Table 1 Transformation Parameters from ITRF2014 to ITRF97 at Epoch 2010.0
项目 Tx/mm Ty/mm Tz/mm D Rx/(″) Ry/(″) Rz/(″) 参数值 7.4 -0.5 -62.8 3.8×10-9 0 0 0.000 26 年变率 0.1 -0.5 -3.3 1.2×10-8 0 0 0.000 22 注：Tx、Ty、Tz分别为框架间3个坐标轴平移量，D为框架间的尺度变化因子，Rx、Ry、Rz分别为框架间3个坐标轴旋转角度 
下载: 导出CSV
表 2 4种模型的内部检核精度比较/(mm·a-1)
Table 2 Comparison of Internal Checking Accuracy of Four Models/(mm·a-1)
模型 MvE MvN SvE SvN 普通克里金法 0.67 0.64 1.07 1.09 反距离加权法 0.90 0.81 1.32 1.22 整体三角网法 0.59 0.62 0.89 1.00 局部无缝三角网法 0.56 0.56 0.87 0.99
下载: 导出CSV
表 3 4种模型的外部检核精度比较/(mm·a-1)
Table 3 Comparison of External Checking Accuracy of Four Models/(mm·a-1)
模型 MvE MvN SvE SvN 普通克里金法 1.27 0.95 1.64 1.44 反距离加权法 1.30 1.00 1.67 1.46 整体三角网法 1.15 0.92 1.52 1.38 局部无缝三角网法 1.11 0.90 1.50 1.35
下载: 导出CSV
-
[1] Argus D F, Gordon R G. No-Net-Rotation Model of Current Plate Velocities Incorporating Plate Motion Model NUVEL-1[J]. Geophysical Research Letters, 1991, 18(11): 2 039-2 042 doi: 10.1029/91GL01532
[2] 吴富梅, 刘光明, 魏子卿. 利用局域欧拉矢量法建立CGCS2000速度场模型[J]. 武汉大学学报·信息科学版, 2012, 37(4): 432-435 http://ch.whu.edu.cn/article/id/169 Wu Fumei, Liu Guangming, Wei Ziqing. Velocity Field Model of CGCS2000 Based on Euler Vector of Local Area[J]. Geomatics and Information Science of Wuhan University, 2012, 37(4): 432-435 http://ch.whu.edu.cn/article/id/169
[3] 姚宜斌, 刘经南, 施闯, 等. ITRF97参考框架下的中国大陆区域地壳板块运动背景场的建立及应用[J]. 武汉大学学报·信息科学版, 2002, 27(4): 363-366 http://ch.whu.edu.cn/article/id/4973 Yao Yibin, Liu Jingnan, Shi Chuang, et al. Construction and Application of Chinese Crustal Motion Background Field Based on ITRF97 Frame[J]. Geomatics and Information Science of Wuhan University, 2002, 27(4): 363-366 http://ch.whu.edu.cn/article/id/4973
[4] 刘经南, 姚宜斌, 施闯. 中国地壳运动整体速度场模型的建立方法研究[J]. 武汉大学学报·信息科学版, 2002, 27(4): 331-336 http://ch.whu.edu.cn/article/id/4967 Liu Jingnan, Yao Yibin, Shi Chuang. Method for Establishing the Speed Field Model of Crustal Movement in China[J]. Geomatics and Information Science of Wuhan University, 2002, 27(4): 331-336 http://ch.whu.edu.cn/article/id/4967
[5] 张勇. 不同数值内插方法建立我国速度场模型[D]. 北京: 中国测绘科学研究院, 2011 Zhang Yong. Different Numerical Interpolation Method to Establish Velocity Model[D]. Beijing: Chinese Academy of Surveying and Mapping, 2011
[6] 任雅奇. 基于GPS数据的中国地壳运动速度场模型的建立[D]. 郑州: 信息工程大学, 2012 Ren Yaqi. Establishment of Chinese Crustal Movement Velocity Model Based on the GPS Data[D]. Zhengzhou: Information Engineering University, 2012
[7] 谢方. 基于2000国家大地坐标系的中国大陆速度场模型建立方法研究[D]. 西安: 长安大学, 2013 Xie Fang. Research on the Modeling Approach of Velocity Field in China Continent Based on China Geodetic Coordinate System 2000[D]. Xi'an: Chang'an University, 2013
[8] 于亮, 朱璇, 陈永祥, 等. 基于CMONOC建立和评估中国大陆地壳运动速度场模型[J]. 大地测量与地球动力学, 2017, 37(7): 704-708 https://www.cnki.com.cn/Article/CJFDTOTAL-DKXB201707007.htm Yu Liang, Zhu Xuan, Chen Yongxiang, et al. The Building and Estimating of Velocity Field in Chinese Mainland Based on CMONOC[J]. Journal of Geodesy and Geodynamics, 2017, 37(7): 704-708 https://www.cnki.com.cn/Article/CJFDTOTAL-DKXB201707007.htm
[9] 苗岳旺, 李军, 朱璇, 等. 基于欧拉矢量的中国大陆局部速度场建模方法[J]. 测绘工程, 2019, 28(1): 26-31 https://www.cnki.com.cn/Article/CJFDTOTAL-CHGC201901006.htm Miao Yuewang, Li Jun, Zhu Xuan, et al. The Method of Local Velocity Field Modeling in Chinese Mainland Based on Euler Vector[J]. Engineering of Surveying and Mapping, 2019, 28(1): 26-31 https://www.cnki.com.cn/Article/CJFDTOTAL-CHGC201901006.htm
[10] 张杰, 范玉磊, 董海政, 等. ITRF框架与CGCS2000坐标转换的研究[J]. 全球定位系统, 2017, 42(6): 79-83 https://www.cnki.com.cn/Article/CJFDTOTAL-QUDW201706014.htm Zhang Jie, Fan Yulei, Dong Haizheng, et al. Research on Coordinate Transformation Between ITRF Frame and CGCS2000[J]. GNSS World of China, 2017, 42(6): 79-83 https://www.cnki.com.cn/Article/CJFDTOTAL-QUDW201706014.htm
[11] Zhao Bin, Huang Yong, Zhang Caihong, et al. Crustal Deformation on the Chinese Mainland During 1998—2014 Based on GPS Data[J]. Geodesy and Geodynamics, 2015, 6(1): 7-15 doi: 10.1016/j.geog.2014.12.006
[12] 张小红, 鄂栋臣. 西南极菲尔德斯海峡断层GPS地壳形变监测网的建立和精度分析[J]. 测绘地理信息, 2000(1): 26-30 https://www.cnki.com.cn/Article/CJFDTOTAL-CHXG200001007.htm Zhang Xiaohong, E Dongchen. The Establishment and Precision Analysis of GPS Crustal Deformation Monitoring Network in the Fildes Strait Fault of Southwest Polar Region[J]. Journal of Geomatics, 2000(1): 26-30 https://www.cnki.com.cn/Article/CJFDTOTAL-CHXG200001007.htm
[13] 瞿伟, 王运生, 张勤, 等. 空间大地测量GPS揭示的汾渭盆地及其邻域现今地壳应变场变化特征[J]. 地球物理学报, 2016, 59(3): 828-839 https://www.cnki.com.cn/Article/CJFDTOTAL-DQWX201603006.htm Qu Wei, Wang Yunsheng, Zhang Qin, et al. Current Crustal Deformation Variation Characteristics of the Fenwei Basin and Its Surrounding Areas Revealed by GPS Data[J]. Chinese Journal of Geophysics, 2016, 59(3): 828-839 https://www.cnki.com.cn/Article/CJFDTOTAL-DQWX201603006.htm
[14] Qu W, Lu Z, Zhang Q, et al. Crustal Deformation and Strain Fields of the Weihe Basin and Surrounding Area of Central China Based on GPS Observations and Kinematic Models[J]. Journal of Geodynamics, 2018, 120: 1-10 doi: 10.1016/j.jog.2018.06.003
[15] Yu J, Tan K, Zhang C, et al. Present-day Crustal Movement of the Chinese Mainland Based on Global Navigation Satellite System Data from 1998 to 2018[J]. Advances in Space Research, 2019, 63(2): 840-856 doi: 10.1016/j.asr.2018.10.001
[16] 徐武平, 邱峰, 徐爱萍. 空间数据插值的自动化方法研究[J]. 武汉大学学报·信息科学版, 2016, 41(4): 498-502 doi: 10.13203/j.whugis20140131 Xu Wuping, Qiu Feng, Xu Aiping. Automatic Method of Kriging Interpolation of Spatial Data[J]. Geomatics and Information Science of Wuhan University, 2016, 41(4): 498-502 doi: 10.13203/j.whugis20140131
[17] 董箭, 彭认灿, 郑义东. 利用局部动态最优Delaunay三角网改进逐点内插算法[J]. 武汉大学学报·信息科学版, 2013, 38(5): 613-617 http://ch.whu.edu.cn/article/id/2648 Dong Jian, Peng Rencan, Deng Yidong. An Improved Algorithm of Point-by-Point Interpolation by Using Local Dynamic Optimal Delaunay Triangulation Network[J]. Geomatics and Information Science of Wuhan University, 2013, 38(5): 613-617 http://ch.whu.edu.cn/article/id/2648
[18] 青文星, 陈伟. Delaunay三角网生成的改进算法[J]. 计算机科学, 2019, 46(S1): 226-229 https://www.cnki.com.cn/Article/CJFDTOTAL-JSJA2019S1047.htm Qing Wenxing, Chen Wei. Delaunay Triangular Mesh Optimization Algorithm[J]. Computer Science, 2019, 46(S1): 226-229 https://www.cnki.com.cn/Article/CJFDTOTAL-JSJA2019S1047.htm
[19] 余杰, 吕品, 郑昌文. Delaunay三角网构建方法比较研究[J]. 中国图象图形学报, 2010, 15(8): 1158-1167 https://www.cnki.com.cn/Article/CJFDTOTAL-ZGTB201008005.htm Yu Jie, lv Pin, Zheng Changwen. A Comparative Research on Methods of Delaunay Triangulation[J]. Joural of Image and Graphics, 2010, 15(8): 1158-1167 https://www.cnki.com.cn/Article/CJFDTOTAL-ZGTB201008005.htm
-
期刊类型引用(12)
1. 成英燕，常春涛，王虎，徐彦田，许长辉. 基于CGCS2000新框架的地物坐标更新策略研究与分析. 测绘通报. 2024(03): 54-62 . 
百度学术
2. 丑凡. 三种速度场模型在CORS基准站时间序列中的对比分析. 测绘技术装备. 2024(01): 30-37 . 百度学术
3. 李宗勋，李宇平，余祖锋. 精密单点定位结果转换至CGCS2000坐标方法. 测绘通报. 2024(S1): 96-100 . 百度学术
4. 苏友能，徐青，朱新铭，邹霓霄，武钰清，刘波. POI点要素空间邻近图质量评价——以Delaunay三角网和缓冲区构建方法为例. 地理与地理信息科学. 2024(03): 14-20 . 百度学术
5. 曲麒霖，秘金钟，徐彦田，王春杨，韩忠倫，李庚辰. 基于监督聚类法的板块格网速度场构建方法研究. 测绘科学. 2024(12): 20-27 . 百度学术
6. 杨国创，尹彤，王金月，高严，孙泽宇. 多元回归格网改正量模型的大区域坐标转换方法. 导航定位学报. 2023(03): 185-190 . 百度学术
7. 杨承午. 中国区域速度场多面函数拟合稀疏建模方法研究. 铁道勘察. 2023(04): 27-34 . 百度学术
8. 党亚民，蒋涛，杨元喜，孙和平，姜卫平，朱建军，薛树强，张小红，蔚保国，罗志才，李星星，肖云，章传银，张宝成，李子申，冯伟，任夏，王虎. 中国大地测量研究进展（2019—2023）. 测绘学报. 2023(09): 1419-1436 . 百度学术
9. Yamin DANG，Hu WANG，Fuping SUN，Guangwei JIANG，Xinhui ZHU，Yingyan CHENG，Qiang YANG，Yingying REN，Jing JIAO. Maintenance of Millimeter-level Geodetic Reference Framework. Journal of Geodesy and Geoinformation Science. 2023(03): 9-18 . 必应学术
10. 赵立志. 基于放射性测量的离子吸附型稀土矿勘探方法设计. 煤化工. 2023(05): 141-144+149 . 百度学术
11. 许文洁，霍亮，沈涛，陈壮，赵博林. 地面沉降时空形变场的构建方法研究. 测绘科学. 2022(01): 204-211 . 百度学术
12. 孙付平，贾彦锋，朱新慧，肖凯，刘婧. 毫米级地球参考框架动态维持技术研究进展. 武汉大学学报(信息科学版). 2022(10): 1688-1700 . 百度学术
其他类型引用(5)
计量
- 文章访问数: 1386
- HTML全文浏览量: 382
- PDF下载量: 78
- 被引次数: 17
