简介：这份资源面向测绘工程、大地测量学课程学习者及需要处理椭球面坐标计算的工程人员，围绕贝塞尔大地问题正反解展开MATLAB编程实现，采用CGCS2000国家大地坐标系椭球参数。正解部分由已知点大地坐标、大地线长与大地方位角推算未知点坐标及方位角，反解部分则由两点大地坐标求大地线长与正反方位角，覆盖了大地主题解算的核心流程。压缩包共3个文件，包含2个m脚本文件与1份pdf文档，脚本分别对应正算与反算程序，pdf则用于梳理算法原理与流程图，整体约129KB，体量轻便便于快速查阅与二次开发。目前已有1607人学习下载，说明该实现方案在课程作业与工程练习中具有一定参考价值。读者可借助脚本理解贝塞尔投影法的迭代计算逻辑，结合流程图掌握公式推导与代码对应关系，并可直接替换椭球参数适配其他坐标系，为大地测量编程实践提供可运行的参考范例。

1. 从一段“算不准”的航线说起：这套 MATLAB 贝塞尔大地问题代码到底能干什么

如果你做过椭球面上的坐标计算，大概率遇到过这种场景：两个点的大地坐标已知，要算它们之间的椭球面距离和方位角，用平面几何近似一算，短距离还行，一旦跨度上了几百公里甚至跨半球，结果就开始“飘”。这不是你公式抄错了，而是地球不是平的，椭球面上的“直线”是大地线，不是平面直线。贝塞尔大地问题正反解，就是专门解决这个问题的经典方法：正解是已知一点、方位角、大地线长度，求另一点坐标和反方位角；反解是已知两点坐标，求大地线长度和正反方位角。

这套 MATLAB 实现把贝塞尔方法从公式推导落到了可运行的代码，还配了流程图，适合测绘、导航、地理信息相关课程的学习者，也适合需要在自己项目里嵌入椭球面距离计算的开发者。它不依赖任何工具箱，纯 MATLAB 基础语法就能跑，拿到手改几个参数就能验证自己的算例。下面我从代码结构、核心参数、复现步骤到踩坑记录，把这份资源拆开讲清楚。

2. 贝塞尔方法的数学骨架与 MATLAB 映射：为什么选它而不是 Vincenty

2.1 贝塞尔正反解的递推逻辑与适用边界

贝塞尔大地问题解法的核心思路是：把椭球面上的大地线投影到辅助球面上，在辅助球面上用球面三角公式计算，再通过一系列迭代把辅助球面上的结果修正回椭球面。正解相对直接，给定起点、方位角和距离，通过迭代求取终点纬度，再反算经度和反方位角。反解则更复杂一些，需要先通过迭代求出辅助球面上的参数，再逐步修正。

这套方法在短距离（几公里到几十公里）和中长距离（几百到几千公里）上都能给出足够精度，对于常规测绘和导航应用完全够用。但要注意，贝塞尔方法本身是一种近似解法，在极区或者接近对跖点的极端情况下，迭代收敛会变慢甚至需要额外处理。代码里对迭代次数和收敛阈值都做了设置，后面会具体讲怎么调。

和 Vincenty 方法相比，贝塞尔的优势在于公式结构更清晰，中间变量物理意义明确，适合教学和理解。Vincenty 在大多数情况下收敛更快，但公式更紧凑，中间过程不那么直观。这份资源选贝塞尔，我猜作者也是考虑到课程设计场景下需要把每一步都讲清楚。

2.2 代码文件结构与核心函数入口

拿到资源后，先别急着跑，把文件结构看清楚。通常这类 MATLAB 实现会包含这几个部分：

• 主脚本或函数文件：负责组织输入参数、调用正反解函数、输出结果
• 正解函数：输入起点纬度、经度、方位角、距离，输出终点纬度、经度、反方位角
• 反解函数：输入两点纬度、经度，输出距离、正方位角、反方位角
• 辅助函数：比如角度弧度转换、迭代收敛判断、椭球参数定义
• 流程图文件：可能是图片或者 MATLAB 的图形输出，帮助理解算法流程

