简介：一个基于TDOA时延估计与广义互相关算法的六麦克风三维声源定位MATLAB脚本，面向声源定位、麦克风阵列信号处理初学者及研究人员。资源利用声波到达时间差构建超定方程，通过广义互相关加权抑制噪声干扰，提升时延估计精度，最终解算声源三维坐标，可迁移至智能语音设备、机器人、监控系统等定位场景。压缩包内仅含一个m脚本文件，大小约2KB，代码精简，围绕信号预处理、互相关时延计算、定位方程求解等主流程展开，便于逐行研读与二次开发。当前已有224人学习下载，适合作为算法验证、课程设计或毕业设计参考。通过这份代码，读者可快速掌握时差定位核心思路，并能结合真实麦克风阵列改造参数，在实际环境中测试定位效果。

1. 六麦克风TDOA时延估计与三维定位：先算对时间，再谈空间

很多人以为六麦克风阵列买回来就能做三维定位，实际上算法链路里最花时间的是TDOA时延估计这一步。这套方案的主线很清楚：先用广义互相关（GCC-PHAT）把声音到达每只麦克风的时间差算准，再把时间差乘以声速换算成距离差，最后解一个三维非线性方程组，得到声源的 x/y/z 坐标。它适合做机器人听觉的工程师——比如人形机器人头部的麦克风阵列、会议室定向拾音——解决的不只是“听到”，而是“知道声音从哪个方向、多远来”。一个反直觉的结论是：时延估计差一个采样点，在 48kHz 采样下空间误差约 0.7cm，在 8kHz 下直接到 4.3cm。所以这类三维定位的地基，不是降噪算法，而是时延差本身的精度。

2. 广义互相关GCC-PHAT做时延估计：公式、代码与可信度判断

TDOA 的第一步是时延估计。六只麦克风两两配对能算出 15 组时延，但实际定位用不了这么多，取 5~6 组高质量的就够。时延估计的质量直接决定后面三维坐标的误差。这一章从原理讲到一份可以直接跑通的 Python 实现，再讲怎么判断一个时延值能不能信。

2.1 时延估计在算什么：互相关、峰值与PHAT加权

互相关的本质是把两路信号滑动对齐，看哪个延迟 τ 下它们最像：

R(τ) = ∫ x₁(t) · x₂(t - τ) dt

R(τ) 的最大值位置就是两路麦克风之间的到达时间差。但直接做时域互相关在语音场景里很容易翻车，因为语音的功率谱很不平坦，共振峰集中的频段能量高，会把互相关峰值“拉”到那些强频段上去，时延就偏了。

广义互相关把计算搬到频域：

R(τ) = ∫ W(f) · X₁(f) · X₂*(f) · e^{j2πfτ} df

其中 X₁(f)X₂*(f) 是互功率谱，W(f) 是加权函数。把 W(f) 取成 1/|X₁(f)X₂*(f)|，就是 PHAT 加权——每一根谱线的幅度都被归一化，只保留相位信息。这正好避开语音谱不平坦的问题：时延信息均匀地来自所有频段，主峰比普通互相关尖锐得多。这也是 GCC-PHAT 成为 TDOA 定位默认算法的原因。

PHAT 的边界在低信噪比下暴露：噪声也会被归一化，等效于把本底噪声放大。所以我一般会在 PHAT 之前先做带通滤波，把 150Hz 以下工频和 4kHz 以上混叠风险频段直接切掉，别让 PHAT 在噪声上“硬算”。

2.2 一份可直接复现的GCC-PHAT代码（Python/numpy）

下面是我常用的 GCC-PHAT 实现，去掉依赖只剩 numpy，核心逻辑全部摊开：

