基于SINS/GNSS的航空矢量重力測(cè)量數(shù)據(jù)處理方法研究

寧津生

（1.武漢大學(xué)測(cè)繪學(xué)院，武漢 430079；2.武漢大學(xué)地球空間環(huán)境與大地測(cè)量教育部重點(diǎn)實(shí)驗(yàn)室，武漢 430079）

基于SINS/GNSS的航空矢量重力測(cè)量數(shù)據(jù)處理方法研究

寧津生1，2

（1.武漢大學(xué)測(cè)繪學(xué)院，武漢 430079；2.武漢大學(xué)地球空間環(huán)境與大地測(cè)量教育部重點(diǎn)實(shí)驗(yàn)室，武漢 430079）

航空矢量重力測(cè)量是未來高效測(cè)定中高頻地球重力場(chǎng)信息的先進(jìn)技術(shù)，與地面重力測(cè)量、海洋重力測(cè)量和衛(wèi)星重力測(cè)量技術(shù)相互補(bǔ)充。本文介紹了基于捷聯(lián)慣性導(dǎo)航系統(tǒng)/全球衛(wèi)星導(dǎo)航系統(tǒng)（SINS/ GNSS）的航空矢量重力測(cè)量的原理，重點(diǎn)討論了航空矢量重力測(cè)量的數(shù)據(jù)預(yù)處理、數(shù)據(jù)歸算、帶限地形影響、向下延拓及大地水準(zhǔn)面確定等問題，為發(fā)展我國(guó)航空矢量重力測(cè)量技術(shù)提供參考。

航空矢量重力測(cè)量；SINS/GNSS數(shù)據(jù)融合；地球重力場(chǎng)；大地水準(zhǔn)面

1 前言

航空重力測(cè)量技術(shù)是以飛機(jī)為載體，快速測(cè)定近地空中重力加速度的重力測(cè)量方法。航空重力測(cè)量最早出現(xiàn)在20世紀(jì)50年代，經(jīng)過數(shù)十年的沉淀與發(fā)展，已成為高效測(cè)定中高頻地球重力場(chǎng)信息的主要手段。在地面重力測(cè)量難以到達(dá)的地區(qū)進(jìn)行航空重力測(cè)量，可以有效填補(bǔ)困難地區(qū)的重力空白，有利于改善重力場(chǎng)數(shù)據(jù)的精度和分辨率及精化局部重力場(chǎng)模型和區(qū)域大地水準(zhǔn)面。航空重力測(cè)量結(jié)合其他技術(shù)手段獲取的高精度地球重力場(chǎng)信息，在國(guó)家經(jīng)濟(jì)建設(shè)、軍事和國(guó)防建設(shè)、地球物理、海洋學(xué)、地球動(dòng)力學(xué)、資源勘探等相關(guān)地球科學(xué)領(lǐng)域具有重要作用。

以美國(guó)、加拿大、丹麥等為代表的發(fā)達(dá)國(guó)家先后開展了航空重力測(cè)量研究，并進(jìn)行了大量的飛行試驗(yàn)[1～3]，我國(guó)也在2002年成功研制了自己的航空標(biāo)量測(cè)量系統(tǒng)Chinese airborne gravimetry system（CHAGS）[4]。眾多研究結(jié)果表明，航空標(biāo)量重力測(cè)量數(shù)據(jù)在波長(zhǎng)5～10 km的尺度上能夠達(dá)到1～3 mGal（1 Gal=0.01 m/s2）的精度，這也意味著航空標(biāo)量重力測(cè)量技術(shù)已經(jīng)進(jìn)入成熟階段。航空矢量重力測(cè)量技術(shù)不僅能獲取航空標(biāo)量重力測(cè)量觀測(cè)值（即重力矢量的垂直分量），也能測(cè)得重力矢量的水平分量，是目前大地測(cè)量領(lǐng)域的研究熱點(diǎn)之一。Jekeli和Garcia[5]驗(yàn)證了全球衛(wèi)星導(dǎo)航系統(tǒng)（GNSS）相位加速度應(yīng)用于航空矢量重力測(cè)量的可行性；Kwon和Jekeli[6，7]等率先啟動(dòng)關(guān)于航空矢量重力測(cè)量的研究與試驗(yàn)并提出了利用慣性導(dǎo)航系統(tǒng)（INS）/GNSS組合系統(tǒng)測(cè)定矢量重力的新方法，所得重力矢量垂直分量精度為3 mGal，水平分量精度為6～8 mGal，分辨率為10 km；Senobari[8]給出了基于捷聯(lián)慣性導(dǎo)航系統(tǒng)（SINS）/GNSS的航空矢量重力測(cè)量的試驗(yàn)結(jié)果；美國(guó)政府也于2012年開始為期10年的GRAV-D計(jì)劃，并將航空矢量重力測(cè)量作為其中的一項(xiàng)關(guān)鍵技術(shù)，預(yù)期確定北美地區(qū)大地水準(zhǔn)面的精度達(dá)到2 cm[9]。目前國(guó)內(nèi)還沒有航空矢量重力測(cè)量系統(tǒng)，在航空矢量重力數(shù)據(jù)處理方面也尚處于起步階段。因此，充分吸收國(guó)內(nèi)外現(xiàn)有研究成果，深入研究航空矢量重力測(cè)量的基礎(chǔ)理論和數(shù)據(jù)處理方法，對(duì)我國(guó)航空矢量重力測(cè)量技術(shù)的發(fā)展具有重要的科學(xué)意義和應(yīng)用價(jià)值。

