高斯投影正反算的直接计算结果是什么

一、基本思想：

高斯投影正算公式就是由大地坐标（L ，B ）求解高斯平面坐标（x ，y ），而高斯投影反算公式则是由高斯平面坐标（x ，y)求解大地坐标（L ，B).

二、计算模型:

基本椭球参数：

椭球长半轴a

椭球扁率f

椭球短半轴:(1)b a f =-

椭球第一偏心率

第 1 页

高斯计_高精闷慎度_高效率_高品质

全面覆盖电离辐射， 非电离辐射及直流磁场的检测。低频电磁场还可在线监测， 多种功能供用户选配

点击立即咨询，了解更多详情

咨询

深圳市柯雷科技开发 广告

：e = 椭球第二偏心率

：e '=高斯投影正算公式：此公式换算的精度为0.001m

64256

442234

22)5861(cos sin 720)495(cos 24cos sin 2l t t B B N l t B simB N l B B N X x ''+-''+''++-''+''⋅'蚂祥敬'+

=ρηηρρ 52224255

32233

)5814185(cos 120)1(cos 6cos l t t t B N l t B

第 2 页

N l B N y ''-++-''+''+-''+''⋅''=ηηρηρρ

其中：角度都为弧度

B 为点的纬度,0l L L ''=-，L 为点的经度,0L 为中央子午线经度；

N 为子午圈曲率半径,1222

(1sin )N a e B -=-；

tan t B =； 222cos e B η'=

180

3600ρπ''=*

其中X 为子午线弧长：

2402464661616sin cos ()(2)sin sin 33X a B B B a a a a a B a B ⎡⎤=--++-+⎢⎥⎣⎦

第 3 页

02468,,,,a a a a a 为基本常量，按如下公式计算：

2004682426844686868

83535281612815722321637816323216128m a m m m m m m a m m m a m m m m a m a ⎧=++++⎪⎪⎪=+++⎪⎪⎪=++⎨⎪⎪=+⎪⎪⎪=⎪⎩

02468,,,,m m m m m 为基本常量，按如下公式计算：

22222020426486379(1)5268

m a e m e m m e m m e m m e m =-====

高斯投影反算公式：此公式换算的宴中精度为0.0001’'。

()()()()22222432465

第 4 页

3

2235

2422250

53922461904572012cos 6cos 5282468120cos f f f f f f f f f f f f f f f f f f f f f f

f f f f f f f

t t B B y t t y

M N M N t y t t y

M N y y l t N B N B y t t t N B L l L ηηηηη=-

+++--++=-+++++++=+

其中: 0L 为中央子午线经度。

第 5 页

f B 为底点纬度，也就是当x X =时的子午线弧长所对应的纬度。按照子午线弧长公式：68240sin 2sin 4sin 6sin82468

a a a a X a B B B B B =-+-+，迭代进行计算； 初始开始时设：10f B X a =

以后每次迭代按下式计算:

10

6824(())()sin 2sin 4sin 6sin82468i

i f f i

i i i i f

f f f f B X X F B a a a a a F B B B B B +=-=-+-+

重复迭代至1i

第 6 页

i

f f B B ε+-<；为止。

1222

(1sin )f f N a e B -=-

3

2222(1)(1sin )f f M a e e B -=-- tan f f t B =；

222cos f

f e B η'=

海福特椭球（1910)我国52年以前基准椭球 a=6378388m b=6356911。9461279m α=0.33670033670

克拉索夫斯基椭球（1940 Krassovsky)

第 7 页

北京54坐标系基准椭球 a=6378245m b=6356863。018773m α=0.33523298692

1975年I 。U 。G.G 推荐椭球（国际大地测量协会1975) 西安80坐标系基准椭球

a=6378140m b=6356755.2881575m α=0.0033528131778

WGS-84椭球（GPS 全球定位系统椭球、17届国际大地测量协会） WGS —84 GPS 基准椭球

a=6378137m b=6356752.3142451m α=0。00335281006247

三、程序代码函数:

