-
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) 式中，Ux为x方向位移；Uy为y方向位移；(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) 式中, eUxi、eUyi分别为对应第i个地表水平方向位移Ux、Uy的误差。
分别对式(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)可知，垂直位移与水平位移观测值对应的系数矩阵B1与B2中不仅存在含有随机误差的坐标数据，同时还存在待估压力源参数；压力源参数的变化必将引起系数矩阵B1与B2的结构性扰动，同时线性化迭代解法如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, y1为n1×1维垂直位移观测值；el1为对应垂直位移观测方程的随机误差；l2=y2-y02, y2为n2×1维水平位移观测值；el2为对应水平位移观测方程的随机误差；B1为n1×4维垂直位移观测方程列满秩系数矩阵, EB1为其误差矩阵；B2为n2×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) 式中，eB1和eB2分别是将矩阵EB1和EB2按列拉直得到的列向量；σ12、σ22分别为垂直及水平观测值的方差分量；Ql1、QB1、Ql2和QB2分别为对应观测向量误差和系数矩阵误差的协因数阵；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-XTⓧIn1]$\left[ \begin{array}{l} {\mathit{\boldsymbol{e}}_{{\mathit{\boldsymbol{l}}_1}}}\\ {\mathit{\boldsymbol{e}}_{{\mathit{\boldsymbol{B}}_1}}} \end{array} \right]$, ē2=[In2-XTⓧIn2]$\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]$=CTC，C为上三角矩阵，则式(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]T，l=[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} $$ -
考虑到向量y01与y02是由近似压力源参数与监测点坐标数据计算得到的近似位移向量，因此，向量y01、y02的协因数阵存在[20]，分别记为Qy01和Qy02。由式(7)、式(9)和式(10)可知，在给定了先验压力源信息θ0后，由θ0求得的第i个监测点的位移fi(θ0)、φxi(θ0)与φyi(θ0)对xi与yi分别求偏导得：
$$ {\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)中，Qy01i、Qy02i分别为协因数阵Qy01和Qy02中对应第i个监测点的子矩阵，联合所有监测点的精度信息便可求得相应Qy01和Qy02，因此，Ql1=Qy1+Qy01、Ql2=Qy2+Qy02，其中，Qy1和Qy2对应水平位移和垂直位移观测值的协因数阵。
同理，对于确定系数矩阵的协因数阵QB1与QB2，由于式(8)与式(11)中B1与B2各元素同样为监测点坐标的函数，则系数矩阵中各元素需分别对xi与yi求偏导，并结合协因数传播律求得各系数阵的协因数阵[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}}}}$；②更新B1、B2、l1和l2；③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点坐标(X，Y)及水平东向、水平北向和垂直方向的形变场结果(VE，VN，VH)[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/km Y/km VH/(mm·a-1) X/km Y/km VE/(mm·a-1) VN/(mm·a-1) IS0 10.30 25.18 0 P0 0.91 6.34 -7.182 2 43.630 1 IS1 9.84 23.25 4.04 P01 10.15 20.35 7.287 0 9.830 8 IS2 9.84 20.25 8.07 P02 18.09 3.22 19.981 0 0.641 6 IS3 8.67 18.77 5.03 P03 3.64 -13.12 -10.951 8 -9.526 7 IS4 7.57 17.07 6.19 P04 -4.13 -5.56 -26.314 7 -20.486 5 IS5 6.44 15.70 9.12 P05 -16.11 1.22 -20.964 7 3.627 0 IS6 5.15 13.66 12.67 P06 -14.94 16.23 -13.269 6 11.158 1 IS7 3.98 12.29 15.01 P07 -12.54 24.46 -2.481 7 1.651 2 IS8 2.80 10.70 19.75 IS9 1.55 9.37 22.71 IS10 1.06 7.97 30.85 IS11 0.89 6.13 47.03 IS12 0.51 4.45 65.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-J 2.385 28 0.388 22 6.879 19 174.186 71 文献[5] -0.862 50 -6.807 00 6.096 90 155.800 00 文献[6] 2.500 00 0.500 00 6.900 00 107.474 37 文献[10] 2.257 40 1.356 00 5.606 00 145.166 15 方案1 2.926 96 -0.812 17 6.859 66 189.693 76 0.99 方案2 2.387 33 0.388 47 6.879 66 174.201 39 0.50 方案3 2.497 07 0.129 21 7.035 25 177.737 39 0.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模型反演实例。由结果分析可知，总体最小二乘联合平差方法由于顾及了系数矩阵的误差影响，且在求解中考虑了不同类观测数据参与联合平差的比重，较之最小二乘联合平差方法更为合理，因而参数估计结果更加可信。数值试验结果也表明了总体最小二乘联合平差方法的可行性。
Application of Total Least Squares Joint Adjustment to Volcano Inversion of Mogi Model
-
摘要: 针对垂直位移与水平位移的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|}$的函数最小化能获得合理的压力源参数估值结果和相对权比大小，具有一定的实际参考价值。Abstract: In this paper, a total least squares joint (TLS-J) adjustment method is proposed to the inversion of Mogi model with vertical and horizontal observational data. The proposed method considers the errors in both observation vector and coefficient matrix of the functional model of joint adjustment problem. Three forms of the minimum discriminate function methods are adopted to determine the weight scaling factor which are used to weigh the vertical and horizontal observation data. In view of the existing ill-posed problems in the joint adjustment, the L-curve method is adopted to determine the ridge parameter. Through practical examples, the total least squares joint method is systematically applied to the inversion of the Mogi model of Changbaishan Tianchi volcano. The results show that the discriminant function as the minimum can obtain the reasonable value of pressure source parameters and the relative weight ratio, which has a certain referential value to practical applications.
-
Key words:
- total least squares /
- joint adjustment /
- weight scaling factor /
- Mogi model /
- ill-posed problem
-
图 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 方案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
表 1 GPS点与水准点坐标及形变值[6]
Table 1. Coordinates of GPS and Leveling Points and Corresponding Displacement[6]
点名 水准点点位观测数据(第1类数据) 点名 GPS点位观测数据(第2类数据) X/km Y/km VH/(mm·a-1) X/km Y/km VE/(mm·a-1) VN/(mm·a-1) IS0 10.30 25.18 0 P0 0.91 6.34 -7.182 2 43.630 1 IS1 9.84 23.25 4.04 P01 10.15 20.35 7.287 0 9.830 8 IS2 9.84 20.25 8.07 P02 18.09 3.22 19.981 0 0.641 6 IS3 8.67 18.77 5.03 P03 3.64 -13.12 -10.951 8 -9.526 7 IS4 7.57 17.07 6.19 P04 -4.13 -5.56 -26.314 7 -20.486 5 IS5 6.44 15.70 9.12 P05 -16.11 1.22 -20.964 7 3.627 0 IS6 5.15 13.66 12.67 P06 -14.94 16.23 -13.269 6 11.158 1 IS7 3.98 12.29 15.01 P07 -12.54 24.46 -2.481 7 1.651 2 IS8 2.80 10.70 19.75 IS9 1.55 9.37 22.71 IS10 1.06 7.97 30.85 IS11 0.89 6.13 47.03 IS12 0.51 4.45 65.79 
下载: 导出CSV
表 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|}$的判别函数最小化法
下载: 导出CSV
表 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-J 2.385 28 0.388 22 6.879 19 174.186 71 文献[5] -0.862 50 -6.807 00 6.096 90 155.800 00 文献[6] 2.500 00 0.500 00 6.900 00 107.474 37 文献[10] 2.257 40 1.356 00 5.606 00 145.166 15 方案1 2.926 96 -0.812 17 6.859 66 189.693 76 0.99 方案2 2.387 33 0.388 47 6.879 66 174.201 39 0.50 方案3 2.497 07 0.129 21 7.035 25 177.737 39 0.88
下载: 导出CSV
-
[1] Mogi K. Relations Between the Eruption of Various Volcanoes and the Deformations of the Ground Surfaces Around Them[J]. Bulletin of the Earthquake Research Institute, University of Tokyo, 1958, 36:99-134 https://www.researchgate.net/publication/247382386_Relations_of_the_Eruptions_of_Various_Volcanoes_and_the_Deformation_of_Ground_Surfaces_around_them [2] 胡亚轩, 王庆良, 崔笃信, 等. Mogi模型在长白山天池火山区的应用[J].地震地质, 2007, 29(1):144-151 doi: 10.3969/j.issn.0253-4967.2007.01.013 Hu Yaxuan, Wang Qingliang, Cui Duxin, et al. Application of Mogi Model at Changbaishan Tianchi Volcano[J]. Seismology and Geology, 2007, 29(1):144-151 doi: 10.3969/j.issn.0253-4967.2007.01.013 [3] 胡亚轩, 王庆良, 崔笃信, 等.三种压力源模型对火山区地面变形的影响[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 Models in Volcanic Area[J]. Seismology Research of Northeast China, 2005, 21(3):33-38 doi: 10.3969/j.issn.1674-8565.2005.03.003 [4] 李克, 刘俊清, 盘晓东, 等. 2000-2007年期间长白山天池火山区地壳变形监测与分析[J].地震地质, 2009, 31(4):639-646 doi: 10.3969/j.issn.0253-4967.2009.04.007 Li Ke, Liu Junqing, Pan Xiaodong, et al. Crust Deformation Monitoring and Research in Tianchi Volcano Area, Changbai Mountains from 2000-2007[J]. Seismology and Geology, 2009, 31(4):639-646 doi: 10.3969/j.issn.0253-4967.2009.04.007 [5] 胡亚轩, 郝明, 王雄, 等. L曲线法在反演火山区压力源参数中的应用[J].大地测量与地球动力学, 2009, 29(2):66-70 http://d.old.wanfangdata.com.cn/Periodical/dkxbydz200902015 Hu Yaxuan, Hao Ming, Wang Xiong, et al. Application of L-curve to the Inversion of Pressure Source Parameters in Volcano Area[J]. Journal of Geodesy and Geodynamics, 2009, 29(2):66-70 http://d.old.wanfangdata.com.cn/Periodical/dkxbydz200902015 [6] 陈国浒. 长白山天池火山形变监测与模拟研究[D]. 北京: 中国地震局地质研究所, 2007 http://www.wanfangdata.com.cn/details/detail.do?_type=degree&id=Y1172335 Chen Guohu. Deformation Monitoring and Simulation Research of Mt. Changbai Tianchi Volcano[D]. Beijing: Institute of Geology, China Earthquake Administration, 2007 http://www.wanfangdata.com.cn/details/detail.do?_type=degree&id=Y1172335 [7] 胡亚轩, 施行觉, 王庆良, 等.腾冲火山区地表垂直形变分析[J].大地测量与地球动力学, 2003, 23(2):37-41 http://d.old.wanfangdata.com.cn/Periodical/dkxbydz200302006 Hu Yaxuan, Shi Xingjue, Wang Qingliang, et al. Analysis of Vertical Deformation in Tengchong Volcano Area[J]. Journal of Geodesy and Geodyna-mics, 2003, 23(2):37-41 http://d.old.wanfangdata.com.cn/Periodical/dkxbydz200302006 [8] 施行觉, 胡亚轩, 毛竹, 等.以垂直形变资料反演腾冲火山区岩浆活动性的初步研究[J].地震研究, 2005, 28(3):256-261 doi: 10.3969/j.issn.1000-0666.2005.03.011 Shi Xingjue, Hu Yaxuan, Mao Zhu, et al. Vertical Deformation Inversion of Magma Activity in Tengchong Volcanic Area[J]. Journal of Seismological Research, 2005, 28(3):256-261 doi: 10.3969/j.issn.1000-0666.2005.03.011 [9] Bifulco I, Raiconi G, Scarpa R. Computer Algebra Software for Least Squares and Total Least Norm Inversion of Geophysical Models[J]. Computers & Geosciences, 2009, 35(7):1427-1438 http://www.wanfangdata.com.cn/details/detail.do?_type=perio&id=JJ0211153394 [10] 王乐洋.基于总体最小二乘的大地测量反演理论及应用研究[J].测绘学报, 2012, 41(4):629 http://www.wanfangdata.com.cn/details/detail.do?_type=perio&id=QK201203765934 Wang Leyang. Research on Theory and Application of Total Least Squares in Geodetic Inversion[J]. Acta Geodaetica et Cartographica Sinica, 2012, 41(4):629 http://www.wanfangdata.com.cn/details/detail.do?_type=perio&id=QK201203765934 [11] 于冬冬. 病态总体最小二乘解算方法及应用[D]. 南昌: 东华理工大学, 2015 http://cdmd.cnki.com.cn/Article/CDMD-10405-1015990441.htm Yu Dongdong. Research on the Ill-Posed Total Least Squares Algorithm and Its Application[D]. Nanchang: East China University of Technology, 2015 http://cdmd.cnki.com.cn/Article/CDMD-10405-1015990441.htm [12] 姚宜斌, 孔建.顾及设计矩阵随机误差的最小二乘组合新解法[J].武汉大学学报·信息科学版, 2014, 39(9):1028-1032 http://ch.whu.edu.cn/CN/Y2014/V39/I9/1028 Yao Yibin, Kong Jian. A New Combined LS Method Considering Random Errors of Design Matrix[J]. Geomatics and Information Science of Wuhan University, 2014, 39(9):1028-1032 http://ch.whu.edu.cn/CN/Y2014/V39/I9/1028 [13] Schaffrin B, Wieser A. On Weighted Total Least Squares Adjustment for Linear Regression[J]. Journal of Geodesy, 2008, 82(7):415-421 doi: 10.1007/s00190-007-0190-9 [14] 王乐洋, 许才军, 鲁铁定.病态加权总体最小二乘平差的岭估计解法[J].武汉大学学报·信息科学版, 2010, 35(11):1346-1350 http://ch.whu.edu.cn/CN/Y2010/V35/I11/1346 Wang Leyang, Xu Caijun, Lu Tieding. Ridge Estimation Method in Ill-Posed Weighted Total Least Squares Adjustment[J]. Geomatics and Information Science of Wuhan University, 2010, 35(11):1346-1350 http://ch.whu.edu.cn/CN/Y2010/V35/I11/1346 [15] 王乐洋, 许才军.附有相对权比的总体最小二乘平差[J].武汉大学学报·信息科学版, 2011, 36(8):887-890 http://ch.whu.edu.cn/CN/Y2011/V36/I8/887 Wang Leyang, Xu Caijun. Total Least Squares Adjustment with Weight Scaling Factor[J]. Geomatics and Information Science of Wuhan University, 2011, 36(8):887-890 http://ch.whu.edu.cn/CN/Y2011/V36/I8/887 [16] 王乐洋, 余航, 陈晓勇. Partial EIV模型的解法[J].测绘学报, 2016, 45(1):22-29 http://d.old.wanfangdata.com.cn/Periodical/chxb201601004 Wang Leyang, Yu Hang, Chen Xiaoyong. An Algorithm for Partial EIV Model[J]. Acta Geodaetica et Cartographica Sinica, 2016, 45(1):22-29 http://d.old.wanfangdata.com.cn/Periodical/chxb201601004 [17] Wang Leyang, Xu Guangyu. Variance Component Estimation for Partial Errors-in-Variables Models[J]. Studia Geophysica et Geodaetica, 2016, 60(1):35-55 doi: 10.1007/s11200-014-0975-2 [18] 王乐洋, 许才军, 温扬茂.利用STLN和InSAR数据反演2008年青海大柴旦Mw6.3级地震断层参数[J].测绘学报, 2013, 42(2):168-176 http://d.old.wanfangdata.com.cn/Periodical/chxb201302002 Wang Leyang, Xu Caijun, Wen Yangmao. Fault Parameters of 2008 Qinghai Dacaidan Mw6.3 Earthquake from STLN Inversion and InSAR Data[J]. Acta Geodaetica et Cartographica Sinica, 2013, 42(2):168-176 http://d.old.wanfangdata.com.cn/Periodical/chxb201302002 [19] 胡川, 陈义.非线性整体最小二乘平差迭代算法[J].测绘学报, 2014, 43(7):668-674 http://www.cnki.com.cn/Article/CJFDTOTAL-CHXB201407003.htm Hu Chuan, Chen Yi. An Iterative Algorithm for Nonlinear Total Least Squares Adjustment[J]. Acta Geodaetica et Cartographica Sinica, 2014, 43(7):668-674 http://www.cnki.com.cn/Article/CJFDTOTAL-CHXB201407003.htm [20] Jazaeri S, Amiri-Simkooei A R. Weighted Total Least Squares for Solving Non-linear Problem:GNSS Point Positioning[J]. Survey Review, 2015, 47(343):265-271 doi: 10.1179/1752270614Y.0000000132 [21] Tikhonov A N. Solution of Incorrectly Formulated Problems and the Regularization Method[J]. Soviet Mathematics Doklady, 1963, 4:1305-1308 http://ci.nii.ac.jp/naid/10004315593 [22] 余航. 总体最小二乘联合平差方法及其应用研究[D]. 南昌: 东华理工大学, 2016 http://cdmd.cnki.com.cn/Article/CDMD-10405-1016758498.htm Yu Hang. Research on the Total Least Squares Joint Adjustment and Its Application[D]. Nanchang: East China University of Technology, 2016 http://cdmd.cnki.com.cn/Article/CDMD-10405-1016758498.htm [23] 王乐洋, 余航.总体最小二乘联合平差[J].武汉大学学报·信息科学版, 2016, 41(12):1683-1689 http://ch.whu.edu.cn/CN/Y2016/V41/I12/1683 Wang Leyang, Yu Hang. Total Least Squares Joint Adjustment[J]. Geomatics and Information Science of Wuhan University, 2016, 41(12):1683-1689 http://ch.whu.edu.cn/CN/Y2016/V41/I12/1683 [24] 王振杰. 大地测量中不适定问题的正则化解法研究[D]. 武汉: 中国科学院测量与地球物理研究所, 2003 http://www.wanfangdata.com.cn/details/detail.do?_type=degree&id=Y568497 Wang Zhenjie. Research on the Regularization Solution of Ill-Posed Problems in Geodesy[D]. Wuhan: Institute of Geodesy and Geophysics, Chinese Aca-demy of Sciences, 2003 http://www.wanfangdata.com.cn/details/detail.do?_type=degree&id=Y568497 [25] 王彬, 高井祥, 刘超, 等. L曲线法在等价权抗差岭估计模型中的应用[J].大地测量与地球动力学, 2012, 32(3):97-101 http://d.old.wanfangdata.com.cn/Periodical/dkxbydz201203021 Wang Bin, Gao Jingxiang, Liu Chao, et al. Application of L-curve Method to Equivalent Weighted Robust Ridge Estimation Model[J]. Journal of Geodesy and Geodynamics, 2012, 32(3):97-101 http://d.old.wanfangdata.com.cn/Periodical/dkxbydz201203021 [26] Dieterich J H, Decker R W. Finite Element Mode-ling of Surface Deformation Associated with Volca-nism[J]. Journal of Geophysical Research, 1975, 80(29):4094-4102 doi: 10.1029/JB080i029p04094 -
点击查看大图
图(6) / 表(3)
计量
- 文章访问数: 1818
- HTML全文浏览量: 112
- PDF下载量: 241
- 被引次数: 0
