中国是世界上地震灾害最严重的国家之一，随着社会经济的发展，地震造成的灾害损失逐步增大，且各种迹象表明我国正在进入新一轮的地震活动周期.如何有效减轻破坏性地震造成的人员伤亡和经济损失，成为我国迫切需要解决的问题.由于地震成因的复杂性，现有科技水平和认识水平下，短时间内还不能做到准确的地震预报，而地震预警则是目前国际上公认的能够减轻地震灾害的有效手段之一(张洪才等，2012).地震预警通常利用几秒钟的P波数据估算地震震级，这种方法估算中小型地震时十分有效，但是在估算强震震级时会出现震级饱和(Kuyuk and Allen, 2013).以3·11日本M_W9.0大地震为例，JMA地震预警系统出现饱和，在震后120 s左右速报的震级仅为M_W8.1，75 min后将震级提升为M_W8.4，3 h后才改为M_W8.8.由于对震级存在低估，导致没有对后续的海啸做出相应的防御措施，最终造成了重大的人员伤亡.

准确的地表位移在反演地震震级方面有非常重要的价值，测定位移的方法主要包括GPS和地震传感器两种，后者包含强震仪和测震仪.但是由于仪器倾斜等原因，强震仪积分得到的位移会出现基线偏移，需要进行事后的基线校正，测震仪在强震过程中还会出现饱和的现象.近几年实时高频GPS技术在地震预警中的应用得到越来越多的重视.与其他方法相比, 实时高频GPS方法最大的优点是在地震发生期间和之后，快速准确地获取地表位移数据，具有量程无限制、误差不累计的优势.同时，利用GPS提供的地表位移数据可以进行震源破裂过程以及断层滑动分布的反演，得到更加准确的矩震级.这将为人们逃生、紧急采取避难措施、政府采取应急救援准备等提供更好的保障.Allen和Ziv(2011)在模拟实时条件下利用EI Mayor-Cucapah地震时1 Hz GPS数据获取了地面震动过程，使用一个移动的时间窗口快速提取地面永久位移，对震级实现快速反演.Colombelli等(2013)在此基础上，提出自适应性方法，允许断层模型在破裂过程中向两端延伸.Li等(2014)对历元差分算法做出改进，并对3·11日本地震进行反演，结果表明实时单站位移解算方法与精密单点定位一样，可以得到可靠的矩震级和断层滑动分布.Ohta等(2012)基于RTK-GPS自动探测3·11日本地震的永久位移，解算断层平面和震级，最终预测海啸高度.Crowell等(2012)基于实时GPS数据，对比了使用断层目录提供的断层模型和快速矩张量解两种反演方法.结果表明，相比传统方法两者都可以更加快速得到强震震级.

本文以汶川地震为例，在模拟的实时条件下使用历元差分算法对高频GPS数据进行处理，提取地震波初至时刻，对震中位置进行反演.由于地震期间GPS数据出现中断，补充使用强震记录，通过二次积分模拟GPS地震波位移特征，探讨使用高频GPS记录中的PGA，PGV和PGD(加速度速度和位移峰值)计算地震烈度的可行性.最后利用强震记录提供的地表位移反演断层破裂过程和矩震级，验证地表位移数据对传统地震预警系统的重要补充作用.

需要说明的是，利用传统的地震学方法确定震源位置的技术已经非常成熟，且相比于GPS技术，地震传感器对地震更加灵敏，结果也更加快速准确.传统方法主要是依据地震动三要素(峰值、频谱、持时)计算地震烈度.强震仪频率可达500 Hz，远高于高频GPS(1 Hz)和超高频GPS(20~50 Hz).因此与GPS相比，强震记录明显包含更加丰富的地震波信息.本文使用高频GPS数据计算震中位置和地震烈度只是希望采用另一种独立的数据源辅助传统方法，提供多一种观测和外部检核条件.

在模拟实时条件下，使用历元差分方法对汶川地震期间13个高频GPS台站数据进行处理，提取各个台站的地震波初至时刻.再利用台站坐标和初至时刻反演震中位置、发震时刻和地震波传播速度.

