企业官网建设流程全解析

做触觉传感器和三维形貌重建的人，多半绕不开 GelSight Sensor 和它输出的 height-map。我最早接触这套传感器时，天真地以为拿到表面法向图之后，高度场无非是沿着像素路径积一下分，拉一条跨平台式的 scanline 就能搞定。直到真的在真实数据上跑了一轮，才发现噪声、边界和累积误差会把结果毁得没法看。后来我把问题改写成泊松方程，配合 DST-II 算子做快速求解，计算量从“跑一次要等到怀疑人生”降到“实时预览毫无压力”。这篇文章就讲讲这套方案的思路、实现和踩坑记录，适合正在做触觉重建、机器人抓取感知，或者想了解快速数值求解如何落到图像处理管线的朋友。

1. 先搞清楚我们要解什么问题

1.1 GelSight 怎么把“触觉”变成“图像”

GelSight 的物理结构并不神秘：一块软凝胶覆盖在透明的支撑板上，凝胶表面喷涂了反光涂层。当物体压上去，凝胶表面跟着物体形状发生微米级的形变，像一面镜子。传感器内部用红、绿、蓝三组 LED 从不同方向照射这面“镜子”，顶上再放一个相机拍下反射光。形变越大、局部倾斜越陡，反射光进入相机的角度就变化越明显，RGB 三个通道的亮度就会呈现规律性的变化。

这背后的光学逻辑是：局部表面法向决定了反光方向，反光方向决定了相机能不能在某个 LED 的照射下看到高光。把三个通道的亮度分别理解成“来自三个方向的光照强度”，就能通过比色法反推出每个像素处表面的倾斜分量。换句话说，GelSight 直接测量的是一个梯度场，而不是高度本身。它输出的原始信号，经过标定后可以转成每个像素处的法向量或者是表面的斜率分量 gx、gy。

这个阶段最重要的标记是：我们要重建的 height-map 不是直接从图像里读出来的，它像一个“积分结果”。传感器给的是导数的证据，我们要从证据里还原原函数。于是问题天然落到数值计算上。

1.2 height-map 的直接积分法为什么不行

我第一次上手时用的是最直观的路径积分：从图像左上角出发，gx 是 x 方向的斜率，gy 是 y 方向的斜率，于是逐像素累加：

h[i][j+1] = h[i][j] + gx[i][j] h[i+1][j] = h[i][j] + gy[i][j]

听起来没有任何问题，可真实数据一跑就露馅了。第一，GelSight 图像里的梯度场带有大量噪声，每个像素的 gx、gy 都有小误差，累加时误差被一路放大，积分久了画面里全是条纹状的漂移。第二，路径选择会直接影响结果：从左上角积分得到的右下角高度，和从右下角反着积分回来的数值完全对不上，因为噪声让梯度场不再是一个“可积”的保守场。这就像你在一个起伏不平的地形图上沿着不同路线量海拔，每条路线量出来的终点高度都不一样，你会怀疑到底是地形有问题还是尺子有问题。

更本质的原因是路径积分只利用了局部信息。每一个点的高度只依赖于它同一条路径上的前一个点，它不知道侧向的梯度值和周围一大片区域传递过来的约束。在存在噪声时，这种“一条路走到黑”的算法几乎是灾难。所以行业里很少直接用路径积分来做 GelSight 的形貌重建，除非只是看一个粗略的剖面趋势。

1.3 泊松方程：把局部梯度场“拼”回全局高度场

既然路径积分不靠谱，那就换成全局优化思路。我们希望找到一个高度场 h，使它的梯度 ∇h 在最小二乘意义下尽量接近传感器测出来的梯度场 g。也就是：

min ∫ || ∇h - g ||² dx dy

对这个目标函数做变分，等价于求解一个泊松方程：

∇²h = div(g)

