蔣 濤，黨亞民，章傳銀，劉站科

1.中國測繪科學(xué)研究院，北京 100830；2.國家測繪地理信息局 第一大地測量隊，陜西 西安 710054

1 引 言

經(jīng)濟(jì)高效的航空重力測量能夠很好地填補(bǔ)衛(wèi)星重力和地面重力測量所得重力場信息的信號盲區(qū)，可覆蓋地面重力測量難以達(dá)到的地區(qū)，數(shù)據(jù)精度在全波長為5～10km的空間分辨率尺度上能夠達(dá)到1～3mGal（1Gal＝1cm／s2）。航空重力測量獲取的是距地面達(dá)數(shù)千米飛行高度處的重力觀測值，大地測量和地球物理等領(lǐng)域需要的是地面或大地水準(zhǔn)面上的重力值，因此需將飛行高度處的重力觀測值向下調(diào)和延拓至地面或大地水準(zhǔn)面上。

在航空重力觀測值向下延拓過程中噪聲會得到放大，可能造成重力場信號的嚴(yán)重失真，若不采用合適的向下延拓方法，將無法得到穩(wěn)定可靠的重力延拓解?，F(xiàn)有的向下延拓方法主要包括逆泊松積分法［1－7］、最小二乘配置（LSC）［8］、快速傅里葉變換（FFT）［9-10］、直接代表法［11］和直接法［12－13］，其中比較常用的是逆泊松積分法，求解泊松積分方程時一般需引入正則化處理。

本文提出了一種基于矩諧分析的航空重力向下延拓新方法，矩諧分析最開始被文獻(xiàn)［14—15］用于區(qū)域地磁場建模，在地磁學(xué)研究中應(yīng)用廣泛。文獻(xiàn)［16］將矩諧分析用于局部重力場逼近，數(shù)值結(jié)果驗證了其可行性。該方法的基本原理是在局部直角坐標(biāo)系中將擾動位或其泛函展開成矩諧級數(shù)，利用飛行高度處的重力異常或重力擾動進(jìn)行矩諧分析以求解表征該區(qū)域內(nèi)重力場的擾動位系數(shù)，再進(jìn)行矩諧綜合計算得到地面或大地水準(zhǔn)面上的重力異常和重力擾動，從而完成重力觀測值從飛行高度處到地面或大地水準(zhǔn)面的向下調(diào)和延拓。本文利用EGM2008重力位模型［17］設(shè)計模擬數(shù)值試驗，將矩諧分析與直接法和基于廣義嶺估計的逆泊松積分法作數(shù)值比較分析，驗證了基于矩諧分析的向下延拓方法的可靠性、精度和穩(wěn)定性。

2 數(shù)學(xué)模型

矩諧分析表達(dá)區(qū)域地球重力場的基本思想是在局部直角坐標(biāo)系下求解拉普拉斯方程，將引力位展開成矩諧級數(shù)，以矩諧系數(shù)的集合來表征引力位。矩諧分析是在局部直角坐標(biāo)系中進(jìn)行的，因此需選取數(shù)據(jù)區(qū)域的中心點(diǎn)作為局部直角坐標(biāo)系的原點(diǎn)，將輸入數(shù)據(jù)點(diǎn)的大地坐標(biāo)或球坐標(biāo)轉(zhuǎn)換到局部直角坐標(biāo)系中。圖1以航空重力測量數(shù)據(jù)為例，給出了局部直角坐標(biāo)系下區(qū)域重力場表達(dá)的幾何關(guān)系。

圖1 局部直角坐標(biāo)系下區(qū)域重力場表達(dá)的幾何關(guān)系Fig.1 Geometry of the regional gravity field representation in local Cartesian coordinate system

在航空重力測量中，GPS動態(tài)定位能夠精密確定飛行高度處的大地高，重力擾動的獲取已不成問題，可采用重力擾動作為重力觀測值。在此直接給出重力擾動的矩諧展開式，推導(dǎo)過程請參見文獻(xiàn)［18］

式中

式（1）中的τnm定義為

