本发明涉及地铁检测及测量方法技术领域，具体来说是一种用于地铁检测及测量的3D扫描方法，包括如下步骤：扫描获得隧道二维点云数据，并同步采集里程数据；将二维点云数据及里程数据相匹配，以获得隧道的三维点云数据；对隧道进行切片，并对切片进行预处理；进行椭圆拟合；本发明同现有技术相比，其优点在于：解决了现有技术在隧道管壁上布设监测点，利用全站仪、断面仪、收敛仪进行监测的方法费时费力的问题，实现了对隧道的整体监测和形变监测，通过本发明方法，对地铁隧道及结构进行检测和测量，能够全面的掌握隧道的整体情况，并留存数字档案，并能实现大范围的隧道收敛测量工作以及隧道轴线测量工作，有效增加了隧道监护测量的效率和准确性。
1.一种用于地铁检测及测量的3D扫描方法，其特征在于所述的方法包括如下步骤：a.扫描获得隧道二维点云数据，并同步采集里程数据；
b.将二维点云数据及里程数据相匹配，以获得隧道的三维点云数据；
e.对隧道进行变形分析，并根据三维点云数据建立隧道平铺图；
所述的步骤a包括用以检验里程数据的里程粗差探测算法，所述的粗差探测算法具体如下：设误差方程为：
式中：V是n维观测值残差向量；A是n×t阶系数矩阵； 是t维参数估值向量；L是n维观测向量，可得最小二乘解为：可求得观测值残差向量V，及其对应的权逆阵为：-1 T -1 T
QV的对角线元素qv1,qv2,...,qvn是观测值残差v1,v2,...,vn的权倒数，由于观测值独立，有qvi＝1-Ai(ATA)-1AiT其中，Ai是观测值系数矩阵A的第i行，观测值残差vi的方差为：标准化残差为：
2.如权利要求1所述的一种用于地铁检测及测量的3D扫描方法，其特征在于所述的步骤a还包括里程粗差剔除方法，所述的里程粗差剔除方法具体如下：a1.根据扫描仪的扫描线数量，将所述的隧道顺序划分为若干环，判断起始里程sDis与终点里程eDis大小，若起始里程sDis>终点里程eDIs,表示从大环推扫至小环；若终点里程eDIs>起始里程sDis,表示从小环推扫至大环；
a2.根据从大环至小环或从小环至大环的顺序，依次读取采集得到的里程信息并匹配每环的起点里程和终点里程；
a3.根据里程数据的起点时间和终点时间，计算扫描总时间，并计算每根扫描线对应的起点时间和终点时间；
a4.对每组数据都进行如下剔除鉴别处理:除了最后和第一个里程，将不在(0.0,max(sDis,eDis)+10)区间内的里程数据删除；
a41.用第一个里程后续的10个里程数据，使用粗差探测算法，估计第一个里程数值大小，判断第一个里程数据是否正确，若不正确，使用估算数据替换，最后一个数据作相同的数据处理；
a42.除了最后一个里程和第一个里程，逐次对其它里程数据判断，判断依据如下：若上一个里程数据不为0.0，下一个里程数据也不为0.0，而该里程数据为0.0，则删除该里程数据；
若上一个里程数据为0.0，下一个里程数据也为0.0，而该里程数据不为0.0，则删除该里程数据；
a43.将所有测段数据拼接成为一个整体,从第二个历元到倒数第二个历元进行速度判断：与前一个历元或后一个历元计算速度，若速度在0-3m/s保存该历元，否则删除；
a44.从第二个历元到倒数第二个历元按照从小到大判断，若历元不满足要求，则删除历元数据。
3.如权利要求1所述的一种用于地铁检测及测量的3D扫描方法，其特征在于所述的对隧道进行切片的方法具体如下：根据采集得到的隧道点云数据的区间长度，以及区间内对应的环数，计算每一环的长度，并将每一环均分为若干切片，并使切片的厚度趋近于10cm，计算每个切片的起点和终点，并依此将每个切片的起点和终点及相应的三维点云数据保存至每个切片的文件数据中。
4.如权利要求1所述的一种用于地铁检测及测量的3D扫描方法，其特征在于所述的对切片进行预处理的方法具体如下：将获得的每个切片的文件数据的数据量整除100kb，获得若干个100kb的文件数据，不足100kb的文件数据视为不存在，记切片文件包含N个100kb的文件数据，对每个切片文件依次进行处理为：若N>9,则将切片文件数据分别存取在其他N/2个文件下，按照顺序读取切片每行文件，并根据行号，依次保留到对应的文件下,即实现数据均匀分配,并删除切片原始文件数据；
5.如权利要求1所述的一种用于地铁检测及测量的3D扫描方法，其特征在于所述的椭圆拟合方法包括对测量点平面坐标的处理方法：设测量点拟合平面方程如下：
式中X＝(x y h)T为任意点坐标，e＝(a b c)T为平面法线单位矢量，d为原点至平面的距离；
取e0＝(1.0 0.0 0.0)T，d0＝0.0，对所有断面点XiT＝(xi yi hi)T列出误差方程：vi＝XiTe+d
建立新的坐标系oo_xxyyhh，原点oo为所有投影点坐标的均值(xm ym hm)T，hh为平面法线方向，xx在平面内，与h轴共面，故：ehh＝(a b c)T
若c≠0，则平面与h轴的交点为 exx为以下向量的单位化：yy与另外两轴正交：
由平面坐标系中坐标求测量坐标系中坐标的转换关系为：
6.如权利要求5所述的一种用于地铁检测及测量的3D扫描方法，其特征在于在得到测量坐标系中坐标后进行椭圆拟合，所述的椭圆拟合方法还包括特征根法椭圆拟合，采用平面坐标系内的平面坐标X＝(xx yy)T，拟合椭圆；
定义椭圆坐标系Oe-x_ey_e，原点Oe为椭圆中心，x_e为长半轴方向，y_e为短半轴方向，y_e方向处在x_e顺时针转过90°的方向；
测量坐标X与椭圆坐标Xe的转换关系表示为：Xe＝X0+R(α)X
认为各测量点坐标不满足椭圆方程的部分为残差，列出误差方程：vi＝XeiTΛXei-1
将Xe＝X0+R(α)X代入误差方程，vi＝(X0+R(α)Xi)TΛ(X0+R(α)Xi)-1vi＝X0TΛX0+2XiTRT(α)ΛX0+XiTRT(α)ΛR(α)Xi-1线形化得：
li＝1-(X0+RXi)TΛ(X0+RXi)＝XeTΛXe对所有测量点列出误差方程，在最小二乘 的条件下组成法方程，求解参数改正数，并迭代至改正数收敛，即得到椭圆的长短半轴a、b，测量坐标与椭圆坐标之间的平移量X0和旋转角α；
拟合时根据残差剔除测量粗差点，若某观测量改正残差大于限差，则视该观测量为粗差点，剔除粗差后根据最小二乘法则重新解算待估参数；
由于隧道切面接近于圆，为了解的稳定，设定α＝90°线性化为：
所述的椭圆拟合方法还包括至焦点距离和不变法椭圆拟合法，定义vi为各点至椭圆的垂直距离：vi＝hi＝f(Xi,a,b,x0,y0,α)参数取初值a0,b0,x00,y00,α0，列出误差方程：
7.如权利要求6所述的一种用于地铁检测及测量的3D扫描方法，其特征在于所述的对隧道进行变形分析的方法包括：e1.对管片进行分块拟合，以确定管片的变形参数：定义θ为管片绕其中点P的旋转角，D为本管片沿OP方向的平移量；
以XeP表示P点的椭圆坐标系坐标，绕P点旋转θ，再沿OP方向平移D后，该管片上各测量点的椭圆坐标Xei变为X'ei：其中αOP为P点在椭圆坐标系内的方位角；
按旋转和平移后的坐标与拟合椭圆之间的偏差列出误差方程：vi＝X'eiTΛX'ei-1
在最小二乘 的条件下组成法方程，求解待估参数，对参数改正数进行迭代至收敛，即可解出旋转角θ和平移量D。
8.如权利要求7所述的一种用于地铁检测及测量的3D扫描方法，其特征在于所述的对隧道进行变形分析的方法还包括：e2.对隧道进行变形分析的方法包括确定管片端点的变形量的方法：在椭圆坐标系内，求出管片的起点坐标、终点坐标及与椭圆的交点X'eS、X'eE；
两交点转过旋转θ角、平移D后，交点坐标为下式中的XeS、XeE：比较X'eS和XeS、X'eE和XeE，即得到该管片的变形量。
9.如权利要求7所述的一种用于地铁检测及测量的3D扫描方法，其特征在于所述的对隧道进行变形分析的方法还包括：e3.确定点的变形量的方法：对面坐标系方位角为α的点i，先取一个距原点10m的临时点，其平面坐标为：将其转换至椭圆坐标：
在椭圆坐标系内，原点与临时点连线上点的坐标为(kxet kyet) ，代入椭圆方程，求得k:从而得到i点的椭圆坐标X'ei，该点是Xei绕中点旋转平移得到的：i点形变量defi为：
一种用于地铁检测及测量的3D扫描方法
[技术领域]
[0001] 本发明涉及地铁检测及测量方法技术领域，具体来说是一种用于地铁检测及测量的3D扫描方法。[背景技术]
[0002] 在当今以计算机技术为依托的信息时代，人们对空间三维信息的需求更加迫切。基于测距、测角的传统工程测量方法，在理论、设备和应用等诸多方面都已相当成熟，新型的全站仪可以完成特征目标的高精度测量，GPS可以全天候、一天24小时精确定位全球任何位置的三维坐标，但它们多用于稀疏目标点的高精度测量。3D激光扫描技术的出现和发展为空间三维信息的获取提供了全新的技术手段，为信息数字化发展提供了必要的生存条件。3D激光扫描技术克服了传统测量技术的局限性，采用非接触主动测量方式直接获取高精度三维数据，能够对任意物体进行扫描，且没有白天和黑夜的限制，快速将现实世界的信息转换成可以处理的数据。
[0003] 而隧道形变监测对隧道施工及运营中的安全具有重要意义，现有的监测方法主要是在隧道管壁上布设监测点，利用全站仪、断面仪、收敛仪等进行监测，但这些方法不仅费时费力，更重要的是难以实现对隧道的整体监测。因此，本申请采用3D扫描技术实现对地铁隧道及结构的检测和测量，全面的掌握隧道的整体情况，并留存数字档案，并实现大范围的隧道收敛测量工作以及隧道轴线测量工作等。[发明内容]
[0004] 本发明的目的在于解决现有技术的不足，提供一种用于地铁检测及测量的3D扫描方法，快速、准确、有效地获取隧道的三维信息，优化数据的处理方法，以得到可靠的隧道变形分析结果。
[0005] 为了实现上述目的，设计一种用于地铁检测及测量的3D扫描方法，包括如下步骤：
[0006] a.扫描获得隧道二维点云数据，并同步采集里程数据；
[0007] b.将二维点云数据及里程数据相匹配，以获得隧道的三维点云数据；
[0008] c.对隧道进行切片，并对切片进行预处理；
[0009] d.进行椭圆拟合；
[0010] e.对隧道进行变形分析，并根据三维点云数据建立隧道平铺图。
[0011] 所述的步骤a包括用以检验里程数据的里程粗差探测算法，所述的粗差探测算法具体如下：
[0012] 设误差方程为：
[0013]
[0014] 式中：V是n维观测值残差向量；A是n×t阶系数矩阵；x是t维参数估值向量；L是n维观测向量。参数最小二乘解为：
[0015]
[0016] 可求得观测值残差向量V，及其对应的权逆阵为：
[0017] QV＝P-1-A(ATA)-1AT
[0018] QV的对角线元素qv1,qv2,...,qvn是观测值残差v1,v2,...,vn的权倒数。由于观测值独立，有
[0019]
[0020] 其中，Ai是观测值A的第i行。观测值vi的方差为：
[0021]
[0022] 标准化残差为：
[0023]
[0024] 当 大于给定的限值时即认定观测值Li存在粗差。
[0025] 所述的步骤a还包括里程粗差剔除算法，所述的里程粗差剔除算法具体如下：
[0026] a1.根据扫描仪的扫描线数量，将所述的隧道顺序划分为若干环，判断起始里程sDis与终点里程eDis大小，若起始里程sDis>终点里程eDIs,表示从大环推扫至小环；若终点里程eDIs>起始里程sDis,表示从小环推扫至大环；
[0027] a2.根据从大环至小环或从小环至大环的顺序，依次读取采集得到的里程信息；
[0028] a3.根据里程数据的起点时间和终点时间，计算扫描总时间，并计算每根扫描线对应的起点时间和终点时间；
[0029] a4.对每组数据都进行如下剔除鉴别处理:
[0030] 除了最后和第一个里程，将不在(0.0,max(sDis,eDis)+10)区间内里程删除；
[0031] a41.用第一个里程后续的10个里程数据，使用粗差探测算法，估计第一个里程数值大小，判断第一个里程数据是否正确，若不正确，使用估算数据替换，最后一个数据作相同的数据处理；
[0032] a42.除了最后一个里程和第一个里程，逐次对其它里程数据判断，判断依据如下：
[0033] 若上一个里程数据不为0.0，下一个里程数据也不为0.0，而该里程数据为0.0，则删除该里程数据；
[0034] 若上一个里程数据为0.0，下一个里程数据也为0.0，而该里程数据不为0.0，则删除该里程数据；
[0035] a43.将所有测段数据拼接成为一个整体,从第二个历元到倒数第二个历元进行速度判断：与前一个历元或后一个历元计算速度，若速度在0-3m/s保存该历元，否则删除；
[0036] a44.从第二个历元到倒数第二个历元按照从小到大判断，若历元不满足要求，则删除历元数据。
[0037] 所述的对隧道进行切片的方法具体如下：根据采集得到的隧道点云数据的区间长度，以及区间内对应的环数，计算每一环的长度，并将每一环均分为若干切片，并使切片的厚度趋近于10cm，计算每个切片的起点和终点，并依此将每个切片的起点和终点及相应的三维点云数据保存至每个切片的文件数据中。
[0038] 所述的对切片进行预处理的方法具体如下：将获得的每个切片的文件数据的数据量整除100kb，获得若干个100kb的文件数据，不足100kb的文件数据视为不存在，记切片文件包含N个100kb的文件数据，对每个切片文件依次进行处理为：
[0039] 若N>9,则将切片文件数据分别存取在其他N/2个文件下，按照顺序读取切片每行文件，并根据行号，依次保留到对应的文件下,即实现数据均匀分配,并删除切片原始文件数据；
[0040] 若N<1,直接将原始数据删除。
[0041] 所述的椭圆拟合方法包括对对测量点平面坐标的处理方法：
[0042] 设测量点拟合平面方程如下：
[0043] XTe+d＝0
[0044] 式中X＝(x y h)T为任意点坐标，e＝(a b c)T为平面法线单位矢量，d为原点至平面的距离；
[0045] 取e0＝(1.0 0.0 0.0)T，d0＝0.0，对所有断面点XiT＝(xi yi hi)T列出误差方程：
[0046] vi＝XiTe+d
[0047]
[0048]
[0049] 则各测量点至平面的距离ti为：
[0050] ti＝XiTe+d
[0051] 所有ti应满足|ti|＜ε；
[0052] 得到各测量点在平面上的投影点坐标为：
[0053] Xi-eti
[0054] 建立新的坐标系oo_xxyyhh，原点oo为所有投影点坐标的均值(xm ym hm)T，hh为平[0055] 面法线方向，xx在平面内，与h轴共面，故：
[0056] ehh＝(a b c)T
[0057] 若c＝0，则，exx＝(0 0 1)T；
[0058] 若c≠0，则平面与h轴的交点为 exx为以下向量的单位化：
[0059]
[0060] yy与另外两轴正交：
[0061] eyy＝exx×ehh
[0062] 因此，平面坐标与测量坐标的关系为：
[0063]
[0064] 式中的：
[0065]
[0066] 由平面坐标系中坐标求测量坐标系中坐标的转换关系为：
[0067]
[0068] 本发明同现有技术相比，其优点在于：解决了现有技术在隧道管壁上布设监测点，利用全站仪、断面仪、收敛仪进行监测的方法费时费力的问题，并且实现了对隧道的整体监测和形变监测，通过本发明方法，对地铁隧道及结构进行检测和测量，能够全面的掌握隧道的整体情况，并留存数字档案，并能实现大范围的隧道收敛测量工作以及隧道轴线测量工作，采用椭圆拟合方法去除了数据中的误差值，易于对隧道情况进行更好的监测，有效增加了隧道监护测量的效率和准确性。[附图说明]
[0069] 图1是一实施例中本发明方法的流程示意图；
[0070] 图2是一实施例中本发明方法中测量点的平面坐标的处理示意图；
[0071] 图3是一实施例中本发明特征根法椭圆拟合的示意图；
[0072] 图4是一实施例中本发明检测管片的变形参数的示意图；
[0073] 图5是一实施例中本发明方法建立隧道平铺图的示意图。[具体实施方式]
[0074] 下面结合附图对本发明作进一步说明，这种方法的原理对本专业的人来说是非常清楚的。应当理解，此处所描述的具体实施例仅仅用以解释本发明，并不用于限定本发明。
[0075] 本发明包括处理器、存储介质和数据采集装置，所述的数据采集装置包括移动小车、3D扫描仪和里程仪，在移动小车上设有3D扫描仪和里程仪，并使移动小车沿地铁轨道向前行进，通过扫描仪和里程仪同步进行数据采集获得相应的二维点云数据和里程数据，将存储于存储介质中，由处理器对数据进行处理和分析，或直接将二维点云数据和里程数据传输给处理器进行运算分析并最终分析得到隧道的变形情况，以对隧道进行安全检测并进行及时、合理的维护，解决现有的利用全站仪、断面仪、收敛仪等进行监测，费时费力难以实现对隧道的整体监测的问题。参见图1，图1左侧为现有技术的一流程实例，图1右侧为本发明方法的一流程示意，本实施方式对数据进行采集和处理的具体步骤如下：
[0076] a.通过扫描仪扫描获得隧道二维点云数据，并在扫描的同时，通过里程仪同步采集里程数据，并将扫描和采集得到的数据传输至存储介质中或直接传输给处理器进行运算分析；
[0077] 扫描仪采用螺旋扫模式，激光镜片能旋转，而扫描仪自身固定，也就是扫描仪只能采集二维轮廓数据且无里程信息。激光镜片每旋转一圈，获得一根扫描线数据。推扫数据由一根根扫描线组合而成。扫描仪默认在停止扫描事件被触发时：直到当下扫描线数据采集结束时，扫描仪停止工作。每根扫描线的点数、质量以及采集速度，对于不同用户输入参数不同。若使用扫描仪用户界面设置扫描参数，要设置分辨率和质量；若通过Faro LS SDK接口函数设置参数，需设置参数分辨率、观测速度、去噪等级。界面和程序参数设置虽然个数不同，但是各组参数之间设置等价。
[0078] 本实施方式的实际扫描过程中，扫描仪的质量选择3X或4X，分辨率选择1/4，1/5或1/8。质量若选择2X或1X，扫描仪的旋转速度极快，点云质量较低。于工程而言，分辨率选择过高或过低，将致使测量数据冗余或稀疏。
[0079] 在扫描仪进入螺旋扫模式后，扫描仪的水平度盘卡死，并默认为0度。地铁隧道情况复杂，在有些转弯路段，铁轨的左右具有一定超高。超高使得小车左右倾斜，而扫描仪会随着轨道倾斜，发生水平方向转动。由于转动的发生，采集数据的质量严重下降，后期数据处理会发现明显错误。提出的解决办法是：使用两根螺丝，左右各一根固定扫描仪，使其不会发生转动。办法虽简单，却很好解决了推扫过程中，扫描仪绕水平度盘转动的问题。
[0080] 数据采集还需解决一个关键问题：推扫速度的大小。推扫太慢容易导致采集点密度过大，推扫时间花费增多且后续数据处理时间递增。经过多次实验发现，扫描线之间间隔越大，每根扫描线点越少，导致切片计算点变少，获得的计算结果精度越低。可结合每条扫描线费时表，计算一定推扫速度下扫描线间隔。最后发现，若选择分辨率1/4、质量为3X，推扫速度低于1.5m/s时，可以计算获得稳定的计算结果。对于其他参数组合，可使用此组参数作为设计依据。
[0081] 里程仪提供扫描时间序列下的里程数据，提供扫描隧道的第三维信息。扫描仪根据采集点时间和里程数据，内插出各点里程数据。里程仪采集数据的频率越高，意味着相邻里程数据时间间隔越小，提供推扫的运动信息越全面，内插出各点里程信息精度越高。
[0082] 里程仪与电脑之间采用串口通讯，数据传输会有极小延迟，但可忽略其对里程结果的影响。当里程仪记录距离过大时，易形成误差累积，应适当距离清零里程。里程仪适宜采用速度为10个/s，若再增大里程仪的频率，则产生的里程数据出错率大大增加。里程仪在采集数据过程中，由于USB与电脑连接不稳定，会出现异常提醒,若不造成程序运行的停止，不会影响里程数据记录。综合考虑，里程仪在采集数据过程中，设置采集速度为10个/s，每隔600个里程数据清零一次。
[0083] 数据采集程序需具有输入扫描仪参数、里程参数及同步采集数据的功能。扫描线数决定了本次螺旋扫需完成的扫描线数，不低于100条。直线段或者曲线段是为了记录扫描测段的形式，起始和终点环号是数据处理需要输入的参数。在推扫过程中，扫描数据的存取方式有两种：存在扫描仪内存卡中或实时传入电脑存储器中，可供选择。扫描仪连接信息反馈了异常及记录情况，异常信息记录了里程仪与电脑串口通讯错误。里程仪清零和接收频率参数，在程序内部直接固定，分别为600个里程数据清零一次和0.1s接受一个里程数据。
[0084] 使用推扫速度、距离递增等已知条件，结合粗差探测的算法获得纯净的里程仪数据。认为推扫速度不会超过3m/s，超过部分可认为是粗差数据。由于里程数据每600个里程数据清零一次，每段起点里程和终点里程数据至关重要，发生错误则导致里程数据大面积偏差，对后续数据处理是致命的。因此，使用粗差探测算法检测，每段首尾里程的正确性。
[0085] 所述的里程粗差探测算法，用以检验里程数据，处理器根据预先设置的里程粗差探测算法程序来检测里程数据中是否存在粗差，所述的粗差探测算法具体如下：
[0086] 设误差方程为：
[0087]
[0088] 式中：V是n维观测值残差向量；x是t维参数估值向量；L是n维观测向量；A是n×t阶系数矩阵，例如，可以认为每次观测历元t对应的里程数据l均可按照如下多项式展开：
[0089] l＝a0+a1t+a2t2+a3t3
[0090] 则系数矩阵A中对应行可表示为：(1,t,t2,t3)。
[0091] 参数最小二乘解为
[0092]
[0093] 式中P是观测值权矩阵，可认为是里程l的权矩阵，此处为单位权矩阵。
[0094] 可求得观测值残差向量V，及其对应的权逆阵为
[0095] QV＝P-1-A(ATA)-1AT
[0096] QV的对角线元素qv1,qv2,...,qvn是观测值残差v1,v2,...,vn的权倒数。由于观测值独立，有
[0097]
[0098] 其中，Ai是观测值系数矩阵A的第i行。观测值残差vi的方差为：
[0099]
[0100] 标准化残差为：
[0101]
[0102] 当观测值没有粗差时， 是服从标准正态分布的随机变量。当 大于给定的限值(例如2或3)时即认定观测值Li存在粗差。
[0103] 若处理器探测发现粗差，则需要对粗差进行剔除，处理器根据预先设置的里程粗差剔除算法程序对粗差进行剔除，里程粗差的剔除方法具体如下：
[0104] a1.根据扫描仪的扫描线数量，将所述的隧道顺序等距划分为若干环，判断起始里程sDis与终点里程eDis大小，若起始里程sDis>终点里程eDIs,表示从大环推扫至小环；若终点里程eDIs>起始里程sDis,表示从小环推扫至大环；
[0105] a2.根据从大环至小环或从小环至大环的顺序，依次读取采集得到的里程信息并匹配每环的起点里程和终点里程；
[0106] 此时，还需判断里程数据的采集时间是否跨越了两天，若跨越了，后续处理里程时间需要判断其具体是哪一天的时间点，以确保数据的时间顺序无误；
[0107] a3.根据里程数据的起点时间和终点时间，计算扫描总时间，并计算每根扫描线对应的起点时间和终点时间；
[0108] 由于每个600个数据里程清零，而长距离推扫，包含多组600里程组。因此，首先计算分组数，对每一组数据开空间，并将每组数据分别存取。
[0109] a4.对每组数据都进行如下剔除鉴别处理:
[0110] 由于里程仪记录的数据中，有很多是错误的记录。若某一个里程数据大于测段的大小，对于超过这一区间的里程数据，肯定是错误的，必须删除，此处人为添加的经验距离范围10，以免误删除。除了最后和第一个里程，将不在(0.0,max(sDis,eDis)+10)区间内的里程数据删除。
[0111] a41.用第一个里程后续的10个里程数据，使用粗差探测算法，估计第一个里程数值大小，判断第一个里程数据是否正确，若不正确，使用估算数据替换，最后一个数据作相同的数据处理。里程数大小估计方法就是通过多项式估计，通过先求解出所述粗差探测算法的多项式中参数，然后仅需要输入里程时刻就可以估计对应时刻下里程数据。
[0112] a42.除了最后一个里程和第一个里程，逐次对其它里程数据判断，判断依据如下：
[0113] 若上一个里程数据不为0.0，下一个里程数据也不为0.0，而该里程数据为0.0，则删除该里程数据；
[0114] 若上一个里程数据为0.0，下一个里程数据也为0.0，而该里程数据不为0.0，则删除该里程数据；
[0115] a43.将所有测段数据拼接成为一个整体,从第二个历元到倒数第二个历元进行速度判断：与前一个历元或后一个历元计算速度，若速度在0-3m/s保存该历元，否则删除，所述的历元是指每个里程数据所对应的时刻。
[0116] a44.从第二个历元到倒数第二个历元按照从小到大判断，若历元不满足要求，则删除历元数据。
[0117] b.剔除里程数据的粗差后，处理器根据预设算法将二维点云数据及里程数据相匹配，以获得隧道的三维点云数据；
[0118] 点云数据是以扫描线为基础的，扫描线数总是一个整数，假设扫描仪分辨率为1/4,质量为3X，那么每根扫描线需花费42ms，里程仪采集一个里程数据需要花费0.1s,鉴于此，只需要依次计算每根扫描仪起点里程以及对应的速度，就可以完成点云数据三维里程的填充，从而通过二维点云数据及里程数据获得三维点云数据。
[0119] c.对隧道进行切片，并对切片进行预处理；
[0120] c1.所述的对隧道进行切片的方法具体如下：处理器根据采集得到的隧道点云数据的区间长度，以及区间内对应的环数，计算每一环的长度，并将每一环均分为若干切片，并使切片的厚度趋近于10cm，计算每个切片的起点和终点，并依此将每个切片的起点和终点及相应的三维点云数据保存至每个切片的文件数据中，所述的文件数据存储于存储介质中。
[0121] c2.所述的对切片进行预处理的方法具体如下：对数据量过大或过小的切片文件数据进行预处理，具体的处理方法是：所述的对切片进行预处理的方法具体如下：将获得的每个切片的文件数据的数据量整除100kb，获得若干个100kb的文件数据，不足100kb的文件数据视为不存在，记切片文件包含N个100kb的文件数据，对每个切片文件依次进行处理为：
[0122] c21.若N>9,则将切片文件数据分别存取在其他N/2个文件下，按照顺序读取切片每行文件，并根据行号，依次保留到对应的文件下,即实现数据均匀分配,并删除切片原始文件数据；
[0123] c22.若N<1,直接将原始数据删除。
[0124] d.分片后，进行椭圆拟合，以去除误差数据，提高后续分析处理的效率和结果的准确性，以更好地对地铁隧道情况进行有效监测。本实施例采用了多种拟合方法，一般二次曲线法剔除了大部分粗差点，并得到拟合参数，特征根法是基于二次曲线法基础上进行，并不是严格按照点到椭圆垂直距离最小列式，但算法高效且与按点到椭圆距离最小列式结果相同，特征值法严格按照点到椭圆的垂直距离最小列式，可获得最小二乘最严密解；至焦点和不变法是从椭圆上任意点至两个焦点的距离和为常数的定义出发，列误差方程求解，可获得椭圆参数解，算法形象生动的表达了椭圆定义，具有显著的创造性。在进行一般二次曲线法拟合后，使用特征根法椭圆拟合、数值导数法椭圆拟合和至焦点距离和不变法椭圆拟合三种方法的至少一种进行椭圆拟合，最好是三种拟合方法全部使用，可提供比较依据。即在使用一般二次曲线法拟合剔除误差点的基础上，再使用特征值法、数值导数法和至焦点和不变法拟合椭圆，若得到的结果在亚毫米级变换，可认为算的结果基本相同。
[0125] d1.参见图2，处理器对三维点云数据进行处理，以得到测量点的平面坐标。
[0126] 设测量点拟合平面方程如下：
[0127] XTe+d＝0
[0128] 式中X＝(x y h)T为任意点坐标，e＝(a b c)T为平面法线单位矢量，d为原点至平面的距离。
[0129] 取近似值e0＝(1.0 0.0 0.0)T，d0＝0.0，对所有断面点XiT＝(xi yi hi)T列出误差方程：
[0130] vi＝XiTe+d
[0131]
[0132]
[0133] 各测量点至平面的距离ti为：
[0134] ti＝XiTe+d
[0135] 所有ti应满足|ti|＜ε。式中，ε为定值，本实施例根据工程经验选择ε为5cm，对于不满足该式的ti进行剔除，不再进行后续拟合。
[0136] 各测量点在平面上的投影点坐标为：
[0137] Xi-eti
[0138] 建立新的坐标系oo_xxyyhh，原点oo为所有投影点坐标的均值(xm ym hm)T，hh为平面法线方向，xx在平面内，与h轴共面。故：
[0139] ehh＝(a b c)T
[0140] 若c＝0，则，
[0141] exx＝(0 0 1)T
[0142] 若c≠0，则平面与h轴的交点为 exx为以下向量的单位化：
[0143]
[0144] yy与另外两轴正交：
[0145] eyy＝exx×ehh
[0146] 因此，平面坐标与测量坐标的关系为：
[0147]
[0148] 式中的：
[0149]
[0150] 由平面坐标系中坐标求测量坐标系中坐标的转换关系为：
[0151]
[0152] d2.通过一般二次曲线标准化法椭圆拟合法剔除粗差点并得到拟合参数：
[0153] 以平面坐标X＝(xx yy)T表示的一般二次曲线方程为：
[0154] a0+a1xx+a2yy+a3xx2+a4xxyy+a5yy2＝0
[0155] 写成矩阵形式：
[0156]
[0157] 为了避免对平面坐标X＝(xx yy)T数值过大，做以下正规化：
[0158]
[0159] 式中的 是坐标分量均值，Δx、Δy分别是坐标分量最大值与最小值之差的一半。
[0160] 正规化后，x'i与y'i在尺度上不一致，为使尺度一致，也可使Δx、Δy取一样的值，如Δx、Δy均取为其中大的一个。坐标分量x'i与y'i的值均落在区间[-1 1]内。将二次曲线方程写为：
[0161] 1+a'1x'+a'2y'+a'3x'2+a'4x'y'+a'5y'2＝0
[0162] 列出误差方程为：
[0163] vi＝x'iδa'1+y'iδa'2+x'i2δa'3+x'iy'iδa'4+y'i2δa'5-l'i
[0164] 式中l'i＝-1-a'1x'i-a'2y'i-a'3x'i2-a'4x'i y'i-a'5y'i2
[0165] 参数a'1,a'2,...,a'5的迭代初值可以取为0，由残差确定单位权重误差，将残差大于三倍单位权重误差的点剔除，迭代求解收敛后，所有点残差均小于三倍单位权重误差。再重新进行计算，直到。再反求(xx yy)T对应的系数。将一般二次曲线方程式写为：
[0166]
[0167] 将X和X'的关系写为：
[0168]
[0169] 式中：
[0170]
[0171] 故有：
[0172]
[0173]
[0174] 式中：
[0175]
[0176] 因此，平面二次曲线方程的系数a0,a1,...,a5为：
[0177]
[0178] 展开为：
[0179]
[0180] d3.参见图3，采用特征根法椭圆拟合法，利用平面坐标系内的平面坐标X＝(xx yy)T，拟合椭圆，剔除粗差点并得到拟合参数；
[0181] 椭圆坐标系原点在椭圆中心，x_e为长半轴方向，y_e为短半轴方向，圆坐标(x_cir y_cir)T是椭圆坐标饶短半轴旋转得到。
[0182] 定义椭圆坐标系Oe-x_ey_e，原点Oe为椭圆中心，x_e为长半轴方向，y_e为短半轴方向，y_e方向处在x_e顺时针转过90°的方向。
[0183] 椭圆在椭圆坐标系中的方程为：
[0184] XeTΛXe＝1
[0185] 式中 为椭圆坐标， a、b为椭圆的长短半轴。
[0186] 测量坐标X与椭圆坐标Xe的转换关系表示为：
[0187] Xe＝X0+R(α)X
[0188] 式中 为平移量， 为旋转矩阵。
[0189] 认为各测量点坐标不满足椭圆方程的部分为残差，列出误差方程：
[0190] vi＝XeiTΛXei-1
[0191] 将Xe＝X0+R(α)X代入误差方程，
[0192] vi＝(X0+R(α)Xi)TΛ(X0+R(α)Xi)-1
[0193] vi＝X0TΛX0+2XiTRT(α)ΛX0+XiTRT(α)ΛR(α)Xi-1
[0194] 线形化得：
[0195]
[0196] 其中各项导数和常数项为：
[0197]
[0198]
[0199]
[0200]
[0201]
[0202] li＝1-(X0+RXi)TΛ(X0+RXi)＝XeTΛXe
[0203] 对所有测量点列出误差方程，在最小二乘 的条件下组成法方程，求解参数改正数，并迭代至改正数收敛，即得到椭圆的长短半轴a、b，测量坐标与椭圆坐标之间的平移量X0和旋转角α；
[0204] 拟合时可以根据残差剔除测量粗差点，若某观测量改正残差大于限差(如限差为3σ)，则视该观测量为粗差点，剔除粗差后根据最小二乘法则重新解算待估参数。
[0205] 由于隧道切面接近于圆，为了解的稳定，加条件：
[0206] α＝90°
[0207] 线性化为：
[0208] δα＝-α0
[0209] 其权为：
[0210] pα为法方程其余四个对角元的均值。
[0211] 测量点平面坐标与椭圆坐标的关系：
[0212]
[0213] 圆坐标与椭圆坐标的关系：
[0214]
[0215] d4.采用数值导数法椭圆拟合剔除粗差点并得到拟合参数，前两个椭圆拟合方法VTpV为代数残差，更合理的是将vi定义为各点至椭圆的垂直距离；
[0216] vi＝hi＝f(Xi,a,b,x0,y0,α)
[0217] 参数取初值a0,b0,x00,y00,α0，
[0218]
[0219]
[0220]
[0221]
[0222]
[0223] 列出误差方程：
[0224]
[0225] 式中：
[0226] li＝-f(a0,b0,x00,y00,α0)，由残差确定单位权重误差，将残差大于三倍单位权重误差的点剔除，直至迭代求解收敛，所有点残差均小于三倍单位权重误差。
[0227] d5.采用至焦点距离和不变法椭圆拟合剔除粗差点并得到拟合参数；
[0228] 椭圆上任意点至两个焦点的距离和为常数，设两个焦点在平面坐标系内的坐标为(x1 y1)和(x2 y2)，距离和为c，列出误差方程：
[0229]
[0230] 参数(x1 y1)和(x2 y2)的初值取为测量点平面坐标均值减去和加上一个小量，即(xm+Δ y1+Δ)和(xm-Δ y1-Δ)，c的初值取为5.5，线性化，列出误差方程，并由残差确定单位权重误差，将残差大于三倍单位权重误差的点剔除，迭代求解收敛后，所有点残差均小于三倍单位权重误差：
[0231]
[0232] 式中：
[0233]
[0234]
[0235]
[0236] 椭圆中心的平面坐标系坐标：
[0237]
[0238] 椭圆长短半轴：
[0239] a＝c/2
[0240]
[0241] 椭圆长半轴方位角α：
[0242]
[0243] 平面坐标与椭圆坐标的转换关系：
[0244]
[0245] 即：
[0246]
[0247] 椭圆中心：
[0248]
[0249] e.处理器根据得到的椭圆拟合参数，对隧道进行变形分析并根据三维点云数据建立隧道平铺图，具体包括：
[0250] e1.对管片进行分块拟合，以检测管片的变形参数；
[0251] 由于每个管片为刚体，需要确定的变形量为每管片两端点相对于椭圆的位移，变形量由θ、D两个参数确定，θ为管片绕其中点P的旋转角，D为本管片沿OP方向的平移量。
[0252] 变形量θ、D的确定方法是经旋转平移后，该管片上所有观测点与椭圆最接近。
[0253] 参见图4：，A、B为某管片的端点，由管片方位角和椭圆心的平面坐标可计算A、B点的平面坐标，及其椭圆坐标，求出OA、OB在椭圆坐标系中的方位角α1、α2。考虑到椭圆构件为刚性材料，整体形变量较小，故将其近似视作圆形，则管片中点P可近似为方位角(α1+α2)/2的射线与椭圆的交点。
[0254] 以XeP表示P点的椭圆坐标系坐标，绕P点旋转θ，再沿OP方向平移D后，该管片上各测量点的椭圆坐标Xei变为X'ei：
[0255]
[0256] 其中αOP为P点在椭圆坐标系内的方位角。
[0257] 按旋转和平移后的坐标与拟合椭圆之间的偏差列出误差方程：
[0258] vi＝X'eiTΛX'ei-1
[0259] 线性化为：
[0260]
[0261] 其中各项导数如下：
[0262]
[0263]
[0264]
[0265]
[0266]
[0267] 在最小二乘 的条件下组成法方程，求解待估参数。对参数改正数进行迭代至收敛，即可解出旋转角θ和平移量D。
[0268] e2.检测管片端点的变形量：在椭圆坐标系内，求出管片的起点坐标、终点坐标及与椭圆的交点X'eS、X'eE。
[0269] 两交点转过旋转θ角、平移D后，交点坐标变为下式中的XeS、XeE：
[0270]
[0271]
[0272] 比较X'eS和XeS、X'eE和XeE，即得到该管片的变形量。
[0273] e3.检测点的变形量：对面坐标系方位角为α的点i，先取一个距原点10m的临时点，其平面坐标为：
[0274]
[0275] 将其转换至椭圆坐标：
[0276]
[0277] 在椭圆坐标系内，原点与临时点连线上点的坐标为(kxet kyet)T，代入椭圆方程，求得k:
[0278]
[0279] 从而得到i点的椭圆坐标X'ei，该点是Xei绕中点旋转平移得到的：
[0280]
[0281] i点形变量defi为：
[0282] defi＝|Xei|-|X'ei
[0283] e4.并根据三维点云数据建立隧道平铺图。
[0284] 图5为隧道扫描某一切片示意图，所述的隧道平铺图包括灰度图或反射率图，以隧道中轴线上一点o为原点，z轴指向隧道顶端，y轴为推扫方向，建立右手系，称o-xyz为中心坐标系。以轨道中线o′x′为横轴，隧道椭圆周方向o′y′将隧道面展开，建立隧道平面坐标系。使用一定大小BMP图像，根据每点坐标，计算像素位置，用点的灰度值或反射值填充,可依次获得隧道灰度图或反射率图。所述的灰度值，是指通过灰度显示材料的颜色，对于不同的材料，显示为不同的颜色。可通过灰度图，直接观看到隧道的病害情况。所述的反射值，是显示反射情况的数值，通过不同材料对于激光的反射情况不同，利用不同的反射情况，可知隧道病害。每一切片均可如此操作，再拼接所有切片，即可获得推扫里程区间灰度图或反射率图。
