参数的先验随机性通常是指参数具有先验期望和先验方差，这些先验信息往往来源于上一个观测和平差任务。经典最小二乘平差只能针对非随机参数，即使参数具有先验期望和方差也无法在平差中顾及其随机性。在广义平差中，如极大验后估计、卡尔曼滤波等，都必须顾及参数的随机性。在参数估计时顾及先验随机信息，可以实现动态数据处理和数据分期处理，并提高数据处理的精度。整体最小二乘(total least squares, TLS)已提出的诸多算法^[1-4]均未考虑参数的先验随机性。2009年Schaffrin^[5]首先提出EIV-REM (errors-in-variables with random effects model)模型，在EIV模型基础上考虑参数随机性，称之为WTLSC加权整体最小二乘配置(weighted total least squares collocation)，同时附加了IID(independent and identically distributed data)条件，即要求数据独立同分布。2012~2013年Snow等^[6-8]研究了EIV-REM模型的相关算法，去掉了IID条件。2013年Fang^[9]研究了随机参数的WTLS参数估计算法。2017年王乐洋等^[10]在经典最小二乘配置的基础上，同时顾及观测向量和趋势项系数矩阵存在随机误差的情形，推导了总体最小二乘配置方法并应用于地壳形变分析。文献[5-9]将全部待估参数作为随机参数，而没有考虑只有部分参数具有先验随机信息的情形。而实践中只有部分参数具有先验随机信息的情况较多，例如三维激光扫描的球形标靶在进行圆曲线拟合时，半径和圆心坐标都是待估参数，其中半径可以有先验随机信息，但圆心坐标由于标靶摆放的任意性而没有先验信息；在最小二乘配置问题中，既有随机参数，又有非随机参数。本文提出部分参数具有先验随机性的加权整体最小二乘平差模型及算法，该模型算法更具一般性，可将整体最小二乘原理从经典平差领域引入到广义平差领域的应用中。同时，给出了函数模型，基于最小二乘准则利用拉格朗日条件极值，通过矩阵合并与分解，推导出待估参数与观测改正数的表达式，给出精度评定公式，通知实例对算法进行验证。

1 部分参数具有随机性的WTLS算法 1.1 函数模型

EIV模型^[11-12]为:

$ \left( {\mathop {\mathit{\boldsymbol{A}}}\limits_{m \times u} + {\mathit{\boldsymbol{E}}_{\mathop A\limits_{m \times u} }}} \right)\mathop {\mathit{\boldsymbol{X}}}\limits_{u \times 1} = \mathop {\mathit{\boldsymbol{L}}}\limits_{m \times 1} + \mathop{ \mathit{\boldsymbol{e}}}\limits_{m \times 1} $

(1a)

式中，系数矩阵A和观测向量L分别含有来源于观测值的误差E_A及e，系数矩阵A的改正矩阵为V_A，观测向量L的改正数为V_L。随机参数X具有先验期望μ_X和先验方差D_X，μ_X含有误差e_X，e_X与e不相关。EIV-REM模型为：

$ \mathop {\mathit{\boldsymbol{X}}}\limits_{u \times 1} = {\mathit{\boldsymbol{\mu }}_{\mathop X\limits_{u \times 1} }} + {\mathit{\boldsymbol{e}}_{\mathop X\limits_{u \times 1} }} $

(1b)

式中，所有u个待估参数X_i都具有先验随机信息。通常称先验期望μ_X为虚拟观测值，e_X为其误差。部分参数具有先验随机性的函数模型为：

$ \mathop {\mathit{\boldsymbol{a}}}\limits_{c \times u} \mathop {\mathit{\boldsymbol{X}}}\limits_{u \times 1} = {\mathit{\boldsymbol{\mu }}_{\mathop X\limits_{c \times 1} }} + {\mathit{\boldsymbol{e}}_{\mathop X\limits_{c \times 1} }} $

(1c)

式中，a为系数矩阵，其元素是单位矩阵的部分或全部行向量。

将式(1a)中误差矩阵E_A及e拉直并改写为观测值改正数向量V，则有：

$ \begin{array}{*{20}{c}} {\mathop {\mathit{\boldsymbol{L}}}\limits_{m \times 1} - \mathop {\mathit{\boldsymbol{A}}}\limits_{m \times u} \mathop {\mathit{\boldsymbol{X}}}\limits_{u \times 1} = \mathop {{\mathit{\boldsymbol{E}}_A}}\limits_{m \times u} \mathop {\mathit{\boldsymbol{X}}}\limits_{u \times 1} - \mathop {\mathit{\boldsymbol{e}}}\limits_{m \times 1} = }\\ {\left[ {\begin{array}{*{20}{c}} {\mathop {{\mathit{\boldsymbol{E}}_A}}\limits_{m \times \left( {u + 1} \right)} }&\mathit{\boldsymbol{e}} \end{array}} \right]\mathop {\left[ {\begin{array}{*{20}{c}} \mathit{\boldsymbol{X}}\\ { - 1} \end{array}} \right]}\limits_{\left( {u + 1} \right) \times 1} = \mathop {\mathit{\boldsymbol{B}}}\limits_{m \times n} \mathop {\mathit{\boldsymbol{V}}}\limits_{n \times 1} } \end{array} $

(2)

式中，$ \mathit{\boldsymbol{B}}\text{=}{{\mathit{\boldsymbol{I}}}_{\mathit{m}}}\otimes \left[ \frac{\partial \mathit{\boldsymbol{f}}}{\partial {{\mathit{\boldsymbol{L}}}_{\mathit{A}}}}-1 \right]$，其中I_m为m阶单位矩阵，f为观测方程，L_A为A矩阵元素中的观测值，求偏导的实质是泰勒级数展开至一次项。部分参数具有随机性的整体最小二乘平差函数模型为：

