秩亏往往意味着线性高斯⁃马尔可夫模型中缺少平差必要的信息，而先验信息可以解决这些信息缺失导致的秩亏问题，提供线性模型平差的基准。秩亏自由网平差^[1]在大地测量中应用广泛，尤其在水准网、GPS基线测量以及大地坐标系框架体系建设^[2]等存在缺少基准或起算数据时。利用先验信息解决秩亏问题也十分常见，但Meissl^{[3, 4]}以及Neimei⁃ er ^[5]提出利用新增的基准确定点存在不确定性，所以应将先验信息表示为带随机性质的先验信息。添加的先验信息作为一种约束需要满足正交方程有解并且不改变残差的值都可以认为是合理的^{[6, 7]}。利用先验信息使正交方程有解时应该将模型表示为混合模型^[8]，或者扩展的高斯⁃马尔可夫模型^[9]，这种混合模型的最小二乘估计通常称作最小二乘配置^{[10, 11]}，并在确定重力场的研究中广泛使用。

先验随机信息的可信度是解决秩亏平差时一个重要的问题，Schaffrin^[12]提出一种稳健配置方法用于观测值信息和随机先验信息更合理的使用。

综上所述，利用随机先验信息克服秩亏问题的方差分量估计在秩亏网平差中的研究十分重要。本文利用最优不变二次无偏估计原理^[13-16]，推导带有先验信息的高斯⁃马尔可夫模型的方差分量估计公式，从而做到对观测信息和随机先验信息进行合理配权，以提高秩亏网中基准的可靠程度。

1 最优不变二次无偏估计方差分量

最优不变二次无偏估计是方差分量估计中广泛应用的估计方法，利用该方法估计随机模型中先验信息部分和观测值部分的方差分量。最优不变二次无偏估计利用观测值二次型$\overset\frown{\boldsymbol{y}}^{{T}} \boldsymbol{D} \overset\frown{\boldsymbol{y}}$估计方差分量，要求方差分量估计满意无偏性、不变性、最小方差，将上述3个条件整合成为一个极值求解问题，估计出满足要求的未知的D矩阵，进而估计方差分量。因为$\overset\frown{\boldsymbol{y}}^{\mathrm{T}} \boldsymbol{D} \overset\frown{\boldsymbol{y}}= \overset\frown{\boldsymbol{y}}^{\mathrm{T}} \boldsymbol{D}^{\mathrm{T}} \overset\frown{\boldsymbol{y}} = \overset\frown{\boldsymbol{y}}^{\mathrm{T}}\left( {\left( {\boldsymbol{D}{\rm{ + }}{{\boldsymbol{D}}^{\rm{T}}}} \right)/2} \right)\overset\frown{\boldsymbol{y}}$所以规定未知的n × n的D矩阵为对称矩阵。

1）无偏性：

$ E\left(\overset\frown{\boldsymbol{y}}^{\mathrm{T}} \boldsymbol{D} \overset\frown{\boldsymbol{y}}\right)=\boldsymbol{p}^{\mathrm{T}} \sigma $

(1)

式中，p =(p_i)，i∈{1，2}，σ₁²tr(DV₁)+σ₂²tr(DV₂)=p₁σ₁²+p₂σ₂²，tr(DV_i)=p_i。根据二次型期望公示将式（10）表示为：

$ E\left(\overset\frown{\boldsymbol{y}}^{T}\boldsymbol D \overset\frown{\boldsymbol{y}}\right)=\sum\limits_{i=1}^{2} \sigma_{i}^{2} \operatorname{tr}\left(\boldsymbol D \boldsymbol V_{i}\right) $

(2)

2）不变性：

$ \boldsymbol D\boldsymbol A_{2}=0 $

(3)

3）最小方差：

$ V\left(\overset\frown{\boldsymbol{y}}^{\mathrm{T}} D \overset\frown{\boldsymbol{y}}\right)=\min $

(4)

当$\overset\frown{\boldsymbol{y}}$满足正态分布时，由最小方差条件以及随机向量的方差性质，详见参考文献[16]，得到：

