錢學(xué)武　 蔡體菁

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

?

旋轉(zhuǎn)加速度計(jì)重力梯度儀數(shù)據(jù)處理方法

錢學(xué)武 蔡體菁

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

為了有效去除旋轉(zhuǎn)加速度計(jì)重力梯度儀輸出信號(hào)中的各種干擾噪聲,提出了一種提取重力梯度信號(hào)的有效方法．首先對(duì)重力梯度儀輸出信號(hào)進(jìn)行故障診斷,然后采用基于Chebyshev最佳一致逼近法設(shè)計(jì)的帶通濾波器對(duì)故障診斷后的信號(hào)進(jìn)行濾波,并對(duì)濾波后的重力梯度信號(hào)進(jìn)行梯度解調(diào),最后采用dmey小波基函數(shù)強(qiáng)制閾值方法對(duì)解調(diào)后的重力梯度信息做進(jìn)一步去噪處理．在重力梯度半物理仿真平臺(tái)上進(jìn)行了仿真試驗(yàn)測(cè)試,結(jié)果表明,所提方法可以有效降低有用信號(hào)均方差,不會(huì)造成數(shù)據(jù)丟失和信號(hào)偏移,提高了重力梯度測(cè)量精度．

重力梯度儀;旋轉(zhuǎn)加速度計(jì);梯度解調(diào);帶通濾波;小波去噪

高精度重力梯度測(cè)量對(duì)于空間科學(xué)、地球科學(xué)、地質(zhì)科學(xué)、能源勘探、慣性導(dǎo)航等領(lǐng)域的研究具有重要意義[1-4]．目前,重力梯度儀主要有旋轉(zhuǎn)加速度計(jì)重力梯度儀、超導(dǎo)重力梯度儀、冷原子重力梯度儀、靜電重力梯度儀,以及基于微機(jī)械加工(MEMS)技術(shù)的重力梯度儀等,其中旋轉(zhuǎn)加速度計(jì)重力梯度儀是目前唯一成功用于機(jī)載/船載動(dòng)基座上的重力梯度測(cè)量儀器[5-8]．

旋轉(zhuǎn)加速度計(jì)重力梯度儀由多種機(jī)械裝置組成,加速度計(jì)性能、加速度計(jì)安裝精度、旋轉(zhuǎn)機(jī)構(gòu)穩(wěn)定性、溫控性能、平臺(tái)穩(wěn)定性、電路系統(tǒng)穩(wěn)定性等因素都會(huì)對(duì)重力梯度信號(hào)造成嚴(yán)重影響[9-10],這些因素會(huì)使重力梯度儀輸出信號(hào)包含各種噪聲,降低了測(cè)量精度,因此研究重力梯度信號(hào)處理方法愈顯重要．國外學(xué)者已對(duì)航空重力梯度儀數(shù)據(jù)處理方法進(jìn)行了相關(guān)研究[11-13],但都是采用基于計(jì)算機(jī)模擬仿真的處理方法,可信度不高．由于傳統(tǒng)數(shù)字濾波器具有平滑效應(yīng),會(huì)對(duì)重力梯度信息起到平滑作用,并且無法濾除重力梯度同頻帶內(nèi)的噪聲,相比而言,小波分析提供了一種自適應(yīng)時(shí)域和頻域同時(shí)局部化的多分辨率分析方法,可以很好地體現(xiàn)出信號(hào)的非平穩(wěn)特性,為信號(hào)處理提供了一種全新途徑[14-16]．本文針對(duì)重力梯度信號(hào)特點(diǎn),提出了一種基于小波分析的去除重力梯度干擾噪聲方法,仿真實(shí)驗(yàn)結(jié)果表明該處理方法可以有效去除重力梯度信號(hào)中的干擾噪聲,提高了重力梯度測(cè)量精度．

1　重力梯度信號(hào)處理方法