使用历元差分方法计算相位变率，采用消电离层组合削弱电离层影响，采用改进的TurboEdit方法进行周跳的探测与修复，添加对流层延迟改正、相对论效应改正、固体潮改正、地球自转和相位缠绕改正，引入台站初始的精确坐标消除几何误差，结合广播星历，实现单站GPS速度的实时解算，对速度进行积分即可得到测站实时位移.在短时间内，由于广播星历、钟误差以及对流层误差等对解算结果不会产生太大影响，积分得到的测站位移精度可以达到cm级(鲁洋为和王振杰，2016).相比于相对定位，该方法不需要设置参考站，在地震应用中具有明显优势.

联合使用长短时平均(STA/LTA)方法和AIC方法提取地震波初至时刻.首先使用STA(信号短时平均值)和LTA(信号长时平均值)之比对地震波初至时刻进行粗提取.STA主要反映信号幅值的瞬时变化，衡量的是地震信号的变化，LTA反映了背景噪声的平均水平(马强，2013).地震波到达时，相应的STA/LTA会明显增加，若超过事先设定的阈值则认为地震波初至.在此基础上，利用AIC准则求解背景噪声与地震信号这两个时间序列的最佳划分点，作为地震波初至时刻.

震源位置是指地震破裂起始点的位置，其在地面上的投影称为震中.震中信息有利于搜索发震断层位置，最终反演断层破裂过程.因此地震发生之后，快速准确的震中信息可以为地震预警系统、地震灾情分析评估和震后救援工作提供重要帮助.理论上利用台站的三维坐标和地震波初至时刻可以反演震源三维坐标.但是由于GPS台站较少，且几乎都在同一个水平面上，受到几何结构的限制无法对震源深度提供有效约束.因此本文忽略震源深度的影响，使用地震波初至时刻，震中和测站之间的球面距离列观测方程，反演震中位置.震中距D_i可表示为：

(1)

其中，B_i和L_i分别为测站的纬度和经度，震中坐标为(B, L)，R是地球半径，D_i表示第i个测站的震中距.

假设地震波沿各个方向的传播速度相等, 设为v, 则根据地震波到达各个观测台站的时间差可以得到如下观测方程(方荣新等，2014)：

(2)

对观测方程线性化，依据最小二乘原理迭代计算即可得到震中坐标(B, L)和地震波速度v.震中位置的初值坐标可以设为任意值，与真实坐标越近收敛速度越快.GPS台站观测到的地震波信号主要是S波，波速一般为3~4 km·s^-1，因此将地震波速初值设为4 km·s^-1.地震的发震时刻T₀可根据震中坐标、测站坐标和地震波速度计算得到，并对多个测站结果取平均削弱粗差影响，公式为：

(3)

式中t_i为第i测站的地震波初至时刻.

选取四川省地震局和中国大陆地壳观测网络的13个高频1 Hz GPS测站的数据，采用历元差分方法获取震时测站三维速度和位移.其中PIXI台站记录了完整的主震过程，选取14: 28: 01为起始时刻，三维速度和位移序列如图 2所示.

图 2

Fig. 2

图 2 PIXI台站速度位移时间序列 Fig. 2 GPS velocity and displacement time series of station PIXI

为验证历元差分算法的精度，使用GAMIT/GLOBK软件中的TRACK运动学分析模块，以LUZH台站为参考站计算PIXI台站位移时间序列，并将两种结果作比较.其中东西向均方根误差为1.4 cm，南北向为2.3 cm，垂直向为2.4 cm.由于位移波形包含中长期低频噪声，且GPS方法垂直方向误差较高，采用测站水平向速度波形提取地震波初至时刻.采用上述STA/LTA方法进行粗提取，短时时间窗口设为2 s，长时时间窗口设为100 s，阈值设为10倍.使用Allen(1978)提出的特征函数放大震相特征，特征函数为：

(4)

式中v表示水平方向合速度.

考虑地震预警对时效性的要求，在粗判结果前取20个采样点，在粗判结果后取10个采样点，在该时间窗口内使用AIC准则精确判断地震波到时.AIC方法计算函数为(Maeda, 1985)：

(5)

式中N表示时间窗口总长度，k表示在时间窗口内逐点滑动的时间点，x表示水平合速度时程，var表示计算时间段内数据的方差.

