中国处于欧亚板块、印度板块和太平洋板块的交汇处，3大板块的共同作用使大陆构造活动十分复杂。中国不是一个完整的大板块，各小板块之间相互作用，地球动力学机理十分复杂。垂直运动是地面点在铅垂线方向的移动量，国内外专家利用大地测量技术获取可靠的地壳垂直运动量化信息，科学合理地解释地球动力学机理问题。随着世界各地频繁出现不同程度的地面沉降灾害，大地测量技术也成为地面沉降灾害监测的重要手段。

精密水准测量是一种经典的垂直运动位移获取技术，自20世纪初起就应用于垂直形变监测^[1-5]。为了获取精确的垂直形变信息，陶本藻^[6]讨论了动态监测网在任何参考基准下参数有约束的参数估计普遍公式；欧吉坤等^[7]利用新丰江水库监测网进一步证明拟稳平差的必要性和优越性；黄立人等^[8]讨论了采用全球导航卫星系统（Global Navigation Satellite System，GNSS）速度作为水准点的速度基准，解决水准点速度缺失问题；王庆良等^[9]采用高差变化速率作为观测值估计水准点速度，避免估计水准点高程；李明峰等^[10]和曾安敏等^[11]将异常误差作为参数，采用函数模型补偿方法估计其大小；龚率等^[12]根据高差观测误差与地形起伏关系确定权值，利用赫尔默特方法估计调整各类权比。以上研究成果优化了水准网平差的函数模型与随机模型，提高了水准网变形分析的精度和可靠性。

地面沉降灾害是一种缓慢性、持久性灾害，具有形成缓慢、持续时间长、影响范围广、成因机制复杂和防治难度大^[13]。山东省是一个沉降比较严重的区域，德州市城区、滨州城区、东营城区、济宁城区、菏泽城区等形成了深层地下水漏斗中心区，成为山东省地面沉降重点防治区^[14]。山东省先后于2006年和2015年开展了全省似大地水准面精化工作，重复观测了山东省二等水准网，这些重复观测资料为研究山东省垂直运动速度提供了重要基础数据。鉴于以上研究基础与数据资源，本文提出一种多期水准网自适应动态平差方法，改进等价权函数，有效抑制异常误差影响，提高多期水准网动态平差精度，利用最小二乘配置方法获得山东省最新基于水准数据的垂直运动速度模型。

1 水准网动态平差

多期水准网动态平差可以采用整体动态平差方法，也可以采用序贯平差或Kalman滤波方法，获取水准点的平均速度或不同历元的速度。水准网动态平差过程中，一般将水准点高程与速度作为待估参数，采用水准高差观测值作为观测量，根据高差观测值的测站数定义权矩阵元素，采用经典间接平差、拟稳平差^[7]等方法获得水准点的垂直运动信息。

水准高差值误差方程为：

$ {V_h}_{_{ij}} = {\rm{ }} - \delta {H_i} - (t - {t_0}){\rm{ }}\delta {v_i} + {\rm{ }}\delta {H_j} + (t - {t_0}){\rm{ }}\delta {v_j} + {l_h}_{_{ij}} $

(1)

式中，H_i、H_j分别为i、j点的正常高；δH_i、δH_j为i、j点的正常高改正数；t为当前时间；t₀为参考历元时间；δv_i、δv_j分别为i、j点的正常高方向垂直运动速度改正数；

$ {l_h}_{_{ij}} = H_j^0 - H_i^0 - {h_{ij}} + (t - {t_0})v_j^0 - (t - {t_0})v_i^0 $

(2)

其中，h_ij为i点至j点的高差观测值；$H_j^0$为j点的正常高近似值；$H_i^0$为i点的正常高近似值；$v_i^0$、$v_j^0$分别为i点、j点的正常高方向垂直运动速度近似值。

误差方程（1）的矩阵表达形式为：

$ V = {\rm{ }}\left[ {\begin{array}{*{20}{c}} 1& \cdots &{ - 1}& \cdots &{(t - {t_0})}& \cdots &{ - (t - {t_0})}&0\\ \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots \\ 0&{\mathop 1\limits_i }& \cdots &{\mathop { - 1}\limits_j }& \cdots &{(t - {t_0})}& \cdots &{ - (t - {t_0})}\\ \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots \\ 0&{ - 1}& \cdots &1& \cdots &{ - (t - {t_0})}& \cdots &{(t - {t_0})} \end{array}} \right]{\rm{ }}\left[ {\begin{array}{*{20}{c}} {\delta {H_1}{\rm{ }}}\\ \vdots \\ {\delta {H_n}{\rm{ }}}\\ {\delta {v_1}}\\ \vdots \\ {{\rm{ }}\delta {v_1}} \end{array}} \right] - \left[ {\begin{array}{*{20}{c}} {\delta {l_{{h_{ik}}}}}\\ \vdots \\ {\delta {l_{h{i_j}}}}\\ \vdots \\ {{\delta _{{l_{hk}}}}} \end{array}} \right]{\rm{ }} = AX - L $

