留言板

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

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

一种基于F-J线性-非线性模型解的迭代最小二乘方法

张彦栋 许才军 汪建军

张彦栋, 许才军, 汪建军. 一种基于F-J线性-非线性模型解的迭代最小二乘方法[J]. 武汉大学学报 ● 信息科学版, 2019, 44(12): 1816-1822. doi: 10.13203/j.whugis20180117
引用本文: 张彦栋, 许才军, 汪建军. 一种基于F-J线性-非线性模型解的迭代最小二乘方法[J]. 武汉大学学报 ● 信息科学版, 2019, 44(12): 1816-1822. doi: 10.13203/j.whugis20180117
ZHANG Yandong, XU Caijun, WANG Jianjun. An Iterative Least Squares Method Based on the F-J Inversion for Linear-Non-Linear Models[J]. Geomatics and Information Science of Wuhan University, 2019, 44(12): 1816-1822. doi: 10.13203/j.whugis20180117
Citation: ZHANG Yandong, XU Caijun, WANG Jianjun. An Iterative Least Squares Method Based on the F-J Inversion for Linear-Non-Linear Models[J]. Geomatics and Information Science of Wuhan University, 2019, 44(12): 1816-1822. doi: 10.13203/j.whugis20180117

一种基于F-J线性-非线性模型解的迭代最小二乘方法

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

国家重点研发计划 2018YFC1503605

国家自然科学基金 41431069

国家自然科学基金 41574002

国家自然科学基金 41774011

详细信息
    作者简介:

    张彦栋, 硕士, 主要研究方向为地球物理大地测量反演算法。ydzhang@whu.edu.cn

    通讯作者: 汪建军, 博士, 副教授, 主要从事地震应力触发研究。jjwang@sgg.whu.edu.cn
  • 中图分类号: P228;P207

An Iterative Least Squares Method Based on the F-J Inversion for Linear-Non-Linear Models

Funds: 

The National Key Research and Development Program of China 2018YFC1503605

the National Natural Science Foundation of China 41431069

the National Natural Science Foundation of China 41574002

the National Natural Science Foundation of China 41774011

