基于优化Tikhonov最小二乘模型的曲面重建

Surface Reconstruction Based on Optimized Tikhonov Least Squares Model

摘要: 曲面3D重建是通过SfS系列方法从灰度图像中获取梯度场后，采用L2范数建立深度函数的能量泛函进行重建，这是一个不适定问题。本研究一方面通过将三角测量获取的粗糙深度图Z0与深度函数进行拟合作为Tikhonov正则项，将其转化为适定问题，并用L曲线选取合适的正则化参数；另一方面，利用先离散再优化的方法将上述问题表示为Sylvester矩阵方程，使用Hessenberg-Schur方法求解得到待重建曲面。数值实验表明，优化后的模型抑制噪声和偏差的效果更优，重建精度更高。

Abstract: Surface 3D reconstruction is to obtain gradient field from gray image by SfS series method, and then reconstruct it by using L2 norm to establish the energy functional of depth function, which is an ill-posed problem. In this study, on the one hand, the rough depth map Z0 obtained by triangulation is fitted with the depth function as Tikhonov regular term, which is transformed into an adaptive problem, and the appropriate regularization parameters are selected with L curve. On the other hand, the above problems are expressed as Sylvester matrix equations by the first discretization and then optimization method, and the reconstructed surfaces are obtained by Hessenberg-Schur method. Numerical experiments show that the optimized model can suppress the noise and devia-tion better, and the reconstruction accuracy is higher.

1. 引言

三维重建是利用测得曲面的二维梯度场或表面法向的矢量数据来恢复物体三维信息的过程，它对许多应用都是必不可少的，如地震成像 [1] 、Hartmann波前检测 [2] 、光栅投影 [3] 和三维精密测量 [4] 等。假

设图像通过正交投影获得，其3D表面为 $z=z\left(x,y\right)$ ，它的梯度场可以表示为 $\nabla z=\left(\frac{\partial z}{\partial x},\frac{\partial z}{\partial y}\right)=\left(z_x,z_y\right)$ 。最

基本的曲面重建方法是基于线积分，并优化局部代价函数，但是它不能添加任何形式的正则项。Frankot和Chellappa [5] 以及Kovesi [6] 他们分别使用复杂的傅里叶变换和Shapelets，这两种方法都是周期性的，并不现实。Horn和Brooks [7] 提出了变分法，此方法目的在于根据重建曲面上的边界条件最小化最小二乘泛函，后来有很多学者也遵循这一思路。Harker [8] 等人将最小化泛函问题表达为一个半正定的Sylvester方程，提高了求解速度。

笔者借鉴Horn和Harker等人的思想，在最小二乘积分重建技术的基础上引入深度图与深度函数拟合后的Tikhonov正则项，利用L曲线选取最佳参数从而优化了Tikhonov最小二乘模型，提高重建精度。在Sylvester方程下使用Hessenberg-Schur (HS)算法进行求解。最后通过数值实验与现有的算法进行对比。

2. 最小二乘泛函曲面重建

2.1. 最小二乘模型

假设图像是通过正交投影获得的，并且表面具有明确的形式 $z=z\left(x,y\right)$ ，将梯度空间参数细化为

$z_x\left(x,y\right)=\frac{\partial z\left(x,y\right)}{\partial x}$ ， $z_y\left(x,y\right)=\frac{\partial z\left(x,y\right)}{\partial y}$ 。获取高度z即解的前提是有一个被噪声破坏的梯度场。假设噪

声是高斯噪声，最小二乘成本函数的连续形式为 [8] ：

$E\left(z\right)={\displaystyle {\int }_c^d{\displaystyle {\int }_a^b{\left({\hat{z}}_x\left(x,y\right)-z_x\left(x,y\right)\right)}^2}}+{\left({\hat{z}}_y\left(x,y\right)-z_y\left(x,y\right)\right)}^2\text{d}x\text{d}y$ (1)

式中 ${\hat{z}}_x$ 和 ${\hat{z}}_y$ 分别代表测量梯度，这个成本函数是函数差的平方项的体积。利用有限差分的思想，此成本函数在矩形网格上等价的离散形式为：

$E\left(Z\right)={\Vert Z_x-{\hat{Z}}_x\Vert }_{\text{F}}^2+{\Vert Z_y-{\hat{Z}}_y\Vert }_{\text{F}}^2$ (2)

