利用垂線偏差計算大地水準面中央區效應的改進方法

李厚樸，邊少鋒

1.海軍工程大學導航工程系，湖北武漢430033；2.海島（礁）測繪技術國家測繪地理信息局重點實驗室，山東 青島266510；3.中國科學院測量與地球物理研究所，湖北武漢430077

利用垂線偏差計算大地水準面中央區效應的改進方法

李厚樸1，2，邊少鋒1，3

1.海軍工程大學導航工程系，湖北武漢430033；2.海島（礁）測繪技術國家測繪地理信息局重點實驗室，山東 青島266510；3.中國科學院測量與地球物理研究所，湖北武漢430077

為提高利用垂線偏差計算大地水準面中央區效應的精度，視中央區為矩形域，將垂線偏差分量表示成雙二次多項式插值形式，引入非奇異變換，推導出中央區效應的計算公式。垂線偏差理論模型下的分析表明傳統公式的誤差與緯度以及垂線偏差子午分量沿經線方向的變化與卯酉分量沿緯線方向的變化之間的比值有關；以中緯度某區域分辨率為2′×2′的垂線偏差數據為背景場進行實際計算，結果表明當中央區為計算點所在的1個網格時，傳統公式與該公式計算得到的中央區效應差值的最大值達數厘米。該公式可為大地水準面中央區效應的高精度計算提供理論依據。

衛星測高；Molodensky公式；非奇異變換；大地水準面；中央區效應

1 引 言

大地水準面表征地球的基本幾何和物理特性，是對定義和建立地理坐標系統起基準作用的一個曲面［1］。高精度高分辨率局部或區域大地水準面為測繪學、地球物理、地球動力學及海洋學等地球科學的研究和應用提供基礎地球空間信息［2－3］。

海域大地水準面的研究和計算是建立大地水準面模型必須解決的問題。衛星測高技術的出現使人們可以在全球范圍內以較高的精度和分辨率觀測海面高。根據豐富的海面高觀測值可以計算出密集的垂線偏差，并且這一過程可以消除地理位置相關的徑向軌道誤差以及長波海面地形等類似系統誤差。因此，利用垂線偏差計算大地水準面的Molodensky公式備受關注，尤其是與快速Fourier變換算法相結合，在大地水準面計算中得到廣泛應用［4－10］。

實際計算中，計算點及其附近區域（中央區）到計算點的理論距離接近于零，導致Molodensky公式中的積分奇異。由于衛星測高數據分辨率的原因，中央區實際上是一個數平方千米乃至數十平方千米的區域，它對大地水準面的貢獻不能不加考慮地忽略。為此，文獻［8］將中央區視為圓域，推導出中央區效應的計算公式。然而，實際計算所用數據通常為網格化分布，由于子午線收斂的影響，中央區更近似于矩形域，因此，將中央區視為圓域的處理方法與數據真實分布情況并不相符，由此產生的誤差在高精度大地水準面計算中能否被忽略值得深入研究。有鑒于此，本文將中央區視為矩形域，引入非奇異變換［11－14］，推導出大地水準面中央區效應的精密計算公式，并對導出公式和傳統公式在實際計算中存在的差異進行分析比較。

2 Molodensky公式及其平面近似形式

根據文獻［6］，直接寫出Molodensky公式

式中，N為大地水準面高；R＝6 371km為地球平均半徑；αQP代表流動點Q至計算點P的方位角；ψPQ為P、Q之間的球面距離；ξQ和ηQ分別代表流動點Q處的垂線偏差在子午圈方向和卯酉圈方向上的分量。

圖1 局部坐標系示意圖Fig.1 The sketch map of the local coordinate system

為對中央區進行計算，首先定義以計算點P為原點的局部切平面直角坐標系，如圖1所示。x軸指向北極，y軸指向東且xy平面過P點與地球表面相切，Q（x，y）為中央區內流動的積分點，l為計算點與流動點的平面距離，且有l＝由圖1可知

3 中央區效應的精密計算

為得到中央區為矩形域時的大地水準面計算公式，需要解決式（3）中的奇異積分計算問題。奇異積分的計算一直是地球物理學中的一個難點問題，制約著重力場泛函的計算精度。文獻［11］提出一組“非奇異變換”，系統地解決了物理大地測量中高程異常、垂線偏差、地形改正以及重力異常梯度等泛函的中央區計算問題［11－14］。這組“非奇異變換”不僅可使中央區直接數值積分，而且可簡化有關的解析推導，值得推廣應用。為此，本文首先將中央區垂線偏差分量表示為雙二次多項式插值形式，之后利用“非奇異變換”對式（3）中的奇異積分進行了處理，推導出大地水準面中央區效應的精密計算公式。

