柯寶貴，張利明，章傳銀，黨亞民

中國測繪科學(xué)研究院，北京 100830

高精度的海域重力場對水下潛器的導(dǎo)航、飛行器的軌道確定及海洋航空重力測量的檢核有重要作用。隨著海洋地質(zhì)調(diào)查資料的增多，不同測高衛(wèi)星的數(shù)據(jù)積累，建立更高精度和分辨率的重力場模型的需求也越來越迫切。融合船載重力數(shù)據(jù)與衛(wèi)星測高數(shù)據(jù)既可以充分發(fā)揮船載數(shù)據(jù)精度高，還可以利用衛(wèi)星測高反演的重力場分布均勻的優(yōu)點(diǎn)。對于多源重力數(shù)據(jù)融合問題，常見的方法有最小二乘配置法，也可采用文獻(xiàn)[1]中的方法，即先利用殘差重力異常修正重力位模型系數(shù)，進(jìn)而融合衛(wèi)星、航空、地面(海面)重力數(shù)據(jù)的迭代計算。其他的方法還有泊松小波徑向基函數(shù)[2]、點(diǎn)質(zhì)量方法[3]及小波分解的數(shù)據(jù)融合方法[4]。點(diǎn)質(zhì)量方法理論模型直接由反映重力場性質(zhì)的模型推導(dǎo)而來。在眾多的利用點(diǎn)質(zhì)量法研究中，研究重點(diǎn)通常是在討論如何克服法方程矩陣求逆導(dǎo)致數(shù)值解算過程失穩(wěn)的問題[5-8]。

本文不研究點(diǎn)質(zhì)量方法的穩(wěn)定性問題，主要利用點(diǎn)質(zhì)量法對于重力場的解析表達(dá)形式，借鑒移去-恢復(fù)技術(shù)的思路，通過點(diǎn)質(zhì)量法擬合船載重力測量數(shù)據(jù)的中低頻信號，將測高重力數(shù)據(jù)中的高頻信號與船載重力數(shù)據(jù)中的高頻信號使用加權(quán)最小曲率方法進(jìn)行格網(wǎng)化，然后采用點(diǎn)質(zhì)量模型恢復(fù)其中重力數(shù)據(jù)的中低頻信號，實(shí)現(xiàn)船載重力數(shù)據(jù)與衛(wèi)星測高數(shù)據(jù)的融合，最后使用國際船載重力數(shù)據(jù)進(jìn)行檢核。

1 計算模型與融合方法

1.1 計算模型

點(diǎn)質(zhì)量模型是

式中,ri、φi、λi表示第i個重力異常；Δgi為觀測點(diǎn)的球面三維坐標(biāo)；Rj、φj、λj表示第j個質(zhì)量源的球面三維坐標(biāo)；M是質(zhì)量源的總數(shù)；lij、ψij分別是第i個重力異常觀測點(diǎn)與第j個質(zhì)量源之間的空間距離及其與地心連線的夾角；G是萬有引力常數(shù)；δMj是第j個質(zhì)量源的質(zhì)量；Rj是Bjerhammar球半徑。點(diǎn)質(zhì)量模型是基于重力場等效原理而建立的，擬合重力異常數(shù)據(jù)時對重力場物理屬性表征更多一些，有別于純數(shù)學(xué)模型的擬合。

點(diǎn)質(zhì)量模型的最小二乘解為

X=(ATPA)-1ATPL

式中，P為觀測值向量L的權(quán)矩陣；L由不同類型的觀測量組成，按定權(quán)公式確定它們的權(quán)比。本文僅使用船載重力數(shù)據(jù)，故P為單位矩陣。

1.2 測高重力異常與船載重力異常融合方法

為了實(shí)現(xiàn)測高重力異常與船載重力異常融合，本文設(shè)計了如下的步驟：

第1步：考慮到船載重力數(shù)據(jù)的精度比測高的重力異常精度要高，用點(diǎn)質(zhì)量法對收集到的全部船載重力測量數(shù)據(jù)(61 226個點(diǎn)值)進(jìn)行擬合。

第2步：利用得到的質(zhì)量源三維坐標(biāo)位置及質(zhì)量大小計算船載重力測線上重力數(shù)據(jù)重力異常gship，并與船載重力異常求差，得到殘差rgs。

第3步：同樣利用得到的質(zhì)量源三維坐標(biāo)位置及質(zhì)量大小計算測高重力格網(wǎng)模型處的重力異常galt，并計算殘差rga。

第4步：分別計算rgs與rga在重合點(diǎn)處的方差，取船載重力殘差方差作為單位權(quán)中誤差，計算測高重力殘差的權(quán)因子。

