3D 高斯泼溅实时渲染深度工程实战：从协方差 EWA 投影、Tile 可微光栅化到自适应密度控制与 HAC 压缩

3D 高斯泼溅实时渲染深度工程实战：从协方差 EWA 投影、Tile 可微光栅化到自适应密度控制与 HAC 压缩

执行摘要：体积渲染让 NeRF 拿到了照片级质量，却也把自己锁死在"每像素几百次 MLP 查询"的性能泥潭里。3D Gaussian Splatting（3DGS）换了一个范式：把场景显式表示为数百万个带各向异性协方差的 3D 高斯椭球，用 可微分光栅化 而不是 ray marching 求积分。这一改动把问题从"采样上百次积分"变成"一次 sort + 一次 forward blend"，1080p 下轻松跑到 100+ FPS。本文从工程角度拆开 3DGS 的主干：协方差参数化与 EWA 投影（为什么必须是 J W Σ Wᵀ Jᵀ + T）、Tile 分级排序（为什么 key 是 tile_id<<32 | depth）、自适应密度控制（split/clone/prune 的真实判定逻辑）、以及真正决定它能否上生产的 压缩（HAC / SOG）与伪影治理（popping / 拉伸 / 飞絮）。文中给出可直接跑的参数化代码与光栅化 kernel 骨架。

一、先算一笔账：为什么必须在 NeRF 之外找路

体积渲染的离散形式是这样的：每条射线在 [t_near, t_far] 上采样 N 个点，每点送一个 MLP 得到 (σ, c)，然后按数值积分累加。

$$C(r) = \sum_{i=1}^{N} T_i (1 - e^{-\sigma_i \delta_i}) c_i, \quad T_i = e^{-\sum_{j<i} \sigma_j \delta_j}$$

问题在于 N 与 MLP 都不能省。就算有了 InstantNGP 的 hash grid 把 MLP 压到几十层算乘法也很重要：1920×1080 分辨率 × 96 采样点 ≈ 2 亿次 MLP 前向。这就是为什么 NeRF 类方法在训练侧越来越快、在 渲染侧 始终进不了实时。

3DGS 的思路反转：不采样，画点。场景是一堆显式的高斯基元，渲染退化为一个排序后的 alpha blending —— 这是 GPU 从 1990 年代就开始优化的 workload。代价是你必须自己发明一套"如何控制这些高斯的分布"的训练机制，这就是后面第三节和第五节的全部内容。

二、数学骨架：高斯的参数化与 EWA 投影

一个 3D 高斯基元由四点决定：均值 μ（中心位置）、协方差矩阵 Σ（椭球的形状与朝向）、不透明度 α、以及视角相关的颜色 c（用球谐 SH 表达）。

G(x) = exp( -1/2 · (x - μ)ᵀ Σ⁻¹ (x - μ) )

Σ 不能直接优化。它是一个对称正定矩阵，直接把 9 个元素交给 SGD，几步就会变成非正定，导致 Σ⁻¹ 爆炸、椭球翻面。标准做法是把它拆成旋转与缩放两部分：

Σ = R S Sᵀ Rᵀ

其中 S = diag(s₁, s₂, s₃) 是各轴缩放向量，R 是旋转矩阵（用单位四元数 q 存储以避免万向锁与过参数化）。优化器只动 q、s、μ、α、SH 系数，Σ 是每帧重建出来的中间量，正定性由构造保证。

import torch

def build_covariance(scale: torch.Tensor, quat: torch.Tensor) -> torch.Tensor:
    """从 (scale[...,3], quat[...,4]) 构造世界空间协方差 Σ = R S Sᵀ Rᵀ。

    scale: 已经过 exp/log 激活，保证恒正
    quat : 单位四元数 (w, x, y, z)
    """
    q = quat / quat.norm(dim=-1, keepdim=True).clamp_min(1e-8)
    w, x, y, z = q.unbind(dim=-1)
    # 四元数 -> 旋转矩阵（行主序）
    R = torch.stack([
        1 - 2*(y*y + z*z),   2*(x*y - w*z),       2*(x*z + w*y),
        2*(x*y + w*z),       1 - 2*(x*x + z*z),   2*(y*z - w*x),
        2*(x*z - w*y),       2*(y*z + w*x),       1 - 2*(x*x + y*y),
    ], dim=-1).reshape(*q.shape[:-1], 3, 3)
    # Σ = R S Sᵀ Rᵀ = (R S)(R S)ᵀ，注意 S 是逐行缩放
    M = R * scale.unsqueeze(-2)          # R · S
    return M @ M.transpose(-1, -2)

