实验3 相机标定实验

组员: 刘欣楠, 翟国辰

班级: 数学强基 2301

电子版密码: zgc-lxn-cvpr3

实验内容

利用棋盘格平面模式, 对定焦相机进行标定 (包括: 相机内参数、相机外参数、相机畸变校正参数).

实验目的

• 1. 掌握针孔相机成像的原理与特点;
• 1. 掌握针孔成像系统中的坐标系 (世界坐标、摄像机坐标、成像平面坐标、图像坐标) 及 各个坐标系之间关系;
• 1. 掌握针孔相机模型中内参数与外参数的物理意义; 掌握相机的经向畸变与切向畸变的特 点及校正方法;
• 1. 掌握基于平面模式的相机标定原理 (张氏标定法) ;
• 1. 统计分析标定误差形成原因, 思考提升标定参数精度的方法.

实验环境与实验器材

硬件: Nvidia 4090laptop, Intel I9-13900HX, 华谷动力相机 WP-YTUT230M, 镜头 WP-5M1614-C, 平面棋盘格 9*12 单个格子边长 100mm.

软件: 华谷动力相机驱动, python (环境如下).

    pip install numpy==1.26.0 scipy==1.11.1 matplotlib==3.8.0 opencv-python==4.8.0.76

实验原理

针孔相机成像模型

针孔相机模型 (Pinhole Camera Model) 是计算机视觉与摄影测量中最基本的成像模型. 其核心思想是通过一个理想的小孔, 将三维空间中的点投影到二维成像平面上, 形成透视投影关系. 假设投影中心为 $O(0,0,0)$, 成像平面位于 $z = f$ 处 ($f$ 为焦距) , 则三维空间中点 $P(X,Y,Z)$ 在成像平面上的投影点 $p(x,y)$ 满足:

$$

\begin{cases} x = f \frac{X}{Z} \ y = f \frac{Y}{Z} \end{cases}

$$

该投影关系体现了“近大远小”的透视特性: 物体越靠近相机, 其像越大. 针孔相机成像属于中心透视投影, 其几何关系如图所示.

此外, 为了消除成像倒置现象, 常在相机前方构造一个虚拟成像平面, 使图像保持正向. 此时投影方程形式保持不变, 但焦距符号相反.

成像系统的多坐标系关系

为便于代数化描述相机成像过程, 通常建立四个基本坐标系:

为了便于从代数上描述摄像机的成像模型, 依据针孔相机成像的具体过程, 建立了四个基本坐标系: - 世界坐标系: 是在3D空间中建立的一个基准坐标系, 描述空间点和摄像机的位置, 如图3-3中建立的世界坐标系 $O_w - X_w Y_w Z_w$; - 摄像机坐标系: 以摄像机中心 $O_c$ 为原点, 以点 $O_c$ 为端点垂直与像平面的射线为主轴 $Z_c$, 以平行与像平面水平线且通过摄像机中心 $O_c$ 的直线为 $X_c$ 轴, 以平行与像平面铅垂线且通过摄像机中心 $O_c$ 的直线为 $Y_c$ 轴, 建立图3-3中的摄像机坐标系 $O_c - X_c Y_c Z_c$; - 成像平面坐标 (也称规范化坐标系) : 以主轴 $Z_c$ 与像平面 (假设为焦距 $f=1$) 的交点为原点 $p$, 该点也称为主点, 以像平面水平线和铅垂线分别为 $x$ 轴和 $y$ 轴, 建立图3-3中的成像平面坐标系 $p - xy$; - 图像坐标: 原点 $o$ 在图像的角上, 平行与 $x$ 轴的 $u$ 轴, 平行于为 $y$ 轴的为 $v$ 轴, 建立下图中的图像坐标系 $o - uv$.

四个坐标系之间的关系: - 世界坐标系通过简单的刚体变换 (旋转和平移) 可以转换成摄像机坐标系; - 摄像机坐标系中的主轴 $Z_c$ 与像平面垂直, 而规范化成像坐标系建立在像平面上; 对摄像机坐标系中的任意一点 $\vec{\mathbf{x}}_c$, 以摄像机中心 $O_c$ 为端点, 与经过点 $\vec{\mathbf{x}}_c$ 构成的射线 $l$ 与像平面的交点 $\mathbf{x}'_{\mathrm{img}}$, 称为 $\vec{\mathbf{x}}_c$ 的像点. 在这个过程中, 3D 空间点 $\vec{\mathbf{x}}_c$ 降维为 2D 平面点 $\mathbf{x}'_{\mathrm{img}}$, 根据透视投影方程, 这两个点在 $x$、$y$ 轴方向上各存在一个尺度缩放因子. 综上摄像机坐标系中的 3D 点到规范化成像坐标系 2D 点的转化关系, 描述了空间 3D 点到成像平面的投影过程. - 成像平面坐标系与图像坐标系建立在像平面上, 区别在于: 成像平面坐标系的原点在主点 $p$ 处, 坐标轴的度量单位为米 ($m$); 而图像坐标系的原点在图像角上, 坐标轴的度量单位为像素个数 (pixel).