由于大地震造成电力系统瘫痪，GPS数据出现中断.使用上述方法只有8个台站探测到了地震波初至时刻.其中QLAI站在北京时间14: 27: 43时在北方向出现一个3 cm左右的误差，使用GAMIT/GLOBK软件中的TRACK运动学分析模块处理，同样存在这个问题，认为是QLAI站数据本身存在缺陷，予以舍弃.受到数据不完整的限制，LESH站没有使用AIC方法.GPS台站分布如图 1所示，以14: 23: 00为0时刻，7个台站坐标和地震波初至时刻见表 1.

图 1

Fig. 1

图 1 GPS台站分布与震中位置示意图 黑色圆点表示台站位置，红色星号表示测震学的震中位置，黑色星号表示GPS数据给出的震中位置. Fig. 1 Distribution of GPS stations and epicentral locations Black dots are GPS stations. Red star is the epicenter constrained by seismic data. Black star is the epicenter constrained by GPS data.

表 1 Table 1

表 1 GPS测站坐标和地震波初至时刻 Table 1 The coordinates of GPS stations and seismic wave arrival times 测站 纬度(°N) 经度(°E) 地震波初至时刻(s) PIXI 30.910 103.757 326 CHDU 30.639 104.064 332 ZHJI 31.006 104.546 347 MYAN 31.440 104.726 351 RENS 30.200 104.103 352 JYAN 30.388 104.545 354 LESH 29.565 103.755 367

表 1 GPS测站坐标和地震波初至时刻 Table 1 The coordinates of GPS stations and seismic wave arrival times

利用测站坐标和地震波初至时刻，利用上述方法迭代反演震中位置，震中坐标为(31.1576°N, 103.4153°E)，发震时刻为14: 28: 12.6.张新东等(2013)使用测震学方法得到的震中坐标为(31.021°N, 103.375°E)，发震时刻为14: 27: 58.8，两个震中位置相差15.7 km，且指向破裂方向.分析产生误差的原因，主要有以下三点：(1)GPS台站全部位于龙门山断裂带下盘，数据分布并不具有完备的方位角；(2)我们采用的是1 Hz的高频GPS数据，最小时间单位为s，而地震波波速一般在3~4 km·s^-1，因而波到时中包含的误差将会对震中位置产生很大的影响；(3)地震仪具有极高的灵敏度，测震学方法给出的震中指向断层破裂的初始位置，这里称作仪器震中，而高频GPS的观测噪声较大，对微弱地表形变不敏感，震中结果更加靠近能量释放的位置，称作宏观震中.两种震中信息的不一致也是产生误差的重要原因之一.值得一提的是理论上可以根据两种不同震中位置判断断层破裂方向，这一点需要日后更多震例提供佐证.

根据野外调查结果，汶川地震烈度分布如图 3所示(吴健等，2008).其中为了使结果更加合理，将烈度边界附近的台站舍去，选用的数据见表 2.利用强震仪数据模拟高频GPS数据特征，通过统计结果表明两种数据的PGA/PGV存在正相关，但是基于同一公式计算的宏观烈度存在较大差异.最后统计了基于GPS方法的PGA、PGV和PGD与宏观地震烈度之间的关系式.由于使用的数据较少，需要日后做进一步研究.

图 3

Fig. 3

图 3 汶川地震宏观烈度与强震台站分布 Fig. 3 Macro intensity isoseismal map of the Wenchuan earthquake and locations of strong motion stations

表 2 Table 2

