大地测量中，高斯投影正反算的问题关键是子午线弧长数学模型的解算和求解底点纬度B_f^[1]。为了提高两者的计算精度、优化计算方法、简化编程或提高运算效率，相关学者提出许多计算方法。文献[2-7]分别引入了递归算法、递推公式、数值积分以及牛顿迭代、插值公式直接求解等方法，使得子午线弧长解算精度优于0.1 mm，底点纬度解算精度达到10^-5(″)量级；文献[8, 9]提出基于第二类椭圆和依不同纬度变量的子午线弧长计算方法，丰富了子午线弧长计算的理论；文献[10]对子午线弧长公式的泰勒级数展开形式进行了简化和优化；文献[11-13]将计算子午线弧长与底点纬度视作求解一阶常微分方程，通过数值计算的方法进行求解，并对各种计算方法的解算精度、解算速度进行了分析与评价；文献[14]利用Excel的“单变量求解”功能实现底点纬度的免编程计算。

本文在综合考虑编程简洁性、运算速度和计算精度的基础上，根据文献[4]的研究思路设计了基于Romberg积分和牛顿迭代的高斯投影坐标正反算算法，并通过C#编程实现了高斯投影坐标正反算。

1 基于Romberg积分和牛顿迭代的高斯投影坐标正反算算法

由已知大地坐标(B，L)计算高斯投影平面坐标(x, y):

$ \begin{array}{l} x = X + \frac{N}{2}\sin B\cos B{l^2} + \frac{N}{{24}}\sin B{\cos ^3}B\left( {5 - {t^2} + 9{\eta ^2} + 4{\eta ^4}} \right){l^4} + \\ \;\;\;\frac{N}{{720}}\sin B{\cos ^5}B\left( {61 - 58{t^2} + {t^4}} \right){l^6} \end{array} $

(1)

$ \begin{array}{l} y = N\cos Bl + \frac{N}{6}{\cos ^3}B\left( {1 - {t^2} + {\eta ^2}} \right){l^3} + \\ \;\frac{N}{{120}}{\cos ^5}B\left( {5 - 18{t^2} + {t^4} + 14{\eta ^2} - 58{\eta ^2}{t^2}} \right){l^5} \end{array} $

(2)

式中，X为自赤道量超至已知点的子午线弧长；N为卯酉圈半径，$n = \frac{{a{\rm{^\circ }}}}{W}$，$W = \sqrt {1 - {e^2}{{\sin }^2}B} $, a°为椭球长半径，e为椭球第一偏心率；t=tan B；η²=e′²cos²B, e′为椭球第二偏心率；l为投影带内中央子午线经差(单位为弧度)。当l < 3.5°时，上式换算精度可精确至0.001 mm。

高斯投影坐标反算公式为：

$ \begin{array}{l} B = {B_f} - \frac{{{t_f}}}{{2{M_f}{N_f}}}{y^2} + \frac{{{t_f}}}{{24{M_f}N_f^3}}\left( {5 + 3t_f^2 + \eta _f^2 - 9\eta _f^2t_f^2} \right){y^4}\\ - \frac{{{t_f}}}{{720{M_f}N_f^5}}\left( {61 + 90t_f^2 + 45t_f^4} \right){y^6} \end{array} $

(3)

$ \begin{array}{l} l = \frac{1}{{{N_f}\cos {B_f}}}y - \frac{1}{{6N_f^3\cos {B_f}}}\left( {1 + 2t_f^2 + \eta _f^2} \right){y^3} + \\ \frac{1}{{120N_f^5\cos {B_f}}}\left( {5 + 28t_f^2 + 24t_f^4 + 6\eta _f^2 + 8\eta _f^2t_f^2} \right){y^5} \end{array} $

(4)

式中，B_f为底点纬度，其余有下标的各量都是B_f的函数。当l < 3.5°时，上式换算精确达0.000 1″；B和l的单位为弧度。

在高斯投影坐标正反算中最为关键的是子午线弧长与底点纬度的计算。

1.1 子午线弧长的计算

子午线弧长数学模型是计算从赤道到任意纬度B的平行圈之间的弧长：

$ X = \int_0^B M {\rm{d}}B = a\left( {1 - {e^2}} \right)\int_0^B {{{\left( {1 - {e^2}{{\sin }^2}B} \right)}^{ - \frac{3}{2}}}} {\rm{d}}B $

(5)

式中，M为子午线曲率半径。

求解式(5)的传统算法^[1]是把被积函数按级数展开，为便于积分将正弦函数的幂函数展开为余弦的倍数函数。积分后经整理为：

$ X = {a_0}B - \frac{{{a_2}}}{2}\sin 2B + \frac{{{a_4}}}{4}\sin 4B - \frac{{{a_6}}}{6}\sin 6B + \frac{{{a_8}}}{8}\sin 8B $

(6)

式中,

