地震同震滑动分布反演平滑因子的确定

王乐洋，赵 雄

1. 东华理工大学测绘工程学院，江西 南昌 330013； 2. 流域生态与地理环境监测国家测绘地理信息局重点实验室，江西 南昌 330013； 3. 江西省数字国土重点实验室，江西 南昌 330013

在地震机制研究领域中，利用大地测量资料研究地震机制是地学领域的前沿研究热点[1]。大地测量数据(如GPS、InSAR数据等)因观测范围广、精度高，反演出的滑动分布结果更为精细、全面，近年来常被国内外学者用于进行地震同震滑动分布反演研究[2-8]。在震源机制研究中，位错模型是反演震源参数的有效模型[9-11]，一般采用弹性半空间矩形位错模型进行地震同震滑动分布反演[12-13]。

在地震同震滑动分布反演过程中，地表形变位移与滑动分布参数之间为线性关系。滑动分布参数与地表形变位移的线性模型依赖于格林函数矩阵的构造，而格林函数矩阵元素通常是确定断层的几何参数后，将断层破裂面进一步离散成若干断层子块，由每个子块的单位走滑与倾滑引起的地表形变位移构成。另外，格林函数矩阵在地震滑动分布反演过程中通常是病态的，一般采用Tikhonov正则化方法[14]解决。

在正则化过程中，平滑因子的确定是关键问题之一，不同的平滑因子会使反演结果大不相同。目前主要方法有广义交叉核实法(generalized cross validation,GCV)[15]、岭迹法[16-17]、L曲线法[18-20]。GCV法理论严密，但有时GCV方法函数绘制的曲线过于平缓，很难定位到最优平滑因子[21]；岭迹法计算简单，但是在确定平滑因子的时候具有一定的主观性；在地震同震滑动分布反演中，目前最常用的平滑因子确定方法为L曲线法，但是L曲线法在确定平滑因子过程中具有以下几个缺点：①无法从L曲线图中根据横纵坐标轴直观读取平滑因子的大小；②L曲线法过于依赖数据拟合精度[22]，求解过程可能不收敛[23]；③计算时间过长。

针对以上问题，本文提出一种确定地震同震滑动分布反演平滑因子的新方法——折中相交曲线法(eclectic intersection curve method，EI)(简称EI曲线法)。通过模拟试验以及拉奎拉与台湾美浓实际震例反演验证EI曲线法确定平滑因子的优势与可行性。

1 EI曲线法确定平滑因子的基本原理

1.1 同震滑动分布基本方程

在同震滑动分布反演中，地表同震位移与断层滑动量之间为线性关系，如式(1)所示

d=Gm+ε

(1)

式中，d表示地表形变数据；G表示格林函数矩阵；m为滑动参数；ε为随机观测误差。

在同震滑动分布反演中，通常将断层面划分成n个矩形单元，基于矩形位错模型计算矩形单元的单位走滑量和单位倾滑量引起地表k个观测点的形变值，构成了k×2n个格林函数矩阵G。将断层面离散化处理，直接用式(1)求解会出现病态问题，从而使滑动分布解不稳定，因为此时格林矩阵G是严重秩亏的。为了避免这种情况，通常对断层面施加一定的平滑约束条件，一方面增加虚拟观测数据解决系数矩阵G的秩亏问题，另一方面约束了相邻断层单元间出现大的梯度变化。

本文采用目前较成功用于断层面平滑约束的拉普拉斯二阶平滑矩阵对断层单元进行约束[24]。根据Tikhonov正则化原理，相应的准则为

(2)

由式(1)、式(2)可得地震滑动分布解为

m=(GTP1G+αTTP2T)-1GTP1d

(3)

1.2 L曲线法的基本原理

图1 L曲线图Fig.1 L curve

1.3 EI曲线法基本原理

针对以上问题，EI曲线法首先利用赫尔默特方差分量估计法[1,25]统一反演过程中各类数据单位权方差，在统一单位权方差过程中需要加入权修正因子[26]a、b进行修正，使σ01=σ02，此时相应的目标函数可以表达如下

