地球物理学报  2019, Vol. 62 Issue (11): 4142-4155   PDF    
利用独立成分分析检测2004年和2012年印度洋地震的重力变化
王陈燕, 游为, 范东明     
西南交通大学地球科学与环境工程学院, 成都 611756
摘要:2004年12月苏门答腊发生MW9.3地震,造成巨大的质量重新分布.利用GRACE卫星月重力场数据计算了研究区域地面1°×1°网格点上的重力变化时间序列,采用主成分分析和独立成分分析两种方法,提取了重力变化的空间与时间特征,结果显示震中两侧区域的重力变化呈两极分布,其中东侧重力下降,西侧重力增加.相较于传统的主成分分析方法,独立成分分析能更好地从原始信号中提取地震的信息,能分解出具有显著阶跃变化的独立成分.除了2004年的苏门答腊大地震外,独立成分分析还分解得到了2012年Andaman地震的特征,与该次地震的空间特征与时间序列基本一致.相较于适合定量分析的多项式拟合方法,独立成分分析更适合大范围区域的定性分析,建议将两种方法相结合,取长补短,从而为GRACE地震监测提供一种更为客观、有效的方法.
关键词: GRACE      苏门答腊大地震      主成分分析      独立成分分析     
Using independent component analysis to detect the gravity changes in the Indian Ocean Earthquakes of 2004 and 2012
WANG ChenYan, YOU Wei, FAN DongMing     
Faculty of Geosciences and Environmental Engineering, Southwest Jiaotong University, Chengdu 611756, China
Abstract: A massive redistribution of mass was caused by the MW9.3 earthquake that occurred in Sumatra in December 2004. By using the gravity data of GRACE satellite, the time series of gravity change on the ground 1°×1°grid points are calculated, and the spatial and temporal characteristics of gravity change are extracted by using Principal Component Analysis and Independent Component Analysis, and the results show that the gravitational change in both sides of the epicenter has the character of polar distribution, the gravity of the up disk in the fault decreases, while the gravity of the down disk increases. Compared with the principal component analysis method, independent component analysis can better extract seismic-related information from the GRACE gravity signal, and can decompose independent components with significant gravity changes. In addition to the 2004 earthquake in Sumatra, the independent component analysis also decomposed the characteristics of the 2012 Andaman earthquake, which is basically consistent with the spatial characteristics and time series of the earthquake. Compared with the polynomial fitting method which is suitable for quantitative analysis, independent component analysis is more suitable for qualitative analysis of large-area regions. The two methods can be combined to provide a more objective and effective method for GRACE seismic monitoring.
Keywords: GRACE    Sumatra earthquake    Principal component analysis    Independent component analysis    
0 引言

Gravity Recovery and Climate Experiment(GRACE)卫星任务由美国宇航局(NASA)和德国空间局(DLR)提出计划并联合研制,通过测量两颗卫星之间的距离变化,恢复时变地球重力场模型,进而推导出地球表面及浅层的质量变化(郑伟等, 2009, 2017; 严畅达和徐亚, 2018),GRACE时变地球重力场的应用研究在地学领域已经取得了不少重大科研成果(Tapley et al., 2004; Han et al., 2006; Rodell et al., 2009; Chen et al., 2009; Matsuo and Heki, 2010; Luo et al., 2012).在众多的地球物理现象中,地震能够造成强烈的地表重力变化,许多学者都从理论上论证了GRACE检测大型地震的能力(Mikhailov et al., 2004; Sun and Okubo, 2004),多次特大型地震也验证了这一点.Han等(2006)首次利用GRACE卫星重力场检测到了2004年苏门答腊大地震引起的重力变化,2010年的智利大地震以及2011年的东日本大地震所造成的重力变化也被GRACE卫星检测到(Han et al., 2010; Heki and Matsuo, 2010; 周新等, 2011; Matsuo and Heki, 2011; Wang et al., 2012; Zhou et al., 2012),以上诸多有关GRACE检测地震的研究验证了GRACE能探测矩震级MW约8.5级以上地震重力变化的能力.

在GRACE的质量变化信号中,相比于陆地水,地震信号本身相对比较弱,为了更好地从中提取所需要的信号,需要用一定的技术手段进行信号分析与增强(严畅达和徐亚,2018),如时间序列的多项式拟合等方法.Li等(2014)系统介绍了利用GRACE数据探测地震的数据处理方法及其局限性,其研究表明,利用多项式拟合等信号分析方法来提取或分离共震和震后信号,可能会导致信号的衰减或失真.而主成分分析(Principal Component Analysis, PCA)可以有效地从这类时空数据中提取出信号的时间序列与空间分布(Preisendorfer, 1988),研究表明,PCA可以提取出60%的8.8级地震引起的同震重力变化(De Viron et al., 2008).周江存等(2013)利用PCA方法计算出GRACE卫星重力结果中的同震重力变化,张克亮等(2014)利用GRACE数据,结合PCA方法分析了2011年发生在日本东部海域的大地震.虽然PCA方法在信号分析的角度具有一定的优势,但其缺陷在于提取出来的主分量未必具有明确的物理含义(邹正波等, 2016),且主成分之间只是彼此不相关,而非相互独立.独立成分分析(Independent Component Analysis, ICA)可以被视为是PCA分解的一种拓展延伸,用以对所得到的成分进行改善,使其更加合理(Richman, 1986),从现有文献来看,目前还没有结合ICA方法提取GRACE数据中地震信号的研究.

