星敏与磁力仪间安装矩阵的一种地面标定方法

刘逸康，任顺清，张高雄，张 祎

(1.哈尔滨工业大学 空间控制与惯性技术研究中心，哈尔滨 150080；2.上海卫星工程研究所，上海 200240)

地磁场提供了丰富的地磁特征[1]，利用卫星地磁测量，能够快速地测量全球的地磁场矢量信息.在此基础上完成地磁场模型库的建立与更新，并广泛应用于地磁导航、地震预报等技术领域[2-3].地磁测量的精度受到测量用磁力仪和星敏感器等的仪器安装误差以及坐标系间传递精度的影响[4-7]，因此提高卫星地磁测量系统内传感器的精度和磁力仪与星敏感器之间安装矩阵的测量精度具有重要的实用价值.

目前对磁力仪及星敏感器之间安装矩阵的标定方法研究相对较少，张艺腾等[8]对航空地磁测量中的磁力仪与姿态仪的联合标定进行了研究，假设磁力仪三轴与姿态仪三轴之间为欧拉转换的关系，通过多次旋转，由最小二乘通过扫描迭代求出最优的旋转欧拉角，最终得到传感器间的安装矩阵.孙闯等[9]对安装矩阵提出一种户外标定的方法，该方法借助一个标准六面体对伸展杆系统进行多个方向的旋转，测量不同姿态下星敏与磁力仪的输出，并根据测出的星敏感器与磁力仪敏感轴坐标系相对地理坐标系的姿态辨识出星敏与三轴磁力仪间的安装矩阵，但该方法达不到精度要求.王凯强等[10]提出了在实验室环境下利用三轴无磁转台提供姿态信息以及恒星模拟器模拟室外观星环境，通过最小二乘辨识出安装矩阵的方法，但用三轴无磁转台在户外标定既不方便也不经济.

本文提出了一种利用单轴立式无磁转台对磁力仪及星敏感器之间安装矩阵在户外进行地面标定的方法，通过无磁转台提供精确的角位置使磁力仪敏感到不同的地磁场分量的激励，不仅辨识出磁力仪的非正交误差，并在此之上完成安装矩阵的标定.

1 地面标定系统介绍

针对三轴磁力仪以及3个星敏感器，如果各个敏感轴的方向矢量能够统一在地理坐标系下表示，即可确定他们间的安装矩阵.当星敏感器在地面观星时，可得到星敏光轴在惯性空间的姿态，再利用当地经纬度信息以及测量时刻的格林尼治恒星时推导出星敏光轴在地理坐标系下的表示.在单个位置下三轴磁力仪测量地磁场无法直接解算出其敏感轴的姿态，本文借助单轴无磁转台以及安装夹具，在多个位置进行测量，推导出初始位置下的磁力仪敏感轴在地理坐标系下的表示.最终通过坐标系转换等数学方法，解算出精确的安装矩阵.根据以上原理，整个地面标定系统如图1所示.

图1 地面标定系统示意

磁力仪和星敏感器由刚性杆连接成固联体，安装在单轴立式无磁转台上，构成了无磁实验平台.地面监测系统包括旋转轴监测系统以及地磁场监测磁强计.无磁转台由于必须采用无磁材料等的原因存在着较大的倾角回转误差，使用旋转轴监测装置对转台的回转误差进行实时监测，并对输出进行补偿.标定时需要提供精确的当地磁场矢量信息，故采用监测用磁强计测量当地磁场的变化.

2 坐标系的建立与误差传递

根据地面标定系统组成，依次建立以下坐标系.地心惯性坐标系OiXiYiZi，Oi为地球的地心，Xi轴和Yi轴在地球的赤道平面内，Xi轴指向春分点，Zi轴指向地球北极，Yi轴的方向由右手定则来决定.

地理坐标系OXYZ，即东北天坐标系g.如果已知格林尼治恒星时t，测量地点的经度λ以及纬度φ，则可求出惯性系i与地理系g间的转换矩阵为

(1)

星敏感器光轴坐标系OpXpYpZp.由3个星敏光轴组成的非正交坐标系.OpXp沿着X星敏光轴，OpYp，OpZp分别沿着其余星敏光轴.