在實(shí)際重力梯度測(cè)量過程中,由于儀器性能和外部環(huán)境的干擾,重力梯度儀輸出信號(hào)中往往含有多種較大幅度干擾噪聲,并且存在數(shù)據(jù)丟失或數(shù)據(jù)錯(cuò)誤等情況．針對(duì)重力梯度信號(hào)特點(diǎn),本文提出一種重力梯度信號(hào)有效處理方法:首先對(duì)重力梯度儀輸出信號(hào)進(jìn)行故障診斷處理,消除故障數(shù)據(jù)點(diǎn)；然后采用數(shù)字帶通濾波器對(duì)故障診斷處理后的重力梯度信號(hào)進(jìn)行濾波處理并進(jìn)行重力梯度解調(diào)；最后采用小波分解方法對(duì)含噪重力梯度信息進(jìn)行去噪處理,得到高精度重力梯度信息．重力梯度信號(hào)處理過程如圖1所示．

圖1　重力梯度信號(hào)處理框圖

1.1FIR數(shù)字帶通濾波器

重力梯度信號(hào)信噪比非常低,為了去除噪聲獲得有用信號(hào),需要設(shè)計(jì)合適的帶通濾波器．現(xiàn)有的數(shù)字濾波器設(shè)計(jì)方法中,在保證濾波器所有性能指標(biāo)滿足的前提下,采用Chebyshev最佳一致逼近方法設(shè)計(jì)的數(shù)字濾波器階數(shù)最小．理論上,FIR濾波器的通帶帶寬越窄,階數(shù)越高,濾波效果越好,但濾波器階數(shù)太高會(huì)出現(xiàn)濾波效果不明顯,甚至出現(xiàn)濾波不穩(wěn)定現(xiàn)象．濾波器階數(shù)可根據(jù)Kaiser[17]提出的經(jīng)驗(yàn)公式進(jìn)行近似估計(jì):

(1)

式中,N為濾波器階數(shù);過渡帶帶寬Δf=(ωs-ωp)/(2π),ωp和ωs分別為濾波器通帶和阻帶截止角頻率;δp和δs分別為濾波器通帶和阻帶的紋波峰值．本文采用Remez算法設(shè)計(jì)濾波器,根據(jù)重力梯度信號(hào)特征,并反復(fù)進(jìn)行實(shí)驗(yàn)測(cè)試驗(yàn)證,最終得到一組濾波效果較好的帶通濾波器參數(shù)．本文選取的濾波器參數(shù)為:中心頻率為0.5Hz,通帶帶寬為0.02Hz,過渡帶帶寬為0.01Hz,濾波器階數(shù)為800．

1.2重力梯度解調(diào)

旋轉(zhuǎn)加速度計(jì)重力梯度儀理想輸出信號(hào)表達(dá)式為

Eout=2RKKI[(Γyy-Γxx)sin2ωt+2Γxycos2ωt]

(2)

式中,Eout為重力梯度儀輸出信號(hào);R為加速度計(jì)質(zhì)量中心到GGI圓盤中心的距離;K為重力梯度信號(hào)放大增益;KI為加速度計(jì)標(biāo)度因數(shù);ω為GGI圓盤旋轉(zhuǎn)角頻率;Γyy,Γxx和Γxy分別為相應(yīng)方向上的重力梯度張量元素．只要對(duì)式(2)在2ω頻率上進(jìn)行幅值解調(diào)就可以得到重力梯度分量Γyy-Γxx和Γxy,解調(diào)公式為

(3)

式中, T為重力梯度儀圓盤旋轉(zhuǎn)周期．

1.3小波去噪

小波分析屬于時(shí)頻分析的一種,它能夠同時(shí)給出信號(hào)的時(shí)域和頻域特征,有效區(qū)分出信號(hào)中的突變成分和干擾噪聲,從而實(shí)現(xiàn)信號(hào)分離和降噪．設(shè)重力梯度信號(hào)為Γ(t),被噪聲污染的重力梯度信號(hào)為G(t),重力梯度信號(hào)模型可表示為

G(t)=Γ(t)+e(t)

(4)

(5)

式中,j,k分別為離散化的縮放因子和平移因子;ψ(t)為小波基函數(shù);ψ*(t)為ψ(t)的共軛．利用式(5)可得到信號(hào)G(t)在第j尺度上的小波系數(shù)．采用Mallat算法對(duì)信號(hào)G(t)進(jìn)行小波分解和信號(hào)重構(gòu)．

