大角度三维基准转换的粗差探测算法

戴鹏洋，王彬; DAI Pengyang; WANG Bin

2025年4月29日 22:26 星期二

网刊加载中。。。

使用Chrome浏览器效果最佳，继续浏览，你可能不会看到最佳的展示效果，

确定继续浏览么?

复制成功，请在其他浏览器进行阅读

大角度三维基准转换的粗差探测算法 PDF

- ORCID：
戴鹏洋 ¹
- ORCID：
王彬 ²
✉

1. 同济大学测绘与地理信息学院，上海 200092； 2. 南京工业大学测绘科学与技术学院，南京 211800

中图分类号： P207

最近更新：2024-03-13

DOI：10.16356/j.1005-2615.2024.01.009

摘要

三维基准转换广泛应用于大地测量、摄影测量、点云配准等领域，求解大角度、任意比例尺的三维基准转换参数的研究有很多。然而，当观测值中含有粗差时，得到的转换参数估值会受到不利影响甚至被严重扭曲。为处理含有粗差的大角度三维基准转换问题，本文首先将大角度三维基准转换问题抽象为具有等式约束的最小二乘问题（Constrained least squares， CLS），推导参数在正交约束条件下的最小二乘解。然后，将灵敏度分析方法应用到CLS问题中，研究残差加权平方和对观测值扰动的局部敏感性，并基于这些敏感度指标构造局部检验统计量，进而推导出一个适用于CLS问题的粗差探测算法。最后，为核实该算法的有效性进行了仿真与实测数据实验。实验结果表明：本文提出的基于灵敏度检验统计量的数据探测算法可以降低粗差的负面影响，得到可靠的参数估值，从而有效解决大角度三维基准转换中的粗差处理问题。

关键词

大角度三维基准转换; 约束最小二乘; 粗差; 数据探测算法

三维基准转换在测绘科学领域广泛应用。空间三维坐标转换的模型中含有7个独立的参数，其中包含3个旋转角度、3个平移参数和1个尺度参数。基准转换的核心在于利用公共点的两套坐标求解转换参数，之后对非公共点进行基准转换。

经典的Bursa‑Wolf模型^［

1］和Molodensky模型^{［参考文献 2

百度学术}2］都是在小角度的基础上进行的简化，用于求取小角度情况下的坐标转换参数。随着测量技术的发展，在摄影测量领域、点云配准等方面都无法避免求解任意旋转角度下的三维基准转换参数，其中不乏大旋转角度的情形。为解决此类问题，文献［3］在七参数的初值处按照泰勒级数展开，仅保留一次项，将误差方程线性化基于最小二乘准则进行参数估计迭代计算控制舍入误差。在进行线性化时不可避免地要进行复杂的偏导数计算，直接解算非常复杂。文献［4］通过将3个平移参数、1个尺度参数、旋转矩阵中的9个方向余弦共13个参数视为未知数，在这13个参数的初始值处做线性化处理，旋转矩阵中的9个方向余弦中仅有3个独立参数，利用旋转矩阵是一个正交矩阵这一先验约束条件通过附有限制条件的间接平差法解算参数。文献［5］用四元数的概念来表示基准参数，在进行计算时不需要线性化和迭代，进一步拓展该算法可以在约束条件下计算任意旋转角度的基于四元数的基准转换参数，但该方法只适用于等精度的条件。文献［6］利用多元模型的矩阵形式，将多点坐标组成矩阵并顾及旋转矩阵的正交性导出了大角度三维基准转换的解析分步解，与传统的迭代解有等效转换结果，同时也避免了复杂的迭代计算过程。然而，以上所提出的算法均未考虑到观测值中含有粗差的情况。当观测值中含有粗差时，最小二乘方法缺乏稳健性，参数估计结果不理想。Lehmann^{［参考文献 7

百度学术}7］指出异常值是可能由观测值含有粗差造成的错误，所以含有粗差的观测值不应参与到平差的解算过程。附有等式约束的最小二乘问题（Constrained least squares， CLS）同样会受到粗差的不利影响，因此有必要研究针对CLS问题的粗差处理方法。基于正态分布线性模型的敏感度分析理论自Cook^{［参考文献 8

百度学术}8］提出之后得到了广泛研究。Guo等^{［参考文献 9

百度学术}9］详细讨论了随机模型扰动引起的LS残差加权平方和的局部敏感性，利用矩阵微分法详细推导了敏感度指标，并基于这些指标进行假设检验从而达到剔除粗差的目的。Wang等^{［参考文献 10

百度学术}10］将广义EIV（Errors‑in‑variable）模型的加权整体最小二乘解重新表达成经典最小二乘问题，依据最小二乘理论推导出检验统计量进行粗差探测。针对转换模型中含有等式约束的问题，Wang等^{［参考文献 11

百度学术}11］进一步提出了具有等式约束的非线性高斯赫尔默特模型的数据探测算法。

