Fisher 信息：相机标定仿真实操

本文最后更新于：2026年8月5日下午

相机标定估出内参后，这些估计准到什么程度？本文用 Fisher 信息算出每个内参的精度下界（CRLB）——不管用什么算法，误差都不可能小于它，而且拍照前就能算。理论推导见《Fisher 信息与精度下界》。

问题描述

标定给出一组内参（焦距、主点、畸变系数），但估得准不准、误差能有多大，通常只能凭重投影残差粗略感受。我们想要一个更根本的量：在既定的拍摄条件下，内参估计的精度极限是多少。

这个极限由 Fisher 信息给出，叫 Cramér–Rao 下界（CRLB）。它只依赖观测方程和噪声，与具体算法无关，因此可以在拍照之前就算出来——用来规划拍多少张、什么角度；拍完也能反过来检验标定是否真的逼近了理论极限。

只关心内参时，第 k 个内参的精度下界是

$$ \mathrm{CRLB}_k=\sqrt{[\mathcal I_\beta^{-1}]_{kk}} $$

其中 $\mathcal I_\beta$ 是扣掉外参之后、只关于内参的信息矩阵。整套推导见理论篇，下面只说怎么把它一步步算出来。

解决思路

高斯噪声下，信息矩阵只看雅可比——每个像素对各参数的偏导有多灵敏：

$$ \mathcal I=\mathbf J^\top\boldsymbol\Sigma^{-1}\mathbf J $$

这一行把"算精度下界"拆成三件事：

雅可比 $\mathbf J$：参数如何影响像素。数值差分就能算，不必手推偏导。
噪声 $\boldsymbol\Sigma$：角点检测的随机误差。各向同性、各通道不相关时退化为对角阵，对角元是单轴噪声方差 σ_a²。
边缘化外参：每张图的位姿也是从同一批数据估出来的，会分走内参的信息，要用 Schur 补扣掉。

三件事自然落成六个步骤：先造出带真值的仿真数据，算雅可比，累加成信息矩阵，边缘化外参，求逆开方得 CRLB，最后用 Monte Carlo 验证估计器是否真达到这个下界。前一步的产出就是后一步的输入。

仿真与真实

这条流水线算的是金标准 CRLB——它在真值处求值，只有仿真算得出（仿真里真值已知）。真实标定没有真值，做法是先标定得到估计 $\hat{\boldsymbol\theta}$，再在估计处算信息（观测信息），不确定度近似为

$$ \mathrm{Cov}(\hat{\boldsymbol\theta})\approx\mathcal I_{\mathrm{obs}}(\hat{\boldsymbol\theta})^{-1} $$

OpenCV 的 stdDeviations 就是这么来的。这之所以可行，靠的是最大似然的一致性：估计会收敛到真值，故观测信息与期望信息在大样本下重合。

那么仿真何必费劲算金标准？恰恰是为了验证：拿真值处的 CRLB 去对照真实标定算法的散布，第 6 步里 std/CRLB ≈ 1 就说明这套方法可信，真实标定里用观测信息法才站得住脚。

无真值时的真实标定方法

上面说真实标定走观测信息——具体怎么算？几乎和仿真一模一样，只把"真值"换成"估计"：θ_true 换成 θ̂，噪声 σ 换成重投影残差估出的 σ̂，其余代码原样。

步	仿真（有真值）	真实标定（无真值）
数据	正向仿真生成	真实角点检测 `y`
求值点	真值 θ	估计 θ̂（`calibrateCamera` 输出）
噪声	`σ_a` 已知	`σ̂_a` = 重投影 RMS / √2
雅可比	`J(θ)` 在真值处	`J(θ̂)`（同一套差分代码）
信息	`(1/σ_a²)JᵀJ` 在真值处	`(1/σ̂_a²)J(θ̂)ᵀJ(θ̂)`
结果	CRLB = `√diag(I_β⁻¹)`	`Cov(θ̂) ≈ I_obs⁻¹`
验证	第 6 步 Monte Carlo	做不了（无真值对照）

不确定度的核心式：

$$ \mathrm{Cov}(\hat{\boldsymbol\theta})\approx\mathcal I_{\mathrm{obs}}(\hat{\boldsymbol\theta})^{-1}=\bigl(\hat\sigma_a^{-2}\,\mathbf J(\hat{\boldsymbol\theta})^\top\mathbf J(\hat{\boldsymbol\theta})\bigr)^{-1} $$

噪声怎么估：calibrateCamera 返回的重投影 RMS 就是径向噪声 std σ̂_det，单轴 σ̂_a = RMS/√2（和仿真里 σ_a = σ_det/√2 同一个换算）。

唯一丢掉的是第 6 步验证。仿真有真值，才能拿金标准 CRLB 去对照估计器散布；真实标定没有真值，这一步天然做不了。但这正是仿真的意义——它替你验证过"在估计处算观测信息"和"在真值处算金标准"几乎重合（std/CRLB ≈ 1），所以真实标定里用观测信息法才站得住。

关心外参时

边缘化外参只是因为我们关心内参。若关心位姿（外参）精度，反过来：

完整协方差：不边缘化，Cov(θ̂) ≈ I⁻¹ 给出所有参数的联合协方差，第 i 张图的位姿不确定度就是 I⁻¹ 里对应那 6 个位姿参数的块。
内参已知、只估位姿（PnP）：参数只剩 6 个位姿，I 是 6×6，Cov = I⁻¹ 直接给位姿精度——PnP 不确定度、手眼标定误差传播的基础。
内参也估、带不确定度：把内参当 nuisance 边缘化掉，得到扣掉内参不确定性后的纯位姿精度，比假装内参已知更诚实。

其实你已经在用

cv2.calibrateCamera 返回的 stdDeviations，内部就是这套：在 θ̂ 处建 J、JᵀJ、按残差方差缩放、求逆取对角。所以真实标定里早就在用观测信息法了，只是 OpenCV 把它打包成那个返回值。这套流水线的额外价值是：把边缘化（Schur 补）和完整协方差块显式算出来——stdDeviations 只给对角，拿不到"扣掉外参后内参的 CRLB"或"某张图位姿的 6×6 协方差"。

什么时候失准

数据太少：θ̂ 离真值远，渐近近似变差 → 多拍、位姿多样。
系统误差：角点检测偏置、镜头模型不匹配（真实畸变 ≠ 假设的 5 系数）→ CRLB 只管随机误差，偏置要另算；有偏时 CRLB 不再是有效下界。
噪声非高斯：CRLB 形式仍对，但渐近有效要打折扣。

配置

具体怎么干，先定一组贴近真实的参数：

图像 2560×1440，焦距 fx = fy = 2200
棋盘 9×6 内角点，方格 25mm
50 张整板入画的图，位置、俯仰、偏航、滚转、远近都不同（多样性是避免参数耦合退化的关键）
角点噪声 σ_det = 0.1px（径向），折合单轴 σ_a = σ_det/√2

文中代码与项目 pipeline/ 下脚本逐字一致，所有数字真跑过。

数据流

step1                step2              step3            step4           step5          step6
生成数据   ──calib_sim.npz──>  算雅可比  ──J.npz──>  算信息  ──I.npz──>  边缘化  ──Ibeta.npz──>  CRLB  ──crlb.npz──┐
                                                                                                                 │
             ┌────────────────── calib_sim.npz（取板、f、σ_a）+ crlb.npz ←──────────────────────────────────────────┘
             ▼
           step6：Monte Carlo 验证  ──mc_stats.npz

每个 .npz 存的是上一步的结果，下一步读它，文件名即内容。

步	脚本	做什么	输入 → 输出
1	`step1_gen_data.py`	生成仿真数据（真值、位姿、f、y）	— → `calib_sim.npz`
2	`step2_jacobian.py`	雅可比 J = ∂f/∂θ	`calib_sim.npz` → `J.npz`
3	`step3_fisher.py`	信息矩阵 I = (1/σ_a²)·JᵀJ	`J.npz` → `I.npz`
4	`step4_marginal.py`	Schur 补边缘化外参 → I_β	`I.npz` → `Ibeta.npz`
5	`step5_crlb.py`	求逆开方 → CRLB	`Ibeta.npz` → `crlb.npz`
6	`step6_validate.py`	Monte Carlo 验证 std/CRLB ≈ 1	`calib_sim.npz`+`crlb.npz` → `mc_stats.npz`

理论符号 ↔ 标定里的真实量

理论符号	标定里的真实量	本例取值
θ	内参（9 个）+ 每图外参	(2200, 2200, 1280, 720, −0.1, 0, 0, 0, 0) + 50 组位姿
y	全部角点像素 (u, v)	50 图 × 54 角 × 2 = 5400 个数
f(θ)	3D 角点投影出的 (u, v)	`projectPoints` 链
J	各角点对参数的偏导	5400 × 309
Σ	角点方差	σ_a²·I = 0.005·I
σ_a / σ_det	单轴 / 径向噪声 std	0.0707 / 0.1 px
I / I_β	全参 / 内参信息矩阵	309² / 9²
CRLB_k	第 k 内参精度下界	fx → 0.154 px

实操流程

第 1 步｜生成测试数据

做什么：正向仿真——已知真值内参和位姿，把 3D 角点投影成"应有像素" f(θ_true)，再加噪声得观测 y。全程没有估计。正向链是板坐标经外参变到相机系、除以深度归一化、再经畸变和内参得到像素。噪声 σ_a = σ_det/√2，因为径向误差平方等于两个轴向误差平方之和。

输入 → 输出：— → data/calib_sim.npz（真值内参、板、外参、img_clean=f、img_noisy=y、σ）。

