地球物理学报  2013, Vol. 56 Issue (9): 3022-3028   PDF    
引用本文
冯威, 黄丁发, 李萌 等. 高频GPS双差残差模型监测强震地表运动. 地球物理学报, 2013, 56(9): 3022-3028, doi: 10.6038/cjg20130915.  
Feng W, Huang D F, Li M, et al. Ground motion monitoring during strong shake with high-rate GPS double-differenced residual model. Chinese J. Geophys. (in Chinese), 2013, 56(9): 3022-3028, doi: 10.6038/cjg20130915.  
高频GPS双差残差模型监测强震地表运动
冯威 , 黄丁发 , 李萌 , 张熙 , 严丽     
西南交通大学测量工程系, 成都 610031
收稿日期 2012-10-25, 2013-01-16 收修定稿
基金项目: 国家高技术研究发展计划(863计划)项目(2012AA12A209);国家自然科学基金项目(40771173)资助
作者简介: 冯威, 男, 1984年生, 讲师, 研究方向:GPS数据处理.E-mail:fengwei99@gmail.com
通讯作者: 黄丁发, 男, 1963年生, 教授, 主要从事GNSS理论与系统开发研究.E-mail:dfhuang@swjtu.edu.cn
摘要: 根据震时地表震动持续时间短的特点, 以及双差残差的相关性, 提出了基于双差残差建模的高频GPS强震位移监测方法.首先对震前若干历元各双差残差进行建模, 地震发生后, 利用模型对各双差残差进行预报, 最后利用预报残差实现短时间尺度内的位移解算.利用一条长约1100 km基线的静态数据和El Mayor-Cucapah 7.2级地震的94个测站数据进行试验分析, 结果表明, 5 min预测时间内, 静态数据的动态位移在N、E、U三个方向的定位中误差分别为6 mm、6 mm和13 mm, 地震数据的解算结果与实际情况有较好的一致性.
关键词: 高频GPS      双差残差模型      强震地面运动      实时动态监测     
Ground motion monitoring during strong shake with high-rate GPS double-differenced residual model
FENG Wei, HUANG Ding-Fa, LI Meng, ZHANG Xi, YAN Li     
Department of Surveying and Mapping, Southwest Jiaotong University, Chengdu 610031, China
Abstract: Due to the short lasting time of the ground motion and the correlation of the double-differenced (DD) residuals, a new high-rate GPS ground motion monitoring method is proposed based on DD residual modeling for strong shake. Firstly, the DD residual is modeled with epochs of data before the earthquake, then the model is used to predict the DD residual after the earthquake start, and finally the displacement is accomplished with the predicted residual in short time scale. 1 Hz GPS data from a static baseline about 1100 km and 94 sites in El Mayor-Cucapah Mw7.2 earthquake is tested, and the experimental results show that the positioning precision of N, E and U can reach 6 mm, 6 mm and 13 mm for the static data, and the results derived from earthquake data coincide with the actual situation well, within 5 predicted minutes..
Key words: High-rate GPS      Double differenced residual model      Strong shake ground motion      Realtime kinematic monitoring     
1 引言

随着GPS接收机技术的发展和高频接收机的普及, 高频GPS技术在强震领域得到了越来越广泛的研究和应用[1-5].且与传统地震仪相比, 由于不受振幅限制, GPS能够有效监测大幅度的地面震动, 其可为地震预警、救援和相关理论研究提供可靠的信息.

为了满足强震地表运动监测的高精度定位和单历元解算的要求, 常用精密单点定位(PPP)模式和差分定位模式来实现高精度的位移解算.由于地震影响的范围较大, 采用PPP技术进行地震的地表位移解算[6]具有一定的优势.但由于PPP技术对卫星钟差的精度要求很高, 而目前还难以实时获取到精密的卫星钟差信息, 因此PPP常用于事后的位移解算.差分定位模式虽然还无需精密卫星钟, 但在地震监测应用中需在附近选择一个参考站, 而选择较远的测站作为参考站时, 数据处理过程变得相对复杂.目前广泛使用的单历元数据处理软件有GAMIT/ GLOBK软件的TRACK模块[7], BERNESE软件的运动学坐标解算模块[8], GIPSY软件的高精度运动学数据处理模块[9], 以及Geodetics公司研制的商业软件RTD[10].虽然这些软件能够实现单历元高精度的位移解算, 但它们难以同时在实时数据处理以及长基线高精度解算方面满足地震地表实时监测的要求.

