直线拟合 - 最小二乘法与 RANSAC 算法

课程视频链接：北京邮电大学_计算机视觉_鲁鹏_第三次课_拟合 和 北京邮电大学_计算机视觉_鲁鹏_第4次课_拟合（RANSAC复习，hough）。

1. 问题描述

如下图所示，设二维平面内存在 $n$ 个点 { ( x i y i ) } i = 1 n \{ \begin{pmatrix} x_i \\ y_i \end{pmatrix} \}_{i = 1}^{n} {(xi​yi​​)}i=1n​，现需构造最优直线模型来拟合该数据集，下面从三种不同的方法来解决这个问题。

对下面所述的各个算法做一个总结图，如下图所示：

2. 最小二乘法

2.1. 传统最小二乘法

最小二乘法(Least Squares, LS)构造最优直线模型 $\hat{y}=mx+b$ 来拟合二维平面内的点，其目标是最小化所有观测点的残差平方和 ${\sum }_{i=1}^n(y_i-{\hat{y}}_i{)}^2$，其中残差 $e_i=y_i-{\hat{y}}_i$ 表示实际值 $y_i$ 与模型预测值 $\widehat{y_i}$ 的垂直偏差，如下图所示：
设 $E={\sum }_{i=1}^n(y_i-m{x}_i-b{)}^2$，$\mathbf{Y}=\left(\begin{array}{c}{y}_1\\ ⋮\\ y_n\end{array}\right)$，$\mathbf{X}=\left(\begin{array}{cc}{x}_1& 1\\ ⋮& ⋮\\ x_n& 1\end{array}\right)$，$\mathbf{B}=\left(\begin{array}{c}m\\ b\end{array}\right)$，则有：
$$E=\Vert \mathbf{Y}-\mathbf{XB}{\Vert }^2\text{\hspace{1em}(1)}$$现在我们需要求解 $\mathbf{B}$ 取什么值时，$E$ 取得最小值。对 $E$ 进行展开有：
$$\begin{array}{rl}E& =(\mathbf{Y}-\mathbf{XB}{)}^T(\mathbf{Y}-\mathbf{XB})=({\mathbf{Y}}^T-{\mathbf{B}}^T{\mathbf{X}}^T)(\mathbf{Y}-\mathbf{XB})\\ & ={\mathbf{Y}}^T\mathbf{Y}-{\mathbf{Y}}^T\mathbf{XB}-{\mathbf{B}}^T{\mathbf{X}}^T\mathbf{Y}+{\mathbf{B}}^T{\mathbf{X}}^T\mathbf{XB}\end{array}$$我们需要计算出 $\frac{\partial E}{\partial \mathbf{B}}$，然后根据 $\frac{\partial E}{\partial \mathbf{B}}=0$ 解出 $\mathbf{B}$。
参考博客：机器学习中的数学理论1：三步搞定矩阵求导，求导步骤如下：

1. 求微分 $dE$。
   $$\begin{array}{rl}dE& =-{\mathbf{Y}}^T\mathbf{X}d\mathbf{B}-d({\mathbf{B}}^T){\mathbf{X}}^T\mathbf{Y}+d({\mathbf{B}}^T){\mathbf{X}}^T\mathbf{XB}+{\mathbf{B}}^T{\mathbf{X}}^T\mathbf{X}d\mathbf{B}\\ & =({\mathbf{B}}^T{\mathbf{X}}^T\mathbf{X}-{\mathbf{Y}}^T\mathbf{X})d\mathbf{B}+(d\mathbf{B}{)}^T({\mathbf{X}}^T\mathbf{XB}-{\mathbf{X}}^T\mathbf{Y})\end{array}$$
2. 计算 $tr(dE)$。
   $$\begin{array}{rl}dE& =tr(dE)=tr[({\mathbf{B}}^T{\mathbf{X}}^T\mathbf{X}-{\mathbf{Y}}^T\mathbf{X})d\mathbf{B}+(d\mathbf{B}{)}^T({\mathbf{X}}^T\mathbf{XB}-{\mathbf{X}}^T\mathbf{Y})]\\ & =tr[{\mathbf{B}}^T{\mathbf{X}}^T\mathbf{X}-{\mathbf{Y}}^T\mathbf{X})d\mathbf{B}]+tr[(d\mathbf{B}{)}^T({\mathbf{X}}^T\mathbf{XB}-{\mathbf{X}}^T\mathbf{Y})]\\ & =tr[({\mathbf{B}}^T{\mathbf{X}}^T\mathbf{X}-{\mathbf{Y}}^T\mathbf{X})d\mathbf{B}]+tr[({\mathbf{X}}^T\mathbf{XB}-{\mathbf{X}}^T\mathbf{Y}{)}^{T}d\mathbf{B}]\\ & =tr[2({\mathbf{B}}^T{\mathbf{X}}^T\mathbf{X}-{\mathbf{Y}}^T\mathbf{X})d\mathbf{B}]\end{array}$$
