翟振和，孫中苗，李迎春，郭飛霄

(西安測繪研究所，陜西 西安 710054)

The Integrative Method of Airborne Disturbing Gravity Data and

Land Gravity Anomaly Data

ZHAI Zhenhe，SUN Zhongmiao，LI Yingchun，GUO Feixiao

?

航空和地面重力測量數(shù)據(jù)的一體化處理方法

翟振和，孫中苗，李迎春，郭飛霄

(西安測繪研究所，陜西 西安 710054)

The Integrative Method of Airborne Disturbing Gravity Data and

Land Gravity Anomaly Data

ZHAI Zhenhe，SUN Zhongmiao，LI Yingchun，GUO Feixiao

摘要：在三維經(jīng)典擾動位協(xié)方差模型基礎(chǔ)上，推導(dǎo)獲得局部重力異常協(xié)方差模型、局部擾動重力協(xié)方差模型，以及擾動重力和重力異常的互協(xié)方差模型，這3個模型在忽略小項后是一致的，即航空擾動重力數(shù)據(jù)與地面重力異常數(shù)據(jù)的聯(lián)合處理采用一個模型即可，大大簡化了處理流程。以最小二乘配置理論及本文推導(dǎo)的協(xié)方差模型為基礎(chǔ)，構(gòu)建了航空擾動重力測線數(shù)據(jù)與地面重力異常數(shù)據(jù)一體化處理的具體方法，利用該方法可推估地面待求點的重力異常數(shù)據(jù)，且處理過程自動完成了航空擾動重力數(shù)據(jù)的向下延拓，同時也大大削弱了航空數(shù)據(jù)中存在的系統(tǒng)偏差影響。在某試驗區(qū)開展了一體化處理試驗，結(jié)果表明，采用一體化處理得到的重力異常數(shù)據(jù)比單獨使用5′分辨率航空重力數(shù)據(jù)推估獲得的數(shù)據(jù)精度提高2.4 mGal，比單獨使用2.5′分辨率地面重力數(shù)據(jù)推估獲得的數(shù)據(jù)精度提高1.2 mGal。

關(guān)鍵詞：航空擾動重力數(shù)據(jù); 地面重力異常數(shù)據(jù)；協(xié)方差模型；一體化處理

一、引言

目前，利用航空重力測量技術(shù)獲取局部區(qū)域的高分辨率重力場信息成為一個非常有效的手段。在航空測量數(shù)據(jù)處理中，傳統(tǒng)的技術(shù)流程是先將觀測得到的航空擾動重力數(shù)據(jù)進(jìn)行濾波、格網(wǎng)化，然后轉(zhuǎn)化為重力異常數(shù)據(jù)，最后向下延拓至對應(yīng)的地面區(qū)域獲得最終的地面重力異常數(shù)據(jù)[1]。上述數(shù)據(jù)處理過程往往會使測量數(shù)據(jù)受到不當(dāng)?shù)奶幚?如采用不準(zhǔn)確的模型)而影響最終數(shù)據(jù)的質(zhì)量(如存在一定的系統(tǒng)誤差)。如果能找到一種直接利用測線擾動重力數(shù)據(jù)進(jìn)行處理的方法則可以在減少處理流程的同時提高產(chǎn)品數(shù)據(jù)的精度。在航空重力測量實施的區(qū)域，地面上或多或少會有已知的重力異常測量點，這些測量點的重力數(shù)據(jù)較之航空測量精度更高，因此，如何將地面已知重力數(shù)據(jù)有效利用起來構(gòu)建航空和地面一體化的處理方法成為國內(nèi)外學(xué)者關(guān)注的重點。最小二乘配置法具有聯(lián)合多種數(shù)據(jù)類型的處理優(yōu)勢[2-5]，然而在運用最小二乘配置法進(jìn)行處理時，多數(shù)學(xué)者都是在固定的協(xié)方差模型及模型參數(shù)下進(jìn)行解算，即所有觀測數(shù)據(jù)都采用同一個先驗信息，而后按照該先驗信息完成數(shù)據(jù)的推估和濾波[6-8]。針對上述問題，本文將在最小二乘配置法基礎(chǔ)上，構(gòu)建一種能夠?qū)娇蘸偷孛嬷亓?shù)據(jù)進(jìn)行一體化處理的方法，該方法應(yīng)顧及不同高度和不同數(shù)據(jù)類型，且兼具格網(wǎng)化和向下延拓等功能。

