《大地测量学基础》里面的高斯投影正反算公式及换带计算的VB或者C语言编程

我有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

工程施工过程中，常常会遇到不同坐标系统间，坐标转换的问题。目前国内常见的转换有以下几种：1，大地坐标（BLH）对平面直角坐标（XYZ）；2，北京54全国80及WGS84坐标系的相互转换；3，任意两空间坐标系的转换。其中第2类可归入第三类中。所谓坐标转换的过程就是转换参数的求解过程。常用的方法有三参数法、四参数法和七参数法。以下对上述三种情况作详细描述如下：

1，大地坐标（BLH）对平面直角坐标（XYZ）

常规的转换应先确定转换参数，即椭球参数、分带标准（3度，6度）和中央子午线的经度。椭球参数就是指平面直角坐标系采用什么样的椭球基准，对应有不同的长短轴及扁率。一般的工程中3度带应用较为广泛。对于中央子午线的确定有两种方法，一是取平面直角坐标系中Y坐标的前两位*3，即可得到对应的中央子午线的经度。如x=3250212m，y=395121123m，则中央子午线的经度=39*3=117度。另一种方法是根据大地坐标经度，如果经度是在155.5~185.5度之间，那么对应的中央子午线的经度=（155.5+185.5）/2=117度，其他情况可以据此3度类推。

另外一些工程采用自身特殊的分带标准，则对应的参数确定不在上述之列。

确定参数之后，可以用软件进行转换，以下提供坐标转换的程序下载。

2，北京54全国80及WGS84坐标系的相互转换

这三个坐标系统是当前国内较为常用的，它们均采用不同的椭球基准。

其中北京54坐标系，属三心坐标系，大地原点在苏联的普而科沃，长轴6378245m，短轴6356863，扁率1/298.3；西安80坐标系，属三心坐标系，大地旅桐原点在陕西省径阳县永乐镇，长轴6378140m，短轴6356755，扁率1/298.25722101；WGS84坐标系，长轴6378137.000m，短轴6356752.314，扁率1/298.257223563。由于采用的椭球基准不一样，并且由于投影的局限性，使的全国各地并不存在一至的转换参数。对于这种转换由于量较大，有条件的话，一般都采用GPS联测已知点，应用GPS软件自动完成坐标的转换。当然若条件不许可，且有足够的重合点，也可以进行人工解算。详细方法见第三类。

3，任意两空间坐标系的转换拆羡坦

由于测量坐标系和施工坐标系采用不同的标准，要进行精确转换，必须知道至少3个重合点（即为在两坐标系中坐标均为已知的点。采用布尔莎模型进行求解。布尔莎公式:

对该公式进行变换等价得到：

解算这七个参数，至少要用到三个已知点（2个坐标系统的坐标都知道），采用间接平差模型进行解算：

其中： V 为残差矩阵

X 为未知七参数

A 为系数矩阵

解之：L 为闭合差

解得派神七参数后，利用布尔莎公式就可以进行未知点的坐标转换了，每输入一组坐标值，就能求出它在新坐标系中的坐标。 但是要想GPS观测成果用于工程或者测绘，还需要将地方直角坐标转换为大地坐标，最后还要转换为平面高斯坐标。

上述方法类同于我们的间接平差，解算起来较复杂，以下提供坐标转换程序，只需输入三个已知点的坐标即可求解出坐标转换的七个参数。如果已知点的数量较多，可以进行参数间的平差运算，则精度更高。

当已知点的数量只有两个时，我们可以采用简单变换法，此法较为方便易行，适于手算，只是精度受到一定的限制。

详细解算方程如下：

式中调x,y和x\'、y\'分别为新旧（或；旧新）网重合点的坐标，a、b、、k为变换参数，显然要解算出a、b、、k，必须至少有两个重合点，列出四个方程。

即可进行通常的参数平差，解求a、x、b、c、d各参数值。将之代人（3）式，可得各拟合点的残差（改正数）代人（2）式，可得待换点的坐标。

求出解算参数之后，可在Excel中，进行其余坐标的转换。

上次笔者用此法进行过80和54坐标的转换，由于当时没有多余的点可供验证和平差，所以转换精度不得而知，但转换之后各点的相对位置不变。估计，实际的转换误差应该是10m量级的。

还有一些情况是先将大地坐标转换 为直角坐标，然后进行相关转换。

A点Y坐标前2位是39，B点Y坐标前2位是29，表示A点属于第39度带，B点属于第29度带，而我国按6度带划分，在12带到23带，按3度带划分，在24带到45带，所以A、B两点均是按3度带高斯投影。那么A点对应的中央子午线是117°，B点对应的中央子午线是87°。由于A点和B点不是相邻的两个带，不能换带计算再求距离，只好先把平面坐标转换成大地经纬度坐标，再按照大地主体解算的方法求出AB间距离。高斯平面坐标换算成大地坐标结果：A（B-32°48′54.23″，L-116°19′57.44″），B（B-44°05′15.95″，L-88°07′14.60″），大地主题反算（具体原理你找本《老禅大地测量学》书籍看看）我选白塞尔方法，解算结散渗果：AB间大冲含脊地球面距离S=2743353.59480769m，大地方位角（从A到B方向）是305°21′09″。

---