錢學(xué)武，蔡體菁

（東南大學(xué) 儀器科學(xué)與工程學(xué)院，南京 210096）

旋轉(zhuǎn)加速度計重力梯度儀標定方法

錢學(xué)武，蔡體菁

（東南大學(xué) 儀器科學(xué)與工程學(xué)院，南京 210096）

提出一種旋轉(zhuǎn)加速度計重力梯度儀標定方法，在重力梯度儀外側(cè)空間確定的四個位置上依次放置一定質(zhì)量的檢測質(zhì)量體，根據(jù)檢測質(zhì)量體在正交位置上引起的重力梯度大小相等，符號相反的關(guān)系確定出重力梯度標度系數(shù)，采用旋轉(zhuǎn)重力梯度儀本體方式確定出重力梯度零位。給出了重力梯度零位和標度系數(shù)計算表達式，在重力梯度半物理仿真系統(tǒng)上進行了仿真試驗驗證，仿真結(jié)果表明，標定后重力梯度測量誤差小于1E。

重力梯度儀；旋轉(zhuǎn)加速度計；檢測質(zhì)量體；梯度標定

高精度重力梯度測量在地球科學(xué)、空間科學(xué)、地質(zhì)科學(xué)、能源勘探以及慣性導(dǎo)航等領(lǐng)域具有非常重要的應(yīng)用價值和廣泛的應(yīng)用前景。重力梯度測量相對于重力測量具有明顯的優(yōu)勢，對于地質(zhì)結(jié)構(gòu)的形狀、尺寸和埋藏深度等信息，重力梯度測量能夠提供更高的解析度，特別適合于探測局部地質(zhì)結(jié)構(gòu)及其細節(jié)[1-4]。目前，國外正在進行研究具有應(yīng)用前景的重力梯度儀主要有旋轉(zhuǎn)加速度計重力梯度儀、超導(dǎo)重力梯度儀、冷原子重力梯度儀等，其中旋轉(zhuǎn)加速度計重力梯度儀技術(shù)最成熟，也是目前唯一成功用于機載/船載動機座并進行商用的重力梯度測量系統(tǒng)，已經(jīng)進行了大量的勘探工作，并取得了良好的效果[5-8]。20世紀90年代國內(nèi)開始研究重力梯度儀[9]，目前有的科研單位已經(jīng)研制出了重力梯度儀實驗室樣機[10-12]。

旋轉(zhuǎn)加速度計重力梯度儀是一種高精密測量儀器，在出廠前或使用前需要對其進行標定。目前，關(guān)于旋轉(zhuǎn)加速度計重力梯度儀的標定方法在公開刊物上沒有報道。

重力梯度儀標定一般是在實驗室環(huán)境條件下進行，在重力梯度儀外側(cè)放置形狀規(guī)則的均質(zhì)檢測質(zhì)量體，通過改變檢測質(zhì)量體與GGI圓盤中心的距離來改變GGI中心處的引力場，變化的引力場被重力梯度敏感器敏感到，然后對重力梯度儀的輸出信號進行重力梯度解調(diào)，就可以獲得距離與GGI測量梯度之間的對應(yīng)關(guān)系，根據(jù)這種對應(yīng)關(guān)系與重力梯度理論值進行對比，可以確定出重力梯度零位和標度系數(shù)，但由于引力梯度對距離敏感，檢測質(zhì)量體距GGI較近時會產(chǎn)生較大畸變，從而導(dǎo)致較大標定誤差。為此，本文提出一種旋轉(zhuǎn)加速度計重力梯度儀部分張量的高精度標定方法。

1 重力梯度標定方法

重力梯度儀主要有全張量和部分張量兩種類型，針對全張量重力梯度儀需要對五個獨立重力梯度元素進行標定，全張量重力梯度儀標定方法與部分張量重力梯度儀標定方法類似，本文只對部分張量重力梯度儀標定方法進行分析。對于部分張量重力梯度儀，需要確定重力梯度分量(Гyy-Гxx)和 Гxy的零位和標度系數(shù)，重力梯度儀誤差模型方程為

