王 燚，姜效典，李德勇

中國海洋大學海洋地球科學學院，山東 青島266100

1 引 言

隨著觀測技術和手段的不斷創(chuàng)新，以及相關科學理論的發(fā)展，目前已形成了海陸空全方位的地球重力場測量觀測體系。重力數(shù)據(jù)在某些地域呈現(xiàn)出多源化狀態(tài)，即在同一區(qū)域內，存在著多種不同手段獲取的重力數(shù)據(jù)。由于測量的手段不同，以及測量技術上的差異等，造成同一區(qū)域內的重力場在不同測量方法下獲取的重力數(shù)據(jù)呈現(xiàn)不一致性，甚至有相互矛盾的數(shù)據(jù)存在。必須將這些具有不同誤差特性的重力數(shù)據(jù)有效地融合在一起，并消除不同數(shù)據(jù)體之間的差別，解出統(tǒng)一、準確、可靠的局部重力基礎數(shù)據(jù)。

數(shù)據(jù)融合算法中，最小二乘配置算法具有天然的優(yōu)勢［1］，因為它能夠方便地將多種不同源的數(shù)據(jù)聯(lián)合起來求解，是多源數(shù)據(jù)融合過程中常用到的方法。但是使用最小二乘配置方法融合數(shù)據(jù)，需要計算多源數(shù)據(jù)之間的協(xié)方差函數(shù)矩陣［2］，在局部重力場中，當觀測數(shù)據(jù)較多時，數(shù)據(jù)運算量大。也有學者采用球諧模型進行多源數(shù)據(jù)的融合研究，通過迭代算法融合重力異常數(shù)據(jù)［3］。球諧函數(shù)雖然是球坐標系下滿足位理論邊界條件的譜函數(shù)，但其特征函數(shù)與局部區(qū)域不對應，不能滿足有限區(qū)域邊值問題的定解條件，因此求解得到的球諧系數(shù)不穩(wěn)定［4］。而且若要擬合分辨率為5′×5′的重力數(shù)據(jù)，根據(jù)奈奎斯特采樣定律，球諧函數(shù)需展開到2159階，總共需要計算4 672 082個球諧系數(shù)，計算量巨大。

文獻［5］成功運用球冠諧理論研究保守力場，指出運用球冠諧展開逼近局部力場的可能性；文獻［6］將球冠諧理論引入海面地形的研究；文獻［7］提出并推導了球冠諧分析法逼近區(qū)域地殼垂直形變場的模型；文獻［8］利用球冠諧分析建立了2005—2010年中國地區(qū)地磁場長期變化模型；文獻［9］利用球冠諧理論建立了我國的局部重力場模型。這些研究工作都取得了比較滿意的結果。由于在球冠諧模型中使用非整階勒讓德函數(shù)，余緯的定義域從［0，π］變換到［0，θ］，因此只需要少量的系數(shù)就可以擬合局部重力場，能夠集中反映局部重力場的高頻特征，并充分發(fā)揮局部地區(qū)重力資料采樣密集的優(yōu)勢。由此本文提出一種基于球冠諧模型，使用抗差嶺估計算法融合多源重力數(shù)據(jù)的方法，以期建立高精度的局域重力場模型。

2 球冠諧模型分析

文獻［9］從 Sturm－Liouvide方程（簡稱 S－L方程）出發(fā)，給出了局部重力場的球冠諧和函數(shù)表達式。在無源區(qū)，重力位V在球冠坐標下的拉普拉斯方程的解為

式中，a是地球平均半徑；r是計算點的地心半徑；θ與λ是球冠坐標下的余緯和經度［10］；k是球冠諧和分析中根的階序號；Kmax是球冠諧和分析最大截斷階數(shù)；是非整階勒讓德函數(shù)（m是整級次，nk（m）是非整階次）與為高斯系數(shù)（也稱為球冠諧系數(shù)）。

