黃謨濤，寧津生，歐陽(yáng)永忠，，鄧凱亮，翟國(guó)君，陸秀平，吳太旗

1.海軍海洋測(cè)繪研究所，天津300061；2.海軍工程大學(xué)導(dǎo)航工程系，湖北 武漢430033；3.武漢大學(xué)測(cè)繪學(xué)院，湖北 武漢430079

1 引 言

自20世紀(jì)80年代末由于高精度衛(wèi)星定位系統(tǒng)的發(fā)展，使得航空重力測(cè)量技術(shù)取得實(shí)質(zhì)性突破以來，航空重力測(cè)量數(shù)據(jù)向下延拓問題一直受到極大的關(guān)注，國(guó)內(nèi)外學(xué)者為此開展了深入研究并提出了許多不同類型的解算方法，具有代表性的主要有梯度法、最小二乘配置法、FFT法、點(diǎn)質(zhì)量方法等［1－3］，但目前在實(shí)際應(yīng)用中，大多仍沿用傳統(tǒng)的球面Poisson積分方法解決此問題，即求逆Poisson積分方法，包括不同形式的迭代求解和非迭代求解法［4－5］。文獻(xiàn)［6］研究了航空重力測(cè)量數(shù)據(jù)向不規(guī)則面（包括地面）延拓的可能性，提出了顧及地形高度但無密度假設(shè)的向下延拓方法。文獻(xiàn)［7—9］先后提出直接代表法和球內(nèi)Dirichlet法。但由于重力場(chǎng)向下延拓計(jì)算在數(shù)學(xué)上屬于不適定反問題，其解算過程存在不穩(wěn)定性是該問題本身固有的一種屬性，即很小的觀測(cè)噪聲也可能引起模型解嚴(yán)重偏離真解［10－11］。因此，無論是早期的逆Poisson積分迭代求解法，還是后期的最小二乘配置法和FFT法，都無法確保向下延拓解算結(jié)果是絕對(duì)穩(wěn)定有效的。針對(duì)此問題，自20世紀(jì)末以來，國(guó)內(nèi)外學(xué)者的注意力一致投向?qū)⒄齽t化方法引入到航空重力測(cè)量向下延拓這個(gè)典型的不適定反問題解算中來，其研究成果可謂層出不窮［12－19］。但必須指出的是，航空重力測(cè)量向下延拓計(jì)算結(jié)果的有效性不僅取決于解算方法的穩(wěn)定性，同時(shí)取決于航空重力觀測(cè)數(shù)據(jù)的噪聲特性［20］。從已有研究成果看，在3km計(jì)算高度上，即使采用有針對(duì)性的正則化處理方法，向下延拓解算精度一般也只能達(dá)到±（3～5）mGal（mGal＝10－5m/s2）的水平［3，9，15，20］。而要想取得這樣的精度，還需要在觀測(cè)高度歸一化、觀測(cè)數(shù)據(jù)網(wǎng)格化、邊緣效應(yīng)處理和模型參數(shù)選擇等多個(gè)方面作出不懈的努力?？傊?，現(xiàn)有的航空重力測(cè)量向下延拓方法仍需要作進(jìn)一步的改進(jìn)和完善。為此，筆者曾從另一角度出發(fā)，提出了一種基于外部數(shù)據(jù)源、獨(dú)立于觀測(cè)數(shù)據(jù)的向下延拓方案［20］。針對(duì)海域重力場(chǎng)的變化特點(diǎn)，文獻(xiàn)［20］提出了直接利用超高階重力場(chǎng)模型計(jì)算延拓改正數(shù)的實(shí)用方法，其解算結(jié)果的穩(wěn)定性不再受到航空重力觀測(cè)噪聲的干擾，突破了向下延拓不適定問題的技術(shù)瓶頸，同時(shí)顯著簡(jiǎn)化了其計(jì)算過程和解算難度，取得了較好的解算效果［20］。本文繼續(xù)沿用這種“差分”思路將該方法拓展應(yīng)用到陸部，提出聯(lián)合使用超高階位模型和地形高數(shù)據(jù)，實(shí)施陸部航空重力向下延拓總改正數(shù)的計(jì)算，以飛行高度面與地面對(duì)應(yīng)點(diǎn)的位模型差分信息表征總改正數(shù)的中長(zhǎng)波分量，以相對(duì)應(yīng)的局部地形改正差分修正量表征總改正數(shù)的中高頻成分，從而實(shí)現(xiàn)航空重力數(shù)據(jù)點(diǎn)對(duì)點(diǎn)向地面的全頻段延拓。在地形變化不同區(qū)域，聯(lián)合使用EGM2008位模型、地面實(shí)測(cè)重力和高分辨率高程數(shù)據(jù)進(jìn)行了實(shí)際數(shù)值計(jì)算、精度評(píng)估和對(duì)比分析，證明該方法是可行有效的，具有推廣應(yīng)用價(jià)值。

2 計(jì)算模型與適用性分析

2.1 利用位模型和地形信息計(jì)算延拓改正數(shù)