本文基于震时地表震动持续时间短的特点, 以及短时间内双差残差的强相关性, 提出了双差残差预报的GPS位移监测方法, 对各双差残差进行建模, 再结合预报残差实现短时间尺度内的位移监测.最后本文将通过长约1100km的静态基线数据和ElMayor-Cucapah7.2级地震时94个测站的数据来验证本文方法的有效性.

2 双差残差建模

长基线GPS双差观测方程可表达成如下形式:

(1)

其中Φ为经天线相位中心、对流层模型等改正后的相位观测值, λ为对应的波长, N为整周模糊度, T为对流层残差, I为电离层延迟, O为轨道误差, M为多路径延迟, εIF为观测噪声, i(i=1, 2)代表不同频率, 代表双差运算.与常用的测站天顶延迟估计加映射函数的模型不同, 对流层延迟先通过经验模型改正[11], 再对各卫星的双差对流层残差分别进行建模.

卫星轨道误差所导致的单差测距误差可表示为[12]:

(2)

式中:O为卫星位置误差向量(m), Xmk为基线向量(km), θ为两个向量的夹角, ρ为卫地距离.国际GNSS服务组织(IGS)提供的超快速精密星历的精度已达到5cm, 根据式(2), 对于1000km长的基线, 超快速精密星历对单差观测值所引入的误差不超过3mm.因此, 使用精密星历计算卫星位置时, 其误差可忽略不计.

式(1)中, 电离层延迟可通过双频消电离层组合消除其影响, 而多路径延迟则可通过改进的恒心日滤波的方法来建模[13-14].经过消电离层和恒心日滤波处理后, 式(1)可写成如下形式:

(3)

式中其中为载波相位的频率, 1和2分别代表GPS的L1和L2载波频率, εIFIF为转换后的噪声.

式(3)中的消电离层双差模糊度不再具有整周特性, 但是在没有周跳发生的情况下, 其是不随时间变化的量.一般情况下, 天顶对流层延迟量约为2.4m, 可以通过模型有效削除其中90%以上的延迟影响. 图 1显示了不同卫星高度角下映射函数的值.当高度角高于25°时(图中的虚线处), 映射函数值随高度角的变化十分缓慢; 当卫星高度角低于20°时, 随着高度角的降低, 映射函数值急剧增大, 因此式(3)中不同卫星高度角的双差对流层残差需采用不同的模型.经过双差处理后, 对流层残差将得到进一步削弱.对流层残差具有极强的时空相关性, 测站天顶延迟变化常是以小时为时间尺度来估计[15], 同时考虑到地震发生时地表震动时间相对较短, 卫星高度角变化也较小, 由映射函数所引起的误差可近似表达成线性函数.结合图 1, 双差残差一般可简化成如下模型:

(4)

图 1 不同卫星高度角的映射函数值 Fig. 1 The values of mapping function of different elevation

式中t为时间, θ为卫星高度角, θ0为不同模型的分界值, 建议值为25°, εr为模型噪声.由于模糊度是一个定值, 将模糊度合并到式(4)后有

(5)

称a和b为双差残差参数, f(t)为双差残差模型.式(3)代入式(5)可得:

(6)

地震发生前, 测站精确坐标已知, 即式(6)中的可视为已知量.通过地震发生前的若干历元, 结合最小二乘法即可解算出双差残差参数a、b.

3 基于双差残差预报的地表运动监测

当地震发生时, 式(6)中的将会发生变化, 可以认为距震中较远(上千公里)的参考站其空间位置保持不动, 的变化是由于离震中较近的测站震动而引起.由于天顶延迟通常是以小时为间隔来估计其变化量, 而地震地表震动时间短, 可忽略其变化量.再结合第2节映射函数相关影响的分析, 可以认为地表震动过程中双差残差不发生突变, 各双差残差模型f(t)保持不变.图 2为4.1节中静态试验部分某两颗卫星的双差残差序列, 可以看出在整个过程中双差残差量变化较为平缓.当然在大气变化异常时也可能会出现异常情况, 此时可以在建模和定位过程中进行适当的判断, 探测是否存在异常.根据上面分析, 式(6)可表达为:

(7)

图 2 不同卫星双差残差序列图 Fig. 2 Time series of the double differenced residual value for different satellites

f0(t)为震前所确定的各双差残差模型f(t)的预报值.当一个历元有n颗卫星(参考星除外)可用时, 可得如下观测方程:

(8)

其中上标代表不同的卫星, 利用最小二乘法就能解算出震时各个历元的地表位移.

4 试验结果分析

为检验本文方法的有效性, 分别采用静态数据和真实地震数据, 利用本文方法进行动态定位.试验过程中, 采用5min的数据来估计各双差残差模型, 并用各模型预报后5 min的双差残差, 利用预报的双差残差来解算测站的动态位移, 解算时采用IGS的预报精密星历.