相對球諧分析，球冠諧分析是在地球的某一球冠部分內進行的球諧分析，因此必須滿足一定的邊界條件［11］。在球冠極點處（經度和余緯θ＝0）與球諧分析一致，但在θ＝θ0處，重力位V要滿足以下兩個條件

式中，f與g是球冠邊界上的函數(shù)，與θ無關。通常情況下，球冠極與地理極是不重合的，因此在進行球冠諧分析前，需要進行坐標變換。首先將各觀測點的大地經緯度（B，L）坐標轉換為地心經度和余緯，然后通過球面三角公式轉換為球冠極點對應的球冠坐標下的經度和余緯［12］。

文獻［5］證明可以通過選取適當?shù)幕瘮?shù)，使得當和成立時，式（2）和式（3）成立。由式（1）可以看到，重力位V與θ有關的函數(shù)只包含在函數(shù)中，因此在給定m和θ0后，只需根據(jù)下列兩式，通過確定根nk，來完成基函數(shù)的選取

滿足式（4）和式（5）的根nk是非整階的，根據(jù)給定的m和θ0計算可以得到兩組根序列，將這兩組根序列按大小排序，從m起始編號，得到根的階序號k。當k－m＝偶數(shù)時，取式（4）的根，k－m＝奇數(shù)時，取式（5）的根。根據(jù)S－L理論，這兩組根構成的基函數(shù)在組內是正交的，但組間并不正交。式（1）中非整階勒讓德函數(shù)（cosθ）在給定nk和m值，可以通過如下遞推公式計算［12］

式中

3 抗差嶺估計算法融合多源重力數(shù)據(jù)

因為重力是重力位沿球冠半徑的梯度［13］，即所以重力異常的展開式是［14］

設有一系列重力異常觀測值G＝（g1，g2，…，gn）T，球冠諧系數(shù)1，2，…，Kmax；m＝0，1，2，…，k；k－m＝偶數(shù)），則將式（9）寫成矩陣的形式

根據(jù)式（10）使用LS估計計算球冠諧系數(shù)時，由于系數(shù)矩陣A中包含復共線性關系，因此LS估計的法方程矩陣是病態(tài)的，存在接近零值的特征根。結果導致LS估計在線性無偏估計中雖然方差最小，但其均方誤差很大，降低了其參數(shù)估值的準確度和穩(wěn)定性。針對這個問題，已經有大量文獻討論使用正則化方法計算方程（10）的解。另一方面，LS估計算法本身不具備抗差性，觀測值的粗差會對估計的結果造成很大干擾，估計的穩(wěn)定性受到質疑［17］。為了克服法方程的病態(tài)性對估計的影響，又能具備抵抗數(shù)據(jù)粗差的影響，本文根據(jù)Tikhonov正則化方法構造最小化目標函數(shù)

解得球冠諧系數(shù)M的抗差估計表達式為

式中，等價權矩陣ˉP是對角陣；α是嶺參數(shù)；Λ是正定阻尼矩陣。

在應用式（9）計算重力異常時，由于只需要k－m＝偶數(shù)的一組根，其基函數(shù)是正交的，可以與球諧分析類似［18］，在式（12）、式（13）中，正定阻尼矩陣Λ設計為一個對角矩陣，其對角元素由下列公式給出

使用抗差算法通過式（13）計算球冠諧系數(shù)M時，需要迭代計算。方法如下：若已經解得第l步的系數(shù)M的估計值，根據(jù)觀測殘差v＝AMG，計算等價權陣ˉP和嶺參數(shù)α，然后根據(jù)式（13）計算第l＋1步的參數(shù)估計，最后計算相鄰兩次參數(shù)估計值的差，當差值小于指定值時，結束迭代計算，解得球冠諧系數(shù)M的估計。

等價權矩陣的計算采用GM估計方法［19］，由觀測值的殘差vi＝AM－G，根據(jù)式（15）計算得到