第5步：對rgs與rga利用加權(quán)最小曲率方法進(jìn)行格網(wǎng)化，得到殘差重力異常模型g1。1974年，文獻(xiàn)[9]首次提出最小曲率網(wǎng)格化方法。經(jīng)過前人的對比分析[10-12]，認(rèn)為該方法的優(yōu)點(diǎn)是所需參數(shù)少、計算工作量小、有利于稀疏數(shù)據(jù)的網(wǎng)格化，且插值結(jié)果能夠保證二階導(dǎo)數(shù)連續(xù)。因此，本文采用加權(quán)最小曲率方法對離散數(shù)據(jù)進(jìn)行格網(wǎng)化。

第6步：利用得到的質(zhì)量源三維坐標(biāo)位置及質(zhì)量大小計算格網(wǎng)點(diǎn)處的重力異常g2。

第7步：將g1與g2相加即為融合的重力異常模型。

第8步：利用船載重力測量數(shù)據(jù)對模型進(jìn)行檢核。

從公開發(fā)表的文獻(xiàn)來看，對于重力場量進(jìn)行融合通常是基于不同重力場量(如重力異常、高程異常和垂線偏差等)，采用點(diǎn)質(zhì)量模型，解算出點(diǎn)質(zhì)量的三維坐標(biāo)位置及質(zhì)量大小來實(shí)現(xiàn)的。本文對此步驟進(jìn)行了改進(jìn)，即在第2步至第7步之間借鑒移去-恢復(fù)的思路，利用點(diǎn)質(zhì)量模型計算測高與船載重力中所包含的相同的中低頻信號，僅對測高與船載重力異常的高頻信號進(jìn)行融合。在移去-恢復(fù)的過程中，未采用全球重力場模型移去-恢復(fù)中低頻的重力場信號，而是基于點(diǎn)質(zhì)量模型反映出的局部重力場的中低頻信號來進(jìn)行移去-恢復(fù)運(yùn)算。這在一定的程度上可避免全球重力場模型在擬合局部重力數(shù)據(jù)時存在系統(tǒng)偏差的影響。

2 實(shí)測數(shù)值計算與結(jié)果分析

2.1 船載重力數(shù)據(jù)

收集的船載重力測量數(shù)據(jù)分布如圖1中點(diǎn)線所示，由國內(nèi)不同部門，在不同年代進(jìn)行測量得到，空間分辨率約10′?？紤]到船載重力測量數(shù)據(jù)處理不是本文的重點(diǎn)，對其精細(xì)處理主要借鑒了文獻(xiàn)[13—17]等提供的方法。該區(qū)域內(nèi)的船載重力數(shù)據(jù)有3種不同分辨率。分別是：測線間距約20 km，測點(diǎn)之間的間距約為3 km；測線間距約11 km，測點(diǎn)之間的間距約為3.7 km；測線間距約6 km，測點(diǎn)之間的間距約為4 km；主要以第1種與第2種數(shù)據(jù)為主。第3種主要分布在海底地形起伏較大區(qū)域。整個研究范圍內(nèi)數(shù)據(jù)分為50個測區(qū)，編號分別為100—149。

采用美國地球物理數(shù)據(jù)中心NGDC(National Geophysical Data Center)網(wǎng)站提供的數(shù)據(jù)作為第三方檢核數(shù)據(jù)，用來對本文融合后的重力異常模型進(jìn)行檢核。其分布情況如圖2中紅色點(diǎn)線所示。有關(guān)NGDC數(shù)據(jù)的介紹可以參考文獻(xiàn)[18—19]。

2.2 點(diǎn)質(zhì)量源的數(shù)量與位置的探討

盡管本文討論的重點(diǎn)不是點(diǎn)質(zhì)量法的穩(wěn)定性問題，但有必要討論點(diǎn)質(zhì)量源的數(shù)量、位置(平面坐標(biāo)與埋深)與原始觀測值偏離情況。為此，根據(jù)GEBCO海底地形模型(general bathymetric chart of the oceans ，分辨率為1′)依據(jù)點(diǎn)質(zhì)量法的正演模型，討論點(diǎn)質(zhì)量源的平面位置、埋深以及質(zhì)量源的數(shù)量。

一般而言，點(diǎn)質(zhì)量反演方法的埋深初值可以根據(jù)數(shù)據(jù)的地球物理先驗信息給定，在反演過程進(jìn)行精化。本文根據(jù)可用的數(shù)據(jù)情況，假定點(diǎn)質(zhì)量源分布于海底，即依據(jù)GEBCO海底地形模型內(nèi)插測線上各點(diǎn)重力異常數(shù)據(jù)的水深值，將該值作為對應(yīng)點(diǎn)處的質(zhì)量埋深初值。