转台轴套坐标系O1X1Y1Z1固联在轴套上.设地理坐标系为单轴无磁转台的参考坐标系，则轴套坐标系相对地理坐标系的变换矩阵为[11]

(2)

式中Δθx0、Δθy0称为主轴轴线的铅垂度，可在试验前通过水平仪精确测量.

主轴坐标系O2X2Y2Z2，主轴坐标系相对于轴套坐标系的姿态矩阵为[11]

rot(z1,γ),

(3)

式中Δθx1(γ)、Δθy1(γ)为转台的倾角回转误差，在试验时由旋转轴监视系统进行实时监控.

安装基面坐标系O3X3Y3Z3，为考虑转台安装基面相对于轴系几何轴线的垂直度Δαx2，Δαy2形成，安装基面坐标系相对主轴坐标系的姿态矩阵为[11]

(4)

磁力仪测量坐标系O4X4Y4Z4.为正交坐标系，磁力仪位于初始位置时，O4X4轴与安装基面O3X3轴相差了一个绕O3Z3轴的初始安装角β，磁力仪测量坐标系相对安装基面坐标系的姿态矩阵为

(5)

由于磁力仪3个敏感轴存在着三轴非正交误差，初始安装磁力仪时，令磁力仪X轴Y轴处于水平状态，此时磁力仪敏感轴O4Z4与磁力仪敏感轴O5Z5平行，现称为第1种安装方式.设磁力仪敏感轴矢量在磁力仪测量坐标系O4X4Y4Z4的表示为

(6)

式中Δθyx、Δθzx、Δθzy分别为3个敏感轴间的垂直度.若3个参数为正，则表示磁力仪敏感轴之间夹角大于90°.由于加工水平与安装工艺的限制，以及磁力仪的内部存在剩磁，磁力仪输出模型中还存在着标度因子kx,ky,kz和零偏bx,by,bz，可在试验之前通过标定方法求出.设初始安装位置如图1所示.结合式(2)～(6)，此时磁力仪的输出表示为[12]

(7)

3 磁力仪误差标定方法

称图1中固联体姿态为第1种安装方式下的初始位置，当单轴无磁转台旋转角度为γ时，采集三轴磁力仪的输出.展开式(7)并忽略二阶小量，有：

Rx1=[kxBx(cosβ+Δθyxsinβ)+kxBy(sinβ-

Δθyxcosβ)+kxBz(-Δθy0cosβ-Δθy1cosβ+

Δθx0sinβ+Δθx1sinβ)]cosγ+[kxBx(-sinβ+

Δθyxcosβ)+kxBy(cosβ+Δθyxsinβ)+

kxBz(Δθx0cosβ+Δθx1cosβ+Δθy0sinβ+

Δθy1sinβ)]sinγ+kxBz(-Δθzx-Δαy2cosβ+

Δαx2sinβ)+bx,

(8)

Ry1=[kyBx(-sinβ)+kyBy(cosβ)+

kyBz(Δθx0cosβ+Δθx1cosβ+Δθy0sinβ+

Δθy1sinβ)]cosγ+[-kyBxcosβ-kyBysinβ+

kyBz(Δθy0cosβ+Δθy1cosβ-Δθx0sinβ-

Δθx1sinβ)]sinγ+kyBz(-Δθzy+Δαx2cosβ+

Δαy2sinβ)+by,

(9)

Rz1=[Δαy2kzBx-Δαx2kzBy]cosγ+[Δαx2kzBx+

Δαy2kzBy]sinγ+kzBx(Δθy0+Δθy1)-

kzBy(Δθx0+Δθx1)+kzBz+bz.

(10)

式(8)～(10)中存在的辨识参数有β，Δθyx，Δθzx，Δθzy，Δαx2和Δαy2.考虑使用谐波分析辨识，此时一周内采样点个数应大于辨识量的两倍.令γ=2πi/n(i=0,1,2,…,n-1)，n=24，即单轴转台等间隔旋转24次，每次转台停稳并采集完磁力仪输出后进行下次旋转.