表 2 强震台站坐标与对应的烈度 Table 2 The coordinates of strong motion stations and their corresponding intensities 测站 纬度(°N) 经度(°E) 烈度 514ZG 29.3 104.8 6 51BXD 30.4 102.8 7 51BXY 30.5 102.9 7 51BXZ 30.5 102.9 7 51CXQ 31.7 105.9 7 51DXY 30.6 103.5 7 51DYB 31.3 104.5 7 51FSB 29.1 104.8 6 51GYS 32.2 105.8 7 51HSD 32.1 103.0 7 51HSL 32.1 103.3 7 51HYJ 29.5 102.6 7 51HYQ 29.6 102.6 6 51HYW 29.3 102.9 6 51HYY 29.7 102.4 6 51JYC 32.0 105.0 8 51JYD 31.8 104.7 8 51JZB 33.3 104.1 6 51JZG 33.1 104.3 7 51JZW 33.0 104.2 7 51JZZ 33.3 103.8 6 51KDG 30.0 101.6 6 51KDT 30.1 102.0 6 51LDD 29.6 102.2 6 51LDJ 29.7 102.2 6 51LDL 29.8 102.2 6 51LDS 29.9 102.2 6 51LSJ 30.2 102.9 7 51LXM 31.6 103.3 8 51LXS 31.5 102.9 7 51MBD 28.8 103.5 6 51MCL 29.0 103.7 6 51MED 31.9 102.2 6 51MEZ 31.9 102.3 6 51MNC 28.6 102.2 6 51MXD 32.0 103.7 8 51MZQ 31.5 104.1 10 51PJD 30.2 103.4 7 51PWM 32.6 104.5 8 51QLY 30.4 103.3 7 51SMC 29.1 102.3 6 51SMK 29.4 102.1 6 51SML 29.0 102.3 6 51SMM 29.3 102.4 6 51SMW 29.4 102.2 6 51SMX 29.3 102.3 6 51SPA 32.5 103.6 7 51SPC 32.8 103.6 6 51SPT 32.6 103.1 6 51TQL 29.9 102.4 6 51WCW 31.0 103.2 9 51XJB 31.0 102.4 7 51XJD 31.0 102.6 7 51YAL 29.9 102.8 6 51YAM 30.1 103.1 7 51YAS 29.9 103.0 6 51YBY 29.0 104.6 6 51YXX 28.7 102.5 6 51YXZ 28.6 102.5 6 51LXT 31.6 103.4 8

表 2 强震台站坐标与对应的烈度 Table 2 The coordinates of strong motion stations and their corresponding intensities

由于地震期间GPS数据出现中断，我们使用中国地震台网中心的强震记录模拟高频GPS数据特征，对利用高频GPS记录中的PGA/PGV计算地震烈度的问题进行探讨.强震记录和高频GPS数据主要有两点不同：(1)地震发生时，强震记录受到仪器倾斜、旋转等因素的影响，由积分得到的地表位移会出现基线偏移，而高频GPS没有这种现象；(2)强震记录一般为200 Hz，高频GPS为1 Hz，采样率不同也会导致两种数据的PGA/PGV存在较大差异.因此首先使用SMBLOC程序(Wang et al., 2011)改正强震记录中的基线偏移，单位为cm·s^-2.然后对基于强震记录得到的位移进行重采样，将采样率改为1 Hz，再对时间求一次和二次导数分别得到速度和加速度数据，以此来模拟震时高频GPS数据.

选取了60个强震台站数据，强震记录中的PGA和PGV分别记为As和Vs，采用上述方法模拟的GPS记录中的PGA，PGV和PGD分别记为Ag，Vg和Dg.对一定烈度内的As和Ag，Vs和Vg取平均值的对数，统计两组对数之间的关系，结果如图 4所示.需要说明的是只使用了汶川地震数据，烈度9和烈度10都只有一个台站.

图 4

Fig. 4

图 4 (a) As和Ag，(b) Vs和Vg之间的统计关系 Fig. 4 Statistical relationships between (a) As and Ag, (b) Vs and Vg

统计公式分别为：

(6)

(7)

由图 4可知，对不同来源的PGA/PGV取对数后，两种数据之间存在明显的正相关，说明高频GPS数据中的PGA/PGV也可以用于计算地震烈度.

1980年，刘恢先在中国地震烈度表中第一次列出了5~10烈度时所对应的PGA和PGV定量参考标准.PGA、PGV与烈度间的关系分别为(中国地震局，2008)：

(8)

(9)

将两种PGA/PGV分别代入上述公式中，在得到的烈度值I_As与I_Ag，I_Vs与I_Vg之间分别求差，以此来评价两种不同来源的PGA/PGV基于同一公式得到的烈度之间的差异.结果如图 5所示.

图 5

Fig. 5

图 5 基于两种不同PGA/PGV的烈度差值 Fig. 5 Intensity differences calculated by two kinds of PGA/PGV

