本发明涉及一种基于整体最小二乘的附加几何约束的点云建模方法，属于三维激光扫描技术领域。
背景技术：
针对三维激光扫描点云数据处理的后期三维建模阶段，国内外学者进行了大量的研究，并取得相应成果。但是，这些方法往往集中于单张曲面拟合研究，自动或者手工分块拟合，而忽略了模型的整体性。这导致拼接之后的模型存在严重问题，即单独拟合后的面与面之间往往不满足已知的几何约束关系，如平行、垂直、固定角度、固定距离等。该问题在一定程度上限制了重建模型精度的提高。而目前涉及尺寸约束、位置约束等特殊建模理论的研究并不广泛。
依据求解方法不同的操作方式，目前存在的几何约束求解方法可以分为基于代数的求解方法、基于几何约束图的方法和基于人工智能的规则推理。基于代数的求解方法，是将几何约束问题直接转化为非线性方程组，然后采用非线性方程组求解领域的方法进行求解。它最主要的优点是具有很强的通用性和直观性。目前，可行的添加几何约束的建模方法之一是应用惩罚函数法将带约束的优化问题转化为无约束优化问题，再使用Levenberg-Marquardt算法求解。该方法存在有效、稳定性好、计算相对简单的特点。但在研究过程中发现，该方法也存在局限性，即最终解算结果受模型初值影响较大。
技术实现要素：
为解决现有技术的不足，本发明的目的在于提供一种基于整体最小二乘的附加几何约束的点云建模方法，解决目前添加几何约束的建模方法并利用Levenberg-Marquardt算法求解的方法中存在的最终解算结果受模型初值影响较大的局限性。
为了实现上述目标，本发明采用如下的技术方案：
一种基于整体最小二乘的附加几何约束的点云建模方法，其特征是，包括以下步骤：
1)将点云模型转换成线性方程组AX＝L，其中A、X、L分别为系数矩阵、模型参数解、常数项，其中，A为s×n矩阵、X为n×1矩阵、L为s×1矩阵；采用整体最小二乘法拟合得到模型参数初值X0，其中[v1,n+1,...,vn+1,n+1]T为增广矩阵[A L]右奇异向量的最后一列；s为观测值个数即点云数，n为待估模型参数个数；
2)确定初始惩罚因子λ0，惩罚因子变化速度v1；
3)确定初始LM算法调节因子w0，调节因子变化速度v2，变量参数k＝1；
4)根据当前变量参数k的取值，确定模型参数Xk-1，惩罚因子λk-1，LM算法调节因子wk-1；
5)计算目标函数G(Xk-1,λk-1)，梯度▽G(Xk-1,λk-1)，海赛矩阵H(Xk-1,λk-1)；
6)计算目标函数可能的下降方向dk-1：
dk-1＝-{H(Xk-1,λk-1)+wk-1diag[H(Xk-1,λk-1)]}-1▽G(Xk-1,λk-1)；
7)判断G(Xk-1+dk-1,λk-1)＞G(Xk-1,λk-1)是否成立：若成立，则令wk-1＝v2wk-1跳转至步骤6)，直至得到更小值的目标函数，若不成立，则跳转至步骤8)；
8)更换模型参数值和LM算法调节因子：令Xk＝Xk-1+dk-1，wk＝wk-1/v2；
9)检查是否满足迭代终止条件||Xk-Xk-1||＜ε1，ε1为阈值，依据经验值选取；若满足，输出Xk作为极小值点，跳转至步骤10)，否则令k＝k+1，惩罚因子不变λk＝λk-1，跳转至步骤5)；
10)判断P(Xk)＜ε2是否成立，ε2为阈值，依据经验值选取；若成立，停止运算，Xk作为原问题的解输出，否则，令k＝k+1，增大惩罚因子，令λk＝v1λk-1，返回步骤4)。
进一步地，所述步骤2)中初始惩罚因子λ0取1，惩罚因子变化速度v1取10。
进一步地，所述步骤3)中初始LM算法调节因子w0取0.01，调节因子变化速度v2取10。
进一步地，所述步骤5)中：
目标函数G(Xk-1,λk-1)的计算公式为：式中，X为模型参数解，m为面的总个数，wi为第i个拟合面的权值，fi(X)为第i个拟合面的误差函数；Ni表示第i个拟合面数据点的总个数，pj表示第i个拟合面的第j个点，di(X,pj)表示第i个拟合面第j个点到该拟合面的距离；t表示约束条件的个数，hp(x)表示第p个约束条件，且hp(x)＝0为多张拟合面之间的所有等式约束条件。
梯度▽G(Xk-1,λk-1)的计算公式为：
海赛矩阵H(Xk-1,λk-1)的计算公式为：
式中，
本发明所达到的有益效果：本方法融入整体最小二乘算法，提高点云建模过程中初始模型参数的精度；融入附加几何约束的点云建模方法，结合惩罚函数法与LM算法，应用惩罚函数法将带约束的优化问题转化为无约束优化问题，再使用LM算法求解；经过上述两步处理，最终模型参数精度大大提高，相较传统点云建模方法，本方案具有精度高、计算相对简单、适用范围广、运算效率高等特点。
附图说明
图1是本方法的流程图。
具体实施方式
下面结合附图对本发明作进一步描述。以下实施例仅用于更加清楚地说明本发明的技术方案，而不能以此来限制本发明的保护范围。
本方法的具体步骤如下：
1)确定初始模型参数X0。
采用整体最小二乘法拟合得到较高精度的模型参数初值X0过程如下：任意点云模型，可写为形为AX＝L的线性方程，其中，A为s×n矩阵、X为n×1矩阵、L为s×1矩阵，s为观测值个数即点云数，n为待估模型参数个数，非线性方程可通过线性化得到。整体最小二乘的求解是通过奇异值分解来实现的。线性方程AX＝L可改写成[A L][XT-1]T＝0，记增广矩阵C＝[A L]，对C进行奇异值分解C＝U∑VT，其中，∑＝diag(σ1,σ2,...,σn,σn+1),σ1≥σ2≥...≥σn≥σn+1≥0。整体最小二乘解可由增广矩阵右奇异向量的最后一列vn+1得到，解为
2)确定初始惩罚因子λ0，惩罚因子变化速度v1。本发明初始惩罚因子λ0取1，惩罚因子变化速度v1取10。
3)确定初始LM算法调节因子w0，调节因子变化速度v2。本发明初始LM算法调节因子w0取0.01，调节因子变化速度v2取10。
4)确定模型参数Xk-1，惩罚因子λk-1，LM算法调节因子wk-1。第一次令k＝1，对应的值由步骤1)-3)可得，其余由以下步骤循环迭代更换旧值得到新值。
5)计算目标函数G(Xk-1,λk-1)，梯度▽G(Xk-1,λk-1)，海赛矩阵H(Xk-1,λk-1)。
目标函数G(Xk-1,λk-1)的计算公式为：式中，X为所求模型的参数，m为面的总个数，wi为第i个拟合面的权值，fi(X)为第i个拟合面的误差函数；Ni表示第i个拟合面数据点的总个数，pj表示第i个拟合面的第j个点，di(X,pj)表示第i个拟合面第j个点到该拟合面的距离；t表示约束条件的个数，hp(x)表示第p个约束条件，且hp(x)＝0为多张拟合面之间的所有等式约束条件。
这一步应用惩罚函数法将带约束的优化问题转化为无约束优化问题，便于以下步骤中利用LM算法对无约束优化问题的求解。
梯度▽G(Xk-1,λk-1)的计算公式为：
海赛矩阵H(Xk-1,λk-1)的计算公式为
式中，
简单来说，梯度函数就是对目标函数一次求导，海塞矩阵就是对目标函数二次求导。
6)计算目标函数可能的下降方向dk-1。
可能的下降方向dk-1的计算公式为dk-1＝-{H(Xk-1,λk-1)+wk-1diag[H(Xk-1,λk-1)]}-1▽G(Xk-1,λk-1)，式中，▽G(Xk-1,λk-1)、H(Xk-1,λk-1)分别为步骤5)计算所得梯度、海塞矩阵，wk-1为步骤4)计算所得LM算法调节因子。
事实上，该步骤中的LM算法其实是最速下降法和牛顿法的混合，在接近最小值之前主要，而在接近最小值时逐步切换为最速下降法。该过程由通过调整惩罚因子λk-1来达到，调节因子变化速度为v2。
7)判断G(Xk-1+dk-1,λk-1)＞G(Xk-1,λk-1)是否成立。
若成立，更新后目标函数值大于原来的目标函数值，则令wk-1＝v2wk-1转步骤6)，继续循环迭代步骤6)～7)，若不是，更新后目标函数值不大于等于原来的目标函数值，转步骤8)。
8)更换模型参数值、LM算法调节因子。计算公式为Xk＝Xk-1+dk-1，wk＝wk-1/v2。
9)检查是否满足迭代终止条件||Xk-Xk-1||＜ε1。其中ε1为阈值，依据经验值选取。若满足，将Xk作为该方程的解，进入步骤10)，否则令k＝k+1，λk＝λk-1即惩罚因子不变，转步骤5)进入循环迭代，这个过程是LM算法求解目标函数有约束优化问题的过程。
10)判断P(Xk)＜ε2是否成立。其中ε2为阈值，依据经验值选取。若成立，说明最终模型参数所建模型较好满足已知几何约束条件，整个循环停止，Xk是最终求得模型参数，否则，令k＝k+1，λk＝v1λk-1即增大惩罚因子，返回步骤4)。
这个过程是利用惩罚函数法的求解过程。
以上所述仅是本发明的优选实施方式，应当指出，对于本技术领域的普通技术人员来说，在不脱离本发明技术原理的前提下，还可以做出若干改进和变形，这些改进和变形也应视为本发明的保护范围。