本文以2004年12月的苏门答腊大地震为例,利用传统的方法对苏门答腊地震的重力变化进行检测,即在变形较大的主震附近选取有代表性的点来研究(Ogawa and Heki, 2007; Chen et al., 2007; Panet et al., 2007),作为PCA和ICA分解结果的参照物.利用PCA和ICA对GRACE重力变化数据进行处理,将分解得到的各个成分与主震附近有代表性的点进行比较,比较其空间特征和时间序列,可以发现,在PCA的前12个主成分中,有两个与2004年的地震有关,但并没有将地震的两个阶跃完全分离开.将PCA的前4个模式进行旋转,实现ICA算法,其对原始信号的分离效果更好,对于某些地球物理实事件的解释更加合理,结果表明ICA比PCA更适合提取代表地震的信号.但PCA和ICA的分解结果表明,2012年Andaman地震对研究区域的影响也比较大,将PCA分解的前12个主成分进行旋转,由此得到的独立成分中可以发现2012年Andaman地震的信号.此外,缩短参与ICA的数据时间段长度以去除掉2004年地震的影响,同时增加参与旋转的模式数量,得到了与2012年地震4个特征点一一对应的独立成分,其空间特征与时间序列都吻合得比较好.以上结果表明,ICA方法可以从GRACE地面重力扰动数据中分离得到2004年和2012年两次地震的信号,得到的独立成分与多项式拟合方法的结果相互验证,相较于常用的多项式拟合方法,ICA方法能更好地对地球物理事件的整体空间特征进行分析.相较于PCA,ICA的结果更加客观,分解得到的结果也更加合理,与实际情况更加符合,ICA可以在不需要先验信息的前提下获得研究区域重力变化的主要空间和时间特征,对于大范围研究区域的定性分析具有一定的优势.在实际的研究过程中,可以与适合定量分析的多项式拟合方法相结合,取长补短.

1 基本原理 1.1 GRACE时变地球重力场计算

利用GRACE数据处理机构GFZ发布的时变地球重力场模型计算地面重力扰动变化的公式为(Wahr et al., 1998; 周江存等, 2009; Zhang et al., 2009)

(1)

式中,GM为万有引力常数和地球质量的乘积,R为地球平均半径,lm分别表示球谐位系数对应的阶和次,Plm(cosθ)是完全规格化的勒让德(Legendre)函数,ΔClm与ΔSlm是每月球谐位系数的变化值,θ是地心余纬,λ是地心经度,Wl, Wm为滤波系数.由于球谐位系数高阶项存在较大误差,阶数大于60的球谐位系数对分析质量变化作用不显著(卢飞, 2015),因此在计算过程中截断阶次的取值为60.值得注意的是仅采用一种滤波并不能完全削弱条带误差的影响,需结合其他方法进行,本文采用Chen等(2007)提出的P3M6的去相关滤波方法,并结合250 km的扇形滤波进行组合滤波.由于GRACE卫星能探测到各种原因引起的地表质量变化,其中最主要的是陆地水,故若将GRACE仅用于探测重大地震的重力变化,需要扣除陆地水对重力变化信号的影响.目前常用的方法有两种,一个是利用陆地水变化的周期性,即利用多项式拟合或年度变化均值来去除,第二种是利用水文模型扣除掉陆地水的影响.两种方法都有其局限性,对于第一种特别是多项式拟合方法,需要了解研究区域的陆地水变化特征,即需要一定的先验信息,而对于第二种方法,不同的水文模型由不同的数据和算法计算得到,引入的水文模型也会带来一定的误差.除了这两种方法,本文提到的ICA方法对于削弱陆地水的影响也有一定的作用,该方法不需要任何先验信息,且不会引入其他的不确定量,可以利用陆地水信号和重力变化信号之间不同的统计信息进行分离,特别是对于发生在陆地上的强震,因为相较于发生在海底的强震,陆地受到陆地水的影响更大,其周期性更明显,对于ICA而言,陆地水与重力变化时间序列的统计特性区别更大,区分效果更好.

1.2 ICA的原理与算法

ICA可被视为PCA的一种扩展.PCA是一种基于二阶统计信息的信号分解方法(Hyvärinen, 1999),可以在减少数据集维数的同时,保留数据集的主要特征.设每个历元含有p个观测点,n个历元的数据可以存储在一个n×p的矩阵X中,则

(2)

式中x1, x2, …, xp均表示n×1维列向量,其中的每一个元素都减去了所对应列向量的平均值,矩阵X的协方差矩阵可以表示为

(3)

可以将矩阵X的PCA分解表示为

(4)

式中,EX的特征向量组成,即经验正交函数,表示空间分布,可以由协方差矩阵C的特征值分解计算得到.P=XE,储存着主成分,代表着时间序列,可以将P视为原始数据X在正交基函数E上的投影.为降低数据集的维度,可以将主要模式保留下来,一般保留前f项主要模式,再对其进行重构,即

(5)

式中PfEf分别包含前f项主成分和经验正交函数,ICA就是将上述PCA分解得到的主成分PC或经验正交函数EOF进行旋转变换,一般而言,如果ICA算法是对PC进行旋转,那么不同成分的时间序列将会尽可能地独立,这样的ICA被称为Temporal ICA(TICA),与此对应,如果ICA算法是对EOF进行旋转,即可称为Spatial ICA(SICA)(Forootan and Kusche, 2012).利用TICA或SICA,可以使得PCA分解的结果在时间或空间上更加合理(Richman, 1986),即

(6)

式中Rf表示正交旋转矩阵,即RfRfT=I.一般由Cardoso和Souloumiac(1993)提出算法计算得到旋转矩阵,该算法又被称为特征矩阵的联合近似对角化(Joint Approximation Diagonalization of Eigen-matrices, JADE).对于PCA和ICA所提取出来的模式,其对应时间序列与空间分布的乘积才是真正的信号,单一地分析空间分布或时间序列是不全面的.

