利用垂线偏差计算大地水准面中央区效应的改进方法

李厚朴，边少锋

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