$ \left\{ {\begin{array}{*{20}{l}} {{a_0} = {m_0} + \frac{{{m_2}}}{2} + \frac{3}{8}{m_4} + \frac{5}{{16}}{m_6} + \frac{{35}}{{128}}{m_8}}\\ {{a_2} = \frac{{{m_2}}}{2} + \frac{{{m_4}}}{2} + \frac{{15}}{{32}}{m_6} + \frac{7}{{16}}{m_8}}\\ {{a_4} = \frac{{{m_4}}}{8} + \frac{3}{{16}}{m_6} + \frac{7}{{32}}{m_8}}\\ {{a_6} = \frac{{{m_6}}}{{32}} + \frac{{{m_8}}}{{16}}}\\ {{a_8} = \frac{{{m_8}}}{{128}}} \end{array}} \right. $

其中，m₀=a(1-e²), m₂=${\frac{3}{{2}}}$e²m₀, m₄=${\frac{5}{{4}}}$e²m₂, m₆=${\frac{7}{{6}}}$e²m₄, m₈=${\frac{9}{{8}}}$e²m₆。

式(6)中，$\frac{{{a_8}}}{8} = \frac{{{m_8}}}{{1024}}$，小于0.000 03 m，故可忽略不计。

为求解式(5)，本文采用Romberg积分方法。Romberg积分将积分区间逐次分半，得到低次复核梯形积分序列，然后运用Richardson外推算法加速收敛，并将得到的相邻两次结果进行比较，直到较差值小于要求的精度时为止^{[15, 16]}。

1.2 底点纬度求解

为求解B_f，用式(5)或式(6)构造关于B_f的一元方程：

$ f\left( {{B_f}} \right) = a\left( {1 - {e^2}} \right)\int_0^{{B_f}} {{{\left( {1 - {e^2}{{\sin }^2}{B_f}} \right)}^{ - \frac{3}{2}}}} {\rm{d}}{B_f} - X $

(7)

$ \begin{array}{*{20}{r}} {f\left( {{B_f}} \right) = {a_0}{B_f} - \frac{{{a_2}}}{2}\sin 2{B_f} + \frac{{{a_4}}}{4}\sin 4{B_f} - }\\ {\frac{{{a_6}}}{6}\sin 6{B_f} + \frac{{{a_8}}}{8}\sin 8{B_f} - X} \end{array} $

(8)

问题转化为求式(7)或式(8)方程的解，利用牛顿迭代算法进行求解^{[6, 15]}。设已知方程f(x)=0有近似根x_k(假定f′(x_k)≠0)，将函数f(x)在x_k展开，有f(x)≈f(x_k)+f′(x_k)(x-x_k)，于是方程f(x)=0可近似地表示为f(x_k)+f′(x_k)(x-x_k)=0，记其根为x_k+1，则x_k+1的计算式为：

$ {x_{k + 1}} = {x_k} - \frac{{f\left( {{x_k}} \right)}}{{{f^\prime }\left( {{x_k}} \right)}}, k = 0, 1, 2, \cdots $

(9)

综上，求解式(7)或式(8)式一元方程解的牛顿迭代式为：${B_{{f_{k + 1}}}} = {B_{{f_k}}} - \frac{{f\left( {{B_{{f_k}}}} \right)}}{{{f^\prime }\left( {{B_{{f_k}}}} \right)}}$，选取大地纬度B的迭代初始值${B_0} = \frac{X}{a}$，迭代至$\left| {{B_{{f_{k + 1}}}} + {B_{{f_k}}}} \right| < \varepsilon $为止。

2 算法实现与实例精度验证

本文所述高斯投影坐标正反算算法的实现过程如下。利用Romberg积分算法求解式(5)，得到子午线弧长X，然后将其代入式(1)、式(2)做高斯投影坐标正算；利用牛顿迭代算法求解式(7)(通过调用上述Romberg积分算法模块计算子午线弧长X)，求得底点纬度并将其代入式(3)、式(4)做高斯投影坐标反算。

利用式(6)求解子午线弧长X，然后将其代入式(1)、式(2)做高斯投影坐标正算；同理，利用式(8)求解底点纬度，并将其代入式(3)、式(4)做高斯投影坐标反算。此为高斯投影坐标正反算的传统算法过程。

为了验证本文所提算法的有效性和正确性，利用C#编程实现算法。并以CGCS2000椭球(a=63 781 37 m, α=1/298.257 222 101)中央子午线为114°的一组模拟数据为例, 将其计算结果与传统算法进行比较，见表 1、表 2。大地坐标的输入及输出形式作如下说明：3°40′40.123 45″其形式为3.404 012 345。

由表 1、表 2可知，高斯投影坐标正算达到1 mm以内的计算精度，高斯投影坐标反算的计算精度达到1×10^-4(″)以内, 说明该算法精度与传统算法相当。

3 结束语

Romberg积分在解算子午线弧长时涉及到的参数较少，省去多项式推导及大量系数的计算，因此该算法代码书写相对简洁，计算精度可扩展性比较强，不考虑多项式余项的截断误差。对于底点纬度求解，本文算法与插值公式直接求解法^[7]相比，省去了复杂数学公式推导和精度估计，更适合编程采用，同时该算法具有收敛速度快和计算精度高的特点，使得大地坐标反算结果与初始值符合性很好。