VB高斯正反算

1、vb下高斯坐标变换的实现曾圣陈伟高斯- 克吕格 (gauss-kruger) 投影，是一种“等角横切圆柱投影”。德国数学家、物理学家、天文学家高斯（carl friedrich auss，1777 一 1855 ）于十九世纪二十年代拟定，后经德国大地测量学家克吕格 （johannes kruger ， 18571928）于 1912 年对投影公式加以补充，故名。设想用一个圆柱横切于球面上投影带的中央经线，按照投影带中央经线投影为直线且长度不变和赤道投影为直线的条件，将中央经线两侧一定经差范围内的球面正形投影于圆柱面. 高斯- 克吕格 (gauss-kruger) 投影是横轴墨卡托投影的变种，

2、高斯- 克吕格投影是“等角横切圆柱投影”，投影后中央经线保持长度不变，即比例系数为1。高斯- 克吕格投影的东伪偏移是500 公里，投影北伪偏移为零。高斯- 克吕格投影正解公式 ：已知（ b,l）值求解 (x,y) ，（原点纬度 0 ，中央经度 l0）上面公式中东纬偏移fe = 500000 米 （本程序中中只设计加500000米的常数，如还要带号的话 fe = 500000+ 带号 * 1000000 米）；高斯 - 克吕格投影比例因子 k0 = 1 高斯- 克吕格投影反解公式 ：已知(x,y) 求解（ b,l），（原点纬度 0 ，中央经度l0）a - 椭球体长半轴b - 椭球体短半轴f -

3、扁率e - 第一偏心率e - 第二偏心率n - 卯酉圈曲率半径r - 子午圈曲率半径b - 纬度， l - 经度，单位弧度 (rad) x - 纵直角坐标 , y - 横直角坐标，单位米 (m) 关于椭球体参数，我国常用的3 个椭球体参数如下：椭球体长半轴 a （米）短半轴 b（米）krassovsky （北京 54采用）6378245 6356863.0188 iag 75 （西安 80采用）6378140 6356755.2882 wgs 84 6378137 6356752.3142 程序实现 : 高斯正算求 x public function x(byval b#, byval l#,

4、 byval l0#) as double dim n#, t#, t2#, m#, m2#, ng2# dim s#, c# x = a1 * b + a2 * sin(2 * b) + a3 * sin(4 * b) + a4 * sin(6 * b)? 子午线弧长s = sin(b) c = cos(b) t = tan(b) t2 = t * t n = a / sqr(1 - e12 * s * s)卯酉圈曲率半径m = c * (l - l0) m2 = m * m ng2 = c * c * e12 / (1 - e12)? x = x + n * t * (0.5 + (5 -

5、 t2 + 9 * ng2 + 4 * ng2 * ng2) / 24# + (61 - 58 * t2 + t2 * t2) * m2 / 720#) * m2) * m2)? end function 高斯正算求 y public function y(byval b#, byval l#, byval l0#) as double dim n#, t#, t2#, m#, m2#, ng2# dim s#, c# s = sin(b) c = cos(b) t = tan(b) t2 = t * t n = a / sqr(1 - e12 * s * s) m = c * (l - l0

6、) m2 = m * m ng2 = c * c * e12 / (1 - e12) y = n * m * (1 + m2 * (1 - t2 + ng2) / 6# + m2 * (5 - 18 * t2 + t2 * t2 + 14 * ng2 - 58 * ng2 * t2) / 120#) y = y + y0 end function 高斯反算求 b（纬度 ) public function b(byval x#, byval y#) as double dim s#, c#, t#, t2#, n#, ng2#, v#, yn# dim preb0#, b0# dim eta#

7、y = y - y0 b0 = x / a1 do preb0 = b0 b0 = (x - (a2 * sin(2 * b0) + a3 * sin(4 * b0) + a4 * sin(6 * b0) / a1 if abs(b0 - preb0) 0.000000001 then exit do loop s = sin(b0) c = cos(b0) t = tan(b0) t2 = t * t n = a / sqr(1 - e12 * s * s) ng2 = c * c * e12 / (1 - e12) v = sqr(1 + ng2) yn = y / n b = b0 -

8、(yn * yn - (5 + 3 * t2 + ng2 - 9 * ng2 * t2) * yn * yn * yn * yn / 12# + (61 + 90 * t2 + 45 * t2 * t2) * yn * yn * yn * yn * yn * yn / 360#) * v * v * t / 2# end function 程序运行界面 : 检测数据 : （1）、已知在北京坐标系下中央子午线： 117度，纬度 b：28 度 32 分 14.5 秒，经度 l：116 度 54 分 12.3秒正算求解北方向 x、东方向 y值。得 : x= 3158054.118 米y= 490547.265 米（2）、已知在北京坐标系下中央子午线： 114度，北方向 x值： 2834561.381 米、东方向 y值：583726.735米反算求解 b、l 值。得 : b= 25 度 36分 56.1 秒l= 114 度 50 分 00.8 秒（3）、已知在北京坐标系下