4.1 静态试验

静态数据采用两个CORS站的观测数据, 基线长约1100km, 采样间隔为1s, 时长为24h, 卫星截止高度角设置为5°.以每10min为一个时段进行分段处理, 并剔除卫星数少于5的时段.

图 3显示了根据预报的双差参数模型解算出的测站动态位移的分布情况, 共约解算了38000个历元的位移.统计结果表明, N、E、U三个方向的中误差分别为6mm, 6mm, 13mm.

图 3 静态基线的动态定位结果分布 Fig. 3 Kinematic positioning results distribution for the static baseline

进一步分析双差残差模型, 图 4显示了两个具有代表性的双差参数的计算值和模型值, 其中, 双差残差模型是根据前300s的数据解算得出, 后300s的模型值为预报值.可以看出, 当卫星高度角较低时, 双差残差变化显著, 设置为常数显然不合理.另外为了防止高次模型因噪声影响而出现震荡现象, 应尽量采用低次的模型.因此本文根据不同的卫星高度角分别用一次模型和常数模型来对双差残差建模.

图 4 不同卫星高度角的双差残差模型 (a)48°; (b)8°. Fig. 4 The double-differenced residual model for different elevations

随着预报时间的增加, 双差残差的相关性将会减弱.对于一般地震来说, 地面震动时候一般不超过5 min, 因此本文将只分析在残差预报的时间达到5min时, 本文方法进行位移监测所能达到的精度.为了削弱随机噪声的影响, 将根据最后20个历元的位移信息取平均来计算最终的位移量.本试验共解算了126个时段的数据, 结果如图 5所示.

图 5 预报5min时的位移偏差 Fig. 5 The offset of the displacement at the 5th predicted minute

可以看出, 水平方向的最大值不超过20 mm, 竖直方向基本不超过40 mm.统计结果表明, N、E方向超过92.0%的位移小于10 mm, U方向超过91.2%的位移小于20mm.

4.2 ElMayor-Cucapah地震位移解算

2010年4月4日, 在墨西哥Baja地区发生了7.2级强烈地震, 主震持续时间超过了40s, 地表破裂长度达120km.本文利用加州实时观测网络(CRTN)的1 Hz GPS数据[16], 通过本文的方法解算该地震各测站的动态位移和地震同震位移场, 参数设置与静态试验相同.

本次试验共解算了94个测站的位移, p090离震中约920km, 将其作为参考站.主震之后的部分位移场图如图 6所示.

图 6 CRTN各测站的同震位移 Fig. 6 Coseismic displacement of CRTN sites

图 7画出了其中离震中较近和较远的几个测站水平方向的实时动态位移信息, 可以看出, 利用本文的方法解算出的GPS时间序列可清楚记录各个测站的地表震动.

图 7 测站的水平位移时间序列 Fig. 7 Times eries of the horizontal displacement

为验证解算结果, 将本方法解算的同震位移与GPSExplorer提供的高精度后处理结果[16]进行对比, 统计了其中66个测站的同震位移的差值, 统计结果如图 8.可以看出水平方向偏差基本小于20mm, 竖直方向偏差基本小于30 mm.统计结果表明, 与后处理结果相比, N、E方向分别有79%和91%的测站相差小于10mm, U方向有82%的测站相差小于20mm.

图 8 本文方法与后处理方法解算的位移之差 Fig. 8 Displacement difference between the proposed method and the post-process method

进一步分析N方向差值最大的4个测站, 偏差从大到小依次为p494、p744、p497和p498.从图 6中可以看出这4个点距震中较近, 且位置相对集中.由于事后处理的结果是根据震后几天的数据解算得出, 其位移场是震后几天的位移场, 而本文解算的是震后几分钟的位移场.两种解算结果在4个测站N、E方向的位移差分别为:-34、-28、-25、-19mm和9、-3、3、0 mm.可以看出它们主要呈现出往南偏的趋势, 这与同震位移的方向大体保持一致, 初步分析两种方法解算结果中较大的差异主要是由于震后滑移所引起.

SOPAC提供的p494和p497两个测站(p744、p498缺少)的快速静态解和本文的解算结果吻合较好, 两测站N方向相差约为8mm和14mm, E方向相差约2mm和16mm, 进一步证明了该方法的有效性.

5 结论

