高斯投影

我们知道椭球面是处理测量问题的基准面，测量中需要把球面元素投影化算到平时使用的平面。我们一般选择高斯投影来处理。

一、高斯投影的概念

设想一个横圆柱体把地球椭球体套在里面,使横轴圆柱的轴心通过球的中心,球面上一根子午线（经线）与横轴圆柱面相切。该子午线在圆柱面上的投影为一直线,赤道面与圆柱面的交线是一条与该子午线投影垂直的直线。将横圆柱面展开成平面,由这两条正交直线就构成高斯-克吕格平面。这个过程就是高斯投影。

二、高斯投影特点

1，中央子午线投影后是一条不变形的直线。

2，无角度变形，图形保持相似（等角或者正形投影）。

3，离中央子午线越远，变形越大。

三，高斯投影分带

按一定经差将地球椭球面划分成若干投影带,这是高斯投影中限制长度变形的最有效方法。分带时既要控制长度变形使其不大于测图误差，又要使带数不致过多以减少换带计算工作，据此原则将地球椭球面沿子午线划分成经差相等的瓜瓣形地带,以便分带投影 。一般有有六度带、三度带两种。显然分带愈多，各带包含的范围愈小，变形也就愈小。由于分带投影后。各投影带有自己的坐标轴和原点。从而形成各自独立的坐标系。这样，在相邻两带的点分别属于两个不同的坐标系，在工程中往往要将不同的坐标系化成同一坐标系，这就要进行相邻带之间的坐标换算。为了减少换带计算，分带不宜过多。

四，坐标换带计算

1，坐标换带的方法

（1）直接换带计算法（根据换带标计算）

（2）间接计算法（根据高斯投影正反算公式计算）

换带计算目前广泛采用高斯投影坐标正反算方法，他适用于任何情况下的换带计算工作。这种方法的程序是：首先将某投影带的已知平面坐标(x1,y1 )，按高斯投影坐标反算公式求得其大地坐标（B,L），然后根据纬度B和对于所选定的中央子午线的经差,按高斯投影坐标正算公式求其在选定的投影带的平面坐标（x2,y2）。

2，坐标换带的应用

（1）平面控制网跨带，平差时需要在同一个带计算。

（2）不同带之间转换（6度转换3度， 6度转换1.5度等）。

（3）跨带测图时，用到另外一带三角点作为控制，需要换带计算。

我有VB的，自己很多年前写的，一直用，但是正算－>反算－>正算后，Y坐标与原来的差了0.5－0.7mm，不知道怎么回事，这两年工作忙也没有时间再深究，但是这样的计算精度做控制足够了，如果楼主或是者是哪位同仁见此贴能顺便把这个问题解决了，咱们就一起进步了！代码如下：

'高斯坐标正算

Private Sub DadiZs()

Dim t As Double, Itp As Double, X0 As Double, N As Double, L0 As Double

Dim V As Double, ll As Double, W As Double, M As Double

Lat = Radian(Lat)

Lon = Radian(Lon)

L0 = Radian(Lo)

If Tq = 0 Then

a = 6378245 '54椭球参数

b = 6356863.01877305

ep = 0.006693421622966

ep1 = 0.006738525414683

f = (a - b) / a

c = a ^ 2 / b

d = b ^ 2 / a

