留言板

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

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

火山Mogi模型反演的总体最小二乘联合平差方法

王乐洋 余航

王乐洋, 余航. 火山Mogi模型反演的总体最小二乘联合平差方法[J]. 武汉大学学报 ● 信息科学版, 2018, 43(9): 1333-1341. doi: 10.13203/j.whugis20160469
引用本文: 王乐洋, 余航. 火山Mogi模型反演的总体最小二乘联合平差方法[J]. 武汉大学学报 ● 信息科学版, 2018, 43(9): 1333-1341. doi: 10.13203/j.whugis20160469
WANG Leyang, YU Hang. Application of Total Least Squares Joint Adjustment to Volcano Inversion of Mogi Model[J]. Geomatics and Information Science of Wuhan University, 2018, 43(9): 1333-1341. doi: 10.13203/j.whugis20160469
Citation: WANG Leyang, YU Hang. Application of Total Least Squares Joint Adjustment to Volcano Inversion of Mogi Model[J]. Geomatics and Information Science of Wuhan University, 2018, 43(9): 1333-1341. doi: 10.13203/j.whugis20160469

火山Mogi模型反演的总体最小二乘联合平差方法

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

国家自然科学基金 41664001

国家自然科学基金 41204003

江西省杰出青-人才资助计划 20162BCB23050

国家重点研发计划 2016YFB0501405

江西省教育厅科技项目 GJJ150595

详细信息
    作者简介:

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

  • 中图分类号: P207

Application of Total Least Squares Joint Adjustment to Volcano Inversion of Mogi Model

Funds: 

The National Natural Science Foundation of China 41664001

The National Natural Science Foundation of China 41204003

Support Program for Outstanding Youth Talents in Jiangxi Province 20162BCB23050

the National Key R & D Program of China 2016YFB0501405

Science and Technology Project of the Education Department of Jiangxi Province GJJ150595

More Information
    Author Bio:

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

图(6) / 表(3)
计量
  • 文章访问数:  1620
  • HTML全文浏览量:  82
  • PDF下载量:  237
  • 被引次数: 0
出版历程
  • 收稿日期:  2017-07-06
  • 刊出日期:  2018-09-05

火山Mogi模型反演的总体最小二乘联合平差方法

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

    国家自然科学基金 41664001

    国家自然科学基金 41204003

    江西省杰出青-人才资助计划 20162BCB23050

    国家重点研发计划 2016YFB0501405

    江西省教育厅科技项目 GJJ150595

    作者简介:

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

  • 中图分类号: P207

摘要: 针对垂直位移与水平位移的Mogi模型,提出采用总体最小二乘联合(total least squares joint,TLS-J)平差方法进行求解。该方法可同时顾及联合平差函数模型中观测向量与系数矩阵的误差项,且采用3种判别函数最小化法确定相对权比,用以权衡垂直位移与水平位移观测数据在联合求解过程中所占的比重。针对平差过程中出现的病态问题,结合L曲线法确定岭参数。通过实际算例,系统研究了总体最小二乘联合平差方法在长白山天池火山Mogi模型反演中的应用。研究结果表明,以判别函数为$\sum\limits_{i=1}^{n1}{\left| {{{\hat{\bar{e}}}}_{1i}} \right|}+\sum\limits_{j=1}^{n2}{\left| {{{\hat{\bar{e}}}}_{2j}} \right|}$的函数最小化能获得合理的压力源参数估值结果和相对权比大小,具有一定的实际参考价值。

English Abstract

