衛(wèi)星測高與船載重力測量數(shù)據(jù)融合的點質(zhì)量擬合法

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

中國測繪科學研究院，北京 100830

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

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

1 計算模型與融合方法

1.1 計算模型

點質(zhì)量模型是

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

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

X=(ATPA)-1ATPL

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

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

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

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

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

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

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

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

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

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

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

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

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

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

收集的船載重力測量數(shù)據(jù)分布如圖1中點線所示，由國內(nèi)不同部門，在不同年代進行測量得到，空間分辨率約10′?？紤]到船載重力測量數(shù)據(jù)處理不是本文的重點，對其精細處理主要借鑒了文獻[13—17]等提供的方法。該區(qū)域內(nèi)的船載重力數(shù)據(jù)有3種不同分辨率。分別是：測線間距約20 km，測點之間的間距約為3 km；測線間距約11 km，測點之間的間距約為3.7 km；測線間距約6 km，測點之間的間距約為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ù)，用來對本文融合后的重力異常模型進行檢核。其分布情況如圖2中紅色點線所示。有關(guān)NGDC數(shù)據(jù)的介紹可以參考文獻[18—19]。

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

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

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

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

表1 不同分辨率反演的點質(zhì)量模型計算的重力異常與對應測區(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ù))反演點質(zhì)量模型，然后利用點質(zhì)量模型計算測線上重力異常，并將原始值與其求差，對差值進行統(tǒng)計((2)—(5)列)；(7)—(11)列為各分區(qū)按照10 km的空間分辨率重采樣后(第(11)列為該區(qū)內(nèi)的數(shù)據(jù)點數(shù))反演點質(zhì)量模型，然后利用點質(zhì)量模型計算測線(第(6)列)上重力異常，并將原始值與其求差，對差值進行統(tǒng)計((7)—(10)列)；(12)—(17)列為各分區(qū)按照15 km的空間分辨率重采樣后(第(17)列為該區(qū)內(nèi)的數(shù)據(jù)點數(shù))反演點質(zhì)量模型，然后利用點質(zhì)量模型計算測線(第(6)列)上重力異常，并將原始值與其求差，對差值進行統(tǒng)計((13)—(17)列)；(18)—(22)列為各分區(qū)按照20 km的空間分辨率重采樣后(第(22)列為該區(qū)內(nèi)的數(shù)據(jù)點數(shù))反演點質(zhì)量模型，然后利用點質(zhì)量模型計算測線(第(6)列)上重力異常，并將原始值與其求差，對差值進行統(tǒng)計((18)—(22)列)。

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

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

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

目前，國際上發(fā)布的兩個重力場模型分別是Sandwell v23.1和DTU13模型。Sandwell v23.1模型于2014年發(fā)布，其中采用部分Cryosat-2及Jason-1的大地測量周期內(nèi)的數(shù)據(jù)，精度達到了1～2 mGal(文獻[20])。DTU13模型由DTU08逐漸發(fā)展而來，在世界范圍內(nèi)得到了廣泛應用，該模型仍在持續(xù)更新中，具體細節(jié)可見文獻[21]。在部分區(qū)域Sandwell v23.1精度要高于DTU13模型，但Sandwell v23.1與DTU13模型均不同程度地融入了船載重力數(shù)據(jù)。為了客觀評價衛(wèi)星測高重力場數(shù)據(jù)與船載重力數(shù)據(jù)融合后構(gòu)建的新重力場模型的精度，本文采用中國測繪科學研究院研制的純衛(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ù)的標準差也比本文收集的數(shù)據(jù)標準差要大。

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所示。標準差在4 mGal左右，但仍然有1～2 mGal的整體偏差。相比于純衛(wèi)星測高的重力異常模型，精度有一定的提升，平均偏差的絕對值有一定增大，標準差也有一定的改善。

表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é) 論

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

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

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

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

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

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