3.1 中央區垂線偏差分量雙二次多項式插值表示

如圖2所示，設中央區為σ：［－a＜x＜a，－b＜y＜b］，其中a＝1，b＝cosφ。中央區共包含4個網格單元，9個網格節點。為充分反映中央區內垂線偏差的變化細節，同時為提高計算精度，將垂線偏差子午分量ξQ、卯酉分量ηQ表示成雙二次多項式插值形式

圖2 垂線偏差雙二次多項式插值表示時的中央區示意圖Fig.2 The sketch map of the innermost area when components of deflections of the vertical are expressed as double quadratic polynomials

為確定出待定系數αij和βij，將式（4）和式（5）改寫為

則有

將圖2中網格節點處的垂線偏差子午分量ξij＝ξ（i，jcosφ）和卯酉分量ηij＝η（i，jcosφ）（i，j＝－1，0，1）作為插值條件代入式（6）和式（7），可得

以上兩式可以簡化為

式中

因此

將式（15）和式（16）分別代入式（8）和式（9）即可確定出待定系數αij和βij。

3.2 中央區效應的精密計算公式

為推導方便，記式（3）中的奇異積分為

式中

將式（4）和式（5）分別代入式（18）和式（19），考慮到奇偶函數的積分性質，可得

由計算點P向右上頂點作連線分右上象限為σ1：［0＜x＜1，0＜y＜bx］和σ2：［0＜x＜b－1y，0＜y＜b］，如圖2所示。對σ1引入非奇異變換［11－14］

對σ2引入非奇異變換［11－14］

則三角形σ1映射為矩形σ′1：［0＜x＜1，0＜k＜b］，三角形σ2映射為矩形σ′2：［0＜λ＜b－1，0＜y＜b］。

將式（22）和式（23）分別代入式（20）和式（21），并注意到兩式中的被積函數均為偶函數，可得

具體的積分值可在計算機代數系統Mathematica［15］下求得

因此，本文方法確定出的式（3）中的奇異積分為

考慮到一般情況下a非單位長度，IM含長度量綱需乘以a，將式（28）代入式（3），可得中央區包含四個網格時，垂線偏差雙二次多項式插值表示下大地水準面中央區效應的計算公式為

4 算例分析

4.1 理論模型下的誤差分析

為分析本文和文獻［8］導出的大地水準面中央區效應計算公式的誤差，取中央區為σ：［－a＜x＜a，－b＜y＜b］，其中a＝1，b＝cosφ，垂線偏差分量取為如下二階泰勒展開式

將式（31）和式（32）代入式（3），考慮到奇偶函數的積分性質，可得大地水準面中央區效應為

引入非奇異變換式（22）和式（23）后，式（33）變換為

N即中央區包含4個網格時的大地水準面中央區效應真值。將由式（8）和式（9）確定的系數αij、βij代入式（29），可得本文方法確定的大地水準面中央區效應N1為

可見該方法確定的N1與非奇異變換后的真值N完全一致。

文獻［8］確定的大地水準面中央區效應NH為

相對誤差βH為

不同緯度φ和比值m處相對誤差βH的具體數值列于表1。為形象地描述βH隨φ、m的變化情況，分別繪制出了m＝－5、φ＝20°時βH的變化曲線，如圖3（a）、（b）所示。分析圖3，并結合表1，可以發現：

表1 不同緯度φ和比值m處的相對誤差βHTab.1 Numerical values of relative errorβHat different latitudesφand ratios m （%）

圖3 相對誤差βH的變化示意圖Fig.3 The sketch map of the characteristic ofβH

（1）m為定值時，βH是關于緯度φ的偶函數。m＝1或φ＝0°時，式（37）分子為零，βH＝0；m＝－1時除φ＝0°以外，βH＝－1。βH不為定值時其絕對值隨緯度φ的升高而增大。

（2）φ為定值時，式（37）分母在m0＝處為零，βH發生跳躍。當m＜m0時，βH為負值，隨m的增大而減??；當m＞m0時，βH亦隨m的增大而減小，并由正值變為負值。

4.2 實際數據下的中央區效應分析

