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

錢學(xué)武　 蔡體菁

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

?

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

錢學(xué)武 蔡體菁

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

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

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

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

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

1　重力梯度信號處理方法

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

圖1　重力梯度信號處理框圖

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

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

(1)

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

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

旋轉(zhuǎn)加速度計重力梯度儀理想輸出信號表達式為

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

(2)

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

(3)

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

1.3小波去噪

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

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

(4)

(5)

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

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

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

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

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

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

圖3　重力梯度信號濾波前后信號波形

圖4　重力梯度信號濾波前后信號頻譜

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

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

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

(b) 重力梯度分量Γxy

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

3　結(jié)語

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

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].物探與化探,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ā)展計劃(863計劃)資助項目(2011AA060501).

10.3969/j.issn.1001-0505.2016.04.006.

U666.1

A

1001-0505(2016)04-0708-05

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