相干函数：从概念到实战，信号关系量化与工程诊断指南

1. 从“相干”说起：一个被误解的工程概念

在信号处理、振动分析乃至光学领域，“相干”这个词出现的频率相当高。但很多刚接触的朋友，包括我当年，都容易把它想得过于玄乎——是不是意味着两个信号“一模一样”？或者“完全同步”？其实不然。我第一次在实验室里对着两台并不同步的振动台数据发愁时，导师就指着屏幕上的两条曲线说：“别光看它们长得像不像，得看它们‘相干’不‘相干’。” 这句话点醒了我。相干函数，本质上是一个统计工具，它回答的核心问题是：在给定的频率上，两个信号之间到底存在多少线性的、因果性的关系？它衡量的是“ predictability”（可预测性）。如果你能通过一个信号完美地预测另一个信号在某个频率成分上的振幅和相位，那么它们在这个频率上就是完全相干的；如果完全无法预测，那就是不相干。

这个概念听起来有点抽象，但它的应用场景却极其广泛且接地气。比如，你想知道汽车在颠簸路面上行驶时，方向盘抖动的“罪魁祸首”是前轮还是后轮的振动？又比如，在一间嘈杂的会议室里，你如何判断两个麦克风录到的噪音是来自同一个声源（比如空调）还是各自独立的干扰？再比如，在结构健康监测中，如何确认桥梁某处的异常振动响应，确实是由特定方向的风载荷引起的，而不是传感器自身故障？这些问题，都可以通过计算和分析相干函数来寻找答案。它就像一个精密的“关系鉴定器”，帮我们从混杂的、看似相关的信号中，剥离出真正有因果联系的部分。对于从事测试、诊断、噪声控制、通信系统设计的工程师和研究人员来说，掌握相干函数是进行有效数据分析的基本功。接下来，我就结合自己踩过的坑和实战经验，把这个工具从里到外拆解清楚。

2. 相干函数的核心原理与数学内涵拆解

要真正用好一个工具，不能只停留在“点按钮出结果”的层面，必须理解其背后的数学逻辑和物理意义。否则，当计算结果出现反常（比如相干值大于1）时，你根本无从判断是数据出了问题，还是自己用错了方法。

2.1 定义式：从互功率谱密度出发

相干函数最经典的定义基于傅里叶变换和功率谱。假设我们有两个平稳随机信号x(t)和y(t)。我们关心的不是它们时域上的具体波形，而是其频域上的关系。

首先，我们定义互功率谱密度：G_xy(f) = E[X*(f) Y(f)]这里，X(f)和Y(f)分别是x(t)和y(t)的傅里叶变换，*表示复共轭，E[]表示期望值（在实操中通常用多次平均来近似）。G_xy(f)是一个复数，它包含了两个信号在频率f上的幅值关系信息和相位差信息。

然后，我们定义自功率谱密度：G_xx(f) = E[X*(f) X(f)]（信号x自身的能量谱）G_yy(f) = E[Y*(f) Y(f)]（信号y自身的能量谱）

那么，常相干函数（Ordinary Coherence Function）γ²_xy(f)定义为：γ²_xy(f) = |G_xy(f)|² / (G_xx(f) * G_yy(f))

这个公式是理解一切的核心。我们来逐项解读：

1. 分子|G_xy(f)|²：这是互功率谱密度幅值的平方。它衡量的是两个信号在频率f上“共同”的功率有多大。如果两个信号在该频率上完全独立（比如一个是信号，一个是纯随机噪声），那么经过多次平均后，它们的互功率谱会趋近于零。
2. 分母G_xx(f) * G_yy(f)：这是两个信号各自功率谱的乘积。它代表了一个“理论上的最大可能共同功率”。可以想象，如果两个信号一模一样，那么|G_xy(f)|²就会等于G_xx(f) * G_yy(f)。
3. 比值γ²_xy(f)：因此，相干函数的值被归一化到 0 和 1 之间。
   • γ²(f) = 1：表示在频率f上，y(t)的能量可以完全由x(t)通过一个线性系统来解释。两者在该频率成分上存在完全确定的线性关系。注意，这并不意味着它们振幅相等，只意味着存在一个固定的复数增益（包含幅值缩放和相位移动）将x(f)映射到y(f)。
   • γ²(f) = 0：表示在频率f上，y(t)与x(t)完全不相关。y(t)在该频率上的能量全部来自噪声、非线性失真或其他与x(t)无关的源。
   • 0 < γ²(f) < 1：这是最常见的情况。意味着y(t)中该频率成分的一部分来源于x(t)的线性传递，另一部分则来源于其他无关的噪声或非线性成分。