import numpy as np def gcc_phat(x1, x2, fs, nfft=None): """ GCC-PHAT 时延估计 返回: tau(秒), sharpness(峰值/次峰比) 约定: tau > 0 表示 x1 比 x2 晚到 """ # 1. nfft 取 2 的幂，保证 FFT 效率 if nfft is None: nfft = int(2 ** np.ceil(np.log2(len(x1)))) # 长度不一致时截短到相同长度 n = min(len(x1), len(x2)) x1, x2 = x1[:n], x2[:n] # 2. 去直流 + 加窗，抑制帧边缘泄漏 x1 = x1 - np.mean(x1) x2 = x2 - np.mean(x2) win = np.hanning(n) x1 = x1 * win x2 = x2 * win # 3. 互功率谱 + PHAT 加权 X1 = np.fft.fft(x1, nfft) X2 = np.fft.fft(x2, nfft) cross = X1 * np.conj(X2) eps = 1e-6 phat = cross / (np.abs(cross) + eps) r = np.fft.ifft(phat, nfft).real # 广义互相关函数 # 4. 峰值索引换算成时延（秒） idx = int(np.argmax(np.abs(r))) if idx >= nfft / 2: tau = (idx - nfft) / fs # 负时延落在右半边 else: tau = idx / fs # 正时延落在左半边 # 5. 主峰/次峰比，反映互相关峰是否尖锐 mags = np.sort(np.abs(r))[::-1] sharpness = mags[0] / (mags[1] + 1e-9) return tau, sharpness

逻辑说明：第 3 步是核心，cross = X₁·X₂* 保留了两路信号在各频点上的相位差；除以幅度加 eps 就是 PHAT 加权。eps 太小会在静音帧上放大数值噪声，太大又会让加权失效，我一般取 1e-6~1e-4。第 4 步的索引换算容易错：ifft 返回的是循环互相关，索引 0 是零延迟，左半边是正延迟，右半边是负延迟，所以 idx 超过 nfft/2 时要先减掉 nfft 再除以 fs。这个正负号约定放在函数头注释里，避免后续定位方程整体镜像。

参数说明：fs 直接决定时延量化粒度，48kHz 下每个采样点约 21μs，对应声程差 0.7cm；16kHz 下则是 2.1cm。nfft 默认取输入长度向上到 2 的幂，工程里我固定用 50ms 帧长、48kHz 采样，nfft 就是 4096。sharpness 返回值后面定位前筛帧要用，一帧互相关峰越尖锐越可信。

2.3 拿到时延先别急着定位：峰值可信度与正负号约定

六麦会得到很多组时延，但每一组的质量参差不齐。我的做法是先用 sharpness 筛一遍：大于 2 的帧直接进定位；1.5~2 的帧标记为弱帧，可以用但要降低权重；小于 1.5 的基本是混响严重或信噪比太低，直接丢。一堆弱帧宁可丢掉，也不要硬塞给定位解算——三维解算对异常值极敏感，一个错误的时延能把坐标带偏几十厘米。

另一个必须统一的是正负号约定。我代码里 tau 表示“x1 相对 x2 延迟多少秒”，tau > 0 就是 x1 晚到。不同资料里的约定可能相反，混用时定位结果会整体镜像翻转，而且很难排查。我习惯统一成以参考麦 0 为基准：tau_i = 到达麦克风 i 的时间 - 到达参考麦 0 的时间，这样后续距离差 d_i = c·tau_i 天然就是“第 i 个麦比参考麦远多少米”，正负永远一致。

3. 六麦克风阵列布局与三维坐标解算：从时延差到x/y/z的完整链路

时延算完，接下来面对两件事：阵列几何布局和方程组解算。麦克风坐标写错，后面算法再对也是白干；解算方程选错初值，结果会收敛到镜像位置。这一章把坐标设计、最小二乘代码和六麦冗余三个问题一次讲透。

3.1 空间混叠约束与六麦坐标设计

麦克风阵列有一个和采样定理对应的空间混叠约束：麦克风间距 d 必须小于最高频率波长的一半，否则相位差会混淆，时延峰跳到相邻周期上。约束式是 d < c/(2·fmax)。按语音 4kHz 上限、声速 343m/s 计算，d 必须小于 4.3cm。很多人上来就把麦间距拉到 10cm，结果 4kHz 附近的高频分量混叠，时延时大时小，定位跟着抖。

