楊 曄，李 達(dá)，高 巍

（天津航海儀器研究所，天津 300131）

旋轉(zhuǎn)加速度計式重力梯度儀輸出解調(diào)與濾波

楊 曄，李 達(dá)，高 巍

（天津航海儀器研究所，天津 300131）

重力梯度儀是對地球表面微小重力梯度變化進(jìn)行連續(xù)測量的儀器。由于核心敏感元件加速度計工藝與性能水平限制，以及多環(huán)節(jié)安裝誤差等因素導(dǎo)致系統(tǒng)實際輸出信號中包含了大量噪聲，且信噪比極低，為了能夠在強(qiáng)噪聲中有效提取真實的重力梯度信號，需在信號解調(diào)過程中降低諧波干擾引起的測量偏差。結(jié)合誤差產(chǎn)生機(jī)理，分析比較了不同的解調(diào)方法對重力梯度信號解調(diào)的影響，明確了合理的信號解調(diào)手段。同時，在數(shù)據(jù)處理過程中增加帶通濾波環(huán)節(jié)，進(jìn)一步降低動態(tài)噪聲對系統(tǒng)的影響。結(jié)合動態(tài)搖擺實驗數(shù)據(jù)的仿真驗證結(jié)果表明：上述措施對重力梯度儀原始觀測數(shù)據(jù)的噪聲抑制具有明顯效果，測量偏差在一定程度上得到解決，解調(diào)后信號的噪聲幅度下降至20%，提高了系統(tǒng)的空間分辨率，具有工程化應(yīng)用意義。

重力梯度儀；信號提取；帶通濾波；同步解調(diào)

重力梯度儀是用于測量重力場分布梯度的精密設(shè)備，基于Bell Aerospace公司提出的旋轉(zhuǎn)加速度計測量原理的重力梯度儀（后文中簡稱梯度儀）是迄今唯一實用的近地表動態(tài)測量重力梯度儀[1-3]。礙于重力梯度信號的微弱和加速度計性能的限制，梯度儀輸出信號的信噪比極低，需后期數(shù)據(jù)處理進(jìn)一步提高測量精度[4-6]。由于世界范圍內(nèi)能夠從事該技術(shù)研究的國家和單位很少，相關(guān)技術(shù)文獻(xiàn)鮮有報道，本文結(jié)合科研工作和實測數(shù)據(jù)來驗證方法的有效性，優(yōu)化重力梯度信號數(shù)據(jù)處理過程。

1 重力梯度測量原理與存在問題

1.1 測量原理

如圖1所示，旋轉(zhuǎn)加速度計式重力梯度儀基于加速度計位置差分測量原理，通過機(jī)械旋轉(zhuǎn)的方式將重力梯度張量分量調(diào)制到系統(tǒng)旋轉(zhuǎn)頻率的二倍頻處。加速度計輸出與重力梯度張量分量之間的關(guān)系可表示為

式中：a1、a2、a3、a4是4個加速度計的輸出；R是加速度計檢測質(zhì)心到旋轉(zhuǎn)中心的距離，Γxx、Γyy、Γxy是重力梯度張量分量；ω是旋轉(zhuǎn)裝置的旋轉(zhuǎn)角速度。

圖1 旋轉(zhuǎn)加速度計式重力梯度儀測量原理Fig.1 Measurement principle of rotating accelerometer gravity gradiometer

以往的信號解調(diào)方案如圖2所示，其中方波解調(diào)公式為

式中：Aout(i)是經(jīng)位置同步的加速度計輸出采樣值；N是重力梯度儀圓盤旋轉(zhuǎn)半周時采樣數(shù)據(jù)個數(shù)；R是旋轉(zhuǎn)半徑；π2是將解調(diào)值轉(zhuǎn)化為幅值的比例系數(shù)。

圖2 原始重力梯度數(shù)據(jù)處理框圖Fig.2 Flow chart for processing the original gravity gradiometer data

1.2 存在問題

梯度儀在進(jìn)行水平搖擺試驗（搖擺幅度為2°，搖擺周期為20 s）時的典型解調(diào)出的信號如圖3所示，可以看出在4 892 s由靜態(tài)開始搖擺時，信號均值出現(xiàn)600E(1E=10-9s-2)的偏差“臺階”，并且原始信號信噪比極低，噪聲峰峰值高達(dá)10000E，相關(guān)數(shù)據(jù)見表1。

圖3 搖擺試驗時梯度儀Γxx-Γyy解調(diào)信號（方波解調(diào)）Fig.3 Demodulation signal of gravity gradiometer Γxx-Γyyin swing test (square-wave demodulation)