2 基本原理

依據(jù)牛頓第二定律，單位質(zhì)點(diǎn)的加速度在數(shù)值上等于質(zhì)點(diǎn)所受作用力的合力。在慣性系i系中質(zhì)點(diǎn)的運(yùn)動(dòng)方程可寫為

ai=fi+gi（1）

式（1）中，ai為慣性加速度，由GNSS觀測(cè)值給出；fi為比力觀測(cè)值，由SINS給出；gi為重力加速度。將重力加速度矢量表示為正常重力矢量γn與擾動(dòng)重力δgn之和，則在導(dǎo)航坐標(biāo)系n系中，航空矢量重力測(cè)量的基本模型可表示為

式（2）中，右邊第三項(xiàng)為Coriolis加速度；Ωnie為地心地固坐標(biāo)系e系相對(duì)于i系的旋轉(zhuǎn)角速度；Ωnen為n系相對(duì)于e系的旋轉(zhuǎn)角速度。對(duì)基于SINS/GNSS的航空矢量重力測(cè)量系統(tǒng)來說，加速度計(jì)和陀螺儀固定在載體上，比力是在載體坐標(biāo)系b系中測(cè)得的，因此需要利用陀螺儀數(shù)據(jù)將比力轉(zhuǎn)換至導(dǎo)航坐標(biāo)系n系，此時(shí)式可寫為

式（3）中，Rnb為b系到n系的旋轉(zhuǎn)矩陣。

3 數(shù)據(jù)處理方法

3.1 數(shù)據(jù)預(yù)處理

航空矢量重力測(cè)量是一個(gè)動(dòng)態(tài)測(cè)量過程，各類觀測(cè)數(shù)據(jù)中均包含了大量的觀測(cè)噪聲，如何從低信噪比的加速度觀測(cè)信息中有效地提取重力信息是航空矢量重力測(cè)量數(shù)據(jù)處理的關(guān)鍵。這一過程又大致可分為數(shù)據(jù)去噪、誤差補(bǔ)償、數(shù)據(jù)融合、低通濾波等，具體數(shù)據(jù)處理流程如圖1所示。在數(shù)據(jù)預(yù)處理過程中對(duì)原始觀測(cè)數(shù)據(jù)進(jìn)行去噪和誤差補(bǔ)償，可以提高數(shù)據(jù)的信噪比，進(jìn)而提高數(shù)據(jù)處理精度。慣性測(cè)量單元（inertial measurement unit，IMU）是航空矢量重力測(cè)量系統(tǒng)的重要組成部分，主要由加速度計(jì)和陀螺儀構(gòu)成。由于IMU自身的誤差特性，在外業(yè)測(cè)量前需對(duì)其進(jìn)行外場(chǎng)標(biāo)定，提高數(shù)據(jù)質(zhì)量，減小系統(tǒng)誤差對(duì)后續(xù)數(shù)據(jù)處理帶來的影響[10]。

圖1 航空矢量重力測(cè)量的數(shù)據(jù)預(yù)處理Fig.1 Data preprocessing of airborne vector gravimetry

初始對(duì)準(zhǔn)是SINS在進(jìn)行各種應(yīng)用前所必需的一項(xiàng)重要工作，計(jì)算出b系與n系或者i系之間的姿態(tài)轉(zhuǎn)移矩陣，其對(duì)準(zhǔn)精度對(duì)后續(xù)的各種應(yīng)用產(chǎn)生直接影響，而Kalman濾波技術(shù)是實(shí)現(xiàn)初始對(duì)準(zhǔn)的一種有效方法。最優(yōu)的Kalman濾波器是建立在函數(shù)模型和噪聲特性均精確已知的條件下，但噪聲統(tǒng)計(jì)特性一般難以獲得，自適應(yīng)Kalman濾波是解決這一應(yīng)用難題的主要方法之一。根據(jù)IMU的自身特性以及作業(yè)場(chǎng)地等因素，加速度計(jì)的噪聲特性通常是未知的，本文提出了基于簡(jiǎn)化自協(xié)方差最小二乘（autocovariance least-squares，ALS）噪聲估計(jì)的SINS靜基座初始對(duì)準(zhǔn)算法（SALS），計(jì)算流程見圖2。

圖2 基于簡(jiǎn)化ALS噪聲估計(jì)的SINS靜基座初始對(duì)準(zhǔn)流程圖Fig.2 Flow chart of the simplified ALS for SINS initial alignment on stationary base

圖3給出了數(shù)值仿真中東、北、天方向姿態(tài)角的初始對(duì)準(zhǔn)誤差，從中可以看出，天方向姿態(tài)角比水平方向姿態(tài)角更易受觀測(cè)噪聲協(xié)方差矩陣的影響，當(dāng)先驗(yàn)觀測(cè)噪聲協(xié)方差矩陣存在較大偏差時(shí)，對(duì)天方向?qū)?zhǔn)結(jié)果有較大影響；本文提出的基于簡(jiǎn)化ALS噪聲估計(jì)的SINS SALS同時(shí)進(jìn)行噪聲估計(jì)和姿態(tài)角修正，因此各方向姿態(tài)角的對(duì)準(zhǔn)精度較常規(guī)Kalman濾波方法有著明顯得提高。靜基座初始對(duì)準(zhǔn)的Kalman濾波模型并非完全可觀，三個(gè)失準(zhǔn)角由于受陀螺常值漂移和加速度計(jì)常值漂移的影響，故均存在估計(jì)偏差。

