田顏鋒 李?yuàn)檴?肖 凡

航空重力測(cè)量水平加速度改正的小波預(yù)處理*

田顏鋒1)李?yuàn)檴?)肖 凡3)

根據(jù)航空重力測(cè)量中水平加速度改正的頻譜特性和小波變換頻帶劃分規(guī)則，導(dǎo)出了水平加速度改正的小波預(yù)濾波尺度。以實(shí)測(cè)數(shù)據(jù)為例，基于信號(hào)小波變換尺度間的相關(guān)性，分別利用Haar小波、db4小波、db6小波和sym6小波對(duì)其進(jìn)行小波預(yù)濾波處理，并與目前常用的FIR低通預(yù)濾波處理結(jié)果進(jìn)行比對(duì)，結(jié)果表明:小波預(yù)濾波處理水平加速度改正的精度較FIR低通預(yù)濾波處理得到的精度高，并且用db系列小波進(jìn)行預(yù)處理得到的結(jié)果精度更加均勻，據(jù)此提出對(duì)水平加速度改正預(yù)處理時(shí)應(yīng)選用db系列小波或者支集較長(zhǎng)的sym8小波。

航空重力測(cè)量;預(yù)濾波;水平加速度改正;厄特弗斯改正;傅里葉變換;小波變換

1 引言

我國(guó)的航空重力測(cè)量系統(tǒng)CHAGS(Chinese Airborne Gravity System)采用基于穩(wěn)定平臺(tái)的航空重力儀，穩(wěn)定平臺(tái)維持其傳感器的垂直指向［1］。當(dāng)穩(wěn)定平臺(tái)保持當(dāng)?shù)厮綍r(shí)，重力儀指向當(dāng)?shù)劂U垂線方向，重力測(cè)量不受水平加速度的影響;當(dāng)穩(wěn)定平臺(tái)不能保持當(dāng)?shù)厮綍r(shí)，重力測(cè)量受到飛行載體水平加速度的影響從而產(chǎn)生水平加速度改正。它是航空重力測(cè)量數(shù)據(jù)處理中必不可少的環(huán)節(jié)，由于改正模型中二次項(xiàng)的出現(xiàn)，改變了整個(gè)模型的線性特征，因此必須對(duì)原始數(shù)據(jù)進(jìn)行預(yù)濾波處理。石磐［2］利用FIR低通濾波器對(duì)水平加速度改正進(jìn)行預(yù)濾波處理，估算出其平均量級(jí)約為3×10-5ms-2;孫中苗［3］論證了計(jì)算水平加速度改正的兩步法和一步法之間的關(guān)系，根據(jù)平臺(tái)傾角的頻譜特性提出水平加速度改正的基本濾波尺度為10 s，利用該尺度得到的重力測(cè)量結(jié)果與地面向上延拓參考值之間的系統(tǒng)差優(yōu)于0.5×10-5ms-2。但是FIR濾波器不能有效抑制混入信號(hào)頻帶內(nèi)部的噪聲分量，由此，本文將基于航空重力測(cè)量減震系統(tǒng)得到水平加速度改正的小波預(yù)處理尺度，并對(duì)其進(jìn)行小波處理。

2 小波變換的基本原理

小波被定義為積分為零且平方可積的函數(shù)，滿足［4-7］

的函數(shù)φ(t)就被稱為基小波，由基小波可以派生一個(gè)小波函數(shù)簇

式中a表示膨脹系數(shù)，b表示平移參數(shù)。在這個(gè)小波基上將函數(shù)f(t)分解即可得到該序列的連續(xù)小波變換

3 水平加速度改正的小波預(yù)處理

3.1 水平加速度改正的計(jì)算模型

計(jì)算水平加速度改正有兩種基本方法:兩步法和一步法。兩步法先計(jì)算穩(wěn)定平臺(tái)的傾角，再計(jì)算水平加速度改正，其數(shù)學(xué)模型為［4］:

其中α和β為平臺(tái)傾角，aE和aN分別是飛機(jī)的東向、北向水平加速度，可由機(jī)載GPS導(dǎo)出，g為測(cè)得的重力值。