基于平面棋盘格的相机标定原理 (张氏标定法)

张氏标定法是张正友博士提出的一种利用平面棋盘格进行相机标定的实用方法, 标定过程仅需要通过相机从不同方向拍摄几张(至少两张)棋盘格图案. 标定原理如下:

1) 棋盘格平面与图像平面之间的单应性 (Homography)

世界坐标系 $O_w - X_w Y_w Z_w$ 建立在棋盘格平面上, $Z_w$ 轴垂直于标板, 即棋盘格平面上的点都有 $Z_w = 0$, 依据针孔相机成像模型, 有:

$$

\mathbf{x} = \begin{bmatrix} u \ v \ 1 \end{bmatrix} \equiv \mathbf{K} \begin{bmatrix} \mathbf{R} & \mathbf{t} \end{bmatrix} \mathbf{X}_w = \mathbf{K} \begin{bmatrix} \mathbf{r}_1 & \mathbf{r}_2 & \mathbf{r}_3 & \mathbf{t} \end{bmatrix} \begin{bmatrix} x_w \ y_w \ 0 \ 1 \end{bmatrix} = \mathbf{K} \begin{bmatrix} \mathbf{r}_1 & \mathbf{r}_2 & \mathbf{t} \end{bmatrix} \begin{bmatrix} x_w \ y_w \ 1 \end{bmatrix} = \mathbf{H} \begin{bmatrix} x_w \ y_w \ 1 \end{bmatrix}

$$

则 $\mathbf{H}$ 为棋盘格平面与图像平面之间的单应性 (Homography) 变换矩阵.

\subsection*{2) 从单应矩阵 $\mathbf{H}$ 中估计相机内参矩阵 $\mathbf{K}$:}

$$

\mathbf{H} = \begin{bmatrix} \mathbf{h}_1 & \mathbf{h}_2 & \mathbf{h}_3 \end{bmatrix} = \lambda \mathbf{K} \begin{bmatrix} \mathbf{r}_1 & \mathbf{r}_2 & \mathbf{t} \end{bmatrix}

$$

$$

\Rightarrow \begin{cases} \mathbf{r}_1 = \dfrac{1}{\lambda} \mathbf{K}^{-1} \mathbf{h}_1 \ \ \mathbf{r}_2 = \dfrac{1}{\lambda} \mathbf{K}^{-1} \mathbf{h}_2 \end{cases} \quad \mathbf{r}_1 \perp \mathbf{r}_2, \; |\mathbf{r}_1| = |\mathbf{r}_2| \quad \Rightarrow \begin{cases} \mathbf{h}_1^T \mathbf{K}^{-T} \mathbf{K}^{-1} \mathbf{h}_2 = 0 \ \ \mathbf{h}_1^T \mathbf{K}^{-T} \mathbf{K}^{-1} \mathbf{h}_1 = \mathbf{h}_2^T \mathbf{K}^{-T} \mathbf{K}^{-1} \mathbf{h}_2 \end{cases}

$$

令

$$

\mathbf{B} = \mathbf{K}^{-T} \mathbf{K}^{-1} = \begin{bmatrix} B_{11} & B_{12} & B_{13} \ B_{12} & B_{22} & B_{23} \ B_{13} & B_{23} & B_{33} \end{bmatrix}, \quad \mathbf{b} = \begin{bmatrix} B_{11}, B_{12}, B_{22}, B_{13}, B_{23}, B_{33} \end{bmatrix}^T

$$

将上述两个等式可写成矩阵相乘的形式 $\mathbf{V} \mathbf{b} = \mathbf{0}$.

通过矩阵的 SVD 分解, 求出向量 $\mathbf{b}$, 再计算出内参矩阵 $\mathbf{K}$ 的各个参数.

3) 估计相机外参矩阵:

$$

