王永彤，朱志剛，郭宗本，徐 超

（北京航天控制儀器研究所，北京100039)

三浮陀螺伺服測試誤差系數(shù)估計方法

王永彤，朱志剛，郭宗本，徐 超

（北京航天控制儀器研究所，北京100039)

為了更好地解決工程實際中三浮陀螺極軸伺服測試時誤差系數(shù)估計問題，采用了Kalman濾波和遞推最小二乘的方法。試驗結果證明，Kalman濾波方法能夠有效地估計出極軸伺服測試時的陀螺誤差系數(shù)。用統(tǒng)計檢驗的方法分析了該模型以及各系數(shù)的顯著性，結果顯示對Dvz和Dzz的估計是顯著的。

三浮陀螺；伺服測試；Kalman濾波；誤差系數(shù)估計；顯著性分析

0 引言

伺服轉臺測試是高精度單自由度機電陀螺測試方法之一，該方法的優(yōu)點是：1)測試精度不受陀螺力矩器精度的限制；2)精度潛力大，由于不直接測量漂移角速度，而是測量漂移角速度隨時間的積分，在一定程度上可以用時間換精度；3)測量過程更接近陀螺在平臺系統(tǒng)中的工作狀態(tài)，對實際應用有更好的指導意義。為提高測試精度、豐富陀螺精度鑒定手段，以及給更高精度機電陀螺奠定測試技術基礎，提高工程化應用，亟待開展陀螺的伺服轉臺測試方法、數(shù)據(jù)處理及誤差系數(shù)標定的研究。

當陀螺輸入軸平行地理極軸做伺服測試時，陀螺誤差系數(shù)的辨識可以采用諧波系數(shù)法、最小二乘估計法和Kalman濾波估計方法［1］。用諧波系數(shù)法時，需要陀螺沿輸入軸相對地理系旋轉至少一周才能解出誤差系數(shù)，由于耗時長、實時性差，難以滿足工程應用。最小二乘估計法使量測估計的精度達到最佳，不使用與估計量有關的動態(tài)信息和統(tǒng)計信息，優(yōu)點是算法簡單，在對被估計量和量測誤差缺乏了解的情況下仍能適用［3］。Kalman濾波是線性最小方差估計，使用時需要知道被估計量驅動噪聲的統(tǒng)計特性、系統(tǒng)狀態(tài)方程，以及量測噪聲的統(tǒng)計特性［3］。最小二乘和Kalman濾波各有優(yōu)缺點，為提高陀螺誤差系數(shù)的估計精度，提升陀螺伺服試驗效率，同時解決估計方法在實際測試中的實用性，還需進行對比研究。本文根據(jù)陀螺漂移的實際特性并結合極軸伺服測試時的特點采用了與主軸角度相關的陀螺誤差模型，并在此基礎上采用Kalman濾波器與最小二乘法對模型參數(shù)進行辨識。隨后，對模型的回歸方程及各系數(shù)進行了顯著性分析。

1 陀螺極軸伺服試驗

三浮陀螺在1g重力場下誤差模型如式（1)所示［1］：

式中，ωd為陀螺的綜合漂移，（°)／h；Df為不敏感比力的漂移角速度，（°)／h；Dx、Dy、Dz為與重力加速度成比例的誤差系數(shù)，（°)／h／g；Dxx、Dzz為與重力加速度平方成比例的誤差系數(shù)，（°)／h／g2；Dxy、Dxz、Dyz為與重力加速度乘積成比例的誤差系數(shù)，（°)／h／g2；gx、gy、gz分別為沿x軸、y軸、z軸的重力加速度分量，g；ωm為陀螺的動態(tài)相關噪聲，（°)／h；ωε為陀螺誤差模型殘差，（°)／h。

x、y、z分別對應陀螺儀的輸入軸IA，輸出軸OA，自轉軸SA。

陀螺伺服測試時常采用極軸伺服試驗考察陀螺的精度［1］。陀螺極軸伺服試驗有兩種安裝方位：陀螺儀輸入軸平行極軸指北和指南。陀螺儀坐標系在地球坐標系（O′XYZ)和地理坐標系（OXeYeZe)中的方位如圖1所示。在這種安裝方式下，陀螺始終敏感地速。重力在陀螺輸入軸上投影為常值，在陀螺自轉軸和輸出軸上投影根據(jù)轉臺主軸轉動的角位置呈現(xiàn)周期性變化。