(4)

式中，a、b为实际观测数据与虚拟观测数据内部权阵P1、P2的权修正因子。

图2 EI曲线图Fig.2 EI curve

综上，利用EI曲线法确定地震滑动分布平滑因子步骤可总结如下：

(1) 利用方差分量估计法统一实际观测数据与虚拟观测数据单位权方差：即使σ01=σ02，并求得a、b。

2 模拟试验

为了验证本文方法的可行性，本文作了以下模拟试验。模拟试验中，模拟断层面几何参数如下：断层面几何中心位置X=0 km，Y=0 km，断层顶深0 km，长度30 km，宽度30 km，走向角60°，倾角45°，滑动角45°。模拟GPS 3方向观测数据如图3所示，根据点位与震中距离设置形变点密度，并给形变点施加观测误差N(0，32mm2)。

在模拟试验中，本文分别利用L曲线法、EI曲线法确定平滑因子，模拟试验中平滑因子α均以0.01为初值、0.001为步长、1为终值绘制L曲线图、EI曲线图(其中，L曲线法确定平滑因子大小为0.179，EI曲线法确定的平滑因子大小为0.095。在利用EI曲线法确定平滑因子过程中，利用方差分量估计法确定实际观测数据与虚拟观测数据内部权阵修正因子大小分别为a=1，b=0.003 7)，并将两种方法确定的平滑因子进行地震滑动分布反演。利用L曲线法、EI曲线法确定平滑因子如图4(a)、(d)所示；两种方法反演地震滑动分布结果以及两种方法反演地震滑动分布与模拟值的差值分布如图4(b)、(e)、(c)、(f)所示。

模拟试验预设参数以及利用L曲线法、EI曲线法反演滑动分布各参数结果见表1。两种方法确定平滑因子反演地震滑动分布部分参数结果与模拟值的差异，以及EI曲线法反演地震滑动分布部分参数结果，较L曲线法精度改进程度见表2。

表1 模拟试验反演地震滑动分布结果

表2模拟试验两种方法反演各参数结果与模拟值的差异

Tab.2Comparisontheresultsofeachparametersandsimulationvaluesoftwomethodsinsimulationexperiments

与模拟值的差值最大滑动量/m平均滑动量/m平均滑动角/(°)L曲线法0.16720.00640.9012EI曲线法0.09990.00340.8826L-EI∗5.23%0.95%0.04%

注：L-EI*表示利用EI曲线法确定平滑因子反演地震滑动分布各参数结果较L曲线法改进程度。

表1给出了利用L曲线法和EI曲线法确定平滑因子反演地震同震滑动分布各参数的结果及两种方法确定平滑因子所用时间；表2给出了两种方法反演地震同震滑动分布各参数结果与模拟值的差异。在各参数反演精度上，从表1、表2可以看出，利用EI曲线法确定的平滑因子反演地震滑动分布结果在最大滑动量、平均滑动量、矩震级等参数反演结果均要优于L曲线法反演结果；在计算时间上，在相同的初值、步长、终值情况下，EI曲线法比L曲线法在计算时间上减少了107.784 3 s，大大地提高了计算效率；在均方根误差上，EI曲线法较L曲线法减小0.1 mm，进一步说明EI曲线法确定平滑因子反演地震滑动分布结果精度要高于L曲线法。

3 实际震例反演

3.1 拉奎拉地震滑动分布反演

2009年4月6日，意大利中部拉奎拉地区发生Mw6.3级地震。据美国地质调查局(USGS)网站发布，地震震中位置为(13.334°E，42.334°N)，震源深度为8.8 km。地震造成1000余人的伤亡，大量的房屋被损毁或破坏[27]。本文从文献[6，28]中获取拉奎拉地震InSAR降轨数据位移1254个，升轨位移数据1282个，并利用文献[29]中通过非线性反演方法反演得到的拉奎拉地震断层的走向角、倾角及断层中心位置等参数来固定断层破裂面。在此基础上将断层破裂面沿断层走向、倾向把断层长度、宽度扩展至30 km、27 km，并将破裂面延伸至地表，将断层面均匀剖分成1.5×1.5 km大小的矩形单元，共得到了360个矩形断层单元。分别采用EI曲线法与L曲线法确定平滑因子进行地震滑动分布反演。

