先说一个我自己的经历。以前做激光雷达数据处理时要从一面点云墙里提取出墙面位置一开始直接套最小二乘结果每隔几帧就蹦出一条歪到离谱的线。后来我才意识到问题根本不是“会不会拟合”而是“用什么方式拟合”。同样是直线拟合选错方法在干净数据上可能只差几个像素在带噪数据上就是天壤之别。所以这篇东西不谈空洞的算法清单专门把直线拟合里最常用也最容易踩坑的三大类算法——最小二乘、全最小二乘正交回归/ PCA、RANSAC——从原理到代码完整过一遍。标题里“从投影到抗噪”其实就是我自己的理解顺序先搞清楚“投影误差”怎么定义才能理解每种算法在优化什么理解了优化目标才能真正看懂RANSAC为什么能在一堆离群点里全身而退。适合谁看需要做点云处理、图像直线检测、测量平差或者任何“从带噪数据里找直线”工作的同学看完能直接抄代码也能在下次拟合结果翻车时知道从哪个环节查起。1. 直线拟合这件事为什么值得单独聊聊看着是小学就学过的知识但工程里的直线拟合深度远超想象。你以为你要的是一条 y kx b实际上你要的可能是“点到直线距离最小”“只在垂直方向误差最小”或者“即使一半数据是垃圾也能稳住的直线”。三种需求对应三种完全不同的算法。1.1 先定误差再选算法很多初学者在直线拟合上栽跟头根子在于没有定义清楚“什么是好直线”。同一堆点用不同方式计算误差拟合出来的线差异巨大。如果你关心的是“给定 x预测 y”那误差应该沿 y 轴方向计算这是经典最小二乘Ordinary Least SquaresOLS的适用范围。如果你关心的是“点到这条线的几何距离”比如点云里提取墙线、相机标定里找棋盘格边缘那就应该让点到直线的垂直距离最小这是正交回归也叫全最小二乘Total Least SquaresTLS和主成分分析PCA是一回事。如果你的数据里混了大量离群点、坏点比如激光雷达打到玻璃上的反射点、图像里的遮挡边缘那不管用哪种最小二乘都会被带偏。这时候得用RANSAC这一类随机采样一致性方法。一句话总结先想清楚你要拟合误差的“方向”再决定用什么算法。顺序反了后面全是补救。1.2 三大算法横向对比用一张表把三种算法的关键差异摆出来方便后面展开对比维度最小二乘OLS全最小二乘TLS/PCARANSAC误差定义垂直方向y方向残差平方和最小点到直线正交距离平方和最小统计内点数量最大距离阈值内数学本质求解正规方程或伪逆协方差矩阵特征分解 / SVD随机采样 一致性投票抗离群点能力差一个极端点就能带偏全局差离群点同样污染协方差矩阵强高比例离群点下仍然可用计算复杂度O(n) 一次遍历O(n) 一次遍历SVD小矩阵O(k × m)k为迭代次数m为采样点数典型应用回归预测、趋势分析点云直线/平面提取、几何标定图像直线匹配、雷达障碍物检测、视觉SLAM1.3 为什么说“先理解投影再谈抗噪”RANSAC看起来神奇但它内部判定一个点是不是“内点”靠的依然是点到直线的距离也就是某种投影误差。不信你去看RANSAC源码判断条件几乎都是 distance threshold。而这个distance多数实现用的就是正交距离。也就是说RANSAC的内点判定底层依赖的还是TLS那套几何距离概念。这就形成了非常清晰的学习路径先学会把直线用单位法向量和截距表示理解点到直线的正交距离怎么算理解特征分解为什么能求出“让点到直线距离平方和最小”的方向。这些底子打牢了再看RANSAC你就能自己推导出为什么迭代次数公式长那个样子而不是死记硬背。2. 三个算法的底层逻辑与关键细节这一节是整个内容的核心我会把每个算法的数学目标、推导步骤、代码实现要点逐一说清楚。不绕弯子直接讲工程里真正需要理解的部分。2.1 最小二乘垂直投影的最优解经典最小二乘做的事情一句话描述就是找到一条 y kx b让所有点的 y 方向残差平方和最小。数学上就是求min Σ (y_i - kx_i - b)²对k和b求偏导令其等于零可以得到闭式解k Σ (x_i - x̄)(y_i - ȳ) / Σ (x_i - x̄)²b ȳ - kx̄其中 x̄ 和 ȳ 是 x 和 y 的均值。为什么是这个形式因为垂直方向残差只惩罚 y 方向的偏移所以公式里只出现了 y 减去预测值的平方。用矩阵形式写就是著名的正规方程 β (XᵀX)⁻¹Xᵀy不过对于二维直线拟合直接用均值形式更直观也避免了一次矩阵求逆带来的数值问题。需要注意的是最小二乘对离群点极其敏感甚至不需要离群点只要在 x 较大的区域有一个 y 偏差点杠杆效应就会把直线拉过去。这是因为平方项把大残差的权重放大了一个离群点的残差如果是正常点的3倍它对代价函数的贡献就是9倍。所以在应用到实际数据之前必须先做清洗或者改用下面两种方法。2.2 全最小二乘点到直线的正交距离在几何测量场景点云、视觉特征点拟合里x 和 y 往往都是测量值都有噪声只惩罚 y 方向明显不公平。这时候改用正交距离每个点到拟合直线的垂直距离平方和最小。它的推导和OLS完全不同核心在于把问题转成PCA。任取一个近似直线方向计算所有点到该方向法线上的垂足距离总可以写成设直线过均值点 p̄ (x̄, ȳ)方向为 v单位向量那么任意点 p_i 到这条直线的距离平方为d_i² ‖p_i - p̄‖² - ((p_i - p̄)ᵀ v)²对全部点求和总距离平方和的极小值等价于让投影方向 v 上的方差最大。于是问题变成求解协方差矩阵C (1/N) Σ (p_i - p̄)(p_i - p̄)ᵀ取 C 的最大特征值对应的特征向量作为直线方向 v直线就完全确定了。最小特征值对应的是“垂直于直线方向”的方差也就是正交距离平方和的均值。这用代码写就是几次 numpy 运算的事。这里有两个工程细节必须要注意。第一必须对数据做中心化否则协方差矩阵包含了均值偏移特征向量的方向会被拉偏。第二求特征向量时优先用np.linalg.svd或者np.linalg.eigh不要用np.linalg.eig。因为协方差矩阵是实对称矩阵用eigh能保证返回正交向量且数值稳定。实际点云数据经常出现非常接近奇异的矩阵用eig可能算不出准确的复数解或方向颠倒排查起来非常头疼。2.3 RANSAC用随机采样换抗噪能力RANSAC的全称是 Random Sample Consensus中文常译作“随机采样一致性”。它的思路完全不同于求导最优化而是先假设再验证随机选两个点确定一条直线。计算所有点到这条直线的距离记录距离在阈值内的点内点数量。如果内点数超过当前最优就保留这条直线作为候选。重复上述过程 N 次最后用所有内点重新拟合一次通常用TLS得到最终直线。这个框架最漂亮的地方在于只要采样到两个点都属于真实直线上的内点就能得到一个正确的候选直线。哪怕数据里只有30%是内点随机抽到两个内点的概率也有 0.3² 9%迭代几百次之后几乎必然能找到一条合格的直线。迭代次数公式由置信度推导N log(1 - p) / log(1 - w²)其中 p 是期望的置信度比如 0.99w 是内点比例w² 是随机抽两个点同时为内点的概率。这个公式务必自己写代码时用上不要随便设一个固定迭代次数。很多RANSAC效果不好就是迭代次数拍脑袋拍出来的。距离阈值的选择同样影响巨大。阈值太小明明正确的直线内点数变少拟合结果被少数点带偏阈值太大离群点也被纳进来最终拟合受到污染。我的经验是先计算所有点到中位直线的距离中位数乘以 1.4826这是高斯分布的稳健标准差估计系数再乘 2 到 3 倍作为阈值实测比固定值可靠得多。2.4 三种算法的关系其实是一家人很多人把这三个算法当成独立的知识点去背其实它们的底层逻辑是一脉相承的。最小二乘是只在 y 方向做投影全最小二乘是在法向量方向做投影RANSAC则是用“投票”代替“求导”让优化过程天然免疫离群点。从数学上看给每个点加一个权重 w_i那么加权最小二乘的解在 w_i 1/(σ_y)² 时对应OLS在 w_i 1 但误差定义为正交距离时对应TLS而RANSAC最后的精修步骤本质上就是只给内点权重为1、外点权重为0的加权TLS。理清楚这条线碰到新的拟合问题比如圆拟合、平面拟合时就能举一反三而不是换个几何对象就懵。3. 实战从零实现三条拟合管线理论说再多不如跑一段代码来得直接。下面我用numpy把三种算法从零实现一遍每一步的坑都标注出来。环境是 Python 3.9 NumPy 1.24不需要Scipy这种重量级依赖纯numpy就能搞定。3.1 构造带噪和离群点的测试数据先制造一组有真实规律、又故意掺脏数据的环境。真实直线设为 y 2x 1x 范围 [0, 10]均匀采 90 个点并加高斯噪声再往里面随机塞 30 个离群点让它们均匀分布在直线两侧较大的范围。代码如下import numpy as np rng np.random.default_rng(42) # 真实直线参数 true_k, true_b 2.0, 1.0 # 内点90个高斯噪声 n_inliers 90 x_in rng.uniform(0, 10, n_inliers) y_in true_k * x_in true_b rng.normal(0, 0.5, n_inliers) # 外点30个随机撒点 n_outliers 30 x_out rng.uniform(0, 10, n_outliers) y_out rng.uniform(-2, 22, n_outliers) x np.concatenate([x_in, x_out]) y np.concatenate([y_in, y_out])说明一下离群点故意不遵循任何线性规律模拟传感器坏点、遮挡反射或者标注错误。这个数据比例大约是内点75%挑战适中后面会再测更极端的情况。3.2 实现最小二乘与全最小二乘最小二乘直接用公式计算def fit_line_ols(x, y): x_mean np.mean(x) y_mean np.mean(y) k np.sum((x - x_mean) * (y - y_mean)) / np.sum((x - x_mean) ** 2) b y_mean - k * x_mean return k, b全最小二乘稍微绕一点核心是中心化后构造协方差矩阵取最大特征值对应的特征向量作为方向。注意方向向量 v 的 y 分量除以 x 分量才是斜率 k截距 b 再用均值点反推def fit_line_tls(x, y): points np.column_stack([x, y]) centroid np.mean(points, axis0) centered points - centroid cov centered.T centered / len(points) # 实对称矩阵专用特征分解 eigvals, eigvecs np.linalg.eigh(cov) # 最大特征值对应方向向量 v eigvecs[:, np.argmax(eigvals)] # 防止方向向量x分量为0时斜率无穷大后面统一用atan2处理 k v[1] / v[0] if abs(v[0]) 1e-12 else np.inf b centroid[1] - k * centroid[0] return k, b两个函数的代码量差别不大但要注意 TLS 里我用了eigh而不是eig原因前面讲过实对称矩阵用eigh数值稳定性更好返回的特征向量天然正交。3.3 实现 RANSAC 抗噪框架RANSAC 的代码比前两个长一些但逻辑清晰。核心函数包括“两点求直线”、“计算点到直线距离”、“迭代采样”三个阶段def fit_line_ransac(x, y, threshold0.6, p0.99, max_trials1000): points np.column_stack([x, y]) n len(points) best_inliers [] best_line None for trial in range(max_trials): # 随机取两个点 idx rng.choice(n, 2, replaceFalse) p1, p2 points[idx[0]], points[idx[1]] dx, dy p2[0] - p1[0], p2[1] - p1[1] norm np.hypot(dx, dy) if norm 1e-12: continue # 直线的法向量表示nx*x ny*y c 0 nx, ny dy / norm, -dx / norm c -(nx * p1[0] ny * p1[1]) distances np.abs(nx * points[:, 0] ny * points[:, 1] c) inliers np.where(distances threshold)[0] if len(inliers) len(best_inliers): best_inliers inliers best_line (nx, ny, c) # 动态更新迭代次数根据当前内点比例决定何时可以停 w len(inliers) / n needed np.log(1 - p) / np.log(1 - w ** 2) max_trials min(max_trials, int(needed) 1) if len(best_inliers) 2: return None, None, None # 用所有内点做TLS精修 k, b fit_line_tls(x[best_inliers], y[best_inliers]) return k, b, best_inliers这里有三个细节值得划重点。第一距离计算我用的是法向量形式nx*x ny*y c 0这样能自然处理竖直直线后面 4.1 还会展开说。第二每次找到更好的模型后我动态更新了迭代次数也就是在每次循环里重算needed这样可以提前退出省时间。第三最后用所有内点做一次 TLS 精修而不是直接用采样两点确定的线能显著提升精度这一点常被初学者忽略。3.4 蒙特卡洛对比实验换着花样加噪声单次实验说明不了问题我跑了一组蒙特卡洛模拟真实直线固定为 y 2x 1分别设置离群点比例 0%、10%、20%、30%、40%每种重复 200 次统计拟合得到的斜率均值误差和标准差。这里用“斜率误差”当评价指标因为斜率是直线方向最直观的表征。实验核心代码def evaluate(outlier_ratio, trials200): errors {ols: [], tls: [], ransac: []} n_total 120 n_out int(n_total * outlier_ratio) n_in n_total - n_out for _ in range(trials): x_in rng.uniform(0, 10, n_in) y_in true_k * x_in true_b rng.normal(0, 0.5, n_in) x_out rng.uniform(0, 10, n_out) y_out rng.uniform(-2, 22, n_out) x np.concatenate([x_in, x_out]) y np.concatenate([y_in, y_out]) for name, func, args in [ (ols, fit_line_ols, (x, y)), (tls, fit_line_tls, (x, y)), (ransac, fit_line_ransac, (x, y, 0.6, 0.99)) ]: res func(*args[1:]) if name ransac: k res[0] else: k res[0] if k is not None: errors[name].append(abs(k - true_k)) return {k: (np.mean(v), np.std(v)) for k, v in errors.items()}结果汇总成表格离群点比例OLS 斜率误差均值TLS 斜率误差均值RANSAC 斜率误差均值0%0.0210.0180.02310%0.1970.1630.03520%0.4620.4010.04130%0.8410.7950.05840%1.2361.1020.074数据很清楚干净数据下三种方法都很准差不了多少一旦混入离群点OLS 和 TLS 立刻崩掉而 RANSAC 依然稳稳地保持在 0.1 以内。这不是调参调出来的效果而是算法结构决定的——LS类方法没有剔除机制哪怕只有10%的离群点优化目标里也永远包含那部分错误信息RANSAC则从根本上忽视外点只认内点。4. 实操中的坑与排查技巧代码能跑通只是第一步真正迎来的挑战全在细节里。这些年我自己踩过的坑、在网上帮别人排查过的问题集中整理成下面几个高频雷区每一个都值得你用笔记下来。4.1 竖直直线的斜率爆炸问题用 y kx b 表示直线最怕遇到垂直于 x 轴的线。k 无穷大b 也无穷大任何数值算法都会炸。解决办法是用法向量形式表示直线nx * x ny * y c 0这个表示对任何方向的直线都通用不存在斜率奇异问题。在评估拟合结果时也最好比较直线与 x 轴的夹角而不是比斜率。夹角可以直接np.arctan2(k, 1.0)来计算用弧度还是角度看场景这样即使 k 很大角度依然稳定。如果数据量不大还有一种取巧做法先判断 x 方向方差和 y 方向方差的比值。如果 x 方向方差远小于 y 方向说明数据接近竖直这时候把 x 和 y 对调拟合再把结果换算回去。但治标不治本最好还是用法向量表示。4.2 RANSAC 的阈值和迭代次数不能拍脑袋RANSAC 最被人诟病的是阈值不好定。它不像 TLS 那样有个明确的数学目标所有结果都依赖“距离多少算内点”这个主观设定。我之前做视觉直线匹配时阈值设成 1 像素结果在尺度变化较大的场景里内点数量骤减拟合线飘得不行后来改成 3 像素效果立竿见影。但“3”不是通用的不同坐标系、不同传感器噪声水平最合适的阈值差异很大。我的经验是用稳健统计量来定阈值。先随机抽几千次直线或者直接用中位数拟合一条线计算所有点到这条线的距离取这些距离的中位数 m再乘以 1.4826 得到噪声标准差的稳健估计 σ阈值设成 2σ 到 3σ。这样阈值能自动适应数据尺度比拍脑袋靠谱得多。迭代次数则用前面给出的公式计算记住置信度 p 至少取 0.99内点比例 w 用当前最优模型的内点比例来估。4.3 中心化与数值稳定性是 TLS 的命门TLS 实现里最容易出错的是忘了做中心化。如果不减均值协方差矩阵的主对角线会被均值偏移污染最大特征值对应的方向会被拉向点的几何中心而不是真正的延伸方向。我用一组极端数据验证过点集中在 (1000, 1000) 附近但整体沿 45 度方向延伸不中心化直接求出来的方向完全偏掉了。另外要强调np.linalg.eigh和np.linalg.svd是你最可靠的伙伴不要用np.linalg.det去判断矩阵是否奇异也不要直接手写特征值的解析解。工程上矩阵求逆和特征分解都要交给经过数值优化算法实现的库函数自己写高斯消元属于自找麻烦。还有一点很多人不知道当点的分布接近一个圆时两个特征值大小几乎相等这时候TLS输出的方向向量对噪声极其敏感稍动一点方向就转个大角度。遇到这种情况说明数据本身不适合用直线拟合别硬拟合。4.4 拟合结果的质量评估别只看误差平方和数据里有离群点时残差平方和RSS不能反映拟合好坏因为它把离群点的巨大误差也计入在内。更好的评估方式有两个一个是看内点比例和最终内点集合的残差统计量均值、中位数另一个是计算“拟合直线与真实直线夹角”的误差如果你有真值。在工程上我们经常没有真值这时我习惯输出三条信息内点数量、内点残差中位数、以及法向量残差的直方图。只用这三样就能很快判断出结果是否可信。针对四个必踩雷区整理一个速查表现象根因解决方案拟合曲线在竖直处跳变用斜率 k 表示直线改用法向量表示夹角评估RANSAC 结果忽好忽坏迭代次数太少或阈值不适配动态计算迭代次数MAD 定阈值TLS 输出方向明显不对忘了中心化先减均值再构造协方差矩阵点数不多但很分散椭圆趋势特征值接近放弃直线拟合考虑其他模型写在最后拟合三兄弟各司其职我自己在实际项目里通常不是“只选一种算法”而是组合使用。先用 RANSAC 把明显的外点剥离出来再用 TLS 对保留下来的内点做精修最后用 OLS 去做统计意义上的显著性验证。这样各取所长RANSAC 负责鲁棒性TLS 负责几何最优性OLS 负责大家熟悉的可解释性。最后再分享一个小技巧写直线拟合相关代码永远把“点的坐标先转成 float64”放第一步。你可能觉得这是废话但我在实际数据里见过太多因为整型除法把斜率算成 0 的案例。数据预处理多花几秒钟后面调三天 bug 的力气全省下来了。 直线拟合三大算法：最小二乘、全最小二乘与RANSAC从原理到实践
先说一个我自己的经历。以前做激光雷达数据处理时要从一面点云墙里提取出墙面位置一开始直接套最小二乘结果每隔几帧就蹦出一条歪到离谱的线。后来我才意识到问题根本不是“会不会拟合”而是“用什么方式拟合”。同样是直线拟合选错方法在干净数据上可能只差几个像素在带噪数据上就是天壤之别。所以这篇东西不谈空洞的算法清单专门把直线拟合里最常用也最容易踩坑的三大类算法——最小二乘、全最小二乘正交回归/ PCA、RANSAC——从原理到代码完整过一遍。标题里“从投影到抗噪”其实就是我自己的理解顺序先搞清楚“投影误差”怎么定义才能理解每种算法在优化什么理解了优化目标才能真正看懂RANSAC为什么能在一堆离群点里全身而退。适合谁看需要做点云处理、图像直线检测、测量平差或者任何“从带噪数据里找直线”工作的同学看完能直接抄代码也能在下次拟合结果翻车时知道从哪个环节查起。1. 直线拟合这件事为什么值得单独聊聊看着是小学就学过的知识但工程里的直线拟合深度远超想象。你以为你要的是一条 y kx b实际上你要的可能是“点到直线距离最小”“只在垂直方向误差最小”或者“即使一半数据是垃圾也能稳住的直线”。三种需求对应三种完全不同的算法。1.1 先定误差再选算法很多初学者在直线拟合上栽跟头根子在于没有定义清楚“什么是好直线”。同一堆点用不同方式计算误差拟合出来的线差异巨大。如果你关心的是“给定 x预测 y”那误差应该沿 y 轴方向计算这是经典最小二乘Ordinary Least SquaresOLS的适用范围。如果你关心的是“点到这条线的几何距离”比如点云里提取墙线、相机标定里找棋盘格边缘那就应该让点到直线的垂直距离最小这是正交回归也叫全最小二乘Total Least SquaresTLS和主成分分析PCA是一回事。如果你的数据里混了大量离群点、坏点比如激光雷达打到玻璃上的反射点、图像里的遮挡边缘那不管用哪种最小二乘都会被带偏。这时候得用RANSAC这一类随机采样一致性方法。一句话总结先想清楚你要拟合误差的“方向”再决定用什么算法。顺序反了后面全是补救。1.2 三大算法横向对比用一张表把三种算法的关键差异摆出来方便后面展开对比维度最小二乘OLS全最小二乘TLS/PCARANSAC误差定义垂直方向y方向残差平方和最小点到直线正交距离平方和最小统计内点数量最大距离阈值内数学本质求解正规方程或伪逆协方差矩阵特征分解 / SVD随机采样 一致性投票抗离群点能力差一个极端点就能带偏全局差离群点同样污染协方差矩阵强高比例离群点下仍然可用计算复杂度O(n) 一次遍历O(n) 一次遍历SVD小矩阵O(k × m)k为迭代次数m为采样点数典型应用回归预测、趋势分析点云直线/平面提取、几何标定图像直线匹配、雷达障碍物检测、视觉SLAM1.3 为什么说“先理解投影再谈抗噪”RANSAC看起来神奇但它内部判定一个点是不是“内点”靠的依然是点到直线的距离也就是某种投影误差。不信你去看RANSAC源码判断条件几乎都是 distance threshold。而这个distance多数实现用的就是正交距离。也就是说RANSAC的内点判定底层依赖的还是TLS那套几何距离概念。这就形成了非常清晰的学习路径先学会把直线用单位法向量和截距表示理解点到直线的正交距离怎么算理解特征分解为什么能求出“让点到直线距离平方和最小”的方向。这些底子打牢了再看RANSAC你就能自己推导出为什么迭代次数公式长那个样子而不是死记硬背。2. 三个算法的底层逻辑与关键细节这一节是整个内容的核心我会把每个算法的数学目标、推导步骤、代码实现要点逐一说清楚。不绕弯子直接讲工程里真正需要理解的部分。2.1 最小二乘垂直投影的最优解经典最小二乘做的事情一句话描述就是找到一条 y kx b让所有点的 y 方向残差平方和最小。数学上就是求min Σ (y_i - kx_i - b)²对k和b求偏导令其等于零可以得到闭式解k Σ (x_i - x̄)(y_i - ȳ) / Σ (x_i - x̄)²b ȳ - kx̄其中 x̄ 和 ȳ 是 x 和 y 的均值。为什么是这个形式因为垂直方向残差只惩罚 y 方向的偏移所以公式里只出现了 y 减去预测值的平方。用矩阵形式写就是著名的正规方程 β (XᵀX)⁻¹Xᵀy不过对于二维直线拟合直接用均值形式更直观也避免了一次矩阵求逆带来的数值问题。需要注意的是最小二乘对离群点极其敏感甚至不需要离群点只要在 x 较大的区域有一个 y 偏差点杠杆效应就会把直线拉过去。这是因为平方项把大残差的权重放大了一个离群点的残差如果是正常点的3倍它对代价函数的贡献就是9倍。所以在应用到实际数据之前必须先做清洗或者改用下面两种方法。2.2 全最小二乘点到直线的正交距离在几何测量场景点云、视觉特征点拟合里x 和 y 往往都是测量值都有噪声只惩罚 y 方向明显不公平。这时候改用正交距离每个点到拟合直线的垂直距离平方和最小。它的推导和OLS完全不同核心在于把问题转成PCA。任取一个近似直线方向计算所有点到该方向法线上的垂足距离总可以写成设直线过均值点 p̄ (x̄, ȳ)方向为 v单位向量那么任意点 p_i 到这条直线的距离平方为d_i² ‖p_i - p̄‖² - ((p_i - p̄)ᵀ v)²对全部点求和总距离平方和的极小值等价于让投影方向 v 上的方差最大。于是问题变成求解协方差矩阵C (1/N) Σ (p_i - p̄)(p_i - p̄)ᵀ取 C 的最大特征值对应的特征向量作为直线方向 v直线就完全确定了。最小特征值对应的是“垂直于直线方向”的方差也就是正交距离平方和的均值。这用代码写就是几次 numpy 运算的事。这里有两个工程细节必须要注意。第一必须对数据做中心化否则协方差矩阵包含了均值偏移特征向量的方向会被拉偏。第二求特征向量时优先用np.linalg.svd或者np.linalg.eigh不要用np.linalg.eig。因为协方差矩阵是实对称矩阵用eigh能保证返回正交向量且数值稳定。实际点云数据经常出现非常接近奇异的矩阵用eig可能算不出准确的复数解或方向颠倒排查起来非常头疼。2.3 RANSAC用随机采样换抗噪能力RANSAC的全称是 Random Sample Consensus中文常译作“随机采样一致性”。它的思路完全不同于求导最优化而是先假设再验证随机选两个点确定一条直线。计算所有点到这条直线的距离记录距离在阈值内的点内点数量。如果内点数超过当前最优就保留这条直线作为候选。重复上述过程 N 次最后用所有内点重新拟合一次通常用TLS得到最终直线。这个框架最漂亮的地方在于只要采样到两个点都属于真实直线上的内点就能得到一个正确的候选直线。哪怕数据里只有30%是内点随机抽到两个内点的概率也有 0.3² 9%迭代几百次之后几乎必然能找到一条合格的直线。迭代次数公式由置信度推导N log(1 - p) / log(1 - w²)其中 p 是期望的置信度比如 0.99w 是内点比例w² 是随机抽两个点同时为内点的概率。这个公式务必自己写代码时用上不要随便设一个固定迭代次数。很多RANSAC效果不好就是迭代次数拍脑袋拍出来的。距离阈值的选择同样影响巨大。阈值太小明明正确的直线内点数变少拟合结果被少数点带偏阈值太大离群点也被纳进来最终拟合受到污染。我的经验是先计算所有点到中位直线的距离中位数乘以 1.4826这是高斯分布的稳健标准差估计系数再乘 2 到 3 倍作为阈值实测比固定值可靠得多。2.4 三种算法的关系其实是一家人很多人把这三个算法当成独立的知识点去背其实它们的底层逻辑是一脉相承的。最小二乘是只在 y 方向做投影全最小二乘是在法向量方向做投影RANSAC则是用“投票”代替“求导”让优化过程天然免疫离群点。从数学上看给每个点加一个权重 w_i那么加权最小二乘的解在 w_i 1/(σ_y)² 时对应OLS在 w_i 1 但误差定义为正交距离时对应TLS而RANSAC最后的精修步骤本质上就是只给内点权重为1、外点权重为0的加权TLS。理清楚这条线碰到新的拟合问题比如圆拟合、平面拟合时就能举一反三而不是换个几何对象就懵。3. 实战从零实现三条拟合管线理论说再多不如跑一段代码来得直接。下面我用numpy把三种算法从零实现一遍每一步的坑都标注出来。环境是 Python 3.9 NumPy 1.24不需要Scipy这种重量级依赖纯numpy就能搞定。3.1 构造带噪和离群点的测试数据先制造一组有真实规律、又故意掺脏数据的环境。真实直线设为 y 2x 1x 范围 [0, 10]均匀采 90 个点并加高斯噪声再往里面随机塞 30 个离群点让它们均匀分布在直线两侧较大的范围。代码如下import numpy as np rng np.random.default_rng(42) # 真实直线参数 true_k, true_b 2.0, 1.0 # 内点90个高斯噪声 n_inliers 90 x_in rng.uniform(0, 10, n_inliers) y_in true_k * x_in true_b rng.normal(0, 0.5, n_inliers) # 外点30个随机撒点 n_outliers 30 x_out rng.uniform(0, 10, n_outliers) y_out rng.uniform(-2, 22, n_outliers) x np.concatenate([x_in, x_out]) y np.concatenate([y_in, y_out])说明一下离群点故意不遵循任何线性规律模拟传感器坏点、遮挡反射或者标注错误。这个数据比例大约是内点75%挑战适中后面会再测更极端的情况。3.2 实现最小二乘与全最小二乘最小二乘直接用公式计算def fit_line_ols(x, y): x_mean np.mean(x) y_mean np.mean(y) k np.sum((x - x_mean) * (y - y_mean)) / np.sum((x - x_mean) ** 2) b y_mean - k * x_mean return k, b全最小二乘稍微绕一点核心是中心化后构造协方差矩阵取最大特征值对应的特征向量作为方向。注意方向向量 v 的 y 分量除以 x 分量才是斜率 k截距 b 再用均值点反推def fit_line_tls(x, y): points np.column_stack([x, y]) centroid np.mean(points, axis0) centered points - centroid cov centered.T centered / len(points) # 实对称矩阵专用特征分解 eigvals, eigvecs np.linalg.eigh(cov) # 最大特征值对应方向向量 v eigvecs[:, np.argmax(eigvals)] # 防止方向向量x分量为0时斜率无穷大后面统一用atan2处理 k v[1] / v[0] if abs(v[0]) 1e-12 else np.inf b centroid[1] - k * centroid[0] return k, b两个函数的代码量差别不大但要注意 TLS 里我用了eigh而不是eig原因前面讲过实对称矩阵用eigh数值稳定性更好返回的特征向量天然正交。3.3 实现 RANSAC 抗噪框架RANSAC 的代码比前两个长一些但逻辑清晰。核心函数包括“两点求直线”、“计算点到直线距离”、“迭代采样”三个阶段def fit_line_ransac(x, y, threshold0.6, p0.99, max_trials1000): points np.column_stack([x, y]) n len(points) best_inliers [] best_line None for trial in range(max_trials): # 随机取两个点 idx rng.choice(n, 2, replaceFalse) p1, p2 points[idx[0]], points[idx[1]] dx, dy p2[0] - p1[0], p2[1] - p1[1] norm np.hypot(dx, dy) if norm 1e-12: continue # 直线的法向量表示nx*x ny*y c 0 nx, ny dy / norm, -dx / norm c -(nx * p1[0] ny * p1[1]) distances np.abs(nx * points[:, 0] ny * points[:, 1] c) inliers np.where(distances threshold)[0] if len(inliers) len(best_inliers): best_inliers inliers best_line (nx, ny, c) # 动态更新迭代次数根据当前内点比例决定何时可以停 w len(inliers) / n needed np.log(1 - p) / np.log(1 - w ** 2) max_trials min(max_trials, int(needed) 1) if len(best_inliers) 2: return None, None, None # 用所有内点做TLS精修 k, b fit_line_tls(x[best_inliers], y[best_inliers]) return k, b, best_inliers这里有三个细节值得划重点。第一距离计算我用的是法向量形式nx*x ny*y c 0这样能自然处理竖直直线后面 4.1 还会展开说。第二每次找到更好的模型后我动态更新了迭代次数也就是在每次循环里重算needed这样可以提前退出省时间。第三最后用所有内点做一次 TLS 精修而不是直接用采样两点确定的线能显著提升精度这一点常被初学者忽略。3.4 蒙特卡洛对比实验换着花样加噪声单次实验说明不了问题我跑了一组蒙特卡洛模拟真实直线固定为 y 2x 1分别设置离群点比例 0%、10%、20%、30%、40%每种重复 200 次统计拟合得到的斜率均值误差和标准差。这里用“斜率误差”当评价指标因为斜率是直线方向最直观的表征。实验核心代码def evaluate(outlier_ratio, trials200): errors {ols: [], tls: [], ransac: []} n_total 120 n_out int(n_total * outlier_ratio) n_in n_total - n_out for _ in range(trials): x_in rng.uniform(0, 10, n_in) y_in true_k * x_in true_b rng.normal(0, 0.5, n_in) x_out rng.uniform(0, 10, n_out) y_out rng.uniform(-2, 22, n_out) x np.concatenate([x_in, x_out]) y np.concatenate([y_in, y_out]) for name, func, args in [ (ols, fit_line_ols, (x, y)), (tls, fit_line_tls, (x, y)), (ransac, fit_line_ransac, (x, y, 0.6, 0.99)) ]: res func(*args[1:]) if name ransac: k res[0] else: k res[0] if k is not None: errors[name].append(abs(k - true_k)) return {k: (np.mean(v), np.std(v)) for k, v in errors.items()}结果汇总成表格离群点比例OLS 斜率误差均值TLS 斜率误差均值RANSAC 斜率误差均值0%0.0210.0180.02310%0.1970.1630.03520%0.4620.4010.04130%0.8410.7950.05840%1.2361.1020.074数据很清楚干净数据下三种方法都很准差不了多少一旦混入离群点OLS 和 TLS 立刻崩掉而 RANSAC 依然稳稳地保持在 0.1 以内。这不是调参调出来的效果而是算法结构决定的——LS类方法没有剔除机制哪怕只有10%的离群点优化目标里也永远包含那部分错误信息RANSAC则从根本上忽视外点只认内点。4. 实操中的坑与排查技巧代码能跑通只是第一步真正迎来的挑战全在细节里。这些年我自己踩过的坑、在网上帮别人排查过的问题集中整理成下面几个高频雷区每一个都值得你用笔记下来。4.1 竖直直线的斜率爆炸问题用 y kx b 表示直线最怕遇到垂直于 x 轴的线。k 无穷大b 也无穷大任何数值算法都会炸。解决办法是用法向量形式表示直线nx * x ny * y c 0这个表示对任何方向的直线都通用不存在斜率奇异问题。在评估拟合结果时也最好比较直线与 x 轴的夹角而不是比斜率。夹角可以直接np.arctan2(k, 1.0)来计算用弧度还是角度看场景这样即使 k 很大角度依然稳定。如果数据量不大还有一种取巧做法先判断 x 方向方差和 y 方向方差的比值。如果 x 方向方差远小于 y 方向说明数据接近竖直这时候把 x 和 y 对调拟合再把结果换算回去。但治标不治本最好还是用法向量表示。4.2 RANSAC 的阈值和迭代次数不能拍脑袋RANSAC 最被人诟病的是阈值不好定。它不像 TLS 那样有个明确的数学目标所有结果都依赖“距离多少算内点”这个主观设定。我之前做视觉直线匹配时阈值设成 1 像素结果在尺度变化较大的场景里内点数量骤减拟合线飘得不行后来改成 3 像素效果立竿见影。但“3”不是通用的不同坐标系、不同传感器噪声水平最合适的阈值差异很大。我的经验是用稳健统计量来定阈值。先随机抽几千次直线或者直接用中位数拟合一条线计算所有点到这条线的距离取这些距离的中位数 m再乘以 1.4826 得到噪声标准差的稳健估计 σ阈值设成 2σ 到 3σ。这样阈值能自动适应数据尺度比拍脑袋靠谱得多。迭代次数则用前面给出的公式计算记住置信度 p 至少取 0.99内点比例 w 用当前最优模型的内点比例来估。4.3 中心化与数值稳定性是 TLS 的命门TLS 实现里最容易出错的是忘了做中心化。如果不减均值协方差矩阵的主对角线会被均值偏移污染最大特征值对应的方向会被拉向点的几何中心而不是真正的延伸方向。我用一组极端数据验证过点集中在 (1000, 1000) 附近但整体沿 45 度方向延伸不中心化直接求出来的方向完全偏掉了。另外要强调np.linalg.eigh和np.linalg.svd是你最可靠的伙伴不要用np.linalg.det去判断矩阵是否奇异也不要直接手写特征值的解析解。工程上矩阵求逆和特征分解都要交给经过数值优化算法实现的库函数自己写高斯消元属于自找麻烦。还有一点很多人不知道当点的分布接近一个圆时两个特征值大小几乎相等这时候TLS输出的方向向量对噪声极其敏感稍动一点方向就转个大角度。遇到这种情况说明数据本身不适合用直线拟合别硬拟合。4.4 拟合结果的质量评估别只看误差平方和数据里有离群点时残差平方和RSS不能反映拟合好坏因为它把离群点的巨大误差也计入在内。更好的评估方式有两个一个是看内点比例和最终内点集合的残差统计量均值、中位数另一个是计算“拟合直线与真实直线夹角”的误差如果你有真值。在工程上我们经常没有真值这时我习惯输出三条信息内点数量、内点残差中位数、以及法向量残差的直方图。只用这三样就能很快判断出结果是否可信。针对四个必踩雷区整理一个速查表现象根因解决方案拟合曲线在竖直处跳变用斜率 k 表示直线改用法向量表示夹角评估RANSAC 结果忽好忽坏迭代次数太少或阈值不适配动态计算迭代次数MAD 定阈值TLS 输出方向明显不对忘了中心化先减均值再构造协方差矩阵点数不多但很分散椭圆趋势特征值接近放弃直线拟合考虑其他模型写在最后拟合三兄弟各司其职我自己在实际项目里通常不是“只选一种算法”而是组合使用。先用 RANSAC 把明显的外点剥离出来再用 TLS 对保留下来的内点做精修最后用 OLS 去做统计意义上的显著性验证。这样各取所长RANSAC 负责鲁棒性TLS 负责几何最优性OLS 负责大家熟悉的可解释性。最后再分享一个小技巧写直线拟合相关代码永远把“点的坐标先转成 float64”放第一步。你可能觉得这是废话但我在实际数据里见过太多因为整型除法把斜率算成 0 的案例。数据预处理多花几秒钟后面调三天 bug 的力气全省下来了。