式中，Dx、Dy分別表示矩形計算區(qū)域在東西（x）和南北（y）方向上的距離（見圖1）。矩諧展開式（1）對應(yīng)的半波長空間分辨率為

式中，y為k階空中重力擾動觀測值向量；A為k×u階設(shè)計矩陣；x為u階待求擾動位系數(shù)向量；e為觀測噪聲，其數(shù)學(xué)期望為0；D｛y｝ 為觀測值的方差—協(xié)方差陣；σ2是先驗方差因子；P是權(quán)陣，對應(yīng)的法方程為

其最小二乘解為

矩諧展開模型（1）采用周期性函數(shù)傅里葉級數(shù)表達(dá)重力場信號，其成立的前提是假定待求信號在計算區(qū)域內(nèi)也是周期性的。顯然有限區(qū)域的重力場信號不具有完全的周期性，用周期函數(shù)表示非周期性信號，在區(qū)域的邊界處，重力位近似等于兩邊界處重力位的平均值，會產(chǎn)生周期延拓邊界效應(yīng)。在矩諧分析中，擴(kuò)展計算區(qū)域的范圍，使計算區(qū)域內(nèi)的待求重力場信號滿足周期性條件，可降低周期延拓邊界效應(yīng)的影響。若數(shù)據(jù)點(diǎn)平面范圍為Dx和Dy（如圖1），將計算區(qū)域的平面范圍擴(kuò)展為Dx＋Δx和Dy＋Δy，使得待估重力場信號在區(qū)間D＋Δ上滿足周期性條件，可顯著降低邊界效應(yīng)，于是式（4）應(yīng)改寫為

航空重力數(shù)據(jù)在向下延拓過程中觀測噪聲會被放大。隨著數(shù)據(jù)格網(wǎng)間隔的減小和延拓高度的增加，設(shè)計矩陣A（式（7））的復(fù)共線性增強(qiáng)，有可能呈病態(tài)甚至奇異，若采用經(jīng)典最小二乘方法求解，即使較小的觀測誤差也會導(dǎo)致解的嚴(yán)重失真甚至錯誤。此外，在矩諧分析中引入擴(kuò)展參數(shù)Δ也會使法方程（式（7））的條件數(shù)增加，加重法方程的病態(tài)性，同樣會影響待求擾動位系數(shù)的求解精度。因此，必要時需引入正則化方法求解病態(tài)法方程，以獲取待求位系數(shù)的最優(yōu)解，可采用Tikhonov正則化方法，則式（6）的正則化解為

式中，I為單位陣；α＞0為正則化參數(shù)，至于最優(yōu)正則化參數(shù)的確定可采用廣義交叉檢驗（GCV），詳細(xì)算法請參見文獻(xiàn)［6］、［19］。

利用飛行高度處的重力擾動觀測值進(jìn)行矩諧分析求得表示該區(qū)域內(nèi)重力場的擾動位系數(shù)后，再按照式（1）進(jìn)行矩諧綜合計算即可得到大地水準(zhǔn)面上的重力擾動，至此便完成了航空重力觀測值的向下延拓。

3 數(shù)值計算與分析

