錢(qián)學(xué)武，蔡體菁

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

旋轉(zhuǎn)加速度計(jì)重力梯度儀標(biāo)定方法

錢(qián)學(xué)武，蔡體菁

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

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

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

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

旋轉(zhuǎn)加速度計(jì)重力梯度儀是一種高精密測(cè)量?jī)x器，在出廠(chǎng)前或使用前需要對(duì)其進(jìn)行標(biāo)定。目前，關(guān)于旋轉(zhuǎn)加速度計(jì)重力梯度儀的標(biāo)定方法在公開(kāi)刊物上沒(méi)有報(bào)道。

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

1 重力梯度標(biāo)定方法

重力梯度儀主要有全張量和部分張量?jī)煞N類(lèi)型，針對(duì)全張量重力梯度儀需要對(duì)五個(gè)獨(dú)立重力梯度元素進(jìn)行標(biāo)定，全張量重力梯度儀標(biāo)定方法與部分張量重力梯度儀標(biāo)定方法類(lèi)似，本文只對(duì)部分張量重力梯度儀標(biāo)定方法進(jìn)行分析。對(duì)于部分張量重力梯度儀，需要確定重力梯度分量(Гyy-Гxx)和 Гxy的零位和標(biāo)度系數(shù)，重力梯度儀誤差模型方程為

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

在進(jìn)行重力梯度測(cè)量時(shí)，重力梯度儀真正測(cè)量的重力梯度信息是GGI周?chē)h(huán)境物體引起的重力梯度異常，因此重力梯度分量計(jì)算表達(dá)式為

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

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

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

2 半物理仿真實(shí)驗(yàn)與結(jié)果分析

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

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

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

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

表2 標(biāo)定前后重力梯度測(cè)量值Tab.2 Gravity gradiometry of GGI before and after calibration

3 結(jié) 論

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

（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]. 中國(guó)慣性技術(shù)學(xué)報(bào), 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] 劉昱, 張宇, 潘軍, 等. 重力梯度測(cè)量中光柵角編碼器安裝誤差的影響及分析[J]. 中國(guó)慣性技術(shù)學(xué)報(bào), 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)中的高精度微弱信號(hào)檢測(cè)電路[J]. 中國(guó)慣性技術(shù)學(xué)報(bào), 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

國(guó)家863計(jì)劃項(xiàng)目（2011AA060501）

錢(qián)學(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