等价条件平差模型的方差-协方差分量最小二乘估计方法

刘志平，朱丹彤，余 航，张克非,2

1. 中国矿业大学环境与测绘学院, 江苏 徐州 221116； 2. 皇家墨尔本理工大学空间科学研究中心，澳大利亚 维多利亚州 墨尔本 3001

随着现代测量技术的发展，测量数据处理的对象已经由单一同类观测数据转变为多源、多类或同类多因素的异方差结构观测数据，在这种情况下随机模型的准确性对参数解估计的统计特性具有重要影响[1-2]。方差-协方差分量估计能够极大地改善随机模型不精确的问题，从而成为国内外学者关注的焦点，并形成了一整套的理论框架[3-8]。近年来主要理论研究包括：①以平差结果之一的残差向量为基本输入量的残差型VCE方法，如Helmert法[1-2,7]、最小二乘方差分量估计法(简称LS-VCE)[9-11]、MINQUE法[12]、BIQUE法[13]、MLE-VCE法[14]、等效残差法[15-16]及EIV模型的VCE方法[4,8]等；②以等价条件闭合差为基本输入量的解析型VCE方法，如VCE-ECM法[17-18]。近年来的主要应用研究有：文献[19—20]研究了常规导线控制网中角度观测值与距离观测值的定权问题；文献[21—22]探讨了GNSS观测值随机模型的精化问题；文献[23—24]研究了坐标转换中不同公共点坐标的精度问题；文献[11,25—26]开展了GNSS站坐标时间序列的噪声估计及评价研究。

大量研究表明，LS-VCE方法[9-10,15]可导出Helmert、MINQUE、MLE-VCE等方法的方差-协方差分量估计公式，具有最优无偏特性。与此同时，基于等价条件平差模型的VCE-ECM法[17-18]以等价条件闭合差作为基本输入量，摒弃了平差值求解或残差估计过程。鉴于此，本文从方差-协方差分量最小二乘估计准则和等价条件平差模型出发，利用等价条件闭合差的二次型构造方差分量估计公式，进而将方差分量估计方程变换为线性Gauss-Markov形式并建立了基本估计方程；然后，顾及矩阵拉直、半拉直算子和Kronecker积的性质进行化简，导出了基于等价条件闭合差的方差-协方差分量最小二乘估计公式(least-square variance-covariance estimation based on the equivalent condition misclosure)，简称为LSV-ECM法。在此基础上，证明了本文LSV-ECM法与残差型VCE方法的等价性，定量分析了LSV-ECM法、Helmert法和LS-VCE法的计算复杂度，以边角网平差和中国区域GNSS站坐标时序建模实例验证本文方法的正确性和高效性。

1 等价条件闭合差的方差分量最小二乘估计

1.1 等价条件平差模型

设概括平差模型为

(1)

式中，A、[BTCT]均为行满秩系数矩阵；W为具有参数的条件方程闭合差；Z为限制条件方程闭合差；DL为观测值L的方差阵；V、x分别为待求的残差与参数向量；下标c、s分别为条件方程个数和限制条件方程个数；n、u分别为观测数和参数个数；平差模型自由度r=(c+s)-u。

(2)

(3)

1.2 方差-协方差分量最小二乘估计方法

假设概括平差模型式(1)中随机模型DL可表示为

(4)

式中，∑(·)表示累加运算；k为方差-协方差分量的个数；Qi为第i个协因数分量矩阵。

(5)

(6)

(7)

(8)

(9)

式中，D为复制矩阵(duplication matrix)[9]，也即半拉直算子与拉直算子的映射矩阵

(10)

分析可知，式(8)—式(10)含有大量的高维矩阵迭代运算(r2阶)，因而需进一步简化矩阵运算。考虑到矩阵拉直算子vec(·)、Kronecker积⊗和矩阵迹算子tr(·)具有如下性质[9]

vec(T)T·(U⊗S)·vec(VT)=tr(UVST)

(11)

(12)

(13)

参量模型条件平差具有参数的条件平差间接平差附有限制的间接平差B=0,C=0C=0A=-I,C=0A=-IQiAQiATHAQiATHTHQiHTHcQiHTcWWHWHWHcW+HsZDWADLATHADLATHTHDLHTHcDLHTc

