An Improved Region-Growing Surface Triangulation Algorithm for Urban Dense Point Cloud
-
摘要: 随着飞行平台与传感器技术的发展,空间信息的获取能力逐渐增强,特别是低空倾斜摄影及其数据处理方法的出现,使得城市真三维信息的获取手段变得更加丰富。针对现有区域生长算法在城市三维密集点云拐角处未能完全生长的问题,设计了点到边(point to edge,PTE)数据结构,提出了遍历PTE数据结构以完善构网的改进策略。实验结果表明,改进算法能对城市真三维密集点云进行较好的表面构网,并且能够解决现有方法在拐角处不能完全生长的问题。Abstract: With the development of flying UAV platforms, sensor technologies, and corresponding data processing methods, low-altitude oblique photography becomes a powerful method to obtain true three-dimensional urban spatial information. In existing region-growing methods however, corners within urban dense point clouds are not fully grown. In this paper, we propose an improvement over the previous mesh-growing algorithms for fast surface reconstruction, a PTE (point to edge) data structure and traversing algorithm. Experimental results show that the improved algorithm is able to handle dense true three-dimensional urban point clouds, and can solve most of the incomplete growing problems in corners, with better performance than the exiting methods.
-
经典平差模型和最小二乘估计理论[1]在大地测量等众多科学研究和工程领域中应用广泛，其中，高斯-马尔科夫模型（Gauss-Markov model, GMM）最为常用，而高斯-赫尔默特模型（GaussHelmert model, GHM）可视为经典平差模型的一般通用形式。在实际应用中，坐标转换、回归模型、数字地面模型和大地测量反演等平差模型的系数矩阵包含随机的观测误差，从而使得GMM扩展为随机系数矩阵的变量含误差（errors-invariables, EIV）模型[2]。文献[3]提出同时顾及观测向量和系数矩阵中随机误差的整体最小二乘（total least squares, TLS）估计算法。TLS的非线性特征导致其受制于计算机技术的发展，直至20世纪80年代，文献[4]将TLS引入数值分析领域并提出奇异值分解算法，TLS才开始广泛应用于各专业领域并取得丰富的研究成果。文献[2]中对TLS进行了改进和扩展；文献[5]从TLS的算法、统计特性和可靠性研究等方面综述了TLS方法的研究进展。当误差相关且精度不等时，采用加权整体最小二乘估计（weighted total least squares, WTLS）方法进行求解。文献[6]研究了基于高斯-牛顿迭代法的WTLS算法，该算法假设权矩阵为特殊情况，得到的解在形式上与最小二乘（least squares, LS）解相同；文献[7]研究了在任意权矩阵的一般情况下的WTLS算法；文献[8]研究了特殊结构下WTLS算法的迭代方法并将其应用于实际场景；文献[9-12]研究了附有等式和不等式约束情况下的WTLS算法；文献[13]研究了稳健WTLS算法。
通过对EIV模型的形式进行变换，文献[14]提出部分EIV(partial EIV, PEIV）模型，提高了系数矩阵仅含部分随机量情况下的计算效率。文献[15]对PEIV模型进行线性化，推导了PEIV模型的LS算法。文献[16-17]从模型的一般性出发，将EIV模型扩展至通用EIV模型，将经典平差的GHM中观测向量的系数矩阵和参数向量的系数矩阵由固定矩阵推广为随机矩阵，涵括随机系数矩阵的各类情况，同时推导了通用EIV模型在任意权矩阵情况下的一般性WTLS算法。
通用EIV模型的非线性使得该算法在估计量较多时计算量大。本文利用非线性平差原理，将通用EIV模型展开后的二阶项纳入平差方程的常数项，从而将其转化为GHM形式，推导出通用EIV模型的线性化整体最小二乘（linearized total least squares, LTLS）算法。相较于WTLS算法，LTLS算法提高了通用EIV模型的计算效率，当参数向量初始值与最优值相差较大时，提升了迭代收敛速度。
1 通用EIV平差模型及其WTLS估计
1.1 通用EIV平差模型
GHM的形式为：
$$ \mathit{\boldsymbol{A}}\left( {\mathit{\boldsymbol{y}} + {\mathit{\boldsymbol{v}}_\mathit{\boldsymbol{y}}}} \right) + \mathit{\boldsymbol{BX}} + \mathit{\boldsymbol{w}} = \mathit{\boldsymbol{0}} $$ (1) 式中，y和vy分别为n×1维观测值向量和观测值改正数向量；X为u×1维参数向量；A为观测值向量对应的f×n维系数矩阵；B为参数向量对应的f×u维系数矩阵；w为f×1维常数向量；在经典平差函数模型的定义中，A和B均不含随机误差，为固定矩阵。
当参数向量的系数矩阵B含随机误差时，GHM（式（1））扩展为经典EIV模型。当观测值向量的系数矩阵A和参数向量的系数矩阵B均含随机误差时，GHM（式（1））扩展为通用EIV平差模型[16]:
$$ \left( {\mathit{\boldsymbol{A}} + {\mathit{\boldsymbol{V}}_\mathit{\boldsymbol{A}}}} \right)\left( {\mathit{\boldsymbol{y}} + {\mathit{\boldsymbol{v}}_\mathit{\boldsymbol{y}}}} \right) + \left( {\mathit{\boldsymbol{B}} + {\mathit{\boldsymbol{V}}_\mathit{\boldsymbol{B}}}} \right)\mathit{\boldsymbol{X}} + \mathit{\boldsymbol{w}} = \mathit{\boldsymbol{0}} $$ (2) 式中，A和VA分别为观测值向量对应的f×n维系数矩阵及其改正数矩阵；B和VB分别为参数向量对应的f×u维系数矩阵及其改正数矩阵。由于A、B和y均为随机矩阵，则通用EIV的随机模型为：
$$ \mathit{\boldsymbol{L}} = \left[ {\begin{array}{*{20}{c}} {{\mathop{\rm vec}\nolimits} (\mathit{\boldsymbol{A}})}\\ {{\mathop{\rm vec}\nolimits} (\mathit{\boldsymbol{B}})}\\ \mathit{\boldsymbol{y}} \end{array}} \right],\mathit{\boldsymbol{v}} = \left[ {\begin{array}{*{20}{c}} {{\mathop{\rm vec}\nolimits} \left( {{\mathit{\boldsymbol{V}}_\mathit{\boldsymbol{A}}}} \right)}\\ {{\mathop{\rm vec}\nolimits} \left( {{\mathit{\boldsymbol{V}}_\mathit{\boldsymbol{B}}}} \right)}\\ {{\mathit{\boldsymbol{v}}_\mathit{\boldsymbol{y}}}} \end{array}} \right] = \left[ {\begin{array}{*{20}{l}} {{\mathit{\boldsymbol{v}}_\mathit{\boldsymbol{A}}}}\\ {{\mathit{\boldsymbol{v}}_\mathit{\boldsymbol{B}}}}\\ {{\mathit{\boldsymbol{v}}_\mathit{\boldsymbol{y}}}} \end{array}} \right] $$ (3) $$ \mathit{\boldsymbol{D}}(\mathit{\boldsymbol{L}}) = \delta _0^2{\mathit{\boldsymbol{P}}^{ - 1}} = \delta _0^2\mathit{\boldsymbol{Q}} = \delta _0^2\left[ {\begin{array}{*{20}{c}} {{\mathit{\boldsymbol{Q}}_\mathit{\boldsymbol{A}}}}&{{\mathit{\boldsymbol{Q}}_{\mathit{\boldsymbol{AB}}}}}&{{\mathit{\boldsymbol{Q}}_{\mathit{\boldsymbol{Ay}}}}}\\ {{\mathit{\boldsymbol{Q}}_{\mathit{\boldsymbol{BA}}}}}&{{\mathit{\boldsymbol{Q}}_\mathit{\boldsymbol{B}}}}&{{\mathit{\boldsymbol{Q}}_{\mathit{\boldsymbol{By}}}}}\\ {{\mathit{\boldsymbol{Q}}_{\mathit{\boldsymbol{yA}}}}}&{{\mathit{\boldsymbol{Q}}_{\mathit{\boldsymbol{yB}}}}}&{{\mathit{\boldsymbol{Q}}_\mathit{\boldsymbol{y}}}} \end{array}} \right] $$ (4) 式中，vec(.)表示将矩阵按列向量化；L和v分别为观测数据的k×1维观测值向量及其改正数向量，包括A、B和y中所有观测值及其改正数，其中k=fn+fu+n; P、Q和D(L)分别为L的权矩阵、协因数矩阵和方差协方差矩阵；δ02为单位权方差。
1.2 WTLS估计
根据TLS准则，通用EIV平差模型的求解可转化为最优化估计问题[16]:
$$ \left\{ {\begin{array}{*{20}{l}} {\min {\mathit{\boldsymbol{v}}^{\rm{T}}}\mathit{\boldsymbol{Pv}}}\\ {s.t.{\rm{ }}\left( {\mathit{\boldsymbol{A}} + {\mathit{\boldsymbol{V}}_\mathit{\boldsymbol{A}}}} \right)\left( {\mathit{\boldsymbol{y}} + {\mathit{\boldsymbol{v}}_\mathit{\boldsymbol{y}}}} \right) + \left( {\mathit{\boldsymbol{B}} + {\mathit{\boldsymbol{V}}_\mathit{\boldsymbol{B}}}} \right)\mathit{\boldsymbol{X}} + \mathit{\boldsymbol{w}} = \mathit{\boldsymbol{0}}} \end{array}} \right. $$ (5) 相应目标函数为：
$$ \begin{array}{*{20}{c}} {\mathit{\Phi }(\mathit{\boldsymbol{r}},\mathit{\boldsymbol{\lambda }},\mathit{\boldsymbol{X}}) = {\mathit{\boldsymbol{v}}^{\rm{T}}}\mathit{\boldsymbol{Pv}} + 2{\mathit{\boldsymbol{\lambda }}^{\rm{T}}}\left( {\mathit{\boldsymbol{Ay}} + \mathit{\boldsymbol{A}}{\mathit{\boldsymbol{v}}_\mathit{\boldsymbol{y}}} + } \right.}\\ {\left. {{\mathit{\boldsymbol{V}}_\mathit{\boldsymbol{A}}}\mathit{\boldsymbol{y}} + {\mathit{\boldsymbol{V}}_\mathit{\boldsymbol{A}}}{\mathit{\boldsymbol{v}}_\mathit{\boldsymbol{y}}} + \mathit{\boldsymbol{BX}} + {\mathit{\boldsymbol{V}}_\mathit{\boldsymbol{B}}}\mathit{\boldsymbol{X}} + \mathit{\boldsymbol{w}}} \right)} \end{array} $$ (6) 将目标函数对估计量分别求偏导并令其等于0，得到非线性方程组：
$$ \frac{{\partial \mathit{\Phi }}}{{\partial \mathit{\boldsymbol{\hat X}}}} = 2\left( {{\mathit{\boldsymbol{B}}^{\rm{T}}}\mathit{\boldsymbol{\hat \lambda }} + \mathit{\boldsymbol{\hat V}}_\mathit{\boldsymbol{B}}^{\rm{T}}\mathit{\boldsymbol{\hat \lambda }}} \right) = \mathit{\boldsymbol{0}} $$ (7) $$ \frac{{\partial \mathit{\Phi }}}{{\partial \mathit{\boldsymbol{\hat v}}}} = 2\left( {\mathit{\boldsymbol{P\hat v}} + {{\mathit{\boldsymbol{\hat C}}}^{\rm{T}}}\mathit{\boldsymbol{\hat \lambda }}} \right) = \mathit{\boldsymbol{0}} $$ (8) $$ \frac{{\partial \mathit{\Phi }}}{{\partial \mathit{\boldsymbol{\hat \lambda }}}} = 2(\mathit{\boldsymbol{Ay}} + \mathit{\boldsymbol{B\hat X}} + \mathit{\boldsymbol{\hat C\hat v}} + \mathit{\boldsymbol{w}}) = \mathit{\boldsymbol{0}} $$ (9) 式中，$ {\boldsymbol{C}}=\left[{\boldsymbol{y}}^{\mathrm{T}} \otimes {\boldsymbol{I}}_{f} \boldsymbol{X}^{\mathrm{T}} \otimes {\boldsymbol{I}}_{f} A+\boldsymbol{V}_{{\boldsymbol{A}}}\right] ; \hat{{\boldsymbol{v}}}, \hat{\boldsymbol{X}} $分别为观测值向量和参数向量的估计值。
根据式（7）~式（9）可导出：
$$ \mathit{\boldsymbol{\hat v}} = - \mathit{\boldsymbol{Q}}{\mathit{\boldsymbol{\hat C}}^{\rm{T}}}\mathit{\boldsymbol{\hat Q}}_\mathit{\boldsymbol{C}}^{ - 1}(\mathit{\boldsymbol{Ay}} + \mathit{\boldsymbol{B\hat X}} + \mathit{\boldsymbol{w}}) $$ (10) $$ \begin{array}{*{20}{c}} {\mathit{\boldsymbol{\hat X}} = - {{\left[ {\left( {{\mathit{\boldsymbol{B}}^{\rm{T}}} + \mathit{\boldsymbol{\hat V}}_\mathit{\boldsymbol{B}}^{\rm{T}}} \right)\mathit{\boldsymbol{\hat Q}}_\mathit{\boldsymbol{C}}^{ - 1}\mathit{\boldsymbol{B}}} \right]}^{ - 1}}\left( {{\mathit{\boldsymbol{B}}^{\rm{T}}} + } \right.}\\ {\left. {\mathit{\boldsymbol{\hat V}}_\mathit{\boldsymbol{B}}^{\rm{T}}} \right)\mathit{\boldsymbol{\hat Q}}_\mathit{\boldsymbol{C}}^{ - 1}(\mathit{\boldsymbol{Ay}} + \mathit{\boldsymbol{w}})} \end{array} $$ (11) 式中，$ \hat{\boldsymbol{Q}}_{C}=\hat{\boldsymbol{C}} \boldsymbol{Q} \hat{\boldsymbol{C}}^{\mathrm{T}} $。
以式（2）的LS解作为初始值，根据式（10）和式（11）进行迭代计算可得通用EIV模型的WTLS最优解。
2 通用EIV模型的线性化估计算法
通用EIV模型是非线性模型，观测值矩阵和系数矩阵均为随机量，WTLS算法的计算量随着待估量个数增多将迅速增加。将式（2）展开，利用非线性函数平差原理[18]将二阶项作为模型误差纳入方程的常数项，从而将通用EIV模型转化为线性的GHM，推导出通用EIV模型的LTLS算法。
令X=X0+x，将式（2）展开得：
$$ \begin{array}{c} \left( {{\mathit{\boldsymbol{V}}_\mathit{\boldsymbol{A}}}\mathit{\boldsymbol{y}} + {\mathit{\boldsymbol{V}}_\mathit{\boldsymbol{B}}}{\mathit{\boldsymbol{X}}_0} + \mathit{\boldsymbol{A}}{\mathit{\boldsymbol{v}}_\mathit{\boldsymbol{y}}}} \right) + \mathit{\boldsymbol{Bx}} + (\mathit{\boldsymbol{w}} + \mathit{\boldsymbol{Ay}} + \\ \left. {\mathit{\boldsymbol{B}}{\mathit{\boldsymbol{X}}_0} + {\mathit{\boldsymbol{V}}_\mathit{\boldsymbol{A}}}{\mathit{\boldsymbol{v}}_\mathit{\boldsymbol{y}}} + {\mathit{\boldsymbol{V}}_\mathit{\boldsymbol{B}}}\mathit{\boldsymbol{x}}} \right) = \mathit{\boldsymbol{0}} \end{array} $$ (12) 式（12）可表示为：
$$ \begin{array}{*{20}{c}} {\left[ {\left( {{\mathit{\boldsymbol{y}}^{\rm{T}}} \otimes {\mathit{\boldsymbol{I}}_\mathit{\boldsymbol{f}}}} \right){\mathop{\rm vec}\nolimits} \left( {{\mathit{\boldsymbol{V}}_\mathit{\boldsymbol{A}}}} \right) + \left( {\mathit{\boldsymbol{X}}_0^{\rm{T}} \otimes {\mathit{\boldsymbol{I}}_\mathit{\boldsymbol{f}}}} \right){\mathop{\rm vec}\nolimits} \left( {{\mathit{\boldsymbol{V}}_\mathit{\boldsymbol{B}}}} \right) + \mathit{\boldsymbol{A}}{\mathit{\boldsymbol{v}}_\mathit{\boldsymbol{y}}}} \right] + }\\ {\mathit{\boldsymbol{Bx}} + \left( {\mathit{\boldsymbol{w}} + \mathit{\boldsymbol{Ay}} + \mathit{\boldsymbol{B}}{\mathit{\boldsymbol{X}}_0} + {\mathit{\boldsymbol{V}}_\mathit{\boldsymbol{A}}}{\mathit{\boldsymbol{v}}_\mathit{\boldsymbol{y}}} + {\mathit{\boldsymbol{V}}_\mathit{\boldsymbol{B}}}\mathit{\boldsymbol{x}}} \right) = \mathit{\boldsymbol{0}}} \end{array} $$ (13) 则通用EIV模型的线性化形式为：
$$ {\mathit{\boldsymbol{A}}_\mathit{\boldsymbol{l}}}{\mathit{\boldsymbol{v}}_\mathit{\boldsymbol{l}}} + \mathit{\boldsymbol{Bx}} + {\mathit{\boldsymbol{w}}_\mathit{\boldsymbol{l}}} = \mathit{\boldsymbol{0}} $$ (14) 式中，$ {\boldsymbol{A}}_{{\boldsymbol{l}}}=\left[\begin{array}{llll} {\boldsymbol{y}}^{\mathrm{T}} \otimes {\boldsymbol{I_{f}}} & \boldsymbol{X}_{0}^{\mathrm{T}} & \otimes {\boldsymbol{I_{f}}} & {\boldsymbol{A}} \end{array}\right] ; \ {\boldsymbol{v_{l}}}= \left[\begin{array}{c} \operatorname{vec}\left(\boldsymbol{V}_{{\boldsymbol{A}}}\right) \\ \operatorname{vec}\left(\boldsymbol{V}_{{\boldsymbol{B}}}\right) \\ \boldsymbol{v}_{{\boldsymbol{y}}} \end{array}\right] ; \boldsymbol{w}_{{\boldsymbol{l}}}=\boldsymbol{w}+\boldsymbol{Ay} + {\boldsymbol{B}} \boldsymbol{X}_{0}+\boldsymbol{V}_{{\boldsymbol{A}}} \boldsymbol{v}_{{\boldsymbol{y}}}+\boldsymbol{V}_{{\boldsymbol{B}}} {\boldsymbol{x}}$为常数项，即将式（2）按泰勒展开略去的二阶项VAvy+VBx作为模型误差纳入常数项。
由式（14）可知，线性化后的通用EIV模型与GHM形式一致，可使用最小二乘法得到观测值向量和参数向量：
$$ \mathit{\boldsymbol{\hat x}} = - \mathit{\boldsymbol{\hat N}}_{\mathit{\boldsymbol{bb}}}^{ - 1}{\mathit{\boldsymbol{B}}^{\rm{T}}}\mathit{\boldsymbol{\hat N}}_{\mathit{\boldsymbol{aa}}}^{ - 1}{\mathit{\boldsymbol{w}}_\mathit{\boldsymbol{l}}} $$ (15) $$ {{\mathit{\boldsymbol{\hat v}}}_\mathit{\boldsymbol{l}}} = - \mathit{\boldsymbol{Q}}{{\mathit{\boldsymbol{\hat A}}}_\mathit{\boldsymbol{l}}}\mathit{\boldsymbol{\hat N}}_{\mathit{\boldsymbol{aa}}}^{ - 1}\left( {\mathit{\boldsymbol{B\hat x}} + {\mathit{\boldsymbol{w}}_\mathit{\boldsymbol{l}}}} \right) $$ (16) 式中，$ \boldsymbol{N}_{{\boldsymbol{a a}}}=\boldsymbol{A}_{{\boldsymbol{l}}} \boldsymbol{Q} \boldsymbol{A}_{{\boldsymbol{l}}}^{\mathrm{T}} ; \boldsymbol{N}_{{\boldsymbol{b b}}}= {\boldsymbol{B}}^{\mathrm{T}} \boldsymbol{N}_{{\boldsymbol{a a}}}^{-1} {\boldsymbol{B}} $; Q为观测值向量的协因数矩阵。
通用EIV模型的LTLS算法步骤如下：
1）将实际模型表示为式（2），将观测数据代入得到A、B和y矩阵，并给出观测值数据的协因数阵Q，包括观测值向量Qy、观测值向量系数矩阵QA、参数向量系数矩阵QB。
2）计算通用EIV模型的LS解作为初始参数解：$ \hat{\boldsymbol{X}}^{0}=-\left({\boldsymbol{B}}^{\mathrm{T}}\left({\boldsymbol{A }}\boldsymbol{Q}_{{\boldsymbol{y}}} {\boldsymbol{A}}^{\mathrm{T}}\right)^{-1} {\boldsymbol{B}}\right)^{-1} {\boldsymbol{B}}^{\mathrm{T}}\left({\boldsymbol{A}} \boldsymbol{Q}_{{\boldsymbol{y}}} {\boldsymbol{A}}^{\mathrm{T}}\right)^{-1} · ({\boldsymbol{A y}}+{\boldsymbol{w}})$, 观测值向量改正数初始值取0。
3）根据式（15）和式（16）进行迭代计算，每次迭代将上一次估计值作为初始值代入新的迭代过程，直至前后两次估计值之差小于设定阈值。
GHM按泰勒级数展开仅包含常数项、一阶项（二阶及以上项全部为零），将二阶项纳入线性化后的常数项，该方法极大地减弱了线性化引起的模型误差。因此，在同样以LS解作为初值的情况下，GHM线性化的LTLS解与WTLS解的收敛性一致。此外，根据文献[19]中EIV模型LS解偏差的研究结果，在当前测量技术手段和观测精度条件下，以有偏的LS解作为初值，能够保证TLS迭代计算收敛，除非出现极特殊情况导致LS初始解严重偏离最优值。
参考GHM，式（2）的LTLS算法估计结果的精度计算式为：
$$ \left\{ {\begin{array}{*{20}{l}} {\hat \delta _0^2 = \mathit{\boldsymbol{v}}_\mathit{\boldsymbol{l}}^{\rm{T}}\mathit{\boldsymbol{P}}{\mathit{\boldsymbol{v}}_\mathit{\boldsymbol{l}}}/\mathit{r}}\\ {\mathit{\boldsymbol{Q}}(\mathit{\boldsymbol{\hat x}}) = {{\left( {{\mathit{\boldsymbol{B}}^{\rm{T}}}{{\left( {{\mathit{\boldsymbol{A}}_\mathit{\boldsymbol{l}}}\mathit{\boldsymbol{QA}}_\mathit{\boldsymbol{l}}^{\rm{T}}} \right)}^{ - 1}}\mathit{\boldsymbol{B}}} \right)}^{ - 1}}}\\ {\mathit{\boldsymbol{Q}}(\mathit{\boldsymbol{\hat L}}) = \mathit{\boldsymbol{Q}} - \mathit{\boldsymbol{QA}}_\mathit{\boldsymbol{l}}^{\rm{T}}\left( {\mathit{\boldsymbol{N}}_\mathit{\boldsymbol{A}}^{ - 1}\left( {\mathit{\boldsymbol{I}} - \mathit{\boldsymbol{B}}{\mathit{\boldsymbol{Q}}_{\mathit{\boldsymbol{\hat x}}}}{\mathit{\boldsymbol{B}}^{\rm{T}}}\mathit{\boldsymbol{N}}_\mathit{\boldsymbol{A}}^{ - 1}} \right){\mathit{\boldsymbol{A}}_\mathit{\boldsymbol{l}}}\mathit{\boldsymbol{Q}}} \right)}\\ {\mathit{\boldsymbol{D}}(\mathit{\boldsymbol{\hat x}}) = \hat \delta _0^2\mathit{\boldsymbol{Q}}(\mathit{\boldsymbol{\hat x}})}\\ {\mathit{\boldsymbol{D}}(\mathit{\boldsymbol{\hat L}}) = \hat \delta _0^2\mathit{\boldsymbol{Q}}(\mathit{\boldsymbol{\hat L}})} \end{array}} \right. $$ (17) 式中，多余观测数r=n-t; n为观测值个数；t为必要观测值个数；r与GHM的多余观测数相同；$ \boldsymbol{N}_{{\boldsymbol{A}}}=\boldsymbol{A}_{{\boldsymbol{l}}} \boldsymbol{P}^{-1} \boldsymbol{A}_{{\boldsymbol{l}}}^{\mathrm{T}} $。
3 LTLS算法实例分析
按照LTLS算法步骤设计实验，比较LTLS算法与WTLS算法的计算结果，验证LTLS算法的正确性、高效性和可行性。实验1设计模拟数据，比较分析单组实验结果和1 000组实验统计结果，验证LTLS算法的正确性；实验2在待估计量数目取不同量级时，比较两种算法的计算时间，验证LTLS算法的高效性；实验3通过实例验证LTLS算法的可行性。
3.1 实验1
在通用EIV模型（式（2））中，设置参数真值X=[5 10]T，系数矩阵A和B中随机量的中误差分别设为0.01和0.02，观测向量y的中误差设为0.03。A、B、y和常数向量w的模拟数据如下：
$$ \mathit{\boldsymbol{A}} = \left[ {\begin{array}{*{20}{c}} {12.469}&{11.096}&{15.872}&{11.725}\\ {8.883}&{10.291}&{2.929}&{3.666}\\ {12.321}&{1.109}&{6.392}&{15.809}\\ {3.551}&{12.104}&{3.867}&{1.257} \end{array}} \right], $$ $$ \mathit{\boldsymbol{B}} = \left[ {\begin{array}{*{20}{l}} {10.410}&{17.544}\\ {18.033}&{15.171}\\ {18.631}&{15.855}\\ {15.671}&{11.878} \end{array}} \right], $$ $$ \mathit{\boldsymbol{y}} = \left[ {\begin{array}{*{20}{l}} {27.543}\\ {20.727}\\ {20.839}\\ {25.033} \end{array}} \right],\mathit{\boldsymbol{w}} = \left[ {\begin{array}{*{20}{c}} { - 1425.323}\\ { - 852.619}\\ { - 1142.913}\\ { - 658.407} \end{array}} \right]。 $$ 采用LTLS算法和WTLS算法的估计结果如表 1所示，LTLS参数解与WTLS参数解完全相等，验证了LTLS算法的正确性。
表 1 参数解及其方差估计值Table 1. Parameter Values and Mean Square Deviations算法 参数解 中误差 $ \hat {\boldsymbol{X}} _1 $ $ \hat {\boldsymbol{X}} _2 $ $ σ _{\hat{\boldsymbol{X}}_1} $ $ σ_ {\hat{\boldsymbol{X}}_2} $ LTLS算法 5.012 551 9.994 964 0.039 9 0.051 9 WTLS算法 5.012 551 9.994 964 0.039 9 0.051 9 为了进一步验证LTLS算法的正确性，首先采用模拟的1 000组数据计算LTLS和WTLS参数解的均值$ {\rm{avg}} (\hat {\boldsymbol{X}} _1) 和 {\rm{avg}} (\hat {\boldsymbol{X}} _2) $，并将参数解均值代入式（17），求得LTLS参数解的协因数阵估值$ {\boldsymbol{Q}}(\hat {\boldsymbol{X}} )$，然后利用参数真值求得精确的参数协因数阵$ \overline{\boldsymbol{Q}}(\hat{\boldsymbol{X}})={\boldsymbol{E}}_{{\boldsymbol{X}}}^{\mathrm{T}} \boldsymbol{E}_{{\boldsymbol{X}}} /(m-n) , {\boldsymbol{E}}_{{\boldsymbol{X}}}=\left[\begin{array}{ll} \hat{{\boldsymbol{X}}}_{1}^{(j)}-5 & \hat{{\boldsymbol{X}}}_{2}^{(j)}-10 \end{array}\right] $，计算结果见表 2。1 000组数据的LTLS参数解和WTLS参数解与参数真值偏差的统计分析见图 1和图 2。
表 2 1000组实验的参数解均值和协因数阵Table 2. Average Parameter Values and Co⁃variance Matrix in 1 000 Experiments算法 参数解均值 协因数阵 $ {\rm{avg}} ( \hat{\boldsymbol{X}}_1) $ $ {\rm{avg}} ( \hat{\boldsymbol{X}}_2) $ $ {\boldsymbol{Q }}( \hat{\boldsymbol{X}} ) $ $ \bar{ \boldsymbol{Q }}( \hat{\boldsymbol{X}} ) $ LTLS算法 5.009 427 9.998 257 0.004 3 −0.005 1 0.004 2 −0.004 9 −0.005 1 0.007 3 −0.004 9 0.007 2 WTLS算法 5.009 427 9.998 257 0.004 3 −0.005 1 0.004 2 −0.004 9 −0.005 1 0.007 3 −0.004 9 0.007 2 ![]()
图 1 1 000组实验数据的LTLS参数解与参数真值偏差统计图Figure 1. Statistical Graph of Deviations Between Truth Values and Parameter Values Soluted by LTLS in 1 000 Experiments![]()
图 2 1 000组实验数据的WTLS参数解与参数真值偏差统计图Figure 2. Statistical Graph of Deviations Between Truth Values and Parameter Values Soluted by WTLS in 1 000 Experiments由表 2、图 1和图 2可以看出，LTLS算法与WTLS算法的统计结果完全一致，说明在每次实验中两种算法所求参数解均一致，验证了LTLS算法的正确性。同时两种计算协因数阵的方法结果非常相近，验证了参数协因数阵一阶近似估计公式（17）的有效性[17]。
3.2 实验2
为分析LTLS算法的计算效率，设计通用EIV模型中待估计量个数在不同的数量级情况，采用LTLS算法和WTLS算法计算100组模拟数据的平均迭代次数N、平均解算时间t和减少比例（LTLS算法较WTLS算法减少的平均解算时间与WTLS算法平均解算时间之比），结果见表 3。
表 3 LTLS算法和WTLS算法计算效率的比较Table 3. Comparison of Computational Efficiency Between LTLS Algorithm and WTLS Algorithm待估量数量 NLTLS NWTLS tLTLS tWTLS 减少比例/% 10 5.18 4.5 0.224 ms 0.191 ms − 100 5.04 6.36 0.681 ms 0.724 ms 5.9 1 000 5 6.75 0.033 s 0.044 s 25.0 10 000 5 7.5 2.681 s 3.909 s 31.4 从表 3可以看出，两种算法每次迭代的平均时间基本一致；当模型估计量数量较少时，两种算法的效率基本相当，随着估计量的数量级逐渐增大，LTLS算法的效率高于WTLS算法。原因在于GHM按泰勒级数展开后，仅包含常数项、一阶项（二阶及以上项全部为零），LTLS算法将二阶项纳入常数项，减小了线性化引起的模型误差，迭代计算收敛更快，迭代次数减少，使得计算效率提高。
3.3 实验3
本文采用的摄影测量实例示意图如图 3所示，由3个地面摄像机S1、S2和S3拍摄两个目标点P1和P2组成，相机主距f=100 mm，距离l1、l2、l3、l4、l5、l6、y1、y2的观测值和中误差见表 4。
表 4 距离观测值及其中误差Table 4. Distance Observations and Standard Deviations统计项 l1/mm l2/mm l3/mm l4/mm l5/mm l6/mm y1/m y2/m 观测值 14.1 16.6 6.1 7.1 22.1 26.3 10.0 8.0 中误差 0.10 0.10 0.10 0.10 0.10 0.10 0.05 0.05 根据图 3可得到误差方程：
$$ \left\{\begin{array}{l} l_{1} x_{2}-f x_{1}=0 \\ l_{2} x_{4}-f x_{3}=0 \\ l_{3} x_{2}+f y_{1}+f x_{1}=0 \\ l_{4} x_{4}+f y_{1}+f x_{3}=0 \\ l_{5} x_{2}-f y_{1}-f y_{2}+f x_{1}=0 \\ l_{6} x_{4}-f y_{1}-f y_{2}+f x_{3}=0 \end{array}\right. $$ (18) 构建通用EIV模型来估计点P1和点P2的坐标，由误差方程可得：
$$ \begin{array}{l} \mathit{\boldsymbol{A}} = \left[ {\begin{array}{*{20}{c}} 0&0\\ 0&0\\ { - f}&0\\ { - f}&0\\ { - f}&{ - f}\\ { - f}&{ - f} \end{array}} \right],\mathit{\boldsymbol{B}} = \left[ {\begin{array}{*{20}{c}} { - f}&{{l_1}}&0&0\\ 0&0&{ - f}&{{l_2}}\\ f&{{l_3}}&0&0\\ 0&0&f&{{l_4}}\\ f&{{l_5}}&0&0\\ 0&0&f&{{l_6}} \end{array}} \right],\\ \mathit{\boldsymbol{y}} = \left[ {\begin{array}{*{20}{l}} {{y_1}}\\ {{y_2}} \end{array}} \right],\mathit{\boldsymbol{x}} = \left[ {\begin{array}{*{20}{l}} {{x_1}}\\ {{x_2}}\\ {{x_3}}\\ {{x_4}} \end{array}} \right],\mathit{\boldsymbol{w}} = \left[ {\begin{array}{*{20}{l}} 0\\ 0\\ 0\\ 0\\ 0\\ 0 \end{array}} \right] \end{array} $$ (19) 由式（19）可知矩阵A不含有随机误差，则该矩阵的改正数矩阵为零矩阵。采用LTLS算法和WTLS算法解算，设置阈值为1×10-8，结果见表 5和表 6。两种方法所得坐标估值和距离观测值估值完全一致，该实例表示为通用EIV模型时系数矩阵中待估计量较少，所以两个算法收敛速度相差不大。
表 5 点P1和点P2的坐标估值/mTable 5. Coordinate Estimates of P1 and P2 /m算法 坐标估值 $ \hat x_1 $ $ \hat x_2 $ $ \hat x_3 $ $ \hat x_4 $ LTLS算法 6.995 056 5 49.715 632 6.981 465 5 41.968 315 9 WTLS算法 6.995 056 5 49.715 632 6.981 465 5 41.968 315 9 表 6 距离观测值估值Table 6. Estimation of Distance Observations算法 $ \hat l_1 $/mm $ \hat l_2 $/mm $ \hat l_3 $/mm $ \hat l_4 $/mm $ \hat l_5 $5/mm $ \hat l_6 $/mm $ \hat y_1 $/m $ \hat y_2 $/m LTLS算法 14.070 1 16.635 1 6.032 4 7.178 4 22.137 7 26.256 7 9.994 1 8.006 8 WTLS算法 14.070 1 16.635 1 6.032 4 7.178 4 22.137 7 26.256 7 9.994 1 8.006 8 4 结语
本文将通用EIV函数模型展开后的二阶项纳入模型的常数项，将通用EIV模型表示为线性形式的GHM，推导出通用EIV模型的线性化整体最小二乘算法和近似精度估计公式。实验结果表明，通用EIV模型的LTLS算法与WTLS算法结果一致，验证了该算法的正确性。此外，LTLS算法估计精度公式和WTLS估计精度公式均为一阶近似精度，因此两种算法参数的估计精度相同。当通用EIV模型的待估量数量较多时，LTLS算法比WTLS计算效率更高，在处理海量数据时更具有优势。
-
表 1 容器包含的所有元素
Table 1 All Elements in Container
容器 元素 pV[v5] v1 v4 v1 v2 v2 v*6 v4 v8 v8 v*9 pv[v6] v*5 v2 v2 v3 v3 v7 v*9 v10 v10 v7 pv[v9] v8 v*5 v*6 v10 v8 v11 v11 v12 v12 v10 
下载: 导出CSV
表 2 数据处理结果统计
Table 2 Statistics of Data Processing Results
区域 点云数量 算法 三角形个数 构网时间/s 漏洞个数 漏洞比例/‰ A 128 160 原始算法2 13 918 10.1 1 115 5.2 改进算法 215 325 10.5 321 1.5 B 322 319 原始算法 635 943 32.4 1 004 1.6 改进算法 638 126 33.5 278 0.4
下载: 导出CSV
-
[1] 穆超,余洁,许磊.基于高分辨率遥感影像的DSM建筑物点的提取研究[J]. 武汉大学学报·信息科学版,2009, 34(4):414-417 Mu Chao, Yu Jie, Xu Lei. Extracting Building Points from the DSM Data Combining the High-Resolution Remote Sensing Image[J]. Geomatics and Information Science of Wuhan University, 2009, 34(4):414-417
[2] 王刃, 徐青, 朱新慧.用多种策略从机载LiDAR数据中提取建筑脚点[J]. 武汉大学学报·信息科学版,2008,33(7):688-691 Wang Ren, Xu Qing. Zhu Xinhui. Picking up Foot Prints of Building from Airborne LiDAR Data with Multi-strategies[J]. Geomatics and Information Science of Wuhan University, 2008, 33(7):688-691
[3] Poullis C, You S. Automatic Creation of Massive Virtual Cities[C]. IEEE Virtual Reality Conference, Louisiana, 2009
[4] 李镇洲,张学之.基于倾斜摄影测量技术快速建立城市3维模型研究[J]. 测绘与空间地理信息,2012,35(4):117-119 Li Zhenzhou, Zhang Xuezhi. Research on the Quick Construction of 3D Model of City Based on Oblique Photogrammetric Technique[J]. Geomatics & Spatial Information Technology, 2012, 35(4):117-119
[5] Musialski P, Wonka P, Aliaga D G, et al. A Survey of Urban Reconstruction[J]. Computer Graphics Forum, 2013, 32(6):146-177
[6] Merchan P, Adan A, Salamanca S, et al. Geometric and Colour Data Fusion for Outdoor 3D Models[J]. Sensors, 2012, 12(6):6893-6919
[7] Kuschk G. Large Scale Urban Reconstruction from Remote Sensing Imagery[J]. International Archives of the Photogrammetry, Remote Sensing and Spatial Information Sciences, 2013, XL-5/W1:139-146
[8] Kuschk G. Model-Free Dense Stereo Reconstruction for Creating Realistic 3D City Models[C]. Urban Remote Sensing Event, IEEE Joint, Brazil, 2013
[9] Canciani M, Falcolini C, Saccone M, et al. From Point Clouds to Architectural Models:Algorithms for Shape Reconstruction[J]. International Archives of the Photogrammetry, Remote Sensing and Spatial Information Sciences, 2013, XL-5/W1:27-34
[10] Boissonnat J D. Geometric Structures for Three Dimensional Shape Representation[J]. ACM Transactions on Graphics, 1984, 3(4):266-286
[11] Kolluri R. Provably Good Moving Least Squares[J]. ACM Transactions on Algorithms, 2005, 4(2):782-795
[12] Kazhdan M, Bolitho M, Hoppe H. Poisson Surface Reconstruction[C]. Eurographics Symposium on Geometry Processing, Switzerland, 2006
[13] Amenta N. A New Voronoi-based Surface Reconstruction Algorithm[C]. ACM Conference on Computer Graphics and Interactive Techniques, Atlanta, 1998
[14] Dey T K, Goswami S. Tight Cocone:A Water-Tight Surface Reconstructor[C]. ACM Symposium on Solid Modeling and Applications, Seattle, 2003
[15] Cohen-Steiner D, Da F. A Greedy Delaunay-based Surface Reconstruction Algorithm[J]. The Visual Computer, 2004, 20(1):4-16
[16] Bernardini F. The Ball-Pivoting Algorithm for Surface Reconstruction[J]. Visualization and Computer Graphics, 1999, 5(4):349-359
[17] Li X, Han C Y, Wee W G. On Surface Reconstruction:A Priority Driven Approach[J]. Computer-Aided Design, 2009, 41(9):626-640
[18] Di Angelo L. A New Mesh-Growing Algorithm for Fast Surface Reconstruction[J]. Computer-Aided Design, 2011, 43(6):639-650
-
期刊类型引用(6)
1. 杨根新，王友昆，谢正明. 基于广义EIV模型的矿区高程异常的无缝推估算法. 工程勘察. 2023(08): 46-51 . 
百度学术
2. Jianjun ZHU，Leyang WANG，Jun HU，Bofeng LI，Haiqiang FU，Yibin YAO. Recent Advances in the Geodesy Data Processing. Journal of Geodesy and Geoinformation Science. 2023(03): 33-45 . 必应学术
3. 翁烨，陈丽，王岩. 线性化通用EIV平差模型的正则化解法. 勘察科学技术. 2023(05): 1-5 . 百度学术
4. 戴中东，孟良，高永攀，项伟. 加权整体最小二乘坐标匹配算法在机场道面测量中的应用. 测绘地理信息. 2022(02): 61-66 . 百度学术
5. 翁烨，邵德盛，甘淑. 线性化通用EIV平差模型的岭估计解法. 全球定位系统. 2022(02): 82-89 . 百度学术
6. 翁烨，邵德盛. 病态加权总体最小二乘的广义岭估计解法. 全球定位系统. 2021(06): 84-89 . 百度学术
其他类型引用(2)
计量
- 文章访问数: 1778
- HTML全文浏览量: 65
- PDF下载量: 429
- 被引次数: 8