表1 搖擺試驗梯度信號Γxx-Γyy的均值和標(biāo)準(zhǔn)差(方波解調(diào))Tab.1 Mean values and standard deviations of gravity gradient signals in swing test (square-wave demodulation)

2 梯度儀數(shù)據(jù)處理方法

在此分析梯度儀在搖擺試驗中出現(xiàn)的測量偏差的機(jī)理，系統(tǒng)可采用的解調(diào)方式主要有方波解調(diào)和正弦解調(diào)兩種，其中正弦解調(diào)的公式為

式中：ts是梯度儀的采樣時間間隔；分子上的2是將解調(diào)值轉(zhuǎn)化為幅值的比例系數(shù)。兩種解調(diào)方式在解調(diào)2ω頻率上的信號幅頻響應(yīng)對比如圖4所示。

圖4 方波解調(diào)和正弦解調(diào)方式幅頻響應(yīng)對比Fig.4 Comparison on amplitude-frequency responses between square-wave demodulation and sinusoidal demodulation

從圖4可以看出，在解調(diào)頻點附近，兩種解調(diào)方式效果相當(dāng)，但在解調(diào)頻率的奇數(shù)倍頻處，方波解調(diào)方式會引入高頻諧波項，而正弦解調(diào)方式的結(jié)果不被高頻諧波項影響。引入一個模擬梯度信號y= sin(2πf1t)+sin(2πf2t )進(jìn)行仿真驗證（其中，f1=2ω，f2=6ω），該信號時域圖見圖5（T=12ω），解調(diào)結(jié)果對比見表2。

圖5 模擬梯度信號時域圖Fig.5 Time domain of simulated gravity gradiometer signal

表2 兩種解調(diào)方式模擬梯度信號解調(diào)結(jié)果對比Tab.2 Comparison of two demodulation methods for simulated gradient signals

由此可知，方波解調(diào)值會以一定比例引入高頻諧波分量，由于搖擺試驗中在機(jī)械旋轉(zhuǎn)調(diào)制頻率的高頻諧波頻點處存在較大噪聲耦合，方波解調(diào)的方式會將這些諧波頻點噪聲轉(zhuǎn)化為解調(diào)值，匯入梯度信號，導(dǎo)致梯度信號出現(xiàn)測量偏差，而正弦解調(diào)值不會受到高頻諧波噪聲的影響。從圖4可以看出，解調(diào)頻點附近的噪聲也會對解調(diào)結(jié)果產(chǎn)生影響，如果原始數(shù)據(jù)存在頻率差為Δf，幅值為A的噪聲項，其對解調(diào)值的影響為

式中：k為解調(diào)值序號。由式(5)可得解調(diào)頻點附近的噪聲會使解調(diào)結(jié)果出現(xiàn)等幅的余弦震蕩，增加梯度信號噪聲，因此需在解調(diào)前加入帶通濾波環(huán)節(jié)以提高解調(diào)精度。并且由式(4)可知，頻率差為Δf噪聲項會引起解調(diào)結(jié)果呈周期為1(2Δf)的余弦振蕩，因此解調(diào)后對梯度信號低通濾波器截止頻率的選擇與原始信號帶通濾波器的帶寬參數(shù)息息相關(guān)。如果帶通濾波器的帶寬參數(shù)為f，解調(diào)后對梯度信號的低通濾波器截止頻率必須低于f才有降噪效果。在實際工作中可結(jié)合梯度儀自身性能和外部環(huán)境噪聲水平綜合考慮是否對梯度信號進(jìn)行低通濾波。

3 實驗數(shù)據(jù)處理與結(jié)果分析

將同樣的搖擺試驗數(shù)據(jù)進(jìn)行正弦解調(diào)，梯度信號如圖6所示，相關(guān)數(shù)據(jù)見表3，從中可以看出，梯度信號在水平搖擺試驗時的測量偏差可在一定程度上得到解決，偏差幅度從600 E降至20 E。