3. 根据 $dE=tr({\frac{\partial E}{\partial \mathbf{B}}}^{T}d\mathbf{B})$，求 $\frac{\partial E}{\partial \mathbf{B}}$。由上述推导可知，
   $${\frac{\partial E}{\partial \mathbf{B}}}^T=2({\mathbf{B}}^T{\mathbf{X}}^T\mathbf{X}-{\mathbf{Y}}^T\mathbf{X})\Rightarrow \frac{\partial E}{\partial \mathbf{B}}=2({\mathbf{X}}^T\mathbf{XB}-{\mathbf{X}}^T\mathbf{Y})$$

从而有：$$\frac{\partial E}{\partial \mathbf{B}}=0\Rightarrow {\mathbf{X}}^T\mathbf{XB}={\mathbf{X}}^T\mathbf{Y}$$因此，非齐次线性方程组 ${\mathbf{X}}^T\mathbf{XB}={\mathbf{X}}^T\mathbf{Y}$ 的解即为所求的 $\mathbf{B}$。
如果 ${\mathbf{X}}^T\mathbf{X}$ 可逆，则有 $\mathbf{B}=({\mathbf{X}}^T\mathbf{X}{)}^{-1}{\mathbf{X}}^T\mathbf{Y}$。

上述普通线性最小二乘法在直线拟合中存在以下局限性：

• 无法表示垂直线。 由于模型 $y=mx+b$ 中的斜率 $m$ 无法取无穷大值，因此该方法无法拟合严格垂直分布的数据点（即 $x$ 值恒定），此时方程无解。
• 不具有旋转不变性。最小二乘拟合的结果依赖于坐标系的方向。这意味着，对同一组数据点进行坐标旋转变换后，使用相同方法拟合得到的直线会不同。例如，当数据点经过旋转变换后，在图像平面（或任何坐标系）中呈现出垂直分布的状态时（如上一条所述），该方法将完全失效。
• 对异常点较敏感。损失函数为残差平方和，异常点的残差平方项过大，会显著扭曲拟合直线的方向

2.2. 总体最小二乘法

总体最小二乘法(Total Least Squares, TLS)​不同于最小二乘法，它基于点到直线的垂线距离​（而非 $y$ 值之差）进行拟合，从而有效地解决了前面所述最小二乘法的不足，如下图所示：

从上图可知，总体最小二乘法的解即为齐次线性方程 ${\mathbf{U}}^T\mathbf{UN}=0$ 的解。令 $\mathbf{A}={\mathbf{U}}^T\mathbf{U}$，这是一个 $2\times 2$ 的实对称矩阵，且半正定（其特征值满足 ${\lambda }_1,{\lambda }_2\ge 0$）。我们从特征值与特征向量的关系($\mathbf{AN}={\lambda }_i\mathbf{N}$)出发，结合实对称矩阵的性质，求解 $\mathbf{AN}=0$。实际上，​方程的解对应于最小特征值（通常为零或接近零）对应的归一化特征向量。​

值得注意的是 ${\Vert \left[\begin{array}{cc}{x}_1-\overline{x}& y_1-\overline{y}\\ ⋮& ⋮\\ x_n-\overline{x}& y_n-\overline{y}\end{array}\right]\left[\begin{array}{c}a\\ b\end{array}\right]\Vert }^2$ 为将所有数据点的中心化坐标 $\left(\begin{array}{c}{x}_i-\overline{x}\\ y_i-\overline{y}\end{array}\right)$ 投影到单位直线法向量 $\left(\begin{array}{c}a\\ b\end{array}\right)$ 上，再求这些投影的平方和。这等价于衡量所有数据点在该法向量方向上的总体偏离程度，如下图所示：

总体最小二乘法解决了最小二乘法无法拟合垂直线的问题，但对异常点依然敏感。数据中的异常点很容易显著扭曲拟合直线的方向。如下图所示：

2.3. 从概率角度分析总体最小二乘法