重力梯度信號(hào)本質(zhì)上是一種非平穩(wěn)低頻隨機(jī)信號(hào),經(jīng)過小波分解后,在某個(gè)分解尺度上的小波系數(shù)具有較大幅值,而噪聲頻帶比重力梯度信號(hào)頻帶大得多．針對(duì)這種情況,在重力梯度數(shù)據(jù)去噪階段采用強(qiáng)制閾值去噪方法．該方法僅對(duì)某個(gè)尺度下的低頻系數(shù)進(jìn)行小波變換來重建信號(hào),直接舍棄其他尺度上的系數(shù)分量,從而達(dá)到信號(hào)降噪的目的．在強(qiáng)制閾值濾波過程中,采用dmey小波基函數(shù),這是由于該小波基具有很好的正則性、緊支撐性和平滑性,并且適合重力梯度信號(hào)的特點(diǎn)．對(duì)重力梯度信號(hào)進(jìn)行10層分解,經(jīng)過比較相鄰層信號(hào)與真實(shí)信號(hào)的變化趨勢(shì)發(fā)現(xiàn),第6層低頻信號(hào)與真實(shí)信號(hào)比較接近,其他層與真實(shí)信號(hào)相差較大,故選擇第6層低頻信號(hào)作為重力梯度信息的最終處理結(jié)果．

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

根據(jù)旋轉(zhuǎn)加速度計(jì)重力梯度儀工作原理,采用高性能計(jì)算機(jī)、可編程高精度電流源、低噪聲電流放大器、低噪聲電壓放大器、多路切換開關(guān)和高精度數(shù)字電壓表等高精密儀器組建了重力梯度半物理仿真系統(tǒng)[18]．利用該仿真系統(tǒng)可以完成信號(hào)產(chǎn)生、信號(hào)放大、信號(hào)濾波、重力梯度解調(diào)等仿真實(shí)驗(yàn)．仿真系統(tǒng)采用LabVIEW軟件開發(fā)．重力梯度半物理仿真系統(tǒng)實(shí)物圖如圖2所示．

圖2　重力梯度半物理仿真系統(tǒng)實(shí)物圖

采用高密度的鎢合金均質(zhì)圓球作為重力梯度引力裝置,利用所開發(fā)的重力梯度半物理仿真系統(tǒng)對(duì)重力梯度信號(hào)進(jìn)行仿真測(cè)試. 實(shí)驗(yàn)參數(shù)如下:圓球密度為18 000kg/m3,質(zhì)量為1 000kg,圓球質(zhì)心坐標(biāo)為(0.8,0,0)m,圓盤旋轉(zhuǎn)周期為4s,加速度計(jì)對(duì)的基線距離為0.2m,加速度計(jì)標(biāo)度因數(shù)為10mA/g,信號(hào)放大增益為5×106,信號(hào)采樣頻率為10Hz．根據(jù)萬有引力定律，可以計(jì)算出在上述參數(shù)下圓球引起的重力梯度分量理論值分別為Γyy-Γxx=-393E,Γxy=0E．將從旋轉(zhuǎn)加速度計(jì)重力梯度儀樣機(jī)上采集到的原始加速度計(jì)數(shù)據(jù)作為重力梯度半物理仿真系統(tǒng)信號(hào)源,對(duì)重力梯度半物理仿真系統(tǒng)輸出信號(hào)進(jìn)行數(shù)據(jù)采集和數(shù)據(jù)故障診斷處理,然后采用帶通濾波器對(duì)故障診斷后的重力梯度信號(hào)進(jìn)行濾波．濾波前后的重力梯度信號(hào)時(shí)域波形如圖3所示,對(duì)應(yīng)的頻譜如圖4所示．

圖3　重力梯度信號(hào)濾波前后信號(hào)波形

圖4　重力梯度信號(hào)濾波前后信號(hào)頻譜

從圖3可看出,重力梯度信號(hào)全部被噪聲淹沒,完全看不到重力梯度信號(hào),濾波后的信號(hào)幅值大幅降低．從圖4可看出,0.5Hz的重力梯度信號(hào)較為明顯,除0.5Hz外的低高頻干擾諧波幅度較大,而帶通濾波后的重力梯度信號(hào)頻譜只含有0.5Hz窄帶信號(hào)和通帶帶寬內(nèi)殘留的少量噪聲信號(hào),0.5Hz外的高幅值低高頻噪聲基本被濾除．利用式(3)對(duì)帶通濾波后的重力梯度信號(hào)進(jìn)行梯度解調(diào),然后對(duì)解調(diào)得到的重力梯度信息進(jìn)行100s平滑濾波、200s平滑濾波和小波強(qiáng)制閾值去噪處理,不同處理方法的結(jié)果如表1所示,波形對(duì)比如圖5所示．