在圖1中，φ為當?shù)鼐暥?；ωe為地球轉速；O′XYZ為地球坐標系；OXeYeZe（OUNW）為地理坐標系，OXe指天，OYe在當?shù)厮矫鎯戎副?；oxyz為陀螺儀坐標系，ox、oy、oz分別為陀螺儀的輸入軸IA、輸出軸OA、自轉軸SA；OX2Y2Z2為伺服轉臺主軸坐標系，OX2、OZ2分別為轉臺主軸、耳軸，轉臺零位時OX2與OXe重合，OY2指向OZe反向，OZ2與OYe重合。

圖1 極軸伺服試驗陀螺儀坐標系方位圖Fig．1 Coordinates in servo test

陀螺極軸伺服試驗時，有如下運動方程：

其中，˙θ為轉臺繞轉臺主軸的旋轉角速度，可由角度差分得到，（°)／h；ωx為地速在陀螺輸入軸上的投影，（°)／h。

陀螺儀輸入軸平行極軸指北的零位方位定義為：陀螺輸入軸IA平行極軸指北，馬達軸SA指東。重力加速度g朝上為正，伺服轉臺沿主軸坐標系OX2軸逆時針轉動為正，轉角為θ，正東方向為角度起始零位。地速在陀螺輸入軸上投影為ωe，重力加速度在陀螺各軸上的投影如下：

其中，gs＝gsinφ，gc＝gcosφ。

陀螺儀輸入軸平行極軸指南的零位方位定義為：陀螺輸入軸IA平行極軸指南，馬達軸SA指東。地速在陀螺輸入軸上投影為-ωe，重力加速度在陀螺各軸上的投影如下：

將式（3)帶入式（1)后采用諧波分析的方法，得到陀螺漂移是一個關于轉臺轉角θ的周期函數(shù)。陀螺漂移可以用式（5)表示，整理后得式（6)，式中各系數(shù)分為與轉臺轉角無關的項、與轉臺轉角一次諧波有關的項、與轉臺轉角二次諧波有關的項。

其中，

ωdN為陀螺輸入軸指北時的綜合漂移。

其中，A0N、A1N、A2N、B1N、B2N如式（8)所示。

同理可得到陀螺輸入軸指南時的諧波系數(shù)方程，A0S、A1S、A2S、B1S、B2S如式（9)所示。陀螺誤差系數(shù)可由諧波系數(shù)通過代數(shù)運算解出，因此辨識誤差系數(shù)的關鍵是得到諧波系數(shù)。

2 伺服測試陀螺誤差系數(shù)估計方法

2．1 Kalman濾波估計方法

三浮陀螺動態(tài)相關噪聲ωm為1階馬氏過程，其相關函數(shù)為：

其中，RN（0)為均方值，α為反相關時間常數(shù)，τ為兩點之間的時間間隔。RN（0)和α可以通過試驗進行估算［2］。

將其白化離散后的表達式［3］為：

其中，T為離散時間；wk為白噪聲序列。

陀螺伺服測試時，陀螺漂移的系統(tǒng)模型和量測模型方程［1，3］分別為：

根據(jù)式（6)、式（7)和式（11)得到Kalman濾波方程中的各項如下：

θk為k時刻伺服轉臺主軸角度。

觀測量Zk為伺服轉臺角速率，由角度差分得出。Vk為伺服轉臺角速率的量測噪聲，角速率的隨機誤差主要由測角誤差造成，設測角誤差的峰峰值為ε，ΔT為動平均時間［1］，角速率讀數(shù)誤差最大值為ε／ΔT，此誤差近似為正態(tài)分布，可視為白噪聲。ΔT選取時需考慮角速率的平滑效果和陀螺動態(tài)誤差的相關時間，這里為2000s。量測方差估計為：

系統(tǒng)噪聲方差為：

由試驗數(shù)據(jù)計算可得，RN（0）＝1.3×10-7［（°)／h］2，α＝2000，T＝1s。

Kalman濾波方程如下［1，3］：

X0和P0初始值均賦0。

2．2 最小二乘估計方法