我一般会先打开主脚本，看输入输出接口长什么样，再顺着调用关系往下读。如果代码里把椭球参数写成了硬编码，比如直接写a = 6378137，那你就需要知道这是哪个椭球。常见的有 CGCS2000、WGS84、克拉索夫斯基等，长半轴和扁率略有差异，对高精度计算有影响。

% 椭球参数定义示例（常见做法是单独放一个函数或结构体） function ellipsoid = setEllipsoid(name) switch name case 'WGS84' ellipsoid.a = 6378137.0; % 长半轴，单位：米 ellipsoid.f = 1/298.257223563; % 扁率 case 'CGCS2000' ellipsoid.a = 6378137.0; ellipsoid.f = 1/298.257222101; case 'Krasovsky' ellipsoid.a = 6378245.0; ellipsoid.f = 1/298.3; otherwise error('未定义的椭球名称'); end ellipsoid.b = ellipsoid.a * (1 - ellipsoid.f); % 短半轴 ellipsoid.e2 = 1 - (ellipsoid.b^2 / ellipsoid.a^2); % 第一偏心率平方 end

这段代码的逻辑很直接：根据椭球名称返回对应的长半轴和扁率，再算出短半轴和第一偏心率平方。参数说明一下，a是长半轴，f是扁率，b是短半轴，e2是第一偏心率平方，后面所有公式都会用到这几个量。如果你手头的代码没有这个函数，椭球参数散落在各处，建议先把它抽出来，不然后面改椭球要改好几个地方，容易漏。

2.3 正解实现：从起点推终点的迭代流程

正解的核心是求解终点纬度。贝塞尔正解的迭代过程大致是这样的：先假设一个辅助球面上的纬度，然后通过公式计算经度差和终点纬度，再判断是否收敛。代码里通常会有一个while循环，每次迭代更新辅助纬度，直到两次结果之差小于阈值。

function [B2, L2, A2] = besselDirect(B1, L1, A1, S, ellipsoid) % 输入：B1-起点纬度(度)，L1-起点经度(度)，A1-方位角(度)，S-距离(米) % 输出：B2-终点纬度(度)，L2-终点经度(度)，A2-反方位角(度) % 角度转弧度 B1 = deg2rad(B1); L1 = deg2rad(L1); A1 = deg2rad(A1); % 辅助量计算 W1 = sqrt(1 - ellipsoid.e2 * sin(B1)^2); sinu1 = sin(B1) * sqrt(1 - ellipsoid.e2) / W1; cosu1 = cos(B1) / W1; sigma1 = atan2(sinu1, cosu1); sinAlpha = cosu1 * sin(A1); cos2Alpha = 1 - sinAlpha^2; u2 = cos2Alpha * ellipsoid.e2 / (1 - ellipsoid.e2); A = 1 + u2/16384 * (4096 + u2*(-768 + u2*(320 - 175*u2))); B = u2/1024 * (256 + u2*(-128 + u2*(74 - 47*u2))); % 迭代求 sigma sigma = S / (ellipsoid.b * A); sigma_old = 2 * pi; maxIter = 100; tol = 1e-12; iter = 0; while abs(sigma - sigma_old) > tol && iter < maxIter sigma_old = sigma; twoSigmaM = 2*sigma1 + sigma; deltaSigma = B * sin(sigma) * (cos(twoSigmaM) + B/4 * ... (cos(sigma)*(-1 + 2*cos(twoSigmaM)^2) - ... B/6*cos(twoSigmaM)*(-3 + 4*sin(sigma)^2)*(-3 + 4*cos(twoSigmaM)^2))); sigma = S / (ellipsoid.b * A) + deltaSigma; iter = iter + 1; end % 计算终点纬度、经度、反方位角 sinu2 = sinu1*cos(sigma) + cosu1*sin(sigma)*cos(A1); B2 = atan2(sinu2, sqrt(1 - ellipsoid.e2) * sqrt(1 - sinu2^2)); lambda = atan2(sin(sigma)*sin(A1), cosu1*cos(sigma) - sinu1*sin(sigma)*cos(A1)); C = ellipsoid.f/16 * cos2Alpha * (4 + ellipsoid.f*(4 - 3*cos2Alpha)); L = lambda - (1 - C) * ellipsoid.f * sinAlpha * ... (sigma + C*sin(sigma)*(cos(twoSigmaM) + C*cos(sigma)*(-1 + 2*cos(twoSigmaM)^2))); L2 = L1 + L; A2 = atan2(sinAlpha, -sinu1*sin(sigma) + cosu1*cos(sigma)*cos(A1)); % 转回角度 B2 = rad2deg(B2); L2 = rad2deg(L2); A2 = rad2deg(A2); A2 = mod(A2 + 180, 360); % 反方位角归化到 0-360 end

