1、第4 3卷 第9期2 0 2 3年9月大 地 测 量 与 地 球 动 力 学J o u r n a l o fG e o d e s ya n dG e o d y n a m i c sV o l.4 3N o.9S e p t.,2 0 2 3收稿日期:2 0 2 2-1 1-1 8项目来源:国家自然科学基金(4 1 7 0 4 0 3 1,4 2 0 6 1 0 7 7);江西省数字国土重点实验室开放基金(D L L J 2 0 2 2 1 3);自然资源部环鄱阳湖区域矿山环境监测与治理重点实验室开放基金(MEM I-2 0 2 1-2 0 2 2-2 9)。第一作者简介:李萌,博士,讲
2、师,主要从事高精度GN S S数据处理和地壳变形研究,E-m a i l:n e m o n 8 1 81 6 3.c o m。D O I:1 0.1 4 0 7 5/j.j g g.2 0 2 3.0 9.0 0 6文章编号:1 6 7 1-5 9 4 2(2 0 2 3)0 9-0 9 0 9-0 5G P S监测慢滑移信号的双曲正切函数建模李 萌1,2 邹小平1,2 袁林果3 吕开云1,2 鲁铁定1,2 卢永刚41 东华理工大学江西省数字国土重点实验室,南昌市广兰大道4 1 8号,3 3 0 0 1 32 东华理工大学自然资源部环鄱阳湖区域矿山环境监测与治理重点实验室,南昌市广兰大道4
3、1 8号,3 3 0 0 1 33 西南交通大学地球科学与环境工程学院,成都市犀安路9 9 9号,6 1 1 7 5 64 中铁武汉勘察设计院有限公司,武汉市关山大道1号,4 3 0 0 7 4摘 要:以日本房总半岛2 0 1 1年、2 0 1 32 0 1 4年、2 0 1 8年3次慢滑移事件(s l o ws l i pe v e n t,S S E)为例,利用自主研发的G P S坐标时序分析(G P Sc o o r d i n a t et i m es e r i e sa n a l y s i s,G T S A)软件批量处理该半岛多个G P S站点2 0 0 92 0 1 9年
4、共1 0a的坐标时序数据,构建包含周期性S S E的坐标时序模型。结果表明:1)双曲正切函数的数学特性与S S E的运动规律十分吻合,表现为缓慢加速-快速滑动-缓慢减速,最后恢复至稳态;2)函数对S S E的精准建模可有效确定S S E发生时间的中心、持续时长和地表位移;3)2 0 1 1年S S E持续时间最长可达5 0d,2 0 1 32 0 1 4年和2 0 1 8年的S S E持续时间最长为2 5d;4)3次S S E的水平位移均朝东南方向,高程位移量均表现为下沉,2 0 1 8年S S E最大水平位移达5.2c m,2 0 1 32 0 1 4年最大高程位移达3.8c m。关键词:G
5、 P S;慢滑移;双曲正切函数;建模;坐标时序中图分类号:P 2 2 8 文献标识码:A 慢滑移事件S S E是一种持续数天至数月的无震滑动事件,发生在地壳内部薄弱带上,会引起地表形变,也可能导致慢地震的发生1。由于不产生地震波,微弱的地表S S E信号通常很难用传统地震仪器直接测量2,因此G P S观测技术成为当前监测S S E的有效手段3。G P S坐标时序分析中容易忽略S S E信号的影响4-5,但S S E信号在构造运动强烈的俯冲带上非常明显。对G P S坐标时序中S S E信号进行准确建模和研究,可以得到S S E发生的时间、大小和位置分布6,有助于进一步了解断层系统及其动力学发展规
6、律,为地震机理研究开辟新途径7。S S E信号属于G P S坐标时序中的瞬态信号,目前已有许多方法应用于G P S坐标时序中瞬态信号的识别与建模。例如利用区域网滤波方法可以过滤出S S E信号,但需要给定G P S坐标时序模型参数或S S E坐标时序,且无法准确估计S S E发生的时间8;还有采用纯粹的数学方法进行检测,如正交函数分解方法9、频谱分析方法1 0等。但上述方法无法对S S E信号进行模型化处理,难以用具体参数解释S S E发生的时间和区域范围。R o u s s e t等1 1首先利用三角函数描述和确定单位振幅滑移的时间,然后结合地球物理模型估计滑移区域;L a r s o n等
7、1 2首先利用双曲正切函数描述S S E的总滑移量和时间参数,然后结合其他地球物理参数(如断层格网的搜索宽度、长度、中心位置等)反演滑移区域范围、深度等信息。但时间参数采用的是概略值,带有一定的主观性和模糊性,因此有必要构建一个明确的S S E信号模型,进一步完善G P S坐标时序模型。本文在自主研发G P S坐标时序分析G T S A软件的基础上,引入双曲正切函数对S S E信号建模。利用该 软 件 批 量 处 理 日 本 房 总 半 岛 多 个G P S站点的坐标时序,分离出S S E信号,并利用S S E信号模型估计S S E发生的时间、地表位移及影响区域。1 G P S坐标时序中S S
8、 E信号模型由于双曲正切函数的几何形态与S S E快速滑移的瞬时运动特征较为相符,因此本文引入双曲正切函数对G P S坐标时序中的S S E信号进行建模。S S E信号模型如下:S(t)=ni=1ui2t a n ht-T0ii -1 (1)大 地 测 量 与 地 球 动 力 学2 0 2 3年9月式中,S(t)为t时刻的S S E坐标时序,n为S S E发生的次数,ui为第i个S S E的地表位移,T0i为第i个S S E发生时间的中心,i为第i个S S E持续的时长。考虑S S E信号后G P S坐标时序模型表述为1 3:x(t)=x0+v t+nj=1bjH(t-tj)+a1l n1+t
9、-te q +mi=1Bis i n(it)+Cic o s(it)+S(t)+(t)(2)式中,x(t)为t时刻测站的G P S原始坐标时序;x0为常数项;v为以年为时间的速度项;bjH(t-tj)为H e a v i s i d e阶跃项,bj为对应的阶跃次数,tj为发生阶跃的时间,阶跃项一般由GN S S测站周边 的 地 震 或 接 收 机 天 线 等 变 更 引 起;a1l n1+t-te q 为震后弛豫项,a1为弛豫项振幅系数,为震后的弛豫时间,te q为大震发生的时刻;Bis i n(it)+Cic o s(it)为以年和半年为单位的周期项,Bi、Ci为周期项振幅,i为年和半年的周
10、期项角速度;(t)为噪声。在考虑S S E信号的G P S坐标时序模型中,G P S原始坐标时序及其对应时刻为已知量,其余各项需通过数据预处理及模型拟合估计方法逐步求得。首先去除阶跃项,然后依次提取和去除常数项、稳态速度项和震后弛豫项,最后进行S S E信号建模。S S E信号建模时需估计3个未知参数:S S E的地表位移、发生时间的中心和持续时长。目前,全球多个俯冲带区域内出现S S E现象:新西兰北岛希库朗吉俯冲带上发现不同特征的S S E1 4;墨西哥俯冲带记录多次S S E及非火山震 颤1 5;日 本 西 南 部 南 开 俯 冲 带 多 次 发 生S S E1 6,沿该俯冲带的日本房总
11、半岛在菲律宾板块断层面曾发生一系列周期性S S E1 7。以房总半岛最近3次S S E为例,对该区域G P S观测的S S E信号进行建模,并围绕S S E模型中的时间和位移参数展开分析。2 S S E信号提取与建模2.1 S S E信号提取自2 0 1 1年日本9.0级地震后,日本房总半岛分别在2 0 1 1年、2 0 1 32 0 1 4年、2 0 1 8年发生3次S S E。研究选取该区域内2 0 0 92 0 1 9年以d为单位的G P S高精度坐标时序,该坐标时序主要由内华达大学采用G I P S Y/O A S I S软件解算得出,地球参考框架为I G S 1 4,E、N方向上的解
12、算精度在1mm以内,U方向上的解算精度在3mm以内,具体解算策略可参考h t t p:g e o d e s y.u n r.e d u/。以测站J 2 2 6的E方向为例,利用自主研发的G T S A软件进行S S E信号的分段提取,第1段为2 0 0 9 2 0 1 1年、第2段为2 0 1 2 2 0 1 3年、第3段为2 0 1 5 2 0 1 7年,结果如图1所示。由图可见,季节性周期项影响在mm级,对S S E信号的影响可忽略不计。因此,数据处理后的G P S坐标时序主要包含S S E坐标时序信号,可用于S S E建模。图1 3次S S E信号提取F i g.1 T h es i
13、g n a l e x t r a c t i o no f t h r e eS S E2.2 S S E信号建模对上述提取的S S E信号进行分段建模,以测站J 2 2 6为例,E、N、U三个方向上的模型拟合效果如图2所示。由图可见,S S E的时间范围约为12个月,且每次S S E发生的时段长短不同,表现为缓慢加速-快速滑动-缓慢减速,最后恢复至匀速稳定状态。受2 0 1 1-0 3-1 1日本宫城MW9.0地震影响,去除震前稳态速度后仍存在速度项,且S S E发生前后每段G P S坐标时序的速度也有差异。因此对G P S坐标时序进行整体建模时,还需考虑速度变化的影响(图3)。图2、3中
14、利用双曲正切函数模型拟合S S E信号的拟合效果最为突出,说明双曲正切函数的数学特性与S S E信号的运动规律相吻合。3 S S E模型的时间与位移参数3.1 时间和位移参数大小S S E的双曲正切函数模型中包含的3个未知参数分别为发生时间的中心T0i、持续时长i和019 第4 3卷第9期李 萌等:G P S监测慢滑移信号的双曲正切函数建模图2 2 0 1 1年、2 0 1 32 0 1 4年、2 0 1 8年S S E信号及建模F i g.2 S S Es i g n a l sa n dm o d e l s i n2 0 1 1,2 0 1 3-2 0 1 4a n d2 0 1 8图3
15、 G P S坐标时序中S S E信号建模F i g.3 S S Es i g n a lm o d e l i n g i nG P Sc o o r d i n a t e t i m es e r i e s地表位移ui。根据上述坐标时序的S S E信号,首先可初步确定S S E发生的时间范围;然后在该时间范围内逐步变动T0i作为已知参数输入,采用最小二乘法估计持续时长和地表位移,并计算模型的拟合优度1 8;最后根据最佳拟合优度确定最优估值T0i、i、ui。经计算发现,房总半岛1 6个测站3次S S E发生时间的中心T0i与持续时长i均存在差异。根据时间中心和持续时长计算测站S S E发生
16、的时间区间为T0i-i/2,T0i+i/2,从1 6个测站S S E模型的时间参数中分别选取起始时间的最小值和终结时间的最大值作为S S E的时间区域。计算结果发现,2 0 1 1年的S S E为5 0d、2 0 1 32 0 1 4年的S S E为2 6d、2 0 1 8年的S S E为3 9d。相比于文献1 7 中的参考时间范围,2 0 1 1年的S S E结束时间提前1 2d,2 0 1 32 0 1 4年的结束时间提前1 0d,2 0 1 8年的开始时间提前1 5d、结束时间提前2 8d;相比于文献1 9 的时间范围,S S E模型的总体时间参数范围缩短了1 01 3d。原因主要是文献研究采用的S S E时间范围为概略参考时间1 7,1 9,包含了无法准确确定的S S E初始缓慢加速和最后缓慢减速阶段的时间范围;而模型化的S S E从测站的几何形态出发,主要包含每个测站快速滑动的时间区间1 9,更能精确凸显S S E的时间特性。计算得出的房总半岛3次S S E的位移参数如表1(单位c m)所示。由表可见,2 0 1 8年的S S E最大水平位移可达5.2 4c m、2 0 1