圖3 初始對(duì)準(zhǔn)誤差Fig.3 Errors of initial alignment

載體的運(yùn)動(dòng)加速度是最重要的改正項(xiàng)，其精度直接影響著航空重力測(cè)量的精度和空間分辨率[11，12]，但載體加速度并非直接觀測(cè)值，它是通過載體位置或速度差分計(jì)算得到的。理論上，已知載體位置的時(shí)間序列，采用數(shù)值差分的方法便可得到載體加速度的時(shí)間序列，但由于數(shù)值差分過程會(huì)放大高頻噪聲，在實(shí)際數(shù)據(jù)處理中，還需進(jìn)行低通濾波，然而當(dāng)噪聲和信號(hào)在同一頻帶上有交疊時(shí)，采用低通濾波方法必然會(huì)造成信號(hào)損失，且不同的低通濾波器及其參數(shù)的選擇也會(huì)對(duì)信號(hào)處理結(jié)果有一定的影響。Kalman濾波方法是一種基于模型的數(shù)字濾波方法，可以用來確定載體的加速度。在進(jìn)行航空重力測(cè)量時(shí)，一般要求飛行過程比較平穩(wěn)，因此可采用常加速度模型對(duì)載體運(yùn)動(dòng)狀態(tài)進(jìn)行建模。由于Kalman濾波模型的噪聲協(xié)方差矩陣特別是狀態(tài)噪聲協(xié)方差矩陣難以獲取，不適當(dāng)?shù)脑肼晠?shù)的選取必然會(huì)對(duì)Kalman濾波的結(jié)果帶來影響。而顧及常加速度模型的狀態(tài)噪聲協(xié)方差矩陣滿足特定結(jié)構(gòu)[13]，本文提出了一種基于ALS噪聲估計(jì)的載體加速度確定算法，其計(jì)算流程見圖4。

圖4 基于ALS噪聲估計(jì)的載體加速度確定算法流程圖Fig.4 Flow chart of the ALS method for acceleration estimation

表1中5種不同方案的狀態(tài)估計(jì)結(jié)果表明：a.采用9點(diǎn)差分法計(jì)算得到的速度和加速度精度最差；b.當(dāng)先驗(yàn)噪聲協(xié)方差存在較大偏差時(shí)，Kalman濾波精度較差，再采用RTS算法進(jìn)行平滑處理，其狀態(tài)估計(jì)精度提高有限；c.先采用ALS進(jìn)行噪聲估計(jì)，再進(jìn)行RTS平滑，其狀態(tài)估計(jì)精度有著明顯提高。圖5說明差分法的差分誤差在高頻部分特別顯著，采用ALS及RTS平滑算法不僅在高頻段對(duì)誤差有較好的抑制作用，在低頻段對(duì)信號(hào)的擬合精度也明顯優(yōu)于差分法。

表1 加速度計(jì)算結(jié)果比較Table 1 Results comparison of acceleration

圖5 加速度的頻譜圖Fig.5 Spectrum of acceleration

INS具有較高的短期穩(wěn)定性，但誤差隨著時(shí)間累積，而GNSS具有較高的精度，但可能存在數(shù)據(jù)空白?；趦烧叩幕パa(bǔ)性，通過差分GNSS得到載體的運(yùn)動(dòng)加速度，將載體加速度與比力及正常重力的差值作為觀測(cè)信息來建立Kalman濾波模型，并對(duì)SINS的各種誤差進(jìn)行估計(jì)，此時(shí)擾動(dòng)重力矢量存在于Kalman濾波的觀測(cè)殘差中，從而可提取到測(cè)線上的擾動(dòng)重力。Kalman濾波是SINS/GNSS數(shù)據(jù)融合時(shí)最為常用的數(shù)據(jù)處理方法，但其有效性有待在實(shí)際數(shù)據(jù)處理中得到驗(yàn)證。

數(shù)據(jù)融合后得到的測(cè)線擾動(dòng)重力中仍含有大量的高頻噪聲，大量研究和試驗(yàn)表明，低通濾波技術(shù)是消除高頻噪聲的有效方法，國(guó)內(nèi)外學(xué)者針對(duì)航空標(biāo)量重力測(cè)量數(shù)據(jù)的低通濾波做了大量的研究[13～18]。本文采用窗函數(shù)法（濾波器1）和最佳一致逼近法（濾波器2）設(shè)計(jì)了適用于航空矢量重力測(cè)量數(shù)據(jù)處理的FIR數(shù)字低通濾波器，結(jié)合GPS靜態(tài)測(cè)量數(shù)據(jù)進(jìn)行了數(shù)值試驗(yàn)。結(jié)果如圖6所示，說明兩類濾波器均能以±1～2 mGal的精度確定垂直加速度，以高于±1 mGal的精度確定水平加速度。

圖6 濾波后GPS測(cè)站的加速度Fig.6 Accelerations of GPS station after lowpass filtering

3.2 數(shù)據(jù)歸算