本文在文献［

4］基础上，将尺度参数与旋转矩阵相乘得到新的旋转矩阵，将其中的9个元素以及3个平移参数视为未知数构建十二参数模型，相比七参数模型简化了线性化的计算过程，加快了计算时间。利用旋转矩阵的正交性推导出新的约束条件，采用附有等式约束的最小二乘法求解未知参数。借鉴文献［11］的思路，推导了CLS模型的粗差探测算法，用于本文采用的十二参数模型中，达到抵抗粗差的目的。

1 三维坐标转换的十二参数模型及其解法

公共点坐标转换的数学模型表示为^［

4］

[\begin{matrix} x_{2 i} \\ y_{2 i} \\ z_{2 i} \end{matrix}] - [\begin{matrix} e_{x 2 i} \\ e_{y 2 i} \\ e_{z 2 i} \end{matrix}] = [\begin{matrix} x_{1 i} \\ y_{1 i} \\ z_{1 i} \end{matrix}] Ξ + [\begin{matrix} Δ x \\ Δ y \\ Δ z \end{matrix}]

（1）

式中： $[x_{1 i} y_{1 i} z_{1 i}]^{T}$ 表示第i个公共点在原坐标系下的观测向量； $[x_{2 i} y_{2 i} z_{2 i}]^{T}$ 表示第i个点在目标坐标系下的观测向量； $[e_{x 2 i} e_{y 2 i} e_{z 2 i}]^{T}$ 代表 $[x_{2 i} y_{2 i} z_{2 i}]^{T}$ 的随机误差向量；等式右端的 ${[Δ x Δ y Δ z]}^{T}$ 代表3个平移参数； $Ξ$ 为3×3阶的矩阵，其中9个元素均视为参数，表示为

Ξ = [\begin{matrix} ξ_{11} & ξ_{12} & ξ_{13} \\ ξ_{21} & ξ_{22} & ξ_{23} \\ ξ_{31} & ξ_{32} & ξ_{33} \end{matrix}]

（2）

参数矩阵 $Ξ = μ M$ ，其中μ代表尺度参数，M为一个正交旋转矩阵，且 $M = M_{3} M_{2} M_{1}$ 。

M_{1} = [\begin{matrix} 1 & 0 & 0 \\ 0 & c o s β_{1} & s i n β_{1} \\ 0 & - s i n β_{1} & c o s β_{1} \end{matrix}]

M_{2} = [\begin{matrix} c o s β_{2} & 0 & - s i n β_{2} \\ 0 & 1 & 0 \\ s i n β_{2} & 0 & c o s β_{2} \end{matrix}]

M_{3} = [\begin{matrix} c o s β_{3} & s i n β_{3} & 0 \\ - s i n β_{3} & c o s β_{3} & 0 \\ 0 & 0 & 1 \end{matrix}]

式中 $β_{1}$ 、 $β_{2}$ 、 $β_{3}$ 表示旋转角度。当这3个角度很小时，旋转矩阵中该角度的余弦值可近似为1，正弦值近似为0，以此来简化计算。对于大角度而言这种简化便大大损失了精度。旋转矩阵 $M$ 为正交矩阵且满足 $M M^{T} = I_{3}$ （ $I_{3}$ 代表三阶单位矩阵）^［

12‑13］，因此，

Ξ Ξ^{T} = μ^{2} I_{3}

。从而可以得到5个独立的条件方程，写成矩阵形式有^{［参考文献 14

百度学术}14］

g (ξ) = [\begin{matrix} ξ_{11}^{2} + ξ_{12}^{2} + ξ_{13}^{2} - ξ_{31}^{2} - ξ_{32}^{2} - ξ_{33}^{2} \\ ξ_{21}^{2} + ξ_{22}^{2} + ξ_{23}^{2} - ξ_{31}^{2} - ξ_{32}^{2} - ξ_{33}^{2} \\ ξ_{11} ξ_{21} + ξ_{12} ξ_{22} + ξ_{13} ξ_{23} \\ ξ_{11} ξ_{31} + ξ_{12} ξ_{32} + ξ_{13} ξ_{33} \\ ξ_{31} ξ_{21} + ξ_{32} ξ_{22} + ξ_{33} ξ_{23} \end{matrix}] = 0

（3）

将式（1）记为 $L_{i} - e_{L_{i}} = A_{i} ξ$ ，其中

L_{i} = [\begin{matrix} x_{2 i} \\ y_{2 i} \\ z_{2 i} \end{matrix}]

，

e_{L_{i}} = [\begin{matrix} e_{x 2 i} \\ e_{y 2 i} \\ e_{z 2 i} \end{matrix}]

