# 第6章 本体感觉与触觉感知

> 原课程: Robots That Learn, UC Berkeley CS294-291, Spring 2026
> 主讲: Jitendra Malik
> 本章基于 Lecture 4B: Proprioception and Tactile Sensing, 阅读论文 Jones (2006) + Gardner (2010)

---

## 一、本章概要

闭上你的眼睛。用右手从口袋里摸出一枚硬币。你的拇指和食指能瞬间判断它的形状（圆形）、纹理（光滑但有凹凸纹路）、温度（低于体温）、硬度（金属的刚性），以及它的边缘曲率。在此过程中，你从未"看"到硬币，但你的大脑构建了关于这个物体的完整表征。

这个看似简单的动作背后，是一套精密的信息处理管道：17,000 条机械感受器纤维在指尖的表皮下排列成精确的阵列，它们的动作电位以毫秒级的精度沿 30,000 条传入神经向上传递，穿过脊髓背柱、延髓楔束核、丘脑腹后外侧核，最终抵达初级体感皮层（S-I）的 Brodmann 3b 区。与此同时，来自肌肉纺锤体和高尔基腱器的本体感觉信号告诉你每根手指的关节角度和施力大小。你的运动皮层发出前向模型预测动作的体感后果，并将这些预测与实际的体感反馈进行实时比较。这就是**主动触觉**（active touch）——不只是"感受"，而是"通过运动来感受"。

对机器人而言，触觉和本体感觉是极其困难的问题。尽管计算机视觉在过去十年取得了革命性突破（ImageNet 精度超越人类、Segment Anything 实现通用分割、扩散模型生成照片级图像），但触觉传感器至今没有一个类似于"摄像头"的通用标准化硬件。GelSight 和 BioTac 等传感器正在填补这一空白，但它们的分辨率、带宽、耐用性和成本距离生物触觉系统仍有数量级的差距。

本章是课程三部曲"生物灵感"部分的收官章，为后续的控制（第7-9章）和灵巧操作（第16章）奠定感知基础。我们将从四个层次解构触觉系统：**（1）生物层次**：机械感受器的类型、分布和神经编码机制；**（2）物理层次**：指尖皮肤的力-位移本构模型和振动传递函数；**（3）工程层次**：GelSight 和 BioTac 等触觉传感器的工作原理与物理模型；**（4）算法层次**：触觉信号处理的 Python 实现。

本章所依赖的核心阅读论文为 Gardner (2010) 的《Touch》综述（全面覆盖机械感受器分类、体感通路、皮层处理）和 Jones & Lederman (2006) 的《Human Hand Function》（指尖力学模型、主动触觉感知、灵巧操作中的体感反馈）。两者共同构成了从"感受器电位"到"有意识的触觉感知"的完整链条。

---

## 二、核心概念与数学基础

### 2.1 四类机械感受器的功能分类

人类无毛皮肤（glabrous skin）中存在四类机械感受器-传入纤维单元，其分类基于两个核心标准：**适应速度**（对持续压痕的响应衰减速率）和**感受野大小**（图 2.1）。

| 类型 | 受体末梢 | 适应特征 | 感受野直径 | 指尖密度 | 最优频率 | 知觉属性 |
|------|---------|---------|-----------|---------|---------|---------|
| **SA-I** | Merkel 细胞-神经突复合体 | 慢适应、不规则放电 | 2-3 mm | ~100/cm² | 0.4-3 Hz | 压力、形状、纹理 |
| **SA-II** | Ruffini 末梢 | 慢适应、规则放电 | ~10 mm | ~10/cm² | 100-500 Hz | 皮肤拉伸、手指位置 |
| **FA-I** | Meissner 小体 | 快适应 | 3-5 mm | ~150/cm² | 3-40 Hz | 颤动(flutter)、滑移检测 |
| **FA-II** | Pacinian 小体 | 快适应、极敏感 | >20 mm | ~350/手指 | 40-500+ Hz | 振动、工具接触 |

**(本章表格资料来源: Gardner 2010 §Primary Afferent Terminals; Johansson & Vallbo 1983; Bolanowski et al. 1988)**

**SA-I（慢适应 I 型）**位于指纹脊的中央，此处皮肤弹性最大。SA-I 感受器是最小的触觉感受器，其放电频率与施加于皮肤的正压力成正比（Goodwin & Wheat 2004; Johnson 2001）。它们的关键功能是**空间编码**——由于单个 Merkel 细胞仅覆盖约 0.5 mm 的区域，SA-I 阵列可以精确解析物体表面的曲率、边缘方向和 Braille 点阵的空间布局。SA-I 传入纤维对垂直于皮肤表面的压痕产生近似线性的响应（直线范围约 0-1500 μm，Blake et al. 1997），这使其成为形状感知的主要神经基质。

**FA-I（快适应 I 型）**由 Meissner 小体介导，位于指纹脊的侧缘，通过胶原纤维与脊缘连接。这一解剖位置使其对切向剪切力极其敏感（Iggo & Andres 1982; Nolano et al. 2003）。当一个被抓握的物体开始滑动时，FA-I 传入纤维以高频率爆发动作电位，触发反射性的握力增加以阻止物体滑落。FA-I 纤维因此被称为"滑移检测器"（Johansson 1996）。它们不响应稳态压力。基于心理物理数据，FA-I 通道在 3-40 Hz 的频率范围内最优激活，产生的知觉体验被称为"颤动"（flutter, Talbot et al. 1968）。

**FA-II（快适应 II 型）**由 Pacinian 小体介导——一种长约 1 mm 的洋葱状多层膜囊结构，位于皮下组织和肠系膜中。Pacinian 小体的膜层结构充当**机械高通滤波器**：当振动频率为 200 Hz 时，膜层刚度不足以及时跟随探针的运动，因此全部膜层同相振动，有效压缩传入神经末梢。而在 20 Hz 的低频刺激下，上层膜被挤压在一起，流体横向移位，底层膜保持静止，神经末梢不受刺激（Bell et al. 1994）。这一物理机制赋予 Pacinian 小体对 100-400 Hz 振动的高灵敏度（最小值约 1 μm 振幅 at 250 Hz），解释了为什么我们能通过手握感知电脑硬盘的马达嗡嗡声。

**SA-II（慢适应 II 型）**与 Ruffini 末梢相关。与 SA-I 不同，SA-II 的感受野约大五倍，但对皮肤形变的敏感度约低六倍。SA-II 的特殊之处在于它对**皮肤侧向拉伸**的方向选择性，尤其对近端方向的拉伸敏感（Birznieks et al. 2001）。由于手背皮肤的松弛性（手指从伸展到完全弯曲，背侧皮肤可被拉伸多达 30 mm），SA-II 传入纤维在手指运动时被可靠地激活——它既非纯粹的"触觉"感受器，也非纯粹的"本体感觉"感受器，而是**模糊了两者的边界**（Edin & Johansson 1995）。

### 2.2 机械感受器传递函数的数学模型

一个完整的机械感受器模型需要描述三条串联的变换链：

**（1）力学前滤波（Mechanical Pre-Filtering）**

皮肤和感受器囊膜构成一个非均匀、各向异性、粘弹性的力学滤波器。设皮肤表面的位移场为 $$\mathbf{u}(\mathbf{x}, t)$$（其中 $$\mathbf{x} \in \mathbb{R}^3$$ 为皮肤表面坐标），感受器神经末梢处的组织应变为 $$\varepsilon(\mathbf{x}_r, t)$$，则传递关系为：

$$P_{\text{mech}}(\mathbf{x}, \omega) = \frac{\tilde{\varepsilon}(\mathbf{x}_r, \omega)}{\tilde{u}(\mathbf{x}_{\text{surf}}, \omega)} = H_{\text{skin}}(\omega) \cdot H_{\text{capsule}}(\omega) \cdot G(\mathbf{x}_{\text{surf}}, \mathbf{x}_r)$$

