平臺式重力儀航空重力測量數據后處理技術研究

閆 方，胡平華，趙 明，姜作喜，羅 鋒，唐江河，劉東斌

(1.北京自動化控制設備研究所，北京 100074；2.中國國土資源航空物探遙感中心，北京 100083)

0 引言

地球重力是地球引力和地球自轉引起的慣性力的合力，是地球的基本物理特征之一。精確的局部地球重力場信息對自然資源勘察、地質科學研究、自然災害監測及預警、地球極區勘探、高精度戰略武器系統研制保障具有重要的意義[1-3]。

航空重力測量是一種以航空器為載體，綜合運用重力敏感器、高精度姿態穩定系統、定位系統，獲取近地面或近海重力場的方法[4]。

結合重力測量的發展歷程，當前還有走行式地面重力測量、海洋船載重力測量和衛星測高重力測量等方法[5]。其中走行式方法借助靜態重力儀在地面進行，這種測量方式耗時長、效率低、人工成本高，同時受山川、湖泊、沼澤等自然條件的制約；船載方式更適合遠洋應用場景，在近海岸區域由于水流波動較大，容易造成測量精度下降；借助衛星測高資料反演得到的重力場只能反映地球重力場的中低頻信息，限制了衛星重力資料的應用領域[6]。而采用航空重力儀系統的機載測量方式，可以克服上述測量手段中的缺點，高效快速地獲取較大區域內的重力資料，同時能夠保證測量的精度與分辨率[7-8]。

對于航空重力數據處理，目前主要有頻域有限脈沖響應(Finite Impulse Response，FIR)低通濾波、時域Kalman濾波以及時頻域小波分解去噪方法[9]。其中Kalman濾波方法對于重力異常模型的精度要求較高，小波方法在重力測量領域缺乏明確的理論依據，因此實際應用都受到了限制。而頻域FIR低通濾波方法物理含義清晰，實際應用廣泛。

除了上述常規手段，國外還有具有特色的濾波方法。Von Frese等提出了波數相關濾波器(Wavenumber Correlation Filter，WCF)，該方法通過比較測量數據與協同數據頻譜的相關程度來實現噪聲分離[10]。Alberts提出了頻域加權的方法，對測量數據的不同頻段賦予不同權值，以此進行降噪處理[11]。Bolotin在時域處理方法的基礎上，將地球重力場的非均勻一致性建模為多狀態隱性馬爾可夫模型，對結果的細節有所提升[12]。

當前，國內對采用加速度計式重力敏感器的三軸穩定平臺重力儀數據后處理方法研究較少。東南大學蔡體菁帶領的研究團隊采用Kalman濾波的方法對平臺重力儀實測數據進行了處理，得到了較好的結果，此外該團隊還提出了巴特沃斯與Kalman平滑濾波級聯的處理方法。

本文采用頻域信息處理的方法，結合平臺重力儀的特點，系統地給出了一套在頻域內有嚴格規范和明確含義的航空重力數據后處理方法。

1 航空重力測量的基本原理與數學模型

航空重力測量圍繞2個基本問題展開，一個是如何在運動的環境下穩定重力敏感器的指向，使其嚴格保持垂向；另一個是如何從重力敏感器輸出的慣性加速度中分離得到重力加速度[13]。以下將結合平臺式航空重力儀的系統組成和重力測量的基本數學模型進行簡述。

1.1 平臺式航空重力儀系統組成

平臺式重力儀系統硬件部分主要由重力儀主機、差分全球定位系統(Differential Global Positioning System，DGPS)、不間斷電源、顯控存儲裝置和減振器組成。

其中重力儀主機的核心裝置為慣性穩定平臺，主要由1個高精度石英撓性加速度計式重力敏感器、2個動力調諧陀螺、2個導航級石英撓性加速度計和3個環架組成。該慣性穩定平臺在水平面內穩定精度優于10″，可以隔離載體角運動對重力敏感器的影響，同時可以跟蹤當地地理坐標系，嚴格保持重力敏感器指向垂向[14]。

差分GPS是航空重力測量系統中另一關鍵部分，由安裝在載體上的移動站和固定在地面的基站組成。原始差分衛星測量數據經過解算可以得到載體高精度的位置與速度信息，用來計算載體垂向方向的運動加速度以及重力相關的改正項[15]。

1.2 航空重力測量的數學原理

航空重力測量的數學模型是由慣性導航原理中比力方程推導而來的。本文研究的是標量重力儀數據后處理，以下給出標量形式重力測量的基本方程

