ARTICLE DETAIL

资讯详情

深耕网站建设与运营推广的一线实战洞察。

自相关、互相关与相干性:从数学定义到工程应用

自相关、互相关与相干性:从数学定义到工程应用 做设备振动监测的朋友拿了一段加速度信号找我说频谱毛刺太多十几万个点里根本找不到轴承故障的特征频率。我让他先把数据做一遍自相关滞后轴上的周期峰值清清楚楚故障间隔就摆在那里。他愣了半天说当年信号处理课上都学过真到现场却想不起来用。这个事让我一直觉得相关性和相干性这套东西属于那种教材里公式写得很整齐、但没几个人真正用熟的工具。这篇东西我就以自相关、互相关、卷积为线索把相关性Correlation和相干性Coherence从数学定义到工程应用完整捋一遍如果你正在做振动分析、声学定位、图像配准或者任何一种从信号里找关系的活儿应该能直接拿走用。1. 三个概念一起聊才不会越学越糊涂1.1 相关性一开始只是在回答像不像先退回到最朴素的场景。你有两组数比如同一台设备上两个测点的振动幅值序列想判断它们是不是一起涨一起落。最简单的办法是算皮尔逊相关系数一个数落在[-1, 1]区间里1表示完全正相关-1表示完全负相关0表示看不出线性关系。这个系数本身没有问题但它把两个序列当成两个整体来比较丢掉了时间结构的信息。实际信号里最常见的情况是两个序列非常相似但存在一个时间错位。同一个声音从声源传到两个麦克风距离不同到达时刻就不同如果你直接把采样点一一对应去求相关系数错开的那一段会把结果拉向0。这时候自然想到一个办法把其中一个序列往前滑动一点再算相似度把同一个错位量下算出的相关系数记录成一个函数。这个函数就是相关函数。根据参与运算的是两个不同序列还是同一个序列又分成互相关和自相关互相关衡量序列x和序列y之间的相似度随滞后量变化的函数。自相关衡量序列x和它自身滞后副本之间的相似度变化。这部分是后面所有计算的起点。很多人一上来就背公式忽略了相关函数本来就是为了解决时间错位的相似度度量这个动机遇到具体问题反而不知道怎么用。1.2 相关函数看的是整个函数不是单个数值理解相关函数的关键是记住它是一条随滞后变化的曲线而不是一个数。横轴是滞后量τ常用单位可以是采样点、秒、毫秒纵轴是相似度。如果两个信号没有可观测的时间错位互相关函数的峰值一般出现在τ0附近如果第二个信号比第一个信号晚到了20个采样点峰值就会往正方向偏移20个点这个偏移量就是时延估计的基础。写成离散形式不带归一化的互相关是Rxy[τ] Σ x[n] · y[n τ]自相关就是令y等于xRxx[τ] Σ x[n] · x[n τ]这个和会随着序列长度增大而变大所以在实际工程里我们一般做归一化处理得到取值在[-1, 1]范围内的归一化相关函数。不要小看这个归一化它能把不同幅度量级的信号拉到同一个比较尺度上否则比较两个幅值相差很大的测点信号时结果很容易被大能量信号主导。1.3 相干性解决的问题不一样哪个频段在相关时域相关函数给出的是一条随滞后变化的曲线它可以告诉你两个信号整体的相似程度但无法告诉你在10Hz附近相关还是在200Hz附近相关。如果两个信号在中频段高度相关、在高频段各自独立时域相关函数会把不同频段的行为混在一起得到一个模糊的结果。相干性Coherence就是专门回答哪个频率成分相关的工具工程上最常用的是幅值平方相干Magnitude Squared Coherence, MSC定义是Cxy(f) |Sxy(f)|² / [ Sxx(f) · Syy(f) ]其中Sxx和Syy是两个信号各自的功率谱密度Sxy是互功率谱密度。从形式上看它就是频域里的相关系数平方取值范围在[0, 1]之间0表示该频率成分完全无关1表示完全线性相关。所以整套知识其实可以分成两层相关函数在时域上刻画像不像、错开多少相干性在频域上刻画哪个频率上相关、相关到什么程度。两者各有位置谁也替代不了谁。2. 卷积与相关就差一翻——数学结构、物理含义与CNN那个假卷积2.1 卷积的定义里有个翻转动作卷积的连续形式是(f * g)(t) ∫ f(τ) · g(t - τ) dτ只要看g的变量t - τ就能发现g要先沿纵轴翻转然后随t滑动。离散卷积写成y[n] Σ x[k] · h[n - k]这里h的序号是n - k也就是说卷积算子本质上是把其中一个序列翻转后再逐点相乘求和。为什么卷积必须翻转这要回到卷积在LTI系统里的物理意义。系统的输出是输入冲击响应h的叠加当前时刻的输出取决于过去的输入所以时间轴的方向在数学上会体现为反向遍历。这个翻转不是形式上的花活而是因果性、叠加积分这些物理约束在数学结构上的投影。2.2 相关没有翻转因为它只关心对齐再看互相关Rxy[τ] Σ x[n] · y[n τ]y没有翻转只是平移了τ。原因是做相关的时候我们关心的是两个时间序列在同一个时间方向上的相似程度就像把两条波形曲线叠在一起左右拖动其中一条看重叠部分像不像。这个过程不需要翻转翻转了反而会把因果顺序弄反。把相关和卷积放在一起看两者之间的数学关系非常直接Rxy[τ] x(-t) * y(t)也就是说互相关可以理解为把x翻转后再与y做卷积。这个关系不只是形式上好看它直接提供了工程上的计算捷径卷积可以用FFT快速实现互相关同样可以用FFT快速实现。2.3 FFT实现为什么卷积和相关都能转成频域乘法卷积定理告诉我们时域卷积等于频域相乘。对离散信号循环卷积可以通过离散傅里叶变换加速y IFFT( FFT(x) · FFT(h) )相关运算同样能转成频域乘法只是其中一项需要取共轭Rxy IFFT( conj(FFT(x)) · FFT(y) )如果信号是实数序列共轭可以忽略直接乘再反变换就行。这里有一个细节使用FFT得到的是循环相关如果序列长度不够循环回绕会造成边缘污染需要做长度填充。最简单的做法是对两个序列补零到长度不小于N L - 1再进入FFT流程。速度上长度为N的序列直接算相关复杂度是O(N²)FFT方法降到O(N log N)。N大到几千甚至几万的时候这个差距是决定性的。当年我处理一段50万点的数据直接循环慢到怀疑人生换FFT方法后几乎瞬间出结果。2.4 深度学习里的卷积其实是互相关现在深度学习很热卷积神经网络CNN里到处是卷积核这个词。但细看CNN的计算方式——把卷积核在输入特征图上滑动对覆盖区域做加权求和整个过程里没有翻转核的步骤。这在信号处理的老规矩里叫互相关不叫卷积。为什么深度学习不翻转核从数学上说卷积和互相关的唯一差别就是翻转。如果用梯度下降训练网络完全可以通过调整核的权值把翻转后的结果学出来翻转与否对模型的表达能力没有影响。所以CNN这个命名其实是历史遗留真正做的是互相关运算。你在做信号处理时如果把这个概念混淆了会搞不懂为什么卷积的交换律在NN里好像不存在——因为互相关本身不满足交换律翻转的缺失让卷积核可交换这一性质消失了。我自己在带新人时有个习惯先让他们算一遍卷积再算一遍相关对比结果再让他们解释CNN为什么不翻转。这一套走通理解立刻上一个台阶不会再出现卷积核要不要翻转的经典迷思。3. 自相关把信号和自己对齐找回藏在噪声里的周期性3.1 自相关的两个关键性质自相关有一些非常重要的性质是实际应用的根本依据。第一自相关在滞后为0处取得最大值。这个由柯西-施瓦茨不等式保证互相关也可能在某个非零滞后处出现较大峰值但自相关一定是自己和自己完全对齐时最大。第二如果信号是周期的自相关也会呈周期性并在滞后等于周期整数倍的位置出现峰。这个性质非常实用周期性信号经过自相关之后周期性不但保留而且周期会以峰值间隔的形式直观呈现出来。第三纯白噪声的自相关趋近于一个冲击函数除了滞后0处有值其他滞后位置几乎为0。这意味着自相关能把淹没在噪声里的周期成分捞出来。3.2 为什么自相关能从噪声里找周期设想一个正弦信号叠加了宽带宽噪声直接看时域波形就是一团毛刺做FFT看幅值谱周期成分的谱线确实存在但可能被噪声底抬高尤其在采样点有限时谱线附近会有明显的方差波动。自相关把时间域里的噪声成分平均掉了只剩信号与自身的关联周期成分会以显著峰的形式出现。工程案例很典型轴承外圈故障会在振动信号里产生周期性冲击冲击间隔等于外圈故障特征频率的倒数。但齿轮啮合、电机电磁噪声、现场其他设备振动都会混进来直接看频谱往往是一堆峰值你不知道哪个是轴承特征频率。对同一段振动信号做自相关滞后轴上的第一个显著峰到原点的距离就是冲击间隔故障特征频率就是它的倒数。这个思路比在频谱里数峰值要稳健得多尤其适合早期微弱故障。3.3 自相关的计算与维纳-辛钦定理维纳-辛钦定理把时域和频域真正连在一起平稳随机信号的自相关函数与功率谱密度互为傅里叶变换对Sxx(f) ∫ Rxx(τ) · e^(-j2πfτ) dτ反过来说只要你有了自相关估计对它做傅里叶变换就能得到功率谱估计。反过来也一样。这个定理还解释了为什么经典功率谱估计里我们会先估计自相关再做变换。实际计算自相关常见三种方法直接按定义循环计算适合N小时做教学演示numpy的correlate函数直接调用适合中等长度FFT方法适合长序列。import numpy as np def autocorr_fft(x): n len(x) # 补零到2n避免循环回绕 X np.fft.rfft(x, n2 * n) corr np.fft.irfft(X * np.conj(X), n2 * n)[:n] return corr def autocorr_direct(x): n len(x) result np.zeros(n) x x - np.mean(x) for lag in range(n): result[lag] np.sum(x[:n - lag] * x[lag:]) if n - lag 0: result[lag] / (n - lag) return result注意直接实现时最后要除以有效长度n - lag否则滞后越大参与计算的样本越少累计值越小会形成一种虚假的衰减趋势。这涉及有偏和无偏估计的问题。除以n的是有偏估计方差更小但会低估除以(n - lag)的是无偏估计能纠正幅度但滞后增大时噪声会变大。实际工程里我倾向于用除以(n - lag)的无偏形式但只看滞后0到N/2的区域太靠后的部分样本太少参考价值有限。3.4 自相关的一个大坑低频趋势会造成假相关自相关对趋势项非常敏感。如果信号里有一个缓慢变化的线性趋势自相关结果会在很长一段滞后上保持较高值形成斜坡式衰减这会让周期峰被淹没或者出现虚假的长周期峰。所以做自相关前先把信号去直流、再去趋势已经是我的固定流程了处理振动、声音、生物电信号都适用。4. 互相关用对齐测时延服务声学定位与图像配准4.1 峰值位置等于时延估计互相关最经典的应用就是时延估计。两个麦克风同时录一个瞬时声音近的那个先到远的那个后到。把两段信号做互相关峰值所在的滞后位置就对应样点数的时延乘上采样周期得到秒再除以声速得到距离差。这是声学定位里TDOA方法的基本盘。代码里的实现非常直接import numpy as np def xcorr(x, y): # 用FFT方法算互相关返回滞后轴和归一化结果 n len(x) len(y) - 1 X np.fft.rfft(x, nn) Y np.fft.rfft(y, nn) corr np.fft.irfft(X * np.conj(Y), nn) lags np.arange(-(len(y) - 1), len(x)) return lags, corr使用时先找到corr的最大值位置再映射到lags数组上得到真实的滞后量。这里要特别注意不同库函数对lag轴的定义不一样numpy的correlate里modefull时返回数组的中间位置是0滞后计算时别把索引当滞直接用。4.2 归一化互相关与相位相关原始互相关的数值受信号能量影响幅度大的信号容易盖过幅度小的信号。为了消除这种影响可以做归一化互相关NCCNCCxy[τ] Rxy[τ] / sqrt( Rxx[0] · Ryy[0] )这个值被压到[-1, 1]之间更接近通常意义的相关系数。如果信号本身是图像或者处理的是频带很宽的信号可以用相位相关phase correlation。相位相关把互功率谱的所有幅度信息归一化掉只保留相位R_phase IFFT( conj(FFT(x)) · FFT(y) / |conj(FFT(x)) · FFT(y)| )这样做的好处是峰值非常尖锐时延估计分辨率更高尤其适合图像平移量的估计。做图像拼接时先做相位相关峰值位置就是两幅图之间的平移量亚像素级别的精度可以通过拟合峰来提高。4.3 窄带信号的多峰值陷阱互相关测时延有一个大坑容易踩如果信号是窄带周期信号比如50Hz正弦波互相关结果会在多个滞后位置都出现峰值间隔正好等于信号周期。这时候你没法从单一峰值判断真实时延因为你看到的可能是周期整数倍偏移后依然很像的假峰。这种情况在振动分析里非常常见设备转速稳定时振动主体是窄带周期成分两个测点信号之间计算互相关会有很多个等间隔的峰哪个都不是唯一答案。解决思路是结合信号的物理特性排除真假比如已知信号来源和测点距离范围把时延约束在一个合理区间在该区间内找峰或者使用宽带激励信号用相位相关等方法改善。知道这个坑比会算互相关更重要因为它会让你的结果稳很多。4.4 互相关与因果不能画等号互相关峰很高只能说明两个信号形态相似且时序对齐良好不能说明一条信号引起另一条。最常见的误导发生在一个共同源被两个传感器接收的场景传感器A和B都主要受到同一个外部振源影响它们的互相关自然接近1但A与B之间并没有任何物理上的单向因果关系。做故障定位或因果推断的时候需要引入系统模型、传递路径分析或者其他先验信息别只靠一个互相关峰下结论。5. 相干性把相关搬进频域看清哪个频段在连动5.1 从功率谱到相干性两个信号各自有功率谱密度描述能量在频率上的分布还有互功率谱密度描述两个信号在不同频率上的共同能量。把互功率谱按幅值平方归一化就得到幅值平方相干MSCCxy(f) |Sxy(f)|² / [ Sxx(f) · Syy(f) ]如果把Sxy、Sxx、Syy都换成某个频率点的复数/实数你会发现这个式子就是把频域的相关系数平方了。0表示该频率处毫不相干1表示完全线性相关。实际估计时Sxx和Syy用Welch方法做分段平均得到Sxy用同样的分段窗口取互谱平均。代码用scipy非常方便from scipy import signal f, Cxy signal.coherence(x, y, fsfs, nperseg1024, noverlap512, windowhann)参数上nperseg决定了频率分辨率noverlap会影响估计方差。经验上重叠50%以上比较稳定。频率分辨率和估计方差是一对矛盾nperseg越大频率分辨越细但分段数量越少方差越大nperseg越小方差越小但频率分辨变宽。具体选多少要看你想研究的谱线相差多少赫兹。5.2 为什么时域相关算完还要看相干性时域相关系数把全部频率成分混成单一指标会掩盖低频相关、高频不相关这样的频率选择性关系。相干性能把这种频率结构展示出来这在机械故障诊断里特别有用。举个例子减速箱输出轴转速40Hz啮合频率320Hz。两个测点如果相干性在40Hz及其谐波处接近1而在其他频率处很低说明两个测点之间的低频旋转激励是高度同步的这个信息单看时域互相关是得不到的。医学信号里也常用相干性看不同脑区的功能连接。两个脑电通道在alpha频段8-12Hz的相干系数高通常被视为这两个脑区在该频段存在功能耦合这也比单纯算一个总体相关系数信息量大得多。5.3 相干性的判读和误区MSC达到多少才算显著相关没有绝对标准和分段数量、窗函数、信号特性都有关系。工程上比较粗糙的经验是MSC大于0.8可以认为该频率成分两个信号之间存在明显的线性关系0.5以下基本不认。但如果你只分段了8次方差会很大0.8以上也有可能是噪声凑出来的。稳妥起见至少分16段以上再把0.8作为关注阈值。另一个经典误区是高相干等于有因果。相干性只能说明两个信号在某个频率上具有线性相关如果存在共同的驱动源即使两个测点之间没有任何信息传递相干性也会很高。例如管道系统两个压力测点都受水泵脉动影响在不同管路位置的测点压力在泵频处可能相干性接近1但你不能由此断定测点1的压力引起了测点2的压力。这个判断错误在工程现场出现过太多次了。5.4 相干性估计的一个实操细节去均值很多人计算相干性前不做去直流导致零频率附近的互谱巨大MSC在接近0Hz区间异常地高看起来好像两个信号在直流上有关系其实只是没去掉零频分量。我的习惯是不管算自相关、互相关还是相干性第一步先把均值减掉如果有线性趋势再用多项式拟合并剔除。这一步虽然看起来基础实际对结果稳定性的贡献比后面任何花哨参数都大。6. 一个完整的仿真实验自相关、互相关、相干性一通跑完6.1 构造仿真信号为了把前面的概念一次串起来我用Python构造一个接近工程场景的信号。设采样率1000Hz时长2秒。源信号是周期为0.08秒即12.5Hz的窄脉冲串模拟轴承故障产生的周期性冲击。再叠加宽带高斯白噪声让信噪比低一点。第二路信号是源信号延迟20个采样点后的版本再叠一部分独立噪声。import numpy as np from scipy import signal fs 1000 t np.arange(0, 2.0, 1/fs) f_bearing 12.5 # 源信号周期性窄脉冲 impulse_train np.zeros_like(t) impulse_train[::int(fs / f_bearing)] 1.0 impulse_signal np.convolve(impulse_train, np.hanning(15), modesame) # 加噪声并加入延迟 x impulse_signal 0.5 * np.random.randn(len(t)) # 第一路信号 delay_samples 20 y np.concatenate([np.zeros(delay_samples), impulse_signal[:-delay_samples]]) 0.5 * np.random.randn(len(t)) # 延迟版 x x - np.mean(x) y y - np.mean(y)6.2 自相关找周期直接算自相关画图后注意看滞后轴上的第一组峰值位置。周期性冲击的自相关会在lag0处最大在lag80点附近出现下一个峰因为周期是0.08秒对应80个采样点。这样就求出了周期再取倒数得到特征频率。corr_xx np.correlate(x, x, modefull)[len(x)-1:] lags np.arange(len(corr_xx)) / fs这里有个常见问题直接用np.correlate做自相关返回范围是[- (n-1), n-1]的完整结果取后半部分就得到滞后大于等于0的自相关。记得先做去均值再做相关否则0滞后附近的直流会把周期峰相对压低。6.3 互相关测时延对x和y做互相关峰值应该出现在lag0.02秒附近也就是20个采样点。直接把峰值所在lag乘上fs就是时延。注意如果使用了modefull滞后轴要从数组索引减去n-1。别看这个映射关系小实际编码里出错率极高。corr_xy np.correlate(x, y, modefull) lags_xy np.arange(-(len(x)-1), len(x)) / fs lag_est lags_xy[np.argmax(corr_xy)] print(f估计时延: {lag_est*1000:.2f} ms)如果信噪比很低互相关峰可能不明显甚至会出现峰位漂移。这时候可以试试把零均值和去趋势再做一遍或者用相位相关。我用相位相关处理过一个实测低信噪比声学信号效果明显好于普通互相关。6.4 相干性验证频段相关性用scipy.signal.coherence算MSC会看到在12.5Hz及其整数倍谐波处MSC接近1而其他频段接近0。这说明两路信号的关联主要集中在这个周期性冲击的频点及其谐波上。这个结果与互相关测出的20点延迟互相印证两路信号确实在周期性成分上高度一致只是存在固定时延。f, Cxy signal.coherence(x, y, fsfs, nperseg256, noverlap128)读数的时候注意只关心有物理意义的频段比如0-200Hz就够了高频段如果只是噪声MSC怎么算都会在0附近抖不用过度解读。6.5 这套流程在真实项目里的执行顺序我自己做振动或者声学数据分析时固定流程是先去除均值和趋势然后看自相关判断信号里有没有周期性结构用互相关估算各测点之间的时延关系最后用相干性验证不同测点之间在哪些频段存在显著的线性关系。这三步走完基本能对信号的空间和时间结构形成一个比较完整的判断比单独拿出一张频谱图靠肉眼猜要可靠得多。7. 踩过几次坑之后的个人习惯最后分享几点从项目里磨出来的操作习惯不一定写在教科书里但能帮你少走弯路。第一互相关的时延估计不要只取最大值就完事。先看互相关峰是锐峰还是圆峰。锐峰说明信号频带较宽时延估计可信圆峰说明信号窄带或者信噪比差峰值位置对噪声非常敏感。遇到圆峰我一般会改做窄带滤波后再算互相关或者直接用相位相关。第二自相关结果出现长滞后位置的高峰时先怀疑趋势项不要急着解读成某个长周期。把信号的线性趋势和直流去掉再跑一遍结果很可能完全变样。第三相干性的频率分辨率和方差必须同时看。某个频率点的MSC高但如果该频段本身能量极低那这个高相干也可能是数字计算在接近零值的两段功率谱之间做出的不稳定结果。判断的时候配合功率谱看如果两个信号在该频段都有明显能量相干性高才有实际意义。第四不管是自相关、互相关还是相干性结论都只能说存在线性关联不要过度推断因果。我见过不少案例把高相关直接等同于物理机制上的因果关系最后在定位故障源时走了弯路。做系统级判断时要把相关结果与设备结构、传递路径、工况信息结合起来才能得到可信的结论。整套方法吃透之后你手里会多出三样极为顺手的工具自相关用于找周期互相关用于测时延相干性用于做频段关联性分析。它们互相配合在振动诊断、声学定位、图像配准、生物信号分析这些领域里都属于绕不开的基本功。希望这篇分享能帮你把概念彻底打通少走我当年走过的弯路。
返回列表