现有文献已证明残差型VCE方法之间的等价性，并形成了一套较为完整的理论框架。通过概括平差因子矩阵R也可以证明本文LSV-ECM方法与残差型VCE方法的等价性。对于概括平差模型，LSV-ECM法与Helmert型通用VCE方法[7](简称通用Helmert法)等价；对于间接平差模型，LSV-ECM法与LS-VCE法[9,11]、基于等效残差的LS-VCE法[15]等价；对于具有参数的条件平差模型，LSV-ECM法与文献[10]相应算法等价。为节省篇幅，这里仅给出LSV-ECM法与通用Helmert法[7]等价关系的证明过程。

(14)

(15)

(16)

分析对比式(15)—式(16)和式(12)—式(13)可知，LSV-ECM方法可以导出通用Helmert法[7]，等价性得证。其次，残差型VCE方法(如Helmert法、LS-VCE法等)均以残差向量为输入量，而本文LSV-ECM法以等价条件闭合差为输入量，后者实现了平差值求解(残差为平差结果之一)与随机模型估计的分离。再次，等价条件闭合差维数低于残差向量维数，且前者为不需求解的已知量，因此计算效率方面明显优于前者。

表2 不同VCE方法的计算复杂度

Tab.2 A comparison of computational complexities of different VCE methods

2 应用结果及分析

2.1 边角网平差

方案1：分别采用通用Helmert法[7]、LS-VCE法[9,11]、本文LSV-ECM法进行方差分量估计和参数平差计算。其中测角和测边单位权中误差的迭代初值均取先验中误差，迭代终止条件为测角和测边的单位权中误差相等。

方案2：在方案1的基础上，剔除含有粗差的12号角度观测值和16号边长观测值[18]，基于余下的16个观测值进行方差分量估计和参数平差计算。

按上述方案分别进行计算，统计不同方法迭代收敛时的方差分量估值、参数估值及其中误差。另外，为避免单次计算的随机误差，将上述方案重复计算100次，并统计LS-VCE法和本文LSV-ECM法相对于通用Helmert法的耗时比值Tratio，所得结果如表3所示。

表3 不同VCE方法的估计结果

从方案1方差分量和待估参数的估计结果看，3种结果完全相同，表明3种方法均能获得一致的估计结果。同时比较3种方法的计算效率，通用Helmert法和LS-VCE的计算效率处于同一水平，本文所提出的LSV-ECM法的计算效率最高，其计算时间约为通用Helmert法的65%，较通用Helmert法和LS-VCE法均有明显提升，从而验证本文方法的有效性。

方案2中同样可以验证本文所提方法的正确性和有效性，不再赘述。同时，方案2在剔除两个粗差观测值以后，两类观测值的单位权方差、参数的估计精度均有较大提升，表明3种方差分量估计方法都无法抑制粗差影响，必须引入质量控制方法进行粗差识别与剔除。

2.2 GNSS站坐标时序建模

为进一步检验本文LSV-ECM法的正确性和高效性，选用中国大陆构造环境监测网络(crustal movement observation network of China，CMONOC)中16个长期线性趋势稳定的GNSS基准站2005.0014(DOY)～2015.0014(DOY)共10年的原始坐标时间序列进行噪声特性分析(方差-协方差分量估计)和三维运动速度估计(平差参数计算)，测站位置分布如图1所示。

GNSS站坐标时序的观测方程和随机模型分别为[11]

y(ti)=a+bti+ccos(2πti)+dsin(2πti)+

ecos(4πti)+fsin(4πti)+vi

(17)

(18)

式中，ti是以年为单位的时间序列历元点；a为常数项；b为线性速度项；c、d组合表示全年性周期运动；e、f组合表示半年性周期运动；vi为残差；σWN、σFN为所求的噪声分量大小；QWN、QFN分别为白噪声和闪烁噪声的协因数阵，具体形式见文献[11]。

图1 CMONOC所选测站分布Fig.1 Geographical distribution of selected stations in the CMONOC network

表4 噪声分量(WN,FN)估计结果

注：高斯白噪声(WN)的单位为mm；闪烁噪声(FN)的单位为mm/a0.25

对比表4中通用Helmert法、LS-VCE法和LSV-ECM法的计算结果可知，3种方法所得的噪声分量结果完全相同，验证了本文方法的正确性。同时，对比3种方法的计算效率可知，通用Helmert法与LS-VCE法基本一致，而本文LSV-ECM法有较大提高，三维方向计算时间比分别为通用Helmert法的74.0%、73.8%、72.1%，表明本文方法的高效性。

