在这里插入图片描述

在足底压力采集、步态分析、康复机器人和外骨骼研究中，热力图只能回答“哪里受力较大”。如果希望继续观察载荷怎样从足跟转移到前足，通常还需要计算压力中心轨迹，即 Center of Pressure，简称 CoP。

不少系统能够直接显示 CoP 曲线，但当我们拿到原始压力矩阵后，仍然会遇到几个工程问题：

  • CoP 是不是最高压力点？
  • 压力值能否直接作为权重？
  • 不同传感单元面积不一致时怎样处理？
  • 刚着地和即将离地时，轨迹为什么容易跳动？
  • 左右脚的 CoP 能否直接画在同一张图上？

本文从数据结构开始，给出一个可以直接运行的 NumPy 实现，并说明实际项目中最容易踩到的坑。

本文代码使用模拟压力数据验证算法，不是受试者实测结果，也不用于医疗诊断。

1. CoP不是“最红点”的轨迹

压力热力图中的最红区域，通常对应局部压力较高的位置。CoP 则是所有有效受力位置共同形成的加权中心。

假设当前帧有多个压力单元，第 i 个单元的位置为 (x_i, y_i)，法向力为 F_i，那么：

x C o P = ∑ i F i x i ∑ i F i x_{CoP}=\frac{\sum_i F_i x_i}{\sum_i F_i} xCoP​=∑i​Fi​∑i​Fi​xi​​

y C o P = ∑ i F i y i ∑ i F i y_{CoP}=\frac{\sum_i F_i y_i}{\sum_i F_i} yCoP​=∑i​Fi​∑i​Fi​yi​​

因此，局部峰值会把 CoP 拉向自己，但 CoP 不一定与最高压力单元重合。

在这里插入图片描述

上图展示了支撑期内 CoP 从足跟向前足移动的基本过程。右下角的示意还说明：最高压力点与受力加权中心可以位于不同位置。

2. 输入数据应该是什么形状？

本文假设压力数据使用三维数组：

pressure.shape == (T, H, W)

其中：

  • T：时间帧数量；
  • H：压力矩阵行数；
  • W：压力矩阵列数；
  • pressure[t, r, c]：第 t 帧、第 r 行、第 c 列的压力值。

另外还需要两个坐标数组：

x_mm.shape == (W,)
y_mm.shape == (H,)

x_mm 表示每一列在鞋垫横向上的物理位置，y_mm 表示每一行在鞋垫前后方向上的物理位置，本文统一使用毫米。

如果硬件不是规则矩阵，而是约 190 个不规则分布的离散感知点，也可以把数据整理成 (T, N)，并为每个感知点提供独立的 (x_i, y_i) 和有效面积。核心加权公式并不会改变。

3. 为什么压力值要先转换成力？

严格来说，CoP 公式的权重应该是力，而不是压力。

当所有感知单元的有效面积完全一致时，面积因子在分子和分母中会相互抵消，直接以压力作为权重可以得到相同坐标。

但如果不同单元的面积不一致，就必须先计算：

F i = P i A i F_i=P_i A_i Fi​=Pi​Ai​

本文输入压力单位为 kPa、面积单位为 mm²。换算关系是：

1 kPa × 1 mm² = 0.001 N

对应代码为：

force_n = pressure_kpa * cell_area_mm2 * 1e-3

忽略不等面积问题，会让面积较大的感知单元和面积较小的单元拥有不合理的相同权重。

4. NumPy实现：逐帧计算单脚CoP

下面的函数支持：

  • (T, H, W) 压力矩阵；
  • 标量或 (H, W) 单元有效面积；
  • 自动处理 NaN、正负无穷和负压力；
  • 低载荷帧输出 NaN；
  • 同时返回 CoP 坐标和估算总法向力。
from __future__ import annotations

import numpy as np