2.1.1 問題的提出

如引言所述，針對(duì)海域重力場(chǎng)的變化特點(diǎn)，筆者已在文獻(xiàn)［20］中提出了直接利用超高階重力場(chǎng)模型計(jì)算海域延拓改正數(shù)的實(shí)用方法。由于EGM2008模型在建模時(shí)已經(jīng)充分利用了當(dāng)時(shí)最好的衛(wèi)星測(cè)高數(shù)據(jù)集［21］，因此在海域使用該模型計(jì)算航空重力測(cè)量向下延拓改正數(shù)，能夠取得較高精度的歸算結(jié)果［20］。但在地形變化比較復(fù)雜的陸部，EGM2008模型的適用性如何還不得而知。文獻(xiàn)［22］利用我國(guó)大陸GPS水準(zhǔn)高程異常和地面平均空間重力異常數(shù)據(jù)，對(duì)EGM2008模型在我國(guó)陸部的符合度進(jìn)行了評(píng)估，本文利用多個(gè)部門在我國(guó)及鄰近海域測(cè)量獲得的1500多萬(wàn)個(gè)船測(cè)重力數(shù)據(jù)，對(duì)EGM2008位模型進(jìn)行了外部檢核，現(xiàn)將兩部分比對(duì)結(jié)果匯總于表1。

表1 EGM2008模型與船測(cè)重力點(diǎn)值和陸部平均重力對(duì)比Tab.1 Differences of EGM2008gravity anomalies and the measurements mGal

從表1比對(duì)結(jié)果可以看出，EGM2008模型對(duì)我國(guó)海域重力場(chǎng)的逼近度較高，已接近于我國(guó)現(xiàn)役航空重力儀的測(cè)量精度水平（±（3～5）mGal）［23－24］，而在表征我國(guó)大陸地區(qū)5′×5′重力異常方面，盡管與EGM96相比，EGM2008模型的逼近度已經(jīng)有了較大幅度的提高，但比對(duì)誤差仍超過±10mGal。此結(jié)果說明，在我國(guó)陸部單純使用位模型計(jì)算延拓改正數(shù)，不會(huì)取得理想的結(jié)果，在地形變化比較劇烈的山區(qū)，情況更是如此。為了改善位模型延拓改正數(shù)的計(jì)算精度，必須考慮在原有基礎(chǔ)上增加陸部重力場(chǎng)高頻信息的作用。地形是一種可見的重力場(chǎng)高頻信息，在局部范圍內(nèi)，空間重力異常與地形的相關(guān)度可達(dá)0.95以上［8］。因此，聯(lián)合利用位模型和地形信息確定延拓改正數(shù)，成為本文研究思路的基本出發(fā)點(diǎn)。

2.1.2 基于差分局部地形改正的延拓歸算模型

由高等物理大地測(cè)量學(xué)得知［5］，Molodensky邊值問題解是通過計(jì)算G1、G2、G3、…各項(xiàng)改正來顧及地球表面形態(tài)變化的。在Bjerhammar邊值理論中，地形質(zhì)量影響則是通過解析延拓手段間接轉(zhuǎn)移到等效重力異?；蛱摂M點(diǎn)質(zhì)量上來的。文獻(xiàn)［4—5］將G1分解為兩項(xiàng)G1＝G11＋G12，同時(shí)證明了G11正好是將地面重力異?！把油亍钡胶Ｆ矫娴目臻g改正，而G12的Stokes積分又相當(dāng)于從海平面“回到”地面的高程異常改正［25］。當(dāng)重力異常與地形高度呈線性相關(guān)時(shí)，可證明局部地形改正C和G1項(xiàng)改正是近似等價(jià)的。以上情況說明，局部地形改正是表征地球重力場(chǎng)高頻信息的一類關(guān)鍵參量，在局部重力場(chǎng)逼近計(jì)算中具有非常重要的作用。重力位模型由于受空間分辨率的限制，高頻含量明顯不足，因此很難在表示地球重力場(chǎng)局部特征方面有所作為，局部地形改正參量正好能夠彌補(bǔ)這方面的缺陷?？罩幸稽c(diǎn)重力異常的高頻分量隨高度增大而衰減，空中測(cè)點(diǎn)與對(duì)應(yīng)地面投影點(diǎn)局部地形改正數(shù)之差異，正是這種高頻信號(hào)衰減的量化體現(xiàn)。據(jù)此，本文提出陸部航空重力測(cè)量向下延拓計(jì)算模型，見式（1）和式（2）

局部地形改正可按照文獻(xiàn)［26］提出的球面公式進(jìn)行嚴(yán)密計(jì)算，但考慮到敘述方便起見，這里仍采用平面近似方式給出相應(yīng)的計(jì)算模型。計(jì)算O點(diǎn)處的局部地形改正數(shù)C0時(shí)，不管是高出計(jì)算點(diǎn)部分（h－h(huán)0＞0），還是低于計(jì)算點(diǎn)部分（h－h(huán)0＜0），其地形改正數(shù)均為正值。C0的統(tǒng)一計(jì)算式為［25］