為滿足地球重力場(chǎng)有關(guān)研究及相關(guān)應(yīng)用的需要，航空矢量重力觀測(cè)值往往被歸算成平均飛行高度面上規(guī)則格網(wǎng)化的擾動(dòng)重力，這一過程包括平均高度面歸算、測(cè)線網(wǎng)平差和格網(wǎng)化數(shù)據(jù)生成。平均高度面歸算本質(zhì)上可理解為實(shí)測(cè)數(shù)據(jù)相對(duì)于平均高度面的向上延拓或向下延拓改正，因此平均高度面歸算可采用重力延拓方法，如譜方法[19]、泊松積分法[20]和徑向改正方法[21]等，當(dāng)實(shí)際航線相對(duì)于平均高度面的變化較小時(shí)，一般只考慮徑向改正方法。測(cè)線網(wǎng)平差實(shí)際上是一個(gè)系統(tǒng)誤差的補(bǔ)償問題，進(jìn)行系統(tǒng)誤差補(bǔ)償?shù)囊话惴椒樽詸z校平差法，但其設(shè)計(jì)矩陣是秩虧的，這時(shí)可選定某一條或幾條精度很高的測(cè)線作為固定測(cè)線，將這些測(cè)線的系統(tǒng)誤差作為虛擬觀測(cè)值，按秩虧自由網(wǎng)平差法進(jìn)行求解[22]。但在很多情況下，因觀測(cè)區(qū)域的地面檢校數(shù)據(jù)不足，觀測(cè)數(shù)據(jù)的精度無法準(zhǔn)確評(píng)定，此時(shí)可采用兩步平差法來進(jìn)行系統(tǒng)誤差補(bǔ)償[23～25]。航空重力測(cè)量數(shù)據(jù)歸算的目的是獲得平均飛行高度面上具有一定空間分辨率、按經(jīng)緯度劃分的規(guī)則格網(wǎng)重力數(shù)據(jù)，格網(wǎng)化方法的好壞對(duì)格網(wǎng)化輸出數(shù)據(jù)的質(zhì)量、精度和可信度有著直接影響。重力數(shù)據(jù)處理中常用的格網(wǎng)化方法主要有加權(quán)平均法、Shepard曲面擬合法、Kriging方法和最小二乘配置法（LSC），這些格網(wǎng)化方法有著各自的優(yōu)缺點(diǎn)，在實(shí)際應(yīng)用中有所差異。

3.3 帶限地形影響

地形起伏反映了地球重力場(chǎng)的不規(guī)則變化，精確確定地形對(duì)地球重力場(chǎng)的影響具有十分重要的意義[26]。第二類Helmert凝集法是目前應(yīng)用最為廣泛的地形歸算方法，分為解析形式和譜形式，前者提供全波段的地形質(zhì)量改正，后者則是將牛頓積分核函數(shù)展開成Legendre多項(xiàng)式的收斂級(jí)數(shù)，可用來計(jì)算帶限地形影響。

航空重力測(cè)量的各類觀測(cè)數(shù)據(jù)中包含了大量的高頻觀測(cè)噪聲，消除這些噪聲的主要手段是使用低通濾波器，濾波后得到的重力信號(hào)為帶限重力信號(hào)[27]。因此，采用解析公式計(jì)算的全頻譜地形直接影響和間接影響對(duì)航空矢量重力觀測(cè)值進(jìn)行改正是不合理的，解決這一問題的方法有兩類：a.將地形改正值也作相應(yīng)的低通濾波處理，使得地形效應(yīng)的頻譜范圍與航空重力觀測(cè)值一致[28]，這種方法需要輸入飛行高度上密集采樣點(diǎn)處的地形質(zhì)量影響值，采用高分辨率的數(shù)字高程模型（digital elevation model，DEM）來進(jìn)行計(jì)算時(shí)工作量相當(dāng)龐大；b.基于第二類Helmert凝集法的譜形式，采用帶限公式來計(jì)算地形質(zhì)量對(duì)引力位和重力擾動(dòng)的直接影響公式和大地水準(zhǔn)面的間接影響[29]。這種方法的優(yōu)點(diǎn)是可依據(jù)重力觀測(cè)值的頻譜范圍來調(diào)整積分核函數(shù)的起始階數(shù)和最高階數(shù)，從而避免了對(duì)地形質(zhì)量影響作低通濾波處理，原理簡(jiǎn)單而且計(jì)算方便。Novák等[29]基于第二類Helmert凝集法的譜形式推導(dǎo)了地形質(zhì)量對(duì)引力位和重力擾動(dòng)的帶限直接影響改正公式以及對(duì)大地水準(zhǔn)面的帶限間接影響改正公式，將基于級(jí)數(shù)核計(jì)算的帶限直接地形影響和間接地形影響與基于解析核計(jì)算并經(jīng)低通濾波處理后得到的帶限直接地形影響和間接地形影響進(jìn)行了比較，驗(yàn)證了算法的有效性，該帶限地形影響改正公式可對(duì)重力擾動(dòng)或重力異常進(jìn)行地形歸算，針對(duì)的是航空標(biāo)量重力測(cè)量觀測(cè)值。航空重力矢量測(cè)量獲取的是擾動(dòng)重力的三個(gè)分量，而針對(duì)水平分量帶限地形影響改正的研究尚未有公開的文獻(xiàn)發(fā)表。本文在Novák研究的基礎(chǔ)上，推導(dǎo)了重力矢量水平分量的帶限直接地形影響改正公式，并在中國(guó)西部山區(qū)選取了兩條沿緯圈方向4 000 m（大地高）航線高度上的航線，采用3″×3″的數(shù)字高程模型基于解析核計(jì)算了水平分量的直接地形影響，經(jīng)低通濾波后與基于帶限公式計(jì)算的地形影響進(jìn)行了比較。數(shù)值計(jì)算結(jié)果如圖7～圖10所示，表明二者吻合較好。