一步法假設(shè)穩(wěn)定平臺(tái)的水平軸嚴(yán)格垂直，令fx和fy分別為平臺(tái)橫軸和縱軸方向的水平加速度，則有［8］

文獻(xiàn)［9］推導(dǎo)出這兩種方法在理論上等價(jià)。

由于式(6)中平方項(xiàng)的存在，產(chǎn)生的系統(tǒng)誤差具有非線性特性，這加大了分析加速度譜的難度。需同時(shí)分析航空重力測(cè)量厄特弗斯改正［2］

式中VE、VN分別是水平速度的東向分量和北向分量，M、N是所在位置子午圈和卯酉圈曲率半徑，h是飛機(jī)航高，ω是地球自轉(zhuǎn)角速度，φ是航點(diǎn)緯度。分析式(6)和式(7)都含有平方項(xiàng)，但是數(shù)據(jù)處理中不同尺度的預(yù)濾波處理對(duì)水平加速度影響較大，而厄特弗斯改正基本不受預(yù)濾波的影響。進(jìn)一步研究水平加速度改正的計(jì)算模型發(fā)現(xiàn):式(6)中所涉及的數(shù)據(jù)不僅包含原始觀測(cè)量fx和fy，也含有計(jì)算得到的載體加速度aE和aN，因此可以初步確定計(jì)算載體加速度時(shí)會(huì)引起數(shù)據(jù)頻譜中噪聲成分的增加。圖1的實(shí)測(cè)數(shù)據(jù)分析也證實(shí)了這個(gè)結(jié)論:速度的能量基本集中在頻譜圖的低頻部分，加速度的能量中含有大量頻譜集中在0.01 Hz附近的噪聲。

圖1 載體速度分量和水平加速度分量的頻譜圖Fig.1 Spectrum diagram of the carrier’s velocity constituent and horizontal acceleration constituent

3.2 小波變換尺度的確定

預(yù)濾波的最佳尺度存在各種假定，較常見的一種假定是認(rèn)為預(yù)濾波的尺度與最終濾波器的尺度相近時(shí)最佳，由經(jīng)驗(yàn)可知，航空重力數(shù)據(jù)濾波處理時(shí)采用的濾波器時(shí)域周期為100 s或者更高［9］，但通過實(shí)驗(yàn)發(fā)現(xiàn)，采用時(shí)域周期為50 s的低通濾波器進(jìn)行預(yù)濾波時(shí)輸出的加速度分量已經(jīng)全部接近零，因此采用過長(zhǎng)的低通濾波器進(jìn)行濾波將濾掉加速度信息的有效部分，可見利用最終的濾波器進(jìn)行預(yù)濾波處理是不可行的。

考慮到航空重力測(cè)量系統(tǒng)的構(gòu)成，水平加速度計(jì)和測(cè)姿裝置固定在穩(wěn)定平臺(tái)內(nèi)部，穩(wěn)定平臺(tái)通過專門的減震系統(tǒng)懸浮在框架之間［10］。保護(hù)穩(wěn)定平臺(tái)的減震系統(tǒng)由減震繩和8個(gè)油壓減震器構(gòu)成，系統(tǒng)外部的高頻噪聲能否傳遞到穩(wěn)定平臺(tái)內(nèi)部在很大程度上取決于減震系統(tǒng)的性能。出于技術(shù)保密的需要，Lamp;R公司并未公布其產(chǎn)品的減震系統(tǒng)性能，我們只能通過相關(guān)數(shù)據(jù)的頻譜特性來(lái)估計(jì)減震系統(tǒng)的減震效果。穩(wěn)定平臺(tái)內(nèi)部固定有水平加速度計(jì)和姿態(tài)測(cè)定裝置，選取記錄的內(nèi)部水平加速度數(shù)據(jù)來(lái)分析減震系統(tǒng)的性能。令fx表示平臺(tái)橫軸加速度計(jì)敏感到的橫向水平加速度，fy表示平臺(tái)縱軸加速度計(jì)敏感到的縱向水平加速度，對(duì)某測(cè)區(qū)所有測(cè)線的fx和fy分量進(jìn)行頻譜分析，結(jié)果如圖2所示，穩(wěn)定平臺(tái)內(nèi)部加速度計(jì)測(cè)得的水平加速度頻譜基本都在0～0.1 Hz之間，據(jù)此建議該航空重力測(cè)量系統(tǒng)應(yīng)選定10 s的預(yù)濾波尺度。