注意：这里用的是γ²（相干系数的平方），很多软件（如MATLAB的mscohere函数）直接输出的是平方后的值。也有些文献和软件输出γ（模值）。在引用或对比数据时，务必弄清楚是平方值还是模值，两者数值和意义差别很大。我个人的习惯是统一使用γ²，因为它直接代表了线性相关功率的百分比。

2.2 物理意义：它到底在衡量什么？

理解公式后，我们可以从三个层面把握其物理意义：

1. 线性度的度量：相干函数接近于1，强烈暗示两个信号在该频率上的关系是线性的。如果系统中存在显著的非线性（如饱和、截止、谐波失真），即使输入输出信号看起来相关，相干函数也会在非线性影响的频率上下降。
2. 噪声影响的度量：相干函数下降，直接表明测量到的输出信号中，混入了与输入信号不相关的噪声。这些噪声可能来自测量环境、传感器本身、或其他未被监控的输入源。
3. 因果关系的指示器（需谨慎）：高相干性通常被解释为存在因果关系（x导致y）的证据。但这并非绝对。如果存在第三个共同源z同时驱动x和y，即使x和y之间没有直接因果联系，它们也可能表现出高相干性。这是解释相干结果时最常见的陷阱之一。

2.3 关键参数：估计精度与平均次数

在实际计算中，我们是用有限长度的采样数据去估计理论上无限长的随机过程的统计特性。这就引入了估计误差。相干函数的估计值Ĝ²_xy(f)是一个随机变量，其统计特性与一个参数强相关：平均次数n_d。

当我们进行谱估计时，通常会将长数据分段（或加窗重叠分段），对每一段计算周期图，然后对所有段的周期图进行平均，从而得到平滑的谱估计。这个平均的段数，就是n_d。

相干函数估计值的偏差和方差近似为：

• 偏差：E[Ĝ²_xy] ≈ γ²_xy + (1 - γ²_xy)/n_d。即使真实相干为0，估计值也可能有1/n_d的正偏差。这意味着，如果平均次数太少，你可能会在不相干的频率上观察到虚假的、非零的相干值。
• 方差：随n_d增大而减小。

实操心得：n_d的选择是精度和频率分辨率的权衡。更多的平均次数（n_d更大）会降低估计方差，使相干曲线更平滑，置信度更高，但会降低频率分辨率（因为每段数据变短了）。我的经验法则是：确保n_d至少大于 20。对于初步分析，我常设置n_d使得最终相干谱的频率分辨率（Δf = fs / N_fft，其中N_fft是每段FFT点数）能满足工程需求，然后尽量增加平均段数。在MATLAB中，使用mscohere(x, y, window, noverlap, nfft, fs)函数时，可以通过调整window长度、noverlap（通常设50%重叠率以增加段数）和nfft来控制这些参数。

3. 实战案例解析：从数据采集到结果解读

光说不练假把式。我们通过一个经典的机械振动传递路径分析的案例，来走一遍完整的流程。这个案例模拟了工业现场常见的问题：一台设备（如泵、风机）的振动传递到了基座和周边结构，我们需要量化源头振动与结构响应之间的关系。

3.1 案例背景与数据模拟

假设我们有一台旋转设备，在某个转速下运行，其转频为f_r = 25 Hz。我们在设备轴承座（输入点x）和远处一个关键的结构测点（输出点y）上分别安装了加速度传感器。我们怀疑结构测点的振动主要来源于设备轴承的振动传递。

为了模拟真实情况，我们构造以下数据：

1. 输入信号x(t)：一个以f_r为主频的正弦波，模拟设备的不平衡力，同时叠加一些宽频带的随机噪声，模拟轴承的滚动体噪声。x(t) = A1 * sin(2π * f_r * t) + A2 * randn(size(t))
2. 传递路径：假设从x到y是一个线性系统，其频率响应函数为H(f)。在f_r处有一个共振峰。
3. 输出信号y(t)：由x(t)经过线性系统H(f)产生，但同时，在结构测点y本身还存在一个独立的环境噪声源（比如远处其他设备的干扰），其频率集中在f_noise = 60 Hz附近。y(t) = IDFT( H(f) * DFT(x(t)) ) + A3 * sin(2π * f_noise * t + φ) + A4 * randn(size(t))其中IDFT/DFT表示逆/离散傅里叶变换，φ是随机相位。

数据采集关键设置：