在“要三维定位”的前提下，我常用的六麦布局是“底面方形四麦 + 中心上方两麦”，坐标如下：

底面方形的 4cm 边长保证了 x/y 平面时延差足够明显；4、5 号麦放在阵列中心正上方形成垂直基线，让 z 轴有独立的时间差来源。三维定位必须有垂直基线，纯平面阵列对高度变化几乎不敏感，这是很多“六麦全放一个平面”的方案到最后 z 轴飘得没法用的原因。均匀圆阵在二维测向里常用，但拿到三维定位里不如“平面方形 + 垂直柱”直观——垂直方向时延信息更容易被 GCC 捕捉，音视频算法里迭代也简单。

3.2 从TDOA方程组到三维解算：最小二乘代码

设声源坐标 s = [x, y, z]，麦克风 i 的位置 pᵢ 已知，那么 TDOA 给出的距离差方程为：

||s - pᵢ|| - ||s - p₀|| = c·τᵢ = dᵢ

这是双曲面方程的几何意义：声源到两个麦克风的距离差固定，轨迹是双曲面。麦克风越多方程越多，但每个方程都有噪声，所以用最小二乘而不是直接解非线性方程组。

import numpy as np from scipy.optimize import least_squares # 与 3.1 表格一致，单位：米 MICS = np.array([ [0.00, 0.00, 0.00], # 0 号：参考麦 [0.04, 0.00, 0.00], [0.04, 0.04, 0.00], [0.00, 0.04, 0.00], [0.02, 0.02, 0.05], [0.02, 0.02, 0.10], ]) M0 = MICS[0] def tdoa_residuals(s, d, c=343.0): """理论距离差 - 实测距离差""" s = np.asarray(s, dtype=np.float64) res = [] for i in range(1, len(MICS)): dist_i = np.linalg.norm(s - MICS[i]) dist_0 = np.linalg.norm(s - M0) res.append((dist_i - dist_0) - d[i]) # d[i] = c * tau_i return res def locate_tdoa(d, x0=None, c=343.0): """ d: 长度为 5 的数组，d[i] 对应第 i 号麦克风与参考麦 0 的声程差 返回: 声源坐标, 残差平方和 """ if x0 is None: x0 = M0.copy() # 用参考麦位置做初值 result = least_squares(tdoa_residuals, x0, args=(d, c), method='lm', max_nfev=500) return result.x, result.cost

逻辑说明：残差函数把每个麦克风对参考麦 0 的理论距离差减去实测距离差，least_squares 负责让残差平方和最小。d 的下标必须和 MICS 对应——d[1] 是 1 号麦与 0 号麦的声程差，d 长度是 5，因为参考麦 0 自己不进方程。数组错位是定位翻车的高频原因，我调试时第一步永远是打印 d 和 MICS 逐行核对。

参数说明：method='lm' 适合这种小规模无约束优化，收敛快；如果后面加了高度边界约束，要换成 method='trf'。max_nfev=500 是保护上限，防止异常时延把优化器拖死。声速 c 我默认 343m/s（20℃），更精确可以用 c ≈ 331.3·√(1+T/273.15)，T 是摄氏温度；温湿度对声速的影响通常在 1% 以内，但做精密测量时温度修正能明显缩短系统偏差。

初值问题值得多说一句：这个方程组不是线性的，存在局部极小。x0 = 参考麦位置在“声源离阵列 1 米以内”时通常能收敛，但声源在阵列侧面或后方时可能收进镜像解。我一般会把四个角麦的坐标轮流当 x0 各算一遍，取 cost 最小的那个解，成本就是多跑几次 LM，效果立竿见影。

3.3 为什么选六麦克风：冗余、容错与几何精度因子

理论上四麦就能解三维坐标：3 个独立 TDOA 方程对 3 个未知数。但真实场景里混响会污染某一路时延，GCC 可能整个判错。六麦多出来的方程价值在于容错：冗余方程让最小二乘可以把残差摊薄，而且解算后看 cost，明显偏大就能判断这一帧不可信，直接丢掉。

