柴 華 王虎彪 武 凜 王 勇

1 中國科學(xué)院測量與地球物理研究所大地測量與地球動力學(xué)國家重點實驗室，武漢市徐東大街340號，430077

利用慣性測量單元（inertial measurement unit，IMU）進行重力測量時，當(dāng)載體靜止時運動加速度為0。此時若不考慮零偏、噪聲等誤差，加速度計的直接觀測量就是重力。這一無載體運動加速度存在的狀態(tài)，可以排除由GNSS系統(tǒng)給出的載體運動加速度觀測誤差的影響，有利于分析IMU 內(nèi)置的高精度石英加速度計對地球重力變化的響應(yīng)及其誤差特性，為考察IMU 的比力觀測質(zhì)量能否達(dá)到重力測量的要求提供有利條件。

本文利用內(nèi)置有高精度石英撓性加速度計的IMU 開展慣性重力測量靜態(tài)試驗研究。通過試驗，對慣性重力觀測數(shù)據(jù)中的誤差進行初步分析與建模，探討在當(dāng)前的硬件條件下實現(xiàn)慣性重力測量的可行性及存在的問題，為后續(xù)的動態(tài)研究提供參考。

1 慣性重力測量試驗原理與誤差處理

慣性重力測量的基本原理可以表示為［1］：

利用式（2）可將載體坐標(biāo)系下的比力觀測fb通過矩陣轉(zhuǎn)換至當(dāng)?shù)厮阶鴺?biāo)系。更一般地，在不考慮噪聲與觀測誤差的前提下，當(dāng)載體靜止時加速度計的讀數(shù)皆由地球重力場引起，該點的重力值g與載體坐標(biāo)系下三軸加速度計的輸出關(guān)系可表示為［2］：

高精度捷聯(lián)式IMU 中通常內(nèi)置有石英撓性加速度計，其誤差構(gòu)成主要包括零偏項、比例因子項和白噪聲等誤差項等。在載體坐標(biāo)系下可以用數(shù)學(xué)模型表達(dá)為［1］：

式中，ba為加速度計的零偏；ka為加速度計比例因子誤差；fb為載體坐標(biāo)系下的比力觀測；wa是高頻白噪聲。零偏ba和加速度計比例因子kafb主要體現(xiàn)為常值零偏、開機零偏和漂移等誤差。常值零偏的主要部分通常已在實驗室中得到標(biāo)定與校正，因此重力測量結(jié)果主要受未補償零偏（即標(biāo)定后隨時間的漂移誤差）、開機零偏、開機后隨機漂移和溫度漂移的影響。白噪聲wa主要表現(xiàn)為加速度計觀測輸出的高頻抖動。

加速度計的誤差構(gòu)成較復(fù)雜，需根據(jù)誤差類型和特性采用不同的手段進行處理。在后續(xù)試驗中，對于高頻噪聲將利用其頻率特性對數(shù)據(jù)進行低通濾波以抑制其影響；對于未得到補償?shù)牧闫?、開機零偏，將通過在重力已知點上的對比觀測試驗來獲?。粚τ陂_機后由于儀器自身原因和觀測環(huán)境造成的漂移，將利用回歸分析方法嘗試進行擬合后補償。

2 試驗與分析

陸地慣性重力測量試驗分為兩次，第一次試驗的目的是通過長時段觀測來研究高頻觀測噪聲的去除方法與開機后漂移的擬合補償方法，并評估該IMU 內(nèi)置加速度計用于重力測量的精度與穩(wěn)定性；第二次試驗在某山區(qū)進行，即在外業(yè)環(huán)境中進行陸地慣性重力靜態(tài)測量的初步試驗。

測量中使用的IMU 置于測量車內(nèi)，表1為試驗中所使用的高精度IMU 內(nèi)置石英加速度計的技術(shù)指標(biāo)。

2.1 車載靜態(tài)測試

將測量車靜止并對IMU 的數(shù)據(jù)進行長時間靜態(tài)采集，IMU 的采樣率為200 Hz。去掉觀測起始階段受人為干擾和溫度影響較大的不穩(wěn)定觀測時段后，采集數(shù)據(jù)的總時長為7 200s。