$ \boldsymbol{V}\left(\overset\frown{\boldsymbol{y}}^{\mathrm{T}} \boldsymbol{D} \overset\frown{\boldsymbol{y}}\right)=2 \operatorname{tr}(\boldsymbol{D} \boldsymbol{\varSigma} \boldsymbol{D} \boldsymbol{\varSigma})+4 \xi_{2}^{\mathrm{T}} \boldsymbol{A}_{2}^{\mathrm{T}} {\boldsymbol{D}} {\boldsymbol{A}}_{2} {\xi}_{2} $

(5)

将不变性式（3）带入式（5）得到：

$ \boldsymbol{V}\left(\overset\frown{\boldsymbol{y}}^{\mathrm{T}} \boldsymbol{D} \overset\frown{\boldsymbol{y}}\right)=2 \operatorname{tr}(\boldsymbol{D} \boldsymbol{\varSigma} \boldsymbol{D} \boldsymbol{\varSigma}) $

(6)

由$\boldsymbol{V}\left(\overset\frown{\boldsymbol{y}}^{\mathrm{T}} \boldsymbol{D} \overset\frown{\boldsymbol{y}}\right)=2 \operatorname{tr}(\boldsymbol{D} \boldsymbol{\varSigma} \boldsymbol{D} \boldsymbol{\varSigma})$做目标函数，DA₂=0，tr(DV_i)=p_i做约束条件得到拉格朗日函数：

$ \begin{gathered} w(\boldsymbol{D})=2 \operatorname{tr}(\boldsymbol{D} \boldsymbol{\varSigma} \boldsymbol{D} \boldsymbol{\varSigma})-4 \operatorname{tr}\left(\boldsymbol{D} \boldsymbol{A}_{2} \boldsymbol{\varLambda}^{\mathrm{T}}\right)- \\ 4 \sum\limits_{i=1}^{2} \lambda_{i}\left(\operatorname{tr}\left(\boldsymbol{D V}_{i}\right)-p_{i}\right) \end{gathered} $

(7)

式中，$\mathop { - 4\mathit{\boldsymbol{ \boldsymbol{\varLambda} }}}\limits_{n \times \left( {m - {m_1}} \right)} $表示DA₂=0约束的拉格朗日乘子矩阵；-4λ₁、-4λ₂表示tr(DV_i)=p_i约束的两个拉格朗日乘子；因为D为对称矩阵所以将D矩阵表达为D=(F+F^T) /2，带入式（7）中，有

$ \begin{gathered} w\left(\left(\boldsymbol{F}+\boldsymbol{F}^{\mathrm{T}}\right) / 2\right)=\operatorname{tr}(\boldsymbol{F} \boldsymbol{\varSigma} \boldsymbol{F} \boldsymbol{\varSigma})+\operatorname{tr}\left(\boldsymbol{F} \boldsymbol{\varSigma} \boldsymbol{F}^{\mathrm{T}} \boldsymbol{\varSigma}\right)- \\ 2 \operatorname{tr}\left(\boldsymbol{F} \boldsymbol{A}_{2} \boldsymbol{\varLambda}^{\mathrm{T}}\right)-2 \operatorname{tr}\left(\boldsymbol{F}^{\mathrm{T}} \boldsymbol{A}_{2} \boldsymbol{\varLambda}^{\mathrm{T}}\right)- \\ 2 \sum\limits_{i=1}^{2} \lambda_{i}\left(\operatorname{tr}\left(\left(\boldsymbol{F}+\boldsymbol{F}^{\mathrm{T}}\right) \boldsymbol{V}_{i}\right)-2 p_{i}\right) \end{gathered} $

由∂tr(ABA^TC)/∂A=(BA^TC)^T+CAB以及∂w((F+F^T)/2)/∂F=0，得到：

$ 2 \boldsymbol{\varSigma} \boldsymbol{F}^{\mathrm{T}} \boldsymbol{\varSigma}+2 \boldsymbol{\varSigma} \boldsymbol{F} \boldsymbol{\varSigma}-2 \boldsymbol{\varLambda} \boldsymbol{A}_{2}^{\mathrm{T}}-2 \boldsymbol{A}_{2} \boldsymbol{\varLambda}^{\mathrm{T}}-4 \sum\limits_{i=1}^{2} \lambda_{i} \boldsymbol{V}_{i}=0 $