嶺參數(shù)α的選取是正則化方法的核心問題，α的取值不同，參數(shù)的估計值可能不同。確定嶺參數(shù)的方法很多［20］，常用的有嶺跡法、L曲線法、GCV法等。其中嶺跡法優(yōu)點是比較直觀，但是具有主觀隨意性；GCV法能夠解析地求得一個最優(yōu)嶺參數(shù)，但是該方法變化過于平緩，難以收斂；L曲線法是比較成熟的一種方法，但是其確定的嶺參數(shù)不是最優(yōu)的，而且在抗差迭代過程中，等價權陣ˉP是變化的，需要在每次迭代過程中都使用L曲線法確定嶺參數(shù)。本文利用Hoerl和Kennard從廣義嶺估計出發(fā)確定普通嶺估計中嶺參數(shù)的思想，使用法方程直接計算嶺參數(shù)［21］

4 數(shù)值計算

為了驗證上述方法的正確性，選取36.95°N—38.05°N，76.95°E—78.05°E范圍作為數(shù)據(jù)融合區(qū)域，球冠半角α約為0.71°。這里位于中國的青藏高原［22－23］，地勢復雜，海拔落差大。北部地殼厚度47～58km，中部地殼厚度67～73km。殼幔密度的差異大，低頻重力場變化顯著。通過http：∥bgi.omp.obs－mip.fr網站獲取該地區(qū)分辨率為2.5′×2.5′衛(wèi)星重力異常網格數(shù)據(jù)，共676個數(shù)據(jù)點。地面重力異常觀測值分辨率為3′×3′網格數(shù)據(jù)，共有395個點。將地面重力異常數(shù)據(jù)分為兩個點集合，其中點集S1包含354個點，用于模型融合試驗。點集S2包含41個點，用于融合模型的外符合精度檢測。精度檢驗標準采用地面重力異常觀測值與模型計算的重力異常值差值的均方根衡量。衛(wèi)星重力數(shù)據(jù)和地面實測重力數(shù)據(jù)的分布如圖1所示。將衛(wèi)星重力異常插值到地面觀測點位置并統(tǒng)計其差值，見表1。

圖1 重力測量數(shù)據(jù)點位分布Fig.1 Distribution of the data points of gravity surveying

表1 衛(wèi)星重力異常與地面重力異常的差值統(tǒng)計Tab.1 Error between the satellite gravity anomaly and the ground gravity anomaly 10－5 m/s2

根據(jù)球冠諧階數(shù)nk和球冠半角α的關系式nk＝90°（Kmax＋0.5）/α－0.5［4］，以及與網格分辨率的關系Δθ＝180°/nk，可取截斷階數(shù)Kmax＝22，這時球冠諧模型能夠擬合的網格分辨率大約為4′×4′。設σA和σT分別是衛(wèi)星重力數(shù)據(jù)和地面重力數(shù)據(jù)的觀測噪聲，若觀測噪聲不相關，則噪聲協(xié)方差矩陣為ΣΔ＝diag｛…σA…，…σT…｝，式（13）中先驗矩陣取觀測數(shù)據(jù)噪聲協(xié)方差矩陣的逆，即在觀測噪聲σA＝8×10－5m/s2和σT＝3×10－5m/s2情況下，設衛(wèi)星重力數(shù)據(jù)和地面重力數(shù)據(jù)都存在粗差，應用式（13）迭代計算得到融合后的球冠諧模型的系數(shù)估計值。最后通過式（10）計算球冠諧模型的重力異常，并使用地面重力數(shù)據(jù)進行檢核，其內外符合精度見表2。比較表1和表2中的數(shù)據(jù)，可以看到融合后的球冠諧模型要比衛(wèi)星重力異常更接近地面重力異常。

表2 模型符合精度Tab.2 Precisions of all models 10－5 m/s2