如上图所示，我们从概率角度出发推导总体最小二乘法的目标函数 $E={\sum }_{i=1}^n(a{x}_i+b{y}_i-d{)}^2$。从概率角度，我们可以认为观测数据点 $\left(\begin{array}{c}x\\ y\end{array}\right)$ 是直线上的真值点 $\left(\begin{array}{c}u\\ v\end{array}\right)$ 叠加上垂直于直线的噪声 $ϵ$ 而产生的，即有：
$$\left(\begin{array}{c}x\\ y\end{array}\right)=\left(\begin{array}{c}u\\ v\end{array}\right)+ϵ\left(\begin{array}{c}a\\ b\end{array}\right)$$其中，$ϵ$ 服从均值为 0、方差为 1 的高斯分布，即 $ϵ\sim N(0,1)$。
观测到所有数据点的联合概率可以表示为：
$$P(x_1,y_1,\dots ,x_n,y_n\mid a,b,d)=\prod _{i=1}^{n}P(x_i,y_i\mid a,b,d)\propto \prod _{i=1}^n\exp \left(-\frac{(a{x}_i+b{y}_i-d{)}^2}{2{\sigma }^2}\right)$$又因为：
$$\prod _{i=1}^n\exp \left(-\frac{(a{x}_i+b{y}_i-d{)}^2}{2{\sigma }^2}\right)=\exp \left(\sum _{i=1}^n-\frac{(a{x}_i+b{y}_i-d{)}^2}{2{\sigma }^2}\right)$$则有：
$$\max P(x_1,y_1,\dots ,x_n,y_n\mid a,b,d)\iff \max \sum _{i=1}^n-\frac{(a{x}_i+b{y}_i-d{)}^2}{2{\sigma }^2}\iff \min \sum _{i=1}^n(a{x}_i+b{y}_i-d{)}^2$$

2.4. 鲁棒估计

​最小二乘法和总体最小二乘法在处理包含异常值的数据时，容易使拟合直线发生方向偏移。这是因为异常点对参数估计造成的干扰过大。​为此，​​ 我们引入鲁棒估计以增强拟合的稳健性。
定义鲁棒函数 $\rho (u;\sigma )=\frac{u^2}{{\sigma }^2+u^2}$，其中 $u=r_i(x_i,\theta )$ 代表参数为 $\theta$ 的数学模型在数据点 $x_i$​ 处的残差值。因此，直线拟合问题可归结为求解模型参数 $\theta$，使得 ${\sum }_{i=1}^n\rho (r_i(x_i,\theta );\sigma )$ 最小化。
我们取不同的 $\sigma$ 值绘制曲线 $\rho (u;\sigma )$，如下图所示：
鲁棒函数 $\rho$ 在残差 $u$ 较小时近似二次函数，而在 $u$ 较大时趋于饱和，从而有效抑制了异常点对直线拟合的影响。
结合上图，分析 $\sigma$ 取值对直线拟合的影响：​

• $\sigma$ 过小​：鲁棒函数 $\rho (u;\sigma )$ 对距离 $u$ 的敏感性显著降低，导致不同残差的点具有相近的贡献权重，无法有效区分内点与外点，拟合结果偏离真实模型
• ​​$\sigma$ 过大​：函数 $\rho (u;\sigma )$ 对较大残差的抑制效果减弱，其行为退化为最小二乘法，外点对拟合结果的影响显著增加，失去鲁棒性

值得注意的是：鲁棒估计是一个非线性优化问题，必须通过迭代算法求解，无法直接获得解析解。最小二乘解可作为迭代优化的有效初始值，显著加速收敛并规避局部最优解。尺度参数 $\sigma$ 建议取值为残差绝对值的中位数的1.5 倍，即：$$\sigma =1.5\times median(|r_i|)$$该准则源自统计学理论(Freedman & Peters, 1984)，通过自适应调整平衡模型对核心数据与异常值的响应特性。

3. RANSAC 算法

3.1. 算法原理

鲁棒估计在少量异常值场景下表现良好，但当数据中存在大量异常值时，其稳定性可能显著下降。随机抽样一致算法(Random Sample Consensus, RANSAC)通过迭代随机采样和模型一致性验证，为高比例异常值的数据提供了解决方案。
RANSAC 算法的步骤可以归纳成如下步骤：