二、局部重力異常的協(xié)方差模型

構(gòu)建能夠反映局部重力場空間變化的重力異常協(xié)方差模型是進(jìn)行一體化處理的基礎(chǔ)工作。根據(jù)文獻(xiàn)[3]的研究成果，采用擾動位的三維協(xié)方差經(jīng)驗?zāi)Ｐ腿缦?/p>

(1)

考慮到局部區(qū)域的應(yīng)用環(huán)境，將式(1)用大地坐標(biāo)替換得到如下形式的局部擾動位協(xié)方差模型

(2)

式中，f、α具體含義同上式；Hi=Hp+Hq+Di，Hp、Hq代表p、q兩點的大地高。

根據(jù)擾動場元間的線性關(guān)系，按照協(xié)方差傳播定律可得p、q兩點重力異常間協(xié)方差為

(3)

式中，CTT為擾動位協(xié)方差模型；ρp、ρq分別是p、q兩點的地心向徑。

(4)

對擾動位協(xié)方差模型的通用形式CTT求導(dǎo)得

(5)

二階導(dǎo)數(shù)為

(6)

將式(5)、式(6)代入式(3)得到局部重力異常協(xié)方差的一般形式為

(7)

其完整形式為

(8)

(9)

(10)

(11)

(12)

模型的完整形式為

(13)

式(9)—式(13)即為簡化后的局部重力異常的協(xié)方差實用模型。

三、局部擾動重力和重力異常的互協(xié)方差模型

擾動重力與擾動位間的關(guān)系為

(14)

考慮到橢球面上大地高方向H與n方向相同，因此

(15)

按照協(xié)方差傳播定律可得空間擾動重力觀測點p和地面重力異常待求點q的互協(xié)方差為

(16)

將局部擾動位協(xié)方差模型CTT代入上式并忽略小項得到互協(xié)方差模型為

(17)

式(17)即擾動重力與重力異?；f(xié)方差模型的完整形式。

同理，局部擾動重力的自協(xié)方差模型按照上述原理推導(dǎo)如下

(18)

將式(18)、式(17)與式(13)對比發(fā)現(xiàn)，擾動重力與重力異常互協(xié)方差模型與重力異常協(xié)方差模型、擾動重力協(xié)方差模型形式上完全一致，這樣在航空擾動重力數(shù)據(jù)和地面重力異常數(shù)據(jù)處理時只采用同一個協(xié)方差模型即可。

四、一體化處理方法

假設(shè)航空擾動重力數(shù)據(jù)觀測值為L1，地面重力異常數(shù)據(jù)觀測值為L2，對應(yīng)的觀測信號為t1、t2，對應(yīng)的觀測噪聲分別為Δ1、Δ2，航空重力數(shù)據(jù)對應(yīng)的系統(tǒng)性參數(shù)為X，此時得到兩組觀測方程分別為

L1=A1X+t1+Δ1

(19)

L2=t2+Δ2

(20)

式中，A1為系數(shù)陣；觀測值的協(xié)方差陣為

(21)

式中，C11、C22表示兩組數(shù)據(jù)的自協(xié)方差；C21、C12表示兩組數(shù)據(jù)的互協(xié)方差。

未測信號S(地面未測點的重力異常)與已測信號t的互協(xié)方差為

(22)

按照經(jīng)典理論可推導(dǎo)獲得未測信號的推估公式[9-11]如下

(23)

(24)

(25)

式中

(26)

(27)

(28)

(29)

系統(tǒng)參數(shù)的估值為

(30)

根據(jù)大量的試驗分析和理論推導(dǎo)，采用以下計算步驟。對于兩類數(shù)據(jù)的自協(xié)方差，分別采用各自的參數(shù)計算即Δg1的自協(xié)方差C11為

(31)

Δg2的自協(xié)方差為

(32)

對于兩類數(shù)據(jù)的互協(xié)方差的計算，首先加入兩個調(diào)整因子Pa、Pb對模型參數(shù)進(jìn)行調(diào)整，為了保證加入調(diào)整因子后的參數(shù)符合實際情況，Pa、Pb應(yīng)滿足如下關(guān)系

Pa+Pb=1

(33)

利用Pa、Pb對模型參數(shù)進(jìn)行調(diào)整得到

Dcom=D1×Pa+D2×Pb