為評價基于矩諧分析的向下延拓方法的可靠性、精度和穩(wěn)定性，利用EGM2008重力位模型設(shè)計模擬數(shù)值試驗，將矩諧分析與直接法［13］和基于廣義嶺估計的逆泊松積分法［6］這兩種下延方法進(jìn)行數(shù)值比較與分析。首先利用EGM2008重力位模型的2～2190階位系數(shù)計算得到飛行高度處的重力擾動，在模擬觀測值中加入數(shù)學(xué)期望為0、標(biāo)準(zhǔn)差為1.5mGal的高斯白噪聲，移去由GOCE衛(wèi)星重力位模型［20］2～120階位系數(shù)計算的參考重力擾動得到殘余重力擾動，將殘余重力擾動向下延拓至大地水準(zhǔn)面得到延拓值，再恢復(fù)參考模型值得到大地水準(zhǔn)面上的重力擾動延拓值。將由EGM2008模型2～2190階位系數(shù)計算的大地水準(zhǔn)面上的重力擾動作為真實值，比較重力擾動延拓值與真實值，可以檢驗向下延拓方法的可靠性與精度并比較其優(yōu)劣。為模擬真實航空重力測量條件，選定4個飛行高度H（正高），分別是2km、3km、4km、5km，數(shù)據(jù)點(diǎn)范圍為緯線［30°，33°］和經(jīng)線［107°，110°］圍成的3°×3°區(qū)域，格網(wǎng)間隔取2.5′×2.5′，共包含5184個數(shù)據(jù)點(diǎn)，重力擾動模擬觀測值統(tǒng)計信息見表1。分別采用矩諧分析、直接法和基于廣義嶺估計的逆泊松積分法，將飛行高度處的重力擾動觀測值向下延拓至大地水準(zhǔn)面。直接法和逆泊松積分法采用EGM96重力位模型計算遠(yuǎn)區(qū)流動點(diǎn)的貢獻(xiàn)，最大階數(shù)截斷至360階，積分球冠區(qū)半徑選為1°。經(jīng)多次對比測試，確定了矩諧分析的最優(yōu)展開階數(shù)和擴(kuò)展參數(shù)，展開階數(shù)取為N＝M＝20，擴(kuò)展參數(shù)取為Δ＝50km，選擇GCV方法確定正則化參數(shù)。

表1 重力擾動模擬觀測值統(tǒng)計表Tab.1 Statistics of simulated gravity disturbance

由于計算區(qū)域外沒有數(shù)據(jù)點(diǎn)分布，直接法、逆泊松積分法和矩諧分析的向下延拓結(jié)果均會受到邊界效應(yīng)影響，矩諧分析中還存在周期延拓邊界效應(yīng)，首先對比這3種下延方法的邊界效應(yīng)。圖2是直接法、基于廣義嶺估計的逆泊松積分法和矩諧分析的向下延拓誤差分布圖，延拓高度為4km，計算點(diǎn)范圍與數(shù)據(jù)點(diǎn)范圍保持一致?？梢钥闯觯谶@3種方法中，矩諧分析的下延結(jié)果受邊界效應(yīng)的影響最小，誤差極值點(diǎn)的數(shù)量和量級是最小的，主要分布在距四周邊界15′以內(nèi)，中心區(qū)域（圖2中虛線圍成的2°×2°矩形）計算點(diǎn)的結(jié)果不受邊界效應(yīng)影響。直接法的邊界效應(yīng)最為嚴(yán)重，誤差極值點(diǎn)的數(shù)量和量級均大大超過其他兩種方法，即使在中心區(qū)域誤差也非常明顯?；趶V義嶺估計的逆泊松積分法的邊界效應(yīng)介于矩諧分析和直接法之間，中心區(qū)域的計算點(diǎn)基本不受邊界效應(yīng)影響。

為消除邊界效應(yīng)的影響，應(yīng)以2°×2°中心區(qū)域的計算結(jié)果來評定3種向下延拓方法的精度，共包含2304個計算點(diǎn)。表2分別給出了取不同延拓高度時直接法、逆泊松積分法和矩諧分析的向下延拓誤差統(tǒng)計信息，向下延拓誤差定義為延拓值與由EGM2008模型計算的真實值之差，其中誤差比定義為

式中，RMS為延拓值相對于真實值的均方根誤差；e為白噪聲的標(biāo)準(zhǔn)差?？梢钥闯?，在不同延拓高度處，矩諧分析所得延拓解的精度都是最高的，RMS隨延拓高度的增大幅度也是最小的，最大誤差比僅為2.12，噪聲放大程度遠(yuǎn)低于直接法，略低于基于廣義嶺估計的逆泊松積分法。當(dāng)延拓高度為5km時，矩諧分析解的RMS為3.183mGal，占直接法所得延拓值RMS（5.665mGal）的56%，略低于逆泊松積分法所得延拓解的RMS（3.245mGal）。圖3是2°×2°中心區(qū)域內(nèi)3種方法的向下延拓誤差分布圖，延拓高度為4km。從中可以看出，矩諧分析的向下延拓誤差均在±10mGal以內(nèi)，大部分計算點(diǎn)處的誤差不超過±3mGal，其空間分布與逆泊松積分法相似，但誤差量級略低。相比之下，直接法的向下延拓效果最差，中心區(qū)域受到邊界效應(yīng)的影響，出現(xiàn)了較多誤差極值點(diǎn)。

