1. 从墙角到算法NLOS 非视域成像究竟在做什么这一讲是我计算成像系列的第 20 个主题也是最容易让人眼前一亮的一个方向NLOS 非视域成像。还没接触过这个词的话一句话说清楚——NLOS 让相机看见视线之外、被遮挡的目标靠的不是透视能力而是墙面上散射回来的光子。你打一束激光到一面普通的白墙上这面墙变成一面中继墙隐藏目标的散射光会再弹回墙面被单光子探测器记录。听起来离谱但最近十年这套东西已经从论文里的概念验证走到了走廊转弯、烟雾遮挡、生物组织浅层成像这些实际场景里。我最初入坑这个方向是被一段演示视频震住的激光逐点扫描墙面探测器记录的是墙上每个点的光子回波时间直方图算法跑完之后墙角后面一个字母牌的形状被清晰地重建出来。那种感觉怎么说呢有点像拿耳朵贴着墙听见了隔壁房间的谈话——不是直接看到但信息就是传过来了。这背后最核心的三板斧就是标题里这三个词瞬态数据transient负责怎么测光锥变换light cone transform负责怎么建模f–k 偏移f–k migration负责怎么快速重建。这套东西适合谁读如果你在做计算成像、光子级视觉、雷达/声呐信号处理或者单纯对隔墙有眼这类问题感兴趣这篇笔记能帮你把三者的关系串起来。我默认你有一点信号处理和傅里叶变换的基础但即便你刚入门我也会把每一步的物理直觉先讲清楚再给公式。1.1 一个看不见的场景怎么成像先建立一个最朴素的问题模型。场景里有一个隐藏物体它完全被遮挡相机看不到它。但场景里存在一面中继墙——墙面的一侧对着相机和激光器另一侧对着隐藏物体。激光脉冲打到墙面某个位置后一部分光子被墙面散射进入隐藏空间照亮隐藏物体隐藏物体再把光子散射回来又一次经过墙面最后被探测器接收。问题在于这些散射回来的光子经历了完全随机、混乱的多次散射为什么还能从中恢复出物体的形状答案藏在时间里。光子从激光器出发打到隐藏物体再回到探测器整个路径长度是已知的——因为光速是常数飞行时间直接对应路径长度。如果把隐藏场景划分成一个个体素那么某个墙点测到的某个时刻的回波强度就等于所有到该墙点距离恰好等于 c·t/2的体素散射贡献之和。测的数据是墙上位置 飞行时间要的答案是隐藏空间的三维体素强度这就是 NLOS 重建的核心矛盾。1.2 谁需要关心这套管线做自动驾驶的关心它因为车辆拐角处可能藏着行人做安防监控的关心它因为遮挡区域往往是案发区域做医学成像的也关心它因为有些组织表面之下藏着病灶但光子穿越深度有限NLOS 提供了一种绕路观察的可能。当然学术上它最大的价值是作为一个典型的反问题测量数据是被时间和空间耦合的正向过程简单反向重建却处处是坑。这篇文章后续的所有内容本质上都是围绕如何又快又稳地解这个反问题展开的。2. 瞬态测量NLOS 的第一手观测数据2.1 光子飞行时间的物理模型NLOS 的原始数据叫瞬态测量英文一般叫 transient measurement。它不像普通相机那样直接得到一张二维图像而是得到一堆时间分辨的直方图。具体来说用波长可调的脉冲激光器通常是可见光或近红外照射墙面上的某个点用单光子雪崩二极管SPAD或条纹相机记录回波。SPAD 的特点是灵敏度极高能记录单个光子的到达时刻配合时间相关单光子计数TCSPC技术重复触发成千上万次脉冲就能累积出一根光子到达时间直方图。横轴是飞行时间 t纵轴是光子计数。假设激光和探测器都对着墙上的同一个点 x这叫共焦配置confocal隐藏场景里有一个体素 p那么光子从 x 出发到 p 再回来飞行时间满足t 2·|p − x| / c其中 c 是光速。也就是说一个墙点 x一根直方图柱对应的不只是一个距离上的物体而是所有到 x 距离相等的物体的总和——这是三维空间里的一个球面。用公式写墙点 x 处测到的瞬态信号 I(x, t) 近似等于I(x, t) ∫ ρ(p) · δ(t − 2|p−x|/c) / |p−x|^4 dpρ(p) 是隐藏目标在 p 点的反射率分母上的 |p−x|^4 来自两程传播的逆平方衰减去程一次、回程一次各贡献 1/R²。注意这里假设墙面和目标表面都是理想朗伯散射体光子方向完全随机化。2.2 中继墙上的扫描策略共焦与非共焦实际采集时你要决定激光和探测器怎么扫描墙面。这里有两类经典配置我分别说一下。共焦配置激光和探测器始终指向墙上的同一点然后这个点在整个墙面上做二维栅格扫描。每一个墙点得到一根直方图最终数据是一个三维数组 (nx, ny, nt)nx、ny 是墙面扫描网格点数nt 是时间直方图的柱数。共焦的好处是几何关系简单——每个测量值对应一个球面反演模型简洁。坏处是慢因为要逐点扫描一个 64×64 的网格要测 4096 个点每个点还要累积大量光子统计。非共焦配置激光和探测器各自独立扫描或者激光固定、探测器扫描。这时候一个探测器墙点 xd、一个激光墙点 xl 和一个隐藏体素 p 满足的路径是折线xl → p → xd飞行时间 t (|p−xl| |p−xd|) / c对应的是椭圆面而不是球面。非共焦的优势是可以灵活部署探测器阵列比如几十个 SPAD 同时测采集速度大幅提升坏处是模型复杂度上升顺带的遮挡、基线距离等问题也更难处理。我在实际课题里两种都试过。共焦适合验证算法调试起来直观非共焦更接近真实部署场景但你的重建算法必须能处理不等权重的多基线段——这会给后面的 f–k 偏移带来额外麻烦后面细说。2.3 数据格式与预处理要点拿到原始数据之后第一件事不是上算法而是先把它弄干净。我总结了一个固定流程暗计数扣除SPAD 即使没有光照也有暗计数dark counts大概几十到几百 Hz 量级先测一段无光照的背景直方图从信号里减掉。环境光本底扣除室内环境光会均匀抬高所有时间柱的计数可以取信号段之外比如时间最早、激光脉冲还没回来的那几柱的平均值作本底。响应校准墙面不是理想朗伯体不同栅格点的反射率可能有差异需要先用一块已知高反射率的目标或直接在墙面前放一面镜子测一遍系统响应做归一化。时间零点校准直方图的时间原点不一定对准激光脉冲发射时刻系统触发延迟会造成整体的时间偏移这个偏移不校准重建出来的物体位置会在深度方向整体漂移而且会先糊后漂。还有一个容易被忽略的点直方图每个柱的宽度是探测器的时间分辨率决定的SPAD 常见在几十皮秒量级对应空间分辨率大约几毫米到一厘米。柱宽太宽会抹掉深度细节太窄则每个柱光子数不够、噪声爆炸。我一般先测一个点目标一块小反光板看系统的脉冲响应宽度再按它的 1/3 左右设置重采样时间间隔既保住分辨率又不至于砸进噪声里。3. 光锥变换把时间维度翻译成空间维度3.1 椭球面反投影的直觉拿到三维瞬态数据 I(x, t) 之后最直觉的重建思路就是反着来对隐藏空间里的每一个体素 p把墙面上所有点 x 在对应时间 t 2|p−x|/c 处的测量值拉出来累加到一起。这就是光锥变换也叫椭球面反投影。为什么叫光锥因为在 (x, y, t) 这个三维观测空间里一个固定的隐藏点 p 对应一个锥面随着墙点 x 变化飞行时间 t 和空间位置的关系画出来是圆锥面。反过来一个固定的测量时刻 t 对应一个空间曲面。光锥变换做的就是把这样的几何约束在时空域里显式地写出来然后把测量值沿锥面倒灌回体素空间。生活化的类比想象你在一个大操场上用一堆麦克风记录一次拍手的声音。每个麦克风听到啪的时间不同因为距离不同。如果你知道声速就可以把每个麦克风的到达时间转化成到拍手点的距离多个距离一交汇拍手的位置就出来了。NLOS 的每个墙点就相当于一个麦克风每个时刻的直方图柱就相当于一次到达时间记录只是这里信号不是声波而是散射光而且存在无数个拍手点隐藏物体的每个体素同时发声。数学上共焦 case 的反投影写出来是这个样子F(p) ∫ W(x, p) · I(x, 2|p−x|/c) dx其中 W(x, p) 是一个权重至少应该包含逆平方衰减的补偿项 |p−x|^4。如果不做任何滤波、直接累加结果会非常糊——因为每个测量值都被摊到了一个球面上球面相交区域是物体但球面本身会在体素空间里留下大片拖尾。所以实际使用的光锥变换不是纯反投影而是类似 CT 里滤波反投影的思路先对时间维做二阶导数或施加一个斜坡滤波器再进行反投影才能把模糊压下去。3.2 离散实现与体素坐标映射真正写代码的时候光锥变换的朴素实现是四重循环外层遍历体素 p内层遍历墙点 x。每个墙点 x 和体素 p 的距离算一次换算成时间索引 t_index直接从三维数组里取对应柱的值乘上权重累加到体素 p 上。伪代码大概长这样# volume: (Nz, Ny, Nx)transient: (Nt, Nyw, Nxw) for iz in range(Nz): for iy in range(Ny): for ix in range(Nx): acc 0.0 p get_voxel_coord(ix, iy, iz, origin, step) for jy in range(Nyw): for jx in range(Nxw): wp get_wall_coord(jx, jy, wall_origin, wall_step) d distance(p, wp) t_idx int(d / c / dt) if 0 t_idx Nt: acc transient[t_idx, jy, jx] * weight(d) volume[iz, iy, ix] acc这段代码能跑但性能很糟糕。假设体素空间 128×128×128墙面扫描 64×64时间柱 512总操作量是 128³ × 64² ≈ 137 亿次纯 Python 跑要几个小时。我实际做的时候做了三层优化第一把墙点循环和体素循环里最内层用 C 扩展或者 JIT 编译比如用 Numba能快一到两个数量级。第二体素坐标到时间索引的映射只跟空间几何有关、跟数据无关可以预先算好一张 lookup 表避免重复计算距离。第三利用对称性——某些配置下墙点关于体素是对称的可以只算一半再镜像累加。3.3 计算量与加速思路这里说一个重要的复杂度对比光锥变换朴素实现的复杂度是 O(Nx·Ny·Nz·Nxw·Nyw)即体素数和墙点数的乘积。64×64 网格、128³ 体素已经是一百多亿次操作。这种量级在离线验证时可以忍但要实时或近实时就完全不现实了。如果只扫一个方向的场景比如目标基本上是平面浮雕可以退化成二维问题复杂度降一个维度快很多。另一种加速思路是查表 GPU 并行把每个体素对应的墙点时间索引做成 (Nz, Ny, Nx, Nyw, Nxw) 的整数索引张量在 GPU 上用 gather 操作批量取数然后做加权累加能比 CPU 快上百倍。我们实验室最终跑 128³ 体素、64² 墙点、512 时间柱的场景用一张中端显卡重建时间从十几分钟压到了十几秒。不过光锥变换最大的问题不是速度而是精度反投影类方法对噪声敏感测量值里任何一个异常高亮的墙点比如反光强烈的表面都会沿着它的球面拖出一条亮线像手电筒扫过一样污染整个体素空间。这就是为什么要研究下一节的 f–k 偏移——它不是简单的反投影而是在频域里找一个全局一致的解。4. f–k 偏移傅里叶域里的快速重建4.1 从地震勘探借来的思路f–k 偏移这个名词听着陌生其实它的老家在地球物理勘探——石油行业用人工地震波探测地下结构地表布置一串检波器记录地震波到达时间和振幅反演地下反射层位置。这个地表记录波场、向下延拓成像的问题和 NLOS墙面记录散射波场、向隐藏空间成像的问题是同一个数学结构。我当初意识到这一点的时候第一反应是计算成像领域很多所谓的新方法其实是把地震学、雷达、超声这些老领域的问题换个物理场景重新讲了一遍。f–k 偏移就是最典型的例子上世纪七八年代地震勘探就用来做叠前偏移的 Stolt 算法到了 NLOS 这里只需要把波速换掉、把边界条件改一改就能直接用来做非视域重建而且做出来的速度比光锥变换快非常多。关键洞察是把测量数据重新解释为一个波场共焦配置下墙面上测到的瞬态信号 I(x, t)可以想象成在墙面这个平面上观察一个向上传播的标量波这个波以有效速度 c/2 传播因为光走了两倍距离去程加回程。隐藏物体里的每一个反射体素等价于在这个波场里的一个源点。于是问题变成了给出一平面上的波场记录反演出源的空间分布——这正是地震偏移的看家本领。4.2 频散关系与频域映射f–k 偏移的具体步骤分四步我按实际实现顺序写第一步对瞬态数据做三维傅里叶变换。沿墙面两个方向做空间 FFT沿时间轴做时间 FFT把数据从 (x, y, t) 域变换到 (kx, ky, ω) 域I~ (kx, ky, ω) FFT₃D[ I(x, y, t) ]第二步利用频散关系dispersion relation把时间频率 ω 映射到深度方向的波数 kz。均匀介质里标量波的频散关系是(kx² ky² kz²) (ω / c_eff)²其中 c_eff c/2 是 NLOS 共焦配置下的有效波速。因此深度的波数为kz sqrt( (ω / c_eff)² − kx² − ky² )注意当 (ω / c_eff)² kx² ky² 时kz 变成虚数对应倏逝波成分直接被舍弃。第三步在频域做 Stolt 插值。I~(kx, ky, ω) 采样点均匀分布在 (kx, ky, ω) 网格上但 ω 和 kz 不是线性关系所以要先把数据插值到均匀的 (kx, ky, kz) 网格上。这一步通常用线性插值或三次插值。第四步对 (kx, ky, kz) 域的三维数组做三维逆 FFT得到体素空间的重建结果F(x, y, z) IFFT₃D[ I~(kx, ky, kz) ]整个流程下来计算复杂度是 O(N³ log N)其中 N 是每个维度的采样数比光锥变换的 O(N⁴) 低一个量级。对于 128³ 输出体素f–k 偏移在 CPU 上也就几秒钟GPU 上亚秒级。这个速度优势直接决定了它适合做近实时系统。4.3 参数选择与伪影控制用 f–k 偏移有几个参数直接决定重建质量我把实际调试经验放在这里。有效波速 c_eff 是最关键的一个。理论上是 c/2但实际中因为墙面和目标散射是非理想朗伯体、探测器响应有时间展宽有效速度会偏离理论值。我习惯的做法是采集场景里放一个已知位置的小目标比如一张白纸折成的角重建后看它的深度位置是否准确不准确就微调 c_eff直到目标落在正确深度。这一步就像相机对焦c_eff 偏大目标深度会整体偏浅偏小则偏深。时间维的采样率也影响大。时间柱太宽ω 域的最高频率低重建的空间分辨率就受限时间柱太窄高频段基本是噪声反而会在 kz 域放大高频抖动。我通常把时间采样率设置为系统脉冲响应半高全宽的 1/2 到 1/3。傅里叶域的窗口函数不能省。直接对 I~(kx, ky, ω) 做插值和逆变换会在物体边缘产生明显的振铃伪影Gibbs 现象。我习惯在 FFT 之前先施加一个 Hanning 或者余弦窗对 kx、ky、ω 三个方向各做一次代价是分辨率略微下降但图像干净得多。还有一个常见的坑低频 DC 分量。墙面返回的直流背景光环境光、暗计数在 ω 域里集中在 ω 0 附近迁移之后会在重建体的边缘形成一圈亮晕。解决方法是在时间 FFT 之前做一次高通滤波直接扣除每个墙点直方图的平均值或者把 ω 接近 0 的那几层直接清零。我踩过这个坑第一次跑 f–k 偏移重建图四边各出现一条亮边排查了半天才发现是 DC 泄漏。5. 从数据到重建一条完整管线的坑与细节5.1 实验采集现场要注意的事纸上谈兵到这里应该聊点现场经验了。做实采的时候有四个容易忽略的环节。第一墙面选型。很多人以为随便找面白墙就行其实墙面的散射特性非常关键。理想的墙应该是高反射率、近似朗伯、双面不能有强镜面成分。实际中的墙面往往有一层呼吸产生的表面纹理或者涂料的散射不是完全各向同性这会直接导致测量值偏离理想模型。我后来养成的习惯是正式采集前先用一面小镜子放在墙面前对比直接反射和墙面散射两种信号的脉冲形状如果墙面散射的脉冲形状拖着长尾说明墙面存在多次散射后面预处理要加一个补偿项。第二扫描网格密度要匹配期望分辨率。墙面扫描步长太粗空间频率 kx、ky 的奈奎斯特上限就低重建分辨率天花板被压低。经验法则墙面步长应小于目标分辨率的 1/2。如果目标识别需要 1 cm 精度墙面试样步长至少取 5 mm再细也不好——因为时间分辨率本身限制了深度方向精度。第三激光功率和探测器计数率的平衡。SPAD 有死时间光子到达率太高会饱和直方图形状被压缩、甚至出现堆积效应导致时间偏移。我在调试时先把激光功率降到很弱观察单柱计数率在峰值处不超过探测器标称饱和率的 30%然后才开始正式扫描。第四时间同步。激光触发、SPAD 门控、扫描平台移动这三者之间要保持严格同步。扫描平台步进电机的振动、触发信号抖动都会在瞬态数据里表现为时间轴的随机扰动。我的做法是每扫描完一行插一个已知距离的参考目标测一次时间零点事后用插值把所有行对齐。5.2 重建参数标定与调试流程一套我反复使用的调试流程是这样先用模拟数据验证管线。用正向模型生成已知场景的瞬态数据加噪声分别跑光锥变换和 f–k 偏移。模拟阶段我发现的问题包括f–k 偏移对墙点视场覆盖不全特别敏感当墙面扫描区域小于隐藏场景的投影范围时频域会出现空洞重建结果有明显的条纹。解决方式是在 kx, ky 平面加一个渐晕窗让边缘频谱平滑衰减到零。然后做实采标定。场景里放三个已知位置的小反光球位置各不相同一个贴近墙面、一个中距、一个远距重建后看三个点的位置误差。这一步可以同时校准 c_eff、时间零点和体素坐标系与墙面的相对位姿。位姿这东西经常被忽略——如果墙面法向量和体素坐标系的 z 轴不完全平行重建结果会整体倾斜肉眼看着就是物体歪了。最后才做实际目标重建。每一步我都保存中间结果原始直方图、滤波后的瞬态、频域幅度谱、最终体素。因为反问题最讨厌的地方是结果错了不知道错在哪一步留中间结果是排查效率最高的习惯。5.3 质量评估与后续处理重建出来的是一个三维体素强度场不是干净的彩色照片。我评估质量的几个指标目标体素的对比度目标区峰值强度与背景均值之比、轮廓清晰度体素图里目标边缘的梯度、以及是否有对称鬼影镜像伪影。对比度至少要超过 3 dB 才算可靠低于这个值基本靠猜。后处理方面最常用的是三维中值滤波和总变差去噪。我的经验是先做一次 3×3×3 中值滤波去掉椒盐状噪声点再用总变差迭代去噪但要控制迭代次数迭代过头会把尖锐边缘也抹平。另外因为体素强度动态范围大显示时用对数变换或百分位截断比如把 99 百分位以上的强度映射到白色比直接线性显示清楚得多。6. 常见问题与排查速查最后把我在实践中遇到的高频问题整理成一张速查表每条都对应一个真实踩过的坑。现象大概率原因排查/解决方向重建图像整体模糊、边缘发虚时间采样过宽或墙点扫描步长过粗细化时间柱减小墙面扫描步长检查是否低于奈奎斯特采样率目标深度位置偏移c_eff 设置不准或时间零点未校准用已知位置目标标定 c_eff检查直方图时间原点与激光触发是否同步图像边缘出现明亮弧线/拖尾反投影类方法对高亮墙点敏感预处理做墙点归一化重建时施加权重窗抑制异常墙点贡献四边出现规则亮边框f–k 偏移里低频 DC 分量泄漏时间 FFT 前扣除每墙点均值或把 ω 低频层清零目标旁边出现对称鬼影非共焦配置下椭圆面交叉歧义或频域重建中 kz 正负分支未处理确认共焦/非共焦模型一致f–k 法中统一取正 kz 分支必要时用先验约束重建结果整体倾斜墙面法向量与体素坐标 z 轴不平行采集前用反光球标定位姿重建时在坐标变换里加入旋转修正强烈的高频闪烁噪声时间维过采样导致高频噪声放大在时间 FFT 前加 Hanning 窗或低通滤波匹配系统带宽目标被拉扁或压缩墙面扫描范围只覆盖了目标投影的一部分扩大扫描区域让隐藏目标的投影完整落在墙面上每一个问题我都可以展开讲半页但这里挑两个最容易反复踩的展开说。第一个是对称鬼影。f–k 偏移本质上是把瞬态数据当波场处理而波场反演天然不区分上行和下行的镜像解。在 NLOS 场景里墙面上的瞬态信号可能对应一个真实的隐藏物体也可能对应它在墙面后面的镜像位置。处理办法是在频散关系里指定 kz 的符号方向——只取向隐藏空间传播的那一支同时在重建后用一个墙面遮挡掩膜把墙面背后的镜像区直接抹掉。这个方法在实际场景里能把鬼影亮度压低一个量级但需要先知道墙面在体素坐标里的确切位置所以位姿标定是前提。第二个是环境光的干扰。白天房间里自然光很强SPAD 会收到大量非信号光子直方图整体抬高重建背景噪声显著增大。我在一个没有遮光窗帘的房间里做实验时经历过对比度从 8 dB 掉到 1.5 dB 的情况。解决方式有三种一是缩短直方图时间窗口只在激光脉冲往返的时间窗内开启 SPAD 门控能把环境光抑制掉一大半二是使用窄带滤光片只让激光波长附近的光进来三是事后扣除本底。三者叠加使用效果最好我最后的配置是门控时间窗 10 nm 窄带滤光片白天也能做采集只是信噪比比暗室还是低不少。从实践角度说NLOS 重建很少有一次跑通、图像完美的情况大多数时候是在标定、调参、去伪影之间反复迭代。这套瞬态数据采集 → 光锥变换建模 → f–k 偏移快速重建 → 后处理去伪影的管线已经是目前学术界和工业界的主流框架先行者用光锥变换把问题讲清楚后来者用 f–k 偏移把速度提上去。如果你理解了每个模块的物理含义和计算角色再去看近几年那些非共焦快速重建多基线段融合实时光子计数 NLOS的新论文会发现它们万变不离其宗都是在这三个环节上做文章。我个人在实际课题里最深的一点体会是NLOS 的方法论成熟度已经足以支撑工程化尝试但它永远受限于物理定律——墙面要够散射、探测器要够快、光子要够多、目标要够反射。遇到重建不出东西的时候别急着怀疑算法先把数据质量从头到尾检查一遍大部分问题都出在采集端的光学设计上而不是重建代码。这个方向后续再往下走我看好两个扩展一个是把深度学习和波动方程方法结合用网络做端到端的伪影抑制另一个是把 fast f–k 偏移做成实时管线真正推向自动驾驶和机器人感知场景。如果你也想动手试我建议从模拟数据开始先把光锥变换和 f–k 偏移在合成数据上都跑通再上实采——这样你踩坑时至少知道坑在数据端还是算法端。