• 采样频率fs = 500 Hz（满足奈奎斯特定律，远高于关心的最高频率）。
• 采样时长T = 20 s，以获得足够的数据进行平均。
• 两个通道严格同步采样，这是计算互谱和相干的前提。如果使用两个独立的采集设备，必须有时钟同步信号。

3.2 计算步骤与MATLAB/Python实操

这里以MATLAB为例，Python的scipy.signal库中coherence函数用法类似。

% 1. 模拟生成数据 fs = 500; % 采样率 (Hz) T = 20; % 总时长 (s) t = 0:1/fs:T-1/fs; % 时间向量 N = length(t); % 数据总点数 f_r = 25; % 转频 (Hz) f_noise = 60; % 干扰噪声频率 (Hz) % 生成输入信号 x A1 = 2; A2 = 0.5; x = A1 * sin(2*pi*f_r*t) + A2 * randn(1, N); % 设计一个简单的模拟传递路径 H(f)（二阶系统，在f_r附近有增益） fn = f_r; zeta = 0.05; % 为了方便，我们在时域用滤波器近似。设计一个带通滤波器，中心在f_r附近。 [b, a] = butter(4, [f_r-2, f_r+2]/(fs/2), 'bandpass'); x_filtered = filter(b, a, x); % 生成输出信号 y：包含滤波后的x、独立的60Hz干扰、附加随机噪声 A3 = 1.5; A4 = 0.3; independent_noise = A3 * sin(2*pi*f_noise*t + rand*2*pi); y = x_filtered + independent_noise + A4 * randn(1, N); % 2. 计算相干函数 window = hann(round(fs*2)); % 窗函数：每段2秒，即1000点 noverlap = round(length(window)/2); % 50%重叠 nfft = max(2^nextpow2(length(window)), 1024); % FFT点数，不小于窗长 [coh, f] = mscohere(x, y, window, noverlap, nfft, fs); % 3. 计算功率谱密度作为对比 [Pxx, f_psd] = pwelch(x, window, noverlap, nfft, fs); [Pyy, f_psd] = pwelch(y, window, noverlap, nfft, fs); % 4. 绘图 figure('Position', [100, 100, 1200, 800]) % 子图1：时域信号片段 subplot(3,1,1) plot(t(1:fs*2), x(1:fs*2), 'b', 'LineWidth', 1); hold on; plot(t(1:fs*2), y(1:fs*2), 'r', 'LineWidth', 1); xlabel('时间 (s)'); ylabel('幅值'); legend('输入 x(t) - 轴承座', '输出 y(t) - 结构点'); title('时域信号对比 (前2秒)'); grid on; % 子图2：自功率谱密度 subplot(3,1,2) semilogy(f_psd, Pxx, 'b', 'LineWidth', 1.5); hold on; semilogy(f_psd, Pyy, 'r', 'LineWidth', 1.5); xlabel('频率 (Hz)'); ylabel('功率谱密度 (dB/Hz)'); legend('G_{xx} - 输入', 'G_{yy} - 输出'); title('自功率谱密度'); xlim([0, 100]); grid on; % 子图3：相干函数 subplot(3,1,3) plot(f, coh, 'k', 'LineWidth', 2); xlabel('频率 (Hz)'); ylabel('相干系数 \gamma^2'); title('常相干函数'); ylim([0, 1.1]); xlim([0, 100]); grid on; hold on; % 标记我们关心的频率点 plot([f_r, f_r], [0, 1], 'b--', 'LineWidth', 1); text(f_r+1, 0.9, sprintf('转频 f_r = %d Hz', f_r), 'Color', 'b'); plot([f_noise, f_noise], [0, 1], 'r--', 'LineWidth', 1); text(f_noise+1, 0.7, sprintf('干扰噪声 %d Hz', f_noise), 'Color', 'r');

3.3 结果解读与工程诊断

运行上述代码后，我们会得到三张图。解读的关键在于第三张图——相干函数图。

1. 在转频f_r = 25 Hz处：

   • 预期：由于输出y在25Hz的成分主要来源于输入x的线性传递（我们用了滤波器模拟），因此这里的相干系数γ²应该非常接近于1。
   • 实际看图：你会观察到在25Hz附近有一个尖锐的峰值，γ²值可能达到0.95以上。这证实了我们的假设：结构测点在25Hz的振动，主要（甚至几乎完全）是由设备轴承座的振动引起的。这是一个强关联证据。

2. 在干扰噪声频率f_noise = 60 Hz处：

   • 预期：输出y中的60Hz成分来源于一个独立的噪声源（independent_noise），与输入x完全无关。因此，在这个频率上，相干系数γ²应该非常低。
   • 实际看图：在60Hz处，γ²会跌落到一个很低的值，可能接近0.1甚至更低（取决于平均次数和噪声强度）。这清晰地告诉我们：结构点在60Hz的振动问题，根源不在我们监控的这台设备上。工程师应该去寻找其他可能的60Hz干扰源（如电源谐波、其他设备的运行频率等）。

