MATLAB与遥感影像处理：做影像镶嵌或拼接时 MATLAB 报“transformation matrix is singular to working precision”，是重叠度不够、镜头畸变太大，还是应该从中间影像开始而不是从第一幅开始？

1. 报错的真实含义：从"矩阵奇异"到"系统病态"

MATLAB 的这条报错并非来自某个专门的"拼接函数"，而是来自线性代数内核。当代码执行 A\b、inv(A)、mldivide 或 fitgeotrans 内部的线性求解时，若矩阵被判定为数值奇异，就会抛出 Matrix is singular to working precision。关键在于"working precision"——它不是说矩阵在数学上严格不可逆，而是说在当前双精度（约 16 位有效数字）下，求逆结果已不可信。

在影像拼接语境中，被求逆的通常是单应矩阵估计的法方程矩阵 AᵀA（或 DLT 中的 AᵀA 的等价形式）。法方程的条件数是原始设计矩阵条件数的平方，这意味着原始矩阵条件数到 10⁸ 量级时，法方程已逼近双精度的可解边界。本文评述：多数工程报错发生在法方程层面而非原始矩阵层面，因此"换一个更稳的求解器"往往比"换一组点"收效更小——病根在数据构型，不在算法实现。

需要区分三类"奇异"：结构性奇异（控制点共线、重合，秩亏）、数值性奇异（构型尚可但坐标量级悬殊，条件数爆炸）、模型性奇异（真实变换不属于单应族，残差结构导致正规方程病态）。三者处置手段不同，混为一谈是排查效率低下的主因。

2. 单应估计的数学骨架与奇异判据

2.1 从对应点到齐次线性系统

设两幅影像上的对应点为 (x, y) 与 (x', y')，平面单应 H 为 3×3 齐次矩阵，满足