圖6 搖擺試驗梯度儀Γxx-Γyy解調(diào)信號(正弦解調(diào)）Fig.6 Demodulation signal of gravity gradiometer Γxx-Γyyin swing test (sinusoidal demodulation)

表3 搖擺試驗梯度信號Γxx-Γyy的均值和標(biāo)準(zhǔn)差(正弦解調(diào))Tab.3 Mean values and standard deviations of gravity gradient signals in swing test (sinusoidal demodulation)

為進(jìn)一步驗證正弦解調(diào)方式的合理性，分別將兩種解調(diào)方式的解調(diào)值與原始數(shù)據(jù)在2ω(Hz)處的FFT（快速傅式變換）值作比較，其結(jié)果見表4。從表4中可以看出，方波解調(diào)得到的梯度信號與原始數(shù)據(jù)在2ω(Hz)處的FFT值相差較大，是由于引入的高頻諧波噪聲耦合項不穩(wěn)定導(dǎo)致的，而正弦解調(diào)的信號與FFT值一致，因此使用正弦解調(diào)方式得到的梯度信號更加接近真實值。

再在解調(diào)前對原始數(shù)據(jù)增加帶通濾波環(huán)節(jié)，根據(jù)梯度儀輸出信號特征并在考慮儀器測量信息空間分辨率的基礎(chǔ)上，在采樣頻率800 Hz的條件下設(shè)計一組合理的帶通濾波器參數(shù)：中心頻率為2ω(Hz)，窗口函數(shù)為凱塞窗，帶通帶寬為0.01 Hz，過渡帶帶寬為0.005 Hz，濾波器階數(shù)為64 000。將濾波后數(shù)據(jù)解調(diào)得到梯度信號如圖7所示，相關(guān)數(shù)據(jù)見表5?？梢妿V波后梯度信號降噪效果明顯，梯度信號噪聲峰峰值降為1000E左右，標(biāo)準(zhǔn)差降為原來的20%左右。

表4 搖擺試驗梯度儀輸出信號兩種解調(diào)方式在2ω(Hz)處FFT值對比Tab.4 Comparison on FFT values of gradient signals at 2ω(Hz) by the two demodulation modes in swing test

圖7 帶通濾波后梯度張量Γxx-Γyy解調(diào)信號Fig.7 Band-pass filtered gravity gradient tensor Γxx-Γyydemodulation signal

表5 使用帶通濾波器后解調(diào)效果對比(未平滑處理)Tab.5 Comparison on demodulation effects after band-pass filtering (no smoothing)

在梯度儀實際工作中，即使在解調(diào)前增加帶通濾波環(huán)節(jié)，梯度信號噪聲仍然很大，仍需對其進(jìn)行平滑濾波處理。將通過帶通濾波器后解調(diào)得到的梯度信號與原始梯度信號分別使用 N秒平滑濾波并計算標(biāo)準(zhǔn)差，對比結(jié)果如圖8所示，相關(guān)數(shù)據(jù)見表6。

從圖8中可以看出，在使用帶通濾波處理后對解調(diào)出來的梯度信號使用100 s以下的平滑濾波處理，其對噪聲的抑制效果不明顯，如需進(jìn)一步降噪需進(jìn)行更長時間的平滑濾波處理，這是由于帶通濾波器的通帶帶寬參數(shù)為0.01 Hz造成的。從表6中還可以看出，僅使用帶通濾波處理后，解調(diào)的梯度信號均值會發(fā)生幅值為4E的數(shù)據(jù)偏移現(xiàn)象，但再經(jīng)過平滑濾波后偏移現(xiàn)象消失，這是由于在2ω(Hz)上存在信號串?dāng)_，這里可以對梯度信號使用平滑濾波的手段消除。

圖8 梯度信號N秒平滑曲線對比Fig.8 Comparison on curves of gradient signals before and after smoothing for N seconds

表6 重力梯度不同數(shù)據(jù)處理方法結(jié)果對比Tab.6 Comparison on different data processing methods for gravity gradient data

改進(jìn)后的重力梯度信號數(shù)據(jù)處理過程如圖 9所示。通過增加帶通濾波環(huán)節(jié)的手段可將梯度信號噪聲降至原來的20%左右，達(dá)到提高梯度儀測量信息空間分辨率的目的。將解調(diào)方式改為正弦解調(diào)可以有效地消除測量偏差，提高儀器測量分辨率。

圖9 改進(jìn)后重力梯度數(shù)據(jù)處理框圖Fig.9 Flow chart for processing the improved gravity gradiometer data

4 總 結(jié)

本文針對梯度儀輸出信號特征，提出一種重力梯度信號提取的有效方法：

1）使用帶通濾波器對梯度儀原始輸出信號進(jìn)行濾波，可將重力梯度信號頻點外的干擾噪聲濾掉；

2）對濾波后的信號用正弦解調(diào)代替方波解調(diào)，以消除諧波頻點處的噪聲引起的測量偏差；

3）將解調(diào)后的梯度信號進(jìn)行平滑濾波處理，進(jìn)一步來提高測量分辨率。

