高斯反算公式详解及代码实现（VB.Net）-CSDN博客

本文链接：https://blog.csdn.net/qq_34451885/article/details/144306242

高斯反算公式是将平面直角坐标（x，y）转换为地理坐标（经度L和纬度B）的数学公式。

已知高斯平面坐标（x，y）、椭球长半轴a、短半轴b（或扁率f），以及中央子午线经度 $L_{0}$ ，

具体公式如下（精度为0.0001″）：

其中，N为卯酉圈曲率半径：

$N=a(1-e^{^{2}}sin^{2}B)^{-\frac{1}{2}}$

注意，B要转为弧度制。

e为椭球体的第一偏心率： $e=\frac{\sqrt{a^{2}-b^{2}}}{a}$

$B_{f}$ 为利用子午线弧长X反算得到的大地纬度，单位为弧度， $B_{f}^{i}$ 为基于迭代法利用子线弧长X反算得到的大地纬度：

其中，X为自赤道至纬度B的子午线弧长，翻我上一篇高斯正算的文章有说过，这里贴结果：

$X=a_{0}B-\frac{a_{2}}{2}sin2B+\frac{a_{4}}{4}sin4B-\frac{a_{6}}{6}sin6B+\frac{a_{8}}{8}sin8B$

其中：

将 $B_{f}^{i}$ 反复迭代直至下式为止，以保证 Bf 精确至0.0001"：

$\left | B_{f}^{i+1} -B_{f}^{i}\right |<1\times 10^{-8}$

由迭代公式得到的 Bf 是以（弧度）为单位。

公式计算出的 L、B 为角度制。

用VB.net实现，代码如下：

 Private Function fansuan(ByVal za As Long, ByVal zf As Double, ByVal zL0 As Double, ByVal zx As Double, ByVal zy As Double, ByRef zL As Double, ByRef zB As Double) As Double
        'za为长半轴，zf为扁率，zl0为中央子午线
        Dim m0 As Double
        Dim m2 As Double
        Dim m4 As Double
        Dim m6 As Double
        Dim m8 As Double
        Dim a0 As Double
        Dim a2 As Double
        Dim a4 As Double
        Dim a6 As Double
        Dim a8 As Double
        Dim ee As Double  '第一偏心率e的平方
        Dim e2 As Double   '第二偏心率
        ee = 2 * zf - zf * zf
        e2 = Math.Sqrt(ee / (1 - ee))
        m0 = za * (1 - ee)
        m2 = 3 * ee * m0 / 2
        m4 = 5 * ee * m2 / 4
        m6 = 7 * ee * m4 / 6
        m8 = 9 * ee * m6 / 8
        a0 = m0 + m2 / 2 + 3 / 8 * m4 + 5 / 16 * m6 + 35 / 128 * m8
        a2 = m2 / 2 + m4 / 2 + 15 / 32 * m6 + 7 / 16 * m8
        a4 = m4 / 8 + 3 / 16 * m6 + 7 / 32 * m8
        a6 = m6 / 32 + m8 / 16
        a8 = m8 / 128
        Dim Bf As Double
        Dim Bf_1 As Double
        Dim Bf_nxt As Double
        Dim Bf_cur As Double
        Dim F_Bf As Double
        Dim Bf_delta As Double
        Bf_1 = zx / a0
        Bf_nxt = Bf_1
        Do
            Bf_cur = Bf_nxt
            F_Bf = -a2 / 2 * Math.Sin(Bf_cur * 2) + a4 / 4 * Math.Sin(Bf_cur * 4) - a6 / 6 * Math.Sin(Bf_cur * 6) + a8 / 8 * Math.Sin(Bf_cur * 8)
            Bf_nxt = (zx - F_Bf) / a0
            Bf_delta = Math.Abs(Bf_nxt - Bf_cur)
        Loop Until Bf_delta < 1 * 10 ^ (-8)
        Bf = Bf_nxt
        Dim yf_pow2 As Double  'ηf的平方
        Dim tf As Double
        Dim tf_pow2 As Double  'tf的平方
        Dim Nf As Double  '卯酉圈曲率半径
        Dim Mf As Double  '子午圈曲率半径
        yf_pow2 = e2 * e2 * Math.Cos(Bf) * Math.Cos(Bf)
        tf = Math.Tan(Bf)
        tf_pow2 = tf * tf
        Nf = za / Math.Sqrt(1 - ee * Math.Sin(Bf) * Math.Sin(Bf))
        Mf = za * (1 - ee) / (Math.Sqrt(1 - ee * Math.Sin(Bf) * Math.Sin(Bf))) ^ 3
        Dim Z As Double
        Z = zy / Nf
        Dim B_hudu As Double
        Dim l_hudu As Double
        B_hudu = Bf - tf / 2 / Mf * zy * Z * (1 - (5 + 3 * tf_pow2 + yf_pow2 - 9 * yf_pow2 * tf_pow2) / 12 * Z * Z + (61 + 90 * tf_pow2 + 45 * tf_pow2 * tf_pow2) / 360 * Z ^ 4)
        l_hudu = Z / Math.Cos(Bf) * (1 - (1 + 2 * tf_pow2 + yf_pow2) / 6 * Z * Z + 1 / 120 * (5 + 28 * tf_pow2 + 24 * tf_pow2 ^ 2 + 6 * yf_pow2 + 8 * yf_pow2 * tf_pow2) * Z ^ 4)
        zB = B_hudu * 180 / PI
        zL = zL0 + l_hudu * 180 / PI
    End Function