def calculate_cop(
    pressure_kpa: np.ndarray,
    x_mm: np.ndarray,
    y_mm: np.ndarray,
    cell_area_mm2: float | np.ndarray = 1.0,
    min_total_force_n: float = 5.0,
) -> tuple[np.ndarray, np.ndarray]:
    """从足底压力矩阵计算每一帧的单脚 CoP。"""
    pressure = np.asarray(pressure_kpa, dtype=float)
    x = np.asarray(x_mm, dtype=float)
    y = np.asarray(y_mm, dtype=float)

    if pressure.ndim != 3:
        raise ValueError("pressure_kpa 必须是 (T, H, W) 三维数组")

    frames, rows, cols = pressure.shape
    if x.shape != (cols,) or y.shape != (rows,):
        raise ValueError("x_mm、y_mm 长度必须与矩阵列数、行数一致")

    area = np.asarray(cell_area_mm2, dtype=float)
    if area.ndim == 0:
        area = np.full((rows, cols), float(area))
    if area.shape != (rows, cols):
        raise ValueError("cell_area_mm2 必须是标量或 (H, W) 数组")
    if np.any(area <= 0):
        raise ValueError("单元有效面积必须大于 0")

    # 无效数值清零；标定后的微小负值不参与受力计算。
    pressure = np.nan_to_num(
        pressure,
        nan=0.0,
        posinf=0.0,
        neginf=0.0,
    )
    pressure = np.maximum(pressure, 0.0)

    # kPa × mm² × 0.001 = N
    force_n = pressure * area[None, :, :] * 1e-3
    total_force_n = force_n.sum(axis=(1, 2))
x_grid, y_grid = np.meshgrid(x, y)
    weighted_x = (force_n * x_grid[None, :, :]).sum(axis=(1, 2))
    weighted_y = (force_n * y_grid[None, :, :]).sum(axis=(1, 2))

    cop_mm = np.full((frames, 2), np.nan, dtype=float)
    valid = total_force_n >= min_total_force_n

    cop_mm[valid, 0] = weighted_x[valid] / total_force_n[valid]
    cop_mm[valid, 1] = weighted_y[valid] / total_force_n[valid]
    return cop_mm, total_force_n

返回的 cop_mm 形状是 (T, 2)：

cop_mm[:, 0] -> 横向CoP，单位mm
cop_mm[:, 1] -> 前后方向CoP，单位mm

无效帧被标记为 NaN，这样在后续统计和绘图时不会被误认为坐标原点。

5. 构造一组模拟步态压力数据

下面用移动的二维高斯分布模拟载荷从足跟向前足转移。它只用于验证代码流程，不代表真实足底压力形态。

def make_demo_pressure():
    frames, rows, cols = 80, 60, 24
    x_mm = np.linspace(-45.0, 45.0, cols)
    y_mm = np.linspace(0.0, 260.0, rows)
    x_grid, y_grid = np.meshgrid(x_mm, y_mm)

    pressure = np.zeros((frames, rows, cols), dtype=float)

    for index, phase in enumerate(np.linspace(0.0, 1.0, frames)):
        center_x = 5.0 * np.sin(phase * np.pi * 1.4) - 2.0
        center_y = 25.0 + 215.0 * phase
        amplitude = 260.0 * np.sin(np.pi * phase) ** 0.35

        main = np.exp(
            -((x_grid - center_x) ** 2 / (2 * 18.0**2))
            -((y_grid - center_y) ** 2 / (2 * 28.0**2))
        )

        secondary = 0.35 * np.exp(
            -((x_grid + 15.0) ** 2 / (2 * 15.0**2))
            -((y_grid - (center_y - 18.0)) ** 2 / (2 * 24.0**2))
        )

        pressure[index] = amplitude * (main + secondary)

    return pressure, x_mm, y_mm

调用计算函数：

pressure, x_mm, y_mm = make_demo_pressure()

cop_mm, total_force_n = calculate_cop(
    pressure,
    x_mm,
    y_mm,
    cell_area_mm2=12.0,
    min_total_force_n=5.0,
)

valid = np.isfinite(cop_mm[:, 0])

print(f"valid frames: {valid.sum()} / {len(valid)}")
print(f"peak total force: {np.nanmax(total_force_n):.2f} N")
print("first valid CoP:", np.round(cop_mm[valid][0], 2))
print("last valid CoP: ", np.round(cop_mm[valid][-1], 2))

本文示例的运行结果为：

valid frames: 78 / 80
peak total force: 708.45 N
first valid CoP (mm): [-3.73  32.90]
last valid CoP (mm):  [-8.03 225.58]

从结果可以看到，CoP 的前后方向坐标由约 33 mm 移动到约 226 mm，符合模拟数据从足跟向前足转移的设定。