由图 5可知，两种数据得到的烈度值存在较大差异.特别是利用PGA计算的结果，5个烈度的差值全部超过3.因此，在使用高频GPS数据计算地震烈度时，需要对原有公式作修改，才能实现两种类型数据的融合.

统计一定烈度下的PGA/PGV观测数据平均值的对数与地震烈度之间的关系，如图 6所示.在烈度为6、7、8时，PGA/PGV与地震烈度存在明显的线性关系，在烈度为9和10时较为离散.主要原因是本文只是用了汶川地震，烈度9和烈度10都只有一个台站，缺乏统计意义.该公式只是一个初步公式，需要日后利用更多震例、更多台站数据做修正.

图 6

Fig. 6

图 6 基于模拟GPS数据的PGA/PGV与宏观烈度之间的关系 Fig. 6 The relationship between macro intensity and PGA/PGV based on simulated GPS data

最终统计公式为：

(10)

(11)

作进一步分析，利用上述方法将强震数据中的位移数据进行重采样，计算各个台站的PGD.采取相同方法对处于同一烈度的PGD取平均值，统计PGD与地震烈度的关系.结果如图 7所示，基本特征与PGA/PGV相同，在烈度为9和10时受到台站数量的限制，结果较为离散，在6、7、8三个烈度时PGD与烈度的线性关系明显.该结果也证实了在台站数量足够多的情况下，实时高频GPS方法提供的PGD也可以作为一种有效数据源计算地震烈度，为传统方法作补充.

图 7

Fig. 7

图 7 基于模拟GPS数据的PGD与宏观烈度之间的关系 Fig. 7 The relationship between macro intensity and PGD based on simulated GPS data

统计公式为：

(12)

另外，根据GPS和强震记录中两种峰值间的定量关系，可以实现二者互相转化.首先使用式(6)和式(7)将GPS峰值转为强震记录中的峰值，再利用式(8)和式(9)计算烈度，结果如图 8所示.

图 8

Fig. 8

图 8 宏观烈度与峰值计算烈度的关系 Fig. 8 The relationship between macro intensity and the intensity calculated by peak values

结果表明利用两种峰值的定量关系，将GPS峰值转换为强震记录峰值同样可以用于计算地震烈度.由于只是用了汶川地震的数据，结果存在系统误差，需要日后更多震例来修正参数.

使用张国宏等(2010)根据野外地质考察研究成果和InSAR形变场特征提出的断层分段模型.但是后续的反演结果表明在主断层西南方端点处仍存在较大的滑动分布，参考中国活动断层目录(邓起东，2007)，将该模型中的主断层向西南方向延长至(30.81°E, 103.25°N)，如图 9所示，以此描述完整的断层滑动分布.

图 9

Fig. 9

图 9 断层模型与强震台站分布 蓝色线表示原有断层模型，红色线表示断层向西南方向延伸的部分.黑色三角表示反演时选取的40个强震台站. Fig. 9 The fault model and locations of strong motion stations Blue lines indicate original fault model. The red line indicates the part of fault extending to the southwest. Black triangles are strong motion stations.

考虑在实际应用中，使用历元差分算法实时获取的地表位移精度，选取三个台站水平位移首次超过5 cm的时刻作为反演的起算时刻.使用跟踪方差方法自动判断事件结束(Melgar et al., 2012).对每个测站的水平位移序列，逐历元计算当前历元之前20个历元的位移序列的方差，定义为跟踪方差.地震发生之后，跟踪方差会不断增大，随着地震波的结束，跟踪方差会逐渐下降至震前状态.根据跟踪方差这一特征设置判断条件，若当前历元的跟踪方差与之前历元中最大跟踪方差的比值小于某个阈值(本文设为0.25)，认为该站点事件结束.对事件结束时刻之后20 s的位移数据取平均值计算同震位移，反演得到的断层滑动分布作为最终反演结果.为了避免站点位移在地震期间存在短暂的平稳时期而导致判断的事件结束时刻过早，这里首先对水平位移最大的站点进行处理，再判断其他测站的结束时刻.若下一个历元搜索到拥有最大水平位移的站点发生变化，则舍弃之前的判断结果，重新计算.