X0 = 111134.8611 * (Lat * 180# / Pi) - (32005.7799 * Sin(Lat) + 133.9238 * (Sin(Lat)) ^ 3 + 0.6973 * (Sin(Lat)) ^ 5 + 0.0039 * (Sin(Lat)) ^ 7) * Cos(Lat)

'X0 = 111134.8611 * (Lat * 180# / Pi) - (32005.7798 * Sin(Lat) + 133.9238 * (Sin(Lat)) ^ 3 + 0.6972 * (Sin(Lat)) ^ 5 + 0.0039 * (Sin(Lat)) ^ 7) * Cos(Lat)

Else

a = 6378140 '75椭球参数

b = 6356755.28815753

ep = 0.006694384999588

ep1 = 0.006739501819473

f = (a - b) / a

c = a ^ 2 / b

d = b ^ 2 / a

X0 = 111133.0047 * (Lat * 180 / Pi) - (32009.8575 * Sin(Lat) + 133.9602 * (Sin(Lat)) ^ 3 + 0.6976 * (Sin(Lat)) ^ 5 + 0.0039 * (Sin(Lat)) ^ 7) * Cos(Lat)

End If

ll = Lon - L0

t = Tan(Lat)

Itp = ep1 * Cos(Lat) ^ 2

W = Sqr(1 - ep * Sin(Lat) ^ 2)

V = Sqr(1 + ep1 * Cos(Lat) ^ 2)

M = c / V ^ 3

N = a / W

'x = X0 + N * t * (Cos(Lat)) ^ 2 * ll ^ 2 / 2 + N * t * (5 - t * t + 9 * Itp + 4 * Itp * Itp) * (Cos(Lat)) ^ 4 * ll ^ 4 / 24 + N * t * (61 - 58 * t ^ 2 + t ^ 4 + 270 * Itp - 330 * t ^ 2 * Itp) * (Cos(Lat)) ^ 6 * ll ^ 6 / 720 + N * t * (1385 - 3111 * t ^ 2 + 543 * t ^ 4 - t ^ 6) * Cos(Lat) ^ 8 * ll ^ 8 / 40320

x = X0 + N * t * (Cos(Lat)) ^ 2 * ll ^ 2 / 2 + N * t * (5 - t * t + 9 * Itp ^ 2 + 4 * Itp ^ 4) * (Cos(Lat)) ^ 4 * ll ^ 4 / 24 + N * t * (61 - 58 * t ^ 2 + t ^ 4 + 270 * Itp ^ 2 - 330 * t ^ 2 * Itp ^ 2) * (Cos(Lat)) ^ 6 * ll ^ 6 / 720 + N * t * (1385 - 3111 * t ^ 2 + 543 * t ^ 4 - t ^ 6) * Cos(Lat) ^ 8 * ll ^ 8 / 40320

y = N * Cos(Lat) * ll + N * (1 - t * t + Itp) * (Cos(Lat)) ^ 3 * ll ^ 3 / 6 + N * (5 - 18 * t * t + t ^ 4 + 14 * Itp - 58 * Itp * t * t) * (Cos(Lat)) ^ 5 * ll ^ 5 / 120 + N * (61 - 479 * t ^ 2 + 179 * t ^ 4 - t ^ 6) * Cos(Lat) ^ 7 * ll ^ 7 / 5040

r = Sin(Lat) * ll + Sin(Lat) * (Cos(Lat)) ^ 2 * ll ^ 3 * (1 + 3 * Itp + 2 * Itp ^ 2) / 3 + Sin(Lat) * (Cos(Lat)) ^ 4 * ll ^ 5 * (2 - t * t) / 15

r = Degree(r)

y = y + 500000#

End Sub

'高斯反算

Private Sub DadiFs()

Dim t As Double, Itp As Double, X0 As Double, Bf As Double, N As Double

Dim v As Double, ll As Double, W As Double, M As Double, L0 As Double

L0 = Radian(Lo)

X0 = x * 0.000001

y = y - 500000#

If Tq = 0 Then

a = 6378245 '54椭球参数

b = 6356863.01877305

ep = 0.006693421622966

ep1 = 0.006738525414683

f = (a - b) / a

c = a ^ 2 / b

d = b ^ 2 / a

If X0 <3 Then

Bf = 9.04353301294 * X0 - 0.00000049604 * X0 ^ 2 - 0.00075310733 * X0 ^ 3 - 0.00000084307 * X0 ^ 4 - 0.00000426055 * X0 ^ 5 - 0.00000010148 * X0 ^ 6

ElseIf X0 <6 Then

Bf = 27.11115372595 + 9.02468257083 * (X0 - 3) - 0.00579740442 * (X0 - 3) ^ 2 - 0.00043532572 * (X0 - 3) ^ 3 + 0.00004857285 * (X0 - 3) ^ 4 + 0.00000215727 * (X0 - 3) ^ 5 - 0.00000019399 * (X0 - 3) ^ 6

End If

Else

a = 6378140 '75椭球参数

b = 6356755.28815753

ep = 0.006694384999588

ep1 = 0.006739501819473

f = (a - b) / a

c = a ^ 2 / b

d = b ^ 2 / a

If X0 <3 Then

Bf = 9.04369066313 * X0 - 0.00000049618 * X0 ^ 2 - 0.00075325505 * X0 ^ 3 - 0.0000008433 * X0 ^ 4 - 0.00000426157 * X0 ^ 5 - 0.0000001015 * X0 ^ 6

ElseIf X0 <6 Then

Bf = 27.11162289465 + 9.02483657729 * (X0 - 3) - 0.00579850656 * (X0 - 3) ^ 2 - 0.00043540029 * (X0 - 3) ^ 3 + 0.00004858357 * (X0 - 3) ^ 4 + 0.00000215769 * (X0 - 3) ^ 5 - 0.00000019404 * (X0 - 3) ^ 6

End If

End If

Bf = Bf * Pi / 180#

t = Tan(Bf)

Itp = ep1 * Cos(Bf) ^ 2

W = Sqr(1 - ep * Sin(Bf) ^ 2)

v = Sqr(1 + ep1 * Cos(Bf) ^ 2)

M = c / v ^ 3

N = a / W

Lat = Bf - 0.5 * v ^ 2 * t * ((y / N) ^ 2 - (5 + 3 * t * t + Itp - 9 * Itp * t * t) * (y / N) ^ 4 / 12 + (61 + 90 * t * t + 45 * t ^ 4) * (y / N) ^ 6 / 360)

ll = ((y / N) - (1 + 2 * t * t + Itp) * (y / N) ^ 3 / 6 + (5 + 28 * t * t + 24 * t ^ 4 + 6 * Itp + 8 * Itp * t * t) * (y / N) ^ 5 / 120) / Cos(Bf)

r = y * t / N - y ^ 3 * t * (1 + t * t - Itp) / (3 * N ^ 3) + y ^ 5 * t * (2 + 5 * t * t + 3 * t ^ 4) / (15 * N ^ 5)

Lat = Degree(Lat)

Lon = Degree(L0 + ll)

r = Degree(r)

End Sub

有了正反算，换带也就完成了！

用到的子程序：

Public Const Pi = 3.14159265358979, p = 206264.806

Public Cktq As String

'角度化弧度

Public Function Radian(a As Double) As Double

Dim Ro As Double

Dim c As Double

Dim Fs As Double

Dim Ib As Integer

Dim Ic As Integer

If a <0 Then a = -a: t = 1

Ro = Pi / 180#

Ib = Int(a)

c = (a - Ib) * 100#

Ic = Int(c + 0.000000000001)

Fs = (c - Ic) * 100#

If t = 1 Then Radian = -(Ib + Ic / 60# + Fs / 3600#) * Ro Else Radian = (Ib + Ic / 60# + Fs / 3600#) * Ro

End Function

'弧度化角度

Public Function Degree(a As Double) As Double

Dim Bo As Double

Dim Fs As Double

Dim Im As Integer

Dim Id As Integer

If a <0 Then a = -a: t = 1

Bo = a

Call DMS(Bo, Id, Im, Fs)

If t = 1 Then Degree = -(Id + Im / 100# + Fs / 10000#) Else Degree = Id + Im / 100# + Fs / 10000#

End Function

Public Sub DMS(a As Double, Id As Integer, Im As Integer, Fs As Double)

Dim Bo As Double

Dim c As Double

c = a

c = 180# / Pi * c

Id = Int(c)

Bo = (c - Id) * 60

Im = Int(Bo)

Fs = (Bo - Im) * 60

End Sub

'取位计算

Public Function Qw(a As Double, Ws As Integer) As Double

Qw = Int(a * 10 ^ Ws + 0.5) / 10 ^ Ws

End Function

另外,站长团上有产品团购,便宜有保证

高斯投影的编程研究目的是为了使程序条理更清晰，数据的复用性更强。根据相关资料显示，高斯投影的编程加深理解高斯—克吕格投影的实质，掌握高斯—克吕格投影直角坐标的运用，理解通用坐标和自然坐标值的关系。

---