留言板

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

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

利用偏差改正的方差分量估计方法确定联合反演相对权比

王乐洋 谷旺旺 赵雄 许光煜 高华

王乐洋, 谷旺旺, 赵雄, 许光煜, 高华. 利用偏差改正的方差分量估计方法确定联合反演相对权比[J]. 武汉大学学报 ● 信息科学版, 2022, 47(4): 508-516. doi: 10.13203/j.whugis20200216
引用本文: 王乐洋, 谷旺旺, 赵雄, 许光煜, 高华. 利用偏差改正的方差分量估计方法确定联合反演相对权比[J]. 武汉大学学报 ● 信息科学版, 2022, 47(4): 508-516. doi: 10.13203/j.whugis20200216
WANG Leyang, GU Wangwang, ZHAO Xiong, XU Guangyu, GAO Hua. Determination of Relative Weight Ratio of Joint Inversion Using Bias-Corrected Variance Component Estimation Method[J]. Geomatics and Information Science of Wuhan University, 2022, 47(4): 508-516. doi: 10.13203/j.whugis20200216
Citation: WANG Leyang, GU Wangwang, ZHAO Xiong, XU Guangyu, GAO Hua. Determination of Relative Weight Ratio of Joint Inversion Using Bias-Corrected Variance Component Estimation Method[J]. Geomatics and Information Science of Wuhan University, 2022, 47(4): 508-516. doi: 10.13203/j.whugis20200216

利用偏差改正的方差分量估计方法确定联合反演相对权比

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

国家自然科学基金 41874001

国家自然科学基金 41664001

国家自然科学基金 42174011

详细信息
    作者简介:

    王乐洋, 博士, 教授, 主要研究方向为大地测量反演及大地测量数据处理。wleyang@163.com

  • 中图分类号: P223

Determination of Relative Weight Ratio of Joint Inversion Using Bias-Corrected Variance Component Estimation Method

Funds: 

The National Natural Science Foundation of China 41874001

The National Natural Science Foundation of China 41664001

The National Natural Science Foundation of China 42174011