為了對比Kalman濾波效果，對三浮陀螺誤差系數(shù)用最小二乘方法進行估計，式（7)表達可以如式（14)所示，Zi為第i次轉臺角速度量測值，Vi為第i次隨機量測噪聲，Hi＝［1 sinθisin2θicosθicos2θi］，為提高數(shù)據(jù)處理效率，采用遞推最小二乘方法［3］估計諧波系數(shù)。

為避免估計誤差跳躍劇烈，＾X0和P0可由初始數(shù)據(jù)經(jīng)最小二乘計算得到［3］，如式（16）

3 試驗分析

3．1 試驗數(shù)據(jù)分析

某陀螺極軸伺服試驗的主軸角速度數(shù)據(jù)如圖2和圖3所示，采樣時間為1s。為獲取更多數(shù)據(jù)進行對比，測試時間較長。其中，輸入軸指北、指南的數(shù)據(jù)樣本長度分別為52.7h和39h，轉臺分別轉過579.3°（1.61圈)和523°（1.45圈)，轉臺繞主軸旋轉一周分別為32.8h和26.9h。用Kalman濾波估計方法，陀螺相關噪聲和諧波系數(shù)的收斂穩(wěn)定時間均在6.5h左右，收斂時間占主軸旋轉一周時間的21％和 26％。辨識結果如圖4和圖5所示。各諧波系數(shù)體現(xiàn)了所包含陀螺誤差系數(shù)的綜合精度，如式（8)、式（9)所示。用遞推最小二乘估計方法辨識結果如圖6和圖7所示。為避免估計誤差起伏劇烈，用前4000個數(shù)據(jù)做最小二乘估計，將其設為遞推初值。各諧波系數(shù)收斂時間約為4.5h左右。將陀螺輸入軸指北、南兩次試驗得到的諧波系數(shù)做代數(shù)運算后，可以得到陀螺各誤差系數(shù)，結果如表1所示。

圖2 陀螺輸入軸指北時轉臺相對主軸角速度Fig．2 The angular velocity of the turntable which axis is parallel to the earth's north polar axis

圖3 陀螺輸入軸指南時轉臺相對主軸角速度Fig．3 The angular velocity of the turntable which axis is parallel to the earth's south polar axis

圖4 陀螺相關噪聲與諧波系數(shù)濾波收斂過程（陀螺輸入軸指北）Fig．4 The identification results of the model（the turntable's axis is parallel to the earth's north polar axis）

圖5 陀螺相關噪聲與諧波系數(shù)濾波收斂過程（陀螺輸入軸指南）Fig．5 The identification results of the model（the turntable's axis is parallel to the earth's north polar axis）

圖6 用遞推最小二乘法估計陀螺諧波系數(shù)收斂過程（陀螺輸入軸指北）Fig．6 The identification results of the model by recursive least?squares（the turntable's axis is parallel to the earth's north polar axis）

圖7 用遞推最小二乘法估計陀螺諧波系數(shù)收斂過程（陀螺輸入軸指南）Fig．7 The identification results of the model recursive least?squares（the turntable's axis is parallel to the earth's south polar axis）

表1 陀螺誤差系數(shù)估計結果（設Dxx＝0）Table 1 The determined coefficients of Gyro（Dxx＝0）

以陀螺輸入軸平行極軸指北的數(shù)據(jù)為例，用Kalman濾波辨識的系數(shù)經(jīng)擬合后與實測漂移相減的殘差如圖8所示，標準差如表2所示。與遞推最小二乘的相比，Kalman濾波擬合后的殘差標準差小，如圖9所示。

圖8 Kalman濾波擬合后的陀螺漂移殘差Fig．8 The drift residuals with Kalman filter

表2 陀螺漂移角速度殘差標準差（陀螺輸入軸指北）［（°）／h］Table 2 The standard deviation of the gyro drift residuals（the turntable's axis is parallel to the earth's north polar axis）

圖9 最小二乘方法擬合后的陀螺漂移殘差Fig．9 The drift residuals with LSM

3．2 顯著性分析

估計出陀螺誤差系數(shù)后，還需對誤差方程及各系數(shù)進行統(tǒng)計檢驗。

假定第j次被測量yj和參數(shù)a0、a1、…、an之間呈線性關系n，j＝1，2，…，N，N＞n，xji為可測的已知變量，ai為n待估計的未知參數(shù)，ej是N個相互獨立服從同一正態(tài)分布N（0，σ2)的隨機變量。有觀測方程Y＝XA＋ε?？偲钇椒胶蜑榛貧w平方和為殘差平方和為可證明ST＝SR＋SE。