式中 ${\hat{Z}}_x$ 和 ${\hat{Z}}_y$ 表示测量得到的梯度矩阵； ${\Vert \cdot \Vert }_{\text{F}}^2$ 表示F范数。由于微分是一个线性算子，导数可以写成一个简单的矩阵乘法，即线性变换：

$\frac{\partial Z}{\partial x}={\hat{Z}}_x\approx Z{D}_x^{\text{T}},\frac{\partial Z}{\partial y}={\hat{Z}}_y\approx D_{y}Z$ (3)

这里“≈”表示最小平方意义上的相等。 $D_x$ 和 $D_y$ 分别是x和y方向的有限差分矩阵。 $D_x^{\text{T}}$ 会对x方向的微分产生影响。因此，在测量的梯度场 ${\hat{Z}}_x$ 和 ${\hat{Z}}_y$ 中，最小二乘意义上最接近的曲面是以下成本函数的最小值：

$E\left(Z\right)={\Vert Z{D}_x^{\text{T}}-{\hat{Z}}_x\Vert }_{\text{F}}^2+{\Vert D_{y}Z-{\hat{Z}}_y\Vert }_{\text{F}}^2$ (4)

2.2. 带有Tikhonov正则项的最小二乘模型

最常见的正则化项是解的范数的一个界限 [9] 。在离散意义上，解的F范数为：

${\rho }_0\left(Z\right)={\Vert Z\Vert }_F^2=\frac{1}{2}\left({\Vert I_{m}Z\Vert }_F^2+{\Vert Z{I}_n\Vert }_F^2\right)$ (5)

它是一个有效的0阶正则化项。用最右端的形式书写，方程(5)就满足Sylvester方程的框架。其中 $I_m$ 和 $I_n$ 是单位矩阵。对于1阶正则化项而言，它可以利用函数沿单位向量 $u={\left[\begin{array}{cc}{u}_x& u_y\end{array}\right]}^{\text{T}}$ 的方向导数推导出曲面上某一点的最大方向导数，进而得到一个限定曲面陡度的正则化项，其离散化为：

${\rho }_1\left(Z\right)={\Vert D_{y}Z\Vert }_{\text{F}}^2+{\Vert Z{D}_x^{\text{T}}\Vert }_{\text{F}}^2$ (6)

类似的可以导出一个二阶正则化项，它是曲面平均曲率的边界：

${\rho }_2\left(Z\right)={\Vert D_y^{2}Z\Vert }_{\text{F}}^2+{\Vert Z{\left(D_x^2\right)}^{\text{T}}\Vert }_{\text{F}}^2$ (7)

对于从梯度重建曲面的问题而言，要最小化的函数取决于2D曲面而不是向量。所以Tikhonov正则化的常用方法不直接适用。为了推导出2D域的适当函数，从方程(4)开始，通过比较方程(5)~(7)，将矩阵 $S_x$ 和 $S_y$ 分别定义为x和y方向上的一般“平滑”算子。因此，Tikhonov最小二乘泛函的一般形式如下：

$E\left(Z\right)={\Vert Z{D}_x^{\text{T}}-{\hat{Z}}_x\Vert }_{\text{F}}^2+{\Vert D_{y}Z-{\hat{Z}}_y\Vert }_{\text{F}}^2+{\lambda }_{}^2\left({\Vert S_y\left(Z-Z_0\right)\Vert }_{\text{F}}^2+{\Vert \left(Z-Z_0\right)S_x^{\text{T}}\Vert }_{\text{F}}^2\right)$ (8)

这里 $Z_0$ 是由三角测量得到的粗糙深度图，将其作为表面的先验估计。然而将先验估计值纳入算法中的事实可能会对光度立体技术在应用过程中产生巨大的影响。下一节将在L曲线框架中考虑正则化参数的选择，以此来优化Tikhonov最小二乘模型。

3. 模型优化与求解

3.1. 正则化参数的选择

L曲线 [10] 是在未知其他噪声特性的情况下选择正则化参数的最简单方法之一。L曲线是 $\left(\rho \left(\lambda \right),\eta \left(\lambda \right)\right)$ 的曲线，其中 ${\rho }^2\left(\lambda \right)$ 是最小平方成本函数：

