从动态规划到DTW：一个Python可视化教程，带你亲手画出时间规整路径图

从动态规划到DTW：一个Python可视化教程，带你亲手画出时间规整路径图

在信号处理和机器学习领域，时间序列的相似性比较是一个基础但极具挑战性的问题。想象一下，当你需要比较两段语音、心电图或股票走势时，简单的逐点对比往往会得到反直觉的结果——这正是动态时间规整(DTW)算法大显身手的地方。

传统欧氏距离在处理时间序列时有个致命弱点：它要求两个序列必须严格对齐。但现实中，相似的波形往往存在时间轴上的非线性变形。比如两个人以不同语速说同一个单词，或者同一首歌曲的快慢版本。DTW通过动态规划的思想，巧妙地解决了这个问题。

本文将带你用Python的Matplotlib库，从零开始可视化DTW的核心计算过程。不同于直接调用现成库的黑箱操作，我们会一步步构建累积距离矩阵，动态展示最优路径的搜索过程，最终绘制出直观的时间规整对齐图。这种"可视化推导"的方式，能让你真正理解DTW如何实现时间轴的弹性匹配。

1. 准备工作与环境配置

首先确保你的Python环境已安装必要的科学计算库。推荐使用Anaconda发行版，它已经集成了我们所需的大部分工具：

conda install numpy matplotlib

为了更直观地展示动态过程，我们还会用到Matplotlib的动画功能。以下是基础导入语句：

import numpy as np import matplotlib.pyplot as plt from matplotlib.animation import FuncAnimation from IPython.display import HTML

让我们创建两个具有时间变形关系的示例序列。这里构造一个正弦波和一个经过非线性拉伸的相同波形：

t = np.linspace(0, 10, 100) seq1 = np.sin(t) # 原始正弦波 # 创建非线性变形的时间轴 warped_time = np.sqrt(t) * 3 seq2 = np.sin(warped_time) # 变形后的正弦波 plt.figure(figsize=(12,4)) plt.plot(t, seq1, label="原始序列") plt.plot(t, seq2, label="变形序列") plt.legend() plt.show()

这个简单的例子展示了时间规整问题的本质——两个序列形态相似，但在时间轴上存在非线性对应关系。直接计算欧氏距离会得到很大的数值，尽管人眼能明显看出它们的相似性。

2. DTW核心算法解析

DTW的核心思想是通过构建累积距离矩阵，寻找两个序列间的最优匹配路径。让我们分解这个过程的每一步。

2.1 构建局部距离矩阵

首先计算两个序列所有点对之间的局部距离。通常使用欧氏距离：

def compute_local_dist(seq1, seq2): return np.array([[abs(x - y) for y in seq2] for x in seq1]) local_dist = compute_local_dist(seq1, seq2)

这个n×m矩阵（n和m是两个序列的长度）中的每个元素代表对应点对的局部不相似度。我们可以用热图直观展示：

plt.figure(figsize=(8,6)) plt.imshow(local_dist, origin='lower', cmap='viridis') plt.colorbar(label="局部距离") plt.xlabel("序列2索引") plt.ylabel("序列1索引") plt.title("局部距离矩阵") plt.show()

2.2 累积距离矩阵的动态构建

DTW的精华在于通过动态规划逐步构建累积距离矩阵。递推公式为：

γ(i,j) = local_dist(i,j) + min(γ(i-1,j), γ(i,j-1), γ(i-1,j-1))

让我们用动画展示这个构建过程：

def init(): im.set_data(np.zeros_like(local_dist)) return [im] def update(frame): i, j = frame if i == 0 and j == 0: gamma[i,j] = local_dist[i,j] elif i == 0: gamma[i,j] = local_dist[i,j] + gamma[i,j-1] elif j == 0: gamma[i,j] = local_dist[i,j] + gamma[i-1,j] else: gamma[i,j] = local_dist[i,j] + min(gamma[i-1,j], gamma[i,j-1], gamma[i-1,j-1]) im.set_array(gamma) return [im] gamma = np.zeros_like(local_dist) fig, ax = plt.subplots(figsize=(8,6)) im = ax.imshow(gamma, origin='lower', cmap='viridis', vmax=local_dist.max()*10) plt.colorbar(im, label="累积距离") plt.xlabel("序列2索引") plt.ylabel("序列1索引") plt.title("累积距离矩阵构建过程") frames = [(i,j) for i in range(len(seq1)) for j in range(len(seq2))] ani = FuncAnimation(fig, update, frames=frames, init_func=init, blit=True, interval=50) HTML(ani.to_jshtml())

这段动画会逐步填充累积距离矩阵，你可以清晰地看到最小值路径是如何形成的。矩阵右下角的值就是两个序列的DTW距离。

3. 回溯最优路径

有了完整的累积距离矩阵后，我们需要从终点(n,m)回溯到起点(0,0)找出最优路径：

def trace_path(gamma): path = [] i, j = gamma.shape[0]-1, gamma.shape[1]-1 path.append((i, j)) while i > 0 or j > 0: if i == 0: j -= 1 elif j == 0: i -= 1 else: min_val = min(gamma[i-1,j], gamma[i,j-1], gamma[i-1,j-1]) if gamma[i-1,j-1] == min_val: i -= 1 j -= 1 elif gamma[i-1,j] == min_val: i -= 1 else: j -= 1 path.append((i, j)) return path[::-1] # 反转使路径从起点开始 path = trace_path(gamma)