本文在任意測區(qū)的測線方向，按照10、15和20 km的分辨率對測線數(shù)據(jù)進(jìn)行了重采樣。將重采樣的數(shù)據(jù)，反演對應(yīng)測區(qū)的點(diǎn)質(zhì)量模型，進(jìn)而由點(diǎn)質(zhì)量模型計算的重力異常與全部測線上重力異常的差異來比較，不同分辨率數(shù)據(jù)的差別。為了文章的簡潔，表1列舉部分測區(qū)采用不同分辨率的數(shù)據(jù)與原始數(shù)據(jù)的差異。

表1 不同分辨率反演的點(diǎn)質(zhì)量模型計算的重力異常與對應(yīng)測區(qū)的重力異常差異

表1中，(1)、(12)列是研究區(qū)內(nèi)數(shù)據(jù)的分區(qū)編號；(2)—(6)列為各分區(qū)采用該區(qū)內(nèi)全部的數(shù)據(jù)作為初始數(shù)據(jù)(第(6)列為該區(qū)內(nèi)的數(shù)據(jù)總數(shù))反演點(diǎn)質(zhì)量模型，然后利用點(diǎn)質(zhì)量模型計算測線上重力異常，并將原始值與其求差，對差值進(jìn)行統(tǒng)計((2)—(5)列)；(7)—(11)列為各分區(qū)按照10 km的空間分辨率重采樣后(第(11)列為該區(qū)內(nèi)的數(shù)據(jù)點(diǎn)數(shù))反演點(diǎn)質(zhì)量模型，然后利用點(diǎn)質(zhì)量模型計算測線(第(6)列)上重力異常，并將原始值與其求差，對差值進(jìn)行統(tǒng)計((7)—(10)列)；(12)—(17)列為各分區(qū)按照15 km的空間分辨率重采樣后(第(17)列為該區(qū)內(nèi)的數(shù)據(jù)點(diǎn)數(shù))反演點(diǎn)質(zhì)量模型，然后利用點(diǎn)質(zhì)量模型計算測線(第(6)列)上重力異常，并將原始值與其求差，對差值進(jìn)行統(tǒng)計((13)—(17)列)；(18)—(22)列為各分區(qū)按照20 km的空間分辨率重采樣后(第(22)列為該區(qū)內(nèi)的數(shù)據(jù)點(diǎn)數(shù))反演點(diǎn)質(zhì)量模型，然后利用點(diǎn)質(zhì)量模型計算測線(第(6)列)上重力異常，并將原始值與其求差，對差值進(jìn)行統(tǒng)計((18)—(22)列)。

從表1可以看出，采用不同分辨率的采樣數(shù)據(jù)，反演的結(jié)果相差不大，也說明點(diǎn)質(zhì)量法模型對數(shù)據(jù)的質(zhì)量不敏感。需要說明的是，采樣數(shù)據(jù)需要包含測線上重力異常的基本特征(最大值、最小值以及變化趨勢)，即數(shù)據(jù)分布不能太稀疏。

為此，本文采用了20 km的分辨率對所有的測線數(shù)據(jù)進(jìn)行重采樣，依據(jù)GEBCO模型的深度值，以及測線的平面坐標(biāo)反演點(diǎn)質(zhì)量源大小，作為后文研究的基礎(chǔ)。

2.3 已有海洋重力場模型的比較

目前，國際上發(fā)布的兩個重力場模型分別是Sandwell v23.1和DTU13模型。Sandwell v23.1模型于2014年發(fā)布，其中采用部分Cryosat-2及Jason-1的大地測量周期內(nèi)的數(shù)據(jù)，精度達(dá)到了1～2 mGal(文獻(xiàn)[20])。DTU13模型由DTU08逐漸發(fā)展而來，在世界范圍內(nèi)得到了廣泛應(yīng)用，該模型仍在持續(xù)更新中，具體細(xì)節(jié)可見文獻(xiàn)[21]。在部分區(qū)域Sandwell v23.1精度要高于DTU13模型，但Sandwell v23.1與DTU13模型均不同程度地融入了船載重力數(shù)據(jù)。為了客觀評價衛(wèi)星測高重力場數(shù)據(jù)與船載重力數(shù)據(jù)融合后構(gòu)建的新重力場模型的精度，本文采用中國測繪科學(xué)研究院研制的純衛(wèi)星測高重力場模型(CASM-ALT)。該模型與本文中研究區(qū)域內(nèi)收集的船載重力異常(非NGDC數(shù)據(jù))的差異如表2中第2行所示。表2第3行為Sandwell v23.1模型與船載數(shù)據(jù)的差異，第4行為DTU13模型與船載數(shù)據(jù)的差異。從表中可以看出，CASM-ALT精度略低于Sandwell v23.1、DTU13模型重力異常。在后文中主要以NGDC數(shù)據(jù)作為檢核。表2中也給出了NGDC船載重力數(shù)據(jù)與各模型的差異。