这里 ∇² 是拉普拉斯算子，div 是散度。右边 div(g) 完全由实测梯度场决定：对 gx 求 x 方向差分，加上 gy 求 y 方向差分。左边是标准二阶微分算子。整个方程的含义非常直观：高度场的“弯曲程度”应该等于实测梯度场的“膨胀程度”，而不再是逐点累加。

这个形式的好处是它在全局范围内分配误差，不会出现路径积分那样的沿路径漂移。即使某个局部梯度值带了噪声，它只会影响该点附近的平滑度，而不会顺着一条线一路污染到图像边缘。用图像处理的术语来说，泊松方程像一个全局低通滤波器，把不可积的梯度场“投射”到可积函数空间里，投影之后剩下的残差就是噪声。

所以剩下的核心问题只有一个：如何快速、稳定地求解 ∇²h = div(g)，尤其是当图像分辨率来到百万像素级别的时候。

2. 为什么选 DST-II 算子而不是硬解矩阵

2.1 从拉普拉斯算子到线性方程组

离散化之后，拉普拉斯算子变成经典的卷积模板。二维情况用五点模板：

∇²h[i][j] ≈ (h[i+1][j] + h[i-1][j] + h[i][j+1] + h[i][j-1] - 4h[i][j]) / Δ²

把整个 height-map 摊成一个一维向量，这个离散拉普拉斯就变成一个巨大的稀疏矩阵 A。泊松方程变成：

A h = b

b 就是散度项 div(g)。如果不考虑效率，直接调一个稀疏线性求解器也能跑。但问题是：当图像尺寸是 1024 x 1024 时，未知量超过一百万，即使 A 是稀疏矩阵，通用求解器也要吃大量内存和时间。在实际项目里，我经常要在一台工控机上做实时重建，这种“通用但昂贵”的算法肯定没法用。

换个角度想：A 是五对角块矩阵，它对应的离散泊松方程是一个结构极其规则的线性系统。对规则系统，数学上早就有比通用求解器快几个数量级的方法——换到频率域去解。

2.2 特征分解思维：换到“频率域”去解方程

回忆一下信号处理里解微分方程的常规操作：对时间信号做傅里叶变换，微分算子就变成乘以 jω，于是微分方程变成一个代数方程。微分算子在频域里是对角的，这个性质简直是为快速求解量身定做的。

离散泊松方程也一样。如果我们能找到一组基，使得离散拉普拉斯矩阵 A 在这组基下变成对角矩阵，那么解方程就变成三步：

1. 把 b 变换到频域：b_hat = transform(b) 2. 在频域做逐点除法：h_hat = b_hat / lambda 3. 逆变换回空间域：h = inverse_transform(h_hat)

整个流程里没有迭代，没有矩阵求逆，只有三次变换加一次逐元素除法。如果变换本身是 O(N log N) 的快速算法，那整体求解速度比通用稀疏求解器快几个数量级，尤其在高分辨率下优势越发明显。

这其实就是“用特征分解理解快速求解”的核心思路。A 的维度再大，只要它能被一组固定基同时对角化，我们就可以绕开显式的矩阵存储和分解。离散余弦变换（DCT）和离散正弦变换（DST）之所以在图像编码里随处可见，本质上也是因为它们能对角化某些特殊边界的二阶微分矩阵。

2.3 DST-II 到底做了什么，边界条件如何对号入座

具体到 GelSight 重建，我为什么不用更常见的 DCT-II，而选了 DST-II？这里有一个特别容易被忽略的细节：边界条件。

在标准的二维泊松重建里，如果高度场四周都是自由边界，也就是法向导数在边缘处为零，最自然的选择是 DCT-II。DCT 的基函数是余弦函数，它自带“边界外推平坦”的 Neumann 特征。但 GelSight 重建 height-map 时，很多场景下我们并不是在整个二维平面上解方程，而是沿着扫描方向一条剖面、一条剖面地重建。比如机器人按压时只关心接触中心区域的高度轮廓，或者我们用滚筒式扫描接触痕迹，每列像素都是一条独立的轮廓线。