$ \left\{ \begin{array}{l} \mathop {\mathit{\boldsymbol{L}}}\limits_{m \times 1} = \mathop {\mathit{\boldsymbol{A}}}\limits_{m \times u} \mathop {\mathit{\boldsymbol{X}}}\limits_{u \times 1} + \mathop {\mathit{\boldsymbol{B}}}\limits_{m \times n} \mathop {\mathit{\boldsymbol{V}}}\limits_{n \times 1} \\ \mathop {{\mathit{\boldsymbol{\mu }}_X}}\limits_{c \times 1} = \mathop {\mathit{\boldsymbol{a}}}\limits_{c \times u} \mathop {\mathit{\boldsymbol{X}}}\limits_{u \times 1} - \mathop {{\mathit{\boldsymbol{V}}_X}}\limits_{c \times 1} \end{array} \right. $

(3)

式中，X是平差问题的u维待估参数向量，L为m阶实际观测列向量，其改正数为m阶列向量V_L，A为m×u阶系数矩阵，其改正数为m×u阶矩阵V_A，V为n维(n≥m)观测改正数向量，B为改正数向量的m×n阶系数矩阵，a为c×u阶系数矩阵，0≤c≤u，a的元素为u阶单位矩阵的一行或多行，c为0表示参数没有先验信息，c为u则所有参数具有先验信息，X的先验数学期望μ_X是c×1阶虚拟观测列向量，其改正数为V_X。合并式(3)得：

$ \left[ {\begin{array}{*{20}{c}} {\mathop {\mathit{\boldsymbol{L}}}\limits_{m \times 1} }\\ {\mathop {{\mathit{\boldsymbol{\mu }}_X}}\limits_{c \times 1} } \end{array}} \right] = \left[ {\begin{array}{*{20}{c}} {\mathop {\mathit{\boldsymbol{A}}}\limits_{m \times u} }\\ {\mathop {\mathit{\boldsymbol{a}}}\limits_{c \times u} } \end{array}} \right]\mathop {\mathit{\boldsymbol{X}}}\limits_{u \times 1} + \left[ {\begin{array}{*{20}{c}} {\mathop {\mathit{\boldsymbol{B}}}\limits_{m \times n} }&{\mathop {\bf{0}}\limits_{m \times c} }\\ {\mathop {\bf{0}}\limits_{c \times n} }&{\mathop { - \mathit{\boldsymbol{I}}}\limits_{c \times c} } \end{array}} \right]\left[ {\begin{array}{*{20}{c}} {\mathop {\mathit{\boldsymbol{V}}}\limits_{n \times 1} }\\ {\mathop {{\mathit{\boldsymbol{V}}_X}}\limits_{c \times 1} } \end{array}} \right] $

(4)

设$ \underset{\left( \mathit{m}+\mathit{c} \right)\times 1}{\mathop{{\mathit{\boldsymbol{\bar{L}}}}}}\,=\left[ \begin{align} & \underset{\mathit{m}\times 1}{\mathop{\mathit{\boldsymbol{L}}}}\, \\ & \underset{\mathit{c}\times 1}{\mathop{{{\mathit{\boldsymbol{ \mu}} \rm{ }}_{\mathit{X}}}}}\, \\ \end{align} \right], $ $ \underset{\left( \mathit{m}+\mathit{c} \right)\times \mathit{u}}{\mathop{{\mathit{\boldsymbol{\bar{A}}}}}}\,=\left[ \begin{align} & \underset{\mathit{m}\times \mathit{u}}{\mathop{\mathit{\boldsymbol{A}}}}\, \\ & \underset{\mathit{c}\times \mathit{u}}{\mathop{\mathit{\boldsymbol{a}}}}\, \\ \end{align} \right], $ $ \underset{\left( \mathit{m}+\mathit{c} \right)\times \left( \mathit{n}+\mathit{c} \right)}{\mathop{{\mathit{\boldsymbol{\bar{B}}}}}}\,=$ $ \left[ \begin{matrix} \underset{\mathit{m}\times \mathit{n}}{\mathop{\mathit{\boldsymbol{B}}}}\, & \underset{\mathit{m}\times \mathit{c}}{\mathop{\bf{0}}}\, \\ \underset{\mathit{c}\times \mathit{n}}{\mathop{\bf{0}}}\, & \underset{\mathit{c}\times \mathit{c}}{\mathop{\rm{-}\mathit{\boldsymbol{I}}}}\, \\ \end{matrix} \right],\underset{\left( \mathit{n}+\mathit{c} \right)\times 1}{\mathop{{\mathit{\boldsymbol{\bar{V}}}}}}\,=\left[ \begin{align} & \underset{\mathit{n}\times 1}{\mathop{\mathit{\boldsymbol{V}}}}\, \\ & \underset{\mathit{c}\times 1}{\mathop{{{\mathit{\boldsymbol{V}}}_{\mathit{X}}}}}\, \\ \end{align} \right] $，则：

$ \mathop {\mathit{\boldsymbol{\bar L}}}\limits_{\left( {m + c} \right) \times 1} = \mathop {\begin{array}{*{20}{c}} {\mathit{\boldsymbol{\bar A}}}&\mathit{\boldsymbol{X}} \end{array}}\limits_{\left( {m + c} \right) \times uu \times 1} + \mathop {\mathit{\boldsymbol{\bar B}}}\limits_{\left( {m + c} \right) \times \left( {n + c} \right)} \mathop {\mathit{\boldsymbol{\bar V}}}\limits_{\left( {n + c} \right) \times 1} $

(5)

随机模型E(V)=0，实际观测值L_V改正数V的方差为D(L_V)=δ₀²Q_VV，权为P_V=Q_VV^-1；X与V的协方差cov(X, V)=0，虚拟观测值μ_X的方差为D(X)=δ₀²Q_XX，其权P_X=Q_XX^-1。定权时采用同一单位权方差δ₀²，权矩阵为：

$ \mathop {\mathit{\boldsymbol{\bar P}}}\limits_{\left( {n + c} \right) \times \left( {n + c} \right)} = \left[ {\begin{array}{*{20}{c}} {{\mathit{\boldsymbol{P}}_V}}&0\\ 0&{{\mathit{\boldsymbol{P}}_X}} \end{array}} \right] $

(6)

1.2 算法推导

目标方程为：

$ {{\mathit{\boldsymbol{\bar V}}}^{\rm{T}}}\mathit{\boldsymbol{\bar P\bar V}} = \min ,{\rm{s}}.\;{\rm{t}}.\;\mathit{\boldsymbol{\bar L}} - \mathit{\boldsymbol{\bar AX}} - \mathit{\boldsymbol{\bar B\bar V}} = 0 $

(7)

构造拉格朗日极值方程：

$ \mathit{\Phi }\left( {\mathit{\boldsymbol{\bar V}},\mathit{\boldsymbol{X}},\mathit{\boldsymbol{K}}} \right) = {{\mathit{\boldsymbol{\bar V}}}^{\rm{T}}}\mathit{\boldsymbol{\bar P\bar V}} + 2{\mathit{\boldsymbol{K}}^{\rm{T}}}\left( {\mathit{\boldsymbol{\bar L}} - \mathit{\boldsymbol{\bar AX}} - \mathit{\boldsymbol{\bar B\bar V}}} \right) $

(8)

矩阵L和AX对V的偏导为0。由式(2)可知：

$ \mathit{\boldsymbol{\bar B\bar V}} = \left[ {\begin{array}{*{20}{c}} {\mathit{\boldsymbol{BV}}}\\ { - {\mathit{\boldsymbol{V}}_X}} \end{array}} \right] = \left[ {\begin{array}{*{20}{c}} {{\mathit{\boldsymbol{V}}_A}\mathit{\boldsymbol{X}} - \mathit{\boldsymbol{e}}}\\ { - {\mathit{\boldsymbol{V}}_X}} \end{array}} \right] $

(9)

故BV对X求偏导得$ \left[ \begin{align} &{{\mathit{\boldsymbol{V}}}_{\mathit{A}}} \\ &\ 0 \\ \end{align} \right]$，设

$ \mathit{\boldsymbol{\tilde A}} = \mathit{\boldsymbol{\bar A}} + \left[ {\begin{array}{*{20}{c}} {{\mathit{\boldsymbol{V}}_A}}\\ 0 \end{array}} \right] = \left[ {\begin{array}{*{20}{c}} {\mathit{\boldsymbol{A}} + {\mathit{\boldsymbol{V}}_A}}\\ \mathit{\boldsymbol{a}} \end{array}} \right] $

(10)

由Φ分别对V、X求偏导得：

即

$ \begin{array}{*{20}{c}} {\frac{{\partial \mathit{\Phi }}}{{\partial \mathit{\boldsymbol{\bar V}}}} = 2{{\mathit{\boldsymbol{\bar V}}}^{\rm{T}}}\mathit{\boldsymbol{\bar P}} - 2{\mathit{\boldsymbol{K}}^{\rm{T}}}\mathit{\boldsymbol{\bar B}} = 0,}\\ {\mathit{\boldsymbol{\bar P\bar V}} - {{\mathit{\boldsymbol{\bar B}}}^{\rm{T}}}\mathit{\boldsymbol{K}} = 0} \end{array} $

(11)

即

$ \begin{array}{*{20}{c}} {\frac{{\partial \mathit{\Phi }}}{{\partial \mathit{\boldsymbol{X}}}} = - 2{\mathit{\boldsymbol{K}}^{\rm{T}}}\left( {\mathit{\boldsymbol{\bar A}} + \left[ {\begin{array}{*{20}{c}} {{\mathit{\boldsymbol{V}}_A}}\\ 0 \end{array}} \right]} \right) = 0,}\\ {{{\mathit{\boldsymbol{\tilde A}}}^{\rm{T}}}\mathit{\boldsymbol{K}} = 0} \end{array} $

(12)

联合式(5)、(11)、(12)可得：

$ \left[ {\begin{array}{*{20}{c}} {\mathop {\mathit{\boldsymbol{\bar P}}}\limits_{\left( {n + c} \right) \times \left( {n + c} \right)} }&{\mathop {\bf{0}}\limits_{\left( {n + c} \right) \times u} }&{\mathop { - {\mathit{\boldsymbol{\bar B}}^{\rm{T}}}}\limits_{\left( {n + c} \right) \times \left( {m + c} \right)} }\\ {\mathop {\bf{0}}\limits_{u \times \left( {n + c} \right)} }&{\mathop {\bf{0}}\limits_{u \times u} }&{\mathop {{{\mathit{\boldsymbol{\tilde A}}}^{\rm{T}}}}\limits_{u \times \left( {m + c} \right)} }\\ {\mathop {\mathit{\boldsymbol{\bar B}}}\limits_{\left( {m + c} \right) \times \left( {n + c} \right)} }&{\mathop {\mathit{\boldsymbol{\bar A}}}\limits_{\left( {m + c} \right) \times u} }&{\mathop {\bf{0}}\limits_{\left( {m + c} \right) \times \left( {m + c} \right)} } \end{array}} \right]\left[ {\begin{array}{*{20}{c}} {\mathop {\mathit{\boldsymbol{\bar V}}}\limits_{\left( {n + c} \right) \times 1} }\\ {\mathop {\mathit{\boldsymbol{X}}}\limits_{u \times 1} }\\ {\mathop {\mathit{\boldsymbol{K}}}\limits_{\left( {m + c} \right) \times 1} } \end{array}} \right] = \left[ {\begin{array}{*{20}{c}} {\mathop {\bf{0}}\limits_{\left( {n + c} \right) \times 1} }\\ {\mathop {\bf{0}}\limits_{u \times 1} }\\ {\mathop {\mathit{\boldsymbol{\bar L}}}\limits_{\left( {m + c} \right) \times 1} } \end{array}} \right] $

(13)

待估参数X和改正数向量V的解为：

$ \left[ {\begin{array}{*{20}{c}} {\mathop {\mathit{\boldsymbol{\hat {\bar V}}}}\limits_{\left( {n + c} \right) \times 1} }\\ {\mathop {\mathit{\boldsymbol{\hat X}}}\limits_{u \times 1} }\\ {\mathop {\mathit{\boldsymbol{\hat K}}}\limits_{\left( {m + c} \right) \times 1} } \end{array}} \right] = {\left[ {\begin{array}{*{20}{c}} {\mathop {\mathit{\boldsymbol{\bar P}}}\limits_{\left( {n + c} \right) \times \left( {n + c} \right)} }&{\mathop {\bf{0}}\limits_{\left( {n + c} \right) \times u} }&{\mathop { - {{\mathit{\boldsymbol{\bar B}}}^{\rm{T}}}}\limits_{\left( {n + c} \right) \times \left( {m + c} \right)} }\\ {\mathop {\bf{0}}\limits_{u \times \left( {n + c} \right)} }&{\mathop {\bf{0}}\limits_{u \times u} }&{\mathop {{{\mathit{\boldsymbol{\tilde A}}}^{\rm{T}}}}\limits_{u \times \left( {m + c} \right)} }\\ {\mathop {\mathit{\boldsymbol{\bar B}}}\limits_{\left( {m + c} \right) \times \left( {n + c} \right)} }&{\mathop {\mathit{\boldsymbol{\bar A}}}\limits_{\left( {m + c} \right) \times u} }&{\mathop {\bf{0}}\limits_{\left( {m + c} \right) \times \left( {m + c} \right)} } \end{array}} \right]^{ - 1}}\left[ {\begin{array}{*{20}{c}} {\mathop {\bf{0}}\limits_{\left( {n + c} \right) \times 1} }\\ {\mathop {\bf{0}}\limits_{u \times 1} }\\ {\mathop {\mathit{\boldsymbol{\bar L}}}\limits_{\left( {m + c} \right) \times 1} } \end{array}} \right] $

(14)

式中，等号两侧都含有V、X，可设定V初值为0，X初值为X₀迭代求解。式(14)需对n+m+2c+u阶矩阵求逆，且V、X、K位于同一矩阵中。由式(11)可得:

$ \mathit{\boldsymbol{\bar V}} = {{\mathit{\boldsymbol{\bar P}}}^{ - 1}}{{\mathit{\boldsymbol{\bar B}}}^{\rm{T}}}\mathit{\boldsymbol{K}} $

(15)

将式(15)代入式(5)得：

$ \mathit{\boldsymbol{\bar L}} = \mathit{\boldsymbol{\bar AX}} + \mathit{\boldsymbol{\bar B}}{{\mathit{\boldsymbol{\bar P}}}^{ - 1}}{{\mathit{\boldsymbol{\bar B}}}^{\rm{T}}}\mathit{\boldsymbol{K}} $

(16)

设$ \mathit{\boldsymbol{\tilde{P}}}\text{=}{{\left( {{\overline{\mathit{\boldsymbol{B}}}}{\overline{\mathit{\boldsymbol{P}}}}^{\text{-}1}}{{\overline{\mathit{\boldsymbol{B}}}}^{\text{T}}} \right)}^{\text{-}1}}$，则：

$ \mathit{\boldsymbol{K}} = {\left( {\mathit{\boldsymbol{\bar B}}{{\mathit{\boldsymbol{\bar P}}}^{ - 1}}{{\mathit{\boldsymbol{\bar B}}}^{\rm{T}}}} \right)^{ - 1}}\left( {\mathit{\boldsymbol{\bar L}} - \mathit{\boldsymbol{\bar AX}}} \right) = \mathit{\boldsymbol{\tilde P}}\left( {\mathit{\boldsymbol{\bar L}} - \mathit{\boldsymbol{\bar AX}}} \right) $

(17)

将式(17)代入式(12)得：

$ {{\mathit{\boldsymbol{\tilde A}}}^{\rm{T}}}\mathit{\boldsymbol{\tilde P}}\left( {\mathit{\boldsymbol{\bar L}} - \mathit{\boldsymbol{\bar AX}}} \right) = 0 $

(18a)

设$ \mathit{\boldsymbol{\tilde{L}}}=\left[ \begin{matrix} \mathit{\boldsymbol{L}}\text{+}{{\mathit{\boldsymbol{V}}}_{\mathit{A}}}\mathit{\boldsymbol{X}} \\ {\mathit{\boldsymbol{\mu }}_{\mathit{X}}} \\ \end{matrix} \right]$，则：

$ {{\mathit{\boldsymbol{\tilde A}}}^{\rm{T}}}\mathit{\boldsymbol{\tilde P}}\left( {\mathit{\boldsymbol{\tilde L}} - \mathit{\boldsymbol{\tilde AX}}} \right) = 0 $

(18b)

由式(18)可得：

$ \mathit{\boldsymbol{\hat X}} = {\left( {{{\mathit{\boldsymbol{\tilde A}}}^{\rm{T}}}\mathit{\boldsymbol{\tilde P\bar A}}} \right)^{ - 1}}\left( {{{\mathit{\boldsymbol{\tilde A}}}^{\rm{T}}}\mathit{\boldsymbol{\tilde P\bar L}}} \right) $

(19a)

或

$ \mathit{\boldsymbol{\hat X}} = {\left( {{{\mathit{\boldsymbol{\tilde A}}}^{\rm{T}}}\mathit{\boldsymbol{\tilde P\tilde A}}} \right)^{ - 1}}\left( {{{\mathit{\boldsymbol{\tilde A}}}^{\rm{T}}}\mathit{\boldsymbol{\tilde P\tilde L}}} \right) $

(19b)

求得X后，由式(17)求得矩阵K，并代入式(15)求得矩阵V为：

$ \mathit{\boldsymbol{\hat {\bar V}}} = {{\mathit{\boldsymbol{\bar P}}}^{ - 1}}{{\mathit{\boldsymbol{\bar B}}}^{\rm{T}}}\mathit{\boldsymbol{\tilde P}}\left( {\mathit{\boldsymbol{\bar L}} - \mathit{\boldsymbol{\bar A\hat X}}} \right) $

(20)

1.3 解的另一形式(a=I时)

展开$ {\mathit{\boldsymbol{\tilde{P}}}}$的表达式：

$ \begin{array}{*{20}{c}} {\mathit{\boldsymbol{\tilde P}} = {{\left( {\mathit{\boldsymbol{\bar B}}{{\mathit{\boldsymbol{\bar P}}}^{ - 1}}{{\mathit{\boldsymbol{\bar B}}}^{\rm{T}}}} \right)}^{ - 1}} = }\\ {{{\left( {\left[ {\begin{array}{*{20}{c}} \mathit{\boldsymbol{B}}&0\\ 0&{ - \mathit{\boldsymbol{I}}} \end{array}} \right]\left[ {\begin{array}{*{20}{c}} {{\mathit{\boldsymbol{Q}}_{VV}}}&0\\ 0&{{\mathit{\boldsymbol{Q}}_{XX}}} \end{array}} \right]\left[ {\begin{array}{*{20}{c}} {{\mathit{\boldsymbol{B}}^{\rm{T}}}}&0\\ 0&{ - \mathit{\boldsymbol{I}}} \end{array}} \right]} \right)}^{ - 1}} = }\\ {\left[ {\begin{array}{*{20}{c}} {{{\left( {\mathit{\boldsymbol{B}}{\mathit{\boldsymbol{Q}}_{VV}}{\mathit{\boldsymbol{B}}^{\rm{T}}}} \right)}^{ - 1}}}&0\\ 0&{\mathit{\boldsymbol{Q}}_{XX}^{ - 1}} \end{array}} \right]} \end{array} $

(21)

结合式(5)、(10)及(19)得：

$ \begin{array}{*{20}{c}} {\mathit{\boldsymbol{\hat X}} = {{\left( {{{\mathit{\boldsymbol{\tilde A}}}^{\rm{T}}}\mathit{\boldsymbol{\tilde P\tilde A}}} \right)}^{ - 1}}\left( {{{\mathit{\boldsymbol{\tilde A}}}^{\rm{T}}}\mathit{\boldsymbol{\tilde P\tilde L}}} \right) = }\\ {{{\left( {\left[ {\begin{array}{*{20}{c}} {{{\left( {\mathit{\boldsymbol{A}} + {\mathit{\boldsymbol{V}}_A}} \right)}^{\rm{T}}}}&{{\mathit{\boldsymbol{a}}^{\rm{T}}}} \end{array}} \right]\left[ {\begin{array}{*{20}{c}} {{{\left( {\mathit{\boldsymbol{B}}{\mathit{\boldsymbol{Q}}_{VV}}{\mathit{\boldsymbol{B}}^{\rm{T}}}} \right)}^{ - 1}}}&0\\ 0&{\mathit{\boldsymbol{Q}}_{XX}^{ - 1}} \end{array}} \right]\left[ {\begin{array}{*{20}{c}} \mathit{\boldsymbol{A}}\\ \mathit{\boldsymbol{a}} \end{array}} \right]} \right)}^{ - 1}} \cdot }\\ {\left( {\left[ {\begin{array}{*{20}{c}} {{{\left( {\mathit{\boldsymbol{A}} + {\mathit{\boldsymbol{V}}_A}} \right)}^{\rm{T}}}}&{{\mathit{\boldsymbol{a}}^{\rm{T}}}} \end{array}} \right]\left[ {\begin{array}{*{20}{c}} {{{\left( {\mathit{\boldsymbol{B}}{\mathit{\boldsymbol{Q}}_{VV}}{\mathit{\boldsymbol{B}}^{\rm{T}}}} \right)}^{ - 1}}}&0\\ 0&{\mathit{\boldsymbol{Q}}_{XX}^{ - 1}} \end{array}} \right]\left[ {\begin{array}{*{20}{c}} \mathit{\boldsymbol{L}}\\ {{\mathit{\boldsymbol{\mu }}_X}} \end{array}} \right]} \right) = }\\ {{{\left[ {{{\left( {\mathit{\boldsymbol{A}} + {\mathit{\boldsymbol{V}}_A}} \right)}^{\rm{T}}}{{\left( {\mathit{\boldsymbol{B}}{\mathit{\boldsymbol{Q}}_{VV}}{\mathit{\boldsymbol{B}}^{\rm{T}}}} \right)}^{ - 1}}\mathit{\boldsymbol{A}} + {\mathit{\boldsymbol{a}}^{\rm{T}}}\mathit{\boldsymbol{Q}}_{XX}^{ - 1}\mathit{\boldsymbol{a}}} \right]}^{ - 1}} \cdot }\\ {\left[ {{{\left( {\mathit{\boldsymbol{A}} + {\mathit{\boldsymbol{V}}_A}} \right)}^{\rm{T}}}{{\left( {\mathit{\boldsymbol{B}}{\mathit{\boldsymbol{Q}}_{VV}}{\mathit{\boldsymbol{B}}^{\rm{T}}}} \right)}^{ - 1}}\mathit{\boldsymbol{L}} + {\mathit{\boldsymbol{a}}^{\rm{T}}}\mathit{\boldsymbol{Q}}_{XX}^{ - 1}{\mathit{\boldsymbol{\mu }}_X}} \right]} \end{array} $

(22)

由于c×u阶矩阵a的秩为c，且c≤u，附加部分随机参数时，a^TQ_XX^-1a秩亏不可逆，仅当a=I时满秩可逆，因此以下公式推导仅基于a=I的情形。

由矩阵反演公式(BD^-1A + C^-1)^-1BD^-1=CB(D + ACB)^-1可得：

$ \begin{array}{*{20}{c}} {{{\left[ {{{\left( {\mathit{\boldsymbol{A}} + {\mathit{\boldsymbol{V}}_A}} \right)}^{\rm{T}}}{{\left( {\mathit{\boldsymbol{B}}{\mathit{\boldsymbol{Q}}_{VV}}{\mathit{\boldsymbol{B}}^{\rm{T}}}} \right)}^{ - 1}}\mathit{\boldsymbol{A}} + \mathit{\boldsymbol{Q}}_{XX}^{ - 1}} \right]}^{ - 1}}}\\ {{{\left( {\mathit{\boldsymbol{A}} + {\mathit{\boldsymbol{V}}_A}} \right)}^{\rm{T}}}{{\left( {\mathit{\boldsymbol{B}}{\mathit{\boldsymbol{Q}}_{VV}}{\mathit{\boldsymbol{B}}^{\rm{T}}}} \right)}^{ - 1}} = {\mathit{\boldsymbol{Q}}_{XX}}{{\left( {\mathit{\boldsymbol{A}} + {\mathit{\boldsymbol{V}}_A}} \right)}^{\rm{T}}} \cdot }\\ {{{\left( {\mathit{\boldsymbol{B}}{\mathit{\boldsymbol{Q}}_{VV}}{\mathit{\boldsymbol{B}}^{\rm{T}}} + \mathit{\boldsymbol{A}}{\mathit{\boldsymbol{Q}}_{XX}}{{\left( {\mathit{\boldsymbol{A}} + {\mathit{\boldsymbol{V}}_A}} \right)}^{\rm{T}}}} \right)}^{ - 1}}} \end{array} $

(23a)

由矩阵反演公式(ACB + D)^-1=D^-1- D^-1A ·(C^-1+ BD^-1A)^-1BD^-1可得：

$ \begin{array}{*{20}{c}} {{{\left[ {{{\left( {\mathit{\boldsymbol{A}} + {\mathit{\boldsymbol{V}}_A}} \right)}^{\rm{T}}}{{\left( {\mathit{\boldsymbol{B}}{\mathit{\boldsymbol{Q}}_{VV}}{\mathit{\boldsymbol{B}}^{\rm{T}}}} \right)}^{ - 1}}\mathit{\boldsymbol{A}} + \mathit{\boldsymbol{Q}}_{XX}^{ - 1}} \right]}^{ - 1}} = }\\ {{\mathit{\boldsymbol{Q}}_{XX}} - {\mathit{\boldsymbol{Q}}_{XX}}{{\left( {\mathit{\boldsymbol{A}} + {\mathit{\boldsymbol{V}}_A}} \right)}^{\rm{T}}} \cdot }\\ {{{\left( {\mathit{\boldsymbol{B}}{\mathit{\boldsymbol{Q}}_{VV}}{\mathit{\boldsymbol{B}}^{\rm{T}}} + \mathit{\boldsymbol{A}}{\mathit{\boldsymbol{Q}}_{XX}}{{\left( {\mathit{\boldsymbol{A}} + {\mathit{\boldsymbol{V}}_A}} \right)}^{\rm{T}}}} \right)}^{ - 1}}\mathit{\boldsymbol{A}}{\mathit{\boldsymbol{Q}}_{XX}}} \end{array} $

(23b)

结合式(22)、(23)得：

$ \begin{array}{*{20}{c}} {\mathit{\boldsymbol{\hat X}}\left| {_{a = I}} \right. = {\mathit{\boldsymbol{\mu }}_X} + {\mathit{\boldsymbol{Q}}_{XX}}{{\left( {\mathit{\boldsymbol{A}} + {\mathit{\boldsymbol{V}}_A}} \right)}^{\rm{T}}}}\\ {{{\left( {\mathit{\boldsymbol{B}}{\mathit{\boldsymbol{Q}}_{VV}}{\mathit{\boldsymbol{B}}^{\rm{T}}} + \mathit{\boldsymbol{A}}{\mathit{\boldsymbol{Q}}_{XX}}{{\left( {\mathit{\boldsymbol{A}} + {\mathit{\boldsymbol{V}}_A}} \right)}^{\rm{T}}}} \right)}^{ - 1}}\left( {\mathit{\boldsymbol{L}} - \mathit{\boldsymbol{A}}{\mathit{\boldsymbol{\mu }}_X}} \right)} \end{array} $

(24a)

$ \begin{array}{*{20}{c}} {\mathit{\boldsymbol{\hat X}}\left| {_{a = I}} \right. = {\mathit{\boldsymbol{\mu }}_X} + {{\left[ {{{\left( {\mathit{\boldsymbol{A}} + {\mathit{\boldsymbol{V}}_A}} \right)}^{\rm{T}}}{{\left( {\mathit{\boldsymbol{B}}{\mathit{\boldsymbol{Q}}_{VV}}{\mathit{\boldsymbol{B}}^{\rm{T}}}} \right)}^{ - 1}}\mathit{\boldsymbol{A + Q}}_{XX}^{ - 1}} \right]}^{ - 1}} \cdot }\\ {{{\left( {\mathit{\boldsymbol{A}} + {\mathit{\boldsymbol{V}}_A}} \right)}^{\rm{T}}}{{\left( {\mathit{\boldsymbol{B}}{\mathit{\boldsymbol{Q}}_{VV}}{\mathit{\boldsymbol{B}}^{\rm{T}}}} \right)}^{ - 1}}\left( {\mathit{\boldsymbol{L}} - \mathit{\boldsymbol{A}}{\mathit{\boldsymbol{\mu }}_X}} \right)} \end{array} $

(24b)

式(24a)适合在Q_XX奇异时采用，式(24b)要求Q_XX^-1存在。当a=I时，式(24)与极大验后估计及卡尔曼滤波公式形式上一致，而式(14)与式(19)更具有一般适用性。

1.4 精度评定

在逐次动态滤波中，参数X和方差D_X都是下一轮计算的基础，准确估计方差D_X与参数X同等重要。整体最小二乘平差计算单位权方差时，V采用矩阵A、观测值L和虚拟观测μ_X中的全部观测改正值。由于增加了c个虚拟观测值，必要观测数是u，因此平差问题的自由度为r=m+c-u，m为观测向量L的阶数。单位权方差δ₀²为：

$ \delta _0^2 = \frac{{{{\mathit{\boldsymbol{\bar V}}}^{\rm{T}}}\mathit{\boldsymbol{\bar P\bar V}}}}{{m + c - u}} $

(25)

式中，V为平差问题中所有观测值的改正数向量，由式(20)计算，阶数为n+c，P是每一个观测值对应的n+c阶权矩阵，见式(6)。

严密的协方差阵计算需考虑WTLS算式的非线性特征，基于误差传播定律，采用迭代计算^[13]或无味变换^[14]来实现，但过程复杂，限制了其应用，以下给出待估参数协方差的近似公式。近似公式忽略了V_A的随机性，由式(19b)根据协因数传播定律，可得平差待估参数X的近似协方差阵$ {{\mathit{\boldsymbol{D}}}_{\mathit{\hat{X}\hat{X}}}}$为：

$ \begin{array}{*{20}{c}} {{\mathit{\boldsymbol{D}}_{\hat X\hat X}} = \delta _0^2{\mathit{\boldsymbol{Q}}_{\hat X\hat X}} = \delta _0^2{{\left( {{{\mathit{\boldsymbol{\tilde A}}}^{\rm{T}}}\mathit{\boldsymbol{\tilde P\tilde A}}} \right)}^{ - 1}} \cdot }\\ {{{\mathit{\boldsymbol{\tilde A}}}^{\rm{T}}}\mathit{\boldsymbol{\tilde P}}{\mathit{\boldsymbol{Q}}_{\tilde L\tilde L}}\mathit{\boldsymbol{\tilde P\tilde A}}{{\left( {{{\mathit{\boldsymbol{\tilde A}}}^{\rm{T}}}\mathit{\boldsymbol{\tilde P\tilde A}}} \right)}^{ - 1}} = \delta _0^2{{\left( {{{\mathit{\boldsymbol{\tilde A}}}^{\rm{T}}}\mathit{\boldsymbol{\tilde P\tilde A}}} \right)}^{ - 1}}} \end{array} $

(26)

1.5 计算步骤

迭代算法步骤如下：1)根据观测方程线性化，形成平差函数模型(5)；2)定权P，根据经典LS计算未知参数初值X₀，将X的初值赋给矩阵A、B中的相应元素，设改正数向量V的初值为0，矩阵$ {\mathit{\boldsymbol{\tilde{A}}}}$的初值为A；3)将B、P、A、$ {\mathit{\boldsymbol{\tilde{A}}}}$、L代入式(14)或式(19)、(20)或式(22)、(24)中，求得未知参数X和改正数V的值；4)把未知参数X和改正数V的新值用于更新算式中对应的更新项；5)重复步骤3)~ 4)，直至相邻两次V和X的值之差小于阈值ε或迭代达到指定次数时结束；6)采用式(25)~(26)近似公式进行精度评定。