/*＊＊**＊＊＊**＊＊高斯投影正算函数***＊＊**＊＊**＊*＊*

第 8 页

输入 : double a ，f 椭球参数，B,L 为大地坐标,L0为中央子午线的经度，单位为弧度，x,y 为高斯平面坐标,y 加上了500000常量

返回:none

＊＊＊**＊*＊*＊*＊＊＊＊*＊＊＊＊*＊*＊*＊*＊**＊＊*＊＊*＊＊*＊＊＊/

void gaosiforward (double a ,double f ，double B ，double L ，double L0，double ＆x ，double ＆y ） ｛

double b , c ，e1， e2； //短半轴，极点处的子午线曲率半径，第一偏心率，第二偏心率

double l ， W ,N ， M ， daihao ；//W 为常用辅助函数，N 为子午圈曲率半径，M 为卯酉圈曲率半径

第 9 页

double X //子午线弧长，高斯投影的坐标

double ruo ， ita , sb , cb ，t

double m [5］，n ［5]；

//计算一些基本常量

{

b =a ＊(1-f )；

e1=sqrt (a ＊a -b ＊b )/a

e2=sqrt （a *a -b ＊b ）/b

c =a ＊a /b ；

m ［0］=a *(1-e1＊e1） m ［1］=3＊（e1＊e1*m [0])/2。0

m[2］=5*(e1＊e1*m［1]）/4。0

第 10 页

m[3]=7*（e1*e1*m［2]）/6.0

m[4］=9*(e1*e1*m[3]）/8。0；

n[0]=m[0］+m［1]/2+3*m［2］/8+5＊m［3]/16+35*m[4]/128

n[1］=m［1］/2+m[2］/2+15*m[3]/32+7＊m[4]/16；

n[2］=m［2]/8+3＊m[3]/16+7*m［4]/32

n［3]=m［3]/32+m[4]/16；

n［4］=m［4]/128 /////by kjh 2014。5。22 把改成了

｝

//由纬度计算子午线弧长

第 11 页

{

X=n［0]*B—sin（B)*cos(B）＊（（n［1]-n[2］+n[3]）+（2*n［2]-（16*n[3］/3.0))＊sin(B)＊sin(B）+16*n[3]*pow(sin(B)，4)/3。0)

}

l=L—L0；//弧度

ita=e2＊cos（B)

sb=sin(B）；

cb=cos(B）；

W=sqrt（1-e1*e1*sb*sb）；

N=a/W；

第 12 页

t=tan(B)

ruo=（180/Pi)＊3600

x=(X+N＊sb*cb*l＊l/2+N＊sb＊cb*cb＊cb*（5—t*t+9＊ita*ita+4＊ita*ita＊ita＊ita）＊l*l*l ＊l/24+N＊sb*cb*cb*cb*cb*cb＊（61-58*t*t+t＊t*t*t)＊l*l＊l*l*l*l/720）

y=（N＊cb*l+N*cb*cb*cb＊（1-t*t+ita＊ita）＊l*l*l/6+N*cb*cb＊cb*cb*cb＊（5—18＊t＊t+t＊t*t＊t+14＊ita＊ita-58＊ita＊ita*t＊t)＊l*l*l＊l*l/120）；

y=y+500000

}

/＊＊＊＊＊＊*＊*＊＊*＊*高斯反算函数***＊****＊＊*****

第 13 页

输入： double a ，f 椭球参数， x，y为高斯平面坐标，L0为中央子午线的经度； B,L为大地坐标，单位为弧度

＊返回:none

＊*＊＊＊＊＊＊*＊＊＊＊**＊*＊*＊**＊＊＊*＊＊＊/

void gaosibackward(double a,double f,double x，double y,double L0,double＆B,double &L)

｛

double b, c,e1, e2//短半轴，极点处的子午线曲率半径，第一偏心率，第二偏心率

double Bf,itaf,tf,Nf，Mf,Wf；

double l

第 14 页

double m[5］,n[5］

y=y-500000

//计算一些基本常量

｛

b=a＊（1—f）；

e1=sqrt(a*a—b*b）/a

e2=sqrt(a*a-b＊b)/b

c=a*a/b

m［0］=a＊（1-e1＊e1)

m[1]=3*（e1*e1*m［0］)/2.0；

m[2]=5*（e1＊e1*m［1］)/4。0；

第 15 页

m[3］=7*(e1＊e1*m[2］）/6。0；

m[4]=9*(e1＊e1＊m［3］)/8.0；

n［0]=m［0］+m［1]/2+3＊m［2]/8+5*m［3]/16+35*m［4］/128；

n[1]=m［1］/2+m[2]/2+15＊m［3]/32+7*m[4］/16

n[2］=m［2］/8+3*m［3］/16+7＊m[4]/32；

n[3]=m［3]/32+m[4］/16

n［4]=m［4]/128；

｝

//计算Bf

第 16 页

｛

double Bf1,Bfi0，Bfi1，FBfi；

Bf1=x/n[0]

Bfi0=Bf1

Bfi1=0

FBfi=0；

int num=0；

do

｛

num=0

FBfi=0.0-n[1］＊sin（2＊Bfi0）/2.0+n[2］*s

第 17 页

in(4*Bfi0）/4。0-n［3]*sin(6*Bfi0）/6。0；

Bfi1=（x—FBfi)/n［0]

if (fabs（Bfi1—Bfi0）>(Pi*pow(10.0,-8）/（36*18）)）

{

num=1；

Bfi0=Bfi1；

｝

｝while（num==1）

Bf=Bfi1；

}

tf=tan（Bf）；

第 18 页

Wf=sqrt（1-e1*e1＊sin（Bf）＊sin（Bf)）

Nf=a/Wf；

Mf=a＊（1—e1*e1)/(Wf＊Wf＊Wf)；

itaf=e2*cos（Bf)

B=Bf-tf＊y＊y/（2*Mf＊Nf）+tf*(5+3*tf＊tf+itaf＊itaf—9＊itaf*itaf*tf＊tf）*pow（y，4)/(24＊Mf＊pow(Nf，3））—tf＊（61+90*tf*tf+45*pow（tf,4）)*pow（y,6）/（720*Mf*pow（Nf，5)）

l=y/（Nf*cos(Bf））-(1+2*tf＊tf+itaf＊itaf)＊pow（y,3)/（6＊pow(Nf,3)＊cos（Bf))+(5+28＊

tf*tf+24＊pow(tf,4）+6*itaf*itaf+8*itaf＊itaf*tf*tf)*pow(y，5）/（120＊pow(Nf，5)*

第 19 页

cos(Bf）)

L=l+L0

｝

2014-5—22

’输入: double a ，f 椭球参数，B，L为大地坐标，L0为中央子午线的经度，单位为弧度，x，y为高斯平面坐标,y加上了常量

Private Function gaosiforward（ByVal a As Double, ByVal f As Double, ByVal B As Double，ByVal L As Double，ByVal L0 As Double）As Double（)

Dim x， y， xy(2) As Double

Dim bb, c, e1， e2 As Double’短半轴，极点

第 20 页

处的子午线曲率半径，第一偏心率，第二偏心率

Dim ll, W, N， M， daihao As Double’W为常用辅助函数，N为子午圈曲率半径，M为卯酉圈曲率半径

Dim xx As Double'子午线弧长，高斯投影的坐标

Dim ruo， ita， sb, cb, t As Double

Dim mm（5), nn（5) As Double

bb = a * （1 — f）

e1 = Math。Sqrt(a ＊ a - bb ＊ bb） / a

e2 = Math.Sqrt（a ＊ a — bb ＊ bb） / bb

c = a ＊ a / bb

第 21 页

mm（0) = a ＊（1 - e1 ＊ e1）

mm(1） = 3 ＊（e1 * e1 ＊ mm(0)） / 2.0

mm(2) = 5 * （e1 ＊ e1 ＊ mm（1)) / 4.0

mm（3) = 7 ＊（e1 ＊ e1 * mm(2)） / 6。0

mm（4） = 9 * （e1 ＊ e1 * mm（3）) / 8.0

nn(0) = mm(0） + mm（1) / 2 + 3 ＊ mm（2) / 8 + 5 * mm（3) / 16 + 35 ＊ mm（4） / 128nn(1) = mm(1) / 2 + mm(2） / 2 + 15 ＊ mm（3) / 32 + 7 ＊ mm（4） / 16

nn（2) = mm(2） / 8 + 3 ＊ mm（3) / 16 + 7 ＊ mm（4) / 32

nn(3) = mm（3) / 32 + mm(4） / 16

nn（4） = mm（4) / 128

第 22 页

xx = nn（0) ＊ B — Sin（B) * Cos（B）＊（（nn(1） - nn(2） + nn(3)) + (2 ＊ nn（2) - (16 ＊ nn(3） / 3。0）） * Sin（B）＊ Sin（B） + 16 * nn(3） * Pow(Sin（B), 4） / 3.0）

ll = L — L0 ’弧度

ita = e2 ＊ Cos（B)

sb = Sin（B)

cb = Cos(B)

W = Sqrt（1 - e1 ＊ e1 * sb * sb）

N = a / W

t = Tan(B）

ruo = (180 / PI） * 3600

第 23 页

x = (xx + N ＊ sb ＊ cb * ll ＊ ll / 2 + N ＊ sb ＊ cb ＊ cb ＊ cb * （5 - t ＊ t + 9 ＊ ita ＊ ita + 4 * ita ＊ ita * ita * ita) * ll ＊ ll * ll * ll / 24 + N * sb ＊ cb ＊ cb * cb ＊ cb ＊ cb * (61 — 58 ＊ t * t + t * t * t * t) * ll * ll ＊ ll * ll ＊ ll * ll / 720)

第 24 页

百度文库

搜索

百度文库10亿海量资料，查找管理一应俱全

打开APP

#网上搜来的

# 高斯坐标转经纬度算法 # B=大地坐标X # C=大地坐标Y  喊誉# IsSix=6度态裤带或3度带

import math

def GetLatLon2(B, C,IsSix):

#带号

D = math.trunc( C/ 1000000)

#中央经线(单位：弧度)

K = 0

if IsSix:

K = D * 6 - 3 #6度带计算

else:

K = D * 3 #3度带计算郑闭段

L = B/(6378245*(1-0.006693421623)*1.0050517739)

M = L +(0.00506237764 * math.sin(2*L)/2-0.00001062451*math.sin(4*L)/4+0.0000002081*math.sin(6*L)/6)/1.0050517739

N = L +(0.00506237764 * math.sin(2*M)/2-0.00001062451*math.sin(4*M)/4+0.0000002081*math.sin(6*M)/6)/1.0050517739

O = L +(0.00506237764 * math.sin(2*N)/2-0.00001062451*math.sin(4*N)/4+0.0000002081*math.sin(6*N)/6)/1.0050517739

P = L +(0.00506237764 * math.sin(2*O)/2-0.00001062451*math.sin(4*O)/4+0.0000002081*math.sin(6*O)/6)/1.0050517739

Q = L +(0.00506237764 * math.sin(2*P)/2-0.00001062451*math.sin(4*P)/4+0.0000002081*math.sin(6*P)/6)/1.0050517739

R = L +(0.00506237764 * math.sin(2*Q)/2-0.00001062451*math.sin(4*Q)/4+0.0000002081*math.sin(6*Q)/6)/1.0050517739

S = math.tan(R)

T = 0.006738525415*(math.cos(R))**2

U = 6378245/math.sqrt(1-0.006693421623*(math.sin(R))**2)

V = 6378245*(1-0.006693421623)/(math.sqrt((1-0.006693421623*(math.sin(R))**2)))**3

W = 5+3*S**2+T-9*T*S**2

X = 61+90*S**2+45*S**4

Y = 1+2*S**2+T**2

Z = 5+28*S**2+24*S**4+6*T+8*T*S**2

Lat= (180/math.pi)*(R-(C-D*1000000-500000)**2*S/(2*V*U)+(C-D*1000000-500000)**4*W/(24*U**3*V)-(C-D*1000000-500000)**6*X/(7200*U**5*V))

Lon= (180/math.pi)*(C-D*1000000-500000)*(1-(C-D*1000000-500000)**2*Y/(6*U**2)+(C-D*1000000-500000)**4*Z/(120*U**4))/(U*math.cos(P))

Lat = Lat

Lon = K + Lon

return (Lon, Lat)