其中 $$\tilde{u}$$ 和 $$\tilde{\varepsilon}$$ 是傅里叶变换后的皮肤表面位移和末梢应变，$$H_{\text{skin}}(\omega)$$ 是皮肤层的频率响应，$$H_{\text{capsule}}(\omega)$$ 是感受器囊膜的频率响应，$$G(\mathbf{x}_{\text{surf}}, \mathbf{x}_r)$$ 是空间格林函数（描述表面某点的位移如何传播到感受器位置）。

对于 **Pacinian 小体**，其多层膜结构可用 Loewenstein & Skalak (1966) 的层合流体力学的简化模型描述：

$$H_{\text{PC}}(s) = \frac{s^2}{(s + \omega_L)(s + \omega_H)}$$

其中 $$s = j\omega$$，$$\omega_L \approx 2\pi \times 40$$ rad/s（低频截止），$$\omega_H \approx 2\pi \times 500$$ rad/s（高频衰减）。这是一个**带通**特性，但低频衰减非常陡峭（因为在低频下各膜层的滑动解耦了力的传递），实际上表现为偏重高频的带通。

**（2）机械电转导（Mechanoelectrical Transduction）**

神经末梢膜上的机械敏感离子通道将膜应变转化为离子电流。目前最被广泛接受的假说是**栓绳门控模型**（tether-gating model, Christensen & Corey 2007）：通道蛋白通过结构蛋白（如胶原纤维）与细胞外基质连接，通过肌动蛋白丝与细胞骨架锚定。当组织形变拉伸栓绳时，通道打开，Na⁺ 和 Ca²⁺ 内流，产生**感受器电位**（receptor potential）。

感受器电位的幅值与施加压力成正比：

$$V_{\text{receptor}}(t) \propto \kappa \cdot F_{\text{applied}}(t) \ast h_{\text{adapt}}(t)$$

其中 $$\kappa$$ 是压电转导系数（单位为 mV/mN），$$h_{\text{adapt}}(t)$$ 是适应脉冲响应函数。对于 RA 型感受器，$$h_{\text{adapt}}(t)$$ 包含一个快速衰减项（时间常数 $$\tau_{\text{fast}} \approx 5-20$$ ms），使得感受器电位在持续刺激下迅速回归基线。

适应动力学的核心微分方程为：

$$\tau_{\text{adapt}} \frac{dV_{\text{receptor}}}{dt} + V_{\text{receptor}} = \kappa \cdot F_{\text{applied}}(t)$$

其中 $$\tau_{\text{adapt}}$$ 区分了 SA ($$\tau_{\text{adapt}} > 1$$ s) 和 RA ($$\tau_{\text{adapt}} \approx 10-50$$ ms) 类型。

**（3）尖峰编码（Spike Encoding）**

当感受器电位超过阈值 $$V_{\text{th}}$$ 时，电压门控钠通道开放，触发动作电位。动作电位的时序——而非幅度——携带信息。这一过程可用**泄漏积分发放模型**（Leaky Integrate-and-Fire, LIF）建模：

$$\tau_m \frac{dV_m}{dt} = -(V_m - V_{\text{rest}}) + R_m I_{\text{receptor}}(t)$$

发放条件：当 $$V_m \geq V_{\text{th}}$$ 时，产生一个尖峰，并重置 $$V_m \to V_{\text{reset}}$$，进入不应期 $$\tau_{\text{ref}}$$。

由此，机械感受器的完整传递函数可写为级联形式：

$$\text{F}(t) \xrightarrow{\text{力学滤波}} \varepsilon(t) \xrightarrow{\text{转导}} V_{\text{receptor}}(t) \xrightarrow{\text{编码}} \{t_1, t_2, \ldots, t_n\}$$

对于 SA-I 感受器，在稳态下的输入-输出关系可由**压缩幂函数**近似（Johnson 2001）：

$$r_{\text{SA-I}}(F) = r_0 + a \cdot (F - F_{\text{th}})^n \cdot \Theta(F - F_{\text{th}})$$

其中典型值 $$n \approx 0.5-0.7$$（亚线性压缩，暗示 SA-I 编码的不是绝对力而是力的相对变化），$$F_{\text{th}} \approx 1-5$$ mN 为力阈值，$$\Theta(\cdot)$$ 为阶跃函数。

### 2.3 神经编码：从脉冲序列到感知

**（a）频率编码 vs. 时序编码**

长期以来，触觉神经编码的研究存在两种竞争性观点。**频率编码**（rate coding）假说认为：感知强度由发放率（spikes per second）编码。**时序编码**（temporal coding）假说认为：精确的尖峰时间（精确到毫秒）携带了比平均发放率更多的信息。

对于 SA-I 感受器，实验证据强烈支持**空间-频率联合编码**：纹理的精细空间结构由哪些 SA-I 感受器被激活（空间编码，space code）来表征，而纹理的粗糙度由 SA-I 和 FA-I 的总体发放率（频率编码，rate code）来表征。Connor et al. (1990) 通过心理物理实验证明，感知到的表面粗糙度 $$R_{\text{perceived}}$$ 与 SA-I 群体发放率的空间变异度高度相关：

$$R_{\text{perceived}} \propto \text{Var}_{\mathbf{x}}\left[r_{\text{SA-I}}(\mathbf{x})\right]$$

其中 $$\mathbf{x}$$ 遍历激活的感受野位置。

**（b）群体编码与空间模式**

单个感受器提供的信息是局部的。物体形状和纹理的整体知觉来自于**群体编码**（population coding）。设在一个指尖上有 $$N \approx 250-300$$ 条机械感受纤维，每条纤维在时刻 $$t$$ 的发放率为 $$r_i(t)$$。刺激模式 $$\mathcal{S}$$（例如一个 Braille 字母）被编码为一个 $$N$$ 维发放率向量：

$$\mathbf{r}(t; \mathcal{S}) = [r_1(t), r_2(t), \ldots, r_N(t)]^\top$$

大脑通过将这一高维向量映射到特定的感知类别来"解码"刺激。在初级体感皮层（S-I area 3b），每个皮层神经元整合来自 300-400 条传入纤维的输入，构建更大的感受野和更抽象的特征检测器（DiCarlo et al. 1998; Sripati et al. 2006）。例如，area 2 中某些神经元专门响应**特定方向**的运动（方向选择性），而另一些则对**边缘朝向**敏感（朝向选择性）（Warren et al. 1986; Gardner 1988）。

**（c）神经信号的信息论测度**

刺激与神经响应之间的互信息（mutual information）提供了编码效率的上界：

$$I(\mathcal{S}; \mathcal{R}) = \sum_{s \in \mathcal{S}} \sum_{r \in \mathcal{R}} P(s, r) \log_2 \frac{P(s, r)}{P(s) P(r)} \quad \text{[bits per stimulus]}$$

Johansson & Flanagan (2009) 估计，单个 FA-I 传入纤维在抓握任务中能传输约 0.1-0.3 bits/spike 的滑移相关信息，而群体编码可将互信息提升到 3-5 bits/事件——足以区分数十种不同的抓握状态。

### 2.4 指尖皮肤的生物力学模型

指尖皮肤的力学响应决定了什么信息被传递给机械感受器。图 2.10（Jones & Lederman 2006, 第22页）描绘了指尖的力-位移关系：在 0-1 N 范围内，皮肤位移量快速增加（指尖软组织变形主导），接触面积迅速扩大；对于 1 N 的力，指尖接触面积已经达到 5 N 时接触面积的约 70%。而在 >1 N 后，组织刚度急剧增加（从约 3.5 N/mm 增加到 20 N/mm），保护远端指骨免受冲击。

Serina et al. (1998) 提出的**椭球膜模型**将指尖建模为充有不可压缩流体的椭球形膜，膜被赋予皮肤的材料属性。该模型成功预测了指尖的力-位移非线性关系和力-接触面积关系。其核心方程为：

