本发明涉及一种用于动基座重力仪的动态精度评估方法，通过将动基座重力仪安装于特殊设计的试验车上，并选择合适的试验路段进行往复跑车测量，通过对比在同一路段上的多次重力异常值，可以实现动基座重力仪内符合精度的准确、有效评估，且方法简单，易于实现，成本低，效率高，此外本发明重力异常数据的处理方法，抗干扰能力强，具有较高的数据处理精度，特别适用于高精度动基座重力仪的重力异常提取，本发明可用于评估动基座重力仪的内外符合测量精度，适用于动基座重力仪的外场跑车测试。
1.一种用于动基座重力仪的动态精度评估方法，其特征在于：包括如下步骤：(1)、在地面架设差分GPS基站系统，所述GPS基站系统包括第一GPS接收机、第一GPS天线、第一记录设备和第一电源；
(2)、将动基座重力仪(1)安装于减震架(2)的减震器(6)上，减震架(2)固定连接在试验车(3)地板上，(3)、安装地面差分GPS移动站系统，所述GPS移动站系统包括第二GPS接收机、第二GPS天线(4)、第二记录设备和第二电源，其中，第二GPS天线(4)固定安装在试验车(3)的顶部，与动基座重力仪(1)的相对位置保持不变；
(4)、对动基座重力仪(1)在校准点进行静态前校，并在选取的试验路段进行试验车(3)的跑车试验，即：使试验车(3)从校准点出发，到达试验路段的另一端后，转弯并沿原路返回至校准点，之后在校准点进行静态后校，完成一次闭环测量；
(5)、重复步骤(4)，形成至少2次闭环测量，每次闭环测量得到一组动基座重力仪(1)记录的试验车(3)的速度增量和角速度增量，并发送给第二记录设备；
(6)、根据所述至少两组试验车(3)的速度增量、角速度增量以及第一记录设备、第二记录设备记录的数据，计算得到至少两组重力异常值；其中第一GPS接收机从第一GPS天线接收数据并发送给第一记录设备，第二GPS接收机从第二GPS天线接收数据并发送给第二记录设备；
(7)、根据所述至少两组重力异常值计算动基座重力仪(1)的内符合精度，并将所述内符合精度与设定的阈值进行比较，若满足阈值要求，则表明动基座重力仪(1)的动态测量精度满足要求；
所述步骤(6)中根据所述至少两组试验车(3)的速度增量、角速度增量以及第一GPS接收机、第二GPS接收机记录的数据，计算得到至少两组重力异常值的具体方法包括如下步骤：(6.1)、将试验车(3)的速度增量、角速度增量投影至载体坐标系下得到速度增量fb和角速度增量(6.2)、进行捷联式海洋重力仪的动基座初始对准，获得从载体坐标系到实际数学平台坐标系的姿态转移矩阵 所述姿态转移矩阵 为对准时间段内最后时刻k时刻的姿态转移矩阵；
(6.3)、根据对准时间段内最后时刻k时刻的姿态转移矩阵 获得导航过程中当前时刻从载体坐标系到实际数学平台坐标系的姿态转移矩阵 东速ve(ti)、北速vn(ti)、纬度lat(ti)和经度lon(ti)，其中ti为当前时刻；
(6.4)、将第一GPS接收机记录的数据设置为基站数据，将第二GPS接收机记录的数据设置为移动站数据，对所述基站数据和移动站数据进行差分处理获得包含载体运动信息的差分GPS数据；
(6.5)、根据捷联惯性系统误差方程，选取状态向量，构建卡尔曼滤波器的系统状态方程，根据步骤(6.3)中当前时刻的东速ve(ti)、北速vn(ti)、纬度lat(ti)和经度lon(ti)，以及步骤(6.4)中差分GPS数据中的当前时刻的东速VE，北速VN、经度λ和纬度L，得到相应东速的差值ve(ti)-VE、北速的差值vn(ti)-VN、经度的差值lon(ti)-λ和纬度的差值lat(ti)-L，作为卡尔曼滤波器的观测量，进行东速误差、北速误差、经度误差、纬度误差和姿态误差的估计；
(6.6)、根据步骤(6.5)中得到的当前时刻的姿态误差修正东北天向比力值，得到修正后的东北天向比力值fn'；
所述步骤(6.6)中修正后的东北天向比力值f '表示如下：
为当前时刻从载体坐标系到实际数学平台坐标系的姿态转移矩阵；
2.根据权利要求1所述的一种用于动基座重力仪的动态精度评估方法，其特征在于：所述步骤(1)中在地面架设差分GPS基站系统时，架设地点选择在空旷的高地，即基站信号不被遮挡或不被环境干扰，并完成至少2小时的监视测量；如果测量范围较大，即测量基线长度超过100公里，则在测线中间位置架设基站，以保证全程GPS信号质量。
3.根据权利要求1所述的一种用于动基座重力仪的动态精度评估方法，其特征在于：所述步骤(1)中在地面架设的差分GPS基站系统还包括支架，其中将第一GPS天线固定在支架上，第一GPS接收机与第一GPS天线通过馈线连接，第一记录设备与第一GPS接收机通过数据线连接，第一电源为第一GPS接收机和第一记录设备供电。
4.根据权利要求1所述的一种用于动基座重力仪的动态精度评估方法，其特征在于：所述步骤(2)中定义动基座重力仪(1)安装的坐标系为载体坐标系oxbybzb，动基座重力仪(1)的三个轴向oxb、oyb和ozb分别指向动基座重力仪(1)的右前上，安装时动基座重力仪(1)的Y轴指向车头方向，即oxb、oyb和ozb指向试验车(3)的右前上。
5.根据权利要求4所述的一种用于动基座重力仪的动态精度评估方法，其特征在于：所述减震架(2)中部开设凹槽，且四端对称位置各安装一个减震器(6)，动基座重力仪(1)放置在所述凹槽中，并将重力仪的安装部位与四个减震器(6)固定连接。
6.根据权利要求1所述的一种用于动基座重力仪的动态精度评估方法，其特征在于：所述步骤(4)中对动基座重力仪(1)进行静态前校、后校的时间不少于20分钟；所述步骤(4)中试验车(3)的跑车试验中，试验车(3)保持匀速行驶，速度为40-45km/h。
7.根据权利要求1所述的一种用于动基座重力仪的动态精度评估方法，其特征在于：还包括动基座重力仪(1)外符合精度的测量，具体为：采用手持式静态重力仪，在测线沿途选定若干测量点，测量各点的相对重力异常值，绘制测线段的重力异常曲线，选取的外符合基准点的间距取为3-5km/个。
8.根据权利要求1所述的一种用于动基座重力仪的动态精度评估方法，其特征在于：所述步骤(6.3)中获得导航过程中当前时刻从载体坐标系到实际数学平台坐标系的姿态转移矩阵 东速ve(ti)、北速vn(ti)、纬度lat(ti)和经度lon(ti)的具体方法如下：根据导航过程中起始时刻，即第k+1时刻的速度增量fb(tk+1)、角速度增量 和所述姿态转移矩阵 获得导航过程中试验车在第k+1时刻的姿态转移矩阵东速ve(tk+1)、北速vn(tk+1)、纬度lat(tk+1)和经度lon(tk+1)，根据试验车在第k+
1时刻的姿态转移矩阵 东速ve(tk+1)、北速vn(tk+1)、纬度lat(tk+1)和经度lon(tk+1)，以及第k+2时刻的速度增量fb(tk+2)、角速度增量 获得的试验车在第k+2时刻的姿态转移矩阵 东速ve(tk+2)、北速vn(tk+2)、纬度lat(tk+2)和经度lon(tk+2)，依次类推，获得试验车在导航过程中当前时刻的姿态转移矩阵 东速ve(ti)、北速vn(ti)、纬度lat(ti)和经度lon(ti)。
9.根据权利要求1所述的一种用于动基座重力仪的动态精度评估方法，其特征在于：所述步骤(6.3)中姿态转移矩阵 通过如下方法获得：其中：
其中： 为载体坐标系到实际数学平台坐标系的姿态转移矩阵，L代表大地纬度，ωie为地球自转角速度，t0为对准时间段内的初始对准时刻，tk为对准时间段内的任意时刻；
为载体坐标系到载体惯性坐标系的姿态转移矩阵，具体表达式为：
式中：q0、q1、q2、q3为对准数据段最后时刻k的四元素；
为载体惯性坐标系到经线地心惯性坐标系的姿态转移矩阵，具体表达式为：其中： g为地球重力值，
分别计算Vi(tk1)和Vi(tk2)的值，tk1和tk2分别是对准时段内的两个时刻；
为当前时刻从载体坐标系到实际数学平台坐标系的姿态转移矩阵；
10.根据权利要求9所述的一种用于动基座重力仪的动态精度评估方法，其特征在于：所述载体坐标系到载体惯性坐标系的姿态转移矩 通过如下方法获得：(3.1)、初始对准时刻t0，载体坐标系到载体惯性坐标系的姿态转移矩阵 表示如下：其中：I为3阶单位矩阵，其对应的初始时刻四元素为Q(t0)＝[1000]；
(3.2)、根据t0时刻的四元素Q(t0)和t1时刻的角速度增量 获得t1时刻的四元素其中， Φ＝|Φ|；
(3.3)、根据t1时刻的四元素Q(t1)和t2时刻的角速度增量 获得t2时刻的四元素Q(t2)，依次类推，获得对准数据段最后时刻k时刻的四元素Q(k)＝[q0 q1 q2 q3]，根据Q(k)计算 如下：
11.根据权利要求1～10之一所述的一种用于动基座重力仪的动态精度评估方法，其特征在于：所述步骤(4)中差分GPS数据包括GPS时间、经度λ、纬度L、海拔高、大地高、东北天速度(VE，VN，VU)、东北天加速度、卫星数、PDOP、HDOP、VDOP、质量数Q以及GPS周。
12.根据权利要求1～10之一所述的一种用于动基座重力仪的动态精度评估方法，其特征在于：所述步骤(5)中根据捷联惯性系统误差方程，选取的状态向量XINS为13阶，具体表示如下：其中：δL为纬度误差；
13.根据权利要求1～10之一所述的一种用于动基座重力仪的动态精度评估方法，其特征在于：对所述步骤(7)中计算的重力异常粗值δg采用数字滤波器进行滤波处理，以提高数据精度。
14.根据权利要求13所述的一种用于动基座重力仪的动态精度评估方法，其特征在于：所述数字滤波采用FIR和IIR低通滤波器，截止频率小于0.01Hz；或者采用正反卡尔曼滤波器。
一种用于动基座重力仪的动态精度评估方法
技术领域
[0001] 本发明涉及一种用于动基座重力仪的动态精度评估方法，属于航空航天高精度惯性仪器的测试技术领域。
背景技术
[0002] 发达国家掌握着航空、海洋、卫星等重力测量技术、仪器、重力场数据处理和应用技术，其尖端技术、仪器对我国实行封锁。我国重力测量技术现状远不能满足国防安全、国家经济建设及地球科学研究对精密重力测量的迫切需求。近年来，国内部分高校和科研院所在这一领域进行了研究，并形成了原理样机。由于国外的技术封锁，并不清楚国外对于仪器精度的评估方案。国内动基座重力仪精度评估多采用海洋或航空测量以交叉点的标准差评估内符合精度。但该方法试验成本太高，在地面试验对于重力测量仪器尚没有好的评估方案。
发明内容
[0003] 本发明的目的在于克服现有技术的上述不足，提供一种用于动基座重力仪的动态精度评估方法，通过将捷联式动基座重力仪安装于特殊设计的试验车上，并选择合适的试验路段进行往复跑车测量，通过对比在同一路段上的多次重力异常值，可以实现捷联式动基座重力仪内符合精度的准确、有效评估，且方法简单，易于实现，成本低，效率高。
[0004] 本发明的上述目的主要是通过如下技术方案予以实现的：
[0005] 一种用于动基座重力仪的动态精度评估方法，包括如下步骤：
[0006] (1)、在地面架设差分GPS基站系统，所述GPS基站系统包括第一GPS接收机、第一GPS天线、第一记录设备和第一电源；
[0007] (2)、将动基座重力仪安装于减震架的减震器上，减震架固定连接在试验车地板上，
[0008] (3)、安装地面差分GPS移动站系统，所述GPS移动站系统包括第二GPS接收机、第二GPS天线、第二记录设备和第二电源，其中，第二GPS天线固定安装在试验车的顶部，与动基座重力仪的相对位置保持不变；
[0009] (4)、对动基座重力仪在校准点进行静态前校，并在选取的试验路段进行试验车的跑车试验，即：使试验车从校准点出发，到达试验路段的另一端后，转弯并沿原路返回至校准点，之后在校准点进行静态后校，完成一次闭环测量；
[0010] (5)、重复步骤(4)，形成至少2次闭环测量，每次闭环测量得到一组动基座重力仪(1)记录的试验车(3)的速度增量和角速度增量，并发送给第二记录设备；
[0011] (6)、根据所述至少两组试验车(3)的速度增量、角速度增量以及第一记录设备、第二记录设备记录的数据，计算得到至少两组重力异常值；其中第一GPS接收机从第一GPS天线接收数据并发送给第一记录设备，第二GPS接收机从第二GPS天线接收数据并发送给第二记录设备；
[0012] (7)、根据所述至少两组重力异常值计算动基座重力仪的内符合精度，并将所述内符合精度与设定的阈值进行比较，若满足阈值要求，则表明动基座重力仪的动态测量精度满足要求。
[0013] 在上述用于动基座重力仪的动态精度评估方法中，步骤(1)中在地面架设差分GPS基站系统时，架设地点选择在空旷的高地，即基站信号不被遮挡或不被环境干扰，并完成至少2小时的监视测量；如果测量范围较大，即测量基线长度超过100公里，则在测线中间位置架设基站，以保证全程GPS信号质量。
[0014] 在上述用于动基座重力仪的动态精度评估方法中，所述步骤(1)中在地面架设的差分GPS基站系统还包括支架，其中将第一GPS天线固定在支架上，第一GPS接收机与第一GPS天线通过馈线连接，第一记录设备与第一GPS接收机通过数据线连接，第一电源为第一GPS接收机和第一记录设备供电。
[0015] 在上述用于动基座重力仪的动态精度评估方法中，所述步骤(2)中定义动基座重力仪安装的坐标系为载体坐标系oxbybzb，动基座重力仪的三个轴向oxb、oyb和ozb分别指向动基座重力仪的右前上，安装时动基座重力仪的Y轴指向车头方向，即oxb、oyb和ozb指向试验车的右前上。
[0016] 在上述用于动基座重力仪的动态精度评估方法中，所述减震架中部开设凹槽，且四端对称位置各安装一个减震器，动基座重力仪放置在所述凹槽中，并将重力仪的安装部位与四个减震器固定连接。
[0017] 在上述用于动基座重力仪的动态精度评估方法中，所述步骤(4)中对动基座重力仪进行静态前校、后校的时间不少于20分钟；所述步骤(4)中试验车的跑车试验中，试验车保持匀速行驶，速度为40-45km/h。
[0018] 在上述用于动基座重力仪的动态精度评估方法中，还包括动基座重力仪外符合精度的测量，具体为：采用手持式静态重力仪，在测线沿途选定若干测量点，测量各点的相对重力异常值，绘制测线段的重力异常曲线，选取的外符合基准点的间距取为3-5km/个。
[0019] 在上述用于动基座重力仪的动态精度评估方法中，所述步骤(6)中根据所述至少两组试验车的速度增量、角速度增量以及第一GPS接收机、第二GPS接收机记录的数据，计算得到至少两组重力异常值的具体方法包括如下步骤：
[0020] (1)、将试验车的速度增量、角速度增量投影至载体坐标系下得到速度增量fb和角速度增量
[0021] (2)、进行捷联式海洋重力仪的动基座初始对准，获得从载体坐标系到实际数学平台坐标系的姿态转移矩阵 所述姿态转移矩阵 为对准时间段内最后时刻k时刻的姿态转移矩阵；
[0022] (3)、根据对准时间段内最后时刻k时刻的姿态转移矩阵 获得导航过程中当前时刻从载体坐标系到实际数学平台坐标系的姿态转移矩阵 东速ve(ti)、北速vn(ti)、纬度lat(ti)和经度lon(ti)，其中ti为当前时刻；
[0023] (4)、将第一GPS接收机记录的数据设置为基站数据，将第二GPS接收机记录的数据设置为移动站数据，对所述基站数据和移动站数据进行差分处理获得包含载体运动信息的差分GPS数据；
[0024] (5)、根据捷联惯性系统误差方程，选取状态向量，构建卡尔曼滤波器的系统状态方程，根据步骤(3)中当前时刻的东速ve(ti)、北速vn(ti)、纬度lat(ti)和经度lon(ti)，以及步骤(4)中差分GPS数据中的当前时刻的东速VE，北速VN、经度λ和纬度L，得到相应东速的差值ve(ti)-VE、北速的差值vn(ti)-VN、经度的差值lon(ti)-λ和纬度的差值lat(ti)-L，作为卡尔曼滤波器的观测量进行东速误差、北速误差、经度误差、纬度误差和姿态误差的估计；
[0025] (6)、根据步骤(5)中得到的当前时刻的姿态误差修正东北天向比力值，得到修正后的东北天向比力值fn'；
[0026] (7)、计算重力异常粗值δg，公式如下：
[0027]
[0028] 其中：gb为前校点处的重力基准值；
[0029] fu为fn'中的天向比力值；
[0030] 为前校点处的天向比力初值；
[0031] au为天向运动加速度；
[0032] δaE为厄特弗斯改正；
[0033] δaF为自由空间改正；
[0034] γ0为正常重力改正；
[0035] δgdrift为零点漂移改正。
[0036] 在上述用于动基座重力仪的动态精度评估方法中，所述步骤(3)中获得导航过程中当前时刻的姿态转移矩阵 东速ve(ti)、北速vn(ti)、纬度lat(ti)和经度lon(ti)的具体方法如下：
[0037] 根据导航过程中起始时刻，即第k+1时刻的速度增量fb(tk+1)、角速度增量和所述姿态转移矩阵 获得导航过程中试验车在第k+1时刻的姿态转移矩阵东速ve(tk+1)、北速vn(tk+1)、纬度lat(tk+1)和经度lon(tk+1)，根据试验车在第k+
1时刻的姿态转移矩阵 东速ve(tk+1)、北速vn(tk+1)、纬度lat(tk+1)和经度lon(tk+1)，以及第k+2时刻的速度增量fb(tk+2)、角速度增量 获得的试验车在第k+2时刻的姿态转移矩阵 东速ve(tk+2)、北速vn(tk+2)、纬度lat(tk+2)和经度lon(tk+2)，依次类推，获得试验车在导航过程中当前时刻的姿态转移矩阵 东速ve(ti)、北速vn(ti)、纬度lat(ti)和经度lon(ti)。
[0038] 在上述用于动基座重力仪的动态精度评估方法中，所述步骤(3)中姿态转移矩阵通过如下方法获得：
[0039]
[0040] 其中：
[0041]
[0042] 其中： 为载体坐标系到实际数学平台坐标系的姿态转移矩阵，L代表大地纬度，ωie为地球自转角速度，t0为对准时间段内的初始对准时刻，tk为对准时间段内的任意时刻；
[0043] 为载体坐标系到载体惯性坐标系的姿态转移矩阵，具体表达式为：
[0044]
[0045] 式中：q0q1q2q3为对准数据段最后时刻k的四元素；
[0046] 为载体惯性坐标系到经线地心惯性坐标系的姿态转移矩阵，具体表达式为：
[0047]
[0048] 其中： g为地球重力值，分别计算Vi(tk1)和Vi(tk2)的值，tk1和tk2分别是对准时段内的两个时刻；
[0049]
[0050] Δtk＝tk-t0，
[0051] 为当前时刻从载体坐标系到实际数学平台坐标系的姿态转移矩阵；
[0052] fb(ti)为当前时刻的速度增量。
[0053] 在上述用于动基座重力仪的动态精度评估方法中，所述载体坐标系到载体惯性坐标系的姿态转移矩 通过如下方法获得：
[0054] (3.1)、初始对准时刻t0，载体坐标系到载体惯性坐标系的姿态转移矩阵 表示如下：
[0055]
[0056] 其中：I为3阶单位矩阵，其对应的初始时刻四元素为Q(t0)＝[1 0 0 0]；
[0057] (3.2)、根据t0时刻的四元素Q(t0)和t1时刻的角速度增量 获得t1时刻的四元素 其中， Φ＝|Φ|；
[0058] (3.3)、根据t1时刻的四元素Q(t1)和t2时刻的角速度增量 获得t2时刻的四元素Q(t2)，依次类推，获得对准数据段最后时刻k时刻的四元素Q(k)＝[q0 q1 q2 q3]，根据Q(k)计算 如下：
[0059]
[0060] 在上述用于动基座重力仪的动态精度评估方法中，所述步骤(4)中差分GPS数据包括GPS时间、经度λ、纬度L、海拔高、大地高、东北天速度(VE，VN，VU)、东北天加速度、卫星数、PDOP、HDOP、VDOP、质量数Q以及GPS周。
[0061] 在上述用于动基座重力仪的动态精度评估方法中，所述步骤(5)中根据捷联惯性系统误差方程，选取的状态向量XINS为13阶，具体表示如下：
[0062]
[0063] 其中：δL为纬度误差；
[0064] δλ为经度误差；
[0065] δve、δvn分别为东速误差和北速误差；
[0066] φe、φn和φu分别为三个姿态误差角；
[0067] εx、εy和εz为激光陀螺的零位；
[0068] 和 为加速度计零位。
[0069] 在上述用于动基座重力仪的动态精度评估方法中，所述步骤(6)中修正后的东北n'天向比力值f 表示如下：
[0070]
[0071] 其中：φ×为反对称阵，
[0072] 为当前时刻从载体坐标系到实际数学平台坐标系的姿态转移矩阵；
[0073] fb(ti)为当前时刻的速度增量；
[0074] ΔT为系统采样间隔时间。
[0075] 在上述用于动基座重力仪的动态精度评估方法中，对所述步骤(7)中计算的重力异常粗值δg采用数字滤波器进行滤波处理，以提高数据精度。
[0076] 在上述用于动基座重力仪的动态精度评估方法中，所述数字滤波采用FIR和IIR低通滤波器，截止频率小于0.01Hz；或者采用正反卡尔曼滤波器。
[0077] 本发明与现有技术相比的优点在于：
[0078] (1)、本发明通过将捷联式动基座重力仪安装于特殊设计的试验车上，并选择合适的试验路段进行往复跑车测量，通过对比在同一路段上的多次重力异常值，可以实现捷联式动基座重力仪内符合精度的准确、有效评估；
[0079] (2)、本发明针对捷联式动基座重力仪的陆地跑车测量环境，设计了专用的试验车工装，能有效衰减由于路面颠簸和环境带来了冲击干扰。
[0080] (3)、本发明在地面环境下，不易形成交叉点，通过往返测线的测量方式获得测量路段的重力异常重复测量值，进而可按经纬度进行内符合评估；
[0081] (4)、本发明对重力异常值的计算方法，抗干扰能力强，具有较高的数据处理精度，可以作为高精度捷联式动基座重力仪异常提取的数据处理方法；
[0082] (5)、本发明对重力异常值的计算过程中，采用kalman滤波器，估计惯性导航解算中的姿态误差，对姿态矩阵和垂向比力分量进行修正，得到更为精确的比力信息，再由差分GPS技术得提供的位置、速度和高度信息计算重力各改正项，最后经低通滤波得到沿测线的重力异常信息，本发明可用于捷联式动基座重力仪的车载重力测量数据处理，抗干扰能力强，具有较高的数据处理精度，特别适用于捷联式动基座重力仪的车载重力测量数据处理；
[0083] (6)、本发明对重力异常值的计算过程中，在初始对准阶段，采用了惯性凝固假设，进行动基座初始对准，可以有效减小由于载体晃动带来的姿态误差，提高初始对准精度；
[0084] (7)、本发明对重力异常值的计算过程中，采用差分GPS技术处理车载试验过程中的GPS数据，获得较高精度的载体运动信息；采用组合导航技术获取天向比力信息，经过改正计算后，通过数字低通方式去除高频噪声得到高精度的重力异常信号。
[0085] (8)、本发明对重力异常值的计算过程中，可以采用数字低通滤波器对计算得到的重力异常粗值去噪处理，截止频率小于0.01Hz，进一步提高了重力异常信号的精度。
[0086] (9)、本发明采用相对式静态重力仪在沿途测线段中选取若干重力测量点，在测线上形成若干离散的测量点作为测线段重力异常的基准值，通过对比跑车测量与静态测量基准点的重力异常值，可以评估仪器的外符合精度。
[0087] (10)、本发明方法实现方便，成本低，并且条件许可的情况下，可以对多台套重力仪同时进行测试，测试效率高。
附图说明
[0088] 图1为本发明动基座重力仪的动态精度评估方法的工作流程图；
[0089] 图2为本发明动基座重力仪在试验车上的安装示意图；
[0090] 图3为本发明重力仪与减震架的安装示意图。
[0091] 图4为本发明动基座重力仪的安装坐标图；
[0092] 图5为本发明重力异常数据处理方法的工作流程图；
[0093] 图6为本发明重力异常数据处理方法中载体坐标系、导航坐标系和实际数学平台坐标系之间关系示意图；
[0094] 图7为本发明重力异常数据处理方法中针对的捷联式海洋重力仪的硬件安装示意图。
具体实施方式
[0095] 下面结合附图和具体实施例对本发明作进一步详细的描述：
[0096] 本发明提出的用于动基座重力仪的动态精度评估方法，采用车载试验方式，沿预定测线往返测量，通过测线重复测量结果评估仪器的内符合精度，还可以通过手持式相对重力仪引入外部基准点的方式评估仪器的外符合精度。其具体实施方式如下所述：
[0097] 步骤(一)、在地面架设差分GPS基站系统，所述GPS基站系统包括第一GPS接收机、第一GPS天线、支架、第一记录设备和第一电源。其中将第一GPS天线固定在支架上，第一GPS接收机与第一GPS天线通过馈线连接，第一记录设备与第一GPS接收机通过数据线连接，第一电源为第一GPS接收机和第一记录设备供电。
[0098] 在地面架设差分GPS基站系统时，架设地点选择在空旷的高地，即基站信号不被遮挡或不被环境干扰，并完成至少2小时的监视测量；如果测量范围较大，即测量基线长度超过100公里，则在测线中间位置架设基站，以保证全程GPS信号质量。
[0099] 步骤(二)、将动基座重力仪1安装于减震架2的减震器6上，减震架2固定连接在试验车3地板上。如图2所示为本发明动基座重力仪在试验车上的安装示意图。如图3所示为本发明重力仪与减震架的安装示意图，由图可知，减震架2中部开设凹槽，且四端对称位置各安装一个减震器6，动基座重力仪1放置在所述凹槽中，并通过螺钉将重力仪的安装部位与四个减震器6固定连接。
[0100] 当路面颠簸时，车辆在行驶过程中会产生垂向运动加速度，该运动会传导至减震架。由于重力仪是通过4个减震器与减震支架固连，因此，能在很大程度上隔绝由路面颠簸导致的垂向运动加速度。
[0101] 定义动基座重力仪1安装的坐标系为载体坐标系oxbybzb，动基座重力仪1的三个轴向oxb、oyb和ozb分别指向动基座重力仪1的右前上，安装时动基座重力仪1的Y轴指向车头方向，即oxb、oyb和ozb指向工程车的右前上。如图4所示为本发明动基座重力仪的安装坐标图。
[0102] 步骤(三)、安装地面差分GPS移动站系统，所述GPS移动站系统包括第二GPS接收机、第二GPS天线4、第二记录设备和第二电源，其中，第二GPS天线4固定安装在试验车3的顶部，与动基座重力仪1的相对位置保持不变；
[0103] 步骤(四)、对仪器进行充分预热后，在试验路段的一端选取合适的校准点，进行不少于20分钟的静态前校。
[0104] 步骤(五)、在选取的试验路段进行试验车3的跑车试验，即使试验车3从校准点出发，到达试验路段的另一端后，转弯并沿原路返回至校准点，之后在校准点进行静态后校，完成一次闭环测量。在试验时，试验车3保持匀速行驶，避免急起急停。例如对于长度为60km的测量基线载体的行驶速度为40-45km/h。
[0105] 步骤(六)、重复步骤(五)，形成至少2次闭环测量，每次闭环测量得到一组动基座重力仪1记录的试验车(3)的速度增量和角速度增量；并发送给第二记录设备；
[0106] 步骤(七)、根据所述至少两组试验车3的速度增量、角速度增量以及第一记录设备、第二记录设备记录的数据，计算得到至少两组重力异常值；其中第一GPS接收机从第一GPS天线接收数据并发送给第一记录设备，第二GPS接收机从第二GPS天线接收数据并发送给第二记录设备。
[0107] 重力异常值具体的计算方法如下：
[0108] 如图5所示为本发明重力异常数据处理方法的工作流程图，本发明重力异常数据处理方法具体包括如下步骤：
[0109] (1)、定义动基座重力仪1安装的坐标系为载体坐标系oxbybzb，动基座重力仪1的三个轴向oxb、oyb和ozb分别指向动基座重力仪1的右前上，安装时动基座重力仪1的Y轴指向车头方向，即oxb、oyb和ozb指向工程车的右前上，如图7所示为本发明所针对的捷联式重力仪的硬件安装示意图。
[0110] 重力仪实时记录载体的速度增量和角速度增量，并投影至载体坐标系下得到速度增量fb和角速度增量
[0111] (2)采用了惯性凝固假设，进行捷联式动基座重力仪的动基座初始对准，获得从载体坐标系到实际数学平台坐标系的姿态转移矩阵 所述姿态转移矩阵 为对准时间段内最后时刻k时刻的姿态转移矩阵。
[0112] 对载体坐标系(b系)到实际数学平台坐标系(n系)的姿态转移矩阵进行初始化设定；重力数据处理中，导航坐标系通常选为地理坐标系。如图6所示为本发明载体坐标系、导航坐标系和实际数学平台坐标系之间关系示意图。
[0113] 该方法用到的坐标系定义如下：
[0114] a)经线地球坐标系e：原点位于地心，oze轴沿地球自转轴方向，oxe轴位于赤道平面内，从地心指向重力仪所在点经线，oye轴在赤道平面内，oxe、oye、oze轴构成右手坐标系。
[0115] b)经线地心惯性坐标系i：在对准起始时刻t0时刻将经线地球坐标系oxeyeze惯性凝固后形成的坐标系。
[0116] c)导航坐标系n’：原点位于捷联式动基座重力仪中心，ox轴指向东(E)，oy轴指向北(N),oz轴指向天(U)。
[0117] d)实际数学平台坐标系n：坐标系Ox1y1z1，近似指向东北天，与理想导航坐标n’系之间存在失准角，例如，水平失准角(φe、φn)为0.005°，方位失准角φu为0.08°。
[0118] e)载体系b系：坐标系Oxbybzb，原点位于载体质心，xb、yb和zb指向载体右前上。
[0119] f)载体惯性坐标系ib0：在对准起始时刻t0时刻将载体坐标系oxbybzb经惯性凝固后的坐标系。
[0120] 动基座初始对准算法中，姿态阵分散成4个矩阵求取。设对准点的纬度为L，则姿态转移矩阵 通过如下方法获得：
[0121]
[0122] 其中：
[0123]
[0124] 其中： 为载体坐标系到实际数学平台坐标系的姿态转移矩阵，L代表大地纬度，ωie为地球自转角速度，t0为对准时间段内的初始对准时刻，tk为对准时间段内的任意时刻；
[0125] 为载体坐标系到载体惯性坐标系的姿态转移矩阵，具体表达式为：
[0126]
[0127] 式中：q0q1q2q3为对准数据段最后时刻k的四元素。
[0128] 载体坐标系到载体惯性坐标系的姿态转移矩 通过如下方法获得：
[0129] (a)、初始对准时刻t0，载体坐标系到载体惯性坐标系的姿态转移矩阵 表示如下：
[0130]
[0131] 其中：I为3阶单位矩阵，其对应的初始时刻四元素为Q(t0)＝[1 0 0 0]；
[0132] (b)、根据t0时刻的四元素Q(t0)和t1时刻的角速度增量 获得t1时刻的四元素 其中， Φ＝|Φ|；
[0133] (c)、根据t1时刻的四元素Q(t1)和t2时刻的角速度增量 获得t2时刻的四元素Q(t2)，依次类推，获得对准数据段最后时刻k时刻的四元素Q(k)＝[q0 q1 q2 q3]，根据Q(k)计算 如下：
[0134]
[0135] 为载体惯性坐标系到经线地心惯性坐标系的姿态转移矩阵，具体表达式为：
[0136]
[0137] 其中： g为地球重力值，分别计算Vi(tk1)和Vi(tk2)的值，tk1和tk2分别是对准时段内的两个时刻；
[0138]
[0139] Δtk＝tk-t0，
[0140] 为当前时刻从载体坐标系到实际数学平台坐标系的姿态转移矩阵；
[0141] fb(ti)为当前时刻的速度增量。
[0142] (3)、根据对准时间段内最后时刻k时刻的姿态转移矩阵 获得载体导航过程中当前时刻从载体坐标系到实际数学平台坐标系的姿态转移矩阵 东速ve
(ti)、北速vn(ti)、纬度lat(ti)和经度lon(ti)，其中ti为当前时刻；具体方法如下：
[0143] 根据导航过程中起始时刻，即第k+1时刻的速度增量fb(tk+1)、角速度增量和所述姿态转移矩阵 获得导航过程中试验车在第k+1时刻的姿态转移矩阵东速ve(tk+1)、北速vn(tk+1)、纬度lat(tk+1)和经度lon(tk+1)；
[0144] 根据载体在第k+1时刻的姿态转移矩阵 东速ve(tk+1)、北速vn(tk+1)、纬度lat(tk+1)和经度lon(tk+1)，以及第k+2时刻的速度增量fb(tk+2)、角速度增量 获得的载体在第k+2时刻的姿态转移矩阵 东速ve(tk+2)、北速vn(tk+2)、纬度lat(tk+2)和经度lon(tk+2)；
[0145] 根据载体在第k+2时刻的姿态转移矩阵 东速ve(tk+2)、北速vn(tk+2)、纬b度lat(tk+2)和经度lon(tk+2)，以及第k+3时刻的加速度增量f(tk+3)、角速度增量 获得的载体在第k+3时刻的姿态转移矩阵 东速ve(tk+3)、北速vn(tk+3)、纬度lat(tk+3)和经度lon(tk+3)。
[0146] 依次类推，获得载体在导航过程中当前时刻的姿态转移矩阵 东速ve(ti)、北速vn(ti)、纬度lat(ti)和经度lon(ti)。按照上述方法可以获得载体在导航过程中各时刻的姿态转移矩阵 东速ve(t)、北速vn(t)、纬度lat(t)和经度lon(t)。
[0147] (4)、将第一GPS接收机记录的数据设置为基站数据，将第二GPS接收机记录的数据设置为移动站数据，对所述基站数据和移动站数据进行差分处理获得包含载体运动信息的差分GPS数据。也就是根据载体测量期间记录的GPS的基站数据和移动站数据，采用waypoint软件中的差分GPS技术，获得差分GPS信息，包括GPS时间、经度λ、纬度L、海拔高、大地高、东北天速度(VE，VN，VU)、东北天加速度、卫星数、PDOP、HDOP、VDOP、质量数Q以及GPS周。获得差分GPS数据后，进行步骤(5)。
[0148] (5)、根据捷联惯性系统误差方程，选取状态向量，构建卡尔曼滤波器的系统状态方程，根据步骤(3)中当前时刻的东速ve(ti)、北速vn(ti)、纬度lat(ti)和经度lon(ti)，以及步骤(4)中差分GPS数据中的当前时刻的东速VE，北速VN、经度λ和纬度L，得到相应东速的差值ve(ti)-VE、北速的差值vn(ti)-VN、经度的差值lon(ti)-λ和纬度的差值lat(ti)-L，作为卡尔曼滤波器的观测量进行东速误差、北速误差、经度误差、纬度误差和姿态误差的估计。
[0149] 根据kalman滤波器估计修正各项参数误差，特别是位置、速度和姿态误差，修正分为开环和闭环修正两种方式。完成后，进入步骤(6)。
[0150] 选取的状态向量XINS为13阶，具体表示如下：
[0151]
[0152] 根据测量环境，忽略部分非主要误差参数，采用了如下的捷联惯性系统的误差方程
[0153]
[0154]
[0155]
[0156]
[0157]
[0158]
[0159]
[0160] 式中，δL为纬度误差；
[0161] δλ为经度误差；
[0162] δve、δvn分别为东、北速误差；
[0163] φe、φn和φu分别为三个姿态误差角，通常情况下，φ为小量；
[0164] εx、εy和εz为激光陀螺的零位；
[0165] 和 为加速度计零位；
[0166] Tij(i＝1,2,3；j＝1,2,3)为姿态阵 的元素。
[0167] (6)、根据步骤(5)中得到的当前时刻姿态误差，修正东北天向比力值，得到修正后的东北天向比力值fn'，完成修正后进入步骤(7)。
[0168]
[0169] 其中：φ×为反对称阵，
[0170] 为当前时刻从载体坐标系到实际数学平台坐标系的姿态转移矩阵；
[0171] fb(ti)为当前时刻的速度增量；
[0172] ΔT为系统采样间隔时间。
[0173] (7)、计算重力异常粗值δg，公式如下：
[0174]
[0175] 式中：
[0176] gb为前校点处的重力基准值，如果没有基准值，则采用前校段的天向比力值代替；
[0177] fu为fn'中的天向比力值；
[0178] 为前校点处的天向比力初值；
[0179] au为天向运动加速度；
[0180] δaE为厄特弗斯改正；
[0181] δaF为自由空间改正；
[0182] γ0为正常重力改正；
[0183] δgdrift为零点漂移改正，零点漂移改正为每个航次重力测量时不同时间在同一点的观测值变化的改正。
[0184] 其中计算重力异常各改正项，包括厄特弗斯改正、天向运动加速度、正常重力改正、自有空间改正、零点漂移改正，各改正项计算公式如下。厄特弗斯：
[0185]
[0186] 天向运动加速度：
[0187] 可以根据GPS提供的伪距、载波相位和多普勒频移观测值，以及其单差、双差组建观测方程，利用最小二乘法，求出载体位置、速度和加速度。然后利用差分GPS数据，采用位置差分、速度差分或者载波相位差分等方法计算天向运动加速度。
[0188] 正常重力：
[0189] γ0＝9.780327(1+0.0053024sin2L-0.0000058sin22L)
[0190] 自由空间：
[0191]
[0192] 零点漂移改正：
[0193] 重力测量的零点漂移率可按线性化近似计算，有
[0194]
[0195] 其中，C为此次测量的零点漂移变化率，f1u和 分别为前后校基准点处的天向比力值，g1b和 分别为前后校基准点处的重力场，t1和t0分别对应f1u和 的观测时间。
[0196] 则零点漂移改正值为
[0197] δgdrift＝C·(t-t0)
[0198] 式中，C为此次测量的零点漂移变化率，Δti为第i个测点的测量时刻与在基准点前校时刻的时间差。
[0199] 完成计算后进入步骤(8)。
[0200] (8)采用数字滤波器对δg进行滤波，得到高精度重力异常信号，采用滤波器可采用FIR和IIR滤波器，截止频率小于0.01Hz；也可采用正反kalman滤波器。
[0201] 步骤(八)、根据所述至少两组重力异常值计算动基座重力仪1的内符合精度，并将所述内符合精度与设定的阈值进行比较，若满足阈值要求(即内符合精度小于或等于阈值)，则表明动基座重力仪1的动态测量精度满足要求。
[0202] 如果具备测量条件，可采用CG-5等手持式静态重力仪，在测线沿途选定若干测量点，测量各点的相对重力异常值，绘制测线段的重力异常曲线。建议重力异常测点的间距5km左右，可用于评估仪器的外符合精度。选取的外符合基准点的间距取为3-5km/个，60km试验路段选择的外基准点数为12-20个。
[0203] 以上所述，仅为本发明最佳的具体实施方式，但本发明的保护范围并不局限于此，任何熟悉本技术领域的技术人员在本发明揭露的技术范围内，可轻易想到的变化或替换，都应涵盖在本发明的保护范围之内。
[0204] 本发明说明书中未作详细描述的内容属于本领域专业技术人员的公知技术。