1. 参数初始化。需要初始化的参数包括：
   • 模型最小样本数 $s$，其中直线拟合需要 2 个点，平面拟合需要 3 个点
   • 误差阈值 $t$，用于判定内点的误差容忍度
   • 期望置信度概率 $p$ 表示算法在随机抽样过程中至少有一次抽到纯内点集​（无异常值）的概率，通常 $p\ge 0.99$
   • 初始迭代次数 $N$，后续可自适应更新
2. 随机抽样。从数据集中随机选择一组最小样本集，称为 “子集”，用于初始化模型。子集的大小由模型类型决定，例如：直线拟合需要至少 2 个点，平面拟合需要至少 3 个点。值得注意的是：选择的样本必须足够小以保证高效性，同时足够大以覆盖模型的基本自由度；
3. 模型估计。利用随机选择的子集计算模型的初始参数，例如，对于直线拟合，通过两点计算斜率和截距，对于平面拟合，通过三点计算平面方程；
4. 验证一致性。计算所有数据点与当前模型的误差（例如欧氏距离、残差等）。根据预设的误差阈值 $t$，将误差小于 $t$ 的点判定为内点，误差大于 $t$ 的点为异常值(outlier)，并且统计当前模型支持的内点总数 $inlier\_count$；
5. 迭代优化。若当前模型的 $inlier\_count$ > 历史最优模型的 $best\_inlier\_count$，更新 $best\_inlier\_count=inlier\_count$，并且暂存当前模型参数及内点集合。
   计算数据集中异常值的比例 $e=\frac{M-inlier\_count}{M}$，再根据概率公式 $N=\frac{\log (1-p)}{\log (1-(1-e{)}^s)}$ 更新最大迭代次数 $N$，重复执行步骤 $2\sim 5$，逐步优化模型参数。
   迭代的终止条件有两个：
   • 达到最大迭代次数 $N$；
   • 内点数 $best\_inlier\_count$ 超过预设阈值；
6. 最终优化​。用最优模型对应的所有内点，通过最小二乘法等重新估计模型参数（如用所有内点拟合直线）。

上述的步骤可视化如下图：

现在我们来分析在不同外点率 $e$ 和抽样点数量 $s$ 下，为确保 $p=0.99$，所需的迭代次数 $N$。如下表所示。由表可知：当 $s=2$ 时，若 $e=10\%$，则 $N\ge 3$，若 $e=50\%$，则 $N\ge 17$。此外，随着 $s$ 增大，所需的 $N$ 也要相应增大。

总得来说，RANSAC 算法通过随机性和统计验证解决了高噪声环境下的模型拟合问题，实际应用中需根据数据特性调整阈值与迭代策略，以平衡精度与效率。

3.2. 概率公式推导

1. 定义变量。定义如下变量：
   • $p$：期望的成功概率，即在 $N$ 次迭代中至少有一次成功抽到全内点样本的概率
   • $e$：数据集中异常值的比例
   • $s$：模型所需最小样本数（如直线拟合需要 2 个点）
   • $N$：所需迭代次数
2. 计算单次抽样失败的概率。假设数据集中内点的比例为 $1-e$，则单次随机抽样中选择到一个内点的概率为 $1-e$。由于抽样是独立的（假设数据集足够大或允许放回抽样），选择 $s$ 个样本全部为内点的概率为 $(1-e{)}^s$，单次抽样失败的概率为 $1-(1-e{)}^s$；
3. 计算 N 次迭代全部失败的概率。在 $N$ 次独立迭代中，所有迭代均失败的概率为 $[1-(1-e{)}^s{]}^N$；
4. 计算成功概率的表达式。要求至少有一次成功的概率为 $p$，即 $1-[1-(1-e{)}^s{]}^N=p$；
5. 解方程求 $N$。将方程整理为 $1-p=[1-(1-e{)}^s{]}^N$，对两边取对数有 $\log (1-p)=N\cdot \log [1-(1-e{)}^s]$，整理可得 $N=\frac{\log (1-p)}{\log \left[1-(1-e{)}^s\right]}$。

3.3. 仿真实验

​我们采用 RANSAC 算法进行二维点集的直线拟合仿真实验。为简化代码，我们固定了最大迭代次数。仿真代码如下：​

import numpy as np
import matplotlib.pyplot as plt


np.random.seed(6)