2 GRACE卫星检测到的苏门答腊地震的重力变化

发生在2004年12月26日的苏门答腊大地震震中位置为(3.317°N, 95.854°E),震源深度28.6 km,地震矩为3.95×1022N·m(周江存等, 2013).以震中为中心、计算一定区域内1°×1°地面网格点上的重力扰动变化,每个网格点上数据起止时间为2003年1月—2014年12月,共144个月的结果,其中缺失了11个月的数据,分别为2003年6月,2011年1月和6月,2012年5月和10月,2013年3月、8月和9月以及2014年2月、7月和12月,对于缺失的月重力场数据,采用整体三次多项式插值的办法补齐.利用式(1)对整个区域(10°S—20°N, 85°E~115°E)内1°×1°网格点上的GRACE时间序列进行处理,得到研究区域内的地面重力扰动变化,根据地震前后两年平均重力场之差得到地震前后的重力场变化(Chen et al., 2007; 王武星等, 2011),如图 1所示,可以发现震中C点东西两侧产生明显的符号相反的重力变化,即二极分布.震中C点东侧马来半岛呈现大范围的重力负变化,A点(6.5°N, 97.5°E)的重力减小最为明显,而震中C点西侧印度洋海域都呈现一定的正变化,B点(1.5°N, 94.5°E)的重力增加最为明显.

图 1 2004年苏门答腊地震前后的重力场变化 A、B两点分别表示重力减少和增加最为明显的地区,C点表示震中.单位:μGal. Fig. 1 Gravity field changes caused by the 2004 Sumatra earthquake The two points A and B respectively indicate the area where gravity is reduced and increased most, and the point C indicates the epicenter. Unit: μGal.

A、B两点自2003年1月至2014年12月的GRACE重力变化时间序列如图 2中黑色曲线所示,图中时间序列标记的点即为地震发生时刻,显示出在地震前、后2个区域出现一次显著的突跳,分别对应重力的下降和增加,A、B两点在地震前均保持稳定的时间特征,同震时都有明显的重力阶跃,并且两点震后的变化趋势相同,均呈上升趋势,该结果与Chen等的结果一致(Chen et al., 2007; 张国庆等, 2015).由地震前后两个月的差值计算A和B两点的重力变化,分别为-13.7 μGal和5.2 μGal,时间序列中有明显的周期变化,表明该时间序列还受到陆地水等周期性变化项的影响,这样的结果作为地震同震重力变化显然是不合理的,对地震前后两段时间序列分别进行多项式拟合,扣除周年变化、半周年变化和161天周期变化项后的变化趋势如图 2中直线表示,由地震前、后两个月拟合值的差值确定出A和B两点的重力变化分别为-11.1 μGal和6.9 μGal.

图 2 观测点A和B的GRACE重力变化时间序列 直线表示扣除周年变化、半周年变化和161天周期变化项后的变化趋势,时间序列中黑点表示地震发生的时刻. Fig. 2 The time series of GRACE gravity change at observation points A and B The straight line indicates the trend after deducting the annual change, the half-anniversary change, and the 161-day cycle change. The black dot in the time series indicates the moment of the earthquake.
3 重力变化的PCA与ICA分解 3.1 重力变化的PCA分解

利用PCA方法对区域内的重力扰动变化时间序列进行处理,根据特征值的大小将主成分进行排序,计算的前12个特征值如图 3所示.可以看出,第1特征值占所有特征值之和的比值为40.9%,远大于其他的特征值,前4个特征值之和占比超过了80%,而前9个特征值之和占比则超过了90%.

图 3 主成分分析中不同模式的特征值频谱分布 两图分别表示不同模式所对应的特征值及特征值累计量所占的百分比,横轴表示对应的模式. Fig. 3 Spectral distribution of eigenvalues of different modes in principal component analysis

提取了前12个主成分,各个主成分空间特征和时间序列如图 4图 5.可以发现,前5个PC时间序列的周年变化都很明显,PC2和PC4都在2004年地震震中附近出现了与图 1类似的二级分布.结合PC1的空间分布和时间序列,可以发现,PC1的信号主要集中在中南半岛地区,其时间序列表现为单纯的周年变化,表现的是中南半岛的陆地水变化情况.在空间模式上,PC2在震中西侧重力正变化区域的信号强度强于东侧重力负变化区域,主要表现了地震中重力正变化区域的变化特征.其时间序列在2003年12月附近,出现了明显的突跳,图中黑色虚线即为归一化后的B点时间序列,可以看出,两个时间序列在趋势和变化幅度上比较接近,但两个在相位上有一定的区别,因为PC2同时反映了震中东西两侧重力变化的情况,结合图 2中A、B两点在时间序列相位上的区别,可以认为,PC2的时间序列是两种不同变化特征相互影响的结果.PC3在震中东北方向有明显的变化特征,且其对应的时间序列也出现了明显的突跳,将图 2中A点时间序列反号处理,并对其进行归一化,如图 5中PC3时间序列中的黑色虚线所示,其与PC3的时间序列相似,可以认为PC3主要表现了地震中重力负变化区域的变化特征.虽然PC4在震中附近出现了二级分布,但其变化最明显的区域集中在中国南海海域,而非苏门答腊一带,且PC4对应的时间序列并没有出现明显的突跳,故可以排除PC4代表地震变化的可能性.

