一、模型研制的目的、必要性和需求分析

1.1 研制的目的

整合现有的大地测量数据，获取中国地震科学实验场地壳运动速度和应变率场模型，为研究实验场其他地震科学问题和地震预测提供形变基础数据和模型支撑。

1.2 必要性和需求分析

地震孕育时间尺度可达数百年甚至上千年，无法为地震仪器观测到，但孕震过程中断层缓慢积累弹性应变，并造成地表变形，反过来，以震间地表形变为约束就可以研究孕震过程中断层的闭锁状态和能量积累速率。因此，以大地测量资料为基础的形变模型是进行地震孕育机理研究和地震预测的重要手段。

形变模型是地震科学实验场的基础模型之一，能够为其他地震科学问题研究提供必要的支撑。在地震科学实验场15个主要科学问题中，与形变模型强相关的有6个，主要表现在：

1）为研究实验场中主要活动断裂应力应变累积过程、地震引起的库仑应力变化提供必要的变形数据约束（科学问题2、9）；

2）为断层结构模型、块体划分、断层分段和级联破裂提供形变依据（科学问题3、4、6）；

3）为地震预测模型提供必需的输入数据（科学问题11）。

二、模型预期指标和技术路线

2.1 预期指标

1）提供地震科学实验场区域高精度的大地测量观测速度场，精度达到1mm/yr；

2）提供地震科学实验场区域速度场模型，精度达到1-2mm/yr；

3）提供地震科学实验场区域应变率场模型，精度达到。

1.2 技术路线

    形变模型研制主要包括大地测量数据的收集与处理、GPS速度场解算、速度场模型建立、应变率场模型建立、结果验证等过程，其技术路线流程如下图所示。

图1 形变模型技术路线图

三、模型数据来源和质量分析

3.1 GPS数据来源（合计517个测站）

1)  “中国地壳运动观测网络”和“中国大陆构造环境监测网络”（简称“陆态网络”一期/二期）1998-2018年间观测的350个区域站和53个连续站，其中部分区域站包括“中国综合地球物理场观测”专项于2010、2012、2014和2016年复测的数据;

2)  973项目“活动地块边界带动力过程与强震预测”布设的30个测站，这些测站受973项目资助观测时间为2005-2007年，并于2015和2017年受国家自然科学基金项目“鲜水河断裂震间形变时空分布特征及地震危险性分析”资助进行了两期复测;

3)  中国地震局地震预测研究所布设的龙门山连续观测网络，合计14个站，采用数据的观测时间约为2008-2018年;

4)  中国地震局地震预测研究所在鲜水河断裂布设的4个跨断层连续观测站，观测时间为2013-2018年;

5)  Kreemer等人于2014年在G-cubed杂志发表的论文中提供的66个GPS测站速度场（Kreemer et al., 2014）。

3.2 质量分析

最终采用的GPS区域站至少观测3期，连续站至少观测2年，均按照中国地壳运动观测网络的相关技术标准进行观测，所有原始观测数据均能满足地壳运动研究需求。

四、模型研制关键技术方法

4.1 GPS速度场解算方法

使用武汉大学研发的PANDA软件（Shi et al., 2015）和精密单点定位模式，对上述GPS原始观测资料进行统一的高精度数据处理，获得了IGS08参考框架下的单日解坐标时间序列。由于中国地震科学实验场及周边地区地震频发，在速度场计算中，我们考虑了周边2004/2012年苏门答腊地震、2001年可可西里地震、2008年汶川地震、2011年日本地震、2015年尼泊尔地震等大震以及实验场区所有Mw5.5级以上地震的影响。地震目录来源于GCMT数据库，地震影响的初始范围采用Herring et al. (2016)提供的经验公式。

在处理中，我们根据地震发生的时间分别计算震前和震后的速度。若震前、震后速度相差小于2mm/yr，则采用所有时段的坐标序列估计同震偏移和震间速度；若震前震后速度差异大于2mm/yr，则直接采用震前速度。由于部分测站只有震后观测，难于判断震后影响的大小，对于这些测站我们直接采用震后速度,并给予标记。对于区域站，采用线性模型估计速度和同震变形；对于连续站，同时估计线性速度、同震变形和年周期及半年周期的季节性形变。在所有解算中，均采用自适应的稳健估计方法探测观测值的粗差，并进行速度和同震变形模型的检验。最后采用欧拉矢量将所有GPS速度转换到欧亚参考框架。

Kreemer et al.（2014）对全球范围内可以获得的GPS原始观测资料进行统一数据处理，同时融合了233项研究工作的GPS测站速度，形成了一个较全面的全球GPS速度场。为了丰富中国地震科学实验场的GPS资料，我们利用“陆态网络”与Kreemer等人发布的公共站点计算欧拉矢量，然后利用该欧拉矢量将Kreemer等人发布的速度场转换到我们采用的参考框架下，从而形成统一参考框架的速度场。