\begin{cases} \mathbf{r}_1 = \dfrac{1}{\lambda} \mathbf{K}^{-1} \mathbf{h}_1 \ \ \mathbf{r}_2 = \dfrac{1}{\lambda} \mathbf{K}^{-1} \mathbf{h}_2 \ \ \mathbf{r}_3 = \mathbf{r}_1 \times \mathbf{r}_2 \end{cases} \quad \begin{cases} \mathbf{t} = \dfrac{1}{\lambda} \mathbf{K}^{-1} \mathbf{h}_3 \ \ \lambda = \left| \mathbf{K}^{-1} \mathbf{h}_1 \right| = \left| \mathbf{K}^{-1} \mathbf{h}_2 \right| \end{cases}

$$

4) 径向畸变估计:

畸变模型: $(x,y)$、$(\tilde{x},\tilde{y})$ 分别表示规范化坐标系下理想和畸变后的坐标点.

$$

\begin{aligned} \tilde{x} &= x + x \left[ k_1(x^2 + y^2) + k_2(x^2 + y^2)^2 + k_3(x^2 + y^2)^3 \right] \ \tilde{y} &= y + y \left[ k_1(x^2 + y^2) + k_2(x^2 + y^2)^2 + k_3(x^2 + y^2)^3 \right] \end{aligned}

$$

$(u,v)$、$(\tilde{u},\tilde{v})$ 分别表示图像像素坐标系下理想和畸变后的坐标点, 经内参矩阵 $\mathbf{K}$ 将上述模型变为:

$$

\begin{aligned} \tilde{u} &= u + (u - u_0) \left[ k_1(x^2 + y^2) + k_2(x^2 + y^2)^2 \right] \ \tilde{v} &= v + (v - v_0) \left[ k_1(x^2 + y^2) + k_2(x^2 + y^2)^2 \right] \end{aligned}

$$

写成矩阵相乘形式:

$$

\underbrace{ \begin{bmatrix} (u - u_0)(x^2 + y^2) & (u - u_0)(x^2 + y^2)^2 \ (v - v_0)(x^2 + y^2) & (v - v_0)(x^2 + y^2)^2 \end{bmatrix} }{\mathbf{D}} \underbrace{ \begin{bmatrix} k_1 \ k_2 \end{bmatrix} } = \underbrace{ \begin{bmatrix} \tilde{u} - u \ \tilde{v} - v \end{bmatrix} }_{\mathbf{d}}}

$$

用线性最小二乘求解得

$$

\mathbf{k} = \left( \mathbf{D}^T \mathbf{D} \right)^{-1} \mathbf{D}^T \mathbf{d}

$$

5) 极大似然估计:

用 $n$ 张棋盘格图片对相机进行标定, 将重投影误差设定为目标函数:

$$

S(\mathbf{K}, \mathbf{k}, \mathbf{R}, \mathbf{t}, \lambda) = \sum_{i=1}^{n} \sum_{j=1}^{m} \left| \mathbf{x}_{ij} - f(\mathbf{K}, \mathbf{k}, \mathbf{R}_i, \mathbf{t}_i, \lambda_i, \mathbf{X}_j) \right|^2

$$

通过非线性最小二乘法, 优化目标函数.

Levenberg-Marquardt-LM 优化法:

目标函数: $ S(\boldsymbol{\beta}) = \sum_{i=1}^{n} \left| \mathbf{y}_i - f(\mathbf{x}_i, \boldsymbol{\beta}) \right|^2 $

一阶 Jacobi 矩阵近似: $ f(\mathbf{x}_i, \boldsymbol{\beta} + \boldsymbol{\delta}) \approx f(\mathbf{x}_i, \boldsymbol{\beta}) + \mathbf{J}_i \boldsymbol{\delta} $, 其中 $ \mathbf{J}_i = \dfrac{\partial f(\mathbf{x}_i, \boldsymbol{\beta})}{\partial \boldsymbol{\beta}} = \dfrac{\partial \mathbf{y}_i}{\partial \boldsymbol{\beta}} $

目标函数 $S$ 的周边信息:

$$