由表4噪声分量的估计结果可知，白噪声分量均远小于闪烁噪声分量，表明有色噪声为中国区域GNSS站坐标时间序列的主要噪声，在参数估计时若直接采用白噪声模型会导致估计结果有偏，并会产生较高的虚假估计精度。另外，对比不同方向的噪声分量估计结果，北方向和东方向的噪声分量相差不大，且量级较小，94%和88%的测站白噪声分量在1 mm以内，50%的测站闪烁噪声分量在3 mm/a0.25以内，而竖直方向的噪声分量远高于水平方向，这与现有结论相一致，约有56%的测站白噪声分量在3 mm以内，31%的测站闪烁噪声分量在9 mm/a0.25以内。

为进一步分析噪声大小与经纬度之间的关系，绘制了白噪声分量和闪烁噪声分量随经纬度的变化情况。由于噪声大小与纬度相关性较强、与经度相关性较弱，为节省篇幅，此处仅讨论噪声大小与纬度变化的关系，见图2。从图2可以看出，水平方向的白噪声和闪烁噪声分量大小较为平稳，竖直方向的噪声大小波动较大。从整体趋势上看，白噪声和闪烁噪声大小随纬度变化的趋势较为明显。具体表现为：噪声大小随纬度增加而逐渐减小，减少过程中有轻微波动，且在中纬度地区有翘尾现象。该现象原因分析可能是GNSS的GDOP值随纬度增大而降低[28]，从而使噪声分量表现出随纬度增加而减小的现象，而且噪声分量在中纬度地区呈现一定的翘尾现象。

此外，为检验方差-协方差分量估计结果的正确性和有效性，以CMONOC的公布数据(http:∥www.cgps.ac.cn/cgs/index.action)作为对比参考值，表5统计了中国区域GNSS站的三维运动速度和不确定度。由该表可知，中国区域不同站点水平方向均呈现东南方向运动趋势，这与亚欧板块的整体运动趋势相一致，但竖直方向运动变化差异性较大，其中约40%的测站呈现下沉趋势。对比本文结果与CMONOC的公布结果，除JIXN、HLAR站N方向，KMIN、GUAN站E方向，CHUN、BJFS站U方向以外，其余计算结果均在2倍中误差范围内。此外，参照文献[29]的计算结果，对比发现除CHUN站U方向外，其余测站的计算结果与本文结果具有较好的相符性，进一步验证了本文方法的正确性。需要指出的是，CHUN站U方向的差异性可能是数据预处理和共模误差提取策略不同所致。

图2 噪声分量大小与纬度的关系Fig.2 The relationship between noise components and latitude

3 结 论

本文基于等价条件平差模型和最小二乘准则，利用等价条件闭合差的二次型构建方差-协方差分量估计方程，并通过矩阵半拉直算子将其变换为线性Gauss-Markov形式，进而顾及矩阵拉直算子、半拉直算子和Kronecker积运算性质，导出了基于等价条件平差模型的方差-协方差分量最小二乘估计公式，简称LSV-ECM法。该方法实现了平差值求解(残差为平差结果之一)与随机模型估计的分离，有效兼顾了等价条件闭合差(已知)和最小二乘的特性。边角网平差实例的结果表明，本文的LSV-ECM法与通用Helmert法、LS-VCE法的结果完全相同，但计算效率更高，验证了该方法与残差型VCE方法的等价性和计算高效性。

分析指出了方差-协方差分量估计方法相较于常规GNSS-MLE法进行GNSS站坐标时间序列噪声分析的优势，并利用LSV-ECM法、通用Helmert法、LS-VCE法计算分析了中国区域16个GNSS站坐标时间序列在WN+FN模型下的噪声估计和站点速度。估计结果表明，有色噪声是中国GNSS站坐标时序的主要噪声，且白噪声和有色噪声分量随纬度增大而减小，速度估计结果与陆态网络的公布结果基本一致，水平方向的整体运动趋势呈东南方向，而竖直方向的运动趋势差异性较大。因此，计算结果也进一步验证了本文LSV-ECM法的正确性和高效性。

表5 站点运动速度和不确定度