通過不同的測量手段獲取的重力觀測值，其數(shù)據(jù)精度是不同的。通常地面靜態(tài)重力測量和船載海洋重力測量的數(shù)據(jù)精度最高，但作業(yè)效率低；通過衛(wèi)星遙測方式可以獲得全球的重力觀測值，但由于受電離層和潮汐效應等因素的影響，數(shù)據(jù)的觀測噪聲較高，而且存在粗差。為了進一步研究數(shù)據(jù)粗差和觀測值噪聲對融合效果的影響，首先將衛(wèi)星數(shù)據(jù)和地面數(shù)據(jù)的觀測噪聲分別設置為①σA＝8×10－5m/s2和σT＝8×10－5m/s2；②σA＝8×10－5m/s2和σT＝5×10－5m/s2；③σA＝8×10－5m/s2和σT＝3×10－5m/s2共3種情況。其次將粗差的影響分為①衛(wèi)星數(shù)據(jù)和地面數(shù)據(jù)都存在粗差；②僅衛(wèi)星數(shù)據(jù)存在粗差，而地面數(shù)據(jù)無粗差，這時只需在式（15）中計算衛(wèi)星數(shù)據(jù)對應的等價權陣，而地面數(shù)據(jù)對應的權陣不變。最后根據(jù)數(shù)據(jù)噪聲和粗差的不同，設計6種計算方案。按上述方案計算的融合球冠諧模型內外符合精度見表3和表4。

表3 各種計算方案的內符合精度Tab.3 Internal precisions of computation result under each scheme 10－5 m/s2

表4 各種計算方案的外符合精度Tab.4 External precisions of computation result under each scheme 10－5 m/s2

從表3和表4的統(tǒng)計結果來看，隨著地面重力數(shù)據(jù)觀測噪聲的逐步下降，融合后的球冠諧模型的精度逐步提高。當衛(wèi)星重力數(shù)據(jù)和地面重力數(shù)據(jù)都存在粗差時，融合后的球冠諧模型受觀測數(shù)據(jù)噪聲的影響較大，當兩者的觀測噪聲一樣時，內符合精度為24.44×10－5m/s2，模型的精度最低。實際上這時融合后球冠諧模型更接近衛(wèi)星重力數(shù)據(jù)，模型與衛(wèi)星重力數(shù)據(jù)的差值均方根為9.4×10－5m/s2。表4中外符合精度反而比表3中的內符合精度高，原因可能是地面數(shù)據(jù)的檢測點多分布于融合區(qū)域中部，而模型與地面觀測數(shù)據(jù)誤差大的點多分布在球冠邊界位置，造成模型的內符合精度低于外符合精度，說明這時球冠諧模型存在一定的邊界效應。

當?shù)孛嬷亓?shù)據(jù)無粗差時，觀測數(shù)據(jù)噪聲對融合后的球冠諧模型的內符合精度影響不大，當兩者的觀測誤差一致時，內符合精度可以達到3.44×10－5m/s2，外符合精度也有6.44×10－5m/s2。當?shù)孛嬗^測數(shù)據(jù)噪聲最小時，模型精度最高，外符合精度可達到3.41×10－5m/s2。此時球冠諧模型重力異常（分辨率2.5′×2.5′）和地面重力異常（分辨率3′×3′）如圖2所示。

5 結 論

通過上面的數(shù)據(jù)分析，發(fā)現(xiàn)數(shù)據(jù)存在粗差以及數(shù)據(jù)的觀測噪聲大小，都對融合后的球冠諧模型精度有一定的影響，從中可以得出以下幾點結論。

（1）當衛(wèi)星和地面重力數(shù)據(jù)都存在粗差時，模型的精度隨著地面數(shù)據(jù)的觀測噪聲降低而顯著提高。當?shù)孛鏀?shù)據(jù)的觀測噪聲相對較小時，模型的精度受數(shù)據(jù)粗差的影響較小。

（2）當?shù)孛嬷亓?shù)據(jù)不存在粗差時，模型的精度雖然也隨著地面數(shù)據(jù)的觀測噪聲降低而提高，但受其影響較小。在兩者觀測噪聲一致時，模型外符合精度也有6.44×10－5m/s2。而當?shù)孛鏀?shù)據(jù)的觀測噪聲進一步減小到3×10－5m/s2時，模型的外符合精度可以達到3.41×10－5m/s2，如圖2所示能夠很好地擬合地面重力數(shù)據(jù)。

