用U曲线法确定地震同震滑动分布反演正则化参数

王乐洋;赵雄

【摘 要】针对地震滑动分布反演中正则化参数选取问题,提出利用U曲线法确定地震滑动分布反演正则化参数.利用U曲线法、L曲线法设计模拟实验,并将两种方法应用到芦山实际震例反演中.模拟实验以及芦山实际震例反演结果表明,利用U曲线法确定正则化参数反演地震滑动分布结果与L曲线法相比具有反演精度高、无需依赖数据拟合精度等优势.

【期刊名称】《大地测量与地球动力学》 【年(卷),期】2018(038)011 【总页数】6页(P1196-1201)

【关键词】同震滑动分布;正则化参数;U曲线法;L曲线法;芦山地震 【作 者】王乐洋;赵雄

【作者单位】东华理工大学测绘工程学院,南昌市广兰大道418号,330013;流域生态与地理环境监测国家测绘地理信息局重点实验室,南昌广兰大道418号,330013;江西省数字国土重点实验室,南昌广兰大道418号,330013;东华理工大学测绘工程学院,南昌市广兰大道418号,330013;流域生态与地理环境监测国家测绘地理信息局重点实验室,南昌广兰大道418号,330013 【正文语种】中 文 【中图分类】P315

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

在Tikhonov正则化过程中，正则化参数的选取尤为重要，不同的正则化参数会使反演结果大不相同。目前常用的方法有广义交叉核实法(generalized cross validation, GCV)[3]、岭迹法[4]、L曲线法[5-6]。GCV法理论上能获得较优的正则化参数，但有时GCV方法函数绘制的曲线过于平缓，很难定位到最优正则化参数[7]；岭迹法计算简单，但在确定正则化参数时存在主观性。在地震同震滑动分布反演中，一般采用数据拟合残差和模型粗糙度之间的折中曲线来确定(由于其形似“L”，故称其为L曲线法)，但L曲线法在确定正则化参数过程中具有以下缺点：无法从L曲线图中根据横坐标轴直观读取正则化参数的大小；L曲线法过于依赖数据拟合精度[8]，求解过程可能不收敛[9]。

针对以上问题，本文提出利用U曲线法确定地震滑动分布反演的正则化参数，通过模拟实验以及芦山实际震例反演，验证U曲线法确定地震滑动分布反演正则化参数的可行性与优势。 1 U曲线法确定正则化参数 1.1 地震滑动分布反演基本方程

地震滑动分布反演中，地震滑动分布位移与地震滑动参数之间关系可表达为： d=Gm+ε (1)

式中，d表示用InSAR观测到的LOS向形变量或GPS观测到的3个方向的形变量；G为系数格林矩阵；m为地震滑动参数；ε为InSAR数据或GPS数据观测误差。

在地震滑动分布反演中，为避免系数格林矩阵的秩亏性同时保证断层面各断层块间的平滑度，本文采用Tikhonov正则化解法求解滑动分布，并引入目前较为成功的用于断层面平滑约束的二阶拉普拉斯平滑算子进行约束，使得： Hm=0 (2)

式中，0视为虚拟观测矩阵，H为拉普拉斯平滑矩阵，R=HTH视为正则化矩阵。 根据Tikhonov提出的正则化理论，得到估计准则为: ‖d-Gm‖2+αΩ(m)= ‖d-Gm‖2+αmTRm=min (3)

式中，Ω(m)为稳定泛函；‖·‖为·的二范数；α为地震滑动分布反演时的正则化参数，在本文中采用U曲线法确定正则化参数α。 1.2 U曲线法基本原理

由式(1)、(2)可得地震滑动参数解为： m=(GTpG+αR)-1GTpd (4)

式(4)表明，地震滑动参数解的大小随正则化参数α的变化而变化，因此选择一个合理的正则化参数在地震滑动分布反演中尤为重要。

U曲线法确定正则化参数与L曲线法类似。L曲线法是根据不同的α值分别得到一组‖d-Gm‖2、‖Hm‖2值，以‖d-Gm‖2为横坐标、‖Hm‖2为纵坐标拟合成一条类似“L”形状的曲线，取L曲线上拐点附近对应的α值作为最优正则化参数[10]。