(3)

式中，A为系数矩阵；X为待估参数，包括正常高及对应速度；L为自由项。

2 水准网自适应动态平差

为了有效抑制变形监测网的异常误差影响，学者们分别从函数模型与随机模型方面提出异常误差有效抑制方法。曾安敏等^[11]分析了监测网的函数模型异常补偿方法，通过估计异常误差的大小，使异常观测值回到正常观测值，从而削弱异常对估计参数的影响；杨元喜等^[15]推导了秩亏平差模型的抗差解和拟稳抗差解，提高水准网平差精度。本文将在固定青岛水准原点的基础上，讨论多期水准网自适应动态平差方法，改进等价权函数形式，等价权函数采用平方形式进一步减小异常观测值的权值，进而削弱异常误差影响。

水准网自适应动态平差的参数估计值为^[15-17]：

$ X = {\rm{ }}{({A^T}\bar P{\rm{ }}A)^{ - 1}}{A^T}\bar P{\rm{ }}L $

(4)

式中，

$ {{\bar P}_{i,i}} = \left\{ {\begin{array}{*{20}{l}} {{\kern 1pt} {\kern 1pt} {P_{i,i}} \times 1.0,{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \left| {\frac{{{V_i}}}{{{\sigma _0}}}} \right| \le {k_0}}\\ {{P_{i,i}} \times \frac{{{k_0}{\sigma _0}}}{{{V_i}}},{\kern 1pt} {\kern 1pt} {\kern 1pt} \left| {\frac{{{V_i}}}{{{\sigma _0}}}} \right|\begin{array}{*{20}{l}} > \end{array}{k_0}} \end{array}} \right. $

(5)

其中，k₀取值范围1.0~1.5；σ₀为残差中位数；V_i为第i个观测值的残差。

为了进一步抑制异常误差影响，等价权函数采用平方形式构造异常观测值的权值，降低异常观测值对估计参数的贡献，即：

$ {{\bar P}_{i,i}} = \left\{ {\begin{array}{*{20}{l}} {{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {P_{i,i}} \times 1.0,{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \left| {\frac{{{V_i}}}{{{\sigma _0}}}} \right| \le {k_0}}\\ {{P_{i,i}} \times {{\left( {\frac{{{k_0}{\sigma _0}}}{{{V_i}}}} \right)}^2},\left| {\frac{{{V_i}}}{{{\sigma _0}}}} \right|\begin{array}{*{20}{l}} > \end{array}{k_0}} \end{array}} \right. $

(6)

3 基于最小二乘配置的垂直运动模型建立

为了进一步描述某区域整体垂直运动趋势，需要在水准网动态平差整体估计水准网点的高程与垂直运动速度的基础上，将离散水准网点垂直运动速度的格网化，从而建立垂直运动模型。常用的垂直运动模型建立方法有Kring、多项式、多面函数、最小二乘配置等。其中，最小二乘配置方法能够削弱垂直运动中的信号影响，提高垂直运动模型的可靠性。

最小二乘配置方法将垂直运动速度分为系统性部分和随机性影响部分，其观测模型为^[18]：

$ v = {A_1}x + {B_1}S + \Delta $

(7)

式中，v为水准点垂直运动速度向量；x为倾向性参数向量，维数为t_n×1；S为信号向量，维数为t_s×1；Δ为观测误差向量，维数为n × 1；A₁、B₁表示系统性参数以及信号参数对应的系数矩阵，维数分别为n×t_n、n×t_s；

$ {({A_1})_i} = \left[ {\begin{array}{*{20}{c}} 1&{{\lambda _i}}&{{\theta _i}}&{{\lambda _i}{\theta _i}}&{\lambda _i^2}&{\theta _i^2}& \cdots \end{array}} \right] $

(8)

$ {({B_1})_i} = \left[ {\begin{array}{*{20}{c}} 0&0&0&0&{\mathop 1\limits_i }&0& \cdots \end{array}} \right] $

(9)

其中，（θ_i，λ_i）为水准点的坐标。

最小二乘配置估计值与信号估计值分别为^[16]：

$ \left\{ \begin{array}{l} \hat x{\rm{ }} = {\rm{ }}{(A_1^T{\rm{ }}\Sigma _L^{ - 1}{A_1})^{ - 1}}A_1^T\Sigma _L^{ - 1}v\\ \hat S{\rm{ }} = {\rm{ }}{\Sigma _S}B_1^T\Sigma _L^{ - 1}\left( {v - {A_1}\hat x} \right)\\ \hat S\prime {\rm{ }} = {\Sigma _{S\prime S}}\Sigma _S^{ - 1}\hat S \end{array} \right. $

(10)

式中，Σ_S'S为已知点与未知点的互协方差矩阵；Σ_S为信号的协方差矩阵；

$ {\Sigma _L} = {B_1}{\Sigma _S}B_1^T + {\Sigma _\Delta } $

(11)

其中，Σ_Δ为观测噪声的协方差矩阵；其他参数含义同前。

4 算例分析

算例采用2006年和2015年山东省似大地水准面精化水准观测数据。2006年华东华中似大水准面精化项目开展了二等水准测量，其路线长度为8 260 km；2015年山东测绘基准体系优化升级项目同样也开展了二等水准测量，其路线长度为12 196 km；两期共有重合点491个（包括青岛水准原点）；两期共有高差观测值1 069个，其中第一期有515个高差观测值，第二期有554个高差观测值。图 1为水准重合点分布图。为了整体反映该区域垂直运动速度，算例首先进行水准网动态平差，获取各水准点垂直运动速度，然后利用最小二乘配置方法实现垂直运动速度模型的格网化。

4.1 水准网动态平差

水准网动态平差基准采用单点基准，起算点为青岛水准原点^[5]，其垂直运动速度为2.2 mm/a。两期二等水准网数据采用4种动态平差方案估计水准点高程与垂直运动速度，这4种方案为：①经典水准网动态平差方法（方案1）；②基于函数模型补偿的水准网动态平差方法^[11]（方案2）；③自适应水准网动态平差方法^[16]（方案3）；④改进的自适应水准网动态平差方法（方案4）。观测残差曲线绘于图 2，观测残差统计于表 1。

由表 1可以看出：

1）4种方案都能获得了水准重合点的高程与速度。由于区域地面沉降不均匀，导致高差观测中存在异常误差，将影响水准点高程与速度估计结果。与方案1相比，方案2~4精度都有提高，方案2采用函数模型异常误差补偿法，估计了异常误差大小，并在相应观测值中扣除异常误差，削弱了异常误差的影响；方案3、方案4采用抗差估计方法降低异常观测值的权比，调整异常数据对估计结果的贡献，提高估计精度；方案4是平差精度最高的方案，该方案进一步修正了等价权函数，有效控制异常误差影响。

2）方案3、方案4采用抗差估计能够抑制异常误差影响，提高水准网平差精度，但某区域的严重不均匀沉降时，可能导致二等水准网观测数据存在异常误差，而水准观测值也满足二等水准网测量规范要求，若删除某些异常观测值，可能会造成部分点无法获取到估计结果，因此，平差时只能降权而不能删除；方案4中采用平方形式进一步缩小异常观测值的权值，有效抑制异常误差影响。

3）方案2与方案3、方案4采用不同方式抑制异常误差影响，都能提高待估参数的可靠性。方案2采用函数模型补偿方法，量化异常误差大小，修正异常观测值到正常观测值，从而提高水准网平差精度；方案3、方案4采用随机模型改进方法抑制异常误差影响，这两种方案通过调整异常观测值的权比，降低了异常观测值对估计参数的影响。

4.2 垂直运动模型建立

利用水准网动态平差获得的491个水准重合点垂直运动速度，采用最小二乘配置方法建立2006—2015年山东省区域垂直运动速度模型，其垂直运动趋势如图 3所示。

2006年、2015年两期二等水准观测数据是山东省最新两期水准观测数据，是研究山东省区域垂直运动研究的重要基础数据资源。利用近十年两期二等水准网观测数据获得了山东省最新垂直运动速度模型。山东省北部、西部垂直运动呈下降趋势，东部、南部呈缓慢上升趋势。东北部、东南部垂直运动呈缓慢的上升趋势，速度在5 mm/a左右；在德州区域垂直运动存在一个持续下降趋势，由北向南速度在—30~—10 mm/a；在菏泽地区垂直运动呈现下降趋势，其平均速度—10 mm/a，且存在一个小沉降漏斗，最大速度近—40 mm/a；黄河三角洲地区垂直运动也呈现下降趋势，且在广饶区域存在一个沉降中心，年变化量达70 mm左右。

5 结语

水准测量是垂直运动研究的经典技术，在各行业中广泛应用。随着山东省多期水准网观测数据的积累，这为区域垂直运动研究提供了重要的基础数据资源。本文利用抗差估计理论抑制了异常误差影响，提高了各水准网点垂直运动速度精度，采用最小二乘配置山东全省范围最新垂直运动速度模型。该模型反映出山东省北部、西部整体上呈下降趋势，且在广饶、德州、菏泽都存在一个沉降漏斗，其中广饶沉降中心的垂直运动速度最大，数值达到—70 mm/a；东部、南部呈缓慢上升趋势，最大上升速度约5 mm/a。由于山东连续运行基准站也运行了10余年，积累的GNSS观测数据也能够反映山东省区域垂直运动趋势，因此，融合GNSS数据和水准数据开展区域垂直运动速度模型建立将是下一步的重点研究工作。

致谢: 感谢山东国土测绘院与陆态网络中心提供的试验观测数据。