(1)

將式(1)中后三項單獨記為δaE，即為厄特弗斯改正項

(2)

由于重力測量一般都是為了獲取被測區域的重力異常值，因此從重力測量值中減去標準重力數值，得到以下計算重力異常的公式

(3)

其中，g0為正常重力場，根據不同的需求采用相應的正常重力模型進行計算。常用的正常重力公式主要有赫爾默特公式和1985年國際正常重力公式。需要注意的是，這些公式計算值為測點在地球旋轉橢球體表面處的正常重力值，而高度每增加1m，重力值就降低約0.3086mGal，因此在航空重力數據處理中使用的正常重力場公式需包含高度校正部分。

以下給出高度校正部分的公式[16]：

1)當飛行海拔高度h≤700m

δgh=0.3086h

(4)

2)當飛行海拔高度h>700m

δgh=0.3086h-0.073h2

(5)

經過上述流程的計算，即可實現重力敏感器輸出中慣性加速度和重力加速度的分離，同時還對重力異常信號進行了必要的改正。在實際航空重力測量中，重力敏感器的輸出還包含由飛機發動機以及空氣湍流造成的隨機振動加速度，這部分噪聲的強度是重力異常的數萬倍。因此，需要利用濾波的方法進一步去除重力異常中的高頻噪聲[17]。

2 航空重力數據后處理流程與關鍵技術

根據上述重力測量的基本原理，結合平臺式重力儀的結構特點，給出如圖1所示的航空重力數據后處理流程。

圖1 航空重力數據處理流程圖Fig.1 Flow chart of gravity measurement data processing

結合上述流程圖，對于平臺式航空重力儀的數據處理，有以下幾個問題需要注意：

1)重力敏感器和差分GPS的輸出頻率分別為100Hz和2Hz，且2個通道獨立工作，因此需要實現通道之間的采樣率統一；

2)在對各個通道內原始測量數據進行處理的過程中，還需考慮每一步操作對信號相位的影響，從而避免在作差環節因相位不同步而引入誤差；

3)在對計算得到的重力異常原始值進行濾波時，同樣需要考慮信號的相位變化，保證解算得到的重力異常值可以歸算到正確的地理位置。

采用低通FIR濾波方法的前提是認為重力異常信號具有低頻特性，因此數據處理的每個環節都需要有明確的頻域含義，避免計算過程對原始測量信號的低頻信息產生破壞，進而保證數據處理的精度。

針對上述問題，提出了以下幾點數據處理的關鍵技術。

2.1 采用降采樣率方法實現通道間采樣率統一

多采樣率數字信號處理方法有著嚴格的頻域理論基礎，在滿足一定條件的情況下可以提升或者降低數據的采樣率。以下簡述該方法的基本思想。

首先介紹降低采樣率的方法。現有原始信號序列x(n)，信號中最高頻率為Fc，原始采樣頻率為Fs。降低采樣率M倍最直接的方法是在序列中每M個點抽取一個點，組成一個新的序列x′(n)。根據奈奎斯特采樣定律，新序列不發生頻譜混疊的必要條件是新的采樣率大于信號內最高頻率的2倍，即

Fs/M≥2Fc

(6)

由于M是可變的，為了保證式(6)始終成立，可以首先對原始信號進行低通濾波以實現抗混疊，然后再進行抽取。其流程圖如圖2所示。

圖2 降采樣流程圖Fig.2 Flow chart of downsampling

對于提升L倍采樣率，首先在每2個相鄰點之間插入(L-1)個零值，然后通過一個低通濾波器，將插入的零值點計算得到相應的抽樣值。該方法在理論上的依據可以參考相關文獻[18]。其實現的流程圖如圖3所示。

圖3 升采樣流程圖Fig.3 Flow chart of upsampling

綜合使用上述方法，可以將通道間采樣統一到相同頻率。本文采用的方案是將重力儀的采樣率直接降低到2Hz，為此設計了截止頻率為1Hz的FIR低通濾波器，其幅頻特性如圖4所示。

圖4 FIR低通濾波器幅頻響應圖(1Hz)Fig.4 Amplitude frequency response curve ofFIR low pass filter (1Hz)

2.2 采用中心差分器求解載體運動加速度

差分GPS可以輸出高精度定位信息，對其中的高度序列做二次差分，可以計算得到載體運動加速度。

理想差分器的頻率響應為