可以看出，L曲线法的精度依赖于‖d-Gm‖2与‖Hm‖2数据之间的拟合程度。U曲线法则是通过定义U(α)函数，根据α的取值得到U(α)-α曲线，曲线左侧近似垂直部分曲率最大点所确定的α值即为最优正则化参数[11-12]。U曲线法定义如下： (5)

式中，其中u为系数格林矩阵G经过奇异值分解后的左侧奇异值矩阵，

σ1≥σ2…≥σr>0为系数格林矩阵G经过奇异值分解后r个从大到小排列的奇异值，β1≥β2…≥βr>0为拉普拉斯二阶平滑矩阵经过奇异值分解后r个从大到小排列的奇异值。 1.3 两种方法比较

由L曲线定义可以得出，L曲线在确定正则化参数的过程中需要依赖于数据拟合残差‖d-Gm‖2以及模型粗糙度‖Hm‖2之间的数据拟合程度，且从L曲线图上无法直观地从横坐标获得最优正则化参数的大小。由文献[11]及U曲线定义可以得出以下结论：U曲线在区间内严格递减，在区间严格递增，此处σ1≥σ2…≥σr>0为系数格林矩阵G经过奇异值分解后r个从大到小排列的奇异值；U曲线在区间存在局部最小值，并且当系数格林矩阵G只含一个非零奇异值时，U(α)函数只含有一个极小值；,。

综合以上分析可以看出，利用U曲线法确定正则化参数较L曲线法具有无需依赖数据拟合度，不会出现因曲线发散而无法确定正则化参数的情况[7]，可以从U曲线图中直观读取正则化参数大小等优势。另外，在同震滑动分布反演中，观测值与位错模型反演值之间的均方根误差(RMS)大小常常用作评判滑动分布反演精度，其可表达为：

(6)

式中，n表示观测值个数，d和c分别表示观测值和模型反演值所构成的矩阵，P表示观测值权阵。 2 模拟实验

在模拟实验中，断层面几何参数如下：断层面几何中心位置X=0 km，Y=0 km，断层中心深度Depth=10 km，长度L=42 km，宽度W=42 km，走向角strike=55°，倾角dip=40°，滑动角φ=60°。根据点位与震中距离设置形变点密度模拟InSAR数据592个如图1所示，并给InSAR形变点施加观测误差N(0,32mm2)。为了使观测数据更贴近实际，本文采取拉奎拉实际地震[13]InSAR观测时的升轨和降轨迹线的方位角以及雷达波入射角，确定InSAR数据LOS向形变观测值与地表形变分量的转换系数。

图1 模拟InSAR数据形变观测点Fig.1 Simulating the deformation observation points of InSAR data

在模拟实验中，本文分别采用L曲线法、U曲线法确定反演过程中的正则化参数(见图2)，其中，L曲线法得到的正则化参数为α1=0.145，U曲线法确定的正则化参数为α2=0.055。两种方法反演地震滑动分布结果及残差分布如图3所示，其中，图3(a)、(b)分别为利用L曲线法、U曲线法反演的滑动分布结果；图3(c)、(d)为两种方法反演的滑动分布结果与真值的差值分布。

在模拟实验中，利用L曲线法、U曲线法反演地震滑动分布各参数结果见表1；部分参数与真值的差值以及U曲线法反演结果较L曲线法改进程度见表2。 图2 用L曲线法、U曲线法确定正则化参数Fig.2 The regularization parameters obtained by using L-curve and U-curve method

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

拟实验反演地震滑动分布结果Tab.1 The inversion of seismic slip distribution results in simulation experiments

方法平滑因子最大滑动量/m平均滑动量/m平均滑动角/(°)矩震级MW地震矩/1019N·m均方根误差/mm预设参数-1.381 40.209 060.000 06.662 41.105 8-L曲线法0.1451.1 60.226 359.073 96.685 61.197 80.536U曲线法0.0551.239 60.218 660.4 96.675 51.156 70.176

表2 模拟实验两种方法反演各参数结果与真值的比较Tab.2 Comparison the results of each parameters and true value of two methods in simulation experiments最大滑动量/m平均滑动量/m平均滑动角/(°)矩震级MW均方根误差/mmL曲线法0.226 80.017 30.926 40.023 20.536U曲线法0.141 80.009 60.4 90.013 10.176L-U*7.36%3.40%0.78%0.15%67.16%

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