S(\boldsymbol{\beta} + \boldsymbol{\delta}) = \sum_{i=1}^{m} \left| \mathbf{y}i - f(\mathbf{x}_i, \boldsymbol{\beta} + \boldsymbol{\delta}) \right|^2 = \sum \right|^2}^{m} \left| \mathbf{y}_i - f(\mathbf{x}_i, \boldsymbol{\beta}) - \mathbf{J}_i \boldsymbol{\delta

$$

使 $S$ 最小时: $ \left( \mathbf{J}^T \mathbf{J} \right) \boldsymbol{\delta} = \mathbf{J}^T \left[ \mathbf{Y} - F(\boldsymbol{\beta}) \right] $

其中

$$

\mathbf{J} = \begin{bmatrix} \dfrac{\partial y_1}{\partial \beta_1} & \dfrac{\partial y_1}{\partial \beta_2} & \cdots & \dfrac{\partial y_1}{\partial \beta_n} \ \dfrac{\partial y_2}{\partial \beta_1} & \dfrac{\partial y_2}{\partial \beta_2} & \cdots & \dfrac{\partial y_2}{\partial \beta_n} \ \vdots & \vdots & \ddots & \vdots \ \dfrac{\partial y_m}{\partial \beta_1} & \dfrac{\partial y_m}{\partial \beta_2} & \cdots & \dfrac{\partial y_m}{\partial \beta_n} \end{bmatrix}

$$

加入阻尼系数: $ \left( \mathbf{J}^T \mathbf{J} + \lambda \mathbf{I} \right) \boldsymbol{\delta} = \mathbf{J}^T \left[ \mathbf{y} - F(\boldsymbol{\beta}) \right] $,

得到: $ \boldsymbol{\delta} = \left( \mathbf{J}^T \mathbf{J} + \lambda \mathbf{I} \right)^{-1} \mathbf{J}^T \left[ \mathbf{y} - F(\boldsymbol{\beta}) \right] $.

实验步骤与结果分析

采集棋盘格图片

使用华谷动力相机 WP-YTUT230M, 镜头 WP-5M1614-C, 平面棋盘格 9*12 单个格子边长 100mm. 共采集了 46 张清晰不同位姿的图片.

思考问题:

1) 理论上只需要 3 张图片, 因为每张图片能提供 2 个关于内参的独立约束.

2) 不同的位姿是为了更好的区分内参与外参, 避免相同位姿导致的参数退化以至于无法求解.

提取棋盘格角点

我们使用 python 语言中的 opencv 包对图片进行角点提取. 以下是两组示例图.

可以注意到, 标注出的角点精度很高.

思考问题:

1) 在 OpenCV 实现中, 回按行顺序返回角点, 我们用同样的方式对棋盘上的角点编号, 可以看到, 在上述两组图片中都以坐下角作为红色起点.

2) 首先二者数量要相同, 这就要求采集的照片一定要包括整个棋盘. 其次棋盘可以清晰确定顺序, 避免因旋转导致的顺序错乱.

求解单应性矩阵

下面是两张图的单应性矩阵示例.

$$

\begin{bmatrix} -9.413411 & 95.165331 & 432.481612 \ -66.974114 & 5.141156 & 877.766191 \ 0.014179 & 0.010800 & 1.000000 \end{bmatrix},\quad \begin{bmatrix} -4.172529 & 95.518424 & 603.945823 \ -77.134468 & 41.120396 & 737.451313 \ 0.016976 & 0.020178 & 1.000000 \end{bmatrix}

$$

估计相机的内参和外参

1) 初始化内参.

从初始标定得到的内参矩阵 $K$：

$$

K = \begin{bmatrix} 3471.651308, & 0.000000, & 982.124233 \ 0.000000, & 3474.180514, & 599.359246 \ 0.000000, & 0.000000, & 1.000000 \end{bmatrix}

$$

说明: - $fx = 3471.651308 \text{ px}$ - $fy = 3474.180514 \text{ px}$ - 主点 $(cx, cy) = (982.124233 \text{ px}, 599.359246 \text{ px})$ - 将像素单位的焦距换算为物理焦距 (使用相机规格的像素间距 pixel\_pitch = 4.8 $\mu$m = 0.0048 mm)： - $\text{pixel\_size} = 4.8 \text{ } \mu\text{m} = 0.0048 \text{ mm}$ - $fx_{\text{mm}} = fx \times \text{pixel\_size} = 3471.651308 \times 0.0048 \approx 16.6639 \text{ mm}$ - $fy_{\text{mm}} = fy \times \text{pixel\_size} = 3474.180514 \times 0.0048 \approx 16.6761 \text{ mm}$ - 传感器尺寸 (由像素尺寸和分辨率算出)： - $\text{sensor\_width} = 1920 \times 0.0048 = 9.2160 \text{ mm}$ - $\text{sensor\_height} = 1200 \times 0.0048 = 5.7600 \text{ mm}$ - 由焦距计算视场 (FOV)： - $FOV_x = 2 \times \text{atan}(\text{sensor\_width} / (2 \times f_{\text{mm}}))$ - 规格值 (用 $f_{\text{spec}} = 16.0 \text{ mm}$): $FOV_{x,\text{spec}} \approx 32.133^\circ$ - 测得 (用 $f_{\text{mm}} \approx 16.6639 \text{ mm}$): $FOV_{x,\text{meas}} \approx 30.915^\circ$ - $FOV_y = 2 \times \text{atan}(\text{sensor\_height} / (2 \times f_{\text{mm}}))$ - 规格值 (用 $f_{\text{spec}} = 16.0 \text{ mm}$): $FOV_{y,\text{spec}} \approx 20.408^\circ$ - 测得 (用 $f_{\text{mm}} \approx 16.6639 \text{ mm}$): $FOV_{y,\text{meas}} \approx 19.597^\circ$