位姿反解（怎么从"板占多大、中心在哪"反推位姿，只看针孔主项）。由板占视场比例 fill 反解距离：针孔下板宽成像像素等于焦距乘板宽除以距离，要板宽占 fill·W 像素，故距离 Z = fx·board_w / (fill·W)，fill 大则板大则离得近。中心放置反解平移 t：把板几何中心投到指定像素，先反针孔算出中心应有的相机坐标，再由 C_c = R·C_w + t 解出 t = 目标相机坐标 − R·板中心，板向四周铺开、容易整块入画。整板入画过滤：只留所有角点都在画面内的位姿——真实标定只看得见画面内角点、只关心画面内畸变，部分出画的板不可用（曾用"3 倍画幅外才挡"导致四成角点出画、虚高精度）。参数独立多样性：位置、俯仰、偏航、滚转、远近用混合基数各自独立遍历，不锁死联动——多样性是打破参数耦合退化的关键。

代码（逐行注释，与 pipeline/step1_gen_data.py 一致）：

"""
Step 1：生成相机标定仿真测试数据
==================================
本脚本只做"正向仿真"：已知真值内参 θ_true 和每张图的位姿，把 3D 棋盘角点
投影成"应有的像素坐标" f(θ_true)，再加噪声，得到"观测像素" y。
全程纯几何 + 加噪，没有任何估计——真值是已知的，这正是仿真的意义。

输出的 npz 文件含（后续步骤都从它读）：
  theta_true / K / dist   真值内参（当 θ_true 用）
  board                   棋盘 3D 角点 (N,3)，Z=0 平面，单位 mm
  rvecs / tvecs           每张图外参（旋转向量、平移）
  img_clean               无噪声投影 f(θ_true)，形状 (n_img, N, 2)
  img_noisy               加噪观测 y = f(θ_true)+噪声，形状 (n_img, N, 2)
  sigma_det / sigma_a     噪声 std（径向 / 单轴）

配置参考 sim 基线（本流水线改用整板入画 + 中心放置 + 多样性采样）；CRLB_fx≈0.154。
运行： python3 pipeline/step1_gen_data.py
"""
import os                                  # 仅用于拼保存路径（让脚本从任意目录运行都能存对地方）
import numpy as np                         # 数值计算（数组、矩阵、随机）
import cv2                                 # OpenCV：projectPoints（投影）、Rodrigues（旋转向量↔矩阵）

rng = np.random.default_rng(0)             # 随机数发生器，种子=0 → 每次运行噪声一样、结果可复现

# 1. 真值内参 θ_true
# 这一组数就是"上帝视角"的真相，后续 CRLB 在它上面算；估计器试图还原的也是它。
theta_true = dict(fx=2200., fy=2200., cx=1280., cy=720.,   # 焦距 fx/fy（像素）、主点 cx/cy（像素）
                  k1=-0.10, k2=0., p1=0., p2=0., k3=0.)    # 畸变系数：k 径向、p 切向（这里只给 k1）
image_size = (2560, 1440)                  # 图像宽×高 (W, H)，像素
dist_keys = ["k1", "k2", "p1", "p2", "k3"] # OpenCV 要求的 5 系数顺序（写错位置会张冠李戴）


def theta_to_cv(th):
    """把内参 dict 转成 OpenCV 要的 (K 3×3 相机矩阵, distCoeffs 畸变向量)。"""
    K = np.array([[th["fx"], 0, th["cx"]],            # 相机内参矩阵 K：[[fx,0,cx],[0,fy,cy],[0,0,1]]
                  [0, th["fy"], th["cy"]],
                  [0, 0, 1.]])
    dist = np.array([th[k] for k in dist_keys])       # 畸变向量按 dist_keys 顺序取值 → [k1,k2,p1,p2,k3]
    return K, dist


K, dist = theta_to_cv(theta_true)         # 算出本例的 K、dist，后面投影要用

# 2. 棋盘 3D 角点
# 标定板的"内部角点"在板自身坐标系下的 3D 坐标（Z=0 平面，单位 mm）。
cols, rows = 9, 6                          # 内角点 列数×行数（9×6 板共 54 个内角点）
sq = 25.0                                  # 每个方格的边长 mm（实物用尺子量）
board = np.zeros((rows * cols, 3))         # 预分配 (54,3)，先全零
for r in range(rows):                      # 逐行
    for c in range(cols):                  # 逐列
        board[r * cols + c] = (c * sq, r * sq, 0.)   # 第 (r,c) 个角点的板坐标 (X,Y,Z=0)，单位 mm

# 3. 正向投影 f
def project(rvec, tvec, X):
    """3D 点 X (N,3) → 像素 (N,2)：经外参→畸变→内参。这一整条链就是正向模型 f(θ_true)。"""
    pts, _ = cv2.projectPoints(X.reshape(-1, 1, 3),       # 3D 点，reshape 成 OpenCV 要的 (N,1,3)
                               rvec.reshape(3, 1), tvec.reshape(3, 1),   # 旋转向量、平移向量（外参）
                               K, dist)                    # 内参矩阵、畸变系数
    return pts.reshape(-1, 2)              # 投影结果压回 (N,2)，每行是该角点的像素 (u,v)


def euler_to_R(pitch, yaw, roll):
    """欧拉角(度) → 3×3 旋转矩阵 R（板坐标系→相机坐标系的姿态）。"""
    p, y, r = map(np.radians, (pitch, yaw, roll))         # 度→弧度
    Rx = np.array([[1, 0, 0], [0, np.cos(p), -np.sin(p)], [0, np.sin(p), np.cos(p)]])   # 绕 X 轴(俯仰)
    Ry = np.array([[np.cos(y), 0, np.sin(y)], [0, 1, 0], [-np.sin(y), 0, np.cos(y)]])   # 绕 Y 轴(偏航)
    Rz = np.array([[np.cos(r), -np.sin(r), 0], [np.sin(r), np.cos(r), 0], [0, 0, 1]])   # 绕 Z 轴(滚转)
    return Rz @ Ry @ Rx                    # 合成：先俯仰、再偏航、再滚转（矩阵乘顺序=作用反序）


# 4. 采样位姿
# 思路：把板的【几何中心】投到画面不同位置（九宫格）、给不同倾斜/滚转、板占视场宽
# 不同比例，得到覆盖多样的 50 张图。距离 Z 由"板宽占视场比例"反解；平移 t 由
# "板中心投到指定像素"反解（中心放置：t = 目标相机坐标 - R·板中心，这样板向四周铺开）。
# 真实约束：只保留【所有角点都在画面内】的位姿——真实标定只能看见画面内的角点、
# 也只关心画面内的畸变，部分出画的板不可用；这同时挡掉了畸变爆裂的荒谬投影。
W, H = image_size                          # 画面宽、高
fx, fy, cx, cy = theta_true["fx"], theta_true["fy"], theta_true["cx"], theta_true["cy"]  # 取出便于书写
n_images = 50                              # 要采多少张图（多图累加信息，见可加性；小板需更多图）
board_w_mm = (cols - 1) * sq               # 板的物理宽度 mm（角点跨度，非方格数）
center_w = np.array([(cols - 1) * sq / 2, (rows - 1) * sq / 2, 0.])   # 板的几何中心（板坐标）
grid = [(x, y) for y in np.linspace(0.3, 0.7, 3) for x in np.linspace(0.3, 0.7, 3)]   # 板中心投到画面的归一化位置（九宫格，位置多样性）
pitches = np.linspace(-30, 30, 3)          # 俯仰角候选（度）：-30, 0, 30
yaws = np.linspace(-30, 30, 3)             # 偏航角候选（度）
rolls = [-45, 0, 45]                       # 滚转角候选（度）
fills = np.linspace(0.35, 0.55, 3)         # 板占视场宽比例（现实大小：板约 1/3~1/2 画面宽；配多图累加信息）

