基于CGCS2000橢球的大地測量實用公式

王愛生，徐 歡，張 棋，魏 猛

基于CGCS2000橢球的大地測量實用公式

王愛生，徐 歡，張 棋，魏 猛

(江蘇師范大學(xué) 測繪學(xué)院，江蘇 徐州 221116)

針對目前已出版的文獻(xiàn)中都沒有給出有關(guān)2000中國大地測量坐標(biāo)系統(tǒng)對應(yīng)橢球(CGCS2000橢球)的實用公式，根據(jù)大地測量學(xué)中有關(guān)橢球計算和高斯投影的基本公式，使用CGCS2000橢球參數(shù)，給出多種常用計算公式的實用公式，包括子午圈曲率半徑、卯酉圈曲率半徑、子午線弧長、底點緯度、白塞爾大地主題解算的A、B、C等系數(shù)、高斯投影等，有些公式形式簡單、使用方便，特別適合工程技術(shù)人員，有些公式有幾種表示方法，適合理論推導(dǎo)和演算。每個公式都附有算例來驗證其正確性，通過這些算例也能了解公式的適用性。

CGCS2000橢球；實用公式；曲率半徑；子午線弧長；底點緯度；白塞爾大地主題；高斯投影

0 引言

在大地測量計算時，會涉及到復(fù)雜的計算公式和長長的橢球參數(shù)。因此，希望有比較簡單的實用公式直接應(yīng)用。以往針對克拉索夫斯基橢球和IAG-75橢球的實用公式在許多文獻(xiàn)中都會找到。例如，針對克拉索夫斯基橢球的主曲率半徑的實用公式[1-3]、針對IAG-75橢球的主曲率半徑計算公式[1-2，4]、針對克拉索夫斯基橢球的子午線弧長實用公式[1-8]、針對IAG-75橢球的子午線弧長實用公式[1-2，4，6-8]、針對克拉索夫斯基橢球的底點緯度實用公式[1，3-5，7-8]、針對IAG-75橢球的底點緯度實用公式[1][4-5，7-8]；針對克拉索夫斯基橢球的白塞爾大地主題計算A、B、C等系數(shù)的實用公式[1-2]、針對IAG-75橢球的白塞爾大地主題計算A、B、C等系數(shù)的實用公式[1-2]、針對克拉索夫斯基橢球的高斯投影計算的實用公式[1-3，5]、針對IAG-75橢球的高斯投影計算實用公式[1-2，5]。從2008-07-01起，我國正式使用2000國家大地坐標(biāo)系(China Geodetic Coordinate System 2000，CGCS2000)，其對應(yīng)的橢球為CGCS2000橢球[9-10]。北斗衛(wèi)星導(dǎo)航系統(tǒng)(BeiDou navigation satellite system，BDS)使用的坐標(biāo)系也是CGCS2000[11]。但是目前針對CGCS2000橢球的實用公式還很少見到，文獻(xiàn)[6]曾給出針對CGCS2000橢球的子午線弧長的實用公式。因此，本文將不加推導(dǎo)給出主曲率半徑、子午線弧長、底點緯度、白塞爾大地主題解算的系數(shù)以及高斯投影正算和反算的實用公式，并通過算例檢核公式的正確性。

1 CGCS2000橢球參數(shù)

表1是CGCS2000橢球參數(shù)以及其它常用的橢球參數(shù)[1，6，10]。表1中a是橢球的長半徑，b是橢球的短半徑，f是橢球的扁率，e2是橢球的第一偏心率，e′2是橢球的第二偏心率，c是橢球的極曲率半徑。我國的1954年北京坐標(biāo)系使用克拉索夫斯基橢球，1980西安坐標(biāo)系使用IAG-75橢球。表中可見，CGCS2000橢球與WGS84橢球僅有微小的差別。

2 主曲率半徑級數(shù)展開公式

主曲率半徑是指子午圈的曲率半徑M和卯酉圈的曲率半徑N，基本公式是：

(1)

2.1 正弦級數(shù)展開式

M=m0+m2sin2B+m4sin4B+m6sin6B+

m8sin8B+m10sin10B

(2)

N=n0+n2sin2B+n4sin4B+n6sin6B+

n8sin8B+n10sin10B

(3)

式(2)中，m0=6 335 439.327 08，m2=63 617.757 70，m4=532.351 81，m6=4.157 73，m8=0.031 31，m10=0.000 23；式(3)中，n0=6 378 137.000 00，n2=21 348.836 46，n4=107.187 92，n6=0.597 96，n8=0.003 50，n10=0.000 02

2.2 余弦級數(shù)展開式

(4)

(5)

2.3 算例

選取不同的緯度B進(jìn)行代入公式(1)、(2)和(4)計算的子午線曲率半徑M，代入公式(1)、(3)和(5)計算的卯酉圈曲率半徑N，結(jié)果如表2。

表1 幾種常用的橢球參數(shù)

表2 主曲率半徑計算

3 子午線弧長與底點緯度