现在让我们可视化这条路径：

plt.figure(figsize=(8,6)) plt.imshow(gamma, origin='lower', cmap='viridis') plt.colorbar(label="累积距离") plt.plot([p[1] for p in path], [p[0] for p in path], 'r', linewidth=2) plt.xlabel("序列2索引") plt.ylabel("序列1索引") plt.title("最优规整路径") plt.show()

红色路径展示了两个序列间的最佳对齐方式。路径的走向反映了时间轴的压缩和拉伸——水平移动表示序列1的一个点对应序列2的多个点（时间拉伸），垂直移动则相反。

4. 时间规整对齐可视化

最直观的理解方式是将两个序列按照最优路径进行对齐展示：

def plot_alignment(seq1, seq2, path): fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(12,8), sharex=True) # 绘制原始序列 ax1.plot(seq1, 'b-', label='序列1') ax1.plot(seq2, 'g-', label='序列2') ax1.legend() ax1.set_title("原始序列对比") # 绘制对齐后的序列 aligned_seq1 = [seq1[i] for i,j in path] aligned_seq2 = [seq2[j] for i,j in path] ax2.plot(aligned_seq1, 'b-', label='序列1(对齐后)') ax2.plot(aligned_seq2, 'g-', label='序列2(对齐后)') # 添加对应关系连线 for step in range(0, len(path), 10): # 每隔10个点画一条线 i, j = path[step] ax1.plot([i, j], [seq1[i], seq2[j]], 'r--', alpha=0.3) ax2.legend() ax2.set_title("规整对齐后的序列") plt.tight_layout() plt.show() plot_alignment(seq1, seq2, path)

上图中，第一个子图展示了原始序列，红色虚线显示了关键点之间的对应关系。第二个子图则展示了按照DTW路径对齐后的序列——你可以看到波形特征点现在完美对齐了。

5. 高级应用与优化技巧

理解了基本原理后，让我们探讨一些实际应用中的高级技巧。

5.1 约束窗口加速计算

完整DTW的复杂度是O(nm)，对于长序列可能很耗时。通过添加约束窗口可以显著加速：

def constrained_dtw(seq1, seq2, window_size=10): n, m = len(seq1), len(seq2) gamma = np.full((n,m), np.inf) gamma[0,0] = abs(seq1[0] - seq2[0]) for i in range(1, n): for j in range(max(1, i-window_size), min(m, i+window_size)): cost = abs(seq1[i] - seq2[j]) gamma[i,j] = cost + min(gamma[i-1,j], gamma[i,j-1], gamma[i-1,j-1]) return gamma[n-1,m-1], gamma

这个版本只计算对角线附近一定范围内的单元格，复杂度降为O(nw)，其中w是窗口大小。

5.2 导数动态时间规整

对于某些应用，序列的形状比绝对值更重要。导数DTW(DDTW)先计算序列的导数：

def derivative(seq): return np.diff(seq, prepend=seq[0]) def ddtw_distance(seq1, seq2): dseq1 = derivative(seq1) dseq2 = derivative(seq2) return dtw_distance(dseq1, dseq2)

这种方法对振幅偏移和线性趋势不敏感，更适合比较形状相似性。

5.3 多维度DTW

对于多维时间序列（如3D运动捕捉数据），只需修改距离计算：

def multivariate_dtw(seq1, seq2): # seq1和seq2形状为(T, D)，D是维度数 local_dist = np.array([[np.linalg.norm(x - y) for y in seq2] for x in seq1]) # 其余部分与标准DTW相同

6. 实战案例：语音信号对齐

让我们用一个真实案例展示DTW的威力。假设我们有两段说"Hello"的录音，语速不同：

# 生成模拟语音信号 def create_voice(word, speed=1.0): t = np.linspace(0, 1, 1000) if word == "hello": sig = np.sin(2*np.pi*50*t) * (t>0.2) * (t<0.8) # 基频 sig += 0.5*np.sin(2*np.pi*120*t) * (t>0.1) * (t<0.9) # 第一共振峰 sig += 0.3*np.sin(2*np.pi*240*t) * (t>0.3) * (t<0.7) # 第二共振峰 return sig[::int(1/speed)] # 通过采样模拟语速变化 normal = create_voice("hello", 1.0) fast = create_voice("hello", 1.5) plt.figure(figsize=(12,4)) plt.plot(normal, label="正常语速") plt.plot(fast, label="快速语音") plt.legend() plt.title("不同语速的语音信号") plt.show()

计算DTW对齐：

distance, gamma = constrained_dtw(normal, fast, window_size=50) path = trace_path(gamma) plt.figure(figsize=(10,8)) plt.imshow(gamma, origin='lower', cmap='viridis', aspect='auto') plt.plot([p[1] for p in path], [p[0] for p in path], 'r', linewidth=2) plt.colorbar(label="累积距离") plt.title("语音信号对齐路径") plt.xlabel("快速语音帧") plt.ylabel("正常语速帧") plt.show()

从路径可以看出，DTW成功找到了非线性对应关系，将快速语音的压缩部分与正常语速的展开部分正确匹配。