圖2 不同方法的向下延拓誤差Fig.2 DWC errors of different methods

圖3 2°×2°中心區(qū)域不同方法的向下延拓誤差Fig.3 DWC errors of different methods in 2°×2°central area

表2 不同飛行高度的向下延拓誤差統(tǒng)計Tab.2 Statistics of the DWC errors at different flight levels

表3 加入不同量級噪聲時的向下延拓誤差統(tǒng)計Tab.3 Statistics of the DWC errors with different noise levels

以上數(shù)值結(jié)果都是在標(biāo)準(zhǔn)差為1.5mGal的高斯白噪聲條件下得到的，為檢驗含較大量級噪聲時3種向下延拓方法的可靠性與精度，分別在重力擾動觀測值中加入標(biāo)準(zhǔn)差為2mGal、2.5mGal和3mGal的高斯白噪聲進(jìn)行向下延拓計算。表3分別給出了加入不同量級噪聲時直接法、逆泊松積分法和矩諧分析的向下延拓誤差統(tǒng)計信息，延拓高度取3km，總體上仍然是矩諧分析表現(xiàn)最優(yōu)，具有最好的精度和穩(wěn)定性。數(shù)值比較結(jié)果表明矩諧分析能夠有效抑制測量噪聲的放大，實現(xiàn)航空重力信號的穩(wěn)定向下延拓，在精度、穩(wěn)定性和邊界效應(yīng)等方面都要優(yōu)于直接法和基于廣義嶺估計的逆泊松積分法。

4 結(jié) 論

矩諧分析通過在局部直角坐標(biāo)系下求解關(guān)于引力位的拉普拉斯方程，建立位系數(shù)與重力觀測值之間的函數(shù)關(guān)系，其推導(dǎo)過程嚴(yán)格滿足地球重力場的位理論，具有明確的物理意義。實際計算時矩諧函數(shù)展開至一定的截斷階次，求得表征區(qū)域重力場結(jié)構(gòu)的有限階次矩諧系數(shù)，在譜域內(nèi)對重力場信號具有一定的平滑作用，在求解關(guān)于擾動位系數(shù)的線性方程組（6）時，又引入了正則化處理，故矩諧分析能夠有效抑制向下延拓過程中觀測噪聲的放大。

直接法和逆泊松積分法屬于空域內(nèi)的向下延拓方法，由于其函數(shù)模型的固有局限性，下延結(jié)果的邊界效應(yīng)比較明顯，特別是直接法的結(jié)果受邊界效應(yīng)影響非常嚴(yán)重。矩諧分析方法同樣存在邊界效應(yīng)，但擴(kuò)展計算區(qū)域的范圍，求解法方程時進(jìn)行正則化處理，可有效控制邊界效應(yīng)的量級和影響范圍，使得中心區(qū)域不致受到邊界效應(yīng)的影響?；贓GM2008重力位模型設(shè)計的空中重力數(shù)據(jù)向下延拓數(shù)值試驗表明：在下延精度、穩(wěn)定性和邊界效應(yīng)等方面，矩諧分析都要優(yōu)于直接法和基于廣義嶺估計的逆泊松積分法，能夠?qū)崿F(xiàn)航空重力數(shù)據(jù)的穩(wěn)定下延，是一種可靠的航空重力向下延拓方法。

［1］ MARTINE C Z.Stability Investigations of a Discrete Downward Continuation Problem for Geoid Determination in the Canadian Rocky Mountains ［J］.Journal of Geodesy，1996，70：805－828.

［2］ VANíCEK P，SUN W，ONG P，et al.Downward Continuation of Helmert’s Gravity［J］.Journal of Geodesy，1996，71：21－34.

