廣義殘差地形模型及其在局部重力場逼近中的應(yīng)用

吳懌昊，羅志才,2,3，周波陽

1. 武漢大學(xué)測繪學(xué)院，湖北 武漢 430079； 2. 武漢大學(xué)地球空間環(huán)境與大地測量教育部重點實驗室，湖北 武漢 430079； 3. 武漢大學(xué)測繪遙感信息工程國家重點實驗室，湖北 武漢 430079； 4. 廣東工業(yè)大學(xué)測繪工程系，廣東 廣州 510006

?

廣義殘差地形模型及其在局部重力場逼近中的應(yīng)用

吳懌昊1，羅志才1,2,3，周波陽4

1. 武漢大學(xué)測繪學(xué)院，湖北 武漢 430079； 2. 武漢大學(xué)地球空間環(huán)境與大地測量教育部重點實驗室，湖北 武漢 430079； 3. 武漢大學(xué)測繪遙感信息工程國家重點實驗室，湖北 武漢 430079； 4. 廣東工業(yè)大學(xué)測繪工程系，廣東 廣州 510006

Foundationsupport：TheNationalNaturalScienceFoundationofChina(Nos.41374023; 41131067; 41174020);TheStateScholarshipFundfromChineseScholarshipCouncil(No.201306270014);TheFundamentalResearchFundsfortheCentralUniversities(No.2042014gf041)

研究了殘差地形模型中的非調(diào)和性問題，比較了基于棱柱體和球冠體的積分模型，提出了基于球冠體積分的廣義殘差地形模型。以泊松小波徑向基函數(shù)為構(gòu)造基函數(shù)，結(jié)合廣義殘差地形模型，融合多源實測重力數(shù)據(jù)構(gòu)建了局部區(qū)域重力場模型。研究結(jié)果表明：基于棱柱體積分的殘差地形模型精度較低，在山區(qū)可能引入毫伽級以上的誤差，建議采用更為接近真實地形表面的球冠體積分模型。相比于原始的殘差地形模型，基于球冠體積分的廣義殘差地形模型能更為精確地逼近局部重力場模型中地形因素引起的高頻效應(yīng)。

殘差地形模型；調(diào)和改正；泊松小波徑向基函數(shù)；似大地水準(zhǔn)面；局部重力場建模

地形起伏反映了局部重力場的高頻擾動信息，精確確定地形影響對于局部重力場的精確逼近具有重要的意義[1-2]。殘差地形模型(RTM)可逼近地形擾動引起的高頻信號，不引入頻譜混疊效應(yīng)，在局部重力場建模中有廣泛的應(yīng)用[3-5]。文獻[6]采用殘差地形模型改善了重力場模型高頻部分的精度，削弱了截斷誤差的影響；文獻[7]聯(lián)合EGM2008重力場模型和殘差地形模型提高了山區(qū)大地水準(zhǔn)面的精度；文獻[8]聯(lián)合EGM2008和SRTM/DTM2006.0殘差地形模型改善了GPS高程轉(zhuǎn)換的精度。殘差地形模型僅考慮平滑面與地形面間的質(zhì)量，平滑面以上的多余質(zhì)量被移除，以下的缺失質(zhì)量被填充。經(jīng)殘差地形改正后，重力場逼近的邊界面變?yōu)樵撈交?。若重力點位于平滑面以下，則該點位于質(zhì)量體內(nèi)部，不滿足重力場逼近中的調(diào)和性條件，即為殘差地形模型中的非調(diào)和性問題[3]。

文獻[3]提出調(diào)和改正處理非調(diào)和性問題。將平滑面內(nèi)部的重力點與平滑面之間的質(zhì)量壓縮為位于該點以下一個無限薄的質(zhì)量層，使得測量點重新位于質(zhì)量體外部，滿足調(diào)和性條件?；谝迫?恢復(fù)法，文獻[3]提出的方法僅在移去階段對平滑面以下的重力點施加調(diào)和改正，在恢復(fù)階段則忽略了引入調(diào)和改正的影響。上述處理方法顯然值得討論，移去階段改變了重力場觀測參量，恢復(fù)階段也需考慮其相應(yīng)的影響。而如何合理地施加調(diào)和改正，綜合考慮調(diào)和改正在移去及恢復(fù)階段對外部重力場影響的研究尚未有成果發(fā)表。此外，文獻[9]從理論上比較基于棱柱體與球冠體積分的殘差地形模型的差別，認(rèn)為廣泛使用棱柱體模型精度較低，可能引入較大誤差，但并未通過實測多源重力數(shù)據(jù)進行數(shù)值分析。本文基于多源重力觀測數(shù)據(jù)，比較棱柱體和球冠體模型的差異，為積分體的選擇提供支持。針對非調(diào)和性問題，基于解析延拓的思想提出廣義殘差地形模型，結(jié)合基于徑向基函數(shù)構(gòu)建局部重力場模型的算例驗證該模型的有效性。

