地壳垂直运动是地壳运动的一种重要表现形式，长期观测成果可直接服务环境预测、防灾减灾，将极大推动大陆动力学、区域应力场以及地壳运动等方面的研究，为国家与地区的基础设施建设提供垂直运动信息。大地测量是定量描述地壳垂直运动的重要方法之一，在国内外得到广泛应用。

文献[1-5]利用全球导航卫星系统（global navigation satellite system，GNSS）、水准等技术获取垂直形变信息，分析了区域地壳垂直运动趋势。GNSS、水准观测数据源于不同观测仪器、误差模型、处理方法、参考基准，这导致GNSS、水准的垂直形变结果可能自身存在误差，且相互之间存在系统偏差，进而影响垂直形变分析的精度。另外水准网一般沿公路、铁路等路线布设，水准网中间存在较大空白区，导致其垂直形变量的可靠性不高。而GNSS控制网布设一般考虑网型结构，按照一定密度布设控制点，这些GNSS控制点能够进一步填补水准网的空白，有效提升垂直形变信息的可靠性。随着多期GNSS、水准观测数据的积累，GNSS水准点具有同时期的GNSS与水准观测数据，这有利于分析大地高变化与正常高变化之间的系统偏差，为实现两类高程变化融合统一提供了基础数据。为了充分发挥不同类型数据的优势，杨元喜等^[6-7]深入研究了大地测量数据融合理论与方法，提出了函数模型、随机模型方面的融合方法，充分发挥大地测量数据对计算结果的贡献。对于GNSS、水准数据融合，文献[8-14]研究了联合平差的速率基准、观测值选取、权重调整、函数模型选取等问题，充分发挥不同类观测值对垂直运动速度的贡献，提高垂直运动速度的估计精度。GNSS与水准数据可以通过联合平差方法实现融合，也可以使用GNSS水准点多期成果，利用拟合逼近方法实现GNSS与水准数据融合。因此，基于拟合逼近原理，本文提出一种自适应最小二乘配置的GNSS与水准融合方法，该方法以水准点获取的正常高变化速度为基础，利用GNSS水准点建立大地高变化速度与正常高变化速度的转换关系，进一步修正GNSS点的大地高变化速度后，再根据融合后的高程变化速度结果，采用反距离加权法建立垂直运动格网模型。

1 自适应最小二乘配置的GNSS水准数据融合

精密水准测量是垂直运动监测的经典技术手段，国内外学者通常将GNSS、星载合成孔径雷达干涉测量（interferometric synthetic aperture radar，InSAR）垂直运动结果与水准垂直运动结果相比较，检验垂直运动成果的可靠性。若建立控制点与GNSS水准点两类高程变化速度之间差异的函数关系，就能实现GNSS与水准两类垂直运动速度的融合。对于函数关系建模，可以采用多面函数^[15]、最小二乘配置^[16]、神经网络^[17]等拟合逼近方法。多面函数能够有效逼近曲面，无法有效抑制系统误差影响。神经网络理论能够削弱不确定因素等影响，但容易陷入极小值，且泛化能力差。最小二乘配置能够较好地解决两类高程变化量之间差异中含有倾向性、随机性因素的问题，控制系统误差等因素影响，提高两类高程变化量之间差异的拟合效果，但协方差矩阵一般根据经验确定。自适应最小二乘配置^{[16, 18]}利用自适应因子调整观测噪声协方差矩阵与信号协方差矩阵的比例，使其更加接近真实关系。因此，将水准正常高变化量作为基础，不考虑垂线偏差造成的大地高与正常高变化速度差异，利用自适应最小二乘配置方法建立GNSS水准点的大地高与正常高速度之间的转换关系，进一步实现非GNSS水准点的垂直运动速度融合。

最小二乘配置的大地高与正常高的速度差异观测量分为固定参数效应和随机效应两部分，前者表示倾向性部分或系统性部分，后者为随机性影响部分。其公式可表示为：

$\mathit{\boldsymbol{L}} = \mathit{\boldsymbol{AX}} + \mathit{\boldsymbol{BS}} + \mathit{\boldsymbol{\varepsilon }}$

(1)

式中，L为大地高变化速度与正常高变化速度之间的差异；A、B分别表示系统性参数以及信号参数对应的系数矩阵，维数分别为n × t_n、n × t_s, n为观测值个数，t_n为系统性参数个数，t_s为信号参数个数；X为倾向性参数向量，维数为t_n × 1；S为信号向量，维数为t_s × 1。