对每条剖面线，边界条件通常是一端固定、一端自由。固定端可以取接触区域的起点，自由端是远离接触的位置。这种混合边界的离散二阶微分矩阵，特征函数正好是 DST-II 的基函数，也就是：

sin(π * (n + 0.5) * (k + 1) / N), n, k = 0, 1, ..., N-1

这个公式不要被吓到，它可以理解为：在像素中心位置采样、但边界上满足一硬一软条件的正弦基。如果你实际做过数值验证就会发现，把离散拉普拉斯矩阵的特征向量画出来，它们和 DST-II 的基函数逐点重合。而 Gradient 和散度在 staggered grid 上的定义方式，也天然让 DST-II 成为对界最准的那个算子。

二维整面重建当然也可以用 DCT-II 一类的方法，但如果要处理的是大量一维剖面，DST-II 是最顺手的方案。用 DST-II 的好处是边界条件和物理模型完全一致，不用人为在图像外侧补一堆假数据来凑 Neumann 边界。

3. 实现流程：从 GelSight 图像到 height-map

3.1 图像采集与梯度场提取

在实际系统里，GelSight 的原始输出是一张彩色图像。以常见的 RGB 三光源配置为例，红、绿、蓝三个通道各自对应一个方向的照射响应。拿到图像后的第一步不是直接算梯度，而是先做预处理。

我习惯先把三个通道拆开。这个操作在不同视觉库里有不同叫法，在 OpenCV 里是split，在类 Halcon 的更底层操作里可以理解为一个通道切片（channel slice）算子，把三通道图按波段切成三个单通道图。拆完通道后要做一次亮度归一化，把白平衡偏差去掉。凝胶表面反射涂层如果不是完全均匀，图像里会出现低频明暗不均，这一步能明显缓解后续梯度场的系统性倾斜。

接下来用一阶微分算子提取梯度。最常用的是 Sobel 算子，它的原理是对图像做带权重的差分，同时起平滑作用。X 方向的 Sobel 核提取水平梯度，Y 方向的 Sobel 核提取垂直梯度。对 GelSight 这种表面纹理较弱的图像，Sobel 十分合适，因为它对噪声有一定抑制，不会像单纯的前向差分那样把像素噪声全部放大。

import cv2 import numpy as np img = cv2.imread("gelsight_sample.png") b, g, r = cv2.split(img) # 用三个通道的亮度关系换算出表面梯度分量 # 具体换算矩阵由 LED 方向标定得到，这里用近似系数示意 gx = 0.7 * (r.astype(np.float32) - g.astype(np.float32)) gy = 0.7 * (g.astype(np.float32) - b.astype(np.float32))

这段代码是一个高度简化的示意。严格的做法是先在标定块上采集已知平面和若干已知曲率目标的光照响应，建立“通道亮度到局部倾斜角度”的查找表或线性映射矩阵。做一次完整标定之后，GelSight 的测量精度可以到微米级。

3.2 构建散度项与边界条件

有了 gx 和 gy 之后，下一步是把它们合成泊松方程的右端项 b = div(g)。离散散度用的是差分格式。为了和 DST-II 的 staggered grid 配合，我通常把梯度定义在半像素位置，散度则用相邻梯度之差来近似：

div(g)[i][j] = (gx[i][j] - gx[i-1][j]) + (gy[i][j] - gy[i-1][j])

边界处则根据具体模型做处理。如果是逐列剖面重建，每一列的起点高度直接固定为 0，这是 Dirichlet 条件；剖面末端认为法向倾角为 0，这是 Neumann 条件。这种一边固定、一边自由的组合，正好和 DST-II 的特征函数匹配。

如果做整面二维重建，边界处理会更繁琐。我的习惯是先把图像边缘外扩 4 个像素，外扩方式用镜像对称，然后再算散度。这样能显著减少重建出的高度图边缘出现的“碗状”扭曲，也就是边界伪影。

3.3 DST-II 快速求解步骤与代码级示意

