GPS定位系列写到第五篇前几篇把坐标系、时间系统这些地基铺完了今天聊一个真正动手的活高斯投影正反算的C程序。干过测绘或者GIS开发的人都知道GPS接收机直接输出的通常是WGS84或CGCS2000的大地坐标也就是B纬度、L经度、H大地高。可是放到工程里无论是施工放样、地形测图还是无人机影像处理最后基本都得落在平面坐标上这个从经纬度到平面坐标的换算绝大部分场景就是高斯投影的活。这篇内容适合两类人一类是刚接触坐标转换、想把原理和代码一起搞明白的初学者另一类是手头正需要写坐标转换模块、想找一份可以放心复用代码的开发者。我会把正算、反算的公式怎么来的、程序怎么搭、实测会踩哪些坑一件件说清楚。1. GPS定位链路里高斯投影卡在哪一环为什么躲不掉1.1 接收机给的是“球面地址”施工要的是“平面坐标”GPS定位最终给到用户手里的通常是一组大地坐标其中纬度B和经度L本质上是球面上的角度量。你可以把经纬度想象成地球表面上的“门牌号”它能把每个点唯一确定下来但你要是拿着两个点的经纬度去量距离、算面积、放样就很痛苦了。地面上两点之间的距离在椭球面上是一条测地线直接算测地线长度需要解比较复杂的微分几何问题而在平面坐标系里两点之间就是勾股定理几行代码就搞定。CAD制图、RTK放样、地形图分幅、土方计算这些全部建立在平面坐标之上。所以GPS定位数据进工程系统之前必须先把大地坐标变成平面坐标。那为什么偏偏要用投影而不是直接找一个“平面坐标系”让GPS输出因为地球是一个不规则的椭球体不可能把椭球表面完美地展开成平面而不产生任何变形。你能做到的就是选择一种投影方式让某种变形尽量小或者让某种特性保持住。测绘领域在我国用得最多的就是高斯投影它把椭球面上的点通过等角投影的方式映射到平面上。1.2 等角投影这个性质决定了它的适用范围高斯投影的全称是“等角横切椭圆柱投影”也叫高斯-克吕格投影。核心性质是投影前后微分图形保持相似角度不发生变形。这句话翻译成人话就是地面上一个小范围的三角形或圆形投影到平面上之后图形的角度关系不变只是整体放缩了一个比例。对测图来说角度不变意味着图上量出来的方位角可以直接用于实地放样这是工程测量最看重的一件事。代价也很明显长度会有变形而且离中央子午线越远长度变形越大。中央子午线本身是长度无变形的线离开它越远变形迅速增加。6度带边缘最大长度变形约0.14%3度带边缘约0.034%。这也是为什么高斯投影一定要分带——把变形控制在允许范围内。理解了这一点你就知道为什么高斯投影程序里那个“中央子午线L0”参数举足轻重它不只是公式里的一个减数而是整个投影变形的基准。1.3 分带方式与中央子午线选错带等于白干我国常用高斯投影分带是6度带和3度带。6度带从0°经线起每隔6°划一带带号N等于经度L除以6取整再加1中央子午线L06N-33度带从1.5°经线起每隔3°划一带带号n等于(L1.5)除以3取整中央子午线L03n。但GPS作业和传统国家基本比例尺测图还不完全一样。很多时候你手里拿到的项目是一个城市、一个工地范围不大完全可以按测区中心经度自定义中央子午线这就是所谓的“任意带”。只要确保测区离中央子午线不太远长度变形就能控制在毫米级甚至更低。因此实际写程序的时候中央子午线最好做成接口参数而不是写死在代码里。另一个高频概念是y坐标的加常数为了避免中央子午线以西出现负坐标我国规定在通用坐标中y值加500km并在前面冠以带号。高斯正反算公式里用的y必须是自然值也就是去掉带号、去掉500km后的数值。这个坑后面我还会专门说。2. 正算公式的“零件”拆解搞懂每一项在干嘛再写代码2.1 输入、输出与那个关键的经差l高斯投影正算的输入是大地坐标B纬度、L经度以及中央子午线经度L0输出是平面坐标x北方向也叫纵坐标、y东方向横坐标。B、L、L0都是弧度制三者之间有一个关键中间量l L - L0l叫经差表示待投影点离开中央子午线的角度距离。整个高斯投影级数公式都是在中央子午线附近展开的所以公式里的每一项都是l的幂次项。l越小收敛越快l越大高阶项越不能省。有些初学者把L直接丢进公式里算出来自然是一堆乱码还有人把L的十进制度数当成弧度用结果差了57.3倍——这类问题后面统一在“高频错误”里讲。2.2 辅助符号表五个量决定整个公式的命运高斯投影公式里有几个基础量必须先算出来它们都和椭球参数以及纬度B有关符号含义计算公式e²第一偏心率的平方e² (a²-b²)/a² 2f-f²e²第二偏心率的平方e² e²/(1-e²)N卯酉圈曲率半径N a/sqrt(1-e²sin²B)M子午圈曲率半径M a(1-e²)/(1-e²sin²B)^(3/2)t纬度正切值t tanBη²辅助量η² e²cos²BX子午线弧长从赤道到纬度B的椭球面弧长这里最容易懵的是N和M为什么不一样。你可以这样理解椭球面上过某一点可以切出无数个法截面不同方向的法截线曲率半径不同。沿子午圈方向的那个曲率半径叫M沿卯酉圈方向与子午圈垂直的法截面的叫N。在赤道上N就是赤道半径a在极点上N和M都趋于同一个值。因为椭球是“扁”的导致这两个半径在不同纬度上表现出不同的变化规律。高斯投影公式里x方向主要涉及M和N的组合y方向涉及cosB和N所以不把这些基础量搞对后面全是白搭。2.3 子午线弧长X正算里x坐标的“起跑线”子午线弧长X表示从赤道沿子午圈量到纬度B的椭球面距离。它只由B决定和经差l无关。展开式常见写成X a(1-e²)[A0·B A2·sin2B A4·sin4B A6·sin6B A8·sin8B]其中系数由e²、e⁴、e⁶、e⁸组合而成A0 1 (3/4)e² (45/64)e⁴ (175/256)e⁶ (11025/16384)e⁸A2 -(3/4)e² - (15/16)e⁴ - (525/512)e⁶ - (2205/2048)e⁸A4 (15/64)e⁴ (105/256)e⁶ (2205/4096)e⁸A6 -(35/512)e⁶ - (315/2048)e⁸A8 (315/16384)e⁸不同教材对系数定义可能略有差异有的把A2写成正号再在X式子里放减号本质是一样的。只要你保证正算和反算用的是同一套X程序内部就是自洽的。但要注意不同椭球参数算出来的X系数不同所以千万不要网上抄一段“经验系数”就固定死一定要让系数从椭球参数动态算出来。高斯正算的x坐标本质上就是“子午线弧长X加上经差l引起的改正项”。这也是为什么中央子午线上的点l0的x坐标正好等于Xy坐标等于0——很好理解投影点落在中央子午线上平面上本来就不该有东西方向的偏移。3. 高斯投影正算的C实现代码骨架和每一步对应哪个公式3.1 结构体设计把椭球参数当“一等公民”写C程序时我习惯把椭球参数定义成独立结构体而不是散落成全局变量。原因很简单同一个工程里很可能既用到CGCS2000又要转换西安80或WGS84坐标椭球参数集中管理切换起来只改一个初始化赋值。坐标结构体则分成“大地坐标”和“投影坐标”两类避免传参时把x、y和B、L搞混。#include stdio.h #include math.h #define PI 3.14159265358979323846 typedef struct { double a; // 长半轴单位米 double f; // 扁率 } Ellipsoid; typedef struct { double B; // 纬度弧度 double L; // 经度弧度 } GeoCoord; typedef struct { double x; // 北坐标米 double y; // 东坐标自然值米 } ProjCoord;这里的y我用“自然值”也就是没有加500km、没有带带号的数值。如果你手里的坐标是“34065789.123”这种带带号的通用值进函数之前必须先把带号去掉、再减掉500km否则反算结果直接飞掉。后面会专门讲这个操作。3.2 正算函数完整代码void GaussForward(const Ellipsoid *ellip, double L0, const GeoCoord *geo, ProjCoord *proj) { double a ellip-a; double f ellip-f; double e2 2.0*f - f*f; double ep2 e2 / (1.0 - e2); double B geo-B; double l geo-L - L0; double sinB sin(B); double cosB cos(B); double cosB2 cosB * cosB; double cosB3 cosB2 * cosB; double cosB5 cosB3 * cosB2; double t tan(B); double t2 t * t; double t4 t2 * t2; double eta2 ep2 * cosB2; double N a / sqrt(1.0 - e2 * sinB * sinB); // 子午线弧长X系数随椭球参数动态计算 double e4 e2 * e2; double e6 e4 * e2; double e8 e4 * e4; double A0 1.0 3.0/4.0*e2 45.0/64.0*e4 175.0/256.0*e6 11025.0/16384.0*e8; double A2 -3.0/4.0*e2 - 15.0/16.0*e4 - 525.0/512.0*e6 - 2205.0/2048.0*e8; double A4 15.0/64.0*e4 105.0/256.0*e6 2205.0/4096.0*e8; double A6 -35.0/512.0*e6 - 315.0/2048.0*e8; double A8 315.0/16384.0*e8; double X a * (1.0 - e2) * (A0*B A2*sin(2.0*B) A4*sin(4.0*B) A6*sin(6.0*B) A8*sin(8.0*B)); double l2 l * l; double l3 l2 * l; double l4 l3 * l; double l5 l4 * l; double l6 l5 * l; proj-x X N * sinB * cosB * l2 / 2.0 N * sinB * cosB3 * (5.0 - t2 9.0*eta2 4.0*eta2*eta2) * l4 / 24.0 N * sinB * cosB5 * (61.0 - 58.0*t2 t4 270.0*eta2 - 330.0*eta2*t2) * l6 / 720.0; proj-y N * cosB * l N * cosB3 * (1.0 - t2 eta2) * l3 / 6.0 N * cosB5 * (5.0 - 18.0*t2 t4 14.0*eta2 - 58.0*eta2*t2) * l5 / 120.0; }3.3 为什么这些级数项不能随意砍很多初学者看到x、y那一大串公式就想精简觉得l是弧度制的小数高次项肯定可以忽略。这个想法在经差很小的时候没错但“很小”是有边界的。以3度带边缘l1.5°为例l≈0.02618radx公式里的l²项算下来接近1000米量级无论如何不能丢l⁴项大约0.1米量级做厘米级精度必须保留l⁶项在3度带内不足毫米但如果你用6度带或任意带让经差达到3°l⁶项就会冲到毫米级。GPS后处理的成果精度通常要求到毫米级所以代码里保留到l⁶是合理的。y方向的l³项对应几十米量级的修正量l⁵项在高纬度也有毫米级贡献同样不能省。这段代码里的每一项都能和上一节的公式对上号x第一项是子午弧长X第二项是l²项的“南北方位总体平移”第三、第四项是更高阶的“东西方向不均匀拉伸带来的南北修正”y第一项是中央子午线附近的“近似水平距离”后面几项则考虑椭球曲率和经差变大后的细节修正。把这些都印在脑子里调试时看量级就能判断某个数字是不是算错了。4. 反算的实现底点纬度迭代和公式展开4.1 反算的核心难点先求底点纬度Bf反算比正算难在“倒着来”。已知平面坐标x、y求大地坐标B、L。y可以直接反推经差l但l的表达式里带N、t这些和纬度有关的量所以得先知道一个中间纬度。这个中间纬度叫底点纬度Bf定义是子午线弧长恰好等于x的那个纬度。也就是说先解方程X(Bf)x。问题在于X(B)是超越函数没有解析反函数只能迭代。常用方法是牛顿迭代因为初值给得好收敛极快。初值可以用Bfx/[a(1-e²)A0]这个初值相当于把子午弧长近似成线性关系。然后每步计算当前X(Bf)和x的差值除以该纬度处的子午圈曲率半径Mf修正量就是差值除以导数。一般四五次迭代就收敛到1e-12弧度对应的弧长误差是微米级。void GaussInverse(const Ellipsoid *ellip, double L0, const ProjCoord *proj, GeoCoord *geo) { double a ellip-a; double f ellip-f; double e2 2.0*f - f*f; double ep2 e2 / (1.0 - e2); double x proj-x; double y proj-y; double e4 e2 * e2; double e6 e4 * e2; double e8 e4 * e4; double A0 1.0 3.0/4.0*e2 45.0/64.0*e4 175.0/256.0*e6 11025.0/16384.0*e8; double A2 -3.0/4.0*e2 - 15.0/16.0*e4 - 525.0/512.0*e6 - 2205.0/2048.0*e8; double A4 15.0/64.0*e4 105.0/256.0*e6 2205.0/4096.0*e8; double A6 -35.0/512.0*e6 - 315.0/2048.0*e8; double A8 315.0/16384.0*e8; // 牛顿迭代求底点纬度 Bf double Bf x / (a * (1.0 - e2) * A0); double d; do { double Xf a * (1.0 - e2) * (A0*Bf A2*sin(2.0*Bf) A4*sin(4.0*Bf) A6*sin(6.0*Bf) A8*sin(8.0*Bf)); double Mf a * (1.0 - e2) / pow(1.0 - e2*sin(Bf)*sin(Bf), 1.5); d (Xf - x) / Mf; Bf - d; } while (fabs(d) 1e-12); double sinBf sin(Bf); double cosBf cos(Bf); double cosBf2 cosBf * cosBf; double tf tan(Bf); double tf2 tf * tf; double tf4 tf2 * tf2; double etaf2 ep2 * cosBf2; double Nf a / sqrt(1.0 - e2 * sinBf * sinBf); double Mf a * (1.0 - e2) / pow(1.0 - e2 * sinBf * sinBf, 1.5); double y2 y * y; double y3 y2 * y; double y4 y3 * y; double y5 y4 * y; double y6 y5 * y; // 反算纬度B和经差l注意符号和量纲 double B Bf - tf / (2.0 * Mf * Nf) * y2 tf / (24.0 * Mf * Nf * Nf * Nf) * (5.0 3.0*tf2 etaf2 - 9.0*etaf2*tf2) * y4 - tf / (720.0 * Mf * Nf * Nf * Nf * Nf * Nf) * (61.0 90.0*tf2 45.0*tf4) * y6; double l y / (Nf * cosBf) - (1.0 2.0*tf2 etaf2) / (6.0 * Nf * Nf * Nf * cosBf) * y3 (5.0 28.0*tf2 24.0*tf4 6.0*etaf2 8.0*etaf2*tf2) / (120.0 * Nf * Nf * Nf * Nf * Nf * cosBf) * y5; geo-B B; geo-L L0 l; }4.2 反算公式里的符号陷阱与量纲检查反算公式看起来和正算“长得像”但符号完全是反着的。纬度B是在Bf的基础上先减一个y²项再加一个y⁴项再减一个y⁶项符号交替出现经差l先减去y³项再加上y⁵项。我最早写这套程序时把第二项的加号抄成了减号做往返验证时怎么都差一个微小角度查了很久才发现是这里出了问题。所以写完反算函数之后第一件事就是拿正算的输出做一次往返测试如果正反算互逆正算出的x、y反算回B、L和原始值差异应该在1e-10量级以下。另一个容易忽视的是l的计算结果需要加回L0得到真正的经度L。有人反算出来L-3度还以为是纬度输错了其实是忘了加L0。这种错误在代码里非常隐蔽因为编译器不会报警程序也“能跑”只有对比已知控制点才暴露。5. 从测试到实跑精度验证和几个“看起来对实际错”的典型案例5.1 往返验证测试程序自洽性的最低要求正反算写完先做一个简单的往返测试。选一个测试点纬度30度、经度120度中央子午线选117度这样经差有3度处于6度带边缘能检验高阶项的作用。然后调用正算得到x、y再调用反算回到经纬度打印差值。int main(void) { // CGCS2000椭球参数 Ellipsoid e {6378137.0, 1.0/298.257222101}; double L0 117.0 * PI / 180.0; // 中央子午线弧度 GeoCoord g, g2; ProjCoord p; g.B 30.0 * PI / 180.0; g.L 120.0 * PI / 180.0; GaussForward(e, L0, g, p); GaussInverse(e, L0, p, g2); printf(x%.6f y%.6f\n, p.x, p.y); printf(dL%.3e dB%.3e\n, g2.L - g.L, g2.B - g.B); return 0; }我在这台机器上跑出来的结果dL和dB都在1e-11弧度附近折算成平面距离是微米量级。这个测试能证明代码内部正反算是一致的但记住往返验证只能证明“程序自己对自己”不能证明输出坐标和外部坐标系统一致。真正要验证外部一致必须拿测区控制点成果来测正算出来的x、y与控制点平面坐标对比误差应在毫米级才算数。5.2 高频错误对照表症状、原因与排查方向症状可能原因检查方法y坐标整体偏移几十公里甚至上百公里中央子午线设错或者带号取错用带号反推中央子午线L06N-3或L03n对比当前设置x、y和已知控制点差几米到几十米椭球参数选错比如把CGCS2000的成果用西安80椭球投影核对椭球a和f检查原始成果属于哪个坐标框架y坐标带一串“34”开头的数字数值特别大通用y值没去掉带号、没减500km通用值转自然值y自然 y通用 - 整数带号×1000000 - 500000纬线看起来像抛物线经线也歪歪扭扭输入经纬度用了十进制度而公式里当成弧度弧度度×π/180进函数前统一转弧度反算经度比实际小了3度或6度整倍数反算结果l忘记加回L0检查LL0l那一步换一个测区后程序立刻失准中央子午线写死没按测区动态设置把L0作为函数参数传递我把这张表贴在工位旁边因为每次帮别人救急排查坐标问题时十个里有八个是第一行和第二行导致的。中央子午线错误通常是“看起来全都对”——高斯坐标长得都差不多六位数的x、七位数的y数字很正常但用RTK去实地放样时才能发现全错。椭球参数错误则更隐蔽因为投影公式结构完全一样只有系数有细微差别成果差异从分米级到几十米都有不对比已知点很难察觉。5.3 高程到底要不要参与投影高斯投影是在椭球面上展开的公式输入只有B和L高程H不参与。GPS直接给的是大地高工程上用的却是正常高两者之间隔着一个高程异常在平原地区一般是几厘米到几十厘米在山区可能达到几米甚至更大。所以如果你是做控制测量或需要高程成果的项目不要指望高斯投影顺带解决高程问题该做高程拟合的还是要单独做高程拟合。投影这一层只需要专心处理平面位置。不少新手把H当成0、或者干脆混进x、y里都是概念没理清。6. 我在实际项目里的一些使用体会这套程序我在不同的数据后处理工具里反复用过有几个体会很直接。第一个体会是程序里的角度量纲必须“一刀切”。我从一开始就定义所有接口都用弧度不在函数内部做度转弧度的动作数据进来的时候统一转好出去的时候统一转回。这样虽然调用方多写几行转换代码但至少不会出现“这个函数内部是度、那个函数内部是弧度”的精神分裂式设计。第二个体会是写完之后立刻用控制点实测不要用往返测试自欺欺人。往返验证能检查自洽性但坐标系是否和外部成果匹配、中央子午线设置是否符合项目要求、椭球参数是否和原始成果一致这些必须靠已知点来锚定。我经手过的一个项目里程序正反算自洽得完美但和甲方控制点一对比x方向整体偏了快一米查到最后是原始CGCS2000成果和WGS84框架之间的差异被忽略了不是高斯投影程序本身的问题。这也是个很有价值的教训坐标转换是一条链投影只是其中一个环节前后环节的框架一致性同样重要。第三个体会是关于扩展。高斯投影和UTM投影原理非常相似UTM其实就是高斯投影加上一个0.9996的比例因子再对带号和中央子午线做一些调整。如果哪天需求从“高斯投影”变成“UTM坐标”你不需要推翻这套代码只要在椭球参数、比例因子、带号处理上做少量扩展大部分框架都能复用。这也是我把椭球参数做成结构体、把L0做成参数的根本原因——工程代码最怕的就是写死。最后再说一个小技巧如果你要在程序里支持多种椭球建议把椭球参数表a和f放在配置文件或者一个集中常量表里。CGCS2000的a是6378137米扁率是1/298.257222101WGS84的a相同扁率是1/298.257223563两者椭球参数差异很小但坐标框架整体差异在中国境内可能达到厘米级甚至更大。所以宁可多写一个配置文件也不要到处散落硬编码否则后期维护的时候你会疯掉。