航空重力測量重復(fù)測線數(shù)據(jù)處理技術(shù)研究

李新宇，周召發(fā)，張志利，常振軍，郝詩文

火箭軍工程大學(xué) 兵器發(fā)射理論和技術(shù)國家重點(diǎn)學(xué)科實(shí)驗(yàn)室，陜西 西安 710025

以飛機(jī)為載體的航空重力測量技術(shù)已成為高效測定地球重力場信息的主要手段[1]。在動態(tài)重力測量中，測量結(jié)果會受到發(fā)動機(jī)振動、載體動態(tài)性、空氣擾動等因素的干擾，測量噪聲廣泛分布，嚴(yán)重影響了重力測量的精度[2]。開展航空重力測量數(shù)據(jù)處理方法研究，提高數(shù)據(jù)處理的精度，具有十分重要的意義。

為了抑制噪聲對實(shí)現(xiàn)重力信息的高精度提取的影響，眾多學(xué)者對航空重力數(shù)據(jù)后處理進(jìn)行了深入研究。頻域有限沖激響應(yīng)（FIR）低通濾波方法在實(shí)際中被廣泛應(yīng)用，可以有效消除高頻噪聲[3]。國外提出了具有特色的濾波方法，如R.R.B.Von Frese 等[4]提出了波數(shù)相關(guān)濾波器（WCF），通過比較測量數(shù)據(jù)與協(xié)同數(shù)據(jù)頻譜的相關(guān)程度來實(shí)現(xiàn)噪聲分離；B.A.Alberts 等[5]提出了頻域加權(quán)的方法，對測量數(shù)據(jù)的不同頻段賦予不同的權(quán)值，以此來進(jìn)行降噪處理；Y.V.Bolotin 等[6]在時域處理方法的基礎(chǔ)上，將地球重力場的非均勻一致性建模為多狀態(tài)隱性馬爾可夫模型，對結(jié)果的細(xì)節(jié)有所提升。國內(nèi)方面，羅鋒等[7]提出使用Kalman平滑算法開展航空重力測量數(shù)據(jù)處理方法研究，在國內(nèi)首次解算出穩(wěn)定平臺式航空重力異常數(shù)據(jù)，且數(shù)據(jù)處理精度達(dá)到國際先進(jìn)水平。參考文獻(xiàn)[8]將經(jīng)驗(yàn)?zāi)B(tài)分解應(yīng)用到航空重力測量數(shù)據(jù)處理中，對已經(jīng)過FIR 低通濾波的航空重力測量數(shù)據(jù)進(jìn)行深化處理，可以很好地消除由飛機(jī)動態(tài)運(yùn)動引起的動態(tài)誤差，進(jìn)一步提高測量精度。針對一直存在的動態(tài)效應(yīng)剩余影響，黃謨濤等[9]基于AIC信息量準(zhǔn)則和互相關(guān)分析方法提出了一種針對各類動態(tài)效應(yīng)剩余誤差的通用補(bǔ)償模型，內(nèi)符合精度從原來的9.35mGal提升到1.01mGal。為了削弱各類誤差源的影響，歐陽永忠等[10]提出了一種相關(guān)分析法和測線網(wǎng)平差兩階段綜合誤差補(bǔ)償方法，實(shí)測數(shù)據(jù)處理結(jié)果驗(yàn)證了該方法的有效性和可靠性。

通過測線網(wǎng)平差求解系統(tǒng)誤差改正量，對各測線重力測量值進(jìn)行補(bǔ)償是處理航空重力測量數(shù)據(jù)的重要方法[11]，但現(xiàn)階段平差方法往往僅適用于規(guī)則的測線網(wǎng)[12]。航空重力測量中常常采用重復(fù)測線的測量方式來檢驗(yàn)航空重力儀動態(tài)測量的重復(fù)一致性，此時測線網(wǎng)平差的方法便不再適用[13]。蔡劭琨等[14]提出了一種航空重力測量重復(fù)測線系統(tǒng)誤差自調(diào)整的方法，也僅僅可以解決兩條重復(fù)測線的數(shù)據(jù)處理，怎樣對多條重復(fù)測線進(jìn)行系統(tǒng)平差還需要繼續(xù)深入研究。

