劉逸康，任順清，張高雄，張 祎

(1.哈爾濱工業(yè)大學 空間控制與慣性技術(shù)研究中心，哈爾濱 150080；2.上海衛(wèi)星工程研究所，上海 200240)

地磁場提供了豐富的地磁特征[1]，利用衛(wèi)星地磁測量，能夠快速地測量全球的地磁場矢量信息.在此基礎(chǔ)上完成地磁場模型庫的建立與更新，并廣泛應(yīng)用于地磁導航、地震預報等技術(shù)領(lǐng)域[2-3].地磁測量的精度受到測量用磁力儀和星敏感器等的儀器安裝誤差以及坐標系間傳遞精度的影響[4-7]，因此提高衛(wèi)星地磁測量系統(tǒng)內(nèi)傳感器的精度和磁力儀與星敏感器之間安裝矩陣的測量精度具有重要的實用價值.

目前對磁力儀及星敏感器之間安裝矩陣的標定方法研究相對較少，張藝騰等[8]對航空地磁測量中的磁力儀與姿態(tài)儀的聯(lián)合標定進行了研究，假設(shè)磁力儀三軸與姿態(tài)儀三軸之間為歐拉轉(zhuǎn)換的關(guān)系，通過多次旋轉(zhuǎn)，由最小二乘通過掃描迭代求出最優(yōu)的旋轉(zhuǎn)歐拉角，最終得到傳感器間的安裝矩陣.孫闖等[9]對安裝矩陣提出一種戶外標定的方法，該方法借助一個標準六面體對伸展桿系統(tǒng)進行多個方向的旋轉(zhuǎn)，測量不同姿態(tài)下星敏與磁力儀的輸出，并根據(jù)測出的星敏感器與磁力儀敏感軸坐標系相對地理坐標系的姿態(tài)辨識出星敏與三軸磁力儀間的安裝矩陣，但該方法達不到精度要求.王凱強等[10]提出了在實驗室環(huán)境下利用三軸無磁轉(zhuǎn)臺提供姿態(tài)信息以及恒星模擬器模擬室外觀星環(huán)境，通過最小二乘辨識出安裝矩陣的方法，但用三軸無磁轉(zhuǎn)臺在戶外標定既不方便也不經(jīng)濟.

本文提出了一種利用單軸立式無磁轉(zhuǎn)臺對磁力儀及星敏感器之間安裝矩陣在戶外進行地面標定的方法，通過無磁轉(zhuǎn)臺提供精確的角位置使磁力儀敏感到不同的地磁場分量的激勵，不僅辨識出磁力儀的非正交誤差，并在此之上完成安裝矩陣的標定.

1 地面標定系統(tǒng)介紹

針對三軸磁力儀以及3個星敏感器，如果各個敏感軸的方向矢量能夠統(tǒng)一在地理坐標系下表示，即可確定他們間的安裝矩陣.當星敏感器在地面觀星時，可得到星敏光軸在慣性空間的姿態(tài)，再利用當?shù)亟?jīng)緯度信息以及測量時刻的格林尼治恒星時推導出星敏光軸在地理坐標系下的表示.在單個位置下三軸磁力儀測量地磁場無法直接解算出其敏感軸的姿態(tài)，本文借助單軸無磁轉(zhuǎn)臺以及安裝夾具，在多個位置進行測量，推導出初始位置下的磁力儀敏感軸在地理坐標系下的表示.最終通過坐標系轉(zhuǎn)換等數(shù)學方法，解算出精確的安裝矩陣.根據(jù)以上原理，整個地面標定系統(tǒng)如圖1所示.

圖1 地面標定系統(tǒng)示意