從表1和圖5可看出,采用平滑濾波處理方法

表1　重力梯度不同處理方法結(jié)果對(duì)比　E

(b) 重力梯度分量Γxy

雖然可以有效去除噪聲,但隨著平滑時(shí)間的延長,重力梯度信息丟失越多,且平滑濾波后的均值發(fā)生了明顯偏移．而采用小波去噪方法不會(huì)出現(xiàn)數(shù)據(jù)丟失和均值偏移現(xiàn)象,并且可以明顯降低重力梯度信息均方差．同時(shí)還發(fā)現(xiàn),重力梯度分量Γyy- Γxx和Γxy均值分別約為-430和-22E,而理論值分別為-393和0E,重力梯度測(cè)量值與理論值相差較大,這主要是由于在0.5Hz上存在信號(hào)串?dāng)_.造成信號(hào)串?dāng)_的原因很多,可通過采用重力梯度標(biāo)定方法進(jìn)行重力梯度補(bǔ)償解決．

3　結(jié)語

本文針對(duì)重力梯度信號(hào)特點(diǎn),提出了一種去除重力梯度干擾噪聲的有效處理方法.首先對(duì)重力梯度信號(hào)進(jìn)行故障診斷處理,然后采用帶通濾波器對(duì)故障診斷后的信號(hào)進(jìn)行濾波,此濾波方法可以把重力梯度信號(hào)頻率外的低高頻干擾噪聲濾除,提高了重力梯度信號(hào)信噪比;接著對(duì)濾波后的重力梯度信號(hào)進(jìn)行梯度解調(diào);最后采用小波去噪方法對(duì)解調(diào)結(jié)果做進(jìn)一步處理．為了驗(yàn)證所提處理方法的有效性,在重力梯度半物理仿真系統(tǒng)上進(jìn)行了仿真實(shí)驗(yàn)測(cè)試．仿真結(jié)果表明,小波去噪方法可以顯著提高重力梯度測(cè)量精度,并且不會(huì)造成數(shù)據(jù)丟失和信號(hào)偏移,為實(shí)現(xiàn)高精度重力梯度測(cè)量提供了一種方法．

References)

[1]Difrancesco D. Advances and challenges in the development and deployment of gravity gradiometer systems[C/OL]//EGM2007InternationalWorkshop. Capri, Italy, 2007. http://www.earthdoc.org/publication/publicationdetails/?publication=41199.

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

[3]Christensen A N, Galder C V, Dransfield M. Improved resolution of fixed-wing airborne gravity gradiometer surveys[C]//2014SEGAnnualMeeting. Denver, Colorado, USA, 2014: 1319-1323. DOI: 10.1190/segam2014-0713.1.

[4]DiFrancesco D, Meyer T, Christensen A, et al. Gravity gradiometry—Today and tomorrow[C]//11thSAGABiennialTechnicalMeetingandExhibition. Swaziland, 2009: 80-83.

[5]Hodges G, Dransfield M H, Shei T C. The Falcon airborne gravity gradiometer for engineering applications[C]//SymposiumontheApplicationofGeophysicstoEngineeringandEnvironmentalProblems2010. Keystone, Colorado, USA, 2010: 443-447. DOI:10.4133/1.3445467.

[6]Christensen A N, Hodges G. HeliFALCON?airborne gravity gradiometer data acquisition in rugged terrain[C]//Proceedingsofthe11thSEGJInternationalSymposium. Yokohama, Japan, 2013:140-145. DOI:10.1190/segj112013-036.

[7]Difrancesco D, Grierson A, Kaputa D, et al. Gravity gradiometer systems—Advances and challenges[J].GeophysicalProspecting, 2009, 57(4): 615-623.