图 4 PCA12个成分的空间模式, 图中黑点表示震中.单位:μGal Fig. 4 The spatial pattern of the 12 components obtained by the PCA and the lack dots in the picture indicate the epicenter. Unit: μGal
图 5 由PCA得到的12个主成分的时间序列 图中黑点表示地震发生的时刻,PC2和PC3中黑色虚线分别表示归一化后B点和A点的时间序列,其中A点时间序列进行了反号处理. Fig. 5 The time series of the 12 components obtained from PCA The black dots in the figure indicates the moment when the earthquake occurred. The black dotted lines in PC2 and PC3 represent the time series after normalization of point B and point A respectively, and the time series of point A is reversed.
3.2 ICA分解

图 1图 2可知,重力正负变化区域的时间序列是不同的,这是两种不同的信号,对于一个合理的信号分解方法,应该能将两种不同的信号完全区分开,无论是空间分布还是时间序列.而在3.1节中PCA的结果中,PC2、PC3和PC4都在重力负变化区域出现了变化的特征,这明显是不合理的,且其时间序列与图 2的两个时间序列相比,也只是比较符合,在相位和趋势上并没有完全一一对应.为了改善PCA分解的结果,对PCA的前四个空间模式进行旋转,即实现了ICA算法,得到的空间模式和时间序列如图 6图 7.

图 6 由ICA得到的4个独立成分的空间模式,图中黑点表示震中 Fig. 6 The spatial pattern of the 4 components obtained by the ICA, and the lack dots in the picture indicate the epicenter
图 7 由ICA得到的4个独立成分的时间序列 图中黑点表示地震发生的时刻,IC2和IC4中黑色虚线分别表示归一化后B点和A点的时间序列,其中A点时间序列进行了反号处理. Fig. 7 The time series of the 4 components obtained from ICA The black dots in the figure indicates the moment when the earthquake occurred. The black dotted lines in IC2 and IC4 represent the time series after normalization of point B and point A respectively, and the time series of point A is reversed.

从空间分布上看,IC1和IC3的信号分别集中在中南半岛和中国南海海域,而其时间序列都表现出明显的周年变化,说明IC1和IC3所代表的信息与地震无关.IC2和IC4的时间序列都出现了明显的突跳,且变化趋势很明显,周期性信号比较弱,将其分别与对应的空间特征相乘,得到这两个IC所代表的重力变化,IC2和IC4分别出现了正值和负值的跳变,与图 2中B、A两点的时间序列相吻合,如图 7中IC2和IC4的黑色虚线所示,而两者信号的空间分布分别位于重力正变化区域和负变化区域,没有出现PCA分解结果中PC2的二极分布.对比IC2和图 2中B点时间序列,发现IC2时间序列的周期项振幅更小些.将PCA与ICA两种方法与地震有关的空间分布进行比较,可以发现,ICA的分解结果更合理一些,将重力负变化区域和正变化区域的空间分布分解为两个独立成分.在其时间序列中,相较于PCA,ICA分解得到的IC2和IC4与图 2所示时间序列更加吻合,除了由地震导致的重力信号突变,周年变化信号要更弱一些,而周年变化信号明显不是由地震造成的.故可以得出结论,相较于PCA,ICA从GRACE数据中提取地震信息的能力要更强一些.

4 GRACE重力信号与2012年的Andaman地震 4.1 Andaman地震重力变化的时空特征

上文已经对比了PCA和ICA两种方法对2004年苏门答腊大地震的分解效果,对比两种方法的前四种模式,可以发现,PC3和IC4的空间模式在2004年震中的西南方向,都出现了相同的空间模式,即两个相对的区域出现了重力增加,而另外两个相对的区域出现了减少的情况,且其对应的时间序列都在2012年前后出现了重力的快速下降.这主要是由2012年4月Andaman地震所造成的(Han et al., 2015),震级为MW8.6和MW8.2,震中位置为(2.29°N, 93.08°E)(薛艳等, 2014).

参照上文的方法分析区域内的地面重力扰动变化,根据地震前后两年平均重力场之差得到地震前后的重力场变化,如图 8所示,研究区域内出现了非常明显的重力变化.震中D点东北侧和西南侧呈现大范围的重力正变化,其中两个区域的最大值分别出现在P2点(3.5°N, 96.5°E)和P3点(0.5°N, 88.5°E),而震中D点西北侧和东南侧都呈现一定的负变化,最大值出现在P1点(6.5°N, 88.5°E)和P4点(2.5°S, 95.5°E).图 9中第一个星型点表示2004年的苏门答腊大地震,而第二个五角星标记表示2012年的Andaman地震,时间序列显示在地震前后四个区域发生了显著的突跳,分别对应重力的下降和增加,该结果与Han等的结果一致(Han et al., 2015),也与上文中PC3和IC4的空间特征相对应,说明从GRACE重力数据中可以分解得到与Andaman地震有关的信息.

图 8 2012年Andaman地震前后的重力场变化 其中P2和P3均表示重力增加的区域,而P1和P4均代表重力减少的区域,D点表示地震的震中.单位:μGal. Fig. 8 Gravity field changes caused by the 2012 Andaman earthquake P2 and P3 both represent the area where gravity increases, while P1 and P4 both represent the area where gravity decreases, and point D represents the epicenter of the earthquake. Unit:μGal.
图 9 观测点P1、P2、P3和P4的GRACE重力变化时间序列 图中第一个黑点标记表示2004年地震发生的时刻,第二个黑色五角星标记表示2012年地震发生时刻.单位:μGal. Fig. 9 The time series of GRACE gravity change at observation points P1、P2、P3 and P4 The first black dots in the figure indicates the moment of the 2004 earthquake, and the second black five-pointed star marks indicates the moment of the 2012 earthquake. Unit: μGal.
4.2 重力变化的ICA分解