6. 为什么低载荷帧必须设置阈值？

在刚着地或即将离地时，总载荷接近零。此时 CoP 公式的分母很小，一个感知点的噪声或零点漂移就可能让坐标大幅跳动。

错误做法是：

cop_x = weighted_x / total_force

然后把无穷大、异常坐标或 (0, 0) 继续连接进轨迹。

更合理的做法是根据系统噪声、受试者体重和实验任务设置有效接触阈值，只在总载荷足够大时计算 CoP：

valid = total_force_n >= min_total_force_n

阈值不能为了“让曲线更好看”而随意调整。正式实验中应记录阈值依据，并在所有受试者和条件中使用一致规则。

7. 不规则离散感知点怎样计算？

很多智能鞋垫不是完整的规则矩阵，而是若干离散感知点。此时可以把每帧数据表示为 (T, N)：

def calculate_cop_points(
    pressure_kpa,
    sensor_x_mm,
    sensor_y_mm,
    sensor_area_mm2,
    min_total_force_n=5.0,
):
    pressure = np.asarray(pressure_kpa, dtype=float)
    pressure = np.maximum(
        np.nan_to_num(pressure, nan=0.0, posinf=0.0, neginf=0.0),
        0.0,
    )

    area = np.asarray(sensor_area_mm2, dtype=float)
    force_n = pressure * area[None, :] * 1e-3
    total_force_n = force_n.sum(axis=1)

    cop = np.full((pressure.shape[0], 2), np.nan)
    valid = total_force_n >= min_total_force_n

    cop[valid, 0] = (
        force_n[valid] * np.asarray(sensor_x_mm)
    ).sum(axis=1) / total_force_n[valid]

    cop[valid, 1] = (
        force_n[valid] * np.asarray(sensor_y_mm)
    ).sum(axis=1) / total_force_n[valid]

    return cop, total_force_n

这里最关键的不是代码形式，而是坐标和有效面积必须来自真实硬件布局，不能用数组下标代替物理坐标。

8. 左右脚坐标为什么经常画反？

左右鞋垫可能使用相同的局部数组方向，也可能使用镜像坐标。

例如，左脚“向内”可能对应 x 增大，右脚“向内”却对应 x 减小。如果直接比较原始 x，会把相同的内侧偏移解释成相反方向。

建议在项目开始时明确：

x轴正方向：统一指向足的内侧，还是统一指向实验室右侧？
y轴正方向：从足跟指向足尖，还是从足尖指向足跟？
原点：足跟中心、鞋垫左下角，还是传感器阵列起点？

用于左右脚足内比较时，可以把右脚坐标镜像到统一的解剖方向；用于计算双脚整体 CoP 时，则需要知道两只脚在同一个实验室坐标系中的位置和朝向。

只有两只鞋垫各自的局部坐标，并不足以恢复人体在地面上的整体 CoP。

9. 平滑CoP轨迹时要注意什么？

原始 CoP 可能带有小幅抖动，适当滤波能够提高可读性，但滤波会改变轨迹长度、速度和转折细节。

常见错误包括：

  • 先把无效帧填成零，再进行低通滤波；
  • 不考虑采样率，直接使用固定窗口长度；
  • 为了消除首尾跳动，过度平滑整个支撑期；
  • 对不同实验条件使用不同滤波参数；
  • 只保存平滑结果，不保留原始轨迹。

工程上建议先完成有效接触分割，再对每个独立支撑期内部进行滤波。报告轨迹长度和速度时，必须说明滤波方法、截止频率或窗口大小。

10. 一条可靠的处理流水线

从智能鞋垫原始压力到 CoP 指标，可以按下面的顺序实现：

读取原始压力与时间戳
        ↓
检查丢包、重复帧和通道数量
        ↓
应用零点与标定参数
        ↓
负值截断、坏点处理、饱和检查
        ↓
压力 × 有效面积 → 单元力
        ↓
根据总载荷划分有效接触阶段
        ↓
逐帧计算CoP
        ↓
统一左右脚坐标与支撑期时间轴
        ↓
必要时进行有记录的滤波
        ↓
计算轨迹长度、前后推进、内外侧偏移和速度

不要从一张插值热力图截图反推 CoP。插值图主要用于显示，原始传感单元数据才是计算输入。