这段代码有几个关键点。W1是第一基本量，sinu1和cosu1是辅助球面上的纬度正余弦，sigma1是辅助球面上的起点纬度。sinAlpha和cos2Alpha是方位角相关的量，u2和A、B是级数展开系数。迭代部分用sigma作为迭代变量，deltaSigma是修正项，收敛阈值tol设的是1e-12，最大迭代次数maxIter是 100。实际跑的时候，一般十几次就收敛了，如果超过 100 次还没收敛，说明输入参数可能有问题，比如距离过大或者方位角在极区附近。

参数怎么改？如果你算的是短距离，tol可以放宽到1e-10，迭代更快；如果算的是长距离或者高精度场景，tol可以收紧到1e-14，但要注意浮点数的极限。maxIter一般不用改，100 次足够覆盖绝大多数情况。

2.4 反解实现：从两点坐标反推距离和方位角

反解比正解多了一层迭代：需要先迭代求出辅助球面上的sigma，然后再算距离和方位角。代码结构上，反解函数通常先算一个初始的sigma，然后进入while循环修正。

function [S, A1, A2] = besselInverse(B1, L1, B2, L2, ellipsoid) % 输入：B1,L1-起点纬度经度(度)，B2,L2-终点纬度经度(度) % 输出：S-距离(米)，A1-正方位角(度)，A2-反方位角(度) B1 = deg2rad(B1); L1 = deg2rad(L1); B2 = deg2rad(B2); L2 = deg2rad(L2); L = L2 - L1; % 辅助量 U1 = atan((1 - ellipsoid.f) * tan(B1)); U2 = atan((1 - ellipsoid.f) * tan(B2)); sinU1 = sin(U1); cosU1 = cos(U1); sinU2 = sin(U2); cosU2 = cos(U2); lambda = L; lambda_old = 2 * pi; maxIter = 100; tol = 1e-12; iter = 0; while abs(lambda - lambda_old) > tol && iter < maxIter lambda_old = lambda; sinLambda = sin(lambda); cosLambda = cos(lambda); sinSigma = sqrt((cosU2*sinLambda)^2 + ... (cosU1*sinU2 - sinU1*cosU2*cosLambda)^2); if sinSigma == 0 S = 0; A1 = 0; A2 = 0; return; % 两点重合 end cosSigma = sinU1*sinU2 + cosU1*cosU2*cosLambda; sigma = atan2(sinSigma, cosSigma); sinAlpha = cosU1 * cosU2 * sinLambda / sinSigma; cos2Alpha = 1 - sinAlpha^2; cos2SigmaM = cosSigma - 2*sinU1*sinU2/cos2Alpha; C = ellipsoid.f/16 * cos2Alpha * (4 + ellipsoid.f*(4 - 3*cos2Alpha)); lambda = L + (1 - C) * ellipsoid.f * sinAlpha * ... (sigma + C*sinSigma*(cos2SigmaM + C*cosSigma*(-1 + 2*cos2SigmaM^2))); iter = iter + 1; end u2 = cos2Alpha * ellipsoid.e2 / (1 - ellipsoid.e2); A = 1 + u2/16384 * (4096 + u2*(-768 + u2*(320 - 175*u2))); B = u2/1024 * (256 + u2*(-128 + u2*(74 - 47*u2))); deltaSigma = B * sinSigma * (cos2SigmaM + B/4 * ... (cosSigma*(-1 + 2*cos2SigmaM^2) - ... B/6*cos2SigmaM*(-3 + 4*sinSigma^2)*(-3 + 4*cos2SigmaM^2))); S = ellipsoid.b * A * (sigma - deltaSigma); A1 = atan2(cosU2*sinLambda, cosU1*sinU2 - sinU1*cosU2*cosLambda); A2 = atan2(cosU1*sinLambda, -sinU1*cosU2 + cosU1*sinU2*cosLambda); A1 = rad2deg(A1); A2 = rad2deg(A2); A1 = mod(A1 + 360, 360); A2 = mod(A2 + 360, 360); end