(34)

Ucom=U1×Pa+U2×Pb

(35)

(36)

(37)

(38)

由于未測信號的先驗信息未知，因此未測信號與已測信號的互協(xié)方差也按照調(diào)整后的模型參數(shù)計算

(39)

(40)

在實際應(yīng)用中，Pa、Pb需要進(jìn)行迭代計算，迭代循環(huán)的截止條件則以每次計算獲得的單位權(quán)中誤差作為評判依據(jù)，當(dāng)前后兩次獲得的單位權(quán)中誤差之差小于某個限值時，循環(huán)結(jié)束，如圖1所示。

五、試驗分析

試驗區(qū)位于我國陸地某區(qū)域，該區(qū)域航空重力測量與陸地重力數(shù)據(jù)分布如圖2所示，整個區(qū)域有分辨率為5′×5′的航空擾動重力數(shù)據(jù)(飛行平均高度3600 m，測線數(shù)據(jù)采樣頻率為1 Hz)和2.5′分辨率的地面重力異常數(shù)據(jù)，地形數(shù)據(jù)采用SRTM 90 m分辨率的地形模型。為了檢驗一體化處理算法，本文假設(shè)參與地面重力異常數(shù)據(jù)ΔgT布設(shè)形式為均勻分布，如圖3所示，其中“★”代表參與一體化處理的地面重力數(shù)據(jù)，其余空格區(qū)域代表需要求解的重力數(shù)據(jù)，同時也用來對一體化處理結(jié)果進(jìn)行檢核。

首先利用航空擾動重力數(shù)據(jù)進(jìn)行格網(wǎng)化，然后利用直接代表法向下延拓[12]獲得地面待求點數(shù)據(jù)。其次利用地面重力數(shù)據(jù)按照經(jīng)典的最小二乘配置法推估該區(qū)域待求點重力數(shù)據(jù)，兩種方法分別與已知的地面重力數(shù)據(jù)進(jìn)行比較，見表1。

圖1　航空和地面重力數(shù)據(jù)一體化處理流程

圖2　航空重力測量與陸地重力數(shù)據(jù)分布

圖3　地面重力數(shù)據(jù)呈均勻分布

最后按照本文提出的一體化處理方法，利用航空擾動重力的測線數(shù)據(jù)求解地面待求點重力數(shù)據(jù)，統(tǒng)計結(jié)果見表2，最終處理后的地面重力異常數(shù)據(jù)如圖4所示。通過表1、表2對比可以看出，一體化處理方法比單獨使用航空和地面重力數(shù)據(jù)進(jìn)行推估的效果都要好，這也說明了一體化處理方法對于聯(lián)合處理航空和地面重力數(shù)據(jù)是非常適用的。

表1　航空和地面重力數(shù)據(jù)分別推估待求點數(shù)據(jù)統(tǒng)計結(jié)果　mGal

表2　利用一體化處理方法得到的結(jié)果　mGal

圖4　航空和地面重力數(shù)據(jù)一體化處理后的三維重力圖

六、結(jié)論

本文深入研究了航空擾動重力數(shù)據(jù)和地面重力異常數(shù)據(jù)的聯(lián)合處理問題，構(gòu)建了一種一體化處理方法，主要結(jié)論和創(chuàng)新點如下：

1) 在三維擾動位協(xié)方差模型基礎(chǔ)上，推導(dǎo)獲得局部重力異常協(xié)方差模型、局部擾動重力協(xié)方差模型，以及擾動重力和重力異常的互協(xié)方差模型。這3個模型在忽略小項后是一致的，即航空擾動重力數(shù)據(jù)與地面重力異常數(shù)據(jù)的聯(lián)合處理采用一個模型即可，大大簡化了處理流程，提高了計算效率。

2) 以最小二乘配置法為基礎(chǔ)，構(gòu)建了航空擾動重力數(shù)據(jù)與地面重力異常數(shù)據(jù)一體化處理的具體方法和實施步驟，利用該方法可推估地面任意待求點的重力數(shù)據(jù)，并可有效估計航空重力數(shù)據(jù)中的系統(tǒng)偏差。