EWA Splatting：从 3D 椭球到屏幕上的 2D 高斯

这是整个算法最精妙的一步，来自 Zwicker 2003 的 EWA surface splatting。我们要的是：把世界空间的 3D 高斯经过透视投影后，在屏幕空间仍然是一个 可以闭式求逆的 2D 高斯。

做法是 用雅可比做局部线性近似。设视图变换为 W，投影映射的局部线性化雅可比为 J，则：

Σ' = J W Σ Wᵀ Jᵀ  +  T

• J 是仿射近似的雅可比，编码了透视"近大远小"的局部形变；
• W 是 view matrix 的左上 3×3 旋转部分；
• T 是刻意加的低通滤波（dilation filter），大小为像素足迹量级。

T 这一项不能省。当某个高斯非常靠近相机或被拉得极扁时，投影后的 2D 协方差可能退化到 小于一个像素，此时像素采样会 miss 这个高斯，训练时出现忽闪忽灭的正负反馈 —— 表现为训练不收敛和刺点伪影。工程上还有一种流行做法是 screen-space 保底：把 Σ' 的最小特征值 clamp 到 0.3 px² 以上（抗锯齿版 3DGS 的核心就是把这个 clamp 做成有理论依据的 Mip 滤波，见第六节）。

def project_mean_cov(mu3d, cov3d, viewmat, focal, principal, img_size, filter_px=0.3):
    """标准 EWA 投影：返回屏幕坐标 μ' 与 2D 协方差 Σ'（ conic 形式）。"""
    t = viewmat @ torch.cat([mu3d, torch.ones_like(mu3d[..., :1])], dim=-1)
    tz = t[..., 2].clamp_min(0.3)                     # 近裁剪：防止除零与反向
    mu2d = torch.stack([
        (t[..., 0] * focal[0] / tz) + principal[0],
        (t[..., 1] * focal[1] / tz) + principal[1],
    ], dim=-1)
    W = viewmat[..., :3, :3]
    # 透视投影的雅可比（对相机空间坐标求导）
    J = torch.zeros(*tz.shape, 2, 3, device=mu3d.device)
    J[..., 0, 0] =  focal[0] / tz
    J[..., 1, 1] =  focal[1] / tz
    J[..., 0, 2] = -t[..., 0] * focal[0] / (tz * tz)
    J[..., 1, 2] = -t[..., 1] * focal[1] / (tz * tz)
    JW = J @ W
    cov2d = JW @ cov3d @ JW.transpose(-1, -2)
    # 低通滤波：沿对角线加像素级方差，等价于最小尺寸扩张
    cov2d[..., 0, 0] += filter_px
    cov2d[..., 1, 1] += filter_px
    return mu2d, cov2d

三、前向渲染：一次 alpha blending 就完事

到这一步渲染主循环极其简单。把 2D 协方差的逆记为 conic，则像素 x 处该高斯的功率为：