More Information
    Author Bio:

    ZHANG Yandong, master, majors in inversion algorithms for geophysical and geodetic data. E-mail:ydzhang@whu.edu.cn

    Corresponding author: WANG Jianjun, PhD, associate professor, specializes in earthquake stress triggering.E-mail:jjwang@sgg.whu.edu.cn
  • 摘要: 基于贝叶斯理论的线性与非线性模型反演方法(Fukuda-Johnson,F-J)已广泛应用于地球物理模型的线性-非线性参数反演。但F-J方法的反演结果可能受马尔可夫链蒙特卡洛采样(Markov chain Monte Carlo,MCMC)经验参数选择的影响,而反复调试合适的经验参数需耗费大量计算时间。对线性与非线性模型进行线性化后,也可以利用迭代最小二乘方法反演,但该方法难以选择合适的初始值。为提高参数反演计算效率和避免参数初值选择影响,提出了一种以F-J方法模型解为初始值的迭代最小二乘方法。该方法只需计算一次F-J方法模型解和有限次最小二乘迭代,既提高了F-J方法的反演效率,又能获得迭代最小二乘全局最优解。针对模拟数据实验和实际数据算例,分别采用F-J方法、随机生成初始值的迭代最小二乘方法和以F-J方法结果为初值的迭代最小二乘方法进行参数反演。结果表明,直接使用F-J方法时,MCMC采样参数会影响反演结果;直接进行迭代最小二乘反演时,初始值选取不当会导致迭代无法收敛到正确的结果;以F-J方法的结果作为迭代最小二乘方法的初始值进行反演,可以充分发挥F-J方法的全局最优性和迭代最小二乘方法计算量小、稳定性好的优势。
  • 图  1  长白山天池火山水准点和GPS点分布图

    Figure  1.  Distribution of the Leveling Benchmarks and GPS Stations Surrounding Changbaishan Tianchi Volcano

    表  1  不同情况下F-J方法求出的参数估计值

    Table  1.   Estimation and Variance of Parameters by Using the F-J Method at Different Conditions

    情况 D/km x0/km y0/km σ12/mm2 σ22/mm2 σ32/(10-12m2∙s-4) ΔV/103 m3
    参考结果 4.007±0.016 0.005±0.013 0.006±0.011 3.906±0.451 4.133±0.187 9.138±0.777 6.007±0.032
    4.007±0.015 0.005±0.009 0.010±0.008 3.910±0.358 4.092±0.144 9.198±0.604 6.096±0.029
    4.353±0.105 0.007±0.029 0.011±0.025 3.797±0.554 4.351±0.290 9.384±0.911 6.302±0.150
    4.026±0.027 0.019±0.018 0.001±0.017 5.292±0.504 4.014±0.207 9.194±0.873 6.344±0.037
    4.029±0.010 0.028±0.003 0.022±0.010 5.380±0.360 3.957±0.141 8.299±0.686 6.548±0.026
    马尔可夫链发散
    下载: 导出CSV

    表  2  迭代最小二乘方法求出的参数估计值

    Table  2.   Results of the Iterative Least Squares Method with Initial Values Resulted from the F-J Method

    统计项目 D/km x0/km y0/km σ12/mm2 σ22/mm2 σ32/(10-12m2∙s-4) ΔV/103 m3
    真值 4 0 0 4 4 9 6
    迭代收敛结果 4.013 0±0.018 2 0.005 7±0.009 4 0.007 6±0.013 6 3.905 7 4.132 4 9.137 7 6.016 2±0.035 0
    下载: 导出CSV

    表  3  实测数据不同条件下F-J方法求出的参数估计值

    Table  3.   Estimation and Variance of Parameters by Using the F-J Method of Experimental Data at Different Conditions

    情况 D/km x0/km y0/km σ12/mm2 σ22/mm2 ΔV/103 m3
    参考结果 7.533 4±0.629 3 2.030 1±0.978 1 0.836 1±0.659 0 6.925 9±2.671 2 6.302 6±2.391 3 22.250 2±3.278 4
    7.554 7±0.548 5 2.029 6±0.880 0 0.839 6±0.514 2 6.921 6±2.421 8 6.299 7±1.915 7 22.249 5±2.792 2
    7.896 0±0.748 2 1.655 7±1.195 6 1.035 7±0.765 3 8.491 7±3.045 4 6.787 5±2.915 2 23.925 5±3.940 4
    7.034 1±0.764 7 2.170 5±1.261 6 0.559 1±0.697 5 5.703 7±2.285 3 4.535 8±2.651 5 21.965 6±3.944 5
    7.508 2±0.531 9 2.559 2±0.779 5 0.748 3±0.472 8 6.980 5±2.566 9 4.633 3±1.916 7 22.574 3±3.056 3
    马尔可夫链发散
    下载: 导出CSV

    表  4  各研究方法结果的参数比较

    Table  4.   Parameters Comparison of Different Research Methods

    方法 D/km x0/km y0/km ΔV/103 m3
    本文F-J方法 7.533 4±0.629 3 2.030 1±0.978 1 0.836 1±0.659 0 22.250 2±3.278 4
    基于F-J方法解的迭代最小二乘方法 7.854 8±0.548 5 1.818 4±1.025 2 1.075 5±1.025 0 23.294 6±3.592 7
    文献[17]方法 6.9 2.5 0.5 5.2
    文献[11]最小二乘方法 10.335 4±1.084 9 3.005 9±0.930 5 3.936 9±1.051 2 22.463 2±1.870 2
    文献[11]TLS方法 5.606 9±0.020 9 2.257 4±0.021 0 1.356 0±0.020 9 12.814 0±0.021 0
    下载: 导出CSV

    表  5  F-J方法求出的参数估计值及误差

    Table  5.   Estimation and Variance of Parameters by Using the F-J Method

    情况 D1/km D2/km σ12/(mm·a-1)2 σ22/(mm·a-1)2 σ32/(mm·a-1)2 1/(mm·a-1) 2/(mm·a-1)
    真值 10 10 1.00 6.25 20.25 25 15
    参考结果 10.004±1.029 9.940±1.579 1.023±0.074 6.450±0.263 20.341±0.590 24.712±1.218 15.151±1.217
    10.004±0.896 9.940±1.316 1.024±0.060 6.450±0.256 20.341±0.496 24.712±1.004 15.151±0.985
    10.591±1.144 8.855±1.816 1.027±0.082 6.449±0.254 20.375±0.595 25.348±1.331 14.552±1.352
    10.829±1.463 8.582±1.826 1.012±0.113 6.737±0.391 20.529±0.963 26.130±1.810 13.779±1.828
    9.561±0.705 10.399±1.036 1.032±0.057 6.457±0.174 20.532±0.384 24.058±0.724 15.774±0.758
    马尔可夫链发散
    下载: 导出CSV

    表  6  以F-J方法的解为初值的迭代最小二乘结果

    Table  6.   Results of the Iterative Least Squares Method with Initial Values Resulted from the F-J Method

    情况 D1/km D2/km σ12/(mm·a-1)2 σ22/(mm·a-1)2 σ32/(mm·a-1)2 1/(mm·a-1) 2/(mm·a-1)
    真值 10 10 1.0 6.25 20.25 25 15
    迭代收敛结果 10.002±1.031 9.948±1.616 1.024±0.076 6.455±0.262 20.349±0.588 24.717±1.228 15.154±1.220
    下载: 导出CSV
  • [1] 许才军, 大地测量联合反演理论和方法研究进展[J].武汉大学学报·信息科学版, 2001, 26(6):555-561 http://ch.whu.edu.cn/CN/abstract/abstract5232.shtml

    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/CN/abstract/abstract5232.shtml
    [2] Fukuda J, Johnson K M. A Fully Bayesian Inversion for Spatial Distribution of Fault Slip with Objective Smoothing[J]. Bulletin of the Seismological Society of America, 2008, 98(3):1128-1146 doi:  10.1785/0120070194
    [3] Fukuda J, Johnson K M. Mixed Linear-Non-Linear Inversion Crustal Deformation Data: Bayesian Inference of Model, Weighting and Regularization Parameters[J]. Geophysical Journal International, 2010, 181(3):1441-1458 http://cn.bing.com/academic/profile?id=6fa2ca1e5794f8891e4f605ac6a940e4&encoded=0&v=paper_preview&mkt=zh-cn
    [4] Du Y, Segall P, Aydin A. Comparison of Various Inversion Techniques as Applied to the Determination of a Geophysical Deformation Model for the 1983 Borah Peak Earthquake[J]. Bulletin of the Seismological Society of America, 1992, 82(4):1840-1866 http://cn.bing.com/academic/profile?id=1a8d759197dcb54712aa44832d9a97af&encoded=0&v=paper_preview&mkt=zh-cn
    [5] Freymueller J, King N E, Segall P. The Co-seismic Slip Distribution of the Landers Earthquake[J]. Bulletin of the Seismological Society of America, 1994, 84(3):646-659 http://cn.bing.com/academic/profile?id=af4c80ad80b916a562356856dff04366&encoded=0&v=paper_preview&mkt=zh-cn
    [6] Jonsson S, Zebker H, Segall P, et al. Fault Slip Distribution of the 1999 Mw 7.1 Hector Mine, California, Earthquake, Estimated from Satellite Radar and GPS Measurements[J]. Bulletin of the Seismological Society of America, 2002, 92(4):1377-1389 doi:  10.1785/0120000922
    [7] Wright T J, Lu Z, Wicks C. Constraining the Slip Distribution and Fault Geometry of the Mw 7.9, 3 November 2002, Denali Fault Earthquake with Interferometric Synthetic Aperture Radar and Global Positioning System Data[J]. Bulletin of the Seismological Society of America, 2004, 94(6B):175-189 doi:  10.1785/0120040623
    [8] Simons M, Fialko Y, Rivera L. Coseismic Deformation from the 1999 Mw 7.1 Hector Mine, California, Earthquake as Inferred from InSAR and GPS Observations[J]. Bulletin of the Seismological Society of America, 2002, 92(4):1390-1402 doi:  10.1785/0120000933
    [9] Xu Caijun, Ding Kaihua, Cai Jianqing, 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
    [10] Metropolis N, Rosenbluth A W, Rosenbluth M N, et al. Equation of State Calculations by Fast Computing Machines[J]. The Journal of Chemical Physics, 1953, 21(6):1087-1092 doi:  10.1063/1.1699114
    [11] 王乐洋.基于总体最小二乘的大地测量反演理论及应用研究[D].武汉: 武汉大学, 2011 http://www.cnki.com.cn/Article/CJFDTotal-CHXB201204029.htm

    Wang Leyang. Research on Theory and Application of Total Least Squares in Geodetic Inversion[D].Wuhan: Wuhan University, 2011 http://www.cnki.com.cn/Article/CJFDTotal-CHXB201204029.htm
    [12] Yamakawa N. On the Strain Produced in a Semi-infinite Elastic Solid by an Interior Source of Stress[J]. Zisin, 1955, 2(8):84-98 http://www.wanfangdata.com.cn/details/detail.do?_type=perio&id=J-STAGE_39860
    [13] Scarpa R, Tilling R I. Monitoring and Mitigation of Volcano Hazards[M]. Berlin: Springer-Verlag, 1996:373-402
    [14] 胡亚轩, 王庆良, 崔笃信, 等.三种压力源模型对火山区地面变形的影响[J].东北地震研究, 2005, 21(3):33-38 doi:  10.3969/j.issn.1674-8565.2005.03.003

    Hu Yaxuan, Wang Qingliang, Cui Duxin, et al. Influences on Surface Deformation by the Three Different Stress Source Models in Volcanic Area[J]. Seismological Research of Northeast China, 2005, 21(3):33-38 doi:  10.3969/j.issn.1674-8565.2005.03.003
    [15] Hagiwara Y. The Mogi Model as a Possible Cause of the Crustal Uplift in the Eastern Part of Izu Peninsula and the Related Gravity Change[J]. Bulletin of the Earthquake Research Institute, 1977, 52:301-309
    [16] 陈国浒.长白山天池火山形变监测与模拟研究[D].北京: 中国地震局地质研究所, 2007 http://cdmd.cnki.com.cn/article/cdmd-85402-2007090274.htm

    Chen Guohu. Deformation Monitoring and Simulation Research of the Tianchi Volcano in the Changbaishan Mountains[D]. Beijing: Institute of Geology, China Earthquake Administration, 2007 http://cdmd.cnki.com.cn/article/cdmd-85402-2007090274.htm
    [17] Savage J C, Burford R O. Geodetic Determination of Relative Plate Motion in Central California[J]. Journal of Geophysical Research, 1973, 78(5):832-845 doi:  10.1029/JB078i005p00832
  • [1] 吕志鹏, 隋立芬.  基于非线性高斯-赫尔默特模型的结构总体最小二乘法 . 武汉大学学报 ● 信息科学版, 2019, 44(12): 1808-1815. doi: 10.13203/j.whugis20180104
    [2] 尹建鹏, 许才军.  利用组合算法反演球体分层模型的断层参数 . 武汉大学学报 ● 信息科学版, 2019, 44(9): 1320-1327. doi: 10.13203/j.whugis20170380
    [3] 王乐洋, 余航.  火山Mogi模型反演的总体最小二乘联合平差方法 . 武汉大学学报 ● 信息科学版, 2018, 43(9): 1333-1341. doi: 10.13203/j.whugis20160469
    [4] 何平 许才军 温扬茂 丁开华 王庆良.  利用PALSAR数据研究长白山火山活动性 . 武汉大学学报 ● 信息科学版, 2015, 40(2): 214-221. doi: 10.13203/j.whugis20130077
    [5] 姚宜斌, 黄书华, 陈家君.  求解自回归模型参数的整体最小二乘新方法 . 武汉大学学报 ● 信息科学版, 2014, 39(12): 1463-1466.
    [6] 宫轶松, 归庆明, 李保利, 边少峰.  动态非线性滤波模型非线性强度的曲率度量及其应用 . 武汉大学学报 ● 信息科学版, 2011, 36(8): 904-908.
    [7] 邱卫宁, 齐公玉, 田丰瑞.  整体最小二乘求解线性模型的改进算法 . 武汉大学学报 ● 信息科学版, 2010, 35(6): 708-710.
    [8] 游为, 范东明, 傅淑娟.  同伦函数与填充函数相结合的非线性最小二乘平差模型 . 武汉大学学报 ● 信息科学版, 2010, 35(2): 185-188.
    [9] 宋迎春, 朱建军, 罗德仁, 曾联斌.  附非负约束平差模型的最小二乘估计 . 武汉大学学报 ● 信息科学版, 2008, 33(9): 907-909.
    [10] 孙红星, 付建红, 袁修孝, 唐卫明.  基于多历元递推最小二乘卡尔曼滤波方法的模糊度解算 . 武汉大学学报 ● 信息科学版, 2008, 33(7): 735-739.
    [11] 宁伟, 陶华学, 卿熙宏.  广义非线性最小二乘测量参数平差的快速差分迭代解算 . 武汉大学学报 ● 信息科学版, 2005, 30(7): 617-620.
    [12] 苗作华, 刘耀林, 王海军.  耕地需求量预测的加权模糊-马尔可夫链模型 . 武汉大学学报 ● 信息科学版, 2005, 30(4): 305-308.
    [13] 张勤, 陶本藻.  基于同伦法的非线性最小二乘平差统一模型 . 武汉大学学报 ● 信息科学版, 2004, 29(8): 708-710.
    [14] 曾文宪, 陶本藻.  三维坐标转换的非线性模型 . 武汉大学学报 ● 信息科学版, 2003, 28(5): 566-568.
    [15] 陶本藻, 张勤.  GPS非线性数据处理的同伦最小二乘模型 . 武汉大学学报 ● 信息科学版, 2003, 28(S1): 115-118.
    [16] 陶本藻.  形变反演模型的非线性平差 . 武汉大学学报 ● 信息科学版, 2001, 26(6): 504-508.
    [17] 王新洲.  非线性模型参数估计的直接解法 . 武汉大学学报 ● 信息科学版, 1999, 24(1): 64-67.
    [18] 王新洲.  非线性模型能否线性化的实用判据 . 武汉大学学报 ● 信息科学版, 1999, 24(2): 145-148.
    [19] 王新洲.  非线性模型线性近似的容许曲率 . 武汉大学学报 ● 信息科学版, 1997, 22(2): 119-121.
    [20] 白亿同.  非线性最小二乘平差迭代解法的收敛性 . 武汉大学学报 ● 信息科学版, 1991, 16(2): 92-95.
  • 加载中