3) 在某試驗區(qū)開展了一體化處理試驗，結(jié)果表明，采用一體化處理獲得的重力異常數(shù)據(jù)比單獨使用5′分辨率航空重力數(shù)據(jù)推估獲得的數(shù)據(jù)精度提高2.4 mGal，比單獨使用2.5′分辨率地面重力數(shù)據(jù)推估獲得的數(shù)據(jù)精度提高1.2 mGal，且一體化處理過程自動完成了航空擾動重力數(shù)據(jù)的向下延拓，同時也大大削弱了航空數(shù)據(jù)中存在的系統(tǒng)偏差。

本文所構(gòu)建的航空擾動重力數(shù)據(jù)和地面重力異常數(shù)據(jù)的一體化處理方法將為未來我國大規(guī)模的航空重力測量及磁力測量數(shù)據(jù)處理提供很好的參考和借鑒。

參考文獻(xiàn)：

[1]孫中苗.航空重力測量理論、方法及應(yīng)用研究[D].鄭州: 中國人民解放軍信息工程大學(xué),2004.

[2]OLESEN A V, ANDERSEN O B, TSCHERNING C C. Merging of Airborne Gravity and Gravity Derived from Satellite Altimetry: Test Cases Along the Coast of Greenland[J]. Studia Geophysica et Geodaetica, 2002, 46(3): 387-394.

[3]STRYKOWSKI G, FORSBERG R. Operational Merging of Satellite, Airborne and Surface Gravity Data by Draping Techniques[C]∥Geodesy on the Move: Gravity, Geoid, Geodynamics and Antarctica. Berlin: Springer Verlag, 1997: 243-248.

[4]BAYOUD F A. Some Investigations on Local Geoid Determination from Airborne Gravity Data[D]. Calgary, Canada:University of Calgary,2001:33-35.

[5]HWANG C, GUO J Y, DENG X L, et al. Coastal Gravity Anomalies from Retracked Geosat/GM Altimetry:Improvement, Limitation and the Role of Airborne Gravity Data[J]. Journal of Geodesy, 2006, 80(4): 204-216.

[6]郝燕玲，成怡，劉繁明, 等. 融合多類型海洋重力數(shù)據(jù)算法仿真研究[J]. 系統(tǒng)仿真學(xué)報, 2007, 19(21): 4897-4899.

[7]KERN M, SCHWARZ K K P P, SNEEUW N. A Study on the Combination of Satellite, Airborne, and Terrestrial Gravity Data[J]. Journal of Geodesy, 2003, 77(3-4): 217-225.

[8]楊元喜, 曾安敏. 大地測量數(shù)據(jù)融合模式及其分析[J]. 武漢大學(xué)學(xué)報: 信息科學(xué)版, 2008, 33(8): 771-774.

[9]楊元喜, 張菊清, 張亮. 基于方差分量估計的擬合推估及其在GIS誤差糾正的應(yīng)用[J]. 測繪學(xué)報, 2008, 37(2): 152-155.

[10]陸仲連. 地球重力場理論與方法[M]. 北京: 解放軍出版社, 1996: 316-320.

[11]鄒賢才.最小二乘配置在物理大地測量中的應(yīng)用[D]. 武漢: 武漢大學(xué),2003.

[12]王興濤,夏哲仁,石磐, 等.航空重力測量數(shù)據(jù)向下延拓方法比較[J].地球物理學(xué)報,2004,47(6):1017-1022.

本刊開通微信公眾平臺服務(wù)賬號

[本刊訊]為了更全面、及時地傳播測繪地理信息學(xué)科相關(guān)研究的最新成果，方便廣大讀者、作者獲取本刊信息，更大地發(fā)揮刊物的作用，本刊自2015年1月起正式開通微信公眾平臺服務(wù)賬號，定期發(fā)布本刊每期論文目錄、重點文章和新聞報道，歡迎廣大讀者、作者關(guān)注。

《測繪通報》微信號：CN11-2246

(本刊編輯部)

引文格式： 翟振和，孫中苗，李迎春，等. 航空和地面重力測量數(shù)據(jù)的一體化處理方法[J].測繪通報，2015(1)：39-43.DOI:10.13474/j.cnki.11-2246.2015.0007

作者簡介：翟振和(1980—)，男，博士，助理研究員，主要從事物理大地測量方面的研究。E-mail： zhaizhenhe1980@163.com

基金項目：國家自然科學(xué)基金(41174017)

收稿日期：2014-07-01

中圖分類號：P223

文獻(xiàn)標(biāo)識碼：B

文章編號：0494-0911(2015)01-0039-05