上文的ICA分解虽然能在IC4中观察到Andaman地震的变化特征,但在IC4的空间特征和时间序列中特征并不是很明显,而在PCA分解得到的PC5、PC6、PC7和PC9中,在该区域都出现了一定强度的信号,且部分重力时间序列在2012年4月附近出现了较为明显的变化.为探究Andaman地震的时空特征,仅利用PCA分解的前四个主成分模式进行旋转,是不合理的,应该利用更多的模式,因为在该研究区域中,2次地震震中位置比较接近,2004年的苏门答腊大地震所造成重力突变最为严重,Andaman地震的影响相对较小,利用较少数量的模式进行旋转,很容易将Andaman地震的信号隐藏在2004年的大地震之下,ICA分解得到的IC很可能就是两次地震的折中,因为用一个信号表现两个地球物理事件,只能体现它们之间共同的部分.

为此,将上文提到的12个模式进行旋转,得到的空间模式和时间序列如图 10图 11.图 10中黑点和白点分别表示2004年和2012年地震震中,而在图 11中两个标记分别表示两次地震发生的时刻.对比两次ICA分解得到的独立成分,两次分解都得到了2004年地震的时空特征,且都出现在比较靠前的位置,说明2004年地震的影响较大,且对应的时间序列都出现了明显的重力变化.而在第二次的ICA分解中,IC8和IC11表现出了2012年地震的特征,其空间分布都在震中附近出现了比较大的变化,而这两个时间序列都在地震发生前后出现了比较明显的重力突跳,其中IC8所对应的是P1所代表的重力减少地区,而IC11对应的是P3和P4所代表的地区,但没有找到P2点对应地区的独立成分,很有可能是代表该点的信号与2004年地震的信号相互影响,折中表示为IC3和IC4.

图 10 由ICA得到的12个独立成分的空间模式 图中黑点表示2004年地震震中,而白点表示2012年地震震中.单位:μGal. Fig. 10 The spatial pattern of the 12 components obtained by the ICA The black dots in the figure indicate the epicenter of the 2004 earthquake, while the white dots indicate the epicenter of the 2012 earthquake. Unit: μGal.
图 11 由ICA得到的12个独立成分的时间序列 图中第一个黑点表示2004年地震发生的时刻,第二个五角星标记表示2012年地震发生时刻,IC3和IC4中黑色曲线分别表示归一化后B点和A点的时间序列,其中A点时间序列进行了反号处理. Fig. 11 The time series of the 12 components obtained from ICA The first black dots in the figure indicates the moment of the 2004 earthquake, and the second black five-pointed stars indicates the moment of the 2012 earthquake. The black curves in IC3 and IC4 represent the time series after normalization of point B and point A respectively, and the time series of point A is reversed.

相较于只旋转4个模式的第一次ICA,第二次ICA分解的效果更好,不仅分解得到了与2004年苏门答腊大地震有关的独立成分,还得到了象征2012年地震的成分.但第二次ICA没有分解出代表P2点特征的独立成分,且P3点和P4点的特征出现在同一个独立成分中,造成这个现象的主要原因在于,2004年苏门答腊大地震的信号太强,其对于2012年Andaman地震的分析而言,是一个很强的噪声源,且两次地震的影响范围过于接近.为了更好地从GRACE重力信号中提取Andaman地震的时空特征,尽可能减小苏门答腊地震的影响,缩小时间序列的范围,利用ICA处理2006年1月至2014年12月的重力扰动变化时间序列,同时将参与旋转的模式数量由12个增加到24个.在得到的24个独立成分中,找到了与图 8中4个特征点相对应的IC,分别为IC4、IC8、IC17和IC22,如图 12图 13所示,图 12中的黑色圆点表示Andaman地震的震中,4个IC的空间模式都在震中附近出现了比较明显的重力扰动变化,相较于图 8中的4个特征点,P1、P2、P3和P4点的空间特征分别与IC22、IC4、IC17和IC8相对应,其时间序列也一一对应,其中图 13中的实线表示4个IC的时间序列,而虚线则表示其对应特征点归一化后的时间序列,两种时间序列还是比较吻合的,其中IC4和IC17反映了研究区域在2012年Andaman地震发生前后重力的增加,而IC8和IC22表现了重力的减少,这与图 8中的4个特征点一一对应,当然,两种时间序列在某些细节处有微小差异,之所以出现这种现象,是由于独立成分时间序列表现的是一定范围内所有点共同的特征,而图 8中特征点的时间序列仅仅代表这个点的特征,两者之间是共性与特性的关系,由此看来,由ICA得到的独立成分时间序列,更能代表研究区域的变化特征.对2006—2014年重力扰动数据的ICA分解结果进行分析,结合独立成分的空间特征与时间序列,可以认为,ICA方法可以较为完美地从GRACE重力扰动变化中提取得到与2012年Andaman地震有关的时空特征.

图 12 由2006—2014年数据得到的与2012年地震有关的独立成分空间模式 图中黑点表示震中.单位:μGal. Fig. 12 The spatial pattern of the independent components related to the 2012 earthquake based on data from 2006 to 2014 The lack dots in the picture indicate the epicenter. Unit: μGal.
图 13 由2006—2014年数据得到的与2012年地震有关的独立成分时间序列 图中的实线表示4个IC的时间序列,而虚线则表示与之对应的特征点归一化后的时间序列. Fig. 13 The time series of the independent components related to the 2012 earthquake based on data from 2006 to 2014 The solid line in the figure represents the time series of 4 ICs, and the dotted line represents the normalized time series of the corresponding feature points.