3.4 向下延拓

航空矢量重力測(cè)量的觀測(cè)值為近地空中的中高頻重力信號(hào)，但地球物理和大地測(cè)量領(lǐng)域的許多應(yīng)用都需要地形表面或大地水準(zhǔn)面上的重力觀測(cè)值。因此，將航線上的重力測(cè)量觀測(cè)值進(jìn)行向下延拓，是航空矢量重力測(cè)量數(shù)據(jù)處理的主要內(nèi)容之一。

向下延拓是一個(gè)不適定的過程，屬于病態(tài)問題，會(huì)放大重力觀測(cè)數(shù)據(jù)中的噪聲，若不采用合適的延拓方法，將會(huì)導(dǎo)致延拓解的嚴(yán)重失真。在現(xiàn)有的向下延拓方法中，逆Possion積分法的應(yīng)用最為廣泛，逆Possion積分法輔以Tikhonov正則化[20]、嶺估計(jì)[30]、廣義嶺估計(jì)[30]、迭代法[31]、濾波[32，33]等正則化方法才能得到穩(wěn)定的延拓解，其關(guān)鍵在于選擇合理的正則化參數(shù)。Novák和Heck[34]認(rèn)為航空重力為帶限信號(hào)，提出了基于帶限Possion核函數(shù)的直接延拓法，從而避免了逆Possion積分法面臨的大型法方程陣求逆問題。LSC兼具向上延拓和向下延拓的雙重功能，在延拓過程中可融入新的重力數(shù)據(jù)且能有效估計(jì)和改善系統(tǒng)偏差[32]，但LSC延拓結(jié)果的好壞很大程度上取決于協(xié)方差函數(shù)模型。解析延拓也是重力信號(hào)向下延拓的一種常用方法，同樣兼具向上延拓和向下延拓的能力，其基本思想是假設(shè)重力信號(hào)的垂直梯度已知，則延拓表達(dá)式可表示為一個(gè)Taylor展開級(jí)數(shù)；解析延拓在空域中的求解十分復(fù)雜，Wu[35]給出了解析延拓在頻域中的求解公式。此外，一些數(shù)學(xué)方法如樣條函數(shù)法、小波算法等的引入，使得重力信號(hào)向下延拓的手段更加豐富。

圖7 沿31.5°N緯線4 km飛行高度處北向水平分量的直接地形影響Fig.7 Direct topographical effect for northern component at 4 km flight height along the parallel of 31.5°N

圖8 沿32.5°N緯線4 km飛行高度處北向水平分量的直接地形影響Fig.8 Direct topographical effect for northern component at 4 km flight height along the parallel of 32.5°N

圖9 沿31.5°N緯線4 km飛行高度處東向水平分量的直接地形影響Fig.9 Direct topographical effect for eastern component at 4 km flight height along the parallel of 31.5°N

以上方法一般用于航空標(biāo)量重力測(cè)量數(shù)據(jù)，即重力擾動(dòng)（重力異常）的向下延拓，而航空矢量重力測(cè)量不僅能獲取重力擾動(dòng)（重力矢量的垂直分量），還能測(cè)出重力矢量的水平分量。理論上水平分量的向下延拓原理與垂直分量類似，但現(xiàn)有公開發(fā)表的文獻(xiàn)中針對(duì)水平分量向下延拓的研究較少，且由于航空矢量重力測(cè)量技術(shù)自身的缺陷，空中獲取的水平分量的精度遠(yuǎn)不如垂直分量，僅為6～8 mGal，如何實(shí)現(xiàn)低信噪比情況下水平分量的向下延拓是一個(gè)需要解決的關(guān)鍵問題。本文引入了頻域輸入輸出系統(tǒng)，在不同的噪聲水平下比較分析了單輸入單輸出系統(tǒng)（single-input single-output system，SISOS）和雙輸入單輸出系統(tǒng)（double-input singleoutput system，DISOS）對(duì)水平分量的向下延拓效果（見圖11和圖12）。研究結(jié)果表明：a.當(dāng)水平分量的精度較高（如1.5 mGal）時(shí)，二者均能實(shí)現(xiàn)水平分量的穩(wěn)定向下延拓；b.當(dāng)水平分量的精度較低（如6 mGal）時(shí)，單輸入單輸出法的延拓結(jié)果較差，而雙輸入單輸出法能實(shí)現(xiàn)水平分量的穩(wěn)定向下延拓。

圖10 沿32.5°N緯線4 km飛行高度處東向水平分量的直接地形影響Fig.10 Direct topographical effect for eastern component at 4 km flight height along the parallel of 32.5°N

圖11 東向分量的精度為1.5 mGal時(shí)延拓誤差的空間分布Fig.11 Distribution of continued errors of estern component with the accuracy of 1.5 mGal