$$F(\delta) = \frac{E h}{R^2} \cdot \delta^2 + \frac{p_0}{R} \cdot \delta \quad \text{(对于小位移)}$$

其中 $$\delta$$ 为压入深度，$$E$$ 为膜的有效杨氏模量，$$h$$ 为膜厚，$$R$$ 为指尖曲率半径，$$p_0$$ 为内压。

对于切向载荷（与抓握中的滑移检测直接相关），指尖的剪切刚度和粘度可用 **Kelvin 模型**描述（Nakazawa et al. 2000）：

$$\tau(t) = G \cdot \gamma(t) + \eta \cdot \frac{d\gamma(t)}{dt}$$

其中 $$\tau$$ 为剪切应力，$$\gamma$$ 为剪切应变，$$G$$ 为剪切刚度（约 10-50 kPa 取决于压力），$$\eta$$ 为剪切粘度。

在抓握轻物体时（剪切力 < 1.2 N），指尖的剪切刚度较大，这意味着由滑移引起的指尖变形较小——这对 FA-I 滑移检测器的信噪比提出了挑战。而当抓握力增大时，指尖压缩程度增加，剪切刚度进一步变化，形成一个高度非线性的力学环境。

### 2.5 本体感觉：位置、运动和力的感知

**本体感觉**（proprioception，Sherrington 1906 年创造该术语）是指对身体各部分相对位置和运动的感觉，其感受器分布在大幅超出皮肤范围的组织中。

**（a）肌肉纺锤体（Muscle Spindle）**

肌肉纺锤体是本体感觉的核心传感器——这是位于肌肉内部的梭形结构（长 4-10 mm），内含 2-14 条梭内纤维，"并联"于梭外肌纤维（即力产生的肌纤维）。每个纺锤体通常有 1-2 个初级感觉末梢（由 Ia 类传入纤维支配）和 1-5 个次级感觉末梢（由 II 类传入纤维支配）。

初级末梢（Ia）对肌肉长度变化**和**变化速度都敏感——它们是速度和方向的检测器。次级末梢（II）的位置灵敏度比初级末梢高约 50%（Edin & Vallbo 1990），但对速度不敏感——它们是**静态位置**的检测器。

纺锤体信号并非"纯净"的长度读数，因为它们有自己的运动神经支配——**γ 运动系统**（fusimotor system）。α-γ 共激活（alpha-gamma coactivation）意味着当运动皮层命令肌肉收缩时，纺锤体的灵敏度被同步调节，以防止在肌肉缩短时纺锤体失去张力。

**（b）高尔基腱器（Golgi Tendon Organ）**

高尔基腱器位于肌肉-肌腱连接处，与 10-20 条梭外肌纤维"串联"。它对肌纤维主动产生的力敏感，对被动拉伸不敏感。这使其成为肌肉张力的专用传感器。值得注意的是，手部蚓状肌（lumbrical muscles）似乎没有腱器（Devanandan et al. 1983）。

**（c）关节感受器与皮肤拉伸**

传统观点认为关节囊感受器参与本体感觉，但证据表明其作用有限。接受人工掌指关节置换手术的风湿性关节炎患者，虽然失去了关节囊感受器，但仍保持正常的关节位置感（Cross & McCloskey 1973）。手背皮肤中的 SA-II 感受器在手指运动中因皮肤拉伸而被可靠激活——单个 SA-II 传入纤维的放电率可在运动引起的皮肤拉伸下高达 50 imp/s，足以编码关节运动信息（Edin & Abbs 1991）。这说明本体感觉是一个**多源融合**的信号：来自肌肉纺锤体、高尔基腱器、关节囊感受器和皮肤拉伸感受器的信号在大脑中被整合。

**（d）随附放电与预测编码**

最后，中枢神经系统不是被动接收本体感觉信号——它使用**随附放电**（corollary discharge）或**传出副本**（efference copy）来预测自身动作的体感后果。在主动触摸中，运动皮层的传出副本被发送到体感区域，允许大脑比较"我预期感受到什么"和"我实际感受到什么"。这就是为什么你无法胳肢自己的神经机制——自生触觉被预测信号精确地减弱了（Flanagan et al. 2003; Johansson & Flanagan 2009）。

### 2.6 工程触觉传感器模型

**（a）GelSight 传感器原理**

GelSight（最初由 MIT 的 Adelson 团队开发）的核心设计包括三个组件：（1）一块透明的弹性凝胶（通常为硅胶），其一面与物体接触；（2）四色（红、绿、蓝、白）LED 阵列从不同角度照明接触表面；（3）内置摄像头通过凝胶观察接触面的形变。当一个物体压在凝胶表面时，表面的三维形变调制了每个点的光照颜色——通过光度立体（photometric stereo）原理，彩色图像可以被转换为高分辨率的深度图。

GelSight 的物理模型从**光度立体视觉**出发。设相机像素 $$(u, v)$$ 处的辐照强度为 $$I_c(u, v)$$，其中 $$c \in \{R, G, B\}$$ 表示颜色通道。LED 颜色编码使每个方向的光照具有不同的光谱成分，从而：

$$I_c(u, v) = \rho(u, v) \sum_{k \in \text{LEDs}} L_k(u, v) \cdot \max(0, \mathbf{n}(u, v) \cdot \mathbf{l}_k) \cdot S_c(k)$$

其中 $$\rho$$ 为凝胶表面反射率，$$\mathbf{n}$$ 为表面法向量，$$\mathbf{l}_k$$ 为第 $$k$$ 个 LED 的方向向量，$$S_c(k)$$ 为第 $$k$$ 个 LED 在颜色通道 $$c$$ 中的发射光谱。通过求解这一逆问题（在已知 LED 方向和光谱的标定条件下），可以从三个颜色通道恢复表面法向图，并通过积分重建深度图。GelSight 可实现约 10-30 μm 的空间分辨率和约 2 μm 的深度精度，足以分辨人类指尖可感知的最小表面特征。

**（b）BioTac 传感器原理**

BioTac（由 SynTouch 开发）采用了仿生学设计，其结构模拟了人类指尖：一个刚性的"指骨"核心上覆盖着弹性"皮肤"，两者之间填充了不可压缩的导电液体。核心上嵌入了 19 个电极，当皮肤因接触而形变时，液体流动改变电极间的阻抗。此外，核心中还有一个水听器（hydrophone）用于检测振动。

BioTac 的物理模型基于**电阻抗层析成像**（Electrical Impedance Tomography, EIT）。皮肤形变改变了核心周围流体的厚度分布，而流体厚度与电极对之间的复阻抗 $$Z_{ij}$$ 之间有如下关系：

$$Z_{ij} \propto \frac{1}{\sigma} \int_{\Gamma_{ij}} \frac{d\ell}{h(\mathbf{x})}$$

其中 $$\sigma$$ 为流体电导率，$$h(\mathbf{x})$$ 为位置 $$\mathbf{x}$$ 处的流体厚度，$$\Gamma_{ij}$$ 为电极 $$i$$ 和 $$j$$ 之间的电流路径。通过测量全部 19 个电极之间的阻抗矩阵，可以重建接触力的分布。BioTac 还能测量绝对温度（通过嵌入的热敏电阻）和振动（通过水听器，可在高达 800 Hz 的频率范围内检测微振动，模仿 FA-II/Pacinian 通道），从而以单一传感器同时实现 SA-I（压力）、FA-I（滑移）和 FA-II（振动）三种机械感受功能的工程仿生。

---

## 三、论文精读

### 3.1 Gardner (2010): Touch

Gardner 的《Touch》发表于 *Encyclopedia of Life Sciences*，是一篇全面覆盖触觉神经生物学的综述，其主要贡献在于将机械感受器的分子/细胞机制、解剖通路、皮层处理和临床病变效应整合为一条完整的因果链。

