星敏與磁力儀間安裝矩陣的一種地面標定方法

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

(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以內.