def generate_data(num_inliers=100, num_outliers=20, slope=0.5, intercept=1.0, noise_level=0.5):
    """ 
        生成带噪声和异常值的数据。
        生成 100 个内点（符合直线 y = 0.5x + 1）和 20 个异常值（随机分布）
        内点包含高斯噪声，异常值包含更大范围的随机噪声
    """
    # 内点
    x_inliers = np.random.rand(num_inliers) * 10
    y_inliers = slope * x_inliers + intercept + np.random.normal(0, noise_level, num_inliers)

    # 异常值（随机噪声）
    x_outliers = np.random.rand(num_outliers) * 10
    y_outliers = np.random.rand(num_outliers) * 10 + np.random.normal(0, 3, num_outliers)
    
    X = np.hstack((x_inliers, x_outliers))
    y = np.hstack((y_inliers, y_outliers))
    return X, y


def ransac_line_fitting(X, y, num_iterations=100, threshold=1.0, min_samples=2):
    best_model = None
    best_inliers = []
    
    for _ in range(num_iterations):
        # 随机选择两个点作为样
        indices = np.random.choice(len(X), size=min_samples, replace=False)
        x_sample = X[indices]
        y_sample = y[indices]

        # 避免除以零
        if x_sample[0] == x_sample[1]: 
            continue

        # 拟合直线模型 y = kx + b
        k = (y_sample[1] - y_sample[0]) / (x_sample[1] - x_sample[0])
        b = y_sample[0] - k * x_sample[0]
        model = (k, b)

        # 计算所有点到当前模型的残差
        residuals = np.abs(k * X + b - y)
        inliers = np.where(residuals < threshold)[0]

        # 更新最优模型
        if len(inliers) > len(best_inliers):
            best_model = model
            best_inliers = inliers

    # RANSAC 算法拟合直线失败
    if len(best_inliers) < 2:
        return None, None

    # # 使用最优模型的所有内点重新拟合（最小二乘法）
    X_best = X[best_inliers]
    y_best = y[best_inliers]
    A = np.vstack([X_best, np.ones(len(X_best))]).T
    k_best, b_best = np.linalg.lstsq(A, y_best, rcond=None)[0]
    return (k_best, b_best), best_inliers


def ols_line_fitting(X, y):
    A = np.vstack([X, np.ones(len(X))]).T
    k_ols, b_ols = np.linalg.lstsq(A, y, rcond=None)[0]
    return (k_ols, b_ols)


def plot_results(X, y, best_model, best_inliers, outliers_indices, ols_model):
    plt.figure(figsize=(10, 6))
    
    # 绘制 inliers 和 outliers
    plt.scatter(X[outliers_indices], y[outliers_indices], label="Outliers", color="red", marker="x")
    plt.scatter(X[best_inliers], y[best_inliers], label="Inliers", color="blue", alpha=0.7)

    # 绘制 RANSAC 拟合的直线
    x_vals = np.array([np.min(X), np.max(X)])
    y_vals_ransac = best_model[0] * x_vals + best_model[1]
    plt.plot(x_vals, y_vals_ransac, color="green", label="RANSAC Model", linewidth=2)

    # 绘制 OLS 拟合的直线
    y_vals_ols = ols_model[0] * x_vals + ols_model[1]
    plt.plot(x_vals, y_vals_ols, color="orange", label="OLS Model", linestyle="--", linewidth=2)

    plt.text(0.05, 0.95, f"RANSAC: y = {best_model[0]:.2f}x + {best_model[1]:.2f}", 
             transform=plt.gca().transAxes, fontsize=12, verticalalignment='top')
    plt.text(0.05, 0.90, f"OLS:     y = {ols_model[0]:.2f}x + {ols_model[1]:.2f}", 
             transform=plt.gca().transAxes, fontsize=12, verticalalignment='top')

    plt.legend()
    plt.title("RANSAC vs OLS Line Fitting with Outliers")
    plt.xlabel("X")
    plt.ylabel("Y")
    plt.grid(True)
    plt.show()


if __name__ == "__main__":
    # 生成数据
    X, y = generate_data()
    
    # RANSAC 拟合
    best_model, best_inliers = ransac_line_fitting(X, y)

    # 获取异常值的索引
    outliers_indices = np.setdiff1d(np.arange(len(X)), best_inliers)

    # 计算 OLS 拟合的直线
    ols_model = ols_line_fitting(X, y)

    if best_model:
        plot_results(X, y, best_model, best_inliers, outliers_indices, ols_model)
    else:
        print("无法找到合适的模型！")

仿真结果如下图所示，可以看出，当异常值较多时，与最小二乘法相比，RANSAC 算法拟合的直线更优。​