从表1、表2可以看出，用U曲线法确定正则化参数反演模拟实验各参数结果较L曲线法反演结果更接近预设参数(真值)，在最大滑动量、平均滑动量上U曲线法反演结果优势更明显。另外，在地震滑动分布反演中均方根误差常常作为评判反演精度的指标，均方根误差越小，说明反演参数大小与真值越接近，反演精度越高。两种方法反演结果的均方根误差大小可以进一步说明U曲线法反演地震滑动分布精度要优于L曲线法。 3 芦山实际震例反演

2013-04-20 08：02 四川芦山发生M7.0地震，震源深度13 km。地震发生后，不同研究机构对此次地震进行快速反演分析，研究地震震源机制和反演地震震源参数，但基于不同的手段所反演出的震源参数结果具有一定的差异。本文采用文献[14]获取的GPS三维形变数据来约束地表形变场，利用文献[1]通过多峰值颗粒群

算法获得的芦山地震断层参数作为反演芦山地震的断层参数。在此基础上将断层破裂面沿断层走向、倾向把断层长度、宽度均扩展至51 km，并将破裂面延伸至地表，将断层面均匀剖分成1.5 km×1.5 km大小的矩形单元，共得到了1 156个矩形断层单元。

分别采用L曲线法、U曲线法确定反演的正则化参数(见图4)，其大小分别为α1=0.198，α2=0.063。

两种方法反演的地震滑动分布结果及残差分布如图5所示，其中，图5(a)、(b)分别为利用L曲线法、U曲线法反演的滑动分布结果；图5(c)为两种方法反演的滑动分布结果差值分布。

两种方法反演芦山地震滑动分布各参数结果见表3；近年来部分学者研究芦山地震各参数结果见表4。

图4 用L曲线法、U曲线法确定正则化参数Fig.4 The regularization parameters obtained by using L-curve and U-curve method

图5 芦山地震滑动分布反演结果及两种方法反演残差分布Fig.5 Slip distribution results of Lushan earthquake and the residuals of two methods表3 用L曲线法、U曲线法确定平滑因子反演芦山地震滑动分布结果Tab.3 The results of slip distribution inversion of Lushan earthquake by L and U curve methods 方法平滑因子最大滑动量/m平均滑动量/m最大值处*的滑动角/(°)矩震级MW地震矩/1019N·m均方根误差/mmL曲线法0.1980.4 30.136 171.279 06.650 71.061 93.6U曲线法0.0630.679 70.159 274.510 16.694 21.242 63.1 注：*最大值处的滑动角是指滑动量最大值处对应的滑动角。

表4 芦山地震断层几何及滑动参数Tab.4 The various source parameters for Lushan earthquake来源走向/(°)倾角/(°)滑动角/(°)经度/(°)纬度/(°)最大滑动量/m长度/m地震矩/1019N·m矩震级MW均匀滑动

206.4744.1180.13102.93830.20.7521.940.726.L曲线法206.4744.1171.28102.93830.20.51.001.066.65U曲线法206.4744.1174.51102.93830.20.6851.001.246.69Jiang等[14]208.0043.0071.00102.93830.2950.6146.000.956.60刘云华等

[15]210.8040.20-103.06830.2380.9431.200.536.45王卫民等[16]205.0038.50---1.59-1.6.70USGS216.0047.0093.00102.88830.308---6.50CENC220.0035.0095.00103.00030.300----