式中：Γc為重力梯度儀實際測量值；K為重力梯度標度系數(shù)；Γt為檢測質(zhì)量體引起的理論重力梯度，即重力梯度精確值；G0為周圍環(huán)境物體引起的重力梯度，Γ0為重力梯度儀儀器本身性能引起的重力梯度值，即為重力梯度零位，重力梯度零位主要由GGI動力學(xué)擾動、線路噪聲竄擾等因素引起，不考慮時變影響，它是一個定值。從式(1)可以得到，G0和 Γ0相互耦合，不易提取出重力梯度零位，而G0主要與環(huán)境物體相對于GGI的位置有關(guān)，只要環(huán)境物體相對位置不改變，其值也不會改變。旋轉(zhuǎn)重力梯度儀本體相當(dāng)于環(huán)境物體相對于GGI位置的改變，若重力梯度儀本體旋轉(zhuǎn)整數(shù)倍周期，那么G0將以周期性呈現(xiàn)，利用G0的周期性和Γ0的定值性，可以確定出重力梯度零位Γ0。

以GGI圓盤中心為原點O，以其所在平面為xy平面建立圓盤坐標系。假設(shè)檢測質(zhì)量體為均質(zhì)長方體，其密度為ρ，檢測質(zhì)量體質(zhì)心坐標為Q(W, D, H)，檢測質(zhì)量體內(nèi)任意一點坐標為P(x, y, z)，檢測質(zhì)量體的長、寬、高分別設(shè)為w、d、h，如圖1所示。

圖1 檢測質(zhì)量體作用于重力梯度儀示意圖Fig.1 Schematic for proof mass’s acting on GGI

根據(jù)萬有引力定律，可得檢測質(zhì)量體對GGI圓盤中心的重力梯度分量計算表達式為：

式中：G為萬有引力常數(shù)。從式(2)可以得到重力梯度與檢測質(zhì)量體的關(guān)系。為了便于說明檢測質(zhì)量體位置對重力梯度的影響，對式(2)變換為：

式中：θ為GGI圓盤中心到檢測質(zhì)量體質(zhì)心矢量在圓盤平面上的投影與圓盤坐標系x軸正方向間的夾角，其中，

從式(3)可以看出，檢測質(zhì)量體相對于重力梯度儀旋轉(zhuǎn)θ角，相當(dāng)于重力梯度儀相對于檢測質(zhì)量體以反方向旋轉(zhuǎn)θ角。設(shè)重力梯度儀本體旋轉(zhuǎn)一周，旋轉(zhuǎn)角度間隔為 10°，在每次旋轉(zhuǎn)位置上測量一次重力梯度值，仿真結(jié)果如圖2所示，從圖2可以看出，通過旋轉(zhuǎn)重力梯度儀本體，環(huán)境物體引起的重力梯度G0以周期性呈現(xiàn)，而重力梯度零位Γ0保持常值，采用這種方法可以分離出環(huán)境梯度 G0和重力梯度零位 Γ0，從而可得重力梯度零位Γ0的計算表達式為：

從式(3)可以進一步得到：當(dāng)θ=0°, 90°, 180°, 360°時，即當(dāng)檢測質(zhì)量體質(zhì)心在x或y坐標軸上時，檢測質(zhì)量體引起的重力梯度分量 Гxy=0E；當(dāng)θ=45°, 135°, 225°, 315°時，即當(dāng)檢測質(zhì)量體質(zhì)心在y=±x直線方程上時，檢測質(zhì)量體引起的重力梯度分量(Гyy-Гxx)=0E；當(dāng)檢測質(zhì)量體質(zhì)心在θ=0°, 90°位置上時，重力梯度分量(Гyy-Гxx)達到最大值，且兩位置處的重力梯度分量(Гyy-Гxx)互為相反數(shù)，因此只要對兩位置處的重力梯度分量(Гyy-Гxx)相減然后取均值，可以消除環(huán)境梯度以及由儀器本身引起的重力梯度零位。

圖2 環(huán)境物體引起的重力梯度Fig.2 Gravity gradient caused by ambient object

