李厚樸，邊少鋒

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

利用垂線偏差計算大地水準(zhǔn)面中央?yún)^(qū)效應(yīng)的改進方法

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

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

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

衛(wèi)星測高；Molodensky公式；非奇異變換；大地水準(zhǔn)面；中央?yún)^(qū)效應(yīng)

1 引 言

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

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

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

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

根據(jù)文獻［6］，直接寫出Molodensky公式

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

圖1 局部坐標(biāo)系示意圖Fig.1 The sketch map of the local coordinate system

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

3 中央?yún)^(qū)效應(yīng)的精密計算

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

3.1 中央?yún)^(qū)垂線偏差分量雙二次多項式插值表示

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

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

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

則有

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

以上兩式可以簡化為

式中

因此

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

3.2 中央?yún)^(qū)效應(yīng)的精密計算公式

為推導(dǎo)方便，記式（3）中的奇異積分為

式中

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

由計算點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），并注意到兩式中的被積函數(shù)均為偶函數(shù)，可得

具體的積分值可在計算機代數(shù)系統(tǒng)Mathematica［15］下求得

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

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

4 算例分析

4.1 理論模型下的誤差分析

為分析本文和文獻［8］導(dǎo)出的大地水準(zhǔn)面中央?yún)^(qū)效應(yīng)計算公式的誤差，取中央?yún)^(qū)為σ：［－a＜x＜a，－b＜y＜b］，其中a＝1，b＝cosφ，垂線偏差分量取為如下二階泰勒展開式

將式（31）和式（32）代入式（3），考慮到奇偶函數(shù)的積分性質(zhì)，可得大地水準(zhǔn)面中央?yún)^(qū)效應(yīng)為

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

N即中央?yún)^(qū)包含4個網(wǎng)格時的大地水準(zhǔn)面中央?yún)^(qū)效應(yīng)真值。將由式（8）和式（9）確定的系數(shù)αij、βij代入式（29），可得本文方法確定的大地水準(zhǔn)面中央?yún)^(qū)效應(yīng)N1為

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

文獻［8］確定的大地水準(zhǔn)面中央?yún)^(qū)效應(yīng)NH為

相對誤差βH為

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

表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是關(guān)于緯度φ的偶函數(shù)。m＝1或φ＝0°時，式（37）分子為零，βH＝0；m＝－1時除φ＝0°以外，βH＝－1。βH不為定值時其絕對值隨緯度φ的升高而增大。

（2）φ為定值時，式（37）分母在m0＝處為零，βH發(fā)生跳躍。當(dāng)m＜m0時，βH為負(fù)值，隨m的增大而減?。划?dāng)m＞m0時，βH亦隨m的增大而減小，并由正值變?yōu)樨?fù)值。

4.2 實際數(shù)據(jù)下的中央?yún)^(qū)效應(yīng)分析

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

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

由表2可以看出，當(dāng)中央?yún)^(qū)包含4個網(wǎng)格時，兩種方法計算得到的中央?yún)^(qū)效應(yīng)差值的均方差和標(biāo)準(zhǔn)差均為1.037cm，最大值高達16.943cm；當(dāng)中央?yún)^(qū)為計算點所在的1個網(wǎng)格時，兩種方法計算得到的中央?yún)^(qū)效應(yīng)差值的均方差和標(biāo)準(zhǔn)差均為0.263cm，最大值為4.288cm，這種差異需要在計算厘米級以及更高精度的大地水準(zhǔn)面時加以考慮。

5 結(jié) 論

為提高利用垂線偏差計算大地水準(zhǔn)面中央?yún)^(qū)效應(yīng)的精度，將中央?yún)^(qū)視為與實際數(shù)據(jù)分布更為相符的矩形域，引入非奇異變換，推導(dǎo)出大地水準(zhǔn)面中央?yún)^(qū)效應(yīng)的精密計算公式，并基于理論模型和實際數(shù)據(jù)分析導(dǎo)出公式和傳統(tǒng)公式計算結(jié)果的差異。研究表明：

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

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

（3）實際數(shù)據(jù)下的中央?yún)^(qū)效應(yīng)分析表明，當(dāng)中央?yún)^(qū)為計算點所在的1個網(wǎng)格時，傳統(tǒng)公式與本文導(dǎo)出公式計算得到的中央?yún)^(qū)效應(yīng)差值的最大值達數(shù)厘米，這種差異需要在計算厘米級以及更高精度的大地水準(zhǔn)面時加以考慮。

［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］ 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.）

［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.（晁定波.關(guān)于我國似大地水準(zhǔn)面的精化及有關(guān)問題［J］.武漢大學(xué)學(xué)報：信息科學(xué)版，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.（陳俊勇，李建成，寧津生，等.中國新一代高精度、高分辨率大地水準(zhǔn)面的研究和實施［J］.武漢大學(xué)學(xué)報：信息科學(xué)版，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.（陳俊勇，李建成，寧津生，等.中國似大地水準(zhǔn)面［J］.測繪學(xué)報，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似大地水準(zhǔn)面的確定［M］.武漢：武漢大學(xué)出版社，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.（李建成，寧津生，陳俊勇，等.我國海域大地水準(zhǔn)面與大陸大地水準(zhǔn)面的拼接研究［J］.武漢大學(xué)學(xué)報：信息科學(xué)版，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.（彭富清，張瑞華，石磐，等.基于衛(wèi)星測高的海域大地水準(zhǔn)面［J］.地球物理學(xué)報，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.（邊少鋒.大地測量邊值問題數(shù)值解法與地球重力場逼近［D］.武漢：武漢測繪科技大學(xué)，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.（邊少鋒，薛芳俠.論地形垂線偏差中央?yún)^(qū)貢獻的計算［J］.測繪學(xué)報，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.（邊少鋒，許江寧.計算機代數(shù)系統(tǒng)與大地測量數(shù)學(xué)分析［M］.北京：國防工業(yè)出版社，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］.測繪學(xué)報，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

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

宋啟凡）

2010－03－03

2011－04－21

李厚樸（1985—），男，博士，講師，研究方向為大地測量和衛(wèi)星導(dǎo)航。

E－mail：lihoupu1985＠126.com