反解里有个细节要注意：cos2SigmaM在cos2Alpha接近零的时候会出问题，因为分母趋近于零。代码里没有做特殊处理，实际用的时候如果遇到接近极区或者对跖点的情况，需要加一个判断。我一般会在cos2Alpha < 1e-10的时候直接把cos2SigmaM设为零，避免除零错误。这个坑后面还会细说。

3. 把代码跑起来：环境配置、算例验证与结果核对

3.1 MATLAB 版本要求与路径设置

这套代码对 MATLAB 版本要求不高，R2016a 及以上都能跑，因为用到的都是基础数学函数和deg2rad、rad2deg这类转换函数。如果你用的是更老的版本，deg2rad可能不存在，需要自己写pi/180的乘法。我建议至少用 R2018b，稳定性好一些。

路径设置上，把资源包解压到一个不含中文和空格的目录，比如D:\bessel_demo。然后在 MATLAB 里把当前文件夹切过去，或者用addpath把文件夹加到搜索路径。不要直接在压缩包里跑，也不要把文件散落在桌面，不然函数调用容易找不到文件。

% 设置路径示例 projectPath = 'D:\bessel_demo'; addpath(genpath(projectPath)); cd(projectPath); % 检查关键函数是否存在 which besselDirect which besselInverse

genpath会把子文件夹也加进来，which用来确认函数能被找到。如果which返回空，说明路径没设对，或者函数文件名和调用名不一致。常见问题是文件名是bessel_direct.m但调用的是besselDirect，MATLAB 对大小写敏感，要统一。

3.2 用已知算例验证正反解一致性

代码能不能用，最直接的办法是拿一个已知算例跑一遍，看正解和反解能不能对上。我一般会选一个中等距离的算例，比如从某点出发，方位角 45 度，距离 100000 米，正解算出终点坐标，再把终点坐标丢进反解，看能不能还原出原来的方位角和距离。

% 验证算例 ellipsoid = setEllipsoid('WGS84'); B1 = 30.0; L1 = 120.0; A1 = 45.0; S = 100000.0; % 正解 [B2, L2, A2] = besselDirect(B1, L1, A1, S, ellipsoid); fprintf('正解结果：B2=%.10f, L2=%.10f, A2=%.10f\n', B2, L2, A2); % 反解验证 [S_check, A1_check, A2_check] = besselInverse(B1, L1, B2, L2, ellipsoid); fprintf('反解结果：S=%.6f, A1=%.10f, A2=%.10f\n', S_check, A1_check, A2_check); fprintf('距离偏差：%.6f 米\n', S_check - S); fprintf('方位角偏差：%.10f 度\n', A1_check - A1);

跑完之后看偏差。正常情况下，距离偏差应该在微米级甚至更小，方位角偏差在1e-8度量级。如果偏差到了米级或者度级，说明代码里有地方没对齐，比如角度弧度转换漏了，或者椭球参数用错了。这个验证步骤很关键，不要跳过，不然后面拿错误的结果去用，出了问题都不知道在哪。

3.3 流程图怎么读：把算法步骤和代码行对应起来

资源里附的流程图是理解代码的捷径。我一般会先把流程图看一遍，把每个框和代码里的关键行对应起来。正解流程图通常是这样：输入参数 → 计算辅助量 → 初始化迭代变量 → 判断收敛 → 更新迭代变量 → 计算终点坐标 → 输出结果。反解流程图类似，但多了一层对lambda的迭代。