(8)

将D=(F+F^T) /2带入式（8）中得到：

$ \boldsymbol{\varSigma} \boldsymbol{D} \boldsymbol{\varSigma}-\frac{1}{2}\left(\boldsymbol{\varLambda} \boldsymbol{A}_{2}^{\mathrm{T}}+\boldsymbol{A}_{2} \boldsymbol{\varLambda}^{\mathrm{T}}\right)-\sum\limits_{i=1}^{2} \lambda_{i} \boldsymbol{V}_{i}=0 $

(9)

R(A₂) 代表由A₂列向量形成的列空间，A₂(A₂^TΣ^-1A₂)^-1A₂^TΣ^-1是到R(A₂)上的广义Σ^-1正交投影算子^[21]，该算子满足条件(I-A₂(A₂^TΣ^-1A₂^-1A₂^TΣ^-1)A₂=0。实际上，任何满足矩阵乘法的矩阵X满足(X (I-A₂(A₂^TΣ^-1A₂)-1 A₂^TΣ^-1))A₂=0。为消去式（9）中ΛA₂^T+A₂Λ^T，(X (I-A₂(A₂^TΣ^-1A₂)^-1A₂^TΣ^-1)) 取为对称矩阵，X=Σ^-1时式（9）为对称矩阵，该对称矩阵记为W。

W=Σ^-1 (I-A₂(A₂^TΣ^-1A₂)^-1A₂^TΣ^-1)，满足WA₂=0和WΣDΣW=D。将式（9）左右乘以W矩阵，结合W满足的上述两个等式，可得矩阵D的估计值，记作$\hat{\boldsymbol D}$：

$ \hat{\boldsymbol D}=\sum\limits_{i=1}^{2} \lambda_{i} \boldsymbol{W} \boldsymbol{V}_{i} \boldsymbol{W} $

(10)

又由tr(DV_i)=p_i得到：

$ \sum\limits_{i=1}^{2} \lambda_{i} \operatorname{tr}\left(\boldsymbol{W} \boldsymbol{V}_{i} \boldsymbol{W} \boldsymbol{V}_{j}\right)=p_{j} $

(11)

$\boldsymbol{p}^{\mathrm{T}} \hat{\boldsymbol{\sigma}}$为p^Tσ的估计值，得到：

$ \boldsymbol{p}^{\mathrm{T}} \hat{\boldsymbol{\sigma}}=\overset\frown{\boldsymbol{y}}^{\mathrm{T}} \hat{\boldsymbol{D}} \overset\frown{\boldsymbol{y}}=\boldsymbol\lambda^{\mathrm{T}} \boldsymbol{q} $

(12)

式中，$\mathit{\boldsymbol{\lambda }} = \left( {{\lambda _i}} \right), \mathit{\boldsymbol{q}} = \left( {{q_i}} \right) = \left( {{{\overset\frown{\boldsymbol{y}}}^{\rm{T}}}\mathit{\boldsymbol{W}}{\mathit{\boldsymbol{V}}_i}\mathit{\boldsymbol{W}}{\overset\frown{\boldsymbol{y}}} } \right)$。记S=(tr (WV_iWV_j))，则由式（11）有S^Tλ=p。因此$\boldsymbol{p}^{\mathrm{T}} \hat{\boldsymbol{\sigma}}$估计值为：

$ \boldsymbol{p}^{\mathrm{T}} \hat{\boldsymbol{\sigma}}=\boldsymbol{p}^{\mathrm{T}} \boldsymbol{S}^{-} \boldsymbol{q} $

(13)

式中，S^-为S的广义逆矩阵，需要证明式（13）得出的方差分量估计式唯一的，也即证明方差分量估计独立于广义逆S^-的选取。利用Cholesky分解将Σ分解为Σ=GG^T，$\mathop H\limits_{{n^2} \times 2} $的列向量为：

$ \boldsymbol H=\left(h_{i}\right)=\operatorname{vec}\left(\boldsymbol{G}^{\mathrm{T}} \boldsymbol{W} \boldsymbol{V}_{i} \boldsymbol{W G}\right) $

(14)

则

$ \boldsymbol{h}_{i}^{\mathrm{T}} \boldsymbol{h}_{j}=\operatorname{tr}\left(\boldsymbol{G}^{\mathrm{T}} \boldsymbol{W} \boldsymbol{V}_{i} \boldsymbol{W} \boldsymbol{G} \boldsymbol{G}^{\mathrm{T}} \boldsymbol{W} \boldsymbol{V}_{j} \boldsymbol{W} \boldsymbol{G}\right) $

(15)

式中，W、Σ满足WΣW=W，矩阵乘积求迹运算满足tr(AB)=tr(BA)，所以h_i^Th_j=tr(WV_iWV_j)，也即H^TH=S。

定义矩阵$\mathit{\boldsymbol{k}} = {\rm{vec}}\left( {{{\boldsymbol{G}}^{-1}}{\overset\frown{\boldsymbol{y}}} {{\hat {\boldsymbol y}}^{\rm{T}}}{{\left( {{{\boldsymbol{G}}^{\rm{T}}}} \right)}^{-1}}} \right)$，$\mathit{\boldsymbol{h}}_i^{\rm{T}}\mathit{\boldsymbol{k = }}{\rm{tr}}\left( {{{\hat {\boldsymbol y}}^{\rm{T}}}\mathit{\boldsymbol{W}}{\mathit{\boldsymbol{V}}_i}\mathit{\boldsymbol{W}}{\overset\frown{\boldsymbol{y}}} } \right) = {\mathit{\boldsymbol{q}}_i}$，也即H^Tk=q。因为H^TH=S以及H^Tk=q，所以q∈R (H) =R (H ^T H)=R (S)。由S^Tλ=p，所以p∈R (S)。存在两向量v、u使得p=Sv和q=Su存在。S为对称矩阵满足S^T=S且广义逆满足AA^-A=A，则$\boldsymbol{p}^{\mathrm{T}} \hat{\boldsymbol{\sigma}}= {\mathit{\boldsymbol{v}}^{\rm{T}}}\mathit{\boldsymbol{Su}}, {\mathit{\boldsymbol{p}}^{\rm{T}}}\mathit{\boldsymbol{\hat \sigma }}$估计值独立于S^-的选择，因此$\boldsymbol{p}^{\mathrm{T}} \hat{\boldsymbol{\sigma}}$的估计唯一。

当S可逆时，则p^T=[1, 0]^T∈R (S) 且p^T=[0, 1]^T∈R (S)，由式（13）可知方差分量的估值唯一确定为：

$ \hat{\boldsymbol\sigma}=\boldsymbol{S}^{-1} \boldsymbol q $

(16)

当S奇异时，则p^T =[1, 0]^T∈R (S) 或者p^T = [0, 1]^T∈R(S)，由式（13）可知方差分量的估值唯一确定为：

$ \hat{\boldsymbol\sigma}=\boldsymbol S^{-}\boldsymbol q $

(17)

利用迭代算法，第i-1次估计求得的${{\mathit{\boldsymbol{\hat \sigma }}}^{i - 1}}$乘以对应的V^i-1得到Vⁱ，再进行第i次估计，直至$\max \left( {\left\| {{{\left( {\hat \sigma _1^2} \right)}^i} - 1} \right\|, \left\| {{{\left( {\hat \sigma _2^2} \right)}^i} - 1} \right\|} \right) < \varepsilon $（ε是给出的极小的正数阈值），确保方差分量每一个估计值都满足迭代终止条件。为了衡量所估计的方差分量的精度，方差分量的方差也需要被估计。

将$\hat{\boldsymbol D} = \sum\limits_{i = 1}^2 {{\lambda _i}\mathit{\boldsymbol{W}}{\mathit{\boldsymbol{V}}_i}\mathit{\boldsymbol{W}}} $带入式（18）中。

$ \begin{array}{l} \ \ \ \ \ \ \ \ \boldsymbol{V}\left(\boldsymbol{p}^{\mathrm{T}} \hat{\boldsymbol{\sigma}}\right)=2 \operatorname{tr}(\hat{\boldsymbol D} \boldsymbol{\varSigma} \hat{\boldsymbol D} \boldsymbol{\varSigma})= \\ 2 \operatorname{tr}\left(\left(\sum\limits_{i=1}^{2} \lambda_{i} \boldsymbol{W} \boldsymbol{V}_{i} \boldsymbol{W}\right) \Sigma\left(\sum\limits_{j=1}^{2} \lambda_{j} \boldsymbol{W} \boldsymbol{V}_{j} \boldsymbol{W}\right) \boldsymbol{\varSigma}\right)= \\ 2\left(\sum\limits_{i=1}^{2} \lambda_{i} \operatorname{vec} \boldsymbol{G}^{\mathrm{T}} \boldsymbol{W} \boldsymbol{V}_{i} \boldsymbol{W} \boldsymbol{G}\right)^{\mathrm{T}}\left(\sum\limits_{j=1}^{2} \lambda_{j} \operatorname{vec} \boldsymbol{G}^{\mathrm{T}} \boldsymbol{W} \boldsymbol{V}_{j} \boldsymbol{W G}\right)= \\ 2(\boldsymbol{H} \boldsymbol{\lambda})^{\mathrm{T}} \boldsymbol{H} \boldsymbol{\lambda}=2 \boldsymbol{\lambda}^{\mathrm{T}} \boldsymbol{S} \boldsymbol{\lambda}=2\left(\boldsymbol{\lambda}^{\mathrm{T}} \boldsymbol{S}\right) \boldsymbol{S}^{-}(\boldsymbol{S} \boldsymbol{\lambda})=2 \boldsymbol{p}^{\mathrm{T}} \boldsymbol{S}^{-} \boldsymbol{p} \end{array} $

(18)

当S可逆时，则p^T=[1, 0]^T∈R(S) 且p^T=[0, 1]^T∈R(S))，由式（17）可知方差分量的方差估值唯一确定为：