More Information
    Author Bio:

    WANG Leyang, PhD, professor, majors in geodetic inversion and geodetic data processing.E-mail: wleyang@163.com

  • 摘要: 在大地测量联合反演中,方差分量估计法用于确定相对权比时并没有考虑大地测量反演的病态性,利用正则化解代替最小二乘解会引入偏差,会造成方差分量估计不准确问题。针对此提出采用偏差改正方差分量估计方法,以消除正则化解引入偏差的影响,并基于残差的偏差改正方差分量估计与方差分量估计法进行模拟实验。结果表明,进行偏差改正后的方差分量估计法能够较好地反演出滑动分布情况,所提方法针对参数进行偏差改正的方差分量估计考虑了迭代初值引入的偏差,理论更为严密。并将所提方法用于Visso地震和Norcia地震反演中,验证了该方法的可靠性、合理性。
  • 图  1  模拟地表观测点

    Figure  1.  Simulated Surface Observation Points

    图  2  模拟实验的滑移分布结果及其残差

    Figure  2.  Slip Distribution Results and Residuals in the Simulation Experiment

    图  3  3种方案反演的Visso地震滑动分布结果

    Figure  3.  Results of the Slip Distribution of the Visso Earthquake for Three Schemes

    图  4  3种方案反演的Norcia地震的滑动分布结果

    Figure  4.  Results of the Slip Distribution of the Norcia Earthquake for Three Schemes

    表  1  模拟实验的3种方案

    Table  1.   Three Schemes of Simulated Experiment

    方案 方法
    方案1 方差分量估计法确定相对权比和正则化参数
    方案2 MSE求解正则化参数,方差分量估计法确定相对权比,迭代求解
    方案3 偏差改正方差分量估计确定相对权比和正则化参数
    下载: 导出CSV

    表  2  模拟实验反演地震滑动分布结果

    Table  2.   Inversion of Seismic Slip Distribution Results in Simulation Experiment

    参数 预设参数值 方案1 方案2 方案3
    相对权比 1∶0.210 8∶0.001 5 1∶0.060 9∶0.001 1∶0.209 3∶0.000 9
    最大滑动量/m 1.444 8 1.398 0 1.409 9 1.421 4
    平均滑动量/m 0.249 2 0.266 6 0.266 8 0.268 1
    最大滑动角/(°) 50 51.065 5 50.961 5 51.041 2
    矩震级(Mw) 6.518 6 6.538 1 6.538 3 6.539 8
    地震矩/(1018 N·m) 6.727 5 7.197 6 7.202 4 7.238 9
    RMS/mm 6.3 6.3 6.3
    下载: 导出CSV

    表  3  Visso地震断层滑动参数

    Table  3.   Slip Parameters for the Visso Earthquake

    方法 相对权比 最大滑动量/m 平均滑动量/m 平均滑动角/(°) 矩震级(Mw) 地震矩/(1018 N·m) RMS/mm
    方案1 1∶2.073 4∶0.002 6 0.605 6 0.194 8 -82.481 5 6.216 2 2.367 3 7.2
    方案2 1∶0.000 83∶0.001 0 0.623 9 0.177 9 -68.742 8 6.189 9 2.161 8 16.2
    方案3 1∶2.888 1∶0.001 0 0.692 6 0.213 7 -85.693 2 6.242 9 2.596 4 6.8
    文献[20] 0.840 0 -67.400 0 6.190 0 2.150 0
    文献[22] 0.900 0 -89.700 0 6.140 0 1.810 0
    文献[23] 0.120 0 -81/-80 6.000 0 1.130 0
    文献[24] 6.160 0 1.970 0
    USGS -89 6.100 0 1.840 0
    GCMT -94 6.100 0 1.610 0
    下载: 导出CSV

    表  4  Norcia地震断层滑动参数

    Table  4.   Slip Parameters for the Norcia Earthquake

    方法 相对权比 最大滑动量/m 平均滑动量/m 平均滑动角/(°) 矩震级(Mw) 地震矩/(1018 N·m) RMS/mm
    方案1 1∶5.268 4∶0.001 1 2.013 6 0.695 2 -87.342 6 6.584 4 8.446 2 13.4
    方案2 1∶0.005 8∶0.01 2.484 4 0.732 6 -84.770 7 6.599 6 8.900 8 22.3
    方案3 1∶4.196 4∶0.000 39 2.267 3 0.707 0 -84.894 7 6.589 3 8.590 2 12.7
    文献[20] 3.440 0 -64.100 0 6.600 0 9.140 0
    文献[22] 2.500 0 6.600 0 8.970 0
    文献[23] 0.420 0 -90.000 0 6.400 0 4.430 0
    文献[24] 6.590 0 8.460 0
    USGS -84.000 0 6.600 0 10.700 0
    GCMT -96.000 0 6.600 0 9.580 0
    下载: 导出CSV
  • [1] 许才军. 大地测量联合反演理论和方法研究进展[J]. 武汉大学学报·信息科学版, 2001, 26(6): 555-561 http://ch.whu.edu.cn/article/id/5232

    Xu Caijun. Progress of Joint Inversion on Geodesy and Geophysics[J]. Geomatics and Information Science of Wuhan University, 2001, 26(6): 555-561 http://ch.whu.edu.cn/article/id/5232
    [2] Li Z C, Wen Y M, Zhang P, et al. Joint Inversion of GPS, Leveling, and InSAR Data for the 2013 Lushan (China) Earthquake and Its Seismic Hazard Implications[J]. Remote Sensing, 2020, 12(4): 715 doi:  10.3390/rs12040715
    [3] Yi L, Xu C J, Wen Y M, et al. Rupture Process of the 2016 Mw 7.8 Ecuador Earthquake from Joint Inversion of InSAR Data and Teleseismic P Waveforms[J]. Tectonophysics, 2018, 722: 163-174 doi:  10.1016/j.tecto.2017.10.028
    [4] Wang S, Xu C J, Wen Y M, et al. Slip Model for the 25 November 2016 Mw 6.6 Aketao Earthquake, Western China, Revealed by Sentinel-1 and ALOS-2 Observations[J]. Remote Sensing, 2017, 9(4): 325 doi:  10.3390/rs9040325
    [5] Yi L, Xu C J, Zhang X, et al. Joint Inversion of GPS, InSAR and Teleseismic Data Sets for the Rupture Process of the 2015 Gorkha, Nepal, Earthquake Using a Generalized ABIC Method[J]. Journal of Asian Earth Sciences, 2017, 148: 121-130 doi:  10.1016/j.jseaes.2017.08.029
    [6] Xu C J, Ding K H, Cai J Q, et al. Methods of Determining Weight Scaling Factors for Geodetic-Geophysical Joint Inversion[J]. Journal of Geodynamics, 2009, 47(1): 39-46 doi:  10.1016/j.jog.2008.06.005
    [7] 方进, 许才军, 温扬茂, 等. 基于方差分量估计的2015年尼泊尔MW 7.8地震同震滑动分布[J]. 地球物理学报, 2019, 62(3): 923-939 https://www.cnki.com.cn/Article/CJFDTOTAL-DQWX201903009.htm

    Fang Jin, Xu Caijun, Wen Yangmao, et al. Coseismic Slip Distribution of 2015 Gorkha (Nepal) MW 7.8 Earthquake Determined Using the Helmert Variance Component Estimation[J]. Chinese Journal of Geophysics, 2019, 62(3): 923-939 https://www.cnki.com.cn/Article/CJFDTOTAL-DQWX201903009.htm
    [8] Fan Q B, Xu C J, 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): 1163-1177 doi:  10.1007/s00190-017-1015-0
    [9] 王乐洋, 温贵森. 火山Mogi模型形变反演的虚拟观测法及方差分量估计[J]. 地球物理学报, 2020, 63(5): 1775-1786 https://www.cnki.com.cn/Article/CJFDTOTAL-DQWX202005005.htm

    Wang Leyang, Wen Guisen. Virtual Observation Method and Variance Component Estimation for Deformation Inversion on the Volcano Mogi Model[J]. Chinese Journal of Geophysics, 2020, 63(5): 1775-1786 https://www.cnki.com.cn/Article/CJFDTOTAL-DQWX202005005.htm
    [10] 王乐洋, 赵雄, 高华. 大地测量地震断层同震滑动分布反演的两步解法[J]. 武汉大学学报·信息科学版, 2019, 44(9): 1265-1273 doi:  10.13203/j.whugis20170382

    Wang Leyang, Zhao Xiong, Gao Hua. A Two-Step Solution Method for the Co-seismic Slip Distribution Inversion of Earthquake Faults in Geodesy[J]. Geomatics and Information Science of Wuhan University, 2019, 44(9): 1265-1273 doi:  10.13203/j.whugis20170382
    [11] Xu P L, Shen Y Z, Fukuda Y, et al. Variance Component Estimation in Linear Inverse Ill-Posed Models[J]. Journal of Geodesy, 2006, 80(2): 69-81 doi:  10.1007/s00190-006-0032-1
    [12] Xu P L. Iterative Generalized Cross-Validation for Fusing Heteroscedastic Data of Inverse Ill-Posed Problems[J]. Geophysical Journal International, 2009, 179(1): 182-200 doi:  10.1111/j.1365-246X.2009.04280.x
    [13] Wang L Y, Zhao X, Gao H. A Method for Determining the Regularization Parameter and the Relative Weight Ratio of the Seismic Slip Distribution with Multi-source Data[J]. Journal of Geodynamics, 2018, 118: 1-10 doi:  10.1016/j.jog.2018.04.005
    [14] 王乐洋, 许光煜, 陈晓勇. 附有相对权比的PEIV模型总体最小二乘平差[J]. 武汉大学学报·信息科学版, 2017, 42(6): 857-863 doi:  10.13203/j.whugis20150001

    Wang Leyang, Xu Guangyu, Chen Xiaoyong. Total Least Squares Adjustment of Partial Errors-in-Variables Model with Weight Scaling Factor[J]. Geomatics and Information Science of Wuhan University, 2017, 42(6): 857-863 doi:  10.13203/j.whugis20150001
    [15] 王乐洋, 余航. 附有相对权比的加权总体最小二乘联合平差方法[J]. 武汉大学学报·信息科学版, 2019, 44(8): 1233-1240 doi:  10.13203/j.whugis20170265

    Wang Leyang, Yu Hang. Weighted Total Least Squares Method for Joint Adjustment Model with Weight Scaling Factor[J]. Geomatics and Information Science of Wuhan University, 2019, 44(8): 1233-1240 doi:  10.13203/j.whugis20170265
    [16] 许才军, 邓长勇, 周力璇. 利用方差分量估计的地震同震滑动分布反演[J]. 武汉大学学报·信息科学版, 2016, 41(1): 37-44 https://www.cnki.com.cn/Article/CJFDTOTAL-WHCH201601007.htm

    Xu Caijun, Deng Changyong, Zhou Lixuan. Coseismic Slip Distribution Inversion Method Based on the Variance Component Estimation[J]. Geomatics and Information Science of Wuhan University, 2016, 41(1): 37-44 https://www.cnki.com.cn/Article/CJFDTOTAL-WHCH201601007.htm
    [17] 王乐洋, 李海燕, 陈汉清. 2013年芦山Ms 7.0级地震断层参数模型反演[J]. 武汉大学学报·信息科学版, 2019, 44(3): 347-354 doi:  10.13203/j.whugis20160485

    Wang Leyang, Li Haiyan, Chen Hanqing. Source Parameters and Slip Distribution Inversion of 2013 Lushan Ms 7.0 Earthquake[J]. Geomatics and Information Science of Wuhan University, 2019, 44(3): 347-354 doi:  10.13203/j.whugis20160485
    [18] Wang L Y, Gu W W. An Optimal Design Method to Determine the Regularization Parameter of Coseismic Slip Distribution Inversion[J]. Geophysical Journal International, 2020, 221(1): 440-450 doi:  10.1093/gji/ggaa030
    [19] Okada Y. Surface Deformation Due to Shear and Tensile Faults in a Half-Space[J]. Bulletin of the Seismological Society of America, 1985, 75(4): 1135-1154 doi:  10.1785/BSSA0750041135
    [20] Wang L Y, Gao H, Feng G C, et al. Source Parameters and Triggering Links of the Earthquake Sequence in Central Italy from 2009 to 2016 Analyzed with GPS and InSAR Data[J]. Tectonophysics, 2018, 744: 285-295 doi:  10.1016/j.tecto.2018.07.013
    [21] 冯万鹏, 李振洪. InSAR资料约束下震源参数的PSO混合算法反演策略[J]. 地球物理学进展, 2010, 25(4): 1189-1196 doi:  10.3969/j.issn.1004-2903.2010.04.007

    Feng Wanpeng, Li Zhenhong. A Novel Hybrid PSO/Simplex Algorithm for Determining Earthquake Source Parameters Using InSAR Data[J]. Progress in Geophysics, 2010, 25(4): 1189-1196 doi:  10.3969/j.issn.1004-2903.2010.04.007
    [22] Xu G Y, Xu C J, Wen Y M, et al. Source Parameters of the 2016—2017 Central Italy Earthquake Sequence from the Sentinel-1, ALOS-2 and GPS Data[J]. Remote Sensing, 2017, 9(11): 1182 doi:  10.3390/rs9111182
    [23] Chiaraluce L, di Stefano R, Tinti E, et al. The 2016 Central Italy Seismic Sequence: A First Look at the Mainshocks, Aftershocks, and Source Models[J]. Seismological Research Letters, 2017, 88(3): 757-771 doi:  10.1785/0220160221
    [24] Cheloni D, de Novellis V, Albano M, et al. Geodetic Model of the 2016 Central Italy Earthquake Sequence Inferred from InSAR and GPS Data[J]. Geophysical Research Letters, 2017, 44(13): 6778-6787 doi:  10.1002/2017GL073580
  • [1] 吕志鹏, 隋立芬.  基于变量投影法的自回归模型方差分量估计 . 武汉大学学报 ● 信息科学版, 2020, 45(2): 205-212. doi: 10.13203/j.whugis20180352
    [2] 王乐洋, 余航.  附有相对权比的加权总体最小二乘联合平差方法 . 武汉大学学报 ● 信息科学版, 2019, 44(8): 1233-1240. doi: 10.13203/j.whugis20170265
    [3] 王乐洋, 余航.  火山Mogi模型反演的总体最小二乘联合平差方法 . 武汉大学学报 ● 信息科学版, 2018, 43(9): 1333-1341. doi: 10.13203/j.whugis20160469
    [4] 王乐洋, 许光煜, 陈晓勇.  附有相对权比的PEIV模型总体最小二乘平差 . 武汉大学学报 ● 信息科学版, 2017, 42(6): 857-863. doi: 10.13203/j.whugis20150001
    [5] 王乐洋, 余航.  总体最小二乘联合平差 . 武汉大学学报 ● 信息科学版, 2016, 41(12): 1683-1689. doi: 10.13203/j.whugis20140670
    [6] 许才军, 邓长勇, 周力璇.  利用方差分量估计的地震同震滑动分布反演 . 武汉大学学报 ● 信息科学版, 2016, 41(1): 37-44. doi: 10.13203/j.whugis20150500
    [7] 张朝玉, 许才军.  具有自适应权比的大地测量联合反演序贯算法及其应用 . 武汉大学学报 ● 信息科学版, 2012, 37(10): 1140-1144.
    [8] 王乐洋, 许才军.  附有相对权比的总体最小二乘平差 . 武汉大学学报 ● 信息科学版, 2011, 36(8): 887-890.
    [9] 许才军, 王乐洋.  大地测量和地震数据联合反演地震震源破裂过程研究进展 . 武汉大学学报 ● 信息科学版, 2010, 35(4): 457-462.
    [10] 丁开华, 许才军.  川滇地区地壳应变场的GPS与地震矩张量联合反演研究 . 武汉大学学报 ● 信息科学版, 2009, 34(3): 265-268.
    [11] 温扬茂, 许才军.  联合GPS与重力资料反演分析川滇地区现今地壳形变 . 武汉大学学报 ● 信息科学版, 2009, 34(5): 568-572.
    [12] 徐天河, 居向明.  基于方差分量估计的CHAMP重力场恢复方法 . 武汉大学学报 ● 信息科学版, 2007, 32(3): 242-245.
    [13] 王新洲, 游扬声, 刘星, 史文中.  GIS叠置图层方差分量的极大似然估计 . 武汉大学学报 ● 信息科学版, 2005, 30(10): 892-896.
    [14] 杨元喜, 徐天河.  基于移动开窗法协方差估计和方差分量估计的自适应滤波 . 武汉大学学报 ● 信息科学版, 2003, 28(6): 714-718.
    [15] 独知行, 刘经南.  利用GPS位移和主应力方向观测资料进行川滇地区边界力的联合反演研究 . 武汉大学学报 ● 信息科学版, 2003, 28(2): 162-166,176.
    [16] 王仲锋.  导线网方差分量估计的综合研究 . 武汉大学学报 ● 信息科学版, 2001, 26(2): 112-117.
    [17] 盛乐山.  三维网平差的方差分量估计 . 武汉大学学报 ● 信息科学版, 1990, 15(1): 91-98.
    [18] 许才军.  卫星网与地面网联合平差中随机模型的估计 . 武汉大学学报 ● 信息科学版, 1989, 14(3): 35-45.
    [19] 吴晓清.  权因子法方差分量估计 . 武汉大学学报 ● 信息科学版, 1989, 14(4): 8-15.
    [20] 张万鹏.  方差分量估计理论及其在边角网平差中的应用 . 武汉大学学报 ● 信息科学版, 1988, 13(1): 9-21.
  • 加载中