同理，檢測質(zhì)量體質(zhì)心在θ=45°, 135°位置上時，重力梯度分量 Гxy達到最大值，且兩位置處的重力梯度分量Гxy互為相反數(shù)。

基于以上分析，給出重力梯度分量(Гyy-Гxx)和Гxy的標定方式：當(dāng)對重力梯度分量(Гyy-Гxx)進行標定時，設(shè)定檢測質(zhì)量體質(zhì)心在x和y坐標軸上；當(dāng)對重力梯度分量Гxy進行標定時，設(shè)定質(zhì)量體質(zhì)心在y=±x直線關(guān)系上，且4個位置上的檢測質(zhì)量體質(zhì)心到GGI圓盤中心的距離相等，檢測質(zhì)量體放置位置如圖3所示。

圖3 檢測質(zhì)量體配置示意圖Fig.3 Setting of proof mass bodies

檢測質(zhì)量體以圖3所示位置配置方式進行放置，根據(jù)重力梯度分量計算式(3)可知，檢測質(zhì)量體在位置1與位置2處引起的重力梯度分量(Гyy-Гxx)大小相等，符號相反；同理，檢測質(zhì)量體在位置3與位置4處引起的重力梯度分量 Гxy大小相等，符號相反，根據(jù)上述關(guān)系給出4個位置檢測質(zhì)量體引起的重力梯度分量測量表達式為

式中：Kyy-xx和Kxy分別表示重力梯度分量(Гyy-Гxx)和Гxy的標度系數(shù)，Gyy-xx和Gxy分別表示由環(huán)境物體引起的重力梯度分量分別表示檢測質(zhì)量體在位置1和位置2實際引起的理論重力梯度分量分別表示檢測質(zhì)量體在位置3和位置4實際引起的理論重力梯度分量Гxy，分別表示重力梯度分量零位。根據(jù)重力梯度零位計算表達式(5)可以得到重力梯度分量(Гyy-Гxx)和Гxy的零位計算表達式為

位置1與位置2處的重力梯度分量(Гyy-Гxx)測量值相減，消除環(huán)境物體和零位引起的重力梯度分量(Гyy-Гxx)，得到如下關(guān)系式：

同理，位置3與位置4處的重力梯度分量Гxy測量值相減，消除環(huán)境物體和零位引起的重力梯度分量Гxy，得到重力梯度分量Гxy的標度系數(shù)計算表達式為

在進行重力梯度測量時，重力梯度儀真正測量的重力梯度信息是GGI周圍環(huán)境物體引起的重力梯度異常，因此重力梯度分量計算表達式為

1）在沒有檢測質(zhì)量體下，旋轉(zhuǎn)重力梯度儀本體整數(shù)倍圈數(shù)，在相同旋轉(zhuǎn)間隔位置處進行重力梯度解調(diào)，利用式(7)計算得到重力梯度零位Γ0；

2）在正交位置處依次放置檢測質(zhì)量體，然后對重力梯度儀輸出信號進行重力梯度解調(diào)，利用式(9)和式(10)計算重力梯度標度系數(shù)K；

3）采用式(11)計算出重力梯度儀真正測量到的重力梯度分量(Гyy-Гxx)和Гxy。

2 半物理仿真實驗與結(jié)果分析