在拉奎拉实际震例中，平滑因子α均以0.05为初值、0.001为步长、1为终值绘制L曲线图、EI曲线图(其中，L曲线法确定平滑因子大小为0.588与文献[28]中用L曲线法确定的平滑因子大小0.6，两者大小相近；EI曲线法确定的平滑因子大小为0.159。在利用EI曲线法确定平滑因子过程中，利用方差分量估计法确定实际观测数据与虚拟观测数据内部权阵修正因子，大小分别为a=1，b=0.002 4)。利用EI曲线法和L曲线法确定平滑因子如图5(a)、(c)所示；两种方法反演拉奎拉地震成果及两种方法反演滑动分布结果差值分布如图5(b)、(d)、(e)所示。

利用L曲线法、EI曲线法确定的平滑因子反演拉奎拉地震滑动分布各参数结果见表3；近年来国内外部分其他学者反演拉奎拉地震滑动分布各参数结果见表4。

表3 用L曲线法、EI曲线法确定平滑因子反演拉奎拉地震滑动分布结果

表4 拉奎拉地震断层滑动参数

由图5(b)、(d)可知，利用L曲线法、EI曲线法确定的平滑因子进行拉奎拉地震反演时，其反演结果均显示拉奎拉同震滑动分布区域主要分布在地下4～15 km，该结果同文献[28]研究成果一致；对比表3、表4可以看出，本文两种方法确定的平滑因子反演拉奎拉地震滑动分布结果均在其他学者研究范围内。对于最大滑动量，利用L曲线法和EI曲线法反演最大滑动量大小分别为1.01、1.02 m，文献[6]联合了不同波长、不同入射角的Enviast和ALOS卫星差分干涉数据，并联立了震区部分GPS三维形变数据，通过弹性三角位错模型，反演得到滑动分布的最大滑动量为1.07 m，对于其最大滑动量略大于本文反演结果，这可能与约束地表的数据源有关。文献[28]利用最小二乘和总体最小二乘反演最大滑动量均为0.95 m，略小于本文两种方法反演结果；对于平均滑动角，本文两种方法反演的结果分别为-101.3°、-103.9°，该研究结果均在其他学者研究范围内，并且与文献[6]反演结果相近；对于矩震级参数，本文两种方法反演结果相近，且在其他学者研究范围内。另外，文献[27]联合ASAR数据反演得到的断层滑动量为0.66 m，给出的矩张量为2.80×1018N·m；文献[31]利用拉奎拉震中附近的GPS观测站数据，研究了拉奎拉地震的同震滑动分布和震后的余震滑动分布，得到的最大滑动量为1.1 m，矩张量为3.90×1018N·m。

3.2 台湾美浓地震滑动分布反演

2016年2月6日3时57分(北京时间)我国台湾省高雄市美浓区发生了Mw6.4(global centroid moment tensor,GCMT)地震(下称美浓地震)。此次地震给台湾南部地区带来巨大的人员伤亡及财产损失。本文采用文献[32]获取美浓GPS站位移数据225个参与计算美浓地震滑动分布反演。利用文献[32]中均匀滑动分布反演的断层走向、倾向、断层深度等参数，作为本文通过多峰值颗粒群算法[29]进行非线性断层参数反演的初始参数大小。通过非线性反演最终确定适合本文数据的美浓地震断层参数走向角为strike=303.6°，倾角为dip=16.2°，震源中心深度depth=12.3 km。将断层长和宽扩展至44、54 km，并将其按照2 km×2 km矩形划分为594个单元。为了避免少量滑动单元出现方向矛盾问题，采用非负约束最小二乘法进行线性滑动分布求解[33]。分别采用EI曲线法与L曲线法确定平滑因子进行地震滑动分布反演。

