# 解决Kinect中的多路径干扰问题：一种逆问题方法
Ayush Bhandari*、Micha Feigin*、Shahram Izadi†、Christoph Rhemann†、Mirko Schmidt†和Ramesh Raskar* 美国马萨诸塞州剑桥市麻省理工学院
电子邮件：{ayush.michaf.raskar}@mit.edu
$^{\dagger}$微软研究院，英国剑桥
电子邮件：{shahram,christoph,mirko}@microsoft.com
摘要—多路径干扰（MPI）是飞行时间（ToF）相机中深度和幅度测量误差的主要来源之一。这一问题近期受到了广泛关注。在本研究中，我们基于频谱估计理论和多频率测量框架探讨了MPI问题。与以往仅考虑最多两条干扰路径的方法不同，我们的模型涵盖了K条干扰路径的一般情况。在理论层面，我们证明对于存在K条光干扰路径的情况，$2K + 1$次频率测量就足以恢复出每条光路径对应的深度和幅度值。我们方法的独特之处在于其算法实现是非迭代的，因此能得到计算效率较高的闭式解。此外，我们首次在市售的微软Kinect for Xbox One上验证了该模型的有效性。
## 一、引言
振幅调制连续波（AMCW）飞行时间（ToF）成像相机[^1]、[^2]能够在每个像素处同时测量幅度和光传播时间（深度），从而获取三维场景信息。图3(e)和3(f)展示了ToF相机生成的此类强度和深度图像示例。这类相机的工作原理是：通过对光源进行振幅调制来发射编码光信号（通常为正弦波），然后测量信号从发射到从场景反射回来的时间延迟（图1(a)），其原理与激光雷达类似。
### A. ToF相机中的多路径干扰
所有现有的ToF相机都基于这样一个假设：每个像素只观测到一条光路径。换而言之，该假设认为场景仅被直接照明，不存在内部反射（在光学和地球物理学领域被称为玻恩近似[^3]）。然而，在许多实际场景中并非如此，比如由于多个物体或拐角导致的场景内部反射、存在窗户等透明物体时，或是发生次表面散射的情况。图1展示了其中的一些实例。
当多条光路径在某一个像素处叠加时，深度测量结果就会被破坏。这就是所谓的多路径干扰问题或混合像素问题（MPI）。它是ToF传感器误差的主要来源之一[^4] [^5] [^6] [^7] [^8] [^9] [^10] [^11] [^12] [^13] [^14] [^15]。为此，几乎所有现有解决方案都只考虑最多两条光干扰路径[^8] [^9] [^10]。K条干扰路径的情况最早在[^4]中被提及，随后在[^5]中也有论述。这两篇论文都采用了稀疏性建模方法，
(a)
(b)
(c)
(d)
Fig. 1: (a) The ToF camera emits a reference signal. Delay in time of arrival of the reflection encodes the depth d. (b) Specular or mirror like and, (c) Semi-transparent reflections cause multiple light paths to mix at the sensor which leads to multipath interference or the mixed-pixel problem (MPP). (d) Case of continuous multipath reflections.
从而导致算法计算量大且需要迭代求解。
### B. 本文的贡献与结构
在本研究中，我们针对K路径干涉问题提出了一种非迭代求解模型，该模型在实现上无需迭代计算，因此与[^4]、[^5]中讨论的基于稀疏性的迭代方案相比，具有更高的计算效率。
在第II节中，我们对ToF图像形成模型进行了总体描述。在我们的逆问题框架下，我们建立了多路径干涉问题的数学模型。正如我们将看到的，对该问题的重新表述表明，多路径干涉本质上与参数谱估计理论密切相关[^16]、[^17]。我们利用已有的研究成果[^6]、[^16]来以闭式形式求解该问题。在第III节中，我们讨论了与微软Kinect One中的多路径干涉消除相关的一些初步结果。最后，我们总结了未来的研究方向。
---
Fig. 2: Time-of-flight imaging system pipeline for solving inverse problems. p is the transmitted signal, h the system (scene) response, r the measured signal and m the correlation measurement.
### II. ToF图像形成模型
ToF成像流程如图2所示。ToF相机的像素级图像形成模型可分为3个步骤。
1) 探测场景：ToF相机属于主动式传感器。ToF相机通过发射形式为$p(t) = p(t + T)$的T周期函数来探测场景。这一功能是通过对通常位于传感器处的光源进行幅度调制来实现的。在当前相机中，该信号为$p(t) = 1 + \cos(\omega_0 t)$、$\omega_0 = 2\pi/T$。
$$ r\left(t\right)=\int_{\Omega}p\left(z\right)h\left(t-z\right)d z\quad\mathrm{(R e f l e c t e d~S i g n a l).} $$
2) 场景交互：探测信号p会与场景相互作用，从而产生反射信号s。这种相互作用可通过场景响应函数h在数学上进行定义。
在标准的多路径情况下，场景交互函数是对称的。而对于更复杂的逆问题，比如荧光寿命成像显微镜（FLIM），该函数可能是非对称的，即$h(z,t)$，此时响应函数可由更通用的弗雷德霍姆积分来描述[^18]。
$$ m\left(t\right)\Big|_{t=k\Delta}=\left(p\otimes r\right)\left(t\right),~k\in\mathbb{Z} $$
3) 互相关运算：ToF锁相传感器会针对某个采样步长$\Delta$，测量参考信号与反射信号的离散化互相关函数。
此处$\otimes$表示相关函数。
### A. 单深度建模
包括Xbox One版Kinect在内的大多数ToF相机，均采用幅度调制的正弦波来探测场景。
$$ p\left(t\right)=1+\cos\left(\omega t\right). $$
在某个深度d处，探测信号会发生反射，同时还会根据当前像素处的场景反照率发生衰减$0 \leq \Gamma \leq 1$。总的飞行时间为$t_d = 2d/c$，其中$c = 3 \times 10^8 m/s$为光速。
简单的计算[^1]、[^2]表明，锁相像素处测得的相关响应为：
$$ m_{\omega}\left(t\right)=\Gamma\left(1+\cos\left(\omega\left(t+t_{d}\right)\right)\right). $$
通过将场景响应函数设为以下形式，也可得出相同的结论：
$$ h\left(t\right)=\Gamma\delta\left(t-t_{d}\right), $$
其中$\delta$表示狄拉克δ函数。
接下来需要从公式(2)中的测量结果m来估算$\{\Gamma, d\}$。大多数ToF相机属于所谓的4抽头设备，它们利用4个相关测量值来在采样步长$\{\widetilde{\Gamma}, \widetilde{d}\}$处计算近似值$\Delta = \pi/2\omega, k = 0\cdots3$ [^2]。
$$ \begin{array}{r l}&{\widetilde{\Gamma}=p_{0}^{-2}\bigvee\left(m_{\omega}\left(3\Delta\right)-m_{\omega}\left(\Delta\right)\right)^{2}+\left(m_{\omega}\left(0\right)-m_{\omega}\left(2\Delta\right)\right)^{2}}\\ &{\widetilde{d}=\frac{c}{2\omega}\mathrm{t a n}^{-1}\left(\frac{m_{\omega}\left(3\Delta\right)-m_{\omega}\left(\Delta\right)}{m_{\omega}\left(0\right)-m_{\omega}\left(2\Delta\right)}\right).}\end{array}\quad(6) $$
需要注意的是，对于任意给定的调制频率$\omega$以及估算出的参数$\{\Gamma, d\}$，我们都可以将其对应为一个复数/相量。
$$ Z_{\omega}=\widetilde{\Gamma}\exp\left(\jmath\left(2\widetilde{d}\omega\right)/c\right)\triangleq\widetilde{\Gamma}\exp\left(\jmath\omega t_{d}\right),\quad Z_{\omega}\in\mathbb{C}. $$
需要注意的一点是，$Z_{\omega}$ 是 $\Gamma\delta(t-t_{d})$ 的傅里叶变换。
现在我们将注意力重新转向多径干扰的情况。
### B. 多径干扰问题
考虑图1(b)和1(c)中讨论的情况。当多次反射的光或光路到达传感器时，深度测量结果会被破坏[^4]、[^8]、[^14]。基于公式(5)建模的系统响应函数为
$$ h_{K}\left(t\right)=\sum_{k=0}^{K-1}\Gamma_{k}\delta\left(t-\frac{2d_{k}}{c}\right)\equiv\sum_{k=0}^{K-1}\Gamma_{k}\delta\left(t-t_{k}\right). $$
我们用$h(t)$来表示$\widehat{h}_{K}(\omega)$的傅里叶变换，即
$$ \begin{array}{r}{h_{K}\left(t\right)\xrightarrow{\mathrm{F o u r i e r}}\widehat{h}_{K}\left(\omega\right)=\sum_{k=0}^{K-1}\Gamma_{k}e^{-\jmath\omega t_{k}}.}\end{array} $$
实际上，对于多径情况，公式(8)中系统响应函数的傅里叶域表示是K个复指数函数的和。利用相量公式(7)，测量结果可以用这个复数值$Z_{\omega} = \widehat{h}_K^*(\omega)$来表示。由此我们得到关于$\{\widetilde{\Gamma}, \widetilde{d}\}$的以下公式
$$ \widetilde{\Gamma}_{\omega}=\left|\sum_{k=0}^{K-1}\Gamma_{k}e^{j\omega t_{k}}\right|\mathrm{~a n d~}\widetilde{\phi}_{\omega}=\angle\left(\sum_{k=0}^{K-1}\Gamma_{k}e^{j\omega t_{k}}\right), $$
其结果是相机无法解析受干扰破坏的测量结果。
### C. 恢复$\{\Gamma_k, d_k\}_{k=0}^{K-1}$的闭式解
我们从测量值$m_{\omega}(t)$开始。由于多径分量与t无关，我们仅考虑在$\omega$中的测量值。
问题描述：给定一个包含N个等间距调制频率下的测量值的向量m，其中
$$ m_{n}=\sum_{k=0}^{K-1}\Gamma_{k}u_{k}^{n},n\leqslant N with u_{k}=\exp\left(jt_{k}\right) $$
提取出$\{\Gamma_k, d_k\}_{k=0}^{K-1}$
这是线谱估计中的经典逆问题[^17]，也是分数阶傅里叶变换背景下所讨论解的一个特例[^16]。
---
我们定义一个洛朗多项式，使其根为$u_{k}$，
$$ \textsf{Q}\left(z\right)=\prod_{\ell=0}^{K-1}\left(1-u_{k}z^{-1}\right)\equiv\underbrace{\sum_{\ell=0}^{K}q_{\ell}z^{-\ell}}_{\mathrm{F a c t o r i z e d~P o l y n o m i a l}},\quad z\in\mathbb{Z}. $$
$$ q_{\ell}\ast m_{\ell}=m_{\ell}\mathsf{Q}\left(u_{k}\right)=0. $$
$$ q_{\ell}{*}m_{\ell}=\sum_{n=0}^{K}q_{n}m_{\ell-n}=\underbrace{\sum_{k=0}^{K-1}\Gamma_{k}u_{k}^{\ell}}_{m_{\ell}}\underbrace{\sum_{n=0}^{K}q_{n}u_{k}^{-n}}_{\mathsf{Q}(u_{k})=0}=0. $$
这是因为$\left.Q\left(z\right)\right|_{z=u_k}=0$（见公式(12)）。因此该多项式的零空间包含了不同路径的解耦相位信息。用向量-矩阵表示法来说，这等价于
$$ \underbrace{\left[\begin{array}{c c c}{m_{K}}&{\cdots}&{m_{0}}\\ {\vdots}&{\ddots}&{\vdots}\\ {m_{2K}}&{\cdots}&{m_{K}}\end{array}\right]}_{}\underbrace{\left[\begin{array}{c}{q_{0}}\\ {\vdots}\\ {q_{K}}\end{array}\right]}_{\mathbf{q}^{K+1}}=\underbrace{\left[\begin{array}{c}{0}\\ {\vdots}\\ {0}\end{array}\right]}_{\mathbf{0}^{K+1}}\Leftrightarrow\mathbf{M}\mathbf{q}=\mathbf{0}, $$
其中$\mathbf{M}$被识别为卷积/托普利茨矩阵，向量$\mathbf{q} \in \mathrm{Null}(\mathbf{M})$是秩亏矩阵$\mathbf{M}$的零空间中的一个向量，其秩为K。我们至少需要$N = 2K + 1$个测量值来构造$\mathbf{M}$，因此需要$2K + 1$个调制频率下的$m_\omega$测量值才能恢复K条光学干涉路径。一旦计算出$\mathbf{q}$，就对其进行因式分解以得到公式(12)中的根$\{u_k\}$。进而可得
$$ \widetilde{d}_{k}=\left(2\jmath\right)^{-1}c\operatorname{l o g}\left(u_{k}\right)\qquad\mathrm{w h e r e}\qquad c=3\times10^{8}m/s. $$
我们定义一个范德蒙矩阵$\mathbf{U}^{N\times K}$，其矩阵元素为
$$ [U]_{n,k}=u_{k}^{n},\quad k=0,\ldots,K-1,\quad n=0,\ldots,N-1. $$
为了计算$\{\Gamma\}_{k}$，我们求解方程组$\mathrm{UT} = \mathbf{m} \Rightarrow \Gamma = \mathbf{U}^{+}\mathbf{m}$，其中$\mathbf{U}^{+}$是$\mathbf{U}$的伪逆。
### D. 备注与讨论
- 相较于以往的工作有所改进。首先，我们的工作推广了[^9]、[^10]中的结果，因为我们的理论结果适用于K条干涉路径的一般情形。例如，当K=2时，需要N=5个频率测量值，这正是[^9]中讨论的结果。此外，与需要同时测量频率和相位的[^9]不同，我们的公式无需进行相位采样。
在 TOF 背景下，利用多个调制频率来分辨多条相互干扰路径 $(K > 2)$ 的思想最早在 [^4] 中提出；然而，在稀疏优化中选用贪心算法会使求解效率低下。在本文中，我们证明 $\{\Gamma_k, d_k\}_{k=0}^{K-1}$ 可以通过闭式解从 $m_\omega$ 中恢复。
- 重建保证。在不存在模型不匹配和噪声的情况下，我们证明需要$N = 2K + 1$次频率测量才能对K条干扰路径进行解混。据我们所知，这一结果在文献中尚未被报道（在TOF成像领域）。此外，针对我们的实验设置，可以利用Cramér–Rao界限[^17]来提供重建保证。
- 实现细节。众所周知，Vandermonde矩阵$U^{+}$在存在噪声时稳定性极差，而过采样是一种有效的应对手段。文献中已提出了多种缓解该问题的方案。为实现实际应用，我们采用Matrix Pencils[^19]来求解该问题。
### 三、实验结果
实验中我们使用了一台市售的微软Xbox One版Kinect TOF相机，其固件经过修改以支持自定义调制频率。Kinect距离场景约1.5米，该场景为一张桌子上的人体模型，背景是一面墙。我们还在相机前方约15厘米处放置了一块透明亚克力板。亚克力板会将部分发射光直接反射回相机，从而产生混合测量数据（即多路径效应）。我们从50 MHz开始，记录Kinect在不同调制频率下的幅度/相位测量值。
实验装置如图3(j)所示。在此示例中，我们设定$K = 2$。52 MHz下的幅度和相位测量结果分别如图3(a)和3(b)所示。我们选取52–72 MHz范围内的21组测量数据进行后续处理。对于图3(j)中的实验，我们的目标是将混合图像分解为两个分量，分别对应两条光路：一条来自亚克力板，另一条来自背景。图3(c)和3(d)展示了由亚克力板产生的幅度/深度配对结果。从图3(c)中可以看出，由于反射性亚克力板产生了镜面反射现象。图3(e)和3(f)则展示了背景部分的贡献；桌面的深度经估算为1.67米。图3(e)中的幅度图像没有受到亚克力板引起的多路径效应的影响，包括镜面反射。
图3(g)和3(g)分别展示了52 MHz和72 MHz下的三维重建结果（幅度信息嵌入在深度数据中）。幅度恢复误差是多路径干扰的明显特征。图3(i)展示了两次反射情况下幅度信息嵌入三维数据的效果。
### 四、结论
我们的理论推导具有通用性，能够涵盖K条光路的一般情况；同时该算法无需迭代，计算效率较高（与近期的相关研究相比）。初步结果证明了该模型的实际可行性以及计算上的便捷性。
## 参考文献
[^1]: M. Hansard, S. Lee, O. Choi, and R. Horaud, Time-of-Flight Cameras: Principles, Methods and Applications. Springer, 2013.
## [^2]: R. Lange, “3d time-of-flight distance measurement with custom solid-state image sensors in cmos/ccd-technology,” Ph.D. dissertation, Department of Electrical Engineering and Computer Science at University of Siegen, 2000
---
Fig. 3: Images (a) and (b) show one of the input amplitude and phase images while (g) and (h) show the embedding of the intensity image onto the depth data. The scene, depicted in (j) is constructed from a mannequin on a table with a transparency between the camera and the scene. As can be seen, there is a strong flashback from the transparency as well as corruption of the depth image due to multi bounce. (c) and (d) show the intensity and depth reconstruction of the transparency while (e) and (f) show that of the scene. (i) Shows the embedding of the intensity image onto the depth reconstruction. Not only have the intensity and depth information been recovered correctly for both the scene and transparency, we have also recovered both amplitude and depth information that was completely occluded by the flashback (highlighted).
[^3]: M. Born and E. Wolf, Principles of Optics: Electromagnetic Theory of Propagation, Interference and Diffraction of Light, 7th ed. Cambridge University Press, 1999.
[^4]: A. Bhandari, A. Kadambi, R. Whyte, C. Barsi, M. Feigin, A. Dorrington, and R. Raskar, “Resolving multipath interference in time-of-flight imaging via modulation frequency diversity and sparse regularization,” Optics Letters, vol. 39, no. 7, 2014.
[^5]: D. Freedman, E. Krupka, Y. Smolin, I. Leichter, and M. Schmidt, “SRA: Fast removal of general multipath for ToF sensors,” arXiv:1403.5919.
[^6]: A. Bhandari, A. Kadambi, and R. Raskar, “Sparse linear operator identification without sparse regularization? Applications to mixed pixel problem in time-of-flight / range imaging,” in Proc. IEEE Int. Conf. Acoust., Speech, and Signal Proc. (ICASSP), 2014, pp. 365–369.
[^7]: A. Kadambi and R. Whyte and A. Bhandari and L. Streeter and C. Barsi and A. Dorrington and R. Raskar, "Coded time of flight cameras: sparse deconvolution to address multipath interference and recover time profiles," ACM Trans. on Graphics (TOG), vol. 32, no. 6, p. 167, 2013.
[^8]: J. P. Godbaz, A. A. Dorrington, and M. J. Cree, “Understanding and ameliorating mixed pixels and multipath interference in amcw lidar,” in TOF Range-Imaging Cameras. Springer, 2013, pp. 91–116.
[^9]: A. Kirmani, A. Benedetti, and P. A. Chou, “Spumic: Simultaneous phase unwrapping and multipath interference cancellation in time-of-flight cameras using spectral methods,” in IEEE Intl. Conf. on Multimedia and Expo (ICME), 2013, pp. 1–6.
[^10]: J. P. Godbaz, M. J. Cree, and A. Dorrington, “Closed-form inverses for the mixed pixel/multipath interference problem in amcw lidar,” in IS&T/SPIE Electronic Imaging, 2012.
[^11]: D. Jiménez, D. Pizarro, M. Mazo, and S. Palazuelos, "Modelling and
“飞行时间相机中多路径干扰的校正”，载于《计算机视觉与模式识别会议（CVPR）》，2012年，第893–900页。
[^12]: A. Jongenelen, D. Bailey, A. Payne, A. Dorrington, and D. Carnegie, "Analysis of errors in ToF range imaging with dual-frequency modulation," IEEE Trans. on Inst. and Meas., vol. 60, no. 5, pp. 1861–1868, 2011.
[^13]: A. A. Dorrington, J. P. Godbaz, M. J. Cree, A. D. Payne, and L. V. Streeter, “Separating true range measurements from multi-path and scattering interference in commercial range cameras,” vol. 7864. SPIE, 2011, pp. 786404–786404–10.
[^14]: S. Fuchs, “Multipath interference compensation in time-of-flight camera images,” in Intl. Conf. Patt. Rec. (ICPR), 2010, pp. 3583–3586.
[^15]: F. Heide, M. B. Hullin, J. Gregson, and W. Heidrich, “Low-budget transient imaging using photonic mixer devices,” in ACM Transactions on Graphics, 2013.
[^16]: A. Bhandari and P. Marziliano, “Sampling and reconstruction of sparse signals in fractional Fourier domain,” IEEE Signal Processing Letters, vol. 17, no. 3, pp. 221–224, 2010.
[^17]: P. Stoica and R. Moses, Introduction to spectral analysis. Prentice Hall Upper Saddle River, NJ, 1997, vol. 89.
[^18]: A. Bhandari, C. Barsi, R. Whyte, A. Kadambi, A. J. Das, A. Dorrington, and R. Raskar, "Coded time-of-flight imaging for calibration free fluorescence lifetime estimation," in Imaging Systems and Applications. Optical Society of America, 2014, pp. IW2C-5.
[^19]: H. Yingbo and T. K. Sarkar, “Matrix pencil method for estimating parameters of exponentially damped/undamped sinusoids in noise.” IEEE Trans. Acoust., Speech, and Signal Proc., pp. 814–824, 1990.