采用上述給出的重力梯度分量零位和標度系數(shù)計算方法，在重力梯度儀半物理仿真分析系統(tǒng)上進行重力梯度標定仿真實驗，仿真參數(shù)如下：加速度計標度因數(shù)為10 mA/g，加速度計對基線距離為0.2 m，檢測質(zhì)量體邊長為0.3 m的均質(zhì)正方體，密度為18 000 kg/m3，檢測質(zhì)量體質(zhì)心到GGI圓盤中心的距離為1 m。根據(jù)重力梯度分量計算式(2)可以計算出檢測質(zhì)量體在位置 1處引起的重力梯度分量(Гyy-Гxx)理論值為-97.53E，在位置3處引起的重力梯度分量Гxy理論值為48.76E。為了模擬環(huán)境物體對重力梯度儀的影響，設(shè)環(huán)境物體邊長為 1 m的均質(zhì)正方體，密度為 2000 kg/m3，環(huán)境物體質(zhì)心位置坐標為(2, 1, 0) m，此環(huán)境物體對GGI引起的重力梯度分量(Гyy-Гxx)和Гxy理論值分別為-21.63E和14.26E。加入2倍頻干擾信號用于模擬重力梯度半物理仿真系統(tǒng)性能引起的重力梯度零位，重力梯度分量(Гyy-Гxx)和 Гxy的零位分別為 200E和10E；檢測質(zhì)量體依次被放置在圖3所示的位置1、2、3、4上，然后分別對重力梯度儀仿真系統(tǒng)輸出信號進行數(shù)據(jù)采集和重力梯度解調(diào)。設(shè)重力梯度儀本體旋轉(zhuǎn)一周，旋轉(zhuǎn)角度間隔為 20°，由環(huán)境物體引起的重力梯度分量輸出信號如圖4所示，從圖4可以看出，重力梯度以近似正弦波形式呈現(xiàn)，這與理論分析結(jié)果一致。采用式(7)計算出半物理仿真系統(tǒng)重力梯度分量(Гyy-Гxx)和 Гxy的零位分別為190.93E和9.98E。檢測質(zhì)量體在不同位置處產(chǎn)生的重力梯度如表1所示，表1中上標t為重力梯度理論值，c為實際測量重力梯度值。根據(jù)表1中的重力梯度分量數(shù)據(jù)，采用式(9)和式(10)計算出重力梯度分量(Гyy-Гxx)和 Гxy的標度系數(shù)分別為0.9578和0.9565。重力梯度零位和標度系數(shù)是重力梯度測量系統(tǒng)固有特性，它只與重力梯度儀本身性能有關(guān)，與外界因素?zé)o關(guān)。

為了檢驗重力梯度零位和標度系數(shù)的正確性，假設(shè)檢驗質(zhì)量體參數(shù)為：質(zhì)量體為均質(zhì)長方體，長、寬和高分別為0.2 m、0.5 m、0.4 m，密度為10 000 kg/m3，其質(zhì)心坐標為(0.8, 0.4, 0) m，此檢驗質(zhì)量體對GGI引起的重力梯度分量(Гyy-Гxx)和 Гxy理論值分別為-68.63E和 39.83E，標定前后重力梯度測量值如表 2所示，從表2的結(jié)果可以得到，重力梯度標定前，由于儀器性能等原因引起的重力梯度被引入到重力梯度信息中，從而造成重力梯度測量誤差，通過對重力梯度標定，標定后重力梯度測量值與重量梯度理論值之間的誤差小于1E，因此提高了重力梯度測量精度。

圖4 半物理仿真系統(tǒng)環(huán)境梯度仿真結(jié)果Fig.4 Gravity gradient caused by ambient object on the Hardware-in-the-loop simulation platform

表1 不同位置處的檢測質(zhì)量體產(chǎn)生的重力梯度Tab.1 Gravity gradient caused by proof mass in different positions

表2 標定前后重力梯度測量值Tab.2 Gravity gradiometry of GGI before and after calibration

3 結(jié) 論

本文從系統(tǒng)級角度提出了一種部分張量旋轉(zhuǎn)加速度計重力梯度儀的標定方法，給出了重力梯度零位和標度系數(shù)的計算表達式以及標定步驟。通過重力梯度半物理仿真系統(tǒng)試驗，仿真結(jié)果表明該方法能夠?qū)崿F(xiàn)重力梯度儀的標定。

（References）：

[1] Dransfield M H, Christensen A N. Performance of airborne gravity gradiometers[J]. The Leading Edge, 2013, 32(8): 908-922.

[2] DiFrancesco D, Grierson A, Kaputa D, et al. Gravity gradiometer systems - advances and challenges[J]. Geophysical Prospecting, 2009, 57(4): 615-623.

[3] Welker T C, Pachter M, Huffman R E. Gravity gradiometer integrated inertial navigation[C]//2013 European Control Conference. Zurich, Switzerland, 2013: 846-851.