在美浓实际震例反演中，平滑因子α均以0.05为初值、0.002为步长、1为终值绘制L曲线图、EI曲线图(其中，L曲线法确定的平滑因子大小为0.445，EI曲线法确定的平滑因子大小为0.328。在利用EI曲线法确定平滑因子过程中，利用方差分量估计确定实际观测数据与虚拟观测数据内部权阵修正因子，大小分别为a=1，b=0.831 1)。利用EI曲线法和L曲线法确定平滑因子如图6(a)、(c)所示；两种方法反演美浓地震结果及两种方法反演结果差值分布如图6(b)、(d)、(e)所示。

图3 GPS 3方向模拟观测点Fig.3 Distribution of three directions simulated points of GPS

图4 模拟试验两种方法反演地震滑动分布结果及解的残差Fig.4 Slip distribution results of simulation experiment and the residuals of two methods

图6 两种方法反演美浓地震滑动分布结果及差值分布Fig.6 The slip distribution inversion results of Meinong earthquake and the residuals of two methods

利用L曲线法、EI曲线法确定的平滑因子反演美浓地震滑动分布各参数结果见表5；近年来国内外部分其他学者反演美浓地震滑动分布各参数结果见表6。

表5 用L曲线法、EI曲线法确定平滑因子反演美浓地震滑动分布结果

表6 美浓地震断层滑动参数

从图6(b)、(d)可以看出，利用L曲线法、EI曲线法确定的平滑因子进行美浓地震反演时，其反演结果均显示美浓地震滑动分布区域主要分布在地下9～14 km，该结果同文献[32]反演结果一致。表6列举了美浓地震的部分研究成果，表5、表6显示对于最大滑动量，本文两种方法反演结果分别为0.453、0.479 m。文献[32]联合InSAR和GPS数据反演美浓地震结果显示最大倾滑量为0.517 m、最大走滑量为0.553 m，略大于本文两种方法反演结果；对于平均滑动角，本文两种方法反演的结果分别为45.32°、43.45°，其中，利用L曲线法反演结果与文献[34]反演结果相近。而EI曲线法反演结果在文献[34]与文献[32]反演结果之间；对于矩震级参数，本文两种方法反演结果相近，且在其他学者研究范围内；文献[35]反演的矩震级结果为Mw6.52，对于其反演结果较其他学者反演结果偏大的原因可能与没有较好地拟合GPS数据有关[32]。

对于本文所做的模拟试验以及拉奎拉与台湾美浓实际震例试验，作出以下几点分析：

(1) 对于两种方法确定平滑因子反演地震滑动分布各参数精度，由模拟试验可知，利用EI曲线法确定平滑因子反演滑动分布的各参数结果与模拟值的差值均小于L曲线法。模拟试验及两个实际震例反演的均方根误差结果表明EI曲线法反演的均方根误差均小于L曲线法，从而进一步说明EI曲线法确定的平滑因子反演地震滑动分布精度要高于L曲线法。

(2) 对于计算效率，在模拟试验、拉奎拉及美浓实际震例中EI曲线法确定平滑因子时间较L曲线法分别减少107.784 3、125.664 1、74.729 4 s。总体上，EI曲线法较L曲线法效率提高2～7倍，故利用EI曲线法确定平滑因子较L曲线法可以有效地提高计算效率。

4 结 论

本文提出的EI曲线法为地震滑动分布反演中平滑因子的确定增加了一条新的途径。本文通过模拟试验及拉奎拉与美浓实际地震反演表明，EI曲线法较L曲线法确定平滑因子具有定位准确、反演精度高、计算效率高，以及EI曲线图横坐标直观反映平滑因子大小、无需依赖于数据拟合度等优势。本文提出的EI曲线法可以应用在地震滑动分布反演中，在其他领域的应用有待于进一步的研究与探索。

致谢：感谢武汉大学温扬茂博士提供意大利拉奎拉地震的InSAR数据。