綜上所述，在局部重力場中，使用多種數(shù)據(jù)源，通過抗差嶺估計迭代算法構建球的冠諧模型，既能克服法矩陣病態(tài)性的影響，又具有抵抗粗差的能力。該模型能夠利用較低分辨率的高質量局部區(qū)域重力數(shù)據(jù)和較高分辨率的含有粗差的衛(wèi)星重力數(shù)據(jù)，構建高精度的局部重力場，在理論和實踐上都有重要意義。

圖2 模型與地面重力異常Fig.2 The model and the ground gravity anomaly

［1］ZHAI Zhenhe，SUN Zhongmiao.The Adaptive Fusion of Multi－source Gravity Data in Bohai Gulf［J］.Acta Geodaetica et Cartographica Sinica，2010，39（5）：444－449.（翟振和，孫中苗.渤海灣多源重力數(shù)據(jù)的自適應融合處理［J］.測繪學報，2010，39（5）：444－449.）

［2］ZHAI Zhenhe，SUN Zhongmiao.Combining Airborne Gravity Data and Terrestrial Gravity Data Based on Least Squares Collocation［J］.Geomatic Science and Engineering，2008，28（1）：10－14.（翟振和，孫中苗.基于最小二乘配置的航空和地面重力數(shù)據(jù)融合處理［J］.測繪科學與工程，2008，28（1）：10－14.）

［3］HAO Yanling，CHENG Yi，LIU Fanming，et al.Simulation of Combination Algorithm for Heterogeneous Marine Gravity Data［J］.Journal of System Simulation，2007，19（21）：4897－4900.（郝燕玲，成怡，劉繁明，等.融合多類型海洋重力數(shù)據(jù)算法仿真研究［J］.系統(tǒng)仿真學報，2007，19（21）：4897－4900.）

［4］WANG Xichen，WANG Guangjie，LI Tonglin，et al.Extracting Different Wavelength Gravity Anomaly by the Spherical Cap Harmonic Analysis Method［J］.Global Geology，1996，15（3）：80－83.（王喜臣，王光杰，李桐林，等.利用球冠諧分析方法提取不同波長重力場異常［J］.世界地質，1996，15（3）：80－83.）

［5］HAINES G V.Spherical Cap Harmonic Analysis［J］.Journal of Geophysical Research，1985，90（B3）：2583－2591.

［6］HWANG C，CHEN S K.Fully Normalized Spherical Cap Harmonics：Application to the Analysis of Sea－level Data from TOPEX/POSEIDON and ERS－1［J］.Geophysical Journal International，1997，129（2）：450－460.

［7］ZHANG Qin，ZHAO Chaoying.The Spherical Cap Harmonic Analysis Method for Crust Vertical Deformation Field Fitting［J］，Acta Geodaetica et Cartographica Sinica，2004，33（1）：39－42.（張勤，趙超英.地殼垂直形變場逼近的球冠諧分析法［J］.測繪學報，2004，33（1）：39－42.）

［8］CHEN Bin，GU Zuowen，GAO Jintian，et al.Analyses of Geomagnetic Field and Its Secular Variation over China for 2005.0Epoch Using Spherical Cap Harmonic Method［J］.Chinese Journal of Geophysics，2011，54（3）：771－779.（陳斌，顧左文，高金田，等.2005.0年代中國地區(qū)地磁場及其長期變化球冠諧和分析［J］.地球物理學報，2011，54（3）：771－779.）

［9］LI Jiancheng.Spherical Cap Harmonic Expansion for Local Gravity Field Representation［J］.Manuscripta Geodaetica，1995（20）：265－277.

［10］YOUNIS G K A，J?GER R，BECKER M.Transformation of Global Spherical Harmonic Models of the Gravity Field to a Local Adjusted Spherical Cap Harmonic Model［J］.Arabian Journal of Geosciences，2013，6（2）：375－381.