，

A_{i} = [\begin{matrix} 1 & 0 & 0 & x_{1 i} & y_{1 i} & z_{1 i} & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 & 0 & 0 & x_{1 i} & y_{1 i} & z_{1 i} & 0 & 0 & 0 \\ 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & x_{1 i} & y_{1 i} & z_{1 i} \end{matrix}]

，

$ξ = {[\begin{matrix} Δ x & Δ y & Δ z & {(v e c (Ξ^{T}))}^{T} \end{matrix}]}^{T}$ 代表待估计的参数向量。其中vec（）代表向量化符号，即把矩阵从左至右，从上至下排列成列向量。将所有公共点的方程排列在一起可以得到

L - e_{L} = A ξ

（4）

式中： $e_{L}$ 代表观测值向量L的随机误差向量， $A$ 代表系数矩阵。

显然，上述模型是关于参数的线性函数。由于独立的参数仅有7个，而总的参数个数却有12个，参数之间无法保持函数独立，求解时需要加入约束条件。但约束条件关于参数是非线性的，因此需要将约束条件方程在参数估值 $ξ^{0}$ 处利用泰勒级数展开得到基于参数改正量的线性函数后再进行平差计算。

可以利用式（4）采用间接平差的方法先获取参数的近似值 $ξ^{0}$ ，再将式（3，4）在 $ξ^{0}$ 处泰勒展开，取一次项可得

l - e_{l} = A δ ξ

（5a）

C δ ξ + W = 0

（5b）

式中： $l = L - A ξ^{0}, W = g (ξ^{0})$ ，e_l代表l的随机误差向量， $Q_{l l} = Q_{L L} = Q$ ，l与L的权阵相同均为P。

E (e_{l}) = E (e_{L}) = 0

（6a）

D (e_{l}) = D (e_{L}) = σ_{0}^{2} Q = σ_{0}^{2} P^{- 1}

（6b）

C = [\begin{matrix} 0 & 0 & 0 & 2 ξ_{11}^{0} & 2 ξ_{12}^{0} & 2 ξ_{13}^{0} & 0 & 0 & 0 & - 2 ξ_{31}^{0} & - 2 ξ_{32}^{0} & - 2 ξ_{33}^{0} \\ 0 & 0 & 0 & 0 & 0 & 0 & 2 ξ_{21}^{0} & 2 ξ_{22}^{0} & 2 ξ_{23}^{0} & - 2 ξ_{31}^{0} & - 2 ξ_{32}^{0} & - 2 ξ_{33}^{0} \\ 0 & 0 & 0 & ξ_{21}^{0} & ξ_{22}^{0} & ξ_{23}^{0} & ξ_{11}^{0} & ξ_{12}^{0} & ξ_{13}^{0} & 0 & 0 & 0 \\ 0 & 0 & 0 & ξ_{31}^{0} & ξ_{32}^{0} & ξ_{33}^{0} & 0 & 0 & 0 & ξ_{11}^{0} & ξ_{13}^{0} & ξ_{13}^{0} \\ 0 & 0 & 0 & 0 & 0 & 0 & ξ_{31}^{0} & ξ_{32}^{0} & ξ_{33}^{0} & ξ_{21}^{0} & ξ_{22}^{0} & ξ_{23}^{0} \end{matrix}]

（7）

为求解上述模型在 $e_{l}^{T} P e_{l} = m i n$ 准则下的最优解，根据条件极值理论，构造拉格朗日函数

Φ = e_{l}^{T} P e_{l} + 2 K^{T} (C δ ξ + W)

（8）

式中K为拉格朗日乘数向量。为求 $Φ$ 的极小值，将 $Φ$ 对 $δ ξ$ 求偏导并令结果为0，则有

\frac{1}{2} \frac{\partial Φ}{\partial ξ} |(\hat{ξ}, \hat{K}) = A^{T} P l - A^{T} P A δ \hat{ξ} + C^{T} \hat{K} = 0

（9a）

\frac{1}{2} \frac{\partial Φ}{\partial K} |(\hat{ξ}, \hat{K}) = C δ \hat{ξ} + W = 0

（9b）

令 $N_{1} = A^{T} P A$ ，则式（9a）两端同时乘以 $C N_{1}^{- 1}$ 后与式（9b）相加可解得

\hat{K} = - N_{2}^{- 1} (C N_{1}^{- 1} A^{T} P l + W)

（10）

式中 $N_{2} = C N_{1}^{- 1} C^{T}$ 。将式（10）代入（9a）并整理便解得所求参数估值为

δ \hat{ξ} = (N_{1}^{- 1} - N_{1}^{- 1} C^{T} N_{2}^{- 1} C N_{1}^{- 1}) A^{T} P l - N_{1}^{- 1} C^{T} N_{2}^{- 1} W