图(1) / 表(6)
计量
  • 文章访问数:  857
  • HTML全文浏览量:  93
  • PDF下载量:  168
  • 被引次数: 0
出版历程
  • 收稿日期:  2018-12-25
  • 刊出日期:  2019-12-05

一种基于F-J线性-非线性模型解的迭代最小二乘方法

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

    国家重点研发计划 2018YFC1503605

    国家自然科学基金 41431069

    国家自然科学基金 41574002

    国家自然科学基金 41774011

    作者简介:

    张彦栋, 硕士, 主要研究方向为地球物理大地测量反演算法。ydzhang@whu.edu.cn

    通讯作者: 汪建军, 博士, 副教授, 主要从事地震应力触发研究。jjwang@sgg.whu.edu.cn
  • 中图分类号: P228;P207

摘要: 基于贝叶斯理论的线性与非线性模型反演方法(Fukuda-Johnson,F-J)已广泛应用于地球物理模型的线性-非线性参数反演。但F-J方法的反演结果可能受马尔可夫链蒙特卡洛采样(Markov chain Monte Carlo,MCMC)经验参数选择的影响,而反复调试合适的经验参数需耗费大量计算时间。对线性与非线性模型进行线性化后,也可以利用迭代最小二乘方法反演,但该方法难以选择合适的初始值。为提高参数反演计算效率和避免参数初值选择影响,提出了一种以F-J方法模型解为初始值的迭代最小二乘方法。该方法只需计算一次F-J方法模型解和有限次最小二乘迭代,既提高了F-J方法的反演效率,又能获得迭代最小二乘全局最优解。针对模拟数据实验和实际数据算例,分别采用F-J方法、随机生成初始值的迭代最小二乘方法和以F-J方法结果为初值的迭代最小二乘方法进行参数反演。结果表明,直接使用F-J方法时,MCMC采样参数会影响反演结果;直接进行迭代最小二乘反演时,初始值选取不当会导致迭代无法收敛到正确的结果;以F-J方法的结果作为迭代最小二乘方法的初始值进行反演,可以充分发挥F-J方法的全局最优性和迭代最小二乘方法计算量小、稳定性好的优势。