磁力儀和星敏感器由剛性桿連接成固聯(lián)體，安裝在單軸立式無磁轉(zhuǎn)臺上，構(gòu)成了無磁實驗平臺.地面監(jiān)測系統(tǒng)包括旋轉(zhuǎn)軸監(jiān)測系統(tǒng)以及地磁場監(jiān)測磁強計.無磁轉(zhuǎn)臺由于必須采用無磁材料等的原因存在著較大的傾角回轉(zhuǎn)誤差，使用旋轉(zhuǎn)軸監(jiān)測裝置對轉(zhuǎn)臺的回轉(zhuǎn)誤差進行實時監(jiān)測，并對輸出進行補償.標定時需要提供精確的當?shù)卮艌鍪噶啃畔?，故采用監(jiān)測用磁強計測量當?shù)卮艌龅淖兓?

2 坐標系的建立與誤差傳遞

根據(jù)地面標定系統(tǒng)組成，依次建立以下坐標系.地心慣性坐標系OiXiYiZi，Oi為地球的地心，Xi軸和Yi軸在地球的赤道平面內(nèi)，Xi軸指向春分點，Zi軸指向地球北極，Yi軸的方向由右手定則來決定.

地理坐標系OXYZ，即東北天坐標系g.如果已知格林尼治恒星時t，測量地點的經(jīng)度λ以及緯度φ，則可求出慣性系i與地理系g間的轉(zhuǎn)換矩陣為

(1)

星敏感器光軸坐標系OpXpYpZp.由3個星敏光軸組成的非正交坐標系.OpXp沿著X星敏光軸，OpYp，OpZp分別沿著其余星敏光軸.

轉(zhuǎn)臺軸套坐標系O1X1Y1Z1固聯(lián)在軸套上.設(shè)地理坐標系為單軸無磁轉(zhuǎn)臺的參考坐標系，則軸套坐標系相對地理坐標系的變換矩陣為[11]

(2)

式中Δθx0、Δθy0稱為主軸軸線的鉛垂度，可在試驗前通過水平儀精確測量.

主軸坐標系O2X2Y2Z2，主軸坐標系相對于軸套坐標系的姿態(tài)矩陣為[11]

rot(z1,γ),

(3)

式中Δθx1(γ)、Δθy1(γ)為轉(zhuǎn)臺的傾角回轉(zhuǎn)誤差，在試驗時由旋轉(zhuǎn)軸監(jiān)視系統(tǒng)進行實時監(jiān)控.

安裝基面坐標系O3X3Y3Z3，為考慮轉(zhuǎn)臺安裝基面相對于軸系幾何軸線的垂直度Δαx2，Δαy2形成，安裝基面坐標系相對主軸坐標系的姿態(tài)矩陣為[11]

(4)

磁力儀測量坐標系O4X4Y4Z4.為正交坐標系，磁力儀位于初始位置時，O4X4軸與安裝基面O3X3軸相差了一個繞O3Z3軸的初始安裝角β，磁力儀測量坐標系相對安裝基面坐標系的姿態(tài)矩陣為

(5)

由于磁力儀3個敏感軸存在著三軸非正交誤差，初始安裝磁力儀時，令磁力儀X軸Y軸處于水平狀態(tài)，此時磁力儀敏感軸O4Z4與磁力儀敏感軸O5Z5平行，現(xiàn)稱為第1種安裝方式.設(shè)磁力儀敏感軸矢量在磁力儀測量坐標系O4X4Y4Z4的表示為

(6)

式中Δθyx、Δθzx、Δθzy分別為3個敏感軸間的垂直度.若3個參數(shù)為正，則表示磁力儀敏感軸之間夾角大于90°.由于加工水平與安裝工藝的限制，以及磁力儀的內(nèi)部存在剩磁，磁力儀輸出模型中還存在著標度因子kx,ky,kz和零偏bx,by,bz，可在試驗之前通過標定方法求出.設(shè)初始安裝位置如圖1所示.結(jié)合式(2)～(6)，此時磁力儀的輸出表示為[12]

(7)

3 磁力儀誤差標定方法