1　廣義殘差地形模型

1.1質(zhì)量體模型

一般而言，積分質(zhì)量體可為棱柱體、楔形體及球冠體等幾何體，只要基于該質(zhì)量體的積分公式具有解析形式即可[9]。使用最為廣泛的質(zhì)量體模型為棱柱體[10-12]，其引力位模型可以表示為

(1)

由于地球曲率的影響，即使采用高分辨率的數(shù)字地形模型，棱柱體之間也會存在空隙，基于棱柱體的積分質(zhì)量元逼近真實的地形模型只是較為粗略的近似，可能會影響到殘差地形模型逼近重力場高頻信號的精度[9]。與棱柱體相比，球冠體模型考慮了地球曲率的影響，能更為準(zhǔn)確地描述地形表面，且球冠體之間沒有空隙，精度較高，其引力位為

(2)

式中，λ1、λ2、φ1、φ2、r1、r2為劃分球冠體的經(jīng)圈、緯圈及地心半徑的積分上下限。

1.2廣義調(diào)和改正模型

重力測量點位于高程平滑面以下時，為滿足調(diào)和性條件，需對其施加調(diào)和改正。假設(shè)重力點與平滑面間的質(zhì)量層可視為一個內(nèi)半徑為R，厚度為Δh(外半徑為R+Δh)，密度為ρ的球?qū)?。此球?qū)釉谄鋬?nèi)部任意點x產(chǎn)生的引力位為

(3)

(4)

(5)

假設(shè)將位于平滑面內(nèi)部的重力點與平滑面間的質(zhì)量層(當(dāng)作球?qū)訉Υ?壓縮為位于重力點以下一無限薄的質(zhì)量層，此時μ=ρΔh。移動球?qū)訒r，對其內(nèi)部點的引力沒有影響。壓縮之后，由于薄層質(zhì)量的影響，引力增加了4πGρΔh，即為調(diào)和改正。對比式(3)和式(4)可知，重力點的引力位變化較小，文獻[3]認(rèn)為施加調(diào)和改正對高程異常的影響也可忽略，因此調(diào)和改正只在移去階段施加于的重力場觀測量，而忽略其在恢復(fù)階段對引力位或高程異常影響。實際情況中，重力點與平滑面之間的質(zhì)量層不能完全當(dāng)作球?qū)訉Υ?，在移去階段對重力觀測值施加調(diào)和改正，必然會對外部引力位產(chǎn)生影響，在恢復(fù)階段需考慮其影響。本文基于解析延拓的思想，將所有位于高程平滑面以下的重力測量點延拓到此平滑面上，綜合考慮調(diào)和改正在移去及恢復(fù)階段的影響，提出改進的殘差地形模型，將其稱為基于廣義調(diào)和改正的殘差地形模型，簡稱為廣義殘差地形模型(generalized residual terrain model)。

假設(shè)P點為殘差地形改正后位于地形質(zhì)量內(nèi)的重力點，Q點為相應(yīng)于P點位于高程平滑面上的點，根據(jù)泰勒展開式，忽略高階展開項的影響，Q點的殘余擾動位為

(6)

基于式(6)，利用移去-恢復(fù)法，Q點的殘余重力觀測值為

Δgres(Q)=Δg(P)-ΔgGM(P)-ΔgRTM(P)-

(7)

(8)

Q點殘余高程異常為

ζres(Q)=ζ(P)-ζGM(P)-ζRTM(P)-

(9)

2　泊松小波基函數(shù)建模