Hd(ejω)=jω

(7)

由式(7)可知，理想差分器的幅頻特性為線性增長。

工程應用中，通常采用在低頻范圍內線性度較好的牛頓中心差分器[19-20]。將差分GPS輸出的高度序列記為x(n)，則牛頓中心差分計算過程為

(8)

中心差分運算可以通過FIR濾波器形式實現，其幅頻特性曲線如圖5所示，相頻特性曲線如圖6所示。

圖5 差分器幅頻特性圖Fig.5 Amplitude frequency response curve of differential filter

圖6 差分器相頻特性圖Fig.6 Phase frequency response curve of differential filter

從幅頻特性圖中可以看出，中心差分器在低頻段內線性度較好。

由相頻特性圖可以看出，其相位具有線性遞減特性。在數字信號處理中，群延遲定義為相位變化與頻率變化的比值，即：

TGroup=-dθ/dω

(9)

群延遲反映了信號中所有頻率成分的時間延遲，對于線性相位的系統，群延遲是常值。借助這一方法，可以精確計算得到中心差分器造成的時延，并在計算過程中予以補償。

2.3 采用零相位濾波器保證重力異常的相位精度

在數據處理的過程中，采用FIR低通濾波器會產生相位延遲，導致輸出信號相對于輸入信號在時間上發生平移。體現在重力數據的處理結果上，會使得濾波解算得到的重力異常位置偏離真實地理位置。因此在設計低通濾波器的過程中，需要采用零相位濾波進行處理[21]。

本文采用Forward-Backward零相位濾波器，其實現方法是：首先將輸入信號x(n)按正向順序通過濾波器，即與數字濾波器沖激響應序列h(n)進行卷積運算

u(n)=x(n)*h(n)

(10)

將得到的結果u(n)進行時間反轉為v(n)

v(n)=u(N-1-n)

(11)

再次通過濾波器得到序列w(n)

w(n)=v(n)*h(n)

(12)

最后將這一序列再次進行時間反轉，即得到精確的零相位結果y(n)

y(n)=w(N-1-n)

(13)

為了進一步說明該方法的零相位特性，以下給出這一過程的頻域描述，其中濾波器的頻域表示為H(ejw)，輸入信號正向通過濾波器后有

U(ejω)=X(ejω)H(ejω)

(14)

然后進行時間反轉

V(ejω)=e-jω(N-1)U(e-jω)

(15)

再次通過濾波器

W(ejω)=V(ejω)H(ejω)

(16)

最后再次將所得結果進行時間反轉，有

Y(ejω)=e-jω(N-1)W(e-jω)=X(ejω)|H(ejω)|2

(17)

由式(17)可以看出，輸入與輸出之間沒有相位變化，可以實現精確的零相位操作。圖7所示為零相位濾波流程圖。

圖7 零相位濾波流程圖Fig.7 Flow chart of zero-phase filter

3 航空重力實測數據處理結果與精度評估

3.1 試驗概況

為了驗證平臺式重力儀的實際工作性能，中國國土資源航空物探遙感中心在海南島附近某海域組織完成了機載航空重力測量試驗。本次試驗采用Cessna208b固定翼飛機，該機型配有自動駕駛儀，重力測線海拔高為600m，飛行速度為60m/s。

飛行試驗共進行了4個架次，其中東西方向在同一測線安排有3個架次，每架次均取得4條有效重復測線；在南北方向有1個架次共4條重復測線。此外，試驗還在東西測線上進行了1次起伏飛行。為更好地評價平臺式重力儀測量效果，試驗提供了同測線GT-2A的事先測量結果。

3.2 航空重力測量結果精度評估方法

航空重力測量數據精度評價指標主要有內符合精度和外符合精度。其中內符合精度是計算同一臺儀器在重復測線測量值的均方根(Root-Mean-Square，RMS)。外符合精度是以其他手段獲取的重力異常值作為標準值來統計的標準差(Standard Deviation，STD)。采用文獻[22]中給出的計算方法，分別給出濾波周期為140s和100s的統計結果。

3.3 濾波周期140s的重復測線結果

采用濾波周期為140s，即截止頻率為0.00714Hz的FIR低通濾波器，可以得到4個架次試驗的測量結果。結果中重力異常的地理分辨率為4.2km。

4個架次的重復測線處理結果及其對應的GT-2A測量結果分別如圖8～圖11所示。

圖10 東西方向第3架次測線及GT-2A結果(140s)Fig.10 Third east-west repeated measurementlines and the result of GT-2A(140s)