［11］WU Zhaocai，LIU Tianyou，GAO Jinyao.Derivative Calculation and Applicationin Spherical Cap Harmonic Analysis of Local Gravity Field［J］.Oceanologia et Limnologia Sinica，2006，37（6）：488－492.（吳招才，劉天佑，高金耀.局部重力場球冠諧分析中的導數(shù)計算及應用［J］.海洋與湖沼，2006，37（6）：488－492.）

［12］ZHAO Jianhu，WANG Shengping，LIU Hui，et al.Study on Establishing Local Geomagnetic Model Using Spherical Cap Harmonic Analysis［J］.Science of Surveying and Mapping，2010，35（1）：51－53.（趙建虎，王勝平，劉輝，等.海洋局域地磁場球冠諧分析建模方法研究［J］.測繪科學，2010，35（1）：51－53.）

［13］WANG Sanjun.The Spherical Cap Harmonic Analysis of GPS Leveling Height［J］.Science of Surveying and Mapping，2008，33（2）：13－14.（王三軍.利用球冠諧模型研究 GPS水準問題［J］.測繪科學，2008，33（2）：13－14.）

［14］CAO Yueling，WANG Jiexian.Application of Spherical Cap Harmonic Analysis to Fit GPS Level Height［J］.Geomatics and Information Science of Wuhan University，2008，33（7）：740－743.（曹月玲，王解先.利用球冠諧分析擬合GPS水準高程［J］.武漢大學學報：信息科學版，2008，33（7）：740－743.）

［15］DESANTIS A.Translated Origin Spherical Cap Harmonic Analysis［J］.Geophysical Journal International，1991，123（6）：253－263.

［16］PENG Fuqing，YU Jinhai.The Characters and Computation of Legendre Function with Non－integral Degree in SCHA［J］.Acta Geodaetica et Cartographica Sinica，2000，29（3）：204－208.（彭富清，于錦海.球冠諧分析中非整階Legendre函數(shù)的性質及其計算［J］.測繪學報，2000，29（3）：204－208.）

［17］WEI Bocheng，LU Guobin，SHI Jianging.Introduction to Statistical Diagnostics［M］.Nanjing：Publishing House of Dongnan University，1991.（韋博成，魯國斌，史建清.統(tǒng)計診斷引論［M］.南京：東南大學出版社，1991.）

［18］KORTE M，HOLME R.Regularization of Spherical Cap Harmonics［J］.Geophysical Journal International，2003，153（1）：253－262.

［19］YANG Yuanxi，SONG Lijie，XU Tianhe.Robust Parameter Estimation for Geodetic Correlated Observations［J］.Acta Geodaetica et Cartographica Sinica，2002，31（2）：95－99.（楊元喜，宋力杰，徐天河.大地測量相關觀測抗差估計理論［J］.測繪學報，2002，31（2）：95－99.）

［20］YANG Yuanxi.Robust Estimation and Its Applications［M］.Beijing：Bayi Publishing House，1993.（楊元喜.抗差估計理論及其應用［M］.北京：八一出版社，1993.）

［21］WANG Guisong.Theory and Application of Liner Model［M］.Hefei：Press of Anhui Education，1987.（王松桂.線性模型理論及其應用［M］.合肥：安徽教育出版社，1987.）

［22］GAO Rui，XIONG Xiaosong，LI Qiusheng，et al.The Moho Depth of Qinghai－Tibet Plateau Revealed by Seismic Detection［J］.Acta Geoscientica Sinica，2009，30（6）：761－773.（高銳，熊小松，李秋生，等.由地震探測揭示的青藏高原莫霍面深度［J］.地球學報，2009，30（6）：761－773.）

［23］TIAN Qianning，GENG Tao，YANG Huiqun，et al.Gravity Anomalies Forward Fitting and Geological Explanation in the Western Qinghai－Tibet Plateau［J］.Geological Bulletin of China，2008，27（12）：2108－2116.（田黔寧，耿濤，楊匯群，等.青藏高原西部重力異常剖面擬合及其地質解釋［J］.地質通報，2008，27（12）：2108－2116.）