十多年前我在做 GIS 坐标转换程序时曾参考过网上流传的一套高斯-克吕格投影正反算 C/C 代码。当时自己也结合《大地测量学基础》理解过公式并使用实际控制点进行过验证。因为计算结果一直满足工程使用精度所以这套代码[CSDN也有引用这段代码]后来也被我用在多个程序中。最近重新翻老代码逐项对照教材和公式检查才发现其中竟然隐藏着几个非常不容易发现的错误。[这里有论文]这些错误有一个共同特点都出现在高阶小量项中所以程序能正常运行结果也“看起来完全正确”普通工程精度下甚至很难发现。网上至今仍能搜索到大量使用相同代码的版本因此把这次检查结果记录下来给后来使用这段代码的朋友一个提醒。一、原高斯正算代码中的第一个错误网上流传版本中有a3(double)15/64*e2*e2 (double)105/256*e2*e2*e2 (double)2205/4096*e2*e2*e2*e2 (double)10359/16384*e2*e2*e2*e2*e2;其中10359/16384应为10395/16384即a3(double)15/64*e2*e2 (double)105/256*e2*e2*e2 (double)2205/4096*e2*e2*e2*e2 (double)10395/16384*e2*e2*e2*e2*e2;这是一个非常像人工抄写时把数字顺序写反的错误。二、第二个错误131072 少了一个“1”原代码a4(double)35/512*e2*e2*e2 (double)315/2048*e2*e2*e2*e2 (double)31185/13072*e2*e2*e2*e2*e2;其中31185/13072正确应为31185/131072修正后a4(double)35/512*e2*e2*e2 (double)315/2048*e2*e2*e2*e2 (double)31185/131072*e2*e2*e2*e2*e2;这个错误更加隐蔽因为只是在一个很长的分母中少了一个数字1。而且这一项属于很高阶的小量所以实际坐标结果受到的影响非常小。这也是为什么这类错误能够在网上流传很多年而没有被轻易发现。三、第三个错误高斯正算 y 坐标中的符号错误原代码yN*m (double)1/6*(1-t*tq2)*N*m*m*m (double)1/120* (5-18*t*tt*t*t*t-14*q2-58*q2*t*t) *N*m*m*m*m*m;这里-14*q2应为14*q2正确公式为yN*m (double)1/6*(1-t*tq2)*N*m*m*m (double)1/120* (5-18*t*tt*t*t*t14*q2-58*q2*t*t) *N*m*m*m*m*m;也就是高斯正算横坐标五次项\[ 5-18t^2t^414\eta^2-58\eta^2t^2 \]这里的 \(q2\) 实际就是\[ \eta^2e^2\cos^2B \]这个符号错误同样处于高阶项因此在正常 3°带、6°带使用范围内对结果影响非常小。但公式本身还是应该改正确。四、高斯反算中还有一个 256 / 252 的错误网上流传版本中常见latitude1 fai - (NN * tan(fai) / R) * ( D * D / 2 - (5 3 * T 10 * C - 4 * C * C - 9 * ee) * D * D * D * D / 24 (61 90 * T 298 * C 45 * T * T - 256 * ee - 3 * C * C) * D * D * D * D * D * D / 720 );其中-256 * ee正确应为-252 * ee所以正确代码为latitude1 fai - (NN * tan(fai) / R) * ( D * D / 2 - (5 3 * T 10 * C - 4 * C * C - 9 * ee) * D * D * D * D / 24 (61 90 * T 298 * C 45 * T * T - 252 * ee - 3 * C * C) * D * D * D * D * D * D / 720 );注意最后仍然是- 3 * C * C不是加号。五、为什么这些错误一直没有被发现原因其实很简单这些错误几乎全部位于高阶展开项。例如椭球第一偏心率平方\[ e^2\approx0.0067 \]到了\[ (e^2)^5 \]已经是一个非常小的量。因此即使系数有一点错误最终对坐标的影响也可能只有亚毫米几十微米甚至更小。我自己当年的程序里还留着一句注释{//经验证可用.现在回头看这句话其实也没有错。因为当时验证的是在实际工程范围内计算结果是否满足使用精度。而不是每一个高阶系数是否都逐项抄写完全正确。这两者是不同的概念。所以“经验证可用”并不等于“代码里的每一个数学系数绝对没有错误”。这也是这次重新检查老程序给我的一个提醒。六、还有两个不是公式本身、但值得修改的问题1.atan(y/x)应改成atan2(y,x)很多老代码中笛卡尔坐标转大地坐标时写的是pcg-longitude atan(pcc-y / pcc-x);更稳妥、正确的写法应该是pcg-longitude atan2(pcc-y, pcc-x);原因是atan(y/x)无法正确区分象限同时在x0时还存在除零问题。而atan2(y,x)能够直接根据 x、y 的符号判断正确象限。一些针对中国区域写的老代码后面会通过180 longitude人为修正象限。这种写法在特定区域可能一直“能用”但并不是一个通用正确的处理方法。2. 一个非常典型的 C 小 bug原代码pDestination-y (iprjno*1000000L);这句实际上什么也没有改变。应该是pDestination-y iprjno * 1000000L;前者只是计算了一个表达式然后把结果丢掉。如果程序平时没有启用带号前缀这个错误同样可能多年都不会暴露。七、修正后的高斯正算核心代码我目前整理后的核心部分如下e2 2*f - f*f; l L - L0; t tan(B); m l * cos(B); N a / sqrt(1 - e2*sin(B)*sin(B)); q2 e2/(1-e2) * cos(B)*cos(B); a1 1 (double)3/4*e2 (double)45/64*e2*e2 (double)175/256*e2*e2*e2 (double)11025/16384*e2*e2*e2*e2 (double)43659/65536*e2*e2*e2*e2*e2; a2 (double)3/4*e2 (double)15/16*e2*e2 (double)525/512*e2*e2*e2 (double)2205/2048*e2*e2*e2*e2 (double)72765/65536*e2*e2*e2*e2*e2; a3 (double)15/64*e2*e2 (double)105/256*e2*e2*e2 (double)2205/4096*e2*e2*e2*e2 (double)10395/16384*e2*e2*e2*e2*e2; a4 (double)35/512*e2*e2*e2 (double)315/2048*e2*e2*e2*e2 (double)31185/131072*e2*e2*e2*e2*e2; b1 a1*a*(1-e2); b2 (double)-1/2*a2*a*(1-e2); b3 (double)1/4*a3*a*(1-e2); b4 (double)-1/6*a4*a*(1-e2); c0 b1; c1 2*b2 4*b3 6*b4; c2 -(8*b3 32*b4); c3 32*b4; s c0*B cos(B) * ( c1*sin(B) c2*sin(B)*sin(B)*sin(B) c3*sin(B)*sin(B)*sin(B)*sin(B)*sin(B) ); x s (double)1/2*N*t*m*m (double)1/24* (5-t*t9*q24*q2*q2) *N*t*m*m*m*m (double)1/720* (61-58*t*tt*t*t*t) *N*t*m*m*m*m*m*m; y N*m (double)1/6* (1-t*tq2) *N*m*m*m (double)1/120* (5-18*t*tt*t*t*t14*q2-58*q2*t*t) *N*m*m*m*m*m; y 500000;八、最后说几句这段代码至少在国内 GIS、GPS、测绘程序圈里流传了很多年。我自己在大约 2007 年做坐标转换程序时就曾经参考过其中一部分代码。当时结合武汉大学出版社《大地测量学基础》理解了高斯投影、椭球、大地坐标与空间直角坐标等关系并进行了实际数据验证。后来因为“经验证可用”很多年也没有再去逐项检查这些高阶系数。直到最近重新翻老代码才把这些小问题逐个找出来。所以把结果留下来。如果你正好在网上搜到类似下面这些代码10359/1638431185/13072-14*q2-256*ee建议检查一下。正确应分别为10395/1638431185/13107214*q2-252*ee这些错误通常不会让你的程序“算错很多”但既然已经知道了就没有理由继续让它们流传下去。也算给这段在互联网上流传了很多年的老代码做一次小小的勘误。—— DZQ2026年10月7日