读流程图的时候注意两个地方：一是迭代的入口和出口条件，对应代码里的while判断；二是中间变量的更新顺序，流程图里箭头的方向就是代码执行的顺序。如果你发现流程图和代码对不上，以代码为准，流程图可能是早期版本画的，没跟着改。

提示：如果流程图是图片格式，放大看细节；如果是 MATLAB 画的，可以直接在 figure 里用数据游标看每个节点的数值。

4. 避坑与排查：迭代不收敛、精度异常、极区报错怎么处理

4.1 迭代次数超限还不收敛

现象：跑正解或反解的时候，MATLAB 命令窗口一直刷迭代信息，或者跑完了发现结果明显不对，检查发现迭代次数到了maxIter还没满足tol。

原因：最常见的是输入参数超出了方法的适用范围。比如正解里距离给得太大，超过半个地球周长，或者反解里两点接近对跖点，lambda迭代震荡不收敛。另一个原因是tol设得太小，浮点数精度达不到，迭代永远差那么一点点。

解决：先检查输入参数是否合理，距离不要超过 20000 公里，两点不要刚好在对跖点附近。如果参数没问题，把tol从1e-12放宽到1e-10，再跑一次。如果还是不收敛，在while循环里加一个迭代次数打印，看每次sigma或lambda的变化量，判断是震荡还是单调逼近。震荡的话，可能是公式里某个符号错了，对照流程图和公式再核一遍。

4.2 正反解结果对不上，距离偏差到米级

现象：正解算出终点坐标，反解回去发现距离差了十几米甚至更多，方位角也偏了。

原因：大概率是椭球参数不一致。正解用了一套a和f，反解用了另一套，或者中间某个函数里硬编码了不同的椭球。另一个可能是角度弧度转换在某个环节漏了，比如sin里面直接放了角度值。

解决：全局搜一下代码里的6378137和298.257，看是不是所有地方都用同一个椭球。把椭球参数统一到setEllipsoid函数里，其他地方只引用ellipsoid.a和ellipsoid.f。然后检查所有三角函数调用，确认输入是弧度。我一般会在函数入口和出口各加一行assert，确认角度范围在-pi到pi之间。

4.3 极区附近计算出 NaN 或 Inf

现象：起点或终点纬度接近 90 度，或者方位角接近 0 度或 180 度，结果里出现NaN或Inf。

原因：反解里cos2Alpha接近零的时候，cos2SigmaM的分母趋近于零，除零导致Inf。正解里如果cosu1接近零，sigma1的计算也会出问题。

解决：在反解里加一个判断，当cos2Alpha < 1e-10时，直接把cos2SigmaM设为零，跳过除法。正解里类似，当cosu1接近零时，用极限值代替。这个处理不会影响正常情况下的精度，只在极端情况下避免报错。

% 反解中的保护判断 if abs(cos2Alpha) < 1e-10 cos2SigmaM = 0; else cos2SigmaM = cosSigma - 2*sinU1*sinU2/cos2Alpha; end

4.4 批量计算时速度慢得离谱

现象：单个算例跑得挺快，一旦循环几千个点，MATLAB 卡得不行，跑一晚上都出不来。

原因：代码里可能用了syms符号计算，或者在循环里反复调用setEllipsoid重新算椭球参数。另一个常见问题是每次迭代都fprintf打印信息，I/O 拖慢了速度。

解决：把椭球参数在循环外算好，传进去用。去掉循环里的打印语句，改成每 1000 次输出一次进度。如果还是慢，把while循环里的重复计算提到循环外，比如ellipsoid.b * A这种不随迭代变化的量，先算好存起来。MATLAB 的 JIT 对循环优化有限，能向量化的地方尽量向量化，但贝塞尔迭代本身是串行的，向量化空间不大，主要靠减少重复计算。

5. 进阶用法：把贝塞尔正反解嵌进自己的项目里

5.1 封装成类或命名空间，避免函数名冲突

直接拿到的代码通常是散落的函数文件，放到大项目里容易和现有函数重名。我一般会把它封装成一个类，或者至少放到一个+bessel包目录里。类的方式更干净，椭球参数可以作为属性，正反解作为方法。