我的经验是：从 15 组可选时延里取哪几组参与解算，比算法本身更影响结果。做法是按 sharpness 排序，只取最锐的 5~6 组构造方程，而不是贪多求全；间距过近的麦克风对远端声源几乎没有时延差，加进去只会稀释方程质量。底面方形四麦之间以及它们和上方柱麦之间的时延通常最锐，优先用。

最后提一个几何层面的限制：声源落在阵列包络内部，或者与阵列平面共面时，几何精度因子很差，微小的时延误差会被空间放大好几倍。工程上的经验做法是让声源待在阵列外侧、和阵列平面成一个角度，而不要把阵列包着声源用——这条约束对三维定位的体感影响非常大。

4. 采集链路与预处理参数：同步采样、近场模型和帧级参数怎么定

TDOA 定位的误差源头往往不在算法，而在前端采集链路。通道时钟不同源、采样率不够、帧参数不合适，都会在时延上留下系统性偏差。这一章讲清楚采集方案和参数选择，再给一套最小预处理流程。

4.1 多通道同步采样：一口时钟、TDM与ES8388级联方案

时延估计要求所有通道共享同一个采样时刻基准。严格说，这个要求等价于“所有通道的采样时钟必须同源”。很多人用四只 USB 麦克风拼成一个阵列，每一只都有自己的晶振和 USB 控制器，这叫异步采集，时延会随时间漂移，完全不能用于 TDOA。

常见做法分两条路线。上位机方案是用一台多通道同步输入的音频接口把六路模拟麦接满，或直接用支持 I2S/PDM 多通道同步采集的板卡，保证所有通道从同一主时钟派生。嵌入式方案里，我一般用多片 ES8388 这类音频 codec 级联：把它配成 TDM slave，三片级联覆盖六通道，MCLK 和 BCLK 统一，六路数据依次塞进同一条 TDM 总线，采样率完全一致。这个方案比 USB 声卡好排查问题，MCLK 有没有、BCLK 有没有，对着示波器一眼就能看出来。

同步自检有一个笨但有效的办法：用同一路正弦波通过 Y 型线同时送入两个输入通道，跑 GCC，时延应该在 ±1 个采样点内。如果稳定偏离某个固定值，是硬件走线延迟；如果逐帧漂移，就是时钟不同源。先把这条过了，再谈下面的参数。

4.2 帧级参数：采样率、帧长、滤波与VAD门控

采样率选 48kHz 是因为时延量化误差直接变成空间误差：1/48000 秒乘以声速 343m/s 约 0.7cm。16kHz 下翻三倍到 2.1cm，三维定位的坐标会明显毛糙。帧长上，50ms 在 48kHz 下是 2400 个采样点，互相关统计量足够；如果做移动声源跟踪，帧长取短端，避免声源在帧内移动导致时延随时间变化。

带通滤波上限 4000Hz 正好压住 3.1 节算出的空间混叠频率。4cm 间距的目标上限约 4.3kHz，把 4kHz 以上切掉，等于从源头上消除了混叠风险频段；下限 150Hz 是为了滤除工频 50Hz 及其谐波，防止电网噪声在 PHAT 加权里被无差别放大。

4.3 近场模型：三维定位为什么不能退化成远场方向余弦

麦克风阵列近场不是一个可选项，而是三维定位的默认前提。TDOA 距离差方程 ||s-pᵢ|| - ||s-p₀|| = dᵢ 本身就是球面波距离差形式，天然支持近场场景。所谓近场，指声源到阵列的距离与阵列孔径可比，波前是球面波而不是平面波。

近场判据是 r < 2D²/λ。按本方案 D 约 0.1m、λ 按 4kHz 的 0.086m 计算，2D²/λ ≈ 0.23m——声源离阵列 23cm 以内，平面波假设彻底失效；即使到 1m 距离，球面波和平面波在距离差上的差异依然有几十毫米到几厘米。所以只要目标写了“三维定位”，老老实实用距离差方程组；只有明确只需要方位角时，才退化成远场方向余弦。那套方向余弦解法在远场雷达里是标准，但室内声音场景声源往往就在一两米内，距离维一丢，“三维”就退化成“二维测向”，z 轴直接失去意义。