最小二乘配置估计值$\mathit{\boldsymbol{\widehat X}}$、信号估计值$\mathit{\boldsymbol{\hat S}}$与未知点信号估计值$ \mathit{\boldsymbol{\hat S'}}$分别为：

$ \left\{ {\begin{array}{*{20}{l}} {\mathit{\boldsymbol{\widehat X}} = {{\left( {{\mathit{\boldsymbol{A}}^{\rm{T}}}\mathit{\boldsymbol{ \boldsymbol{\varSigma} }}_\mathit{\boldsymbol{L}}^{ - 1}\mathit{\boldsymbol{A}}} \right)}^{ - 1}}{\mathit{\boldsymbol{A}}^{\rm{T}}}\mathit{\boldsymbol{ \boldsymbol{\varSigma} }}_\mathit{\boldsymbol{L}}^{ - 1}\mathit{\boldsymbol{L}}}\\ {\mathit{\boldsymbol{\widehat S}} = {\mathit{\boldsymbol{ \boldsymbol{\varSigma} }}_\mathit{\boldsymbol{S}}}{\mathit{\boldsymbol{B}}^{\rm{T}}}\mathit{\boldsymbol{ \boldsymbol{\varSigma} }}_\mathit{\boldsymbol{L}}^{ - 1}(\mathit{\boldsymbol{L}} - \mathit{\boldsymbol{A\hat X}})}\\ {{{\mathit{\boldsymbol{\widehat S}}}^\prime } = {\mathit{\boldsymbol{ \boldsymbol{\varSigma} }}_{\mathit{\boldsymbol{S'S}}}}\mathit{\boldsymbol{ \boldsymbol{\varSigma} }}_\mathit{\boldsymbol{S}}^{ - 1}\mathit{\boldsymbol{\hat S}}} \end{array}} \right. $

(2)

式中，Σ_L为观测值协方差矩阵；Σ_S为信号的协方差矩阵；Σ_S′S为已知点与未知点的互协方差矩阵。

$ {\mathit{\boldsymbol{ \boldsymbol{\varSigma} }}_\mathit{\boldsymbol{L}}} = \mathit{\boldsymbol{B}}{\mathit{\boldsymbol{ \boldsymbol{\varSigma} }}_\mathit{\boldsymbol{S}}}{\mathit{\boldsymbol{B}}^{\rm{T}}} + {\mathit{\boldsymbol{ \boldsymbol{\varSigma} }}_\mathit{\boldsymbol{\varepsilon }} } $

(3)

式中，Σ_ε为观测噪声的协方差矩阵。

自适应最小二乘配置的自适应因子α^{[16, 18-19]}为：

$ \alpha=\frac{\hat{\sigma}_{0 S}^{2}(i)}{\hat{\sigma}_{0 \varepsilon}^{2}(i)} $

(4)

式中，$ \hat{\sigma}_{0 \varepsilon}^{2}(i) $为观测噪声的第i次验后单位权方差；$ \hat{\sigma}_{0 S}^{2}(i) $为信号的第i次验后单位权方差。自适应因子α具体计算方法可以参考文献[16, 18-19]。

观测噪声协方差矩阵$\mathit{\boldsymbol{ \boldsymbol{\widehat \varSigma} }}_\varepsilon ^{i + 1} $为：

$ \mathit{\boldsymbol{ \boldsymbol{\widehat \varSigma} }}_\varepsilon ^{i + 1} = \frac{1}{\alpha }\mathit{\boldsymbol{ \boldsymbol{\varSigma} }}_\varepsilon ^i $

(5)

如此反复迭代，直至α趋于1。

2 垂直运动模型的反距离加权建立方法

目前，我国局部区域出现了严重的非均匀沉降，这也给精确描绘区域垂直运动规律增加了难度。若采用某曲面函数整体拟合大区域垂直运动，可能造成垂直运动剧烈区域存在较大偏差。因此，为了更加可靠地反映垂直运动信息，本文采用水准点以及融合后的GNSS点垂直运动速度作为拟合量，充分利用垂直运动剧烈特征点变化信息，使用反距离加权法实现格网化，建立真实可靠的高分辨率垂直运动格网模型。

反距离加权法^[20-21]是一种整体计算方法，根据待定格网点与控制点之间距离确定权重。该方法依据各点与待定点的距离平方大小确定权值，具体公式如下：

$ P_{k}=\left\{\begin{array}{ll} {\frac{1}{R^{2}}, } & {S_{k} \leqslant R} \\ {\frac{1}{S_{k}^{2}}, } & {S_{k}>R} \end{array}\right. $

(6)

式中，R为搜索半径；S_k为第k点与待定点之间的距离。

格网点垂直运动速度V_i，j为：

$ {V_{i, j}} = \frac{{\sum\limits_{k = 1}^N {{P_k}} V_k^\prime }}{{\sum\limits_{k = 1}^N {{P_k}} }} $

(7)

式中，N为控制点总个数；V_k′为第k点垂直运动速度；i、j分别为格网点的行、列号。

3 算例分析

算例采用2006、2015年山东省似大地水准面精化GNSS和水准观测数据。水准网观测按照二等水准网观测纲要分别于2004—2006年、2015年进行观测，GNSS控制网按照C级网观测纲要也于2004—2006年、2015年同步进行数据采集。水准点速度是两期似大地水准面精化水准网成果之差除以时段长度，两期水准网平差都采用青岛原点为基准点，且假设青岛原点垂直运动速度为0，山东省区域两期水准网重合点共有490个。两期GNSS控制网重合153个GNSS点，且都为GNSS水准点，水准点与GNSS点分布如图 1所示。为了实现GNSS水准高程变化速度有效融合，采用如下步骤进行处理：（1）使用GAMIT/GLOBK软件处理，选取相同的起算基准点，获得GNSS点在国际地球参考框架ITRF2008下不同历元的成果，计算GNSS控制点的垂直运动速度。（2）采用水准网动态平差软件获取水准网重合点垂直运动速度。（3）选取大地高变化与正常高变化一致的GNSS水准点，建立正常高变化速度与大地高变化速度的函数关系，实现正常高变化速度与大地高变化速度的融合。（4）利用融合后的高程变化速度建立垂直运动模型。

3.1 GNSS水准数据融合

为了验证算法可靠性，算例从153个GNSS点中选取60个点作为融合点（见图 2），选取原则如下：（1）大地高变化速度与正常高变化速度方向尽量一致。（2）大地高变化速度与正常高变化速度之间差异的绝对值不大于1.5 mm/a。对于GNSS、水准高程变化融合，从60个融合点中，选取40个点作为拟合点，20个点作为检核点，将大地高变化速度与正常高变化速度差异作为拟合量，采用3种方案进行计算，方案1：多面函数；方案2：最小二乘配置；方案3：自适应最小二乘配置。

拟合点与检核点误差统计见表 1，拟合点残差见图 3，检核点残差见图 4。

由表 1、图 3、图 4的计算结果可以看出，3种方案都能有效减小GNSS高程变化速度与水准高程变化速度之间的差异，实现GNSS高程变化速度与水准高程变化速度之间的有效融合。自适应最小二乘配置比多面函数方法拟合精度与检核精度有所提高。另外，在数据准备过程中，要考虑GNSS与水准的起算基准、有色噪声等因素影响。

为了建立高精度区域垂直运动模型，本文实验利用了60个融合点作为拟合点，采用自适应最小二乘配置建立融合关系，获取剩余93个GNSS点的垂直运动速度。

3.2 垂直运动格网模型建立

利用490个水准点及转换后的93个GNSS点，即583个高程变化点，采用反距离加权方法获取垂直运动格网模型，其中搜索半径R为1 km，山东省区域垂直运动趋势如图 5所示。

从图 5中可以看出：

1）近十年来，山东垂直运动总体上呈西降东升、北降南升的趋势。西北、西南部地区沉降严重，整体呈下降趋势，垂直运动速度平均20.0 mm/a；德州、菏泽、广饶地区出现3个沉降中心，其中广饶地区沉降最严重，十年平均沉降速度70.0 mm/a左右，这主要是由于地下水开采造成的；而在胶东半岛呈缓慢的上升趋势，垂直运动速度约5.0 mm/a。

2）采用反距离加权方法构建垂直运动速度模型，最大限度地发挥距离较近GNSS点与水准点对格网点的贡献，提高垂直运动模型的真实性与精度。

3）水准网一般沿公路布设，其包围区域内为空白区。分布均匀的GNSS网点能够有效补充水准网空白区，提高垂直运动建模所需的控制点分辨率，提升垂直运动模型的可靠性。

4 结语

GNSS与水准是地壳垂直运动速度研究的重要技术。利用自适应最小二乘配置方法建立GNSS水准点的大地高变化速度与正常高变化速度之间的融合模型，可实现GNSS控制点的大地高变化速度与正常高变化速度的统一，有效发挥GNSS高时间分辨率和水准的高精度优势。本文利用反距离加权法建立了山东省高精度、现势性强的垂直运动速度模型，可以为省内基础工程建设以及地面沉降灾害防治等提供决策依据。山东省垂直运动趋势整体上西部、北部下降，其平均下降速度达20.0 mm/a左右，另外在广饶、德州、菏泽区域存在几个沉降中心，特别是广饶地区沉降速度达到70.0 mm/a左右，而在东部、东南部垂直运动趋势为上升，这些不均匀地面沉降对该区域内的基础设施建设造成了严重威胁，大大降低了高程基准成果的可靠性和现势性。随着各省市连续运行参考站（continuously operating reference stations，CORS）系统的进一步完善，利用多期水准与连续GNSS观测数据，削弱各类数据之间的系统偏差，采用自适应融合方法实现多源数据有效融合，建立国家、区域地壳垂直运动速度模型，维护高程基准，是当前大地测量的重要研究方向。

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