圖2 試驗(yàn)區(qū)域所有測(cè)線橫向加速度分量的頻譜圖Fig.2 Spectrum diagram of horizontal acceleration component of all flight path in the test area

穩(wěn)定平臺(tái)的兩根水平軸互相垂直，加速度計(jì)敏感到的水平加速度f(wàn)x和fy與GPS測(cè)量所得的載體水平加速度aE和aN存在如下關(guān)系［8］

其中α、β表示穩(wěn)定平臺(tái)的傾角。從式(8)可以看出:平臺(tái)敏感到的水平加速度f(wàn)x和fy是載體水平加速度aE和aN的線性組合，根據(jù)傅里葉變換的性質(zhì)，fx和fy的頻譜也是載體水平加速度aE和aN頻譜的線性組合，即:在理論上fx、fy、aE和aN具有相同的頻譜構(gòu)成，由此可以確定aE、aN、fx和fy的預(yù)濾波尺度相同。

根據(jù)小波分析的思想，如果采集得到的是頻率范圍為［0，Ω］的信號(hào)，那么對(duì)該信號(hào)進(jìn)行一次小波分解可以將原信號(hào)分成低頻部分和高頻部分，它們的頻率范圍分別為［0，Ω/2］、［Ω/2，Ω］。再次對(duì)低頻部分進(jìn)行分解，可以得到頻帶在［0，Ω/4］、［Ω/4，Ω/2］中的兩個(gè)信號(hào)，以此類推可以將信號(hào)不停地分解至最頂層。當(dāng)然，也可以采用小波包的方法對(duì)信號(hào)的高頻部分進(jìn)行同樣的分解，這樣就可以完成濾波和去噪的工作［11］。航空重力測(cè)量觀測(cè)數(shù)據(jù)的采樣間隔是1 s，對(duì)應(yīng)的小波分析頻率為0.5 Hz，進(jìn)行三層的小波變換基本符合水平加速度改正的頻率要求。

3.3 水平加速度改正的小波預(yù)處理

小波濾波方法對(duì)噪聲的性質(zhì)有3個(gè)基本假設(shè)［12］:一是噪聲均勻低分布在所有系數(shù)中;二是數(shù)據(jù)中的噪聲部分經(jīng)過小波變換后大多數(shù)為零或者近似為零;三是噪聲水平不能太高。航空重力測(cè)量中的噪聲可以看做隨機(jī)的Gauss白噪聲，符合第一個(gè)條件。載體的加速度分量aE小波分解提取后如圖3所示。

從圖3可以看出，原始數(shù)據(jù)經(jīng)小波變換后的細(xì)節(jié)部分基本為0。結(jié)合圖1和圖2，數(shù)據(jù)的噪聲水平不高，符合小波濾波方法對(duì)噪聲的3個(gè)基本假設(shè)。因此可以利用小波濾波方法對(duì)加速度數(shù)據(jù)進(jìn)行濾波處理。

選取Harr小波、db4小波、db6小波和sym4小波對(duì)水平加速度改正進(jìn)行預(yù)濾波處理，它們的函數(shù)如圖4所示。

為驗(yàn)證水平加速度改正的小波濾波處理效果，對(duì)某地區(qū)航空重力測(cè)量數(shù)據(jù)進(jìn)行處理分析。該測(cè)區(qū)屬陸海交界區(qū)域，載體平均飛行速度為60 ms-1，航高約500 m，運(yùn)載平臺(tái)為國(guó)產(chǎn)航測(cè)機(jī)?？紤]到目前航空重力測(cè)量數(shù)據(jù)濾波處理普遍采用FIR低通濾波器，故此處計(jì)算不同預(yù)濾波尺度下該測(cè)區(qū)某條測(cè)線水平加速度改正后將其結(jié)果與10 s預(yù)濾波結(jié)果進(jìn)行對(duì)比分析。圖5展示了不同尺度預(yù)濾波處理對(duì)水平加速度改正最終濾波結(jié)果的影響，從圖5可以看出，水平加速度改正經(jīng)小波預(yù)濾波處理后的最終結(jié)果與FIR預(yù)處理的精度相當(dāng)。表1列出不同方法處理了某條測(cè)線對(duì)最終結(jié)果影響的數(shù)據(jù)統(tǒng)計(jì)特性。