**核心论点**：触觉不是单一的感觉模态，而是由**四类机械感受器-传入纤维**并行编码的信息流在大脑中经过**串行和并行处理**形成的统一知觉。外周感受器将触觉刺激分解为独立的特征通道（压力、运动、振动、拉伸），而初级体感皮层的四个细胞构筑区（Brodmann 3a, 3b, 1, 2）以层级方式逐步整合这些特征，最终在顶叶后部构建出物体的多模态表征（形状、大小、纹理、重量、位置）。

**第1部分：外周机械转导**。Gardner 详细描述了触觉转导的生物物理原理：皮肤压痕或拉伸通过直接门控阳离子通道诱发感受器电位，其幅度与压力大小和施加速度成正比。这一核心机制将机械能转化为电信号。Pacinian 小体的层合结构作为一个机械高通滤波器被重点分析——其物理机制（不同频率下层合膜的运动耦合/解耦）优雅地解释了为什么 FA-II 传入纤维对 200 Hz 振动具有约 1 μm 的惊人灵敏度，但对 <40 Hz 的缓变刺激几乎不响应。

**第2部分：体感通路的解剖组织**。沿着两个正交轴组织：一是**拓扑轴**（somatotopy，身体的邻域关系在大脑中被保留），二是**功能轴**（不同感觉亚模态在独立的解剖束和细胞列中保持分离，直至最高层的多模态整合）。背柱-内侧丘系通路将触觉和本体感觉信息从脊髓传输到丘脑的腹后外侧核，再投射到 S-I 皮层。

**第3部分：皮层处理层级**。S-I 区域 3b 和 1 接收皮肤信息，3a 和 2 接收本体感觉信息。在每个层级，感受野逐渐变大（从指尖小点 → 整个指腹 → 多个相邻手指），神经元从"忠实地复现外周活动"过渡到"检测抽象特征"。Area 2 的某些神经元同时接收触觉和本体感觉输入——这是"手的感觉信息和运动信息第一次在单细胞水平上汇聚"。

**第4部分：体感皮层的可塑性**。也许是最引人注目的发现：截肢者常报告在缺失的肢体上产生幻触觉——因为原本代表那个肢体的皮层区域并未死亡，而是被邻近的受体输入（如面部皮肤或残肢）重新支配。同时，职业小提琴家的功能性磁共振成像（fMRI）显示其按弦手的手指皮层表征区显著大于对照组。

### 3.2 Jones & Lederman (2006): Human Hand Function

Jones & Lederman 的《Human Hand Function》为本章提供了两个关键的理论支柱：**指尖力学模型**和**主动触觉的理论框架**。

**手指运动-感觉连续体**。该书的开创性贡献之一是将手部功能概念化为一条感觉-运动连续体，包含四个类别：（1）**触觉感知**（Tactile Sensing）——手静止、被动接受触觉刺激；（2）**主动触觉感知**（Active Haptic Sensing）——手自主运动，通过探索物体来提取信息；（3）**抓握**（Prehension）——伸手抓取物体；（4）**非抓握性技能运动**——如弹钢琴或打字。这一框架比传统的"被动触觉 = 无运动，主动触觉 = 有运动"的二分法更精确地描述了手的感知和运动功能如何连续过渡。

**第3章：手部功能的神经生理学**（本书引用的核心材料）。Jones & Lederman 提供了对四类机械感受器-传入纤维单元的详细量化描述，包括（1）手指各区域的密度分布图（Vallbo & Johansson 1984 的经典数据，本书第28页），（2）每类传入纤维的适应性特征和频率响应（第27页），（3）指尖的力-位移关系和力-接触面积关系（第22页）。尤其有价值的是他们收录了 Johansson et al. (1982) 的生理学频率响应数据与 Bolanowski et al. (1988) 的心理物理学振动阈值函数之间的对比（第30页）——这张图直观地展示了四个触觉"通道"的存在，每个通道由一类机械感受器介导，并且生理学响应和心理物理学的阈值函数曲线形状高度对应。此图是我们后续触觉信号处理代码实现中**通道建模**的核心依据。

**肌肉感受器与本体感觉**。第3章还全面覆盖了肌肉纺锤体（Ia 和 II 类传入纤维）、高尔基腱器（Ib 类传入纤维）和关节感受器的结构和功能。关键发现包括：（1）手部肌肉纺锤体在安静休息位时发放率极低，甚至在肌纤维松弛时静默，这与猫的纺锤体行为完全不同，提示人类的 γ 运动系统在自然体位下对纺锤体灵敏度的调节可能比猫更弱（Vallbo 1974）；（2）当关节被主动保持在某一位置时，手指伸肌的初级和次级纺锤体传入纤维的发放率无显著差异（图 3.9，Hulliger et al. 1982），表明 γ 运动活动的变化恰好补偿了肌肉长度的变化，使得纺锤体信号本身不携带静态位置信息——位置感知必须来自中枢对传出副本和纺锤体信号的联合解码。

**指尖模型**。Serina et al. (1998) 对指尖的力-位移和力-接触面积关系的定量测量数据表明，指尖组织在低力（<1 N）区间充当高效的触觉传感器——微小的力变化导致显著的接触面积增加，最大化机械感受器的激活范围。这解释了为什么人类在精细触觉辨别任务中倾向于使用刚好足够的接触力（通常远小于 1 N）。

---

## 四、算法与代码实现

以下代码实现一套**触觉信号处理管线**：从模拟四类机械感受器的原始脉冲响应，到尖峰序列的分析、群体编码、以及触觉刺激的分类。代码使用 Python 3 + NumPy + SciPy + Matplotlib，所有核心函数均标注参考文献来源。

### 4.1 机械感受器模拟 (tactile_model.py)