由图5(a)、(b)可以看出，本文两种方法反演的芦山地震结果均显示芦山地震主震区分布在地下4～20 km深度范围内。由表3可以看出，由于正则化参数选取方法不同，使反演芦山地震震源各参数结果具有一定的差异。表4给出了近年来多位学者基于不同数据、不同方法反演的芦山地震各震源参数结果。综合表3、表4可以看出，本文确定的断层位置与文献[14]确定的位置以及CENC、USGS机构利用地震波确定的结果较一致；对于最大滑动量，本文两种方法反演结果分别为0. m、0.68 m，略小于文献[1]、文献[15]反演的结果。对于最大滑动量，文献[16]综合地震震源机制及断层的地质构造资料确定出断层的位置、走向、倾角等参数，反演出最大滑动量为1.59 m；文献[14]反演得到的最大滑动量为0.61 m，结果也小于其通过格网搜索法均匀反演得到的0.70 m。最大滑动量差别比较明显，主要的原因可能与约束的地震数据有关。对于矩震级参数，本文两种方法反演的矩震级分别为MW6.65、MW6.69，对应的矩张量分别为1.06×1019 N·m、1.24×1019 N·m，这两种方法反演结果均在其他学者研究范围内。均匀滑动反演的矩震级和矩张量分别为MW6.和0.72×1019 N·m。均匀滑动将断层面看作均匀整体，所以可能遗漏了细节滑动，导致整体矩张量偏小。文献[16]反演的矩震级和矩张量分别为MW6.7和1.×1019 N·m，相对其他几位学者反演结果更大；对于均方根误差，U曲线法反演结果要小于L曲线法。由以上分析可以看出，本

文两种方法反演的芦山地震滑动分布结果均在其他学者研究范围内；由于原始数据来源、研究方法以及预设的断层参数的不同，导致最后反演得出的各断层参数具有一定的差异。综上，利用U曲线法确定实际地震正则化参数具有一定的可行性与优势。

对于本文所做的模拟实验以及芦山实际地震实验，可以看出U曲线法确定正则化参数较L曲线法的优势为：1)对于反演精度，从模拟实验中可以看出，利用U曲线法确定正则化参数反演的滑动分布各参数结果与真值的差值均小于L曲线法，即在模拟实验中U曲线法确定的正则化参数反演精度要高于L曲线法；从模拟实验以及芦山实际震例反演的均方根误差参数可以看出，利用U曲线法反演结果的均方根误差均小于L曲线法，其原因是利用U曲线法确定的正则化参数要小于L曲线法，由图2(b)、图4(b)及式(6)可知，正则化参数越小，对应的拟合残差越小，使反演的均方根误差也相对变小。但正则化参数α的大小也不能无限变小，只能在一定的范围内，否则无法起到正则化的作用，不能削弱系数格林矩阵病态问题，从而使反演结果不稳定。从本文模拟实验以及实际震例实验可以看出，U曲线法确定的平滑因子能很好地控制在α的合适范围内，使反演结果平滑度相当，从而进一步表明在地震滑动分布反演中利用U曲线法确定正则化参数较L曲线法更为合适。2)从绘制的L曲线图与U曲线图可以看出，由于U曲线图横坐标即为平滑因子α的大小，故从U曲线图上可以直观地读取平滑因子的大小，而从L曲线图中却不能直观读取α大小；L曲线图绘制的时候需要依赖‖d-Gm‖2和‖Hm‖2之间的数据拟合精度，U曲线法则不需要。 4 结 语

本文提出的U曲线法为地震滑动分布反演中正则化参数的确定增加了一条新的途径。本文通过模拟实验以及芦山实际地震反演表明，U曲线法较L曲线法确定正则化参数具有反演精度高、U曲线图横坐标直观反映正则化参数大小、无需依赖于

数据拟合度等优势。另外，对于U曲线法确定的正则化参数要小于L曲线法，这可能与L曲线法过度依赖于‖d-Gm‖2与‖Hm‖2之间的数据拟合程度，使绘制的L曲线有时收敛性差有关，至于具体原因有待于进一步探讨与研究。 参考文献

【相关文献】

[1] 李海燕. 震源位错模型参数反演方法研究[D]. 南昌：东华理工大学,2016(Li Haiyan. Research on Inversion Method for Hypocenter Parameter Base on Dislocation Model[D]. Nanchang：East China University of Technology, 2016)

[2] Tikhonov A N. Regularization of Incorrectly Posed Problems[J]. Sov Math, 1963, 4(1):1 624-1 627

[3] Golub G H, Heath M, Wahba G. Generalized Cross-Validation as a Method for Choosing a Good Ridge Parameter[J]. Technometrics, 1979, 21(2):215-223 [4] 王振杰, 欧吉坤. 用L-曲线法确定岭估计中的岭参数[J]. 武汉大学学报:信息科学版,

