一种基于主子惯导结合的行进间炮口振动测量方法

1.本发明涉及炮口振动测量技术领域，特别是一种基于主子惯导结合的行进间炮口振动测量方法。


背景技术：

2.火炮行进作业时，身管炮口处的振动会影响火炮的射击精度，炮口振动测量研究是如今一个难题。炮口振动研究最早开展于国外，komnkov建立的数学方程表述了身管振动与射击精度之间的关系，证明身管振动是影响射击精度的主要因素之一；苏联火炮设计专家奥尔洛夫与马利科夫等人在1982年给出了估算身管固有频率的计算方法。国内于会杰、王德石等人将身管简化成梁类的结构模型，在此基础上利用梁的振动理论，研究身管振动响应；徐达等人采用线激光散斑场光电检测方法，研究了身管振动信号与输出电流数学关系，提出了一种基于身管表面线激光散斑效应的炮口振动测量的方法，实现炮口振幅测量。
3.在多数研究中，以位移量变化参数开展研究多采用非接触式测量，不足之处即是非接触式测量仪器安装不方便，采样数据较为困难，仪器校准也较为麻烦。


技术实现要素：

4.本发明的目的在于提出一种安装简便、稳定性高、抗干扰能力强、准确性高的基于主子惯导结合的行进间炮口振动测量方法。
5.实现本发明目的的技术解决方案为：一种基于主子惯导结合的行进间炮口振动测量方法，包括以下步骤：
6.步骤1、根据火炮平台，设计主、子惯导安装位置，包括主惯导和mems子惯导；
7.步骤2、根据火炮行进间主惯导、mems子惯导传感器输出，设计基于陀螺仪角增量更新的姿态更新算法；
8.步骤3、在传递对准阶段，根据火炮身管特性进行外杆臂误差与挠曲形变的动态建模，对杆臂效应及挠曲形变进行补偿；
9.步骤4、将mems子惯导解算出的速度、姿态，相对于主惯导解算出的速度、姿态偏差量作为mems子惯导的量测信息，建立以速度加姿态匹配方式的传递对准卡尔曼滤波器；
10.步骤5、在导航更新阶段，mems子惯导利用传递对准结束时刻的对准误差，作为导航解算的初始误差，进行导航状态更新，完成炮口的振动测量。
11.本发明与现有技术相比，其显著优点为：(1)在炮塔处安装高精度主惯导，在身管炮口处安装mems子惯导的接触式测量方式，相对传统光学等非接触式测量方式，安装方便、稳定性高，避免了传统光学方式的易受干扰、安装困难、仪器校准繁琐等缺点；(2)在初始对准阶段，采用高精度主惯导对mems子惯导进行传递对准，属于惯导精对准范畴，可以保证mems子惯导有准确的初始惯导信息，提高了炮口振动测量的可靠性和准确性；(3)在导航更新阶段，子惯导利用传递对准结束时刻的对准误差，作为导航解算的初始误差，利用姿态更新算法不断迭代解算，可以保证mems子惯导的位姿输出准确性，提高了炮口振动测量方法
的实时性和可靠性。
附图说明
12.图1是本发明主子惯导安装的结构示意图。
13.图2是本发明基于主子惯导结合的行进间炮口振动测量方法的流程图。
14.图3是本发明中主子惯导杆臂的原理示意图。
15.图4是本发明中kalman滤波的流程示意图。
16.图5是本发明中kalman滤波计算回路与增益计算回路的流程示意图。
具体实施方式
17.本发明一种基于主子惯导结合的行进间炮口振动测量方法，包括以下步骤：
18.步骤1、根据火炮平台，设计主、子惯导安装位置，包括主惯导和mems子惯导；
19.步骤2、根据火炮行进间主惯导、mems子惯导传感器输出，设计基于陀螺仪角增量更新的姿态更新算法；
20.步骤3、在传递对准阶段，根据火炮身管特性进行外杆臂误差与挠曲形变的动态建模，对杆臂效应及挠曲形变进行补偿；
21.步骤4、将mems子惯导解算出的速度、姿态，相对于主惯导解算出的速度、姿态偏差量作为mems子惯导的量测信息，建立以速度加姿态匹配方式的传递对准卡尔曼滤波器；
22.步骤5、在导航更新阶段，mems子惯导利用传递对准结束时刻的对准误差，作为导航解算的初始误差，进行导航状态更新，完成炮口的振动测量。
23.进一步地，步骤1中，根据火炮平台，设计主、子惯导安装位置，包括主惯导和mems子惯导，具体如下：
24.所述主惯导固连于炮塔，保持水平状态，炮塔进行航向运动带动主惯导一同运动；
25.所述mems子惯导固连于火炮身管前端，用于测量火炮身管前端炮口的位姿。
26.进一步地，步骤2中根据火炮行进间主惯导、mems子惯导传感器输出，设计基于陀螺仪角增量更新的姿态更新算法，具体如下：
27.步骤2.1、令惯导静止放置，获取加速度计数据，用重力分量关系获得俯仰角θ和滚转角γ，利用地磁传感器获得导航坐标系下的偏航角完成惯导系统的初始对准；
28.步骤2.2、设置初始零点，利用四元数以及陀螺仪的角增量，对惯导载体的姿态进行更新。
29.进一步地，所述步骤2.1具体如下：
30.步骤2.1.1、利用加速度三轴分量估算俯仰角θ和滚转角γ：
[0031][0032][0033]
其中，为载体坐标系下第k次加速度计的x、y、z三轴输出值；
[0034]
步骤2.1.2、利用地磁传感器获得载体坐标系下的磁场强度结合获得的俯仰角θ和滚转角γ，经过转换得导航坐标系下得磁场强度进一步计算得到偏航角为：
[0035][0036]
其中，θk、γk分别表示第k次解算出的俯仰角、滚转角；
[0037][0038]
其中，分别表示导航坐标系下x、y、z三轴对应的磁场强度；
[0039]
步骤2.1.3、在静止条件下，利用三个欧拉角获得从载体坐标系到导航坐标系的方向余弦矩阵，完成惯导的初始对准。
[0040]
进一步地，所述步骤2.2，具体如下：
[0041]
步骤2.2.1、由陀螺仪的更新数据计算得到角增量δ：
[0042][0043]
式中，ω
x
、ωy、ωz表示仪输出的三轴角速度标量值，tm为采样时间；
[0044]
步骤2.2.2、利用角增量对四元数进行更新：
[0045][0046]
式中，q
1|k+1
、q
2|k+1
、q
3|k+1
、q
4|k+1
分别为k+1时刻四元数的q1、q2、q3、q4值，q
1|k
、q
2|k
、q
3|k
、q
4|k
分别为k时刻四元数的q1、q2、q3、q4值；
[0047]
步骤2.2.3、四元数进行归一化处理：
[0048][0049]
其中，q1'
|k+1
、q'
2|k+1
、q'
3|k+1
、q'
4|k+1
分别为q
1|k+1
、q
2|k+1
、q
3|k+1
、q
4|k+1
归一化处理后
的值；进行归一化使得四元数变为单位四元数，四个值的平方和为1；
[0050]
步骤2.2.4、根据单位四元数得到方向余弦矩阵，公式为：
[0051][0052]
步骤2.2.5、根据方向余弦矩阵得到欧拉角如下：
[0053][0054]
步骤2.2.6、根据方向余弦矩阵信息和比力信息，得到当前加速度三轴的分量，区分三轴分量能够对重力分量进行补偿，除去重力分量获得导航坐标系下的运动加速度，对运动加速度进行牛顿力学积分得到速度和位置，公式为：
[0055][0056]
其中，为根据补偿重力分量后的运动加速度计输出值，表示从k-1时刻到k时刻，从载体坐标系b系到导航坐标系n系的坐标变换矩阵。
[0057]
进一步地，步骤3中的外杆臂误差补偿，具体如下：
[0058]
定义惯性坐标系为oixiyizi，载体坐标系为obxbybzb，ob是载体的摇摆中心，子惯导的加速度计安装在载体坐标系中的固定点p，为载体坐标系原点的位置矢量，为p点相对于惯性坐标系原点的位置矢量，为p点相对于载体坐标系原点的位置矢量，则p点相对于惯性坐标系的线加速度表达式为：
[0059][0060]
上式后面两项就是由于杆臂效应引起的惯导敏感到的杆臂加速度基本表达式；
[0061]
其中，表示陀螺仪输出的在惯性系下的角速度；
[0062]
定义表示杆臂加速度，则杆臂效应误差的基本方程为：
[0063][0064]
其中，表示杆臂加速度，为切向加速度，为向心加速度。
[0065]
进一步地，步骤3中所述的挠曲变形的动态建模，具体如下：
[0066]
设定动态变形角为λ(t)，动态变形角速率为wf(t)，即有：展开为：
[0067][0068]
其中，w
fx
(t)、w
fy
(t)、w
fz
(t)分别表示x，y，z三轴动态变形角速率，(t)分别表示x，y，z三轴动态变形角速率，分别为x，y，z三轴动态变形角的一阶导数；
[0069]
动态变形角速率运动方程：
[0070][0071]
其中，令i＝x,y,z，βi＝2.146/τi，τi为相关时间；wi(t)是白噪声，如果是有色噪声，则需要进行白化处理，白噪声wi(t)的方差满足(t)的方差满足表示三个弹性变形角λ
x
,λy,λz的方差。
[0072]
进一步地，步骤4中的mems子惯导解算出的速度、姿态，相对于主惯导解算出的速度、姿态偏差量，具体如下：
[0073]
主、子惯导的速度差微分方程为：
[0074][0075]
式中，表示主子惯导的速度差，为子惯导系统的姿态矩阵，φ
×
表示主、子惯导的初始安装误差角对应的反对称矩阵，ξ表示主、子惯导的安装误差角估计值；表示主惯导在子惯导坐标系下感受到的比力；表示地球在导航坐标系下的自转角速率；表示子惯导加速度计误差；
[0076]
姿态量测方程为：
[0077][0078]
式中，φ表示主、子惯导解算的姿态偏差量；z
φ
表示主、子惯导去掉静态变形角与动态变形角影响的姿态差；μ表示静态变形角，包括安装误差；θ表示动态变形角，即挠曲变形。
[0079]
进一步地，步骤4中所述的速度加姿态匹配方式的卡尔曼滤波模型，具体如下：
[0080]
(1)状态一步预测
[0081][0082]
其中，为k时刻状态的一步预测；φ
k,k-1
为k-1时刻到k时刻状态一步转移矩阵；为k-1时刻状态的最优估计；
[0083]
(2)状态一步预测均方误差阵
[0084][0085]
其中，p
k|k-1
为状态一步预测均方误差阵；p
k-1
为k-1时刻状态最优估计的均方误差
阵；(3)滤波增益
[0086][0087]
其中，系数矩阵kk称为滤波增益；为k时刻的系统量测矩阵的转置；pk为k时刻状态最优估计的均方误差阵；
[0088]
(4)状态估计
[0089][0090]
其中，为k时刻状态的最优估计；zk为k时刻的量测；hk为k时刻的系统量测矩阵；
[0091]
(5)状态估计均方误差阵
[0092]
pk＝(i-k
khk
)p
k|k-1
ꢀꢀ
(21)。
[0093][0094]
进一步地，所述步骤5具体如下：
[0095]
传递对准结束后，子惯导从状态估计量中获取当前陀螺仪的零漂及加速度计的零偏，以及主、子惯导的失准角；子惯导利用获取的三个欧拉角得到从载体坐标系到导航坐标系的方向余弦矩阵，结合步骤2中的姿态更新算法，实现子惯导的实时位置姿态测量；
[0096]
由于子惯导与炮口固连，所以子惯导的实时位置姿态信息即为炮口的实时振动情况。
[0097]
下面结合附图和具体实施例，对本发明做进一步的详细说明。
[0098]
实施例
[0099]
结合图1，本发明一种基于主子惯导结合的行进间炮口振动测量方法，包括以下步骤：
[0100]
步骤1、根据火炮平台，设计主、子惯导安装位置，包括高精度主惯导1和mems子惯导2；
[0101]
步骤2、根据火炮行进间主、子惯导传感器输出，设计基于陀螺仪角增量更新的姿态更新算法；
[0102]
步骤3、在传递对准阶段，根据火炮身管特性进行外杆臂误差与挠曲形变的动态建模，对杆臂效应及挠曲形变进行补偿；
[0103]
步骤4、将mems子惯导解算出的速度、姿态相对于高精度主惯导解算出的速度、姿态偏差量作为子惯导的量测信息，建立以速度加姿态匹配方式的传递对准卡尔曼滤波器；
[0104]
步骤5、在导航更新阶段，子惯导利用传递对准结束时刻的对准误差，作为导航解算的初始误差，进行导航状态更新，完成炮口的振动测量。
[0105]
进一步地，步骤1所述的根据火炮平台，设计主、子惯导安装位置，包括高精度主惯导1和mems子惯导2，具体如下：
[0106]
所述高精度抗高过载惯导1固连于炮塔4，保持水平状态，炮塔4进行航向运动带动高精度抗高过载惯导1一同运动；
[0107]
所述mems子惯导2固连于火炮身管3前端，用于测量火炮身管3前端炮口的位姿。
[0108]
进一步地，步骤2所述的根据火炮行进间主、子惯导传感器输出，设计基于陀螺仪角增量更新的姿态更新算法，具体如下：
[0109]
在不同的导航系统中，对坐标系的选取也各不相同，特别是对导航坐标系选取。默
认ox
nynzn
作为导航坐标系。
[0110]
(1)地理坐标系
[0111]
地理坐标系，原点位于载体质心，其中一个坐标轴沿当地地理垂线的方向，另外两轴在当地水平面内分别沿当地经线和纬线的切线方向。选取东北天(enu)地理坐标系，即:x轴指向东，y轴指向北，z轴垂直于当地水平面，沿当地垂线向上。
[0112]
(2)载体坐标系
[0113]
载体坐标系oxbybzb，(原)点与载体质心重合。xb沿载体横轴向右，yb沿载体纵轴向前，zb沿载体竖轴向上，又称“右前上”坐标系。
[0114]
(3)导航坐标系
[0115]
由于涉及多个坐标系系统，需要将多种坐标系转换到同一坐标系进行求解，用ox
nynzn
表示。
[0116]
(4)主惯导坐标系
[0117]
主惯导坐标系选取与载体坐标系一致，用oxmymzm表示。
[0118]
(5)子惯导坐标系
[0119]
子惯导坐标系用ox
syszs
表示，惯导安装在炮口位置。
[0120]
步骤2.1、首先需要对主、子惯导进行初始对准，令惯导静止放置，获取加速度计数据。用重力分量关系获得俯仰角θ和滚转角γ，利用地磁传感器获得导航坐标系下的偏航角完成惯导系统的初始对准。
[0121]
步骤2.1.1、利用加速度三轴分量估算俯仰角θ和滚转角γ：
[0122][0123]
其中，为载体坐标系下第k次加速度计的x、y、z三轴输出值。
[0124]
步骤2.1.2、结合获得的俯仰角θ和滚转角γ，将载体坐标系下的磁场强度转换至导航坐标系下
[0125][0126]
计算得到偏航角为：
[0127][0128]
步骤2.1.3、在静止条件下，三个欧拉角可以获得从载体坐标系到导航坐标系的方向余弦矩阵，完成惯导的初始对准。
[0129]
步骤2.2所述的设置初始零点，利用四元数以及陀螺仪的角增量，对惯导载体的姿态进行更新；
[0130]
步骤2.2.1、由陀螺仪的更新数据计算得到角增量δ：
[0131][0132]
ω
x
、ωy、ωz表示仪输出的三轴角速度标量值，tm为采样时间；
[0133]
步骤2.2.2、利用角增量对四元数进行更新：
[0134][0135]
式中，q
1|k+1
为k+1时刻四元数的q1值，q
1|k
为k时刻四元数的q1值，其他同理；步骤2.2.3、四元数进行归一化处理：
[0136][0137]
进行归一化使得四元数变为单位四元数，四个值的平方和为1；
[0138]
步骤2.2.4、根据单位四元数得到方向余弦矩阵，公式为：
[0139][0140]
步骤2.2.5、根据方向余弦矩阵得到欧拉角如下：
[0141][0142]
步骤2.2.6、根据上述矩阵信息再加上比力信息可以得到当前加速度三轴的分量，区分三轴分量可以对重力分量进行补偿，除去重力分量可以获得导航坐标系下的运动加速度，对运动加速度进行牛顿力学积分可以得到速度和位置，公式为：
[0143][0144]
其中，为根据补偿重力分量后的运动加速度计输出值。
[0145]
至此，可以获得炮塔、炮口的瞬时姿态、速度以及位置信息。
[0146]
进一步地，步骤3所述的在传递对准阶段，根据火炮身管特性进行外杆臂误差与挠曲形变的动态建模，对杆臂效应及挠曲形变进行补偿，具体如下：
[0147]
由于主、子惯导的安装位置不同，其各自敏感的信息存在一定的差别。为提高传递对准的性能，在进行信息匹配时，需要采取一定的措施对杆臂效应、挠曲变形等引起的误差进行补偿，尽可能真实的反映各个误差状态同观测量之间的关系。在杆臂效应中，外杆臂误差远远大于内杆臂误差影响，所以现忽略内杆臂误差，仅对外杆臂误差进行。
[0148]
(1)外杆臂误差补偿
[0149]
在实际的运用中，认为主、子惯导可能均未安装在载体的摇摆中心，甚至距摇摆中心有较长的距离，因此主、子惯导的加速度计敏感的加速度不同。由杆臂引起的误差在传递匹配过程中需要实时补偿。
[0150]
如图3所示，定义惯性坐标系为oixiyizi，载体坐标系为obxbybzb，并认为ob是载体的摇摆中心，在多数的研究中被认为是载体的重心，一般根据设计的载荷分布情况，求出重心位置，并且认为重心是固定的，且常假设主惯导安装位置与ob重合，子惯导的加速度计安装在载体坐标系中的固定点p。
[0151]
图4中，为载体坐标系原点的位置矢量，为p点相对于惯性坐标系原点的位置矢量，为p点相对于载体坐标系原点的位置矢量。
[0152]
显然它们有如下关系：
[0153][0154]
将上式两边对时间求微分可以得到：
[0155][0156]
再对上式求相对时间的微分可以得到：
[0157][0158]
根据矢量微分的相对微分原理，可以得到：
[0159][0160]
其中，表示p点相对于载体坐标系的运动线加速度。
[0161]
同理可以得到：
[0162][0163]
将上式组合可以得到p点相对于惯性坐标系的线加速度表达式：
[0164][0165]
在研究杆臂效应时一般认为载体是刚性结构，当存在挠性变形时需要采用其它的方法对挠性变形造成的误差加以补偿，设定p点相对于载体坐标系是固定的，所以：
[0166][0167][0168]
线加速度表达式可以进一步简化为：
[0169][0170]
理想情况下安装点应该在载体的摇摆中心，即
rp
＝0，这样就不存在杆臂效应。根据实际情况，在传递匹配过程中主、子惯导均未安装在载体的摇摆中心，杆臂效应引起的误差通常不可忽略。上式后面两项就是由于杆臂效应引起的惯导敏感到的杆臂加速度基本表达式。定义表示杆臂加速度，则杆臂效应误差的基本方程为：
[0171][0172]
由于：
[0173][0174]
因此，进一步整理可得：
[0175][0176]
其中，为切向加速度，为向心加速度。
[0177]
(2)挠曲变形的动态建模
[0178]
载体变形可以分为两类：一类是静态变形，它并不是绝对不变的，只是变化周期比较长而已，因此也称其为准静态变形；另一类就检测而言，变换较快，也即动态挠性变形。
[0179]
在近年对机体的研究中，通常把相应载体的动态结构变形视为马尔科夫过程，因为载体的动态变形是随机扰动干扰的随机变量，本发明中拟采用二阶马尔科夫过程作为载体动态变形的模型，并且认为各个轴的动态变形过程是独立的，即干扰形成的纵摇的噪声和干扰形成的横滚的噪声是相互独立的。动态变形的二阶马尔科夫过程如下：
[0180]
设定动态变形角为λ(t)，它是白噪声激励的二阶马尔科夫过程，另取动态变形角速率为
wf(t)
，即有：展开为：
[0181][0182]
动态变形角速率运动方程：
[0183][0184]
其中，β
x
＝2.146/τi(i＝x,y,z)，τi为相关时间，可以视具体的载体情况而定，wi(t)一般认为是具有一定方差的白噪声，如果是有色噪声，则需要进行白化处理，其方差满足：足：表示三个弹性变形角λ
x
,λy,λz的方差。
[0185]
进一步地，步骤4所述的将mems子惯导解算出的速度、姿态相对于高精度主惯导解算出的速度、姿态偏差量作为子惯导的量测信息，建立以速度加姿态匹配方式的传递对准卡尔曼滤波器，具体如下：
[0186]
(1)速度加姿态匹配的误差模型
[0187]
子惯导中各项速度矢量关系为：
[0188][0189]
由惯导基本方程可知：
[0190][0191]
将速度矢量公式左右两边对时间求导，代入上式，并向导航坐标系投影可知：
[0192][0193]
其中：其中：为子惯导系统的姿态矩阵。
[0194]
进一步整理可得：
[0195][0196]
子惯导相对i系的加速度为：
[0197][0198]
进一步可知：
[0199][0200]
同理可知主惯导中：
[0201][0202]
上述两式中都含有项，不过是在不同的主、子惯导惯导系统中求得的，分别表示不同的加速度信息：
[0203]
子惯导中：
[0204][0205]
主惯导中：
[0206][0207]
组合上式可知：
[0208][0209][0210]
上述两式作差，令可知：
[0211][0212]
且：
[0213][0214][0215]
进一步可得：
[0216][0217]
其中：是子惯导在s中感受到的比力
[0218][0219]
其中：
[0220]
子惯导敏感到的主惯导刚性转动引起的加速度；
[0221]
子惯导敏感到的挠性加速度；
[0222]
子惯导加速度计误差。
[0223]
进一步得到：
[0224]
[0225]
忽略高阶小量并且考虑整理得到：
[0226][0227]
假设杆臂效应和挠性运动引起的误差项已经得到补偿，整理可得：
[0228][0229]
用子惯导的姿态矩阵代替主惯导的姿态矩阵最终得到主、子惯导的速度差微分方程：
[0230][0231]
设定静态变形角(包括安装误差)为(简记μ)，动态变形角(即挠曲变形)为(简记θ)，此时由于存在静态变形角及动态变形角，应改为：
[0232][0233]
式中，[μ
×
]与[θ
×
]分别表示静态变形角与动态变形角的反对称阵。上式忽略二阶小项，可近似为：
[0234][0235]
进一步得考虑静态变形角与动态变形角时的姿态差方程为：
[0236][0237]
整理后得到姿态量测方程为：
[0238][0239]
式中，φ表示主、子惯导解算的姿态偏差量，z
φ
表示主、子惯导去掉静态变形角与动态变形角影响的姿态差。
[0240]
(2)速度加姿态匹配方式的卡尔曼滤波模型
[0241]
卡尔曼滤波算法本质上利用从初始时刻到当前时刻的所有测量信息，采用迭代方法，不需要保存之前的测量量，同时利用状态和量测方程。在利用量测方程的基础上，将状态方程也纳入滤波算法中，同时利用状态自身的变化规律和测量量，最大程度地提高炮口惯导估计精度。
[0242]
kalman滤波主要采用离散递推的表达形式，它的状态空间模型如下：
[0243][0244]
其中，xk为系统在k时刻的状态向量，也就是要求的被估计状态量；φ
k,k-1
为k-1时刻到k时刻的状态转移矩阵；γ
k-1
为系统噪声分配矩阵，它表示由k-1时刻到k时刻的各个系统噪声分别影响k时刻各个状态量的程度；w
k-1
代表k-1时刻的系统噪声矢量；zk表示的是系统在k时刻的量测矢量；hk为k时刻的系统量测矩阵；vk表示k时刻的量测噪声矢量。同时，根据卡尔曼滤波的要求，wk和vk为相互独立的零均值高斯白噪声矢量序列，满足：
[0245][0246]
上式为卡尔曼滤波过程中噪声需要满足的条件，同时由于qk为系统的噪声方差阵，而系统的某个状态量可能并不存在噪声，因此qk为非负定矩阵；rk为量测噪声方差阵，而每个量测量都含有噪声，因此rk为正定矩阵。
[0247]
卡尔曼滤波算法可以由以下几个方程进行表示：
[0248]
(1)状态一步预测
[0249][0250]
(2)状态一步预测均方误差阵
[0251][0252]
(3)滤波增益
[0253][0254]
(4)状态估计
[0255][0256]
(5)状态估计均方误差阵
[0257]
pk＝(i-k
khk
)p
k|k-1
ꢀꢀꢀ
(75)
[0258]
设定静态变形角(包括安装误差)为(简记μ)，动态变形角(即挠曲变形)为(简记θ)，动态变形角速率(简记ω)。动态变形角θ采用二阶markov模型：
[0259][0260]
即：
[0261][0262]
其中，var(w)＝q＝4β3σ2，σ为动态变形角θ的方差，τ为相关时间，w,q为激励噪声及其方差。
[0263]
选状态量其中ε为陀螺仪零漂,为加速度计零偏，状态方程展开为
[0264][0265]
量测方程z＝hx+v展开为
[0266][0267]
那么，系统转换矩阵为
[0268][0269]
观测矩阵为
[0270][0271]
其中，为子惯导确立的姿态阵。
[0272]
进一步地，步骤5所述在导航更新阶段，子惯导利用传递对准结束时刻的对准误差，作为导航解算的初始误差，进行导航状态更新，完成炮口的振动测量，具体如下：
[0273]
传递对准结束后，子惯导能从状态估计量中获取当前陀螺仪的零漂及加速度计的零偏，以及主、子惯导的失准角等。子惯导利用获取的三个欧拉角得到从载体坐标系到导航坐标系的方向余弦矩阵，结合步骤2中的姿态更新算法可以实现子惯导的实时位置姿态测量。
[0274]
由于子惯导与炮口固连，所以子惯导的实时位置姿态信息即为炮口的实时振动情况，至此完成炮口的振动测量。
[0275]
综上所述，本发明在炮塔处安装高精度主惯导，在身管炮口处安装mems子惯导的接触式测量方式，相对传统光学等非接触式测量方式，安装方便、稳定性高，避免了传统光学方式的易受干扰、安装困难、仪器校准繁琐等缺点；在初始对准阶段，采用高精度主惯导对mems子惯导进行传递对准，属于惯导精对准范畴，可以保证mems子惯导有准确的初始惯导信息，提高了炮口振动测量的可靠性和准确性；在导航更新阶段，子惯导利用传递对准结束时刻的对准误差，作为导航解算的初始误差，利用姿态更新算法不断迭代解算，可以保证mems子惯导的位姿输出准确性，提高了炮口振动测量方法的实时性和可靠性。