2 实例分析 2.1 直线拟合实验1

随机参数X的先验期望与先验方差由经典LS平差计算。先验期望与方差为：

$ {\mathit{\boldsymbol{\mu }}_X} = {\left[ {\begin{array}{*{20}{c}} { - 0.610\;812\;957}&{6.100\;109\;317} \end{array}} \right]^{\rm{T}}} $

$ {\mathit{\boldsymbol{D}}_X} = {10^S} \cdot \left[ {\begin{array}{*{20}{c}} {0.003\;886\;394\;5}&{ - 0.026\;036\;202\;9}\\ { - 0.026\;036\;202\;9}&{0.179\;826\;418\;9} \end{array}} \right] $

其中，S=-10~10。由直线方程y=kx+b对观测值加改正数得y+v_y=k(x+v_x)+b，由此构造如式(5)的平差数学模型，即

$ y = \left[ {\begin{array}{*{20}{c}} x&1 \end{array}} \right]\left[ {\begin{array}{*{20}{c}} k\\ b \end{array}} \right] + \left[ {\begin{array}{*{20}{c}} k&{ - 1} \end{array}} \right]\left[ {\begin{array}{*{20}{c}} {{v_x}}\\ {{v_y}} \end{array}} \right] $

(27)

设R=[k -1]，B=I_i⊗R为i阶单位矩阵I_i与R的Kronecker张量积，收敛阈值ε设为10^-16。实验数据见表 1^[15]，k、b参数估计见表 2，精度评定见表 3。