2004,29(3):235-238(Wang Zhenjie, Ou Jikun. Determination of Ridge Parameters in Ridge Estimation by L- Curve Method[J]. Geomatics and Information Science of Wuhan University, 2004, 29(3):235-238)

[5] Hansen P C. Analysis of Discrete Ⅲ-Posed Problems by Means of the L-Curve[J]. SIAM Review, 1992, 34(4):561-580

[6] Fan Q, Xu C, Yi L, et al. Implication of Adaptive Smoothness Constraint and Helmert Variance Component Estimation in Seismic Slip Inversion[J]. Journal of Geodesy, 2017,91(10): 1 163-1 177

[7] 鲁洋为, 王振杰. 用U曲线法确定岭估计中的岭参数[J]. 导航定位学报, 2015(3):132-134(Lu Yangwei, Wang Zhenjie. Determination of Ridge Parameters in Ridge Estimate by U-Curve Method[J]. Journal of Navigation and Positioning, 2015(3): 132-134)

[8] 袁强强, 沈焕锋, 李平湘,等.自适应正则化多幅影像超分辨率重建[J]. 中国图像图形学报,2010, 15(12):1 720-1 727(Yuan Qiangqiang, Shen Huanfeng, Li Pingxiang, et al. Adaptive Regularization Multi Image Super Resolution Reconstruction[J]. Journal of Image and Graphics, 2010, 15(12):1 720-1 727)

[9] Xu P L. Truncated SVD Methods for Discrete Linear Ⅲ-Posed Problems[J]. Geophysical Journal International, 2010, 135(2):505-514

[10] 向阳, 于鹏, 陈晓,等. 大地电磁反演中改进的自适应正则化因子选取[J]. 同济大学学报:自然科学版, 2013, 41(9):1 429-1 434 (Xiang Yang, Yu Peng, Chen Xiao, et al. An Improved Adaptive Regularized Parameter Selection in Magnetotelluric Inversion[J]. Journal of Tongji University:Natural Science, 2013, 41(9):1 429-1 434)

[11] Krawczyk-Stańdo D, Rudnicki M. Regularization Parameter Selection in Discrete Ill-Posed Problems-The Use of the U-Curve[J]. International Journal of Applied Mathematics & Computer Science, 2007，17 (2):157-1

[12] Arnrich S, Brüer P, Klauck M, et al. Suitability of L- and U-Curve Methods for Calculating Reasonable Adsorption Energy Distributions from Nitrogen Adsorption Isotherms[J]. Adsorption Science & Technology, 2014, 32(7):521-534

[13] 王乐洋, 李海燕, 温扬茂,等. 地震同震滑动分布反演的总体最小二乘方法[J]. 测绘学报, 2017, 46(3):307-315(Wang Leyang, Li Haiyan, Wen Yangmao. Total Least Squares Method Inversion for Coseismic Slip Distribution[J]. Acta Geodaetica et Cartographica Sinica, 2017, 46(3): 307-315)

[14] Jiang Z S, Wang M, Wang Y, et al.GPS Constrained Coseismic Source and Slip

Distribution of the 2013 MW6.6 Lushan, China, Earthquake and Its Tectonic Implications[J]. Geophysical Research Letters, 2014, 41(2):407-413

[15] 刘云华, 汪驰升, 单新建,等. 芦山MS7.0级地震InSAR形变观测及震源参数反演[J]. 地球物理学报, 2014, 57(8):2 495-2 506(Liu Yunhua, Wang Chisheng, Shan Xinjian, et al. Result of SAR Differential Interferometry for the Coseismic Deformation and Source Parameter of the MS7.0 Lushan Earthquake[J]. Chinese Journal of Geophysics, 2014, 57(8):2 495-2 506) [16] 王卫民, 郝金来, 姚振兴. 2013年四川芦山地震震源破裂过程反演初步结果[J]. 地球物理学报, 2013, 56(4):1 412-1 417(Wang Weimin, Hao Jinlai, Yao Zhengxing. Preliminary Result for Rupture Process of 2013,Lushan Earthquake, Sichuan, China[J].Chinese Journal of Geophysics, 2013, 56 (4):1 412-1 417)