根据地震地表震动时间短的特点、以及短时间尺度内双差残差具有极强的相关性, 提出了基于双差残差预报的强震地表运动实时监测方法.与传统的基线解算方法不同, 该方法不单独解算模糊度信息, 而是直接对各双差残差建模, 并利用残差模型预报实现上千公里基线的解算.试验结果表明, 在5min的预报时间内, 1Hz的静态数据该方法N、E、U三个方向的中误差分别为6mm、6mm和13mm, ElMayor-Cucapah地震数据的解算结果也与实际情况有较好的一致性, 试验结果证明了本文方法的可行性.本方法只需要预报精密星历, 且计算效率高, 在有实时观测数据的情况下, 可实现强震时测站位移的实时监测, 这对地震预警和地震快速救援具有一定的实用和研究价值.

需要注意的是双差参数预报的精度会随着预报时间的增加而降低, 本文方法不适用于长时间的地表运动缓慢的监测.本文初步测试了5 min预报时间内方法的有效性, 对于更长时间的位移监测, 方法的有效性尚需进一步研究.此外, 周跳对本方法的影响不容忽视, 其将影响残差模型和后面的位移解算.

参考文献
[1] Shi C, Lou Y D, Zhang H P, et al. Seismic deformation of the Mw8.0 Wenchuan earthquake from high-rate GPS observations. Advances in Space Research , 2010, 46(2): 228-235. DOI:10.1016/j.asr.2010.03.006
[2] Elósegui P, Davis J L, Oberlander D, et al. Accuracy of high-rate GPS for seismology. Geophys. Res. Lett. , 2006, 33(11). DOI:10.1029/2006G1_026065
[3] Langbein J, Bock Y. High-rate real-time GPS network at Parkfield: Utility for detecting fault slip and seismic displacements. Geophys. Res. Lett. , 2004, 31(15). DOI:10.1029/2003GL019408
[4] 殷海涛, 张培震, 甘卫军, 等. 高频GPS测定的汶川Ms8.0级地震震时近场地表变形过程. 科学通报 , 2010, 55(23): 2529–2634. Yin H T, Zhang P Z, Gan W J, et al. Near-field surface movement during the Wenchuan Ms8.0 earthquake measured by high-rate GPS. Chinese Science Bulletin (in Chinese) , 2010, 55(23): 2529-2634. DOI:10.1007/s11434-010-4026-2
[5] Crowell B W, Bock Y, Squibb M B. Demonstration of earthquake early warning using total displacement waveforms from real-time GPS networks. Seismological Research Letters , 2009, 80(5): 772-782. DOI:10.1785/gssrl.80.5.772
[6] 方荣新, 施闯, 徐培亮, 等. GPS地震仪: PANDA软件测试结果与验证. 武汉大学学报(信息科学版) , 2011, 36(4): 453–456. Fang R X, Shi C, Xu P L, et al. GPS seismometer: PANDA software testing results and validation. Geomatics and Information Science of Wuhan University (in Chinese) , 2011, 36(4): 453-456.
[7] http://geoweb.mit.edu/~tah/track_example/
[8] Bock Y, Nikolaidis R, de Jonge P J, et al. Instantaneous geodetic positioning at medium distances with the global positioning system. J. Geophys. Res. , 2000, 105(B12): 28233-28253.
[9] Jet Propulsion Laboratory California Institute of Technology. GIPSY5.0 Release Note.https://gipsy-oasis.jpl.nasa.gov/gipsy/docs/Release_Nates_5.0.pdf, 2008
[10] http://www.geodetics.com/products/software/rtd_pro.php
[11] Saastamoinen J. Atmospheric correction for the troposphere and stratosphere in radio ranging of satellites, in the use of artificial satellites for geodesy. Geophysics Monograph , 1972, 15(16): 247-251.
[12] Han S W. Carrier phase-based long-range GPS kinematic positioning. Sydney: The University of New South Wales , 1997.
[13] Choi K, Bilich A, Larson K M, et al. Modified sidereal filtering: Implication for high-rate GPS positioning. Geophys. Res. Lett. , 2004, 31. DOI:10.1029/2004GL021621
[14] Niell A E. Global mapping functions for the atmosphere delay at radio wavelengths. J. Geophys. Res. , 1996, 101(B2): 3227-3246. DOI:10.1029/95JB03048
[15] 熊永良, 黄丁发, 丁晓利, 等. 虚拟参考站技术中对流层误差建模方法研究. 测绘学报 , 2006, 35(2): 118–121. Xiong Y L, Huang D F, Ding X L, et al. Research on the modeling of tropospheric delay in virtual reference station. Acta Geodaetica et Cartographica Sinica (in Chinese) , 2006, 35(2): 118-121.
[16] http://geoapp03.ucsd.edu/gridsphere/gridsphere?cid=El+Mayor+Cucapah