［3］ SUN W，VANíCEK P.On Some Problems of the Downward Continuation of the 5′×5′Mean Helmert Gravity Disturbance［J］.Journal of Geodesy，1998，72：411－420.

［4］ 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.（王興濤，石磐，朱非洲.航空重力測量數(shù)據(jù)向下延拓的正則化算法及其譜分解 ［J］.測繪學(xué)報，2004，33（1）：33－38.）

［5］ 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.（顧勇為，歸慶明.航空重力測量數(shù)據(jù)向下延拓基于信噪比的正則化方法的研究［J］.測繪學(xué)報，2010，39（5）：458－464.）

［6］ JIANG Tao，LI Jiancheng，WANG Zhengtao，et al.The 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］.測繪學(xué)報，2011，40（6）：684－689.）

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

［8］ FORSBERG R，KENYON S.Evaluation and Downward Continuation of Airborne Gravity Data：The Greenland Example［C］∥Proceedings of International Symposium Kinematic Systems in Geodesy，Geomatics and Navigation.Banff：［s.n.］；1994：531－538.

［9］ BáLHA T，HIRSCH M，KELLER W，et al.Application of a Spherical FFT Approach in Airborne Gravimetry［J］.Journal of Geodesy，1996，70：663－672.

［10］ HWANG C，HSIAO Y，SHIH H，et al.Geodetic and Geophysical Results from a Taiwan Airborne Gravity Survey：Data Reduction and Accuracy Assessment［J］.Journal of Geophysical Research，2007，112（2）：101－110.

［11］ 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.（石磐，王興濤.空中測量地面平均重力異常的頻域分析 ［J］.測繪學(xué)報，1995，24（4）：301－308.）

［12］ 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］.測繪學(xué)報，1999，28（3）：195－198.）

［13］ NOVáK P，HECK B.Downward Continuation and Geoid Determination Based on Band－limited Airborne Gravity Data［J］.Journal of Geodesy，2002，76：269－278.

［14］ ALLDREDGE L R.Rectangular Harmonic Analysis Applied to the Geomagnetic Field ［J］.Journal of Geophysical Research，1981，86（B4）：3021－3026.

［15］ LI Mingming，HUANG Xianlin，LU Hongqian，et al.Modeling of High Accuracy Local Geomagnetic Field Based on Rectangular Harmonic Analysis［J］.Journal of Astronautics，2010，31（7）：1730－1736.（李明明，黃顯林，盧鴻謙，等.基于矩諧分析的高精度局部地磁場建模研究［J］.宇航學(xué)報，2010，31（7）：1730－1736.）

［16］ BIAN Shaofeng.Numerical Solution for Geodetic Boundary Value Problem and the Earth′s Gravity Field Approximation［D］.Wuhan：Wuhan Technical University of Surveying and Mapping，1992.（邊少鋒.大地測量邊值問題數(shù)值解法和地球重力場逼近 ［D］.武漢：武漢測繪科技大學(xué)，1992.）

［17］ PAVLIS N K，HOLMES A，KENYON S C，et al.An Earth Gravitational Model to Degree 2160：EGM08［C］∥Proceedings of the 2008General Assembly of the European Geoscience Union.Vienna：［s.n.］，2008.

［18］ JIANG Tao.Regioanl Geoid Determination Using Airborne Gravimetry Data［D］.Wuhan：Wuhan University，2011.（蔣濤.利用航空重力測量數(shù)據(jù)確定區(qū)域大地水準(zhǔn)面［D］.武漢：武漢大學(xué)，2011.）

［19］ KUSCHE J，KLEES R.Regularization of Gravity Field Estimation from Satellite Gravity Gradients［J］.Journal of Geodesy，2002，76：359－368.

［20］ BRUINSMA S L，MARTY J C，BALMINO G，et al.GOCE Gravity Field Recovery by Means of the Direct Numerical Method［C］∥Proceedings of the ESA Living Planet Symposium.Bergen：［s.n.］，2010.