如果變量y與xi之間沒有線性關系，只需檢驗假設H0：a1＝0，…，ai＝0，這可以通過比較SR和SE來實現(xiàn)。可以證明X滿秩和假設H0成立時，SR／σ2～χ2（n)，SE／σ2～χ2（N-n-1)，SR和SE相互獨立，從而：

對給定的一組數(shù)據(jù)算得F＞Fa，則在顯著水平α下，認為回歸方程是有顯著意義的。反之，則認為該線性回歸方程沒有顯著意義。

此外，還需對每個變量進行考察，如果某個變量對y作用不明顯，在模型中它前面的系數(shù)ai可取值為0。檢驗因子xi是否顯著可以等價于是服從正態(tài)分布的隨機變量上的第i個對角元與SE相互獨立，在假設H0∶ai＝ 0條件下，計算統(tǒng)計量：

如果F＞Fa，則認為在顯著水平α下回歸系數(shù)顯著，反之則認為不顯著。

在估計陀螺誤差系數(shù)時，根據(jù)陀螺極軸伺服試驗的特點，XTX的最大值為5，因此將陀螺誤差模型轉化為式（14)。對該線性回歸模型及各系數(shù)進行顯著性檢驗。其中，y＝ωd，a0、a1、…、an分別為A0、A1、A2、B1、B2，Vi近似為正態(tài)分布。

陀螺輸入軸指北時，SEN＝0.044，SRN＝1.27× 105，N＝129727，n＝5，經(jīng)計算得回歸方程統(tǒng)計量FN＝7.5×106，在顯著水平α＝0.01下，F(xiàn)0.01（5，129721)＝3.02，因此認為回歸方程是顯著的。

再檢驗各誤差系數(shù)。統(tǒng)計量F的下角標代表陀螺各系數(shù)，N指陀螺輸入軸指北。對角元素pii計算為經(jīng)計算得FA0N＝5.8×1012，F(xiàn)A1N＝3.4×1011，F(xiàn)A2N＝2.8×104，F(xiàn)B1N＝3.2×106，F(xiàn)B2N＝9.5×104。在顯著水平α＝0.01下，F(xiàn)0.01（1，129721)＝6.63。因此，各系數(shù)是顯著的。當陀螺輸入軸指南時，也能得出同樣的結論。

4 結論

本文在三浮陀螺做極軸伺服測試時，用Kalman濾波器和遞推最小二乘方法估計陀螺誤差系數(shù)，并對最小二乘回歸方程以及各諧波系數(shù)進行了顯著性分析，試驗結果表明：

1)所采用的Kalman濾波器是一致漸進收斂的，相對遞推最小二乘方法，用Kalman濾波估計法的陀螺漂移殘差標準差小。

2)在諧波系數(shù)估計的收斂時間方面，遞推最小二乘法比Kalman濾波收斂的快，但遞推最小二乘的收斂過程對各向量的初值較敏感，數(shù)據(jù)需要預先處理后才能進入迭代運算，而Kalman濾波對系統(tǒng)初值無特殊要求。

3)統(tǒng)計檢驗表明，陀螺誤差模型的諧波系數(shù)回歸方程以及估計的各諧波系數(shù)是顯著的。由此可以判斷陀螺誤差系數(shù)Dyz和Dzz顯著，但是尚不能給出數(shù)值較小的Dxy、Dxz和Dy是否顯著的結論。

當三浮陀螺進行伺服測試時，在模型及參數(shù)準確的前提下，用Kalman濾波的方法估計誤差系數(shù)是一種更好的選擇。需要注意的是，當測試陀螺長時間精度時，為了防止新量測值對估計值修正作用下降，陳舊量測值修正作用相對上升，可以根據(jù)最優(yōu)算法結合實際情況采用漸消記憶或限定記憶的方法辨識誤差模型。極軸伺服測試得到的陀螺某些誤差系數(shù)數(shù)值較小，一些二次項也無法用統(tǒng)計檢驗進行顯著性分析，建議對這些誤差項在高g下進行辨識。