選定52°N～58°N，2°W～8°W海域作為試算區，由最新的EGM2008地球重力場模型［16］計算得到了該區域2′×2′分辨率的垂線偏差數據，分別利用本文和文獻［8］導出公式計算該區域的中央區效應。記中央區為4′×4′即包含4個網格時本文和文獻［8］導出公式的計算結果分別為N1、NH；中央區為計算點所在的1個網格時本文和文獻［8］導出公式的計算結果分別為N2、N′H。兩種方法計算結果之間的比較情況如表2所示。

表2 兩種方法計算結果之間的比較Tab.2 Comparisons of the innermost effects calculated by the two methods respectively cm

由表2可以看出，當中央區包含4個網格時，兩種方法計算得到的中央區效應差值的均方差和標準差均為1.037cm，最大值高達16.943cm；當中央區為計算點所在的1個網格時，兩種方法計算得到的中央區效應差值的均方差和標準差均為0.263cm，最大值為4.288cm，這種差異需要在計算厘米級以及更高精度的大地水準面時加以考慮。

5 結 論

為提高利用垂線偏差計算大地水準面中央區效應的精度，將中央區視為與實際數據分布更為相符的矩形域，引入非奇異變換，推導出大地水準面中央區效應的精密計算公式，并基于理論模型和實際數據分析導出公式和傳統公式計算結果的差異。研究表明：

（1）將中央區垂線偏差表示成雙二次多項式插值形式，引入非奇異變換后，Molodensky公式平面近似式中的奇異積分可變換為直接可積的非奇異積分。非奇異變換在解決地球物理奇異問題中具有重要的應用價值，值得推廣使用。

（2）在給定的垂線偏差理論模型下，該公式的計算結果與真值完全一致，傳統公式誤差則與緯度以及垂線偏差子午分量沿經線方向的變化和卯酉分量沿緯線方向的變化之間的比值有關。

（3）實際數據下的中央區效應分析表明，當中央區為計算點所在的1個網格時，傳統公式與本文導出公式計算得到的中央區效應差值的最大值達數厘米，這種差異需要在計算厘米級以及更高精度的大地水準面時加以考慮。

［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.（陳俊勇，李建成，寧津生，等.我國大陸高精度、高分辨率大地水準面的研究和實施［J］.測繪學報，2001，30（2）：95－100.）

［2］ 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高分辨率厘米級高精度大地水準面的確定［J］.測繪學報，2003，32（2）：102－107.）

［3］ CHAO Dingbo.Refinement of Quasi－geoid in China and Relevant Problems［J］.Geomatics and Information Science of Wuhan University，2003，28（sup）：111－114.（晁定波.關于我國似大地水準面的精化及有關問題［J］.武漢大學學報：信息科學版，2003，28（增刊）：111－114.）

［4］ CHEN Junyong，LI Jiancheng，NING Jinsheng，et al.A New Chinese Geoid with High Resolution and High Accurac［J］.Geomatics and Information Science of Wuhan University，2001，26（4）：283－289.（陳俊勇，李建成，寧津生，等.中國新一代高精度、高分辨率大地水準面的研究和實施［J］.武漢大學學報：信息科學版，2001，26（4）：283－289.）

［5］ 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.（陳俊勇，李建成，寧津生，等.中國似大地水準面［J］.測繪學報，2002，31（增刊）：1－6.）

［6］ LI Jiancheng，CHEN Junyong，NING Jinsheng，et al.The Approximation Theory of the Earth’s Gravity Field and Determination of the 2000Quasi Geoid in China［M］.Wuhan：Publishing House of Wuhan University，2003.（李建成，陳俊勇，寧津生，等.地球重力場逼近理論與中國2000似大地水準面的確定［M］.武漢：武漢大學出版社，2003.）

［7］ LI Jiancheng，NING Jinsheng，CHEN Junyong，et al.Problem of Piecing Together Between China Mainlandgeoid and Altimetry－derived Geoid in China Sea Area［J］.Geomatics and Information Science of Wuhan University，2003，28（5）：542－546.（李建成，寧津生，陳俊勇，等.我國海域大地水準面與大陸大地水準面的拼接研究［J］.武漢大學學報：信息科學版，2003，28（5）：542－546.）

［8］ HWANG C.Inverse Vening Meinesz Formula and Deflection－geoid Formula：Applications to the Predictions of Gravity and Geoid over the South China Sea［J］.Journal of Geodesy，1998，72：304－312.