（11）

在得到参数估值后残差可以由下式计算

{\hat{e}}_{l} = l - A δ \hat{ξ} = R l + h_{0}

（12）

其中

R = I - A (N_{1}^{- 1} - N_{1}^{- 1} C^{T} N_{2}^{- 1} C N_{1}^{- 1}) A^{T} P

（13）

h_{0} = A N_{1}^{- 1} C^{T} N_{2}^{- 1} W

（14）

2 敏感度分析及附有等式约束的最小二乘数据探测算法

在先验单位权方差因子 $σ_{0}^{2}$ 已知的情况下，要检验 ${\hat{σ}}_{0}^{2}$ 与 $σ_{0}^{2}$ 的一致性，其中 ${\hat{σ}}_{0}^{2} = \frac{{\hat{e}}_{L}^{T} P {\hat{e}}_{L}}{r} = \frac{{\hat{e}}_{L}^{T} P {\hat{e}}_{L}}{n - u + s}$ 代表验后单位权方差。原假设H₀： ${\hat{σ}}_{0}^{2}$ = $σ_{0}^{2}$ ，备择假设H₁： ${\hat{σ}}_{0}^{2}$ ≠ $σ_{0}^{2}$ 。构造统计量： $χ^{2} = (n - u + s) \frac{{\hat{σ}}_{0}^{2}}{σ_{0}^{2}} ~ χ^{2} (n - u + s)$ ，自由度为n $-$ u+s，其中， $n$ 为观测值个数， $u$ 为参数个数， $s$ 为约束方程个数。若满足 $χ_{α / 2}^{2} (n - u + s) < χ^{2} < χ_{1 - α / 2}^{2} (n - u + s)$ ， $α$ 为显著水平，则接受原假设H₀，否则，接受备择假设H₁，拒绝H₀，这意味着观测值中可能含有粗差。为了进一步检测异常点的具体位置，需要对每个观测值构造局部检验统计量^［

15］。

最小二乘法的残差加权平方和Ω在最小二乘平差中起到关键作用。Ω对观测或随机模型扰动的局部敏感性可作为粗差探测重要的指标^［

16］。不过上述研究是基于标准最小二乘法的，为了构造CL的粗差检验统计量，将敏感度分析方法运用到上述三维基准转换的CLS问题。利用矩阵微分方法推导出观测值扰动的残差加权平方和的局部敏感度指标，可以依据这些敏感度指标构造检验统计量进行数据探测。

易知式（12）中R和h₀具有以下性质

R R = R R^{T} = Q^{- 1} = Q^{- 1} R R h_{0} = h_{0}

（15）

残差的加权平方和对 $l$ 进行微分可得

\frac{\partial Ω}{\partial l} = \frac{\partial {\hat{e}}_{l}^{T}}{\partial l} \cdot \frac{\partial Ω}{\partial {\hat{e}}_{l}} = 2 R^{T} Q^{- 1} {\hat{e}}_{l}

（16）

结合式（15）可知

\frac{\partial Ω}{\partial l} = 2 Q^{- 1} R (R l + h_{0}) = 2 Q^{- 1} (R l + h_{0}) = 2 Q^{- 1} {\hat{e}}_{l}

（17）

因此，Ω对第i个观测值的微分为

\frac{\partial Ω}{\partial l_{i}} = 2 c_{i}^{T} Q^{- 1} {\hat{e}}_{l}

（18）

式中 $c_{i}$ 为一个列向量，这个列向量中的第i个位置上的元素是1，其余位置上的元素为0。可将式（18）作为评价残差的加权平方和Ω对第i个观测值的扰动灵敏度的指标。

由协方差传播率及式（16）的关系可得

\begin{array}{l} v a r (c_{i}^{T} Q^{- 1} {\hat{e}}_{l}) = σ_{0}^{2} c_{i}^{T} Q^{- 1} R Q R^{T} Q^{- 1} c_{i} = \\ σ_{0}^{2} c_{i}^{T} Q^{- 1} R c_{i} \end{array}

（19）

式中随机变量服从正态分布，可以利用数理统计中的知识构造如下统计量

l_{i, N} = \frac{c_{i}^{T} Q^{- 1} {\hat{e}}_{l}}{σ_{0} \sqrt[]{c_{i}^{T} Q^{- 1} R c_{i}}} ~ N (0,1)

（20）

式中 $N (0,1)$ 代表标准正态分布。当满足 $|l_{i, N}| > N_{\frac{α}{2}} (0,1)$ 时，说明第i个观测值的扰动对Ω的影响较大，这个观测值中很有可能存在粗差或者在进行平差时被粗差“污染”。将此观测值组成的方差剔除之后再利用剩下的方程进行平差。