表 2结果表明，2种算法的参数估计和单位权方差结果一致，但本文算法的迭代次数N较少。计算δ₀²时，V采用了所有x_i、y_i和μ_X改正数，平差自由度为10。当S=-10时，先验方差几乎为0，先验期望权很大，参数估计结果就是μ_X；当S=10时，先验方差很大，先验期望权接近0，μ_X失去作用，参数估计结果就是WTLS估计的结果；当S=0时，参数估计结果介于两者之间。

当X不具有随机性，即无先验期望与方差时(c=0)，单位权方差δ₀²=1.483 294 149 3，略小于表 2中S=0时的δ₀²=1.500 705 449 8。由表 3可知，参数X的先验信息参与平差(c=2)提高了参数估计的精度。

2.2 直线拟合实验2

本实例数据同表 1，将其中x_i坐标的权进行放大p_xi→10^S· p_xi，S=11~0。先验数据为:

$ {\mathit{\boldsymbol{\mu }}_X} = {\left[ {\begin{array}{*{20}{c}} { - 0.610\;812\;95\;7}&{6.100\;109\;317} \end{array}} \right]^{\rm{T}}}, $

$ \begin{array}{*{20}{c}} {{\mathit{\boldsymbol{D}}_X} = {{10}^{12}} \cdot }\\ {\left[ {\begin{array}{*{20}{c}} {0.003\;886\;394\;5}&{ - 0.026\;036\;202\;9}\\ { - 0.026\;036\;202\;9}&{0.179\;826\;418\;9} \end{array}} \right]} \end{array} $