```python
#!/usr/bin/env python3
"""
tactile_model.py — 四类机械感受器的计算模型
参考文献:
  - Gardner (2010) Touch. Encyclopedia of Life Sciences.
  - Johansson & Vallbo (1983) Trends in Neuroscience 6:27-32.
  - Bolanowski et al. (1988) J. Acoust. Soc. Am. 84(5):1680-1694.
  - Jones & Lederman (2006) Human Hand Function, Ch.3.
  - Johnson (2001) Current Opinion in Neurobiology 11:455-461.
"""

import numpy as np
from scipy.signal import lfilter, butter
from dataclasses import dataclass
from typing import Tuple, List, Optional

# ============================================================
# 2.1 力学前滤波: 皮肤 + 感受器囊膜的频率响应
# 资料来源: Gardner 2010 §Biophysics of sensory transduction;
#           Bell et al. (1994) 对 Pacinian 小体的膜层模型
# ============================================================

@dataclass
class MechanoreceptorParams:
    """四类机械感受器的模型参数 [UNSPECIFIED: 参数取文献中的典型值]
    参考: Johansson et al. (1982); Bolanowski et al. (1988) Fig.5"""
    name: str
    f_low: float       # 低频截止 (Hz)
    f_high: float      # 高频截止 (Hz)
    tau_adapt: float   # 适应时间常数 (s) — 区分 SA/RA
    gain: float        # 转导增益 (imp/s/mN)
    threshold: float   # 力阈值 (mN)
    spont_rate: float  # 自发发放率 (imp/s)
    rf_diameter: float # 感受野直径 (mm)

# 四类标准参数 (来源: Johnson 2001 Table 1; Bolanowski et al. 1988 Fig.5)
# 某些参数为 [UNSPECIFIED] — 论文未给出精确值, 取文献中合理范围的中值
MECHANORECEPTOR_PARAMS = {
    "SA-I": MechanoreceptorParams(
        name="SA-I (Merkel)", f_low=0.3, f_high=10.0,
        tau_adapt=5.0,    # 慢适应: 秒级时间常数
        gain=15.0,        # [UNSPECIFIED] 精确增益, 取 Johnson 2001 的中间估计
        threshold=1.5, spont_rate=0.5, rf_diameter=2.5
    ),
    "SA-II": MechanoreceptorParams(
        name="SA-II (Ruffini)", f_low=0.3, f_high=5.0,
        tau_adapt=3.0,
        gain=4.0,         # SA-II 敏感度约为 SA-I 的 1/6 (Jones & Lederman 2006 p.31)
        threshold=5.0, spont_rate=8.0, rf_diameter=10.0
    ),
    "FA-I": MechanoreceptorParams(
        name="FA-I (Meissner)", f_low=3.0, f_high=60.0,
        tau_adapt=0.025,  # 快适应: 25 ms 时间常数 [UNSPECIFIED: 精确值]
        gain=12.0,
        threshold=2.0, spont_rate=0.2, rf_diameter=3.5
    ),
    "FA-II": MechanoreceptorParams(
        name="FA-II (Pacinian)", f_low=40.0, f_high=500.0,
        tau_adapt=0.005,  # 极快适应: 5 ms [UNSPECIFIED: 精确值]
        gain=50.0,        # Pacinian 最敏感, 增益最大
        threshold=0.1, spont_rate=0.0, rf_diameter=25.0
    ),
}


def mechanical_filter(
    force: np.ndarray,
    dt: float,
    params: MechanoreceptorParams
) -> np.ndarray:
    """力学前滤波: 模拟皮肤-囊膜的频率选择性
    资料来源: Gardner 2010 §Pacinian corpuscle is designed to detect vibration

    Pacinian 小体的多层膜结构充当机械高通滤波器. 此函数用 Butterworth
    带通滤波器近似所有四类的力学频率响应.

    Args:
        force: 施加的力信号 (mN), shape (n_samples,)
        dt: 采样间隔 (s)
        params: 感受器参数

    Returns:
        滤波后的应变等效信号 (无量纲)
    """
    nyquist = 0.5 / dt
    low = params.f_low / nyquist
    high = params.f_high / nyquist

    # [UNSPECIFIED] 滤波器阶数: 实际皮肤是连续介质, 阶数更高
    # 但二阶 Butterworth 已能复现频率响应的基本形状
    b, a = butter(2, [low, high], btype='band')
    filtered = lfilter(b, a, force)
    filtered = np.maximum(filtered, 0.0)  # 半波整流: 模拟通道的整流特性
    return filtered


# ============================================================
# 2.2 转导 + 适应: 感受器电位到发放率的转换
# 资料来源: Gardner 2010 §Biophysics of sensory transduction
# ============================================================

def receptor_potential(
    strain: np.ndarray,
    dt: float,
    params: MechanoreceptorParams
) -> np.ndarray:
    r"""感受器电位模型: 一阶适应动力学
    资料来源: Gardner 2010 §Biophysics of sensory transduction;
              Eq. tau_adapt * dV/dt + V = kappa * F(t)

    微分方程离散化 (Euler 方法):

        V[k] = (tau * V[k-1] + dt * gain * strain[k]) / (tau + dt)

    RA 类型的时间常数 tau_adapt < 50ms => 快速衰减 (适应)
    SA 类型的时间常数 tau_adapt > 1s  => 缓慢衰减 (持续响应)
    """
    tau = params.tau_adapt
    n = len(strain)
    V = np.zeros(n)

    for k in range(1, n):
        # Euler 前向离散化: tau * (V[k] - V[k-1])/dt + V[k-1] = gain * strain[k]
        V[k] = (tau * V[k-1] + dt * params.gain * strain[k]) / (tau + dt)

    V = np.maximum(V, 0.0)  # 整流: 感受器电位不能为负
    return V


# ============================================================
# 2.3 尖峰编码: Leaky Integrate-and-Fire (LIF) 神经元模型
# 资料来源: 标准计算神经科学模型,
#           与 Johansson & Flanagan 2009 的发放模式一致
# ============================================================

def lif_spike_train(
    receptor_pot: np.ndarray,
    dt: float,
    tau_m: float = 0.010,     # 膜时间常数 10 ms [UNSPECIFIED]
    V_th: float = 1.0,         # 发放阈值 (归一化)
    V_reset: float = 0.0,      # 重置电位
    tau_ref: float = 0.002,    # 不应期 2 ms
    spont_rate: float = 0.0,   # 自发发放率 (imp/s)
) -> Tuple[np.ndarray, np.ndarray]:
    """Leaky Integrate-and-Fire 尖峰生成
    资料来源: 标准计算神经科学 LIF 模型;
              发放性质与 Gardner 2010 Fig.3 中的响应一致

    Returns:
        spike_times: 尖峰时间数组 (s)
        spike_train: 二进制尖峰序列 (0/1), shape (n_samples,)
    """
    n = len(receptor_pot)
    V_m = np.zeros(n)
    spike_train = np.zeros(n, dtype=np.int8)

    ref_counter = 0
    t_since_last_spike = 0.0

    # 自发发放: Poisson 过程
    if spont_rate > 0:
        # [UNSPECIFIED] 精确的随机发放模型 — 使用 Poisson 近似
        spont_prob = spont_rate * dt
        spont_spikes = np.random.rand(n) < spont_prob
    else:
        spont_spikes = np.zeros(n, dtype=bool)

    for k in range(n):
        if ref_counter > 0:
            ref_counter -= 1
            V_m[k] = V_reset
            continue

        # LIF 动力学: tau_m * dV/dt = -(V - V_rest) + I_receptor
        # 离散化 (Euler): V[k] = V[k-1] + (dt/tau_m) * (-V[k-1] + I[k])
        I_syn = receptor_pot[k]
        V_m[k] = V_m[max(k-1, 0)] + (dt / tau_m) * (-V_m[max(k-1, 0)] + I_syn)

        if spont_spikes[k]:
            V_m[k] = V_th  # 强制发放

        if V_m[k] >= V_th:
            spike_train[k] = 1
            V_m[k] = V_reset
            ref_counter = int(tau_ref / dt)

    spike_times = np.arange(n) * dt
    spike_times = spike_times[spike_train == 1]

    return spike_times, spike_train


# ============================================================
# 2.4 完整的机械感受器管线
# ============================================================

def simulate_mechanoreceptor(
    force: np.ndarray,
    dt: float,
    receptor_type: str,
    params: Optional[MechanoreceptorParams] = None,
) -> Tuple[np.ndarray, np.ndarray, np.ndarray]:
    """完整的机械感受器仿真管线

    变换链: Force → Mech Filter → Transduction → Spike Encoding

    Args:
        force: 力信号 (mN), shape (n_samples,)
        dt: 采样间隔 (s)
        receptor_type: 感受器类型 ("SA-I", "SA-II", "FA-I", "FA-II")
        params: 感受器参数 (若为 None, 使用默认参数)

    Returns:
        strain: 滤波后的应变信号
        V: 感受器电位
        spike_train: 二进制尖峰序列
    """
    if params is None:
        params = MECHANORECEPTOR_PARAMS[receptor_type]
    else:
        params = params

    # Step 1: 力学前滤波
    strain = mechanical_filter(force, dt, params)

    # Step 2: 转导 + 适应动力学
    V = receptor_potential(strain, dt, params)

    # Step 3: 尖峰编码
    _, spike_train = lif_spike_train(
        V, dt, spont_rate=params.spont_rate
    )

    return strain, V, spike_train


# ============================================================
# 示例: 对比四类感受器对同一刺激的响应
# ============================================================

if __name__ == "__main__":
    import matplotlib.pyplot as plt

    # 刺激: 持续 1s, dt = 0.1ms
    dt = 0.0001
    t = np.arange(0, 1.0, dt)
    n = len(t)

    # 复合力信号: 0.2s 开始压入, 0.5s 释放, 叠加 200Hz 微振动
    force = np.zeros(n)
    force[int(0.2/dt):int(0.5/dt)] = 10.0  # 10 mN 持续压力
    # 叠加 200 Hz, 振幅 0.5 mN 的振动 (模拟纹理扫描)
    vibration = 0.5 * np.sin(2 * np.pi * 200 * t)
    force += vibration

    fig, axes = plt.subplots(5, 1, figsize=(12, 10), sharex=True)
    axes[0].plot(t, force)
    axes[0].set_ylabel("Force (mN)")
    axes[0].set_title("Composite Stimulus (10 mN pressure + 200 Hz/0.5 mN vibration)")

    for idx, rtype in enumerate(["SA-I", "SA-II", "FA-I", "FA-II"]):
        strain, V, spike_train = simulate_mechanoreceptor(force, dt, rtype)
        ax = axes[idx + 1]
        # 显示尖峰光栅图: 每个尖峰画一个点
        spike_idx = np.where(spike_train > 0)[0]
        spike_t = t[spike_idx]
        # 降采样应变信号以加快渲染
        downsample = 100
        ax.plot(t[::downsample], strain[::downsample],
                alpha=0.5, label=f'{rtype} strain')
        if len(spike_t) > 0:
            ax.plot(spike_t, np.ones_like(spike_t) * np.max(strain) * 0.9,
                    'r|', markersize=10, label=f'{rtype} spikes')
        ax.set_ylabel(rtype)
        ax.legend(loc='upper right', fontsize=7)

    axes[-1].set_xlabel("Time (s)")
    plt.tight_layout()
    plt.savefig("mechanoreceptor_response.png", dpi=150)
    print("Figure saved to mechanoreceptor_response.png")
    print(f"SA-I  spikes: {np.sum(spike_train_for_SAI := simulate_mechanoreceptor(force, dt, 'SA-I')[2])}")
    # Note: 上面的 inline 赋值仅用于演示; 实际运行时会为每种类型独立计算
```