在不知道先验单位权方差的情况下，可以采用平差计算之后得到的验后单位权方差公式推导新的统计量，表示如下

l_{i, T a u} = \frac{l_{i, N}}{\sqrt[]{{\hat{σ}}_{0}^{2} / σ_{0}^{2}}} = \frac{c_{i}^{T} Q^{- 1} {\hat{e}}_{l}}{{\hat{σ}}_{0} \sqrt[]{c_{i}^{T} Q^{- 1} R c_{i}}} ~ τ (n - u + s)

（21）

式中 $τ (n - u + s)$ 代表自由度为n $-$ u+s的τ分布^［

18‑20］。它和t分布存在如下关系

τ_{n - u + s} = \frac{t_{n - u + s - 1} \sqrt[]{n - u + s}}{\sqrt[]{t_{n - u + s - 1}^{2} + n - u + s - 1}}

（22）

式中 $t_{n - u + s - 1}^{}$ 代表自由度为n $-$ u+s $-$ 1的t分布。

类似地，当满足 $|l_{i, T a u}| > τ_{\frac{α}{2}} (n - u + s)$ 时，就说明第i个观测值的扰动对Ω的影响较大，这个观测值中很有可能存在粗差或者在进行平差时被粗差“污染”。将此观测值组成的方差剔除之后再利用剩下的“干净”的方程进行平差计算。

检验异常值的程序通常称为数据探测^［

17］。如果数据中有很多异常值，需要循环进行数据探测程序^{［参考文献 21‑22}21‑22］。如果单位权方差已知时，有

|l_{j, N}| = \underset{1 \leq i \leq n}{m a x} ((|l_{i, N}|) > N_{\frac{α}{2}} (0,1))

（23）

或者单位权方差未知时，有

|l_{j, T a u}| = \underset{1 \leq i \leq n}{m a x} (|l_{i, T a u}|) > τ_{\frac{α}{2}} (n - u + s)

（24）

这意味着Ω对l中的第j个观测值的扰动灵敏度较高，这个观测值可以认为是含有异常值的观测值，应该将其移除后用剩下的观测值重新进行平差计算，直到所有的统计量均通过检验即完成数据探测过程。

3 实验分析

3.1　仿真实验

假设在某一空间中均匀分布25个点。它们在原始坐标系下的坐标真值已知。选取其中15个点作为公共点，剩下10个点作为检查点。两个坐标系之间的转换参数可以设置如下：∆x=100，∆y= $-$ 100，∆z=50，μ=1.2。旋转矩阵由下列旋转角度β₁=1.0 rad，β₂=1.5 rad，β₃= $-$ 0.5 rad计算得到。由此便可以利用比例参数以及旋转角度得到式（1）中的Ξ。再利用式（1）便可以计算出这25个点在目标坐标系中的坐标转换值，将其视为真值。

$M A T L A B$ 中的normrnd函数可以生成随机误差。每次实验在公共点的原始坐标与目标坐标值中均加入由该函数生成独立的服从均值为0、方差为0.01的正态分布的随机误差；另外设置3个粗差，位置随机产生，其大小的绝对值介于先验标准差的5~30倍之间。重复进行1 000次模拟实验。

假设σ₀=0.01。具体方案表示如下：

（1）不在观测值中添加粗差，由公共点坐标利用附有等式约束的最小二乘法计算十二参数估值。

（2）在观测值中添加粗差，用附有等式约束的最小二乘法以及约束最小二乘的数据探测算法分别计算十二参数的估值，在单位权中误差的先验值已知的情况下，数据探测算法中可以使用式（20）中的统计量 $l_{i, N}$ （显著性水平 $α = 0.05$ ）进行粗差检验。

通过上述两方案得到参数估值之后，对于检查点，便可以利用之前提出的坐标转换公式计算它们在目标坐标系中的值 $[{\hat{x}}_{_{2 i}} {\hat{y}}_{2 i} {\hat{z}}_{2 i}]^{T}$ 。假设 $[x_{_{2 i}} y_{2 i} z_{2 i}]^{T}$ 为这些点在目标坐标系中的真实坐标，则转换坐标在x、y和z方向上的均方根误差分别为

R M S E_{x} = \sqrt[]{\sum_{i = 1}^{10} ({\hat{x}}_{_{2 i}} - x_{_{2 i}})^{2} / 10}

（25）

R M S E_{y} = \sqrt[]{\sum_{i = 1}^{10} ({\hat{y}}_{_{2 i}} - y_{_{2 i}})^{2} / 10}

（26）

R M S E_{x} = \sqrt[]{\sum_{i = 1}^{10} ({\hat{z}}_{_{2 i}} - z_{_{2 i}})^{2} / 10}