［1］胡恒章.陀螺儀漂移測試原理及其實驗技術［M］.北京：國防工業(yè)出版社，1981. HU Heng?zhang.Test theory and technology of gyro drift［M］.Beijing：National Defense Industry Press，1981.

［2］王廣雄.控制系統(tǒng)設計［M］.北京：中國宇航出版社，1992. WANG Guang?xiong.Control system design［M］.Beijing：China Aerospace Press，1992.

［3］秦永元，張洪鉞，汪叔華.卡爾曼濾波與組合導航原理［M］.西安：西北工業(yè)大學出版社，1998. QIN Yong?yuan，ZHANG Hong?yue，WANG Shu?hua. Kalman filter and integrated navigation theory［M］.Xi'an：Northwestern Polytechnical University Press，1998.

［4］魏宗康，夏剛，朱宗榮，等.陀螺儀漂移誤差參數(shù)估計的方案設計［J］.導航與控制，2003，2（1)：59?62. WEI Zong?kang，XIA Gang，ZHU Zong?rong，et al.Cali?bration scheme of gyro's drift errors［J］.Navigation and Control，2003，2（1)：59?62.

［5］張克志，田蔚風，錢峰，等.靜電陀螺靜態(tài)漂移誤差模型系數(shù)的UKF標定方法［J］.上海交通大學學報，2008，42（8)：1387?1391. ZHANG Ke?zhi，TIAN Wei?feng，QIAN Feng，et al.Un?scented Kalman filter calibration method of electrostatically suspended gyroscope drift［J］.Journal of Shanghai Jiaotong University，2008，42（8)：1387?1391.

［6］吳廣玉.系統(tǒng)辨識與自適應控制（下)［M］.哈爾濱：哈爾濱工業(yè)大學出版社，1987. WU Guang?yu.System identification and self?adaptation control（2nd)［M］.Harbin：Harbin Institute of Technology Press，1987.

［7］安學彬，萬彥輝.三浮陀螺儀隨機漂移模型研究［J］.壓電與聲光，2009，31（5)：643?645. AN Xue?bin，WAN Yan?hui.Research on random drift model of three?floated gyro［J］.Piezoelectrics＆Acous?tooptics，2009，31（5)：643?645.

［8］李璐，孫純祥.三浮陀螺隨機漂移數(shù)據(jù)建模［J］.壓電與聲光，2011，33（3)：378?381. LI Lu，SUN Chun?xiang.Modeling of random drift data of the liquid floated gyro［J］.Piezoelectrics＆Acoustooptics，2011，33（3)：378?381.

［9］何曉群，劉文卿.應用回歸分析［M］.北京：中國人民大學出版社，2015. HE Xiao?qun，LIU Wen?qing，Applied regression analysis［M］.Beijing：China Renmin University Press，2015.

［10］劉建波，魏宗康.石英加速度計誤差系數(shù)顯著性分析［J］.中國慣性技術學報，2011，19（5)：615?620. LIU Jian?bo，WEI Zong?kang.Significance analysis of QFPA's error model's coefficients［J］.Journal of Chinese Inertial Technology，2011，19（5)：615?620.

［11］嚴恭敏，李四海，秦永元.慣性儀器測試與數(shù)據(jù)分析［M］.北京：國防工業(yè)出版社，2015. YAN Gong?min，LI Si?hai，QIN Yong?yuan.Inertial de?vice test and data analysis［M］.Beijing：National Defense Industry Press，2015.

Calibration of SDFG Coefficient Based on Servo Test

WANG Yong?tong，ZHU Zhi?gang，GUO Zong?ben，XU Chao
（Beijing Institute of Aerospace Control Devices，Beijing 100039)

When single?degree?of?freedom rate?integrating gyro（SDFG)in servo test，to calibrate the coefficient of SDFG whose axis of rotation is nominally parallel to the earth's polar axis，Kalman filter and recursive least?squares method were adap?ted in this paper.The data show that the Kalman filter method is efficient.Subsequently，the regression equation and coefficients were analyzed by statistical test，and the results verified that the calibrated coefficient Dvzand Dzzis creditable.

single?degree?of?freedom rate?integrating gyro;servo test;Kalman filter;error coefficient estimation; prominence analysis

U666.16＋3

A

1674?5558（2017)05?01254

10.3969／j.issn.1674?5558.2017.01.018

王永彤，女，高級工程師，研究方向為導航、制導與控制。

2016?03?16