因此，本文提出了基于經(jīng)驗(yàn)?zāi)B(tài)分解（EMD）和半系統(tǒng)誤差調(diào)整的航空重力測量重復(fù)測線數(shù)據(jù)處理方法。首先，推導(dǎo)了航空重力測量的基本原理，并詳細(xì)分析了測量噪聲的來源及誤差特性。其次，提出了基于EMD的誤差分離方法，對航空重力測量數(shù)據(jù)的深化處理進(jìn)一步提高動態(tài)重力測量的精度。最重要的是，基于半系統(tǒng)誤差調(diào)整模型提出了適用于多條重復(fù)測線的系統(tǒng)平差方法，填補(bǔ)了該領(lǐng)域的研究空白。最后，通過對實(shí)測數(shù)據(jù)的處理，驗(yàn)證了本文所提數(shù)據(jù)處理方法的精度和有效性。

1 航空重力測量數(shù)學(xué)模型

航空重力測量的兩個基本問題是如何保持重力敏感器的穩(wěn)定指向及如何從慣性加速度中分離出重力加速度[15]。平臺式重力儀基于跟蹤當(dāng)?shù)氐乩硐档姆€(wěn)定平臺，可以隔絕載體角運(yùn)動的影響，可以有效解決第一個問題。第二個問題的解決依賴于差分GPS，由其可以獲取高精度的位置和速度信息，用以計算載體垂向加速度以及重力相關(guān)的改正項(xiàng)，進(jìn)而分離得到重力信息。

1.1 數(shù)學(xué)原理

式（1）所述的比力方程是動態(tài)重力測量和慣性導(dǎo)航系統(tǒng)共同的理論基礎(chǔ)

式中，δg為重力異常；f為重力敏感器測量得到的比力；h為載體高度信息，h?為載體垂向加速度；γ為正常重力矢量的垂向分量；ωie為地球自轉(zhuǎn)角速度；VE，VN分別為載體相對于地球的東向和北向速度；φ為載體緯度信息；RN，RM分別為地球參考橢球的卯酉圈和子午圈的曲率半徑；后三項(xiàng)合記為厄特弗斯改正項(xiàng)δaE。

除厄特弗斯改正項(xiàng)外，航空重力測量還涉及正常重力改正、空間改正、偏心改正、水平加速度改正及垂直加速度改正等多個改正項(xiàng)。

1.2 噪聲誤差源分析

以飛機(jī)為載體的航空動態(tài)重力測量數(shù)據(jù)受到大量具有不同頻譜特性且分布廣泛的噪聲影響，本節(jié)依據(jù)噪聲特性進(jìn)行分類簡述。

1.2.1 高頻噪聲

在實(shí)際航空重力測量中，重力敏感器的輸出包含由飛機(jī)發(fā)動機(jī)以及空氣湍流造成的隨機(jī)振動加速度，這部分噪聲的強(qiáng)度是重力異常的數(shù)萬倍，屬于高頻噪聲。因此，需要利用低通濾波的方法進(jìn)一步去除重力異常信息中的高頻噪聲。

1.2.2 未知特性噪聲

測量條件中各種隨機(jī)因素產(chǎn)生的偶然誤差具有偶然性和隨機(jī)性，由載體動態(tài)性、空氣擾動等因素產(chǎn)生的噪聲也無法明確其誤差特性?；诮?jīng)驗(yàn)?zāi)B(tài)分解的數(shù)據(jù)深化處理是消除這些誤差的有效方法。

1.2.3 低頻誤差

受測量條件中某些特定因素如長周期傳感器誤差的系統(tǒng)性影響而產(chǎn)生的低頻誤差被稱為系統(tǒng)誤差，其與重力信息在同一頻段，無法通過低通濾波予以消除。系統(tǒng)誤差對觀測成果具有累積作用，因此在測量過程中常使用平差方法予以消除。

2 基于EMD的誤差分離

構(gòu)建補(bǔ)償模型是抑制重力測量動態(tài)效應(yīng)影響的一種方法，但模型構(gòu)建的科學(xué)性及模型參數(shù)的準(zhǔn)確性嚴(yán)重制約上述方法的有效性?；跀?shù)據(jù)處理的誤差分離是另一種有效的抑制方法。本節(jié)提出EMD 對已經(jīng)低通濾波過后的航空重力測量數(shù)據(jù)進(jìn)行深化處理，削弱未知特性噪聲對重力測量結(jié)果的影響，進(jìn)一步提高動態(tài)重力測量的精度。

2.1 EMD基本原理

EMD 方法認(rèn)為，任何復(fù)雜的時間序列x(t)都可以分解為一組從高頻到低頻的本征模態(tài)函數(shù)（IMF）和一個殘余信號r(t)，這些本征模態(tài)函數(shù)可以很好地反映信號在任何時間局部的頻率特征[16]。分解結(jié)果如式（4）所示