（27）

从而可得点位均方根误差为

\begin{array}{l} R M S E_{P} = \\ \sqrt[]{(R M S E_{x})^{2} + (R M S E_{y})^{2} + (R M S E_{z})^{2}} \end{array}

（29）

计算得到点位均方根误差如图1~3所示，具体的计算结果见表1。

图1 方案1所得点位均方根误差序列

Fig.1 Sequences of RMSE_P obtained by Scheme 1

图2 方案2所得点位均方根误差对比

Fig.2 Comparison of sequences of RMSE_P obtained by Scheme 2

图3 方案2所得结果x、y、z方向误差序列

Fig.3 Sequences of RMSE_P obtained by Scheme 2 in directions x，y，z

表1 均方根误差序列的统计信息

Table 1 Statistics for RMSE sequences

RMSE	参数	数据探测算法	CLS算法
RMSE_x	平均值	0.008 8	0.030 4
	最大值	0.040 3	0.105 4
	标准差	0.006 7	0.023 1
RMSE_y	平均值	0.008 4	0.033 7
	最大值	0.039 0	0.122 0
	标准差	0.006 5	0.025 0
RMSE_z	平均值	0.007 9	0.026 5
	最大值	0.030 5	0.084 4
	标准差	0.005 9	0.018 7
RMSE_P	平均值	0.016 7	0.061 9
	最大值	0.042 3	0.123 1
	标准差	0.007 3	0.021 2

根据以上实验结果，可以发现：当观测值中不含有粗差时，附有等式约束的最小二乘算法可以得出较好的转换参数值。然而当粗差存在时该算法会受到影响，所得结果严重扭曲，导致转换坐标的点位均方根误差增大了近3倍。

数据探测算法抵抗粗差的效果很明显，均方根误差值序列的平均值、最大值以及标准差明显减少很多，点位均方根误差的平均值、最大值以及标准差分别是CLS算法的26.9%、34.4%、34.3%，可见数据探测方法可以有效解决观测含有粗差的问题，并且得到精度较高的转换参数。

本文还进行了另外两组实验，设置的粗差个数分别为4个、5个，粗差的绝对值介于3~20倍先验标准差之间，这两组实验结果均表明数据探测算法可以准确的剔除粗差，得到精度更高的转换参数。

3.2　实例数据

使用两组实际测量的数据进行实验，验证本文提出的方法的性能与精度。

第1个例子中采用文献［

23］提供的实测三维坐标值（表2），这组数据也为文献［10⁃11，14］所采用，使用本文方法的参数结算结果见表3。第2个例子采用文献［24］中的数据，使用本文方法的参数解算结果见表4。

表2 公共点三维坐标值的实测值^{[参考文献 23

百度学术}23]

Table 2 Coordinates of the common points^{[参考文献 23

百度学术}23]

点号	源坐标系			目标坐标系
点号	x₁/m	y₁/m	z₁/m	x₂/m	y₂/m	z₂/m
80601	5 234 251.250	905 003.201 1	3 518 869.674	5 233 991.482	905 003.106 4	3 519 305.459
32127	5 218 851.932	919 148.974 9	3 537 928.348	5 218 595.021	919 152.324 4	3 538 363.627
80600	5 220 818.669	772 128.361 3	3 569 828.606	5 220 565.466	772 130.563 0	3 570 253.01
32136	5 148 067.252	803 912.306 0	3 668 491.426	5 147 806.722	803 921.322 3	3 668 928.371
80598	5 081 676.230	771 786.812 2	3 765 023.787	5 081 410.788	771 799.425 6	3 765 460.689
80597	5 022 479.060	955 283.548 7	3 801 754.143	5 022 218.176	955 297.253 8	3 802 185.975

表3 CLS方法及其数据探测算法的解算结果（算例1）

Table 3 Calculation results obtained by CLS and the corresponding data‑snooping algorithm （Example 1）

参数	CLS算法		数据探测算法
参数	参数估值	标准差	参数估值	标准差
∆x/m	-293.362 9	82.233 2	-226.745 1	26.081 3
∆y/m	40.798 1	157.556 1	27.190 8	53.898 6
∆z/m	354.730 3	85.386 2	394.981 9	26.148 4
ξ₁₁	1.000 010 667	1.209×10^-5	0.999 998 914	3.705×10^-6
ξ₁₂	0.000 021 228	2.144×10^-5	0.000 023 882	7.256×10^-6
ξ₁₃	-0.000 010 763	1.380×10^-5	-0.000 013 443	4.348×10^-6
ξ₂₁	-0.000 021 228	2.144×10^-5	-0.000 023 882	7.256×10^-6
ξ₂₂	1.000 010 667	1.209×10^-5	0.999 998 913	3.705×10^-6
ξ₂₃	0.000 018 196	1.755×10^-5	0.000 028 442	5.768×10^-6
ξ₃₁	0.000 010 763	1.380×10^-5	0.000 0134 44	4.348×10^-6
ξ₃₂	-0.000 018 196	1.755×10^-5	-0.000 028 441	5.768×10^-6
ξ₃₃	1.000 010 667	1.209×10^-5	0.999 998 914	3.705×10^-6

