Lucas-Kanade光流法梯度计算疑问：简易方法为何效果更优？

问题背景 #

实现逐像素版Lucas-Kanade光流法时，用邻域4像素的简易方式计算梯度（dx、dy、dt）能得到符合预期的光流结果，但改用X/Y方向Sobel算子计算一阶导数后，光流结果极小且完全不符合图像实际运动方向，这一现象的核心原因可从两种梯度计算方式的差异入手分析。

两种梯度计算逻辑对比 #

简易邻域计算 #

这种方式本质是在2x2小邻域内，同时对前后帧的空间差分做时间维度平均，对时间差分做空间维度平均，自带平滑抗噪效果：

for (int y = 0; y < src.GetHeight() - 1; y++) {
    for (int x = 0; x < src.GetWidth() - 1; x++) {
        float src00 = src.GetPixelStart(x, y)[0];
        float src10 = src.GetPixelStart(x + 1, y)[0];
        float src01 = src.GetPixelStart(x, y + 1)[0];
        float src11 = src.GetPixelStart(x + 1, y + 1)[0];

        float dst00 = dst.GetPixelStart(x, y)[0];
        float dst10 = dst.GetPixelStart(x + 1, y)[0];
        float dst01 = dst.GetPixelStart(x, y + 1)[0];
        float dst11 = dst.GetPixelStart(x + 1, y + 1)[0];

        float srcVal = src10 - src00 + src11 - src01;                           
        float dstVal = dst10 - dst00 + dst11 - dst01;

        dFx.SetValue(0.25 * (srcVal + dstVal), 0, x, y);

        //------

        srcVal = src01 - src00 + src11 - src10;                
        dstVal = dst01 - dst00 + dst11 - dst10;

        dFy.SetValue(0.25 * (srcVal + dstVal), 0, x, y);

        //------

        dFt.SetValue(0.25 * (dst00 + dst10 + dst01 + dst11 - src00 - src10 - src01 - src11), 0, x, y);
    }
}

• 空间梯度（dx/dy）是前后帧各自邻域差分的平均，相当于对时间维度做平滑，降低单帧噪声的影响
• 时间梯度（dt）是2x2邻域内前后帧像素差的平均，通过空间平滑减少单像素噪声干扰

Sobel算子计算 #

当前实现直接对前后帧分别求Sobel梯度后相加，时间梯度采用单像素差分：

Image2d<float> dFx = (derivatorA.dx + derivatorB.dx);
Image2d<float> dFy = (derivatorA.dy + derivatorB.dy);
Image2d<float> dFt = (dst - src);

结果异常的核心原因 #

1. 噪声放大与矩阵奇异问题
   Sobel是3x3卷积核，会提取更大范围的空间梯度，但如果原始图像存在噪声，Sobel的差分计算会显著放大噪声，导致Lucas-Kanade求解时依赖的$\nabla\nabla^T$矩阵接近奇异。奇异矩阵的逆运算会产生数值不稳定的结果，最终得到极小的错误光流值。

2. 时间梯度未做空间平滑
   简易计算中dt经过2x2邻域平均平滑，而Sobel版本的dt是单像素差分，没有空间平滑，时间梯度的噪声极大。光流的基本约束方程$\nabla F \cdot [u, v]^T + F_t = 0$中，噪声大的$F_t$会直接导致解偏离真实运动方向。

3. 梯度计算的物理意义偏差
   Lucas-Kanade所需的是时空联合梯度，即$\nabla F = [\frac{\partial F}{\partial x}, \frac{\partial F}{\partial y}, \frac{\partial F}{\partial t}]$，其中$F(x,y,t)$是连续的时空图像。当前Sobel实现将前后帧的空间梯度直接相加（$\frac{\partial I_1}{\partial x} + \frac{\partial I_0}{\partial x}$），这和理论上的时空梯度（如$\frac{\partial (I_1 - I_0)}{\partial x}$）物理意义不符，导致梯度输入本身错误。

4. 缺少高斯平滑预处理
   Sobel算子通常需要配合高斯平滑使用（即高斯-差分组合），直接用Sobel求导会对高频噪声极度敏感。而简易计算的2x2邻域本身自带微弱的低通平滑效果，抗噪性远优于直接使用Sobel。

修复建议 #

• 先对前后帧分别做高斯平滑，再计算Sobel梯度，降低噪声干扰
• 修正梯度计算逻辑：计算差分图像$F_t = dst - src$后，对$F_t$求空间梯度，或采用时空联合的差分方式（如$\frac{\partial F}{\partial x} = \frac{I_1(x+1,y) - I_1(x-1,y) + I_0(x+1,y) - I_0(x-1,y)}{4}$）
• 对时间梯度$dt$也进行邻域平滑（如3x3窗口平均）
• 求解线性方程组时加入正则项（如添加小的对角矩阵），避免矩阵奇异导致的不稳定解

内容的提问来源于stack exchange，提问作者Martin Perry