平差问题的自由度为10，参数估计及单位权方差估计结果见表 4。

表 4显示，2种方法参数估计结果一致，结果符合常理，当先验方差数值非常大时先验权降低，先验期望几乎失去作用。当S=11时，x_i坐标的权非常大，可以认为x_i坐标几乎没有误差，由于y_i坐标的权未改变，此时得到的解与视x_i坐标无误差的经典最小二乘解一致；当S=0时，x_i坐标的权未放大，平差结果与未附加随机参数的WTLS平差结果一致。

2.3 圆曲线拟合实验

本例数据来源于文献[16]，用该数据拟合最佳圆曲线，求圆心坐标与圆半径，观测坐标见表 5。

由点到圆心的距离与半径之差平方和最小，得目标方程为：

$ \begin{array}{*{20}{c}} {{d^2} = }\\ {{{\left\| {\sqrt {{{\left( {x - {x_0}} \right)}^2} + {{\left( {y - {y_0}} \right)}^2}} - r} \right\|}_{{x_0},{y_0},r}} = \min } \end{array} $

(28)

泰勒公式展开得：

$ \begin{array}{*{20}{c}} {\sqrt {{{\left( {x - x_0^0} \right)}^2} + {{\left( {y - y_0^0} \right)}^2}} - {r_0} + }\\ {\frac{{\left( {x - x_0^0} \right)\left( {{v_x} - {\delta _x}} \right) + \left( {y - y_0^0} \right)\left( {{v_y} - {\delta _y}} \right)}}{{\sqrt {{{\left( {x - x_0^0} \right)}^2} + {{\left( {y - y_0^0} \right)}^2}} }} - {\delta _r} = 0} \end{array} $