${\rho }^2\left(\lambda \right)={\Vert Z\left(\lambda \right)D_x^{\text{T}}-{\hat{Z}}_x\Vert }_{\text{F}}^2+{\Vert D_{y}Z\left(\lambda \right)-{\hat{Z}}_y\Vert }_{\text{F}}^2$ (9)

${\eta }^2\left(\lambda \right)$ 是标准形式的正则化项：

${\eta }^2\left(\lambda \right)={\Vert Z\left(\lambda \right)\Vert }_{\text{F}}^2$ (10)

因此，L曲线是最小二乘残差和正则化项大小相互作用的可视化结果。将x方向和y方向导数算子的奇异值分解(SVD)表示为

$D_x=U_x{L}_x{V}_x^{\text{T}}$ (11)

$D_y=U_y{L}_y{V}_y^{\text{T}}$ (12)

一旦计算了导数矩阵的SVD，L曲线上的点可以计算为：

${\rho }^2\left(\lambda \right)={\Vert M\left(\lambda \right)L_x-Q\Vert }_{\text{F}}^2+{\Vert L_{y}M\left(\lambda \right)-P\Vert }_{\text{F}}^2$ (13)

以及

${\eta }^2\left(\lambda \right)={\Vert M\left(\lambda \right)\Vert }_{\text{F}}^2$ (14)

上述式子中 $M=V_y^{\text{T}}Z{V}_x$ ， $P=U_y^{\text{T}}\hat{Z}{}_{y}V{}_x$ ， $Q=V_y^{\text{T}}\hat{Z}{}_{x}U{}_x$ 。根据F范数在正交变换下的不变性，由于 $L_x$ 和 $L_y$ 是对角矩阵，这些计算成本比SVD的计算成本更小，因此可以通过计算L曲线上的几个点来确定适当的正则化参数，将最佳参数值带入方程(8)中，从而达到优化模型的目的。

3.2. 建立Sylvester方程求解成本函数

若待重建图像有 $m\times n$ 个采样点，则方程(8)中系数矩阵 $D_x^{\text{T}}$ 和 $D_y$ 的大小为 $2mn\times mn$ 。虽然 $D_x^{\text{T}}$ 和 $D_y$ 是稀疏矩阵，空间复杂度较高，但Harker使用Sylvesyer方程将系数矩阵大小降低为 $m\times n$ 。从而本研究将成本函数公式(8)中的参数利用L曲线优化后对 $Z$ 进行微分并使得它的值为0，得到该函数相应的正规方程：

$\left(D_y^{\text{T}}{D}_y+{\lambda }^2{S}_y^{\text{T}}{S}_y\right)Z+Z\left(D_x^{\text{T}}{D}_x+{\lambda }^2{S}_x^{\text{T}}{S}_x\right)-D_y^{\text{T}}{\hat{Z}}_y-{\hat{Z}}_x{D}_x-{\lambda }^2{S}_y^{\text{T}}{S}_y{Z}_0-Z_0\left({\lambda }^2{S}_x^{\text{T}}{S}_x\right)=0$ (15)

此方程是未知曲面 $Z$ 中的Sylvester方程。对于曲面重建问题，本研究考虑“平滑”算子 $S_x=I_n$ ， $S_y=I_m$ 在这种情方程(15)可以改写为：

$\left(D_y^{\text{T}}{D}_y+{\lambda }^2{I}_m\right)Z+Z\left(D_x^{\text{T}}{D}_x+{\lambda }^2{I}_n\right)-D_y^{\text{T}}{\hat{Z}}_y-{\hat{Z}}_x{D}_x-2{\lambda }^2{Z}_0=0$ (16)

方程(16)也是一个Sylvester方程，但由于 $D_y^{\text{T}}{D}_y$ 和 $D_x^{\text{T}}{D}_x$ 有相同的特征值，所以此方程没有唯一解。要想求得唯一解，可以采取HS算法 [8] 求解。首先使用Householder变换将矩阵 $D_x$ 和 $D_y$ 分解为上三角形式的矩阵(11)和(12)，Householder矩阵 $V_y$ 和 $V_x$ 是正交的，不会影响解的数值精度，显然有 $V_y{V}_y^{\text{T}}=I$ ， $V_x{V}_x^{\text{T}}=I$ 。方程(16)可以化简为：