11. 常见错误清单

错误1：直接寻找每帧最大值坐标

最大值坐标是局部峰值位置，不是 CoP。

错误2：不同面积单元使用相同权重

应先由压力和有效面积计算力。

错误3：总载荷接近零时继续计算

会产生严重跳点，应使用有依据的接触阈值。

错误4：使用数组下标作为毫米坐标

只能得到像素坐标，无法与不同鞋码或设备比较。

错误5：左右脚未统一方向

内外侧结论可能被完全反转。

错误6：把单脚CoP当成双脚整体CoP

后者还需要左右脚在共同空间中的位置与载荷关系。

错误7：看到轨迹偏移就给出诊断

CoP 需要与足底分区压力、步速、动作、鞋类、IMU或动作捕捉及临床信息结合解释。

12. 星枢智步鞋垫在多模态步态研究中的使用思路

星枢智步鞋垫（BioLink-X Insole Plus / Pro）的科研方案，将约 190 个足底压力感知点、双 IMU 运动感知以及弯曲／形变感知纳入同一采集体系，并支持最高 100 Hz 同步采集与科研软件开放。

在 CoP 分析中，压力矩阵可用于计算载荷中心的空间移动；同步的 IMU 和弯曲数据则可以帮助研究者进一步判断：轨迹变化发生在哪个运动阶段、足部姿态是否同时变化、前足推进和弯曲过程是否对应。

销售方说明：星枢智步鞋垫由威海启航信息技术有限公司销售。

查看威海启航信息技术有限公司销售的星枢智步鞋垫产品信息与科研方案

正式选型或二次开发时，建议重点确认：原始通道坐标与有效面积、标定参数、时间戳、双足同步方式、数据导出格式以及 CoP 的具体计算规则。具体配置、接口和性能以正式技术参数表及实际 Demo 为准。

13. 总结

从压力矩阵计算 CoP，公式本身并不复杂：对所有有效受力位置进行加权平均即可。

真正决定结果是否可信的，是公式之外的工程细节：

  • 压力是否经过正确标定；
  • 单元有效面积是否一致；
  • 物理坐标是否准确；
  • 低载荷帧如何处理；
  • 左右脚坐标是否统一；
  • 滤波参数是否有记录；
  • 单脚 CoP 与整体 CoP 是否被正确区分。

把这些条件说明白，CoP 才不是软件界面里的一条装饰线，而是可以复核、比较和继续分析的步态指标。

完整示例代码见同目录文件 cop_numpy.py。

参考资料

  1. Hu X, Zhao J, Peng D, Qu X. Estimation of Foot Plantar Center of Pressure Trajectories with Low-Cost Instrumented Insoles Using an Individual-Specific Nonlinear Model. Sensors, 2018.
  2. Wang D, Cai P, Mao Z. The configuration of plantar pressure sensing cells for wearable measurement of COP coordinates. BioMedical Engineering OnLine, 2016.
  3. Cudejko T, Button K, Al-Amri M. Wireless pressure insoles for measuring ground reaction forces and trajectories of the centre of pressure during functional activities. Scientific Reports, 2023.
  4. Warmerdam E, Burger LM, Mergen DF, et al. The walking surface influences vertical ground reaction force and centre of pressure data obtained with pressure-sensing insoles. Frontiers in Digital Health, 2024.

推荐标签： Python、NumPy、智能鞋垫、足底压力、步态分析、生物力学、传感器

说明：本文用于算法实现、科研设备选型和实验设计交流，不构成医疗诊断或治疗建议。配图为原理示意，不代表真实受试者数据。

# 矩阵 # numpy # 线性代数
Logo
DAMO开发者矩阵

DAMO开发者矩阵，由阿里巴巴达摩院和中国互联网协会联合发起，致力于探讨最前沿的技术趋势与应用成果，搭建高质量的交流与分享平台，推动技术创新与产业应用链接，围绕“人工智能与新型计算”构建开放共享的开发者生态。

更多推荐

cover

VLAct 详细总结

avatar DAMO开发者矩阵
cover

Sensodrive｜SensoJoint 力控关节模组：连续体机器人绳张力闭环的实验平台方案

avatar DAMO开发者矩阵
cover

VLAct 数据清洗总结 (对应代码)

avatar DAMO开发者矩阵