全部栏目

首页 旅游资讯 线路攻略 景点大全 国内游 境外游 美食特产
您的当前位置：首页高斯投影正反算公式

高斯投影正反算公式

2021-09-18 来源：乌哈旅游


高斯投影坐标正反算

一、基本思想：

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

二、计算模型：

基本椭球参数：

椭球长半轴a

椭球扁率f

椭球短半轴：ba(1f)

a2b2ea椭球第一偏心率 ：

椭球第二偏心率 ：

ea2b2b 高斯投影正算公式：此公式换算的精度为0.001m

xXNN232244sinBcosBlsimBcosB(5t94)l22244N5246sinBcosB(6158tt)l7206

yNN3223cosBlcosB(1t)l63N5242225cosB(518tt1458t)l1205 其中：角度都为弧度

B为点的纬度，lLL0，L为点的经度，L0为中央子午线经度；

N为子午圈曲率半径，Na(1esinB)；

2212ttanB；

2e2cos2B

1803600

1616Xa0BsinBcosB(a2a4a6)(2a4a6)sin2Ba6sin4B33 其中X为子午线弧长：

a0,a2,a4,a6,a8为基本常量，按如下公式计算：

m23535ammmm804602816128am2m415m7m682223216m437m6m8a481632m6m8a63216am88128

m0,m2,m4,m6,m8为基本常量，按如下公式计算：

379m0a(1e2);m2e2m0;m45e2m2;m6e2m4;m8e2m6268；

高斯投影反算公式：此公式换算的精度为0.0001’’.

tf2MfNftf24MfN3fBBftfy253t2f2242f9ftfy720MfN5f46y6190t2f45tfyyy32l12t2ff3NfcosBf6NfcosBfy5422528t22f24tf6f8ftf5120NfcosBfLlL0

其中：

L0为中央子午线经度。

Bf为底点纬度，也就是当xX时的子午线弧长所对应的纬度。按照子午线弧长公式：

Xaa20B2sin2Ba44sin4Ba6a6sin6B88sin8B，迭代进行计算；

初始开始时设：

B1fXa0

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

Bi1fX(XF(Bif))a0F(Bif)a22sin2Bia4aaf4sin4Bif66sin6Bif88sin8Bif

重复迭代至

Bi1fBif为止。

Nfa(1e2sin2B1f)2；

M222fa(1e)(1esinBf)32

tftanBf；

222fecosBf

海福特椭球(1910) 椭球

a=6378388m b=6356911.9461279m α=0.33670033670我国52年以前基准

克拉索夫斯基椭球(1940 Krassovsky) 北京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

三、程序代码函数：

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

输入 : 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为卯酉圈曲率半径

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;

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 把改成了

}

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

{

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;

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;

}

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

输入 : 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;

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;

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

{

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]*sin(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);

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)*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 '短半轴，极点处的子午线曲率半径，第一偏心率，第二偏心率

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

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) / 128

nn(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

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

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)

y = (N * cb * ll + N * cb * cb * cb * (1 - t * t + ita * ita) * ll * ll * ll / 6 + N * cb * cb * cb * cb * cb * (5 - 18 * t * t + t * t * t * t + 14 * ita * ita - 58 * ita * ita * t * t) * ll * ll * ll * ll * ll / 120)

y = y + 500000

xy(0) = x

xy(1) = y

Return xy

End Function

Private Function gaosibackward(ByVal a As Double, ByVal f As Double, ByVal x As Double, ByVal y As Double, ByVal L0 As Double) As Double()

Dim b, l, bl(2) As Double

Dim bb, c, e1, e2 As Double '短半轴，极点处的子午线曲率半径，第一偏心率，第二偏心率

Dim Bf, itaf, tf, Nf, Mf, Wf As Double

Dim ll As Double

Dim m(5), n(5) As Double

y = y - 500000

bb = a * (1 - f)

e1 = Sqrt(a * a - bb * bb) / a

e2 = Sqrt(a * a - bb * bb) / bb

c = a * a / bb

m(0) = a * (1 - e1 * e1)

m(1) = 3 * (e1 * e1 * m(0)) / 2.0

m(2) = 5 * (e1 * e1 * m(1)) / 4.0

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

Dim Bf1, Bfi0, Bfi1, FBfi As Double

Bf1 = x / n(0)

Bfi0 = Bf1

Bfi1 = 0

FBfi = 0

Dim num As Integer = 0

Do

num = 0

FBfi = 0.0 - n(1) * Sin(2 * Bfi0) / 2.0 + n(2) * Sin(4 * Bfi0) / 4.0 - n(3) * Sin(6 * Bfi0) / 6.0

Bfi1 = (x - FBfi) / n(0)

If (Abs(Bfi1 - Bfi0) > (PI * Pow(10.0, -8) / (36 * 18))) Then

num = 1

Bfi0 = Bfi1

End If

Loop While num = 1

Bf = Bfi1

tf = Tan(Bf)

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))

ll = 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) * Cos(Bf))

l = ll + L0

bl(0) = hdtod(b)

bl(1) = hdtod(l)

Return bl

End Function

'度转换为弧度

Private Function dtohd(ByVal d As Double) As Double

Dim hd As Double

hd = d * PI / 180

Return hd

End Function

'弧度转换为度

Private Function hdtod(ByVal hd As Double) As Double

Dim d As Double

d = hd * 180 / PI

Return d

End Function

Private Function dfmtod(ByVal dfm As Double) As Double

Dim dfm2, dfm3(), du, fen, miao, miao1, miao2 As String

Dim duf As Double

dfm2 = dfm.ToString

dfm3 = dfm2.Split(\".\")

du = dfm3(0)

fen = dfm3(1).Substring(0, 2)

miao1 = dfm3(1).Substring(2, 2)

miao2 = dfm3(1).Substring(4)

miao = miao1 & \".\" & miao2

duf = Convert.ToInt32(du) + Convert.ToInt32(fen) / 60.0 + Convert.ToDouble(miao) / 3600

Return duf

End Function

'国家坐标系椭球参数2000,84,54,80

a2 = 6378137

f2 = 1 / 298.257222101

a84 = a2

f84 = 1 / 298.257223563

a54 = 6378245

f54 = 0.33523298692

a80 = 6378140

f80 = 0.0033528131778

因篇幅问题不能全部显示，请点此查看更多更全内容

查看全文