简介这套由Matlab编写的附和导线网平差程序面向测绘工程专业学生与测量数据处理人员用于导线网观测数据的坐标推算与精度评定支持界面交互导入数据、运行计算并直接查看结果可满足课程设计、毕业设计及实际工程复算等场景需要。压缩包共28个文件以m源码、fig图形界面、txt数据与说明文档为主另有prj工程文件和c接口文件核心算法与界面逻辑分文件存放便于按需学习和改造。包内附带示例数据、结果输出表与运行说明文档下载解压后按说明配置即可在Matlab中运行适合新手及有开发经验的测量人员快速上手。程序经测试校正、运行稳定如需深入理解附合导线闭合差分配、点位误差计算等原理可结合源码逐步研读。目前已有924人浏览学习全包仅72KB是一份精简实用的测量平差参考。1. 外业数据到严密坐标一条附和导线为什么值得写程序外业测量回来手里一沓观测记录转折角闭合差差几秒坐标闭合差差几公分手工近似平差推一上午最后的精度评定还算不出来。这是很多测量内外业人员在没有专用软件时的真实状态。这里讲的是一个基于 MATLAB 实现的附和导线平差程序读入转折角与边长观测值在两端已知点和已知方位角的约束下用间接平差解出全部待定点的坐标并输出单位权中误差、点位中误差和误差椭圆。程序不需要专用工具箱只用矩阵分解和几个基础函数看懂误差方程组装方式后从单条附和导线扩展成导线网也只是数据结构层面的替换。适合测绘工程技术人员、基坑与管线监测数据处理人员以及做测量程序设计课题的学生。2. 附和导线平差的数学模型观测方程、误差方程与权的设置2.1 附和导线的观测值与未知参数多余观测从哪来一条两端定向的附和导线通常有已知起点坐标、已知终点坐标以及起算方位角和终边方位角。中间有n个待定点实测观测值是转折角含两个连接角和边长。以一条只有 2 个待定点的导线为例转折角观测 3 个、边长观测 3 个总计 6 个观测值待定坐标参数是 2 个点 × 2 4 个所以多余观测数 r 6 − 4 2。这 2 个多余观测在手工平差里对应角度闭合差和坐标闭合差两项检核在严密平差中它们不会像近似平差那样先分配角度闭合差、再分配坐标闭合差而是通过最小二乘同时被满足。多余观测数越大越能暴露观测值里的粗差精度评定也越有统计意义。闭合导线、结点导线网的 r 会比单条附和导线大这也是“导线网平差”在数据质量评价上优于单导线手工平差的根本原因。2.2 为什么选间接平差而不是条件平差测量平差有两条经典路线条件平差和间接平差参数平差。条件平差把多余观测数 r 作为条件方程个数每一个条件式描述观测值之间必须满足的几何关系它的优点是未知数少但每个网型的条件方程都要单独列程序化时很难用同一套模板遍历生成。间接平差把所有待定点坐标直接作为参数每个观测值只要写成“观测值 计算值 改正数”的形式就产生一行误差方程。对程序来说逻辑简化为“遍历观测值 → 组装一行系数 → 累加法方程”。单条附和导线和复杂导线网在这个框架下没有本质区别只是参数数量和观测行数变多。所以我一般用间接平差代码量和出错率都更低。2.3 方向角与边长的误差方程先把系数推导写清楚附和导线里边长观测值和坐标的关系是两点间欧氏距离转折角观测值和坐标的关系是两个坐标方位角之差。这两个关系在近似坐标处线性化得到误差方程 v A·dx − l。先看方向角误差方程。约定 X 指北、Y 指东坐标方位角 T 从北方向顺时针起算。对边 ij 有T atan2(ΔY, ΔX)其中 ΔX Xj − XiΔY Yj − Yi。偏导数为参数偏导数 ∂T/∂(参数)XiΔY / S²Yi−ΔX / S²Xj−ΔY / S²YjΔX / S²边长 S 的偏导数是另一组系数参数偏导数 ∂S/∂(参数)Xi−ΔX / SYi−ΔY / SXjΔX / SYjΔY / S在实际代码里我会用函数直接返回这些系数避免每次手推function [az, coef, paramList] az_coeff(Pi, Pj, paramIdx) % 计算坐标方位角及其对两点坐标的偏导数 % paramIdx(i) 是 Pi 点 X 参数编号paramIdx(i)1 是 Y 参数编号 dx Pj(1) - Pi(1); dy Pj(2) - Pi(2); S2 dx*dx dy*dy; S sqrt(S2); az atan2(dy, dx); % 系数顺序: xi, yi, xj, yj coef [dy/S2, -dx/S2, -dy/S2, dx/S2]; paramList [paramIdx(1), paramIdx(1)1, ... paramIdx(2), paramIdx(2)1]; end转折角观测值通常以秒为单位。坐标改正数以米为单位方向角变化量以弧度为单位两者差一个尺度因子 ρ 206264.806″/rad。因此角度误差方程的系数行要整体乘以 ρ常数项 l 用秒表示。这是实现平差程序最容易漏掉的一步不乘 ρ角度行的数值会被边长行淹没法方程数值条件变差结果可能仍然“能算”但精度评定会失真。2.4 角度权与边长权先验精度如何进法方程权的定义是观测值方差倒数的相对比例。取单位权中误差 σ₀角度观测中误差 σβ边长观测中误差 σS则Pβ σ₀² / σβ²PS σ₀² / σS²。σ₀ 取角度先验中误差或 1 都可以因为法方程右边 W 和左边 N 同时被缩放解的坐标不变它影响的是单位权中误差的绝对大小进而不影响相对精度指标。实际项目中σβ 用全站仪标称测角精度比如 ±2″σS 用标称测距精度比如 ±(2mm 2ppm·S)。代码里按每条边逐条计算% 逐边计算边长权 sigmaS 0.002 2e-6 * S_meter(i); % 单位 m P(i) sigma0^2 / sigmaS^2;提示角度误差方程里的 ρ 是弧度与秒的换算不是权的一部分。权公式里的 σβ 仍以秒为单位二者独立设置不要混在一起乘。权对坐标成果的影响通常只有毫米级但对误差椭圆的方向和大小影响显著。如果手头没有检定精度至少按仪器标称值设置比把所有观测值都设成等权更接近真实情况。3. MATLAB 实现数据表结构、近似坐标推算与法方程求解3.1 观测文件用三张表组织点表、边表、角表程序的第一步是确定数据交换格式。我习惯用三个 CSV 文件分别描述控制点、边长观测值和转角观测值文件列说明points.csvid, x, y, knownid 为字符串known 1 表示已知点0 表示待定点edges.csvfrom, to, dist_m边的起点、终点和平距单位米angles.csvstation, back, fore, angle_dms, turn测站点、后视点、前视点、水平角观测值、左角/右角读文件用readtable并给每个点建立从 id 到行号的映射避免后续用字符串找索引pts readtable(points.csv, TextType, string); % 关键把点号映射到矩阵行号 idMap containers.Map(pts.id, 1:height(pts));angle_dms是度分秒格式比如165.3024表示 165°30′24″。转换公式为function deg dms2deg(dms) d floor(dms / 10000); % 度 m floor(mod(dms, 10000) / 100); % 分 s mod(dms, 100); % 秒 deg d m/60 s/3600; end用这个函数而不是直接除以 180 再乘 pi可以避免 60 进制写成 100 进制这类低级错误。观测值进入误差方程前统一转为弧度但角度闭合差的报表输出保留秒方便和外业手簿核对。3.2 近似坐标递推把度分秒观测值变成初始坐标平差要在近似坐标处线性化所以先按导线路线从已知点推一遍近似坐标。左角递推公式为T_next wrap(T_back β_left − π)X_next X_current S·cos(T_next)Y_next Y_current S·sin(T_next)方向角要标准化到 (−π, π]写一个统一函数function T wrap2pi(T) T mod(T pi, 2*pi) - pi; end递推循环按观测顺序进行。已知点本身不作为误差方程参数但参与误差方程系数计算。近似坐标不需要高精度厘米级就够如果从手簿人工录入也可以用坐标正算替代。3.3 误差方程组装关键在参数编号与已知点处理这是整个程序的核心。所有待定点按顺序分配参数编号第 i 个待定点的 X 参数编号为 2i−1Y 为 2i。已知点不占参数编号误差方程里与已知点相关的系数直接乘已知改正数 0等于并入常数项 l。角度观测值的一行误差方程由两个方向角误差方程相减得到。左角 β T(前视方向) − T(后视方向)代码组织如下function [ai, li] angle_row(st, bk, fr, obsSec, coord, paramIdx, rho) % st: 测站点; bk: 后视点; fr: 前视点 % 前视方向误差方程系数 [az_sf, c_sf, p_sf] az_coeff(coord(st,:), coord(fr,:), paramIdx); % 后视方向误差方程系数 [az_bs, c_bs, p_bs] az_coeff(coord(bk,:), coord(st,:), paramIdx); % 计算值: 左角 T(前视) - T(后视) beta_calc wrap2pi(az_sf - az_bs); li obsSec - rad2sec(beta_calc); % 秒 % 系数: 前视方向偏导 减 后视方向偏导再乘 rho nParams max([p_sf, p_bs]); ai zeros(1, nParams); ai(p_sf) ai(p_sf) rho * c_sf; ai(p_bs) ai(p_bs) - rho * c_bs; end右角观测值把相减顺序反过来即可也就是 β T(后视) − T(前视)。我当时处理右角的方法是在读表时加一个布尔量组装前统一转成“等效左角”这样误差方程主体只写一遍。边长观测值的一行更直接function [ai, li] dist_row(fr, to, obsDist, coord, paramIdx) dx coord(to,1) - coord(fr,1); dy coord(to,2) - coord(fr,2); S sqrt(dx*dx dy*dy); % 对起点的偏导为负对终点的偏导为正 coef [-dx/S, -dy/S, dx/S, dy/S]; paramList [paramIdx(fr), paramIdx(fr)1, ... paramIdx(to), paramIdx(to)1]; li obsDist - S; nParams max(paramList); ai zeros(1, nParams); ai(paramList) coef; end把所有角度行和边长行纵向拼接成 A、l再按观测值先验精度组对角矩阵 P最后解法方程N A * P * A; W A * P * l; % 病态检查条件数超过 1e12 时需要回头查参数编号和已知点配置 if cond(N) 1e12 warning(法方程接近病态请检查已知点与观测值配置); end dx N \ W;这里用N \ W而不是inv(N) * W。对 2n×2n 的法方程反斜杠走 LU 分解数值更稳代码也更短。3.4 解算与迭代N\W 代替 inv(N)收敛条件看 dx线性化误差方程在近似坐标处展开如果近似坐标有偏差理论上需要重新线性化。对附和导线近似坐标由观测值递推通常误差在秒级和毫米级一次解算已经足够。更稳妥的做法是迭代到改正数小到可忽略for iter 1:5 % 在当前 coord 下重新组装 A、l [A, l, P] build_system(pts, edges, angles, coord, paramIdx); dx (A*P*A) \ (A*P*l); coord update_coord(coord, dx, paramIdx); if norm(dx) 1e-4 break; end end收敛阈值取 1e−4 米即 0.1 毫米。这样既避免第一次模型线性化误差影响成果又不会在亚毫米级别空转。最后一次迭代的 A、l、P 保留下来供精度评定使用——这一点容易被忽略如果用第一次的 A 和最后一次的坐标组合算改正数v A·dx − l 会出现不一致后续统计量全部失真。4. 精度评定单位权中误差、点位中误差与误差椭圆输出4.1 单位权中误差vTPv 与自由度的对应关系平差结束后先计算改正数向量 v A·dx − l再统计观测值内部符合程度v A * dx - l; r size(A,1) - size(A,2); % 自由度: 观测数 - 参数数 sigma02 (v * P * v) / r; % 单位权方差 sigma0 sqrt(sigma02);自由度 r 是 6 个观测值减去 4 个坐标参数的例子就是 2。如果 r 小于等于 0精度评定没有统计意义程序应当直接报错而不是输出一组看起来正常的中误差。单位权中误差的数值应与先验 σ₀ 同量级若明显偏大说明先验精度设置过于乐观或观测值中存在粗差。4.2 协因数阵与点位中误差Qxx 对角线的物理含义坐标参数的协因数阵是 Qxx N⁻¹点位协方差阵为 Σ σ₀²·Qxx。第 i 个待定点的平面位置属于它的两行两列子块Qxx inv(A * P * A); sigma0_2 sigma02; for k 1:nUnknown Qk Qxx(2*k-1 : 2*k, 2*k-1 : 2*k) * sigma0_2; mx sqrt(Qk(1,1)); % X 方向中误差单位 m my sqrt(Qk(2,2)); % Y 方向中误差 mp sqrt(Qk(1,1) Qk(2,2)); % 点位中误差 end这里必须把单位权方差乘进去。直接把 Qxx 对角线开方当作业中误差是国内教学代码里最常见的错误。4.3 误差椭圆特征值分解不是可选项点位中误差只给出一个标量误差椭圆则描述该点在各个方向上的不确定性。对每个点的协方差子矩阵做特征值分解[V, D] eig(Qk); [ev, ord] sort(diag(D), descend); % 特征值降序 a sqrt(ev(1)); % 长半轴 b sqrt(ev(2)); % 短半轴 theta atan2(V(2, ord(1)), V(1, ord(1))) * 180/pi; % 长轴方位角特征值对应协方差矩阵主方向特征向量给出长轴方向。误差椭圆长半轴的方向就是该点位误差最大的方向往往垂直于导线边方向这在布网设计和控制点选取时很有用。画图不需要额外工具箱用 MATLAB 基础画图函数即可先按极角生成椭圆点列再乘以长短半轴并旋转到方位角叠加到导线图上。MATLAB 画图语法本身很简单核心是把误差椭圆放到每个待定点上而不是单独画一个形状。4.4 把精度结果写成报告点号、误差椭圆与相对点位精度除单点精度外相邻待定点的相对精度也对放样和监测有实际意义。相对点位误差要从 Qxx 中同时取两点的行和列Qpair Qxx([2*i-1, 2*i, 2*j-1, 2*j], ... [2*i-1, 2*i, 2*j-1, 2*j]) * sigma0_2; dX [1, 0, -1, 0]; dY [0, 1, 0, -1]; mRelX sqrt(dX * Qpair * dX); % 两点间 X 方向相对中误差 mRelY sqrt(dY * Qpair * dY);最终把结果导出为 CSV便于直接用于报告审阅点号X(m)Y(m)mx(mm)my(mm)点位中误差(mm)长半轴(mm)短半轴(mm)长轴方位角(°)P11101.23452003.45672.13.23.84.01.967.3这个表格的数值只是示意输出。实际程序中用writetable输出即可点号保持字符串类型坐标保留到毫米位中误差保留到 0.1 毫米。5. 验证与排错用模拟数据确认程序可信再上真实外业5.1 模拟自检用“真值坐标反算观测值”做回归平差程序最隐蔽的错误往往不是语法而是系数符号反了、已知点编号错位、角度方向反了。这些错误在真实数据上表现出的闭合差可能很小让人误以为结果正确。我拿到任何新写的平差程序第一件事就是做模拟自检。构造一个简单附和导线先给定一组“真值坐标”用坐标反算出每条边长和每个转折角的真值给观测值加上毫米级随机噪声再让程序平差把平差坐标对比真值坐标。流程只有四步% 第一步给定真值坐标 truth [1000, 2000; 1101.234, 2003.456; ...]; % 第二步由坐标反算边长与方向角生成观测值真值 for e each edge obsDist norm(coord(to,:) - coord(fr,:)); end % 第三步叠加已知大小的随机噪声 rng(2025); % 固定种子保证可复现 obsDist obsDist randn(size(obsDist)) * 0.003; % 第四步平差后对比 mismatch coord_result - truth; assert(max(abs(mismatch(:))) 0.02);固定随机种子rng(2025)是关键。这样每一次修改程序后重新跑得到的结果完全可比较而不是因为随机噪声不同而误判程序出了问题。把模拟自检和真实数据平差分开模拟数据负责验证程序逻辑真实数据负责验证外业质量。两者混在一起时出现异常很难定位是观测值质量差还是程序有 bug。5.2 MATLAB 里最容易踩的五个坑按出现频率排序这些坑我基本都踩过一遍度分秒转十进制的进制错误。165.3024是 165°30′24″不是 165.3024°。转换函数里要按度、分、秒分别取整不能用mod(dms, 1)之类的方式。atan2参数顺序写反。坐标方位角要写成atan2(dy, dx)其中 dx 是北方向坐标差。写成atan2(dx, dy)会把东西和南北互换结果在方向角接近 0° 和 90° 时很难一眼看出来。已知点的系数塞进了参数列。已知点改正数为 0系数要么不进入 A要么在进入前乘 0。最直观的检查法方程 N 的维数应该等于2 * 待定点数而不是2 * 总点数。权矩阵只给角度不给边长或者量纲不一致。角度权要和角度秒值配套边长权要和米值配套单位不统一时误差椭圆的长轴方向会明显偏离导线边的垂直方向。迭代过程中 A、l 没有同步更新。只更新坐标不重新组装 A改正数序列会在迭代后出现不合理的跳跃精度统计全是错的。5.3 从单导线到简单导线网只改数据结构不改平差核心标题里带“导线网平差”实现层面和单条附和导线的差异没有想象中那么大。单导线是“一条路线顺序编号”导线网则是多个结点和闭合环共享同一套参数编号空间。把待定点放进一个全局列表每个点只分配一次参数编号所有角度和边长观测值按同一个规则逐行加入 A 矩阵误差方程的逻辑完全复用。唯一要额外处理的是起算数据已知方位角在单导线里作为方向起算值导线网中需要明确哪些边带有已知方位角避免重复使用或漏用。参数编号统一后单导线的angle_row和dist_row可以直接在导线网中使用。这也是我强调“参数编号与 id 解耦”的原因导线网的代码和单导线的代码共用同一个build_system差别只在数据表内容。写一个模拟数据、固定随机种子的回归脚本每次改动后跑一遍并与上次结果比对坐标差异这个方法比任何人工检查都更能兜住程序改动的风险。本文还有配套的精品资源点击获取