表4 CLS方法及其数据探测算法的解算结果（算例2）

Table 4 Calculation results obtained by CLS and the corresponding data‑snooping algorithm （Example 2）

参数	CLS算法		数据探测算法
参数	参数估值	标准差	参数估值	标准差
∆x/m	-4.706 1	0.011 6	-4.684 9	0.001 0
∆y/m	1.889 9	0.011 3	1.906 8	0.001 2
∆z/m	-2.705 2×10^-5	0.020 0	0.003 3	0.001 7
ξ₁₁	-0.276 504 732	8.175 51×10^-4	-0.270 746 440	6.881 48×10^-5
ξ₁₂	-0.959 605 627	7.850 82×10^-4	-0.963 204 676	7.245 09×10^-5
ξ₁₃	0.002 161 123	1.230 09×10^-3	0.001 149 131	9.622 63×10^-5
ξ₂₁	0.959 605 703	7.853 32×10^-4	0.963 204 784	7.247 14×10^-5
ξ₂₂	-0.276 498 367	8.169 75×10^-4	-0.270 744 778	6.877 74×10^-5
ξ₂₃	0.002 836 273	2.030 76×10^-3	0.001 418 579 5	1.733 75×10^-4
ξ₃₁	-0.002 127 027	1.920 10×10^-3	-0.001 054 697	1.617 70×10^-4
ξ₃₂	0.002 861 932	1.339 85×10^-3	0.001 490 128	1.077 86×10^-4
ξ₃₃	0.998 643 972	7.822 20×10^-4	1.000 532 294	6.866 03×10^-5
$\hat{σ}$ ₀/m	0.024 1		0.001 8

本节分别采用CLS算法及其数据探测算法根据两组实际测量的数据求解基准转换的十二参数，在数据探测过程中由于先验单位权方差因子未知，无法采用正态分布统计量，应选择式（21）推导的l_i_，Tau（α=0.05）统计量。第1组数据的计算结果以及精度信息（参数的标准差）如表3所示。由表3中数据可见，当含有粗差的观测值被剔除之后，求得的转换参数尤其是3个平移参数有很大的变化。相比CLS算法，数据探测算法所得参数的标准差下降了67%左右，得到了更准确、精度更高的转换参数。对比第2组数据的计算结果，由于CLS算法无法抵抗粗差，所得的验后单位权误差较大，达到2.4 cm。而数据探测算法所得的验后单位权误差仅为1.8 mm，与实验仪器的标称精度1~2 mm符合。各参数的标准差也有显著的改善。在计算过程中，第7、8、5、6个观测方程相继被识别为含粗差的方程并剔除，均为2号点和3号点对应的观测方程，这与测量过程中2号、3号标靶球被小幅移动有关，因此数据探测算法可以提高参数估计的稳健性。

4 结论

本文将大角度三维基准转换问题抽象成附有等式约束的最小二乘问题，采用十二参数模型来求解大角度下的三维基准转换参数，推导了参数在正交矩阵条件约束下的最小二乘解。当观测值中含有粗差时所得结果会严重失真。为解决这一问题，本文将敏感度分析方法运用到上述三维基准转换的CLS问题，利用矩阵微分方法推导出观测值扰动的残差加权平方和的局部敏感度指标。在先验单位权方程已知和未知的情况下构造了两个不同的检验统计量用于数据探测。模拟实验与实测数据实验结果均表明所提出的算法可以有效地抵抗观测粗差的不利影响，参数的准确性与精度得到明显的提高，从而有效地解决了大角度三维基准转换中的粗差处理问题。

参考文献

WOLF H. Geometric connection and re-orientation of three-dimensional triangulation nets［J］.Bulletin Géodésique，1963， 37（2）： 165-169. [百度学术]

MOLODENSKY M S， EREMEEV V F， YURKINA M I. Methods for the study of the external gravitational field and figure of the earth［M］. Jerusalem： Israeli Programme for the Translation of Scientific Publications， 1962. [百度学术]

姚宜斌，黄承猛，李程春，等.一种适用于大角度的三维坐标转换参数求解算法［J］. 武汉大学学报：信息科学版， 2012， 37（3）： 253-256. [百度学术]

YAO Yibin， HUANG Chengmeng， LI Chengchun， et al. A new algorithm for solution of transformation parameters of big rotation angle’s 3D coordinate［J］. Geomatics & Information Science of Wuhan University， 2012， 37（3）： 253-256. [百度学术]