### 4.2 触觉信号处理与特征提取 (tactile_processing.py)

```python
#!/usr/bin/env python3
"""
tactile_processing.py — 触觉信号处理与群体编码分析
参考文献:
  - Connor et al. (1990) J. Neurosci. 10:3823-3836 (粗糙度神经编码)
  - Johansson & Flanagan (2009) Nat. Rev. Neurosci. 10:345-359
  - Sripati et al. (2006) J. Neurosci. 26:2101-2114 (时空感受野)
  - Weber et al. (2013) "GelSight: High-Resolution Robot Tactile Sensors"
"""

import numpy as np
from scipy import signal
from scipy.ndimage import gaussian_filter
from typing import Tuple, Optional


# ============================================================
# 3.1 皮层感受野模型: 空间-时间感受野 (STRF)
# 资料来源: Sripati et al. 2006, Fig.7 — 3b 区神经元的 STRFs
# ============================================================

def spatial_temporal_rf(
    grid_size: int = 32,
    sigma_center: float = 3.0,
    sigma_surround: float = 8.0,
    tau_exc: float = 0.015,
    tau_inh: float = 0.050,
    dt: float = 0.001,
    kernel_duration: float = 0.2
) -> Tuple[np.ndarray, np.ndarray]:
    r"""生成模拟的皮层神经元空间-时间感受野 (STRF)

    STRF 建模为 DoG (Difference of Gaussians) 空间核 ×
    双指数时间核 (兴奋性+抑制性):

        RF_space(x,y) = A_c * exp(-r^2/2σ_c^2) - A_s * exp(-r^2/2σ_s^2)
        RF_time(t)   = exp(-t/τ_exc) - k * exp(-t/τ_inh)

    资料来源: Sripati et al. 2006, 该研究从外周传入纤维和 3b 区
              神经元的联合记录中估计了 STRFs; DoG 模型复现了
              中心-周围拮抗结构.

    Returns:
        spatial_rf: (grid_size, grid_size) 空间核
        temporal_rf: (n_time_bins,) 时间核
    """
    # 空间 DoG 核
    x = np.linspace(-grid_size/2, grid_size/2, grid_size)
    y = np.linspace(-grid_size/2, grid_size/2, grid_size)
    xx, yy = np.meshgrid(x, y)
    rr = np.sqrt(xx**2 + yy**2)
    A_c = 1.0 / (2 * np.pi * sigma_center**2)
    A_s = 0.4 / (2 * np.pi * sigma_surround**2)  # 周围强度为中心强度的 40%
    spatial_rf = (A_c * np.exp(-rr**2 / (2 * sigma_center**2)) -
                  A_s * np.exp(-rr**2 / (2 * sigma_surround**2)))
    spatial_rf /= np.sum(np.abs(spatial_rf))  # 归一化

    # 时间双指数核
    n_t = int(kernel_duration / dt)
    t_kernel = np.arange(n_t) * dt
    k_inh = 0.6  # 抑制-兴奋比 [UNSPECIFIED: 精确值随神经元而异]
    temporal_rf = (np.exp(-t_kernel / tau_exc) -
                   k_inh * np.exp(-t_kernel / tau_inh))
    # 因果性: t < 0 时核为零 (已由 np.arange 确保)

    return spatial_rf, temporal_rf


# ============================================================
# 3.2 群体编码: 从触觉图像到发放率向量
# 资料来源: Johansson & Flanagan (2009) 的群体编码框架;
#           Connor et al. (1990) 的粗糙度编码模型
# ============================================================

def population_rate_code(
    tactile_image: np.ndarray,
    rf_params: dict,
    receptor_density: float = 1.0,     # 感受器密度 (units/mm²)
    image_scale: float = 0.1,           # 像素到 mm 的缩放因子
) -> np.ndarray:
    r"""触觉图像的群体发放率编码

    模拟 N 个 SA-I 感受器阵列对触觉图像的响应.
    每个感受器的响应由图像 patch 与其空间感受野的内积确定,
    并经过幂函数压缩 (Johnson 2001).

    粗糙度的群体编码定义为发放率在空间上的方差 (Connor et al. 1990):

        R_perceived ∝ Var_x[r_SA-I(x)]

    Args:
        tactile_image: (H, W) 触觉图像 (接触力分布, mN)
        rf_params: dict with 'sigma' 键用于控制感受野大小
        receptor_density: 每 mm² 的感受器数量 (指尖≈100)
        image_scale: 每个像素对应的 mm 数

    Returns:
        rate_array: (n_receptors_h, n_receptors_w) 发放率阵列 (imp/s)
    """
    H, W = tactile_image.shape
    sigma_pix = rf_params.get('sigma', 1.5) / image_scale

    # 对触觉图像进行高斯平滑以模拟感受野的有限空间分辨率
    # 平滑核的 sigma = 感受野半径 (≈1.25 mm for SA-I) / image_scale
    smoothed = gaussian_filter(tactile_image, sigma=sigma_pix)

    # 将平滑后的图像下采样为感受器网格
    # 感受器间距 = 1/sqrt(density) mm, 对 SA-I 约 0.1 mm
    spacing_mm = 1.0 / np.sqrt(receptor_density)
    spacing_pix = max(1, int(spacing_mm / image_scale))

    # 子采样: 每个感受器取感受野中心位置的值
    rate_array = smoothed[::spacing_pix, ::spacing_pix]

    # 幂函数压缩: SA-I 的亚线性响应 (n ≈ 0.5-0.7, Johnson 2001)
    # r(F) = r_0 + a * F^n  for F > F_th, else r_0
    F_th = 1.0  # mN 阈值
    n_power = 0.6
    a_gain = 15.0

    rate_array = np.where(
        rate_array > F_th,
        a_gain * np.power(rate_array - F_th, n_power),
        0.0
    )

    return rate_array


def roughness_index(rate_array: np.ndarray) -> float:
    r"""粗糙度指数: SA-I 群体发放率的空间变异度 (Connor et al. 1990)

    $$R_{perceived} \propto \frac{1}{N} \sum_{i=1}^{N} (r_i - \bar{r})^2$$

    """
    return np.var(rate_array)


# ============================================================
# 3.3 GelSight 仿真: 从深度图到触觉图像
# 资料来源: Weber et al. (2013) 对 GelSight 原理的描述;
#           [UNSPECIFIED] 光度立体逆问题的精确参数依赖于具体标定
# ============================================================

def simulate_gelsight(
    height_map: np.ndarray,
    led_directions: Optional[np.ndarray] = None,
    pixel_scale: float = 0.01  # mm/pixel
) -> np.ndarray:
    r"""模拟 GelSight 传感器的 RGB 输出

    从高度图 (height_map) 计算表面法向, 然后使用简化的
    光度立体模型生成模拟的 GelSight RGB 图像.

    GelSight 的核心: 每个像素的法向量编码了接触面的三维几何.

    $$I_c(u,v) = \rho \cdot \sum_k L_k \cdot \max(0, \mathbf{n} \cdot \mathbf{l}_k) \cdot S_c(k)$$

    Args:
        height_map: (H, W) 高度图 (mm), 正值表示凸起
        led_directions: (N, 3) LED 方向向量数组
        pixel_scale: 每个像素对应的 mm 数 (用于梯度计算)

    Returns:
        gelsight_rgb: (H, W, 3) 模拟的 GelSight RGB 图像 (0-255)
    """
    H, W = height_map.shape

    # 通过中心差分计算梯度
    gy, gx = np.gradient(height_map, pixel_scale)

    # 法向量: n = (-gx, -gy, 1) / ||(-gx, -gy, 1)||
    nz = np.ones_like(height_map)
    norm = np.sqrt(gx**2 + gy**2 + 1.0)
    nx = -gx / norm
    ny = -gy / norm
    nz = nz / norm

    # 默认 LED 配置: 4 个 LED, 方向分别为上、下、左、右, 各倾斜 30°
    # [UNSPECIFIED] 精确的 LED 倾角和色彩分配 — 基于 Weber et al. 2013 的典型配置
    if led_directions is None:
        theta = np.radians(30.0)  # 倾斜角
        led_directions = np.array([
            [np.sin(theta), 0, np.cos(theta)],          # 红色: 右
            [-np.sin(theta), 0, np.cos(theta)],         # 绿色: 左
            [0, np.sin(theta), np.cos(theta)],          # 蓝色: 上
            [0, -np.sin(theta), np.cos(theta)],         # 白色: 下 (混合为蓝色通道)
        ])

    # 颜色分配
    color_weights = np.array([
        [1.0, 0.0, 0.0],  # 红色
        [0.0, 1.0, 0.0],  # 绿色
        [0.0, 0.0, 1.0],  # 蓝色
        [0.3, 0.3, 0.4],  # 白色: 混合 RGB
    ])

    rgb = np.zeros((H, W, 3))
    for k in range(len(led_directions)):
        lk = led_directions[k] / np.linalg.norm(led_directions[k])
        # Lambertian 反射: I = max(0, n·l)
        cos_theta = np.maximum(0.0,
                               nx * lk[0] + ny * lk[1] + nz * lk[2])
        for c in range(3):
            rgb[:, :, c] += cos_theta * color_weights[k, c]

    # 归一化到 [0, 255]
    rgb = np.clip(rgb / np.max(rgb + 1e-8) * 255.0, 0, 255).astype(np.uint8)
    return rgb


# ============================================================
# 3.4 振动触觉频率分析 (微纹理感知)
# 资料来源: Bensmaia & Hollins (2003) "The vibrations of texture"
#           Johansson & Flanagan (2009) §Tactile signals
# ============================================================

def vibratory_texture_analysis(
    acceleration_signal: np.ndarray,
    fs: float = 2000.0,
    scanning_velocity: float = 50.0,  # mm/s
) -> dict:
    r"""振动触觉信号中的纹理信息提取

    当指尖扫描纹理表面时, 皮肤加速度中的频谱峰值
    对应于表面周期特征的频率: f_peak = v_scan / λ_spatial
    (Bensmaia & Hollins 2003).

    此函数模拟从加速度信号中提取空间纹理周期.

    Args:
        acceleration_signal: 皮肤加速度信号 (m/s²)
        fs: 采样率 (Hz)
        scanning_velocity: 扫描速度 (mm/s)

    Returns:
        dict with keys:
          - 'frequencies': 频率轴 (Hz)
          - 'psd': 功率谱密度
          - 'spatial_periods': 空间周期 (mm) = v_scan / frequency
          - 'dominant_spatial_period': 主导空间周期 (mm)
    """
    n = len(acceleration_signal)
    freqs = np.fft.rfftfreq(n, d=1.0/fs)
    fft_vals = np.fft.rfft(acceleration_signal)
    psd = np.abs(fft_vals)**2 / n

    # 排除 DC 和 Nyquist
    valid = (freqs > 0) & (freqs < fs/2)
    freqs_valid = freqs[valid]
    psd_valid = psd[valid]

    # 频率 → 空间周期的转换: spatial_period = velocity / frequency
    spatial_periods = scanning_velocity / (freqs_valid + 1e-8)

    # 主导空间周期 (在合理范围内: 0.1 mm — 10 mm)
    spatial_mask = (spatial_periods >= 0.1) & (spatial_periods <= 10.0)
    if np.any(spatial_mask):
        dominant_idx = np.argmax(psd_valid[spatial_mask])
        dominant_spatial_period = spatial_periods[spatial_mask][dominant_idx]
    else:
        dominant_spatial_period = np.nan

    return {
        'frequencies': freqs_valid,
        'psd': psd_valid,
        'spatial_periods': spatial_periods,
        'dominant_spatial_period': dominant_spatial_period,
    }


# ============================================================
# 示例: 生成触觉图像并使用群体编码分析
# ============================================================

if __name__ == "__main__":
    import matplotlib.pyplot as plt

    # 生成模拟的 GelSight 高度图: 一个 Braille 点阵 (2x3 凸点)
    H, W = 128, 128
    x = np.linspace(-6, 6, W)
    y = np.linspace(-6, 6, H)
    xx, yy = np.meshgrid(x, y)
    height = np.zeros((H, W))
    # 放置 6 个 Braille 点 (间距 2.3 mm, 符合 Braille 标准)
    for dx, dy in [(-2.3, -1.15), (-2.3, 1.15), (0, -1.15),
                    (0, 1.15), (2.3, -1.15), (2.3, 1.15)]:
        dist = np.sqrt((xx - dx)**2 + (yy - dy)**2)
        height += 0.5 * np.exp(-dist**2 / (2 * 0.5**2))  # 每个点高 0.5 mm

    # 生成模拟 GelSight RGB 图像
    gelsight_img = simulate_gelsight(height)

    # 转换为力等效图像 (简化为 height * 弹性刚度)
    # [UNSPECIFIED] GelSight 凝胶的弹性模量 — 使用典型硅胶值
    gel_stiffness = 5.0  # mN/mm 每单位高度
    force_map = height * gel_stiffness

    # 群体编码
    rate_array = population_rate_code(
        force_map,
        rf_params={'sigma': 1.25},  # SA-I 感受野半径 ≈ 1.25 mm
        receptor_density=100.0,      # 指尖 SA-I 密度 ≈ 100/cm²
        image_scale=12.0/128         # 图像宽度 12 mm / 128 pixels
    )
    roughness = roughness_index(rate_array)

    # 可视化
    fig, axes = plt.subplots(1, 4, figsize=(16, 4))
    axes[0].imshow(height, cmap='gray')
    axes[0].set_title("Height Map (Braille Dots)")
    axes[0].axis('off')

    axes[1].imshow(gelsight_img)
    axes[1].set_title("Simulated GelSight RGB")
    axes[1].axis('off')

    axes[2].imshow(force_map, cmap='hot')
    axes[2].set_title("Force Map (mN)")
    axes[2].axis('off')

    axes[3].imshow(rate_array, cmap='viridis')
    axes[3].set_title(f"SA-I Population Rate\nRoughness = {roughness:.1f}")
    axes[3].axis('off')

    plt.tight_layout()
    plt.savefig("tactile_population_code.png", dpi=150)
    print("Figure saved to tactile_population_code.png")
    print(f"Roughness index: {roughness:.2f}")
```