圖12 東向分量的精度為6 mGal時(shí)延拓誤差的空間分布Fig.12 Distribution of continued errors of estern component with the accuracy of 6 mGal

4 航空矢量重力應(yīng)用于大地水準(zhǔn)面的確定

利用航空矢量重力測(cè)量數(shù)據(jù)確定大地水準(zhǔn)面可歸結(jié)為求解航空矢量重力測(cè)量邊值問題。與傳統(tǒng)的重力學(xué)邊值問題不同，航空矢量重力測(cè)量邊值問題的邊界面為飛行高度面，可描述為[36]

式（4）中，Δ為拉普拉斯算子；T為擾動(dòng)位；H為航線高度。由式（4）可知，擾動(dòng)重力的水平分量和垂直分量均可用于精化區(qū)域大地水準(zhǔn)面。垂直分量在全波長(zhǎng)為5～10 km的分辨率尺度上能達(dá)到1～3 mGal的精度，已能滿足區(qū)域大地水準(zhǔn)面確定的要求。

利用垂直分量確定大地水準(zhǔn)面的方法主要有兩類。第一類為空域積分法，它可分為兩步法和一步法[37]。兩步法包括兩個(gè)計(jì)算步驟：a.求解位理論Dirichlet邊值問題，將飛行高度面上的垂直分量向下調(diào)和延拓至大地水準(zhǔn)面；b.求解位理論Neumann邊值問題，按照Hotine積分公式由大地水準(zhǔn)面上的垂直分量得到大地水準(zhǔn)面上的擾動(dòng)位，再由Bruns公式得到大地水準(zhǔn)面高。一步法是一種比兩步法更穩(wěn)定、精度更高的空域解法，它將向下延拓和求解大地水準(zhǔn)面進(jìn)行了綜合。由于航空重力觀測(cè)值為帶限信號(hào)，一步法采用頻譜范圍與觀測(cè)值一致的帶限積分核函數(shù)，由廣義Hotine積分公式結(jié)合Bruns公式即可得到大地水準(zhǔn)面高。第二類為譜方法。通過某種形式的逼近，Hotine積分、Possion積分可化為空域卷積形式。研究表明，采用快速傅里葉變換（fast Fourier transform，F(xiàn)FT）將Hotine積分和Possion積分轉(zhuǎn)換至頻域進(jìn)行計(jì)算可使積分法的計(jì)算速度得到極大的提升[38～40]。此外，還可采用球冠諧分析或矩諧分析將重力觀測(cè)數(shù)據(jù)表達(dá)為局域重力場(chǎng)的諧展開式。球冠諧分析[41]和矩諧分析[42]最早應(yīng)用于地磁場(chǎng)表達(dá)，Li等[43]首次將球冠諧分析引入?yún)^(qū)域重力場(chǎng)研究中，邊少峰[44]和蔣濤[30]則將矩諧分析用于局部重力場(chǎng)表達(dá)。對(duì)于大范圍的區(qū)域重力場(chǎng)逼近，球冠諧分析是一種較為理想的方法，而矩諧分析則適用于區(qū)域重力場(chǎng)精細(xì)結(jié)構(gòu)的逼近。

類似于天文水準(zhǔn)原理，可采用逐測(cè)線剖面積分將擾動(dòng)重力的水平分量轉(zhuǎn)化為航線上的擾動(dòng)位，將擾動(dòng)位向下延拓至大地水準(zhǔn)面后，再采用Bruns公式計(jì)算大地水準(zhǔn)面起伏[45]。值得注意的是，垂直分量確定的是絕對(duì)大地水準(zhǔn)面，而水平分量確定的是相對(duì)大地水準(zhǔn)面。

在聯(lián)合重力矢量的三個(gè)分量精化大地水準(zhǔn)面方面，隨機(jī)解法、頻域輸入輸出法、多分辨率配置法、譜組合法等方法可將航線的重力矢量轉(zhuǎn)化為航線上的擾動(dòng)位，將擾動(dòng)位向下延拓至大地水準(zhǔn)面，再采用Bruns公式計(jì)算大地水準(zhǔn)面起伏。目前航空矢量重力測(cè)量水平分量的精度水平遠(yuǎn)不如垂直分量，需要分析不同信噪比情況下的水平分量對(duì)三分量聯(lián)合求解大地水準(zhǔn)面的影響，進(jìn)一步研究聯(lián)合航空重力矢量的三個(gè)分量確定大地水準(zhǔn)面的融合方法。

5 結(jié)語(yǔ)

基于SINS/GNSS的航空矢量重力測(cè)量是未來高效測(cè)定中高頻地球重力場(chǎng)信息的先進(jìn)技術(shù)，聯(lián)合地面重力測(cè)量、海洋重力測(cè)量和衛(wèi)星重力測(cè)量技術(shù)，可以獲取全球高分辨率高精度地球重力場(chǎng)信息，從而滿足國(guó)家經(jīng)濟(jì)建設(shè)、軍事和國(guó)防建設(shè)、地球物理、海洋學(xué)、地球動(dòng)力學(xué)、資源勘探等相關(guān)地球科學(xué)領(lǐng)域的需求。本文在國(guó)內(nèi)外現(xiàn)有研究成果的基礎(chǔ)上，討論了航空矢量重力測(cè)量的數(shù)據(jù)預(yù)處理、數(shù)據(jù)歸算、帶限地形影響、向下延拓及大地水準(zhǔn)面確定等問題。

