李姍姍,吳曉平,張傳定,歐陽永忠
1.信息工程大學(xué)測繪學(xué)院,河南鄭州450052;2.海軍海洋測繪研究所,天津300061
顧及地形與完全球面布格異常梯度項改正的區(qū)域似大地水準(zhǔn)面精化
李姍姍1,吳曉平1,張傳定1,歐陽永忠2
1.信息工程大學(xué)測繪學(xué)院,河南鄭州450052;2.海軍海洋測繪研究所,天津300061
推證顧及地形與完全球面布格異常梯度改正的完全到一階項的物理大地測量邊值問題的嚴(yán)密解式,并在某試驗區(qū)綜合利用地形、重力、GPS/水準(zhǔn)等數(shù)據(jù)進(jìn)行區(qū)域似大地水準(zhǔn)面的計算與檢驗。通過對高程異常計算絕對與相對精度的比較分析,結(jié)果表明,完全球面布格異常梯度改正項對高程異常的影響能夠達(dá)到厘米的量級。因此,提高區(qū)域似大地水準(zhǔn)面的建模精度,尤其是在地形起伏較大的區(qū)域,除需顧及地形改正項影響外,還應(yīng)考慮完全球面布格異常梯度改正項對高程異常的影響。
似大地水準(zhǔn)面;地形改正;完全球面布格異常梯度;奇異積分;高程異常
似大地水準(zhǔn)面是我國高程系統(tǒng)的一個基準(zhǔn)面,目前可提供使用的2000中國重力場與似大地水準(zhǔn)面模型(CGGM2000),其絕對精度在我國東部地區(qū)優(yōu)于±0.2m,西部地區(qū)優(yōu)于±0.3m,并保持了陸地和近海區(qū)域重力場與大地水準(zhǔn)面的整體性。然而隨著衛(wèi)星、航空、海洋、陸地重力測量技術(shù)的發(fā)展與成熟,地面重力數(shù)據(jù)覆蓋的密度與均勻度越來越高,因此為了進(jìn)一步提高區(qū)域大地水準(zhǔn)面建模的絕對與相對精度,以滿足工程上的應(yīng)用需求,有必要對基于Molodensky理論求解地球重力場的嚴(yán)密理論、方法以及實現(xiàn)進(jìn)行深入研究,從而才有可能實現(xiàn)“GPS+似大地水準(zhǔn)面模型”技術(shù),快速為用戶提供厘米級精度的正常高。
似大地水準(zhǔn)面精化問題一直以來都是我國物理大地測量學(xué)者致力研究的熱點。文獻(xiàn)[1—3]利用地球重力場“可疊加性”的特性以及我國GPS/水準(zhǔn)、重力、地形以及全球重力場模型,采用均衡理論進(jìn)行重力歸算與內(nèi)插外推,采用球面一維FFT計算技術(shù)、基于Stokes和顧及G1項的Molodenskey公式,并顧及橢球改正,建立了我國分米級的(似)大地水準(zhǔn)面模型;文獻(xiàn)[4]利用高精度重力、地形、高階地球重力場模型和GPS/水準(zhǔn)數(shù)據(jù),采用嚴(yán)密的球面核函數(shù)和1D-FFT技術(shù),建立了我國深圳市1km分辨率的厘米級大地水準(zhǔn)面模型;文獻(xiàn)[5]研究探討了在已知厘米級精度地球重力場位系數(shù)模型、地形密度分布模型以及平均海面高模型的前提下,基于虛擬壓縮恢復(fù)理論實現(xiàn)全球1°×1°厘米級精度大地水準(zhǔn)面的方法;文獻(xiàn)[6]研究分析了基于二維平面和二維球面FFT計算大地水準(zhǔn)面算法的特點以及差異,并利用改進(jìn)后的數(shù)學(xué)模型對我國海陸大地水準(zhǔn)面進(jìn)行了計算與精度評定;文獻(xiàn)[7]提出了為滿足高精度局域大地水準(zhǔn)面的建模要求,在不同地形區(qū)域?qū)PS水準(zhǔn)網(wǎng)精度與密度的布設(shè)要求。
在實際基于高分辨率、高精度的重力數(shù)據(jù)、高分辨率數(shù)字地面模型(digital terrain model,DTM)以及全球重力場模型計算高程異常時,一般采用移去-恢復(fù)技術(shù)[8],且認(rèn)為Molodenskey一階解等價于地形改正項影響,假設(shè)的前提是重力異常與地形成線性關(guān)系[9-10]。然而在地形起伏較大的區(qū)域,重力異常與地形之間并不嚴(yán)格成線性關(guān)系。針對這一問題,本文在上述研究的基礎(chǔ)上,在顧及一階項的Molodenskey理論計算重力似大地水準(zhǔn)面方面作相應(yīng)研究,推導(dǎo)建立其嚴(yán)密解式,并利用某試驗區(qū)數(shù)據(jù)進(jìn)行計算分析與精度評定。
球近似意義下,Δg與外部擾動位T的關(guān)系為[11]
式中,Σ為地球自然表面,其上的重力異常值為Δg。若以似大地水準(zhǔn)面(或大地水準(zhǔn)面)為界將擾動位分為兩個部分,分別以Ti表示其內(nèi)部擾動質(zhì)量的位,Te表示外部地形質(zhì)量的位。相應(yīng)的,Σ面上的重力異常也可分解為兩個部分Δgi和Δge,即
對于外部地形物質(zhì)位,可表示為
若只取通常的球近似,也就是說,用一個半徑為R的球來代替似大地水準(zhǔn)面,則上式積分區(qū)域為似大地水準(zhǔn)面之外(R≤r≤R+h);dτ是體積單元,設(shè)點P是計算地形位的點,則l是dτ至P點的距離;ρ是地形質(zhì)量的密度,式中已假定為常數(shù)2.67gcm-3,略去h/R≤0.14%的影響,則地形物質(zhì)的位就變?yōu)?/p>
式中,對于dσ的積分,表示對整個立體角積分,而
它表示積分dτ的海拔高(球近似意義的海面也是指r=R的球面)。將式(5)分解為兩個部分,完成積分可將地形質(zhì)量的位表示為
地形物質(zhì)的垂直引力是地形物質(zhì)位的徑向方向?qū)?shù)的負(fù)值[12],即
又
代入式(8),得
并有
與式(5)比較,也可以將其分解為相應(yīng)的兩項,分別完成積分后,得
式中,C為球面近似地形改正[13-15]
將式(12)代入式(10),得
式(14)代入式(3),得
由式(2)和式(15),得
為完全球面布格異常。由于完全球面布格異常不僅扣除了地形物質(zhì)引力的影響,而且也扣除了地形位的影響;與平面布格異常相比,進(jìn)一步消除了地形有關(guān)的高頻部分[16],因而其變化更加平緩,值得在實踐中采用。
剩下的問題是如何確定計算點處的(Ti)A。對于線性地面邊值問題,一般都是作為固定邊值問題處理,即過計算點A的重力位水準(zhǔn)面W=WA是一個物理曲面,同時也可作為一個幾何曲面來看待。在減去地形物質(zhì)產(chǎn)生的重力異常Δge之后,則可視似大地水準(zhǔn)面外已沒有質(zhì)量存在。此時Ti在似大地水準(zhǔn)面外是調(diào)和的,并且該擾動位對應(yīng)的重力異常在Σ面上的值Δgi也是已知的。這樣的邊值問題可解釋為重力在Σ面上進(jìn)行測量,確定Σ上的擾動位Ti的問題。其解可依照Moritz的解析延拓理論[8],將Σ面上P點處的(Δgi)P解析延拓到移動地形物質(zhì)前的過A點的水準(zhǔn)面W=WA上(這時可視為幾何曲面,是固定界面)的P0點,則
式中
擾動位Ti在A點處的值為
式中
將式(7)和式(19)相加直接寫出聯(lián)合地形和地面重力異常數(shù)據(jù)確定地面擾動位的基本公式
式中
式(22)這一解式,各項物理意義明晰,但是每個量的計算都是一個全球積分,在進(jìn)行實際數(shù)據(jù)處理時,除了地形高數(shù)據(jù)和重力異常數(shù)據(jù)密度不匹配會引進(jìn)很大的附加誤差外,還存在布格改正直接、間接影響量級很大,相互抵償后也會引進(jìn)計算誤差。為此,利用球近似下恒等式[17]
取χ0=fρh,并顧及式(7)、(15)和(23),則可得到
利用式(25)的第2個等式,可將式(22)簡化為
至此,可以看出,基于球面布格異常所得球面地形改正解的主項是前述空間異常與高程成線性關(guān)系時的地形改正解,再加上空間異常與高程異常非線性部分的改正項和所有地面邊值問題解中都要涉及的高階改正項,這里為δTe+δTi。
最后略去式(26)的高階項,由布隆斯公式得到聯(lián)合地形和重力異常數(shù)據(jù)確定似大地水準(zhǔn)面的公式
在實際計算應(yīng)用過程中,一般將積分區(qū)域分為內(nèi)區(qū)σ0和外區(qū)σ-σ0影響,并視內(nèi)區(qū)近似為平面,因此式(20)可表示為
為解決內(nèi)區(qū)中的奇異積分,在計算點附近中心引入重力異常雙二次曲面擬合模型[18-20]
則中心區(qū)非奇異計算公式為
由于外區(qū)影響與距離的立方成反比,因此外區(qū)積分只需取至一定的積分半徑區(qū)域范圍,此時積分變?yōu)槔奂忧蠛偷男问?/p>
選定某試驗區(qū),該區(qū)域重力異常與高程數(shù)值變化大小如圖1、2所示,數(shù)據(jù)分辨率均為1′×1′,按照上述公式計算得到了該試驗區(qū)完全球面布格異常垂向梯度,結(jié)果如圖3所示,并計算了該區(qū)域重力異常的垂向梯度,如圖4,兩者數(shù)值統(tǒng)計比較結(jié)果如表1所示。
圖1 試驗區(qū)重力異常Fig.1 Gravity anomalies in the experiment area
圖2 試驗區(qū)高程Fig.2 Heights in the experiment area
圖3 完全球面布格異常垂向梯度Fig.3 Vertical gradient of complete spherical Bouguer anomalies
圖4 重力異常垂向梯度Fig.4 Vertical gradient of gravity anomalies
表1 重力異常與完全球面布格異常的垂向梯度的統(tǒng)計Tab.1 Statistics of vertical gradients of gravity anomalies and complete spherical Bouguer anomalies E
從圖3、圖4可以看出,完全球面布格異常垂向梯度L(ΔgB)較之于重力異常垂向梯度L(Δg)的變化要平緩得多,表1的統(tǒng)計結(jié)果也反映出重力異常垂向梯度L(Δg)比完全球面布格異常垂向梯度L(ΔgB)的變化幅度大一個數(shù)量級。顯然,這是由于重力異常與地形的相關(guān)性較強,而扣除了地形相關(guān)部分所得的完全球面布格異常的垂向梯度的量級較小,便于實際數(shù)值的內(nèi)插、外推以及歸算等應(yīng)用計算。進(jìn)一步將完全球面布格異常及其垂向梯度與完全平面布格異常及垂向梯度進(jìn)行比較,結(jié)果如圖5、圖6所示。
圖5 兩者異常差值Fig.5 Difference between the two types of Bouguer anomalies
圖6 兩者垂向梯度差值Fig.6 Vertical gradient difference between the two types of Bouguer anomalies
表2 兩者及其垂向梯度差值的統(tǒng)計Tab.2 Statistics of difference and vertical gradient difference between the two types of Bouguer anomalies
從表2兩者差值的統(tǒng)計結(jié)果可以看出,完全球面布格異常與完全平面布格異常相比,量級上的變化大概在幾個毫伽左右,垂向梯度數(shù)值上的變化并不明顯,用完全平面布格異常梯度替代完全球面布格異常梯度項對大地水準(zhǔn)面所產(chǎn)生的影響在目前測量水平的基礎(chǔ)上可以忽略。然而實現(xiàn)完全球面布格異常以及梯度項計算,一是從理論上而言更趨于嚴(yán)密,二是隨著計算機(jī)性能水平的提高,計算速度已不是問題,尋求更準(zhǔn)確、更精細(xì)地確定地球重力場元將是未來的發(fā)展趨勢。
地形改正C在線性邊值問解中的作用已被認(rèn)可,本節(jié)主要探討完全球面布格異常梯度改正g1項(g1=-(h-h(huán)A)L(ΔgB))對高程異常的影響,并與地形改正C的作用相比較。該試驗區(qū)地形改正C和g1項對高程異常影響的計算結(jié)果如圖7、圖8所示,統(tǒng)計結(jié)果如表3所示。
圖7 地形改正對高程異常的影響Fig.7 Effect on height anomaly by terrain correction
圖8 g1項對高程異常的影響Fig.8 Effect on height anomaly by g1term
表3 地形改正與g1項對高程異常的影響Tab.3 Effect on height anomaly by terrain correction and g1term cm
從表3中可以看出,g1項對高程異常的影響也可達(dá)到幾個厘米。這些數(shù)值結(jié)果也說明邊值問題的地形改正項與完全球面布格異常所引起的梯度改正項在某些區(qū)域能達(dá)到同等的量級。若同時顧及這兩項改正,可達(dá)到取至一階項解析延拓解的精度。
進(jìn)一步探討該試驗區(qū)完全球面布格異常與完全平面布格異常梯度改正項對高程異常計算結(jié)果的影響,從圖9和表4可以看出兩者對高程異常計算結(jié)果影響的量級基本相當(dāng),但實踐中編程計算速度并無明顯區(qū)別,因此從理論的嚴(yán)密程度來看,建議在計算高程異常時采用完全球面布格異常梯度改正項。
圖9 完全球面與完全平面布格異常梯度改正項對高程異常計算結(jié)果影響的差值Fig.9 Difference between the height anomaly calculation results considering complete spherical and complete plane Bouguer anomaly’s gradient terms separately
表4 完全球面布格異常梯度改正項與完全平面布格異常梯度改正項對高程異常計算結(jié)果的差值統(tǒng)計Tab.4 Statistics of difference between the height anomaly calculation results considering complete spherical and complete plane Bouguer anomaly’s gradient terms separately cm
采用EGM2008地球重力場模型,基于移去-恢復(fù)技術(shù),計算了圖1試驗區(qū)1′×1′的重力高程異常值(未加完全球面布格異常梯度改正項計算的格網(wǎng)高程異常值,加了完全球面布格異常梯度改正項計算的格網(wǎng)高程異常值)。為了評定兩者的計算精度,利用該試驗區(qū)范圍內(nèi)的20個GPS水準(zhǔn)點作為真值,按插值方法由格網(wǎng)高程異常值分別推求出20個點的高程異常計算值,二者之差作為真誤差,得到它們的中誤差如表5所示。
表5 重力高程異常的計算精度Tab.5 Calculation precision of gravity height anomaly m
從表5中誤差的數(shù)值可以看出,顧及完全球面布格異常梯度改正項,能提高高程異常的計算精度。此外統(tǒng)計40km左右區(qū)域范圍內(nèi)的相對高程異常差,以GPS水準(zhǔn)兩點之間的高程異常差作為真值,并分別計算加入g1項與不加入g1項的高程異常差,與GPS水準(zhǔn)高程異常差求差進(jìn)行精度評定,統(tǒng)計結(jié)果如表6所示。
表6 重力高程異常差的計算精度Tab.6 Calculation precision of gravity height anomaly difference m
顯然顧及地形改正及其完全球面布格異常梯度改正項的影響,高程異常相對計算精度能夠達(dá)到厘米級水平。
區(qū)域厘米級似大地水準(zhǔn)面的建立是物理大地測量的重要工作之一。隨著測量技術(shù)的日益發(fā)展,重力及地形數(shù)據(jù)的覆蓋日趨密集與均勻,因而對重力場數(shù)據(jù)處理方法也要求盡可能嚴(yán)密,并能得以工程化實現(xiàn)。本文推導(dǎo)了完全到一階項的聯(lián)合重力和地形數(shù)據(jù)確定高程異常的解式,并在某試驗區(qū)進(jìn)行了計算比較,數(shù)值結(jié)果表明:① 完全球面布格異常垂向梯度較之于空間異常垂向梯度的數(shù)值量級要小得多,說明了完全球面布格異常變化平緩,更適宜于實踐中重力場的內(nèi)插、外推以及重力數(shù)據(jù)歸算等工作;② 顧及地形改正及其完全球面布格異常梯度改正項的影響,與僅考慮地形改正項相比能提高高程異常建模的絕對精度,同時相對高程異常的計算精度能夠達(dá)到厘米級水平。
[1] CHEN Junyong,LI Jiancheng,NING Jinsheng,et al.On a High Resolution and High Accuracy Geoid in China Mainland[J].Acta Geodaetica et Cartographica Sinica,2001,30(2):95-100.(陳俊勇,李建成,寧津生,等.我國大陸高精度、高分辨率大地水準(zhǔn)面的研究和實施[J].測繪學(xué)報,2001,30(2):95-100.)
[2] CHEN Junyong,LI Jiancheng,NING Jinsheng,et al.On a Chinese New Quasi Geoid[J].Acta Geodaetica et Cartographica Sinica,2002,31(Sup.):1-6.(陳俊勇,李建成,寧津生,等.中國似大地水準(zhǔn)面[J].測繪學(xué)報,2002,31(增刊):1-6.)
[3] CHEN Junyong,LI Jiancheng,NING Jinsheng,et al.A New Chinese Geoid with High Resolution and High Accuracy[J].Geomatics and Information Science of Wuhan University,2001,26(4):283-289.(陳俊勇,李建成,寧津生,等.中國新一代高精度、高分辨率大地水準(zhǔn)面的研究和實施[J].武漢大學(xué)學(xué)報:信息科學(xué)版,2001,26(4):283-289.)
[4] NING Jinsheng,LUO Zhicai,YANG Zhanji,et al.Determination of Shenzhen Geoid with 1km Resolution and Centimeter Accuracy[J].Acta Geodaetica et Cartographica Sinica,2003,32(2):102-107.(寧津生,羅志才,楊沾吉,等.深圳市1km高分辨率厘米級高精度大地水準(zhǔn)面的確定[J].測繪學(xué)報,2003,32(2):102-107.)
[5] CHAO Dingbo,SHEN Wenbin,WANG Zhengtao.Investigations of the Possibility and Method of Determining Global Centimeter-level Geoid[J].Acta Geodaetica et Cartographica Sinica,2007,36(4):370-376.(晁定波,申文斌,王正濤.確定全球厘米級精度大地水準(zhǔn)面的可能性和方法探討[J].測繪學(xué)報,2007,36(4):370-376.)
[6] HUANG Motao,ZHAI Guojun,GUAN Zheng,et al.Comments on FFT Technique in Spectral Geoid Determination[J].Acta Geodaetica et Cartographica Sinica,2000,29(2):124-131.(黃謨濤,翟國君,管錚,等.利用FFT技術(shù)計算大地水準(zhǔn)面高若干問題研究[J].測繪學(xué)報,2000,29(2):124-131.)
[7] CHEN Junyong.Needs for GPS Leveling Network and Gravimetry in a High Accuracy Local Area Geoid[J].Acta Geodaetica et Cartographica Sinica,2001,30(3):189-191.(陳俊勇.高精度局域大地水準(zhǔn)面對布測GPS水準(zhǔn)和重力的要求[J].測繪學(xué)報,2001,30(3):189-191.)
[8] HOFMANN-WELLENHOF B,MORITZ H.Physical Geodesy[M].Graz,Wien,New York:Springer,2005.
[9] GUAN Zelin,GUAN Zheng,HUANG Motao,et al.Theory and Method of Reginonal Gravity Field Approximation[M].Beijing:Surveying and Mapping Press,1997.(管澤霖,管錚,黃謨濤,等.局部重力場逼近理論和方法[M].北京:測繪出版社,1997.)
[10] GUO Junyi.Foundation of Physical Geodesy[M].Wuhan:Wuhan Surveying and Mapping Techonology University Press,1994.(郭俊義.物理大地測量學(xué)基礎(chǔ)[M].武漢:武漢測繪科技大學(xué)出版社,1994.)
[11] LU Zhonglian.Theory and Method of the Earth’s Gravity Field[M].Beijing:PLA Publishing House,1996.(陸仲連.地球重力場的理論與方法[M].北京:解放軍出版社,1996.)
[12] EMELJANOV N V,KANTER A A.A Method to Compute Inclination Functions and Their Derivatives[J].Manuscripta Geodaetica,1989,14(2):77-83.
[13] FORSBERG R.Gravity Field Terrain Effect Computations by FFT[J].Journal of Geodesy,1984,59(4):342-360.
[14] SCHWARZ K P,SIDERIS M G,F(xiàn)ORSBERG R.The Use of FFT Techniques in Physical Geodesy[J].Geophysical Journal International,1990,100(3):485-514.
[15] SIDERIS M G.A Fast Fourier Transform Method of Computing Terrain Corrections[J].Manuscript Geodaetica,1985,10(1):66-73.
[16] LI Shanshan,WU Xiaoping,ZHANG Chuanding,et al.Calculation and Analysis of the New Statistical Character Parameters of Gravity Field in China[J].Chinese Journal of Geophysics,2010,53(5):1099-1108.(李姍姍,吳曉平,張傳定,等.我國重力場新的統(tǒng)計特征參數(shù)的計算分析[J].地球物理學(xué)報,2010,53(5):1099-1108.)
[17] LI Jiancheng,CHEN Junyong,NING Jinsheng,et al.Approximation Theory of the Earth Gravity Field and Determination of the Chinese Gravity Geoid Model 2000[M].Wuhan:Wuhan University Press,2003.(李建成,陳俊勇,寧津生,等.地球重力場逼近理論與中國2000似大地水準(zhǔn)面的確定[M].武漢:武漢大學(xué)出版社,2003.)
[18] ZHU Changqing.Method of Calculation and Application in Surveying and Mapping[M].Beijing:Surveying and Mapping Press,1988.(朱長青.計算方法及其在測繪中的應(yīng)用[M].北京:測繪出版社,1988.)
[19] WU Zongmin.Model,Method and Theory of Scattering Data Fitting[M].Beijing:Science Press,2008.(吳宗敏.散亂數(shù)據(jù)擬合的模型、方法和理論[M].北京:科學(xué)出版社,2008.)
[20] LI Shanshan,WU Xiaoping,ZHAO Dongming.Coons Curved Surface Reconstruction Method of Marine Gravity Anomaly Map for Navigation[J].Acta Geodaetica et Cartographica Sinica,2010,39(5):508-515.(李姍姍,吳曉平,趙東明.導(dǎo)航用海洋重力異常圖的孔斯曲面重構(gòu)方法[J].測繪學(xué)報,2010,39(5):508-515.)
E-mail:zzy_lily@sina.com
Regional Quasi-Geoid Refining Considering Corrections of Terrain and Complete Spherical Bouguer Anomaly’s Gradient Term
LI Shanshan1,WU Xiaoping1,ZHANG Chuanding1,OUYANG Yongzhong2
1.Institute of Surveying and Mapping,Information Engineering University,Zhengzhou 450052,China;2.Naval Research Institute of Surveying and Mapping,Tianjin 300061,China
The rigorous formulae of the boundary value problem(BVP)of physical geodesy complete to the first-order solution were derived considering corrections of terrain and complete spherical Bouguer anomaly’s gradient term.The formulae were tested through computation of a regional quasi-geoid based on terrain,gravity and GPS/leveling observations within a selected experiment area.Comparison and analysis of the calculated height anomalies in term of absolute and relative accuracy show that the effect of complete spherical Bouguer anomaly’s gradient term on height anomaly could reach centimeter level.Therefore,besides the effect of terrain corrections,the effect of complete spherical Bouguer anomaly’s gradient term on height anomaly should also be taken into account in order to improve the accuracy of regional quasi-geoid modeling especially in areas with rolling terrain.
quasi-geoid;terrain correction;complete spherical Bouguer anomaly’s gradient;singular integral;height anomaly
LI Shanshan(1970—),female,PhD,professor,majors in physical geodesy.
LI Shanshan,WU Xiaoping,ZHANG Chuanding,et al.Regional Quasi-Geoid Refining Considering Corrections of Terrain and Complete Spherical Bouguer Anomaly’s Gradient Term[J].Acta Geodaetica et Cartographica Sinica,2012,41(4):510-516.(李姍姍,吳曉平,張傳定,等.顧及地形與完全球面布格異常梯度項改正的區(qū)域似大地水準(zhǔn)面精化[J].測繪學(xué)報,2012,41(4):510-516.)
P223
A
1001-1595(2012)04-0510-07
國家自然科學(xué)基金(41174026);國家重大科學(xué)儀器設(shè)備開發(fā)專項(2011YQ12004503)
宋啟凡)
2011-11-22
2012-02-08
李姍姍(1970—),女,博士,教授,研究方向為物理大地測量。