或?qū)憺?/p>

式中，R為地球平均半徑；hx和hy分別代表計(jì)算點(diǎn)在經(jīng)圈和緯圈方向上的地形梯度。

計(jì)算空中測(cè)點(diǎn)P的局部地形改正數(shù)Cp時(shí)，仍然選擇與其相對(duì)應(yīng)的地面投影點(diǎn)O作為地形高程面，此時(shí)高出O點(diǎn)高程面部分的地形改正數(shù)為負(fù)值，低于O點(diǎn)高程面部分的地形改正數(shù)為正值。Cp的統(tǒng)一計(jì)算式為

在實(shí)際應(yīng)用中，一般使用數(shù)值積分完成式（3）和式（6）的計(jì)算，式（4）更適合使用FFT技術(shù)進(jìn)行解算，但當(dāng)使用特別高分辨率的地形數(shù)據(jù)（如3″×3″甚至1″×1″）時(shí)，由于不能確保，式（3）無法近似過渡到式（4），此時(shí)需對(duì)核函數(shù)作非奇變換［27］，因此必須注意，使用式（4）是有前提條件的?？罩袦y(cè)點(diǎn)和地面點(diǎn)局部地形改正計(jì)算各參量之間的幾何關(guān)系見圖1。

圖1 局部地形改正參量幾何關(guān)系Fig.1 Geometry relations of parameters in topographic correction computations

不難看出，本文提出的延拓計(jì)算模型即式（1），同樣適用于陸海交界區(qū)域的重力向下延拓解算。當(dāng)計(jì)算點(diǎn)從陸區(qū)向海域延伸時(shí)，局部地形差分改正數(shù)δCp0將逐步減??；當(dāng)計(jì)算點(diǎn)距離陸地地形超出局部地形改正積分半徑時(shí)，差分改正數(shù)δCp0為零。

2.2 基于差分層間地形改正的延拓歸算模型

文獻(xiàn)［7］曾依據(jù)空中一點(diǎn)重力異常主要與其相應(yīng)地面一定區(qū)域內(nèi)的重力異常大小相關(guān)聯(lián)這一事實(shí)，提出直接使用空中測(cè)點(diǎn)重力異常代替相應(yīng)地面方塊平均重力異常的延拓歸算方案，簡(jiǎn)稱直接代表法。文獻(xiàn)［8］在此基礎(chǔ)上，提出了利用地形信息改善直接代表法的改進(jìn)型延拓歸算方案，其基本思想是：首先通過使用空中測(cè)點(diǎn)重力異常代表相應(yīng)地面方塊平均重力異常的方法，確定地面投影點(diǎn)重力異常的主項(xiàng)，然后利用地形信息補(bǔ)充該主項(xiàng)的高頻分量，最終按式（7）完成延拓歸算［8］

式中，HT代表地面點(diǎn)高程；代表最佳相關(guān)區(qū)平均高程，可利用高分辨率DEM數(shù)據(jù)按逐次趨近方法計(jì)算得到；b為布格系數(shù)，取b＝2πGρ0＝0.111 9mGal/m，其他符號(hào)意義同前。

文獻(xiàn)［8］的數(shù)值試驗(yàn)結(jié)果表明，利用式（7）作為延拓歸算模型的主要缺陷是，延拓計(jì)算值與地面觀測(cè)值之間存在一個(gè)比較明顯的系統(tǒng)性偏差。這顯然跟直接使用空中測(cè)點(diǎn)重力異常作為地面點(diǎn)重力異常的主項(xiàng)存在不確定性有關(guān)。實(shí)際上，隨著飛行高度的增大，空中測(cè)點(diǎn)重力異常不僅高頻分量在衰減，中低頻分量也在衰減。因此，在和之間，除了由地形質(zhì)量引起的高頻分量差異外，還應(yīng)包括由高度差引起的中低頻分量差異。本文在式（1）中特別引入位模型延拓改正數(shù)項(xiàng)，其作用就是為了補(bǔ)償由高度差引起的中低頻分量衰減。為了對(duì)比分析不同方法的計(jì)算效果，這里特意列出現(xiàn)行國(guó)家軍用標(biāo)準(zhǔn)選用的直接代表法模型即式（7）［28］，并在后面的數(shù)值檢核中一同給出其計(jì)算結(jié)果。因式（7）右端第二項(xiàng)實(shí)為地面投影點(diǎn)與平均高程點(diǎn)的層間改正互差，故將式（7）稱為基于差分層間地形改正的延拓歸算模型。

2.3 新方法特點(diǎn)分析

