留言板

尊敬的读者、作者、审稿人, 关于本刊的投稿、审稿、编辑和出版的任何问题, 您可以本页添加留言。我们将尽快给您答复。谢谢您的支持!

姓名
邮箱
手机号码
标题
留言内容
验证码

基于GNSS观测网络的断层滑移时空反演

徐克科 伍吉仓 雷伟伟

徐克科, 伍吉仓, 雷伟伟. 基于GNSS观测网络的断层滑移时空反演[J]. 武汉大学学报 ● 信息科学版, 2016, 41(10): 1391-1397. doi: 10.13203/j.whugis20140400
引用本文: 徐克科, 伍吉仓, 雷伟伟. 基于GNSS观测网络的断层滑移时空反演[J]. 武汉大学学报 ● 信息科学版, 2016, 41(10): 1391-1397. doi: 10.13203/j.whugis20140400
XU Keke, WU Jicang, LEI Weiwei. Spatio-Temporal Inversion of Fault Slip Based on GNSS Observation Network[J]. Geomatics and Information Science of Wuhan University, 2016, 41(10): 1391-1397. doi: 10.13203/j.whugis20140400
Citation: XU Keke, WU Jicang, LEI Weiwei. Spatio-Temporal Inversion of Fault Slip Based on GNSS Observation Network[J]. Geomatics and Information Science of Wuhan University, 2016, 41(10): 1391-1397. doi: 10.13203/j.whugis20140400

基于GNSS观测网络的断层滑移时空反演

doi: 10.13203/j.whugis20140400
基金项目: 

国家重点基础研究发展计划 No. 2013CB733304

国家自然科学基金 No. 41404023

详细信息
    作者简介:

    徐克科,博士,主要从事GNSS数据处理与地壳形变分析研究。12xkk@tongji.edu.cn

    通讯作者: 伍吉仓,博士,教授。jcwu@tongji.edu.cn
  • 中图分类号: P228.41;P542

Spatio-Temporal Inversion of Fault Slip Based on GNSS Observation Network

Funds: 

The National Program on Key Basic Reseach Project of China No. 2013CB733304

National Natural Science Foundation of China No. 41404023

More Information
图(8) / 表(1)
计量
  • 文章访问数:  2847
  • HTML全文浏览量:  62
  • PDF下载量:  286
  • 被引次数: 0
出版历程
  • 收稿日期:  2014-10-10
  • 刊出日期:  2016-10-05

基于GNSS观测网络的断层滑移时空反演

doi: 10.13203/j.whugis20140400
    基金项目:

    国家重点基础研究发展计划 No. 2013CB733304

    国家自然科学基金 No. 41404023

    作者简介:

    徐克科,博士,主要从事GNSS数据处理与地壳形变分析研究。12xkk@tongji.edu.cn

    通讯作者: 伍吉仓,博士,教授。jcwu@tongji.edu.cn
  • 中图分类号: P228.41;P542

摘要: 利用GNSS网络位移时空序列,基于弹性位错理论,构建了断层滑移时空分布的动态卡尔曼滤波反演模型。考虑断层面的非均匀滑动,将断层面细分为多个子断层,获取了较精细的滑移空间分布,并顾及了先验信息和拉普拉斯平滑约束。鉴于断层滑移引起的地表形变具有高空间相关性的特点,利用整个GNSS观测网络数据一起参与反演,有效分离了空间不相关的噪声。实验结果表明,当断层形变位移量和噪声水平相当,且点位分布间隔沿走向和倾向至少与子断层长宽等同时,均能反演得到正确的断层滑移时空分布。若信噪比不变,测站分布密度继续增大时,对反演效果提高并不显著,但能够容忍较低信噪比的观测数据。

English Abstract