p(x) = exp( -1/2 · dᵀ Σ'⁻¹ d ),   d = x - μ'
α(x) = o · p(x)

按 深度从近到远 做 volume rendering 的离散化：

C = Σᵢ [ cᵢ · αᵢ · Π_{j<i} (1 - αⱼ) ]

其中 T = Π(1-αⱼ) 是累积透过率（transmittance），可以在一个寄存器里前向串行更新，并在 T < 1/255 时 提前终止 —— 这条 early-ray-termination 是性能的一半来源。

注意这里既没有立方体遍历，也没有二次 ray marching：全部 全部像素值来自一次 GPU 排序。这就是为什么 3DGS 能在一个 4050 显卡上跑 4K@60。

四、Tile 可微光栅化：一切都变成一次 key sort

难点从"渲染"转移到了"排序"。全逐像素 排序是 O(N log N)，而且无法并行。原始 3DGS 的做法是 分而治之 + 单次排序：

1. 屏幕划成 16×16（可配置）的 tile；
2. 每个高斯投影后算屏幕 AABB 向外扩一圈半径（3σ 覆盖 99% 概率质量），与它相交的每个 tile 生成 一个 instance；
3. key = (tile_id << 32) | packed_depth（32 位 tile id + 32 位浮点深度的可排序整数编码）；
4. 对所有 instance 做 一次 GPU radix sort（实践中用 OneSweep 单次遍历算法）；
5. 为每个 tile 记录它在排序后数组中的 [start, end) 区间（用一个 atomic + prefix sum 即可）；
6. 每个 tile 启动一个 thread block（典型 256 线程），把高斯元数据 load 进 shared memory，所有像素线程在这批共享数据上并行累加。

// tile 光栅化 kernel 骨架（每个 block 负责一个 16x16 tile）
__global__ void render_tile(const int2 dim,
                            const uint2* __restrict__ tile_ranges,
                            const uint32_t* __restrict__ sorted_keys,
                            const float2*  __restrict__ means2d,
                            const float3*  __restrict__ conic_and_opacity,
                            const float*   __restrict__ sh,
                            float*         __restrict__ out_rgb) {
    uint32_t tile_id = blockIdx.y * gridDim.x + blockIdx.x;
    uint2 range = tile_ranges[tile_id];
    uint32_t tid = threadIdx.x;
    // 协作把 [range.x, range.y) 区间的高斯载入 shared memory
    __shared__ float2 s_mean[256];
    __shared__ float4 s_conic_opacity[256];
    for (uint32_t base = range.x; base < range.y; base += blockDim.x) {
        uint32_t idx = base + tid;
        if (idx < range.y) {
            uint32_t g = sorted_keys[idx] & 0xFFFFFFFFu;   // 低 32 位是高斯索引
            s_mean[tid]          = means2d[g];
            s_conic_opacity[tid] = make_float4(conic_and_opacity[g].x,
                                               conic_and_opacity[g].y,
                                               conic_and_opacity[g].z, 0.f);
        }
        __syncthreads();

        uint32_t px = blockIdx.x * 16 + (tid % 16);
        uint32_t py = blockIdx.y * 16 + (tid / 16);
        if (px >= dim.x || py >= dim.y) continue;
        float T = 1.0f;                 // 累积透过率
        float3 C = make_float3(0.f, 0.f, 0.f);
        uint32_t n = min(range.y - base, blockDim.x);
        for (uint32_t j = 0; j < n; ++j) {
            if (T < 1.0f / 255.0f) break;            // early ray termination
            float2 d = make_float2(px - s_mean[j].x, py - s_mean[j].y);
            float power = -0.5f * (s_conic_opacity[j].x * d.x * d.x +
                                   s_conic_opacity[j].z * d.y * d.y)
                          - s_conic_opacity[j].y * d.x * d.y;
            if (power > 0.f) continue;               // 3σ 之外的像素直接跳过
            float alpha = min(0.99f, s_conic_opacity[j].w * expf(power));
            if (alpha < 1.0f / 255.0f) continue;
            C += /* SH 求值 */ evaluate_sh(sh, d) * alpha * T;
            T *= (1.0f - alpha);
        }
        if (T > 1e-4f) out_rgb[py * dim.x + px] = /* 背景 */ C;
        __syncthreads();
    }
}

工程关键点：

• tile size 的权衡：16×16 是经验最优。太小 → tile 数量暴涨，每个高斯的 instance 复制变多，radix sort 的 key 数量爆炸；太大 → 每个 tile 里高斯太多，逐像素线程要遍历的长尾变长，而且大部分像素对这个高斯的功率早已为 0，等于白算。
• 排序开销不可忽略：百万级高斯在高分辨率 下会产生 200 万 ~ 500 万 个 instance，sort 常常占到单帧时间的 30%~50%。这也是衍生工作（StopThePop、SortFree 高斯 rasterizer）的发力点。
• 反向传播必须重跑：中间量太多无法存下来，反向 pass 会按逆序重新累乘 T 并链式展开 —— 这也是显存里总是要屯一份"每个像素最后命中的 N 个高斯"的原因，训练峰值显存明显高于推理。

五、训练循环与自适应密度控制：算法的灵魂

如果只有前向 rasterizer，3DGS 就退化成一个"点云渲染器"。真正让它拿到 SOTA 的是 Adaptive Density Control —— 在训练过程中动态地分裂、克隆、剪枝高斯。

初始化来自 SfM（通常是 COLMAP）的稀疏点云，颜色用点的颜色，协方差初始化为"到最近 3 个邻居的平均距离的等比缩放"（保证初始椭球既不过大也不退化）。

训练目标是一个混合损失：

L = (1 - λ) · L1  +  λ · D-SSIM      (λ 通常 0.2)

纯 L1 会让结果偏糊，纯 SSIM 会让高频闪光；混合是标配。

然后每 100 个 iteration 执行一次密度控制。它的核心洞察是：位置梯度大的地方，说明这个区域还没被足够的高斯表达，需要加密。

@torch.no_grad()
def adaptive_density_control(gaussians, grad_accum, denom, step,
                             tau_pos=0.0002, scene_extent=1.0, max_screen_size=20):
    """每 100 步调用一次：先筛候选，再决定 clone 还是 split。

    返回需要删除的索引与新增的高斯，由调用方合并到模型中。
    """
    # 1) 累积的位置梯度（是对 N 张视角/mipmap 累积的平均）
    grad = grad_accum / denom.clamp_min(1.0)
    is_grad_high = grad.norm(dim=-1) > tau_pos

    # 2) scale 太大 -> 需要 split；太小（瓶颈尺寸不足）-> 需要 clone
    scales = gaussians.scales                      # 已激活后的实际尺度
    is_large = scales.max(dim=-1).values > (0.01 * scene_extent * max_screen_size / 100)

    to_split = is_grad_high & is_large             # 大 + 高梯度 -> 分裂
    to_clone = is_grad_high & ~is_large & (gaussians.opacity > 0.05)

    new_gaussians = []
    # split：以原高斯作为 PDF 采样两个新点，缩放除以 phi = 1.6
    for i in torch.where(to_split)[0]:
            samples = sample_from_gaussian(gaussians.mean[i], gaussians.cov[i], n=2)
        for s in samples:
            ng = gaussians[i].clone()
            ng.mean = s
            ng.scales = ng.scales / 1.6            # phi 超参，1.6 是论文值
            new_gaussians.append(ng)
    # clone：直接复制一个同残差兄弟
    for i in torch.where(to_clone)[0]:
        new_gaussians.append(gaussians[i].clone())

    # 3) prune：每 3000 步做一次 opacity reset + 清理低不透明度/过大的高斯
    low_opacity = gaussians.opacity < 0.005
    too_big = scales.max(dim=-1).values > (0.1 * scene_extent)
    to_prune = low_opacity | too_big

    return new_gaussians, torch.where(to_prune)[0]

为什么是 split 与 clone 两个分支？ 这是很有工程美学的设计：

• clone：这个高斯太小了、位置梯度大 → 说明这里细节还没填满，复制一个兄弟让两人各覆盖 一半区域，覆盖密度直接翻倍。
• split：这个高斯已经很大了却还在拉扯 → 说明它同时承担了多个几何/纹理的表达，把它裂成两个更小的子高斯，缩放除以 φ=1.6 保证总体积不暴涨但分辨率提升。
• prune + opacity reset：训练中总会有一些高斯跑向"透明"或"无限大"，它们不贡献画面却吃掉算力。定期 reset opacity 让它们重新竞争上岗，不行就删。

这套机制让高斯总数从初始的几千动态长到 100 万 ~ 500 万，而且长出来的分布是自适应的：复杂几何区域密，天空区域稀。

六、伪影治理：popping、拉伸、混叠是三个不同的问题

生产环境中用户最先投诉的从来不是 PSNR，而是 动起来的时候画面在抽。3DGS 有三类典型伪影，根因完全不同：

1. Popping（深度排序跳变）：tile 排序是 per-tile 的，相邻 tile 可能对同一组高斯给出不同顺序；当相机移动使某个 tile 的排序发生翻转时，整块的混合结果突变。StopThePop 的解法是引入 层级 tile 排序 + 交错（stochastic interleave）更新：让排序变化在时间上错开，把一次大跳变拆成多次微小变化，视觉上就察觉不到了。代价是每 N 帧要重做一次部分排序。
2. Dilation / Erosion（放大缩小时的体积胀缩）：根因是第二节里的低通滤波 T 是按 训练时 的像素足迹标的。训练用 1080p、部署时用户放大到 4K 或用不同 FOV，采样率变了，高斯的投影面积却没变 —— 于是要么出现飞絮（dilation，变胖），要么出现空洞（erosion，变瘦）。Mip-Splatting 的做法最简单有效：训练时对多频（不同分辨率）做混合训练，把 Σ' 的低通按实际采样率算，让模型学会"在同一个世界里同时做多个分辨率的正确渲染"。
3. Floaters（飞絮）：场景中飘浮的、不成表面的模糊团。根因是 view-dependent 的 SH 让高斯可以"用变色的方式解释任何视角"，从而不必真正贴合几何。解法是加几正则化：

• 各向异性 loss：惩罚椭球最薄轴，鼓励高斯变扁成 surfels（2DGS、GOF 就是把这个做彻底，直接把高斯基元换成 2D disk/surfel）；
• depth-normal 一致性：从像素导数推出的法线，与单目法线先验或 SfM 深度做一致性损失。

七、规模与压缩：从 500 MB 走到 10 MB

每个高斯的原始参数量：

300 万个高斯 ≈ 700 MB。这显然不能上线。三条主流路线：

• 量化 + 码本：先用 k-means 对 scale / rotation / SH DC 做 vector quantization 得到码本，其余高斯只存 index + residual，再过一遍 GPCC（几何点云压缩）做熵编码。典型压缩比 10×~20×。
• HAC（Hash-grid Assisted Context，CVPR'24 亮点）：压缩的核心其实是提升"下一个符号"的可预测性。HAC 用一个极紧凑的 hash grid 为每个高斯预测属性上下文（自回归），把剩余信息量压得很低，相比原始 3DGS 可达 ~75× 压缩同时画质几乎不降。这是目前质量/体积比最好的路线。
• SOG / Self-Organizing Gaussians：更工程化 ——先用 PLIC/网格排序把百万高斯重排成空间局部相关的 1D 序列，然后 把属性写成 2D 图像（每格对应一个高斯），用 WebP/PNG/JPEG-XL 编解码器直接压。好处是可以挂到标准图像 pipeline 上做流式传输与渐进加载，浏览器端解压成本极低。代价是不如 HAC 极致。

实践建议：如果目标是 Web / 移动端，优先 SOG 路线（顺带解决了"没有 compute shader 也能渲染"的问题）；如果是 离线归档 / CDN 分发 且对质量敏感，用 HAC。

八、生产环境真实会咬人的六个坑

1. 初始化失败是最大的单点故障。COLMAP 在无纹理墙面、镜面、重复纹理、弱光照下会直接崩。现代做法是换前馈式重建（DUSt3R / MASt3R / VGGT）先拿相机位姿与稠密点云，再进 3DGS。这个替换带来的稳定性提升远比调超参大。
2. 训练显存峰值远大于你算出来的参数量。因为反向 pass 需要保存每个像素的最后 N 个高斯。百万级高斯 + 1080p + batch 8 轻松吃到 40 GB+。对策：用小 batch + 累积、用 sparse Adam（只更新可见高斯）、必要时混合精度训练。
3. 旋转必须用四元数并每步归一化，不要用欧拉角，也不要不归一化 —— 漂移后椭球会翻转。
4. sort 的时间比 render 更容易成为瓶颈，尤其在 4K 或者含有大面积天空（几十个巨大高斯覆盖整个屏幕 → instance 数爆炸）的场景。给超大尺寸的高斯加一个"按面积公平拆分"的上限处理。
5. 不要指望单个 3DGS 直接渲染边界清晰的 3D 资产。它是 view-dependent 的场景表示，不是几何网格。若要导出 mesh 复用（例如进 Unity / Blender），需要额外做 surfels-to-mesh（GOF、SuGaR 这类），这一步的质量损失要预先在方案里告知。
6. 评估别只看 PSNR。PSNR 对糊出来的模糊很宽容。生产上要同时看 SSIM、LPIPS，以及最重要的 —— 一个"沿 unseen 轨迹自由飞行"的人工检查视频。飞一遍不该 pop 的地方有没有 pop、不该有 blob 的地方有没有 blob，比任何标量指标都准。

九、结论与判断标准

回到底层立项：3DGS 的本质贡献，是把 NeRF 的 "采样 → 逐点查网络 → 数值积分" 换成了 "投影 → 排序 → 前向混合"。前者是 MLP bound，后者是 memory bandwidth bound，而在现代 GPU 上，后者永远比前者快一到两个数量级。它所有的复杂度（split/clone/prune、EWA 投影里那个 +T、tile 排序里那个 key packing）都服务于一个目标：让这个过程保持闭式可微，同时不让排序吃掉全部收益。

要用不用它，问三个问题：你的场景是静态还是动态？（动态场景需要额外的形变场或逐帧重建，成本陡增）；你的部署终端有没有 compute shader？（没有 → 走 SOG + instanced quad 路线）；你能不能接受 view-dependent 表示而不是标准几何？（要导出资产就得加后处理）。三个都想清楚了，3DGS 大概率是当下质量-时延曲线上最靠外的那一个点。

最后一句实在话：3DGS 不是魔法，它是一个把"采样"换成"光栅化"的工程权衡。理解了 Σ' = J W Σ Wᵀ Jᵀ + T 里那个为什么非加不可的 T，以及为什么 key 是 tile_id << 32 | depth，你就不会被"某某 renderer 更快"的营销话术带偏 —— 你会直接去看它的排序策略和梯度流是怎么设计的。