4.2 速度场和应变率场模型计算方法

    将实验场区域划分为精细的三角形网格（网格大小约为15-30公里），通过球面函数建立网格节点速度与三角形内观测点速度的函数模型，依据最小二乘法则反演节点的速度。通过解算得到的网格节点的速度值，代入球面函数模型中能够获得研究区域任一点的速度值，从而建立起空间连续分布的速度场。在速度场的基础上，采用Savage et al.(2001)球面应变率计算方法，建立实验场区域的应变率场模型。详细方法介绍见Wang and Wright (2012)。

五、模型完成情况与成果

5.1 GPS速度场

目前共获得了中国地震科学实验场517个测站的GPS速度场，如图2所示，具体数据见附件。

图2.（A）GPS速度场及其受地震影响的程度。其中红色箭头表示震前与震后速度差异较大，采用震前速度；蓝色箭头表示震前与震后速度差异较小，采用所有数据估计的速度；绿色箭头表示只有震后观测；紫色箭头表示只有震前观测。（B）GPS速度场观测值的期数（区域站）或年数（连续站）。

5.2 形变模型

依据上述GPS观测速度，计算得到实验场区域的速度场模型（图3A-3B）和应变率模型（图3C），具体数据见附件。

图3. (A-B)速度场模型的东向和北向分量，其中圆圈里的颜色代表GPS观测速度的相应分量值。（C）应变率张量的第二不变量，其中蓝点为历史地震。

六、模型验证（测试）与精度评价

在Zheng et al. (2017)基础上，一方面增加了973、龙门山、鲜水河等项目的站点，提高了测站密度；另一方面，“陆态网络”站点新增了2017和2018年观测，部分站点增加了“中国综合地球物理场观测”专项于2010、2012、2014、2016的观测，提高了观测结果的精度和可靠性（图4）。站点的观测期数如图2B所示。为了保证结果的可靠性，最终的GPS速度场中所有站点至少观测3期。比较发现，Zheng et al. (2017)有38个站点与本结果差异达到3mm/yr以上，这些站点的可靠性得到显著改善。最终结果中，连续站的中误差约为0.4mm/yr，流动站的中误差约为1mm/yr。各测站东向和北向速度中误差见附件。

在建立速度场模型中，格网节点的东向和北向速度中误差均为1.4mm/yr，应变率场的中误差平均值为，约为应变率平均值的15%。速度场和应变率场模型的精度如图5所示。GPS东向和北向速度残差（RMS）分别为1.7mm/yr和1.0mm/yr。

图5. (A-B)速度场模型的东向和北向分量中误差，其中圆圈里的颜色代表GPS观测速度的相应分量的中误差。（C）应变率张量的第二不变量的中误差。

为了进一步验证速度场模型的精度，我们从观测值中随机抽取50个测站作为交叉检验，而用剩余的测站建立速度场模型。结果显示，这50个测站的东向和北向速度观测值与模型预测值之差分别为-0.3±1.8mm/yr和0.2±1.0mm/yr（图6）。

图6. 随机抽取的50个测站空间分布及其（A）东向和（B）北向速度预测残差

七、模型使用说明

形变模型包括三部分：GPS观测速度场、GPS速度场模型、应变率场模型。其中GPS观测速度场为文本文件(cses_gps19.dat)，其格式参见文件中的说明,GPS速度场模型(ve.grd, vn.grd)和应变率场模型（2nd_strain.grd）均为GMT的GRID格式，能够用GMT或其他GIS软件打开。

八、应用案例

研究团队2017年在JGR杂志上发表的GPS速度场（Zheng et al., 2017）两年共被论文引用37次。上述最新版本的实验场形变模型为首次提供，尚未有应用案例。

参考文献：

[1]   Herring, T. A., Melbourne, T. I., Murray, M. H., Floyd, M. A., Szeliga, W. M., King, R. W., et al. (2016). Plate boundary observatory and related networks: GPS data analysis methods and geodetic products. Reviews of Geophysics, 54, 759–808. https://doi.org/10.1002/2016RG000529

[2]   Kreemer, C., G. Blewitt, E. C. Klein (2014). A geodetic plate motion and Global Strain Rate Model. Geochemistry, Geophysics, Geosystems, 15, 3849-3889.

[3]   Savage, J. C., Gan, W., & Svarc, J. L. (2001). Strain accumulation and rotation in the Eastern California Shear Zone. Journal of Geophysical Research: Solid Earth, 106(10), 21,995-22,007.

[4]   Shi, C., M. Li, Q. Zhao, C. Wang, H. Zhang (2015). WHU Analysis Center Technical Report 2014, IGS TECHNICAL REPORT 2014, 81–88.

[5]   Zheng, G., H. Wang, T. J. Wright, Y. Lou, R. Zhang, W. Zhang, C. Shi, N. Wei (2017). Crustal deformation in the India-Eurasia collision zone from 25 years of GPS measurements. Journal of Geophysical Research: Solid Earth, 122, 9290–9312.

[6]   Wang H., Wright T.J. Satellite geodetic imaging reveals internal deformation of western Tibet, Geophysical Research Letters, 2012, 39, L07303, doi:10.1029/2012GL051222.