rvecs, tvecs = [], []                      # 收集合格位姿：旋转向量、平移向量
i = 0                                      # 候选计数（过采样，凑够 50 个为止）
while len(rvecs) < n_images and i < n_images * 40:   # 凑够 50 张，或候选用尽就停
    # 混合基数独立遍历：让 位置/俯仰/偏航/滚转/远近 各自独立变化（全交叉组合），
    # 而非锁死联动——多样性是打破参数耦合退化的关键。
    fi = i % 3                              # 远近（fill）变得最快
    ri = (i // 3) % 3                       # 滚转
    yi = (i // 9) % 3                       # 偏航
    pi = (i // 27) % 3                      # 俯仰
    gi = (i // 81) % len(grid)              # 位置变得最慢
    gx, gy = grid[gi]                       # 板中心投到画面的归一化位置
    pitch, yaw, roll, fill = pitches[pi], yaws[yi], rolls[ri], fills[fi]  # 各自独立取
    Z = board_w_mm * fx / (fill * W)       # 由"板宽占视场 fill"反解相机到板的距离 Z(mm)：板宽像素 = board_w*fx/Z = fill*W
    R = euler_to_R(pitch, yaw, roll)       # 由欧拉角算旋转矩阵
    # 把板【中心】投到 (gx*W, gy*H)：先反针孔算出中心应有的相机坐标，再 t = 目标 - R·中心
    target_cx = (gx * W - cx) * Z / fx     # 中心的相机系 X 坐标（使其投到 gx*W）
    target_cy = (gy * H - cy) * Z / fy     # 中心的相机系 Y 坐标（使其投到 gy*H）
    t = np.array([target_cx, target_cy, Z]) - R @ center_w   # 平移：目标相机坐标 - R·板中心
    rvec, _ = cv2.Rodrigues(R); rvec = rvec.flatten()   # 旋转矩阵 → OpenCV 用的旋转向量(3,)
    u = project(rvec, t, board)            # 用这组位姿把整块板投到像素，看落点合不合理
    # 真实约束：所有角点必须在画面 [0,W]×[0,H] 内（任一出画即丢弃该位姿）
    if (np.all(np.isfinite(u))             # 投影没有 NaN/inf
            and u[:, 0].min() >= 0 and u[:, 0].max() <= W    # u 方向全在画面内
            and u[:, 1].min() >= 0 and u[:, 1].max() <= H):  # v 方向全在画面内
        rvecs.append(rvec); tvecs.append(t)   # 合格：整板在画面内，收下这组位姿
    i += 1                                 # 试下一组候选
if len(rvecs) < n_images:                  # 没凑够说明参数太紧（fill 太大 / 倾斜太大 / grid 太靠边）
    print(f"[warn] 只采到 {len(rvecs)}/{n_images} 个整板入画的位姿，考虑降 fill/倾斜 或收紧 grid")
rvecs = np.array(rvecs); tvecs = np.array(tvecs)   # 列表 → 数组，形状 (n_img, 3)
n_img = len(rvecs)                         # 实际采到的图数（应为 50）

# 5. 无噪声投影 f(θ_true)
# 对每张图，用它的外参把整块板投到像素，得到"应有坐标"——这就是 f(θ_true)。
img_clean = np.stack([project(rvecs[j], tvecs[j], board) for j in range(n_img)])   # 堆叠成 (n_img, 54, 2)

# 6. 加噪 → 观测 y
# 真实角点检测器有随机误差，用高斯噪声模拟。注意单轴 σ_a 与径向 σ_det 的换算。
sigma_det = 0.1                            # 径向噪声 std（像素），与重投影 RMS 同刻度；可由静态重复实测得到
sigma_a = sigma_det / np.sqrt(2)           # 单轴噪声 std（u、v 各自的抖动）
# 径向误差² = Δu² + Δv² = σ_a² + σ_a² = 2σ_a² = σ_det²，故 σ_a = σ_det/√2
noise = rng.normal(0, sigma_a, img_clean.shape)   # 抽同形状的高斯噪声，std=σ_a
img_noisy = img_clean + noise              # 观测 = 应有坐标 + 噪声 = 角点检测器会给出的 y

# 7. 保存
HERE = os.path.dirname(os.path.abspath(__file__))   # 本脚本所在目录（不受运行时 cwd 影响）
DATA = os.path.join(HERE, "data")                   # 数据子目录 pipeline/data
os.makedirs(DATA, exist_ok=True)                    # 不存在就建（已存在不报错）
out = os.path.join(DATA, "calib_sim.npz")           # 输出文件路径
theta_keys = ["fx", "fy", "cx", "cy"] + dist_keys   # 内参的展开顺序，存下来便于后续按名取值
np.savez(out,                                     # 一次性存成 npz（后续 np.load 读）
         theta_true=np.array([theta_true[k] for k in theta_keys]),   # 真值内参向量(9,)
         theta_keys=np.array(theta_keys),          # 内参名称，与上一行一一对应
         K=K, dist=dist,                           # OpenCV 形式的内参矩阵、畸变
         board=board, rvecs=rvecs, tvecs=tvecs,    # 板 3D 角点、每图外参
         img_clean=img_clean, img_noisy=img_noisy, # f(θ_true) 与 观测 y
         sigma_det=sigma_det, sigma_a=sigma_a,     # 噪声 std（两种刻度）
         image_size=np.array(image_size),           # 画面尺寸
         pattern=np.array([cols, rows]), square_mm=sq)   # 棋盘规格

# 8. 打印:生成内容与计算方式
print("=" * 64)                            # 分隔线
print(f"图像 {W}×{H}，棋盘 {cols}×{rows}（{sq:g}mm），{n_img} 张图，每图 {rows*cols} 个角点")  # 配置概览
print(f"真值内参 fx=fy={theta_true['fx']:.0f}, cx={cx:.0f}, cy={cy:.0f}, k1={theta_true['k1']}")  # 真值
print(f"噪声 σ_det={sigma_det}px → 单轴 σ_a={sigma_a:.4f}px  (σ_a = σ_det/√2)")  # 噪声换算
print(f"无噪声投影 img_clean 形状 {img_clean.shape} = (图数, 角点数, 2)  ← f(θ_true)")  # f 的形状
print(f"加噪观测   img_noisy 形状 {img_noisy.shape}                  ← y")   # y 的形状
print("-" * 64)                            # 分隔线
# 追踪一个边缘角点走完正向链，证明 img_clean 就是这么一步步算出来的
j, idx = 0, rows * cols - 1                # 选第 0 张图、最后一个角点（板角 = 大半径边缘角点 #53）
R0, _ = cv2.Rodrigues(rvecs[j])            # 第 0 张图的旋转向量 → 旋转矩阵
Xw = board[idx]                            # 该角点的板坐标 X_w
Xc = R0 @ Xw + tvecs[j]                    # 外参变到相机系：X_c = R·X_w + t
xn, yn = Xc[0] / Xc[2], Xc[1] / Xc[2]      # 归一化坐标 (x,y) = (X_c/Z_c, Y_c/Z_c)
print(f"[追踪] 第 {j} 张图 角点 #{idx}（板坐标 X_w = {Xw}）")           # 起点的板坐标
print(f"   外参变到相机系  X_c = ({Xc[0]:.1f}, {Xc[1]:.1f}, {Xc[2]:.1f})  [R·X_w + t]")  # 相机系坐标
print(f"   归一化          (x, y) = ({xn:.4f}, {yn:.4f})             [X_c/Z_c]")  # 归一化
print(f"   无噪声像素      u = ({img_clean[j, idx, 0]:.2f}, {img_clean[j, idx, 1]:.2f})  [经畸变+内参]")  # 像素=img_clean
print(f"   加噪观测        y = ({img_noisy[j, idx, 0]:.2f}, {img_noisy[j, idx, 1]:.2f})")  # 加噪后
print("-" * 64)                            # 分隔线
res = img_noisy - img_clean                # 残差 = 观测 − 应有 = 噪声本身
print(f"残差 (y - f) 的 std：u 轴 {res[..., 0].std():.4f}, v 轴 {res[..., 1].std():.4f}  （应 ≈ σ_a={sigma_a:.4f}）")  # 单轴 std 应回到 σ_a
print(f"径向残差 std = √(σ_a²+σ_a²) ≈ {np.sqrt((res**2).sum(-1).mean()):.4f}  （应 ≈ σ_det={sigma_det}）")  # 径向 std 应回到 σ_det
print(f"已保存 → {out}")                    # 保存路径
print("=" * 64)                            # 分隔线

运行结果解读。整板入画 50/50；100% 在画面内 表示真实约束满足，50 张图全部整块在画面内。[追踪] 角点 #53 跟一个角点走完正向链，证明 img_clean 就这么一步步算出来：板坐标 (200, 125, 0) → 相机系 (−3.7, −98.3, 514.0) → 归一化 (−0.0073, −0.1912) → 像素 (1264.1, 300.8)。残差 std ≈ 0.0708/0.0701 表示观测减应有等于噪声本身，其 std 应回到单轴 σ_a = 0.0707，径向 std 回到 σ_det = 0.1，这是自检：噪声加对了。

输出（下一步的输入）calib_sim.npz：

键	形状	含义	谁读它
`theta_true`	(9,)	真值内参 [fx, fy, cx, cy, k1, k2, p1, p2, k3]	第 2 步（雅可比求值点）、第 6 步（对照真值）
`board`	(54,3)	棋盘 3D 角点（Z=0，mm）	第 2 步（投影）、第 6 步（objectPoints）
`rvecs`,`tvecs`	(50,3)×2	每图外参（旋转、平移）	第 2 步（雅可比的外参列）
`img_clean`	(50,54,2)	无噪声投影 f(θ_true)	第 6 步（加噪基准）
`img_noisy`	(50,54,2)	加噪观测 y	第 6 步（喂给 calibrateCamera）
`sigma_a`	标量	单轴噪声 std	第 3 步（信息除 σ_a²）、第 6 步（加噪）
`K`,`dist`,`image_size`,`pattern`	—	内参矩阵、畸变、画面、棋盘规格	（备用）

CRLB 用真值算、不需要观测 y（img_noisy）；y 只在第 6 步验证时用。

第 2 步｜雅可比

做什么：算每个角点 (u, v) 对全部参数（9 内参 + 每图 6 外参 = 9 + 6×50 = 309）的偏导，用数值中心差分（扰动 ±h、重投影、差分），堆成 (5400, 309)；再用解析偏导 ∂u/∂k3 = fx·x·r⁶ 交叉验证。参数布局 θ = [fx, fy, cx, cy, k1, k2, p1, p2, k3, 每图 (r0, r1, r2, t0, t1, t2)]。

输入 → 输出：calib_sim.npz → J.npz（雅可比 + 参数名布局）。

"""
Step 2：算雅可比 J = ∂f/∂θ
==========================
读第 1 步的 calib_sim.npz（真值内参、板、位姿），用数值差分算完整雅可比：
每个角点的 (u,v) 对全部参数（9 内参 + 每图 6 外参）的偏导，堆成 (图数×角点×2, 9+6×图数)。
再用角点 #53 的解析偏导 ∂u/∂k3 做交叉验证。结果存 J.npz 给第 3 步。

输入： pipeline/data/calib_sim.npz   （第 1 步产出）
输出： pipeline/data/J.npz           （含 J、参数名与布局）
运行： python3 pipeline/step2_jacobian.py
"""
import os
import numpy as np
import cv2

HERE = os.path.dirname(os.path.abspath(__file__))
DATA = os.path.join(HERE, "data")

# 0. 读第 1 步数据
d = np.load(os.path.join(DATA, "calib_sim.npz"), allow_pickle=True)
theta = d["theta_true"]                  # 真值内参向量(9,) = [fx,fy,cx,cy,k1,k2,p1,p2,k3]
theta_keys = list(d["theta_keys"])       # 对应名称
board = d["board"]                       # 棋盘 3D 角点 (54,3)
rvecs = d["rvecs"]; tvecs = d["tvecs"]   # 外参 (n_img,3)
n_img = rvecs.shape[0]
n_pts = board.shape[0]                   # 54
print(f"读入：{n_img} 张图，每图 {n_pts} 角点，内参 {theta_keys}")

# 1. 参数向量布局
# θ = [9 个内参] + [每图 6 个外参(rvec3, tvec3)]，共 9 + 6*n_img 个参数。
intr_names = ["fx", "fy", "cx", "cy", "k1", "k2", "p1", "p2", "k3"]   # 前 9 个
extr_names = [f"img{i}_{t}" for i in range(n_img) for t in ("r0", "r1", "r2", "t0", "t1", "t2")]
param_names = intr_names + extr_names
n_intr = len(intr_names)                 # 9
n_param = len(param_names)               # 9 内参 + 6×图数 外参

# 拼出真值参数向量 params（数值差分在它附近扰动）
params = np.empty(n_param)
params[:n_intr] = theta                  # 前 9：内参
for i in range(n_img):                   # 后 6×图数：外参
    params[n_intr + 6*i:n_intr + 6*i + 3] = rvecs[i]      # 旋转向量
    params[n_intr + 6*i + 3:n_intr + 6*i + 6] = tvecs[i]  # 平移


def f_all(p):
    """参数向量 p(参数数,) → 全部角点像素 (n_img, n_pts, 2)。即 f(θ)。"""
    fx, fy, cx, cy = p[0], p[1], p[2], p[3]
    dist = p[4:9]                        # k1,k2,p1,p2,k3
    K = np.array([[fx, 0, cx], [0, fy, cy], [0, 0, 1.]])
    out = np.empty((n_img, n_pts, 2))
    for i in range(n_img):
        rv = p[n_intr + 6*i:n_intr + 6*i + 3]
        tv = p[n_intr + 6*i + 3:n_intr + 6*i + 6]
        pts, _ = cv2.projectPoints(board.reshape(-1, 1, 3),
                                   rv.reshape(3, 1), tv.reshape(3, 1), K, dist)
        out[i] = pts.reshape(-1, 2)
    return out


# 2. 数值雅可比（中心差分）
# J 行 = (图, 角点, u/v)，共 n_img×n_pts×2 行；列 = 9+6×n_img 个参数。
J = np.empty((n_img * n_pts * 2, n_param))
for k in range(n_param):
    h = 1e-6 * max(1.0, abs(params[k]))              # 自适应步长（大参数用大步长）
    pp = params.copy(); pp[k] += h
    pm = params.copy(); pm[k] -= h
    J[:, k] = (f_all(pp) - f_all(pm)).reshape(-1) / (2 * h)   # 中心差分 → 该参数那一列
print(f"雅可比 J 形状：{J.shape} = (图数×角点×2, 参数数)")

# 3. 行索引工具 + 交叉验证
def row(i, c, coord):
    """第 i 张图、第 c 个角点、coord(0=u,1=v) 在 J 中的行号。"""
    return (i * n_pts + c) * 2 + coord

k3_col = intr_names.index("k3")          # k3 在参数里的列号 = 8

# 角点 #53（板角）的 ∂u/∂k3，数值 vs 解析 fx·x·r⁶（证明雅可比算对了）
j_img, c_edge = 0, n_pts - 1             # 第 0 张图、最后一个角点 #53
num_du_dk3 = J[row(j_img, c_edge, 0), k3_col]

# 解析值：取该角点在第 0 张图的归一化坐标与 r，算 fx·x·r⁶
R0, _ = cv2.Rodrigues(rvecs[j_img])
Xc = R0 @ board[c_edge] + tvecs[j_img]
x, y = Xc[0] / Xc[2], Xc[1] / Xc[2]
r2 = x * x + y * y
ana_du_dk3 = theta[0] * x * r2**3        # fx · x · r⁶（r⁶=(r²)³）
print(f"[验证] 角点 #53 的 ∂u/∂k3：数值 {num_du_dk3:.4f}  vs  解析 fx·x·r⁶={ana_du_dk3:.4f}")

# 第 0 张图里【离主点最近】的角点，∂u/∂k3 应很小（半径小→畸变灵敏度低）
Xc_all = (R0 @ board.T).T + tvecs[j_img]              # 全部角点的相机坐标 (54,3)
rad = np.sqrt((Xc_all[:, 0] / Xc_all[:, 2])**2 + (Xc_all[:, 1] / Xc_all[:, 2])**2)  # 归一化半径
c_near = int(np.argmin(rad))
print(f"[验证] 第0张图离主点最近的角点 #{c_near}（r={rad[c_near]:.3f}）的 ∂u/∂k3："
      f"{J[row(j_img, c_near, 0), k3_col]:.3e}（半径小→畸变灵敏度低，应≈0）")

# 4. 保存
out = os.path.join(DATA, "J.npz")
np.savez(out, J=J, param_names=np.array(param_names),
         n_intr=n_intr, n_img=n_img, n_pts=n_pts)
print(f"已保存 → {out}")

运行结果解读。J 形状 (5400, 309)，5400 行 = 50 图 × 54 角 × 2(u/v)，309 列 = 9 内参 + 50×6 外参。角点 #53 ∂u/∂k3 数值 −0.0008 = 解析 −0.0008 表示数值雅可比与解析偏导完全一致，雅可比算对了（自检）。离主点最近角点 ∂u/∂k3 ≈ −2.8e-3，半径小则畸变灵敏度低，符合 fx·x·r⁶ 在小半径处趋零的预期。

输出（下一步的输入）J.npz：

键	形状	含义	谁读它
`J`	(5400,309)	雅可比 J = ∂f/∂θ	第 3 步（算 JᵀJ）
`param_names`	(309,)	参数布局名（前 9 内参、后外参）	第 4 步（按名取对角、分块）
`n_intr`	9	内参数（分块边界）	第 4 步（A/B/C 分块）
`n_img`,`n_pts`	50,54	图数、角点数	（备用）

深入：数值雅可比的数学原理（中心差分）

这一步用数值微分付算雅可比，不手推偏导。原理拆成五点。

雅可比是什么。f 把参数 θ（309）映射成像素（5400）。雅可比的第 (i, j) 元素是第 i 个观测对第 j 个参数的全偏导，整体是 5400×309 的矩阵。第 j 列就是"只动参数 j，所有像素怎么变"。

中心差分公式。导数定义取不了 h → 0，用对称扰动近似：

$$ \frac{\partial\mathbf f}{\partial\theta_k}\approx\frac{\mathbf f(\boldsymbol\theta+h\,\mathbf e_k)-\mathbf f(\boldsymbol\theta-h\,\mathbf e_k)}{2h} $$

代码对每个参数 k（雅可比的一列）扰动 +h 算一次 f、扰动 −h 算一次 f、差除 2h，得到该参数对所有观测的偏导（5400 个数 = 整列）。逐列构建：309 列 × 2 次 f_all = 618 次投影。

为什么中心差分不用前向，看泰勒展开。前向差分 [f(θ+h) − f(θ)]/h 的展开：

$$ f(\theta+h)=f(\theta)+h f'(\theta)+\frac{h^2}{2}f''(\theta)+\frac{h^3}{6}f'''(\theta)+\cdots $$ $$ \frac{f(\theta+h)-f(\theta)}{h}=f'(\theta)+\frac{h}{2}f''(\theta)+\cdots $$

首项误差正比于 h，是一阶精度 O(h)。中心差分 [f(θ+h) − f(θ−h)]/(2h)，把两边都展开：

$$ f(\theta+h)=f(\theta)+h f'+\frac{h^2}{2}f''+\frac{h^3}{6}f'''+\frac{h^4}{24}f^{(4)}+\cdots $$ $$ f(\theta-h)=f(\theta)-h f'+\frac{h^2}{2}f''-\frac{h^3}{6}f'''+\frac{h^4}{24}f^{(4)}-\cdots $$

相减时，偶数阶项（f, f’‘, f⁽⁴⁾, …）同号被抵消，奇数阶项（f’, f’‘’, …）异号翻倍：

$$ f(\theta+h)-f(\theta-h)=2h f'(\theta)+\frac{h^3}{3}f'''(\theta)+\cdots $$ $$ \frac{f(\theta+h)-f(\theta-h)}{2h}=f'(\theta)+\frac{h^2}{6}f'''(\theta)+\cdots $$

首项误差正比于 h²，是二阶精度。对称性让一阶误差抵消了。h = 10⁻⁶ 时前向约 10⁻⁶、中心约 10⁻¹²，差六个数量级。

步长 h = 10⁻⁶·max(1, |θ_k|) 由误差极小化定。两种误差在抢 h：截断误差（泰勒舍掉的项）正比于 h²，要 h 小；舍入误差（浮点精度约 2.2×10⁻¹⁶，两个接近的数相减丢有效位）正比于 ε/h，要 h 大。总误差 E(h) ≈ h² + ε/h，对 h 求导令零：

$$ \frac{dE}{dh}=2h-\frac{\varepsilon}{h^2}=0\;\Longrightarrow\;h^3=\frac{\varepsilon}{2}\;\Longrightarrow\;h\sim\varepsilon^{1/3}\approx6\times10^{-6} $$

取 10⁻⁶ 接近最优。max(1, |θ|) 自适应：大参数（fx = 2200）用 h = 2.2×10⁻³（相对扰动 10⁻⁶），小参数（k3 = 0）用 10⁻⁶。

数值与解析各有取舍。解析（手推 ∂u/∂fx 等）精确但只覆盖内参、要逐个推；数值万能（含外参、任意模型）、不记公式但慢。本步用数值、用解析（角点 #53 的 ∂u/∂k3）交叉验证。

第 3 步｜信息矩阵

做什么：把雅可比按 I = JᵀΣ⁻¹J 累加。各向同性独立噪声时 Σ = σ_a²·I，故 I = (1/σ_a²)·JᵀJ——这就是可加性：观测越多，偏导平方和越大，信息越多。

输入 → 输出：J.npz + calib_sim.npz（取 σ_a）→ I.npz（309×309）。

"""
Step 3：信息矩阵 I = Jᵀ Σ⁻¹ J
================================
读第 2 步的雅可比 J 和第 1 步的噪声 σ_a，累加成 Fisher 信息矩阵。
各向同性独立噪声时 Σ = σ_a²·I，故 I = (1/σ_a²)·JᵀJ。结果存 I.npz 给第 4 步。

输入： pipeline/data/J.npz            （第 2 步的雅可比）
       pipeline/data/calib_sim.npz    （取 σ_a）
输出： pipeline/data/I.npz            （信息矩阵）
运行： python3 pipeline/step3_fisher.py
"""
import os
import numpy as np

HERE = os.path.dirname(os.path.abspath(__file__))
DATA = os.path.join(HERE, "data")

# 0. 读输入
Jd = np.load(os.path.join(DATA, "J.npz"), allow_pickle=True)
J = Jd["J"]                              # 雅可比 (n_img×n_pts×2, 参数数)
param_names = list(Jd["param_names"])
n_intr = int(Jd["n_intr"])

cd = np.load(os.path.join(DATA, "calib_sim.npz"), allow_pickle=True)
sigma_a = float(cd["sigma_a"])          # 单轴噪声 std
print(f"读入：J {J.shape}，σ_a={sigma_a:.6f}（σ_a²={sigma_a**2:.6f}）")

# 1. 信息矩阵 I = (1/σ_a²) JᵀJ
# JᵀJ[i,j] = Σ_m J[m,i]·J[m,j]：对角(i=j)=参数i的偏导平方和(自信息)，非对角(i≠j)=i,j的交叉积(耦合)；除 σ_a² 折入噪声刻度。
JtJ = J.T @ J                            # 无噪声刻度的信息（JᵀJ）
I = JtJ / (sigma_a ** 2)                 # Fisher 信息矩阵
print(f"信息矩阵 I 形状：{I.shape}")

# 2. 验证关键对角元
def diag(name):
    return I[param_names.index(name), param_names.index(name)]

# 验证关键对角元（值随配置变，这里只打印、不硬编码期望）
print(f"I[k3,k3] = {diag('k3'):.4e}")
print(f"I[fx,fx] = {diag('fx'):.1f}  (未边缘化)")

# 3. 保存
out = os.path.join(DATA, "I.npz")
np.savez(out, I=I, param_names=np.array(param_names), n_intr=n_intr, sigma_a=sigma_a)
print(f"已保存 → {out}")

运行结果解读。I 形状 (309, 309) 是方阵，行列都是 309 个参数。I[k3,k3] = 1.54e7 是 k3 的信息：所有角点偏导平方和除以 σ_a²。I[fx,fx] = 25362 是 fx 的信息，但这是未边缘化外参时的表观值，第 4 步会缩水。

JᵀJ 怎么把每个参数的"所有观测偏导平方和"累起来，看一行。JᵀJ 的第 (i, j) 元素是对所有观测 m 求和 J[m,i]·J[m,j]。对角元（i = j）是参数 i 各偏导的平方和，即该参数的自信息；非对角元（i ≠ j）是参数 i、j 偏导的交叉积，刻画两个参数的耦合。整个矩阵除以 σ_a² 折入噪声刻度。

输出（下一步的输入）I.npz：

键	形状	含义	谁读它
`I`	(309,309)	信息矩阵 I = (1/σ_a²)·JᵀJ	第 4 步（分块边缘化）
`param_names`	(309,)	参数布局	第 4 步（分块边界 n_intr）
`n_intr`	9	内参数（前 9 行列是内参）	第 4 步（A/B/C 分块）

第 4 步｜边缘化外参

做什么：外参（nuisance）也从同一批数据估、带不确定度，会"吃掉"内参表观信息。按内参、外参分块，Schur 补扣掉外参，得 9×9 内参信息：

$$ \mathcal I=\begin{pmatrix}\mathbf A&\mathbf B\\ \mathbf B^\top&\mathbf C\end{pmatrix},\qquad \mathcal I_\beta=\mathbf A-\mathbf B\mathbf C^{-1}\mathbf B^\top $$

其中 C⁻¹Bᵀ 用 np.linalg.solve 解方程，比直接求逆稳。

输入 → 输出：I.npz → Ibeta.npz（9×9）。

"""
Step 4：边缘化外参 → 内参信息矩阵 I_β
======================================
读第 3 步的信息矩阵 I，按 内参(9)/外参(6×图数) 分块，用 Schur 补把外参这个
nuisance 参数的不确定性扣掉，得到只关于内参的 9×9 信息矩阵 I_β。
为什么必须做：外参也是从同一批数据估的、带不确定度，会"吃掉"内参的表观信息。
结果存 Ibeta.npz 给第 5 步。

输入： pipeline/data/I.npz   （第 3 步的信息矩阵）
输出： pipeline/data/Ibeta.npz （9×9 内参信息矩阵）
运行： python3 pipeline/step4_marginal.py
"""
import os
import numpy as np

HERE = os.path.dirname(os.path.abspath(__file__))
DATA = os.path.join(HERE, "data")

# 0. 读输入
Id = np.load(os.path.join(DATA, "I.npz"), allow_pickle=True)
I = Id["I"]                              # 信息矩阵
param_names = list(Id["param_names"])
n_intr = int(Id["n_intr"])               # 9
intr_names = param_names[:n_intr]
print(f"读入：I {I.shape}，内参数 {n_intr}（{intr_names}）")

# 1. 分块
# I = [[A, B], [Bᵀ, C]]：A=内参×内参，C=外参×外参，B=内参×外参的耦合。
A = I[:n_intr, :n_intr]                  # (9,9)
B = I[:n_intr, n_intr:]                  # 内参×外参（耦合）
C = I[n_intr:, n_intr:]                  # 外参×外参

# 2. Schur 补：I_β = A - B·C⁻¹·Bᵀ
# C⁻¹·Bᵀ 用解线性方程（比直接求逆更稳）：解 C·X = Bᵀ，得 X = C⁻¹Bᵀ
CinvBT = np.linalg.solve(C, B.T)         # = C⁻¹ Bᵀ
Ibeta = A - B @ CinvBT                   # (9,9) 内参信息矩阵（扣掉外参不确定性后）
print(f"内参信息矩阵 I_β 形状：{Ibeta.shape}")

# 3. 验证：fx 信息边缘化前后对比
fx = intr_names.index("fx")
before = I[fx, fx]                       # 未边缘化时 fx 的信息（来自全 I 的对角）
after = Ibeta[fx, fx]                    # 边缘化后
print(f"I_β[fx,fx]：{before:.0f} → {after:.0f}（缩 {before/after:.1f} 倍）")

# 4. 保存
out = os.path.join(DATA, "Ibeta.npz")
np.savez(out, Ibeta=Ibeta, intr_names=np.array(intr_names))
print(f"已保存 → {out}")

运行结果解读。I_β 形状 (9, 9) 只剩内参的信息矩阵，外参 nuisance 已扣掉。I_β[fx,fx] 25362 → 532（缩 47.7 倍） 表示边缘化掉 50 组外参后，fx 的表观信息被外参不确定性吃掉绝大部分；不边缘化、假装外参已知，CRLB 会假性乐观约 √48 ≈ 7 倍。

I_β 是个对称阵、且有负数，正常吗。正常。对称来自 JᵀΣ⁻¹J 的构造（J 的列两两内积，可交换）；非对角元为负表示该对内参负相关（如 fx 与 cx：沿对角线方向的偏导一正一负，交叉积为负）；矩阵整体仍正定，对任意非零向量 v 都有 vᵀI_βv > 0，这是 CRLB 能求逆的前提。负的耦合元不意味着信息为负，只表示参数间互相牵制。

输出（下一步的输入）Ibeta.npz：

键	形状	含义	谁读它
`Ibeta`	(9,9)	内参信息矩阵（边缘化外参后）I_β	第 5 步（求逆开方得 CRLB）
`intr_names`	(9,)	内参名 [fx, fy, cx, cy, k1, k2, p1, p2, k3]	第 5 步（按名打印/对齐）

深入:`np.linalg.solve(C, B.T)` 的作用与不用 `inv` 的原因

代码里 np.linalg.solve(C, B.T) 解的是线性方程组 C·X = Bᵀ，求 X。

C（300×300）是外参信息块，Bᵀ（300×9）是耦合块转置。结果 X = C⁻¹Bᵀ（300×9），然后 B·X = BC⁻¹Bᵀ（9×9），就是从 A 扣掉的渗透量。

两种写法数学等价，但 solve 更好：

	`solve(C, B.T)`	`inv(C) @ B.T`
内部	LU 分解 + 前代/回代，不显式构造 C⁻¹	先求完整 300×300 逆矩阵，再乘
稳定性	更稳（带部分选主元，避免逆放大舍入）	C 病态时显式逆误差大
速度	快（只解方程）	慢（多算完整逆）

类比解 3x = 6。inv 法先算 1/3 = 0.333… 再乘 6（多一步中间量，精度损失）；solve 法直接 6/3 = 2（不经过 0.333…，更准）。能解方程就别显式求逆，是数值线性代数的标准实践。

第 5 步｜CRLB

做什么：对 I_β 求逆、取对角、开方，得到每个内参的精度下界：

$$ \mathrm{CRLB}_k=\sqrt{[\mathcal I_\beta^{-1}]_{kk}} $$

输入 → 输出：Ibeta.npz → crlb.npz。

"""
Step 5：CRLB = √diag(inv(I_β))
===============================
读第 4 步的内参信息矩阵 I_β，求逆、取对角、开方，得到每个内参的精度下界
（Cramér–Rao 下界）：任何无偏估计算法，该参数的 1σ 误差都不可能小于它。
结果存 crlb.npz 给第 6 步做对照。

输入： pipeline/data/Ibeta.npz   （第 4 步的 9×9 内参信息矩阵）
输出： pipeline/data/crlb.npz    （每个内参的精度下界）
运行： python3 pipeline/step5_crlb.py
"""
import os
import numpy as np

HERE = os.path.dirname(os.path.abspath(__file__))
DATA = os.path.join(HERE, "data")

# 0. 读输入
Bd = np.load(os.path.join(DATA, "Ibeta.npz"), allow_pickle=True)
Ibeta = Bd["Ibeta"]                      # (9,9)
intr_names = list(Bd["intr_names"])
print(f"读入：I_β {Ibeta.shape}（内参 {intr_names}）")

# 1. 求精度下界
# Cov(θ̂) ≥ I_β⁻¹；第 k 个参数的方差下界 = [I_β⁻¹]_kk，std 下界 = √之。
cov_lb = np.linalg.inv(Ibeta)            # 内参协方差矩阵的下界 (9,9)
crlb = np.sqrt(np.diag(cov_lb))          # 每个内参的 std 下界 (9,)

# 2. 打印 + 验证
print("内参精度下界 CRLB（1σ，像素 / 无量纲）：")
for name, v in zip(intr_names, crlb):
    print(f"   {name:>3} : {v:.4g}")

fx = intr_names.index("fx")
print(f"\nCRLB_fx = {crlb[fx]:.4f} px")
# 关键：CRLB 来自完整逆，不是对角元倒数（因内参强耦合）
diag_recip = 1.0 / np.sqrt(Ibeta[fx, fx])
print(f"1/√I_β[fx,fx] = {diag_recip:.4f}  ← 错（假性偏小），CRLB 须用完整逆 → {crlb[fx]:.4f}")

# 3. 保存
out = os.path.join(DATA, "crlb.npz")
np.savez(out, crlb=crlb, cov_lb=cov_lb, intr_names=np.array(intr_names))
print(f"\n已保存 → {out}")

运行结果：

1
2
3

  fx: 0.1542   fy: 0.1579   cx: 0.2548   cy: 0.1998
  k1: 3.25e-4  k2: 1.91e-3  p1: 2.20e-5  p2: 2.87e-5  k3: 3.49e-3
CRLB_fx = 0.154 px；1/√I_β[fx,fx] = 0.0434（错，须用完整逆）

为什么是对角线之逆、不是逆对角线，这是最容易搞错的地方。直觉上"信息越大精度越高"，让人以为 CRLB_k = 1/√I_β[k,k]。用 2×2 的例子看清：

$$ \mathcal I_\beta=\begin{pmatrix}a&b\\b&c\end{pmatrix},\qquad \mathcal I_\beta^{-1}=\frac{1}{ac-b^2}\begin{pmatrix}c&-b\\-b&a\end{pmatrix} $$

第 (1,1) 元素是 c/(ac − b²)，而不是 1/a。只有 b = 0（无耦合、对角阵）时才等于 1/a。耦合 b ≠ 0 时通过行列式 ac − b² 改变了每个参数的方差下界。标定实例：1/√I_β[fx,fx] = 0.043px（错，假性偏小），√[I_β⁻¹][fx,fx] = 0.154px（对，完整逆），差 3.5 倍。

深入:条件数与数值稳定(求逆出现 NaN 的原因)

求 CRLB 要对 I_β 求逆。条件数衡量矩阵有多病态：

$$ \mathrm{cond}(\mathcal I_\beta)=\frac{\lambda_{\max}}{\lambda_{\min}}\quad(\text{最大/最小特征值；np.linalg.cond 默认即此}) $$

对 Fisher 信息，λ_max 是最易辨识方向的信息（如畸变 k3），λ_min 是最难辨识方向的信息（如 fx、cy 耦合方向）。

为什么条件数大会让求逆出 NaN。求逆时 λ_min 方向的方差是 1/λ_min（很大），而浮点精度 ε ≈ 2.2×10⁻¹⁶，相对 λ_max 的舍入"地板"是 ε·λ_max。一旦 λ_min < ε·λ_max，λ_min 落在地板之下，算出来可能接近零甚至为负，求逆爆炸，再开方就是 NaN。

配置	λ_max	λ_min	cond	ε·λ_max	结果
旧（角点出画面、参数锁死联动）	~10¹⁴	~10⁻⁴	~10¹⁸	~2×10⁻²	NaN
本例（整板入画、独立多样性）	8.1×10⁹	15.3	5.3×10⁸	~1.8×10⁻⁶	正常

本例三件事让条件数可控：整板入画压低 λ_max（不再有出画面的大半径角点把畸变撑到 10¹⁴）；参数独立多样性打破 fx、cy 耦合退化、抬高 λ_min；多图累加各方向信息。

运行结果解读。每行是一个内参的 CRLB（1σ 精度下界，像素/无量纲）：焦距 fx 0.154、fy 0.158，主点 cx 0.255、cy 0.200，畸变系数更小（k1 3.2e-4、k3 3.5e-3），畸变估得比焦距更准。CRLB_fx = 0.154 表示任何算法估 fx 的 1σ 误差都不可能小于 0.154px，这是精度天花板。1/√I_β[fx,fx] = 0.0434（错） 是警示：CRLB 来自完整逆、不是对角元倒数，差 3.5 倍，因内参强耦合。

输出（下一步的输入）crlb.npz：

键	形状	含义	谁读它
`crlb`	(9,)	每个内参精度下界 √[I_β⁻¹]_kk（金标准，真值处算）	第 6 步（对照实测 std）
`cov_lb`	(9,9)	内参协方差下界 I_β⁻¹	（备用：完整协方差）
`intr_names`	(9,)	内参名	第 6 步（按名对齐）

第 6 步｜Monte Carlo 验证

做什么：检验真实标定算法（cv2.calibrateCamera，高斯噪声下等价 MLE）的散布是否真达到第 5 步的金标准 CRLB。固定几何，换噪声种子重复 N 次——每次加新噪声、跑标定、收下内参；算每个内参跨多次的 std，除以 CRLB，应接近 1。这正是"仿真有真值"的价值：用真值处的 CRLB 验证估计器（也即验证真实标定里要用的观测信息法）靠谱。

输入 → 输出：calib_sim.npz（取板、f、σ_a）+ crlb.npz → mc_stats.npz。

"""
Step 6：Monte Carlo 验证 std / CRLB ≈ 1
=========================================
前 5 步在真值上算出了精度下界 CRLB。这一步反过来检验：真实标定算法
（cv2.calibrateCamera，高斯噪声下等价于 MLE）的散布是否真的达到这个下界。
做法：固定几何（板、位姿），换不同噪声种子重复 N 次——每次给角点加新噪声、
跑一遍 calibrateCamera，收下估出的内参；最后算每个内参跨多次估计的 std，
除以 CRLB。≈1 → 估计器最优（达到下界）；≫1 → 低效或退化。

输入： pipeline/data/calib_sim.npz   （取板、位姿、img_clean=f、σ_a）
       pipeline/data/crlb.npz        （第 5 步的精度下界）
输出： pipeline/data/mc_stats.npz    （每个内参的 std、std/CRLB）
运行： python3 pipeline/step6_validate.py
"""
import os
import numpy as np
import cv2

HERE = os.path.dirname(os.path.abspath(__file__))
DATA = os.path.join(HERE, "data")
N_RUNS = 60                               # Monte Carlo 次数（越多 std 估得越稳，但更慢）

# 0. 读输入
cd = np.load(os.path.join(DATA, "calib_sim.npz"), allow_pickle=True)
board = cd["board"].astype(np.float64)            # (54,3)
img_clean = cd["img_clean"]                       # (n_img,54,2) = f(θ_true)
sigma_a = float(cd["sigma_a"])
W, H = int(cd["image_size"][0]), int(cd["image_size"][1])
n_img = img_clean.shape[0]
n_pts = board.shape[0]

Cd = np.load(os.path.join(DATA, "crlb.npz"), allow_pickle=True)
crlb = Cd["crlb"]                                 # (9,) 第 5 步的精度下界
intr_names = list(Cd["intr_names"])
print(f"读入：{n_img} 图，每图 {n_pts} 角点，σ_a={sigma_a:.4f}，CRLB 已就绪。Monte Carlo {N_RUNS} 次。")

# cv2 要的输入格式：每图一份 objectPoints(54,1,3 float32) 与 imagePoints(54,1,2 float32)
obj = [np.ascontiguousarray(board.astype(np.float32).reshape(-1, 1, 3)) for _ in range(n_img)]

# 1. Monte Carlo：换噪声种子重复标定
est = []                                  # 收每次估计出的 9 个内参 [fx,fy,cx,cy,k1,k2,p1,p2,k3]
for s in range(N_RUNS):
    rng = np.random.default_rng(1000 + s)             # 每次换种子 → 新噪声
    noisy = img_clean + rng.normal(0, sigma_a, img_clean.shape)
    img = [np.ascontiguousarray(noisy[i].astype(np.float32).reshape(-1, 1, 2)) for i in range(n_img)]
    K0 = np.zeros((3, 3)); K0[2, 2] = 1.0; dist0 = np.zeros(5)
    # 不给内参先验（纯标定），让 cv2 自己初估；估 5 系数畸变（与生成模型一致）
    ok, K_est, dist_est, r_, t_ = cv2.calibrateCamera(obj, img, (W, H), K0, dist0,
                                                      flags=0)
    if not ok:
        continue
    p = [K_est[0, 0], K_est[1, 1], K_est[0, 2], K_est[1, 2]] + list(dist_est[:5])
    est.append(p)
est = np.array(est)                       # (n_ok, 9)
print(f"成功 {est.shape[0]}/{N_RUNS} 次。")

# 2. 实测统计：bias（估计均值−真值，判无偏）+ std vs CRLB（判有效）
mean_est = est.mean(axis=0)                # 多次估计的均值
theta_true = cd["theta_true"]
bias = mean_est - theta_true               # bias = 均值 − 真值（应≈0）
std = est.std(axis=0)                      # 跨多次估计的 std
ratio = std / crlb                         # std/CRLB（应≈1 → 有效）
bias_ratio = bias / crlb                   # bias/CRLB（应≪1 → 无偏）
print("\n内参   真值        估计均值      bias        bias/CRLB   实测std    CRLB      std/CRLB")
for k, name in enumerate(intr_names):
    print(f"  {name:>3}  {theta_true[k]:>10.4g}  {mean_est[k]:>10.4g}  {bias[k]:>+10.3e}  {bias_ratio[k]:>8.2f}   {std[k]:.4g}   {crlb[k]:.4g}   {ratio[k]:.2f}")
print(f"\nbias/CRLB 最大 = {np.abs(bias_ratio).max():.2f}（应≪1 → 无偏）；std/CRLB 均值 ≈ {ratio.mean():.2f}（应≈1 → 有效）")

# 3. 保存
out = os.path.join(DATA, "mc_stats.npz")
np.savez(out, est=est, mean_est=mean_est, bias=bias, std=std, crlb=crlb,
         ratio=ratio, bias_ratio=bias_ratio, intr_names=np.array(intr_names))
print(f"已保存 → {out}")

运行结果（完整输出）：

成功 60/60 次。

内参   真值        估计均值      bias        bias/CRLB   实测std    CRLB      std/CRLB
   fx        2200        2200  +6.400e-03      0.04   0.1255   0.1542   0.81
   fy        2200        2200  +1.352e-02      0.09   0.1287   0.1579   0.81
   cx        1280        1280  -3.122e-02     -0.12   0.2539   0.2548   1.00
   cy         720         720  -1.462e-02     -0.07   0.1989   0.1998   1.00
   k1        -0.1        -0.1  -1.098e-05     -0.03   0.0003282   0.0003248   1.01
   k2           0   4.603e-05  +4.603e-05      0.02   0.001933   0.001905   1.01
   p1           0  -8.597e-07  -8.597e-07     -0.04   2.345e-05   2.203e-05   1.06
   p2           0  -2.517e-06  -2.517e-06     -0.09   2.489e-05   2.866e-05   0.87
   k3           0  -2.172e-05  -2.172e-05     -0.01   0.003633   0.00349   1.04

bias/CRLB 最大 = 0.12（应≪1 → 无偏）；std/CRLB 均值 ≈ 0.96（应≈1 → 有效）

逐列含义：

列	含义	作用
真值	仿真里设的真值（上帝知道、估计器不知道）	对照基准
估计均值	60 次标定结果取平均	应接近真值
bias	估计均值 − 真值	应接近 0，正负表示偏高偏低
bias/CRLB	bias 除以 CRLB	应远小于 1（偏差远小于随机散布，判无偏）
实测 std	60 次估计的散布（标准差）	应接近 CRLB
CRLB	第 5 步算出的金标准下界	理论最优 std
std/CRLB	实测散布 ÷ 理论下界	应接近 1（达到下界，判有效）

fx、fy（焦距）bias 极小（+0.006、+0.014px，不到散布的 4%–9%），无偏；std/CRLB = 0.81，实测散布略小于理论下界，可能因 60 次样本估 std 有轻微下偏，或高斯下界稍保守（渐近成立、有限样本有波动）。cx、cy（主点）bias/CRLB = −0.12、−0.07（最大的一项绝对值 0.12），偏差仍远小于散布（bias 0.03px 对 CRLB 0.25px），可接受；std/CRLB = 1.00，主点估计恰好达到下界。k1（一阶径向畸变）bias/CRLB = −0.03、std/CRLB = 1.01，几乎完美。k2、p1、p2、k3（高阶畸变，真值 = 0）bias/CRLB 绝对值都不超过 0.09，无偏；std/CRLB 在 0.87–1.06，有效。

总结判据：

判据	公式	标准	本例结果	结论
无偏性		bias/CRLB		远小于 1
有效性	std/CRLB	接近 1	均值 0.96	有效

无偏加有效，是最优估计器。这验证了两件事：整条 Fisher 信息到雅可比到信息到边缘化到 CRLB 的流水线，推导和代码都对；真实标定里用的观测信息法（在估计处算、而非真值处）也能给出正确的不确定度——仿真替你验证过了。

输出 mc_stats.npz：

键	形状	含义
`est`	(60,9)	每次 Monte Carlo 估出的 9 个内参
`mean_est`	(9,)	跨 60 次的估计均值
`bias`	(9,)	bias = 估计均值 − 真值（应接近 0）
`bias_ratio`	(9,)	bias/CRLB（应远小于 1，判无偏）
`std`	(9,)	跨 60 次的实测 std
`crlb`	(9,)	金标准（复制自第 5 步）
`ratio`	(9,)	std/crlb（应接近 1，判有效）

结果一览

结果	数值	说明什么	出处
fx 的 CRLB	0.154 px	任何算法估 fx 的 1σ 下限，精度天花板，拍照前可算（仿真里）；真实标定用观测信息近似此值	CRLB
边缘化缩减	25362 → 532（48 倍）	外参 nuisance 吃掉 fx 表观信息；不边缘化 CRLB 假性乐观约 7 倍	Schur 补
CRLB ≠ 对角倒数	0.154 ≠ 0.043	内参强耦合，精度来自完整逆	CRLB
条件数	cond(I_β) = 5.3×10⁸	整板入画加多样性把 cond 压到可控；过大则求逆 NaN	第 5 步
std / CRLB ≈ 1	Monte Carlo 实测	估计器达下界（渐近有效），验证观测信息法可靠	MLE 渐近

踩过的坑

仿真用真值、真实用估计。CRLB 在真值处算（期望信息，金标准，仅仿真可得）；真实标定在估计 θ̂ 处算观测信息、取其逆作不确定度（渐近等价）。

整板必须入画。真实标定只看得见画面内角点，部分出画会虚高精度。

位姿参数要独立多样。锁死联动会让 fx、cy 退化、条件数爆、CRLB 出 NaN。用混合基数独立遍历。

σ_a = σ_det/√2。公式里用 σ_a²，否则偏 √2 倍。径向误差平方是两个轴向平方和，故 σ_det² = 2σ_a²。

外参必须边缘化。Schur 掉 nuisance，否则内参 CRLB 假性偏小。

CRLB 不是对角元倒数。内参强耦合，来自完整逆，差可达数倍。

条件数大时 CRLB 不可信。λ_min < ε·λ_max 会让求逆出 NaN；靠入画（压 λ_max）、多样性（抬 λ_min）、多图控制。

cv2 要 float32。calibrateCamera 的点须 float32 加 ascontiguousarray，float64 会出错。

σ_det 要实测。静态重复法：相机加板静止采 100 帧，取角点径向 RMS。

不能只拍正面的原因

标定时常听到的建议"板要倾斜着拍"不是玄学——只拍正面（板正对相机，pitch、yaw≈0）会让焦距、主点、畸变几乎定不准。做个对照：把第 1 步的 pitch、yaw 从 ±30° 改成 0，其余不变，重跑整条流水线。

指标	基线（±30°）	只拍正面（=0）	放大
fx 的 CRLB	0.154 px	4.2×10⁵ px	2.7×10⁶
cx 的 CRLB	0.255 px	1.0×10⁵ px	4.1×10⁵
k1 的 CRLB	3.2×10⁻⁴	38	1.2×10⁵
p2 的 CRLB	2.9×10⁻⁵	4.75	1.7×10⁵
I_β 的条件数	5.3×10⁸	∞（矩阵不定）	—

正面拍的退化极其彻底：边缘化后的内参信息矩阵 I_β 出现负特征值（最小特征值从 15 掉到小于 0），条件数直接 ∞，求逆失去意义。焦距的 CRLB 从 0.15px 飙到 4×10⁵px，差两百多万倍；径向畸变 k1 从 3×10⁻⁴ 飙到 38，切向 p2 从 3×10⁻⁵ 飙到 4.75。表中没列的 cy、p2 显示为 0 也不是精度高，而是矩阵已经不定、方差算出负值被截断了——纯属数值崩溃。

原因还是信息矩阵：正面拍时板的透视退化为相似变换，几组参数对像素的影响变得共线，信息矩阵出现近零特征方向——

焦距 f ↔ 深度 t_z（尺度模糊）：正面板的表观大小 = f·板宽/Z，只能定 f/Z 这个比值，f 和 Z 可同比例缩放。
主点 c_x/c_y 与板位置、焦距混在一起。
畸变：径向要大半径角点（正面居中板半径小，测不到），切向 p1/p2 本质描述镜头偏心，几乎必须有倾斜或不对称才能定。

倾斜板带来透视压缩（分开 f 和距离）、把角点扫过大半径并制造不对称（分开畸变），把这些方向重新撑开。第 1 步采样里 pitch/yaw ±30°、roll ±45°、九宫格位置、不同 fill 的混合遍历，正是为此——这也是为什么"位姿参数要独立多样"那条坑会致命。

退化的发现与定位

退化数据最阴险的地方在于：重投影误差往往很低。欠定的参数能把这批退化数据拟合得很好（拟合好 ≠ 参数对），正面拍 50 张图标完，RMS 可能还是 0.1px，看着很美，但 fx 其实可以是 2200 也可以是 40 万。所以 RMS 是必要不充分条件，绝不能只靠它判好坏。

发现方法(从省事到深入)

看 stdDeviations（OpenCV 白给）。calibrateCamera 返回的 stdDeviations 就是 √diag(I_obs⁻¹)，即每个参数的不确定度。fx 的 stdDev 是 0.15px → 正常；是几万 px → 直接报警。比 RMS 靠谱得多。
子集复现性（不用任何理论）。把图分两半各标定一次，比内参；两次 fx 差很多 → 欠定。leave-one-out（每次抽掉一张重标）同理，fx 来回跳 → 约束脆弱。这是没真值也能查的穷人验证。
条件数 cond(I_β)。自己建信息矩阵（本流水线），cond 飙到 1e14+ 或 ∞ → 退化，直接刻画病态程度。
留出图测试。用正面图标定，去重投影一张倾斜的图，误差爆 → 模型根本没学到真畸变/真焦距。

定位步骤:先审输入,再特征分解

最该先做的是审输入位姿分布：每张图的板法向（从 rvec 算）、板在画面的位置、板大小，各自方差≈0 的那一维就是退化维。正面拍的法向全指向相机，一眼就出。

更精确的是特征分解 I_β。先把两个角色分清：

特征值（标量 λ）：这条方向有多不可辨识，越小越退化，≈0 就是定不出来。
特征向量（9 维向量）：活在参数空间里，第 k 个分量对应第 k 个参数（第 1/2 位 = fx/fy，第 3/4 位 = cx/cy，后 5 位 = 畸变）。读哪几个分量绝对值大，就是哪几个参数。

把"为什么"说透——信息矩阵是对数似然在最优点处的曲率，特征值就是沿该方向的曲率：

$$ \mathcal I_\beta=\mathbf V\boldsymbol\Lambda\mathbf V^\top,\qquad [\mathcal I_\beta^{-1}]_{kk}=\sum_j \frac{V_{kj}^{2}}{\lambda_j} $$

λ≈0 的方向似然是平的（扰动参数、像素几乎不变）→ 不可辨识；上面第二个式子说明，任何近零 λ_j 都会让 1/λ_j 爆炸，把所有在 v_j 上有分量的参数的方差下界一起抬高——这就是退化会"连累"一堆参数 CRLB 飙升的根因。

用前面「只拍正面」那组退化数据，I_β 有 3 个近零特征值，特征向量这样读：

近零 λ	大分量	读作
3.3×10⁻⁹	cy=−0.88, cx=−0.48	主点
1.9×10⁻⁷	cx=+0.87, cy=−0.47	主点（正交组合）
3.6×10⁻⁷	fx=+0.70, fy=+0.70	焦距

两条主点 + 一条焦距，即"正面拍 → 主点和焦距不可辨识"，和直觉完全吻合。分量间的符号比值也有含义：fx、fy 同号等大表示二者同步缩放（典型的 f/距离尺度退化）；cx、cy 的组合给出主点沿某朝向移动。

定位后的修复方法

特征向量告诉你缺哪一维多样性，对症补采集：

主点方向退化 → 多拍板偏离中心（含画面角落）的图；
焦距方向退化 → 多拍不同倾斜、不同距离的图；
畸变方向退化 → 多拍大半径（板占满画面）+ 倾斜的图。

补齐后 cond(I_β)、stdDeviations 会一起回落，就是修好了的信号。

采集前的准备

正经采集前，核心是把会改变内参的变量锁死、系统误差源提前排除——而且可以用仿真在拍之前就预估精度够不够。这份准备清单决定了整轮采集能不能用。

锁死相机的几何相关设置

内参是"这台相机在这个设置下"的属性，设置一变就作废：

手动对焦并锁定：标定中绝不能让自动对焦跑——焦距变了，每张图的内参都不同，标定无意义。在工作距离对实。
固定焦距（变焦别动 zoom），定焦最好。
固定光圈（中段 f/5.6~f/8 通常最锐、像差最小）。
关掉机内畸变/镜头校正（数字变焦、畸变校正 profile 都关），否则标定的是校正后的图，不是镜头本身。
用实际处理的分辨率拍：fx/fy/cx/cy 是像素单位、和分辨率绑定，别满分辨率标定再降采样。
固定 ISO/曝光选低噪声档，相机热稳定后再开拍。

标定板

平：翘曲是 CRLB 看不见的系统偏差。用铝板/玻璃背衬的刚性板，别用软纸或易翘覆膜，拿直尺验。
方格尺寸量准：数显卡尺量多个方格取平均。它决定物理尺度（fx 像素值不受影响，但外参和物理解读依赖它）。
大小匹配 FOV/工作距离：板要能在工作距离占画面约 1/3~1/2；广角要大板，长焦要小板。
清洁、无反光、黑白高对比。

环境与拍摄

匀光漫射光，避免反光和强阴影；防抖（三脚架或高速快门，运动模糊直接抬高 σ_a）；每个曝光内板和相机静止。

先测噪声底，再预估精度（最值钱的一步）

整套理论最实在的落地——拍之前就能知道够不够：

静态重复测 σ_a：相机和板都不动，连拍约 100 帧，检测角点，算角点位置跨帧的 std → 得径向 σ_det，单轴 σ_a = σ_det/√2。这就是你的噪声底。
仿真预估 CRLB：把计划好的板、规划的位姿覆盖、测到的 σ_a 喂进本文流水线（第 1~5 步），算出预期 CRLB。满足下游精度 → 放心拍；不满足 → 先调（加图、加多样性、降噪声）再拍，省下一整轮废采集。

这一步把"实验设计"变成可计算的：拍前估、拍后验，是 Fisher 信息的核心用途，也是本文仿真流水线最实在的价值。

规划覆盖、流程就绪、记元数据

列位姿清单，确保倾斜全方向、位置含角落、远近多样、滚转都有（对应那张"多样性↔参数方向"表）。
管线先测通：先拍 1~2 张跑检测+标定，确认角点检测稳、代码能跑、现场能即时查 RMS/stdDeviations。
记元数据：相机型号、镜头、标称焦距、光圈、对焦、分辨率、工作距离、板型号与实测方格尺寸、光照。

可靠标定的标准流程

前面用仿真验证了方法、用退化实验看清了坑。把结论收拢成一条可照着做的标准流程：多样的采集 + 恰当的模型 + 自检诊断闭环，目标是让 I_β 良态、stdDeviations 小到够用、且无系统偏差。Fisher 信息把每一步都变成可量化、可验证的事。

采集：多样覆盖，让 I_β 良态

这是决定成败、也最难补救的一步。本质是在参数空间里给每个特征方向喂够信息——信息是各图累加、并按噪声缩放的：

$$ \mathcal I=\sum_{\text{图}}\frac{1}{\sigma_a^{2}}\,\mathbf J^\top\mathbf J $$

所以"图越多越好"只在新方向上成立（雷同图在已覆盖方向上重复累加，是浪费），而 σ_a 越小（图越清晰）每个角点贡献越大。每维多样性对应它撑起的参数方向：

采集维度	撑起的参数方向	缺了会怎样
板多向倾斜（pitch/yaw 正负都覆盖）	焦距（破 f↔Z 尺度）、主点、切向畸变	焦距/主点退化（实验里 ±30° 够，=0 崩）
板偏离中心、拍到画面角落	主点、大半径径向畸变	主点和 k1/k2/k3 定不准
不同距离（板占画面 1/3~1/2）	分开焦距与畸变（角点扫不同半径）	焦距畸变耦合
板自身滚转	解耦面内参数	部分参数强相关
覆盖全视场（尤其边缘）	径向畸变（要大半径角点）	高阶畸变测不到

清单：倾斜全方向、位置铺满含角落、远近多样、滚转有；图清晰对焦光照足；板要平、方格尺寸量准（翘曲或错尺寸是 CRLB 看不见的系统偏差）。数量服从多样性——20 张充分错开胜过 100 张雷同正面；实际 15–30 张高质量全覆盖通常够。

检测与求解

findChessboardCorners + cornerSubPix 亚像素细化；检测差的图（模糊、反光、遮挡）剔除——它们抬高 σ 还可能引偏。求解用 cv2.calibrateCamera，畸变模型别过度参数化：镜头若 5 系数够用就别上 8/14，多出来的弱约束参数会变成新的近零特征方向，反而让 I_β 病态。按镜头实际畸变量级选阶数。

验证诊断（闭环，别跳过）

把"看着不错"变成"确实可靠"，就是上一节的诊断套件，按代价从低到高：重投影 RMS 接近 σ_a（亚像素，好检测 < 0.3px，但必要不充分）→ 逐图误差看分布、剔除离群 → stdDeviations 每个内参都小到合理（fx 的不能是几万，最省事最该盯）→ 子集复现性 / leave-one-out 看稳定性 → cond(I_β) + 特征分解看退化与缺维 → 物理交叉校验（fx ≈ 焦距(mm)/像元(mm)、主点 ≈ 画面中心）挡系统偏差 → 留出图重投影测泛化。

精度预算与定稿

内参要连同 stdDeviations 一起报告，不能只给一个数。按下游任务反查：需要 fx 准到 0.1px 而 stdDev 是 0.15px → 这组不够，回去补数据；够用才定稿。Fisher 信息把"精度预算"变成拍前能估、拍后能验的量。

几条原理

多样性 > 数量：信息只在新的参数方向上累加。
条件数是命脉：目标是 I_β 良态，λ_min 远高于 ε·λ_max 地板。
噪声是地板：σ_a 决定单角点信息，清晰图 + 亚像素检测最划算。
CRLB 只管随机误差：系统偏差靠物理控制和交叉校验。
闭环验证：每改采集或模型，都用 stdDev + 特征分解 + 物理校验重新确认。

最后记住：真实流程里 Fisher 信息是在估计 θ̂ 处算的（观测信息），渐近等价于真值处的金标准。本文那条仿真流水线的价值，正是先用真值验证了"观测信息法可靠"（std/CRLB≈1），才敢把它用在真实数据上——仿真不是玩具，是真实流程的计量校准。那套真实诊断的完整 10 步实操（含 k 折、独立测试、数据一致性补回"没有真值就无法验证"这步）见《Fisher 信息：相机标定真实数据实操》。

参考资料

文章链接：
https://www.zywvvd.com/notes/study/camera-imaging/fisher-information-calib/fisher-information-calib/

“觉得不错的话，给点打赏吧 ୧(๑•̀⌄•́๑)૭”

微信支付

支付宝支付

Camera

#Camera