---

## 五、核心 Takeaway 与延伸阅读

### 五个核心洞察

**1. 触觉是"四个通道"的并行编码，不是单一感觉。**
四类机械感受器（SA-I、SA-II、FA-I、FA-II）各自提取物理刺激的不同方面：SA-I 编码空间细结构，FA-I 编码皮肤运动，FA-II 编码远程振动，SA-II 编码皮肤拉伸。四个并行通道将触觉刺激分解为独立的特征流，大脑通过重新整合这些流来构建统一的触觉。这一组织原则意味着：机器人触觉系统不应追求"万能传感器"，而应设计多通道、互补的传感器阵列，每个通道专门化处理特定特征。

**2. 感受器的力学前滤波决定了它"看到"什么信息。**
触觉不是从力到神经脉冲的直接映射。皮肤的粘弹性、感受器囊膜的结构力学（Pacinian 小体的层合结构是最经典的例子）在物理层面完成了信息的预滤波——这些物理滤波器决定了什么频率、什么方向的刺激能够到达神经末梢。理解这些力学前滤波器，是设计有效触觉传感器（如 GelSight 的硅胶膜选择）的物理基础。

**3. 主动触觉是感觉-运动闭环的控制过程，不是被动的信息接收。**
触觉（被动）和主动触觉（主动运动中获取的体感信息）在功能上存在本质区别。主动触觉涉及：运动命令（传出副本）→ 预测体感后果 → 与实际体感反馈比较 → 误差用于更新运动策略。这一闭环解释了为什么人类在抓握时可以自动适应物体的重量和摩擦系数变化，也映射到机器人抓取控制的前馈+反馈结构。