表1 試驗中采用的高精度撓性加速度計的技術(shù)指標(biāo)Tab.1 The technical index of the accelerometers used in experiments

針對慣性元件輸出數(shù)據(jù)的特性，可以采用頻譜分析法和回歸分析法進行處理，前者適用于存在確定周期的數(shù)據(jù)而后者適用于非周期性的數(shù)據(jù)。圖1為載體坐標(biāo)系三軸的加速度讀數(shù)時間序列（截取時長為5s，共1 000 個歷元）和利用式（3）計算得到的重力觀測值的時間序列（截取時長為5s）。從圖中可以觀察到，由于加速度計各類誤差，特別是高頻觀測噪聲的存在，使單歷元重力觀測結(jié)果變化劇烈。因此在試驗一中，首先將根據(jù)高頻噪聲的特性，研究減小其對靜態(tài)重力觀測影響的方法，之后將對觀測結(jié)果中可能存在的趨勢項（長周期項）進行回歸分析，以支持后續(xù)重力觀測的補償，提高測量精度。

圖1 載體坐標(biāo)系X、Y、Z 軸加速度及單歷元計算出的重力加速度Fig.1 The body frame acceleration in X，Yand Zaxis and the calculated gravity acceleration in single epoch

2.1.1 重力觀測數(shù)據(jù)的低通濾波與頻譜分析

大量噪聲的存在使重力觀測值的時間序列無法直接通過式（3）獲得穩(wěn)定的觀測結(jié)果。重力信號通常位于低頻段，根據(jù)重力信號與觀測噪聲的特點，可以利用低通濾波器抑制高頻噪聲，恢復(fù)重力信號［3］。在重力測量領(lǐng)域，常用的低通濾波器包括巴特沃斯IIR 低通濾波器和有限沖擊響應(yīng)FIR 低通濾波器［4－5］，這些濾波技術(shù)通常被應(yīng)用于基于海空重力儀的動態(tài)重力測量數(shù)據(jù)處理中。本文采用滑動平均濾波的方法對高頻噪聲進行處理，式（5）為滑動平均濾波的基本原理［6］：

式中，g［i］為第i個歷元滑動平均后的重力觀測值；M為觀測歷元個數(shù)，它決定平滑的時間。

為選取最合適的平滑時間，采用不同平滑尺度下相對于整段觀測的標(biāo)準(zhǔn)差來評定平滑質(zhì)量，即將觀測數(shù)據(jù)分別按1s、2s、5s、10s、30s、60s、300s為長度進行濾波，結(jié)果見表2。

表2 平滑時間與精度Tab.2 Smoothing length and precision

從表2看出，隨著平滑時間的增加，觀測重力值標(biāo)準(zhǔn)差顯著減小，濾波結(jié)果越平滑，經(jīng)過30s平均后的重力觀測值已趨于穩(wěn)定。權(quán)衡觀測效率與觀測精度，也考慮后續(xù)試驗的需要，選擇30s作為濾波平滑時間長度。

分別將7 200s的原始重力觀測值及經(jīng)過30 s濾波平滑的重力觀測值進行快速傅立葉變換，得到它們的功率譜密度（圖2）。可以看出，噪聲在整個頻段均有分布，且在77.5 Hz處的功率譜密度最大，即該頻率的高頻噪聲對測量結(jié)果影響最為顯著，但經(jīng)30s平滑低通濾波后高頻噪聲的影響已得到有效抑制。

圖2 原始重力觀測值的功率譜密度（a）與經(jīng)30s平滑后重力觀測值的功率譜密度（b）Fig.2 The power spectrum density comparison between original calculated gravity acceleration（a）and the gravity acceleration after 30ssmoothing（b）

2.1.2 重力觀測數(shù)據(jù)的回歸分析

慣性元器件的漂移將使靜態(tài)重力的觀測結(jié)果隨時間變化，本節(jié)主要考慮儀器溫度與觀測環(huán)境的不穩(wěn)定性引起的漂移及其對重力觀測的影響。采用與上一組試驗相同的數(shù)據(jù)，將經(jīng)30s濾波平滑后靜態(tài)重力觀測的時間序列繪制于圖3，可以看出重力觀測值存在明顯的上升趨勢。

