0 引言

Mogi模型(Mogi，1958)是日本学者Kiyoo Mogi依据Yamakawa(1955)理论公式，将岩浆压力源与火山地表形变间建立联系的一种模型，该模型在模拟和解释大量的火山区形变有重要的应用意义(Hu et al., 2003; 胡亚轩等, 2004, 2007a, 2009; 王乐洋, 2012; 余航, 2016)；在长白山天池火山及腾冲火山的研究中，Mogi模型得到了广泛的应用(施行觉等, 2005; 陈国浒, 2007).获得精确的地表形变数据对研究火山活动及火山喷发的预测有重大的意义.目前，随着观测手段技术的提高，GPS全天候的观测及InSAR技术的发展，使得对数据的获取更加轻松.在数据解算方面，最小二乘是Mogi模型压力源参数反演的主要方法(Hu et al., 2003; 胡亚轩等, 2004, 2007b, 2009; 施行觉等, 2005; 陈国浒, 2007)；然而，Mogi模型的系数矩阵含有存在误差的坐标数据及含有待求的压力源参数，其理论上属于非线性变量含误差模型(王乐洋，2012)；目前，总体最小二乘(Total Least Squares, TLS)是变量误差模型(Errors-in-Variables, EIV)常用的解算方法，该方法在国内外已经得到了广泛的研究(Schaffrin and Wieser, 2008; Fang, 2011; 王乐洋, 2012; 王乐洋和许才军, 2013; 曾文宪, 2013)，将算法应用于实践是当前要做的工作.Bifulco等(2009)采用结构总体最小二乘的L1范数，结合1982—1985年的垂直形变数据计算得到了意大利南部Campi Flegrei火山压力源参数；王乐洋(2012)针对垂直观测数据及水平观测数据，使用总体最小二乘联合平差方法反演了长白山天池火山压力源参数；于冬冬(2015)研究了Mogi模型的病态总体最小二乘方法，探讨了虚拟观测法、谱修正迭代法和共轭梯度法在Mogi模型反演中的应用；余航(2016)分析了两类数据在平差中占的比重关系，采用附有权比因子的加权总体最小二乘联合平差方法来反演长白山火山天池火山压力源参数，该方法只考虑两类观测，在病态问题的处理上没有顾及正则化部分的权比关系.

Mogi模型垂直位移单一反演、水平位移单一反演或垂直位移与水平位移联合反演法方程都呈现出严重病态问题，病态问题(袁振超等, 2009; 葛旭明和伍吉仓, 2012, 2013; 朱建军等, 2012; 顾勇为等, 2016)的处理一直是研究的热点.在实际应用中经常遇到这一问题，如航空重力向下延拓、GPS快速定位、大地测量反演；处理病态问题的方法也较多，如正则化方法、岭估计法、截断奇异值法等，处理病态问题的关键是正则化参数的选取，且需要根据实际情况选择合适的方法.朱建军等(2012)针对准则带参数的情况提出了虚拟观测法，虚拟观测法是将附有先验信息的参数作为观测方程与实际观测方程进行联合平差，此时可以使用方差分量估计方法确定病态问题的正则化参数；冯光财等(2007)使用方差分量估计方法确定病态问题的正则化参数并取得了较好的结果；王乐洋和于冬冬(2014)推导了病态总体最小二乘的虚拟观测法，但是基于单位权阵推导的公式，不具有统一性；Koch和Kusche(2002)在使用卫星数据确定重力场时，考虑了卫星与梯度计测量两种类型的数据，将未知参数的先验信息作为观测方程进行联合解算，使用方差因子法确定正则化参数并取得了较好的结果，该方法使用的是G-M(Gauss-Markov)模型.方差分量估计方法(赵俊和郭建锋, 2013; Amiri-Simkooei et al., 2009; Amiri-Simkooei, 2013; Mahboub, 2014)是针对验前平差随机模型不准确而提出的验后估计方法，同时，该方法也可用来确定不同类数据之间的权比，许才军等(2016)和Fan等(2017)分析了在地震反演的方差分量方法并取得了较好的结果.在地球科学反演问题中，联合反演是一个具有很强优势的方法，将各资料联合解算可以进行优势互补、加快数据收敛及得到更稳定的结果(Ji et al., 2004)；在地震破裂过程反演中，远场的体波数据有良好的时间分辨率，因为其较高的数据采样率，但是空间分辨率不足；GPS位移数据和InSAR数据是近场地表形变数据，反演得到的结果空间分辨率高；将大地测量数据作为地震波数据的约束条件，进行联合反演可以得到更佳的结果(Zhang et al., 2012; Li et al., 2016).在本文中，GPS数据因为其垂直方向分辨率不高，因此增加水准测量的垂直方向数据，增强数据完整性；而约束条件在本文中也作为观测方程与观测数据进行联合解算，可以减弱病态性的影响，增强参数的稳定性.