每個IMF分量都滿足兩個條件：一是整個信號上極值點(diǎn)數(shù)目與零點(diǎn)數(shù)目相差為1；二是極大值與極小值分別構(gòu)成的包絡(luò)線均值為0。EMD的分解步驟為：（1）確定原序列x(t)中的極大值和極小值點(diǎn)，并采用三次樣條分別對極大值和極小值點(diǎn)進(jìn)行插值，構(gòu)造上下包絡(luò)線xmax(t)和xmin(t)，可得上下包絡(luò)線的均值為m(t)；（2）計算原時間序列與包絡(luò)線均值序列的差值為h(t) =x(t)-m(t)；（3）若h(t)滿足IMF的兩個基本條件，則將h(t)設(shè)定為一個IMF分量，并求出剩余分量r(t) =x(t)-h(t)，繼續(xù)執(zhí)行步驟（4）；若不滿足，將h(t)替代x(t)返回執(zhí)行步驟（1）～（3）；（4）將r(t)作為新的序列重復(fù)執(zhí)行步驟（1）～（4）；循環(huán)執(zhí)行直到提取的最后一個IMF分量或者余項(xiàng)小于預(yù)先設(shè)定的閾值或已成為單調(diào)函數(shù)即可結(jié)束算法。

2.2 誤差分離

基于EMD對航空重力測量數(shù)據(jù)進(jìn)行深化處理的誤差分離，首先對低通濾波后去除高頻噪聲的重力數(shù)據(jù)進(jìn)行EMD分解，而后根據(jù)重力數(shù)據(jù)與IMF分量的相關(guān)性進(jìn)行篩選并重構(gòu)得到去噪后的信號。誤差分離的流程如圖1所示。

圖1 基于EMD的誤差分離流程Fig.1 Error separation flowchart based on EMD

若IMF 分量與重力數(shù)據(jù)的相關(guān)系數(shù)r（t）> 0.7，則認(rèn)定為強(qiáng)相關(guān)；若0.2

3 重復(fù)測線系統(tǒng)誤差調(diào)整

航空重力測量中各條重復(fù)測線在同一測點(diǎn)上的重力異常值理論上應(yīng)該是相同的，但由于受到系統(tǒng)誤差的影響，實(shí)際測量值并不相等，存在測點(diǎn)不符值。本節(jié)基于測區(qū)半系統(tǒng)差調(diào)整的算法提出適用于重復(fù)測線的平差方法，對測點(diǎn)數(shù)據(jù)進(jìn)行補(bǔ)償，抑制系統(tǒng)誤差的影響。

n條重復(fù)測線上各選取m個測點(diǎn)進(jìn)行平差計算，重復(fù)測線上重力測量值分別為g1j，g2j，…，gnj，其中j表示某一測點(diǎn)，且j∈[1,m]。將上述測量值用以下模型表示。

式中，g0是j測點(diǎn)上的重力真值；εij和Δij分別表示測線1～n上的系統(tǒng)誤差和隨機(jī)誤差。對n條測線上的重力測量值計算均方誤差，其中測線i的均方誤差計算式為

選取均方誤差最小的測線上的測量值為重復(fù)測線的重力基準(zhǔn)值gˉij。因此可得n條重復(fù)測線j測點(diǎn)上的重力不符值為

除去基準(zhǔn)測線外的（n-1）條重復(fù)測線測點(diǎn)不符值的均值作為該測點(diǎn)的總體不符值，并以該值作為平差的指標(biāo)，迭代計算至測點(diǎn)總體不符值代數(shù)和為0 時結(jié)束調(diào)整，并按照該方法遍歷重復(fù)測線上所有測點(diǎn)即可完成重復(fù)測線的系統(tǒng)平差處理。

4 實(shí)測數(shù)據(jù)分析

在某次航空重力測量中4條重復(fù)測線上獲得的重力異常結(jié)果如圖2所示，所示結(jié)果已通過相同的FIR低通濾波器將高頻噪聲濾除。由圖2 可以看出，重力異常數(shù)據(jù)低通濾波后相對平滑，高頻噪聲濾除效果較好，但仍然存在明顯的振蕩和測點(diǎn)不符值，計算得到內(nèi)符合精度為2.62mGal。

圖2 重復(fù)測線重力異常結(jié)果圖Fig.2 Gravity anomaly result map of repeated survey lines

4.1 EMD分解

前文已經(jīng)介紹了EMD分解的原理，下面首先對重復(fù)測線1 的重力異常結(jié)果進(jìn)行EMD 分解，得到如圖3 所示的分解結(jié)果。進(jìn)行相關(guān)性分析后選擇相關(guān)系數(shù)大于0.7 的IMF分量對重力異常信息進(jìn)行重構(gòu)，得到EMD去噪前后重力異常結(jié)果對比圖，如圖4 所示，可知EMD 深化處理后重力異常結(jié)果中的振蕩得到了很好的分離。