林鹏，刘超，高井祥，等.附加约束的三维基准转换的高斯-赫尔默特模型［J］.中国矿业大学学报， 2017， 46（5）： 1152-1158. [百度学术]

LIN Peng， LIU Chao， GAO Jingxiang， et al. Three-dimensional datum transformation based on Gauss-Helmet model with constraints［J］. Journal of China University of Mining & Techonology，2017，46（5）： 1152-1158. [百度学术]

SHEN Y Z， CHEN Y， ZHENG D H. A quaternion-based geodetic datum transformation algorithm［J］. Journal of Geodesy， 2006，80（5）： 233-239. [百度学术]

李博峰，黄善琪.大角度三维基准转换的解析封闭解［J］. 测绘学报，2016，45（3）： 267-273. [百度学术]

LI Bofeng， HUANG Shanqi. Analytical close-form solutions for three-dimensional datum transformation with big rotation angles［J］. Acta Geodaetica et Cartographica Sinica， 2016， 45（3）： 267-273. [百度学术]

LEHMANN R. On the formulation of the alternative hypothesis for geodetic outlier detection［J］. Journal of Geodesy， 2013， 87（4）： 373-386. [百度学术]

COOK R D. Detection of influential observation in linear regression［J］. Technometrics， 1977， 19（1）： 15-18. [百度学术]

GUO J F， OU J K， WANG H T. Robust estimation for correlated observations： Two local sensitivity-based down weighting strategies［J］. Journal of Geodesy， 2010， 84（4）： 243-250. [百度学术]

WANG B， YU J， LIU C， et al. Data snooping algorithm for universal 3D similarity transformation based on generalized EIV model［J］. Measurement， 2018， 119： 56-62. [百度学术]

WANG B， FANG X， LIU C， et al. Data snooping for the equality constrained nonlinear Gauss‑Helmert model using sensitivity analysis［J］. Journal of Surveying Engineering， 2020， 146（4）： 04020015. [百度学术]

WANG B， LI J， LIU C， et al. Generalized total least squares prediction algorithm for universal 3D similarity transformation［J］.Advances in Space Research， 2017， 59（3）： 815-823. [百度学术]

XU T， CHANG G， WANG Q， et al. Analytical 3D rotation estimation using vector measurements with full variance-covariance matrix［J］. Measurement，2017，98： 131-138. [百度学术]

FANG X. Weighted total least-squares with constraints： A universal formula for geodetic symmetrical transformations［J］. Journal of Geodesy， 2015， 89（5）： 459-469. [百度学术]

TEUNISSEN P J G. Distributional theory for the DIA method［J］. Journal of Geodesy， 2017， 92（2）： 59-80. [百度学术]

GUO J F， OU J K， WANG H T. Quasi-accurate detection of outliers for correlated observations［J］. Journal of Surveying Engineering， 2007， 133（3）： 129-133. [百度学术]

BAARDA W. A testing procedure for use in geodetic networks［J］. Publication on Geodesy， New Series， 1968， 2（5）： 1-42. [百度学术]

POPE A J. The statistics of residuals and the detection of outliers： NOAA Technical Report NOS NGS； 65-1［R］. ［S.l.］： NOAA， 1976. [百度学术]

LEHMANN R. Improved critical values for extreme normalized and studentized residuals in Gauss‑Markov models［J］. Journal of Geodesy， 2012， 86（12）：1137-1146. [百度学术]

BASELGA S. Critical limitation in use of τ test for gross error detection［J］. Journal of Surveying Engineering， 2007， 133（2）： 52-55. [百度学术]

AMIRI-SIMKOOEI A R， JAZAERI S. Data-snooping procedure applied to errors-in-variables models［J］. Studia Geophysicaet Geodaetica，2013，57（3）：426-441. [百度学术]

ROFATTO V F， MATSUOKA M T， KLEIN I， et al. A half-century of Baarda’s concept of reliability： A review， new perspectives， and applications［J］. Survey Review， 2018，52（1）： 261-277. [百度学术]

FELUS Y A， BURTCH R C. On symmetrical three-dimensional datum conversion［J］. GPS Solut，2009， 13： 65-74. [百度学术]

王彬.广义整体最小二乘的拓展理论及其在测量数据处理中的应用研究［D］.武汉：武汉大学，2017. [百度学术]

WANG Bin. Study on the extended theories of generalized total least squares and their applications in surveying data processing［D］. Wuhan， China：Wuhan University，2017. [百度学术]

请使用 Firefox、Chrome、IE10、IE11、360极速模式、搜狗极速模式、QQ极速模式等浏览器，其他浏览器不建议使用!