本文方法的最大特點(diǎn)是，向下延拓改正數(shù)計(jì)算過程完全獨(dú)立于航空重力測(cè)量觀測(cè)數(shù)據(jù)，不受觀測(cè)噪聲的影響，更重要的是，不受向下延拓不適定反問題解算過程固有的不穩(wěn)定性影響。如前所述，由于空中一點(diǎn)重力異常的高頻分量隨高度增大而衰減，因此相比較而言，空中重力場(chǎng)變化要比地面重力場(chǎng)平緩是不言而喻的。而從理論上講，要想通過求解傳統(tǒng)的逆Poisson積分方程，單純由變化相對(duì)平緩的空中重力異常信息恢復(fù)包含更多高頻分量的地面重力異常，幾乎是不可能的。即使能夠獲取一部分高頻信息，那也是觀測(cè)噪聲被放大后的虛假的重力異常信息。因此從這個(gè)角度上講，傳統(tǒng)延拓方法一直存在理論上的缺陷。本文方法的獨(dú)特之處正是完全避開了傳統(tǒng)方法的弊端，首先利用超高階地球位模型恢復(fù)延拓改正數(shù)的中長(zhǎng)波部分，然后利用地形信息恢復(fù)地面重力場(chǎng)的高頻分量，從而實(shí)現(xiàn)航空重力測(cè)量數(shù)據(jù)向地面的全頻延拓。

3 精度分析與數(shù)值檢核

3.1 精度分析