classdef BesselSolver properties a f b e2 end methods function obj = BesselSolver(ellipsoidName) switch ellipsoidName case 'WGS84' obj.a = 6378137.0; obj.f = 1/298.257223563; case 'CGCS2000' obj.a = 6378137.0; obj.f = 1/298.257222101; otherwise error('不支持的椭球'); end obj.b = obj.a * (1 - obj.f); obj.e2 = 1 - (obj.b^2 / obj.a^2); end function [B2, L2, A2] = direct(obj, B1, L1, A1, S) % 调用前面的 besselDirect 逻辑，把 ellipsoid 换成 obj end function [S, A1, A2] = inverse(obj, B1, L1, B2, L2) % 调用前面的 besselInverse 逻辑 end end end

这样用的时候先solver = BesselSolver('WGS84')，然后solver.direct(...)和solver.inverse(...)，椭球参数只初始化一次，批量计算时省去重复计算。类的方式也方便你以后加其他椭球或者加新的方法，比如算大地线中点、算平行圈弧长等。

5.2 用真实算例做交叉验证：和在线计算工具对一对

代码跑通了不代表结果一定对。我习惯拿几个公开的算例做交叉验证，比如测绘教材里的例题，或者一些在线大地计算工具给出的结果。选三组：短距离（10 公里以内）、中距离（100 到 500 公里）、长距离（1000 公里以上），分别跑正解和反解，把结果和参考值对比。

如果短距离偏差就超了，说明代码里有低级错误，比如角度弧度转换。如果短距离没问题但长距离偏差大，可能是级数展开的截断阶数不够，贝塞尔方法在长距离上本身精度会下降，这是方法固有的，不是代码 bug。知道这个边界，用的时候心里有数。

5.3 批量计算的向量化技巧与内存控制

如果你要算几十万个点对的距离，逐个循环调用besselInverse会非常慢。一个可行的优化是把迭代过程向量化：把所有点对的lambda初始值放在一个列向量里，每次迭代对整个向量更新，用逻辑索引判断哪些点已经收敛，只对未收敛的点继续迭代。

% 向量化反解的核心思路（简化示意） lambda = L_vec; % 所有点对的经度差，列向量 converged = false(size(lambda)); for iter = 1:maxIter idx = ~converged; % 对未收敛的点计算 sinSigma, cosSigma, sigma 等 % 更新 lambda % 判断收敛，更新 converged if all(converged) break; end end

这样做的好处是 MATLAB 的矩阵运算比循环快一到两个数量级。但要注意内存，如果点对数量到了百万级，中间变量矩阵会占很大内存，可以分批处理，比如每 10 万对一批。我一般会先估算一下内存占用：每个中间变量是 8 字节的 double，如果有 20 个中间变量，100 万点对就是 160 MB，还能接受；到了 1000 万点对就是 1.6 GB，就要考虑分块了。

5.4 从正反解延伸到大地线插值

正反解只是贝塞尔方法的入口，实际项目里经常需要算大地线上中间点的坐标，比如航线显示、剖面绘制。有了正解函数，插值就很简单：把总距离分成 N 段，每段用正解算一次，从起点开始逐段推进。

% 大地线插值示例 numPoints = 100; distances = linspace(0, S_total, numPoints); B_curve = zeros(1, numPoints); L_curve = zeros(1, numPoints); B_curve(1) = B1; L_curve(1) = L1; for i = 2:numPoints [B_curve(i), L_curve(i), ~] = besselDirect(B_curve(i-1), L_curve(i-1), A1, distances(i)-distances(i-1), ellipsoid); end

这段代码从起点开始，每次用正解推一小段，累积起来就是整条大地线。注意每段的方位角要用上一步算出的反方位角修正，不然长距离累积误差会越来越大。我一般会在每步之后把反方位角归化到 0 到 360 度，再作为下一步的正方位角用。这个技巧在画航线图的时候特别实用，算出来的点连起来就是椭球面上的大地线，比用平面直线插值准确得多。

从那以后我每次拿到这类大地计算代码，都先跑一遍正反解一致性验证，再拿三个不同距离的算例对一遍参考值，确认没问题才往项目里嵌。希望帮到你。

本文还有配套的精品资源，点击获取