徐克科, 伍吉仓, 雷伟伟. 基于GNSS观测网络的断层滑移时空反演[J]. 武汉大学学报 ● 信息科学版, 2016, 41(10): 1391-1397. doi: 10.13203/j.whugis20140400
引用本文: 徐克科, 伍吉仓, 雷伟伟. 基于GNSS观测网络的断层滑移时空反演[J]. 武汉大学学报 ● 信息科学版, 2016, 41(10): 1391-1397. doi: 10.13203/j.whugis20140400
XU Keke, WU Jicang, LEI Weiwei. Spatio-Temporal Inversion of Fault Slip Based on GNSS Observation Network[J]. Geomatics and Information Science of Wuhan University, 2016, 41(10): 1391-1397. doi: 10.13203/j.whugis20140400
Citation: XU Keke, WU Jicang, LEI Weiwei. Spatio-Temporal Inversion of Fault Slip Based on GNSS Observation Network[J]. Geomatics and Information Science of Wuhan University, 2016, 41(10): 1391-1397. doi: 10.13203/j.whugis20140400
  • 地震孕震过程是一个长期复杂的缓慢过程。而地震发生时断层破裂所释放的能量只是其中一部分,很大一部分能量在常规地震前或后,以断层无震蠕滑的形式释放。慢地震和无震蠕滑是伴随活动断层地震应力成核的重要过程,每一次断层蠕滑都可能会转移部分应力到其上部的锁定区域,使其应力承载在滑移时刻的增加大大高于平时的增加,地震危险性也会大大高于平时[1-2]。因此,动态反演断层滑移时空分布特征及其演变过程,正确确定活动断层的初始应力场和边界条件以及判断断层锁定、蠕滑转换的周期性滑移过程,均对地震危险性评估至关重要。

    连续运行的GNSS网络提供了时空采样分辨率越来越高的数据,在断层形变检测和反演方面显示出了巨大优势。Segall最早提出了GNSS网滤波反演的思想[3]。McGuire研究了GNSS网络扩展滤波反演,并成功获取1999年加拿大Cascadia慢地震滑移分布[4]。王武星和方颖等人对GNSS网滤波方法也给予了初步探讨,并进行了相应试验[5-6]。Kositsky等人提出了主成分反演方法,并获取了2005年8.6级Nias地震震后余滑分布[7]。Radiguet等人基于线性位错模型用15个GPS台站反演了2006年墨西哥慢地震滑移时空分布[8]。陈光齐等人利用GNSS分析了震前应变积累及其时空特征和同震变形特征的差异[9]。丁开华等人利用GNSS时序分析,采用单纯形法求得了汶川地震震后松弛时间约为38 d,并基于粘弹性模型反演了流变学参数[10]。目前,利用GNSS数据不仅仅局限于计算地壳平均运动速度和反演同震滑移、震间复位错等静态模式方面,诸多学者开始致力于震间和震后断层微动态形变异常检测和断层滑移动态反演等方面的研究[11-12]。如何构建稳健的断层滑移时空分布反演模型,至少需要多大的信噪比和GNSS台站分布密度才能正确反演断层滑移时空分布等问题需要进一步研究和考证。

    本文基于OKADA弹性位错理论,利用断层分布区域的GNSS网络时空序列,顾及先验信息和平滑约束,构建了GNSS网络卡尔曼滤波断层滑移时空反演模型。通过大量模拟实验,分析了不同信噪比和不同台站分布密度情况下的反演效果,得出了正确反演断层滑移时空分布所需要的最低信噪比和最优的台站分布密度。

    • GNSS台站位移时间序列可表示为:

      (1)

      式中,t为时间;a0为初始位置;b为线性趋势项;第三项为阶跃项,用来修正地震或天线改变在时间序列所造成的落差,可利用差分方法确定阶跃发生的时刻进行消除;第四项为周年和半周年周期项;第五项代表断层滑移所造成的地表位移,其中G表示OKADA弹性位错模型格林函数[13]s为断层滑移量;δcme为空间共模误差,可结合PCA(principal component analysis)和KLE(Karhunen-Loeve expansion)方法去除[14]εgnss为观测误差,为避免解算时某些误差未被模型化而导致误差低估,需对其协方差阵Σgnss乘以一比例因子σ2,即变为σ2Σgnss

    • 细分断层时,未知数个数将远大于观测值个数。此时,附加不同约束条件,保持解的稳定至关重要。为避免滑动分布解的振荡,避免相邻子断层滑动量在大小和方向上存在显著差异,通常采用拉普拉斯平滑约束,设置相邻断层间滑动量的梯度为最小[15],拉普拉斯二阶差分算子为:

      (2)

      式中,s(i,j)表示位于第i行、第j列的子断层上的滑动量;Δx、Δy分别表示相邻子断层沿走向和倾向的距离。按式(2) 分别对所有断层单元求相应的二阶差分算子,观测方程的形式为:

      (3)

      式中,左边表示虚拟观测值,一般假定为零;Hn为所有子断层滑移拉普拉斯二阶差分算子;s为所有子断层走向和倾向滑移大小;γ为描述断层面上的滑移在空间上的平滑程度,值越小越平滑,反之滑移变化大;I为单位阵。

      通常,通过地球物理勘探可探知断层的某些先验信息,如滑移速度、最大滑移量和滑移特征等,可采用顾及被估参数的先验信息约束条件。设s0为待估参数s的先验值,ε+为参数先验误差。观测方程的形式为:

      (4)

      式中,β用于制约断层滑移大小,使滑移及滑移率始终控制在一定范围内变动。

    • 根据位置和速度系统噪声的谱密度矩阵[16],可以得到滑移和滑移率的状态方程,其转移矩阵为:

      (5)

      过程噪声矩阵为:

      (6)

      式中,α为速度谱密度,用来控制滑移和滑移率在时间域的平滑程度。

      在上述观测方程及卡尔曼滤波状态模型中,共有4个超参数需要估计,分别为σ、γ、βα,可采用极大似然方法进行估计[3]。因极大似然估计比较耗时,本文把超参数转换至观测方程中,利用扩展卡尔曼滤波方法随同模型参数一并求解[4]。转换后,观测方程(1) 、(3) 、(4) 依次变为:

      (9)

      设第k历元总的观测值为dk,则包括GNSSN、E、U向位移dgnss和约束条件虚拟观测值d、d+。设待估参数为${\hat{x}}$,则包括每个历元的断层滑移s和滑移率、GNSS位移时序中的初值a0、线性项b、周期项系数a1、a2、b1、b2和超参数σ、γ、β、α。设观测值误差项为ε,包括εgnss、ε∇、ε+,其协方差阵可表示为:

      (10)

      则总观测方程可表示为:

      (11)

      在整个时间序列中,设待估参数中除断层滑移s和滑移率${\hat{s}}$外,其它10个参数均为恒定参数,不随时间改变。因此,其转移矩阵为单位阵I,过程噪声矩阵为零,则整个待估参数总转移矩阵和过程噪声矩阵分别为:

      (12)

      正向状态预测向量为:

      (13)

      相应协方差矩阵为:

      (14)

      k历元新息向量可表示为:

      (15)

      其协方差阵为:

      (16)

      根据卡尔曼滤波递推公式,参数估值及其协方差阵可求解得:

      (17)
      (18)
      (19)

      根据观测值的预报残差向量来进行卡尔曼滤波反演结果的精度评定,求其单位权方差为:

      (20)

      其中,n为观测值个数,则第k历元参数估值的协方差阵为:

      (21)

      同样,回代状态预测向量为:

      (22)

      相应协方差矩阵为:

      (23)

      回代卡尔曼滤波参数估计与精度评定同正向。

    • 目标断层按图 1所示,沿走向和倾向划分为2.5 km×1 km大小的子断层,共8×5个方格,每一方格代表一子断层在平面的投影,星型符号表示测站。断层模型的具体参数如表 1所示。设置断层为非均匀滑动,滑移特征为逆冲兼走滑,模拟300 d呈指数函数变化的滑移时间序列。其中,每隔30 d的滑移分布见图 2。设置观测台站均匀分布在断层区域,由OKADA弹性位错模型正演300 d地表站点N、E、U向的位移时间序列。在不同密度的测站分布和不同信噪比的情况下进行反演分析。

      图  1  断层划分与测站分布

      Figure 1.  Fault Division and Distribution of Observation Stations

      表 1  断层参数

      Table 1.  Fault Parameters

      下边缘中心位置(N,E)/km 方位角/(°) 长/km 宽/km 下边缘深/km 倾角/(°)
      (50,50) 90 40 5 5 70

      图  2  每隔30天的实际滑移分布

      Figure 2.  Actual Slip Distribution Every 30 Days

      试验1 设置观测站点沿走向和倾向分布间隔约为2倍子断层长和宽的距离,观测站数为5×4个,分布如图 1中的红色星号。根据噪声功率谱密度,由Fakenet模拟软件[17]加入和实际位移大小相当的噪声,白噪声和有色噪声水平之比约为1:1,合成模拟的地表位移序列见图 3

      图  3  模拟噪声与合成地表位移序列(试验1)

      Figure 3.  Simulation Noises and Synthetic Displacement Series(Test 1)

      利用表 1的断层参数,由§1构建的GNSS网卡尔曼滤波反演模型进行断层滑移时空分布反演。反演的所有子断层每隔30 d滑移分布见图 4。在时间分布上,反演的最大滑移量为350 mm,而实际为400 mm。从时间分布上可以得出与实际断层滑移一致的演变特征。在空间分布上,滑移量最大的位置位于东向52 km、北向50 km,实际滑移量最大的地方为东向50 km、北向50 km;整体空间分布与实际断层滑移分布较一致。若噪声水平不变,减少站点分布密度,反演效果较差。若点位分布密度保持不变,放大噪声倍数,反演效果较差。由此得出,当断层滑移贡献位移与噪声水平相当时,测站分布密度沿走向和倾向至少为2倍子断层长和宽时,顾及平滑约束和先验信息约束的网络卡尔曼滤波反演能够得到正确的断层滑移时空分布。

      图  4  试验1中每隔30天的反演滑移分布

      Figure 4.  Distribution of Inversion Slip,Every 30 Days,Test 1

      试2二 设置噪声水平同试验1,增加观测站点的分布密度为原来的2倍,即沿走向和倾向的分布间隔约等同子断层长和宽,见图 5。反演的每隔30 d滑移分布见图 6。由图 6可知,在空间分布上,整体比模拟试验1略有提高;时间分布同模拟试验1,反演所有子断层滑移量总体提高不明显。由此得出,当断层滑移贡献位移与噪声水平相当时,测站分布密度沿走向和倾向为2倍子断层长和宽时,已能得到正确的断层滑移时空分布,若增加测站分布密度,提高并不显著。

      图  5  测站点位分布

      Figure 5.  Distribution of Observation Stations

      图  6  试验2中每隔30天的反演滑移分布

      Figure 6.  Distribution of Inversion Slip,Every 30 Days,Test 2

      试验3 点位分布密度同试验2,噪声放大一倍;其合成地表位移时间序列如图 7所示;反演的每隔30 d断层滑移分布见图 8。从反演结果看,仍能正确得到断层滑移时空分布及其演化特征。由此得出,当测站分布密度增大时,能够容忍GNSS位移序列较低信噪比(SNR<1) ,有效提高反演信噪比。

      图  7  模拟噪声与合成地表位移序列(试验3)

      Figure 7.  Simulation Noises and Displacement Series(Test 3)

      图  8  试验3中每隔30天的反演滑移分布

      Figure 8.  Distribution of Inversion Slip,Every 30 Days,Test 3

    • 本文基于GNSS网络位移时空序列构建了顾及平滑约束和先验信息约束的卡尔曼滤波反演模型。试验验证了通过滤波的方法能够有效分离空间不相关的噪声和空间高相关的断层形变,准确估计断层滑移时空分布特征。通过试验得出,当断层滑移贡献位移与噪声水平相当时,测站分布密度沿走向和倾向为2倍子断层长和宽时,已能反映出正确的断层滑移时空分布,若增加测站分布密度,反演结果提高并不显著。但增加测站分布密度,能有效提高反演信噪比,甚至适用于信噪比低于1的情况。由此,结合GNSS位移解算精度和所要检测的断层滑移大小与分布,可选择最优的GNSS站点布设方案。或根据现有的GNSS站点分布,选择最优的子断层划分。

      在反演前需要对GNSS数据进行预处理,如剔除粗差、补齐缺失数据、修正阶跃项、去除共模误差等。当观测方程较少或参数之间相关时,为保持解的稳定性,要尽量减少参数个数。本文反演模型中没有考虑随机游动和GNSS坐标参考框架误差。在保持参数求解稳定的情况下,改正更多的误差项则反演效果将会更好。因为唯有尽可能去除各种相关噪声,才有可能检测慢地震或无震蠕滑等更微弱的断层形变时空分布信息,获取更为准确的断层滑移特征及演变过程。

参考文献 (17)

目录

    /

    返回文章
    返回