由式（1）知，本文提出的陸部航空重力測(cè)量向下延拓計(jì)算模型精度，主要取決于位模型延拓改正數(shù)和局部地形差分改正數(shù)δCp0的計(jì)算精度。筆者曾對(duì)位模型延拓改正數(shù)的計(jì)算精度做過理論估算［20］，其中，在5km延拓高度上，由位模型系數(shù)誤差引起的延拓改正數(shù)計(jì)算誤差不超過2.5mGal，由位模型階次截?cái)嘁鸬难油馗恼龜?shù)計(jì)算誤差不超過2.8mGal。因在式（1）中已經(jīng)引入局部地形差分改正數(shù)，來補(bǔ)償位模型缺失的高頻分量，故在新模型中不需考慮后一項(xiàng)位模型截?cái)嗾`差的影響。又由文獻(xiàn)［25—27］分析結(jié)果知，在具備高精度和高分辨率數(shù)字地形模型條件下，局部地形差分改正數(shù)的計(jì)算誤差完全可以控制在1mGal以內(nèi)。因此，從理論上講，新模型中兩項(xiàng)延拓改正數(shù)的總體估算精度應(yīng)優(yōu)于3mGal。但必須指出的是，新模型的實(shí)際計(jì)算精度還受多種不確定因素的影響，包括位模型在陸部逼近度的非均勻性及數(shù)字地形模型精度和分辨率水平的不一致性等。

3.2 數(shù)值檢核

3.2.1 檢核方法與試驗(yàn)數(shù)據(jù)

檢驗(yàn)向下延拓計(jì)算模型精度最有效的方法是，直接使用航空和地面實(shí)際重力觀測(cè)數(shù)據(jù)求互差作為基準(zhǔn)值，將由新模型計(jì)算得到的延拓改正數(shù)與基準(zhǔn)值進(jìn)行比較，由此可得到新模型計(jì)算精度的評(píng)價(jià)指標(biāo)。遺憾的是，陸部同時(shí)擁有高精度航空和地面實(shí)際重力觀測(cè)數(shù)據(jù)的區(qū)域并不多見，目前筆者還缺乏這方面的可靠資料。為此，本文改用向上與向下延拓比對(duì)方法對(duì)新模型計(jì)算精度進(jìn)行外部檢核，即首先利用地面網(wǎng)格重力和地形高數(shù)據(jù)，通過向上延拓方法計(jì)算一定高度面上的空間重力異常，進(jìn)而求取地面與計(jì)算高度面重力異常的互差，將此互差值分別與單獨(dú)由位模型計(jì)算得到的改正數(shù)（即文獻(xiàn)［20］中的式（8））及由式（1）、式（7）計(jì)算得到的延拓改正數(shù)進(jìn)行逐一比較，即可獲得相應(yīng)延拓歸算模型的精度評(píng)價(jià)。由于向上延拓過程是穩(wěn)定可靠的，因此其解算結(jié)果可作為檢核向下延拓計(jì)算結(jié)果的基準(zhǔn)值。這里考慮到向上延拓的Poisson積分方程是球面公式，當(dāng)?shù)匦纹鸱^大時(shí)，利用球面向上延拓公式推算低空重力異常精度難以得到保證［4］。為此，本文改用虛擬點(diǎn)質(zhì)量方法進(jìn)行地面重力異常向上延拓計(jì)算，具體解算模型可參見文獻(xiàn)［25］和［27］。因點(diǎn)質(zhì)量模型以最簡(jiǎn)單的方式顧及了地形效應(yīng)，故向上延拓解算結(jié)果更接近實(shí)際［27］。

本文選用的試驗(yàn)區(qū)包括兩個(gè)區(qū)塊，一個(gè)位于我國(guó)大陸中部，區(qū)塊大小為3°×3°（稱為區(qū)塊1），屬于地形變化相對(duì)比較平緩的中等山區(qū)；另一個(gè)位于美國(guó)大陸南部，區(qū)塊大小為4°×4°（稱為區(qū)塊2），屬于地形變化比較劇烈的大山區(qū)。區(qū)塊1和2分別擁有1′×1′和2′×2′網(wǎng)格地面重力異常數(shù)據(jù)，同時(shí)擁有30″×30″地形高數(shù)據(jù)。兩個(gè)區(qū)塊所對(duì)應(yīng)的重力和地形數(shù)據(jù)變化特征如表2統(tǒng)計(jì)結(jié)果所示。

表2 兩個(gè)試驗(yàn)區(qū)塊重力和地形數(shù)據(jù)變化特征Tab.2 Statistics of gravity anomalies and terrain data in the test areas

3.2.2 區(qū)塊1檢核結(jié)果

試驗(yàn)區(qū)塊1的平均地形高為358m，最大高程為1698m。據(jù)此，這里取飛行高度（相對(duì)橢球面高）hp＝3km，進(jìn)行地面重力異常Δg地向上延拓計(jì)算，首先解算得到3km高度面上的空間重力異常Δgh，進(jìn)而計(jì)算相對(duì)應(yīng)的延拓改正數(shù)基準(zhǔn)值δΔg地＝Δg地－Δgh。又利用EGM2008位模型和地面30″×30″地形高數(shù)據(jù)，分別按文獻(xiàn)［20］中的式（8）計(jì)算位模型延拓改正數(shù)δΔg位，按式（1）計(jì)算位模型和差分局部地形改正的延拓總改正數(shù)δΔg位局＝δΔg位－δC，按式（7）計(jì)算差分層間地形改正的延拓改正數(shù)。將基準(zhǔn)值δΔg地分別同以上3組計(jì)算值δΔg位、δΔg位局和δΔg層作比較，可得到相應(yīng)的精度評(píng)估信息。表3為地面和3km高度面局部地形改正計(jì)算值統(tǒng)計(jì)結(jié)果，具體比對(duì)結(jié)果見表4。

表3 區(qū)塊1地面點(diǎn)與3km高度點(diǎn)局部地形改正計(jì)算值統(tǒng)計(jì)結(jié)果Tab.3 Statistics of local topographic corrections at ground and 3km altitude in area 1st mGal

表4 第1區(qū)塊3種延拓方案計(jì)算結(jié)果與地面實(shí)測(cè)數(shù)據(jù)計(jì)算延拓改正數(shù)的比較Tab.4 Comparisons between the continuation corrections from ground gravity and three models in area 1st mGal

為了說明單獨(dú)使用位模型計(jì)算延拓改正數(shù)的精度并不完全取決于位模型自身的絕對(duì)精度［20］，表4第2行同時(shí)給出了同一區(qū)域位模型重力異常與地面實(shí)測(cè)網(wǎng)格重力的比對(duì)結(jié)果。從表4可以看出，位模型重力異常與地面觀測(cè)重力比對(duì)的系統(tǒng)性偏差超過1mGal、均方根誤差超過9mGal，但由位模型計(jì)算延拓改正數(shù)的精度仍達(dá)到了6mGal，加上局部地形改正差分延拓改正數(shù)后的比對(duì)精度則提升到了3mGal，系統(tǒng)性偏差減小到0.5mGal。而利用層間地形改正模型進(jìn)行延拓計(jì)算的比對(duì)效果并不十分理想，均方根誤差接近6mGal。這個(gè)結(jié)果說明，將本文提出的差分模式應(yīng)用于陸部航空重力測(cè)量數(shù)據(jù)延拓處理，能夠取得較高精度的延拓解算結(jié)果。

3.2.3 區(qū)塊2檢核結(jié)果

為了進(jìn)一步了解各種計(jì)算方案在地形變化劇烈山區(qū)的適用性，這里繼續(xù)在表2所示的區(qū)塊2進(jìn)行與表4相類同的數(shù)值檢核。其中，考慮到區(qū)塊2的平均地形高為1987m，最大高程為4116m，實(shí)施地面重力異常向上延拓計(jì)算時(shí)，特別把飛行平均高度改取為hp＝5km，其他參量保持不變。表5為地面和5km高度面局部地形改正計(jì)算值統(tǒng)計(jì)結(jié)果，具體比對(duì)結(jié)果見表6。圖2同時(shí)給出了（δΔg地－δΔg位）和（δΔg地－δΔg位局）的等值線圖，從圖2可清晰看出，單獨(dú)由位模型計(jì)算得到的延拓改正數(shù)誤差（見圖2（a））明顯與地形變化幅度大小有關(guān)；由位模型和局部地形改正差分得到的延拓改正數(shù)與地面實(shí)測(cè)數(shù)據(jù)解算結(jié)果則具有很好的一致性（見圖2（b）），誤差分布比較均勻，不存在明顯的邊緣效應(yīng)。

從表6結(jié)果可以看出，即使在地形變化劇烈的大山區(qū)，由位模型計(jì)算延拓改正數(shù)的精度仍達(dá)到了5.71mGal，加上局部地形改正差分延拓改正數(shù)后的比對(duì)精度則提升到了2.97mGal，系統(tǒng)性偏差很小。而利用層間地形改正模型計(jì)算延拓改正數(shù)的效果明顯失效。對(duì)比表4和表6計(jì)算結(jié)果不難看出，地形變化劇烈大山區(qū)（即區(qū)塊2）的比對(duì)效果還略好于地形變化相對(duì)平緩的區(qū)塊1，這個(gè)結(jié)果也從另一個(gè)側(cè)面說明，EGM2008模型在美國(guó)本土的符合度要明顯高于在我國(guó)的符合度，與該模型在建模時(shí)已經(jīng)使用了美國(guó)本土比較高精度和高分辨率重力及地形數(shù)據(jù)這一事實(shí)相符。

表5 區(qū)塊2地面點(diǎn)與5km高度點(diǎn)局部地形改正計(jì)算值統(tǒng)計(jì)結(jié)果Tab.5 Statistics of local topographic corrections at ground and 5km altitude in area 2nd mGal

表6 第2區(qū)塊3種延拓方案計(jì)算結(jié)果與地面實(shí)測(cè)數(shù)據(jù)計(jì)算延拓改正數(shù)的比較Tab.6 Comparisons between the continuation corrections from ground gravity and three models in area 2nd mGal

圖2 地面數(shù)據(jù)延拓改正數(shù)與模型計(jì)算值互差等值線圖Fig.2 Maps of continuation correction differences between the true and model values

3.2.4 正則化方法檢核結(jié)果

為了從另一個(gè)側(cè)面說明本文新方法的有效性，這里進(jìn)一步采用文獻(xiàn)［19］介紹的逆Poisson積分正則化方法，將由點(diǎn)質(zhì)量模型向上延拓計(jì)算得到的3km高度（對(duì)應(yīng)區(qū)塊1）和5km高度（對(duì)應(yīng)區(qū)塊2）上的重力異常，分別向下延拓到兩個(gè)區(qū)塊的平均地形高度，然后將向下延拓計(jì)算結(jié)果與地面已知重力異常作比較，可得到向下延拓計(jì)算模型的精度評(píng)價(jià)指標(biāo)，其結(jié)果見表7。

表7 兩區(qū)塊逆Poisson積分正則化方法向下延拓比對(duì)統(tǒng)計(jì)結(jié)果Tab.7 Comparisons of the continuation corrections from ground gravity and regularization model in test areas mGal

由于這里采用的空中重力異常是由點(diǎn)質(zhì)量模型向上延拓計(jì)算得到的，因此它與地面已知重力異常之間形成閉環(huán)關(guān)系，即相當(dāng)于不存在通常意義下的觀測(cè)誤差，向下延拓計(jì)算結(jié)果與已知值互差所反映的應(yīng)當(dāng)是純粹的計(jì)算模型誤差。對(duì)比表7和表4、表6計(jì)算結(jié)果可以看出，即使沒有觀測(cè)噪聲的干擾，逆Poisson積分正則化方法的延拓計(jì)算比對(duì)效果，無論是系統(tǒng)偏差還是均方根誤差指標(biāo)，都不及本文提出的新方法。此外，從計(jì)算過程看，新方法不需要對(duì)觀測(cè)重力做網(wǎng)格化處理，能夠?qū)崿F(xiàn)點(diǎn)對(duì)點(diǎn)延拓計(jì)算也是傳統(tǒng)逆積分方法無法比擬的。

4 結(jié)論與建議

前面的理論分析和數(shù)值檢核結(jié)果表明，利用EGM2008模型和地形信息計(jì)算延拓總改正數(shù)，將陸部航空重力數(shù)據(jù)向地面延拓是可行有效的，具有很好的推廣應(yīng)用價(jià)值。

（1）聯(lián)合使用超高階位模型和地形信息進(jìn)行陸部航空重力數(shù)據(jù)向下延拓解算，在5km高度下，其理論估計(jì)精度優(yōu)于±3mGal；在地形變化劇烈的大山區(qū)，實(shí)際比對(duì)精度達(dá)±3mGal。此結(jié)果與當(dāng)前國(guó)際航空重力測(cè)量精度水平相當(dāng)，因此使用新方法不會(huì)顯著降低航空重力測(cè)量成果質(zhì)量。

（2）利用本文方法計(jì)算向下延拓改正數(shù)，不僅簡(jiǎn)化了向下延拓計(jì)算過程，更重要的是可獲得穩(wěn)定可靠的解算結(jié)果。同時(shí)相對(duì)降低了航空重力測(cè)量測(cè)線布設(shè)要求，不需要對(duì)觀測(cè)數(shù)據(jù)作高度歸一化和網(wǎng)格化處理，可直接對(duì)不同高度的測(cè)點(diǎn)進(jìn)行點(diǎn)對(duì)點(diǎn)計(jì)算。

（3）本文方法計(jì)算效果優(yōu)于傳統(tǒng)的直接代表法和逆積分正則化方法，建議盡早對(duì)現(xiàn)行的國(guó)家軍用標(biāo)準(zhǔn)相關(guān)內(nèi)容作必要的修改。

［1］MUELLER F，MAYER－GUERR T.Comparison of Downward Continuation Methods of Airborne Gravimetry Data［C］∥International Association of Geodesy Symposia Volume 128：A Window on the Future of Geodesy.Berlin：Springer，2003：254－258.

［2］TZIAVOS I N，ANDRITSANOS V D，F(xiàn)ORSBERG R，et al.Numerical Investigation of Downward Continuation Methods for Airborne Gravity Data［C］∥International Association of Geodesy Symposia Volume 129：Gravity，Geoid and Space Missions.Berlin：Springer，2005：119－124.

［3］WANG Xingtao，XIA Zheren，SHI Pan，et al.A Comparison of Different Downward Continuation Methods for Airborne Gravity Data［J］.Chinese Journal of Geophysics，2004，47（6）：1017－1022.（王興濤，夏哲仁，石磐，等.航空重力測(cè)量數(shù)據(jù)向下延拓方法比較［J］.地球物理學(xué)報(bào)，2004，47（6）：1017－1022.）

［4］HEISKANEN W A，MORITZ H.Physical Geodesy［M］.San Francisco：Freeman W H and Company，1967.

［5］MORITZ H.Advanced Physical Geodesy［M］.Karlsruhe：Herbert Wichmann Verlag，1980.

［6］JEKELI C.The Downward Continuation of Aerial Gravimetric Data without Density Hypothesis［J］.Bulletin Géodésique，1987，61（4）：319－329.

［7］SHI Pan，WANG Xingtao.The Frequence Domain Analysis for the Determination of Terrestrial Mean Gravity Anomaly Using Airborne Gravimetry［J］.Acta Geodaetica et Cartographica Sinica，1995，24（4）：301－308.（石磐，王興濤.空中測(cè)量地面平均重力異常的頻域分析［J］.測(cè)繪學(xué)報(bào)，1995，24（4）：301－308.）

［8］SHI Pan，WANG Xingtao.Determination of the Terrain Surface Gravity Field Using Airborne Gravimetry and DEM［J］.Acta Geodaetica et Cartographica Sinica，1997，26（2）：117－121.（石磐，王興濤.利用航空重力測(cè)量和DEM確定地面重力場(chǎng)［J］.測(cè)繪學(xué)報(bào)，1997，26（2）：117－121.）

［9］SHI Pan，SUN Zhongmiao.The Solution to the Problem of the Spherical Interior Dirichlet and Its Application［J］.Acta Geodaetica et Cartographica Sinica，1999，28（3）：195－198.（石磐，孫中苗.球內(nèi)Dirichlet問題解及其應(yīng)用［J］.測(cè)繪學(xué)報(bào)，1999，28（3）：195－198.）

［10］KELLER W，HIRSCH M.Downward Continuation versus Free－air Reduction in Airborne Gravimetry［C］∥International Association of Geodesy Symposium Volume 112：Geodesy and Physics of the Earth.Heidelberg，Berlin：Springer Verlag，1992，266－270.

［11］WANG Yanfei.Computational Methods for Inverse Problems and Their Applications［M］.Beijing：Higher Education Press，2007.（王彥飛.反演問題的計(jì)算方法及其應(yīng)用［M］.北京：高等教育出版社，2007.）

［12］KERN M.An Analysis of the Combination and Downward Continuation of Satellite，Airborne and Terrestrial Gravity Data［D］.Calgary：University of Calgary，2003.

［13］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，2007，112（B4），DOI：10.1029/2005JB004220.

［14］ALBERTS B.Regional Gravity Field Modeling Using Airborne Gravimetry Data［C］∥Publications on Geodesy 70.Delft：Netherlands Geodetic Commission，2009.

［15］WANG Xingtao，SHI Pan，ZHU Feizhou.Regularization Methods and Spectral Decomposition for the Downward Continuation of Airborne Gravity Data［J］.Acta Geodaetica et Cartographica Sinica，2004，33（1）：33－38.（王興濤，石磐，朱非洲.航空重力測(cè)量數(shù)據(jù)向下延拓的正則化算法及其譜分解［J］.測(cè)繪學(xué)報(bào)，2004，33（1）：33－38.）

［16］GU Yongwei，GUI Qingming.Study of Regularization Based on Signal－to－noise Index in Airborne Gravity Downward to the Earth Surface［J］.Acta Geodaetica et Cartographica Sinica，2010，39（5）：458－464.（顧勇為，歸慶明.航空重力測(cè)量數(shù)據(jù)向下延拓基于信噪比的正則化方法的研究［J］.測(cè)繪學(xué)報(bào)，2010，39（5）：458－464.）

［17］JIANG Tao，LI Jiancheng，WANG Zhengtao，et al.Solution of Ill－posed Problem in Downward Continuation of Airborne Gravity［J］.Acta Geodaetica et Cartographica Sinica，2011，40（6）：684－689.（蔣濤，李建成，王正濤，等.航空重力向下延拓病態(tài)問題的求解［J］.測(cè)繪學(xué)報(bào)，2011，40（6）：684－689.）

［18］DENG Kailiang，HUANG Motao，BAO Jingyang，et al.Tikhonov Two－parameter Regulation Algorithm in Downward Continuation of Airborne Gravity Data［J］.Acta Geodaetica et Cartographica Sinica，2011，40（6）：690－696.（鄧凱亮，黃謨濤，暴景陽(yáng)，等.向下延拓航空重力數(shù)據(jù)的Tikhonov雙參數(shù)正則化法［J］.測(cè)繪學(xué)報(bào)，2011，40（6）：690－696.）

［19］WU Taiqi，DENG Kailiang，HUANG Motao，et al.An Improved Singular Values Decomposition Method for Illposed Problem［J］.Geomatics and Information Science of Wuhan University，2011，36（8）：900－903.（吳太旗，鄧凱亮，黃謨濤，等.一種改進(jìn)的不適定問題奇異值分解法［J］.武漢大學(xué)學(xué)報(bào)：信息科學(xué)版，2011，36（8）：900－903.）

［20］HUANG Motao，OUYANG Yongzhong，LIU Min，et al.Practical Methods for the Downward Continuation of Airborne Gravity Data in the Sea Area［J］.Geomatics and Information Science of Wuhan University，2014，39（10）：1147－1152.（黃謨濤，歐陽(yáng)永忠，劉敏，等.海域航空重力測(cè)量數(shù)據(jù)向下延拓的實(shí)用方法［J］.武漢大學(xué)學(xué)報(bào)：信息科學(xué)版，2014，39（10）：1147－1152.）

［21］PAVLIS N K，HOLMES S A，KENYON S C，et al.The Development and Evaluation of the Earth Gravitational Model 2008（EGM2008）［J］.Journal of Geophysical Research，2012，117（B4），DOI： 10.1029/2011JB008916.

［22］ZHANG Chuanyin，GUO Chunxi，CHEN Junyong，et al.EGM 2008and Its Application Analysis in Chinese Mainland［J］.Acta Geodaetica et Cartographica Sinica，2009，38（4）：283－289.（章傳銀，郭春喜，陳俊勇，等.EGM2008地球重力場(chǎng)模型在中國(guó)大陸適用性分析［J］.測(cè)繪學(xué)報(bào)，2009，38（4）：283－289.）

［23］SUN Zhongmiao，ZHAI Zhenhe，LI Yingchun，et al.Concurrent Flight Test of LaCoste ＆ Romberg（LCR）Airborne GravimeterⅡ andⅠ System［J］.Journal of Geodesy and Geodynamics，2012，32（2）：24－27.（孫中苗，翟振和，李迎春，等.LCRⅡ型和Ⅰ型航空重力儀的同機(jī)飛行試驗(yàn)［J］.大地測(cè)量與地球動(dòng)力學(xué)，2012，32（2）：24－27.）

［24］OUYANG Yongzhong，DENG Kailiang，LU Xiuping，et al.Tests of Multi－type Airborne Gravimeters and Data Analysis［J］.Hydrographic Surveying and Charting，2013，33（4）：6－11.（歐陽(yáng)永忠，鄧凱亮，陸秀平，等.多型航空重力儀同機(jī)測(cè)試及其數(shù)據(jù)分析［J］.海洋測(cè)繪，2013，33（4）：6－11.）

［25］LI Jiancheng，CHEN Junyong，NING Jinsheng，et al.Theory of the Earth’s Gravity Field Approximation and Determination of China Quasi－geoid 2000［M］.Wuhan：Wuhan University Press，2003.（李建成，陳俊勇，寧津生，等.地球重力場(chǎng)逼近理論與中國(guó)2000似大地水準(zhǔn)面的確定［M］.武漢：武漢大學(xué)出版社，2003.）

［26］ZHANG Chuanyin，CHAO Dingbo，DING Jian，et al.Precision Topographical Effects for any Kond of Field Quantities for any Altitude［J］.Acta Geodaetica et Cartographica Sinica，2009，38（1）：28－34.（章傳銀，晁定波，丁劍，等.球近似下地球外空間任意類型場(chǎng)元的地形影響［J］.測(cè)繪學(xué)報(bào)，2009，38（1）：28－34.）

［27］HUANG Motao，ZHAI Guojun，GUAN Zheng，et al.The Determination and Application of Marine Gravity Field［M］.Beijing：Surveying and Mapping Press，2005.（黃謨濤，翟國(guó)君，管錚，等.海洋重力場(chǎng)測(cè)定及其應(yīng)用［M］.北京：測(cè)繪出版社，2005.）

［28］HEADQUARTERS OF GENERAL EQUIPMENT OF PLA.GJB 6561—2008，Rules for Operations of Airborne Gravimetry［S］.Beijing：Military Standard Press of the Headquarters of General Equipment，2008.（中國(guó)人民解放軍總裝備部.GJB 6561—2008，航空重力測(cè)量作業(yè)規(guī)范［S］.北京：總裝備部，2008.）