[4] Hodges G, Dransfield M H, Shei T C. The falcon airborne gravity gradiometer for engineering applications[C]//23rdSymposium on the Application of Geophysics to Engineering and Environmental Problems 2010. Keystone, USA, 2010: 443-447.

[5] Uzun S. Estimating parameters of subsurface structures from airborne gravity gradiometry data using a Monte-Carlo optimization method[R]. Geodetic Science, The Ohio State University, Columbus, Ohio, USA, 2013.

[6] Rummel R, Yi W, Stummer C. GOCE gravitational gradiometry[J]. Journal of Geodesy, 2011, 85(11): 777-790.

[7] Roberts D, Chowdhury P R, Lowe S J, et al. Airborne gravity gradiometer surveying of petroleum systems under Lake Tanganyika, Tanzania[C]//ASEG-PESA 2015. Perth, Australia, 2015: 1-5.

[8] Christensen A N, Galder C V, Dransfield M. Improved resolution of fixed-wing airborne gravity gradiometer surveys[C]//2014 SEG Annual Meeting. Denver, Colorado, USA, 2014: 1319-1323.

[9] 蔡體菁, 周百令. 重力梯度儀的現(xiàn)狀和前景[J]. 中國慣性技術(shù)學(xué)報, 1999, 7(1): 39-42. Cai Ti-jing, Zhou Bai-ling. Status and prospects of gravity gradiometers[J]. Journal of Chinese Inertial Technology, 1999, 7(1): 39-42.

[10] 劉昱, 張宇, 潘軍, 等. 重力梯度測量中光柵角編碼器安裝誤差的影響及分析[J]. 中國慣性技術(shù)學(xué)報, 2015, 23(1): 6-8. Liu Yu, Zhang Yu, Pan Jun, et al. Influence and analysis of optical angular encoder’s installation error in gravity gradiometry[J]. Journal of Chinese Inertial Technology, 2015, 23(1): 6-8.

[11] 李海兵, 楊慧, 馬存尊, 等. 用于重力梯度儀系統(tǒng)中的高精度微弱信號檢測電路[J]. 中國慣性技術(shù)學(xué)報, 2013, 21(5): 581-584. Li Hai-bing, Yang Hui, Ma Cun-zun, et al. High precision and week signal detection circuits of gravity gradiometer system[J]. Journal of Chinese Inertial Technology, 2013, 21(5): 581-584.

[12] Tu L C, Wang Z W, Liu J Q, et al. Implementation of the scale factor balance on two pairs of quartz-flexure capacitive accelerometers by trimming bias voltage[J]. Review of Scientific Instruments, 2014, 85: 095108-1-7.

Calibration method for gravity gradiometer of rotating accelerometer

QIAN Xue-wu, CAI Ti-jing
(School of Instrument Science and Engineering, Southeast University, Nanjing 210096, China)

A calibration method for rotating accelerometer gravity gradiometer instrument (GGI) is presented. The proof mass bodies are sequentially placed at four positions outside the GGI. The gravity gradient scale coefficient is determined based on the relationship at perpendicular position of gravity gradient component, i.e. the gradients caused by proof masses at orthogonal positions are the same in magnitude and opposite in sign. The gravity gradient bias is determined by the method of rotating the GGI body. The calculation expression of gravity gradient bias and scale coefficient are deduced in details. A validation test of calibration method is performed on the Hardware-in-the-loop simulation platform. Simulation results show that the measurement error of the calibrated gravity gradient is less than 1E.

gravity gradiometer instrument; rotating accelerometer; proof mass; gradient calibration

U666.1

：A

2016-04-27；

：2016-07-08

國家863計劃項目（2011AA060501）

錢學(xué)武（1981—），男，博士研究生，從事精密儀器研究。E-mail: njqxwu@163.com

聯(lián) 系 人：蔡體菁（1961—），男，教授，博士生導(dǎo)師。E-mail: caitij@seu.edu.cn

1005-6734(2016)04-0427-04

10.13695/j.cnki.12-1222/o3.2016.04.002