表2 CASM-ALT 、Sandwell v23.1、DTU13模型與船載重力數(shù)據(jù)的偏差

從表2可以看出，本文收集的數(shù)據(jù)與NGDC的數(shù)據(jù)有一定的系統(tǒng)差。另外由于研究區(qū)域位于馬里亞納海溝東側(cè) ，該處海底地形起伏較為劇烈(圖1中的背景)，且NGDC數(shù)據(jù)的標(biāo)準(zhǔn)差也比本文收集的數(shù)據(jù)標(biāo)準(zhǔn)差要大。

2.4 融合的重力異常模型檢核

基于上節(jié)的步驟，得到融合衛(wèi)星測高與船載重力數(shù)據(jù)的重力異常模型。圖2中紅線表示NGDC的國際船載重力測量數(shù)據(jù)?？紤]到區(qū)域A、區(qū)域B和區(qū)域C的海底有起伏的海山，高頻信號豐富，選擇這幾個區(qū)域作為檢核區(qū)域更能反映融合后重力場模型的精度。

檢核區(qū)域A、B、C的海底地形及NGDC船載重力檢核數(shù)據(jù)的分布如圖3—圖5所示。相對于整個研究區(qū)域而言，所選的檢核區(qū)域是海底地形起伏較大的區(qū)域。對融合后的重力異常檢核結(jié)果如表3所示。標(biāo)準(zhǔn)差在4 mGal左右，但仍然有1～2 mGal的整體偏差。相比于純衛(wèi)星測高的重力異常模型，精度有一定的提升，平均偏差的絕對值有一定增大，標(biāo)準(zhǔn)差也有一定的改善。

表3 檢核區(qū)內(nèi)融合后的重力異常模型與NGDC測線重力的差異

圖1 研究區(qū)內(nèi)船載重力數(shù)據(jù)的分布Fig.1 The distribution of shipborne gravity anomaly

圖2 NGDC國際船載重力測量數(shù)據(jù)分布及檢核區(qū)的分布Fig 2 The distribution of the trackline gravity anomaly of NGDC and the validation regions

圖3 檢核區(qū)A的海底地形及船載重力數(shù)據(jù)分布Fig.3 Bathymetry and shipborne gravity data distribution of region A

圖4 檢核區(qū)B的海底地形及船載重力數(shù)據(jù)分布Fig.4 Bathymetry and shipborne gravity data distribution of region B

圖5 檢核區(qū)C的海底地形及船載重力數(shù)據(jù)分布Fig.5 Bathymetry and shipborne gravity data distribution of region C

3 結(jié) 論

基于點(diǎn)質(zhì)量法擬合船載重力異常數(shù)據(jù)的中、低頻信號，在衛(wèi)星測高重力異常中將其移去。對船載重力數(shù)據(jù)的高頻信號及衛(wèi)星測高的高頻信號進(jìn)行加權(quán)最小曲率格網(wǎng)化，基于點(diǎn)質(zhì)量模型恢復(fù)中低頻信號，實(shí)現(xiàn)二者的融合。經(jīng)國際船載重力測量數(shù)據(jù)檢核，研究結(jié)果表明：

(1) 文中采用的融合步驟相比以往常用的數(shù)據(jù)融合步驟進(jìn)行了改進(jìn)，即引進(jìn)了基于點(diǎn)質(zhì)量模型移去-恢復(fù)技術(shù)的融合方法。

(2) 采用GEBCO的海底地形先驗?zāi)Ｐ?，指定點(diǎn)質(zhì)量的平面坐標(biāo)位置，減少反演參數(shù)，回避了點(diǎn)質(zhì)量模型法方程求逆失穩(wěn)的問題。

(3) 在移去-恢復(fù)的過程中，沒有使用全球重力場模型作為參考，而是依據(jù)點(diǎn)質(zhì)量模型參數(shù)所具有的局部重力場特征，計算區(qū)域重力異常的中低頻分量，這有助于避免引入全球重力場模型的誤差。

(4) 對不同平臺下獲取的重力異常的高頻信息進(jìn)行融合，所得的融合后重力場模型的精度有一定提高，成果標(biāo)準(zhǔn)差約為4 mGal，反映了本文的技術(shù)思路及方法具有一定的可行性與可靠性。

(5) 構(gòu)建的融合重力場模型依然與NGDC重力異常數(shù)據(jù)存在1～2 mGal平均偏差，可能是因為重力數(shù)據(jù)之間的系統(tǒng)偏差或者由于不同年代進(jìn)行重力測量時測量儀器之間存在偏差造成的，其實(shí)際成因擬在后續(xù)研究做進(jìn)一步分析。