柴 華 王虎彪 武 凜 王 勇

1 中國科學(xué)院測(cè)量與地球物理研究所大地測(cè)量與地球動(dòng)力學(xué)國家重點(diǎn)實(shí)驗(yàn)室，武漢市徐東大街340號(hào)，430077

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

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

1 慣性重力測(cè)量試驗(yàn)原理與誤差處理

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

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

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

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

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

2 試驗(yàn)與分析

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

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

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

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

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

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

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

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

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

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

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

表2 平滑時(shí)間與精度Tab.2 Smoothing length and precision

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

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

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

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

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

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

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

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

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

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

2.2 車載零速修正點(diǎn)上的靜態(tài)測(cè)試

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

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

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

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

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

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

圖4 ZUPT 點(diǎn)上的重力值對(duì)比Fig.4 The gravity value comparison on ZUPT points

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

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

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

3 結(jié) 語

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

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

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

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

［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］.長(zhǎng)沙：國防科技大學(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的航空重力測(cè)量方法研究［D］.長(zhǎng)沙：國防科技大學(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è)計(jì)及其在航空重力測(cè)量中的應(yīng)用［J］.地球物理學(xué)報(bào)，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ìn)行慣性重力測(cè)量的誤差分析［J］.大地測(cè)量與地球動(dòng)力學(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）