4.4 给一套最小预处理与单帧时延流程

import numpy as np from scipy import signal def process_frame(frame, fs, rms_threshold=0.01): """ frame: (6, samples) 短时帧 返回: d(声程差数组), sharpness；静音帧返回 None """ # 1. 静音检测：声压不够直接丢帧，省后续算力 rms = np.sqrt(np.mean(frame**2, axis=1)) if np.max(rms) < rms_threshold: return None # 2. 150~4000Hz 带通滤波，使用 filtfilt 保持零相位 b, a = signal.butter(4, [150/(fs/2), 4000/(fs/2)], btype='band') frame_f = signal.filtfilt(b, a, frame, axis=1) # 3. 去直流、加 hann 窗 frame_f = frame_f - np.mean(frame_f, axis=1, keepdims=True) win = np.hanning(frame.shape[1]) frame_win = frame_f * win[None, :] # 4. 以 0 号麦为参考，逐对跑 GCC-PHAT d = np.zeros(5) c = 343.0 for i in range(1, 6): tau, sharp = gcc_phat(frame_win[0], frame_win[i], fs) d[i-1] = c * tau # 声程差，单位：米 return d, sharp

逻辑说明：VAD 放在滤波之前，可以在静音帧上省掉滤波和 FFT 的开销。filtfilt 是零相位滤波，不会给信号引入与频率相关的群延迟；如果图省事用普通 butter 滤波，相位会随频率弯曲，时延整体偏移，这类误差在 GCC 里特别隐蔽，现象是“低音时延准、高音时延偏”。win[None, :] 是把 (samples,) 的窗函数广播到 (6, samples) 数组上。d[i-1] 对应第 i 号麦相对参考麦 0 的声程差，直接喂给 3.2 节的 locate_tdoa。

5. 六麦克风TDOA定位避坑指南：4个让三维坐标翻车的实测问题

这一章记几条血泪经验。每次看到有人调了两周 GCC 还在怀疑参数，多半是踩了下面四个坑里的某一个。

5.1 混响把互相关峰值拉偏

现象：吸音室测试定位很好，搬到办公室后开始跳变；时延数组偶发偏离大半毫秒，定位结果突然跑偏 20cm。互相关主峰还在，但次峰已经和主峰差不多高。

原因：反射声叠加在直达声上，互相关函数出现多个与房间反射路径相关的旁瓣。GCC-PHAT 的混响鲁棒性比普通互相关好很多，但混响能量过强时，主峰仍然可能被某个反射峰压制。

解决：第一，把频带压缩到 150~4000Hz，只让有效语音能量进入 GCC；第二，按 2.3 节的 sharpness 阈值丢弃弱帧，弱帧硬解只会带偏整体结果；第三，如果还不够，改用平方根归一化的 Sqrt-PHAT，或做自适应频带加权，降低强混响频段的贡献；第四，工程上有一个后悔药：连续 N 帧定位结果都在同一坐标附近异常跳变时，直接保持上一帧有效位置，等互相关峰恢复尖锐再更新。这个策略不好看，但现场体验提升明显。

5.2 采样时钟不同源：时延漂移是线性的

现象：固定声源不动，上午测误差 1cm，下午误差 7cm；把误差随时间画出来，接近一条直线，斜率稳定。

原因：两路通道的采样率存在几十 ppm 偏差。以 40ppm 为例，每秒产生 40μs 的时延误差，按声速 343m/s 折算成距离就是每秒累积 1.4cm。这是随运行时间线性增长的误差，降噪、滤波一概无效。

解决：换同一时钟域的同步采集方案是最彻底的。如果硬件已经定死，退而求其次是每次定位前用固定位置参考声源重新标定时钟偏差，然后用软件重采样把各通道对齐到同一时间轴。重采样前必须先测出偏差量，方法就是 4.1 节那个“一路正弦同时送两通道”的自检流程。