English Abstract

张彦栋, 许才军, 汪建军. 一种基于F-J线性-非线性模型解的迭代最小二乘方法[J]. 武汉大学学报 ● 信息科学版, 2019, 44(12): 1816-1822. doi: 10.13203/j.whugis20180117
引用本文: 张彦栋, 许才军, 汪建军. 一种基于F-J线性-非线性模型解的迭代最小二乘方法[J]. 武汉大学学报 ● 信息科学版, 2019, 44(12): 1816-1822. doi: 10.13203/j.whugis20180117
ZHANG Yandong, XU Caijun, WANG Jianjun. An Iterative Least Squares Method Based on the F-J Inversion for Linear-Non-Linear Models[J]. Geomatics and Information Science of Wuhan University, 2019, 44(12): 1816-1822. doi: 10.13203/j.whugis20180117
Citation: ZHANG Yandong, XU Caijun, WANG Jianjun. An Iterative Least Squares Method Based on the F-J Inversion for Linear-Non-Linear Models[J]. Geomatics and Information Science of Wuhan University, 2019, 44(12): 1816-1822. doi: 10.13203/j.whugis20180117
  • 随着观测技术的不断发展,各种大地测量观测数据和地球物理数据为地球科学研究提供了丰富的源数据,多类数据联合反演成为大地测量和地球物理反演的趋势[1]。联合反演模型中,部分参数与观测值呈线性关系(线性参数),另一部分呈非线性关系(非线性参数)。线性-非线性模型联合反演方法通常要估计模型中的线性参数和非线性参数,也要估计多类观测数据的权重参数和约束条件中的正则化参数(二者统称为超参数)[2-3]

    传统反演方法是先固定权重参数和正则化参数,然后对非线性参数做线性化处理,最后采用最小二乘方法反演[4-6]。该方法的正则化参数和权重参数分别采用L曲线法[4, 6-7]和试错法[8]确定。这种反演方法的不足为:①需要选择非线性参数的近似值,如果选取不合适会影响反演结果;②权重参数和正则化参数的选择具有主观性。因此,实际计算中通常采用方差分量估计,经迭代计算求出权重参数和正则化参数的最佳估计值[9]。但是,当其初值选取不合适时,会导致反演解迭代不收敛或者收敛到局部最小。

    文献[2]提出了一种基于贝叶斯理论的反演方法,称为全贝叶斯方法。该方法首先推导出线性参数、非线性参数、权重参数以及正则化参数的概率密度分布,然后通过马尔可夫链蒙特卡洛(Markov chain Monte Carlo,MCMC)采样对概率密度函数采样,计算参数的估计值及其精度。全贝叶斯方法求出的概率密度函数非常复杂,无法直接求解参数的期望和方差;且该方法采用MCMC采样的同时估计线性参数和非线性参数,极大地增加了MCMC采样计算量。基于此,Fukuda和Johnson[3]进一步改进了全贝叶斯方法,采用最小二乘估计线性参数,降低了MCMC采样计算量,该方法简称F-J方法。但该方法结果受MCMC采样经验参数选择影响,且MCMC采样计算量仍较大。基于此,本文提出一种基于F-J线性-非线性模型解的迭代最小二乘方法。该方法结合了F-J方法解的全局最优性和迭代最小二乘方法的计算高效性。

    • 贝叶斯理论将待求参数视作随机变量且采用概率分布进行描述,贝叶斯公式的概率密度函数如下:

      $$ p(\mathit{\boldsymbol{z}}|\mathit{\boldsymbol{d}}) = \frac{{p(\mathit{\boldsymbol{d}}|\mathit{\boldsymbol{z}})p\left( \mathit{\boldsymbol{z}} \right)}}{{p\left( \mathit{\boldsymbol{d}} \right)}} $$ (1)

      式中,z是模型中待求的线性参数(如断层的滑动位错量参数)、非线性参数(如断层的走向和倾角参数)或超参数(如断层滑动分布反演中的平滑参数);d是观测数据;p(z|d)是在观测数据d已知条件下参数z的后验概率密度;p(d|z)是参数z已知条件下观测数据的概率密度;p(z)代表参数的先验概率密度;p(d)是归一化常数。在线性-非线性混合反演问题中,变量z包含非线性参数m、观测数据的权重参数$\sigma _1^2, \sigma _2^2 \cdots \sigma _n^2$和约束条件中的正则化参数$\alpha _1^2, \alpha _2^2 \cdots \alpha _J^2$。

      在非线性参数已知的情况下,线性参数服从正态分布,而z不服从正态分布。可以使用MCMC方法对z的概率密度函数进行采样,并根据每一个采样点计算对应的线性参数,估计所有参数的值和精度。

    • F-J方法中,MCMC采样由Metropolis算法[10]实现。在该算法中,记非线性参数x=$\mathit{\boldsymbol{x}} = \left( {\mathit{\boldsymbol{m}}, \sigma _1^2, \sigma _2^2 \cdots \sigma _n^2, \alpha _1^2, \alpha _2^2 \cdots \alpha _J^2} \right)$的验后分布为p(x|d)。采用收敛的Markov随机游走过程实现该验后概率分布采样。

      x(i)为第i次循环时x的取值,利用Markov链随机游走过程对当前状态加上扰动得到x(i+1)的候选值x',其被接受的概率为:

      $${P_{{\rm{accept}}}} = {\rm{min}}\left( {1, \frac{{p(\mathit{\boldsymbol{x}}'|\mathit{\boldsymbol{d}})}}{{p({\mathit{\boldsymbol{x}}^{\left( i \right)}}|\mathit{\boldsymbol{d}})}}} \right)$$ (2)

      Metropolis算法中,采样点数量、随机游走步长、舍弃的最初采样点数量和选取采样点的间隔等参数选取都会影响到反演结果。因此,基于F-J方法的反演结果与采样参数的选择有关。

    • 将对函数模型中的非线性参数m线性化,进行最小二乘计算和精度估计。观测值的线性-非线性模型为:

      $${\mathit{\boldsymbol{d}}_k} = {\mathit{\boldsymbol{G}}_k}\left( \mathit{\boldsymbol{m}} \right)\mathit{\boldsymbol{s}} + {\varepsilon _k}, k = 1, 2 \cdots K$$ (3)

      式中,dk是第k类观测数据;mMN维非线性参数;sML维线性参数;Gk(m)是Nk×ML维矩阵;假定εk服从正态分布,其均值为0,方差为σk2Σk

      线性参数的先验约束为:

      $${\mathit{\boldsymbol{D}}_j}\mathit{\boldsymbol{s}} = {\mathit{\boldsymbol{\delta }}_j}, j = 1, 2 \cdots J$$ (4)

      式中,DjML×ML的矩阵;δjML维的向量,δj服从正态分布,其均值0,协方差矩阵为αj2IML,其中IMLML×ML的单位矩阵;αj2是比例因子。

      最小二乘方法的目标函数为:

      $$\begin{array}{l} \mathit{\Phi }\left( {\mathit{\boldsymbol{m}}, \mathit{\boldsymbol{s}}} \right) = {\sum\limits_{k = 1} }^K \frac{1}{{\sigma _k^2}}{[{\mathit{\boldsymbol{d}}_k} - {\mathit{\boldsymbol{G}}_k}\left( \mathit{\boldsymbol{m}} \right)\mathit{\boldsymbol{s}}]^{\rm{T}}}\mathit{\boldsymbol{ \boldsymbol{\varSigma} }}_k^{ - 1}\\ \;\;\;\;\;\;\left[ {{\mathit{\boldsymbol{d}}_k} - {\mathit{\boldsymbol{G}}_k}\left( \mathit{\boldsymbol{m}} \right)\mathit{\boldsymbol{s}}} \right] + {\sum\limits_{j = 1} }^J \frac{1}{{\alpha _j^2}}{\left\| {{\mathit{\boldsymbol{D}}_j}\mathit{\boldsymbol{s}}} \right\|^2} = {\rm{min}} \end{array}$$ (5)

      G(m)中,每个元素都是m的函数,即:

      $$\mathit{\boldsymbol{G}}\left( \mathit{\boldsymbol{m}} \right) = \left( {\begin{array}{*{20}{c}} {{g_{11}}\left( \mathit{\boldsymbol{m}} \right)}&{{g_{12}}\left( \mathit{\boldsymbol{m}} \right)}& \cdots &{{g_{1{M_L}}}\left( \mathit{\boldsymbol{m}} \right)}\\ {{g_{21}}\left( \mathit{\boldsymbol{m}} \right)}&{{g_{22}}\left( \mathit{\boldsymbol{m}} \right)}& \cdots &{{g_{2{M_L}}}\left( \mathit{\boldsymbol{m}} \right)}\\ \vdots & \vdots &{}& \vdots \\ {{g_{N1}}\left( \mathit{\boldsymbol{m}} \right)}&{{g_{N2}}\left( \mathit{\boldsymbol{m}} \right)}& \cdots &{{g_{N{M_L}}}\left( \mathit{\boldsymbol{m}} \right)} \end{array}} \right)$$

      观测数据d中的第i个元素对m中的第j个元素的偏导数为$\frac{{\partial {d_i}}}{{\partial {d_j}}} = {\sum\limits_{k = 1} }^{{M_L}} \frac{{\partial {g_{ik}}\left( \mathit{\boldsymbol{m}} \right)}}{{\partial {d_j}}}{s_k}$,从而得到线性化后的函数模型。

      由于线性-非线性模型中的非线性参数存在初始值误差,只有近似值足够准确时才能忽略线性化的误差,故而需要进行迭代最小二乘计算。在每次迭代后,将最小二乘求得的改正数加到前一次最小二乘计算的估计值上,将方差分量的估计值分别乘到前一次得到的方差因子上。经过若干次迭代计算后,如果改正数足够小且各个方差分量的估计值足够接近1,则认为迭代收敛,从而得到最终的反演结果。

    • 为避免F-J方法和迭代最小二乘方法的不足,本文采用基于F-J方法模型解的迭代最小二乘方法。其具体思想是:首先利用F-J方法和一组采样参数同时反演模型的非线性参数和线性参数;然后基于反演参数解,利用迭代最小二乘方法反演模型的最优解。由于F-J方法提供了较好的初值,迭代最小二乘方法能较快地迭代收敛到全局最优解。该组合方法既避免了F-J方法中概率密度采样费时低效和反演结果依赖采样参数的问题,也避免了单独基于最小二乘方法的迭代反演解可能陷入局部最优的问题。

    • Mogi模型在火山区地表形变和重力变化研究中已得到广泛的应用[11]。火山岩浆活动引起的地表垂直位移的表达式为[12]

      $${\rm{\Delta }}h = \frac{{3D{\rm{\Delta }}V}}{{4{\rm{ \mathsf{ π} }}{{({D^2} + {r^2})}^{3/2}}}}$$ (6)

      式中,Δh为膨胀源引起的垂直位移;D为膨胀源的深度;r为监测点到膨胀源之间的地表径向位移;ΔV为体积增量。

      火山岩浆活动引起的地表水平位移[13-14]为:

      $$\left\{ {\begin{array}{*{20}{l}} {{U_x} = \frac{{3{\rm{\Delta }}V\left( {x - {x_0}} \right)}}{{4{\rm{ \mathsf{ π} }}{{({D^2} + {r^2})}^{3/2}}}}}\\ {{U_y} = \frac{{3{\rm{\Delta }}V\left( {y - {y_0}} \right)}}{{4{\rm{ \mathsf{ π} }}{{({D^2} + {r^2})}^{3/2}}}}} \end{array}} \right.$$ (7)

      式中,UxUy分别为x方向和y方向的位移;r为监测点到膨胀源之间的地表径向位移;D为膨胀源的深度;ΔV为体积增量;(x0, y0)为岩浆囊中心在地表的投影坐标;其余变量同式(6)。

      火山岩浆活动引起的地表重力变化为[15]

      $${\rm{\Delta }}g = \frac{{G\rho D{\rm{\Delta }}V}}{{{{({D^2} + {r^2})}^{3/2}}}}$$ (8)

      式中,Δg为重力变化;G为万有引力常量;ρ为膨胀体内异常体的密度;其余变量同式(6)。

      根据监测点的垂直位移、水平位移和重力变化,可以联合反演D、ΔV和(x0, y0)。Mogi模型含有线性-非线性参数,其中ΔV为线性参数,D和(x0, y0)为非线性参数。

    • 模拟观测区域为东西向和南北向均长20 km的正方形,其中心点坐标为(0, 0),相邻垂直形变观测点间隔为1 km,由此得到在区域内均匀分布的21×21=441个观测点。设置Mogi模型的参数值如下:D=4 km,(x0, y0)为(0,0)km,ΔV= 6×103 m3。根据Mogi模型分别正演计算由膨胀源引起的地表垂直位移、地表水平位移和地表重力变化(在计算重力变化时,取ρ=3.38×103 kg/m3),并对3种数据分别加上均值为0、标准差分别为σ1=2 mm、σ2=2 mm、σ3=3×10-6 m/s2的误差,得到3种不同模拟观测值。

      首先用F-J方法进行反演。MCMC采样参数参考值选取如下:采样点的间隔设为每10个取一个;采样点次数为100万次;舍弃最初采样点数为10 000个。仅改变一个MCMC采样参数参考值的方法包括以下几种情况:①减小选取采样点的间隔,即由参考结果中原始采样点每10个取一个改为每5个取一个;②减少采样点数量,即由100万次减少到1万次;③舍弃更少的最初采样点,即由参考结果中舍弃前面10 000个采样点改为只舍弃前面1 000个采样点;④增大随机游走步长为3倍;⑤增大随机游走步长为6倍。基于这些MCMC采样参数,反演的Mogi模型线性和非线性参数如表 1所示。可以看出,采样参数的选取不仅影响后验估计参数的值,也影响其精度。其中,改变MCMC采样经验参数中的采样点数、舍弃初始采样点数和随机游走步长对反演参数结果影响比较明显。

      表 1  不同情况下F-J方法求出的参数估计值

      Table 1.  Estimation and Variance of Parameters by Using the F-J Method at Different Conditions

      情况 D/km x0/km y0/km σ12/mm2 σ22/mm2 σ32/(10-12m2∙s-4) ΔV/103 m3
      参考结果 4.007±0.016 0.005±0.013 0.006±0.011 3.906±0.451 4.133±0.187 9.138±0.777 6.007±0.032
      4.007±0.015 0.005±0.009 0.010±0.008 3.910±0.358 4.092±0.144 9.198±0.604 6.096±0.029
      4.353±0.105 0.007±0.029 0.011±0.025 3.797±0.554 4.351±0.290 9.384±0.911 6.302±0.150
      4.026±0.027 0.019±0.018 0.001±0.017 5.292±0.504 4.014±0.207 9.194±0.873 6.344±0.037
      4.029±0.010 0.028±0.003 0.022±0.010 5.380±0.360 3.957±0.141 8.299±0.686 6.548±0.026
      马尔可夫链发散

      本文同样采用直接进行迭代最小二乘计算的方法,初始值在一定的范围内随机选取,参数Dx0y0σ12σ22σ32和ΔV的选取范围依次为[0,10]、[-5,5]、[-5,5]、[1,16]、[1,16]、[1,20]和[1,15]。实验中共随机生成了30组迭代初始值,经实验发现,有16组初始值不能收敛,5组收敛到错误结果,9组收敛到正确结果。可见迭代最小二乘的结果与初值选择有非常大的关系。再分别采用以上5组F-J方法计算得到的结果为初始值进行最小二乘计算,每组都收敛到相同的结果,如表 2所示。

      表 2  迭代最小二乘方法求出的参数估计值

      Table 2.  Results of the Iterative Least Squares Method with Initial Values Resulted from the F-J Method

      统计项目 D/km x0/km y0/km σ12/mm2 σ22/mm2 σ32/(10-12m2∙s-4) ΔV/103 m3
      真值 4 0 0 4 4 9 6
      迭代收敛结果 4.013 0±0.018 2 0.005 7±0.009 4 0.007 6±0.013 6 3.905 7 4.132 4 9.137 7 6.016 2±0.035 0

      以上模拟数据实验表明,在MCMC采样参数不导致随机游走过程发散的前提下,基于F-J方法反演解的迭代最小二乘方法可以得到具有相同收敛结果的线性-非线性参数解。

    • 长白山天池火山区形变监测的大地测量观测主要有布设于长白山北坡、长约25 km的13个水准点组成的水准监测剖面和8个全球定位系统(Global Positioning System,GPS)台站组成的GPS观测网(图 1)。

      图  1  长白山天池火山水准点和GPS点分布图

      Figure 1.  Distribution of the Leveling Benchmarks and GPS Stations Surrounding Changbaishan Tianchi Volcano

      长白山天池火山的水准观测GPS水平观测资料来源于文献[16]。水准监测剖面在2002—2005年进行了4期水准观测,其中2002年9月至2003年9月的垂直位移量最大。水平位移监测是由8个GPS台站完成的。为了与2002年9月至2003年9月的垂直位移进行联合反演,本文采用2002—2003年的数据。

      首先采用F-J方法。选取不同的MCMC采样参数得到的F-J方法反演结果如表 3所示,表 3中各种情况与表 1中一致。其中,参考结果是经过多次实验后得到的选取了合适的采样参数的反演结果,其他结果与参考结果相比调整了一个参数(参数调整方法与前述模拟实验相同)。表 3同样表明,改变MCMC采样参数中的采样点数、舍弃初始采样点数和随机游走步长对反演参数结果影响比较明显。

      表 3  实测数据不同条件下F-J方法求出的参数估计值

      Table 3.  Estimation and Variance of Parameters by Using the F-J Method of Experimental Data at Different Conditions

      情况 D/km x0/km y0/km σ12/mm2 σ22/mm2 ΔV/103 m3
      参考结果 7.533 4±0.629 3 2.030 1±0.978 1 0.836 1±0.659 0 6.925 9±2.671 2 6.302 6±2.391 3 22.250 2±3.278 4
      7.554 7±0.548 5 2.029 6±0.880 0 0.839 6±0.514 2 6.921 6±2.421 8 6.299 7±1.915 7 22.249 5±2.792 2
      7.896 0±0.748 2 1.655 7±1.195 6 1.035 7±0.765 3 8.491 7±3.045 4 6.787 5±2.915 2 23.925 5±3.940 4
      7.034 1±0.764 7 2.170 5±1.261 6 0.559 1±0.697 5 5.703 7±2.285 3 4.535 8±2.651 5 21.965 6±3.944 5
      7.508 2±0.531 9 2.559 2±0.779 5 0.748 3±0.472 8 6.980 5±2.566 9 4.633 3±1.916 7 22.574 3±3.056 3
      马尔可夫链发散

      然后使用迭代最小二乘计算方法反演,初始值在一定的范围内随机选取。参数Dx0y0σ12σ22、ΔV的取值范围分别为[0,15]、[-5,5]、[-5,5]、[1,16]、[1,16]和[10,30]。实验中共随机生成了30组迭代初始值,同样发现,有23组初始值不能使迭代收敛,4组收敛到错误的结果,3组收敛到了正确的结果。因此,在可能的取值范围内随机取值,通常不能保证迭代最小二乘方法收敛到正确的结果。

      本文同时分别采用以上5组F-J方法计算得到的结果为初始值进行迭代最小二乘反演,每组都收敛到相同的结果,则参数D迭代收敛结果为(7.854 8±0.548 5) km;x0=(1.818 4±1.025 2) km;y0=(1.075 5±1.025 0) km;σ12=7.145 4 mm2σ22=6.624 8 mm2;ΔV=(23.294 6±3.592 7)×103 m3

      文献[17]利用2002—2003年的水准和GPS数据,基于Mogi模型反演了长白山地区火山的岩浆囊参数。文献[11]也使用相同的数据分别进行了最小二乘和总体最小二乘(total least square,TLS)反演,其结果与本文的结果对比如表 4所示。

      表 4  各研究方法结果的参数比较

      Table 4.  Parameters Comparison of Different Research Methods

      方法 D/km x0/km y0/km ΔV/103 m3
      本文F-J方法 7.533 4±0.629 3 2.030 1±0.978 1 0.836 1±0.659 0 22.250 2±3.278 4
      基于F-J方法解的迭代最小二乘方法 7.854 8±0.548 5 1.818 4±1.025 2 1.075 5±1.025 0 23.294 6±3.592 7
      文献[17]方法 6.9 2.5 0.5 5.2
      文献[11]最小二乘方法 10.335 4±1.084 9 3.005 9±0.930 5 3.936 9±1.051 2 22.463 2±1.870 2
      文献[11]TLS方法 5.606 9±0.020 9 2.257 4±0.021 0 1.356 0±0.020 9 12.814 0±0.021 0

      表 4表明不同研究得到的参数值之间存在一定的差异,可能是所采用的反演方法不同造成的。相较于文献[17],文献[11]采用TLS方法得到的结果与本文基于F-J方法解的迭代最小二乘方法的岩浆囊部分参数(岩浆囊中心位置(x0, y0))更为接近。其原因可能是这两类方法都顾及了模型参数先验误差,但相较于TLS解,以F-J方法的解作为迭代最小二乘反演的初值更能改善最小二乘的反演结果。

    • 为验证本文方法,根据文献[17]进行了模拟实验。首先采用不同的MCMC采样参数进行F-J方法反演,结果如表 5所示。表 5中情况①~③与表 1同,情况④表示1.5倍游走步长,情况⑤表示3倍游走步长。

      表 5  F-J方法求出的参数估计值及误差

      Table 5.  Estimation and Variance of Parameters by Using the F-J Method

      情况 D1/km D2/km σ12/(mm·a-1)2 σ22/(mm·a-1)2 σ32/(mm·a-1)2 1/(mm·a-1) 2/(mm·a-1)
      真值 10 10 1.00 6.25 20.25 25 15
      参考结果 10.004±1.029 9.940±1.579 1.023±0.074 6.450±0.263 20.341±0.590 24.712±1.218 15.151±1.217
      10.004±0.896 9.940±1.316 1.024±0.060 6.450±0.256 20.341±0.496 24.712±1.004 15.151±0.985
      10.591±1.144 8.855±1.816 1.027±0.082 6.449±0.254 20.375±0.595 25.348±1.331 14.552±1.352
      10.829±1.463 8.582±1.826 1.012±0.113 6.737±0.391 20.529±0.963 26.130±1.810 13.779±1.828
      9.561±0.705 10.399±1.036 1.032±0.057 6.457±0.174 20.532±0.384 24.058±0.724 15.774±0.758
      马尔可夫链发散

      然后使用迭代最小二乘计算,初始值在一定的范围内随机选取。参数D1D2σ12σ22σ3212的取值范围分别为[0,20]、[0,20]、[1,5]、[3,10]、[10,30]、[0,40]和[0,40]。实验中共随机生成了30组迭代初始值,经实验发现有18组初始值不能收敛,4组收敛错误,8组收敛正确。

      最后,分别采用以上5组F-J方法计算得到的结果为初始值进行最小二乘计算,每组都收敛到相同的结果,如表 6所示。

      表 6  以F-J方法的解为初值的迭代最小二乘结果

      Table 6.  Results of the Iterative Least Squares Method with Initial Values Resulted from the F-J Method

      情况 D1/km D2/km σ12/(mm·a-1)2 σ22/(mm·a-1)2 σ32/(mm·a-1)2 1/(mm·a-1) 2/(mm·a-1)
      真值 10 10 1.0 6.25 20.25 25 15
      迭代收敛结果 10.002±1.031 9.948±1.616 1.024±0.076 6.455±0.262 20.349±0.588 24.717±1.228 15.154±1.220
    • 本文使用F-J方法需要非常大的计算量,在计算时使用Intel的Xeon E3-1231 v3处理器,完成100万次采样需要20 h左右。而迭代最小二乘方法收敛只需几秒,基于F-J线性-非线性模型解的迭代最小二乘方法则需约20 h,极大地提升了计算效率。F-J方法属于随机反演方法,具有全局最优性。但F-J方法在进行MCMC采样时,采样经验参数会显著影响反演结果,而实际计算中通常很难直接确定合适的采样参数。如果单纯使用F-J方法,需要多次调整各个采样参数,耗时较长。

      最小二乘法原理简单,计算量远少于F-J方法,但在初始值选取不当的情况下,通常会收敛到局部最优解或者迭代不收敛。本文算例中,当F-J方法的采样参数选取不合适但马尔可夫链没有发散时,以F-J方法反演结果作为迭代最小二乘方法的初始值可以得到一致收敛结果。

      在实际反演中,几乎没有文献分析讨论F-J方法中的采样参数对反演结果的影响。本文提出的基于F-J方法解的迭代最小二乘方法可以有效解决此问题。如果无法选择合适的待估参数的近似值,则可以先用F-J方法进行初步反演。只要MCMC采样不发散,就能得到近似的反演结果。以此初步反演结果作为迭代最小二乘方法的初始值,可以得到一致收敛的全局稳定最优解。

    • F-J方法可以从理论上求出模型参数的概率分布,然后通过MCMC采样估计参数并评定精度,然而MCMC采样参数的选择不当会影响F-J方法的效果。F-J反演方法需要试探性地选择MCMC参数,计算量非常大。F-J方法的初步反演结果可以作为迭代最小二乘方法的初始值。基于F-J线性-非线性模型解的迭代最小二乘方法,可以充分利用F-J方法的全局最优特性和最小二乘方法计算量少的优点,在线性-非线性模型参数反演方面有很好的应用前景。

参考文献 (17)

目录

    /

    返回文章
    返回