平差随机模型的验后估计
概述
众所周知,一个平差问题必须首先建立该平差问题的数学模型。平差的数学模型包括函数 模型和随机模型两类。描述平差问题中观测量与观测量之间、观测量与未知参数之间相互关系 的函数表达式,称平差函数模型。随机模型是描述观测误差 \(\bm{\Delta}\) 的一些随机特征,在 平差中主要是 \(\bm{\Delta}\) 的数学期望和方差,具有 \[E(\bm{\Delta})=\bm{0} \tag{3-1-1}\] 和 \[D(\bm{\Delta})=\sigma^{2}\bm{Q}=\sigma^{2}\bm{P}^{-1} \tag{3-1-2}\] (3-1-1)式表明观测误差中不含系统误差和粗差,是一般情况下最小二乘平差的要求。 (3-1-2)式是平差时定权的根据。
注意 (3-1-2) 式只确定了权的相对比例:方差阵是乘积 \(\sigma_{0}^{2}\bm{P}^{-1}\),若把权阵 \(\bm{P}\) 整体乘以常数 \(k\)、同时把单位权方差 \(\sigma_{0}^{2}\) 换成 \(\sigma_{0}^{2}/k\),方差阵不变。因此"定权"本质上定的是各类观测之间的权比,\(\sigma_{0}^{2}\) 只是吸收共同比例因子的标量。这正是后文赫尔默特迭代中 (3-2-22) 式的常数 \(c\) 可以任取(例如令某类观测的权保持为 1,见例 3-2-1 中 \(P_{\beta}=1\))而不影响平差结果的根本原因。
平差前,随机模型要已知 \(D(\bm{\Delta})\),称为验前方差。只有精确地已知验前方差 \(D(\bm{\Delta})\) 才能精确地定权,所以随机模型的估计,就是验前方差 \(D(\bm{\Delta})\) 的估计,也就是观测值权的估计。
过去很长的时间,平差都在单一的同类观测量中进行,例如测角网平差、水准网平差。定权 可从定义式(3-1-2)出发,采用测量平差中常用方法定权,例如,水准高差按路线长度倒数定 权等。随着平差对象从单一同类观测量扩展为不同类的多种观测量,一般,它们的验前方差又不 能都已知,如何能精确地估计它们的方差,达到精确地定权就需要深入研究了。所以,国内外测 量界把平差随机模型的估计作为主要课题进行研究,取得了丰富的成果。
对不同类的观测量,一般采用经验公式定权,即根据仪器出厂标明的标称精度估算各自的 方差,然后再按定义式(3-1-2)定权。例如在边角同测的控制网中,测距仪给出的测边中误差 标称精度公式为 \[\sigma_{s_i}=a+bs_i\] 测角中误差为 \(\sigma_{\beta}\)(按规范),以 \(\sigma_{\beta}\) 和 \(\sigma_{s_i}\) 为测边和 测角的验前方差定权,得 \[P_{\beta}=1,\qquad P_{s_i}=\frac{\sigma_{\beta}^{2}}{\sigma_{s_i}^{2}} \qquad\text{(单位:}('')^{2}/\mathrm{cm}^{2}\text{)}\]
此处隐去了单位权方差的选取:按定义式(3-1-2)令角度观测的权为 \(1\)(即取 \(\sigma_0=\sigma_{\beta}\)),故 \(P_{\beta}=1\), \(P_{s_i}=\sigma_{\beta}^{2}/\sigma_{s_i}^{2}\)。
在卫星网与地面网、重力网与水准网的联合平差,摄影测量与大地测量数据联合处理中, 也可按上述经验公式的方法定权。
这种估计验前方差确定各类观测量权的方法,实践证明,在许多情况下是不够精确的。为 了提高方差估计的精度,20 世纪 70 年代开始出现了用验后的方法估计各类观测量的方差,然 后定权,我们称为平差随机模型的验后估计法。
随机模型的验后估计,其基本思想是,先对各类观测量定初权,进行预平差,利用预平差后 得到的信息,主要是各类观测值的改正数 \(\bm{V}\),依据一定的原则对各类观测量的验前方差 和协方差作出估计,依此定权。实践已经证明这种定权方法的优越性,并已在实际工作中广泛应 用。
验前定权与验后估计是一条"闭环":验前定权靠仪器标称精度这类经验值,验后估计则先按初权做一次预平差,再用预平差得到的改正数 \(\bm{V}\) 反推各类观测的真实方差,据此重新定权、重新平差。可以把 \(\bm{V}\) 想成"数据对初权不满意的申诉"——某类观测的残差平方和 \(\bm{V}^{T}\bm{P}\bm{V}\) 偏大,说明该类的权定得偏大(即低估了其方差),反之亦然。赫尔默特方差分量估计做的就是把这句"申诉"翻译成方差的修正量;本章 3-2 节给出了翻译所需的全部数学工具。
本章首先介绍方差估计法,即赫尔默特估计法,并介绍该法在实际计算中的一些简化计算公 式。接着介绍方差、协方差估计法。然后介绍二次无偏估计法,主要介绍 C. R. Rao 于 1970 年提出的最小范数二次无偏估计(MINQUE)法和 K. R. Koch 于 1980 年提出的最优不变二次 无偏估计(BIQUE)法。最后介绍方差分量估计中的精度评定以及方差分量估计在测量实践中的 应用。
本章的"函数模型/随机模型"划分,对应《最优估计基础》第1章"期望和方差"(1.2 节)与"多维随机变量"(1.4 节)中的统计基础:函数模型给出 \(E(\bm{L})=\bm{B}\bm{X}\)(该书第2章的高斯-马尔可夫模型),随机模型给出 \(D(\bm{\Delta})=\sigma_{0}^{2}\bm{Q}\),即该书第1章对观测噪声的统计描述;最小二乘平差下单位权方差估值 \(\hat{\sigma}_{0}^{2}=\bm{V}^{T}\bm{P}\bm{V}/(n-t)\) 的统计含义,见该书第2章"最小二乘估计"(2.2 节)。
赫尔默特方差估计法
间接平差时的方差分量估计
利用预平差的改正数 \(\bm{V}\),按验后估计各类观测量验前方差的方法,其思想最早是由赫尔默特(1924)提出的。若各类观测量之间相互独立,即观测量的方差阵是拟对角型矩阵,称为方差估计,或称方差分量估计。以下介绍由 Welsch(1978)推证的赫尔默特方差估计严密公式。
间接平差的基本公式为
函数模型: \[\underset{n\times 1}{\bm{L}}=\underset{n\times t}{\bm{B}}\,\underset{t\times 1}{\bm{X}}+\underset{n\times 1}{\bm{\Delta}} \tag{3-2-1}\]
随机模型: \[\begin{aligned} &E(\bm{L})=\bm{B}\widetilde{\bm{X}},\quad E(\bm{\Delta})=\bm{0},\\ &D(\bm{L})=\sigma_{0}^{2}\bm{P}^{-1},\quad D(\bm{\Delta})=D(\bm{L})=\sigma_{0}^{2}\bm{P}^{-1} \end{aligned} \tag{3-2-2}\]
误差方程: \[\bm{V}=\bm{B}\hat{\bm{X}}-\bm{L} \tag{3-2-3}\]
法方程及其解为 \[\bm{N}\hat{\bm{X}}=\bm{W},\qquad \hat{\bm{X}}=\bm{N}^{-1}\bm{W} \tag{3-2-4}\] 式中 \(\bm{N}=\bm{B}^{T}\bm{P}\bm{B}\),\(\bm{W}=\bm{B}^{T}\bm{P}\bm{L}\)。
现设在 \(\bm{L}\) 中包含有两类相互独立的观测值 \(\underset{n_{1}\times 1}{\bm{L}_{1}}\) 和 \(\underset{n_{2}\times 1}{\bm{L}_{2}}\),它们的权阵分别为 \(\underset{n_{1}\times n_{1}}{\bm{P}_{1}}\) 和 \(\underset{n_{2}\times n_{2}}{\bm{P}_{2}}\),并且 \(\bm{P}_{12}=\bm{0}\),它们的误差方程分别为 \[\left.\begin{aligned} \bm{V}_{1}&=\bm{B}_{1}\hat{\bm{X}}-\bm{L}_{1}\\ \bm{V}_{2}&=\bm{B}_{2}\hat{\bm{X}}-\bm{L}_{2} \end{aligned}\right\} \tag{3-2-5}\] 且有下列关系式: \[\begin{gathered} \bm{L}=\begin{bmatrix}\bm{L}_{1}\\ \bm{L}_{2}\end{bmatrix},\quad \bm{V}=\begin{bmatrix}\bm{V}_{1}\\ \bm{V}_{2}\end{bmatrix},\quad \bm{B}=\begin{bmatrix}\bm{B}_{1}\\ \bm{B}_{2}\end{bmatrix},\quad \bm{P}=\begin{bmatrix}\bm{P}_{1}&\bm{0}\\ \bm{0}&\bm{P}_{2}\end{bmatrix}\\ \bm{N}=\bm{B}^{T}\bm{P}\bm{B}=\bm{B}_{1}^{T}\bm{P}_{1}\bm{B}_{1}+\bm{B}_{2}^{T}\bm{P}_{2}\bm{B}_{2}=\bm{N}_{1}+\bm{N}_{2}\\ \bm{W}=\bm{B}^{T}\bm{P}\bm{L}=\bm{B}_{1}^{T}\bm{P}_{1}\bm{L}_{1}+\bm{B}_{2}^{T}\bm{P}_{2}\bm{L}_{2}=\bm{W}_{1}+\bm{W}_{2} \end{gathered} \tag{3-2-6}\]
一般地说,第一次平差时给定的两类观测值的权 \(\bm{P}_{1}\) 和 \(\bm{P}_{2}\) 是不恰当的,或者说它们所对应的单位权方差不相等,令其分别为 \(\sigma_{0_{1}}^{2}\) 和 \(\sigma_{0_{2}}^{2}\),则有 \[\left.\begin{aligned} D(\bm{L}_{1})&=\sigma_{0_{1}}^{2}\bm{P}_{1}^{-1}\\ D(\bm{L}_{2})&=\sigma_{0_{2}}^{2}\bm{P}_{2}^{-1} \end{aligned}\right\} \tag{3-2-7}\]
方差分量估计的目的是利用各次平差后各类改正数的平方和 \(\bm{V}_{1}^{T}\bm{P}_{1}\bm{V}_{1}\) 及 \(\bm{V}_{2}^{T}\bm{P}_{2}\bm{V}_{2}\) 来估计 \(\sigma_{0_{1}}^{2}\) 及 \(\sigma_{0_{2}}^{2}\)。为此,必须建立残差平方和与 \(\sigma_{0_{1}}^{2}\) 及 \(\sigma_{0_{2}}^{2}\) 之间的关系式。
赫尔默特估计的迭代直觉:先按初权做一次平差,得到各类观测各自的残差平方和 \(\bm{V}_{1}^{T}\bm{P}_{1}\bm{V}_{1}\)、\(\bm{V}_{2}^{T}\bm{P}_{2}\bm{V}_{2}\)。若初权定得恰当,两类观测应共享同一个单位权方差;若它们各自的单位权方差 \(\sigma_{0_{1}}^{2}\neq\sigma_{0_{2}}^{2}\),就说明权比定错了。关键一步是 (3-2-8) 式:\(\bm{V}_{i}^{T}\bm{P}_{i}\bm{V}_{i}\) 的数学期望是 \(\sigma_{0_{1}}^{2}\)、\(\sigma_{0_{2}}^{2}\) 的线性函数,把期望符号去掉、解方程组 \(\bm{S}\hat{\bm{\theta}}=\bm{W}_{\theta}\),就得到 \(\sigma_{0_{i}}^{2}\) 的估值;再按 (3-2-22) 式重新定权、重新平差,反复迭代直到各 \(\sigma_{0_{i}}^{2}\) 相等。整个流程就是"平差 \(\to\) 用改正数估计方差 \(\to\) 按新方差定权 \(\to\) 再平差"的闭环。
因为对于数学期望为 \(\bm{\eta}\),方差阵为 \(\bm{\Sigma}\) 的随机向量 \(\bm{Y}\),其二次型 \(\bm{Y}^{T}\bm{M}\bm{Y}\)(\(\bm{M}\) 为任一对称可逆阵)的数学期望为 \[E(\bm{Y}^{T}\bm{M}\bm{Y})=\operatorname{tr}(\bm{M}\bm{\Sigma})+\bm{\eta}^{T}\bm{M}\bm{\eta} \tag{3-2-8}\] 而改正数 \(\bm{V}\) 的期望为零,即有 \[E(\bm{V}_{1})=\bm{0} \tag{3-2-9}\] 所以 \[E(\bm{V}_{1}^{T}\bm{P}_{1}\bm{V}_{1})=\operatorname{tr}(\bm{P}_{1}D(\bm{V}_{1})) \tag{3-2-10}\] 式中 \(D(\bm{V}_{1})\) 为改正数 \(\bm{V}_{1}\) 的方差。
(3-2-8) 式二次型期望公式 \(E(\bm{Y}^{T}\bm{M}\bm{Y})=\operatorname{tr}(\bm{M}\bm{\Sigma})+\bm{\eta}^{T}\bm{M}\bm{\eta}\) 的来历:由 \(E(\bm{Y}\bm{Y}^{T})=\bm{\Sigma}+\bm{\eta}\bm{\eta}^{T}\) 与"标量的迹等于它自己、期望与迹可交换"即可推出(严格的证明见 3-5 节 (3-5-3) 式)。对改正数恒有 \(E(\bm{V})=\bm{0}\)((3-2-9) 式),故 (3-2-10) 式 \(E(\bm{V}_{1}^{T}\bm{P}_{1}\bm{V}_{1})=\operatorname{tr}(\bm{P}_{1}D(\bm{V}_{1}))\)——二次型期望只剩第一项 \(\operatorname{tr}(\bm{M}\bm{\Sigma})\),问题的核心于是转为求 \(D(\bm{V}_{1})\),即 (3-2-11)、(3-2-12) 两式所做的事。
由 (3-2-5) 式可知: \[\begin{aligned} \bm{V}_{1}&=\bm{B}_{1}\hat{\bm{X}}-\bm{L}_{1} =\bm{B}_{1}\bm{N}^{-1}\bm{W}-\bm{L}_{1} =\bm{B}_{1}\bm{N}^{-1}(\bm{W}_{1}+\bm{W}_{2})-\bm{L}_{1}\\ &=(\bm{B}_{1}\bm{N}^{-1}\bm{B}_{1}^{T}\bm{P}_{1}-\bm{E})\bm{L}_{1} +\bm{B}_{1}\bm{N}^{-1}\bm{B}_{2}^{T}\bm{P}_{2}\bm{L}_{2} \end{aligned} \tag{3-2-11}\] 故 \(\bm{V}_{1}\) 的方差为 \[\begin{aligned} D(\bm{V}_{1})={}&(\bm{B}_{1}\bm{N}^{-1}\bm{B}_{1}^{T}\bm{P}_{1}-\bm{E})\,D(\bm{L}_{1})\,(\bm{B}_{1}\bm{N}^{-1}\bm{B}_{1}^{T}\bm{P}_{1}-\bm{E})^{T}+{}\\ &\bm{B}_{1}\bm{N}^{-1}\bm{B}_{2}^{T}\bm{P}_{2}D(\bm{L}_{2})\bm{P}_{2}\bm{B}_{2}\bm{N}^{-1}\bm{B}_{1}^{T} \end{aligned}\] 将上式展开,并顾及 \(D(\bm{L}_{1})=\sigma_{0_{1}}^{2}\bm{P}_{1}^{-1}\),\(D(\bm{L}_{2})=\sigma_{0_{2}}^{2}\bm{P}_{2}^{-1}\),则上式可整理得: \[\begin{aligned} D(\bm{V}_{1})={}&\sigma_{0_{1}}^{2}\left(\bm{B}_{1}\bm{N}^{-1}\bm{N}_{1}\bm{N}^{-1}\bm{B}_{1}^{T} -2\bm{B}_{1}\bm{N}^{-1}\bm{B}_{1}^{T}+\bm{P}_{1}^{-1}\right)+{}\\ &\sigma_{0_{2}}^{2}\left(\bm{B}_{1}\bm{N}^{-1}\bm{N}_{2}\bm{N}^{-1}\bm{B}_{1}^{T}\right) \end{aligned} \tag{3-2-12}\]
由 \(D(\bm{V}_1)\) 展开到 (3-2-12) 式的关键是:交叉项合并为 \(-2\bm{B}_1\bm{N}^{-1}\bm{B}_1^{T}\sigma_{0_1}^2\),并利用了 \(\bm{P}_1\bm{B}_1\bm{N}^{-1}\bm{B}_1^{T}\bm{P}_1\bm{P}_1^{-1}=\bm{B}_1\bm{N}^{-1}\bm{B}_1^{T}\) 与 \(\bm{B}_1^{T}\bm{P}_1\bm{B}_1=\bm{N}_1\),把含 \(\bm{B}_1\bm{N}^{-1}\bm{B}_1^{T}\bm{P}_1\) 的平方项化为 \(\bm{B}_1\bm{N}^{-1}\bm{N}_1\bm{N}^{-1}\bm{B}_1^{T}\)。
将上式代入 (3-2-10) 式,得: \[\begin{aligned} E(\bm{V}_{1}^{T}\bm{P}_{1}\bm{V}_{1})={}&\operatorname{tr}(\bm{P}_{1}D(\bm{V}_{1}))\\ ={}&\sigma_{0_{1}}^{2}\operatorname{tr}\{\bm{P}_{1}\bm{P}_{1}^{-1} -2\bm{P}_{1}\bm{B}_{1}\bm{N}^{-1}\bm{B}_{1}^{T} +\bm{P}_{1}\bm{B}_{1}\bm{N}^{-1}\bm{N}_{1}\bm{N}^{-1}\bm{B}_{1}^{T}\}+{}\\ &\sigma_{0_{2}}^{2}\operatorname{tr}\{\bm{P}_{1}\bm{B}_{1}\bm{N}^{-1}\bm{N}_{2}\bm{N}^{-1}\bm{B}_{1}^{T}\}\\ ={}&\sigma_{0_{1}}^{2}\operatorname{tr}\{\underset{n_{1}\times n_{1}}{\bm{E}} -2\bm{N}^{-1}\bm{B}_{1}^{T}\bm{P}_{1}\bm{B}_{1} +\bm{N}^{-1}\bm{N}_{1}\bm{N}^{-1}\bm{B}_{1}^{T}\bm{P}_{1}\bm{B}_{1}\}+{}\\ &\sigma_{0_{2}}^{2}\operatorname{tr}\{\bm{N}^{-1}\bm{N}_{2}\bm{N}^{-1}\bm{B}_{1}^{T}\bm{P}_{1}\bm{B}_{1}\}\\ ={}&\left\{n_{1}-2\operatorname{tr}(\bm{N}^{-1}\bm{N}_{1}) +\operatorname{tr}(\bm{N}^{-1}\bm{N}_{1}\bm{N}^{-1}\bm{N}_{1})\right\}\sigma_{0_{1}}^{2}+{}\\ &\operatorname{tr}(\bm{N}^{-1}\bm{N}_{1}\bm{N}^{-1}\bm{N}_{2})\sigma_{0_{2}}^{2} \end{aligned} \tag{3-2-13}\]
上式第二步到第三步用了迹的循环性质 \(\operatorname{tr}(\bm{A}\bm{B})=\operatorname{tr}(\bm{B}\bm{A})\),把 \(\bm{P}_1\bm{B}_1\) 移到乘积末尾,使括号内化为 \(\bm{N}^{-1}\bm{N}_1\) 的形式。
把 (3-2-13) 的三项来源再拆细一层:第一项 \(\operatorname{tr}(\bm{P}_{1}\bm{P}_{1}^{-1})=n_{1}\) 来自 (3-2-12) 中 \(\sigma_{0_{1}}^{2}\) 对应的 \(\bm{P}_{1}^{-1}\) 项;交叉项 \(-2\bm{B}_{1}\bm{N}^{-1}\bm{B}_{1}^{T}\) 左乘 \(\bm{P}_{1}\) 后取迹得 \(-2\operatorname{tr}(\bm{N}^{-1}\bm{B}_{1}^{T}\bm{P}_{1}\bm{B}_{1})=-2\operatorname{tr}(\bm{N}^{-1}\bm{N}_{1})\);平方项 \(\bm{B}_{1}\bm{N}^{-1}\bm{N}_{1}\bm{N}^{-1}\bm{B}_{1}^{T}\) 左乘 \(\bm{P}_{1}\) 取迹得 \(\operatorname{tr}(\bm{N}^{-1}\bm{N}_{1}\bm{N}^{-1}\bm{N}_{1})\)。而第 2 类观测对第 1 类残差平方和的贡献只有正项 \(\operatorname{tr}(\bm{N}^{-1}\bm{N}_{1}\bm{N}^{-1}\bm{N}_{2})\sigma_{0_{2}}^{2}\),不含与自身单位权方差相关的 \(-2\) 型项。
同理可得: \[\begin{aligned} E(\bm{V}_{2}^{T}\bm{P}_{2}\bm{V}_{2})={}&\operatorname{tr}(\bm{N}^{-1}\bm{N}_{1}\bm{N}^{-1}\bm{N}_{2})\sigma_{0_{1}}^{2}+{}\\ &\left\{n_{2}-2\operatorname{tr}(\bm{N}^{-1}\bm{N}_{2})+\operatorname{tr}(\bm{N}^{-1}\bm{N}_{2})^{2}\right\}\sigma_{0_{2}}^{2} \end{aligned} \tag{3-2-14}\]
注意 (3-2-13)、(3-2-14) 中 \(\operatorname{tr}(\bm{N}^{-1}\bm{N}_{i})^{2}\) 的记法:本书用它表示 \(\operatorname{tr}(\bm{N}^{-1}\bm{N}_{i}\bm{N}^{-1}\bm{N}_{i})\)(对照 (3-2-13) 第三项写全的形式),并不是"迹的平方" \([\operatorname{tr}(\bm{N}^{-1}\bm{N}_{i})]^{2}\)。另注意 \(\bm{S}\) 的非对角元 \(\operatorname{tr}(\bm{N}^{-1}\bm{N}_{1}\bm{N}^{-1}\bm{N}_{2})\) 一般不为零——两类观测通过共同的参数解 \(\hat{\bm{X}}\) 耦合在一起,这正是方差分量必须联立方程组求解、而不能逐类单独估计的原因。若偷懒把 \(\bm{V}_{i}^{T}\bm{P}_{i}\bm{V}_{i}\) 简单地除以 \(n_{i}\)((3-2-40) 式),就丢掉了这些耦合项,得到的是有偏估计。
在 (3-2-13) 和 (3-2-14) 两式中,将数学期望的符号去掉,改成由平差得到的计算值 \(\bm{V}_{1}^{T}\bm{P}_{1}\bm{V}_{1}\) 和 \(\bm{V}_{2}^{T}\bm{P}_{2}\bm{V}_{2}\),则求出的 \(\sigma_{0_{1}}^{2}\) 和 \(\sigma_{0_{2}}^{2}\) 也应改为估值 \(\hat{\sigma}_{0_{1}}^{2}\) 和 \(\hat{\sigma}_{0_{2}}^{2}\)。将上两式写成矩阵形式为 \[\underset{2\times 2}{\bm{S}}\,\underset{2\times 1}{\hat{\bm{\theta}}}=\underset{2\times 1}{\bm{W}_{\theta}} \tag{3-2-15}\] 式中 \[\bm{S}=\begin{bmatrix} n_{1}-2\operatorname{tr}(\bm{N}^{-1}\bm{N}_{1})+\operatorname{tr}(\bm{N}^{-1}\bm{N}_{1})^{2} & \operatorname{tr}(\bm{N}^{-1}\bm{N}_{1}\bm{N}^{-1}\bm{N}_{2})\\ (\text{对\quad 称}) & n_{2}-2\operatorname{tr}(\bm{N}^{-1}\bm{N}_{2})+\operatorname{tr}(\bm{N}^{-1}\bm{N}_{2})^{2} \end{bmatrix} \tag{3-2-16}\] \[\hat{\bm{\theta}}=\begin{bmatrix}\hat{\sigma}_{0_{1}}^{2}&\hat{\sigma}_{0_{2}}^{2}\end{bmatrix}^{T},\qquad \bm{W}_{\theta}=\begin{bmatrix}\bm{V}_{1}^{T}\bm{P}_{1}\bm{V}_{1}&\bm{V}_{2}^{T}\bm{P}_{2}\bm{V}_{2}\end{bmatrix}^{T}\]
公式 (3-2-15)、(3-2-16) 即为两类观测值按间接平差时的赫尔默特估算公式。由 (3-2-15) 式知,被估参数与方程的个数相同,一般说来,有唯一解,即 \[\hat{\bm{\theta}}=\bm{S}^{-1}\bm{W}_{\theta}\]
将两类观测值扩展到 \(m\) 类观测值的一般情况,则对应的公式如下: \[\underset{n_{i}\times 1}{\bm{V}_{i}}=\underset{n_{i}\times t}{\bm{B}_{i}}\,\underset{t\times 1}{\hat{\bm{X}}}-\underset{n_{i}\times 1}{\bm{L}_{i}},\quad (i=1,2,\cdots,m) \tag{3-2-17}\] \[D(\bm{L}_{i})=\sigma_{0_{i}}^{2}\bm{P}_{i}^{-1} \tag{3-2-18}\] 将 \[\hat{\bm{X}}=\bm{N}^{-1}\bm{W}=\bm{N}^{-1}\sum_{j=1}^{m}\bm{W}_{j}\] 代入 (3-2-17) 式,并整理集项得: \[\begin{aligned} \bm{V}_{i}&=\bm{B}_{i}\hat{\bm{X}}-\bm{L}_{i} =\bm{B}_{i}\bm{N}^{-1}\sum_{j=1}^{m}\bm{W}_{j}-\bm{L}_{i}\\ &=(\bm{B}_{i}\bm{N}^{-1}\bm{B}_{i}^{T}\bm{P}_{i}-\bm{E})\bm{L}_{i} +\bm{B}_{i}\bm{N}^{-1}\sum_{\substack{j=1,\\ j\neq i}}^{m}\bm{B}_{j}^{T}\bm{P}_{j}\bm{L}_{j} \end{aligned}\] 由协方差传播律得: \[D(\bm{V}_{i})=\left(\bm{P}_{i}^{-1}+\bm{B}_{i}\bm{N}^{-1}\bm{N}_{i}\bm{N}^{-1}\bm{B}_{i}^{T} -2\bm{B}_{i}\bm{N}^{-1}\bm{B}_{i}^{T}\right)\sigma_{0_{i}}^{2} +\sum_{\substack{j=1,\\ j\neq i}}^{m}\left\{(\bm{B}_{i}\bm{N}^{-1}\bm{N}_{j}\bm{N}^{-1}\bm{B}_{i}^{T})\,\sigma_{0_{j}}^{2}\right\}\] 根据二次型的期望定理:\(E(\bm{Y}^{T}\bm{M}\bm{Y})=\operatorname{tr}(\bm{M}\bm{\Sigma})+\bm{\eta}^{T}\bm{M}\bm{\eta}\),并顾及 \(\bm{Y}=\bm{V}_{i}\),\(\bm{M}=\bm{P}_{i}\),\(\bm{\Sigma}=D(\bm{V}_{i})\) 及 \(\bm{\eta}=E(\bm{V}_{i})=\bm{0}\),则有 \[\begin{aligned} E(\bm{V}_{i}^{T}\bm{P}_{i}\bm{V}_{i})={}&\operatorname{tr}(\bm{P}_{i}\bm{P}_{i}^{-1} +\bm{P}_{i}\bm{B}_{i}\bm{N}^{-1}\bm{N}_{i}\bm{N}^{-1}\bm{B}_{i} -2\bm{P}_{i}\bm{B}_{i}\bm{N}^{-1}\bm{B}_{i}^{T})\,\sigma_{0_{i}}^{2}+{}\\ &\sum_{\substack{j=1,\\ j\neq i}}^{m}\left\{\operatorname{tr}(\bm{P}_{i}\bm{B}_{i}\bm{N}^{-1}\bm{N}_{j}\bm{N}^{-1}\bm{B}_{i}^{T})\,\sigma_{0_{j}}^{2}\right\} \end{aligned}\] 所以 \[\begin{aligned} E(\bm{V}_{i}^{T}\bm{P}_{i}\bm{V}_{i})={}&\left(n_{i}-2\operatorname{tr}(\bm{N}^{-1}\bm{N}_{i}) +\operatorname{tr}(\bm{N}^{-1}\bm{N}_{i})^{2}\right)\sigma_{0_{i}}^{2}+{}\\ &\sum_{\substack{j=1,\\ j\neq i}}^{m}\left\{\operatorname{tr}(\bm{N}^{-1}\bm{N}_{i}\bm{N}^{-1}\bm{N}_{j})\,\sigma_{0_{j}}^{2}\right\} \end{aligned} \tag{3-2-19}\] 在上式中,被估参数 \(m\) 个。将上式写成矩阵形式,即得 \(m\) 类观测值的赫尔默特估计公式: \[\underset{m\times m}{\bm{S}}\,\underset{m\times 1}{\hat{\bm{\theta}}}=\underset{m\times 1}{\bm{W}_{\theta}} \tag{3-2-20}\] 式中 \[\resizebox{\textwidth}{!}{$ \bm{S}=\begin{bmatrix} n_{1}-2\operatorname{tr}(\bm{N}^{-1}\bm{N}_{1})+\operatorname{tr}(\bm{N}^{-1}\bm{N}_{1})^{2} & \operatorname{tr}(\bm{N}^{-1}\bm{N}_{1}\bm{N}^{-1}\bm{N}_{2}) & \cdots & \operatorname{tr}(\bm{N}^{-1}\bm{N}_{1}\bm{N}^{-1}\bm{N}_{m})\\ & n_{2}-2\operatorname{tr}(\bm{N}^{-1}\bm{N}_{2})+\operatorname{tr}(\bm{N}^{-1}\bm{N}_{2})^{2} & \cdots & \operatorname{tr}(\bm{N}^{-1}\bm{N}_{2}\bm{N}^{-1}\bm{N}_{m})\\ &&\ddots&\vdots\\ (\text{对\qquad\qquad\qquad 称})&&& n_{m}-2\operatorname{tr}(\bm{N}^{-1}\bm{N}_{m})+\operatorname{tr}(\bm{N}^{-1}\bm{N}_{m})^{2} \end{bmatrix}$}\] \[\hat{\bm{\theta}}=\begin{bmatrix}\hat{\sigma}_{0_{1}}^{2}&\hat{\sigma}_{0_{2}}^{2}&\cdots&\hat{\sigma}_{0_{m}}^{2}\end{bmatrix}^{T},\qquad \bm{W}_{\theta}=\begin{bmatrix}\bm{V}_{1}^{T}\bm{P}_{1}\bm{V}_{1}&\bm{V}_{2}^{T}\bm{P}_{2}\bm{V}_{2}&\cdots&\bm{V}_{m}^{T}\bm{P}_{m}\bm{V}_{m}\end{bmatrix}^{T}\] 其解为 \[\hat{\bm{\theta}}=\bm{S}^{-1}\bm{W}_{\theta} \tag{3-2-21}\]
方差分量估计的迭代计算步骤如下:
(1) 将观测值按等级或按不同观测来源分类,并进行验前权估计,即确定各类观测值的权的初值 \(\bm{P}_{1},\bm{P}_{2},\cdots,\bm{P}_{m}\);
(2) 进行第一次平差,求得 \(\bm{V}_{i}^{T}\bm{P}_{i}\bm{V}_{i}\);
(3) 按 (3-2-20) 式进行第一次方差分量估计,求得各类观测值单位权方差的第一次估值 \(\hat{\sigma}_{0_{i}}^{2}\),再依下列定权: \[\hat{\bm{P}}_{i}=\frac{c}{\hat{\sigma}_{0_{i}}^{2}\bm{P}_{i}^{-1}} \tag{3-2-22}\] 式中 \(c\) 为任一常数,一般是选 \(\hat{\sigma}_{0_{i}}^{2}\) 中的某一个值;
(4) 反复进行第二项和第三项,即进行:平差—方差分量估计—定权后再平差,直至 \[\hat{\sigma}_{0_{1}}^{2}=\hat{\sigma}_{0_{2}}^{2}=\cdots=\hat{\sigma}_{0_{m}}^{2}\] 为止,或通过必要的检验认为各类单位权方差之比等于 1 为止。
定权公式 (3-2-22) \(\hat{\bm{P}}_{i}=c/\left(\hat{\sigma}_{0_{i}}^{2}\bm{P}_{i}^{-1}\right)\) 的直观含义:由于 \(D(\bm{L}_{i})=\sigma_{0_{i}}^{2}\bm{P}_{i}^{-1}\),希望换用公共单位权方差 \(c\) 表示所有类的方差,于是令 \(\hat{\bm{P}}_{i}^{-1}=\hat{\sigma}_{0_{i}}^{2}\bm{P}_{i}^{-1}/c\)。\(c\) 只决定权重共同的比例,通常取某一类的 \(\hat{\sigma}_{0_{i}}^{2}\) 使该类权保持为 1(例 3-2-1 中 \(P_{\beta}=1\)),便于与先验精度对比。迭代的收敛判据 \(\hat{\sigma}_{0_{1}}^{2}=\cdots=\hat{\sigma}_{0_{m}}^{2}\) 是相对判据——各方差分量收敛到同一公共值即可,该公共值本身仍由 \(c\) 定标。例 3-2-1 的数值可以看到:三次迭代后单位权方差比从 \(1:0.96\) 收敛到 \(1:0.99\),权的变化(\(0.56\to0.59\to0.60\))也相应趋于稳定。
(3-2-22) 式中常数 \(c\) 的选取只改变各权的共同比例因子,不影响平差结果;通常取某一类的 \(\hat{\sigma}_{0_i}^2\),使该类的权保持为 1(见例 3-2-1 中 \(P_\beta=1\))。
单位权方差估值 \(\hat{\sigma}_{0}^{2}=\bm{V}^{T}\bm{P}\bm{V}/r\)、法方程 \(\bm{N}\hat{\bm{X}}=\bm{W}\) 与协因数阵 \(\bm{Q}_{\hat{x}}=(\bm{B}^{T}\bm{P}\bm{B})^{-1}\),均与《最优估计基础》第2章"最小二乘估计"(2.2 节)完全对应;把整体平差解按"平差—定权"循环改写的做法,与该书第2章"递推最小二乘估计"(2.3 节)把整体解改写为递推形式是同一类"化整为零"的思路。本例边角网中用到的单位权方差、法方程解及 \(\bm{N}^{-1}\)、\(\bm{N}_{1}\)、\(\bm{N}_{2}\) 分块关系,都可回查该书第2章相关各节。
例 3-2-1
有边角网如图 3-1,\(A\)、\(B\)、\(C\) 为已知点,\(P_{1}\)、\(P_{2}\) 为待定点,网中观测了 12 个角度和 6 个边长,起算数据和观测值分别列于表 3-1 和表 3-2。根据经验,测角中误差为 \(\pm1.5''\),边长测量中误差为 \(\pm2.0\,\mathrm{cm}\),试按间接平差法进行赫尔默特估计,并求:
(1) 角度、边长观测值的方差估值;
(2) 待定点坐标的平差值及其方差估值。
| 点号 | 坐标 (m) | 坐标方位角 | 边长 | ||
|---|---|---|---|---|---|
| 2-3(lr)4-5(lr)6-6 | \(X\) | \(Y\) | \({}^{\circ}\) | \({}^{\prime}\quad{}^{\prime\prime}\) | (m) |
| \(A\) | 4899.846 | 130.812 | 14 | 00 | 4001.117 |
| \(B\) | 8781.945 | 1099.443 | 123 | 10 | 7734.443 |
| \(C\) | 4548.795 | 7572.622 | |||
| 编号 | 观测角 | 编号 | 观测角 | 编号 | 观测边 | ||||
|---|---|---|---|---|---|---|---|---|---|
| 2-4(lr)6-8(lr)10-10 | \({}^{\circ}\) | \({}^{\prime}\) | \({}^{\prime\prime}\) | \({}^{\circ}\) | \({}^{\prime}\) | \({}^{\prime\prime}\) | (m) | ||
| 1 | 84 | 07 | 38.2 | 7 | 74 | 18 | 16.8 | 13 | 2463.94 |
| 2 | 37 | 46 | 34.9 | 8 | 77 | 27 | 59.1 | 14 | 3414.61 |
| 3 | 58 | 05 | 44.1 | 9 | 28 | 13 | 43.2 | 15 | 5216.16 |
| 4 | 33 | 03 | 03.2 | 10 | 55 | 21 | 09.9 | 16 | 6042.94 |
| 5 | 126 | 01 | 55.7 | 11 | 72 | 22 | 25.8 | 17 | 5085.08 |
| 6 | 20 | 55 | 02.3 | 12 | 52 | 16 | 20.5 | 18 | 5014.99 |
解(1) 根据先验方差(\(m_{\beta}=\pm1.5''\),\(m_{s}=\pm2\,\mathrm{cm}\))进行第一次定权,即 \[P_{\beta}=\frac{\sigma_{0}^{2}}{\sigma_{\beta}^{2}}=1,\ (\text{无量纲})\] \[P_{s}=\frac{\sigma_{0}^{2}}{\sigma_{s}^{2}}=\frac{1.5^{2}}{2^{2}}=0.56\,(('')^{2}/\mathrm{cm}^{2})\]
列误差方程。待定点坐标的近似值取为 \[X_{1}^{0}=5656.89\,\mathrm{m},\ Y_{1}^{0}=2475.56\,\mathrm{m}\] \[X_{2}^{0}=663.90\,\mathrm{m},\ Y_{2}^{0}=2943.91\,\mathrm{m}\] 误差方程的系数和常数项列于表 3-3。
| 序 | \(a\) | \(b\) | \(c\) | \(d\) | 常数项 (\(-l\)) |
|---|---|---|---|---|---|
| 1 | 0.5532 | \(-0.8100\) | 0 | 0 | 0.18 |
| 2 | 0.2434 | 0.5528 | 0 | 0 | \(-0.53\) |
| 3 | \(-0.7966\) | 0.2572 | 0 | 0 | 3.15 |
| 4 | \(-0.2434\) | \(-0.5528\) | 0 | 0 | 0.23 |
| 5 | 0.6298 | 0.6368 | 0 | 0 | \(-2.44\) |
| 6 | \(-0.3864\) | \(-0.0840\) | 0 | 0 | 1.01 |
| 7 | 0.7966 | \(-0.2572\) | \(-0.2244\) | \(-0.3379\) | 2.68 |
| 8 | \(-0.8350\) | \(-0.1523\) | 0.0384 | 0.4095 | \(-4.58\) |
| 9 | 0.0384 | 0.4095 | 0.1860 | \(-0.0716\) | 2.80 |
| 10 | \(-0.0384\) | \(-0.4095\) | 0.2998 | 0.1901 | \(-3.10\) |
| 11 | \(-0.3480\) | 0.3255 | \(-0.0384\) | \(-0.4095\) | 8.04 |
| 12 | 0.3864 | 0.0840 | \(-0.2614\) | 0.2194 | \(-1.14\) |
| 13 | 0.3072 | 0.9516 | 0 | 0 | \(-0.84\) |
| 14 | \(-0.9152\) | 0.4030 | 0 | 0 | 1.54 |
| 15 | 0.2124 | \(-0.9772\) | 0 | 0 | \(-3.93\) |
| 16 | 0 | 0 | \(-0.6429\) | \(-0.7660\) | 2.15 |
| 17 | 0 | 0 | \(-0.8330\) | 0.5532 | \(-12.58\) |
| 18 | 0.9956 | \(-0.0934\) | \(-0.9956\) | 0.0934 | \(-8.21\) |
第一次平差的法方程是 \[\begin{bmatrix} 4.3124 & -0.2886 & -0.8579 & -0.3418\\ & 3.4214 & 0.0229 & -0.2024\\ & & 1.4212 & 0.0592\\ (\text{对\quad 称}) & & & 1.0438 \end{bmatrix} \hat{\bm{x}}+ \begin{bmatrix} -7.5529\\ 6.0311\\ 8.4751\\ -12.3622 \end{bmatrix} =\bm{0}\] 它的解列于表 3-4。
第一次平差后的改正数列于表 3-5。
| 迭代次数 | 1 | 2 | 3 | |
|---|---|---|---|---|
| 坐 | \(\hat{x}_{1}\) | 1.5844 | 1.5719 | 1.5664 |
| 标 | \(\hat{y}_{1}\) | \(-0.8516\) | \(-0.8698\) | \(-0.8756\) |
| 改 | \(\hat{x}_{2}\) | \(-5.5103\) | \(-5.5852\) | \(-5.6086\) |
| 正数 (cm) | \(\hat{y}_{2}\) | 12.5115 | 12.4387 | 12.4159 |
| 序号 | 1 | 2 | 3 | |||
|---|---|---|---|---|---|---|
| 2-3(lr)4-5(lr)6-7 | \(P_{i}\) | \(V_{i}\) | \(P_{i}\) | \(V_{i}\) | \(P_{i}\) | \(V_{i}\) |
| 1 | 1.749 | 1.754 | 1.756 | |||
| 2 | \(-0.614\) | \(-0.628\) | \(-0.633\) | |||
| 3 | 1.665 | 1.674 | 1.677 | |||
| 4 | 0.314 | 0.328 | 0.333 | |||
| 5 | \(-1.981\) | \(-2.004\) | \(-2.011\) | |||
| 6 | 1 | 0.467 | 1 | 0.476 | 1 | 0.478 |
| 7 | 1.174 | 1.206 | 1.217 | |||
| 8 | \(-0.866\) | \(-0.881\) | \(-0.886\) | |||
| 9 | 0.592 | 0.575 | 0.569 | |||
| 10 | \(-2.086\) | \(-2.114\) | \(-2.123\) | |||
| 11 | 2.298 | 2.331 | 2.341 | |||
| 12 | 3.588 | 3.583 | 3.582 | |||
| 13 | \(-1.162\) | \(-1.185\) | \(-1.192\) | |||
| 14 | \(-0.258\) | \(-0.249\) | \(-0.246\) | |||
| 15 | 0.56 | \(-2.760\) | 0.59 | \(-2.746\) | 0.60 | \(-2.742\) |
| 16 | \(-3.891\) | \(-3.787\) | \(-3.754\) | |||
| 17 | \(-1.069\) | \(-1.040\) | \(-1.040\) | |||
| 18 | 0.107 | 0.159 | 0.175 | |||
按 (3-2-16) 式进行赫尔默特估计。其中 \[\bm{N}^{-1}=\begin{bmatrix} 0.2724 & 0.0269 & 0.1604 & 0.0853\\ & 0.2984 & 0.0087 & 0.0662\\ & & 0.8000 & 0.0088\\ (\text{对\qquad\qquad 称}) & & & 0.9983 \end{bmatrix}\] \[\bm{N}_{1}=\begin{bmatrix} 3.2102 & -0.0774 & -0.3028 & -0.3939\\ & 2.2837 & -0.0292 & -0.1975\\ & & 0.2461 & 0.0936\\ (\text{对\qquad\qquad 称}) & & & 0.5390 \end{bmatrix}\] \[\bm{N}_{2}=\begin{bmatrix} 1.1022 & -0.2112 & -0.5551 & 0.0521\\ & 1.1377 & 0.0521 & -0.0049\\ & & 1.1751 & -0.0344\\ (\text{对\qquad\qquad 称}) & & & 0.5048 \end{bmatrix}\]
方程 (3-2-16) 式的系数和常数项及其解算结果列于表 3-6。
| 迭代次数 | 系数矩阵 \(\bm{S}\) | \(\bm{V}_{i}^{T}\bm{P}_{i}\bm{V}_{i}\) | \(\hat{\sigma}_{0_{\beta}}^{2}\,(''^{2})\) | \(\hat{\sigma}_{0_{s}}^{2}\) | \(\hat{\sigma}_{0_{\beta}}^{2}:\hat{\sigma}_{0_{s}}^{2}\) | |
|---|---|---|---|---|---|---|
| 1 | \(\begin{bmatrix}9.1388 & 0.7640\\ 0.7640 & 3.3333\end{bmatrix}\) | 35.4339 | 14.1840 | 3.5904 | 3.4323 | \(1:0.96\) |
| 2 | \(\begin{bmatrix}9.1756 & 0.7670\\ 0.7670 & 3.2903\end{bmatrix}\) | 35.9250 | 14.4359 | 3.6190 | 3.5438 | \(1:0.98\) |
| 3 | \(\begin{bmatrix}9.1875 & 0.7681\\ 0.7681 & 3.2763\end{bmatrix}\) | 36.0880 | 14.5228 | 3.6285 | 3.5820 | \(1:0.99\) |
根据第一次估算出的两类观测值的单位权方差 \(\hat{\sigma}_{0_{\beta}}^{2}\) 和 \(\hat{\sigma}_{0_{s}}^{2}\),计算角度和边长观测值的方差估值,其计算公式为 \[\hat{\sigma}_{\beta}^{2}=\hat{\sigma}_{0_{\beta}}^{2}P_{\beta}^{-1},\qquad \hat{\sigma}_{s}^{2}=\hat{\sigma}_{0_{s}}^{2}P_{s}^{-1} \tag{3-2-23}\] 即 \[\hat{\sigma}_{\beta}^{2}=\hat{\sigma}_{0_{\beta}}^{2}=3.5904\,(''^{2})\] \[\hat{\sigma}_{s}^{2}=3.4323\times(0.56)^{-1}=6.1291\,(\mathrm{cm}^{2})\]
(2) 根据第一次估算出的方差 \(\hat{\sigma}_{\beta}^{2}\) 和 \(\hat{\sigma}_{s}^{2}\) 进行第二次定权。则有 \[P_{\beta}=\frac{\sigma_{0}^{2}}{\hat{\sigma}_{\beta}^{2}}=1\] \[P_{s}=\frac{\sigma_{0}^{2}}{\hat{\sigma}_{s}^{2}}=\frac{3.5904}{6.1291}=0.59\] 第二次平差的法方程是 \[\begin{bmatrix} 4.3715 & -0.2999 & -0.8876 & -0.3390\\ & 3.4823 & 0.0257 & -0.2026\\ & & 1.4842 & 0.0574\\ (\text{对\quad 称}) & & & 1.0709 \end{bmatrix} \hat{\bm{x}}+ \begin{bmatrix} -7.8732\\ 6.1640\\ 8.9932\\ -12.6434 \end{bmatrix} =\bm{0}\] 它的解仍列于表 3-4 中。
第二次平差后的改正数仍列于表 3-5 中。
按 (3-2-16) 式进行第二次赫尔默特估计,其中 \(\bm{N}_{1}\) 矩阵与第一次估计时的 \(\bm{N}_{1}\) 相同,而 \[\bm{N}^{-1}=\begin{bmatrix} 0.2688 & 0.0267 & 0.1571 & 0.0817\\ & 0.2931 & 0.0085 & 0.0635\\ & & 0.7672 & 0.0102\\ (\text{对\qquad\qquad 称}) & & & 0.9711 \end{bmatrix}\] \[\bm{N}_{2}=\begin{bmatrix} 1.1613 & -0.2225 & -0.5848 & 0.0549\\ & 1.1986 & 0.0549 & -0.0051\\ & & 1.2381 & -0.0362\\ (\text{对\qquad\qquad 称}) & & & 0.5319 \end{bmatrix}\]
第二次的估计结果仍列于表 3-6。
根据第二次估算出的 \(\hat{\sigma}_{0_{\beta}}^{2}\),\(\hat{\sigma}_{0_{s}}^{2}\),\(\hat{\sigma}_{\beta}^{2}\) 和 \(\hat{\sigma}_{s}^{2}\): \[\hat{\sigma}_{\beta}^{2}=\hat{\sigma}_{0_{\beta}}^{2}P_{\beta}^{-1}=3.6190\,(''^{2})\] \[\hat{\sigma}_{s}^{2}=\hat{\sigma}_{0_{s}}^{2}P_{s}^{-1}=3.5438\times(0.59)^{-1}=6.0064\,(\mathrm{cm}^{2})\]
(3) 重复 (2) 的计算程序,进行第三次定权、第三次平差及第三次赫尔默特估计。其中: \[P_{\beta}=1,\qquad P_{s}=0.60\] 法方程是: \[\begin{bmatrix} 4.3912 & -0.3036 & -0.8975 & -0.3381\\ & 3.5027 & 0.0266 & -0.2027\\ & & 1.5052 & 0.0568\\ (\text{对\quad 称}) & & & 1.0799 \end{bmatrix} \hat{\bm{x}}+ \begin{bmatrix} -7.9800\\ 6.2083\\ 9.1659\\ -12.7372 \end{bmatrix} =\bm{0}\] \(\bm{N}_{1}\) 矩阵仍与第一次估计时的 \(\bm{N}_{1}\) 相同,而 \(\bm{N}^{-1}\) 和 \(\bm{N}_{2}\) 分别为 \[\bm{N}^{-1}=\begin{bmatrix} 0.2677 & 0.0267 & 0.1561 & 0.0806\\ & 0.2914 & 0.0084 & 0.0626\\ & & 0.7569 & 0.0106\\ (\text{对\qquad\qquad 称}) & & & 0.9624 \end{bmatrix}\] \[\bm{N}_{2}=\begin{bmatrix} 1.1810 & -0.2262 & -0.5947 & 0.0558\\ & 1.2190 & 0.0558 & -0.0052\\ & & 1.2591 & -0.0368\\ (\text{对\qquad\qquad 称}) & & & 0.5409 \end{bmatrix}\]
由表 3-6 可知,经过三次迭代计算后,角度和边长的单位权方差已趋于一致,因此,迭代计算结束。现将三次赫尔默特估计法的主要数据汇集于表 3-7 中。
| 迭代次数 | 观测值的精度估值 | 单位权方差之比 | |
|---|---|---|---|
| 2-3(lr)4-4 | \(\hat{\sigma}_{\beta}\,('')\) | \(\hat{\sigma}_{s}\,(\mathrm{cm})\) | \(\hat{\sigma}_{0_{\beta}}^{2}:\hat{\sigma}_{0_{s}}^{2}\) |
| 0 | 1.50 | 2.00 | (根据经验估计的精度) |
| 1 | 1.89 | 2.48 | \(1:0.96\) |
| 2 | 1.90 | 2.45 | \(1:0.98\) |
| 3 | 1.91 | 2.44 | \(1:0.99\) |
(4) 待定点的最后坐标及其精度指标的计算。坐标计算式为 \[\hat{\bm{X}}=\bm{X}_{0}+\hat{\bm{x}}\] 式中的 \(\hat{\bm{x}}\) 为第三次平差的结果。精度计算式为 \[m_{x_{i}}=m_{0}\sqrt{Q_{x_{i}}}\] \[m_{y_{i}}=m_{0}\sqrt{Q_{y_{i}}}\] 式中 \(m_{0}\) 与 \(Q_{i}\) 均采用第三次平差的结果。其计算数值列于表 3-8 中。
| 点号 | 最后坐标 | 精度指标 | ||
|---|---|---|---|---|
| 2-3(lr)4-5 | \(\hat{X}\,(\mathrm{m})\) | \(\hat{Y}\,(\mathrm{m})\) | \(m_{x}\) (cm) | \(m_{y}\) (cm) |
| \(P_{1}\) | 5656.906 | 2475.551 | 0.98 | 1.03 |
| \(P_{2}\) | 663.844 | 2944.034 | 1.65 | 1.87 |
其他经典平差方法的方差分量估计
(1) 条件平差时的方差分量估计。
条件平差的函数模型: \[\underset{r\times n}{\bm{A}}\,\underset{n\times 1}{\bm{L}}+\underset{r\times n}{\bm{A}}\,\underset{n\times 1}{\bm{\Delta}}+\underset{r\times 1}{\bm{A}_{0}}=\bm{0} \tag{3-2-24}\] 随机模型: \[\begin{aligned} &E(\bm{L})=\widetilde{\bm{L}},\ E(\bm{\Delta})=\bm{0}\\ &D(\bm{L})=\sigma_{0}^{2}\bm{P}^{-1},\ D(\bm{\Delta})=D(\bm{L})=\sigma_{0}^{2}\bm{P}^{-1} \end{aligned} \tag{3-2-25}\] 条件方程: \[\bm{A}\bm{V}-\bm{f}=\bm{0} \tag{3-2-26}\] 式中 \[-\bm{f}=\bm{A}\bm{L}+\bm{A}_{0}\] 法方程: \[\bm{N}\bm{K}-\bm{f}=\bm{0} \tag{3-2-27}\] 式中 \[\bm{N}=\bm{A}\bm{P}^{-1}\bm{A}^{T}\] 法方程的解: \[\bm{K}=\bm{N}^{-1}\bm{f} \tag{3-2-28}\] 改正数方程: \[\bm{V}=\bm{P}^{-1}\bm{A}^{T}\bm{K}\]
设 \(\bm{L}\) 中仅含有两类相互独立的观测值 \(\underset{n_{1}\times 1}{\bm{L}_{1}}\) 和 \(\underset{n_{2}\times 1}{\bm{L}_{2}}\),且它们的权阵分别为 \(\bm{P}_{1}\) 和 \(\bm{P}_{2}\)(\(\bm{P}_{12}=\bm{0}\)),则条件方程为 \[\bm{A}_{1}\bm{V}_{1}+\bm{A}_{2}\bm{V}_{2}-\bm{f}=\bm{0} \tag{3-2-29}\] 式中 \[-\bm{f}=\bm{A}_{1}\bm{L}_{1}+\bm{A}_{2}\bm{L}_{2}+\bm{A}_{0}\] 且下列关系式成立: \[\begin{gathered} \bm{L}=\begin{bmatrix}\bm{L}_{1}\\ \bm{L}_{2}\end{bmatrix},\quad \bm{V}=\begin{bmatrix}\bm{V}_{1}\\ \bm{V}_{2}\end{bmatrix},\quad \bm{A}=\begin{bmatrix}\bm{A}_{1}&\bm{A}_{2}\end{bmatrix},\quad \bm{P}=\begin{bmatrix}\bm{P}_{1}&\bm{0}\\ \bm{0}&\bm{P}_{2}\end{bmatrix}\\ \bm{N}=\bm{A}\bm{P}^{-1}\bm{A}^{T}=\bm{A}_{1}\bm{P}_{1}^{-1}\bm{A}_{1}^{T}+\bm{A}_{2}\bm{P}_{2}^{-1}\bm{A}_{2}^{T}=\bm{N}_{1}+\bm{N}_{2} \end{gathered}\]
若第一次平差时给定的两类观测值的权 \(\bm{P}_{1}\) 和 \(\bm{P}_{2}\) 不恰当,则可设 \[D(\bm{L}_{1})=\sigma_{0_{1}}^{2}\bm{P}_{1}^{-1},\qquad D(\bm{L}_{2})=\sigma_{0_{2}}^{2}\bm{P}_{2}^{-1}\] 现在的任务是:利用 \(\bm{V}_{1}^{T}\bm{P}_{1}\bm{V}_{1}\) 和 \(\bm{V}_{2}^{T}\bm{P}_{2}\bm{V}_{2}\) 估计两类观测值的单位权方差因子 \(\sigma_{0_{1}}^{2}\) 和 \(\sigma_{0_{2}}^{2}\)。
因 \[\bm{V}_{1}=\bm{P}_{1}^{-1}\bm{A}_{1}^{T}\bm{K}=\bm{P}_{1}^{-1}\bm{A}_{1}^{T}\bm{N}^{-1}\bm{f} =-\bm{P}_{1}^{-1}\bm{A}_{1}^{T}\bm{N}^{-1}(\bm{A}_{1}\bm{L}_{1}+\bm{A}_{2}\bm{L}_{2}+\bm{A}_{0})\] 所以有 \[E(\bm{V}_{1})=\bm{0}\] \[\begin{aligned} D(\bm{V}_{1})&=\bm{P}_{1}^{-1}\bm{A}_{1}^{T}\bm{N}^{-1} \begin{bmatrix}\bm{A}_{1}&\bm{A}_{2}\end{bmatrix} \begin{bmatrix}\sigma_{0_{1}}^{2}\bm{P}_{1}^{-1}&\bm{0}\\ \bm{0}&\sigma_{0_{2}}^{2}\bm{P}_{2}^{-1}\end{bmatrix} \begin{bmatrix}\bm{A}_{1}^{T}\\ \bm{A}_{2}^{T}\end{bmatrix} \bm{N}^{-1}\bm{A}_{1}\bm{P}_{1}^{-1}\\ &=\sigma_{0_{1}}^{2}\bm{P}_{1}^{-1}\bm{A}_{1}^{T}\bm{N}^{-1}\bm{N}_{1}\bm{N}^{-1}\bm{A}_{1}\bm{P}_{1}^{-1} +\sigma_{0_{2}}^{2}\bm{P}_{1}^{-1}\bm{A}_{1}^{T}\bm{N}^{-1}\bm{N}_{2}\bm{N}^{-1}\bm{A}_{1}\bm{P}_{1}^{-1} \end{aligned}\] 将上两式代入 (3-2-8) 式,则得 \[\begin{aligned} E(\bm{V}_{1}^{T}\bm{P}_{1}\bm{V}_{1})={}&\operatorname{tr}(\bm{P}_{1}D(\bm{V}_{1}))\\ ={}&\operatorname{tr}(\bm{P}_{1}\bm{P}_{1}^{-1}\bm{A}_{1}^{T}\bm{N}^{-1}\bm{N}_{1}\bm{N}^{-1}\bm{A}_{1}\bm{P}_{1}^{-1})\,\sigma_{0_{1}}^{2}+{}\\ &\operatorname{tr}(\bm{P}_{1}\bm{P}_{1}^{-1}\bm{A}_{1}^{T}\bm{N}^{-1}\bm{N}_{2}\bm{N}^{-1}\bm{A}_{1}\bm{P}_{1}^{-1})\,\sigma_{0_{2}}^{2}\\ ={}&\operatorname{tr}(\bm{N}^{-1}\bm{N}_{1}\bm{N}^{-1}\bm{N}_{1})\,\sigma_{0_{1}}^{2} +\operatorname{tr}(\bm{N}^{-1}\bm{N}_{1}\bm{N}^{-1}\bm{N}_{2})\,\sigma_{0_{2}}^{2} \end{aligned} \tag{3-2-30}\] 同理可得 \[E(\bm{V}_{2}^{T}\bm{P}_{2}\bm{V}_{2}) =\operatorname{tr}(\bm{N}^{-1}\bm{N}_{1}\bm{N}^{-1}\bm{N}_{2})\,\sigma_{0_{1}}^{2} +\operatorname{tr}(\bm{N}^{-1}\bm{N}_{2})^{2}\,\sigma_{0_{2}}^{2} \tag{3-2-31}\] 由 (3-2-30) 和 (3-2-31) 两式,即得估算两类观测值单位权方差因子的公式为 \[\underset{2\times 2}{\bm{S}_{A}}\,\underset{2\times 1}{\hat{\bm{\theta}}}=\underset{2\times 1}{\bm{W}_{\theta}} \tag{3-2-32}\] 式中 \[\bm{S}_{A}=\begin{bmatrix} \operatorname{tr}(\bm{N}^{-1}\bm{N}_{1})^{2} & \operatorname{tr}(\bm{N}^{-1}\bm{N}_{1}\bm{N}^{-1}\bm{N}_{2})\\ (\text{对\quad 称}) & \operatorname{tr}(\bm{N}^{-1}\bm{N}_{2})^{2} \end{bmatrix} \tag{3-2-33}\] \[\hat{\bm{\theta}}=\begin{bmatrix}\hat{\sigma}_{0_{1}}^{2}&\hat{\sigma}_{0_{2}}^{2}\end{bmatrix},\qquad \bm{W}_{\theta}=\begin{bmatrix}\bm{V}_{1}^{T}\bm{P}_{1}\bm{V}_{1}&\bm{V}_{2}^{T}\bm{P}_{2}\bm{V}_{2}\end{bmatrix} \tag{3-2-34}\]
(2) 附有参数的条件平差时的方差分量估计。
在上述条件平差的方差分量估计公式 (3-2-33) 中,只要将 \(\bm{N}^{-1}\) 改写联系数 \(\bm{K}\) 的协因数阵 \(\bm{Q}_{K}\),即得附有参数的条件平差时的方差分量估计。即其估计公式为 (3-2-32) 式、(3-2-34) 式及 \[\bm{S}=\begin{bmatrix} \operatorname{tr}(\bm{Q}_{K}\bm{N}_{1})^{2} & \operatorname{tr}(\bm{Q}_{K}\bm{N}_{1}\bm{Q}_{K}\bm{N}_{2})\\ (\text{对称}) & \operatorname{tr}(\bm{Q}_{K}\bm{N}_{2})^{2} \end{bmatrix} \tag{3-2-35}\]
(3) 附有条件的间接平差时的方差分量估计。
由 (3-2-16) 式改写为 \[\bm{S}=\begin{bmatrix} n_{1}-2\operatorname{tr}(\bm{Q}_{\hat{x}}\bm{N}_{1})+\operatorname{tr}(\bm{Q}_{\hat{x}}\bm{N}_{1})^{2} & \operatorname{tr}(\bm{Q}_{\hat{x}}\bm{N}_{1}\bm{Q}_{\hat{x}}\bm{N}_{2})\\ (\text{对称}) & n_{2}-2\operatorname{tr}(\bm{Q}_{\hat{x}}\bm{N}_{2})+\operatorname{tr}(\bm{Q}_{\hat{x}}\bm{N}_{2})^{2} \end{bmatrix} \tag{3-2-36}\] 式中 \(\bm{Q}_{\hat{x}}\) 为未知参数估值 \(\hat{\bm{x}}\) 的协因数阵。(3-2-32)、(3-2-34) 式及 (3-2-36) 式即为其方差分量估计公式。
秩亏自由网平差时方差分量估计公式
由间接平差公式可得 \[\begin{gathered} \underset{2\times 2}{\bm{S}}\,\underset{2\times 1}{\hat{\bm{\theta}}}=\underset{2\times 1}{\bm{W}_{\theta}}\\ \bm{S}=\begin{bmatrix} n_{1}-2\operatorname{tr}(\bm{N}^{+}\bm{N}_{1})+\operatorname{tr}(\bm{N}^{+}\bm{N}_{1})^{2} & \operatorname{tr}(\bm{N}^{+}\bm{N}_{1}\bm{N}^{+}\bm{N}_{2})\\ (\text{对称}) & n_{2}-2\operatorname{tr}(\bm{N}^{+}\bm{N}_{2})+\operatorname{tr}(\bm{N}^{+}\bm{N}_{2})^{2} \end{bmatrix}\\ \hat{\bm{\theta}}=\begin{bmatrix}\hat{\sigma}_{0_{1}}^{2}&\hat{\sigma}_{0_{2}}^{2}\end{bmatrix},\qquad \bm{W}_{\theta}=\begin{bmatrix}\bm{V}_{1}^{T}\bm{P}_{1}\bm{V}_{1}&\bm{V}_{2}^{T}\bm{P}_{2}\bm{V}_{2}\end{bmatrix} \end{gathered} \tag{3-2-37}\]
滤波配置的方差分量估计公式
在滤波和配置问题中,如果按 (2-5-53) 式或按 (2-5-52) 式求得的 \(\hat{\sigma}_{0}^{2}\) 不等于 1,而 \(\hat{\sigma}_{0}^{2}-1\) 是一个较大的数值,则说明所给定的 \(D_{\Delta}\) 和 \(D_{X}\) 不合适,应当重新求定。为此,也可以按赫尔默特方法来求定 \(\bm{Q}_{\Delta}\) 和 \(\bm{Q}_{X}\) 所对应的单位权方差估值 \(\hat{\sigma}_{\Delta}^{2}\) 和 \(\hat{\sigma}_{X}^{2}\),则有 \[\left.\begin{aligned} \bm{V}^{T}\bm{P}_{\Delta}\bm{V}={}&\{n-2\operatorname{tr}(\bm{N}^{-1}\bm{N}_{1}) +\operatorname{tr}(\bm{N}^{-1}\bm{N}_{1}\bm{N}^{-1}\bm{N}_{1})\}\,\hat{\sigma}_{\Delta}^{2}+{}\\ &\operatorname{tr}(\bm{N}^{-1}\bm{N}_{1}\bm{N}^{-1}\bm{N}_{2})\,\hat{\sigma}_{X}^{2}\\ \bm{V}_{X}^{T}\bm{P}_{X}\bm{V}_{X}={}&\operatorname{tr}(\bm{N}^{-1}\bm{N}_{1}\bm{N}^{-1}\bm{N}_{2})\,\hat{\sigma}_{\Delta}^{2} +\{t_{1}-2\operatorname{tr}(\bm{N}^{-1}\bm{N}_{2})+{}\\ &\operatorname{tr}(\bm{N}^{-1}\bm{N}_{2}\bm{N}^{-1}\bm{N}_{2})\}\,\hat{\sigma}_{X}^{2} \end{aligned}\right\} \tag{3-2-38}\] 对于滤波,(3-2-28) 式中的 \(\bm{N}\)、\(\bm{N}_{1}\)、\(\bm{N}_{2}\) 为 \[\left.\begin{aligned} \bm{N}&=\bm{B}^{T}\bm{P}_{\Delta}\bm{B}+\bm{P}_{\bar{X}}\\ \bm{N}_{1}&=\bm{B}^{T}\bm{P}_{\Delta}\bm{B}\\ \bm{N}_{2}&=\bm{P}_{\bar{X}} \end{aligned}\right\}\] 对于配置,有 \[\left.\begin{aligned} \bm{N}&=\begin{bmatrix} \bm{B}^{T}\bm{P}_{\Delta}\bm{B}+\bm{P}_{\bar{X}} & \bm{B}^{T}\bm{P}_{\Delta}\bm{G}\\ \bm{G}^{T}\bm{P}_{\Delta}\bm{B} & \bm{G}^{T}\bm{P}_{\Delta}\bm{G} \end{bmatrix}\\ \bm{N}_{1}&=\begin{bmatrix} \bm{B}^{T}\bm{P}_{\Delta}\bm{B} & \bm{B}^{T}\bm{P}_{\Delta}\bm{G}\\ \bm{G}^{T}\bm{P}_{\Delta}\bm{B} & \bm{G}^{T}\bm{P}_{\Delta}\bm{G} \end{bmatrix}\\ \bm{N}_{2}&=\begin{bmatrix} \bm{P}_{\bar{X}} & \bm{0}\\ \bm{0} & \bm{0} \end{bmatrix} \end{aligned}\right\} \tag{3-2-39}\]
上述法方程中未出现推估的未测点信号 \(\bm{X}^{\prime}\),这是因为方差分量估计与推估无关。
赫尔默特方差分量估计的简化公式
上面已经导出了方差分量的严密估算公式,从统计的角度讲,它们均具有无偏的良好特性。但对于一个较大型的控制网而言,计算工作量很大,这主要表现在几个矩阵的乘法运算和求迹运算。即使有了大型的电子计算机,要想同时存储逆矩阵 \(\bm{N}^{-1}\) 及所有的子块矩阵 \(\bm{N}_{i}\,(i=1,2,\cdots,m)\),就难以实现,更不必说进行矩阵运算了。为此,需要寻求近似的、快速的计算方法。
在严密公式 (3-2-20) 中,略去求迹部分,则有 \[\bm{V}_{i}^{T}\bm{P}_{i}\bm{V}_{i}=n_{i}\hat{\sigma}_{0_{i}}^{2}\] 即各类观测值的单位权方差的估值为 \[\hat{\sigma}_{0_{i}}^{2}=\frac{\bm{V}_{i}^{T}\bm{P}_{i}\bm{V}_{i}}{n_{i}} \tag{3-2-40}\] 对 (3-2-40) 式求和便知,它是有偏估计,上式为赫尔默特近似公式。
在严密公式 (3-2-19) 中,假定 \(\sigma_{0_{1}}^{2}=\sigma_{0_{2}}^{2}=\cdots=\sigma_{0_{i}}^{2}=\cdots=\sigma_{0_{m}}^{2}=\sigma_{0_{i}}^{2}\,(\neq\sigma_{0}^{2})\),则得 \[\begin{aligned} E(\bm{V}_{i}^{T}\bm{P}_{i}\bm{V}_{i})&=\{n_{i}-2\operatorname{tr}(\bm{N}^{-1}\bm{N}_{i}) +\operatorname{tr}(\bm{N}^{-1}\bm{N}_{i}\bm{N}^{-1}\sum_{j=1}^{m}\bm{N}_{j})\}\,\sigma_{0_{i}}^{2}\\ &=\left(n_{i}-2\operatorname{tr}(\bm{N}^{-1}\bm{N}_{i})+\operatorname{tr}(\bm{N}^{-1}\bm{N}_{i}\bm{N}^{-1}\bm{N})\right)\sigma_{0_{i}}^{2}\\ &=\left(n_{i}-\operatorname{tr}(\bm{N}^{-1}\bm{N}_{i})\right)\sigma_{0_{i}}^{2} \end{aligned}\] 由此得简化公式为 \[\bm{V}_{i}^{T}\bm{P}_{i}\bm{V}_{i}=\left(n_{i}-\operatorname{tr}(\bm{N}^{-1}\bm{N}_{i})\right)\hat{\sigma}_{0_{i}}^{2} \tag{3-2-41}\] 上式也可通过下述途径导出:
由间接平差的基本公式可知,改正数 \(\bm{V}\) 的协因数阵为 \[\bm{Q}_{V}=\bm{Q}-\bm{B}\bm{N}^{-1}\bm{B}^{T} \tag{3-2-42}\] 即有 \[\bm{Q}_{V_{i}}=\bm{Q}_{i}-\bm{B}_{i}\bm{N}^{-1}\bm{B}_{i}^{T} \tag{3-2-43}\] \[D(\bm{V}_{i})=\left(\bm{Q}_{i}-\bm{B}_{i}\bm{N}^{-1}\bm{B}_{i}^{T}\right)\sigma_{0_{i}}^{2}\] 所以 \[E(\bm{V}_{i}^{T}\bm{P}_{i}\bm{V}_{i}) =\operatorname{tr}(\bm{P}_{i}\bm{Q}_{i}-\bm{P}_{i}\bm{B}_{i}\bm{N}^{-1}\bm{B}_{i}^{T})\,\sigma_{0_{i}}^{2} =\left(n_{i}-\operatorname{tr}(\bm{N}^{-1}\bm{N}_{i})\right)\sigma_{0_{i}}^{2} \tag{3-2-44}\] 上式与 (3-2-41) 式完全一致,因此,由上式得到的简化公式应为 \[\hat{\sigma}_{0_{i}}^{2}=\frac{\bm{V}_{i}^{T}\bm{P}_{i}\bm{V}_{i}}{n_{i}-\operatorname{tr}(\bm{N}^{-1}\bm{N}_{i})} \tag{3-2-45}\]
(3-2-45) 是在"各类单位权方差相等"这一假定下从 (3-2-19) 导出的简化公式:当 \(\sigma_{0_{1}}^{2}=\cdots=\sigma_{0_{m}}^{2}\) 时,(3-2-19) 中所有 \(\operatorname{tr}(\bm{N}^{-1}\bm{N}_{i}\bm{N}^{-1}\bm{N}_{j})\) 型交叉项与主项合并为 \(\left(n_{i}-\operatorname{tr}(\bm{N}^{-1}\bm{N}_{i})\right)\sigma_{0_{i}}^{2}\),交叉耦合项才恰好消失。若初权严重不合理、各 \(\sigma_{0_{i}}^{2}\) 相差很大,简化公式仍有偏差;实践中常用严密公式 (3-2-20) 先迭代若干步、待权比稳定后再换简化公式加速。注意 (3-2-40)(分母 \(n_{i}\))与 (3-2-45)(分母 \(r_{i}\))在同一轮迭代中会给出不同结果,收敛到的权比也不同。
由 (3-2-43) 式不难看出,上式中的 \[n_{i}-\operatorname{tr}(\bm{N}^{-1}\bm{N}_{i})=\operatorname{tr}(\bm{P}\bm{Q}_{V})_{i}\] 因此,一般称 \((n_{i}-\operatorname{tr}(\bm{N}^{-1}\bm{N}_{i}))\) 为第 \(i\) 类观测值 \(\bm{L}_{i}\) 的多余观测分量,若令 \[\operatorname{tr}(\bm{P}\bm{Q}_{V})_{i}=n_{i}-\operatorname{tr}(\bm{N}^{-1}\bm{N}_{i})=r_{i} \tag{3-2-46}\] 则 (3-2-45) 式可写成: \[\hat{\sigma}_{0_{i}}^{2}=\frac{\bm{V}_{i}^{T}\bm{P}_{i}\bm{V}_{i}}{r_{i}} \tag{3-2-47}\]
多余观测分量 \(r_{i}=n_{i}-\operatorname{tr}(\bm{N}^{-1}\bm{N}_{i})\) 的直觉:它是第 \(i\) 类观测中"真正冗余"的个数。对全部类求和,利用 \(\sum_{i}\bm{N}_{i}=\bm{N}\) 与迹的性质得 \[\sum_{i=1}^{m}r_{i}=n-\operatorname{tr}\!\left(\bm{N}^{-1}\textstyle\sum_{i}\bm{N}_{i}\right)=n-\operatorname{tr}(\bm{N}^{-1}\bm{N})=n-t=r,\] 恰好回到总的多余观测数。\(r_{i}\) 大,说明第 \(i\) 类观测多、对参数的约束弱(\(\operatorname{tr}(\bm{N}^{-1}\bm{N}_{i})\) 小),其残差平方和越"可信"。这也是 3-5 节精度评定中"各类观测的 \(r_{i}\) 大致相等时,方差分量估值的可靠度才一致"(例 3-5-3)的由来。
一般情况下,在大型控制网平差中,法方程系数矩阵 \(\bm{N}\) 的非主对角元相对较小。因此,略去法方程系数阵 \(\bm{N}\) 和子块矩阵 \(\bm{N}_{i}\) 中的非主对角元仍利用严密公式进行方差估计也可得到较好的近似。此法由周江文给出(文献 [9])。
令略去非主对角元后的矩阵为 \[\begin{aligned} \bm{N}&=\bm{B}^{T}\bm{P}\bm{B}=\operatorname{diag}\begin{bmatrix}\bm{N}_{a}&\bm{N}_{b}&\cdots&\bm{N}_{t}\end{bmatrix}\\ \bm{N}_{i}&=\bm{B}_{i}^{T}\bm{P}_{i}\bm{B}_{i}=\operatorname{diag}\begin{bmatrix}\bm{N}_{a_{i}}&\bm{N}_{b_{i}}&\cdots&\bm{N}_{t_{i}}\end{bmatrix} \end{aligned}\] 则有 \[\begin{aligned} \operatorname{tr}(\bm{N}^{-1}\bm{N}_{i})&=\left[\frac{\bm{N}_{k_{i}}}{\bm{N}_{k}}\right]\quad(k=a,b,\cdots,t)\\ \operatorname{tr}(\bm{N}^{-1}\bm{N}_{i})^{2}&=\left[\frac{\bm{N}_{k_{i}}^{2}}{\bm{N}_{k}^{2}}\right]\\ \operatorname{tr}(\bm{N}^{-1}\bm{N}_{i}\bm{N}^{-1}\bm{N}_{j})&=\left[\frac{\bm{N}_{k_{i}}\bm{N}_{k_{j}}}{\bm{N}_{k}^{2}}\right] \end{aligned}\]
方差-协方差分量估计
为推导公式简便,本节将首先导出两类相关观测值 \(\bm{L}\) 的方差-协方差估计公式,然后再将其扩展到 \(m\) 类相关观测值的一般情况。
间接平差的基本公式仍为 (3-2-1) (3-2-4),其中相关观测值的相关权矩阵设为 \(\bm{P}_{L}\),再设在 \(\bm{L}\) 中包含有两类(或精度不等)相关观测值 \(\bm{L}_{1}\) 和 \(\bm{L}_{2}\),一般说来,第一次平差时给定两类观测值的权阵是不恰当的,也可以说它们所对应的单位权方差与单位权协方差不相等,令其分别为 \(\sigma_{0_{1}}^{2}\),\(\sigma_{0_{12}}^{2}\),\(\sigma_{0_{2}}^{2}\),则观测值 \(\bm{L}\) 的方差阵为 \[D(\bm{L})=\begin{bmatrix} D(\bm{L}_{1}) & D(\bm{L}_{1},\bm{L}_{2})\\ D(\bm{L}_{2},\bm{L}_{1}) & D(\bm{L}_{2}) \end{bmatrix} =\begin{bmatrix} \sigma_{0_{1}}^{2}\bm{Q}_{11} & \sigma_{0_{12}}^{2}\bm{Q}_{12}\\ \sigma_{0_{12}}^{2}\bm{Q}_{12}^{T} & \sigma_{0_{2}}^{2}\bm{Q}_{22} \end{bmatrix} \tag{3-3-1}\] 式中,自协因数阵和互协因数阵 \(\bm{Q}_{11}\),\(\bm{Q}_{12}\),\(\bm{Q}_{22}\) 是已知的,并组成对称正定阵,而 \(\sigma_{0_{1}}^{2}\)、\(\sigma_{0_{12}}^{2}\) 和 \(\sigma_{0_{2}}^{2}\) 是待定量。所谓方差-协方差估计,就是通过平差得到的残差平方和 \(\bm{V}_{1}^{T}\bm{P}_{11}\bm{V}_{1}\),\(\bm{V}_{1}^{T}\bm{P}_{12}\bm{V}_{2}+\bm{V}_{2}^{T}\bm{P}_{12}^{T}\bm{V}_{1}\) 和 \(\bm{V}_{2}^{T}\bm{P}_{22}\bm{V}_{2}\) 去估计 \(\sigma_{0_{1}}^{2}\),\(\sigma_{0_{12}}^{2}\) 和 \(\sigma_{0_{2}}^{2}\)。为叙述方便,令待估量用向量 \(\bm{\theta}\) 表示出,即有 \[\bm{\theta}=\begin{bmatrix}\sigma_{0_{1}}^{2}&\sigma_{0_{12}}^{2}&\sigma_{0_{2}}^{2}\end{bmatrix}^{T}\] 为与向量的维数号一致,上式也可表示为 \[\bm{\theta}=\begin{bmatrix}\sigma_{0_{1}}^{2}&\sigma_{0_{12}}^{2}&\sigma_{0_{2}}^{2}\end{bmatrix}^{T} =\begin{bmatrix}\theta_{1}&\theta_{2}&\theta_{3}\end{bmatrix}^{T} \tag{3-3-2}\] 下面推导向量 \(\bm{\theta}\) 的估求公式。
为平差时定权的需要,应取 \(\bm{\theta}\) 的初值,一般取 \[\bm{\theta}(0)=\begin{bmatrix}1&1&1\end{bmatrix}^{T} \tag{3-3-3}\] 此时观测值 \(\bm{L}\) 的方差阵为 \[D_{0}(\bm{L})=\begin{bmatrix} D_{0}(\bm{L}_{1}) & D_{0}(\bm{L}_{1},\bm{L}_{2})\\ D_{0}(\bm{L}_{2},\bm{L}_{1}) & D_{0}(\bm{L}_{2}) \end{bmatrix} =\begin{bmatrix} \bm{Q}_{11} & \bm{Q}_{12}\\ \bm{Q}_{12}^{T} & \bm{Q}_{22} \end{bmatrix} \tag{3-3-4}\] 上式即为观测值 \(\bm{L}\) 的协因数阵,由此得观测值的权阵为 \[\bm{P}_{0}(\bm{L})=D_{0}^{-1}(\bm{L}) =\begin{bmatrix}\bm{Q}_{11}&\bm{Q}_{12}\\ \bm{Q}_{12}^{T}&\bm{Q}_{22}\end{bmatrix}^{-1} =\begin{bmatrix}\bm{P}_{11}&\bm{P}_{12}\\ \bm{P}_{12}^{T}&\bm{P}_{22}\end{bmatrix} \tag{3-3-5}\] 为推导公式之方便,现采用如下符号: \[\begin{gathered} \widetilde{\bm{Q}}_{1}=\begin{bmatrix}\bm{Q}_{11}&\bm{0}\\ \bm{0}&\bm{0}\end{bmatrix},\quad \widetilde{\bm{Q}}_{2}=\begin{bmatrix}\bm{0}&\bm{Q}_{12}\\ \bm{Q}_{12}^{T}&\bm{0}\end{bmatrix},\quad \widetilde{\bm{Q}}_{3}=\begin{bmatrix}\bm{0}&\bm{0}\\ \bm{0}&\bm{Q}_{22}\end{bmatrix}\\ \widetilde{\bm{P}}_{1}=\begin{bmatrix}\bm{P}_{11}&\bm{0}\\ \bm{0}&\bm{0}\end{bmatrix},\quad \widetilde{\bm{P}}_{2}=\begin{bmatrix}\bm{0}&\bm{P}_{12}\\ \bm{P}_{12}^{T}&\bm{0}\end{bmatrix},\quad \widetilde{\bm{P}}_{3}=\begin{bmatrix}\bm{0}&\bm{0}\\ \bm{0}&\bm{P}_{22}\end{bmatrix} \end{gathered}\] 则 (3-3-1)、(3-3-4) 和 (3-3-5) 三式可表示为 \[\left.\begin{aligned} D(\bm{L})&=\sum_{i=1}^{3}\widetilde{\bm{Q}}_{i}\theta_{i}\\ D_{0}(\bm{L})&=\sum_{i=1}^{3}\widetilde{\bm{Q}}_{i}\\ \bm{P}_{0}(\bm{L})&=\sum_{i=1}^{3}\widetilde{\bm{P}}_{i} \end{aligned}\right\} \tag{3-3-6}\]
把 3-2 节的赫尔默特方差估计推广到相关观测:3-2 假设两类相互独立(\(\bm{P}_{12}=\bm{0}\)),本节去掉这个假设,观测方差阵中出现互协方差块 \(\sigma_{0_{12}}^{2}\bm{Q}_{12}\),待估量由 \(m\) 个 \(\sigma_{0_{i}}^{2}\) 增加到 \(m(m+1)/2\) 个(各类方差 \(\sigma_{0_{i}}^{2}\) 加两类间协方差 \(\sigma_{0_{ij}}^{2}\))。为统一起见,(3-3-6) 把 \(\bm{Q}_{11}\)、\(\bm{Q}_{12}\)、\(\bm{Q}_{22}\) 分别嵌成对角/反对角矩阵 \(\widetilde{\bm{Q}}_{i}\),使 \(D(\bm{L})=\sum_{i}\widetilde{\bm{Q}}_{i}\theta_{i}\) 成为待估量 \(\bm{\theta}\) 的线性结构——后面 (3-3-12) 至 (3-3-22) 的全部推导只对这一线性结构做运算,这就是 3-3 与 3-2 公式形式相近却更一般的原因。
由此可得未知参数 \(\bm{X}\) 在 \(\bm{P}_{0}(\bm{L})\) 下的最小二乘估值为 \[\hat{\bm{X}}=(\bm{B}^{T}\bm{P}_{0}(\bm{L})\bm{B})^{-1}\bm{B}^{T}\bm{P}_{0}(\bm{L})\bm{L} \tag{3-3-7}\] 在不致混淆的情况下,上式仍记为 \[\hat{\bm{X}}=\bm{N}^{-1}\bm{B}^{T}\bm{P}_{L}\bm{L} \tag{3-3-8}\]
观测值改正数 \(\bm{V}\) 在 \(\bm{P}_{0}(\bm{L})\) 下的计算公式为 \[\bm{V}=\bm{B}\hat{\bm{X}}-\bm{L}=\bm{B}\bm{N}^{-1}\bm{B}^{T}\bm{P}_{L}\bm{L}-\bm{L} =(\bm{B}\bm{N}^{-1}\bm{B}^{T}\bm{P}_{L}-\bm{E})\bm{L}=\bm{R}\bm{L} \tag{3-3-9}\] 式中 \[\bm{R}=\bm{B}\bm{N}^{-1}\bm{B}^{T}\bm{P}_{L}-\bm{E}\] 值得指出的是,上式中的 \(\bm{N}\)、\(\bm{P}_{L}\) 均是在初值 (3-3-3) 下的计算值,因此,\(\bm{R}\) 也应是在上述初值下的计算值,严格地讲,应记为 \(\bm{R}(0)\)。
顾及等式 \(\bm{L}-\bm{B}\bm{X}=\bm{\Delta}\),则有 \[\begin{aligned} \bm{V}&=\bm{B}\hat{\bm{X}}-\bm{L}=\bm{B}\bm{X}-\bm{B}\bm{X}+\bm{B}\hat{\bm{X}}-\bm{L}\\ &=(\bm{B}\hat{\bm{X}}-\bm{B}\bm{X})-(\bm{L}-\bm{B}\bm{X})\\ &=(\bm{B}\bm{N}^{-1}\bm{B}^{T}\bm{P}_{L}\bm{L}-\bm{B}\bm{N}^{-1}\bm{B}^{T}\bm{P}_{L}\bm{B}\bm{X})-(\bm{L}-\bm{B}\bm{X})\\ &=\bm{B}\bm{N}^{-1}\bm{B}^{T}\bm{P}_{L}(\bm{L}-\bm{B}\bm{X})-(\bm{L}-\bm{B}\bm{X})\\ &=(\bm{B}\bm{N}^{-1}\bm{B}^{T}\bm{P}_{L}-\bm{E})\,(\bm{L}-\bm{B}\bm{X})\\ &=\bm{R}\bm{\Delta} \end{aligned} \tag{3-3-10}\] 其二次型为 \[\bm{V}^{T}\bm{P}_{L}\bm{V}=\bm{\Delta}^{T}\bm{R}^{T}\bm{P}_{L}\bm{R}\bm{\Delta} \tag{3-3-11}\] 因 \[\begin{aligned} \bm{V}^{T}\bm{P}_{L}\bm{V}&=\operatorname{tr}(\bm{V}^{T}\bm{P}_{L}\bm{V})\\ &=\operatorname{tr}(\bm{\Delta}^{T}\bm{R}^{T}\bm{P}_{L}\bm{R}\bm{\Delta})\\ &=\operatorname{tr}(\bm{R}^{T}\bm{P}_{L}\bm{R}\bm{\Delta}\bm{\Delta}^{T}) \end{aligned}\] 对上式取期望值,并顾及 \(E(\bm{\Delta}\bm{\Delta}^{T})=D(\bm{\Delta})=D(\bm{L})\),则得 \[\begin{aligned} E(\bm{V}^{T}\bm{P}_{L}\bm{V})&=\operatorname{tr}\{E(\bm{R}^{T}\bm{P}_{L}\bm{R}\bm{\Delta}\bm{\Delta}^{T})\}\\ &=\operatorname{tr}\{\bm{R}^{T}\bm{P}_{L}\bm{R}E(\bm{\Delta}\bm{\Delta}^{T})\} =\operatorname{tr}(\bm{R}^{T}\bm{P}_{L}\bm{R}D(\bm{\Delta}))\\ &=\operatorname{tr}(\bm{R}^{T}\bm{P}_{L}\bm{R}D(\bm{L})) \end{aligned} \tag{3-3-12}\]
从 \(\bm{V}=\bm{R}\bm{L}\) 到 (3-3-12) 的关键步骤:(3-3-10) 先借 \(\bm{L}-\bm{B}\bm{X}=\bm{\Delta}\) 把 \(\bm{V}\) 用真误差表达为 \(\bm{V}=\bm{R}\bm{\Delta}\),于是 \(\bm{V}^{T}\bm{P}_{L}\bm{V}\) 是 \(\bm{\Delta}\) 的二次型;它是标量,故等于自身的迹,把 \(\bm{\Delta}\bm{\Delta}^{T}\) 从迹中提出后再取期望,用 \(E(\bm{\Delta}\bm{\Delta}^{T})=D(\bm{\Delta})=D(\bm{L})\) 即得 \(E(\bm{V}^{T}\bm{P}_{L}\bm{V})=\operatorname{tr}(\bm{R}^{T}\bm{P}_{L}\bm{R}D(\bm{L}))\)。与 3-2 逐类处理 \(\bm{V}_{i}^{T}\bm{P}_{i}\bm{V}_{i}\) 不同,这里先整体后拆分:(3-3-13) 把 \(\bm{P}_{L}=\sum\widetilde{\bm{P}}_{i}\) 与 \(D(\bm{L})=\sum\widetilde{\bm{Q}}_{j}\theta_{j}\) 对齐展开,即得系数阵 \(\bm{T}_{ij}=\operatorname{tr}(\bm{R}^{T}\widetilde{\bm{P}}_{i}\bm{R}\widetilde{\bm{Q}}_{j})\) 和右端 \(\bm{W}_{V_{i}}=\bm{V}^{T}\widetilde{\bm{P}}_{i}\bm{V}\)。
由 (3-3-1) 式知,所需估计的方差-协方差因子 \(\sigma_{0_{1}}^{2}\)、\(\sigma_{0_{12}}^{2}\) 和 \(\sigma_{0_{2}}^{2}\) 已包含在 \(D(\bm{L})\) 中。将上式左边的期望符号去掉,则可得方差-协方差因子的估值 \(\hat{\sigma}_{0_{1}}^{2}\)、\(\hat{\sigma}_{0_{12}}^{2}\) 和 \(\hat{\sigma}_{0_{2}}^{2}\)。
为把 (3-3-12) 式表示成方程组的形式,可将 (3-3-6) 式代入 (3-3-12) 式的左边,得 \[E(\bm{V}^{T}\bm{P}_{L}\bm{V})=\sum_{i=1}^{3}E(\bm{V}^{T}\widetilde{\bm{P}}_{i}\bm{V})\] 代入 (3-3-12) 式的右边,得 \[\operatorname{tr}(\bm{R}^{T}\bm{P}_{L}\bm{R}D(\bm{L})) =\sum_{i=1}^{3}\sum_{j=1}^{3}\operatorname{tr}(\bm{R}^{T}\widetilde{\bm{P}}_{i}\bm{R}\widetilde{\bm{Q}}_{j}\theta_{j})\] 因此有等式 \[\sum_{i=1}^{3}E(\bm{V}^{T}\widetilde{\bm{P}}_{i}\bm{V}) =\sum_{i=1}^{3}\sum_{j=1}^{3}\operatorname{tr}(\bm{R}^{T}\widetilde{\bm{P}}_{i}\bm{R}\widetilde{\bm{Q}}_{j}\theta_{j}) \tag{3-3-13}\] 很显然,上式为一具有 3 个未知数的线性方程组,去掉数学期望符号,并将它写成矩阵形式: \[\begin{bmatrix} \operatorname{tr}(\bm{R}^{T}\widetilde{\bm{P}}_{1}\bm{R}\widetilde{\bm{Q}}_{1}) & \operatorname{tr}(\bm{R}^{T}\widetilde{\bm{P}}_{1}\bm{R}\widetilde{\bm{Q}}_{2}) & \operatorname{tr}(\bm{R}^{T}\widetilde{\bm{P}}_{1}\bm{R}\widetilde{\bm{Q}}_{3})\\ \operatorname{tr}(\bm{R}^{T}\widetilde{\bm{P}}_{2}\bm{R}\widetilde{\bm{Q}}_{1}) & \operatorname{tr}(\bm{R}^{T}\widetilde{\bm{P}}_{2}\bm{R}\widetilde{\bm{Q}}_{2}) & \operatorname{tr}(\bm{R}^{T}\widetilde{\bm{P}}_{2}\bm{R}\widetilde{\bm{Q}}_{3})\\ \operatorname{tr}(\bm{R}^{T}\widetilde{\bm{P}}_{3}\bm{R}\widetilde{\bm{Q}}_{1}) & \operatorname{tr}(\bm{R}^{T}\widetilde{\bm{P}}_{3}\bm{R}\widetilde{\bm{Q}}_{2}) & \operatorname{tr}(\bm{R}^{T}\widetilde{\bm{P}}_{3}\bm{R}\widetilde{\bm{Q}}_{3}) \end{bmatrix} \begin{bmatrix} \hat{\sigma}_{0_{1}}^{2}\\ \hat{\sigma}_{0_{12}}^{2}\\ \hat{\sigma}_{0_{2}}^{2} \end{bmatrix} = \begin{bmatrix} \bm{V}^{T}\widetilde{\bm{P}}_{1}\bm{V}\\ \bm{V}^{T}\widetilde{\bm{P}}_{2}\bm{V}\\ \bm{V}^{T}\widetilde{\bm{P}}_{3}\bm{V} \end{bmatrix} \tag{3-3-14}\] 顾及 \[\bm{V}=\begin{bmatrix}\bm{V}_{1}\\ \bm{V}_{2}\end{bmatrix}\] 及符号 \(\widetilde{\bm{P}}_{i}\) 的含义,上式也可表示为 \[\begin{bmatrix} \operatorname{tr}(\bm{R}^{T}\widetilde{\bm{P}}_{1}\bm{R}\widetilde{\bm{Q}}_{1}) & \operatorname{tr}(\bm{R}^{T}\widetilde{\bm{P}}_{1}\bm{R}\widetilde{\bm{Q}}_{2}) & \operatorname{tr}(\bm{R}^{T}\widetilde{\bm{P}}_{1}\bm{R}\widetilde{\bm{Q}}_{3})\\ \operatorname{tr}(\bm{R}^{T}\widetilde{\bm{P}}_{2}\bm{R}\widetilde{\bm{Q}}_{1}) & \operatorname{tr}(\bm{R}^{T}\widetilde{\bm{P}}_{2}\bm{R}\widetilde{\bm{Q}}_{2}) & \operatorname{tr}(\bm{R}^{T}\widetilde{\bm{P}}_{2}\bm{R}\widetilde{\bm{Q}}_{3})\\ \operatorname{tr}(\bm{R}^{T}\widetilde{\bm{P}}_{3}\bm{R}\widetilde{\bm{Q}}_{1}) & \operatorname{tr}(\bm{R}^{T}\widetilde{\bm{P}}_{3}\bm{R}\widetilde{\bm{Q}}_{2}) & \operatorname{tr}(\bm{R}^{T}\widetilde{\bm{P}}_{3}\bm{R}\widetilde{\bm{Q}}_{3}) \end{bmatrix} \begin{bmatrix} \hat{\sigma}_{0_{1}}^{2}\\ \hat{\sigma}_{0_{12}}^{2}\\ \hat{\sigma}_{0_{2}}^{2} \end{bmatrix} = \begin{bmatrix} \bm{V}_{1}^{T}\bm{P}_{11}\bm{V}_{1}\\ 2\bm{V}_{1}^{T}\bm{P}_{12}\bm{V}_{2}\\ \bm{V}_{2}^{T}\bm{P}_{22}\bm{V}_{2} \end{bmatrix} \tag{3-3-15}\] 上两式是用改正数的二次型,即用残差的平方和估求 \(\bm{\theta}\) 的计算公式,若将 (3-3-9) 式,即将 \(\bm{V}=\bm{R}\bm{L}\) 代入 (3-3-14) 式,可得用观测值 \(\bm{L}\) 的二次型估求 \(\bm{\theta}\) 的计算公式: \[\begin{bmatrix} \operatorname{tr}(\bm{R}^{T}\widetilde{\bm{P}}_{1}\bm{R}\widetilde{\bm{Q}}_{1}) & \operatorname{tr}(\bm{R}^{T}\widetilde{\bm{P}}_{1}\bm{R}\widetilde{\bm{Q}}_{2}) & \operatorname{tr}(\bm{R}^{T}\widetilde{\bm{P}}_{1}\bm{R}\widetilde{\bm{Q}}_{3})\\ \operatorname{tr}(\bm{R}^{T}\widetilde{\bm{P}}_{2}\bm{R}\widetilde{\bm{Q}}_{1}) & \operatorname{tr}(\bm{R}^{T}\widetilde{\bm{P}}_{2}\bm{R}\widetilde{\bm{Q}}_{2}) & \operatorname{tr}(\bm{R}^{T}\widetilde{\bm{P}}_{2}\bm{R}\widetilde{\bm{Q}}_{3})\\ \operatorname{tr}(\bm{R}^{T}\widetilde{\bm{P}}_{3}\bm{R}\widetilde{\bm{Q}}_{1}) & \operatorname{tr}(\bm{R}^{T}\widetilde{\bm{P}}_{3}\bm{R}\widetilde{\bm{Q}}_{2}) & \operatorname{tr}(\bm{R}^{T}\widetilde{\bm{P}}_{3}\bm{R}\widetilde{\bm{Q}}_{3}) \end{bmatrix} \begin{bmatrix} \hat{\sigma}_{0_{1}}^{2}\\ \hat{\sigma}_{0_{12}}^{2}\\ \hat{\sigma}_{0_{2}}^{2} \end{bmatrix} = \begin{bmatrix} \bm{L}^{T}\bm{R}^{T}\widetilde{\bm{P}}_{1}\bm{R}\bm{L}\\ \bm{L}^{T}\bm{R}^{T}\widetilde{\bm{P}}_{2}\bm{R}\bm{L}\\ \bm{L}^{T}\bm{R}^{T}\widetilde{\bm{P}}_{3}\bm{R}\bm{L} \end{bmatrix} \tag{3-3-16}\]
下面将两类观测值的方差-协方差的估计公式扩展到 \(m\) 类观测值的一般情况,此时,观测值 \(\bm{L}\) 的方差阵为 \[\resizebox{\textwidth}{!}{$ D(\bm{L})=\begin{bmatrix} D(\bm{L}_{1}) & D(\bm{L}_{1},\bm{L}_{2}) & \cdots & D(\bm{L}_{1},\bm{L}_{m})\\ D(\bm{L}_{2},\bm{L}_{1}) & D(\bm{L}_{2}) & \cdots & D(\bm{L}_{2},\bm{L}_{m})\\ \vdots & \vdots & & \vdots\\ D(\bm{L}_{m},\bm{L}_{1}) & D(\bm{L}_{m},\bm{L}_{2}) & \cdots & D(\bm{L}_{m}) \end{bmatrix} =\begin{bmatrix} \sigma_{0_{1}}^{2}\bm{Q}_{11} & \sigma_{0_{12}}^{2}\bm{Q}_{12} & \cdots & \sigma_{0_{1m}}^{2}\bm{Q}_{1m}\\ \sigma_{0_{12}}^{2}\bm{Q}_{12}^{T} & \sigma_{0_{2}}^{2}\bm{Q}_{22} & \cdots & \sigma_{0_{2m}}^{2}\bm{Q}_{2m}\\ \vdots & \vdots & & \vdots\\ \sigma_{0_{1m}}^{2}\bm{Q}_{1m}^{T} & \sigma_{0_{2m}}^{2}\bm{Q}_{2m}^{T} & \cdots & \sigma_{0_{m}}^{2}\bm{Q}_{mm} \end{bmatrix}$}\] 式中待估量共有 \(m(m+1)/2\) 个,令 \[k=m(m+1)/2\] 其中有 \(m\) 个各类观测值的单位权方差因子 \(\sigma_{0_{i}}^{2}\),以及 \(m(m-1)/2\) 个两类观测值之间的单位权协方差因子 \(\sigma_{0_{ij}}^{2}\)。将待估量用向量 \(\bm{\theta}\) 表示,并记为 \[\begin{aligned} \underset{k\times 1}{\bm{\theta}}&=\begin{bmatrix}\sigma_{0_{1}}^{2}&\sigma_{0_{12}}^{2}&\sigma_{0_{2}}^{2}&\sigma_{0_{13}}^{2}&\sigma_{0_{23}}^{2}&\sigma_{0_{3}}^{2}&\cdots&\sigma_{0_{(m-1)m}}^{2}&\sigma_{0_{m}}^{2}\end{bmatrix}^{T}\\ &=\begin{bmatrix}\theta_{1}&\theta_{2}&\theta_{3}&\cdots&\theta_{k-1}&\theta_{k}\end{bmatrix}^{T} \end{aligned}\] (3-3-6) 式相应地扩展为 \[\left.\begin{aligned} D(\bm{L})&=\sum_{i=1}^{k}\widetilde{\bm{Q}}_{i}\theta_{i}\\ D_{0}(\bm{L})&=\sum_{i=1}^{k}\widetilde{\bm{Q}}_{i}\\ \bm{P}_{0}(\bm{L})&=\sum_{i=1}^{k}\widetilde{\bm{P}}_{i} \end{aligned}\right\} \tag{3-3-17}\] (3-3-13) 式则变为 \[\sum_{i=1}^{k}E(\bm{V}^{T}\widetilde{\bm{P}}_{i}\bm{V}) =\sum_{i=1}^{k}\sum_{j=1}^{k}\operatorname{tr}(\bm{R}^{T}\widetilde{\bm{P}}_{i}\bm{R}\widetilde{\bm{Q}}_{j}\theta_{j}) \tag{3-3-18}\] 去掉期望符号,并将上式写成矩阵形式,即得 \(m\) 类观测值的方差-协方差估计公式: \[\begin{bmatrix} \operatorname{tr}(\bm{R}^{T}\widetilde{\bm{P}}_{1}\bm{R}\widetilde{\bm{Q}}_{1}) & \operatorname{tr}(\bm{R}^{T}\widetilde{\bm{P}}_{1}\bm{R}\widetilde{\bm{Q}}_{2}) & \cdots & \operatorname{tr}(\bm{R}^{T}\widetilde{\bm{P}}_{1}\bm{R}\widetilde{\bm{Q}}_{k})\\ \operatorname{tr}(\bm{R}^{T}\widetilde{\bm{P}}_{2}\bm{R}\widetilde{\bm{Q}}_{1}) & \operatorname{tr}(\bm{R}^{T}\widetilde{\bm{P}}_{2}\bm{R}\widetilde{\bm{Q}}_{2}) & \cdots & \operatorname{tr}(\bm{R}^{T}\widetilde{\bm{P}}_{2}\bm{R}\widetilde{\bm{Q}}_{k})\\ \vdots & \vdots & & \vdots\\ \operatorname{tr}(\bm{R}^{T}\widetilde{\bm{P}}_{k}\bm{R}\widetilde{\bm{Q}}_{1}) & \operatorname{tr}(\bm{R}^{T}\widetilde{\bm{P}}_{k}\bm{R}\widetilde{\bm{Q}}_{2}) & \cdots & \operatorname{tr}(\bm{R}^{T}\widetilde{\bm{P}}_{k}\bm{R}\widetilde{\bm{Q}}_{k}) \end{bmatrix} \begin{bmatrix} \hat{\theta}_{1}\\ \hat{\theta}_{2}\\ \vdots\\ \hat{\theta}_{k} \end{bmatrix} = \begin{bmatrix} \bm{V}^{T}\widetilde{\bm{P}}_{1}\bm{V}\\ \bm{V}^{T}\widetilde{\bm{P}}_{2}\bm{V}\\ \vdots\\ \bm{V}^{T}\widetilde{\bm{P}}_{k}\bm{V} \end{bmatrix} \tag{3-3-19}\] 或写为 \[\begin{bmatrix} \operatorname{tr}(\bm{R}^{T}\widetilde{\bm{P}}_{1}\bm{R}\widetilde{\bm{Q}}_{1}) & \operatorname{tr}(\bm{R}^{T}\widetilde{\bm{P}}_{1}\bm{R}\widetilde{\bm{Q}}_{2}) & \cdots & \operatorname{tr}(\bm{R}^{T}\widetilde{\bm{P}}_{1}\bm{R}\widetilde{\bm{Q}}_{k})\\ \operatorname{tr}(\bm{R}^{T}\widetilde{\bm{P}}_{2}\bm{R}\widetilde{\bm{Q}}_{1}) & \operatorname{tr}(\bm{R}^{T}\widetilde{\bm{P}}_{2}\bm{R}\widetilde{\bm{Q}}_{2}) & \cdots & \operatorname{tr}(\bm{R}^{T}\widetilde{\bm{P}}_{2}\bm{R}\widetilde{\bm{Q}}_{k})\\ \vdots & \vdots & & \vdots\\ \operatorname{tr}(\bm{R}^{T}\widetilde{\bm{P}}_{k}\bm{R}\widetilde{\bm{Q}}_{1}) & \operatorname{tr}(\bm{R}^{T}\widetilde{\bm{P}}_{k}\bm{R}\widetilde{\bm{Q}}_{2}) & \cdots & \operatorname{tr}(\bm{R}^{T}\widetilde{\bm{P}}_{k}\bm{R}\widetilde{\bm{Q}}_{k}) \end{bmatrix} \begin{bmatrix} \hat{\theta}_{1}\\ \hat{\theta}_{2}\\ \vdots\\ \hat{\theta}_{k} \end{bmatrix} = \begin{bmatrix} \bm{L}^{T}\bm{R}^{T}\widetilde{\bm{P}}_{1}\bm{R}\bm{L}\\ \bm{L}^{T}\bm{R}^{T}\widetilde{\bm{P}}_{2}\bm{R}\bm{L}\\ \vdots\\ \bm{L}^{T}\bm{R}^{T}\widetilde{\bm{P}}_{k}\bm{R}\bm{L} \end{bmatrix} \tag{3-3-20}\] 简记为 \[\underset{k\times k}{\bm{T}}\,\underset{k\times 1}{\hat{\bm{\theta}}}=\underset{k\times 1}{\bm{W}_{V}} \tag{3-3-21}\] 和 \[\underset{k\times k}{\bm{T}}\,\underset{k\times 1}{\hat{\bm{\theta}}}=\underset{k\times 1}{\bm{W}_{L}} \tag{3-3-22}\] 式中 \[\begin{gathered} \bm{T}_{ij}=\operatorname{tr}(\bm{R}^{T}\widetilde{\bm{P}}_{i}\bm{R}\widetilde{\bm{Q}}_{j})\\ \bm{W}_{V_{i}}=\bm{V}^{T}\widetilde{\bm{P}}_{i}\bm{V}\\ \bm{W}_{L_{i}}=\bm{L}^{T}\bm{R}^{T}\widetilde{\bm{P}}_{i}\bm{R}\bm{L}\\ (i,j=1,2,\cdots,k) \end{gathered}\] 常称 (3-3-21) 或 (3-3-22) 式为赫尔默特型方差-协方差的估计公式。此公式由 Graferend 等(1980)导出。
至于 (3-3-21) 式的解,则取决于系数矩阵 \(\bm{T}\) 的性质。当 \(\bm{T}\) 为满秩矩阵,即它的秩为 \[\operatorname{rk}(\bm{T})=k\] 时,方程有唯一解 \[\hat{\bm{\theta}}=\bm{T}^{-1}\bm{W}_{V} \tag{3-3-23}\] 或 \[\hat{\bm{\theta}}=\bm{T}^{-1}\bm{W}_{L} \tag{3-3-24}\]
当 \(\bm{T}\) 为降秩矩阵,即当 \[\operatorname{rk}(\bm{T})<k\] 时,方程 (3-3-21) 式有唯一的最小范数解 \[\hat{\bm{\theta}}=\bm{T}^{+}\bm{W}_{V} \tag{3-3-25}\] 或 \[\hat{\bm{\theta}}=\bm{T}^{+}\bm{W}_{L} \tag{3-3-26}\] 式中,\(\bm{T}^{+}\) 为矩阵 \(\bm{T}\) 的最小二乘最小范数逆。
方差-协方差估计的具体计算步骤如下:
(1) 将观测值按等级或按不同观测来源分类,并选取单位权方差和单位权协方差因子的初值,即确定 \(\bm{\theta}(0)\),然后定权 \(\bm{P}_{0}(\bm{L})\);
(2) 进行第一次平差,求得 \(\bm{V}^{T}\widetilde{\bm{P}}_{i}\bm{V}\);
(3) 按 (3-3-21) 式求得单位权方差因子 \(\sigma_{0_{i}}^{2}\) 和单位权协方差因子 \(\sigma_{0_{ij}}^{2}\) 的第一次估值 \(\hat{\sigma}_{0_{i}}^{2}\) 和 \(\hat{\sigma}_{0_{ij}}^{2}\),将它们代入 (3-3-17) 式,估计观测值 \(\bm{L}\) 的方差阵,并以此重新定权;
(4) 反复进行第 2 项和第 3 项,直至 \[\hat{\theta}_{1}=\hat{\theta}_{2}=\cdots=\hat{\theta}_{k}\] 为止。
进行方差-协方差估计应注意以下几个问题:
(1) 必须有足够多的多余观测值,以使估值具有良好的统计意义。即方差-协方差估计应在大子样的条件下进行。
(2) 为避免方程 (3-3-21) 出现病态或降秩,应选择适当的 \(m\) 值,即不要将观测值分类过多。
(3) 在估计过程中,单位权方差因子 \(\hat{\sigma}_{0_{i}}^{2}\) 有可能出现负值,即出现了负方差,此时,应分析并找出其原因,然后再重新分类估计。
本节列出的注意事项都是实操教训:(1) 方差分量是二次型估计,统计上需要大子样,多余观测数 \(r\) 太小则估值不可靠;(2) 分类 \(m\) 过多会使 \(\bm{T}\) 阵病态甚至降秩,此时 (3-3-21) 无唯一解,只能改用最小范数解 \(\hat{\bm{\theta}}=\bm{T}^{+}\bm{W}_{V}\)((3-3-25) 式),且 \(\bm{T}\) 越接近奇异,解对 \(\bm{W}_{V}\) 的误差越敏感;(3) 出现负方差估值说明模型或分类有问题——例如某些类间的相关性被忽略、某类的多余观测分量太小——不能强行采用,应重新分类估计。特别地,若 \(k=m(m+1)/2\) 接近甚至超过多余观测数 \(r\),方程组本身就不具备足够的统计信息。
赫尔默特方差估计的严密公式是方差-协方差估计公式的一个特例,即在 \(\bm{Q}_{ij}=\bm{0}\) 条件下的一种特殊情况。
本节末句"赫尔默特方差估计的严密公式是方差-协方差估计公式在 \(\bm{Q}_{ij}=\bm{0}\) 时的特例",正好把 3-2 与 3-3 串成一条线;相关观测的协方差传播、互协方差阵的概念,见《最优估计基础》第1章"期望和方差"(1.2 节)与"多维随机变量"(1.4 节)。例 3-3-1 中 \(m=1\) 退化为 \(\hat{\sigma}_{0}^{2}=\bm{V}^{T}\bm{P}\bm{V}/(n-t)\) 的结论,则是检验任何方差分量估计公式正确性的通用手段。
例 3-3-1
当 \(m=1\) 时,试由 (3-3-19) 式导出间接平差时,单位权方差的估值公式 \[\hat{\sigma}_{0}^{2}=\frac{\bm{V}^{T}\bm{P}\bm{V}}{n-t}\] 式中 \(t\) 为必要观测数。
当 \(m=1\) 时,则 (3-3-19) 式变为 \[\operatorname{tr}(\bm{R}^{T}\widetilde{\bm{P}}_{1}\bm{R}\widetilde{\bm{Q}}_{1})\,\hat{\sigma}_{0_{1}}^{2} =\bm{V}^{T}\widetilde{\bm{P}}_{1}\bm{V} \tag{3-3-27}\] 此时: \[\begin{gathered} \hat{\sigma}_{0_{1}}^{2}=\hat{\sigma}_{0}^{2},\quad \widetilde{\bm{P}}_{1}=\bm{P}_{1}=\bm{P}=\bm{Q}^{-1}\\ \widetilde{\bm{Q}}=\bm{Q}_{1}=\bm{Q}=\bm{P}^{-1},\quad \bm{R}=\bm{B}\bm{N}^{-1}\bm{B}^{T}\bm{P}-\bm{E} \end{gathered}\] 则有 \[\begin{aligned} \bm{R}^{T}\widetilde{\bm{P}}_{1}\bm{R} &=(\bm{B}\bm{N}^{-1}\bm{B}^{T}\bm{P}-\bm{E})^{T}\bm{P}(\bm{B}\bm{N}^{-1}\bm{B}^{T}\bm{P}-\bm{E})\\ &=\bm{P}-\bm{P}\bm{B}\bm{N}^{-1}\bm{B}^{T}\bm{P} \end{aligned}\] \[\bm{V}^{T}\widetilde{\bm{P}}_{1}\bm{V}_{1}=\bm{V}^{T}\bm{P}\bm{V}\] \[\begin{aligned} \operatorname{tr}(\bm{R}^{T}\widetilde{\bm{P}}_{1}\bm{R}\widetilde{\bm{Q}}_{1}) &=\operatorname{tr}\{(\bm{P}-\bm{P}\bm{B}\bm{N}^{-1}\bm{B}^{T}\bm{P})\,\bm{P}^{-1}\}\\ &=\operatorname{tr}(\underset{n\times n}{\bm{E}}-\bm{P}\bm{B}\bm{N}^{-1}\bm{B}^{T}) =\operatorname{tr}(\bm{E})-\operatorname{tr}(\bm{P}\bm{B}\bm{N}^{-1}\bm{B}^{T}) \end{aligned}\] \[=\operatorname{tr}(\bm{E})-\operatorname{tr}(\bm{B}^{T}\bm{P}\bm{B}\bm{N}^{-1}) =\operatorname{tr}(\underset{n\times n}{\bm{E}})-\operatorname{tr}(\underset{t\times t}{\bm{E}})=n-t \tag{3-3-28}\] 于是 (3-3-27) 式为 \[(n-t)\,\hat{\sigma}_{0}^{2}=\bm{V}^{T}\bm{P}\bm{V}\] 即 \[\hat{\sigma}_{0}^{2}=\frac{\bm{V}^{T}\bm{P}\bm{V}}{n-t}\] 此式说明,当 \(m=1\) 时,以上方差分量估计就是单位权方差估计。
(3-3-28) 的关键两步:先用迹的循环性质 \(\operatorname{tr}(\bm{P}\bm{B}\bm{N}^{-1}\bm{B}^{T})=\operatorname{tr}(\bm{B}^{T}\bm{P}\bm{B}\bm{N}^{-1})\),再由 \(\bm{N}=\bm{B}^{T}\bm{P}\bm{B}\) 得 \(\bm{B}^{T}\bm{P}\bm{B}\bm{N}^{-1}=\bm{E}_{t\times t}\),故迹为 \(n-t\)(即多余观测数)。
例 3-3-1 是全书公式的"退化检验":令 \(m=1\) 时 \(\widetilde{\bm{P}}_{1}=\bm{P}=\bm{Q}^{-1}\)、\(\widetilde{\bm{Q}}_{1}=\bm{Q}\),系数为 \[\operatorname{tr}(\bm{R}^{T}\widetilde{\bm{P}}_{1}\bm{R}\widetilde{\bm{Q}}_{1}) =\operatorname{tr}\{(\bm{P}-\bm{P}\bm{B}\bm{N}^{-1}\bm{B}^{T}\bm{P})\,\bm{P}^{-1}\} =\operatorname{tr}(\bm{E}-\bm{P}\bm{B}\bm{N}^{-1}\bm{B}^{T})=n-t,\] 其中 \(\bm{R}^{T}\bm{P}\bm{R}=\bm{P}-\bm{P}\bm{B}\bm{N}^{-1}\bm{B}^{T}\bm{P}\) 展开时用到 \(\bm{P}\bm{B}\bm{N}^{-1}\bm{B}^{T}\bm{P}\bm{B}\bm{N}^{-1}\bm{B}^{T}\bm{P}=\bm{P}\bm{B}\bm{N}^{-1}\bm{B}^{T}\bm{P}\)(即 \(\bm{N}=\bm{B}^{T}\bm{P}\bm{B}\))。于是 (3-3-27) 化为 \((n-t)\hat{\sigma}_{0}^{2}=\bm{V}^{T}\bm{P}\bm{V}\),与经典间接平差的单位权方差公式完全一致——今后检验任何新导出的方差分量公式,都可先令 \(m=1\) 看能否退回此式。
二次无偏估计法
本节介绍两种二次无偏估计法:一是 Rao 于 1970 年提出的最小范数二次无偏(简记为 MINQUE 法),二是 Koch 于 1980 年提出的最优不变二次无偏估计法(简记为 BIQUE 法)。
二次无偏估计法的主要目的是:如果存在不同类观测值,那么该法可估计各类观测值的方差分量,即估计各类观测值的单位权方差因子 \(\sigma_{0_{i}}^{2}\);如果同类观测值的方差存在着不同因素的不同影响,则该法也可估计这些不同因素的方差分量。例如,在距离或高差测量中,人们通常把它们的观测误差分为与距离有关的一部分和与距离无关的一部分。
二次无偏估计的基本途径:先提出估计应具有的性质,然后把满足这些性质所加的条件构成一个极值问题,在 MINQUE 中,是所谓的最小范数(也可叫最小迹)问题;在 BIQUE 中,是所谓的最小方差问题,求极值问题的解,即可得到所要的估计结果。
最小范数二次无偏估计法(MINQUE 法)
设间接平差的数学模型为
函数模型: \[\underset{n\times 1}{\bm{L}}=\underset{n\times t}{\bm{B}}\,\underset{t\times 1}{\bm{X}}+\underset{n\times 1}{\bm{\Delta}},\quad \operatorname{rk}(\bm{B})=t \tag{3-4-1}\] 随机模型: \[\left.\begin{aligned} &E(\bm{L})=\bm{B}\widetilde{\bm{X}},\ E(\bm{\Delta})=\bm{0}\\ &D(\bm{\Delta})=D(\bm{L})=\sigma_{0}^{2}\bm{P}^{-1} \end{aligned}\right\} \tag{3-4-2}\] 为推导公式更具一般性,现设观测误差向量 \(\bm{\Delta}\) 具有如下形式: \[\underset{n\times 1}{\bm{\Delta}}=\underset{n\times n_{1}}{\bm{H}_{1}}\,\underset{n_{1}\times 1}{\bm{\xi}_{1}} +\underset{n\times n_{2}}{\bm{H}_{2}}\,\underset{n_{2}\times 1}{\bm{\xi}_{2}} +\cdots+\underset{n\times n_{m}}{\bm{H}_{m}}\,\underset{n_{m}\times 1}{\bm{\xi}_{m}}=\bm{H}\bm{\xi} \tag{3-4-3}\] 式中 \[\bm{H}=\begin{bmatrix}\bm{H}_{1}&\bm{H}_{2}&\cdots&\bm{H}_{m}\end{bmatrix}\] \[\bm{\xi}^{T}=\begin{bmatrix}\bm{\xi}_{1}^{T}&\bm{\xi}_{2}^{T}&\cdots&\bm{\xi}_{m}^{T}\end{bmatrix}\] \(\bm{\xi}_{i}\) 是随机误差分量,\(\bm{H}_{i}\) 为已知的 \(n\times n_{i}\) 系数矩阵,且 \(E(\bm{\xi}_{i})=\bm{0}\),\(D(\bm{\xi}_{i})=\sigma_{0_{i}}^{2}\bm{E}_{i}\,(i=1,2,\cdots,m)\),\(D(\bm{\xi}_{i},\bm{\xi}_{j})=\bm{0}\,(i\neq j)\)。由此 \[D(\bm{L})=D(\bm{\Delta})=\sum_{i=1}^{m}\bm{H}_{i}D(\bm{\xi}_{i})\bm{H}_{i}^{T} =\sum_{i=1}^{m}\sigma_{0_{i}}^{2}\bm{H}_{i}\bm{H}_{i}^{T} =\sum_{i=1}^{m}\sigma_{0_{i}}^{2}\bm{T}_{i} \tag{3-4-4}\] 式中 \[\bm{T}_{i}=\bm{H}_{i}\bm{H}_{i}^{T}\] 在上式中,\(\bm{H}_{i}\) 是已知的,所以 \(\bm{T}_{i}\) 可以计算。而单位权方差分量 \(\sigma_{0_{i}}^{2}\) 是未知的,为书写方便,设 \[\underset{m\times 1}{\bm{\theta}}=\begin{bmatrix}\sigma_{0_{1}}^{2}&\sigma_{0_{2}}^{2}&\cdots&\sigma_{0_{m}}^{2}\end{bmatrix}^{T} =\begin{bmatrix}\theta_{1}&\theta_{2}&\cdots&\theta_{m}\end{bmatrix}^{T}\] 再设上述方差分量的任意线性函数为 \[\varOmega=\alpha_{1}\sigma_{0_{1}}^{2}+\cdots+\alpha_{m}\sigma_{0_{m}}^{2} =\sum_{i=1}^{m}\alpha_{i}\sigma_{0_{i}}^{2}=\bm{\alpha}^{T}\bm{\theta} \tag{3-4-5}\] 式中 \[\underset{m\times 1}{\bm{\alpha}}=\begin{bmatrix}\alpha_{1}&\alpha_{2}&\cdots&\alpha_{m}\end{bmatrix}^{T}\] 是已知的 \(m\) 维向量(任意)。现选取观测值向量 \(\bm{L}\) 的某个二次型去估计 \(\varOmega\),即取 \[\widetilde{\varOmega}=\bm{L}^{T}\bm{M}\bm{L} \tag{3-4-6}\] 式中,二次型对称矩阵 \(\bm{M}\) 待定。在 MINQUE 法中,要求选出的待定矩阵 \(\bm{M}\) 能使估计量 \(\widetilde{\varOmega}\) 具有如下三个性质:(1) 不变性;(2) 无偏性;(3) 最小范数条件。
(1) 不变性
所谓不变性,是指二次估计 \(\bm{L}^{T}\bm{M}\bm{L}\) 与未知参数 \(\bm{X}\) 的选择无关。下面说明,仅当对称矩阵 \(\bm{M}\) 满足条件 \[\underset{n\times n}{\bm{M}}\,\underset{n\times t}{\bm{B}}=\bm{0} \tag{3-4-7}\] 时,则二次估计 \(\bm{L}^{T}\bm{M}\bm{L}\) 是不变的。
设 \[\bm{X}=\bm{X}_{0}+\delta\bm{X}\] 作为未知参数,则观测方程变为 \[\bm{L}=\bm{B}\bm{X}+\bm{\Delta}=\bm{B}\bm{X}_{0}+\bm{B}\delta\bm{X}+\bm{\Delta}\] 即有 \[\bm{l}=\bm{B}\delta\bm{X}+\bm{\Delta}\] 式中 \[\bm{l}=\bm{L}-\bm{B}\bm{X}_{0}\] 此时,\(\bm{\alpha}^{T}\bm{\theta}\) 的估计应为 \(\bm{l}^{T}\bm{M}\bm{l}\),即 \((\bm{L}-\bm{B}\bm{X}_{0})^{T}\bm{M}(\bm{L}-\bm{B}\bm{X}_{0})\),我们要求它对所有的 \(\bm{X}_{0}\) 恒等于 \(\bm{L}^{T}\bm{M}\bm{L}\),这个要求是合理的。因为现在估计的是 \(m\) 个方差分量,它们的估计量应该对于均值 \(E(\bm{L})=\bm{B}\widetilde{\bm{X}}\) 是不变的,由于 \[(\bm{L}-\bm{B}\bm{X}_{0})^{T}\bm{M}(\bm{L}-\bm{B}\bm{X}_{0}) =\bm{L}^{T}\bm{M}\bm{L}-2\bm{L}^{T}\bm{M}\bm{B}\bm{X}_{0}+\bm{X}_{0}^{T}\bm{B}^{T}\bm{M}\bm{B}\bm{X}_{0}\] 易知,当 \(\bm{M}\bm{B}=\bm{0}\) 时,上式右边的后两项均为零,则 \[(\bm{L}-\bm{B}\bm{X}_{0})^{T}\bm{M}(\bm{L}-\bm{B}\bm{X}_{0})=\bm{L}^{T}\bm{M}\bm{L} \tag{3-4-8}\] 由上式可知,在 \(\bm{M}\bm{B}=\bm{0}\) 的条件下,则对任意的 \(\bm{X}_{0}\),二次型 \(\bm{L}^{T}\bm{M}\bm{L}\) 的值不变。特别是当取 \[\bm{X}_{0}=\hat{\bm{X}}\] 式中 \(\hat{\bm{X}}\) 为 \(\bm{X}\) 的最小二乘估值,此时 \[\bm{L}-\bm{B}\bm{X}_{0}=\bm{L}-\bm{B}\hat{\bm{X}}=-(\bm{B}\hat{\bm{X}}-\bm{L})=-\bm{V}\] 则 \[(\bm{L}-\bm{B}\bm{X}_{0})^{T}\bm{M}(\bm{L}-\bm{B}\bm{X}_{0}) =(\bm{L}-\bm{B}\hat{\bm{X}})^{T}\bm{M}(\bm{L}-\bm{B}\hat{\bm{X}})=\bm{V}^{T}\bm{M}\bm{V} \tag{3-4-9}\] 顾及 (3-4-8) 式,则 \[\bm{V}^{T}\bm{M}\bm{V}=\bm{L}^{T}\bm{M}\bm{L} \tag{3-4-10}\] 由上式可知,一个观测值向量 \(\bm{L}\) 的二次型可以用它的改正数向量 \(\bm{V}\) 的二次型来表达。
(2) 无偏性
这里指的是在满足不变性条件下的无偏性。即二次估计 \(\bm{L}^{T}\bm{M}\bm{L}\) 应满足等式: \[E(\bm{L}^{T}\bm{M}\bm{L})=\bm{\alpha}^{T}\bm{\theta}=\sum_{i=1}^{m}\alpha_{i}\sigma_{0_{i}}^{2} \tag{3-4-11}\] 应用二次型的期望公式,则得 \[E(\bm{L}^{T}\bm{M}\bm{L})=\operatorname{tr}(\bm{M}D(\bm{L}))+E^{T}(\bm{L})\bm{M}E(\bm{L}) \tag{3-4-12}\] 将 \[E(\bm{L})=\bm{B}\widetilde{\bm{X}}\] \[D(\bm{L})=\sum_{i=1}^{m}\sigma_{0_{i}}^{2}\bm{T}_{i}\] 代入 (3-4-12) 式,则得 \[E(\bm{L}^{T}\bm{M}\bm{L})=\sum_{i=1}^{m}\left\{\sigma_{0_{i}}^{2}\operatorname{tr}(\bm{M}\bm{T}_{i})\right\} +\widetilde{\bm{X}}^{T}\bm{B}^{T}\bm{M}\bm{B}\widetilde{\bm{X}}\] 由不变性条件可知,上式中的 \(\bm{M}\bm{B}=\bm{0}\),故 \[E(\bm{L}^{T}\bm{M}\bm{L})=\sum_{i=1}^{m}\left\{\sigma_{0_{i}}^{2}\operatorname{tr}(\bm{M}\bm{T}_{i})\right\} \tag{3-4-13}\] 由 (3-4-11) 和 (3-4-13) 两式可知,欲使 \(\bm{L}^{T}\bm{M}\bm{L}\) 为 \(\sum\limits_{i=1}^{m}\alpha_{i}\sigma_{0_{i}}^{2}\) 的无偏估值,则 \(\bm{M}\) 还必须满足条件 \[\alpha_{i}=\operatorname{tr}(\bm{M}\bm{T}_{i})\quad(i=1,2,\cdots,m) \tag{3-4-14}\] 换言之,当待定阵 \(\bm{M}\) 满足条件 \(\alpha_{i}=\operatorname{tr}(\bm{M}\bm{T}_{i})\) 时,则 \(\bm{L}^{T}\bm{M}\bm{L}\) 为 \(\sum\limits_{i=1}^{m}\alpha_{i}\sigma_{0_{i}}^{2}\) 的无偏估值。顾及 (3-4-10) 式,则在相同条件下,\(\bm{V}^{T}\bm{M}\bm{V}\) 也是 \(\sum\limits_{i=1}^{m}\alpha_{i}\sigma_{0_{i}}^{2}\) 的无偏估值。
(3) 最小范数条件
这里指的是在不变和无偏条件下的最小范数条件,也叫最小范数准则,在这个准则下,即可具体确定待定矩阵 \(\bm{M}\)。
现假设随机变量 \(\underset{n_{i}\times 1}{\bm{\xi}_{i}}\) 是已知的,则 \(\sigma_{0_{i}}^{2}\) 的理论估值应为 \[\sigma_{0_{i}}^{2}=\frac{\bm{\xi}_{i}^{T}\bm{\xi}_{i}}{n_{i}} \tag{3-4-15}\] 则 \(\varOmega=\bm{\alpha}^{T}\bm{\theta}\) 的理论估值为 \[\begin{aligned} \varOmega=\bm{\alpha}^{T}\bm{\theta}&=\sum_{i=1}^{m}\alpha_{i}\sigma_{0_{i}}^{2} =\left(\frac{\alpha_{1}}{n_{1}}\right)\bm{\xi}_{1}^{T}\bm{\xi}_{1}+\cdots+\left(\frac{\alpha_{m}}{n_{m}}\right)\bm{\xi}_{m}^{T}\bm{\xi}_{m}\\ &=\begin{bmatrix}\bm{\xi}_{1}^{T}&\cdots&\bm{\xi}_{m}^{T}\end{bmatrix} \begin{bmatrix} \dfrac{\alpha_{1}}{n_{1}}\bm{E}_{1}&&&\\ &\ddots&&\\ &&\ddots&\\ &&&\dfrac{\alpha_{m}}{n_{m}}\bm{E}_{m} \end{bmatrix} \begin{bmatrix}\bm{\xi}_{1}\\ \vdots\\ \vdots\\ \bm{\xi}_{m}\end{bmatrix} =\bm{\xi}^{T}\bm{R}\bm{\xi} \end{aligned} \tag{3-4-16}\] 式中 \[\bm{R}=\operatorname{diag}\begin{bmatrix}\dfrac{\alpha_{1}}{n_{1}}\bm{E}_{1}&\cdots&\dfrac{\alpha_{m}}{n_{m}}\bm{E}_{m}\end{bmatrix}\] 前已提及,二次无偏估计是选取观测值向量的某个二次型 \(\bm{L}^{T}\bm{M}\bm{L}\) 去估计 \(\varOmega\)。在满足不变性的条件下,即在 \(\bm{M}\bm{B}=\bm{0}\) 的条件下,并顾及 \(\bm{\Delta}=\bm{H}\bm{\xi}\),则 \(\varOmega=\bm{\alpha}^{T}\bm{\theta}\) 的实际估值为 \[\hat{\varOmega}=\bm{L}^{T}\bm{M}\bm{L}=(\bm{B}\bm{X}+\bm{\Delta})^{T}\bm{M}(\bm{B}\bm{X}+\bm{\Delta}) =\bm{\Delta}^{T}\bm{M}\bm{\Delta}=\bm{\xi}^{T}\bm{H}^{T}\bm{M}\bm{H}\bm{\xi} \tag{3-4-17}\] 由 (3-4-16) 和 (3-4-17) 两式可知,\(\varOmega\) 的实际估值与理论估值之差为 \[\hat{\varOmega}-\varOmega=\bm{\alpha}^{T}\hat{\bm{\theta}}-\bm{\alpha}^{T}\bm{\theta} =\bm{\xi}^{T}(\bm{H}^{T}\bm{M}\bm{H}-\bm{R})\bm{\xi} \tag{3-4-18}\] 欲使 \(\hat{\varOmega}=\bm{L}^{T}\bm{M}\bm{L}\) 为一个好的估计,自然应要求实际估值与理论估值之差为最小,即要求 \(\bm{H}^{T}\bm{M}\bm{H}-\bm{R}\) 在某种意义下达到最小。最小范数二次无偏估计则选择其欧氏范数(二次范数)为最小,即适当地选择某一矩阵 \(\bm{M}\),使得 \[\|\bm{H}^{T}\bm{M}\bm{H}-\bm{R}\|^{2}=\min\] 由于欧氏范数就是矩阵中各元素的平方和,因此, \[\|\bm{H}^{T}\bm{M}\bm{H}-\bm{R}\|^{2} =\operatorname{tr}(\bm{H}^{T}\bm{M}\bm{H}\bm{H}^{T}\bm{M}\bm{H}) -2\operatorname{tr}(\bm{H}^{T}\bm{M}\bm{H}\bm{R})+\operatorname{tr}(\bm{R}^{2}) \tag{3-4-19}\] 式中 \[\operatorname{tr}(\bm{H}^{T}\bm{M}\bm{H}\bm{H}^{T}\bm{M}\bm{H}) =\operatorname{tr}(\bm{M}\bm{H}\bm{H}^{T}\bm{M}\bm{H}\bm{H}^{T})=\operatorname{tr}(\bm{M}\bm{T}\bm{M}\bm{T})\] \[\bm{T}=\bm{H}\bm{H}^{T}=\sum_{i=1}^{m}\bm{H}_{i}\bm{H}_{i}^{T}=\sum_{i=1}^{m}\bm{T}_{i}\ (\text{正定})\] \[\begin{aligned} \operatorname{tr}(\bm{H}^{T}\bm{M}\bm{H}\bm{R})&=\operatorname{tr}(\bm{M}\bm{H}\bm{R}\bm{H}^{T})\\ &=\operatorname{tr}\left(\bm{M}\frac{\alpha_{1}}{n_{1}}\bm{H}_{1}\bm{H}_{1}^{T}+\cdots+\bm{M}\frac{\alpha_{m}}{n_{m}}\bm{H}_{m}\bm{H}_{m}^{T}\right)\\ &=\frac{\alpha_{1}}{n_{1}}\operatorname{tr}(\bm{M}\bm{T}_{1})+\cdots+\frac{\alpha_{m}}{n_{m}}\operatorname{tr}(\bm{M}\bm{T}_{m})\\ &=\sum_{i=1}^{m}\frac{\alpha_{i}^{2}}{n_{i}}=\operatorname{tr}(\bm{R}^{2}) \end{aligned}\] 将这些式子代入 (3-4-19) 式,则得 \[\|\bm{H}^{T}\bm{M}\bm{H}-\bm{R}\|^{2}=\operatorname{tr}(\bm{M}\bm{T}\bm{M}\bm{T})-\operatorname{tr}(\bm{R}^{2}) \tag{3-4-20}\] 由 (3-4-16) 式可知,当向量 \(\bm{\alpha}\) 给定时,\(\operatorname{tr}(\bm{R}^{2})\) 是一常数,即与 \(\bm{M}\) 无关。所以,最小范数条件 \(\|\bm{H}^{T}\bm{M}\bm{H}-\bm{R}\|^{2}=\min\) 等价于 \[\operatorname{tr}(\bm{M}\bm{T}\bm{M}\bm{T})=\min \tag{3-4-21}\] 由以上三个性质的讨论可知,所谓最小范数二次无偏估计,其实质是当矩阵 \(\bm{M}\) 的选择,能同时满足 (3-4-7)、(3-4-14) 和 (3-4-21) 三式的条件下,利用二次型 \(\bm{L}^{T}\bm{M}\bm{L}\)(或 \(\bm{V}^{T}\bm{M}\bm{V}\))去估计方差分量的线性函数 \(\sum\limits_{i=1}^{m}\alpha_{i}\sigma_{0_{i}}^{2}\)。
MINQUE 的三条性质各有分工。不变性条件 \(\bm{M}\bm{B}=\bm{0}\) 使二次型 \(\bm{L}^{T}\bm{M}\bm{L}\) 与参数 \(\bm{X}\) 的取值无关((3-4-8) (3-4-10) 式),从而可以写成改正数的二次型 \(\bm{V}^{T}\bm{M}\bm{V}\)——若某类观测的均值未知(如自由网情形),估计方差不该受它影响;无偏性条件 \(\operatorname{tr}(\bm{M}\bm{T}_{i})=\alpha_{i}\)((3-4-14) 式)保证估计量在期望意义下对准 \(\bm{\alpha}^{T}\bm{\theta}\);最小范数条件 \(\operatorname{tr}(\bm{M}\bm{T}\bm{M}\bm{T})=\min\) 则挑出"实际估计与理论估计范数差最小"的那个 \(\bm{M}\)。三条性质依次对应估计量设计的三个维度:对均值稳健(不变)、准确(无偏)、高效(范数最小)。
采用数学命题的形式,即可叙述为:若矩阵 \(\bm{M}\) 是下述极值问题的解 \[\left\{\begin{aligned} &\text{迹最小:}\ \operatorname{tr}(\bm{M}\bm{T}\bm{M}\bm{T})=\min\\ &\text{且满足:}\ \bm{M}\bm{B}=\bm{0}\\ &\qquad\qquad\ \operatorname{tr}(\bm{M}\bm{T}_{i})=\alpha_{i}\quad(i=1,2,\cdots,m) \end{aligned}\right.\] 则称二次型 \(\bm{L}^{T}\bm{M}\bm{L}\) 为 \(\sum\limits_{i=1}^{m}\alpha_{i}\sigma_{0_{i}}^{2}\) 的最小范数二次无偏估计。为此,作拉格朗日函数: \[\varPhi(\bm{M})=2\operatorname{tr}(\bm{M}\bm{T}\bm{M}\bm{T})-4\operatorname{tr}(\bm{M}\bm{B}\bm{K}^{T}) -4\sum_{i=1}^{m}\lambda_{i}(\operatorname{tr}(\bm{M}\bm{T}_{i})-\alpha_{i})\] 式中,\(\bm{K}\) 是对应于条件 \(\bm{M}\bm{B}=\bm{0}\) 的 \(n\times t\) 联系数矩阵,\(\lambda_{i}\) 是对应于条件 (3-4-14) 式的 \(m\) 个联系数。由 \[\partial\varPhi(\bm{M})/\partial\bm{M}=\bm{0}\] 得: \[\bm{T}\bm{M}\bm{T}-\bm{K}\bm{B}^{T}-\sum_{i=1}^{m}\lambda_{i}\bm{T}_{i}=\bm{0} \tag{3-4-22}\]
由 \(\partial\varPhi/\partial\bm{M}=\bm{0}\) 得 (3-4-22) 用了矩阵迹求导公式 \(\partial\operatorname{tr}(\bm{M}\bm{T}\bm{M}\bm{T})/\partial\bm{M}=2\bm{T}\bm{M}\bm{T}\)(\(\bm{M}\)、\(\bm{T}\) 对称)与 \(\partial\operatorname{tr}(\bm{M}\bm{B}\bm{K}^{T})/\partial\bm{M}=\bm{K}\bm{B}^{T}\),代入 \(\varPhi\) 的系数即得。
把 (3-4-22) 的来路再捋一遍:目标 \(\operatorname{tr}(\bm{M}\bm{T}\bm{M}\bm{T})=\min\) 本质上是对 \(\bm{M}\) 作矩阵 Frobenius 范数最小化(\(\operatorname{tr}(\bm{A}^{T}\bm{A})=\|\bm{A}\|^{2}\)),再把 (3-4-7) 式(\(\bm{M}\bm{B}=\bm{0}\))与 (3-4-14) 式(\(\operatorname{tr}(\bm{M}\bm{T}_{i})=\alpha_{i}\))作为约束用拉格朗日乘子吸收。注意乘子的形状由约束决定:\(\bm{M}\bm{B}=\bm{0}\) 是矩阵等式,对应 \(n\times t\) 联系数矩阵 \(\bm{K}\);\(\operatorname{tr}(\bm{M}\bm{T}_{i})=\alpha_{i}\) 是标量等式,对应标量 \(\lambda_{i}\)。求导后先由 (3-4-23) 解出 \(\bm{M}\),再把约束代回确定 \(\bm{K}\) 与 \(\lambda_{i}\),这是 (3-4-23) (3-4-28) 的整体逻辑。
可见,MINQUE 法实质上可归结为求下列方程组的解: \[\begin{gathered} \bm{T}\bm{M}\bm{T}-\bm{K}\bm{B}^{T}-\sum_{i=1}^{m}\lambda_{i}\bm{T}_{i}=\bm{0}\\ \bm{M}\bm{B}=\bm{0}\\ \operatorname{tr}(\bm{M}\bm{T}_{i})=\alpha_{i} \end{gathered}\] 由 (3-4-22) 式得: \[\bm{M}=\bm{T}^{-1}\bm{K}\bm{B}^{T}\bm{T}^{-1}+\bm{T}^{-1}\left(\sum_{i=1}^{m}\lambda_{i}\bm{T}_{i}\right)\bm{T}^{-1} \tag{3-4-23}\] 再代入 (3-4-7) 式得: \[\bm{T}^{-1}\bm{K}\bm{B}^{T}\bm{T}^{-1}\bm{B}+\bm{T}^{-1}\left(\sum_{i=1}^{m}\lambda_{i}\bm{T}_{i}\right)\bm{T}^{-1}\bm{B}=\bm{0}\] 由上式即可解得联系数 \(\bm{K}\) 为 \[\bm{K}=-\left(\sum_{i=1}^{m}\lambda_{i}\bm{T}_{i}\right)\bm{T}^{-1}\bm{B}(\bm{B}^{T}\bm{T}^{-1}\bm{B})^{-1} \tag{3-4-24}\] 将上式再代入 (3-4-23) 式,即得: \[\begin{aligned} \bm{M}&=-\bm{T}^{-1}\left(\sum_{i=1}^{m}\lambda_{i}\bm{T}_{i}\right)\bm{T}^{-1}\bm{B}(\bm{B}^{T}\bm{T}^{-1}\bm{B})^{-1}\bm{B}^{T}\bm{T}^{-1} +\bm{T}^{-1}\left(\sum_{i=1}^{m}\lambda_{i}\bm{T}_{i}\right)\bm{T}^{-1}\\ &=\bm{T}^{-1}\left(\sum_{i=1}^{m}\lambda_{i}\bm{T}_{i}\right)\left(\bm{T}^{-1}-\bm{T}^{-1}\bm{B}(\bm{B}^{T}\bm{T}^{-1}\bm{B})^{-1}\bm{B}^{T}\bm{T}^{-1}\right)\\ &=\bm{T}^{-1}\left(\sum_{i=1}^{m}\lambda_{i}\bm{T}_{i}\right)\bm{C} \end{aligned} \tag{3-4-25}\] 式中 \[\bm{C}=\bm{T}^{-1}-\bm{T}^{-1}\bm{B}(\bm{B}^{T}\bm{T}^{-1}\bm{B})^{-1}\bm{B}^{T}\bm{T}^{-1}\] 由于 \(\bm{T}=\sum\limits_{i=1}^{m}\bm{T}_{i}\),\(D(\bm{L})=\sum\limits_{i=1}^{m}\sigma_{0_{i}}^{2}\bm{T}_{i}\),因此,当取 \(\sigma_{0_{i}}^{2}\) 的初值均等于 1 时,观测值向量的方差阵的初值为 \(D(\bm{L})=\bm{T}=\bm{P}^{-1}\)。此时,\(\bm{C}\) 又可表示为 \[\bm{C}=\bm{P}-\bm{P}\bm{B}(\bm{B}^{T}\bm{P}\bm{B})^{-1}\bm{B}^{T}\bm{P}=\bm{P}-\bm{P}\bm{B}\bm{N}^{-1}\bm{B}^{T}\bm{P}\] 将上式右乘 \(\bm{T}\bm{M}\),并顾及 \(\bm{P}\bm{T}=\bm{E}\) 和 \(\bm{M}\bm{B}=\bm{0}\),则得: \[\bm{C}\bm{T}\bm{M}=\bm{P}\bm{T}\bm{M}-\bm{P}\bm{B}\bm{N}^{-1}\bm{B}^{T}\bm{P}\bm{T}\bm{M} =\bm{M}-\bm{P}\bm{B}\bm{N}^{-1}\bm{B}^{T}\bm{M}=\bm{M}\] 再由 (3-4-25) 式,左乘 \(\bm{C}\bm{T}\) 后即得待定矩阵 \(\bm{M}\) 的解为 \[\bm{M}=\bm{C}\bm{T}\bm{M}=\bm{C}\left(\sum_{i=1}^{m}\lambda_{i}\bm{T}_{i}\right)\bm{C} =\sum_{i=1}^{m}\lambda_{i}\bm{C}\bm{T}_{i}\bm{C} \tag{3-4-26}\] 将上式代入 (3-4-14) 式,即得联系数 \(\lambda_{i}\) 应满足如下方程: \[\sum_{i=1}^{m}\lambda_{i}\operatorname{tr}(\bm{C}\bm{T}_{i}\bm{C}\bm{T}_{j})=\alpha_{j}\quad(i=1,2,\cdots,m) \tag{3-4-27}\] 写成矩阵形式,则为 \[\underset{m\times m}{\bm{S}}\,\underset{m\times 1}{\bm{\lambda}}=\underset{m\times 1}{\bm{\alpha}} \tag{3-4-28}\] 当 \(\bm{S}\) 满秩时,则得联系数 \(\bm{\lambda}\) 的解为 \[\bm{\lambda}=\bm{S}^{-1}\bm{\alpha}\] 式中,系数矩阵 \(\bm{S}\) 的具体表达式为 \[\begin{aligned} \underset{m\times m}{\bm{S}}&=\operatorname{tr}(\bm{C}\bm{T}_{i}\bm{C}\bm{T}_{j})\\ &=\begin{bmatrix} \operatorname{tr}(\bm{C}\bm{T}_{1}\bm{C}\bm{T}_{1}) & \operatorname{tr}(\bm{C}\bm{T}_{1}\bm{C}\bm{T}_{2}) & \cdots & \operatorname{tr}(\bm{C}\bm{T}_{1}\bm{C}\bm{T}_{m})\\ \operatorname{tr}(\bm{C}\bm{T}_{2}\bm{C}\bm{T}_{1}) & \operatorname{tr}(\bm{C}\bm{T}_{2}\bm{C}\bm{T}_{2}) & \cdots & \operatorname{tr}(\bm{C}\bm{T}_{2}\bm{C}\bm{T}_{m})\\ \vdots & \vdots & & \vdots\\ \operatorname{tr}(\bm{C}\bm{T}_{m}\bm{C}\bm{T}_{1}) & \operatorname{tr}(\bm{C}\bm{T}_{m}\bm{C}\bm{T}_{2}) & \cdots & \operatorname{tr}(\bm{C}\bm{T}_{m}\bm{C}\bm{T}_{m}) \end{bmatrix} \end{aligned} \tag{3-4-29}\] 可见 \(\bm{S}\) 为对称方阵。
将 (3-4-26) 式代入 (3-4-11) 式,并顾及 (3-4-10) 式,即得 \(\bm{\alpha}^{T}\bm{\theta}\) 的最小范数二次无偏估值为 \[\begin{aligned} \bm{\alpha}^{T}\hat{\bm{\theta}}&=\bm{L}^{T}\bm{M}\bm{L}=\bm{V}^{T}\bm{M}\bm{V}\\ &=\bm{V}^{T}\left(\sum_{i=1}^{m}\lambda_{i}\bm{C}\bm{T}_{i}\bm{C}\right)\bm{V} =\sum_{i=1}^{m}\lambda_{i}\bm{V}^{T}\bm{C}\bm{T}_{i}\bm{C}\bm{V}\\ &=\begin{bmatrix}\lambda_{1}&\lambda_{2}&\cdots&\lambda_{m}\end{bmatrix} \begin{bmatrix} \bm{V}^{T}\bm{C}\bm{T}_{1}\bm{C}\bm{V}\\ \bm{V}^{T}\bm{C}\bm{T}_{2}\bm{C}\bm{V}\\ \vdots\\ \bm{V}^{T}\bm{C}\bm{T}_{m}\bm{C}\bm{V} \end{bmatrix} =\bm{\lambda}^{T}\bm{W}_{\theta} \end{aligned} \tag{3-4-30}\] 式中 \[\bm{W}_{\theta}=\begin{bmatrix}\bm{V}^{T}\bm{C}\bm{T}_{1}\bm{C}\bm{V}&\bm{V}^{T}\bm{C}\bm{T}_{2}\bm{C}\bm{V}&\cdots&\bm{V}^{T}\bm{C}\bm{T}_{m}\bm{C}\bm{V}\end{bmatrix}^{T} \tag{3-4-31}\] 将 (3-4-29) 式代入 (3-4-30) 式,即得 \[\bm{\alpha}^{T}\hat{\bm{\theta}}=\bm{\alpha}^{T}\bm{S}^{-1}\bm{W}_{\theta}\] 由于 \(\bm{\alpha}\) 的任意性,因此,当分别取 \(\bm{\alpha}^{T}=\begin{bmatrix}1&0&\cdots&0\end{bmatrix}\),\(\bm{\alpha}^{T}=\begin{bmatrix}0&1&0&\cdots&0\end{bmatrix}\),\(\cdots\),\(\bm{\alpha}^{T}=\begin{bmatrix}0&0&\cdots&1\end{bmatrix}\) 时,即得: \[\hat{\bm{\theta}}=\bm{S}^{-1}\bm{W}_{\theta} \tag{3-4-32}\] 故单位权方差分量的估值 \(\hat{\bm{\theta}}\) 即为下述线性方程组的解: \[\bm{S}\hat{\bm{\theta}}=\bm{W}_{\theta} \tag{3-4-33}\]
最优不变二次无偏估计法(BIQUE 法)
当平差函数模型与随机模型均与 (3-4-1) (3-4-4) 式相同的情况下,BIQUE 法要求方差分量的任意线性函数 \(\bm{\alpha}^{T}\bm{\theta}\) 的二次估计 \(\bm{L}^{T}\bm{M}\bm{L}\) 应具有如下三个性质:(1) 不变性;(2) 无偏性;(3) 最小方差条件,即 \(\operatorname{var}(\bm{L}^{T}\bm{M}\bm{L})=\min\)。
可见,第一、二两个性质与 MINQUE 法相同,当然,由这两个性质推出的条件理应与 (3-4-7)、(3-4-14) 式相同。下面介绍第三个条件,即介绍 BIQUE 法中的最小方差条件,在一定的条件下与 MINQUE 法中的最小范数条件是一致的。也就是要证明:在 BIQUE 中,由最小方差条件也可导出 (3-4-21) 式: \[\operatorname{tr}(\bm{M}\bm{T}\bm{M}\bm{T})=\min\] 设观测值向量 \(\bm{L}\) 服从正态分布,即 \[\bm{L}\sim N(\bm{\mu}_{L},\bm{D}_{L})\] 式中 \[\bm{\mu}_{L}=E(\bm{L})=\bm{B}\widetilde{\bm{X}}\] \[\bm{D}_{L}=D(\bm{L})=\sigma_{0}^{2}\bm{P}^{-1}=\sum_{i=1}^{m}\sigma_{0_{i}}^{2}\bm{T}_{i}\] 根据期望为 \(\bm{\eta}\),方差阵为 \(\bm{\Sigma}\) 的正态随机向量 \(\bm{Y}\) 的二次型 \(\bm{Y}^{T}\bm{A}\bm{Y}\)(\(\bm{A}\) 为任一对称矩阵)的方差公式 \[\operatorname{var}(\bm{Y}^{T}\bm{A}\bm{Y})=2\operatorname{tr}(\bm{A}\bm{\Sigma}\bm{A}\bm{\Sigma})+4\bm{\eta}^{T}\bm{A}\bm{\Sigma}\bm{A}\bm{\eta}\] 并顾及 \(\bm{M}\bm{B}=\bm{0}\),可得二次型 \(\bm{L}^{T}\bm{M}\bm{L}\) 的方差为 \[\begin{aligned} \operatorname{var}(\bm{L}^{T}\bm{M}\bm{L})&=2\operatorname{tr}(\bm{M}\bm{D}_{L}\bm{M}\bm{D}_{L})+4\bm{\mu}_{L}^{T}\bm{M}\bm{D}_{L}\bm{M}\bm{\mu}_{L}\\ &=2\operatorname{tr}(\bm{M}\bm{D}_{L}\bm{M}\bm{D}_{L})+4\widetilde{\bm{X}}^{T}\bm{B}^{T}\bm{M}\bm{D}_{L}\bm{M}\bm{B}\widetilde{\bm{X}}\\ &=2\operatorname{tr}(\bm{M}\bm{D}_{L}\bm{M}\bm{D}_{L}) \end{aligned}\] 可见,欲使 \[\operatorname{var}(\bm{L}^{T}\bm{M}\bm{L})=\min\] 即应要求 \[\operatorname{tr}(\bm{M}\bm{D}_{L}\bm{M}\bm{D}_{L})=\min \tag{3-4-34}\] 如果存在某一对称矩阵 \(\bm{M}\) 使得上式成立,且又满足条件 \(\bm{M}\bm{B}=\bm{0}\) 和 \(\operatorname{tr}(\bm{M}\bm{T}_{i})=\alpha_{i}\),那么,二次型 \(\bm{L}^{T}\bm{M}\bm{L}\) 就是 \(\bm{\alpha}^{T}\bm{\theta}\) 的最优不变二次无偏估计量。
为使 (3-4-34) 式成立,首先要给定 \(\bm{D}_{L}\) 中的未知方差分量 \(\sigma_{0_{i}}^{2}\) 的初值。习惯做法是令 \(\sigma_{0_{i}}^{2}\) 的初值均为 1,此时观测值方差的初值即其近似值为 \[\bm{D}_{L}=D(\bm{L})=\sum_{i=1}^{m}\bm{T}_{i}=\bm{T}=\bm{P}^{-1}\] 将上式代入 (3-4-34) 式,即得 \[\operatorname{tr}(\bm{M}\bm{T}\bm{M}\bm{T})=\min\] 这表明,最小方差条件,在一定的前提下与最小范数条件是一致的。因此,由 BIQUE 导出的方差分量估计公式,理应与 MINQUE 法一致。
值得指出的是:第一,在 BIQUE 中,假设了 \(\bm{L}\) 服从正态分布,而在 MINQUE 中,则无需作上述假定。第二,在 BIQUE 中,当各方差分量的初值均等于 1 时,\(\operatorname{tr}(\bm{M}\bm{T}\bm{M}\bm{T})=\min\) 才能成立,这就意味着,其方差分量的估计公式与其近似值有关。因此,严格地说,这种估计方法应是局部最优二次无偏估计法。
BIQUE 与 MINQUE 的等价是有前提的,两点缺一不可:其一,BIQUE 要求 \(\bm{L}\) 服从正态分布,因为最小方差条件 \(\operatorname{var}(\bm{L}^{T}\bm{M}\bm{L})=\min\) 用到了正态二次型的方差公式(见 3-5 节 (3-5-4) 式:\(\operatorname{var}(\bm{Y}^{T}\bm{A}\bm{Y})=2\operatorname{tr}(\bm{A}\bm{\Sigma}\bm{A}\bm{\Sigma})+4\bm{\eta}^{T}\bm{A}\bm{\Sigma}\bm{A}\bm{\eta}\)),而 MINQUE 不需要任何分布假定;其二,最小方差条件要化成 \(\operatorname{tr}(\bm{M}\bm{T}\bm{M}\bm{T})=\min\),必须令各方差分量初值 \(\sigma_{0_{i}}^{2}=1\),使 \(D_{L}\) 的初值恰为 \(\bm{T}=\bm{P}^{-1}\),初值一变等价性即破。因此书上特意说"在一定的前提下两者一致"——严格地讲 BIQUE 只是局部最优。
最小范数二次无偏估计与赫尔默特估计
下面将从最小范数二次无偏估计的计算公式出发,导出赫尔默特估计的严密计算公式,即证明在特定的条件下,(3-2-20) 式与 (3-4-33) 式是完全一致的。
当观测值 \(\underset{n\times 1}{\bm{L}}\) 中包含有 \(m\) 类相互独立的观测分量时,则有 \[\underset{1\times n}{\bm{L}^{T}}=\begin{bmatrix}\bm{L}_{1}^{T}&\bm{L}_{2}^{T}&\cdots&\bm{L}_{m}^{T}\end{bmatrix}\] 观测误差向量为 \[\underset{1\times n}{\bm{\Delta}^{T}}=\begin{bmatrix}\bm{\Delta}_{1}^{T}&\bm{\Delta}_{2}^{T}&\cdots&\bm{\Delta}_{m}^{T}\end{bmatrix}\] 对照 (3-4-3) 式,此时 \(\bm{\Delta}\) 的具体形式可视为 \[\underset{n\times 1}{\bm{\Delta}}=\begin{bmatrix}\bm{\Delta}_{1}\\ \bm{\Delta}_{2}\\ \vdots\\ \bm{\Delta}_{m}\end{bmatrix} =\underset{n\times n}{\bm{H}}\,\underset{n\times 1}{\bm{\xi}} =\begin{bmatrix} \overline{\bm{H}}_{1}&&&\\ &\overline{\bm{H}}_{2}&&\\ &&\ddots&\\ &&&\overline{\bm{H}}_{m} \end{bmatrix} \begin{bmatrix}\bm{\xi}_{1}\\ \bm{\xi}_{2}\\ \vdots\\ \bm{\xi}_{m}\end{bmatrix}\] 即 \[\underset{n_{i}\times 1}{\bm{\Delta}_{i}}=\underset{n_{i}\times n_{i}}{\overline{\bm{H}}_{i}}\,\underset{n_{i}\times 1}{\bm{\xi}_{i}} \tag{3-4-35}\] 顾及 \(D(\bm{\xi}_{i})=\sigma_{0_{i}}^{2}\bm{E}_{i}\),则 \[D(\bm{\Delta}_{i})=\sigma_{0_{i}}^{2}\overline{\bm{H}}_{i}\overline{\bm{H}}_{i}^{T} \tag{3-4-36}\] 即有 \[\underset{n_{i}\times n_{i}}{\bm{Q}_{i}}=\underset{n_{i}\times n_{i}}{\overline{\bm{H}}_{i}}\quad\underset{n_{i}\times n_{i}}{\overline{\bm{H}}_{i}^{T}} \tag{3-4-37}\] 再将 (3-4-35) 式与 (3-4-3) 式比较可知,此时 \[\underset{n\times n_{i}}{\bm{H}_{i}}=\begin{bmatrix}\bm{0}\\ \vdots\\ \overline{\bm{H}}_{i}\\ \vdots\\ \bm{0}\end{bmatrix}\] 故 \[\underset{n\times n}{\bm{T}_{i}}=\bm{H}_{i}\bm{H}_{i}^{T} =\begin{bmatrix} \bm{0}&&&\\ &\ddots&&\\ &&\overline{\bm{H}}_{i}\overline{\bm{H}}_{i}^{T}&\\ &&&\ddots\\ &&&&\bm{0} \end{bmatrix} =\begin{bmatrix} \bm{0}&&&\\ &\ddots&&\\ &&\bm{Q}_{i}&\\ &&&\ddots\\ &&&&\bm{0} \end{bmatrix} =\underset{n\times n}{\widetilde{\bm{Q}}_{i}} \tag{3-4-38}\] \[\begin{gathered} \bm{T}=\bm{H}\bm{H}^{T}=\sum_{i=1}^{m}\bm{H}_{i}\bm{H}_{i}^{T}=\sum_{i=1}^{m}\bm{T}_{i} =\sum_{i=1}^{m}\widetilde{\bm{Q}}_{i}=\bm{Q}=\bm{P}^{-1}\\ \bm{T}^{-1}=\bm{Q}^{-1}=\bm{P} \end{gathered} \tag{3-4-39}\] 将上式代入 (3-4-25) 式中的 \(\bm{C}\),则得: \[\begin{aligned} \bm{C}&=\bm{T}^{-1}-\bm{T}^{-1}\bm{B}(\bm{B}^{T}\bm{T}^{-1}\bm{B})^{-1}\bm{B}^{T}\bm{T}^{-1} =\bm{P}-\bm{P}\bm{B}(\bm{B}^{T}\bm{P}\bm{B})^{-1}\bm{B}^{T}\bm{P}\\ &=\bm{P}-\bm{P}\bm{B}\bm{N}^{-1}\bm{B}^{T}\bm{P} \end{aligned} \tag{3-4-40}\] 将 (3-4-38) 和 (3-4-40) 两式代入 (3-4-31) 式,则得: \[\begin{aligned} \bm{W}_{i}&=\bm{V}^{T}\bm{C}\bm{T}_{i}\bm{C}\bm{V}\\ &=\bm{V}^{T}(\bm{P}-\bm{P}\bm{B}\bm{N}^{-1}\bm{B}^{T}\bm{P})\,\bm{T}_{i}\,(\bm{P}-\bm{P}\bm{B}\bm{N}^{-1}\bm{B}^{T}\bm{P})\,\bm{V}\\ &=\bm{V}^{T}\bm{P}\bm{T}_{i}\bm{P}\bm{V}-2\bm{V}^{T}\bm{P}\bm{T}_{i}\bm{P}\bm{B}\bm{N}^{-1}\bm{B}^{T}\bm{P}\bm{V} +\bm{V}^{T}\bm{P}\bm{B}\bm{N}^{-1}\bm{B}^{T}\bm{P}\bm{T}_{i}\bm{P}\bm{B}\bm{N}^{-1}\bm{B}^{T}\bm{P}\bm{V} \end{aligned}\] 顾及到 \[\begin{gathered} \bm{P}\bm{T}_{i}\bm{P}=\operatorname{diag}\begin{bmatrix}\bm{0}&\bm{0}&\cdots&\bm{P}_{i}&\bm{0}&\cdots&\bm{0}\end{bmatrix}\\ \bm{B}^{T}\bm{P}\bm{V}=\bm{0},\quad \underset{1\times n}{\bm{V}^{T}}=\begin{bmatrix}\bm{V}_{1}^{T}&\bm{V}_{2}^{T}&\cdots&\bm{V}_{m}^{T}\end{bmatrix} \end{gathered}\] 则 \[\bm{W}_{i}=\bm{V}^{T}\bm{P}\bm{T}_{i}\bm{P}\bm{V}=\bm{V}_{i}^{T}\bm{P}_{i}\bm{V}_{i} \tag{3-4-41}\]
从 MINQUE 一般公式化简到 (3-4-41),点睛之处有二。其一:\(\bm{W}_{i}=\bm{V}^{T}\bm{C}\bm{T}_{i}\bm{C}\bm{V}\) 展开后本有三项,后两项都含 \(\bm{B}^{T}\bm{P}\bm{V}\),而由法方程知 \(\bm{B}^{T}\bm{P}\bm{V}=\bm{0}\),两项同时消失,只剩 \(\bm{V}^{T}\bm{P}\bm{T}_{i}\bm{P}\bm{V}\);又 \(\bm{P}\bm{T}_{i}\bm{P}\) 只在第 \(i\) 个子块非零,故 \(\bm{W}_{i}=\bm{V}_{i}^{T}\bm{P}_{i}\bm{V}_{i}\)——这就是"独立观测下 MINQUE 的右端退化为残差平方和"的代数原因。其二:系数阵化简中 \(\bm{C}\bm{T}_{i}=\widetilde{\bm{E}}_{i}-\bm{P}\bm{B}\bm{N}^{-1}\widetilde{\bm{B}}_{i}^{T}\)(\(\widetilde{\bm{E}}_{i}\)、\(\widetilde{\bm{B}}_{i}\) 只在第 \(i\) 块非零),展开 \(S_{ii}\) 得 \(n_{i}-2\operatorname{tr}(\bm{N}^{-1}\bm{N}_{i})+\operatorname{tr}(\bm{N}^{-1}\bm{N}_{i})^{2}\),展开 \(S_{ij}\) 得 \(\operatorname{tr}(\bm{N}^{-1}\bm{N}_{i}\bm{N}^{-1}\bm{N}_{j})\),与 (3-2-20) 逐项吻合。
由 (3-4-29) 式可知,系数矩阵 \(\bm{S}\) 中的主对角线元素和非主对角线元素分别为 \[\left.\begin{aligned} \bm{S}_{ii}&=\operatorname{tr}(\bm{C}\bm{T}_{i}\bm{C}\bm{T}_{i})\quad(i=1,2,\cdots,m)\\ \bm{S}_{ij}&=\operatorname{tr}(\bm{C}\bm{T}_{i}\bm{C}\bm{T}_{j})\quad(i\neq j) \end{aligned}\right\} \tag{3-4-42}\] 顾及 (3-4-38) 和 (3-4-40) 两式,则有 \[\bm{C}\bm{T}_{i}=(\bm{P}-\bm{P}\bm{B}\bm{N}^{-1}\bm{B}^{T}\bm{P})\,\widetilde{\bm{Q}}_{i} =\bm{P}\widetilde{\bm{Q}}_{i}-\bm{P}\bm{B}\bm{N}^{-1}\bm{B}^{T}\bm{P}\widetilde{\bm{Q}}_{i} =\widetilde{\bm{E}}_{i}-\bm{P}\bm{B}\bm{N}^{-1}\widetilde{\bm{B}}_{i}^{T}\] 式中 \[\underset{n\times n}{\widetilde{\bm{E}}_{i}}=\bm{P}\widetilde{\bm{Q}}_{i} =\begin{bmatrix} \bm{0}&&&\\ &\ddots&&\\ &&\bm{E}_{i}&\\ &&&\ddots\\ &&&&\bm{0} \end{bmatrix} \longleftarrow\text{第 $i$ 个子块}\] \[\begin{aligned} \underset{t\times n}{\widetilde{\bm{B}}_{i}^{T}}=\bm{B}^{T}\bm{P}\widetilde{\bm{Q}}_{i} &=\begin{bmatrix}\bm{B}_{1}^{T}&\cdots&\bm{B}_{i}^{T}&\cdots&\bm{B}_{m}^{T}\end{bmatrix} \begin{bmatrix} \bm{0}&&&\\ &\ddots&&\\ &&\bm{E}_{i}&\\ &&&\ddots\\ &&&&\bm{0} \end{bmatrix}\\ &=\begin{bmatrix}\bm{0}&\bm{0}&\cdots&\underset{t\times n_{i}}{\bm{B}_{i}^{T}}&\bm{0}&\cdots&\bm{0}\end{bmatrix}\\ &\qquad\qquad\qquad\uparrow\\ &\qquad\quad\text{第 $i$ 个子块} \end{aligned}\] 同理 \[\bm{C}\bm{T}_{j}=(\bm{P}-\bm{P}\bm{B}\bm{N}^{-1}\bm{B}^{T}\bm{P})\,\widetilde{\bm{Q}}_{j} =\bm{P}\widetilde{\bm{Q}}_{j}-\bm{P}\bm{B}\bm{N}^{-1}\bm{B}^{T}\bm{P}\widetilde{\bm{Q}}_{j} =\widetilde{\bm{E}}_{j}-\bm{P}\bm{B}\bm{N}^{-1}\widetilde{\bm{B}}_{j}^{T}\] 式中 \[\underset{n\times n}{\widetilde{\bm{E}}_{j}}=\bm{P}\widetilde{\bm{Q}}_{j} =\begin{bmatrix} \bm{0}&&&\\ &\ddots&&\\ &&\bm{E}_{j}&\\ &&&\ddots\\ &&&&\bm{0} \end{bmatrix} \longleftarrow\text{第 $j$ 个子块}\] \[\underset{t\times n}{\widetilde{\bm{B}}_{j}^{T}}=\bm{B}^{T}\bm{P}\widetilde{\bm{Q}}_{j} =\begin{bmatrix}\bm{0}&\bm{0}&\cdots&\underset{t\times n_{j}}{\bm{B}_{j}^{T}}&\bm{0}&\cdots&\bm{0}\end{bmatrix}\] 将以上诸式代入 (3-4-42) 式,则得: \[\begin{aligned} \bm{S}_{ii}&=\operatorname{tr}\{(\bm{C}\bm{T}_{i})^{2}\} =\operatorname{tr}\{(\widetilde{\bm{E}}_{i}-\bm{P}\bm{B}\bm{N}^{-1}\widetilde{\bm{B}}_{i}^{T})^{2}\}\\ &=\operatorname{tr}(\widetilde{\bm{E}}_{i})-2\operatorname{tr}(\bm{P}\bm{B}\bm{N}^{-1}\widetilde{\bm{B}}_{i}^{T}) +\operatorname{tr}(\bm{P}\bm{B}\bm{N}^{-1}\widetilde{\bm{B}}_{i}^{T}\bm{P}\bm{B}\bm{N}^{-1}\widetilde{\bm{B}}_{i}^{T})\\ &=n_{i}-2\operatorname{tr}(\bm{N}^{-1}\bm{B}_{i}^{T}\bm{P}_{i}\bm{B}_{i}) +\operatorname{tr}(\bm{N}^{-1}\bm{B}_{i}^{T}\bm{P}_{i}\bm{B}_{i}\bm{N}^{-1}\bm{B}_{i}^{T}\bm{P}_{i}\bm{B}_{i})\\ &=n_{i}-2\operatorname{tr}(\bm{N}^{-1}\bm{N}_{i})+\operatorname{tr}(\bm{N}^{-1}\bm{N}_{i})^{2} \end{aligned}\] \[\begin{aligned} \bm{S}_{ij}&=\operatorname{tr}(\bm{C}\bm{T}_{i}\bm{C}\bm{T}_{j}) =\operatorname{tr}\{(\widetilde{\bm{E}}_{i}-\bm{P}\bm{B}\bm{N}^{-1}\widetilde{\bm{B}}_{i}^{T})\, (\widetilde{\bm{E}}_{j}-\bm{P}\bm{B}\bm{N}^{-1}\widetilde{\bm{B}}_{j}^{T})\}\\ &=\operatorname{tr}(\bm{P}\bm{B}\bm{N}^{-1}\widetilde{\bm{B}}_{i}^{T}\bm{P}\bm{B}\bm{N}^{-1}\widetilde{\bm{B}}_{j}^{T}) =\operatorname{tr}(\bm{N}^{-1}\bm{B}_{i}^{T}\bm{P}_{i}\bm{B}_{i}\bm{N}^{-1}\bm{B}_{j}^{T}\bm{P}_{j}\bm{B}_{j})\\ &=\operatorname{tr}(\bm{N}^{-1}\bm{N}_{i}\bm{N}^{-1}\bm{N}_{j}) \end{aligned}\] 将 (3-4-41) 式和上式代入方程 (3-4-33) 式,则得: \[\resizebox{0.98\textwidth}{!}{$ \begin{bmatrix} n_{1}-2\operatorname{tr}(\bm{N}^{-1}\bm{N}_{1})+\operatorname{tr}(\bm{N}^{-1}\bm{N}_{1})^{2} & \operatorname{tr}(\bm{N}^{-1}\bm{N}_{1}\bm{N}^{-1}\bm{N}_{2}) & \cdots & \operatorname{tr}(\bm{N}^{-1}\bm{N}_{1}\bm{N}^{-1}\bm{N}_{m})\\ & n_{2}-2\operatorname{tr}(\bm{N}^{-1}\bm{N}_{2})+\operatorname{tr}(\bm{N}^{-1}\bm{N}_{2})^{2} & \cdots & \operatorname{tr}(\bm{N}^{-1}\bm{N}_{2}\bm{N}^{-1}\bm{N}_{m})\\ &\text{对称}&& \vdots\\ &&& n_{m}-2\operatorname{tr}(\bm{N}^{-1}\bm{N}_{m})+\operatorname{tr}(\bm{N}^{-1}\bm{N}_{m})^{2} \end{bmatrix} \cdot\hat{\bm{\theta}} = \begin{bmatrix} \bm{V}_{1}^{T}\bm{P}_{1}\bm{V}_{1}\\ \bm{V}_{2}^{T}\bm{P}_{2}\bm{V}_{2}\\ \vdots\\ \bm{V}_{m}^{T}\bm{P}_{m}\bm{V}_{m} \end{bmatrix}$} \tag{3-4-43}\] 将上式与 (3-2-20) 式比较可知,两者完全一致。这就是说,赫尔默特估计公式可以由最小范数二次无偏估计公式导出。
(3-4-43) 与 (3-2-20) 完全一致,说明 MINQUE 在"各类观测相互独立"这一特殊情形下退化为赫尔默特估计——这是本章两条推导路线(3-2 节的验后估计路线与 3-4 节的二次无偏估计路线)的汇合点。MINQUE/BIQUE 所追求的不变性、无偏性、最小方差,与《最优估计基础》第2章中估计量的性质一脉相承(无偏性、有效性,见 2.2 节"最小二乘估计"与 2.6 节"最小方差估计");正态二次型方差公式的出处,即本书 3-5 节与本套装书第1章"多维随机变量"(1.4 节)的多维正态部分。
方差分量估计中的精度评定
前面几节,我们已经介绍了方差分量估计的多种方法,如赫尔默特(Helmert)法,最小范数二次无偏估计(MINQUE)法,以及最优不变二次无偏估计(BIQUE)法。本节将主要讨论方差分量估计中的精度评定问题。
方差分量估计公式多用改正数 \(\bm{V}\) 的二次型来表达,即方差分量因子的估值 \(\hat{\sigma}_{0_{i}}^{2}\) 是观测值改正数二次型的函数。为讨论问题方便,有必要简述一下二次型的期望、方差和协方差公式,并给予一些证明。
二次型的期望、方差和协方差公式
设 \(E(\bm{X})=\bm{\mu}_{x}\),\(E(\bm{Y})=\bm{\mu}_{y}\),\(\operatorname{var}(\bm{X})=\bm{D}_{X}\),\(\operatorname{cov}(\bm{Y},\bm{X})=\bm{D}_{YX}\),则有 \[\left.\begin{aligned} E(\bm{X}^{T}\bm{A}\bm{X})&=\operatorname{tr}(\bm{A}\bm{D}_{X})+\bm{\mu}_{x}^{T}\bm{A}\bm{\mu}_{x}\\ E(\bm{X}^{T}\bm{A}\bm{Y})&=\operatorname{tr}(\bm{A}\bm{D}_{YX})+\bm{\mu}_{x}^{T}\bm{A}\bm{\mu}_{y} \end{aligned}\right\} \tag{3-5-1}\] 式中 \(\bm{A}\) 为任意的对称可逆阵。
当 \(\bm{X}\sim N(\bm{0},\bm{D}_{X})\) 时,则有 \[\left.\begin{aligned} E(\bm{X}^{T}\bm{A}\bm{X})&=\operatorname{tr}(\bm{A}\bm{D}_{X})\\ E(\bm{X}^{T}\bm{A}\bm{Y})&=\operatorname{tr}(\bm{A}\bm{D}_{YX}) \end{aligned}\right\} \tag{3-5-2}\] 下面证明 (3-5-1) 中第一式。由 \(\bm{X}\) 的方差定义式知 \[\bm{D}_{X}=E\{(\bm{X}-\bm{\mu}_{x})(\bm{X}-\bm{\mu}_{x})^{T}\}=E(\bm{X}\bm{X}^{T})-\bm{\mu}_{x}\bm{\mu}_{x}^{T}\] 所以有 \[\begin{aligned} E(\bm{X}\bm{X}^{T})&=\bm{D}_{X}+\bm{\mu}_{x}\bm{\mu}_{x}^{T}\\ E(\bm{X}^{T}\bm{A}\bm{X})&=E\{\operatorname{tr}(\bm{X}^{T}\bm{A}\bm{X})\} =E\{\operatorname{tr}(\bm{X}\bm{X}^{T}\bm{A})\} =\operatorname{tr}\{E(\bm{X}\bm{X}^{T})\bm{A}\}\\ &=\operatorname{tr}\{(\bm{D}_{X}+\bm{\mu}_{x}\bm{\mu}_{x}^{T})\bm{A}\} =\operatorname{tr}(\bm{A}\bm{D}_{X})+\bm{\mu}_{x}^{T}\bm{A}\bm{\mu}_{x} \end{aligned} \tag{3-5-3}\]
补 (3-5-3) 证明中"标量等于自身迹"的一步:\(\bm{X}^{T}\bm{A}\bm{X}\) 是 \(1\times 1\) 标量,故 \(\bm{X}^{T}\bm{A}\bm{X}=\operatorname{tr}(\bm{X}^{T}\bm{A}\bm{X})\);用迹的循环性质 \(\operatorname{tr}(\bm{X}^{T}\bm{A}\bm{X})=\operatorname{tr}(\bm{X}\bm{X}^{T}\bm{A})\) 把随机部分移到一起,取期望后与迹交换,再代入 \(E(\bm{X}\bm{X}^{T})=\bm{D}_{X}+\bm{\mu}_{x}\bm{\mu}_{x}^{T}\) 即得 (3-5-1) 第一式。注意第二项 \(\bm{\mu}_{x}^{T}\bm{A}\bm{\mu}_{x}\) 是二次型而非迹;当 \(\bm{\mu}_{x}=\bm{0}\)(如改正数 \(\bm{V}\))时该项消失,就退回 (3-5-2) 式。证明中中心化 \(\bm{e}=\bm{X}-\bm{\mu}_{x}\)((3-5-8) 式)的手法,在 (3-5-7) 式方差推导中还要再用一次。
二次型的方差和协方差公式为:当 \(\bm{X}\sim N(\bm{\mu}_{x},\bm{D}_{X})\) 时,则有 \[\left.\begin{aligned} \operatorname{var}(\bm{X}^{T}\bm{A}\bm{X})&=2\operatorname{tr}(\bm{A}\bm{D}_{X}\bm{A}\bm{D}_{X}) +4\bm{\mu}_{x}^{T}\bm{A}\bm{D}_{X}\bm{A}\bm{\mu}_{x}\\ \operatorname{cov}(\bm{X}^{T}\bm{A}\bm{X},\bm{X}^{T}\bm{B}\bm{X})&=2\operatorname{tr}(\bm{A}\bm{D}_{X}\bm{B}\bm{D}_{X}) +4\bm{\mu}_{x}^{T}\bm{A}\bm{D}_{X}\bm{B}\bm{\mu}_{x} \end{aligned}\right\} \tag{3-5-4}\] 特别是当 \(\bm{X}\sim N(\bm{0},\bm{D}_{X})\) 时,则有 \[\left.\begin{aligned} \operatorname{var}(\bm{X}^{T}\bm{A}\bm{X})&=2\operatorname{tr}(\bm{A}\bm{D}_{X}\bm{A}\bm{D}_{X})\\ \operatorname{cov}(\bm{X}^{T}\bm{A}\bm{X},\bm{X}^{T}\bm{B}\bm{X})&=2\operatorname{tr}(\bm{A}\bm{D}_{X}\bm{B}\bm{D}_{X}) \end{aligned}\right\} \tag{3-5-5}\] 当 \(\bm{X}\sim N(\bm{0},\bm{D}_{X})\),\(\bm{Y}\sim N(\bm{0},\bm{D}_{Y})\),且 \(\operatorname{cov}(\bm{X},\bm{Y})=\bm{D}_{XY}\) 时,则有 \[\operatorname{cov}(\bm{X}^{T}\bm{A}\bm{X},\bm{Y}^{T}\bm{B}\bm{Y}) =2\operatorname{tr}(\bm{A}\bm{D}_{XY}\bm{B}\bm{D}_{YX}) \tag{3-5-6}\] 以上诸式中的 \(\bm{A}\)、\(\bm{B}\) 均为任意的对称可逆阵。
下面仅证明 (3-5-4) 式中第一式。由一维随机变量方差的定义,并顾及 (3-5-1) 式中第一式,即得 \[\begin{aligned} \operatorname{var}(\bm{X}^{T}\bm{A}\bm{X})&=E\{(\bm{X}^{T}\bm{A}\bm{X}-E(\bm{X}^{T}\bm{A}\bm{X}))^{2}\}\\ &=E\{(\bm{X}^{T}\bm{A}\bm{X}-\operatorname{tr}(\bm{A}\bm{D}_{X})-\bm{\mu}_{x}^{T}\bm{A}\bm{\mu}_{x})^{2}\} \end{aligned} \tag{3-5-7}\] 为证明简洁,现将 \(\bm{X}\) 中心化,设 \[\bm{e}=\bm{X}-\bm{\mu}_{x} \tag{3-5-8}\] \[\bm{X}=\bm{e}+\bm{\mu}_{x}\] 即可得: \[\bm{X}^{T}\bm{A}\bm{X}=\bm{e}^{T}\bm{A}\bm{e}+2\bm{\mu}_{x}^{T}\bm{A}\bm{e}+\bm{\mu}_{x}^{T}\bm{A}\bm{\mu}_{x} \tag{3-5-9}\] \[E(\bm{e})=\bm{0},\quad D(\bm{e})=\bm{D}_{X}=E(\bm{e}\bm{e}^{T}) \tag{3-5-10}\] 即 \(\bm{e}\sim N(\bm{0},\bm{D}_{X})\)。
将 (3-5-9) 式代入 (3-5-7) 式,并顾及 (3-5-10) 式,即得: \[\begin{aligned} \operatorname{var}(\bm{X}^{T}\bm{A}\bm{X}) ={}&E\{(\bm{e}^{T}\bm{A}\bm{e}+2\bm{\mu}_{x}^{T}\bm{A}\bm{e}-\operatorname{tr}(\bm{A}\bm{D}_{X}))^{2}\}\\ ={}&E\{\bm{e}^{T}\bm{A}\bm{e}\bm{e}^{T}\bm{A}\bm{e} +4\bm{\mu}_{x}^{T}\bm{A}\bm{e}\bm{e}^{T}\bm{A}\bm{\mu}_{x} +\operatorname{tr}(\bm{A}\bm{D}_{X})\operatorname{tr}(\bm{A}\bm{D}_{X})+{}\\ &4\bm{\mu}_{x}^{T}\bm{A}\bm{e}\bm{e}^{T}\bm{A}\bm{e} -2\bm{e}^{T}\bm{A}\bm{e}\operatorname{tr}(\bm{A}\bm{D}_{X}) -4\bm{\mu}_{x}^{T}\bm{A}\bm{e}\operatorname{tr}(\bm{A}\bm{D}_{X})\}\\ ={}&E(\bm{e}^{T}\bm{A}\bm{e}\bm{e}^{T}\bm{A}\bm{e}) +4\bm{\mu}_{x}^{T}\bm{A}E(\bm{e}\bm{e}^{T})\bm{A}\bm{\mu}_{x} +\operatorname{tr}(\bm{A}\bm{D}_{X})\operatorname{tr}(\bm{A}\bm{D}_{X})-{}\\ &2E(\bm{e}^{T}\bm{A}\bm{e})\operatorname{tr}(\bm{A}\bm{D}_{X})\\ ={}&\operatorname{tr}(\bm{A}\bm{D}_{X})\operatorname{tr}(\bm{A}\bm{D}_{X}) +2\operatorname{tr}(\bm{A}\bm{D}_{X}\bm{A}\bm{D}_{X}) +4\bm{\mu}_{x}^{T}\bm{A}\bm{D}_{X}\bm{A}\bm{\mu}_{x}+{}\\ &\operatorname{tr}(\bm{A}\bm{D}_{X})\operatorname{tr}(\bm{A}\bm{D}_{X}) -2\operatorname{tr}(\bm{A}\bm{D}_{X})\operatorname{tr}(\bm{A}\bm{D}_{X})\\ ={}&2\operatorname{tr}(\bm{A}\bm{D}_{X}\bm{A}\bm{D}_{X}) +4\bm{\mu}_{x}^{T}\bm{A}\bm{D}_{X}\bm{A}\bm{\mu}_{x} \end{aligned}\]
单位权方差估值 \(\hat{\sigma}_{0}^{2}\) 的精度
在最小二乘平差问题中,不管采用何种平差方法,单位权方差的估值公式均可表示为 \[\hat{\sigma}_{0}^{2}=\frac{\bm{V}^{T}\bm{P}\bm{V}}{d_{f}}=\frac{\bm{V}^{T}\bm{P}\bm{V}}{r} \tag{3-5-11}\] 式中,\(\bm{V}^{T}\bm{P}\bm{V}\) 是观测值改正数向量 \(\bm{V}\) 关于权阵 \(\bm{P}\) 的二次型,\(d_{f}\) 为平差问题的自由度,即多余观测数 \(r\)。
依广义传播律,并顾及 \(E(\bm{V})=\bm{0}\) 和 (3-5-5) 中第一式,则由 (3-5-11) 式可得估值 \(\hat{\sigma}_{0}^{2}\) 的方差为 \[\operatorname{var}(\hat{\sigma}_{0}^{2})=\operatorname{var}(\bm{V}^{T}\bm{P}\bm{V}/r) =\frac{1}{r^{2}}\operatorname{var}(\bm{V}^{T}\bm{P}\bm{V}) =\frac{2\sigma_{0}^{4}}{r^{2}}\operatorname{tr}(\bm{P}\bm{Q}_{V}\bm{P}\bm{Q}_{V}) \tag{3-5-12}\] 众所周知,在间接平差中,若设误差方程为 \[\underset{n\times 1}{\bm{V}}=\underset{n\times t}{\bm{B}}\,\underset{t\times 1}{\hat{\bm{X}}}-\underset{n\times 1}{\bm{L}} \tag{3-5-13}\] 则 \[\bm{Q}_{V}=\bm{Q}-\bm{B}(\bm{B}^{T}\bm{P}\bm{B})^{-1}\bm{B}^{T} \tag{3-5-14}\] 所以 \[\bm{P}\bm{Q}_{V}=\bm{E}-\bm{P}\bm{B}(\bm{B}^{T}\bm{P}\bm{B})^{-1}\bm{B}^{T} \tag{3-5-15}\] 易知,\(\bm{P}\bm{Q}_{V}\) 为幂等阵,即 \((\bm{P}\bm{Q}_{V})^{2}=\bm{P}\bm{Q}_{V}\),则有 \[\operatorname{tr}(\bm{P}\bm{Q}_{V}\bm{P}\bm{Q}_{V})=\operatorname{tr}(\bm{P}\bm{Q}_{V}) =\operatorname{tr}(\bm{E}-\bm{P}\bm{B}(\bm{B}^{T}\bm{P}\bm{B})^{-1}\bm{B}^{T}) =\operatorname{tr}(\underset{n\times n}{\bm{E}})-\operatorname{tr}(\underset{t\times t}{\bm{E}})=n-t=r \tag{3-5-16}\] 将上式代入 (3-5-12) 式,即得 \(\hat{\sigma}_{0}^{2}\) 的方差为 \[\operatorname{var}(\hat{\sigma}_{0}^{2})=2\sigma_{0}^{4}/r \tag{3-5-17}\] (3-5-17) 式说明,随着多余观测数的增多,单位权方差的估值 \(\hat{\sigma}_{0}^{2}\) 就愈加可靠。
(3-5-17) 的关键是 (3-5-16) 中 \(\bm{P}\bm{Q}_{V}\) 的幂等性:\((\bm{P}\bm{Q}_{V})^{2}=\bm{P}\bm{Q}_{V}\),于是 \(\operatorname{tr}(\bm{P}\bm{Q}_{V}\bm{P}\bm{Q}_{V})=\operatorname{tr}(\bm{P}\bm{Q}_{V})=\operatorname{tr}(\bm{E}-\bm{P}\bm{B}\bm{N}^{-1}\bm{B}^{T})=n-\operatorname{tr}(\bm{B}^{T}\bm{P}\bm{B}\bm{N}^{-1})=n-t=r\)。幂等性来自 \(\bm{P}\bm{Q}_{V}\) 是残差空间上的投影阵,投影阵的秩等于迹(第 2 章秩运算规则 (6))。代回 (3-5-12) 得 \(\operatorname{var}(\hat{\sigma}_{0}^{2})=\dfrac{2\sigma_{0}^{4}}{r^{2}}\cdot r=\dfrac{2\sigma_{0}^{4}}{r}\),标准差按 \(1/\sqrt{r}\) 衰减——这就是"多余观测越多,单位权方差估值越可靠"的定量表达。
Helmert 型方差分量估计中的精度评定
对两类独立观测值采用间接平差,则 Helmert 型方差分量估计模型为 \[\bm{S}\hat{\bm{\theta}}=\bm{W}_{\theta} \tag{3-5-18}\] 它的解一般可写为 \[\hat{\bm{\theta}}=\bm{S}^{-1}\bm{W}_{\theta}\] 式中 \[\left.\begin{aligned} \bm{S}&=\begin{bmatrix} n_{1}-2\operatorname{tr}(\bm{N}^{-1}\bm{N}_{1})+\operatorname{tr}(\bm{N}^{-1}\bm{N}_{1})^{2} & \operatorname{tr}(\bm{N}^{-1}\bm{N}_{1}\bm{N}^{-1}\bm{N}_{2})\\ \text{对\qquad 称} & n_{2}-2\operatorname{tr}(\bm{N}^{-1}\bm{N}_{2})+\operatorname{tr}(\bm{N}^{-1}\bm{N}_{2})^{2} \end{bmatrix}\\ &=\begin{bmatrix} \bm{S}_{11} & \bm{S}_{12}\\ \text{对称} & \bm{S}_{22} \end{bmatrix}\\ \hat{\bm{\theta}}&=\begin{bmatrix}\hat{\sigma}_{0_{1}}^{2}&\hat{\sigma}_{0_{2}}^{2}\end{bmatrix}^{T},\quad \bm{W}_{\theta}=\begin{bmatrix}\bm{V}_{1}^{T}\bm{P}_{1}\bm{V}_{1}&\bm{V}_{2}^{T}\bm{P}_{2}\bm{V}_{2}\end{bmatrix}^{T} \end{aligned}\right\} \tag{3-5-19}\] 由 (3-5-18) 式,依广义传播律可得方差分量估值 \(\hat{\bm{\theta}}\) 的方差阵为 \[\operatorname{var}(\hat{\bm{\theta}})=\bm{S}^{-1}\operatorname{var}(\bm{W}_{\theta})\bm{S}^{-1} \tag{3-5-20}\] 由 (3-5-19) 式知,为了推求 \(\operatorname{var}(\bm{W}_{\theta})\),只需分别求出 \(\operatorname{var}(\bm{V}_{1}^{T}\bm{P}_{1}\bm{V}_{1})\)、\(\operatorname{var}(\bm{V}_{2}^{T}\bm{P}_{2}\bm{V}_{2})\) 和 \(\operatorname{cov}(\bm{V}_{1}^{T}\bm{P}_{1}\bm{V}_{1},\bm{V}_{2}^{T}\bm{P}_{2}\bm{V}_{2})\) 即可。
由于 \(\bm{V}_{1}\sim N(\bm{0},\sigma_{0}^{2}\bm{Q}_{V_{1}})\),当 \(\bm{V}_{1}=\bm{B}_{1}\hat{\bm{X}}-\bm{L}_{1}\) 时,则有 \[\begin{gathered} \bm{Q}_{V_{1}}=\bm{Q}_{1}-\bm{B}_{1}\bm{N}^{-1}\bm{B}_{1}^{T}\\ \bm{P}_{1}\bm{Q}_{V_{1}}=\bm{E}-\bm{P}_{1}\bm{B}_{1}\bm{N}^{-1}\bm{B}_{1}^{T}\\ \bm{P}_{1}\bm{Q}_{V_{1}}\bm{P}_{1}\bm{Q}_{V_{1}} =\bm{E}-2\bm{P}_{1}\bm{B}_{1}\bm{N}^{-1}\bm{B}_{1}^{T} +\bm{P}_{1}\bm{B}_{1}\bm{N}^{-1}\bm{N}_{1}\bm{N}^{-1}\bm{B}_{1}^{T} \end{gathered}\] 依 (3-5-5) 式,则得: \[\begin{aligned} \operatorname{var}(\bm{V}_{1}^{T}\bm{P}_{1}\bm{V}_{1}) &=2\sigma_{0}^{4}\operatorname{tr}(\bm{P}_{1}\bm{Q}_{V_{1}}\bm{P}_{1}\bm{Q}_{V_{1}})\\ &=2\sigma_{0}^{4}\operatorname{tr}(\underset{n_{1}\times n_{1}}{\bm{E}} -2\bm{P}_{1}\bm{B}_{1}\bm{N}^{-1}\bm{B}_{1}^{T} +\bm{P}_{1}\bm{B}_{1}\bm{N}^{-1}\bm{N}_{1}\bm{N}^{-1}\bm{B}_{1}^{T})\\ &=2\sigma_{0}^{4}\left(n_{1}-2\operatorname{tr}(\bm{N}^{-1}\bm{N}_{1}) +\operatorname{tr}(\bm{N}^{-1}\bm{N}_{1})^{2}\right) =2\sigma_{0}^{4}\bm{S}_{11} \end{aligned} \tag{3-5-21}\] 同理可得: \[\operatorname{var}(\bm{V}_{2}^{T}\bm{P}_{2}\bm{V}_{2})=2\sigma_{0}^{4}\bm{S}_{22} \tag{3-5-22}\] 又因 \(\bm{V}_{1}\sim N(\bm{0},\sigma_{0}^{2}\bm{Q}_{V_{1}})\),\(\bm{V}_{2}\sim N(\bm{0},\sigma_{0}^{2}\bm{Q}_{V_{2}})\),当 \(\bm{V}_{1}=\bm{B}_{1}\hat{\bm{X}}-\bm{L}_{1}\),\(\bm{V}_{2}=\bm{B}_{2}\hat{\bm{X}}-\bm{L}_{2}\) 时,则有 \[\begin{gathered} \bm{Q}_{V_{1}V_{2}}=-\bm{B}_{1}\bm{N}^{-1}\bm{B}_{2}^{T},\quad \bm{Q}_{V_{2}V_{1}}=-\bm{B}_{2}\bm{N}^{-1}\bm{B}_{1}^{T}\\ \bm{P}_{1}\bm{Q}_{V_{1}V_{2}}\bm{P}_{2}\bm{Q}_{V_{2}V_{1}} =\bm{P}_{1}\bm{B}_{1}\bm{N}^{-1}\bm{N}_{2}\bm{N}^{-1}\bm{B}_{1}^{T} \end{gathered}\] 由 (3-5-6) 式知 \[\begin{aligned} &\operatorname{cov}(\bm{V}_{1}^{T}\bm{P}_{1}\bm{V}_{1},\bm{V}_{2}^{T}\bm{P}_{2}\bm{V}_{2}) =2\sigma_{0}^{4}\operatorname{tr}(\bm{P}_{1}\bm{Q}_{V_{1}V_{2}}\bm{P}_{2}\bm{Q}_{V_{2}V_{1}})\\ ={}&2\sigma_{0}^{4}\operatorname{tr}(\bm{P}_{1}\bm{B}_{1}\bm{N}^{-1}\bm{N}_{2}\bm{N}^{-1}\bm{B}_{1}^{T}) =2\sigma_{0}^{4}\operatorname{tr}(\bm{N}^{-1}\bm{N}_{1}\bm{N}^{-1}\bm{N}_{2})\\ ={}&2\sigma_{0}^{4}\bm{S}_{12} \end{aligned} \tag{3-5-23}\] 由此即可写出 \[\operatorname{var}(\bm{W}_{\theta}) =2\sigma_{0}^{4}\begin{bmatrix} \bm{S}_{11} & \bm{S}_{12}\\ \text{对称} & \bm{S}_{22} \end{bmatrix} =2\sigma_{0}^{4}\bm{S} \tag{3-5-24}\] 再将上式代入 (3-5-20) 式,即得: \[\operatorname{var}(\hat{\bm{\theta}})=2\sigma_{0}^{4}\bm{S}^{-1}\bm{S}\bm{S}^{-1}=2\sigma_{0}^{4}\bm{S}^{-1} \tag{3-5-25}\]
Helmert 型方差分量估值的精度公式,与最小二乘参数估计如出一辙:参数估值 \(\hat{\bm{X}}\) 的协因数阵是 \(\bm{N}^{-1}\),方差分量估值 \(\hat{\bm{\theta}}\) 的方差阵则是 \(2\sigma_{0}^{4}\bm{S}^{-1}\)——系数阵 \(\bm{S}\) 同时扮演了"法方程系数"与"协因数阵来源"的角色。原因在 (3-5-24):\(\bm{W}_{\theta}\) 的方差阵恰好是 \(2\sigma_{0}^{4}\bm{S}\),这是 \(\bm{S}\) 中 \(\operatorname{tr}\) 项的第二层含义——它刻画了残差二次型 \(\bm{V}_{1}^{T}\bm{P}_{1}\bm{V}_{1}\) 与 \(\bm{V}_{2}^{T}\bm{P}_{2}\bm{V}_{2}\) 之间的统计耦合。多余观测越多、各类观测对参数的约束越充分,\(\bm{S}\) 的"量级"越大,\(\hat{\bm{\theta}}\) 越精;反之 \(r\) 不足时 \(\bm{S}\) 近奇异,\(\bm{S}^{-1}\) 放大误差,精度公式本身就不可靠。
现设 \[\bm{S}^{-1}=\begin{bmatrix} \bm{S}_{11} & \bm{S}_{12}\\ \text{对称} & \bm{S}_{22} \end{bmatrix}^{-1} =\begin{bmatrix} \overline{\bm{S}}_{11} & \overline{\bm{S}}_{12}\\ \text{对称} & \overline{\bm{S}}_{22} \end{bmatrix}\] 并顾及 \(\hat{\bm{\theta}}=\begin{bmatrix}\hat{\sigma}_{0_{1}}^{2}&\hat{\sigma}_{0_{2}}^{2}\end{bmatrix}^{T}\),则得方差分量估值 \(\hat{\sigma}_{0_{i}}^{2}\) 的方差公式为 \[\operatorname{var}(\hat{\sigma}_{0_{i}}^{2})=2\sigma_{0}^{4}\overline{\bm{S}}_{ii} \tag{3-5-26}\] 式中,\(\overline{\bm{S}}_{ii}\) 为逆矩阵 \(\bm{S}^{-1}\) 中的第 \(i\) 个主对角线元素。
MINQUE 和 BIQUE 方差分量估计中的精度评定
在最小范数二次无偏估计与最优不变二次无偏估计中,当需估计 \(m\) 个方差分量时,其估计模型为 \[\bm{S}\hat{\bm{\theta}}=\bm{W}_{\theta} \tag{3-5-27}\] 它的解一般可写为 \[\hat{\bm{\theta}}=\bm{S}^{-1}\bm{W}_{\theta}\] 式中 \[\left.\begin{aligned} \bm{S}&=\begin{bmatrix} \operatorname{tr}(\bm{C}\bm{T}_{1}\bm{C}\bm{T}_{1}) & \operatorname{tr}(\bm{C}\bm{T}_{1}\bm{C}\bm{T}_{2}) & \cdots & \operatorname{tr}(\bm{C}\bm{T}_{1}\bm{C}\bm{T}_{m})\\ & \operatorname{tr}(\bm{C}\bm{T}_{2}\bm{C}\bm{T}_{2}) & \cdots & \operatorname{tr}(\bm{C}\bm{T}_{2}\bm{C}\bm{T}_{m})\\ \text{对\qquad 称} && \vdots\\ &&& \operatorname{tr}(\bm{C}\bm{T}_{m}\bm{C}\bm{T}_{m}) \end{bmatrix} =(\bm{S}_{ij})_{m\times m}\\ \hat{\bm{\theta}}&=\begin{bmatrix}\hat{\sigma}_{0_{1}}^{2}&\hat{\sigma}_{0_{2}}^{2}&\cdots&\hat{\sigma}_{0_{m}}^{2}\end{bmatrix}^{T}\\ \bm{W}_{\theta}&=\begin{bmatrix}\bm{V}^{T}\bm{C}\bm{T}_{1}\bm{C}\bm{V}&\bm{V}^{T}\bm{C}\bm{T}_{2}\bm{C}\bm{V}&\cdots&\bm{V}^{T}\bm{C}\bm{T}_{m}\bm{C}\bm{V}\end{bmatrix}^{T} \end{aligned}\right\} \tag{3-5-28}\] 由 (3-5-27) 式得: \[\operatorname{var}(\hat{\bm{\theta}})=\bm{S}^{-1}\operatorname{var}(\bm{W}_{\theta})\bm{S}^{-1} \tag{3-5-29}\] 下面推导 \(\operatorname{var}(\bm{W}_{\theta})\) 的表达式。因 \(\bm{V}\sim N(\bm{0},\sigma_{0}^{2}\bm{Q}_{V})\),所以,由 (3-5-5) 式(此时 \(\bm{A}=\bm{C}\bm{T}_{i}\bm{C}\))得: \[\operatorname{var}(\bm{V}^{T}\bm{C}\bm{T}_{i}\bm{C}\bm{V}) =2\sigma_{0}^{4}\operatorname{tr}(\bm{C}\bm{T}_{i}\bm{C}\bm{Q}_{V}\bm{C}\bm{T}_{i}\bm{C}\bm{Q}_{V}) =2\sigma_{0}^{4}\operatorname{tr}(\bm{C}\bm{Q}_{V}\bm{C}\bm{T}_{i}\bm{C}\bm{Q}_{V}\bm{C}\bm{T}_{i}) \tag{3-5-30}\] 设平差函数模型为 \(\bm{L}=\bm{B}\bm{X}+\bm{\Delta}\),误差方程为 \(\bm{V}=\bm{B}\hat{\bm{X}}-\bm{L}\) 时,则有 \(\bm{C}=\bm{P}-\bm{P}\bm{B}\bm{N}^{-1}\bm{B}^{T}\bm{P}\),\(\bm{Q}_{V}=\bm{Q}-\bm{B}\bm{N}^{-1}\bm{B}^{T}\),此时 \[\begin{aligned} \bm{C}\bm{Q}_{V}\bm{C}&=(\bm{P}-\bm{P}\bm{B}\bm{N}^{-1}\bm{B}^{T}\bm{P})\,(\bm{Q}-\bm{B}\bm{N}^{-1}\bm{B}^{T})\,\bm{C} =(\bm{E}-\bm{P}\bm{B}\bm{N}^{-1}\bm{B}^{T})\,\bm{C}\\ &=\bm{C}-\bm{P}\bm{B}\bm{N}^{-1}\bm{B}^{T}\bm{C} =\bm{C}-\bm{P}\bm{B}\bm{N}^{-1}\bm{B}^{T}\bm{P}+\bm{P}\bm{B}\bm{N}^{-1}\bm{N}\bm{N}^{-1}\bm{B}^{T}\bm{P} =\bm{C} \end{aligned} \tag{3-5-31}\] 将上式代入 (3-5-30) 式,即得: \[\operatorname{var}(\bm{V}^{T}\bm{C}\bm{T}_{i}\bm{C}\bm{V}) =2\sigma_{0}^{4}\operatorname{tr}(\bm{C}\bm{T}_{i}\bm{C}\bm{T}_{i}) =2\sigma_{0}^{4}\bm{S}_{ii}\quad(i=1,2,\cdots,m) \tag{3-5-32}\] 再由 (3-5-5) 式(此时 \(\bm{A}=\bm{C}\bm{T}_{i}\bm{C}\),\(\bm{B}=\bm{C}\bm{T}_{j}\bm{C}\)),并顾及 (3-5-31) 式,则得: \[\begin{aligned} &\operatorname{cov}(\bm{V}^{T}\bm{C}\bm{T}_{i}\bm{C}\bm{V},\bm{V}^{T}\bm{C}\bm{T}_{j}\bm{C}\bm{V})\\ ={}&2\sigma_{0}^{4}\operatorname{tr}(\bm{C}\bm{T}_{i}\bm{C}\bm{Q}_{V}\bm{C}\bm{T}_{j}\bm{C}\bm{Q}_{V}) =2\sigma_{0}^{4}\operatorname{tr}(\bm{C}\bm{Q}_{V}\bm{C}\bm{T}_{i}\bm{C}\bm{Q}_{V}\bm{C}\bm{T}_{j})\\ ={}&2\sigma_{0}^{4}\operatorname{tr}(\bm{C}\bm{T}_{i}\bm{C}\bm{T}_{j}) =2\sigma_{0}^{4}\bm{S}_{ij}\quad(i\neq j=1,2,\cdots,m) \end{aligned} \tag{3-5-33}\] 有了以上两式,则无需详细推导,就可得以下诸式: \[\operatorname{var}(\bm{W}_{\theta}) =2\sigma_{0}^{4}\begin{bmatrix} \bm{S}_{11} & \bm{S}_{12} & \cdots & \bm{S}_{1m}\\ & \bm{S}_{22} & \cdots & \bm{S}_{2m}\\ \text{对\quad 称} && \vdots\\ &&& \bm{S}_{mm} \end{bmatrix} =2\sigma_{0}^{4}\bm{S}\] \[\operatorname{var}(\hat{\bm{\theta}})=2\sigma_{0}^{4}\bm{S}^{-1} \tag{3-5-34}\] 式中,\(\overline{\bm{S}}_{ii}\) 为逆矩阵 \(\bm{S}^{-1}\) 中的第 \(i\) 个主对角线元素。
(3-5-34) 中的 \(\sigma_{0}^{4}\) 是未知真值,实用中只能代入估计值 \(\hat{\sigma}_{0}^{4}\),故精度公式本身带近似性;且 \(\hat{\bm{\theta}}=\bm{S}^{-1}\bm{W}_{\theta}\) 是 \(\bm{W}_{\theta}\) 的线性函数,当 \(\bm{S}\) 近奇异时 \(\bm{S}^{-1}\) 会放大 \(\bm{W}_{\theta}\) 的随机误差。量纲上,若 \(\sigma_{0}\) 是单位权中误差(量纲如 \(('')\)),则 \(\hat{\sigma}_{0_{i}}^{2}\) 的方差量纲含 \(\sigma_{0}^{4}\)。例 3-5-1 表 3-11 中 \(2.90\)、\(8.14\) 两值与 \(2\hat{\sigma}_{0}^{4}\overline{\bm{S}}_{ii}\) 基本吻合(取第三次迭代的 \(\hat{\sigma}_{0}=\pm1.90''\):\(2\times1.90^{4}\times0.1110\approx2.90\),\(2\times1.90^{4}\times0.3113\approx8.14\)),单位沿用原书记法。MINQUE/BIQUE 的 \(\bm{S}\) 阵由 \(\operatorname{tr}(\bm{C}\bm{T}_{i}\bm{C}\bm{T}_{j})\) 组成((3-5-28)),与赫尔默特型 \(\bm{S}\) 形式不同,但 (3-4-43) 已证独立观测下两者一致,故精度公式同样适用。
算例
[例 3-5-1]图 3-2 所示的边角网,其中角度为平差元素;
[例 3-5-2]图 3-3 所示的边角网,其中方向为平差元素;
[例 3-5-3]图 3-4 所示的水准网,\(\mathrm{I}\)、\(\mathrm{II}\) 等联合平差。
对上述三个算例分别进行了平差计算和方差分量估计,以确定各个平差网中两类观测值的合理权比,最后还进行了方差分量估计的精度评定计算。
表 3-9 列出了三个网的基本情况;
表 3-10 列出了三个网的方差分量估计的主要结果;
表 3-11 列出了三个网的方差分量估计的精度计算结果。
有了验后中误差,可得验后权。对 [例 3-5-1],则有:\(P_{\beta}=1\),\(P_{S}=0.61\);对 [例 3-5-2],则有:\(P_{a}=1\),\(P_{S}=5.3/(0.5+4\times 10^{-4}S_{i})^{2}\);对 [例 3-5-3],则有 \(P_{\mathrm{I}}=100/S_{i}\),\(P_{\mathrm{II}}=100/1.6S_{i}\)。
| 算例 | 单位权中误差估值 (\(\hat{\sigma}_{0}\)) | \(\bm{S}^{-1}\) | \(\hat{\operatorname{var}}(\hat{\sigma}_{0_{i}}^{2})\) |
|---|---|---|---|
| 例 3-5-1 | \(\pm1.90''\) | \(\begin{bmatrix}0.1110 & -0.0260\\ -0.0260 & 0.3113\end{bmatrix}\) | \(\begin{aligned}2.90\,(''^{2})\\ 8.14\,(''^{2})\end{aligned}\) |
| 例 3-5-2 | \(\pm1.15''\) | \(\begin{bmatrix}0.0434 & -0.0133\\ -0.0133 & 0.1472\end{bmatrix}\) | \(\begin{aligned}0.15\,(''^{2})\\ 0.51\,(''^{2})\end{aligned}\) |
| 例 3-5-3 | \(\pm7.51\,\mathrm{mm/100km}\) | \(\begin{bmatrix}0.1679 & -0.0294\\ -0.0294 & 0.1303\end{bmatrix}\) | \(\begin{aligned}1068.17\,(\mathrm{mm}^{4})\\ 828.96\,(\mathrm{mm}^{4})\end{aligned}\) |
对算例结果稍作分析,即可看出:
(1) 当网形较大,多余观测数较多时,各类观测值的方差分量估值才具有较高的精度。如在 [例 3-5-2] 中 \(r=35\),方差分量估值的中误差仅为方差分量自身数值的 \(1/3\) \(1/2\)。这表明:在自由度 \(d_{f}\) 较大时,方差分量的估值才是可靠的。当然,由此而定出的两类观测值的权比也是可靠的。
(2) 在同一平差问题中,当各类观测值的多余观测分量相差较大时(例 3-5-1,例 3-5-2),其方差分量估值的精度也相差较大。这表明,欲使各类观测值的方差分量估值具有同等的可靠度,则在进行方差分量估计时,应尽可能地使各类观测值的多余观测分量大体保持相等(如例 3-5-3)。
从估计的角度看,测量平差与方差分量估计均属参数估计问题。所不同的是:测量平差是估求未知参数 \(\bm{X}\) 的最佳估值 \(\hat{\bm{X}}\) 和观测值的平差值 \(\hat{\bm{L}}\),即估求平差函数模型中的参数;而方差分量估计是估求平差随机模型中的未知参数,即估求各类观测值的方差分量 \(\sigma_{0_{i}}^{2}\)。因此,严格地说,欲使观测数据处理得更加完善,必须是平差计算与方差分量估计同时进行。
众所周知,未知参数的估值 \(\hat{\bm{X}}\) 总可表示为观测值 \(\bm{L}\) 的一次函数,而方差分量的估值 \(\hat{\bm{\theta}}\) 则是观测值 \(\bm{L}\) 的二次函数。由于 \(\bm{L}\) 是随机向量,因此,未知参数的估值 \(\hat{\bm{X}}\) 和方差分量的估值 \(\hat{\bm{\theta}}\) 也是随机向量,正如估求 \(\hat{\bm{X}}\) 的方差一样,也必须估求 \(\hat{\bm{\theta}}\) 的方差。
在单位权方差一定的条件下,方差分量估值 \(\hat{\bm{\theta}}\) 的方差均与平差问题的自由度有关。当 \(d_{f}\) 较大时,方差分量的估值才是可靠的;换言之,由方差分量的估值重新确定各类观测值的权比也是可靠的。这就告诉我们:只有在数据处理的精度要求较高,平差网形较大,且多余观测数较多的情况下,进行方差分量估计才是必要的和有利的。