圖11 南北方向架次測線及GT-2A結果(140s)Fig.11 South-north repeated measurementlines and the result of GT-2A(140s)

上述處理曲線的內符合和外符合精度結果統計如表1所示。

將3個架次東西方向12條重復測線的結果綜合到一起進行分析，并與其對應的GT-2A測量結果進行對比，其結果如圖12所示。

表1 四架次重復測線內、外符合精度(140s)

3.4 濾波周期100s的重復測線結果

采用濾波周期為100s，即截止頻率為0.01Hz的FIR低通濾波器，可以得到4個架次試驗的測量結果。結果中重力異常的地理分辨率為3km。

4個架次的重復測線處理結果及其對應的GT-2A測量結果分別如圖13～圖16所示。

圖13 東西方向第1架次測線及GT-2A結果(100s)Fig.13 First east-west repeated measurementlines and the result of GT-2A(100s)

圖14 東西方向第2架次測線及GT-2A結果(100s)Fig.14 Second east-west repeated measurementlines and the result of GT-2A(100s)

圖15 東西方向第3架次測線及GT-2A結果(100s)Fig.15 Third east-west repeated measurementlines and the result of GT-2A(100s)

圖16 南北方向架次測線及GT-2A結果(100s)Fig.16 South-north repeated measurementlines and the result of GT-2A(100s)

將3個架次東西方向12條重復測線的結果綜合到一起進行分析，并與其對應的GT-2A測量結果進行對比，其結果如圖17所示。

圖17 東西方向所有測線及GT-2A結果(100s)Fig.17 All east-west repeated measurementlines and the result of GT-2A(100s)

上述處理曲線的內符合和外符合精度結果統計如表2所示。

表2 四架次重復測線內、外符合情況(100s)

3.5 起伏飛行條件下的測試結果

本次試驗在東西方向第2架次飛行中進行了2次起伏飛行，最大起伏高度為100m，垂向速度為2m/s左右。起伏飛行的高度與垂向速度變化曲線如圖18所示。

圖18 起伏飛行高度和垂向速度曲線圖Fig.18 Figure of height and vertical velocity underrise-and-fall flight

分別給出140s和100s濾波周期下，起伏飛行測得的重力異常與GT-2A平穩飛行測量結果及同架次其他4條重復線的均值結果的對比圖，分別如圖19和圖20所示。

圖19 起伏飛行測量結果與重復線均值及GT-2A結果比較(140s)Fig.19 Comparison of result under rise-and-fall flight conditionwith mean of repeated lines and the result of GT-2A(140s)

圖20 起伏飛行測量結果與重復線均值及GT-2A結果比較(100s)Fig.20 Comparison of result under rise-and-fall flight conditionwith mean of repeated lines and the result of GT-2A(100s)

在140s濾波周期下，起伏測量結果與GT-2A的標準差為0.81mGal，與同濾波尺度下重復測線均值的標準差為0.60mGal。

在100s濾波周期下，起伏測量結果與GT-2A的標準差為1.11mGal，與同濾波尺度下重復測線均值的標準差為1.03mGal。

上述指標及具體曲線表明，平臺式重力儀具備一定的起伏飛行能力。

4 結論

本文針對平臺式重力儀航空測量數據后處理的要求，結合平臺式重力儀的系統原理和結構特點，設計了一套數據后處理的流程，并對其中的關鍵處理環節進行了頻域分析。通過對飛行試驗實測數據進行處理，得到如下結論：

1)試驗中4個架次的重復線內符合精度達到：12條東西測線0.43mGal/140s和0.84mGal/100s，4條南北測線0.39mGal/140s和0.79mGal/100s，表明平臺式重力儀性能穩定可靠；

2)以GT-2A單次測量結果(其本身也存在一定誤差)為標準值統計外符合精度，12條東西測線為0.72mGal/140s和0.98mGal/100s，4條南北測線為1.41mGal/140s和1.53mGal/100s，東西測線精度較好，南北測線精度較差。經分析，原因可能是GT-2A的此次南北測線測量結果存在較大誤差，后續將針對此問題進行重點驗證；

3)在起伏高度為100m，垂向速度為2m/s的條件下，重力異常的內符合及外符合精度沒有顯著變大，表明平臺式重力儀具備一定的起伏飛行能力；

4)通過對以上實測數據的處理，驗證了本文設計的航空重力數據后處理方法的正確性。