1 引 言
近40年来，国内外学者持续作了大量的方差-协方差分量估计(VCE)理论与方法研究，至今仍然是多类观测值融合或多因素同类观测值平差[1, 2]以及大地测量反演[3, 4]的难点与开放课题。当前，VCE代表性成果主要包括：基于二次型期望公式导出的Helmert法[5, 6]；基于二次估计无偏性、不变性与最小范数准则导出的MINQUE法[7, 8]；基于二次估计无偏性、不变性和最小方差准则导出的BIQUE法[9]；基于二次估计无偏性、非负性和最小方差准则导出的BQUNE法及其改进[10, 11]；基于极大似然估计原理导出的MLE法[12, 13]；基于最小二乘准则的LS-VCE法[14]以及不同VCE方法之间等价性证明等研究[15]。此外，结合测量平差质量控制理论，发展了稳健BIQUE法[16]、含有粗差或系统误差的稳健VCE法[17]、处理病态模型的VCE法[18]、基于Lp估计的稳健VCE法[19]以及基于等效残差积的稳健VCE法[20]等成果。
现有VCE方法的共同特点是：以残差向量作为随机模型分量估计的基本输入量，故本文称之为残差型VCE方法。但是，残差向量是平差结果之一，使得残差型VCE方法必须迭代求解。因此，残差型VCE方法无法避免迭代算法的收敛性[21]与迭代停止条件选取的派生问题。目前，广泛采用的迭代停止条件是各类观测值单位权方差相等。然而，研究发现该停止条件的本质是方差一致性检验统计量χ2等于其数学期望[22]，即与模型自由度r相等。换言之，迭代处理后的χ2检验统计量不再具有随机性，而是在任何情况下均能通过平差模型检验的确定值。因此，残差型VCE方法改变了χ2检验统计量的统计特征。此外，残差型VCE方法大多数是基于间接平差模型或概括平差模型，而很少利用条件平差模型的特殊性。在条件平差模型中，条件闭合差虽是残差向量的线性函数，但为先验已知量。因此，若能以条件闭合差的某种表达作为方差分量估计的基本输入量，就可望导出VCE解析估计方法。与此同时，残差型VCE方法的固有缺陷及其派生问题便可迎刃而解。
2 VCE-ECM方法 2.1 等价条件平差模型设概括平差模型为[1]
式中，A为行满秩系数矩阵；[BT CT]为行满秩系数矩阵；W=－(AL+Bx0+A0)表示具有参数的条件方程闭合差；Z=－(Cx0+C0)表示限制条件方程闭合差；DL表示观测值L的方差阵。V、x分别为待求的残差与参数向量。下标c、s分别表示条件方程个数和限制条件方程个数；n、u分别表示观测数和参数个数；平差模型自由度r=c+s－u。
利用正交投影矩阵H消去参数向量，即
式中，H由零空间算子null(·)得到[23]
式(1)等式两边均左乘H，根据式(2)可将概括平差模型化为等价的条件平差模型
式中，A=HcA；=HcW+HsZ。其中，Hc由H前c列构成，Hs由H后s列形成。
基于最小二乘准则VTDL－1V=min，残差向量V可由W表示如下
式中，Na=ADLAT。
为叙述方便，文中将称为等价条件闭合差向量(equivalent condition misclosure,ECM)，其方差阵DW
平差模型的正确性一般采用方差一致性原理进行检验[1]。不论何种平差模型，若观测值方差阵正定，则由式(5)可导出以W表示的χ2检验统计量
式(7)表明等价条件闭合差可方便地用于χ2统计量计算。在一定的显著水平上，若χ2(r)统计量落在了接受域，则认为平差模型正确合理，否则认为平差模型误差(仅讨论随机模型DL)超出允许范围。由于参数向量x的LS估值与观测值方差阵DL可分别表示为观测值L的线性函数与二次型[1]，欲使观测数据处理更加完善，必须同时进行参数估计与方差估计(精度评定)。鉴此，下文讨论利用等价条件闭合差W进行方差分量估计的方法。
2.2 可逆方差分量模型对于k类观测值，其方差阵可用分块结构矩阵表示如下
式中，σ0ij2表示观测值的单位权方差-协方差因子；qij表示已知的观测值协因数阵。
目前，大多数研究均直接基于方差因子σ0ij2构造方差分量模型[1]
式中，qij称为协因数分量矩阵。qij的维数、矩阵分块结构与DL一致，且其中ij与ji位置分别为qij,qijT，其他位置为零矩阵。
为分析方差分量模型式(9)存在的不足，以两类观测值为例，将式(9)重新表达为
分析式(9)～(10)可发现，qij具有对称性和线性独立性[14]，qij(i=j)的行列式为零，矩阵不可逆。在协因数分量模型构造研究中，其对称性条件可以保证方差阵对称，其线性独立性条件可确保方差因子法矩阵可逆。然而，线性独立性条件并不是方差因子法矩阵可逆的必要条件，却导致了qij(i=j)的不可逆。鉴于此，为使构造的方差分量矩阵同时满足对称性和可逆性，对式(9)右端第一项进行可逆线性变换
式中
由于qij（i≠j）对称可逆，将式(11)第3项中的σ0ij2,qij(i≠j)分别按对应顺序编入第1项αi,Qi(k＜i≤m)，因此基于式(11)可得新的方差分量模型
式中，m=k(k+1)/2。
相应的，以两类观测值对式(12)说明如下
式中，α1=σ0112－σ0222/2,α2=σ0112+σ0222/2,α3=σ0122，且有σ0112=α1+α2,σ0222=α2－α1。
分析式(11)～(13)可知，对于i∈[1,k]，αi为方差因子σ0ij2或其可逆线性组合，Qi为协因数分量qij或其可逆线性组合。此外，由于Rankl=k，则αi可唯一地确定σ0ij2，Qi对称和可逆性质分别保证了方差阵对称性和方差因子法矩阵可逆性。因此，所建立的新模型称为可逆方差分量模型，αi、Qi可相应地称为可逆方差因子、可逆协因数分量。
2.3 基于的VCE解析估计首先，给出下文推导过程中用到的二次型期望公式[1]：若有服从任一分布的n维随机向量Y，其数学期望为μ，方差阵为Σ，则Y的任一二次型的数学期望可表达为
式中，Q为一个n阶对称可逆阵，tr(·)表示迹算子。
根据式(6)与式(12)，等价条件闭合差的可逆方差分量模型为
式中，Qi=AQiAT。
由于Qi对称可逆，则Qi对称可逆。此外，由E(V)=0，可得E(W)=E(AV)=0。因此，将式(14)中的Q替换为Qi、Y替换为等价条件闭合差W，可得
进一步地，将式(15)结果代入式(16)右端，并去除左端期望符号，整理可得
与式(17)同理，考虑m个可逆方差因子，可导出方差分量模型的解估计
分析上述推导过程可知，式(11)～(13)未对αi的取值范围作要求，表明新方法适用于各类随机模型估计研究。其次，式(18)由二次型期望公式导出，其输入量为已知的可逆协因数阵Qi和等价条件闭合差W，表明新方法具有无偏性且无需迭代求解，从而克服了残差型VCE方法的固有缺陷及派生问题。再次，只需变换Qi,W计算式，便可得到某特定平差模型的方差分量解析估计式(见表 1)。因此，将该方法称为基于等价条件闭合差W的VCE无偏解析方法，简称VCE-ECM法。
| 模型 | 条件平差 | 具参数的条件平差 | 间接平差 | 附限制的间接平差 |
| 参量 | B=0,C=0 | C=0 | A=－I,C=0 | A=－I |
| Qi | AQiAT | HAQiATHT | HQiHT | HcQiHTc |
| W | W | HW | HW | HcW+HsZ |
| χ2 | WT(ADLAT)－1W | W(HAQiDLATHT)－1WT | W(HDLHT)－1WT | W(HcDLHcT)－1WT |
算例1[1]：有一边角网，A、B、C是已知点，P1、P2是待定点，网中观测了12个角度，编号为1～12，观测了6条边长，编号为13～18。采用间接平差模型V=Bx－L，先验测角中误差为1.5″，测边中误差为2 cm。记后验测角方差为σβ2，测边方差为σs2，观测值均相互独立，由式(10)构造可逆方差分量模型。为验证所提出新方法的正确有效性，两个设计方案如下(计算结果见表 2)。
| 方法 | 残差型VCE法(Helmert方法不同迭代次数结果) | VCE-ECM法 | |||
| 方案A | 0 | 1 | 8 | - | |
| σβ2 | 2.250 | 3.581 | 3.638 | 3.840 | |
| σs2 | 4.000 | 6.168 | 5.929 | 5.230 | |
| χ2 | 22.080 | 13.998 | 13.999 | 13.971 | |
| x1、y1、σp1 | 1.588、-0.853、1.132 | 1.577、-0.864、1.422 | 1.558、-0.883、1.421 | 1.499、-0.945、1.418 | |
| x2、y2、σp2 | -5.518、12.506、2.008 | -5.563、12.461、2.508 | -5.639、12.385、2.483 | -5.872、12.138、2.401 | |
| 方案B | 0 | 1 | 4 | - | |
| σβ2 | 2.250 | 0.751 | 0.743 | 0.777 | |
| σs2 | 4.000 | 0.395 | 0.330 | 0.228 | |
| χ2 | 3.065 | 10.598 | 10.999 | 11.238 | |
| x1、y1、σp1 | 2.796、0.238、1.236 | 2.552、0.237、0.587 | 2.509、0.238、0.561 | 2.398、0.241、0.513 | |
| x2、y2、σp2 | -3.579、1.779、2.577 | -3.854、17.352、1.221 | -3.909、17.262、1.177 | -4.059、17.008、1.107 | |
方案A：分别采用Helmert法、VCE-ECM法进行方差分量估计与参数平差计算。Helmert法以先验方差为初始值，迭代停止条件为测角与测边单位权方差相等。
方案B：在方案A基础上，剔除含有粗差的12号角度观测值和15、16号边长观测值[19]，基于余下的15个观测值分别进行随机模型估计与参数平差。
在方案A中，比较第0次(先验值)与第2次、第4次Helmert法迭代结果，坐标改正参数平差值及中误差几乎没有变化，但观测值方差与χ2检验统计量却作了明显调整。当Helmert法经过8步迭代收敛时，χ2检验统计量等于14(平差模型自由度)，表明残差型VCE方法中的迭代停止条件的隐性约束为χ2=r。同时，比较两方法的坐标点位中误差可知，VCE-ECM方法略优于残差型VCE方法。由此可见，残差型VCE方法迭代停止条件与LS估计准则的内在要求不一致，改变了χ2检验统计量的统计特征。与此相反，VCE-ECM方法并不以χ2=r为隐性约束，而是在平差计算之前便获得测角、测边方差分量的无偏解析估计，后续平差结果也保留了χ2检验统计量的统计性质。
方案B同样可以得出，VCE-ECM方法完全克服了残差型VCE方法在计算效率与χ2统计量特征两方面的固有缺陷。比较两方案可知，在剔除了3个含有粗差的观测值后，方案B的参数精度较方案A提高了2～3倍。分析认为VCE可用观测值的二次型表示，使得VCE较参数估计对粗差更加不具稳健性，表明在实际应用中应引入必要的质量控制方法。
算例2：为了进一步验证VCE-ECM方法的正确性与有效性，设计仿真算例对VCE-ECM与MINQUE方法进行统计比较。试验模型为Vi=[I10⊗150×1]x－Li,(i=1,2)，其中，观测向量L1～N(0,1),L2～N(0,1.5),σL1,L2=0.5；I10表示10阶单位阵，150×1表示50维列向量。分别按VCE-ECM与MINQUE方法试验1000次，结果见表 3和图 1。
| 估计量 | σL12 | σL22 | σL1,L2 | tr(D ) | MSE||||2
| |
| 均值 | MINQUE均值 | 0.999 591 | 1.498 210 | 0.500 400 | 0.166 620 | 0.151 562 |
| ECM均值 | 0.999 072 | 1.498 937 | 0.500 505 | 0.166 597 | 0.151 556 | |
| MINQUE置信区间 | ±0.015 495 | ±0.024 559 | ±0.014 231 | ±0.002 765 | 0.105 841 | |
| ECM置信区间 | ±0.016 632 | ±0.024 968 | ±0.014 250 | ±0.002 741 | 0.105 740 |
|
| 图 1 不同方法的χ2估计结果比较 Fig. 1 Result comparisons of χ2 using different methods |
比较表 3中的σL12、σL22、σL1,L2可知，两方法所得方差分量估计的均值之差不超过10－3；在置信水平5%条件下，ECM方法较MINQUE方法的置信区间略宽(其单侧宽度之差不超过10－3)，表明ECM方法虽在理论上不具有最小方差性，但在统计意义下两方法所得方差-协方差分量估计值的差异甚微。
比较表 3中的tr(D)、MSE及其置信区间，ECM方法略优于MINQUE方法，进一步表明ECM方法与MINQUE方法几乎相当。结合图 1可以看出，对于MINQUE方法，χ2=r在1000次试验中恒成立；对于ECM方法，χ2值在以r值为中心的一定范围内波动。因此，仿真结果进一步表明，MINQUE迭代解法与不合理的隐性约束χ2=r等价[22]，改变了χ2检验统计量的随机特性；而ECM解析法为解决上述问题提供了一种有效的新思路。
以Helmert、MINQUE、BIQUE及MLE等为代表的方差-协方差分量估计已形成了一个普遍性理论框架，其方法特点是以残差向量为基本输入量的迭代解法。文中详细分析了方差-协方差迭代估计的固有缺陷及派生问题，引入零空间算子将概括平差模型转换为等价条件平差模型，并通过构造可逆方差分量模型提出了以等价条件闭合差(ECM)为基本输入量、以解析法为估计形式的方差-协方差分量估计新方法。计算结果与分析表明，VCE-ECM方法由于无需迭代求解而保留了χ2检验统计量的统计性质，且与现有VCE方法的方差-协方差分量估计值在统计意义上无明显差异。需指出的是，方差因子非负性约束、可逆方差分量模型构造方法对VCE-ECM方差分量估计的影响有待进一步研究。
| [1] | CUI Xizhang, YU Zongchou, TAO Benzao, et al. Generalized Surveying Adjustment (New Edition)[M]. Wuhan: Wuhan University Publishing House, 2005.(崔希璋, 於宗俦, 陶本藻, 等. 广义测量平差(新版)[M]. 武汉: 武汉大学出版社, 2005.) |
| [2] | AMIRI-SIMKOOEI A R.Least-squares Variance Component Estimation Theory and GPS Applications [D].Delft:Delft University of Technology, 2007. |
| [3] | DU Zhixing, OU Jikun, JIN Fengxiang, et al. Optimization Inversion of the Relative Weight Ratio λ in the Joint Inversion Models[J]. Acta Geodaetica et Cartographica Sinica, 2003,32(1):15-19.(独知行, 欧吉坤, 靳奉祥, 等. 联合反演模型中相对权比的优化反演[J]. 测绘学报, 2003,32(1):15-19.) |
| [4] | XU Caijun, WANG Leyang. Progress of Joint Inversion of Geodetic and Seismological Data for Seismic Source Rupture Process[J]. Geomatics and Information Science of Wuhan University, 2010,35(4):457-462.(许才军, 王乐洋. 大地测量和地震数据联合反演地震震源破裂过程研究进展[J]. 武汉大学学报:信息科学版, 2010,35(4):457-462.) |
| [5] | YU Zongchou. The General Formulas of Helmert Type for Estimation Variance and Covariance Components[J]. Journal of Wuhan Technical University of Surveying and Mapping, 1991,16(2):8-17.(於宗俦. Helmert型方差—协方差分量估计的通用公式[J]. 武汉测绘科技大学学报, 1991,16(2):8-17.) |
| [6] | LIU Changjian, CHAI Hongzhou, WU Hongju, et al. Extended Formulae of Helmert Type for Estimating Variance Components[J]. Acta Geodetica et Cartographica Sinica, 2008,37(1):1-4.(刘长建, 柴洪洲, 吴洪举, 等. 扩展的Helmert型方差分量估计公式[J]. 测绘学报, 2008,37(1):1-4.) |
| [7] | WANG Zhizhong, ZHU Jianjun. A Universal Formula of MINQUE of Variance Components[J]. Journal of Central South University of Technology, 2001,32(4):433-436.(王志忠, 朱建军. 方差分量的MINQUE通用公式[J]. 中南工业大学学报:自然科学版, 2001,32(4):433-436.) |
| [8] | PENG Junhuan, SHI Yun, LI Shuhui,et al. MINQUE of Variance-Covariance Components in Linear Gauss-Markov Models [J]. Journal of Surveying Engineering, 2011, 137(4):129-139. |
| [9] | CROCETTO N, GATTI M, RUSSO P. Simplified Formulae for the BIQUE Estimation of Variance Components in Disjunctive Observation Groups [J]. Journal of Geodesy, 2000, 74(6): 447-457. |
| [10] | SJBERG L E. Non-negative Variance Component Estimation in the Gauss-Helmert Adjustment Model[J]. Manuscripta Geodaetica, 1984(9): 247-280. |
| [11] | ESHAGH M, SJBERG L E. The Modified Best Quadratic Unbiased Non-Negative Estimator (MBQUNE) of Variance Components [J]. Studia Geophysica et Geodaetica, 2008, 52(3): 305-320. |
| [12] | YU Zongchou. A Universal Formula of Maximum Likelihood Estimation of Variance-covariance Components[J]. Acta Geodaetica et Cartographica Sinica, 1994,23(1):6-13.(於宗俦. 方差-协方差分量极大似然估计的通用公式[J]. 测绘学报, 1994,23(1):6-13.) |
| [13] | YU Zongchou. A Universal Formula of Maximum Likelihood Estimation of Variance-covariance Components [J]. Journal of Geodesy, 1996, 70(4): 233-240. |
| [14] | TEUNISSEN P J G, AMIRI-SIMKOOEI A R. Least-squares Variance Component Estimation [J]. Journal of Geodesy, 2008,82(2):65-82. |
| [15] | YU Zongchou. The Uniformity Theory of Estimation of Variance-covariance Components[J]. Acta Geodaetica et Cartographica Sinica, 1991,20(3):161-171.( 於宗俦. 方差—协方差分量估计的统一理论[J]. 测绘学报, 1991,20(3):161-171.) |
| [16] | WANG Xinzhou, LIU Jingnan, TAO Benzao. Robust Best Invariant Quadratic Unbiased Estimation[J]. Journal of Wuhan Technical University of Surveying and Mapping, 1995,20(3):240-245.(王新洲, 刘经南, 陶本藻. 稳健最优不变二次无偏估计[J]. 武汉测绘科技大学学报, 1995,20(3):240-245.) |
| [17] | LIU Changjian, MA Gaofeng, QIAO Shubo. Variance Component Estimation with Gross Errors or Systematic Errors[J]. Bulletin of Surveying and Mapping, 2005(10):7-8.(刘长建, 马高峰, 乔书波. 含有粗差或系统误差的方差分量估计[J]. 测绘通报, 2005(10):7-8.) |
| [18] | XU Peiliang, SHEN Yunzhong, FUKUDA Y, et al. Variance Component Estimation in Linear Inverse Ill-posed Models [J]. Journal of Geodesy, 2006,80(2):69-81. |
| [19] | ZHANG Yali, PENG Junhuan, ZHOU Jinguo, et al. The Estimate of Variance Components Based on M-Estimate Residuals[J]. Journal of Chongqing Jianzhu University, 2007,29(6):62-66.( 张亚利, 彭军还, 周金国, 等. 基于M残差的方差分量估计[J]. 重庆建筑大学学报, 2007,29(6):62-66.) |
| [20] | LI Bofeng, SHEN Yunzhong. Equivalent Residual Product Based Outlier Detection for Variance and Covariance Component Estimation[J]. Acta Geodaetica et Cartographica Sinica, 2011,40(1):10-14.(李博峰, 沈云中. 基于等效残差积探测粗差的方差-协方差分量估计[J]. 测绘学报, 2011,40(1):10-14) |
| [21] | SHANG Yanliang. Analysis and Improvement for the Convergence Performance in the Estimation of Variance Components[J]. Science of Surveying and Mapping, 2010,35(6):129-130.(尚艳亮. 方差分量估计迭代算法收敛性分析与改进[J]. 测绘科学, 2010,35(6):129-130.) |
| [22] | LIU Changjian, CHAI Hongzhou, WU Hongju, et al. Pilot Study on the Precondition of Variance Components Estimation[J]. Science of Surveying and Mapping, 2009,34(3):64-67.(刘长建, 柴洪洲, 吴洪举, 等. 方差分量估计前提初探[J]. 测绘科学, 2009,34(3):64-67) |
| [23] | LU Guifu, WANG Yong, ZOU Jian. A Fast Implementation of Null Space Based Algorithm[J]. Journal of Xi'an Jiaotong University, 2012,46(2):59-63.(卢桂馥, 王勇, 邹健. 一种快速的零空间算法[J]. 西安交通大学学报, 2012,46(2):59-63.) |