(29)

设$ {{\mathit{d}}_{\text{0}}}=\sqrt{{{\left( x-x_{0}^{0} \right)}^{2}}+{{\left( y-y_{0}^{0} \right)}^{2}}}$，构造如式(30)的函数模型：

$ \begin{array}{*{20}{c}} {{d_0} - {r_0} = \left[ {\begin{array}{*{20}{c}} {\frac{{x - x_0^0}}{{{d_0}}}}&{\frac{{y - y_0^0}}{{{d_0}}}}&1 \end{array}} \right]\left[ {\begin{array}{*{20}{c}} {{\delta _x}}\\ {{\delta _y}}\\ {{\delta _r}} \end{array}} \right] + }\\ {\left[ {\begin{array}{*{20}{c}} {\frac{{x_0^0 - x}}{{{d_0}}}}&{\frac{{y_0^0 - y}}{{{d_0}}}} \end{array}} \right]\left[ {\begin{array}{*{20}{c}} {{v_x}}\\ {{v_y}} \end{array}} \right]} \end{array} $

(30)

设l=d₀-r₀为观测值，经典LS平差时视l为观测值，改正数为v_l；而WTLS平差时，视X与Y坐标为观测值，改正数为v_x、v_y。由式(30)可知，l中包含了观测值X、Y和未知参数，v_l中融合了改正数v_x、v_y，并非只对Y改正而认为X没有误差。因此本例在X、Y等权的情形下，LS解与WTLS解等价，只是2种方法得到的改正数不同。WTLS解更具一般性，更适合处理X、Y不等权的情形。

