1 引言
脉冲修正弹药是当前世界各国智能弹药发展的一个重要方向。由于采用了弹道闭环控制，脉冲修正弹药的命中精度大幅度提高。使得其不仅具备攻击点-面目标的能力，甚至具备攻击点目标的能力，大大拓展了作战能力和综合作战效能[1]。
其中，弹丸姿态解算直接影响到脉冲修正弹的控制精度，是脉冲修正弹技术体系中一项关键的技术。随着地磁理论的不断发展和完善，技术也在不断进步，基于磁传感器测量姿态角的方法已得到广泛运用，该方法具有高可靠性、误差也不会随时间累积的优点[2]。
国内外众多专家和学者对利用地磁信息进行制导弹药姿态解算的方法进行了研究，取得了一定的成果。本文在这些研究成果的基础上[3-10]，提出了一套利用发射点地磁信息和弹上传感器获取的地磁信息来解算弹丸滚转角的算法，并对解算精度进行了验证。
2 数学模型
2.1 地磁场要素
通常用地磁感应强度T及其分量描述地磁场的特征[11]。以发射点O为原点建立坐标系，X轴沿着地理子午线，向北为正；Y轴沿着纬度方向，向东为正；Z轴垂直于平面，向下为正。并将地磁感应强度T投影在此坐标系下，得到各个分量。如图1所示。
图1 地磁场感应强度示意图
Fig.1 Geomagnetic field elements
图1中各个参数意义如下：T为磁感应强度；HX为北向强度，磁感应强度T在X轴上的投影；HY为东向强度，磁感应强度T在Y轴上的投影；HZ为垂直强度，磁感应强度T在Z轴上的投影；H为水平强度，磁感应强度T在水平面上的投影；D为磁偏角，磁子午面与地理子午面的夹角，且规定东偏为正，反之为负；I为磁倾角，磁感应强度T与水平面之间的夹角，且在北半球指向地平线之下，此时磁倾角为正；反之在南半球指向上，磁倾角为负。
以上7个量统称为地磁要素，他们之间的关系如式(1)：
(1)
2.2 坐标转换
设地磁场的磁感应强度T在地面坐标系Oxyz三个轴上的分量为：Hx、Hy、Hz。如图2所示。
图2 磁感应强度在地面坐标系的分量图
Fig.2 Component diagram of magnetic induction intensity in ground coordinate system
定义地面坐标系的Ox轴偏离地理北方向的角度η为射向角。射向角以北向为基准，北偏东为正。则磁感应强度矢量T在地面坐标系Oxyz三个轴上的分量可表示为
若给定弹丸的俯仰角θ和偏航角φ，磁感应强度矢量T在准弹体坐标系Ox4y4z4三个轴上的分量可表示为
2.3 俯仰角、偏航角更新方法
弹丸飞行过程中，需要实时获取俯仰角和偏航角数据。若采用雷达跟踪弹丸，可以得到弹丸在雷达坐标系下的斜距R、高低角α和方位角β。雷达坐标系以发射点为原点，发射点与弹丸连线的长度为斜距R，连线与水平面的夹角为高低角α，连线在水平面上方，则高低角α为正，反之为负。连线在水平面上的投影与北向的夹角为方位角β，方位角β以北向为基准，北偏东为正。
以发射点为原点，北向为X轴，东向为Y轴，垂直向上为Z轴，建立直角坐标系。将雷达坐标系下的坐标参数转化为直角坐标系上的坐标参数。如式(2)所示：
(2)
弹丸速度矢量在直角坐标系上表示为
(3)
设η为射向角，射向角以北向为基准，北偏东为正。则弹丸速度矢量在地面坐标系上表示为
假设弹丸速度矢量和弹轴重合，俯仰角和偏航角可用下列公式计算：
(4)
(5)
在某些工程应用条件下，弹丸速度矢量和弹轴之间夹角较小，利用式(4)和式(5)计算出来的弹丸俯仰角和偏航角在一定误差允许范围内可以被接受。
2.4 滚转角解算
在弹体横截面内，准弹体坐标系和弹体坐标系上的磁感应强度分量如图3所示。
图3 磁场感应强度分解示意图
Fig.3 Exploded view of magnetic field
定义弹体横截面内磁感应强度矢量H与y4轴的夹角φ0为滚转基准角。计算公式为
(6)
定义弹体横截面内磁感应强度矢量H与y1轴的夹角为α0。计算公式为
(7)
则滚转角γ为
(8)
3 数值仿真
利用Matlab搭建数值仿真模型，进行该滚转角解算算法的数值仿真验证。采用变步长ode45解算器，最大步长限制为0.01 s，仿真时间设置为0～10 s，相对误差限设置为0.001。
3.1 输入参数
在该阶段的工作中，只对该算法本身解算滚转角的可行性进行验证，因此可将俯仰角和偏航角作为定值，给定射向角，给定作者工作单位实验室附近的地磁要素HX、HY、HZ，各项输入参数值如表1所示。
表1 输入参数值
Table 1 Value of input parameter
参数值参数值HX/nT31 727射向角/(°)90HY/nT-2 218 偏航角/(°)0HZ/nT42 206俯仰角/(°)0
设置初始滚转角为0，转速为1 r/s，转向为左旋，则磁传感器磁轴输出值Hy1、Hz1的理论值为[12]
Hy1=52 800cos(2πt+3.786)
Hz1=52 800sin(2πt+3.786)
3.2 结果及分析
在3.1节给出的条件下，弹丸滚转角的准确值可用下式计算。公式形式为
γ=(-2π×t)%(2π)
(9)
式(9)中，%为取模运算。
将弹丸滚转角的准确值曲线与算法解算得到的计算值曲线绘制在同一坐标下，如图4所示。准确值结果用黑线表示，离散点标记符号为黑色菱形；计算值结果用蓝线表示，离散点标记符号为蓝色正方形。从图4可知，两条曲线完全重合，说明算法解算得到的滚转角计算值和滚转角的准确值完全一致，验证了算法的可行性。
图4 滚转角准确值和算法计算值曲线
Fig.4 Comparison of the accurate value of the roll angle and the calculated value of the algorithm
4 半实物仿真
利用手摇无磁转台加LABVIEW开发环境组建的半实物仿真平台进行该滚转角解算算法的半实物仿真验证。相比数值仿真利用地磁信息理论值进行滚转角解算的方式，半实物仿真验证能够通过地磁传感器获取的真实地磁信息进行滚转角解算，对算法在实际应用情况下的解算精度进行验证。
4.1 试验器材
所需试验器材：手摇无磁转台、三轴磁通门地磁传感器、数据采集装置。试验器材如图5所示。
图5 试验器材图片
Fig.5 Test equipment
4.2 试验方法
试验方法和步骤如下：
1) 将三轴地磁传感器装配进钛合金弹体(选用钛合金弹体是想要消除弹体对地磁信号的影响)并灌封；
2) 将数学公式运用LABVIEW编制成计算程序并制作传感器输出数据和滚转角解算结果的显示界面；
3) 通过螺纹连接弹体和手摇无磁转台，将地磁传感器的导线和传感器数据采集装置连接，将传感器y轴处于铅锤面时的转盘刻度值当成0°；
4) 调整转台俯仰角和射向角并固定，如图6所示。
图6 转台示意图
Fig.6 Schematic diagram of turntable placement
5) 输入初始参数，在该阶段的工作中，只对该算法解算滚转角的可行性和精度进行验证，因此可将俯仰角和偏航角作为定值。输入初始参数值如表2所示。
表2 输入初始参数值
Table 2 Value of input parameters
参数值参数值HX/nT31 727射向角/(°)90HY/nT-2 218 偏航角/(°)0HZ/nT42 206俯仰角/(°)0、30、60
6) 启动LABVIEW数据处理平台，进行传感器输出的采集。待数据处理平台显示曲线稳定之后，转动手摇转台，通过数据处理平台界面读取滚转角的输出值；
7) 连续转动转台，观察滚转角变化曲线，通过与数值仿真结果的对比验证算法在实物平台上运行并且进行滚转角解算的可行性。
8) 离散的转动转台，记录数据处理平台输出的滚转角和手摇无磁转台上的刻度值，并将刻度值作为滚转角的准确值。通过对比，验证算法解算精度。
4.3 试验结果
4.3.1 干扰影响分析
由于实验室电磁干扰和弹体的影响，传感器测量的地磁信号值和理论值必然存在偏差。在实验室环境下和室外环境下，将传感器获得的地磁信号值和理论值绘制在同一坐标下以研究电磁干扰和弹体对传感器测量值的影响[13-15]，如图7和图8所示。
图7 实验室地磁信号曲线(90°射向，0°射角)
Fig.7 Indoor laboratory (90 degree firing direction, 0 degree firing angle)
图8 室外环境地磁信号曲线(90°射向，0°射角)
Fig.8 Outdoor environment (90° shooting direction, 0° shooting angle)
观察两图发现：图7中测量值曲线和理论值曲线相差较大，图8中测量值曲线和理论值曲线相差较小。并且从图7和图8的对比中可知，选用钛合金弹体基本达到了预期目的，在室外电磁干扰较弱的情况下，传感器测量值和理论值非常接近。而在室内环境下，由于电磁干扰较强，传感器测量值和理论值的偏差则更大一些。
4.3.2 可行性验证
连续向左转动转台6圈，在LABVIEW数据处理平台界面显示了滚转角变化曲线，如图9。上方为地磁传感器输出值的变化曲线，弹体y轴和z轴方向的地磁强度按正弦规律变化，弹体x轴方向的地磁强度基本不变。下方为算法解算得到的滚转角变化曲线，波形为直角三角形，呈现循环下降的态势。和图4比较可知，半实物仿真结果和数值计算结果相吻合，验证了算法的可行性。
图9 滚转角变化曲线
Fig.9 Roll angle change curve
4.3.3 精度验证
分别将弹丸俯仰角设置为0°、30°、60°，然后将传感器y轴转至铅锤面内，并将此时的转盘刻度值当做0°。以此为基准，旋转1圈，每旋转30°记录1次算法解算结果。数据统计结果如表3所示。
表3 不同弹丸俯仰角算法解算数据统计结果
Table 3 Data statistics
转盘刻度值/(°)弹丸俯仰角0°30°60°00.571 70.061 947359.0723031.069 529.658 328.538 96060.680 759.289 958.170 69091.577 190.326 589.035 2120123.052122.121120.839150153.433152.401151.965180182.521182.369182.678210212.161211.790212.882240242.355242.566242.279270272.616272.778272.593300301.989301.407300.024330330.872330.291329.4823600.802 8359.342358.567
从表3可知，算法解算得到的滚转角精度在3°以内。若是在室外环境下，外界电磁干扰更弱，同时给定的输入参数更靠近准确值，那么预期将会获得更好的解算精度。
5 结论
1) 设计了一套非常简洁的可用于解算弹丸滚转角的解算算法，并进行了数值仿真验证和半实物仿真验证。该算法配合磁通门传感器解算弹丸滚转角具有良好的精度。
2) 利用弹道测量数据进行弹丸俯仰角、偏航角的近似计算，在弹丸速度矢量和弹轴重合或者偏差较小的理想条件下，具有一定的可行性。
3) 通过试验证明了钛合金弹体对于磁场的影响微乎其微。在实际应用过程中无需再通过标定消除弹体对弹内传感器测量值的干扰。
[1] 陈胜政,杨波,杨栓虎,等.脉冲修正弹药技术研究[J].兵器装备工程学报,2017(02):60-64.
Chen S Z,Yang B,Yang S H,et al.Research on impulse correction munitions technology[J].Journal of Ordnance Equipment Engineering,2017(02):60-64.
[2] 任先贞.基于磁传感器组合的旋转弹姿态测试研究[D].太原：中北大学,2017.
Ren X Z.Research on attitude measurement of rotating projectile based on magnetic sensor[D].Taiyuan:North University Of China,2017.
[3] 马秀丽.基于磁阻/陀螺组合测量弹体飞行姿态的数据融合方法研究[D].南京：南京理工大学,2012.
Ma X L.Investigation on data fusion method based on magnetometer/gyro combine measurement of rojectile’s flying attitudes[D].Nanjing:Nanjing University of sciences & Technology,2012.
[4] 张莺莺,张晓明,薛羽阳.地磁旋转弹滚转角及角速率的估计方法[J].中国测试,2021,47(04):38-43.
Zhang Y Y,Zhang X M,Xue Y Y.An estimation method of roll angle and rate of spinning projectile based on geomagnetic information[J].China Measurement & Test,2021,47(04):38-43.
[5] Liu H,Dong H B,GeJ,et al.Research on a secondary tuning algorithm based on SVD & STFT for FID signal[J].Measurement Science and Technology,2016,27(10):0957-0233.
[6] Yan X L,Chen G G,Tian X L.Two-step adaptive augmented unscented kalman filter for roll angles of spinning missiles based on magnetometer measurements[J].Measurement & Control,2018,51(05):73-82.
[7] Li L B,Li M X,Jiang L X,et al.Angular rate estimation and damping control of satellite with magnetometer data[J].Optik,2019,180(11):1049-1055.
[8] Maley J M.Efficient attitude estimation for a spin-stabilized projectile[J].Journal of Guidance Control and Dynamics,2016,36(02):339-350.
[9] Ma H L,Xu S J.Magnetometer-only attitude and angular velocity filtering estimation for attitude changing spacecraft[J].Acta Astronautica,2014,102(05):89-102.
[10] Zhou Y K,Zhang X Y,XiaoW.Spinning projectile’s angular measurement using crest and through data of a geomagnetic sensor[J].Measurement Science and Technology,2018,29(09):1-9.
[11] 徐华威.MEMS惯性/地磁组合测姿技术研究[D].太原：中北大学,2018.
Xu H W.Research on MEMS inertial/geomagnetic combination posture measurement technology[D].Taiyuan:North University of China,2018.
[11] Zhang C.Progress in time synchronization technology for wireless sensor networks[J].Journal of Chongqing Technology and Business University(Natural Science Edition),2019,36(06):88-94.
[13] Klingbeil L，Eling C，Florian Z，et al.Magnetic field sensor calibration for attitude determination[J].Journal of Applied Geodesy，2014，8(02):97-108.
[14] Renaudin V，Afzal M H，Lachapelle G.Complete triaxis magnetometer calibration in the magnetic domain[J].Journal of Sensors，2010(01):1-10.
[15] 李冰,雷泷杰,陈超.基于椭圆拟合的双轴磁传感器标定方法[J].探测与控制学报,2020,42(03):20-23.
Li B,Lei L J,Chen C.A calibration method of biaxial magnetic sensor based on ellipse fitting[J].Journal of Sensors，2020,42(03):20-23.