航空重力測量水平加速度改正的小波預處理*

田顏鋒 李姍姍 肖 凡

航空重力測量水平加速度改正的小波預處理*

田顏鋒1)李姍姍2)肖 凡3)

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

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

1 引言

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

2 小波變換的基本原理

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

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

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

3 水平加速度改正的小波預處理

3.1 水平加速度改正的計算模型

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

其中α和β為平臺傾角，aE和aN分別是飛機的東向、北向水平加速度，可由機載GPS導出，g為測得的重力值。

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

文獻［9］推導出這兩種方法在理論上等價。

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

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

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

3.2 小波變換尺度的確定

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

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

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

穩(wěn)定平臺的兩根水平軸互相垂直，加速度計敏感到的水平加速度fx和fy與GPS測量所得的載體水平加速度aE和aN存在如下關系［8］

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

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

3.3 水平加速度改正的小波預處理

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

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

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

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

圖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 不同小波預濾波對水平加速度改正最終結果的影響Fig.5 The influence of the different wavelets pre-filters to the final results of the horizontal acceleration correction

表1 與10 s FIR預濾波相比，小波預濾波對最終結果影響的數(shù)據(jù)統(tǒng)計特性(單位：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 精度評定

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

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

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

4 結論

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

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

2 石磐，孫中苗，肖云.航空重力測量中水平加速度改正的計算與頻域分析［J］.武漢大學學報(信息科學版)，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］.測繪科學技術學報，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小波分析與應用(第2版)［M］.北京:國防工業(yè)出版社，2007.(Gao Cheng.Wavelet analysis and its application(2nd)［M］.Beijing:National Defence Industry Press，2007)

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

7 孫延奎.小波分析及其應用［M］.北京:機械工業(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 孫中苗，等.航空重力測量數(shù)據(jù)的濾波與處理［J］.地球物理學進展，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 程正興.小波分析與應用實例［M］.西安:西安交通大學出版社，2006.(Cheng Zhengxing.Examples of wavelet analysis and its application［M］.Xi’an:Xi’an Traffic University Press，2006)

12 潘泉，等.小波濾波方法及應用［M］.北京:清華大學出版社，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部隊，南京 210049 2)解放軍信息工程大學測繪學院，鄭州 450052 3)解放軍61365部隊，天津300140)

1671-5942(2012)02-0115-05

2011-12-04

田顏鋒，男，1981年生，碩士，主要研究方向為物理大地測量.E-mail:yanftian@126.com

A