王乐洋, 余航. 火山Mogi模型反演的总体最小二乘联合平差方法[J]. 武汉大学学报 ● 信息科学版, 2018, 43(9): 1333-1341. doi: 10.13203/j.whugis20160469
引用本文: 王乐洋, 余航. 火山Mogi模型反演的总体最小二乘联合平差方法[J]. 武汉大学学报 ● 信息科学版, 2018, 43(9): 1333-1341. doi: 10.13203/j.whugis20160469
WANG Leyang, YU Hang. Application of Total Least Squares Joint Adjustment to Volcano Inversion of Mogi Model[J]. Geomatics and Information Science of Wuhan University, 2018, 43(9): 1333-1341. doi: 10.13203/j.whugis20160469
Citation: WANG Leyang, YU Hang. Application of Total Least Squares Joint Adjustment to Volcano Inversion of Mogi Model[J]. Geomatics and Information Science of Wuhan University, 2018, 43(9): 1333-1341. doi: 10.13203/j.whugis20160469
  • Mogi模型是日本学者Kiyoo Mogi依据Yamakawa理论公式,将岩浆压力源与火山地表形变建立联系的一种模型[1]。应用该模型的前提假设是:将岩浆压力源置于弹性半空间中,当源半径远小于源的深度时,火山岩浆的活动量可以看作是等效球状源。该模型已用来模拟和解释大量的火山区地表变形,应用火山区的形变测量数据来反演岩浆压力源参数,对火山的危险性评价有很大的意义[2]。在中国,Mogi模型主要应用于长白山天池火山[3-6]以及腾冲火山[7-8]研究。目前, 针对火山Mogi模型的压力源参数反演主要采用最小二乘方法,并取得了丰富的成果。然而,对于Mogi模型而言,其系数矩阵不仅含有待求的压力源参数,同时还存在含有误差的坐标数据,其实属于非线性的变量含误差模型[9-11]。近年来,总体最小二乘平差方法研究较多[12-17],而针对Mogi模型反演方面的文献相对较少。文献[9]采用结构总体最小二乘(structured total least squares,STLS)的L1范数,并结合1982-1985年的垂直形变数据计算了意大利南部Campi Flegrei火山压力源参数;文献[10]采用总体最小二乘联合平差方法反演了长白山天池火山压力源参数,但文献并未顾及联合反演中垂直形变数据与水平形变数据之间的权比关系,即视各类数据等权比。针对求解Mogi模型存在系数矩阵病态的问题,文献[11]在顾及系数矩阵误差的情况下,采用虚拟观测法、谱修正迭代法和共轭梯度法探讨了其在火山Mogi模型中的应用。

    本文采用总体最小二乘联合平差(total least squares joint, TLS-J)方法,研究该方法在长白山天池火山Mogi模型反演中的应用。针对该模型非线性的特点,对Mogi模型进行线性化,给出观测向量与系数矩阵协因数阵的计算公式。平差时,由于模型法矩阵病态性严重,故采用L曲线法确定岭估计中的岭参数,以抑制平差过程中出现的病态性。采用单压力源模型,运用判别函数最小化法来确定垂直形变数据与水平形变数据间的权比,给出了针对火山Mogi模型加权总体最小二乘联合平差的迭代计算公式。

    • 由于岩浆房内外的压力不同,存在一定程度的压力差,从而会导致岩浆房体积发生变化。由膨胀压力源引起地表垂直位移和水平位移与等效压力源参数之间的关系式为[1]

      $$ \Delta h = \frac{{3K + 4\mu }}{{2{\rm{ \mathsf{ π} }}\left( {3K + \mu } \right)}} \cdot \frac{D}{{{{\left( {{D^2} + {r^2}} \right)}^{3/2}}}}\Delta V $$ (1)
      $$ \Delta r = \frac{{3K + 4\mu }}{{2{\rm{ \mathsf{ π} }}\left( {3K + \mu } \right)}} \cdot \frac{r}{{{{\left( {{D^2} + {r^2}} \right)}^{3/2}}}}\Delta V $$ (2)

      式中,Δh为地表垂直位移;Δr为地表水平位移;K为体积弹性模量;D为源的中心深度; ΔV为岩浆房体增量; μ为剪切模量;r为地表径向距离,且当r=0时垂直位移取得最大值。

      当地壳为泊松介质时,且顾及岩浆房体积的膨胀是半径为R的等效球体,式(1)与式(2)可简化为[7]

      $$ \Delta h = \frac{{{R^3}D}}{{{{\left( {{D^2} + {r^2}} \right)}^{3/2}}}} $$ (3)
      $$ \Delta r = \frac{{{R^3}r}}{{{{\left( {{D^2} + {r^2}} \right)}^{3/2}}}} $$ (4)

      在实际的火山形变监测中,观测到的水平位移往往是东向分量与北向分量的形式,如GPS得到的形变速度场,因此在水平位移反演时常采用如下模型:

      $$ {U_x} = \frac{{{R^3}\left( {x - {x_0}} \right)}}{{{{\left( {{D^2} + {r^2}} \right)}^{3/2}}}} $$ (5)
      $$ {U_y} = \frac{{{R^3}\left( {y - {y_0}} \right)}}{{{{\left( {{D^2} + {r^2}} \right)}^{3/2}}}} $$ (6)

      式中,Uxx方向位移;Uyy方向位移;(x0, y0)为压力源在平面上的投影坐标。

      i个监测点(xi, yi)到源之间的地表径向距离为r2=(xi-x0)2+(yi-y0)2。若有n个监测点,由式(3)可得:

      $$ \begin{array}{*{20}{c}} {\Delta {h_i} + {e_i} = \frac{{{R^3}D}}{{{{\left( {{D^2} + {{\left( {{x_i} - {x_0}} \right)}^2} + {{\left( {{y_i} - {y_0}} \right)}^2}} \right)}^{3/2}}}} = }\\ {{f_i}\left( \mathit{\boldsymbol{\theta }} \right)} \end{array} $$ (7)

      式中,ei为对应第i个地表垂直位移观测数据Δhi的误差; i=1, 2…nθ=[R  D  x0  y0]T

      对式(7)利用泰勒级数在初值θ0处展开,并舍去二次项,得如下矩阵形式[10]

      $$ \begin{array}{*{20}{c}} {\underbrace {\left[ {\begin{array}{*{20}{c}} {\Delta {h_1}}\\ \vdots \\ {\Delta {h_n}} \end{array}} \right]}_{{\mathit{\boldsymbol{y}}_1}} + \underbrace {\left[ {\begin{array}{*{20}{c}} {{e_1}}\\ \vdots \\ {{e_n}} \end{array}} \right]}_{{\mathit{\boldsymbol{e}}_{{\mathit{\boldsymbol{y}}_{_1}}}}} = }\\ {\underbrace {\left[ {\begin{array}{*{20}{c}} {\frac{{\partial {f_1}}}{{\partial R}}}&{\frac{{\partial {f_1}}}{{\partial D}}}&{\frac{{\partial {f_1}}}{{\partial {x_0}}}}&{\frac{{\partial {f_1}}}{{\partial {y_0}}}}\\ \vdots&\vdots&\vdots&\vdots \\ {\frac{{\partial {f_n}}}{{\partial R}}}&{\frac{{\partial {f_n}}}{{\partial D}}}&{\frac{{\partial {f_n}}}{{\partial {x_0}}}}&{\frac{{\partial {f_n}}}{{\partial {y_0}}}} \end{array}} \right]}_{{\mathit{\boldsymbol{B}}_1}}\left[ {\begin{array}{*{20}{c}} {{\rm{d}}R}\\ {{\rm{d}}D}\\ {{\rm{d}}{x_0}}\\ {{\rm{d}}{y_0}} \end{array}} \right] + \underbrace {\left[ {\begin{array}{*{20}{c}} {{f_1}\left( {{\mathit{\boldsymbol{\theta }}_0}} \right)}\\ \vdots \\ {{f_n}\left( {{\mathit{\boldsymbol{\theta }}_0}} \right)} \end{array}} \right]}_{{\mathit{\boldsymbol{y}}_{01}}}} \end{array} $$ (8)

      同理,对于水平位移式(5)和式(6),若有n个监测点,则有:

      $$ \begin{array}{*{20}{c}} {{U_{{x_i}}} + {e_{{U_{{x_i}}}}} = }\\ {\frac{{{R^3}\left( {x - {x_0}} \right)}}{{{{\left( {{D^2} + {{\left( {{x_i} - {x_0}} \right)}^2} + {{\left( {{y_i} - {y_0}} \right)}^2}} \right)}^{3/2}}}} = {\varphi _{{x_i}}}\left( \mathit{\boldsymbol{\theta }} \right)} \end{array} $$ (9)
      $$ \begin{array}{*{20}{c}} {{U_{{y_i}}} + {e_{{U_{{y_i}}}}} = }\\ {\frac{{{R^3}\left( {y - {y_0}} \right)}}{{{{\left( {{D^2} + {{\left( {{x_i} - {x_0}} \right)}^2} + {{\left( {{y_i} - {y_0}} \right)}^2}} \right)}^{3/2}}}} = {\varphi _{{y_i}}}\left( \mathit{\boldsymbol{\theta }} \right)} \end{array} $$ (10)

      式中, eUxieUyi分别为对应第i个地表水平方向位移UxUy的误差。

      分别对式(9)、式(10)进行泰勒展开至一次项,写成矩阵形式[10]

      $$ \begin{array}{*{20}{c}} {\underbrace {\left[ {\begin{array}{*{20}{c}} {{U_{{x_1}}}}\\ {{U_{{y_1}}}}\\ \vdots \\ {{U_{{x_n}}}}\\ {{U_{{y_n}}}} \end{array}} \right]}_{{\mathit{\boldsymbol{y}}_2}} + \underbrace {\left[ {\begin{array}{*{20}{c}} {{e_{{U_{{x_1}}}}}}\\ {{e_{{U_{{y_1}}}}}}\\ \vdots \\ {{e_{{U_{{x_n}}}}}}\\ {{e_{{U_{{y_n}}}}}} \end{array}} \right]}_{{\mathit{\boldsymbol{e}}_{{\mathit{\boldsymbol{y}}_{_2}}}}} = }\\ {\underbrace {\left[ {\begin{array}{*{20}{c}} {\frac{{\partial {\varphi _{{x_1}}}}}{{\partial R}}}&{\frac{{\partial {\varphi _{{x_1}}}}}{{\partial D}}}&{\frac{{\partial {\varphi _{{x_1}}}}}{{\partial {x_0}}}}&{\frac{{\partial {\varphi _{{x_1}}}}}{{\partial {y_0}}}}\\ {\frac{{\partial {\varphi _{{y_1}}}}}{{\partial R}}}&{\frac{{\partial {\varphi _{{y_1}}}}}{{\partial D}}}&{\frac{{\partial {\varphi _{{y_1}}}}}{{\partial {x_0}}}}&{\frac{{\partial {\varphi _{{y_1}}}}}{{\partial {y_0}}}}\\ \vdots &{}&{}& \vdots \\ {\frac{{\partial {\varphi _{{x_n}}}}}{{\partial R}}}&{\frac{{\partial {\varphi _{{x_n}}}}}{{\partial D}}}&{\frac{{\partial {\varphi _{{x_n}}}}}{{\partial {x_0}}}}&{\frac{{\partial {\varphi _{{x_n}}}}}{{\partial {y_0}}}}\\ {\frac{{\partial {\varphi _{{y_n}}}}}{{\partial R}}}&{\frac{{\partial {\varphi _{{y_n}}}}}{{\partial D}}}&{\frac{{\partial {\varphi _{{y_n}}}}}{{\partial {x_0}}}}&{\frac{{\partial {\varphi _{{y_n}}}}}{{\partial {y_0}}}} \end{array}} \right]}_{{\mathit{\boldsymbol{B}}_2}}\left[ {\begin{array}{*{20}{c}} {{\rm{d}}R}\\ {{\rm{d}}D}\\ {{\rm{d}}{x_0}}\\ {{\rm{d}}{y_0}} \end{array}} \right] + \underbrace {\left[ \begin{array}{l} {\varphi _{{x_1}}}\left( {{\mathit{\boldsymbol{\theta }}_0}} \right)\\ {\varphi _{{y_1}}}\left( {{\mathit{\boldsymbol{\theta }}_0}} \right)\\ \vdots \\ {\varphi _{{x_n}}}\left( {{\mathit{\boldsymbol{\theta }}_0}} \right)\\ {\varphi _{{y_n}}}\left( {{\mathit{\boldsymbol{\theta }}_0}} \right) \end{array} \right]}_{{\mathit{\boldsymbol{y}}_{02}}}} \end{array} $$ (11)
    • 由线性化后的式(8)和式(11)可知,垂直位移与水平位移观测值对应的系数矩阵B1B2中不仅存在含有随机误差的坐标数据,同时还存在待估压力源参数;压力源参数的变化必将引起系数矩阵B1B2的结构性扰动,同时线性化迭代解法如Gauss-Newton法和Levenbeg-Marquardt法需要求解系数矩阵的雅可比矩阵,线性化必然使得求导得到的系数矩阵含有误差,应将此误差与观测值误差同时考虑[10]。因此,要获得Mogi模型火山压力源参数,实则是求解非线性总体最小二乘联合平差问题[18-20]。不同于已有的非线性总体最小二乘平差方法,火山Mogi模型的联合求解问题除了需顾及系数矩阵的扰动影响和相对权比的确定外,还需顾及火山形变反演过程中出现病态或奇异的情况。Mogi模型线性化后系数矩阵的条件数较大,一般在106量级[5],可通过岭估计方法来处理法矩阵求逆不稳定的问题。

      总体最小二乘联合平差的函数模型为:

      $$ \left\{ \begin{array}{l} \left( {{\mathit{\boldsymbol{B}}_1} + {\mathit{\boldsymbol{E}}_{{\mathit{\boldsymbol{B}}_1}}}} \right)\mathit{\boldsymbol{X}} = {\mathit{\boldsymbol{l}}_1} + {\mathit{\boldsymbol{e}}_{{\mathit{\boldsymbol{l}}_1}}}\\ \left( {{\mathit{\boldsymbol{B}}_2} + {\mathit{\boldsymbol{E}}_{{\mathit{\boldsymbol{B}}_2}}}} \right)\mathit{\boldsymbol{X}} = {\mathit{\boldsymbol{l}}_2} + {\mathit{\boldsymbol{e}}_{{\mathit{\boldsymbol{l}}_2}}} \end{array} \right. $$ (12)

      式中,l1=y1-y01, y1n1×1维垂直位移观测值;el1为对应垂直位移观测方程的随机误差;l2=y2-y02, y2n2×1维水平位移观测值;el2为对应水平位移观测方程的随机误差;B1n1×4维垂直位移观测方程列满秩系数矩阵, EB1为其误差矩阵;B2n2×4维水平位移观测方程列满秩系数矩阵,EB2为其误差矩阵;X=[dR  dD  dx0  dy0]T为4×1待估参数向量。

      联合平差的随机模型为:

      $$ \left\{ \begin{array}{l} \left[ {\begin{array}{*{20}{c}} {{\mathit{\boldsymbol{e}}_{{\mathit{\boldsymbol{l}}_1}}}}\\ {{\mathit{\boldsymbol{e}}_{{\mathit{\boldsymbol{B}}_1}}}} \end{array}} \right] = \left[ {\begin{array}{*{20}{c}} {{\mathit{\boldsymbol{e}}_{{\mathit{\boldsymbol{l}}_1}}}}\\ {{\rm{vec}}\left( {{\mathit{\boldsymbol{E}}_{{\mathit{\boldsymbol{B}}_1}}}} \right)} \end{array}} \right] \sim N\left( {\left[ {\begin{array}{*{20}{c}} {\bf{0}}\\ {\bf{0}} \end{array}} \right],\sigma _1^2\left[ {\begin{array}{*{20}{c}} {{\mathit{\boldsymbol{Q}}_{{\mathit{\boldsymbol{l}}_1}}}}&{\bf{0}}\\ {\bf{0}}&{{\mathit{\boldsymbol{Q}}_{{\mathit{\boldsymbol{B}}_1}}}} \end{array}} \right]} \right)\\ \left[ {\begin{array}{*{20}{c}} {{\mathit{\boldsymbol{e}}_{{\mathit{\boldsymbol{l}}_2}}}}\\ {{\mathit{\boldsymbol{e}}_{{\mathit{\boldsymbol{B}}_2}}}} \end{array}} \right] = \left[ {\begin{array}{*{20}{c}} {{\mathit{\boldsymbol{e}}_{{\mathit{\boldsymbol{l}}_2}}}}\\ {{\rm{vec}}\left( {{\mathit{\boldsymbol{E}}_{{\mathit{\boldsymbol{B}}_2}}}} \right)} \end{array}} \right] \sim N\left( {\left[ {\begin{array}{*{20}{c}} {\bf{0}}\\ {\bf{0}} \end{array}} \right],\sigma _2^2\left[ {\begin{array}{*{20}{c}} {{\mathit{\boldsymbol{Q}}_{{\mathit{\boldsymbol{l}}_2}}}}&{\bf{0}}\\ {\bf{0}}&{{\mathit{\boldsymbol{Q}}_{{\mathit{\boldsymbol{B}}_2}}}} \end{array}} \right]} \right) \end{array} \right. $$ (13)

      式中,eB1eB2分别是将矩阵EB1EB2按列拉直得到的列向量;σ12σ22分别为垂直及水平观测值的方差分量;Ql1QB1Ql2QB2分别为对应观测向量误差和系数矩阵误差的协因数阵;vec(·)为矩阵拉直运算。

      进行移项变换,函数模型式(12)等价于:

      $$ \left\{ \begin{array}{l} {{\mathit{\boldsymbol{\bar e}}}_1} = {\mathit{\boldsymbol{B}}_1}\mathit{\boldsymbol{X}} - {\mathit{\boldsymbol{l}}_1}\\ {{\mathit{\boldsymbol{\bar e}}}_2} = {\mathit{\boldsymbol{B}}_2}\mathit{\boldsymbol{X}} - {\mathit{\boldsymbol{l}}_2} \end{array} \right. $$ (14)

      ē1=[In1-XTIn1]$\left[ \begin{array}{l} {\mathit{\boldsymbol{e}}_{{\mathit{\boldsymbol{l}}_1}}}\\ {\mathit{\boldsymbol{e}}_{{\mathit{\boldsymbol{B}}_1}}} \end{array} \right]$, ē2=[In2-XTIn2]$\left[ \begin{array}{l} {\mathit{\boldsymbol{e}}_{{\mathit{\boldsymbol{l}}_2}}}\\ {\mathit{\boldsymbol{e}}_{{\mathit{\boldsymbol{B}}_2}}} \end{array} \right]$,其中,ⓧ为Kronecker积。根据式(13),由协因数传播律容易得到ē1ē2的协因数阵为:

      $$ \begin{array}{*{20}{c}} {{\mathit{\boldsymbol{Q}}_{{{\mathit{\boldsymbol{\bar e}}}_1}}} = }\\ {\left[ {{\mathit{\boldsymbol{I}}_{{\mathit{\boldsymbol{n}}_1}}} - {\mathit{\boldsymbol{X}}^{\rm{T}}} \otimes {\mathit{\boldsymbol{I}}_{{\mathit{\boldsymbol{n}}_1}}}} \right]\left[ {\begin{array}{*{20}{c}} {{\mathit{\boldsymbol{Q}}_{{\mathit{\boldsymbol{l}}_1}}}}&{\bf{0}}\\ {\bf{0}}&{{\mathit{\boldsymbol{Q}}_{{\mathit{\boldsymbol{B}}_1}}}} \end{array}} \right]\left[ {\begin{array}{*{20}{c}} {{\mathit{\boldsymbol{I}}_{{\mathit{\boldsymbol{n}}_1}}}}\\ { - {\mathit{\boldsymbol{X}}^{\rm{T}}} \otimes {\mathit{\boldsymbol{I}}_{{\mathit{\boldsymbol{n}}_1}}}} \end{array}} \right] = }\\ {{\mathit{\boldsymbol{Q}}_{{\mathit{\boldsymbol{l}}_1}}} + \left( {{\mathit{\boldsymbol{X}}^{\rm{T}}} \otimes {\mathit{\boldsymbol{I}}_{{\mathit{\boldsymbol{n}}_1}}}} \right){\mathit{\boldsymbol{Q}}_{{\mathit{\boldsymbol{B}}_1}}}\left( {\mathit{\boldsymbol{X}} \otimes {\mathit{\boldsymbol{I}}_{{\mathit{\boldsymbol{n}}_1}}}} \right)} \end{array} $$ (15)
      $$ \begin{array}{*{20}{c}} {{\mathit{\boldsymbol{Q}}_{{{\mathit{\boldsymbol{\bar e}}}_2}}} = }\\ {\left[ {{\mathit{\boldsymbol{I}}_{{\mathit{\boldsymbol{n}}_2}}} - {\mathit{\boldsymbol{X}}^{\rm{T}}} \otimes {\mathit{\boldsymbol{I}}_{{\mathit{\boldsymbol{n}}_2}}}} \right]\left[ {\begin{array}{*{20}{c}} {{\mathit{\boldsymbol{Q}}_{{\mathit{\boldsymbol{l}}_2}}}}&{\bf{0}}\\ {\bf{0}}&{{\mathit{\boldsymbol{Q}}_{{\mathit{\boldsymbol{B}}_2}}}} \end{array}} \right]\left[ {\begin{array}{*{20}{c}} {{\mathit{\boldsymbol{I}}_{{\mathit{\boldsymbol{n}}_2}}}}\\ { - {\mathit{\boldsymbol{X}}^{\rm{T}}} \otimes {\mathit{\boldsymbol{I}}_{{\mathit{\boldsymbol{n}}_2}}}} \end{array}} \right] = }\\ {{\mathit{\boldsymbol{Q}}_{{\mathit{\boldsymbol{l}}_2}}} + \left( {{\mathit{\boldsymbol{X}}^{\rm{T}}} \otimes {\mathit{\boldsymbol{I}}_{{\mathit{\boldsymbol{n}}_2}}}} \right){\mathit{\boldsymbol{Q}}_{{\mathit{\boldsymbol{B}}_2}}}\left( {\mathit{\boldsymbol{X}} \otimes {\mathit{\boldsymbol{I}}_{{\mathit{\boldsymbol{n}}_2}}}} \right)} \end{array} $$ (16)

      由Tikhonov正则化理论[21],Mogi模型反演的总体最小二乘联合平差估计准则为[22]

      $$ \mathit{\Phi } = \lambda \mathit{\boldsymbol{\hat {\bar e}}}_1^{\rm{T}}\mathit{\boldsymbol{Q}}_{{{\mathit{\boldsymbol{\hat {\bar e}}}}_1}}^{ - 1}{{\mathit{\boldsymbol{\hat {\bar e}}}}_1} + \left( {1 - \lambda } \right)\mathit{\boldsymbol{\hat {\bar e}}}_2^{\rm{T}}Q_{{{\mathit{\boldsymbol{\hat {\bar e}}}}_2}}^{ - 1}{{\mathit{\boldsymbol{\hat {\bar e}}}}_2} + \alpha {{\mathit{\boldsymbol{\hat X}}}^{\rm{T}}}\mathit{\boldsymbol{\hat X}} = \min $$ (17)

      式中,α为岭参数;λ为相对权比(0 < λ < 1);ē1的平差值为${{{\mathit{\boldsymbol{\hat{\bar{e}}}}}}_{1}}$;ē2的平差值为${{{\mathit{\boldsymbol{\hat{\bar{e}}}}}}_{2}}$;X的平差值为${\mathit{\boldsymbol{\hat X}}}$; ${\mathit{\boldsymbol{Q}}_{{{{\mathit{\boldsymbol{\hat{\bar{e}}}}}}_{1}}}}$和${\mathit{\boldsymbol{Q}}_{{{{\mathit{\boldsymbol{\hat{\bar{e}}}}}}_{2}}}}$为将${\mathit{\boldsymbol{\hat X}}}$代入${\mathit{\boldsymbol{Q}}_{{{{\mathit{\boldsymbol{\hat{\bar{e}}}}}}_{1}}}}$、${\mathit{\boldsymbol{Q}}_{{{{\mathit{\boldsymbol{\hat{\bar{e}}}}}}_{2}}}}$中所得。

      由式(17)对${\mathit{\boldsymbol{\hat X}}}$求导,并令其为零得[22-23]

      $$ \begin{array}{*{20}{c}} {\mathit{\boldsymbol{\hat X}} = {{\left( {\lambda {{\mathit{\boldsymbol{\bar B}}}_1} + \left( {1 - \lambda } \right){{\mathit{\boldsymbol{\bar B}}}_2} + \alpha {\mathit{\boldsymbol{I}}_{4 \times 4}}} \right)}^{ - 1}} \times }\\ {\left( {\lambda \mathit{\boldsymbol{\bar B}}_1^{\rm{T}}Q_{{{\mathit{\boldsymbol{\hat {\bar e}}}}_1}}^{ - 1}{{\mathit{\boldsymbol{\bar l}}}_1} + \left( {1 - \lambda } \right)\mathit{\boldsymbol{\bar B}}_2^{\rm{T}}Q_{{{\mathit{\boldsymbol{\hat {\bar e}}}}_2}}^{ - 1}{{\mathit{\boldsymbol{\bar l}}}_2}} \right)} \end{array} $$ (18)

      式中,${{\mathit{\boldsymbol{\bar l}}}_1} = {\mathit{\boldsymbol{l}}_1} + {{\mathit{\boldsymbol{\hat E}}}_{{\mathit{\boldsymbol{B}}_1}}}\mathit{\boldsymbol{\hat X}}$;${{\mathit{\boldsymbol{\bar l}}}_2} = {\mathit{\boldsymbol{l}}_2} + {{\mathit{\boldsymbol{\hat E}}}_{{\mathit{\boldsymbol{B}}_2}}}\mathit{\boldsymbol{\hat X}}$。

      已有的L曲线方法是在权阵为单位阵的情况下推导的[24],若权阵不是单位阵,则可对其进行单位化[24-25]。由于矩阵$\left[ {\begin{array}{*{20}{c}} {\lambda \mathit{\boldsymbol{Q}}_{{{{\mathit{\boldsymbol{\hat{\bar{e}}}}}}_{1}}}^{ - 1}}&{\bf{0}}\\ {\bf{0}}&{\left( {1 - \lambda } \right)\mathit{\boldsymbol{Q}}_{{{{\mathit{\boldsymbol{\hat{\bar{e}}}}}}_{2}}}^{ - 1}} \end{array}} \right]$正定,可采用Cholesky分解进行单位化,得$\left[ {\begin{array}{*{20}{c}} {\lambda \mathit{\boldsymbol{Q}}_{{{{\mathit{\boldsymbol{\hat{\bar{e}}}}}}_{1}}}^{ - 1}}&{\bf{0}}\\ {\bf{0}}&{\left( {1 - \lambda } \right)\mathit{\boldsymbol{Q}}_{{{{\mathit{\boldsymbol{\hat{\bar{e}}}}}}_{2}}}^{ - 1}} \end{array}} \right]$=CTCC为上三角矩阵,则式(18)可写为:

      $$ \mathit{\boldsymbol{\hat X}} = {\left( {{{\left( {\mathit{\boldsymbol{C\bar B}}} \right)}^{\rm{T}}}\left( {\mathit{\boldsymbol{C\bar B}}} \right) + \alpha {\mathit{\boldsymbol{I}}_{4 \times 4}}} \right)^{ - 1}}{\left( {\mathit{\boldsymbol{C\bar B}}} \right)^{\rm{T}}}\mathit{\boldsymbol{C\bar l}} $$ (19)

      式中,B=[B1T  B2T]Tl=[l1T  l2T]T。令D=CB, W=Cl即实现了单位化,进而再使用L曲线法即确定岭参数。

    • 选定一个步长,让权比λ遍历整个取值区间,结合§2.1推导,将每个权比λ的取值代入式(19)中计算,可得到对应联合平差的参数估值。将参数估值代入式(14),并采用条件平差法可求得等效残差向量的估值${{{\mathit{\boldsymbol{\hat{\bar{e}}}}}}_{1}}$、${{{\mathit{\boldsymbol{\hat{\bar{e}}}}}}_{2}}$[15, 23],进而由选定的判别函数计算对应权比λ下的判别函数。最后取判别函数值最小时对应的相对权比和此相对权比下的参数估值作为最终结果。本文的火山水平与垂直位移Mogi模型联合反演实例中,选取了如下3种判别函数${\mathit{\hat{\bar{\Phi }}}}$:

      $$ \mathit{\hat {\bar \Phi }} = \mathit{\boldsymbol{\hat {\bar e}}}_1^{\rm{T}}Q_{{{\mathit{\boldsymbol{\hat {\bar e}}}}_1}}^{ - 1}{{\mathit{\boldsymbol{\hat {\bar e}}}}_1} + \mathit{\boldsymbol{\hat {\bar e}}}_2^{\rm{T}}Q_{{{\mathit{\boldsymbol{\hat {\bar e}}}}_2}}^{ - 1}{{\mathit{\boldsymbol{\hat {\bar e}}}}_2} $$
      $$ \mathit{\hat {\bar \Phi }} = \sum\limits_{i = 1}^{{n_1}} {\left| {{{\mathit{\boldsymbol{\hat {\bar e}}}}_{1i}}} \right|} + \sum\limits_{j = 1}^{{n_2}} {\left| {{{\mathit{\boldsymbol{\hat {\bar e}}}}_{2j}}} \right|} $$
      $$ \mathit{\hat {\bar \Phi }} = \lambda \mathit{\boldsymbol{\hat {\bar e}}}_1^{\rm{T}}Q_{{{\mathit{\boldsymbol{\hat {\bar e}}}}_1}}^{ - 1}{{\mathit{\boldsymbol{\hat {\bar e}}}}_1} + \left( {1 - \lambda } \right)\mathit{\boldsymbol{\hat {\bar e}}}_2^{\rm{T}}Q_{{{\mathit{\boldsymbol{\hat {\bar e}}}}_2}}^{ - 1}{{\mathit{\boldsymbol{\hat {\bar e}}}}_2} $$
    • 考虑到向量y01y02是由近似压力源参数与监测点坐标数据计算得到的近似位移向量,因此,向量y01y02的协因数阵存在[20],分别记为Qy01Qy02。由式(7)、式(9)和式(10)可知,在给定了先验压力源信息θ0后,由θ0求得的第i个监测点的位移fi(θ0)、φxi(θ0)与φyi(θ0)对xiyi分别求偏导得:

      $$ {\rm{d}}{f_{i,1}} = \frac{{\partial {f_i}\left( {{\mathit{\boldsymbol{\theta }}_0}} \right)}}{{\partial {x_i}}}\left| { = - \frac{{3D{R^3}\left( {{{\hat x}_i} - {x_0}} \right)}}{{{{\left( {{D^2} + \hat r_i^2} \right)}^{5/2}}}}} \right. $$ (20)
      $$ {\rm{d}}{f_{i,2}} = \frac{{\partial {f_i}\left( {{\mathit{\boldsymbol{\theta }}_0}} \right)}}{{\partial {y_i}}}\left| { = - \frac{{3D{R^3}\left( {{{\hat y}_i} - {y_0}} \right)}}{{{{\left( {{D^2} + \hat r_i^2} \right)}^{5/2}}}}} \right. $$ (21)
      $$ \begin{array}{*{20}{c}} {{\rm{d}}{\varphi _{{x_i},1}} = \frac{{\partial {\varphi _{{x_i}}}\left( {{\mathit{\boldsymbol{\theta }}_0}} \right)}}{{\partial {x_i}}}\left| = \right.}\\ {\frac{{{R^3}}}{{{{\left( {{D^2} + \hat r_i^2} \right)}^{3/2}}}} - \frac{{3{R^3}{{\left( {{{\hat x}_i} - {x_0}} \right)}^2}}}{{{{\left( {{D^2} + \hat r_i^2} \right)}^{5/2}}}}} \end{array} $$ (22)
      $$ {\rm{d}}{\varphi _{{x_i},2}} = \frac{{\partial {\varphi _{{x_i}}}\left( {{\mathit{\boldsymbol{\theta }}_0}} \right)}}{{\partial {y_i}}}\left| { = - \frac{{3{R^3}\left( {{{\hat x}_i} - {x_0}} \right)\left( {{{\hat y}_i} - {y_0}} \right)}}{{{{\left( {{D^2} + \hat r_i^2} \right)}^{5/2}}}}} \right. $$ (23)
      $$ {\rm{d}}{\varphi _{{y_i},1}} = \frac{{\partial {\varphi _{{y_i}}}\left( {{\mathit{\boldsymbol{\theta }}_0}} \right)}}{{\partial {x_i}}}\left| { = - \frac{{3{R^3}\left( {{{\hat x}_i} - {x_0}} \right)\left( {{{\hat y}_i} - {y_0}} \right)}}{{{{\left( {{D^2} + \hat r_i^2} \right)}^{5/2}}}}} \right. $$ (24)
      $$ \begin{array}{*{20}{c}} {{\rm{d}}{\varphi _{{y_i},2}} = \frac{{\partial {\varphi _{{y_i}}}\left( {{\mathit{\boldsymbol{\theta }}_0}} \right)}}{{\partial {y_i}}}\left| = \right.}\\ {\frac{{{R^3}}}{{{{\left( {{D^2} + \hat r_i^2} \right)}^{3/2}}}} - \frac{{3{R^3}{{\left( {{{\hat y}_i} - {y_0}} \right)}^2}}}{{{{\left( {{D^2} + \hat r_i^2} \right)}^{5/2}}}}} \end{array} $$ (25)

      式(20)~(25)中,$\left( {{{\hat x}_i}, {{\hat y}_i}} \right)$为第i个监测点的坐标估计值,可通过传统测量或者卫星定位的方式获得,且${{\hat r}_i} = \sqrt {{{\left( {{{\hat x}_i} - {x_0}} \right)}^2} + {{\left( {{{\hat y}_i} - {y_0}} \right)}^2}} $。

      若第i个监测点坐标的平面精度为(σxiσyi),则可通过协因数传播律求得对应第i个监测点垂直位移与水平位移计算值的协因数阵:

      $$ \mathit{\boldsymbol{Q}}_{{\mathit{\boldsymbol{y}}_{01}}}^i = \left[ {\begin{array}{*{20}{c}} {{\rm{d}}{f_{i,1}}}&{{\rm{d}}{f_{i,2}}} \end{array}} \right]\left[ {\begin{array}{*{20}{c}} {\sigma _{{x_i}}^2}&0\\ 0&{\sigma _{{y_i}}^2} \end{array}} \right]{\left[ {\begin{array}{*{20}{c}} {{\rm{d}}{f_{i,1}}}&{{\rm{d}}{f_{i,2}}} \end{array}} \right]^{\rm{T}}} $$ (26)
      $$ \begin{array}{*{20}{c}} {\mathit{\boldsymbol{Q}}_{{\mathit{\boldsymbol{y}}_{02}}}^i = }\\ {\left[ {\begin{array}{*{20}{c}} {{\rm{d}}{\varphi _{{x_i},1}}}&{{\rm{d}}{\varphi _{{x_i},2}}}\\ {{\rm{d}}{\varphi _{{y_i},1}}}&{{\rm{d}}{\varphi _{{y_i},2}}} \end{array}} \right]\left[ {\begin{array}{*{20}{c}} {\sigma _{{x_i}}^2}&0\\ 0&{\sigma _{{y_i}}^2} \end{array}} \right]{{\left[ {\begin{array}{*{20}{c}} {{\rm{d}}{\varphi _{{x_i},1}}}&{{\rm{d}}{\varphi _{{x_i},2}}}\\ {{\rm{d}}{\varphi _{{y_i},1}}}&{{\rm{d}}{\varphi _{{y_i},2}}} \end{array}} \right]}^{\rm{T}}}} \end{array} $$ (27)

      式(26)、(27)中,Qy01iQy02i分别为协因数阵Qy01Qy02中对应第i个监测点的子矩阵,联合所有监测点的精度信息便可求得相应Qy01Qy02,因此,Ql1=Qy1+Qy01Ql2=Qy2+Qy02,其中,Qy1Qy2对应水平位移和垂直位移观测值的协因数阵。

      同理,对于确定系数矩阵的协因数阵QB1QB2,由于式(8)与式(11)中B1B2各元素同样为监测点坐标的函数,则系数矩阵中各元素需分别对xiyi求偏导,并结合协因数传播律求得各系数阵的协因数阵[22]

    • 外循环:

      1)给定压力源参数的先验信息θ0,迭代计数i=0;

      2) λ=0.01;

      3) 设置迭代计数i=i+1;内循环:①计算${\mathit{\boldsymbol{Q}}_{{{{\mathit{\boldsymbol{\hat{\bar{e}}}}}}_{1}}}}$和${\mathit{\boldsymbol{Q}}_{{{{\mathit{\boldsymbol{\hat{\bar{e}}}}}}_{2}}}}$;②更新B1B2l1l2;③Cholesky分解得$\left[ {\begin{array}{*{20}{c}} {\lambda \mathit{\boldsymbol{Q}}_{{{{\mathit{\boldsymbol{\hat{\bar{e}}}}}}_{1}}}^{ - 1}}&{\bf{0}}\\ {\bf{0}}&{\left( {1 - \lambda } \right)\mathit{\boldsymbol{Q}}_{{{{\mathit{\boldsymbol{\hat{\bar{e}}}}}}_{2}}}^{ - 1}} \end{array}} \right]$=CTC;④采用L曲线法确定岭参数;⑤由式(19)计算${\mathit{\boldsymbol{\hat X}}}$;⑥重复步骤①~⑤,直至‖${\mathit{\boldsymbol{\hat X}}}$‖ < δ为止(本文δ=1×10-8);

      4) λ=λ+0.01;

      5) 计算判别函数,并存储内循环步骤⑤中求得的压力源参数;

      6) 重复步骤3)~5),直至i=99时迭代结束;

      7) 将判别函数取最小时对应的压力源参数与相对权比值作为最终平差结果。

    • 长白山天池火山区域(如图 1所示)的垂直与水平位移观测资料来源于文献[6],表 1给出了以火山中心点经纬度(128.05°,42.00°)作为笛卡尔坐标系原点,经坐标转换后得到的各水准点、GPS点坐标(XY)及水平东向、水平北向和垂直方向的形变场结果(VEVNVH)[6]。相关数据也在文献[9]中得到过相关研究,其中,长白山水准观测为位于长白山天池至止道白河镇之间的一条水准测线,长白山火山GPS观测网由8个站点组成(见图 2)。

      图  1  长白山天池火山研究区域

      Figure 1.  Study Area of Changbaishan Tianchi Volcano

      表 1  GPS点与水准点坐标及形变值[6]

      Table 1.  Coordinates of GPS and Leveling Points and Corresponding Displacement[6]

      点名水准点点位观测数据(第1类数据)点名GPS点位观测数据(第2类数据)
      X/kmY/kmVH/(mm·a-1)X/kmY/kmVE/(mm·a-1)VN/(mm·a-1)
      IS010.3025.180P00.916.34-7.182 243.630 1
      IS19.8423.254.04P0110.1520.357.287 09.830 8
      IS29.8420.258.07P0218.093.2219.981 00.641 6
      IS38.6718.775.03P033.64-13.12-10.951 8-9.526 7
      IS47.5717.076.19P04-4.13-5.56-26.314 7-20.486 5
      IS56.4415.709.12P05-16.111.22-20.964 73.627 0
      IS65.1513.6612.67P06-14.9416.23-13.269 611.158 1
      IS73.9812.2915.01P07-12.5424.46-2.481 71.651 2
      IS82.8010.7019.75
      IS91.559.3722.71
      IS101.067.9730.85
      IS110.896.1347.03
      IS120.514.4565.79

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

      Figure 2.  Distribution Map of GPS and Leveling Points

      地表垂直位移对火山岩浆房压力中心位置比较敏感,而对岩浆压力源形状变化不太敏感,因此,为识别源的形状和深度,宜同时对垂直形变和水平形变进行反演[26]。先利用形变观测资料进行单压力源模型反演,以文献[6]中GPS和水准得到的单压力源模型参数平差结果作为最小二乘联合(least square joint, LS-J)平差的初始值,依据LS-J平差结果作为本文总体最小二乘联合平差的初始值(最小二乘联合平差采用等权比,并采用L曲线法确定岭参数)。总体最小二乘联合平差计算过程中,令各坐标点等精度。采用的不同判别函数如表 2所示。

      表 2  各方案所采用的判别函数

      Table 2.  Discriminate Function of Each Scheme

      方案方法
      1$\mathit{\hat{\bar{\Phi }}}=\lambda \mathit{\boldsymbol{\hat{\bar{e}}}}_{1}^{\rm{T}}\mathit{\boldsymbol{Q}}_{{{{\mathit{\boldsymbol{\hat{\bar{e}}}}}}_{1}}}^{-1}{{{\mathit{\boldsymbol{\hat{\bar{e}}}}}}_{1}}+\left( 1-\lambda \right)\mathit{\boldsymbol{\hat{\bar{e}}}}_{2}^{\rm{T}}\mathit{\boldsymbol{Q}}_{{{{\mathit{\boldsymbol{\hat{\bar{e}}}}}}_{2}}}^{-1}{{{\mathit{\boldsymbol{\hat{\bar{e}}}}}}_{2}}$的判别函数最小化法
      2$\mathit{\hat{\bar{\Phi }}}=\mathit{\boldsymbol{\hat{\bar{e}}}}_{1}^{\rm{T}}\mathit{\boldsymbol{Q}}_{{{{\mathit{\boldsymbol{\hat{\bar{e}}}}}}_{1}}}^{-1}{{{\mathit{\boldsymbol{\hat{\bar{e}}}}}}_{1}}+\mathit{\boldsymbol{\hat{\bar{e}}}}_{2}^{\rm{T}}\mathit{\boldsymbol{Q}}_{{{{\mathit{\boldsymbol{\hat{\bar{e}}}}}}_{2}}}^{-1}{{{\mathit{\boldsymbol{\hat{\bar{e}}}}}}_{2}}$的判别函数最小化法
      3$\mathit{\hat{\bar{\Phi }}}=\sum\limits_{i=1}^{n1}{\left| {{{\hat{\bar{e}}}}_{1i}} \right|}+\sum\limits_{j=1}^{n2}{\left| {{{\hat{\bar{e}}}}_{2j}} \right|}$的判别函数最小化法

      为了更好地与已有文献进行对比,本文将文献[5]、文献[6]和文献[10]的结果与本文计算结果列于表 3。由表 3可知,采用不同判别函数的总体最小二乘联合平差方法(方案1~3)与最小二乘联合平差结果有一定的差别。在压力源中心坐标参数结果$\left( {{{\hat x}_0}, {{\hat y}_0}} \right)$方面,方案3结果与文献[6]、文献[10]较为接近,由于文献[5]与本文采用的坐标系不一致,本文方法是在将大地坐标投影到平面坐标后再转换到以火山压力源中心为原点的平面坐标中,因此没有可比性。然而,从压力源深度参数结果${\hat D}$与等效半径参数结果${\hat R}$来看,方案3与文献[5]的结果接近。无论是总体最小二乘联合平差还是最小二乘联合平差,压力源深度结果均在4~8 km范围内,说明了总体最小二乘联合平差结果是有效的[10]

      表 3  长白山天池火山Mogi模型压力源参数反演结果

      Table 3.  Inversion Results of Pressure Source Parameters of Mogi Model of Changbaishan Tianchishan Volcano

      方法${{{\hat x}_0}}$/kmŷ0/km${\hat D}$/km${\hat R}$/kmλ
      LS-J2.385 280.388 226.879 19174.186 71
      文献[5]-0.862 50-6.807 006.096 90155.800 00
      文献[6]2.500 000.500 006.900 00107.474 37
      文献[10]2.257 401.356 005.606 00145.166 15
      方案12.926 96-0.812 176.859 66189.693 760.99
      方案22.387 330.388 476.879 66174.201 390.50
      方案32.497 070.129 217.035 25177.737 390.88

      采用方案1的判别函数确定的权比为0.99(图 3(a)),如果不加0 < λ < 1的限制,采用方案1将退化为垂直位移数据的单一反演;方案2确定的权比大小为0.50(图 3(b)),表明水平数据与垂直位移数据在联合求解过程中等权比,没有起到权衡各类数据的作用;方案1与方案2的判别函数计算结果与文献[22-23]所得结论一致。采用方案3的判别函数最小化法最终确定的相对权比大小为0.88(图 3(c)),说明垂直位移观测数据在联合平差中的比重较之水平位移的要大,这是因为联合平差中精度较高、数据量较大的那类数据一般占据较高的比重。文献[6]数据中,由于垂直与水平位移的精度相当,但垂直位移观测数据相对水平位移数据要多,垂直位移在联合平差中的比重较大是合理的。图 4反映了总体最小二乘联合平差过程中采用L曲线法求得各次迭代的条件数,其值在94.5~97.5之间,法方程病态性得到了很好的改善。

      图  3  判别函数与相对权比关系图

      Figure 3.  Relationship Between Discriminate Function and Weight Scaling Factor

      图  4  方案3 TLS-J方法法矩阵条件数

      Figure 4.  Condition Numbers of Coefficient Matrix by Using Scheme 3 of TLS-J

      图 5图 6为将方案3总体最小二乘联合平差方法与最小二乘联合平差方法的压力源参数结果进行正演得到的垂直位移与水平位移图。由图 5图 6可知,方案3结果较之最小二乘联合平差结果整体上较好地拟合于观测数据。在图 5中IS2点及图 6中的P03点,两种方法得到的拟合结果相对于观测值的拟合情况较差,这与文献[5]、文献[6]和文献[10]研究成果一致。以图 6中P03点为例,方案3总体最小二乘联合平差结果与最小二乘联合平差结果在拟合观测值方向上是一致的,但均与观测水平位移观测值方向相差较大。这是因为P03点的GPS水平位移观测存在异常[3, 10],可能受到背景场影响,计算时可视为粗差处理[10]

      图  5  方案3 TLS-J方法与LS-J方法正演垂直位移

      Figure 5.  Vertical Displacement of Inversion Results by TLS-J and LS-J of Scheme 3

      图  6  方案3 TLS-J方法与LS-J方法正演水平位移

      Figure 6.  Horizontal Displacement of Inversion Results by TLS-J and LS-J of Scheme 3

    • 本文主要介绍了垂直位移和水平位移的Mogi模型线性化方法,结合线性化后垂直位移和水平位移观测方程,推导了总体最小二乘方法中观测向量与系数矩阵协因数阵的计算公式;针对采用垂直位移与水平位移的火山Mogi模型反演压力源参数过程中出现的病态情况,给出了病态总体最小二乘联合平差处理方法,并应用于长白山天池火山Mogi模型反演实例。由结果分析可知,总体最小二乘联合平差方法由于顾及了系数矩阵的误差影响,且在求解中考虑了不同类观测数据参与联合平差的比重,较之最小二乘联合平差方法更为合理,因而参数估计结果更加可信。数值试验结果也表明了总体最小二乘联合平差方法的可行性。

参考文献 (26)

目录

    /

    返回文章
    返回