$\left(V_y{L}_y^2{V}_y^{\text{T}}+{\lambda }^2{I}_m\right)Z+Z\left(V_x{L}_x^2{V}_x^{\text{T}}+{\lambda }^2{I}_n\right)-V_y{L}_y{U}_y^{\text{T}}{\hat{Z}}_y-{\hat{Z}}_x{U}_x{L}_x{V}_x^{\text{T}}-2{\lambda }^2{Z}_0=0$ (17)

此方程左乘 $V_y^{\text{T}}$ 右乘 $V_x$ 可得：

$\left(L_y^2+{\lambda }^2\right)V_y^{\text{T}}Z{V}_x+V_y^{\text{T}}Z{V}_x\left(L_x^2+{\lambda }^2\right)-L_y{U}_y^{\text{T}}{\hat{Z}}_y{V}_x-V_y^{\text{T}}{\hat{Z}}_x{U}_x{L}_x-M_0=0$ (18)

其中 $M_0=2{\lambda }^2{V}_y^{\text{T}}{Z}_0{V}_x$ 。

令 $M=V_y^{\text{T}}Z{V}_x$ ， $A=L_y^2+{\lambda }^2$ ， $B=L_x^2+{\lambda }^2$ ， $C=-L_y{U}_y^{\text{T}}{\hat{Z}}_y{V}_x-V_y^{\text{T}}{\hat{Z}}_x{U}_x{L}_x-M_0$ ，方程(18)可以写为：

$AM+MB+C=0$ (19)

此方程可以写为以下形式：

$\left[\begin{array}{cc}0& 0^{\text{T}}\\ 0& A\end{array}\right]\left[\begin{array}{cc}{m}_{00}& m_{01}^{\text{T}}\\ m_{10}& M_{11}\end{array}\right]+\left[\begin{array}{cc}{m}_{00}& m_{01}^{\text{T}}\\ m_{10}& M_{11}\end{array}\right]\left[\begin{array}{cc}0& 0^{\text{T}}\\ 0& B\end{array}\right]+\left[\begin{array}{cc}0& c_{01}^{\text{T}}\\ c_{10}& C\end{array}\right]=\left[\begin{array}{cc}0& 0^{\text{T}}\\ 0& 0\end{array}\right]$ (20)

它等价于下面方程组：

$\{\begin{array}{l}{m}_{00}=0\\ m_{01}^{\text{T}}B+c_{01}^{\text{T}}=0^{\text{T}}\\ A{m}_{10}+c_{10}=0\\ A{M}_{11}+M_{11}B+C=0\end{array}$ (21)

以 $M$ 中元素为未知量进行方程组(21)的求解，将矩阵重新组合，得到待重构的曲面 $Z$ 为：

$Z=V_{y}M{V}_x^{\text{T}}$ (22)

4. 实验结果与分析

为了将更新参数后的结果与现有的方法进行对比，本研究进行了以下实验：

1) 用理想数据合成一个待重建的地形图。利用蒙特卡洛测试，比较成本函数的误差，并将本研究中的重建结果图与现有方法结果图进行对比，分析方法的鲁棒性；

2) 采用蒙特卡洛测试，就不同尺寸图像的计算时间而言，将现有方法的结果与本研究中的结果进行对比。

本研究选取了尺寸大小为 $m\times n=128\times 128$ 的合成地形图，其类似于MATLAB中的“峰值”检测，测试曲面及其梯度场如图1所示。

(a) 测试曲面Z (b) 梯度场Z_x (c) 梯度场Z_y

Figure 1. Topographic map and its gradient field to be reconstructed

图1. 待重建的地形图及其梯度场

将全局最小二乘(Global Least Squares, GLS)以及Dirichlet边界结果与本研究得到的结果进行比较。选择GLS的原因在于，一方面它和本研究一样，都是在Sylvester方程框架下进行求解，另一方面，当梯度场被高斯噪声破坏后，GLS是一个标准的解决方案；选择Dirichlet边界的原因在于，确定曲面的边界可以减少自由度的数量，从而减少未知数的数量，在这种情况下，Dirichlet问题的解是唯一的，同样可以利用Sylvester方程求解，并且对异常值具有极强的鲁棒性。

根据现有文献的方法，将 $\sigma =0.05$ 的高斯噪声添加到梯度场中，当受高斯噪声影响时，解析梯度不再是可积的。图2展示了图1中测试曲面Z被高斯噪声破坏时的梯度。