圖3 水平加速度分量aE的db4小波分解Fig.3 db4 wavelet decomposition of the horizontal acceleration component aE

圖4 本文選取的小波函數(shù)Fig.4 Wavelet functions chosen in this article

圖5 不同小波預(yù)濾波對(duì)水平加速度改正最終結(jié)果的影響Fig.5 The influence of the different wavelets pre-filters to the final results of the horizontal acceleration correction

表1 與10 s FIR預(yù)濾波相比，小波預(yù)濾波對(duì)最終結(jié)果影響的數(shù)據(jù)統(tǒng)計(jì)特性(單位：10-5ms-2)Tab.1 Statistic characteristic of the wavelet pre-filter impact on the final results compared with 10 s FIR pre-filter(unit：10-5ms-2)

3.4 精度評(píng)定

從表1可以看出，水平加速度改正經(jīng)小波濾波和FIR低通濾波預(yù)處理后的差值小于 1×10-5ms-2。由于航空重力測(cè)量中加速度改正的真值無(wú)法確定，因此無(wú)法從表1判斷出各濾波方法的優(yōu)劣。選取四種方法中的任一個(gè)預(yù)濾波結(jié)果作為真值，再加入均值為零、方差為0.02的Gauss白噪聲，分別用FIR低通濾波和不同的小波對(duì)其進(jìn)行濾波處理，并與不加噪聲的數(shù)據(jù)進(jìn)行對(duì)比，其終結(jié)果如表2。

表2 不同預(yù)濾波方法對(duì)最終結(jié)果的影響(單位：10-5ms-2)Tab.2 Influence of different wavelet pre-filter on the final results(unit：10-5ms-2)

從表2可以看出，與FIR低通預(yù)濾波相比，小波預(yù)濾波的精度有不同程度的提高，特別是sym8小波和db小波簇更是將水平加速度改正的誤差控制在0.2×10-5ms-2內(nèi)。遺憾的是航空重力測(cè)量數(shù)據(jù)的信噪比低，目前尚未找到行之有效的小波濾波方法，因此將小波濾波用于最終的濾波處理還需進(jìn)一步的研究。

4 結(jié)論

傳統(tǒng)的FIR低通濾波首先用Fourier變換將信號(hào)變換到頻域，然后再進(jìn)行低通濾波處理。當(dāng)信號(hào)的有效頻譜中混有噪聲時(shí)，F(xiàn)IR的濾波效果較差，因此FIR低通濾波器不能抑制混入信號(hào)頻帶內(nèi)部的噪聲分量。采用Haar小波、db小波簇和sym小波簇對(duì)水平加速度改正進(jìn)行預(yù)濾波處理，可以得出:1)小波變換能夠抑制混入有效頻帶范圍內(nèi)的噪聲，進(jìn)一步提高了數(shù)據(jù)處理的精度;2)小波變換實(shí)現(xiàn)簡(jiǎn)單，對(duì)邊界數(shù)據(jù)的影響較小。表2的精度分析表明，在不產(chǎn)生邊界問題的前提下，sym8小波和 Daubechies小波簇精度最高，適用于水平加速度改正的預(yù)濾波處理。

1 孫中苗.航空重力測(cè)量理論、方法及應(yīng)用［D］.鄭州:信息工程大學(xué)測(cè)繪學(xué)院，2004.(Sun Zhongmiao.Theory，methods and application of airborne gravimetry［D］.Zhengzhou: Institute of Surveying and Mapping，2004)

2 石磐，孫中苗，肖云.航空重力測(cè)量中水平加速度改正的計(jì)算與頻域分析［J］.武漢大學(xué)學(xué)報(bào)(信息科學(xué)版)，2001，26(6):549-553.(Shi Pan，Sun Zhongmiao and Xiao Yun.Calculation and spectra analysis of horizontal acceleration correction(HACC)for airborne gravity［J］.Geomatics and Information Science of Wuhan University，2001，26(6): 549-553)