图(4) / 表(4)
计量
  • 文章访问数:  549
  • HTML全文浏览量:  72
  • PDF下载量:  47
  • 被引次数: 0
出版历程
  • 收稿日期:  2020-05-07
  • 刊出日期:  2022-04-05

利用偏差改正的方差分量估计方法确定联合反演相对权比

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

    国家自然科学基金 41874001

    国家自然科学基金 41664001

    国家自然科学基金 42174011

    作者简介:

    王乐洋, 博士, 教授, 主要研究方向为大地测量反演及大地测量数据处理。wleyang@163.com

  • 中图分类号: P223

摘要: 在大地测量联合反演中,方差分量估计法用于确定相对权比时并没有考虑大地测量反演的病态性,利用正则化解代替最小二乘解会引入偏差,会造成方差分量估计不准确问题。针对此提出采用偏差改正方差分量估计方法,以消除正则化解引入偏差的影响,并基于残差的偏差改正方差分量估计与方差分量估计法进行模拟实验。结果表明,进行偏差改正后的方差分量估计法能够较好地反演出滑动分布情况,所提方法针对参数进行偏差改正的方差分量估计考虑了迭代初值引入的偏差,理论更为严密。并将所提方法用于Visso地震和Norcia地震反演中,验证了该方法的可靠性、合理性。