技术特征：
1.一种基于主子惯导结合的行进间炮口振动测量方法，其特征在于，包括以下步骤：步骤1、根据火炮平台，设计主、子惯导安装位置，包括主惯导和mems子惯导；步骤2、根据火炮行进间主惯导、mems子惯导传感器输出，设计基于陀螺仪角增量更新的姿态更新算法；步骤3、在传递对准阶段，根据火炮身管特性进行外杆臂误差与挠曲形变的动态建模，对杆臂效应及挠曲形变进行补偿；步骤4、将mems子惯导解算出的速度、姿态，相对于主惯导解算出的速度、姿态偏差量作为mems子惯导的量测信息，建立以速度加姿态匹配方式的传递对准卡尔曼滤波器；步骤5、在导航更新阶段，mems子惯导利用传递对准结束时刻的对准误差，作为导航解算的初始误差，进行导航状态更新，完成炮口的振动测量。2.根据权利要求1所述的基于主子惯导结合的行进间炮口振动测量方法，其特征在于，步骤1中，根据火炮平台，设计主、子惯导安装位置，包括主惯导和mems子惯导，具体如下：所述主惯导固连于炮塔，保持水平状态，炮塔进行航向运动带动主惯导一同运动；所述mems子惯导固连于火炮身管前端，用于测量火炮身管前端炮口的位姿。3.根据权利要求1所述的基于主子惯导结合的行进间炮口振动测量方法，其特征在于，步骤2中根据火炮行进间主惯导、mems子惯导传感器输出，设计基于陀螺仪角增量更新的姿态更新算法，具体如下：步骤2.1、令惯导静止放置，获取加速度计数据，用重力分量关系获得俯仰角θ和滚转角γ，利用地磁传感器获得导航坐标系下的偏航角完成惯导系统的初始对准；步骤2.2、设置初始零点，利用四元数以及陀螺仪的角增量，对惯导载体的姿态进行更新。4.根据权利要求3所述的基于主子惯导结合的行进间炮口振动测量方法，其特征在于，所述步骤2.1具体如下：步骤2.1.1、利用加速度三轴分量估算俯仰角θ和滚转角γ：步骤2.1.1、利用加速度三轴分量估算俯仰角θ和滚转角γ：其中，为载体坐标系下第k次加速度计的x、y、z三轴输出值；步骤2.1.2、利用地磁传感器获得载体坐标系下的磁场强度结合获得的俯仰角θ和滚转角γ，经过转换得导航坐标系下得磁场强度进一步计算得到偏航角为：其中，θ
k
、γ
k
分别表示第k次解算出的俯仰角、滚转角；
其中，分别表示导航坐标系下x、y、z三轴对应的磁场强度；步骤2.1.3、在静止条件下，利用三个欧拉角获得从载体坐标系到导航坐标系的方向余弦矩阵，完成惯导的初始对准。5.根据权利要求3所述的基于主子惯导结合的行进间炮口振动测量方法，其特征在于，所述步骤2.2，具体如下：步骤2.2.1、由陀螺仪的更新数据计算得到角增量δ：式中，ω
x
、ω
y
、ω
z
表示仪输出的三轴角速度标量值，t
m
为采样时间；步骤2.2.2、利用角增量对四元数进行更新：式中，q
1|k+1
、q
2|k+1
、q
3|k+1
、q
4|k+1
分别为k+1时刻四元数的q1、q2、q3、q4值，q
1|k
、q
2|k
、q
3|k
、q
4|k
分别为k时刻四元数的q1、q2、q3、q4值；步骤2.2.3、四元数进行归一化处理：其中，q1'
|k+1
、q'
2|k+1
、q'
3|k+1
、q'
4|k+1
分别为q
1|k+1
、q
2|k+1
、q
3|k+1
、q
4|k+1
归一化处理后的值；进行归一化使得四元数变为单位四元数，四个值的平方和为1；步骤2.2.4、根据单位四元数得到方向余弦矩阵，公式为：步骤2.2.5、根据方向余弦矩阵得到欧拉角如下：
步骤2.2.6、根据方向余弦矩阵信息和比力信息，得到当前加速度三轴的分量，区分三轴分量能够对重力分量进行补偿，除去重力分量获得导航坐标系下的运动加速度，对运动加速度进行牛顿力学积分得到速度和位置，公式为：其中，为根据补偿重力分量后的运动加速度计输出值，表示从k-1时刻到k时刻，从载体坐标系b系到导航坐标系n系的坐标变换矩阵。6.根据权力要求1所述的基于主子惯导结合的行进间炮口振动测量方法，其特征在于，步骤3中的外杆臂误差补偿，具体如下：定义惯性坐标系为o
i
x
i
y
i
z
i
，载体坐标系为o
b
x
b
y
b
z
b
，o
b
是载体的摇摆中心，子惯导的加速度计安装在载体坐标系中的固定点p，为载体坐标系原点的位置矢量，为p点相对于惯性坐标系原点的位置矢量，为p点相对于载体坐标系原点的位置矢量，则p点相对于惯性坐标系的线加速度表达式为：上式后面两项就是由于杆臂效应引起的惯导敏感到的杆臂加速度基本表达式；其中，表示陀螺仪输出的在惯性系下的角速度；定义表示杆臂加速度，则杆臂效应误差的基本方程为：其中，表示杆臂加速度，为切向加速度，为向心加速度。7.根据权利要求1所述的基于主子惯导结合的行进间炮口振动测量方法，其特征在于，步骤3中所述的挠曲变形的动态建模，具体如下：设定动态变形角为λ(t)，动态变形角速率为w
f
(t)，即有：展开为：
其中，分别表示x，y，z三轴动态变形角速率，分别表示x，y，z三轴动态变形角速率，分别为x，y，z三轴动态变形角的一阶导数；动态变形角速率运动方程：其中，令i＝x,y,z，β
i
＝2.146/τ
i
，τ
i
为相关时间；w
i
(t)是白噪声，如果是有色噪声，则需要进行白化处理，白噪声w
i
(t)的方差满足表示三个弹性变形角λ
x
,λ
y
,λ
z
的方差。8.根据权利要求1所述的基于主子惯导结合的行进间炮口振动测量方法，其特征在于，步骤4中的mems子惯导解算出的速度、姿态，相对于主惯导解算出的速度、姿态偏差量，具体如下：主、子惯导的速度差微分方程为：式中，表示主子惯导的速度差，为子惯导系统的姿态矩阵，φ
×
表示主、子惯导的初始安装误差角对应的反对称矩阵，ξ表示主、子惯导的安装误差角估计值；表示主惯导在子惯导坐标系下感受到的比力；表示地球在导航坐标系下的自转角速率；表示子惯导加速度计误差；姿态量测方程为：式中，φ表示主、子惯导解算的姿态偏差量；z
φ
表示主、子惯导去掉静态变形角与动态变形角影响的姿态差；μ表示静态变形角，包括安装误差；θ表示动态变形角，即挠曲变形。9.根据权利要求1所述的基于主子惯导结合的行进间炮口振动测量方法，其特征在于，步骤4中所述的速度加姿态匹配方式的卡尔曼滤波模型，具体如下：(1)状态一步预测其中，为k时刻状态的一步预测；φ
k,k-1
为k-1时刻到k时刻状态一步转移矩阵；为k-1时刻状态的最优估计；(2)状态一步预测均方误差阵其中，p
k|k-1
为状态一步预测均方误差阵；p
k-1
为k-1时刻状态最优估计的均方误差阵；(3)滤波增益其中，系数矩阵k
k
称为滤波增益；为k时刻的系统量测矩阵的转置；p
k
为k时刻状态最优估计的均方误差阵；
(4)状态估计其中，为k时刻状态的最优估计；z
k
为k时刻的量测；h
k
为k时刻的系统量测矩阵；(5)状态估计均方误差阵p
k
＝(i-k
k
h
k
)p
k|k-1
(21)。10.根据权利要求1所述的基于主子惯导结合的行进间炮口振动测量方法，其特征在于，所述步骤5具体如下：传递对准结束后，子惯导从状态估计量中获取当前陀螺仪的零漂及加速度计的零偏，以及主、子惯导的失准角；子惯导利用获取的三个欧拉角得到从载体坐标系到导航坐标系的方向余弦矩阵，结合步骤2中的姿态更新算法，实现子惯导的实时位置姿态测量；由于子惯导与炮口固连，所以子惯导的实时位置姿态信息即为炮口的实时振动情况。

技术总结
本发明公开了一种基于主子惯导结合的行进间炮口振动测量方法，该方法为：在炮塔上安装抗高过载高精度惯导，在身管的炮口处安装低精度的MEMS惯导；在初始阶段，炮塔处的高精度主惯导与炮口处的MEMS子惯导完成传递对准；在传递对准阶段，根据火炮身管特性进行外杆臂误差与挠曲形变的动态建模，对杆臂效应及挠曲形变进行补偿；建立以速度加姿态匹配方式的传递对准卡尔曼滤波器；在导航更新阶段，子惯导利用传递对准结束时刻的对准误差，作为导航解算的初始误差，进行导航状态更新，完成炮口的振动测量。本发明具有安装简便、稳定性高、抗干扰能力强、准确性高的优点。准确性高的优点。准确性高的优点。


技术研发人员：朱建良 田超亚 王立雅 薄煜明 吴盘龙 王军
受保护的技术使用者：南京理工大学
技术研发日：2022.12.08
技术公布日：2023/7/19