3 孫中苗，等.Lamp;R航空重力儀的水平加速度改正［J］.測(cè)繪科學(xué)技術(shù)學(xué)報(bào)，2007，24(4):259-262.(Sun Zhongmiao，et al.Horizontal acceleration correction for the airborne gravity［J］.Journal of Geomatics Science and Technology，2007，24(4):259-262)

4 Donald B P and Andrew T W.Wavelet methods for time series analysis［M］.Cambridge:Cambridge University Publishing Company，2006.

5 高成.Matlab小波分析與應(yīng)用(第2版)［M］.北京:國(guó)防工業(yè)出版社，2007.(Gao Cheng.Wavelet analysis and its application(2nd)［M］.Beijing:National Defence Industry Press，2007)

6 耿則勛.小波變換理論及在遙感影像壓縮中的應(yīng)用［M］.北京:測(cè)繪出版社，2002.(Gen Zexun.Wavelet transform and its application in image compression of remote sensing［M］.Beijing:Surveying and Mapping Publication，2002)

7 孫延奎.小波分析及其應(yīng)用［M］.北京:機(jī)械工業(yè)出版社，2005.(Sun Yankui.Wavelet analysis and its application［M］.Beijing:Machine Industry Press，2005)

8 Peters M F and Brozena J M.Methods to improve existing shipboard gravimeters for airborne gravimetry［A］.Presented at IAG Symposium on Airborne Gravity Field Determination［C］.Boulder，Colorado，1995，39-45.

9 孫中苗，等.航空重力測(cè)量數(shù)據(jù)的濾波與處理［J］.地球物理學(xué)進(jìn)展，2004，19(1):119-124.(Sun Zhongmiao，et al.Filtering and processing for the airborne gravity data［J］.Progress in Geophysics，2004，19(1):119-124)

10 LaCosteamp;Romberg Company.Model“S”air-sea dynamic gravity meter system II instruction manual［S］.Texas: LRC，2003.

11 程正興.小波分析與應(yīng)用實(shí)例［M］.西安:西安交通大學(xué)出版社，2006.(Cheng Zhengxing.Examples of wavelet analysis and its application［M］.Xi’an:Xi’an Traffic University Press，2006)

12 潘泉，等.小波濾波方法及應(yīng)用［M］.北京:清華大學(xué)出版社，2005.(Pan Quan，et al.Wavelet filtering method and its application［M］.Beijing:Tsinghua University Press，2005)

WAVELET PRE-PROCESSING FOR HORIZONTAL ACCELERATION OF AERIAL GRAVITY MEASUREMENT

Tian Yanfeng1)，Li Shanshan2)and Xiao Fan3)

(1)PLA 73603 Troops，Nanjing 210049 2)Institute of Surveying and Mapping，Information Engineering University，Zhengzhou 450052 3)PLA 61365 Troops，Tianjin 300140)

According to the spectrum characteristic of the horizontal acceleration correction in the aerial gravity measurement and the partition rules of wavelet transformation frequency band，the wavelet pre-filter scale of the horizontal acceleration correction is deduced.On the basis of the correlation among the wavelet transform scales of the signal，as an example，some surveying data were pre-processed by Haar wavelet、db4 wavelet、db6 wavelet and sym6 wavelet，and the results are compared with the results from FIR low pass filter in common use.The comparison shows that the precision of the wavelet pre-filter is higher than the FIR low pass filter in the pre-processing of the horizontal acceleration correction and it is more symmetrical as pre-processed by db serial wavelets.Thus it is proposed that the Daubechies wavelets or the sym8 wavelet with much longer branch are more useful in the pre-processing of the horizontal acceleration correction.

aerial gravity measurement;pre-filter;horizontal acceleration correction;E?tv?s correction;Fourier transform;wavelet transform

(1)解放軍73603部隊(duì)，南京 210049 2)解放軍信息工程大學(xué)測(cè)繪學(xué)院，鄭州 450052 3)解放軍61365部隊(duì)，天津300140)

1671-5942(2012)02-0115-05

2011-12-04

田顏鋒，男，1981年生，碩士，主要研究方向?yàn)槲锢泶蟮販y(cè)量.E-mail:yanftian@126.com

A