English Abstract

王乐洋, 谷旺旺, 赵雄, 许光煜, 高华. 利用偏差改正的方差分量估计方法确定联合反演相对权比[J]. 武汉大学学报 ● 信息科学版, 2022, 47(4): 508-516. doi: 10.13203/j.whugis20200216
引用本文: 王乐洋, 谷旺旺, 赵雄, 许光煜, 高华. 利用偏差改正的方差分量估计方法确定联合反演相对权比[J]. 武汉大学学报 ● 信息科学版, 2022, 47(4): 508-516. doi: 10.13203/j.whugis20200216
WANG Leyang, GU Wangwang, ZHAO Xiong, XU Guangyu, GAO Hua. Determination of Relative Weight Ratio of Joint Inversion Using Bias-Corrected Variance Component Estimation Method[J]. Geomatics and Information Science of Wuhan University, 2022, 47(4): 508-516. doi: 10.13203/j.whugis20200216
Citation: WANG Leyang, GU Wangwang, ZHAO Xiong, XU Guangyu, GAO Hua. Determination of Relative Weight Ratio of Joint Inversion Using Bias-Corrected Variance Component Estimation Method[J]. Geomatics and Information Science of Wuhan University, 2022, 47(4): 508-516. doi: 10.13203/j.whugis20200216
  • 大地测量反演问题是利用大地测量手段探测地震致灾机理的关键。随着观测技术的发展,用于反演的数据类型逐渐丰富,其中全球定位系统(global positioning system,GPS)数据具有高时间分辨率,而合成孔径干涉雷达(interferometric synthetic aperture radar,InSAR)数据具有高空间分辨率,为了提高反演的精度,通常联合GPS和InSAR两类数据进行大地测量联合反演[1-5]。而大地测量联合反演的关键在于相对权比的确定,数据的权重体现了该类数据对于反演的贡献程度,因此一个相对合理的权比能够提高反演的精度。

    文献[6]总结了目前较为常用的定权方法,将相对权比作为待求参数统一求解,根据先验信息确定相对权比,当相对权比为1∶1时进行等权处理,并采用赫尔默特方差分量估计法[7-9]定权;文献[10]提出了确定联合反演中相对权比的两步法,并讨论了不同的目标函数对相对权比的影响;文献[11]研究了病态问题中方差分量估计的适用性,考虑到参数的偏差可能会大于参数本身,提出了对残差进行偏差改正的方差分量估计来确定相对权比,并利用均方误差(mean square error,MSE)方法确定正则化参数,通过迭代计算的方式,减弱正则化参数和方差分量之间的相互影响,准确估计出方差分量与正则化参数;文献[12]采用广义交叉验证法确定正则化参数,再利用方差分量估计法确定相对权比,并通过迭代的方式确定最优权比与正则化参数;文献[13]采用判别函数最小化法[14-15]确定联合反演中的相对权比;文献[16]提出利用方差分量估计法同时确定相对权比与正则化参数,拓展了方差分量估计法的应用范围。在大地测量联合反演中,方差分量估计法仍然是目前最为常用的确定相对权比的方法。而在确定大地测量联合反演相对权比的过程中,通常没有考虑大地测量反演问题的病态性,直接利用正则化解代替最小二乘解进行方差分量估计,由于正则化解引入了偏差,此时正则化解的期望值并不为零。因此,针对以上问题,本文采用偏差改正方差分量估计方法同时确定联合反演的相对权比和正则化参数。

    • 在同震滑动分布反演中,地表同震形变量与断层滑动量之间为线性关系[17-18],可表示为:

      d=Gm+ε ]]>

      式中,d表示观测向量(如InSAR视线向位移、GPS位移场等);G表示格林函数矩阵;m为滑动参数;ε为InSAR数据或GPS数据等观测值的随机观测误差。

      对于多类相互独立的地表形变观测数据,则式(1)可以表示为:

      d1=G1m+ε1d2=G2m+ε2dn=Gnm+εn ]]>

      在同震滑动分布反演中,为了得到更精细的滑动,通常将断层面划分为n个矩形单元,基于矩形位错模型计算矩形单元的单位走滑量和单位倾滑量引起地表k个观测点的形变值,构成k×2n的格林函数矩阵G。断层面离散化处理后,会导致G严重秩亏,此时式(1)的法矩阵病态,从而使滑动分布解不稳定。为了避免这种情况,可对断层的位错施加一定的约束,一方面增加虚拟观测数据避免系数矩阵的秩亏性,另一方面抑制相邻断层单元出现大的梯度变化[17]

      本文采用拉普拉斯约束对断层单元进行平滑约束,使:

      Hm=0 ]]>

      式中,H为拉普拉斯二阶平滑矩阵。联合式(2)和式(3),得到滑动分布联合反演公式为:

      d1d2dn0T= G1G2GnHTm+ε ]]>
      P=diag1σ12D1-11σ22D2-11σn2Dn-11α2I ]]>

      式中,P为联合反演的权阵;Di为各类数据的方差阵;1/σ121/σ221/σn2表示各类观测数据的相对权重;1/α2表示模型平滑度与数据拟合度之间的相对权比。式(5)的最小二乘解为:

      mα=(G1TP1G1+G2TP2G2++GnTPnGn+HTPn+1H)-1(G1TP1d1+ G2TP2d2++GnTPndn) ]]>

      采用方差分量估计法对式(5)进行迭代计算,可以求解各类观测值的方差分量,从而确定联合反演的相对权比与正则化参数。然而在滑动分布反演中,利用方差分量估计方法确定相对权比时需要初始值进行迭代,而未加入平滑约束的滑动分布解是不可靠的,故将有偏的滑动分布解作为初始值进行迭代运算。而此时需要考虑其偏差的影响,且已有文献并没考虑滑动分布解是有偏的,因此本文利用偏差改正的方差分量估计进行求解。这里仅列出GPS数据与InSAR数据联合反演的公式,误差方程可表示为:

      v1=G1m-d1v2=G2m-d2v3=Hm-0 ]]>

      所求方差为:

      D(v1)=(G1N-1N1N-1G1T-2G1N-1G1T+P1-1)σ12 +G1N-1N2N-1G1Tσ22+G1N-1N3N-1G1Tσ32 ]]>

      其中,

      N1=G1TP1G1N2=G2TP2G2N3=HTP3H]]>
      N=N1+N2+N3]]>

      根据方差分量估计公式,可得:

      E(v1TP1v1)=(n1+tr(N-1N1N-1N1)-2tr(N-1N1))σ12+tr(N-1N2N-1N1)σ22+tr(N-1N3N-1N1)σ32 ]]>

      式中,E()表示期望运算;tr表示矩阵的求迹运算。

      同理可得到E(v2TP2v2)E(v3TP3v3)的表达式,将数学期望的符号去掉,则可以求出方差分量的估值,其矩阵形式为:

      Sθ̂=Wθ ]]>

      式中,

      θ̂=σ̂12σ̂22α2T]]>
      Wθ=v1TP1v1v2TP2v2v3TP3v3T]]>
      S=Atr(N-1N2N-1N1)tr(N-1N3N-1N1)Btr(N-1N3N-1N2)C]]>

      其中,A=n1+tr(N-1N1N-1N1)-2tr(N-1N1)B=n2+tr(N-1N2N-1N2)-2tr(N-1N2)C=n3+tr(N-1N3N-1N3)-2tr(N-1N3)

      滑动分布解的偏差为:

      β=-N-1N3m̂ ]]>

      式中,m̂为无偏估计值,由于其无法获取,故采用正则化解替代,从而得到滑动分布解的偏差计算如下:

      βα=-N-1N3mα ]]>

      由于每次迭代计算时,所代入的滑动分布解都是有偏的,因此需要进行偏差改正,改正后的滑动分布解公式为:

      m=mα-βα ]]>

      将偏差改正的滑动分布解作为初值进行迭代,每次迭代后都进行偏差改正,从而得到偏差改正的方差分量估计。

      文献[11]提出对参数的残差进行偏差改正,其中参数表达式为:

      mλ=(G1TP1G1+G2TP2G2+λHTH)-1 (G1TP1d1+G2TP2d2) ]]>

      得到参数残差的偏差为:

      E(v1)= -λG1Nλ-1RmλNλ-1=(G1TP1G1+G2TP2G2+λHTH) ]]>

      式中,λ表示正则化参数;R=HTH为正则化矩阵。

      根据方差分量估计公式,可得:

      E(v1TP1v1)=tr(P1D(v1))+ E(v1T)P1E(v1)=(n1-2tr(Nλ-1N1)+tr(Nλ-1N1Nλ-1N1))σ12+tr(Nλ-1N2Nλ-1N1)σ22+λ2mλTRTNλ-1N1Nλ-1Rmλ ]]>

      同理可得E(v2TP2v2)表达式。将式(16)数学期望的符号去掉,则可以求出方差分量的估值,得到其矩阵形式为:

      S1θ̂1+M1=W1 ]]>

      式中,

      S1=Dtr(Nλ-1N1Nλ-1N2)E]]>

      其中,D=n1-2tr(Nλ-1N1)+tr(Nλ-1N1)2E=n2-2tr(Nλ-1N2)+tr(Nλ-1N2)2

      θ̂1=σ12σ22T]]>
      M1=λ2mλTRTNλ-1N1Nλ-1Rmλλ2mλTRTNλ-1N2Nλ-1Rmλ]]>
      W1=v1TP1v1v2TP2v2T]]>

      考虑到正则化参数会放大方差分量,采用MSE法求解正则化参数,并通过迭代求解得到最后的正则化参数和方差分量,迭代步骤如下:(1)采用MSE法确定滑动分布反演的正则化参数;(2)利用获得的正则化参数求得滑动分布解的有偏估计量;(3)计算出残差的偏差,对残差进行偏差改正,利用偏差改正后方差分量估计求得各类数据的方差分量;(4)根据确定的方差分量得到具有不同权比的滑动分布解,以此作为初值重复步骤(1)~(3),直到均方误差收敛,迭代停止。

      为了更好地比较不同方差分量估计的适用性和合理性,采用3种方案进行模拟实验,如表 1所示。

      表 1  模拟实验的3种方案

      Table 1.  Three Schemes of Simulated Experiment

      方案 方法
      方案1 方差分量估计法确定相对权比和正则化参数
      方案2 MSE求解正则化参数,方差分量估计法确定相对权比,迭代求解
      方案3 偏差改正方差分量估计确定相对权比和正则化参数
    • 为了验证基于偏差改正方差的分量估计方法的合理性与优势,本文进行了模拟实验,模拟一个双主滑动区的断层进行滑动分布反演。设置断层几何中心位置为X=0 km,Y=0 km,长度为30 km,宽度为30 km,走向角为50°,倾角为45°,滑动角为50°,断层顶深为0 km。首先模拟了断层的真实滑动分布情况,然后利用均匀弹性半空间地表位移与矩形位错公式[19]计算了每个断层单元的滑动对地表417个站台的地表位移之和,并给每个台站的位移施加观测误差N(0,32 mm2),最终得到模拟形变数据。采用相同的方法获得592个InSAR形变位移,考虑到InSAR数据质量较差,对InSAR形变位移施加误差N(0,12 cm2),得到InSAR模拟形变数据。GPS与InSAR形变数据如图 1所示。在滑动分布反演中,将断层离散成20×20个1.5 km×1.5 km的断层单元,并利用偏差改正方差分量估计进行滑动分布反演。

      图  1  模拟地表观测点

      Figure 1.  Simulated Surface Observation Points

      在模拟实验中,分别采用表 1中的3种方案确定多类数据的相对权比,方案1、2、3确定的GPS、InSAR、正则化参数的相对权比分别为1∶0.210 8∶0.001 5、1∶0.060 9∶0.001、1∶0.209 3∶0.000 9。图2(a)2(b)2(c)分别为采用方案1、2、3反演得到的滑动分布结果;图2(d)2(e)2(f)分别为方案1、2、3反演得到的滑动分布结果与真实滑移的差值分布。

      图  2  模拟实验的滑移分布结果及其残差

      Figure 2.  Slip Distribution Results and Residuals in the Simulation Experiment

      图 2可以看出,3种方案都能很好地反演出两个主滑动区,浅部滑动区的滑动能够较好地反演出来,而深部滑动区滑动量存在明显的低估,这与数据难以分辨出地下深层的滑动有关。对于两个主滑动区之间的区域,滑动量均被高估,方案3较方案2更接近模拟滑移,方案1结果最差;对于深部滑动区,方案1难以较好地反演出实际滑移,而方案2较方案1有所改进,方案3相比方案1、2,反演结果更接近于模拟滑动情况。

      表 2统计了模拟实验中3种方案确定的相对权比以及反演得到的各参数结果和模拟值。从表 2可以看出,方案1反演的最大滑动量最小,方案2反演的最大滑动量要略高于方案1,方案3结果略高于方案2,与滑动分布图结果一致。对于最大滑动角,3个方案反演结果差别不大,而方案2结果更接近实际滑动角。对比3个方案反演的地震矩,方案1至方案3依次增大,这与最大滑动量情况一致。均方根(root mean square,RMS)误差作为大地测量反演中的一种衡量反演结果指标,3种方案反演结果均为6.3 mm,从RMS角度分析,3个方案的反演精度相当,但方案2和方案3考虑了反演中偏差的影响,理论更为严密。

      表 2  模拟实验反演地震滑动分布结果

      Table 2.  Inversion of Seismic Slip Distribution Results in Simulation Experiment

      参数 预设参数值 方案1 方案2 方案3
      相对权比 1∶0.210 8∶0.001 5 1∶0.060 9∶0.001 1∶0.209 3∶0.000 9
      最大滑动量/m 1.444 8 1.398 0 1.409 9 1.421 4
      平均滑动量/m 0.249 2 0.266 6 0.266 8 0.268 1
      最大滑动角/(°) 50 51.065 5 50.961 5 51.041 2
      矩震级(Mw) 6.518 6 6.538 1 6.538 3 6.539 8
      地震矩/(1018 N·m) 6.727 5 7.197 6 7.202 4 7.238 9
      RMS/mm 6.3 6.3 6.3
    • 2016-10-26意大利发生Mw 6.1地震,根据美国地震调查局(United States Geological Survey,USGS)网站发布的震源机制,其震中位置为(13.067°E,42.956°N),震源深度为10 km。127个Visso地震GPS台站同震形变数据来自网站http://ring.gm.ingv.it/,814个InSAR升轨位移数据来自文献[20],利用获取的GPS和InSAR数据进行联合反演。首先采用多粒子群算法 [21]进行非线性反演,得到意大利地震断层的走向角为170.3°,倾角为40°,以及断层中心位置等参数,并根据断层参数来构建断层。将断层沿走向/倾向把断层长度、宽度分别拓展至27 km、15 km,并将断层面延伸至地表。将断层面均匀剖分成1 km × 1 km大小的矩形单元,共得到405个矩形单元。分别采用表 1的3种方案确定GPS和InSAR的相对权比以及正则化参数。

      在Visso地震反演实验中,方案1、2、3得到的相对权比分别为1∶2.073 4∶0.002 6、1∶0.000 83∶0.001、1∶2.888 1∶0.001。图 3为利用3种方案反演Visso地震滑动分布结果,其结果与近年来部分学者反演Visso地震滑动分布各参数结果见表 3。其中,文献[22]联合了Sentinel-1和ALOS-2差分干涉数据,及震区GPS三维形变数据,反演得到滑动分布的最大滑动量为0.900 0 m;文献[20]利用Sentinel-1多景影像数据和127个GPS台站形变数据进行联合反演,得到最大滑动量为0.840 0 m。

      图  3  3种方案反演的Visso地震滑动分布结果

      Figure 3.  Results of the Slip Distribution of the Visso Earthquake for Three Schemes

      表 3  Visso地震断层滑动参数

      Table 3.  Slip Parameters for the Visso Earthquake

      方法 相对权比 最大滑动量/m 平均滑动量/m 平均滑动角/(°) 矩震级(Mw) 地震矩/(1018 N·m) RMS/mm
      方案1 1∶2.073 4∶0.002 6 0.605 6 0.194 8 -82.481 5 6.216 2 2.367 3 7.2
      方案2 1∶0.000 83∶0.001 0 0.623 9 0.177 9 -68.742 8 6.189 9 2.161 8 16.2
      方案3 1∶2.888 1∶0.001 0 0.692 6 0.213 7 -85.693 2 6.242 9 2.596 4 6.8
      文献[20] 0.840 0 -67.400 0 6.190 0 2.150 0
      文献[22] 0.900 0 -89.700 0 6.140 0 1.810 0
      文献[23] 0.120 0 -81/-80 6.000 0 1.130 0
      文献[24] 6.160 0 1.970 0
      USGS -89 6.100 0 1.840 0
      GCMT -94 6.100 0 1.610 0

      图 3可知,采用本文3种方案反演Visso地震时,其反演结果均显示Visso地震同震滑动分布区域主要分布在地下0~6 km,其最大滑动区域为4 km,与文献[20]研究结果较一致,方案1~3反演的最大滑动量逐渐增加,但相差不大,且滑动已延伸至地表,这与现场勘查的地表破裂情况吻合。由表 3可以看出,方案1、2、3确定的最大滑动量分别为0.605 6 m、0.623 9 m和0.692 6 m,均低于文献[2022]的结果。这可能是由于两者数据源差造成的,也表明了多类数据能够更好地约束断层,提高反演精度。结合文献[2022]结果发现,方案3较方案1、2能够更好地反演出最大滑动量。方案1、2、3反演的平均滑动角分别为-82.481 5°、-68.742 8°和-85.693 2°,均在各学者的反演结果范围,且方案2反演结果与文献[20]较为接近。方案1、2、3反演矩震级分别为Mw 6.216 2、Mw 6.189 9和Mw 6.242 9,其中方案2结果与文献[20]较为一致,而方案1和3均高于其他文献结果。方案1、2、3反演的RMS分别为7.2 mm、16.2 mm、6.8 mm,方案2结果较差,方案3要略优于方案1。综合以上分析,方案3在方案1的基础上进行偏差改正,两者反演结果较为一致,且方案3反演精度更高。而方案2采用了MSE法求解正则化参数,仅利用方差分量估计求解相对权比,反演结果与方案1和3存在差异,反演参数与文献[20]较为一致。

    • 2016-10-30意大利继Visso地震之后又发生了Mw 6.6地震,根据USGS网站发布的震源机制,其震中位置为(13.096°E,42.862°N),震源深度为8 km。本文采用的Norcia地震114个GPS台站同震形变数据来自网站http://ring.gm.ingv.it/,采用的1 820个InSAR升轨数据来自文献[20],并利用获得的GPS和InSAR进行同震滑动联合反演。同样采用§3.1中的非线性反演算法获得Norcia地震断层的走向角、倾角及断层中心位置参数。构建好断层之后,将断层沿走向/倾向把断层长度、宽度分别拓展至27 km、15 km,并将破裂面延伸至地表;将断层面均匀剖分成1 km×1 km大小的矩形单元,共得到405个矩形单元。分别采用表 3的方案确定GPS和InSAR的相对权比以及正则化参数。

      在Norcia地震反演实验中,方案1、2、3确定的相对权比分别为1∶5.268 4∶0.001 1、1∶0.005 8∶0.01和1∶4.196 4∶0.000 39。图 4为利用3种方案反演Norcia地震滑动分布结果。

      图  4  3种方案反演的Norcia地震的滑动分布结果

      Figure 4.  Results of the Slip Distribution of the Norcia Earthquake for Three Schemes

      图 4可知,利用3种方案反演Norcia地震时,其反演结果均显示Norcia地震同震滑动分布区域主要分布在地下0~6 km,且地震破裂到地表,与文献[20]结果一致;3种方案反演结果较为一致,方案2反演的最大滑动量最大,其次为方案3和方案1。

      本文统计了3种方案反演Norcia地震滑动分布各参数以及近年来部分学者反演Norcia地震滑动分布各参数,结果见表 4。由表 4可以看出,方案1、2、3反演的最大滑动量分别为2.013 6 m、2.484 4 m和2.267 3 m,均小于文献[20]结果,其中方案2与文献[22]结果接近,可能是因为融合了多类数据,能够更好地约束断层。方案1、2、3反演的平均滑动角分别为-87.342 6°、-84.770 7°和-84.894 7°,其中方案2、3反演的结果与USGS较为接近。方案1、2、3反演的矩震级分别为Mw 6.584 4、Mw 6.599 6、Mw 6.589 3,反演的地震矩分别为8.446 2×1018 N·m、8.900 8×1018 N·m、8.590 2×1018 N·m,略低于USGS公布的结果,方案1、3与文献[24]反演结果较为一致,方案2与文献[2022]较为接近。方案1、2、3反演的RMS分别为13.4 mm、22.3 mm、12.7 mm,从RMS角度分析,方案2结果较差,方案3略优于方案1。

      表 4  Norcia地震断层滑动参数

      Table 4.  Slip Parameters for the Norcia Earthquake

      方法 相对权比 最大滑动量/m 平均滑动量/m 平均滑动角/(°) 矩震级(Mw) 地震矩/(1018 N·m) RMS/mm
      方案1 1∶5.268 4∶0.001 1 2.013 6 0.695 2 -87.342 6 6.584 4 8.446 2 13.4
      方案2 1∶0.005 8∶0.01 2.484 4 0.732 6 -84.770 7 6.599 6 8.900 8 22.3
      方案3 1∶4.196 4∶0.000 39 2.267 3 0.707 0 -84.894 7 6.589 3 8.590 2 12.7
      文献[20] 3.440 0 -64.100 0 6.600 0 9.140 0
      文献[22] 2.500 0 6.600 0 8.970 0
      文献[23] 0.420 0 -90.000 0 6.400 0 4.430 0
      文献[24] 6.590 0 8.460 0
      USGS -84.000 0 6.600 0 10.700 0
      GCMT -96.000 0 6.600 0 9.580 0
    • 本文采用了3种方案对模拟实验、Visso地震和Norcia地震进行滑动分布反演,并对3种方案进行对比分析,综合实验结果发现,本文3种方案均能很好地反演出合理的相对权比与正则化参数。但无论是对参数进行偏差改正还是对残差进行偏差改正,都能够有效地提高最大滑动量,因为两者均考虑到反演中偏差的影响,但两种方案针对不同的实际情况,改进的程度有所不同。方案1、3反演的滑动分布较为一致,而方案2反演结果在主滑动区较为一致,在周围区域存在差别,这可能与方案2采用MSE法求解正则化参数相关。对于3种方案反演的矩震级,在模拟实验中,3种方案都能较准确地反演出实际的矩震级,而在本文所采用的两个实际地震中,方案2的反演结果更为准确,Visso地震反演实验中,方案1、3的结果要略高于其他文献结果。对于RMS值,方案3的结果明显优于其他两种方案,在两个实际地震中,方案2的RMS值明显高于方案1、3,由此也可以说明方案1、3反演的滑动分布更为合理。

    • 本文分别采用了方差分量估计、对残差进行偏差改正的方差分量估计和对参数进行偏差改正方差分量估计法3种方法确定联合反演的相对权比和正则化参数,验证了3种方法的合理性与可行性,且进行偏差改正后的结果优于未进行偏差改正的方差分量估计法。针对不同偏差改正的方式在不同地震中的应用,其反演结果存在差异。针对残差进行偏差改正的初衷在于参数的偏差可能大于参数本身,但通过本文的实验表明,对参数进行偏差改正时,其反演结果改正较为显著,这也表明了偏差改正方差分量估计法在滑动分布反演中具有一定优势。

参考文献 (24)

目录

    /

    返回文章
    返回