简介使用MATLAB实现的贝塞尔大地问题正反解计算程序基于CGCS2000国家大地坐标系椭球参数面向测绘工程与大地测量学课程学习者及科研人员。资源包含两个MATLAB函数文件和一个PDF说明文档分别实现正算由已知点大地坐标、大地线长与大地方位角求解未知点坐标及方位角反算由两点大地坐标求解大地线长与正反大地方位角并附算法流程图及使用说明。压缩包仅129KB共3个文件结构精简可直接在MATLAB中运行适合课程设计、作业参考或工程验证。目前已有1607人学习程序注释清晰便于理解贝塞尔投影逼近法的核心思想及其在椭球大地问题中的应用也可作为二次开发或算法改进的起点。1. 贝塞尔大地问题正反解的MATLAB编程为什么课程作业都会卡在这一题用 MATLAB 做椭球大地测量学课程里的贝塞尔大地问题正反解编程最容易翻车的不是公式背不熟而是球面映射和角度象限处理写错。正算已知起点纬度、方位角和长度求终点反算已知两点求距离和方位角听起来只是三角公式但一旦加上椭球扁率、归化纬度、辅助球面这三个关键词很多人的程序就开始输出南辕北辙的结果。这份资源就是把贝塞尔正反解从原理到可运行代码完整拆开包含正算函数、反算函数、测试脚本和流程图适合课程设计、测绘实习以及需要自己从零实现大地线解算的从业者。照着走一遍你能清楚知道每一步在算什么以及结果出错时该往哪里查。2. 贝塞尔法原理从椭球到辅助球面的三步映射与参数定义2.1 椭球到辅助球面归化纬度与球面方位角的关系贝塞尔法的核心思想是把椭球面上的计算搬到辅助球面上绕开椭球积分的复杂性。第一步要把大地纬度转换成归化纬度u公式是u atan((1 - f) * tan(B));这里f是椭球扁率B是大地纬度u是归化纬度。为什么要用归化纬度因为把椭球面沿法线方向投影到球面上时归化纬度能保持经度方向的一致性使得辅助球面的经度差可以直接当作椭球面的经度差这是贝塞尔法能简化的关键。第二步是求辅助球面上的方位角alpha1。椭球面上的大地方位角A1和球面方位角alpha1并不一样需要通过起点的大地纬度和归化纬度换算alpha1 asin(cos(B1) * sin(A1) / cos(u1));这个式子是从微分几何关系推导出来的。实际编程时要注意式子右边可能因为浮点误差略大于 1导致asin出现复数结果。我一般会加上一行保护ratio cos(B1) * sin(A1) / cos(u1); ratio min(max(ratio, -1), 1); alpha1 asin(ratio);参数说明B1、A1都要转成弧度u1是由B1计算得到的归化纬度alpha1是辅助球面上的大圆方位角后续所有球面三角公式都基于它展开。2.2 球面三角模型正算已知弧长求终点反算已知两点求弧长在辅助球面上大地线变成大圆弧问题就退化成球面三角。正算时已知起点归化纬度u1、球面方位角alpha1和球心角sigma求终点的归化纬度u2和经度差dLsin_u2 sin(u1) * cos(sigma) cos(u1) * sin(sigma) * cos(alpha1); u2 asin(sin_u2); dL atan2(sin(sigma) * sin(alpha1), ... cos(u1) * cos(sigma) - sin(u1) * sin(sigma) * cos(alpha1));第一行是球面三角形中的余弦定理变形用来算纬度第二行用四象限反正切atan2算经度差这样能避免atan把角度限制在[-pi/2, pi/2]的象限问题。反算时已知两端点的归化纬度u1、u2和经度差dL球心角sigma用球面余弦定理直接求cos_sigma sin(u1) * sin(u2) cos(u1) * cos(u2) * cos(dL); sigma acos(max(min(cos_sigma, 1), -1));然后求球面方位角alpha1和alpha2alpha1 atan2(cos(u2) * sin(dL), ... cos(u1) * sin(u2) - sin(u1) * cos(u2) * cos(dL)); alpha2 atan2(cos(u1) * sin(dL), ... -sin(u1) * cos(u2) cos(u1) * sin(u2) * cos(dL));这里的alpha2是从终点往起点方向看的球面方位角后续还需要修正成椭球面上的反方位角。球面模型的好处是只用余弦定理和四象限反正切没有迭代逻辑非常清晰。但有一个前提sigma和椭球面上的实际大地线长度S不是简单线性关系这是下一节要处理的问题。2.3 长度映射为什么不能直接用球面弧长乘半径辅助球面取长半轴a作为半径如果把毫秒级测线直接按sigma S / a来算短距离可能误差不大但一旦测线超过几十公里终点纬度就会偏掉。原因在于椭球面的曲率半径随纬度变化同样的球心角对应的椭球面弧长并不是常数。常见的做法是用测线平均纬度处的平均曲率半径Rm来折算e2 f * (2 - f); Bm (B1 B2) / 2; Rm a * sqrt(1 - e2) / (1 - e2 * sin(Bm)^2); S_approx Rm * sigma;这是一次近似。更严格的做法是在正算时把sigma当作未知量用迭代或者二分法求解使S_approx逼近期望的S。资源包里的完整版本就是用二分法迭代sigma而教学简化版直接采用起点平均曲率半径折算跑通流程没问题精度要求高时建议换成迭代版。参数说明e2是第一偏心率的平方由扁率f得到sin(Bm)要确保Bm是弧度Rm的单位与a一致所以S_approx的单位是米。这一段是整个贝塞尔法里最容易出问题的地方很多教程只给球面公式不提长度修正导致初学者按球面算完却对不上真实距离。3. 正算实现MATLAB函数设计、迭代求解与数值验证3.1 正算函数接口参数顺序、单位约定与格式说明设计函数的第一件事是定接口推荐把所有角度统一成弧度长度统一成米避免函数内部来回转换。我用的是下面这种签名function [B2, L2, A2] geodesic_direct(B1, L1, A1, S, a, f)输入参数B1、L1是起点大地纬度和经度弧度A1是起点大地方位角弧度从北方向顺时针量S是大地线长度米a是椭球长半轴米f是扁率无量纲。输出为终点纬度、终点经度和终点反方位角都是弧度。之所以把a和f都作为参数而不是函数内部写死是为了方便切换坐标系。做课程作业时习惯用克拉索夫斯基椭球参数是a 6378245f 1/298.3做 GPS 相关用 WGS84参数是a 6378137f 1/298.257223563。接口里带有这两个参数后换椭球只需要改调用处不需要改函数体。单位约定宁可啰嗦也要写在注释里。我见过不少同学把角度从度直接传给函数结果sin(30)被当成弧度算终点跑到南极去了。所以在函数第一行注释里明确写“所有角度均为弧度长度单位与 a 保持一致”会省很多调试时间。3.2 核心迭代从初值到二分法求球心角简化版直接用sigma S / R0一次算到底但完整版推荐用二分法迭代。原因是sigma和S的关系不是线性的尤其当测线穿过不同纬度带时用起点曲率半径折算会累计误差。二分法实现简单且稳定function [B2, L2, A2] geodesic_direct_iter(B1, L1, A1, S, a, f) e2 f * (2 - f); u1 atan((1 - f) * tan(B1)); ratio cos(B1) * sin(A1) / cos(u1); ratio min(max(ratio, -1), 1); alpha1 asin(ratio); lo 0; hi pi; for k 1:60 sigma (lo hi) / 2; % 球面三角解算中点 sin_u2 sin(u1) * cos(sigma) cos(u1) * sin(sigma) * cos(alpha1); u2 asin(sin_u2); dL atan2(sin(sigma) * sin(alpha1), ... cos(u1) * cos(sigma) - sin(u1) * sin(sigma) * cos(alpha1)); B2 atan(tan(u2) / (1 - f)); L2 L1 dL; % 用实测平均曲率半径折算当前 sigma 对应的长度 Bm (B1 B2) / 2; Rm a * sqrt(1 - e2) / (1 - e2 * sin(Bm)^2); S_calc Rm * sigma; if S_calc S hi sigma; else lo sigma; end end sigma (lo hi) / 2; % 最后一遍重新算终点 sin_u2 sin(u1) * cos(sigma) cos(u1) * sin(sigma) * cos(alpha1); u2 asin(sin_u2); dL atan2(sin(sigma) * sin(alpha1), ... cos(u1) * cos(sigma) - sin(u1) * sin(sigma) * cos(alpha1)); B2 atan(tan(u2) / (1 - f)); L2 L1 dL; % 球面反方位角再作偏量修正 alpha2 asin(sin(alpha1) * cos(u1) / cos(u2)); A2 alpha2 (A1 - alpha1); end逻辑说明二分法的上下界是0和pi因为球心角不可能超过半圈。每次用中点sigma算出临时终点再算该终点对应的椭球长度S_calc和真实S比较大小来缩小区间。循环 60 次足够把精度收敛到接近机器精度。参数说明lo、hi分别是球心角的搜索区间Bm是起终点平均纬度用于计算平均曲率半径S_calc是当前sigma对应的长度估计。这个迭代思路比牛顿法更省心不用求导也不会因为初值太差发散。有一点要提醒A2的修正方式A2 alpha2 (A1 - alpha1)是短距离近似。对于超过几百公里的测线椭球面方位角和球面方位角的差值并不是常数完整版资源里用了一个小迭代来重新求这个差值。如果只做课程演示这个近似足够。3.3 数值验证已知测线数据对照与误差检查写完函数不能直接交差至少要拿一组已知数据做闭合验证。我用过一组典型测线起点北纬 30°、东经 0°方位角 45°长度 100 kmWGS84 椭球下终点约在北纬 30.6448°、东经 0.6397°。你可以用下面这段代码调用a 6378137; f 1 / 298.257223563; B1 deg2rad(30); L1 deg2rad(0); A1 deg2rad(45); S 100e3; [B2, L2, A2] geodesic_direct_iter(B1, L1, A1, S, a, f); fprintf(终点纬度: %.6f\n, rad2deg(B2)); fprintf(终点经度: %.6f\n, rad2deg(L2)); fprintf(反方位角: %.6f\n, rad2deg(A2));如果输出和参考值相差在 1e-6 度以内说明函数基本没问题。注意这里的参考值是我用更高精度的软件工具算出来的不保证与教材附表完全一致但闭合差验证更实用把正算得到的终点和反方位角再拿去做反算反算出的距离和起始方位角应该能回到原来的值。这也是课程设计里老师最常查的一点。验证时如果发现经度差方向反了先检查atan2的第二个参数是不是cos(u1) * cos(sigma) - sin(u1) * sin(sigma) * cos(alpha1)正负号写反会直接导致经度差取补角。这种问题在入门阶段很常见后面避坑章会专门展开。4. 反算实现球面三角形求距离与方位角的完整代码4.1 反算的球面三角形解算余弦定理与四象限方位角反算的输入是两点经纬度输出是大地线长度S、起点方位角A1和终点反方位角A2。核心还是先转归化纬度再用球面余弦定理求球心角。球心角sigma由两点归化纬度u1、u2和经度差dL决定dL L2 - L1; cos_sigma sin(u1) * sin(u2) cos(u1) * cos(u2) * cos(dL); sigma acos(max(min(cos_sigma, 1), -1));求方位角时要特别注意方位角的定义从北方向顺时针量范围是[0, 2*pi]。MATLAB 的atan2输出范围是[-pi, pi]所以算完后要加一圈转换成 0~360 度否则后续画图或输出会很不方便。起点球面方位角的标准公式是alpha1 atan2(cos(u2) * sin(dL), ... cos(u1) * sin(u2) - sin(u1) * cos(u2) * cos(dL));分子是终点纬度方向的经度贡献分母是起点到终点的纬度变化项。这个形式在几乎所有测绘教材里都有但很多人把分母里的sin(u1) * sin(u2) - cos(u1) * cos(u2)抄错导致方位角差 180°。终点反方位角alpha2公式alpha2 atan2(cos(u1) * sin(dL), ... -sin(u1) * cos(u2) cos(u1) * sin(u2) * cos(dL));注意第二项第一个符号是负号因为反方位角是从终点看起点三角形顶点换了。4.2 反算代码从经纬度到距离方位角的完整实现下面给出一个可运行的反算函数采用与正算相同的简化思路function [S, A1, A2] geodesic_inverse_simple(B1, L1, B2, L2, a, f) % 贝塞尔法反算简化版 % 输入两点大地纬度、经度弧度椭球长半轴a扁率f % 输出大地线长度S米起算方位角A1弧度反方位角A2弧度 u1 atan((1 - f) * tan(B1)); u2 atan((1 - f) * tan(B2)); dL L2 - L1; % 球心角 cos_sigma sin(u1) * sin(u2) cos(u1) * cos(u2) * cos(dL); cos_sigma min(max(cos_sigma, -1), 1); sigma acos(cos_sigma); % 球面方位角 alpha1 atan2(cos(u2) * sin(dL), ... cos(u1) * sin(u2) - sin(u1) * cos(u2) * cos(dL)); alpha2 atan2(cos(u1) * sin(dL), ... -sin(u1) * cos(u2) cos(u1) * sin(u2) * cos(dL)); % 距离用平均曲率半径折算 e2 f * (2 - f); Bm (B1 B2) / 2; Rm a * sqrt(1 - e2) / (1 - e2 * sin(Bm)^2); S Rm * sigma; % 简化修正用正算的偏量近似 A1 alpha1; A2 alpha2 (alpha1 - alpha2); % 此处仅为演示完整版有独立修正 end逻辑说明前两行把大地纬度换成归化纬度后面所有计算都在球面上进行。sigma的余弦值如果因为浮点误差超出[-1,1]不加保护的话acos会返回复数这是常见问题所以先做了min/max夹逼。距离S用的是中纬度平均曲率半径这和正算的二分法标准是一致的。方位角修正方面简化版直接把球面方位角当成椭球方位角。真正严格的贝塞尔法里椭球方位角和球面方位角的差会随弧长变化完整版资源中用一个迭代修正函数来消除这个误差。课程作业如果老师不深究简化版能跑通如果要求高请使用资源包中的完整实现。4.3 正反互逆测试同一测线两种计算结果的闭合差正算和反算必须对得上这是我最看重的测试。做法很简单先用正算得到终点和反方位角再把这组结果代入反算看得到的距离和起始方位角是否回到原始值。% 正算 [B2, L2, A2] geodesic_direct_iter(B1, L1, A1, S, a, f); % 反算 [S_inv, A1_inv, A2_inv] geodesic_inverse_simple(B1, L1, B2, L2, a, f); fprintf(距离闭合差: %.6f 米\n, S - S_inv); fprintf(方位角闭合差: %.6f 度\n, rad2deg(A1 - A1_inv));理想情况下距离闭合差应在毫米量级方位角闭合差应在 1e-8 度量级。如果用的是简化版正反算闭合差可能会到几十米或 0.01 度这时不要慌说明是长度修正没到位。我见过有人拿着简化版的反算结果去猛调迭代参数最后发现是正算的sigma没收敛白白折腾一下午。互逆测试是检验程序最好的试金石。如果闭合差达不到要求优先怀疑两点一是球心角sigma的求解是否用了迭代二是椭球方位角修正量是否考虑了弧长变化。这两处正确后闭合差自然就降下来了。5. 避坑与排查四个让正反解结果不对劲的常见原因5.1 角度单位错乱弧度与度混用导致的结果偏移现象输入纬度 30输出终点纬度变成 -39 或者跑到地图上完全不对的位置但公式看起来没毛病。原因MATLAB 的三角函数默认输入是弧度而很多课程数据习惯给度。sin(30)实际是sin(30 弧度)不是sin(30°)。贝塞尔公式里经度差dL和方位角都是弧度量一旦混用所有球面三角结果全部乱套。解决在函数入口统一处理。推荐禁止在函数内部做单位切换强制要求调用方传弧度。调用处用deg2rad()转换B1 deg2rad(30); L1 deg2rad(120); A1 deg2rad(60);如果非要支持度输入可以在函数开头加一个标志位判断例如if ang_unit deg B1 deg2rad(B1); A1 deg2rad(A1); end从那以后我每次写这种程序都会在注释第一行写明单位约定避免第二天自己都忘掉。5.2 方位角象限错误atan2用成atan的经典翻车现象计算出的反方位角与正确值相差 180° 或 90°有时候经度方向也反了。原因atan(y/x)无法区分角度在第几象限atan2(y, x)通过两个参数的符号判断象限。很多人写球面方位角公式时图省事写成atan(...)结果在北半球某区域角度正确跨象限后全面翻车。解决所有方位角计算一律用atan2并且保持y参数在前、x参数在后。以起点球面方位角为例alpha1 atan2(cos(u2) * sin(dL), ... cos(u1) * sin(u2) - sin(u1) * cos(u2) * cos(dL));第一个参数是“分子”第二个是“分母”。另外atan2返回范围是[-pi, pi]如果想输出[0, 2*pi]加一行if alpha1 0 alpha1 alpha1 2 * pi; end这一步看似多余但能避免画方位角图时出现负数角度。5.3 迭代不收敛初值距真解太远或容差设错现象正算迭代循环跑满 60 次后S_calc依然和真实S差很远或者sigma不再变化但长度对不上。原因二分法的上下界设得太窄或者长度公式里用了错误的e2。最常见的是把e2 f * (2 - f)写成e2 f^2这会导致平均曲率半径偏大sigma被估算得过小迭代始终向着错误目标压缩区间。解决先单独打印e2和Rm检查数值。WGS84 的e2约等于0.00669437999如果算出0.000011左右说明e2公式错误。另外检查hi上界是否给了pi不要为了加快速度人为缩小上界。还有一点平均曲率半径公式的分母是1 - e2 * sin(Bm)^2如果漏掉平方项长度会系统性偏大几十米。血泪经验迭代不收敛时先查常数再查公式最后才是算法。5.4 互逆性校验失败正算结果再反算对不上现象正算完用反算去还原距离差了几十米方位角差了几百分之一度。看起来“差不多”但拿给导师看被打回。原因简化版正算用sigma S / R0一次折算而反算用S Rm * sigma再折算两次使用的曲率半径所在纬度不同。正算迭代用的是最终中点纬度反算也用了中点纬度理论上应该一致但简化版正算没有迭代导致sigma本身就偏了。解决统一使用完整的迭代正算函数确保S_calc与真实S在毫米量级闭合。互逆测试脚本要写成一个独立的.m文件每次修改代码后自动跑一遍。我一般把闭合差阈值设为 1e-6 米超了就停止执行并报错。6. 从单函数到程序包流程图绘制与封装验收技巧6.1 流程图用哪一层抽象表达正反解流程贝塞尔正反解程序包的流程图不需要画到每个公式画到模块级别就够了。典型流程是输入椭球参数和源数据 - 归化纬度转换 - 选择正算或反算分支 - 球面三角解算 - 长度与方位角修正 - 输出终点或距离方位角。如果你用 MATLAB 的annotation绘制可以把每一步做成一个注释块plot([0 10 10 0 0], [0 0 8 8 0], k); % 外框 text(2, 4, 输入B1,L1,A1,S);实际上更推荐用纸笔画好再截图放进文档MATLAB 画流程图排版很受罪。资源包里的流程图是用通用绘图工具做的直接插入文档改成自己的课程名就能用。6.2 程序包封装把几个函数打包成工具箱的检查清单正式交付前至少要做三件事。第一函数名统一命名正算geodesic_direct反算geodesic_inverse不要出现test1、ceshi这类名字。第二加一个run_demo.m里面包含正算、反算、互逆测试三个调用示例别人拿到手不用看注释也能跑通。第三写一行说明文档写明椭球参数如何切换、角度单位是弧度、所有长度用米。封装完成后用一组独立数据做最终验收。我会用一条 5000 公里的大地线做压测比如从北纬 20° 到南纬 20°横跨多个纬度带。这种长测线能暴露短距离下看不出的方位角修正问题。从那以后我每次写完贝塞尔正反解都强制走一遍长距离互逆闭合差测试先算出毫米级误差再谈交付。希望帮到你。本文还有配套的精品资源点击获取