航空矢量重力測(cè)量技術(shù)的實(shí)用化還需要解決一系列關(guān)鍵問題，如高精度SINS的自主研制、SINS/ GNSS的數(shù)據(jù)融合方法及滿足工程化應(yīng)用的數(shù)據(jù)處理軟件、SINS和GNSS的誤差校正、地殼橫向密度擾動(dòng)對(duì)向下延拓和精化大地水準(zhǔn)面的影響、航空矢量重力及其他類型重力場(chǎng)信息的融合處理等。

[1]Wei M，Schwarz K P.Flight test results from a strapdown airborne gravity system[J].Journal of Geodesy，1998，72：323-332.

[2]Verdun J，Klingelé E，Bayer R，et al.The alpine Swiss-French airborne gravity survey[J].Geophysical Journal International，2003，152：8-19.

[3]Forsberg R，Olesen A，Vest A，et al.Gravity field improvements in the north atlantic region[C]//Proc.GOCE Workshop，ESA-ESRIN，2004.

[4]夏哲仁，石 磐，孫中苗，等.航空重力測(cè)量系統(tǒng)CHAGS[J].測(cè)繪學(xué)報(bào)，2004，33（3）：216-220.

[5]Jekeli C，Garcia R.GNSS phase accelerations for movingbase vector gravimetry[J].Journal of Geodesy，1997，71：630-639.

[6]Kwon J.Airborne vector gravimetry using GNSS/INS[D].USA：The Ohio State University，2000.

[7]Kwon J，Jekeli C.A new approach for airborne vector gravimetry using GNSS/INS[J].Journal of Geodesy，2001，74（10）：690-700.

[8]Senobari M S.New results in airborne vector gravimetry using strapdown INS/GNSS[J].J Geoid，2010，84：277-291.

[9]Smith D.The GRAV-D project：Gravity for the redefinition of the American Vertical Datum[OL].NOAA website：http：//www. ngs.noaa.gov/GRAV-D/pubs/GRAV-D_v2007_12_19.pdf，2007.

[10]Lee J G，Park C G，Park H W.Multiposition alignment of strapdown inertial navigation system[J].Aerospace and Electronic Systems，IEEE Transactions on，1993，29（4）：1323-1328.

[11]Kwon J H，Jekeli C.A new approach for airborne vector gravimetry using GNSS/INS[J].Journal of Geodesy，2001，74（10）：690-700.

[12]孫中苗.航空重力測(cè)量理論，方法及應(yīng)用研究[D].鄭州：解放軍信息工程大學(xué)，2004.

[13]Singer R A.Estimating optimal tracking filter performance for manned maneuvering targets[J].IEEE Transactions on Aerospace and Electronic Systems，1970（4）：473-483.

[14]Hammada Y.A comparison of filtering techniques for airborne gravimetry[D].Calgary：Department of Geomatics Engineering at University of Calgary，1996.

[15]孫中苗.航空重力測(cè)量中FIR低通濾波器的設(shè)計(jì)與應(yīng)用[J].解放軍測(cè)繪學(xué)院學(xué)報(bào)，1996，13（4）：247-250.

[16]柳林濤，許厚澤.航空重力測(cè)量數(shù)據(jù)的小波濾波處理[J].地球物理學(xué)報(bào)，2004，47（3）：490-494.

[17]孫中苗，翟振和.航空重力測(cè)量數(shù)據(jù)的小波閾值濾波[J].武漢大學(xué)學(xué)報(bào)（信息科學(xué)版），2009，34（10）：1222-1225.

[18]蔡劭琨，吳美平，張開東.航空重力測(cè)量中FIR低通濾波器的比較[J].物探與化探，2010，34（1）：74-78.

[19]蔣 濤，黨亞明，章傳銀，等.基于矩諧分析的航空重力向下延拓[J].測(cè)繪學(xué)報(bào)，2013，42（4）：475-480.

[20]王興濤，石 磐，朱非洲.航空重力測(cè)量數(shù)據(jù)向下延拓的正則化算法及其譜分解[J].測(cè)繪學(xué)報(bào)，2004，33（1）：33-38.

[21]鐘 波，劉華亮，羅志才，等.衛(wèi)星重力梯度數(shù)據(jù)的歸算與格網(wǎng)化處理[J].大地測(cè)量與地球動(dòng)力學(xué)，2011，31（3）：79-84.

[22]Hwang C，Hsiao Y S，Shih H C.Data reduction in scalar airborne gravimetry：Theory，software and case study in Taiwan [J].Computers&Geosciences，2006，32（10）：1573-1584.

[23]黃謨濤，翟國(guó)君，歐陽(yáng)永忠，等.海洋重力測(cè)量誤差補(bǔ)償兩步處理法[J].武漢大學(xué)學(xué)報(bào)（信息科學(xué)版），2002，27（3）：251-255.

[24]周波陽(yáng)，羅志才，林 旭，等.航空重力測(cè)量測(cè)線網(wǎng)平差中的粗差處理[J].大地測(cè)量與地球動(dòng)力學(xué)，2012，32（2）：110-114.

[25]周波陽(yáng)，羅志才，林 旭.航空重力測(cè)量測(cè)線網(wǎng)平差中的粗差處理方法[C]//中國(guó)地球物理學(xué)會(huì)第二十七屆年會(huì)論文集.合肥：中國(guó)科技大學(xué)出版社，2011.

