用实现测量坐标转换系统

1、坐标转换系统（VB ）东华理工大学Theory北京54高斯坐标转换西安80高斯坐标转换系统1.0版-匸Fj-Si-Ifi冋需点赴理打幵处共点戲爵njFwmunm厂UF-一MhinBOJS点处理* J吒址标系I i *a0D坐棣系高聊坐标k高斷坐标F高厮生际y酬h1 1桃h多点处理结果二甲央子午卿经度R尹幵弐址空矗|町-54坐様点导Xi-yTH0413762520 59409497513 333MG403TB2M9 14123媳7760 964館34-换3TT822< 127IB49fcT30 43803343TT9033 3791996706 3191434G37377C236.3593

2、0495$S3 0451733.鬲3778505 76197496250 39T0933d3117197 8295987 93TC34537T77S7 35212136603 96T393<G333T7B420.6632149M5S 7232334.«323776392 19456<97141, J24SB34的1377513 fC&t帕7蛙7羽嗨I34<303T7B09O 110004 歯KE $312134377356. &1256498753.50420<T 1专t曽qgi rffvnncrrw? 8仲BJ-54 加 i inSO 草点

3、it理监20换系核心源代码:'高斯坐标转换成大地坐标过程Public Sub GausReverse(a As Double, f As Double, x() As Double, y() As Double, RB() AsDouble, RL() As Double, k As In teger)Dim i As In teger, fxb As Double, fxbl As Double, fybl As Doublee = Sqr(2 * f - f A 2)C = a / Sqr(1 - e a 2)e2 = e / Sqr(1 - e a 2)beitaO = 1 -

4、(3 / 4) * e2 a 2 + (45 / 64) * e2 a 4 - (175 / 256) * e2 a 6 + (11025 / 16384) * e2 a 8beita2 = beitaO - 1beita4 = (15 / 32) * e2 a 4 - (175 / 384) * e2 a 6 + (3675 / 8192) * e2 a 8beita6 = (-35 / 96) * e2 a 6 + (735 / 2048) * e2 a 8beita8 = (315 / 1024) * e2 a 8Fori = 1 To kBO = x(i) / (C * beitaO)