[8]Dransfield M H, Christensen A N. Performance of airborne gravity gradiometers[J].TheLeadingEdge, 2013, 32(8):908-922. DOI:10.1190/tle32080908.1.

[9]Jekeli C. The gravity gradiometer survey system (GGSS)[J].EosTransactionsAmericanGeophysicalUnion, 1988, 69(8):105, 116-117.

[10]Jekeli C. Statistical analysis of moving-base gravimetry and gravity gradiometry[R]. Columbus,Ohio,USA: Geodetic and GeoInfomation Science, The Ohio State University, 2003.

[11]Cesar J, Lyrio S O. Wavelet denoising of gravity gradiometry data[C]//SEGInternationalExpositionand71stAnnualMeeting. San Antonio, USA, 2001:1474-1477. DOI:10.1190/1.1816384.

[12]Carlos D U, Braga M A, Galbiatti H F, et al. Airborne gravity gradiometry-data processing and interpretation [J].RevistaBrasileiraDeGeofǐsica, 2013, 31(3):427-453.

[13]Christensen A N, Dransfield M H, Galder C V. Noise and repeatability of airborne gravity gradiometry [J].FirstBreak, 2015, 33:55-63.

[14]Barnes G, Lumley J. Processing gravity gradient data [J].Geophysics, 2011, 76(2):133-147.

[15]Jekeli C, Abt T L. The statistical performance of the matched filter for anomaly detection using gravity gradients[R]. Columbus, Ohio, USA: Division of Geodetic Science, Ohio State University, 2010.

[16]de Oliveira Lyrio J C S, Tenorio L, Li Y, et al. Efficient automatic denoising of gravity gradiometry data[J].Geophysics, 2004, 69(3):772-782. DOI:10.1190/1.1759463.

[17]Kaiser J F. Nonrecursive digital filter design using the Io-sinh window function[C]//IEEEInternationalSymposiumonCircuits&Systems. San Francisco, CA, USA, 1974:20-23.

[18]蔡體菁,錢學(xué)武,丁昊.旋轉(zhuǎn)加速度計(jì)重力梯度儀重力梯度信號(hào)仿真[J].物探與化探,2015,39(S1):76-79.

Cai Tijing, Qian Xuewu, Ding Hao. Signal simulation of gravity gradiometer of rotating accelerometer[J].GeophysicalandGeochemicalExploration, 2015, 39(S1):76-79. (in Chinese)

DOI:10.3969/j.issn.1001-0505.2016.04.007

Data processing method for rotating accelerometer gravity gradiometer

Qian Xuewu Cai Tijing

(School of Instrument Science and Engineering, Southeast University, Nanjing 210096, China)

In order to remove the noise in the output signals of the rotating accelerometer gravity gradiometer instrument (GGI), an effective method for extracting the gravity gradient signal is presented. First, the fault diagnosis is carried out to preprocess the output signals of the GGI. Then a band-passing filter based on the Chebyshev optimal uniform approximation method is used to filtrate the diagnosed signals, and the output signals of the band-pass filter are demodulated in the next step. Finally, the gravity gradient is further denoised by the compulsory threshold method based on the demy wavelet function. The proposed method is tested on the gravity gradient hard-in-the-loop simulation platform. The simulation results show that the method can effectively reduce the mean square deviation of the useful signals and avoid data loss or signal offset. Therefore, the method can improve the precision of the gravity gradiometer.

gravity gradiometer; rotating accelerometer; gradient demodulation; band-passing filtering; wavelet de-noising

10.3969/j.issn.1001-0505.2016.04.006

2015-12-28.作者簡介: 錢學(xué)武(1981—)，男，博士生；蔡體菁(聯(lián)系人)，男，博士，教授，博士生導(dǎo)師，caitij@seu.edu.cn.

國家高技術(shù)研究發(fā)展計(jì)劃(863計(jì)劃)資助項(xiàng)目(2011AA060501).

10.3969/j.issn.1001-0505.2016.04.006.

U666.1

A

1001-0505(2016)04-0708-05

引用本文: 錢學(xué)武,蔡體菁.旋轉(zhuǎn)加速度計(jì)重力梯度儀數(shù)據(jù)處理方法[J].東南大學(xué)學(xué)報(bào)(自然科學(xué)版),2016,46(4):708-712.