综上所述，本文考虑火山Mogi模型反演中水平位移与垂直位移反演时系数矩阵含有误差的情况，并计算观测值与系数矩阵中随机元素的协方差矩阵，以加权形式对公式进行推导，最后获得更为一般的参数估值表达式；在水平位移与垂直位移的联合反演中，考虑反演时的病态性，尝试使用虚拟观测法来减弱模型的病态性，区别于王乐洋和余航(2018)的算法，本文是进行三类数据的联合平差，并尝试通过方差分量估计方法平衡三类数据的权比因子，收敛速度较快，推导了更一般的表达式及迭代算法.

1 Mogi模型 1.1 垂直位移与水平位移的Mogi模型

由膨胀压力源引起地表垂直位移和水平位移与等效压力源参数之间的表达式为(Yamakawa, 1955)：

(1)

(2)

式中，Δh表示为地表垂直位移，Δr表示为地表水平位移，K表示为体积弹性模量，D表示为源的中心深度，μ表示为剪切模量，r表示为地表径向距离，当r=0时垂直位移取得最大值.

Hu等(2003)对Mogi模型进行简化处理，当地壳为泊松介质时，取泊松比υ=0.25，此时体积模量K与剪切模量μ有关系式为(余航, 2016)：

(3)

因此式(1)和式(2)垂直位移和水平位移的表达式可以简化为

(4)

(5)

若将岩浆房的膨胀量化为半径为R的等效球体(王乐洋, 2012; 余航, 2016)：

(6)

式(4)和(5)等价于：

(7)

(8)

胡亚轩等(2007b)在水平位移反演中考虑到水平位移是东向分量和北向分量的形式，将水平位移改写为

(9)

(10)

由式(7)可知，根据地表垂直位移Δh反演压力源的待反演参数为：源的深度D、等效半径R及源在平面上投影的坐标(x₀, y₀).第i个监测点到源之间的地表径向距离为r²=(x_i-x₀)²+(y_i-y₀)².若有n个监测点，则有：

(11)

式中，i=1, 2, …, n，θ=[R D x₀ y₀]^T为待估参数.

对式(11)在参数初值θ₀处展开并舍去二次项得(王乐洋, 2012; 余航, 2016)：

(12)

式中，, 且中各元素可表示为

(13)

式(11)可线性化为

(14)

写成矩阵形式为

(15)

式(15)即为垂直位移反演压力源参数的平差模型.针对水平位移式(9)和(10)，若有n个监测点，则有：

(16)

将式(16)在参数初值θ₀处展开得到：

(17)

式中，，且中各元素为

(18)

中各元素为

(19)

n个监测点不同水平位移分量可线性化为

(20)

写成矩阵形式为

(21)

式(21)即为水平位移反演压力源参数的平差模型.

1.2 协因数阵计算

由式(15)和式(21)的平差模型可知，垂直位移和水平位移观测值对应的系数矩阵A₁和A₂中不仅存在含有随机误差的坐标数据，还存在待估压力源参数θ；压力源参数θ的变化必然引起系数矩阵的结构性扰动，在线性化后必然使得求导得到的系数矩阵含有误差.此时既考虑观测向量误差又考虑系数矩阵的误差平差方法即为总体最小二乘，相应的针对火山联合反演的即为非线性总体最小二乘联合平差问题(王乐洋, 2012; 余航, 2016)，记向量y₀₁、y₀₂的协因数阵为Q_y01、Q_y02，由θ₀计算求得的第i个监测点的位移f_i(θ₀)、φ_xi(θ₀)和φ_yi(θ₀)可改写为