在满足IID数据及权矩阵为单位矩阵等权的情况下，Schaffrin等^[17]的算法结果及本文WTLS参数估计结果见表 6。

表 6显示，2种方法的参数估计与精度评定结果一致，但本文算法在对初值的要求、迭代次数、收敛阈值等方面更具优势。由于每个点到拟合圆偏离的距离不同，即它们是不等权的，对X、Y坐标附加权，权由点到拟合圆的距离远近确定，权数据见表 7。

当X与Y的权相等时，由权倒数传播律可计算出l=d₀-r₀的权为表 7中X或Y坐标对应的权；反之X与Y的权不相等时则不成立。实验分3种情况：不加先验信息；加1个先验信息r；加3个先验信息x₀、y₀、r。先验信息均采用LS近似解。设3次加权实验中的初值均为x₀⁰=5，y₀⁰=5，r₀=5，收敛阈值均为1×10¹⁰，采用式(19b)计算，结果见表 8。

若仅半径r为随机参数，μ_r=4.714，D_r=4×10⁸，圆心坐标没有先验信息，则：

$ {\mu _r} - {r_0} = \left[ {\begin{array}{*{20}{c}} 0&0&1 \end{array}} \right]\left[ {\begin{array}{*{20}{c}} {{\delta _x}}\\ {{\delta _y}}\\ {{\delta _r}} \end{array}} \right] - {v_r} $