参数对比:

`latex \begin{tabular}{|l|l|l|l|l|}

参数 & 规格值 & 测量值 & 绝对差 & 相对差 \

fx(px) & 3333.33 px & 3471.651 px & +138.32 px & +4.15\% \

fy(px) & 3333.33 px & 3474.181 px & +140.85 px & +4.23\% \

f (mm) & 16.00 mm & 16.6639 mm & +0.6639 mm & +4.15\% \

cx (px) & 960.0 px & 982.124 px & +22.12 px & +1.152\% \

cy (px) & 600.0 px & 599.359 px & -0.64 px & -0.053\% \

FOV_x & 32.133° & 30.915° & -1.218° & -3.79\% \

MSE & --- & 0.463979 & --- & --- \

\end{tabular} `

2) 计算外参.

rvec = (-0.1623, -0.8116, -1.6377)

tvec = (-5.9215, 2.9969, 37.3366)

rot_deg = 105.14

trans_norm = 37.92

分析标定误差

    \begin{tabular}{|c|c|c|c|}

        **标定图像数量** & **组别** & **重投影误差 (px)** & **主点坐标 [cx cy]** \\

        \multirow{4}{*}{3}
        & 第一组 & 1.5836 & [981.05\ \ 635.81] \\
        & 第二组 & 1.3101 & [883.89\ \ 626.09]\\
        & 第三组 & 0.4751 & [1006.99\ \ 605.07] \\
        & 随机 50 次均值 & 0.8243 & [970.88\ \ 603.85] $\pm$ [54.79\ \ 48.79]\\

        \multirow{4}{*}{7}
        & 第一组 & 0.4551 & [1010.55\ \ 538.60] \\
        & 第二组 & 0.3844 & [976.00\ \ 557.00] \\
        & 第三组 & 0.6414 & [1008.33\ \ 583.30] \\
        & 随机 50 次均值 & 0.49 & [979.58\ \ 600.85] $\pm$ [38.03\ \ 27.42]\\

    \end{tabular}

如表格所示, 图片数量的增加可以显著提高标定实验的稳定性, 但重投影误差并不一定会显著减少, 说明数量不是唯一因素. 结合之前对图片采集的限制, 位姿的多样性也极大的影响着标定的精度.

误差形成的原因可能是存在标定偏移大的图片干扰, 数据采集多样性不足.

提升方法:

1. 剔除显著差异过大的图片数据.

2. 从旋转, 俯仰, 距离等多位姿情况采集数据, 避免只是单纯重复相似的位姿.

3. 选择不同子集进行多次重复模拟以确定稳定性.

结论与讨论

通过本次实验, 成功完成了从棋盘格角点提取、相机标定、单应矩阵求解到相机内参与畸变参数估计的完整标定流程. 利用多组标定图像实现了相机的几何建模与参数求解. 验证了张正友标定法在常规实验条件下的有效性.

当标定图像数量较少 (如 3 张)时, 结果受角点检测精度与图像噪声影响较大, 导致 MSE 误差及主点位置波动明显. 当图像数量增加至 7 张或以上后, 参数估计结果趋于稳定, 表明冗余观测对降低随机误差起到了关键作用.

主点位置略偏离图像中心 (约 $10\sim30$ 像素), 可能由镜头装配偏心或相机几何不对准引起, 此类偏差在实际工业相机中属于正常现象. 焦距估计值与理论焦距 $f=16\,\text{mm}$ 误差较小, 验证了模型假设的正确性.

本实验使我们对相机成像模型及标定算法有了更深入的理解. 通过实践操作, 认识到理论模型虽理想化, 但实际数据受光照, 噪声, 镜头畸变等因素影响明显; 标定精度不仅依赖算法本身, 更依赖实验图像质量及拍摄几何布置.