(22)

由f_i(θ₀)对x_i和y_i求偏导得：

(23)

φ_xi(θ₀)与φ_yi(θ₀)分别对x_i和y_i求偏导得：

(24)

若第i个监测点坐标的平面精度(σ_xi, σ_yi)已知，则可以通过协方差传播律求得对应观测值的协因数阵：

(25)

式中，Q_y01ⁱ、Q_y02ⁱ分别表示为协因数阵Q_y01和Q_y02中对应的第i个子矩阵.相应的，针对系数矩阵的协因数阵，可参考余航(2016).

2 Mogi模型反演的虚拟观测法

将线性化后的垂直位移模型及水平位移分别写成Partial EIV(Partial Errors-in-Variables)模型(Xu et al., 2012)形式：垂直位移函数模型为

(26)

垂直位移随机模型为

(27)

水平位移函数模型为

(28)

水平位移随机模型为

(29)

考虑到联合解算中的病态性，根据朱建军等(2012)，虚拟观测方程及随机模型为

(30)

(31)

在本文中，将水平位移与垂直位移联合得到随机模型为

(32)

令：

(33)

则可以建立平差准则为

(34)

式中，

令：

(35)

则可以得到：

(36)

Koch和Kusche(2002)在解算GPS数据时考虑两类观测，并增加了虚拟观测方程，使用方差因子法确定正则化参数；在本文中，考虑火山Mogi模型联合反演的病态性，使用虚拟观测法，则随机模型变为

(37)

并令，则构造平差准则为

(38)

式中，

对平差准则进行求偏导计算并令其等于0得到：

(39)

进而可得到：

(40)

将式(39)第四式代入第二式得：

(41)

式中，从而可以得到参数估值的表达式为

(42)

3 Mogi模型反演的方差分量估计

从式(42)可以看出，最终参数估值的表达式为三类数据联合的标准化最小二乘形式，因此，从式(42)可以得到：

(43)

即：

(44)

此时等同于最小二乘间接平差，其随机模型可以表示为

(45)

式中，Q₁⁰、Q₂⁰、Q_y3⁰表示为初始的协因数阵.

若三类数据的方差分量已知，且令Q₁⁰=P₁^-1，Q₂⁰=P₂^-1，Q_y3⁰=P₃^-1，则(42)式可以表示为

(46)

式中，表示为权比.

对式(42)进行Helmert方差分量估计，即存在：

(47)

(48)

(49)

在迭代过程中，假设第一类观测值的权阵合理，通过改变第二类和第三类数据的权来达到修正参数估值的效果，此时两个权比因子的形式为

(50)

式中，k表示方差分量估计迭代中收敛时的迭代次数.

使用方差分量估计方法反演Mogi模型的压力源参数的迭代流程见图 1.

4 算例分析 4.1 火山Mogi模型模拟算例

由Mogi模型正演得到地表垂直位移，[x, y]∈[-20, 20] km×[-20, 20] km，相邻垂直形变观测点间隔为2 km，在该区域中正演得到地表水平位移，得到的图见图 2；根据式(5)，考虑地壳为泊松介质且取泊松比为0.25，设置压力源参数真值为：中心坐标为x₀=0 km，y₀=0 km，源深D为6 km，等效半径R为150 m；在实际观测中水准点或GPS监测点并不密集，因此，在该区域中随机选择选点作为监测点，对垂直位移模拟均值为0，标准差为1 mm的随机误差，对水平位移模拟均值为0，标准差为2 mm的随机误差.

现分别使用垂直位移的总体最小二乘单一反演(TLS-V)、水平位移的总体最小二乘单一反演(TLS-H)、总体最小二乘联合反演(TLS-J)和本文方法进行解算，得到的压力源参数结果及与真值的差值范数‖Δβ‖见表 1.

从表 1可以看出，总体最小二乘联合反演可以得到比总体最小二乘单一反演更佳的参数估值，进一步说明了王乐洋(2012)方法的有效性，而本章方法将改善病态问题的虚拟观测方程作为一类数据进行方差分量估计，一方面赋予了正则化参数实际意义，另一方面进行方差分量估计可以得到更合理的参数估值；此时进行方差分量估计则变为三类数据的联合平差，通过方差分量估计方法来确定不同类数据的权比，从而达到更合理的参数估值.