3. 在其他频段（例如宽频随机噪声部分）：

   • 预期：输入x和输出y都包含宽带随机噪声（randn部分），这些噪声成分彼此独立。因此，在这些频段，相干系数会围绕一个较低的值波动。
   • 实际看图：你会看到一条在低值附近（例如0.1-0.3）波动的“基线”。这个基线的水平与公式中的1/n_d偏差有关。平均次数n_d越大，这条基线就越贴近0。

工程诊断结论： 通过这次相干分析，我们可以向项目团队出具明确的报告：

• 主要传递路径确认：结构在25Hz的共振问题，根源是设备的不平衡力。建议的解决措施是对设备进行动平衡校正，而非盲目加固结构。
• 排除无关干扰：60Hz的振动与当前设备无关，避免了错误的维修方向。下一步应排查厂区内其他可能产生60Hz振动的设备或电气问题。
• 量化贡献度：相干系数γ²的值（如0.97）可以近似理解为“该频率成分上，输出能量中来源于输入线性贡献的比例”。这为振动贡献量排序和设定治理目标提供了数据支撑。

4. 深入应用：偏相干与多重相干——剥离复杂系统中的关联

常相干函数在处理单一输入、单一输出，且噪声干扰不大的情况下非常有效。但现实工程问题往往更复杂：多个输入源同时作用于一个输出，且这些输入之间可能还存在关联。例如，汽车行驶时，方向盘振动可能同时来源于发动机、轮胎和传动轴，而这些源之间本身可能相关（发动机振动会传递到车架再影响其他部分）。这时，常相干函数就会“失灵”或“误导”，因为它无法区分一个输出成分是来自输入A的直接作用，还是通过输入B间接产生的作用。

这就需要引入更高级的工具：偏相干函数和多重相干函数。

4.1 偏相干函数：排除其他输入影响后的“纯”关系

偏相干函数衡量的是，在排除掉其他所有已知输入信号的影响之后，某个特定输入x_i与输出y之间剩余的线性相关程度。

定义：考虑有n个输入x_1, x_2, ..., x_n和一个输出y。输入x_i与输出y的偏相干函数γ²_yxi·z，其中z代表除x_i之外的所有其他输入集合。其计算需要用到条件谱密度的概念，通常通过求解一个线性方程组（与计算多重回归残差相关）来获得。

物理意义：它回答了“如果我能完美地测量并扣除掉x_1, x_2, ..., x_{i-1}, x_{i+1}, ..., x_n对y的贡献，那么剩下的y里面，还有多少是和x_i线性相关的？” 这就像在多人讨论中，让你只听其中一个人的发言，而自动屏蔽掉其他人的声音和他们对发言者的影响。

应用场景：

• 汽车NVH诊断：在方向盘上布置多个加速度计作为输入（发动机悬置点、副车架连接点、转向柱），以方向盘振动为输出。计算发动机悬置点与方向盘的偏相干，可以排除通过副车架传递路径的间接影响，直接评估发动机振动对方向盘的“纯净”贡献。
• 建筑风振分析：多个方向的风压作为输入，建筑顶部的加速度作为输出。计算某个方向风压的偏相干，可以评估该方向风荷载对响应的独立贡献，排除其他方向风荷载的耦合效应。

实操心得：计算偏相干需要所有输入输出信号的同步测量数据，且对数据的质量和线性假设要求更高。在MATLAB中，可以使用cpsd（互功率谱密度）和pwelch函数构建谱矩阵，然后通过矩阵运算求解。也可以借助系统辨识工具箱中的函数。关键步骤是确保所有输入信号之间的互谱也被准确估计，因为偏相干计算的核心就是扣除这部分关联的影响。

4.2 多重相干函数：所有已知输入对输出的总解释力

如果说偏相干是“分拆责任”，那么多重相干就是“汇总功劳”。

定义：多重相干函数γ²_yx衡量的是，所有n个已知输入信号x_1, x_2, ..., x_n共同对输出y的线性影响程度。它表示y的总功率中，能够被这组输入信号的线性组合所解释的比例。

计算：γ²_yx(f) = (G_xy^H * G_xx^{-1} * G_xy) / G_yy其中，G_xy是输入与输出的互谱密度列向量，G_xx是输入之间的互谱密度矩阵，^H表示共轭转置，^{-1}表示矩阵求逆。