[26]羅志才，寧津生，晁定波.衛(wèi)星重力梯度向下延拓的譜方法[J].測(cè)繪學(xué)報(bào)，1997，26（2）：168-175.

[27]Novák P，Kern M，Schwarz K P.Numerical studies on the harmonic downward continuation of band-limited airborne gravity [J].Studia Geophysica et Geodaetica，2001，45（4）：327-345.

[28]Bayoud F A，Sideris M G.Two different methodologies for geoid determination from ground and airborne gravity data[J]. Geophysical Journal International，2003，155（3）：914-922.

[29]Novák P，Kern M，Schwarz K P，et al.Evaluation of band-limited topographical effects in airborne gravimetry[J].Journal of Geodesy，2003，76（11-12）：597-604.

[30]蔣 濤.利用航空重力測(cè)量數(shù)據(jù)確定區(qū)域大地水準(zhǔn)面[D].武漢：武漢大學(xué)，2011.

[31]Serpas J.Local and regional geoid determination from airborne vector gravimetry[D].USA：The Ohio State University，2003.

[32]Hwang C，Hsiao Y S，Shih H C，et al.Geodetic and geophysical results from a Taiwan airborne gravity survey：Data reduction and accuracy assessment[J].Journal of Geophysical Research：Solid Earth（1978-2012），2007，112（B4）.

[33]周波陽(yáng)，羅志才，許 闖，等.航空重力測(cè)量數(shù)據(jù)向下延拓最小二乘配置法和逆Possion積分法的比較[C]//中國(guó)地球物理.合肥：中國(guó)科技大學(xué)出版社，2012.

[34]Novák P，Heck B.Downward continuation and geoid determination based on band-limited airborne gravity data[J].Journal of Geodesy，2002，76（5）：269-278.

[35]Wu L.Spectral methods for post processing of airborne vector gravity data[D].Canada：University of Calgary，1996.

[36]Schwarz K P，Li Z.An introduction to airborne gravimetry and its boundary value problems[M]//Geodetic boundary value problems in view of the one centimeter geoid.Germany：SpringerBerlin Heidelberg，1997：312-358.

[37]Novák P.Optional model for geoid determination from airborne gravity[J].Studia geophysica et geodaetica，2003，47：1-36.

[38]Schwarz K P，Sideris M G，F(xiàn)orsberg R.The use of FFT techniques in physical geodesy[J].Geophysical Journal International，1990，100（3）：485-514.

[39]Sideris M G，She B B.A new，high-resolution geoid for Canada and part of the US by the 1D-FFT method[J].Bulletin Géodésique，1995，69（2）：92-108.

[40]周波陽(yáng)，羅志才，林 旭，等.航空重力數(shù)據(jù)向下延拓FFT快速算法比較[J].大地測(cè)量與地球動(dòng)力學(xué)，2013，33（1）：64-68.

[41]Haines G V.Spherical cap harmonic analysis[J].Journal of Geophysical Research：Solid Earth（1978-2012），1985，90（B3）：2583-2591.

[42]Alldredge L R.Rectangular harmonic analysis applied to the geomagnetic field[J].Journal of Geophysical Research，1981，86（B4）：3021-3026.

[43]Li J C，Chao D B，Ning J S.Spherical cap harmonic expansion for local gravity field representation[J].Manuscripta Geodaetica，1995，20：265-277.

[44]邊少峰.大地測(cè)量邊值問題數(shù)值解法與地球重力場(chǎng)逼近[D].武漢：武漢測(cè)繪科技大學(xué)，1992.

[45]Jekeli C，Kwon J H.Geoid profile determination by direct integration of GNSS inertial navigation system vector gravimetry[J].Journal of Geophysical Research：Solid Earth（1978-2012），2002，107（B10）：ETG 3-1-ETG 3-10.

Research on the data processing methods of airborne vector gravimetry using SINS/GNSS

Ning Jinsheng1，2
（1.School of Geodesy and Geomatics，Wuhan University，Wuhan 430079，China；2.Key Laboratory of Geospace Environment and Geodesy，Ministry of Education，Wuhan University，Wuhan 430079，China）

Airborne vector gravimetry is an advanced and efficient technology to determine high frequency information of earth’s gravity field in the future，and is a complement to ground gravimetry，marine gravimetry and satellite gravimetry.The principle of airborne vector gravimetry using SINS（strapdown inertial navigation system）/GNSS（global navigation satellite systems）is introduced in the paper，and then the data preprocessing，data reduction，band-limited topographical effect and downward continuation are discussed，as well as the geoid determination from airborne vector gravimetry.It provides the support for the development of airborne vector gravimetry in our country.

airborne vector gravimetry；SINS/GNSS data fusion；earth gravity field；geoid

P227

A

1009-1742（2014）03-0004-10

2013-12-03

國(guó)家自然科學(xué)基金項(xiàng)目（41174062）；地球空間環(huán)境與大地測(cè)量教育部重點(diǎn)實(shí)驗(yàn)室開放基金項(xiàng)目（12-02-05，12-02-09）

寧津生，1932年出生，安徽桐城縣人，中國(guó)工程院院士，教授，主要從事物理大地測(cè)量的教學(xué)和科研工作；E-mail：jsning@sgg.whu.edu.cn