針對該數(shù)據(jù)進行回歸分析，可構(gòu)造多項式函數(shù)：

若采用線性擬合方法，多項式函數(shù)可簡化為：

式中，斜率p1代表加速度計的漂移速率，與加速度計的材質(zhì)、制造工藝和觀測環(huán)境溫度有關(guān)；截距l(xiāng)為擬合起始時刻的重力觀測值，該值與該點真實重力值之間的誤差主要由未得到補償?shù)某Ｖ盗闫兔看伍_機的非固定零偏引起。

對經(jīng)過30s濾波后的數(shù)據(jù)采用最小二乘法求出其最優(yōu)線性擬合多項式：

由式（8）可知，開機后石英撓性加速度計的隨機游走、環(huán)境變化引起的漂移使7 200s（1 440 000個歷元）內(nèi)重力觀測值增大約4.0mGal。

2.2 車載零速修正點上的靜態(tài)測試

零速修正（ZUPT）能有效減小慣導(dǎo)系統(tǒng)中誤差隨時間的積累，該過程需要載體維持靜止?fàn)顟B(tài)數(shù)10s。試驗一已經(jīng)證明，通過30s的濾波平滑可以使測試中采用的IMU 得到標(biāo)準(zhǔn)差小于1.5 mGal的靜態(tài)重力觀測值，因此在慣性導(dǎo)航中可以利用加速度計的觀測數(shù)據(jù)獲得零速修正點上的重力值。

試驗二在某山區(qū)一條長約48km 的測線上進行。IMU 被安置于動態(tài)測量車上，用于采集ZUPT 點上的重力值。試驗中部分零速修正點的重力值已提前通過相對重力測量獲得。測線的大地高從123m 變化至390m，根據(jù)垂向重力梯度改正公式：

式中，h為高程變化。試驗中大地高變化為267 m，對應(yīng)的重力變化理論值約為82.4mGal，較大的重力變化有利于測量成果的檢驗。

考慮到石英撓性加速度計的不穩(wěn)定性會使相隔較短時間的加速度計的零偏值發(fā)生變化，且慣性測量系統(tǒng)每次開機引起的零偏量亦不相同，故在試驗二中由IMU 重力測量的初始偏移量bias0需在重力已知點上進行觀測來確定。由于試驗一與試驗二采用同型號的IMU 且都是在冬季進行，氣溫較低，外部環(huán)境溫度近乎一致。在溫度與外界環(huán)境對測量結(jié)果的影響大致相同的前提下，斜率p1將沿用試驗一中獲得的參數(shù)。因此，相應(yīng)補償量的計算公式為：

式中，t0為確定bias0時刻的歷元數(shù)，ti為任一零速修正點上的歷元數(shù)。假設(shè)載體運動對漂移特性無顯著影響，最后得到零速修正點上補償后的重力觀測值為：

將測線上所有ZUPT 點的重力觀測值經(jīng)式（10）、（11）改正后作圖4。圖中的空心點及其連線表示零速修正點上的IMU 觀測重力值，離散的實心點表示對應(yīng)零速修正點上的已知重力值。為排除儀器開機時升溫對測量的影響，試驗將選擇第2個重力值已知的ZUPT 點作為初始重力已知點，經(jīng)計算得到該點起始重力觀測誤差，即初始偏移量bias0為12.3mGal。從圖4看出，在測線上大地高升高的同時，IMU 給出的重力觀測值減小，基于慣性重力測量結(jié)果正確地反映了重力的變化趨勢。

圖4 ZUPT 點上的重力值對比Fig.4 The gravity value comparison on ZUPT points

表3為觀測值已知的零速修正點上重力測量精度統(tǒng)計。點1上補償后誤差高達(dá)－20.4mGal，經(jīng)分析主要是受開機時IMU 內(nèi)部升溫影響所致。由于試驗在環(huán)境溫度較低的冬天進行，開機后IMU 內(nèi)部溫度變化較大，使石英撓性加速度計的觀測受到嚴(yán)重影響，故在試驗中沒有將該點作為初始的重力已知點。

