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

基于CGCS2000橢球的大地測(cè)量實(shí)用公式

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

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

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

CGCS2000橢球；實(shí)用公式；曲率半徑；子午線弧長(zhǎng)；底點(diǎn)緯度；白塞爾大地主題；高斯投影

0 引言

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

1 CGCS2000橢球參數(shù)

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

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

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

(1)

2.1 正弦級(jí)數(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 余弦級(jí)數(shù)展開式

(4)

(5)

2.3 算例

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

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

表2 主曲率半徑計(jì)算

3 子午線弧長(zhǎng)與底點(diǎn)緯度

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

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

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

16.832 6sin4B-0.022sin6B

(6)

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

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

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

133.959 8sin3B+0.697 5sin5B)cosB

(7)

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

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

135.366 9cos3B+0.709 5cos5B)sinB

(8)

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

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

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

10-10sin8β

(9)

3.5 用余弦升冪多項(xiàng)式表示的底點(diǎn)緯度公式

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)可計(jì)算子午線弧長(zhǎng)，將計(jì)算出的子午線弧長(zhǎng)代入式(9)和式(10)，計(jì)算出的底點(diǎn)緯度應(yīng)該與已知的緯度相同。計(jì)算結(jié)果顯示在表3中。

表3 子午線弧長(zhǎng)與底點(diǎn)緯度計(jì)算

4 白塞爾大地主題解算

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

4.1 正算時(shí)的A、B、C

(11)

4.2 正算時(shí)的α、β

(12)

4.3 反算時(shí)的A、B″、C″

(13)

4.4 反算時(shí)的α、β′

(14)

4.5 算例

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

表4 白塞爾大地主題解算

5 高斯投影正反算公式

5.1 正算公式

(15)

式(15)中

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

5.2 反算公式

(16)

式(16)中

(17)

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

5.3 算例

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

表5 高斯投影正反算

6 結(jié)束語

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

[11]中國(guó)衛(wèi)星導(dǎo)航系統(tǒng)管理辦公室.北斗衛(wèi)星導(dǎo)航系統(tǒng)空間信號(hào)接口控制文件公開服務(wù)信號(hào)(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橢球的大地測(cè)量實(shí)用公式[J].導(dǎo)航定位學(xué)報(bào),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—)，男，山西孝義人，教授，主要從事大地測(cè)量的研究和教學(xué)。

P228

A

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