稱圖1中固聯(lián)體姿態(tài)為第1種安裝方式下的初始位置，當單軸無磁轉(zhuǎn)臺旋轉(zhuǎn)角度為γ時，采集三軸磁力儀的輸出.展開式(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)中存在的辨識參數(shù)有β，Δθyx，Δθzx，Δθzy，Δαx2和Δαy2.考慮使用諧波分析辨識，此時一周內(nèi)采樣點個數(shù)應(yīng)大于辨識量的兩倍.令γ=2πi/n(i=0,1,2,…,n-1)，n=24，即單軸轉(zhuǎn)臺等間隔旋轉(zhuǎn)24次，每次轉(zhuǎn)臺停穩(wěn)并采集完磁力儀輸出后進行下次旋轉(zhuǎn).

首先應(yīng)辨識初始旋轉(zhuǎn)角β，因只有在β已知的時才能獲取磁力儀敏感軸在地理坐標系下的表示.利用采集到的磁力儀的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)變?yōu)?/p>

(14)

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

(15)

精確值滿足β=β′+Δβ，其中Δβ為小量.將其代回式(9)，并將等式右邊的Δθx0、Δθy0、Δθx1和Δθy1項補償至Ry1，進行整理，并重新進行諧波分析，設(shè)補償后的磁力儀輸出為Fy，有:

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

Δβsinβ′),

(16)

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

Δβcosβ′).

(17)

聯(lián)立式(16)、(17)，求出Δβ的值為

(18)

至此完成了初始旋轉(zhuǎn)角β的辨識.在此基礎(chǔ)上，將等式右邊的Δθx0、Δθy0、Δθx1和Δθy1項補償至Rx1、Rz1，分別對磁力儀的X軸和Z軸輸出進行諧波分析.設(shè)補償后的輸出分別為Fx和Fz.

(19)

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

Δθyxcosβ),

(20)

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

Δθyxsinβ),

(21)

(22)

(23)

結(jié)合式(20)、(21)求出Δθyx，聯(lián)立式(22)、(23)求出Δαx2和Δαy2，分別為：

(24)

(25)

考慮式(11)、(19)，可以看出Δθzy和Δθzx項與磁力儀的零偏誤差發(fā)生混疊，辨識精度會受到嚴重影響，因此不能直接辨識.至此在第1次安裝方式下，可以得到此時的初始旋轉(zhuǎn)角β，磁力儀水平方向兩軸之間的非正交誤差Δθyx，和安裝基面垂直度Δαx2和Δαy2.Δθzx，Δθzy由于不能精確標定，需要改變安裝方式進行測量.

重新調(diào)整固聯(lián)體的安裝方式，當磁力儀X軸和Z軸位于水平，Y軸豎直時，稱為第2種安裝方式；當磁力儀Y軸和Z軸位于水平，X軸豎直時，稱為第3種安裝方式.按照第1種安裝方式的步驟進行處理，此時不同安裝方式下的初始旋轉(zhuǎn)角將發(fā)生變化，進行諧波分析，可以求出各個辨識參數(shù)的值.對于第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個星敏感器均能夠觀星.設(shè)某星敏感器視場內(nèi)能觀測到的某顆恒星的赤經(jīng)為α，赤緯為δ.則在慣性坐標系下的單位恒星矢量為[13]

(27)

星敏感器觀測該恒星，得到該恒星發(fā)出的光映射在該星敏感器測量坐標系下的坐標:

(28)

星敏感器在同一時刻能夠觀測多個恒星，得到了多組矢量對Ai和Si.根據(jù)坐標系間變換關(guān)系，有

(29)

(30)

由于3顆星敏安裝工藝的限制，三軸之間存在著垂直度誤差，將星敏光軸進行正交化，建立星敏感器正交坐標系OXtYtZt，其與星敏感器光軸坐標系的關(guān)系為

(31)