由上文可以得出结论,对于ICA分解,参与旋转的模式越多越好,模式越多表明信号来源越丰富,分解的效果也会越好.但模式越多,ICA算法对计算机的要求也越高,同时,所得到的模式也就越多,对于结果的分解也就越复杂.所以,在实际计算过程中,也不是说模式的数量越多越好,应该综合考虑,达到一个综合最佳的效果,一般而言,至少要包含90%的信息,这样进行的ICA分解才能得到比较合理的结果.

同时,需要确定的是,与常用的多项式拟合方法进行对比,PCA、ICA等基于统计信息的信号分解方法只是一种数学工具,单纯从信号统计特性的角度进行计算,有时可能会得到无法解释的分解结果,并不适合于定量计算,其优点在于定性分析,即不通过任何先验信息,通过信号分解得到大范围研究区域主要的变化特征,对研究区域有个大致的了解.至于精确的数值计算,还是需要利用多项式拟合等办法.在有关研究中,可以将两种方法结合起来,先利用ICA对原始数据进行定性分析,得到一定的先验信息,利用ICA分解得到的空间模式,找到需要重点研究的区域,再结合其对应的时间序列,得到相对准确的拟合模型,最后利用多项式拟合方法进行定量的分析,使分析的结果更加客观和准确.

5 讨论与结论

在一定程度上说,由GRACE球谐位系数计算得到的地面重力扰动是一种多维时空数据,可以利用ICA等信号分解方法对其进行分析,基于参与旋转的参数不同,ICA可以分为SICA和TICA,两者的区别在于分别是对PCA分解得到的空间模式和时间序列进行旋转,分别使其空间模式和时间序列尽可能的独立,本文所采用的是SICA,如果应用于其他的地震的分析中,需要结合分析的目的,如果想重点研究地震的空间特征,可以采用SICA,如果想分析地震的时间序列变化,可以采用TICA,当然,也可以将两种方法结合起来.

ICA方法可以有效地从GRACE观测结果中获得地震重力变化信号,以2004年12月的苏门答腊大地震为例,采用PCA和ICA两种信号分解方法,提取重力变化,并与变形较大的主震附近特征点的时间序列进行比较.结果表明,相较于PCA,ICA从原始信号中提取地震信号的能力更强,在提取信号比例较大的情况下,可以比较完美地获取与地震有关的模式,能分解出具有显著阶跃变化的独立成分,并且将震中附近不同变化特征的信号分解为不同的独立成分.增加参与旋转的PCA模式数量,ICA分解还提取到了与2012年Andaman地震有关的信号,而将与2004年苏门答腊大地震有关的月份删除掉,可以提升ICA关于Andaman地震的分解效果.由此可知,影响ICA分解效果的因素很多,如采用数据时间段的长短、参与旋转的模式数量以及计算重力扰动变化时的滤波方法等,都会对ICA的分解结果造成一定的影响.以采用数据时间段的长短为例,由于所采用时间段不同,4.2节中第一次ICA分解得到了2004年和2012年两次地震的信息,但第二次分解仅得到了2012年地震的信息,这是因为后者通过时间段的选择,消除了2004年地震的影响,可以这样认为,采用的数据时间段不同,其原始数据统计信息会发生一定的改变,但对时间序列数据主要的时空特征体现不会有影响,因为ICA的分解结果都只会表现原始数据中最主要的时空特征.同理,不同的滤波方法会对计算得到的地面重力扰动变化数值造成一定的影响,如果所选择的方法滤波效果不是很好,那么ICA分解结果也会受到比较明显的影响,故首先需要选取合适的滤波方法,得到较为准确的地面重力扰动变化数值,之后才进行独立成分分析.对于不同滤波方法获得的结果,独立主成分分析方法会导致一定的差异性,但差异的大小取决于滤波方法的差异.虽然影响ICA的最终结果的因素很多,有一点必须要明确,即ICA的分解结果只会表现原始数据中最主要的时空特征.

得到较为准确的地面重力扰动变化数值,除了利用合适的滤波方法及信号分解方法外,对各种误差和因素的分析也是必不可少的,对于2004年苏门答腊地震这样的“洋-陆”板块碰撞型地震,从定量分析的角度上看,“海水改正”是必不可少的,但本文的目的在于探究独立成分分析能否从GRACE数据中提取出代表地震的信号,属于定性分析,且“海水改正”对本文ICA分解的影响很难量化,故没有考虑“海水改正”,这还有待后续研究.

本文利用ICA方法可以比较完美地分离2004年和2012年两次地震的信号,由GRACE数据计算得到的地表重力扰动属于多维时空数据,利用ICA可以从时间与空间两方面进行分析,且得到的独立成分时空关系彼此一一对应,地震所造成的影响可以通过独立成分的时间与空间特征两方面表现出来,其空间特征会在震中附近出现比较明显的变化特征,而在其对应的时间序列中,也会在地震发生前后出现比较明显的重力突变,两者缺一不可.本文正是利用了两次地震差异较大的同震信号,从时空特征入手,实现了两次地震的分离.从某种角度上说,PCA和ICA等信号分解方法只是一种数学手段,没有利用先验信息,只基于原始数据统计信息得到的模式不一定能完全用自然事件进行解释,这也是这类信号分解方法的缺陷之一.相较于传统的多项式拟合方法,ICA方法可以得到信号的整体空间特征,有利于从宏观的角度对整个地球物理事件进行定性的分析,这也是常用的多项式拟合方法所不具备的特点.可以将ICA方法与传统的多项式拟合方法相结合,利用ICA进行定性分析,得到研究区域的时空变化特征,选取重点区域,利用多项式拟合进行定量分析,从而更好地识别像地震这样的地球物理事件,让分析结果更为客观、有效.