物理意义：γ²_yx(f)的值越接近1，说明在频率f上，输出y的能量几乎全部可以由我们测量到的这组输入{x_i}来解释。如果存在显著的、未被测量的输入源或强烈的非线性，多重相干就会小于1。

应用场景：

• 传递路径分析（TPA）完备性检查：在汽车NVH的TPA中，我们测量了所有主要的力输入点（如发动机悬置、排气吊挂、轮胎接地点等），并预测了车内噪声。计算预测噪声与实际测量噪声的多重相干，如果在整个频段上都接近1，说明我们的模型（传递路径和力输入）是完备的；如果在某些频段明显低于1，则提示我们可能遗漏了重要的传递路径或输入源。
• 传感器布局优化：在大型结构健康监测中，布置大量传感器不现实。可以在候选位置布点，计算不同传感器子集与关键响应点的多重相干。选择那些能以最少传感器数量使多重相干在关键频段保持高值的布点方案，实现成本与监测效果的平衡。

案例启示：在一次风机塔筒振动分析中，我们同时测量了机舱三个方向的振动（x1, x2, x3）和塔筒顶部的振动（y）。常相干显示x1与y在某个频率有高相干。但计算偏相干时，排除了x2和x3的影响后，x1与y在该频率的偏相干变得很低。而x2与y的偏相干却很高。最终发现，是x2方向的振动（机舱横向摆动）是主要根源，x1方向的高常相干只是因为机舱各方向振动存在耦合，它“蹭”了x2的关联性。这个例子深刻说明了在复杂多源系统中，直接使用常相干可能得出错误结论，必须借助偏相干进行深入剖析。

5. 常见陷阱、误区与高级注意事项

即使理解了原理，在实际操作中依然会踩坑。下面是我总结的几个关键注意事项和高级技巧。

5.1 相干值低的可能原因（不仅仅是“无关”）

看到相干值低，不要立刻下结论说“这两个信号没关系”。需要系统排查：

5.2 相干值高的解读陷阱（不一定等于“因果”）

看到相干值高，尤其是接近1，也要保持警惕：

• 共同源问题：如前所述，如果存在第三个源z同时驱动x和y，x和y会表现出高相干，但并无直接因果。例如，电网电压波动（z）可能同时导致实验室的电源（x）和照明（y）闪烁，计算x和y的相干会很高。
• 反馈回路：在闭环控制系统中，输出y会反馈回来影响输入x。此时的高相干反映的是闭环系统的整体特性，不能简单理解为x单向引起y。
• 过拟合与虚假相干：如果数据长度极短，或频率分辨率极低，可能会偶然产生高相干峰值。一定要结合物理意义判断：这个频率峰值是否对应系统已知的固有频率、转频或其倍频？

5.3 实操中的高级技巧

1. 置信区间估计：相干函数的估计值是有统计波动的。可以基于平均次数n_d和置信水平（如95%），计算出相干估计值的置信区间。如果计算出的相干值落在“零相干”的置信区间内，那么即使它不为零，也可能没有统计显著性。一些高级分析软件会提供此功能。
2. 结合频率响应函数（FRF）使用：相干函数常与FRF（H1或H2估计）一起使用。H1估计假设所有噪声都在输出端，H2估计假设所有噪声都在输入端。观察不同频率点上的相干值，可以帮助判断哪种FRF估计更可靠（相干高的频段，H1和H2估计结果一致；相干低的频段，两者差异大）。
3. 时变系统的分析：对于非平稳信号（如启动过程、变速运行），传统的基于平稳假设的相干函数不适用。可以采用短时傅里叶变换（STFT）或小波变换，计算“时频相干”或“小波相干”，来观察相干性随时间的变化。这在旋转机械的启停机分析中非常有用。
4. 循环平稳信号的相干：对于旋转机械等产生的循环平稳信号，其统计特性随时间周期性变化。直接计算普通相干会丢失信息。可以采用循环相干函数，它能揭示与循环频率相关的耦合特性，对于诊断齿轮、轴承的局部故障更为敏感。

相干函数是一个强大而精细的工具，它像一把尺子，能量化信号间线性关系的强度。但尺子本身不会思考，解读测量结果、洞察数据背后的物理机制，始终是工程师最重要的任务。从明确分析目标、精心设计测试方案、规范采集数据，到选择合适的参数进行计算，最后结合物理背景审慎解读，每一步都至关重要。避免陷入“唯数据论”的陷阱，让相干函数成为连接测试数据与工程真相的可靠桥梁，这才是掌握它的真正意义。