首先应辨识初始旋转角β，因只有在β已知的时才能获取磁力仪敏感轴在地理坐标系下的表示.利用采集到的磁力仪的Y轴输出对式(9)进行谐波分析，有：

(11)

kyBz(Δθx0cosβ+Δθx1cosβ+Δθy0sinβ+

Δθy1sinβ)-kyBxsinβ+kyBycosβ,

(12)

kyBz(Δθy0cosβ+Δθy1cosβ-Δθx0sinβ-

Δθx1sinβ)-kyBxcosβ-kyBysinβ.

(13)

由于Δθx0、Δθy0、Δθx1和Δθy1为小量，式(12)、(13)变为

(14)

解式(14)，求出β的初值β′为

(15)

精确值满足β=β′+Δβ，其中Δβ为小量.将其代回式(9)，并将等式右边的Δθx0、Δθy0、Δθx1和Δθy1项补偿至Ry1，进行整理，并重新进行谐波分析，设补偿后的磁力仪输出为Fy，有:

-kyBx(sinβ′+Δβcosβ′)+kyBy(cosβ′-

Δβsinβ′),

(16)

kyBx(-cosβ′+Δβsinβ′)-kyBy(sinβ′+

Δβcosβ′).

(17)

联立式(16)、(17)，求出Δβ的值为

(18)

至此完成了初始旋转角β的辨识.在此基础上，将等式右边的Δθx0、Δθy0、Δθx1和Δθy1项补偿至Rx1、Rz1，分别对磁力仪的X轴和Z轴输出进行谐波分析.设补偿后的输出分别为Fx和Fz.

(19)

kxBx(cosβ+Δθyxsinβ)+kxBy(sinβ-

Δθyxcosβ),

(20)

kxBx(-sinβ+Δθyxcosβ)+kxBy(cosβ+

Δθyxsinβ),

(21)

(22)

(23)

结合式(20)、(21)求出Δθyx，联立式(22)、(23)求出Δαx2和Δαy2，分别为：

(24)

(25)

考虑式(11)、(19)，可以看出Δθzy和Δθzx项与磁力仪的零偏误差发生混叠，辨识精度会受到严重影响，因此不能直接辨识.至此在第1次安装方式下，可以得到此时的初始旋转角β，磁力仪水平方向两轴之间的非正交误差Δθyx，和安装基面垂直度Δαx2和Δαy2.Δθzx，Δθzy由于不能精确标定，需要改变安装方式进行测量.

重新调整固联体的安装方式，当磁力仪X轴和Z轴位于水平，Y轴竖直时，称为第2种安装方式；当磁力仪Y轴和Z轴位于水平，X轴竖直时，称为第3种安装方式.按照第1种安装方式的步骤进行处理，此时不同安装方式下的初始旋转角将发生变化，进行谐波分析，可以求出各个辨识参数的值.对于第2种安装方式，可以求出水平面两轴垂直度Δθzx的值，对于第3种安装方式，可以求出水平面两轴垂直度Δθzy的值.

在第1次安装方式的初始位置下，磁力仪敏感轴坐标系在地理坐标系下的表示为

(26)

其中:

A=sinβ(Δθx0+Δθx1(0))-cosβ(Δθy0+Δθy1(0))+

Δαx2sinβ-Δαy2cosβ,

B=cosβ(Δθx0+Δθx1(0))+sinβ(Δθy0+Δθy1(0))+

Δαx2cosβ+Δαy2sinβ.

4 安装矩阵标定方法

磁力仪测量磁场的同时，利用星敏感器进行观星.在第1次安装方式的初始位置，假定3个星敏感器均能够观星.设某星敏感器视场内能观测到的某颗恒星的赤经为α，赤纬为δ.则在惯性坐标系下的单位恒星矢量为[13]

(27)

星敏感器观测该恒星，得到该恒星发出的光映射在该星敏感器测量坐标系下的坐标:

(28)

星敏感器在同一时刻能够观测多个恒星，得到了多组矢量对Ai和Si.根据坐标系间变换关系，有

(29)

(30)