致谢  感谢GFZ所提供的GRACE卫星数据.
References
Cardoso J F, Souloumiac A. 1993. Blind beamforming for non-Gaussian signals. IEE Proceedings F(Radar and Signal Processing), 140(6): 362-370. DOI:10.1049/ip-f-2.1993.0054
Chen J L, Wilson C R, Tapley B D, et al. 2007. GRACE detects coseismic and postseismic deformation from the Sumatra-Andaman earthquake. Geophysical Research Letters, 34(13): L13302. DOI:10.1029/2007GL030356
Chen J L, Wilson C R, Blankenship D, et al. 2009. Accelerated antarctic ice loss from satellite gravity measurements. Nature Geoscience, 2(12): 859-862. DOI:10.1038/ngeo694
De Viron O, Panet I, Mikhailov V, et al. 2008. Retrieving earthquake signature in grace gravity solutions. Geophysical Journal International, 174(1): 14-20. DOI:10.1111/j.1365-246X.2008.03807.x
Forootan E, Kusche J. 2012. Separation of global time-variable gravity signals into maximally independent components. Journal of Geodesy, 86(7): 477-497. DOI:10.1007/s00190-011-0532-5
Han S C, Shum C K, Bevis M, et al. 2006. Crustal dilatation observed by GRACE after the 2004 Sumatra-Andaman earthquake. Science, 313(5787): 658-662. DOI:10.1126/science.1128661
Han S C, Sauber J, Luthcke S. 2010. Regional gravity decrease after the 2010 Maule (Chile) earthquake indicates large-scale mass redistribution. Geophysical Research Letters, 37(23): L23307. DOI:10.1029/2010GL045449
Han S C, Sauber J, Pollitz F. 2015. Coseismic compression/dilatation and viscoelastic uplift/subsidence following the 2012 Indian Ocean earthquakes quantified from satellite gravity observations. Geophysical Research Letters, 42(10): 3764-3772. DOI:10.1002/2015GL063819
Heki K, Matsuo K. 2010. Coseismic gravity changes of the 2010 earthquake in Central Chile from satellite gravimetry. Geophysical Research Letters, 37(24): L24306. DOI:10.1029/2010GL045335
Hyvärinen A. 1999. Survey on independent component analysis. Neural Computing Surveys, 2: 94-128.
Li J, Chen J L, Zhang Z Z. 2014. Seismologic applications of GRACE time-variable gravity measurements. Earthquake Science, 27(2): 229-245. DOI:10.1007/s11589-014-0072-1
Lu F. 2015. Global mass change analysis in recent ten years by GRACE RL05 data[Master's thesis] (in Chinese). Chengdu: Southwest Jiaotong University.
Luo Z C, Li Q, Zhang K, et al. 2012. Trend of mass change in the antarctic ice sheet recovered from the grace temporal gravity field. Science China Earth Sciences, 55(1): 76-82.
Matsuo K, Heki K. 2010. Time-variable ice loss in asian high mountains from satellite gravimetry. Earth and Planetary Science Letters, 290(1-2): 30-36. DOI:10.1016/j.epsl.2009.11.053
Matsuo K, Heki K. 2011. Coseismic gravity changes of the 2011 Tohoku-Oki earthquake from satellite gravimetry. Geophysical Research Letters, 38(7): L00G12. DOI:10.1029/2011GL049018
Mikhailov V, Tikhotsky S, Diament M, et al. 2004. Can tectonic processes be recovered from new gravity satellite data?. Earth and Planetary Science Letters, 228(3-4): 281-297. DOI:10.1016/j.epsl.2004.09.035
Ogawa R, Heki K. 2007. Slow postseismic recovery of geoid depression formed by the 2004 Sumatra-Andaman Earthquake by mantle water diffusion. Geophysical Research Letters, 34(6): L06313. DOI:10.1029/2007GL029340
Panet I, Mikhailov V, Diament M, et al. 2007. Coseismic and post-seismic signatures of the Sumatra 2004 December and 2005 March earthquakes in GRACE satellite gravity. Geophysical Journal International, 171(1): 177-190. DOI:10.1111/j.1365-246X.2007.03525.x
Preisendorfer R W. 1988. Principal Component Analysis in Meteorology and Oceanography. Amsterdam: Elsevier.
Richman M. B. 1986. Rotation of principal components. Journal of Climatology, 6(3): 293-335.
Rodell M, Velicogna I, Famiglietti J S. 2009. Satellite-based estimates of groundwater depletion in India. Nature, 460(7258): 999-1002. DOI:10.1038/nature08238
Sun W K, Okubo S. 2004. Coseismic deformations detectable by satellite gravity missions:A case study of Alaska (1964, 2002) and Hokkaido (2003) earthquakes in the spectral domain. Journal of Geophysical Research, 109(B4): B04405. DOI:10.1029/2003JB002554
Tapley B D, Bettadpur S, Ries J C, et al. 2004. GRACE measurements of mass variability in the earth system. Science, 305(5683): 503-505. DOI:10.1126/science.1099192
Wahr J, Molenaar M, Bryan F. 1998. Time variability of the Earth's gravity field:Hydrological and oceanic effects and their possible detection using GRACE. Journal of Geophysical Research, 103(B12): 30205-30229. DOI:10.1029/98JB02844
Wang L, Shum C K, Simons F J, et al. 2012. Coseismic slip of the 2010 MW8. 8 Great Maule, Chile, earthquake quantified by the inversion of GRACE observations. Earth and Planetary Science Letters, 335-336: 167-179. DOI:10.1016/j.epsl.2012.04.044
Wang W X, Shi Y L, Sun W K, et al. 2011. Viscous lithospheric structure beneath Sumatra inferred from post-seismic gravity changes detected by GRACE. Sci. China Earth Sci., 54(8): 1257-1267. DOI:10.1007/s11430-011-4217-y
Xue Y, Cheng J, Liu J, et al. 2014. Geodynamic mechanism of the 2012 Sumatra MW8.6 Earthquake and relationships with other great shocks in surrounding areas. Earth Science-Journal of China University of Geosciences (in Chinese), 39(4): 481-491. DOI:10.3799/dqkx.2014.046
Yan C D, Xu Y. 2018. Progress of application of GRACE satellite gravity in research of earthquakes. Progress in Geophysics (in Chinese), 33(3): 1005-1012. DOI:10.6038/pg2018BB0429
Zhang G Q, Fu G Y, Zhou X. 2015. The evolution process of the gravitational field after the Sumatra MW9.3 earthquake from GRACE RL05 data. Journal of Geodesy and Geodynamics (in Chinese), 35(2): 303-308.
Zhang K L, Gan W J, Zhou X. 2014. Detection of coseismic changes of great earthquakes in grace time-variable gravity field with empirical orthogonal functions:A case study of the MW9.0 Tohoku-Oki earthquake. Seismology and Geology (in Chinese), 36(3): 763-774.
Zhang Z Z, Chao B F, Lu Y, et al. 2009. An effective filtering for GRACE time-variable gravity:Fan filter. Geophysical Research Letter, 36(17): L17311. DOI:10.1029/2009GL039459
Zheng W, Xu H Z, Zhong M, et al. 2009. Effective processing of measured data from GRACE key payloads and accurate determination of Earth's gravitational field. Chinese J. Geophys. (in Chinese), 52(8): 1966-1975. DOI:10.3969/j.issn.0001-5733.2009.08.003
Zheng W, Xu H Z, Li Z W, et al. 2017. Precise establishment of the next-generation Earth gravity field model from HIP-3S based on combination of inline and pendulum satellite formations. Chinese J. Geophys. (in Chinese), 60(8): 3051-3061. DOI:10.6038/cjg20170813
Zhou J C, Sun H P, Xu J Q. 2009. Validating global hydrological models by ground and space gravimetry. Chinese Sci. Bull., 54(9): 1534-1542.
Zhou J C, Sun H P, Xu J Q. 2013. Coseismic gravity signals detection from GRACE results by EOF method. Journal of Geodesy and Geodynamics (in Chinese), 33(3): 25-29.
Zhou X, Sun W K, Fu G Y. 2011. Gravity satellite GRACE detects coseismic gravity changes caused by 2010 Chile MW8.8 earthquake. Chinese J. Geophys. (in Chinese), 54(7): 1745-1749. DOI:10.3969/j.issn.0001-5733.2011.07.007
Zhou X, Sun W K, Zhao B, et al. 2012. Geodetic observations detecting coseismic displacements and gravity changes caused by the MW=9.0 Tohoku-Oki earthquake. Journal of Geophysical Research, 117(B5): B05408. DOI:10.1029/2011JB008849
Zou Z B, Li H, Wu Y L, et al. 2016. Spatial and temporal characteristics of long-term satellite gravity change in the epicenter of MW9.0 Japan earthquake and its surrounding regions. Acta Seismologica Sinica (in Chinese), 38(3): 417-428.
卢飞. 2015.基于GRACE RL05数据的近十年全球质量变化分析[硕士论文].成都: 西南交通大学. http://d.wanfangdata.com.cn/Thesis_Y2813915.aspx
王武星, 石耀霖, 孙文科, 等. 2011. 利用GRACE观测数据研究苏门答腊区域的黏滞性结构. 中国科学:地球科学, 41(6): 773-783.
薛艳, 程佳, 刘杰, 等. 2014. 2012年苏门答腊MW8.6地震成因及与周边大震关系. 地球科学(中国地质大学学报), 39(4): 481-491.
严畅达, 徐亚. 2018. GRACE卫星重力在地震研究中的应用进展. 地球物理学进展, 33(3): 1005-1012. DOI:10.6038/pg2018BB0429
张国庆, 付广裕, 周新. 2015. 利用GRACE卫星数据提取苏门答腊地震同震和震后重力变化. 大地测量与地球动力学, 35(2): 303-308.
张克亮, 甘卫军, 周新. 2014. GRACE卫星重力场同震变化的经验正交函数分解:以日本MW9.0地震为例. 地震地质, 36(3): 763-774. DOI:10.3969/j.issn.0253-4967.2014.03.017
郑伟, 许厚泽, 钟敏, 等. 2009. GRACE卫星关键载荷实测数据的有效处理和地球重力场的精确解算. 地球物理学报, 52(8): 1966-1975. DOI:10.3969/j.issn.0001-5733.2009.08.003
郑伟, 许厚泽, 李钊伟, 等. 2017. 联合串行式和钟摆式卫星编队精确建立下一代HIP-3S地球重力场模型. 地球物理学报, 60(8): 3051-3061. DOI:10.6038/cjg20170813
周江存, 孙和平, 徐建桥. 2009. 用地表和空间重力测量验证全球水储量变化模型. 科学通报, 54(9): 1282-1289.
周江存, 孙和平, 徐建桥. 2013. EOF方法检测GRACE卫星重力结果中的同震重力变化. 大地测量与地球动力学, 33(3): 25-29.
周新, 孙文科, 付广裕. 2011. 重力卫星GRACE检测出2010年智利MW8.8地震的同震重力变化. 地球物理学报, 54(7): 1745-1749. DOI:10.3969/j.issn.0001-5733.2011.07.007
邹正波, 李辉, 吴云龙, 等. 2016. 日本MW9.0地震震区及其周缘2002-2015年卫星重力变化时空特征. 地震学报, 38(3): 417-428.