表3 重力測量精度統(tǒng)計Tab.3 The precision statistic of gravity measurement

除起始點1外，補償后的IMU 觀測重力值誤差均較補償前顯著減小，大多數(shù)點上的重力觀測精度由補償前的12～20 mGal提升至優(yōu)于2 mGal，說明本文試驗采用的補償方法有效，同時也證明利用IMU 開展陸地mGal級重力測量的可能性。從表3中可看出，經(jīng)過誤差補償后點3、4精度最優(yōu)，誤差小于1mGal；點5、6精度次之，誤差分別為1.8 和－1.8 mGal；點7 誤差約6.1 mGal，懷疑與車輛振動等相對復(fù)雜的外部觀測環(huán)境有關(guān)。

3 結(jié) 語

本文利用高精度的慣性測量設(shè)備開展測定重力場的初步靜態(tài)測試，證明了基于試驗型號的高精度捷聯(lián)式IMU 內(nèi)置加速度計能夠?qū)崿F(xiàn)mGal級的重力觀測精度，同時也為后續(xù)試驗提供了一些經(jīng)驗：

1）試驗中儀器開機后IMU 內(nèi)部的升溫對重力測量有較大影響，已超出當(dāng)前高精度加速度計溫度模型補償?shù)哪芰Ψ秶９试诤罄m(xù)的試驗中，一方面有必要在測量前對儀器進行預(yù)熱，以減小開機后的溫度效應(yīng)；另一方面，可從硬件著手對系統(tǒng)溫度進行控制，例如為石英撓性加速度計提供恒溫環(huán)境，以獲得穩(wěn)定的觀測值。

2）在相同的硬件和相似的溫度環(huán)境下，開機后重力觀測值的漂移具有一定的復(fù)現(xiàn)性，因此可提前測定并用于后續(xù)觀測的補償。在動態(tài)試驗中，可以考慮采用GNSS與SINS聯(lián)合數(shù)據(jù)處理的方式，實時地對IMU 給出的比力觀測誤差進行求定與補償。

致謝：對德國達(dá)姆斯塔特工業(yè)大學(xué)（TU Darmstadt）物理與衛(wèi)星大地測量研究所（PSGD）提供的試驗機會與幫助表示感謝。

［1］Jekeli C.Inertial Navigation Systems with Geodetic Applications［M］.Berlin：Walter de Gruyter，2001

［2］Gerlach C，Dorobantu R，Rothacher M.Results of a Combined INS／GPS Experiment for Geodetic Application［J］.Navigation，2005，53（212）：31－47

［3］董緒榮，張守信，華仲春.GPS／INS組合導(dǎo)航定位及其應(yīng)用［M］.長沙：國防科技大學(xué)出版社，1998（Dong Xurong，Zhang Shouxin，Hua Zhongchun.GPS／INS Integrated System Navigation，Positioning and Its Applications［M］.Changsha：Publishing Company of National University of Defense Technology，1998）

［4］張開東.基于SINS／DGPS的航空重力測量方法研究［D］.長沙：國防科技大學(xué)，2007（Zhang Kaidong.Research on the Methods of Airborne Gravimetry Based on SINS／DGPS［D］.Changsha：National University of Defense Technology，2007）

［5］孫中苗，夏哲仁.FIR 低通差分器的設(shè)計及其在航空重力測量中的應(yīng)用［J］.地球物理學(xué)報，2000，43（6）：850－855（Sun Zhongmiao，Xia Zheren.Design of FIR Lowpass Differentiator and Its Applications in Airborne Gravimetry［J］.Chinese J Geophys，2000，43（6）：850－855）

［6］Steven W S.Digital Signal Processing：A Practical Guide for Engineers and Scientists［M］.Boston：Newnes，2003

［7］柴華，王勇，王虎彪，等.GNSS／SINS組合進行慣性重力測量的誤差分析［J］.大地測量與地球動力學(xué)，2011，31（6）：73－78（Chai Hua，Wang Yong，Wang Hubiao，et al.Error Analysis for Inertial Gravimetry by Use of GNSS／SINS Combination［J］.Journal of Geodesy and Geodynamics，2011，31（6）：73－78）