由于3颗星敏安装工艺的限制，三轴之间存在着垂直度误差，将星敏光轴进行正交化，建立星敏感器正交坐标系OXtYtZt，其与星敏感器光轴坐标系的关系为

(31)

由于星敏感器与磁力仪通过夹具体固联，二者的夹角不会发生变化，磁力仪三轴垂直度以及星敏感器三轴垂直度也不会发生变化，所以由三轴磁力仪至星敏感器的正交安装矩阵表示为

(32)

在实际的卫星地磁测量中，通过磁力仪测量磁场，并瞬时测得星敏相对惯性空间的姿态，以及此时的经纬度和格林尼治恒星时，结合式(32)的正交安装矩阵以及辨识出的磁力仪三轴垂直度误差Δθyx，Δθzx，Δθzy， 星敏感器3个光轴的垂直度，即可解算出该位置下的地磁场矢量信息，为地磁场信息库的建立打下基础.

5 仿真验证

为了验证本文基于误差模型提出的标定方法的正确性，利用Monte Carlo方法验证.假定测量位置地磁场强度为三分量Bx为-3 275.9 nT，By为27 102.2 nT，Bz为-48 085.0 nT.各个参数初值见表1.

表1 误差参数初值

以第1次安装方式为例，选择矢量磁力仪测量噪声为σ=1.00 nT的零均值高斯白噪声，单轴转台每次安装方式等间隔旋转24次，利用Monte Carlo方法设计仿真算法辨识误差参数，共辨识100次，将计算结果绘制成曲线，如图2所示.

图2 误差参数辨识曲线

针对得到的辨识结果，计算出各个参数的残差εi(β)，εi(Δθyx)，εi(Δαx2)和εi(Δαy2).残差的计算公式和残差的标准差公式分别为：

(33)

利用式(33)求出辨识参数的不确定度为:σ(β)=2.29″，σ(Δθyx)=3.17″，σ(Δαx2)=2.20″，σ(Δαy2)=2.20″.

由于3次安装方式下的条件相同，可以认为Δθzx，Δθzy的辨识精度与Δθyx相同.利用不确定度合成公式，结合求出的各个辨识参数的不确定度即可求出磁力仪敏感轴的辨识值以及辨识精度.忽略小量，有：

(34)

(35)

(36)

假定经过长时间观星后，星敏感器光轴指向精度Δα为1″，设X星敏感器的光轴矢量为

(37)

按照不确定度的分配规律，则此时的不确定度为:

(38)

当星敏观测时刻惯性坐标系到地理坐标系的转换矩阵为

(39)

于是X星敏光轴矢量在地理系下表示和不确定度为:

(40)

(41)

σ(F·G)2=(g1σ(f1))2+(g2σ(f2))2+(g3σ(f3))2+

(f1σ(g1))2+(f2σ(g2))2+(f3σ(g3))2.

(42)

σ(L11)2=(xx|cosβ|σ(β))2+(xy|sinβ|σ(β))2+

(xzσ(Δαx2))2+(cosβσ(xx))2+

(sinβσ(xy))2+(Aσ(xz))2.

(43)

假设观测地的经纬度信息为北纬40.9°，东经115.5°，取对应的格林尼治恒星时为218.63°.利用式(1)得到惯性坐标系到地理坐标系的转换矩阵:

(44)

3个星敏感器的光轴矢量分别为:

(45)

(46)

(47)

按照上述方法计算，最终求出安装矩阵及各个元素不确定度为:

(48)

(49)

6 结 论

1)在建立误差模型的基础上，利用谐波分析法可以辨识得到包括磁力仪三轴垂直度、初始安装角、安装误差等误差参数.

2)根据辨识结果求出磁力仪测量坐标系在地理坐标系下的表示，结合星敏感器观星确定星敏感器在惯性空间姿态并转化到地理坐标系，最终推导了磁力仪测量坐标系到星敏感器正交坐标系的安装矩阵计算公式.

3)通过仿真验证，在1 nT的磁力仪测量噪声以及1″的星敏测量精度下，求出安装矩阵内各个元素的精度在1.14×10-5rad以内.