圖3 reline1重力異常結(jié)果EMD分解結(jié)果圖Fig.3 EMD decomposition result diagram of reline1 gravity anomaly results

圖4 reline1重力異常數(shù)據(jù)深化處理前后對比Fig.4 Comparison between reline1 gravity anomaly data before and after deepening processing

使用同樣的方法對其他測線進(jìn)行EMD分解并重構(gòu)，得到EMD深化處理后的重復(fù)測線重力異常測量結(jié)果，如圖5所示。重復(fù)測線數(shù)據(jù)EMD處理后消除了原有明顯的振蕩，計算可得內(nèi)符合精度為0.75mGal，但測點(diǎn)不符值依然存在，有待進(jìn)一步平差處理。

圖5 EMD去噪后的重復(fù)測線重力異常結(jié)果圖Fig.5 Gravity anomaly results of repeated survey lines after EMD denoising

4.2 重復(fù)測線平差調(diào)整

在重復(fù)測線上等間隔選取756 個測點(diǎn)，按照前文設(shè)計的重復(fù)測線平差算法進(jìn)行系統(tǒng)誤差調(diào)整。本文設(shè)計的系統(tǒng)平差算法精度與平差次數(shù)呈正相關(guān)，兩者的相關(guān)關(guān)系如圖6所示。由圖6 可知，迭代計算至第16 次后測點(diǎn)不符值精度已收斂，756 個測點(diǎn)不符值代數(shù)和以及標(biāo)準(zhǔn)差基本為0，此時認(rèn)為平差結(jié)束。

圖6 平差算法精度與平差次數(shù)關(guān)系圖Fig.6 Relationship between adjustment algorithm accuracy and adjustment times

平差結(jié)束后的重復(fù)測線重力異常結(jié)果如圖7 所示，相較平差前的結(jié)果，單個測點(diǎn)的不符值仍然存在，似乎表示該算法的優(yōu)化作用不是特別明顯，但通過計算可知，此時重復(fù)測線的內(nèi)符合精度提升至0.28mGal，相較平差前的0.75mGal已有大幅提升。這是因?yàn)楸疚乃岬闹貜?fù)測線平差方法是針對總體測線而言的，因此重復(fù)測線總體的內(nèi)符合精度提升較為明顯。

圖7 平差后重復(fù)測線重力異常結(jié)果圖Fig.7 Gravity anomaly results of repeated survey lines after adjustment

圖8給出了平差前后測點(diǎn)不符值的分布情況。顯而易見，平差結(jié)束后測點(diǎn)不符值誤差明顯減小，系統(tǒng)誤差得到有效補(bǔ)償。且測點(diǎn)不符值分布大致滿足標(biāo)準(zhǔn)正態(tài)分布，說明了各測點(diǎn)之間相互獨(dú)立。

圖8 平差前后測點(diǎn)不符值分布圖Fig.8 Distribution map of inconsistent values of measurement points before and after adjustment

4.3 結(jié)果分析

基于EMD分解有效分離了信號與誤差，本文提出的重復(fù)測線平差方法在總體上對系統(tǒng)誤差進(jìn)行調(diào)整，進(jìn)一步提高了總體測線的內(nèi)符合精度。重復(fù)測線內(nèi)符合精度統(tǒng)計見表1，由此可知本文的研究內(nèi)容顯著提升了數(shù)據(jù)處理的精度。

表1 重復(fù)測線內(nèi)符合精度統(tǒng)計Table 1 Statistical table of inner coincidence accuracy ofrepeated lines

5 結(jié)論

本文針對航空動態(tài)重力測量噪聲分布廣泛的問題，對重復(fù)測線的誤差分離展開研究，提出了基于EMD和半系統(tǒng)誤差調(diào)整的航空重力測量重復(fù)測線數(shù)據(jù)處理方法。實(shí)測數(shù)據(jù)的處理結(jié)果表明，EMD 的深化處理將內(nèi)符合精度從2.62mGal 提升至0.75mGal，有效分離出噪聲的影響；經(jīng)平差處理后抑制系統(tǒng)誤差的影響，內(nèi)符合精度進(jìn)一步提升至0.28mGal，精度依次提升了71.37%和62.67%。本文研究內(nèi)容對進(jìn)一步推動航空動態(tài)重力測量精度提升和工程應(yīng)用化具有積極意義。