CN104613966B - 一种地籍测量事后数据处理方法 - Google Patents
一种地籍测量事后数据处理方法 Download PDFInfo
- Publication number
- CN104613966B CN104613966B CN201510051138.9A CN201510051138A CN104613966B CN 104613966 B CN104613966 B CN 104613966B CN 201510051138 A CN201510051138 A CN 201510051138A CN 104613966 B CN104613966 B CN 104613966B
- Authority
- CN
- China
- Prior art keywords
- mrow
- data
- measurement unit
- msub
- inertial measurement
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Active
Links
Classifications
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01C—MEASURING DISTANCES, LEVELS OR BEARINGS; SURVEYING; NAVIGATION; GYROSCOPIC INSTRUMENTS; PHOTOGRAMMETRY OR VIDEOGRAMMETRY
- G01C21/00—Navigation; Navigational instruments not provided for in groups G01C1/00 - G01C19/00
- G01C21/20—Instruments for performing navigational calculations
Landscapes
- Engineering & Computer Science (AREA)
- Radar, Positioning & Navigation (AREA)
- Remote Sensing (AREA)
- Automation & Control Theory (AREA)
- Physics & Mathematics (AREA)
- General Physics & Mathematics (AREA)
- Navigation (AREA)
Abstract
本发明提出了一种地籍测量事后数据处理方法，依次通过对现场采集数据中的惯性测量单元数据进行数据预处理，去除误差；对经去除误差后的惯性测量单元数据进行分类，并对分类的数据分别采用卡尔曼固定区间最优平滑、纯捷联解算以及逆向导航解算进行处理，有效滤除了惯性测量单元数据中各种的有害误差成分，并实现地籍测量中待测点的高精度定位信息的获取。
Description
技术领域
本发明属于数据处理领域，涉及测量、测绘技术和捷联惯性定位定向技术，具体涉及一种地籍测量事后数据处理方法。
背景技术
我国经济社会的快速发展，土地利用类型、方式，以及权属变化大、范围广，要求地籍测量工作速度快、效率高、精度高。
随着GPS及全站仪等多技术结合的综合测量方法精度和效率的大幅提高，促进了地籍测量技术的进步，但GPS定位技术存在的环境依赖度高和全站仪存在的低效率等缺陷仍然无法满足地籍测量的需求。捷联惯性定位系统具有自主性强、抗干扰能力强、环境依赖度低及结构简单易维护等特点，近年来在地籍测量领域越来越受到研究人员的重视。
根据地籍测量低动态的测量特点和高精度的测量要求，将捷联惯性定位技术应用于地籍测量中。通过充分利用地籍测量过程“随停随测随走”的特点，以及可以在测量过程中人为引入辅助观测信息的优势，可以对捷联惯性定位系统定位误差进行控制。然而，由于受到惯性传感器精度限制，以及实际测量现场各种干扰因素的影响，通过在线解算，获得待测点的高精度定位信息，是十分困难的。
因此本发明提出了一种地籍测量事后数据处理方法，通过各种有效的事后数据处理手段，弥补惯性定位技术自身的不足，获取待测点的高精度定位信息。
发明内容
针对现有技术捷联惯性定位系统通过在线解算难以获得待测点的高精度定位信息的不足，本发提出了一种地籍测量事后数据处理方法，利用地籍测量允许进行事后数据处理的优势，引入数据预处理、卡尔曼固定区间最优平滑滤波、纯捷联解算和逆向导航解算相关技术，实现地籍测量中待测点的高精度定位信息的获取。
为实现上述目的，本发明采用如下技术方案：
本发明地籍测量事后数据处理方法，包括以下步骤：
步骤1、对现场采集数据中的惯性测量单元数据进行数据预处理，去除误差；
步骤2：对经步骤1去除误差后的惯性测量单元数据进行分类，将对应于坐标已知点、位置已知点、调平修正点、零速修正点以及待测点上静态测量过程时的惯性测量单元数据记为第一类；将对应于待测点上静态测量过程，将对应于所述五点两两之间的行进间动态测量过程时的惯性测量单元数据记为第二类；对应于待测点上静态测量过程，同时与其紧邻的上一个静态测量过程为坐标已知点上的静态测量过程时的惯性测量单元数据记为第三类；
步骤3：对经步骤2分类的数据分别作如下处理：第一类采用卡尔曼固定区间最优平滑对现场采集数据进行处理获得陀螺仪和加速度计的误差估计值，并根据陀螺仪和加速度的误差估计值分别对惯性测量单元数据中陀螺仪测得的惯性测量单元的角速度数据和加速度计测得的惯性测量单元的线加速度数据进行进一步误差补偿；第二类采用纯捷联解算对现场采集数据进行处理得到整个行进间动态测量过程惯性测量单元的位置、速度和姿态值；第三类采用逆向导航解算对坐标已知点与待测点之间行进间动态测量过程的现场采集数据进行处理获得坐标已知点和待测点之间的行进间动态测量过程惯性测量单元的位置值，并将所述逆向导航解算获得的位置值校正待测点处惯性测量单元的位置值。
进一步地，所述步骤1中数据预处理为采用时序建模与卡尔曼滤波相结合方法去除惯性测量单元数据误差；所述的时序建模包括野值剔除、低通滤波、数据统计分析、模型结构选取、模型参数估计和建立误差模型；所述的卡尔曼滤波根据上述时序建模建立的误差模型，对陀螺仪和加速度计的误差进行估计和补偿。
进一步地，所述第一类采用卡尔曼固定区间最优平滑对现场采集数据进行处理，具体为：
1)建立与静态测量过程相对应的卡尔曼滤波观测模型；
2)利用卡尔曼滤波观测模型对现场采集数据进行正向卡尔曼滤波；
3)对正向卡尔曼滤波的结果进行逆向平滑滤波；
4)获得陀螺仪和加速度计的误差估计值；
5)在滤波完成后，根据陀螺仪和加速度的误差估计值分别对惯性测量单元数据中陀螺仪测得的惯性测量单元的角速度数据和加速度计测得的惯性测量单元的线加速度数据进行进一步误差补偿。
进一步地，第二类采用纯捷联解算对现场采集数据进行处理，选择当地地理坐标系作为导航坐标系，选择四元数法进行姿态更新，具体为：
1)进行惯性测量单元的速度、位置和姿态进行初始化，并根据姿态初始值计算姿态矩阵初始值和四元数初始值；
2)按时间顺序读取现场采集数据中的一组惯性测量单元数据；
3)利用姿态矩阵对加速度计测得的下线加速度数据进行坐标变换；
4)根据比力方程，利用坐标变换的得到的线速度数据，进行速度更新，得到惯性测量单元速度最新值；
5)利用惯性测量单元速度最新值进行位置更新，得到惯性测量单元位置最新值；
6)利用惯性测量单元速度最新值和位置最新值进行地球角速度和位置角速度的计算，具体如下式：
其中：ωie为地球自转角速度；L为当地地理纬度；VE、VN分别为惯性测量单元东向、北向速度分量；RN、RM分别为当地卯酉圈曲率半径和子午圈曲率半径；h为当地高度分量；
7)利用6)中计算得到的地球角速度、位置角速度以及姿态矩阵和陀螺仪测得的角速度数据进行姿态角速度的计算，具体如下式：
其中：为惯性定位系统中陀螺仪直接输出的角速度值；为姿态矩阵；为地球角速度；为位置角速度；
8)利用7)中计算得到的姿态角速度进行四元数的更新，得到四元数最新值；
9)利用8)中计算得到的四元数最新值进行姿态矩阵更新，得到姿态矩阵最新值；
10)根据姿态与姿态矩阵的对应关系，计算得到惯性测量单元的姿态最新值；
11)重复步骤2)～10)，用所得到的位置、速度、姿态、姿态矩阵、四元数最新值进行递推计算，直到所有现场采集数据处理完毕，得到整个行进间动态测量过程惯性测量单元的位置、速度和姿态值。
进一步地，第三类采用逆向导航解算对坐标已知点与待测点之间行进间动态测量过程的现场采集数据进行处理，具体为：
1)将纯捷联解算的位置终值作为逆向导航解算的位置初始值，将纯解算的姿态终值作为逆向导航解算的姿态初始值，将纯捷联解算的速度终值取反，作为逆向导航解算速度初始值；
2)逆向读取现场采集数据中惯性测量单元数据，将陀螺仪测得的角速度值和地球自转角速度值符号取反，进行逆向导航解算，得到坐标已知点和待测点之间的行进间动态测量过程惯性测量单元的位置、速度和姿态值；
3)根据上述逆向导航解算获得的位置值，校正待测点处惯性测量单元的位置值。
与现有技术相比，本发明的有益效果是：
1)通过对惯性测量单元数据进行数据预处理，有效滤除了惯性测量单元数据中各种的有害误差成分，提高了陀螺仪和加速度计的使用精度。
2)通过卡尔曼固定区间最优平滑数据处理环节，通过增加反向的平滑滤波环节，弥补了单一采用卡尔曼滤波数据利用率低，滤波估计结果平稳性和精确性不高的缺陷。
3)通过纯捷联解算数据处理环节对五种点两两之间的行进间动态测量过程的定位信息进行更新，保证了地籍测量整个过程结果的完整性，并为进行逆向导航解算提供基础数据。
4)通过逆向导航解算数据处理环节，充分利用坐标已知点处卡尔曼固定区间最优平滑处理后的结果，基于坐标已知点和待测点之间的数据正向纯捷联解算的结果对惯性测量单元数据进行逆向导航解算，将逆向导航解算获得的位置值校正待测点处惯性测量单元的位置误差，有效提高了待测点处的位置测量精度。
5)地籍测量事后数据处理方法总体上提高了地籍测量的定位精度，有效解决了在线解算不能给出待测点高精度定位信息的问题。
附图说明：
图1为地籍测量事后数据处理总体流程图；
图2为数据预处理流程图；
图3为卡尔曼固定区间最优平滑流程图；
图4为纯捷联解算流程图；
图5为逆向导航解算流程图。
具体实施方式：
下面结合附图，对本发明的具体实施方式做进一步说明。
如图1-5所示，本发明地籍测量事后数据处理方法，包括以下几个步骤：
步骤1、对现场采集数据中的惯性测量单元数据进行数据预处理，去除误差。
步骤2：对经步骤1去除误差后的惯性测量单元数据进行分类，将对应于坐标已知点、位置已知点、调平修正点、零速修正点以及待测点上静态测量过程时的惯性测量单元数据记为第一类；将对应于待测点上静态测量过程，将对应于所述五点两两之间的行进间动态测量过程时的惯性测量单元数据记为第二类；对应于待测点上静态测量过程，同时与其紧邻的上一个静态测量过程为坐标已知点上的静态测量过程时的惯性测量单元数据记为第三类。
步骤3：对经步骤2分类的数据分别作如下处理：第一类采用卡尔曼固定区间最优平滑对现场采集数据进行处理获得陀螺仪和加速度计的误差估计值，并根据陀螺仪和加速度的误差估计值分别对惯性测量单元数据中陀螺仪测得的惯性测量单元的角速度数据和加速度计测得的惯性测量单元的线加速度数据进行进一步误差补偿；第二类采用纯捷联解算对现场采集数据进行处理得到整个行进间动态测量过程惯性测量单元的位置、速度和姿态值；第三类采用逆向导航解算对坐标已知点与待测点之间行进间动态测量过程的现场采集数据进行处理获得坐标已知点和待测点之间的行进间动态测量过程惯性测量单元的位置值，并将所述逆向导航解算获得的位置值校正待测点处惯性测量单元的位置值。
现场采集数据包括惯性测量单元数据以及辅助信息，所述的辅助信息包括地籍测量过程中外界提供的惯性测量单元的位置、速度、姿态观测信息，其实际内容因其对应的惯性测量单元数据不同而不同，具体为：
当现场采集数据对应于坐标已知点上的静态测量过程时，其辅助信息为坐标已知点提供的惯性测量单元位置、速度和姿态观测信息；
当现场采集数据对应于位置已知点上静态测量过程时，其辅助信息为位置已知点提供的惯性测量单元位置、速度观测信息，此时姿态数据无效；
当现场采集数据对应于调平修正点上静态测量过程时，其辅助信息为调平修正点提供的惯性测量单元速度和姿态观测信息，辅助信息的位置数据无效；
当现场采集数据对应于零速修正点上静态测量过程时，其辅助信息为零速修正点提供的惯性测量单元速度观测信息，辅助信息的位置、姿态数据无效；
当现场采集数据对应于待测点上静态测量过程时，其辅助信息为待测点提供的惯性测量单元速度观测信息，辅助信息的位置、姿态数据无效；
当现场采集数据对应于所述五点两两之间的行进间动态测量过程时，辅助信息的位置、速度、姿态数据均无效。
惯性测量单元经过测试和标定后可补偿陀螺仪和加速度计的确定性误差，如标定因数误差、安装误差等，地籍测量中陀螺仪的随机漂移误差和加速度计的随机偏置误差成为影响系统精度的主要因素。所述步骤1中数据预处理采用时序建模与卡尔曼滤波相结合方法去除惯性测量单元数据误差；所述的时序建模包括野值剔除、低通滤波、数据统计分析、模型结构选取、模型参数估计和建立误差模型；所述的卡尔曼滤波根据上述时序建模建立的误差模型，对陀螺仪和加速度计的误差进行估计和补偿。
如图2所示，以陀螺仪的数据预处理为例，数据预处理具体实施过程为：
1)、野值剔除：陀螺仪进行长时间的数据采集时，不可避免的会存在外界环境如振动或者电路电压不稳定等干扰对陀螺仪测试输出数据产生影响，在数据值上体现为异于常值，称为野值。主要表现为数据图中有尖峰出现或明显豁口出现。数据预处理，首先剔除这些野值，避免对后续的数据建模产生不利影响。
采用基于拉依达准则的奇异数据滤波法剔除测量数据中的野值。拉依达准则的具体描述是：当测量次数N足够多且测量服从正态分布时，在各次测量值中，若某次测量值Xi所对应的剩余误差Vi>3σ，则认为该Xi为野值，予以剔除。
其中：剩余误差定义为为N次测量值的均值；σ为N次测量值的标准差。
2)、低通滤波：在地籍测量系统中，由于惯性测量单元的振幅摇摆频率不会大于10Hz，因此，高于此频率的信号均可视为干扰噪声去除，这对消除实际测量时的车辆经过、施工等振动干扰十分有用。剔除野之后，对陀螺仪输出的数字信号进行低通滤波，进而提高信噪比。
陀螺仪输出信号依照其工作环境，所测量的角速度不可能发生突变。故其滤波输出不仅与过去一段时刻的输入有关，与过去一段时间的滤波输出亦有关系，满足无限脉冲响应(IIR)滤波器条件。综合考虑，采用巴特沃斯低通滤波器，设定滤波器的参数为：通带截止频率10Hz，阻带截止频率为15Hz，通带、阻带波动值分别为1dB和30dB。
3)、数据统计分析：实际测得的陀螺仪随机漂移序列一般为非平稳随机过程，在剔除野值并进行低通滤波后，需要先对序列趋向性检验，具体为：对序列利用Daniel检验法检验序列是否为平稳序列，如果检验为非平稳序列，则需要进行平稳化处理，处理方法是先利用最小二乘法低阶拟合序列的趋势项，去除趋势项得到残差，利用Daniel检验法继续检验残差是否为平稳序列，如果检验仍为非平稳序列，则提高最小二乘法的阶数继续拟合序列的趋势项，直到残差序列为平稳序列，完成序列的趋向性检验。
陀螺仪随机漂移经过趋向性检验并提取线性趋势项后，漂移数据中仍可能含有周期项，须进一步识别并提取隐含周期项。采用周期图分析和Fisher检验方法对去掉趋势项的陀螺仪漂移信号的隐含周期项进行识别。
4)、模型结构选取和模型参数估计：实际应用中，陀螺仪随机漂移通常采用ARMA和AR模型逼近随机序列，由于陀螺仪随机漂移模型阶次都比较低，一般不超过2到3阶，因此可以在模型参数数目等于3的范围内寻找合适的模型。并且对于最小实现的线性系统，ARMA模型的自回归阶数不低于滑动平均部分的阶次，故实际需要选择的只有AR(1)、AR(2)、AR(3)、ARMA(1,1)和ARMA(2,1)共5种模型。对于AR模型和ARMA模型的参数估计，可调用Matlab工具中相应的时序建模函数完成。
5)、建立误差模型：根据步骤4确定的模型结构和模型参数，确定陀螺仪随机漂移误差模型。
6)、传感器误差的估计与补偿：建立陀螺仪的随机漂移误差模型后，采用卡尔曼滤波进行滤波进行误差估计。以AR(3)模型为例，该模型可以表示为：
其中：n表示n时刻；x(n)表示n时刻系统状态；分别为x(n-1)、x(n-2)、x(n-3)对应的自回归系数；{ξ(n)}为零均值白噪声序列。
将上式改写为状态空间模型为：
Yk＝[0 0 1]Xk+νk＝HkXk+νk
其中：k表示k时刻；Φk表示k时刻状态一步转移矩阵；Γk为k时刻过程噪声输入阵；ξk为k时刻过程噪声；Yk为k时刻滤波观测数据；Hk为k时刻的量测矩阵；νk为k时刻测量噪声。
滤波初值设置为X0、P0，设{ξk}与{νk}的方差分别为Qk和Rk，则可构造如下的卡尔曼滤波器：
其中：k表示k时刻；Φk表示k时刻状态一步转移矩阵；Γk为k时刻过程噪声输入阵；Hk+1为k+1时刻的量测矩阵；为k+1时刻的状态一步预测；为k+1时刻的状态估计；Pk+1,k为k+1时刻一部预测误差方差阵；Gk+1为k+1时刻滤波增益矩阵；Pk+1,k+1为k+1时刻估计误差方差阵。
估计出陀螺仪的随机误差后，对其进行补偿，减小陀螺仪的漂移，提高陀螺仪的使用精度。
加速度计的零偏建模方法和补偿过程和陀螺仪相同。
如图3所示，第一类采用卡尔曼固定区间最优平滑对现场采集数据进行处理，具体实施过程为：
1)、建立与静态测量子过程相对应的卡尔曼滤波观测模型，具体方法描述如下：
在建立捷联惯性定位系统误差模型过程中，将系统状态向量X取为：
X＝[δVe δVn φe φn φu δL δλ ▽x ▽y εx εy εz]T
其中：δVe和δVn分别是东向和北向速度误差；φe、φn和φu分别为两个水平失准角和方位失准角；δL和δλ分别是纬度和经度误差；▽x、▽y和▽z分别为东向、北向加速度计偏置；εx、εy和εz分别为东向、北向及方位陀螺仪漂移。
定义捷联姿态矩阵为：
其中：C11、C12、C13、C21、C22、C23、C31、C32、C33为姿态矩阵中对应元素；H为航向角；P为纵摇角；R为横摇角。
a、当信息源点为坐标已知点，外部可以提供精确的位置、姿态和水平速度信息，观测方程如下：
其中：δH、δP、δR分别为航向角、俯仰角、横滚角姿态观测量；δVe、δVn为东向和北向速度观测量；δL、δλ分别为纬度、经度位置观测量；δH、δP、δR构成观测矢量Z1，H1为对应于位置观测量的观测矩阵，ν为对应观测噪声矢量；δVe、δVn构成观测矢量Z2，H2为对应于位置观测量的观测矩阵，κ为对应观测噪声矢量；δL、δλ构成观测矢量Z3，H3为对应于位置观测量的观测矩阵，ζ为对应观测噪声矢量；Z1、Z2、Z3构成坐标已知点处观测矢量Zl；X为系统状态矢量。
令
其中：HSINS、PSINS和RSINS为捷联惯性定位系统解算得到的姿态角；Hplatform、Pplatform和Rplatform为地籍测量辅助设备给出的参考姿态角；δH、δP、δR分别为航向角、俯仰角、横滚角姿态观测量，三者构成坐标已知点处姿态观测矢量Z1，ν为对应观测噪声矢量；X为系统状态矢量；H1为姿态观测量对应的观测矩阵，其具体形式如下：
其中：C11、C12、C13、C21、C22、C23、C31、C32、C33为姿态矩阵Cb n中对应元素。
令
其中：为捷联惯性定位系统解算得到的水平速度；δVe、δVn为东向和北向速度观测量，二者构成坐标已知点处速度观测矢量Z1，κ为对应观测噪声矢量；X为系统状态矢量；H2为速度观测量对应的观测矩阵，其具体形式如下：
令
其中：Lsins、λsins为捷联惯性定位系统解算得到的纬度、经度值；LGPS、λGPS是坐标已知点处给出的精确纬度和经度值；δL、δλ构成观测矢量Z3，ζ为对应观测噪声矢量；X为系统状态矢量；H3为位置观测量对应的观测矩阵，其具体形式为：
b、当信息源点为调平修正点，该点处无法得到精确的位置参考信息，所以相对坐标已知点，只有三个姿态角和两个水平速度五个观测量，观测方程如下：
其中：δH、δP、δR分别为航向角、俯仰角、横滚角姿态观测量；δVe、δVn为东向和北向速度观测量；δH、δP、δR、δL、δλ构成调平修正点处观测矢量Zf；X为系统状态矢量；ν为对应于姿态观测量的观测噪声矢量；κ为对应于速度观测量的观测噪声矢量；H1为对应于姿态观测量的观测矩阵分量，H2为对应于速度观测量的观测矩阵分量，二者具体形式如下：
其中：C11、C12、C13、C21、C22、C23、C31、C32、C33为姿态矩阵中对应元素。
c、当信息源点为位置已知点，以此点的水平速度信息和纬度、经度信息作为外观测量，观测方程如下：
其中：δVe、δVn为东向和北向速度观测量；δL、δλ分别为纬度、经度位置观测量；δVe、δVn、δL、δλ构成位置已知点处观测矢量Zλ；X为系统状态矢量；κ为对应于速度观测量的观测噪声矢量；ζ为对应于位置观测量的观测噪声矢量；H2为对应于速度观测量的观测矩阵分量，H3为对应于位置观测量的观测矩阵分量，二者具体形式如下：
d、当信息源点为零速修正点，利用两个水平轴向的速度误差作为卡尔曼滤波器外观测量，得到如下观测方程：
其中：为捷联解算得到的东向和北向速度；δVe、δVn为东向和北向速度观测量，二者构成零速修正点处观测矢量Zv；X为系统状态矢量；κ为对应于速度观测量的观测噪声矢量；H2为观测矩阵，其具体形式如下：
2)、利用卡尔曼滤波观测模型对现场采集数据进行正向卡尔曼滤波，具体方法为：
将滤波初值设置为X0、P0，设{ξk}与{νk}的方差分别为Qk和Rk，构造如下的卡尔曼滤波器：
其中：k表示k时刻；Φk表示k时刻状态一步转移矩阵；Γk为k时刻过程噪声输入阵；Hk+1为k+1时刻的量测矩阵；为k+1时刻的状态一步预测；为k+1时刻的状态估计；Pk+1,k为k+1时刻一部预测误差方差阵；Gk+1为k+1时刻滤波增益矩阵；Pk+1,k+1为k+1时刻估计误差方差阵。
利用卡尔曼滤波估计误差状态Xk，并存储各时刻的有关滤波信息，供逆向平滑滤波算法使用。
3)、利用正向卡尔曼滤波的结果，进行逆向平滑滤波，逆向平滑基本方程为：
其边界条件为：当k＝N-1时，Pk+1/N＝PN/N。
其中：k＝N-1,N-2,...,1,0；和Pk/N分别为tk时刻状态的平滑值和平滑误差方差阵；和Pk+1/k分别是tk时刻的滤波估计值、滤波估计误差均方差阵和一步预测误差均方差阵；Aa(k)为平滑增益矩阵，它与由滤波器确定的估计误差均方阵Pk/k和一步预测误差协方差阵Pk+1/k有关，而与平滑器的平滑误差协方差阵Pk/N无关。固定区间平滑算法从k＝N开始，计算到k＝0，由反向逐步递推，依次得到
4)、获得陀螺仪和加速度计的误差估计值；
5)、在滤波完成后，根据陀螺仪和加速度的误差估计值分别对惯性测量单元数据中陀螺仪测得的惯性测量单元的角速度数据和加速度计测得的惯性测量单元的线加速度数据进行进一步误差补偿。
如图4所示，地籍测量中，选择当地地理坐标系作为导航坐标系，选择四元数法进行姿态更新，第二类采用纯捷联解算对现场采集数据进行处理具体实施过程为：
1)、进行惯性测量单元的速度、位置和姿态进行初始化，并根据姿态初始值与姿态矩阵初始值和四元数初始值的关系计算姿态矩阵初始值和四元数初始值；
2)、按时间顺序读取现场采集数据中的一组惯性测量单元数据；
3)、利用姿态矩阵对加速度计测得的下线加速度数据进行坐标变换，具体如下式：
将加速度计测量的载体系上的比力fb转换为导航坐标系上的比力fn
其中，fb为加速度计测量的线加速度在载体坐标系中表示，fn为加速度计测量的线加速度在导航坐标系中表示，为载体坐标系到导航坐标系的姿态变换矩阵。
4)、根据比力方程，利用坐标变换的得到的线速度数据，进行速度更新，得到惯性测量单元速度最新值，具体如下式：
其中：Vn为载体三维速度矢量；为载体坐标系到导航坐标系的姿态变换矩阵；fb为三轴加速度计比力数据构成的矢量；为地球速度；为位置速度；gn为载体所在位置重力矢量在导航坐标系投影。
通过解该微分方程来求得惯性测量单元速度最新值。
5)、利用惯性测量单元速度最新值进行位置更新，得到惯性测量单元位置最新值，具体如下式：
其中：λ为当地地理经度；L为当地地理纬度；VE、VN分别为载体东向、北向速度分量；RN、RM分别为当地卯酉圈曲率半径和子午圈曲率半径；h当地高度分量。
过解该微分方程来求得惯性测量单元位置最新值。
6)、利用惯性测量单元速度最新值和位置最新值进行地球角速度和位置角速度的计算，具体如下式：
其中：ωie为地球自转角速度；L为当地地理纬度；VE、VN分别为惯性测量单元东向、北向速度分量；RN、RM分别为当地卯酉圈曲率半径和子午圈曲率半径；h为当地高度分量；
7)、利用6)中计算得到的地球角速度、位置角速度以及姿态矩阵和陀螺仪测得的角速度数据进行姿态角速度的计算，具体如下式：
其中：为惯性定位系统中陀螺仪直接输出的角速度值；为姿态矩阵；为地球角速度；为位置角速度；
8)、利用姿态角速度进行四元数的更新，得到四元数最新值，具体如下式：
设载体坐标系相对导航坐标系的转动四元数为：
Q＝q0+q1i+q2j+q3k
其中：q0、q1、q2、q3为四个元素；i、j、k为导航坐标系的基。
Q的即时修正通过求解下面的四元数微分方程来实现：
其中：分别为姿态速度在载体坐标系上x、y、z三个轴向上的投影。
为确保“数学平台”中导航坐标系的3个基的模均为1，需要对四元数周期性的进行归一化处理，以欧几里得范数最小为准则的四元数最佳归一化由下式获得：
其中：Q为归一化前后的四元数值；i、j、k为导航坐标系的基； 为四元数的四个分量；为四元数Q的四个分量，进行若干次即时修正后进行归一化处理。
9)、利用四元数最新值进行姿态矩阵更新，得到姿态矩阵最新值；
其中：其中姿态矩阵，q0、q1、q2、q3为四元数的四个分量。
10)、根据姿态角与姿态矩阵的关系，计算得到惯性测量单元的姿态最新值；
11)、重复步骤2～10，用所得到的位置、速度、姿态、姿态矩阵、四元数最新值进行递推计算，直到所有现场采集数据处理完毕，得到整个行进间动态测量过程惯性测量单元的位置、速度和姿态值。
如图5所示，第三类采用逆向导航解算对坐标已知点与待测点之间行进间动态测量过程的现场采集数据进行处理，具体实施过程为：
1)、将纯捷联解算的位置终值作为逆向导航解算的位置初始值，将纯捷联解算的姿态终值作为逆向导航解算的姿态初始值，将纯捷联解算的速度终值取反，作为逆向导航解算速度初始值。
2)、逆向读取现场采集数据中惯性测量单元数据，将陀螺仪测得的角速度值和地球自转角速度值符号取反，进行逆向导航解算，得到坐标已知点和待测点之间的行进间动态测量过程惯性测量单元的位置、速度和姿态值。
3)、根据上述逆向导航解算获得的位置值，校正待测点处惯性测量单元的位置值。
将所有解算结果保存到.txt或.dat格式的数据文件，供地籍测量后续分析使用。
地籍测量的最终目的在于给出待测点的高精度定位信息。本发明，通过引入数据预处理、卡尔曼固定区间最优平滑、纯捷联解算、逆向导航解算四个核心环节对地籍测量的现场采集数据进行了事后数据处理，最终获得了待测点的高精度定位信息。通过本发明提出的地籍测量事后数据处理方法，弥补惯性定位技术自身的不足，有效解决了在线解算不能给出待测点的高精度定位信息的问题。
Claims (2)
1.一种地籍测量事后数据处理方法，其特征在于，包括以下步骤：
步骤1、对现场采集数据中的惯性测量单元数据进行数据预处理，去除误差；
步骤2：对经步骤1去除误差后的惯性测量单元数据进行分类，将对应于坐标已知点、位置已知点、调平修正点、零速修正点以及待测点上静态测量过程时的惯性测量单元数据记为第一类；将对应于待测点上静态测量过程，将对应于所述五点两两之间的行进间动态测量过程时的惯性测量单元数据记为第二类；对应于待测点上静态测量过程，同时与其紧邻的上一个静态测量过程为坐标已知点上的静态测量过程时的惯性测量单元数据记为第三类；
步骤3：对经步骤2分类的数据分别作如下处理：第一类采用卡尔曼固定区间最优平滑对现场采集数据进行处理获得陀螺仪和加速度计的误差估计值，并根据陀螺仪和加速度的误差估计值分别对惯性测量单元数据中陀螺仪测得的惯性测量单元的角速度数据和加速度计测得的惯性测量单元的线加速度数据进行进一步误差补偿；第二类采用纯捷联解算对现场采集数据进行处理得到整个行进间动态测量过程惯性测量单元的位置、速度和姿态值；第三类采用逆向导航解算对坐标已知点与待测点之间行进间动态测量过程的现场采集数据进行处理获得坐标已知点和待测点之间的行进间动态测量过程惯性测量单元的位置值，并将所述逆向导航解算获得的位置值校正待测点处惯性测量单元的位置值；
所述第一类采用卡尔曼固定区间最优平滑对现场采集数据进行处理，具体为：
1)建立与静态测量过程相对应的卡尔曼滤波观测模型；
2)利用卡尔曼滤波观测模型对现场采集数据进行正向卡尔曼滤波；
3)对正向卡尔曼滤波的结果进行逆向平滑滤波；
4)获得陀螺仪和加速度计的误差估计值；
5)在滤波完成后，根据陀螺仪和加速度的误差估计值分别对惯性测量单元数据中陀螺仪测得的惯性测量单元的角速度数据和加速度计测得的惯性测量单元的线加速度数据进行进一步误差补偿；
第二类采用纯捷联解算对现场采集数据进行处理，选择当地地理坐标系作为导航坐标系，选择四元数法进行姿态更新，具体为：
1)进行惯性测量单元的速度、位置和姿态进行初始化，并根据姿态初始值计算姿态矩阵初始值和四元数初始值；
2)按时间顺序读取现场采集数据中的一组惯性测量单元数据；
3)利用姿态矩阵对加速度计测得的下线加速度数据进行坐标变换；
4)根据比力方程，利用坐标变换的得到的线速度数据，进行速度更新，得到惯性测量单元速度最新值；
5)利用惯性测量单元速度最新值进行位置更新，得到惯性测量单元位置最新值；
6)利用惯性测量单元速度最新值和位置最新值进行地球角速度和位置角速度的计算，具体如下式：
<mrow>
<msubsup>
<mi>&omega;</mi>
<mrow>
<mi>i</mi>
<mi>e</mi>
</mrow>
<mi>n</mi>
</msubsup>
<mo>=</mo>
<mfenced open = "[" close = "]">
<mtable>
<mtr>
<mtd>
<mn>0</mn>
</mtd>
</mtr>
<mtr>
<mtd>
<msub>
<mi>&omega;</mi>
<mrow>
<mi>i</mi>
<mi>e</mi>
</mrow>
</msub>
<mi>cos</mi>
<mi>L</mi>
</mtd>
</mtr>
<mtr>
<mtd>
<msub>
<mi>&omega;</mi>
<mrow>
<mi>i</mi>
<mi>e</mi>
</mrow>
</msub>
<mi>sin</mi>
<mi>L</mi>
</mtd>
</mtr>
</mtable>
</mfenced>
</mrow>
1
<mrow>
<msubsup>
<mi>&omega;</mi>
<mrow>
<mi>e</mi>
<mi>n</mi>
</mrow>
<mi>n</mi>
</msubsup>
<mo>=</mo>
<mfenced open = "[" close = "]">
<mtable>
<mtr>
<mtd>
<mfrac>
<mrow>
<mo>-</mo>
<msub>
<mi>V</mi>
<mi>N</mi>
</msub>
</mrow>
<mrow>
<msub>
<mi>R</mi>
<mi>M</mi>
</msub>
<mo>+</mo>
<mi>h</mi>
</mrow>
</mfrac>
</mtd>
</mtr>
<mtr>
<mtd>
<mfrac>
<msub>
<mi>V</mi>
<mi>E</mi>
</msub>
<mrow>
<msub>
<mi>R</mi>
<mi>N</mi>
</msub>
<mo>+</mo>
<mi>h</mi>
</mrow>
</mfrac>
</mtd>
</mtr>
<mtr>
<mtd>
<mfrac>
<msub>
<mi>V</mi>
<mi>E</mi>
</msub>
<mrow>
<msub>
<mi>R</mi>
<mi>N</mi>
</msub>
<mo>+</mo>
<mi>h</mi>
</mrow>
</mfrac>
<mi>tan</mi>
<mi>L</mi>
</mtd>
</mtr>
</mtable>
</mfenced>
</mrow>
<mrow>
<msubsup>
<mi>&omega;</mi>
<mrow>
<mi>i</mi>
<mi>e</mi>
</mrow>
<mi>n</mi>
</msubsup>
<mo>+</mo>
<msubsup>
<mi>&omega;</mi>
<mrow>
<mi>e</mi>
<mi>n</mi>
</mrow>
<mi>n</mi>
</msubsup>
<mo>=</mo>
<mfenced open = "[" close = "]">
<mtable>
<mtr>
<mtd>
<mfrac>
<mrow>
<mo>-</mo>
<msub>
<mi>V</mi>
<mi>N</mi>
</msub>
</mrow>
<mrow>
<msub>
<mi>R</mi>
<mi>M</mi>
</msub>
<mo>+</mo>
<mi>h</mi>
</mrow>
</mfrac>
</mtd>
</mtr>
<mtr>
<mtd>
<msub>
<mi>&omega;</mi>
<mrow>
<mi>i</mi>
<mi>e</mi>
</mrow>
</msub>
<mi>cos</mi>
<mi>L</mi>
<mo>+</mo>
<mfrac>
<msub>
<mi>V</mi>
<mi>E</mi>
</msub>
<mrow>
<msub>
<mi>R</mi>
<mi>N</mi>
</msub>
<mo>+</mo>
<mi>h</mi>
</mrow>
</mfrac>
</mtd>
</mtr>
<mtr>
<mtd>
<msub>
<mi>&omega;</mi>
<mrow>
<mi>i</mi>
<mi>e</mi>
</mrow>
</msub>
<mi>sin</mi>
<mi>L</mi>
<mo>+</mo>
<mfrac>
<msub>
<mi>V</mi>
<mi>E</mi>
</msub>
<mrow>
<msub>
<mi>R</mi>
<mi>N</mi>
</msub>
<mo>+</mo>
<mi>h</mi>
</mrow>
</mfrac>
<mi>tan</mi>
<mi>L</mi>
</mtd>
</mtr>
</mtable>
</mfenced>
</mrow>
其中：ωie为地球自转角速度；L为当地地理纬度；VE、VN分别为惯性测量单元东向、北向速度分量；RN、RM分别为当地卯酉圈曲率半径和子午圈曲率半径；h为当地高度分量；
7)利用6)中计算得到的地球角速度、位置角速度以及姿态矩阵和陀螺仪测得的角速度数据进行姿态角速度的计算，具体如下式：
<mrow>
<msubsup>
<mi>&omega;</mi>
<mrow>
<mi>n</mi>
<mi>b</mi>
</mrow>
<mi>b</mi>
</msubsup>
<mo>=</mo>
<msubsup>
<mi>&omega;</mi>
<mrow>
<mi>i</mi>
<mi>b</mi>
</mrow>
<mi>b</mi>
</msubsup>
<mo>-</mo>
<msubsup>
<mi>C</mi>
<mi>n</mi>
<mi>b</mi>
</msubsup>
<mrow>
<mo>(</mo>
<msubsup>
<mi>&omega;</mi>
<mrow>
<mi>i</mi>
<mi>e</mi>
</mrow>
<mi>n</mi>
</msubsup>
<mo>+</mo>
<msubsup>
<mi>&omega;</mi>
<mrow>
<mi>e</mi>
<mi>n</mi>
</mrow>
<mi>n</mi>
</msubsup>
<mo>)</mo>
</mrow>
</mrow>
其中：为惯性定位系统中陀螺仪直接输出的角速度值；为姿态矩阵；为地球角速度；为位置角速度；
8)利用7)中计算得到的姿态角速度进行四元数的更新，得到四元数最新值；
9)利用8)中计算得到的四元数最新值进行姿态矩阵更新，得到姿态矩阵最新值；
10)根据姿态与姿态矩阵的对应关系，计算得到惯性测量单元的姿态最新值；
11)重复步骤2)～10)，用所得到的位置、速度、姿态、姿态矩阵、四元数最新值进行递推计算，直到所有现场采集数据处理完毕，得到整个行进间动态测量过程惯性测量单元的位置、速度和姿态值；
第三类采用逆向导航解算对坐标已知点与待测点之间行进间动态测量过程的现场采集数据进行处理，具体为：
1)将纯捷联解算的位置终值作为逆向导航解算的位置初始值，将纯解算的姿态终值作为逆向导航解算的姿态初始值，将纯捷联解算的速度终值取反，作为逆向导航解算速度初始值；
2)逆向读取现场采集数据中惯性测量单元数据，将陀螺仪测得的角速度值和地球自转角速度值符号取反，进行逆向导航解算，得到坐标已知点和待测点之间的行进间动态测量过程惯性测量单元的位置、速度和姿态值；
3)根据上述逆向导航解算获得的位置值，校正待测点处惯性测量单元的位置值。
2.根据权利要求1所述的地籍测量事后数据处理方法，其特征在于，所述步骤1中数据预处理为采用时序建模与卡尔曼滤波相结合方法去除惯性测量单元数据误差；所述的时序建模包括野值剔除、低通滤波、数据统计分析、模型结构选取、模型参数估计和建立误差模型；所述的卡尔曼滤波根据上述时序建模建立的误差模型，对陀螺仪和加速度计的误差进行估计和补偿。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201510051138.9A CN104613966B (zh) | 2015-01-30 | 2015-01-30 | 一种地籍测量事后数据处理方法 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201510051138.9A CN104613966B (zh) | 2015-01-30 | 2015-01-30 | 一种地籍测量事后数据处理方法 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN104613966A CN104613966A (zh) | 2015-05-13 |
| CN104613966B true CN104613966B (zh) | 2017-11-14 |
Family
ID=53148525
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN201510051138.9A Active CN104613966B (zh) | 2015-01-30 | 2015-01-30 | 一种地籍测量事后数据处理方法 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN104613966B (zh) |
Families Citing this family (6)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN105607105A (zh) * | 2016-03-07 | 2016-05-25 | 东南大学 | 一种面向不动产实地测量的惯性定位方法 |
| CN109211234B (zh) * | 2017-06-30 | 2022-03-29 | 湖南格纳微信息科技有限公司 | 一种矿下惯性测绘方法及装置 |
| CN109725213B (zh) * | 2018-12-12 | 2020-05-22 | 江南大学 | 基于逆向卡尔曼滤波器的Buck变换器故障检测方法 |
| CN110703288A (zh) * | 2019-10-15 | 2020-01-17 | 河海大学 | 一种适用于天线内置卫星定位设备的信号质量评估方法 |
| CN110730502B (zh) * | 2019-10-23 | 2020-11-03 | 珠海优特电力科技股份有限公司 | 一种定位方法及装置 |
| CN111193495B (zh) * | 2019-12-12 | 2023-06-02 | 浙江工业大学 | 一种工件校直数据的滤波处理方法 |
Citations (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN102033235A (zh) * | 2010-10-15 | 2011-04-27 | 东南大学 | 基于gps/sins城镇地籍图根点快速测量及事后数据处理方法 |
Family Cites Families (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US7541974B2 (en) * | 2005-12-15 | 2009-06-02 | Trimble Navigation Limited | Managed traverse system and method to acquire accurate survey data in absence of precise GPS data |
-
2015
- 2015-01-30 CN CN201510051138.9A patent/CN104613966B/zh active Active
Patent Citations (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN102033235A (zh) * | 2010-10-15 | 2011-04-27 | 东南大学 | 基于gps/sins城镇地籍图根点快速测量及事后数据处理方法 |
Non-Patent Citations (1)
| Title |
|---|
| 基于地籍测量的惯性系统初始对准方法研究;樊征等;《测绘与空间地理信息》;20140930;第37卷(第9期);46-49 * |
Also Published As
| Publication number | Publication date |
|---|---|
| CN104613966A (zh) | 2015-05-13 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN104613966B (zh) | 一种地籍测量事后数据处理方法 | |
| CN105509739B (zh) | 采用固定区间crts平滑的ins/uwb紧组合导航系统及方法 | |
| CN105737823B (zh) | 一种基于五阶ckf的gps/sins/cns组合导航方法 | |
| CN102706366B (zh) | 一种基于地球自转角速率约束的sins初始对准方法 | |
| CN104655131B (zh) | 基于istssrckf的惯性导航初始对准方法 | |
| CN109000642A (zh) | 一种改进的强跟踪容积卡尔曼滤波组合导航方法 | |
| CN107796391A (zh) | 一种捷联惯性导航系统/视觉里程计组合导航方法 | |
| CN101706284B (zh) | 提高船用光纤陀螺捷联惯导系统定位精度的方法 | |
| Georgy et al. | Vehicle navigator using a mixture particle filter for inertial sensors/odometer/map data/GPS integration | |
| CN109931955B (zh) | 基于状态相关李群滤波的捷联惯性导航系统初始对准方法 | |
| CN105865452B (zh) | 一种基于间接卡尔曼滤波的移动平台位姿估计方法 | |
| CN101949703A (zh) | 一种捷联惯性/卫星组合导航滤波方法 | |
| CN104061934A (zh) | 基于惯性传感器的行人室内位置跟踪方法 | |
| CN106840211A (zh) | 一种基于kf和stupf组合滤波的sins大方位失准角初始对准方法 | |
| CN107270893A (zh) | 面向不动产测量的杆臂、时间不同步误差估计与补偿方法 | |
| CN101216321A (zh) | 捷联惯性导航系统的快速精对准方法 | |
| CN111982106A (zh) | 导航方法、装置、存储介质及电子装置 | |
| CN103379619A (zh) | 一种定位方法和系统 | |
| CN103674064B (zh) | 捷联惯性导航系统的初始标定方法 | |
| CN107702712A (zh) | 基于惯性测量双层wlan指纹库的室内行人组合定位方法 | |
| CN105806363A (zh) | 基于srqkf的sins/dvl水下大失准角对准方法 | |
| CN101246012A (zh) | 一种基于鲁棒耗散滤波的组合导航方法 | |
| CN109507706B (zh) | 一种gps信号丢失的预测定位方法 | |
| CN104880201A (zh) | Mems陀螺自动标定方法 | |
| Luo et al. | A position loci-based in-motion initial alignment method for low-cost attitude and heading reference system |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| C06 | Publication | ||
| PB01 | Publication | ||
| C10 | Entry into substantive examination | ||
| SE01 | Entry into force of request for substantive examination | ||
| GR01 | Patent grant | ||
| GR01 | Patent grant |