★ FEATURED ARTICLE 简介导线网平差程序是一套面向武汉大学测绘学院平差课程设计的C# WPF实现完整覆盖导线网近似值计算、角度与边长观测值处理、平差值迭代、坐标改正数求解以及验后单位权中误差、坐标参数方差阵、点位中误差和边角观测值改正数等精度评定指标并可绘制导线网图进行可视化展示。RAR资源包内共80个文件、约5.98MB核心代码为16个C#源码文件另有17个txt数据与结果文本、3个xaml界面文件、6个dll运行库以及项目配置文件还附带两组测试数据和郭辉《导线网间接平差算法分析与实现》PDF论文方便对照验证。目前已有3790人学习下载特别适合正在准备武汉大学测绘学院平差课设、需要快速拿到可运行程序并理解间接平差实现细节的本科生。代码结构清晰Adjustment.cs与Data.cs为计算核心类既可整体运行也可提取迁移到新的WPF或WinForm工程自行构建视图再按按钮功能调用相应函数完成课设内置的DisplayWindow.xaml也给出了显示窗口的参考实现能帮助读者节省算法框架搭建时间将重点放在测量平差原理与数据验证上。1. 导线网平差程序为什么实习代码写了两遍还是差一大截导线网平差程序是测绘工程专业实习里绕不开的一件硬活很多同学拿着教材里的公式手算了一遍又一遍最后结果和老师给的答案差出好几厘米只能归因于“观测数据的玄学误差”。实际上导线网平差的原理并不复杂把水平角和边长观测值组织成误差方程用间接平差方法解一个线性方程组几百行源代码就能跑通整个流程。这篇文章按我从零写实习程序的思路来讲先说明为什么选间接平差而不是条件平差再落到观测文件怎么设计、误差方程怎么组、法方程怎么解最后把新手最容易翻车的五个坑挑出来。适合正在赶平差实习报告的你也适合想用程序替代手算的现场工程师。2. 导线网平差的数学模型为什么选间接平差而不是条件平差2.1 导线网里到底有哪些观测值和未知数导线网平差处理的观测值主要两类水平角和边长。水平角来自全站仪测回观测边长来自测距或全站仪斜距换算后的平距。待求的未知数很简单就是每个待定点的平面坐标改正数一个点有 (dx, dy) 两个分量。已知点坐标是起算数据不进入未知数列表。设观测值总数为 n必要观测数为 t多余观测数 r n - t。导线网里 r 必须大于等于 0等于 0 时就是导线计算没有平差可言大于 0 时才有“用多余观测约束误差”的空间。实习里常见的闭合导线、附合导线都有 2 到 5 个多余观测所以平差程序的意义就是把这一部分冗余信息转化为精度提升。这里有个容易被忽略的前提观测值和未知数不是一一对应的。一个水平角观测值同时关联三个点的坐标改正数——测站点、后视点、前视点一条边长观测值关联两个点。写程序时如果只按“一个观测值一行数据”来索引后面组误差方程矩阵时就会漏系数。2.2 条件平差与间接平差为什么程序代码只选间接平差平差的经典模型有两套。条件平差要列出闭合条件方程比如闭合导线内角和条件、坐标增量闭合条件再把观测值改正数投影到条件方程上求解。对于单个闭合导线条件方程数目少手算很合适。但导线网一旦有多个节点、多条闭合环条件方程怎么列全就成了最大的负担列漏一条整个程序就悄悄偏掉而且这种偏错很难在结果上第一时间察觉。间接平差参数平差的思路完全不同直接选待定点坐标为参数对每一个观测值单独列一个误差方程形式统一为“观测值 改正数 观测值的函数”。函数是非线性的所以要在近似坐标处线性化最后解线性化的法方程。代码上就是一个 for 循环遍历所有观测值逐行填充矩阵 A 和常数项 l完全不需要人工判断网形。实习代码我一般强烈建议选间接平差没有第二种选择。两种模型解出来的平差结果是等价的精度评定结果也在数值上一致。但代码实现的复杂度和调试难度差异极大。条件平差的未知数是观测值改正数间接平差的未知数是坐标改正数后者更直观算完直接加到坐标上就是平差坐标。2.3 定权策略测角中误差与测距中误差怎么换算成权平差里的“权”代表观测值的相对可信度权越大该观测值在解算中的话语权越大。定权不当是实习程序结果偏差的最大来源之一。常见做法是以测角中误差为基准把测边观测值的权换算成与角度相同的量纲级别。设测角中误差为 σβ单位弧度测边中误差为 σs单位米则角度观测值的权 p_angle 1 / σβ^2边长观测值的权 p_dist 1 / σs^2。直接这么算出来的权在数值上可能差很多个数量级但不影响解算结果因为法方程会通过矩阵运算自动平衡。重要的是权之间的比例关系要对比如角度中误差 1 秒边长中误差 2 毫米那么边长权的相对值约为角度权的 (1″ / 2mm)^2换算到统一单位后差约 3.6×10^7 倍。在代码里权阵 P 通常按对角阵组织对角线元素按观测值顺序排列。这样法方程 N A^T P A 和 U A^T P l 的计算就是两次矩阵乘法代码非常整齐。如果实习只要求模拟数据可以把角度权设为 1边长权按比例系数换算后填入这就是后面验证程序时最关键的输入参数之一。3. 观测数据与已知点的组织先把文件格式设计成不会翻车的样子3.1 观测文件格式设计一行一条观测值谁都能看懂写平差程序第一步不是写算法而是设计观测文件的读写格式。很多同学把观测值硬编码放在源代码里换一组数据就要改代码重新编译实习报告里也讲不清楚。更合理的做法是把已知点、待定点、观测值拆成两份文本文件程序启动时读入。我习惯用如下格式已知点文件每行三列点名、X、Y观测文件每行五列观测类型、测站、后视、前视/目标、观测值。角度用“度.分秒”编码边长直接用米。这里要注意一旦定下格式就别中途改动否则后续数据处理脚本全部要跟着改。def read_known_points(path): known {} with open(path, r, encodingutf-8) as f: for line in f: line line.strip() if not line or line.startswith(#) or line.startswith(//): continue parts line.split(,) name, x, y parts[0].strip(), float(parts[1]), float(parts[2]) known[name] [x, y] return known def read_observations(path): obs [] with open(path, r, encodingutf-8) as f: for line in f: line line.strip() if not line or line.startswith(#) or line.startswith(//): continue parts line.split(,) otype parts[0].strip() if otype angle: station, back, fore, value parts[1].strip(), parts[2].strip(), parts[3].strip(), float(parts[4]) obs.append((angle, station, back, fore, value)) elif otype dist: station, target, value parts[1].strip(), parts[2].strip(), float(parts[3]) obs.append((dist, station, target, value)) return obs这段代码的作用是把文本文件转成内存里的字典和列表。known 字典的键是点名值是两个坐标分量obs 列表里每个元素是一个元组第一种是角度观测测站、后视、前视、观测值第二种是边长观测测站、目标、观测值。读文件时的#和//注释支持是为了方便在数据文件里做点位说明和分段标记排查问题时能快速定位是哪一段数据出错。3.2 待定点近似坐标的来源支导线推算与概略坐标文件间接平差是非线性问题线性化的过程必须在近似坐标处展开所以每个待定点都要有一个初始坐标。近似坐标从哪来实习里最常见的是支导线推算从一个已知点出发依次用角度和边长推算出每个待定点的坐标。推算过程不复杂先由已知边方位角和观测角推导线边方位角再用边长乘以方位角的正余弦得到坐标增量逐点累加。这里最怕的是方位角推算到负值或大于 360 度时不作归一化导致后续三角函数符号全部错乱。我一般会在每步方位角计算后做一次az az % (2 * math.pi)的归零处理。如果不想写支导线推算也可以直接在已知点文件里附一个待定点概略坐标文件用全站仪现场测的坐标或地图上量的坐标都行。近似坐标的精度要求不苛刻边长误差在几厘米内通常都能在迭代中收敛但极端情况下——比如待定点近似坐标偏离真实位置超过几米——第一轮迭代的改正数会巨大此时必须做迭代解算而不是只解一次这一点后面单独讲。3.3 点号到未知数序号的索引映射错误的高发地带平差程序里最大的坑不是数学推导而是点号和矩阵列号的对应关系。未知数向量 x 的长度是 2 乘以待定点个数第 k 个待定点的 dx 在 x[2k] 位置dy 在 x[2k1] 位置。如果观测文件里点名顺序和这个索引对不上误差方程的系数就填到了错误的列上法方程解出来是一堆没有意义的数。我建议在建观测系列之前先遍历所有观测值涉及的点名把待定点单独收集起来建立一个字典映射。def build_index_mapping(known, obs): idx {} for name in known: idx[name] -1 # 已知点不参与未知数求解索引记为-1 unknown_list [] for ob in obs: for name in ob[1:4] if ob[0] angle else ob[1:3]: if name not in known and name not in idx: idx[name] len(unknown_list) unknown_list.append(name) return idx, unknown_list这段逻辑是先给已知点赋 -1表示“该点在误差方程中系数为 0”再遍历观测值提取所有未出现在已知点字典里的点名按出现顺序分配未知数序号。注意这里遍历观测值的顺序就是索引分配顺序因此观测文件里点名的出现顺序会直接影响最终坐标解算结果的输出顺序。解决方法是平差结束后再按点名排序输出不要信任索引顺序就是输出顺序。4. 误差方程与法方程求解的代码实现从观测值到平差结果的一次跑通4.1 角度观测值的误差方程方位角导数的推导与代码角度观测值 β 是前视方向方位角减去后视方向方位角β T_sf - T_sb。平差时待定点坐标有微小改正数方位角 T 也会随之变化。对 T 做 Taylor 展开得到线性化的方向导数。先定义从测站 fr 到目标点 to 的方位角系数def direction_coeff(fr, to, coords): dx coords[to][0] - coords[fr][0] dy coords[to][1] - coords[fr][1] s math.hypot(dx, dy) az math.atan2(dy, dx) return -math.sin(az) / s, math.cos(az) / s返回值第一个是 dT/dx第二个是 dT/dy。为什么 dT/dx -sin(az)/s因为方位角 T atan2(dy, dx)对 dx 求偏导得到 -dy / s^2 -sin(T)/s对 dy 求偏导得到 dx / s^2 cos(T)/s。这两项的物理意义是目标点横向和纵向微小移动对方位角的影响量。角度误差方程的系数行组装如下def angle_equation(station, back, fore, beta_obs, coords, idx): a_sb, b_sb direction_coeff(station, back, coords) a_sf, b_sf direction_coeff(station, fore, coords) az_sb math.atan2(coords[back][1] - coords[station][1], coords[back][0] - coords[station][0]) az_sf math.atan2(coords[fore][1] - coords[station][1], coords[fore][0] - coords[station][0]) beta0 az_sf - az_sb l beta_obs - beta0 if l math.pi: l - 2 * math.pi if l -math.pi: l 2 * math.pi n_unk len(idx) row np.zeros(2 * n_unk) if idx[station] 0: row[2 * idx[station]] a_sb - a_sf row[2 * idx[station] 1] b_sb - b_sf if idx[back] 0: row[2 * idx[back]] - a_sb row[2 * idx[back] 1] - b_sb if idx[fore] 0: row[2 * idx[fore]] a_sf row[2 * idx[fore] 1] b_sf return row, l这里 l 是观测值减计算值方向是“观测值 - 计算值”对应误差方程 V A x - l 的符号约定。角度差可能超过 ±π所以要做弧度归一化否则方程组里的常数项出现 2π 级别的跳变解算出来的改正数会把坐标推出十万八千里。系数填列时测站点本身也是未知数时它的系数是两个方向导数之差逻辑推导正确之后写代码就一行赋值的事。4.2 边长观测值的误差方程线性化的陷阱边长观测值 S_ij 是两点间欧氏距离。设近似边长 S0线性化后得到V cos(T) * (dx_j - dx_i) sin(T) * (dy_j - dy_i) - (S0 - S_obs)其中 T 是近似方位角cos(T) 和 sin(T) 就是坐标增量除以边长。这个式子直观上比角度方程简单但初学者容易犯一个错把 cos(T) 和 sin(T) 的符号写反或者漏掉测站点的负号。def dist_equation(station, target, s_obs, coords, idx): dx coords[target][0] - coords[station][0] dy coords[target][1] - coords[station][1] s0 math.hypot(dx, dy) cos_t dx / s0 sin_t dy / s0 l s0 - s_obs n_unk len(idx) row np.zeros(2 * n_unk) if idx[station] 0: row[2 * idx[station]] - cos_t row[2 * idx[station] 1] - sin_t if idx[target] 0: row[2 * idx[target]] cos_t row[2 * idx[target] 1] sin_t return row, l这段代码的关键在于测站点的系数永远是负号目标点的系数永远是正号不管这两个点谁是已知点谁是待定点。idx 里已知点是 -1所以系数自动跳过。边长观测值单位是米常数项 l 也是米和角度观测值的弧度制混在一个法方程里时完全依靠权阵 P 来平衡量纲差异这里没有其他技巧。4.3 组法方程并解算N矩阵不要用求逆把所有观测值的行向量拼成矩阵 A常数项拼成向量 l再结合权阵 P法方程就是N A^T P AU A^T P lN x U写成 Python 就是几行矩阵运算n_obs len(obs) n_unk len(unknown_list) A np.zeros((n_obs, 2 * n_unk)) l np.zeros(n_obs) P np.zeros((n_obs, n_obs)) for k, ob in enumerate(obs): if ob[0] angle: A[k], l[k] angle_equation(ob[1], ob[2], ob[3], ob[4], coords, idx) P[k, k] 1.0 / (sigma_beta ** 2) else: A[k], l[k] dist_equation(ob[1], ob[2], ob[3], coords, idx) P[k, k] 1.0 / (sigma_dist ** 2) N A.T P A U A.T P l x np.linalg.solve(N, U)N 的维度是 2n_unk 阶方阵。实习网通常只有几个到十几个待定点矩阵阶数不超过几十np.linalg.solve 足够。这里有一个血泪经验不要用 np.linalg.inv(N) 求逆再乘 U虽然结果一样但数值稳定性差而且当 N 接近奇异时不会给出任何警告solve 至少会抛异常提醒你注意。实测数据里法方程的条件数经常很大用求逆容易把微小误差放大到坐标改正数上。4.4 精度评定单位权中误差与点位精度计算解出坐标改正数后误差方程残差 V A x - l。平差的标准做法是计算单位权中误差σ0 sqrt(V^T P V / (n - t))其中 n 是观测值总数t 是未知数个数。V^T P V 不必先算 V 再相乘可以用恒等式 V^T P V l^T P l - U^T x比先算残差再乘快很多而且避免了二次舍入误差。V A x - l vpv V P V sigma0 math.sqrt(vpv / (n_obs - n_unk)) Qxx np.linalg.inv(N) point_precision {} for name in unknown_list: k idx[name] q Qxx[2 * k, 2 * k] Qxx[2 * k 1, 2 * k 1] point_precision[name] sigma0 * math.sqrt(q)point_precision 里存的就是每个待定点的点位中误差单位和坐标一致米。实习报告里要求的点位误差精度表本质上就是把 Qxx 对角线元素开方再乘单位权中误差。注意 Qxx 是协因数阵不是协方差阵必须乘 σ0^2 才是方差前面代码直接对 Qxx 开方再乘 σ0顺序不能反过来。迭代问题也在这里处理间接平差只做一次线性化是不严格的改正数较大时必须把解出的 x 加到 coords 上重新组 A 和 l 再解直到 max(|x|) 小于某个阈值比如 1e-4 米。上面这段代码是单次解算版本实际使用时要包一个 5 到 10 次的循环判断是否收敛后再跳出。5. 导线网平差程序避坑指南5个让你的结果对不上的常见问题5.1 度分秒转换错误导致的“系统性偏差”现象平差结果和老师给的坐标相差固定的一小段距离而且所有点位偏移方向一致像是整个网被平移了几厘米。原因观测员记录的水平角不是十进制度而是“度.分秒”编码。如果直接把 245.30123 当成 245.30123 度传进 math.atan2 的导数计算方位角偏出去一大截误差方程的常数项全部带上一个系统偏差。解决写一个度分秒转十进制度数的小函数在读数时统一转换。注意编码规则是度占整数部分分占两位秒占两位秒可以带小数。def dms2deg(dms): d int(dms) minutes int((dms - d) * 100) seconds ((dms - d) * 100 - minutes) * 100 return d minutes / 60.0 seconds / 3600.0转换之后还要再乘 math.pi / 180 才是弧度这个步骤放进读观测文件的地方不要在误差方程里重复转否则代码里到处是 math.radians 会让排查变得很痛苦。5.2 错误定权把所有观测值当等权现象解算过程不报错坐标结果也在合理范围但精度评定的单位权中误差严重偏大点位误差大得离谱。原因权阵 P 被设成了单位阵也就是角度观测值和边长观测值被同等对待。水平角单位是弧度、边长单位是米数值上差了几个数量级法方程实际是在用非常畸形的比例加权等价于完全信任角度、忽略边长或反过来。解决按 2.3 节的办法给出测角中误差和测边中误差后构造对角权阵。实习中如果拿不准中误差就先用仪器标称精度角度取 1 秒换算成弧度约 4.85e-6边长取 2 毫米。这两者的权比自动会调整到合理的尺度精度评定结果也才有参考价值。5.3 漏掉测站点改正数误差方程少了一半系数现象平差结果的坐标变化量异常巨大每个点都偏移几十厘米但观测值残差 V 看起来又很正常。原因角度误差方程里测站点如果也是待定点它的坐标改正数同时影响后视和前视两个方位角算下来测站点上应有 (a_sb - a_sf, b_sb - b_sf) 的系数。写代码时只填了后视和前视点的系数漏掉测站点。这会让法方程错误地认为测站点完全可靠把误差全部往相邻点上推。解决回头检查 4.1 节 angle_equation 函数重点看测站点系数分支是否被 idx[station] 0 正确保护。如果测站点是已知点跳过是对的如果是待定点必须填入两方向导数之差。这条也是我在实习报告里反复强调的“误差方程少了一半系数”。5.4 改正数太大却不迭代线性化失效现象第一轮解算出来坐标改正数达到几十米甚至上百米直接加到近似坐标上第二轮迭代直接发散。原因间接平差的误差方程是线性化近似的只有在近似坐标接近真值时才能保证收敛半径。如果近似坐标是支导线推算但推算过程中有一站方位角传错或者起始数据本来就有误第一次解算的改正数会特别大线性化失效。解决在迭代循环里加一个判断如果 max(|x|) 5 米说明近似坐标可能有问题优先检查已知点和观测值是否读对而不是强行迭代。正常收敛的导线网第一次解算的坐标改正数应该在毫米到厘米量级超过 0.1 米就要回头查数据。迭代终止条件我习惯用 1e-4 米配合最多 10 次的循环上限既不漏收敛也不需要无意义空转。5.5 法方程奇异或病态网形不连接时的处理现象程序报 LinAlgError: Singular matrix或者解出来的坐标出现正负几十米的异常值没有任何合理解释。原因导线网里某个待定点没有被足够的观测值连接起来比如只有一个角度观测、没有边长观测导致法方程秩亏。更隐蔽的是两个未知数完全被同一条边锁定自由度为 0矩阵奇异。解决先检查观测数据覆盖情况确认每个待定点至少与两个其他点有观测联系。临时处理可以在该点上增加一个边长观测或者固定约束但这会改变网形结构必须在报告中注明。如果只是个别实测点位条件差可以考虑把该点从平差中剔除用支导线计算结果代替。要彻底理解这个问题还需要补一点秩亏自由网的知识但实习层面做到发现并定位问题就够了。6. 验证程序正确性与误差椭圆的小进阶6.1 用模拟数据验证平差程序的正确性调试平差程序最可靠的办法不是拿实习数据硬试而是先自己造一个数据网络。选两个已知点 A、B再预设几个待定点的“真坐标”用这些坐标正算出一组理论角度和理论边长当作观测值喂给程序。如果程序正确平差解出的坐标应当与预设真坐标完全一致误差等于零。在此基础上叠加微小随机误差到观测值上平差结果应该在一个合理的误差范围内贴近真值。def generate_test_network(): known {A: [0.0, 0.0], B: [500.0, 0.0]} true_coords { P1: [150.0, 120.0], P2: [350.0, 180.0], P3: [420.0, -50.0], } obs [] for station in [A, B, P1, P2, P3]: for target in [B, P1, P2, P3]: if station target: continue if station in known and target in known: continue dx true_coords[target][0] - true_coords[station][0] if station in true_coords else \ true_coords[target][0] - known[station][0] dy ... # 生成边长观测值和方位角再两两组合成角度观测 ... return known, true_coords, obs上面代码里省略了 dx、dy 的完整计算实际就是把站心坐标和目标坐标相减调用 math.hypot 拿距离调用 math.atan2 拿方位角再把两个方位角相减得到水平角观测值。这个模拟网的价值在于真值已知任何程序逻辑错误都会表现为坐标残差超限不需要依赖老师给的答案去“碰运气”。6.2 误差椭圆的进阶应用从Qxx里读什么平差程序输出坐标和点位中误差之后实习报告往往还要求画误差椭圆。误差椭圆的长半轴方向和大小完全由协因数阵 Qxx 对应点的 2×2 子矩阵决定。取待定点 P 对应的 Qxx 子矩阵它的两个特征值开方后乘单位权中误差就是椭圆长半轴和短半轴特征向量给出长轴方向。特征值大说明该点在一个方向上观测约束弱点位不确定性主要集中在这个方向。我实习时养成的习惯是每次平差完先做事后验算把观测值重新代回解算后的坐标算一遍残差看最大残差是否在观测精度的 3 倍范围以内。如果某个角度残差超过 5 秒大概率是观测值里有粗差或者程序中该行误差方程填错。直到了解这个习惯才感觉导线网平差程序真正从“交作业的代码”变成了“能拿来独立分析测量数据的工作”。希望帮到你。本文还有配套的精品资源点击获取
阅读完成 · 觉得有帮助？