泊松小波徑向基函數(shù)在空間域和頻率域都有較好的局部化特性，可融合多源重力場數(shù)據(jù)，適合局部重力場建模[14-18]?；谝迫?恢復(fù)法，殘余擾動位Tres可表示為有限個泊松小波基函數(shù)之和[17]

(10)

多源重力場觀測數(shù)據(jù)可表示為擾動位的泛函，重力擾動δg、重力異常Δg及高程異常ζ在球面近似條件下分別與擾動位存在如下函數(shù)關(guān)系

(11)

結(jié)合式(10)，對于某一類觀測值可以建立觀測方程如下

(12)

式中，Lp表示第p類的重力場信息觀測值；Δp表示觀測誤差；fp表示此類觀測值與擾動位之間的泛函關(guān)系；J表示觀測數(shù)據(jù)種類的個數(shù)。將式(12)改寫為誤差方程的形式如下

Vp=ApX-lp

(13)

式中，Ap表示mp×K設(shè)計矩陣；X表示K×1基函數(shù)的未知參數(shù)向量；lp表示mp×1誤差方程的常數(shù)項向量；Vp表示mp×1觀測值殘差向量，mp表示該類重力場觀測值的個數(shù)?？偟恼`差方程為

V=AX-l

(14)

式中

(15)

假定不同類型的觀測值之間互不相關(guān)，觀測數(shù)據(jù)的方差-協(xié)方差陣為

(16)

利用最小二乘原理，基函數(shù)的未知參數(shù)向量X的估值可表示為

(17)

式中

(18)

(19)

采用方差分量估計的方法對各類觀測值進行合理地定權(quán)，各類觀測值的單位權(quán)方差因子估計

如下[19]

(20)

3　數(shù)值計算與分析

3.1采用數(shù)據(jù)

基于Eurodem、SRTM以及GEBCO3種數(shù)字地形模型構(gòu)建了局部區(qū)域高精度、高分辨率(2″×2″)陸海統(tǒng)一的數(shù)字地形模型，見圖1(a)。重力場模型采用DGM1S模型，其球諧展開階數(shù)達到250階[20]。采用移動平均法構(gòu)建殘差地形模型中的高程平滑面[10-11]，如圖1(b)所示，其空間分辨率與DGM1S相一致。收集了覆蓋整個荷蘭、比利時、英國及部分德國、法國、丹麥、挪威和北海區(qū)域的多源重力數(shù)據(jù)。通過交叉點平差完成了船載、航空重力數(shù)據(jù)中系統(tǒng)偏差的校正；利用閾值法和hampel濾波剔除了多源數(shù)據(jù)中的粗差；采用低通濾波削弱了船載、航空重力數(shù)據(jù)中的高頻噪聲的影響，并將各類重力數(shù)據(jù)歸算到同一參考框架(ETRS89)及垂直基準(zhǔn)(EVRF2007)。基于移去-恢復(fù)法，利用DGM1S模型和殘余地形模型移去重力場的長波和短波信號，計算的殘余重力觀測數(shù)據(jù)見圖2，其統(tǒng)計信息見表1。同時收集了上述國家部分區(qū)域的高精度的GPS水準(zhǔn)數(shù)據(jù)，并將其歸算到統(tǒng)一的參考框架及垂直基準(zhǔn)。

圖1　局部區(qū)域的數(shù)字地形模型(a)及高程平滑面(b)Fig.1　Digital terrain model (a) and mean elevation surface (b)

圖2　殘余重力觀測數(shù)據(jù)Fig.2　Residual gravity observations

表1　殘余重力觀測數(shù)據(jù)統(tǒng)計信息

3.2積分體模型的選擇

基于不同積分體模型計算的殘差地形改正的差異如圖3所示，結(jié)合表2可知，基于棱柱體與球冠體積分的殘差地形改正存在差異，其陸地、船載重力異常和航空重力擾動地形改正差異的標(biāo)準(zhǔn)差分別為1.156mGal、0.390mGal及0.078mGal。以陸地重力異常為例，兩者的差異主要集中在英國北部、挪威南部及德國和法國的部分多山區(qū)域，其量級達到毫伽級。與球冠體積分模型相比，棱柱體模型的求解精度相對較低，特別在地勢起伏較大的多山區(qū)域，忽略地球曲率的影響可能會引入較大的誤差。后續(xù)計算中采用基于球冠體積分的殘差地形模型用于高頻重力場信息的逼近。

表2　基于球冠體和棱柱體模型計算的殘差地形改正的差異的統(tǒng)計信息

3.3基于廣義殘差地形模型的局部重力場逼近

基于泊松小波基函數(shù)，結(jié)合廣義殘差地形模型，融合陸地、船載重力異常及航空重力擾動構(gòu)建局部重力場模型。試算區(qū)域覆蓋荷蘭、比利時、英國、愛爾蘭及部分挪威、德國、法國和北海區(qū)域，其緯度范圍為47°—63°，經(jīng)度范圍為-11°—13°，基于等間隔格網(wǎng)構(gòu)建基函數(shù)的網(wǎng)絡(luò)[17]，采用約為55 000個泊松小波基函數(shù)，其空間分辨率約為7.1km。圖4(a)表示基于廣義殘差地形模型計算的重力觀測數(shù)據(jù)的調(diào)和改正，其信號主要集中在挪威南部、英國北部及德國西南部的山脈區(qū)域。在地形起伏劇烈的多山區(qū)域，反映地形中長波段信息的高程平滑面與真實地形起伏差距較大。若重力觀測點位于高程平滑面以下，需施加的調(diào)和改正較大，局部范圍內(nèi)可達數(shù)十毫伽。圖4(b)顯示了調(diào)和改正在恢復(fù)階段對高程異常的影響。在地形較為平坦的地區(qū)，其影響僅達到毫米量級。而在地勢起伏較大的英國北部及挪威南部山區(qū)，調(diào)和改正的影響可達厘米甚至分米級，忽略其影響會引入較大的誤差。表3、表4分別給出了基于廣義殘差地形模型和原始模型構(gòu)建的似大地水準(zhǔn)面的檢核結(jié)果。基于兩種模型構(gòu)建的似大地水準(zhǔn)面在地形起伏較為平緩的荷蘭、比利時及德國區(qū)域精度相當(dāng)。而在地形起伏較大挪威南部和英國北部山區(qū)，結(jié)合廣義殘差地形模型構(gòu)建的似大地水準(zhǔn)面的精度分別提高了約2.1cm和1.7mm。

圖3　基于球冠體和棱柱體積分模型計算的殘差地形改正的差異Fig.3　The difference of the RTM correction computed from prism and tesseroid model

圖4　基于廣義殘差地形模型解算的調(diào)和改正對重力觀測數(shù)據(jù)和高程異常的影響Fig.4　The effect of harmonic corrections on gravity observations and height anomaly computed from the generalized RTM

表3　基于廣義殘差地形模型構(gòu)建的局部似大地水準(zhǔn)面的精度評價

4　結(jié)　論

本文研究了殘差地形模型中的非調(diào)和性問題，比較了基于棱柱體和球冠體的積分模型，并基于解析延拓的思想提出了基于球冠體積分的廣義殘差地形模型。利用徑向基函數(shù)方法，結(jié)合廣義殘差地形模型，融合實測的陸地、船載重力異常及航空重力擾動數(shù)據(jù)構(gòu)建了局部區(qū)域陸海統(tǒng)一的似大地水準(zhǔn)面模型。結(jié)果表明：與球冠體積分模型相比，棱柱體模型的求解精度相對較低，特別在地勢起伏較大的多山區(qū)域，忽略地球曲率的影響可能會引入毫伽級以上的誤差。在地形較為平坦的區(qū)域，調(diào)和改正對高程異常產(chǎn)生的影響僅達到毫米量級，在厘米級似大地水準(zhǔn)面的構(gòu)建中可忽略其影響。而在地形起伏較大的山區(qū)，忽略調(diào)和改正的影響可能引入厘米甚至分米級的誤差。廣義殘差地形模型綜合考慮了調(diào)和改正在移去及恢復(fù)階段的影響，能更為精確地逼近地形因素引起的高頻效應(yīng)。與原始的殘差地形模型相比，結(jié)合廣義殘差地形模型構(gòu)建的似大地水準(zhǔn)面的精度在地形起伏較大的英國及挪威山區(qū)分別提高了1.7mm和2.1cm。

表4　基于原始殘差地形模型構(gòu)建的局部似大地水準(zhǔn)面的精度評價

[1]羅志才, 陳永奇, 寧津生. 地形對確定高精度局部大地水準(zhǔn)面的影響[J]. 武漢大學(xué)學(xué)報(信息科學(xué)版), 2003, 28(3): 340-344.

LUOZhicai,CHENYongqi,NINGJinsheng.EffectofTerrainontheDeterminationofHighPreciseLocalGravimetricGeoid[J].GeomaticsandInformationScienceofWuhanUniversity, 2003, 28(3): 340-344.

[2]章傳銀, 晁定波, 丁劍, 等. 厘米級高程異常地形影響的算法及特征分析[J]. 測繪學(xué)報, 2006, 35(4): 308-314.ZHANGChuanyin,CHAODingbo,DINGJian,etal.ArithmeticandCharactersAnalysisofTerrainEffectsforcmOrderPrecisionHeightAnomaly[J].ActaGeodaeticaetCartographicaSinica, 2006, 35(4): 308-314.

[3]FORSBERGR,TSCHERNINGCC.TheUseofHeightDatainGravityFieldApproximationbyCollocation[J].JournalofGeophysicalResearch, 1981, 86(B9): 7843-7854.

[4]OMANGOCD,FORSBERGR.HowtoHandleTopographyinPracticalGeoidDetermination:ThreeExamples[J].JournalofGeodesy, 2000, 74(6): 458-466.

[5]HIRTC,MARTIU,BüRKIB,etal.AssessmentofEGM2008inEuropeUsingAccurateAstrogeodeticVerticalDeflectionsandOmissionErrorEstimatesfromSRTM/DTM2006.0ResidualTerrainModelData[J].JournalofGeophysicalResearch, 2010, 115(B10): 196-211.

[6]HIRTC.RTMGravityForward-modelingUsingTopography/BathymetryDatatoImproveHigh-degreeGlobalGeopotentialModelsintheCoastalZone[J].MarineGeodesy, 2013, 36(2): 183-202.

[7]HIRTC,FEATHERSTONEWE,MARTIU.CombiningEGM2008andSRTM/DTM2006.0ResidualTerrainModelDatatoImproveQuasigeoidComputationsinMountainousAreasDevoidofGravityData[J].JournalofGeodesy, 2010, 84(9): 557-567.

[8]張興福, 劉成. 綜合EGM2008模型和SRTM/DTM2006.0剩余地形模型的GPS高程轉(zhuǎn)換方法[J]. 測繪學(xué)報, 2012, 41(1): 25-32.

ZHANGXingfu,LIUCheng.TheApproachofGPSHeightTransformationBasedonEGM2008andSRTM/DTM2006.0ResidualTerrainModel[J].ActaGeodaeticaetCartographicaSinica, 2012, 41(1): 25-32.

[9]HECKB,SEITZK.AComparisonofTheTesseroid,PrismandPoint-massApproachesforMassReductionsinGravityFieldModelling[J].JournalofGeodesy, 2007, 81(2): 121-136.

[10]HIRTC.PredictionofVerticalDeflectionsfromHigh-degreeSphericalHarmonicSynthesisandResidualTerrainModelData[J].JournalofGeodesy, 2010, 84(3): 179-190.

[11]TZIAVOSIN,VERGOSGS,GRIGORIADISVN.InvestigationofTopographicReductionsandAliasingEffectsonGravityandtheGeoidoverGreeceBasedonVariousDigitalTerrainModels[J].SurveysinGeophysics, 2010, 31(1): 23-67.

[12]NAGYD,PAPPG,BENEDEKJ.TheGravitationalPotentialandItsDerivativesforthePrism[J].JournalofGeodesy, 2000, 74(7): 552-560.

[13]寧津生, 李建成, 晁定波. 由地面重力數(shù)據(jù)確定擾動位徑向二階梯度[J]. 武漢測繪科技大學(xué)學(xué)報, 1996, 21(3): 201-206.NINGJinsheng,LIJiancheng,CHAODingbo.Determinationofthe2ndRadialGradientofDisturbingPotentialUsingTerrestrialGravityData[J].JournalofWuhanTechnicalUniversityofSurveyingandMapping, 1996, 21(3): 201-206.

[14]ANTUNESC,PAILR,CATALOJ.PointMassMethodAppliedtotheRegionalGravimetricDeterminationoftheGeoid[J].StudiaGeophysicaetGeodaetica, 2003, 47(3): 495-509.

[15]BENTELK,SCHMIDTM,GERLACHC.DifferentRadialBasisFunctionsandTheirApplicabilityforRegionalGravityFieldRepresentationontheSphere[J].GEM-InternationalJournalonGeomathematics, 2013, 4(1): 67-96.

[16]BENTELK,SCHMIDTM,DENBYCR.ArtifactsinRegionalGravityRepresentationswithSphericalRadialBasisFunctions[J].JournalofGeodeticScience, 2013, 3(3): 173-187.

[17]KLEESR,TENZERR,PRUTKINI,etal.AData-drivenApproachtoLocalGravityFieldModellingUsingSphericalRadialBasisFunctions[J].JournalofGeodesy, 2008, 82(8): 457-471.

[18]TENZERR,KLEESR.TheChoiceoftheSphericalRadialBasisFunctionsinLocalGravityFieldModeling[J].StudiaGeophysicaetGeodaetica, 2008, 52(3): 287-304.

[19]KUSCHEJ.NoiseVarianceEstimationandOptimalWeightDeterminationforGOCEGravityRecovery[J].AdvancesinGeosciences, 2003, 1: 81-85.

[20]FARAHANIHH,DITMARP,KLEESR,etal.TheStaticGravityFieldModelDGM-1SfromGRACEandGOCEData:Computation,ValidationandanAnalysisofGOCEMission’sAddedValue[J].JournalofGeodesy, 2013, 87(9): 843-867.

(責(zé)任編輯：宋啟凡)

GeneralizedResidualTerrainModelandItsApplicationinRegionalGravityFieldApproximation

WUYihao1,LUOZhicai1,2,3,ZHOUBoyang4

1.SchoolofGeodesyandGeomatics,WuhanUniversity,Wuhan430079,China; 2.KeyLaboratoryofGeospaceEnvironmentandGeodesy,MinistryofEducation,WuhanUniversity,Wuhan430079,China; 3.StateKeyLaboratoryofInformationEngineeringinSurveying,MappingandRemoteSensing,WuhanUniversity,Wuhan430079,China; 4.DepartmentofSurveyingandMapping,GuangdongUniversityofTechnology,Guangzhou510006,China

Westudiedtheso-callednon-harmonicprobleminresidualterrainmodel(RTM)andcomparedtheRTMcorrectionsbasedonprismsandtesseroids.Moreover,weproposedthegeneralizedRTMbasedontesseroids,togetherwithwhichtheregionalgravityfieldwasmodeledbasedonPoissonwaveletsbasisfunctionbyusingheterogeneousgravitydatasets.TheresultsshowthattheRTMcorrectionbasedonprismintegralhasapoorperformanceinmountainousregions,whichmayintroduceerrorswiththemagnitudeofseveralmGal.Thus,wesuggestusingRTMbasedontesseroids,whichleadtoabetterapproximationofthetopography.Comparedtooriginalresidualterrainmodel，thegeneralizedoneleadsabetterapproximationoftheregionalgravityfiledatthehigh-frequencypartcausedbylocaltopographicalvariation.

residualterrainmodel;harmoniccorrection;Poissonwaveletsradialbasisfunction;quasi-geoid;regionalgravityfieldmodeling

WUYihao,LUOZhicai,ZHOUBoyang.GeneralizedResidualTerrainModelandItsApplicationinRegionalGravityFieldApproximation[J].ActaGeodaeticaetCartographicaSinica,2016,45(5):531-537.DOI:10.11947/j.AGCS.2016.20150219.

P223

A

1001-1595(2016)05-0531-07

國家自然科學(xué)基金(41374023; 41131067; 41174020)；國家留學(xué)基金委基金(201306270014)；中央高?；究蒲袠I(yè)務(wù)費專項(2042014gf041)

引文格式：吳懌昊，羅志才，周波陽.廣義殘差地形模型及其在局部重力場逼近中的應(yīng)用[J].測繪學(xué)報，2016,45(5)：531-537.DOI:10.11947/j.AGCS.2016.20150219.