4.2 长白山天池火山Mogi模型反演

长白山天池火山区域位于中国东北大陆东缘、中—朝古老克拉通北缘附件和西太平洋板块俯冲带前沿，且同时位于北东、北西及东西走向的断层的交汇部位.长白山天池火山区域的垂直与水平位移观测资料来源于陈国浒(2007)，见表 2.长白山火山水准测量按一等水准测量规范执行，共有13个水准点组成，2002—2005年进行了四期水准测量，每公里偶然中误差ΔM=0.13 mm，成果质量可信；GPS观测网由8个站点组成，2002年和2003年的两期观测，因同时利用了8台接收机进行同步观测，且观测时间段长度大于48 h，结果可信度较高.本文主要采用这些数据进行反演火山压力源参数.相关数据也在王乐洋(2012)和余航(2016)中得到过研究.

针对表 2所提供的数据，建立相应的模型及求得协因数阵，以最小二乘联合反演结果作为初值，根据本文推导的方法进行解算，得到的结果见表 3.

使用Helmert方差分量估计确定正则化参数时，三个方差分量估值随迭代次数变化如表 4，变化图见图 3.从表 4可以看出，方差分量估计经过10次迭代即可收敛，三个单位权方差因子估值最后结果都为1，而权比因子ω和μ的结果分别为0.007和3.5386×10^-10，此时μ的值即为正则化方法中的正则化参数.正则化参数量级达到了10^-10，这与Koch和Kusche(2002)得到的结果是相似的，这是因数数据结构的原因.根据表 3得到的压力源参数估值结果正演得到的垂直与水平位移见表 5.

根据表 4和图 3可知，方差分量估计方法确定正则化参数时只需要迭代10次，有效的减少了计算时间.虚拟观测法作为准则带参数的通用方法，在Mogi模型反演中，使用虚拟观测法来减弱联合反演中病态性对反演参数的影响；而针对病态问题中正则化参数的选择，本文将正则化参数反应在虚拟观测值的权上，并转换为数据间的权比选取，通过方差分量估计方法来确定不同类数据之间的权比关系，理论上更加严密.

从表 5可以看出，通过方差分量估计方法反演得到的压力源参数正演得到的垂直位移和水平位移与原始的垂直、水平位移观测值接近；同样，在GPS点P03上，正演得到的东方向与北方向的位移与原始观测值较大，这与胡亚轩等(2007b)和王乐洋(2012)指出的P03点的GPS水平位移观测存在异常是相同的，计算时可做粗差处理.王乐洋(2012)指出，总体最小二乘联合反演火山Mogi模型的结果要优于最小二乘联合反演，但经过线性化过程后，模型对初值比较敏感，算法还需要进一步研究；本文在总体最小二乘联合反演的基础上，将Mogi模型的垂直位移和水平位移作为两类观测值，考虑观测值先验权值的不准确及病态性，增加虚拟观测值作为第三类观测值，进行三类观测值的联合平差并进行方差分量估计，得到修正权比因子下改进的压力源参数估值.

5 结论

本文针对Mogi模型垂直位移与水平位移联合反演中的病态问题，使用准则带参数的通用方法-虚拟观测来处理该问题，在已有方法的基础上进行改进，将参数先验信息作为一类观测数据，与垂直位移和水平位移观测方程联合进行三类数据的联合平差，推导了与已有方法相似的参数估值表达式，得到更为一般的形式.当考虑三类数据的联合平差时，使用方差分量估计方法确定病态中的正则化参数(即数据间的权比).方差分量估计方法是通过修正观测值随机模型之间在平差时的占比关系从而来达到平衡参数估值的作用，且该方法收敛速度较快，方差分量估计只迭代了10次.由火山Mogi模型的模拟算例和长白山实际数据的计算结果可以看出，当模型病态时，使用虚拟观测法进行三类观测值的联合平差并采用方差分量估计获得不同类数据的之间的权比可以得到修改后的参数估值，在实际中值得借鉴.