求解阶段用 DST-II 变换。Python 生态里没有像 FFT 那样开箱即用的标准 DST-II 封装，但可以用 FFTW 的 RODFT10 变换类型，或者自己在 NumPy 里手写一个基于 FFT 的快速实现。

一维剖面重建的核心逻辑如下：

def dst2_1d(x, axis=0): # DST-II 可以借助 FFT 实现，这里为可读性直接用矩阵乘法示意 N = x.shape[axis] n = np.arange(N) k = np.arange(N).reshape(-1, 1) basis = np.sin(np.pi * (n + 0.5) * (k + 1) / N) return basis @ x # 假设 gx 是某一行剖面的 x 方向梯度 N = gx.shape[0] b = np.zeros(N) b[1:-1] = gx[1:] - gx[:-1] # 内部散度 b[0] = gx[0] # 起点固定 h[0] = 0，散度单独处理 # 变换到频域 b_hat = dst2_1d(b) # 频率域逐点除以特征值 # 特征值与离散拉普拉斯矩阵和边界条件绑定，实际建议数值验证后使用 k = np.arange(N) lam = 2 * np.cos(np.pi * (k + 1) / N) - 2 lam[0] = 1.0 # 避免除零，直流分量在固定起点场景下无意义 h_hat = b_hat / lam # 逆变换（DST-III 或 DST-II 的逆，注意归一化系数） h = np.zeros(N) # 实际工程中建议用 FFTW 的 RODFT10 和 RODFT01 配对 # 这里的矩阵版本省略了归一化，作示意用

实际工程项目里我不会建议手写这个变换。成熟的方案是：

FFTW 的 RODFT10 做正变换，RODFT01 做逆变换

FFTW 里这两个变换类型就是标准 DST-II 和它的逆。如果你用的是基于 FFTW 封装的 Python 库，例如pyfftw，可以直接调用。配合多线程，1080p 图像的逐列剖面重建能在几十毫秒内完成。

二维整面重建也可以照搬这个思路：先在 x 方向做一次 DST-II，再在 y 方向做一次 DST-II，散度项用二维差分构造，最后在频域做一次二维逐元素除法。边界条件需要和变换类型严格匹配，实际操作时我用一维剖面版本做了一个 small test，再扩展到二维。

3.4 可视化与灰度拉伸

重建出的 height-map 是浮点数据，数值范围往往非常小。比如一个深度 50 微米的微小压痕，它的高度值在 -2.5e-5 到 2.5e-5 之间浮动，直接换算成 8 位图就是一片漆黑。这时候必须要做灰度拉伸。

在 Halcon 里，我经常用灰度值拉伸算子配合min_max_gray统计图像中的极值，再把整个深度范围线性映射到 0 到 255。OpenCV 的话可以用normalize配合 percentile 截断，把 0.5% 到 99.5% 的数值映射到全灰度范围，防止个别极端噪声点拖垮整体对比度。

除了可视化，灰度拉伸还能辅助后续的特征提取。比如做缺陷检测时，经过拉伸后微小划痕变成清晰的亮暗纹理，再用阈值分割或者边缘检测就方便多了。

4. 实操避坑与常见问题

4.1 边界伪影与条带噪声

第一个坑是边界伪影。用 DCT/DST 类变换求解泊松方程时，边界条件如果和变换的特征结构不匹配，会看到重建结果在图像边缘出现明显的“拱起”或“下陷”。尤其是二维重建时，如果四周都当自由边界处理，但实际接触区域外 GelSight 的梯度场并不为零，重建出的高度图边缘就会像碗边一样翘起来。

我常用的解决办法是外部扩展。在进入散度计算前，先把原始梯度图向外侧镜像扩展 8 到 16 个像素，求解完成后再把扩展部分裁切掉。这样外扩区域承担了边界效应，目标区域的边缘就干净很多。条带噪声则通常来自光源不均匀或者反射涂层缺陷，处理思路是先对 gx、gy 做沿扫描方向的中值滤波，再用灰度拉伸增强对比。