选取40个强震台站的数据如图 9所示，利用SMBLOC程序去除积分后产生的基线偏移，得到各个测站震时地表位移时间序列.位移时间序列当中既有永久位移，也包含动态震荡的影响.这里引入一个移动的20 s平均窗口对位移序列作平滑，快速地提取永久位移，实时得到永久位移时间序列，最后使用SDM程序(Wang et al., 2008)实现断层滑动分布的准实时反演.采用上述方法得到起始和终止时刻分别为震后20 s和震后159 s，在震后179 s时输出最终震级.断层破裂过程如图 10所示，对应时刻的震级见表 3.

图 10

Fig. 10

图 10 断层滑动分布随时间的变化 从上到下分别表示震后20 s、30 s、40 s、60 s、70 s、80 s、159 s、179 s和强震数据完全静止时断层滑动分布情况. Fig. 10 Changes of fault slip distribution with time From top to bottom indicate the fault slip distributions at 20 s, 30 s, 40 s, 60 s, 70 s, 80 s, 159 s, and 179 s after the earthquake and the final motionless, respectively.

表 3 Table 3

表 3 震级随时间变化情况 Table 3 Changes of magnitude with time OT(s) 20 30 40 60 70 80 159 179 数据静止 MW 7.0 7.3 7.6 7.7 7.8 7.8 7.8 7.8 7.8 拟合度 83.61% 88.3% 95.98% 97.66% 97.27% 96.12% 93.79% 93.85% 93.83% 注：OT代表地震发生的起始时间.

表 3 震级随时间变化情况 Table 3 Changes of magnitude with time

由图 10可知，地震发生之后，主断层由西南向东北方向破裂，一共产生了四个明显的滑动区，与Zhang等(2012)研究成果相符合.震后30 s时，形成了第一个滑动区.此后破裂继续向东北方向延伸，到震后60 s时，又出现了两个明显的滑动区，在震后80 s时，滑动量继续增加，特别是第三个滑动区，滑动范围持续扩大，且出现了第四个滑动区.截止到跟踪方差方法提供终止时刻，即震后159 s，第四个滑动区继续扩展，且与最终断层滑动分布基本相同.将最终断层滑动分布与宏观烈度作对比，发现两者存在明显的对应关系.第一个和第三个滑动区滑动量最高，刚好对应烈度为11的区域，第二个滑动区的滑动量和范围相对较小，对应的烈度为10级.而第四个滑动集中区滑动量最小，其对应的区域烈度也最小为9级.震级方面，震后20 s提供了初始震级为7.0，随后地震波达到的台站逐步增加，到震后70 s达到7.8并稳定，在震后179 s作为最终震级输出.

高频GPS实时获取的地表位移数据可以为现有地震预警系统作补充，特别是在大地震(M_W＞7)发生时不会出现震级饱和.本文以汶川地震为例，研究分析了利用高频GPS快速确定震中、地震烈度，以及准实时反演滑动分布、震级的方法.得到如下认识：

(1) 使用7个GPS台站的坐标和长短时平均与AIC方法联合判断的地震波初至时刻反演了震中位置，坐标为(31.1576°N，103.4153°E)，以张新东使用测震学方法得到的震中结果为参考，两者相差15.7 km.误差产生的原因主要包含三个方面：测站分布不均匀，不具有完备的方位角；波到时包含的误差也会使计算结果与实际位置产生较大差距；两种方法确定的震相不一致导致两种震中指向不同位置；

(2) 利用强震数据模拟了高频GPS数据特征，探讨了基于高频GPS的PGA、PGV和PGD计算地震烈度的可行性，以汶川地震为例给出了一个初步经验公式.由于采用的是模拟数据且只有一次震例，该公式需要在今后的研究中进一步完善；

(3) 准实时反演断层破裂过程显示，汶川地震主断层由西南向东北方向破裂，共有4个明显的断层滑动量集中区域，与前人研究成果相吻合.在震后20 s提供初始震级M_W7.0，并逐步增大最终稳定为M_W7.8.结果表明地表位移数据在地震预警，灾情评估和灾后救援工作等方面具有重要应用潜力.但是我们也必须注意到位移数据对中小地震不敏感，因而在实际应用中，应该将大地测量与测震学方法相结合.地震发生之后，先使用测震学方法实现地震预警，若震级超过某一阈值，再使用大地测量方法更新预警信息，提高预警系统的准确性.