［9］ HWANG C，HSU H Y，JANG R J.Global Mean Sea Surface and Marine Gravity Anomaly from Multi－satellite Altimetry：Applications of Deflection－geoid and Inverse Vening Meinesz Formulae［J］.Journal of Geodesy，2002，76：407－418.

［10］ PENG Fuqing，ZHANG Ruihua，SHI Pan，et al.Marine Geoid from Satellite Altimeter Data［J］.Chinese Journal of Geophysics，2003，46（4）：462－466.（彭富清，張瑞華，石磐，等.基于衛星測高的海域大地水準面［J］.地球物理學報，2003，46（4）：462－466.）

［11］ BIAN Shaofeng.Numerical Solution for Geodetic Boundary Value Problem and the Earth’s Gravity Field Approximation［D］.Wuhan：Wuhan Technical University of Surveying and Mapping，1992：96－98.（邊少鋒.大地測量邊值問題數值解法與地球重力場逼近［D］.武漢：武漢測繪科技大學，1992：96－98.）

［12］ BIAN S F，SUN H Q.The Expression of Common Singular Integrals in Physical Geodesy［J］.Manuscripta Geodaetica，1994，19：62－69.

［13］ BIAN S F.Some Cubature Formulas for Singular Integrals in Physical Geodesy［J］.Journal of Geodesy，1997，71：443－453.

［14］ BIAN Shaofeng，XUE Fangxia.Discussion on the Calculation of Plumb Line Deviation［J］.Acta Geodaetica et Cartographica Sinica，1997，26（1）：33－36.（邊少鋒，薛芳俠.論地形垂線偏差中央區貢獻的計算［J］.測繪學報，1997，26（1）：33－36.）

［15］ BIAN Shaofeng，XU Jiangning.Computer Algebra System and Mathematical Analysis in Geodesy［M］.Beijing：National Defense Industy Press，2004：168－177.（邊少鋒，許江寧.計算機代數系統與大地測量數學分析［M］.北京：國防工業出版社，2004：168－177.）

［16］ ZHANG Chuanyin，GUO Chunxi，CHEN Junyong，et al.EGM2008and Its Application Analysis in Chinese Mailand［J］.Acta Geodaetica et Cartographica Sinica，2009，38（4）：283－289.（章傳銀，郭春喜，陳俊勇，等.EGM2008地球重力場模型在中國大陸適用性分析［J］.測繪學報，2009，38（4）：283－289.）

The Improved Method of Calculating the Geoid Innermost Area Effects Using Deflections of the Vertical

LI Houpu1，2，BIAN Shaofeng1，3
1.Department of Navigation，Naval University of Engineering，Wuhan430033，China；2.Key Laboratory of Surveying and Mapping Technology on Island and Reef，National Administration of Surveying，Mapping and Geoinformation，Qingdao 266510，China；3.Institute of Geodesy and Geophysics，Chinese Academy of Sciences，Wuhan 430077，China

In order to improve the precision of the innermost area effects in geoid calculation using deflections of the vertical，components of deflections of the vertical were expressed as double quadratic polynomials regarding the innermost area as a rectangular one，and the formulae to calculate the innermost area effects were derived after the non－singular transformation was introduced.The analysis based on the theoretical model of deflections of the vertical shows that errors of traditional formulae depend on the latitude and the ratio of the change of the northsouth deflection component along the meridian to that of the west－east deflection component along the parallel.A practical calculation was done based on deflections of the vertical data with a resolution of 2′×2′in a middle latitude area.The results indicate that the maximal differences of the innermost area effects calculated by traditional and this formulae are several centimeters when the innermost area is only agrid including the calculation point.The formulae could provide theoretical basis for the calculation of the geoid innermost area effects with high precision.Key words：satellite altimetry；Molodensky formula；non－singular transformation；geoid；innermost area effects

LI Houpu（1985—），male，PhD，lecturer，majors in geodesy and satellite navigation.

1001－1595（2011）06－0730－06

P223

A

國家自然科學基金（40774002；40904018；41071295）；海島（礁）測繪技術國家測繪地理信息局重點實驗室開放研究基金（2010B04）；海軍工程大學自然科學基金（HGDYDJJ009）

宋啟凡）

2010－03－03

2011－04－21

李厚樸（1985—），男，博士，講師，研究方向為大地測量和衛星導航。

E－mail：lihoupu1985＠126.com