[x']   [h11 h12 h13] [x]
[y'] = [h21 h22 h23] [y]
[w ]   [h31 h32 h33] [1]

消去尺度因子 w 后，每对点贡献两个线性方程。n 对点堆叠得到 A h = 0，其中 A 为 2n×9。标准解法是对 A 做 SVD，取最小奇异值对应的右奇异向量；若使用归一化 DLT，则先对坐标做 Hartley 归一化再求解，最后反归一化。本文评述：MATLAB 的 fitgeotrans 在 'projective' 模式下内部已做一定归一化，但若用户自行实现最小二乘或调用 estimateGeometricTransform2D 之外的路径，归一化缺失是最常见的数值性奇异来源。

2.2 奇异判据：条件数与最小奇异值

对设计矩阵 A（2n×9）做 SVD 得奇异值 σ₁ ≥ … ≥ σ₉。工程上可用的判据有三条：其一，σ₉/σ₁ 低于约 10⁻⁸ 时，双精度下解不可信；其二，法方程矩阵 AᵀA 的条件数超过约 10¹⁶ 即触及双精度极限；其三，最小奇异值绝对值小于 eps * σ₁ * 9 量级时，MATLAB 的秩判定会将其视为零。本文评述：这三条判据在文献中常被混用，但适用对象不同——第一条针对原始矩阵、第二条针对法方程、第三条针对 MATLAB 内部秩判定。排查时应先明确报错发生在哪一层，否则容易在错误的层面反复调参。

3. 诱因一：重叠度不足如何退化为构型退化

3.1 重叠度不是标量，而是空间分布

"重叠度不够"在口语中常被理解为一个百分比，但在单应估计中真正起作用的是控制点在重叠区内的空间分布。即便重叠面积达到 40%，若匹配点全部挤在一条窄带或一个角落，设计矩阵仍会因构型退化而接近秩亏。遥感影像的条带式采集尤其容易出现这种"面积够、构型差"的情形：相邻航带的重叠区呈狭长矩形，匹配点沿航向高度相关，横向约束薄弱。

从几何上看，单应有 4 个自由度（8 参数，尺度归一后），要求控制点在平面上至少覆盖两个方向。若点集近似共线，则垂直于该方向的平移与旋转分量无法被约束，对应的奇异值趋近于零。本文评述：把重叠度问题转化为"点集协方差矩阵的最小特征值"问题，比单纯追求重叠百分比更贴近病根。工程上可计算匹配点坐标的 2×2 协方差矩阵，若最小特征值与最大特征值之比低于约 0.05，即应视为构型不良。

3.2 弱纹理与重复纹理的放大效应

遥感影像中的大面积水体、裸地、规则农田网格，会产生大量低区分度或周期性重复的匹配点。这类点并非随机分布，而是聚集在少数高对比结构上，进一步恶化构型。近年基于深度学习的局部特征（如 SuperPoint、LoFTR 系列）在弱纹理区有所改善，但在周期性纹理上仍存在误匹配聚集问题。本文评述：学习型匹配提升的是"点是否存在"，并不自动保证"点是否分散"，因此后端的几何筛选与空间均衡采样仍然不可省略。

4. 诱因二：镜头畸变与模型失配的耦合

4.1 畸变不是"让矩阵奇异"的直接原因

严格说，镜头畸变本身不会让单应矩阵奇异。畸变的作用是破坏单应模型的成立前提：单应描述的是平面到平面的射影变换，而带畸变的成像过程是非线性的。当畸变显著时，用单应拟合全部匹配点会产生系统性残差，RANSAC 为容纳这些残差不得不放宽阈值，进而纳入更多外点，最终使内点集的空间分布被污染、设计矩阵条件数恶化。本文评述：畸变是通过"污染内点集"这一间接通道触发奇异的，因此单纯增加匹配点数量无法缓解，反而可能加剧。

4.2 边界条件：何时畸变可以忽略

对于视场角较小、畸变已标定的航空或卫星推扫影像，残余畸变在亚像素量级，单应模型通常足够。对于大视场无人机影像、鱼眼或广角镜头，径向畸变在像幅边缘可达数十像素，此时应先在单幅影像上做畸变校正，再做拼接。工程判据可参考：若在像幅边缘的畸变残差超过匹配点定位精度的 2 至 3 倍，就应把畸变校正前置。本文评述：这一判据是经验性的，但方向明确——畸变校正的成本远低于在病态系统上反复调参的成本。

5. 诱因三：拼接顺序、误差累积与基准选择

5.1 "从第一幅开始"意味着什么

顺序拼接的本质是链式变换：第 k 幅影像的坐标通过 H₁→₂ · H₂→₃ · … · H₍ₖ₋₁₎→ₖ 累积映射到参考帧。以第一幅为基准时，误差沿链条单调累积，末端影像的几何偏差最大；同时，链式相乘会放大中间某一环的病态——若某一对影像的单应估计条件数偏高，其误差会向后传递并在后续配准中表现为"怎么调都对不齐"。本文评述：从中间影像开始并非玄学，其合理性在于把误差累积路径从"单向长链"改为"双向短链"，使最大累积长度减半，从而降低末端偏差与病态传播概率。

5.2 顺序策略的适用边界

需要指出，改变拼接顺序只能缓解误差累积，不能修复单对影像之间的结构性奇异。若某一对影像本身重叠不足或构型退化，无论从哪幅开始都会在该环节报错。因此顺序策略应被视为"二级优化"，而非"一级修复"。工程上更稳健的做法是引入全局平差：把所有影像的变换参数放入一个联合优化问题，用光束法平差或位姿图优化统一求解，从根本上避免链式累积。本文评述：对于数十幅以上的镶嵌任务，全局平差的收益远大于任何顺序技巧；顺序策略更适合影像数量少、算力受限的快速出图场景。

6. 四步处置路径：诊断、归一化、筛选、平差

6.1 第一步：定位报错层级

先确认报错来自哪一行。若是用户自写的 A\b，直接检查 A 的条件数；若是 fitgeotrans 或 estimateGeometricTransform2D 抛出，则应检查输入点对的数量、分布与是否含重复点。重复点会导致设计矩阵出现完全相同的行，秩亏是必然的。本文评述：诊断阶段最常见的浪费是把"点太少"与"点太挤"混为一谈，前者靠增加数量解决，后者必须靠空间均衡采样解决。

6.2 第二步：坐标归一化

Hartley 归一化的标准做法是：计算点集质心并平移到原点，再按平均距离缩放到 √2。对遥感影像而言，像素坐标动辄数千，若直接进入法方程，x² 项可达 10⁷ 量级，与常数项 1 相差七个数量级，条件数必然爆炸。归一化后各坐标量级接近 1，条件数通常可下降数个数量级。本文评述：归一化是性价比最高的一步，几乎零成本，应作为所有自写单应估计的默认前置步骤。

6.3 第三步：几何筛选与空间均衡

在 RANSAC 得到初步内点后，不要直接全部送入最小二乘。建议做两件事：其一，剔除坐标重复或距离过近的点对（例如像素距离小于 2 像素）；其二，对重叠区做网格划分，每个网格保留若干最强匹配，强制点集在空间上分散。本文评述：空间均衡会牺牲一部分匹配数量，但换来的是条件数下降与几何约束增强，在病态场景中通常是净收益。

6.4 第四步：分块与全局平差

当影像数量多、场景大时，单应模型本身可能不再适用。此时应转向分块拼接加全局平差的架构：先按空间邻接关系分块，块内用单应或仿射，块间用少量控制点做相似变换或多项式校正，最后用全局优化统一残差。本文评述：这一架构在开源摄影测量工具链中已有成熟实践，MATLAB 侧可借助优化工具箱实现简化版位姿图优化，核心是把"逐对估计"改为"联合估计"。

7. 代码实现：MATLAB 侧的可复现流程

7.1 条件数诊断片段

function [c, s] = diagnoseHomography(pts1, pts2)
% pts1, pts2: N×2 对应点
% 构造 DLT 设计矩阵（未归一化，用于暴露尺度问题）
N = size(pts1,1);
A = zeros(2*N, 9);
for i = 1:N
    x = pts1(i,1); y = pts1(i,2);
    u = pts2(i,1); v = pts2(i,2);
    A(2*i-1,:) = [-x -y -1  0  0  0  u*x u*y u];
    A(2*i  ,:) = [ 0  0  0 -x -y -1  v*x v*y v];
end
s = svd(A);
c = s(1) / max(s(end), eps);
fprintf('cond(A)=%.3e, sigma_min=%.3e\n', c, s(end));
end

该片段刻意不做归一化，目的是让尺度问题显性化。若 cond(A) 超过 10⁸，再对坐标做 Hartley 归一化后重算，对比条件数变化即可判断病根是否在尺度。

7.2 归一化 DLT 与均衡采样

function H = normalizedDLT(pts1, pts2)
% Hartley 归一化 + DLT
[T1, p1] = hartleyNormalize(pts1);
[T2, p2] = hartleyNormalize(pts2);
A = buildDLT(p1, p2);          % 同 7.1 的构造
[~,~,V] = svd(A, 'econ');
Hn = reshape(V(:,end), 3, 3)';
H  = T2 \ Hn * T1;
H  = H / H(3,3);
end

function [T, p] = hartleyNormalize(pts)
c = mean(pts, 1);
d = mean(sqrt(sum((pts - c).^2, 2)));
s = sqrt(2) / max(d, eps);
T = [s 0 -s*c(1); 0 s -s*c(2); 0 0 1];
p = (T * [pts, ones(size(pts,1),1)]')';
p = p(:,1:2);
end

均衡采样可在 RANSAC 内点集上执行：把重叠区划分为约 8×8 网格，每格保留匹配质量最高的 3 至 5 对点，再合并。本文评述：网格数不宜过细，否则每格样本不足反而引入噪声；经验上以每格至少保留 3 对点为下限。

7.3 顺序策略的工程实现

若坚持链式拼接，建议以"中间影像为根、向两侧扩展"的树形结构组织变换，而非线性链。实现上可先估计所有相邻对之间的单应，再以根节点为基准做一次广度优先的变换复合，最后可选地做一轮全局残差最小化。本文评述：树形结构把最大链长从 N 降到约 N/2，对误差累积的抑制是可量化的，但实现复杂度高于线性链，适合中等规模任务。

8. 关键数据与案例对照

下表汇总了不同构型与处理策略下单应估计条件数的典型量级。表中数值为模拟数据，由笔者在合成点集上按不同分布与尺度生成，用于说明趋势而非代表特定数据集。

表中趋势与经典数值分析结论一致：法方程条件数约为设计矩阵条件数的平方，归一化可显著降低尺度引起的病态，而构型退化（共线、团簇、重复）无法仅靠归一化修复。本文评述：这解释了为何"先归一化再看是否仍报错"是一个高效的二分诊断——若归一化后条件数降到 10⁴ 量级，病根在尺度；若仍在 10⁸ 以上，病根在构型。

9. 对比与局限：单应模型何时该被放弃

9.1 单应、仿射与多项式模型的适用边界

单应适用于平面场景或相机纯旋转情形；仿射是单应在无穷远平面下的近似，参数更少、更稳健，但对透视效应敏感；多项式校正（如二次、三次）能吸收部分地形起伏与镜头畸变，但外推能力差、控制点需求高。本文评述：模型选择应服务于场景几何，而非反过来用模型去硬套数据。当残差呈现明显的非射影结构（例如随像幅位置系统性弯曲）时，继续增加单应估计的匹配点只会让病态更隐蔽。

9.2 学习型方法与全局优化的进展

近三年，基于检测器-描述子联合学习与端到端匹配的方法在跨时相、跨模态遥感配准中进展明显，部分工作报告在弱纹理区的匹配召回率较传统方法有可观提升。但需注意，这些指标多在特定数据集上取得，迁移到新传感器、新季节时性能会下降。本文评述：学习型匹配解决的是"对应点从哪来"，并不直接解决"对应点如何组织成良态系统"，因此几何筛选、归一化与全局平差仍是不可或缺的后端。工程上更现实的路线是"学习型前端 + 稳健几何后端"，而非单点押注。

9.3 本文分析的局限

本文以平面单应为主线，未深入讨论地形起伏显著时的正射校正流程，也未对具体深度学习模型的网络结构与训练细节展开。文中条件数量级为合成数据趋势，实际项目应以自身数据的诊断结果为准。本文评述：这些局限不影响主线的因果判断，但在落地时应结合具体传感器模型与精度要求做补充验证。

10. 结论与展望

回到最初的问题：重叠度不够、镜头畸变太大、还是应从中间影像开始？本文的结论是——三者不在同一因果层级。重叠度通过控制点构型直接决定设计矩阵是否秩亏或近秩亏；畸变通过模型失配污染内点集，间接恶化构型；拼接顺序只影响误差累积路径，不能修复单对影像的结构性奇异。三者共同作用于法方程条件数这一通道，因此排查应按"先诊断层级、再归一化、再几何筛选、最后考虑全局平差"的顺序推进。

展望而言，随着全局平差与学习型匹配的成熟，单应估计的"病态焦虑"会部分让位于"系统建模"问题：与其在病态矩阵上反复调参，不如从一开始就采用联合优化框架，让几何约束在全局层面被正确表达。本文评述：工具会变，但"先看条件数、再谈算法"的诊断习惯不会过时。

主要参考文献

1. Hartley R, Zisserman A. Multiple View Geometry in Computer Vision. 2nd ed. Cambridge University Press, 2004.（经典单应估计与归一化 DLT 理论来源）
2. Fischler M A, Bolles R C. Random Sample Consensus. Communications of the ACM, 1981, 24(6): 381–395.（RANSAC 原始文献）
3. Lowe D G. Distinctive Image Features from Scale-Invariant Keypoints. IJCV, 2004, 60(2): 91–110.（SIFT 原始文献）
4. DeTone D, Malisiewicz T, Rabinovich A. SuperPoint: Self-Supervised Interest Point Detection and Description. CVPR Workshops, 2018.（学习型特征代表工作）
5. Sun J, Shen Z, Wang Y, et al. LoFTR: Detector-Free Local Feature Matching with Transformers. CVPR, 2021.（无检测器匹配代表工作）
6. Triggs B, McLauchlan P F, Hartley R I, et al. Bundle Adjustment — A Modern Synthesis. Vision Algorithms: Theory and Practice, 2000: 298–372.（全局平差经典综述）
7. Schönberger J L, Frahm J M. Structure-from-Motion Revisited. CVPR, 2016.（COLMAP 核心方法）
8. Zhang Z. Parameter Estimation Techniques: A Tutorial with Application to Conic Fitting. Image and Vision Computing, 1997, 15(1): 59–76.（数值病态与参数估计教程）
9. Golub G H, Van Loan C F. Matrix Computations. 4th ed. Johns Hopkins University Press, 2013.（条件数与奇异值分解数值基础）

说明：文中条件数量级表格为合成点集上的模拟数据，用于说明趋势；涉及数据集的处理均以原始文献描述为准。参考文献总数逾 60 篇，此处仅列主要 9 篇，其余在正文论述中按主题对应，引用时请以原始文献为准。