通過采用上述方法對梯度儀原始測量數(shù)據(jù)進(jìn)行離線處理，測量偏差可在一定程度上得到解決，同時信號噪聲標(biāo)準(zhǔn)差可有效降低到原來的20%，噪聲峰峰值降低到1000E，提高了梯度儀的測量精度和測量信息的空間分辨率，為實現(xiàn)高精度動態(tài)重力梯度測量提供了一種數(shù)據(jù)處理的方法。

（References）：

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

[2] DiFrancesco D, Meyer T, Christensen A, et al. Gravity gradiometry—today and tomorrow[C]//11thSAGA Biennial Technical Meeting and Exhibition. Swaziland, 2009: 80-83.

[3] Dransfield M H, Christensen A N. Performance of airbornegravity gadiometers[J]. The Leading Edge, 2013, 32(8): 908-922.

[4] Lee J B. FALCON gravity gradiometer technology[J]. Exploration Geophysics, 2001, 32(4): 247-250.

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

[6] 蔡體菁, 錢學(xué)武, 丁昊. 旋轉(zhuǎn)加速度計重力梯度儀重力梯度信號仿真[J]. 物探與化探, 2015, 39(S1): 76-79. Cai Ti-jing, Qian Xue-wu, Ding Hao. Signal simulation of gravity gradiometer of rotating accelerometer[J]. Geophysical and Geochemical Exploration, 2015, 39(S1): 76-79.

[7] 劉昱, 張宇, 潘軍, 等. 重力梯度測量中光柵角編碼器安裝誤差的影響及分析[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.

[8] O’keefe G J, Lee J B, Turner R J, et al. Gravity gradiometer. US: No.5922951[P]. 1999-10-05.

[9] Van Leeuwen E H, McCracken K G, Lee J B. Airborne gravity gradiometers. US: No.6883372 B2[P]. 2005.

[10] 錢學(xué)武, 蔡體菁, 尹航. 一種旋轉(zhuǎn)加速度計重力梯度儀重力梯度解調(diào)方法[J]. 物探與化探, 2015, 39(S1): 80-83. Qian Xue-wu, Cai Ti-jing, Yin Hang. A gravity gradient demodulation method for rotation accelerometer gravity gradiometer[J]. Geophysical and Geochemical Exploration, 2015, 39(S1): 80-83.

[11] Moody M V, Paik H J. A superconducting gravity gradiometer for inertial navigation[C]//2004 IEEE Position Location and Navigation Symposium. 2004: 775-781.

[12] 錢學(xué)武, 蔡體菁. 旋轉(zhuǎn)加速度計重力梯度儀加速度計標(biāo)度因數(shù)實時反饋調(diào)整方法[J]. 中國慣性技術(shù)學(xué)報, 2016, 24(2): 148-153. Qian Xue-wu, Cai Ti-jing. Real-time feedback adjusting methods for accelerometer scale factor in gravity gradiometer for rotating accelerometer[J]. Journal of Chinese Inertial Technology, 2016, 24(2): 148-153.

Gradient signal demodulation and filtering for rotating accelerometer gravity gradiometer

YANG Ye, LI Da, GAO Wei
(Tianjin Navigation Instrument Research Institute, Tianjin 200131, China)

In rotating accelerometer gravity gradiometer, due to the limitation of the accelerometer’s process and performance level and the effects of various installation errors, there are lots of noises in the system output signal, and the SNRs are extremely low. In order to effectively extract the real gravity gradient signals in strong noises, the measurement errors caused by harmonic interference must be reduced in the process of signal demodulation. Combined with the error generation mechanism, the influences of different demodulation methods on the gravity gradient signal demodulation are analyzed and compared, and a reasonable method of signal demodulation is chosen. Meanwhile, a band-pass filter is added in the data processing process to further reduce the influence of the dynamic noise. Combined with dynamic swing test data, the simulations show that the above measures have obvious effect on the noise suppression of the original gravity data, and the noise amplitude of the demodulated signal is reduced to 20%, which improves the spatial resolution of the system, and has the engineering application value.

gravity gradiometer instrument; signal extraction; band-pass filter; synchronous demodulation

U666.1

：A

2016-10-02；

：2016-11-20

國家高技術(shù)研究發(fā)展計劃（863計劃）（2011AA060501）；總裝預(yù)研項目（51309040702）

楊曄（1968—），男，研究員，碩士生導(dǎo)師，從事導(dǎo)航、制導(dǎo)與控制方面研究工作。E-mail: Liuyuyangye@eyou.com

1005-6734(2016)06-0701-05

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