5.3 阵列坐标误差被TDOA方程放大

现象：sharpness 很高、时延看起来全对，但定位结果整体偏向某个方向 3~5cm，怎么调算法都纹丝不动。

原因：麦克风实际安装位置和代码 MICS 数组里写的坐标差了 1~2mm。近距离声源下，坐标误差会沿着双曲面方程梯度被放大，尤其是声源靠近阵列时，毫米级误差变成厘米级偏差，这就是几何精度因子在起作用。

解决：用激光测距或游标卡尺量麦克风振膜中心位置，把真实坐标写进 MICS 数组；不要信外壳开孔位置，振膜中心才是声学中心。阵列固定在刚性板上之后，应该再复测一次坐标。这一步成本最低，但最容易被跳过——很多人调算法调了两周，最后发现是坐标差了 1.5mm。

5.4 近场当远场算：z轴飘了20cm

现象：x/y 方向定位凑合能用，z 轴高度误差特别大，而且声源低于或高于阵列平面时，误差进一步放大。

原因：算法里用了远场平面波/方向余弦假设，把距离差方程近似成了角度问题。近场场景下，被近似丢掉的球面波修正项恰好包含高度信息，z 轴自然失准。

解决：回到 4.3 节的近场模型，用距离差方程组解 s。另一个隐蔽因素在垂直基线本身：本方案垂直方向有效基线只有 5~10cm，时延误差 1 个采样点（48kHz 下 21μs）对应约 7mm 声程差，反映在 z 轴上可能被放大到数倍。所以垂直麦对的时延要选 sharpness 最高的帧来解，再在 z 轴输出上加一节平滑滤波，轨迹会稳得多。

6. 验证闭环与进阶技巧：用chirp信号量化定位误差

先做端到端验证：放一个扬声器在已知坐标，播放 0.5 秒、200Hz 到 4kHz 的对数扫频 chirp，跑完整个链路，重复 10 次。置信坐标均值误差是系统误差（阵列坐标、声速、硬件延迟），标准差是随机误差（互相关噪声、风噪）。我给内部定的合格线是均值偏差小于 5cm、标准差小于 3cm；超差时按“声速温度修正 → 复测麦克风坐标 → 检查通道同步”的顺序排查，这三样占了系统误差的大头。

import numpy as np from scipy import signal fs = 48000 t = np.arange(int(0.5 * fs)) / fs # 对数扫频 chirp，能量从 200Hz 平移到 4kHz chirp = signal.chirp(t, f0=200, f1=4000, method='logarithmic') chirp = chirp * np.hanning(len(t))

chirp 信号的好处是频率连续覆盖目标频段，互相关峰比语音更尖锐，验证时延链路时不容易被共振峰干扰。注意播放器和功放会引入固定延迟，这部分会整体平移时延，在验证时先测出来、在距离差里减掉，再谈绝对坐标精度。

再往后做移动声源跟踪，我习惯在定位结果上套一个位置-速度模型卡尔曼滤波，把每帧独立解算的坐标变成平滑轨迹；同时用 3.2 节的 cost 做异常值门控，cost 突增时让卡尔曼忽略该帧。多目标场景可以加 RANSAC 式的多次解算投票，但算力消耗要提前评估。

最后分享一个我自己的调试习惯：先仿真、后实采。用已知时延加高斯噪声的合成波形先跑通 GCC 和最小二乘，把算法侧的 bug 全部清干净，再接入真实阵列；装完阵列先做同步自检，最后放 chirp 做端到端验收。这个顺序能避免一个问题——当定位不准时，你分不清是算法错了、坐标错了还是时钟错了。说实话，我在这个方向上栽过最大的跟头，就是跳过同步自检直接去调滤波参数。后来养成了先同步自检、再仿真、后实采的顺序，定位问题的排查时间基本能砍一半。这套六麦克风 TDOA 三维定位的链路不算新，但它对细节的敏感度极高，希望这些参数和避坑经验帮到你。

本文还有配套的精品资源，点击获取