1. 引言
地震同震滑动分布反演问题是不适定的(ill-posed),其反演结果具有非唯一性。为了增强反演过程的稳定性、减少反演结果的非唯一性,需要加入先验约束条件进行正则化反演。根据不同的先验约束条件,模型约束目标函数主要有最小模型、最平坦模型和最光滑模型,分别对应模型参数的平方和最小、模型参数导数的平方和最小和模型参数二阶导数的平方和最小。Masterlark等通过滑动量的二阶导数来约束断层间的滑动,实现了走向滑动量或倾向滑动量的光滑约束 [1] ;Tong等采用滑动量的一阶导数最小进行约束,实则是一种最平坦模型约束形式 [2] ;Maerten等构造了三角形单元剖分的尺度依赖伞算子(scale-dependent umbrella operator),通过结合点源位错理论,实现了复杂断层面构造的同震滑动分布反演。另外,即使在反演中加入了先验光滑约束条件,但由于GPS测量数据的分辨率仍然无法保证断层面的滑动分布模型的合理性 [3] ;此时,除了光滑约束外,需要借助于更强的约束;Barnhart和Lohman尝试了固定滑动角,再利用最小二乘法进行求解 [4] ;Pollitz等将滑动角固定到纯走滑或纯倾滑的情况,只反演滑动量 [5] ;Fukuda和Johnson将断层面滑动分布反演问题中线性参数和非线性参数分离,混合使用线性与非线性反演方法,但反演计算效率较低 [6] ;Ziv等采用非负最小二乘求解了同震反演过程中的超定问题,对滑动量进行了非负约束 [7] 。因此,先验光滑约束是一种弱约束,而滑动角约束则是一种强制约束。
在正则化反演计算过程中,正则化因子选取的合理与否直接影响着反演结果的分辨率和反演过程的稳定性,决定着反演求解过程成功和失败的命运。正则化因子决定了反演的主要拟合对象,过大则主要拟合先验约束信息,过小则主要拟合实际观测数据,因此,正则化因子的不同取值对反演结果有重要影响。Yabuki和Matsu’ura [8] 、Fukuda和Johnson [9] 通过贝叶斯信息准则评判正则化因子,同时计算出滑动量分布以及数据相对权比,但计算时间较长;Hreinsdóttir等利用交叉验证准则确定正则化参数,然而这种估值准则的实施过程较为复杂,并牵涉到矩阵的分解操作,而且当观测数据存在相关噪声时,这一方法求得的正则化参数往往并不理想 [10] ;Barnhart和Lohman利用模型重构矩阵形成的扩散函数选取正则化因子,该策略不需要曲线出现拐点,且对噪声数据反演非常灵敏 [4] ;Jiang等利用加权拟合差与粗糙度的关系确定正则化因子,两者之间曲率变化最大的点即为最佳正则化因子 [11] ;许才军等提出了基于方差分量估计方法确定正则化因子,该方法不仅可以确定不同数据集的相对权重,而且还能得到正则化因子的值,突破了方差分量估计方法仅用于数据定权的应用范围 [12] 。
本文将利用光滑约束的正则化反演方法,计算断层面上的滑动量分布,而其中的Green函数采用Okada提出的矩形有限断层的同震位移计算方法 [13] [14] 。反演计算过程中,构建断层面非均匀离散的滑动量分布光滑约束模型,设计一种快速稳定的正则化因子求取方法,并采用非负最小二乘进行求解。文中首先通过理论滑动分布模型的最小二乘与非负最小二乘反演对比分析,验证反演方法的有效性;然后利用噪声数据反演的数值试验,验证反演算法的稳定性;最后通过对2005年Nias地震的同震滑动分布正则化反演,探讨数据分布和断层构造的差异、以及不同反演计算方法对重构同震滑动分布的影响。
2. 模型与算法
2.1. 正则化方法
考虑到断层面同震滑动分布反演问题是不适定的(ill-posed),其反演结果具有非唯一性,也就是不同滑动量分布模型所对应的响应数据与观测数据具有同样的精度拟合。为了改善解的稳定性和非唯一性问题,通常是引入Tikhonov的正则化思想 [15] :
(1)
式中,
为总目标函数;
为模型参数向量,即断层面同震滑动量,包括走向滑动量和倾向滑动量;
为正则化因子;
为观测数据与预测数据之差的平方和(即数据目标函数);
为稳定器(即模型约束目标函数)。
因此,断层面同震滑动分布反演问题的总目标函数可表示为
(2)
式中:
为模型参数空间到地表形变空间的格林函数;
为观测到的地表形变数据向量;
为数据权重系数矩阵,假设观测数据的噪声满足高斯正态分布、观测点个数为
,则有
;
为模型权系数矩阵,通常以微分、积分或矩阵形式给出。
将上述目标函数对
求导并令其等于0,可得线性方程组:
(3)
则可得断层面滑动分布反演问题的最小二乘解:
(4)
为了避免最小二乘反演可能得到不合理的滑动分布,可以对滑动量进行非负约束,即
。这里,我们采用非负最小二乘法进行反演求解,反演所得的走向滑动只能是左滑或右滑 [16] [17] 。
2.2. 光滑约束模型构建
反演过程中,将断层面离散化为均匀或非均匀的子断层,通常需要对子断层上的滑动量施加一定的光滑约束条件,从而避免相邻子断层之间的走向滑动量或倾向滑动量的大小存在着显著差异。本文采用二阶有限差分算子,其目的是相邻断层片之间滑动量的二阶导数值最小 [18] 。对于非均匀离散化断层的走向滑动量s或倾向滑动量s,如图1所示,其光滑约束的拉普拉斯二阶差分算子可以表示为:
(5)
式中:
为第i行、第j列子断层面的走向滑动量或倾向滑动量;
和
为走向方向上相邻子断层中心点的距离;
和
为倾向方向上子断层中心点的距离。
因此,光滑矩阵
可以构建为:
(6)
这里的nc为倾向方向上断层离散化的单元数。于是,同震滑动面分布反演计算过程中的光滑约束矩
![](//html.hanspub.org/file/6-1770499x35_hanspub.png)
Figure 1. Schematic figure of fault patches for the smoothed coseismic slip
图1. 同震滑动量光滑约束的断层片示意图
阵可写为
。
2.3. 正则化因子求取
采用Tikhonov正则化方法解决反演问题的不适定性时,式(1)中的正则化因子
将在反演过程中起着关键性的作用 [19] [20] [21] 。当
时,由于反演方程的奇异性,必然会增大解的方差,导致反演过程的不稳定,最终会导致反演失败。当
时,解的光滑性和反演过程的稳定性最好,但是以严重降低解的分辨率为代价。因此,正则化因子的选择是否合理掌握着反演求解过程的成功与失败。本文采用一种改进的L曲线法,其目的是快速稳定的寻找到最佳正则化因子。
将走向滑动量和倾向滑动量的光滑约束矩阵
转换为主对角线元素值为−1的形式,并设
,取正则化因子值的最小值与最大值分别为:
(7)
和
(8)
式中:M为模型参数个数。
确定了正则化因子的取值区间后,令
和
,则反演过程中的最佳正则化因子计算公式可以表述为:
(9)
式中:
和
分别表示
和
的一阶导数;
和
分别表示
和
的二阶导数。正则化因子的搜索过程如图2所示,最佳正则化因子是曲线曲率最大的地方。
3. 理论模型反演试算
半无限均匀弹性空间中存着在长度为15 km,宽度为8 km,走向角为37˚,倾角为60˚的断层,其上顶面距离地表2 km。将断层划分为
的子断层,并在各断层面上加入一定的滑动量(见图3(a)、图3(b)和图3(c)),滑动角设定为63.4°,其中最大走向滑动量和最大倾向滑动量分别为10 m和20 m,总的最大滑动量为22.4 m。
根据上述理论断层分布和半无限均匀弹性空间的位错模型,在地表观测位移中加入5%的高斯分布噪声,可以计算出均匀分布的441个站点(见图4)的三维同震位移。图3中,白色箭头代表地表水平位移,其最大值约为1.79 m;等值线底图代表垂直位移,其平均下沉的最大值约为4.02 m。
反演计算过程中,实际计算的模型参数个数为
。图3(d)~图3(f)为最小二乘反演计算结果,分别对应的是走向滑动分布、倾向滑动分布和总体滑动分布,其分布特征与理论模型的滑动分布吻合得很好,但走向滑动量和倾向滑动量都出现了负值,这与理论模型的实际滑动情况不符。图3(g)~图3(i)为非负最小二乘反演计算结果,分别对应的是走向滑动分布、倾向滑动分布和总体滑动分布,其最大走向滑动量值为10.2 m,最大倾向滑动量值为21.1 m,总的最大滑动量值为23.1 m,且分布特征同样与理论模型的滑动分布吻合得很好,并改善了出现滑动量负值的计算结果。
![](//html.hanspub.org/file/6-1770499x57_hanspub.png)
Figure 2. The best regularization parameter used in the inversion
图2. 反演过程中的最佳正则化因子
![](//html.hanspub.org/file/6-1770499x58_hanspub.png)
Figure 3. Coseismic slip distribution of a synthetic model and inversion results
图3. 理论模型的同震滑动分布及反演结果图
将模型反演计算的预测数据与实际的观测数据进行对比分析,如图4所示,最小二乘反演与非负最小二乘反演的预测数据都与实际观测数据吻合,这也说明了正则化反演方法的准确性。利用反演得到的断层面上的滑动分布,取
,可以计算地震发生时释放的地震矩为
,比理论模型实际的地震矩
稍大。
表1给出了光滑约束正则化反演方法对噪声数据反演的数值试验结果,其中:第一列数据给出了理论模型中观测数据加入噪声的百分比,第二列数据给出了反演计算的噪声水平,第三列数据显示了实际噪声的数量级(其中
为添加到观测数据中的高斯噪声)。可以看出,估算噪声与实际噪声的数量级吻合的很好。观测数据的误差越大,反演结果的可信度就越低。但值得注意的是,由于观测误差的随机性,反演参数的误差并不一定大于观测误差,验证了反演算法的有效性和稳定性。
4. 2005年Nias地震
Nias地震是2005年3月28日发生在印度尼西亚的尼亚斯岛附近海域的8.6级强烈地震,其震中位置为(2.074˚N, 97.013˚E),震源深度为30 km。该地震同震形变观测数据选用Kreemer等提供的解算结果 [22] ,表2给出了10个GPS观测站点的三维同震位移。通过数据可以看出,地表水平移动的最大值约4.5 m,平均下沉的最大值约2.9 m。
反演计算过程中,断层面参数采用Konca等给出的单段断层模型:走向角为325˚,倾向角为10˚,断层面长度和宽度分别为416 km和320 km,且断层面顶部距离地表3.21 km [23] 。将整个断层面划分为
的子断层,则反演计算的模型参数个数为
。利用光滑约束正则化反演所得的滑动量分布情况如图5所示,其中图5(a)揭示的走向滑动量最大值为3.6 m、图5(b)揭示的倾向滑动量最大值为12.6 m、图5(c)揭示的总体滑动量最大值为12.8 m,这与Konca等给出的反演结果基本相似。另外,
![](//html.hanspub.org/file/6-1770499x65_hanspub.png)
Figure 4. Coseismic displacement of a synthetic model and inversion results; (a)Observed data; (b)LS inversion; (c)NNLS inversion
图4. 理论模型的三维同震位移及反演结果图;(a) 观测数据;(b) 最小二乘反演;(c) 非负最小二乘反演
![](Images/Table_Tmp.jpg)
Table 1. Numerical experiments with regularized inversion for noise data
表1. 噪声数据反演的数值试验
![](Images/Table_Tmp.jpg)
Table 2. Coseismic displacement of GPS sites for Nias earthquake
表2. Nias地震GPS观测站点的同震位移
![](//html.hanspub.org/file/6-1770499x68_hanspub.png)
![](//html.hanspub.org/file/6-1770499x69_hanspub.png)
Figure 5. Coseismic slip distribution of the 2005 Nias earthquake; (a) Strike-slip distribution; (b) Dip-slip distribution; (c) Total-slip distribution
图5. 2005年Nias地震的同震滑动分布特征;(a) 走滑分布;(b) 倾滑分布;(c) 总体滑动分布
走向滑动量和倾向滑动量出现极大值的区域不一致,且两者的分布特征不同,这说明反演计算的各个子断层的滑动角是变化的。
通过反演得到的滑动量分布来计算观测点上的GPS站点形变值,图5给出了观测点上形变预测值和观测值的比较,其中五角星是USGS公布的震中位置,箭头代表观测点的水平形变值。可以看出反演得到的断层面上的滑动分布在地表观测点上引起的形变预测值与观测值拟合得很好。
利用反演得到的断层面的滑动分布,可以计算地震发生时释放的地震矩为
,与Konca等得到的地震矩
接近。反演得到的Nias地震震级为
,与USGS公布的结果一致。
5. 结论与讨论
本文实现了同震滑动分布的正则化反演算法,构建了断层面非均匀离散的走向滑动量与倾向滑动量的光滑约束模型,并在正则化因子计算过程中提出了一种改进的L曲线计算方法。通过理论滑动分布模型的最小二乘与非负最小二乘反演对比分析,以及实际震例的非负最小二乘反演求解,均证实了正则化反演方法是合理且有效的。通过噪声数据反演的数值试验,验证了正则化反演算法的有效性和稳定性。
在2005年Nias地震的同震滑动分布反演计算过程中,断层面参数采用Konca等给出的单段断层模型,而反演所得的滑动量分布情况类似于Konca等获得的反演结果,且滑动模型对应地震矩和地震震级同样吻合。这说明在同震形变反演过程中,断层构造一致时,不同反演计算方法所得的同震滑动分布情况大同小异。
对于实际地震的反演计算,特别是俯冲型地震,滑动角一般为
(一般情况下
),而本文的非负最小二乘只是把滑动角设定为
或
,这与真实情况存在一定的偏差。因此,下一步研究重点将对线性反演问题进行非线性处理,把反演计算的模型参数转换为滑动量和滑动角,进行非线性迭代反演。
基金项目
国家自然科学基金项目(41674080)和国家级大学生创新创业训练计划项目(201710533236)联合资助。