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

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

0 引言

地球重力是地球引力和地球自轉(zhuǎn)引起的慣性力的合力，是地球的基本物理特征之一。精確的局部地球重力場信息對自然資源勘察、地質(zhì)科學(xué)研究、自然災(zāi)害監(jiān)測及預(yù)警、地球極區(qū)勘探、高精度戰(zhàn)略武器系統(tǒng)研制保障具有重要的意義[1-3]。

航空重力測量是一種以航空器為載體，綜合運(yùn)用重力敏感器、高精度姿態(tài)穩(wěn)定系統(tǒng)、定位系統(tǒng)，獲取近地面或近海重力場的方法[4]。

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

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

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

當(dāng)前，國內(nèi)對采用加速度計式重力敏感器的三軸穩(wěn)定平臺重力儀數(shù)據(jù)后處理方法研究較少。東南大學(xué)蔡體菁帶領(lǐng)的研究團(tuán)隊采用Kalman濾波的方法對平臺重力儀實(shí)測數(shù)據(jù)進(jìn)行了處理，得到了較好的結(jié)果，此外該團(tuán)隊還提出了巴特沃斯與Kalman平滑濾波級聯(lián)的處理方法。

本文采用頻域信息處理的方法，結(jié)合平臺重力儀的特點(diǎn)，系統(tǒng)地給出了一套在頻域內(nèi)有嚴(yán)格規(guī)范和明確含義的航空重力數(shù)據(jù)后處理方法。

1 航空重力測量的基本原理與數(shù)學(xué)模型

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

1.1 平臺式航空重力儀系統(tǒng)組成

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

其中重力儀主機(jī)的核心裝置為慣性穩(wěn)定平臺，主要由1個高精度石英撓性加速度計式重力敏感器、2個動力調(diào)諧陀螺、2個導(dǎo)航級石英撓性加速度計和3個環(huán)架組成。該慣性穩(wěn)定平臺在水平面內(nèi)穩(wěn)定精度優(yōu)于10″，可以隔離載體角運(yùn)動對重力敏感器的影響，同時可以跟蹤當(dāng)?shù)氐乩碜鴺?biāo)系，嚴(yán)格保持重力敏感器指向垂向[14]。

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

1.2 航空重力測量的數(shù)學(xué)原理

航空重力測量的數(shù)學(xué)模型是由慣性導(dǎo)航原理中比力方程推導(dǎo)而來的。本文研究的是標(biāo)量重力儀數(shù)據(jù)后處理，以下給出標(biāo)量形式重力測量的基本方程

(1)

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

(2)

由于重力測量一般都是為了獲取被測區(qū)域的重力異常值，因此從重力測量值中減去標(biāo)準(zhǔn)重力數(shù)值，得到以下計算重力異常的公式

(3)

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

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

1)當(dāng)飛行海拔高度h≤700m

δgh=0.3086h

(4)

2)當(dāng)飛行海拔高度h>700m

δgh=0.3086h-0.073h2

(5)

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

2 航空重力數(shù)據(jù)后處理流程與關(guān)鍵技術(shù)

根據(jù)上述重力測量的基本原理，結(jié)合平臺式重力儀的結(jié)構(gòu)特點(diǎn)，給出如圖1所示的航空重力數(shù)據(jù)后處理流程。

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

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

1)重力敏感器和差分GPS的輸出頻率分別為100Hz和2Hz，且2個通道獨(dú)立工作，因此需要實(shí)現(xiàn)通道之間的采樣率統(tǒng)一；

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

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

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

針對上述問題，提出了以下幾點(diǎn)數(shù)據(jù)處理的關(guān)鍵技術(shù)。

2.1 采用降采樣率方法實(shí)現(xiàn)通道間采樣率統(tǒng)一

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

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

Fs/M≥2Fc

(6)

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

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

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

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

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

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

2.2 采用中心差分器求解載體運(yùn)動加速度

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

理想差分器的頻率響應(yīng)為

Hd(ejω)=jω

(7)

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

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

(8)

中心差分運(yùn)算可以通過FIR濾波器形式實(shí)現(xiàn)，其幅頻特性曲線如圖5所示，相頻特性曲線如圖6所示。

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

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

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

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

TGroup=-dθ/dω

(9)

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

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

在數(shù)據(jù)處理的過程中，采用FIR低通濾波器會產(chǎn)生相位延遲，導(dǎo)致輸出信號相對于輸入信號在時間上發(fā)生平移。體現(xiàn)在重力數(shù)據(jù)的處理結(jié)果上，會使得濾波解算得到的重力異常位置偏離真實(shí)地理位置。因此在設(shè)計低通濾波器的過程中，需要采用零相位濾波進(jìn)行處理[21]。

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

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

(10)

將得到的結(jié)果u(n)進(jìn)行時間反轉(zhuǎn)為v(n)

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

(11)

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

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

(12)

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

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

(13)

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

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

(14)

然后進(jìn)行時間反轉(zhuǎn)

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

(15)

再次通過濾波器

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

(16)

最后再次將所得結(jié)果進(jìn)行時間反轉(zhuǎn)，有

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

(17)

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

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

3 航空重力實(shí)測數(shù)據(jù)處理結(jié)果與精度評估

3.1 試驗概況

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

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

3.2 航空重力測量結(jié)果精度評估方法

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

3.3 濾波周期140s的重復(fù)測線結(jié)果

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

4個架次的重復(fù)測線處理結(jié)果及其對應(yīng)的GT-2A測量結(jié)果分別如圖8～圖11所示。

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

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

上述處理曲線的內(nèi)符合和外符合精度結(jié)果統(tǒng)計如表1所示。

將3個架次東西方向12條重復(fù)測線的結(jié)果綜合到一起進(jìn)行分析，并與其對應(yīng)的GT-2A測量結(jié)果進(jìn)行對比，其結(jié)果如圖12所示。

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

3.4 濾波周期100s的重復(fù)測線結(jié)果

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

4個架次的重復(fù)測線處理結(jié)果及其對應(yīng)的GT-2A測量結(jié)果分別如圖13～圖16所示。

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

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

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

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

將3個架次東西方向12條重復(fù)測線的結(jié)果綜合到一起進(jìn)行分析，并與其對應(yīng)的GT-2A測量結(jié)果進(jìn)行對比，其結(jié)果如圖17所示。

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

上述處理曲線的內(nèi)符合和外符合精度結(jié)果統(tǒng)計如表2所示。

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

3.5 起伏飛行條件下的測試結(jié)果

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

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

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

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

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

在140s濾波周期下，起伏測量結(jié)果與GT-2A的標(biāo)準(zhǔn)差為0.81mGal，與同濾波尺度下重復(fù)測線均值的標(biāo)準(zhǔn)差為0.60mGal。

在100s濾波周期下，起伏測量結(jié)果與GT-2A的標(biāo)準(zhǔn)差為1.11mGal，與同濾波尺度下重復(fù)測線均值的標(biāo)準(zhǔn)差為1.03mGal。

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

4 結(jié)論

本文針對平臺式重力儀航空測量數(shù)據(jù)后處理的要求，結(jié)合平臺式重力儀的系統(tǒng)原理和結(jié)構(gòu)特點(diǎn)，設(shè)計了一套數(shù)據(jù)后處理的流程，并對其中的關(guān)鍵處理環(huán)節(jié)進(jìn)行了頻域分析。通過對飛行試驗實(shí)測數(shù)據(jù)進(jìn)行處理，得到如下結(jié)論：

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

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

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

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