子午線弧長是在子午線上緯度從0到B的曲線距離，用X表示，以m為單位。底點緯度是子午線弧長對應(yīng)的緯度，用Bf表示，底點緯度的計算其實是已知子午線弧長計算對應(yīng)的緯度。

3.1 表示成正弦的倍數(shù)函數(shù)的子午線弧長公式

X=111 132.952 54B°-16 038.508 69sin2B+

16.832 6sin4B-0.022sin6B

(6)

文獻(xiàn)[6]給出的子午線弧長的表達(dá)方式與上式相同，但是取到sin8B項。同時，各個系數(shù)的有效數(shù)字取位與上式不同。

3.2 表示成正弦的n次冪和余弦的乘積的子午線弧長公式

X=111 132.952 54B°-(32 009.818 6sinB+

133.959 8sin3B+0.697 5sin5B)cosB

(7)

3.3 表示成余弦的n次冪與正弦的乘積的子午線弧長公式

X=111 132.952 54B°-(32 144.480 0cosB-

135.366 9cos3B+0.709 5cos5B)sinB

(8)

3.4 用正弦倍數(shù)函數(shù)表示的底點緯度公式

Bf=β+2.518 826 589×10-3sin2β+3.701 005×

10-6sin4β+7.447×10-9sin6β+1.1×

10-10sin8β

(9)

3.5 用余弦升冪多項式表示的底點緯度公式

Bf=β+(50 228 929.6+(293 785.7+(2 171.0+

141.4cos2β)cos2β)cos2β)×10-10cosβsinβ

(10)

3.6 算例

已知緯度，根據(jù)式(6)、式(7)和式(8)可計算子午線弧長，將計算出的子午線弧長代入式(9)和式(10)，計算出的底點緯度應(yīng)該與已知的緯度相同。計算結(jié)果顯示在表3中。

表3 子午線弧長與底點緯度計算

4 白塞爾大地主題解算

大地主題解算包括大地主題正解和反解。設(shè)在橢球面上有長度為S的大地線，其兩端點為P1和P2，如果已知P1的大地坐標(biāo)(B1，L1)和在P1點處的大地方位角A1，計算P2的大地坐標(biāo)(B2，L2)和在P2點處的大地方位角A2，稱為大地主題正解；如果已知(B1，L1)和(B2，L2)，計算S、A1和A2，稱為大地主題反解。有關(guān)白塞爾大地主題解算的所有公式都在文獻(xiàn)[1]中，其中涉及到的幾個系數(shù)與橢球參數(shù)有關(guān)。下面給出的系數(shù)與文獻(xiàn)[1]中的系數(shù)是對應(yīng)的，式中的有關(guān)符號的含義也與文獻(xiàn)[1]相同。

4.1 正算時的A、B、C

(11)

4.2 正算時的α、β

(12)

4.3 反算時的A、B″、C″

(13)

4.4 反算時的α、β′

(14)

4.5 算例

已知B1=35°，L1=114°，A1=25°，S分別等于1 000 m、10 000 m等，代入貝塞爾正解公式計算(B2，L2)和A2，將正解的結(jié)果(B2，L2)連同已知的(B1，L1)代入反算公式，計算出來的S(表中為S′)和A1(表中為A1′)應(yīng)與已知S和A1相同，并且正解的A2與反解的A2(表中為A2′)相同。

表4 白塞爾大地主題解算

5 高斯投影正反算公式

5.1 正算公式

(15)

式(15)中

式(15)中，B和L的單位是弧度，L0是中央子午線經(jīng)度。N的計算式其實就是式(5)，只不過取的項數(shù)和數(shù)字的位數(shù)不同而已。當(dāng)然也可以使用式(3)來計算。

5.2 反算公式

(16)

式(16)中

(17)

式(17)中，B、L的單位均為弧度，Bf的計算式與式(9)相同。當(dāng)然也可以使用式(10)。

5.3 算例

由B和L計算x和y稱為高斯投影正算，由x和y計算B和L稱為高斯投影反算。將正算的結(jié)果x和y代入反算公式，計算出來的B(表5中的B′)和L(表5中的L′)應(yīng)該與初始的B和L相同。

表5 高斯投影正反算

6 結(jié)束語

本文給出的幾個公式，通過數(shù)值試驗，證明是正確的。

1)子午線曲率半徑正弦級數(shù)展開式和余弦級數(shù)展開式都取至sin或cos的10次方項，可精確到0.05 mm；卯酉圈曲率半徑的正弦級數(shù)展開式取至sin的10次方項，可精確到0.05 mm，余弦級數(shù)展開到cos的10次方項，可精確到0.1 mm。

2)子午線弧長的三個公式中，式(6)和式(7)在緯度在30°左右時，僅相差0.1 mm，在45°左右時相差0.3 mm，在60°以上時可能相差0.7 mm，而式(8)與式(6)和式(7)相差較大，最大達(dá)0.9 mm。因此推薦使用式(6)和式(7)。