由于星敏感器與磁力儀通過夾具體固聯(lián)，二者的夾角不會發(fā)生變化，磁力儀三軸垂直度以及星敏感器三軸垂直度也不會發(fā)生變化，所以由三軸磁力儀至星敏感器的正交安裝矩陣表示為

(32)

在實際的衛(wèi)星地磁測量中，通過磁力儀測量磁場，并瞬時測得星敏相對慣性空間的姿態(tài)，以及此時的經(jīng)緯度和格林尼治恒星時，結(jié)合式(32)的正交安裝矩陣以及辨識出的磁力儀三軸垂直度誤差Δθyx，Δθzx，Δθzy， 星敏感器3個光軸的垂直度，即可解算出該位置下的地磁場矢量信息，為地磁場信息庫的建立打下基礎(chǔ).

5 仿真驗證

為了驗證本文基于誤差模型提出的標定方法的正確性，利用Monte Carlo方法驗證.假定測量位置地磁場強度為三分量Bx為-3 275.9 nT，By為27 102.2 nT，Bz為-48 085.0 nT.各個參數(shù)初值見表1.

表1 誤差參數(shù)初值

以第1次安裝方式為例，選擇矢量磁力儀測量噪聲為σ=1.00 nT的零均值高斯白噪聲，單軸轉(zhuǎn)臺每次安裝方式等間隔旋轉(zhuǎn)24次，利用Monte Carlo方法設(shè)計仿真算法辨識誤差參數(shù)，共辨識100次，將計算結(jié)果繪制成曲線，如圖2所示.

圖2 誤差參數(shù)辨識曲線

針對得到的辨識結(jié)果，計算出各個參數(shù)的殘差εi(β)，εi(Δθyx)，εi(Δαx2)和εi(Δαy2).殘差的計算公式和殘差的標準差公式分別為：

(33)

利用式(33)求出辨識參數(shù)的不確定度為:σ(β)=2.29″，σ(Δθyx)=3.17″，σ(Δαx2)=2.20″，σ(Δαy2)=2.20″.

由于3次安裝方式下的條件相同，可以認為Δθzx，Δθzy的辨識精度與Δθyx相同.利用不確定度合成公式，結(jié)合求出的各個辨識參數(shù)的不確定度即可求出磁力儀敏感軸的辨識值以及辨識精度.忽略小量，有：

(34)

(35)

(36)

假定經(jīng)過長時間觀星后，星敏感器光軸指向精度Δα為1″，設(shè)X星敏感器的光軸矢量為

(37)

按照不確定度的分配規(guī)律，則此時的不確定度為:

(38)

當星敏觀測時刻慣性坐標系到地理坐標系的轉(zhuǎn)換矩陣為

(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)

假設(shè)觀測地的經(jīng)緯度信息為北緯40.9°，東經(jīng)115.5°，取對應(yīng)的格林尼治恒星時為218.63°.利用式(1)得到慣性坐標系到地理坐標系的轉(zhuǎn)換矩陣:

(44)

3個星敏感器的光軸矢量分別為:

(45)

(46)

(47)

按照上述方法計算，最終求出安裝矩陣及各個元素不確定度為:

(48)

(49)

6 結(jié) 論

1)在建立誤差模型的基礎(chǔ)上，利用諧波分析法可以辨識得到包括磁力儀三軸垂直度、初始安裝角、安裝誤差等誤差參數(shù).

2)根據(jù)辨識結(jié)果求出磁力儀測量坐標系在地理坐標系下的表示，結(jié)合星敏感器觀星確定星敏感器在慣性空間姿態(tài)并轉(zhuǎn)化到地理坐標系，最終推導了磁力儀測量坐標系到星敏感器正交坐標系的安裝矩陣計算公式.

3)通過仿真驗證，在1 nT的磁力儀測量噪聲以及1″的星敏測量精度下，求出安裝矩陣內(nèi)各個元素的精度在1.14×10-5rad以內(nèi).