**4. 本体感觉是多模态融合的结果，而非单一感受器的输出。**
没有任何单一类型的感受器（甚至包括肌肉纺锤体）能独自编码精确的关节位置或运动速度。本体感觉是肌肉纺锤体（长度+速度）、高尔基腱器（力）、关节囊感受器（极限位）、以及皮肤 SA-II（拉伸）信息的贝叶斯融合，被随附放电（对自生动作的预期）进一步调制。

**5. GelSight 和 BioTac 代表了触觉传感器设计的两种互补范式。**
GelSight 遵循"光学范式"：用高分辨率 CCD 相机通过硅胶-光照界面获得丰富的空间几何信息。BioTac 遵循"仿生范式"：用人造皮肤、导电液体、电极阵列和水听器来模拟四类机械感受功能。两者各有优势：GelSight 在形状和纹理分辨率上远优于 BioTac，但 BioTac 能同时提供力分布、振动和温度——更接近人类手指的多模态特性。下一代触觉传感器需要融合这两种范式。

### 三个开放问题

1. **如何在机器人触觉系统中实现自适应增益控制？** 人类触觉系统通过 γ 运动系统动态调节肌肉纺锤体的灵敏度，并通过下行抑制通路控制丘脑和脑干核团对触觉信号的增益。机器人目前缺乏类似的"注意力驱动的触觉灵敏度调节"机制。

2. **如何设计可泛化的触觉表征？** 人类可以在完全不同的物体（橡胶球 vs. 玻璃杯）之间迁移抓握技能——因为触觉表征抽象了材料的特定属性。当前的机器人触觉表征大多是传感器特定的低维特征向量，缺乏跨物体、跨任务的可迁移性。

3. **主动触觉的探索策略如何学习？** 生物学设定了触觉探索的最优策略（例如，判断纹理时用横向扫描，判断硬度时用按压，判断形状时用轮廓追踪），但机器人如何在强化学习框架下自主发现这些探索策略仍然是一个开放问题。

### 推荐延伸阅读

- **Mountcastle, V.B. (2005).** *The Sensory Hand: Neural Mechanisms of Somatic Sensation.* Harvard University Press. — 触觉神经科学的"圣经"级著作。
- **Johansson, R.S. & Flanagan, J.R. (2009).** "Coding and use of tactile signals from the fingertips in object manipulation tasks." *Nature Reviews Neuroscience* 10:345-359. — 灵巧操作中触觉编码的权威综述。
- **Dahiya, R.S. et al. (2010).** "Tactile sensing — from humans to humanoids." *IEEE Trans. Robotics* 26(1):1-20. — 从生物触觉到机器人触觉传感器的系统工程视角。
- **Yuan, W., Dong, S., & Adelson, E.H. (2017).** "GelSight: High-Resolution Robot Tactile Sensors for Estimating Geometry and Force." *Sensors* 17(12):2762. — GelSight 原理与应用的详细介绍。
- **Fishel, J.A. & Loeb, G.E. (2012).** "Bayesian exploration for intelligent identification of textures." *Frontiers in Neurorobotics* 6:4. — BioTac 传感器上的贝叶斯纹理识别算法。

---

## 六、练习与思考

### 基础练习

**6.1** 根据 Gardner (2010) 的综述，四类机械感受器分别检测物理刺激的哪些不同方面？如果一名患者失去了所有 FA-I 传入纤维（但保留了 SA-I、SA-II 和 FA-II），他/她在以下任务中会遇到什么困难：(a) 阅读 Braille，(b) 检测手中酒杯开始滑落，(c) 感知电脑风扇的振动？

**6.2** 推导指尖皮肤的力-位移本构关系。假设指尖为半径为 R 的半球，皮肤为厚度 h 的线弹性膜（杨氏模量 E），内部充满不可压缩流体。证明：在小位移 δ << R 条件下，$$F \propto \delta^2$$（提示：使用 Hertz 接触理论的修正版本）。结合 Serina et al. (1998) 的实验数据，讨论在 δ > 1 mm 时偏离二次关系的物理原因。

**6.3** 编程练习：使用本章提供的 `tactile_model.py` 代码，模拟四类感受器对以下刺激的响应：(a) 5 mN 的阶跃压力（模拟静态抓握），(b) 振幅 0.2 mN、频率 10 Hz 的正弦力（模拟缓慢纹理扫描），(c) 振幅 0.2 mN、频率 200 Hz 的正弦力（模拟快速振动）。记录每种感受器的尖峰数，并解释为什么不同感受器对不同刺激的响应差异如此大。

### 进阶思考

**6.4** 在 GelSight 传感器中，光度立体从三个颜色通道恢复表面法向。考虑一个简化的一维情况：表面高度为 $$h(x)$$，法向量倾角为 $$\theta(x) = \arctan(dh/dx)$$，两个 LED 分别从角度 $$+\alpha$$ 和 $$-\alpha$$ 照明。证明：在 Lambertian 反射假设下，强度比 $$R = I_+/I_-$$ 唯一地决定了 $$\theta$$。推导 $$h(x)$$ 从 $$R(x)$$ 的积分重建公式。

**6.5** Connor et al. (1990) 提出感知粗糙度与 SA-I 群体发放率的空间方差成正比。批判性地分析：为什么是方差，而不是均值？如果一个纹理的表面特征间距小于单个 SA-I 感受野的直径（即超出了空间分辨率极限），方差还能可靠编码粗糙度吗？在什么条件下 FA-I 通道可能提供补充信息？

**6.6** 设计一个完整的机器人触觉系统（包含传感器选择、信号处理管线、特征提取和控制策略），使得一个多指机器手能完成以下任务：(a) 从一堆相似的物体中（如不同尺寸的螺丝）通过触觉识别并拾取目标物体，(b) 在抓握过程中检测滑移并实时调整抓握力，(c) 通过触觉识别表面材料（金属/塑料/木材）。请列出每一步的信号流，并讨论哪些生物机械感受功能你需要用哪种工程传感器来实现。


---

<!-- chapter-nav -->
<div style="display:flex; justify-content:space-between; align-items:center; padding:1em 0;">
  <div><a href="ch05_人手与机器人手.md">← 第05章 人手与灵巧操作</a></div>
  <div><a href="index.md">↑ 目录</a></div>
  <div><a href="ch07_运动控制的发展视角.md">第07章 运动控制的发展视角 →</a></div>
</div>
