多源重力數(shù)據(jù)球冠諧模型抗差融合法

王 燚，姜效典，李德勇

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

1 引 言

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

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

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

2 球冠諧模型分析

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

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

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

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

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

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

式中

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

因?yàn)橹亓κ侵亓ξ谎厍蚬诎霃降奶荻龋?3］，即所以重力異常的展開式是［14］

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

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

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

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

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

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

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

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

4 數(shù)值計(jì)算

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

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

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

根據(jù)球冠諧階數(shù)nk和球冠半角α的關(guān)系式nk＝90°（Kmax＋0.5）/α－0.5［4］，以及與網(wǎng)格分辨率的關(guān)系Δθ＝180°/nk，可取截?cái)嚯A數(shù)Kmax＝22，這時(shí)球冠諧模型能夠擬合的網(wǎng)格分辨率大約為4′×4′。設(shè)σA和σT分別是衛(wèi)星重力數(shù)據(jù)和地面重力數(shù)據(jù)的觀測噪聲，若觀測噪聲不相關(guān)，則噪聲協(xié)方差矩陣為ΣΔ＝diag｛…σA…，…σT…｝，式（13）中先驗(yàn)矩陣取觀測數(shù)據(jù)噪聲協(xié)方差矩陣的逆，即在觀測噪聲σA＝8×10－5m/s2和σT＝3×10－5m/s2情況下，設(shè)衛(wèi)星重力數(shù)據(jù)和地面重力數(shù)據(jù)都存在粗差，應(yīng)用式（13）迭代計(jì)算得到融合后的球冠諧模型的系數(shù)估計(jì)值。最后通過式（10）計(jì)算球冠諧模型的重力異常，并使用地面重力數(shù)據(jù)進(jìn)行檢核，其內(nèi)外符合精度見表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)星遙測方式可以獲得全球的重力觀測值，但由于受電離層和潮汐效應(yīng)等因素的影響，數(shù)據(jù)的觀測噪聲較高，而且存在粗差。為了進(jìn)一步研究數(shù)據(jù)粗差和觀測值噪聲對融合效果的影響，首先將衛(wèi)星數(shù)據(jù)和地面數(shù)據(jù)的觀測噪聲分別設(shè)置為①σ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ù)無粗差，這時(shí)只需在式（15）中計(jì)算衛(wèi)星數(shù)據(jù)對應(yīng)的等價(jià)權(quán)陣，而地面數(shù)據(jù)對應(yīng)的權(quán)陣不變。最后根據(jù)數(shù)據(jù)噪聲和粗差的不同，設(shè)計(jì)6種計(jì)算方案。按上述方案計(jì)算的融合球冠諧模型內(nèi)外符合精度見表3和表4。

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

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

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

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

5 結(jié) 論

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

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

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

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

圖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ù)的自適應(yīng)融合處理［J］.測繪學(xué)報(bào)，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］.測繪科學(xué)與工程，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)仿真學(xué)報(bào)，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］.世界地質(zhì)，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］.測繪學(xué)報(bào)，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］.地球物理學(xué)報(bào)，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.（吳招才，劉天佑，高金耀.局部重力場球冠諧分析中的導(dǎo)數(shù)計(jì)算及應(yīng)用［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］.測繪科學(xué)，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水準(zhǔn)問題［J］.測繪科學(xué)，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水準(zhǔn)高程［J］.武漢大學(xué)學(xué)報(bào)：信息科學(xué)版，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ù)的性質(zhì)及其計(jì)算［J］.測繪學(xué)報(bào)，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)計(jì)診斷引論［M］.南京：東南大學(xué)出版社，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.（楊元喜，宋力杰，徐天河.大地測量相關(guān)觀測抗差估計(jì)理論［J］.測繪學(xué)報(bào)，2002，31（2）：95－99.）

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

［21］WANG Guisong.Theory and Application of Liner Model［M］.Hefei：Press of Anhui Education，1987.（王松桂.線性模型理論及其應(yīng)用［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］.地球?qū)W報(bào)，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.（田黔寧，耿濤，楊匯群，等.青藏高原西部重力異常剖面擬合及其地質(zhì)解釋［J］.地質(zhì)通報(bào)，2008，27（12）：2108－2116.）