4.2 标定偏差与光源不均

GelSight 的梯度场是从三通道亮度换算过来的，换算矩阵的准确度直接决定 height-map 的精度。我在早期试验时发现重建出的球面轮廓总是带一个莫名其妙的“锥形偏置”，排查了很久才发现是红色 LED 的亮度衰减比蓝色 LED 快，导致换算系数在高亮度区域出现了非线性偏移。

解决方法是做多点标定：用一堆已知曲率半径的标准球或标准斜块，在不同亮度区间采集样本，拟合分段线性或多项式映射。另外，光源老化会导致标定参数漂移，所以项目上线后最好定期重新标定。如果环境光干扰严重，可以考虑在传感器外面加遮光罩，或者在算法里加一个暗帧扣除：盖上不透光盖板拍一张全黑图，后续每帧都减去这张暗帧，能有效消除固定模式噪声。

4.3 高度场漂移与倾斜修正

即使用了泊松方程，重建出的高度场有时还是会带整体倾斜项。原因是传感器本身的相机光轴和凝胶表面法向不完全平行，导致梯度场里混入了一个恒定偏置。

解决方式是在重建后做一个平面拟合，把拟合得到的倾斜平面从高度场里减去：

def remove_tilt(h): H, W = h.shape y, x = np.mgrid[0:H, 0:W] A = np.stack([x.ravel(), y.ravel(), np.ones_like(x.ravel())], axis=1) coef, _, _, _ = np.linalg.lstsq(A, h.ravel(), rcond=None) tilt = (coef[0] * x + coef[1] * y + coef[2]).reshape(H, W) return h - tilt

这个操作在离线分析和实时处理中都很常见。如果高度场还伴随低频波浪，一般是因为凝胶表面自身不是绝对平面，这时候可以用高通滤波或者多帧平均背景相减来处理。

4.4 性能与硬件实现：算子优化方向

整套 GelSight 重建流程实际是由一连串算子组成的：拆分通道、灰度归一化、Sobel 梯度提取、散度计算、DST-II 正变换、频域除法、DST-II 逆变换、灰度拉伸。在桌面处理器上跑，1080p 分辨率下面向实时毫无压力，但一旦要把这套能力部署到机器人夹爪、嵌入式视觉模块里，硬件压力马上就来了。

大量的算子顺序执行，每一遍都涉及对整幅图像的内存读写，瓶颈通常不在算力而在内存带宽。所以优化第一优先级是算子融合：把 Sobel 梯度提取和散度计算合并成一次卷积，省一趟内存读写；把拆分通道和白平衡合并到相机原始数据到 RGB 的 ISP 流程里；频域除法和逆变换的前半部分可以在同一个循环里完成。

我做过一个嵌入式平台上的 DST-II 实现，发现手写蝶形运算时主要性能杀手是三角函数查表带来的 cache miss。把角度预计算成 LUT 后整体算子快了近 30%。如果再往底层走，可以用定点整数近似替代浮点运算，不过要注意动态范围，建议先仿真验证误差。DST-II 本身是一个高度规整的移位加乘结构，很适合在 GPU 或 DSP 上用 SIMD 指令加速。这也是现在算子开发工作里很常见的优化方向：从数学形态不变的算法等价变换开始，再针对目标硬件做手脚。

最后再分享一个小技巧。DST-II 的特征值在不同库里的排列顺序和归一化系数差异很大，第一次写代码时不要直接套论文公式，先用一个已知的小矩阵做数值验证：随手构造一个小尺寸测试高度图，算出它的梯度场，再用你的 DST-II 求解器重建回来，看误差是不是在浮点精度范围内。我每次换平台重写这个算法都会先跑一遍这个小测试，能省掉后面排查一堆莫名其妙问题的功夫。这套方法后续还可以扩展，比如用在柔性触觉阵列的非接触标定、压痕深度统计、或者机器人滑移检测里的微几何变化追踪上，核心都是同一套“梯度场到高度场”的快速求解逻辑。