5、a1 = a * Cos(BO) / Sqr(1 - e a 2 * (Si n(B0) a 2)lO = y(i) / a1Dofxb = 0fxbl = 0fybl = 0t = Tan (BO)yita = e2 * Cos(BO)n = a / Sqr(1 - e a 2 * (Si n(B0) a 2)a2 = (1 / 2) * n * Sin(B0) * Cos(BO)a3 = (1 / 6) * n * (Cos(BO) a 3 * (1 - t a 2 + yita a 2)a4 = (1 / 24) * n * Sin (BO) * (Cos(BO) a 3 * (5 -

6、t a 2 + 9 * yita a 2 +4*yita a 4)a5 = (1 / 120) * n * (Cos(BO) a 5 * (5 - 18 * t a 2 + t a 4 + 14 * yita a2-58 *yita a 2 * t a 2)a6 = (1 / 720) * n * Sin (BO) * (Cos(BO) a 5 * (61 - 58 * t a 2 + t a 4)fxb = fxb + (C * beita6 + C * beita8 * (Cos(BO) a 2) * (Cos(BO) a 2fxb = (fxb + C * beita4) * (Cos(

7、BO) A 2fxb = (fxb + C * beita2) * Sin (BO) * Cos(B0)fxbl = a2 * 10 a 2 + a4 * 10 a 4 + a6 * 10 a 6fybl = a3 * 10 a 3 + a5 * 10 a 5RB(i) = (x(i) - fxb - fxbl) / (C * beita0)al = (a * Cos(RB(i) / Sqr(1 - e a 2 * (Si n( RB(i) a 2)RL(i) = (y(i) - fyb1) / alIf Abs(RB(i) - B0) <= 0.0000000001 And Abs(R

8、L(i) - 10) <= 0.0000000001 ThenRL(i) = zrl + 10Exit DoE1seB0 = RB(i)10 = RL(i)End IfLoopNext iEnd Sub'大地坐标B , L转换为高斯坐标 x,y的过程Public Sub BLHGaus(RB() As Double, RL() As Double, GX() As Double, GY() As Double, a AsDouble, f As Double, k As In teger)Dim 10 As Double, fxb As Double, gxbl As Doubl

9、e, fybl As Doublee = Sqr(2 * f - f a 2)C = a / Sqr(1 - e a 2)e2 = e / Sqr(1 - e a 2)beita0 = 1 - (3 / 4) * e2 a 2 + (45 / 64) * e2 a 4 - (175 / 256) * e2 a 6 + (11025 / 16384) * e2 a 8beita2 = beita0 - 1beita4 = (15 / 32) * e2 a 4 - (175 / 384) * e2 a 6 + (3675 / 8192) * e2 a 8beita6 = (-35 / 96) *

10、e2 a 6 + (735 / 2048) * e2 a 8beita8 = (315 / 1024) * e2 A 8Fori = 1 To k10 = RL(i) - zrln = a / Sqr(1 - e a 2 * (Si n(RB(i) a 2)t = Tan (RB(i)yita = e2 * Cos(RB(i)al = (a* Cos(RB(i) / Sqr(1 - e a 2 * (Si n( RB(i) a 2)a2 = (1/ 2) * n * Sin(RB(i) * Cos(RB(i)a3 = (1/ 6) * n * (Cos(RB(i) a 3 * (1 - t a

11、 2 + yita a 2)a4 = (1/ 24) * n * Sin (RB(i) * (Cos(RB(i) a 3 * (5 - t a2 +9 *yitaa2 +4 * yitaa4)a5 = (1/ 120) * n * (Cos(RB(i) a 5 * (5 - 18 * t a 2 + ta4 +14* yitaa 2- 58 *yitaa 2 * t a2)a6 = (1 / 720) * n * Sin( RB(i) * (Cos(RB(i) a 5 * (61 - 58 * t a 2 + t a 4)fxb = 0fxbl = 0fybl = 0fxb = fxb + (

12、C * beita6 + C * beita8 * (Cos(RB(i) a 2) * (Cos(RB(i) a 2fxb = (fxb + C * beita4) * (Cos(RB(i) a 2fxb = (fxb + C * beita2) * Sin(RB(i) * Cos(RB(i)fxbl = a2 * 10 a 2 + a4 * 10 a 4 + a6 * 10 a 6fybl = a3 * 10 a 3 + a5 * 10 a 5GX(i) = C * beita0 * RB(i) + fxb + fxblGY(i) = a1 * l0 + fyblNext iEnd Sub&

13、#39;三维直角坐标 XYZ转换成大地坐标 BLH过程Public Sub BLHXYZ1(SX() As Double, SY() As Double, SZ() As Double, RB() As Double, RL()As Double, RH() As Double, k As In teger, a As Double, f As Double)Dim i As In tegerDim NO As Double, H0 As Double, BO As Double, sb As Double, Ni As Doublee = Sqr(2 * f - f A 2)C = a /

14、Sqr(1 - e a 2)e2 = e / Sqr(1 - e a 2)sb = a * Sqr(1 - e a 2)pi = 4 * Atn(1)Fori = 1 To kRL(i) = Atn (Abs(SY(i) / SX(i)If SY(i) > 0 And SX(i) < 0 The nRL(i) = pi - RL(i)End IfNO = aHO = Sqr(SX(i) a 2 + SY(i) a 2 + SZ(i) a 2) - Sqr(a * sb)BO = Atn(SZ(i) / (Sqr(SX(i) a 2 + SY(i) a 2) * (1 - e a 2

15、 * NO / (NO + HO)DoNi = a / Sqr(1 - e a 2 * (Si n( BO) a 2)RH(i) = Sqr(SX(i) a 2 + SY(i) a 2) / Cos(BO) - NiRB(i) = Atn( SZ(i) / (Sqr(SX(i) a 2 + SY(i) a 2) * (1 - e a 2 * Ni / (Ni + RH(i)If Abs(RB(i) - BO) < O.OOOOOOOOO1 And Abs(RH(i) - HO) < O.OOOOOOOOO1 ThenExit DoElseBO = RB(i)HO = RH(i)En

16、d IfLoopNext iEnd Sub'最小二乘法求解七参数布尔萨模型Public Sub SloveBuersa(XD() As Double, XG() As Double, R() As Double, k As Integer)Dim a0 As Double, a1 As Double, a2 As DoubleFori = 1 To ka0 = XG(i, 1)a1 = XG(i, 2)a2 = XG(i, 3)Call Cjuzhe n(aO, a1, a2, G()Forj = 1 To 3For s = 1 To 7EC(3 * (i - 1) + j, s) =

17、 G(j, s)Next sNext jNext iFori = 1 To kForj = 1 To 3LC(3 * (i - 1) + j) = XD(i, j) - XG(i, j)Next jNext iFor i = 1 To 3 * kForj = 1 To 7ET(j, i) = EC(i, j)Next jNext iFor i = 1 To 7Forj = 1 To 7CTC(i, j) = 0CTC(i, j) = CTC(i, j) + ET(i, s) * EC(s, j)Next sNext jNext iFor i = 1 To 7Forj = 1 To 7CTC1(

18、i, j) = CTC(i, j)Next jNext iCall Comm.Reverse(CTC(), 7)For i = 1 To 7CTL(i) = 0For j = 1 To 3 * kCTL(i) = CTL(i) + ET(i, j) * LC(j)Next jNext iFor i = 1 To 7R(i) = 0Forj = 1 To 7R(i) = R(i) + CTC(i, j) * CTL(j)Next jNext iEnd Sub'布尔萨模型系数矩阵一部分Public Sub Cjuzhe n(x As Double, y As Double, z As Do

19、uble, D() As Double)Fori = 1 To 3D(i, i) = 1Next iD(1,4) = xD(1,5) = 0D(1,6) = -zD(1,7) = yD(2, 4) = yD(2, 5) = zD(2, 6) = 0D(2, 7) = -xD(3, 4) = zD(3, 5) = -yD(3, 6) = xD(3, 7) = 0End Sub'使用七参数求解新坐标系下的坐标Public Sub UseBuersa(XD() As Double, XG() As Double, EC() As Double, R() As Double, k AsIn t

20、eger)Dim a1 As Double, a2 As Double, a3 As Double, CR(1000) As DoubleFori = 1 To ka1 = XG(i, 1)a2 = XG(i, 2)a3 = XG(i, 3)Call Cjuzhen(a1, a2, a3, G()Forj = 1 To 3For s = 1 To 7EC(3 * (i - 1) + j, s) = G(j, s)Next sNext jNext iFor i = 1 To 3 * kCR(i) = 0CR(i) = CR(i) + EC(i, j) * R(j)Next jNext iFori

21、 = 1 To kForj = 1 To 3XD(i, j) = XG(i, j) + CR(3 * (i - 1) + j)Next jNext iEnd Sub'弧度化为度分秒'弧度化成角度Public Sub RuJiao(ByVal rudu As Double, jiaodu As Double)Dim ja As In teger, jb As In teger, jc As Doublepi = 4 * Atn (1)jiaodu = rudu * 180 / pija = Fix(jiaodu)jb = Fix(jiaodu - ja) * 60)jc = (j

22、iaodu - ja) * 60 - jb) * 60jiaodu = ja + jb / 100 + jc / 10000End Sub矩阵求逆Public Sub Reverse(Ba, n%)Dim k%, K1%, j%, i%Dim C As Double, Aa(100, 200)For i = 1 To nAa(i, j) = Ba(i, j)Next jNext iFor i = 1 To nForj = 1 To nIf i = j The nAa(i, j + n) = 1ElseAa(i, j + n) = 0End IfNext jNext iFor k = 1 To nForj = k T