3)用底點緯度公式(9)計算出來的緯度與已知的緯度最大相差0.000 02″，對應(yīng)的球面距離約為0.6 mm，而用式(10)計算出來的緯度與已知緯度最大相差0.000 03″，對應(yīng)的球面距離約為0.9 mm，因此，在計算底點緯度時推薦使用式(9)。

4)利用本文給出的系數(shù)進(jìn)行白塞爾大地主題解算時，當(dāng)S在1 000 km以下時，正反算檢核的效果很好，反算的大地線長度與已知的大地線長度最大相差0.014 m，而方位角最大相差0.005 1″。即使S達(dá)到10 000 km，長度相差僅為0.222 m，方位角相差僅為0.117 1″。

5)將高斯投影正算結(jié)果代入反算公式，當(dāng)距中央子午線小于5 km時，與原始值只差0.000 01″，距中央子午線300 km時，只差0.000 04″。由此說明高斯投影正反算實用公式能夠精確到1 mm。

[1] 孔祥元，郭際明，劉宗泉.大地測量學(xué)基礎(chǔ)[M].武漢：武漢大學(xué)出版社，2010.

[2] 張華海，王寶山，趙長勝，等.應(yīng)用大地測量學(xué)[M].徐州：中國礦業(yè)大學(xué)出版社，2011.

[3] 張鳳舉，邢永昌.礦區(qū)控制測量(上冊)[M].北京：煤炭工業(yè)出版社，1987.

[4] 黨亞民，成英燕，薛樹強(qiáng).大地坐標(biāo)系統(tǒng)及其應(yīng)用[M].北京：測繪出版社，2010.

[5] 武漢測繪學(xué)院控制測量教研室，同濟(jì)大學(xué)大地測量教研室.控制測量學(xué)(下冊)[M].北京：測繪出版社，1988.

[6] 呂志平，喬書波.大地測量學(xué)基礎(chǔ)[M].北京：測繪出版社，2010.

[7] 郭群長，李仲勤，李輝.大地測量學(xué)概論[M].成都：西南交通大學(xué)出版社，2009.

[8] 李玉寶，沈?qū)W標(biāo)，吳向陽.控制測量學(xué)[M].南京：東南大學(xué)出版社，2013.

[9] 陳俊勇.中國現(xiàn)代大地基準(zhǔn)——中國大地坐標(biāo)系統(tǒng)2000(CGCS2000)及其框架[J].測繪學(xué)報，2008，37(8)：269-271.

[10]程鵬飛，文漢江，成英燕，等.2000國家大地坐標(biāo)系橢球參數(shù)與GRS80和WGS84的比較[J].測繪學(xué)報，2009，38(3)：189-194.

[11]中國衛(wèi)星導(dǎo)航系統(tǒng)管理辦公室.北斗衛(wèi)星導(dǎo)航系統(tǒng)空間信號接口控制文件公開服務(wù)信號(2.0 版)[EB/OL].(2013-12-01)[2015-04-21].http：//www.beidou.gov.cn/attach/2013/12/26/2013122604a521b35b7f4a54b44cfbbc8abd74a8.pdf.

Practical Formulas for Geodetic Surveying Based on CGCS2000 Ellipsoid

WANGAi-sheng，XUHuan，ZHANGQi，WEIMeng

(School of Geomatics and Geodesy，Jiangsu Normal University，Xuzhou 221116，China)

Because published literature did not give practical formula corresponding to Chinese Geodetic Coordinate System 2000 ellipsoid(CGCS2000 ellipsoid)，according to the basic principle of the ellipsoid calculation and Gauss projection in geodesy，using CGCS2000 ellipsoid parameters，many practical calculation formula are presented.These formulas consist of curvature radius of the meridian，curvature radius of the prime vertical，meridian arc length，latitude of pedal，Bessel formula coefficients A，B，C etc for solution of geodetic problem and Guass project.Some formulas are simple and easy to use，especially for engineering technicians.Some formulas have several representations，suitable for theoretical derivation and calculation.Each formula is accompanied by an example to verify its correctness，and can also be used to understand the applicability of the formula.

CGCS2000 ellipsoid；practical formula；radius of curvature；meridian arc length；latitude of pedal；Bessel solution of geodetic problem；Guass project

王愛生，徐歡，張棋，等.基于CGCS2000橢球的大地測量實用公式[J].導(dǎo)航定位學(xué)報,2015,3(3):105-109+131 (.WANG Ai-sheng,XU Huan,ZHANG Qi,et al.Practical Formulas for Geodetic Surveying Based on CGCS2000 Ellipsoid[J].Journal of Navigation and Positioning,2015,3(3):105-109+131.)

10.16547/j.cnki.10-1096.20150321.

2015-05-18

江蘇師范大學(xué)科研基金(09XLR18)。

王愛生(1965—)，男，山西孝義人，教授，主要從事大地測量的研究和教學(xué)。

P228

A

2095-4999(2015)-03-0105-05