李厚樸，邊少鋒

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

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

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

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

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

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

1 引 言

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

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

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

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

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

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

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

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

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

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

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

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

圖2 垂線偏差雙二次多項(xiàng)式插值表示時(shí)的中央?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é)點(diǎn)處的垂線偏差子午分量ξij＝ξ（i，jcosφ）和卯酉分量ηij＝η（i，jcosφ）（i，j＝－1，0，1）作為插值條件代入式（6）和式（7），可得

以上兩式可以簡(jiǎn)化為

式中

因此

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

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

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

式中

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

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

對(duì)σ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ù)，可得

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

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

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

4 算例分析

4.1 理論模型下的誤差分析

為分析本文和文獻(xiàn)［8］導(dǎo)出的大地水準(zhǔn)面中央?yún)^(qū)效應(yīng)計(jì)算公式的誤差，取中央?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個(gè)網(wǎng)格時(shí)的大地水準(zhǔn)面中央?yún)^(qū)效應(yīng)真值。將由式（8）和式（9）確定的系數(shù)αij、βij代入式（29），可得本文方法確定的大地水準(zhǔn)面中央?yún)^(qū)效應(yīng)N1為

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

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

相對(duì)誤差βH為

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

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

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

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

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

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

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

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

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

5 結(jié) 論

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

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

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

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

［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.（陳俊勇，李建成，寧津生，等.我國(guó)大陸高精度、高分辨率大地水準(zhǔn)面的研究和實(shí)施［J］.測(cè)繪學(xué)報(bào)，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í)高精度大地水準(zhǔn)面的確定［J］.測(cè)繪學(xué)報(bào)，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)于我國(guó)似大地水準(zhǔn)面的精化及有關(guān)問題［J］.武漢大學(xué)學(xué)報(bào)：信息科學(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.（陳俊勇，李建成，寧津生，等.中國(guó)新一代高精度、高分辨率大地水準(zhǔn)面的研究和實(shí)施［J］.武漢大學(xué)學(xué)報(bào)：信息科學(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.（陳俊勇，李建成，寧津生，等.中國(guó)似大地水準(zhǔn)面［J］.測(cè)繪學(xué)報(bào)，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.（李建成，陳俊勇，寧津生，等.地球重力場(chǎng)逼近理論與中國(guó)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.（李建成，寧津生，陳俊勇，等.我國(guó)海域大地水準(zhǔn)面與大陸大地水準(zhǔn)面的拼接研究［J］.武漢大學(xué)學(xué)報(bào)：信息科學(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)星測(cè)高的海域大地水準(zhǔn)面［J］.地球物理學(xué)報(bào)，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.（邊少鋒.大地測(cè)量邊值問題數(shù)值解法與地球重力場(chǎng)逼近［D］.武漢：武漢測(cè)繪科技大學(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ū)貢獻(xiàn)的計(jì)算［J］.測(cè)繪學(xué)報(bào)，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.（邊少鋒，許江寧.計(jì)算機(jī)代數(shù)系統(tǒng)與大地測(cè)量數(shù)學(xué)分析［M］.北京：國(guó)防工業(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地球重力場(chǎng)模型在中國(guó)大陸適用性分析［J］.測(cè)繪學(xué)報(bào)，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

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

宋啟凡）

2010－03－03

2011－04－21

李厚樸（1985—），男，博士，講師，研究方向?yàn)榇蟮販y(cè)量和衛(wèi)星導(dǎo)航。

E－mail：lihoupu1985＠126.com