$ \boldsymbol V(\hat{\boldsymbol\sigma})=2 \boldsymbol{S}^{-1} $

(19)

当S奇异时且p^T=[1, 0]^T∈R(S)成立时，由式（19）可知方差分量的方差估值唯一确定为$V\left( {\hat \sigma _1^2} \right) = 2{s_{11}}$；S奇异且p^T=[0, 1]^T∈R(S)成立时，由式（19）可知方差分量的方差估值唯一确定为$V\left( {\hat \sigma _2^2} \right)=2{s_{22}}$。其中S^-1=(s_ij)。

2 实例验证

如图 1所示，秩亏水准网中现有两部分数据，一部分为观测数据，另一部分为由之前平差成果中获得的先验信息数据。A~E为水准点。高差测量数据如表 1所示，其中高差观测值方差计算公式为：σ_i²=(0.003 m) ² × D_i（D_i表示相应高程差对应测量时的路程长度）。A点先验高程为5. 016m，B点先验高程为6. 017m，A点和B点先验高程的方差协方差矩阵为diag (1.5² 2²)（mm²）。

由表 2可知，根据误差向量估计值，算出方差估计值为${{\hat \sigma }_0}$=1.886 3 mm。

计算中发现自第二次迭代后均收敛于1，于是在第二次计算完成后结束迭代计算。如表 3、表 4所示，得到最终估计的方差分量为5. 452 3×10^-6和3. 949 5×10^-7，两部分方差分量的比值为：0. 072 4。相比于最小二乘估计单位权方差，方差分量估计更加准确的描述了带有先验信息的高斯⁃马尔可夫的随机模型，这也使得先验信息中的绝对基准更加合理可靠地应用在秩亏平差中。

3 结束语

新增的先验信息常用作确定基准，但新增的先验信息带有随机性质，这意味着先验信息的可靠程度无法衡量。先验信息在平差模型中的不确定性需要寻求一种合理可靠的方法来度量先验信息和观测值，对两部分信息进行合理地分配权重，需要得到平差模型中准确的验后随机模型。本文利用最优不变二次无偏估计先验信息和观测值两部分的方差分量以及方差分量方差，得到验后随机模型，准确地描述两部分的权重。然后可以通过合理配权，从而进一步提高先验信息在解决秩亏网平差时的可靠性。