为了使得从被噪声破坏的梯度场中重建曲面的效果更好，本研究利用L曲线算法，选取100个要计算的L曲线上的点数，从中选取正则化参数 $\lambda$ 的最佳值。 $\lambda$ 的选择结果如图3中(a)所示。从图中可以看出，当 $\lambda =0.2833$ 时为最佳值，优化模型后重建的曲面如图3中(b)所示。图3中的(c)和(d)分别为GLS以及Dirichlet边界重建曲面的效果图，可以看出利用L曲线以及Dirichlet边界重建的曲面抑制噪声以及偏差的效果更好，但是本研究的方法在边界上要比Dirichlet边界更稳定更平滑；GLS以及本研究方法整体效果更优，虽然本研究方法在左侧和右侧边缘处高度值恢复有误差，边缘处有少量锯齿状的不稳定值存在，但是比GLS重建表面更光滑。

(a) 被噪声破坏后的梯度场Z_x (b) 被噪声破坏后的梯度场Z_y

Figure 2. Gradient of test surface Z damaged by Gaussian noise

图2. 测试曲面Z被高斯噪声破坏时的梯度

(a) L曲线 (b) Tikhonov (L曲线) (c) GLS (d) Dirichlet边界

Figure 3. L curve and the reconstructed surface renderings under the three methods

图3. L曲线以及三种方法下的重建曲面效果图

为了更进一步分析三种方法的鲁棒性，采用均方误差(MSE)以及相对均方误差(RMSE)来衡量。

$\text{MSE}:=\frac{1}{mn}{\displaystyle \underset{i=1}{\overset{m}{\sum }}{\displaystyle \underset{j=1}{\overset{n}{\sum }}{\left|{\hat{z}}_{i,j}-z_{i,j}\right|}^2}}$ ， $\text{RMSE}=\frac{\text{MSE}}{{\left|{\displaystyle \underset{i=1}{\overset{m}{\sum }}{\displaystyle \underset{j=1}{\overset{n}{\sum }}{z}_{i,j}}}\right|}^2}$ (23)

其中， ${\hat{z}}_{i,j}$ 表示图像像素点 $\left(i,j\right)$ 处测得的高度， $z_{i,j}$ 表示真实高度。所得MSE和RMSE如表1所示，分别展示了三种方法下的重建误差。

Table 1. Reconstruction error of each algorithm under noise σ = 0.05 and surface size of 128 × 128

表1. 在噪声 $\sigma =0.05$ 下表面尺寸为128 × 128下各算法的重建误差

从表1可以看出，本研究所采用方法比另外两种方法下的误差显著降低，其主要原因在于本研究中的方法通过惩罚曲面的平均曲率来对成本函数起作用，因此是一个低通型滤波器，可以达到消除噪声的目的，鲁棒性增强；GLS受异常值扰动的影响较大，因此鲁棒性弱于本研究中的方法；Dirichlet边界在图像边缘处鲁棒性较弱。

不同尺寸图像的计算时间比较的结果将在表2中展示。

Table 2. The computation time of the algorithm

表2. 算法的计算时间

表2结果显示，通过L曲线优化后的模型在小尺寸图像的重建上有着较强的优势，比另外两种方法快了将近0.1 s；在中等尺寸和大尺寸图像的重建中弱于GLS；虽然大尺寸图像上Dirichlet边界优势显著，但本研究仍具有较好的应用价值。

5. 结语

笔者通过引入L曲线来选取参数 $\lambda$ 的最优值，从而优化了Tikhonov最小二乘模型。利用有限差分思想将模型进行离散，使得离散后的解满足Sylvester方程。利用HS算法求解降低计算时间。仿真实验表明，在理想情况下，优化模型后的方法能够有效抑制噪声和偏差，在小尺寸图像上比五点差分格式的GLS耗时更短，误差更小；优化模型后重建出的三维图精度更高，更接近原始曲面，整体效果更优。未来考虑具有边界条件的Tikhonov正则化，设计出更快速有效的数值方法，进一步缩短重建时间，提高重建精度和质量。

基金项目

2022年湖南省教育厅科学研究项目(22A0368)；吉首大学校级科研项目(Jdy22012)。

NOTES

^*通讯作者。