(31a)

若3个参数均为随机参数，μ_X=[4.739 2.983 4.714]^T，P_X=diag([2.5 2.5 2.5])，则：

$ \left[ {\begin{array}{*{20}{c}} {{\mu _x}}\\ {{\mu _y}}\\ {{\mu _r}} \end{array}} \right] - \left[ {\begin{array}{*{20}{c}} {x_0^0}\\ {y_0^0}\\ {{r_0}} \end{array}} \right] = \left[ {\begin{array}{*{20}{c}} 1&0&0\\ 0&1&0\\ 0&0&1 \end{array}} \right]\left[ {\begin{array}{*{20}{c}} {{\delta _x}}\\ {{\delta _y}}\\ {{\delta _r}} \end{array}} \right] - \left[ {\begin{array}{*{20}{c}} {{v_{{x_0}}}}\\ {{v_{{y_0}}}}\\ {{v_r}} \end{array}} \right] $

(31b)

表 8显示，参数具有先验随机性可提高待估参数的精度，加1个先验信息r时设置其方差非常小(4×10^-8)，估计结果的精度甚至优于附加3个先验信息的情形，说明先验信息精度越高则待估参数的精度越高。

3 结语

EIV-REM模型只能处理全部未知参数都具有先验随机信息的情形，本文提出部分参数具有先验随机性的加权整体最小二乘平差函数模型及其算法，推导了部分参数具有先验随机性的加权整体最小二乘平差参数估计和精度评定近似公式。分析表明，EIV-REM模型算法是本文模型算法在a=I时的特例，与EIV-REM模型相比，本文模型算法的优势在于其普适性，可处理部分或全部参数具有先验随机性的情形。不同算例解算结果显示，本文算法结果在特殊情形下可与加权整体最小二乘及经典最小二乘等算法一致，附加的先验信息精度越高则估计结果精度越高，算法可靠，结果正确，且迭代收敛速度比EIV-REM模型算法更优。